equationEquationequations \creflabelformatequation#2(#1)#3 \crefnamepluraleqEquationsequations \creflabelformatpluraleq#2(#1)#3 \crefnamesectionSectionsections \crefnamechapterChapterchapters \crefnamefigureFigurefigures \crefnametableTabletables \crefnamesubfigureFigurefigures \creflabelformatsubfigure#2#1#3
Calculation of interface curvature with the level-set 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] NO-7491 Trondheim, Norway
[-1.0mm] e–mail: firstname.lastname@example.org
Summary The level-set method is a popular method for interface capturing. One of the advantages of the level-set method is that the curvature and the normal vector of the interface can be readily calculated from the level-set function. However, in cases where the level-set method is used to capture topological changes, the standard discretization techniques for the curvature and the normal vector do not work properly. This is because they are affected by the discontinuities of the signed-distance function half-way between two interfaces. This article addresses the calculation of normal vectors and curvatures with the level-set method for such cases. It presents a discretization scheme that is relatively easy to implement in to an existing code. The improved discretization scheme is compared with a standard discretization scheme, first for a case with no flow, then for a case where two drops collide in a shear flow. The results show that the improved discretization yields more robust calculations in areas where topological changes are imminent.
The level-set method was introduced by Osher and Sethian . 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. .
This article addresses the calculation of interface geometries with the level-set 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 level-set function. These calculations are typically done with standard finite-difference methods. Since the level-set function is chosen to be a signed-distance 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 level-set method, see \creffig:cdroplets. The derivative of the level-set 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 finite-difference methods. For convenience the areas where the derivative of the level-set function is not defined will hereafter be refered to as kinks.
To the authors knowledge, this issue was first described in , where the level-set 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 . This was further expanded to include the normal vectors .
An alternative method to avoid the kinks is presented in , where a level-set extraction technique is presented. This technique uses an extraction algorithm to reconstruct separate level-set functions for each distinct body.
Accurate calculation of the curvature is important in many applications, in particular in curvature-driven 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  use a coupled level-set and volume-of-fluid method based on a fixed Eulerian grid, and they use a height function to calculate the curvatures. In  a refined level-set grid method is used to study two-phase flows on structured and unstructured grids for the flow solver. An interface-projected curvature-evaluation method is presented to achieve converging calculation of the curvature. In  they adopt a discontinuous Galerkin method and a pressure-stabilized finite-element method to solve the level-set equation and the Navier-Stokes equations, respectively. They develop a least-squares approach to calculate the normal vector and the curvature accurately, as opposed to using a direct derivation of the level-set function. This method is used in , where they show impressive results for simulation of turbulent atomization.
This article applies the level-set method to incompressible two-phase flow in two dimensions. The direction difference described in  is used to calculate the normal vectors, and a curvature discretization is presented which is based on the geometry-aware discretization given in . 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 level-set 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.
2 Governing equations
2.1 Navier-Stokes equations for two-phase flow
Consider a two-phase domain where and denote the regions occupied by the respective phases. The domain is divided by an interface as illustrated in \creffig:domain. The governing equations for incompressible and immiscible two-phase flow in the domain with an interface force on the interface can be stated as
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 Navier-Stokes equations for incompressible two-phase 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,
where is the tangent vector along the interface and denotes the jump across an interface, that is
2.2 Level-set method
The interface is captured with the zero level set of the level-set function , which is prescribed as a signed-distance function. That is, the interface is given by
and for any ,
The interface is updated by solving an advection equation for ,
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
to steady state, c.f. . Here is pseudo-time and is a smeared sign function which is equal to zero at the interface,
When \crefeq:lseq is solved numerically, the level-set function loses its signed-distance property due to numerical dissipation. The level-set function is therefore reinitialized regularly by solving
to steady state as proposed in . Here is the level-set function that needs to be reinitialized.
One of the advantages of the level-set method is that normal vectors and curvatures can be readily calculated from the level-set function, i.e.
3 Numerical methods
The Navier-Stokes equations, (1) and (2), 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 second-order central difference scheme, except for the convective terms which are discretized by a fifth-order WENO scheme. The temporal discretization is done with explicit strong stability-preserving Runge-Kutta (SSP RK) schemes, see . A three-stage third-order SSP-RK method is used for the Navier-Stokes equations (1) and (2), and a four-stage second-order SSP-RK method is used for the level-set equations (10), (11) and (13).
The method presented in  is used to improve the computational speed. The method is often called the narrow-band method, since the level-set 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 Ghost-Fluid 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 right-hand side of the Poisson equation for the pressure. Consider a one-dimensional case where and where the interface lies between and . In this case,
where is the general right-hand side value and is the curvature at the interface. The sign of the added term depends on the sign of . See  for more details on how the GFM is used for the Navier-Stokes equations and  for a description on how to use the GFM for a variable-coefficient Poisson equation.
The normal vector and the curvature defined by \crefeq:norm,eq:curv are typically discretized by the second-order 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 \crefeq:pressure-poisson,
If the level-set method is used to capture non-trivial geometries, the level-set function will in general contain areas where it is not smooth, i.e. kinks. This is depicted in \creffig:level-set-function, which shows a level-set function in a one-dimensional 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 level-set function. If the level-set 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 narrow-band approach.
4 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 level-set 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.
4.1 The normal vector
A discretization scheme is presented in  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 one-sided differences based on the values of a quality function,
which is approximated with central differences. The quality function effectively detects the areas where the level-set function differs from the signed-distance 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,
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  for more details.
The direction difference is now defined as
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 one-sided 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.
4.2 The curvature
The curvature is calculated with a discretization that is based on the improved geometry-aware curvature discretization presented by Macklin and Lowengrub . This is a method where the curvature is calculated at the interfaces directly with the use of a least-squares curve parametrization of the interface. The curve parametrization is used to create a local level-set function from which the curvature is calculated using standard discretization techniques. The local level-set 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 \crefeq:pressure-poisson directly, whereas the present method calculates and with the improved curvature discretization and then use linear interpolation as described in \crefeq:curvature-gamma 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 Navier-Stokes equations. An example of such a case is when the curvature is used to model interfacial flows with surfactant .
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 least-square parametrization used in  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 level-set function based on the parametrization .
Use the standard discretization on the local level-set function to calculate the curvature.
5 Details of the curvature algorithm
5.1 Locating the closest interface
A breadth-first search is used to to identify the closest interface, see \creffig:first-point-2. 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,
In the diagonal cases the interface point is found with bilinear interpolation along the diagonal. This leads to
where is the solution of
The values depend on the grid cell. For instance, when searching along the diagonal between and the values will be
5.2 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 left-hand 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 marching-squares
In some rare cases the algorithm must handle the ambiguous case depicted in \creffig:marching-square-amb. In these cases there are four interface crossing-points 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 level-set method is not able to resolve the interfaces on a sub-cell 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.
5.3 Curve fitting
Cubic Hermite splines are used to fit a curve to the set of points
Let the curve parametrization be denoted for . A cubic spline is a parametrization where
and each interpolant is a third-order 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
for , where is the curve tangents and and are Hermite basis polynomials,
The choice of the tangents is non-unique, 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 level-set 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 . 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,
for . Next the tangents are initialized as the average of the secants at every point,
for . The curve tangents at the endpoints are set to and . Finally let pass from 1 through and set where , and where .
5.4 Local level-set function
The local level-set 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 level-set function, .
A precise definition of is
where is the signed-distance 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
When is composed of cubic polynomials as is the case for cubic Hermite splines, the computation of the distance requires the solution of several fifth-order 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 Newton-Raphson and bisection method is used to refine them. The correct sign is found by solving
where is the tangent vector of .
6 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 second-order central differences are used as the standard discretization.
6.1 A static disc above a rectangle
The first case is a simple and static test-case where a disc of radius is positioned at a distance above a rectangle, see \creffig:circleandline. Only the level-set function and the geometrical quantities are considered. This means that none of the governing equations are solved (\crefeq:masseq,eq:momeq,eq:lseq,eq:lsvelext,eq:lsreinit). 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 .
fig:simple-normal-vectors 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.
fig:simple-curvatures 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 \crefeq:curvature-gamma. The effect of these errors will become more clear in the next case.
6.2 Drop collision in shear flow
The second case considers drop collision in shear flow, see \creffig:drop-in-shear. 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 curve-fitting based discretization is properly revealed.
The purpose of this case is to study the behaviour of the level-set 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
In the following results the Reynolds and the Capillary numbers were set to
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 .
fig:ycf84 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. \crefeq:pressure-poisson. \creffig:ycf84-pressure shows the pressure field at . It can be seen that the pressure field for the standard discretization is distorted in the thin-film 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 thin-film region and lower at the edges. The pressure change induces a flow out of the region which is more as expected.
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 . The curvature was discretized with a scheme that is based on the geometry-aware discretization presented in . 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.
This work was financed through the Enabling Low-Emission 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.
- The marching-squares algorithm is an equivalent two-dimensional formulation of the well known marching-cubes algorithm presented in . The algorithm was mainly developed for use in computer graphics.
- D.Adalsteinsson and J. A.Sethian A fast level set method for propagating interfaces Journal of Computational Physics, vol.118, 269–277, 1995.
- 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.
- F. N.Fritsch and R. E.Carlson Monotone piecewise cubic interpolation SIAM Journal of Numerical Analysis, vol.17(2), 238–246, 1980.
- S.Gottlieb, C. W.Shu and E.Tadmor Strong stability-preserving high-order time discretization methods SIAM Review, vol.43, 89–112, 2001.
- 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 82-471-7318-2.
- M.Herrmann A balanced force refined level set grid method for two-phase flows on unstructured flow solver grids Journal of Computational Physics, vol.227, 2674–2706, 2008.
- 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.
- S.Lang Algebra Graduate texts in mathematics. Springer, 2002.
- 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.
- 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.
- P.Macklin and J.Lowengrub Evolving interfaces via gradients of geometry-dependent interior poisson problems: Application to tumor growth Journal of Computational Physics, vol.203, 191–220, 2005.
- P.Macklin and J.Lowengrub An improved geometry-aware curvature discretziation for level set methods: Application to tumor growth Journal of Computational Physics, vol.215, 392–401, 2006.
- 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.
- 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.
- S.Osher and R. P.Fedkiw The Level-Set Method and Dynamic Implicit Surfaces Springer, 2003.
- S.Osher and J. A.Sethian Fronts propagating with curvature dependent speed: Algorithms based on Hamilton-Jacobi formulations Journal of Computational Physics, vol.79, 12–49, 1988.
- H.Prautzsh, W.Boehm and M.Paluszny Bézier and B-spline Techniques Springer, 2002.
- D.Salac and W.Lu A local semi-implicit level-set method for interface motion Journal of Scientific Computing, vol.35, 330–349, 2008.
- J. A.Sethian and P.Smereka Level set methods for fluid interfaces Annual Review of Fluid Mechanics, vol.35, 341–372, 2003.
- M.Sussman, P.Smereka and S.Osher A level set approach for computing solutions to incompressible two-phase flow Journal of Computational Physics, vol.114, 146–159, 1994.
- B.Waerden, E.Artin and E.Noether Algebra Number v. 1 in Algebra. Springer-Verlag, 2003.
- Z.Wang and A. Y.Tong A sharp surface tension modeling method for two-phase incompressible interfacial flows International Journal for Numerical Methods in Fluids, vol.64, 709–732, 2010.
- 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.
- 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.