### 3.3.7 Link to other discretisation spaces

As pointed out in the introduction, the developed spaces can be linked to the $H(\text{div})$-conforming virtual discretisation setting presented in [100] and outlined in Section 3.2.4, where the considered space reads

$$ V_{2,k}^{face}(K) = \{ v \in H(\text{div}, K) \cap H(\text{curl}, K) \quad \text{s.t.} \quad v \cdot n|_f \in \mathbb{P}_k(f) \forall f \in \partial K, \\ \nabla(\nabla \cdot v) \in \nabla \mathbb{P}_{k-1}(K), \text{ and } \nabla \times v \in \mathbb{P}_{k-1}(K) \}. $$

Indeed, restricting the setting on the boundary to polynomial functions, that is, introducing

$$ \tilde{V}_{2,k}(K) = \{ v \in H(\text{div}, K) \cap H(\text{curl}, K) \quad \text{s.t.} \ v \cdot n|_f \in \mathbb{P}_k(f) \quad \text{and} \ v|_f \in \times_{i=1}^d \mathbb{Q}_{\zeta_i(k+1, k, \dots, k)}(f) \quad \forall f \in \partial K, \\ \nabla(\nabla \cdot v) \in \nabla \mathbb{P}_{k-1}(K), \text{ and } \nabla \times v \in \mathbb{P}_{k-1}(K) \}. $$

it can in particular be shown the following inclusion.

**Proposition 3.3.1** Relation between $\mathbb{H}_k(K)$ and $\tilde{V}_{2,k}(K)$

Let us consider the $k^{\text{th}}$ order space $\mathbb{H}_k(K)$, $k \in \mathbb{N} \setminus \{0\}$, constructed from the refinement pattern $(l_1, l_2) = (0, k)$ and $(m_1, m_2) = (k, -1)$, that is, define

$$ \mathbb{H}_k(K) = \{ u \in H^1(K), u|_{\partial K} \in \mathcal{H}_0(\partial K), \Delta u \in \mathbb{Q}_k(K) \}^d \\ + x \{ u \in H^1(K), u|_{\partial K} \in \mathcal{H}_k(\partial K), \Delta u \in \mathbb{Q}_{[-1]}(K) \}. $$

Then,

$$ v \in \tilde{V}_{2,k} \Rightarrow v \in \mathbb{H}_k(K), $$

**Proof.**

Let us show that for any $k \in \mathbb{N} \setminus \{0\}$, any element of $\tilde{V}_{2,k}(K)$ can be recast as an element of $\mathbb{H}_k(K)$ constructed from the coefficients $(l_1, l_2) = (0, k)$ and $(m_1, m_2) = (k - 1, -1)$. We start by showing that $\mathbb{H}_k(K)$ is more general in terms of regularity, before deriving the inclusion relations by using from the block structure of the investigated spaces.

* (Regularity) Any element $v \in \tilde{V}_{2,k}(K)$ belongs to $H^1(K)$. Indeed, $v \in H(\text{div}, K) \cap H(\text{rot}, K) \subset (H^1(K))^d$. Furthermore, as any element $w$ of $\mathbb{H}_k(K)$ writes $w = u_1 + x u_2$ with $u_1, u_2$ both having for regularity $H^1(K)$, $w \in H_{\text{loc}}^1(K)$. And since $K$ is compact and bounded, $H_{\text{loc}}^1(K) = H^1(K)$. Thus, the regularity asked for any element $v$ to be in $\tilde{V}_{2,k}(K)$ is stricter than the one asked for any element $w$ to be in $\mathbb{H}_k(K)$.

* (Inclusion on the boundary) Any $v \in \tilde{V}_{2,k}(K)$ is a polynomial and satisfies $v \cdot n|_f \in \mathbb{P}_f(K)$ on every face $K$ of the polygon $K$. But $v \cdot n|_f$ is nothing else than the linear combination $\sum_{i=1}^d v_i n_i$ of the coordinate-wise functions $v_i$ with the normal's coefficients $n_i$. Thus, each polynomial $v_i$ has no choice
but to live in the space $\mathbb{P}_{k, \dots, k, k+1, k, \dots, k}(K)$, where the $k+1$ is in the $i^{\mathrm{th}}$ position. The space $\times_{i=1}^d \mathbb{P}_{k, \dots, k, k+1, k, \dots, k}(K)$ is indeed the smallest space that is polynomial (required by the linear combination) and that contains all the functions $v$ such that $v \cdot n|_f \in \mathbb{P}_f(K)$. Note that allowing a higher degree in the $i^{\mathrm{th}}$ variable is required as on each face $f \in \partial K$, $x \cdot n \equiv c$ for some constant $c$. And $\times_{i=1}^d \mathbb{P}_{k, \dots, k, k+1, k, \dots, k}(K) \subset \mathbb{P}_0(K) + x \mathbb{P}_k(K)$, which is exactly the structure of the $\mathbb{H}_k(K)$ space on the boundary.

• (Inclusion within the element) For any $v \in \tilde{V}_{2,k}(K)$, it holds

$$ \left\{ \begin{array}{l} \nabla(\nabla \cdot v) \in \nabla(\mathbb{P}_{k-1}(K)) \\ \nabla \times v \in \mathbb{P}_{k-1}(K) \end{array} \right. $$

Thus, it comes $\underbrace{\nabla(\nabla \cdot v)}_{\in \nabla(\mathbb{P}_{k-1}(K))} - \underbrace{\nabla \times (\nabla \times v)}_{\in \nabla(\mathbb{P}_{k-2}(K))} \in \nabla(\mathbb{P}_{k-1}(K))$, which implies $\nabla^2 v \in \nabla(\mathbb{P}_{k-1}(K))$ and

$$ \nabla^2 v = \begin{pmatrix} \sum_{i=1}^d \frac{\partial^2 v_1}{\partial x_i^2} \\ \vdots \\ \sum_{i=1}^d \frac{\partial^2 v_d}{\partial x_i^2} \end{pmatrix} \in \begin{pmatrix} \mathbb{P}_{k-2, k-1, \dots, k-1}(K) \\ \vdots \\ \mathbb{P}_{k-1, k-1, \dots, k-2}(K) \end{pmatrix} $$

Thus, we have naturally

$$ \nabla^2 v \in \{u, \begin{pmatrix} \Delta u_1 \\ \vdots \\ \Delta u_d \end{pmatrix}, \Delta u_i \in \mathbb{Q}_{k-1}(K) \}. $$

• (Wrap-up) Let $v \in V_{2,k}(K)$. Then $v$ can be decomposed as

$$ v = (\hat{v}_1 + \hat{v}_2)|_{\tilde{K}} + (\bar{v}_1 + \bar{v}_2)|_{\partial K}, $$

where

$$ \left\{ \begin{array}{ll} \bar{v}_1 \in (\mathbb{P}_0(\partial K))^d & \text{represents the constant part of } v \text{ on the boundary} \\ \bar{v}_2 \in (\mathbb{P}_k(\partial K) \setminus \mathbb{P}_0(\partial K))^d & \text{represents the higher parts of } v \text{ on the boundary} \end{array} \right. $$

and

$$ \left\{ \begin{array}{l} \hat{v}_1 = v|_{\tilde{K}} - p_{[1], \dots, [1]} \\ \bar{v}_2 = p_{[1], \dots, [1]} \end{array} \right. $$

for any $p_{[1], \dots, [1]}$ belonging to $x\mathbb{P}_0$. Doing so, we get

$$ \left\{ \begin{array}{l} \bar{v}_{1|\partial K} \in (\mathbb{P}_0(\partial K))^d \\ \hat{v}_{1|\tilde{K}} = v|_{\tilde{K}} - p_{[1], \dots, [1]} \end{array} \right. $$
for the lower boundary characterisation, while for the higher one it comes

$$
\left\{
\begin{array}{l}
\bar{v}_{2|\partial K} \quad (\mathbb{P}_k(\partial K) \setminus \mathbb{P}_0(\partial K))^d \\
\hat{v}_{2|\check{K}} = p_{[1], \dots, [1]}
\end{array}
\right.
$$

There, we get straightforwardly from the previous paragraphs that

$$
\left\{
\begin{array}{l}
v_1 \in H^1(K) \\
v_1 \cdot n_{|\partial K} \in \mathbb{P}_0(\partial K) \subset \mathbb{Q}_0(\partial K) \\
\Delta v_{1|\check{K}} = \underbrace{v_{|\check{K}}}_{\in \nabla(\mathbb{P}_{k-1}(K))} - \underbrace{\Delta(p)}_{=0} \in \mathbb{P}_{k-1}(K) \subset \mathbb{Q}_{k-1}(K),
\end{array}
\right.
$$

and

$$
\left\{
\begin{array}{l}
v_2 \in H^1(K) \\
v_2 \cdot n_{|\partial K} = v_{|\partial K} \cdot n - v_1 \cdot n_{|\partial K} \quad \in \mathbb{P}_k(\partial K) \subset \mathbb{Q}_k(\partial K) \\
\Delta v_{2|\check{K}} = \Delta(x c) = x \Delta(c) = 0 \quad \in \mathbb{P}_{[-1]}(K) \subset \mathbb{Q}_{[-1]}(K)
\end{array}
\right.
$$

for some constant $c \in \mathbb{P}_0(K)$. Writing $v_2$ as $v_2 = xw$ on the boundary with $w \in x \mathbb{P}_{k, \dots, k}(\partial K) \subset \times_{i=1}^d \mathbb{P}_{k, \dots, k, k+1, k, \dots, k}(\partial K) \setminus \mathbb{P}_0(\partial K)$, it comes further

$$
\left\{
\begin{array}{l}
w \in H^1(K) \\
w \cdot n_{|\partial K} \quad \in \mathbb{P}_k(\partial K) \subset \mathbb{Q}_k(\partial K) \\
\Delta w_{|\check{K}} = 0 \quad \in \mathbb{P}_{-1}(K) \subset \mathbb{Q}_{[-1]}(K)
\end{array}
\right.
$$

Taking without loss of generality $c = 1$, considering $v = v_1 + xw$ and setting
$(l_1, l_2) = (0, k)$, $(m_1, m_2) = (k - 1, -1)$, we get:

$v_1 \in \mathcal{A}_k, \quad v_2 = xw \in x \mathcal{B}_k,$

and therefore,

$v = v_1 + v_2 \in \mathbb{H}_k(K),$

concluding the proof.

Even though this inclusion holds for any $k \in \mathbb{N} \setminus \{0\}$ and defines $\tilde{V}_k$ as a discretisation subspace of $\mathbb{H}_k(K)$, the selected degrees of freedom in the framework of [100] are different from those we considered in Section 3.3.4.1, providing a less refined and less flexible element. Indeed, the Virtual element setting [100] allows less freedom in the inner characterisation than the proposed framework. Furthermore, contrarily to the more general setting presented in this section, the normal component on the boundary belongs at least to $\mathbb{P}_1(K)$ for any admissible space. A lowest order element guaranteeing a constant value for the discretised functions' normal component on the boundaries thus cannot be natively defined in the conforming Virtual elements framework [100], but is subject to a slight modification (see e.g. [112]).
### 3.3.8 Numerical results

Let us now confirm the properties obtained along the theoretical developments by exploring the basis functions associated to the main element $Ia$ and its variant $Ib$ presented in the previous sections, for both the general and reduced frameworks. We particularly focus on the normal component $\varphi \cdot n|_{\partial K}$ of representative basis functions $\varphi$ on the boundary of the element $K$. Note here that in practice, those basis functions have been constructed by tuning a natural basis of the space $H_k(K)$ towards selected sets of degrees of freedom through a transfer matrix.

As an example, we consider the non-convex nine-edges polygon presented in Figure 3.7, on which the elements are built. In all the results, the polynomial projectors used in the definition of the degrees of freedom were chosen as Hermite polynomials, as we have observed experimentally that this improves the conditioning of the linear system.

<figure><img src="image_43.png" /><figcaption>**Figure 3.7.** Left: Considered polygon. Right: Normal component of a representative internal basis functions plotted on every edge.</figcaption></figure>

#### 3.3.8.1 Elements constructed in a general setting

Let us start by considering the space and elements as described in the Section 3.3.6.1. There, we start by investigating the behaviour of the internal basis functions, before discussing the results obtained on the boundary for the normal ones.

**Investigating the first configuration of degrees of freedom** Having built the elements based on the degrees of freedom configuration $I$, we represented the normal component of the associated internal basis functions in the right side of Figure 3.7. As desired, the normal component of basis functions corresponding to internal degrees of freedom vanish on the boundary, as shown on the right graph.
In order to study the behaviour of the normal basis functions on the boundaries, we then have considered the case $k = 2$, where we expect five basis function to have a quadratic normal component. We have plotted in the leftmost side of Figure 3.8 the normal component of one of the normal basis functions associated to the element $Ib$, and can first observe that its support is contained on one single edge.

<figure><img src="image_44.png" /><figcaption>**Figure 3.8.** Left: normal component of one representative of the normal basis functions for the element $Ib$ and $k = 2$ along the edges. Middle: normal component of all the functions generated from the edge number 5, plotted on the edge number 5. Right: as middle, for the element $Ia$.</figcaption></figure>

We then have plotted all the basis functions associated to the fifth edge on Figure 3.8, for both the configurations $Ib$ (middle) and $Ia$ (right). There, one can first notice that their normal components, given in the middle graph of Figure 3.8, are polynomial of degree $k \le 2$, that together generate $\mathbb{P}_2(\mathbb{R})$.

<figure><img src="image_45.png" /><figcaption>**Figure 3.9.** Degeneration of a degenerating normal basis functions' representative in the case $k = 2$, for the element $Ib$. Left: normal component on all the boundaries. Right: internal behaviour of the basis functions.</figcaption></figure>

Observing further, it appears that the normal component of two basis are vanishing, that is $\varphi \cdot n = \varphi_1 \cdot n_{5,x} + \varphi_2 \cdot n_{5,y} = 0$. This generates this straight line equal to zero in the graph. Indeed, those two basis functions characterise the coordinate-wise freedom $\varphi_1 \cdot n_{5,x} \neq 0$ and $\varphi_2 \cdot n_{5,y} \neq 0$. This additional freedom is not reflected through the global term $\varphi \cdot n$ as addressed when discussing the
impact of the degrees of freedom's choice in Section 3.3.4.1. This comes from the fact that only three basis functions are required to generate $\mathbb{P}_2(\mathbb{R})$, where the global component $\varphi \cdot n$ lives. The two components $\varphi_1 n_{5,x}$ and $\varphi_2 n_{5,y}$ of the vector $\varphi \circ n$ are compensating themselves. Those functions are nevertheless regular within the polygon $K$ and not identically vanishing on $K$, as shown in Figure 3.9, where the left graph represents the value of the normal component $\varphi \cdot n$ on the boundary and where the right graph represents the components $\varphi_1$ and $\varphi_2$ on the element $K$. They can therefore be reclassified as internal basis functions. Note that this can be suppressed when using the reduced setting, as one can observe below in Section 3.3.8.2.

<figure><img src="image_46.png" /><figcaption>**Figure 3.10.** Scaling of the non-vanishing basis functions generated from the edge number 5 when $k=0$. From left to right: $Ia$ and $Ib$.</figcaption></figure>

Finally, one can consider the scaling of the basis functions by plotting the normal basis functions corresponding to the configurations $Ia$ and $Ib$ in the lowest order case, i.e. for $k=0$ (see Figure 3.10). There, only the configuration $Ib$-using a pointwise value- scales to one. The fully moment-based configuration $Ia$ scales to another constant that depends on the edge's length and orientation with respect to the axes. This example emphasises that the configuration $Ib$ leads to basis functions which share similar properties like the Raviart – Thomas elements.

**Investigating the second configuration of degrees of freedom** For the sake of completeness, let us represent the shape of the normal component of the normal basis functions obtained with the configurations $IIa$ and $IIb$ in Figure 3.11. There, contrarily to the configuration $I$, only one function is degenerating as an internal one (in $IIb$, two quadratic functions have an identical normal component, making them not distinct in the second picture). Indeed, the focus is made on the coordinate-wise determination of the normal component. As we are in two dimensions, two functions are required to determine the coordinate-wise constants. The similarity with the configuration $I$ lies in the scaling. Indeed, only the type $b$ using pointwise value forces the scaling of the dual basis functions to one in the lowest order case.
<figure><img src="image_47.png" /><figcaption>**Figure 3.11.** Normal basis functions issued from one same edge. From left to right: ($k = 2, IIa$), ($k = 2, IIb$), ($k = 0, IIa$), ($k = 0, IIb$).</figcaption></figure>

### 3.3.8.2 Elements constructed in a reduced setting

As a last example, let us derive some results obtained for the reduced element $Ib$, offering a complete parallel with the Raviart – Thomas elements on the boundary by suppressing the further coordinate-wise liberty granted by the general setting. The internal basis functions being unchanged from the general setting, they are not represented.

<figure><img src="image_48.png" /><figcaption>**Figure 3.12.** Top left: regularity of the components of one representative of the basis functions for the reduced element $Ib$ and $k = 2$ within the element. Top right: its normal component along the boundary. Bottom: normal component of all the functions generated from the edge number 5, plot on the edge number 5. From left to right: $k = 0, k = 1, k = 2$.</figcaption></figure>

In this context, one can observe on the bottom of Figure 3.12 that there is no more degenerating normal basis functions. Therefore, all normal basis functions are acting globally to characterise the polynomial behaviour of functions of the reduced $H_k(K)$ space on the boundary. Furthermore, one can observe that the scaling of the basis functions corresponding to lowest order element, as well as the amplitude of the basis functions describing the higher order ones, make the discretisation framework reliable.
### 3.3.9 Critical assessment of the proposed approach

Motivated by defining a Flux Reconstruction scheme on general polytopes [52], we have developed a new $H(\operatorname{div})$-conforming discretisation framework that can be set up on any polytope, not necessarily convex. As a final note, let us shortly wrap up the properties and observations applicable to the constructed elements, before going through perspectives of the presented work.

**Limitations** First of all, let us mention that as for any $H(\operatorname{div})$-conforming element, the use of the developed elements from its basis functions through a reference shape can be delicate. Indeed, preserving the properties of the normal components through geometrical transformations require conforming mappings, that can be extremely delicate, especially in the particular case of non-convex polytopal geometries. As a consequence, the use of the elements as formulated by their degrees of freedom should be favoured, not requiring to express any basis function nor to use a reference element.

More specific to the configuration $\mathcal{H}$ of the developed elements, it can also be observed that the conditioning of the determination matrix is numerically sensitive towards the orientation of the element, impacting the overall robustness. Even though this drawback is easily avoidable by small local modifications of the projectors, a more comprehensive solution would be desirable.

**Strengths** On a positive note, let us point out that merging the geometrical flexibility of the Virtual elements setting with the boundary properties of the Raviart-Thomas elements, the introduced framework is vectorial and allows a lot of flexibility in the definition of the degrees of freedom. In particular, the choices of discretisation quality and degrees of freedom on the boundary are independent from the ones made within the element. The discretised quantities benefit from an extensive coordinate-wise freedom. Hence, upon the choice made while selecting the degrees of freedom, some dual normal basis functions may be reclassified into internal ones. As a consequence, allowing a complete parallel with the Raviart-Thomas setting on the boundary from the lowest order on requires to construct straightforwardly a reduced space, along with reduced elements.

Detailing a particular example of a discretisation framework through a series of spaces and the definition of particular elements, it could be observed that in both general and reduced frameworks the type of used degrees of freedom (pointwise values or moments) impacts the shape and scaling of the dual basis functions. Such behaviour could in particular be observed in the lowest order case of the given examples, where only the dual basis functions of the configuration of type "b" scale to one. This implicit access to the shape of the basis functions from the degrees of freedom's expressions therefore al-
lows us to instil the desired properties of basis functions directly at the dis-
cretisation level. The degrees of freedom expressions not being as rigid as
in the Raviart – Thomas framework, one can then develop one’s own $H(\text{div})$–
conforming element depending on further features of interest, typically im-
posing non-negative basis functions or a specific scaling.

**Perspectives** To conclude, let us point out that the results presented here are already of theoretical interest, guaranteeing the considerations about Flux Reconstruction schemes on general polytopes hold, and that the conjecture about the correction functions made in [52] is correct.

One could now consider to explore the existence and definition of projectors
from the introduced spaces onto polynomial ones, hence completing the anal-
ysis of the developed elements and widening their application range. There, a
possibility is to follow the line drawn by the Virtual Elements methods, where
projectors have been constructed in [100] and seems to be generalisable to the
presented developments. Then, the projectors defined, estimates could be de-
rived and the introduced discretisation spaces and elements used in a more
practical context, as by example in a finite element framework.
<figure><img src="image_49.png" /></figure>
A conservative and continuity
preserving correction procedure
for Level Set methods in evolv-
ing interfaces

Another facet of conservation problems lies in the scope of two-phase flow modelling, field investigating the motions of two fluids that coexist in a same domain. There, when the dynamic of each fluid is driven by a set of conservative hyperbolic equations, the mass quantity of each phase is preserved through time. Nevertheless, the spatial distribution of each phase is possibly going through heavy distortions and topological changes.

An accurate approximation of the solution to those problems is therefore
delicate, as it requires a time-dependent spatial representation of the inter-
face that does not prevent the comprehensive numerical scheme to be mass
conservative. In particular, at the discrete level, the local conservation of the
partial masses usually contradicts with a sharp resolution of the interface.
As a consequence, the few numerical methods achieving both usually have to
waive interface continuity, yielding another major issue.

Focusing on numerical schemes based on unstructured meshes, we present here a corrective numerical approach that enforces the mass conservation at the discrete level while preserving the cross-cell continuity of the interface. After specifying the considered problem in Section 4.1, we present in Section 4.2 the three fundamental works from which stem our investigations; two classical approaches and a first conservative method combining the strengths of each classical approach together. We then improve this first bridge by deriving a novel conservative and interface-continuous approach in Section 4.3.
## Chapter contents

<table><tr><td>4.1</td><td>Problem statement, modelling, and approximation strategy</td><td>169</td></tr><tr><td>4.1.1</td><td>Modelling strategy</td><td>170</td></tr><tr><td>4.1.2</td><td>Conservation principle</td><td>174</td></tr><tr><td>4.1.3</td><td>Dedicated approaches and numerical schemes</td><td>177</td></tr><tr><td>4.1.4</td><td>Scope of the presented approach</td><td>183</td></tr><tr><td>4.2</td><td>Description of three fundamental methods</td><td>184</td></tr><tr><td>4.2.1</td><td>The Level Set approach</td><td>184</td></tr><tr><td>4.2.2</td><td>The Volume of Fluid method</td><td>186</td></tr><tr><td>4.2.3</td><td>The THINC-LS approach</td><td>188</td></tr><tr><td>4.3</td><td>A continuity preserving procedure</td><td>192</td></tr><tr><td>4.3.1</td><td>The general idea in a nutshell</td><td>192</td></tr><tr><td>4.3.2</td><td>An intuitive continuous interface correction</td><td>193</td></tr><tr><td>4.3.3</td><td>Extensions</td><td>205</td></tr><tr><td>4.3.4</td><td>Numerical approximations</td><td>219</td></tr><tr><td>4.3.5</td><td>Numerical results</td><td>230</td></tr><tr><td>4.3.6</td><td>Critical assessment of the proposed approach</td><td>242</td></tr></table>
## 4.1 Problem statement, modelling, and numerical approximation strategy

Two-phase flow modelling investigates the time-space behaviour of two different entities coexisting in a same domain $\Omega$, called phases. Considered in their largest definition, those phases are characterised by the entities' nature, that is, the combination of its chemical composition (e.g. oil, water) and physical state (liquid, gaseous, solid). By example, with this definition, aqueous vapour-water, oil-water and ice-oil interaction are all considered as two-phase problems. Each of the phases is further described by a set of variables $Y$ representing its dynamical state, usually related to density, pressure and velocity.

The set $Y$ is evolved accordingly to a dynamic that may be specific to each phase and that can be impacted by external conditions imposed on the domain as e.g. imposed velocity, pressure or heat flux. Moreover, the two phases may interact themselves through their interface $\Lambda$, a time-dependent hypersurface of the domain $\Omega \subset \mathbb{R}^d, d \in \mathbb{N}$, where both phases are spatially in contact. As the two phases split the domain $\Omega$, one usually considers only one reference phase lying on a subset $\Omega_0 \subset \Omega$, inferring the spatial description of the second phase by the complementary relation $\Omega_1 = \Omega \setminus \Omega_0$.

<figure><img src="image_50.png" /><figcaption>**Figure 4.1.** Representation of the considered domain, phases and interface, the interface not including the domain's boundaries.</figcaption></figure>

Approximating the time-space behaviour of the two phases simply reduces to retrieving both the dynamical state and the location of the interface, the latter dictating by itself the spatial distribution of the two phases. However, as the interface location is a direct consequence of the fluid motion, its shape can become tremendously complicated and even change topology across time as soon as the velocity is not spatially constant. In some applications, it is possible to know *a priori* if those features will develop depending on the domain shape and the intrinsic properties of each phase. Typically, the Reynolds
number will indicate if the fluid is prone to turbulence, and the interplay be-
tween the gas mass flow and the interfacial liquid velocity can foresee the flow
regime and therefore the interface's topology type. A typical classification for
horizontal pipes is given as an example in Figure 4.2.

<figure><img src="image_51.png" /><figcaption>**Figure 4.2.** Prominent flow behaviours encountered in horizontal pipes. (a)-(h): Different flow regimes. (1)-(3): Different flow natures identified when pulsing in the second phase into the reference one.</figcaption></figure>

Depending on its spatial distribution, each phase is either called dispersed or continuous. The type of spatial cohabitation of the two phases within the considered domain give rise to different modelling strategies.

**4.1.1 Modelling strategy**

Ideally, the physical model should describe accurately the phases properties
and dynamic, grasp the topologically changing location of the subdomains and
interface, and consider quantity exchange between phases as heat and mass
transfer, if and when relevant.

In practice, considering all of those
properties at once requires a thorough
physical model that is delicate to discre-
tise, when even possible. Therefore, it
is natural to select directly a modelling
approach that is tailored to the simu-
lation's finality and the foreseen flow
regime. Indeed, if for example simulat-
ing the behaviour of breaking waves or
annular dynamics requires an accurate
spatial description of the interface, the
motion prediction of a dispersed pollu-
tion cloud does not require an accurate
boundary representation of each small
pollution particle.

<figure><img src="image_52.png" /></figure>
At the discrete level, the representation of continuous phases is classically done in the Eulerian framework by fixing a reference frame and observing the phases passing through so-called control volumes, small portions $C$ of the domain $\Omega$ where the state quantities are discretised, monitored and updated. Disperse phases may also be immersed in this framework when the spatial distribution of the interface or its properties are important to the dynamic (see Figure 4.3a). In contrast, when the individual particles dynamic prevails and that the surface effects and interface deformations are negligible it may be suitable to use a Lagrangian framework [113, 114], where the dispersed entities are reduced to points to follow in time (see Figure 4.3b). An intermediate perspective called interface tracking [115, 116, 117], mainly suited for cases where the interface distortions are limited, consists in discretising both phases in an Eulerian framework while explicitly keeping track of the interface by following its Lagrangian markers (see Figure 4.3c).

<figure><img src="image_53.png" /><figcaption>**Figure 4.3.** Depiction of classical two-phase flows modelling strategies.</figcaption></figure>

**Model and assumptions** As we aim here to spatially describe a highly distorting interface that separates two phases whose connected components are of non-negligible size, we consider a fully Eulerian framework.

Within this scope, and to ease the presentation of the corrective technique, we assume that the two phases can be described identically in both subdomains by a time-space dependent state vector $Y: \mathbb{R}_+ \times \Omega \to \mathbb{R}^q, q \in \mathbb{N}$, usually related to the pressure $p$, velocity $v$, and internal energy $e$. We further assume that the phases dynamic is driven by the same set of hyperbolic equations (4.1), referred to as the system dynamics.

$$ egin{cases} \frac{\partial Y}{\partial t} + \nabla \cdot \mathcal{F}(Y) = 0 \\ Y|_{\partial \Omega} = g, Y(0, \cdot) \in \mathcal{D}, \text{ given.} \end{cases} \qquad (4.1) $$

There, the flux $\mathcal{F}: \mathcal{D} \to \mathbb{R}^q$ describes the considered physical model, as e.g. the Euler or the shallow water equations. The boundary conditions, here simply exemplified as Dirichlet with for second member $g: \mathbb{R}_+ \times \partial\Omega \to \mathbb{R}^q$, specify the impact of the domain's boundary on the dynamic. This set of hyperbolic equations is then closed by the equation of state (4.2), describing the relationship between the state variables accordingly to phase-wise properties.

$$ EoS(Y, \psi) = 0 \qquad (4.2) $$

Those phase-wise properties are selected through the time-dependent value of $\psi$, the interface representation defined by the Lagrangian marker (4.3) that
provides a fluid selector. This representation is initiated by a function $\psi_0: \Omega \to \mathbb{R}$ associated to the initial subdomains and evolved either by (4.4a) or (4.4b) accordingly to the fluid motion, driven by the underlying velocity field $v$.

$$ \psi: \mathbb{R}_+ \times \Omega \to \mathbb{R} \quad (4.3) \quad \left\{ \begin{array}{l} \displaystyle \frac{\partial \psi}{\partial t} + v \cdot \nabla \psi = 0 \quad (4.4a) \\ \displaystyle \psi(0, \cdot) = \psi_0 \end{array} \right. \quad \left\{ \begin{array}{l} \displaystyle \frac{\partial \psi}{\partial t} + \nabla (v \psi) = 0 \quad (4.4b) \\ \displaystyle \psi(0, \cdot) = \psi_0 \end{array} \right. $$

The two-phase state and spatial dynamics is thus completely described in the Eulerian–Eulerian framework by the comprehensive set of equations (4.1–4.4), defining a so-named *one-fluid model* in the literature.

Yet, the initial value problem (4.4a) does not consider the nature of the phases nor their properties. On the same line, the system dynamics (4.1) does not consider any interfacial effect, as no $\psi$-dependent source term representing a possible heat exchange or phase interaction is involved. The interface is thus considered as passive, simply transported as a free surface by the underlying velocity field and making in particular the phase selection of importance only in the equation (4.2).

**Dynamics of the discretised model** When considered continuously, all the above-mentioned dynamics are interlinked and evolve simultaneously. At the discrete level however, each of the equations is discretised in both space and time, making the interplay effective only at specific points in the time space. A partially sequential approximation of the dynamic has therefore to be introduced, as depicted in Figure 4.4.

<figure><img src="image_54.png" /><figcaption>Figure 4.4. Semisequential update of the time-discretised model.</figcaption></figure>
At each time step, the space dependent phase properties are inferred from the interface representation accordingly to the value of the fluid selector. The equation (4.2) thence closes the system (4.1), and the state is evolved to the next time step. Simultaneously, the interface location is updated from the fluid selector values accordingly to the motion of the interface representation (4.4a), motion induced by the velocity field that is directly obtained from the initial state vector. Due to this step-wise information exchange, a reliable numerical solution can only be obtained by designing a numerical scheme that guarantees the discrete state conservation principles while providing an interface location that is in agreement with the mass quantity of each phase.

**Discretisation issues** Driven by conservation laws where the total mixture's mass is preserved despite a possibly reckless velocity field acting on the phases (see Section 4.1.2), the solutions to the problem (4.1–4.4a) may experience heavy distortions and topological changes of their subdomains and interfaces repeatedly in time (see Figure 4.5 for an illustration).

<figure><img src="image_55.png" /><figcaption>Figure 4.5. Topological changes that the interface may undergo.</figcaption></figure>

The discrete spatial representation of the interface has therefore to be topologically flexible, robust, and accurate while preserving the continuity across control volumes. Indeed, as the interface location has a direct impact on the state equation and thus on the upcoming fluid motion, failing in either of those properties will notably yield topological concerns, spurious deformations, and more dramatically front speeds and phase locations inaccuracies fuelling conservation issues.
### 4.1.2 Conservation principle

In order to develop a correction procedure that turns any comprehensive numerical scheme into a mass-conservative one that is free from the above numerical issues, let us first derive the discrete mass-conservation principle. More than specifying the scheme property we aim to enforce, it will also allow us to identify the reason of extra numerical difficulties and trade-offs.

**Continuous two-phase mass conservation principle** The time dynamic of a flow is in particular ruled by the mass conservation principle acting on the full domain $\Omega$ and directly expressed in the definition of $Y$ and $\mathcal{F}$ whenever a system of conservation laws is investigated. We saw in the preliminary Section 2.1.2.1 that this conservation principle reads in its integral form

$$ \int_{\Omega} \frac{\partial \rho}{\partial t} dx + \int_{\Omega} \nabla \cdot (\rho v) dx = 0, \qquad (4.5) $$

where $\rho$ represents the quantity responsible for the mixture's mass (e.g. the fluid's density or height) and the term $\int_{\Omega} \nabla \cdot (\rho v) dx$ expresses the mass transfer, here reducing to a possibly net mass income from the boundary.

In the particular case of immiscible two-phase flows, where no mixing nor chemical transformation occur, this principle also holds for each of the phases themselves. Each of the phase-wise total masses

$$ M_i: \mathbb{R}_+ \rightarrow \mathbb{R} \\ t \mapsto \int_{\Omega_i(t)} \rho \, dx = \int_{\Omega} \rho \, 1_{\Omega_i(t)} \, dx, \quad i \in [0, 1], \qquad (4.6) $$

is there constant through time and satisfies the relation

$$ \frac{dM_i}{dt} = \int_{\Omega} \frac{\partial(\rho 1_{\Omega_i(t)})}{\partial t} dx + \int_{\Omega} \nabla \cdot (\rho v 1_{\Omega_i(t)}) dx. \qquad (4.7) $$

At the continuous level, those phase-wise relations are directly enforced by (4.5) due to the gradient linearity and the disjoint union relation $\Omega = \Omega_0 \cup \Omega_1$, as illustrated in Figure 4.6. Indeed, the covering of $\Omega$ by the two phases yield $\int_{\Omega} \rho dx = \int_{\Omega} \rho 1_{\Omega_0(t)} dx + \int_{\Omega} \rho 1_{\Omega_1(t)} dx$ for any $t \in \mathbb{R}_+$, splitting (4.5) accordingly.

<figure><img src="image_56.png" /><figcaption>**Figure 4.6.** Illustration of the continuous conservation principle.</figcaption></figure>
**Spatially discrete two-phase mass conservation principle** At the spatially discrete level, the same mass conservation principle applies directly on any time-independent control volume $C \subset \Omega$, where $\Omega$ represents the computational domain approximating $\Omega$. There, the change in the partial mixture's mass $m: t \mapsto \int_C \rho dx$ is simply given by

$$ \frac{dm}{dt} = \int_C \frac{\partial \rho}{\partial t} dx + \int_C \nabla \cdot (\rho v) dx, \qquad (4.8) $$

and the term $\int_C \nabla \cdot (\rho v) dx$ represents the quantity of mass received from and sent to neighbouring control volumes (see Figure 4.6). This transfer can be seen more easily by applying the Stokes theorem to the relation (4.8), yielding

$$ \frac{dm}{dt} = \int_C \frac{\partial \rho}{\partial t} dx + \int_{\partial C} \rho v \cdot n dx, \qquad (4.9) $$

expression emphasising each transfer across the control volume's boundary. In addition, as in the continuous case when considering immiscible phases, each of the phase-wise partial masses defined over any control volume, i.e.

$$ m_{i,C}: \mathbb{R}_+ \to \mathbb{R} \\ t \mapsto \int_{C \cap \Omega_i(t)}^{\rho} x \, d\mathbf{x} = \int_C^{\rho \mathbb{1}_{\Omega_i(t)}} x \, d\mathbf{x}, \quad i \in [0, 1], \qquad (4.10) $$

is conserved through time and follows the spatially discrete principle

$$ \frac{dm_{i,C}}{dt} = \int_C \frac{\partial(\rho \mathbb{1}_{\Omega_i(t)})}{\partial t} dx + \int_C \nabla \cdot (\rho v \mathbb{1}_{\Omega_i(t)}) dx, \quad i \in [0, 1]. \qquad (4.11) $$

Considering then a set of control volumes $T^h = \{C\}$ that covers $\Omega$, the sum of all the control volume-wise conservation relations

$$ \sum_{C \in T^h} \frac{dm_{i,C}}{dt} = \sum_{C \in T^h} \int_C \frac{\partial \rho \mathbb{1}_{\Omega_i(t)}}{\partial t} dx + \sum_{C \in T^h} \int_{\partial C \cap \Omega_i(t)} \rho v \cdot n dx \qquad (4.12) $$

reduces to (4.7), and the reconstructed phase-wise masses $M_i, i \in [0, 1]$ are preserved (see Figure 4.7 for an illustration). However, this reconstruction is only possible provided an exact representation of the interface in the control volumes, whose resolving is extremely delicate numerically.

<figure><img src="image_57.png" /><figcaption>**Figure 4.7.** Illustration of the discrete conservation principle.</figcaption></figure>
Indeed, any finite dimensional discretisation space defined over a control volume does not necessarily allow an exact representation of the interface

location within the cell, furthermore continuously in time. Hence, the approximated interface location does not necessarily preserve the distribution of both phases within the control volume accordingly to the mass conservation principle, neither in the phases' proportions nor in their spatial layout. It is thus very challenging to retrieve a numerical solution for which (4.12) reduces to (4.7), as the definition of

the dedicated numerical schemes outlined in Section 4.1.3.1 shows.

<figure><img src="image_58.png" /><figcaption>Figure 4.8. Interface approximation within a control volume.</figcaption></figure>

**Fully discrete two-phase mass conservation principle** Discretising further in time, the conservation statement (4.11) has to be understood with respect to a step-wise time evolution of the solution. We therefore consider the cumulative mass change, dictated in any control volume $C$ over each complete time step $\Delta t = t^{n+1} - t^n$ by the time-space discrete conservation relations

$$ \int_C [\rho \mathbb{1}_{\Omega_i(t)}]_{t^n}^{t^{n+1}} dx + \int_{t^n}^{t^{n+1}} \int_C \rho \mathbb{1}_{\Omega_i(t)} v \cdot n dx dt = 0 \quad \forall C \subset \Omega, \forall i \in [0, 1]. \quad (4.13) $$

When applied to a numerical solution, the left hand side value is representative of the scheme's conservativeness as it quantifies the cumulative mass loss, indicating the scheme defects even when the errors due to the interface location are compensating within the considered time step.

However, the numerical monitoring of (4.13) is itself subject to the approximation of both time and space integrals. Thus, on top of the numerical difficulties encountered when designing the numerical scheme, quadrature defects will impact the mass loss quantification and thus the validation of the scheme. In particular, the point distribution of the space quadrature

may not suit the fluid profile, as illustrated in Figure 4.9. Furthermore, the accuracy of the time quadrature is bound to the number of possibly subtimesteps and their distribution in the time interval. Thus, the choice of the discretisation space in which the fluid selector is expressed also impacts how well the scheme's performance can be assessed, making its design even harder.

<figure><img src="image_59.png" /><figcaption>Figure 4.9. Impact of the quadrature points on the evaluated mass.</figcaption></figure>
**Note** When considering incompressible flows, not only the mass but also the phases' volumes are preserved. Thus, instead of monitoring (4.6), one can simply consider the values of $|\Omega_0|$ and $|\Omega_1|$, and instead of monitoring (4.13) one can simply check that $\int_{t^n}^{t^{n+1}} \int_{\partial C \cap \Omega_i(t)} v \cdot n \, dx \, dt = 0$. ▲

### 4.1.3 Dedicated approaches and numerical schemes

Due to the high complexity of two-phase flow problems that manifest both in their behaviour and modelling up to the discrete level, the development of dedicated numerical schemes is relatively recent, especially in the case of two-phase flows driven by hyperbolic systems of equations. Indeed, hitherto unresolved discontinuities, stability and conservation issues disallowing any trust in the solution were already encountered in the single-phase case, field remaining an active research area in itself.

In order to motivate the choice of our starting points detailed in Section 4.2, let us outline the main numerical achievements in tackling two-phase flow problems that are modelled by (4.1–4.4) in the Eulerian–Eulerian framework. Within this scope, all the mentioned methods are of *interface capturing* type, a class of schemes where the interface is inferred from the flow state without being specifically tracked. There, as observed when deriving the mass conservation principle in Section 4.1.2, the discrete preservation of the phases’ mass is often in contradiction with an accurate interface representation.

Most of the interface capturing approaches therefore only focus either on preserving mass conservation or on resolving accurately a sharp interface location. Hybrid approaches attempting to resolve both are newfangled and still suffer from substantial drawbacks. Though not being able to be exhaustive, we aim to provide a progressive view of resolved issues and current approximation challenges in those three directions.

#### 4.1.3.1 Approaches focusing on mass conservation

Approaches focusing on mass conservation usually see the domain $\Omega$ as a testellation of small fluid portions on which physical quantities are investigated. Initially built as thorough physical models treating the interface as a thick interaction layer, they have been then simplified to the case of thin passive interfaces, paving timidly the way for a sharp resolution of their location.

**Diffuse interface** One of the pioneering works focusing on mass conservation in two-phase flows is the diffuse interface model introduced by Anderson, McFadden and Wheeler [118], inspired by the physical theory of phase-transitions [119, 120]. In this approach, the fluid selector $\psi$ is defined by a continuous phase-representative parameter constructed from averaged molec-
ular quantities on small fluid portions, exploiting typically the molecular density field. Considered as the energy-free surface of the fluid selector, the interface is then advected as to minimise the total energy coming from elementary intermolecular interactions while preserving the equilibrium profile across itself [121] by observing the energy-free surface rule [122]. This dynamic is usually enforced by augmenting (4.4b) with the Cahn-Hilliard equation [123].

This model is physically accurate as long as the size of each fluid portion crossed by the interface matches the region where microscopic physical interactions occur. The interface has thus a physically meaningful width, giving its name to the diffuse interface model. Besides, the conservation is guaranteed whenever a consistent numerical scheme is used in the model's discretisation. Dramatical numerical behaviours are however encountered in the cases of very thin interfaces and topological changes, due to emerging singularities.

**Phase field** Often assimilated with the diffuse interface in the literature, the phase-field method has been introduced quasi-simultaneously by Jacqmin in [124] and differs only in the way the interface is represented. Yet, as an emergence from the modelling of free boundary problems [125], the method is restricted to flows whose tangential stress across the interface can be neglected, i.e. to the free surface case. There, the molecular interactions across the interface and their width are negligible, and the physical meaning of the interface can be consequently dropped. The interface is thus considered as mathematically sharp and its representation reduced to a phase field, an ad-hoc indicator function $\psi$ satisfying $\psi = 1$ in the reference phase and $\psi = -1$ otherwise. The transition between the two states has then to be artificially smoothened to allow numerical computations. As the function's shape does not impact the model dynamics, it is usually taken as

$$ \begin{gathered} \psi: \Omega \rightarrow [-1, 1] \\ x \mapsto -\tanh(\beta x), \end{gathered} \qquad (4.14) $$

where $\beta$ acts as a sharpness control parameter and determines the artificial interface width, adapted upon the numerical schemes needs.

Unlike the diffuse interface method, the interface representation is not bound to any physical quantities and therefore handles natively topological changes, making it popular for modelling turbulent flows [126], binary alloys [127] and physical phase change [128]. However, going towards the sharp interface limit by reducing the artificial width of the interface still yields instabilities. Furthermore, another severe drawback relates to conservation, where regardless of the employed numerical scheme the partial phase-wise masses are not necessarily conserved. Indeed, due to the minimisation driven coupled convection-advection dynamic of $\psi$, any perturbed energy equilibrium of the interface is restored by artificially diffusing into the other phase, yield-
ing partial mass leakages that shrinks small domains and grows already large ones. As a consequence, the discrete mass conservation is only achieved at the global level (4.8), not down to the control-volume wise level in the sense of (4.11). To reduce the impact of this issue, profile corrected [129] and flux corrected [130] approaches have been investigated, without succeeding yet to provide a fully conservative framework.

**Free surface capturing** Another approach tackling two-phase flows separated by a passive interface is the free surface capturing method introduced by Kelecy, Pletcher [131], Pan and Chang [132], here again quasi-simultaneously with the diffuse interface method. It also adopts its perspective, but instead of augmenting the equation (4.4b) to incorporate the interface dynamic, the free interface is directly described in the system (4.1–4.4b) as a contact discontinuity in the density field. The conservation is there immediate by definition of the enforced dynamics, provided that a consistent numerical scheme is used.

So defined, the free surface capturing method can handle complex surface representation, deformations, and to some extent topology changes. Very popular in modelling dam-break and tsunami problems, its numerical discretisation uses however high resolution Riemann solvers coupled with the artificial compressibility method, being thus very demanding for non-linear problems.

**Volume of Fluid** A very popular strategy among approaches that focus on mass conservation is to consider the interface representation $\psi$ cell-wise as a volume fraction $\bar{H} \in [0, 1]$ representing the proportion of the reference phase within the mixture contained in a given control volume. The motion of the interface is then simply described as the advection of its representation by the underlying fluid velocity field (4.4b). Introduced by Hirt and Nichols in [133] under the name *Volume of Fluid*, this modelling strategy is flexible with respect to the mesh geometry and is naturally mass conservative, provided that a consistent numerical scheme is used in the approximation of (4.1–4.4b).

This numerical approximation has been investigated from several perspectives. A first direction has been to use Riemann solvers [134], linearized shortly after in [135]. In the same line, Godunov type schemes were considered in [136, 137] and later extended to (W)ENO schemes, still investigated nowadays [138, 139]. Parallelly, the development of flux vector [140, 141] and flux difference [142, 143, 144] splittings schemes yielded accurate simulations and have been extended to higher dimensions [145]. Further efforts to reduce oscillations have been made in [146], and positivity preservation has been achieved in [147]. Later, hybrid flux splitting methods have been developed [148, 149, 150, 151] in order to reduce the computational cost and restrictions caused by the need of the Jacobian, and the need of eigenstructure has been removed in [152].
Other Finite Volume methods also aimed at retrieving an interface location that is consistent with the physics for the turbulent case [153, 154, 155], paying an extra attention to the flux discretisation. More recently, Discontinuous Galerkin-like collocation methods have been investigated in [156], paving the way to entropy stability [18, 157, 158].

Parallelly, to offset the fact that the Volume of Fluid method provides the interface location only as a layer of interfacial cells, reconstruction techniques have been developed to infer it more precisely within each interfacial control volume. A first method, *SLIC* [159], defined the interface as a line parallel to one of the control volume's boundaries, and has been improved by the *PLIC* reconstruction [160] where a linear interface is determined with respect to the gradient field of the volume fraction distribution. In the same vein, spline and quadratic representations have been later investigated in [161] and [162]. More involved reconstruction techniques were also developed, as e.g. the Moment Of Fluid method that minimises the phase-wise volume defects in each control volume [163], or the physics-aware *OST* [164] and *PROST* [165] methods, taking into account the surface tension. However, all of those techniques appear to be unsatisfactory as they do not ensure the continuity of the interface across cells nor an actual physical relevance of their shape and location.

### 4.1.3.2 Approaches focusing on the interface location

Resolving accurately the interface location is usually achieved by carefully designing the interface representation so that it can be precisely followed or tracked in time. The most accurate methods in this direction are mainly of Lagrangian type, as e.g. the Sharp Interface [166, 167], Ghost Fluid [168], Immerse Boundary [169] and more recently the cut cells [170, 171] methods, falling therefore out of our scope.

Among the few interface capturing methods that focus on resolving precisely the interface location, accurate techniques making use of Normalised Variable Diagram [172] and Flux Corrected Transport [173] were developed. Providing a sharp description of the interface down to discrete level and making a step towards conservation, those methods however encounter untoward stability problems. Less accurate but more reliable, the widely used Level Set method [174, 175, 176] dedicates itself to the passive interface case. There, the interface is defined implicitly by the zero level of a $d+1$ dimensional surface that is evolved in time via its advection by the underlying fluid velocity field, obeying the equation (4.4a). Not requiring any reconstruction and being suited for representing any interface regardless its topology, this method does not however control the interface shape by any dedicated physical constraint and is far from being conservative.
Another perspective in improving the interface resolution is the use of
mesh-adaptation techniques [177, 178, 179, 180] and hybrid-meshes [181],
where the interface is sharpened through mesh refinement rather than by
elaborating suitable discretisation spaces. Despite allowing the application
of specific local treatments, the width of the interface then usually remains
linked with the interfacial cell widths, and a fine resolution can only be ob-
tained at great computational cost.

4.1.3.3 Hybrid methods

Aiming to retrieve conservation while specifying the interface location, hybrid methods borrowing the strengths of several schemes that belong to either of the above classes have been investigated.

Using coupling techniques, those methods usually adjust the definition of
numerical fluxes in order to enforce the conservativeness of the interface rep-
resentation, exploiting natively conserved quantities. A natural approach in
this direction is to combine the Level Set and Volume of Fluid methods, ini-
tiated in [182] and popularised with the CLSVOF method [183] by making
use of operator splitting techniques. The need of interface reconstruction has
been later removed in [184], and flexibility with respect to the mesh type has
been improved in [185] and [186]. Likewise, a coupled Level Set and Moment
of Fluid method, CLSMOF [187], has also been proposed. All of those methods
are by definition fully conservative and provide a precise location of the inter-
face. However, as modifying the numerical fluxes is subject to discretisation
and geometrical restrictions, the topological flexibility of the Level Set method
does not completely transfer to the coupled methods, and issues in highly dis-
torted control volumes and topological changes are typically encountered.

To overcome those restrictions, more recent corrective approaches seem
promising. The location of the interface is there corrected after its standalone
evolution, within each control volume and accordingly to local conservation
constraints. The coupling is therefore not simultaneous in discrete time, and
no geometrical restriction is entailed by the correction process. In practice,
works defining such procedures originally made use of projection methods
[188] or predictive-corrective approaches [189]. More recently, a promising
procedure coupling the Level Set and Volume of Fluid methods by using an
intermediate interface reconstruction inspired from the phase-field approach
has been proposed in [190] and developed in [20]. Both conservative and geo-
metrically flexible, there is however no cross-cell interface continuity and no
physical relevance of the corrected interface location.
4.1.4 Scope of the presented approach

In this chapter, we propose to construct a continuous and conservative repre-
sentation of the interface that separates two flows whose motion is ruled by
the system (4.1–4.4). As in this model the interface is assumed to be passive
and is possibly subject to large deformations, we consider a fully Eulerian
framework and focus on the particular class of interface capturing schemes.

**Background choices and starting point**

Observed while deriving the discrete conservation principle in Section 4.1.2 and noticeable in the historical developments of numerical schemes, shortly reviewed in Section 4.1.3, continuous and conservative representations of the interface are not natively compatible. Thorough representations are therefore not retrieved by stand-alone models, but rather by hybrid approaches that bring the strengths of separate monotasked numerical schemes together.

We aim to step further in this direction by improving existing hybrid ap-
proaches. To allow an independent numerical discretisation of each scheme
component and to prevent any topological issue in the interface representa-
tion, we focus on corrective approaches that do not restrict the use or features
of any scheme component. Each of those components is then selected upon
the properties we aim to instil in the comprehensive hybrid scheme. Here, we
consider a level set representation for determining the location of a possibly
topologically changing interface, and select the Volume of Fluid approach for
providing a set of phase-wise discrete mass conservation constraints. Given
those three choices, the hybrid method THINC-LS [20] appears as a natural,
flexible, and powerful framework. We therefore consider it as our starting
point, combining the Level Set and Volume of Fluid approaches in a correc-
tive way so that the discrete conservation principle (4.11) holds.

**Objectives and achievements** The two main drawbacks of the *THINC-LS* approach are the lack of cross-cell continuity and physical meaning of the retrieved interface. We eliminate here the former by redefining the discrete representation of the initial interface and modifying the expression of the *THINC* correction term. Fully conservative, this technique can be applied to any numerical scheme assessing (4.1–4.4a), provided that a pointwise reconstruction of the solution is available.

**Outline** In order to ease the presentation of the proposed technique, we start in Section 4.2 by describing more in details the three fundamental works from which stem our investigations. We then describe the core of our improved method in Section 4.3, presenting an intuitive design of an admissible correction function in Section 4.3.2. Various extensions are then given in Section 4.3.3, followed in Section 4.3.5 by the numerical details required in practical applications. Essential numerical results are finally showed in Section 4.3.5.
## 4.2 Description of three fundamental methods

As the presented technique consists in an improvement of the hybrid THINC-LS approach, we first describe the essentials of this approach. We start by describing each of the involved components, that are the Level Set and the Volume of Fluid methods, presented respectively in Section 4.2.1 and 4.2.2. The classical THINC-LS approach is then introduced in Section 4.2.3.

### 4.2.1 The Level Set approach

Introduced in particular by Dervieux, Thomasset, Osher, Sethian, Sussmann and Fedkiw [174, 175, 176], the Level Set approach is an interface-capturing method that relies on an implicit description of the interface location. Providing a continuous phase-selector to the system dynamics, this interface representation allows a pointwise selection of the phase properties to consider. However, despite allowing an easy spatial representation of the interface, this method is prone to conservation errors.

**Implicit geometry representation** In the Level Set approach, the interface is represented implicitly as the zero level of a $d+1$ dimensional evolving surface, parametrised by a function $\phi: \mathbb{R}_+ \times \Omega \to \mathbb{R}$. The space-time-dependent sign of the function selects the phase, a positive sign representing a given reference phase and a negative one its complementary in $\Omega$. Hence, assuming without loss of generality that the reference phase is occupying the subset $\Omega_0$, the level set function $\phi$ satisfies the following property.

$$ \phi: \mathbb{R}_+ \times \Omega \to \mathbb{R} \\ (t,x) \mapsto \begin{cases} \phi(t,x) > 0 & \text{if } x \in \Omega_0(t) \\ \phi(t,x) = 0 & \text{if } x \in \Lambda(t) \\ \phi(t,x) < 0 & \text{if } x \in \Omega_1(t) \end{cases} \quad (4.15) $$

Even though the shape of the function $\phi$ is a-priori irrelevant, one usually seeks $||\nabla\phi||=1$. One thus generally defines the initial level set $\phi_0: \Omega \to \mathbb{R}$ as the signed Euclidean distance function to the interface, that is,

$$ \forall x \in \Omega, \quad \phi(0,x) = \phi_0(x) := \min_{y \in \Lambda(0)} ||x-y||_2. \quad (4.16) $$

<figure><img src="image_60.png" /><figcaption>**Figure 4.10.** Topological flexibility in representing interfaces by a level set.</figcaption></figure>

Note that only one level set is defined over the entire domain $\Omega$, even if the interface is made of several components. The level set function is therefore led to represent interfaces having various topologies and being possibly highly distorted. Seemly challenging, this flexibility is in fact handled naturally. Indeed, any implicit representation always relies on the same parametrisation nature of the level set function, depending only on its given expression.
**Spatial evolution of the interface** The interface location is evolved by advecting its level set representation accordingly to the velocity field $v$ derived from the state variables. Likewise, the state evolution depends on the time-dependent interface location. Therefore, the comprehensive dynamic is given by augmenting the system (4.1–4.2) by

$$ \frac{\partial \phi}{\partial t} + v \cdot \nabla \phi = 0, \qquad (4.17) $$

coming down to specifying the expressions of (4.3) and (4.4a) of the general picture given in Section 4.1.1. There, $\phi$ indicates to (4.1)-(4.2) the phase to consider at any point of the domain $\Omega$, selecting in particular the relevant equation of state and physical properties. Especially, when considering any point-value based discretisation it allows to select possibly different phases for the discretisation points of a same interfacial cell.

**Method strengths and drawbacks** By definition, the Level Set method provides an interface representation that is suitable to any time-dependent topology and allows a pointwise selection of the relevant phase over the domain $\Omega$. However, the level set function suffers from the development of singularities as its gradient becomes far from $|\nabla \phi|=1$ through its advection, even for smooth initial conditions (see Figure 4.11 for an illustration).

<figure><img src="image_61.png" /><figcaption>Figure 4.11. Development of singularities.</figcaption></figure>

As particular consequence, the interface representation may become sloppy when $\phi$ becomes too flat near the zero level set. On the computational side, spurious interfaces may even develop due to singularities or possible numerical diffusion (see Figure 4.12 for an illustration). To prevent those pitfalls, a redistancing technique can be regularly employed. Commonly consisting in evolving the Eikonal equation

$$ \left\{ \begin{array}{l} \frac{\partial \tilde{\phi}}{\partial \tau} + \operatorname{sign}(\tilde{\phi}_0)(||\nabla \tilde{\phi}||_2 - 1) = 0 \\ \tilde{\phi}(0, \cdot) = \phi(t, \cdot) \end{array} \right. \qquad (4.18) $$

until a space-dependent pseudo-final time $\tau_{F,x}$ for which $|\nabla\phi(\tau_{F,x}, x)|=1$, the level set is then reinitialized to $\phi(t,x)=\tilde{\phi}(\tau_F(x), x)$. The interface therefore stays sharp all through its evolution [191].
An other main issue is yielded by the sole discretization of (4.17). Indeed, as the level set field is practically known from the degrees of freedom only, the shape of the interface within each cell depends on the point-value reconstruc- tion technique associated to the discretisation setting. Therefore, the mass conservation principle (4.7) does not natively transfer to the discrete level.

<figure><img src="image_62.png" /><figcaption>Figure 4.12. Phase determination issue for flattened or perturbed level set.</figcaption></figure>

4.2.2 The Volume of Fluid method

Primarily focusing on mass conservation, the Volume of Fluid method tackles two-phase problems without describing specifically the interface location. In- stead of considering two fluids as distinct entities, it appraises the properties of a single mixed phase characterised within each cell by a volume fraction.

**Phase representation** This mixture is continuously defined over Ω by an indicator function H featuring a given reference phase, assumed here without loss of generality to be lying on Ω₀ ⊂ Ω, i.e.

$$
\begin{equation}
\begin{aligned}
H: \Omega &\mapsto \{0, 1\} \\
x &\mapsto \mathbb{1}_{x \in \Omega_0}.
\end{aligned}
\tag{4.19}
\end{equation}
$$

The nature of the mixture is then given cell-
wise accordingly to the average value of H on
the considered cell Cᵢ, referred to as volume of
fluid, denoted H̄ᵢ and defined by

$$
\bar{H}_i := \int_{C_i} H(x) \, dx = |\Omega_0 \cap C_i|. \quad (4.20)
$$

The associated field is then set as the piece-
wise constant function $\bar{H}$ expressed as

$$
\begin{equation}
\begin{aligned}
\bar{H} : \Omega &\mapsto [0, 1] \\
x &\mapsto \bar{H}_i, \quad \text{i s.t. } x \in C_i,
\end{aligned}
\tag{4.21}
\end{equation}
$$

and the interface is seen as the layer of cells Cᵢ
whose associated value Ĥᵢ are neither 0 nor 1.

<figure><img src="image_63.png" /><figcaption>Figure 4.13. Mixture (top) and interface (bottom) definition for given layouts of two original phases and control volumes.</figcaption></figure>

<figure><img src="image_64.png" /></figure>

The space-dependent mixture's properties to be considered in the system dynamics are then derived cell-wise from the volume of fluid values (4.20) and the respective properties of the two distinct phases themselves.
**Mixture properties** By definition of the volume of fluid value, the mixture's properties in any cell $C_i$ for which $\bar{H}_i=1$ are the ones of the reference phase. Respectively, the properties of the other phase are considered wherever $\bar{H}_i=0$. In the case of an interfacial cell $C_i$, where $0 < \bar{H}_i < 1$, the mixture's properties are derived for each cell from the value of $\bar{H}_i$ and the properties of the two phases themselves. In accordance with physical principles, the specific density $\rho$ would typically be taken as the weighted average $\rho = \rho_1 \bar{H}_i + \rho_2(1 - \bar{H}_i)$, while other properties as viscosity may involve a non-linear expression.

**Figure 4.14.** Properties of the single mixture defined from the two original phases.

<figure><img src="image_65.png" /></figure>

**Phase and state evolution** The space-dependent nature of the mixture is evolved cell-wise by advecting the discrete volume of fluid quantities (4.20) accordingly to the velocity field $v$ obtained from the state variables, that is,

$$ \frac{\partial \bar{H}}{\partial t} + \nabla \cdot (v \bar{H}) = 0. \qquad (4.22) $$

Similarly as in the Level Set method, the state (and thus the velocity field) is obtained by evolving the augmented system (4.1–4.2)–(4.22), the space-time dependent mixed phase properties of the mixed phase being only considered in the equation of state. The interface location is then reassessed accordingly to the updated interfacial control volume, as illustrated in Figure 4.15.

<figure><img src="image_66.png" /><figcaption>**Figure 4.15.** Redefinition of the interface location with respect to the updated volume of fluid values. Left: Previous state. Right: Updated state.</figcaption></figure>

**Method strengths and drawbacks** Advecting volume fractions from a velocity field, the Volume of Fluid method is mass conservative in essence. However, the volume of fluid field is piecewise constant, thus discontinuous, and the location of the interface is not exactly known within interfacial cells. Therefore, should the interface location be known, a reconstruction technique has to be employed. Although a variety of techniques exists, as e.g. [159, 160], they are usually restricted to low order and yield discontinuous interfaces.

A more recent technique, THINC-LS, retrieves a spatially determined interface while preserving mass conservation by combining the Volume of Fluid technique with the Level Set method via the THINC reconstruction.
### 4.2.3 The THINC-LS approach

The THINC reconstruction is a local correction procedure that is applied after each time step to the level set field. Using a phase-field function generated out of the level set values by a tangent hyperbolic profile, the interface correction is determined within each interfacial control volume with respect to conservative volume of fluid constraints. When considered as part of the THINC-LS hybrid scheme [20], it modifies the discrete model dynamics presented in Section 4.1.1 by incorporating the correction step after the individual interface and state evolutions, as depicted in Figure 4.16.

<figure><img src="image_67.png" /><figcaption>Figure 4.16. Semisequential update of the time-discretised model as performed by the THINC hybrid approach.</figcaption></figure>

**Time step** Being a hybrid approach, each time step comprehends the evolution of three entities; the state, the level set representation $\phi$ of the interface, and the volume of fluid representation $H$ of the interface. Each of the above evolutions is assumed to be performed on a same mesh $T^h = \{C_i\}_i$ accordingly to the system (4.1–4.2) for the state and to the dynamics (4.17) and (4.22) for the respective interface representations. More precisely, the discrete level set and volume of fluid fields are expressed at the time $t^n$ by

$$ \phi^n: \Omega \to \mathbb{R} \\ \text{and} \\ x \mapsto \phi_i^n(x), \quad \text{i s.t. } x \in C_i \\ x \mapsto \bar{H}_i^n, \quad \text{i s.t. } x \in C_i, $$

where $\phi_i^n$ is the local reconstruction of the level set within the cell $C_i$ and $\bar{H}_i^n$ is the value of the piecewise constant volume of fluid within the cell $C_i$.
The update of both quantities from the discrete time $t^n$ to $t^{n+1}$ is then computed independently for each control volume $C_i$ through the approximations of

$$ \frac{d\bar{H}_i}{dt} + \nabla \cdot (v \bar{H}_i) = 0 \quad (4.23) \quad \text{and} \quad \frac{\partial \phi_i}{\partial t} + v \cdot \nabla \phi_i = 0, \quad (4.24) $$

determining the degrees of freedom values characterising the level set $\phi_i^{n+1}$ and volume of fluid $\bar{H}_i^{n+1}$ discrete fields. There, the numerical method assessing (4.24) can a-priori be chosen freely among stable and consistent schemes, while the discretisation of (4.23) requires to be further conservative. Once evolved, the level set expression $\phi_i^{n+1}$ is then corrected with respect to the volume of fluid value of $\bar{H}_i^{n+1}$ in each interfacial control volume $C_i \in \mathcal{T}^h$.

**Correction procedure** After each time step, the level set expression $\phi_i^{n+1}$ is corrected in each cell $C_i$ with respect to the volume of fluid value $\bar{H}_i^{n+1}$ considered here as a constraint. This correction consists in six steps, detailed and depicted in Algorithm 1.

**Method strengths and drawbacks** As a corrective approach, the THINC-LS correction procedure preserves the topology flexibility of the level set and conserves the partial masses down to the local control volumes. However, by construction of the local fluid selector $\psi_i$ and by definition of the correction term $c_i$, the expression of the corrected level set yields a discontinuous interface. Furthermore, the correction being only an artificial shift up or down of the initial level set, the shape of the obtained interface has no physical meaning and is subject to the choices of the polynomial locator function's order and of its determination constraints.

A less visible and yet critical drawback lies in the nature of the corrected level set's expression. Indeed, the corrected level set can be expressed in the original level set discretisation framework if and only if the shape of the initial interface locator and of its offsets can. Therefore, whenever one likes to preserve the nature of the discrete representation through the correction procedure, a special attention has to be paid to the compatibility between the shape of the interface locator and the discretisation of (4.24). In this respect and to allow any interface locator, the classical THINC approach [20] evolves the level set with a semi-Lagrangian method. One could also be more tailored and choose an expression of $\psi_i$ so that for any $c_i \in \mathbb{R}$, $\psi_i + c_i$ can be represented in the level set discretisation associated to the cell $C_i$. However, this issue cannot be resolved in a general scope not specifying the underlying level set advection scheme.

In what follows, we propose to improve the THINC-LS approach by enforcing the interface cross-cell continuity by simply redefining the expressions of the locator and corrector functions.
1. Redefine the implicit representation of the interface location within $C_i$ to an arbitrary order $k \in \mathbb{N}$ representation by constructing a constrained polynomial $\psi_i \in \mathbb{Q}_k(\mathbb{R}^d)$ out of the local discrete level set field $\phi_i$. Choosing in particular the constraints as

$$ \frac{\partial^d \psi_i(x_m)}{\partial d_1 x_1 \partial d_2 x_2} = \frac{\partial^d \phi_i(x_m)}{\partial d_1 x_1 \partial d_2 x_2} $$

for $k \in \mathbb{N}$, $d_1, d_2 \in [0, d]$ such that $d_1 + d_2 = d$, and for $\{x_m\}_m \in [1, M]$ a set of $M \in \mathbb{N}$ evaluation points contained in $C_i$, the values of $\psi_i$ and its derivatives at any point $x \in \{x_m\}_m$ are consistent with the initial local level set reconstruction $\phi_i$. Furthermore, the so constrained interpolation function $\psi_i$, referred to in the following as *interface locator*, is of maximum steepness near the interface.

2. Add a corrective constant $c_i \in \mathbb{R}$ to the expression of $\psi_i$, whose value will be later adjusted with respect to the local volume of fluid constraint.

3. Construct a smoothed indicator function $\mathcal{H}_i: C_i \to [0, 1]$ tagging the phases within $C_i$ through the tangent hyperbolic phase-field profile

$$ \mathcal{H}_i = \frac{1}{2} (1 + \tanh(\epsilon\psi_i + \epsilon c_i)) $$

referred to in the following as *volume of fluid estimator*. There, $2\epsilon \in \mathbb{R}_+^*$ represents a chosen artificial interface thickness parameter.

4. Retrieve the *evaluated volume of fluid value* $\bar{H}_i \in [0, 1]$, representing the total quantity of the reference phase in $C_i$ as dictated by the corrected interface location represented by $\psi_i + c_i$, i.e.

$$ \bar{H}_i = \int_C \mathcal{H}_i(x) dx. $$

5. Determine the value of $c_i$ by matching the evaluated volume of fluid value with the actual volume of fluid constraint $\bar{H}_i$, that is,

Find $c_i$ such that $\bar{H}_i = \bar{H}_i$.

6. Update the local expression of the level set $\phi_i$ in $C_i$ by using the reciprocal *THINC* function, i.e.

$$ \phi_i(x) = \frac{1}{\epsilon} \tanh^{-1}(2\mathcal{H}_i(x) - 1) = \psi_i(x) + c_i. \quad (4.25) $$

<figure><img src="image_68.png" /></figure>

<figure><img src="image_69.png" /></figure>

<figure><img src="image_70.png" /></figure>

<figure><img src="image_71.png" /></figure>

**Algorithm 1.**

First four steps of the correction procedure in the one-dimensional case.
## 4.3 A continuity preserving procedure

The discontinuity problem encountered in [20] and observed in Section 4.2.3 is only due to the definition of the cell-wise locator function $\psi_i$ and to the nature of its correction $c_i$. By redefining their expressions, we propose a correction procedure that ensures both mass conservation and interface continuity.

After sketching the introduced idea in Section 4.3.1, we present in Section 4.3.2 intuitive definitions that can be used with continuous discretisations of (4.24) and (4.23) considering primal meshes of convex elements. They are then extended to the cases of non-convex control volumes, staggered settings, and discontinuous discretisations in Section 4.3.3. A set of admissibility conditions for designing general expressions of locator and corrector functions is then provided in Section 4.3.3.4, while possible numerical approximations are proposed in Section 4.3.4. Lastly, numerical results are presented in Section 4.3.5, followed in Section 4.3.6 by a critical assessment of the method.

### 4.3.1 The general idea in a nutshell

As in the *THINC-LS* approach, the conservative location of the interface is determined by an extraneous correction enforcing volume of fluid constraints. This correction is applied on the level set field $\phi$ in each interfacial cell of the volume of fluid discretisation, that is, within each cell $C_i$ where $\bar{H}_i \in (0, 1)$.

There, the initial non-conservative interface is represented as in [192, 193] by a parameter-dependent function that is constrained by locally reconstructed point-values of the level set and its derivatives. However, contrarily to [192], this reconstruction has not for objective to provide an arbitrary high order representation, but to enforce cross-cell continuity and to allow flexibility in the schemes design. Indeed, when different schemes are used to tackle the equations (4.23) and (4.24), the spatial discretisations of the volume of fluid and level set fields are not necessarily identical. The control volume $C_i$ may therefore not match any of the cells associated to the level set discretisation, and no local point-value reconstruction of the level set may have a support matching any of the given volume of fluid constraints. As the initial discretisation may be further discontinuous, a cross-cell continuous interface locator function $\psi_i$, local reconstruction of the level set field $\phi$ over $C_i$, is therefore constructed before applying the correction procedure on the local control volume.

The expression of this initial locator function is defined as to enforce the interface cross-cell continuity while only requiring a compact stencil for determining its parameters. Mass conservation is then achieved by adding a correction term that is tuned towards the associated volume of fluid constraint $\bar{H}_i$. There, we consider a parameter-dependent function $\psi_i^\Delta$ that vanishes on the control volume boundaries, hence preserving any pre-existing cross-cell interface continuity. Thus, as the initial interface locator is cross-cell continuous, the obtained corrected interface is conservative and globally continuous.
### 4.3.2 An intuitive continuous interface correction

In order to ease the presentation of the proposed technique, we first consider the two spatial dimensional case where both the discrete level set and volume of fluid fields are resolved on a same primal mesh of convex elements. There, we apply step by step the THINC-LS procedure and detail the proposed improvements, impacting only the second and third steps.

#### 4.3.2.1 Discretisation setting

**Notations** We consider $\Omega \subset \mathbb{R}^2$ the study domain where the spatial points $x \in \Omega$ have for coordinates $x = (x_1, x_2)$. Let then $\mathcal{T}^h$ be a primal mesh consisting of convex elements $K \subset \Omega$ and $\Omega := \cup_{K \in \mathcal{T}^h} \{K\}$ represent the computational domain. Each element $K$ has for boundary $\partial K$, a union of $n_K$ faces generically denoted by $f$. Moreover, during the correction process, $K$ is identified to a control volume $C_i$ for some given index $i \in \mathbb{N}$. We thus denote by $\Xi_i$ the set of faces $f \in \partial K$ that are crossed by the interface $\Lambda$ and refer to any face $f \in \Xi_i$ by the term “interfacial face”.

**Known information layout** We assume that the equation (4.24) is tackled with a primal scheme that offers a continuous representation of the level set field over $\Omega$ and provides a local pointwise reconstruction of the level set within each element $K$. That is, we assume that there exists a local level set function

$$ \phi^K: K \to \mathbb{R} \\ x \mapsto \phi^K(x) $$

determined from local degrees of freedom associated to the element $K$, and that the global level set function

$$ \phi: \Omega \to \mathbb{R} \\ x \mapsto \phi^K(x), \quad K \text{ s.t. } x \in K $$

is continuous across cells. We further assume that the volume of fluid value is given for each primal element $K$ and ruled by a discretely conservative evolution of (4.23), providing conservative constraints on the associated control volume $C_i := K$. A visual description of the spatial layout of the state, volume of fluid and level set information is provided in Figure 4.17.

<figure><img src="image_72.png" /><figcaption>(a) Used spatial notations</figcaption></figure>

<figure><img src="image_73.png" /><figcaption>(b) Variational layout</figcaption></figure>

**Figure 4.17.** Description of the spatially known information.
### 4.3.2.2 The modified correction procedure

As in this setting the representation of the level set is continuous over the primal mesh, yielding a cross-cell continuous interface location, we do not require any reconstruction of a local interface locator. The level set field is therefore directly corrected in each primal control volume $C_i := K$ by a local correction function $\psi_{i,c}^{\Delta}: C_i \to \mathbb{R}$. The corrected level set field then reads

$$ \psi: \Omega \to \mathbb{R} \\ x \mapsto \phi^K(x) + \psi_{i,c}^{\Delta}(x), \quad i \text{ s.t. } C_i = K \text{ and } x \in K, \qquad (4.26) $$

where the parameters of each correction function $\psi_{i,c}^{\Delta}$ are determined with respect to the local volume of fluid constraint $\bar{H}_i$. The adaptation of the THINC-LS procedure we present here therefore only lies in the expression of the locator function $\psi_i$ and the corrector function $\psi_{i,c}^{\Delta}$, impacting only the second and third steps. We however derive the complete procedure once again for the sake of self-sufficiency of the detailed algorithm.

Having defined all the prerequisites, let us now explicit the continuous and conservative correction procedure. We consider the correction step performed at a given discrete time $t^{n+1}$ and assume that the discrete level set $\phi^n$ and volume of fluid $\bar{H}^n$ fields are known and have been evolved to $\phi^{n+1}$ and $\bar{H}^{n+1}$ through primal discretisations of

$$ \frac{d\bar{H}_i}{dt} + \nabla \cdot (v\bar{H}_i) = 0 \quad (4.23) \quad \text{and} \quad \frac{\partial \phi_i}{\partial t} + v \cdot \nabla \phi_i = 0, \quad (4.24) $$

respectively. In other words, we explicit the contents of the orange box in the Figure 4.16 so that the correction preserves the interface continuity.

1. Select an interfacial element $K$ within which to correct the interface location and define the control volume as $C_i := K$. Retrieve then the pointwise reconstruction of the local level set $\phi^K$ as well as the average volume of fluid value $\bar{H}_i$ on $C_i$.

2. Define the locator function $\psi_i$ as being equal to the local level set reconstruction $\phi^K$.

3. Define the corrector function $\psi_{i,c}^{\Delta}: C_i \to \mathbb{R}$ as the polynomial defined as the pondered product of the faces' Cartesian equations, i.e.

$$ \psi_{i,c}^{\Delta}(x) = c \prod_{j \in [1, n_K]} (d_j x_2 - a_j x_1 - b_j). \qquad (4.27) $$

There, $a_j, d_j$ are the slope and $b_j$ the offset coefficients of the line extending the face $f_j \in \partial K$, and $c \in \mathbb{R}$ is a constant.

Those two definitions alone are sufficient to ensure a cross-cell continuous representation of the interface. Indeed, the interface locator is in particular cross-
cell continuous by assumption, and the corrector function is defined as to be continuous and vanishing on the cell boundary $\partial K$. The crossing locations of the interface on the cell boundary are thus preserved regardless the value of the constant $c$, that can be adjusted against any volume of fluid constraint without further assumption (see Section 4.3.2.3 for an existence statement). To this aim, an auxiliary volume of fluid estimator is constructed.

4. Set an artificial interface thickness $2\epsilon \in \mathbb{R}_+^*$.

5. Define a volume of fluid estimator $\mathcal{H}_{i,c}$ for the control volume $C_i$ by applying the THINC reconstruction to the corrected locator function $\psi_i + \psi_{i,c}^\Delta$ pondered by the interface thickness. Namely, define

$$ \mathcal{H}_{i,c} = \frac{1}{2} (1 + \tanh(\epsilon\psi_i + \epsilon\psi_{i,c}^\Delta)). \quad (4.28) $$

Following this definition, the integral of $\mathcal{H}_{i,c}$ over $C_i$ provides an estimated volume of fluid value that corresponds to the corrected interface location. Knowing in parallel the volume of fluid constraint associated to $C_i$ is therefore enough to enforce mass conservation, as one can match the evaluated volume of fluid value with the constraint $\bar{H}_i$.

6. Apply the conservative constraint associated with the cell $C_i$ to the volume of fluid estimator and resolve the parameter $c \in \mathbb{R}$ of the correction function accordingly. In other words,

$$ \text{Find } c \in \mathbb{R} \text{ such that } \frac{1}{|C_i|} \int_{C_i} \mathcal{H}_{i,c}(x) \, dx = \bar{H}_i. \quad (4.29) $$

The solution to (4.29) provides the parameter $c \in \mathbb{R}$ of the correction function $\psi_{i,c}^\Delta$ to be used within the cell $C_i$, and therefore the conservative locator function $\psi_i + \psi_{i,c}^\Delta$. There is then only left to update the values of the level set field $\phi$ by using the reciprocal volume of fluid estimator.

7. Redefine the level set function on the primal elements $K$ from the reciprocal volume of fluid estimator, namely,

$$ \begin{aligned} \forall x \in K, \quad \phi(x) &= \frac{1}{\epsilon} \tanh^{-1}(2\mathcal{H}_{i,c}(x) - 1) \\ &= \psi_i(x) + \psi_{i,c}^\Delta(x). \end{aligned} \quad (4.30) $$

At the end of the procedure, the mass conservation is enforced by the value of $c \in \mathbb{R}$, while the cross-cell interface continuity is enforced by the expression of $\phi^K$ and preserved by the definition of the correction functions $\psi_{i,c}^\Delta$. Furthermore, the interface thickness $2\epsilon$ that has been used in the correction procedure is not involved anymore in the corrected level set evaluation.
### 4.3.2.3 Analysis

The feasibility and relevance of the proposed improvements to the THINC-LS approach rely entirely on the shape given to the corrector function $\psi_{i,c}^\Delta$. Indeed, its sole definition is responsible for the existence of a corrected interface location that lies within the control volume $C_i$ and further satisfies the relation (4.29). We show that the proposed shape (4.27) fulfils those two essentials and subsequently derive an additional property that will be beneficial to the numerical approximation of the correction procedure.

**Existence of a conservative corrected interface location** By definition of the correction algorithm, the obtention of a conservative interface is conditional on the existence of a solution to (4.29).

The free parameter of the corrector function should therefore have a direct impact on the location of the interface within the control volume, deforming the zero level set of $\psi_i + \psi_{i,c}^\Delta$ without creating any artificial boundary crossing or modifying existing ones. In addition, the set of corrected interface locations that can be obtained through the parameter range must contain at least one interface that splits the control volume into areas whose proportion corresponds to $\bar{H}_i$ for any

$\bar{H}_i \in (0, 1)$. Together, those two conditions reduce to requiring the surjectivity of the evaluated volume of fluid function

<figure><img src="image_74.png" /><figcaption>**Figure 4.18.** Local correction of the interface while preserving the cross-cell continuity.</figcaption></figure>

$$ \bar{H}_i: \mathbb{R} \to \mathbb{R} \\ c \mapsto \frac{1}{|C_i|} \int_{C_i} \frac{1}{2} (1 + \tanh(\epsilon\psi_i(x) + \epsilon\psi_{i,c}^\Delta(x))) dx, \quad (4.31) $$

property that is ensured for the proposed function (4.27) by Proposition 4.3.1.

**Proposition 4.3.1 Feasibility of the volume of fluid constraint**

Let $C_i$ be the control volume to which is associated the volume of fluid constraint $\bar{H}_i \in (0, 1)$ and consider $c$ the parameter of the correction function $\psi_{i,c}^\Delta$ given by (4.27). Then,

$$ \forall \bar{H}_i \in (0, 1), \exists c \in \mathbb{R} \text{ s.t. } \frac{1}{|C_i|} \int_{C_i} \mathcal{H}_{i,c}(x) dx = \bar{H}_i. $$

The presented proof, based on the intermediate value theorem, is fairly simple but requires specific properties of the proposed corrector and volume of fluid estimator functions. Hence, we start by specifying them in the two following lemmas before presenting surjectivity arguments in Proof of Proposition 4.3.1.
**Lemma 4.3.1** *Properties of the correction function*

For any control volume $C_i \subset \Omega$ whose any of its faces $f_j \in \partial K$ is parametrised by the slope $a_j$, $d_j \in \mathbb{R}$ and offset $b_j \in \mathbb{R}$ coefficients, and for any constant $c \in \mathbb{R}^*$, the function

$$ \psi_{i,c}^{\Delta}: C_i \to \mathbb{R} \\ x \mapsto c \prod_{j \in [1, n_K]} (d_j x_2 - a_j x_1 - b_j), $$

• is never vanishing strictly within $C_i$,

• does not change sign within the control volume $C_i$,

• is bounded within $C_i$ and achieves its maximum absolute value strictly within $C_i$, i.e. $\max_{x \in C_i}\{|\psi_{i,c}^{\Delta}(x)|\} = \max_{x \in \tilde{C}_i}\{|\psi_{i,c}^{\Delta}(x)|\} := |c\gamma_i| \in \mathbb{R}$ with $|\gamma_i| > 0$.

**Proof.**

• Let us start by showing that the correction function $\psi_{i,c}^{\Delta}$ does not change sign strictly within the control volume $C_i$. For this purpose, we can directly use the definition by observing that $\psi_{i,c}^{\Delta}$ is a polynomial belonging to $\mathbb{Q}_{n_K}(K) = \mathbb{Q}_{n_K}(C_i)$ and consists in a weighted product of the factors $(d_j x_2 - a_j x_1 - b_j)$ for $j \in [1, n_K]$. Therefrom, one can apply the factor theorem and assert that whenever $c \neq 0$, the zeros of $\psi_{i,c}^{\Delta}$ correspond exactly to the union of those of the factors. The kernel of $\psi_{i,c}^{\Delta}$ is thus given by $\cup\{x = (d_j x_2, a_j x_1 + b_j), j \in [1, n_K]\}$, describing the union of the lines extending the faces of $C_i$. As here the control volume $C_i$ is assumed to be convex and of positive area, we can directly infer that none of those line intersect with the interior of the cell. Consequently, whenever $c \neq 0$, the function $\psi_{i,c}^{\Delta}$ never vanishes within any convex control volume $C_i$.

• The property of $\psi_{i,c}^{\Delta}$ to have a constant sign strictly within $C_i$ is an immediate consequence of the above statement. Indeed, as the function $\psi_{i,c}^{\Delta}$ is a continuous polynomial never vanishing in the interior of $C_i$, it follows by the contraposition of the intermediate value theorem that the function $\psi_{i,c}^{\Delta}$ never changes sign strictly within the control volume $C_i$.

• The maximum of $|\psi_{i,c}^{\Delta}|$ being bounded and achieved strictly within $C_i$ is equally straightforward. Indeed, as the control volume $C_i$ is compact and $\psi_{i,c}^{\Delta}: C_i \to \mathbb{R}$ continuous, the function is bounded on $C_i$. Furthermore, we saw above that $\psi_{i,c}^{\Delta}$ is never vanishing strictly within $C_i$ whenever $c \neq 0$. The maximum absolute value of $\psi_{i,c}^{\Delta}$ is therefore positive and real valued. Moreover, as by definition $\psi_{i,c}^{\Delta}|_{\partial C_i} = 0$, the maximum absolute value is necessarily achieved strictly within $C_i$. Thus, $\max_{x \in C_i}\{|\psi_{i,c}^{\Delta}(x)|\} = \max_{x \in \tilde{C}_i}\{|\psi_{i,c}^{\Delta}(x)|\} := |c\gamma_i| > 0$. ■

The obtained properties of the corrector function engender those of the estimator function, crucial to proving the well-posedness of the problem (4.29).
**Note** In order to ease the analysis of the estimator function (4.28), let us now break the expression of the corrector function (4.27) into its two natural parts, the first one depending only on the free parameter and the second depending only on the shape of the control volume. Expressly, we write the corrector function as $\psi_{i,c}^{\Delta} := c\omega_{i}^{\Delta}$ where $c \in \mathbb{R}$ and $\omega_{i}^{\Delta}$ is defined by the parameter-free function

$$ \omega_i^{\Delta}: C_i \longrightarrow \mathbb{R} \\ x \longmapsto \prod_{j \in [1, n_K]} (d_j x_2 - a_j x_1 - b_j), \qquad (4.32) $$

$d_j$, $a_j$ and $b_j$ being predetermined by the control volume $C_i$ as the slope and offset coefficients of the line extending the face $f_j \in \partial K$. So defined, note that the function $\omega_i^{\Delta}$ itself still satisfies the properties derived in Lemma 4.3.1 when considering the constant $c$ as equal to 1. ▲

**Lemma 4.3.2** *Properties of the volume of fluid estimator function*

For any control volume $C_i \subset \Omega$, any constant $c \in \mathbb{R}$, and any interface semi-width $\epsilon \in \mathbb{R}_+$, the volume of fluid estimator function

$$ \mathcal{H}_{i,c}: C_i \longrightarrow \mathbb{R} \\ x \longmapsto \frac{1}{2} (1 + \tanh(\epsilon\psi_i(x) + \epsilon c\omega_i^{\Delta}(x))) $$

* is non-negative, i.e. $\forall x \in C_i, \mathcal{H}_{i,c}(x) \ge 0$,

* is bounded, and more particularly $\max\{\mathcal{H}_{i,c}(x), x \in C_i\} \le 1$,

* is continuous and Lebesgue integrable, i.e. $\mathcal{H}_{i,c} \in L^1(C_i) \cap C^0(C_i)$.

**Proof.**

* (First two points) By definition of the level set locator and corrector functions, one has that $\psi_i$ and $\psi_{i,c}^{\Delta}$ are real-valued at any point $x \in C_i$. Hence, we get by definition of the hyperbolic tangent that

$$ -1 \le \tanh(\epsilon\psi_i(x) + \epsilon c\omega_i^{\Delta}(x)) \le 1 \quad \forall x \in C_i, \forall \epsilon \in \mathbb{R}_+ \text{ and } \forall c \in \mathbb{R}. $$

Reconstructing the expression of $\mathcal{H}_i$ by elementary operations, we retrieve

$$ 0 \le \frac{1}{2} (1 + \tanh(\epsilon\psi_i(x) + \epsilon c\psi_i^{\Delta}(x))) \le 1 \quad \forall x \in C_i, \forall \epsilon \in \mathbb{R}_+ \text{ and } \forall c \in \mathbb{R}, $$

bounding $\mathcal{H}_{i,c}$ on $C_i$ and implying straightforwardly the two first properties.

* (Last point) To begin with, we notice that the function $\mathcal{H}_{i,c}$ is continuous for any $c \in \mathbb{R}$ as being defined as the addition and composition of continuous functions. From there, one can simply use the above point and state that as $\mathcal{H}_{i,c}$ is bounded and $C_i \in \Omega$ compact, $\int_{C_i} \mathcal{H}_{i,c}(x)dx$ is real valued. Thus, $\mathcal{H}_{i,c}$ is further Lebesgue integrable and $\mathcal{H}_{i,c} \in C^0 \cap L^1(C_i)$. ■

Given those regularity properties, we can now proceed to proving the surjectivity of the function $\bar{\mathcal{H}}_i$ when considering the corrector function (4.27).
**Proof of Proposition 4.3.1**

In order to prove the surjectivity of the function $\bar{\mathcal{H}}_i$ defined by (4.31), it is enough to show that $\bar{\mathcal{H}}_i$ is continuous and satisfies

$$ \left\{ \begin{array}{ll} \lim_{c \to -\infty} \bar{\mathcal{H}}_i(c) = 0 & \text{and} \quad \lim_{c \to +\infty} \bar{\mathcal{H}}_i(c) = 1 \quad \text{when } \omega_i^\Delta \ge 0 \\ \\ \lim_{c \to -\infty} \bar{\mathcal{H}}_i(c) = 1 & \text{and} \quad \lim_{c \to +\infty} \bar{\mathcal{H}}_i(c) = 0 \quad \text{when } \omega_i^\Delta \le 0, \end{array} \right. $$

the surjectivity of $\bar{\mathcal{H}}_i$ then following from the intermediate value theorem.

• The continuity of $\bar{\mathcal{H}}_i: \mathbb{R} \to \mathbb{R}$ is an immediate consequence of the continuity of its integrand $\mathcal{H}_{i,c}$ given by Lemma 4.3.2. Indeed, the function $\mathcal{H}_{i,c}$ is itself continuous as being formed by additions and compositions of continuous functions. Besides, as its definition domain $C_i \subset \mathbb{R}^2$ is compact, it is further absolutely continuous. Thus, its integral over $C_i$ is continuous with respect to the parameter $c$, implying by definition the continuity of $\bar{\mathcal{H}}_i$.

• Let us now investigate the limits of $\bar{\mathcal{H}}_i$ as $c \to \pm\infty$ when $\omega_i^\Delta$ is non-negative. Analogous, the case $\omega_i^\Delta \le 0$ will not be detailed here for the sake of concision.

▷ We first derive the limit of $\bar{\mathcal{H}}_i$ as $c \to +\infty$ by investigating the expression of the integrand $\mathcal{H}_{i,c}$ and using its properties for obtaining integral estimates depending on the value of $c$. More precisely, we know by Lemma 4.3.2 that for any constant $c \in \mathbb{R}$ the function $\mathcal{H}_{i,c}$ is Lebesgue integrable, non-negative and bounded by 1. Hence, we can apply Fatou's lemma and assert that

$$ \frac{1}{|C_i|} \int_{C_i} \lim_{c \to +\infty} \frac{1}{2} \left( 1 + \tanh\left(\epsilon\psi_i(x) + \epsilon c \omega_i^\Delta(x)\right) \right) dx \le \lim_{c \to +\infty} \frac{1}{|C_i|} \int_{C_i} \mathcal{H}_{i,c}(x) dx, $$

where we have made explicit the expression of $\mathcal{H}_{i,c}$ in the left hand side. There, using the non-negativity assumption over $\omega_i^\Delta$ and recalling from Lemma 4.3.1 that in this case $\omega_i^\Delta > 0$ strictly within $C_i$, we obtain

$$ \lim_{c \to +\infty} \frac{1}{2} \left( 1 + \tanh\left(\epsilon\psi_i(x) + \epsilon c \omega_i^\Delta(x)\right) \right) = 1, \quad \forall \epsilon \in \mathbb{R}_+, \forall x \in \mathring{C}_i. $$

Indeed, the functions $\omega_i^\Delta$ and $\psi_i$ are continuous and therefore bounded on their compact support $C_i$. They can thus be treated as constants, and the limit follows by definition of the hyperbolic tangent. The lower bound then reads

$$ 1 \le \lim_{c \to +\infty} \frac{1}{|C_i|} \int_{C_i} \mathcal{H}_{i,c}(x) dx. $$

In turn, the upper bound is simply obtained by plugging the global maximum of $\mathcal{H}_{i,c}$, given by Lemma 4.3.2 and reading

$$ \forall c \in \mathbb{R}, \quad \max_{x \in C_i} \{\mathcal{H}_{i,c}(x)\} \le 1, $$
in the integrand's expression. Namely,

$$ 
\begin{aligned}
\lim_{c \to +\infty} \frac{1}{|C_i|} \int_{C_i} \mathcal{H}_{i,c}(x) \, dx &\le \lim_{c \to +\infty} \frac{1}{|C_i|} \int_{C_i} \max_{x \in C_i}\{H_{i,c}(x)\} \, dx \\
&\le \lim_{c \to +\infty} \frac{1}{|C_i|} \int_{C_i} 1 \, dx = 1,
\end{aligned}
$$

allowing us to conclude that indeed $\lim_{c \to +\infty} \bar{\mathcal{H}}_i = 1$ by the squeeze theorem.

▷ Let us now derive the limit of $\bar{\mathcal{H}}_i$ as $c_i \to -\infty$. In this case, we can bound directly the integral from above by using only the expression of $\mathcal{H}_{i,c}$ itself. Indeed, by assumption for $\psi_i$ and definition for $\omega_i^\Delta$, those functions are continuous. Therefore, as $C_i$ has a compact support their values are bounded by two constants $\alpha_\psi \in \mathbb{R}$ and $\alpha_\omega \in \mathbb{R}$. In other words,

$$ \forall x \in C_i, \quad |\psi_i(x)| \le \alpha_\psi \quad \text{and} \quad |\omega_i^\Delta(x)| \le \alpha_\omega. $$

The hyperbolic tangent function being strictly increasing, one can then bound the value of $\mathcal{H}_{i,c}$ at any point $x \in C_i$ by writing

$$ \mathcal{H}_{i,c}(x) = \frac{1}{2}(1 + \tanh(\epsilon\psi_i(x) + \epsilon c\omega_i^\Delta(x))) \le \frac{1}{2}(1 + \tanh(\epsilon\alpha_\psi + \epsilon c\alpha_\omega)), $$

yielding naturally by non-negativity of the volume of fluid estimator $\mathcal{H}_{i,c}$

$$ \frac{1}{C_i} \int_{C_i} \frac{1}{2} (1 + \tanh(\epsilon \psi_i(x) + \epsilon c \omega_i^\Delta(x))) dx \le \frac{1}{2} (1 + \tanh(\epsilon \alpha_\psi + \epsilon c \alpha_\omega)). $$

Considering the limit of both terms as $c \to -\infty$, we then obtain

$$ \lim_{c \to -\infty} \frac{1}{C_i} \int_{C_i} \frac{1}{2} (1 + \tanh(\epsilon \psi_i(x) + \epsilon c \omega_i^\Delta(x))) dx \le \lim_{c \to -\infty} \frac{1}{2} (1 + \tanh(\epsilon \alpha_\psi + \epsilon c \alpha_\omega)) $$

and thus retrieve by definition of the hyperbolic tangent the upper bound

$$ \lim_{c \to -\infty} \frac{1}{|C_i|} \int_{C_i} \mathcal{H}_{i,c}(x) \, dx \le 0. $$

Recalling further that $\mathcal{H}_{i,c} \ge 0$ for any $c \in \mathbb{R}$, we observe that $\bar{\mathcal{H}}_i$ is necessarily non-negative and apply the squeeze theorem, yielding $\lim_{c \to -\infty} \bar{\mathcal{H}}_i(c) = 0$.

▷ Summarising the above two cases, we can state that whenever $\omega_i^\Delta \ge 0$,

$$ \lim_{c_i \to -\infty} \bar{\mathcal{H}}_i = 0 \quad \text{and} \quad \lim_{c_i \to +\infty} \bar{\mathcal{H}}_i = 1. $$

The case of $\psi_i^\Delta$ turns to be analogous by writing $c = -c$ and $\omega_i^\Delta = -\omega_i^\Delta$. There, we derive similarly $\lim_{c_i \to -\infty} \bar{\mathcal{H}}_i = 1$ and $\lim_{c_i \to +\infty} \bar{\mathcal{H}}_i = 0$, concluding the proof. ■

Now proven, Proposition 4.3.1 provides the existence of a solution to the problem (4.29), and therefore the existence of a corrected interface for any value of the conservative constraint $\bar{\mathcal{H}}_i$. If our approach can thus be validated in principle,
ciple, we still have to attest of its robustness. We therefore derive the strict monotonicity of the corrector function $\bar{H}_i$, ensuring the uniqueness of the corrected interface and the continuous dependency of its associated parameter $c$ on the volume of fluid constraint $\bar{H}_i$.

**Lemma 4.3.3** *Monotonicity of the volume of fluid estimator*

The function $\bar{H}_i: \mathbb{R} \to \mathbb{R}$ actually takes values in $(0, 1)$, and

* is strictly increasing whenever $\omega_i^\Delta \ge 0$,

* is strictly decreasing whenever $\omega_i^\Delta \le 0$.

*Proof.*

* (Monotonicity) Let us first derive the monotonicity of $\bar{H}_i$ by studying the sign of its derivative. As $\bar{H}_i$ is defined by the integral of $H_{i,c}$ over $C_i$, we can directly apply the Leibnitz integral rule and make use of the integrand's properties. There,

$$ \begin{aligned} \frac{d}{dc} \left( \frac{1}{|C_i|} \int_{C_i} \mathcal{H}_{i,c}(x) dx \right) &= \frac{1}{2|C_i|} \int_{C_i} \frac{\partial}{\partial c} \left( 1 + \tanh(\epsilon \psi_i(x) + \epsilon c \omega_i^\Delta(x)) \right) dx \\ &= \frac{1}{2|C_i|} \int_{C_i} \omega_i^\Delta(x) \sech^2(\epsilon \psi_i(x) + \epsilon c \omega_i^\Delta(x)) dx. \quad (4.33) \end{aligned} $$

Thus, as we saw in Lemma 4.3.1 that the sign of $\omega_i^\Delta$ is constant within $C_i$, bounds ensuring monotonicity can be derived accordingly to the sign of $\omega_i^\Delta$.

▷ Let us first consider the case where $\omega_i^\Delta$ is non-negative. To begin with, one can observe that by definition of the hyperbolic secant and the fact that $\psi_i$ and $\omega_i^\Delta$ are real valued, $\sech^2(\epsilon \psi_i + c\epsilon \omega_i^\Delta) \in (0, 1)$ and is thus greater than zero. Hence, $C_i$ being compact there exists a constant $\eta_1 \in \mathbb{R}_+^*$ so that

$$ \max_{x \in C_i} \{\sech^2(\epsilon \psi_i(x) + \epsilon c \omega_i^\Delta(x))\} \ge \eta_1 > 0, \quad (4.34) $$

and we have that

$$ \frac{d}{dc} \left( \frac{1}{|C_i|} \int_{C_i} \mathcal{H}_{i,c}(x) dx \right) \ge \frac{\eta_1}{2|C_i|} \int_{C_i} \omega_i dx. \quad (4.35) $$

Furthermore, we know by Lemma 4.3.1 that the function $\omega_i^\Delta$ is never changing sign nor vanishing strictly within $C_i$, and that $\max\{|\omega_i^\Delta(x)|, x \in \mathring{C}_i\} = |\gamma_i|$ for some $\gamma_i \in \mathbb{R}$. As we further assume that $\omega_i^\Delta \ge 0$, we can infer that $\gamma_i > 0$ and reduce (4.35) to

$$ \frac{d}{dc} \left( \frac{1}{|C_i|} \int_{C_i} \mathcal{H}_{i,c}(x) dx \right) \ge \eta_1 \gamma_i > 0, \quad (4.36) $$

guaranteeing that $\bar{H}_i$ is strictly increasing along with the value of $c$.
▷ The case were $\omega_i^\Delta \le 0$ can be treated analogously. Indeed, arguing with the identical arguments as for obtaining (4.34), we get the existence of a constant $\eta_2 \in \mathbb{R}^*$ that verifies

$$ \max_{x \in C_i} \{\mathrm{sech}^2(\epsilon \psi_i(x) + \epsilon c \omega_i^\Delta(x))\} \le \eta_2 < 1. $$

Using this estimate directly in (4.33), we retrieve

$$ \frac{d}{dc} \left( \frac{1}{|C_i|} \int_{C_i} \mathcal{H}_{i,c}(x) dx \right) \le \frac{\eta_2}{2|C_i|} \int_{C_i} \omega_i dx. \quad (4.37) $$

Thus, as here the function $\omega_i^\Delta$ is negatively valued everywhere within the $C_i$, control volume that is a subset of $\Omega$ having a non-null measure, we get

$$ \frac{d}{dc} \left( \frac{1}{|C_i|} \int_{C_i} \mathcal{H}_{i,c}(x) dx \right) < 0, \quad (4.38) $$

guaranteeing that $\bar{\mathcal{H}}_i$ is strictly decreasing along with the value of $c$.

• (Image space) By the above point, we know that the function $\bar{\mathcal{H}}_i$ is strictly monotonous over $\mathbb{R}$. Therefore, recalling the limit bounds obtained in Proof of Proposition 4.3.1, the image space of $\bar{\mathcal{H}}_i$ is reduced to $(0, 1)$. ■

As a consequence of Lemma 4.3.3, there always exists a unique corrector function associated to a given constraint $\bar{\mathcal{H}}_i \in (0, 1)$. Not requiring any condition on the initial interface locator other than its continuity, the presented technique is therefore able to correct any initial interface towards any volume of fluid constraint within any control volume, regardless the amount of correction needed to achieve conservation. Moreover, the presented technique is robust towards possibly small numerical perturbations, the core problem of the correction procedure (4.29) being well-posed.

**Proposition 4.3.1** **Well-posedness**

Considering the corrector function (4.27), the problem (4.29) is well-posed.

**Proof.**

The proof is immediate given the properties of (4.27). Indeed, we saw in Proposition 4.3.1 that the function $\bar{\mathcal{H}}_i: \mathbb{R} \to (0, 1)$ is continuous and in Lemma 4.3.3 that it is strictly monotonous. Therefore, it is bijective and the problem (4.29) admits a unique solution depending continuously on the value of the volume of fluid constraint $\bar{\mathcal{H}}_i$. ■

The well-posedness of the problem (4.27) states in particular that the value of the constant $c$ adapts continuously to the constraint value $\bar{\mathcal{H}}_i$. Therefore, as the corrector function $\psi_{i,c}^\Delta$ is also continuous, the spatial location of its zero level set changes smoothly, and we can investigate the spatial properties of the obtained corrected interfaces.
**Spatial properties of the corrected interface** Let us now drive our interest towards the spatial properties of the interface given as the zero level set of the corrected function $\psi_i + \psi_{i,c}^{\Delta}$. We aim there to justify the cross-cell continuity of its zero level set and get a glimpse of its topological behaviour.

To begin with, let us start by confirming that the main goal of achieving a conservative and cross-cell continuous interface has been reached. To this aim, we summarise the above discussions in the following proposition.

**Proposition 4.3.2** *Interface continuity and conservation principle*

Given the conservative volume of fluid constraint $\bar{H}_i$ associated to any control volume $C_i$, the locally corrected interface defined as the zero level set of the function $\psi_i + \psi_{i,c}^{\Delta}$ is conservative and cross-cell continuous.

*Proof.*

* Let us first focus on the conservation principle. According to the course of the correction procedure, the compliance of the corrected interface location with the discrete conservation principle is directly ensured at the last step. Indeed, as we assume here that the underlying numerical scheme evolving (4.23) is conservative, we have by Proposition 4.3.1 that (4.29) always admits a solution on any control volume $C_i$ yielding to the definition of a conservative interface location on $C_i$. Furthermore, as we are considering a fully primal setting, any control volume $C_i$ match some given element $K \in \mathcal{T}^h$ on which the initial level set field has been determined. The interface locator function $\psi_i$ is therefore built and modified on each primal cell $K = C_i$ to become consistent with the conservative volume fraction. Hence, the mass conservation is achieved down to the cell level of the level set discretisation, and (4.11) is fulfilled.

* As shortly stated when describing the correction procedure, the cross-cell continuity is achieved by construction. Indeed, assuming that the initial interface locator function $\psi_i$ is already cross-cell continuous, defining the interface corrector $\psi_{i,c}^{\Delta}$ as a continuous function that vanishes on the boundary is enough to preserve all the intersection points of the interface with the boundary of the control volume. The cross-cell continuity of the interface is therefore preserved through all the correction procedure. ■

While preserving the continuity of the interface within the cell and across its boundaries, the proposed shape of the interface corrector makes the procedure entirely local, not impacting the interface's representation in neighbouring cells. In particular, the above statement implies that no spurious interface crossing point will be created on the boundary of the control volume.
**Corollary 4.3.1** *Behaviour of the local interface within its cell*

Each connected component of a locally corrected interface is continuous
within its supporting cell, and no spurious boundary crossing is created
during the correction procedure.

**Proof.** • (Continuity) Already cross-cell continuous by construction, the corrected interface is also continuous within the cell. Indeed, any primal cell $K$ is here directly identified to some control volume $C_i$, a bounded subset of $\mathbb{R}^2$. There, the locally defined interface locator and corrector functions are continuous. Thus, the zero level set of their sum has for connected components subsets of $C_i$ that are either closed and contained in the cell or crossing it from one subset of its boundary to another. Especially, in the case where the zero level set of the initial locator function contains only one-dimensional curves, classical scope of this work, the corrected interface has for connected components single curves that are either closed or having their end points on the cell's boundary. Thus, each connected component of the interface is continuous within $C_i$.

• (No spurious boundary crossing) The corrected interface not crossing the boundary at more, less or different locations than the original ones is again a consequence of the chosen shape of the corrector function $\psi_{i,c}^{\Delta}$. Indeed, vanishing on the boundary being continuous and not changing sign within $C_i$, the number of boundary crossing points remains the same.

Although not creating any spurious boundary crossing point, we saw that the
correction procedure is able to correct any continuous initial interface locator
function for any volume of fluid constraint, regardless its shape. In particular,
initial locator functions not having any single zero level set within the cell
can still be corrected. There, the corrector function is enforcing the conserv-
vative constraint by creating an interface whose connected components are
contained within the control volume, thus not creating any boundary crossing
point. On the same vein, the sign-constant polynomial nature of the corrector
function yields specific patterns of the corrected interface location that depends
both on the topology of the initial interface and on the needed correction
amount. Not yet rigorously proven, those observations are simply reported as
numerical results in Section 4.3.5.

Lastly, let us point out once again that this procedure holds for any contin-
uous shape of the initial interface locator and does not depend explicitly on the
underlying numerical schemes. More than defining a conservative and contin-
uous interface, the presented approach is therefore flexible towards the mesh
geometry and discretisation setting, only requiring the existence of control
volumes that cover the computational domain $\Omega$ and to which are associated
conservation constraints.
### 4.3.3 Extensions

The flexibility of the THINC-LS approach combined with the local definition of the corrector function allows us to extend the presented approach to more general cases. We go through crucial extensions, ranging from the drop of the continuity requirement on the initial level set field to non-convex control volumes and staggered discretisations. As an opening, we present a set of sufficient conditions for defining continuity-preserving corrector functions, paving the way for designs suited e.g. to specific physical or geometrical constraints.

#### 4.3.3.1 Extension to non-convex control volumes

In order to preserve the interface continuity, the presented modification to the THINC-LS procedure requires in particular the correction function $\psi_{i,c}^{\Delta}$ to have a constant sign within the control volume $C_i$ and to vanish on its boundary. Even though the proposed expression (4.27) always fulfils the second condition, being part of the vanishing ideal of $\partial C_i$, the first property does not hold for non-convex control volumes. Indeed, as the corrector function $\psi_{i,c}^{\Delta} \in Q_{n_K}(C_i)$

is explicitly built so that the lines extending the faces of the control volume provide its zeros, its construction over any non-convex control volume would result in the existence of zero lines in its interior. The regularity of the polynomial would then imply a sign change of $\psi_{i,c}^{\Delta}$ within the cell, failing to meet the requirements.

<figure><img src="image_75.png" /><figcaption>**Figure 4.19.** Maximum volume fraction that can be reached by a corrected interface bluit directly from (4.27) on a non-convex cell.</figcaption></figure>

As particular consequence, a collection of polygonal subcells where $\psi_{i,c}^{\Delta}$ has constant sign will be created. Thus, even if on each of those subcell the function $\psi_{i,c}^{\Delta}$ meets the Lemma 4.3.1, the range of volume fraction values that can be achieved through the correction process would be upper bounded by the cumulative area of the subcells that support the initial interface. The existence of a solution to the problem (4.29) could therefore not be guaranteed anymore in case of large corrections (see an illustration in Figure 4.19).

To overcome this issue, it is enough to define a bijective mapping $\zeta: C^{ref} \to C_i$ that deforms a convex polygon of reference $C^{ref}$ to the considered control volume $C_i$ so that the transformation of any path contained in $C^{ref}$ lies in $C_i$. However, the design of such transformation usually requires smooth continuous bijective mappings and is computationally costly. Thus, as a proof of concept only, we consider $\zeta$ as the smooth harmonic map proposed in [194], applied between the regular polygon having the same number of faces as the considered control volume $C_i$ and $C_i$ itself (see an illustration in Figure 4.20).
<figure><img src="image_76.png" /><figcaption>**Figure 4.20.** Illustration of a continuous mapping ensuring that the path transformed from any arc contained in $C_i$ is contained in the reference shape, and vis-versa.</figcaption></figure>

More precisely, assuming that the control volume is a generic n-gon, we can define the reference shape $C^{ref}$ as the regular polygon having for vertices

$$ \{e^{i \frac{2\pi}{n+1}k}\}_{k \in [1, n+1]} $$

in polar coordinates. The corrector function can then be defined as

$$ \psi_{i,c}^{\Delta}: C_i \to \mathbb{R} \\ x \mapsto \psi_i^{\Delta,ref}(\zeta^{-1}(x)), \qquad (4.39) $$

where $\zeta$ is given component-wise as in [194] through

$$ \left\{ \begin{array}{l} \Delta \zeta_l = 0 \\ \zeta_l|_{\partial C_i} = \mathcal{L}_l, \end{array} \right. \quad l \in [1, d], $$

the boundary term $\mathcal{L}: \partial C^{ref} \to \partial C_i$ being any bijective boundary mapping for which a possible choice is described in Section 4.3.5.1. The corrector function $\psi_{i,c}^{\Delta,ref}$

built on the regular polygon $C^{ref}$ can then be defined as previously by

<figure><img src="image_77.png" /><figcaption>**Figure 4.21.** Convex reference polygons for $n=3, 5, 6$.</figcaption></figure>

$$ \psi_{i,c}^{\Delta,ref}: C^{ref} \to \mathbb{R} \\ x \mapsto c \prod_{j \in [1, n_K]} (d_j x_2 - a_j x_1 - b_j), $$

where $a_j, d_j$ are the leading and $b_j$ the offset coefficients of the line extending the face $f_j \in \partial K$, and satisfies the set of properties described in Lemma 4.3.1.

Thereby, the same continuous and conservative interface correction procedure as the one detailed in Section 4.3.2 can be applied. Indeed, the mapping $\zeta$ is smooth, bijective and all the trajectories contained within the reference polygon $C^{ref}$ have their image contained in the non-convex control volume $C_i$. The properties of the reference corrector function $\psi_{i,c}^{\Delta,ref}$ thus naturally transfer to $\psi_{i,c}^{\Delta}$, and the well-posedness of the problem (4.29) constructed on the expression (4.39) of $\psi_{i,c}^{\Delta}$ can be established. A conservative and continuous interface entirely contained in $C_i$ can therefore be retrieved by the proposed approach, even when the conservative volume of fluid constraint is given on non-convex control volumes.
### 4.3.3.2 Extensions to discontinuous and high order representations of the initial level set field

So far, as defined in Section 4.3.2 and Section 4.3.3.1, the correction procedure only preserves the cross-cell continuity of the interface, not bringing it in when the interface is initially discontinuous. Therefore, dropping the assumption on the continuous initial level set field requires to change the definition of the initial locator function $\psi_i$ to a local reconstruction of the level set enforcing cross-cell continuity.

A straightforward definition of a local reconstruction that interpolates the level set in such a way relies on a simple adaptation of the arbitrary high order THINC-LS approach [192]. Indeed, there, the interface locator is defined as a polynomial constrained by the level set field, providing an arbitrary high order interface representation within the considered cell. Following this line, we simply modify the shape of the interpolation polynomial in order to enforce by construction the continuity of the interface across the cell boundary.

To this aim, we design an interface locator $\psi_i$ whose expression on each interfacial face is independent of the actual control volume $C_i$ that supports it. Aiming in addition to keep the determination stencil local to $C_i$, we define $\psi_i$ as a polynomial whose expression on any interfacial face depends on that face's parametrisation alone. Namely, denoting $\Xi_i$ the set of all interfacial faces of the control volume $C_i$, the proposed interface locator reads

$$ \psi_i: C_i \to \mathbb{R} \\ x \mapsto \sum_{j \in \Xi_i} p_k^j(\gamma_j(x)) \prod_{\substack{i \in \Xi_i \\ i \neq j}} (d_i x_2 - a_i x_1 - b_i), \quad (4.40) $$

where $p_k^j \in Q_k(\mathbb{R}^{d-1})$ is a polynomial of degree at most $k \in \mathbb{N}$, $a_i$, $d_i$ are the leading and $b_i$ the offset coefficients of the line extending the face $i$, and $\gamma_i$ is any of its parametrisation. Naming $x_A$ and $x_B$ the endpoints of a given face $j$, one may typically consider the distance based representation

$$ \gamma_j(x) = \frac{\|x - x_A\|_2}{\|x_A - x_B\|_2}, \quad (4.41) $$

defining $p_k$ as a radial function centred in $x_A$ whenever a robust representation is desirable, or simply

$$ \left\{ \begin{array}{ll} \gamma_j(x) = x_2 & \text{if } d_j = 0 \\ \gamma_j(x) = \frac{a_j}{d_j} x_1 + \frac{b_j}{d_j} & \text{otherwise} \end{array} \right. \quad (4.42) $$

when one prefers to consider an expression that is similar to the definition of the corrector function $\psi_{i,c}^\Delta$. Either way, the product term of (4.40) cancels the contribution of polynomial taking care of the representation on any other interfacial face. Furthermore, the number of interfacial faces defines the minimum polynomial order, while the degree of $p_k$ prescribes the available freedom for fitting the polynomial to the possibly discontinuous initial level set field $\phi$.
**Choice of the constraints characterising the locator function**
Once the shape (4.40) of the locator function $\psi_i$ has been selected, characterising its expression on the control volume $C_i$ upon the initial level set field $\phi$ comes down to finding the coefficients of each polynomial $p_k^j, j \in \Xi_i$. To this aim, and in order to ensure cross-cell continuity of the locator function, we constrain the polynomial accordingly to the local level set reconstructions $\phi^K$ whose support contains any of the faces $f \in \partial C_i$. More precisely, we consider the local set of reconstructions

$$ \{\phi^K: K \to \mathbb{R}, K \in \mathcal{T}^h \text{ s.t. } K \cap \partial C_i \neq \emptyset\} \qquad (4.43) $$

and select constraints that ensure both the compactness and relevance of the characterisation. We therefore focus specifically on constraints acting on the faces $j$ that force the zero level set of the locator function to be consistent with the boundary crossings of the initial interface as described by the possibly discontinuous level set field. In practice, we simply use the same set of point-value constraints as the one proposed in [20], set up however at a finite set of evaluation points $\{x_m\}_m$ that belongs to the cell boundary $\partial C_i$. Namely, denoting $\mathcal{K} := \{K \text{ s.t. } x_m \in \mathcal{K}\}$ and assuming a discontinuous level set discretisation, we impose the constraints

$$ \frac{\partial^d \psi_i}{\partial^{l_1} x_1 \partial^{l_2} x_2}(x_m) = \frac{1}{|\mathcal{K}|} \sum_{K \in \mathcal{K}} \frac{\partial^d \phi^K}{\partial^{l_1} x_1 \partial^{l_2} x_2}(x_m), \quad \text{where} \quad \substack{l_1, l_2 \in [0, d], \\ \text{s.t. } l_1+l_2=d} \qquad (4.44) $$

and where the set $\{x_m\}_m$ is defined so that the number of evaluation points distributed on the face $j \in \partial C_i$ matches the freedom amount of the polynomial $p_k^j$, that is e.g. $k+1$ in two dimensions. Alternatively, the behaviour of the locator function can be prescribed more directly by

$$ \frac{\partial \psi_i}{\partial x_l}(x_m) = \frac{1}{|\mathcal{K}|} \sum_{K \in \mathcal{K}} \frac{\partial \phi}{\partial x_l}(x_m), \quad l \in [1, d] \quad \text{and} \quad \psi_i(x_m) = \frac{1}{|\mathcal{K}|} \sum_{K \in \mathcal{K}} \phi(x_m), \qquad (4.45) $$

even though this set of constraints may yield unwanted behaviour of the polynomial within $C_i$ in case of distorted initial level set field, as direct point-values constraints are prone to yielding large derivative values and thus to creating spurious polynomial roots within the control volume.

When the underlying level set discretisation setting does not provide a pointwise reconstruction on $\partial C_i$ but rather face-averaged quantities, one can easily turn the above constraints into moments by testing their expressions against given functions and integrating over the control volume. On the same line, whenever the degrees of freedom of the initial level set discretisation are pointwise values or face-supported moments, one could also apply their expression directly to the locator function and use their values as constraints.
A note on high order representation and characterisation flexibility

Following the characterisation of the interface locator presented in the above paragraph when assuming the existence of a pointwise reconstruction of the initial level set field, allowing us to consider the pointwise quantities (4.44), the number of constraints that can be set up on each boundary is not re- restricted. The order of the polynomials $p_k^j$ is therefore up to one's choice, pro- viding an arbitrary high order ($\ge n_K$) representation of the interface that is similar to the one proposed in [192]. However, the conditioning of the charac- terisation system may worsen with the order increase, as it can be observed in the shape of the characterisation matrix shape derived in Section 4.3.4.2.

Another consequence of the chosen determination stencil is the lack of flex-
ibility in the polynomial representation, impacting the range of local features
that the interface locator (4.43) is able to represent. Indeed, as the actual
shape of the interpolation polynomial is directed by the face parametrisations
and constrained solely by boundary quantities, the interface locator would
smooth out local features that may have been initially detected by a high
order level set representation. In particular, any closed component of the in-
terface whose location is contained within the control volume and could be
inferred by a high order initial level set field would be completely ignored
or mishandled, as its existence would not be reflected in any boundary con-
straint and thus not directly known by any of the interpolation polynomials
$p_k^j$. Indeed, when not explicitly enforced, the local and closed zero level set of
the interface locator (4.43) is only a by-product of its variations, which may be
misplaced or even not existing if the local extrema are of the same sign. There-
fore, when a high order representation of the level set is available, one may
improve the definition of the locator function by refining the inner cell rep-
resentation, adding by example a polynomial that vanishes on the boundary
and whose dimension of the minimal space containing it matches the order
of the initial level set representation within the control volume. When this
refinement appears to be too costly, it is alternatively possible to perform the
interface correction on the subcells defining the variational grid, on which the
level set representation would be reduced to first order.

From a geometrical perspective, the given shape of the locator function is
very flexible towards the geometry of the control volume, and does not retain
the cell affiliation of the level set constraints in its expression. The initial
approximation of the interface location can thus be performed e.g. on non-
convex dual control volumes, even when the initial level set field is discretised
on convex primal elements. Therefore, as the corrector function (4.39) can
equally be defined for any control volume shape, one can straightforwardly
extend the correction procedure to any fully dual discretisation of the level
set and volume of fluid fields by considering (4.45) as the interface locator.
More challenging, let us now extend the procedure to the staggered setting.
### 4.3.3.3 Extensions to the staggered setting

Let us outline the correction procedure for staggered settings where the level set field is assumed to be reconstructed from simplicial primal elements while the volume of fluid constraints are given on the possibly non-convex dual cells, the reverse case being addressable in a similar way. As in the primal approach, the presented procedure can be set up regardless the schemes used for tackling (4.24) and (4.23). In particular, any discontinuous primal representation of the level set field is admissible, as well as any representation of the volume of fluid field from which an average value on the dual control cells can be determined.

**Notations** We consider here $\mathcal{J}^h$ the dual mesh of $\mathcal{T}^h$ and denote by $C_i$ its $i^{th}$ cell, centred at $x_i$. Any portion $C_i \cap K$ is called subcell. The boundary of $C_i$ is denoted $\partial C_i$, any of its face $f$. The set of faces crossed by the interface $\Lambda$ is further denoted $\Xi_i$ and any of those faces $j$ are referred to as “interfacial face”.

**Known information layout** We assume a local state and level set field to be defined on each primal element $K$, upon possibly a reconstruction from associated degrees of freedom values. We also assume $\bar{H}_i$, the average volume of fluid value, to be known or retrievable over $C_i$ without approximation error other than the one of (4.23)’s discretisation. The spatial layout of the known information is pictured in Figure 4.22b.

<figure><img src="image_78.png" /><figcaption>(a) Used spatial notations</figcaption></figure>

<figure><img src="image_79.png" /><figcaption>(b) Variational layout</figcaption></figure>

<caption>**Figure 4.22.** Spatially known information in a given staggered setting.</caption>

**Correction procedure** As presented in the fully primal approach, the level set values are updated at the end of the correction procedure that follows the guidelines of [20] and enforces conservation from the local volume of fluid constraint given for each dual cell $C_i$. Simultaneously, the interface continuity is obtained by the sole definition of the corrected locator function

$$ \psi: \Omega \to \mathbb{R} \\ x \mapsto \psi_i(x) + \psi_{i,c}^{\Delta}(x), \quad i \text{ s.t. } x \in C_i. \qquad (4.46) $$

More precisely, the continuity across each dual cell is enforced by the definition (4.40) of each $\psi_i$, and preserved using the corrector functions $\psi_i^{\Delta}$ ex-
pressed by (4.39), vanishing on $\partial C_i$. Paying attention to the interplay between the cells of different discretisations, the complete procedure takes eight steps.

1. Select the dual cell $C_i$ within which to construct the corrected interface and retrieve the set of primal elements $K$ that overlap with $C_i$, that is, define $\mathcal{K} = \{K \text{ s.t. } C_i \cap K \neq \emptyset\}$.

2. Retrieve the average volume of fluid value $\bar{H}_i$ on $C_i$.

3. Define a cross-cell continuous interface locator function $\psi_{i,c}$ on $C_i$ in the shape (4.40), and retrieve its parameters by setting constraints taken from the sets (4.44) and (4.45) based on the local level set reconstructions $\phi_i^K$ for all $K \in \mathcal{K}$. Alternatively, use directly the sets of degrees of freedom associated to the primal elements discretisation whose support are restricted to interfacial faces.

This third step can be simplified when the discretisation of (4.24) provides a continuous representation of the level set field over the computational domain by taking benefit of the staggered setting. Indeed, defining $\psi_i$ on each subcell $K \cap C_i, K \in \mathcal{K}$ as the primal reconstruction function over $K$ simply reads

$$ \psi_i: C_i \to \mathbb{R} \\ x \mapsto \phi^K(x), \quad K \text{ s.t. } x \in C_i \cap K. \tag{4.47} $$

There, the local reconstruction from the degrees of freedom to the continuous representation $\phi^K$ on each element $K$ contributing to $C_i$ brings directly the continuity across $\partial C_i$, each segment of the boundary being contained in a single primal element $K$ (see Figure 4.22b). Furthermore, the continuity within $C_i$ is equally immediate as the primal setting is continuous across primal elements. Hence, no additional determination constraint is required and one has just left to define a continuity-preserving corrector function.

4. Define a correction function $\psi_{i,c}^\Delta$ on $C_i$ in the shape (4.39), having a free parameter whose value impacts the zero level set of $\psi_i + \psi_{i,c}^\Delta$ and will be tuned against the volume of fluid constraint.

5. Set an artificial interface thickness $2\epsilon \in \mathbb{R}_+^*$.

6. Define the volume of fluid estimator $\mathcal{H}_i$ by applying the THINC reconstruction to the corrected locator function $\psi_i + \psi_{i,c}^\Delta$ and considering the interface thickness $2\epsilon$, that is,

$$ \mathcal{H}_{i,c} = \frac{1}{2} (1 + \tanh(\epsilon \psi_i + \epsilon \psi_{i,c}^\Delta)). \tag{4.48} $$
7. Set the conservative constraint on the cell $C_i$ and determine the pa- rameter of the corrector function $\psi_{i,c}^{\Delta}$ accordingly, that is,

$$ \text{Find } c \text{ s.t. } \frac{1}{|C_i|} \int_{C_i} \mathcal{H}_{i,c}(x) dx = \bar{H}_i. \qquad (4.49) $$

Note that as the properties of $\psi_{i,c}^{\Delta}$ are independent of the shape of the dual cell, the well-posedness of this problem is identical to the fully primal setting and will not be specifically discussed. However, the existence of a solution to (4.49) is here a-priori not enough to guarantee the interface continuity. Indeed, one needs to ensure that the solution found on the dual cell transfers to the primal setting, property ensured only by the last step.

8. Redefine the level set function on the primal elements $K$ from the dual reciprocal volume of fluid estimator

$$ \begin{aligned} \forall x \in K \cap C_i, \quad \phi(x) &= \frac{1}{\epsilon} \tanh^{-1}(2\mathcal{H}_{i,c}(x) - 1) \qquad (4.50) \\ &= \psi_i(x) + \psi_{i,c}^{\Delta}(x). \end{aligned} $$

By definition of the functions $\psi_i$ and $\psi_{i,c}^{\Delta}$, the interface is continuous within and across any dual cell $C$, as depicted in Figure 4.23a. Therefore, as the corrected level set values (4.50) are directly built on each subcell from the corrected locator function $\psi_i + \psi_{i,c}^{\Delta}$, the interface stays continuous and conservative when considered from the primal representation (see Figure 4.23b for an illustration).

<figure><img src="image_80.png" /><figcaption>Figure 4.23. Correction procedure and evaluation of the level set values.</figcaption></figure>
### 4.3.3.4 Extensions to general forms of interface functions

Even though we could already extend the proposed corrector and locator functions to more general discretisation settings as discontinuous or staggered, their proposed shapes remain intuitive and are not able to take into account specific interface behaviours within the inner cell, as e.g. those resulting from a local fluid state or obeying to a physical principle such as surface tension. We therefore aim to further generalise the approach by deriving sufficient conditions on the shape of the interface locator and corrector functions to yield a continuous and conservative corrected interface by construction.

**Admissibility conditions on locator functions** The correction procedure requires a continuous interface locator function whose expression can be modified locally on each control volume of the volume of fluid discretisation. Thus, whenever the discretisation of the level set field is not performed on the same mesh, is discontinuous, or when a higher representation order is aimed at, a cross-cell continuous local reconstruction of the interface locator is required.

In order to derive admissible definitions of reconstruction functions, we consider interfacial control volumes {$C_i$}$_i$ either belonging to $\mathcal{T}^h$ or $\mathcal{J}^h$ within which to design a local locator function $\psi_i$, and aim to approximate continuously the possibly non-conservative initial interface as its zero level set. This goal in mind, one should therefore define those local functions so that their sum defines over $\Omega$ a continuous implicit interface representation of the form

$$ \psi: \Omega \to \mathbb{R} \\ x \mapsto \psi_i(x), \quad i \text{ s.t. } x \in C_i. $$

Furthermore, in view of the correction procedure, one may consider locator functions $\psi_i$ that are exactly supported on $C_i$ in order to allow a compact correction and prevent multiply defined values. Asking $\psi$ to determine implicitly a continuous interface would there reduce to asking the cross-cell and within continuities of the local interface represented in each cell $C_i$ by the local locator function $\psi_i$. Hence, one should shape the functions {$\psi_i$}$_i$ so that the boundary crossing locations of its zero level set on a given face coincide with those determined by the neighbouring local reconstructions.

Lastly, in addition to being continuous, the interface should delimit the actual subdomains of the problem (4.1–4.2). Each local reconstruction should therefore be constrained by the initial level set field $\phi$, the location of the interface’s boundary crossings approximating its zeros in order to preserve the interfacial nature of the selected faces under the reconstructed representation $\psi$. One can therefore derive sufficient conditions over the interface locator shape by asking it to vanish on the boundary at the same locations as $\phi$. All in all, the set of admissibility conditions then reads as follows.
Admissibility conditions 4.3.1 Locator function

A function $\psi_i$ can be used as an interface locator function on a cell $C_i$ when

* $\psi_i$ is supported on $C_i$,

* $\psi_i$ is continuous within $C_i$ and $\psi$ is continuous across $\Xi_i$,

* the zeros' location of $\psi_i$ on $\partial C_i$ coincide with those of the interface representations $\phi^K$ on $C_i \cap K$ for each element $K$ of the initial level set discretisation that overlaps with $C_i$.

Additionally, one may ask the construction of $\psi_i$ to stay compact, that is,

* $\psi_i$ is determined from level set values contained in or on $C_i$.

The first and second properties ensure that the expression of $\psi_i$ is compact and continuous within and across the control volume boundaries. The third property guarantees that the interface represented by $\psi_i$ coincides on $\partial C_i$ with the local representations of the initial level set field whose support intersect with $C_i$, where the term coincide has to be understood in the error-minimisation sense whenever the initial level set field $\phi$ is discontinuous and yields different interface locations on the boundary. In practice, one could consider e.g. finding the interface locator $\psi_i$ whose zero level set minimises the minimum homotopy area between the different interface locations.

Having then defined an interface locator $\psi_i$ fulfilling Admissibility conditions 4.3.1, its zero level set naturally provides an admissible interface contained within the cell $C_i$, interpolating on its boundary the one determined by the non-conservative initial level set field. Within the control volume however, the Admissibility conditions 4.3.1 do not require the interface locator to be specifically constrained, as its inner behaviour does not prevent to obtain a cross-cell continuous and conservative interface trough the correction procedure presented in Section 4.3.2, as already observed in Section 4.3.3.3. One is therefore free to design its inner shape prior to its correction by using any continuous function supported on $C_i$, aiming e.g. to approximate the initial zero level set as accurately as possible, to enforce a given representation order or even to consider physical principles such as surface tension. However, as those tunings necessarily increase the dimension of the space where the locator function live, one has to design enough independent determination constraints in order to guarantee unisolvence, or be keen on retrieving a rough initial approximation through an over or underdetermined system. In addition, one should prevent the creation of spurious interfaces by a misdistribution of constraints, possibly yielding to oscillating interpolation functions whenever high order polynomials or harmonic functions are involved.
**Note** For simplifying the derivation of the interface locator expression $\psi$, one can strengthen the Admissibility conditions 4.3.1 by asking continuity of the local functions $\psi_i$ themselves instead of their zero level set.

**Admissibility conditions on corrector function** Now that we could determine sufficient conditions for the locator function $\psi_i$ to be cross-cell continuous, let us derive a set of conditions for the corrector function $\psi_i^\Delta$ to be admissible. According to the correction procedure, it is enough for $\psi_i^\Delta$ to be a parameter-dependent function that has enough freedom for deforming the zero level of $\psi_i + \psi_i^\Delta$ towards the conservation constraint (4.49) while preserving the continuity of $\psi$. Considered together, those requirements reduce to the following Admissibility conditions 4.3.1.

### Admissibility conditions 4.3.2 Corrector function

A function $\psi_i^\Delta$ can be used as an interface corrector on a cell $C_i$ when the following properties hold.

* $\psi_i^\Delta$ is supported on $C_i$

* $\psi_i^\Delta$ is continuous within $C_i$ and preserves the interface location of $\psi$ across any interfacial face $j \in \Xi_i$

* $\psi_i^\Delta$ has at least one free parameter

* The parameters range of $\psi_i^\Delta$ allows the zero level set of $\psi_i + \psi_i^\Delta$ to define local subdomains given any area proportion constraint in (0, 1)

Impacting directly the correction process, the two last properties ensure the feasibility of (4.49) and therefore the existence of a corrector function for any value of the volume of fluid constraint. Note however that here again no requirement is set on the shape of the interface within the control volume other than respecting the interface boundary crossing points and splitting the control volume so that the volume of fluid constraint is matched. It can therefore be freely modified within $C$ in order to meet the conservation constraint, in particular without asking its implicit representation to yield uniqueness of the solution to the problem (4.49). Even though this freedom amount may seem unrealistic, there do exist several phase distributions within $C_i$ that matches the same volume of fluid constraint $\bar{H}_i$ (see Figure 4.24 for an illustration). Narrowing down the possible spatial layouts comes with enforcing physical properties or inner level set constraints, either directly in the design of the implicit interface corrector function or by imposing further determination constraints within the cell.
More importantly, on the top of preserving any pre-existing continuity of $\psi$ over the computational domain $\Omega$, the two first constraints ensure that the discrete conservation properties enforced on each control volume associated to the volume of fluid discretisation transfer to the possibly different mesh of the initial level set discretisation. Indeed, the last step of the correction procedure updates the level set values regardless of its original discretisation, rebuilding the field through the evaluation of the continuous interface locator function (4.26) that is here considered as a reference. The updated values of the level set thus become grid-independent, transferring the discrete conservation principle to any partition of $\Omega$.

<figure><img src="image_81.png" /><figcaption>**Figure 4.24.** Several interface distortions defining cells partitions that yield a same volume fraction in a given convex control volume.</figcaption></figure>

**A short note on further extensions** Encompassing a wide range of expressions both for the locator and corrector functions, the Admissibility conditions 4.3.1 and Admissibility conditions 4.3.2 open a lot of possibilities for further extensions, paving in particular the way to refining the interface location within the cells. In particular, one can improve the physical relevance of the interface by adding inner constraints, and achieve higher interface regularity by dropping the optional requirement for the locator and corrector functions to be determined compactly within the cell. Generalisation to higher dimensions are also possible, being especially straightforward when considering the proposed interface locator (4.40) and corrector (4.39) functions.

The continuous and conservative correction approach being now well defined, let us specify the numerical approximations required to turn the theoretical procedure into a practical technique.
### 4.3.4 Numerical approximations

The practical definition of the presented correction procedure comes with a numerical method that determines the prescribed coefficients of $\psi_i$, the definition of a quadrature scheme on the control volumes $C_i$ for discretising the conservation problem (4.29), and a numerical technique to solve it. The following subsections detail a possible numerical method for each of those aspects.

#### 4.3.4.1 Quadrature rule over any control volume

Let us first define a quadrature rule over any control volume $C_i$ by taking advantage of its polygonal shape. Its convexity not being necessarily known a-priori, we start by tessellating the control volume by simplices, applying then any classical quadrature rule on each of the obtained simplicial subcells.

In practice, any control volume $C_i$ can be seen as a mesh element and its shape can be assimilated to a simple polygon. Geometrically speaking, one can therefore tile $C_i$ with simplicial subdomains simply by using a Delaunay triangulation [195]. Gaussian quadrature points are then typically considered in each subcell, and the quadrature scheme on $C_i$ can be expressed for any integrable integrand $f: C_i \to \mathbb{R}$ by

$$ \frac{1}{|C_i|} \int_{C_i} f(x) dx = \sum_{T \in C_i} \sum_{g=1}^{G} w_g f(x_g) + o(h^{G+1}), $$

where $\{(x_g, w_g)\}_{g \in [1, G]}$ is the set of points and weights of the degree $G \in \mathbb{N}$ Gaussian quadrature rule associated to the triangle $T$ of the simplicial tessellation (see e.g a

possible quadrature points distribution in Figure 4.25). When the computational power is not an issue, one can significantly improve the overall quadrature points distributions in $C_i$ by optimising the aspect ratio of each subcell, using e.g. a ray shooting algorithm [196] in place of the Delaunay triangulation.

<figure><img src="image_82.png" /><figcaption>Figure 4.25. Quadrature points distribution and impact of the geometrical regularity of the subcells.</figcaption></figure>

#### 4.3.4.2 Determination of the locator function coefficients

Let us now focus on the determination of the initial interface locator's coefficients under a set of constraints derived from the level set field. Even though the presented technique holds for any linear setting, we focus for the sake of the explanation on the particular interface locator (4.40) built from the face parametrisation (4.41) on a control volume $C_i$, and characterised by the sets of level set constraints (4.44) and (4.45).
In this context, we denote by $m = |\Xi_i|$ the number of interfacial faces $j \in \Xi_i$ and consider $m$ polynomials $\{p_k^j \in \mathbb{Q}_k(\mathbb{R})\}_{j \in [1,m]}$ supported on $C_i$. Defined by the expressions (4.40) and (4.41), the initial locator function $\psi_i$ then reads

$$ \psi_i: C_i \to \mathbb{R} \qquad T_j: C_i \to \mathbb{R} \\ x \mapsto \sum_{j \in \Xi_i} (p_k^j(\gamma_j(x))T_j(x)), \quad \text{where} \quad x \mapsto \prod_{\substack{i \in \Xi_i \\ i \neq j}} \Gamma_i(x). \tag{4.51} $$

There, as detailed in Section 4.3.3.2, $\gamma_j = \frac{\|x-x_A\|}{\|x_B-x_A\|}$ and $\Gamma_i = d_j x_2 - a_i x_1 - b_i$ represent face parametrisations, where $x_A, x_B$ are the end points of the face $j \in \Xi_i$, and $d_j, a_j$ are the leading and $b_j$ the offset coefficients of the line extending the face $j \in \Xi_i$. So defined, the coefficients $\{d_j, a_j, b_j\}_{j \in [1,m]}$ are known by the geometrical definition of $C_i$ and no freedom is given to the expression of $\Gamma_i$. One has then to determine $m(k+1)$ unknowns, each of them being only involved in the expression of a given polynomial $p \in \{p_k^j\}_j$.

In order to specify those unknowns with respect to the constraints (4.44) and (4.45), a straightforward Lagrangian interpolation may be used. However, as the conservation is anyhow retrieved afterwards, there is no need to interpolate exactly the level set values or its derivatives. Therefore, we rather consider as in [20] a least square approximation method, flexible towards the number of constraints applied to the initial locator function, and derive here the associated determination system.

**Derivatives expressions** Let us develop the expression of the derivatives encountered in the constraints (4.44) and (4.45). Starting by $T_j$, and denoting indiscriminately $a_r = a_j, b_r = b_j, d_r = d_j$ the coefficients of the Cartesian equation of the line extending the face $j \in \Xi_i$ indexed by $r \in \mathbb{N}$, we retrieve

$$ \begin{aligned} T_j &= \prod_{\substack{i \in \Xi_i \\ i \neq j}} (d_j x_2 - \underbrace{(a_i x_1 + b_i)}_{:= z_i}) \\ &= \sum_{k=1}^{m+1} (-1)^{m-k+1} \left(\prod_{r=1}^{k} d_r\right) x_2^k \sum_{\substack{0 \le t_1 < t_2 < \dots \\ t_{m-k} \le t_{m-k+1}}} \left(\prod_{j=1}^{m-k+1} z_{t_j}\right) \\ &= \sum_{k=1}^{m+1} (-1)^{m-k+1} \left(\prod_{r=1}^{k} d_r\right) x_2^k \sum_{\substack{0 \le t_1 < t_2 < \dots \\ t_{m-k} \le t_{m-k+1}}} \left(\sum_{p=1}^{m-k+2} \left(\prod_{r=1}^{p} a_r\right) x_1^p \sum_{\substack{0 \le l_1 < l_2 < \dots \\ l_{m+1-k-p} \le l_{m-k-p+2}}} \left(\prod_{q=1}^{m-k-p+2} b_{l_{l_q}}\right)\right), \end{aligned} $$

allowing us to specify the coefficients $T_j^{k,p}$ attached to each monomial

$$ T_j^{k,p} := \sum_{\substack{0 \le t_1 < t_2 < \dots < t_{m-k} < l_{m+1-k-p} \\ l_{m+1-k} \le t_{m+1} \le l_{m+2-k-p}}} \sum_{r=1}^{k} \left(\prod_{r=1}^{k} d_r\right) \left(\prod_{r=1}^{p} a_r\right) \left(\prod_{q=1}^{m+2-k-p} b_{l_{l_q}}\right), \tag{4.52} $$
and to express the expression (4.53) as

$$T_j = \sum_{k=1}^{m+1} (-1)^{m+1-k} \sum_{p=1}^{m+2-k} T_f^{k,p} x_1^p x_2^k. \quad (4.53)$$

The interface locator $\psi_i$ thus reduces to

$$\psi_i(x) = \sum_{j \in \Xi_i} [p_j^l(\gamma_i(x)) T_j(x)], \quad (4.54)$$

and knowing that $p_k^l(\gamma_f(x)) = \sum_{k=1} A_k^l(\gamma_f(x))^k$, for $A_k^l \in \mathbb{R}$, we get

$$\psi_i(x) = \sum_{j \in \Xi_i} \sum_{k=1} A_k^j \underbrace{(\gamma_f(x))^k T_j(x)}_{:=C_j^k}. \quad (4.55)$$

The partial derivative of the interface locator in the $x_1$ direction then reads

$$\frac{\partial \psi}{\partial x_1}(x) = \sum_{j \in \Xi_i} \sum_{k=1} A_k^j \underbrace{\left( k \gamma_j^{k-1}(x) \frac{\partial \gamma_j}{\partial x_1}(x) T_j(x) + k \gamma_j^k(x) \frac{\partial T_j}{\partial x_1}(x) \right)}_{:=D_{x_1} C_j^k} \quad (4.56)$$

its second order partial derivative in the $x_1$ direction thus writing

$$\frac{\partial^2 \psi}{\partial x_1^2}(x) = \sum_{j \in \Xi_i} \sum_{k=1} A_k^j \left( k T_j(x) (\gamma - 1)^{k-2}(x) \left(\frac{\partial \gamma_j}{\partial x_1}\right)^2(x) + k \gamma_j^{k-1}(x) \frac{\partial^2 \gamma_j}{\partial x_1^2}(x) \right. \\ \left. + 2k \frac{\partial T_j}{\partial x_1}(x) \gamma_j^{k-1}(x) \frac{\partial \gamma_j}{\partial x_1}(x) + \gamma_j^k(x) \left(\frac{\partial T_j}{\partial x_1}\right)^2(x) \right). \quad (4.57)$$

The first and second order partial derivative along $x_2$ are given analogously, while the second order partial derivative along $x_1x_2$ reads however as follows.

$$\frac{\partial^2 \psi}{\partial x_2 \partial x_1}(x) = \sum_{j \in \Xi_i} \sum_{k=1} A_k^j \left( T_j(x) \left( k(k-1)\gamma_j^{k-2}(x) \frac{\partial \gamma_j}{\partial x_2}(x) \frac{\partial \gamma_j}{\partial x_1}(x) + k^2\gamma_j^{k-1}(x) \frac{\partial \gamma_j}{\partial x_2 \partial x_1}(x) \right) \right. \\ \left. + \frac{\partial T_j}{\partial x_2}(x) T_j(x) k \frac{\partial \gamma_j}{\partial x_1}(x) \gamma_j^{k-1}(x) \right. \\ \left. + \frac{\partial T_j}{\partial x_1}(x) k \frac{\partial \gamma_j}{\partial x_2}(x) \gamma_j^{k-1}(x) \right. \\ \left. + \underbrace{\frac{\partial^2 T_j}{\partial x_2 \partial x_1}(x) \gamma_j^k(x)}_{:=D_{x_2 x_1} C_j^k} \right). \quad (4.58)$$

There is then just left to express the derivatives of the face parametrisations $\Gamma_i$ and $\gamma_j$, given for reference in Table 4.1 without further computation details.
<table><caption>**Table 4.1.** Derivatives of the face parametrisations.</caption><thead><tr><th colspan="2">Derivatives of the face-selecting polynomial T<sub>j</sub></th></tr></thead><tbody><tr><td>First order</td><td><div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding: 0 5px;">∂T<sub>j</sub> /∂x<sub>1</sub>(x) = &sum;<sub>k=1</sub><sup>m+1</sup> (-1)<sup>m+1-k</sup> &sum;<sub>p=1</sub><sup>m+2-k</sup> p T<sub>j</sub><sup>k,p</sup> x<sub>1</sub><sup>p-1</sup> x<sub>2</sub><sup>k</sup></div></td></tr><tr><td>Second order</td><td><div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding: 0 5px;">∂<sup>2</sup>T<sub>j</sub> /∂x<sub>1</sub><sup>2</sup>(x) = &sum;<sub>k=1</sub><sup>m+1</sup> (-1)<sup>m+1-k</sup> &sum;<sub>p=1</sub><sup>m+2-k</sup> p(p-1) T<sub>j</sub><sup>k,p</sup> x<sub>1</sub><sup>p-2</sup> x<sub>2</sub><sup>k</sup></div></td></tr><tr><td></td><td><div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding: 0 5px;">∂<sup>2</sup>T<sub>j</sub> /∂x<sub>2</sub>x<sub>1</sub>(x) = &sum;<sub>k=1</sub><sup>m+1</sup> (-1)<sup>m+1-k</sup> &sum;<sub>p=1</sub><sup>m+2-k</sup> pk T<sub>j</sub><sup>k,p</sup> x<sub>1</sub><sup>p-1</sup> x<sub>2</sub><sup>k-1</sup></div></td></tr><tr><th colspan="2">Derivatives of the face parametrisation γ<sub>j</sub></th></tr><tr><td>First order</td><td><div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding: 0 5px;">∂γ<sub>j</sub> /∂x<sub>1</sub>(x) = <div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding: 0 5px;">(x<sub>1</sub>-x<sub>A</sub>)</div></div> / ||x-x<sub>A</sub>||||x<sub>B</sub>-x<sub>A</sub>||</div></td></tr><tr><td>Second order</td><td><div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding: 0 5px;">∂<sup>2</sup>γ<sub>j</sub> /∂x<sub>1</sub><sup>2</sup>(x) = <div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding: 0 5px;">||x-x<sub>A</sub>||<sup>2</sup>-x<sub>1</sub>(x<sub>1</sub>-x<sub>A</sub>)</div></div> / ||x-x<sub>A</sub>||<sup>3</sup>||x<sub>B</sub>-x<sub>A</sub>||</div></td></tr><tr><td></td><td><div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding: 0 5px;">∂<sup>2</sup>γ<sub>j</sub> /∂x<sub>1</sub>x<sub>2</sub>(x) = <div style="display: inline-block; vertical-align: middle; text-align: center;"><div style="border-bottom: 1px solid black; padding: 0 5px;">-(x<sub>2</sub>-x<sub>A</sub>)(x<sub>1</sub>-x<sub>A</sub>)</div></div> / ||x-x<sub>A</sub>||<sup>3</sup>||x<sub>B</sub>-x<sub>A</sub>||</div></td></tr></tbody></table>

**Derivation of the determination system** Using the expressions (4.56), (4.57) and (4.58), let us now express the constraints selected within the sets (4.44) and (4.45) over the expanded expression of $\psi$ (4.53) as a linear system. To this aim, we assume that the constraints are set on a finite number $M$ of sampling points $\{x_m\}_{m \in [1, M]}$ within $C_i$, and that the value of their right hand side can be directly obtained from any local reconstruction of the level set field $\phi$. In such a setting, the determination system over $\psi_i$ on the cell $C$ then reads

$$XA = \Psi, \qquad (4.59)$$

where $X$ is the matrix containing the evaluated monomials of (4.53) at the sampling points $\{x_m\}_m$, writing

$$X = \begin{pmatrix}
\begin{array}{c}
\text{Point value constraints} \\
\begin{array}{c}
C_1^1(x_1) \quad \dots \quad C_1^k(x_1) \\
\vdots \\
C_1^1(x_M) \quad \dots \quad C_1^k(x_M)
\end{array}
\end{array} & C_2^1(x_1) \quad \dots \quad C_2^k(x_1) \quad \dots \\
\begin{array}{c}
\text{Quantities associated to the first edge, evaluated at } \{x_i, i \in [1, M]\} \\
\begin{array}{c}
C_1^1(x_M) \quad \dots \quad C_1^k(x_M)
\end{array}
\end{array} & C_2^1(x_M) \quad \dots \quad C_2^k(x_M) \quad \dots \\
\begin{array}{c}
\text{Point value considering first order derivatives} \\
\begin{array}{c}
D_{x_1} C_1^1(x_i) \quad \dots \quad D_{x_1} C_1^k(x_i) \\
\vdots \\
D_{x_2} C_1^1(x_i) \quad \dots \quad D_{x_2} C_1^k(x_i)
\end{array}
\end{array} & D_{x_1} C_2^1(x_i) \quad \dots \quad D_{x_1} C_2^k(x_i) \quad \dots \\
\begin{array}{c}
\vdots \\
\vdots \\
\vdots
\end{array} & D_{x_2} C_2^1(x_i) \quad \dots \quad D_{x_2} C_2^k(x_i) \quad \dots \\
\begin{array}{c}
D_{x_1 x_1} C_1^1(x_i) \quad \dots \quad D_{x_1 x_1} C_1^k(x_i) \\
\vdots \\
D_{x_2 x_2} C_1^1(x_i) \quad \dots \quad D_{x_2 x_2} C_1^k(x_i)
\end{array} & D_{x_2 x_2} C_2^1(x_i) \quad \dots \quad D_{x_2 x_2} C_2^k(x_i) \quad \dots \\
\begin{array}{c}
\vdots \\
\vdots \\
\vdots
\end{array} & D_{x_1 x_2} C_2^1(x_i) \quad \dots \quad D_{x_1 x_2} C_2^k(x_i) \quad \dots \\
\begin{array}{c}
\vdots \\
\vdots \\
\vdots
\end{array} & D_{x_1 x_2} C_1^1(x_i) \quad \dots \quad D_{x_1 x_2} C_1^k(x_i)
\end{array}
\end{pmatrix}, \quad (4.60)$$
and A is the vector of coefficients of each polynomial $p_k^j$ given by

$$A = \left( A_1^1 \cdots A_1^k \left[ A_2^1 \cdots A_2^k \cdots \cdots A_m^1 \cdots A_m^k \right] \right)^T \quad (4.61)$$

On the right-hand-side, $\Psi$ represents the level set values or its derivative at the evaluated points and reads

$$\Psi = \left( \phi(x_1) \cdots \phi(x_M) \cdots \frac{\partial \phi}{\partial x_1}(x_i) \cdots \frac{\partial \phi}{\partial x_2}(x_i) \cdots \frac{\partial \phi}{\partial x_1^2}(x_i) \cdots \frac{\partial \phi}{\partial x_2^2}(x_i) \cdots \frac{\partial \phi}{\partial x_1 x_2}(x_i) \cdots \right)^T \quad (4.62)$$

So defined, the system (4.59) can be either under or over-determined depending on the number of selected constraints and wished polynomial order on each interfacial face. In those cases, a Least-Square approximation method may be used to retrieve the vector A, providing a rough locator function $\psi_i$ in the form of (4.55).

**A note on the constraint selection** As the solution of the system (4.59) is to be retrieved with a Least-Square approach, one does not explicitly require the matrix X to be invertible and is therefore free to choose the type of constraints to apply on the expression of $\psi_i$, even though prescribing constraints regardless of unisolvence considerations may significantly reduce the accuracy of the initial interface approximation. In the proposed setting however, the unisolvence of (4.59) is guaranteed by the definitions of the interface locator (4.40) and constraints expressions (4.44) and (4.45). Indeed, as derived in the first paragraph of this section, the left hand side of the constraints expressions applied to (4.40) at a given point $x_m \in C_i$ are linearly independent. As a consequence, whenever as many constraints as unknowns are selected, the solution of (4.59) is an actual polynomial interpolation, not only a best polynomial fit.

Regarding the constraint selection, one can first notice that the number of constraints that are applicable to the interface locator is unbounded, as the sets (4.44) and (4.45) define pointwise constraints. It is therefore enough to select as many points $x_m$ as needed in order to cope with the dimension of the space the interface locator is discretised in, knowing that the set of constraints (4.44) and (4.45) provide several constraints per discrete point (e.g. up to three in two dimensions). However, the number and spatial distribution of the selected points impact the conditioning number of the determination matrix, making their choice determinant for an accurate representation of the initial interface. On a last note, let us point out that those observations transfer to moment-type constraints, where the point selection have to be replaced by the definition of a set of test functions.
### 4.3.4.3 Numerical techniques to enforce conservation

Having approximated the quadrature scheme and determined the initial locator function, let us focus on a numerical technique tackling (4.49) in the particular case of the interface locator and corrector shapes introduced in Section 4.3.2. In essence, we aim to derive a numerical procedure that approximates the solution $c \in \mathbb{R}$ to the problem

$$ \bar{H}_i := \frac{1}{|C_i|} \int_{C_i} \frac{1}{2} (1 + \tanh(\epsilon\psi_i + \epsilon\psi_{i,c}^{\Delta})) dx = \bar{H}_i, \quad (4.49) $$

where the interface locator $\psi_i$ and corrector $\psi_{i,c}^{\Delta}$ are respectively given by the expressions (4.40) and (4.39), the constraint $\bar{H}_i$ being directly retrieved on the control volume $C_i$ from a conservative approximation of (4.23).

**Problem discretisation**

Considering $C_i$ as any simple polygon, we start with the discretisation of the left-hand side by directly making use of the quadrature rule derived in Section 4.3.4.1. More precisely, given the set of quadrature points and weights $\{(x_g, w_g)\}_{g \in [1, G]}$ distributed on the control volume $C_i$, the problem (4.49) approximates to

$$ \bar{H}_i(c) := \frac{1}{2|C_i|} \sum_{g=1}^{G} w_g (1 + \tanh(\epsilon\psi_i(x_g) + \epsilon\psi_{i,c}^{\Delta}(x_g))) = \bar{H}_i. \quad (4.63) $$

There, the expressions of $\psi$ and $\psi^{\Delta}$ are continuous and explicit within $C_i$, allowing their direct evaluation. The right hand side term $\bar{H}_i$ being a predefined constant, the discretisation of (4.49) thus immediately follows.

**Numerical approximation of a solution**

As seen in Section 4.3.2.3, the shape of the corrector function allows to correct any initial interface locator by providing a unique value of the correction parameter $c$ for which the conservative constraint $\bar{H}_i$ is achieved, without imposing any restriction on the shape of the initial locator function other than its continuity. In addition, as defined in (4.63), the approximated volume of fluid estimator function $\bar{H}_i$ is valued in $(0, 1)$ and, up to quadrature errors, is strictly monotonous on $\mathbb{R}$. Any root finding method applying on open intervals is therefore suitable, and we choose here the simple secant method for the sake of generality, not relying on any knowledge or approximation of $\bar{H}_i$'s derivative. Namely, we iterate

$$ \begin{cases} c_0 = c \\ c_1 = -c \\ c_{n+1} = c_n - \bar{H}_i(c_n) \frac{c_n - c_{n-1}}{\bar{H}_i(c_n) - \bar{H}_i(c_{n-1})} \end{cases} $$

over $n \in \mathbb{N}$ until the tolerance $|c_{n+1} - c_n| < 1e-8$ is reached from an initial guess $c \in \mathbb{R}$, chosen as to avoid initial computational instabilities and round off errors in the right hand side of the iteration rule.
**A note on numerical defects** Relying on several approximations, the practical realisation of the discrete conservation principle through the correction procedure remains up to the accuracy of involved numerical methods and the amount of computational errors.

In particular, as defined in [20] and adapted in Section 4.3.2, the correction procedure is sensitive to the quadrature scheme, both in terms of its accuracy and layout of its sampling points. Especially, as already mentioned in Section 4.1.2 and depicted in Figure 4.9, the discrete volume of fluid estimator $\bar{H}_i$ may yield different values for cells within which the local subdomains occupy a same volume but have a different spatial distribution, the layout of their interface not necessarily encompassing weight-equivalent quadrature points. The distribution of the quadrature points $x_g$ may therefore prevent the continuous deformation of the interface from being accurately described, in turn impeding the value $\bar{H}_i$ to adapt smoothly with the correction amount. In particular, the approximated function $\bar{H}_i$ may lose the strict monotonicity of $\bar{H}_i$ in case of subdomains expanding or deforming along with $c$, as depicted in Figure 4.26. In extreme cases, the numerical method used to tackle (4.63) may be impacted, as e.g. when a local extremum is created or when the approximated volume of fluid estimator function $\bar{H}_i$ adopts a constant behaviour in a vicinity of its zero.

<figure><img src="image_83.png" /><figcaption>**Figure 4.26.** Expanding (left) or distorting (right) interfaces for which the continuous expansion of the generated subdomain of reference is not represented through the value of the approximated volume of fluid estimator $\bar{H}_i$.</figcaption></figure>

Another delicate point when considering the numerical application of the presented technique relates to the artificial width of the interface, defined as the steepness parameter of the volume of fluid estimator. Indeed, its value should never be too high for the root finding algorithm to operate, but still be high enough for the volume of fluid evaluation to remain meaningful. If this concern has already been taken into account in the classical *THINC-LS* approach [20], one can here be more tactful and consider one steepness parameter for each of the two functions, corrector and locator, each parameter depending on the function's general shape and the natural steepness of its variations around the interface.

Having described all the involved numerical approximations, let us now assess the versatility of the presented correction approach as well as its relevance in the numerical results Section 4.3.5.
4.3.5 Numerical results

Well established theoretically, let us assess the practical relevance of the pre-
sented approach by testing the corrector function ψ Δ i,c as expressed by (4.3.3.1)
for various control volume shapes. We start in Section 4.3.5.2 by investigating
the spatial layout of corrected interfaces, before exemplifying the Proposition
4.3.1 in Section 4.3.5.3 by testing the behaviour of the evaluated volume of
fluid quantity H̅ᵢ against the correction parameter. Lastly, we explore the lim-
itations of the proposed approach in Section 4.3.5.4, investigating the impact
that highly non-linear locator functions have on the shape of the corrected
interface, and focusing on huge correction amounts.

**4.3.5.1 Implementation of involved numerical techniques**

As already detailed in Section 4.3.4, the correction procedure requires a collection of numerical techniques to be carried out in practice. Here, the presented results were obtained by using the composite quadrature rule presented in Section 4.3.4.1 for computing the evaluated volume of fluid Hᵢ(c), and by using the secant method for finding the correction parameter, as presented in Section 4.3.4.3.

Furthermore, all the results involving non convex control volume are based
on the smooth bijective map described in Section 4.3.3.1, using as boundary
term the collection of the *n* affine transformations mapping the *i*<sup>th</sup> face of the
control volume to the *i*<sup>th</sup> face of the reference polygonal shape by designating
the plane origin as its invariant. The two Poisson problems have there been
practically solved by a Fenics [197] solver using a second order Finite Element
Continuous Galerkin method on a simplicial mesh of the control volume.

<figure><img src="image_84.png" /><figcaption>**Figure 4.27.** Smooth bijective map between a non-convex octagon (left) and its reference shape (right), showing the inner distortions.</figcaption></figure>

Lastly, note that the initial locator function satisfying |∇ψᵢ| ≡ 1 have been ob-
tained by redistancing a possibly highly non-linear function that generates
the zero level set of interest. Already used right after the level set evolution
when considering a two-phase flow problem, the practical implementation is
borrowed from the *LESARD* code [198] and involves a Gradient descent pro-
cedure applied to the Lax-Hopf formulation of the problem (4.18).
4.3.5.2 Spatial behaviour of the interface

Let us start by investigating the dependence of the corrected interface shape
on the nature of the control volume and the regularity of the initial interface
locator. We represented in Figure 4.28 the considered control volume, the zero
level of the initial interface locator as a plain line and the corrected interfaces
in dashed line for several values of the correction parameter *c*.

<figure><img src="image_85.png" /><figcaption>(a) Corrected interfaces for an initially linear locator function</figcaption></figure>

<figure><img src="image_86.png" /><figcaption>(b) Corrected interfaces for an initially quintic locator function</figcaption></figure>

<figure><img src="image_87.png" /><figcaption>(c) Corrected interfaces for an initially linear locator function not having zeros within the control volume nor on its boundary</figcaption></figure>

**Figure 4.28.** Representation of the initial interface (plain line) and corrected interfaces (dashed lines) for various values of the correction parameter *c* and control volume shapes.

First of all, we can observe that as stated in Corollary 4.3.1, the corrected
interface is contained within the control volume, not modifying the number of
interface boundary crossings or changing their locations. Valid for any shape
of the initial locator function, we further observe that the spatial deformation
of the initial interface is smooth with respect to the correction parameter.
In light of those behaviours, one can further wonder whether topological
changes of an already existing interface can occur through its correction, as it
already appears in the bottom row of Figure 4.28 that a new interface can be
forced into the control volume when no interface pre-existed. We thus investi-
gate the particular behaviour of multi-component interfaces, whose test cases
are reported in Figure 4.29.

<figure><img src="image_88.png" /><figcaption>**Figure 4.29.** Layout of an initially multi-component interface depending on the correction parameter *c*. Top: Degree 6 locator function having an initial zero level set made of three distinct connected components. Bottom: Degree 5 locator function having an initial zero level set made of three distinct connected components, two contained within the control volume and one crossing it.</figcaption></figure>

There, we are indeed facing topological changes of the interface whenever the
correction amount becomes extensive. Especially, in both presented cases the
originally disconnected subdomains merge as the evaluated volume of fluid
value gets high. Reversely, when its value decreases, connected components
that are contained within the control volume simply vanish, while crossing
interfaces are preserved and tend to run alongside the boundary.

As those topological changes depend on the space-dependent relationship
between the variations of the locator and corrector, they may seem dramati-
cal. Yet, the changes in the spatial layout of the corrected interface remain
smooth with respect to the correction parameter, and the subdomains propor-
tions are monotonously adjusting. Therefore, the topological changes can be
seen as a simple by-product of the problem (4.49)'s well-posedness when fac-
ing an oscillating locator function whose variations locally overgrows those
of the corrector function. Indeed, as one is able to reach any volume of fluid
constraint and thus any subdomains' proportion through the correction proce-
dure, the deformation, merge or creation of subdomains is a necessary feature.
The location and number of created or vanished components then depends on
the shape of the initial locator function, as observed later in particular exam-
ples given in Section 4.3.5.4.
### 4.3.5.3 Correction procedure

As we saw that the shape of the numerically corrected interface is smoothly modified with respect to the parameter $c$ and is consistent with the above analysis, let us now investigate the properties of the discrete volume of fluid estimator function $\bar{H}_i$, representing in the bottom row of Figure 4.30 the evaluated volume of fluid value $\bar{H}_i$ against the correction parameter $c$.

<figure><img src="image_89.png" /><figcaption>(a) Control volume, initial (solid) and corrected (dashed) interfaces</figcaption></figure>

<figure><img src="image_90.png" /><figcaption>(b) Value of $\bar{H}_i$ against the correction parameter $c$</figcaption></figure>

**Figure 4.30.** Behaviour of the volume of fluid estimator function $\bar{H}_i$ and of the corrected interface location for an initially linear locator function. Top: Associated control volume and example of interfaces that are generated from values of $c$ following a logarithmic scale. Bottom: Value of the volume of fluid estimator against the parameter $c$.

First and foremost, one can observe that regardless the shape of the control volume, the continuous and monotonous behaviour of $\bar{H}_i$ obtained in Lemma 4.3.3 is preserved at the discrete level, and that any volume of fluid constraint in (0, 1) can indeed be reached. However, even though the convexity does not impact the qualitative behaviour of the estimator function, the order of magnitude of the parameter $c$ required for spanning a same range of constraint values $\bar{H}_i$ is considerably larger in the non-convex case due to the influence of the bijective mapping.

Secondly, one can assert that the observed behaviour of $\bar{H}_i$ is independent on the shape of the initial locator function, confirming here too the analysis performed in Section 4.3.2.3. Indeed, originally presented in Figure 4.30 for the case of an affine locator function, exact same observations can be made by analysing the results obtained for any initial locator function, as e.g. when considering the quintic case reported in Figure 4.31.
There, one can further confirm that the increasing or decreasing property of $\overline{H}_i$ relates only to the Cartesian parametrisation of each polygon face, and not to the shape of the initial locator function by simply comparing the results of the simplicial and convex shapes given for the same initial interface locator.

<figure><img src="image_91.png" /><figcaption>**Figure 4.31.** Behaviour of the volume of fluid estimator function $H_i$ and of the corrected interface location for an initially quintic locator function. Top: Associated control volume and example of interfaces generated from values of $c$ following a logarithmic scale. Bottom: Value of the volume of fluid estimator against the parameter $c$.</figcaption></figure>

Given those test cases, one can infer that the well-posedness property derived in Proposition 4.3.1 transfers to the discrete level for any shape of control volume and initial interface locator. However, paying attention to the shape of the volume of fluid estimator itself, one can foresee a significant numerical sensitivity of the corrected interface location towards possibly computational errors of the volume of fluid constraint. Indeed, the estimator function being very sharp around $c=0$, most of the correction cases are achieved in a narrow parameter range. Thus, as the spatial location of the interface evolves smoothly with respect to the parameter $c$, a small numerical defect can greatly impact the interface location. To emphasise this behaviour, corrected interfaces have been represented as dashed lines in the top row of Figure 4.30 and Figure 4.31 for parameters $c$ sampled in $\mathbb{R}$ in a logarithmic scale. One can notice that the regular stretching out of the interface towards the faces of the polygon is non linear, the vicinity of the original interface location being covered by interfaces that have been corrected by a small magnitude parameter. Furthermore, this non-linearity is locally amplified in the non convex case through the deformations induced by the continuous mapping.
On the bright side, the steepness of the estimator function is independent of the artificial interface thickness parameter $\epsilon$ and of the initial locator function, only depending on the number of faces of the polygon through the definition (4.39) of the corrector function. One can then *a-priori* apply a damping parameter on the corrector function $\psi_{i,c}^{\Delta}$, reducing the impact of numerical defects when solving the problem (4.49) for intermediate constraint values.

<figure><img src="image_92.png" /><figcaption>**Figure 4.32.** Corrected interfaces in several juxtaposed control volumes within which high or low volume of fluid constraints are prescribed. Top: (a)-(c) Bubbly shape interface spanning four elements. (d) Single component interface going through topological change across the control volumes. Bottom: (e)-(g) Single component interface spanning four elements. (h) No initial interface, forced into by the conservation constraint.</figcaption></figure>

Now that we assessed the legitimacy of the correction function behaviour, let us turn to practical applications and exemplify the cross cell continuity of the corrected interface. We represented in Figure 4.32 examples of interface corrections for various volume of fluid values and control volume shapes where one can see that in any case the cross cell continuity of the interface is preserved and the conservative volume of fluid constraint realised. Especially, in the top right case involving a two-components initial interface, all the interface crossings are preserved and the created subdomains are meaningful both within each control volume and in the full domain. However as already mentioned in the theoretical section, in case of a large correction amount the interface aspect is mainly determined by the correction function shape and is therefore not physically relevant.

Lastly, let us mention here again that those observations hold regardless the shape and convexity of the control volume, making the correction procedure applicable to any primal, dual, and hybrid meshes.
4.3.5.4 Impact of the interface locator function’s shape

Overlooked until now as being negligible with respect to the overall behaviour
of the interface corrector during the correction process, let us now focus on the
locator function and its impact on the spatial location of the corrected inter-
face. In particular, let us point out that although the interface is smoothly
deforming from the initial interface towards the control volume boundaries
along with the correction parameter, local extrema of the locator function
could cause spatial artefacts in the process, as represented in Figure 4.33.

<figure><img src="image_93.png" /><figcaption>**Figure 4.33.** Corrected interfaces for various parameters *c* sampled in a logarithmic scale, and given volume of fluid constraints.</figcaption></figure>

Even though those artefacts do not impact the behaviour of the estimator
function, unphysical bubbles and unwanted heavy distortions may be created
at locations that only depend on the interplay between the locator and correc-
tor functions, and are not not explicitly manageable without imposing further
constraints.

<figure><img src="image_94.png" /><figcaption>**Figure 4.34.** Interface correction of a locator function having huge variations amplitudes. (a): Initial locator function. (b): Shape of the corrected interfaces for various parameters *c* sampled on $\mathbb{R}$ in logarithmic scale. (c): volume of fluid estimator function.</figcaption></figure>

To prevent those effects while keeping a possibly high order representation of
the interface, one can always redistance the initial locator function to obtain
smoother volume of fluid estimator function and corrected interface shapes,
as it can be seen by comparing Figure 4.35 with Figure 4.36.
<figure><img src="image_95.png" /><figcaption>**Figure 4.35.** Interface correction based on an initial redistanced locator function. (a): Initial locator function. (b): Shape of the corrected interfaces for various parameters *c* sampled on $\mathbb{R}$ in logarithmic scale. (c): volume of fluid estimator.</figcaption></figure>

The initial locator function so redistanced, no spurious connected component is created during the correction process and the regularity of the interface shapes is improved. The stability of the correction procedure and the quality of the interface shape are consequently enhanced.

<figure><img src="image_96.png" /></figure>

<figure><img src="image_97.png" /><figcaption>**Figure 4.36.** Corrected interface shapes obtained for various conservation constraints. Left: $\bar{H}_i = 0.2$. Middle: $\bar{H}_i = 0.5$. Right: $\bar{H}_i = 0.8$.</figcaption></figure>

Especially, and as depicted in Figure 4.36, the corrected interface based on the redistanced interface locator tend to keep a better shape memory of the initial interface, improving the reliability of the corrected interface.
Lastly, let us note that a binding issue related to the shape of the interface
locator arises in the case of interfaces that are made of several connected
components initially crossing the control volume. There, the connectivity of
the crossing points may not be preserved whenever the shape of the control
volume yields a corrector function that parametrises an hyperbolic surface.
Indeed, for large corrections, the shape of interface corrector dominates the
one of the initial interface locator, and the hyperbolic nature of the corrector
surface will force the point connectivity to swap whenever the corrected level
set lies across a saddle point (see e.g. Figure 4.37).

<figure><img src="image_98.png" /><figcaption>**Figure 4.37.** Correction of a multi-component interface yielding a change in the boundary crossing points' connectivity. Plain line: initial interface. Dashed line: corrected interface for various correction parameters.</figcaption></figure>

Unavoidable without imposing constraints on the corrector function within
the cell, this pitfall simply highlights that by using the definition (4.55) as
part of the correction procedure, one corrects the interface location unques-
tionably to the path taken within the cell, thus not having any shape control.
As a consequence, applications requiring a high amount of correction and fre-
quent multi-component interfaces within cells would require an augmented
definition of the proposed corrector function in order to encompass adequate
constraints within the control volume.
### 4.3.6 Critical assessment of the proposed approach

As proven in Section 4.3.2.3 and observed in Section 4.3.5, the adaptation of the conservative THINC-LS correction procedure [193, 20] presented in Section 4.3.2 allows us to preserve the interface continuity of the interface by a simple modification of the corrector function's definition. As a final note, let us shortly wrap up the derived properties and observations, before going through possible perspectives of the presented work.

**Limitations** First of all, let us mention that by using the THINC-LS framework, we inherit the intrinsic drawbacks of the general correction procedure. If some are purely numerical, as achieving conservation only up to the accuracy of the employed quadrature scheme, others are structural and more problematic. Especially, as already pointed out back in Section 4.2.3, a representation of the corrected interface (4.25) in discretisation setting of the initial level set scheme does not necessarily exist, making the time evolution of the conservative level set particularly delicate. In the same vein, one may also have to apply more regularly a redistancing technique due to a possible recklessness of the locator function.

On the top of those intrinsic disadvantages, two limitations that are specific to the choice (4.39) of the corrector function are encountered. Requiring careful considerations, the first drawback consists in the trouble of preserving the connectivity of boundary crossing points of multicomponent interfaces in the case of very large corrections. Less cumbersome, the second relates to the possible need to ponder the corrector function in cells having a high number of faces, in order to lessen its steepness and improve the computational accuracy of the resolved correction parameter.

**Strengths** On the bright side, let us mention again that the proposed corrector function can be constructed on any polygonal control volume, allowing to correct any locator function to any volume of fluid constraint while preserving the continuity of the interface. In particular, when the control volume associated to the conservative constraint does not match with any mesh element of the level set discretisation, one can locally reconstruct the interpolation (4.46) of the level set field, whose shape and selected constraints enforce cross-cell continuity. The comprehensive correction procedure is thus free from any geometrical or variational constraint and can be immediately adapted to fully primal, fully dual or even staggered settings.

Applied on each control volume where the local volume of fluid constraint is contained in (0, 1), the mass conservation principle is then achieved at the discrete level, and the global continuity of the interface naturally extends the conservation principle to the phase-wise level both to any other discrete mesh and on the full computational domain $\Omega$.
In addition, the proposed definition of the corrector function enforces a
monotonic behaviour of the local volume of fluid estimator function with re-
spect to the correction parameter, property that is independent of both the
expression of the locator function and the shape of the control volume where
the correction is performed. Solving (4.49) thus becomes computationally eas-
ier and more stable than in the classical THINC-LS technique. Lastly, this
independence further ensures that any conservative constraint can always be
achieved, unlike in [20] where the feasibility of (4.49) depends on the shape
of the initial locator function. As a by product, it also allows us to safely re-
define any high order initial interface locator by redistancing it, preventing
the apparition of spurious connected components in the shape of the corrected
interface within the control volume.

**Perspectives** Let us first point out that being flexible, the proposed setting can be directly extended to curvilinear element by defining additional bijec-tive mappings, and to higher dimensions without a single modification. Fur-ther improvements can also be straightforwardly obtained, as e.g. by extend-ing the definition of the corrector function so as to include physically relevant constraints controlling the interface shape within the control volume, already initiated in Section 4.3.3.4. Among the same lines, one could also consider enhancing the cross-cell regularity of the interface whenever the locality of the determination stencil can be dropped.

On a more technical note, one could also reduce the computational sensi-
tivity on the volume of fluid constraint induced by the numerical resolution
of (4.49) by determining directly the corrected interface, considering the vol-
ume of fluid constraint as the value of a moment-type degree of freedom to
be observed by the reconstructed locator function. The correction procedure
would then be transformed into a constrained fit of a parameter-dependent
the locator function, more flexible and numerically reliable.
<figure><img src="image_99.png" /></figure>
<figure><img src="image_100.png" /></figure>
Bibliography

[1] H. Poincaré, "Sur les équations aux dérivées partielles de la physique mathématique," *American Journal of Mathematics*, pp. 211–294, 1890.

[2] D. S. Ross, "Ion etching: An application of the mathematical theory of hyperbolic conservation laws," *Journal of the Electrochemical Society*, vol. 135, no. 5, p. 1235, 1988.

[3] A. KasaHara, "Computational aspects of numerical models for weather prediction and climate simulation," in *Methods in Computational Physics: Advances in Research and Applications*, vol. 17, Elsevier, 1977, pp. 1–66.

[4] L. Alvarez-Vázquez, A. Martínez, M. Vázquez-Méndez, and M. Vilar, "An application of optimal control theory to river pollution remediation," *Applied Numerical Mathematics*, vol. 59, no. 5, pp. 845–858, 2009.

[5] M. R. Wiitala, "A dynamic programming approach to determining optimal forest wildfire initial attack responses," *Fire Economics, Planning, and Policy: Bottom Lines*, pp. 115–123, 1999.

[6] K. Lipnikov, G. Manzini, and M. Shashkov, "Mimetic finite difference method.," English, *Journal of Computational Physics*, vol. 257, 2014.

[7] L. B. Da Veiga, F. Brezzi, A. Cangiani, G. Manzini, M. L. D., and A. Russo, "Basic principles of virtual element methods," *Mathematical Models and Methods in Applied Sciences*, vol. 23, no. 01, pp. 199–214, 2013.

[8] F. Bonaldi, D. Di Pietro, G. Geymonat, and F. Krasucki, "A hybrid high-order method for Kirchhoff-Love plate bending problems.," *ESAIM: Mathematical Modelling and Numerical Analysis*, vol. 52, 2018.

[9] J. Aghili, D. Di Pietro, and B. Ruffini, "An *hp*-hybrid high-order method for variable diffusion on general meshes.," English, *Computational Methods in Applied Mathematics*, vol. 17, no. 3, pp. 359–376, 2017.

[10] L. B. Da Veiga, F. Brezzi, L. D. Marini, and A. Russo, "Mixed virtual element methods for general second order elliptic problems on polygonal meshes.," English, *ESAIM: Mathematical Modelling and Numerical Analysis*, vol. 50, no. 3, pp. 727–747, 2016.
[11] C. Talischi, A. Pereira, G. H. Paulino, I. F. M. Menezes, and M. S. Carvalho, "Polygonal finite elements for incompressible fluid flow," *International Journal for Numerical Methods in Fluids*, vol. 74, no. 2, pp. 134–151, 2013.

[12] L. Botti, D. Di Pietro, and J. Droniou, "A hybrid high-order method for the incompressible Navier-Stokes equations based on Temam's device," English, *Journal of Computational Physics*, vol. 376, pp. 786–816, 2019.

[13] R. Loubère, P. H. Maire, and M. Shashkov, "ReALE: a Reconnection Arbitrary-Lagrangian-Eulerian method in cylindrical geometry," English, *Computers and Fluids*, vol. 46, no. 1, pp. 59–69, 2011.

[14] S. K. Godunov, "A difference method for numerical calculation of discontinuous solutions of the equations of hydrodynamics," *Matematicheskii Sbornik*, vol. 89, no. 3, pp. 271–306, 1959.

[15] A. Harten, B. Engquist, S. Osher, and S. R. Chakravarthy, "Uniformly high order accurate essentially non-oscillatory schemes, III," in *Upwind and high-resolution schemes*, Springer, 1987, pp. 218–290.

[16] A. Harten, P. D. Lax, and B. v. Leer, "On upstream differencing and Godunov-type schemes for hyperbolic conservation laws," *SIAM review*, vol. 25, no. 1, pp. 35–61, 1983.

[17] G.-S. Jiang and C.-W. Shu, "Efficient implementation of weighted ENO schemes," *Journal of computational physics*, vol. 126, no. 1, pp. 202–228, 1996.

[18] T. Chen and C.-W. Shu, "Entropy stable high order discontinuous Galerkin methods with suitable quadrature rules for hyperbolic conservation laws," *Journal of Computational Physics*, vol. 345, pp. 427–461, 2017.

[19] H. T. Huynh, "A flux reconstruction approach to high-order schemes including discontinuous Galerkin methods," in *18th AIAA Computational Fluid Dynamics Conference*, 2007, p. 4079.

[20] R. Kumar, L. Cheng, B. Xie, and F. Xiao, "THINC-scaling scheme that unifies VOF and level set methods," *arXiv preprint arXiv:1908.04529*, 2019.

[21] E. Godlewski and P.-A. Raviart, *Numerical approximation of hyperbolic systems of conservation laws*. Springer Science & Business Media, 2013, vol. 118.

[22] S. N. Kružkov, "First order quasilinear equations in several independent variables," *Mathematics of the USSR-Sbornik*, vol. 10, no. 2, p. 217, 1970.

[23] A. Bressan, "An ill posed Cauchy problem for a hyperbolic system in two space dimensions," *Rendiconti del Seminario Matematico della Università di Padova*, vol. 110, pp. 103–117, 2003.
[24] S. Bianchini and A. Bressan, “Vanishing viscosity solutions of nonlinear hyperbolic systems,” *Annals of mathematics*, pp. 223–342, 2005.

[25] L. Ambrosio, F. Bouchut, and C. D. Lellis, “Well-Posedness for a Class of Hyperbolic Systems of Conservation Laws in Several Space Dimensions,” *Communications in Partial Differential Equations*, vol. 29, no. 9-10, pp. 1635–1651, 2005.

[26] C. M. Dafermos, “Generalized characteristics and the structure of solutions of hyperbolic conservation laws,” *Indiana University Mathematics Journal*, vol. 26, no. 6, pp. 1097–1119, 1977.

[27] E. Conway and J. Smoller, “Global solutions of the Cauchy problem for quasi-linear first-order equations in several space variables,” *Communications on Pure and Applied Mathematics*, vol. 19, no. 1, pp. 95–105, 1966.

[28] C. De Lellis, *Notes on hyperbolic systems of conservation laws and transport equations*. University of Zurich, Institute of Mathematics, 2006.

[29] J. Glimm, “Solutions in the large for nonlinear hyperbolic systems of equations,” *Communications on pure and applied mathematics*, vol. 18, no. 4, pp. 697–715, 1965.

[30] P. D. Lax, “Hyperbolic systems of conservation laws II,” *Communications on pure and applied mathematics*, vol. 10, no. 4, pp. 537–566, 1957.

[31] T. Nishida, “Global solution for an initial boundary value problem of a quasilinear hyperbolic system,” *Proceedings of the Japan Academy*, vol. 44, no. 7, pp. 642–646, 1968.

[32] T.-P. Liu, “The Riemann problem for general systems of conservation laws,” *Journal of Differential Equations*, vol. 18, no. 1, pp. 218–234, 1975.

[33] O. A. Oleinik, “Discontinuous solutions of non-linear differential equations,” *Uspekhi Matematicheskikh Nauk*, vol. 12, no. 3, pp. 3–73, 1957.

[34] B. L. Keyfitz and H. C. Kranzer, “Existence and uniqueness of entropy solutions to the Riemann problem for hyperbolic systems of two nonlinear conservation laws,” *Journal of Differential Equations*, vol. 27, no. 3, pp. 444–476, 1978.

[35] T.-P. Liu and K. Zumbrun, “On nonlinear stability of general undercompressive viscous shock waves,” *Communications in mathematical physics*, vol. 174, no. 2, pp. 319–345, 1995.

[36] J. Goodman, “Nonlinear asymptotic stability of viscous shock profiles for conservation laws,” *Archive for Rational Mechanics and Analysis*, vol. 95, no. 4, pp. 325–344, 1986.

[37] C. M. Dafermos, “The entropy rate admissibility criterion for solutions of hyperbolic conservation laws,” *Journal of Differential Equations*, vol. 14, no. 2, pp. 202–212, 1973.
[38] —, "Quasilinear hyperbolic systems that result from conservation laws(with application to wave propagation)," *Nonlinear waves*. (A 74-20717 07-23) *Ithaca, N. Y.*, Cornell University Press, 1974, pp. 82–102, 1974.

[39] P. D. Lax, "Hyperbolic systems of conservation laws in several space variables," New York Univ., NY (USA). Courant Mathematics and Computing Lab., Tech. Rep., 1985.

[40] Y. Wang and C. Wei, "Completely linear degeneracy for quasilinear hyperbolic systems," *Journal of Mathematical Analysis and Applications*, vol. 465, no. 1, pp. 318–330, 2018.

[41] C. De Lellis and L. Székelyhidi, "On admissibility criteria for weak solutions of the Euler equations," *Archive for rational mechanics and analysis*, vol. 195, no. 1, pp. 225–260, 2010.

[42] M. V. Kondratieva, A. B. Levin, A. V. Mikhalev, and E. V. Pankratiev, *Differential and Difference Dimension Polynomials*. Springer Netherlands, 1999, vol. 461.

[43] J.-M. Thomas, "Méthode des éléments finis hybrides duaux pour les problèmes elliptiques du second ordre," *Revue française d'automatique, informatique, recherche opérationnelle*. Analyse numérique, vol. 10, 1976.

[44] S. T. Miller and R. B. Haber, "A spacetime discontinuous Galerkin method for hyperbolic heat conduction," *Computer Methods in Applied Mechanics and Engineering*, vol. 198, no. 2, pp. 194–209, 2008.

[45] W. Dörfler, S. Findeisen, and C. Wieners, "Space-time discontinuous Galerkin discretizations for linear first-order hyperbolic evolution systems," *Computational Methods in Applied Mathematics*, vol. 16, no. 3, pp. 409–428, 2016.

[46] R. L. Hardy, "Multiquadric equations of topography and other irregular surfaces," *Journal of geophysical research*, vol. 76, no. 8, pp. 1905–1915, 1971.

[47] E. J. Kansa, "Multiquadrics - A scattered data approximation scheme with applications to computational fluid-dynamics - II solutions to parabolic, hyperbolic and elliptic partial differential equations," *Computers & mathematics with applications*, vol. 19, no. 8-9, pp. 147–161, 1990.

[48] D. Levy, G. Puppo, and G. Russo, "Central WENO schemes for hyperbolic systems of conservation laws," *ESAIM: Mathematical Modelling and Numerical Analysis*, vol. 33, no. 3, pp. 547–571, 1999.

[49] R. M. Beam and R. F. Warming, "An implicit finite-difference algorithm for hyperbolic systems in conservation-law form," *Journal of computational physics*, vol. 22, no. 1, pp. 87–110, 1976.

[50] R. Abgrall, "High order schemes for hyperbolic problems using globally continuous approximation and avoiding mass matrices," *Journal of Scientific Computing*, vol. 73, no. 2, pp. 461–494, 2017.
[51] P. Lax, “Systems of conservation laws,” Los Alamos National Lab., Tech. Rep., 1959.

[52] R. Abgrall, E. Le Mélédo, and P. Öffner, “On the Connection between Residual Distribution Schemes and Flux Reconstruction,” 2018.

[53] W. H. Reed and T. R. Hill, “Triangular mesh methods for the neutron transport equation,” Los Alamos Scientific Lab., Tech. Rep., 1973.

[54] B. Cockburn and C.-W. Shu, “The Runge-Kutta local projection-discontinuous-Galerkin finite element method for scalar conservation laws,” *ESAIM: Mathematical Modelling and Numerical Analysis*, vol. 25, 1991.

[55] A. Jameson, “A proof of the stability of the spectral difference method for all orders of accuracy,” *Journal of Scientific Computing*, vol. 45, no. 1-3, pp. 348–358, 2010.

[56] P. Castonguay, P. E. Vincent, and A. Jameson, “A new class of high-order energy stable flux reconstruction schemes for triangular elements,” *Journal of Scientific Computing*, vol. 51, no. 1, pp. 224–256, 2012.

[57] G. J. Gassner, “A skew-symmetric discontinuous Galerkin spectral element discretization and its relation to SBP-SAT finite difference methods,” *SIAM Journal on Scientific Computing*, vol. 35, no. 3, A1233–A1253, 2013.

[58] G. J. Gassner and D. A. Kopriva, “A comparison of the dispersion and dissipation errors of Gauss and Gauss-Lobatto discontinuous Galerkin spectral element methods,” *SIAM Journal on Scientific Computing*, vol. 33, no. 5, pp. 2560–2579, 2011.

[59] F. D. Witherden and P. E. Vincent, “An analysis of solution point coordinates for flux reconstruction schemes on triangular elements,” *Journal of Scientific Computing*, vol. 61, no. 2, pp. 398–423, 2014.

[60] Y. Abe, I. Morinaka, T. Haga, T. Nonomura, H. Shibata, and K. Miyaji, “Stable, non-dissipative, and conservative flux-reconstruction schemes in split forms,” *Journal of Computational Physics*, vol. 353, pp. 193–227, 2018.

[61] H. Ranocha, P. Öffner, and T. Sonar, “Summation-by-parts operators for correction procedure via reconstruction,” *Journal of Computational Physics*, vol. 311, pp. 299–328, 2016.

[62] V. Ervin, “Computational bases for RTk and BDMk on triangles,” *Computers & Mathematics with Applications*, vol. 64, no. 8, pp. 2765–2774, 2012.

[63] R. Abgrall, É. Le Mélédo, and P. Öffner, “A class of finite dimensional spaces and $H$-(div) conformal elements on general polytopes,” 2019.

[64] S. Clain, S. Diot, and R. Loubère, “A high-order finite volume method for systems of conservation laws—Multi-dimensional Optimal Order Detection (MOOD),” *Journal of computational Physics*, vol. 230, no. 10, pp. 4028–4050, 2011.