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. . 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   and Cut-cell  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  . 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 . The idea here is to impose a no-slip condition directly on the boundary using a mirrored flow over the boundary. In  , 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 , 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 , this kind of approach seems to be more accurate than  .
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)  requires the addition of a penalty term in the conservation equations, such that:
where denotes the penalty parameter which tends to 0. Hence, in the original equation becomes negligible and is imposed. In   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 . In  , Sarthou et al. have discretized the volume penalty term with a second order using implicit interpolations as in . 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  and SMP method  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   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  and widely described in . 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  or Navier-Stokes equations . In , 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  and Navier-Stokes .
The Ghost Fluid Method, originally developed by Fedkiw et al.  , 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    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 , 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  to discretize additional equations. In , 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:
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:
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 :
Figure 2. Definition of the discretization kernels for the AIIB method.
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  , the authors used a global methodology partly based on  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  . 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   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   and the immersed interface methods ), 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 :
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:
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:
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:
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 :
To our knowledge, it has been verified that for schemes close to the one presented here (including ghost-like schemes  and IBM-Direct Forcing schemes ) 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  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 .
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 :
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
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:
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:
which is equivalent to build a constraint imposing at with a interpolation:
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  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:
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 :
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
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:
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:
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  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  method (when the penalty parameter is sufficiently small) and the DF-IB method . 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:
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:
when interpolations are used. With , we obtain:
which is the first constraint to be imposed.
2.4.2. The Flux Constraint
Following the same idea and using interpolations,
for the case presented in Figure 3 left. The normal is pointing from the + to the − sides. Using (21), we obtain:
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 , 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.  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 , and an iterative BiCGSTAB solver , preconditioned under a ILUK method . 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:
with the analytical solution.
The discrete error is defined as:
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:
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.
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:
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:
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    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
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.  proposed changing the interpolation by using a new node which is the image of through the interface. In  , 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 . 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  which are:
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   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.
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
 Glowinski, R., Pan, T.W., Hesla, T.I. and Joseph, D.D. (1999) A Distributed Lagrange Multiplier/Fictitious Domain Method for Particulate Flows. International Journal of Multiphase Flow, 25, 755-794.
 Johansen, H. and Colella, P. (1998) A Cartesian Grid Embedded Boundary Method for Poisson’s Equation on Irregular Domains. Journal of Computational Physics, 147, 60-85.
 McCorquodale, P., Colella, P. and Johansen, H. (2001) A Cartesian Grid Embedded Boundary Method for the Heat Equation on Irregular Domains. Journal of Computational Physics, 173, 620-635.
 Ye, T., Mittal, R., Udaykumar, H.S. and Shyy, W. (1999) An Accurate Cartesian Grid Method for Viscous Incompressible Flows with Complex Immersed Boundaries. Journal of Computational Physics, 156, 209-240.
 Fadlun, E.A., Verzicco, R., Orlandi, P. and Mohd-Yusof, J. (2000) Combined Immersed-Boundary Finite-Difference Methods for Three-Dimensional Complex Flow Simulations. Journal of Computational Physics, 161, 35-60.
 Tseng, Y.-H. and Ferziger, J.H. (2003) A Ghost-Cell Immersed Boundary Method for Flow in Complex Geometry. Journal of Computational Physics, 192, 593-623.
 Sarthou, A., Vincent, S., Caltagirone, J.-P. and Angot, P. (2008) Eulerian-Lagrangian Grid Coupling and Penalty Methods for the Simulation of Multiphase Flows Interacting with Complex Objects. International Journal for Numerical Methods in Fluids, 56, 1093-1099.
 Angot, P., Bruneau, C.-H. and Fabrie, P. (1999) A Penalization Method to Take into Account Obstacles in Incompressible Viscous Flows. Numerische Mathematik, 81, 497-520. https://doi.org/10.1007/s002110050401
 Khadra, K., Angot, P., Parneix, S. and Caltagirone, J.-P. (2000) Fictitious Domain Approach for Numerical Modelling of Navier-Stokes Equations. International Journal for Numerical Methods in Fluids, 34, 651-684.
 Lacanette, D., Vincent, S., Sarthou, A., Malaurent, P. and Caltagirone, J.P. (2009) An Eulerian/Lagrangian Method for the Numerical Simulation of Incompressible Convection Flows Interacting with Complex Obstacles: Application to the Natural Convection in the Lascaux Cave. International Journal of Heat and Mass Transfer, 52, 2528-2542.
 Ramière, I., Angot, P. and Belliard, M. (2007) A General Fictitious Domain Method with Immersed Jumps and Multilevel Nested Structured Meshes. Journal of Computational Physics, 225, 1347-1387.
 Gibou, F., Fedkiw, R.P., Cheng, L.-T. and Kang, M. (2002) A Second-Order-Accurate Symmetric Discretization of the Poisson Equation on Irregular Domains. Journal of Computational Physics, 176, 205-227.
 Gibou, F. and Fedkiw, R. (2005) A Fourth Order Accurate Discretization for the Laplace and Heat Equations on Arbitrary Domains, with Applications to the Stefan Problem. Journal of Computational Physics, 202, 577-601.
 Leveque, R.J. and Li, Z. (1994) The Immersed Interface Method for Elliptic Equations with Discontinuous Coefficients and Singular Sources. SIAM Journal of Numerical Analysis, 31, 1001-1025.
 Hou, T.Y., Li, Z.L., Osher, S. and Zhao, H. (1997) A Hybrid Method for Moving Interface Problems with Application to the Hele-Shaw Flow. Journal of Computational Physics, 134, 236-252.
 Lai, M.-C. and Tseng, H-.C. (2008) A Simple Implementation of the Immersed Interface Methods for Stokes Flows with Singular Forces. Computers and Fluids, 37, 99-106.
 Li, Z.L., Ito, K. and Lai, M.-C. (2007) An Augmented Approach for Stokes Equations with a Discontinuous Viscosity and Singular Forces. Computers & Fluids, 36, 622-635.
 Ito, K., Lai, M.-C. and Li, Z.L. (2009) A Well-Conditioned Augmented System for Solving Navier-Stokes Equations in Irregular Domains. Journal of Computational Physics, 228, 2616-2628.
 Fedkiw, R.P., Aslam, T., Merriman, B. and Osher, S. (1999) A Non-Oscillatory Eulerian Approach to Interfaces in Multimaterial Flows (The Ghost Fluid Method). Journal of Computational Physics, 152, 457-492.
 Liu, X.-D., Fedkiw, R.P. and Kang, M. (2000) A Boundary Condition Capturing Method for Poisson’s Equation on Irregular Domains. Journal of Computational Physics, 160, 151-178.
 Zhou, Y.C., Zhao, S., Feig, M. and Wei, G.W. (2006) High Order Matched Interface and Boundary Method for Elliptic Equations with Discontinuous Coefficients and Singular Sources. Journal of Computational Physics, 213, 1-30.
 Zhou, Y.C. and Wei, G.W. (2006) On the Fictitious-Domain and Interpolation Formulations of the Matched Interface and Boundary (MIB) Method. Journal of Computational Physics, 219, 228-246.
 Yu, S.N., Zhou, Y.C. and Wei, G.W. (2007) Matched Interface and Boundary (MIB) Method for Elliptic Problems with Sharp-Edged Interfaces. Journal of Computational Physics, 224, 729-756.
 Cisternino, M. and Weynans, L. (2012) A Parallel Second Order Cartesian Method for Elliptic Interface Problems. Communications in Computational Physics, 12, 1562-1587.
 Sussman, M. and Fatemi, E. (1999) An Efficient, Interface-Preserving Level Set Redistancing Algorithm and Its Application to Interfacial Incompressible Fluid Flow. SIAM Journal of Scientific Computing, 20, 1165-1191.
 Shin, S. and Juric, D. (2002) Modelling Three-Dimensional Multiphase Flow Using a Level Contour Reconstruction Method for Front Tracking without Connectivity. Journal of Computational Physics, 180, 427-470.
 Sarthou, A., Vincent, S. and Caltagirone, J.P. (2010) A Second-Order Curvilinear to Cartesian Transformation of Immersed Interfaces and Boundaries. Application to Fictitious Domains and Multiphase Flows. Computers & Fluids, 46, 422-428.
 Svard, M. and Nordstrom, J. (2006) On the Order of Accuracy for Difference Approximations of Initial-Boundary Value Problems. Journal of Computational Physics, 352, 219-333.
 Schenk, O. and Gartner, K. (2004) Solving Unsymmetric Sparse Systems of Linear Equations with Pardiso. Journal of Future Generation Computing Systems, 20, 475-487.
 Gustafsson, I. (1978) On First and Second Order Symmetric Factorization Methods for the Solution of Elliptic Difference Equations. Research Report, Department of Computer Sciences, Chalmersuniversity of Technology and the University of Gateborg, Gateborg.
 Saad, Y. and Schultz, M.H. (1986) Gmres: A Generalized Minimal Residual Algorithm for Solving Nonsymmetric Linear Systems. SIAM Journal on Scientific and Statistical Computing, 7, 856-869.
 Randrianarivelo, T.N., Pianet, G., Vincent, S. and Caltagirone, J.-P. (2005) Numerical Modelling of Solid Particle Motion Using a New Penalty Method. International Journal for Numerical Methods in Fluids, 47, 1245-1251.
 Stéphane, V., De Motta, J.C.B., Sarthou, A., Estivalezes, J.-L., Simonin, O. and Climent, E. (2012) A Lagrangian VOF Tensorial Penalty Method for the DNS of Resolved Particle-Laden Flows. Journal of Computational Physics, 256, 582-614.