Conservative Level Set

For cases involving multiple phases of viscous fluid, the different phases are tracked by a set of functions of the form $φ_{i} (x, t)$ . The function $φ_{i} (x, t)$ denotes the regions of the domain occupied by phase $i$ at an arbitrary spatial point $x$ and time $t$ .

The location of each phase can be described by the distance function $φ_{i} (x, t)$ according to:

Figure 1. EQUATION_DISPLAY

φ_{i} (x, t) = {\begin{matrix} 1, x \in Ω_{i} \\ 0, x \notin Ω_{i} \end{matrix}

(1074)

$Ω_{i}$ is the domain occupied by the $i$ th phase. This scalar field for each phase of the $N$ -phase system is then advected with the flow velocity field v according to:

Figure 2. EQUATION_DISPLAY

\frac{\partial φ_{i}}{\partial t} + \nabla\cdot (v φ_{i}) = 0, i = 1, 2, 3 \dots N

(1075)

The numerical discretization of the above transport equations always results in an error in $φ_{i} (x, t)$ that ultimately brings about signiﬁcant volume conservation errors, so that further attention and techniques are required to improve volume conservation. To address this issue, Olsson and Kreiss [273] and Olsson et al. [274] developed a conservative level set (CLS) function that is capable of implicitly tracking the interface. In this case, this function is deﬁned as:

Figure 3. EQUATION_DISPLAY

ϕ_{i} (x, t) = \frac{1}{2} (\tanh [\frac{φ_{i} (x, t)}{2 ϵ}] + 1)

(1076)

where $ϕ_{i} (x, t)$ is the conservative level set function and $ϵ$ determines the interface thickness. With this proper mapping of $φ_{i}$ , the interface can now be represented by the 0.5 iso-contour of $ϕ_{i} (x, t)$ . Furthermore, to obtain accurate values for $ϕ_{i} (x, t)$ , Olsson et al. [274] proposed a two-step advection and reinitialization scheme for a two-phase flow system where the following pure advection is considered:

Figure 4. EQUATION_DISPLAY

\frac{\partial ϕ}{\partial t} + \nabla\cdot (v ϕ) = 0

(1077)

The resulting value of the above equation is then considered as the initial guess for the re-initialization step, which can be written as:

Figure 5. EQUATION_DISPLAY

\frac{\partial ϕ}{\partial τ} + \nabla\cdot (ϕ [1 - ϕ] n - ϵ (\nabla ϕ \cdot n) n) = 0

(1078)

where the first term (first-order derivative in $τ$ ,which is pseudo-time, not the physical time) is solely introduced to guarantee the well-posedness of the above equation. The above scheme can be extended to an $N$ -phase system by adopting the strategy recently proposed by Howard and Tartakovsky [249] to avoid the possibility of formation of overlaps and voids between the phases. This can be achieved by introducing Lagrange multipliers to the right-hand side of the re-initialization equations such that the following constraints are satisfied throughout the domain:

Figure 6. EQUATION_DISPLAY

\begin{matrix} \sum_{i = 1}^{N} ϕ_{i}^{0} (x) = 1 \\ \sum_{i = 1}^{N} \frac{\partial ϕ_{i}^{n} (x)}{\partial τ} = 0 \end{matrix}

(1079)

where $ϕ_{i}^{n}$ is the phase indicator function for the $i$ th phase at time $n$ . The Lagrange multipliers can be determined analytically and thus the two-step conservative level-set method can be extended to an $N$ -phase system according to:

Figure 7. EQUATION_DISPLAY

\begin{matrix} \frac{\partial ϕ_{i}}{\partial t} + v \cdot \nabla ϕ_{i} = 0, i = 1, 2 \dots N \\ \frac{\partial ϕ_{i}}{\partial τ} + \nabla\cdot f_{i} = 0, i = 1, 2 \dots N \end{matrix}

(1080)

where, in the above, the flow is incompressible (that is, $\nabla\cdot v = 0$ ) and $f_{i}$ is:

Figure 8. EQUATION_DISPLAY

f_{i} = ϕ_{i} (1 - ϕ_{i}) n_{i} - ϵ (\nabla ϕ_{i} \cdot n_{i}) n_{i} - \frac{ϕ_{i}^{p}}{\sum_{i = 1}^{N} ϕ_{i}^{p}} \sum_{i = 1}^{N} ϕ_{j} (1 - ϕ_{j}) n_{j}

(1081)

where $p$ , the Lagrange multiplier moment, is a positive integer with a default value of 2. To obtain a variational form, first multiply the above equation by a weight function ${\tilde{ϕ}}_{i}$ and then integrate over the entire domain:

Figure 9. EQUATION_DISPLAY

\begin{matrix} \int_{Ω} {\tilde{ϕ}}_{i} \frac{\partial ϕ_{i}}{\partial t} d Ω + \int_{Ω} {\tilde{ϕ}}_{i} v \cdot \nabla ϕ_{i} d Ω + \sum_{e = 1}^{n_{e l}} τ_{supg} v \cdot \nabla {\tilde{ϕ}}_{i} [\frac{\partial ϕ_{i}}{\partial t} + v \cdot \nabla ϕ_{i}] d Ω_{e} = 0, \forall {\tilde{ϕ}}_{i} \in H_{0}^{1} \\ \int_{Ω} {\tilde{ϕ}}_{i} \frac{\partial ϕ_{i}}{\partial t} d Ω + \int_{Ω} {\tilde{ϕ}}_{i} v \cdot f_{i} d Ω = 0, \forall {\tilde{ϕ}}_{i} \in H_{0}^{1} \end{matrix}

(1082)

where

$Ω$ is volume and $Ω_{e}$ is the volume of the $e$ th element.
$n_{e l}$ is the number of elements.
$H_{0}^{1} (Ω)$ a Sobolev space.
This space consists of square integrable functions for which the first derivatives are also square integrable, with the magnitude zero at the boundaries with essential (Dirichlet) boundary conditions.

Notice that the third term in the Eqn. (1082) is the conventional upwinding scheme (SUPG) used to handle the advection term. For the second equation, integration by parts gives:

Figure 10. EQUATION_DISPLAY

\int_{Ω} {\tilde{ϕ}}_{i} \frac{\partial ϕ_{i}}{\partial τ} d Ω - \int_{Ω} \nabla {\tilde{ϕ}}_{i} \cdot f_{i} d Ω + \int_{\partial Ω} {\tilde{ϕ}}_{i} f_{i} \cdot e_{n} d \partial Ω = 0, \forall {\tilde{ϕ}}_{i} \in H_{0}^{1}

(1083)

where $e_{n}$ is the unit normal vector pointing outward to the boundary $\partial Ω$ enclosing the domain $Ω$ . The third (boundary) term in the above equation can be completely neglected on the wall where there is zero flux of the field variable. By contrast, a Dirichlet boundary condition for the phase field variable ( $ϕ_{i} = 0$ or $ϕ_{i} = 1$ ) for mass and velocity inlets can be imposed.