DropletSMOKE++: a comprehensive multiphase CFD framework for the evaporation of multidimensional fuel droplets

This paper aims at presenting the DropletSMOKE++ solver, a comprehensive multidimensional computational framework for the evaporation of fuel droplets, under the influence of a gravity field and an external fluid flow. The Volume Of Fluid (VOF) methodology is adopted to dynamically track the interface, coupled with the solution of energy and species equations. The evaporation rate is directly evaluated based on the vapor concentration gradient at the phase boundary, with no need of semi-empirical evaporation sub-models. The strong surface tension forces often prevent to model small droplets evaporation, because of the presence of parasitic currents. In this work we by-pass the problem, eliminating surface tension and introducing a centripetal force toward the center of the droplet. This expedient represents a major novelty of this work, which allows to numerically hang a droplet on a fiber in normal gravity conditions without modeling surface tension. Parasitic currents are completely suppressed, allowing to accurately model the evaporation process whatever the droplet size. DropletSMOKE++ shows an excellent agreement with the experimental data in a wide range of operating conditions, for various fuels and initial droplet diameters, both in natural and forced convection. The comparison with the same cases modeled in microgravity conditions highlights the impact of an external fluid flow on the evaporation mechanism, especially at high pressures. Non-ideal thermodynamics for phase-equilibrium is included to correctly capture evaporation rates at high pressures, otherwise not well predicted by an ideal gas assumption. Finally, the presence of flow circulation in the liquid phase is discussed, as well as its influence on the internal temperature field. DropletSMOKE++ will be released as an open-source code, open to contributions from the scientific community.


Introduction
The high energy density of liquid fuels is nowadays exploited in many engineering devices such as diesel engines and industrial burners. Spray injection systems are widely used in order to disperse the liquid fuel in an oxidizing environment and the control of this process is currently an active area of research. In particular, the study of a single isolated droplet allows to neglect complex interaction phenomena (coalescence, breakup, etc.) among droplets and to focus on the evaporation and combustion mechanisms. In the last decades the numerical modeling of droplets evaporation and combustion has considerably improved. Starting from the works of Spalding [1] and Godsave [2] concerning the well known "d 2 law", derived under the hypothesis of a perfectly spherical, motionless and constant liquid temperature droplet, Abramzon and Sirignano [3,4] started to account for a convective flow, providing numerous correlations for the heat and mass transfer outside and within the droplet in presence of a relative gas motion. In this context, Dwyer et al. analyzed the droplet dynamics in high T fields [5,6]. In particular, the transient surface heating has been investigated by Law [7], assuming a constant spatial liquid temperature, while Kotake and Okazaki [8] started to analyze the influence of the liquid internal circulation on the vaporization rate. More recently, Sazhin et al. analyzed the heating and evaporation process [9,10], as well as the effect of thermal radiation [11].
With the growing available computational power, the direct numerical solution of the equations governing the evolution of isolated fuel droplets became affordable, allowing the detailed 1D modeling of spherical fuel droplets in microgravity conditions [12,13,14]. Microgravity conditions can be reproduced conducting the experiments in free-falling chambers [15,16,17], both for suspended and free droplets. Recently, some experiments have been conducted aboard the International Space Station (ISS) for the study of evaporation and multi-stage combustion phenomena.
The absence of gravity allows to adopt a simple 1D mathematical description and to focus on physical aspects such as differential species diffusion, non-ideal thermodynamics and combustion kinetics. In particular, the implementation of detailed mechanisms for gas-phase chemistry in 1D models has paved the way for a better understanding of low-T chemistry [18,19,20] , ignition and extinction phenomena [21,22].
However, most of experiments on droplets are conducted under the influence of a gravity field, since experimental devices are much cheaper and the operating conditions are closer to the real ones where droplets are commonly adopted. 1D models are intrinsically not able to predict physical phenomena such as buoyancy, droplet deformation, relative gas motion and internal circulation, which are always present in real conditions. In particular, the external flow field (both induced and forced) around the droplet strongly modifies the evaporation rate at the interface, while the liquid circulation governs the internal heat transfer.
From the modeling point of view, the description of these phenomena requires at least a 2D multiphase model, able to predict the anisotropic deformation of the droplet and to provide both liquid and gas velocity fields. This must be coupled with a reliable evaporation model to describe the droplet shrinking, the diffusion of vapor towards the gas-phase and its further transport, both by the external convection and by the induced Stefan flow. The energy transfer between the two phases is also necessary, as well as proper correlations for the fluid properties.
Among the several methods available for multiphase flows, the Volume Of Fluid (VOF) methodology [23] is widely known for its simplicity, robustness and especially for the excellent mass conservation properties. Its application for evaporation and condensation phenomena has been investigated in recent years [24,25,26,27]. However, most of these works rely on phase-change models based on experimental data, kinetic theories and semi-empirical correlations. A more general evaporation model based on the interface concentration gradient has been implemented by Shlottke [28], neglecting however the detailed characterization of thermodynamic equilibrium at the interface and without a comparison of the results with experimental data. More recently, Gatha et al. [29,30] used a VOF technique to model droplet evaporation and combustion, limiting however the investigation to microgravity conditions. A major numerical issue of VOF method concerns surface tension. Small droplets dynamics involve strong surface tension forces, which are extremely difficult to model due to the presence of the so-called parasitic currents [31]. These currents tend to numerically deform and eventually destroy the droplet, due to inaccuracies on the evaluation of the surface curvature. Some techniques [32,33] have been developed to mitigate parasitic currents for particular cases (e.g. rising bubbles, inviscid static droplets, capillary oscillations), but small droplets evaporation still suffers from this issue. The rapid size reduction of the droplet makes surface tension to be more and more dominant, further amplifying the problem. Even if some works concerning droplets vaporization with VOF methodology are available in literature [34,35,36], the problem of parasitic currents is often not even mentioned. These works are based on commercial CFD codes (mainly Ansys FLUENT R and COMSOL R ), making very difficult to reproduce the results with open-source codes since the adopted numerical algorithms are not provided in detail. Concerning the numerical framework, promising results have been reached by codes such as FS3D [37] and Gerris [38], whereas the OpenFOAM R environment is still lacking in comprehensive solvers for the detailed analysis of droplets evaporation.
This work aims at presenting an open-source general computational tool (called DropletSMOKE++) for the modeling of isolated droplets evaporation under convective conditions in the presence of a gravity field, overcoming the aforementioned limits of the currently available solvers. In particular it includes: • An accurate interface tracking method, based on a geometrical advection, coupled with the Navier-Stokes equations, numerically solved for both liquid and gas phase.
• A reliable evaporation model based on the surface mass flux, coupled with a species equation, with no need of semi-analytical approaches or correlations based on mass-heat transfer dimensionless numbers.
• The detailed description of the vapor-liquid equilibrium thermodynamics, including fugacity coefficients and Poynting correction calculation using a cubic Equation of State.
• A numerical technique which suppresses parasitic currents and by-passes the problem of curvature evaluation, introducing an external centripetal force to mimic surface tension effects.
The DropletSMOKE++ code is fully implemented within the open-source OpenFOAM R framework to manage the computational mesh and the discretization of governing equations. The numerical code will be freely available and open for contributions from the scientific community. The code version used for this work is provided in the supplementary material.
The paper organization includes a thorough description of the mathematical model, as well as the numerical methodology adopted. The centripetal force substituting surface tension is discussed in detail, including its numerical implementation, its effect on the droplet shape and how it impacts evaporation.
A quantitative assessment of the code has been performed, comparing the DropletSMOKE++ results with a reference finite-differences based solver for droplets evaporation in microgravity conditions [13]. Afterwards, the gravity field has been introduced, allowing the analysis of the natural and forced convection regime.
An extensive validation with available experimental data has been carried on for n-heptane, ndecane and n-hexadecane droplets in a very wide range of operating conditions, both for natural and forced convection. Experimental data from various authors include the normalized equivalent diameter decay and some data on the surface temperature behavior, which have been compared with the DropletSMOKE++ numerical results. Equivalent numerical cases (with the same initial conditions) have been simulated also in microgravity conditions in order to highlight the main differences between the models and the intrinsic inability of microgravity solvers to correctly predict the experimental data. Finally, the presence of internal circulation in the liquid phase is briefly analyzed.

Interface advection
The VOF methodology is often referred to as a "one-fluid" approach, where the two phases are treated as a single fluid whose properties vary abruptly at the phase boundary. A scalar marker function α represents the liquid volumetric fraction, varying from value 0 in the gas-phase to value 1 in the liquid phase. The α advection equation in the most general form is: The source termṁ represents the evaporation/condensation contribution to the liquid. Rewriting the equation in function of α: This equation is solved in two steps, treating the advection and the source terms in a segregated approach, similarly to what happens in an operator-splitting technique. The interface tracking is solved first, using the isoAdvector library developed by Roenby and Jasak [39] : isoAdvector performs a geometric advection of the interface, whose quality is superior ( Figure   1) to the MULES (Multidimensional Universal Limiter with Explicit Solution) compressive scheme by Weller [41] usually used in OpenFOAM R multiphase solvers. The source terms are included in a second step: where the first term accounts for the evaporation, while the second term describes the liquid density variation, and eventually the droplet expansion due to external heating.

Pressure and velocity fields
A single Navier-Stokes equation is solved for both phases: It is convenient to separate the hydrostatic pressure ρg · r from the total pressure p, in order to simplify the definition of boundary conditions for pressure: where p rgh = p − ρg · r is the dynamic pressure and r is the position vector.
The continuity equation (needed for the pressure-velocity coupling) accounts for the density variations and for the abrupt change of density at the interface occurring during the evaporation, generating the Stefan flow: In case of constant density and no evaporation Equation 7 reduces to the continuity equation for incompressible flows ∇ · v = 0.

Temperature and species fields
Additionally, energy and species equations are included: The last term in Equation 8 accounts for the interface cooling during evaporation. The diffusion fluxes j d,i are discussed in the next section.

Evaporation model
In this work, the mass evaporation flux is calculated based on the interface concentration gradient of the vapor. Its evaluation is not trivial, since it strongly depends on how the diffusion coefficients are computed. The species diffusion coefficients D i in the mixture are evaluated based on the binary diffusion coefficients and the mole fractions: which means that the diffusion velocities v d,i must be evaluated based on the mole fraction gradient of the species [42]: The diffusion fluxes j d,i are defined as: This expression of the diffusion flux is also used in Equation 9.
A convective flux j c,i is generated by the interfacial density change during the evaporation process, as stated in Equation 7 : The total evaporating flux j i is the sum of diffusive and convective fluxes: ).
The local surface area per unit volume is usually computed from the α function as |∇α|, since its volume integral is the surface area by definition: The evaporation rate for species i is then: It is important to point out that |∇α| is calculated at the interface cells (where |∇α| = 0), while j i must be computed at the gas cells adjacent to the interface cells. Referring to Figure 2 (a), a specific algorithm has been developed to associate every cell at the interface (represented in dotted blue) to the closest cell in the gas-phase (represented in red) along the interface normal (black arrow). Indicating with subscript int a generic interface cell and with adj the closest adjacent gas-phase cell with respect to int, we have: which has to be evaluated for every cell at the interface.
The total evaporation rateṁ is the sum over the liquid (evaporating) species: When solving Equation 4 it is important to ensure a sharp interface. However, it is possible for α to reach negative values in certain cells, especially ifṁ is high enough. In order to avoid this, m is multiplied by a proper function of α, whose aim is to forceṁ to be zero when α approaches zero. Differently from previous works [27], the evaporating flux is evaluated as follows: where the constant K is introduced to ensure the conservation of the integral mass flux Vṁ dV .
It is simply recovered equaling Equations 18 and 19 after a volume integration: (20) We have chosen the function f (α) as follows: which is a redistribution ofṁ over the interface, as shown in Figure 2. This function goes to 0 and to 1 much faster, which helps to maintain a sharp interface after Equation 4 solution.

Thermodynamic equilibrium
The characterization of equilibrium thermodynamics at the phase boundary is crucial for a correct prediction of experimental data, especially for non-ideal mixtures and high pressure cases.
In this work a detailed thermodynamics is included, based on Peng and Robinson equation of state [43] for the evaluation of the compressibility factor Z and the fugacity coefficients φ i .
The vapor-liquid equilibrium relation for a two-phase systems is [44]: where p 0 i (T ) is the vapor pressure of species i, φ i is the gas-phase fugacity coefficient for the pure species andφ i is the gas-phase mixture fugacity coefficient. The exponential term represents the Poynting correction, while x i and y i are the liquid and gas mole fractions of species i. The gaseous mole fraction y i is evaluated in a segregated approach, from equation 22: This saturation mole fraction is assigned to the whole liquid phase and then advected in the gas phase through Equation 9, following the approach described in [45].
The gas density is calculated by: where Z is provided by the Peng-Robinson cubic equation of state.
The molar volume v L in the Poynting correction is simply evaluated as: 2.6. The problem of surface tension: introduction of the centripetal force field The droplet-gas interaction is characterized by high density ratios between the two phases and strong surface tension forces, particularly enhanced in small droplet sizes. Even though VOF methods are widely recognized for the great efficiencies in describing topologically complex interfaces, the accurate representation of surface tension is still a major problem in multiphase flows [33].
The continuous representation of surface tension force f s by Brackbill [31] requires an accurate evaluation of the interface curvature κ: where curvature κ is evaluated as: The α gradient evaluation is numerically challenging, because of the discontinuous nature of the α function. A direct implementation of Equation 26 creates a numerical imbalance between the surface tension force f s and the pressure gradient ∇p in the Navier-Stokes equation. As a consequence, unphysical flows and pressure spikes are generated around the interface, whose magnitude is proportional to σ and inversely proportional to the droplet size. For the typical droplet diameters of the experiments (usually no more than 1-1.5 mm) these flows can strongly deform the interface and eventually break the droplet apart. Several methods have been proposed in this context, such as balanced force algorithms [46] to reduce the numerical imbalance in surface tension force evaluation, height functions [46] to compute the interface curvature directly from the α field and artificial viscosity methods [47]. A sophisticated combination of quad/octree spatial discretization, height functions and balanced algorithms has been implemented by Popinet [38] in Gerris solver.
Concerning the OpenFOAM R environment, Albadawi [48] proposed a coupled VOF-Level Set to mitigate the parasitic currents, while Raeini [49] introduced smoothing factors for κ field combined with the introduction of a capillary pressure equation and a curtailing of the α marker. Even if these latter methods significantly improve the modeling of surface tension forces, they have been validated on bubbles growth, stationary droplets and capillary flows, whereas no test case can be found on small droplets dynamics. Their application on the cases of our interest did not provide satisfactory results, mainly because the droplet shrinking due to the evaporation exponentially amplified the problem: the lower was the droplet size, the stronger were the parasitic currents around the interface, even if the initial stationary droplet was stabilized enough.
In this work, in order to avoid these problems, parasitic currents have been suppressed directly from their source imposing zero surface tension. In presence of a gravity field surface tension is essential to hang the droplet on a fiber or a wire. Therefore, its elimination required the introduction of an additional force to compensate the droplet weight.
A small spherical fiber is initially introduced inside the liquid droplet. A centripetal force directed toward the center of the fiber have been imposed, defined as: where ξ is a scalar potential field defined as: where D f is the fiber radius (D f = 0.05 mm in this work), r is the distance from the fiber center and ξ 0 the intensity of the potential field. The α term in Equation 28 forces f m to be applied to the sole liquid phase and makes f m to be proportional to the droplet volume, as it is for gravity forces.
The result is a force field directed towards the droplet center and only applied to the liquid phase, similarly to what happens in a magnetic attraction of a ferrofluid. The Navier-Stokes equation becomes: It is important to notice that the combination of a geometrical advection (i.e. isoAdvector) and the absence of surface tension σ allows to track the interface maintaining a very high resolution, without any additional correction (such as filtering kernels, smoothing factors etc.) usually needed in the VOF methodology because of the strong gradients involved. The detailed numerical description and the implementation of the centripetal force f m are presented in the next section, as well as a sensitivity analysis on the value ξ 0 to be used in the simulations. The effect of the initial droplet shape on the evaporation will be also discussed.

Fluid properties
The evaluation of transport and thermodynamics properties is based on the OpenSMOKE++ library [50]. Gas properties (diffusion coefficients, thermal conductivity, heat capacities and viscosity) are based on the kinetic theory of gases, while liquid properties (vapor pressure, density, conductivity, heat capacity, viscosity and vaporization heat) are evaluated based on the correlations available in the Yaws database [51].
Mixture properties to be used in the governing equations can be then computed. For a generic property χ: with χ = ρ, µ, C p , k, D i .

Description of the DropletSMOKE++ solver
The DropletSMOKE++ code is embedded within the OpenFOAM R framework which allows to manage the spatial discretization of the governing equations on arbitrary geometries. The PIM-PLE algorithm [40], a combination between SIMPLE (Semi-Implicit Method for Pressure-Linked Equations) and PISO (Pressure Implicit Splitting of Operators), is used to manage the pressurevelocity coupling, computing a velocity field which satisfies both momentum and continuity equation through an iterative procedure.   Finally, Gauss linear upwind scheme is used for spatial discretization of convective terms, while an orthogonal correction is adopted for Laplacian terms.

Computational mesh
The computational mesh used for all the numerical cases proposed in this work has been built with the commercial CFD code Ansys FLUENT R v17.2. The geometry is 2D and axisymmetric, representing a slice of a cylinder having a radius W = 5 mm and a height H = 30 mm. At 10 mm from the base a spherical fiber (diameter D = 0.1 mm) is introduced, which is needed as a source for the force field f m previously described. The liquid droplet will be place around this small spherical fiber. The geometry is meshed in a non-structured way (Figure 4), with a particular attention to the gas-liquid boundary: the mesh is structured and orthogonal on the cylinder boundaries. The region including the liquid droplet and the gas in its proximity is meshed with a concentric pattern, getting finer while approaching the inner fiber. Therefore, the presence of a spherical boundary is also necessary for the mesh construction.
In this work, the mesh is composed by 70,000 cells, with a maximum mesh non-orthogonality equal to 56.6 and a maximum Skewness of 1.3. It has been verified that the numerical solution does not change doubling the number of cells, proving a complete mesh independence.    [13]. The related plots are reported in Figure 6.

Boundary conditions
The definition of proper boundary conditions is fundamental to reach reliable results. The geometry has four boundaries, shown in Figure

Initial conditions
Referring to the mesh presented in Figure 4, the liquid droplet must be positioned around the spherical fiber. In Figure 5 the initial conditions for the evaporation of a cold n-heptane droplet in a hot environment are presented as an example.

Code parallelization
The numerical cases presented in this work are run on a multi-processor (Intel Xeon X5675, 3.07 GHz) machine. The Domain Decomposition Method is used to split the mesh into sub-domains and allocate them to separate processors. The DropletSMOKE++ code can then run in parallel mode, with communication between processors with MPI communication protocol, allowing a significant reduction of the computational time. The optimal number of processors is found to be around 12, which corresponds to an average of 6000 cells for each processor and a speed-up performance ∼ 6 times higher with respect to the serial mode.
The simulations required from 2 to 20 hours, depending on the operating conditions.

Validation of microgravity cases
In this section DropletSMOKE++ is validated in microgravity conditions (imposing g = 0), against the numerical results of the solver developed by Cuoci et al. [13]. This latter code has been validated in a wide range of operating conditions over the last 10 years [22,21] and represents a reliable reference for comparison. This validation is done in order to assess the validity of the equations (in terms of diameter decay, temperature profiles and vaporization velocity) before the activation of the gravity field. Obviously, the centripetal force f m = 0 in these cases.
The 4 numerical cases analyzed are listed in Table 2. Three liquids have been chosen (ndecane, n-heptane and water) at different temperature and pressure conditions. Numerical results are summarized in Figure 6.
For each case (lines) three plots are reported (columns): the temporal evolution of the droplet size, the liquid surface temperature and the vaporization velocity. More precisely, this latter is the convective velocity of the gas generated at the phase boundary.

Results and discussion
For all cases the squared diameter, after the initial transient period, reduces linearly with time, as the d 2 law states. A wet-bulb temperature is reached by the interface and maintained along the whole evaporation process, due to the equilibrium between the incoming heat flux and the outcoming vaporization enthalpy. The transient period depends on the initial conditions. In cases 1-2-4 (initially isothermal) the interface is immediately cooled, diminishing the vaporization flux.
In case 3 the liquid is firstly heated (reducing ρ L and increasing the droplet size) providing an increasing vaporization flux.   The agreement between the models can be considered satisfactory, even if small differences exist.
In particular the surface temperature tends to be slightly over-predicted, especially in the final steps of the evaporation. This is probably due to the spherical boundary of the fiber in the DropletSMOKE++ simulations (Figure 4), which is absent in the 1D model. The presence of a small sphere inside the liquid provides a lower droplet mass if compared to a 1D droplet with the same initial diameter. As can be seen in Figure 6 its effect is completely negligible for most of the simulation (the "missing" mass is less than 0.01% of the initial total mass). While the droplet surface approaches the fiber, this effect becomes important: the incoming heating flux is now distributed over a lower liquid mass, providing a slightly higher equilibrium temperature.
This will not be a problem when modeling evaporation under convection, since fibers are always present in experiments to hang liquid droplets.

Implementation of the centripetal force f m
In order to sustain the droplet against the gravity field, a centripetal force f m (Equation 28) is introduced, based on a potential ξ defined by Equation 29, represented in Figure 7 (a) and previously described in the mathematical model. The liquid droplet is then fixed around the fiber (Figure 7 b). Inside the droplet a small pressure gradient is established (from ∼100 to ∼500 Pa, depending on the value ξ 0 ), similarly to what happens for a liquid column subjected to a gravity field (Figure 7 c). This is actually consistent with surface tension forces, which provide an internal pressure higher than the external one (capillary pressure).

Control of the droplet shape
A droplet suspended under the effect of a gravity field reaches a steady-state shape imposed by the equilibrium between the liquid weight, the pressure gradient and the surface tension forces: In this work surface tension in not considered, and the centripetal force f m governs the equilibrium shape of the droplet: The intensity of f m can be controlled by the value of the constant ξ 0 , which defines the potential field ξ (Equation 29). Figure 8 shows the equilibrium shape assumed by a n-heptane droplet (1 mm diameter, ambient T ) for different values of ξ 0 , from 0.35 to 2 m 2 s 2 . For values ξ 0 < 0.35 m 2 s 2 the droplet falls down, while for ξ 0 > 2 m 2 s 2 the droplet shape remains spherical. These threshold values have been found to be nearly the same for different liquids.
The effect is very similar to surface tension forces: for low values of ξ 0 (or surface tension σ) the droplet tends to assume an elongated shape, because of the dominant gravity forces. Increasing ξ 0 (or surface tension σ), the droplet approaches a spherical geometry. In this work ξ 0 represents a degree of freedom of the problem, which we can saturate imposing a ξ 0 value that better describes the droplet shape of the experimental cases.

Effect of the droplet shape on evaporation
In order to set a proper initial shape of the droplet, a sensitivity analysis on the evaporation process has been carried out. Two numerical simulations of an evaporating n-heptane droplet have been run for the threshold values ξ 0 = 0.35 m 2 s 2 and ξ 0 = 2 m 2 s 2 , (Figures 8 a, e). The initial droplet diameter D 0 = 1 mm, the liquid temperature is T L = 300 K, the ambient temperature is T G = 364 K and the pressure is p = 20 bar.
From Figure 9 we can see that the droplet shape does not have a major impact on the evaporation process. The surface area increases for lower values of ξ 0 , explaining the slightly larger evaporation rates (Figure 9 c), while the surface temperature (Figure 9 b) is almost insensitive to ξ 0 .
We expect the experimental cases of our interest to behave somewhere in-between these two cases, closer to the upper limit described by the black lines in Figure 9 In order to verify this hypothesis, a simple analysis has been conducted on the droplet shapes usually involved in experiments. 10 pictures of suspended droplets available in literature from experimental works [52,53,54,55,56,57,58] have been collected (Figure 10 b) and organized based on two dimensionless numbers, the sphericity ψ and the Eötvös number Eo, defined as: where D y and D x are the maximum lengths of the droplet, in vertical and horizontal directions.
The sphericity ψ is a measure of how the droplet shape approaches that of a perfect sphere (which has ψ = 1). The Eötvös number Eo is the ratio between gravitational and surface tension forces.
A further classification is introduced based on the tipe of supporting fiber (vertical, horizontal or cross fiber). The 10 experimental cases are reported in Figure 10  The Eötvös number is below 1 for all cases, indicating that surface tension forces always dominate over gravity. The sphericity 0.86 < ψ < 1.12 indicates nearly spherical shapes for all the droplets, especially when cross fibers are used (black points). As expected, the horizontal fibers tend to slightly deform the droplet along D x direction, providing ψ < 1 (white points) while the vertical ones along D y , providing ψ > 1 (orange points). Figure 10 clearly shows that most of experimental cases lay between the two lines ψ = 1 (Figure 8 e) and ψ = 1.2 (Figure 8 a). A rig-orous modeling of the surface topology would require different computational meshes to describe the large variety of fibers used to support the droplets (vertical, horizontal, thermocouples, cross fibers...), as well as precise data regarding the dynamic contact angle between the solid and the different liquids. In order to simplify the problem and based on Figures 9 and 10, we assume the evaporation regime not to be significantly influenced by small deviations of the droplet shape from the spherical one, which probably better describes on average the different experimental cases analyzed in this work. For these reasons we decided to adopt an intermediate value ξ 0 = 1 m 2 s 2 (Figure 8 d).

Description of the experimental cases
The numerical model is now validated against 11 sets of experimental data in natural convection regime, taken from 4 papers in literature from 1997 to 2018. All the experimental cases are summarized in Table 3.
Ghassemi [52] used a fine quartz fiber of 0.125 mm of diameter (similarly to the spherical fiber used in our simulations) to hang and evaporate small n-heptane droplets in a hot nitrogen environment.
Nomura et al. [59] examined the evaporation of cold n-heptane droplets in a hot air environment from 300 K to 450 K at various pressures (from 5 atm to 50 atm). The experimental data are taken from the paper of Gogos et al. [60], who proposed a very complex axysimmetric model to model them. An aluminum cross-fiber support frame has been utilized by Verwey et al. [61], in order to maintain a spherical shape of the droplet. N-decane and n-heptane evaporation has been investigated. Finally, n-hexadecane droplets evaporation was examined through a horizontal fiber by Han et al. [55]. Where needed, electric furnaces have been used to heat up the environment and nitrogen to pressurize the vessel. High-speed video cameras have been adopted to follow the droplet lifetime. For all cases, the experimental data provided concern the droplet "equivalent" diameter decay. Han et al. [55] also provided some data on the mean droplet temperature by using a thermocouple as a fiber.   To our knowledge, the cases from Verwey et al. [61] and Han et al. [55] have never been modeled.
Experimental data from Ghassemi et al. [52] and Nomura et al. [59] have only been approached so far either with simplified semi-analytical methods, or specific correlations based on mass-heat transfer dimensionless numbers for the evaporation sub-model [62,63].
In principle, the DropletSMOKE++ code can handle the thermal perturbation of the fiber on the droplet. This could be done either imposing a heat flux of the sphere boundary ( Figure 5) which mimics its presence, or directly meshing the fiber geometry and solve the multi-region heat transfer between the fluid and the solid. Because of the small size of the fibers used in the experiments and the relatively low temperatures involved (since we are not modeling combustion processes), we found the evaporation process to be marginally affected by the thermal perturbation of the fiber, deciding to neglect it and postponing the problem to future works.

Numerical simulations: results and discussion
For each one of the cases presented in Table 3, an equivalent simulation (imposing the same initial conditions) has been performed with the microgravity solver of Cuoci et al. [13] in order to highlight the impact of gravity on the numerical simulation and better support the need of a multidimensional CFD model.
In Figure 11 the numerical results of the simulation of case 1 of Ghassemi et al. [52] are reported as a reference example for a qualitative analysis. The liquid volume fraction (α), the n-heptane  Table 3.
mass fraction field, the temperature and the vectorial velocity field are reported respectively in The droplet is initially at an ambient temperature T L = 300 K and it is progressively heated by the external environment at T G = 773 K (Figure 11 c). A small amount of vapor is released from the droplet surface (Figure 11 b). The gas around the droplet is cooler than T G and the n-heptane vapor has a higher density. This density gradient induces a downward flow field (Figure 11 d) around the droplet, generating a laminar natural convection regime. After a transient period, an equilibrium temperature is reached since the equilibrium mass fraction vapor in Figure 11 (b) reaches a steady state value. Figure 12 reports the numerical results of the cases from Ghassemi et al. [52] (cases 1, 2), compared against the available experimental data. DropletSMOKE++ can predict the experimental data with a reasonable accuracy for both cases. The linear behavior of the squared diameter can be easily recognized, as well as its initial increase due to the droplet heating shown in Figure 11.
The numerical simulations of the cases from Nomura et al. [59] (cases 4, 5, 6) are reported in (c) Figure 13: Model prediction of the experimental data of Nomura et al. [59], cases 3, 4, 5 of Table 3. The reason is straightforward: the Grashof number increases with the gas density (with a power 2), which is proportional to pressure: Therefore, the presence of natural convection has a stronger impact in high pressures cases, if compared with the same cases in microgravity (where Gr = 0).
In Case 5 (Figure 13 c) we slightly overpredict the droplet lifetime, especially if compared to the other cases from the same author. It is worth noticing that this experimental case has been also modeled by Gogos et al. [60], using a complex non-VOF axisymmetric model and assuming a perfectly spherical droplet. His model generally shows good agreement with experiments, but the results of this specific case show deviation from the experimental data very similar to ours.
The cases from Verwey et al. [61] are carried out in mild conditions, characterized by a low vaporization rate due to the low temperatures and moderately high pressures involved. This is clear from Figure 14, where the evaporation times are much longer with respect to the previous cases. Cases 8 and 9 are particularly interesting, since the microgravity solver predicts a quasi steady-state condition. An evaporation flux is still present, but it is extremely low due to the purely diffusive regime involved (the Stefan flow in this conditions is completely negligible). It  n-Hexadecane droplets evaporation from Han et al. [55] have been modeled ( Figure 15). Experimental data on the mean droplet temperature are also available. There is a good agreement with the experimental results. The droplet temperature is slightly underestimated, especially towards the end of the simulation. Experimentally, this is due to the exposure of the thermocouple to the hot gas when the droplet is almost entirely consumed, or to the additional heat flux on the as will be further discussed later. This provides an average liquid temperature which is higher if compared to the 1D model. Moreover, the temperature distribution predicted by DropletSMOKE++ is more homogeneous, while stronger T gradients characterize the 1D model.
Whereas the evaporation rate and the droplet lifetime are strongly influenced by natural convection, in terms of droplet surface temperature there is no significant deviation between the 1D model and DropletSMOKE++. More precisely, the initial droplet heating is faster (since evaporation is here negligible), but the steady state wet-bulb temperature is the same. The experimental data confirm this fact (Figure 15 b, e).
In steady state conditions, there is an equilibrium between the heat flux on the droplet surface  Figure 16 (b) where ρ i,s = from which it is easy to derive h M in function of h T : substituting h M in Equation 37: The resolution of this equation provides the equilibrium temperature T . Equation 40 clearly indicates that the wet-bulb temperature does not depend on the external fluid flow, but only on the properties of the fluids and on the bulk temperature T G .

Description of the experimental cases
To conclude this work, two cases in forced convection have been simulated and compared with experimental results taken from Yang et al. [64]. According to our knowledge, this is the first time that these cases are numerically modeled. n-heptane and n-hexadecane droplets have been evaporated in an upward hot gas flow generated by a heating coil. The experimental data are given for different supporting fibers diameters. The data chosen for comparison refer to a fiber diameter of 0.15 mm, which is the closest to the fiber dimension used in this work (D = 0.1 mm).
The experimental cases are summarized in Table 4.

Numerical results and discussion
The boundary conditions presented in Table 1 impose a fixed upward velocity at the inlet boundary, as well as a fixed temperature value. In Figure 16 (c) we can see the upward velocity field and the relative effect on the vapor cloud around the droplet and on the temperature field.
The agreement with the experimental data (Figure 16 a, b) is satisfactory, with a slight overestimation of the initial heat up in case 2.
From Figure 16 (c), it is clear that the liquid phase can be considered almost uniform in these conditions. This is further analyzed in Figure 17: the flow field streamlines are reported for the liquid phase, colored by the velocity magnitude (Figure 17 a) and temperature (Figure 17 b) for case 2 (Table 4) at time t = 3 s. The external convective flow induces a shear stress on the surface providing a strong circulation inside the droplet, with a maximum magnitude of ∼ 2 cm/s. The resulting toroidal structure is well known in literature, and it is usually approximated as an inviscid Hill's spherical vortex [4,65] to derive sub-models of droplet heating. It has been demonstrated that internal circulation needs to be included in modeling droplet vaporization [8], and this is usually accomplished with effective conductivity models [3]. Since the entire flow field is directly solved in this work, the internal temperature field can be predicted with no need of such additional sub-models. The maximum temperature difference inside the liquid phase is only ∼ 3 K (Figure 17 b). As expected, a more intense heat transfer occurs at the bottom of the droplet, where the flow field impacts the liquid.
Finally, it is important to point out that the suppression of surface tension forces adopted in this work could in principle affect the internal motion in the liquid phase, because of the different the tangential shear stress at the interface. However, our comparative analyses in this sense, not reported in this work, show that the results obtained with DropletSMOKE++ are in a reasonable agreement both with experimental [66,67] and numerical [36,68] works not based on the VOF methodology. These investigations are still preliminary, but indicate that even neglecting surface tension the internal circulation phenomenon can be quantitatively and qualitatively captured with a reasonable accuracy.

Conclusions
In this work we presented a multidimensional CFD numerical framework for the evaporation of isolated fuel droplets. The solver (called DropletSMOKE++) accounts for the presence of external flow field and for the influence of a gravity field, allowing to numerically model the evaporation of droplets both in natural and forced convection. The code is embedded within the well-known OpenFOAM R framework, open for contributions from the scientific community.
The core of the code is based on the VOF (Volume Of Fluid) methodology, widely known for its simplicity and mass conservation properties. The interface is geometrically advected and subjected to source terms to account both for evaporation and liquid expansion. This main solver is coupled with an evaporation model, which directly computes the vaporization rate from the interface molar fraction gradient, and with energy and species equations.
The main novelty of this work is the surface tension treatment. The parasitic currents, still an open problem for the modeling of small droplets dynamics, have been completely suppressed by-passing the modeling of surface tension effects. A centripetal force directed towards the droplet center has been introduced to stabilize the droplet position against the gravity force. This expedient allowed to accurately model the evaporation process, whatever the droplet size. By tuning this centripetal force, it is possible to modify the initial droplet shape to better match the experimental conditions (where the droplets are nearly spherical).
The DropletSMOKE++ simulations has been compared with a well-known 1D microgravity solver, in order to verify the correct implementation of the governing equations, obtaining a very good agreement between the models. Afterwards an extensive validation against experimental data has been carried out, both for natural and forced convection cases available in literature, over a very wide range of temperature, pressure and initial diameter conditions.
The agreement with experimental data was excellent and allowed to analyze the impact of an external flow field on the evaporation process, in particularly compared to the same cases modeled with a 1D model. The droplet lifetimes are largely overpredicted by the 1D model, while the CFD simulations carried out with DropletSMOKE++ are able to correctly reproduce the data. This difference between the models is particularly enhanced at high pressures and in slow evaporating conditions (low temperature and high pressure).
Finally, the effect of the liquid internal circulation has been highlighted. In particular a much more uniform temperature field is provided by the internal heat transfer, strongly enhanced by the vortical structure of the internal flow field.
Future works will concern the introduction of a gas-phase kinetic mechanism to simulate the combustion process in presence of a convective flow, as well as the investigation of the liquid phase chemistry.

VOF Volume Of Fluid
Greek letters