Abstract
Asphalt mixture is a particulate composite material consisting of aggregate, mastic, and air voids. The computed tomography (CT) image-based finite element approach is used as an effective method to simulate micromechanical response of asphalt mixture. For finite element analysis, the accuracy of the finite results is determined by the size of the finite element. In this paper, a voxel-based three-dimensional (3D) digital reconstruction model of asphalt mixture with the CT images after being processed was proposed. In this 3D model, the aggregate phase was considered as elastic materials while the asphalt mastic phase was considered as linear viscoelastic material. Four micromechanical digital models were generated, whose voxel sizes were 0.5 mm, 0.67 mm, 1.0 mm, and 2.0 mm, respectively. The four digital models were used to conduct uniaxial creep test for predicting creep stiffness modulus to investigate the effect of voxel size. Simulation results showed that the voxel sizes had a significant effect on creep stiffness modulus. For the creep simulation test, the most appropriate voxel size whose creep stiffness modulus changes within 2.5% is 1.0 mm with regard to time steps, computational time, aggregate, and mastic shape representations.
1. Introduction
1.1. Background
Asphalt mixture is the most widely used road construction material in the paving industry. 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 does not provide information that can be used to explain why certain mixtures perform better than others. And also this approach is cost intensive, high-resource consuming, and time consuming. To overcome these shortcomings, some microstructural methods and numerical techniques are developed to gain insight into the performance of asphalt mixture.
The composition of asphalt mixture can be considered to be included in 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 mechanical behavior of asphalt mixture is complex due to the heterogeneity of the asphalt composite material. Several research studies have investigated the internal properties and performance of asphalt mixture at microstructural scale using the nondestructive X-ray CT technique [1–9]. There have also 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 [10–14]. The majority of these microstructure numerical samples for asphalt mixture face challenges with regard to accurate geometry of different phases. The geometry accuracy and volume composition of aggregate and mastic phases are directly affected by the resolutions of CT image for asphalt mixture. But most of the previous image-based finite element approaches are limited to low precision numerical models in order to reduce the time steps and computational time [15–17]. However, it is not sufficient enough to consider the quantitative effects of finite element dimensions on asphalt mixture simulation results.
1.2. Objective
The objective of this study is to investigate size effects of finite element model for three-dimensional microstructural modeling of asphalt mixture. Four image-based microstructure models with different finite element dimensions for asphalt mixture were established in order to quantify the size effects of finite element model for three-dimensional microstructural modeling. Firstly, the nondestructive industrial X-ray CT technique was used to capture the internal microstructure of the asphalt mixture specimen. The grayscale thresholds for dividing aggregate, asphalt mastic, and air voids were chosen based on the Otsu method. A morphological multiscale algorithm was applied to segment the coarse aggregate adhesion images. Then, four voxel-based 3D digital reconstruction models of the asphalt mixture specimen is constructed based on the X-ray CT images. After that, the parameters of the asphalt mastic phase are obtained using the dynamic shear rheometer (DSR) test. Finally, the image-based numerical specimens are applied to conduct the creep simulations in order to investigate size effects of finite elements on creep stiffness.
2. Three-Dimensional Microstructural Reconstruction
2.1. CT Image Processing
At present, there is no image analysis software especially for material segmentation of CT slice images of asphalt mixture. The usual method involves segmentation of the internal materials according to the gray distribution characteristics of asphalt mixture. The most widely used approach is the Otsu method [18–20].
Because the difference of material density is not very obvious in the areas where the particles are in close contact with each other, cohesion of aggregate particles occurs when using the Otsu method to segment the aggregate and asphalt mastic phases. So, the adhesive aggregates are identified as individual particles by computer in the subsequent mechanical simulation, resulting in a calculation result which is inconsistent with the actual situation. Therefore, before constructing the numerical model, adhesive aggregates in the images must be separated. At present, coarse aggregate adhesive separation is mainly based on manual processing, which is convenient for single image processing [21–23]. However, it is not conducive for rapid and efficient image processing if the number of images is too large. Therefore, an automatic separation algorithm for adhesive aggregates is necessary.
This section presents a morphological multiscale separation algorithm for adhesive aggregates based on different structural element radius. The algorithm first classifies the binary CT image in MATLAB before eroding in order to prevent smaller aggregates from being filtered out in the subsequent segmentation process. Then, structural element with a minimum size is applied to perform ultimate erosion on all the labeled binary images. In the ultimate erosion process, if two or more divided objects appeared in the binary image, then they are classified as possible adhesive aggregates for the next analysis step. If two or more split objects are not observed, then the object could be considered as independent aggregate. Structural elements with different sizes are then selected as the ultimate erosion structural elements for the morphological algorithm of the adhesion particle binary map. According to the adhesive separation map based on morphological algorithm under different structural elements, the separation number with the largest number of division lines is taken as the actual separation of the binary image, and the separation image corresponding to the smallest structural element is taken as the actual separation image. Finally, all the separated marked images are superimposed to form a complete glue separation image. The algorithm flow chart is shown in Figure 1, which was implemented in the MATLAB.

2.2. Numerical Model Reconstructing Method
The three-dimensional microstructural reconstruction of asphalt mixture based on voxel is a reconstruction of complete asphalt mixture specimen from a series of asphalt mixture CT images. This reconstruction is a reverse solution process. For an M × N binary image, the basic constituent unit is a pixel. The expansion of a pixel in three dimensions is a voxel, as shown in Figure 2. A voxel is the most basic unit in the reconstruction of three-dimensional numerical samples.

After CT scanning of the asphalt mixture specimen, a series of continuous CT images can be obtained. These continuous pseudocolor CT images were converted into a matrix of gray images, and the adjacent images were superimposed on each other, thereby expanding into three-dimensional space to form a three-dimensional gray unit space network, as shown in Figure 3.

(a)

(b)
For finite element mesh models, the element, node number, and node coordinates must be specified. For a series of consecutive CT images, the unit and node number of the gray pixel matrix of the first CT image were defined as shown in Figure 4(a), and the unit and node number of the gray pixel matrix of no. k CT image were defined as shown in Figures 4(b) . Figure 5 shows a diagram of an initial three-dimensional numerical model obtained in accordance with the above algorithm.

(a)

(b)

Figure 5 is a numerical model containing the background and mixture. The background, void, asphalt mastic, and aggregate of the CT image were segmented according to the relevant method in Sections 2.1 and 2.2. The segmentation process is shown in Figure 6.

2.3. Reconstruction of Numerical Model of AC-16 Asphalt Mixture
In this paper, an AC-16 asphalt mixture test specimen was used for analysis. According to the above numerical modeling method, four microstructure numerical models with voxel lengths of 0.5 mm, 0.67 mm, 1 mm, and 2 mm, respectively, were reconstructed for AC-16 asphalt mixture. The numerical models of four different voxel lengths for AC-16 asphalt mixture are shown in Figure 7. The eight-node 3D solid reduction integration elements (C3D8R) were adopted in constructing the mesh for the four types of numerical models. Numbers of units and nodes are shown in Table 1.

(a)

(b)

(c)

(d)
3. Material Parameters
Material properties are required for numerical modeling before simulation. The coarse aggregate can be considered to be linear elastic material due to its small deformation. In this paper, modulus and Poisson’s ratio of the coarse aggregate were assumed to be 25 GPa and 0.25, respectively. Asphalt mastic was set as a linear viscoelastic material, and the generalized Maxwell model is often used as its constitutive model [24–26]. The material properties of asphalt mastic in this numerical model adopted the 8-element generalized Maxwell model, and the parameters were obtained by dynamic shear rheological test. Initially, the dynamic shear rheometer (DSR) test was performed at five temperatures (15, 40, 50, 60, and 70°C), and the shear strain level was controlled within 0.15%. Then, the Williams–Landel–Ferry (WLF) equation was applied to shift the measured shear modulus at test temperature to the reference temperature (25°C in this study). The curve of the real part of the complex modulus can be described by the sigmoidal mathematical model:where is the real part of the complex modulus; , , , and are the parameters of the equation; and is the loading time at the reference temperature.
The real part of the complex modulus data was fitted by equation (1), as shown in Figure 8. The generalized Maxwell model parameters can be obtained (Table 2) according to the Prony series representation of the real part of the complex modulus as shown in equation (2) using the collocation method.where is the number of elements of the generalized Maxwell model; is the equilibrium modulus; is the relaxation strength of no. ; and is the shear relaxation time.

Given that the viscoelastic material was set to the generalized Maxwell model, which was the line viscoelastic constitutive, the load magnitude did not affect the analysis results. Therefore, the uniaxial creep simulation test was carried out with a dead load of 80 kPa. The displacement-based boundary conditions were applied to the asphalt mixture digital sample. The test specimen was completely fixed at one end to prevent the rigid body from being displaced, and the load was applied to the other end of the test specimen. For simplicity, the rigid body was used to simulate the loading plates at top and bottom of the specimen. A rigid body is a collection of nodes and elements whose motion is governed by the motion of a reference node. In this study, the dead load was imposed to the reference node, and the rigid body distributed the concentrated load on the top surface of the asphalt mixture digital model. The viscoelastic analysis step was used in this paper. The total analysis step time was 10 and the incremental step time was 0.1.
4. Numerical Simulation and Result Analysis
4.1. Computational Time Analysis
As can be seen in Table 1, the length of voxel became smaller; the number of elements of the model and the number of nodes increased. Taking a model of about 1.7 million elements as an example, the minimum memory required for computer access is approximately 244 GB. So, it is difficult for a PC or a workstation computer to conduct the simulation. Therefore, using supercomputer for massive parallel simulation is necessary. Simulation and analysis of the four numerical models in this study were based on the supercomputer Milkyway-2 in the National Supercomputer Center in Guangzhou. In this paper, massive parallel simulation experiments were conducted on the numerical models with a total of 5 nodes having cumulative capacity of 120 processors and 600 GB of memory.
The computational time of each numerical model is shown in Table 3. The computational time consumed by different voxel length models varied greatly, and the maximum calculation time was 225.92 min. Thus, the computational time can be greatly reduced by increasing the number of applied nodes.
4.2. Creep Stiffness Modulus Analysis
Creep stiffness modulus was used as the result analysis index in this paper. The creep stiffness modulus of the four numerical models is shown in Figure 9. When the voxel lengths were 0.5, 0.67, and 1 mm, the creep stiffness modulus curves were very close.

Since Figure 9 cannot clearly demonstrate differences among the four numerical models’ creep stiffness modulus, stiffness ratios were used to quantify their differences. The ratio of the creep stiffness modulus for each numerical model shown in equation (3) is defined as a fraction of the reference creep stiffness modulus (creep stiffness modulus with a voxel length of 0.5 mm). This is taken as a variation in the evaluation of each creep stiffness modulus curve, as shown in Figure 10.

Figure 10 shows that a longer length of the voxel results in a higher rate of variation. When the length of the voxel was 0.67 mm, the ratio of the creep stiffness modulus from the first second to the tenth second was 0.55%–1.86%. When the length of the voxel was 1 mm, the ratio of the creep stiffness modulus from the first second to the tenth second was 2.14%–4.15%. When the length of the voxel was 2 mm, the ratio of the creep stiffness modulus from the first second to the tenth second was 6.56%–15.58%. Thus, when the accuracy of the CT image of the asphalt mixture reached a certain level, that is, when the length of the voxel is less than or equal to 1 mm, the ratio of the creep stiffness modulus within 10 s is less than 5%, that is, the length of the voxel had no significant effect on the creep test results. When the length of the voxel was 2 mm, the creep stiffness modulus curve was different from the other dimensions, with the maximum change of 15.58%. A smaller voxel size provides more accurate angular information of the aggregate. Voxel size has a great influence on the three-dimensional visualization of the asphalt mixture and the simulation results at local mesolevel. However, remarkably small voxel size showed no significant influence on the creep result. To reduce the computational time of simulation, a numerical model reconstructed with a pixel size of 1 mm can be recommended as an analytical model when performing the creep numerical simulation test.
5. Conclusions and Future Work
5.1. Summary and Conclusions
The size effect of finite elements in the three-dimensional microstructural numerical model of asphalt mixture was studied in this paper. Firstly, the industrial X-ray CT technique was used to obtain the CT image of the internal composition of the asphalt mixture. The Otsu method was used to segment the CT image material phases, and a morphological multiscale segmentation algorithm was proposed to separate the adhesive aggregates after material classification. Finally, a three-dimensional numerical model reconstruction method for asphalt mixture based on voxel was proposed from a series of discrete asphalt mixture CT images. Four microstructural numerical models with voxel lengths of 0.5, 0.67, 1, and 2 mm were reconstructed. Uniaxial creep virtual test was implemented on the numerical model with four voxel sizes using the supercomputer Milkyway-2 in the National Supercomputer Center in Guangzhou to analyze the influence of the size effect on the simulation results, and the following conclusions were obtained:(1)A morphological multiscale algorithm for ultimate erosion of coarse aggregate adhesion image of asphalt mixture was proposed based on structural elements with different sizes. The algorithm can effectively separate adhesive aggregate particles.(2)A voxel-based three-dimensional numerical model of AC-16 asphalt mixture reconstruction method was proposed. Numerical model was successfully reconstructed based on CT images of asphalt mixture processed by the Otsu method and morphological multiscale algorithm.(3)The numerical models with four different voxel sizes varied greatly in terms of the computational time when performing the creep numerical simulation test, and the maximum computational time was 225.92 min. Thus, more cumulative capacity of processors and memory can remarkably reduce the computational time.(4)When the length of the voxel was less than 1 mm, the creep stiffness modulus changed within 2.5%. This finding suggests that when the voxel size is less than 1 mm, the influence of the voxel size on creep stiffness modulus is not significant. Moreover, voxel size mainly affects the computational time and the angular information of the aggregate. A numerical model reconstructed with 1 mm voxel size is recommended as the analytical model.
5.2. Recommendations for Future Work
(1)Many researchers have developed the low-resolution digital samples to simulate the behavior of asphalt mixture subjected to specific loading conditions, and this approach has been verified as a relatively reliable method to conduct the virtual test compared with the experimental test results. It is vital to validate the high-resolution simulation results with a series of experimental tests, and it will be the focus in our future papers.(2)The aggregate-mastic interface is a very challenging task in the asphalt mixture virtual simulation. For simplicity, the mastic phase was assumed to be perfectly bonded to the aggregate phase ignoring the internal cohesive behavior of the two phases in this paper. In order to obtain a more accurate digital sample, the interfacial conditions will be considered as part of our ongoing modeling efforts.(3)Asphalt mixture is a highly heterogeneous material consisting of aggregate, mastic, and air voids. The asphalt mastic has very complex temperature-, time-, and rate-dependent three-dimensional viscoelastic property. In this paper, the mastic is considered as a one-dimensional linear viscoelastic material. The underway paper will be focused on a three-dimensional viscoelastic constitutive model developed for asphalt mastic.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
This study was supported by the Science and Technology Program of Guangzhou (no. 201804010231) and National Natural Science Foundation of China (no. 51878193). The authors thank all those who contributed in the X-ray CT image technique and the experimental part of this study. In addition, the authors acknowledge the National Supercomputer Center in Guangzhou (NSCC) for providing computing resources useful in conducting the research reported in this paper.