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.

## I. INTRODUCTION

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*-C_{7}H_{16}) 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 Bellan^{15} proposed a basic theory for modeling single-droplet vaporization under sub/supercritical conditions, and Zhu and Aggarwal^{16} 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*-C_{12}H_{26})] 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 modeling^{20–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 (N_{2}) jets under different N_{2} 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 (CH_{3}) and methylene (CH_{2}) 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 Qiao^{32} 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*-C_{12}H_{26} 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 model^{39} 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 Qiao^{41} 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.

## II. MOLECULAR DYNAMICS SIMULATION

### A. Potential model

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:

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

where *σ*_{i,j} and *ϵ*_{i,j} are the size and energy parameters, respectively, and *r*_{c} 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,

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

where *k*_{b} 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:

where *k*_{θ} is the bending coefficient and $\theta 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:

where *k*_{1}, *k*_{2}, and *k*_{3} are the model coefficients.

The Rivera and Alejandre model^{43} was adopted to describe the interatomic potential of the N_{2} 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 N_{2} atoms was fixed at 0.108 97 nm.

### B. Hybrid SKS model

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*-C_{7}H_{16}. However, the SKS model still fails to quantitatively reproduce the physical properties of longer-chain alkanes (e.g., *n*-C_{12}H_{26}). 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,

where *A*_{n} 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.

Potential parameters . | OPLS . | SKS . | TraPPE . | NERD . | SKS-5 . |
---|---|---|---|---|---|

$\u03f5CH3kcal/mol$ | 0.175 | 0.2265 | 0.195 | 0.207 | 0.2265 |

$\u03f5CH2kcal/mol$ | 0.118 | 0.0934 | 0.091 | 0.091 | 0.0934 |

$\sigma CH3nm$ | 0.3905 | 0.393 | 0.375 | 0.391 | 0.393 |

$\sigma CH2nm$ | 0.3905 | 0.393 | 0.395 | 0.393 | 0.393 |

$kbkcal/mol\u22c5nm2$ | 52012 | 19 177 | |||

$rieqnm$ | 0.1526 | 0.154 | |||

$k\theta kcal/mol\u22c5rad2$ | 124.15 | 124.15 | |||

$\theta ieq\xb0$ | 112 | 114 | |||

$kn/Ankcal/mol$ | $k1=1.411k2=\u22120.271k3=3.144$ | $A1=4.825A2=0.162A3=\u22126.224A4=\u22120.325A5=\u22120.503$ |

Potential parameters . | OPLS . | SKS . | TraPPE . | NERD . | SKS-5 . |
---|---|---|---|---|---|

$\u03f5CH3kcal/mol$ | 0.175 | 0.2265 | 0.195 | 0.207 | 0.2265 |

$\u03f5CH2kcal/mol$ | 0.118 | 0.0934 | 0.091 | 0.091 | 0.0934 |

$\sigma CH3nm$ | 0.3905 | 0.393 | 0.375 | 0.391 | 0.393 |

$\sigma CH2nm$ | 0.3905 | 0.393 | 0.395 | 0.393 | 0.393 |

$kbkcal/mol\u22c5nm2$ | 52012 | 19 177 | |||

$rieqnm$ | 0.1526 | 0.154 | |||

$k\theta kcal/mol\u22c5rad2$ | 124.15 | 124.15 | |||

$\theta ieq\xb0$ | 112 | 114 | |||

$kn/Ankcal/mol$ | $k1=1.411k2=\u22120.271k3=3.144$ | $A1=4.825A2=0.162A3=\u22126.224A4=\u22120.325A5=\u22120.503$ |

### C. Simulation configuration

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 *N*_{s} 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 *k*th slice could be determined by the following equations:

where $Nk$ is the atom number in the *k*th slice; $Vk$ is the corresponding slice volume calculated as $Vk=LxLyLzNs$, where *L*_{x}, *L*_{y}, and *L*_{z} are the box lengths in the *x*, *y*, and *z* directions, respectively; and *m*_{i}, *v*_{i}, *r*_{i}, and *F*_{i} are the mass, velocity, radius vector, and force of the *i*th atom, respectively.

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:

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.

Configuration . | Box length (nm^{3})
. | Number of molecules . | Temperature range (K) . |
---|---|---|---|

Thermodynamic | |||

(n-C_{7}H_{16}) | 30 ×8× 8 | 2048 | 293–533 |

Thermodynamic (n-C_{12}H_{26}) | 40 ×7.7× 7.7 | 1536 | 450–650 |

Transport (n-C_{7}H_{16}) | 5 ×5× 5 | 500 | 273–368 |

Transport (n-C_{12}H_{26}) | 6 ×6× 6 | 576 | 308–364 |

Configuration . | Box length (nm^{3})
. | Number of molecules . | Temperature range (K) . |
---|---|---|---|

Thermodynamic | |||

(n-C_{7}H_{16}) | 30 ×8× 8 | 2048 | 293–533 |

Thermodynamic (n-C_{12}H_{26}) | 40 ×7.7× 7.7 | 1536 | 450–650 |

Transport (n-C_{7}H_{16}) | 5 ×5× 5 | 500 | 273–368 |

Transport (n-C_{12}H_{26}) | 6 ×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*-C_{7}H_{16} and *n*-C_{12}H_{26} in different mole ratios was placed at the center of a three-dimensional cuboid box with dimensions of 265 × 100 × 100 nm^{3}, and N_{2} 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}

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 A^{19}). Meanwhile, the N_{2} 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*-C_{7}H_{16} and *n*-C_{12}H_{26}. 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.

## III. RESULTS AND DISCUSSION

### A. Comparison of different potential models

This section describes the evaluation of three physical properties (the VLE curve, saturated vapor pressure, and self-diffusion coefficient) of *n*-C_{7}H_{16} and *n*-C_{12}H_{26} to test the performance of the five potential models.

The saturated density distribution curves of *n*-C_{7}H_{16} 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*-C_{7}H_{16}), 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.

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*-C_{12}H_{26}.

The predicted and measured^{47,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.

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 data^{49,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.

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 C_{7}–C_{12}, and more fundamental tests are needed beyond this range of alkane chain length.

### B. Vaporization characteristics under typical sub/supercritical conditions

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*-C_{12}H_{26} 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.

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.

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*-C_{7}H_{16}/50% *n*-C_{12}H_{26} 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*-C_{7}H_{16}) was obvious (i.e., the evaporation rate is much higher) compared with the heavy-end component (*n*-C_{12}H_{26}) 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.

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.

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}

where *ρ*_{l} and *ρ*_{v} are the liquid- and gas-phase densities, respectively; *x*_{0} is the liquid phase center; and *d* is the gas–liquid interface thickness.

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.

### C. Characteristics of the supercritical transition under different ambient conditions

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*-C_{7}H_{16} means pure *n*-C_{12}H_{26}). In a similar study, Mo and Qiao^{32} introduced a dimensionless transition time (*τ*_{T}) defined as *τ*_{T} = *t*_{T}/*t*_{L} 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 *t*_{T}. Mo and Qiao^{32} 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.

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.

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.

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 Qiao^{32} 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.

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 Qiao^{32} 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:

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*-C_{7}H_{16} 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.

## IV. CONCLUSIONS

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 (C_{7}–C_{12}).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.

## ACKNOWLEDGMENTS

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).

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

## DATA AVAILABILITY

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

## REFERENCES

_{10}–C

_{28}range,”