Molecular Simulation of Minerals-Asphalt Interfacial Interaction

The interfacial interaction between asphalt binder and mineral aggregate makes different components have different diffusion behavior. It influences the performance of interface and consequently that of the mix. In this research, the models of four asphalt components (asphaltene, resin, aromatics and saturate) and five minerals were constructed individually, and then the Al2O3-asphalt interface model was constructed by adding the asphalt layer and mineral layer together. The interfacial behavior at molecular scale was simulated by setting boundary conditions, optimizing the structure and canonical ensemble. The mean square displacement (MSD) and diffusion coefficient of particles were selected as indicators to study the diffusion of asphalt components on the surface of Al2O3. The results show that when the temperature was 65 ◦C, asphalt binder showed more viscosity, the diffusion speed of asphalt components ranked according to their molecular mass, which was saturate > aromatics > resin > asphaltene. At 25 ◦C and 165 ◦C, the resin had the fastest diffusion speed, which was nearly twice of asphaltene. The interaction between asphalt components and Al2O3 mineral surface can make the diffusion of asphalt components independent on temperature.


Introduction
Molecular dynamics simulation technology is a computer modeling approach.This approach employs classical mechanics model of molecular interactions as a foundation.With the development of computer technology, molecular dynamics (MD) simulation is widely being used in the fields of material science, biology, medicine and physical chemistry.In this research, molecular dynamics simulation technology was used to study the interfacial adsorption between asphalt binder and mineral aggregate.This is significantly important to understand the interfacial mechanism between asphalt binder and mineral aggregate at nanoscale [1][2][3][4][5].
The nano-mechanics method has recently emerged as a powerful computational approach to model from nano to macro microstructure deformation and failure behaviours in materials, including asphaltic materials.Bhasin et al. [6] used the molecular simulation techniques to investigate the correlation of chain length and chain branching to self-diffusivity of binder molecules.The findings reported by them were consistent with observations reported in previous studies and expanded on the understanding of the relationship between molecular architecture, self-diffusivity, and self-healing properties of asphalt binders.Sun et al. [7] used molecular dynamics simulation to examine the hypothesis of healing mechanism and to evaluate the self-healing capability of neat and Styrene Butadiene Styrene (SBS) modified asphalt binders.Densities and glass transition temperatures of neat and SBS modified asphalt binders were obtained by them based on the MD simulation results.They created an artificial crack in MD model to simulate the diffusion of molecules across a crack interface.They proposed three indices calculated to evaluate the self-healing capability of asphalt binders, including diffusion coefficient, activation energy, and pre-exponential factor.They also compared the values of these three indices obtained from MD results with the values derived from fatigue-rest-fatigue test, to assess the reliability of above indices as self-healing performance indicators.Xu and Wang [8] used molecular dynamics simulation to investigate the effect of oxidative aging on asphalt binder.They used mean square displacement (MSD) and radial distribution function (RDF) of molecules and diffusion coefficient to analyze molecular structures of virgin and aged asphalt.They found that the molecular structure of asphalt was in agreement with the colloidal theory.Aging weakened the nano-aggregation behavior of asphaltene molecules and reduces the translational mobility of asphalt molecules.Compared to virgin asphalt, aged asphalt had the higher activation energy barrier for self-healing.Inclusion of water molecules in the asphalt-aggregate interface caused degradation of work of adhesion and this effect was more significant for aged asphalt.
Some researchers used MD method to interfacially investigate the interaction between asphalt binder and mineral aggregate.Lu and Wang [9] used the molecular mechanics method to calculate the quartz bulk elastic constants, e.g.stiffness matrix, shear modulus, Young's modulus and Poisson's ratio.They also used the MD simulations to model the deformation and failure behavior of asphalt-rock interfaces when subject to uniaxial tension.A 3D asphalt-quartz interface structure model was proposed in terms of density, position and thickness.The interfacial atom trajectories were visualized to represent the properties in simulations that characterized those of the asphalt-rock interface under uniaxial tension at nanoscale.Interfacial debonding characteristics or the adhesive failure was implemented with a large-scale MD simulation technology.They found that highly anisotropic elastic properties of a quartz structure will appear from atomistic scale and tensile strength of the asphalt-quartz interface system was controlled by the stress state at the asphalt-rock interface layer.Asphalt-rock interface adhesive failure appeared to be ductile at freezing environmental temperature and low strain rate.Later, Lu and Wang [10] reviewed the advances of nano-mechanics modelling method.They introduced the recent advances of their research on the molecular origin of deformation and failure processes of asphalt-aggregate interfaces, including studies of bulk mechanical properties of asphalt and mineral aggregate, as well as asphalt-aggregate interface properties under tensile or confined shear loading.Their nanoscale modelling work focused on understanding the mechanical quantities of hierarchical asphalt/mineral materials to develop atomic scale structure-property relations.The large-scale atomistic simulation provided a powerful approach in studying mechanical properties of asphalt concrete materials.They also concluded a discussion of the significance of nanoscale structure-properties relationship of the interface structure under various loading conditions and water affinities of specific asphalt-aggregate pairs.Li et al. [11] used MD to establish the molecular models of asphaltene and four oxides, and calculated the interfacial energy.They found that the interfacial energy between asphaltene and oxide reaches the maximum value at 25 • C and 80 • C and the minimum value at 40 • C. Dong et al. [12] used the MD simulation and the peak force quantitative nano-mechanics atomic force microscopy (PFQN-AFM) to characterize the nanostructure of the asphalt-aggregate interface.They found that that the intermolecular gaps can be filled by naphthene aromatics, and molecular micelles can be lubricated by saturates.Xu and Wang [13,14] used MD to study thermodynamic and cohesive properties of asphalt binder, such as density, solubility parameter, cohesive energy density, and surface free energy.They analyzed the adhesion properties by calculating the interaction energy and the work of adhesion at asphalt-aggregate interface for the first time.They also performed tensile simulations to obtain the stress-separation responses to analyze the interfacial mechanical behavior.They found that the interface failure was mainly adhesive failure although large air voids were formed in the bulk asphalt as the loading rate decreases to a certain level.The objective of this research is to build an effective interfacial model between asphalt binder and minerals and study the factors that influence interfacial interaction.

The Construction of Asphalt Model
The molecular structure of four asphalt components selected in this research can be seen in Reference [15].The molecular weight of different asphalt component can be calculated, and their mass ratio was selected according to the Reference [16].The relative atomic mass of C, H, O, N, S are 12, 1, 16, 14, 32 respectively.The weight of element equals atomic weight multiplying atomic number.The molecular weight can be obtained by accumulating all elements weights.The calculated results are 2048 for asphalt molecular, 824 for resin, 618 for aromatic and 310 for saturate.Accordingly, the ratio of the four asphalt components can be calculated as follows: asphaltenes:resin:aromatics:saturate = 1:3:7:5.

The Construction of Mineral Model
Other researchers have already found that SiO 2 , Al 2 O 3 , CaO, MgO and Fe 2 O 3 are the most important components affecting the performance of aggregate, particularly for the tensile strength and adhesion behavior [17,18].Therefore, this research mainly focused on the crystal structure of Al 2 O 3 , and analyzed the interfacial diffusion behavior between it and asphalt binder.
(1) Construction of crystal unit cell Crystal unit cell is the smallest unit in crystal.Its size and shape can be characterized by the edge length and crossing angle.The lattice constants used in this research can be found in Reference [15].
(2) Interception and optimization of lattice planes In order to describe the structure of crystal unit cell clearly, a plane was first intercepted.For example, the (001) lattice plane of Al 2 O 3 was intercepted and its thickness was set as 10 Å.The energy minimization procedure was used to optimize the structure, and then the supercell was constructed, followed by adding a vacuum layer with 10 Å thickness.The finally constructed Al 2 O 3 mineral surface model is shown in Figure 1.
The objective of this research is to build an effective interfacial model between asphalt binder and minerals and study the factors that influence interfacial interaction.

The Construction of Asphalt Model
The molecular structure of four asphalt components selected in this research can be seen in Reference [15].The molecular weight of different asphalt component can be calculated, and their mass ratio was selected according to the Reference [16].The relative atomic mass of C, H, O, N, S are 12, 1, 16, 14, 32 respectively.The weight of element equals atomic weight multiplying atomic number.The molecular weight can be obtained by accumulating all elements weights.The calculated results are 2048 for asphalt molecular, 824 for resin, 618 for aromatic and 310 for saturate.Accordingly, the ratio of the four asphalt components can be calculated as follows: asphaltenes:resin:aromatics:saturate = 1:3:7:5.

The Construction of Mineral Model
Other researchers have already found that SiO2, Al2O3, CaO, MgO and Fe2O3 are the most important components affecting the performance of aggregate, particularly for the tensile strength and adhesion behavior [17,18].Therefore, this research mainly focused on the crystal structure of Al2O3, and analyzed the interfacial diffusion behavior between it and asphalt binder.
(1) Construction of crystal unit cell Crystal unit cell is the smallest unit in crystal.Its size and shape can be characterized by the edge length and crossing angle.The lattice constants used in this research can be found in Reference [15].
(2) Interception and optimization of lattice planes In order to describe the structure of crystal unit cell clearly, a plane was first intercepted.For example, the (001) lattice plane of Al2O3 was intercepted and its thickness was set as 10 Å .The energy minimization procedure was used to optimize the structure, and then the supercell was constructed, followed by adding a vacuum layer with 10 Å thickness.The finally constructed Al2O3 mineral surface model is shown in Figure 1.

The Construction of Interface Model
The interface model was obtained by overlaying the following three layers.The first layer was mineral crystal surface.The second layer was asphalt layer.The third layer was a 30 Å vacuum layer.The constructed models have a size 35 Å × 35 Å × 75 Å .The details and verification on effectiveness of the models can be found in the following section.

Molecular Dynamic Simulation
Condensed-phase optimized molecular potentials for atomistic simulation studies (COMPASS) force field was selected in this research.Its potential energy function is shown in Equation ( 1).
where, i and j are different atoms; q is the atomic charge (C); r is the distant between different atoms (m); ε is the potential well depth (eV).More detailes about structure optimization, boundary conditions and molecular dynamic simulation can be found in Reference [19].Particularly this research selected a single cell for the molecular dynamics simulation and a periodic boundary condition was used.A 12.5 Å Cut-off radius was used to calculate the non-bonded remote force.

Verification of Model Effectiveness
The verification of effectiveness of asphalt model could be found in Guo et al.'s previous work [15].In this paper, the effectiveness of mineral crystal models is discussed.
In this research, the elastic constants of constructed structure models were calculated as follows: (i) the surface model of mineral crystals was input by Materials Studio software (version 8.0, BIOVIA, San Diego, CA, USA); (ii) the Calculation function in the Forcite module was selected; (iii) the mechanical properties was selected in the task window.The calculated accuracy was Medium, and the calculated method was Constant Strain.The calculation principle was as follows: (i) First, select a value for strain tensor according to the principle that the relationship between strain amplitude, elastic constant and stress tensor should be the simplest.(ii) According to the relationship between strain tensor and lattice matrix vector, the lattice matrix vector after strain was obtained.(iii) The fractional coordinates of atoms between different strains were selected as the atom position, and the stress tensor was calculated after optimizing the atom position.(iv) The elastic constant was obtained by analyzing the relationship between stress tensor, strain amplitude and elastic constant.
The elastic constants of different mineral crystal models in this research were calculated as follows.
The elastic stiffness constant SiO 2 crystal structure model: The physical properties of different minerals calculated in this research are shown in Table 1.It can be seen from Table 1 that the calculated shear modulus of SiO 2 was 46.7GPa, which was similar with the quartz's shear modulus 44.5 GPa reported in Reference [20].The Young's modulus of SiO 2 were all 130.2 GPa along three directions, which were larger than the tested value 76-97 GPa reported in Reference [20].It was due to that the purity of the practical material was poorer than the simulated material and the practical material had more defects.The calculated shear modulus of CaO was 76.8 GPa, larger than the tested value 35 GPa.The Young's modulus of CaO were all 193.9 GPa along three directions, which were larger than the tested value 72.4-88.2GPa reported in Reference [20].Huang et al. found that the Young's modulus of kaolin mainly comprised of SiO 2 and Al 2 O 3 was between 31 GPa and 70 GPa [21].This was smaller than the calculated results in this research.The difference between the simulated results and tested results was due to: (i) the simulated mineral was not exactly the same with the practical mineral, and the latter had more impurities; (ii) the practical mineral structure had more crystal defects compared to the simulated minerals.Overall, the simulated values in this research are acceptable compared to the practical tested values.To some extent, they can reflect the practical situation.

Diffusion of Asphalt Components on the Surface of Al 2 O 3
Mean square displacement (MSD) was used to study the diffusion trend of four asphalt components on the surface of Al 2 O 3 .MSD can be calculated as follows: where, r i (t) is the displacement of particle i at t time (m); r i (0) is the displacement of particle i at the initial time (m); <> is the average of all particles.
The bigger the slope of MSD curve, the faster the migration of the particle.The diffusion coefficient of particles can be calculated according to Equations (5).
where, N is the number of diffusion particles in the system; D is the diffusion coefficient of particle (m 2 /s).If the differential in Equation ( 5) can be approximated as the slope of correlation between MSD and time, the diffusion coefficient of particles can be revised as follows: where, a is the slope of the curve between MSD and time.
The relationship between MSD of four asphalt components on the surface of Al 2 O 3 and time are shown in Figures 2 and 3. where, a is the slope of the curve between MSD and time.
The relationship between MSD of four asphalt components on the surface of Al2O3 and time are shown in Figures 2 and 3  It can be seen from Figure 2 that the diffusion speed of saturate is the fastest on the surface of Al2O3.It was not dependent on the temperature.At 25 °C , resin had a similar diffusion speed with saturate, obviously faster than asphaltene and aromatic.At 65 °C , asphaltene, resion and aromatic had similar diffusion speed.At 165 °C , the resin diffused slower than saturate, followed by aromatic.The asphaltene diffused the slowest.The diffusion speeds of different asphalt components are consistent with the ranking of the relative molecular mass.It indicates that the interaction between asphalt components and SiO2 at 165 °C became too weak to influence the diffusion of asphalt components.Figure 3 indicates that the diffusion speeds of asphaltene, resin and aromatic on surface of Al2O3 were not sensitive to temperature.For saturate, the diffusion speed at 165 °C was obviously higher than the diffusion speeds at 25 °C and 65 °C .This is due to that the higher temperature can result in the faster diffusion.The slope of curves in Figures 2 and 3 can be obtained by fitting the relationship between MSD and time.The diffusion coefficients of four asphalt components at different temperature can be calculated according to Equation (6).The results are shown in Figure 4.    Figure 3 indicates that the diffusion speeds of asphaltene, resin and aromatic on surface of Al2O3 were not sensitive to temperature.For saturate, the diffusion speed at 165 °C was obviously higher than the diffusion speeds at 25 °C and 65 °C .This is due to that the higher temperature can result in the faster diffusion.The slope of curves in Figures 2 and 3 can be obtained by fitting the relationship between MSD and time.The diffusion coefficients of four asphalt components at different temperature can be calculated according to Equation (6).The results are shown in Figure 4.It can be seen from Figure 4 that at 65 • C, around the softening point of asphalt with penetration 60-100, the diffusion speed of asphalt components ranked according to their molecular weight, that is the larger the molecular weight, the faster the diffusion would be.At 25 • C and 165 • C, the rule also applied except resin.This may due to that the resin was adsorbed around the asphaltene.The interaction between different components made the resin diffused faster.By comparing the diffusion speeds at different temperatures, it can be known that the diffusion of four asphalt components were not sensitive to the temperature.This was due to the interaction between asphalt components and Al 2 O 3 mineral crystal.
seen from the above matrix of elastic stiffness constant of different minerals that the elastic stiffness constants of SiO 2 , CaO and MgO were same at three principal directions.It indicates that the structure of SiO 2 , CaO and MgO were symmetry, and the hardness were same at the a, b and c edge-length directions.However, the C 33 of Al 2 O 3 was 488.3 GPa, a little bit larger than C 11 = C 22 = 479.1 GPa.This indicates that it was more difficult for Al 2 O 3 crystal to deform along the c direction than the a and b direction.The C 33 of Fe 2 O 3 was 185.1GPa, smaller than C 11 = C 22 = 533.8GPa.It means that it was much easier for Fe 2 O 3 crystal to deform along the c direction than the a and b direction. .

Figure 2 .
Figure 2. The relationship between mean square displacement (MSD) of four asphalt conponents on the surface of Al2O3 and time: (a) Temperature is 25 °C; (b) Temperature is 65 °C; (c) Temperature is 165 °C.

Figure 2 .
Figure 2. The relationship between mean square displacement (MSD) of four asphalt conponents on the surface of Al 2 O 3 and time: (a) Temperature is 25 • C; (b) Temperature is 65 • C; (c) Temperature is 165 • C.

Figure 4 .
Figure 4. Diffusion coefficient of four asphalt components on the surface of Al2O3 at different temperatures.

Figure 3
Figure3indicates that the diffusion speeds of asphaltene, resin and aromatic on surface of Al 2 O 3 were not sensitive to temperature.For saturate, the diffusion speed at 165 • C was obviously higher than the diffusion speeds at 25 • C and 65 • C.This is due to that the higher temperature can result in the faster diffusion.The slope of curves in Figures2 and 3can be obtained by fitting the relationship between MSD and time.The diffusion coefficients of four asphalt components at different temperature can be calculated according to Equation(6).The results are shown in Figure4.

Figure 4 .
Figure 4. Diffusion coefficient of four asphalt components on the surface of Al2O3 at different temperatures.Figure 4. Diffusion coefficient of four asphalt components on the surface of Al 2 O 3 at different temperatures.

Figure 4 .
Figure 4. Diffusion coefficient of four asphalt components on the surface of Al2O3 at different temperatures.Figure 4. Diffusion coefficient of four asphalt components on the surface of Al 2 O 3 at different temperatures.
is the potential energy function of bond length; ∑ E coulomb is the Interaction energy caused by coulomb electrostatic force (kcal/mol); E vdw is the interaction energy caused by Van der Waals force (kcal/mol).E coulomb and E vdw can be determined by Equations (

Table 1 .
The physical parameters of crystal structure model of different mineral.