Abstract
This paper presents an image-based micromechanical modeling approach to predict the viscoelastic behavior of asphalt mixture. An improved image analysis technique based on the OTSU thresholding operation was employed to reduce the beam hardening effect in X-ray CT images. We developed a voxel-based 3D digital reconstruction model of asphalt mixture with the CT images after being processed. In this 3D model, the aggregate phase and air void were considered as elastic materials while the asphalt mastic phase was considered as linear viscoelastic material. The viscoelastic constitutive model of asphalt mastic was implemented in a finite element code using the ABAQUS user material subroutine (UMAT). An experimental procedure for determining the parameters of the viscoelastic constitutive model at a given temperature was proposed. To examine the capability of the model and the accuracy of the parameter, comparisons between the numerical predictions and the observed laboratory results of bending and compression tests were conducted. Finally, the verified digital sample of asphalt mixture was used to predict the asphalt mixture viscoelastic behavior under dynamic loading and creep-recovery loading. Simulation results showed that the presented image-based digital sample may be appropriate for predicting the mechanical behavior of asphalt mixture when all the mechanical properties for different phases became available.
1. Introduction
Asphalt mixture is the most widely used road construction material in the paving industry. The mechanical behavior of the asphalt mixture is complex due to the heterogeneity of the asphalt composite material and the time-dependent viscoelastic binder. In the aim of optimizing materials design, it is necessary to study the behavior of this composite material. The traditional trial and error approach in industrial practice has focused on empirical experiments that develop correlations between the macro phenomena and the material characteristics. However, this traditional approach is cost-intensive, high-resource-consuming, and time-consuming. The developed correlations are not sufficient to gain insight into the performance of asphalt mixture. To overcome these difficulties, numerical techniques are developed to gain some insight into the problems.
Many numerical attempts have been made to investigate the internal behavior of asphalt mixture or to estimate the test indices that can represent the macro phenomena under certain conditions [1–5]. But these researches consider the asphalt mixture as being homogenous material with equivalent properties. When using these homogenous numerical models, the influences of the microstructure such as angularity of the aggregates, percentage of the air voids, and distributions of the asphalt mastic are not taken into account. The mechanical behavior of asphalt mixture strongly depends on the particular components including aggregates, asphalt mastic, and air voids. In order to gain accurate predictions of the mechanical behavior for asphalt mixture, particular attention has to be paid to the study of its microstructure and to the mechanical behaviors of its different phases.
A number of researchers have paid attention to employing the discrete element model (DEM) in studying the micromechanical behaviors of asphalt mixtures. Chang and Meegoda [6] used a micromechanical model based on the discrete element method (DEM) to investigate the viscoelastic behavior of the internal structure of HMA. Collop et al. [7, 8] applied the discrete element method (DEM) to generate a three-dimensional (3D) idealized mixture model for simulating the behavior of a highly idealized asphalt mixture under uniaxial and triaxial compressive creep tests. Wu et al. [9] used the discrete element method (DEM) to simulate constant strain rate compressive tests for an idealized asphalt mixture comprising approximately single-sized sand particles. Cai et al. [10, 11] used the discrete element method (DEM) to simulate the uniaxial compression tests of a realistic asphalt mixture under constant strain rate. In these studies, the produced numerical samples are either highly idealized model or randomly generated model. The most important advantages of these discrete element samples are the ability to handle complex and changing contact geometries and large displacements. However, these discrete element samples were still not adequate enough to describe the complicated microstructure of asphalt mixture with highly irregular shaped components. Another disadvantage of these discrete element samples is the relatively high computational cost.
To overcome this shortcoming, many researchers developed micromechanical models to investigate the behavior of asphalt mixture by using the finite element method (FEM). Dai [12, 13] developed a two-dimensional (2D) microstructure finite element model of the asphalt mixture generated from the scanned image to predict the viscoelastic properties of the asphalt mixtures, where elastic properties were assigned to aggregates and viscoelastic properties were assigned to mastic phase. Kim and Lutif [14] proposed a computational micromechanics modeling approach to predict damage-dependent constitutive behavior of asphalt mixtures. Arshadi and Bahia [15] developed the two-dimensional (2D) images of mastic and mortar scales artificially and used them to characterize the properties of those scales. And then, the 2D scanned images of asphalt mixtures are utilized to study the asphalt mixture behavior under uniaxial creep and recovery loading. But the majority of the image-based microstructure finite element models is limited to two-dimensional samples. Lacking of three-dimensional information in FE framework makes them less able to represent complicated morphology of asphalt mixture when used for simulating.
The nondestructive X-ray CT technique is used as an effective method to capture images of the internal microstructure of asphalt mixture. There have been a number of recent attempts to use X-ray CT images to reconstruct a three-dimensional (3D) finite element specimen to predict the mechanical properties of asphalt mixture [16–19]. The image-based finite element method requires the image analysis techniques in order to segment the CT slices consisting of different compositions, the mechanical properties of different phases, and the use of mechanical testing to identify the material parameters. The majority of these microstructure numerical samples for asphalt mixture faces challenges with regard to accurate geometry of different phases and precise description of asphalt mastic behavior.
This paper presents an image-based micromechanical modeling approach for simulating the viscoelastic response of asphalt mixture at a given temperature. In order to enhance the segmentation of different phases in the CT image, we developed an improved image analysis technique based on the OTSU thresholding operation for use. Then, the three-dimensional viscoelastic constitutive model is applied to better describe the mechanical behavior of asphalt mastic and implemented in a finite element code using the ABAQUS user material subroutine (UMAT). After that, we propose an experimental procedure for identifying and verifying the viscoelastic parameters at a given temperature in the linear viscoelasticity range. 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. Microstructural Reconstruction Based on X-Ray CT Images
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 nondestructive industrial X-ray CT technique. Then, the grayscale thresholds for dividing aggregate, matrix, and air voids were chosen based on the improved OTSU method. Additionally, a pixel-based 3D image model of the specimen was constructed. In Sections 2.1 and 2.2, the reconstructed process of the microstructure model is described in detail.
2.1. Microstructure Acquisitions and Digital Image Processing
X-ray CT is a commonly used nondestructive technique for characterizing the internal microstructure of asphalt mixture. The basic image-forming principle of the industrial X-ray CT is that suppose the X-ray attenuation coefficient of the ideal material is . If the original X-ray strength is , the attenuated X-ray strength is . For homogeneous material, the attenuated X-ray strength follows the Bill Exponential Law described below:
For heterogeneous materials, the attenuated X-ray strength can be calculated by dividing the object into small units. So we can get the following equation:
The original X-ray strength and the attenuated X-ray strength can be easily measured and the right side of (2) can be calculated by numerical method. In this study, the YXLON Compact-225 CT X-ray scanner is used to obtain the detailed microscopic structure of the asphalt mixture specimen.
Generally, asphalt mixture is a three-phase composite material including aggregate, matrix, and air voids. Scanned CT images of asphalt mixture have different grayscale intensities between 0 and 255, where denser materials have higher intensity. The OTSU algorithm is a popular thresholding method in adaptive optimal threshold selection for image segmentation.
Beam hardening effect is a common phenomenon in X-ray CT images where the edge is brighter than the center as presented in Figure 1(a). Figure 1(b) shows the phase-segmented HMA mixture (AC-20) image using the OTSU method and it can be noticed that the OTSU thresholding operation only captures the aggregates at the edges of the image. To reduce the beam hardening effect in the CT image, an improved OTSU method developed by Liu and Li [20] is used for this purpose. In this method, the CT image is decomposed into a series of circular subimages with 50% overlap between each subimage given in Figure 1(c). Then, the OTSU thresholding operation is applied for each circular subimage to segment the three phases. Figure 1(d) shows the aggregate phase segmented from the original CT image using the improved OTSU method and the phase-segmented image indicated with different colors can be seen in Figure 1(e).
(a)
(b)
(c)
(d)
(e)
2.2. Numerical Model Generation
A binary image is represented by an logical matrix where pixel values are 1 (true) or 0 (false). The voxel defined in Figure 2 is generated by expanding the pixel into three-dimensional space (3D). A voxel-based 3D digital reconstruction model of asphalt mixture is constructed when every pixel in the consecutive CT images is converted into voxel as shown in Figures 3(a) and 3(b).
(a)
(b)
(c)
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 voxel-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, and the last number of each line is followed by the next line; the position of corresponding element number in the first image differs from that of the adjacent image by , while the position of corresponding node number in the first image differs from that of the adjacent image by . The element and node information is generated with the MATLAB programming and written into an input file of BAQUS for numerical simulations. Figure 3(c) presents the voxel-based numerical model of asphalt mixture.
3. Modeling of Asphalt Mastic
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. The challenge in the modeling of micromechanical finite element model for asphalt concrete includes the time, temperature, and rate-dependent behavior of the asphalt mastic [16]. In this section, the three-dimensional viscoelastic constitutive model is used to represent the behavior of asphalt mastic phase. Then, the incremental formulations of the constitutive model implemented in a finite element code are developed in details.
3.1. Viscoelastic Constitutive Model
For linear viscoelastic materials, including asphalt mastic, the stress–strain constitutive relation is expressed by convolution integrals. In the case of strain response at constant stress, the convolution relations is explained as follows for one-dimensional problems:where the is the instantaneous elastic compliance, is the reduced time, and is the transient compliance. It is given bywhere is the th coefficient of the Prony series and is the th retardation time.
For stress response at constant strain, the convolution relations can be represented as follows:where the and are the instantaneous elastic modulus and the transient modulus, respectively; is the th coefficient of the Prony series and is the th relaxation time.
For three-dimensional problems, (3) can be decomposed into deviatoric and volumetric components, such thatwhere and are the deviatoric strain and volumetric strain, respectively; is the total strain and is the Kronecker delta; and are the instantaneous effective elastic shear and bulk compliances, respectively; and are the transient shear compliance and bulk compliance, respectively.
Similarly, (5) can be decomposed into deviatoric and volumetric components, such thatwhere and are the deviatoric stress and volumetric stress, respectively; , , and are the total stress, the instantaneous effective elastic shear modulus, and bulk modulus, respectively; and are the transient shear modulus and bulk modulus, respectively.
3.2. Numerical Implementation for the Constitutive Model
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 three-dimensional incremental deviatoric and volumetric formulations can be derived with some algebraic manipulations.
3.2.1. Stress-Based Incremental Formulations
For stress-based incremental deviatoric and volumetric formulations, the results are expressed aswhere and are the incremental shear and bulk strains, respectively; the variables and are the shear and volumetric hereditary integrals for every Prony series term at previous time , respectively; the hereditary integrals are updated at the end of every converged time increment; for the initial increment, the variables and are set to and , respectively; is the total incremental strain.
3.2.2. Displacement-Based Incremental Formulations
Obviously, a similar derivation procedure may be carried out for the case of stress response at constant strain. The incremental formulations are given bywhere and are the incremental shear and bulk stress, respectively; for the initial increment, the variables and are set to and , respectively; is the total incremental stress.
The stress-based and displacement-based three-dimensional numerical constitutive models are implemented within the FE code using FORTRAN language. The ABAQUS user material subroutine (UMAT) is applied for this purpose.
4. Identification and Verification
4.1. Identification of the Material Parameters
According to the elastic–viscoelastic correspondence principle, each of the tensile creep compliance, the shear compliance, and the bulk compliance can be obtained from the other two using Laplace transform. The tensile creep compliance and the shear compliance can be easily determined directly from uniaxial tensile tests and torsion tests, respectively. Then, the bulk compliance can be calculated from the tensile creep compliance and the shear compliance . A similar set of equations may be formulated for the bulk modulus when the tensile modulus and the shear compliance are available. For displacement-based linear viscoelastic constitutive model, an indirect method verified by previous research [21] is applied to determine the fundamental relaxation modulus and from the known compliance function.
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), 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 tensile tests and torsion tests were applied to determine the material parameters in three-dimensional viscoelastic constitutive model at a temperature of 20°C. In the theory of linear viscoelasticity, the strain response to any applied stress is independent of the stress magnitude. This characteristic can be adopted to the static creep test by monitoring the creep compliance as stress increases. If the creep compliance curves vary a little subjected to a range of loading levels, linear viscoelasticity holds. The tensile creep data and shear creep data obtained from uniaxial tensile tests and torsion tests, respectively, were employed in determining Prony series coefficients of creep compliances. Then, the Prony-based series coefficients of relaxation modulus can be converted by the method proposed in previous study [21]. The results of Prony series coefficients are shown in Table 1.
4.2. Verifications of the Incremental Constitutive Model
The capability of the numerical constitutive model can be examined by comparing the numerical predictions with the observed laboratory tests. In this paper, two basic loading paths (bending loading and compression loading) were utilized to conduct a series of displacement-based tests at a temperature of 20°C to verify the capability of the incremental constitutive model.
Beams with a dimension of 30 mm × 30 mm × 120 mm in length, width, and height were applied to conduct bending flexural tests. Two primary modes of displacement-based loading bending tests were performed on the beams. For the verification purpose, three-dimensional numerical constitutive model mentioned above was introduced to reproduce the bending tests. The comparisons between the predictions and the corresponding experimental results are given in Figure 4.
(a)
(b)
From Figure 4, it can be seen that the numerical predictions exhibit consistent trend with the experimental results. The differences of the vertical reaction force between the numerical predictions and measured results for the two primary modes of displacement-based loading are 2.1% and 2.2%, respectively. Obviously, a rather good agreement between the simulations and experimental measurements can be obtained.
As for compression loading, cylindrical asphalt mastic specimens having a diameter and height of 100 × 100 mm were utilized for testing. The same shapes of loading were used for compression responses to arbitrary displacement-based loading step history of mastic asphalt. The experimental results and the testing results are given in Figure 5.
(a)
(b)
It can be observed that the predicted results obtained from the numerical constitutive model show good agreement with the corresponding experimental data. This suggests that once the appropriate viscoelastic constitutive parameters for asphalt mastic at a constant temperature are available, the numerical constitutive model is capable of describing the material response under compression loading.
5. Numerical Simulations
The purpose of this section is to conduct the micromechanical simulations of the digital sample reconstructed from X-ray CT slices to predict the viscoelastic properties under dynamic loading modes. Asphalt mixture specimen with a dimension of 100 mm diameter by 200 mm height is usually recommended for the dynamic modulus tests and the creep tests to minimize edge effects (AASHTO TP 62-03, AASHTO TP-70). Usually, the thickness of a typical asphalt surface layer is far less than 150 mm. The indirect tensile test (IDT) used on field cores is the most effective method to evaluate the mechanical properties of the existing pavement. In this study, an asphalt mixture specimen with a dimension of 100 mm diameter by 30 mm height was utilized for the numerical simulations. The specimen was scanned by the X-ray CT device with a resolution of 1024 × 1240. All the CT slices were converted to low-resolution images to reduce the number of elements. The numerical model of the asphalt mixture was generated from the image pixels using the MATLAB code as shown in Figure 3(c).
The eight-node 3D solid integration elements (C3D8) with a unit thickness were used in constructing the mesh. The aggregates contained a total of 35,611 elements, the matrix phase had 44,181 elements, and the remaining elements were part of the air voids inclusions with a total number of 179 elements. The aggregates and the steel loading strip are considered as a linear elastic material and the modulus of elasticity and Poisson’s ratio for the aggregate and the steel loading strip are assumed to be 25 GPa and 0.25, 80 GPa and 0.2, respectively. The air void will be included in this digital sample with an elastic modulus of 0.5 MPa and Poisson’s ratio of 0.3. The matrix is considered as a linear viscoelastic material. Parameters for the viscoelastic constitutive of the matrix are determined in Section 4.1 and the numerical viscoelastic constitutive model is implemented within the FE code using FORTRAN. The ABAQUS user material subroutine (UMAT) is applied for this purpose.
5.1. The Dynamic Test
In the Mechanistic Empirical Pavement Design Guide (MEPDG), the dynamic modulus of asphalt mixture is used as an important input parameter to characterize the temperature and frequency dependent behavior for pavement design and construction [22]. The standard dynamic modulus determination procedure consists of the uniaxial partial sinusoidal compressive test and the indirect partial sinusoidal tension test, while the stress state of the indirect tension test specimen subjected to vertical loading is very similar to that of the field. It is apparent that dynamic modulus measured from the indirect tension test can better characterize the in-situ behavior of asphalt mixture.
The dynamic modulus simulations in indirect tension mode under six partial sinusoidal cycle loading frequencies (0.1, 0.5, 1, 5, 10, and 25 Hz) were conducted to show the utility of the developed microstructure digital model of asphalt mixture containing the three main phases. The displacement-based loading was imposed to the top of the steel loading strip, and the steel loading strip distributed the applied load on the top surface of the numerical model. In order to ensure that the asphalt mixture behaves as a linear viscoelastic material, the vertical strain was confined to 0.01%. The displacement-based incremental constitutive model for asphalt mastic derived in Section 3.2.2 would be incorporated in ABAQUS user material subroutine (UMAT) to model the effective asphalt mastic behavior.
For frequencies of 0.1, 0.5, and 1 Hz, six displacement-based loading cycles were used and ten displacement loading cycles were applied for the frequencies of 5, 10, and 25 Hz. To balance the computational cost and the smoothness of the stress/strain-time response curves, 20 time increments (computation points) were applied for each loading cycle. In the dynamic modulus simulation, the reacting force and the vertical displacement of the whole model were recorded. Stress is defined as the reacting force divided by the cross section area of the digital specimen. Strain is the change in the vertical deformation of the digital specimen divided by the initial height. The dynamic modulus is calculated by dividing the peak stress amplitude with the peak strain amplitude:where is the dynamic modulus, is the peak-to-peak stress amplitude, and is the peak-to-peak strain amplitude.
In this study, the averages of the reacting force and the vertical displacement of the last two loading cycles were used to calculate the dynamic modulus. The stress/strain-time responds under different loading frequencies were presented in Figure 6.
(a) Frequency = 0.1 Hz
(b) Frequency = 0.5 Hz
(c) Frequency = 1 Hz
(d) Frequency = 5 Hz
(e) Frequency = 10 Hz
(f) Frequency = 25 Hz
It is clearly shown that the overall stress responds under different loading frequencies drop as the loading time increases, and the peak strain lags behind the peak loading under the same loading cycle, which demonstrates that the digital sample behaves as a viscoelastic material under cycle loading. The dynamic modulus was calculated with (11) as shown in Figure 7. It shows that the dynamic modulus increases with the loading frequencies. The results indicate that the developed digital microstructure sample is capable to predict the macroscopic viscoelastic behaviors of asphalt mixture with experimental material parameters.
5.2. Repeated Creep-Recovery Test
A typical vehicle load in flexible pavement is a cyclic load with loading and unloading periods. So in the laboratory tests, repeated creep-recovery test is more representative to what the asphalt mixture experiences under traffic loading. Figure 8 shows the stress input history for the repeated creep-recovery test simulation with a 6-second loading duration and a 2-second unloading duration.
Repeated creep-recovery test simulation is conducted using the same microstructure digital model of asphalt mixture utilized for dynamic simulations. The vertical stress applied to the top of the steel loading strip is 1.2 MPa due to the high stiffness for AC-20 at a temperature of 20°C. To balance the computational cost and the smoothness of the strain-time response curves, 40 time increments (computation points) were applied for each loading cycle. The force-based incremental constitutive model for matrix derived in Section 3.2.1 would be incorporated in ABAQUS user material subroutine (UMAT) to model the effective asphalt mixture behavior.
Figure 9 shows the vertical displacement and strain contour at 60th time step, respectively. It is found in Figure 9(a) that the vertical displacement decreases from the top to the bottom but not left-right symmetric due to the heterogeneity of the asphalt mixture digital model. Figure 9(b) presents that the local vertical strain distributions are largely in the state of compression. The vertical compressive strain mainly occurs in the vicinity of the vertical central axis. This is because the aggregates around the central axis will be extruded into a skeleton to withstand the compressive loading. It also can be seen from Figure 9(b) that the matrix phase suffers greater deformation due to the large difference in stiffness between aggregate and matrix.
(a)
(b)
The corresponding strain-time response under vertical stress for 64 time steps is presented in Figure 10. It can clearly be noticed that the primary and secondary creep regions can be observed from this figure. However, damage was not considered in this study, and the tertiary creep region does not occur as the numbers of loading cycle increases. This trend is in close agreement with experimental data regarding the linear viscoelastic behavior of asphalt mixture under repeated loading in the previous studies [23, 24]. It is clear that the three-dimensional (3D) image-based numerical model and the verified incremental constitutive model can be used to predict the mechanical response of asphalt mixture under repeated creep loading.
6. Conclusions
In this paper, the nondestructive industrial X-ray CT technique was employed to capture the internal microstructure of asphalt mixture. The CT images were segmented with the improved OTSU method to reduce the beam hardening effect which is a common phenomenon in X-ray CT images. The image-based three-dimensional (3D) FE model including aggregates, asphalt mastic, and air voids was developed for numerical predictions. In this 3D model, the aggregate phase and air void were considered as elastic materials while the asphalt mastic phase was considered as linear viscoelastic material. We derived the displacement-based and stress-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 material parameters in the three-dimensional viscoelastic constitutive model were determined by applying the uniaxial tensile tests and torsion tests. Furthermore, the capability of the incremental constitutive model was examined by comparing the numerical predictions with the observed laboratory tests at a temperature of 20°C. Finally, the image-based three-dimensional (3D) FE mode incorporated with viscoelastic asphalt mastic phase and elastic aggregates and air voids was used for the dynamic and the repeated creep-recovery test simulations.
The improved OTSU method presented in this study can better identify the CT images segmentation and the developed image-based three-dimensional (3D) FE model can provide geometric information for the aggregates, air voids, and asphalt mastic. Comparison between the numerical predictions and the experimental results shows that the incremental constitutive model has considerable promise for predicting the asphalt mastic mechanical properties under two basic loading paths. It can be concluded that once the appropriate viscoelastic constitutive parameters for asphalt mastic at a constant temperature are available, the numerical constitutive model is capable of describing the material response. Simulation results of the dynamic test and the creep test presented in this study showed that the 3D finite element models provided reasonable predictions of the complex modulus and creep characteristics. It is feasible to utilize the proposed approach to conduct numerical simulations for mechanical responses of the asphalt mixture. Our future work will focus on extending this approach to conduct triaxial tests of microstructural asphalt mixture under complex loading path, as well as further generalizing the high-resolution numerical model including aggregate-mastic interface.
Competing Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
This publication is supported by Project (51508109), research on dynamic mechanical behavior of asphalt mixture and key parameters prediction based on energy method sponsored by National Natural Science Foundation of China. The authors thank all those who contributed in the experimental part of this study.