Many unfitted mesh (fictitious domain, immersed boundary, cut cell) discretization methods have been proposed in the literature to enforce interface conditions for applications ranging from the Poisson equation to multiphase flow models and fluid-structure interaction problems. In this work, we focus on weak formulations in which interfacial terms are represented by surface or volume integrals. In level set methods osher ; sethian_review , the interface is defined as a manifold on which an approximate signed distance function (SDF) becomes equal to zero. The SDF property of advected level sets is usually preserved using a postprocessing procedure known as redistancing. The monolithic conservative version proposed in mono incorporates redistancing and mass correction into a nonlinear transport equation for the level set function. A variety of unfitted finite element methods (FEM) are available for numerical integration over sharp, diffuse, or surrogate interfaces. A sharp interface method enriches the local finite element space to accommodate boundary/jump conditions on an embedded surface. For example, this strategy is used in XFEM xfem , unfitted Nitsche FEM hansbo2 ; wadbro , and CutFEM cutfem ; hansbo2014 approaches. Integration over sharp interfaces has the potential drawback that small cut cells may cause ill-conditioning or severe time step restrictions. Possible remedies to this problem include the use of stabilization terms hansbo2014 ; zahedi , ghost penalties burman ; olsh , and (explicit-)implicit schemes for time integration may .
Diffuse interface methods replace surface integrals by volume integrals depending on regularized Dirac delta functions kublik2016 ; kublik2018 ; teigen2011 ; teigen2009 . This family of unfitted FEM is closely related to classical immersed boundary methods mittal ; peskin2002 and well suited, e.g., for implementation of surface tension effects in two-phase flow models CSF ; hysing . The construction of approximate delta functions for level set methods was addressed in engquist ; smereka ; towers ; zahedi . Müller et al. mueller2013 devised an elegant alternative, which uses the divergence theorem and divergence-free basis functions to reduce integration over an embedded interface to that over a fitted boundary. Another promising new approach is the shifted boundary method (SBM) sbm ; sbm1 ; sbm2 , which belongs to the family of surrogate interface approximations. To avoid the need for dealing with small cut cells, the SBM applies surrogate boundary conditions on common edges/faces of cut and uncut cells. The extrapolation of data to shifted boundaries is performed using Taylor expansions. Interface jump conditions are treated similarly sbm . The SBM approach is conceptually simple and backed by theoretical analysis atallah ; sbm1 . A potential drawback is high complexity of closest-point projection algorithms for constructing a map to the true interface.
In the present paper, we introduce a diffuse interface version of SBM, in which surface integration is performed using an algorithm that fits into the framework developed in kublik2016
. Approximating surface integrals by volume integrals, we use level set functions to construct regularized delta functions, approximate global normals, and closest-point mappings. The extrapolation of interface data into the quadrature points for volume integration involves the following steps: (i) find the closest point on the interface; (ii) calculate the quantities of interest at that point; (iii) perform constant extension of normal derivatives (for diffusive fluxes) and linear extension of solution values (for inviscid fluxes). In Step (i), we evaluate the approximate normal at the quadrature point, construct a straight line parallel to the resulting vector, and perform a line search to find the nearest root of the level set function. In Step (ii), normal derivatives are calculated using the given interface value and the interpolated value at the internal edge of the diffuse interface. If tangential derivatives are needed for linear extrapolation in Step (iii), they are approximated using an interpolation stencil centered at the internal point. The way in which the interface data is projected into quadrature points distinguishes our algorithm from shifted boundary methods and extrapolation-based approaches proposed inkublik2016 ; kublik2018 .
The design of our diffuse level set method is motivated by embedded domain problems in which the interface motion is driven by gradients of concentration fields hogea ; chen ; teigen2011 ; teigen2009 ; vermolen2007 . Level set algorithms for such problems require a suitable extension of the interface velocity sethian ; chen ; utz . We perform it in the same way as extrapolation of normal fluxes. The level set function is evolved as in mono ; monons . The support of extended interface data is restricted to a narrow band around the interface. The cost of integrating extrapolated quantities is minimized in this way. For numerical experiments, we design a set of two-dimensional elliptic, parabolic, and hyperbolic test problems with closed-form analytical solutions. The interface is fixed in the elliptic and hyperbolic test cases but is moving in the parabolic scenario. We report and discuss the results of grid convergence studies for our new benchmarks.
2 Sharp interface problem
Many problems of practical interest require numerical solution of conservation laws in evolving domains. Let be a scalar quantity depending on the space location and time . As a representative general model, we consider the parabolic equation
where is an inviscid flux and is a diffusion coefficient. The time-dependent domain is enclosed by an interface and embedded in a fictitious domain with a fixed Lipschitz boundary . The boundary and initial conditions that we impose are given by
The instantaneous position of the interface is determined by a level set function such that is positive in and negative in . This assumption implies that
The vector fields represent extended unit outward normals to . If is a signed distance function satisfying the Eikonal equation , then .
Level set methods osher ; sethian_review evolve using a velocity field defined in . In incompressible two-phase flow models, is obtained by solving the Navier-Stokes equations with piecewise-constant density and viscosity hysing ; monons . In mass and heat transfer models, may represent an extension of a normal velocity depending on concentration gradients hogea ; vermolen2007 . The evolution of is governed by
The exact solution of this linear transport equation satisfies the nonlinear conservation law
is the Heaviside function. However, approximate solutions of (5) may fail to satisfy a discrete form of the volume-of-fluid equation (6). This drawback of the level set approach can be cured using various mass correction procedures. The monolithic level set method that we present in the next section is conservative by construction and does not require any postprocessing (such as redistancing, which is commonly employed to approximately preserve the distance function property).
where is a space of test functions, is the unit outward normal, is the Dirichlet boundary data, and is an approximate normal derivative (to be defined below).
Using the Heaviside function defined by (7), volume integration over the time-dependent embedded domain can be replaced by that over the fixed fictitious domain as follows:
Note that we have added a ghost penalty term depending on a positive function and an extension such that in . Discrete counterparts of such terms are used in unfitted finite element methods for stabilization and regularization purposes burman ; olsh . Introducing the Dirac delta function , formulation (8) can be interpreted as a weak form of (cf. diffuse ; teigen2011 ; teigen2009 )
where the Dirichlet boundary condition (2) is taken into account using . The ghost penalty term imposes the Dirichlet extension condition in the external subdomain and vanishes in , where by definition. The existence of continuous extension operators from a bounded domain to is guaranteed by the theoretical result referred to in olsh . In practice, we extrapolate into to construct (see Section 5.4).
Neumann-type ghost penalties of the form can be defined using a vector field such that in . This version penalizes the weak residual of
If in , a harmonic extension gorb of into is defined by the solution of this boundary value problem. In Section 5.4, we construct a discrete counterpart of using extrapolation. Other kinds of discrete ghost penalties for fictitious domain methods can be found in olsh .
3 Diffuse interface problem
To avoid numerical difficulties associated with volume integration of discontinuous functions and surface integration over sharp evolving interfaces, we approximate (8) by (cf. hysing ; kublik2016 ; zahedi )
where is an extension of the interface flux into . The functions and represent regularized approximations to and , respectively. The design and analysis of such approximations have received significant attention in the literature during the last two decades engquist ; kublik2016 ; kublik2018 ; smereka ; towers ; zahedi . Many diffuse interface methods teigen2011 ; teigen2009 and level set algorithms hysing ; mono ; monons use regularized Heaviside and/or delta functions. The novelty of our unfitted finite element method lies in the construction of , which we discuss in the next section.
where is the Gauss error function.
To evolve in moving interface scenarios, we use the monolithic conservative level set method proposed in mono . The strong form of the underlying nonlinear evolution equation reads
where is a smoothed sign function and is the Eikonal flux. The term depending on a penalty parameter regularizes (12) and enforces the approximate distance function property of . The finite element method for solving (12) uses the mixed weak form mono
The vector that appears in (13) is the unit outward normal to the fixed boundary . If the level set function satisfies the Eikonal equation , then for . A nonvanishing small value of the regularization parameter is used in (14) to avoid ill-posedness at critical points, i.e., in the limit . For a detailed presentation, we refer the interested reader to mono ; monons .
4 Finite element discretization
We discretize (10), (13), and (14) in space using the continuous Galerkin method and a conforming mesh of linear () or multilinear () Lagrange finite elements. The corresponding global basis functions are denoted by . They span a finite-dimensional space and have the property that for each vertex of . The indices of elements that meet at are stored in the integer set . The global indices of nodes belonging to are stored in the integer set . The computational stencil of node is the index set .
The finite element approximations to can be written as
The standard Galerkin discretization of our system yields the semi-discrete problem
where . The diffuse interface thickness is given by . Thus the sharp interface formulation is recovered in the limit . The penalty functions and are chosen to be piecewise constant. Their values in depend on the local mesh size . If necessary, the semi-discrete problem for can be stabilized, e.g., using algebraic flux correction tools afc1 . No stabilization is needed for (16) because attains values in by definition. The coefficients of can be approximated using the explicit lumped-mass formula (mono, , eq. (20))
Further implementation details and parameter settings for (16),(17) can be found in mono ; monons . In Sections 5.1–5.4, we present the algorithms for calculating and . The calculation of an extension velocity for the discretized phase field equation (16) is discussed in Section 5.5.
An appropriate choice of a time integrator for the spatial semi-discretization (15) is problem dependent. In Section 6, we use Heun’s explicit Runge–Kutta method in the hyperbolic test and the implicit Crank–Nicolson scheme in the parabolic one. In the elliptic test, we march the solution of the corresponding parabolic problem to a steady state using a pseudo-time stepping method.
5 Extrapolation using level sets
The main highlight of our diffuse interface method is a simple new algorithm for calculating the fluxes , ghost penalty functions , and extension velocities . We define
using continuous extensions of pointwise solution values, normal derivatives, and normal velocities defined on the zero level set . Of course, the extensions also depend on the mesh size but we omit the subscript for brevity.
5.1 Closest-point search
Let be an approximate signed distance function and an approximation to the extended normal . In our diffuse level set method, and are given by (16) and (17), respectively. However, any other approximation or analytical formula may be used instead. To calculate the extensions and/or for a quadrature point , we first need to find the closest that is connected to by a line parallel to . The point
is usually a good approximation. Indeed, if is an exact SDF, then and, therefore, is the closest point. If is not exact, can be found using a line search along
Note that and for . Let for . Using this parametrization, the exact location of the interface point can be found as follows:
Choose the smallest such that .
Use the bisection method to find the root of .
This way to zoom in on requires repeated evaluation of the level set function at random locations on the straight line . The mesh cells containing the interpolation points are easy to find on uniform meshes. However, the cost of searching becomes large if the mesh is unstructured. Moreover, the bisection method may require many iterations to achieve the desired accuracy.
To locate the point efficiently on general meshes, we use the following search algorithm:
For let with be the next intersection of with the boundary of a mesh cell . Exit the loop when becomes negative for .
Solve a linear/quadratic equation to find the root of .
The sketch in Fig. 1(a) illustrates the construction of for a point on a triangular mesh. For belonging to , it is easy to find and an adjacent cell into which the interface navigation vector is pointing at . By definition, the points and belong to the same mesh cell. Since is linear (for elements) or quadratic (for elements) along the line segment connecting the two points, the root can be easily found using a closed-form formula. The monotone sequence can also be constructed efficiently, especially if the quadrature point lies in a narrow band around .
5.2 Gradient reconstruction
Having determined the location of , we reconstruct the gradient at that point using finite difference approximations on a uniform stencil aligned with the unit outward normal vector
as shown in Fig. 1(a). The straight line through intersects the inner edge of the diffuse interface region at the point . The three-point interpolation stencil
is constructed using a unit vector orthogonal to . In the 3D case, the five-point stencil
is used for interpolation purposes. Note that the cross product is a unit vector orthogonal to and . Hence, the vectors and span the shifted tangential plane at .
If a Neumann (rather than Dirichlet) boundary condition is prescribed on , we use it to define the interface value directly. Otherwise, we use the Dirichlet value and the interpolated value to construct the finite difference approximation
to the normal derivative at . Tangential derivatives are approximated similarly using interpolated values at the points belonging to the interpolation stencil ( in 2D, in 3D). Since the interface thickness parameter is proportional to the mesh size , the reconstructed gradient represents a consistent approximation to with approaching as .
5.3 Flux extrapolation
To evaluate , as defined by (18), at the point , we need the values of and . Since we are using a or finite element approximation, constant extrapolation
is sufficient for the normal derivative. The linear extrapolation formula (cf. sbm1 )
is used for the argument of the inviscid flux. Note that extrapolation is performed along the vector , which is generally not parallel to the normal (see Fig. 1(b)). The gradient to be used in (21) is uniquely defined by the reconstructed normal and tangential derivatives. Its calculation requires a coordinate transformation from the local coordinate system aligned with to the standard Cartesian reference frame. The corresponding transformation formulas can be found, e.g., in vectorlim . No transformations need to be performed for the normal derivatives to be used in (20).
5.4 Ghost penalties
The best way to define a ghost penalty function for (15) depends on the problem at hand and on the available data. A simple and natural Dirichlet extension is given by
Ghost penalties of Neumann type (as mentioned in Remark 2) can be defined, e.g., using
This definition corresponds to constant extrapolation of into . The dimensions of local penalization parameters should be chosen in accordance with the type of the extension (Dirichlet: , Neumann: ). Detailed analysis of ghost penalties is beyond the scope of the present work. However, we envisage that it can be performed following burman ; olsh .
Implicit treatment of ghost penalties changes the sparsity pattern of the finite element matrix. For that reason, we implement implicit schemes / steady-state solvers using fixed-point iterations
based on a split form of the fully discrete problem without penalization (using ). One iteration may suffice if the problem is time dependent and the time step is small.
5.5 Extension velocities
If the evolution of is driven by interfacial phenomena and only the normal velocity is provided by the mathematical model (as in hogea ; vermolen2007 ), a velocity field for evolving the level set function can be constructed using (18) with a suitably defined extension velocity . The value can be determined using constant or linear extrapolation, as in (20) and (21), respectively. For example, suppose that is proportional to a (reconstructed) normal derivative of a concentration field . Then may be used to define . This way to calculate is an efficient alternative to methods that perform a constant extension of
in the normal direction by solving a partial differential equation (cf.sethian ; chen ; utz ).
5.6 Damping functions
If the regularized delta function has a global support, extrapolation needs to be performed at each quadrature point in the fictitious domain. In this case, the cost of numerical integration and matrix assembly can be drastically reduced using damping functions such as
Note that the integrals of and vanish also in subdomains where . Thus an efficient narrow-band implementation is possible.
In a similar vein, the damped extension velocity can be used in (16) to eliminate the surface integral and reduce volume integration to a narrow band. Even if no extension is required because the model provides a vector field defined in , this velocity may be multiplied by to speed up computations and evolve in a domain without inlets (thus eliminating the need to prescribe unknown inflow values in standard implementations of the level set method).
6 Test problems and results
In this section, we apply the diffuse level set method to three new test problems. We design them to have closed-form polynomial exact solutions, which are independent of time even if the equation is time-dependent and the interface is moving. Problems of elliptic, parabolic, and hyperbolic type are defined in Sections 6.1–6.3. The numerical results presented in these sections were calculated using the finite element interpolant of the exact signed distance function . In the elliptic case, the interface is a fixed embedded boundary. The objective of our study is to quantify the error caused by weak imposition and the proposed numerical treatment of interface conditions. In particular, ghost penalties of Dirichlet and Neumann type, as well as the narrow-band damping functions (24), are tested in this experiment. The parabolic and hyperbolic cases are designed to assess additional errors caused by the motion of the interface and by convective fluxes across fixed embedded boundaries, respectively. In Section 6.4, we solve the parabolic problem using (13) and (14) to evolve the level set function numerically. The normal velocity of the interface is extended into using extrapolation with damping (as described in Sections 5.5 and 5.6). This experiment illustrates that all components of our diffuse level set algorithm fit together and produce reasonable results in 2D.
Computations for all test cases are performed on uniform quadrilateral meshes using and a C++
implementation of the methods under investigation in the open-source finite element library MFEMmfem . The numerical solutions are visualized using the open-source software GlVis glvis .
6.1 Elliptic test
The elliptic scenario corresponds to the steady-state limit of (1) with and . We consider the fictitious domain formulation of the boundary value problem
The domain is enclosed by and embedded into . The boundary data and analytical solution are given by
In our numerical experiments, we use the constant penalization parameter for the Dirichlet extension and for the Neumann extension .
The error of the finite element approximation in the embedded domain is approximated by
To keep the integration domain constant on all refinement levels, we use the interface thickness of the finest mesh in this formula.
Table 1 summarizes the results of grid convergence studies for the diffuse level set algorithm using Gaussian regularization (11) and four kinds of ghost penalties. The experimental order of convergence (EOC) is approximately 2 for all types of penalization that we compare in this study (Dirichlet and Neumann, without and with multiplication by the damping functions defined by (24)). Hence, the restriction of level set extrapolation to the narrow band around the embedded boundary has no adverse effect on the accuracy of the proposed method. The numerical solutions obtained using damped ghost penalties and the mesh size are shown in Fig. 2.
|full DGP||EOC||damped DGP||EOC||full NGP||EOC||damped NGP||EOC|
6.2 Parabolic test
The parabolic test problem is a time-dependent version of the elliptic case. We solve (1) with and . The initial-boundary value problem is formulated as follows:
The moving boundary is an expanding circle of radius centered at . The final time is chosen to be short enough for to stay embedded in . As in the elliptic test, the exact solution and the Dirichlet boundary data are given by formula (25).
In this example, and in the following ones, we use Dirichlet ghost penalties (DGP). The penalty parameter is defined as before. Time discretization is performed using the second-order accurate Crank-Nicolson scheme and the time step for the parabolic test. The errors and EOCs are presented in Table 2. Second-order convergence behavior is observed again in this experiment. The damped DGP result corresponding to the mesh size is shown in Fig. 3.
|full DGP||EOC||damped DGP||EOC|
6.3 Hyperbolic test
The third test case is the hyperbolic () limit of (1) with the linear flux , where
is a divergence-free rotating velocity field. We solve the solid body rotation problem
where is the inflow part of and is embedded into . Computations are stopped at the final time . The exact solution
is used to define the boundary data at the inlet when it comes to constructing .
In this example, the numerical flux is given by the extended upwind approximation
where and is calculated using (22) with replaced by
The parameter settings for the space discretization are chosen as before. Time integration is performed using Heun’s method, a second-order explicit Runge-Kutta scheme which is known to be strong stability preserving (SSP) under a CFL-type condition. We use the time step corresponding to the maximum CFL number in our computations. The results are presented in Table 3 and Fig. 4. We observe nearly first-order convergence, except on the finest grid, on which the EOC drops to 0.77. The suboptimal convergence behavior and the occurrence of small oscillations in numerical solutions to the solid body rotation problem are not surprising because we discretize a pure advection equation using the continuous Galerkin method without any stabilization.
|full DGP||EOC||damped DGP||EOC|
6.4 Level set advection
In the numerical studies above, we used the finite element interpolant of the exact signed distance function for extrapolation purposes. In the final test, we use the monolithic conservative level set method mono to compute a numerical approximation for the parabolic test problem formulated in Section 6.2. The constant normal velocity of the interface is approximated by
where is defined in terms of similarly to (19). That is, we use