Four-component united-atom model of bitumen
We propose a four-component molecular model of bitumen. The model includes realistic chemical constituents and introduces a coarse graining level that suppresses the highest frequency modes. Molecular dynamics simulations of the model are be carried out using Graphic-Processor-Units based software in time spans in order of microseconds, and this enables the study of slow relaxation processes characterizing bitumen. This paper focuses on the high-temperature dynamics as expressed through the mean-square displacement, the stress autocorrelation function, and rotational relaxation. The diffusivity of the individual molecules changes little as a function of temperature and reveals distinct dynamical time scales as a result of the different constituents in the system. Different time scales are also observed for the rotational relaxation. The stress autocorrelation function features a slow non-exponential decay for all temperatures studied. From the stress autocorrelation function, the shear viscosity and shear modulus are evaluated at the highest temperature, showing a viscous response at frequencies below 100 MHz. The model predictions of viscosity and diffusivities are compared to experimental data, giving reasonable agreement. The model shows that the asphaltene, resin and resinous oil tend to form nano-aggregates. The characteristic dynamical relaxation time of these aggregates is different from the homogeneously distributed parts of the system, leading to strong dynamical heterogeneity.
Refined bitumen (or asphalt) is a residual product of the refinery process of crude oil, and is highly viscous at room temperatures. The chemical composition of bitumen is very complex in that it is composed of up to 10 relatively large molecules wiehe:1996; rogel:ef:2003; artok:ef:1999; murgich:ef:1996 with molar masses above 200 g/mol. Bitumen is not chemically unique in the sense that its composition depends on the crude oil source, the age of the bitumen barth:1962, and possible chemical modifiers zhang:ef:2008. The functional chemistry includes heteroatoms like sulfur and nitrogen, cyclic structures, aromatics compounds and saturated hydrocarbons barth:1962, and these enter larger and complicated molecular structures. A correct and detailed chemical categorization is not feasible, and bitumen is often characterized simply by its rheological properties such as bulk and shear moduli and penetration depth reed:2003. A fundamental understanding on the molecular level is, however, necessary if one wishes to understand and control the mechanical properties.
The aim of the present work is to do just that, i.e., we wish to study dynamical properties of bitumen through molecular modelling. A further impetus of this work lies on the definition and characterization of a “Cooee-bitumen model” with the aim to bring new insights on how to reduce the rolling resistance generated between the tyre and the pavement schmidt:2012. Molecular simulations has previously been used to investigate bitumen properties, for example, the density and structure in asphaltene systems rogel:ef:2003 and molecular alignment murgich:ef:1996. Recently, Zhang and Greenfield (Z & G) published a series of papers zhang:ef:2007; zhang:ef:2007:2; zhang:jcp:2007; zhang:ef:2008; zhang:jcp:2010, wherein they study structural, thermodynamical and dynamical properties of different bitumen models with and without polymer additives. To this end the authors used atomistic molecular dynamics simulations. Their bitumen models are based on a tertiary mixture where each chemical component represents a distinct constituent : high-weight asphaltenes, aromatic maltenes (resin) and saturated maltenes (-alkanes). In their simulations they applied an all-atom force field that includes the interactions between both hydrogen and carbon and hydrogen and sulfur, as well as the hydrogen-carbon bonds. While this allows for standard and well tested force field parametrizations and gives very detailed description of the system, it also includes very fast modes which are unlikely to couple to the long relaxation times characterizing viscous liquids like bitumen. The inclusion of these fast modes necessitates a small integration time step in the numerical intregration of Newton’s equations of motion, and Z & G were able to simulate a time span in the order of 10 nsec., even using parallelized molecular dynamics software zhang:jcp:2010. For the aromatic maltene constituent Z & G used dimethylnaphtalene or closely related molecular structures. This relatively low weight molecule has a boiling point of around 278 C lide:1991, which is much lower than that of the other constituents and dimethylnaphtalene is therefore likely to be removed from the crude oil during the refinery process at approximately 500 C, see Ref. barth:1962.
We propose a new four-component bitumen model, based on the Hubbard-Stanfield classification hubbard:1948; rostler:1965 The model is coarse grained such that the hydrogen atoms are not explicitly included, but only implicitly through parametrization of the particle interactions. This allows for larger time steps and, more importantly, it reduces significantly the number of bond-, angle- and pair-forces one needs to evaluate in each time step. The molecular dynamics simulations are carried out using the Graphical-Processor-Unit (GPU) based molecular dynamics software package RUMD bailey:2012. In this way time spans in the order of sec. are accessible using a single graphics card. We are thus able to study much slower relaxation phenomena in bitumen compared to what was previously possible.
In this paper we focus on the high-temperature dynamical properties of the model and identify the temperature where the relaxation time becomes very long. Specifically, we investigate characteristic relaxation times for collective, chemical and single molecule quantities. From this we show that the dynamical properties of the model are quite different from those of Z & G’s model, in that the dynamics are much slower and changes only little in the high temperature regime. For relatively low temperatures, the dynamics is characterized by complex glass-like behavior, where the relaxation times widely exceed the accessible simulation times.
The paper is organized as follows. Section II describes the parameters chosen in the model and the simulation method. Section III contains results and discussion of diffusivity, viscous properties and rotational dynamics of the model bitumen. Finally, a summary of the paper is provided in section IV.
Ii Model and parametrization
As stated in the introduction, a detailed chemical model for bitumen is not possible. We do therefore not target any specific bitumen, but aim at capturing the characteristic properties of a typical bitumen. We conjecture that an explicit inclusion of the hydrogen atoms is not important for our modelling purpose, even if it may be useful for other specific purposes (see for example Ref. toxvaerd:1990). We therefore coarse grain the chemical structures which enables us to investigate longer time dynamics.
The bitumen constituents are characterized using different separation techniques rostler:1959; barth:1962; reed:2003. We will here use a general classification scheme by Hubbard and Stanfield hubbard:1948; rostler:1965, see Fig. 1. We let unsaturated hydrocarbons and the aromatic/cyclic compounds without any heteroatoms like sulfur fall into the catagory resinous oils barth:1962. Also, we interpret the class non-resinous oils to be paraffin wax, i.e. staturated hydrocarbons. This motivates a four-constituent model composed of (i) asphaltene, (ii) resin, (iii) resinous oil and (iv) saturated hydrocarbon. The model also follow the SARA classification (Saturates-Aromatics-Resins-Asphaltenes)astm:2007; although the unsaturated hydrocarbons and the aromatic/cyclic compounds fall under the catagory ’Aromatics’ in the SARA classification, we prefer to use the term resinous oil to follow Hubbard and Stanfield.
A single molecule is used to represent each constituent. Based on NMR studies Artok artok:ef:1999 proposed an asphaltene structure, shown in Fig. 2, which is adopted here. This structure was also used by Z & G which they call asphaltene 1. For the resinous oil and resin, we use the structures by Rossini et al. rossini:1953 and Murgich et al. murgich:ef:1996, respectively. In Ref. zhang:ef:2007, Z & G argued that the saturated hydrocarbons can be modelled by docosane (-C), since this molecule represents the average chain length found by Storm et al. storm:ef:1994. It is also used here. All four molecular structures are shown in Fig. 2.
We represent the methyl (CH), methylene (CH), and methine (CH) groups by the same Lennard-Jones particle. This is then a united atomic unit (UAU) with mass 13.3 g/mol and Lennard-Jones parameters = 3.75 Å and are estimated from the OPLS force field jorgensen:jacs:1988. The sulfur atoms are represented by the same Lennard-Jones interaction, but with a mass of 32 g/mol.
The force field is given through the total potential energy of the system
The force field thus includes pair interactions as well as interactions due to (flexible) bonds, angles and dihedral angles. The bonds, angles and dihedral angles are specified via bonds between the UAUs. UAUs belonging to the same bond, angle and dihedral angle do not interact via the pair interaction. From the force field it can bee seen that the model has a large parameter space. To reduce the number of free parameters, we assume that all bonds have the same length Å. This particular value of is chosen to be the average value of all the bonds in the system ( Å) computed from energy minimization of the individual molecules avogadro. This is, of course, a simplification since energy minimization of, for example, the asphaltene gives the single bonded units a bond length of around 1.50-1.54 Å and the double bonded units a bond length 1.38-1.41 Å. These values again depend on the specific molecule and amount to very many different bond length, so we choose an average value for simplicity. The same argument lies behind the choices of the angles and dihedral angles, but we here use two distinct angles and three distinct dihedral angles. We use the Generalized Amber Force Field (GAFF) wang:jcc:2004 parametrization for bond and angle parameters, see Table 1. For the saturated hydrocarbons and the side chains in the asphaltene and resin, we use the Ryckaert-Bellemans parameters catlow:1989. These parameters are optimized with respect to liquid butane, but they capture the important cis, trans and gauche configurations reasonably well. We are thereby left with a single free parameter, namely the zero-force coefficient for dihedral angles and 0 degrees in the ring structures.
The parameter is estimated from one structural criterion, namely the alignment between the aromatic rings in a one-component system of 1,7-dimethylnaphtalene (also used by Z & G). This particular criterion is chosen because of the simple molecular structure of 1,7-dimethylnaphtalene. Let the vector be defined as the cross product , where and are two vectors parallel with two bonds in one of the ring structure of the molecule, as shown in Fig. 3. Thus is a vector normal to the ring structure. Likewise, the vector is defined as normal to the second ring structure. The angle between and then defines the ring alignment angle, .
Figure 3 shows the histogram of the alignment angle for 1,7-dimethylnaphtalene at mass density 994 kg/m and temperature K. The values kcal/mol and kcal/mol were chosen for the and degrees zero-force dihedral angles, respectively, as they lead to a good agreement between the all-atom model used by Z & G and the united-atom model used in the present work. Table 1 lists all parameters used in the model.
ii.2 Simulation Method
In all simulations we used asphaltene, resin, resinous oil and docosane molecules, giving a total of 3114 UAUs. The mass fraction corresponds closely to Z & G’s mixture 1 zhang:ef:2007, for which dimethylnaphtalene represents both the resin and resinous oil, but we used twice as many molecules. Initially, the molecules center-of-mass are arranged on a simple low-density lattice and the system is compressed (while integrating the equation of motion) to the desired density. We use a target density of kg/m, and the quench is carried out for a temperature of K. The system was coupled to a Nosé-Hoover thermostat nose:mp:1984; hoover:pra:1985 and the equations of motion were integrated using a leap-frog integrator scheme frenkel:1996. After the compression the volume is kept fixed in order to perform NVT ensemble simulations. A first production run is carried out at K over a total time span of 0.85 sec. During this run the system reaches equilibrium (see the discussion below). The final configuration from this state point is used to reduce the temperature from K to and 452 K at a rate of 2.2 K per nanosecond. In this way we study the systems at four different temperatures. For K production runs of 1.7 sec. were performed. For all temperatures the production runs were divided into 10 intervals of same time span (0.085 sec. and 0.17 sec. for = 603 K and 528 K, respectively), which enables us to follow the equilibration of the system. We only use the last eight sample intervals for data analysis in order to avoid the initial fast relaxations in the system.
In order to highlight the different dynamics we also studied the system at a relatively low temperature and far from equilibrium. To this end a dilute and high-temperature system is rapidly quenched to the target density and K. This system is also simulated over a time span of 1.7 sec.
To determine the integration time step, the energy drift was investigated by performing a series of NVE simulations of the bitumen mixture at the desired density and for different temperatures over one million time steps. For low temperatures and time steps larger than 3.4 femtoseconds, we noticed a slight systematic increase in the energy fluctuations (quantified by the energy standard error). To be conservative, we chose a time step of 1.7 femtoseconds for 528 K and 0.85 femtoseconds for T=603 . These time steps are around 35-70 % larger than those used by Z & G zhang:jcp:2007, and are made possible due to the coarse graining. This enables us to simulate a time span of up to 1.7 sec. using one billion time steps.
Iii Results and discussion
This study analyzes (i) the average diffusivity of the four different chemical components, (ii) the collective viscous properties through the stress (or equivalently the pressure tensor) autocorrelation function and (iii) the orientational dynamics of single resinous oil and asphaltene molecules.
The diffusivity of the molecules is evaluated through the mean-square displacement of the center-of-mass of each molecule. If the center-of-mass position of molecule is denoted , the mean-square displacement of molecules of a certain type is defined as hansen:2006: , where index runs over molecules of that type. In the normal diffusive regime
where is the diffusion coefficient. In Fig. 4, the mean square displacements of the four components is plotted for 452 K and 603 K.
It is seen that the molecules are in the normal diffusive regime at both temperatures, for times longer than 1 nanosecond. The corresponding diffusion coefficients are plotted in Fig. 5 for different temperatures. Upreti and Mehrotra upreti:2002 reported diffusivities in the order of 0.1 - 1.0 m/s for small weight gases in bitumen over a temperature range of 20-200 C, which is in reasonable agreement with our results in view of the fact that we study higher weight molecules. From Fig. 5 it follows that , where is the diffusivity for docosane, for resinous oil, for resin and for asphaltene. This, of course, agrees with the fact that the mobility decreases with increasing molecular size.
Both asphaltene and docosane have diffusion coefficients that are lower than the values reported by Z & G zhang:jcp:2007, who found m/s and m/s at K. This discrepancy is likely due to the current model replacing 1,7-dimethylnaptalene with the resin and resinous oil compounds. These relatively large structures hinder the motion of the compounds significantly compared to the smaller 1,7-dimethylnaphtalene. This may also explain another difference between the two models: we find a large difference in diffusivity between the docosane and asphaltene () whereas Z & G’s model shows a difference of only a factor of 5.
The mean-square displacement for a system rapidly quenched to the temperature K was also studied and is plotted in Fig. 6. In order to highlight the details and as opposed to what was done in Fig. 4, no sample averaging is done here. The docosane molecules still exhibit an approximate diffusive motion whereas the larger resin and resinous oil molecules perform vibrational motion around a given position (arrested), and then on average make a rapid movement (jump). The resin molecules eventually enter an approximately diffusive regime. The arrest and jump motions are both random events, which are also observed for simple models of glasses berthier:rmp:2011; edinger:jcp:2012.
This behavior shows the inherent heterogeneous dynamics present in the bitumen model where the asphaltenes are immobile, on the time scales shown in the figure, whereas the docosanes approaches normal diffusivity.
iii.2 Viscous properties
To study the viscosity we use the Irving-Kirkwood irving:1951 expression for the pressure tensor
where is the momentum of molecule , its mass, is the center-of-mass displacement vector between molecules and and is the total force acting between molecules and . The force between molecules is the sum of the atomic forces todd:ms:2007, i.e. , where atom is a constituent atom in molecule and in molecule . is the system volume. The pressure tensor defined above is in general not symmetric todd:ms:2007, and the traceless symmetric part can be extracted using . From this we can define the shear stress autocorrelation function
where runs over the , and pressure tensor elements and denotes an ensemble average. The ensemble average is considered equivalent to a time average. Each sample interval is divided into sample windows for with time spans of around 2 nsec. The product is computed for a time , nsec. and for different initial times within the specific window. The average is performed over the initial times within a period and over the different periods, leading to the stress autocorrelation function for time going from zero to 2 nsec. Comparing the results of the 8 successive samples in a production run enables us to probe the equilibration process.
In Fig. 7 the stress autocorrelation function is plotted for K and for =452 K.
For higher temperatures, (Fig. 7 (a)), the correlation function is fully decayed (note the logarithmic scale of the abscissae) and time invariant within statistical error, as shown by the equivalent results obtained for different samples. For low temperatures (Fig. 7 (b)), it is seen that the autocorrelation function does not relax fully and that the 8 successive samples in a production run does not lead to the same result, i.e. the relaxation time characterizing collective dynamics is larger than the sample window for the stress autocorrelation function of 0.17 sec. The 2 nsec. period used here is approximately ten times longer than what was used by Z & G to study their bitumen model (including equilibrium properties) at temperatures 298-443 K zhang:jcp:2007.
For highly viscous liquids the relaxation processes are often quantified by non-exponential functions cavagna:2009. In particular, relaxation at later times is frequently described by a stretched exponential function
where is the characteristic relaxation time and . Figure 7 shows the best fits of Eq. (5) to the averaged data (red line). It is seen that the stretched exponential form fits data well even for small times for K. For this temperature the coefficient is around 0.20 which may be interpreted as an indication of strong spatial dynamical heterogeneity cavagna:2009. We discuss this point further below.
The frequency-dependent viscosity is given by the stress autocorrelation function as follows hansen:2006
This expression can, of course, only be used for temperatures where the stress autocorrelation function is fully relaxed and time invariant; in Fig. 8 (a) we have plotted the running integral indicating that this is the case within statistical uncertainty for =603 K.
Figure 8 (b) displays for =603 K; it is observed that the fluid response is purely viscous for MHz. The imaginary part of the viscosity exhibits a peak at around 0.5 GHz, which corresponds to a characteristic relaxation time of approximately 0.2 nsec. The frequency-dependent shear modulus is plotted in Fig. 8 (c). Two distinct regimes are clearly observed; at low frequencies which is the zero-frequency behavior and at high frequencies . This later frequency dependency is different from models for simple polymers bird:1987, so great care should be taken when applying the theoretical predictions from polymer science to the field of bitumen.
From Fig. 8 (b) we can extract the zero-frequency viscosity for = 603 K. In order to compare the predicted viscosity with experimental data available for lower temperatures we extrapolate these to = 603 K. In Fig. 9 this is done for the standard SV11274-70/100 bitumen, which is measured in this work, see Ref. erik:SV for a short description of the measurement technique. Only data for the three highest temperatures have been used for the fitting procedures. Data were fitted to (i) first order polynomial (pure Arrhenius behavior) (ii) the Vogel-Fulcher-Tamman (VFT) function where and are fitting parameters, and (iii) a power law (shown). The power law is found to fit data best. The model predicts the viscosity reasonably well for SV11274-70/100 bitumen; but not for PG 64-22 bitumen. As we mentioned above, it is not the purpose of the model to fit a specific bitumen, and the fact that the model cannot capture the dynamical properties of different bitumens is expected. Note that we used the stress autocorrelation function to calculate the viscosity, which is only possible for high temperatures. Other approximative methods bird:1987; zhang:jcp:2007 may be used to estimate the viscosity; however it is not clear whether these methods can be applied to bitumen mixtures.
iii.3 Rotational dynamics
Above it was observed that while the molecular mean-square displacements indicate a well-defined diffusive regime for all temperatures studied, the stress autocorrelation functions showed a not fully relaxed system for times around 0.17 sec. To better understand this behavior of the model we will study the rotational dynamics for resinous oil and asphaltene. Rotation is quantified through the dynamics of the normalized orientational vector normal to one of the ring structures. As in Sect. II, the vector is defined as the cross product of two vectors pointing along a chemical bond as illustrated in Fig. 2, where the numbers indicate the atom indices involved in the definition of .
In Fig. 10, projections of the orientational vector are plotted as it evolves over time. It is evident that the smaller resinous oil goes through all possible orientations (to a good approximation) for all temperatures within the 0.85 to 1.75 sec. sample time. This is not the case for the larger asphaltene. Especially for lower temperatures the orientation of this molecule is confined to a small subset of orientations. Naturally, this subset increases as a function of time.
The first-order rotational correlation function, hansen:2006, is plotted in Fig. 11. The rotational relaxation appears to be governed by multiple relaxation times because of the nonsymmetric molecular structures. From it is possible to define a characteristic relaxation time . For example, we obtain nsec. for the resinous oil at 452 K. Again this is evidence of a much slower dynamics compared to Z & G’s model where nsec. for the asphaltene at 443.15 K. In combination with Fig. 4 it can be concluded that while, on average, the asphaltene molecules enter the diffusive regime, the single molecules keep their orientation relatively fixed. Such a decoupling between the rotational and diffusional dynamics is known from many experiments on viscous liquids, see Ref. edinger:jcp:2012 and references therein.
The snapshot shown in Fig. 12 highlights this picture; the asphaltene molecules tend to align and form nano-aggregates, but occasionally they leave the aggregate and then later attach at a new position or aggregate. This dynamical behavior reveals that the diffusivity of the asphaltenes is determined by the mobility of just a small fraction of the molecules at any instant in time. The same remark applies to resin and resinous oil that also participate in the nano-aggregate formation. However, resin and resinous oil evolve at a different time scale. This long-time fixed orientation and clustering naturally leads to a very slowly decaying stress autocorrelation function. We conclude that the dynamics of the model is dominated by dynamical heterogeneities.
In this paper we have proposed a new molecular model of bitumen. Compared to Z & G’s model, we have replaced the dimethylnaphtalene with two higher-weight molecular structures representing resin and resinous oil molecules. The new model is a united-atomic-unit model, which suppresses the high-frequency modes and enables long simulation times. Using GPU-based software we are able to simulate time spans in order of sec., rather than nsec. zhang:ef:2007:2. This makes it possible to study additional slow relaxation processes characterizing bitumen.
For the density and temperatures studied here the model is characterized by slow dynamics, which changes little over a relatively large temperature range. The model predictions are in reasonable agreement with experimental data for diffusion and viscosity. Better agreement with specific bitumen solutions can be obtained by model re-parametrization and by choosing the state point carefully.
The model captures the heterogeneous dynamics since the asphaltene, resin and resinous oil form nano-aggregates. Murgich et al. murgich:ef:1996 and Z & G zhang:ef:2007:2 have reported that the asphaltene molecules tend to align, which supports this picture. The fact that asphaltenes from nano-aggregates is also known from experimental work yen:ac:1961; mullins:ef:2012 The characteristic dynamical relaxation time of these localized aggregates are different from the more homogeneously distributed parts of the system leading to a high degree of dynamical heterogeneity. This is also evident in the stretched exponential fit of the stress autocorrelation function. In Fig. 5 one observes separate time scales, the slowest one being due to the dynamics of larger asphaltene molecules. This indicates that the asphaltene nano-aggregate dynamics may be studied in more details on a longer time scale using coarse graining techniques like Brownian dynamics. This may form the foundation for understanding the dynamical heterogeneity of bitumen.