Acta Polytechnica Vol. 52 No. 6/2012 Hierarchical Modeling of Mastic Asphalt in Layered Road Structures Based on the Mori-Tanaka Method Richard Valenta, Michal Šejnoha Department of Mechanics, Faculty of Civil Engineering, Czech Technical University in Prague Thákurova 7, 166 29 Praha 6, Czech Republic Corresponding author: sejnom@fsv.cvut.cz Abstract We present an application of the Mori-Tanaka micromechanical model for a description of the highly nonlinear behavior of asphalt mixtures. This method is expected to replace an expensive finite element-based fully-coupled multi-scale analysis while still providing useful information about local fields on the meso-scale that are not predictable by strictly macroscopic simulations. Drawing on our recent results from extensive experimental and also numerical investigations this paper concentrates on principal limitations of the Mori-Tanaka method, typical of all two-point averaging schemes, when applied to material systems prone to evolving highly localized deformation patterns such as a network of shear bands. The inability of the Mori-Tanaka method to properly capture the correct stress transfer between phases with increasing compliance of the matrix phase is remedied here by introducing a damage like parameter into the local constitutive equation of reinforcements (stones) to control an amount of stress taken by this phase. A deficiency of the Mori-Tanaka method in the prediction of creep response is also mentioned particularly in the light of large scale simulations. A comparison with the application of macroscopic homogenized constitutive model for an asphalt mixture is also presented. 1 Introduction Asphalt mixtures represent in general a highly het- erogeneous material with a complex microstructure consisting at minimum of a mastic binder, aggregates and voids. Our paper integrates various aspects of micromechanical modeling into a relatively simple, yet reliable and efficient computational scheme. The solution strategy relies on both uncoupled and cou- pled multi-scale homogenization approach. It com- bines advanced simulation-based homogenization of a certain representative volume element (RVE), e.g. the so-called statistically equivalent periodic unit cell (SEPUC), and classical micromechanics-based averag- ing techniques such as the Mori-Tanaka (MT) method, in a search for a reliable macroscopic constitutive law enabling a computationally efficient analysis of full- scale structures. Taking advantage of rapidly develop- ing imaging techniques the proposed approach follows the prominent bottom-up class of hierarchical model- ing promoting synergy of small scale experiments and detailed numerical simulations throughout individual scales to accomplish two basic guidelines, which on the one hand, considerably improve the predictive capability of numerical analysis and, on the other, provide significant cost savings in the material design. These are: • Take into account details of the microstructure of the analyzed heterogeneous material. • Avoid expensive and often intricate large-scale laboratory measurements. Implementing this task inevitably calls for the con- cept of a Virtual Testing Tool (VTT), which is cur- rently at the forefront of engineering interest. Details are available, e.g. in [6, 7]. When limiting our attention to Mastic Asphalt mix- tures, which are typically used in traffic arteries of substantial importance, the fraction of voids becomes negligible. This allows us to consider a two-phase com- posite consisting of irregular stones randomly spaced in a binding matrix. A color image of such a material system is plotted in Fig. 1(a). Its binary represen- tation is then also readily available, see Fig. 1(b). The literature offers several distinct routes that take advantage of this kind of representation. Following our previous works [9, 10], the approach adopted here combines various aspects of microme- chanical modeling into a relatively simple, yet reli- able and efficient computational scheme. While still incorporating the bottom-up uncoupled multi-scale homogenization strategy, the effective properties on individual scales are found as the volume averages of the local stress and strain fields developed inside SEPUC [11]. One particular example of SEPUC, representing the real microstructure of the asphalt mixture in Fig. 1(a), is shown in Fig. 1(c). Given such micromechanical models on every scale supported by associated experimental work, we may attempt to derive a macroscopic constitutive model that describes the homogenized response of Mastic Asphalt mixture (MAm) to general loading actions. The particular scales are shown in Fig. 2. It has been 48 Acta Polytechnica Vol. 52 No. 6/2012 (a) Original colour image (b) Simplified binary image (c) SEPUC ⇒ ⇒ Figure 1: Representation of asphalt mixture. Figure 2: Three distinct scales of Mastic Asphalt mixture. assumed that at each computational level the homog- enized response can be described by the nonlinear viscoelastic generalized Leonov (GL) model [7, 8]. Note that the necessary laboratory measurements are required only at the level of mastic asphalt. The remaining scales can be treated numerically through a set of virtual experiments performed on artificial computational models. While on the mortar scale, representing the scale of small stones removed from the original binary image, a simple hexagonal arrange- ment of stones is assumed, an advanced statistically equivalent representation of aggregates in terms of SEPUC is considered on the meso-scale. Such a unit cell is typically constructed by comparing the mate- rial statistics, e.g. the two point probability function, of the most appropriate representative of the real microstructure and the periodic unit cell [6, 7, 11]. This upscaling procedure consistent with Fig. 2 was applied in [7] to deliver the necessary parameters of the homogenized GL model. The essential results are summarized in Section 2 to show the potential of the homogenized model in practical applications, and also to motivate subsequent discussions. It is seen that reliable estimates of the macro- scopic response of real (large scale) structures can be obtained with the homogenized constitutive model. However, this approach lacks knowledge of the lo- cal mechanisms (the local stress and strain fields), which essentially drive the macroscopic response. On the other hand, a detailed fully coupled multi-scale analysis, typically employing the finite element dis- cretization on every scale, may prove computationally unfeasible. In this regard, the Mori-Tanaka averag- ing scheme appears to be a reasonable alternative for bridging the two above mentioned approaches. In the present context, the Mori-Tanaka method is introduced in a detailed finite element analysis at the level of material point to assess the time-dependent macroscopic response of an asphalt layer in a multi- layered road structure. It thus substitutes the macro- scopic constitutive model, which may not be available in general. Inasmuch as we expect that the local constitutive laws of individual phases are known, the derivation of macroscopic stresses and an instanta- neous homogenized stiffness matrix of asphalt is pro- vided by a very efficient MT homogenization step, see Section 3. At minimum, two-scale analysis is still maintained, but intensive finite element simulations at the level of SEPUC are avoided. The choice of the Mori-Tanaka method is first sup- ported by a detailed micromechanical analysis of the real microstructure of a Mastic Asphalt mixture show- ing reasonable agreement with a detailed finite ele- ment analysis of SEPUC. Its application in the nonlin- ear analysis of the layered system with subsoil defor- mation governed by one of the available constitutive models for soils is presented next in Section 4. A re- view of the essential findings is provided in Section 5. 2 Virtual testing tool This section summarizes essential steps for calibrating the macroscopic constitutive law in the framework of the Virtual Testing Tool outlined in detail in [7]. 49 Acta Polytechnica Vol. 52 No. 6/2012 (a) (b) Figure 3: (a) The HAAKE MARS dynamic shear rheometer; (b) loading scheme: plate-plate shear. The concept of VTT assumes that the experimental program is carried out only at the level of individual constituents — stones and mastic. At larger scales, a numerical simulation of virtual tests is performed. 2.1 Small scale experiments Although the mastic phase itself is a composite consist- ing of a filler and a bituminous binder, it is supposed in the present study to be well represented by a tem- perature and rate dependent homogeneous isotropic material. Since limiting our attention to moderate and elevated temperatures exceeding 0◦C the GL model is used to simulate the experimentally observed non- linear viscoelastic behavior of bituminous matrices. For the details of the proposed experimental program to calibrate the model parameters, and also for the actual material parameters of the model used in all the calculations the reader is referred to [6]. In short, the standard approach for model calibra- tion from a set of constant shear strain rate experi- ments accompanied by several creep tests at various temperatures [8] failed due to the inability of available experimental devices to provide the required viscous flow at zero elastic strain increment. A dynamic shear rheometer, displayed in Fig. 3, was therefore used to −100 0 100 −100 0 100 0 0.2 0.4 x1 x2 S m m (a) −100 −50 0 50 100 −100 0 100 0 0.2 0.4 x1 x2 S m m (b) Figure 4: (a) Target two-point probability function; (b) optimal two-point probability function. derive the required frequency and standard creep char- acteristics of the mastic phase. As shown in Fig. 3(b), this device operates on the plate–plate shear principle and allows the stress as well as the strain loading con- ditions to be applied either in simple rotation (shear) or in oscillatory mode at the prescribed temperature. The measured dynamic shear moduli obtained for the selected range of frequencies were used to con- struct the time-dependent compliance master curve of mastic asphalt, which in turn was used to calibrate the material parameters of the individual Maxwell units of the GL model. 2.2 Virtual numerical experiments At larger scales (mortar scale, asphalt mixture), the material parameters of the GL model were found through homogenization. The concept of first order homogenization of periodic fields was adopted to re- visit the original experimental program [8] numerically. Here, the key point is to formulate a reliable computa- tional model given in terms of SEPUC. For simplicity, we consider the scale of the asphalt mixture only. The analysis starts from the reference binary image of the asphalt concrete shown in Fig. 1(b), charac- terized by the two-point probability function of the 50 Acta Polytechnica Vol. 52 No. 6/2012 (a) (b) (c) (d) Figure 5: Examples of SEPUCs corresponding to the binary image in Fig. 1(b) and distributions of local shear strain due to the applied shear strain rate Ėxy: (a,b) SEPUC 37; (c,d) SEPUC 48. matrix phase shown in Fig. 4(a). Note that this im- age has undergone several image processing steps, described in detail in [6], that aim at preserving meso- scopic aggregates while filtering out very fine aggre- gate fractions. SEPUC, see Fig. 1(c), is defined by the following parameters: number of aggregates having elliptical shape, size, position, orientation and aspect ratio of the axes of individual ellipses. The size of the stones is derived on the basis of the cumulative distribution function [7]. For example, if 10 stones are selected for SEPUC, then the smallest stone corresponds to an average size of 10% of all stones determined from the cumulative distribution function. The next stone then reflects the size of the subsequent 10% of stones, etc. By adjusting the parameters of the aggregates to optimally approximate the target two-point probabil- ity function, we arrive at the idealized representation shown in Fig. 4(b). Since they were derived using a soft computing [4] strategy, the resulting SEPUCs may differ slightly in their microstructural details, though from the selected statistical descriptor point of view they are identical. Two particular examples appear in Fig. 5 together with the distributions of the local shear strains caused by the applied macroscopically uniform shear strain rate. Clearly, the localized shear bands in Figs. 5(b,c) are not identical as they follow the actual geomet- rical details of the corresponding SEPUC. It there- fore remains to confirm that all cells provide at least the “same” macroscopic response that is decisive for Figure 6: Macroscopic response for various SEPUCs, T = 20◦C, Ėxy = 10−4. the calibration of the homogenized macroscopic GL model. To answer this question, we refer the reader to Fig. 6, which depicts the macroscopic, homogenized, shear stress-strain curves attributed to the applied shear strain rate Ėxy = 10−4. Although no “perfect match” is observed, the difference in the estimated load bearing capacities does not exceed 10%. For comparison, we also plot the results provided by the Mori-Tanaka method. Although assuming only spherical inclusions for simplicity, these results still identify essential limitations of two-point averaging schemes — the inability to capture the localization phenomena observed in composites with a highly non- linear response of the binder phase [8]. A possible remedy, the main objective of this paper, will be offered in Section 3. 2.3 Model verification from large scale laboratory measurements We close this section by confirming the applicability of the derived macroscopic constitutive model through a numerical simulation of a large-scale tracking wheel test. Since details of the experimental measurements are available in [6], we limit our attention here to the simulation part. To that end, a two-dimensional plane strain model of an asphalt sample with the same cross-section as used in the laboratory measurements was created to perform the numerical analysis. It is imagined that during experiment, which was carried out according to EN 12 697-22, the wheel passing in the direction normal to the cross-section loads a 50mm wide zone (b) of the cross-section for a very short time of 1.2 s (one cycle of two passages lasts approximately 2.4 s). The idea is to model wheel tracking as short impulses of pressure applied to the loaded zone of the cross 51 Acta Polytechnica Vol. 52 No. 6/2012 Figure 7: Numerical model scheme of wheel tracking. section. Details of the computational model together with the applied load history appear in Fig. 7. By analogy with the actual experiment, the temperature was set equal to 40°C. Judging from the geometry of the structure and the deformed shape of the tested sample the actual width of contact area s to determine the applied pressure value was set equal to 20 mm which amounted to pressure p = 350 kPa. The dura- tion of pressure impulse ti = 0.2 s was calculated from the wheel speed and the length of contact area s. To reduce the computational time only one half of the structure was considered due to symmetry conditions. To minimize the number of required time steps, the actual loading diagram was simplified, and was assumed to attain a simple triangular shape as evident from Fig. 7. To reproduce the actual experiment, which was stopped after 2300 cycles of wheel tracking, the num- ber of impulses in Fig. 7 was set to 4600 (one impulse represents one passage of the wheel). In addition, a recovery period lasting one day was simulated. The results of the numerical calculation are plotted in Fig. 8. In Fig. 8(a) the gray curve represents the passing period while the black curve represents the recovery period. Note that only 2 hours of the recovery period are shown, for clarity of presentation. It cab be seen that the rut depth remains almost constant after this time. Its value after 24 hours reached 8.1 mm. The shaded region represents all short time impulses of the loading period. To better illustrate the loading se- quence we zoomed the first five and the last five passes now clearly visible in Figs. 8(c) and (d), respectively. Evidently, the rate of the increase in rutting slows down during each pass and appears almost negligible in Fig. 8(d), though the rut depth is still increasing as seen in Fig. 8(a). This Figure will be further discussed in the next section when we address the Mori-Tanaka method in conjunction with the creep response. A comparison between the experimental data and the numerical simulation is depicted in Fig. 8(b). Only a depression associated with the peak loads of passage of the wheel is displayed. Although experimentally observed response seems slightly more compliant, the (a) (b) (c) (d) Figure 8: Wheel tracking: (a) numerical simulation — loading and recovery, (b) comparison of maximum rut depth from numerical and experimental results, (c) first 5 and (d) last 5 passes from simulation. 52 Acta Polytechnica Vol. 52 No. 6/2012 numerical predictions appear rather appealing in the light of the simplified loading diagram. Note that more realistic short periods of creep (constant pres- sure between the loading and unloading branch) were excluded from the numerical simulation as they would cause the analysis computationally unfeasible. On the other hand, this may have a great impact on the final results to eventually arrive at even better agreement. In every case, the presented results show applicability of the proposed version of the GL model in the analysis of real engineering problems and con- firm the potential of the proposed two step uncoupled homogenization scheme as one particular step of the underlying Virtual testing tool. Recall, however, that no information about the lo- cal fields is available when purely macro-scale analysis is assumed. If we wish to keep the computational cost low, while still obtaining certain knowledge about the phase response, the Mori-Tanaka method discussed in the remaining part of this paper may become par- ticularly attractive. 3 Augmented Mori-Tanaka method In the preceding section the Mori-Tanaka method was quoted as a potential candidate for replacing de- manding finite element computations in hierarchical modeling of multi-layered road structures. However, direct application of this method in its original for- mat [1, 2] is often precluded by the prediction of an excessively stiff response in comparison to SEPUC analysis as has been already seen in Fig. 6. Thus in the present section, the local constitutive equations are reformulated to address: • Gradually decreasing the value of the matrix viscoelastic modulus. • The effect of strain localization already observed with a detailed FE analysis of SEPUC. 3.1 Theoretical formulation If admitting only piecewise uniform variation of the stress and strain fields within individual phases, the stress increments become ∆σs = ωLs∆εs = L̂s∆εs (1) ∆σm = Ltm∆εm + ∆λm, (2) where Lr is the stiffness tensor of the phase r = s,m (subscript s,m being assigned to stones and matrix, respectively) and the increment of the eigenstress ∆λm represents contributions caused by the creep strains developed in the binder phase thus leaving the aggregates to remain linearly elastic. Superscript (t) denotes an instantaneous quantity. The damage like parameter ω in Eq. (1) is introduced to reduce the stress carried by the stones when localized strains are developed in the matrix. It is imagined that such a behavior, which cannot be addressed by standard formulation of the Mori-Tanaka method, can be phys- ically associated with a gradual debonding of the aggregates allowing the matrix to take all the stress when the composite deformation is realized solely through a network of localized shear bands in the matrix. Note that ω = 0 corresponds in this case to a fully debonded material. The following representation is proposed ω = [ τ teq Nτ0 / sinh ( τ teq Nτ0 )]M , (3) where τeq is the current equivalent deviatoric stress andM,N are model parameters. In the present study, these were found by comparing the MT predictions with the results provided by the homogenized GL model, recall Section 2, under the strain control con- dition for the prescribed macroscopic shear strain rate Ėxy. The local strain increments can be estimated from Dvorak’s Transformation Field Analysis [2] to get ∆εs = Ats∆E +Dtsm∆µm, (4) ∆εm = Atm∆E +Dtmm∆µm, (5) where Atr,Dtrm, r = s,m, are the mechanical strain localization tensors and the transformation influence functions, respectively. The phase eigenstrain µm is introduced to account for the nonlinear viscoelastic behavior of the binder. Note that for a two-phase medium, the transformation influence functions are directly available in terms of the strain localization tensors [2] as Dtsm = ( I −Ats ) ( Ltm −Ls )−1 Ltm, Dtmm = ( I −Atm ) ( Ltm −Ls )−1 Ltm. (6) Also note that the strain concentration tensors and the transformation influence functions satisfy the fol- lowing universal connections [2] N∑ r=1 crAr = I, (7) N∑ r=1 crDrs = 0, (8) N∑ r=1 Dsr = I −As, (9) N∑ r=1 DsrMr = 0, (10) N∑ r=1 crLrDrsMs = cs ( I −AT s ) , (11) 53 Acta Polytechnica Vol. 52 No. 6/2012 N∑ r=1 N∑ s=1 crLrDrsMs = 0, (12) where N denotes the number of phases (N = 2 in our particular case) andMr = L−1 r is the compliance matrix of a given phase r. To finalize the list of needed expressions, we provide the specific forms of the strain concentration factors Atm = [ cmI + csT t s ]−1 , (13) Ats = T tsA t m, (14) where the partial strain concentration factor T ts is given by T ts = [ I − St(Ltm)−1 (Ltm −Ls)]−1 . (15) Closed form solutions for the Eshelby tensor St are provided in the literature for certain types of hetero- geneities, see e.g. [5]. Next, writing the macroscopic constitutive law ∆Σ = Lhom∆E + ∆Λ, (16) and realizing that ∆Σ = cs∆σs + cm∆σm, (17) yields the macroscopic stiffness matrix and the macro- scopic eigenstress [2] Lhom = csL̂sA t s + cmL t mA t m, (18) ∆Λ = cs ( L̂sD t sm ) ∆µm +cm [ Ltm ( Dtmm − I )] ∆µm (19) = − [ cs ( L̂sD t smM t m ) + cm ( LtmD t mmM t m )] ︸ ︷︷ ︸ 6=cs(I−AT s), recall Eq. (11) ∆λm. Note that the connection given by Eq. (11) is no longer valid due to the presence of the modified matrix L̂s. 3.2 Numerical predictions A comparison of the resulting macroscopic predictions including the results from various SEPUCs are avail- able in Fig. 9(a) for the uniform macroscopic shear strain rate Ėxy. The evolution of damage parameter ω for optimal values of M = 0.3, N = 5 is plotted for illustration in Fig. 9(c). Considerable improvement in the behavior of the MT method is evident partic- ularly when compared with the original predictions presented in Fig. 6. Another strong motivation for mastering the MT method is the possibility to estimate at least the local phase averages if it is used in place of the macroscopic homogenized constitutive law in a multi-scale like analysis. For illustration, these estimates are com- pared separately with predictions from SEPUC 37 in Fig. 9(e). Having in mind that the parameters of the damage model were fitted against the homogenized response, the agreement of the local fields is quite satisfactory. A closer match might be expected if deriving the damage model parameters directly from the local matrix stress averages provided by one of the corresponding periodic unit cells. It is also advisable to test the proposed homoge- nized GL model against more complex loading condi- tions. As an example we consider the simple tension represented by the prescribed macroscopic stress rate Σ̇xx. The available homogenized and local distribu- tions appear in Figs. 9(b,d,f). The results indicate several observations not evident for pure shear load- ing. For the most important observation, we draw the readers attention to Fig. 9(b), which shows a rel- atively large deviation of the macroscopic predictions delivered by individual SEPUCs. This suggests that the assumed size of SEPUC is insufficient for SEPUC to be considered a unique RVE since small changes in a microstructure cause quite severe differences in macroscopic predictions. Fig. 9(d) demonstrates the evolution of damage parameter ω. It is not surprising to observe a very slow evolution of this parameter for this type of loading since only a portion of the total load accrues to a deviatoric stress. Note that Ed = √ 2eijeij is the equivalent deviatoric strain. Finally, we inspected the predictive capability of the MT method for loading conditions dominating the creep or relaxation response. For illustration, two other unit cells plotted in Fig. 10(a,b) were selected to carry out the comparative finite elements based simulations. The creep response is examined in Fig. 10(c). There is an almost perfect match between the results found from the homogenized GL model and from the finite element analysis of SEPUC No. 43, since this partic- ular RVE was used in [7] to derive the homogenized master curve in full from a set of creep experiments. It is not so much the deviation between the two unit cells but, rather, the predictions pertinent to the Mori-Tanaka method associated with the interval of constant load (creep response) that are disturbing. The inability of the method to predict the correct creep behavior is quite pronounced. However, the loading branches and also the unloading breaches controlled by the damage parameter ω are captured relatively well. Somewhat better agreement is pro- vided by the strain control loading visible in Fig. 10(d). Nevertheless, the differences between local fields are still quite severe. Possible explanation to this poor performance of the Mori-Tanaka method is provide in [6]. Further evidence for the inability of the two-point averaging variant of the Mori-Tanaka method to pre- dict creep behavior correctly is evident from Fig. 11, which compares the results from the tracking wheel 54 Acta Polytechnica Vol. 52 No. 6/2012 (a) (b) 0.00 0.05 0.10 0.15 0.20 Strain E xy [-] 0.60 0.70 0.80 0.90 1.00 ω [- ] Strain control E xy 0.00 0.01 0.02 0.03 0.04 Strain E d [-] 0.94 0.96 0.98 1.00 ω [- ] Stress control Σ xx (c) (d) (e) (f) Figure 9: (a,b) Macroscopic response for various SEPUCs, the homogenized GL model and the MT method; (c,d) variation of damage parameter ω; (e,f) phase averages of local fields from SEPUC 37 and the MT method; (a,c,e) T = 20°C, Ėxy = 10−4 s−1; (b,d,f) T = 20°C, Σ̇xx = 0.025 kPa. 55 Acta Polytechnica Vol. 52 No. 6/2012 (a) (b) (c) (d) Figure 10: (a) SEPUC No. 6; (b) SEPUC No. 43: macroscopic response for reference temperature T = 40°C; (c) stress control: loading–creep–unloading– recovery; (d) strain control: loading–relaxation– unloading–recovery. experiment derived by the homogenized GL model and by the Mori-Tanaka method. 4 Numerical simulation of a layered road structure Both the homogenized GL model and the MT method were implemented in the GEO FEM [3] finite element code to allow for the modeling of real engineering structures. It is not the purpose to analyze one par- ticular structure but rather to assess the feasibility 0 1 2 3 4 5 6 Time [s] -0.20 -0.15 -0.10 -0.05 0.00 0.05 R u t d e p th [ m m ] Leonov hom Mori-Tanaka First 5 passings Wheel tracking (reduced load) Figure 11: Numerical model scheme of wheel track- ing. Comparison between the predictions based on the homogenized GL model and the Mori-Tanaka method. Figure 12: Scheme of flexible pavement structure. of the two approaches when combined with inelastic deformation of the subsoil. A simplified three-layer structure was therefore considered. The soil profile also evident from Fig. 12 consists of one layer of asphalt 100 mm in thickness and two lay- ers of subsoil 400 mm in thickness. The upper subsoil layer represents a heavily dense gravel overlaying a layer of sand. The two subsoil layers were represented by the Drucker-Prager constitutive model with the material data listed in Table 1. Unsuitability of the present asphalt mixture to be used in road applications, see [7], required adjusting the material parameters of the GL model on the mortar scale to comply with 0◦C temperature in order to arrive at a converged solution even for a relatively short duration of loading. The wheel load F = 5 kN, corresponding to a stand- ing personal vehicle, was assumed to be distributed over a 200× 200mm surface area, which corresponds to the prescribed surface pressure p = 125 kPa, see Fig. 13(a) also showing the boundary conditions and the finite element mesh. We observed that even for 0°C temperature the asphalt material showed an enor- mous susceptibility to creep causing the loaded area to 56 Acta Polytechnica Vol. 52 No. 6/2012 Soil Young modulus Poisson number Weight Cohesion Friction angle E [MPa] ν [-] γ [kN/m3] c [kPa] ϕ [-] Gravel 50 0.2 21 5 41.5 Sand 30 0.3 17.5 5 31.5 Table 1: Material properties of subsoil layers (a) (b) (c) (d) Figure 13: (a)–(b) Displacement dz [mm]; (c)–(d) equivalent plastic strain Epld [%]. be essentially squashed in all the way to the bottom of the asphalt layer within just a few seconds. Although a converged solution was achieved, the results are not realistic owing to the assumed small strain theory. The results that correspond to almost instantaneous loading (the total load was applied in 1 s) are shown in Figs. 13(a,c), whereas Figs. 13(b,d) were found after just 3 s for permanent loading reaching the peak load. In particular, variation of the vertical displace- ment appears in Figs. 13(a,b). The distribution of the equivalent deviatoric plastic strain in the subsoil layers is displayed in Figs. 13(c,d). An increase in the localized plastic deformation due to creep even for such a short duration of permanent loading is evident. Also note a considerable increase in vertical displacement from 2 to 5 mm. Bear in mind, however, some limitations of the MT method for predicting creep. Nevertheless, similar results were presented in [6] using the homogenized GL model. Although it fails to predict unwanted rutting caused by long-term permanent loading, the analysis still indi- cates that the proposed homogenized GL model or the MT method can be used for large scale applications. 5 Summary Although research interests on flexible pavements have been quite intense in the past two decades, the field is still very much under development and will cer- tainly witness considerable activity in the coming decade particularly in connection with hierarchical modeling and micromechanics. Within this frame- work, the present work provides theoretical tools for formulating a macroscopic constitutive law reflecting the confluence of threads coming from experimental work, image analysis, statistical mechanics and tra- ditional disciplines of micromechanics and first order computational homogenization. Section 2 reviewed the totally uncoupled multi-scale 57 Acta Polytechnica Vol. 52 No. 6/2012 modeling approach for the prediction of macroscopic homogenized constitutive model that will enable inex- pensive analyses of real world large scale structures. If we are interested not only in the bulk response but also in the local phase behavior, the stress update procedure can be performed in the framework of the Mori-Tanaka method described in Section 3. This es- sentially resembles a two-scale analysis, but in terms of efficiency this method is superior to considerably more time consuming finite element simulations at the level of SEPUC. The results presented here favor this approach but at the same time suggest caution in applications where strong creep behavior is expected. Nevertheless, a number of points still deserve further attention. The effect of microstructure anisotropy, the influence of the first stress invariant and the three-dimensional character of asphalt mix- tures are just a few issues which need to be addressed, yielding flexible pavements a fertile field of future research. Acknowledgements Financial support provided by MSM 684077003 is gratefully acknowledged. References [1] Y. Benveniste. A new approach to the application of Mori-Tanaka theory in composite materials. Mechanics of Materials 6:147–157, 1987. [2] G. J. Dvorak, Y. Benveniste. On transformation strains and uniform fields in multiphase elas- tic media. Proceedings of the Royal Society of London Series A - Mathematical, Physical and Engineering Sciences 437(1900):291–310, 1992. [3] Fine-Ltd. GEO 5 manual. http://www.fine. cz/, 2012. [4] A. Kučerová. Identification of nonlinear mechan- ical model parameters based on softcomputing methods. Ph.D. thesis, Ecole Normale Supérieure de Cachan, Laboratoire de Mécanique et Tech- nologie, 2007. [5] T. Mura. Micromechanics of Defects in Solids. No. 3 in Mechanics of elastic and inelastic solids. Kluwer Academic Publishers, second revised edn., 1987. [6] R. Valenta. Modeling of Asphalt Mixtures. Ph.D. thesis, Czech Technical University in Prague, 2011. [7] R. Valenta, M. Šejnoha, J. Zeman. Macroscopic constitutive law for mastic asphalt mixtures from multiscale modeling. International Journal for Multiscale Computational Engineering 8(1):131– 149, 2010. [8] M. Šejnoha, R. Valenta, J. Zeman. Nonlinear viscoelastic analysis of statistically homogeneous random composites. International Journal for Multiscale Computational Engineering 2(4):645– 673, 2004. [9] M. Šejnoha, J. Zeman. Overall viscoelastic re- sponse of random fibrous composites with sta- tistically quasi uniform distribution of reinforce- ments. Computer Methods in Applied Mechanics and Engineering 191(44):5027–5044, 2002. [10] J. Zeman, M. Šejnoha. Numerical evaluation of effective properties of graphite fiber tow im- pregnated by polymer matrix. Journal of the Mechanics and Physics of Solids 49(1):69–90, 2001. [11] J. Zeman, M. Šejnoha. From random microstruc- tures to representative volume elements. Mod- elling and Simulation in Materials Science and Engineering 15(4):S325–S335, 2007. 58