A Cartesian Grid Method for Modeling Charge Distribution on Interfaces via Augmented Technique

In the electrostatic field computations, second-order elliptic interface problems with nonhomogeneous interface jump conditions need to be solved. In realistic applications, often the total electric quantity on the interface is given. However, the charge distribution on the interface corresponding to the nonhomogeneous interface jump condition is unknown. This paper proposes a Cartesian grid method for solving the interface problemwith the given total electric quantity on the interface.The proposedmethod employs both the immersed finite element with the nonhomogeneous interface jump condition and the augmented technique. Numerical experiments are presented to show the accuracy and efficiency of the proposed method.

Note that other boundary conditions can be treated using standard techniques.The function  + (, ) represents the external charge density of the part outside of the conductor.The solution (, )| Ω ± =  ± (, ) is the potential.Once the potential  is obtained, the electrostatic field can be computed by  = −∇(, ).If the charge distribution (, ) on the interface Γ is given, then the potential (, ) can be solved efficiently by the immersed finite element (IFE) method for nonhomogeneous interface jump conditions (see [2,3]).The significance of the IFE method [4][5][6][7] is the use of structured meshes which are independent of the interface, 2 Mathematical Problems in Engineering such as Cartesian meshes.The IFE method modifies the basis function on interface elements according to the interface conditions to capture the jumps of the exact solution.The weak form and the degrees of freedom remain the same as if there was no interface.If the coefficient is a constant without jumps, then the stiffness matrix is the same as that obtained by traditional finite element for the problem without interfaces.And only the right-hand side needs to be modified according to the interface conditions.The elliptic interface problem can be solved efficiently by the IFE method with the given nonhomogeneous interface jump (, ).However, for the problem discussed in this paper, only the total electric quantity on the interface is known, not the charge density distribution on the interface which is related to the nonhomogeneous interface jump condition of this problem.In [1], the authors proposed an iterative IFE method for this interface problem.An IFE method with nonhomogeneous interface jump conditions and a standard finite element method with ghost nodes are combined to get a "Prediction-Correction-Prediction" iteration.Numerical examples in [1] show that the iterative method is convergent and can solve this problem efficiently.Note that for partial differential equations there are many other methods in the literature [8][9][10][11].In this paper, we present a new Cartesian grid method based on the IFE method and the augmented technique [12,13].By introducing the jump of the normal derivative of the exact solution as an augmented variable, we can get an efficient discretization in which the fast Fourier transform-(FFT-) based fast Poisson solver can be applied.The augmented variable is chosen such that the nonhomogeneous interface jump condition and the total electric quantity are satisfied.In the numerical method, the augmented variable is solved by using the GMRES iteration.Compared with the iterative IFE method proposed in [1], the advantage of our Cartesian grid method is that the FFT-based fast Poisson solver can be used.Numerical experiments are also provided in this paper to show the performance of the proposed method.
The rest of the paper is organized as follows.In Section 2, we describe the augmented technique for the interface problem with given electric quantity on the interface.We choose an augmented variable and rewrite the interface problem to a new one for which the FFT-based fast Poisson solver can be applied.In Section 3, we briefly recall the IFE method for the nonhomogeneous interface jump conditions, where the augmented variable is assumed to be given.In Section 4, the constraint of the total electric quantity is enforced and some implementation details about the GMRES iteration are described.Finally, some numerical examples are provided in Section 5 to show the accuracy and efficiency of the proposed method.

Augmented Technique for Given Electric Quantity on Interfaces
By introducing an augmented variable  = ∇ + ⋅n−∇ − ⋅n and using the fact that the coefficient  + is a constant, the original problem can be written as with It is obvious that the solution of ( 5)-( 6) is dependent on the augmented variable .Thus, we denote the solution of ( 5)-( 6) by   (, ).Let the solution of the original model problem ( 1)-( 4) be  * (, ) and define Then  * (, ) satisfies ( 5)-( 6) with  ≡  * .In other words,   * (, ) =  * (, ).From the original model problem, the augmented variable  should be chosen such that the solution of ( 5)-( 6)   (, ) satisfies This is the constraint for the choice of the augmented variable.When the conductor Ω − is in electrical equilibrium, the electrical potential is a constant and electric field  = −∇ − = 0 inside the conductor; that is, ∇ − ⋅ n = 0 on the interface Γ. Hence, we have  = − + ∇ + ⋅ n and ∫ Γ ∇ + ⋅ n  = −/ + .
In the continuous case, the augmented method is to find the solution (, ) of ( 5)-( 6) and the augmented variable  is constrained by To present the numerical method, first we partition the domain into uniform rectangles with mesh size (ℎ  , ℎ  ); then we obtain the triangulation T ℎ with mesh size ℎ = √(ℎ  ) 2 + (ℎ  ) 2 by cutting the rectangles along one of diagonals in the same direction.We call an element  an interface element if Γ intersects ; otherwise, we call  a noninterface element.The sets of all interface elements and noninterface elements are denoted by T int ℎ and T non ℎ , respectively.The interface Γ is approximated by Γ ℎ , the union of the line segments connecting the intersections of the interface and the edges of elements.That is, where x Γ  is the intersection of the interface and edges of elements.Let Ω − ℎ be the domain with Γ ℎ as its boundary, and There is a small region around the interface Γ, whose area is of order (ℎ 3 ).And the distance between Γ and Γ ℎ satisfies (Γ, Γ ℎ ) ≤ ℎ 2 .In the discrete case, the augmented variable is piecewise constant and is defined at line segments Γ  ℎ on interface elements; that is, In other words, on the interface Γ ℎ , the augmented variable  is approximated by the piecewise constant function: In the numerical method, if the vector  (or  ℎ in the function form) is given, then we can use the IFE method which will be described in the next section to get the discrete solution  ℎ (, ).

Immersed Finite Element for Nonhomogeneous Interface Jump Conditions
In the IFE method, the interface jump conditions [] Γ = 0 and [∇ ⋅ n] Γ =  are used to construct the discrete trial function space  ℎ .For any  ℎ ∈  ℎ , the finite dimensional function  ℎ is piecewise linear on each element and is broken along Γ ℎ to satisfy [ ℎ ] = 0 and [ ℎ ] Γ =  ℎ on Γ ℎ .For the test function space, we use the standard P 1 conforming finite element space  ℎ ∈  1 0 (Ω).It will be shown later that the IFE function  ℎ ∈  ℎ can be decomposed as where V ℎ ∈  ℎ and   ℎ is a piecewise linear function that is nonzero only on noninterface elements.Note that the function   ℎ depends on the augmented variable .Given the augmented variable , the IFE method for (5)-( 6) is to find The discrete solution is  ℎ =   ℎ +   ℎ .From (16), it is obvious that the stiffness matrix is the same as that obtained by traditional finite element for the problem without interfaces; only the right-hand side needs to be modified.Thus, we can take advantage of the fast Poisson solver to solve the system of equations efficiently when the augmented variable  is known.

Construction of the Function 𝜙 𝐽
ℎ .First, we describe the function  ℎ in the space  ℎ in detail.On a non-interface element,  ℎ is the standard linear function and the degrees of freedom are functional values on the vertices of the element.On interface elements, for example, on  (see Figure 1 for an illustration), the function  ℎ ∈  ℎ is constructed as the following piecewise linear function: The coefficients  ± ,  ± ,  ± are chosen such that where the vector   represents the degrees of freedom.Obviously, the function  ℎ can be decomposed as where V ℎ ∈  ℎ and   ℎ satisfy with coefficients ã± , b± , c± chosen such that (21)

Enforce the Total Electric Quantity on the Interface
In the discrete case, constraints ( 9) and ( 10) are replaced by In matrix-vector form, the discretization ( 16) and ( 22) can be written as For ( 23), we have the vector form   Φ +    = −/ + , where  and  are vectors.To enforce this constraint, we augmented the linear system with the equation and a Lagrangian multiplier  to get We use the GMRES iteration method to solve the augmented variable  and the Lagrangian multiplier  first and then to solve Φ by using one more fast Poisson solver.We refer the readers to Section 6.1.2 in [14] for the details about the GMRES iteration.

Numerical Experiments
In this section, we present some examples to show the accuracy and the efficiency of the proposed numerical method.First, we consider the following example in which the exact solution is given.This example is taken from [1].
Example 1.Consider a conductor Ω − = {(, ) ∈ R 2 :  2 +  2 <  2 } that is placed in an externally applied field  0 along the  direction.The external charge density (, ) = 0.It is easy to verify that the exact solution is and on the interface We choose  0 = 1,  = 0.3, and  = 1 in this example.Not only the potential  ℎ (, ) but also the electric field  ℎ (, ) = −∇ ℎ (, ) are computed by using the proposed method.We choose the square domain  the computational domain.We first partition the domain into  2 congruent squares, and then we obtain the triangulation T ℎ by cutting the squares along one of diagonals in the same direction.We compute errors in  2 for the numerical potential  ℎ and the numerical electric field  ℎ and estimate the convergence rate by using Numerical results reported in Table 1 show that the proposed method achieves first-order convergence in the  2 norm for both the potential and the electric field.The numerical charge distribution  ℎ (, ) on the interface, the numerical potential  ℎ (, ), and the numerical electric field  ℎ (, ) obtained by the proposed method with  = 64 are plotted in Figures 2, 3, and 4, respectively.We choose  = 5 and the total electric quantity  = 2 in this example.Note that there is no explicit analytical solution for this example.Numerical results are plotted in Figures 5, 6, and 7. Note that the external charge is negative in a small area and the total electric quantity is positive.In order to maintain the electric field  = 0 inside the conductor, the positive surface charges should move towards the external negative charges. Figure 5 shows that the positive charges gather on the right side of the circle.Example 3. Consider a conductor with complicated boundary where We choose  = 5,  = 0.1, and  0 = 0.5.The external charge density is set to be (, ) = 0. Dirichlet boundary condition

Conclusion
We presented a Cartesian grid method for solving the surface charge distribution problem in electromagnetism.The advantage of the method is that the used mesh does not need to be aligned with the interface.This is quite convenient for the problem with complex interfaces.The proposed method employs both the immersed finite element and the augmented technique.The GMRES iteration and FFT-based fast Poisson solver are used to solve the discrete systems.Numerical examples are provided to show the accuracy and efficiency of the proposed method.

Figure 1 :
Figure 1: An example of the interface element.

Figure 2 :
Figure 2: Exact charge distribution (, ) (blue) and the numerical charge distribution  ℎ (, ) on the interface (red) obtained by the proposed method with  = 64 for Example 1.

Figure 5 :
Figure 5: Numerical charge distribution  ℎ (, ) on the interface (red) obtained by the proposed method with  = 256 for Example 2.

Figure 8 :
Figure 8: Numerical charge distribution  ℎ (, ) on the interface (red) obtained by the proposed method with  = 256 for Example 3.

𝜙 = 1
is applied to the left boundary and Neumann boundary condition / = 0 is applied on the rest of boundaries.We choose  = 10,  = 20, and Ω = [−1, 1] × [−1, 1] in this example.Similar numerical results are plotted in Figures8, 9, and 10.From these figures, we can see that the surface charge is distributed according to the potential  = 1 on the left boundary to maintain the electric field  ℎ ≈ 0 inside the conductor.

Table 1 :
A grid refinement analysis for Example 1 for the proposed method.