Hierarchical Modeling of Mastic Asphalt in Layered Road Structures Based on the Mori-Tanaka Method

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.


Introduction
Asphalt mixtures represent in general a highly heterogeneous 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 coupled multi-scale homogenization approach.It combines 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 averaging techniques such as the Mori-Tanaka (MT) method, in a search for a reliable macroscopic constitutive law enabling a computationally efficient analysis of fullscale structures.Taking advantage of rapidly developing imaging techniques the proposed approach follows the prominent bottom-up class of hierarchical modeling 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 concept of a Virtual Testing Tool (VTT), which is currently at the forefront of engineering interest.Details are available, e.g. in [6,7].
When limiting our attention to Mastic Asphalt mixtures, which are typically used in traffic arteries of substantial importance, the fraction of voids becomes negligible.This allows us to consider a two-phase composite 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 representation 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 micromechanical modeling into a relatively simple, yet reliable 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  assumed that at each computational level the homogenized 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 arrangement 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 material 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 macroscopic response of real (large scale) structures can be obtained with the homogenized constitutive model.However, this approach lacks knowledge of the local 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 discretization on every scale, may prove computationally unfeasible.In this regard, the Mori-Tanaka averaging 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 multilayered road structure.It thus substitutes the macroscopic 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 instantaneous homogenized stiffness matrix of asphalt is provided 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 supported by a detailed micromechanical analysis of the real microstructure of a Mastic Asphalt mixture showing reasonable agreement with a detailed finite element analysis of SEPUC.Its application in the nonlinear analysis of the layered system with subsoil deformation governed by one of the available constitutive models for soils is presented next in Section 4. A review of the essential findings is provided in Section 5.

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].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.

Small scale experiments
Although the mastic phase itself is a composite consisting of a filler and a bituminous binder, it is supposed in the present study to be well represented by a temperature 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 nonlinear 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 calibration from a set of constant shear strain rate experiments 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 derive the required frequency and standard creep characteristics 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 conditions 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 construct 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.

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 revisit the original experimental program [8] numerically.Here, the key point is to formulate a reliable computational 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 matrix phase shown in Fig. 4(a).Note that this image has undergone several image processing steps, described in detail in [6], that aim at preserving mesoscopic aggregates while filtering out very fine aggregate 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 probability 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 geometrical details of the corresponding SEPUC.It therefore remains to confirm that all cells provide at least the "same" macroscopic response that is decisive for 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 nonlinear response of the binder phase [8].A possible remedy, the main objective of this paper, will be offered in Section 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 50 mm 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 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 duration of pressure impulse t i = 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 number 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 sequence 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 numerical predictions appear rather appealing in the light of the simplified loading diagram.Note that more realistic short periods of creep (constant pressure 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 confirm 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 local 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 particularly attractive.

Augmented Mori-Tanaka method
In the preceding section the Mori-Tanaka method was quoted as a potential candidate for replacing demanding finite element computations in hierarchical modeling of multi-layered road structures.However, direct application of this method in its original format [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.

Theoretical formulation
If admitting only piecewise uniform variation of the stress and strain fields within individual phases, the stress increments become where L r 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 physically 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 where τ eq is the current equivalent deviatoric stress and M, 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 condition for the prescribed macroscopic shear strain rate Ėxy .The local strain increments can be estimated from Dvorak's Transformation Field Analysis [2] to get where A t r , D t rm , 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 Also note that the strain concentration tensors and the transformation influence functions satisfy the following universal connections [2] N r=1 where N denotes the number of phases (N = 2 in our particular case) and M r = 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 where the partial strain concentration factor T t s is given by Closed form solutions for the Eshelby tensor S t are provided in the literature for certain types of heterogeneities, see e.g.[5].
Next, writing the macroscopic constitutive law and realizing that yields the macroscopic stiffness matrix and the macroscopic eigenstress [2] ∆λ m .
Note that the connection given by Eq. ( 11) is no longer valid due to the presence of the modified matrix L s .

Numerical predictions
A comparison of the resulting macroscopic predictions including the results from various SEPUCs are available 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 particularly 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 compared 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 homogenized GL model against more complex loading conditions.As an example we consider the simple tension represented by the prescribed macroscopic stress rate Σxx .The available homogenized and local distributions appear in Figs.9(b,d,f).The results indicate several observations not evident for pure shear loading.For the most important observation, we draw the readers attention to Fig. 9(b), which shows a relatively 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 E d = 2e ij e ij 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 particular 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 provided 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 predict creep behavior correctly is evident from Fig. 11, which compares the results from the tracking wheel

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 particular structure but rather to assess the feasibility  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 layers 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 standing personal vehicle, was assumed to be distributed over a 200 × 200 mm 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 enormous susceptibility to creep causing the loaded area to  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 displacement 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 indicates that the proposed homogenized GL model or the MT method can be used for large scale applications.

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 certainly witness considerable activity in the coming decade particularly in connection with hierarchical modeling and micromechanics.Within this framework, 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 traditional disciplines of micromechanics and first order computational homogenization.Section 2 reviewed the totally uncoupled multi-scale modeling approach for the prediction of macroscopic homogenized constitutive model that will enable inexpensive 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 essentially 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 mixtures are just a few issues which need to be addressed, yielding flexible pavements a fertile field of future research.

Figure 2 :
Figure 2: Three distinct scales of Mastic Asphalt mixture.

Figure 5 :
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.

Figure 7 :
Figure 7: Numerical model scheme of wheel tracking.

Figure 8 :
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.

Figure 11 :
Figure 11: Numerical model scheme of wheel tracking.Comparison between the predictions based on the homogenized GL model and the Mori-Tanaka method.

Table 1 :
Material properties of subsoil layers