1 Introduction
The advection of the phase interface plays a central role in the simulation of immiscible, and in our case incompressible, twophase 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 heightfunctions (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 YoungLaplace 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 timedependent 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 LHFbased 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 timedependent problems.
Throughout this paper we will denote the positive timestep 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
(1) 
where the timestep satisfies the CFL timestep restriction
(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(3) 
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
(4) 
The liquid volume and its control volume fraction are then defined as
(5) 
for some control volume . The liquid first moment and centroid are similarly denoted by and respectively.
2.1 Lagrangian remapping
figure [2] [0.48][][t] [0.48][][t]
We let denote the flow map of the velocity field over the time interval , that is, solves the initial value problem
(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
(7) 
By the invertibility of it holds that and therefore
(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
(9) 
Equation (9) is key to understanding how a Lagrangian remappingbased geometric VoF method can be constructed, since such a method approximates each of the terms on the righthand 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
(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.
figure [3] [0.31][][t] [0.31][][t] [0.31][][t]
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 secondorder accurate Heun method, which, as we will see, is sufficiently accurate for this purpose. This results in the following approximate liquid volume
(11) 
The method described here is very similar to the LagrangianEulerian 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(12) 
Using the symmetric difference we define the reconstruction error as
(13) 
Here denotes the liquid neighbourhood centred around the control volume
(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
(15) 
and the integration error
(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
(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
(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))
(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
(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 )
(23) 
∎
The result of theorem 1 is limited to a single timestep. 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 timestep convergence result, the stability of the interface reconstruction method w.r.t. perturbations in the reference moments (i.e. errors from the previous timestep) 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 secondorder 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 timedependent 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 heightfunction
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 coordinate axes^{1}^{1}1For 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
(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 secondorder 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
(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 wellresolved 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
(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 secondorder accuracy.
The fully approximated curvature is then given by
(28) 
where the secondorder finite difference approximations to the first and second derivatives of some function are defined as
(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 heightfunction (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 secondorder 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
(30) 
Proof.
Let the approximation error of the averaged LHF be defined as
(31) 
where by the assumed smoothness of . By the midpoint rule we find that . It follows that
(32) 
by the differentiability of .
Similarly, we let the approximation error due to the approximation of the volume fractions be defined as
(33) 
where we have defined
(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
(35) 
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.
4 Optimisationbased 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 optimisationbased reconstruction methods. That is, the reconstruction step for each control volume can be written as an optimisation problem
(36) 
where is the search space consisting of linear interfaces and is a cost function. More precisely, the search space is defined as
(37) 
where is the unit 1sphere 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 halfspace if is linear) is denoted by
(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
(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
(40) 
4.1 The least squares VoF interface reconstruction algorithm
The least squares VoF interface reconstruction algorithm (LVIRA) (puckett1991volume) uses the following cost function
(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 minimum^{2}^{2}2For 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
(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 solutions^{3}^{3}3For 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
(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)
(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
figure [3] [0.305][][t] [0.305][][t] [0.305][][t]
Steps 1 (Centroid advection).
The liquid first moment of the intersection of the approximate preimage with the approximate liquid domain is computed
(45) 
The centroid of this intersection is advected as a point particle using the linearly interpolated velocity field (in space and time)
(46) 
where for simplicity in presentation we consider here the forward Euler method for approximation of the timeintegration 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
(47) 
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
(49)  
(50)  
(51) 
It follows that if
(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
(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
(54) 
In Dyadechko2005 the following estimate is given for the approximation error made in the first moment when the MoF method is used
(55) 
which assumes that the interface that is being approximated is a curve with a bounded (from below) radius of curvature. After a single timestep however, the interface has become piecewise smooth, and for the general nonsmooth case the authors of (Dyadechko2005) provide the pessimistic estimate
(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 timestep, we can not prove consistency of the firstmoment when the MoF cost function is used.
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
(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 one^{4}^{4}4The ‘unrestricted’ search space includes a tangential shift and is given by
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.
(eq. 37)  (eq. 58)  

(eq. 41)  (E)LVIRA (puckett1991volume; Pilliod2004)  PLVIRA  PROST (Renardy2002) 
(eq. 42)  MoF (Dyadechko2005)  PMoF  N/A 
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
(60) 
Similarly we compute the first moment
(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 righthand side of eq. 9 require the computation of either the zeroth or first moment of the intersection of some polygon with a parabola^{5}^{5}5The 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
(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
(63) 
where again corresponds to an interface of negative curvature.
figure [3] [0.31][][c] [0.31][][c] [0.31][][c]
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).
Comments
There are no comments yet.