Specificity Switching Pathways in Thermal and Mass Evaporation of Multicomponent Hydrocarbon Droplets: A Mesoscopic Observation

For well over one century, the Hertz–Knudsen equation has established the relationship between thermal – mass transfer coefficients through a liquid – vapour interface and evaporation rate. These coefficients, however, have been often separately estimated for one-component equilibrium systems and their simultaneous influences on evaporation rate of fuel droplets in multicomponent systems have yet to be investigated at the atomic level. Here we first apply atomistic simulation techniques and quantum/statistical mechanics methods to understand how thermal and mass evaporation effects are controlled kinetically/thermodynamically. We then present a new development of a hybrid method of quantum transition state theory/improved kinetic gas theory, for multicomponent hydrocarbon systems to investigate how concerted-distinct conformational changes of hydrocarbons at the interface affect the evaporation rate. The results of this work provide an important physical concept in fundamental understanding of atomistic pathways in topological interface transitions of chain molecules, resolving an open problem in kinetics of fuel droplets evaporation.

can model transient "internal rotations" during phase transition as introduced in ref. 14. The general expression developed in ref. 10 and applied in ref. 9 for estimation of β m is suggested to be valid if "isotropic" requirements at the interface are established 10 . In the vicinity of n-dodecane droplets surfaces the "anisotropic" effects have been suggested using ab initio molecular dynamics (AIMD) simulation 15 and dynamic reaction coordinate (DRC) analysis, however 16 . Classical MD simulation results were consistent with a general expression for estimation of β m 10 but reliability of GAFF 10 and OPLS force field 11 becomes particularly questionable at "interface" for molecules with multi-conformers. Support for this is to simulate surface tensions of organic molecules which have been calculated using GAFF and OPLS force fields within 10-20% of experimental values at room temperature 29,30 . Additionally, the aforementioned non-reactive FFs 10,11 and NERD force field 31 have predicted n-alkanes molecular orientation along the surface differently which are not in agreement with experimental measurement by vibrational sum frequency spectroscopy (VSFS) 32 . Therefore, this general expression for estimation of evaporation coefficient 10 is likely to fail for modeling conformational changes at the interfacial layers 14 .
In this article, thermal (β T ) and mass evaporation (β m ) coefficients and evaporation rate (γ) are respectively calculated using novel transient reactive molecular dynamics simulations, the statistical associating fluid theory (SAFT) and "quantum transition state theory/improved kinetic gas theory" (QTST/IKGT) 14 . We apply a non-equilibrium MD simulation technique using ReaxFF 33 and SAFT equation of state 34 with quantum chemical calculations 35 to fundamentally understand how interfacial flows in n-dodecane droplets affect β T and β m kinetically/thermodynamically. We also develop and present the results of a new version of QTST/IKGT for multicomponent hydrocarbon systems which have been inspired by the "discrete" Boltzmann method. These findings provide an important physical concept: dynamic coupling between liquid and gas phases during the evaporation, which should be useful in better understanding the simultaneous influences of thermal and mass transfer on the evaporation rate of multicomponent fuel droplets at the atomic level.

Results and Discussion
A non-equilibrium MD simulation technique using ReaxFF is proposed to determine thermal evaporation coefficient of n-dodecane, a representative of n-alkanes in Diesel fuel. The reliability of reactive force field ReaxFF has also been assessed and compared for modeling the evaporation of hydrocarbons 15 using the quantum , re-equilibration stages are presented at initial temperature 400 K. The temperature jump is observed at the interface during the re-equilibration stages due to cooling effects. chemical calculations (DFT and PM7 methods) and available experimental data on the determination of bond energies, Gibbs free energies of internal molecular dynamics of a set of n-dodecane conformers and collision energies of attacking molecules with the surface of the droplet. It was shown that ReaxFF performs better than semi-empirical quantum chemistry PM7 method in terms of both cost and accuracy of calculations of the evaporation of n-dodecane. Therefore, the bond energy bond order approach of ReaxFF is applied to study thermal effects induced over the interfacial flows during the evaporation process of n-dodecane. The aim of this simulation is to investigate whether the thermal coefficient values are temperature dependent similar to mass evaporation/condensation coefficient, while we examine the interfacial temperature discontinuities. The thermal evaporation coefficient is defined as: where T i , T g and T l refer to, respectively, the effective temperatures in the interfacial layer, gas and liquid phases assuming a semi-spherical droplet evaporates into vacuum without any movement (see Methods). The results of this simulation will give us important information for better understanding the energy transfer mechanisms in the initial transient stage of the evaporation process. The vacuum conditions have already been investigated experimentally, theoretically and computationally on simple fluids [36][37][38][39] . The time evolutions of various average molecular energies and corresponding temperatures were obtained at various stages of droplet heating and evaporation, identified as heating, perturbation, and re-equilibration. The values of temperatures are shown in Fig. 1. The initial droplet heating was set up to take place during 1000 fs using the Berendsen's thermostat 40 . During this period the average temperature of the droplet reached 400 K. At the   [9][10][11] , temperature gradient at interface could be modelled using the ReaxFF method in a three dimensional model.
later times some oscillations of droplet temperature were observed, with the liquid temperatures being almost always below the interface temperature. This stage is called the heating stage (see Fig. 1a). At the next stage, the system was perturbed using various coupling time constants (τ T ) as described in Methods. When the interface is strongly coupled to the thermostat using τ T = 1 fs and T 0 = 400 K, the temperature is controlled by the velocity rescaling algorithm used in the Berendsen's thermostat. But when τ T = 100 ps is specified for liquid phase, the temperature is maintained only through thermal effects induced by conformational changes since the thermostat does not function due to very slow rate of change of kinetic energy, therefore canonical ensemble (NVT) is essentially converted to the micro-canonical one (NVE). Application of these two different coupling time constants on a liquid drop leads to the control of liquid and interface temperatures in two different ways causing an oscillation of liquid phase temperature around the droplet temperature and showing transient transfer of energy mostly between liquid and interface since internal rotations or torsions do not change the centre of mass of molecules (see Fig. 1b). The averaged liquid temperature value was determined to be higher than interface one during the perturbation stage (T l = 399.83 K against to T i = 398.72 K). This discrepancy can be explained by the fact that molecules at the surface with high energy leave the drop leading to cooling effect caused at the interface. Moreover, as already mentioned our analysis in estimation of temperature has been done on molecules that stay in their relevant sub-systems (liquid or interface) during the evaporation process. After imposing these non-equilibrium conditions, the formation of some nano-bubbles of 1-3 nm in diameter was observed in the liquid phase. This illustrates how inversion of heat energy affects the structure of the liquid phase. These nano-bubbly flows into the droplet gradually disappeared when the sub-systems (gas, interface and liquid) reached the quasi-equilibrium state and molecules in the liquid phase could show expected behaviours again (compare structures in Fig. 2 after stages of b and c). As seen in Fig. 1c, the liquid phase has temperature higher than the droplet and even the interface during the re-equilibration in which system will be simulated using coupling τ T = 100 ps implying minimal perturbation effects caused by thermostat. This is related to the fact that the directions of transfer of heat and mass are not the same during evaporation leading to higher temperatures in the liquid phase relative to the interface. Gas temperature during the evaporation drops about an order of magnitude and reaches the saturation state as the energy transfers from the gas phase to the interface and then into the liquid phase in a stepwise manner. As shown in Table 1, values of the evaporation coefficient are identical at temperatures 350 and 400 K with a time constant of 2.3 × 10 −4 ps −1 and we can expect that those do not change dramatically at higher temperatures as well.
Mass evaporation coefficient can be derived in terms of thermodynamic potentials and SAFT molecular based equation of state 34 . SAFT can be applied for predicting interfacial layer thickness of fluids and it incorporates the effects of chain length, molecular association and other interactions such as long-range dipolar forces and dispersions. While the interfacial layer effects were not explicitly modelled in refs 9, 13, 15 and 16, we consider these effects in this study by setting up an equation including interfacial width, δ. A standard state has to be defined for the evaporation/condensation process and with this thickness the relationship between the free energy of evaporation/condensation (ΔG g↔int ) and the coefficient β m becomes: where <ΔG g↔int > presents the average difference values of Gibbs free energy of conformers in the interfacial layer and gas phase (see Methods). Taking the interfacial layer effects and relevant correction terms into consideration, the same results are obtained as reported in ref. 9. One question arises concerning whether or not adding the interfacial layer using SAFT has had no effects on evaporation/condensation coefficient of n-dodecane. The answer is no since SAFT, which is a thermodynamic-based approach, cannot model properly transient processes such as the internal rotations in chain molecules 14 . While this molecular theory can provide useful interfacial properties, it cannot describe the interface at an atomic level. More specifically, in all diffuse interface models the existence of interfacial width is inherent and once it reduces to a length scale which is small in comparison with the macroscopic length scale associated with the motion of the two bulk fluids, these models are related to the free-boundary problems 41 . We believe that these sorts of equations are fundamentally unable to track "thermal effects" induced with "transition states" over the interfacial layers during the evaporation process 42 . We do not think that classical diffuse interface models can capture "quasi-equilibrium" transition states and internal molecular dynamics effects in complex molecules which have multi-structural effects. The internal rotations in multi-conformers cannot be modelled based on classical and harmonic models and therefore anharmonicity effects (conformational changes and the coupling between torsions and vibrational modes) should be considered based on quantum mechanics theory and a suitable statistical mechanics method in which the atoms in molecules (AIM) motions are taken into account. In order to understand simultaneous relationship between thermal and mass evaporation with evaporation rate in multicomponent fuel droplet hydrocarbons, we have applied an extension of the quantum transition state theory/improved kinetic gas theory (QTST/IKGT) 14 . The evaporation flux is first predicted based on the assumption that single molecular events occurring during the evaporation of individual components from a multicomponent liquid phase are independently and identically distributed; and then we generalize the total solution evaporation flux for a c-component system as a summation of individual component evaporation fluxes which are in equilibrium in gas and liquid phases. These expressions can be easily applied to mixtures with any number of chemical components (and not just binary mixtures, as is the case for this study because of the available experimental data):  conformers i and A jk represents the gas/vapour molecules or clusters/droplets accessible surface area. ∆ ↔ G g l c presents the average difference values of Gibbs free energy of each component between liquid and gas phases.
is the activation Gibbs free energy induced by internal rotations in each conformer including zero-point energy. The m j and r j present the mass and radii of gas/vapor molecules or clusters/droplets colliding with other clusters/droplets and gas/vapour molecules with the mass m k and radii of r k .
We distinguish the quasi-equilibrium phenomena induced by the "internal rotations" dynamics relevant to thermal evaporation effects over the interfacial layers from equilibrium mass evaporation/condensation occurring between the gas and liquid phases. For the sake of simplicity, the model used in Fig. 3 includes two active site loops for two-component systems. Although, in reality there may be a large number of different conformers R i and P i (i = 1,2,3,….), we examine nano-confinement mechanistic hypothesis in which two conformers of each component are confined across the interface to be actively involved in phase transitions (see Fig. 3). We will also only consider the case of an ideal liquid mixture with incompressible liquid components and an ideal vapour mixture with each vapour component treated as an ideal gas. The expressions for Gibbs free energy of each component and their mixtures in liquid and vapour phases are given by: 43 Figure 4. Evaporation rates of a binary fuel droplet. The fits show that QTST/IKGT reproduces temperatureand pressure-dependent evaporation rate in binary fuel droplet with 1.2 mm diameter. The (un) circles and solid (dash) lines respectively represent experimental measurements and results obtained by our model -with the parameters given in Table 2. The fitted data present effects of (a) temperature and (b) pressure on evaporation rate of a mixture of 50% n-heptane and 50% n-hexadecane in liquid phase and at six different mole fractions in the gas phase.
where y Ci and x Ci refer to mole fractions of i th component C in the liquid and vapour phases, respectively. In this study these components are C 7 and C 16 . Substituting equations (7-10) in KGT-based equation (6), we obtain the final expression for predicting evaporation flux for each component in a c-component mixture. Equation 3 can also be rearranged to predict evaporation rates for multicomponent liquid mixtures. Although details of the model presented here are novel for better understanding of new mechanistic pathways in evaporation of multicomponent fuel droplets, we note that it shares similarities with previous models applied for kinetics modeling of mono-component hydrocarbon droplets evaporation 9,14 . We do not claim that QTST/IKGT based on this nano-confinement hypothesis should be taken as the correct kinetic model for each hydrocarbon. Nevertheless, the results in Fig. 4a,b show that QTST/IKGT establish temperature and pressure dependence of the evaporation rate of binary fuels, as long as two "equilibrium" conformational changes in liquid and gas phases are cyclically switched on and off with two other "quasi-equilibrium" transition states at the liquid-gas interface of each component (those are controlled by k 1 (k′ 1 ) and k 2 (k′ 2 ) -see Fig. 3). We also note that neglecting conformational effects in hydrocarbons in which evaporation rate is treated with a temperature-dependent term has been proposed by Elliott and War 18,23 . We used equation (3) Table 2 and results in Fig. 4 for two coupled conformers of each component in each phase. Indeed the hybrid QTST/IKGT method is inspired by discrete methods such as the Lattice Boltzmann methodology [45][46][47][48][49] and Lattice Boltzmann simulations [50][51][52] . The current hybrid methodology is explained by jiggling and wiggling of atoms in a few discretised conformers in the very vicinity of fuel droplets surfaces (both in pure 14 and binary fuels). Their energies have been very well quantized and there is therefore no continuity (see Table 2). While the conversions of these conformers in the gas and liquid phases are taking place easily, their conformational changes at the interface need to pass specific pathways which are switched on-off. The conformational changes in conformers, collision rate effects and equilibrium vapour concentrations of the components in the gas phase play key roles to make ready these nano-pathways for the dynamic coupling between gas and liquid phases which are really important in the phase transitions (e.g. evaporation). Following from Fig. 4a, the model provides the results (red and green solid lines) which are fitted very well with the experimental evaporation rate (shown with star and cub symbols) respect to the temperatures 670-970 K once the X C16 and X C7 are respectively equal to 0.08 and 0.04 at pressure 0.5 MPa and 0.04 and 0.005 at pressure 0.1 MPa. Deviations are significant when these equilibrium vapour concentrations are a little change. The same scenario took place when the pressure dependency of evaporation rates of C7 and C16 was studied (see Fig. 4b).

Conclusion
In summary, the results of QTST/IKGT which have been inspired by the "discrete" Boltzmann method provide a new important physical concept for understanding dynamic coupling between liquid and gas phases during evaporation of multicomponent fuel droplets. This QTST/IKGT level of thermal and mass transport description in the vicinity of evaporating and condensing fuel droplets indicates two concerted-distinct hydrocarbon topologies in each component for coupling thermal -mass evaporation upon phase transitions. Moreover, we also gain further physical insight into the pathways followed by switching on and off mechanisms at the interface via internal rotations -this insight was previously lacking for multicomponent systems. These pathways are very sensitive to the collision effects, and conformational changes and equilibrium vapour concentrations next to interface. Moreover, the approach presented herein is anticipated to lead to a more refined QTST/IKGT method for reactive multicomponent interfacial transport as simple adsorption -desorption of long chain molecules on a substrate can induce not only conformational changes, but also spontaneous breaking of covalent carbon -carbon bonds 53 .

Methods
Thermal Evaporation Coefficient and Reactive MD simulations. In our approach the droplet (see Fig. 5) was first minimised and subsequently pre-equilibrated to desired temperatures of 350 and 400 K. The Berendsen's thermostat 40 controlled the kinetic energy of the system by scaling the velocities. A Velocity-Verlet algorithm was used to integrate the equations of motion. After equilibrating the systems, the interface layers were strongly coupled with thermostat (with relaxation time τ T = 1 fs) while the rest of the system was weakly coupled with τ T = 100 ps. The "coupling time constant", τ, was used to estimate the time evolution of temperatures based on this equation: where τ = 2τ T C V /(N f k B ), C V is the specific heat capacity at constant volume, k B represents the Boltzmann constant, and N f is the number of degrees of freedom of the system. The time constants, by which systems are allowed to reach the quasi-equilibrium state in micro-canonical conditions (NVE), clarified the β T , for which energy transformations were considered via the interface in a non-steady way and exchanged suddenly. This method allowed us to study gradients of temperature during the evaporation/condensation processes in the vicinity of the liquid-gas interface. We used the Amsterdam Density Functional (ADF) package 54 for all ReaxFF simulations. The temperature in the system under consideration is estimated based on the analysis of various parts of the system (e.g. interface, gas, liquid, drop) separately. The analysis of the interface has been performed only for the molecules which stay in the drop during the whole simulation, ignoring the molecules which leave the droplet. The average energy of gas molecules was obtained based on gas (vapour) temperature which was determined from the conservation of energy: g total total drop drop g where subscript 'drop' refers to the sum of the interfacial layer and liquid phase as shown in Fig. 5. The number of evaporated molecules (N g ) was estimated based on a cut-off distance at which molecules belong to the drop or to the gas phase. It was set to 0.5 nm as inferred from the pair correlation function (g(r)) of n-dodecane (see Fig. 6).

Mass Evaporation Coefficient and Quantum/Statistical Mechanics Methods.
We have first estimated interfacial width, δ, which was unknown in the equation (2), using the following equation: [55][56][57][58] c where σ is a temperature-independent diameter parameter of the methylene and methyl functional groups in n-dodecane conformers, which is assumed to be 3.93 × 10 −10 m; a = 1.16 m and υ = 0.5 are constants and T c = 658.15 K is the critical temperature for n-dodecane. We apply the multistructural statistical thermodynamic Figure 5. A schematic view of a nano-droplet. The droplet has a diameter of 10 nm (96900 atoms); the liquid phase is surrounded by the interfacial layer of thickness of about 1.7 nm when the system is heated up to 400 K. The location of the Gibbs dividing surface that corresponds to the area where the density is equal to 0.5 (ρ liq + ρ vap ) is used to estimate thickness of interfacial layer (see equation (19) for more detail).
Scientific RepoRts | 7: 5001 | DOI:10.1038/s41598-017-05160-z method 59 alongside density functional theory to calculate the Gibbs free energies of n-Dodecane conformers in the gas phase (G g (T)); 60 where − Q g MS T represents the multi-structural partition functions in the gas phase in which rotational, vibration, conformational and torsional effects have been taken into account based on the following formulae: where k B is the Boltzmann constant, and U i is the energy of the i th conformer, N is the number of conformers and φ i,τ is a factor that takes account of torsional potential anharmonicity. Q rot,i is a classical expression for the rotational partition function for conformer i; where F and ω i,l indicate the number of degrees of freedom for vibration modes and vibration frequency of the l th mode of the i th conformer, respectively. To calculate the Gibbs free energies of each conformer at the interface (G i (T)) we employ a modified version of continuum solvation model SMD 35,61 in which some correction terms in temperature dependence of interfacial density and surface tension have been taken into account. SMD is based on the solute electron density, the dielectric constant and the atomic surface tension. The temperature dependence of the surface tension is included using the following formula: 62 c n where B and n are constants: B = 80.1946*10 −3 kcal/(mol*Å 2 ), n = 1.3325, and T c is the critical temperature of n-dodecane. The temperature dependency of interfacial density of n-dodecane is also computed with the self-consistent reaction field (SCRF) method, implemented in the Gaussian 09 suite 63 . The interfacial density, ρ(z), can be expressed as a hyperbolic tangent function: where superscripts l and g denote liquid and gas phases, respectively, and z 0 is the position of the Gibbs dividing surface. The saturated densities of liquid and gas at temperatures T = 298.15 K to 648.15 K are taken from the NIST 64 . Since the translational motions are suppressed at the surface of liquid and all SMD calculations have also been performed based on existence of a conformer in the cavity, the pressure corrections also need to be taken into account using; Figure 6. Pair correlation function for n-dodecane at 350 K and 400 K as function of r (distance of centre of mass of the molecules).
Scientific RepoRts | 7: 5001 | DOI:10.1038/s41598-017-05160-z where p is the pressure, τ = T c /T, δ = ρ/ρ c , ρ and ρ c = 1.33 mol/dm 3 are the density and critical density of n-dodecane, respectively, and A is the Helmholtz free energy: 65 where p is the pressure in the centre of interfacial layer and ρ is the experimental interfacial density of n-dodecane changing from 372.8 kg/m 3 at 298.15 K to 117.5 kg/m 3 at 648.15 K 64 . The constants n 1 , n 2 ,…n 12 are given in Table 3.
The Gibbs free energy of the ensemble of conformers at the interface was determined by the formula; where subscripts int and g refer to the interface and gas phase.  Table 3. Constants for the Helmholtz free energy of n-dodecane at the interfacial layer 65 .