Molecular dynamics (MD) simulation is a powerful tool to reveal the microscopic characteristics of supercritical transitions. However, the accuracy of MD depends strongly on the potential model that describes the interaction forces between atoms. In this study, four commonly used potential models for long-chain n-alkanes in MD simulations are evaluated, and a hybrid model is introduced. The vaporization and phase-transition characteristics of n-alkane blended fuels with different mole fractions are then explored under a wide variety of ambient conditions by using the hybrid model. Compared to the commonly used potentials, the hybrid model shows higher accuracy for predicting the thermodynamic and transport properties. In subcritical environments, vaporization belongs to typical two-phase evaporation with a sharp gas–liquid interface. The preferential evaporation of the light-end component is obvious, and the evaporation rate of the heavy-end component is maximized after the light-end component is consumed. Under supercritical conditions, the interface dissolves rapidly, the evaporation rates for both the light- and heavy-end components increase simultaneously, and both components coexist throughout the evaporation process. Based on the maximum potential energy and evaporation rate, a new criterion for the supercritical transition is proposed. The dimensionless transition time, which reflects the proportion of the sub/supercritical stage within the lifetime, is nearly independent of the ambient temperature and fuel composition; instead, it mainly depends on the ambient pressure. Finally, an empirical formula is obtained by curve-fitting to describe the variation in the dimensionless transition time with ambient pressure.

Fuel injection is a critical factor affecting the ignition and combustion characteristics of power generation equipment such as rockets, diesel engines, and gas turbines. To achieve more efficient and cleaner combustion, higher injection and boosting pressures are needed. These pressures may be higher than the critical points of the fuels, resulting in a complicated situation in which both subcritical and trans/supercritical injections exist simultaneously.1 A deeper understanding of these phenomena is crucial for developing the corresponding models to guide engine performance optimization.

Substantial efforts have been made to better understand the supercritical vaporization characteristics and transition mechanisms behind trans/supercritical injections. Existing research has mainly focused on single droplets and jets/sprays using optical measurements,2–7 computational fluid dynamics (CFD),8,9 and theoretical analyses.10–12 Morin et al.13 studied the effects of ambient temperature and pressure on the evaporation rates of n-heptane (n-C7H16) and n-decane in consideration of the critical conditions of the liquid fuel in a high-pressure and high-temperature gasification facility. Nomura et al.14 investigated the vaporization characteristics of a single n-hexadecane droplet in supercritical environments under microgravity conditions. Harstad and Bellan15 proposed a basic theory for modeling single-droplet vaporization under sub/supercritical conditions, and Zhu and Aggarwal16 highlighted the effects of the equation of state (EOS) on single-droplet evaporation and the minimum pressure required for supercritical transition.

The injection of cryogenic liquid propellants under various supercritical conditions has been extensively studied by experiments,3,17,18 demonstrating that a supercritical jet is close in appearance to a turbulent gas jet with significant diffuse behaviors. Recent studies have examined the supercritical injection of diesel-like fuels under diesel engine-relevant conditions. Manin et al.19 experimentally studied the atomization and mixing processes of Spray A [i.e., n-dodecane (n-C12H26)] under sub/supercritical conditions in a constant-volume chamber by using long-distance microscopy and diffused back-illumination. Crua et al.4 investigated the transition of microscopic droplets at the end of diesel injections from the subcritical condition to the diffuse mixing regime as functions of the critical temperature and pressure of different fuels.

Trans/supercritical injection modeling20–24 is an effective method that can complement experimental studies. To accurately predict the complex behaviors, a real-fluid EOS is needed. However, a single universal model that can well reproduce the sub/supercritical injection regime is difficult to construct due to the liquid–gas interface issue (i.e., the interface disappears during the transition from the subcritical state to the supercritical state) and the spurious pressure oscillations caused by the non-linearity of the thermodynamic properties crossing the Widom line. Therefore, separate modeling approaches are usually selected for the subcritical and trans/supercritical injection scenarios, and diffuse interface models have been successfully applied to trans/supercritical injection. The characteristics of nitrogen (N2) jets under different N2 supercritical conditions were reported by Zhang et al.25 based on large-eddy simulations. For multi-component systems, the vapor–liquid equilibrium (VLE) must also be evaluated to test the phase stability and determine the local phase fractions for evaporation modeling in CFD simulations.24,26,27

Although the above macroscopic experimental and numerical studies provide valuable insights into the trans/supercritical transition, the transition mechanism is not fully understood due to the limitations of the measurements and some assumptions adopted in the numerical simulations. Alternatively, molecular dynamics (MD) simulation is a powerful tool to reveal the microscale supercritical vaporization characteristics and transition mechanism because MD does not introduce any assumptions about the investigated physical process,28–31 unlike continuum-based numerical methods. However, the accuracy of MD depends strongly on the potential model that describes the interaction forces between atoms. Commonly used potential models for hydrocarbons include all-atom (AA) and united-atom (UA) models. AA models treat each atom as one site, while UA models take each methyl (CH3) and methylene (CH2) group as two different sites. Therefore, UA models are widely adopted in MD simulations of hydrocarbons with large numbers of carbon atoms due to the much higher computational efficiency. However, comparative studies of different UA models in terms of their accuracy in predicting hydrocarbon thermodynamic and transport properties in MD simulations are limited.

Recently, MD simulations have attracted considerable attention for exploring hydrocarbon fuel vaporization and phase-change characteristics. Mo and Qiao32 studied the supercritical transition characteristics of three pure n-alkanes in nitrogen environments using the SKS potential model.33 Xiao et al.34 investigated the sub/supercritical vaporization characteristics of pure n-C12H26 by using the Optimized Potentials for Liquid Simulations (OPLS) UA model.35 Wang et al.36 studied the evaporation and phase transition modes of hydrocarbon droplets in supercritical environments by using the Transferable Potentials for Phase Equilibria (TraPPE) UA model.

Considering the multi-component nature of practical fuels, Zhang et al.37 explored the gas–liquid interfacial characteristics of binary n-alkane fuels and tested the VLE assumption adopted in the continuum-based vaporization models based on the OPLS UA model. Gong et al.38 investigated the transition characteristics of a three-component hydrocarbon fuel by using the TraPPE UA model39 and developed a transition criterion based on the quantitative Voronoi tessellation. Wang et al.40 studied the influence of multi-component droplets on the relationship between the supercritical transition and evaporation by using the OPLS model. By using the SKS model, Chakraborty and Qiao41 investigated the multi-component characteristics of fuel and ambient gas during supercritical evaporation. However, the comprehensive effects of ambient temperature, ambient pressure, and fuel composition on the supercritical vaporization and transition characteristics under a wide range of ambient conditions have not been well studied.

Motivated by the above considerations, we examined the ability of UA potential models commonly used in MD simulations to reproduce the thermodynamic and transport properties of n-alkanes. We then introduced a hybrid model to better predict the above physical properties related to fuel vaporization and evaluated the vaporization and phase transition characteristics of n-alkane blended fuels with different mole fractions under a wide range of ambient temperature (T) and pressure (P) conditions.

We considered four commonly used UA models for hydrocarbons: OPLS,35 TraPPE,39 NERD,42 and SKS.33 These models have the following similar form but with different model coefficients:

Utotal=Ui,j+Ub+Uθ+Uψ,
(1)

where Ui,j, Ub, Uθ, and Uψ are the non-bonded, bond stretching, bond angle bending, and bond angle torsion potentials, respectively. The non-bonded potential Ui,j of atoms i and j can be described by the Lennard-Jones 12–6 potential,

Ui,j=4ϵi,jσi,jri,j12σi,jri,j6σi,jrc12σi,jrc6,
(2)

where σi,j and ϵi,j are the size and energy parameters, respectively, and rc is the cut-off distance beyond which the potentials between different sites are neglected. To determine the Lennard-Jones parameters for the interactions between different kinds of atoms, we adopted the Lorentz–Berthelot mixing rule,

ϵI,j=ϵiϵj,σi,j=12σi+σj.
(3)

The following harmonic potential was adopted to describe the bond stretching between two bonded atoms:

Ub=12kbririeq,
(4)

where kb is the stretching coefficient and rieq is the equilibrium bond distance.

Bond angle bending exists among three neighboring atoms, which can be described by another harmonic potential as follows:

Uθ=12kθθiθieq,
(5)

where kθ is the bending coefficient and θieq is the equilibrium angle.

Bond torsion, which takes place among four neighboring atoms in the same molecule, describes the vibration of the dihedral formed by these atoms. The following potential is employed:

Uψ=k11+cosψ/2+k21+2cosψ/2+k31+3cosψ/2,
(6)

where k1, k2, and k3 are the model coefficients.

The Rivera and Alejandre model43 was adopted to describe the interatomic potential of the N2 molecule. The non-bonded interactions can also be characterized by the Lennard-Jones 12–6 potential shown in Eq. (2) with parameters of σN = 3.31 Å and ɛN = 0.072 kcal/mol. The distance between two bonded N2 atoms was fixed at 0.108 97 nm.

Based on our previous simulation results, among the four UA models, the SKS model has the highest accuracy in predicting the critical point and thermodynamic properties of n-C7H16. However, the SKS model still fails to quantitatively reproduce the physical properties of longer-chain alkanes (e.g., n-C12H26). Replacing the three-parameter bound torsion model shown in Eq. (6) with the following five-parameter model taken from Refs. 44 and 45 can considerably improve the results predicted by the SKS model,

Uψ=n=15Ancosn1ψ,
(7)

where An are model coefficients. This hybrid model is hereafter referred to as SKS-5; model validation and comparison are found in Sec. III A. In addition, detailed information on the model coefficients for all the UA models used in this study can be found in Table I.

TABLE I.

Coefficients for different UA models.

Potential parametersOPLSSKSTraPPENERDSKS-5
ϵCH3kcal/mol 0.175 0.2265 0.195 0.207 0.2265 
ϵCH2kcal/mol 0.118 0.0934 0.091 0.091 0.0934 
σCH3nm 0.3905 0.393 0.375 0.391 0.393 
σCH2nm 0.3905 0.393 0.395 0.393 0.393 
kbkcal/molnm2 52012 19 177 
rieqnm 0.1526 0.154 
kθkcal/molrad2 124.15 124.15 
θieq° 112 114 
kn/Ankcal/mol k1=1.411k2=0.271k3=3.144 A1=4.825A2=0.162A3=6.224A4=0.325A5=0.503 
Potential parametersOPLSSKSTraPPENERDSKS-5
ϵCH3kcal/mol 0.175 0.2265 0.195 0.207 0.2265 
ϵCH2kcal/mol 0.118 0.0934 0.091 0.091 0.0934 
σCH3nm 0.3905 0.393 0.375 0.391 0.393 
σCH2nm 0.3905 0.393 0.395 0.393 0.393 
kbkcal/molnm2 52012 19 177 
rieqnm 0.1526 0.154 
kθkcal/molrad2 124.15 124.15 
θieq° 112 114 
kn/Ankcal/mol k1=1.411k2=0.271k3=3.144 A1=4.825A2=0.162A3=6.224A4=0.325A5=0.503 

The simulations were carried out by using the Larger-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) code. The Velocity–Verlet algorithm was employed to integrate the equations of atomic motion with a time step of 2 fs in all cases. The initial atom velocity was assigned randomly based on the Gaussian distribution. To eliminate the boundary effect caused by the scaling limitation of the simulation system, periodic boundary conditions were adopted in the three directions of the simulation box. More specific information on the physical properties and vaporization calculations can be found below.

1. Physical properties

A cuboid box was used for the gas–liquid–gas configuration to calculate the saturation thermodynamic properties, as shown in Fig. 1(a). The liquid phase was placed at the center of the box, while the two ends were the gas regions when equilibrium was reached. The simulation box was equally divided into Ns slices along the x direction (i.e., perpendicular to the interface), and the spatial distribution of thermal properties such as temperature, pressure, and density could be obtained through the microscopic properties of atomic number, velocity, and mass. The temperature and pressure in the kth slice could be determined by the following equations:

Tk=13NkkBi=1Nkmivi2,Pk=NkkBTVk+i=1NkriFi3Vk,
(8)

where Nk is the atom number in the kth slice; Vk is the corresponding slice volume calculated as Vk=LxLyLzNs, where Lx, Ly, and Lz are the box lengths in the x, y, and z directions, respectively; and mi, vi, ri, and Fi are the mass, velocity, radius vector, and force of the ith atom, respectively.

FIG. 1.

Simulation configurations for calculating the thermodynamic and transport properties: (a) gas–liquid–gas configuration for calculating saturated thermodynamic properties and (b) pure liquid-phase configuration for determining the self-diffusion coefficient.

FIG. 1.

Simulation configurations for calculating the thermodynamic and transport properties: (a) gas–liquid–gas configuration for calculating saturated thermodynamic properties and (b) pure liquid-phase configuration for determining the self-diffusion coefficient.

Close modal

The cubic box configuration shown in Fig. 1(b) was employed to determine the self-diffusion coefficient. The molecular motion can be approximated as a random walk when the time is sufficiently long, and the mean square displacement (MSD) is linear with time. The self-diffusion coefficient D can be calculated from the relationship between MSD and time according to the following Einstein formula:

D=limtrtr026t.
(9)

The canonical (NVT) ensemble was employed to keep the temperature of the system constant throughout the simulation process when calculating the saturation thermodynamic properties using the configuration shown in Fig. 1(a). However, the isothermal–isobaric (NPT) ensemble was run for 1 ns to stabilize the pressure of the entire system to the normal pressure, and the NVT ensemble was then run for 2 ns to determine the self-diffusion coefficient using the configuration shown in Fig. 1(b). Other simulation settings can be found in Table II. When the four indices of temperature, pressure, total energy, and potential energy of the entire system fluctuated around a fixed value, the system was considered to have reached the equilibrium state, and the preceding 1 ns was selected as the statistical time to calculate the corresponding physical properties.

TABLE II.

Mesh information for the single-droplet transition case.

ConfigurationBox length (nm3)Number of moleculesTemperature range (K)
Thermodynamic 
(n-C7H1630 ×8× 8 2048 293–533 
Thermodynamic (n-C12H2640 ×7.7× 7.7 1536 450–650 
Transport (n-C7H165 ×5× 5 500 273–368 
Transport (n-C12H266 ×6× 6 576 308–364 
ConfigurationBox length (nm3)Number of moleculesTemperature range (K)
Thermodynamic 
(n-C7H1630 ×8× 8 2048 293–533 
Thermodynamic (n-C12H2640 ×7.7× 7.7 1536 450–650 
Transport (n-C7H165 ×5× 5 500 273–368 
Transport (n-C12H266 ×6× 6 576 308–364 

2. Vaporization

To investigate the fuel vaporization characteristics under a wide range of ambient conditions, the liquid phase composed of n-C7H16 and n-C12H26 in different mole ratios was placed at the center of a three-dimensional cuboid box with dimensions of 265 × 100 × 100 nm3, and N2 molecules were located on the two sides of the box, as shown in Fig. 2. The spatial scale used in this study was capable of eliminating the scaling effect. Detailed information on the effect of the spatial scale on the computational results can be found in our previous work.37 

FIG. 2.

Configuration of the simulation system for fuel vaporization in a nitrogen environment.

FIG. 2.

Configuration of the simulation system for fuel vaporization in a nitrogen environment.

Close modal

Initially, the liquid phase was placed in the middle of the simulation box and equilibrated by using the NVT ensemble with a constant liquid temperature of 363 K (standard fuel injection temperature for Spray A19). Meanwhile, the N2 molecule was also equilibrated by using the NVT ensemble at the target temperature and pressure. Once the liquid and gas phases were in equilibrium, the two phases were combined, and the micro-canonical (NVE) ensemble was employed for the entire domain to reproduce the natural process of heat and mass transfer. For the heating/deletion regions located at the two ends along the length direction, additional Langevin thermostats were adopted to accurately maintain the target ambient temperature. Meanwhile, to reproduce the true evaporation process, the fuel molecules were deleted when they entered the heating regions (i.e., fuel vaporizes into an infinite space); thus, the subsequent evaporation process was not affected by the vapor phase of the fuel.

The supercritical vaporization and transition characteristics were explored under a wide range of ambient conditions. The ambient temperatures ranged from 800 to 1200 K, which are higher than the critical temperatures of n-C7H16 and n-C12H26. The ambient pressures ranged from 10 to 120 bars, covering the critical pressures of both components. The total calculation time for the vaporization case was set as 100 ns for all conditions.

This section describes the evaluation of three physical properties (the VLE curve, saturated vapor pressure, and self-diffusion coefficient) of n-C7H16 and n-C12H26 to test the performance of the five potential models.

The saturated density distribution curves of n-C7H16 at different temperatures using the configuration shown in Fig. 1(a) are shown in Fig. 3. When the entire system is equilibrated, the liquid, gas, and gas–liquid interface can be identified based on the density distribution. The five UA models had good performance below 473 K. Increasing the temperature led to a lower liquid density and a higher gas saturated density, as expected. However, when the temperature was close to the critical temperature (i.e., 540 K for n-C7H16), only the OPLS and SKS-5 models performed well in predicting the two-phase interfaces. The density gradient predicted by the original SKS model disappeared at ∼533 K (i.e., the SKS model underestimated the critical temperature), while it disappeared at ∼513 K for the TraPPE and NERD models.

FIG. 3.

Saturated density distribution curves of n-C7H16 predicted by different potential models: (a) OPLS, (b) SKS, (c) TraPPE, (d) NERD, and (e) SKS-5.

FIG. 3.

Saturated density distribution curves of n-C7H16 predicted by different potential models: (a) OPLS, (b) SKS, (c) TraPPE, (d) NERD, and (e) SKS-5.

Close modal

The VLE data taken from the National Institute of Standards and Technology (NIST)46 were adopted to verify the two-phase saturation density calculated by the five potential models in Fig. 4. The predictions obtained by using the TraPPE and NERD models deviated significantly from the experimental data in terms of the saturation density and critical point, consistent with the results shown in Fig. 3. The OPLS model overestimated the saturated liquid-phase density and underestimated the gas-phase value. Compared to the above three models, the SKS model better reproduced the two-phase saturation densities. However, the SKS model still underestimated the critical temperature. In contrast, the SKS-5 model, which combines the original SKS model and the five-parameter bound torsion sub-model, showed better agreement with the measured two-phase saturation density and critical point, especially for n-C12H26.

FIG. 4.

VLE curves predicted by different potential models for (a) n-C7H16 and (b) n-C12H26.

FIG. 4.

VLE curves predicted by different potential models for (a) n-C7H16 and (b) n-C12H26.

Close modal

The predicted and measured47,48 saturated vapor pressures for different species are compared in Fig. 5. The saturated vapor pressure was obtained as the statistical average of the pressure in the saturated gas-phase zone under the VLE state. As expected, a higher saturated gas-phase density resulted in a higher saturated vapor pressure. Consistent with the trend in the saturated gas-phase density, the predictions obtained by using the SKS-5 model agreed better with the experimental results and the NIST data than the predictions of the other models.

FIG. 5.

Comparison of the predicted and measured47,48 saturated vapor pressures for (a) n-C7H16 and (b) n-C12H26.

FIG. 5.

Comparison of the predicted and measured47,48 saturated vapor pressures for (a) n-C7H16 and (b) n-C12H26.

Close modal

We further evaluated the self-diffusion coefficient to test the ability of the potential models to reproduce the transport properties; the comparison is presented in Fig. 6. Unlike the previous results, the OPLS model showed the best agreement with the experimental data49,50 followed by the SKS-5 model. Since the potential model describes the interaction force between atoms, a greater mutual attraction between atoms leads to a higher liquid-phase density and lower gas-phase density under the VLE state. Therefore, the higher predicted liquid-phase saturated density shown in Fig. 4 indicates stronger mutual attraction and lower displacement of diffusion movement when calculating D. In summary, SKS-5, the hybrid SKS model, showed good overall accuracy for predicting the thermodynamic and transport properties of long-chain alkanes. Thus, the SKS-5 model was used in subsequent calculations.

FIG. 6.

Comparison of the predicted and measured49,50 self-diffusion coefficients for (a) n-C7H16 and (b) n-C12H26.

FIG. 6.

Comparison of the predicted and measured49,50 self-diffusion coefficients for (a) n-C7H16 and (b) n-C12H26.

Close modal

By introducing the five-parameter torsion sub-model, the bond torsion potential among four neighboring atoms can be better reproduced compared to the three-parameter torsion sub-model, especially for long-chain alkanes with complicated bond torsions. However, it must be noted that a simple UA model cannot provide accurate predictions for both short and long alkanes.42 Thus, the improved model is applicable to C7–C12, and more fundamental tests are needed beyond this range of alkane chain length.

In this section, we explore the fuel vaporization characteristics in a wide range of sub/supercritical environments in detail. Figure 7 shows the evolution of fuel molecule number, surface tension, and evaporation rate for pure n-C12H26 at T = 1200 K and P = 30 bars. The evaporation rate was determined based on the number of fuel molecules passing through the cross-section of the heating/deletion area per unit time. Initially, the liquid phase was in the subcritical stage with the presence of surface tension. The evaporation rate increased in this stage and reached the maximum value as the surface tension disappeared. Subsequently, the evaporation decreased in the diffusion-dominated supercritical stage, primarily due to the decreased density gradient caused by the disappearance of the two-phase interface. Interestingly, the total potential energy of the entire system reached its peak close to the supercritical transition. The potential energy reflects the magnitude of the intermolecular forces, which is closely related to the atom number density. The rapid decrease in potential energy after the peak value indicates the rapid transition of the system from the liquid-like state to the gas-like state. More information about the supercritical transition can be found in Sec. III C.

FIG. 7.

Pure n-C12H26 vaporization process at T = 1200 K and P = 30 bars for (a) evolution of the number of fuel molecules, surface tension, and evaporation rate and (b) typical vaporization moments.

FIG. 7.

Pure n-C12H26 vaporization process at T = 1200 K and P = 30 bars for (a) evolution of the number of fuel molecules, surface tension, and evaporation rate and (b) typical vaporization moments.

Close modal

Figure 8(a) compares the peak evaporation rates obtained under different conditions. In general, a higher ambient temperature led to a higher peak evaporation rate, as expected. However, the effect of ambient pressure on the peak evaporation rate was complicated, and the peak evaporation rate was also related to the ambient temperature. A higher proportion of the light-end component in the blended fuel resulted in more pronounced preferential evaporation in the early stage of evaporation, leading to a lower peak evaporation rate [Fig. 8(b)]. As the ambient pressure increased, the difference in the peak evaporation rate between different blending ratios decreased since evaporation led to the diffusion-dominated stage, and the evaporation rate of each species was enhanced simultaneously.

FIG. 8.

Comparison of the evaporation rates under various conditions: (a) maximum evaporation rate and (b) evaporation rate at 10 ns.

FIG. 8.

Comparison of the evaporation rates under various conditions: (a) maximum evaporation rate and (b) evaporation rate at 10 ns.

Close modal

To further demonstrate the evaporation characteristics of two-component fuels under different ambient temperatures and pressures, four typical cases were chosen for detailed analysis. Figure 9 illustrates the evaporation behaviors of a 50% n-C7H16/50% n-C12H26 blended fuel under ambient pressures of 10 and 120 bars and an ambient temperature of 800 K. As shown in Fig. 9(a), when the ambient pressure was much lower than the critical pressure of the fuel, the preferential vaporization of the light-end component (n-C7H16) was obvious (i.e., the evaporation rate is much higher) compared with the heavy-end component (n-C12H26) in the initial stage of evaporation. Furthermore, the evaporation rate of the light-end component was relatively constant without an obvious peak. As the liquid-phase temperature increased continuously, the evaporation rate of the heavy-end component increased significantly and reached a peak value after the light-end component is consumed. In this case, the evaporation belonged to subcritical-dominated evaporation since the most of the lifetime was in thev subcritical stage.

FIG. 9.

Evolution of the number of fuel molecules, evaporation rate, and mole fraction of each species for the 50% n-C7H16/50% n-C12H26 blended fuel at T = 800 K and (a) P = 10 bars and (b) P = 120 bars.

FIG. 9.

Evolution of the number of fuel molecules, evaporation rate, and mole fraction of each species for the 50% n-C7H16/50% n-C12H26 blended fuel at T = 800 K and (a) P = 10 bars and (b) P = 120 bars.

Close modal

When the ambient pressure was increased to 120 bars, the evaporation switched to supercritical-dominated evaporation with a much shorter subcritical stage, as shown in Fig. 9(b). In the initial subcritical stage, the evaporation rates for both components increased simultaneously, although the evaporation rate for the light-end component was still a little higher than that of the heavy-end component in this stage. Unlike the subcritical-dominated evaporation (where only the heavy-end component showed unimodal evaporation behavior), the evaporation rates of both the light- and heavy-end components showed unimodal characteristics, and both components coexisted throughout the evaporation process. However, the maximum evaporation rate around the transition time decreased at high ambient pressure [120 bars; also refer to Fig. 8(a)]. The greatly enhanced solubility of nitrogen in the liquid phase under high ambient pressure significantly altered the structure of the gas–liquid interface, which promoted heat transfer within the interface and resulted in a higher mass transfer rate. However, on the other hand, the higher ambient pressure resulted in higher diffusion resistance (i.e., lower diffusion capability), which reduced the evaporation rate (or, more accurately, the diffusion mixing rate) in the supercritical stage.

The results obtained at the higher ambient temperature of 1200 K are shown in Fig. 10. The general trend in the vaporization behavior with ambient pressure is nearly the same as the trend observed at the ambient temperature of 800 K [Fig. 9]. As expected, higher ambient temperature led to a shorter lifetime, especially under higher ambient pressure, due to the increased heat and mass transfer rates. The “supercritical transition rate” was also improved since the transition time decreased with increasing ambient temperature. However, it was interesting to find that the ratios of both the sub- and supercritical stage time periods to the total lifetime were not sensitive to the ambient temperature, whereas they were significantly affected by the ambient pressure. This indicates that the rates of change in the supercritical transition rate and the evaporation rate caused by the variation in ambient temperature were nearly identical. In other words, the supercritical transition characteristics were only influenced by the ambient pressure in the range of ambient temperature investigated in this study; more detailed discussion on this phenomenon can be found in Sec. III C.

FIG. 10.

Evolution of the number of fuel molecules, evaporation rate, and mole fraction of each species for the 50% n-C7H16/50% n-C12H26 blended fuel at T = 1200 K and (a) P = 10 bars and (b) P = 120 bars.

FIG. 10.

Evolution of the number of fuel molecules, evaporation rate, and mole fraction of each species for the 50% n-C7H16/50% n-C12H26 blended fuel at T = 1200 K and (a) P = 10 bars and (b) P = 120 bars.

Close modal

Finally, we discuss the evolution of the liquid/gas phase densities and gas–liquid interface thickness during the evaporation process under typical sub/supercritical conditions. This information can be obtained by curve-fitting the density distribution (Fig. 11) by using the following formula suggested by Thompson et al.:51 

ρx=12ρl+ρv12ρlρvtan2xx0d,
(10)

where ρl and ρv are the liquid- and gas-phase densities, respectively; x0 is the liquid phase center; and d is the gas–liquid interface thickness.

FIG. 11.

Curve-fitting of the density distribution at 24 ns for pure n-C12H26 at T = 800 K and P = 10 bars.

FIG. 11.

Curve-fitting of the density distribution at 24 ns for pure n-C12H26 at T = 800 K and P = 10 bars.

Close modal

The evolution of the gas/liquid-phase density and gas–liquid interface thickness with varying ambient pressure at 800 K is shown in Fig. 12. Under the subcritical environment condition (i.e., P = 10 bars), a sharp gas–liquid interface with a high density gradient across the interface was observed; this thin two-phase interface persisted throughout the lifetime of evaporation. When the pressure was increased to 120 bars, the gas–liquid interface began to dissolve, which reduced the density gradient and increased the interface thickness. After the supercritical transition at around 13 ns, the gas–liquid interface thickness increased rapidly, the gas–liquid interface promptly disappeared, and the entire system rapidly transitioned from the liquid-like state to the gas-like state. In addition, fuel blends with a higher proportion of light-end components led to a lower liquid-phase density and higher gas-phase density as a result of the higher preferential evaporation rate, as shown in Fig. 8. However, the blending ratio had little influence on the gas–liquid interface thickness, especially in the subcritical stage.

FIG. 12.

Evolution of the gas–liquid interface characteristics under different ambient pressures at T = 800 K for (a) gas/liquid-phase density and (b) gas–liquid interface thickness.

FIG. 12.

Evolution of the gas–liquid interface characteristics under different ambient pressures at T = 800 K for (a) gas/liquid-phase density and (b) gas–liquid interface thickness.

Close modal

In this section, we describe a detailed investigation of the supercritical transition characteristics under a wide variety of ambient conditions for different blended fuels (0% n-C7H16 means pure n-C12H26). In a similar study, Mo and Qiao32 introduced a dimensionless transition time (τT) defined as τT = tT/tL to quantitatively describe the proportion of the sub/supercritical stage illustrated in Fig. 7 out of the total lifetime (e.g., τT = 0.2 means 20% of the lifetime is in the subcritical stage, while the remaining 80% is in the supercritical stage). We employed this same parameter in this study; however, we updated the criterion used to determine the transition time tT. Mo and Qiao32 determined the transition time as the point of zero surface tension based on the surface tension history. However, the accurate calculation of surface tension at the atomic level by using MD simulations under supercritical conditions is quite challenging.52 Alternatively, we propose a new criterion based on the maximum potential energy and evaporation rate, which can accurately reflect the thermodynamic state of the entire system. Figure 13 shows the evolution of the potential energy and evaporation rate under different conditions. The maximum potential energy clearly occurred almost simultaneously with the point of the maximum evaporation rate. The average time at which the potential energy and evaporation rate reach their maximum values was chosen as the transition time in this study.

FIG. 13.

Evolution of the potential energy and evaporation rate under different ambient pressures for (a) 0% n-C7H16 at T = 800 K, (b) 0% n-C7H16 at T = 1200 K, (c) 25% n-C7H16 at T = 800 K, and (d) 25% n-C7H16 at T = 1200 K.

FIG. 13.

Evolution of the potential energy and evaporation rate under different ambient pressures for (a) 0% n-C7H16 at T = 800 K, (b) 0% n-C7H16 at T = 1200 K, (c) 25% n-C7H16 at T = 800 K, and (d) 25% n-C7H16 at T = 1200 K.

Close modal

To validate the new transition criterion, the evolutions of the gas–liquid interface on the temperature–mole fraction (Tx) VLE diagram under different ambient conditions are shown in Fig. 14. The left and right curves are, respectively, the gas and liquid equilibrium lines, and the area between the two curves is the gas–liquid interface area, in which liquid and gas are in equilibrium states. In the subcritical stage, the Tx line will cross the two-phase region; when the supercritical transition occurs, the Tx line will be close to the critical point; afterward, it is only located in the single-phase region. It can be seen that the transition time obtained by the new transition criterion is very close to the critical point, proving that the new transition criterion is reliable.

FIG. 14.

The evolution of gas–liquid interface on the temperature–mole fraction VLE diagram for pure n-C12H26 at (a) T = 800 K and P = 90 bars and (b) T = 1200 K and P = 90 bars.

FIG. 14.

The evolution of gas–liquid interface on the temperature–mole fraction VLE diagram for pure n-C12H26 at (a) T = 800 K and P = 90 bars and (b) T = 1200 K and P = 90 bars.

Close modal

Figures 15(a) and 15(b), respectively, show the distributions of liquid-phase lifetime and supercritical transition time under different ambient temperatures and pressures. Increasing the ambient temperature effectively increased the heat transfer from the surrounding environment to the interior of the liquid phase, resulting in a shorter lifetime and transition time. However, the effect of ambient pressure on the lifetime was much more complicated. When the ambient pressure was increased from 10 to 30 bars, the liquid phase lifetime was reduced due to the enhanced evaporation rate; however, further increasing the pressure led to a longer lifetime, consistent with the trend in evaporation rate shown in Fig. 8. For transition time, higher pressure yields shorter results, while it has a limited effect when further increasing the ambient pressure. Overall, the effects of the ambient pressure on the lifetime and transition time were weakened as the ambient temperature increased from 800 to 1200 K.

FIG. 15.

Comparison of the vaporization-relevant times under a range of ambient temperature and pressure conditions: (a) liquid-phase lifetime and (b) supercritical transition time.

FIG. 15.

Comparison of the vaporization-relevant times under a range of ambient temperature and pressure conditions: (a) liquid-phase lifetime and (b) supercritical transition time.

Close modal

Finally, Fig. 16 shows the distribution of τT by combining the results shown in Figs. 15(a) and 15(b). Surprisingly, Fig. 16 shows that τT is nearly independent of the ambient temperature and fuel composition; instead, τT mainly varies with the ambient pressure. Increasing the ambient pressure effectively reduced τT, leading to more supercritical-dominated vaporization. However, Mo and Qiao32 reported that τT depends strongly on the ambient temperature and pressure along with the fuel type and that the pressure required for the transition increases with the number of carbon atoms in n-alkanes.

FIG. 16.

Dimensionless supercritical transition time as a function of pressure.

FIG. 16.

Dimensionless supercritical transition time as a function of pressure.

Close modal

The above discrepancy might be explained in two ways. The first explanation involves the simulation configuration for vaporization, as shown in Fig. 2. In this study, the evaporated molecules were deleted at two sides of the box to reproduce the real evaporation process (i.e., fuel evaporating into an infinite space). Thus, the evaporated fuel molecules did not affect the subsequent evaporation process. However, Mo and Qiao32 did not provide information on the fate of the evaporated fuel molecules; thus, the evaporated fuel molecules may have condensed when they moved back to the liquid phase, thereby influencing the subsequent evaporation process. The second explanation involves the transition criterion used to determine the transition time. As discussed above, it is difficult to accurately calculate the surface tension under supercritical conditions by using MD simulations. According to our preliminary results, the calculated surface tension fluctuated greatly, especially at the near-zero position, which can easily result in inaccurate transition times. The transition criterion used in this study is more reasonable and can more accurately determine the supercritical transition time based on the evolution of potential energy and evaporation rate, as illustrated in Figs. 13 and 14.

The following formula for transition time was obtained by curve-fitting the results shown in Fig. 16:

τT=a+bPc,
(11)

where P is the ambient pressure (unit: bar); and a, b, and c are the curve-fitting constants with values of a = 1.100 57, b = 0.019 07, and c = −1.189 10. To further validate this formula, additional simulations were conducted at P = 200 bars for pure n-C7H16 under ambient temperatures of 800 and 1200 K, resulting in a predicted τT of 0.17. The corresponding τT obtained by using Eq. (11) is 0.15, which is close to the value predicted by MD. Nevertheless, further study is required, especially for other blended fuels such as alcohol/n–alkane or benzene/n–alkane mixtures.

In this study, we evaluated four commonly used UA potential models for n-alkanes in MD simulations (OPLS, TraPPE, NERD, and SKS) and introduced a hybrid model (SKS-5). Based on the hybrid model, we explored the vaporization and phase transition characteristics of n-alkane blended fuels under a wide range of ambient conditions. The main conclusions are summarized as follows:

  • The hybrid potential model, SKS-5, which combines the original SKS model and the five-parameter torsion sub-model, showed overall good accuracy in predicting the thermodynamic and transport properties of long-chain n-alkanes (C7–C12).

  • Higher ambient temperature led to a higher evaporation rate and a shorter lifetime, as expected, while the effect of the ambient pressure was much more complicated. When the ambient pressure was increased from 10 to 30 bars, the lifetime was reduced due to the enhanced evaporation rate; however, further increasing the pressure led to a longer lifetime. The effect of fuel composition on the evaporation characteristics became weaker as the ambient pressure increased.

  • In subcritical environments, a sharp gas–liquid interface with a high density gradient persisted throughout the evaporation process. The fuel composition had little influence on the interfacial properties (e.g., interface thickness). The vaporization belonged to typical two-phase evaporation. For blended fuels, the preferential vaporization of the light-end component was obvious, and the evaporation rate of the heavy-end component reached a peak after the light-end component was totally consumed.

  • Under supercritical conditions, the gas–liquid interface still existed before the supercritical transition. After the transition, the interface dissolved rapidly in association with increased interface thickness and finally disappeared. Evaporation in this case belonged to diffusion-dominated vaporization. The evaporation rates for both the light- and heavy-end components increased simultaneously, and both components coexisted throughout the evaporation process.

  • A new supercritical transition criterion based on the maximum potential energy and evaporation rate was proposed. Based on this new criterion, the dimensionless transition time, which reflects the proportion of the sub/supercritical stage within the lifetime, was found to be nearly independent of the ambient temperature and fuel composition; instead, it mainly varied with the ambient pressure. Finally, a formula was obtained by curve-fitting to describe the variation in the dimensionless transition time with ambient pressure.

This work was supported by the National Natural Science Foundation of China (Grant No. 52106154) and the Natural Science Foundation of Jiangsu Province, China (Grant No. BK20190855).

The authors have no conflicts to disclose.

The data that support the findings of this study are available from the corresponding authors upon reasonable request.

1.
F.
Poursadegh
,
J. S.
Lacey
,
M. J.
Brear
, and
R. L.
Gordon
, “On the fuel spray transition to dense fluid mixing at reciprocating engine conditions,”
Energy Fuels
31
,
6445
6454
(
2017
).
2.
M.
Oschwald
,
J. J.
Smith
,
R.
Branam
,
J.
Hussong
,
A.
Schik
,
B.
Chehroudi
, and
D.
Talley
, “Injection of fluids into supercritical environments,”
Combust. Sci. Technol.
178
,
49
100
(
2006
).
3.
W.
Mayer
and
H.
Tamura
, “Propellant injection in a liquid oxygen/gaseous hydrogen rocket engine,”
J. Propul. Power
12
,
1137
1147
(
1996
).
4.
C.
Crua
,
J.
Manin
, and
L. M.
Pickett
, “On the transcritical mixing of fuels at diesel engine conditions,”
Fuel
208
,
535
548
(
2017
).
5.
Z.
Falgout
,
M.
Rahm
,
D.
Sedarsky
, and
M.
Linne
, “Gas/fuel jet interfaces under high pressures and temperatures,”
Fuel
168
,
14
21
(
2016
).
6.
M.
Wensing
,
T.
Vogel
, and
G.
Götz
, “Transition of diesel spray to a supercritical state under engine conditions,”
Int. J. Engine Res.
17
,
108
119
(
2016
).
7.
Z.
Falgout
,
M.
Rahm
,
Z.
Wang
, and
M.
Linne
, “Evidence for supercritical mixing layers in the ECN Spray A,”
Proc. Combust. Inst.
35
,
1579
1586
(
2015
).
8.
F.
Ries
and
A.
Sadiki
, “Supercritical and transcritical turbulent injection processes: Consistency of numerical modeling,”
Atomization Sprays
31
,
37
(
2021
).
9.
P.
Koukouvinis
,
C.
Rodriguez
,
J.
Hwang
,
I.
Karathanassis
,
M.
Gavaises
, and
L.
Pickett
, “Machine learning and transcritical sprays: A demonstration study of their potential in ECN Spray-A,”
Int. J. Engine Res.
(published onine,
2021
).
10.
R. N.
Dahms
and
J. C.
Oefelein
, “Liquid jet breakup regimes at supercritical pressures,”
Combust. Flame
162
,
3648
3657
(
2015
).
11.
R. N.
Dahms
and
J. C.
Oefelein
, “On the transition between two-phase and single-phase interface dynamics in multicomponent fluids at supercritical pressures,”
Phys. Fluids
25
,
092103
(
2013
).
12.
R. N.
Dahms
,
J.
Manin
,
L. M.
Pickett
, and
J. C.
Oefelein
, “Understanding high-pressure gas-liquid interface phenomena in diesel engines,”
Proc. Combust. Inst.
34
,
1667
1675
(
2013
).
13.
C.
Morin
,
C.
Chauveau
,
P.
Dagaut
,
I.
Gökalp
, and
M.
Cathonnet
, “Vaporization and oxidation of liquid fuel droplets at high temperature and high pressure: application to N-alkanes and vegetable oil methyl esters,”
Combust. Sci. Technol.
176
,
499
529
(
2004
).
14.
H.
Nomura
,
T.
Murakoshi
,
Y.
Suganuma
,
Y.
Ujiie
,
N.
Hashimoto
, and
H.
Nishida
, “Microgravity experiments of fuel droplet evaporation in sub- and supercritical environments,”
Proc. Combust. Inst.
36
,
2425
2432
(
2017
).
15.
K.
Harstad
and
J.
Bellan
, “An all-pressure fluid drop model applied to a binary mixture: Heptane in nitrogen,”
Int. J. Multiphase Flow
26
,
1675
1706
(
2000
).
16.
G. S.
Zhu
and
S. K.
Aggarwal
, “Transient supercritical droplet evaporation with emphasis on the effects of equation of state,”
Int. J. Heat Mass Transfer
43
,
1157
1171
(
2000
).
17.
R.
Branam
and
W.
Mayer
, “Length scales in cryogenic injection at supercritical pressure,”
Exp Fluids
33
,
422
428
(
2002
).
18.
B.
Chehroudi
,
D.
Talley
, and
E.
Coy
, “Visual characteristics and initial growth rates of round cryogenic jets at subcritical and supercritical pressures,”
Phys. Fluids
14
,
850
861
(
2002
).
19.
J.
Manin
,
M.
Bardi
,
L. M.
Pickett
,
R. N.
Dahms
, and
J. C.
Oefelein
, “Microscopic investigation of the atomization and mixing processes of diesel sprays injected into high pressure and temperature environments,”
Fuel
134
,
531
543
(
2014
).
20.
S.
Jafari
,
H.
Gaballa
,
C.
Habchi
, and
J.-C.
de Hemptinne
, “Towards understanding the structure of subcritical and transcritical liquid–gas interfaces using a tabulated real fluid modeling approach,”
Energies
14
,
5621
(
2021
).
21.
B. M.
Ningegowda
,
F. N. Z.
Rahantamialisoa
,
A.
Pandal
,
H.
Jasak
,
H. G.
Im
, and
M.
Battistoni
, “Numerical modeling of transcritical and supercritical fuel injections using a multi-component two-phase flow model,”
Energies
13
,
5676
(
2020
).
22.
P.
Koukouvinis
,
A.
Vidal-Roncero
,
C.
Rodriguez
,
M.
Gavaises
, and
L.
Pickett
, “High pressure/high temperature multiphase simulations of dodecane injection to nitrogen: Application on ECN Spray-A,”
Fuel
275
,
117871
(
2020
).
23.
H.
Müller
,
C. A.
Niedermeier
,
J.
Matheis
,
M.
Pfitzner
, and
S.
Hickel
, “Large-eddy simulation of nitrogen injection at trans- and supercritical conditions,”
Phys. Fluids
28
,
015102
(
2016
).
24.
L.
Qiu
and
R. D.
Reitz
, “An investigation of thermodynamic states during high-pressure fuel injection using equilibrium thermodynamics,”
Int. J. Multiphase Flow
72
,
24
38
(
2015
).
25.
J.
Zhang
,
X.
Zhang
,
T.
Wang
, and
X.
Hou
, “A numerical study on jet characteristics under different supercritical conditions for engine applications,”
Appl. Energy
252
,
113428
(
2019
).
26.
P.
Yi
,
S.
Yang
,
C.
Habchi
, and
R.
Lugo
, “A multicomponent real-fluid fully compressible four-equation model for two-phase flow with phase change,”
Phys. Fluids
31
,
026102
(
2019
).
27.
S.
Yang
,
P.
Yi
, and
C.
Habchi
, “Real-fluid injection modeling and LES simulation of the ECN Spray A injector using a fully compressible two-phase flow approach,”
Int. J. Multiphase Flow
122
,
103145
(
2020
).
28.
B.
Xu
,
Y.
Zhu
,
H.
Jin
,
Y.
Guo
, and
J.
Fan
, “Transcritical transition of the fluid around the interface,”
Phys. Fluids
33
,
122106
(
2021
).
29.
H.
Yamaguchi
and
G.
Kikugawa
, “Molecular dynamics study on flow structure inside a thermal transpiration flow field,”
Phys. Fluids
33
,
012005
(
2021
).
30.
A.
Tokunaga
and
T.
Tsuruta
, “Nonequilibrium molecular dynamics study on energy accommodation coefficient on condensing liquid surface—Molecular boundary conditions for heat and mass transfer,”
Phys. Fluids
32
,
112011
(
2020
).
31.
M.
Foroutan
,
F.
Esmaeilian
, and
M. T.
Rad
, “The change in the wetting regime of a nanodroplet on a substrate with varying wettability: A molecular dynamics investigation,”
Phys. Fluids
33
,
032017
(
2021
).
32.
G.
Mo
and
L.
Qiao
, “A molecular dynamics investigation of n-alkanes vaporizing into nitrogen: Transition from subcritical to supercritical,”
Combust. Flame
176
,
60
71
(
2017
).
33.
B.
Smit
,
S.
Karaborni
, and
J. I.
Siepmann
, “Computer simulations of vapor–liquid phase equilibria of n-alkanes,”
J. Chem. Phys.
102
,
2126
2140
(
1995
).
34.
G.
Xiao
,
K. H.
Luo
,
X.
Ma
, and
S.
Shuai
, “A molecular dynamics study of fuel droplet evaporation in sub- and supercritical conditions,”
Proc. Combust. Inst.
37
,
3219
3227
(
2019
).
35.
W. L.
Jorgensen
,
J. D.
Madura
, and
C. J.
Swenson
, “Optimized intermolecular potential functions for liquid hydrocarbons,”
J. Am. Chem. Soc.
106
,
6638
6646
(
1984
).
36.
Z.
Wang
,
L.
Zhou
,
G.
Shu
, and
H.
Wei
, “Droplet evaporation and phase transition modes in supercritical environment by molecular dynamic simulation,”
Phys. Fluids
33
,
062001
(
2021
).
37.
Y.
Zhang
,
M.
Jia
,
P.
Yi
,
Y.
Chang
,
Z.
He
, and
Q.
Wang
, “A molecular dynamics study of binary-component n-alkane fuel vaporization characteristics at sub/supercritical nitrogen environments,”
Proc. Combust. Inst.
38
,
3303
3312
(
2021
).
38.
Y.
Gong
,
G.
Xiao
,
X.
Ma
,
K. H.
Luo
,
S.
Shuai
, and
H.
Xu
, “Phase transitions of multi-component fuel droplets under sub- and supercritical conditions,”
Fuel
287
,
119516
(
2021
).
39.
M. G.
Martin
and
J. I.
Siepmann
, “Transferable potentials for phase equilibria. 1. United-atom description of n-alkanes,”
J. Phys. Chem. B
102
,
2569
2577
(
1998
).
40.
Z.
Wang
,
W.
Zhao
,
L.
Zhou
,
G.
Shu
, and
H.
Wei
, “A molecular dynamic study of evaporation/supercritical-transition inter-relationship and multicomponents interaction for alkane/alcohol droplets,”
Phys. Fluids
34
,
022002
(
2022
).
41.
S.
Chakraborty
and
L.
Qiao
, “Molecular investigation of sub-to-supercritical transition of hydrocarbon mixtures: Multi-component effect,”
Int. J. Heat Mass Transfer
145
,
118629
(
2019
).
42.
S. K.
Nath
,
F. A.
Escobedo
, and
J. J.
de Pablo
, “On the simulation of vapor–liquid equilibria for alkanes,”
J. Chem. Phys.
108
,
9905
9911
(
1998
).
43.
J. L.
Rivera
and
J.
Alejandre
, “Vapor–liquid equilibrium simulations of nitrogen and n-alkane binary mixtures,”
Colloids Surf., A
207
,
223
228
(
2002
).
44.
V.
Kalyanasundaram
,
D. E.
Spearot
, and
A. P.
Malshe
, “Molecular dynamics simulation of nanoconfinement induced organization of n-decane,”
Langmuir
25
,
7553
7560
(
2009
).
45.
P.
Padilla
and
S.
Toxvaerd
, “Self-diffusion in n-alkane fluid models,”
J. Chem. Phys.
94
,
5650
5654
(
1991
).
46.
P. J.
Linstrom
and
W. G.
Mallard
, “The NIST chemistry webbook: A chemical data resource on the internet,”
J. Chem. Eng. Data
46
,
1059
1063
(
2001
).
47.
L. A.
Weber
, “Vapor pressure of heptane from the triple point to the critical point,”
J. Chem. Eng. Data
45
,
173
176
(
2000
).
48.
D. L.
Morgan
and
R.
Kobayashi
, “Direct vapor pressure measurements of ten n-alkanes m the C10–C28 range,”
Fluid Phase Equilib.
97
,
211
242
(
1994
).
49.
E.
Fishman
, “Self-diffusion in liquid normal pentane and normal heptane,”
J. Phys. Chem.
59
,
469
472
(
1955
).
50.
K. R.
Harris
,
B.
Ganbold
, and
W. S.
Price
, “Viscous calibration liquids for self-diffusion measurements,”
J. Chem. Eng. Data
60
,
3506
3517
(
2015
).
51.
S. M.
Thompson
,
K. E.
Gubbins
,
J. P. R. B.
Walton
,
R. A. R.
Chantry
, and
J. S.
Rowlinson
, “A molecular dynamics study of liquid drops,”
J. Chem. Phys.
81
,
530
542
(
1984
).
52.
A.
Ghoufi
,
P.
Malfreyt
, and
D. J.
Tildesley
, “Computer modelling of the surface tension of the gas–liquid and liquid–liquid interface,”
Chem. Soc. Rev.
45
,
1387
1409
(
2016
).