Calculation of interface curvature with the levelset method
Karl Yngve Lervåg
[0mm] Norwegian University of Science and Technology
[1.0mm] Department of Energy and Process Engineering
[1.0mm] Kolbjørn Hejes veg 2
[1.0mm] NO7491 Trondheim, Norway
[1.0mm] e–mail: karl.y.lervag@ntnu.no
[1.0mm]
Introduction
The levelset method was introduced by Osher and Sethian [16]. It is designed to implicitly track moving interfaces through an isocontour of a function defined in the entire domain. In particular, it is designed for problems in multiple spatial dimensions in which the topology of the evolving interface changes during the course of events, c.f. [19].
This article addresses the calculation of interface geometries with the levelset method. This method allows us to calculate the normal vector and the curvature of an interface directly as the first and second derivatives of the levelset function. These calculations are typically done with standard finitedifference methods. Since the levelset function is chosen to be a signeddistance function, in most cases it will have areas where it is not smooth. Consider for instance two colliding droplets where the interfaces are captured with the levelset method, see Figure 1. The derivative of the levelset function will not be defined at the points outside the droplets that have an equal distance to both droplets. When the droplets are in near contact, this discontinuity in the derivative will lead to significant errors when calculating the interface geometries with standard finitedifference methods. For convenience the areas where the derivative of the levelset function is not defined will hereafter be refered to as kinks.
To the authors knowledge, this issue was first described in [11], where the levelset method was used to model tumor growth. Here Macklin and Lowengrub presented a direction difference to treat the discretization across kinks for the normal vector and the curvature. They later presented an improved method where curve fitting was used to calculate the curvatures [12]. This was further expanded to include the normal vectors [13].
An alternative method to avoid the kinks is presented in [18], where a levelset extraction technique is presented. This technique uses an extraction algorithm to reconstruct separate levelset functions for each distinct body.
Accurate calculation of the curvature is important in many applications, in particular in curvaturedriven flows. There are several examples in the literature of methods that improve the accuracy of the curvature calculations, but that do not consider the problem with the discretization across the kinks. The authors in [22] use a coupled levelset and volumeoffluid method based on a fixed Eulerian grid, and they use a height function to calculate the curvatures. In [6] a refined levelset grid method is used to study twophase flows on structured and unstructured grids for the flow solver. An interfaceprojected curvatureevaluation method is presented to achieve converging calculation of the curvature. In [14] they adopt a discontinuous Galerkin method and a pressurestabilized finiteelement method to solve the levelset equation and the NavierStokes equations, respectively. They develop a leastsquares approach to calculate the normal vector and the curvature accurately, as opposed to using a direct derivation of the levelset function. This method is used in [2], where they show impressive results for simulation of turbulent atomization.
This article applies the levelset method to incompressible twophase flow in two dimensions. The direction difference described in [11] is used to calculate the normal vectors, and a curvature discretization is presented which is based on the geometryaware discretization given in [12]. The main advantage of the present scheme is that is is relatively easy to implement, since it requires very little change to a typical implementation of the levelset method.
The article starts by briefly describing the governing equations. It continues with a description of the numerical methods that are used. Then the discretization schemes for the normal vector and the curvature are presented, followed by a detailed description of the curvature discretization. Next a comparison of the improved discretization and the standard discretization is made, first on static interfaces in near contact, then on two drops colliding in a shear flow. Finally some concluding remarks are made.
Governing equations
NavierStokes equations for twophase flow
Consider a twophase domain where and denote the regions occupied by the respective phases. The domain is divided by an interface as illustrated in Figure 2. The governing equations for incompressible and immiscible twophase flow in the domain with an interface force on the interface can be stated as
(0)  
(0) 
where is the velocity vector, is the pressure, is the specific body force, is the coefficient of surface tension, is the curvature, is the normal vector which points to , is the Dirac Delta function, is a parametrization of the interface, is the density and is the viscosity. These equations are often called the NavierStokes equations for incompressible twophase flow.
It is assumed that the density and the viscosity are constant in each phase, but they may be discontinuous across the interface. The interface force and the discontinuities in the density and the viscosity lead to a set of interface conditions,
(0)  
(0)  
(0)  
(0) 
where is the tangent vector along the interface and denotes the jump across an interface, that is
(0) 
See [7, 5] for more details and a derivation of the interface conditions.
Levelset method
The interface is captured with the zero level set of the levelset function , which is prescribed as a signeddistance function. That is, the interface is given by
(0) 
and for any ,
(0) 
The interface is updated by solving an advection equation for ,
(0) 
where is the velocity at the interface extended to the entire domain. The interface velocity is extended from the interface to the domain by solving
(0) 
to steady state, c.f. [24]. Here is pseudotime and is a smeared sign function which is equal to zero at the interface,
(0) 
When Levelset method is solved numerically, the levelset function loses its signeddistance property due to numerical dissipation. The levelset function is therefore reinitialized regularly by solving
(0) 
to steady state as proposed in [20]. Here is the levelset function that needs to be reinitialized.
One of the advantages of the levelset method is that normal vectors and curvatures can be readily calculated from the levelset function, i.e.
(0)  
(0) 
Numerical methods
The NavierStokes equations, (NavierStokes equations for twophase flow) and (NavierStokes equations for twophase flow), are solved by a projection method on a staggered grid as described in [5, Chapter 5.1.1]. The spatial terms are discretized by the secondorder central difference scheme, except for the convective terms which are discretized by a fifthorder WENO scheme. The temporal discretization is done with explicit strong stabilitypreserving RungeKutta (SSP RK) schemes, see [4]. A threestage thirdorder SSPRK method is used for the NavierStokes equations (NavierStokes equations for twophase flow) and (NavierStokes equations for twophase flow), and a fourstage secondorder SSPRK method is used for the levelset equations (Levelset method), (Levelset method) and (Levelset method).
The method presented in [1] is used to improve the computational speed. The method is often called the narrowband method, since the levelset function is only updated in a narrow band across the interface at each time step.
The interface conditions are treated in a sharp fashion with the GhostFluid Method (GFM), which incorporates the discontinuities into the discretization stencils by altering the stencils close to the interfaces. For instance, the GFM requires that a term is added to the stencil on the righthand side of the Poisson equation for the pressure. Consider a onedimensional case where and where the interface lies between and . In this case,
(0) 
where is the general righthand side value and is the curvature at the interface. The sign of the added term depends on the sign of . See [7] for more details on how the GFM is used for the NavierStokes equations and [9] for a description on how to use the GFM for a variablecoefficient Poisson equation.
The normal vector and the curvature defined by Levelset method and Levelset method are typically discretized by the secondorder central difference scheme, c.f [7, 19, 23]. The curvatures are calculated on the grid nodes and then interpolated with simple linear interpolation to the interface, e.g. for in Numerical methods,
(0) 
If the levelset method is used to capture nontrivial geometries, the levelset function will in general contain areas where it is not smooth, i.e. kinks. This is depicted in Figure 3, which shows a levelset function in a onedimensional domain that captures two interfaces, one on each side of the grid point . The kink at will lead to potentially large errors with the standard discretization both for the curvature and the normal vector. The errors in the curvature will lead to erroneous pressure jumps at the interfaces, and the errors in the normal vector affects both the discretized interface conditions and the advection of the levelset function. If the levelset method is used to study for example coalescence and breakup of drops, these errors may severly affect the simulations.
It should be noted that the kinks that appear far from any interfaces are handled by ensuring that the denominators do not become zero, as explained in [15, Sections 2.3 to 2.4]. This works fine, since only the values of the curvature at the grid nodes adjacent to any interface are used. Also, the normal vector only needs to be accurate close to the interface due to the narrowband approach.
Improved discretization of geometrical quantities
The previous section explained why it is necessary to develop new discretization schemes for the normal vector and the curvature that can handle kinks in the levelset function. This section will give a brief presentation of a better discretization scheme for the normal vector and an overview of an algorithm to calculate the curvature.
The normal vector
A discretization scheme is presented in [11] which uses a quality function to ensure that the differences never cross kinks. The basic strategy is to use a combination of central differences and onesided differences based on the values of a quality function,
(0) 
which is approximated with central differences. The quality function effectively detects the areas where the levelset function differs from the signeddistance function. Let and , then can be used to detect kinks. The parameter is tuned such that the quality function will detect all the kinks. The value is used in the present work.
The quality function is used to define a direction function,
(0) 
where
(0) 
is defined in a similar manner. If or is undetermined, is chosen as the vector normal to . It is normalized, and the sign is chosen such that it points in the direction of best quality. See [11] for more details.
The direction difference is now defined as
(0) 
where is a piecewise smooth function. The normal vector is calculated using the direction difference on , which is equivalent to using central differences in smooth areas and onesided differences in areas close to the kinks. This method makes sure that the differences do not cross any kinks, and the normal vector can be accurately calculated even close to a kink.
The curvature
The curvature is calculated with a discretization that is based on the improved geometryaware curvature discretization presented by Macklin and Lowengrub [12]. This is a method where the curvature is calculated at the interfaces directly with the use of a leastsquares curve parametrization of the interface. The curve parametrization is used to create a local levelset function from which the curvature is calculated using standard discretization techniques. The local levelset function only depends on one interface and is therefore free of kinks.
The main difference between the present method and that of Macklin and Lowengrub is that they calculate the curvature at the interface directly, whereas the present method instead calculates the curvature at the grid nodes. In other words, Macklin and Lowengrub calculate in Numerical methods directly, whereas the present method calculates and with the improved curvature discretization and then use linear interpolation as described in Numerical methods to find . The main motivation behind this difference is that the present method does not require a significant change to any existing code. Thus it is relatively straightforward to implement the present method even when the curvature is needed for more than the Capillary force term in the NavierStokes equations. An example of such a case is when the curvature is used to model interfacial flows with surfactant [23].
An important consequence of the previously explained difference is that it becomes more important to have an accurate representation of the interface. The curvature discretization presented here uses monotone cubic Hermite splines to parametrize the curve. The leastsquare parametrization used in [12] is only accurate very close to the point where the curvature needs to be calculated. The Hermite spline is more accurate along the entire interface representation.
The algorithm to calculate the curvature at can be summarized as follows. The details are explained in the next section.

If , where and , then it is safe to use the standard discretization. Otherwise continue to the next step.

Locate the closest interface, .

Find a set of points .

Create a parametrization of the points .

Calculate a local levelset function based on the parametrization .

Use the standard discretization on the local levelset function to calculate the curvature.
Details of the curvature algorithm
Locating the closest interface
A breadthfirst search is used to to identify the closest interface, see Figure 4. Let denote the starting point and denote the desired point on the closest interface. The search iterates over all the eight edges from to its neighbours and tries to locate an interface which is identified by a change of sign of . If more than one interface is found, is chosen to be the point that is closest to . If no interfaces are located the search continues at the next depth. The search continues in this manner until an interface is found. Note that this algorithm does not in general return the point on the interface which is closest to .
The crossing points between the grid edges and the interfaces are found with linear and bilinear interpolation. E.g. if an interface crosses the edge between and at , the interface point is found by linear interpolation,
(0) 
where
(0) 
In the diagonal cases the interface point is found with bilinear interpolation along the diagonal. This leads to
(0) 
where is the solution of
(0) 
The values depend on the grid cell. For instance, when searching along the diagonal between and the values will be
(0)  
(0)  
(0) 
Searching for points on an interface
When an interface and a corresponding point on the interface are found, the next step is to find a set of points on the same interface. The points should be ordered such that when traversing the points from to , the phase with is on the lefthand side. Note that the ordering of the points may be done after all the points are found. Three criteria are used when searching for new points:

The points are located on the grid edges.

The distance between and for all is greater than a given threshold .

The points that are closest to are selected, where is the initial point where the curvature is to be calculated.
Let be given. To find a new point on , a variant of the marchingsquares algorithm^{1}^{1}1The marchingsquares algorithm is an equivalent twodimensional formulation of the well known marchingcubes algorithm presented in [10]. The algorithm was mainly developed for use in computer graphics. is used. Given and a search direction which is either clockwise or counter clockwise, the algorithm searches for all the points where an interface crosses the edges of the mesh rectangle . In most cases there will be two such points and is one of them. is then selected based on the search direction. If , the search is continued in the adjacent mesh rectangle. The search process is depicted in Figure 5(a).
In some rare cases the algorithm must handle the ambiguous case depicted in Figure 5(b). In these cases there are four interface crossingpoints and two solutions. Either solution is valid, and it is not possible to say which solution is better. The current implementation selects the first solution that it finds, which will be in all practical sense a random choice. Note that the ambiguous cases only occur when two interfaces cross a single grid cell. The ambiguity comes from the fact that the levelset method is not able to resolve the interfaces on a subcell resolution.
It was found that points where necessary in order to ensure that the closest points on the interface with respect to the different grid points are captured with the spline parametrization.
Curve fitting
Cubic Hermite splines are used to fit a curve to the set of points
(0) 
Let the curve parametrization be denoted for . A cubic spline is a parametrization where
(0) 
where
(0) 
and each interpolant is a thirdorder polynomial. A Hermite spline is a spline where each interpolant is in Hermite form, see [17, Chapter 4.5]. The interpolants are created by solving the equations
(0) 
for , where is the curve tangents and and are Hermite basis polynomials,
(0) 
The choice of the tangents is nonunique, and there are several possible options for a cubic Hermite spline.
It is essential that the spline is properly oriented. This is because we require to find both the distance and the position of a point on the grid relative to the spline in order to define a local levelset function. The orientation of the spline is defined such that when increases, is to the left.
To ensure that our curve is properly oriented, the tangents are chosen as described in [3]. This will ensure monotonicity for each component as long as the input data is monotone. The tangents are modified as follows. First the slopes of the secant lines between successive points are computed,
(0) 
for . Next the tangents are initialized as the average of the secants at every point,
(0) 
for . The curve tangents at the endpoints are set to and . Finally let pass from 1 through and set where , and where .
Local levelset function
The local levelset function, here denoted as , is calculated at the grid points surrounding and including . The curvature is then calculated with the standard discretization stencil where is used instead of the global levelset function, .
A precise definition of is
(0) 
where is the signeddistance function, which is negative in phase one and positive in phase two. This function is calculated by first finding the minimum distance between and and then deciding the correct sign. The minimum distance is found by minimizing the norm
(0) 
When is composed of cubic polynomials as is the case for cubic Hermite splines, the computation of the distance requires the solution of several fifthorder polynomial equations. Sturm’s method (see [21, Section 11.3] or [8, Chapter XI,§2]) is employed to locate and bracket the solutions and a combined NewtonRaphson and bisection method is used to refine them. The correct sign is found by solving
(0) 
where is the tangent vector of .
Verification and testing
This section presents results of calculating normal vectors and curvatures with the improved discretization schemes. The results are compared with the standard discretization. Note that in both the following cases the standard secondorder central differences are used as the standard discretization.
A static disc above a rectangle
The first case is a simple and static testcase where a disc of radius is positioned at a distance above a rectangle, see Figure 6. Only the levelset function and the geometrical quantities are considered. This means that none of the governing equations are solved (Levelset method, Levelset method, Levelset method, NavierStokes equations for twophase flow and NavierStokes equations for twophase flow). When is small, the kinks along the dotted line will affect the discretization stencils as has been explained.
The parameters for this case is and . The domain is , and the straight line is positioned at . The grid size is .
Figure 7 shows a comparison of the calculated normal vectors. The standard discretization is depicted with red vectors and the direction difference is depicted with green vectors. The figure shows that the standard discretization yields much less accurate results along the kink region than the direction difference.
Figure 8 shows a comparison of the calculated curvatures. Note that the curvature is only calculated at grid points adjacent to the interfaces. At the grid points where it is not calculated, it is set to zero. The figure shows that the standard discretization leads to large errors in the calculated curvatures in the areas that are close to two interfaces. In particular note that the sign of the curvature becomes wrong. The analytic curvature for this case is , and the curvature spikes seen for the standard discretization is in the order of . These spikes will lead to large errors in the pressure jumps through Numerical methods. The effect of these errors will become more clear in the next case.
Drop collision in shear flow
The second case considers drop collision in shear flow, see Figure 9. The drops both have a radius , and they are initially placed at a distance apart in a shear flow where the flow velocity changes linearly from at the bottom wall to at the top wall. The computational domain is , and the grid size is . The size of the grid is chosen to be relatively coarse, such that the difference between the standard discretization of the curvature and the curvefitting based discretization is properly revealed.
The purpose of this case is to study the behaviour of the levelset method, in particular the calculation of the curvatures, when the drops are in close proximity. It is therefore a natural simplification to only consider the case where the density difference and the viscosity difference of the phases are zero, i.e. there is no jump in density or viscosity across the interface.
The flow is governed by the Reynolds number and the Capillary number, which in the current case can be defined by
(0)  
(0) 
In the following results the Reynolds and the Capillary numbers were set to
(0) 
The choice was made such that the drops would not be severely deformed in the shear flow. The radius of the drops was , and the distance from the center line to the drop centers was .
Figure 10 shows a comparison of the interface evolution and the curvature between the standard discretization and the improved discretization. The top and bottom rows show the evolution for the standard discretization and the improved discretization, respectively. The kinks between the drops lead to curvature spikes with the standard discretization, whereas the improved discretization calculates the curvature along the kink in a much more reliable manner. The curvature spikes are seen to prevent coalescence. This is due to the effect they have on the pressure field.
The errors in the curvature with the standard discretization lead to an erroneous pressure field between the drops that prevents coalescence, c.f. Numerical methods. Figure 11 shows the pressure field at . It can be seen that the pressure field for the standard discretization is distorted in the thinfilm region. This distortion in the pressure leads to a flow in the film region which suppresses coalescence. The corresponding result for the improved method shows that the pressure is not distorted. It is high in the center of the thinfilm region and lower at the edges. The pressure change induces a flow out of the region which is more as expected.
Conclusions
This article has implemented improved discretization schemes for the normal vector and the curvature of the interface between two phases. The normal vector was discretized by the direction difference which is presented in [11]. The curvature was discretized with a scheme that is based on the geometryaware discretization presented in [12]. The main advantage of the present discretization method for the curvature is that it is relatively straightforward to implement in to an existing code since it does not require a change of the existing framework.
The implementation of the curvature discretization have been described in detail. The improved schemes are compared with the standard discretization in two different cases. The first case is a direct comparison of the schemes for a case with no flow. The second case compares the evolution of two drops colliding in shear flow. Both tests demonstrate that the standard discretization of the normal vector and the curvature leads to erroneous behaviour at the kink locations. The second case shows that this behaviour prevents coalescence from occurring due to an erroneous pressure field. The curvature spikes at the kink regions are not observed with the improved discretization schemes, and coalescence is achieved for the second case.
Acknowledgements
This work was financed through the Enabling LowEmission LNG Systems project, and the author acknowledge the contributions of GDF SUEZ, Statoil and the Petromaks programme of the Research Council of Norway (193062/S60).
The author acknowledges Bernhard Müller (NTNU) and Svend Tollak Munkejord (SINTEF Energy Research) for valuable feedback on the manuscript. The author also acknowledges Leif Amund Lie and Eirik Svanes for several good discussions.
References
 [1] D.Adalsteinsson and J. A.Sethian A fast level set method for propagating interfaces Journal of Computational Physics, vol.118, 269–277, 1995.
 [2] O.Desjardins, V.Moureau and H.Pitsch An accurate conservative level set/ghost fluid method for simulating turbulent atomization Journal of Computational Physics, vol.227, 8395–8416, 2008.
 [3] F. N.Fritsch and R. E.Carlson Monotone piecewise cubic interpolation SIAM Journal of Numerical Analysis, vol.17(2), 238–246, 1980.
 [4] S.Gottlieb, C. W.Shu and E.Tadmor Strong stabilitypreserving highorder time discretization methods SIAM Review, vol.43, 89–112, 2001.
 [5] E. B.Hansen Numerical Simulation of Droplet Dynamics in the Presence of an Electric Field Doctoral thesis, Norwegian University of Science and Technology, Department of Energy and Process Engineering, Trondheim, Nov. 2005 ISBN 8247173182.
 [6] M.Herrmann A balanced force refined level set grid method for twophase flows on unstructured flow solver grids Journal of Computational Physics, vol.227, 2674–2706, 2008.
 [7] M.Kang, R. P.Fedkiw and X.D.Liu A boundary condition capturing method for multiphase incompressible flow Journal of Scientific Computing, vol.15(3), 323–360, 2000.
 [8] S.Lang Algebra Graduate texts in mathematics. Springer, 2002.
 [9] X.D.Liu, R. P.Fedkiw and M.Kang A boundary condition capturing method for poisson’s equation on irregular domains Journal of Computational Physics, vol.160, 151–178, 2000.
 [10] W. E.Lorensen and H. E.Cline Marching cubes: A high resolution 3D surface construction algorithm Computer Graphics, vol.21(4), 163–169, July 1987.
 [11] P.Macklin and J.Lowengrub Evolving interfaces via gradients of geometrydependent interior poisson problems: Application to tumor growth Journal of Computational Physics, vol.203, 191–220, 2005.
 [12] P.Macklin and J.Lowengrub An improved geometryaware curvature discretziation for level set methods: Application to tumor growth Journal of Computational Physics, vol.215, 392–401, 2006.
 [13] P.Macklin and J. S.Lowengrub A new ghost cell/level set method for moving boundary problems: Application to tumor growth Journal of Scientific Computing, vol.35, 266–299, 2008.
 [14] E.Marchandise, P.Geuzaine, N.Chevaugeon and J.F.Remacle A stabilized finite element method using a discontinuous level set approach for the computation of bubble dynamics Journal of Computational Physics, vol.225, 949–974, 2007.
 [15] S.Osher and R. P.Fedkiw The LevelSet Method and Dynamic Implicit Surfaces Springer, 2003.
 [16] S.Osher and J. A.Sethian Fronts propagating with curvature dependent speed: Algorithms based on HamiltonJacobi formulations Journal of Computational Physics, vol.79, 12–49, 1988.
 [17] H.Prautzsh, W.Boehm and M.Paluszny Bézier and Bspline Techniques Springer, 2002.
 [18] D.Salac and W.Lu A local semiimplicit levelset method for interface motion Journal of Scientific Computing, vol.35, 330–349, 2008.
 [19] J. A.Sethian and P.Smereka Level set methods for fluid interfaces Annual Review of Fluid Mechanics, vol.35, 341–372, 2003.
 [20] M.Sussman, P.Smereka and S.Osher A level set approach for computing solutions to incompressible twophase flow Journal of Computational Physics, vol.114, 146–159, 1994.
 [21] B.Waerden, E.Artin and E.Noether Algebra Number v. 1 in Algebra. SpringerVerlag, 2003.
 [22] Z.Wang and A. Y.Tong A sharp surface tension modeling method for twophase incompressible interfacial flows International Journal for Numerical Methods in Fluids, vol.64, 709–732, 2010.
 [23] J.J.Xu, Z.Li, J.Lowengrub and H.Zhau A level set method for interfacial flows with surfactants Journal of Computational Physics, vol.212(2), 590–616, March 2006.
 [24] H. K.Zhao, T.Chan, B.Merriman and S.Osher A variational level set approach to multiphase motion Journal of Computational Physics, vol.127, 179–195, 1996.