Finite Element Methods for Microelectromechanical Systems

(1)

IT 09 060

Examensarbete 30 hp November 2009

Finite Element Methods for Microelectromechanical Systems

Yang Zeng

Institutionen för informationsteknologi

Department of Information Technology

(2)

(3)

Teknisk- naturvetenskaplig fakultet UTH-enheten

Besöksadress:

Ångströmlaboratoriet Lägerhyddsvägen 1 Hus 4, Plan 0 Postadress:

Box 536 751 21 Uppsala Telefon:

018 – 471 30 03 Telefax:

018 – 471 30 00 Hemsida:

http://www.teknat.uu.se/student

Abstract

Finite Element Methods for Microelectromechanical Systems

Yang Zeng

The stationary Joule heating problem is a crucial multiphysical problem for many microelectromechanical (MEMS) applications. In our paper, we derive a finite element method for this problem and introduce iterative solution-techniques to compute the numerical simulation. Further we construct an adaptive algorithm for mesh refinement based on a posteriori error estimation.

Finally, we present two numerical tests: convergences analysis of different iterative methods for distinct materials which are classified by electrical conductivities, and a test of the new adaptive refinement algorithm. All the numerical implementations have been done in MATLAB.

Tryckt av: Reprocentralen ITC IT 09 060

Examinator: Anders Jansson Ämnesgranskare: Maya Neytcheva Handledare: Axel Målqvist

(4)

(5)

List of Figures

3.1 A domain Ω and boundary ∂Ω, ∂Ω = Γ

_D

∪ Γ

_N

. . . . 9

3.2 The example of area coordinates . . . . 14

4.1 Refinement algorithms: Rivara and Regular . . . . 30

5.1 U-shape . . . . 33

5.2 Electrical conductivity of a metal . . . . 34

5.3 Electrical conductivity of a semiconductor . . . . 35

5.4 Electrical conductivity of a superconductor . . . . 36

5.5 Numbers of iterations for the metal with different ω in Jacobi and GS methods . . . . 38

5.6 Numbers of iterations for the semiconductor with different ω in Jacobi and GS methods . . . . 39

5.7 Numbers of iterations for the semiconductor with different ω in Jacobi and GS methods . . . . 39

5.8 Final refined mesh on the metal . . . . 40

5.9 Finite element approximation U on the metal . . . . 40

5.10 Finite element approximation Φ on the metal . . . . 41

5.11 Relative error energy norm for the metal . . . . 41

5.12 Final refined mesh on the semiconductor . . . . 42

5.13 Finite element approximation U on the semiconductor . . . . 42

5.14 Finite element approximation Φ on the semicondctor . . . . . 43

5.15 Relative error energy norm for the semiconductor . . . . 43

(8)

(9)

Acknowledgements

First of all, I would like to express my warmest gratitude to my supervisor, Axel M˚ alqvist, for his patience and instructive suggestions on the writing of this thesis. Without his invaluable help and generous encouragement, I would not accomplish my thesis.

Besides, I wish to thank the department of information technology of Uppsala University for giving me the opportunity to study here. I am also greatly indebted to the professors and teachers in Uppsala University, who have instructed and helped me a lot in the past two years.

The last but not the least, my thanks would go to my beloved family

and friends for their loving considerations and great confidence in me all

through these years.

(10)

(11)

Chapter 1

Introduction

Partial differential equations are used to model physical phenomena. In many important applications, several different physical processes are active at the same time. One such example is the design of microelectromechanical systems (MEMS).

MEMS are e.g. used to build sensors on micrometer scale. Here electric potential, heat transfer and mechanical stresses are coupled in a system of non-linear elliptic partial differential equations. The stationary Joule heating problem is a model problem for this application. A voltage is applied at the boundary of a device. A current that flows through the device is produced, which leads to heating of the material. The equations describing electrostatical potential and the temperature are coupled. The heat equation is driven by the electrical current and the electric conductivity depends on the temperature, which means we get coupling in both directions.

Many industrial codes are optimized for solving single physics problems, such as the two individuals equations in the Joule heating application. In order to take advantage of this infrastructure, engineers typically couple to- gether such optimized single physics solvers when solving coupled problems.

One approach is to iterate between two problems in a Jacobi or Gauss-Seidel fashion with a given initial guess. However, it is hard to predict whether the iteration will converge or not.

In this paper, we will study numerical simulation of the Joule heating

(12)

problem using the Finite Element Method, then construct an adaptive al-

gorithm for mesh refinement based on a posteriori error estimates as well as

analyze convergence of different iterative methods for distinct materials in

MEMS devices.

(13)

Chapter 2

Physical Background

2.1 Single Physics Problems

2.1.1 Heat equation

The heat equation is an important partial differential equation which de- scribes the distribution of heat (or variation in temperature) in a given region over time. The stationary heat equation reads

−∇ · k∇u = f (2.1)

where k is the thermal conductivity, u is the temperature and f is a given heat source.

2.1.2 Potential equation

Electrostatic phenomena arises from the forces that electric charges exert on each other. One of the cornerstones of electrostatics is the posing and solving of problems described by the Poisson equation

−∇ · σ(u)∇φ = g (2.2)

where σ(u) is the electrical conductivity, which is strongly dependent on

temperature, φ is the electric potential and g is a given function.

(14)

2.2 Multiphysics Problems

Joule heating is generated by the resistance of materials to electric current and presents in any electric conductor. Therefore it is crucial in many MEMS applications. Let u be the absolute temperature and φ be the electric potential in a solid electric conductor represented by a bounded domain Ω. Under steady conditions, the stationary Joule heating problem consists of the following nonlinear elliptic system

−∇ · k∇u = σ(u) · |∇φ|

²

, in Ω (2.3)

−∇ · σ(u)∇φ = f, in Ω (2.4)

with some suitable boundary conditions, where k and σ(u) are the thermal and electrical conductivities respectively which we suppose to be given positive functions.

A number of different materials can be used in MEMS technology. Silicon is the material used to create most integrated circuits used in consumer elec- tronic in the modern world. Because of availabilities of cheap high-quanlity and incorporating electronics functionality, silicon has been exploited for a wide variety of MEMS applications. Metals can also be used to create MEMS elements because of their high reliability. The devices based on so called high-temperature superconductivity have been paid more and more attention in past two decades. These different materials exhibit distinct physical properties, but we are more interested in their electrical conductivities in our paper.

2.3 Electrical Conductivities

2.3.1 Common definition

When an electrical potential difference is placed across a conductor, its mov-

able charges flow, giving rise to an electrical current. The conductivity σ is

defined as the ratio of the current density J to the electrical field strength

(15)

E

J = σE. (2.5)

Conductivity is the reciprocal of electrical resistivity ρ, and has SI unit of siemens per meter S · m

⁻¹

σ = 1

ρ . (2.6)

2.3.2 Classifications

According to the electrical conductivity, materials can be classified into the following categories:

• A conductor such as metal has high conductivity and a low resistivity;

• An insulator like glass or a vacuum has low conductivity and a high resistivity;

• The conductivity of a semiconductor is generally intermediate, but varies widely under different conditions, such as exposure of the material to electric fields or specific frequencies of light and, most important, with temperature and composition of the semiconductor material;

• Superconductivity occurs at extremely low temperature (not far from absolute zero), and materials have been found to exhibit very high electrical conductivity in this phenomenon.

Based on this classifications, we will test convergences of iterative solution- techniques when we are solving the Joule heating problem for three sorts of materials: metals, semiconductors and superconductors.

2.3.3 Temperature dependency

Electrical conductivity is strongly dependent on temperature. This depen-

dence is often expressed using a conductivity-vs-temperature graph. Let us

define σ

⁰

to be the electrical conductivity at a standard temperature u

⁰

and

(16)

α to be the temperature compensation slope. Based on the physical property of metals such that electrical conductivity decreases with the increasing temperature, we can describe this kind of conductivity-vs-temperature graph as

σ(u) =







σ

⁰

u ≤ u

⁰

σ⁰

1+α(u−u⁰)

u > u

⁰

(2.7)

where σ(u) is electrical conductivity at temperature u.

In semiconductors, electrical conductivity increases as temperature is increasing, hence its dependence can be written as

σ(u) =







σ

⁰

u ≤ u

⁰

σ

⁰

1 + α(u − u

⁰

)

u > u

⁰

(2.8)

where σ(u), σ

⁰

, α and u

⁰

have the same explanations as equation (2.7).

Finally let us consider superconductors. As we know, one of the most important physical properties is that when temperature is decreased to a very low value (not far from absolute zero), the electrical resistance of a superconductor lowers to exact zero. Electrical conductivity is the inverse of electrical resistivity, i.e. it will be extremely large in such situation. Based on this property, we can define the conductivity-vs-temperature graph for superconductors as

σ(u) =







∞ u < u

⁰

σ⁰

α(u−u⁰)

u ≥ u

⁰

(2.9)

where σ(u), σ

⁰

, α and u

⁰

are the same as equations (2.7) and (2.8).

(17)

Chapter 3

Mathematical Model and FEM

3.1 Preliminaries

First of all, let us settle some definitions and notations that will be frequently used in the paper. The scalar product (·, ·) is the ordinary L

²

= L

²

(Ω) product and k · k is the corresponding norm. A triangulation, or mesh, K of Ω is a set {K} of triangles K such that Ω = S

K∈K

K and intersection of two triangles is either a triangle edge, a triangle corner, or empty. A closed polyline ∂K of ∂Ω is a set {E} of edges E such that ∂Ω = S

E∈∂K

E. In other words, {E} is a set of sides of triangles which are on the boundary.

The boundary of a triangle K is denoted by ∂K. If an edge is shared by two triangles, we call it an interior edge, otherwise it will be on the boundary of Ω, i.e. ∂K ∪ ∂Ω = E ∈ ∂K. We denote the longest edge of triangle K by h

_K

and the length of edge E by h

_E

.

Let K be a triangle with nodes at the corners N

₁

= (x

¹₁

, x

¹₂

), N

₂

= (x

²₁

, x

²₂

), and N

₃

= (x

³₁

, x

³₂

) and let P

₁

(K) denote the vector space of linear polynomials defined on K

P

₁

(K) = {v : v(x

1

, x

2

) = c

0

+ c

1

x

1

+ c

2

x

2

, c

0

, c

1

, c

2

∈ R}. (3.1)

(18)

Let V

h

be the vector space of all continuous piecewise linear polynomials V

_h

= {v : v ∈ C(Ω), v|

_K

∈ P

₁

(K) ∀K ∈ K} (3.2) where C(Ω) denotes the space of all continuous functions on Ω and P

1

(K) is the space of linear polynomials on K as defined by (3.1). We denote the basis functions ϕ

j

(N

i

) ⊂ V

_h

, such that

ϕ

_j

(N

_i

) =







0 i = j 1 i 6= j

for i, j = 1, 2, . . . , N . Using these basis functions, all v ∈ V

_h

can be written as

v(x

₁

, x

₂

) =

N

X

i=1

α

_i

ϕ

_i

(x

₁

, x

₂

) (3.3) where the coefficients α

i

are the nodal values of the function v, that is

α

i

= v(N

i

), i = 1, 2, . . . , N. (3.4) Given a continuous function f ∈ C(K) on a triangle K with nodes at N

_i

= (x

ⁱ₁

, x

ⁱ₂

), i = 1, 2, 3, we define the interpolant πf ∈ P

₁

(K) of f as

πf =

3

X

i=1

f (N

_i

)ϕ

_i

.

We let Df and D

²

f be defined by Df = (| ∂f

∂x |

²

+ | ∂f

∂y |

²

)

^1/2

, D

²

f = (| ∂

²

f

∂x

²

|

²

+ 2| ∂

²

f

∂x∂y |

²

+ | ∂

²

f

∂y

²

|

²

)

^1/2

.

Using these notations, the following proposition for the interpolation error estimates has been shown, i.e. in [2].

Proposition 3.1. The following estimates hold X

K

kf − πf k

²_L2(K)

≤ C X

K

h

²_K

kDf k

²_L2(K)

X

K

kD(f − πf )k

²_L2(K)

≤ C X

K

kDf k

²_L2(K)

where C is a constant.

(19)

3.2 The Model Problem

Let us define a model problem for the stationary Joule heating problem with mixed boundary conditions. It can be written as

−∇ · k∇u = σ(u) · |∇φ|

²

in Ω (3.5) n · k∇u = κ

₁

(g

₁

− u) on ∂Ω (3.6)

−∇ · σ(u)∇φ = f in Ω (3.7)

φ = g

₂

on Γ

_D

(3.8)

n · σ(u)∇φ = κ

2

(g

3

− φ) on Γ

_N

(3.9) where g

1

, g

2

, g

3

and f are given functions, Ω is a two-dimensional domain with the boundary ∂Ω, ∂Ω = Γ

D

∪ Γ

_N

, and n is the outward normal of ∂Ω, as shown in Figure 3.1. Next we will derive a finite element method for this

Figure 3.1: A domain Ω and boundary ∂Ω, ∂Ω = Γ

_D

∪ Γ

_N

.

model.

(20)

3.3 Finite Element Method

3.3.1 Heat equation

In this section, we consider the stationary heat equation (3.5) – (3.6). We assume for now that the right hand side of (3.5) is a given function. To distinguish between u and φ in the left hand sides of (3.5) and (3.7), we denote it by σ(u

^∗

) · |∇φ

^∗

|

²

.

Let us introduce a space

V = {v : kvk + k∇vk < ∞}. (3.10) Multiplying (3.5) by a test function v ∈ V , and integrating on both sides, then using Green’s formula, we get

− Z

Ω

(∇ · k∇u) · vdx = Z

Ω

k∇u · ∇vdx − Z

∂Ω

n · k∇u · vds

= Z

Ω

k∇u · ∇vdx + Z

∂Ω

κ

1

(u − g

1

) · vds

= Z

Ω

σ(u

^∗

) · |∇φ

^∗

|

²

· vdx

since n·k∇u = κ

₁

(g

₁

−u) on ∂Ω. Thus we obtain the variational formulation for the problem (3.5) – (3.6): Find u ∈ V such that

Z

Ω

k∇u · ∇vdx + Z

∂Ω

κ

1

u · vds (3.11)

= Z

Ω

σ(u

^∗

) · |∇φ

^∗

|

²

· vdx + Z

∂Ω

κ

1

g

1

· vds, ∀v ∈ V.

Based on this form, we define a finite element method: Find U ∈ V

_h

⊂ V , such that

Z

Ω

k∇U · ∇vdx + Z

∂Ω

κ

1

U · vds (3.12)

= Z

Ω

σ(u

^∗

) · |∇φ

^∗

|

²

· vdx + Z

∂Ω

κ

1

g

1

· vds, ∀v ∈ V

_h

From the variational formulation and the finite element method, we get

the following theorem:

(21)

Theorem 3.2. For the problem (3.5) – (3.6) with a given exact right hand side of (3.5) which has been denoted as σ(u

^∗

)·|∇φ

^∗

|

²

, the following Galerkin orthogonality property holds

Z

Ω

k∇e · ∇vdx + Z

∂Ω

κ

1

evds = 0, ∀v ∈ V

_h

(3.13) where e = u−U is the error between the exact and the finite element solution.

Proof. Subtracting two equations (3.11) and (3.12), then using the fact that V

h

⊂ V immediately proves this claim.

3.3.2 Potential equation

Let us seek a solution to the potential problem which consisting of (3.7) – (3.9). Assume σ(u) is a given function which is the same as in the right hand side of heat equation (3.5), hence we can rewrite it as σ(u

^∗

).

Let us introduce a space

V

_g₂_,D

= {v : kvk + k∇vk < ∞, v|

_Γ_D

= g

₂

}. (3.14) Multiplying (3.7) by a test function v ∈ V

0,D

, and integrating on both sides, then using Green’s formula, we get

− Z

Ω

(∇ · σ(u

^∗

)∇φ) · vdx = Z

Ω

σ(u

^∗

)∇φ · ∇vdx − Z

∂Ω

n · σ(u

^∗

)∇φ · vds

= Z

Ω

σ(u

^∗

)∇φ · ∇vdx + Z

ΓN

κ

₂

(φ − g

₃

) · vds

= Z

Ω

f · vdx

since n · σ(u

^∗

)∇φ = κ

₂

(g

₃

− φ) on Γ

_N

and v = 0 on Γ

_D

. Thus the variational problem reads as follows: Find φ ∈ V

g2,D

such that

Z

Ω

σ(u

^∗

)∇φ · ∇vdx+

Z

ΓN

κ

2

φ · vds = Z

Ω

f · vdx+

Z

ΓN

κ

2

g

3

· vds, ∀v ∈ V

_0,D

.

(3.15)

We assume g

₂

to be piecewise polynomial and continuous on the boundary

Γ

D

. That means, there is a function Φ

g2

∈ V

_h,D

such that Φ

g2

= g

2

on Γ

D

.

(22)

We introduce the affine subspace V

h,g2,D

= {v ∈ V

h

: v|

ΓD

= g

2

}, then the finite element method reads: Find Φ ∈ V

h,g2,D

, such that

Z

Ω

σ(u

^∗

)∇Φ · ∇vdx+

Z

ΓN

κ

2

Φ · vds = Z

Ω

f · vdx+

Z

ΓN

κ

2

g

3

· vds, ∀v ∈ V

_h,0,D

. (3.16) From the variational formulation and the finite element method, we obtain the following theorem:

Theorem 3.3. For the problem defined by (3.7) – (3.9) with a given function σ(u

^∗

), the following Galerkin orthogonality property holds

Z

Ω

σ(u

^∗

)∇e · ∇vdx + Z

ΓN

κ

2

evds = 0, ∀v ∈ V

_h,0,D

(3.17) where e = φ−Φ is the error between the exact and the finite element solution.

Proof. Subtracting two equations (3.15) and (3.16), then using the fact that V

_h,0,D

⊂ V

_0,D

immediately proves this claim.

To derive an equation for Φ, we will use a technique presented in [1]. We write Φ in the form

Φ = Φ

0

+ Φ

g2

(3.18)

where Φ

_g₂

is any fixed function in V

_h,g₂_,D

and Φ

₀

= 0 on Γ

_D

and thus Φ

₀

∈ V

_h,0,D

. This construction of Φ will satisfy the boundary conditions because of Φ

_g₂

= g

₂

on Γ

_D

. Since Φ

_g₂

is known it remains to determine Φ

₀

, we get the equation: Find Φ

₀

∈ V

_h,0,D

, such that

Z

Ω

σ(u

^∗

)∇Φ

0

· ∇vdx + Z

ΓN

κ

2

Φ

0

· vds

= Z

Ω

f · vdx + Z

ΓN

κ

2

g

3

· vds − Z

Ω

σ(u

^∗

)∇Φ

g2

· ∇vdx − Z

ΓN

κ

2

Φ

g2

· vds

(3.19)

for ∀v ∈ V

_h,0,D

. This is a problem of the same kind as above but with a

modified right hand side. One can prove that the solution Φ = Φ

₀

+ Φ

_g₂

is

independent of the particular choice of the function Φ

_g₂

. In practice Φ

_g₂

is

often chosen to be zero at all interior nodes plus all the nodes on Γ

_N

. And

more details about this technique will be illustrated in the later section.

(23)

Next, let us derive linear systems resulting from these two finite element problems. However, we need to introduce quadrature rules which will be utilized in this paper.

3.4 Gauss Quadrature Rule

In numerical analysis, a quadrature rule is an approximation of the definite integral of a function, usually has the form of a weighted sum of function values at specified points within the domain of integration. Assume f is a given function, a general quadrature rule on a triangle K or a edge E takes the form

Z

K

f dx ≈ X

j

w

_j

f (q

_j

) or

Z

E

f ds ≈ X

j

w

j

f (q

j

)

where q

j

is the set of quadrature points within K or on E.

3.4.1 Area coordinates

To explain the interpolation functions with higher degree, we will introduce the definition so-call area coordinates in [5].

Def inition 3.4 (Area coordinates). For triangular elements, it is possible to construct three non-dimensional coordinates L

_i

(i = 1, 2, 3), which vary in a direction normal to the sides directly opposite each node. The coordinates are defined such that

L

i

= A

i

A i = 1, 2, 3 (3.20)

A =

3

X

i=1

A

i

(3.21)

where A

_i

is the area of the triangle formed by nodes j and k (j, k = 1, 2, 3)

and arbitrary point P in the element, and A is the total area of the element.

(24)

Figure 3.2: The example of area coordinates

For example, assume A

1

is the area of the triangle which is formed by nodes N

2

and N

3

and point P , as showed in Figure 3.2. The point P is at a distance of s from the side connecting nodes N

2

and N

3

. We have

A

1

= b · s 2 A = b · h

2 where h is the distance from the node 1 to the side connecting nodes N

₂

and N

₃

and b is the length of this side. Hence,

L

₁

= A

₁

A = s

h .

Clearly, L

₁

is zero on side N

₂

− N

₃

(hence, zero at nodes N

₂

and N

₃

) and

has a value of unity at node N

₁

. In other words, L

₁

is the finite element

basis function associated with node N

₁

. Similarly, L

₂

and L

₃

are the basis

functions associated with nodes N

2

and N

3

.

(25)

3.4.2 Gauss quadrature rule on triangles

On an arbitrary triangle, using the Gauss quadrature rule, we can obtain the following approximation

Z

K

f (x)dx ≈ |K| ·

N

X

i=1

w

_i

f (G

_i

)

where N is the number of the Gaussian points; G

i

is the i

^th

Gaussian point;

w

i

is the weight of i

^th

Gaussian point; |K| is the area of the triangle. Assume this triangle is defined by three nodes (x

¹₁

, x

¹₂

), (x

²₁

, x

²₂

) and (x

³₁

, x

³₂

), then G

i

can be calculated by the formula

G

_i

=





x

¹₁

x

¹₂

x

¹₃

x

²₁

x

²₂

x

²₃



 ·





 L

ⁱ₁

L

ⁱ₂

L

ⁱ₃







where L

ⁱ₁

, L

ⁱ₂

and L

ⁱ₃

are the area coordinates of point G

_i

in the reference triangle. They can be can be found in Table A.1 in Appendix.

3.4.3 Gauss quadrature rule on edges

On the boundary, it is possible to construct two non-dimensional coordinates L

1

and L

2

. If the edge E lies between two boundary nodes (x

¹₁

, x

¹₂

) and (x

²₁

, x

²₂

), then we use the Gauss quadrature rule to obtain the following approximation

Z

E

f (x)dx ≈ |h

_E

| ·

N

X

i=1

w

_i

f (G

_i

)

where N is the number of the Gaussian points; G

_i

is the i

^th

Gaussian point;

w

_i

is the weight of i

^th

Gaussian point; |h

_E

| is the length of the edge. Further, G

_i

is defined by

G

i

=





x

¹₁

x

¹₂

x

²₁

x

²₂



 ·



 L

ⁱ₁

L

ⁱ₂





where L

ⁱ₁

and L

ⁱ₂

are the area coordinates of point G

_i

on the reference edge.

They can be given in Table A.2 in Appendix.

(26)

3.5 Implementation Details

3.5.1 Derivation of the linear system for the heat equation

Let {ϕ

_i

}

^N_i=1

be basis functions of V

_h

which are defined on the mesh K. Then we write U as a linear combination of the basis functions

U =

N

X

j=1

ζ

j

ϕ

j

(3.22)

with unknown coefficients ζ

j

where j = 1, 2, . . . , N . Inserting (3.22) into (3.12), we get

Z

Ω

k(∇

N

X

j=1

ζ

j

ϕ

j

) · ∇ϕ

i

dx + Z

∂Ω

κ

1

(

N

X

j=1

ζ

j

ϕ

j

)ϕ

i

ds (3.23)

= Z

Ω

σ(u

^∗

) · |∇φ

^∗

|

²

ϕ

i

dx + Z

∂Ω

κ

1

g

1

ϕ

i

ds for i = 1, 2, . . . , N . Introducing the notations

a

¹_i,j

= Z

Ω

k∇ϕ

_i

· ∇ϕ

_j

dx, (3.24)

m

¹_i,j

= Z

∂Ω

κ

₁

ϕ

_i

ϕ

_j

ds, (3.25)

b

¹_i

= Z

Ω

σ(u

^∗

) · |∇φ

^∗

|

²

ϕ

i

dx, (3.26) r

¹_i

=

Z

∂Ω

κ

1

g

1

ϕ

i

ds, (3.27)

for i, j = 1, 2, . . . , N . Then we get

N

X

j=1

(a

¹_i,j

+ m

¹_i,j

)ζ

_j

= b

¹_i

+ r

_i¹

(3.28)

for i = 1, 2, . . . , N , which is a linear system for the coefficients ζ

j

. In the matrix form we write as

(A

¹

+ M

¹

)ζ = b

¹

+ r

¹

(3.29)

with N × N matrices A

¹

and M

¹

, N × 1 vectors b

¹

, r

¹

are defined above

respectively.

(27)

3.5.2 Derivation of the linear system for the potential equation

Let {ϕ

_i

}

^N_i=1

be basis functions of V

_h

which are defined on the mesh K as well. Then we write Φ as a linear combination of the basis functions

Φ =

N

X

j=1

ξ

_j

ϕ

_j

(3.30)

with unknown coefficients ξ

_j

where j = 1, 2, . . . , N . Inserting (3.30) into equation (3.16), we get

Z

Ω

σ(u

^∗

)(∇

N

X

j=1

ξ

j

ϕ

j

) · ∇ϕ

i

dx + Z

ΓN

κ

2

(

N

X

j=1

ξ

j

ϕ

j

)ϕ

i

ds (3.31)

= Z

Ω

f ϕ

i

dx + Z

ΓN

κg

3

ϕ

i

ds for i = 1, 2, . . . , N . Introducing the notations

a

²_i,j

= Z

Ω

σ(u

^∗

)∇ϕ

_i

· ∇ϕ

_j

dx, (3.32) m

²_i,j

=

Z

ΓN

κ

₂

ϕ

_i

ϕ

_j

ds, (3.33)

b

²_i

= Z

Ω

f ϕ

_i

dx, (3.34)

r

_i²

= Z

Γ_N

κ

2

g

3

ϕ

i

ds, (3.35)

for i, j = 1, 2, . . . , N . Then we get

N

X

j=1

(a

²_i,j

+ m

²_i,j

)ξ

_j

= b

²_i

+ r

²_i

(3.36)

for i = 1, 2, . . . , N , which is a linear systems for the coefficients ξ

_j

. In the matrix form we write as

(A

²

+ M

²

)ξ = b

²

+ r

²

(3.37) with N × N matrixes A

²

and M

²

, N × 1 vectors b

²

and r

²

are defined above respectively.

Based on (3.18) and (3.19), let us assume that the first S nodes are

interior nodes including the nodes on Γ

N

, while the last N − S nodes are

(28)

boundary nodes on Γ

D

for the linear system of the potential equation (3.37).

These boundary nodes are fixed since the nodal values of Φ should be g

2

. Using these nodes numbering, we can partition the linear system (A

²

+ M

²

)ξ = b

²

+ r

²

into the following form





A

²_0,0

+ M

_0,0²

A

²_0,g₂

+ M

_0,g² ₂

A

²_g₂_,0

+ M

_g²₂_,0

A

²_g₂_,g₂

+ M

_g²₂_,g₂



 ·



 ξ

₀

ξ

_g₂



 =





b

²₀

+ r

²₀

b

²_g₂

+ r

²_g₂





where A

²_0,0

and M

_0,0²

are the upper left S × S block of A

²

and M

²

, while A

²_g₂_,g₂

and M

_g²₂_,g₂

are the lower right (N − S) × (N − S) block of A

²

and M

²

. Rearranging the first S equations of this linear system, we have the S × S linear system

(A

²_0,0

+ M

_0,0²

)ξ

0

= (b

²₀

+ r

²₀

) − (A

²_0,g₂

+ M

_0,g² ₂

)ξ

g2

(3.38) from which the unknown interior nodal values as well as the nodal values on Γ

_N

of Φ can be determined.

3.5.3 Assembling of the stiffness matrices

Recall the notations (3.24) and (3.32), the local 3 × 3 stiffness matrices are given by

A

^1,K_i,j

= Z

K

k∇ϕ

_i

∇ϕ

_j

dx, A

^2,K_i,j

=

Z

K

σ(u

^∗

)∇ϕ

_i

∇ϕ

_j

dx, for i, j = 1, 2, 3.

Consider a triangle K with the nodes (x

¹₁

, x

¹₂

), (x

²₁

, x

²₂

) and (x

³₁

, x

³₂

). To each node N

i

(i = 1, 2, 3), there is a hat function ϕ

i

associated, which takes on the value 1 at node N

_i

and 0 at other nodes. Each hat function is a linear function on K so it takes the form

ϕ

i

= a

i

+ b

i

x

i

+ c

i

x

2

, i = 1, 2, 3

where the coefficients a

i

, b

i

and c

i

are determined by the following linear

(29)

systems







1 x

¹₁

x

¹₂

1 x

²₁

x

²₂

1 x

³₁

x

³₂







· 



 a

₁

b

₁

c

₁







=





 1 0 0







= e

1

,







1 x

¹₁

x

¹₂

1 x

²₁

x

²₂

1 x

³₁

x

³₂







· 



 a

₂

b

₂

c

₂







=





 0 1 0







= e

2

,







1 x

¹₁

x

¹₂

1 x

²₁

x

²₂

1 x

³₁

x

³₂







· 



 a

₃

b

₃

c

₃







=





 0 0 1







= e

3

.

The gradient of ϕ

i

is just the constant vector ∇ϕ

i

= [b

i

c

i

]. Then using the Gauss quadrature rule, we get

A

^1,K_i,j

= Z

K

k∇ϕ

i

∇ϕ

_j

dx

= (b

_i

b

_j

+ c

_i

c

_j

) Z

K

kdx

≈ (b

_i

b

j

+ c

i

c

j

) · |K| ·

Q

X

l=1

(w

_l

· k(G

_l

)),

A

^2,K_i,j

= Z

K

σ(u

^∗

)∇ϕ

i

∇ϕ

_j

dx

≈ (b

_i

b

j

+ c

i

c

j

) · |K| ·

Q

X

l=1

(w

_l

· σ(u

^∗_l

))

where i, j = 1, 2, 3; G

_l

is the l

^th

Gaussian point; w

_l

is the weight of l

^th

Gaussian point; L

^l_i

is the i

^th

area coordinate of the l

^th

Gaussian point defined

by (3.20); σ(u

^∗_l

) is the value of σ(u

^∗

) on the l

^th

Gauss point; Q is the number

of Gaussian points.

(30)

3.5.4 Assembling of the load vectors

Recall the notations (3.26) and (3.34), on each element K, we get local 3 × 1 element vectors b

¹_K

and b

²_K

with entries

b

^1,K_i

= Z

K

σ(u

^∗

) · |∇φ

^∗

|

²

ϕ

i

dx, b

^2,K_i

=

Z

K

f ϕ

i

dx

for i = 1, 2, 3. Using the Gauss quadrature rule to compute these integrals, we obtain

b

^1,K_i

≈ |K| ·

Q

X

l=1

(w

_l

· σ(u

^∗_l

) · |∇φ

^∗_l

|

²

· L

^l_i

)

b

^2,K_i

≈ |K| ·

Q

X

l=1

(w

_l

· f (G

_l

) · L

^l_i

),

where i, j = 1, 2, 3; G

_l

is the l

^th

Gaussian point; w

_l

is the weight of l

^th

Gaussian point; L

^l_i

is the i

^th

area coordinate of the l

^th

Gaussian point defined by (3.20); σ(u

^∗_l

) is the value of σ(u

^∗

) on the l

^th

Gauss point; |∇φ

^∗_l

|

²

is the value of |∇φ

^∗

|

²

on the l

^th

Gauss point; Q is the number of the Gaussian points.

3.5.5 Assembling of the boundary contributions

Two nodes of a triangle K lie along the domain ∂Ω, then the edge between them will contribute to matrices entries M

_i,j¹

, M

_i,j²

and vectors entries r

¹_i

, r

²_i

M

_i,j^1,E

= Z

E

κ

1

ϕ

i

ϕ

j

ds ≈ κ

1

|h

_E

| ·

Q

X

l=1

(w

l

· L

^l_i

· L

^l_j

)

r

^1,E_i

= Z

E

κ

1

g

1

ϕ

i

ds ≈ κ

1

|h

_E

| ·

Q

X

l=1

(w

l

· g

₁

(G

l

) · L

^l_i

)

M

_i,j^2,E

= Z

E

κ

₂

ϕ

_i

ϕ

_j

ds ≈ κ

₂

|h

_E

| ·

Q

X

l=1

(w

_l

· L

^l_i

· L

^l_j

)

r

^2,E_i

= Z

E

κ

₂

g

₃

ϕ

_i

ds ≈ κ

₂

|h

_E

| ·

Q

X

l=1

(w

_l

· g

₃

(G

_l

) · L

^l_i

)

(31)

where i, j = 1, 2, 3; G

l

is the l

^th

Gaussian point; w

l

is the weight of l

^th

Gaussian point; L

^l_i

is the i

^th

area coordinate of the l

^th

Gaussian point defined by (3.20); Q is the number of the Gaussian points.

3.6 Iterative Methods for non-linear problems

So far we have derived two linear systems of two single physical problems which are based on an assumption that σ(u

^∗

) · |∇φ

^∗

|

²

and σ(u

^∗

) are two given functions. Let us write the relationship between u

^∗

, φ

^∗

and U , Φ in the following way

U = f

₁

(u

^∗

, φ

^∗

), (3.39)

Φ = f

2

(u

^∗

). (3.40)

However, our model problem is a multiphysical problem involving both heat and potential equations. Hence u

^∗

and φ

^∗

must be the exact solutions u and φ defined by (3.5) – (3.9). Then as the numerical solutions of u and φ in the model problem, (3.39) and (3.40) can be represented as

U = f

₁

(U, Φ), (3.41)

Φ = f

2

(U ). (3.42)

To solve these linear systems, we need to introduce iterative methods in our model.

An iterative method attempts to solve a problem by finding successive approximations to the solution starting from an initial guess. If an equation can be put into the form F (X) = X (here, X could be a vector and contains several elements), a solution X is an attractive fixed point of the function F , then one may begin with a point X

⁽¹⁾

in the basin of attraction of X.

Let X

⁽ⁱ⁺¹⁾

= F (X

ⁱ

) for i ≥ 1, and the sequence {X

⁽ⁱ⁾

}

_i≥1

will converge to the solution X.

Examples of iterative methods are Jacobi method, Gauss-Seidel (GS)

method and Successive Over-relaxation (SOR) method.

(32)

• Jacobi method: is an iterative technique that solves present values X by using previous values X in the right hand side. This method can be written as X

⁽ⁱ⁺¹⁾

= F (X

⁽ⁱ⁾

);

• GS method: the computation of x

⁽ⁱ⁺¹⁾_j

uses the elements of X

⁽ⁱ⁺¹⁾

that have already been computed (denoted as X

∗⁽ⁱ⁺¹⁾

) and the elements of X

⁽ⁱ⁾

that have to be advanced to iteration i + 1 (denoted as X

∗∗⁽ⁱ⁾

).

This method can be written as X

⁽ⁱ⁺¹⁾

= F (X

∗⁽ⁱ⁺¹⁾

, X

∗∗⁽ⁱ⁾

);

• SOR method: is a variant of the GS method, resulting in faster convergence. It introduces a relaxation factor ω, which is a constant and greater than 0. This method can be written as X

⁽ⁱ⁺¹⁾

= (1 − ω)X

⁽ⁱ⁾

+ ωF (X

∗⁽ⁱ⁺¹⁾

, X

∗∗⁽ⁱ⁾

). When ω = 1, it is the GS method. Since the similar method can be used for any slowly converging iterative process, we can use it to improve Jacobi method as well.

In this paper, assume we start from two guessing values U

⁰

and Φ

⁰

, then these iterative methods would be implemented between heat and potential equations as:

• Jacobi method : {U

⁽ⁱ⁾

, Φ

⁽ⁱ⁾

} = {f

₁

(U

⁽ⁱ⁻¹⁾

, Φ

⁽ⁱ⁻¹⁾

), f

₂

(U

⁽ⁱ⁻¹⁾

)};

• GS method : {U

⁽ⁱ⁾

, Φ

⁽ⁱ⁾

} = {f

₁

(U

⁽ⁱ⁻¹⁾

, Φ

⁽ⁱ⁻¹⁾

), f

2

(U

⁽ⁱ⁾

)};

• SOR method : {U

⁽ⁱ⁾

, Φ

⁽ⁱ⁾

} = {(1 − ω)U

⁽ⁱ⁻¹⁾

+ ωf

₁

(U

⁽ⁱ⁻¹⁾

, Φ

⁽ⁱ⁻¹⁾

), (1 − ω)Φ

⁽ⁱ⁻¹⁾

+ ωf

₂

(U

⁽ⁱ⁾

)}

where i ≥ 1 and all procedures stop when a certain tolerance is reaching.

(33)

Chapter 4

Adaptive Finite Element Method

4.1 A Posteriori Error Estimates

Let us revisit the model problem

−∇ · k∇u = σ(u) · |∇φ|

²

in Ω, (4.1) n · k∇u = κ

1

(g

1

− u) on ∂Ω, (4.2)

−∇ · σ(u)∇φ = f in Ω, (4.3)

φ = g

2

on Γ

D

, (4.4)

n · σ(u)∇φ = κ

₂

(g

₃

− φ) on Γ

_N

. (4.5) For the heat equation defined by (4.1) – (4.2), we have a posteriori estimate:

Theorem 4.1. For the finite element approximation U of the exact solution u to (4.1) and (4.2) with a given right hand side of (4.1), σ(u

^∗

) · |∇φ

^∗

|

²

, the following a posteriori error estimate holds

k∇(u − U )k

²_L2(Ω)

+ ku − U k

²_L2(∂Ω)

≤ C X

K∈K

ρ

²_K

(U, u

^∗

, φ

^∗

) + X

E∈∂K

ρ

²_E

(U )

(4.6)

where C is a constant, the element residual in the interior domain ρ

K

(U, u

^∗

, φ

^∗

)

(34)

is defined by

ρ

_K

(U, u

^∗

, φ

^∗

) = h

_K

kσ(u

^∗

) · |∇φ

^∗

|

²

+ ∇ · k∇U k

_L2(K)

+ 1

2 h

^1/2_K

k[n · k∇U ]k

_L2(∂K\∂Ω)

(4.7) and the element residual on the boundary ρ

_E

(U ) is defined by

ρ

_E

(U ) = h

^1/2_K

kκ

₁

(g

₁

− U ) − n · k∇U k

_L2(E∩∂Ω)

. (4.8) Here [n · k∇U ] denotes the jump in the k times normal derivative of U at an interior edge ∂K

₁

∩ ∂K

₂

, i.e.

[n · k∇U ]|

_∂K₁∩∂K2

= n

1

· k∇U

₁

+ n

2

· k∇U

₂

(4.9) with U

i

= U |

Ki

and n

i

is the exterior unit normal of K

i

.

Proof. Let e = u − U be the error. If k, κ

₁

≥ α > 0, we have α · k∇ek

²_L2(Ω)

+ kek

²_L2(∂Ω)

≤ Z

Ω

k∇e · ∇edx + Z

∂Ω

κ

₁

e · eds

= Z

Ω

k∇e · ∇(e − πe)dx + Z

∂Ω

κ

₁

e · (e − πe)ds

where we have used the Galerkin orthogonality (see Theorem 3.2) to subtract the interpolant πe. Splitting this into a sum over the elements and using the Green’s formula, or integration by parts, we further have

α · k∇ek

²_L2(Ω)

+ kek

²_L2(∂Ω)

≤ X

K∈K

Z

K

k∇e · ∇(e − πe)dx + X

E∈∂K

Z

E

κ

₁

e · (e − πe)ds

= X

K∈K

(−

Z

K

∇ · k∇e · (e − πe)dx + Z

∂K

n · k∇e · (e − πe)ds)

+ X

E∈∂K

Z

E

κ

₁

e · (e − πe)ds.

First of all, let us consider the element residuals in the interior domain.

Finite Element Methods for Microelectromechanical Systems