Homogenization of the elastic-viscoplastic damage behavior of asphalt mixtures based on the mesomechanical Mori–Tanaka method

To improve the efficiency and accuracy in characterizing the nonlinear behavior of multi-phased asphalt mixtures and to further facilitate asphalt pavement design based on the nonlinear properties, this study proposes a homogenized method for the elastic-viscoplastic-damage constitutive model. The Drucker–Prager-type yield surfaces and non-association flow rule were employed to describe the viscoplastic strain of asphalt materials. The continuum damage mechanics (CDM) were incorporated to characterize the evolutions of internal micro-cracks or micro-voids in structures. The mesomechanical Mori–Tanaka (MT) method was used to yield the homogenized nonlinear strain components within multi-phase structures. The proposed constitutive model was then implemented in finite element (FE) simulations based on a self-developed subroutine. Several case studies were conducted, in which composite structures with one inclusion were simulated under constant stress and strain loading rate. Amongst the composite structures, the inclusions with various volume contents and shapes were taken into account. In addition, the influence of air voids in structures was considered by defining the zero stiffness for inclusions. The results indicated that the nonlinear behavior of composite with single aggregate or air void can be effectively represented using the proposed method. Furthermore, by coupling the homogenized nonlinear constitutive relation into the locally homogenous model from previous work, not only was the viscoplastic-damage behavior of the composite with multiple inclusions effectively demonstrated by the definition of the nonlinear material, but the internal heterogeneous feature was also precisely represented by the local cells. Therefore, the proposed homogenized method can effectively predict the viscoplastic and damage behavior of asphalt mixtures with various inclusions by simply specifying the parameters in the FE simulations.


Introduction
Asphalt mixtures account for a large proportion of constructed roads due to their high flexibility and load-bearing capacity. In engineering practice, asphalt mixtures exhibit nonlinear viscoplastic and damage behavior under traffic loadings, which pose great challenges to pavement analyses and designs. In decades, categories of stress-strain relations have been developed for describing the mechanism of asphalt mixtures according to the phenomenological test results [1][2][3][4]. However, such relations were based on conducting laborious experiments; moreover, asphalt mixtures were regarded as homogenous structures instead of composite structures. Therefore, characterizations for asphalt mixtures have been mostly untenable or limited.
To account for the internal structure of asphalt mixtures, mesomechanical theories are presented and widely employed in homogenizing the material performance of composites by introducing the "Eshelby tensor" [5]. During this process, inclusions are assumed to be "dispersed" in the matrix and their contribution to strain components are considered using the local and global strain concentration tensors. Based on the Eshelby tensor, the dilute approach assumes that inclusions are embedded in an infinity space, where the distance 1 3 between two particles is far enough [5][6][7]. However, in the so-called dilute method, the interactions between inclusions are greatly underestimated. For asphalt mixtures comprising multiple inclusions, such a dilute method is inappropriate. To this end, the Mori-Tanaka (MT) method [8] was proposed to better predict the global performance of composite structures with "dense" distributed inclusions. The MT method calculates the average strain for the matrix before embedding each particle, and hence considers the interactions between multiple inclusions. Furthermore, the self-consistent and generalized self-consistent methods [9][10][11] were developed to accurately predict the elasticity of heterogeneous structures using the implicit equations. The abovementioned research facilitated the mechanical characterization of composite materials considering their internal constituents. Nevertheless, the mesomechanical homogenizations yield the global parameters of composites comprising multi-phase under linear configurations.
For the nonlinear performance of structures, such as plastic or damage behavior, several research efforts have been conducted. Mercier and Molinari [12] proposed a homogenization scheme based on the interaction law [13], in which the disordered materials were characterized by two self-consistent schemes, and the composite materials were presented by a MT model. Compressible elasticity and anisotropy of the materials were taken into account in a general framework. Wu et al. [14] presented an incremental secant mean-field homogenization (MFH) for describing the elastoplastic behavior of composites. This approach estimated the residual strain in each phase during the unloading process, and then predicted the strain increments in each phase using the developed secant operators. Zecevic and Lebensohn [15] developed three homogenization approaches for predicting the elasto-viscoplastic deformation of polycrystals, which were suitable for the compression of copper, tension of stainless steel, and compression of magnesium. Wang et al. [16] formulated a multiscale incremental constitutive model to represent the breakage mechanism for frozen soil and gave a more reasonable theoretical function for the volumetric breakage ratio by thermodynamic theory. Such an analytical solution for plasticity homogenization enhanced the capability in characterizing nonlinear behavior of multi-phase structures. However, in engineering practice, the analytical results can rarely be applied for realistic structure prediction and design.
Numerical simulation allows constructing physical models of complex structures, and further discretizing them into finite or discrete elements. By solving the partial differential equations within finite elements or the interaction force between discrete elements, the stress-strain responses of the entire structure can be precisely demonstrated, which greatly reduces the labor consumptions in analyzing the mechanical performance of structures in laboratory or field tests. In addition, varieties of constitutive behaviors can be specified in the model parameters to represent different materials. For example, commercialized finite element (FE) programs, such as ABAQUS, enable users to defined any type of material behavior in a self-developed subroutine by computer programming [17,18]. Based on the subroutines, any kind of constitutive material model can be specified for the use within FE simulations. Subsequently, plastic homogenizations were realized in simulations. Doghri and Friebel [19] improved the Hill-type incremental formulation [20] for elastoplastic composite and linked it to the FE program by a DIGIMAT/UMAT module. Bouhamed et al. [21] endeavored to homogenize the elastoplastic functionally graded material (FGM) through the single point incremental forming (SPIF) process and implemented the homogenization equations into ABAQUS using the subroutines (VUMAT and VUSFLD). Ogierman [22] presented a two-stage hybrid homogenization approach combining the mean-field approach and FE simulation, which guarantees high accuracy and time efficiency for characterizing the elastic-plastic functionally graded composite. Therefore, it is reasonable and necessary to introduce the numerical simulations into the viscoplastic and damage characterization to demonstrate the nonlinear behavior of asphalt materials and further predict the performance of asphalt pavements with high efficiency and effectiveness.
For the nonlinear characterization of asphalt materials, numerous studies have been completed towards the viscoplastic and damage mechanism investigations. Zhu et al. [23][24][25] did lots of work in developing the viscoelasticviscoplastic damage constitutive model of asphalt mixture based on the thermodynamic theory. In their research, the Mohr-Coulomb yield criterion was used to represent the pressure dependency of asphalt mixtures. To fulfill the non-associated flow rule of asphalt materials, the Drucker Prager model was employed as the plastic potential, which can be used to define the viscoplastic strain increment rate. The damage evolution law in their work was proposed according to [26], which was coupled with the viscoplastic internal variable. Additionally, a research group at Texas A&M University developed abundant nonlinear analyses of asphalt mixtures [27][28][29][30][31][32][33][34]. Unlike Zhu's model, they used a Drucker-Prager-type equation to describe both yield surface and viscoplastic potential function but employed different parameters for the two functions to represent the non-associated flow rule for asphalt mixtures. In addition, continuum damage mechanics (CDM) were used to represent the internal damage evolutions including micro-cracks or micro-voids. The abovementioned two constitutive models both regarded the asphalt mixtures as homogeneous materials and applied laboratory tests to curve-fit the viscoelastic and viscoplastic strains. However, they failed to predict the nonlinear behavior of mixtures at smaller scales, where the internal structures of asphalt mixtures depend on the randomly distributed aggregates, and hence the nonlinear behavior cannot be well considered in the asphalt mixture design and analysis.
The previously introduced mesomechanical approaches and viscoplastic constitutive relation provided comprehensive background for homogenizing the nonlinear behavior of composite materials. Nevertheless, such homogenized viscoplastic characterizations for asphaltic materials have rarely been reported. To this end, the mesomechanical based homogenization approach was conducted towards the elastic-viscoplastic-damage constitutive relation of asphalt mixtures. To simplify the simulation process, the asphalt matrix was defined as elastic material before satisfying the yield criteria, and the Drucker-Prager-type equations with different model constants were used to define the yield surface and the viscoplastic potential function, respectively. The simple two-phase composite structure was developed to stand for a small unit of an asphalt composite, where the matrix and inclusions were respectively defined as elasticviscoplastic and linear elastic, as can be observed in Figure 1. Furthermore, the CDM was involved to represent the internal damage evolutions within asphalt mixtures. The mesomechanical MT method was used to homogenize the abovementioned elastic-viscoplastic-damage material properties of the composite structure. The constitutive model was finally packaged in ABAQUS with an UMAT subroutine. To verify the proposed approach, the composite structure and the homogenous structure were both developed in FE simulation. The composite structure consisted of an elasticviscoplastic-damage asphalt matrix and elastic inclusions, whereas the homogenous structures represented the composite material with homogenized elastic-viscoplastic-damage properties. Therefore, for numerical nonlinear simulation on asphalt mixtures, various inclusions can be easily considered by specifying different model parameters. In the end, the proposed constitutive relation was coupled with the locally homogeneous model from the previous researches, and the nonlinear behavior of different asphalt mixtures can be efficiently and effectively predicted at global and local scales.

Notation
In this study, the boldface symbols denote tensors, matrices, and vectors. The component-wise can be expressed as, Einstein's summation convention is used in the tensor operations, The second and fourth-order identity tensors can be expressed as, where δ ij is the Kronecker's symbol,

Homogenization
The general stress-strain constitutive relationship of the elastic materials at mesoscale can be formulated as, where σ(x), L(x), and ε(x) refer to the stress tensor, stiffness tensor, and strain tensor of materials at mesoscale coordinates, respectively.
Assuming that a representative volume element (RVE) withstands homogeneous strain E on its boundary, the strain tensor at mesoscale can be represented as, where A(x) is the so-called global strain concentration tensor, which is related to the mechanical properties and geometrical features of the RVE. Therefore, Eq. (10) can be transformed as, By homogenizing the stress tensor at the volume of RVE, the homogeneous stress tensor Σ can be represented as, where L hom is the homogeneous stiffness tensor of RVE, and it can be obtained by, For r-th phase of the composite structure (matrix or inclusion), the homogenizations of L(x) and A(x) can be calculated to determine its homogeneous strain tensor, which is denoted as L r and A r . Therefore, the homogeneous stiffness tensor of RVE can be formulated as, where c r is the volume content of the r-th phase in RVE; ∑ r c r A r =I.
To calculate the global strain concentration tensor A r , a relationship in strains between matrix and inclusion phases is considered: where r and 0 are homogeneous strain in r-th inclusion and matrix, respectively; G r is the local strain concentration tensor. Furthermore, the globally homogeneous strain can be formulated as, Therefore, the relationships in strains between global and local scales can be determined as, Substituting Eqs. (16)- (19) into Eq. (11), and the global concentration tensor A r can be identified as, Finally, the homogeneous stiffness tensor of the composite structure can be expressed as, In the MT mesomechanical method, the local concentration tensor G r can be explicitly determined by the Eshelby tensor S r of the r-th inclusion [35], where the Eshelby tensor S r can be identified by the geometrical feature of the inclusion. Therefore, as soon as the shape and volume of the inclusion is identified, the homogeneous stiffness of the composite can be easily calculated.

Continuum damage mechanics
To represent the micro damage evolution in structures, the CDM was introduced and applied by [36], which defines the damage variable ϕ as, where A and Ā are the apparent area and undamaged area carrying the load. After micro damages (including microcracks or micro-voids) initiate in structures, the undamaged area Ā , also called the effective area, would decrease from A to zero, which causes the damage variable ϕ to increase from 0 to 1. The damage variable is widely used to account for the damages in the stress-strain constitutive relations. Based on the CDM theory, [37,38] defined the effective stress in the damaged structures, which accurately characterizes the degradation in strength and stiffness because of the damage evolution, where ̄ refers to the effective stress tensor under the undamaged configurations; σ is the nominal Cauchy stress tensor in the damaged configuration.
In the subsequent sections, the effective stress components were used in the constitutive model formulation since they derive more elastic and viscoplastic deformations. Therefore, in the following equations, the superimposed "-" was added above the stress components to represent the effective configuration. In the numerical simulations, a relationship should be defined between stress σ and strain ε. According to [39], strain components in nominal and effective configurations have little difference under small deformation and isotropic damage assumptions. Therefore, the strain equivalence hypothesis was adopted in this study, which assumes that the nominal strain tensors are equal to the corresponding effective strain tensors. The superimposed "-" was also added to the strain components in the subsequent formulations.

Elastic-viscoplastic constitutive model
In this study, a composite structure consisting of the elasticviscoplastic matrix phase and the linear elastic inclusion phase was considered. For the elastic inclusion phase, the stress-strain relation can be described by Hook's law. For the matrix phase, the total strain can be decomposed into the recoverable (elastic) part and the irrecoverable (viscoplastic) part concerning the stress state, where is the total strain tensor; e and vp are elastic strain and viscoplastic strain tensors, respectively. In this study, the viscoplastic strain rate can be identified by Perzyna's model [30,40], where Γ is the viscosity parameter, f is the yield function, and ϕ(f) represents the overstress function of f, ϕ(f) = (f/σ 0 ) N ; σ 0 and N are material constants; 〈〉 is the Macaulay brackets, in which 〈x〉 = x when x > 0, and 〈x〉 = 0 when x < = 0; the superimposed dot (•) indicates the derivation to time. The Drucker-Prager yield surface was employed to take into account the confinement, deviatoric stress, and dilative deformation, which are presented in I 1 -τ plane, is the deviatoric stress tensor; α is material constant relating to the internal friction of materials; δ ij is the Kronecker delta; Ī 1 =̄k k is the first stress invariant; (̄v p e ) is the isotropic hardening function, which can be expressed as [41], where κ 0 , κ 1, and κ 2 are material constants; vp e is the scalar viscoplastic strain, which represents the internal state variable of materials. The rate of vp e can be determined as [42], The non-associated viscoplastic flow rule was used herein to avoid overestimating the dilative behavior [43,44], which defines that the direction of viscoplastic strain increment is not normal to the yield surface. To simplify the simulation, the Drucker-Prager-type model was used as the viscoplastic potential function with different coefficients from the yield surface, where β is the model parameter and smaller than α.

Damage evolution
The evolution law for the damage variable refers to [30]. To simply the simulation process, the evolution of the damage variable in this study was assumed as temperature-independent, which gives, where Γ ϕ is the reference damage viscosity parameter; Y 0 is the reference damage force; q is the stress-dependency parameter; k is the material parameter; eff are the effective strain parameters, which can be determined by eff = √ ij ij ; − Y is the effective damage force, which represents the energy

Numerical implementations
To guarantee the convergence of the simulation and meanwhile satisfy the consistency condition, Euler's method was used to update the stress state based on the strain increments [45][46][47].
At the moment of t n , all the necessary components of the structures, ij t n , ij t n , vp ij t n and vp e t n , have been determined from the previous steps. In addition, the total strain increment Δ̄i j t n+1 and time increment Δt are given by the numerical program. According to Euler's method, the basic strain and stress increments component-wise in the numerical simulation are listed below, However, for composite materials consisting of elastic inclusions and viscoplastic matrix, the viscoplastic strain increment would only exist in the matrix. Therefore, according to Eq. (17), the total elastic strain increment can be written as, where Δ e/m ij t n+1 and Δ e/r ij t n+1 are elastic strain increments of matrix and r-th inclusion, respectively. Based on Eqs. (18) and (19), their elastic strain increments can be respectively represented as, where Δ̄m ij t n+1 and Δ̄r ij t n+1 are total strain increments in the matrix and r-th inclusion; Δ̄v p ij t n+1 is the viscoplastic strain increment that is only in the matrix phase. Hence, Eq. (36) and Eq. (35) can be re-written as, According to the MT mesomechanical theory, the homogenous stiffness tensor of composite structures can be determined as, where L m ijkl and L r ijkl are stiffness tensor of matrix and r-th inclusion, respectively; G r ijkl and S r ijkl are the local concentration tensor and Eshelby tensor of r-th inclusion.
Therefore, the fundamental variable for updating the stress state depends on calculating the viscoplastic strain increment. From Eq. (26), the viscoplastic strain increment can be determined by, To simplify the simulation process, the viscoplastic multiplier Δγ vp was introduced herein, Therefore, the viscoplastic strain and scalar viscoplastic strain increments can be formulated as, The differentiation of g to stress components can be written as, (39) At the time moment of t n+1 , strain increment Δ ̄i j (t n+1 ) will be firstly assumed as elastic, and the trial stress components can be calculated as Hence, the trial yield surface can be calculated by the trial stress components, Afterwards, a judgment should be made towards the trial yield surface. If f tr < = 0, no viscoplastic strain would be induced, and the stress state of the structure would be updated by the trial stress; if f tr > 0, the calibrated stress components would be calculated to consider the viscoplastic strain increment, Substituting Eq. (49) into Eq. (51), the calibrated stress should be, Therefore, the stress state can be solved as soon as the viscoplastic strain increment is determined. From Eq. (45), the viscoplastic strain increment can be identified by the viscoplastic multiplier Δγ vp . The calculation of Δγ vp is presented in the following parts.
Firstly, the deviatoric viscoplastic strain increment to Δγ vp can be formulated based on Eq. (45), The calibrated deviatoric stress component can be derived as, where G hom is the homogeneous shear modulus. The second deviatoric stress invariant after calibration can be expressed as, The calibrated first stress invariant is, where K hom is the homogeneous bulk modulus. The relationship between K hom , G hom , and L hom ijkl can be expressed as L hom ijkl = K hom δ ij δ kl + G hom (δ ik δ jl + δ il δ jk -2/3δ ij δ kl ). Meanwhile, the isotropic hardening function can be updated by Δγ vp , Based on the calibrated stress component, the yield surface can be modified as, The consistency condition of viscoplasticity can be defined according to [48,49], and hence a dynamic yield surface can be determined by Eq. (44) and Eq. (59), The Newton-Raphson (N-R) iteration method is employed to calculate Δγ vp (t n+1 ), in which the differential of χ to Δγ vp is expressed as, Substituting Eqs. (55)-(57) into Eq. (61), the Eq. (61) can be expressed as, T h e " ± " i n E q . ( 6 2 ) i s d et e r m i n e d by � Before the N-R iteration, an initial value of Δγ vp (t n+1 ) should be defined. At the k + 1 iteration, the viscoplastic multiplier can be calculated as, where Δγ vp (t n+1 ) can be finally determined when χ k+1 is smaller than the tolerance error. Afterwards, the viscoplastic strain increment can be expressed by Eq. (45). Therefore, the elastic strain increment and stress increment can be consequently updated.
For the damage simulation, it is noteworthy that the damages only occur in the matrix phase within composite structures, and the damage variable can be updated by, where, where ̄m eff is the effective strain parameter in the matrix phase, which can be given by, Subsequently, the nominal stress can be calculated by Eq. (24).
In the numerical simulation, the derivative of stress increment to strain increment should be calculated to derive the Jacobian matrix in UMAT. The differential of stress increment is expressed as, According to Eq. (45), the differential of the viscoplastic strain increment can be expressed as, Substituting Eq. (68) into Eq. (67), the differential of the stress component can be expressed as where, where ∂ 2 g/(∂∆σ ij ∂∆σ kl ) can be formulated as, The plastic consistency condition is used herein to explicitly represent dΔγ vp (t n+1 ), which defines the differentiation, Therefore, Eq. (72) gives, Substituting Eq. (73) into Eq. (69), and finally determined dΔγ vp (t n+1 ) as, Substituting dΔγ vp (t n+1 ) into Eq. (69), and Δ̄∕ Δ̄ can be expressed as, where df 1 /dσ ij (t n+1 ) and df 2 /dΔγ vp (t n+1 ) can be calculated according to Eq. (72), Therefore, the Jacobian matrix can be determined, which can be further used for numerical simulations. The entire (76) Fig. 2 Flow chart of homogeneous elastic-viscoplastic simulation process of the present homogeneous elastic-viscoplastic damage numerical simulation is presented in Fig. 2,

Finite element simulations
The simple inclusion-matrix composite model was developed in the commercial finite element program Abaqus (Version 2017). Within the model (see Fig. 3), the inclusion was defined as elastic sphere or ellipsoid with different sizes; the matrix phase was defined as elastic-viscoplastic damage material, where the model parameters were specified by the UMAT subroutine using the Fortran programming. Meanwhile, a homogeneous model with the same size as the composite model was established, in which the material parameters were homogenized by the proposed method. Furthermore, the locally homogeneous model was developed to account for the multiple inclusions. The material parameters were listed in Table 1, where the elastic properties of the inclusions referred to typical values; the viscoplastic and damage properties of the matrix were determined according to [30]. To simplify the homogenization process and reduce the influence of viscosity behavior before the yield surface, the elastic properties of the matrix were defined according to the viscoelastic Prony series in [30], where the equilibrium modulus of the asphalt material was used as the elastic Young's modulus for the matrix phase.
In this study, the stress-control and strain-control loading modes with constant loading rates were applied to the models, and the inclusions with different volume contents and different shapes were taken into consideration. In addition, the air voids in the model were regarded as inclusions with zero stiffness and further considered in the simulations.

Single spherical inclusion
The simulation results of spherical inclusions with different volume contents were exhibited in Fig. 4. The results present that the proposed method is capable of describing the non-linear viscoplastic damage behavior of composite materials with homogeneous parameters. For a composite structure with smaller inclusion (0.5% volume content), the stress-strain relations of composite structure are close to the matrix without inclusions, which indicates that the influence of inclusion with low volume content can be neglected in the nonlinear performance. With the increase of the inclusion contents, the stiffness of the structure is increased, and the stress-strain curves presented by the proposed homogeneous model fit the composite structure very well. In addition, the difference of results from the two loading modes became more significant with larger inclusions. For composites comprising large inclusion (27% volume content), the simulation process becomes unstable, especially under strain-controlled loading mode. Such a phenomenon can be explained by the discontinuous strain behavior within the structure.
To find a possible explanation for the unstable simulation on structures with larger inclusions, the contour plots of scalar viscoplastic strain ̄v p e were extracted and presented in Fig. 5, in which the volume contents of inclusion are 3% and 27%. It can be observed that, compared with 27% inclusion, the viscoplastic strain distributions in the composite structure with 3% inclusion are more similar to the corresponding homogeneous structure. Due to the embedded elastic inclusion, remarkable strain concentration behavior can be observed surrounding the inclusion,   and the strain concentration behaves more significantly in a composite structure with larger inclusion (volume content of 27%), which leads to the unstable stress-strain results and even causes convergence problems. After homogenizing the viscoplastic properties by the proposed method, the strain concentration behavior was avoided due to the uniform model parameters. Therefore, the proposed model can effectively represent the homogenous viscoplastic behavior of composite at the global scale; meanwhile, the simulation efficiency of the nonlinear viscoplastic simulation can be greatly improved. However, such homogenizations underestimated the interfacial behavior between matrix and inclusions, which makes it difficult to deeply investigate the stress-strain behavior of composite in smaller scales, especially for multiple inclusions. To address this issue, the previously developed locally homogenous model was introduced in the following part to represent the asphalt mixture structure with multiple inclusions, which can characterize the homogeneous plastic behavior of asphalt composite in the local parts based on its internal structures.
In the current simulations, the internal damage variables within structures were too small (around 1E-8) and their influence on the mechanical responses can be neglected. However, to exhibit the effectiveness of the developed method in characterizing the damage behavior of composite structures, the contour plots of damage variable ϕ distributions in composite and homogenous structures with 3% inclusion were extracted and presented in Fig. 6. The values of damage variables show great agreement between the composite structure and homogenized structure, which indicates that the proposed method is capable of describing the nonlinear damage evolution in composite structures.

Single ellipsoid inclusion
The stress-strain curves of the composite consisting of ellipsoid inclusion with an aspect ratio of 2 are exhibited in Fig. 7, which represents that the proposed method can effectively characterize the viscoplastic behavior of ellipsoidbased composites. However, more simulations on composites with higher volume contents and larger aspect ratios of inclusions are not possible because the strain concentration behavior causes the convergence problems, which will be presented in detail in the following parts. Figure 8 presents the internal effective viscoplastic strain distributions caused by inclusions with different aspect ratios. The higher aspect ratio of ellipsoids induces more viscoplastic strain concentrations near the tip of the inclusions, which would further reduce the simulation efficiency. Therefore, the composite structure model is not appropriate to address the nonlinear viscoplastic behavior of composites with flat and elongated particles. Homogenizing the viscoplastic behavior of composites with such inclusions would affect the accuracy of the nonlinear predictions, and hence the ellipsoidal inclusions with large aspect ratios were neglected in the proposed approach. In fact, it is reasonable to conduct this process because the realistic engineering practice generally remove such inclusions with an aspect ratio larger than 2 to guarantee the performance of the structure.

Air void
To prove the capability of the proposed model in addressing the influence of air voids, structures consisting of different air void contents were developed and simulated. The air void contents in this simulation were defined as 0.5%, 2.83%, and 3%, which was based on the requirements on the dense-graded asphalt mixtures (AC) in engineering practice (around 4%). The stress-strain curves under two different loading modes are exhibited in Fig. 9. The results from the homogeneous structures are close to the composite structures, which proves that the proposed method can

Multiple inclusions
As mentioned above, the homogenization process failed to consider the interfacial behavior and interactions between different phases, which reduces the effectiveness of the proposed method in addressing the nonlinear behavior especially for asphalt mixtures with multiple inclusions. Realistic composite materials are comprised of randomly distributed inclusions, which pose great changes to the nonlinear characterization. To this end, a composite structure with random inclusions was developed for viscoplastic simulation. Correspondingly, two different homogeneous models were established, which, respectively were the globally homogeneous model and the locally homogeneous model [35]. The three models are presented in Fig. 10. The composite structure contains matrix phases and inclusion phases. To simplify the simulation, the inclusion phases represent stones with elastic properties. In the globally homogeneous model, the material parameters were homogenized by considering all the inclusions in the composite structure. In the locally homogeneous model, the space was divided into local cells based on the locations and volumes of the inclusions. The local separation process was conducted using the weighted Voronoi tessellation, according to the authors' previous work [50]. In The simulation results under stress and strain loading modes were presented in Fig. 11. The stress-strain curves from the three models are similar to each other, which proves the effectiveness of the proposed homogenization method in considering the multiple inclusions. It should be noted that the simulation on the composite structure stopped before reaching the final point in the stress-controlled loading mode, whereas the other two models finished their simulations steadily. This phenomenon indicates that the proposed method can describe the nonlinear behavior of asphalt mixture more efficiently. Furthermore, the scalar viscoplastic strain distributions within the three models were extracted and presented in Fig. 12. Compared to the globally homogeneous structure, the locally homogeneous model can better demonstrate the internal heterogeneous stress-strain distribution and meanwhile increase the simulation efficiency. Moreover, by incorporating the locally homogeneous models from previous work, the internal heterogeneous structure of composite, such as the locations and sizes of the inclusions, can be effectively characterized. The abovementioned results indicate that the proposed homogenization method can demonstrate the viscoplastic performance of composite materials with both high accuracy and fewer convergence problems. Therefore, the nonlinear and heterogeneous performance of composite materials can be effectively presented.

Conclusions and outlook
This study presents a homogenization method for characterizing the elastic-viscoplastic and damage behavior of asphalt composites based on the mesomechanical theories. The Drucker-Prager model and non-associated flow rule were employed to determine the viscoplastic strain components. Before the viscoplastic strain was induced, linear elasticity was defined to simplify the simulation and to reduce the influence of the viscoelasticity of asphaltic materials. The CDM was used to represent the internal micro-cracks and micro-voids of structures. The MT method was used to homogenize the mechanical properties of the composite, and the nonlinear constitutive model was defined in FE simulation by the user-material UMAT subroutine. Using this model, the elastic-viscoplastic damage behavior of asphalt composite with various inclusions can be simulated by modifying the model parameters instead of laborious model developments. To verify the effectiveness of the proposed method, case studies with different inclusions were simulated under stress-control and strain-control loading modes. The conclusions are listed as following: (1) The proposed method shows a powerful ability in homogenizing the viscoplastic behavior of composite structures. For the lower contents of inclusions, the simulation results from homogeneous models and composite models were almost equal to each other. For the higher contents of inclusions, simulations on the composite structures were unstable and even failed with convergence problems, which is caused by the strain concentration behavior near the inclusions. Such phenomenon reduced the simulation efficiency and increase the risk of convergence problems in characterizing the nonlinear behavior of composite materials. However, the proposed homogenization method can accurately represent the viscoplastic behavior of composites without inducing too many strain concentrations, and further can characterize the plastic performance more efficiently. (2) The damage evolutions in the models were too small to affect the mechanical responses; however, the contour plots of damage variable distribution indicate that the proposed method can effectively consider the damage evolutions of composites. (3) Inclusions with various shapes were effectively considered in this study. Specifically, spherical inclusions and ellipsoidal inclusions with a smaller aspect ratio (smaller than 2) were accurately represented in the homogeneous viscoplastic model. For ellipsoidal inclusions with a larger aspect ratio (great than 2), remarkable strain concentrations would be induced near the sharp-tips of the inclusion in the composite structures and thus reduce the simulation efficiency. (4) The air voids can be considered in the proposed method by defining as the inclusions with zero stiffness. The simulation results from the proposed model and the composite structure showed little difference when the air void contents were small (smaller than 4%). (5) The locally homogeneous model from previous research was introduced. By cooperating with the proposed homogenization method, the locally homogenous model can represent the nonlinear viscoplastic behavior of asphaltic composite with high efficiency and accuracy; meanwhile, the heterogeneous internal structure feature of the composites can be effectively considered.
The fundamental theories and corresponding case studies for viscoplastic damage homogenization were presented in this study. However, asphalt composites exhibit various internal structures and stand complex loading conditions. In future research, the nonlinear behavior of the large-scale asphalt mixture model with randomly distributed inclusions should be simulated together with the locally homogenous model. Moreover, the viscoelastic-viscoplastic and damage material properties of asphalt mortar and asphalt mixture will be determined by laboratory tests. According to the test data, the effectiveness of the proposed method in realistic engineering can be further proved. Therefore, the simulation results can be used for predicting the nonlinear performance of the real asphalt mixture structures, and further for facilitating the asphalt pavement design.