Property Analysis of Exfoliated Graphite Nanoplatelets Modified Asphalt Model Using Molecular Dynamics (MD) Method

Abstract

This Molecular Dynamics (MD) simulation paper presents a physical property comparison study between exfoliated graphite nanoplatelets (xGNP) modified and control asphalt models, including density, glass transition temperature, viscosity and thermal conductivity. The three-component control asphalt model consists of asphaltenes, aromatics, and saturates based on previous references. The xGNP asphalt model was built by incorporating an xGNP and control asphalt model and controlling mass ratios to represent the laboratory prepared samples. The Amber Cornell Extension Force Field (ACEFF) was used with assigned molecular electro-static potential (ESP) charge from NWChem analysis. After optimization and ensemble relaxation, the properties of the control and xGNP modified asphalt models were computed and analyzed using the MD method. The MD simulated results have a similar trend as the test results. The property analysis showed that: (1) the density of the xGNP modified model is higher than that of the control model; (2) the glass transition temperature of the xGNP modified model is closer to the laboratory data of the Strategic Highway Research Program (SHRP) asphalt binders than that of the control model; (3) the viscosities of the xGNP modified model at different temperatures are higher than those of the control model, and it coincides with the trend in the laboratory data; (4) the thermal conductivities of the xGNP modified asphalt model are higher than those of the control asphalt model at different temperatures, and it is consistent with the trend in the laboratory data.

1. Introduction

1.1. Asphalt Material

Asphalt is a byproduct of petroleum refinement and is also widely applied to many fields such as transportation, recreation, building construction, etc. Around 90% of asphalt consists of carbon and hydrogen. Based on the Corbett method, asphalt can be separated into four components: asphaltenes, saturates, napthene aromatics, and polar aromatics. Asphalt is composed of asphaltenes, paraffins, first acidiffins, second acidiffins, and nitrogen bases using the Rostler method [1]. These different types of molecules in asphalt interact with each other and affect the chemo-physical properties of asphalt [1].

Due to their special properties, nanomaterials have been introduced and used in different fields to enhance composite materials [2]. Some of their special properties include the dominance of interfacial phenomena, size and quantum effects, etc. [3]. Nanomaterials are used for electronics, agriculture, construction, food and medical technologies. Also, different types of nanoclay have been widely used in the modification of asphalt. The test results show that the layered structure of nanoclay improved the high-temperature performance of asphalt and the resistance to rutting and fatigue cracking [4,5]. The nanosilica material was also used and added to the asphalt matrix to improve performance. The micro images of nanosilica modified asphalt were observed, and the test results indicate that the resistance to permanent deformation in the modified asphalt improved [6]. The literature shows that graphite was used for the modification of asphalt, and the addition of graphite improved the electrical property of asphalt [7,8]. In this study, the multi-layer graphite sheets were applied to modify the asphalt in consideration of the high thermal stability, self-lubrication, and high electrical conductivity of multi-layer graphite sheets [9,10,11,12]. This is also the motivation to use the material for modification of the asphalt model.

1.2. Molecular Dynamics (MD) Method

Molecular Dynamics (MD), originating in theoretical physics, was applied widely in materials science [13,14]. MD is a kind of computer program to simulate the movements of atoms in materials, and the atoms and molecules interact for a designated time based on the Newton’s law of motion. The trajectories of atoms and molecules are monitored, and the energy of the system is computed. In physics, MD was used to examine physical properties [13,14]. The evolution of dynamics in a single molecule is used to determine the macro properties of the system. The “statistical mechanics by number” and “Laplace’s vision of Newtonian mechanics” were also used to describe the molecular dynamics. The simulation size and total duration were selected so that the calculation can be finished within a reasonable amount of time [15]. The Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) [16] and the Monte Carlo for Complex Chemical Systems (MCCCS) program [17] were commonly used for MD simulations. The computation algorithm of the MD simulation is shown in Figure 1. In addition, compared to other methods (such as Finite Element Method (FEM) and Discrete Element Method (DEM)), the MD method helps address specific problems or principles of atoms or materials, and also, a specific property can be studied by altering specific contributions. Moreover, the material behaviors or response can be analyzed under extreme conditions on a nanoscale.

Recently, three components of mixtures (asphaltenes, aromatics and saturates) were chosen to represent the asphalt using MD simulation. In this reference model, 1,7-dimethylnaphthalene (C₁₂H₁₂) and docosane (n-C₂₂H₄₆) were represented as naphthene aromatics and saturates, respectively [18]. Two kinds of asphaltene structures were used, and the density of each component was calculated using the MD method. From the simulation results, there are still many differences between the test and the simulation data in the calculations of the glass transition temperature and bulk modulus [18]. A new asphalt model with four components (asphaltenes, polar aromatics, naphthene aromatics and saturates) was created, and the density and thermal expansion coefficient of the asphalt model were calculated. The MD simulation results agreed with the laboratory data [19]. In addition, polystyrene was added to this asphalt model for polymer modification analysis. The radial distribution functions g(r) of components of the asphalt model were computed [19]. The MD simulation was also recently applied to study the self-diffusivity properties of asphalt binders. The effect of healing on the fatigue performance of binders was studied, and the MD model of an asphalt binder was created to predict the healing effect. The self-diffusion effect caused the binder molecules to flow across the crack interface. The correlation between the length and branching of molecules and self-diffusion of asphalt molecules was investigated [20]. The relationship between the asphalt and aggregate was established and the asphalt-quartz structure model of the interface was used in the system. The Consistent-Valence Force Field (CVFF_aug) was considered in this simulation to characterize the inter-atom interactions [21]. Based on a literature review, few researchers used MD to simulate modified asphalt and studied the physical properties of modified asphalt models in engineering disciplines. Graphite has a high thermal stability, high electrical conductivity, good self-lubricating and dry lubricating properties. The graphite was used to modify asphalt binder, and thermal conductivity and anti-aging properties improved after the addition of graphite in the asphalt binder [22,23]. Due to these improvements, in this study, the common multi-layer graphite model was adopted, and the components of the control asphalt model was composed of asphaltenes, saturates and aromatics. The Amber Cornell Extension Force Field (ACEFF) was developed and used to simulate the asphalt modified with exfoliated multi-layered graphite nanoplatelets (xGNP).

1.3. Objectives and Sections

The objectives of this study are to use the MD method to simulate and compare the properties of the xGNP modified and control asphalt model. The three-component control asphalt model consists of asphaltenes, aromatics, and saturates based on previous work [24]. The MD simulation and optimization methods were described in Section 2, as well as the force field. The common multi-layer graphite model was incorporated in the xGNP modified asphalt model to represent the xGNP modifier. The Amber Cornell Extension Force Field and Electrostatic Potential (ESP) charges were assigned to the components of the control and modified asphalt models, as described in Section 3. Different properties of the asphalt models were computed including the density, glass transition temperature (Tg), viscosity and thermal conductivity. The MD simulation data of the control and modified models and their laboratory results were compared in Section 4, Section 5 and Section 6.

2. Force Field and Optimization Methods

2.1. Classical MD Simulation Methods

Different ensembles can be used in the MD method, including the Microcanonical ensemble (NVE ensemble), Canonical ensemble (NVT ensemble), Isothermal-isobaric ensemble (NPT ensemble), Isoenthalpic-isobaric ensemble (NPH ensemble), and Generalized ensembles [25]. For instance, in the NVE ensemble, the number of moles (N), volume (V) and energy (E) in the isolated system are not be changed. The system experiences the adiabatic process and no heat exchange would occur.

2.2. Force Field

A force field, presented as parameters of mathematical functions in molecular mechanics, was used to describe the energies of the atoms. Force field parameters and functions are obtained from experimental tests and quantum mechanical calculations. Many force fields were developed and introduced by researchers, such as the Chemistry at HARvard Molecular Mechanics (CHARMM) Force Field [26], Assisted Model Building and Energy Refinement (AMBER) Force Field [27], Condensed-phase Optimized Molecular Potential for Atomistic Simulation Studies (COMPASS) Force Field [28], Optimized Potential for Liquid Simulation (OPLS) Force Field [22] and DREIDING Force Field [29]. In this study, the Amber Cornell Extension Force Field (ACEFF) was used to define the movement in the molecular system based on the Amber Cornell Force Field [30], and the experimental parameters in this force field were referenced from the General Amber Force Field (GAFF) [31]. The formula is shown in Equation (1).

$$E_{total}={{\displaystyle \sum }}_{bonds}{K}_r{\left(r-r_{eq}\right)}^2+{{\displaystyle \sum }}_{angles}{K}_{\theta }{\left(\theta -{\theta }_{eq}\right)}^2+{{\displaystyle \sum }}_{dihedrals}\frac{V_n}{2}\left[1+\cos \left(n\phi -\gamma \right)\right]+{{\displaystyle \sum }}_{i<j}[\frac{A_{ij}}{R_{ij}^{12}}-\frac{B_{ij}}{R_{ij}^6}+\frac{q_i{q}_j}{ϵR_{ij}}]$$

(1)

where $r_{eq}$ and ${\theta }_{eq}$ are the equilibrium structural data from an X-ray test; $K_r$ is the force coefficient determined by linear interpolation between the single and double bond values; $K_{\theta }$ is the force coefficient from vibrational analysis of a simple sp2 atom; n is the multiplicity for dihedrals; $\mathsf{\gamma }$ is the phase angle for the torsional angle parameters; A, B and q are the non-bonded potentials between atom pairs; $R_{ij}$ is the distance between the atoms; and $ϵ$ is the well depth for van der Waals energy.

2.3. Optimization Methods

When the molecular systems are built, energy optimization and data smoothing are required to optimize the molecular systems and output results, respectively. These procedures help the researchers understand more about the systems. The following methods were used in this study.

(1) Conjugate Gradient Method

The conjugate gradient method is a kind of iterative algorithm to solve the partial differential equations. The solutions for unconstrained optimization problems like energy minimizations were also developed by Hestenes and Stiefel [32]. The formula is shown in Equation (2). The iterative method was essential and required for energy optimization (lowest energy) of large system, and an initial guess at the solution was the start of the iteration approach. The iterative method does not provide the exact solution, but can improve the approximation after a certain number of iterations. This function was also restricted by computational resources.

$$f\left(x\right)=\frac{1}{2}{x}^{T}Ax-x^{T}b,x\in R^n(\text{iterative method})$$

(2)

where A is symmetric, positive and real; b is a known coefficient; and vectors n and T are non-zero.

(2) Particle-Particle-Particle-Mesh Method (PPPM)

The Particle-Particle-Particle-Mesh Method (PPPM or P³M) is used to compute long-rang electrostatic force, which can be divided into two parts: short- and long-range interparticle forces [33]. Short-range interactions are computed from the particle-particle (PP) calculation, and the long-range interactions are processed by the particle-mesh method [34]. The formula is shown in Equation (3). During the energy or force calculation in the system, the particles normally are forced to occupy a low spatial resolution, and it may cause errors in the results. The P³M method is designed to calculate potential through a direct sum for particles. The P³M method was adopted to consider the speed and accuracy of simulations in this study.

$$F_{ij}=F_{ij}^{sr}+F_{ij}^m$$

(3)

where $F_{ij}$ is the interparticle forces in the system; $F_{ij}^{sr}$ is the rapidly varying short-range component; $F_{ij}^m$ is the slowly varying component.

3. Model Generation

Based on previous work of the authors [24] and reference [18], three components were used in the control asphalt including the asphaltenes, aromatics, and saturates at a ratio of 5:27:41. The ratio of asphaltenes, aromatics and saturates was cited form the reference [18], and this ratio (asphaltenes, aromatics and saturates at 5:27:41) in the asphalt model was based on the asphaltene mass fraction (20.7 wt%) and alkane/aromatic carbon ratio (59.6 wt%, 19.7 wt%; wt% is weight percentage) [35]. The 1,7-Dimethylnaphthalene [36] (Figure 2a) and docosane [37] (Figure 2b) were used to represent the aromatics and saturates, respectively. 1,7-Dimethylnaphthalene (C₁₂H₁₂) was reported by Zhang [18] and Groenzin [36] based on the ratio of alkane and aromatic in the asphalt. The docosane (C₂₂H₄₆) was reported by Zhang [18] and Kowalwski [37] based on the properties (melting and boiling points) of docosane and saturates. The asphaltene structure (Figure 2c) is from the references [18,38], C₆₄H₅₂S₂. In addition, the ESP was calculated and assigned to the components by NWChem. The control asphalt model was built through the compression of three components with NPT ensemble, and running time is around 1 ns. The start density of the control asphalt system is around 0.1 g/cm³, and the target density is 1.0 g/cm³, which is similar to that of the real asphalt. The low start density is good for relocations of atoms or components of systems during the energy optimization. Based on the previous property calculation of the asphalt model with the same components [18,24], the properties of the asphalt model are similar to those of the real asphalt tested in the laboratory. In this study, the ACEFF and ESP were assigned in the MD systems, and it is expected that the improvement of properties will be observed during the calculations. Therefore, the control asphalt model with ACEFF was generated to represent the control asphalt PG 58-28.

In the laboratory tests, the modifier, xGNP graphene nanoplatelets, used in this study is produced by XG Sciences Inc., and its micro-image (Figure 2d) was obtained by the field emission scanning electron microscope (FE-SEM). The distance between the graphene layers is around 3.35 Å, and the mole mass content of xGNP nanoplatelets in the modified asphalt is 2% by the weight of the control asphalt. During the preparation of xGNP modified asphalt in the laboratory, 2% xGNP multi-layer graphite particles were slowly added in the asphalt matrix at the temperature of 145 ºC. The modified asphalt was sheared in the high shear machine for two hours to ensure that particles were well dispersed. Similarly, in the simulation test, the common multi-layer graphite model (Figure 2e) was used to represent the xGNP nanoplatelets, and 2% xGNP nanoplatelets by the weight of control model were randomly added to the control model. Mass mentioned in this study is based on 1 mole of the simulation box. The xGNP model with four layers was placed in the control asphalt model, and NPT ensemble was employed to compress the modified system. The xGNP modified asphalt model was generated and is shown in Figure 2f. The composition of the modified asphalt system is shown in

Table 1. The composition of xGNP modified asphalt system.

Modified Asphalt Model	Mass (g/mol)	Sum Formula	Number of Atoms per Molecule	Number of Bonds per Molecule	Number of Molecules	Total Mass (g)	Mass Fraction (%)
Asphaltene	885.23	C64H52S2	118	132	20	17,704.59	20.25
1,7-Dimethylnaphthalene	156.22	C12H12	24	25	108	16,872.16	19.30
Docosane	310.6	C22H46	68	67	164	50,938.50	58.30
xGNP	472.53	C38H16	54	65	4	1890.13	2.15
Modified asphalt model	-	-	-	-	296	87,405.39	-

. Different optimization methods mentioned above, the conjugated gradient method and the PPPM method, were adopted during the energy optimizations. The optimized system with the lowest energy was stable for calculations.

4. Physical Properties of the Control and xGNP Modified Asphalt Models

4.1. Density

When the control asphalt and modified asphalt models were generated using molecular dynamics, the densities of these models were computed at the conditions of room temperature and standard atmosphere pressure to evaluate the molecular model. The LAMMPS and optimization methods mentioned above were used to conduct the experimental MD simulations. The NPT ensemble simulations were employed to compress or relax the unit cell. The temperatures, pressures, and densities of the control and xGNP modified asphalt models are shown in Figure 3.

Figure 3a,b show the temperatures and pressures in the control and modified asphalt systems under the NPT ensemble. The systems were run more than 1 ns to be stable and optimized. Some of the results with the simulation steps are shown in Figure 3. The data were fitted by a Savitzky-Golay filter [39] (a kind of generalized moving average) with a span (a parameter to control the average) of 10%. The temperatures of the control and xGNP modified asphalt are close to 298 K during these steps, and the pressures fluctuate around 1 atmosphere (atm). Meanwhile, it is interesting to note that the data fluctuation of the control model is obviously greater than that of the modified asphalt model due to the difference in the molecular number. The fluctuated data range varies with a $1/\sqrt{N}$ variation based on the baseline of the moving average, where N is the number of molecules in the MD system [40]. Figure 3c shows densities of the control and xGNP modified asphalt systems, and the data was fitted by a Savitzky-Golay filter [39] with a span of 10%. The densities of the xGNP modified asphalt system are larger than those of the control asphalt model. It is reasonable that the addition of xGNP particles in the control system increases the density of the modified system. It is apparent that the density data amplitude of the control system is larger than that of the xGNP modified asphalt model, as well as the stability, due to their being more molecules. In addition, the densities of the control and xGNP modified asphalt systems are also close to the laboratory (0.95 g/cm³–1.05 g/cm³) [18,40,41] and reference data [18]. Therefore, the modified asphalt model can be deemed reasonable by obtaining a comparable mass density. Figure 3d displays the density curve of the control and xGNP modified asphalt with different temperatures, which range from 233.15 K to 443.15 K. The densities of the asphalt models decrease with the increase in temperature, and the density of xGNP modified system is slightly greater than that of the control model under different temperatures. Meanwhile, the density trends of the control and xGNP modified models are similar to that of the reference model [18]. In addition, the data fluctuation amplitude of xGNP modified model is smaller than that of the control system.

4.2. Glass Transition Temperature

The glass transition temperatures (Tg) of materials are influenced by their components, and the addition of a new component in the material results in a difference in the glass transition temperature. The formula [42] of the glass transition temperatures of composite materials is shown in Equation (4). The asphalt transfers from the viscoelastic state to a brittle one. The internal stress increases and thermal energy are insufficient below the glass transition temperature in the materials. Therefore, the glass transition temperature is an important parameter or property for amorphous materials, and it should be low for a good low-temperature performance [42]. The glass transition temperature is the temperature where two asymptotes intersect on the specific volume-temperature curve. In the laboratory, the glass transition temperature of materials can be tested by differential scanning calorimetry (DSC). In this MD simulation study, the control and xGNP modified asphalt systems were relaxed under the NPT ensemble with a temperature range of 233.15 K–443.15 K. The specific volumes of these systems were calculated and the glass transition temperatures of the models were determined. The MD simulation results of the control and xGNP modified asphalt models are shown in Figure 4.

$$\frac{1}{T_g}={\displaystyle \sum }\frac{w_i}{T_{gi}}$$

(4)

where T_g and T_gi are the glass transition temperature of the composite material and its component, respectively; and w_i is the mass fraction of the component.

Figure 4a or Figure 4b shows the specific volumes of the control system or xGNP modified asphalt model, respectively. Figure 4c shows the specific volumes of both control and xGNP modified asphalt models. The specific volumes increase with the increasing temperatures of the models. The glass transition temperature of the control asphalt model is around 300 K [24]. As shown in the Figure, it is deduced that the Tg of the xGNP modified asphalt system is around 250 K. Based on the laboratory data in the reference [43], the Tg of asphalt ranges from 223 K to 303 K. The Tg of the modified asphalt model is within the laboratory data range, and it is also better than the reference data [18] (298 K–358 K) of the asphalt model. In order to get better thermal relaxation in the asphalt, a low Tg of asphalt is expected. Figure 4d shows the comparisons of glass transition temperatures of the references and MD simulations. The glass transition temperatures of the Strategic Highway Research Program (SHRP) asphalt binders are all around 250 K [42], including SHRP asphalt AAA-1, AAB-1, AAC-1, AAD-1, AAF-1, AAG-1, AAK-1 and AAM-1. The Tg of xGNP modified asphalt model is close to the glass transition temperatures of SHRP asphalt binders. Therefore, the Tg of xGNP modified asphalt is reasonable and better than the results of the control [24] and reference asphalt models [18].

5. Rheological Properties of the Control and xGNP Modified Asphalt Models

5.1. Viscosity Measurement Method and Results

Dynamic shear viscosity is an important parameter of fluids to measure the resistance to gradual deformation induced by shear stress. If the shear speed caused by the external force is appropriate, the fluid particles move parallel to the particles sheared. The speed varies linearly from the sheared layers to different layers. The resistance between these layers is caused by friction. The formula to calculate viscosity is shown in Equation (5).

$$\eta =\frac{Fy}{Au}$$

(5)

where $F$ is the applied force; A is the area of the plate; $\eta$ is the dynamic shear viscosity; and u/y is the shear gradient.

During the construction of asphalt pavement, the viscosity determines the mixing and compaction temperatures, which relates to the pump ability, mix ability and workability of asphalt. Based on the American Society for Testing and Materials (ASTM) D4402, the Brookfield DV-II plus viscometer (Figure 5a) was selected to test the viscosity of asphalt in the laboratory. The viscosity test results are shown in Figure 5b.

From the test data, the viscosities of the xGNP modified asphalt binder (light blue line, Figure 5b) are higher than those of the control asphalt binder (red line). It indicates that multi-layer xGNP particles increase the viscosity of the modified asphalt binder. In addition, the viscosity decreases with the increase in temperatures of the asphalt binder. Exponential trends were also observed in the test data. The MD simulation for viscosity of the asphalt binder model was discussed in the following Section 5.2.

5.2. MD Viscosity Aimulation Methods and Results

In the MD experimental simulation, there are four common methods to calculate the dynamic shear viscosity in the MD systems [16]: (1) a non-equilibrium MD (NEMD): the unit cell is sheared by “fix deform” and the temperature is controlled; (2) a NEMD: the viscosity is computed through the velocity and pressure in the systems; (3) a reverse non-equilibrium MD (rNEMD): the momentum flux is transferred between different layers in the unit cell through the Muller-Plathe algorithm; (4) equilibrium MD (EMD): the Green-Kubo (GK) formula is used to compute the viscosity, and continuous momentum flows are applied in the unit cell.

In this MD study, the Muller-Plathe method was used to calculate the viscosities of the control and xGNP modified asphalt systems. The unit cell of asphalt models was split into 20 layers. The viscosity calculation is shown in Equation (6). During the calculation of viscosity in the MD simulation, unit conversion is needed. The viscosity unit in the MD simulation is gram/mol/angstrom/femtosecond, but the unit in the laboratory test is kilogram/meter/second. Avogadro’s constant is used to convert from microscopic to macroscopic states. In addition, the momentums transferred in the control and xGNP modified systems at the temperature of 443 K are shown in Figure 6a. The viscosities of the control and xGNP modified asphalt systems at a temperature of 423 K are shown in Figure 6b,c. The MD and laboratory viscosity results of the control and xGNP modified asphalt systems at different temperatures are shown in Section 5.4.

$$j_z\left(p_x\right)=-\eta \frac{\partial {\nu }_x}{{\partial }_z}\text{ and }j\left(p_x\right)=\frac{P_x}{2tA}$$

(6)

where $\eta$ is the dynamic shear viscosity; $\frac{\partial {\nu }_x}{{\partial }_z}$ is the shear rate; $j\left(p_x\right)$ is the input momentum flux; $P_x$ is the input momentum; $t$ is the simulation time; $A=L_x{L}_y$, $L_x$ is the length of unit cell in the x direction; and $L_y$ is the length of unit cell in the y direction.

Figure 6a displays the momentum transferred in the control and xGNP modified asphalt binder systems at a temperature of 443 K, as well as some of the results from different temperatures for repeatability. The momentum transferred in the xGNP modified asphalt binder model is more than that of the control model. The temperatures in the control asphalt binder model fluctuate more than those of the xGNP modified asphalt binder model due to there being fewer molecules. The temperatures are also around 443 K, and do not have large variations. It indicates that more molecules in the MD system lead to less data variation and a stable structure of materials. It coincides with the conclusion that the vibrating range of MD data is within a $1/\sqrt{N}$ variation [40] (N is the number of molecules in the system). Figure 6b,c demonstrates the viscosities of the control and xGNP modified asphalt binder systems at a temperature of 423 K, and it is also part of data analysis under different temperatures. The temperature fluctuation of the control model is larger than that of xGNP modified asphalt binder model due to fewer molecules compared to xGNP modified asphalt binder system. The temperatures are also centered at 423 K with a 20 K variation. It coincides with the temperature setting of simulations. Furthermore, the statistical analysis of viscosities of the control and xGNP modified asphalt models was performed to analyze the distributions of viscosity data in MD simulation.

5.3. Statistical Analysis for Viscosities of the Control and xGNP Modified Asphalt Models

The variation of data in the viscosity calculation is observed in the last section. The statistical analysis was used to better understand the data distribution, and it is also good for describing and reproducing the data. It is well known that many experimental or observational data arising in engineering is shaped by a lognormal distribution due to its various appealing properties. If a random variable x follows a lognormal distribution, the random variable Y = log(X) is distributed as a normal. The probability density function (PDF) of a lognormal distribution with parameters μ and σ is given by

$$f\left(x\right)=\frac{1}{x\sigma \sqrt{2\pi }}{e}^{-\frac{{\left(logx-\mu \right)}^2}{2{\sigma }^2}}$$

(7)

where $x>0,-\infty <\mu <\infty ,$ and $\sigma >0$. In this study, we consider the lognormal distribution, because it provides heavier tails compared with the normal one and is thus more flexible to experimental data when studying robustness to outliers. Due to large variations of data from 7.375894 to 53,239.325870, we consider the more appropriate fits based on data ranging from 0 to 300, and from 0 to 500. The parameter estimates of μ and σ with their standard errors in parenthesis are presented in

Table 2. Parameter estimates of μ and σ with their standard errors.

Type	$\mathit{\mu }$	$\mathit{\sigma }$
xGNP (<500)	3.81542596 (0.01393144)	0.85901818 (0.00985102)
Control (<500)	3.01951678 (0.01529361)	0.95801986 (0.01081422)

Note: xGNP: viscosity data of the xGNP modified asphalt model; Control: viscosity data of the control modified asphalt model.

. The histograms with the fitted lognormal distributions are depicted in the top two figures (Figure 7). It can be seen from the two figures that with different choices of truncations, the lognormal distribution provides more flexible fits to xGNP data. A similar conclusion can also been drawn for the fitness of the lognormal distribution to control data (Figure 7). Consequently, we may conclude that the data departure from the lognormality is acceptable or slight, the lognormal distribution is a more robust and flexible model, allowing a better fit as shown above.

5.4. Comparison of Viscosity Predictions of the Control and xGNP Modified Asphalt Models

Based on the data analysis from the MD simulation and laboratory data, the viscosity data between the control and xGNP modified asphalt binder models was compared and analyzed, as well as the experimental results. The exponential regressions were used to fit the viscosity data of the control and xGNP modified asphalt binder models. The comparisons between the control and xGNP modified asphalt binder models were conducted, as shown in Figure 8, including the reference data.

Figure 8 shows the viscosities of the control and xGNP modified asphalt binder models. The viscosities of MD simulations were averaged from the calculations under each separate temperature. The MD simulation viscosities of the xGNP modified asphalt binder model (red line in Figure 8) are also higher than those of the control asphalt binder model. It is similar to the trend in the laboratory data. The exponential trends are observed to be fitted for the MD simulation results due to the same trend in the laboratory data [44]. In this figure, it is obvious that the viscosities of the control and xGNP modified asphalt binder models are higher than those of the reference models [19] (0.65 cp and 1.35 cp at 443.15 K) using Green-Kubo and Einstein (Ein) EMD methods. The relatively flat line is observed in the viscosity data of the control asphalt binder model, and the viscosity results (92.47 cp and 80.13 cp) at 423.15 K and 443.15 K are close to the laboratory data (155.0 cp at 423.15 K and 95.0 cp at 438.15 K in Figure 5b) at two different temperatures for the control asphalt binder model, respectively. It is noticed that there is some improvement between viscosities of the control asphalt binder model (Figure 8) and the asphalt model in the reference [24]. Furthermore, for the xGNP modified asphalt model, the viscosities (452.42 cp, 195.98 cp and 108.96 cp) at 403.15 K, 423.15 K, and 443.15 K, respectively, are very close to the laboratory data (530.0 cp, 270.0 cp and 122.5 cp). The viscosity simulation results of the xGNP modified asphalt binder system are better than those of the control asphalt binder model, and the trend in viscosity of the xGNP modified asphalt binder model is similar to that of the laboratory data. With regard to this improvement, it is caused by the increase in the molecular number in the xGNP modified asphalt model compared to the control asphalt model. It is also an improvement to use the Muller-Plathe method to calculate the viscosity of the asphalt binder model, as well as the optimization methods used in the simulations. It is expected that more molecules in the MD asphalt system improve the accuracy of data prediction. Therefore, the viscosity calculation of MD simulations in the xGNP modified asphalt system with the Amber Cornell Extension Force Field provides a better prediction using the Muller-Plathe method compared to the results of the references [19] and the control asphalt model.

6. Thermal Property of the Control and xGNP Modified Asphalt Models

6.1. Thermal Conductivity Measurement Methods and Results

The thermal conductivity is a kind of measure of materials to transmit the heat energy in a diffusive manner based on Fourier’s law. The materials with a high thermal conductivity are applied to the heat sink, and the materials with a low thermal conductivity are manufactured for thermal insulation. In the laboratory, the xGNP modified asphalt was mixed with ultrasonic stirring during the process of high shear so that the xGNP particles can be homogenously dispersed in the asphalt matrix. The thermal properties’ analyzer (KD2 Pro) was employed to measure the thermal conductivity of asphalt based on the transient line heat source method [45]. The asphalt was placed in the glass tube as shown in Figure 9a,b, and the single needle TR-1, with a 2.4 mm diameter and 60 mm length, was used to test the thermal conductivity. During the heating and cooling processes, the temperature-time relationship is monitored by the sensor located in the needle. The thermal conductivity is calculated with the parameters from the fitted curve for temperature-time. The formula for thermal conductivity is shown in Equation (8). The thermal conductivities of asphalt under different temperatures in the laboratory are shown in Figure 9c.

$$\begin{gathered} T=m_0+m_{2}t+m_{3}lnt(\text{Heating process}) \\ T=m_1+m_{2}t+m_{3}ln\frac{t}{\left(t-t_h\right)}(\text{Cooling process})\text{ and }k=\frac{q}{4\pi m_3} \end{gathered}$$

(8)

where $m_0$ and $m_1$ are the ambient temperatures in the heating and cooling processes, respectively; $m_2$ is the rate of drift of the background temperature; $m_3$ is the slope of a line relating temperature rise to the logarithm of temperature; q is the heat input and $k$ is the thermal conductivity.

Figure 9a,b display the KD2 Pro thermal conductivity apparatus and the temperature control chamber used in the laboratory, respectively. Figure 9c demonstrates the thermal conductivity results tested in the laboratory. The thermal conductivity of the control asphalt binder at room temperature is 0.148 W/m·K, and it is close to the reference data of asphalt [22] (0.170 W/m·K). The thermal conductivity of the xGNP modified asphalt binder at room temperature is 0.226 W/m·K, and it is near to the reference data [22] from 0.396 W/m·K to 0.934 W/m·K with different amounts of graphite (different types from the multi-layer graphite used in this study). Through the test results, the addition of xGNP particle increases the thermal conductivity of the asphalt binder. It is likely that high thermal conductivity of xGNP particles (around 3000 W/m·K) enhances the thermal transfer in the asphalt binder matrix based on the thermal conductivity data of xGNP particles from xgsciences.com [46]. It can also be expected that the light absorption is improved after the addition of xGNP particles in the asphalt binder.

6.2. MD Simulation Methods and Results

In the MD simulations, there are four methods to compute the thermal conductivity for MD systems: (1) NEMD: energy is applied to the hot region, and an equal amount of energy is subtracted from the cold area in the simulation cell. The heat flux is monitored between different temperature layers; (2) NEMD: Energy is added or subtracted in two regions, and the temperature difference of the intermediate region is monitored; (3) rNEMD: The kinetic energies of two atoms in different layers are swapped, and the temperature gradient is monitored; (4) EMD: the heat flux can be computed from the fluctuations of per-atom potential and kinetic energies, as well as the stress tensor. It is common in NEMD (non-equilibrium MD) for calculating the thermal conductivity of systems to impose the temperature gradient and the responded heat flux is measured. However, a reverse non-equilibrium MD (rNEMD) algorithm is used in the Muller-Plathe method [47]. The heat flux is applied in the system and the temperature gradient is measured. When the heat flux is imposed in the simulation cell, which is divided into N slabs (N is an even number, 20 was used in this study) with identical thickness. Energy transfer (Figure 10a) is produced from hot to cold slabs through the z-direction and it causes the temperature difference (Figure 10b) between these two slabs. Velocity exchange occurs in two particles, and the energy conservation is satisfied. The formulas for thermal conductivity and heat flux are shown in Equation (9). In addition, the Avogadro constant was used to complete the unit conversion due to different scales from microscopic to macroscopic states based on Equation (9). The mass/volume effect in the MD simulation was also considered in the unit conversion. The results of thermal conductivity of the xGNP modified and control asphalt binder models are presented in Figure 10c.

$$J=-\lambda \nabla T\text{ and }\lambda =-\frac{{{\displaystyle \sum }}_{transfers}\frac{m}{2}(v_{hot}^2-v_{cold}^2)}{2t{L}_x{L}_y\langle \partial T/\partial z\rangle }$$

(9)

where $\nabla T$ is the temperature gradient (scalar) in the simulation cell; J is the energy transferred (scalar) through the surface of layers; $\lambda$ is the thermal conductivity; t is the simulation time; $v_{hot}$ is the velocity of the hot particle; $v_{cold}$ is the velocity of the cold particle; m is the identical mass of particles; $L_x$ is the length of the simulation box in the x-direction; and $L_y$ is the length of the simulation box in the y-direction; and $\partial T/\partial z$ is the temperature gradient in the z-direction.

Figure 10a shows the cumulative energy input in different molecular binder systems at a temperature of 298.15 K. The control asphalt binder system has a relatively low energy input in contrast to the energy of the xGNP modified asphalt binder system due to the small number of molecules and low volume in the control binder system. The energy of the xGNP modified asphalt binder system is four times greater than that of the control asphalt binder system, and it is the same as the mass ratio of the xGNP modified asphalt binder system to the control binder system. The temperature variation of the xGNP modified asphalt binder model is also lower than that of the control binder system due to the large number of molecules. Figure 8d shows the temperature difference in different asphalt binder systems after the input of the heat flux at a temperature of 298.15 K. The variation in the temperature difference of the xGNP modified binder model is less than that of the control binder model, as well as the temperature variation in the MD simulation. It is likely that more molecules in the system result in better stress and heating responses and produce a stable system. The cumulative energy and temperature difference of the control and xGNP modified asphalt binder models at different temperatures were calculated using MD simulations, and based on Equation (9), the thermal conductivity results of the asphalt binder models are shown in Figure 11 (next section).

6.3. Comparison of Thermal Conductivity of the Control and xGNP Modified Asphalt Binder Models

Figure 11 demonstrates the thermal conductivity results through MD simulations. The addition of xGNP model in the asphalt binder model increases the thermal conductivity of the modified asphalt binder model. It is consistent with the laboratory data. It is reasonable that the thermal conductivity of the control and xGNP modified asphalt models at room temperature is around 0.275 W/m·K and 1.146 W/m·K, respectively, compared to the reference data of graphite modified asphalt binders [22] from 0.396 W/m·K to 0.934 W/m·K. There is an insignificant difference between the laboratory data and MD simulation results. The xGNP particles in the control binder improves the thermal conductivity of the modified asphalt binder from the experimental data. The same trend of thermal conductivity is also observed in the data from the MD simulation after the addition of the multi-layer graphite xGNP model in the control asphalt binder model. The thermal conductivities of the control and xGNP modified asphalt binders increase with the increase in temperatures of the experimental tests, and the thermal conductivities of the control and xGNP modified asphalt binder models also increase by increasing the temperatures of the systems. However, there are minor differences between the experimental data and MD simulation results. There may be a few reasons for this: (1) in the preparation of the samples and laboratory testing, the xGNP particles in the modified asphalt were not perfectly dispersed in the tested area due to the mixing method and not due to operational errors, and this causes inhomogeneous heating of the modified asphalt during testing; (2) the test area for thermal conductivity is relatively small; (3) the multi-layer graphite xGNP model does not fully represent the xGNP particles in the asphalt binder matrix for the calculation of thermal conductivity, and there are some improvements needed for models of xGNP particles and asphalt binder. After the analysis of laboratory and MD data, it is confirmed that xGNP particles can improve the thermal conductivity of asphalt, and the multi-layer graphite xGNP model can also enhance the thermal conductivity of the control asphalt binder model. The trend in temperature versus thermal conductivity of the MD simulation results is the same as the trend in the experimental data.

7. Discussion and Conclusions

The MD model of multi-layer graphite xGNP nanoplatelets was created and used for the investigation of the effect of modification on the control asphalt binder model. The control asphalt binder model was composed of three components: asphaltenes, aromatics, and saturates at a certain ratio. The xGNP modified asphalt binder model was generated from the addition of the xGNP model in the control asphalt binder model. The conjugate gradient method and PPPM were used for energy optimization, and the Savitzky-Golay filter was used to smooth data. The Amber Cornell Extension Force Field and ESP charge were used in these asphalt models, and the physical properties of the MD binder models were calculated including density, the glass transition temperature, viscosity, and thermal conductivity. The following conclusions may be made.

(1) The densities of these asphalt binder models were computed, and the addition of the multi-layer xGNP model increased the density of the xGNP modified binder model compared to that of the control binder model. The molecular number in MD systems significantly affects the data variation for density calculation. The density of MD asphalt binder systems decreases with the increase in temperatures.

(2) The glass transition temperature of the xGNP modified asphalt binder model is around 250 K, and it is better than the results of the reference, 298 K–358 K [18]. This glass transition temperature is better than previous results (around 300 K [24]) for the control asphalt binder model, because it is the same as the glass transition temperature of SHRP asphalt binders, around 250 K, from laboratory results [42].

(3) The Muller-Plathe method was used to calculate the viscosity of the control and xGNP modified asphalt binder models. The 20 layers in the MD asphalt models were separated for this calculation. The addition of xGNP particles in the control asphalt binder matrix improves viscosities of the modified asphalt at different temperatures, and the same effect of multi-layer xGNP models in the control asphalt binder model was observed. Compared to the experimental viscosities of the xGNP modified asphalt binder, the viscosities of the MD simulation is close to the experimental results at the temperatures of 403 K, 423 K, and 443 K. The relationship between viscosities and temperatures in the data of the MD simulations is also the same as that of the laboratory results.

(4) The experimental data shows that the xGNP particles in the control asphalt binder increase the thermal conductivity of the modified binder at room temperature. During the calculation of thermal conductivity, the Muller-Plathe method was used in these MD simulations, and the multi-layer xGNP model in the control binder model also improves the thermal conductivity of the control binder model at room temperature. The thermal conductivities of the control and xGNP modified asphalt binders increase with increasing temperatures, and the same trend is observed in the data of MD simulations.

Therefore, the multi-layer xGNP graphite particles in the asphalt binder can improve viscosity and thermal conductivity of the asphalt binder, and the xGNP model in the control asphalt binder model can also enhance the density, glass transition temperature, viscosity and thermal conductivity of the control binder model. It is obvious that the same trend of experimental data and MD results is observed during the testing and MD calculations of different properties of asphalt. It is likely that the xGNP particles can be utilized and generalized for pavement construction and heat sinks. The contributions of this paper include (1) the use of the xGNP graphite particles to enhance the performance of the asphalt binder; (2) the generation of the xGNP model for the modification of the asphalt model; (3) the application of the Muller-Plathe method to compute the thermal conductivity of the asphalt models; and (4) the use of the correlation analysis to reveal the linear relationship in MD simulation data. In addition, more properties of the xGNP modified asphalt binder and its models will be tested and calculated for future research.