RHEOLOGICAL AND THERMODYNAMIC PROPERTIES OF MODEL ASPHALT USING FFT AND GROUP CONTRIBUTION METHODS

Asphalt is an amorphous material whose mechanical performance relies on viscoelastic responses to applied strain or stress. High quality pavement can enhance the commute experience and also can reduce the risk of accidents as well as preventing air pollution. Furthermore, it also can reduce the need for expensive reconstruction and maintenance of roads. Having a broad understanding of asphalt chemistry is the key to understand and improve mechanical properties of this material. Although asphalt is a widely seen material, our knowledge about its molecular structure and properties is limited. Molecular simulations of asphalts can be exploited to infer how the actions of individual molecules contribute to the nanoscale mechanical behavior of a model system. In this work, chemical composition and its effect on the viscoelastic properties of asphalts have been investigated by computing complex modulus from molecular dynamics simulation results for two different model asphalts (Zhang, L., & Greenfield, M. L. (2008). Energy Fuels, 22(5), 3363–3375 [ZG08] and Li, D. D., & Greenfield, M. L. (2014). J. Chem. Phys., 140(3), 034507 [LG14]) whose compositions each resemble the Strategic Highway Research Program AAA-1 asphalt in different ways. Results from equilibrium molecular dynamics simulations these have been interpreted by converting the stress relaxation modulus G(t) to the complex modulus and phase angle delta. Complex modulus at different temperatures have been calculated using fast Fourier transform to study the effect of temperature on viscoelastic properties of asphalt models. Because of inherent noise, any comparison was challenging. To remove the noise, signal processing techniques have been exploited. Signal processing techniques enhanced the clarity of the results and enabled removing the noise from the modulus results. The LG14 system contains larger molecules, and its results have shown a very good agreement with the low and high frequency scaling limits of the Maxwell model within the frequency ranges spanned by the molecular dynamics simulations, while results for the ZG08 model asphalt only follow the high frequency scaling limits of the Maxwell model. A Black plot or van Gurp-Palman plot of complex modulus vs. phase angle for the LG14 system suggests some overlap among results at different temperatures for less high frequencies, with an interdependence consistent with the empirical Christensen-Anderson-Marasteanu model. Both model asphalts are thermorheologically complex at very high frequencies, where they show a loss peak that appears to be independent of temperature and density. While molecular simulation can provide relatively accurate results, it is computationally expensive. To provide a faster approach to determine the relationship between the asphalt chemistry and macroscale properties, an equation of state (EOS) can be used. Most equations of state need critical properties which are not always available for all compounds. Since most molecules within the ZG08 and LG14 model asphalts are hypothetical molecules and no experimental data are available for them, a group contribution method by Nannoolal et al. (Nannoolal, Y., et al. (2004). Fluid Phase Equilib., 226, 45–63; Nannoolal, Y., et al. (2007). Fluid Phase Equilib., 252(1-2), 1–27; Nannoolal, Y., et al. (2008). Fluid Phase Equilib., 269(1-2), 117–133) has been applied for all molecules in both model asphalts to determine their critical properties and acentric factors. These provide the required properties to use in a common cubic EOS, e.g., Peng-Robinson (PR), to predict phase behavior of multicomponent systems. The results for densities and thermal expansion have been compared to experimental results. Since cubic equations of state do not predict liquid density precisely, COSTALD volume translation has been applied to improve density predictions. The results from PR EOS were quite good but the results from molecular dynamics were better. This fact should not devalue the EOS results because an EOS calculation is much faster and takes a fraction of one step of molecular simulation time. The result from PR without COSTALD volume translation for thermal expansion was even better than results from molecular dynamics. By knowing asphalt phase behavior, one can predict its performance in the real world. For example, rutting of asphalt is usually related to the waxy molecules within asphalt mixture. Precipitation of waxes can cause the cracking of asphalt pavement as well as increase asphalt mixture viscosity during compaction of asphalt pavement. To predict this phenomenon, precipitation of squalane was investigated. The chemical potential of pure squalane (as a waxy compound) and squalane in the multicomponent system have been estimated. To calculate the Gibbs free energy required to study phase behavior of asphalt, the critical properties of molecules were calculated using the Nannoolal et al. correlations. The chemical potential of pure squalane somewhat below its melting point was lower than its chemical potential in the multicomponent system, which suggests that squalane will precipitate below its melting point. The temperature that the wax precipitates can help to choose the right asphalt binder for a specific location. In total, the approaches used here provide ways to model how differences in asphalt can affect thermodynamic and mechanical properties. Signal processing can help remove the noise from data which leads to a better understanding and interpretation of the results. Having equation of state parameters for a system will help to predict phase stability of the asphalt systems at different conditions.

Equilib., 269(1-2), 117-133) has been applied for all molecules in both model asphalts to determine their critical properties and acentric factors. These provide the required properties to use in a common cubic EOS, e.g., Peng-Robinson (PR), to predict phase behavior of multicomponent systems. The results for densities and thermal expansion have been compared to experimental results. Since cubic equations of state do not predict liquid density precisely, COSTALD volume translation has been applied to improve density predictions. The results from PR EOS were quite good but the results from molecular dynamics were better. This fact should not devalue the EOS results because an EOS calculation is much faster and takes a fraction of one step of molecular simulation time. The result from PR without COSTALD volume translation for thermal expansion was even better than results from molecular dynamics. By knowing asphalt phase behavior, one can predict its performance in the real world. For example, rutting of asphalt is usually related to the waxy molecules within asphalt mixture. Precipitation of waxes can cause the cracking of asphalt pavement as well as increase asphalt mixture viscosity during compaction of asphalt pavement. To predict this phenomenon, precipitation of squalane was investigated. The chemical potential of pure squalane (as a waxy compound) and squalane in the multicomponent system have been estimated. To calculate the Gibbs free energy required to study phase behavior of asphalt, the critical properties of molecules were calculated using the Nannoolal et al. correlations. The chemical potential of pure squalane somewhat below its melting point was lower than its chemical potential in the multicomponent system, which suggests that squalane will precipitate below its melting point. The temperature that the wax precipitates can help to choose the right asphalt binder for a specific location.
In total, the approaches used here provide ways to model how differences in asphalt can affect thermodynamic and mechanical properties. Signal processing can help remove the noise from data which leads to a better understanding and interpretation of the results. Having equation of state parameters for a system will help to predict phase stability of the asphalt systems at different conditions. First ps of the stress relaxation modulus of both AAA-1 model asphalts. Temperature suffixes and þ indicate the ZG08 and ZG08p systems; [1,2] others refer to the LG14 system. [3] . . . . 90  [1,3] Squares indicate experimental data [5] for a

Frequency Analysis of Stress Relaxation Dynamics in Model Asphalts
The following manuscript has been published in The Journal of Chemical Physics.

Abstract
Asphalt is an amorphous or semi-crystalline material whose mechanical performance relies on viscoelastic responses to applied strain or stress. Chemical composition and its effect on the viscoelastic properties of model asphalts have been investigated here by computing complex modulus from molecular dynamics simulation results for two different model asphalts whose compositions each re- thermorheologically complex at very high frequencies, where they show a loss peak that appears to be independent of temperature and density.

Introduction
Asphalt has been a widely used material since ancient eras. Roberts et al. [1] cite the use of asphalt in shipbuilding in Sumeria as long ago as 6000 B.C, for example. Despite this long history, a thorough description of how asphalt chemistry directly impacts its mechanical properties has been incomplete. Progress in understanding asphalt [2,3,4] and asphaltene [5] chemistry is enabling advances that link asphalt chemistry to physics and mechanics. Achieving such a deep understanding of asphalt mechanical properties can help to understand pavement mechanical behavior and to improve pavement performance and durability. [6] The aim of this paper is to calculate viscoelastic properties of proposed model asphalt systems via frequency analysis of stress relaxation dynamics. Stress correlations from equilibrium molecular dynamics simulations [7,8] using two different model systems [9,10] that each represent the atomic composition of asphalt AAA-1 of the Strategic Highway Research Program (SHRP) [11] have been investigated by calculating storage and loss moduli from stress relaxation modulus data. Results are compared to simulation results of others and to experimental results available in the literature for real asphalts.
The SHRP revolutionized the mechanical tests that are employed to characterize asphalts and bitumens. [12] Tests such as viscometry and penetration were found to be inadequate for estimating asphalt performance. Instead, thermal and rheological tests such as complex modulus and viscosity were found to provide a more thorough understanding of asphalt characteristics, including estimates of asphalt mechanical performance and elasticity. [1] While magnitude and temperature dependence of complex modulus and viscosity describe many asphalt mechanical properties, [12] asphalts with comparable complex modulus under certain conditions show different performance in "real world" pavement applications because of independent behaviors of other mechanical properties. [1,6,13] Another challenge of asphalt is the composition of its components. Because different asphalts come from different crude oil sources, they have different chemical components and different properties. [1] Asphalt generally contains approximately 10 5 to 10 6 different molecule types. [14] It can be separated into at least three different fractions -asphaltene, resin, and maltene -that range from least non-polar and most viscous to most non-polar and least viscous. Asphaltene molecules are considered to have a great role in the viscosity of bitumen. [1,6] Asphalt nanoscale structure also affects its properties. [4] Aging asphalt changes its properties further. [2] Asphalts have been characterized since the SHRP by using complex modulus.
Storage modulus, G ′ (ω), is a measure of energy stored and recovered per cycle of sinusoidal deformation and can be related to stress relaxation modulus G(t) by [15] G e is the equilibrium tensile modulus and can be calculated by G e = lim t→∞ G(t), and it equals zero for viscoelastic liquids. [15] The angular frequency ω is defined using frequency as ω = 2πf . Loss modulus, G ′′ (ω), is a measure of energy lost per cycle of sinusoidal deformation and can be calculated by [15] G ′ (ω) and G ′′ (ω) can be merged to define complex modulus G * (ω) = G ′ (ω) + iG ′′ (ω), [15] which in asphalt rheology is represented by its magnitude and phase angle, [15] In an oscillatory shear experiment, |G * (ω)| corresponds to the ratio of stress and strain amplitudes. Phase angle δ(ω) corresponds to the difference in phase between stress and strain oscillations, and it represents the ratio between viscous and elastic responses during the shearing process. [15] Complex modulus and δ can also indicate glass transition (T g ) and secondary molecular relaxation. [15] Dynamic rheology tests are conducted so deformation data are in the linear viscoelastic region. Often the results can be expressed well by master curves that use the time-temperature superposition principle (TTSP). However many distress mechanisms such as failure, cracking, moisture damage, and tertiary creep are nonlinear. [16] Extra relaxation mechanisms that reflect morphological alterations can lead to TTSP not applying to low frequencies or high temperatures. [16] To use TTSP successfully, all contributing viscoelastic mechanisms must follow the same temperature dependence. If TTSP is applicable to a material, the material is called thermorheologically simple. [17,18] The Black (or van Gurp-Palman) plot, [19] which encompasses phase angle and complex modulus, can be used to test the applicability of TTSP.
In general, a relaxation log-log plot of G ′ (ω) and G ′′ (ω) vs. frequency has four typical regions: terminal zone, plateau zone, transition zone, and glassy zone. [15] The terminal zone corresponds to flow at low frequencies (i.e. long times). [15] One general characteristic of the plateau zone is a relatively small change in moduli over the frequency region. Another is that G ′ is always at least as large as G ′′ , with phase angle passing through a minimum. This property is characteristic of networks, such as entangled or cross-linked polymers, and may be missing if a network does not exist over the measurement time scale. [15] Furthermore, in polymer modified asphalts, the width of this plateau is a measure of compatibility of an elastomeric polymer and the asphalt. [16] Moduli increase by several order of magnitude in the transition zone from rubber-like to glass-like and remain higher as a function of frequency. Another characteristic of the transition zone is that at a given frequency, storage modulus rises rapidly with decreasing temperature, and the lower the chosen frequency, the steeper the rise. [15] A large maximum in phase angle occurs when the transition from rubber-like to glass-like consistency is crossed. [15] The temperature of the midpoint T M of the primary transition can be roughly specified by the inflection in G ′ . In the glassy region, the viscoelastic properties reflect local molecular motions that include rotation around backbone bonds, configurational rearrangement of side chains, and rotations of terminal groups. [15] Viscoelastic relaxation corresponds to altering the distributions of torsional fluctuations, bending and stretching of bonds, and of intermolecular arrangments over the time scale of mechanical deformation. [15] Local motions are studied by high-frequency viscoelastic measurements above T g rather than by low-frequency below T g . These two types of studies usually give closely related but not quite identical information. Above T g , the local structure changes with temperature far more than below T g , leading to a larger thermal expansion coefficient. As a result, the temperature dependence of the magnitude of relaxation may be different above and below T g . Normally, at higher temperatures molecules need less time to relax. [15] Molecular simulations have previously been applied to isolated asphaltenes and to bitumens, as recently reviewed. [20] Boek et al. [21,22] [28] have investigated the nanoaggregation of asphaltene in toluene and in heptane as well as asphaltene-asphaltene potential of mean force from molecu-lar dynamics simulations. Jian et al. simulated the effects of different side chain lengths for a model asphaltene with a 9-ring core and two ester substituents. [29] Our group proposed combinations of compounds to represent asphalt systems, focusing on asphalt as a system rather than asphaltenes as a solubility class. Initial work employed ternary systems. [30] A system with a larger number of more polar molecules was proposed next; [9] its molecules were later found to be too small. [7] A new set of yet larger molecules was proposed recently. [31,10] A coarse-grained four-component "Cooee bitumen" system was recently proposed by Hansen and co-workers. [32,33] Stress relaxations up to nanosecond time scales were shown for simulations at 452 and 603 K, with a long-time limit reached at the higher temperature. [32] Chemical aging was then simulated via changes in asphaltene and resin content. [33] The long durations of these simulations enabled independent estimates of stress relaxation modulus over widely separated time periods.
The work presented here provides further interpretations of our prior molecular dynamics simulations. [7,8]

Methods
Molecular simulations [34] evaluate statistical mechanical averages numerically, given a force field that represents intermolecular potential energy among molecules.
The work described here entails further post-processing analysis of results from molecular dynamics simulations that were described in earlier papers. [9,7,10,8] There, the optimized potentials for liquid simulations (OPLS) [35,36,37] all-atom force field was chosen because its parameters have been optimized for liquid state properties and because it works sufficiently well [38] for aliphatic and aromatic hydrocarbons as well as sulfur compounds. In those prior simulations, [9,31,10] a fixed number of molecules, volume, and temperature (NVT ensemble) and the Nose-Hoover [39,40] thermostat in molecular dynamics (MD) were used, with a time step ∆t of 1 fs and fluctuations sampled for at least 4 ns. The volume was equal to the average volume that resulted from simulations at constant temperature and pressure (1 atm). A single simulation was run at each state point because in our prior work we obtained similar results from multiple independent simulations.
Two model asphalt systems simulated in prior work had been devised so their chemical compositions resembled those of an asphalt termed "AAA-1" in the SHRP.
Simulation results from those two different model asphalt systems are investigated here. The earlier of these model asphalts [9] used relatively small molecules to represent the different SARA (saturate, aromatic, resin, asphaltene) fractions of real AAA-1 asphalt. In this work, we refer to that model asphalt system [9] as ZG08. It contained 18 linear C 22 alkane chains as a saturate; these were consistent with the size range invoked by Kowalewski et al. [25] 100 ethyltetralin molecules were used as a naphthene aromatic. 113 molecules of one-to three-ring aromatics based on thiophene, quinoline, and benzothiophene represented resins. The latter compounds were suggested by Masson and Lacasse [41] and from degradation studies by Strausz. [42] 9 molecules that used a structure proposed by Groenzin and Mullins [43] represented asphaltenes. The ZG08 system matched the atomic composition reported in experiments on real AAA-1 asphalt, [11] but it led to a predicted viscosity that was too low and relaxation times that were too fast, [7] compared to experimental data. A polymer modified version of ZG08 was also simulated, in which the composition was augmented by one molecule of polystyrene, which was selected to represent the more polar fraction of the styrene-butadiene-styrene (SBS) co-polymers that are used in practice to modify asphalts. [44] In this work we refer to that polymer modified system as ZG08p.
In a newer model asphalt system [10] representative of real AAA-1 asphalt, the asphaltene structure of Groenzin and Mullins was replaced with three types of asphaltene molecules proposed by Mullins[45] that are considered more consistent with experimental characterizations and quantum mechanics studies. Some side chain positions were changed to alleviate high energy "pentane effect" overlaps. [46] Five types of polar aromatics compounds were chosen that were consistent with geochemistry studies. [47,48,49,50,51] Two types of naphthene aromatics were taken from previously proposed asphalt models. [52,53] Saturates were chosen to be consistent with characterizations by Netzel and Rovani; [54] a notable change from our prior models is the use of branched and naphthenic compounds (squalane and hopane [55]), because the cited measurements [54] reported that 97 to 99 wt% of alkanes are branched or cyclic in real AAA-1 asphalt. In this work we refer to this system as LG14. It contains 72 molecules among the 12 molecule types. Full details and discussion of its composition are available elsewhere. [10] Time-dependent stress relaxation modulus data were taken from stress fluctuations calculated in the prior molecular dynamics simulations. [7,8,31] There they were calculated by in which V is volume, k B is the Boltzmann constant, T represents temperature, and σ st is a symmetric traceless stress tensor that was calculated using the molecular virial. [56] The stress tensor σ is averaged with its transpose and has the pressure subtracted from each diagonal component; the pressure equals one third of the trace of σ. Then summations over directions (x, y, z), indicated by indices u and v, average over all 9 components of the symmetric, traceless stress tensor. This Green-Kubo method [34] obtains the stress relaxation modulus from spontaneous fluctuations in stress under zero shear conditions, rather than by imposing a prescribed time-dependent strain.
The time integral of Eq. (20) was proposed [57,58] for use in estimating zero shear viscosity from shear fluctuations. Symmetrizing the stress tensor, removing the trace, and summing all components was found to yield a more accurate and quickly converged viscosity estimate, compared to averaging a single off-diagonal component of the stress tensor. [57,58] We assumed here that fluctuations of the stress tensor can be quantified by fluctuations of its symmetric, traceless form.
Equations (21) and (22) with G e = 0 indicate that numerical integration, enables discrete Fourier cosine and Fourier sine transforms to be applied to calculate storage and loss moduli. The time equals the number of equal time intervals, t = n∆t, and the corresponding angular frequencies are ω = 2π(n/N )/∆t. The maximum time available for G(t) equals N ∆t and is of order nanoseconds, which leads to minimum frequencies of order 10 −3 rad/ps.
Since Eq 24 provides a large amount of equally spaced numerical data, the transforms can be performed using a fast Fourier transform (FFT), which requires of order (N log 2 N ) operations. [59] The fastest Fourier transform in the west (FFTW) package has been exploited in this paper. Its detailed algorithm and method are described in the package documentation [60] and in the literature. [61] An FFT calculates a discrete Fourier transform [61], which by definition is periodic, [59,62] meaning that the FFT algorithm implicitly treats G(t) data as being periodic over the time domain considered. Thus G(t) data were augmented for all times t by G(−t) = −G(t) for storage modulus (Fourier sine transform, odd function) and by G(−t) = G(t) for loss modulus (Fourier cosine transform, even function).
Several facts must be taken into account while using FFT and interpreting data. The fewer MD results for longer time differences in G(t) lead to a poor signalto-noise ratio at low frequencies, though for a given set of time data, results from FFT for lower frequencies are more accurate than for higher frequencies. Vratsanos and Farris [63] suggest that the frequency range used in FFT should be ten times larger than the highest frequency result that can be used reliably. Considering the Nyquist-Shannon sampling theorem [64] and the MD time step of ∆t = 1 fs, the highest reliable frequency range is ω < 10 2 THz. Data beyond the first 10% (ω > 10 2 rad/ps) of the FFT output are not shown for this reason.
Next, to smooth the results and to avoid small fluctuations, a moving average filter was applied to all frequency domain results. The moving average was conducted over a lead and lag of 100 data points. These parameters were chosen such that high frequency noise has been removed from the FFT outputs. Whenever not enough points were available for the lead or lag direction, the extent of the lead and lag were both reduced to the available number of points. Applying a moving average on the frequency domain was found to be a better method that led to reasonable and smooth results, while applying a moving average to G(t) data led to complex modulus results that still contained harmonic noise.
where G i is the i-th stress modulus term and τ i is the i-th relaxation time, for a set of m discrete relaxations. The storage and loss moduli for the Maxwell model can be written as [65] For low frequencies, Eqs. (8) and (9) with slopes of +2 and +1 on log-log scales. For high frequencies, Eqs. (8) and (9) simplify to G ′ (ω) = m i=1 G i and G ′′ (ω) = (1/ω) m i=1 G i /τ i , with slopes of 0 and -1.

Results and Discussion
Stress relaxations for the LG14, ZG08, and ZG08p systems have been calcu- show the stress correlations fluctuate rapidly in a wide range of ±4 × 10 7 Pa after 400 ps. Ideally the G(t) reaches a plateau of zero for a viscoelastic liquid, but converging to this limit at all times is beyond the capability of the duration of our molecular dynamics simulations. Integrals of G(t) lead to viscosity and reach welldefined plateaus at high enough temperatures, as shown in prior work [9,8,66] and discussed in Fig. 11 below. Insets in Figs. 2 and 3 indicate results obtained by averaging G(t) over ±10 ps. This reduces noise significantly and illustrates the trends and fluctuations about an eventual zero limit.

Figures 2 and 3 illustrate that the stress relaxation modulus of the LG14
system is more than twice as large as that of the ZG08 and ZG08p systems at the same temperature (443 K). Figures B.1 and 3 show that incorporating polymer into the ZG08 model asphalt increased its stress relaxation modulus by a factor of up to 1.5. Furthermore, the higher the temperature for LG14, the smaller the relaxation modulus and the shorter the relaxation time.
To interpret zero-shear rheological behavior of model asphalts, Fourier sine Temperature suffixes and þ indicate the ZG08 and ZG08p systems; [9,7] others refer to the LG14 system. [8] and cosine transforms have been applied on the G(t) data to calculate storage and loss moduli. To examine thermorheological simplicity, complex modulus and phase angle results are shown together as a Black (or van Gurp-Palmen) plot in Fig. 4.
Because of large temperature differences, a vertical shift of |G * | has been applied by considering T = 533 K as reference. Over a to 3 × 10 8 Pa, phase angle decreases with increased modulus for each temperature, with a few exceptions. The results below 3 × 10 8 Pa show some coarse overlap among different temperatures. Results at higher complex moduli show a peak in phase angle (a loss peak) that does not overlap. In contrast, better overlap for this peak is found when applying a "reverse vertical shift", which corresponds to removing the volume and temperature factors in Eq. (20).
The origin of this peak is a subject of follow-up work. Because the curves at different temperatures do not overlap completely, the LG14 system does not follow thermorheological simplicity over the entire frequency range of the simulation. In other words, this model asphalt is thermorheologically complex and TTSP does not Gurp-Palman) plots for two asphalts described in the literature. da Silva et al. [16] modeled their rheology measurements of AC-1 and AC-2 bitumen systems at 313 K by using the Christensen-Anderson-Marasteanu (CAM) rheological model, [16,67] which relates complex modulus and phase angle to frequency as where ν and w are parameters. While bitumens AC-1 and AC-2 differ from the AAA-1 asphalt of SHRP, the intent of this comparison is to assess if the simulation results follow trends known to occur for real asphalts. Phase angle-complex modulus curves calculated over frequencies of 10 1 to 10 8 rad/s are shown in the figure, with a temperature and density shift applied to |G * | using a thermal expansion The G g parameter needed to be increased by a factor of 10 so the numerical results would match those shown in the figures of da Silva et al. [16]; other parameters for the CAM model were taken directly from their paper. The simulation results at frequencies below ∼1 rad/ps follow the qualitative trend of the CAM model fits. The CAM model does not incorporate the capability to display a loss peak, such as the peak found in the simulations. Figure 5 shows the averaged frequency-dependent moduli and phase angle of the LG14 system at 400 K. Over a frequency range of 7×10 −4 rad/ps to 1.5 rad/ps, storage modulus is larger than loss modulus and phase angle is passing through a minimum, which suggests that this is the plateau zone of the LG14 system. The fact that storage modulus is changing slowly and smoothly in this frequency range supports this idea as well. The presence of a plateau between 10 −2 and 10 0 rad/ps suggests that the model system acts entangled over simulation time scales of order ∼ 1/ω ∼ 100 ps at 400 K. Such entanglement is consistent with the decrease of rotation autocorrelation functions from 1 to 0.8 over ∼ 100 ps. [8] Figure B.6 illustrates moduli for both ZG08 and LG14 systems at 443 K.  . Black (or van Gurp-Palman) plot for the LG14 system for different temperatures. Lines indicate a qualitative comparison with CAM model fits reported [16] for AC-1 and AC-2 asphalt systems and shifted to 533 K.
Across the whole frequency range, the complex modulus magnitude is greater for the LG14 system than the ZG08, which consists of smaller molecules. These results show that using larger molecules leads to higher moduli and thus a stiffer model asphalt. The storage and loss moduli in the LG14 system are comparable, while in the ZG08 system the loss modulus accounts for almost all of the complex modulus (G * and G ′′ curves nearly overlap). We interpret these as simulation artifacts of a poor signal to noise ratio for this nearly viscous material. Analogous findings of a negative loss modulus for a nearly elastic material have been reported. [63] All moduli and phase angle of the LG14 system at 533 K are illustrated in Fig. 8. Over a frequency range of 9 × 10 −4 to 3 rad/ps, we find comparable relative increases in storage and loss moduli and a minimum in phase angle. Loss modulus is larger than storage modulus over this range, indicating more viscous than elastic behavior. A zero-shear viscosity can be estimated at this temperature. [8] To have a more thorough investigation of temperature effects on mechanical properties, the storage modulus, loss modulus, and |G * |/ω are separately illus- to 10. Furthermore, the ZG08 system at 443 K has a bigger viscous to elastic ratio than the LG14 system for frequencies below ∼10 rad/ps, as indicated by its larger phase angle. Its smaller molecules lead to a lower loss modulus (Fig. B.10) and much lower storage modulus (Fig. B.9). Storage moduli at all different temperatures converge in the transition zone (5 to 20 rad/ps) and then they remain more than 2 orders of magnitude larger than at lower frequencies. This trend occurs for loss moduli and |G * |/ω as well.  new results for all systems. |G * |/ω for the LG14 system at 400 and 443 K continues to increase as frequency decreases, which indicates that the limiting values of the zero-shear viscosities at low frequency have not yet been reached. |G * |/ω shows a slow rise with decreasing frequency at 533 K, with a semiquantitative approach toward the Green-Kubo value calculated previously. [8] The ZG08 and ZG08p systems shows good agreement at lower frequencies with the prior Green-Kubo results. [9] Also, the magnitude of complex moduli of the LG14 system at 400 and 443 K are larger by more than a factor of 10 at lower frequencies compared to the ZG08 and ZG08p systems, and they remain larger for higher frequencies as well. This again shows an effect of the larger molecules used for the LG14 system. The corresponding plot for the four component model of Hansen et al. showed convergence for frequencies below 10 9 Hz at 603 K. Experimental data for the zero-shear viscosity at temperatures within 10 degrees of 400 and 443 K for two asphalt systems described in the literature [69,32] are shown on the y-axis of Fig. 11. Differences between the two systems reinforce that "asphalt" refers to a class of materials for which their range of mechanical properties depends on the chemical constitution, molecular weight, and microstructure. The measured viscosities are higher than those calculated in the simulations at the least high frequencies, which is consistent with the rising trend with decreasing frequency that is seen in the simulation results. The frequency at which the rise converts to a plateau, as seen in the range −1.5 <∼ log ω < −0.2 for the 533 K system, cannot be predicted directly.
To show another effect of adding polymer on rheological properties of the model asphalt, Fig. 12 illustrates a Black (or van Gurp-Palman) plot for the ZG08 and ZG08p systems at 443 K. In the simulation, adding polymer to the asphalt increases the viscous contributions to the modulus at higher frequencies. Over 2 × 10 7 < |G * | < 5 × 10 8 Pa, the systems show significant overlap. At lower complex modulus (lower frequencies), the elastic contributrion is stronger in the system with the polymer.

Conclusions
Molecular simulations of model asphalts provide a method to understand how molecular motions contribute to viscoelastic response to applied strain or stress.
Storage and loss moduli were calculated using stress relaxation modulus G(t) data from previous equilibrium molecular dynamics simulations of model asphalt sys- tems ZG08, ZG08p, [9,7] and LG14. [ The following manuscript is prepared for submission to The Journal of Chemical Physics.

Abstract
Signal processing and its application in removing inherent noise in a correlation function has been investigated. Signal is a property that depends on at least one independent variable (e.g., time) and that carries information about the behavior or attributes of some phenomenon. Correlation function data for stress relaxation modulus obtained by equilibrium molecular dynamics simulations were

Introduction
Using data in the presence of large fluctuations is always challenging. The Usually noise is the unwanted part of input data or any unwanted modification of original data. Noise is an unavoidable part of our everyday lives and we have to deal with it everyday. Every measurement has some degree of uncertainty or noise in it. [1] Since 17th century, some techniques have been developed to remove these noise from our data. These techniques became the foundation of a science we call "signal processing" today. [2] Signal is a property that depends on at least one independent variable (e.g., time) and that carries information about the behavior or attributes of some phenomenon. A signal is not necessarily electrical and can be anything like an acoustic wave, thermometer output, pH-meter output, or even an electromagnetic wave. [3] To use already established techniques of signal processing, one can consider these data as signals. The act of doing some operations on a signal to produce another signal is known as signal processing. [3] Several methods can be applied to remove the noise from a given signal. Noise where A(ω) and B(ω) are Fourier cosine and sine transforms respectively, [12] A and L is the data acquisition time duration. The discovery of Fourier series and Fourier transform (FT) had a great influence on mathematics as well as our understanding about functions and integration theory. [13] Fourier transform is a powerful tool to calculate various problems including calculating complex integrals. [14] Fourier transform by definition can be applied on well-defined functions. In many cases we are interested in applying Fourier transform on a series of data that is not considered as a well-defined function, like equally spaced numerical data. This situation happens in many scientific and engineering applications like telecommunications, simulation problems, and data series. [15,16,17,12] These necessities led mathematicians to develop the discrete Fourier transform (DFT). [12] For discrete times, (13) and (14) lead to N is the total number of sampled data.
For large N , the calculation of Eqs. (15) and (16)  FFT codes calculate a discrete Fourier transform, [16] which by definition is periodic. [12,20] This means that the FFT algorithm implicitly treats data as being periodic over the time domain considered. This assumption causes the information to repeat beyond the given period. [19] The details of the theory and algorithm can be found elsewhere. [15,16,17,12] Fourier transform of a real function will give a complex function. Having just the magnitude of the transformed function in frequency cannot uniquely reconstruct the sinusoidal data in time domain. The other complementary information is called phase angle δ. Phase angle can be calculated using both an imaginary (β) and real part (α) of a complex number (ζ = α + iβ in general), [12] δ(ω) = tan −1 (β/α) The phase angle is usually represented in a 0 to 2π interval; this system is called "modulo 2π phase." [19] Understanding signals is not always straightforward. Several facts must be taken into account while using FFT and interpreting data. Sometimes errors are in the results, but it is claimed that more often the errors are simply misinterpretation of the results. [19] The rate at which data are being sampled determines how well those data to a poor signal-to-noise ratio at low frequencies, though for a given set of time samples, results from FFT for lower frequencies are more accurate than for higher frequencies. [19,21] Vratsanos and Farris [22] suggest that the frequency range used in FFT should be ten times larger than the highest frequency result that can be used reliably.
Fundamental frequency is the lowest frequency that can be reached by a given data acquisition duration, [12] This value is very important and will dictate the resolution in frequency domain, [19] Each frequency axis value is an integer multiple of this fundamental frequency and is referred to as a "harmonic." Since ∆f has a reciprocal relation with the length of sampling duration, increasing L can increase the resolution in the frequency domain.
The signal can be sampled in finite time duration. This interval can be selected by multiplying the signal by a function which is zero outside of the interested interval; this function is called a window function or simply a window. [18] Limiting the integral to a short interval or window is the same as multiplying the signal by a square pulse having a width equal to the interval of integration. [19] Limiting the length of the data acquisition time interval of the signal causes smearing of frequency domain spectral components to adjacent frequencies. This phenomenon is known as spectral leakage or simply leakage. [19,18]  Complex modulus is used to characterized asphalts by the SHRP. The complex modulus has two parts: storage modulus and loss modulus. Storage modulus, , is a measure of energy stored and recovered per cycle of sinusoidal deformation. Storage modulus can be calculated by Fourier transform from stress relaxation modulus (G(t)). The stress relaxation modulus can be related to stress tensor by V is volume, k B is the Boltzmann constant, T represents temperature, and σ st is a symmetric traceless stress tensor that was calculated using the molecular virial, [29] which incorporates the forces between molecule centers of mass and the kinetic energy of each molecule. The stress tensor σ is averaged with its transpose and has the pressure subtracted from each diagonal component; the pressure equals one third of the trace of σ. Then summations over directions (x, y, z), indicated by indices u and v, average over all 9 components of the symmetric, traceless stress tensor. The storage modulus is related to stress relaxation modulus G(t) by [30] G G e is the equilibrium tensile modulus and can be calculated by G e = lim t→∞ G(t). [30] For viscoelastic liquids, G e = 0. The angular frequency ω is defined using frequency as ω = 2πf . Loss modulus, G ′′ (ω), is a measure of energy lost per cycle of sinusoidal deformation and can be calculated by [30] G Phase angle (δ(ω)) corresponds to the difference in phase between stress and strain oscillations, and it defines as the ratio between viscous and elastic responses during the shearing process. [30] Using Eq. (17) the phase angle is defined as tan δ = G ′′ /G. [30] Calculating Eqs. (21) and (22) by definition is very tedious. Instead of using numerical integration to calculate G ′ (ω) and G ′′ (ω) one can use Fourier transform.
Equations (21) and (22) are in the same form of Eqs. (15) and (16) respectively; instead of using numerical integration to calculate G ′ (ω) and G ′′ (ω) one can use FFT to accelerate the calculations.

Methods
Time-dependent stress relaxation modulus data of LG14 were taken from stress fluctuations calculated in the prior constant volume molecular dynamics simulations. [9,31] This Green-Kubo method [32] obtains the stress relaxation modulus from spontaneous fluctuations in stress under zero shear conditions, rather than by imposing a prescribed time-dependent strain.
Numerical simulations have sampled continuous G(t) data to obtain discrete results; therefore the digital signal processing must be applied here instead of analog signal processing to interpret the results. G(t) data are the input correlation function that is used as our input signal here. Equations (21) and (22) with G e = 0 indicate that numerical integration, Since G(t) does not goes to zero in longer times, as for liquid viscoelastic liquids, a window function has been applied to G(t) data. The following window function has been applied χ(n∆t) = 1 n∆t ≤ a ps exp[−ζ(n∆t − a)/(b − a)] n∆t > a ps (25) ζ is decaying factor and considered to be 10. a is the start point of damping and considered to be 1 ps. This a value has been chosen to have the least effect on data for which we have more confidence about them. b is the sampling duration and it is 3500 ps. ζ is chosen in this way to have the least damping effect while forcing the window function to go to zero at its end.
Applying this window function leads to some suppression of frequency response magnitude. To restore the magnitude of frequency response, the results have to be divided by coherent gain. [33] Coherent gain of a given window function is calculated by [34] One of the most common filters to remove random noise is the moving average filter. The moving average is the fastest digital filter available. [17] Moving average is the optimal filter for removing random noise while retaining a sharp step response. As the number of points involved in the moving average increases, the noise becomes lower, and the edges of sharp step response becomes smeared. The moving average filter provides the least noise possible while retaining a given edge sharpness and it is optimal for removing random noise while retaining a sharp step response. It can reduce the noise by square-root of number of points used in the averaging. [17] To smooth G(t) data, a moving average filter was applied to G(t) data beyond the first 1ps. The moving average was conducted over a lead and lag of 0.1ps.
Whenever there were not enough points available for the lead or lag direction, the extent of the lead and lag were both reduced to the available number of points.
The fastest Fourier transform in the west (FFTW) package has been exploited in this paper to accelerate the calculations. It is based on FFT to reduce the computational cost. The detailed algorithm and method are described in the package documentation [35] and in the literature. [16] Using FFTW, the G(t) data have been used to calculate G ′ (ω) and G ′′ (ω).
To calculate storage modulus, G(t) was extended as an odd function for negative time span to use Fourier sine transform. An even function generalization has been applied to calculate loss modulus using Fourier cosine transform. Using these two moduli and Eqs. (17) and (23)

Results
Using Eq. (20), stress relaxation modulus has been calculated for LG14 system. The stress relaxation tensor has been averaged over all their elements to have one scalar stress relaxation at each time. The first nanosecond of stress relaxations for the LG14 system has been plotted in Figs. 1 and 2 in the previous paper. [11] The fluctuations at longer times were in the range of ± 8 × 10 8 Pa after 400 ps.
These high amplitude fluctuations in the G(t) data after 400 ps will introduce spectral leakage in FFT results. To reduce the increasing random fluctuations at longer relaxation times, which is caused by a decreasing number of available time separations in the signal averaging of Eq. (20), a moving average of ± 0.1 ps has been applied to all relaxation times beyond 1 ps. This moving average can reduce the spectral leakage in FFT results. The stress relaxation data beyond 1 ps for LG14 system after applying moving average on it is shown in Fig. 13. Since the data within first picosecond of stress relaxation data were the most reliable part, no operation has been applied to this part. Stress relaxation in Fig. 13 does not go to zero at longer time. This cause spectral leakage in frequency domain. Furthermore, from the physics of viscoelastic liquid it is expected that G(t) goes to zero at longer times. To make this decay happen, the window function of Eq. (25) has been applied here. This window function has the benefit of not changing the first picosecond of stress relaxation data while damping the data beyond it in a way that has the least change on data which have more certainties. The data at longer times, which have less certainties, will change more. Figure 14 shows the effect of this window function (Eq. (25)) on G(t) signal.
The applied window function did not change the first 10 ps of G(t) data that much. Between 10 ps and 100 ps, the noise damping effect of the applied window starts to be noticeable. Beyond 100 ps, the noise have been damped a lot and the transformed signal reached 0 Pa after 1000 ps. This window function forced the G(t) to meet the 0 Pa limit of stress relaxation modulus of viscoelastic liquids. Applying the window function on G(t) data along with applying moving average on moduli results and removing insignificant components led to additional improvements in the results. Figure 18 shows the results of these techniques. In this case, the severe fluctuations in storage modulus as seen in Fig. 17 in frequency range of 3 to 30 rad/ps have disappeared.
The best improvement was after applying the window function and moving average on G(t) data plus removing insignificant components from frequency domain data and applying moving average on all moduli while sampling at highest rate. Figure 19 shows the result that was obtained after applying all these techniques.
Phase angle has no high amplitude fluctuation and moduli are smooth.
To show the effect of sampling rate on the frequency domain results, the same techniques that were applied in Fig. 19 have been applied to G(t) and all moduli data with a lower sampling rate of 1 sample per 5 fs, compared to 1 sample per 1 fs in previous case. The results are shown in Fig. 20. Using a lower sampling rate induces some low amplitude fluctuations especially in storage modulus and consequently phase angle and particularly at frequencies above 1 rad/ps. Figure 19 is clear enough to investigate asphalt rheological properties at 533 K.
The fact that storage modulus is less than loss modulus for frequencies below 200 rad/ps shows that asphalt is more viscous than elastic at this temperature.
For frequencies above 200 rad/ps, the elastic components are more dominant. The less noisy data should also enable investigating time-temperature superposition more clearly.

Conclusions
Noise removal can improve the signal to noise ratio and lead to more clear  Figure 15. Storage modulus and loss modulus that correspond to raw G(t) data for LG14 at 533 K.

MANUSCRIPT 3 Equation of State Modeling of Bitumen Thermodynamic Properties
The following manuscript is prepared for submission to the journal of Industrial & Engineering Chemistry Research. was investigated. In a specific temperature, the chemical potential of pure squalane below its melting point was lower than its chemical potential in the multicomponent system which suggests that squalane will precipitate below its melting point. The temperature in which an asphalt model precipitate is different for each models.
This confirms that the chemistry of asphalt can effect its mechanical properties.

Introduction
Bitumen is a widely seen material nowadays but it is not a new discovery. It has been used since ancient era. [1,2] It has been used in shipbuilding in Sumeria around 6000 B.C. [2] and as cement or adhesive material in sculptures by Sumerians and Babylonians. [1] The term "asphalt" is traced back to Babylonians and means "secure" or "firm", which implies its application as an adhesive material. [1] The first record of using asphalt in pavement was in Babylon around 625 B.C. [3] This material has been used without a thorough understanding of its composition and properties. Understanding chemistry [4,5,6,7], viscoelastic properties, [8,9,10,11,12,13,4,14,15,16,17,18] and thermodynamics [11,19,20,21,22,23,24] of phase behavior asphalts and asphaltenes are routes toward improving its performance and durability. Recently the chemistry of asphalt and its impact on its properties have been investigated by several authors. [5,6,4,7] The aim of this paper is to provide required parameters for corresponding states [25] based equations of state (EOS) to enable predicting thermodynamic properties of bitumen. As one example, the density and thermal expansion will be calculated without molecular calculations or experiments. A second example will demonstrate the driving force for precipitation of wax from the bitumen, a problem known from experiments. Waxes are undesirable high molecular weight alkanes in petroleum and petrochemical industries.
The quality and performance of asphalt depends on crude oil source, and therefore depends on chemical composition of asphalt. [11] Naphthenic based crude oils are known to yield high quality asphalt while paraffinic crude oils are known to show poor performance on roads. [26] The major point of this paper is to demonstrate a multiscale method that can account for the effects of chemistry while modeling asphalt properties.
The first step in molecular simulation of any system is to define the molecules in the system. Asphalt is a mixture of 10 5 to 10 6 different molecule types. [27] Doing molecular studies on a system containing this huge number of molecule types is laborious if not impossible. Instead, we use a small set of components that represents the asphalt properties. Several representative models have been introduced. [28,29,30,31,32] Several researchers using advanced EOS like SAFT [33,34] have been carried out on asphaltene and asphalt recently. [21,19,22,20,23,35,24] While these methods are supposed to be more accurate, they need lots of experimental data to estimate unknown parameters. Ting et al. [19] have predicted stability boundaries of asphaltene in reservoir conditions. Gonzalez et al. [20,22]  The Strategic Highway Research Program (SHRP) introduced and named different asphalt binders for research studies. The temperature dependent complex modulus and viscosity were used by SHRP to characterize different binders. [40,11] Although magnitude and temperature dependence of complex modulus and viscosity describe many asphalt mechanical properties, [40] asphalts with comparable complex modulus show different performance in pavement applications because of independent behaviors of other mechanical properties. [2,11,41] Our group proposed several combinations of compounds that represent asphalt as a system rather than asphaltenes as a solubility class. Initial work employed ternary systems. [28] A system with a larger number of more polar molecules [29] was proposed next to represent atomic composition of AAA-1 asphalts in SHRP classification. Here, it has been refered as ZG08. Its molecules were later found to be too small. [42] A new set of yet larger molecules was proposed recently [43,30] to represent SHRP asphalts AAA-1, AAK-1, and AAM-1 (labeled here as LG14A, LG14K, LG14M) as composition variations of the same compounds. Details that motivate these model representations of AAA-1, AAK-1, and AAM-1 can be found elsewhere. [29,30] The critical properties of several hundreds of general compounds are available in the literature. [44,45] But the critical properties of interest are not always available in the literature for all the molecules we need for modeling asphalt. Although experimental methods are available for measuring critical properties of materials, [46,47] measuring these properties are sometimes costly [48] or require special techniques for unstable materials. [49,50,51,52,53] Sometimes experi-mental methods are not even possible because some molecules in the models are hypothetical and are nontrivial to produce as pure compounds.
An alternative way of finding unknown critical properties is group contribution (GC) methods, which have been shown to be a great alternative. [54,55,56] The underlying theory behind GC methods is that the properties of a compound can be calculated from the properties of its fragments (groups) [57]. Full discussion about different methods and their parameters can be found elsewhere. [44,45,58,59,60,61,62,54,55,56] A relationship between critical temperature and normal boiling point has been studied for more than 100 years and was first reported by Guldberg [63] in 1890.

This relation has been modified several times and almost all GC methods follow
Guldberg's relation, [63] T c = 1.5 T b (27) where T c and T b are in Kelvin.
A GC method for critical properties has been introduced by Riedel, [57] and since then several modifications [64,65,66,67,54,55,56] have been introduced to improve the accuracy and reliability of prior approaches. Joback and Reid [64] Tables 3 and 4 respectively.
The details and discussion of these model asphalts are available elsewhere. [42,30]

Methods
The relationship between temperature (T ), pressure (P ), and molar volume a is the energy parameter and is a function of mole fractions and critical properties, [69] a = n i n j x i x j a ij (29) where a ij is [69] a ij = (1 − k ij )(a i a j ) 1/2 (30) k ij is a binary interaction parameter between species i and j. Here we set all k ij = 0 for simplicity. a i and a j are pure component parameters and can be calculated by critical properties [69] a i = 0.45724 R is the universal gas constant, T c is critical temperature, P c is critical pressure, T r = T /T c , and ω is acentric factor. [69] α i is calculated from reduced temperature by [69] κ i is a function of acentric factor only, [69] b is the mixture volume parameter, [69,70] b i is the pure component volume parameter which can be estimated by [69] b i = 0.07780 Since most of the compounds in the models are hypothetical compounds, their critical properties are unknown. To estimate unknown critical properties, the method of Nannoolal et al. [54,55,56] has been used. This method was chosen because of its generality, reliability, and accuracy. This method is a three-order method: [54,55,56] 1. The first-order is just like any other basic GC method. The first order contribution is calculated by summing the occurrence frequency of any group multiplied by its contribution value.
2. The second-order contribution accounts for contributions from some special functional groups as well as isomers.
3. The third-order accounts for interactions between different groups.
We have implemented the method of Nannoolal et al. [54,55,56]  A corresponding states related parameter that is required by most cubic equations of state is the acentric factor, ω, which is a measure of the non-sphericity of the molecules. [71,72,25] While this parameter was not provided directly in papers by Nannoolal et al., they have provided an expression for saturation pressure (P * ) that we have exploited to calculate acentric factor. Their expression for saturation pressure in atm is [56] log P * = (4.1012 + dB) where T b is the normal boiling point temperature, T rb = T /T b , and GI is a third-order group interaction contribution, which is [54] (38) m is the total number of interacting groups and n is the number of atoms in molecule excluding hydrogen atoms. The notation C i−j indicates the contribution from interactions between group i and group j; it is not a subtraction. We calculated acentric factor using its definition, [71,25] ω ≡ −1.0 − log P * r | Tr=0.7 (39) where P * r = P * /P c . The temperature 0.7T c required for acentric factor can be written using the expression from Nannoolal et al. for critical temperature [55] as Substituting T rb,ω in Eq. (36) and using the definition of acentric factor with the expression from Nannoolal et al. for critical pressure leads to M is the molecular weight in g/mol. The C i and GI terms in equations 37 and 41 are different, and all C i and C i−j have separate values that can be found else- where. [54,55,56] Unfortunately, most cubic equations of state do not predict liquid volume accurately. [73,74,75] Several methods have been developed [76,77,78,79,80,81,73,82] that improve liquid density results without altering phase equilibria.
In this paper, the COSTALD [81] method has been used for calculating sat- V • is the characteristic volume for each pure component. If no experimental value available for the characteristic volume, the following expression can be used, [81] Here ∆ H f usion is in units of kJ/kgmol. N i and C i are the occurrence frequency of group i and the group contribution of group i respectively. C i parameter can be found in the Joback and Reid paper. [64] The solid phase free energy for a pure component has been estimated away from T m by assuming that ∆H and ∆S are independent of temperature over a moderate range. This leads to The subscripts S and L indicate pure solid and liquid phases, with liquid phase properties obtained from the equation of state. The pure component chemical potential equalsḠ i (intensive free energy). Squalane was chosen from the compounds in LG14 as a representative component for wax precipitation. Its chemical potential in the multicomponent bitumen has been calculated using the usual EOS approach.

Results
Published data exist for some of components involved in this study. [83,84,85,53] The estimated critical properties and acentric factor of them have been listed in modeling ZG08 and LG14 model asphalts.
The estimated critical properties and acentric factors of components in ZG08 model asphalt are summarized in Table 3 and of LG14 model asphalts in Table   4.  Table 4. Molecule types, number of molecules, and thermodynamic properties for compounds in LG14. [30] The units for properties are the same as Table 3.
bles 3 and 4 and with all k ij set to zero in Aspen HYSYS. Figure 21 shows the temperature-dependent density predictions for both ZG08 and LG14 asphalt models. Experimental density data and thermal expansion coefficient [86] for SHRP AAA-1 asphalt binder are shown in this figure. The density of AAA-1 asphalt at different temperatures were previously [29,30] Table 5 along with the experimental and MD results for AAA-1 asphalt.
For ZG08 model asphalt, COSTALD correlation can provide the most accurate thermal expansion prediction. The PR EOS prediction for thermal expansion is better than MD prediction but it is not as good as COSTALD prediction.
COSTALD and MD predictions for LG14 model asphalts are equally good while LG14-PR LG14-COSTALD LG14-MD Experimental ZG08-PR ZG08-COSTALD ZG08-MD Figure 21. Density vs. temperature for AAA-1 asphalt binder determined for ZG08 and LG14 models using MD (solid line), PR EOS (dashed), and COSTALD (dotted) compared to experimental data (dot-dashed).  Table 5. Thermal expansion of AAA-1 asphalt binder and its absolute relative error with respect to experimental data.      Figure 24. Comparison between ∆H/RT from MD and EOS calculations ZG08 model asphalt is 814 K and its critical pressure is 31 atm. The mixture critical temperature for AAA-1 is 955 K and the mixture critical pressure is 11 atm. Figure 25 and Tables 3 and 4 can be used to conclude that ZG08 molecules are more volatile in asphalt processing condition. The ZG08 model asphalt bubble point at 1 atm is 500 K, and bitumen processing usually extracts volatile species at 800 K. [90] This is a further indication that the molecules in the ZG08 model asphalt are too small and volatile and are not suitable for being used specially in distillation column conditions and asphalt processing conditions.

Conclusions
The critical properties and acentric factors of all components involved in both ZG08 and LG14 model asphalts have been estimated using the methods of Nannoolal et al. [54,55,56] The densities and the thermal expansions of both model asphalts were calculated in Aspen HYSYS 13.2 using these estimated properties and the PR EOS.
LG14 model asphalt predicted densities and thermal expansion coefficient more accurate with respect to experimental data. The calculated densities from MD is better than results from PR EOS corrected by COSTALD LG14-BPP LG14-DPP Figure 25. T-P envelope of modeled AAA-1. BPP stands for bubble point pressure and DPP stands for dew point pressure.
correlation. But the calculation cost is which smaller using the EOS. The PR EOS is not performing very well in predicting densities, which is common in cubic equations of state. On the other hand, PR EOS without volume translation is predicting thermal expansion of asphalts much better than any other methods.
Comparing the squalane chemical potential in the multicomponent system to that of pure liquid or solid squalane shows that the chemistry of an asphalt binder can affect the temperature at which waxes can precipitate.  where N i is the frequency of the group i in the molecule, C i is the group contribution of group i, and n is the number of atoms in molecule excluding hydrogen atoms. GI is the group interaction contribution (third-order), which is [1]

List of References
m is the total number of interacting groups. The notation C i−j indicates the contribution from interactions between group i and group j; it is not a subtraction.
While T b does not appear directly in the EOS, its value is required for critical temperature calculations.
The relation between critical temperature and normal boiling point is known for more than 100 years and was first reported by Guldberg [2] in 1890. This relation has been modified several times and almost all GC methods follow Guldberg's relation, [2] T c = 1. The critical pressure (P c ) is another parameter which is used by equations of state based on "principle of corresponding states." [4] This value can be calculated by [3] P c = M −0.14041 which P c is in kP a and M is the molecular weight in g/mol.
critical volume can be calculated based on Nannoolal method as [3] V c = i N i C i + GI n −0.2266 + 86.1539 (A.6) The dimension for critical volume is 10 −6 m 3 mol −1 .

APPENDIX B
Improved results for chapter 1 This appendix is showing the updated results of Chapter 1 after applying some signal processing techniques to improve the frequency domain results.