A new high-order immersed interface method for solving elliptic equations with imbedded interface of discontinuity

Abstract

This paper presents a new high-order immersed interface method for elliptic equations with imbedded interface of discontinuity. Compared with the original second-order immersed interface method of [R.J. LeVeque, Z. Li. The immersed interface method for elliptic equations with discontinuous coefficients and singular sources. SIAM J. Numer. Anal. 31 (1994) 1001–25], the new method achieves arbitrarily high-order accuracy for derivatives at an irregular grid point by imposing only two physical jump conditions together with a wider set of grid stencils. The new interface difference formulas are expressed in a general explicit form so that they can be applied to different multi-dimensional problems without any modification. The new interface algorithms of up to O(h⁴) accuracy have been derived and tested on several one and two-dimensional elliptic equations with imbedded interface. Compared to the standard second-order immersed interface method, the test results show that the new fourth-order immersed interface method leads to a significant improvement in accuracy of the numerical solutions. The proposed method has potential advantages in the application to two-phase flow because of its high-order accuracy and simplicity in applications.

Introduction

Recently, there has been strong interest in developing numerical methods for computing multi-phase flow with unsteady interface. These methods have many practical applications, such as the simulation of the dynamics of gas bubbles in a liquid [11], drop deformation and breakup in viscous flow [44], free surface flow [51], [41], and the breakup of a liquid jet emanating into another fluid [30].

Compared with single-phase numerical methods, algorithms for two-phase flow simulation face additional difficulties related to the interface treatment. Firstly, the shape of the interface can be complex, and can undergo change, merge and breakup during the course of the simulation. Consequently, it is difficult to use body-fitted unsteady grid to fit the evolving interface. A fixed Cartesian grid, where the interface can cut through the grid lines, is often used. In a fixed grid, the interface can be treated by, among others, the volume-of-fluid method, the front tracking method [50], [14], the level set method [42], [43], [35], [36], and the boundary element method [39]. Each of these methods has its own advantages and disadvantages. The volume-of-fluid method is simple and robust. It can maintain a conservation of bubble or droplet volumes. But it is relatively inaccurate in tracking the interface. The front tracking method can track the interface with relatively high accuracy. But it is difficult to use the method to model the connectivity of the interface undergoing complex changes. The level set method, on the other hand, can easily handle the connectivity of complex interface by using a level set function to track the location and movement of the interface.

Secondly, flow variables and their derivatives can be discontinuous across the interface. Specific jump conditions at the interface depend on the physical property of the problems, the unsteadiness of the interface, and the geometric characteristics of the interface. Consequently, special treatments are necessary for computing flow equations at grid points adjacent to the interface (i.e. irregular points). One of the popular methods in treating interface discontinuity is the immersed boundary method (IBM) originally developed by Peskin (see review by Peskin [37]) for simulating blood flow in the heart. The basic idea of the immersed boundary method is to model the interface by adding a delta-function source term to the Navier–Stokes equations. The resulting equations are then discretized by a standard finite difference method in a fixed Cartesian (or non-Cartesian) grid. The singular delta function is regularized by an approximate smooth function spanning a few grid cells. The immersed boundary method has been incorporated in the front tracking method [50] and the level set method [45], [46], [6] in the interface treatment.

The immersed boundary method, however, is only first order accurate in computing two-phase flow with discontinuous solutions across the interface, even though higher-order approximation to the delta function can be achieved for problems with smooth solutions [4], [22], [49], [15], [9]. Beyer and LeVeque [4] studied the approximation to the delta function by a smooth function for the one-dimensional heat equation. The accuracy was measured by a discrete moment condition. They found that it is possible to achieve second-order accuracy by carefully choosing the discrete representation of the delta function. For some problems, however, it is necessary to add a correction term to the approximation to the delta function in order to maintain second-order accuracy. Griffith and Peskin [15] showed that higher-order convergence rates can only be achieved for sufficiently smooth problems. Tornberg and Engquist [49] studied the numerical approximations of singular source terms in differential equations. Specifically, regularization methods for the delta function were analyzed. They showed that fourth-order convergence can be achieved away from the singularity, when a fourth-order difference formula of the ordinary differential operator is coupled to a regularization function with moment order 4. In general, any delta function regularization produces O(h) errors in the neighborhood of the singularity. Consequently, the interface is “smeared” in a numerical solution computed by the immersed boundary method [31].

An alternative to the immersed boundary method is the “sharp-interface” methods which achieve uniformly second (or higher) order accuracy by incorporating the jump conditions into the finite difference formulas. The immersed interface method (IIM) introduced by LeVeque and Li [24] is one of these methods. A similar idea was used earlier by Mayo [34] for the fast solution of the Poisson’s and the biharmonic equations. In presenting their original IIM method, LeVeque and Li [24] considered finite difference methods for the following elliptic equation:

$$\nabla ·(\beta (x)\nabla u(x))+\kappa (x)u(x)=f(x)$$

The equation is defined in a simple region with a uniform Cartesian grid. Fig. 1 shows a schematic of a two-dimensional grid. There is an irregular surface Γ, which may cut across the grid lines, in the computational domain. Across the interface, β, κ, and f may be discontinuous, and along it f may have a delta function singularity. In the derivation of finite difference formulas, the computational grid points are classified into two categories depending on their relative locations with respect to Γ: regular points away from Γ and irregular points adjacent to Γ. A globally O(h²) accuracy is achieved by using the conventional O(h²) central scheme for the regular points and a locally O(h) scheme for the irregular points. In the one-dimensional case, a finite difference formula of O(h) accuracy at an irregular point uses a three-point grid stencil together with an additional correction term. A Taylor series expansion at the interface is used to obtain a set of linear equations for the undetermined coefficients and the correction term. The linear equations are often problem dependent, and they need to be solved numerically every time they are used in the simulation. In order to reach a locally O(h) approximation, the correction term requires jump conditions of up to the second derivatives, i.e.

$$[u]\text{,}[\beta u_x]\text{,}\text{and}[\beta u_{\mathit{xx}}]$$

where [ ] denotes the jump in variables across the interface.

Since the publication of the original immersed interface method [24], there have been many further developments and analysis in various aspects of the immersed interface methods [26], [27], [55], [56], [18], [12], [28], [7], [1], [8], [2]. Among these developments, Wiegmann and Bube [55], [56] developed an explicit-jump immersed interface method for the special cases where the explicit jump conditions of variables and derivatives ([u], [u_x], [u_xx], etc.) are known. Though in a simpler explicit form, the explicit-jump immersed interface method is not applicable to the general jump conditions given by (2).

The immersed interface methods have been applied to the Stokes flow with elastic boundaries or surface tension [25], Hele–Shaw flow [17], incompressible flow based on the Navier–Stokes equations with singular source terms [29], [5], [23], and nonlinear problems in magneto-rheological fluids [19]. Despite these applications, the immersed interface methods are often difficult to apply to complex two or three-dimensional two-phase flow problems. In order to maintain a second-order accuracy, it is necessary to obtain jump conditions at the interface for flow variables and their first and second derivatives. For the Navier–Stokes equations with an interface of discontinuity, it is easy to derive the physical jump conditions for flow variables and their first derivatives across the interface. But it is difficult to obtain jump conditions for the second or higher-order derivatives. In order to develop third or higher-order immersed interface methods, it is necessary to obtain jump conditions for the third and higher derivatives. In addition, the finite difference formulas of the original immersed interface method need to be re-derived for different problems. The coefficients and the correction terms in the finite difference formulas at irregular points cannot be obtained explicitly. They are often computed numerically by solving a matrix equation. The repeated computations for the coefficients and correction terms can be computationally expensive.

Linnick and Fasel [31] presented a high-order immersed interface method for simulating unsteady incompressible flow in an irregular domain. Their method is an extension of the explicit-jump immersed interface method of Wiegmann and Bube [56]. Instead of using analytical jump conditions, they compute the jump conditions for higher derivatives numerically. A fourth-order compact scheme was successfully tested for computing incompressible flow over a cylinder. This method, however, is not applicable to two-phase flow with moving interface with general jump conditions, such as those of Eq. (2). Piraux and Lombard [38], [32] presented another sharp-interface method for numerical computations of interface for wave equations. In order to discretize derivatives at irregular points, a set of modified variables across the interface are computed by using the original variables on both sides of the interface and a set of jump conditions for variables and derivatives. This method also requires the knowledge of jump conditions of high-order derivatives.

Another “sharp-interface” method is the ghost fluid method of Fedkiw et al. [10], [21], [13]. The basic idea is to extrapolate variables on one side of the interface into the “ghost cells” on the other side. Gibou and Fedkiw [13] introduced an O(h⁴) accurate finite difference discretization for the Laplace and heat equations on irregular domain. However, the ghost fluid method is only first order accurate for two-phase flow simulation. Helenbrook et al. [16] presented a second order interface method with ghost cells for incompressible flow with surface discontinuity. The limitation of the method is that it is most suitable for inviscid flow only. It should also be mentioned that sharp interface Cartesian grid methods have been developed for flow with moving solid bodies [54], [53], [52], [3], [20], [33], [40]. These methods can be up to second order accurate. The treatment of irregular points near the solid-fluid interface is mainly based on local polynomial extrapolations.

Based on the brief review above, it is desirable for a high-order immersed interface method to have the following properties:

• Only two physical jump conditions of the variables and their first derivatives should be needed in the second or higher-order immersed interface methods.

• Finite difference formulas at irregular points should be expressed in a general explicit form (without the need to compute matrix equations repeatedly) so that they can be applied to different problems without any modification.

To reach these goals, this paper presents a new set of high-order immersed interface methods. They can be arbitrarily high-order accurate, and require only jump conditions for variables and their first derivatives. The new methods also have the advantage that the finite difference formulas at irregular points are derived in a general explicit form. Compared with the original IIM method, one of the disadvantages of the new methods is that they lead to wider grid stencils. They also do not recover the original finite difference expressions in regular grid points when there is no jump at the interface. In the case of no discontinuity at the interface, however, the current methods are equivalent to a local high-order spline approximation at the interface. While the main purpose of the current approach is to develop higher-order interface treatments, at second order approximation, the current method provides an alternative approach to original second-order immersed interface method.

The derivation, analysis, and test results of the new methods are presented in following sections. Though the motivation of the present work is to apply the methods to multi-phase flow simulation [47], [48], the new high-order immersed interface methods are presented in this paper for elliptic equations in the form of Eq. (1) with imbedded interface of discontinuity only. Nevertheless, the method has potential advantages in the application to two-phase flow because of its high-order accuracy and simplicity in applications by requiring only the physical jump conditions for variables and their first derivatives are needed in the finite difference formulas. The derivation of jump conditions involving the second or higher-order derivatives can be difficult for two-phase flow problem involving the Navier–Stokes equations.

During the review of this paper, one of the reviewers pointed out a paper by Zhou et al. [57] on a new high-order matched interface and boundary (MIB) method for solving elliptic equations with discontinuous coefficients and singular sources on Cartesian grids. The paper was in press by the Journal of Computational Physics at the time of making the revision of this paper. The MIB method is based on the use of fictitious points to achieve high-order accuracy. To construct higher-order interface schemes, the MIB method bypasses the major challenge of implementing high-order jump conditions by repeatedly enforcing the lowest order jump conditions. In treating straight, regular interfaces, MIB schemes up to 16th-order were constructed. For more general elliptic problems with curved, irregular interfaces and boundary, up to 6th-order MIB schemes were demonstrated. The approach presented in the current paper is similar to Zhou et al.’s approach in that both are high-order methods using lower order jump conditions only. The two methods are different in that Zhou’s methods employ fictitious points to achieve high-order accuracy, which is similar to the idea of Ghost Fluid Methods. In addition, the MIB formulas are not explicitly derived. Instead, they are computed by a computer program. The current methods, on the other hand, do not use fictitious points and express the finite difference formulas at irregular points in a generally applicable and explicit form.