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

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:

1, , 1, , 1 1, 1 1, 1

+ + - ++ +-

fff


**<sup>n</sup>** (56)

+- +- - - <sup>=</sup> - - D D (57)

 f

1 1 1

+ + -

*i i j jj*

, , 2 2 , , , ,, ,,

> 11 11 22 22 , , , , ,, ,,

> > *i j*

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

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

a

**Fn n** = = (58)

a

Ñ

[ ] *s s*

sk d sk

*xi j xi j yi j xi j*

<sup>1</sup> 1 11 <sup>2</sup> 2 22

*nn nn x y*

*xi j i j xi j yi j*

= =+ **n**

*i j ij i j ij i j i j*

2 , <sup>2</sup>

*div div x x y yy*

*div <sup>n</sup> div div* <sup>+</sup> + + ++

1 1 2 2

> 1 2

1 2

+

,

k

*i j*

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

results.

method as:

*i j*

, , *xi j*

, , , ,

+ +

ff

388 Numerical Simulation - From Brain Imaging to Turbulent Flows

*xi j yi j*

In the four‐point method, free surface (as illustrated in **Figure 14**) is approximated using a continuous function *f* ( ⋅ )∈*C* <sup>2</sup> (*x*1, *x*4) (a set of functions with continuous second derivation) so that the distance between the function and the points is minimized according to Eq. (62), and variations of curvature are bounded based on Eq. (63). In this case, the radius curvature can be calculated as in Eq. (64).

$$\inf\_{f(\cdot) \in c^2(\mathbf{x}\_1, \mathbf{x}\_4)} I\left(f\left(\cdot\right)\right) = \sum\_{i=1}^4 \left| f\left(\mathbf{x}\_i\right) - \mathbf{y}\_i \right| \tag{62}$$

$$\begin{cases} \left| \kappa \left( \mathbf{x} \right) - C\_1 \right| \le \varepsilon & \mathbf{x} \in \left( \mathbf{x}\_1, \mathbf{x}\_2 \right] \\ \left| \kappa \left( \mathbf{x} \right) - C\_2 \right| \le \varepsilon & \mathbf{x} \in \left( \mathbf{x}\_2, \mathbf{x}\_3 \right] \\ \left| \kappa \left( \mathbf{x} \right) - C\_3 \right| \le \varepsilon & \mathbf{x} \in \left( \mathbf{x}\_3, \mathbf{x}\_4 \right] \\ C\_i \ge 0 & i = 1, 2, 3 \end{cases} \tag{63}$$

$$\kappa \left( x \right) = \frac{\left| \ddot{f} \left( x \right) \right|}{\sqrt{1 + \dot{f}^2 \left( x \right)}} \tag{64}$$

**Figure 14.** Free surface modeling in FPM.

where *ε* is a small arbitrary given number.

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

$$\inf I\left(P\_n\left(\cdot\right)\right) = \sum\_{i=1}^4 \left| P\_n\left(x\_i\right) - y\_i \right| \tag{65}$$

$$\kappa\_n(\mathbf{x}) = \frac{\left| \ddot{P}\_n(\mathbf{x}) \right|}{\sqrt{1 + \dot{P}\_n^{-2}(\mathbf{x})}} \tag{66}$$

where *κ* is the curvature of the free surface.

It is supposed that Q is the set of *f* ( ⋅ )∈*C* <sup>2</sup> (*x*1, *x*4) such that Eq. (62) is feasible and Q(n) is the set of *Pn*(.) such that Eq. (64) is feasible. Then, the following theorem proves that the se‐ quence of solutions for Eqs. (65) and (66) converges to the solution of Eqs. (62) and (64) as *n* goes towards infinity.

**Theorem 1**: if *<sup>η</sup>* <sup>=</sup>*in f <sup>Q</sup>I*( *<sup>f</sup>* (.)) and *η*(*n*)=*in f Qn I*(*Pn*(.)) then *η* =*limn*→*∞η*(*n*).

() ( ) <sup>2</sup> ( ) ( ) ( ) 1 4

, <sup>1</sup> inf *i i f C xx <sup>i</sup> If fx y* × Î <sup>=</sup>

, . .

<sup>ì</sup> -£ Î <sup>ï</sup> ïï -£ Î <sup>í</sup> ï -£ Î

*S t*

390 Numerical Simulation - From Brain Imaging to Turbulent Flows

**Figure 14.** Free surface modeling in FPM.

where *ε* is a small arbitrary given number.

where *κ* is the curvature of the free surface.

with unknown constant coefficients. Therefore, we have:

k

k

k

k=

ï

0 1,2,3 *<sup>i</sup>*

( ) ( )

1 *f x <sup>x</sup>*

+

&&

*C i*

ïî ³ =

4

1 1 2 2 2 3 3 3 4

> ( ) 2

*f x*

In this method, the desired function is approximated using an n‐degree polynomial function

( ) ( ) ( ) <sup>4</sup> 1 inf *<sup>n</sup> ni i i IP P x y* =

( ) ( )

*P x <sup>x</sup>*

*n*

k= ( ) <sup>2</sup> 1 *n*

+ &&

*P x*

*n*

,

,

( ) ( ] ( ) ( ] ( ) ( ]

*x C x xx*

*x C x xx x C x xx*

 e

> e

 e

×= - å (62)

& (64)

×= - å (65)

& (66)

(63)

**Proof**. It is obvious that *Q*(1)⊂*Q*(2)⊂...⊂*Q*, then *η*(1)≥*η*(2)≥...≥*η*. Therefore, {*η*(*n*)} is a non‐ increasing and bounded sequence. Then it converges to a number called *ξ*. Set *<sup>W</sup>* <sup>=</sup>∪ *n*=1 *∞ Q*(*n*); therefore, inf*Qn I*(*Pn*( ⋅ )) =*ξ*. Since *W* ⊂*Q*, then *ξ* ≥*η*. By the properties of infimum, for every *ε* >0, there exists *f* ( ⋅ )∈*Q* such that:

$$
\mathfrak{n} < I\left(f\left(\cdot\right)\right) < \mathfrak{n} + \mathfrak{s} \tag{67}
$$

As we have *f* ( ⋅ )∈*C* <sup>2</sup> (*x*1, *<sup>x</sup>*4), there is a set of polynomials such that {*Pn*( <sup>⋅</sup> )}, {*P*˙ *<sup>n</sup>*( <sup>⋅</sup> )}, {*<sup>P</sup>* ¨ *<sup>n</sup>*( ⋅ )} are uniformly convergent. Therefore, there exists a natural number such that for every *n* ≥ *N* we have:

$$\left\| P\_n(\cdot) - f\_n(\cdot) \right\|\_\diamond < \delta \tag{68}$$

$$\left\|\dot{P}\_{\boldsymbol{n}}(\cdot) - \dot{f}\_{\boldsymbol{n}}(\cdot)\right\|\_{\infty} < \delta \tag{69}$$

$$\left\| \ddot{P}\_{\pi}(\cdot) - \ddot{f}\_{\pi}(\cdot) \right\|\_{\circ} < \delta \tag{70}$$

Now, it is claimed that there is an *N*<sup>1</sup> ≥ *N* such that for *i* =1, 2, 3, 4, the relation |*κn*(*x*)−*Ci* | ≤*ε*; ∀ *x* ∈(*xi* , *xi*+1 is true. Since otherwise for every *n* ≥ *N* , there is a |*κn*(*x*)−*Ci* | ≤*ε*; ∀ *x* ∈(*xi* , *xi*+1 . Therefore, lim*<sup>n</sup>*→*<sup>∞</sup>* |*κn*(*x*)−*Ci* | >*ε* which contradicts the assump‐ tion of *f* ( ⋅ )∈*Q*. Thus, *PN* ( ⋅ )∈*Q*(*N*1)⊂*W* ⊂*Q*. Based on what was mentioned, | *I*(*PN*<sup>1</sup> ( <sup>⋅</sup> )) <sup>−</sup> *<sup>I</sup>*( *<sup>f</sup>* ( <sup>⋅</sup> ))| <sup>&</sup>lt;*ε* or *I*(*PN*<sup>1</sup> ( ⋅ )) < *I*( *f* ( ⋅ )) + *ε* <*η* + 2*ε* or *η* ≤*ξ* + 2*ε*; therefore, *ξ* =*η*, or lim*n*→*∞η*(*n*)=*η*.

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

$$\inf I\left(P\_n(\cdot)\right)$$

$$\begin{aligned} \left\| \int\_{x\_1}^{x\_2} \left\| \kappa\_n(\mathbf{x}) - C\_1 \right\| - \varepsilon + \left\| \kappa\_n\left(\mathbf{x}\right) - C\_1 \right\| - \varepsilon \right\| \mathrm{d}\mathbf{x} &= \mathbf{0} \\ \left\| \int\_{x\_1}^{x\_2} \left\| \kappa\_n\left(\mathbf{x}\right) - C\_2 \right\| - \varepsilon + \left\| \kappa\_n\left(\mathbf{x}\right) - C\_2 \right\| - \varepsilon \right\| \mathrm{d}\mathbf{x} &= \mathbf{0} \\ \left\| \int\_{x\_1}^{x\_2} \left( \kappa\_n\left(\mathbf{x}\right) - C\_3 \right) - \varepsilon + \left\| \kappa\_n\left(\mathbf{x}\right) - C\_3 \right\| - \varepsilon \right\| \mathrm{d}\mathbf{x} &= \mathbf{0} \\ C\_i \succeq \mathbf{0} & i = 1, 2, 3, n = 1, 2, \dots \end{aligned} \tag{71}$$

or

inf ( ) ( )× *<sup>n</sup> I P*

$$S.I = \begin{bmatrix} \frac{h\_1}{2} \left[ E\_{1u} \left( \mathbf{x}\_1 \right) + E\_{1u} \left( \mathbf{x}\_1 + h \right) + \dots + E\_{1u} \left( \mathbf{x}\_1 + m\_1 h \right) \right] = 0\\ \frac{h\_2}{2} \left[ E\_{2u} \left( \mathbf{x}\_2 \right) + E\_{2u} \left( \mathbf{x}\_2 + h \right) + \dots + E\_{2u} \left( \mathbf{x}\_2 + m\_2 h \right) \right] = 0\\ \frac{h\_3}{2} \left[ E\_{3u} \left( \mathbf{x}\_3 \right) + E\_{3u} \left( \mathbf{x}\_3 + h \right) + \dots + E\_{3u} \left( \mathbf{x}\_3 + m\_3 h \right) \right] = 0\\ C\_i \ge 0 & i = 1, 2, 3. \qquad n = 1, 2, \dots \end{cases} \tag{72}$$

such that *Ein*(*x*)= ∥*κn*(*x*)−*Ci* | −*ε* + ∥*κn*(*x*)−*Ci* | −*ε* ∥.

Now, the intervals *x*1, *x*<sup>2</sup> , *x*2, *x*3 and *x*3, *x*4 are divided into three equal sections, *m*1, *m*2, *m*3, respectively. This means that *h*<sup>1</sup> <sup>=</sup> *<sup>x</sup>*<sup>2</sup> <sup>−</sup> *<sup>x</sup>*<sup>1</sup> *<sup>m</sup>*<sup>1</sup> , *h*<sup>2</sup> <sup>=</sup> *<sup>x</sup>*<sup>3</sup> <sup>−</sup> *<sup>x</sup>*<sup>2</sup> *<sup>m</sup>*2 and *h*<sup>3</sup> <sup>=</sup> *<sup>x</sup>*<sup>4</sup> <sup>−</sup> *<sup>x</sup>*<sup>3</sup> *<sup>m</sup>*<sup>3</sup> . Thus, using a numerical integration method such as trapezoidal rule, the problem above can be reformulated as:

$$\inf \sum\_{i=1}^{4} \left| P\_u \left( x\_i \right) - y\_i \right| $$

$$S.t = \begin{bmatrix} \frac{h\_1}{2} \left[ E\_{1u} \left( \mathbf{x}\_1 \right) + E\_{1u} \left( \mathbf{x}\_1 + h \right) + \dots + E\_{1u} \left( \mathbf{x}\_1 + m\_1 h \right) \right] = \mathbf{0} \\\\ \frac{h\_2}{2} \left[ E\_{2u} \left( \mathbf{x}\_2 \right) + E\_{2u} \left( \mathbf{x}\_2 + h \right) + \dots + E\_{2u} \left( \mathbf{x}\_2 + m\_2 h \right) \right] = \mathbf{0} \\\\ \frac{h\_3}{2} \left[ E\_{3u} \left( \mathbf{x}\_3 \right) + E\_{3u} \left( \mathbf{x}\_3 + h \right) + \dots + E\_{3u} \left( \mathbf{x}\_3 + m\_3 h \right) \right] = \mathbf{0} \\\\ C\_i \ge \mathbf{0} & i = 1, 2, 3. \qquad n = 1, 2, \dots$$

(74)

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

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

( ) ( )

 ek

*n n*

*n n*

*n n*

1 1

*x C x C dx*

0

 e

> e

 e

... 0

... 0

*<sup>m</sup>*2 and *h*<sup>3</sup> <sup>=</sup> *<sup>x</sup>*<sup>4</sup> <sup>−</sup> *<sup>x</sup>*<sup>3</sup>

... 0

... 0

0

(71)

(72)

(73)

*<sup>m</sup>*<sup>3</sup> . Thus, using a numerical

2 2

*x C x C dx*

<sup>0</sup> . .

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

ï - -+ - - =

ï - -+ - - =

3 3

*x C x C dx*

0 1,2,3, 1,2,...

inf ( ) ( )× *<sup>n</sup> I P*

() ( ) ( )

11 11 11 1

<sup>ì</sup> é + + + + + ù= <sup>ï</sup> ë û <sup>ï</sup>

*<sup>h</sup> E x E x h E x mh*

*n n n*

() ( ) ( )

22 22 22 2

... <sup>0</sup> . <sup>2</sup>

<sup>ï</sup> é + + + + + ù= ë û <sup>=</sup> <sup>í</sup>

*<sup>h</sup> E x E x h E x mh*

*<sup>h</sup> E x E x h E x mh*

*n n n*

*n n n*

*C in*

î ³ ==

33 33 33 3

<sup>ï</sup> é + + + + + ù= ë û <sup>ï</sup>

0 1,2,3. 1,2,...

() ( ) ( )

Now, the intervals *x*1, *x*<sup>2</sup> , *x*2, *x*3 and *x*3, *x*4 are divided into three equal sections, *m*1, *m*2, *m*3,

( ) <sup>4</sup> 1 inf =

() ( ) ( )

11 11 11 1

<sup>ì</sup> é + + + + + ù= <sup>ï</sup> ë û <sup>ï</sup>

*<sup>h</sup> E x E x h E x mh*

() ( ) ( )

22 22 22 2

... <sup>0</sup> . <sup>2</sup>

<sup>ï</sup> é + + + + + ù= ë û <sup>=</sup> <sup>í</sup>

*<sup>h</sup> E x E x h E x mh*

*<sup>h</sup> E x E x h E x mh*

*n n n*

*n n n*

*C in*

î ³ ==

33 33 33 3

<sup>ï</sup> é + + + + + ù= ë û <sup>ï</sup>

0 1,2,3. 1,2,...

() ( ) ( )

å *ni i* -

*Px y*

*<sup>m</sup>*<sup>1</sup> , *h*<sup>2</sup> <sup>=</sup> *<sup>x</sup>*<sup>3</sup> <sup>−</sup> *<sup>x</sup>*<sup>2</sup>

integration method such as trapezoidal rule, the problem above can be reformulated as:

*i*

*n n n*

2

k

k

k

*x*

ì

ï ï ï

1 3

*x x*

ò

2 4

*x x*

ò

í ï

*S t*

392 Numerical Simulation - From Brain Imaging to Turbulent Flows

or

3

*x i*

ò

ï ï

1

2

2

ï

ï

ï

respectively. This means that *h*<sup>1</sup> <sup>=</sup> *<sup>x</sup>*<sup>2</sup> <sup>−</sup> *<sup>x</sup>*<sup>1</sup>

*S t*

3

2

*i*

such that *Ein*(*x*)= ∥*κn*(*x*)−*Ci* | −*ε* + ∥*κn*(*x*)−*Ci* | −*ε* ∥.

1

2

2

ï

ï

ï

*S t*

3

2

*i*

( ) ( )

 ek

( ) ( )

 ek

*C in*

î ³ ==

In this method as illustrated in **Figure 15**, the main goal is to find a function as *f* ( ⋅ )∈*C* <sup>2</sup> (*x*1, *x*4) within *x* ∈ *xj* , *x <sup>j</sup>*+1 such that the distance between the function *f* (*x*) and the line *L <sup>i</sup>* (*x*) which connects the given points (*xi* , *yi* ) and (*xi*+1, *yi*+1) is minimized for *i* =1, 2, 3. Thus, variation of curvature is bounded according to Eq. (74):

$$\begin{aligned} \text{inf } I\left(f\left(\cdot\right)\right) &= \sum\_{i=1}^{j} \int\_{x\_i}^{x\_{i+1}} \left| f\left(x\right) - L\_i\left(x\right) \right| dx \\\\ \text{s.t. } & \left| \kappa\left(x\right) - C\_1 \right| \le \varepsilon \qquad x \in \left(x\_1, x\_2\right] \\\\ \text{s.t. } & \left| \kappa\left(x\right) - C\_2 \right| \le \varepsilon \qquad x \in \left(x\_2, x\_3\right] \\\\ & \left| \kappa\left(x\right) - C\_3 \right| \le \varepsilon \qquad x \in \left(x\_3, x\_4\right] \\\\ & \text{s.t. } \varepsilon \qquad x \in \left(x\_3, x\_4\right] \end{aligned}$$

0 1,2,3 *<sup>i</sup>*

*C i*

ïî ³ =

**Figure 15.** Free surface modeling by TLM.

Similar to what was discussed in the four‐point method, in this method the function *f* ( ⋅ ) can be replaced with an n‐degree polynomial, *Pn*( ⋅ ) as below:

$$\inf I\left(P\_{\pi}\left(\cdot\right)\right) = \sum\_{i=1}^{3} \int\_{x\_{i}}^{x\_{i+1}} \left| P\_{\pi}\left(x\right) - L\_{i}\left(x\right) \right| d\mathbf{x}$$

$$\begin{aligned} \left| \begin{array}{cc} \left| \kappa\_{n}(\mathbf{x}) - C\_{1} \right| \leq \varepsilon & \mathbf{x} \in \left[ \mathbf{x}\_{1}, \mathbf{x}\_{2} \right] \\ \left| \kappa\_{n}(\mathbf{x}) - C\_{2} \right| \leq \varepsilon & \mathbf{x} \in \left( \mathbf{x}\_{2}, \mathbf{x}\_{3} \right] \\ \left| \kappa\_{n}(\mathbf{x}) - C\_{3} \right| \leq \varepsilon & \mathbf{x} \in \left( \mathbf{x}\_{3}, \mathbf{x}\_{4} \right] \\ C\_{i} \geq 0 & i = 1, 2, 3, \quad n = 1, 2, \dots \end{array} \end{aligned} \tag{75}$$

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

$$\inf I\left(P\_{\boldsymbol{\pi}}\left(\cdot\right)\right) = \sum\_{i=1}^{3} \frac{h\_i}{2} \begin{pmatrix} \left|P\_{\boldsymbol{\pi}}\left(\mathbf{x}\_{\boldsymbol{1}}\right) - L\_i\left(\mathbf{x}\_{\boldsymbol{1}}\right)\right| + 2\left|P\_{\boldsymbol{\pi}}\left(\mathbf{x}\_{\boldsymbol{1}} + h\right) - L\_i\left(\mathbf{x}\_{\boldsymbol{1}} + h\right)\right| + \dots \\ + 2\left|P\_{\boldsymbol{\pi}}\left(\mathbf{x}\_{\boldsymbol{1}} + (m\_i - 1)h\right) - L\_i\left(\mathbf{x}\_{\boldsymbol{1}} + (m\_i - 1)h\right)\right| + \\ \left|P\_{\boldsymbol{\pi}}\left(\mathbf{x}\_{\boldsymbol{1}} + m\_i h\right) - L\_i\left(\mathbf{x}\_{\boldsymbol{1}} + m\_i h\right)\right| \end{pmatrix}$$

$$S.S.t = \begin{bmatrix} \frac{h\_1}{2} \left[ E\_{1s} \left( \mathbf{x}\_1 \right) + E\_{1s} \left( \mathbf{x}\_1 + h \right) + \dots + E\_{1s} \left( \mathbf{x}\_1 + m\_1 h \right) \right] = 0\\ \frac{h\_2}{2} \left[ E\_{2s} \left( \mathbf{x}\_2 \right) + E\_{2s} \left( \mathbf{x}\_2 + h \right) + \dots + E\_{2s} \left( \mathbf{x}\_2 + m\_2 h \right) \right] = 0\\ \frac{h\_3}{2} \left[ E\_{3s} \left( \mathbf{x}\_3 \right) + E\_{3s} \left( \mathbf{x}\_3 + h \right) + \dots + E\_{3s} \left( \mathbf{x}\_3 + m\_3 h \right) \right] = 0\\ C\_i \ge 0 & i = 1, 2, 3. \qquad n = 1, 2, \dots \end{cases} \tag{76}$$

The steps of using the above equations are as follows:

Step 1: Read *ε*, and set *n* =1.

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

Step 5: If | *In* − *In*−<sup>1</sup> | >*ε*1 then go to step 4, else *In* is the final answer.
