Using a Molecular Dynamics Simulation to Investigate Asphalt Nano-Cracking under External Loading Conditions

Recent research shows that macro-scale cracking in asphalt binder may originate from its intrinsic defects at the nano-scale. In this paper, a molecular dynamics (MD) simulation was conducted to evaluate the nucleation of natural defects in asphalt. The asphalt microstructure was modeled using an ensemble of three different types of molecules to represent a constituent species: asphaltenes, naphthene aromatics and saturates, where the weight proportion of 20:60:20 was used to create an asphalt-like ensemble of molecules. Tension force was then applied on the molecular boundaries to study the crack initiation and propagation. It was discovered that the natural distribution of atoms at microscale would affect the intrinsic defects in asphalt and further influence crack initiation and propagation in asphalt.


Background and Introduction
Asphalt pavement is one of the main highway forms throughout the world. With the continuous development of the world economy, there have been tremendous numbers of vehicles passing on the pavements every day, causing significant irreversible damage to pavement. It is therefore necessary to study asphalt material performance, since it is directly relevant to its service life, where structure is one of the most important factors that affects its mechanical performance. Recent studies show that the asphalt material microstructure may affect its properties at the macro-scale [1,2].
There have been many experimental approaches on the asphalt material microstructure due to the development of modern testing instruments. Mikhailenko et al. [3], used scanning electron microscopy (SEM) to observe the microstructure of asphalt materials. Fourier transform infrared spectroscopy (FTIR) can be used to study the distribution and contents of asphalt components [4][5][6]. Gel permeation chromatography (GPC) can determine the molecular weight and distribution of asphalt rapidly and reliably [7]. Differential scanning calorimetry (DSC) can be used to analyze the difference of the thermal stability of asphalt [8,9]. Atomic force microscope (AFM) can be used to observe three-dimensional (3D) spatial images of asphalt [10][11][12]. These new experimental methods can be used to further investigate the effect of modifiers on the modification of asphalt materials at nanometer scales [13,14]. For example, SEM can be used to observe and compare the microstructures of different nanometer modified asphalts [15], and AFM, FTIR, GPC and DSC can be used to further analyze the modify mechanism of nano-modifier materials.
There also have been many studies on asphalt material properties using the simulation, where molecular dynamics (MD) simulation is one of the effective methods to investigate the microstructure and micro-properties of asphalt materials [16][17][18]. In the aspect of asphalt oxidation aging, MD simulation can be employed to investigate internal chemical and molecular mechanical property changes [19]. Molecular Dynamics simulation can also be used to determine the thermodynamic properties of asphaltene before and after oxidative aging [20]. In other aspects, molecular dynamics simulation can also be used to study molecular movements and microstructure changes of epoxy modified asphalt in the process of healing [21]. Xu et al. [22], conducted AFM experiments and MD simulation to study the adhesion between asphalt and aggregate. The mechanical properties of asphalt and aggregate can be further investigated by MD simulation at the nano-scale [23], and by investigating mechanical properties of the asphalt-aggregate interface [24].
Although there has been significant progress in asphalt material microstructure and micro-properties, research on the fundamental mechanism of asphalt microstructure cracking using simulations and experiment approaches are very few. In this study, MD simulation was conducted to study micro-crack initiation, and propagation of typical asphalt molecular structures. The ultimate goal was to try to understand the fundamental mechanism of asphalt molecular structure cracking.

Molecular Dynamics Simulations
In this study, the MD simulation was employed to investigate the mechanical behavior of asphalt molecules.

Asphalt Molecular Structure Construction
Asphalt is often characterized chemically by its hydrocarbon class composition, e.g., three main constituent species including asphaltene, naphthene aromatics and saturates. The three-fraction asphalt model was employed to establish the molecular structure following Zhang and Greenfield [25], where three different types of molecules are used to represent the corresponding constituent species, and the three are then formed together to create an asphalt-like ensemble of molecules. In this model, 1,7-dimethylnapthalene was used to represent naphthene aromatics, and n-C 22 (n-doccosane) molecules were used to represent saturates, with the mass fraction of asphaltenes, naphthene aromatics, and saturates being approximately 20:60:20, respectively.
The condensed-phase optimized molecular potentials for atomistic simulation studies (COMPASS) force field was used for our simulation, considering that this force field could accurately describe the inorganic material properties. The initial molecular structure was constructed by using the commercial molecular computation software Material Studio 7.0 (BIOVIA Company, San Diego, CA, USA). The construction function in amorphous cell tools of Material Studio was used to construct the representative volume element in asphalt binder. The target density of the final configurations was set to 1 g/cm 3 . Table 1 shows the detailed configuration of the asphalt-like ensemble. Figure 1 shows the constructed molecular structure, where the edge length (Å) of the lattice is 50.9 × 50.9 × 50.9. The atoms in each molecule in the asphalt-like amorphous cell were then subjected to conjugate gradient energy minimization within 10,000 fs.  Our molecular dynamics simulation was split into two parts: constructing the model in Materials Studio due to convenience, and applying the loading conditions in the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) software. After construction of the asphalt molecular model in Material Studio, it was exported in a PDB-format file and imported into the VMD (Visual Molecular Dynamics) software (University of Illinois at Urbana-Champaign, Urbana-Champaign, IL, USA). By using the "Topo writelammpsdata asphalt.data full" command, the asphalt data file can be read by the LAMMPS software.
The following potential function was used in our LAMMPS simulation for convenience: where E is the potential, is the depth of a potential well that reflects the attraction between two atoms, is the distance between atoms when the action potential is equal to 0, r is the distance between two atoms, and rc is the cutoff.

Simulation Results
The relaxation process of asphalt molecules is shown in Figure 2. The relaxation was used to reach the energy minimization state of the asphalt molecular structure. It can be seen that the relaxation process would influence the structure of asphalt molecules, where agglomeration and densification occurred. Our molecular dynamics simulation was split into two parts: constructing the model in Materials Studio due to convenience, and applying the loading conditions in the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) software. After construction of the asphalt molecular model in Material Studio, it was exported in a PDB-format file and imported into the VMD (Visual Molecular Dynamics) software (University of Illinois at Urbana-Champaign, Urbana-Champaign, IL, USA). By using the "Topo writelammpsdata asphalt.data full" command, the asphalt data file can be read by the LAMMPS software.
The following potential function was used in our LAMMPS simulation for convenience: where E is the potential, ε is the depth of a potential well that reflects the attraction between two atoms, σ is the distance between atoms when the action potential is equal to 0, r is the distance between two atoms, and r c is the cutoff.

Simulation Results
The relaxation process of asphalt molecules is shown in Figure 2. The relaxation was used to reach the energy minimization state of the asphalt molecular structure. It can be seen that the relaxation process would influence the structure of asphalt molecules, where agglomeration and densification occurred. In the study, molecular dynamics simulation were exerted on the asphalt layer, where "region" and "group" commands were applied to set different regions in the asphalt molecule ensemble. The In the study, molecular dynamics simulation were exerted on the asphalt layer, where "region" and "group" commands were applied to set different regions in the asphalt molecule ensemble. The velocity command was used to set the velocities of different groups at 0.3 Å/fs along the y axis to simulate the state of tension. The ensemble temperature was set to room temperature, which is 298.15 K. For every 100 time steps, the simulation temperature was adjusted. For asphalt molecular ensemble, the energy minimization process was firstly conducted to achieve the stability of the system, and then the simulation of the stretching process was conducted. The time step was set to 0.003 fs, and a total operation time with 6000 steps was conducted. In the OVITO (Open Visualization Tool) software, asphalt molecules were sliced (Figure 3), where the internal development of asphalt molecules under tension could be observed. Figure 3a shows the complete whole microstructure and Figure 3b shows the half-sliced part. It was found that, along with the applying of continuous tension loading, obvious gaps could be found in the internal asphalt molecule structure (Figure 4), where these locations may be the weak regions and prone to easily crack.
Appl. Sci. 2017, 7, 770 5 of 10 velocity command was used to set the velocities of different groups at 0.3 Å/fs along the y axis to simulate the state of tension. The ensemble temperature was set to room temperature, which is 298.15 K. For every 100 time steps, the simulation temperature was adjusted. For asphalt molecular ensemble, the energy minimization process was firstly conducted to achieve the stability of the system, and then the simulation of the stretching process was conducted. The time step was set to 0.003 fs, and a total operation time with 6000 steps was conducted. In the OVITO (Open Visualization Tool) software, asphalt molecules were sliced (Figure 3), where the internal development of asphalt molecules under tension could be observed. Figure 3a shows the complete whole microstructure and Figure 3b shows the half-sliced part. It was found that, along with the applying of continuous tension loading, obvious gaps could be found in the internal asphalt molecule structure (Figure 4), where these locations may be the weak regions and prone to easily crack.  The radial distribution function (RDF) is an important parameter to describe molecular structure. The essence of the radial distribution function gab (r) is the atomic ratio of local density and average density in the system [26]. When A and B atoms are distributed uniformly in the spherical shell, the probability is calculated as follows:  velocity command was used to set the velocities of different groups at 0.3 Å/fs along the y axis to simulate the state of tension. The ensemble temperature was set to room temperature, which is 298.15 K. For every 100 time steps, the simulation temperature was adjusted. For asphalt molecular ensemble, the energy minimization process was firstly conducted to achieve the stability of the system, and then the simulation of the stretching process was conducted. The time step was set to 0.003 fs, and a total operation time with 6000 steps was conducted. In the OVITO (Open Visualization Tool) software, asphalt molecules were sliced (Figure 3), where the internal development of asphalt molecules under tension could be observed. Figure 3a shows the complete whole microstructure and Figure 3b shows the half-sliced part. It was found that, along with the applying of continuous tension loading, obvious gaps could be found in the internal asphalt molecule structure (Figure 4), where these locations may be the weak regions and prone to easily crack.  The radial distribution function (RDF) is an important parameter to describe molecular structure. The essence of the radial distribution function gab (r) is the atomic ratio of local density and average density in the system [26]. When A and B atoms are distributed uniformly in the spherical shell, the probability is calculated as follows: The radial distribution function (RDF) is an important parameter to describe molecular structure. The essence of the radial distribution function g ab (r) is the atomic ratio of local density and average density in the system [26]. When A and B atoms are distributed uniformly in the spherical shell, the probability is calculated as follows: where V is the total volume of the simulation system, n A and n B are the number of A and B atoms, ρ A and ρ B are the average density of atoms in the simulation system, and P A and P B are the uniform distribution probabilities of A and B atoms. By calculating the radial distribution function between different atoms, the local distribution of different atoms could be studied, as shown in Figure 5. It is observed that there existed a small difference in RDF before and after tension loading. For the asphalt microstructure before tension loading, when the distance was 2.375 Å, a maximum value of 16.5250 occurred; for the asphalt microstructure after the tension loading, when the distance was 2.375 Å, maximum value of 17.039 occurred. The trends of the curves of the two curves were completely consistent, indicating that most of the atoms in the asphalt microstructure did not significantly change their locations while a small portion changed the locations. It was further concluded that the applied tension loading would not only result in the whole asphalt microstructure change, but also in some of the specific locations. where V is the total volume of the simulation system, and are the number of A and B atoms, and are the average density of atoms in the simulation system, and PA and PB are the uniform distribution probabilities of A and B atoms.
By calculating the radial distribution function between different atoms, the local distribution of different atoms could be studied, as shown in Figure 5. It is observed that there existed a small difference in RDF before and after tension loading. For the asphalt microstructure before tension loading, when the distance was 2.375 Å, a maximum value of 16.5250 occurred; for the asphalt microstructure after the tension loading, when the distance was 2.375 Å, maximum value of 17.039 occurred. The trends of the curves of the two curves were completely consistent, indicating that most of the atoms in the asphalt microstructure did not significantly change their locations while a small portion changed the locations. It was further concluded that the applied tension loading would not only result in the whole asphalt microstructure change, but also in some of the specific locations. The application of tension loadings on the asphalt molecular structure was done by setting a constant speed in its y direction (vertical direction) of different regions, resulting in the stretching of the asphalt microstructure. Figure 6 shows the morphological changes of asphalt microstructures. Note that different color represents different atom types on the top and three different divisions on the bottoms. It was noticed that by applying the tension force, the original "perfect" structure generally stretched, and a "hole", i.e., natural defect, occurred in the microstructure, which could be seen as the initiation of the microcrack. It was expected that with an increase in the tension loading, the inside small hole propagated and finally resulted in the failure of the asphalt microstructure [27].  The application of tension loadings on the asphalt molecular structure was done by setting a constant speed in its y direction (vertical direction) of different regions, resulting in the stretching of the asphalt microstructure. Figure 6 shows the morphological changes of asphalt microstructures. Note that different color represents different atom types on the top and three different divisions on the bottoms. It was noticed that by applying the tension force, the original "perfect" structure generally stretched, and a "hole", i.e., natural defect, occurred in the microstructure, which could be seen as the initiation of the microcrack. It was expected that with an increase in the tension loading, the inside small hole propagated and finally resulted in the failure of the asphalt microstructure [27].
In order to analyze the temperature effects on the asphalt molecular, four kinds of different temperatures were chosen, including −40, 0, 25 and 60 • C, to analyze the temperature influence on the morphology and properties of the asphalt molecules during tension. The time step was 0.003 fs, with a total of 5000 steps. The asphalt molecular morphology with different temperture at 5000 timesteps are shown in Figure 7. It was seen that temperature seemed do not have a significant influence on the asphalt microstructure evolution. It might be caused by the small molecule amount in the model. Compared with the internal molecular component distribution, temperature does not have a significant effect on the molecular structure for this computation scale. To further study the temperature effects on asphalt microstructure evolution quantitatively and in more detail, future studies will be conducted on a supercomputer.
The application of tension loadings on the asphalt molecular structure was done by setting a constant speed in its y direction (vertical direction) of different regions, resulting in the stretching of the asphalt microstructure. Figure 6 shows the morphological changes of asphalt microstructures. Note that different color represents different atom types on the top and three different divisions on the bottoms. It was noticed that by applying the tension force, the original "perfect" structure generally stretched, and a "hole", i.e., natural defect, occurred in the microstructure, which could be seen as the initiation of the microcrack. It was expected that with an increase in the tension loading, the inside small hole propagated and finally resulted in the failure of the asphalt microstructure [27].  In order to analyze the temperature effects on the asphalt molecular, four kinds of different temperatures were chosen, including −40, 0, 25 and 60 °C, to analyze the temperature influence on the morphology and properties of the asphalt molecules during tension. The time step was 0.003 fs, with a total of 5000 steps. The asphalt molecular morphology with different temperture at 5000 timesteps are shown in Figure 7. It was seen that temperature seemed do not have a significant influence on the asphalt microstructure evolution. It might be caused by the small molecule amount in the model. Compared with the internal molecular component distribution, temperature does not have a significant effect on the molecular structure for this computation scale. To further study the temperature effects on asphalt microstructure evolution quantitatively and in more detail, future studies will be conducted on a supercomputer.
Through the analysis, it was found that there were "holes" in the asphalt molecules by stretching, with these holes existing due to the distribution of molecules, creating the original defects. The results showed that the original defects naturally existed in the asphalt molecules, resulting in different cracking performance. Note that since this was a preliminary study on the nucleation of natural defects in asphalt, we did not study this phenomenon quantitatively. In future studies, a parameter such as "total distance moved" by the atoms from some reference points may be proposed to describe this phenomenon. In order to analyze the temperature effects on the asphalt molecular, four kinds of different temperatures were chosen, including −40, 0, 25 and 60 °C, to analyze the temperature influence on the morphology and properties of the asphalt molecules during tension. The time step was 0.003 fs, with a total of 5000 steps. The asphalt molecular morphology with different temperture at 5000 timesteps are shown in Figure 7. It was seen that temperature seemed do not have a significant influence on the asphalt microstructure evolution. It might be caused by the small molecule amount in the model. Compared with the internal molecular component distribution, temperature does not have a significant effect on the molecular structure for this computation scale. To further study the temperature effects on asphalt microstructure evolution quantitatively and in more detail, future studies will be conducted on a supercomputer.
Through the analysis, it was found that there were "holes" in the asphalt molecules by stretching, with these holes existing due to the distribution of molecules, creating the original defects. The results showed that the original defects naturally existed in the asphalt molecules, resulting in different cracking performance. Note that since this was a preliminary study on the nucleation of natural defects in asphalt, we did not study this phenomenon quantitatively. In future studies, a parameter such as "total distance moved" by the atoms from some reference points may be proposed to describe this phenomenon. The strain and stress developments over time for the asphalt microstructure subject to external loading at 25 °C are shown in Figure 8. Figure 8a shows that the strain increased almost linearly, as we applied a constant strain deformation rate. Figure 8b shows that the stress gradually decreased and remained almost constant after 1300 timesteps, indicating that due to the occurrence of the "internal hole", the asphalt microstructure bore less stress, and when the final structure was stabilized, the stress became stable. Through the analysis, it was found that there were "holes" in the asphalt molecules by stretching, with these holes existing due to the distribution of molecules, creating the original defects. The results showed that the original defects naturally existed in the asphalt molecules, resulting in different cracking performance. Note that since this was a preliminary study on the nucleation of natural defects in asphalt, we did not study this phenomenon quantitatively. In future studies, a parameter such as "total distance moved" by the atoms from some reference points may be proposed to describe this phenomenon.
The strain and stress developments over time for the asphalt microstructure subject to external loading at 25 • C are shown in Figure 8. Figure 8a shows that the strain increased almost linearly, as we applied a constant strain deformation rate. Figure 8b shows that the stress gradually decreased and remained almost constant after 1300 timesteps, indicating that due to the occurrence of the "internal hole", the asphalt microstructure bore less stress, and when the final structure was stabilized, the stress became stable. Appl. Sci. 2017, 7, 770 8 of 10 It should be mentioned that our current research results seem to be a bit different from the observations in previous studies: 1. Menapace et al. (2016) discovered that temperature will significantly affect the "Bumble-bee" structure in asphalt samples at microscale [28]; 2. Jahangir et al. (2015) claimed that the tension loading would result in the change in asphalt microstructure [29]. There might be two reasons for the differences: 1. We currently only apply a very small tension loading. Therefore, compared with the internal molecule component distribution, the latter will have a larger effect on the molecular structure; 2. The current molecular model is established with relatively small molecule amount due to the limitation of our computer ability. It is expected that a much larger molecular model with millions of molecules may resemble the results by the previous research.

Conclusions
In this study, the MD simulation was conducted to model the asphalt molecular structure. The following conclusions were made: 1. It was discovered that by applying the tension force, the original "perfect" structure generally stretched and a "hole", i.e., natural defect, occurred in the microstructure, which can be seen as the initiation of the microcrack. It was expected that with an increase of tension loading, the inside small hole would propagate and finally result in the failure of the asphalt microstructure. 2. It was observed that there existed a difference in RDF before and after tension loading. For the asphalt microstructure before the tension loading, when the distance was 2.375 Å, a maximum value of 16.5250 occurred; for the asphalt microstructure after the tension loading, when the distance was 2.375 Å, maximum value of 17.039 occurred. The trends of the curves for the two curves were completely consistent, indicating that many of the atoms in the asphalt microstructure did not change their locations, while a small portions changed.