**4. Calculating surface tension**

2 1

 aa

*D DD*

230

*f DD D*

% %% %

a

ï

 aa

a

384 Numerical Simulation - From Brain Imaging to Turbulent Flows

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

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

relation with mass conservation.

**3.7. Integrated methods**

1 2

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

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

1 1

î < <

<sup>ì</sup> < > <sup>ï</sup> =- + < < <sup>í</sup>

% %%

2

(44)

*D*

a

%

 a

0, 1

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 meth‐ od 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 equa‐ tions, 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 perpen‐ dicular 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:

$$\mathbf{F}\_s = \sigma \int\_S \kappa \mathbf{n} \,\delta\left(\mathbf{x} - \mathbf{x}\_s\right) dS \tag{45}$$

where *σ*, *κ*, **n**, *δ*(**x**−**x***s*) and **X***<sup>s</sup>* are surface tension force, curvature of the surface, normal unit vector, Dirac delta function, and the position of a given base point on the interface *S*, respec‐ tively. This equation has been discretized for numerical methods of two‐phase fluids.

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

$$\mathbf{F}\_s = \sigma \mathbf{\kappa} \mathbf{n} \,\boldsymbol{\delta}\_s = -\nabla \cdot \mathbf{T}, \quad \mathbf{T} = -\sigma \left(\mathbf{I} - \mathbf{n} \otimes \mathbf{n}\right) \boldsymbol{\delta}\_s \tag{46}$$

where **I** and **T** are unit tensor and tangential tension tensor of the interface respectively.

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**):

**Figure 12.** Cells on the interface in contact with left and right cells of the free surface [38].

$$P\_L = P\_{1L}\frac{I\_L}{\Delta \mathbf{y}} + P\_{2L}\frac{\Delta \mathbf{y} - I\_L}{\Delta \mathbf{y}} = P\_{1L}H\_L + P\_{2L}\left(\mathbf{l} - H\_L\right) = P\_{2L} + H\_L\left(P\_{1L} - P\_{2L}\right) \tag{47}$$

where *PL* is the mean pressure on the left side of the given cell. In the same way, mean pressure can be calculated for other faces of the cell as follows:

$$P\_R = P\_{2R} + H\_R \left(P\_{1R} - P\_{2R}\right) \tag{48}$$

$$P\_T = P\_{2T} + H\_T \left(P\_{1T} - P\_{2T}\right) \tag{49}$$

$$P\_B = P\_{2B} + H\_B \left(P\_{1B} - P\_{2B}\right) \tag{50}$$

where *H* is a dimensionless number which shows the position of free surface in different directions. For faces completely immersed in the main fluid, H is equal to 1, and for those completely in the secondary fluid, *H* is zero. For other cases, *H* varies between zero and unity. Note that for 3D models, it is just needed to replace edge faces with surfaces faces.

On the other hand, the change in pressure *Ps* in every point of the interface is calculated as:

$$P\_s = P\_1 - P\_2 = \sigma \mathbf{x} \tag{51}$$

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

$$P\_k = P\_{2k} + H\_k \sigma \mathbf{\dot{x}} \tag{52}$$

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:

$$\mathbf{F}\_s = H \sigma \mathbf{\kappa} \mathbf{m} \tag{53}$$

where **F***s* is the surface tension force vector.

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**):

386 Numerical Simulation - From Brain Imaging to Turbulent Flows

**Figure 12.** Cells on the interface in contact with left and right cells of the free surface [38].

*y y*

can be calculated for other faces of the cell as follows:

D -

1 2 1 2 () ( ) 2 12 1 *L L LL L LL L L L L L L l yl P P P PH P H P H P P*

where *PL* is the mean pressure on the left side of the given cell. In the same way, mean pressure

where *H* is a dimensionless number which shows the position of free surface in different directions. For faces completely immersed in the main fluid, H is equal to 1, and for those completely in the secondary fluid, *H* is zero. For other cases, *H* varies between zero and unity.

Note that for 3D models, it is just needed to replace edge faces with surfaces faces.

D D (47)

*P P HP P R R RR R* =+ - 2 12 ( ) (48)

*P P HP P T T TT T* =+ - 2 12 ( ) (49)

*P P HP P B B BB B* =+ - 2 12 ( ) (50)

= + = + - =+ -

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 equiva‐ lent 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 *y*(*x*) is fitted to the nine neighboring cells of the inter‐ face. In this way, summation of the volume of fluid of three cells located in a column creates a value for *y*. By fitting parabolic curve to the values of three neighboring columns and then two times differentiating of these values, one can define the curvature of the surface. This method, however, suffers from low accuracy and some limiting conditions. Another method is presented by Chorin [39] in which a circle is defined based on trial and error in order to satisfy the nine‐cell set in the best possible way. Thereafter, curvature of the surface can be defined on the basis of the calculated circle. However the main problem of this method is its dependency on the several times of trials and errors in order to calculate the best way of estimating the circle, the thing that makes it practically inefficient. Ashgriz and Poo [40] proposed another method in which a parabolic curve is fitted to the 9 or 25 neighboring cells. This method is more accurate than the previous one; however, its applications are very limited to specific cases.

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

$$\mathbf{x} = -\nabla \cdot \mathbf{n} = -\nabla \cdot \frac{\nabla \phi}{|\nabla \phi|}\tag{54}$$

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:

$$\text{div}\_{\mathbf{x},i\prec\mathbf{y},j} = 2\frac{\phi\_{i+1,j} - \phi\_{i,j}}{\Delta\mathbf{x}\_{i+1} + \Delta\mathbf{x}\_{i}}, \quad \text{div}\_{\mathbf{y},i\prec\mathbf{y},j} = \frac{\phi\_{i+1,j} - \phi\_{i,j-1} + \phi\_{i+1,j+1} - \phi\_{i+1,j-1}}{\Delta\mathbf{y}\_{j+1} + 2\Delta\mathbf{y}\_{j} + \Delta\mathbf{y}\_{j-1}}\tag{55}$$

$$\ln \mathbf{n}\_{\boldsymbol{x},i+\boldsymbol{\xi},j} = \frac{\text{div}\_{\boldsymbol{x},i+\boldsymbol{\xi},j}}{\left| \mathbf{n}\_{i+\boldsymbol{\xi},j} \right|}, \quad \left| \mathbf{n}\_{i+\boldsymbol{\xi},j} \right| = \sqrt{\text{div}\_{\boldsymbol{x},i+\boldsymbol{\xi},j}^2 + \text{div}\_{\boldsymbol{\nu},i+\boldsymbol{\xi},j}^2} \tag{56}$$

$$\mathbf{x}\_{i,j} = -\frac{\mathbf{n}\_{\mathbf{x},i+\mathbf{y}}\mathbf{y}\_{i,j} - \mathbf{n}\_{\mathbf{x},i-\mathbf{y}}\mathbf{y}\_{i,j}}{\Delta \mathbf{x}\_i} - \frac{\mathbf{n}\_{\mathbf{y},i,j+\mathbf{y}}\mathbf{y}\_i - \mathbf{n}\_{\mathbf{x},i,j-\mathbf{y}}\mathbf{y}\_i}{\Delta \mathbf{y}\_j} \tag{57}$$

**Figure 13.** Position of the normal vectors of the surface in the cell faces.

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:

$$\mathbf{F}\_s = \sigma \mathbf{\kappa} \mathbf{n} \,\boldsymbol{\delta}\_s = \sigma \mathbf{\kappa} \mathbf{n} \frac{\left| \nabla \overline{\boldsymbol{\alpha}} \right|}{\left[ \boldsymbol{a} \right]} \tag{58}$$

where *δs* is the delta Dirac function which is infinite on the interface and zero otherwise, and its integration is equal to unity. The bar sign in Eq. (58) shows a smoothed (or filtered) value of volume fraction. The bracket sign shows the difference between maximum and minimum of volume of the fluid fraction. In CSF method, this function is estimated using |∇*α*¯ | / *α* . Some other references use the following equation to improve the results' accuracy [41]:

$$\mathbf{F}\_s = \sigma \kappa \mathbf{n} \,\boldsymbol{\delta}\_s = \sigma \kappa \mathbf{n} \,\frac{\left| \nabla \overline{a} \right|}{\left[ a \right]} \frac{\rho}{\left[ \rho \right]} \tag{59}$$

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:

$$F\_{\mathbf{x}\_{i},i+\underline{\mathbf{y}}\_{i},j} = 2\sigma\kappa \frac{\overline{\alpha}\_{i+1,j} - \overline{\alpha}\_{i,j}}{\Delta\mathbf{x}\_{i+1} + \Delta\mathbf{x}\_{i}} \frac{\rho\_{i+\underline{\mathbf{y}}\_{i},j}}{\rho\_{\max} - \rho\_{\min}}\tag{60}$$

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

$$F\_{x\_{i},i} \chi\_{\mathbb{1}^{J}} = 2\sigma\kappa \frac{\overline{\alpha}\_{i+1,j} - \overline{\alpha}\_{i,j}}{\Delta x\_{i+1} + \Delta x\_{i}} H\_{i+\underline{\mathsf{X}},j} \tag{61}$$

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