Volume of fluid method

History

The volume of fluid method is based on earlier Marker-and-cell (MAC) methods. First accounts of what is now known as VOF have been given by Noh & Woodward in 1976, where fraction function C (see below) appeared, although the first publication in a Journal was by Hirt and Nichols in 1981. Since VOF method surpassed MAC by lowering computer storage requirements, it quickly became popular. Early applications include Torrey et al. from Los Alamos, who created VOF codes for NASA (1985,1987). First implementations of VOF suffered from imperfect interface description, which was later remedied by introducing a Piecewise-Linear Interface Calculation (PLIC) scheme. Using VOF with PLIC is a contemporary standard, used in number of computer codes, such as FLOW-3D, Gerris (software), ANSYS Fluent and STAR-CCM.

Overview

The method is based on the idea of a so-called fraction function C . It is a scalar function, defined as the integral of a fluid's characteristic function in the control volume, namely the volume of a computational grid cell. The volume fraction of each fluid is tracked through every cell in the computational grid, while all fluids share a single set of momentum equations. When a cell is empty with no traced fluid inside, the value of C is zero; when the cell is full, C = 1 ; and when there is a fluid interface in the cell, 0 < C < 1 . C is a discontinuous function, its value jumps from 0 to 1 when the argument moves into interior of traced phase. The normal direction of the fluid interface is found where the value of C changes most rapidly. With this method, the free-surface is not defined sharply, instead it is distributed over the height of a cell. Thus, in order to attain accurate results, local grid refinements have to be done. The refinement criterion is simple, cells with 0 < C < 1 have to be refined. A method for this, known as the marker and micro-cell method, has been developed by Raad and his colleagues in 1997.

The evolution of the m -th fluid in a system on n fluids is governed by the transport equation (actually the same equation that has to be fulfilled by the level set method distance function ϕ ):

with the following constraint

i.e., the volume of the fluids is constant. For each cell, properties such as density ρ are calculated by a volume fraction average of all fluids in the cell

These properties are then used to solve a single momentum equation through the domain, and the attained velocity field is shared among the fluids.

The VOF method is computationally friendly, as it introduces only one additional equation and thus requires minimal storage. The method is also characterized by its capability of dealing with highly non-linear problems in which the free-surface experiences sharp topological changes. By using the VOF method, one also evades the use of complicated mesh deformation algorithms used by surface-tracking methods. The major difficulty associated with the method is the smearing of the free-surface. This problem originates from excessive diffusion of the transport equation.

Discretization

To avoid smearing of the free-surface, the transport equation has to be solved without excessive diffusion. Thus, the success of a VOF method depends heavily on the scheme used for the advection of the C field. Whereas a first order upwind scheme smears the interface, a downwind scheme of the same order will cause a false distribution problem which will cause erratic behavior in case of the flow is not oriented along a grid line. As these lower order schemes are inaccurate, and higher order schemes are unstable and induce oscillations, it has been necessary to develop schemes which keep the free-surface sharp while also producing monotonic profiles for C. Over the years, a multitude of different methods for treating the advection have been developed. In the original VOF-article by Hirt, a donor-acceptor scheme was employed. This scheme formed a basis for the compressive differencing schemes.

The different methods for treating VOF can be roughly divided into three categories, namely the donor-acceptor formulation, higher order differencing schemes and line techniques. The donor-acceptor scheme is based on two fundamental criteria, namely the boundedness criterion and the availability criterion. The first one states that the value of C has to be bounded between zero and one. The latter criterion ensures that the amount of fluid convected over a face during a time step is less than or equal to the amount available in the donor cell, i.e., the cell from which the fluid is flowing to the acceptor cell. In his original work, Hirt treated this with a blended scheme consisting of controlled downwinding and upwind differencing.

In the higher order differencing schemes, as the name suggests, the convective transport equation is discretized with higher order or blended differencing schemes. Such methods include the Compressive Interface Capturing Scheme for Arbitrary Meshes (CICSAM) and High Resolution Interface Capturing (HRIC) scheme, which are both based on the Normalized Variable Diagram (NVD) by Leonard.

Line techniques circumvent the problems associated with the discretization of the transport equation by not tracking the interface in a cell explicitly. Instead, the fluid distribution in a cell an interface is obtained by using the volume fraction distribution of neighbouring cells. The Simple Line Interface Calculation (SLIC) by Noh and Woodward from 1976 uses a simple geometry to reconstruct the interface. In each cell the interface is approximated as a line parallel to one of the coordinate axes and assumes different fluid configurations for the horizontal and vertical movements respectively. A widely used technique today is the Piecewise Linear Interface Calculation by Youngs. PLIC is based on the idea that the interface can be represented as a line in R² or a plane in R³; in the latter case we may describe the interface by:

where n is a vector normal to the interface. Components of the normal are found i.e. by using the Finite Difference method or its combination with least squares optimization. The free term α is then found (analytically or by approximation) by enforcing mass conservation within computational cell. Once the description of the interface is established, the advection equation of C is solved using geometrical techniques such as finding the flux of C between grid cells, or advecting the endpoints of interface using discrete values of fluid velocity.