**3.6. Higher order differencing schemes**

Another method to reconstruct the interface between two fluids is to discretize the convec‐ tion 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 diffu‐ sion 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:

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

$$
\tilde{\alpha}\_{\dot{\gamma}} = \gamma\_{\dot{\gamma}} \tilde{\alpha}\_{\dot{\gamma}\_{\text{CC}}} + (1 - \gamma\_{\dot{\gamma}}) \tilde{\alpha}\_{\dot{\gamma}\_{\text{I}\dot{\gamma}}} \tag{25}
$$

where *α*˜ *<sup>f</sup> CBC* is an index for transport bound which defines a bound as follows:

$$
\tilde{\alpha}\_{f\_{\rm conc}} = \begin{cases}
\min\left\{1, \frac{\tilde{\alpha}\_D}{c\_{f\_{\rm c}}}\right\} & \text{when } 0 \le \tilde{\alpha}\_D \le 1 \\
\tilde{\alpha}\_D & \text{when } \tilde{\alpha}\_D \le 0, \ \tilde{\alpha}\_D > 1
\end{cases} \tag{26}
$$

where *α*˜ *<sup>D</sup>* is defined based on the following equation:

min min , , max , ( ) *x y* max ( ) *x y n nn n nn* = = (24)

When unit normal vector of a surface is defined, the true position of the interface can be easily

Another method to reconstruct the interface between two fluids is to discretize the convec‐ tion 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

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 diffu‐ sion 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

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:

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

*α*lim,1 and *α*lim,2 are shown in **Figure 9**.

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

**3.6. Higher order differencing schemes**

simulation.

interface as follows:

determined using volume of fluid in each cell.

378 Numerical Simulation - From Brain Imaging to Turbulent Flows

$$
\tilde{\alpha}\_D = \frac{\alpha\_D - \alpha\_U}{\alpha\_A - \alpha\_U} \tag{27}
$$

$$
\tilde{\boldsymbol{\alpha}}\_{\boldsymbol{\tilde{\alpha}}\_{\boldsymbol{\tilde{\alpha}}}} = \begin{cases}
\min\left\{\frac{8\boldsymbol{c}\_{\boldsymbol{\tilde{\alpha}}}\tilde{\boldsymbol{\alpha}}\_{\boldsymbol{\tilde{\alpha}}} + (\mathbb{I} - \boldsymbol{c}\_{\boldsymbol{\gamma}})(6\tilde{\boldsymbol{\alpha}}\_{\boldsymbol{\tilde{\alpha}}} + 3)}{8}, \tilde{\boldsymbol{\alpha}}\_{\boldsymbol{\tilde{\alpha}}\_{\boldsymbol{\tilde{\alpha}}}}\right\} & \text{when } 0 \le \tilde{\boldsymbol{\alpha}}\_{\boldsymbol{\tilde{\alpha}}} \le 1 \\
\tilde{\boldsymbol{\alpha}}\_{\boldsymbol{\tilde{\alpha}}} & \text{when } \quad \tilde{\boldsymbol{\alpha}}\_{\boldsymbol{\tilde{\alpha}}} < 0, \ \tilde{\boldsymbol{\alpha}}\_{\boldsymbol{\tilde{\alpha}}} > 1
\end{cases} \tag{28}
$$

More details on determining *α*˜ *<sup>f</sup>* can be derived based on the combined schemes such as Normalized Variable Diagram (NVD) and discussed by Ubbink and Issa [26]. Accordingly, a local boundedness is defined for *α*˜ *<sup>f</sup> UQ*. In Eq. (25), *γf* is a weight factor which is defined as:

$$\gamma\_{\gamma} = \min \left\{ k\_{\gamma} \frac{\cos(2\theta\_{\gamma}) + 1}{2}, 1 \right\} \tag{29}$$

where *kγ* is a constant, usually set to unity, and *θf* is the angle between free surface normal vector, (∇*α*)*D*, and the vector *d* ⇀ *<sup>f</sup>* such that it connects the center of the adjacent cells based on **Figure 10**, and defined as:

$$\theta\_f = \arccos \left| \frac{(\nabla \alpha)\_D \vec{d}\_f}{\left| (\nabla \alpha)\_D \right| \left| \vec{d}\_f \right|} \right| \tag{30}$$

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.

**Figure 10.** Calculation of the upwind value for an arbitrary mesh.

#### *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* is calculated using a weighting factor β<sup>f</sup> as follows:

$$\alpha\_f = \left(1 - \beta\_f\right)\alpha\_D + \beta\_f\alpha\_A \tag{31}$$

where

$$\mathcal{B}\_f = \frac{\tilde{\mathcal{a}}\_f - \tilde{\mathcal{a}}\_D}{1 - \tilde{\mathcal{a}}\_D} \tag{32}$$

**α˜f** and **α˜D** can be calculated by Eqs. (25) and (27), respectively.

#### *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, *γ<sup>f</sup>* , of Eq. (29) can be redefined as:

$$\gamma\_f = \min\left(\left(\theta\_f\right)^m, 1\right) \tag{33}$$

where *θf* was previously defined in Eq. (30). For *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*

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

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,

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* is calculated using a

*f f D fA* =- + (31)

% (32)

considering the angle between the interface and the direction of motion [27].

( ) 1

 b a ba

> 1 *f D*

a

*D*

a <sup>=</sup> - % %

a

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

*f*

b-

a

and **α˜D** can be calculated by Eqs. (25) and (27), respectively.

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

method, the weighing factor, *γ<sup>f</sup>* , of Eq. (29) can be redefined as:

linearization, which makes it less accurate for 3D modeling reconstructions.

**Figure 10.** Calculation of the upwind value for an arbitrary mesh.

380 Numerical Simulation - From Brain Imaging to Turbulent Flows

as follows:

*3.6.2 THOR scheme*

weighting factor β<sup>f</sup>

where

**α˜f**

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]:

$$
\tilde{\boldsymbol{\alpha}}\_{\prime}^{\*\*} = \begin{cases}
\tilde{\boldsymbol{\alpha}}\_{D} & \tilde{\boldsymbol{\alpha}}\_{D} < 0 \\
2\tilde{\boldsymbol{\alpha}}\_{D} & 0 < \tilde{\boldsymbol{\alpha}}\_{D} < 0.5 \\
1 & 0.5 < \tilde{\boldsymbol{\alpha}}\_{D} < 1 \\
\tilde{\boldsymbol{\alpha}}\_{D} & 1 < \tilde{\boldsymbol{\alpha}}\_{D}
\end{cases} \tag{34}
$$

$$
\tilde{\alpha}\_{f}^{\*} = \begin{cases}
\tilde{\alpha}\_{f}^{\*\*} & c\_{f} < 0.3 \\
\tilde{\alpha}\_{D}^{\*} + \left(\tilde{\alpha}\_{f}^{\*\*} - \tilde{\alpha}\_{D}\right) \frac{0.7 - c\_{f}}{0.4} & 0.3 < c\_{f} < 0.7 \\
\tilde{\alpha}\_{D}^{\*} & 0.7 < c\_{f}
\end{cases} \tag{35}
$$

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

$$
\tilde{\alpha}\_{\uparrow} = \tilde{\alpha}\_{\uparrow}^\* \sqrt{\cos \theta} + \tilde{\alpha}\_{\diamond} (1 - \sqrt{\cos \theta}) \tag{36}
$$

In this way, *αf* is discretized based on the neighboring cells.

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:

$$\tilde{\boldsymbol{\alpha}}\_{f,\,\,FRCS}^{\*} = \begin{cases} \mathbf{3}\tilde{\boldsymbol{\alpha}}\_{\,\,D}\tilde{\boldsymbol{\alpha}}\_{\,\,D} \le \bigvee\_{\mathcal{S}} \\ \tilde{\boldsymbol{\alpha}}\_{\,\,D} + \frac{1}{4}\bigvee\_{\mathcal{S}} < \tilde{\boldsymbol{\alpha}}\_{\,\,D} < \tilde{\mathcal{Y}}\_{\,\,\,A} \\ 1 \bigvee\_{\mathcal{A}} < \tilde{\boldsymbol{\alpha}}\_{\,\,D} < 1 \\ \tilde{\boldsymbol{\alpha}}\_{\,\,D}\tilde{\boldsymbol{\alpha}}\_{\,\,D} \le 0 \,\,or\,\tilde{\boldsymbol{\alpha}}\_{\,\,D} \ge 1 \end{cases} \tag{37}$$

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:

$$
\tilde{\alpha}\_{f} = \tilde{\alpha}\_{f}^{\*} + \left(\tilde{\alpha}\_{f}^{\*} - \tilde{\alpha}\_{D}\right) \frac{0.7 - Cn}{0.7 - 0.3} \tag{38}
$$

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 cos<sup>4</sup> *θ* term is designed instead of cos*θ* in Eq. (36) as follows:

$$
\tilde{\alpha}\_f = \tilde{\alpha}\_{f,\text{sup}}^\* \cos^4 \theta + \tilde{\alpha}\_{f,\text{static}}^\* \left(1 - \cos^4 \theta\right) \tag{39}
$$

where *α*˜\* *f* , *sup* and *α*˜\* *<sup>f</sup>* , *stoic* are calculated as follows:

$$
\tilde{\boldsymbol{\tilde{\alpha}}}\_{f,\text{sup}}^{\*} = \begin{cases}
\tilde{\boldsymbol{\tilde{\alpha}}}\_{D} \tilde{\boldsymbol{\tilde{\alpha}}}\_{D} \leq 0 \\
10 < \tilde{\boldsymbol{\alpha}}\_{D} < 1 \\
\tilde{\boldsymbol{\alpha}}\_{D} \tilde{\boldsymbol{\alpha}}\_{D} \geq 1
\end{cases} \qquad \tilde{\boldsymbol{\tilde{\alpha}}}\_{f,\text{sup}}^{\*} = \begin{cases}
\tilde{\boldsymbol{\tilde{\alpha}}}\_{D} \tilde{\boldsymbol{\tilde{\alpha}}}\_{D} \leq 0 \\
\frac{1}{2} + \frac{1}{2} \tilde{\boldsymbol{\tilde{\alpha}}}\_{D} 0 < \tilde{\boldsymbol{\alpha}}\_{D} < \frac{1}{2} \\
\frac{3}{8} + \frac{3}{4} \tilde{\boldsymbol{\tilde{\alpha}}}\_{D} \frac{1}{2} < \tilde{\boldsymbol{\tilde{\alpha}}}\_{D} < \frac{5}{6} \\
\frac{1}{8} < \tilde{\boldsymbol{\tilde{\alpha}}}\_{D} < 1 \\
\frac{1}{6} < \tilde{\boldsymbol{\tilde{\alpha}}}\_{D} < 1 \\
\tilde{\boldsymbol{\tilde{\alpha}}}\_{D} \tilde{\boldsymbol{\tilde{\alpha}}}\_{D} \geq 1
\end{cases} \tag{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*

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

> <sup>1</sup> <sup>3</sup> <sup>8</sup> 1 1 3 8 4 4

*D D*

% %

a a

<sup>ì</sup> £ <sup>ï</sup>

a

% % %

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

increases. The Courant number, *Cn*, in HRIC method can be written as follows:

 aa

( ) \* \* 0.7 0.7 0.3 *f f fD*

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

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 cos<sup>4</sup>

( ) \* 4\* 4

*cos*

q

(39)

a

, , 1 *f f sup f stoic* %% % = +-

*<sup>f</sup>* , *stoic* are calculated as follows:

ï

aa

*D D*

%% %

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

<sup>ï</sup> + << <sup>=</sup> <sup>í</sup> <sup>ï</sup> < < <sup>ï</sup> <sup>ï</sup> £ ³ <sup>î</sup>

31 1 <sup>4</sup>

a

%

*D DD D or*

 a

0 1

*Cn*


(37)

*θ*

 a

\* ,

a

382 Numerical Simulation - From Brain Imaging to Turbulent Flows

*f FBICS*

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

aa

for CICSAM when the *Cn* is equal to unity.

term is designed instead of cos*θ* in Eq. (36) as follows:

where *α*˜\*

*f* , *sup*

and *α*˜\*

a *cos* qa

scheme as:

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

$$\frac{\partial \mathcal{\gamma}}{\partial t} + \nabla \cdot \left( \mathbf{U} \boldsymbol{\gamma} \right) + \nabla \cdot \left( \mathbf{U}' \boldsymbol{\gamma} \left( \mathbf{l} - \boldsymbol{\gamma} \right) \right) = \mathbf{0} \tag{41}$$

where **U***<sup>r</sup>* is a velocity field for compressing the free surface. This artificial term is activated only in the presence of free surface because of having the term *γ*(1−*γ*). The solution to this equation is bounded from zero to unity with the inter‐gamma scheme. Eq. (40) can be rewritten as:

$$
\frac{\partial \hat{\boldsymbol{\gamma}}}{\partial t} + \nabla \cdot \left( \boldsymbol{\phi} \left[ \boldsymbol{\mathcal{I}} \boldsymbol{\mathcal{I}} \right]\_{f(\boldsymbol{\phi}, \boldsymbol{S})} \right) + \nabla \cdot \left( \boldsymbol{\phi}\_{\boldsymbol{\phi}} \left[ \boldsymbol{\mathcal{I}} \right]\_{f(\boldsymbol{\phi}\_{\boldsymbol{\phi}}, \boldsymbol{S})} \right) = \mathbf{0} \tag{42}
$$

where *<sup>ϕ</sup>* <sup>=</sup>**S**⋅**<sup>U</sup>** *<sup>f</sup>* is the volume flux and *ϕrb* =(1−*γ*) *f* (−*ϕ <sup>r</sup>* ,*<sup>S</sup>* )*ϕ<sup>r</sup>* . Eq. (43) is used to calculate ∅*<sup>r</sup>* as:

$$\phi^{r} = K\_c n^\* \max \left\{ \frac{\left| \boldsymbol{n}^\* \phi \right|}{\left| \boldsymbol{S} \right|^2} \right\} \tag{43}$$

where *n* \* is the free surface normal unit vector and *Kc* is an adjustable coefficient, with an appropriate value of 1.5, which defines the compression rate of free surface. The inter‐gamma scheme is also similar to CICSAM and is based on a donor‐acceptor equation and normal‐ ized variable diagram, NVD. The formulation is as:

$$
\tilde{\alpha}\_{f} = \begin{cases}
\tilde{\alpha}\_{\mathcal{D}} & \tilde{\alpha}\_{\mathcal{D}} < 0, \ \tilde{\alpha}\_{\mathcal{D}} > 1 \\
1 & \mathcal{Y}\_{2} < \tilde{\alpha}\_{\mathcal{D}} < 1
\end{cases} \tag{44}
$$

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

**Figure 11.** 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 recon‐ struct 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 satisfy‐ ing 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, integra‐ tion 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.
