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

Authors: Haifeng Ji; Qian Zhang; Bin Zhang
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

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 problem with the given total electric quantity on the interface. The proposed method 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.


---

## Body

## 1. Introduction

We consider an isolated conductor that is placed near other charges or in an external electric field. For simplicity, we take a rectangular domainΩ
∈
R
2 as the computational domain. The conductor occupying Ω
- is included in the domain Ω. The boundary of Ω
- is denoted by Γ which is also called an interface in this paper. We assume that Γ
∩
∂
Ω
≠
∅. Let Ω
+
=
Ω
\
Ω
-, let n be the unit normal vector of Γ pointing from Ω
- to Ω
+, and let ε
+ be the permittivity of Ω
+. The modeling of the electrostatic field leads to the interface problem [1],(1)
-
∇
·
∇
ϕ
-
=
0
,
i
n
Ω
-
,
-
∇
·
ε
+
∇
ϕ
+
=
f
+
x
,
y
,
i
n
Ω
+
,together with the interface jump conditions on the interface Γ,(2)
ϕ
+
-
ϕ
-
=
0
,
ε
+
∇
ϕ
+
·
n
-
∇
ϕ
-
·
n
=
-
q
x
,
y
,
x
,
y
∈
Γ
,where the interfacial function q
(
x
,
y
) represents the charge density on the interface Γ and is also unknown; however, the total electric quantity(3)
Q
=
∫
Γ
q
x
,
y
d
sis given. For simplicity, we consider a homogeneous Dirichlet boundary condition:(4)
ϕ
x
,
y
=
0
,
x
,
y
∈
∂
Ω
.Note that other boundary conditions can be treated using standard techniques. The function f
+
(
x
,
y
) represents the external charge density of the part outside of the conductor. The solution ϕ
(
x
,
y
)
|
Ω
±
=
ϕ
±
(
x
,
y
) is the potential. Once the potential ϕ is obtained, the electrostatic field can be computed by E
=
-
∇
ϕ
(
x
,
y
). If the charge distribution q
(
x
,
y
) on the interface Γ is given, then the potential ϕ
(
x
,
y
) 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–7] is the use of structured meshes which are independent of the interface, 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 q
(
x
,
y
).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–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 Section2, 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.


## 2. Augmented Technique for Given Electric Quantity on Interfaces

By introducing an augmented variableg
=
∇
ϕ
+
·
n
-
∇
ϕ
-
·
n and using the fact that the coefficient ε
+ is a constant, the original problem can be written as(5)
-
Δ
ϕ
=
f
x
,
y
,
i
n
Ω
-
∪
Ω
+
,
ϕ
Γ
=
ϕ
+
-
ϕ
-
=
0
,
o
n
Γ
,
∇
ϕ
·
n
Γ
=
∇
ϕ
+
·
n
-
∇
ϕ
-
·
n
=
g
,
o
n
Γ
,
ϕ
x
,
y
=
0
,
o
n
∂
Ω
,with(6)
f
x
,
y
=
0
i
n
Ω
-
,
f
+
x
,
y
ε
+
i
n
Ω
+
.It is obvious that the solution of (5)-(6) is dependent on the augmented variable g. Thus, we denote the solution of (5)-(6) by ϕ
g
(
x
,
y
). Let the solution of the original model problem (1)-(4) be ϕ
∗
(
x
,
y
) and define (7)
g
∗
=
∇
ϕ
∗
Ω
+
·
n
-
∇
ϕ
∗
Ω
-
·
n
,
o
n
Γ
.Then ϕ
∗
(
x
,
y
) satisfies (5)-(6) with g
≡
g
∗. In other words, ϕ
g
∗
(
x
,
y
)
=
ϕ
∗
(
x
,
y
). From the original model problem, the augmented variable g should be chosen such that the solution of (5)-(6) ϕ
g
(
x
,
y
) satisfies(8)
ε
+
∇
ϕ
g
Ω
+
·
n
-
∇
ϕ
g
Ω
+
·
n
=
-
q
x
,
y
,
x
,
y
∈
Γ
,
Q
=
∫
Γ
q
x
,
y
d
s
. 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 E
=
-
∇
ϕ
-
=
0 inside the conductor; that is, ∇
ϕ
-
·
n
=
0 on the interface Γ. Hence, we have q
=
-
ε
+
∇
ϕ
+
·
n and ∫
Γ
∇
ϕ
+
·
n
d
s
=
-
Q
/
ε
+.In the continuous case, the augmented method is to find the solutionϕ
(
x
,
y
) of (5)-(6) and the augmented variable g is constrained by(9)
∇
ϕ
-
·
n
=
0
,
x
,
y
∈
Γ
,
(10)
∫
Γ
∇
ϕ
+
·
n
d
s
=
-
Q
ε
+
.To present the numerical method, first we partition the domain into uniform rectangles with mesh size (h
x, h
y); then we obtain the triangulation T
h with mesh size h
=
(
h
x
)
2
+
(
h
y
)
2 by cutting the rectangles along one of diagonals in the same direction. We call an element T an interface element if Γ intersects T; otherwise, we call T a noninterface element. The sets of all interface elements and noninterface elements are denoted by T
h
i
n
t and T
h
n
o
n, respectively. The interface Γ is approximated by Γ
h, the union of the line segments connecting the intersections of the interface and the edges of elements. That is,(11)
Γ
h
=
⋃
i
=
1
N
Γ
Γ
h
i
,
Γ
h
i
=
x
i
Γ
x
i
+
1
Γ
¯
,
x
1
Γ
=
x
N
Γ
+
1
Γ
,where x
i
Γ is the intersection of the interface and edges of elements. Let Ω
h
- be the domain with Γ
h as its boundary, and Ω
h
+
=
Ω
\
Ω
h
-. There is a small region around the interface Γ,(12)
Γ
r
=
Ω
-
Ω
+
∩
Ω
h
+
-
Ω
-
∩
Ω
h
-
, whose area is of order O
(
h
3
). And the distance between Γ and Γ
h satisfies d
(
Γ
,
Γ
h
)
≤
C
h
2.In the discrete case, the augmented variable is piecewise constant and is defined at line segmentsΓ
h
i on interface elements; that is,(13)
G
=
G
1
,
G
2
,
…
,
G
N
Γ
T
,
G
i
=
g
Γ
h
i
,
i
=
1,2
,
…
,
N
Γ
. In other words, on the interface Γ
h, the augmented variable g is approximated by the piecewise constant function:(14)
g
h
x
,
y
=
G
i
,
i
f
x
,
y
∈
Γ
h
i
.In the numerical method, if the vectorG (or g
h 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 ϕ
h
(
x
,
y
).


## 3. Immersed Finite Element for Nonhomogeneous Interface Jump Conditions

In the IFE method, the interface jump conditions[
ϕ
]
Γ
=
0 and [
∇
ϕ
·
n
]
Γ
=
g are used to construct the discrete trial function space S
h. For any w
h
∈
S
h, the finite dimensional function w
h is piecewise linear on each element and is broken along Γ
h to satisfy [
w
h
]
=
0 and [
w
h
]
Γ
=
g
h on Γ
h. For the test function space, we use the standard P
1 conforming finite element space V
h
∈
H
0
1
(
Ω
). It will be shown later that the IFE function w
h
∈
S
h can be decomposed as(15)
w
h
=
v
h
+
ϕ
h
J
,where v
h
∈
V
h and ϕ
h
J is a piecewise linear function that is nonzero only on noninterface elements. Note that the function ϕ
h
J depends on the augmented variable G.Given the augmented variableG, the IFE method for (5)-(6) is to find ϕ
h
L
∈
V
h such that(16)
∫
Ω
∇
ϕ
h
L
·
∇
v
h
d
x
d
y
=
∫
Ω
+
f
+
v
h
d
x
d
y
+
∫
Γ
h
q
h
v
h
d
s
-
∑
T
∈
T
h
i
n
t
∫
T
∇
ϕ
h
J
·
∇
v
h
d
x
d
y
,
∀
v
h
∈
V
h
.The discrete solution is ϕ
h
=
ϕ
h
L
+
ϕ
h
J. 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 G is known.


### 3.1. Construction of the Functionϕ
h
J

First, we describe the functionw
h in the space S
h in detail. On a non-interface element, w
h 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 ∆
A
B
C (see Figure 1 for an illustration), the function w
h
∈
S
h is constructed as the following piecewise linear function:(17)
w
h
x
=
w
h
+
x
=
a
+
+
b
+
x
+
c
+
y
,
x
=
x
,
y
∈
T
+
,
w
h
-
x
=
a
-
+
b
-
x
+
c
-
y
,
x
=
x
,
y
∈
T
-
.The coefficients a
±
,
b
±
,
c
± are chosen such that(18)
w
h
A
=
V
1
,
w
h
B
=
V
2
,
w
h
C
=
V
2
,
w
h
+
D
-
w
h
-
D
=
0
,
w
h
+
E
-
w
h
-
E
=
0
,
∇
w
h
+
·
n
D
E
¯
-
∇
w
h
-
·
n
D
E
¯
=
g
h
D
E
¯
, where the vector V
i represents the degrees of freedom. Obviously, the function w
h can be decomposed as(19)
w
h
=
v
h
+
ϕ
h
J
,where v
h
∈
V
h and ϕ
h
J satisfy(20)
ϕ
h
J
x
=
a
~
+
+
b
~
+
x
+
c
~
+
y
,
x
=
x
,
y
∈
T
+
,
a
~
-
+
b
~
-
x
+
c
~
-
y
,
x
=
x
,
y
∈
T
-
,with coefficients a
~
±
,
b
~
±
,
c
~
± chosen such that(21)
ϕ
h
J
A
=
0
,
ϕ
h
J
B
=
0
,
ϕ
h
J
C
=
0
,
ϕ
h
J
D
=
0
,
ϕ
h
J
E
=
0
,
∇
ϕ
h
J
Ω
+
·
n
D
E
¯
-
∇
ϕ
h
J
Ω
-
·
n
D
E
¯
=
g
h
D
E
¯
.Figure 1
An example of the interface element.


## 3.1. Construction of the Functionϕ
h
J

First, we describe the functionw
h in the space S
h in detail. On a non-interface element, w
h 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 ∆
A
B
C (see Figure 1 for an illustration), the function w
h
∈
S
h is constructed as the following piecewise linear function:(17)
w
h
x
=
w
h
+
x
=
a
+
+
b
+
x
+
c
+
y
,
x
=
x
,
y
∈
T
+
,
w
h
-
x
=
a
-
+
b
-
x
+
c
-
y
,
x
=
x
,
y
∈
T
-
.The coefficients a
±
,
b
±
,
c
± are chosen such that(18)
w
h
A
=
V
1
,
w
h
B
=
V
2
,
w
h
C
=
V
2
,
w
h
+
D
-
w
h
-
D
=
0
,
w
h
+
E
-
w
h
-
E
=
0
,
∇
w
h
+
·
n
D
E
¯
-
∇
w
h
-
·
n
D
E
¯
=
g
h
D
E
¯
, where the vector V
i represents the degrees of freedom. Obviously, the function w
h can be decomposed as(19)
w
h
=
v
h
+
ϕ
h
J
,where v
h
∈
V
h and ϕ
h
J satisfy(20)
ϕ
h
J
x
=
a
~
+
+
b
~
+
x
+
c
~
+
y
,
x
=
x
,
y
∈
T
+
,
a
~
-
+
b
~
-
x
+
c
~
-
y
,
x
=
x
,
y
∈
T
-
,with coefficients a
~
±
,
b
~
±
,
c
~
± chosen such that(21)
ϕ
h
J
A
=
0
,
ϕ
h
J
B
=
0
,
ϕ
h
J
C
=
0
,
ϕ
h
J
D
=
0
,
ϕ
h
J
E
=
0
,
∇
ϕ
h
J
Ω
+
·
n
D
E
¯
-
∇
ϕ
h
J
Ω
-
·
n
D
E
¯
=
g
h
D
E
¯
.Figure 1
An example of the interface element.


## 4. Enforce the Total Electric Quantity on the Interface

In the discrete case, constraints (9) and (10) are replaced by(22)
∇
ϕ
h
-
·
n
Γ
h
=
0
,
x
,
y
∈
Γ
h
,
(23)
∫
Γ
h
∇
ϕ
h
+
·
n
Γ
h
d
s
=
-
Q
ε
+
.In matrix-vector form, the discretization (16) and (22) can be written as(24)
A
B
E
K
Φ
G
=
F
1
F
2
.For (23), we have the vector form a
T
Φ
+
b
T
G
=
-
Q
/
ε
+, where a and b are vectors. To enforce this constraint, we augmented the linear system with the equation and a Lagrangian multiplier μ to get(25)
A
B
a
E
K
b
a
T
b
T
0
Φ
G
μ
=
F
1
F
2
-
Q
ε
+
.We use the GMRES iteration method to solve the augmented variable G 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.


## 5. 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Ω
-
=
{
(
x
,
y
)
∈
R
2
:
x
2
+
y
2
<
r
2
} that is placed in an externally applied field E
0 along the x direction. The external charge density f
(
x
,
y
)
=
0. It is easy to verify that the exact solution is(26)
ϕ
+
x
,
y
=
-
x
E
0
+
r
2
x
E
0
x
2
+
y
2
,
ϕ
-
x
,
y
=
0
, and on the interface(27)
q
x
,
y
=
2
x
E
0
ε
+
r
,
Q
=
0
. We choose E
0
=
1, r
=
0.3, and ε
=
1 in this example. Not only the potential ϕ
h
(
x
,
y
) but also the electric field E
h
(
x
,
y
)
=
-
∇
ϕ
h
(
x
,
y
) are computed by using the proposed method. We choose the square domain Ω
=
[
-
1,1
]
×
[
-
1,1
] as the computational domain. We first partition the domain into N
2 congruent squares, and then we obtain the triangulation T
h by cutting the squares along one of diagonals in the same direction. We compute errors in L
2 for the numerical potential ϕ
h and the numerical electric field E
h and estimate the convergence rate by using(28)
R
a
t
e
=
l
o
g
2
e
2
h
L
2
e
h
L
2
,
e
h
=
ϕ
h
-
ϕ
o
r
e
h
=
E
h
-
E
. Numerical results reported in Table 1 show that the proposed method achieves first-order convergence in the L
2 norm for both the potential and the electric field. The numerical charge distribution q
h
(
x
,
y
) on the interface, the numerical potential ϕ
h
(
x
,
y
), and the numerical electric field E
h
(
x
,
y
) obtained by the proposed method with N
=
64 are plotted in Figures 2, 3, and 4, respectively.Table 1
A grid refinement analysis for Example1 for the proposed method.
N
ϕ
h
-
ϕ
L
2
Rate
E
h
-
E
L
2
Rate
16
0.59718
D
-
01
0.18051
D
+
00
32
0.13765
D
-
01
2.11712
0.79071
D
-
01
1.19087
64
0.11122
D
-
01
0.30758
0.44268
D
-
01
0.83689
128
0.44936
D
-
02
1.30749
0.21733
D
-
01
1.02640
256
0.26655
D
-
02
0.75350
0.11267
D
-
01
0.94776
512
0.11016
D
-
02
1.27473
0.53901
D
-
02
1.06372Figure 2
Exact charge distributionq
(
x
,
y
) (blue) and the numerical charge distribution q
h
(
x
,
y
) on the interface (red) obtained by the proposed method with N
=
64 for Example 1.Figure 3
Numerical potentialϕ
h
(
x
,
y
) obtained by the proposed method with N
=
64 for Example 1.Figure 4
Numerical electric fieldE
h
(
x
,
y
) obtained by the proposed method with N
=
64 for Example 1.Example 2.
ConsiderΩ
-
=
{
(
x
,
y
)
∈
R
2
:
(
x
-
x
0
)
2
+
(
y
-
y
0
)
2
<
r
0
2
} with x
0
=
-
1, y
0
=
0, and r
0
=
1.6. The computational domain is Ω
=
[
-
10,10
]
×
[
-
10,10
]. The external charge density is set to be f
(
x
,
y
)
=
-
1 inside the rectangular [
4.75,5.25
]
×
[
-
0.25,0.25
] and f
(
x
,
y
)
=
0 everywhere else. Dirichlet boundary condition ϕ
=
0 is applied to the left boundary and Neumann boundary condition ∂
ϕ
/
∂
n
=
0 is applied on the right, bottom, and top boundaries. We choose ε
=
5 and the total electric quantity Q
=
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 E
=
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.Figure 5
Numerical charge distributionq
h
(
x
,
y
) on the interface (red) obtained by the proposed method with N
=
256 for Example 2.Figure 6
Numerical potentialϕ
h
(
x
,
y
) obtained by the proposed method with N
=
256 for Example 2.Figure 7
Numerical electric fieldE
h
(
x
,
y
) obtained by the proposed method with N
=
256 for Example 2.Example 3.
Consider a conductor with complicated boundary(29)
Ω
-
=
x
,
y
∈
R
2
:
x
2
+
y
2
-
r
0
+
λ
s
i
n
ω
θ
≤
0
, where (30)
θ
=
a
r
c
c
o
s
x
r
,
i
f
x
≥
0
,
2
π
-
a
r
c
c
o
s
y
r
,
i
f
y
<
0
. We choose ω
=
5, λ
=
0.1, and r
0
=
0.5. The external charge density is set to be f
(
x
,
y
)
=
0. Dirichlet boundary condition ϕ
=
1 is applied to the left boundary and Neumann boundary condition ∂
ϕ
/
∂
n
=
0 is applied on the rest of boundaries. We choose ε
=
10, Q
=
20, and Ω
=
[
-
1,1
]
×
[
-
1,1
] in this example. Similar numerical results are plotted in Figures 8, 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 E
h
≈
0 inside the conductor.Figure 8
Numerical charge distributionq
h
(
x
,
y
) on the interface (red) obtained by the proposed method with N
=
256 for Example 3.Figure 9
Numerical potentialϕ
h
(
x
,
y
) obtained by the proposed method with N
=
256 for Example 3.Figure 10
Numerical electric fieldE
h
(
x
,
y
) obtained by the proposed method with N
=
256 for Example 3.


## 6. 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.



---

*Source: 7521273-2018-05-14.xml*