Modulus simulation of asphalt binder models using Molecular Dynamics (MD) method

Construction and Building Materials 162 (2018) 430–441

Contents lists available at ScienceDirect

Construction and Building Materials journal homepage: www.elsevier.com/locate/conbuildmat

Modulus simulation of asphalt binder models using Molecular Dynamics (MD) method Hui Yao a, Qingli Dai a, Zhanping You a,⇑, Andreas Bick b, Min Wang c a

Department of Civil and Environmental Engineering, Michigan Technological University, Houghton, MI 49931, USA Scienomics SARL, 16 rue de l’Arcade, 75008 Paris, France c Department of Mathematical Sciences, Michigan Technological University, Houghton, MI 49931, USA b

h i g h l i g h t s
Simulate asphalt binder modified with xGNP using the MD method.
Analyze different modulus properties of these asphalt binder models.
Temperature-modulus trends of models were comparable to those of the laboratory data.

a r t i c l e

i n f o

Article history: Received 4 December 2016 Received in revised form 8 September 2017 Accepted 18 September 2017

Keywords: Asphalt Molecular Dynamics (MD) Bulk modulus Young’s modulus Shear modulus

a b s t r a c t The objectives of this study are, 1) to simulate asphalt binder modified with exfoliated multi-layered graphite nanoplatelets (xGNP) using the Molecular Dynamics (MD) method, and 2) to analyze different modulus properties of these asphalt binder models compared with those of the control asphalt binder model. The multi-layered graphene model was used to represent the xGNP particles, which were used to modify the control asphalt in the laboratory. The three-component control asphalt binder model was used as in the authors’ previous study. The xGNP modified asphalt binder model was built by incorporating the xGNP model and control asphalt binder model and controlling mass ratios to represent the laboratory prepared samples. After the xGNP modified asphalt binder model was generated, the densities of the control and xGNP modified asphalt binder models were computed and verified. Mechanical properties of these models were simulated and calculated in MD simulations using procedures similar to those in the experiments, which include the bulk modulus, Young’s modulus, shear modulus and Poisson’s ratio. The simulation results indicate that the temperature-modulus trends of these asphalt binder models were comparable to those of the laboratory data. The MD simulation data were larger than the laboratory results due to limitations of the current MD simulation, which are discussed in this study. In addition, Poisson’s ratios calculated from the MD simulations coincided with the laboratory results. Ó 2017 Elsevier Ltd. All rights reserved.

1. Introduction Nowadays, asphalt pavement is the main type of pavement used around the world, and the performance, properties, and composition of asphalt are studied by researchers. Asphalt is composed of carbon, hydrogen and some trace elements [1–4]. The performance of asphalt also depends on the source and ambient temperatures. Different modifiers have been applied to the original asphalt in order to obtain better pavement performance. Rubber

⇑ Corresponding author. E-mail addresses: [email protected] (H. Yao), [email protected] (Q. Dai), [email protected] (Z. You), [email protected] (A. Bick), minwang@mtu. edu (M. Wang). https://doi.org/10.1016/j.conbuildmat.2017.09.106 0950-0618/Ó 2017 Elsevier Ltd. All rights reserved.

particles [5–10] and fibers of different sizes [11–13] have been widely utilized to enhance the high-temperature performance of asphalt binders and mixtures and, recently, nanomaterials have been used [12,14–16]. The modification mechanism of these modifiers was also explored and the microstructure of the modified asphalt was examined by researchers [5,14,17–21]. However, the fundamental understanding of the modification and interaction between the modifiers and asphalt are still unclear on a nanoscale in the modified asphalt materials. Molecular dynamics (MD), originating from physics, is widely applied in different research areas. It is a computer program used to simulate the materials and their behaviors on the atomic scale [22]. The physical movements and interactions of atoms and molecules are simulated for a specific period of time. The trajectories of

431

H. Yao et al. / Construction and Building Materials 162 (2018) 430–441

atoms and molecules are recorded based on the solution of Newton’s equation of motion for the interacting particles, where forces and energies between the particles are calculated by the force field [23]. At present, the molecular dynamics (MD) method is a new tool for researchers to understand the physical properties of materials and functions of macromolecules on the atomic scale. MD describes the details of each particle/atom motion over a period of time, and it helps address some specific problems or principles of atoms or materials. It is easier than doing the experiments in the laboratory. In MD simulation, a specific property of materials can be solely studied by altering the specific contributions, and it can be used for fundamental research or functional predictions of a new material. There are many differences between the MD simulations and other numerical tools: 1) inputs in the MD model are descriptions of interatomic and intermolecular interactions, which are different from other numerical tools. Most inputs from other tools are modulus, stiffness, or other parameters; 2) no assumptions are required for processes or mechanisms, whereas other tools may need assumptions for mechanisms; 3) detailed molecular or atomic level information is provided while other tools cannot provide atomic information; 4) the results are more reliable compared to other numerical tools. The asphalt material was simulated with three different components [3] and four different components [24,25]. With three kinds of components, two different asphaltenes, saturates and naphthene aromatics were first used to build the three-component molecular model of asphalt [3]. Different properties of the asphalt binder models, such as density, thermal expansion, isothermal compressibility, and bulk modulus, were calculated and evaluated with the All-atom Optimized Parameters for Liquid Simulation (OPLS-aa) force field [3]. When using four kinds of components, an asphalt binder model, which is similar to the asphalt AAA-1 of the Strategic Highway Research Program (SHRP) [24], was generated with seven components. The composition of this model consists of asphaltene, polar aromatic, naphthene aromatic, and saturate [24]. The Hansen solubility parameters of components in the asphalt binder models were calculated, and the radial distribution function and intramolecular orientation of components were also analyzed [24]. The new compositions of asphalt models were proposed to represent the AAA-1, AAK-1, and AAM-1 asphalt with twelve molecules, and the partial charges and force field parameters were determined through quantum mechanics. The Hansen solubility parameters were used to understand the components in asphalt models. The densities of the new models are better than those of prior models [25]. Based on these references, the feasibility and validity of the MD application in asphalt materials were demonstrated, however, the predicted accuracy and algorithms or methods used for analysis need to be improved. The objectives of this study are as follows: apply the Molecular Dynamics (MD) method to simulate and calculate the moduli of the control and xGNP modified asphalt binder models and compare these values with the laboratory results. The modulus properties include the bulk modulus, Young’s modulus and shear modulus. In this study, the control asphalt binder model was composed of asphaltenes, naphthene aromatics, and saturates, based on the Corbett method [2]. The graphite model with the multi-layer molecular structure was used to represent exfoliated multi-layered graphite nanoplatelet (xGNP) particles as a modifier for this study. The modified asphalt binder model was built by incorporating the xGNP model and control asphalt binder model and controlling the mass ratios to represent the laboratory prepared samples. Mechanical modulus properties of the asphalt binder models were analyzed with different MD simulation methods: 1) compute the bulk modulus of the asphalt binder models through compressive and expansive volumetric strain-pressure relations; 2) simulate the Young’s modulus of these models by applying the

compressive and tensile strain in one direction at different temperatures, and then, calculate stress/strain ratio in that direction; 3) calculate the shear modulus of the models by applying a shear strain in the XY direction at different temperatures, and then, compute shear stress/strain ratios; and 4) analyze Poisson’s ratio of the models using the relationship between Poisson’s ratio and the moduli (bulk modulus and shear modulus). 2. Force field The force field, a kind of energy function, provides the parameters for every type of atom in the MD simulation to calculate the energy of the system. The parameters of the energy functions can be derived from the physical and chemical experiments, as well as from the calculation results of quantum mechanics. The basic functional form of energy in MD simulations includes bonded and nonbonded terms. The bonded terms for interactions are described by covalent bonds, and nonbonded terms are represented by long-range electrostatic and van der Waals forces. The Amber Cornell Extension Force Field (ACEFF) was used based on the Amber Cornell Force filed [26] in this study, and some experimental parameters were obtained from the General Amber Force Field (GAFF) [27]. The formula for energy calculation of ACEFF is shown in Eq. (1). The NWChem [28] was used to optimize the geometry and calculate the electrostatic potential (ESP) charge for each component of the asphalt binder model. When assigning the ESP charge to atoms, the density functional theory (DFT) module was adopted for optimization and charge fit.

Etotal ¼

X bonds

þ

K r ðr  req Þ2 þ

X

K h ðh  heq Þ2

angles

" # X Vn X Aij Bij qi qj  þ ½1 þ cosðnu  cÞ þ 12 2 R6ij Rij i
ð1Þ

where r eq or heq is the equilibrium of the structural data; K r is the force coefficient between double bonds; K h is the force coefficient from vibrational analysis; n is the multiplicity for dihedrals; c is the phase angle for the torsional angle parameters; A, B and q are the non-bonded potentials; qi and qj are charges of two atoms, and the charges were calculated and obtained from NWChem analysis; Rij is the distance between atoms; and  is the depth of the van der Waals energy well. 3. Energy optimization and asphalt binder models 3.1. Model optimization methods The xGNP modified asphalt binder model (Fig. 1b) was optimized by the conjugate gradient method (a kind of algorithm to solve optimization problems) and the Particle-Particle-ParticleMesh Method (PPPM) (an accurate and efficient method to calculate interactions in MD simulation) [29]. The conjugate gradient method was used to relocate the positions of atoms and reallocate the velocities of atoms in asphalt binder models [30,31]. Based on the equation and calculation of force fields, the positions of atom relates to the energy calculation and the velocities of atoms relate to the moment of the atoms. The conjugate gradient method was used to optimize the positions and velocities for energy and parameter optimizations. The conjugate gradient method is often used as an iterative method (Eq. (2)) for large systems to resolve unconstrained optimization problems, such as energy optimizations in MD simulations, which cannot be optimized by direct methods. The PPPM (a kind of Ewald summation) was used to optimize the system energy by mapping charges to a threedimensional mesh and using three-dimensional fast Fourier

432

H. Yao et al. / Construction and Building Materials 162 (2018) 430–441

xGNP modified asphalt model Exfoliated Graphite Nanoplatelets xGNP model (Layer Distance: 3.35 Ǻ)

W

H L

b

a

Fig. 1. The molecular model of multi-layer graphene sheets and xGNP modified asphalt model, (a) the molecular model of exfoliated graphite nanoplatelets (xGNP model); (b) xGNP modified asphalt model. (Simulation box dimension: 66.5 Å (Length)  69.1 Å (Width)  34.7 Å (Height)).

transforms (FFTs) to solve Poisson’s equation on the mesh, then interpolating the electric fields back to atoms [29]. In PPPM calculation, the interparticle forces can be divided into two parts: short-range and slowly varying (Eq. (3)) [29].

f ðxÞ ¼

1 T x Ax  xT b; x 2 Rn ðiterative methodÞ 2

ð2Þ

where A is symmetric, positive and real; b is a known coefficient; and vectors u and v are non-zero. m F ij ¼ F sr ij þ F ij

ð3Þ

where F ij is the interparticle forces in the system; F sr ij is the rapidly varying short-range component; F m ij component.

is the slowly varying

3.2. Control asphalt binder model The control asphalt binder model was composed of three components, asphaltene, naphthene aromatic, and saturate at a ratio of 5:27:41 [3]. The structure of the asphaltene (C64H52S2) was introduced from the references [1,3,32]. Based on the similarity of properties of docosane and saturate, docosane (nAC22H46) [3,33] was used to represent the saturate. The 1, 7- dimethylnaphthalene (C12H12) [3,34] was used to represent the naphthene aromatic due to the similar ratio of alkane and aromatic. The control asphalt binder model with ACEFF and ESP charge was built at a low start density of 0.1 g/cm3. The initial positions of molecules were randomly placed in the simulation cell with the initial configurations from references. Energy minimization (conjugate gradient method) was conducted on molecules and binder models with initial configurations. The velocities of molecules or atoms were randomly assigned from a Maxwell-Boltzmann distribution with no net angular momentum at an initial state. The target density was around 1.0 g/cm3. The NPT (isothermal-isobaric ensemble) was applied to the control asphalt binder model to compress the model for geometry and energy optimizations. The composition of the control asphalt binder model is shown in Table 1. 3.3. xGNP modified asphalt binder model In the laboratory, the xGNP graphene nanoplatelets were used to improve the performance of the asphalt binder, and the PG

grade of the control asphalt binder was PG 58–28. The 2% xGNP particles were slowly added to the control asphalt binder, and the modified asphalt binder was mixed in the high shear machine for around two hours at a temperature of 450 K. In MD simulation, the model of multi-layer graphene sheets (Fig. 1a) was used to represent the xGNP particles (nanoparticle size is relatively small), and this xGNP model was added to the control asphalt binder model at the temperature of 450 K and under 1 atm of pressure. The 2% (by mass) xGNP model was dispersed in the matrix of the control asphalt binder model, and more molecules of the modified asphalt binder model were used to maintain the mass ratio of the xGNP model in the control asphalt binder model. The percentage of the xGNP model in the modified asphalt binder model is based on the mass of the control asphalt binder model. It is commonly used in the references [35–37]. The total mass of the previous control model was 21378.81 g, and the mass of the xGNP model was 1890.13 g. Mass mentioned here is based on 1 mol of the simulation box. The mass of the control model expanded to 85515.24 g (four times the previous one) in order to fully satisfy the required percentage of xGNP particles in the control model. The 2% xGNP particles were used in the modification of asphalt for laboratory tests. The mass of the xGNP modified asphalt binder model was around 87405.39 g. It is also beneficial to analyze the variation in data due to mass in MD simulations. The composition of the xGNP modified asphalt binder model is shown in Table 2. 4. Density calculation of control and modified asphalt binder models Before the modulus calculations, the densities of the asphalt binder models were computed and verified (Fig. 2). The NPT (isothermal-isobaric ensemble) simulation was employed to compress and relax the asphalt binder models. The Savitzky-Golay filter was employed to smooth the data trend. It uses linear least squares and low-degree polynomial regression to fit the successive subsets of adjacent data points. Based on the fitting trends, the densities of the control and xGNP modified asphalt binder models with ACEFF were around 0.903 g/cm3 and 0.911 g/cm3, respectively. The densities of these asphalt binder models were stable after around 1 ns. The addition of xGNP model in the control asphalt binder model increased the density of the asphalt binder model. It is likely that more molecules were added in the model

433

H. Yao et al. / Construction and Building Materials 162 (2018) 430–441 Table 1 The composition of the control asphalt binder model. Control asphalt binder model

Mass (g/mol)

Sum Formula

Number of molecules

Total Mass (g)

Mass fraction (%)

Asphaltene 1, 7- Dimethylnaphthalene Docosane Control asphalt binder model

885.23 156.22 310.6

C64H52S2 C12H12 C22H46

5 27 41 73

4426.15 4217.94 12734.6 21378.81

20.685 19.734 59.581

Table 2 The composition of the xGNP modified asphalt binder model. xGNP modified asphalt binder model

Mass (g/mol)

Sum Formula

Number of molecules

Total Mass (g)

Mass fraction (%)

Asphaltene 1, 7- Dimethylnaphthalene Docosane xGNP xGNP modified asphalt binder model

885.23 156.22 310.6 472.53

C64H52S2 C12H12 C22H46 C38H16

20 108 164 4 296

17704.59 16872.16 50938.50 1890.13 87405.39

20.25 19.30 58.30 2.15

Density 0.93

0.92

0.91

Density(g/cm3)

0.9

0.89

0.88

0.87 Density of xGNP modified asphalt Fitted density curve for xGNP modified asphalt Density of control asphalt Fitted density curve for control asphalt

0.86

0.85

1

2

3

4

5

6

7

8

Step

9

10 10

5

Fig. 2. Densities of the control and xGNP modified asphalt models.

without chemical reactions, improving the density of the modified model.

5. MD simulation methods for mechanical modulus prediction Different moduli were calculated for the asphalt binder models using MD simulations, such as the bulk modulus, Young’s modulus and shear modulus. The strain or displacement controlled testing methods were applied in the MD simulation of the models. After energy optimizations and relaxation of the NPT ensemble, the modulus calculations were carried out, and the results of MD simulations were as follows. Based on the references regarding the relationship between stress and strain [38,39], it is impossible to reproduce the laboratory strain rates through MD simulations under current computational resources. It is difficult to realize the experimental deformation mechanisms using an appropriate rate in MD simulations. Even though the MD simulation data are close to the experimental results, different deformation mechanisms are observed under different strain rates in MD simulations.

5.1. Bulk modulus The bulk modulus, a kind of numerical constant, is a measure of the ability of a material to withstand the changes in the volume under compression or tension in three directions. The bulk modulus is defined by Eq. (2). In the laboratory, there are four approaches to test the bulk modulus: two separate tests (uniaxial and shear), a long cylindrical test in tension and compression, confined compression test, and poker-chip geometry in tension and compression. There are two methods to calculate the bulk modulus in MD simulation, one is to directly use the applied strain and responding stress (the first calculation method), and another one is to use the regression function to compute the bulk modulus (the second regression method). The average volumetric stressstrain relationships (the first method) in these asphalt binder models were applied in this study. In the MD simulation, the small volumetric strains were applied to the control and xGNP modified asphalt binder models, and the compressive and expansive volumetric strains were applied using Eq. (3) [40]. The Nose/Hoover thermostat [41,42] and SLLOD algo-

434

H. Yao et al. / Construction and Building Materials 162 (2018) 430–441

rithm (Lee-Edwards boundary conditions) [43] were employed for these asphalt binder models when performing NVT (canonical ensemble) simulations in these models. NVT, a kind of statistical ensemble, represents the probable distribution of states in a mechanical system at thermal equilibrium. The probability P is assigned to each microscopic configuration through the following FE

equation P ¼ e kT (E is the total energy of the configuration of the system, k is Boltzmann’s constant, T is the temperature of the system, and F is the free energy). The positions and velocities of atoms were updated during the simulation time, and pressure changes were also updated and recorded. These models were fixed by NVT simulations for at least 2 ns with a simulation step of 1 fs. The pressure data over each simulation step were recorded. The applied compressive and expansive volumetric strains were calculated with changes in the simulation box edge length (shown in Eq. (4)). During the MD simulations, the compressive and expansive volumetric stress were calculated using Eq. (5). The average values of the bulk modulus for these asphalt binder models were computed by averaging the stress/strain ratio (Eq. (6)) over simulation time. The MD simulation was conducted at different temperatures to obtain temperature-modulus trends.

K ¼ V

dP dV

ð4Þ

exx ¼ eyy ¼ ezz ¼ 0:005

ð5Þ

D ¼ ð1  exx Þð1  eyy Þð1  ezz Þ  1

ð6Þ

1 3

rh ¼ ðrxx þ ryy þ rzz Þ K¼

rh D

ð7Þ ð8Þ

where P is the pressure, V is the volume, exx , eyy , and ezz are the small strains in all three coordinate directions; rxx , ryy , and rzz are the volume averaged stresses in all three coordinate directions; D is the overall dilatation of the simulation box; rh is the averaged stress in the MD simulation; and K is the bulk modulus in the MD system. After a few ns NVT simulations, the corresponding stress results in three directions of these asphalt binder models were calculated, and one of the stress results (under different temperatures, three directions, and compressive and expansive strains) was presented. The stress results of the control and xGNP modified asphalt binder models under compressive strain at the temperature of 298.15 K are shown in Fig. 3a. The average volumetric stress-strain relationships in these asphalt binder models are shown in Fig. 3b based on the second method. Consequently, the bulk moduli of these asphalt binder models at different temperatures are shown in Fig. 3c after averaging based on the first method. The displacement controlled simulation was applied to represent the laboratory nano-indentation tests. Fig. 3a shows the stress magnitudes of the control and xGNP modified asphalt binder models under compressive strain at a temperature of 298.15 K. The data amplitudes of asphalt binder models in different directions are different and relatively stable around ±108 Pa. The data of the xGNP modified asphalt binder model fluctuates less than that of the control asphalt binder model due to different numbers of molecules or atoms. Fig. 3b shows the relationships between volumetric stresses and strains of these asphalt binder models and their regression lines for the asphalt binder models under compressive strain at the temperature of 298.15 K after averaging every 100 steps. The linear regression is fit for a stress-strain relationship in the modulus calculations. The changes in volumetric strains range from 0 to 0.015 in two asphalt binder models, and they are consistent with the applied strains. The predicted magnitude of the bulk

modulus with applied compressive or expansive strain in this MD simulation is consistent, and it is around 109 Pa. Fig. 3c displays the bulk moduli of the control and xGNP modified asphalt binder models at different temperatures based on the first method. The moduli of the control and xGNP modified asphalt binder models decrease by increasing the system temperatures. The moduli of the xGNP modified asphalt binder model are higher than those of the control asphalt binder model at different temperatures. The addition of the xGNP model in the control asphalt binder model increases the bulk modulus. Based on the reference [44], the bulk modulus of the asphalt binder was tested in the laboratory using the poker-chip test, and it is around 3 MPa at the temperature of 293 K. There is still a difference between the simulation and laboratory data. In this MD simulation, the actual strain was applied to the asphalt binder model, and it is similar to the laboratory test procedure. The time frame of the applied strain is around 1 fs, and it is really short compared to physical experimental tests. The simulation time cannot fully account for the time-dependent asphalt material properties. The three components asphalt binder model may not represent the complex asphalt materials. Also, a limited number of molecules could be another reason for this simulation data. These will lead to the difference between the predicted and laboratory tested bulk modulus. However, the presented MD simulation method provided a foundation for the nanoscale property prediction. These simulations can be further improved by addressing some of these limitations for closer prediction. 5.2. Young’s modulus Young’s modulus is a stiffness measurement for elastic materials and it is the ratio of the responding stress to the (responding) strain along an axis under Hooke’s law [45]. The Young’s modulus is different from the bulk modulus and shear modulus, sometimes called the modulus of elasticity. The asphalt can be treated as elastic material under low temperatures because most materials exhibit nearly Hookean characteristics for a small strain or stress. Considering that, the test temperature in this MD simulation is lower than room temperature. In the laboratory, the nanoindentation test was used to evaluate the Young’s modulus of asphalt. Both the loading and depth control modes can be used in the test. Considering the influence of the indenter tip, the reduced Young’s modulus was introduced, and the Young’s modulus Es of the sample was calculated using Eq. (7). The reduced modulus was calculated from the slope of the load-indentation curve when unloading using Eq. (8) [46]. The Young’s moduli of PG 70– 22 and PG 76–28 asphalt binders were around 0.76 MPa and 5.22 MPa at room temperature, respectively [46].

1 1  v 2s 1  v 2i ¼  Er Es Ei

ð9Þ

pffiffiffiffi Er ¼

p S

2

pffiffiffi A

ð10Þ

In the above equations, Er is the reduced modulus; Es is the Young’s modulus of the tested sample; Ei is the Young’s modulus of the indenter tip; v s is Poisson’s ratio of the tested sample; v i is Poisson’s ratio of the indenter tip; S is the initial unloading stiffness, and A is the contact area during the test. During the MD simulation, small strains (expansive and compressive 5‰) were applied to these asphalt binder models in one direction. The NPH (isenthalpic-isobaric ensemble) simulation, which maintains constant enthalpy, pressure and number of particles in the systems, was performed on the asphalt models with temperature control (temp/rescale). The enthalpy (H) of a system can be defined as H = E + PV: E is the internal energy of the system,

435

H. Yao et al. / Construction and Building Materials 162 (2018) 430–441

8

10 8

10 8

5 4

Average Volumetric Stress(Pa)

6 4

Stress(Pa)

2 0 -2 Stress in X direction-control asphalt model Stress in Y direction-control asphalt model Stress in Z direction-control asphalt model Stress in X direction-xGNP modified asphalt model Stress in Y direction-xGNP modified asphalt model Stress in Z direction-xGNP modified asphalt model

-4 -6 -8 280

Average Volumetric Stress vs. Average Volumetric Strain Average volumetric stress vs. average volumetric strain-control asphalt model Linear regression-control asphalt model Average volumetric stress vs. average volumetric strain-xGNP modified asphalt model Linear regression-xGNP modified asphalt model

3 2 1 0 -1 -2 -3

for control asphalt model: y=1.609 e9*x+1.879 e7 for xGNP modified asphalt model: y=2.191 e9*x+6.436 e6

-4 0

0.002

0.004

0.006

0.008

0.01

0.012

0.014

0.016

Average Volumetric Strain 285

290

295

300

305

310

315

Temperature (K)

(a) Stress level of the control and xGNP modified asphalt models under compressive strain at the temperature of 298.15 K (1.2×105data)

3.5

10

9

(b) Average volumetric stress-strain relationship between the control and xGNP modified asphalt models under compressive strain (Strain rate in each direction: 0.005fs-1) at the temperature of 298.15 K. Note: The two regression lines are very close, one is a black dotted line for the control asphalt model, and the other is a red dashed line for the xGNP modified asphalt model. The regression function is designed to determine the bulk moduli of these asphalt models.

Bulk Modulus vs. Temperature

Bulk Modulus (Pa)

3

2.5

xGNP modified asphalt model Linear regression-xGNP modified asphalt model Control asphalt model

2

Linear regression-control asphalt model

1.5

1 265

270

275

280

285

290

295

300

Temperature (K)

(c) Bulk moduli of the control and xGNP modified asphalt models at different temperatures Fig. 3. The modulus calculations for the control and xGNP modified asphalt.

which contains kinetic and potential energies; P is the pressure of the system; and V is the volume of the system. The enthalpy change is equal to the energy flow of heat at a constant pressure [47] (DH = DE + PDV = qp, qp is the heat at a constant pressure). Meanwhile, the uniaxial tensile or compressive stress was applied in one direction and no stress was applied in the other directions. The position and velocities of atoms in the systems were updated at different temperatures, and the pressure or stress was recorded during the application of strain. The MD simulation time was at least 1 ns, and the output data was averaged every 100 timesteps (simulation timestep, 1 fs). The averaged pressures in the control and xGNP modified asphalt binder models were calculated under the expansive and compressive strains in independently different directions at different temperatures. The pressure results of the control and xGNP modified asphalt binder models in MD simulations are shown in Fig. 4a (under expansive strain) and Fig. 4b (under compressive strain). In Fig. 4a, the label

‘‘Stress in X direction – control asphalt model at 298.15 K” represents the stress in the control asphalt model when the expansive strain was applied in the x direction at the temperature of 298.15 K, and no strain or stress in the other directions. In Fig. 4b, the label ‘‘Stress in X direction – control asphalt model at 298.15 K” represents the stress in the control asphalt model when the compressive strain was applied in the x direction at the temperature of 298.15 K, and no strain or stress in the other directions. The relationship between the applied strain and stress in these models and the linear regression are shown in Fig. 4c (under expansive strain) and Fig. 4d (under compressive strain) based on the second regression method, and the data was regressed and calculated by MATLAB. In Fig. 4c, the two regression lines are close, one is a black dotted line for the control asphalt model and the other is a red dashed line for the xGNP modified asphalt model. The regression analysis was used to calculate the Young’s modulus of these asphalt models under expansive strain. In Fig. 4d, the two regres-

436

H. Yao et al. / Construction and Building Materials 162 (2018) 430–441

6

Stress

10 8

Stress in X direction-control asphalt model at 298.15K Stress in Y direction-control asphalt model at 298.15K Stress in Z direction-control asphalt model at 298.15K Stress in X direction-xGNP modified asphalt model at 298.15K Stress in Y direction-xGNP modified asphalt model at 298.15K Stress in Z direction-xGNP modified asphalt model at 298.15K

4

Stress (Pa)

2

0

-2

-4

-6 285

290

295

300

305

310

Temperature (K)

(a) Stresses of the control and xGNP modified asphalt models under expansive strain 6

Stress

10 8

Stress vs. Strain

108

6

Linear regression-control asphalt model: y = -2.291 e9 *x - 9.383 e6

4

2

2

Stress (Pa)

Stress (Pa)

Linear regression-xGNP modified asphalt model: y = 1.02 e9 *x - 3.289 e6

4

0

-2

-4

-2

-4 Stress in X direction-control asphalt model at 298.15K Stress in Y direction-control asphalt model at 298.15K Stress in Z direction-control asphalt model at 298.15K Stress in X direction-xGNP modified asphalt model at 298.15K Stress in Y direction-xGNP modified asphalt model at 298.15K Stress in Z direction-xGNP modified asphalt model at 298.15K

-6

-8 285

0

Stress vs. Strain-control asphalt model Linear regression-control asphalt model Stress vs. Strain-xGNP modified asphalt model Linear regression-xGNP modified asphalt model

-6

-8 290

295

300

305

310

-5

-6

-4

Temperature (K)

(b) Stresses of the control and xGNP modified asphalt models under compressive strain

-2

-3

-1

0

Strain

10

-3

(d) Stress-strain relationship of the control and xGNP modified asphalt models under compressive strain in the z direction at the temperature of 298.15 K and the linear regression.

4

10

Youngs moudulus vs. Temperature

9

Youngs modulus (Pa)

3.5

3

MD results-control asphalt model Linear regression-control asphalt model MD results-xGNP modified asphalt model Linear regression-xGNP modified asphalt model

2.5

2

1.5

1

0.5 265

270

275

280

285

290

295

300

Temperature (K)

(c) Stress-strain relationship of the control and xGNP modified asphalt -1 models under expansive strain (strain rate: 0.005fs ) in the z direction at the temperature of 298.15 K, and the linear regression.

(e) Young’s modulus versus temperature of the control and xGNP modified asphalt models

Fig. 4. Young’s modulus calculation of the control and xGNP modified asphalt models.

sion lines are close, one is a black dotted line for the control asphalt model and the other is a red dashed line for the xGNP modified asphalt model. The regression analysis is employed to compute the Young’s modulus of these asphalt models under compressive

strain. The Young’s moduli of these models at different temperatures are shown in Fig. 4e. Fig. 4a shows the stresses of these asphalt binder models in six different simulations at the temperature of 298.15 K in three dif-

437

H. Yao et al. / Construction and Building Materials 162 (2018) 430–441

direction since it is hard to relax the stress in a largely confined system. 5.3. Shear modulus The shear modulus is designed to describe the response of the materials to shear stress and measure the stiffness of materials based on the generalized Hooke’s law. The shear modulus can be calculated from Eq. (12). The shear modulus of asphalt is related to the field performance of asphalt pavement. In the laboratory, the shear modulus of the asphalt binders is tested by a Modular Compact Rheometer (MCR) (Fig. 5a) with different plates (Fig. 5b) for different temperatures, and the temperatures can range from 253 K to 355 K. The round plate with a diameter of 25 mm was selected to test the complex shear modulus using the oscillatory mode. The testing data of the control (PG 58–28) and the modified asphalt binder with 2% xGNP particles under the condition of the temperatures (319.15 K–343.15 K) and frequency of 45 Hz are shown in Fig. 5c. In the MD simulation, the small shear strains (compressive and expansive 5‰) were applied to the asphalt binder models in the XY direction, and it was similar to the procedure used in the laboratory test. The applied strain between the laboratory and simulation was different, because the asphalt binder model was assumed to be a purely elastic material in this study. The NVT ensemble was used with the SLLOD algorithm [43] when the shear strain (Eq. (12)) was applied to the models. The positions and velocities of atoms in the asphalt binder models were updated every step, and the pressures and temperatures were recorded. The stress results of these asphalt binder models under expansive and compressive strains are shown in Fig. 6a and b. The shear moduli of these asphalt binder models are shown in Fig. 6c based on the first calculation method.

b

a

(a) Modular Compact Rheometer (MCR); (b)Accessories for testing plates with different diameters and temperatures 1.0E+06

log Shear Modulus (Pa)

ferent directions when the expansive strain was applied. Fig. 4b displays the stresses of these asphalt binder models in the six different simulations at the temperature of 298.15 K in three different directions when the compressive strain was applied. It is noted that the stresses of the xGNP modified asphalt binder model have less fluctuations and more concentrated data than those of the control asphalt binder model due to the larger number of molecules and atoms. The stresses from different directions (x, y, and z) in the control asphalt binder model are of the same magnitude and range. It indicates that there is no difference in stresses among three different directions when applying the same strain to the asphalt binder model. The same situation occurs in the xGNP modified asphalt binder model. Fig. 4c shows the relationship between the tensile stress and strain of the control (green squares) and xGNP (blue circles) asphalt binder models in the z direction at the temperature of 298.15 K under expansive strain. Fig. 4d displays the relationship between the stress and strain of these models in the z direction at the temperature of 298.15 K under compressive strain. The amplitude of variation in the stress-strain data of the xGNP modified asphalt binder model (blue circles) is less than that of the control asphalt model (green squares), and the numbers of atoms and molecules may definitely influence the variation in results. The regression function in Fig. 4c and d was also used to calculate the Young’s modulus of the control (a black dotted line) and xGNP (a red dashed line) asphalt binder models based on the linear relationship between stress and strain. The linearity of stress-strain data was tested in the statistical section of this study. Fig. 4e demonstrates the Young’s modulus of the control (a red solid line) and xGNP (a blue dotted line) asphalt binder models at different temperatures. The Young’s modulus of these models was averaged from the results of compressive or expansive strain in different directions. The Young’s modulus of these models decreased with increasing ambient temperatures in the models. Compared to the laboratory testing results (0.76 MPa for PG 70–22 asphalt and 5.22 MPa for PG 76–28 asphalt) of Young’s modulus at room temperature [46], there is a wide gap between the experimental and simulation data (around 109 Pa). The compressive or expansive strain was applied in these asphalt binder models for a short time, a few thousand time steps. It is reasonable that the small deformation at a very short time results in a large stress in an assumed elastic model at low temperatures. Meanwhile, the simulation or interaction time and numbers of atoms or molar mass in these models are also the determining factors for the calculation of the Young’s modulus. It is obvious that the Young’s moduli of the xGNP modified asphalt binder model are lower than those of the control models at different temperatures, and the Young’s modulus of the xGNP modified asphalt tends to approach the experimental results relative to the control asphalt model. It is probably due to a large number of atoms and molecules in the xGNP modified asphalt binder model compared to that in the control asphalt model. More molecules could bring more relaxation possibilities, from space, conformations or configurations of models, as well as similar deformation mechanisms in the laboratory. Also, the NPH ensemble was used for this simulation, the volume and one or two dimensions of simulation box changed (DH = DE + DPV), as well as internal energy, when the strain was applied during MD simulations under constant enthalpy. The internal energy relates to the conformation of the MD system. When compressive or expansive strain was applied for a short time, more atoms or molecular numbers in the xGNP modified asphalt binder model helped relax the responding stress caused by the small strain applied in one direction when the other directions were not confined. If the simulation cell is confined or the strains are applied in three directions, the responding stress would be higher than the stress on the condition of no confinement or strain in one

1.0E+05

y = 3E+21e-0.117x R² = 0.9993

y = 5E+21e-0.12x R² = 0.9996 1.0E+04

Shear modulus-xGNP modified asphalt binder Shear modulus-Control asphalt binder Expon. (Shear modulus-xGNP modified asphalt binder) Expon. (Shear modulus-Control asphalt binder)

1.0E+03 315

320

325

330

335

340

345

Temperature (K)

(c) Shear moduli of the control (PG 58-28) and modified asphalt binders with 2% xGNP particles Fig. 5. The modular compact rheometer (MCR) for the shear modulus of asphalt binders.

438

G¼

H. Yao et al. / Construction and Building Materials 162 (2018) 430–441

sxy and cxy ¼ 0:005 cxy

of the control and xGNP modified asphalt binder models are around 1.4 atm and 15 atm, respectively, due to the applied strain. The pressure in these MD systems is relatively concentrated, and median of the temperature is around 337.15 K. It is also observed in Fig. 6b where stresses of these models under compressive strain in the XY direction are at a temperature of 337.15 K. Fig.6c shows the shear moduli of the control and xGNP modified asphalt binder models under different temperatures after averaging. In these MD simulations, the simulation data of the shear modulus of these asphalt binder models decreases with the increase in test temperature, and the exponential regression is well fitted to the simulation results. The trends of the shear modulus of these asphalt binder models agree well with those of the experimental data. Additionally, the magnitude of the shear modulus is around 108 Pa, and it is larger than the laboratory results. The shear modulus of the control asphalt binder model in the MD simulation is also higher than that of the xGNP modified asphalt binder model. There are several reasons for the large shear modulus: (1) the

ð11Þ

where G is the shear modulus of the sample; sxy is the shear stress in the XY direction; and cxy is the shear strain in the XY direction. Fig. 5c shows the shear moduli of the control and xGNP modified asphalt binders at different temperatures. The laboratory data of the shear modulus was regressed and fitted by the exponential curve, and the value of R squared is around 0.99. The shear modulus of the control and xGNP modified asphalt binder decreases with increasing test temperatures. The magnitude of the shear modulus of unaged asphalt is around 105 Pa. However, the magnitude of the shear modulus of aged asphalt can reach more than about 107 Pa [15]. Fig. 6a displays the stresses of the control and xGNP modified asphalt models under expansive shear strain in the XY direction at the temperature of 337.15 K. The variation in the control asphalt binder model is larger than that of the xGNP modified asphalt binder, and it is due to fewer molecules/atoms. The average pressures

5 4

10

Stress

8

4 Stress - control asphalt model (positive strain in XY direction) Stress - xGNP modified asphalt model (positive strain in XY direction)

Stress

10 8

3

3 2 2

Stress(Pa)

Stress(Pa)

1 1 0 -1

0

-1

-2

-2 -3

-3

-4

-4

-5 315

320

325

330

335

340

345

350

355

360

-5 315

Stress - control asphalt model (negative strain in XY direction) Stress - xGNP modified asphalt model (negative strain in XY direction)

320

325

330

Temperature (K)

(a) Stresses of the control and xGNP modified asphalt models under positive shear strain in the XY direction at a temperature of 337.15 K

10 8

Shear Modulus - control asphalt model (XY Direction) Shear Modulus - xGNP modifed asphalt model (XY Direction) Exponential Regression - Control asphalt model Exponential Regression - xGNP modified asphalt model

10 7 270

280

340

345

350

355

(b) Stresses of the control and xGNP modified asphalt models under negative shear strain in the XY direction at a temperature of 337.15 K

Shear Modulus vs. Temperature

10 9

Shear Modulus (Pa)

335

Temperature (K)

290

300

310

320

330

340

350

Temperature (K)

(c) Shear moduli of the control and xGNP modified asphalt models under different temperatures (log scale in y-axis) Fig. 6. Shear modulus calculation of the control and xGNP modified asphalt models.

H. Yao et al. / Construction and Building Materials 162 (2018) 430–441

interaction time for applying the strain to these models is really short, several femtoseconds, and the running time of these simulations, a few nanoseconds, is also very short compared to the realtime scale in laboratory tests. It is a reasonable response for elastic materials in asphalt binder models. (2) Asphalt is a viscoelastic material at medium and high temperatures. The response of the asphalt models is elastic, and the results tend to approach the laboratory results when the test temperature is relatively low. (3) The molecular numbers of these asphalt models are limited, and from the trend, the more molecules or atoms in the model, the closer the results are to the laboratory data. More molecules could provide more conformation and configurations for the models, as well as the same deformation mechanism as that in the laboratory. (4) It is possible that the NVT ensemble with SLLOD algorithm needs to be improved for the case of the asphalt binder. In addition, when the strain was applied to these asphalt models, the xGNP modified asphalt binder model expanded and relaxed the strain and responding stress effects due to more molecules and atoms compared to the control asphalt binder model. This is also the explanation as to why the shear modulus of xGNP modified asphalt is less than that of the control asphalt model. It indicates that if there are more molecules/atoms in the model, the modulus can approach the experimental results.

When materials are compressed in one direction, the materials normally expand in the other two directions; this phenomenon is usually called the ‘‘Poisson effect”. Poisson’s ratio is used to measure the Poisson effect, and it can be calculated using Eq. (13). The Poisson’s ratio of isotropic and linear elastic materials ranges from 1.0 to 0.5 due to the expansive values of these moduli. Poisson’s ratio can also be calculated from the bulk modulus and shear modulus after determining these values in the asphalt binder models, and the formula for Poisson’s ratio is shown in Eq. (14) [40,44]. The results of Poisson’s ratio of these asphalt binder models are shown in Fig. 7.

v

detrans deaxial

ð12Þ

3K  2G ¼ 6K þ 2G

ð13Þ

where, v is Poisson’s ratio; detrans is the transverse strain; deaxial is the axial strain; K is the bulk modulus and G is the shear modulus.

Poissons ratio vs. Temperature 0.5

Poissons ratio

0.45

0.4

0.35

0.3

0.25 265

Poissons ratio - control asphalt model Linear regression - control asphalt model Poissons ratio - xGNP modified asphalt model Linear regression - xGNP modified asphalt model

270

Fig. 7 shows the results of Poisson’s ratio of the control and xGNP modified asphalt binder models at a temperature range from 268.15 K to 298.15 K. Poisson’s ratio of the control asphalt binder model is around 0.35, and Poisson’s ratio of the xGNP modified asphalt model is around 0.45 at temperatures from 268.15 K to 298.15 K. It falls within the area of theoretical values. The variation of Poisson’s ratio of the control asphalt binder model is larger than that of the xGNP modified asphalt, and it is likely due to the fewer number of molecules or atoms in the control asphalt model. Considering different loading times and temperatures, the Poisson’s ratio of a typical asphalt binder was from 0.35 to 0.5 based on the reference [48]. Recently, the Poisson’s ratio of two types of asphalt binders was tested at 293.15 K, and it starts from around 0.25 to 0.5 under different loading times [44]. It can be seen that Poisson’s ratios of the control and xGNP modified asphalt binder models are reasonable and agree well with the reference and experimental data. It indicates that the responses of asphalt models are acceptable when applying strain in these asphalt binder models and calculating moduli in MD simulations. 6. Statistical analysis for confident linear stress-strain relations of MD simulation 6.1. Pearson’s product-moment correlation

5.4. Poisson’s ratio

v ¼

439

275

280

285

290

295

300

Temperature (K) Fig. 7. Poisson’s ratios of the control and xGNP modified asphalt models at different temperatures.

The regression method was used to calculate the Young’s modulus (the data from Fig. 4c), and the correlation analysis was used to measure and interpret the strength of the relationship between two variables (strain and stress). The Pearson’s product-moment correlation coefficients were used for evaluation, whose values ranged from 1 to +1. This correlation is commonly used for measuring the linear relationship between two variables [49]. For example, the Pearson correlation coefficient, ‘r’, between two variables X i and Y i (i = 1, 2, . . ., n) can be computed using Eq. (9). Based on this reference [49], the parameter ‘r’, is a measure of the strength of the linear relationship between two variables. A value of 0 indicates that there is no association between the two variables, and a value less than 0 indicates a negative association; ‘r’ has a t distribution with n-2 degrees of freedom.

Pn   i¼1 ðX i  XÞðY i  YÞ ffi r ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi P Pn  2 n ðY i  YÞ  2 i¼1 ðX i  XÞ i¼1

ð14Þ

 and Y  are the sample means of the X i and Y i values, respecwhere X tively; X i and Y i are the stress and strain, respectively, and n is the sample size. The statistics t-test was used for the population correlation coefficient, q, to check whether or not a linear relationship exists between the two variables. More specifically, we tested the null hypothesis, Eq. (10). The Pearson correlation coefficient, r, between strain and stress was obtained from Eq. (9), r = 0.02868321. To test the above hypothesis, the correlation test statistic was equivalent to Eq. (11).

H0 : q ¼ 0 against Ha : q–0

ð15Þ

pffiffiffiffiffiffiffiffiffiffiffiffi r n2 t ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffi ¼ 2:8692 1  r2

ð16Þ

Based on the statistical results, the corresponding p-value from the t-statistic is 0.004124 (<0.01) with degrees of freedom df = 9998 comparing the test statistic to a t-distribution with 9998 degrees of freedom. We can thus reject the null hypothesis. There is sufficient statistical evidence at the a = 0.01 significance level to conclude that there is a significant linear relationship between strain and stress in MD simulations.

H. Yao et al. / Construction and Building Materials 162 (2018) 430–441

Stress (pa)

-1.5e+08

a

5.0e+07

440

0.000

0.001

0.002

0.003

0.004

0.005

Strain

1 -2 -1 0

Stress (pa)

2

b

-1.5

-1.0

-0.5

0.0

0.5

1.0

1.5

Strain Fig. 8. The scatter diagram with the linear regression line for the relationship between stress and strain (data from Fig. 4c in the manuscript).

The control and xGNP modified asphalt binders were simulated in this study using the Molecular Dynamics (MD) method, and different properties of these asphalt binder models were simulated and computed. The multi-layer graphene model was generated to represent the xGNP particles, which were used to modify the control asphalt in the experiments. The xGNP modified asphalt binder model was prepared using the same conditions and procedures from the experiments. After calculating different properties, the following summaries can be drawn.

the increase in temperature. The trends of the bulk modulus in the MD simulations agree well with those of the reference and laboratory data. The difference between the simulation and experimental results was observed due to the time scale and limited molecules. 2) The test process of the Young’s modulus in the asphalt binder models was simulated using the MD method, and the Young’s modulus of the models was calculated based on the second regression method. The Young’s moduli of these models increase linearly with decreasing test temperature. The trend in Young’s modulus in MD simulations is similar to the predicted trend in the experiments. The Young’s moduli of the control and xGNP modified asphalt binder models are larger than laboratory data, and it is caused by several factors: time steps, simulation time, molecular numbers and calculation methods. 3) The shear modulus of the asphalt binder models was calculated by applying a small shear strain to the models in one direction (XY direction) based on the first method mentioned above. The shear moduli of these asphalt models decrease exponentially with increasing test temperature. The trends in the MD simulation data of these models fit the regressed curves of the experimental results. The Young’s moduli of these models are higher than the experimental results. The large modulus of the models is induced by the time scale and limited molecules in MD simulations. 4) Poisson’s ratios of the asphalt binder models are also obtained from different modulus results. The range of Poisson’s ratio is the same as that of the experimental data. It implies that the stress response of asphalt binder models is reasonable after applying the small strain to these models.

1) The density of the xGNP modified asphalt binder model is higher than that of the control asphalt, and it indicates that the addition of xGNP model increases the density of the asphalt binder model. The test process of the bulk modulus in the laboratory was simulated using the MD algorithm, and the bulk moduli of the control and xGNP modified asphalt binder models were calculated with two methods. There is no difference between the results of the two methods. The bulk moduli of the models linearly decrease with

In this paper, the modulus properties of the control and xGNP modified asphalt binder models were analyzed and compared. The trends in different moduli of the asphalt binder models were the same as those of the experimental results, and good predictions of Poisson’s ratio in the asphalt binder models were observed. The contributions of this study are that: 1) obtaining the Young’s modulus of the asphalt binder models using the NPH ensemble; 2) utilizing NVT and SLLOD algorithm for the calculation of bulk and shear moduli of the asphalt binder models with the Lee-Edwards

6.2. Analysis after data standardization and normalization Both stress and stain were standardized and normalized by subtracting their mean and then dividing the value by their standard deviation. Then, the possible atypical observations and outliers were removed from the dataset if their normalized values were either larger than 2.5 or smaller than 2.5. This technique helped reduce the effects of outliers and provide the probable states that contribute more to the average. The scatter diagram with the linear regression line is plotted in Fig. 8. The upper panel (Fig. 8a) is the plot based on the original data, and the lower panel (Fig. 8b) is based on the normalized data after removing possible outliers. The p-value for the model adequacy based on the original data is 0.01393 and the p-value based on the standardized data after removing outliers is 0.02969. Both p-values show that the linear model is adequate for the relationship between stress and strain at a 5% significance level. This analysis further strengthens the linear relationship after removing atypical values to put more weight into the more probable state. 7. Conclusions

H. Yao et al. / Construction and Building Materials 162 (2018) 430–441

boundary conditions; and 3) applying a correlation analysis to reveal the rules for producing data in MD simulation. Even though a large difference between the simulation and experimental data was observed, which is likely due to time frame and complexity as well as the elastic response of the model materials, the research method is useful for the simulation of other materials. In addition, future research by our group would focus on the self-healing properties of the asphalt models and reduce the differences between the simulation and laboratory results. Acknowledgements The authors appreciate the financial support of the U.S. National Science Foundation (NSF) under grant 1300286. The computational studies were performed using Scienomics MAPS software suite (MAPS, Version 3.4, Scienomics, Paris, France, 2014) and a computer cluster (Superior Research Center) at Michigan Technological University. Any opinion, finding, and conclusions expressed in this paper are those of the authors and do not necessarily represent the view of any organization. References [1] H. Yao, Q. Dai, Z. You, Molecular dynamics simulation of physicochemical properties of the asphalt model, Fuel 164 (2016) 83–93. [2] Superpave Series No 1(SP-1), Superpave Performance Graded Asphalt Binder Specification and Testing, Asphalt Institute, Lexington KY, USA, 2003. [3] L. Zhang, M.L. Greenfield, Analyzing properties of model asphalts using molecular simulation, Energy Fuels 21 (3) (2007) 1712–1716. [4] Jennings PWSHRPNRC, Binder characterization and evaluation by nuclear magnetic resonance spectroscopy, Strategic Highway Research Program, National Research Council, SHRP-A-335, Washington, DC, 1993. [5] X. Shu, B. Huang, Recycling of waste tire rubber in asphalt and portland cement concrete: an overview, Constr. Build. Mater. 67 (2014) 217–224. [6] Presti D. Lo, Recycled Tyre Rubber modified bitumens for road asphalt mixtures: a literature review, Constr. Build. Mater. 49 (2013) 863–881. [7] C. Thodesen, K. Shatanawi, S. Amirkhanian, Effect of crumb rubber characteristics on crumb rubber modified (CRM) binder viscosity, Constr. Build. Mater. 23 (1) (2009) 295–303. [8] J. Shen, S. Amirkhanian, F. Xiao, B. Tang, Influence of surface area and size of crumb rubber on high temperature properties of crumb rubber modified binders, Constr. Build. Mater. 23 (1) (2009) 304–310. [9] G.D. Airey, T.M. Singleton, A.C. Collop, Properties of polymer modified bitumen after rubber-bitumen interaction, J. Mater. Civil Eng. 14 (4) (2002) 344–354. [10] S. Kim, S.W. Loh, H. Zhai, H.U. Bahia, Advanced characterization of crumb rubber-modified asphalts, using protocols developed for complex binders, Transp. Res. Rec.: Transp. Res. Board (2001) 15–24. [11] H. Yao, Z. You, L. Li, S.W. Goh, J. Mills-Beale, X. Shi, et al., Evaluation of asphalt blended with low percentage of carbon micro-fiber and nanoclay, J. Test. Eval. 41 (2) (2013) 278–288. [12] M.J. Khattak, A. Khattab, H.R. Rizvi, in: Mechanistic Characteristics of Asphalt Binder and Asphalt Matrix Modified with Nano-Fibers, ASCE, Dallas, TX, 2011, p. 492. [13] Z. Ge, H. Wang, Q. Zhang, C. Xiong, Glass fiber reinforced asphalt membrane for interlayer bonding between asphalt overlay and concrete pavement, Constr. Build. Mater. 101 (Part 1) (2015) 918–925. [14] H. Yao, Z. You, L. Li, C.H. Lee, D. Wingard, Y.K. Yap, et al., Rheological properties and chemical bonding of asphalt modified with nanosilica, J. Mater. Civil Eng. 25 (11) (2013) 1619–1630. [15] H. Yao, Z. You, L. Li, X. Shi, S.W. Goh, J. Mills-Beale, et al., Performance of asphalt binder blended with non-modified and polymer-modified nanoclay, Constr. Build. Mater. 35 (2012) 159–170. [16] H. Yao, Z. You, L. Li, S.W. Goh, C.H. Lee, Y.K. Yap, et al., Rheological properties and chemical analysis of nanoclay and carbon microfiber modified asphalt with Fourier transform infrared spectroscopy, Constr. Build. Mater. 38 (2013) 327–337. [17] H. Wang, Z. You, J. Mills-Beale, P. Hao, Laboratory evaluation on high temperature viscosity and low temperature stiffness of asphalt binder with high percent scrap tire rubber, Constr. Build. Mater. 26 (1) (2012) 583–590. [18] M. Khattak, A. Khattab, P. Zhang, H. Rizvi, T. Pesacreta, Microstructure and fracture morphology of carbon nano-fiber modified asphalt and hot mix asphalt mixtures, Mater. Struct. 46 (12) (2013) 2045–2057.

441

[19] I. Al-Qadi, H. Wang, Impact of wide-base tires on pavements, Transp. Res. Rec.: J. Transp. Res. Board 2304 (2012) 169–176. [20] H. Li, J. Harvey, D. Jones, Multi-dimensional transient temperature simulation and back-calculation for thermal properties of building materials, Build. Environ. 59 (2013) 501–516. [21] M.R. Mohd Hasan, B. Colbert, Z. You, A. Jamshidi, P.A. Heiden, M.O. Hamzah, A simple treatment of electronic-waste plastics to produce asphalt binder additives with improved properties, Constr. Build. Mater. 110 (2016) 79–88. [22] M. Karplus, J.A. McCammon, Molecular dynamics simulations of biomolecules, Nat. Struct. Mol. Biol. 9 (9) (2002) 646–652. [23] H. Yao, Q. Dai, Z. You, A. Bick, M. Wang, S. Guo, Property analysis of exfoliated graphite nanoplatelets modified asphalt model using molecular dynamics (MD) method, Appl. Sci. 7 (1) (2017) 43. [24] L. Zhang, M.L. Greenfield, Effects of polymer modification on properties and microstructure of model asphalt systems, Energy Fuels 22 (5) (2008) 3363– 3375. [25] D.D. Li, M.L. Greenfield, Chemical compositions of improved model asphalt systems for molecular simulations, Fuel 115 (2014) 347–356. [26] W.D. Cornell, P. Cieplak, C.I. Bayly, I.R. Gould, K.M. Merz, D.M. Ferguson, et al., A second generation force field for the simulation of proteins, nucleic acids, and organic molecules, J. Am. Chem. Soc. 117 (19) (1995) 5179–5197. [27] J. Wang, R.M. Wolf, J.W. Caldwell, P.A. Kollman, D.A. Case, Development and testing of a general amber force field, J. Comput. Chem. 25 (9) (2004) 1157– 1174. [28] M. Valiev, E.J. Bylaska, N. Govind, K. Kowalski, T.P. Straatsma, H.J.J. Van Dam, et al., NWChem: a comprehensive and scalable open-source solution for large scale molecular simulations, Comput. Phys. Commun. 181 (9) (2010) 1477– 1489. [29] R.W. Hockney, J.W. Eastwood, Computer Simulation Using Particles, CRC Press, 1988. [30] W.H. Press, Numerical Recipes in C: The Art of Scientific Computing, Cambridge University Press, 1992. [31] M. Hazewinkel, Encyclopaedia of Mathematics (set), Springer, Netherlands, 1994. [32] L. Artok, Y. Su, Y. Hirose, M. Hosokawa, S. Murata, M. Nomura, Structure and reactivity of petroleum-derived asphalteney, Energy Fuels 13 (2) (1999) 287– 296. [33] I. Kowalewski, M. Vandenbroucke, A.Y. Huc, M.J. Taylor, J.L. Faulon, Preliminary results on molecular modeling of asphaltenes using structure elucidation programs in conjunction with molecular simulation programs, Energy Fuels 10 (1) (1996) 97–107. [34] H. Groenzin, O.C. Mullins, Molecular size and structure of asphaltenes from various sources, Energy Fuels 14 (3) (2000) 677–684. [35] D. Sun, T. Lin, X. Zhu, Y. Tian, F. Liu, Indices for self-healing performance assessments based on molecular dynamics simulation of asphalt binders, Comput. Mater. Sci. 114 (2016) 86–93. [36] Y. Ding, B. Huang, X. Shu, Y. Zhang, M.E. Woods, Use of molecular dynamics to investigate diffusion between virgin and aged asphalt binders, Fuel 174 (2016) 267–273. [37] F. Khabaz, R. Khare, Glass transition and molecular mobility in styrenebutadiene rubber modified asphalt, J. Phys. Chem. B. 119 (44) (2015) 14261– 14269. [38] B.D. Jensen, K.E. Wise, G.M. Odegard, The effect of time step, thermostat, and strain rate on ReaxFF simulations of mechanical failure in diamond, graphene, and carbon nanotube, J. Comput. Chem. 36 (21) (2015) 1587–1596. [39] H. Wang, E. Lin, G. Xu, Molecular dynamics simulation of asphalt-aggregate interface adhesion strength with moisture effect, Int. J. Pavement Eng. 1–10 (2015). [40] A. Bandyopadhyay, Molecular Modeling of EPON 862-DETDA Polymer, Michigan Technological University, 2012. [41] S. Nosé, A unified formulation of the constant temperature molecular dynamics methods, J. Chem. Phys. 81 (1) (1984) 511–519. [42] W.G. Hoover, Canonical dynamics: equilibrium phase-space distributions, Phys. Rev. A. 31 (3) (1985) 1695–1697. [43] D.J. Evans, G.P. Morriss, Nonlinear-response theory for steady planar Couette flow, Phys. Rev. A. 30 (3) (1984) 1528–1530. [44] A. Motamed, A. Bhasin, K. Liechti, Using the poker-chip test for determining the bulk modulus of asphalt binders, Mech. Time-Depend Mater. 18 (1) (2014) 197–215. [45] A.D. McNaught, A. Wilkinson, IUPAC Compendium of Chemical Terminology (Gold Book), Second ed., Oxford Blackwell Scientific Publications, 1997. [46] R.A. Tarefder, A.M. Zaman, W. Uddin, Determining hardness and elastic modulus of asphalt by nanoindentation, Int. J. Geome. 10 (3) (2010) 106–116. [47] S.S. Zumdahl, S.A. Zumdahl, Chemistry, Cengage Learning, 2008. [48] H.I. Ling, A. Smyth, R. Betti, Poromechanics IV, DEStech Publications Incorporated, 2009. [49] J.T. McClave, T.T. Sincich, Statistics, Pearson Education, 2012.