A molecular dynamics study of evaporation mode transition of hydrocarbon fuels under supercritical conditions

The mode transition of evaporation for single-and multi-component hydrocarbon fuels is investigated at the molecular level. This study scrutinizes ﬁrst the subcritical and supercritical evaporation of n-hexadecane droplets and liquid ﬁlms by molecular dynamics (MD) simulations. The mode regime map of n-hexadecane droplets is obtained. Then the mode transition of evaporation of a three-component droplet and a six-component droplet is studied. A critical dimensionless number τ 0.9P of 0.5 based on the average displacement increment (ADI) of fuel atoms is used to identify the evaporation mode transition of fuels with any type and number of components. It is found that in the diffusion mode of evaporation, the entropy becomes the dominant factor in the evaporation of fuels, and the disorder of the fuel molecules increases signiﬁcantly compared with that in the classic evaporation mode. Compared with the case of the quiescent droplet, with increasing relative velocity between the droplet and the ambient gas, the mode transition becomes easier, although this is a non-linear process. Fuel droplets and liquid ﬁlms with different initial sizes are investigated to understand the size effect. In addition, for the same ambient temperature and pressure, the smaller the normalized speciﬁc heat transfer surface area of the fuel is, the easier the mode transition of evaporation is. A correlation was proposed to compare the possibility of mode transition of evaporation for single-and multi-component fuels. © 2022 The Author(s). Published by Elsevier Inc. on behalf of The Combustion Institute. This is an open access article under the CC BY license ( http://creativecommons.org/licenses/by/4.0/ )


Introduction
Combustion of hydrocarbon fuels under ever high temperature and ever high pressure conditions is an important measure to improve the efficiency of power devices such as liquid rocket engines and internal combustion engines.Such temperature and pressure often exceed the critical point of the injected fuel so that supercritical combustion takes place in the combustion chamber [1] .Due to the drastic heat and mass transfer, the evaporation mode of fuels may undergo a transition from classical evaporation to diffusive mixing via the transitional mixing [2] .Therefore, understanding the evaporation mode transition of fuels and its effect on the subsequent air-fuel mixing and combustion processes has motivated many research efforts [3][4][5][6] .
Subcritical and supercritical evaporation has significant differences [7] .There is a sharp vapor-liquid interface under subcritical conditions, and the surface tension plays an important role.But in the supercritical evaporation, the vapor-liquid interface almost disappears.Several reviews have discussed the theory and modeling of supercritical droplets and jets [3][4][5][6] .According to these, there are still unanswered key questions about the physics of the mode transition of evaporation for liquid fuels.CFD simulations have been extensively attempted to answer these questions [8 , 9] .However, the dramatic changes in physical properties of the fluids near its thermodynamic critical point may invalidate the conventional CFD approaches [10] , especially for complex evaporation systems with multi-fuel components [11] .In particular, the basic assumption of vapor-liquid equilibrium (VLE) adopted by conventional CFD has been shown to be invalid under supercritical and multicomponent conditions [12] .Similar issues have been encountered in the theoretical analysis of supercritical evaporation.Dahms et al. [13] applied the linear gradient theory and proposed a Knudsen-number (Kn) criterion to delineate the spray regime diagram.When Kn is less than 0.1, the mixing of fuels was thought to be in the continuum regime.However, such theoretical analyses are usually based on the steady-state hypothesis and fail to account for the unsteady nature of mode transition of evaporation [14] .Due to the limitations of experimental devices [2] , there are limited experimental data available for the validation of models [2 , 15 , 16] .In particular, limited by the intrinsic optical effects under high-pressure conditions, the response of imaging system, and image processing techniques for moving objects [2] , it is hard to experimentally investigate the dynamic details in the mode transition from evaporation to diffusion.In recent years, MD simulations have started to make an impact in this field [17 , 18] .MD is based on Newton's second law and does not require other thermodynamic assumptions other than the potential function, and the simulation results are statistically treated [19] .MD has unique advantages for the study of interfacial transition [20] , considering that the thickness of the vapor-liquid interface is usually only a few nanometers [13] , which is hard to resolve based on current experimental techniques.Moreover, MD does not need to make any assumptions about the interface, and therefore can easily handle the disappearance of the nanoscale interface, which is beyond the conventional CFD [20] .Therefore, MD is ideally suited for investigating the physical nature of the evaporation and mode transition.
It also should be noted that MD has two noticeable limitations.Firstly, MD simulations are computationally expensive and require supercomputing resources, so it is not yet possible to directly simulate the evaporation of macroscopic droplets using MD [17] .Generally speaking, the spatial and temporal scales of such studies are currently limited to a few hundred nanometers and nanoseconds [10] , respectively, and the number of molecules studied is limited to millions [20] .In this context, it should be stated how relevant the MD results on evaporation are to the evaporation behaviors occurring in real combustion chambers [21] .In fact, researchers have been trying to answer this question since MD was applied to evaporation studies, including this work.Following the MD work of Little [22] , Kaltz [23] conducted the MD scaling for droplets at a micron level using the Long-Micci method and investigated the validity of extrapolating MD results to macroscopic regimes by scaling the intermolecular potential parameters under subcritical and supercritical conditions.Qiao et al. [14] demonstrated that MD results could be used to approximate the transition from subcritical to supercritical evaporation occurring in macroscopic experiments.Recently, as reported [24] , atomistic simulations have revealed fundamental mechanisms that also exist in macroscopic phenomena, which complements macroscopic experiments.Secondly, the accuracy and reliability of MD results depend on the interatomic potentials used [21] .However, as the interatomic potentials are determined by the fundamental interactions between atoms/molecules, there are far less empirical assumptions and uncertainties in interatomic potentials than in macroscopic models.Thus, the MD simulation results are physics-based rather than a priori [21] .Several potential models have been successfully developed for complex hydrocarbon molecules [14] , which can be divided into two types: the All-Atom Model (AAM) and the United Atom Model (UAM).As reported [25] , the computational time for MD simulations is roughly proportional to the square of the total number of interaction sites.The AAM regards each atom as one site and is suitable for systems with a small number of molecules and a higher requirement for computational accuracy [25] .On the other hand, the UAM's treatment of each pair of CH 3 and CH 2 as one site, which keeps a good balance between computational accuracy and efficiency [10 , 25 , 26] , is widely used for large molecule systems, such as hydrocarbons.Jorgensen et al. [25] compared the simulation results based on the AAM and the UAM.With suitable parameters, the AAM and the UAM could both give highly accurate descriptions of fluids [25] .Despite these limitations, the physics obtained in MD simulations will be helpful for understanding the mechanism of mode transition of evaporation at the atomic level [10 , 20 , 24 , 26] .MD is a complement to present experimental and CFD techniques [24 , 27] .As reported [21 , 27 , 28] , the fluid property data obtained in MD simulations can provide a reliable benchmark for CFD to employ the correct equation of state (EoS) and mixing rules of components under subcritical and supercritical conditions.In addition, the knowledge of droplet evaporation will be expanded by such MD simulations, which is valuable for modern combustion devices [17] .In summary, MD is valuable for further understanding the physics of evaporation [17 , 20 , 24 , 26] , examining assumptions in conventional CFD [14] and providing benchmarks for CFD models [27 , 28] .
Although there have been progresses in understanding mode transition of evaporation using MD [10 , 14 , 29] , the effect of the relative velocity between the droplet (or liquid film) and the ambient gas on the mode transition of evaporation under supercritical conditions has not been sufficiently studied.In fact, in a real engine combustion chamber, the relative velocity between the fuel and the ambient gas is very large [16] .In optical experiments, this velocity effect is hard to investigate because of the difficulties in tracking the boundary of moving droplets under supercritical conditions, although there are some useful observations about the atomization regimes in subcritical environments [2] .Based on the theoretical analysis of the time-scale of evaporation and breakup, Poursadegh et al. [16] found that dense mixing could occur at lower ambient pressures and temperatures when the injected fuel velocity became larger.Ray et al. [30] investigated the effects of convection velocity on the evaporation of fuel droplets and found that the higher velocity gradient near the droplet interface promoted the evaporation with ambient pressures ranging from 1 bar to 55 bar.Recently, Liu et al. [29] studied the effect of forced convection on the preferential evaporation of light components in multi-component fuel droplets in CFD simulations and found that there were two different controlling regimes at different convective intensities.However, this study was focused on cases at subcritical pressures, leaving unanswered questions for supercritical conditions.Given these situations, it is essential to conduct further investigations of the relative velocity effects on the mode transition of evaporation.
As mentioned, MD can only resolve the evaporation and mode transition of droplets/films at nanoscales.However, the droplets generated by jets after atomization in the combustion chamber of an actual engine are usually micron-sized.Therefore, it is necessary to investigate the size effect to bridge MD and experimental studies.In fact, existing studies have inconsistent views on this issue.Qiao et al. [14] investigated the size effect of liquid films and found that the results of nanoscale films were size-independent.They attributed this to interfacial dynamics and thought that the interfacial elements were responsible for the phase transition [14] .But in another MD research, Xiao et al. [10] studied the supercritical evaporation of n-dodecane droplets with different initial diameters (from 21 nm to 43 nm), and found that when the initial droplet size increased, the supercritical transition could occur at lower ambient pressures.Recently similar conclusions were conducted by Wang et al. [31] .This was attributed to the fact that larger droplets had relatively longer evaporation lifetime to achieve the mode transition.Clearly, this topic needs further investigations.
Actual fuels are usually complex mixtures with hundreds of components [32] .However, most existing MD studies have only considered single-component fuels [20] and two-component fuels [12 , 33-35 ].Chakraborty et al. [34] studied the evaporation of a two-component fuel liquid film in supercritical nitrogen, and found that under supercritical conditions, light components and heavy ones dominated the evaporation simultaneously, which was different from the preferential evaporation under subcritical conditions.Zhang et al. [12] studied the evaporation of a two-component liquid film composed of n-heptane and n-dodecane in subcritical and supercritical conditions.They concluded that compared with the mixing ratio of fuel components, the ambient temperature and pressure had greater impacts on the mode transition of evaporation [12] .In a rare example, Gong et al. [20 , 26] studied molecular interactions between individual components in multi-component fuel droplets and conducted Voronoi tessellation analyses for fuels.However, many key questions remain unanswered concerning the evaporation behaviors of multi-component fuels.
Motivated by the above analysis, this study scrutinizes the evaporation of n-hexadecane droplets and liquid films in subcritical and supercritical conditions by MD simulations.The objectives of this work are to study the mode transition mechanism of evaporation for hydrocarbon fuels and propose a criterion to quantitatively identify the mode transition of evaporation to fill the gaps of the CFD modeling and macroscopic experiments.In Section 2 , the research methods are described.The results and discussion are presented in Section 3 .The evaporation mode transition of a threecomponent droplet and a six-component droplet is studied in detail.A critical dimensionless number τ 0.9P , based on the average displacement increment (ADI) of fuel atoms, is proposed to identify the mode transition of evaporation of fuels.Based on this criterion, the mode regime map of n-hexadecane droplets is obtained.The mixing characteristics of fuels in different evaporation modes are discussed.Relative velocity effects are also investigated.Fuel droplets and liquid films with varied initial sizes are simulated by MD to understand the size effect.Finally, conclusions are drawn in Section 4 , with a focus on multi-component effects.

Interatomic potentials
All MD simulations in this paper were performed using the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS).Potential functions are the basis of MD simulations and describe all intermolecular interactions [34] .As mentioned, the UAM's treatment of each pair of CH 3 and CH 2 as one site, which keeps a good balance between computational accuracy and efficiency [10 , 25 , 26] , is widely used for large molecule systems, such as hydrocarbons.According to previous reports [36 , 37] , compared with other UAMs [12 , 38-40 ], such as the Optimized Potentials for Liquid Simulations United Atom (OPLS-UA) [12] model, the Transferable Potentials for Phase Equilibria United Atom (TraPPE-UA) model is more suitable for the simulation of long-chain alkane molecules and has been widely used in previous studies [10 , 20 , 24 , 26 , 41] .This study focuses on fuel molecules with longer chain lengths, so the TraPPE-UA model was chosen for simulations here.The 12-6 Lennard-Jones (LJ) potential was chosen to describe the non-bonded interaction between pseudo-atoms separated by more than three bonds or belonging to two different molecules: where U LJ denotes the non-bonded LJ interaction potential, r i j is the distance between the two interacting pseudo-atoms, ε i j is the energy parameter and σ i j is the size parameter.The threecomponent fuel investigated [26] [20 , 26] .All the LJ parameters, the bondstretching, bond bending and bond torsion parameters for alkanes simulated here can be found in authors' previous publications [20 , 26] .The LJ parameters for all unlike atoms were determined by the standard Lorentz-Berthelot combining rules [10] , ε i j = ε ii • ε j j and σ i j = ( σ ii + σ j j ) / 2 .Aromatic molecules were treated as rigid bodies [20] .In addition, the Coulombic pairwise interaction between charged molecules was described by: U Coul r i j = 1 4 πε 0 q i q j r i j (2) where U Coul represents the Coulombic pairwise interaction, ε 0 denotes the permittivity of free space.q i and q j denote the charges on the interactional charged atoms.Referring to previous studies [20] , the coulomb interactions of fuel molecules and linear diatomic molecules were ignored in the simulation.The potential energy parameters of molecules of ambient gasses are listed in The experiments by Crua et al. [2] are representative work in the field of supercritical evaporation of fuels.The fuel in their experiments evaporated in the combustion exhaust gas [2] .However, nitrogen is the most common background gas for fuel evaporations in MD simulations, including this study.In order to investigate the difference of evaporation of fuels in the nitrogen and the combustion exhaust gas, the evaporation of n-hexadecane in the combustion exhaust gas was also studied.Referring to the experiments by Crua et al. [2] , the composition of the combustion exhaust gas was determined by the complete combustion of a mixture consisting of acetylene (C 2 H 2 ), oxygen (O 2 ) and nitrogen (N 2 ) in certain ratios: where X is the dilution ratio.For convenience, nitrogen and combustion exhaust gas were marked as "A1" and "A2", respectively.When X = 8 in Eq. ( 3) , the composition of A2 was determined by the exhaust gas produced by the combustion of C 2 H 2 .Thus, A2 was composed of 87 mol% nitrogen, 8.7 mol% carbon dioxide, and 4.3 mol% water, similar to the composition of the combustion exhaust gas in the experiment [2] (the mole fractions of the three gasses were 89.7 mol%, 6.5 mol%, and 3.8 mol%, respectively).
In an earlier publication of the authors [26] , the potential energy function of n-hexadecane was verified by comparisons of VLE data and evaporation rate constants with experiments [44 , 45] .Moreover, the present potential models for multi-component fuels have also been validated by the distillation calculation and evaporation curves of fuels [20 , 26] .Results [20 , 26] have been shown to be in good agreement with experimental data [46 , 47] .

Simulation configurations
In this study, twelve initial configurations of fuel vaporization were constructed, including ten initial configurations for nhexadecane, one for the three-component droplet and one for the six-component droplet.The initial configuration of the nhexadecane evaporation system is shown in Fig. 1 , which includes droplets and films.Spray combustion is widely used in highpressure combustion devices such as diesel engines and rocket en- gines [10 , 17 , 45] .Many researchers in this field have developed fundamental theories about fuel droplet evaporation [3][4][5][6] , treating it as a single entity or contained in the spray [17] .It is of great significance to study the mode transition of evaporation of droplets in such combustion applications [10] .However, studying the configuration of liquid films is also necessary for research in this field.In the MD study by Qiao et al. [14] , the evaporation of the liquid film is size-independent.The configuration of the liquid film seems to be more suitable for studying the dynamics of the gasliquid interface [12 , 48] .Besides, the configuration of the liquid film is also widely used to verify intermolecular potentials in MD by comparing VLE data obtained with experiments [12 , 14 , 26] .Combining studies of droplets and liquid films will be helpful for answering the size effect of the transition of evaporation modes and bridging MD to practical combustion applications.
In the quiescent case for a droplet shown in Fig. 1 a, a single suspended n-hexadecane droplet was located in the center of a cubic simulation box, with a nitrogen ambient surrounding it.Before the simulations of fuel evaporation, the fuel droplet and the nitrogen ambient were simulated separately to respectively reach their own thermodynamic equilibrium state using the canonical ensemble (NVT, which means constant atom number N, constant volume V, and constant temperature T).After that, the two were combined together for the simulations of fuel evaporation.The initial diameter for the n-hexadecane droplet was 26.0 nm.The boundary of the droplet or the film was defined as the contour surface where the fluid density was equal to the average of the maximum and the minimum densities of the evaporation system.The diameter of the droplet was defined as that of a sphere of the same volume as the droplet.The evaporation lifetime was defined based on the temporal variation of the square diameter (D 2 ) of the droplet [10] .
To avoid the overlap of molecules, the nitrogen molecules which initially were located in the same region with the droplet in the simulation box were deleted.The initial temperature of the fuel for all cases was set to 363 K, which was near that of the fuel in real engines before injection [10] .The size of the simulation domain was 80 nm × 80 nm × 80 nm .Periodic boundary conditions were applied in all three dimensions.The simulations for fuel evaporation were conducted using the micro-canonical ensemble (NVE, which means constant atom number N, constant volume V, and constant energy E).The region outside of the sphere with a diameter of 79.8 nm was named "thermostat region", where molecular velocities were rescaled every time step using a speed reset method [20] .A fuel molecule would be deleted when it reached the thermostat region, for simulating the evaporation of the fuel taking place in an infinite space [10] .
Figure 1 b shows the configuration of the non-quiescent case.The initializations of quiescent and non-quiescent cases were the same except few details.To study relative velocity effects, in the non-quiescent configuration, the center of mass of the droplet was set to move along the x-direction at a constant speed.However, the ambient gas had no macroscopic translational speed.Correspondingly, the thermostat region was different.As shown in Fig. 1 b, the region for evaporation was a rectangular block with a size of 80 nm × 57.68 nm × 57.68 nm to keep the volume of thermostat region in this configuration the same as that in the quiescent case.Referring to the experiments by Crua et al. [2] , three initial relative velocities of droplets were investigated, which were 0 m/s, 30 m/s and 100 m/s, labeled "V0","V1" and "V2" respectively for convenience.It should be noted that, in the simulations here, referring to the previous study [23 , 24 , 49] , the droplet is given desired relative velocity initially.The ambient gas is stationary, that is, it has no initial bulk velocity.The relative velocity between the droplet and the ambient gas will decrease slowly [49] .Referring to previous MD studies on moving droplets [24 , 49 , 50] , the evaporation simulations here also use the NVE ensemble in the nonquiescent case.The kinetic energy of the bulk flow is also included in the energy term.Figure 1 c shows the initial evaporation configuration of the liquid film.The initialization of this configuration was similar to that of the droplet in Fig. 1 a.The size of the cuboid simulation box was 15 nm × 200 nm × 15 nm.The center of mass of the n-hexadecane liquid film was located in the center of the simulation domain, and the liquid film was surrounded by nitrogen molecules.Two separated thermostat regions with a size of 15 nm × 40 nm × 15 nm were located in two sides of the simulation box in the y direction.Similarly, fuel molecules would be removed when they reached the thermostat region.Periodic boundary conditions were used in all three directions.As a result, the current configuration could simulate the evaporation of an infinitearea fuel film in an open space [41] .Different initial configuration settings, such as the size of a simulation domain and the thermostat setting, were tested before the formal simulations presented.Five replica runs were performed to reduce or eliminate statistical errors.For n-hexadecane, seven initial configurations of droplets or liquid films with different initial sizes were constructed to study the size effect of MD, as detailed in Table 2 .It is worth mentioning that the initial size for droplets is the diameter and that for liquid films is the thickness.It should be noted that MD is a first-principle method that describes the fundamental interactions between atoms/molecules, without any assumptions about the system's macroscopic behaviours.Therefore, the pros of MD simulations are that MD results (1) are more reliable than conventional CFD (with prior assumptions of behaviours such as subcritical or supercritical evaporation), and (2) can reveal dynamic evaporation behaviours and provide quantitative information that are difficult to obtain by experiments.The cons of MD simulations are the small droplet sizes attainable.It is thus important to investigate what evaporation properties are size-dependent and what properties are size-independent.And if some properties are sizedependent, can scaling laws be obtained to bridge the size differences?This is an ongoing investigation but some conclusions can already be drawn.For example, the regression rates [22 , 23] and squares of droplet diameters [10 , 26] are size-dependent, but with the method of scaling, MD results can achieve satisfactory agreement with experimental values.More comparisons of distillation curves and evaporation rate constants from MD simulations with experiments could be found in [10 , 26] .As reported [10 , 26] , although the evaporation rate is size-dependent [22 , 23] , the evaporation rate constant obtained by MD simulations can agree with the macroscopic experiments [45] under subcritical conditions.Fortunately, some properties are size-independent [17] , and MD re-sults can be directly mapped onto macroscopic systems.More discussions about the relevance of MD results on evaporation to those in real combustion chambers, and the advantages and limitations of MD for investigating mode transition of evaporation could be found in the Introduction section.
The temperature and pressure for the various cases are initial values.As mentioned, the simulations for fuel evaporation were performed using the NVE ensemble.The velocities of the molecules located in the thermostat region were rescaled every time step [20] according to where E k indicates the total kinetic energy of the N t atoms in the thermostat region, and T indicates the target ambient temperature.As a consequence, the temperature of this region could be kept at a constant target value [10] .The calculations of temperature in MD simulations are based on the average kinetic energy of molecules [34] : where N is the number of atoms in the target atom group, m i and v i are respectively the mass and velocity of atom i .i indicates the identifier of the atoms.In the case with relative velocities, the temperature of a group of atoms is calculated after subtracting out the center-of-mass velocity of the group caused by the bulk movement, temperature only depends on the random portion of the molecular velocities, that is molecular thermal velocity [23] .
After the center-of-mass velocity has been subtracted from each atom, temperature is calculated by Eq. ( 5) .The number of atoms in the simulation domain varies with time due to evaporation, and the number of atoms N and their associated degrees of freedom in the calculation group are recalculated each time the temperature is computed to ensure that the temperature is properly normalized.The removal of the center-of-mass velocity is essentially computing the temperature after a "bias" has been removed from the velocity of the atoms.Moreover, in the thermostat region, this bias will be subtracted from each atom, thermostatting of the remaining thermal velocity will be performed, and the bias will be added back in to keep a constant temperature value.The initial ambient pressure was determined by a combination of the initial ambient temperature and the initial number of nitrogen molecules in the simulation box in NVT simulations.However, the ambient pressure was determined by choosing a suitable initial number of nitrogen molecules in the system in NVE simulations [10 , 20] .The number of nitrogen molecules ranges from tens of thousands to millions in different cases, which depends on the ambient condition and simulation domain size [10 , 20] .As mentioned, a fuel molecule will be removed when it gets to the thermostat region to simulate the evaporation of the droplet taking place in an infinite space [10] .It is worth mentioning that the droplet accounts for less than 4% of the volume of the simulation box, so the fluctuations of ambient pressure could be negligible [10 , 20 , 26] .For the case of liquid films, the conditions are similar [41] .
Figure 2 shows the initial configurations of multicomponent droplets.The detailed size parameters, such as the size of thermostat regions, have been marked on Fig. 2 .The initialization of configurations of multi-component droplets was the same as that of the single-component droplet in Fig. 1 a.The critical pressure and temperature of the fuel have important effects on its evaporation.For n-hexadecane, the three-component fuel and the sixcomponent fuel studied here, the estimated critical temperatures are 1.4 MPa, 1.89 MPa and 2.45 MPa, and the estimated critical temperatures are 722 K, 697 K and 686 K, respectively [20 , 26] .Taking future operating conditions of rocket engines into account [14] , the simulated ambient pressure and ambient temperature for fuels were in the range of 4-36 MPa and 750-3600 K, respectively, which were much higher than those investigated in similar reports [10 , 12 , 14 , 20 , 26] .Moreover, evaporation of droplets [10 , 17 , 20 , 26 , 31] or liquid films [12 , 14 , 34] was investigated in isolation in previous similar studies.In fact, those studies have not yet reached a unified conclusion on the size effect [10 , 14 , 31] , leaving unanswered questions on the mode transition of evaporation.This paper studies droplets and liquid films simultaneously, and proposes a size factor in Section 3.4 to unify the findings of the two configurations, which provides new atomic-level insights into the size effect of evaporation and further bridges MD to actual combustion applications.The time step was 1.0 femtosecond for "S" cases shown in Table 2 .The time step was 2.0 femtosecond in other cases.The total time steps for a case ranged from 10 0 0,0 0 0 to 960 0,0 0 0 determined by evaporation lifetimes.

Mode transition mechanism of evaporation for hydrocarbon fuels
It is well known that ambient temperature and pressure have a significant impact on evaporation modes of fuels.This paper investigated the profiles of the average displacement (AD) per fuel atom and ADI per fuel atom with time under varied ambient conditions.The results of configuration M were discussed here as an example.In case M, the diameter of the droplet was 26.0 nm.For comparisons, the displacement and its increment of fuel atoms were normalized here.The mixing time of fuels was also normalized.The dimensionless AD and ADI were obtained via dividing the displacement and its increment at each moment by the radius of the droplet (13.0 nm).The dimensionless mixing time, indicated by τ , was obtained via dividing the mixing time by the evaporation lifetime.The equations of AD and ADI are as follows: (6) where t n represents the present moment and t n −1 represents the previous moment.The sampling time interval ( t n − t n −1 ) is 0.05 ns.k represents the identifier of the fuel atoms and m represents the maximum identifier of the fuel atoms (the total number of fuel atoms) remaining in the simulation domain at time t n .AD ( t n ) represents AD at time t n , and the coordinate of displacement vector of the fuel atom k at time As stated in Eq. ( 6) , for fuel atom k , the reference position for calculating the displacement magnitude here is its initial position before the evaporation.In other words, the initial coordinate of the displacement vector of fuel atom k is (0, 0, 0).ADI ( t n ) represents ADI at time t n .R 0 is the initial droplet radius.AD R 0 ( t n ) represents the dimensionless AD at time t n .AD I R 0 ( t n ) represents the dimensionless ADI at time t n .
Figure 3 reveals the AD and ADI profiles of fuel atoms under different ambient conditions.As shown in Fig. 3 a, at low ambient pressures and temperatures (8 MPa@750 K, case1), the AD per fuel atom gradually increased with time, and the time-dependent profile exhibited a concave shape, meaning that the AD increased faster and faster with time.Correspondingly, in Fig. 3 b, the ADI per fuel atom gradually increased with time until near the end of the mixing process.However, at high pressures and temperatures (24 MPa@1350 K, case2), the time-dependent profile of AD was concave initially and then convex, meaning that the profile of growth rate of AD with time was unimodal, which was proved in Fig. 3 b.Obviously, two different mechanisms dominated evaporation in the case1 and case2 respectively.
In case1, the attractive forces of fuel molecules were always strong and only a small amount of ambient gas dissolved into the droplets via the droplet surface [20] .For the fuel, the vapor and the liquid phases coexisted, the interaction energy of the fuel molecules was dominant in this case, and the distribution of the fuel molecules was ordered.With continuous gaseous heat and mass transfer, the fuel was transformed from liquid phase to gas phase via the sharp phase interface in evaporation mode.It is worth mentioning the behaviors of liquid phase dominated the profile of ADI, because the liquid phase was much denser than the gas phase.After that, the evaporated fuel molecules via the vaporliquid interface diffused into the ambient gas.
In case 2, since it took time for heat and mass transfer, the surface tension of the droplet still played an important role initially, the liquid fuel initially transitioned into the gas phase in evaporation mode, and the ADI of the fuel atoms increased with time.However, with the dramatic transfer of heat and mass under high ambient temperatures and pressures, a large number of ambient gas molecules dissolved into the fuel droplets, the fuel density decreased rapidly [20] .Meanwhile, the attractive force between fuel molecules was greatly reduced, the surface tension declined rapidly until it disappeared.The fuel firstly was transformed from the liquid phase to the liquid-like phase, further to the gaslike phase, and finally to the gas phase [31 , 33] .In this process, the molecular disorder rose due to the difference in concentration of the species and the diffusion mode was dominant.For the fuel, the vapor and liquid phases could not coexist, the effect of entropy was decisive.The detailed discussion of entropy will be carried out based on Fig. 8 .Once the vapor-liquid interface disappeared, the mode transition of evaporation was fully completed.As the concentration difference of species decreased with time, the ADI of the fuel atoms gradually decreased with time after reaching a peak.
It is crucial to achieve a quantitative distinction between the two evaporation modes for the experiment and modeling of supercritical evaporation.Considering the fluctuation of MD results, a dimensionless number τ 0.9P was used here which was the smaller dimensionless mixing time corresponding to 90% of the peak of the time-dependent ADI profile for fuel atoms.If τ 0.9P is less than 0.5, it is defined as the diffusion mode.If not, it is in evaporation mode, as shown in Fig. 3 b.
As mentioned, the size of the droplets studied here is only tens of nanometers, and the corresponding evaporation lifetime is only a few to tens of nanoseconds, which is similar to previous MD studies [10 , 14] .Considering scale differences, the corresponding time in MD simulations is usually normalized by, for example, the evaporation lifetime [10 , 12 , 14 , 24 , 51] and a scale factor [17 , 21 , 22] .Kaltz [23] used the Long-Micci method to scale the time in MD simulations and investigated the possibility of extrapolating MD results to macroscopic systems.In this context, such timescale analyses used are valid and suitable for MD simulations.Referring to previous studies [10 , 12 , 14 , 24] , the time in this paper is normalized by the evaporation lifetime, and the criterion τ 0.9P here is based on the dimensionless normalized time.
Figure 4 illustrates the mode transition map on the P-T diagram for n-hexadecane droplets (case M) in the nitrogen ambient.As mentioned, MD simulations have inherent randomness.To overcome this, we have 5 replica runs with different initial configurations for each operating condition and the presented results are based on the averages of multiple runs [24] .Referring to previous MD work [24] , for each presented case, the repeatability of results is good and the statistical deviations are small.For example, for the case of evaporation under 16 MPa and 900 K in A1, τ 0 .9P = 0 .49 ± 0 .026 based on 5 runs.A fitted transition line of τ 0.9P = 0.5 (indicated by the red solid line) was obtained based on simulation points (T: 750-1700 K, P: 6-20 MPa), and the corresponding regression model using the least squares method was as follows: R 2 for this fitting model was 0.9708.As the ambient temperature or pressure increases, τ 0.9P gradually decreases, and the dominant mode gradually transitions from classic evaporation to diffusion.Therefore, the operating conditions above the transition line are dominated by diffusion, and those below are dominated by evaporation.The results here were compared with those from others, as shown in Fig. 4 .The experimental data in Fig. 4 were from Crua et al. [2] .Based on the morphological evolution of microscopic droplets in optical experiments, Crua et al. [2] provided the criteria to clarify the mixing regimes of fuels.In the experiments, the droplets in the regime of classical evaporation followed the classical evaporation process, and the surface tension played a dominant role [2] .For the droplets in the regime of transitional mixing, initially they followed the classic evaporation.But then due to the rapid reduction of intermolecular forces, the apparent deformation of the droplet was observed, and the mixing of the droplet was accelerated [2] .As for the regime of diffusion mixing, the evaporation of droplets initially was controlled by the surface tension.But soon, the droplets lost spherical shapes, and the classical two-phase flow transitioned to the single-phase mixing [2] .The regression models were obtained to distinguish the different mixing regimes [2] based on the experimental data using the least squares method, as shown in Fig. 4 .The transition line colored by blue was redrawn based on the MD simulations by Qiao et al. [14] .
In their study [14] , based on the critical mixing point from the VLE theory, they proposed a threshold dimensionless transition time of 0.35 to separate the subcritical and supercritical regimes, which was also plotted in Fig. 4 .
As shown in Fig. 4 , the transition line with τ 0.9P = 0.5 was generally in good agreement with the diffusion mixing line obtained by experiments [2] .The following discussion focused on the dif-ferences in the results of several methods.From Fig. 4 , the transition line obtained here can satisfactorily distinguish the experimental data [2] between transitional mixing and diffusive mixing, with acceptable differences.Specifically, the transition line for diffusion here was slightly higher than that obtained by the experiment [2] .In other words, the minimum ambient temperature for transition at constant ambient pressure was slightly higher or the minimum ambient pressure for transition at constant ambient temperature was slightly higher in this study.This consistency indicates that although there are huge scale differences between experiments and MD simulations, atomistic simulations can reveal fundamental mechanisms that also exist in macroscopic phenomena and complement macroscopic experiments, which is consistent with the results on multi-component fuels [24] .The possible reasons were as follows.First, the following discussion of the effect of relative velocity (in Fig. 9 ) showed that when a moving droplet evaporated in the ambient gas, compared to the quiescent case considered here, the minimum ambient temperature and pressure for transition became smaller.The theoretical analyses by Poursadegh et al. [16] led to similar conclusions.Second, the transition line shown in Fig. 4 was obtained by the evaporation of quiescent n-hexadecane droplets with a diameter of 26.0 nm in the nitrogen (for case M, A1(V0)).In optical experiments [2] , however, much larger micron-sized droplets were observed.According to the discussion of the size effect below (in Fig. 10 ), with larger size of the droplet, SSA H decreased, and the minimum ambient temperature and ambient pressure for the transition tended to decrease.The studies by Xiao et al. [10] and Wang et al. [31] also provided support for this conclusion.Third, the ambient gas was the nitrogen here, while it was the combustion exhaust gas in the experiment [2] .According to the following discussion on the multicomponent effects (in Fig. 11 ), compared with the nitrogen, the transition of the evaporation mode in the combustion exhaust gas was a little easier, which meant a smaller minimum ambient temperature or pressure for transition.However, judging from the difference in results, based on the ADI criterion, the results for nanoscale droplets slightly larger than here will be in better agreement with the experimental values, although size, relative velocity and multicomponent effects still exist.
A discussion about the limitations of experimental results themselves is instructive.Crua et al. [2] investigated the injection of three single-component fuels into high-temperature and highpressure ambient gasses using high-speed long-range microscopy.However, such optical experiments have intrinsic limitations [2] , such as optical effects under high-pressure conditions, the response of imaging system, and spatial resolution, and it is impossible to investigate actual states of interfacial structures [52] .In addition, such qualitative observations do not provide quantitative details such as droplet temperature profiles, and do not directly reveal the physics of mode transition of evaporation [10] .Crua et al. [2] addressed some of these limitations by using high-resolution microscopy at the University of Brighton [53] and high-speed microscopy at Sandia National Laboratories [15 , 52] .And, three criteria were proposed to account for optical uncertainty when classifying their experimental results [2] .With these limitations in mind, experiments can provide independent reference data for comparison with numerical simulations [2] , such as the MD studies by Qiao et al. [14] .
Moreover, the transition line colored by blue from Qiao et al. [14] was closer to the experimental transitional mixing line [2] , however, the transition line obtained here was closer to the experimental diffusion mixing line [2] .This may be due to the following reasons.The transition line of Qiao et al. [14] was based on the MD simulations of liquid films, while that in this paper was based on those of droplets.According to the discussion below on the size effect (in Fig. 10 ), SSA H of the liquid film studied by Qiao et al. varied from 0.05 (1/20) to 0.1 (1/10), which was smaller than that of the droplet studied here (0.231).Nonetheless, Qiao et al. [14] concluded that the transition line based on liquid films was size-independent.
In fact, previous studies have proposed the criteria for transition from evaporation to diffusion for liquid fuels.The critical locus criterion has been widely used, relying on the calculation of the critical point of a mixture composed of fuels and ambient gasses, whose base is the VLE theory [54] .Moreover, this approach relies on numerous thermodynamic assumptions and estimations, and fails to resolve complex evaporation systems with multi-components, because the critical points of them are hard to be available [11] .As mentioned, the Knudsen number criterion by Dahms et al. [13] is also proposed to distinguish the mixing regime of fuels.Nevertheless, the transition from classic evaporation to diffusive mixing is assumed to be in steady state in this theory, failing to resolve the dynamics of mode transition under supercritical conditions [14] .Qiao et al. [14 , 34] have tried to indicate the subcritical-to-supercritical transitions based on the calculation of surface tension in MD simulations.However, it has been hard and even impossible to achieve reliable calculation of surface tension under supercritical conditions [55] .In summary, a reliable and versatile criterion is urgently needed for supercritical and complex multi-component conditions to distinguish the evaporation mode of liquid fuels under certain conditions.The criterion of τ 0.9P here is a new attempt to achieve this goal.It should be noted that this criterion is based on the profile of ADI of fuel atoms, which is available for the fuel-ambient gas mixing system composed of any type and number of components, even under highly supercritical conditions.

Mixing characteristics of fuel droplets in different evaporation modes
Figure 5 reveals typical snapshots of the molecular distributions of n-hexadecane droplets in different evaporation modes.In Fig. 5 a (the non-transition case), the droplet always had a distinct vaporliquid interface and remained spherical during most of the mixing time.However, as shown in Fig. 5 b, at higher ambient temperatures and pressures, the mode transition of evaporation occurred.Before the transition ( Fig. 5 b1), the droplets remained spherical with a sharp interface.When mode transition occurred ( Fig. 5 b2), the spherical boundary was already difficult to recognize.After the transition ( Fig. 5 b3), the droplet was deformed and the original vapor-liquid interface was replaced by an expanding mixing layer full of fuel and nitrogen molecules with densities of continuous changes.
Figure 6 reveals the evaporation histories of n-hexadecane droplets in different evaporation modes.In evaporation mode (4 MPa@900 K), the profile of fuel molecular number in the droplet with time was first convex and then concave.The evaporation rate of fuel molecules gradually increased at first, and then decreased continuously in the second half of evaporation.In diffusion mode (16 MPa@1500 K), however, the profile of that with time generally was concave.The evaporation rate of fuel molecules reached the maximum value soon after the start of mixing, and then decreased continuously.Figure 6 b shows the histories of the number of newly evaporated n-hexadecane molecules leaving from the droplet, corresponding to the evaporation rate defined by the number of fuel molecules in Fig. 6 a.As shown in Fig. 6 b, in evaporation mode, the maximum evaporation rate occurred in the second half of evaporation.In contrast, in diffusion mode, the evaporation rate reached a peak soon after the start of mixing.In evaporation mode, the density gradient near the droplet interface was larger.Initially, due to the large temperature difference between the droplet and the ambient gas, the droplet temperature increased rapidly with heat  24MPa@1350 K).Here τ stands for the dimensionless mixing time.In Figure 5a, the mode transition does not occur.Figure 5(b1), (b2) and (b3) show the situation before transition, transition and after transition, respectively.Ambient nitrogen molecules are not shown here.transferring from the gas phase to the liquid phase via the interface [10] .Meanwhile, the evaporation rate became higher, during which the droplet experienced an initial heat-up period [26 , 34] .With time, the temperature difference between the droplet and the ambient gas decreased, and the droplet became smaller.The evaporation rate of the droplet decreased after reaching a peak.In diffusion mode, with higher ambient pressures and temperatures, a lot of nitrogen molecules dissolved into the droplet [20] , the density gradient between the liquid core and the ambient gas decreased [14] , and the evaporation rate peaked rapidly and then decreased.Meanwhile, the evaporation rate was determined by the combined effect of diffusion coefficient, mass fraction gradient and temperature gradient [14] .The temporal variation of the square diameter D 2 of the n-hexadecane droplet in evaporation mode (4 MPa@900 K) is shown in Fig. 7 .It can be found that the evaporation of droplets follows the classical D 2 law after an initial heat-up period, which is consistent with previous studies [10 , 26] .
The free energy, a function of internal energy, temperature and entropy, is a crucial parameter for measuring the state of a thermodynamic system.Specifically, the competition between energy and  entropy determines the state of a thermodynamic system, while the temperature determines the relative weight between them.The time-dependent profiles of the mean molecular kinetic energy, pair potential energy, and intermolecular interactions for fuel evaporation system like here have been presented in our earlier publication [20] and will not be discussed.The focus here is on the history of entropy of evaporation systems.The pair entropy fingerprint for each atom was calculated to investigate the degree of disorder of the mixing system.The calculation method of the entropy fingerprint of each atom was derived from Piaggi et al. [56] and Nettleton et al. [57] .It should be noted that the pair entropy, which is always negative, uses k B as its unit.The lower pair entropy stands for the more ordered environment.Figure 8 shows the histories of the sum of pair entropy for each atom and that of its corresponding increment.The mixture was composed of n-hexadecane and nitrogen.As shown in Fig. 8 a, although there were some differences, the profile of entropy change of the mixture was basically the same as that of the alkane, which meant the entropy change of the alkane was dominant.There were significant differences in the profiles of entropy for evaporation and diffusion modes.In evaporation mode (4 MPa@900 K), the time-dependent profile of entropy of alkanes nearly always was concave and eventually became convex.The entropy of alkanes increased gradually with mixing, whose growth rate rose to a peak near the end of evaporation and then declined, which was indicated by the increment profile in Fig. 8 b.In diffusion mode (16 MPa@1500 K), however, the timedependent profile of entropy of alkanes nearly always was convex except concave during a short period initially, whose growth rate rose to a peak in the first half of evaporation and then continuously dropped.Special attention should be drawn to the fact that the physics of the entropy change for the mode transition was consistent with that of evaporation histories shown in Fig. 6 and that of AD and ADI in Fig. 3 .This consistency provides a theoretical support and a further evidence for the ADI criterion.In addition, from Fig. 8 a, the pair entropy for alkane atoms in diffusion mode was greater than that in evaporation mode at almost all times.In diffusion mode, the ambient temperature and pressure were high enough relative to the critical point of the fuel.Meanwhile, the intermolecular interactions of fuel molecules declined dramatically [20] , the order of molecules of liquid fuel was hardly maintained, and the vapor-liquid interface no longer existed.The entropy which caused molecular disorders, instead of energy which maintained molecular orders, became the dominant factor in the mixing process.As a result, the disorder of the fuel molecules increased significantly compared with that in the two-phase evaporation.
It should be mentioned that in Figs. 3 , 6 and 8 , the concavity or convexity of the curves corresponds to different patterns in the slope change of the curves.In fact, Figs. 3 , 6 and 8 illustrate the physics of different evaporation modes and the mode transition from the perspectives of AD, evaporation rate, and pair entropy, which reveal different aspects of the same physical process.Therefore, the physics revealed by these curves are consistent.It was reported [2 , 10 , 14] that even under very high ambient temperatures and pressures, the fuel droplet did not achieve the mode transition of evaporation instantly, and the surface tension did not disappear immediately, because heat and mass transfer took time.As mentioned here, the competition between energy and entropy determines the state of a thermodynamic system, which is also applicable for evaporation and mode transition, so the shape of the curve is determined by a combination effect of energy and entropy.Under this condition, the concavity and convexity of curves and their connection point (inflection point) reflect exactly which factor plays a dominant role at that time.In evaporation mode, the surface tension was relatively strong and the interaction energy of the liquid fuel was dominant compared with its entropy change.The order of fuel molecules in the liquid phase was more or less kept.With continuous heat and mass transfer, the evaporation was accelerated, which could be revealed by concavity of the profiles (The slope of profiles increased).In diffusion mode, however, massive ambient gas molecules entered into the fuel droplet and fuel density decreased dramatically [20] .The attractions between fuel molecules, which was responsible for the molecular order, dramatically declined.The original order of fuel molecules in the liquid phase was lost because the entropy change caused by the mixing became dominant.The evaporation was first accelerated and then decelerated due to the reduced density gradient, which could be revealed by the convexity of the profiles (The slope of profiles decreased).

Relative velocity effects
As mentioned, effects of relative velocity on the mode transition of evaporation of fuels have not been well understood.For this purpose, in addition to the quiescent case (V0), we investigated the evaporation of n-hexadecane droplets with initial relative velocities of 30 m/s (V1) and 100 m/s (V2) in the supercritical nitrogen.Figure 9 shows the relationship between fuel temperature and molar fraction of the fuel for different droplet relative velocities at an ambient pressure of 8 MPa.The temperature was calculated from the average kinetic energy of the molecules [14] .The estimated critical temperature and the critical molar fraction of fuels for the mixing system of n-hexadecane and nitrogen at 8 MPa were shown in Fig. 9 , whose calculation method could be found in the authors' previous report [26] .When the profile of temperature-molar fraction passes through or is above the critical point of the mixing system, supercritical evaporation of fuel can occur [14] .As shown in Fig. 9 , under an ambient pressure of 8 MPa, for the quiescent case, when the ambient temperature ranged from 750 K to 1150 K, the mode transition of evaporation did not occur.As the relative velocity increased, the temperature-molar fraction curve shifted toward the region with higher temperature and mole fraction of fuels in Fig. 9 , although this was a non-linear process, which meant the mode transition became easier.In particular, the mode transition occurred when droplets had relative velocities of 30 m/s and even higher 100 m/s.However, this did not happen in the quiescent case under the same ambient condition.This may be due to the promoted gaseous heat and mass transfer [29] .This finding provided a direct evidence for the results of Poursadegh et al. [16] by theoretical analyses.Further work is needed in this area.
Another discussion is about the possible slip effect.In a direct injection engine, such as a diesel engine, the droplets produced by the fuel spray are usually micron-sized [2] .In this condition, it is appropriate to use the continuum assumption and an approximation of no-slip condition is valid [58] .However, the evaporation of the nanoscale droplet here may conform to the slip regime based on the Knudsen number [59] , where the no-slip condition is invalid.As reported [60] , the slip velocity may affect heat and mass transfer and thus the evaporation rate.Therefore, it is necessary to evaluate the effect of slip velocity on MD results, especially for the cases with relative velocities.In contrast, the relative droplet velocity simulated here is smaller, similar to that in dilute regions of sprays in diesel engines, where the droplet-gas slip velocity is small [60] .Kaltz [23] studied the evaporation of nanoscale droplets with a relative velocity of 500 m/s under a high pressure environment of 20 MPa, and found that the slip between nanoscale droplets and ambient gasses was almost negligible.This relative velocity is much larger than the velocity studied in this paper, which means that the slip velocity in this paper is also negligible.Michaelides et al. [61] found that for moving nanoscale particles, a velocity slip at the interface did not significantly affect the overall Nusselt number.Masri et al. [62] found that droplets of diameter less than 10 μm in spray had similar velocities to the carrier flow.Similar findings have also been reported by Mastorakos et al. [63] .Nomura et al. [45] mentioned that the fuel droplets in the spray combustor were so fine that even in high pressure environments, the effect of convection on droplet evaporation was negligible.In fact, regardless of droplet sizes or slip velocities, ambient temperature and pressure play a dominant role in the mode transition of evaporation [2 , 10 , 14 , 24] .The study here supports the argument that slip has little effect on the mode transition of evaporation, which is supported by the fact that the mode transition map obtained by MD can agree satisfactorily with macroscopic experimental results, consistent with conclusions in [10 , 14 , 24] .On the other hand, MD is based on Newton's second law and does not make any additional thermodynamic assumptions about the simulated physical process other than interatomic potentials.The simulation results are physics-based rather than a priori, which means that even if the slip effect is present here, it will be naturally included in the results.Given that the slip velocities and droplet sizes in this paper are not large enough, future work is needed in this area, which is beyond the scope of the present study.

Size effects
To study the size effect, the evaporation of several fuel droplets and liquid films with different initial sizes was investigated.Referring to previous studies on size effects of MD [10 , 14 , 31] , the diameters of the droplets were selected as 20.4 nm, 26.0 nm,30.5 nm and 73.0 nm, respectively.The thicknesses of the liquid films were 16.9 nm, 25.1 nm and 40.0 nm, respectively.With varying ambient temperatures and pressures, the dimensionless transition time τ 0.9P for every case was calculated, as shown in Fig. 10 a.With increasing ambient temperature and pressure, τ 0.9P for fuels with different initial sizes decreased, and the mode transition of evaporation became easier.In addition, τ 0.9P of several sizes of fuels had the same size relationship for investigated conditions.SSA H of several initial sizes of fuels was studied, which was defined as the ratio of the surface area available for heat and mass transfer to the total volume of the fuel.When it comes to the size effect, the discussion is actually on the same fuel.From the MD point of view, the evaporation is derived from the movement of fuel molecules, and the atomic diameter is a key parameter, which could be indicated by the size parameter σ mentioned before.For fuel atoms, the size parameter σ is indicated by the weighted average of that of all the united atoms in fuel molecules.In this case, σ for the hexadecane is 3.925 nm [ 26 ].The dimensionless size factor S F is defined by multiplying SSA H by σ , which is shown in Fig. 10 b.In this condition, S F can be named as a normalized specific heat transfer surface area.It is well known that the evaporation of fuels is strongly dependent on the heat and mass transfer processes under certain conditions.Hence, there is no doubt that S F defined here will play an important role in this process.As shown in Fig. 10 , for the same ambient temperature and pressure, the smaller S F of the fuel, the easier the transition.A plausible reason for this was that the smaller S F of the fuel, the slower the evaporation of the fuel, the longer the evaporation lifetime, and the more sufficient time for the fuel to achieve the transition.In the studies of droplets by Xiao et al. [10] and Wang et al. [31] , it was believed that larger droplets were easier to transition because of their longer evaporation lifetime compared to their transition time.Similar conclusions have been reached in studies of droplets from continuum-based simulations [64] .However, in the study on liquid films by Qiao et al. [14] , the supercritical transition was size-independent.Note that their initial configurations for fuel evaporation were different.In the studies of Xiao et al. [10] and Wang et al. [31] , the fuel molecules would be removed when they reached the thermostat region, which was not the case in the study of Qiao et al. [14] .S F proposed here can give unified analyses for the results of droplets or films with different initial sizes.Although size effects exist, as shown in Fig. 4 , the results for medium-sized droplets (for case M) were already in good agreement with the experimental data [2] .Based on the above facts, al-  though the mode transition maps obtained by droplets or liquid films with different initial sizes were slightly different, the results of droplets or liquid films with sizes of tens of nanometers studied here were acceptable for practical applications.Moreover, in Ref. [10] , the droplet size effects on evaporation of single-component fuels have been investigated, and no qualitative differences are observed but there are gradual quantitative changes.Similar conclusions have been drawn in Ref. [31] .Further work is needed in this area, provided computational resources are made available.

Multicomponent effects
Although there have been some attempts [12 , 33 , 34] for the different multi-component evaporation systems, if subtle differences of the evaporation process are ignored, there has not yet a consensus concerning which factor plays a decisive role in the mode transition of evaporation.In other words, which factor is primarily responsible for the difference in the mode transition of evaporation for a mixing system composed of different components?Following the author's previous publications [20 , 26] , this work was a new attempt to answer this question.Four systems were set up to study the multi-component effect, that is, a sixcomponent fuel droplet (A1), a three-component droplet (A1), nhexadecane droplet (A1), n-hexadecane droplet (A2), respectively.As mentioned, for n-hexadecane, the three-component fuel and the six-component fuel studied here, the estimated critical temperatures are 1.4 MPa, 1.89 MPa and 2.45 MPa, and the estimated critical temperatures are 722 K, 697 K and 686 K, respectively [20 , 26] .Figure 11 illustrates τ 0.9P for several evaporation systems under the same ambient temperature and pressure condi- tions.When the ambient gas was nitrogen, in the order of sixcomponent fuel, three-component fuel and single-component fuel, τ 0.9P decreased in turn, and the mode transition of evaporation became easier.When the fuel was n-hexadecane, due to the ambient gas molecules changing from A1 to A2, τ 0.9P slightly decreased, and the transition also became easier, which is consistent with a previous study on multicomponent fuels [24] .
Figure 12 exhibits the mixing regime of Crua et al. [2] for several single-component fuels (n-heptane, n-dodecane and nhexadecane) which evaporated in the exhaust gas.The critical temperatures of n-heptane, n-dodecane and n-hexadecane were 540 K, 658 K and 722 K; and the critical pressures were 2.74 MPa, 1.82 MPa and 1.4 MPa, respectively.The reduced ambient temperature Tr and the reduced ambient pressure Pr were both calculated by dividing ambient values by critical values of fuel.Based on their experimental data, Crua et al. [2] proposed a regression model, whose expression was Tr * Pr. ^ 0.5, to unify the results of different single-com ponent fuels.If Tr * Pr. ^ 0.5 was lar ger, it meant that the mixing regime of the fuel tended to be the diffusive mixing.As shown in Fig. 12 , Tr * Pr ^ 0.5 of the three fuels under the same ambient temperature and pressure was calculated for investigating the possibility for mode transition of the three fuels under the same ambient conditions, using that of n-hexadecane as a reference.As shown, in the order of heptane, n-dodecane and n-hexadecane, Tr * Pr ^ 0.5 increased in turn, which meant that the mode transition became easier.In other words, for the three investigated fuels, under the same ambient temperature and pressure, n-hexadecane had the largest Tr * Pr. ^ 0.5, and was the easiest to achieve the diffusive mixing.Based on the above discussions, an interesting consistency between multi-component fuels and singlecomponent fuels was found.It is generally believed that the critical pressure and temperature of the fuel have important effects on its evaporation.Here, Pc1 and Pc2 indicate critical pressures of Fuel 1 and Fuel 2, respectively.And Tc1 and Tc2 indicate critical temperatures of Fuel 1 and Fuel 2, respectively.For the three fuels investigated (including two multi-component fuels and one singlecomponent fuel), when the fuel with a larger critical pressure is regarded as Fuel 1, if Pc1/Pc2 > Tc1/Tc2, then Fuel 2 is easier to achieve the mode transition of evaporation, compared with Fuel 1.This correlation implies that pressure seems to play a more important role than temperature for mode transition of evaporation.It is noted that the minimum ambient temperature (750 K) studied here is greater than the critical temperatures of three fuels investigated (722 K, 697 K and 686 K, respectively).As reported by Zhang et al. [12] , compared with the fuel composition, the ambient pressure and temperature were more crucial factors for determining the transition time.Qiao et al. [65] found that at an ambient pressure fixed to the critical pressure of n-heptane, when the ambient temperature was higher than the critical temperature of n-heptane, the required transition time from evaporation to diffusion decreased extremely slowly with increasing temperature.On the other hand, according to Gong et al. [26] , for single-and multi-component fuels, at an ambient temperature higher than the critical temperatures of fuels, when the ambient pressure was far higher than the critical pressure of the fuel, the fuel evaporation lifetime decreased very slowly with increasing pressure, while that decreased rapidly with increasing temperature.Similar results are reported by Xiao et al. [10] and He et al. [66] .The above discussions show that the evaporation lifetime has different sensitivities to the changes of pressures and temperatures.Since τ 0.9P is a dimensionless time via dividing the absolute time by the evaporation lifetime, it may lead to the fact that the effect of pressure on mode transition of evaporation indicated τ 0.9P is greater than that of temperature.As reported [10 , 14] , having the time normalized by the evaporation lifetime is a common method in this field.In fact, experimental results of the single-component fuels by Crua et al. [2] are also consistent with the above correlation.

Conclusions
The evaporation processes of n-hexadecane droplets in the nitrogen or exhaust gas ambient were studied using MD simulations.Moreover, three-component droplet and six-component droplet evaporation processes were also studied.The ambient pressure ranged from 4 MPa to 36 MPa and the ambient temperature ranged from 750 K to 3600 K, which covered ambient conditions for both subcritical and supercritical evaporations.Significant findings of this study included: 1) The profiles of the average displacement (AD) and average displacement increment (ADI) per fuel atom with time under varied ambient conditions were investigated in detail.At low ambient pressures and temperatures, the ADI per fuel atom gradually increased with time until near the end of the mixing process.However, at high pressures and high temperatures, the time-dependent profile of the ADI had a peak soon after the start of evaporation.Considering scale differences, the time in this paper is normalized by the evaporation lifetime.
A dimensionless number τ 0.9P was proposed here which was the smaller dimensionless mixing time corresponding to 90% of the peak of the time-dependent ADI profile for fuel atoms.A critical dimensionless number τ 0.9P of 0.5 was used to identify the mode transition of evaporation of fuels with any type and number of components.If τ 0.9P was less than 0.5, it indicates the diffusion mode.Otherwise, it signals the evaporation mode.The comparison between mode transition map from experiments and MD simulations revealed that based on the ADI criterion, the results from droplets with diameters of tens of nanometers could be in satisfactory agreement with those from experiments of micro-sized droplets.This consistency indicates that although there are huge scale differences between experiments and MD simulations, atomistic simulations can reveal fundamental mechanisms that also exist in macroscopic phenomena and complement macroscopic experiments.The molecular physics obtained in MD simulations will be helpful for understanding the mechanism of mode transition of evaporation.
In addition, the knowledge of droplet evaporation will be expanded by such MD simulations.
2) The pair entropy fingerprint for each atom was calculated to investigate the degree of disorder of the mixing system.The entropy change of the alkane dominated that of the mixture.In evaporation mode, the time-dependent profile of entropy of alkanes nearly always was concave and eventually became convex.The entropy of alkanes increased gradually with mixing, whose growth rate peaked near the end of evaporation and then declined.Nonetheless, in diffusion mode, the timedependent profile of entropy of alkanes nearly always was convex but concave during a short initial period, whose growth rate peaked before half of evaporation and then continuously dropped.The physics of the entropy change for the mode transition was consistent with that of evaporation histories and that of AD and ADI.The competition between energy and entropy determines the state of a thermodynamic system, which is also applicable for evaporation and mode transition, so the shape of such curves is determined by a combination effect of energy and entropy.Under this condition, the concavity and convexity of curves and their connection point (inflection point) reflect exactly which factor plays a dominant role at that time.
In the diffusion mode, the entropy became the dominant factor in the mixing process, and the disorder of the fuel molecules increased significantly compared with that in the two-phase evaporation.Moreover, MD can reveal dynamic evaporation behaviours and provide quantitative information that are difficult to obtain by experiments.For some size-dependent evaporation properties, it is promising to obtain scaling laws to bridge the size differences.For some size-independent evaporation properties, MD results can be directly mapped onto macroscopic systems.3) Compared with the case of a quiescent droplet, the mode transition became easier with increasing relative velocity between the droplet and the ambient gas, although this was a non-linear process.This was due to the promoted gaseous heat and mass transfer and stagnation effects.Moreover, the slip effect possibly exists during the evaporation of nanoscale droplets, especially for moving droplets.However, even if considering real spray combustors, such as diesel engines, it is estimated that the slip effect on evaporation is really limited and will be naturally included in the MD results.For the same ambient temperature and pressure, the smaller the normalized specific heat transfer surface area (S F ) of the fuel, the easier the mode transition.The possible explanation for this is that the smaller SSA H of the fuel, the slower the evaporation of the fuel, the longer the evaporation lifetime, and the more sufficient time for the fuel to achieve the transition.S F proposed here could give unified analyses for droplets or films with different initial sizes, which will be helpful for answering the size effect on the mode transition of evaporation and bridging MD to practical combustion applications.A correlation was proposed to compare the possibility of mode transition of evaporation for single-and multi-component fuels.When the ambient temperature exceeds the critical temperature of fuels, pressure seems to play a more important role than temperature for mode transition of evaporation of fuels.
Future work is needed in this area.Considering the large size difference between MD simulations and actual engineering applications, MD simulations with larger sizes should be performed, provided computational resources are made available.More accurate intermolecular potentials with lower computational costs should be developed further for complex hydrocarbons.Moreover, effects of relative velocity on mode transition of evaporation should also be further investigated in detail.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 1 .
Fig. 1.Initial configurations of the n-hexadecane droplet/film: (a) Quiescent n-hexadecane droplet (for case M), (b) Non-quiescent n-hexadecane droplet (for case M) and (c) n-Hexadecane film (for case F2).Purple particles indicate n-hexadecane molecules and yellow particles indicate nitrogen molecules.Molecules of ambient gasses surround the fuel, which are not shown here in (a) and (b).

Fig. 2 .
Fig. 2. Initial configurations of multicomponent fuel droplets (front view): (a) The three-component droplet and (b) The six-component droplet.Nitrogen molecules surround the droplet, which are not shown here.

Fig. 7 .
Fig. 7.The temporal variation of the square diameter D 2 of the n-hexadecane droplet.The dotted blue line indicates the fitting of D 2 profile [10] .

Fig. 8 .
Fig.8.Temporal variations of the sum of pair entropy and that of its increment for each atom (for case M, A1(V0)).

Fig. 9 .
Fig. 9. Temperature as a function of molar fraction of n-hexadecane in the case of different droplet relative velocities (for case M, A1).

Fig. 10 .
Fig. 10.Size effects on the mode transition of evaporation of fuels.Here D indicates the diameter for droplets or the thickness for films.

Table 2
Simulation details of n-hexadecane droplet/films.