Volume of Fluid (VOF) History

The Volume of Fluid (VOF) method is based on the earlier Marker-and-Cell (MAC) method [1]. The MAC method used marker particles to locate where fluid exists in a fixed Eulerian grid. MAC was the first computational method to simulate the dynamics of an incompressible fluid with a free surface. The use of marker particles to track fluid is computationally expensive, particularly in three dimensions, and introduces computational noise because the change in grid element properties (such as mass) undergo discrete changes when a particle moves from one grid element to another. A variety of attempts were made to replace marker particles with interface tracking schemes, but these mostly failed since fluid masses often breakup or coalesce, which leads to the creation and destruction of interface surfaces.

The idea of a volume of fluid method that tracks fluid volume instead of fluid surfaces arose from studies of two-phase (water and steam) problems in which it is customary to use a volume of steam variable. The steam volume fraction is a continuous variable, recording the amount of steam in a mixture of water and steam. Extending this volume concept to a discontinuous variable to locate free surfaces in an incompressible fluid (e.g., a unit value in liquid and zero elsewhere) was first demonstrated in a 1975 publication, “Methods for Calculating Multi –Dimensional, Transient Free Surface Flows Past Bodies” [2], by Nichols and Hirt.

Making Computational Sense

The VOF concept makes computational sense because flow models generally use only one numerical value in each grid element to store dependent variables such as pressure, density, temperature, etc. So why should more than one variable be needed to define the fluid distribution within an element? If, for example, the fluid in an element is distributed in more than one blob, then more dependent variables would be needed for each of the blobs. From this point of view, recording only the fluid volume in an element makes sense. The problem, however, is the supposed discontinuous nature of the volume fraction variable. Tracking the movement of a discontinuous fluid interface through a Eulerian grid requires more information.

The problem of discontinuous fluid

This problem has been addressed by many people in a large body of publications. Nearly all of the proposed methods rely on some type of approximation based on an examination of the volume fractions in adjacent grid elements. For example, in one-dimensional flow, an exact method is easy to derive. Suppose liquid is flowing along a one-dimensional duct with a sharp interface separating the liquid from a gas. In grid elements upstream of the interface, the volume fractions are equal to 1, while downstream of the interface the volume fractions are equal to 0. In the grid element containing the interface, which has a volume fraction value somewhere between 0 and 1, it is easy to locate the interface within that cell because liquid must be located at the side of the cell connected to its neighbor containing a volume fraction of 1. The interface is then located a distance downstream from the cell edge connected to its liquid neighbor by the product of the volume fraction multiplied by the size of the cell. This location can then be used when advecting the fluid in such a way that the interface remains a sharp discontinuity. Unfortunately, in two or three dimensions such a simple method of locating the interface within a grid element does not exist.

One method proposed for advecting discontinuous fluid interfaces was presented in the 1980 Los Alamos Scientific Laboratory report, “SOLA-VOF: A Solution Algorithm for Transient Fluid Flow with Multiple Free Boundaries,” [3] by Nichols, Hirt and Hotchkiss, and in a 1981 publication, “Volume of Fluid (VOF) Method for the Dynamics of Free Boundaries,” [4] by Hirt and Nichols. Early applications of this program, mostly for light-water-reactor safety studies, can be found in [5] and [6].

VOF Variations

Many variations of the VOF method have been reported in the literature, but most of them do not follow the method used in the original article [4]. In particular, the original VOF method solved the Navier-Stokes equations for fluid dynamics only in the incompressible liquid and not in the surrounding gas. Instead, the fluid free surface was treated by boundary conditions and the list of grid elements that contained fluid was continually updated. The gas regions in the original model were assumed to have low densities whose momentum could be ignored, and to have gas pressures that were spatially uniform. The alternative used by most other VOF models is to use two-fluid simulations in order to avoid setting boundary conditions on interfaces. This option requires considerably more computational resources than the original method due to having to solve for the gas dynamics. Most two-fluid models also ignore the possibility of a velocity “slip” existing between gas and liquid at an interface. By ignoring the existence of a slip and moving the interface at some average velocity for the gas/liquid mixture it is possible that serious errors may occur.

Modeling Fluid Advection

Another point not always appreciated by developers of alternative VOF methods is the equation that is modeled for the advection of the VOF fluid fraction quantity F. The original method [4] used a conservative transport equation for F,

$latex \frac{{\partial F}}{{\partial t}}+\nabla \bullet \left( {F\vec{u}} \right)=0$

Many alternative VOF formulations, such as those using the level-set method for advection, use the non-conservative transport equation,

$latex \frac{{\partial F}}{{\partial t}}+\vec{u}\bullet \nabla F=0$

The advantage of the conservative method is that it offers one simple check on the accuracy of incompressibility in a simulation, since the fluid volume, which should not change, is easily computed and displayed. 

The TruVOF Solution

Among popular commercial codes that are available, only FLOW-3D is based on the original, one-fluid model referenced in [4]. Of course, many improvements to this software have been made over its lifetime, including many models for a variety of physical processes such as heat transfer, surface tension, phase changes, moving obstacles and fluid-structure interactions.

References

1. F.H. Harlow and J. E. Welch, “Numerical Calculation of Time-Dependent Viscous Incompressible Flow,” Phys. Fluids 8, 2182 (1965); J. E. Welch, F. H. Harlow, J. P. Shannon, and B. J. Daly, “THE MAC METHOD: A Computing Technique for Solving Viscous, Incompressible, Transient Fluid-Flow Problems Involving Free Surfaces,” Los Alamos Scientific Laboratory report LA-3425 (March 1966).
2. B.D. Nichols and C. W. Hirt, “Methods for Calculating Multi-Dimensional, Transient Free Surface Flows Past Bodies,” Proc. Of the First International Conference on Numerical Ship Hydrodynamics, Gaithersburg, Maryland, October 20-23, 1975.
3. B.D. Nichols, C. W. Hirt and R. S. Hotchkiss, “SOLA-VOF: A Solution Algorithm for Transient Fluid Flow with Multiple Free Boundaries,” Los Alamos Scientific Laboratory report LA-8355 (August 1980).
4. C.W. Hirt and B. D. Nichols, “Volume of Fluid (VOF) Method for the Dynamics of Free Boundaries,” Jour. Computational Physics, 39, 201 (1981).
5. B.D. Nichols and C. W. Hirt, Numerical Simulation of BWR Vent Clearing Hydrodynamics,” Proc. 1978 Annual Meeting ANS, San Diego, CA; Nuc. Sci. Eng. 73, 196 (1980).
6. C.W. Hirt and B. D. Nichols, “A Computational Method for Free Surface Hydrodynamics,” ASME 1980 Pressure Vessels and Piping Conf. San Francisco, CA (August 1980) Jour, Pressure Vessel Technology, 103, 136 (1981)