Numerical simulation of droplet evaporation using interface tracking methods: phase change modelling

The present work gives a description and comparison of recent methods for the direct numerical simulation of droplet vaporization using a two-scalar approach for energy and species equations. Those methods require a choice of modelling for the phase change process. The comparison relies on the work of Palmore at al. [1] in the context of Volume of fluid and on the phase change procedure of Rueda Villegas et al. [2] in the Level set framework. Both approaches have been implemented and coupled to the same two-phase low Mach solver. A quantitative analysis is given on the canonical Stefan flow problem where analytical solutions are available. First the 1D planar Stefan problem is performed to avoid any errors associated to topology. Then, the 3D spherical Stefan problem is tested to evaluate both methodologies on a multidimensional case where the choice of interface representation is prevalent. Finally, a 2D droplet convected in a quiescent ambient gas is presented. This last test case aims to demonstrate robustness of the solver on a more demanding test case implying convection, interface deformation and non homogeneous vaporization.


Introduction
Direct numerical simulation (DNS) of droplet evaporation has been a growing subject of interest in the last decades with the emergence of multiple DNS solvers and various associated numerical methods. Two-phase flow simulation is a very challenging subject because of the singular pressure jump and flow properties at the interface. This is even more challenging when dealing with phase change, with the apparition of new interface jumps, in velocity and temperature gradients. The sharp treatment of the interface seems to be a promising method to deal with such discontinuities. Numerous recent numerical approaches have shown very encouraging results using either Volume of fluid (VOF) [1], Level set (LS) [2] or Front tracking [3]. However, some numerical aspects of droplet evaporation simulations remain open subjects. For instance, the choice of an interface tracking method conditions the numerical methods associated to phase change. To handle discontinuity and compute quantities at the interface, the literature often relies on the PLIC reconstruction [4] in a VOF framework while in a LS framework, the distance function is used as in the Ghost Fluid Method (GFM) [11]. This work aims to compare these 2 paradigms in the context of general vaporization studies. In Sec. 2, the physical problem of phase change in a low Mach context is presented. A brief solver description is given in Sec. 3 while numerical methods associated to phase change modelling are detailed in Sec. 4 and 5. Finally a comparison of VOF and LS solvers is provided in Sec. 6 on planar and spherical Stefan flow problems and on a 2D convected droplet.

Set of conservation equations
The two fluids are described by a set of conservation equations, the condensable fluid is assumed to be monocomponent, and the gas part of the fluid does not react with the bulk gas. In this work, the following system of conservation equations is considered with ρ the density of the fluid, u the velocity of the fluid, S = µ ∇u + ∇u − 2 3 (∇ · u) I , p the pressure, g the acceleration of gravity, T the temperature, c p the specific heat, k the thermal conduction, Y the species fraction of condensable fluid in the bulk gas and D m the species diffusivity.

Interface representation
The interface evolution is deduced from an interface tracking method leading to the evolution of c = f the volume fraction or c = φ the signed distance function with u Γ the interface velocity. To solve equation (5), u Γ is deduced from mass conservation across the interface with u l and u g the liquid and gas velocities respectively and n Γ the interface normal pointing outward the liquid phase.
This leads to u Γ = u l −ṁ ρ l n Γ and u Γ = u g −ṁ ρ g n Γ .
The set of conservation equation have associated jumps at the interface. For a given quantity Φ, it is defined as withṁ the vaporization rate and n Γ the interface normal, σ the surface tension (here considered constant), κ the curvature, L vap the specific latent heat boiling, T Γ and T sat the interface and saturation temperature respectively.

Definition of the vaporization rate
The evaporation rate has two definitions based on the heat and mass flux jumps at the interface (9) and (10). This corresponds to different regimes of vaporizatioṅ Then, the computation of the evaporation rateṁ is either based on equation (11) or (12). The first choice can lead to non-physical results if one imposes T Γ = T sat , more specifically when the interface temperature is very far from the boiling point. The second choice suffers from lack of numerical boundness in the saturation limit (the denominator 1 − Y Γ tends to zero), in fact, at saturation the vaporization is only driven by the heat flux at the interfaceṁ = [k∇T · n Γ ] Γ L vap and equation (12) is not defined. A methodology is proposed in [5] to address this issue by switching between both vaporization regimes when a critical temperature close to T sat is reached. Another way to deal with this issue is to always use equation (11) forṁ computation and to provide a closure for T Γ which can be very different from T sat as in [1] and [2].
Closure for Y Γ and T Γ Considering that the pressure is at saturation at the interface, T Γ and Y Γ are related through the Clausius-Clayperon relation with W vap and W g the molar masses of vapour and ambient gas respectively and R the gas constant. r P defines the ratio between saturation and ambient pressures.

Global solver description
A brief survey of the solver discretization used in this work is given without considering the phase-change modelling. Our solver is similar to the one presented by Palmore et al. in [1] for momentum and scalar transport. Here, it is either coupled with VOF or LS.
The VOF scheme used in this work is the dimensional-splitting method of [6] which allows to conserve mass up to machine precision while keeping a fairly simple and efficient implementation. For LS, the standard method of [7] is used here with 2 reinitialization iterations after each time step. This is enough to conserve the distance properties of φ. The interface velocity is It has been previously noticed [8] that using a liquid velocity with divergence-free extensions improves significantly mass predictions in the context of evaporation. Here, u l is obtained through the procedure proposed in [1] where liquid velocity is extrapolated in the gas and projected to its divergence-free form.
The momentum equation (2) is solved with a one-field representation of the velocity using a classical projection method [9] to ensure mass conservation. The prediction step uses a massmomentum consistent methodology for the convective term [10]. The convective fluxes are computed using WENO5 and degenerate to first order upwind when the stencil implies points from both liquid and gas. The diffusion term follows second order discretization. The pressure jump condition is prescribed using GFM [11] with a curvature computed from the algorithm proposed in [12] in VOF or from finite differences in LS. The velocity jump is imposed as in [1] for VOF and the procedure described in [2] is used for LS.
For jump condition treatment in the energy and mass species equations, one approach is to discretize jump conditions as source terms [13] in the conservation equations (3) and (4). However, this leads to a mixing of phase quantities on two or three cells around the interface. In [14], the issue is handled by splitting the temperature into a liquid part and a gas part. Both quantities are transported independently and coupled through the Dirichlet boundary condition T Γ at the interface. In the same manner, the species equation can be considered only in the gas phase with a boundary condition at the interface. Here, two paradigms are presented for handling the two-scalar energy and species equations. The first one makes use of the PLIC reconstruction to solve scalar at their phase barycentre in a VOF framework [1] while the other one divides the domain between two phases based on the sign of the distance function in a LS framework [2]. In both cases, this leads to the resolution of 3 scalar equations T l , T g and  Fig. 1.
The scalar transport discretization is based on the mass consistent transport of [1]. The convective term is discretized using BQUICK and degenerate to first order upwind close to the interface. The diffusion treatment is detailed in Sec. 4 and 5.

Full phase change procedure for the VOF solver
The procedure used in this work for the VOF solver is the one presented in [1] : • T l and T g are extrapolated in the other phase linearly using the PDE approach of [15] to obtain T ghost l and T ghost g •ṁ is computed with equation (11) with gradients computed from linearly extrapolated and values in a stencil around the considered cell • T Γ and Y Γ are computed using the iterative process of [1] imposing both equations (11) and (12) with the constraint (13) discretely in a given cell • The diffusion problems of T l and T g are solved with the embedded Dirichlet boundary condition T Γ while Y g is solved with the embedded Dirichlet boundary condition Y Γ . The method uses the PLIC representation of Fig. 1 to build a symmetric linear system.
This combination leads to a consistent methodology as Dirichlet boundary conditions T Γ and Y Γ both depend on T l , T g and Y g . Note that the discretization can lead to infinitesimal cells which require to treat the diffusion implicitly to avoid any stability problem.

Full phase change procedure for the LS solver
The procedure used in this work for the LS solver is the GFTSBE of [2] : • T l and T g are extrapolated in the other phase linearly using the PDE approach of [15] to obtain T ghost l and T ghost g •ṁ is computed with equation (11) with gradients computed from linearly extrapolated and values in a stencil around the considered cell • Y Γ is computed directly from Y g with the first order approximation Y Γ ≈ Y g and T Γ is obtained by inverting equation (13) • The diffusion problems of T l and T g are solved with the immersed Dirichlet boundary condition T Γ following [16] while Y g is solved with the embedded Robin boundary condition with the approach of [17] This combination is consistent too as the Dirichlet boundary condition T Γ only depends on Y g while the Robin boundary condition only depends onṁ which is a function of T l and T g . Note that the discretization can lead to infinitesimal distances which require to treat the diffusion implicitly to avoid any stability problem.

Test cases
Validations are performed using the VOF and LS strategies described above. The Stefan flow part aims to give a quantitative assessment of the phase change procedure as analytical solutions are available and allows to study the mesh convergence with respect to the exact solution of the problem. On the other hand, the convected droplet test case demonstrates robustness and feasibility of the solver on more complex configurations.

Stefan flow
Considering a 1D planar Stefan flow allows to put aside all the interface tracking and topological issues of curved interfaces. As shown in Fig. 2, air/water vapor gas and liquid water are separated by an interface with gas phase at the left and liquid phase at the right in a domain L = 1 mm. The hot gases will cause vaporization of the liquid water and the interface will evolve to the right through phase change with a free outlet condition at the right of the domain. As no convection effects are present in the test case, the accuracy of scalar diffusion and interface quantities computation are isolated. The parameters are T ∞ = 323.15 K and Y ∞ = 0.2, this gives T Γ = 296.12 K and Y Γ = 0.22106. The displacement of the front can be derived analytically and is given as x Γ (t) = 2γ k g C p,g /ρ g t with γ = 0.1157 a constant deduced from thermodynamic properties (all details of these derivations are in [1]). Finally, the simulation is performed between t = 0.01 and 0.1s, this gives an initial interface position x Γ = 0.0732 mm. The spherical problem introduces normal gradients which are not aligned with the mesh, it is then important to observe the associated errors compared to a 1D problem. This test case is performed in 3D with an expected diameter evolution following the d2 law. The derivations in [2] with T ∞ = 700 K and Y ∞ = 0.2 leads to interface quantities T Γ = 294.92 K and Y Γ = In Fig. 3 (left), a first order convergence of x Γ for both approaches is presented with a better accuracy for VOF. This first order trend is expected as the vaporization rate is computed from differentiation of T g and T g which are at best second order because of the embedded boundary approach. LS is less accurate because of the computation of T Γ and Y Γ which is inherently first order because of the approximation Y Γ ≈ Y g . This is discussed in [1] where a second order extrapolation is proposed and shows huge improvement in accuracy. The improvement has been tested here and corresponds to the LS-extrap line which shows a better accuracy than the VOF procedure.
For the multidimensional test case, LS approach is unable to predict the regression correctly for low resolutions. This is explained by the reinitialisation step which implies a mass loss which is more important than the regression due to phase change. For completeness, a plot without evaporation has been added on the right of Fig. 3 to show mass conservation convergence of the LS method. It is obvious that the test cases with N/D = 4 and 8 cells in the diameter have a regression driven by the reinitialization. However, for better resolutions, the slope is very close to the expected regression. The fact that the N/D = 16 case is closer than the 32 one is explained again by the mass loss through reinitialization which counter-balance the under-predicted vaporization rate. On the other side, the VOF procedure shows a converging trend but it seems that the asymptotic slope is lower than the exact one. For all meshes, the vaporization rate is under predicted with huge errors for the under resolved test cases. This behaviour was also observed in [1] with a convergence rate lower than one.

Convected droplet
This test case has already been presented in [5], a water droplet with D 0 = 5 · 10 −4 m is placed at the top (2.5D 0 , 19D 0 ) of a domain [0, 5D 0 ] × [0, 20D 0 ] with an initial velocity of 1 m/s. The initial solution is discontinuous with a uniform temperature in the liquid T l = 323 K, and a uniform temperature in the gas T g = 873 K which is only composed of air Y g = 0. The left and right boundary conditions are taken periodic while the upper boundary is an inlet with quiescent air and the lower boundary is a free outlet. Hence, this last case shows the capability of the solver to handle transient vaporization regimes, deformations of the droplet due to the velocity and the surface tension and scalar convection effects on the phase change process. The test case is performed on a 64 × 256 mesh which gives N/D = 12.8 during t = 0.009s. It can be seen in Fig. 4 that both methods are able to handle a transient regime from the discontinuous initial solution to a profile which is characteristic of a convected droplet evaporating. The temperature remains very high at the head of the droplet leading to a strong vaporization rate while at the trail, a smoother evolution is observed leading to a vaporization rate closer to the static evaporation problem. Note that the main difference observed here is the velocity of the droplet which is faster for VOF than for LS. This can be explained by the difference in ρ u definition in the momentum equation which is more accurate and consistent for VOF as it is based on the PLIC reconstruction, hence it leads to a better momentum conservation.

Conclusions
This work compares VOF and LS strategies adapted to the same low Mach solver. The Stefan flow test cases presented here show that VOF has a regression following d2 law even for the most under-resolved test cases. However, the error in the regression slope is important for small resolutions. The LS method seems to give very accurate results when the mesh is fine enough to limit the mass lost through reinitialization. In both cases, a minimal resolution is required to give satisfactory results. The dynamic test case enlightens the importance of the interface tracking method choice. They lead to a different conservation of momentum which has an impact on the vaporization process. This problem arises in every evaporation problem as the error in momentum increases with the density ratio which is of 3 order of magnitude in the target applications. Both VOF and LS shows a good qualitative behaviour for this last case while a discrepancy between both method can be observed in the last snapshots of the simulation when momentum conservation errors are accumulated. For a general vaporization problem, the choice of an approach will depend on the application case and the reachable resolution. DNS of evaporation needs to capture the thermal boundary layer which can be very thin in high convection configurations. The mesh associated to such constraints will naturally lead to a good resolution of the droplet where LS can be very interesting as it is accurate and easy to implement in a low Mach solver. However, in evaporating spray simulations, some structures hold on a very limited number of cells and in that case, both VOF and LS lead to huge errors in regression predictions. VOF will always under predict the regression rate while LS will highly over predict the regression because of the mass lost through reinitialization. Note that this last aspect could be improved with more sophisticated mass-preserving reinitialization schemes.