# A Consistent Immersed Finite Element Method for the Interface Elasticity Problems

Authors: Sangwon Jin; Do Y. Kwak; Daehyeon Kyeong
Journal: Advances in Mathematical Physics
 (2016)

---

## Abstract

We propose a new scheme for elasticity problems having discontinuity in the coefficients. In the previous work (Kwak et al., 2014), the authors suggested a method for solving such problems by finite element method using nonfitted grids. The proposed method is based on theP
1-nonconforming finite element methods with stabilizing terms. In this work, we modify the method by adding the consistency terms, so that the estimates of consistency terms are not necessary. We show optimal error estimates in H
1 and divergence norms under minimal assumptions. Various numerical experiments also show optimal rates of convergence.


---

## Body

## 1. Introduction

Linear elasticity equations are important governing equations in continuum mechanics, they describe how solid objects are deformed when external forces are applied on them. In particular, elasticity interface problems are important in various fields such as solid mechanics, material sciences, and biological sciences. There are various numerical methods to solve these problems such as finite difference methods, finite volume methods, and finite element methods; see [1–6] and references therein. But the problems involving composite materials lead to the discontinuity in the coefficients of the governing equations. In this case, there are two types of numerical methods from the point of view of mesh generation. In the first type, people solve the problems by using fitted grid which is created to align with the interface [7, 8]; another type is to use an unfitted (uniform) grid which is generated independently of the interface [9–15].Fitted grid method is well known and the most widely used accurate method but it has weaknesses for the time dependent problems which include moving interface. In moving interface problems, we need to regenerate the appropriate mesh for the interface which changes according to time. Therefore, it is not adequate to apply to the moving interface problems since this mesh generation has considerable computational costs.However, using a fixed grid has an obvious advantage that we can use the mesh in the previous time step, in the case when the interface changes over time. Therefore, it is suitable for moving interface problems. There are two classes of methods which belong to this type. One is the extended finite element method (XFEM) and the other is the immersed finite element method (IFEM). In the XFEM, they use enrichment basis functions in addition to the standard finite element basis. Therefore, they have more degrees of freedom than the standard FEM. In IFEM, the degrees of freedom are the same as the standard FEM basis; instead we modify the shape functions so as to satisfy the interface conditions along the interior interface. Recently, Lin et al. [9] have suggested a numerical method for solving elasticity problem with an interface using rotated Q
1-nonconforming finite element on uniform grids and Kwak et al. [10] proved the optimal error estimate for P
1-nonconforming finite element on triangular grids under an some extra regularity that the stress component belongs to H
2.In this paper, we modify the numerical scheme studied in [10] for the elasticity problem with an interface. We add consistency terms in the bilinear form to avoid the consistency term estimate. As a result, the extra regularity assumption which is necessary to estimate the consistency term error in the previous work [10] can be avoided. We prove optimal error estimates by the standard framework of finite element error analysis, which includes proving the coercivity and boundedness of the bilinear form.The rest of our paper is organized as follows. In Section2, we introduce a problem with interface along which Laplace-Young condition holds. In Section 3, we review the vector basis functions introduced in [10] based on the P
1 nonconforming elements satisfying the interface conditions and discretize the problem using (uniform) triangular grids. In Section 4, we prove the approximation property of our finite element space and the coercivity of our variational form. Next, we prove optimal H
1 and divergence norm error estimates. Finally, numerical experiments are presented in Section 5 for various Lamé constants and for shape interfaces.


## 2. Preliminaries

LetΩ be a connected, convex polygonal domain in R
2 which is divided into two subdomains Ω
+ and Ω
- by C
2 interface Γ
=
∂
Ω
+
∩
∂
Ω
-; see Figure 1. We assume the subdomains Ω
+ and Ω
- are occupied by two different elastic materials. For a differentiable function v
=
(
v
1
,
v
2
) and a tensor τ
=
τ
11
τ
12
τ
21
τ
22, we let (1)
∇
v
=
∂
v
1
∂
x
∂
v
1
∂
y
∂
v
2
∂
x
∂
v
2
∂
y
,
div
⁡
τ
=
∂
τ
11
∂
x
+
∂
τ
12
∂
y
∂
τ
21
∂
x
+
∂
τ
22
∂
y
.Then the displacement u
=
(
u
1
,
u
2
) of the elastic body under an external force satisfies the Navier-Lamé equation as follows:(2)
-
div
⁡
σ
u
=
f
i
n
Ω
s
,
s
=
+
,
-
,
(3)
u
Γ
=
0
,
(4)
σ
u
·
n
Γ
=
0
,
(5)
u
=
0
o
n
∂
Ω
,where(6)
σ
u
=
2
μ
ϵ
u
+
λ
t
r
ϵ
u
δ
,
ϵ
u
=
1
2
∇
u
+
∇
u
Tare the stress tensor and the strain tensor, respectively, n is outward unit normal vector, δ is the identity tensor, and f
∈
(
L
2
(
Ω
)
)
2 is the external force. Here (7)
λ
=
E
ν
1
+
ν
1
-
2
ν
,
μ
=
E
2
1
+
νare the Lamé constants, satisfying 0
<
μ
1
<
μ
<
μ
2 and 0
<
λ
<
∞, and E is Young’s modulus and ν is the Poisson ratio. When the parameter λ
→
∞, this equation describes the behavior of nearly incompressible material. Since the material properties are different in each region, we set Lamé constants μ
=
μ
s
,
λ
=
λ
s
o
n
Ω
s for s
=
+
,
-.Figure 1
A domainΩ with interface.Multiplyingv
∈
(
H
0
1
(
Ω
)
)
2 and applying Green’s identity in each domain Ω
s, we obtain(8)
∫
Ω
s
2
μ
s
ϵ
u
:
ϵ
v
d
x
+
∫
Ω
s
λ
s
div
⁡
u
div
⁡
v
d
x
-
∫
∂
Ω
s
σ
u
n
·
v
d
s
=
∫
Ω
s
f
·
v
d
x
,where (9)
ϵ
u
:
ϵ
v
=
∑
i
,
j
=
1
2
ϵ
i
j
u
ϵ
i
j
v
.Summing over s
=
+
,
- and applying the interior traction condition (4), we obtain the following weak form:(10)
a
u
,
v
=
f
,
v
,where(11)
a
u
,
v
=
∫
Ω
2
μ
ϵ
u
:
ϵ
v
d
x
+
∫
Ω
λ
div
⁡
u
div
⁡
v
d
x
.As usual, (
·
,
·
) denotes the L
2
(
Ω
) inner product. Then we have the following result [5, 13].Theorem 1.
There exists a unique solutionu
∈
(
H
0
1
(
Ω
)
)
2 of (2)–(5) satisfying and u
∈
(
H
2
(
Ω
s
)
)
2, s
=
+
,
-. Here, H
1
(
Ω
)
,
H
2
(
Ω
s
), and so forth are usual Sobolev spaces on respective domains and H
0
1
(
Ω
) is a subspace of H
1
(
Ω
) functions having zero trace.


## 3. A Consistent IFEM Based on Crouzeix-Raviart (CR) Element

In this section, we propose a new IFEM for elasticity equation with interface, by modifying the method suggested in [10]. More specifically, we use the modified CR element [16] and add consistency terms similar to those of the discontinuous Galerkin (DG) formulation [17], but our methods have the same degrees of freedom as those in [10], which has the same degrees of freedom as the standard CR FEM. Here we need to clarify some notations. The bracket [
·
] means the jump across the interface: (12)
u
Γ
≔
u
Ω
+
-
u
Ω
-
.Let {
T
h
} be a given quasi-uniform triangulations of Ω by the triangles of maximum diameter h. We allow the grid to be cut by the interface. Let the collection of all the interior edges of T
∈
T
h be denoted by E. For every e
∈
E, there are two elements T
1 and T
2 sharing e as a common edge. Let n
T
i, i
=
1,2 be the unit outward normal vector to the boundary of T
i. We choose a direction of the normal vector, say n
e
=
n
T
1 for each edge e, and fix it once and for all. For functions v defined on T
1
∪
T
2, we let {
·
} denote the average across e; that is,(13)
v
e
=
1
2
v
T
1
+
v
T
2
.We call an elementT
∈
T
h aninterface element if the interface Γ passes through the interior of T; otherwise, we call it anoninterface element. Let T
h
∗ be the collection of all interface elements. We assume the following situations which are easily satisfied when h is small enough:(i)
The interface intersects the edges of an element at no more than two points.
(ii)
The interface intersects each edge at most once, except possibly it passes through two vertices.The main idea of the IFEM for elasticity problem is to use two pieces of linear shape functions (vector form) on an interface element to construct a new shape function that satisfy the Laplace-Young condition. We set, fori
=
1,2
,
…
,
6,(14)
ϕ
^
i
x
,
y
=
ϕ
^
i
+
x
,
y
=
ϕ
^
i
1
+
ϕ
^
i
2
+
=
a
1
+
+
b
1
+
x
+
c
1
+
y
a
2
+
+
b
2
+
x
+
c
2
+
y
,
x
,
y
∈
T
+
,
ϕ
^
i
-
x
,
y
=
ϕ
^
i
1
-
ϕ
^
i
2
-
=
a
1
-
+
b
1
-
x
+
c
1
-
y
a
2
-
+
b
2
-
x
+
c
2
-
y
,
x
,
y
∈
T
-and require these functions that satisfy the nodal value conditions (edge average), continuity, and jump conditions along the interface (see Figure 2):(15)
ϕ
^
i
1
¯
∣
e
j
=
δ
i
j
,
j
=
1,2
,
3
,
ϕ
^
i
2
¯
∣
e
j
=
δ
i
-
3
j
,
j
=
1,2
,
3
,
ϕ
^
i
D
=
0
,
ϕ
^
i
E
=
0
,
σ
ϕ
^
i
·
n
D
E
¯
=
0
.It is shown in [10] that these twelve conditions uniquely determine the basis functions ϕ
^
i, i
=
1
,
…
,
6.Figure 2
A typical interface triangle.We denote byN
^
h
(
T
) the space of functions generated by ϕ
^
i, i
=
1
,
…
,
6 constructed above. Using this local finite element space, we define the globalimmersed finite element space
N
^
h
(
Ω
) by (16)
N
^
h
Ω
=
ϕ
^
∈
N
^
h
T
i
f
T
∈
T
h
∗
,
ϕ
^
∈
N
h
T
i
f
T
∉
T
h
∗
;
∫
e
ϕ
^
∂
T
1
d
s
=
∫
e
ϕ
^
∂
T
2
d
s
;
i
f
T
1
a
n
d
T
2
s
h
a
r
e
a
n
e
d
g
e
e
∫
∂
T
∩
∂
Ω
ϕ
^
d
s
=
0
.First we note the identity(17)
a
b
-
c
d
=
1
2
a
+
c
b
-
d
+
a
-
c
b
+
d
.Multiplying (2) by v
∈
(
H
1
(
T
)
)
2, integration by parts and using (17), we see(18)
∑
T
∈
T
h
∫
T
2
μ
ϵ
u
:
ϵ
v
d
x
+
∫
T
λ
div
⁡
u
div
⁡
v
d
x
-
∑
e
∈
E
∫
e
σ
u
·
n
e
e
v
e
d
s
=
∫
Ω
f
·
v
d
x
.We now propose an IFEM scheme for (2)–(5) based on the form (18).CRIFEM-C. Find u
h
∈
N
^
h
(
Ω
) such that(19)
a
h
u
h
,
v
h
=
f
,
v
h
,
∀
v
h
∈
N
^
h
Ω
,where(20)
a
h
u
,
v
≔
∑
T
∈
T
h
∫
T
2
μ
ϵ
u
:
ϵ
v
d
x
+
∫
T
λ
div
⁡
u
div
⁡
v
d
x
+
∑
e
∈
E
∫
e
τ
h
u
v
d
s
(21)
-
∑
e
∈
E
∫
e
σ
u
·
n
·
v
d
s
-
ϵ
∑
e
∈
E
∫
e
σ
v
·
n
·
u
d
s
,
∀
u
,
v
∈
H
h
Ω
.The variations of CRIFEM-C scheme are motivated to IIPG, SIPG, and NIPG of DG methods (ϵ
=
0, ϵ
=
1, and ϵ
=
-
1, resp.). In this paper we prove the case ϵ
=
1 only. Other cases can be similarly proved.If the last two terms in (21) are not present, it is the same as the CRIFEM studied in [10]. The difference between our new scheme and CRIFEM lies in the consistency terms.


## 4. Error Analysis

We introduce necessary function spaces and norms. Letm
≥
0 be an integer. For any domain D, we let H
m
(
D
) be the usual Sobolev space with (semi)norms denoted by |
·
|
m
,
D and ·
m
,
D. For m
=
1,2 and any domain D
=
T
(
∈

T
h
) or D
=
Ω, let (22)
H
~
m
D
2
≔
u
=
u
1
,
u
2
∈
H
m
-
1
D
2
:
u
D
∩
Ω
s
∈
H
m
D
∩
Ω
s
2
,
s
=
+
,
-with norms (23)
u
H
~
m
D
2
≔
u
m
,
D
∩
Ω
+
2
+
u
m
,
D
∩
Ω
-
2
,
u
H
~
m
D
2
≔
u
m
-
1
,
D
2
+
u
H
~
m
D
2
.When a finite element triangulation{
T
h
} is involved, the norm u
H
~
m
(
D
)
2 is understood as piecewise: (
∑
T
∈
T
h
‍
u
H
~
m
(
T
)
2
)
1
/
2 and denote it by u
m
,
h. Let H
h
(
Ω
)
≔
(
H
0
1
(
Ω
)
)
2
+
N
^
h
(
Ω
). We need subspaces of (
H
~
2
(
T
)
)
2 and (
H
~
2
(
Ω
)
)
2 satisfying the jump conditions: (24)
H
~
Γ
2
T
2
≔
u
∈
H
~
2
T
2
,
σ
u
·
n
Γ
∩
T
=
0
,
H
~
Γ
2
Ω
2
≔
u
∈
H
0
1
Ω
2
,
u
T
∈
H
~
Γ
2
T
2
,
∀
T
∈
T
h
.At last, we introduce the following mesh dependent energy-like norms(25)
v
a
h
2
≔
∑
T
∈
T
h
v
a
,
T
2
+
∑
e
∈
E
∫
e
h
σ
v
·
n
2
d
s
+
∫
e
τ
h
v
2
d
s
,where (26)
v
a
,
T
2
=
∫
T
2
μ
ϵ
v
:
ϵ
v
d
x
+
∫
T
λ
div
⁡
v
2
d
x
.Throughout the paper, the constantsC
,
C
0
,
C
1, and so forth are generic constants independent of the mesh size h and functions u
,
v but may depend on the problem data μ
,
λ
,
f
, and Ω and are not necessarily the same on each occurrence.


### 4.1. Approximation Property ofN
^
h
(
T
)

We need an interpolation operator: for anyu
∈
(
H
1
(
T
)
)
2, we define I
h
u
∈
N
^
h
(
T
) using the average of u on each edge of T by (27)
∫
e
i
I
h
u
d
s
=
∫
e
i
u
d
s
,
i
=
1,2
,
3and call I
h
u theinterpolant of u in N
^
h
(
T
). We then define I
h
u for u
∈
(
H
1
(
Ω
)
)
2 by (
I
h
u
)
|
T
=
I
h
(
u
|
T
). The following result is shown in [10].Proposition 2.
For anyu
∈
(
H
~
Γ
2
(
Ω
)
)
2, there exists a constant C
>
0 such that for m
=
0,1
(28)
u
-
I
h
u
m
,
h
+
m
·
λ
div
⁡
u
-
I
h
u
L
2
Ω
≤
C
h
2
-
m
u
H
~
2
Ω
+
m
·
λ
M
div
⁡
u
H
~
1
Ω
,
u
-
I
h
u
m
,
h
≤
C
h
2
-
m
u
H
~
2
Ω
,where λ
M
≔
m
a
x
s
=
+
,
-
λ
s.Proposition 3.
Letu be the solution of (2). We have(29)
u
-
I
h
u
a
h
≤
C
h
u
H
~
2
Ω
+
λ
M
div
⁡
u
H
~
1
Ω
,for some constant C
>
0.Proof.
Since the estimates of other terms are given in [10], it suffices to show that the term ∑
e
∈
E
‍
∫
e
‍
h
|
{
σ
(
u
-
I
h
u
)
·
n
}
|
2
d
s is bounded by u
-
I
h
u
1
,
h. It can be proved by splitting σ
(
u
-
I
h
u
)
·
n into normal and tangential component; see [18] for details. The conclusion will follow from Proposition 2.To prove the coercivity of our bilinear form, we need the following lemmas.Lemma 4 (Korn’s inequality [5, 19]).
LetD
⊂
R
n, n
=
2,3 be a simply connected domain. There exists constant C
(
D
)
>
0 such that(30)
v
H
1
D
2
≤
C
ϵ
v
L
2
D
2
+
Q
v
L
2
D
2
,
∀
v
∈
H
1
D
n
,where Q
(
v
)
≔
v
-
1
/
|
D
|
∫
D
‍
v
d
x.Corollary 5.
There exists constantC
^
0
>
0 such that (31)
v
h
1
,
T
≤
C
^
0
v
h
a
,
T
,
∀
v
h
∈
N
^
T
.Proof.
By Lemma4 and approximation property, we have (32)
v
h
1
,
T
2
≤
C
ϵ
v
h
L
2
T
2
+
Q
v
h
L
2
T
2
≤
C
v
h
a
,
T
2
+
C
h
v
h
1
,
T
2
.Hence, we get the result.Lemma 6.
There exists constantC
^
1
>
0 such that (33)
h
σ
v
h
·
n
0
,
e
2
≤
C
^
1
∇
v
h
0
,
T
2
,
∀
v
h
∈
N
^
T
,where e
⊂
∂
T is an edge of T.Proof.
It can be obtained by the similar technique in [18].Now we showa
h-form is bounded and coercive.Lemma 7.
If we chooseτ large enough, there exists constant C
0
,
C
1
>
0 such that the following holds:(34)
a
h
u
,
v
≤
C
0
u
a
h
v
a
h
,
∀
u
,
v
∈
H
h
,
C
1
v
h
a
h
2
≤
a
h
v
h
,
v
h
,
∀
v
h
∈
N
^
h
.Proof.
Boundedness is trivial by definition ofa
h-form. Using Corollary 5 and using Lemma 6 we have (35)
∑
e
∈
E
∫
e
σ
v
h
·
n
v
h
d
s
≤
∑
e
∈
E
h
σ
v
h
n
0
,
e
2
1
/
2
∑
e
∈
E
h
-
1
v
h
0
,
e
2
1
/
2
≤
C
^
1
∑
T
∈
T
h
∇
v
h
0
,
T
2
1
/
2
∑
e
∈
E
h
-
1
v
h
0
,
e
2
1
/
2
≤
C
^
0
C
^
1
∑
T
∈
T
h
v
h
a
,
T
2
1
/
2
∑
e
∈
E
h
-
1
v
h
0
,
e
2
1
/
2
≤
γ
2
∑
T
∈
T
h
v
h
a
,
T
2
+
C
^
0
C
^
1
2
γ
∑
e
∈
E
h
-
1
v
h
0
,
e
2for every γ
>
0. Hence, we have (36)
a
h
v
h
,
v
h
=
∑
T
∈
T
h
v
h
a
,
T
2
+
∑
e
∈
E
∫
e
τ
h
v
h
2
d
s
-
2
∑
e
∈
E
∫
e
σ
v
h
·
n
v
h
d
s
≥
1
-
γ
∑
T
∈
T
h
v
h
a
,
T
2
+
τ
h
-
C
^
0
C
^
1
h
γ
∑
e
∈
E
∫
e
v
h
2
d
s
≥
1
-
γ
1
+
C
^
1
∑
T
∈
T
h
v
h
a
,
T
2
+
τ
h
1
-
C
^
0
C
^
1
τ
γ
∑
e
∈
E
∫
e
v
h
2
d
s
+
1
-
γ
C
^
1
1
+
C
^
1
∑
T
∈
T
h
v
h
a
,
T
2
≥
1
-
γ
1
+
C
^
1
∑
T
∈
T
h
v
h
a
,
T
2
+
τ
h
1
-
C
^
0
C
^
1
τ
γ
∑
e
∈
E
∫
e
v
h
2
d
s
+
1
-
γ
1
+
C
^
1
C
^
0
∑
e
∈
E
∫
e
h
σ
v
h
·
n
2
d
ssince (37)
∑
e
∈
E
∫
e
h
σ
v
h
·
n
2
d
s
≤
C
^
1
∑
T
∈
T
h
∇
v
h
0
,
T
2
≤
C
^
0
C
^
1
∑
T
∈
T
h
v
h
a
,
T
2
.
If we chooseγ
<
1 and τ large enough so that (
1
-
C
^
0
C
^
1
/
τ
γ
)
>
0, then with C
c
≔
min
⁡
(
1
-
γ
/
1
+
C
^
1
,
1
-
C
^
0
C
^
1
/
τ
γ
,
(
1
-
γ
)
/
(
1
+
C
^
1
)
C
^
0
), we have (38)
a
h
v
h
,
v
h
≥
C
c
v
h
a
h
2
.Now we are ready to show our main result.Theorem 8.
Letu (resp., u
h) be the solution of (2) (resp., (19)). Then we have (39)
u
-
u
h
a
h
≤
C
h
u
H
~
2
Ω
+
λ
M
div
⁡
u
H
~
1
Ω
.Proof.
By triangular inequality, we have(40)
u
-
u
h
a
h
≤
u
h
-
I
h
u
a
h
+
u
-
I
h
u
a
h
.From Lemma 7, it follows that (41)
c
u
h
-
I
h
u
a
h
2
≤
a
h
u
h
-
I
h
u
,
u
h
-
I
h
u
=
a
h
u
-
I
h
u
,
u
h
-
I
h
u
≤
C
0
u
h
-
I
h
u
a
h
u
-
I
h
u
a
h
.So we have (42)
u
h
-
I
h
u
a
h
≤
C
u
-
I
h
u
a
h
.By Proposition 3, we have (43)
u
-
u
h
a
h
≤
C
u
-
I
h
u
a
h
≤
C
h
u
H
~
2
Ω
+
λ
M
div
⁡
u
H
~
1
Ω
.Remark 9.
If, in addition, the following relation(44)
2
μ
u
H
~
2
Ω
+
λ
div
⁡
u
H
~
1
Ω
≤
C
f
L
2
Ωholds, then the result of Theorem 8 becomes (45)
u
-
u
h
a
h
≤
C
h
f
L
2
Ω
.It means that the our estimate holds uniformly when λ
→
∞. Furthermore, one can show the following L
2 estimate by the standard duality argument: (46)
u
-
u
h
L
2
Ω
≤
C
h
2
f
L
2
Ω
.


## 4.1. Approximation Property ofN
^
h
(
T
)

We need an interpolation operator: for anyu
∈
(
H
1
(
T
)
)
2, we define I
h
u
∈
N
^
h
(
T
) using the average of u on each edge of T by (27)
∫
e
i
I
h
u
d
s
=
∫
e
i
u
d
s
,
i
=
1,2
,
3and call I
h
u theinterpolant of u in N
^
h
(
T
). We then define I
h
u for u
∈
(
H
1
(
Ω
)
)
2 by (
I
h
u
)
|
T
=
I
h
(
u
|
T
). The following result is shown in [10].Proposition 2.
For anyu
∈
(
H
~
Γ
2
(
Ω
)
)
2, there exists a constant C
>
0 such that for m
=
0,1
(28)
u
-
I
h
u
m
,
h
+
m
·
λ
div
⁡
u
-
I
h
u
L
2
Ω
≤
C
h
2
-
m
u
H
~
2
Ω
+
m
·
λ
M
div
⁡
u
H
~
1
Ω
,
u
-
I
h
u
m
,
h
≤
C
h
2
-
m
u
H
~
2
Ω
,where λ
M
≔
m
a
x
s
=
+
,
-
λ
s.Proposition 3.
Letu be the solution of (2). We have(29)
u
-
I
h
u
a
h
≤
C
h
u
H
~
2
Ω
+
λ
M
div
⁡
u
H
~
1
Ω
,for some constant C
>
0.Proof.
Since the estimates of other terms are given in [10], it suffices to show that the term ∑
e
∈
E
‍
∫
e
‍
h
|
{
σ
(
u
-
I
h
u
)
·
n
}
|
2
d
s is bounded by u
-
I
h
u
1
,
h. It can be proved by splitting σ
(
u
-
I
h
u
)
·
n into normal and tangential component; see [18] for details. The conclusion will follow from Proposition 2.To prove the coercivity of our bilinear form, we need the following lemmas.Lemma 4 (Korn’s inequality [5, 19]).
LetD
⊂
R
n, n
=
2,3 be a simply connected domain. There exists constant C
(
D
)
>
0 such that(30)
v
H
1
D
2
≤
C
ϵ
v
L
2
D
2
+
Q
v
L
2
D
2
,
∀
v
∈
H
1
D
n
,where Q
(
v
)
≔
v
-
1
/
|
D
|
∫
D
‍
v
d
x.Corollary 5.
There exists constantC
^
0
>
0 such that (31)
v
h
1
,
T
≤
C
^
0
v
h
a
,
T
,
∀
v
h
∈
N
^
T
.Proof.
By Lemma4 and approximation property, we have (32)
v
h
1
,
T
2
≤
C
ϵ
v
h
L
2
T
2
+
Q
v
h
L
2
T
2
≤
C
v
h
a
,
T
2
+
C
h
v
h
1
,
T
2
.Hence, we get the result.Lemma 6.
There exists constantC
^
1
>
0 such that (33)
h
σ
v
h
·
n
0
,
e
2
≤
C
^
1
∇
v
h
0
,
T
2
,
∀
v
h
∈
N
^
T
,where e
⊂
∂
T is an edge of T.Proof.
It can be obtained by the similar technique in [18].Now we showa
h-form is bounded and coercive.Lemma 7.
If we chooseτ large enough, there exists constant C
0
,
C
1
>
0 such that the following holds:(34)
a
h
u
,
v
≤
C
0
u
a
h
v
a
h
,
∀
u
,
v
∈
H
h
,
C
1
v
h
a
h
2
≤
a
h
v
h
,
v
h
,
∀
v
h
∈
N
^
h
.Proof.
Boundedness is trivial by definition ofa
h-form. Using Corollary 5 and using Lemma 6 we have (35)
∑
e
∈
E
∫
e
σ
v
h
·
n
v
h
d
s
≤
∑
e
∈
E
h
σ
v
h
n
0
,
e
2
1
/
2
∑
e
∈
E
h
-
1
v
h
0
,
e
2
1
/
2
≤
C
^
1
∑
T
∈
T
h
∇
v
h
0
,
T
2
1
/
2
∑
e
∈
E
h
-
1
v
h
0
,
e
2
1
/
2
≤
C
^
0
C
^
1
∑
T
∈
T
h
v
h
a
,
T
2
1
/
2
∑
e
∈
E
h
-
1
v
h
0
,
e
2
1
/
2
≤
γ
2
∑
T
∈
T
h
v
h
a
,
T
2
+
C
^
0
C
^
1
2
γ
∑
e
∈
E
h
-
1
v
h
0
,
e
2for every γ
>
0. Hence, we have (36)
a
h
v
h
,
v
h
=
∑
T
∈
T
h
v
h
a
,
T
2
+
∑
e
∈
E
∫
e
τ
h
v
h
2
d
s
-
2
∑
e
∈
E
∫
e
σ
v
h
·
n
v
h
d
s
≥
1
-
γ
∑
T
∈
T
h
v
h
a
,
T
2
+
τ
h
-
C
^
0
C
^
1
h
γ
∑
e
∈
E
∫
e
v
h
2
d
s
≥
1
-
γ
1
+
C
^
1
∑
T
∈
T
h
v
h
a
,
T
2
+
τ
h
1
-
C
^
0
C
^
1
τ
γ
∑
e
∈
E
∫
e
v
h
2
d
s
+
1
-
γ
C
^
1
1
+
C
^
1
∑
T
∈
T
h
v
h
a
,
T
2
≥
1
-
γ
1
+
C
^
1
∑
T
∈
T
h
v
h
a
,
T
2
+
τ
h
1
-
C
^
0
C
^
1
τ
γ
∑
e
∈
E
∫
e
v
h
2
d
s
+
1
-
γ
1
+
C
^
1
C
^
0
∑
e
∈
E
∫
e
h
σ
v
h
·
n
2
d
ssince (37)
∑
e
∈
E
∫
e
h
σ
v
h
·
n
2
d
s
≤
C
^
1
∑
T
∈
T
h
∇
v
h
0
,
T
2
≤
C
^
0
C
^
1
∑
T
∈
T
h
v
h
a
,
T
2
.
If we chooseγ
<
1 and τ large enough so that (
1
-
C
^
0
C
^
1
/
τ
γ
)
>
0, then with C
c
≔
min
⁡
(
1
-
γ
/
1
+
C
^
1
,
1
-
C
^
0
C
^
1
/
τ
γ
,
(
1
-
γ
)
/
(
1
+
C
^
1
)
C
^
0
), we have (38)
a
h
v
h
,
v
h
≥
C
c
v
h
a
h
2
.Now we are ready to show our main result.Theorem 8.
Letu (resp., u
h) be the solution of (2) (resp., (19)). Then we have (39)
u
-
u
h
a
h
≤
C
h
u
H
~
2
Ω
+
λ
M
div
⁡
u
H
~
1
Ω
.Proof.
By triangular inequality, we have(40)
u
-
u
h
a
h
≤
u
h
-
I
h
u
a
h
+
u
-
I
h
u
a
h
.From Lemma 7, it follows that (41)
c
u
h
-
I
h
u
a
h
2
≤
a
h
u
h
-
I
h
u
,
u
h
-
I
h
u
=
a
h
u
-
I
h
u
,
u
h
-
I
h
u
≤
C
0
u
h
-
I
h
u
a
h
u
-
I
h
u
a
h
.So we have (42)
u
h
-
I
h
u
a
h
≤
C
u
-
I
h
u
a
h
.By Proposition 3, we have (43)
u
-
u
h
a
h
≤
C
u
-
I
h
u
a
h
≤
C
h
u
H
~
2
Ω
+
λ
M
div
⁡
u
H
~
1
Ω
.Remark 9.
If, in addition, the following relation(44)
2
μ
u
H
~
2
Ω
+
λ
div
⁡
u
H
~
1
Ω
≤
C
f
L
2
Ωholds, then the result of Theorem 8 becomes (45)
u
-
u
h
a
h
≤
C
h
f
L
2
Ω
.It means that the our estimate holds uniformly when λ
→
∞. Furthermore, one can show the following L
2 estimate by the standard duality argument: (46)
u
-
u
h
L
2
Ω
≤
C
h
2
f
L
2
Ω
.


## 5. Numerical Results

In this section we present numerical examples. The domain isΩ
=
(
-
1,1
)
×
(
-
1,1
). The interface is the zero set of L
(
x
,
y
)
=
x
2
+
y
2
-
r
0
2. Let Ω
+
=
Ω
∩
{
x
,
y
∣
L
(
x
,
y
)
>
0
}, Ω
-
=
Ω
∩
{
(
x
,
y
)
∣
L
(
x
,
y
)
<
0
}. The exact solution is chosen as (47)
u
=
1
μ
x
2
+
y
2
-
r
0
2
x
,
1
μ
x
2
+
y
2
-
r
0
2
ywith various values of μ and λ. For numerical simulation we partition the domain into uniform right triangles having size h
=
2
-
k, k
=
3,4
,
….Example 1.
In this example, we test two sets of parameters and radii of the interface.(1)
We chooseμ
-
=
1
,
μ
+
=
100
,
λ
=
5
μ, and r
0
=
0.36.
(2)
We chooseμ
-
=
1
,
μ
+
=
10
,
λ
=
5
μ, and r
0
=
0.48.
Tables1 and 2 show the convergence behavior of our numerical schemes for both examples. In both cases, we see the optimal order of convergence in L
2, H
1, and divergence norms. x-components of the solution are plotted in Figures 3 and 4.Table 1
μ
-
=
1, μ
+
=
100, and λ = 5μ.
1
/
h
u
-
u
h
0
Order
u
-
u
h
1
,
h
Order
div
⁡
u
-
div
⁡
u
h
0
Order
8
1.271
e
-
3
4.230
e
-
2
3.950
e
-
2
16
3.096
e
-
4
2.038
1.979
e
-
2
1.119
1.834
e
-
2
1.107
32
6.461
e
-
5
2.261
9.780
e
-
3
1.017
9.374
e
-
3
0.968
64
1.534
e
-
5
2.075
4.899
e
-
3
0.977
4.790
e
-
3
0.969
128
3.636
e
-
6
2.076
2.432
e
-
3
1.011
2.412
e
-
3
0.990
256
8.882
e
-
7
2.033
1.214
e
-
3
1.003
1.210
e
-
3
0.995Table 2
μ
-
=
1, μ
+
=
10, and λ = 5μ.
1
/
h
u
-
u
h
0
Order
u
-
u
h
1
,
h
Order
div
⁡
u
-
div
⁡
u
h
0
Order
8
1.729
e
-
3
7.374
e
-
2
7.260
e
-
2
16
4.147
e
-
4
2.060
3.786
e
-
2
0.962
3.774
e
-
2
0.944
32
1.030
e
-
4
2.009
1.938
e
-
2
0.966
1.934
e
-
2
0.965
64
2.565
e
-
5
2.007
9.734
e
-
3
0.993
9.744
e
-
3
0.989
128
6.384
e
-
6
2.006
4.886
e
-
3
0.994
4.897
e
-
3
0.993
256
1.590
e
-
6
2.005
2.447
e
-
3
0.998
2.454
e
-
3
0.997Figure 3
x-component, μ
-
=
1
,
μ
+
=
100
,
λ
=
5
μ.Figure 4
x-component, μ
-
=
1
,
μ
+
=
10
,
λ
=
5
μ.Example 2 (nearly incompressible case).
(1) We letμ
-
=
1, μ
+
=
10, λ
=
100
μ, ν
=
0.495, and r
0
=
0.7.
(2) We letμ
-
=
1, μ
+
=
10, λ
=
1000
μ, ν
=
0.4995, and r
0
=
0.6.
Tables3 and 4 show the convergence behavior. In both cases, we see the optimal order of convergence in L
2, H
1, and divergence norms. No locking phenomena occurs in both cases. Again x-components of the solution are plotted in Figures 5 and 6.Table 3
μ
-
=
1, μ
+
=
10, and λ = 100μ.
1
/
h
u
-
u
h
0
Order
u
-
u
h
1
,
h
Order
div
⁡
u
-
div
⁡
u
h
0
Order
8
7.748
e
-
3
2.156
e
-
1
1.457
e
-
1
16
2.432
e
-
3
1.672
1.126
e
-
1
0.937
7.363
e
-
2
0.985
32
7.363
e
-
4
1.724
5.533
e
-
2
1.025
3.724
e
-
2
0.983
64
2.121
e
-
4
1.796
2.735
e
-
2
1.016
1.876
e
-
2
0.989
128
5.475
e
-
5
1.954
1.349
e
-
2
1.020
9.417
e
-
3
0.994
256
1.409
e
-
5
1.959
6.689
e
-
3
1.012
4.719
e
-
3
0.997Table 4
μ
-
=
1, μ
+
=
10, and λ = 1000μ.
1
/
h
u
-
u
h
0
Order
u
-
u
h
1
,
h
Order
div
⁡
u
-
div
⁡
u
h
0
Order
8
7.799
e
-
2
1.636
e
-
0
1.059
e
-
1
16
2.418
e
-
2
1.689
9.117
e
-
1
0.844
5.514
e
-
2
0.942
32
6.848
e
-
2
1.820
4.521
e
-
1
1.012
2.817
e
-
2
0.969
64
1.860
e
-
3
1.880
2.249
e
-
1
1.007
1.417
e
-
2
0.991
128
4.848
e
-
4
1.940
1.112
e
-
1
1.016
7.110
e
-
3
0.995
256
1.237
e
-
4
1.971
5.536
e
-
2
1.006
3.563
e
-
3
0.997Figure 5
x-component, μ
-
=
1
,
μ
+
=
10
,
λ
=
100
μ.Figure 6
x-component, μ
-
=
1
,
μ
+
=
10
,
λ
=
1000
μ.Example 3 (ellipse interface case).
Next we consider examples with elliptic shaped interface. The domain is the same as above, and the interface is represented byL
(
x
,
y
)
=
x
2
/
4
+
y
2
-
r
0
2
=
0. The exact solution is chosen as (48)
u
=
1
μ
x
2
4
+
y
2
-
r
0
2
x
,
1
μ
x
2
4
+
y
2
-
r
0
2
ywith various values of μ and λ. (1)
We letμ
-
=
1
,
μ
+
=
10
,
λ
=
5
μ
,
r
0
=
0.4.
(2)
We letμ
-
=
1
,
μ
+
=
100
,
λ
=
5
μ
,
r
0
=
0.3.
Tables5 and 6 show the convergence behavior. We observe similar optimal convergence rates for all norms. Figures 7 and 8 show the x-components of the solution.Table 5
μ
-
=
1, μ
+
=
10, λ = 5μ, and elliptical interface.
1
/
h
u
-
u
h
0
Order
u
-
u
h
1
,
h
Order
div
⁡
u
-
div
⁡
u
h
0
Order
8
1.791
e
-
3
6.431
e
-
2
5.906
e
-
2
16
4.543
e
-
4
1.979
3.105
e
-
2
1.050
2.902
e
-
2
1.025
32
1.139
e
-
4
1.996
1.568
e
-
2
0.986
1.474
e
-
2
0.977
64
2.785
e
-
5
2.032
7.908
e
-
3
0.987
7.452
e
-
3
0.984
128
6.969
e
-
6
1.999
3.970
e
-
3
0.994
3.747
e
-
3
0.992
256
1.740
e
-
6
2.002
1.988
e
-
3
0.998
1.878
e
-
3
0.997Table 6
μ
-
=
1, μ
+
=
100, λ
=
5
μ, and elliptical interface.
1
/
h
u
-
u
h
0
Order
u
-
u
h
1
,
h
Order
div
⁡
u
-
div
⁡
u
h
0
Order
8
1.299e − 3
3.532e − 2
3.108e − 2
16
4.055e − 4
1.680
1.681e − 2
1.072
1.389e − 2
1.162
32
8.503e − 5
2.254
8.272e − 3
1.023
7.251e − 3
0.937
64
1.848e − 5
2.202
4.089e − 3
1.016
3.713e − 3
0.966
128
4.517e − 6
2.032
2.041e − 3
1.002
1.869e − 3
0.990
256
1.116e − 6
2.017
1.020e − 3
1.001
9.385e − 4
0.994Figure 7
x-component, μ
-
=
1
,
μ
+
=
10
,
λ
=
5
μ elliptical interface.Figure 8
x-component, μ
-
=
1
,
μ
+
=
100
,
λ
=
5
μ elliptical interface.Example 4 (unknown solution).
This last example computes a problem with unknown solution. We chooseμ
-
=
1, μ
+
=
100, ν
-
=
0.28, ν
+
=
0.4, r
0
=
0.3, and F
=
-
11
/
4
-
λ
/
μ
x
,
-
29
/
4
-
λ
/
μ
y with the same elliptical interface as in the previous example.
Figure9 shows the x-component of the computed solution.Figure 9
x-component, μ
-
=
1
,
μ
+
=
100
,
ν
-
=
0.28
,
ν
+
=
0.4
,
r
0
=
0.3 unknown solution.


## 6. Conclusion

We presented a numerical scheme using a uniform grid for the elasticity problem with an interface. This is a modification of the scheme studied in [10] by adding consistency terms. This scheme is very useful in computing problems with moving interface since we use a grid independent of the interface. We proved an optimal order of convergence in H
1-norm with minimal assumption. Numerical results show the efficiency of our method.



---

*Source: 3292487-2016-05-09.xml*