Investigation on Microstructural Damage Properties of Asphalt Mixture Using Linear and Damage-Coupled Viscoelastic Model

Abstract

This paper presents an image-based micromechanical modeling approach for simulating the damage-couple viscoelastic response of asphalt mixture. Details of the numerical damage-couple viscoelastic constitutive formulation implemented in a finite element code are presented and illustrated by using the ABAQUS user material subroutine (UMAT). Then, an experimental procedure based on the Linear Amplitude Sweep test for obtaining the viscoelastic and damage parameters at a given temperature was conducted. An improved morphological multi-scale algorithm was employed to segment the adhesive coarse aggregate images. We developed a pixel-based digital reconstruction model of asphalt mixture with the X-ray CT image after being processed. Finally, the image-based FE model incorporated with damage-coupled viscoelastic asphalt mastic phase and elastic aggregates was used for the compressive test simulations successfully in this study. Simulation results showed that the damaged simulation results have a larger stress distribution compared with the undamaged simulation due to the irregularity of the coarse aggregates. The von Mises stress distribution is smaller as the loading time increases due to the viscoelastic behavior of asphalt mastic. It can also provide insight on the damaged mechanisms and the possible location in asphalt mixture where microscopic cracking would most likely occur.

1. Introduction

Asphalt mixture is the most widely used road building material in the pavement industry. The mechanical behavior of the asphalt mixture is complex due to the heterogeneity of the asphalt composite material and the time- and temperature-dependent viscoelastic binder. In the aim of optimizing materials design, it is necessary to study the behavior of this composite material.

The mechanical response of rheology-like materials can be well-predicted using the linear viscoelastic constitutive model expressed by convolution integrals [1,2,3,4,5,6]. In fact, the linear viscoelastic theory is widely used to predict the fundamental properties of asphalt mixture. For example, Keshavarzi et al. [7] proposed a new method for studying the viscoelastic behavior of asphalt mixture in uniaxial tension. Gudmarsson et al. [8] compared the complex modulus and complex Poisson’s ratio determined through modal testing with cyclic tension–compression measured complex moduli and complex Poisson’s ratio of asphalt mixture to further apply the simpler and more economic modal testing for characterizing the linear viscoelastic properties of asphalt mixture. Pszczola et al. [9] performed the bending beam creep test at temperatures from −20 °C to +10 °C in order to evaluate the low-temperature performance of asphalt mixtures. Furthermore, García et al. [10], Zhang et al. [11], and Huang et al. [12] have developed numerical schemes for implementation of linear viscoelastic constitutive model in finite element (FE) codes to simulate the mechanical properties of asphalt mixture under certain conditions. However, different material’s response models, based on visco-elastic-plastic equations implemented on finite element codes [13,14,15,16], rather than with the discrete element method [17,18] can be found in the literature.

Asphalt mixture can be considered to be consist of three phases: the mastic phase (asphalt cement with a filler that is smaller than 2.36 mm in size), the aggregate phase; and the air voids. The asphalt mastic phase is a typical viscoelastic material which gives the asphalt mixture its rheological characteristics. Recently, many researchers have applied advanced techniques to capture the internal microstructure of asphalt mixtures in order to understand the properties of asphalt mixture from a microscopic perspective [19,20,21,22,23,24,25,26]. Since the asphalt mastic is composed of micro-cracks, micro-voids, and defects that result in strength loss, local cracks in the asphalt mixture would start to form near existing cracks and, thus, would lead to the mixture failure as loading is increased. The damage behavior of asphalt mixture is extremely complex due to its microstructure. Several theoretical attempts have been made to investigate the damage behavior of asphalt mastic under certain conditions. Dai [27,28] proposed a rate-independent damage-couple viscoelastic model by replacing two nonlinear parameters in Schapery nonlinear viscoelastic model to predict the microstructural damage properties of asphalt mixture. Kim et al. [29] employed a computational micromechanics damage-induced model for predicting damage behavior of asphalt mixtures. Zeng et al. [30,31] proposed a two-phase microstructural numerical simulation framework to analyze damage behaviors of heterogeneous asphalt mixture under low temperature. In this framework, the asphalt mastic damage is considered to be only relative to the deviatoric stress. In these studies, the damage density is treated as a constant at each Gaussian integral points for each incremental step. In fact, asphalt mixture is a particulate composite material consisting of aggregate, mastic, and air voids. The mechanical behavior of asphalt mixture strongly depends on the particular components. During the damage evolution process in heterogeneous materials, the amount of damage in each region of the material is related to the current strain. The damage density distribution of heterogeneous materials under loading is discrete in incremental step.

In order to predict the micromechanical response of asphalt mixture at a constant temperature, very few damage models have been coupled to linear viscoelasticity and damage constitutive equations considering the discrete damage density for each incremental step. This paper presents an image-based micromechanical modeling approach for simulating the damage-couple viscoelastic response of asphalt mixture at a given temperature. We propose a damage-coupled viscoelastic model based on the Linear Amplitude Sweep (LAS) test and then implement the constitutive model in a finite element code using the ABAQUS user material subroutine (UMAT). In this constitutive model, the damage density is defined as a function of the current strain. After that, we develop an experimental procedure for obtaining the viscoelastic and damage parameters at a given temperature in the linear viscoelasticity range. Moreover, a pixel-based method is used to generate the microstructural digital sample of the asphalt mixture based on X-ray CT image. To overcome the adhesive coarse aggregate images, an improved morphological multi-scale algorithm based on several different sizes of structural elements is introduced. Finally, the image-based numerical sample combined with the ABAQUS user material subroutine will be applied to predict the viscoelastic behavior of a real microstructure asphalt mixture specimen.

2. Theoretical Background

Asphalt mixture is a particulate composite material consisting of coarse aggregate, asphalt mastic, and air voids. The coarse aggregate phase is basically considered as elastic material in nature. The asphalt mastic phase is a typical viscoelastic material which gives the asphalt mixture its rheological characteristics. In this section, the linear viscoelastic constitutive model used to represent the behavior of asphalt mastic phase is introduced. Then, a displacement-based damage-coupled viscoelastic model for asphalt mastic is proposed as a function of the current strain. Finally, the incremental formulations of the constitutive model implemented in a finite element code are developed in details.

2.1. Linear Viscoelastic Model

For linear viscoelastic materials, including asphalt mastic, the stress–strain constitutive relation is expressed by convolution integrals. In the case of stress response at constant strain, the convolution relations is explained in Equation (1) for one-dimensional problems:

$$\sigma \left(t\right)=E_{\infty }\epsilon +{{\displaystyle \int }}_0^t\Delta E\left({\phi }^t-{\phi }^{\tau }\right)\frac{d\left(\epsilon \right)}{d\tau }d\tau$$

(1)

where the $E_{\infty }$ is the instantaneous elastic modulus, ${\phi }^{\tau }$ is the reduced time and ΔE is the transient modulus. It is given by Equation (2):

$$\Delta E\left(\phi \right)={\displaystyle \sum }_{n=1}^N{E}_n\exp \left(-{\rho }_n\phi \right)$$

(2)

where $E_n$ is the n-th coefficient of the Prony series and ${\rho }_n$ is the n-th relaxation time.

For 2D and 3D problems, Equation (1) can be decomposed into deviatoric and volumetric components in Equation (3) to Equation (5):

$$S_{ij}^t=2G_0{e}_{ij}^t+2{{\displaystyle \int }}_0^t\Delta G\left({\phi }^t-{\phi }^{\tau }\right)\frac{d\left(e_{ij}^{\tau }\right)}{d\tau }d\tau$$

(3)

$${\sigma }_{kk}^t=3K_0{\epsilon }_{kk}^t+3{{\displaystyle \int }}_0^t\Delta K\left({\phi }^t-{\phi }^{\tau }\right)\frac{d\left({\epsilon }_{kk}^{\tau }\right)}{d\tau }d\tau$$

(4)

$${\sigma }_{ij}^t=S_{ij}^t+\frac{1}{3}{\sigma }_{kk}^t{\delta }_{ij}$$

(5)

where $S_{ij}$ and ${\sigma }_{ij}$ are the deviatoric stress tensor and volumetric stress tensor, respectively; ${\sigma }_{kk}^t$, $G_0$, and $K_0$ are the total stress tensor, the instantaneous effective elastic shear modulus, and bulk modulus, respectively; $\Delta G$ and $\Delta K$ are the transient shear modulus and bulk modulus, respectively. ${\delta }_{ij}$ is the Kronecker delta.

In most cases, the asphalt mastic Poison’s ratio μ is assumed to be time independent [32,33,34]. So all the material properties in Equations (3) and (4) can be determined from the data of uniaxial tensile tests in Equation (6):

$$\Delta G_0=\frac{E_0}{2\left(1+\mu \right)},K_0=\frac{E_0}{3\left(1-2\mu \right)},\Delta G\left(t\right)=\frac{\Delta E\left(t\right)}{2\left(1+\mu \right)},\Delta K\left(t\right)=\frac{\Delta E\left(t\right)}{3\left(1-2\mu \right)}$$

(6)

2.2. Damage-Coupled Viscoelastic Model

Kachanove [35] proposed the concept of continuum damage mechanics to describe the progressive decay of materials, where Rabotnov [36] introduced the continuity ζ to define the reduction of area due to micro-damages in Equation (7):

$$\varnothing =1-\zeta =\frac{A-\overline{A}}{A}$$

(7)

where A is the damaged (apparent) area and $\overline{A}$ is the real area (intact or undamaged area) carrying the load. In other words, $\overline{A}$ is the resulted effective area after micro-damages (micro-cracks and micro-voids) are removed from the damaged area A. $\varnothing$ is the so-called damage density or damage variable.

For cases of the isotropic (scalar) damage, the effective stress is given by Equation (8):

$${\overline{\sigma }}_{ij}^t=\frac{{\sigma }_{ij}^t}{1-\varnothing \left(t\right)}=\frac{S_{ij}^t+\frac{1}{3}{\sigma }_{kk}^t{\delta }_{ij}}{1-\varnothing \left(t\right)}$$

(8)

where ${\overline{\sigma }}_{ij}^t$ is the effective stress tensor in the effective (undamaged) configuration.

Generally, the damage density can be a function of stress, strain, temperature, and damage history. However, a damage law considering various influencing factors is not appropriate for predicting the damage growth in asphalt mixture since it is very difficult to express these influences quantitatively in one equation due to their mutual influence. Rabotnov [36] proposed a damage evolution law in terms of strain under exponent format, where he eliminated other influence factors from the evolution law. Belloni et al. [37] proposed an approximation for the damage variable relying that the strain is the most important factor. Dai et al. [27] proposed a unified damage-couple viscoelastic model by using the Schapery model, where the damage variables are functions of maximum strain. Darabi et al. [32] proposed the first approximation of the damage evolution law as an exponential form of the total effective strain.

As for asphalt binder damage evolution law, the Linear Amplitude Sweep (LAS) test is proposed by Johnson to specify the fatigue resistance [38,39]. The LAS test is a cyclic sinusoidal torsion test that uses increasing loading amplitude to accelerate damage. The damage accumulation (D) in the asphalt binder sample can be calculated using in Equation (9) to in Equation (12) [40]:

$$D\left(t\right)\cong {\displaystyle \sum }_{i=1}^N{\left[\pi {\gamma }_0^2\left(C_{i-1}-C_i\right)\right]}^{\frac{\alpha }{1+\alpha }}{\left(t_{i-1}-t_i\right)}^{\frac{1}{1+\alpha }}$$

(9)

$$G\prime \left(\omega \right)\cong \left|G^{*}\right|\left(\omega \right)\times cos \delta \left(\omega \right)$$

(10)

$$log G\prime \left(\omega \right)=m\left(log \omega \right)+b$$

(11)

$$\alpha =1/m$$

(12)

where C_i ≅ |G∗|(i)/|G∗|_initial which is the complex shear modulus for a given data point divided by its initial undamaged value and C_i-1 is the complex shear modulus divided by its initial undamaged value at a former data point.${\gamma }_0$ is the applied shear strain amplitude at each loading cycle and α is the exponent determining energy release rate. $\left|G^{*}\right|\left(\omega \right)$ and $\delta \left(\omega \right)$ are the dynamic modulus and phase angle at frequency ω, respectively. m is the slope of a best-fit straight line plotting with log ω on the horizontal axis and log $G\prime \left(\omega \right)$ on the vertical axis.

Many numerical attempts have been made to investigate the damage behavior of asphalt mixture under certain conditions. However in these numerical models, the asphalt mixture is considered as being a homogenous material with equivalent properties. Thus, the damage density is treated as a constant at each Gaussian integral points for each incremental step. Asphalt mixture is a particulate composite material consisting of aggregate, mastic, and air voids. The mechanical behavior of asphalt mixture strongly depends on the particular components. During the damage evolution process in heterogeneous materials, the amount of damage in each region of the material is related to the current strain. Thus, the damage density distribution of heterogeneous materials under loading is discrete in incremental step.

Motivated and guided by the aforementioned damage evolution laws and the LAS test, in this study, the damage density of heterogeneous materials such as asphalt mixture under loading is proposed as a function of the current strain. For the incremental formulation, the damage evolution law can be calculated using Equation (13), such that:

$$\varnothing \left(t\right)=\varnothing \left({\epsilon }_{t-\Delta t}\right)+\stackrel{•}{\varnothing }\left({\epsilon }_t\right)\Delta t$$

(13)

where $\varnothing \left(t\right)$ is the damage density at time t. $\varnothing \left({\epsilon }_{t-\Delta t}\right)$ is the damage density when the strain value is ɛ at time t−Δt. $\stackrel{•}{\varnothing }\left({\epsilon }_t\right)$ is the rate of damage density at time t with a strain value of ɛ.

According to the LAS test, a normalized damage density versus strain curve is produced as shown in Figure 1.

For a given strain at the current incremental step, the rate of the damage density can be calculated using the normalized damage density versus strain curve. Details of constructing the normalized damage density versus strain curve will be presented in Section 3.2, which is the main focus of this paper.

2.3. Numerical Implementation

The finite element method (FEM) is actually an incremental approach for numerical analysis. Current stress and strain at integration points of each element at every time increment are obtained from the stress and strain over the previous loading history. So the incremental deviatoric and volumetric formulations can be derived with some algebraic manipulations.

Substitution of Equations (2) and (6) into Equation (3), we can obtain the deviatoric stress tensor for cases of the isotropic damage. The results are expressed as Equations (14) and (15):

$${\overline{S}}_{\mathrm{ij}}\left(t\right)=\frac{{2\mathrm{G}}_e{e}_{\mathrm{ij}}\left(t\right)+2{{\displaystyle \int }}_0^t{{\displaystyle \sum }}_{n=1}^N{G}_n{e}^{-\left(t-\zeta \right){/\rho }_n}\frac{{\mathrm{de}}_{\mathrm{ij}}}{\mathrm{d}\mathsf{\zeta }}\mathrm{d}\mathsf{\zeta }}{1-\varnothing \left(t\right)}=\frac{{2\mathrm{G}}_e{e}_{\mathrm{ij}}\left(t\right)+2{{\displaystyle \sum }}_{n=1}^N{G}_n{q}_{\mathrm{ij},\mathrm{n}}^t}{1-\varnothing \left(t\right)}$$

(14)

$$q_{ij,n}^t={{\displaystyle \int }}_0^t{e}^{-\left(t-\zeta \right)/{\rho }_n}\frac{d{e}_{ij}\left(\zeta \right)}{d\zeta }d\zeta$$

(15)

Substitution of Equations (2) and (6) into Equation (4), we can obtain the volumetric stress tensor for cases of the isotropic damage. The results are expressed as Equations (16) and (17):

$${\overline{\sigma }}_{kk}\left(t\right)=\frac{3K_e{\epsilon }_{kk}\left(t\right)+3{{\displaystyle \sum }}_{n=1}^N{K}_n{q}_{kk,n}^t}{1-\varnothing \left(t\right)}$$

(16)

$$q_{kk,n}^t={{\displaystyle \int }}_0^t{e}^{-\left(t-\zeta \right)/{\rho }_n}\frac{d{\epsilon }_{kk}\left(\zeta \right)}{d\zeta }d\zeta$$

(17)

where ${\overline{S}}_{\mathrm{ij}}\left(t\right)$ and ${\overline{\sigma }}_{kk}\left(t\right)$ are the effective deviatoric and volumetric stress tensors, respectively. $q_{ij,n}^t$ and $q_{kk,n}^t$ are the shear and volumetric hereditary integrals, respectively. The hereditary integrals are updated at the end of every converged time increment.

For implementation in the finite element method, the incremental shear and bulk stress tensors for cases of the isotropic damage are obtained by using Equations (14) and (16) as follows:

$$\begin{array}{cc}\hfill \Delta {\overline{S}}_{ij}\left(t\right)={\overline{S}}_{ij}\left(t\right)& -{\overline{S}}_{ij}\left(t-\Delta t\right)\hfill \\ & =\frac{2G_e{e}_{ij}\left(t\right)+2{{\displaystyle \sum }}_{n=1}^N{G}_n{q}_{ij,n}^t}{1-\varnothing \left(t\right)}-\frac{2G_e{e}_{ij}\left(t-\Delta t\right)+2{{\displaystyle \sum }}_{n=1}^N{G}_n{q}_{ij,n}^{t-\Delta t}}{1-\varnothing \left(t-\Delta t\right)}\hfill \end{array}$$

(18)

$$\begin{array}{cc}\hfill \Delta {\overline{\sigma }}_{kk}\left(t\right)={\overline{\sigma }}_{kk}\left(t\right)& -{\overline{\sigma }}_{kk}\left(t-\Delta t\right)\hfill \\ & =\frac{3K_e{\epsilon }_{kk}\left(t\right)+3{{\displaystyle \sum }}_{n=1}^N{K}_n{q}_{kk,n}^t}{1-\varnothing \left(t\right)}-\frac{3K_e{\epsilon }_{kk}\left(t-\Delta t\right)+3{{\displaystyle \sum }}_{n=1}^N{K}_n{q}_{kk,n}^{t-\Delta t}}{1-\varnothing \left(t-\Delta t\right)}\hfill \end{array}$$

(19)

It is assumed that the damage density $\varnothing \left(t\right)$ equals to $\varnothing \left(t-\Delta t\right)$ at the beginning of each time increment. Then, Equations (18) and (19) can be reduced to Equations (20) and (21), respectively, in Equation (20) to Equation (22):

$$\Delta {\overline{S}}_{ij}\left(t\right)=\frac{\left[2G_e+2{{\displaystyle \sum }}_{n=1}^N{G}_n\frac{{\rho }_n\left(1-e^{-\Delta t/{\rho }_n}\right)}{\Delta t}\right]\Delta e_{ij}\left(t\right)+2{{\displaystyle \sum }}_{n=1}^N{G}_n\left(e^{-\Delta t/{\rho }_n}-1\right)q_{ij,n}^{t-\Delta t}}{1-\varnothing \left(t-\Delta t\right)}$$

(20)

$$\Delta {\overline{\sigma }}_{kk}\left(t\right)=\frac{\left[3K_e+3{{\displaystyle \sum }}_{n=1}^N{K}_n\frac{{\rho }_n\left(1-e^{-\Delta t/{\rho }_n}\right)}{\Delta t}\right]\Delta {\epsilon }_{kk}\left(t\right)+3{{\displaystyle \sum }}_{n=1}^N{K}_n\left(e^{-\Delta t/{\rho }_n}-1\right)q_{kk,n}^{t-\Delta t}}{1-\varnothing \left(t-\Delta t\right)}$$

(21)

$$\Delta {\overline{\sigma }}_{ij}=\Delta {\overline{S}}_{ij}+\frac{1}{3}\Delta {\overline{\sigma }}_{kk}{\delta }_{ij}$$

(22)

3. Parameter Determination

3.1. Identification of Viscoelastic Model Parameters

Asphalt mastic comprises of fine aggregates and asphalt binder. The asphalt binder content in asphalt mastic is the same as the full HMA mixture (AC-20 modified asphalt mixture), excluding the binder absorbed by coarse aggregates (larger than 2.36 mm). The finished experimental beams cut from a cylindrical specimen have a dimension of 10 mm × 10 mm × 50 mm in length, width, and height, respectively. Uniaxial creep tensile tests were conducted to determine the creep compliance of asphalt mastic at a temperature of 20 °C. For the displacement-based linear viscoelastic constitutive model, an indirect method verified by our previous research [41] is applied to determine the fundamental relaxation modulus E₀ and ΔE(t) from the known compliance function. Then, all the material parameters in Equation (6) can be calculated. The results of Prony series coefficients are shown in

Table 1. Prony series coefficients of relaxation modulus for asphalt mastic.

i	ρn/s	En/MPa
1	104	0.1806
2	103	1.9356
3	102	8.959
4	101	56.961
5	100	72.7733
6	10−1	24.5185
7	10−2	124.834
8	10−3	434.302
E0 = 0.5972

. The Poisson’s ratio μ is assumed to be a constant of 0.4.

3.2. Identification of Damage Parameters

The Linear Amplitude Sweep (LAS) test has been approved by AASHTO as a standard method estimating damage tolerance of asphalt binders. The LAS test is a cyclic torsion test conducted in the dynamic shear rheometer (DSR) that uses increasing loading amplitudes to accelerate damage. In order to evaluate the damage of asphalt mastic, the same testing method in LAS test is applied.

The same sample used to determine the viscoelastic model parameters in Section 3.1 is prepared for the LAS test. Two types of testing are performed in succession. The first one is the frequency sweep designed to obtain information on the rheological properties, and the second one is the amplitude sweep intended to measure the damage characteristics of the material. The Kinexus lab+ rotational rheometer produced by Kinexus Rotational Rheometer (Produced by Malvern Panalytical Ltd., Royston, United Kingdom) is applied to conduct the two types of testing.

The frequency sweep test is performed at the selected temperature and applies oscillatory shear loading at constant amplitude over a range of loading frequencies to determine the damage parameter α in Equation (12). For this test method, the frequency sweep test employs an applied load of 0.01% strain over a range of frequencies from 0.01 Hz to 30 Hz at 20 °C. Data is sampled at the following 13 frequencies selected from the controller software: 0.01 Hz, 2.57 Hz, 5.14 Hz, 7.71 Hz, 10.27 Hz, 12.83 Hz, 15.39 Hz, 17.95 Hz, 20.52 Hz, 23.08 Hz, 25.66 Hz, 28.22 Hz, and 30.78 Hz. Complex shear modulus |G*| and phase angle δ are recorded at each frequency as shown in Figure 2.

In order to determine the damage parameter α in Equation (12) from frequency sweep test, data for the dynamic modulus|G∗|(ω) and phase angle δ(ω)for each frequency is converted to storage modulus G’(ω) by Equation (10). Then, A best-fit straight line is applied to a plot with log ω on the horizontal axis and log G’ (ω) on the vertical axis using Equation (11) as shown in Figure 3. Thus, the value obtained for m is recorded and the value of α is obtained by Equation (12). In this test, the value of α is: α = 1/m = 1/0.3319 = 3.013.

The second test is run at temperature 20 °C using oscillatory shear in strain control mode at a frequency of 5 Hz. The loading scheme consists of a continuous oscillatory strain sweep. Strain is increased linearly from 0 to 1% over the course of 25,000 cycles of loading for a total test time of 5000 s. Peak shear strain and peak shear stress are recorded every 5 load cycles (1 s), along with phase angle δ and complex shear modulus |G*|.

Figure 4 shows the complex shear stress and shear modulus extracted from amplitude sweep test for asphalt mastic at temperature 20 °C. The stress–strain diagram shows that the complex shear stress grows with a higher rate until the peak point as the strain increases, where the peak stress occurs at 0.7% strain level. And then the stress decreases after the stress–strain peak point. The shear modulus versus strain diagram shows that the shear modulus decreases with a higher rate until 0.7% strain level. After this region the shear modulus decreases with a constant slope.

In the pre-peak region, micro-cracking evolves and accumulates as the strain increases. In the post-peak region, the micro-cracking leads to a separation or debonding after the stress–strain peak point. In linear viscodamage theory domain, the failure behavior is not considered. Therefore, we assume that the maximum value of the rate of the damage density occurs at the peak point of the stress–strain diagram and considered the inflection point of stress-strain diagram as the strain at which the maximum nominal stress occurs.

The normalized damage density versus strain curve is produced as shown in Figure 5 using Equation (13). The best fitting equation is obtained in Equation (23):

$$\varnothing \left(t\right)=1.5571{\epsilon }_t^2+0.4213{\epsilon }_t-0.0265$$

(23)

Thus, for a given strain at the current incremental step, the rate of the damage density can be calculated, given as Equation (24):

$$\stackrel{•}{\varnothing }\left(t\right)=3.1142{\epsilon }_t+0.4213$$

(24)

Thus, the damage evolution law for asphalt mastic can be calculated using Equations (23) and (24) when the current strain is given. For the initial increment, the rate of the damage density equals to 0.4213.

4. Microstructural FE Model Construction

An asphalt mixture is a particulate composite material consisting of aggregate, mastic, and air voids. In this section, a cylindrical HMA mixture (AC-20) lab specimen was prepared for capturing the internal microstructure with the non-destructive industrial X-ray CT technique. Then, the grayscale thresholds for dividing aggregate, matrix, and air voids were chosen based on the OTSU method. Furthermore, an improved morphological multiscale algorithm was applied to segment the coarse aggregate adhesion images. Additionally, a pixel-based 2D image model of the specimen was constructed. In the following sub-sections, the reconstructed process of the microstructure model is described in details.

4.1. Microstructure Acquisitions

X-ray CT is an effective method to capture the internal microstructure of asphalt mixture. Industrial CT systems consist of the following basic components: a radiation source, radiation detector, collimator, data-collecting system, mechanical specimen-scanning system, computer system (hardware and software), auxiliary system (auxiliary power supply and auxiliary security system), and so on. In this study, the YXLON Compact-225 CT X-ray scanner is used to obtain the detailed microscopic structure of the asphalt concrete specimen.

4.2. Digital Image Processing

Scanned CT images of asphalt mixture have different grayscale intensities between 0 and 255, where denser materials have higher intensity. The pixel grayscale intensities of aggregate are close to 255. The pixel grayscale intensities of air voids and background are close to 0. The pixel grayscale intensities of the matrix are between that of aggregates and air voids. The OTSU method was utilized for segmentation of the AC micro-structure. Figure 6 shows the phase-segmented AC image.

Particles that are connected each other is a common phenomenon in the segmentation process of X-ray CT slices which will increase the difficulty of numerical modeling [42,43,44] since the finite element code will treat the connected particles as an isolated one. To overcome this difficulty, an improved morphological multiscale algorithm based on several different sizes of structural elements was introduced to segment the coarse aggregate adhesion images. Details of this segmentation algorithm can be found in our previous research [45]. The segmented image can be seen in Figure 7.

4.3. Digital Samples Generation

A binary image is represented by an m × n logical matrix where pixel values are 1 (true) or 0 (false). In order to input the element and node information into the finite element software, such as ABAQUS, the node and element numbering rules for generating the pixel-based numerical model are defined as follows: the node and element numbering sequences start from the lower left corner of the matrix, and then go to the right side of each line, the last number of each line is followed by the next line. The node and element numbering diagram of the CT image pixel matrix is shown in Figure 8. The element and node information is generated with the MATLAB programming and written into an input file of ABAQUS for numerical simulations. The four-node bilinear reduced integration with hourglass control elements (CPS4R) with a unit length were used in constructing the mesh. The aggregates contained a total of 318,918 elements, the mastic phase had 143,654 elements, and the remaining elements were part of the air voids inclusions with a total number of 1,209,643 elements. The air void was excluded in this digital sample to reduce the number of elements. Consolidation constraint boundary conditions were applied to the bottom of the asphalt mixture digital sample and the displacement-based loading was imposed to the top of the digital specimen. Figure 9 presents the pixel-based numerical model of asphalt mixture.

5. Numerical Simulation and Discussion

The purpose of this section is to conduct the micromechanical simulations of the digital sample reconstructed from X-ray CT slice to predict the viscoelastic responses with/without damage behavior. The aggregates are considered as a linear elastic material and the modulus of elasticity and Poisson’s ratio for the aggregates are assumed to be 25 GPa and 0.25, respectively. The mastic phase is considered as a linear viscoelastic material. Parameters for the viscoelastic constitutive and damage of the mastic phase are determined in Section 3.1 and Section 3.2. The numerical viscoelastic constitutive model is implemented within the FE code using FORTRAN. The ABAQUS user material subroutine (UMAT) is applied for this purpose.

In order to investigate the viscoelastic behavior of asphalt mixture with and without damage, four different compressive strain levels (10 μɛ, 20 μɛ, 50 μɛ, and 70 μɛ) were applied to the 2D digital sample generated in Section 5.2 at 20 °C.

5.1. The von Mises Stress in the Whole Specimen

von Mises Stress is a yield criterion and is commonly called Mises equivalent stress, which follows the fourth strength theory of material mechanics. The von Mises stress statistical distribution of asphalt mastic for three time increments (step 1, step 100, and step 400) in the whole specimen under different loading levels is presented in Figure 10.

As shown in Figure 10, the von Mises stress curves with and without damage for different step in the whole specimen have the same trend. Comparisons between the damaged and undamaged von Mises stress curves in step 1, step 100, and step 400 for different compressive strain levels shows that the damaged simulation results have a larger stress distribution. This might be caused by the damage density which is less than 1 in Equation (20) leading to a higher stress level than that of no damage. As expected, the von Mises stress distribution is smaller as the loading time increases due to the viscoelastic behavior of asphalt mastic. All this indicates that the microscopic cracking characteristic of asphalt mastic can be demonstrated by the damage-coupled viscoelastic model used in this paper.

5.2. The von Mises Stress of Three Selected Locations

Three critical locations in the mastic phase were selected for the von Mises stress analysis. Location 1 and location 3 are in the area where two aggregates are very close to each other while location 2 is in the region with a larger mastic distribution. All these selected locations can be seen in Figure 11.

The von Mises stress vs. loading time within these three locations were obtained as shown in Figure 12. It is found that the von Mises stress curves with and without damage under various loading levels show a significantly consistency in different locations. The von Mises stress at location 1 and location 3 shows a higher magnitude than that at location 2. This is because the regions of location 1 and location 3 are located in between two irregular coarse aggregates which shall cause stress concentration. This is also why the tiny cracks would most likely appear at these positions in asphalt mixture pavement.

As indicated in Figure 12, the damaged von Mises stress at three selected locations show a higher magnitude than that of undamaged von Mises stress due to the damage density. The variation increases obviously with the increasing of compressive loading strain levels and the variation decreases as the loading time increases.

Since Figure 12 cannot clearly demonstrate differences among the four compressive strain levels (10 μɛ, 20 μɛ, 50 μɛ, and 70 μɛ), stress variation ratios were used to quantify their differences. The stress variation ratio for each compressive strain level is defined as relative change in every incremental time step:

$$variation ratio=\frac{the damage stress-the undamaged stress}{the undamaged stress}$$

The calculated variation ratios vs. incremental time step were plotted in Figure 13. From the graph, it is observed that the stress variation ratio increases faster in the first 40 s. This is because the initial stress would mainly release in this period. The variation ratio increases obviously with the increasing of compressive loading strain levels due to the strain-based damage density. This finding is also consistent with the finding that from analysis of stress variation under different loading strain. However, the variation ratio may have very different curve shapes when the loading level is high. This is mainly due to the fact that the asphalt mixture is a particulate composite material. This may cause the inhomogeneous stress field when the load is applied to the specimen which leads to different damage density at every Gauss integral point.

6. Conclusions

Development of structural mechanistic models that can effectively predict the performance of pavements under realistic loading and boundary conditions is highly desirable. Micromechanical models can predict material properties based upon the properties of the individual constituents such as the asphalt mastic and aggregate. These structural models need to incorporate constitutive material models that describe the behavior of asphalt mixtures. Thus, in the aim of optimizing materials design, it is necessary to study the micromechanical damage response of this composite material.

In this paper, a damage-coupled viscoelastic model of asphalt mastic was proposed based on the LAS test. We derived the displacement-based incremental formulations of the viscoelastic constitutive models for asphalt mastic with some algebraic manipulations and then implemented the incremental formulations in a finite element code. The viscoelastic parameters in the constitutive model were determined by applying the uniaxial tensile tests and torsion tests. Then, the LAS test was applied to obtain the damage parameters in the constitutive model. Furthermore, the nondestructive industrial X-ray CT technique was employed to capture the internal microstructure of asphalt mixture. The X-ray CT image was processed using segmentation techniques to better identify and model the mixture components. A pixel-based 2D image digital model consisting of coarse aggregate and asphalt mastic was constructed. Finally, the image-based FE model incorporated with damage-coupled viscoelastic asphalt mastic phase and elastic aggregates was used for the compressive test simulations successfully in this study.

The improved morphological multiscale algorithm presented in this study can better segment the coarse aggregate adhesion images and the image-based FE mode can provide geometric information for the aggregates, air voids, and asphalt mastic. The compressive test simulation results show that the damaged simulation results have a larger stress distribution compared with the undamaged simulation due to the irregularity of the coarse aggregates. The von Mises stress distribution is smaller as the loading time increases due to the viscoelastic behavior of asphalt mastic. It can also provide insight on the damaged mechanisms and the possible location in asphalt mixture where microscopic cracking would most likely occur.

It can be concluded that the presented damage-coupled viscoelastic constitutive and the image-based microstructure digital sample are capable of describing the damaged response of asphalt mixture at a constant temperature. It is feasible to utilize the proposed constitutive and digital sample to conduct numerical simulations for damaged mechanical responses of asphalt mixtures. Our future work will focus on extending this approach to conduct three-dimensional digital sample tests of microstructural asphalt mixture under complex loading path, as well as further generalizing the high-resolution numerical model including aggregate-mastic interface.