3.2.20. Volume of Fluid Model

Beta Capability

The volume-of-fluid capability is a beta feature. The numerics for laminar isothermal flow are well tested, but advanced features such as non-isothermal reacting and evaporating simulations are more experimental.

The volume-of-fluid (VOF) capability allows simulation of two-phase systems. The phases are tracked with a conserved scalar, $\alpha$ , which is the volume fraction of fluid in a given control volume. This scalar is $1$ in the liquid, $0$ in the gas, and between $0$ and $1$ in the transition region.

There are a number of numerical challenges associated with multi-phase modeling using a diffuse interface. Many of these challenges exist regardless of whether one uses a level set or a volume-of-fluid approach, or some hybrid of the two.

3.2.20.1. Governing Equation

The basic advection equation for the volume fraction ( $\alpha$ ) in VOF equations is

(3.595) $\frac{\partial \alpha}{\partial t} + \vec{u} \cdot \nabla \alpha = S_{\alpha}$

Expanding the conservative form of the convection term with the chain rule as

(3.596) $\nabla \cdot \left(\vec{u} \alpha\right) = \vec{u} \cdot \nabla \alpha + \alpha \left( \nabla \cdot \vec{u} \right)$

we can express the governing equation in corrected conservative form as

(3.597) $\frac{\partial \alpha}{\partial t} + \nabla \cdot \left(\vec{u} \alpha\right) = \alpha \left( \nabla \cdot \vec{u} \right) + S_{\alpha}$

where the right hand side term is $0$ for an incompressible flow with no phase change. The standard advection operator is typically too diffusive for practical use, so one of two approaches are typically used:

A geometric advection operator [84, 85, 86, 87]

An additional interface compression term [88, 89]

Several varieties of geometric advection operators have been proposed in the literature [84, 85, 86, 87], many of which are only applicable for 2D quadrilateral meshes. Those that can be extended to higher dimensions often do so at considerable complexity. For a more detailed discussion of the different options, see [90]. The form implemented here uses the interface compression approach to modify the governing equation to be

(3.598) $\frac{\partial \alpha}{\partial t} + \nabla \cdot \left(\vec{u} \alpha\right) + \nabla \cdot \left( \vec{u_c} \alpha (1 - \alpha) \right) = \alpha \left( \nabla \cdot \vec{u} \right) + S_{\alpha}$

where $\vec{u_c}$ is the compressive velocity. A common form for this is

(3.599) $\vec{u_c} = C_{\alpha} |\vec{u}| \vec{n}$

where $C_{\alpha}$ is a constant and $\vec{n}$ is the normal vector at the phase interface.

When phase change is present, the source term $S_{\alpha}$ is non-zero. The volume source of liquid due to evaporation is

(3.600) $S_{\alpha} = -\frac{\dot{m}_{evap}}{\rho_L}$

The evaporation term is calculated using the method described by Hardt and Wondra [91] which converts an evaporative mass flux to a volume source using the interphase area, $A$ , and the mass flux predicted by the evaporation model, $m^{''}_{evap}$ , as

(3.601) $\dot{m}_{evap} = m^{''}_{evap} A$

(3.602) $A = W |\nabla \alpha| \frac{\int{|\nabla \alpha| dV}}{\int{W|\nabla \alpha| dV}}$

(3.603) $W = \alpha$

One variation of this that we explored is $W = (1-\alpha)^2$ , which is well suited for a diffuse interface approach. We find that the form by Hardt and Wondra is sometimes destabilizing to the diffuse interface and leads to an artificially high surface regression rate. For a very simple 1D problem or a sharp interface method, the $W = (1-\alpha)^2$ weighting scheme places all the evaporation outside the liquid and does not correctly predict the surface regression rate. However, it can be better for a diffuse interface representation since evaporation acts as a stabilizing liquid sink to counteract interface diffusion. When using a sharp interface method or a 1D problem where there is no interface diffusion, the weighting of Hardt and Wondra ( $W=\alpha$ ) should be used.

Different evaporation models can be used to define $m^{''}_{evap}$ . The model used by Hardt and Wondra uses deviation from the saturation temperature to define the rate as

(3.604) $m^{''}_{evap} = K(T-T_{sat})$

(3.605) $K = \frac{2\beta}{2-\beta} \frac{h_{fg}}{\sqrt{2\pi R}} \frac{\rho_g}{T_{sat}^{1.5}}$

This form is only physically applicable to either fully saturated gasses, or interface temperatures above the boiling point since it does not take into account the gas composition. As such, this model only allows positive evaporation rates (so no condensation) in order to produce realistic behavior during heat-up and away from the hot regions of the simulation. When using a multi-species gas, the more complete Hertz-Knudsen-Langmuir model can be used, where the rate is driven by the difference between the partial pressure of the vapor and the saturation pressure. This model will be added to the available evaporation models once the multi-species gas treatment is complete.

3.2.20.2. Interface Reconstruction

When the phase interface is represented by a volume-of-fluid variable, numerical diffusion in the normal advection operator causes the interface (cells with volumes of fluid between $0$ and $1$ ) to become wider over time. Although the total volume of liquid is conserved, the volume inside the $\alpha=0.5$ isosurface may change as a result, and the effect of the higher density liquid will be felt well outside the correct liquid interface location. Techniques to counter-act this diffusion in volume-of-fluid include using compressive advection operators in the VOF equation, using a separate anti-diffusive sharpening step, using a geometric reconstruction of the interface to define the advection, and using sharpening and liquid relocation schemes.

With a level set approach, the interface location is defined by the $\phi=0$ contour. By definition, the interface does not become diffuse, but there are no guarantees of global conservation of liquid either. Techniques to address this typically involve doing a volume-conserving redistancing operation to the level set field. Such operations can preserve the total volume of liquid, but may not preserve the local shape of the liquid-gas interface.

In both VOF and level set approaches, calculation of the interface curvature is required to determine surface tension forces. This is done by taking the divergence of the interface normal vector ( $\vec{n}$ ).

(3.606) $\kappa = -\nabla \cdot \vec{n}$

It should be immediately apparent that $\kappa$ will contain not only discretization error from the divergence operator, but also discretization error from the calculation of the normal vector, which is itself a gradient of some other quantity (call it $\alpha_s$ here, it could be the level set field, the volume of fluid field, or a smoothed version of either).

(3.607) $\vec{n} = \frac{\nabla \alpha_s}{|\nabla \alpha_s| + \epsilon}$

The simple approach to use $\alpha_s=\alpha$ results in considerable noise in the interface normal vector calculation, particularly in regions of high curvature relative to mesh sizes. In the coupled level set VOF approach (CLSVOF), $\alpha_s$ is the companion level set distance function that is either advected along with the volume fraction or redistanced from the volume fraction contour at $0.5$ (so $\alpha_s = \phi$ ).

In Fuego we have implemented both a CLSVOF approach and a diffusive approach, where $\alpha_s$ is a diffusively smoothed version of the volume fraction field, $\alpha$ , using a user-specified Fourier number ( $Fo$ ) and number of iterations.

(3.608) $\alpha_s^0 = \alpha$

(3.609) $\alpha_s^{n+1} = \alpha_s^n + Fo \Delta x^2 \nabla^2 \alpha_s^n$

In addition, recall that the analytical curvature for a spherical drop is a function of radius: $\kappa=2/r$ . This means that if $\kappa$ is calculated in cells that span an interface of finite, non-zero width, there will also be a spatial variation of $\kappa$ across that interface. This means that although in reality the surface tension force is applied sharply at the interface, in a diffuse interface scheme it will be applied over a finite interface thickness of several cells.

This noise and spatial variation in the curvature translates directly into noise in the surface tension force and is the sole cause of the so-called “parasitic currents” that plague diffuse interface methods. Many authors have demonstrated that when using a uniform, analytically prescribed curvature the parasitic currents in the classical static drop test drop to machine precision levels if the forces are otherwise properly balanced.

There are various techniques for improving the curvature calculation. Calculation of the normal vector when the VOF field is sharp is very noisy, so VOF schemes often use a smoothed field or higher order least squares technique. Coupled level set VOF (CLSVOF) carries a level set field along to use for calculating this normal vector on a smooth field (the level set field).

Finally, once the interface has been advected and sharpened, and the curvature has been properly calculated, the surface tension force must be integrated into the momentum and continuity equations in such a manner as to not induce spurious velocities. This is known as the “balanced force” technique, and involves including a surface tension force vector consistently with the pressure gradient.

3.2.20.3. Continuity Equation

In order to handle the high density ratios in multiphase simulations, the form of the continuity equation is modified by dividing both sides of the equation by density. Interpolation of nodal quantities to subcontrol surfaces is an important consideration to note. Quantities interpolated to faces are indicated with an $f$ subscript in the following equations, while nodal quantities are given no subscript.

The regular form of the pressure Poisson equation in discrete form is

(3.610) $\sum_f \left( \tau_f \nabla p_f \cdot A_f \right) = \sum_f \left(\left(\rho \vec{u}\right)_f \cdot A_f + \left(\tau \nabla p\right)_f \cdot A_f \right)$

which includes the projected nodal gradient of pressure $\nabla p$ . In the prior form, the projected nodal gradient $\nabla p$ was calculated using the divergence theorem as

(3.611) $\nabla p = \frac{1}{V}\sum_f{p_f A_f}$

For the VOF implementation we use the method described by Francois et al. [92] to calculate a projected nodal density-scaled force term ( $F$ ) based on calculation of forces at faces, which includes both pressure gradient and surface tension forces. The balanced force term is

(3.612) $\vec{F} = \nabla p - \kappa \sigma \nabla \alpha$

and the projection from subcontrol faces to nodes is done using the method by Francois as

(3.613) $F_{j} = \frac{\sum_f \left(\nabla p_{f,j} - (\sigma \kappa)_f \nabla \alpha_{f,j} \right) / \rho_f}{\sum_f |A_{f,j}|}$

This term is updated after the solution to the pressure equation, when updated values for pressure and volume-of-fluid have both been calculated in order to keep them synchronized. This nodal force term is included as a source term in the momentum equation (negative, scaled by volume, and multiplied by nodal density) and can be added to the fourth order pressure projection scheme in discrete form as

(3.614) $\sum_f \left( \left(\frac{\tau}{\rho}\right)_f \nabla p_f \cdot A_f \right) = \sum_f \left( \vec{u}_f \cdot A_f + \left(\tau \vec{F}\right)_f \cdot A_f + \left(\frac{\tau}{\rho}\right)_f (\sigma \kappa)_f \nabla \alpha_f \cdot A_f \right)$

where $\vec{F}$ is the projected nodal density-weighted balanced force, $\kappa$ is the curvature, and $\sigma$ is the surface tension.

The advantage of this projection scheme over a simpler approach where one might define the nodal force term (as done by Lin et al. [93]) as

(3.615) $F = \frac{1}{V} \left( \sum_f p_f A_f - \sigma \kappa \sum_f \alpha_f A_f \right)$

is that in such a scheme the evaluation of curvature ( $\kappa$ ) occurs at a different location than the gradient of the volume-of-fluid term ( $\nabla \alpha$ ). In the presence of gradients in $\kappa$ –which are guaranteed to occur in real problems–this introduces a force imbalance that can cause significant spurious non-physical currents in the fluid. Spatial inaccuracies in curvature are the primary cause of spurious currents in VOF simulations – and a baseline test for any balanced force implementation is to see how a static drop with a prescribed (not calculated) value for curvature behaves. In a properly balanced force implementation the resulting parasitic currents should be basically zero (on the order of machine precision). However, a weakness of relying solely on a prescribed curvature test is that it fails to show the force imbalance that arises in the projected nodal force calculation due to interpolation of the curvature. The form published by Francois et al. [92] and implemented in Sierra/Fuego not only recovers machine-precision parasitic currents with the prescribed curvature test, but also keeps them acceptably low when using calculated curvatures.

In addition, when there is phase change there is an additional source term on the right hand side of this equation, which is

(3.616) $S = \dot{m}_{evap}\left(\frac{1}{\rho_g} - \frac{1}{\rho_L}\right)$

3.2.20.4. Property Evaluation

When using VOF, you are required to provide three material models. There is one model for the liquid, one for the gas, and one top level model. The top level model simply specifies the names for the materials to use for gas and liquid, as well as specifying properties like surface tension. Properties needed for flow calculations are combined from the liquid and gas properties using the following rules

A simple averaging is used for non-specific quantities,

(3.617) $\phi = \alpha \phi_L + (1 -\alpha) \phi_g$

where $\phi$ is thermal conductivity ( $k$ ), density ( $\rho$ ), viscosity ( $\mu$ ), absorption ( $a$ ), or density-pressure-derivative ( $\frac{\partial \rho}{\partial p}$ ).

For specific quantities, such as specific heat ( $C_p$ ) and enthalpy ( $h$ ), a density-weighted average is used

(3.618) $\phi = \frac{\alpha \rho_L \phi_L + (1 -\alpha) \rho_g \phi_g}{\rho}$

For properties that are only relevant for the gas phase, the property value is simply copied from the gas phase

(3.619) $\phi = \phi_g$

Examples include species enthalpy, mass diffusivity, and molecular weight.