## 1. Introduction

Simulations of free surface flows have progressed rapidly over the last decade, and it is now possible to simulate the motion of complicated waves and their interactions with structures considering even deformable bubbles in turbulent flows. In the continuum mechanics, there are two methods to express the motion in the environment. The first description is the Eulerian approach. In this method, attention is paid to a special volume in the space. A mesh remains fixed in the Eulerian method and fluid regions change in shape and location on the mesh. It uses a fixed grid system which is not transformed during the solution procedure. The fluid is studied while passing this volume and continuously replaced in time. Therefore, this method is not appropriate for formulation of basic equations of fluid movement. The Eulerian method has some limitations. For example, when the portion of the perimeter to the area of a zone of fluid is large, the error of this method is increased. In the Eulerian method, it is not possible to decompose the equation on the boundaries with the same precision of inner region of fluid and accordingly, the finer mesh should be used near the boundaries. Therefore, when the free surface of a discontinuous region is modeled by this method, finer grid should be employed in order to achieve more precise results, specifically if this surface has large deflections. This is crucial when the portion of the area to the perimeter of a zone is low, for example on phase of a multiphase fluid. In this case, using finer mesh could increase the portion of the number of the inner elements to the boundary elements, which in turn, increases the precision of the numerical solution. The main superiority of the Eulerian description is the possibility of modeling of complicated surfaces. For example, the collapse of a column of a fluid could be modeled in the Eulerian grid which is shown in **Figure 1**.

In Lagrangian method, the flow field of the considered fluid is covered by a mesh moving with the fluid. The fluid boundaries always coincide with the grid boundaries and the fluid inside each cell of the grid always remains in that computational cell. Although this method is not applicable to flows undergoing large distortions, where meshes can be twisted into unacceptable shapes, but its advantage is the ease with which it handles free surfaces and interfaces, which makes it applicable to a wide variety of problems. For example, the grid shown in **Figure 2** is Lagrangian in the vertical coordinate. For free surface problems, if the free surface movement or the tangential acceleration gradient in the perpendicular direction to its surface is not large, the Lagrangian method can be used to simulate free surfaces. The grid lines are located on the free surface and move with it. Therefore, there is no need for any special boundary condition in this location [1].

## 2. Governing equations

Governing equations for a compressible viscous fluid flow with no phase change are as follows:

In these equations, **u**, **F**_{s} are density, velocity vector, time, total pressure, kinematic viscosity, gravity acceleration, and body forces, respectively. Body forces include forces due to surface tension in the interface. Here, properties of a fluid such as density and viscosity are included in the equations. However, it should be kept in mind that the information changes from one fluid to another. Thus for mesh‐based numerical methods, new properties based on fluid properties of both materials should be considered for Eq. (2) in the cell containing the free surface, and the governing equations should be rewritten in the following form:

where indices 1 and 2 show first and second fluids properties, and

(7) |

This phase indicator function is the fluid property or volume fraction, which moves with it and can be derived as follows:

This function can be used to calculate the fluid properties in each phase as a weight function. In order to use a set of governing equations using the weight function, each fluid property should be calculated based on the volume occupied by this fluid in the surface cell as expressed in Eqs. (9) and (10) [2]:

Free surfaces considered here are those on which discontinuities exist in one or more variables. This has been the challenge for researchers to omit or reduce this problem as much as possible. The transient state as well as phenomena such as surface tension, changing of fluid phase and Kelvin‐Helmholtz instability makes numerical simulation of such problems cumbersome. It is expected that methods used to simulate interface of fluids have a number of characteristics. These include mass conservation, simulating the interface as thin as possible, being able to reproduce complicated topologies, generalization of expansion to 3D problems, and being able to model surface phenomena and be computationally efficient.

## 3. Free surface modeling methods

There are different methods to simulate free surface flow, each of which has its own advantages and disadvantages:

### 3.1. Donor‐acceptor method

The main idea of donor‐acceptor approach is that the value of volume fraction in downwind cell, the acceptor cell, is used for anticipation of transferring fluid in each time step. The problem in this approach is that using downwind cell in calculations may lead to unreal situations which are values out of zero and unity domain in surface cells. **Figure 3a** shows this method with the first fluid with gray color and volume of fluid equals to unity. It could be seen that using donor‐acceptor approach with downwind differencing scheme results in values greater than unity in donor cell. It is because the second fluid in the acceptor cell is greater than the value needed in the donor cell. Similarly in **Figure 3b**, using downwind differencing scheme leads to negative values for volume of fluid, which is because the needed fluid in acceptor cell is more than what is in the donor cell [3].

In order to be assured that volume of fluid is between zero and one, the amount of fluid or volume of fluid in donor cell should be used to regulate the estimated fluid transferring between two adjacent cells [5].

One drawback of donor‐acceptor method is that this method changes any finite gradient into step, and consequently increases the slope of the surface model in the direction of flow. This problem was alleviated by proposing a method to consider the slope of interface for flux transferring in adjacent cells by Hirt and Nichols [6]. For this purpose, a donor‐acceptor equation was proposed so that it could detect the direction of the flow in interface and then define the upwind and downwind cells accordingly. Thereafter, this model was expanded for 3D domains by Torrey et al. [7]. The Surfer method is one version of volume of fluid which deals with merging and fragmenting of interfaces in multiphase flows [8].

The volume of fluid method is one of the most popular methods for anticipation of interfaces, and many researches have been conducted based on this method including dam break, Rayleigh‐Taylor instability, wave generation and bubble movement [6, 9–12]. This method was modified in 2008 to get more accurate results by considering diagonal changes in fluxes of adjacent cells for structured grid domains [13, 14].

### 3.2. The Hirt‐Nichols method

The volume of fluid (VOF) method was first proposed by Hirt and Nichols [6]. In this method, similar to the SLIC method, free surfaces can be reconstructed based on parallel lines with respect to one of the principal coordinates of the system. However, nine neighboring cells are considered for flux changes and defining the normal vector in a desired cell. Then, free surface is considered as either a horizontal or a vertical line in cell with respect to the relative normal vector components. **Figure 4** shows the actual free surface and what was simulated by Hirt‐Nichols method.

Upwind fluxes are used for fluxes parallel to the reconstructed interface, while donor‐acceptor fluxes are used for those fluxes normal to it. For instance according to **Figure 5a**, the interface in the cell *x* coordinate in the face **Figure 5b**), and this cell is considered a downwind cell for the cell

where

The “min” operator has been designed to ensure the fluid leaving the cell

### 3.3. Flux Corrected Transport (FCT) method

The FCT method is based on the idea to present a formulation which combines the upwind and downwind fluxes. This formulation aimed to leave out upwind numerical diffusion and instability of downwind scheme [15]. Idea of neighboring fluxes based on higher order translate scheme was first proposed by Boris and Book [16] and then developed by Zalesak [17] to multidimensional.

In this method calculations consist of some steps. First, an intermediate value of volume of fluid, **Figure 6** shows schematically the solution for a 1D governing equation of fluid volume fraction for cell *i* as:

where

Thereafter, an anti‐diffusive flux is needed to be defined (*F ^{L}*) in order to correct the diffusion of the previous step. This is the difference between upwind and downwind fluxes as:

To make this stable, a correction factor,

### 3.4. Youngs’ method

This method was first proposed by Youngs in 1982 [18]. It was then developed by Rudman [19] with more details. In this method, at first the slope of the interface position is estimated. Then, the free surface is defined as a straight line with the slope of

Assuming that

Using components of the normal unit vector, the angle

The angle

It is possible to set **Figure 7**.

What is behind this conclusion is as follows:

Four side fractions **Table 1**. In this table, positive value is set for velocities towards the outer edges of a cell, and there is no flux calculation for negative velocities into the cell.

### 3.5. Piecewise Linear Interface Calculation (PLIC) method

To solve fluid volume transfer equation with FDM or FVM, diffusion error in interface reconstruction occurs. This leads to poor modeling of free surfaces, specifically in the interface of two adjacent fluids with large density difference. PLIC is one of the methods to reconstruct the interface between fluids with second‐order accuracy [20]. It can increase the accuracy of transferred flux estimation and geometric fluid distribution in each cell. In this method, unit normal vectors of the surface are calculated based on the volume fraction of fluid using Youngs’ least square method as:

where

(22) |

where **Figure 8** [21]. All other situations can be achieved with a mirror reflection of the first quarter with respect to the x and y axes and bisectors between them. The exact position of the free surface is determined defining surface unit normal vector using volume fraction of fluid in each cell. To do this, extreme values of

in which

**Figure 9**.

When unit normal vector of a surface is defined, the true position of the interface can be easily determined using volume of fluid in each cell.

### 3.6. Higher order differencing schemes

Another method to reconstruct the interface between two fluids is to discretize the convection term using higher order differencing schemes or blended differencing scheme. The accuracy of less/non‐diffusive schemes and compressive schemes was compared by Davis [22]. Less/non‐diffusive schemes prevent the interface profile from being diffused. Compressive schemes not only prevent the interface from being diffused, but also omit any diffusion in the neighboring of the interface. Thus, they are considered as powerful tools for thin interface simulation.

Ghobadian [23] applied the higher order scheme proposed by Van Leer [24]. However, his results showed that this scheme has poor ability in terms of removing diffusion. Therefore, he proposed solutions for decreasing numerical diffusion. Other methods for omitting diffusion proposed by Pericleous and Chen [25] proved to be associated with interface diffusion. Although first‐order upwind or downwind schemes lead to diffusion, higher order methods result in numerical fluctuations in the interface. There are other methods for reducing the interface as follows:

#### 3.6.1. Compressive Interface Capturing Scheme for Arbitrary Meshes (CICSAM) scheme

The CICSAM scheme, presented by Ubbink, is a combined method to reduce the diffusion problem in interface modeling. This method imposes some limitations on the fluid fraction value. It is obvious that the value of a fluid in a cell should be constant in the absence of a source. The CICSAM approach presents an equation for free surface volume fraction as:

where

where

More details on determining

where **Figure 10**, and defined as:

The CICSAM method well satisfies the bounds defined within it, and can be accurately reconstruct the free surface. The basis of the method, however, is on the 1D equations and linearization, which makes it less accurate for 3D modeling reconstructions.

#### 3.6.2 THOR scheme

This scheme is based on the CICSAM and switches smoothly between the upper bound of the universal limiter and ULTIMATE‐QUICK, a combination of the universal limiter and QUICK, considering the angle between the interface and the direction of motion [27].

Analogous to CICSAM, this scheme is an algebraic advection scheme for the interface, which is designed for the implicit time advancing algorithm. In this method _{f} as follows:

where

#### 3.6.3. Higher Resolution Artificial Compressive (HiRAC) scheme

HiRAC scheme is another modification of the CICSAM method [28]. This newly proposed method tries to improve the computational efficiency and maintain the accuracy. In this method, the weighing factor,

where *m*=2, the new formulation reduces to the weighting function of Ubbink and Issa [26]. As *m* increases, the interpolation becomes more biased towards the diffusive higher resolution scheme. It is shown that *m*=2 provides a good balance between the compressive and diffusive higher resolution schemes.

#### 3.6.4. High Resolution Interface Capturing (HRIC) scheme

This method is somehow similar to CICSAM, which benefits from a combined interpolation scheme. In HRIC method, the difference between two upwind schemes is calculated based on the normal vector angle of the free surface as [29]:

The portion of each of the two terms in the above equations can be defined as:

In this way,

It should be noted that an improved scheme of HRIC, called Flux-Blending Interface-Capturing Scheme (FBICS), has been recently proposed. In this method, analogous to CICSAM and HRIC, the difference between two upwind schemes is calculated based on the normal vector angle of the free surface. Based on FBICS, Eq. (33) can be reformulated to obtain a more accurate scheme as:

Some other modifications are also proposed by Tsui et al. [30].

#### 3.6.5. Switching Technique for Advection and Capturing of Surfaces (STACS) method

One of the drawbacks of HRIC and CICSAM schemes is high Courant numbers. Both methods lack a proper switching strategy to accurately model the interface when Courant number increases. The Courant number, *Cn*, in HRIC method can be written as follows:

which is suitable for Courant numbers between 0.3 and 0.7. For a *Cn* below 0.3 the scheme is not modified, while for a Courant number above 0.7 the upwind scheme is used. This is true for CICSAM when the *Cn* is equal to unity.

STACS method has been proposed to improve the accuracy and stability of the results specifically in high Courant numbers by Darwish and Moukalled [31]. It uses an implicit transient discretization, i.e. no transient bounding is applied, and in order to minimize the stepping behavior of HRIC scheme, a modification is proposed. In this method, applying

where

(40) |

This enables a rapid but smooth switching strategy that works very well, especially where the normal to the free‐surface face is not along the grid direction.

#### 3.6.6. Inter‐gamma scheme

In this method, presented by Jasak and Weller [32], free surface compression is modeled using additional compressive artificial terms.

where

where

where

**Figure 11** shows the NVD for inter‐gamma scheme.

### 3.7. Integrated methods

As mentioned before, volume of fluid is among the most popular methods in free surface modeling. Having in mind that this method is based on defining a discontinuous function, the color function, there is not a unique form for free surface. Therefore, it is required to reconstruct the free surface using volume of fraction function. In one hand, VOF method satisfies the conservation of mass while it is unable to calculate free surface parameters including curvature radius and normal unit vector directly. On the other hand, in level set methods as the distance function is smooth, the surface geometry can be easily calculated, while satisfying the conservation of mass is very demanding. In order to resolve the problems of level set methods, a number of different researches have been conducted. For example, higher order schemes were proposed to improve the conservation of continuity equation by Peng et al. [33]. Adaptive mesh refinement techniques were also proposed to increase the accuracy of the local mesh consistency. In 2009, an integrated method known as hybrid Particle Level Set (PLS) was proposed to improve the accuracy of the results. However, the problem still remained in relation with mass conservation.

In order to take the advantages of both methods and eliminate their disadvantages, integration of volume of fluid and level set methods was proposed in a new scheme known as coupled level set and volume-of-fluid (CLSVOF) method to model two‐phase incompressible flows by Sussman and Pucket [34]. It should be noted that although accurate, this method cannot be easily employed, because these two methods, VOF and level set, should be individually solved and their effects need to be coupled based on the reconstructed interface.

## 4. Calculating surface tension

Defining the pressure difference inserted on the surface of two fluids with different densities and tension stresses is one of the most demanding problems in fluid mechanics. One method to do this is the Pressure Calculation based on the Interface Location (PCIL) method which is presented here. Surface tension, that changes the value of variables in momentum equations, imposes a discontinuity at the position of the interface between two fluids [21].

Stress from surface tension inserts a force upon the interface. The resultant force is perpendicular to the surface and its curvature is dependent on the geometry of the surface. Surface tension can be considered in two ways. In the first approach, it is considered as a boundary condition in the equations for the surface. This needs using an iterative method for true approximation of pressure, which in result, increases the time and cost of calculation and consequently makes it inefficient. In order to address this problem, some other methods have been proposed in which the precise calculation of interface position is not necessary. In these methods, the direct force of surface tension has been replaced with the body force in the momentum equation. The Continuum Surface Force (CSF) method is a base method for calculation of body forces of fluid surface tension [2]. The body forces can be considered to act smoothly on a narrow strip of cells in interface zone. In this method the surface stresses are replaced with the body forces which are calculated as:

where

Another approach based on CSF method was proposed by Torrey et al. [7] called Continuum Surface Stress (CSS), in which body forces of CSF method were replaced by tension tensors of surface tension based on the following equation:

where

It should be mentioned that employing CSF or CSS methods has some drawbacks. For instance, spurious velocities of the thinner fluid near the interface is one the reported problems.

A number of researches have been performed in order to resolve the problem of spurious velocities [35, 36]. In one approach, using virtual particles moving along with the surface could improve the results [37]. One of the latest methods presented in this field is PCIL. This method shows that having more precise border cells and calculating their associated pressure based on the momentum equation can lead to significant reduction of bothersome flows near this region. PCIL is a simple and efficient method of calculating free surfaces. The total pressure on the left side of the cell can be calculated as (see **Figure 12**):

where

where

On the other hand, the change in pressure

Accordingly, the above equation can be reformulated for pressure in the *K*th face of every common cell as follows:

where the second term introduces the normal force of the surface tension per unit area of the interface. This can be presented in the vector form as:

where **F**_{s} is the surface tension force vector.

One of the most fundamental steps to perform surface tension calculations is defining the curvature of the interface. Defining this curvature is not so demanding as long as the precise position of the interface is known. However, using volumetric tracing methods and equivalent alternatives representing the interface position make the estimation of the curvature cumbersome.

The method of volume of fluid presented by Hirt and Nichols [6] is one of the earliest methods in this field. In this method, a curve

The value of the surface curvature can be defined based on the distance function as:

To discretize the above equation, it is required to first calculate the normal vector of the surface based on **Figure 13** and the following relations, and consequently estimate the curvature:

One of the advantages of this approach, the level‐set method, is using a distance function which is smooth and uniform, so that it increases the simplicity of the calculation and accuracy of the results.

The value of the body force from surface tension of the cell faces can be calculated in CSF method as:

where

in which the ratio of the densities is inserted in order to reduce spurious velocities of the thinner fluid. The discretized version of Eq. (58) can be obtained as follows:

Based on what was discussed for PCIL method, the following relation is adopted to the present method:

It can be seen that in this equation, the ratio of densities is replaced by the variable *H*.

## 5. Parametric method for calculation of curvature of free surfaces

This newly proposed method is based on two sub‐models, the Four-Point Method (FPM) and the Three-Line Method (TLM). In the former sub‐model, a curve is fitted to the intersection of the points of grid lines for central and two neighboring cells, while the latter fits a curve to the free surface so that the distance between the curve and its linear interface approximation is minimized [42].

### 5.1. The Four-Point Method (FPM)

In the four‐point method, free surface (as illustrated in **Figure 14**) is approximated using a continuous function

where

In this method, the desired function is approximated using an n‐degree polynomial function with unknown constant coefficients. Therefore, we have:

where

It is supposed that Q is the set of

**Theorem 1**: if

**Proof**. It is obvious that

As we have

Now, it is claimed that there is an

Thus, our aim is to solve Eqs. (63) and (64), and one can write them in the following forms:

(71) |

or

such that

Now, the intervals

(73) |

This is a nonlinear set of equations and can be easily solved using Matlab or Lingo software.

### 5.2. The Three-Line Method (TLM)

In this method as illustrated in **Figure 15**, the main goal is to find a function as

Similar to what was discussed in the four‐point method, in this method the function

**Theorem 2**: Sequence of the solution of Eq. (75) converges to the solution of Eq. (74).

**Proof**: The method of proof of this theorem is similar to the previous theorem. In the same approach of FPM, the following problem is achieved:

(76) |

The steps of using the above equations are as follows:

Step 1: Read

Step 2: Solve Eq. (73) or (76) in the FPM or the TLM, respectively.

Step 3: If the previous step is infeasible, set *n*=*n*+1, and go to step 2, else set the value of target function in

Step 4: Set *n*=*n*+1, and solve Eq. (73) or (76) in the FPM or the TLM, respectively.

Then set the value of target function in

Step 5: If

## 6. Conclusions

In this chapter volume of fluid (VOF) scheme was introduced. This is one of the most effective methods employed in the simulation of two fluid flows interfaces with dramatic changes in density and viscosity. . These interfaces are represented implicitly by the values of a color function which is the fluid volume fraction. The advantage of the method is its ability to deal with arbitrarily shaped interfaces and to cope with large deformations, as well as interface rupture and coalescence in a natural way. In VOF the mass is rigorously conserved, provided the discretization is conservative. However, advecting the interface without diffusing, dispersing, or wrinkling is a big issue. This can either be performed algebraically, in schemes such as CICSAM or geometrically, in schemes such as PLIC. Herein, the viscous fluid governing equations which are Navier‐Stokes coupled with VOF equation were presented. Then the most popular VOF schemes such as donor‐acceptor, Hirt‐Nichols, FCT, Youngs, and PLIC were explained. CICSAM, HiRAC, HRIC, STACS, and some other up‐to‐date proposed methods were introduced and the accuracy and time calculation of each method were evaluated. Moreover, surface tension modeling and parametric study of interfaces were discussed. The author hopes this brief presentation of the VOF method will be beneficial for scientists and students in their further researches and will help them to massively and continuously expand this very challenging field of fluid mechanics.