The Algebraic Immersed Interface and Boundary Method for Elliptic Equations with Jump Conditions ()
1. Introduction and General Motivations
Simulating flows and heat transfer interacting with complex objects on Cartesian structured grids requires an efficient coupling between such grids and the corresponding numerical methods and complex shape interfaces. Such a coupling is often performed thanks to fictitious domain methods, where the computational domain does not match the physical domain. The advantages of this approach are numerous. A second-order accurate discretization of the spatial operators is simple to obtain; grid generation is trivial, and furthermore there is no need to remesh the discretization grid in the case of moving or deformable boundaries. Concerning this last point, fictitious domain methods can be useful even on unstructured grids: Eulerian fixed unstructured grids can fit immobile obstacles, (e.g. a stator of an aircraft motor) while mobile objects (a rotor) are treated with fictitious domain methods. Two particular classes of problems can be drawn: the immersed boundary problems and the immersed interface problems. The first deals with complex boundaries, such as flow past objects, where no attention has to be paid to the solution inside the obstacles. The immersed interface problems consider subdomains delimited by interfaces, and the solution is required in both sides of the interface. As particular conditions, such as jump conditions, can be required on the interface, this second class of problems is often more difficult to treat.
Let us consider a domain of interest
, typically a fluid domain, which is embedded inside a computational domain of simple shape
, d being the spatial dimension of the problem. The auxiliary domain
, typically a solid particle or an obstacle, is such that as shown in Figure 1:
where
is an immersed interface, and
is the unit outward normal vector to
on
. Let us consider the following model scalar immersed boundary problem in
with a Dirichlet boundary condition (BC) on the interface
:
A boundary condition is also required on the other part of the boundary
so that the whole problem is well-posed.
A first approach dealing with immersed boundaries is the distributed Lagrange Multiplier method proposed by Glowinski et al. [1]. Lagrange multipliers are introduced into the weak formulation of the initial elliptic equation to ensure the immersed boundary condition.
Figure 1. Definition of the subdomains and the interface.
Cartesian grid [2] [3] and Cut-cell [4] methods use a structured grid in the whole domain except near obstacles where unstructured cells are created from structured cells. These methods are hard to implement due to the numerous different space configurations of the intersections between cells and objects. Furthermore, the existence of small cells can induce solver troubles.
The immersed boundary method (IBM) was initially presented by Peskin [5] [6]. Fictitious boundaries are taken into account through a singular source term defined only near the boundaries. As the source term is weighted with a discrete Dirac function smoothed on a non-zero support, the interface influence is spread over some grid cells. This method is first-order in space and explicit. Another class of IBM, the direct-forcing (DF) method, was initially proposed by Mohd-Yusof [7]. The idea here is to impose a no-slip condition directly on the boundary using a mirrored flow over the boundary. In [8] [9], the correct boundary velocity is obtained by interpolating the solution on the boundary and far from the boundary on grid points in the near vicinity of the interface. In [10], Tseng et al. use the same principle but extrapolate the solution in ghost cells outside the domain. This approach can be seen as a generalization of the mirror boundary conditions used in Cartesian staggered grids to impose a velocity Dirichlet condition on pressure nodes. As discussed in [11], this kind of approach seems to be more accurate than [8] [9].
The penalty methods for fictitious domains consist in adding specific terms in the conservation equations to play with the order of magnitude of existing physical contributions so as to obtain at the same time and with the same set of equations two different physical properties. The volume penalty method (VPM) [12] requires the addition of a penalty term
in the conservation equations, such that:
(1)
where
denotes the penalty parameter which tends to 0. Hence, in
the original equation becomes negligible and
is imposed. In [13] [14] authors add a Darcy term
to the Navier-Stokes (NS) equations where
is the dynamic viscosity and K the permeability. In the fluid medium,
so the Darcy term is then negligible and the original set of NS equations is retrieved. In the solid medium,
and consequently the NS equations tend to
. Classical discretizations of the penalty terms are of first order only since they consider the projected shape of the interface on the Eulerian grid to define the penalty parameters [15]. In [16] [17], Sarthou et al. have discretized the volume penalty term with a second order using implicit interpolations as in [10]. This method is called the sub-mesh penalty (SMP) method and has been applied to both elliptic and NS equations.
Applied to problem
, the ghost cell immersed boundary method [10] and SMP method [16] used the cells in
in the vicinity of the interface to enhance the accuracy of the solution in
.
Another approach which considers the extension of the solution is considered in [18] [19] by Gibou and Fedkiw. Ghost nodes and simple interpolations are considered, but contrary to the SMP and the IBM-DF methods, only 1D interpolations are used and the operators are rediscretized “by-hand”.
Let us now consider a model immersed interface problem with jump interface conditions:
A first class of method is the immersed interface methods (IIM) initially introduced by LeVeque and Li [20] and widely described in [21]. This group of methods uses Taylor series expansion of the solution at discretization points in the vicinity of
to modify the discrete operators at these points. Much work has been devoted to the immersed interface method and its numerous applications, such as moving interfaces [22] or Navier-Stokes equations [23]. In [24], Li uses an augmented variables approach. Additional variables and interface equations are added to the initial linear system. The new variables are the values of jumps at some interface points. This method has been extended to the incompressible Stokes [25] and Navier-Stokes [26].
The Ghost Fluid Method, originally developed by Fedkiw et al. [27] [28], introduces ghost nodes where the solution is extended from one side of the interface to the other side. As for IIM, the operator discretization must be modified “by-hand”. Zhou et al. overcome this drawback with the matched interface and boundary (MIB) method [29] [30] [31] by using interface conditions to express the solution at ghost nodes with respect to the solution on physical nodes. Hence, the discretization is automatically performed whatever the discretization scheme. Contrary to [24], the additional equations for these two last methods are not written at “random” points of the interface but at the intersections between the Eulerian grid and the immersed interface. Furthermore, simple Lagrange polynomials are used whereas a more complicated weighted least squares approach is used in [24] to discretize additional equations. In [32], Cisternino and Weynans propose a quite simple method with additional unknowns located at the interface. Interfaces conditions are discretized at these points and are added to the final linear system.
The method presented in this work solves elliptic problems using an augmented method coupled with an auxiliary unknown approach. Its main contribution is its simplicity that eases its implementation, and its generality that allows the method to be applied to various boundary and interface conditions. Contrary to ghost nodes, auxiliary unknowns are present in the linear system and the discretization of the interface conditions is independant from the discretization of the initial equation.
Compact interpolations are used to discretize the additional interface constraints. The method is simple to implement even for interfaces of complex shapes, i.e. not described by analytical equations. Except for the discretization of interface conditions, all operations are automatically performed with algebraic modification or directly by the “black-box” matrix solver. This new method is called the algebraic immersed interface and boundary (AIIB) method. In Section 2, the method is presented for immersed boundary problems. Then, the method is extended to immersed interface problems with known solution on the interface. Finally, the method is applied to immersed interfaces with transmission and jump conditions. A special attention is paid to the management of the discretized interface, especially the way to project it onto the Eulerian grid using a fast ray-casting method. In Section 3, validation tests and convergence studies are presented. Conclusions and perspectives are finally drawn in Section 4.
2. The Algebraic Immersed Interface and Boundary Method
The AIIB method is now presented. The method is first formulated for immersed boundary problems when a Dirichlet or a Neumann boundary condition is required. The method is then extended to simple immersed interface problems where the solution is a priori known on the interface. Finally, an extension to jump and transmission conditions is described.
2.1. Definitions and Notations
Our objective is to numerically impose the adequate conditions on the interface
for problems similar to
and
. These conditions will be discretized in space on an Eulerian structured mesh covering
. As the discretization of the interface or boundary conditions requires interpolations, the following interpolations in 2D:
and
are used. In 3D, we use
and
. An additional interpolation,
, is also possible to be chosen for 2D and 3D problems. The superscript is the dimension of the interpolation while the subscript is the order of spatial accuracy.
The computational domain
is approximated with a curvilinear mesh
composed of
(
in 3D) cell-centered finite volumes (
) for
,
being the set of indexes of the Eulerian orthogonal curvilinear structured mesh. Let
be the vector coordinates of the center of each volume
. In 2D, the horizontal and vertical mesh steps are respectively
and
. This grid is used to discretize the conservation equations. A dual grid is introduced for the management of the AIIB method. The grid lines of this dual cell-vertex mesh are defined by the network of the cell centers
. The volumes of the dual mesh are denoted by (
). The Eulerian unknowns are noted
which are the approximated values of
, i.e. the solution at the cell centers
.
The discrete interface
, hereafter called the Lagrangian mesh, is given by a discretization of the original interface
. It is described by a piecewise linear approximation of
:
,
being the set of indexes of the Lagrangian mesh and K being the cardinal of
. Typically,
are segments in 2D and triangles in 3D. The vertices of each face
are denoted by
for
and the set of all vertices is
,
being the set of indexes of the vertices. The intersection points between the grid lines of the Eulerian dual mesh and the faces
of the Lagrangian mesh are denoted by
(see Figure 2). Our objective is to discretize Dirichlet, Neumann, transmission and jump conditions at these interface points to build a general fictitious domain approach. This method is expected to reach a global second-order spatial accuracy.
We shall use the following Eulerian volume functions in order to implicitly locate
:
• The Heaviside function
, defined as:
(2)
This function is built with a point in solid method presented below. The function
will be used to perform fictitious domain algorithms and to build a level-set function.
• The level-set function
, with:
(3)
and
. The unsigned distance is computed geometrically. The sign is directly obtained with the discrete Heaviside function
.
• The colour phase functions C, which is the ratio of a given phase in a control volume. We denote
the phase ratio in the control volume centered in
. This function is approximated from the
function by using the formula proposed by Sussman and Fatemi [33]:
Figure 2. Definition of the discretization kernels for the AIIB method.
(4)
New sets of Eulerian points
are defined near the interface so that each one has a neighbor
verifying
(with
and
), i.e. the segment
is cut by
. These Eulerian “interface” points are also sorted according to their location inside
or
. Two sets
and
are thus obtained, where
and
.
For each
,
or
, we associate two unknowns: the physical one denoted as
and the auxiliary one
.
2.2. Projection of the Lagrangian Shape on the Eulerian Grid
The generation of the Lagrangian mesh of the interface is achieved using a computer graphics software. Specific algorithms have been developed to project this Lagrangian grid onto the Eulerian physical grid. In order to obtain the discrete Heaviside function
, one has to determine which Eulerian points are inside the domain
defined by a Lagrangian surface. Such a surface must be closed and not self-intersecting. In [11] [14], the authors used a global methodology partly based on [34] where
is obtained thanks to a PDE. This method suffers from a lack of accuracy and robustness. A Ray-casting method based on the Jordan curve theorem is more adapted and is used in the present work. The principle is to cast a ray from each Eulerian point to infinity and to test the number of intersections between the ray and the Lagrangian mesh. If the number of intersections is odd, the Eulerian point is inside the object, and outside otherwise. The Ray-casting method can be enhanced by classifying elements of the Lagrangian mesh with an octree sub-structure which recursively subdivides the space in boxes. If a ray does not intersect a box, it does not intersect the triangles inside the box. A fast and simple optimization is to test if a given point is in the box bounding the Lagrangian mesh. An improvement of the Ray-casting algorithm, the Thread Ray-casting, can be found in [35] [36]. Instead of casting rays in any direction, rays are all aligned on a principal direction of the grid such as each ray passes through a line of points
. Hence, only one ray is cast for each line of points allowing the complexity of the algorithm to be reduced by one order.
The Lagrangian points
of
,
are required to couple the Lagrangian surface and the Eulerian grid used to solve the conservation equations. These points can be obtained with two methods. A geometrical computation of the intersections gives the most accurate result. If not optimized the computational cost of this method is not always negligible for some cases.
Using the Level-set function is a faster but less accurate way to obtain the intersection points. Let us consider two Eulerian points
and
. We denote by
and
the unsigned distances between Eulerian points and the interface
. Then,
.
Algorithmic problems can be encountered if the Lagrangian mesh is too complex compared to the Eulerian mesh. In some cases, more than one segment of the Lagrangian mesh can intersect a same grid segment
, and more than one intersecting points can then be found between
and
with the geometric method. In this case, only one intersecting point is considered.
Concerning the use of the Level-set, this function is a projection of the shape on a discrete grid. The local curvature of the projected shape is thus limited by the accuracy of the Eulerian grid. Consequently, no more than one intersecting point can be found between
and
with the Level-set.
2.3. The AIIB for Dirichlet and Neumann Conditions
2.3.1. General Principle
Once the shape informations are available on the Eulerian grid, the problem discretization has to be modified to take into account the fictitious domain delimited by an immersed boundary or an immersed interface. The sub-mesh penalty (SMP) method [11] [16] was originally designed to treat immersed boundary problems. It could be extended to treat immersed interface problems by symmetrization of the algorithm with introduction of auxiliary unknowns as in the AIIB method. This new method is an enhancement of the SMP method which is also able to solve immersed interface problems. The main idea of the AIIB method is to embed an interface into a given domain by modifying the final matrix only. As no modification of the discretization of the operators is required (contrary to [18] [19] and the immersed interface methods [20]), the AIIB method is thus simple to implement.
Let
be a model problem discretized in the whole domain
as
where A is a square matrix of order m, u the solution vector and b a source term. The basic idea of the AIIB method is to add new unknowns and equations to the initial linear system so as to take into account additional interface constraints. The new unknowns, so-called the auxiliary or fictitious unknowns and labeled with *, are defined as being the extrapolation of the solution from one side of the interface to the other, and are used to discretize the interface conditions. Hence, the original problem
becomes
, with
a square matrix of order
, with n the number of auxiliary constraints related to the interface conditions. The solution
is decomposed such as
and the source term as
. The interface constraints are discretized with a
block matrix C and the source term
. The discretizations of the constraints for various interface and boundary conditions are given in Sections 2.3.2 and 2.4.
According to the interface conditions, the regularity of the solution on the interface is often lower than in the rest of the domain. Hence, the discretization of operators with a stencil cutting the interface can induce a great loss of accuracy. The first idea is to consider unknowns
(resp.
) as the extension of the solution in
(resp.
). The initial algebraic link between unknowns from both sides of the interface is cut, and the new link over the interface is obtained thanks to auxiliary unknowns. Practically, matrix coefficients must be modified to take into account the new connectivities. Let
be a coefficient of A at row I, column J and
the new coefficient in
. If
and
,
.
Hence, there is no more direct link between the initial unknowns in
and
. However, the initial value of
is reused to express the new relation between the original unknown
and its new neighbor
. Consequently,
, where
is the index corresponding to
.
This is exactly the way how we proceed for the practical algorithm. However, this modification can be expressed algebraically with permutation and mask matrices as follows.
We define the two following mask matrices
of dimensions
and
of size
:
(5)
(6)
The matrices
and
are defined such as
,
if
, else
. Similarly
if
else
. Finally, the connectivities are changed using the permutation matrices
and
:
is defined to switch row I with row J if
,
and
to switch row I with row J if
,
. Hence, the new problem matrix is now defined by:
(7)
The new problem is
with
written with 4 blocks of various sizes:
,
,
,
. The matrix
is thus the modification of the initial matrix A by setting to zero the coefficient
if
, and
and
are the two sub-matrices of the matrix C. The problem can be written as:
(8)
The entire problem can then be solved to obtain
. However,
being the auxiliary solution is not required to be computed explicitly. Hence, the Schur complement method can be used to calculate the solution for the physical unknowns only. The final problem is now:
(9)
The opportunity of such a reduction will be discussed later.
2.3.2. AIIB Algorithm for a Scalar Equation with Dirichlet Boundary Conditions
For the sake of clarity, let us first describe in 2D the AIIB method for the model scalar problem
with a Dirichlet boundary condition on the interface
. For this version of the AIIB algorithm,
is the domain of interest and auxiliary unknowns are created in
only. Let us consider a point
. At location
, two unknowns coexist: a physical one
and an auxiliary one
. We first describe the case when
has only one neighbor
in
. The Lagrangian point
is the intersection between
and
(Figure 3 right). Then, the solution
at the interface is approximated by the
interpolation between the Eulerian unknowns
and
:
(10)
To our knowledge, it has been verified that for schemes close to the one presented here (including ghost-like schemes [19] and IBM-Direct Forcing schemes [10]) only a linear interpolation is required to reach a second order of accuracy as a Dirichlet BC is imposed here. It has been theoretically demonstrated in [37] that under the assumption of pointwise stability of the scheme the discretization at the boundary can be two orders less accurate than in the rest of the domain without reducing the global convergence order. It is important to notice this does not apply for immersed interface conditions [32].
If now
has a second neighbor
in
, the intersection
between
and
is considered with
. We choose
, a new point of
between
and
(see Figure 3 left). The solution
is then imposed using a
-interpolation of the values
,
and
:
(11)
A
interpolation of
,
,
and
can be also used by extending the interpolation stencil with the point
which is the fourth point of the cell of the dual mesh defined by
,
and
(see Figure 3 left). As a third
(a) (b)
Figure 3. Example of selection of points for Dirichlet (a) and Neumann (b) constraints.
choice, two independent
linear 1D interpolations can be used (one for each direction) for an almost equivalent result. It produces:
(12)
In this case, two auxiliary unknowns are created.
A simple choice for
is the barycenter between
and
where
. A summation of the two constraints from (12) gives:
(13)
which is equivalent to build a constraint imposing
at
with a
interpolation:
(14)
Hence, an easy general implementation consists in summing the constraints corresponding to each direction, no matter the number of neighbors of
. If the elements
of
used to define
and
are not the same, the barycenter
of these two points is not necessarily on
, especially for interfaces of strong curvature. However, the distance
between
and
varies like
and so this additional error does not spoil the second-order precision of our discretization. The convergence of this additional error is numerically tested in Section (3.3.1). If the curvature of
is small enough relatively to the Eulerian mesh, i.e. if the Eulerian mesh is sufficiently fine,
almost never has a third or a fourth neighbor in
. However, if this case appears, a simple constraint
is used with
being an average of
at the neighbor intersection points. In any case, by decreasing the Eulerian mesh step h, the number of points
having more than two neighbors in
also decreases.
Concerning the C submatrix, it is filled with the coefficients
. One line corresponds to the discretization of the boundary constraint related to one auxiliary node. Each value of the second member
is the boundary value
related on auxiliary point.
Hence, the present method is suitable to impose a Dirichlet boundary condition on
for
, when the solution in
is of no interest. The solution
for
is an extrapolation of the solution in
in order to satisfy the boundary condition on
and thus is non-physical. Hence, the solution at the nodes of
far from the interface does not impact on the solution in
as the unknowns
for
have no more algebraic relations with the nodes
for
. Nevertheless, if nothing specific is performed, a solution with no interest is computed by default in
. A prescribed solution inside
can be obtained easily with a volume penalty method such as VPM [13] where initial equations are penalized to impose a desired analytical solution. However, a desired solution in
can also be obtained by modifying the solution obtained by solver. Then, the computational cost of the approach can be reduced by switching the solving of
off, or by totally removing these nodes in the solving matrix.
2.3.3. Symmetric Version for Dirichlet Interface Conditions
The next step is to allow for multiple Dirichlet boundary conditions on both sides of the immersed interface. Thin objects could be treated with this approach. The problem is now:
(15)
The problem (15) requires for each point
a physical unknown
as well as an auxiliary unknown
on both sides of the interface.
Practically, the algebraic modification part of the AIIB algorithm for a Dirichlet BC is applied a first time with
as domain of interest, and auxiliary unknowns are created near
in
. As a second step, the Heaviside function is modified as
and the algorithm is applied a second time. Now,
is the domain of interest and auxiliary unknowns are created near
in
.
Then, both problems are solved at the same time with the resolution of the linear system.
2.3.4. AIIB Algorithm for a Scalar Equation with Neumann Boundary Conditions
Let us now consider the following model scalar problem with a Neumann BC on the interface
:
(16)
The principle is about the same as for Dirichlet BC, and the same interpolations, once derived, can be used to approximate the quantity
. Hence, at any point
on
we use
(17)
For
, we get
whereas for
,
is obtained which means that the normal gradient is approximated by a constant over the whole support.
For example, in the configuration of Figure 3(a), with
, we have:
(18)
The diagonal coefficient of the row related to
in
is
. The case where
leads to numerical instabilities. However, it can happens for this configuration case only if
and
have the same sign, which will not happen if the interface is sufficiently smooth. For a smooth enough interface,
and
. A smoothed interface can be simply obtained using the normal vector of the segment
implies that the signs of
and
are always different so the diagonal coefficient is always dominant. The same property occurs for the other stencil configurations.
When
has only one neighbor
in
, the
and
interpolations degenerate to
interpolations which suit for Dirichlet BC. For Neumann BC, this loss of dimension no longer allows the interface orientation to be accurately taken into account, as one of the components of the normal unit vector disappears from the interfacial constraint. Hence, a third point
in
is caught to build
interpolations (see Figure 3 right). This point is a neighbor of
and is taken as
. In 2D, two choices generally appear, and the point being so that the angle
is in
is taken.
2.3.5. Algebraic Elimination Using the Schur Complement
The Schur complement method allows an algebraic reduction to be performed. For a Dirichlet or Neumann BC, each constraint is written in such a way that only one auxiliary unknown is needed:
(19)
where
is the source term. In this case, the matrix
in (8) is diagonal and thus the Schur complement
is easy to calculate. Practically, when the algebraic reduction is made,
is built directly by the suitable modification of A without considering the extended matrix
. The part
is then added to
whereas
is added to b. As it will be subsequently demonstrated, the algebraic reduction decreases the computational cost of the solver by 10% - 20%.
If only
interpolations are used with the algebraic elimination, the matrix obtained with this method is similar to the one obtained in [18] for a Dirichlet problem. However, in this last paper, the auxiliary unknowns are taken into account before the discretization of the operator which requires additional calculations for each discretization scheme.
If
interpolations are used, the computed solution in
are very similar (within machine error if the solver perfectly converges) to those that can be obtained with the SMP [16] method (when the penalty parameter is sufficiently small) and the DF-IB method [10]. These methods discretize equivalent boundary constraints, but the way the constraint is applied to the initial discretization of the equation for the nodes
differs. The SMP method uses a penalty term and the DF-IB method uses terms of opposite signs to erase some part of the initial equation. The discretization matrix obtained with both methods is not equivalent to the one obtained with the AIIB method, with or without algebraic reduction. With algebraic reduction, the discretization for the nodes
is modified, and without algebraic reduction, both auxiliary and physical unknowns coexist at
. The accuracy of these methods will be discussed in the next section.
The present algorithm seems simpler, as the standard discretization of the operators is automatically modified in an algebraic manner. So, various discretization schemes of the spatial operators can be used. However, the discretization of an operator at
can only use in
the fictitious unknowns and not the physical ones. Hence, the only limitation concerns the stencil of these operators which have to be limited, if centered, to three points by direction.
2.4. The AIIB for Immersed Interface Problems with Jump Conditions
With the symmetric method described in (2.3.3), the problem can be solved on both sides of the interface when explicit Dirichlet BC are imposed. For many problems, the solution is not a priori known on the interface and some jump transmission conditions on the interface
are required. Let us now consider the problem:
where the interface conditions are:
(20)
(21)
The notation
denotes the jump of a quantity over the interface
. In the symmetric version of the AIIB method, a given intersection point
, is associated with two auxiliary unknowns, one on each side of the interface. Hence, the interface constraints (20) and (21) of
can be imposed at each intersection point
by using the two auxiliary unknowns. For example, the row I of the matrix
with
, can be used to impose the constraint (20) and the line J of the matrix with
, is then used to impose the constraint (21).
2.4.1. The Solution Constraint
The symmetrized AIIB methods for Dirichlet BC reads:
(22)
when
interpolations are used. With
, we obtain:
(23)
which is the first constraint to be imposed.
2.4.2. The Flux Constraint
Following the same idea and using
interpolations,
(24)
for the case presented in Figure 3 left. The normal is pointing from the + to the − sides. Using (21), we obtain:
(25)
which is the second constraint to be imposed. With such an interpolation, the solution gradient is constant over the whole stencil. As demonstrated later, the second-order accuracy can be reached on Cartesian grids when
.
Three auxiliary unknowns are thus involved in the discretizations (23) and (25). The auxiliary unknown
is also involved in the discretization of (20) and (21) at another intersection point on
. Hence, the whole system
is closed.
One can notice that in [32], the same kind of augmented system is considered. Contrary to the AIIB method, no auxiliary nodes are used and the spatial discretization at the grid points at the vicinity of the interface has to be modified “by-hand”.
2.4.3. Algebraic Reduction
Since we need more than one auxiliary unknown to discretize each constraint, the matrix
is not diagonal and a solver has to be used to compute
. This matrix is not diagonally dominant. For instance, if the discretized solution constraint (23) is inserted in
for the row related to
, column
has a coefficient
while the coefficient of the column
has a coefficient
. In many situations,
is greater that
, and the situation can be critical if
. The equivalent behaviour can happen with the discretized flux constraint (25). This point is discussed later in Section 3.4.
For the matched interface and boundary (MIB) method, Zhou et al. [29] use a different discretization of the interface conditions which allows an easy algebraic reduction which is directly performed row by row.
The algebraic reduction for the immersed interface problems has not been yet implemented. However, the standard discretization of the AIIB method requires a more compact stencil than for the MIB method and then the need for a reduction seems to be less critical. More generally, as auxiliary points are in the vicinity of a surface immersed in a volume, their count is one order below the total count of grid points. It suggests that the number of unknowns added to the initial matrix is low compared to the total number of unknowns.
3. Numerical Results for Scalar Problems
Elliptic equations are discretized using the standard second-order centered Laplacian. For all problems, similar results have been obtained with a PARDISO direct solver [38], and an iterative BiCGSTAB solver [39], preconditioned under a ILUK method [40]. Unless otherwise mentioned, a numerical domain
is used for every simulation. Two discrete errors are used.
The discrete relative
error in a domain
is defined as:
(26)
(27)
with
the analytical solution.
The discrete
error is defined as:
(28)
without AIIB method, the present code only reach a first-order spatial error convergence rate for cases with irregular interfaces. With the AIIB method, the code is expected to reach a second-order convergence rate as linear interpolations are used to approximate the solution at the vicinity of the interface.
Only
is taken into account for the immersed boundary problems. For all cases,
is the boundary of
as illustrated in Figure 1. If not explicitly written,
and
have a default value of 0. For all cases, the analytical solution is imposed for the cells at the boundary of the numerical domain.
3.1. Immersed Boundary Problems
3.1.1. Problem 1
The homogenous 2D Laplace equation is solved. The interface
is a centered circle of radius
with a Dirichlet condition of
. An analytical solution accounts for the presence of a second circle with a radius
and
is imposed on the boundary conditions. The analytical solution is:
(29)
Accuracy tests are performed with
,
and
interpolations. Figure 4 shows the solution and the error map for a 32 × 32 mesh with
interpolations. The same results are always obtained with and without algebraic reduction. Figure 5 shows the convergence of the error for the
and
norms. For all interpolations, the convergence slopes are approximately 2 for the relative
error. For the
error, the slopes are about 1.8. The
interpolation is the more accurate, followed by the
interpolation although it uses more auxiliary points (but a smaller stencil). However, the differences of accuracy between the different interpolations remain small. The performances of the ILUK-BiCG-Stab
Figure 4. Solution and error map for problem 1 on a 32 × 32 grid.
Figure 5. Curves of errors for Section 3.1.1.
solver are now benchmarked for the three interpolations with and without algebraic reduction and for the SMP method. Table 1 shows the computational times of the matrix inversions (average time in seconds for 25 matrix inversions) and Table 2 shows the time ratio between the standard and the reduced matrix. Except for the
interpolation on the 1024 × 1024 mesh, the differences between the two methods seem to decrease with the size of the matrix. In fact, as interfaces are
manifolds, the number of intersection points does not increase as fast as the Eulerian points. Hence, the ratio between the size of a reduced and a complete matrix tends to 1. The computational time for the SMP method is quite similar to the one obtained with the AIIB method and algebraic reduction. Figure 6 shows the convergence of the ILUK-BiCG-Stab solver for
Table 1. Computational times in seconds for problem 1. Tests are performed with three different interpolations with (red) and without (std) algebraic reduction, and compared to the SMP method.
Figure 6. Residual against iterations of ILUK solver for problem 1 with a 256 × 256 mesh.
the seven interpolations with a 256 × 256 mesh. The same behavior has been observed for other grid resolutions. The type of interpolation does not significantly impact on solver performances. As expected, the number of solver iterations has to be increased to reach a given residual when the number of computational nodes is also increased.
As a conclusion on the interpolation type,
interpolation seems to be the most desirable as it is generally slightly easier to implement, and requires less memory than a
interpolation. To explain that, one can suppose that the boundary constraint forces the points of the stencil to follow a polynomial solution. As the real solution is not polynomial, adding points to the stencil increases the zone where the analytical solution is more difficult to match. For the next problems,
will be used for constraints on the solution.
3.1.2. Problem 2
The 3D equation
is solved. The solution is
. The solution is imposed on an immersed centered sphere of radius 0.2. Results of the numerical accuracy test with the spherical inner boundary are presented in Figure 7. The average slope for the
norm is 2.33 and increases for the denser meshes. Even if the method has a convergence order similar to the order of the polynomial solution, computer error accuracy cannot be expected because the immersed boundary
is an approximation of the initial interface
. The analytical solution
Table 2. Ratio of computational times for reduced and standard matrices for Section 3.1.1.
Figure 7. Curves of errors for problem 2.
(the value at the boundary) is imposed on the boundary, but the boundary location is approximated so an error which has the order of the discretization error of the interface is introduced.
3.1.3. Problem 3
The 3D equation
is solved. The solution is
. Results with and without an immersed interface (a circle of radius
) with the analytical solution imposed on it. The results are presented in Figure 8. For the
norm, the second order is regularly obtained. For the
norm, the convergence order tends to a second order with the size of the mesh. As can be noticed by comparing results with and without the AIIB method, this last one does not spoil the convergence order of the code, and the presence of the immersed interface with an analytical solution imposed on
improves the accuracy. For both cases, the numerical solution tends to a second order in space.
3.1.4. Problem 4
The 2D equation
is solved. The analytical solution is imposed on the boundaries of the domain and a Neuman BC is imposed on a centered circle of radius
. As can be seen in Figure 9, the global convergence has an average slope of 1.10 and can be explained by the simplicity of the discretization of the Neumann BC.
Figure 8. Curves of errors for problem 3.
Figure 9. Curves of errors for problem 4.
3.2. Immersed Interface Problems
3.2.1. Problem 5
The 2D problem
with
and
is solved. As the equation remains the same in both domains, this problem can be solved without immersed interface method. The analytical solution is
. As can be expected with our second-order code, computer error is reached for all meshes with or without AIIB method. The difference with problem 2, where the solution is a second-order polynomial too, is that instead of imposing the solution on an approximated shape, we impose a solution jump. If the numerical nodes stay in the correct side of the interface (where they should be if the interface were not approximated), the location of the jump does not influence the numerical solution.
In Figure 10 the same result is obtained with a solution jump
on a circular interface such as
for
(
) and
otherwise (
).
An equivalent quality of result is obtained (see Figure 11) with
such as:
Figure 10. The solution and the error for problem 5 with a 32 × 32 mesh.
Figure 11. The solution and the error for problem 5 with a 64 × 64 mesh.
(30)
with
. The small stencil of the method allows interfaces with relatively strong curvatures to be used.
3.2.2. Problem 6
For this case,
is now considered with a discontinuous coefficient a such as
in
and
in
, involving the following analytical solution:
(31)
Accuracy tests are first performed with the interface almost passing by some grid points (the mesh is referred as centered). The interface does not strictly lies on these points, as the shape is shifted by a value of 10−10. This configuration is difficult as the interpolations degenerates. Accuracy tests are then performed with a numerical domain of size 1.0001 (called shifted mesh as the center of
does not match the center of the numerical domain). In this configuration, the interface never passes by a grid point. The numerical errors on the solution are presented in Figure 12. For the centered mesh, the slope is 1.86 for the
and
errors. For the even series, where no geometrical singularity is present, the slope for both errors is 2.04.
Figure 13 shows the solution and the
relative error for a 32 × 32 mesh. As the analytical solution is imposed on the numerical boundary, the error is principally located in the interior subdomain.
3.2.3. Problem 7
The homogenous 2D Laplace equation is considered with the following analytical solution:
(32)
Figure 12. Curves of errors for the centered and shifted meshes for problem 6.
Figure 13. The solution and the error for problem 6 with a 33 × 33 mesh.
where
and
are delimited by
a centered circle of radius 0.5. The jumps coefficients are
,
and
.
Figure 14 shows that the convergence for both
and
error are of first order only. Figure 15 shows the numerical solution (which is not so different from the analytical solution) and the error map for a 32 × 32 mesh.
Hence, the drawback of the compacity of the discretization is a first-order convergence in space only for cases with flux jump. As for the Neumann BC, the simplicity of the discretization of the flux condition seems to be the cause of this low order. Contrary to the discretization of the solution jump, the precise position of the interface is not present which induces a first order of convergence. Generally, other approaches [20] [29] [32] use the precise position of the interface (at the expanse of a larger stencil) and can reach higher orders with all flux conditions.
3.3. Shape Management
3.3.1. Convergence
We measure the sensibility of the method with the accuracy of the discretization of the immersed interface. Problem 1 is solved on 32 × 32 and 128 × 128 meshes. Figure 16 shows the accuracy of the solution with respect to the number of points used to discretize the interface which is here a circle. The reference solutions (Figure 5) have been computed with an analytical circle. As can be seen, a second order in space is globally obtained. The numerical solutions of reference for the 32 × 32 and 128 × 128 meshes are different but the sensitivity of the error to the number of points in the Lagrangian mesh is almost the same.
3.3.2. The Stanford Bunny
This last case demonstrates how a method with a second-order convergence rate
Figure 14. Convergence of the
relative error and the
error for problem 7.
Figure 15. The solution and the
relative error for problem 7 with a 32 × 32 mesh.
on spatial accuracy enhances the representation of the boundary condition compared to a first-order method such as the VPM which imposes a prescribed solution at the center of the control volume as shown in Equation (1).
The homogenous Laplace problem with a Dirichlet BC
is solved on a
mesh bounding an obstacle of complex shape (the Stanford bunny). The extension of the solution in
is used for the post treatment. Thus, all
, are replaced by
. Then, the iso-surface
gives an idea of the approximation of the boundary condition. Figure 17 shows the iso-surface for a first order method. As can be seen, the shape of the obstacle endures a
Figure 16. Convergence of the error with respect to the number of elements forming the Lagrangian shape for problem 1.
Figure 17. Iso-surface
for the Stanford bunny with a first order method.
rasterization effect as the solution is imposed in the entire control volumes. Figure 18 shows the iso-surface for the second order AIIB method. The improvement brought by this approach is straightforward. The graphic processing has been applied to both methods. Figure 19 shows a slice of the solution passing through the bunny. As can be seen, overshoots are present inside the shape which corresponds to the auxiliary values allowing the correct solution at the Lagrangian interface points to be obtained.
3.4. Some Remarks about the Solvers
The kind of interpolation function used and the position of the interface have an impact on the final discretization matrix
, especially on its conditioning. Let us consider an intersection
of
between two points
and
. A Dirichlet BC
is imposed on it. The constraint constructed with a
interpolation is
, with
. Hence,
tends to 0 when
tends to
. As the matrix loses its diagonal dominance, solver
Figure 18. Iso-surface
for the Stanford bunny with a second order method.
Figure 19. Iso-surface
and a slice of the solution.
problems can be encountered. Tseng et al. [10] proposed changing the interpolation by using a new node which is the image of
through the interface. In [18] [19], authors pointed out this problem and suggest to slightly move the interface to a neighboring point (in our case
) if
is too close to
.
In our approach, a simple solution is to shift the interface such as the interface matches
. The loss of accuracy depends on the threshold of distance between
and
under which the solver is unstable. In our case, test 3.2.2 suggests that a threshold below the error of the discretization of
into the piecewise linear interface
can be taken.
In the case of interface shift, for the Dirichlet BC, an unknown
is created, and the equation at
is simply
. For the Neumann BC, the standard interpolation is written at
with
and its neighbor unknowns in
.
For the transmission conditions (20)-(21), if
and
, no auxiliary unknown is created and the standard finite-volume centered discretization is used. However, for this case, or for
and
, our implementation using ILUK preconditioner or a PARDISO direct solver does not necessarily require such methods, even if
. Correct solutions are always obtained.
4. Discussion and Conclusions
A new immersed interface method, the algebraic immersed interface and boundary (AIIB) method, which uses algebraic manipulations and compact stencil discretizations, has been presented. This method is able to treat elliptic equations with discontinuous coefficients and solution jumps over complex interfaces.
For the immersed boundary problems with a Dirichlet BC, the method has shown a second order of convergence in space for various kinds of interpolations. An algebraic reduction has been applied to accelerate the convergence of the solver. For the Neumann BC, a first order has been obtained. For the immersed interfaces, a second order of convergence in space is obtained when the jump of the normal flux is zero, even if the equation has discontinuous coefficients. For a non-zero flux jump, the method reaches a first order.
The main advantage of the method is its simplicity and its ease of implementation. First, the method does not require to modify the original discretization of the equation and works only at the level of the discretization matrix. Second, the method has a compact stencil and interfaces with high curvature can be treated with no particular attention.
Future work could be devoted to increase the accuracy of the method when the jump of the normal flux is not zero. It is the main drawback of the method compared to the IIM, MIB method or [32]. A challenge will be to keep a compact stencil. A study of the matrix conditioning would be important too as a strong solver is required for the present method. Another interesting point would be to couple the method with alternative interface conditions such as the Jump Embedded Boundary Conditions proposed in [41] which are:
(33)
(34)
where
denotes the arithmetic mean of the traces of a quantity on both sides of the interface;
,
, h and g are scalar values which can be chosen in order to obtain Dirichlet, Neumann, Fourier transmission conditions on the interface.
As the method acts on an algebraic level, it can be extended to other implicit approaches where the solution is obtained through the resolution of a linear system. It would be interesting to apply this method to finite element approaches. The interface conditions currently written for finite volumes could be discretized in a finite element manner.
At last, a long-term goal is to extend the method to the Navier-Stokes equations with immersed interfaces. To perform fluid/structure coupling, the implicit tensorial penalty method [42] [43] can be considered. With this approach, the solid medium is treated as a fluid with a high viscosity. At the fluid/solid interface, average physical quantities are imposed. Such strategy is generally less accurate than methods using polynomial interpolations so a coupling with the AIIB method is desirable.
Acknowledgements
The authors greatly thank the reviewers for their accurate and pertinent questions and remarks that helped to improve this article. The authors also thank Lisl Weynans for the helpful discussions.
Definition of the Parameters