# Parabolic interface reconstruction for 2D volume of fluid methods

For capillary driven flow the interface curvature is essential in the modelling of surface tension via the imposition of the Young-Laplace jump condition. We show that traditional geometric volume of fluid (VoF) methods, that are based on a piecewise linear approximation of the interface, do not lead to an interface curvature which is convergent under mesh refinement in time-dependent problems. Instead, we propose to use a piecewise parabolic approximation of the interface, resulting in a class of piecewise parabolic interface calculation (PPIC) methods. In particular, we introduce the parabolic LVIRA and MoF methods, PLVIRA and PMoF, respectively. We show that a Lagrangian remapping method is sufficiently accurate for the advection of such a parabolic interface. It is numerically demonstrated that the newly proposed PPIC methods result in an increase of reconstruction accuracy by one order, convergence of the interface curvature in time-dependent advection problems and Weber number independent convergence of a droplet translation problem, where the advection method is coupled to a two-phase Navier–Stokes solver.

There are no comments yet.

## Authors

• 1 publication
• 1 publication
12/23/2021

### Analysis of the diffuse interface method for the Stokes-Darcy coupled problem

We consider the interaction between a free flowing fluid and a porous me...
04/23/2021

### Consistent and symmetry preserving data-driven interface reconstruction for the level-set method

Recently, machine learning has been used to substitute parts of conventi...
06/23/2020

### Analytic Solution to the Piecewise Linear Interface Construction Problem and its Application in Curvature Calculation for Volume-of-Fluid Simulation Codes

The plane-cube intersection problem has been around in literature since ...
06/26/2020

### NPLIC: A Machine Learning Approach to Piecewise Linear Interface Construction

Volume of fluid (VOF) methods are extensively used to track fluid interf...
06/11/2021

### Projection-based resolved interface mixed-dimension method for embedded tubular network systems

We present a flexible discretization technique for computational models ...
08/05/2021

### a Decision-Tree based Moment-of-Fluid (DTMOF) Method in 3D rectangular hexahedrons

The moment-of-fluid (MOF) method is an extension of the volume-of-fluid ...
09/15/2017

### Generalized Internal Boundaries (GIB)

Representing large-scale motions and topological changes in the finite v...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

The advection of the phase interface plays a central role in the simulation of immiscible, and in our case incompressible, two-phase flow. Therefore, much research has been performed towards the development of accurate, efficient and robust interface advection methods. Explicit modelling of the interface (Tryggvason2001) results in a highly accurate and continuous representation, for which the interface resolution is essentially unrelated to the resolution of the mesh. Changes in topology can be taken into account (Bo2011), but this is nontrivial.

The interface can also be represented implicitly, either using a level set or a volume fraction function, resulting in level set and volume of fluid (VoF) methods respectively. Using level sets (Osher1988; Gibou2017) results in an efficient method for which changes in topology are automatically taken into account. Level set methods however do not inherently conserve mass and require artificial redistancing to ensure that the level set function remains a distance function (Sussman1994). On the contrary, VoF methods (see for instance the works of Hirt1981; Youngs1982; Rider1997; Puckett1997; Rudman1998; Harvie2001; Lopez2004; Weymouth2010) can inherently conserve mass and result, in particular if geometric VoF methods are considered, in a sharper interface. Geometric VoF methods rely on the geometric reconstruction of the interface which results in a sharp interface representation (see e.g. fig. 4) thus reflecting the assumed immiscibility of the fluids. On the other hand, algebraic VoF methods lack a geometric reconstruction of the interface, and result in a diffuse interface. In this paper we will focus our attention on geometric VoF methods that are based on Lagrangian remapping; this is discussed in section 2.

Level set methods are sometimes preferred for their ‘smoother’ interface representation, thus resulting in the ability to easily approximate the local interface geometry, such as the interface normal vector as well as the interface curvature. In our experience however, a sufficiently accurate volume fraction field equally provides the ability to compute the interface normal vector as well as the interface curvature at high accuracy, the latter using local height-functions (LHFs). What we exactly mean by ‘sufficiently accurate’ is discussed in

section 3.

For capillary driven flow the interface curvature is essential in the modelling of surface tension via the imposition of the Young-Laplace jump condition. When traditional geometric reconstruction methods, that are based on a piecewise linear approximation of the interface (see fig. 4 and section 4), are used, we find that the curvature from LHFs does not converge under mesh refinement for time-dependent problems. Nonconvergence of the interface curvature results in spurious and unphysical currents, as discussed in (Magnini2016). A theoretical understanding is provided in sections 3 and 2 where we analyse the accuracy of the advection method and of the LHF-based curvature respectively. These results are confirmed numerically in sections 6.3 and 6.2. Based on this observed lack of convergence, we propose to use a piecewise parabolic approximation of the interface instead, as shown in fig. 4 and discussed in section 5. This essentially includes the curvature in the interface advection, and in section 6 we will demonstrate that this leads to convergence of the curvature, also for time-dependent problems.

Throughout this paper we will denote the positive time-step by , and the maximum diameter of all control volumes by . Furthermore, a numerical approximation to some quantity is denoted by where we never explicitly indicate the dependence on . The statement ‘ converges’ is understood as converges to under mesh refinement

 limh→0˜y(δ,h)=y, (1)

where the time-step satisfies the CFL time-step restriction

 δ≤hU, (2)

such that implies .

## 2 Advection methods based on Lagrangian remapping

In this section we discuss the advection method that will be used to approximately track the interface , and show that this advection method is indeed sufficiently accurate for the advection of a parabolic interface approximation, as we propose in section 5.

We denote the computational domain by , the interface between the two fluids by and the corresponding liquid domain by .

The zeroth and first moment of some set

will be denoted by

respectively (the first moment is introduced for use in the MoF reconstruction method discussed in section 4.2). Given the zeroth and first moment, the centroid is given by

 C(A):=M1(A)M0(A). (4)

The liquid volume and its control volume fraction are then defined as

 Ml0,c:=M0(c∩Ωl),¯χc(t):=Ml0,cM0(c), (5)

for some control volume . The liquid first moment and centroid are similarly denoted by and respectively.

### 2.1 Lagrangian remapping

We let denote the flow map of the velocity field over the time interval , that is, solves the initial value problem

 ddtx=u(t,x(t)),x(t(n))=x0, (6)

resulting in . The action of the flow map naturally extends to the application on sets of positions, such as . The preimage of the control volume under the flow map, which we denote by , is the set of points which end up inside under the flow map , and will henceforth be referred to as ‘the preimage of ’. See also fig. 7.

Note that since the interface is advected with the velocity field , it follows that

 I(n+1)=ΨδI(n)⇒Ωl,(n+1)=ΨδΩl,(n). (7)

By the invertibility of it holds that and therefore

 c∩Ωl,(n+1)=Ψδ(Ψ−δc∩Ωl,(n)). (8)

Computing the zeroth moment of eq. 8, and by making use of the fact that the velocity field is divergence free (which implies that the flow map is area preserving, i.e. ) results in

 Ml,(n+1)0,c=M0(Ψ−δc∩Ωl,(n)). (9)

Equation (9) is key to understanding how a Lagrangian remapping-based geometric VoF method can be constructed, since such a method approximates each of the terms on the right-hand side of eq. 9:

• The liquid domain is approximated per control volume. This will be discussed in sections 5 and 4 for the piecewise linear and piecewise parabolic approximation of the interface respectively.

• The preimage is approximated by a polygon whose vertices are approximated using numerical integration, leading to two approximation errors. This is what we will discuss next.

### 2.2 Approximate Lagrangian remapping

We denote the exact preimage of the control volume by

 Pc:=Ψ−δc. (10)

The approximate preimage is constructed by making two approximations. First we approximate the preimage by a polygonal representation denoted by , which is defined by connecting the corners of the exact preimage by straight line segments, as is shown in fig. 11.

Secondly, we approximate the flow map using a numerical integration method. This means we approximately integrate along pathlines, where the velocity field

is now linearly interpolated (in space and time) from a staggered velocity field. For the time integration we use the second-order accurate Heun method, which, as we will see, is sufficiently accurate for this purpose. This results in the following approximate liquid volume

 ˜Ml0,c=M0(˜Pc∩˜Ωlc). (11)

The method described here is very similar to the Lagrangian-Eulerian advection scheme (LEAS) presented in Zinjala2015.

### 2.3 Error analysis

We will now estimate the error in the volume fractions, due to the two approximations made in the preimage as well as the approximation of the liquid domain. To this end we denote the symmetric difference of two sets

by

 A△B:=(A∪B)∖(A∩B). (12)

Using the symmetric difference we define the reconstruction error as

 ERec0:=M0(˜Ωlc△Ωlc). (13)

Here denotes the liquid neighbourhood centred around the control volume

 Ωlc=Ωl∩⋃c′∈C(c)c′, (14)

where is the set of control volumes which share at least one vertex with (i.e. a neighbourhood of control volumes if the mesh is rectilinear). Moreover we define the polygonal representation error

 ERep0:=M0(PRepc△Pc), (15)

and the integration error

 EInt0:=M0(˜Pc△PRepc). (16)

The three errors are illustrated in fig. 11.

The following lemma, which is based on Zhang2013, shows how the error resulting from a Lagrangian remapping method can be bounded by the previously introduced errors.

###### Lemma 1 (Error decomposition).

The approximation error of a Lagrangian remapping method can be bounded by

 ∣∣˜Ml0,c−Ml0,c∣∣≤E% Rec0+ERep0+EInt0. (17)

Furthermore, lemma 2 provides estimates for the approximation errors of the preimage. Proofs are found in appendix A.

###### Lemma 2 (Remapping error estimates).

The representation and integration errors are given by

 ERep0=O(hδ(h+δ)2),EInt0=O(hδ(h+δ)2+hδq+1), (18)

where is the order of accuracy of the time integration method.

Combining both lemmas allows us to prove the following consistency result for the Lagrangian remapping method described in section 2.2.

###### Theorem 1 (Consistency of a Lagrangian remapping method).

The following single time step consistency result holds if (as is the case under a CFL time step restriction)

 (19)

where denotes the local order of accuracy of the liquid domain approximation (e.g. for the MoF method (Dyadechko2005))

 ERec0=O(hr), (20)

and denotes the order of accuracy of the time integration method (e.g. for Heun’s method).

###### Proof.

The -norm of some scalar function , which is defined on the grid, is defined as its maximal value in absolute sense

 ∥α∥L∞:=maxc∈C|αc|, (21)

and thus

 (22)

which follows from lemma 1. Then, using the result of lemma 2 and the assumed local order of accuracy of the liquid domain approximation, we find that (assuming )

 ∥∥˜¯χ−¯χ∥∥L∞≤O(hr)+O(h4)+O(h4+hq+2)h2=O(hr−2)+O(h2)+O(h2+hq). (23)

The result of theorem 1 is limited to a single time-step. When the Lagrangian remapping method is not coupled to a Navier–Stokes solver, then the only possible source of error propagation is via the interface reconstruction. Hence for generalising theorem 1 to a multiple time-step convergence result, the stability of the interface reconstruction method w.r.t. perturbations in the reference moments (i.e. errors from the previous time-step) must be analysed. Thus far we are not aware of any such results.

From theorem 1 we find that the accuracy of the volume fractions is limited by the reconstruction accuracy whenever a piecewise linear approximation of the interface is used, for which (we numerically demonstrate this in section 6.1 and it is claimed in (Dyadechko2005) for the MoF method). An increase of the reconstruction accuracy to , would yield an improvement of the -error of the volume fractions from first- to second-order accuracy, while using the same definition of the approximate preimage. This is achieved in section 5 where we consider the piecewise parabolic reconstruction of the interface.

## 3 Curvature convergence

We now turn to the observed lack of curvature convergence in time-dependent problems, as discussed in section 1. For simplicity in presentation we assume (only in this section) a uniform rectilinear mesh, and denote each control volume by an index . The control volume centroid is now denoted by .

### 3.1 Curvature computation using a local height-function

The curvature is computed for each interface control volume using a LHF resulting from a principal normal direction which is aligned with one of the co-ordinate axes111For unstructured meshes LHFs can also be used, see Owkes2015., consider for example fig. 14 where the principal normal direction coincides with the positive -direction. Given such a LHF, which we denote by , the curvature can be computed as follows

 κ(y)=H[2](y)[1+(H[1](y))2]−3/2, (24)

where denote the first and second derivative of , respectively.

Rather than using the LHF directly, the curvature is instead computed using the averaged LHF, which is defined as

 (25)

which by the midpoint rule is a second-order accurate approximation of  (quarteroni2010). The reason for using the averaged LHF rather than the LHF itself, is that for rectilinear meshes the former can easily be computed by summing over adjacent volume fractions

 ¯Hj:=¯H(yj)=1hNH∑k=−NHh2¯χi+k,j−(NH+12)h, (26)

where we have omitted the dependence on the index in order to simplify the notation. See also fig. 14. The parameter should be sufficiently large to be able to determine the validity of the averaged LHF: the bottom and top control volumes must be entirely full and empty respectively () and the volume fractions must be monotonically decreasing ( for ). If the principal normal direction coincides with the main direction of the interface normal, then for a sufficiently well-resolved interface these conditions are satisfied for , as is the case in fig. 14 and used by Afkhami2007.

The averaged LHF can be approximated by summing over adjacent and approximate volume fractions

 ˜¯Hj:=1hNH∑k=−NHh2˜¯χi+k,j−(NH+12)h. (27)

Provided that the volume fractions are -th order accurate, this results in a -st order accurate approximation to . Recall that according to theorem 1, where is the order of accuracy of the liquid domain approximation and is the order of accuracy of the time integration method. Hence for traditional piecewise linear interface reconstruction methods, for which , we find that the averaged LHF is approximated at second-order accuracy.

The fully approximated curvature is then given by

 ˜κj:=˜¯H˜[2](y)[1+(˜¯H˜[1](y))2]−3/2, (28)

where the second-order finite difference approximations to the first and second derivatives of some function are defined as

 f˜[1]j:=fj+1−fj−12h,f˜[2]j:=fj+1−2fj+fj−1h2. (29)

Hence three consecutive values of the LHF are needed for the LHF based approximation of the interface curvature. Whenever this is not possible we combine LHFs from different principal normal directions, as proposed in (Popinet2009), resulting in the generalised height-function (GHF) method.

### 3.2 Error analysis

The following lemma provides estimates for the approximation errors resulting from first and second derivatives of the LHF when second-order finite differences are used with the approximated and averaged LHFs.

###### Lemma 3 (Derivative error estimates).

Assuming that and that -th order accurate volume fractions are used to approximate the averaged LHFs, we find

 ˜¯H˜[λ]j=H[λ]j+O(h2)+O(hp+1−λ),λ=1,2. (30)
###### Proof.

Let the approximation error of the averaged LHF be defined as

 ¯e(y;h):=¯H(y)−H(y), (31)

where by the assumed smoothness of . By the midpoint rule we find that . It follows that

 ¯H˜[λ]j=H˜[λ]j+¯e˜[λ](yj;h)=[H[λ]j+O(h2)]+[¯e[λ](x;h)+O(h4)]=H[λ]j+O(h2),λ=1,2, (32)

by the differentiability of .

Similarly, we let the approximation error due to the approximation of the volume fractions be defined as

 ˜¯e(y;h):=˜¯H(y)−¯H(y), (33)

where we have defined

 ˜¯H(y):=1h∫y+h/2y−h/2˜H(τ)dτ. (34)

Here is the approximate LHF based on the piecewise approximation of the interface, see also fig. 14. By the assumed -th order accuracy of the volume fractions, we find that . Contrary to , we can not show differentiability of w.r.t. due to the discontinuity of at and (see also fig. 14 where the error is clearly non differentiable at ). This implies that approximation errors in the volume fractions will be amplified due to the division by and in the approximation of the first and second derivative respectively

 ˜¯H˜[λ]j=¯H˜[λ]j+˜¯e˜[λ](yj;h)=¯H˜[λ]j+O(hp+1−λ),λ=1,2. (35)

Combining eqs. 35 and 32 yields the desired result. ∎

Lemma 3 implies that when a piecewise linear approximation of the interface is used, for which , we find that the volume fractions will be first order accurate () and thus the approximated second derivative is inconsistent, resulting in a lack of curvature convergence. Convergence of the second derivative of the LHF requires a more accurate liquid domain approximation for which (which results in ), in turn this would then also lead to a convergent curvature, as we numerically demonstrate in sections 6.3 and 6.2.

Obtaining a more accurate liquid domain approximation for which will be discussed in section 5 where we introduce piecewise parabolic approximations of the interface as generalisations of the piecewise linear approximations of the interface that are discussed in section 4.

## 4 Optimisation-based PLIC methods

We will present the piecewise linear interface calculation (PLIC) methods, that are based on the optimisation of some cost function over some search space (e.g. the space of linear interfaces ), in a rather abstract way that straightforwardly allows us to generalise any such PLIC method to a parabolic interface reconstruction method. This generalisation does not alter the cost function , but rather replaces the search space with a larger search space which includes parabolic interfaces, which is the topic of section 5.

We consider two PLIC methods in the context of optimisation-based reconstruction methods. That is, the reconstruction step for each control volume can be written as an optimisation problem

 p∗=argminp∈P1f(p), (36)

where is the search space consisting of linear interfaces and is a cost function. More precisely, the search space is defined as

 P1:={p(x)=ηT(x−xc)−s(η;Ml0,c)∣η∈S1}, (37)

where is the unit 1-sphere and denotes the centroid of the control volume . The shift is uniquely defined by requiring that the reconstructed liquid volume equals , and therefore the only remaining unknown is the interface normal .

The resulting interface inside the control volume is defined as the zero level set of . Similarly, the subset of for which (i.e. defining a half-space if is linear) is denoted by

 l(p):={x∈R2∣p(x)≤0}, (38)

and is defined such that defines the subset of that contains the liquid phase. Given the optimal level set for each control volume , we define the approximate liquid domain as

 ˜Ωl:=⋃c∈Cc∩l(p∗c), (39)

where is the set of all control volumes. See also fig. 4.

Having defined the search space, all that remains is the definition of appropriate cost functions. Note that such a cost function aims at approximating the actual error in terms of the area of the symmetric difference

 f△(p):=M0((c∩l(p))△(c∩Ωl))M0(c). (40)

### 4.1 The least squares VoF interface reconstruction algorithm

The least squares VoF interface reconstruction algorithm (LVIRA) (puckett1991volume) uses the following cost function

 fL2(p):=⎡⎢⎣∑c′∈C(c)M0(c′)⎛⎝Ml0,c′−M0(c′∩l(p))M0(c′)⎞⎠2⎤⎥⎦1/2, (41)

where denotes the set of control volumes which share at least one vertex with . That is, the -norm of the difference between the actual and reconstructed (by extending the interface to ) volume fractions is considered in a neighbourhood around the control volume . See also fig. 18.

This means that if the volume fractions originate from a globally linear interface, then the reconstruction will be exact.

The cost function does not always have a unique global minimum222For example, if for a uniform mesh the volume fraction field for has the symmetry , then multiple optimal solutions exist to eq. 36., but for a resolved interface which can be expressed as a LHF (with any principal normal direction), we find that such problems do not arise.

An efficient LVIRA (ELVIRA) method was proposed (Pilliod2004) for rectilinear meshes, where six trial normals are considered, and the one for which is smallest is then selected. This means that no numerical optimisation is required. It should be noted however that the resulting normal may not be optimal. Nevertheless, it can still be shown (Pilliod2004) that the reconstruction will be exact if the volume fractions originate from a globally linear interface.

### 4.2 The moment of fluid reconstruction algorithm

The moment of fluid (MoF) (Dyadechko2005) method proposes to include a reference first moment for each control volume to determine the interface normal. Given such a reference first moment, the interface normal is determined by requiring the reconstructed first moment to be as close as possible to this reference first moment, this is achieved via optimisation of the following cost function

 fM1(p):=|Ml,∗1,c−M1(c∩l(p))|2, (42)

where denotes the Euclidean norm. Contrary to other PLIC methods, the MoF method uses only information (the zeroth and first moment) from the control volume itself, resulting in increased accuracy in particular when the interface is only nearly resolved (). Of course, if the reference moments and originate from a linear interface, then there must exist a normal angle and thus for which and hence the MoF method reconstructs linear interfaces exactly.

For the MoF method it is known that not all reference moments result in stability w.r.t. perturbations of the reference first moment (Dyadechko2005), just as there are reference moments which result in multiple solutions333For example, if the reference centroid coincides with the centroid of a square control volume while , then for any normal direction the same value of the cost function is found for any other normal direction provided that and . to the optimisation eq. 36. The set of reference moments for which multiple solutions exist has measure zero, and therefore this does not pose a problem in practice (Dyadechko2005).

In two spatial dimensions an efficient variant of the MoF method was proposed (Lemoine2017) for rectilinear meshes, which uses three trial normals. The construction of trial normals is such that the trial normal for which is minimal is indeed the optimal one and therefore this efficient variant is merely a new solution strategy for the MoF method.

### 4.3 Advection of the first moment

The MoF reconstruction method requires a reference first moment, and therefore the first moment must also be advected. The advection equation of the first moment can be derived by computing the first moment of eq. 8

 Ml,(n+1)1,c=M1(Ψδ(Ψ−δc∩Ωl,(n))). (43)

Hence after the preimage of the control volume is intersected with the liquid domain at , this intersection must be advected forward in time before computing the first moment. We choose to approximate this last step by instead advecting only the centroid of this intersection forward in time (as a point particle), resulting in the following approximation (Dyadechko2005)

 Cl,(n+1)c=Ψδ(C(Ψ−δc∩Ωl,(n)))+O(δh2). (44)

#### 4.3.1 Approximate advection of the first moment

Numerically the advection of the first moment is done in three steps, as illustrated in figs. 22, 22 and 22

###### Steps 1 (Centroid advection).

The liquid first moment of the intersection of the approximate preimage with the approximate liquid domain is computed

 ˜Ml,∗∗1,c=M1(˜Pc∩˜Ωl,(n)c). (45)

The centroid of this intersection is advected as a point particle using the linearly interpolated velocity field (in space and time)

 ˜Cl,∗c=˜Cl,∗∗c+δu(n)c⇒˜Ml,∗1,c=˜Ml,∗∗1,c+˜Ml,(n+1)0,cδu(n)c, (46)

where for simplicity in presentation we consider here the forward Euler method for approximation of the time-integration in eq. 44.

The interface is reconstructed by minimising the MoF cost function (eq. 42) using as the reference first moment. The first moment corresponding to the reconstructed interface, which is defined as the zero level set of some , is denoted by

 ˜Ml,(n+1)1,c=M1(c∩l(p∗)). (47)

Note that for the advection of the zeroth moment, as given by eq. 11, we considered only step 1, since steps 1 and 1 do not alter the zeroth moment.

#### 4.3.2 Error analysis

Each of the steps introduces several errors. For the error analysis of step 1 we note that for any

 (48)

Hence if we define the first moments and relative to the centroid of the control volume , then we can directly use the results of section 2 to find the estimates

 ∣∣ERec1∣∣2 :=∣∣M1(˜Ωlc△Ωlc)∣∣2=O(hr+1) (49) ∣∣ERep1∣∣2 (50) ∣∣EInt1∣∣2 (51)

It follows that if

 ∣∣∣˜Ml,∗∗1,c−Ml,∗∗1,c∣∣∣2≤∣∣M1((˜Pc∩˜Ωlc)△(Pc∩Ωlc))∣∣2=O(hr+1)+O(h5)+O(hq+3). (52)

Note that the value of the first moment itself is .

Step 1 introduces the approximation error due to advecting the centroid as a point particle (eq. 44) as well as a time integration error, where we again assume that a method of order is used. In total this results in the following error estimate for the first two steps

 ∣∣∣˜Ml,∗1,c−Ml,∗1,c∣∣∣2=O(hr+1)reconstruction+O(h5+hq+3)preimage approximation+O(h5)\lx@cref{creftype~refnum}{eqn:plic:mome% nt_advection:point_particle}+O(hq+3)time % integration. (53)

Finally, in step 1 the interface is reconstructed according to the MoF cost function (eq. 42), resulting in some optimal for which the error made in the third step is minimised

 ∣∣∣˜Ml,∗1,c−˜Ml,(n+1)1,c∣∣∣2=∣∣∣˜Ml,∗1,c−M1(c∩l(p∗))∣∣∣2=fM1(p∗). (54)

In Dyadechko2005 the following estimate is given for the approximation error made in the first moment when the MoF method is used

 fM1(p∗)=O(h5),p∗∈P1, (55)

which assumes that the interface that is being approximated is a curve with a bounded (from below) radius of curvature. After a single time-step however, the interface has become piecewise smooth, and for the general non-smooth case the authors of (Dyadechko2005) provide the pessimistic estimate

 fM1(p∗)=O(h3), (56)

which, since the value of the first moment itself is also , corresponds to a relative error of . It follows that, even for a single time-step, we can not prove consistency of the first-moment when the MoF cost function is used.

The numerical results shown in fig. 38 however suggest that

 fM1(p∗)=O(h4), (57)

when a PLIC reconstruction is used. Hence the non-smooth estimate eq. 56 is indeed found to be too pessimistic.

## 5 Parabolic reconstruction methods

The PLIC methods introduced in section 4 result in a local order of accuracy of the liquid domain approximation of (this is claimed in (Dyadechko2005) for the MoF method and observed numerically in section 6.1). Therefore, using the result of theorem 1 we find that the volume fractions are at most first order accurate, which using lemma 3 with implies that the second derivative of the LHF, and therefore the curvature, does not converge. To this end we now consider the generalisation of the reconstruction methods discussed in section 4 to parabolic reconstruction, which will increase the local order of accuracy of the liquid domain approximation to .

### 5.1 Parabolic search spaces

We define the restricted parabolic search space as

 P2:={p(x)=ηT(x−xc)−s(η,κ;Ml0,c)+κ2(τT(x−xc))2∣η∈S1,κ∈R}, (58)

where is the interface tangent and denotes the curvature. The search space is restricted in the sense that the parabola can be defined relative to any point , but we instead fix it around , this reduces dimensionality of by one444The ‘unrestricted’ search space includes a tangential shift and is given by

(59)
which corresponds to . . Furthermore we denote by the subspace of where the curvature is known a priori (we let the curvature be approximated by the GHF method (Popinet2009)), this reduces the dimensionality of the search space by one once more.

At this point we have two cost functions, and (eqs. 42 and 41), as discussed in section 4, as well as three search spaces: and . In table 1 we show the names of the methods for each of the possible combinations that can be formed.

Reconstruction methods based on or will be referred to as piecewise parabolic interface calculation (PPIC) methods.

### 5.2 PLVIRA and PROST

The parabolic reconstruction of surface tension (PROST) method (Renardy2002) uses the cost function and the search space in order to accurately compute the interface curvature for use in their numerical surface tension model. They report the elimination of spurious currents thanks to a balanced surface tension formulation as well as the accurate curvature obtained from the parabolic reconstruction.

We are however quite satisfied with the curvature obtained from the GHF method (Popinet2009), and therefore propose to use this curvature in combination with the search space. This method is referred to as the parabolic LVIRA (PLVIRA) method.

In fig. 25 we show the cost function for both search spaces as well as the optimal angles and also the error measured in terms of the symmetric difference eq. 40. We note that approximates quite well. In fig. 25 we show the resulting reconstructed interfaces for each of the optimal angles. As expected, the parabolic interface yields a significant increase in approximation accuracy of the interface.

### 5.3 PMoF

Parabolic reconstruction using the MoF cost function

is less straightforward. The reason for this is that the amount of information that is available is rather small: we have two degrees of freedom (

) and only the first moment as reference data (recall that the volume fraction uniquely defines the shift ). The following example illustrates that the MoF cost function cannot be used in combination with the search space .

###### Example 1.

Consider an interface defined by the level set in the unit control volume , where . The liquid volume fraction can then be computed as

 Ml0,c=∫12−12∫s−κ2x2−12dydx=12+s−κ24. (60)

Similarly we compute the first moment

 M1(cl)=∫12−12∫s−κ2x2−12xdydx=12[0κ2720+Ml0,c(Ml0,c−1)], (61)

where we substituted . It follows that for this example the sign of the curvature does not affect the first moment, resulting in at least two distinct solutions to the optimisation problem eq. 36.

Hence for using the MoF method with a parabolic reconstruction we must either decrease the dimensionality of the search space, or add more information such as higher order moments or some information from neighbouring control volumes.

We choose for the former option: the parabolic MoF (PMoF) method refers to the combination of and . Note that we could also interpret this as the latter option, since we include information from neighbouring control volumes indirectly via the GHF curvature.

In fig. 25 we show the cost function for both search spaces as well as the optimal angles and also the error measured in terms of the symmetric difference. The resulting reconstructed interfaces are shown in fig. 25. We find, for this example, that even though the PMoF method yields similar accuracy in terms of the cost function when compared to the MoF method, the interface approximation is significantly improved to the point where it can hardly be visually distinguished from the exact interface, which is a circular arc.

### 5.4 An algorithm for intersecting a polygon with a parabola

The evaluation of both cost functions, as well as the computation of the right-hand side of eq. 9 require the computation of either the zeroth or first moment of the intersection of some polygon with a parabola555The PROST method (Renardy2002) computes this intersection approximately. However since we consider 2D only, the intersection can still rather easily be computed exactly. defined as the zero level set of some . That is, we want to compute

 M0(P∩l(p)),M1(P∩l(p)). (62)

The description of the algorithm is limited to to the case . This is not restrictive, since if corresponds to an interface with positive curvature then we instead compute the moments using

 M0(P∩l(p))=M0(P)−M0(P∩(R2∖l(p)))=M0(P)−M0(P∩l(−p)), (63)

where again corresponds to an interface of negative curvature.

Let the vertices of the polygon be denoted by , for . Using the level set function we can determine for each vertex whether it is below () or above () the parabola, and therefore decide whether the vertex is retained or removed respectively. For each edge we can then determine to which of the following cases the edge belongs (see also fig. 29).