A volume of fluid framework for interface-resolved simulations of vaporizing liquid-gas flows

TL;DR

This paper presents a comprehensive computational framework for simulating vaporizing liquid-gas flows, accurately capturing interface dynamics and phase change processes using advanced numerical techniques.

Contribution

It introduces a novel, divergence-free velocity extrapolation and sharp interface treatments within a volume of fluid framework for vaporization simulations.

Findings

Framework achieves first-order convergence in accuracy.

Successfully simulates unsteady droplet vaporization.

Handles large density ratios and complex interface phenomena.

Abstract

This work demonstrates a computational framework for simulating vaporizing, liquid-gas flows. It is developed for the general vaporization problem which solves the vaporization rate based as from the local thermodynamic equilibrium of the liquid-gas system. This includes the commonly studied vaporization regimes of film boiling and isothermal evaporation. The framework is built upon a Cartesian grid solver for low-Mach, turbulent flows which has been modified to handle multiphase flows with large density ratios. Interface transport is performed using an unsplit volume of fluid solver. A novel, divergence-free extrapolation technique is used to create a velocity field that is suitable for interface transport. Sharp treatments are used for the vapor mass fractions and temperature fields. The pressure Poisson equation is treated using the Ghost Fluid Method. Interface equilibrium at the interface is computed using the Clausius-Clapeyron relation, and is coupled to the flow solver using a monotone, unconditionally stable scheme. It will be shown that correct prediction of the interface properties is fundamental to accurate simulations of the vaporization process. The convergence and accuracy of the proposed numerical framework is verified against solutions in one, two, and three dimensions. The simulations recover first order convergence under temporal and spatial refinement for the general vaporization problem. The work is concluded with a demonstration of unsteady vaporization of a droplet at intermediate Reynolds number.

Tables

Table 1. Table 1 : Fluid properties for Neumann’s problem.

Property	Units	Gas	Liquid
$\rho$	$\text{kg}/{\text{m}}^3$	$1$	$1000$
$\mu$	$\text{kg}/\left(\text{m}\cdot \text{s}\right)$	$1\times {10}^{-5}$	$1\times {10}^{-2}$
$C_p$	$\text{J}/\left(\text{kg}\cdot \text{K}\right)$	$1000$	$1000$
$k$	$\text{W}/\left(\text{m}\cdot \text{K}\right)$	$1\times {10}^{-2}$	$1\times {10}^{-1}$
$D$	${\text{m}}^2/\text{s}$	$1\times {10}^{-5}$	—
$M$	$\text{kg}/\text{mol}$	$0.018$	—
$\sigma$	$\text{N}/\text{m}$	—	$0.01$
$L_V$	$\text{J}/\text{kg}$	—	$1\times {10}^6$
$T_{boil}$	K	—	$373.15$

Table 2. (a) L 1 subscript 𝐿 1 L_{1} error in gas temperature.

$C_{\lambda }$	MB	REA2
$10$	$1.36e-1$	$1.36e-1$
$25$	$1.65e-1$	$1.66e-1$
$50$	$1.89e-1$	$1.95e-1$
$100$	$2.06e-1$	$2.21e-1$
$200$	$2.22e-1$	$7.25e-1$

Table 3. (a) L 1 subscript 𝐿 1 L_{1} error in gas temperature.

$C_{\lambda }$	MB	REA2
$10$	$1.36e-1$	$1.36e-1$
$25$	$1.66e-1$	$1.66e-1$
$50$	$1.91e-1$	$1.94e-1$
$100$	$2.12e-1$	$2.19e-1$
$200$	$2.27e-1$	$2.39e-1$

Table 4. Table 4 : Fluid properties for d 2 superscript 𝑑 2 d^{2} law problem.

Properties	Units	Gas	Acetone
$\rho$	$\text{kg}/{\text{m}}^3$	$1$	$700$
$\mu$	$\text{kg}/\left(\text{m}\cdot \text{s}\right)$	$1\times {10}^{-5}$	$3.26\times {10}^{-4}$
$C_p$	$\text{J}/\left(\text{kg}\cdot \text{K}\right)$	$1000$	$2000$
$k$	$\text{W}/\left(\text{m}\cdot K\right)$	$5.2\times {10}^{-2}$	$1.61\times {10}^{-1}$
$D$	${\text{m}}^2/\text{s}$	—	$5.2\times {10}^{-5}$
$M$	$\text{kg}/\text{mol}$	$0.029$	$0.058$
$\sigma$	$\text{N}/\text{m}$	—	$0.0237$
$L_V$	$\text{J}/\text{kg}$	—	$5.18\times {10}^5$
$T_{boil}$	K	—	$329$

Equations

$$\frac{\partial \left(\rho \mathbf{u}\right)}{\partial t}+\nabla \cdot \left(\rho \mathbf{u}\otimes \mathbf{u}\right)=-\nabla p+\nabla \cdot \left(\mu \mathcal{S}\right),\text{ where}$$

$$\mathcal{S}=\left(\nabla \mathbf{u}+{\nabla \mathbf{u}}^{\top }-\frac{2}{3}\left(\nabla \cdot \mathbf{u}\right)\mathcal{I}\right),$$

$$\frac{\partial \rho }{\partial t}+\nabla \cdot \left(\rho \mathbf{u}\right)=0,$$

$$\frac{\partial \left(\rho Y\right)}{\partial t}+\nabla \cdot \left(\rho Y\mathbf{u}\right)=\nabla \cdot \left(\rho D\nabla Y\right),$$

$$\frac{\partial \left(\rho C_{p}T\right)}{\partial t}+\nabla \cdot \left(\rho C_{p}T\mathbf{u}\right)=\nabla \cdot \left(k\nabla T\right)+\frac{Dp}{Dt}+\frac{\mu }{2}\mathcal{S}:\mathcal{S},$$

$$\frac{Dp}{Dt}=\frac{\partial p}{\partial t}+\mathbf{u}\cdot \nabla p,$$

$$\mathcal{S}:\mathcal{S}=S_{ij}{S}_{ij},$$

$$\frac{\partial \left(\rho C_{p}T\right)}{\partial t}+\nabla \cdot \left(\rho C_{p}T\mathbf{u}\right)=\nabla \cdot \left(k\nabla T\right).$$

$$\displaystyle f(\bm{x},t)=\left\{\begin{array}[]{cc}1&\text{if }\bm{x}\in L,\\ 0&\text{otherwise.}\end{array}\right.$$

$${\mathbf{n}}_{\Gamma }\cdot {\rho }_G\left({\mathbf{u}}_G-{\mathbf{u}}_S\right)={\mathbf{n}}_{\Gamma }\cdot {\rho }_L\left({\mathbf{u}}_L-{\mathbf{u}}_S\right)=\dot{m}.$$

$$\frac{d}{dt}\underset{L}{\int }{\rho }_{L}dV=-\underset{\Gamma }{\int }\dot{m}dS,$$

$$\frac{\partial {\rho }_L\alpha }{\partial t}+\nabla \cdot \left({\rho }_L\alpha \mathbf{u}\right)=-\dot{m}{\delta }_{\Gamma },$$

$${\delta }_{\Gamma }=\underset{\Gamma \cap \Omega }{\int }dS/\underset{\Omega }{\int }dV.$$

$${\left[\varphi \right]}_{\Gamma }={\varphi }_G-{\varphi }_L,$$

$${\left[{\mathbf{n}}_{\Gamma }\cdot \mathbf{u}\right]}_{\Gamma }=\dot{m}{\left[\frac{1}{\rho }\right]}_{\Gamma }.$$

$${\left[p\right]}_{\Gamma }=-\sigma \kappa -{\dot{m}}^2{\left[\frac{1}{\rho }\right]}_{\Gamma }.$$

$$\dot{m}=\frac{{\left[{\mathbf{n}}_{\Gamma }\cdot k\nabla T\right]}_{\Gamma }}{\left[h\right]},$$

$$\dot{m}=\frac{{\mathbf{n}}_{\Gamma }\cdot {\rho }_{G}D\nabla Y}{Y^{\Gamma }-1},$$

$$X$$

$$Y$$

$$\rho ={\rho }_L\alpha +{\rho }_G(1-\alpha ),$$

$$\mu ={\mu }_L\alpha +{\mu }_G(1-\alpha ).$$

$$\frac{\partial \left(\rho \varphi \right)}{\partial t}=-\nabla \cdot \left(\rho \varphi \mathbf{u}\right)+\nabla \cdot \left(\rho \xi \nabla \varphi \right)+S(\mathbf{x}),$$

$$\frac{{\left(\rho \varphi \right)}^{n+1,k+1}-{\left(\rho \varphi \right)}^n}{\Delta t}$$

$${\varphi }^{n+1/2,k}$$

$$\frac{{\left(\rho \varphi \right)}^{n+1,k+1}-{\left(\rho \varphi \right)}^n}{\Delta t}$$

$${\varphi }^{n+1/2,k}$$

$$\frac{\partial \varphi }{\partial t}=\frac{\partial }{\partial x}\left(\xi \frac{\partial \varphi }{\partial x}\right),$$

$$\frac{{\partial }^2\varphi }{\partial x^2}=\frac{{\varphi }_{i+1}-2{\varphi }_i+{\varphi }_{i-1}}{\Delta x^2},$$

$$-C{u}_{i+1}^{n+1,k+1}+(1+2C)u_i^{n+1,k+1}-C{u}_{i-1}^{n+1,k+1}=u_i^n,$$

$$u^{n+1/2,k}$$

$$\frac{{\rho }_u^{n+1,k+1}-{\rho }_u^n}{\Delta t}+\frac{\partial \left({\rho }_u^n{u}^{n+1/2,k}\right)}{\partial x}$$

$$\frac{{\rho }_u^{n+1,k+1}{u}^{n+1,k+1}-{\rho }_u^n{u}^n}{\Delta t}+\frac{\partial \left({\rho }_u^n{u}^n{u}^{n+1/2,k}\right)}{\partial x}$$

Full text

It will be shown that correct prediction of the interface properties is fundamental to accurate simulations of the vaporization process. The convergence and accuracy of the proposed numerical framework is verified against solutions in one, two, and three dimensions. The simulations recover first order convergence under temporal and spatial refinement for the general vaporization problem. The work is concluded with a demonstration of unsteady vaporization of a droplet at intermediate Reynolds number.

1 Introduction

Several common engineering systems involve the flow of liquids and gases. For many of these flows, it is often the case that the primary objective of these systems is to exchange heat or mass between the two phases. Some common systems include heat exchangers, spray coolers, bubble column reactors, and spray combustors. In these systems heat and mass transfer occurs primarily through the process of phase change. As such, understanding phase change is an important step in better understanding the dynamics of these engineering systems. This paper focuses on liquid-to-gas phase change, here called vaporization. For many of these problems, it is appropriate to conceptualize the flow as consisting of a continuous carrier phase with dispersed droplets. Several reviews have been written on the theory of the vaporization of droplets [7, 8, 9]. These studies were focused primarily on either vaporization for idealized flows. The study of the interaction of vaporization and turbulent flows remains an open problem.

Recently, experiments have been performed studying vaporization and combustion of droplets in turbulent flows (e.g., the review by Birouk and Gokalp [10]), however a fundamental understanding of the mechanism of vaporization-turbulence interaction is lacking. In particular, the question of which scales of turbulence are most important to the interaction is yet to be answered. Direct numerical simulation (DNS) can provide an alternative avenue to study this problem.

In this context, DNS are taken to be those simulations that solve the dynamics of the flow from first principles, i.e., directly from the laws for conservation of mass, momentum, and energy. Vaporization is represented as a set of coupled matching conditions at the liquid-gas interface which ensure the conservative transfer of mass, momentum, and energy across the interface. According to Sirignano [11], a multiphase flow simulation can be considered a true DNS when it resolves both the gas film around the droplet and the internal droplet flow. Unfortunately, due to the sharp discontinuities in pressure, velocity, and density at the interface, the stable numerical simulation of this problem is difficult.

Only recently have numerical methods capable of performing resolved simulations of vaporization have been discussed [12, 13, 14, 15, 16, 17, 5, 1, 18]. The aforementioned numerical frameworks have been used to study film boiling, vaporizing bubbles rising in liquid, and the classical $d^{2}$ law. However, the studies are performed under several assumptions regarding the regime of vaporization and the flow conditions.

First, simulations are typically performed within one of two limiting cases of vaporization: an evaporation limit wherein the vaporization rate is assumed to be limited by the local values in the vapor species [14, 16]; and a boiling limit wherein the vaporization rate is assumed to be limited by the local values in temperature [12, 13, 15, 17, 18]. Few authors have studied vaporization without the imposition of these constraints [16, 5, 1]. All of the frameworks mentioned so far used incompressible solvers, and assumed constant thermophysical properties. Furthermore, most of the studies were performed in two dimensions which inherently neglects many important physical processes in turbulent flows.

The literature also contains very little discussion of the convergence and stability of the numerical methods used in the simulations. The current work demonstrates a numerical framework for vaporization in three dimensions, and develops a suite of tests to study the stability and convergence of the flow solver. In developing the framework, it will be shown that traditional semi-implicit treatments of scalar transport can lead to stability issues for large timesteps.

The work is organized as follows. Section 2 gives the physical description of the problem and outlines the mathematical model of the flow. Section 3 details the algorithms used in the work. Section 4 verifies the numerical algorithm against known analytical solutions in one dimension. Section 5 performs a verification study in two dimensions, and it demonstrates the temporal accuracy and stability of the solver. Section 6 verifies the solver against the well know $d^{2}$ solution in three dimensions. Finally, demonstrations of unsteady droplet vaporization (Section 7) are shown.

2 Physical Description

This section details the mathematical equations used to describe the vaporization process. The process is described by the conservation equations for low-Mach flow, i.e., conservation of mass, momentum, energy, and chemical species in each phase. Thermodynamic equilibrium at the interface is resolved using the Clausius-Clapeyron relation. As an initial assumption, a single component fuel is used, and the fuel vapor does not react with the surrounding gas. Both the liquid and the gas are assumed to be of constant density. The exposition of the governing equations follows.

The motion of a Newtonian fluid is governed by the Navier-Stokes equations,

[TABLE]

which is statement of the conservation of momentum. Here, $\rho$ and $\mu$ are the fluid density and dynamic viscosity; $\bm{u}$ is the velocity; and $p$ is the pressure. The identity tensor is represented with $\bm{\mathcal{I}}$. Equation 1 must be coupled with the continuity equation, here expressed as

[TABLE]

which is an expression of the conservation of mass.

In the gas phase, conservation of mass must be supplemented with an equation for conservation of chemical species,

[TABLE]

where the vapor mass fraction and its diffusivity are represented using $Y$ and $D$, respectively. As the problem is nonreactive, remaining inert gas can be expressed as $1-Y$.

The equation for conservation of energy can be expressed as

[TABLE]

where $T$ is the fluid temperature, $C_{p}$ is the specific heat under constant pressure, and Fourier’s law is used for conduction ($k$ is the thermal conductivity). The notation, $D/Dt$ represents the Lagrangian derivative, i.e.,

[TABLE]

and double tensor contraction expands as

[TABLE]

using Einstein summation convention.

However, the full energy equation is not used in this work. Consistent with the low-Mach assumption, we neglect those terms of the energy equation that are quadratic in the Mach number (or equivalently the velocity), yielding

[TABLE]

In deriving Eq. 7 it was assumed that the pressure scales as the dynamic pressure, $\frac{1}{2}\rho|\bm{u}|^{2}$. The resulting equation neglects all pressure and viscous work on the system and, accordingly, acoustic effects.

2.1 Phase Change Strategy

The conservation equations Eqs. 1, LABEL:, 2, 3 and 7 do not offer a complete description of the flow for multiphase problems. A sample domain is shown in Fig. 1(a). To completely describe the system, the interface location and interface matching conditions must be known. Label the liquid and gas portions of the domain using $L$ and $G$, respectively. The liquid-gas interface is then defined as the intersection of these two regions, and is referred to with $\Gamma$. On the interface, an outward facing normal, $\bm{n}_{\Gamma}$, is defined pointing away from the liquid. An vaporization mass flux, $\dot{m}$, forms normal to the interface, and removes an amount of mass $\Delta M$ from the liquid (Fig. 1(b)).

The notation $\cap$ will be used to refer to the intersection of two sets, so that $\Gamma=L\cap G$. The notation $\cup$ will be used for the union of two sets, so that the entire domain is given by $L\cup G$. In Fig. 1, the liquid region is represented using a dark color, the gas is represented using white, and the interface is represented as a thick line. This convention is used throughout the remainder of this work.

Because only two phases exist in this problem, $\Gamma$ can be identified using the liquid indicator function,

[TABLE]

The interface lies where $f$ passes between [math] and $1$, or equivalently when $\nabla f$ is nonzero. In practice, numerical schemes transport the liquid volume fraction, $\alpha$, which is the volume average of $f$ in each computational cell.

2.2 Interfacial Conditions

Consider the local velocity of interface motion, $\bm{u_{S}}$. It follows from the principle of conservation of mass that,

[TABLE]

This expresses the fact that the flux of mass into the gas from the interface is the same as the flux from the liquid into the interface. This flux is denoted $\dot{m}$. It satisfies

[TABLE]

hence, $\dot{m}$ is the local vaporization rate per unit area. (By convention, integrals of type $dS$ refer to surface integrals and of type $dV$ refer to volume integrals.) This also implies a transport equation for the liquid volume fraction in the form,

[TABLE]

which is derived from substituting the definitions of $f$ and $\alpha$ into Eq. 12 and using the Reynolds transport theorem. Here $\delta_{\Gamma}$ is the interfacial surface area density which can be defined for a region $\Omega$ as,

[TABLE]

Equation 11 can also be expressed as an interfacial compatibility condition between the velocities of the two phases. Defining a jump notation as

[TABLE]

the interface condition, Eq. 11 can be written as

[TABLE]

For the momentum equation, the interfacial condition is the jump in the pressure across the interface given by

[TABLE]

The term $\sigma\kappa$ is the pressure caused by surface tension, where $\sigma$ is the surface tension coefficient and $\kappa=\nabla\cdot\bm{n}_{\Gamma}$ is twice the mean curvature.

The interfacial conditions for the remaining equations are more complex. The gradients of the vapor mass fraction and temperature are related implicitly through the definition of $\dot{m}$. The value of $\dot{m}$ varies locally on the interface, and its value reflects the local imbalance in thermal energy and species mass across the interface,

[TABLE]

and,

[TABLE]

which are derived form Eqs. 7 and 3, respectively. Since both Eqs. 18 and 19 must be satisfied simultaneously, it follows that the gradients of vapor mass fraction and temperature are not independent at the interface. In Eq. 18, the specific enthalpy, $h$, in each phase is approximated as $h=h_{sat}+C_{p}\left(T-T_{sat}\right)$. This leads to $\left[h\right]=L_{V}+\left[C_{p}\right]\left(T_{\Gamma}-T_{sat}\right)$, where $L_{V}$ is the specific latent heat of boiling. Accordingly, $\left[h\right]$ is an effective latent heat, which varies as a function of temperature. Because Eq. 18 is derived from the thermal energy equation, Eq. 7, it implicitly presumes that thermal energy transfer at the interface dominates mechanical energy transfer. However, this presumption is consistent with the thermodynamic assumptions used in the work.

The relationship between $T$ and $Y$ is further complicated by the requirement of the thermodynamic equilibrium at the interface between the fuel and its vapor. This implies a functional relationship between the vapor mass fraction and the temperature fields. This relationship can be described by the Clausius-Clapeyron relations,

[TABLE]

Here $X$ is the vapor mole fraction and $W_{V}$ and $W_{A}$ are the vapor and ambient gas molar masses, respectively. Equation 20 is derived assuming the gas behaves ideally, hence the use of the gas constant, $R$. The saturation temperature, $T_{sat}$ are assumed to be known and are fixed in time and space. The temperature across the interface is assumed to be continuous, i.e., $\left.T_{L}\right|_{\Gamma}=\left.T_{G}\right|_{\Gamma}=T_{\Gamma}$. Equations 18, 19, 20a and 20b must be satisfied simultaneously to correctly resolve the vaporization.

3 Numerical Approach

This section details the algorithms used to perform direct numerical simulations of vaporizing liquid-gas flows. The algorithms are embedded into NGA, a computational code for solving low-Mach, multiphase flows [2]. NGA uses a mixture of finite volume and finite difference schemes to solve the conservation equations (Section 2) on a staggered grid. Two time integration schemes are used which are detailed in Section 3.2. A detailed exposition of the methods used to solve the vaporization problem follows.

The numerical methods used in this work will be detailed in semi-discrete form. For time evolving fields, superscripts are used. When necessary, secondary superscripts are used to indicate sub-iterations, e.g., $f^{n,k}$ occurs at time level $n$ at the $k$ sub-iteration. By convention, $f^{n+1,0}=f^{n,m}$ where $m$ is the total number of sub-iteration used in each timestep. The elapsed time between time levels is represented as $\Delta t$. For spatially varying fields, subscripts are used to identify the position, e.g., $f_{i,j,k}^{n}\equiv f(\bm{x},t)$, where $\bm{x}$ is the centroid of the grid cell indexed with $(i,j,k)$ and $t=n\Delta t$. The distance between the grid cell centroids $(i,j,k)$ and $(i+1,j,k)$ is given by $\Delta x$. The distances $\Delta y$ and $\Delta z$ are defined analogously. For brevity, superscript and subscript notations will be omitted when the reference is unambiguous.

3.1 Sharp Representation of Liquid-Gas Flows

Before discussing the numerical methods themselves, it is worthwhile to discuss the strategy used to represent physical quantities and field variables in each grid cell. Following, Anumolu and Trujillo [18], a strategy is said to be sharp when

the interface position is known, 2. 2.

the thermophysical properties at the interface are discontinuous, and 3. 3.

fluxes are evaluated using only information from one phase.

Item 1 requires the use of an interface tracking scheme with explicit interface reconstruction. Such a scheme is used in this work (detailed in Section 3.7). With precise knowledge of the interface position, Items 2 and 3 reduce to a statement on how information is averaged between the two phases. This work uses a mixture of both sharp and non-sharp representations.

A two-phase momentum equation is formed by taking the volume average of the Navier-Stokes equations, Eq. 1, evaluated in each phase. Using $\alpha$, the effective density and viscosity in a computational cell can be defined as

[TABLE]

and

[TABLE]

The form of the momentum equations remains the same, simply replacing the effective density and viscosity for their single phase counterparts. A two-phase continuity equation is constructed similarly. Throughout the remainder of this work, the effective density and viscosity will be referred to using $\rho$ and $\mu$ while the phase-specific values will always use subscripts $L$ and $G$. The resulting representation of $\rho$, $\mu$, and $\bm{u}$ is not sharp, since values from both sides of the interface are used in the computation of these quantities. Even so, this approach is preferred as it is consistent with the finite volume methodology used in this work.

In contrast, the temperature is represented using two fields, $T_{L}$ and $T_{G}$, which are defined within the liquid and gas, respectively. This approach follows Ma and Bothe [5] who demonstrated that such a representation produced more accurate temperature profiles for vaporization simulations. This representation of the temperature field is sharp, since it uses only information from one side of the interface. The vapor mass fraction field is defined only in the gas, hence it is also treated sharply. For grid cells that contain both liquid and gas, data are computed at the phase barycenters $\bm{x}_{L}$ and $\bm{x}_{G}$, respectively. For example, in the cell $(i)$ the value of $T_{L}$ is computed at $\bm{x}_{L,i}$ while the value of $T_{G}$ is computed at $\bm{x}_{G,i}$.

3.2 Time Integration

In general, the transport equations used in this work can be rewritten for a variable $\phi$ in the form

[TABLE]

where $\xi$ is some appropriately defined diffusion coefficient and $S$ represents sources to the equation. Selecting a time integration involves finding appropriate temporal discretization of these operators.

Two time integration schemes are used in this work. The default scheme is an iterative semi-implicit scheme given in [19] which is second order accurate. This time integration scheme corresponds to solving the semi-discrete equation

[TABLE]

with an arbitrary discretization of the spatial operators. It can be proven that this scheme is linearly stable [19]. The second discretization scheme used in this work is a

[TABLE]

The major difference between the two scheme is the time discretization of the diffusive fluxes, which are handled fully implicitly in the later expression. The reason for this is as follows. Consider a simplified transport equation consisting only of the diffusive component of the original equation,

[TABLE]

which has been taken in one dimension, for ease of exposition. (Equation 28 exploits the constant density assumption by dividing density from both sides of the equation.) When discretized temporally using Eq. 26 and spatially using the standard second order finite volume discretization

[TABLE]

this yields,

[TABLE]

where $C=\xi\Delta t/{\Delta x}^{2}$. This equation corresponds to a monotone scheme (see the proof in Section 10.1) and has the property of being not only linearly stable, but also non-oscillatory. Since vaporization is controlled by diffusive mechanisms via Eqs. 18 and 19, the non-oscillatory solution of diffusion is important. This idea is further explored in Section 5.2.

3.3 Momentum Equation

The momentum equation is solved using a mass-momentum consistent advection strategy [3]. The scheme is known to be stable for flows with large density ratios such as those encountered in a fuel-air system. The basic idea of the scheme is to update both $\rho$ and $\rho\bm{u}$ in time such that their discrete solution implies the discrete transport equation for $\bm{u}$, via $\bm{u}^{n+1}=\left(\rho\bm{u}\right)^{n+1}/\rho^{n+1}$. A overview of the solution is shown below for a one dimensional system. Due to the staggered grid approach, velocities are solved on a mesh that is offset from the primary mesh (Fig. 2). A new set of effective densities and viscosities, $\rho_{u}$ and $\mu_{u}$, must be defined using $\alpha_{u}$, the local fraction of liquid volume in the $u$ velocity cell. This local volume fraction is computed by constructing the $u$ velocity cell, and computing the liquid volume contained therein (Fig. 2). An interface reconstruction scheme must be used to compute the exact volume inside of the $u$ velocity cell. The scheme used in this work is discussed in Section 3.7. At each time step, $\alpha^{n}_{u}$ is computed from the cell centered field, $\alpha^{n}$ at the current time level. The resulting algorithm to update $u$ is

[TABLE]

where $RHS$ represents all non-advective terms in the Navier-Stokes equations. First Eq. 31b is solved for $\rho^{n+1,k+1}_{u}$, and that value is used in Eq. 31c to update the velocity. The key aspect of this treatment is that the same discrete representation of the flux operator $\partial/{\partial x}$ is used for advection in Eqs. 31b and 31c. Although any discretization may be used for this operator in general, care must be taken to ensure the stability of the scheme. In this work the flux operator is discretized using first order upwind in interface containing cells, and a second order finite volume discretization[2] otherwise. The specifics of the discretization including the definition of $RHS$ are given in [3].

3.4 Pressure Poisson Equation

For low-Mach flows, the continuity equation is enforced through the pressure. In this work, Chorin splitting [20] is used to couple the momentum and continuity equations. First, an intermediate velocity is computed by solving the Navier-Stokes equations, Eq. 1, with zero pressure. This field is referred to as ${\bm{u}}^{NS}$, the Navier-Stokes velocity. Next the velocity field is updated with a pressure projection according to

[TABLE]

Taking the divergence of Eq. 32 yields the Poisson equation for pressure,

[TABLE]

Conceptually, Equations 32 and 33 are similar to the projection technique used in purely incompressible flows. In that context, the velocity is projected onto its divergence-free part using a Helmholtz decomposition. However, in this work, it is not assumed that $\bm{u}^{n+1}$ is incompressible, so the right hand side of Eq. 33 is nontrivial. By rearranging Eq. 2, we see

[TABLE]

and this value is used in Eq. 33. For brevity, the notation

[TABLE]

is adopted.

Here it should be noted the abuse of notation used in writing operators such as $\nabla\cdot\bm{u}$ in the preceding discussions. The flow solver solves the equations of motion in the finite volume sense, meaning that the flow variables are updated in terms of fluxes through control volumes. This approach is consistent with directly solving the integral form of the conservation equations. However, in this work flow variables are not necessarily differentiable, requiring a careful definition of the various operators. As such, operators must be understood in the following sense. For some vector, $\bm{f}$, define the operator,

[TABLE]

which is the net flux per unit volume out of $\Omega$. For sufficiently smooth vector fields this limit yields

[TABLE]

which follows from the Gauss-Ostrogradsky (Divergence) Theorem. In this work, all flow variables are continuous except at the interface, so relation Eq. 37 applies almost everywhere. Throughout this work, the suggestive notation $\nabla\cdot\bm{f}$ has been used with the meaning $\text{div}\left(f\right)$.

The form of Eq. 37 as a surface integral divided by an volume integral suggests that the surface density, $\delta_{\Gamma}$, may be important in evaluating the divergence operator, div, for discontinuous vectors. This can be seen by directly computing $\text{div}\left(\bm{u}\right)$ near the interface. Define a volume $\Omega_{\Gamma}$ by extending a small distance from the interface in both directions (i.e. in directions of both $\bm{n}_{\Gamma}$ and $-\bm{n}_{\Gamma}$). In such a volume $\text{div}\left(\bm{u}\right)_{\Gamma}$ satisfies,

[TABLE]

In deriving Eq. 38 it was important to assume that there is no jump in tangential velocity, such that the net flux through $\partial\Omega_{\Gamma}$ equals the net flux through $\Gamma$.

The overall value of $\text{div}\left(\bm{u}\right)$ in a discrete computational cell contains three components: a component from the fully gas portion of the cell, a component from the fully liquid component of the cell, and an interfacial component as detailed in Eq. 38. Using this information, and expanding $\left[\bm{n}_{\Gamma}\cdot\bm{u}\right]_{\Gamma}$ using Eq. 34, the pressure Poisson equation can be written as

[TABLE]

where Eq. 35 is applied separately in each phase, so that

[TABLE]

$H$ can change due to variations in the chemical composition, temperature, and thermodynamic pressure of the fluid. As a first step, the current work explicitly neglects changes in density due to chemical composition and temperature. This simplifies the exposition of the system thermodynamics considerably. Consistent with the low-Mach assumption, the thermodynamic pressure is assumed to vary only temporally. This implies that each phase’s value of $H$ is also a function of time only. Finally, we explicitly assume that the liquid is fully incompressible, implying that $H$ can be nonzero only in the gas.

By assumption, nonzero $H$ correspond to a build up of pressure in the computational domain. Since sources due to chemical composition and temperature are neglected, this build up is due entirely to the evaporative flow. In the presence of outflow conditions, this flow is able to freely exit the domain. Accordingly, vaporization does not cause any pressure buildup, and $H=0$. For closed boundary conditions, such as periodic boundary conditions, the flow cannot exit the domain, and a gradual rise in pressure occurs. For low-Mach flows, a consistent choice for this source is

[TABLE]

This choice ensures that the compatibility condition of the pressure Poisson equation is satisfied, i.e,

[TABLE]

Satisfaction of this condition is sufficient to guarantee existence and uniqueness of the pressure correction, $\nabla p$ [21].

The gas density is then updated as

[TABLE]

which can be considered to be a low-Mach correction to the Poisson equation that accounts for gas phase compressibility.

The pressure Poisson equation, Eq. 39, is solved using a technique based upon the Ghost Fluid Method [6]. An illustrative discretization in one dimension is defined below. The pressure Laplacian operator is discretized as,

[TABLE]

where the assumption of constant grid spacing has been used. The velocity divergence is written as

[TABLE]

The singular source term involving $\delta$ requires a special treatment. Here this term is discretized as

[TABLE]

the value $Q$ is one in interface containing cells and zero otherwise. The pressure jump, Eq. 17 is discretized using the standard GFM technique. As a single pressure field is solved for the entire domain, this discretization is not sharp, in the sense of [18]. However, this approach was chosen to be consistent with the one field approach to the velocity field. Furthermore, in practice, this technique only distributes the pressure jump over a region of one or two grid cells around the interface.

3.5 Scalar Transport Equation

Transport of the temperature and vapor mass fraction fields is handled sharply to avoid artificial smearing of the profiles. One of the necessary aspects of a sharp representation is that data for one phase are not computed using values from the other phase. For this purpose, several auxiliary variables are defined that will be necessary to discuss in the transport equations. At each time step, a gas volume fraction $\alpha_{G}^{n}=1-\alpha^{n}$ and secondary liquid volume fraction $\alpha_{L}^{n}=\alpha^{n}$ are constructed. To aid in the advection of scalars, a continuous velocity field is constructed for both the liquid and gas fields, labeled as ${\bm{u}}_{L}$ and ${\bm{u}}_{G}$, respectively. The fields are initialized using the true velocity field at the latest time level, ${\bm{u}}^{n+1,k}$, and the data are extrapolated using the zero normal gradient technique of Aslam [22]. In the technique, a PDE is solved in pseudotime, $t_{ps}$, until steady state, and the steady state solution verifies the zero normal gradient property. The equation to extend an arbitrary field, $\phi$, from the liquid into the gas can be expressed as

[TABLE]

where

[TABLE]

A similar process is used to extrapolate from gas to liquid. The technique results in an extrapolated field that is constant along the direction of $\bm{n}_{\Gamma}$. In [22], a general technique for extrapolation to arbitrary polynomial order is discussed. However, for the current work it was found by experience that constant extrapolation yields superior numerical stability.

The transport of scalars uses a variation of the mass consistent scheme [3], outlined below for the vapor mass fraction

[TABLE]

First, Eq. 51b is solved for $\rho_{G}^{n+1}\alpha_{G}^{n+1,k+1}$ and then this value is used in Eq. 51c to update $Y$. For advection, the BQUICK algorithm [23] is used. The algorithm modifies the QUICK advection routine, by locally switching to an simple upwind scheme wherever the computed $Y^{n+1}$ becomes unbounded. For the present work, a second switching condition is used. Whenever the QUICK stencil crosses the interface, the simple upwind scheme is used. To use the BQUICK algorithm, reasonable bounds must be chosen for each scalar. For the mass fraction field, the maximum and minimum values of $1$ and [math], respectively, are chosen. For the temperature fields, the maximum of the initial domain temperature, $T_{M}$, and [math] are chosen. This choice is appropriate if there are no sources of heat, in which case $T_{M}$ will always be the highest domain temperature physically possible. The overall accuracy of this scheme is first order, because the simple upwind scheme is first order accurate. In practice, higher accuracy has been observed [23].

The treatment of diffusion in interface-containing cells requires special care. Figure 3 demonstrates a region near the interface in one dimension. The interface does not occur inside either of cells $i$ and $i+1$, so the flux between the cells is computed using the second order finite volume discretization,

[TABLE]

A different treatment is used to compute the flux whenever one (or both) of the two cells contains the interface, such is the case with cells $i$ and $i-1$. In this case, the flux is computed using the relative distances between the appropriate phase barycenters. For example, the discretization of the flux for the vapor mass fraction is

[TABLE]

Here $x_{i}$ is the $x$ axis component of $\bm{x}_{G,i}$. The value $\beta_{i-1/2}$ is the fraction of gas area on the face between cells $i$ and $i-1$. The value $D_{i-1/2}$ is linearly interpolated between $D_{i}$ and $D_{i-1}$. Note that the VOF scheme does not represent the interface smoothly between cells, and accordingly, the gas area fraction may be different on either side of a cell face. In such a case, $\beta_{i-1/2}$ is taken to be the minimum of the values.

In cell $i-1$, an additional diffusive flux from the interface into the cell must be computed which satisfies

[TABLE]

or in terms of the surface area density,

[TABLE]

The evaluation of the flux at the interface is the focus of Section 3.6. Although not detailed here $T_{G}$ and $T_{L}$ are solved in a manner analogous to the vapor mass fraction field.

3.6 Thermodynamic Equilibrium at the Interface

Vaporization involves the flux of mass, momentum, and energy across the interface. The imbalance of these fluxes is directly tied to the definition of $\dot{m}$ in Eqs. 17, 16, 18 and 19. This relation is further complicated by the requirement of the system to satisfy thermodynamic equilibrium, i.e., Eq. 20. Because of this complexity, the solution of this nonlinear system is often neglected in the literature by assuming a priori a limiting vaporization regime such as boiling conditions. There are, however, a few strategies in the literature capable of solving the general problem.

The first strategy used in this work is similar to a technique used in Ma and Bothe [5]. In their work, they assume an error function profile for the temperature field, and substitute the thermodynamic relations to yield a nonlinear system that can be solved for the interface temperature. Such an assumption is appropriate for low $\dot{m}$ vaporization (see Section 4.1), but fails to be appropriate for situations with appreciable convection. Instead, the current work constructs the system directly from the discretization of the governing equations. Recalling that gas and liquid temperatures (as well as mass fraction) are stored in their respective phase barycenter, specialized computational stencils must be use to accurately compute the gradient if the stencil includes the interface. Although several discretizations are possible, a general form of such a stencil is

[TABLE]

where the summation is taken over computational cells near the interface and is restricted to the gas side. A liquid stencil can be defined similarly. Substituting these expressions into Eq. 18 yields

[TABLE]

where $w_{G,\Gamma}$ and $w_{L,\Gamma}$ refer to the gas-side and liquid-side coefficient to the interface-containing cell. Notably, knowledge of $\dot{m}$ and the current temperature field allows the expression of the interface temperature in the form $T_{\Gamma}={\mathcal{T}}(\dot{m})$.

The system is closed by noting that discretizing Eq. 19 yields a relationship of the form $\dot{m}=\dot{\mathcal{M}}(Y_{\Gamma})$, and the Clausius-Clapeyron relation, Eq. 20, requires $Y_{\Gamma}=\mathcal{Y}(T_{\Gamma})$. The entire system may thus be written as

[TABLE]

The solution of Equation 58 can be solved numerically using any root finding routine. A simple regula falsi solver was used in this work. The method outlined above solves the system with discrete consistency in the sense that each equation is satisfied simultaneously in each grid cell. This method for computing interface equilibrium will be henceforth referred to as MB method.

Rueda-Villegas et al. [1] forgo the solution of the nonlinear problem, Eq. 58. Instead, they use the approximation $Y_{\Gamma}\approx Y(\bm{x})$, and then invert the Clausius-Clapeyron relation to find $T_{\Gamma}$. This method is not discretely consistent. To correct for this, they do not directly compute $\bm{n}_{\Gamma}\cdot\nabla Y$ on the interface and instead substitute $\frac{1-Y_{\Gamma}}{\rho_{G}D}\dot{m}$ where this quantity is needed. Note that the temperature field is used to compute $\dot{m}$ through Equation 18. This method is abbreviated as the REA method. The use of the local value of $Y$ for the interface value introduces a first order error into the system. In the current work, we introduce a modified version of the REA technique that instead extrapolates the $Y$ field to the interface with second order accuracy. The importance of this higher order accuracy will be discussed in Section 4.1. This modified technique is referred to as REA2.

In this work, two stencils are used to discretize the vaporization fluxes at the interface. The pressure Poisson equation, Eq. 33, and the VOF transport equation, Eq. 13, use $\dot{m}$ but have no direct knowledge of the scalar fields. The scalar transport equations, in turn, are effectively independent of each other, only being coupled implicitly through the boundary conditions, $T_{\Gamma}$ and $Y_{\Gamma}$. This leads to the following splitting. Inside of scalar transport equations, the vaporization fluxes are implemented using a first order accurate gradient approximation. For the gas temperature, this takes the form,

[TABLE]

where $\bm{x}_{G}$ and $\bm{x}_{\Gamma}$ are the gas phase barycenter and interface barycenter within the computational cell. Note that Eq. 59 is discretized using the $n+1$ time level, so that its use corresponds to a fully implicit time integration. Furthermore, this discretization is monotone and therefore does not introduce oscillations to the solution of the scalar fields. This discretization is important. Since the interface moves freely through the domain, $\|\bm{x}_{G}-\bm{x}_{\Gamma}\|$ can approach zero making the interfacial flux term be of arbitrary magnitude. For interface containing cells, this term usually dominates the scalar transport phenomena, so its accurate resolution is paramount.

For the pressure Poisson equation and VOF transport equation it was found that a higher order treatment of $\dot{m}$ was necessary to achieve a convergent solution. First, each scalar is extended across the interface using linear extrapolation [22]. Then a least squares fit to the temperature profile in each phase is computed over each cell and its $5^{dim}-1$ neighbors, in $dim=3$ dimensions. Gradients are calculated from the least square fit, and finally Eq. 18 is used to compute $\dot{m}$. For example, in one dimension the cell $\Omega_{k}$ is fit with the functional form

[TABLE]

yielding the system

[TABLE]

which is solved in the least squares sense, i.e., $\bm{c}=({\bm{\mathcal{L}}}^{\top}{\bm{\mathcal{L}}})^{-1}{\bm{\mathcal{L}}}^{\top}\bm{T}$. It follows from Eq. 60 that the coefficients $c_{i}$ are the least squares fits for the temperature profile at the cell centroid. In practice, this technique does not lead to any issues with stability, due to the weak coupling of Eqs. 33 and 13 with the scalar transport equations. The technique provides another interesting property. As noted in [18], knowledge of $\dot{m}$ is needed not only in interface-containing cells, but also in a narrow band of cells around the interface. The treatment of $\dot{m}$ given in our work provides a natural framework to extend $\dot{m}$ across into such a band, namely, by computing Eq. 18 locally in each cell using the extended temperature field.

3.7 Interface Transport and Divergence-free Velocity Extrapolation

Interface transport is accomplished using a second order accurate, unsplit, geometric volume of fluid method [4]. The volume fraction is transported using the extrapolated liquid velocity, ${\bm{u}}_{L}$. The VOF method used in this work ensures consistency between the discrete velocity divergence and the flux volumes used in the transport step. An implication of this feature is that errors in divergence of the extrapolated velocity field lead to erroneous variations in the liquid volume. To remove this source of error, the extrapolation is completed in two steps. First, the velocity is extrapolated using constant extrapolation [22]. The resulting velocity is labeled ${\bm{u}}_{L}^{*}$. Next, ${\bm{u}}_{L}^{*}$ is projected onto its divergence-free part using

[TABLE]

where W is a potential derived from the Helmholtz-type equation,

[TABLE]

The coefficient $a$ is zero for all liquid-containing cells and cells that are within three grid cells of the interface, which reduces it to a Poisson equation in this region. For all other cells, $a$ is set to an arbitrary constant, here given by $\frac{1}{{\Delta t}^{2}}$. The advantage of solving Eq. 63 over a Poisson equation is that the given technique does not need to satisfy the compatibility relation

[TABLE]

Equation 64 will not hold, in general, for the extrapolated velocity field.

The volume of fluid advection equation, Eq. 13, can be solved in a modified form. The interface recession term, $\dot{m}\left/\rho_{L}\delta_{\Gamma}\right.$, is rewritten in terms of a recession velocity, $\dot{m}/\rho_{L}\bm{n}_{\Gamma}$. Furthermore, since $\nabla\cdot{\bm{u}}_{L}=0$ everywhere that $\alpha$ is nonzero, the resulting equation simplifies to

[TABLE]

In practice, the velocity $\left({\bm{u}}_{L}-\frac{\dot{m}}{\rho_{L}}\bm{n}_{\Gamma}\right)$ is used as the velocity seen by the VOF solver.

The volume fraction field is used to compute all relevant geometric information for the solver. The interface is reconstructed from $\alpha$ using the Piecewise Linear Interface Calculation (PLIC) method with ELVIRA [24] which also defines the interface normal $\bm{n}_{\Gamma}$. PLIC represents the interface in each cell as a simple polygon which separates a liquid polyhedron from a gas polyhedron. The representation of the geometry in terms of polytopes allows for a simple algorithm to compute the area (volume) of each polygon (polyhedron) by subdividing it into triangles (tetrahedra) [4]. The curvature $\kappa$ is computed using the mesh decoupled height function approach [25]. A diagram showing the flow of a single time step of the solver is shown in Fig. 4.

4 Verification in One Dimension

A suite of numerical tests is performed to study the accuracy of the code. First, spatial convergence is demonstrated against a well-studied heat transfer limited case. Second, this case is modified to study vaporization without the constraint of heat transfer limitation. All simulations are required to satisfy a CFL stability restraint, which is a numerical limitation on the possible time step that can be used in a simulation. The CFL numbers depend on the physical problem being studied as well as the numerical discretization of the governing equation. The CFL numbers for advection, viscous diffusion, mass diffusion, temperature diffusion, and surface tension are

[TABLE]

where the notations $u_{M}=\sum\limits_{i\in 1,2,3}|u_{i}|$and $\Delta_{m}=\min\left(\Delta x,\Delta y,\Delta z\right)$ are used. The thermal diffusivity, $\lambda=k/(\rho C_{p})$, has been introduced for convenience. Equations 67 and 69 should be understood to apply to both the liquid and the gas values of $\mu$ and $\lambda$. For the numerical schemes used in this work, the conditions $C_{\sigma}<1$ and $C_{u}<1$ must be satisfied, whereas the relations Eqs. 67, 68 and 69 are relaxed due to the use of implicit solvers.

4.1 Spatial Convergence

Due to the inherent complexities of the vaporizing problem, there are very few analytical solutions to verify against. One problem that does have such a solution is Neumann’s problem, whose solution is given in the paper by Hardt and Wondra [17] The vaporization problem is solved in one dimension with an assumption of a slow moving interface and negligible convective effects. A schematic is shown in Fig. 5. A region of hot gas heats the liquid, causing vaporization. Dirichlet conditions are specified for the temperature and mass fraction of the gas at the left boundary whereas the liquid is able to flow freely from the right boundary. The domain is $1$ mm in length.

Under the assumptions given here, the energy equation, Eq. 7, reduces to the unsteady diffusion equation with time varying boundary. Define the initial and current interface location as $x_{0}$ and $x_{\Gamma}$, respectively. The solution is given in terms of the error function as

[TABLE]

In Eq. 71, the interface position has been assumed to vary as $x_{\Gamma}=x_{0}+2\ell\sqrt{\lambda t}$. The diffusion layer value, $\ell$, must be solved from thermodynamic variables as

[TABLE]

In their work, Hardt and Wondra solve only for the gas temperature field, however, the same analysis can be performed for the mass fraction field, yielding

[TABLE]

In general, an equation similar to equation Eq. 72 can be written for $\ell_{Y}$, however for the unity Lewis number case, the uniqueness of the interface position implies $\ell_{Y}=\ell$. By assumption, the liquid temperature field remains at the interface temperature, $T_{\Gamma}$. Equations 71 and 73 can be combined with the condition of equality of the two expression for $\dot{m}$, Eqs. 18 and 19, to yield an analytical relation between the interface temperature and mass fraction,

[TABLE]

Equation 74 should be combined with the Clausius-Clapeyron relation, Eq. 20 to close the system.

The parameters used in this study mimic the properties of a vaporizing water-air system. They are listed in Table 1. The fluid and thermal properties listed are the same as those used in the work of Hardt and Wondra [17]. A unity Lewis number is imposed, and for consistency with their work, the vapor was given the molar mass of water. Two test cases are performed with these parameters by varying $T_{\infty}$ and $Y_{\infty}$.

Case 1: Film Boiling

The first test, Case 1, corresponds exactly to that used in the paper of Hardt and Wondra, which uses $T_{\infty}=383.15\text{ K}$ and $Y_{\infty}=1$, leading to $T_{\Gamma}=373.15\text{ K}$ and $Y_{\Gamma}=1$. To aid in the imposition of the boundary conditions, $x_{0}$ is chosen to be inside the computational domain at $x_{0}=0.02L=2\times 10^{-5}\text{ m}$. We chose to initialize the simulation with the analytical solution at time $t=0.1\text{ s}$. Simulations are performed of vaporization using the given parameters between times $t=0.\text{ s}$ and $t=0.2\text{ s}$. We initialize the liquid field to the interface temperature, $T_{\Gamma}$. Both the MB technique and REA technique are used to resolve interface equilibrium. A grid convergence study is performed using $50$, $75$, $100$, $150$, and $200$ grid cells in the domain. All simulations are performed for a fixed $C_{\lambda}=C_{D}=2$. The analytical solution at $t=0.3\text{ s}$ is shown in Fig. 5 for the 50 cell grid. It can be seen that both the MB and REA simulations match well with the expected solution. We define error norms for the gas temperature as

[TABLE]

The vapor mass fraction and liquid temperature fields are omitted, because they remain within machine precision of their initial values for all cases, as expected. The $L_{1}$ and $L_{\infty}$ errors in gas temperature are plotted in Figs. 6(a) and 6(b). It can be seen that approximately first order convergence is noted, which is consistent with the overall scheme accuracy. Furthermore, the two schemes produce almost identical results.

Case 2: General Vaporization

While valuable as a comparison to previous work, the standard Neumann’s problem does not test the accuracy of the dynamic interface computation, in a situation with nontrivial mass fraction field. A second case, Case 2, was performed using the same grid and fluid parameters with $T_{\infty}=323.15\text{ K}$ and $Y_{\infty}=0.2$. This yields $T_{\Gamma}=296.163\text{ K}$ and $Y_{\Gamma}=0.221022$. The simulation was performed between times $t=0.01$ and $0.1\text{ s}$. A grid convergence study is performed with $C_{\lambda}=C_{D}=2$. For a more complete analysis of this field, additional error measurements are taken. Errors for the vapor mass fraction, the liquid temperature, and liquid velocity are defined in manner analogously to Eqs. 76 and 77. The gas temperature and vapor mass fraction are referenced to the analytical solution, Eqs. 71 and 73. The liquid temperature is referenced to initial value, which is predicted to be constant in time. The velocity error is computed relative to the following solution. Since the problem is one dimensional, the continuity equation predicts a piecewise constant velocity, with a discontinuity only at the interface. The liquid velocity can then be predicted from the definition of $\dot{m}$

[TABLE]

where $u_{L}=-\bm{n}_{\Gamma}\cdot\bm{u}_{L}$, $u_{G}=-\bm{n}_{\Gamma}\cdot\bm{u}_{G}$, $u_{S}=-\bm{n}_{\Gamma}\cdot\bm{u}_{S}$. Note that Eq. 78 has been expressed such that positive speeds corresponds to motion in the positive $x$ direction. This expression rearranges into,

[TABLE]

Since $x_{\Gamma}=x_{0}+2\ell\sqrt{\lambda t}$, it follows

[TABLE]

With this information, the liquid velocity can be expressed as

[TABLE]

where $u_{G}=0$ has been substituted.

The gas temperature converges with first order accuracy in both $L_{1}$ and $L_{\infty}$ sense for all three techniques (Figs. 7(a) and 7(b)). However, it is interesting to note that errors using REA are consistently higher than either of the other two techniques, and errors using REA2 are consistently the lowest. A similar trend is seen for the vapor mass fraction field in Figs. 7(c) and 7(d), although greater variance about the trend is seen for the REA technique. The liquid temperature field (Figs. 7(e) and 7(f)) demonstrates a different trend. While first order convergence is seen using the REA and REA2 techniques, no convergence is seen using MB. Indeed, MB appears to divergence under grid refinement. However, the overall liquid temperature error associated with the MB technique is orders of magnitude lower than that of the REA or REA2 techniques. The liquid velocity error is demonstrated in Fig. 7(g). Since the liquid velocity is spatially constant, there is no distinction between $L_{1}$ and $L_{\infty}$ errors. This property has been verified, but the redundant plot is not shown here. A trend of approximately linear convergence is displayed, although the trend appears to have some variation.

5 Vaporization of a Curved Interface: Convergence and Stability

The solutions to the vaporization problems presented so far have been analytical solutions in one dimension. These solutions, while interesting for verification, do not allow for a more full test of the capabilities of the solver. One of the goals of this work was to develop a robust numerical solver for simulating vaporizing multiphase flows. However solutions in 1D cannot fully test this. Many of the most challenging aspects of multiphase flow simulation are due to non-alignment between the mesh and flow structures; however, in one dimension all flow structures are aligned with the mesh. To understand what influence, if any, mesh alignment has on the flow solver simulations are performed in two dimensions with a complex, curved interface. Since no analytical solution is known to this problem, errors are approximated by comparing simulation data to those of the same simulation run on a fine grid. This section also serves to test the numerical stability of the solver, as a temporal convergence study is performed. Results are compared for both the semi-implicit and fully implicit numerical time integrator.

The interface shape and initial conditions are depicted in Fig. 8. The interface shape is parameterized as a function of the $y$ position and is given by

[TABLE]

where $w_{h}$, $w_{a}$, and $w_{n}$ are arbitrary constants here chosen as $\frac{2}{3}$, $\frac{1}{6}$, and $5$, respectively. The grid length in the $x$ and $y$ directions are represented by $L_{x}$ and $L_{y}$, respectively and are equal $L_{x}=L_{y}=1e-3$. The initial conditions used in this problem are chosen arbitrarily, but have been inspired by those of the previous section. The condition may be written as

[TABLE]

The fluid properties are the same as the previous cases and are listed in Table 1.

5.1 Solution Convergence

Simulations are performed for grids of size $25\times 25$, $50\times 50$, $75\times 75$, and $100\times 100$, and the data are compared to a simulation with grid size $200\times 200$ for accuracy. Simulation data are taken at time $t=0.05s$. Both the MB and REA2 techniques have been used for this study, and both technique perform well for grid convergence. Images of the convergence of the interface position are shown in Fig. 9. The profile seems to sharpen under grid refinement. The center most hump appears taller, whereas the outer humps appear to decrease in height consistent with the principle of mass conservation. The interface converges rapidly towards the fine mesh profile. Figure 10 depicts the convergence of the gas temperature, vapor mass fraction, liquid temperature, and liquid velocity in the domain. A first order scaling is recovered for both $L_{1}$ and $L_{\infty}$ errors in the gas fields. For the liquid temperature field, first order convergence is seen in $L_{1}$, but slower convergence is noted for $L_{\infty}$. The liquid velocity field appears to converge at a very slow rate in $L_{1}$, but diverges in $L_{\infty}$. This is not entirely surprising. From the previous section, we noted that the liquid velocity is directly proportional to $\dot{m}$ and therefore directly proportional to the scalar gradient at the interface. Since the scalar field is itself only first order accurate, its gradient would be expected to lose an order of accuracy, recovering only zeroth order accuracy. This explains the lack of convergence in the liquid velocity field. Since the velocity controls the liquid motion and the rate of heat advection into and out of the interface, it is reasonable that this would have an impact on the accuracy in the liquid temperature field as well. However, since the gas is stationary, only the boundary effect is felt by the gas, decreasing the impact of the liquid velocity on its accuracy.

It is important to note the significance of the trends in Fig. 10. First, it can be seen that the effective accuracy of the simulation in two dimensions is somewhat lower than what was recovered in one dimension. This is unsurprising because of the greater geometric complexity of the interface in two dimensions. The resolution of this complex interface shape becomes a limiting factor in the resolution of the fluid flow. Second, it can be seen that the error effect is particularly strong when measured in the $L_{\infty}$ sense. It is common within the literature to report errors in terms of smooth error norms such as $L_{1}$. This form of measurement can be misleading when numerical errors are not distributed uniformly, and that is the case for these simulations where errors are largest near the interface. This calls into question the use of smooth error norms for numerical analysis of flows with sharp features or discontinuities.

While the current work fails to convergence in $L_{\infty}$, it does converge in the more traditional $L_{1}$ sense. Consistent with the trends in the literature, the remainder of this paper will assume that the recovered $L_{1}$ convergence is sufficient to perform exploratory studies of vaporizing multiphase flows. However, future work will focus on improving the accuracy of the flow solver, so that improved convergence will can be recovered.

5.2 Temporal Stability and Convergence

The temporal convergence and stability of the solver is tested by performing simulations on the $100\times 100$ grid using various time steps. The time steps are chosen to yield $C_{D}=C_{\lambda}\in\{10,25,50,100,200\}$. The two time integration schemes are compared using both the MB and REA2 techniques. The data are summarized in Tables 2 and 3.

Using the semi-implicit time integrator, similar trends are seen in errors for gas temperature and vapor mass fraction using both MB and REA2 techniques. Both $L_{1}$ and $L_{\infty}$ errors decrease with $C_{\lambda}$. There is a noticeable difference in the error of the liquid temperature and velocity, as has been seen in the previous sections. It is interesting to note that neither of these field appears to converge under temporal refinement (decreasing $C_{\lambda}$) according to the $L_{1}$ metric. Surprisingly, the liquid velocity converges in both cases under the $L_{\infty}$ metric. It should also be noted that for the REA2 technique, there is spike in solution error at $C_{\lambda}=200$. It will later be demonstrated that this is due to an issue with solution monotonicity.

Using the fully implicit time integrator a few trends can be noted. The MB technique seems to recover superior temporal convergence using the fully implicit solver than the semi-implicit solver. Ignoring the $C_{\lambda}=200$ case of REA2 technique, their errors are largely the same using both time integrators. This is not a surprise. While the semi-implicit solver has higher formal accuracy than the fully implicitly integrator, this property is no use when the leading error term in the flow is dominated by the first order treatment of the interfacial diffusion flux Section 3.6.

5.2.1 Monotonicity and Stability

The data suggested both MB and REA2 schemes recovered better stability and convergence properties using the fully implicit time integration scheme. This follows from the fact that this scheme is a monotone scheme, and therefore disallows the development of numerical oscillation developing in the flow. Figure 13 demonstrates snapshots of the vapor mass fraction and $\dot{m}$ fields at simulation end using the REA2 technique and $C_{\lambda}=200$. Using the semi-implicit time integration a region of high vapor mass fraction appears near the interface leading to a non-monotonic variation of the field. This is not seen in the fully implicit time integrator. As a result, the $\dot{m}$ fields look very different. The fully implicit scheme demonstrates smooth variation of $\dot{m}$ , whereas the semi-implicit scheme demonstrates chunking. The effect can also be seen by monitoring the temporal behavior of $\dot{m}$ in Fig. 14(a). While $\dot{m}$ varies smoothly for the fully implicit scheme, it is oscillatory for the semi-implicit scheme. This translates to rapid changes in the slope of the normalized volume (Fig. 14(b)), and the eventual divergence of the solution.

6 Vaporization of an Isolated Spherical Droplet

The previous solutions, while interesting for verification, have little practical value. A somewhat more physically relevant problem is the problem of vaporizing a spherical droplet in a quiescent flow. Assuming the problem is quasi-steady (temporal derivatives of the scalars and momentum are negligible) and the gas obeys the continuity equation, $\nabla\cdot\bm{u}=0$, the solution to the temperature and mass fraction fields can be found analytically, and is given in any text on the subject [11]. The remarkable aspect of the solution is that the analytical result predicts that the squared diameter of the droplet ($d^{2}$) varies linearly in time. This result is also observed experimentally for both vaporization and combustion [7].

We follow the derivation of Rueda-Villegas et al. [1] who obtain the solution

[TABLE]

In Eqs. 86 and 87, $Q$ is a constant related to $\dot{m}$. The values $\dot{m}$, $Y_{\Gamma}$, and $T_{\Gamma}$ are related through

[TABLE]

which is derived form Eqs. 86 and 87. Equation 88 is more commonly written in the form

[TABLE]

where $B_{M}=\frac{Y_{\Gamma}-Y_{\infty}}{1-Y_{\Gamma}}$ is the mass transfer number. This relation can also be used to compute the vaporization rate per unit area, yielding

[TABLE]

Note the explicit dependence of the solution on the current droplet diameter, $d$. Equations 88 and 89 must be solved simultaneously with the Clausius-Clapeyron relation, Eq. 20, to complete the system.

In this section, parameters are chosen to imitate the properties of acetone. It was pointed out by Raessi et al. [26], that for practical simulation of liquid vaporization using real physical properties, the condition $C_{\sigma}<1$ often becomes the limiting CFL condition. Since surface tension is not important to the mechanics of the $d^{2}$ law, the surface tension was reduced to zero, allowing for the use of a larger numerical time step. Similarly, the use of the monotone, fully implicit scalar diffusion solver frees the simulation from time scales $C_{D}$, $C_{\lambda,L}$, and $C_{\lambda,G}$. Only the advection CFL condition need be satisfied, here forced to be $C_{u}<1$ for all simulations.

Some discussion must be devoted to the implementation of the boundary conditions. The simulation is run on a periodic, cubic domain using a Cartesian solver in three dimensions. The spherical symmetry implied by the far field boundary condition is approximated by centering the domain on the droplet center of mass. Tests are performed to characterize the influence of the boundary condition on the evaporation rate (see Section 10.2). It was found that at sufficiently early times in the simulation, boundary effects were negligible.

A further modification is required to deal with the boundary conditions. Experience has shown that small errors in the liquid velocity field can lead to motion of the droplet away from the center of the domain, breaking the implied symmetry of the boundary conditions. To remove this source of error from this suite of tests, $-\bm{n}_{\Gamma}\dot{m}$ is used as the velocity for VOF transport. This effectively enforces the condition, $\bm{u}_{L}=\bm{0}$. However, the reconstructed $\bm{u}_{L}$ is still used (rather than $\bm{0}$) for the liquid temperature advection.

The chosen fluid properties are listed in Table 4. The actual fluid surface tension is listed, even though the value of $\sigma=0$ is used in these simulations. The far field temperature is chosen to be $700\text{ K}$ and the far field vapor mass fraction is [math]. This results in interface temperature and vapor mass fraction of $294.92\text{ K}$ and $0.43993$, respectively. The vapor mass fraction and gas temperature fields are initialized with the analytical solution, Eqs. 86 and 87, and the liquid temperature field is initialized simply as $T_{\Gamma}$. Note that in contrast to the simulations of previous sections, this section demonstrates a large range of values in both the temperature and vapor mass fraction fields. The initialization is demonstrated in Figs. 15(a) and 15(b).

A grid convergence study was performed using a domain of size $8\times 10^{-4}$m which corresponds to exactly eight droplet diameters. The convergence study was performed using $N=32$, $64$, $128$, and $256$ grid cells in each of the three grid directions. This corresponds to resolutions of $4$, $8$, $16$, and $32$ cells per droplet diameter. The study focuses on the prediction of the $d^{2}$ law vaporization time constant, because this parameter is one of the few values that is easily measured in both experimental and computational studies. Following Sirignano [11] we write

[TABLE]

where $d_{0}$ is the initial droplet diameter, and $\tau$ is the vaporization time constant. From Eq. 92 it can be seen that $\tau$ is the time necessary to fully vaporize a droplet. This value can be computed as

[TABLE]

Using the parameters of this study, the time constant is computed to be $\tau=0.029126\text{ s}$ [11]. Simulations are run until time $t=1.5\times 10^{-3}\text{ s}$ when the ratio $d^{2}/d_{0}^{2}$ is computed, and Eq. 92 is used to extract a time constant. This choice of simulation end time corresponds to $t/\tau\approx 0.01$, which is before boundary effects are important.

Figure 16(a) shows the temporal variation of the effective droplet square diameter, computed as

[TABLE]

where $V_{0}$ and $V$ are the initial and current droplet volumes, respectively. It can be seen that a linear trend in $d^{2}$ is recovered for all cases, however, the slope of the line depends on the mesh resolution. As the grid is refined, $\tau$ appears to approach its expected value as shown in Fig. 16(b), however, the rate appears to be less than first order accurate. This is not entirely surprising. As pointed out in Section 5, the evaporation rate is related to the derivative of the scalar fields, and therefore is expected to be less accurate than the scalar field which is only first order.

7 Vaporizing Droplet in Uniform Flow

A final test is performed of a vaporizing droplet in uniform flow. There is no known analytical solution to this problem. A droplet is placed at the center of a gaseous domain with inflow velocity $U=40$ m/s. Unlike the simulation of the previous section, the assumption $\bm{u}_{L}=0$ is not imposed, so the droplet is able to move and deform freely in space. The same fluid properties are used as the $d^{2}$ law case, Table 4. The solution is initialized with a $d^{2}$ law initialization. The parameters chosen lead to a flow Reynolds number $Re=\rho_{G}Ud/\mu_{G}=400$ and gas Schmidt and Prandtl numbers as $Sc=\mu_{G}/\left(\rho_{G}D\right)=0.192\text{ and }Pr=\mu_{G}/\left(\rho_{G}\lambda_{G}\right)=0.192$, respectively. The surface tension is not assumed to be zero, and is given as $\sigma=0.0237$ N/m. This leads to a Weber number $We=\frac{\rho_{G}U^{2}d}{\sigma}\approx 7$.

The simulation used $128$ grid cells in each direction, which was chosen to resolve the gas film thickness, $\delta$. Following the paper by Sirignano and Abramzon [27], the following correlation for gas film thickness can be derived,

[TABLE]

For the parameters in this study, Eq. 95 yields $\delta\approx 0.157d$ or $\delta/\Delta\approx 2.51$. Equation 95 is noted to underpredict the film thickness [27], so the use of Eq. 95 can be considered a upper bound on the required resolution.

A series of visualizations is shown in Fig. 17. Figure 17(a) demonstrates the initial condition. As the flow progresses, aerodynamic forces cause the droplet to flatten and move slightly in the domain (Fig. 17(b)). Figure 17(c) shows the wake becoming asymmetrical and the formation of a large vortex under the droplet. After some time, the droplet wake becomes chaotic due to the interaction of several shed vortices, Fig. 17(d). The temporal variation in vaporization rate is shown in Fig. 18. The vaporization rate is greatly enhanced by the flow when compared to the $Re=0$ limit (the $d^{2}$ law). The vaporization rate appears to increase in time until the attached vortices begin to shed, indicating strong coupling between vortex formation and vaporization.

8 Conclusion

A numerical solver capable of simulating droplet vaporization in a complex flow field is demonstrated. A discretely conservative volume of fluid scheme is coupled with appropriate sources based on local thermodynamic equilibrium of the vapor mass fraction and temperature fields. The evaporation sources are coupled to the scalar transport equation using an unconditionally stable, monotone scheme. The convergence and stability properties of the solver are explored, and it shown that first order convergence is seen for gas temperature and vapor mass fraction, consistent with the schemes used. Also consistent with the scheme accuracy is the fact that the liquid velocity seems to converge poorly. Future work will require the use of higher order schemes for the scalar transport, to avoid this issue. Finally, it is shown that when used with a fully implicit time integrator the solver is capable of running stably at high time steps (for example $C_{\lambda}=200$) with non-oscillatory solutions. This is an improvement over the common Crank-Nicolson time integrator, and the improvement comes without any sacrifice of accuracy.

As true test of the robustness of the numerical solver, three dimensional simulations at high Reynolds number are performed. A vaporizing droplet in uniform flow at $Re=400$ is demonstrated. Because of the high speed flow, there is significant deformation of the interface, so no analytical solution is known for this problem. It is noted that vaporization is enhanced by the flow, which is consistent with the literature on non-deforming vaporizing droplets [27].

9 Acknowledgments

This material is based upon work supported by the National Science Foundation Graduate Research Fellowship under Grant No. DGE-1650441 and by the Alfred P. Sloan Foundation. This work used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant number ACI-1548562. The XSEDE resource used was Stampede2 at the Texas Advanced Computing Center (TACC) at The University of Texas at Austin. The authors also acknowledge Advanced Research Computing at Virginia Tech for providing computational resources and technical support that have contributed to the results reported within this paper.

10 Appendix

10.1 Monotone Schemes

Harten [28] in the development of total variation diminishing (TVD) schemes defines three types of schemes. Monotonicity preserving schemes, TVD schemes, and monotone schemes. Monotonicity preserving schemes are those that do not introduce new extrema to the solution, nor do they change the value of the current extrema. TVD schemes are defined in terms of a mathematical property called the total variation, and they are shown to be a class of monotonicity preserving schemes. The final category, monotone schemes, are a robust class of TVD schemes, and therefore also preserve monotonicity. Following Harten, define an update operator such that

[TABLE]

where the vector $\bm{\phi}=[\phi_{1},\phi_{2},\ldots]^{\top}$. A monotone scheme is one such that the coefficients of $\bm{\mathcal{L}}$ are all positive. Comparing Eq. 30 to Eq. 96, the chosen time-space integration scheme is monotone when the matrix,

[TABLE]

is invertible, and the coefficients of ${\bm{\mathcal{A}}}^{-1}$ are all positive. This property follows directly from the fact that $\bm{\mathcal{A}}$ is a nonsignular M matrix[29]. That is, $\bm{\mathcal{A}}$ is an M matrix, and it is strictly diagonally dominant with positive diagonal.

10.2 Farfield Boundary

The $d^{2}$ law is the vaporization of an isolated droplet in a quiescent field, which calls for an infinite domain. In this work, the domain is approximated using a triply periodic cube, with a fixed droplet in the center. The domain is made to be much larger than the droplet so that the effects of the boundary will be minimal. Still, as the finite sized domain fills with vapor, the effects of the boundary condition will be felt. A study was performed to look at the impact of domain size (measured in droplet diameters) on the evaporation rate. The droplet diameter, fluid properties, and initial conditions were made to be consistent with those in Section 6. For this study, the grid resolution is fixed so that the number of grid cells per droplet diameter is always $8$. The domain sizes were chosen to be integral multiples of the droplet diameter, viz., $4$, $8$, and $16$ diameters. The results of this study can be seen in Fig. 19. At the beginning of the simulation, all domains give the same trend in $d^{2}$. As the simulation progresses, boundary effects become important. This manifests itself in a gradual plateau of the droplet square diameter as a function of time. Simulations with a greater ratio of domain size to droplet diameter are able to continue longer without this plateau. As a reference, a straight line is shown by extrapolating the slope from the initial condition to the first time step. This reference was chosen instead of the analytical result, $1/\tau$, since the effective vaporization rate, and therefore the slope, changes under mesh refinement. For the given parameters, a domain eight times the droplet diameter leaves an error of only about $1.3\%$ in the slope at $t=0.01\tau$, which is the standard used in Section 6.

11 References