**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.

**Figure 4.** Free surface (a) actual surface and (b) reconstructed surface based on Hirt‐Nichols method [6].

**Figure 5.** Hirt‐Nichols scheme (a) actual surface and (b) reconstructed surface.

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 inter‐ face in the cell (*i*, *j*) is considered to have positive celerity direction with respect to *x* coordi‐ nate in the face *i* + <sup>1</sup> 2 in donor‐acceptor method. Therefore, the reconstructed surface in the (*i*, *j*) cell is vertical (**Figure 5b**), and this cell is considered a downwind cell for the cell (*i* + 1, *j*), (*αi*, *<sup>j</sup>* >*αi*+1, *<sup>j</sup>* ). According to donor‐acceptor method, transferred flux from the face (*i* + <sup>1</sup> 2 , *j*) can be calculated as follows:

$$F\_{i \ast \underline{\mathcal{K}}, j} = \delta \circ \left( \min \left[ \alpha\_{i,j} \delta \mathbf{x}, U\_{i \ast \underline{\mathcal{K}}, j} \alpha\_{i+1,j} \delta t + \max \left( 0, U\_{i \ast \underline{\mathcal{K}}, j} \left( 1 - \alpha\_{i+1,j} \right) \delta t - \left( 1 - \alpha\_{i,j} \right) \delta \mathbf{x} \right) \right] \right) \tag{11}$$

where *αi*, *<sup>j</sup> δx* is the maximum fluid available for exiting the cell (*i*, *j*); *U i*+ 1 2 , *j αi*+1, *<sup>j</sup> δt* is the esti‐ mation of downwind flux for the volume of the fluid, *α*; *U i*+ 1 2 , *j* (1−*αi*+1, *<sup>j</sup>* )*δt* is the estimation of downwind flux resulted from the convey of void portion of the cell (*i*, *j*); and (1−*αi*, *<sup>j</sup>* )*δx* is the maximum void portion which can exit from the cell (*i*, *j*).

The "min" operator has been designed to ensure the fluid leaving the cell (*i*, *j*) is not more than the calculated available fluid in it from the previous time step. As the fluid in a cell transfers, so does the whole void space in the cell. Thus, the "max" operator has been designed in order to assure that amount of void exits the cell is bounded by what was in it calculated from the previous step. In this scheme, the combination of downwind and upwind fluxes has been considered in such a way that not only the solution stability is guaranteed, but also avoids the numerical diffusion.

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

for 3D domains by Torrey et al. [7]. The Surfer method is one version of volume of fluid which

The volume of fluid method is one of the most popular methods for anticipation of interfa‐ ces, 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

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‐

**Figure 4.** Free surface (a) actual surface and (b) reconstructed surface based on Hirt‐Nichols method [6].

**Figure 5.** Hirt‐Nichols scheme (a) actual surface and (b) reconstructed surface.

deals with merging and fragmenting of interfaces in multiphase flows [8].

adjacent cells for structured grid domains [13, 14].

370 Numerical Simulation - From Brain Imaging to Turbulent Flows

**3.2. The Hirt‐Nichols method**

Nichols 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, *α* \* , must be defined based on a first‐order scheme. **Figure 6** shows schematically the solution for a 1D governing equation of fluid volume fraction for cell *i* as:

$$\alpha\_i^\* = \alpha\_i^\* - \frac{1}{\delta \mathbf{x}} \left( F\_{i \ast \underline{\mathbf{x}}}^L - F\_{i \cdot \underline{\mathbf{x}}}^L \right) \tag{12}$$

**Figure 6.** Three adjacent typical cells in FCT method.

where *<sup>F</sup> <sup>L</sup>* is the flux in the downwind cell, and *Ff* is defined in the face *f* as:

$$F\_f = U\_f \delta t a\_f \tag{13}$$

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

$$F\_{i \ast}^{\mathcal{A}} = F\_{i \ast}^{\mathcal{H}} - F\_{i \ast}^{\mathcal{L}} \tag{14}$$

To make this stable, a correction factor, *q*, is needed to modify the fluxes values. Finally, a value for fluid fraction in next time step is defined as:

$$\alpha\_{i}^{n+1} = \alpha\_{i}^{\*} - \frac{\left(q\_{i+\underline{\lambda}^{\*}}F\_{i+\underline{\lambda}^{\*}}^{A} - q\_{i-\underline{\lambda}^{\*}}F\_{i-\underline{\lambda}^{\*}}^{A}\right)}{\delta\mathbf{x}}\tag{15}$$

#### **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 *β* in each cell of the numerical domain. The position of this line segment in each cell is defined such that the area reconstruct‐ ed from the line and the perimeter of the cell is equal to the amount of volume of fluid, *α*. The geometry of the polygon from this reconstruction is used to calculate the flux transferred from the cell faces.

Assuming that *αi*, *<sup>j</sup>* is predefined in every cell, the first step is to calculate first‐order upwind fluxes. Then, the Youngs exiting fluxes of every cell can be calculated by considering the values in each cell. To do so, the angle *β*, between free surface and x‐coordinate, must be calculated. Different methods can be used for calculation of *β*. One of them is first using the gradient function of fluid volume fraction for defining unit normal vector of the free surface, and then calculating *β* [20]. The method of defining the normal vector, however, can affect the accura‐ cy of the final results. This formulation for uniform grid is as follows:

**Figure 7.** Four possible positions for free surface in Youngs' method [19].

$$\mathbf{n}\_{i,j}^{\times} = \frac{\mathbf{a}\_{i+1,j+1} + \mathbf{2}\mathbf{a}\_{i+1,j} + \mathbf{a}\_{i+1,j-1} - \mathbf{a}\_{i-1,j+1} - \mathbf{2}\mathbf{a}\_{i-1,j} - \mathbf{a}\_{i-1,j-1}}{\delta \mathbf{x}} \tag{16}$$

$$n\_{i,j}^y = \frac{\alpha\_{i+1,j+1} + 2\alpha\_{i,j+1} + \alpha\_{i-1,j+1} - \alpha\_{i+1,j-1} - 2\alpha\_{i,j-1} - \alpha\_{i-1,j-1}}{\delta y} \tag{17}$$

Using components of the normal unit vector, the angle *β* can be calculated as:

$$\beta = \tan^{-1} \left( \frac{-n\_{\varepsilon}}{n\_{\gamma}} \right) \quad \left( -\pi < \beta \le \pi \right) \tag{18}$$

The angle *γ* is also defined as:

**Figure 6.** Three adjacent typical cells in FCT method.

372 Numerical Simulation - From Brain Imaging to Turbulent Flows

for fluid fraction in next time step is defined as:

**3.4. Youngs' method**

the cell faces.

Assuming that *αi*, *<sup>j</sup>*

where *<sup>F</sup> <sup>L</sup>* is the flux in the downwind cell, and *Ff* is defined in the face *f* as:

*F Ut f ff* = d a

of the previous step. This is the difference between upwind and downwind fluxes as:

1 11 2 22

( 11 11 ) 1 \* 22 22

*ii ii n*

*i i*

 a

a

Thereafter, an anti‐diffusive flux is needed to be defined (*FL*) in order to correct the diffusion

To make this stable, a correction factor, *q*, is needed to modify the fluxes values. Finally, a value

*qF qF x*

d

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 *β* in each cell of the numerical domain. The position of this line segment in each cell is defined such that the area reconstruct‐ ed from the line and the perimeter of the cell is equal to the amount of volume of fluid, *α*. The geometry of the polygon from this reconstruction is used to calculate the flux transferred from

fluxes. Then, the Youngs exiting fluxes of every cell can be calculated by considering the values in each cell. To do so, the angle *β*, between free surface and x‐coordinate, must be calculated. Different methods can be used for calculation of *β*. One of them is first using the gradient

*A A*

(13)

*A HL FFF i ii* + ++ = - (14)

++ -- <sup>+</sup> - = - (15)

is predefined in every cell, the first step is to calculate first‐order upwind

$$\gamma = \tan^{-1} \left( \frac{\delta \chi}{\delta \chi} \tan \beta \right) \quad \left( 0 \le \gamma \le \frac{\pi}{2} \right) \tag{19}$$

It is possible to set 0 ≤*γ* ≤90 by rotating the cell. Therefore, there are only four possible positions for the free surface, which are depicted in **Figure 7**.

What is behind this conclusion is as follows:

( ) ( ) 1 2 1 2 1 2 1 2 4 tan 1 tan cot 1 cot *if if Case I elseif Case II else Case IV else if Case I elseif Case III else Case IV* g p a g a g a g a g < £ £ - £ £ -



( )

 g

1 2

1 tan

£ -

1 2

*Case I elseif*

*Case II*

a

*else*

a

*else if* tan

 g

4

< £

a

*if if*

374 Numerical Simulation - From Brain Imaging to Turbulent Flows

g p

> 1 2

*Case I elseif*

£

*Case III*

a

*Case IV*

*else*

**Case I Case II**

*st* <sup>0</sup> <sup>0</sup>

*sb* 2*α*cot*γ* <sup>1</sup>

*if Ut* >0 *if Utδt* ≤(1−*sr*)*δy Ft* =0 *else*

*Ft* <sup>=</sup> <sup>1</sup>

*Fr* =*αδxδy*

<sup>2</sup> *Urδt*(2<sup>−</sup> *Urδ<sup>t</sup>*

*sbδx*)*srδ<sup>y</sup>*

*if Ur* >0 *if Urδt* ≥*sbδx*

*else*

*Fr* <sup>=</sup> <sup>1</sup>

*sr* 2*α*tan*γ α* + <sup>1</sup>

*sl* <sup>0</sup> *α* − <sup>1</sup>

<sup>2</sup> *Utδt* −(1−*sr*)*δy* 2cot*γ*

*Case IV*

cot

 g

£ -

( )

 g

<sup>2</sup> tan*γ*

<sup>2</sup> tan*γ*

*Ft* =0

*Ft* <sup>=</sup> <sup>1</sup>

*else*

*if Utδt* ≤(1−*sr*)*δy*

*elseif Utδt* ≤(1−*sl*)*δy*

*Fr* <sup>=</sup>*Urδt*(*srδ<sup>y</sup>* <sup>−</sup> <sup>1</sup>

<sup>2</sup> *Utδt* −(1−*sr*)*δy* 2cot*γ*

<sup>2</sup> *Urδt*tan*γ*)

*Ft* =*Utδtδx* −(1−*α*)*δxδy*

1 2

1 cot


$$F\_b = \mathcal{U}\_b \delta t \delta \chi - \frac{1}{2} \{\mathcal{U}\_b \delta t - \mathbf{s}\_l \delta y\}^2 \text{coty}$$


**Table 1.** Calculation of exiting flux in Youngs' method.

Four side fractions (*sl* , *sb*, *sr*, *st*) for up, right, down, and left faces can be calculated with the selection of the free surface position in a cell. Thereafter, flow fluxes can be geometrically computed for each face (*Fl* , *Fb*, *Fr*, *Ft*) based on these side fractions. More details are present‐ ed in **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 inter‐ face 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:

$$\vec{m}\_{i,j} = -\left(\frac{\nabla \alpha\_{i,j}}{\sqrt{\left(\nabla \alpha\_{x,i,j}\right)^2 + \left(\nabla \alpha\_{y,i,j}\right)^2}}\right) \tag{20}$$

where (∇*α*) is defined as:

$$\left(\nabla a\right)\_{i,j} = \frac{1}{4} \Big[ \left(\nabla a\right)\_{i-\underset{K^{\circ}\cdot\boldsymbol{\aleph}^{\circ}\cdot\boldsymbol{K}}} + \left(\nabla a\right)\_{i-\underset{K^{\circ}\cdot\boldsymbol{\aleph}^{\circ}\cdot\boldsymbol{K}}} + \left(\nabla a\right)\_{i+\underset{K^{\circ}\cdot\boldsymbol{\aleph}^{\circ}\cdot\boldsymbol{K}}} + \left(\nabla a\right)\_{i+\underset{K^{\circ}\cdot\boldsymbol{\aleph}^{\circ}\cdot\boldsymbol{K}}} \right] \tag{21}$$

#### Free Surface Flow Simulation Using VOF Method http://dx.doi.org/10.5772/64161 377

$$\begin{aligned} \left(\nabla\alpha\right)\_{i+\overset{j}{\mathbf{X}},j+\overset{i}{\mathbf{X}}} &= \frac{1}{2} \left(\frac{\alpha\_{i+1,j}-\alpha\_{i,j}}{\mathbf{x}\_{i+1,j}-\mathbf{x}\_{i,j}} + \frac{\alpha\_{i+1,j+1}-\alpha\_{i,j+1}}{\mathbf{x}\_{i+1,j+1}-\mathbf{x}\_{i,j+1}}\right)\mathbf{\overset{\leftharpoonup}{i}} + \\ &\frac{1}{2} \Big(\frac{\alpha\_{i,j+1}-\alpha\_{i,j}}{\mathbf{y}\_{i,j+1}-\mathbf{y}\_{i,j}} + \frac{\alpha\_{i+1,j+1}-\alpha\_{i+1,j}}{\mathbf{y}\_{i+1,j+1}-\mathbf{y}\_{i+1,j}}\Big)\mathbf{\overset{\leftharpoonup}{j}} \end{aligned} \tag{22}$$

**Figure 8.** Different positioning of the interface for (0≤θ≤π/ 4).

**Case III Case IV**

and there is no flux calculation for negative velocities into the cell.

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

*i j*

*n*

*if Ul*

*else*

*Fl* =*Ul*

*Fl* =*Ul*

, *sb*, *sr*, *st*) for up, right, down, and left faces can be calculated with the

, *Fb*, *Fr*, *Ft*) based on these side fractions. More details are present‐

selection of the free surface position in a cell. Thereafter, flow fluxes can be geometrically

ed in **Table 1**. In this table, positive value is set for velocities towards the outer edges of a cell,

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 inter‐ face 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

> ( )( ) , , 2 2

a

( ) ( ) 1 1 ( ) 1 1 ( ) 1 1 ( ) 1 1 2 2 2 2 2 2 2 2 , ,,,,

 a

4 *i j i j i j i j i j*


a

æ ö ç ÷ <sup>Ñ</sup> = - Ñ +Ñ è ø

, , , , *i j*

*xi j yi j*

 a

> a

ë û

<sup>r</sup> (20)

 a (21)

*δt* ≥(1−*st*)*δx*

*δt*(*sl δy* + <sup>1</sup> <sup>2</sup> *Ul*

*δtδy* −(1−*α*)*δxδy*

*δt*tan*γ*)

*δt* −(1−*sb*)*δx* 2tan*γ*

*δtδy* −(1−*α*)*δxδy*

*if Ul* >0 *if Ul*

*δt* ≤*sbδx*

*δt* ≤*stδx*

376 Numerical Simulation - From Brain Imaging to Turbulent Flows

**Table 1.** Calculation of exiting flux in Youngs' method.

using Youngs' least square method as:

where (∇*α*) is defined as:

1

aa

*Fl* =0 *elseif Ul*

> *Fl* <sup>=</sup> <sup>1</sup> <sup>2</sup> *Ul*

*Fl* =*Ul*

*else*

Four side fractions (*sl*

computed for each face (*Fl*

where *θ* is the angle between the normal vector and the horizontal coordinate varies be‐ tween zero and 2*π*. For *θ* in the first one‐eighth space (0≤*θ* ≤*π* / 4), eight different conditions are possible for the position of free surfaces as illustrated in **Figure 8** [21]. All other situa‐ tions 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 defin‐ ing surface unit normal vector using volume fraction of fluid in each cell. To do this, ex‐ treme values of *α*lim,1 and *α*lim,2 are determined as:

$$
\alpha\_{\rm lim,1} = \frac{n\_{\rm min}}{2n\_{\rm max}} \quad , \quad \alpha\_{\rm lim,2} = 1 - \alpha\_{\rm lim,1} \tag{23}
$$

in which *v*(*i* + 1, *j* −<sup>1</sup> 2 ) and *n*max are calculated as:

$$n\_{\min} = \min\left( \left| n\_x \right|, \left| n\_y \right| \right) \quad , \quad n\_{\max} = \max\left( \left| n\_x \right|, \left| n\_y \right| \right) \tag{24}$$

*α*lim,1 and *α*lim,2 are shown in **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.

**Figure 9.** Various positions of an interface in a cell.
