‰ Elliptic Equations

7

P ARTIAL D IFFERENTIAL E QUATIONS

‰ Introduction

‰ Elliptic Equations

‰ Parabolic Equations

‰ Hyperbolic Equations

‰ Finite Element Method — An Introduction

‰ Boundary Element Method — An Introduction

7.1 Introduction

• Partial Differential Equations (PDE) are differential equations which have at least two independent variables, e.g.

3

4 ₂

2 2

2 + =

∂ + ∂

∂

∂ u

y xy u x

u second order & linear

y x x

u x

u =

∂

∂ + ∂

⎟⎟⎠

⎜⎜ ⎞

⎝

⎛

∂

∂

2 3 3

2 2

6 third order & nonlinear

• For a second order linear two-dimensional equation, a general equation is 0

, , ,

2 ,

2 2

2

2 ⎟⎟=

⎠

⎜⎜ ⎞

⎝

⎛

∂

∂

∂ + ∂

∂ + ∂

∂

∂ + ∂

∂

∂

y u x u u y x y D

C u y x B u x

A u (7.1)

which can be divided into three types.

TABLE 7.1 Types of second order linear PDEs AC

B² −4 Type Examples

< 0 Elliptic ^{2 0}

2 2 2

∂ = +∂

∂

∂

y T x

T

Laplace equation

= 0 Parabolic ²

2

x k T t T

∂

= ∂

∂

∂

Heat conduction equation

> 0 Hyperbolic ²

2 2 2

2 1

x y c t

y

∂

= ∂

∂

∂

Wave equation

7.2 Elliptic Equations

• Elliptic PDEs are generally related to steady-state problems with diffusivity having boundary conditions, e.g. the Laplace equation.

• Consider a steady-state 2-D heat conduction problem:

Δz

Δx

Δy q(x)

q(y)

q(x + Δx) q(y + Δy)

x y

FIGURE 7.1 Thin plate having a thickness Δz with temperatures at the boundary

Taking q(x) as the heat flux [J/m²⋅s], the heat flow over the element can be written in a time interval of Δt as:

( )

( )

(

)

(

)

( ) ( )

( ) ( ( ) ( ) )

( ) ( ) ( ) ( )

0 0 Δ =

Δ + + −

Δ Δ +

−

= Δ Δ +

− +

Δ Δ +

−

y y y q y q x

x x q x q

x y y q y q y x x q x q

For Δx,Δy→0:

=0

∂

−∂

∂

−∂

y q x

qx y

(7.2)

• In conduction heat transfer, the relationship between heat flux and temperature is given by the Fourier Law:

y C T k x q

C T k

q_{x y}

∂

− ∂

∂ =

− ∂

= ρ , ρ (7.3)

where k is heat diffusivity [m²/s], ρ is density [kg/m³] and C is heat capasity [J/kg°C]. Combining Eq. (7.2) with Eq. (7.3) gives the Laplace equation:

2 0

2 2

2 =

∂ +∂

∂

∂

y T x

T (7.4)

• If a heat source or heat generation of Q(x,y) is present in the domain, a Poisson equation can be formed:

( )

2 2 2

2 + =

∂ +∂

∂

∂ Q x y

y T x

T (7.5)

• To solve Eq. (7.4), use the central differencing:

( )

( )

1 , 2 2

2 , 1 ,

, 1 2

2

2 2

y T T T

y T

x T T T

x T

j i j i j

i

j i j i j i

Δ +

= −

∂

∂

Δ +

= −

∂

∂

− +

− +

Hence, Eq. (7.4) can be written in an algebraic form:

( ) ( )

2 2

2

1 , , 1 , 2

, 1 ,

,

1 =

Δ + + −

Δ +

− _{− + −}

+

y T T T

x T T

T_{i j i j i j i j i j i j}

Taking Δx = Δy gives

0 4 _,

1 , 1 , , 1 ,

1 + ₋ + ₊ + ₋ − =

+ j i j i j i j i j

i T T T T

T (7.6)

0,0 x y

m+1,0 0,n+1

i,j

i,j+1 i−1,j

i,j−1 i+1,j

FIGURE 7.2 Finite difference grid for solving the Laplace equation

• Eq. (7.6) needs boundary conditions (BC), which may be in the form of:

1. Fixed value — Dirichlet (see Fig. 7.3), e.g. Fixed temperature

2. First derivative, or gradient —Neumann, e.g. Fix heat flux or insulated

75°C

T11 T21 T31

T12 T22 T32

T13 T23 T33

100°C

50°C

0°C

FIGURE 7.3 Grid/Node representatioan and BC for Fig. 7.1

For node (1,1) — T01 = 75°C and T10 = 0°C:

0 4 ₁₁

10 12 01

21+T +T +T − T =

T

75 4T₁₁−T₁₂ −T₂₁ =

Hence, the complete system (can be assembled into a matrix equation) is:

150 4

100 4

175 4

50 4

0 4

75 4

50 4

0 4

75 4

33 23

32

33 23

13 22

23 13

12

33 32

22 31

23 32

22 12

21

13 22

12 11

32 31

21

22 31

21 11

12 21

11

= +

−

−

=

− +

−

−

=

− +

−

=

− +

−

−

=

−

− +

−

−

=

−

− +

−

=

− +

−

=

−

− +

−

=

−

−

T T

T

T T

T T

T T

T

T T

T T

T T

T T

T

T T

T T

T T

T

T T

T T

T T

T

• If the Gauss-Seidel method is used, at each node (i,j):

4

1 , 1 , , 1 ,

1 ,

− +

−

+ + + +

= ^{i j i j i j i j}

j i

T T

T

T T (7.7)

which can be solve iteratively via relaxation approach:

( )

new , new

,_{j i j} 1 _{i j}

i T T

T =λ + −λ (7.8)

where λ is the relaxation parameter 1≤ λ ≤2, and is subjected to the termination criterion as followed:

%

new 100

, old , new ,

, − ×

=

j i

j i j i

a T

T T

j

ε i

Example 7.1

Solve the problem in Fig. 7.3 using the Gauss-Seidel method using the termination criterion ε_a ≤1% dan the relaxation parameter λ = 1.5.

Solution

Let all the initial values be 0°C. For the first iteration:

75 . 4 18

0 0 75 0 4

10 12 01 21

11 + + + =

+ = +

=T +T T T T

(

) (

)

5 .

new 1

11 = + − =

T

03125 .

4 7

0 0 125 . 28 0 4

20 22 11 31

21 + + + =

+ = +

=T +T T T T

(

) (

)

.

new 1

21 = + − =

T

13672 .

4 15

0 0 54688 .

10 50 4

30 32 21 41

31 + + + = + + + =

=T T T T

T

(

) (

)

5 .

new 1

31 = + − =

T

and, for other nodes:

99554 .

96

18579 .

34 46900

. 74

45703 .

18 12696

. 80

67188 .

38

33 23 23

22 13

12

=

=

=

=

=

=

T T T

T T

T

Error for node (1,1) — for the first iteration, all errors are 100%:

% 100 125 100

. 28

0 125 . 28

11 = − × =

εa

For the second iteration:

68736 .

67

86833 .

71

60108 .

28 95872

. 87

63333 .

61

35718 .

22 21973

. 75

95288 .

57

51953 .

32

33 32 31

23 22 21

13 12 11

=

=

=

=

=

=

=

=

=

T T T

T T T

T T T

The process is repeated until the ninth iteration in which the termination criterion is fulfilled (εa < 1%):

71050 .

69

33999 .

52

88506 .

33 06402

. 76

11238 .

56

29755 .

33 58718

. 78

21152 .

63

00061 .

43

33 32 31

23 22 21

13 12 11

=

=

=

=

=

=

=

=

=

T T T

T T T

T T T

`

7.3 Parabolic Equations

• Parabolic PDEs are generally related to transient problems with diffusivity, e.g. the 1-D heat conduction equation.

• For a transient problem, there are three approaches:

1. Explicit method, 2. Implicit method,

3. Semi-implicit method — the Crank-Nicolson method.

• Consider a transient 1-D heat conduction problem:

input − output = storage

( )

(

)

x

Hot Cold

FIGURE 7.4 A rod with different temperature at its ends

Dividing with the volume ΔxΔyΔz and the time interval Δt:

( ) ( )

t C T x

x x q x q

Δ

= Δ Δ

Δ +

− ρ

In the limits of Δx,Δt →0:

t C T x q

∂

= ∂

∂

−∂ ρ _(7.9)

By using the Fourier law, Eq. (7.3), Eq. (7.9) becomes

t T x

k T

∂

= ∂

∂

∂

2 2

(7.10)

0,0 x t

m+1,0 i,l

i,l+1 i−1,l

i,l−1 i+1,l

FIGURE 7.5 Finite difference grid for the heat conduction equation

• For the explicit method, the right-hand side of Eq. (7.9) can be discretised via central difference as:

( )

1 2

2 2

x T T T

x

T _i^l _i^l _i^l Δ

+

= −

∂

∂ _{+ −}

and, the left-hand side can be discretised via forward difference as:

t T T

t

T _i^l _i^l Δ

= −

∂

∂ ⁺¹

Hence, Eq. (7.10) can be written in the algebraic form:

( )

kT

l i l i l i l i l

i

Δ

= − Δ

+

− ₋ ⁺

+ 1

2 1

1 2

(7.11)

(

)

l i l

i l

i l

i T T T T

T ⁺¹ = +λ ₊₁ −2 + _{−1 (7.12)}

where λ =kΔt (Δx)² .

x_i

Time differencing grid Space differencing grid

xi−1 xi+1

t^l t^l+1

FIGURE 7.6 Computational grid for the explicit method

Example 7.2

Use the explicit method to determine the temperature distribution for a slender rod having a length of 10 cm. At time t= 0, the temperature of the rod is 20°C and the boundary conditions are T(0) = 100°C and T(10 cm) = 50°C. Use the conduction coefficient k = 0.835 cm²/s, the time interval Δt= 0.5 s and the step size Δx = 2 cm.

Solution Calculate λ:

( ) ( )( )

104375 .

2 0 5 . 0 835 . 0

2

2 = =

Δ

= Δ x

t λ k

At t= 0:

0 20

4 0 3 0 2 0

1 =T =T =T =

T

Use Eq. (7.12). At t= 0.5 s:

[

]

1044 . 0

1 20

1 = + − + =

T₂¹ = 20+0.1044

[

]

[

]

[

]

1044 . 0

1 20

4 = + − + =

T

At t= 1 s:

[

]

1044 . 0 352 .

2 28

1 = + − + =

T₂² =20+0.1044

[

]

[

]

[

]

1044 . 0 132 .

2 23

4 = + − + =

T

The calculation can be continued to produce the result as in Fig. 7.7.

`

0 20 40 60 80 100

0 2 4 6 8 10 x (cm)

T (°C)

t = 0

t = 15 s t = 10 s

t = 6 s t = 3 s t = 20 s

t → ∞

FIGURE 7.7 Result for Ex. 7.2

• For the explicit method, a more accurate result can be obtained if Δx and Δt approach zeros. However, stability of the results is only obtained if:

( )

k t x

2 2 1 2 1

or Δ

≤ Δ

λ ≤

If the stability condition is not fulfilled, the results are contaminated by oscillation.

• To prevent oscillation, the implicit method can be used, where the second order spatial derivative is approximated at time l+1:

( )

1 1 1 1

1 2

2 2

x T T

T x

T _i^l _i^l _i^l Δ

+

= −

∂

∂ ₊^{+ +} ₋⁺

x_i

Time differencing grid Space differencing grid

x_i−1 x_i+1

t^l t^l+1

FIGURE 7.8 Computational grid for the implicit method

Hence, Eq. (7.10) becomes

( )

T T

k T

l i l i l

i l i l

i

Δ

= − Δ

+

− ⁺ ₋^{+ +}

+

+ 1

2

1 1 1 1

1 2

(7.13)

( )

l

i T T T

T + + − =

−λ ₋⁺₁¹ 1 2λ ⁺¹ λ ₊⁺₁¹ _(7.14)

Eq. (7.14) leads to n simultaneous linear equations having n unknowns.

Example 7.3

Repeat Ex. 7.2 using the implicit method.

Solution

From Ex. 7.2, λ = 0.104375. From Eq. (7.14):

( )

[ ]

( )

[ ]

M

44 . 30 1044

. 0 1044

. 0 2 1

20 1044

. 0 1044

. 0 2 1 ) 100 ( 1044 . 0

1 2 1

1

1 2 1

1

=

− +

=

− +

+

−

T T

T T

Hence, the following system of linear equations can be formed:

⎪⎪

⎭

⎪⎪

⎬

⎫

⎪⎪

⎩

⎪⎪

⎨

⎧

=

⎪⎪

⎭

⎪⎪

⎬

⎫

⎪⎪

⎩

⎪⎪

⎨

⎧

⎥⎥

⎥⎥

⎦

⎤

⎢⎢

⎢⎢

⎣

⎡

−

−

−

−

−

−

22 . 25 20 20

44 . 30

2088 . 1 1044 . 0 0

0

1044 . 0 2088 . 1 1044 . 0 0

0 1044 . 0 2088 . 1 1044 . 0

0 0

1044 . 0 2088 . 1

1 4

1 3

1 2

1 1

T T T T

[

]

= T

`

• The combination of the explicit and implicit approaches produces the semi- implicit method, and one of its kind is the Crank-Nicolson method:

x_i

Time differencing grid Space differencing grid

x_i−1 x_i+1

t^l t^l+1

FIGURE 7.9 Computational grid for the Crank-Nicolson method

In this method, the finite difference term for the spatial derivative is

( ) ( )

⎢⎣

⎡

Δ + + −

Δ +

= −

∂

∂ _{+ − +}^{+ +} ₋⁺

2

1 1 1 1

1 2

1 1

2

2 2 2

2 1

x T T

T x

T T T

x

T _i^l _i^l _i^l _i^l _i^l _i^l

(7.15) Hence, Eq. (7.10) becomes

( )

( )

l

i T T T T T

T₋⁺₁¹ +21+ ⁺¹ − ₊⁺₁¹ = ₋₁ +21− + ₊₁

−λ λ λ λ λ λ _(7.16)

Example 7.4

Repeat Ex. 7.2 using the Crank-Nicolsonmethod.

Solution

From Ex. 7.2, λ = 0.104375. From Eq. (7.16):

7866 . 58 10437

. 0 20874

. 2

) 20 ( 10437 .

0 20 ) 10437 .

0 1 ( 2 ) 100 ( 10437 .

0

10437 .

0 )

10437 .

0 1 ( 2 ) 100 ( 10437 .

0

1 2 1

1

1 2 1

1

=

−

+

− +

=

− +

+

−

T T

T T

By considering other nodes, the following system of linear equations can be formed:

⎪⎪

⎭

⎪⎪

⎬

⎫

⎪⎪

⎩

⎪⎪

⎨

⎧

=

⎪⎪

⎭

⎪⎪

⎬

⎫

⎪⎪

⎩

⎪⎪

⎨

⎧

⎥⎥

⎥⎥

⎦

⎤

⎢⎢

⎢⎢

⎣

⎡

−

−

−

−

−

−

3496 . 48 40 40

7866 . 58

2087 . 2 1044 . 0 0

0

1044 . 0 2087 . 2 1044 . 0 0

0 1044 . 0 2087 . 2 1044 . 0

0 0

1044 . 0 2087 . 2

1 4

1 3

1 2

1 1

T T T T

[

]

= T

`

7.4 Hyperbolic Equations

• Hyperbolic PDEs are generally related to transient problems with convection, e.g. the 1-D wave equation.

• Consider a 1-D wave equation, which is a hyperbolic PDE:

=0

∂ + ∂

∂

∂

x a u t

u (7.17)

• One of the methods is the MacCormack’s technique, which is an explicit finite-difference technique and is second-order-accurate in both space and time. By using the Taylor series:

t t u u

u_i^{t t} _i^t ⎟ Δ

⎠

⎜ ⎞

⎝

⎛

∂ + ∂

Δ =

+

avg

(7.18)

• This method consists of two steps: predictor and corrector. In predictor step, use forward difference in the right-hand side:

⎟⎟⎠

⎜⎜ ⎞

⎝

⎛ Δ

− −

⎟ =

⎠

⎜ ⎞

⎝

⎛

∂

∂ ₊

x u a u

t

u ^t _i^t _i^t

i

1 (7.19)

Thus, from the Taylor series, the predicted value of u is:

t t u u

u

t

i t

i t t

i ⎟ Δ

⎠

⎜ ⎞

⎝

⎛

∂ + ∂

Δ =

+ (7.20)

• In corrector step, by replacing the spatial derivatives with rearward differences:

⎟⎟⎠

⎜⎜ ⎞

⎝

⎛

Δ

− −

⎟⎟ =

⎠

⎜⎜ ⎞

⎝

⎛

∂

∂ ^{+Δ +Δ} ₋^+Δ x

u a u

t

u ^{t t} _i^{t t} _i^{t t}

i

1 (7.21)

The average of time derivative can be obtained using

⎥⎥

⎦

⎤

⎢⎢

⎣

⎡

⎟⎟⎠

⎜⎜ ⎞

⎝

⎛

∂ + ∂

⎟⎠

⎜ ⎞

⎝

⎛

∂

= ∂

⎟⎠

⎜ ⎞

⎝

⎛

∂

∂ ^t+Δt

i t

i t

u t

u t

u

2 1

avg

(7.22)

• Hence the final, “corrected” value at time t+Δt is:

t t u u

u_i^{t t} _i^t ⎟ Δ

⎠

⎜ ⎞

⎝

⎛

∂ + ∂

Δ =

+

avg

• The accuracy of the solution for a hyperbolic PDE is dependent on truncation and round off errors, and the term representing it is called artificial viscosity ₂¹ aΔx

(

)

• The effect of artificial viscosity leads to numerical dissipation, which is originated by the even-order derivatives in the truncated term, but it improves stability.

FIGURE 7.10 Numerical dissipation: (a) t = 0, (b) t > 0

• Another opposite effect is known as numerical dispersion, which is originated by the odd-order derivatives in the truncated term and causes

‘wiggles’.

FIGURE 7.11 Numerical dispersion: (a) t = 0, (b) t > 0

7.5 Finite Element Method ⎯ An Introduction

• The Finite Element Method (FEM) is a computer aided mathematical technique used to obtain an approximate numerical solution of a response of a physical system which is subjected to an external loading.

• By using this technique, the computational domain which is theoretically a continuum, is being discretised in form of simple geometries.

• The mesh is the computational domain which is an assembly of discrete elemental blocks known as finite elements, and the vertices defining the elements are called nodes.

• Governing equation is employed at each element to form a set of algebraic equations ⎯local system.

• Local equations assembled to form a global system which is the solved to yield a vector of variables.

Node & element Sample problem Mesh

y

x

z

x,y node

element

3-D 2-D

1-D

FIGURE 7.10 Examples of elements and their applications

Geometry, material properties, mesh,BC, IC

Discretisation of structures into elements

Formation of local stiffness matrix [k]

Formation of global stiffness matrix [K]

Boundary conditions as constraints to the model

Load conditions to the model

Solving the system of Linear equation

[K]{a}= {F}

Calculaton of strains, stresses and forces Boundary

conditions

Loads

Graphical visualisation print-out & files

FIGURE 7.11 Algorithm for a force analysis using FEM

• Consider a 1-D steady state heat conduction:

( )

2

2 + =

∂

∂ Q x

x

k T (7.17)

Eq. (7.17) needs appropriate boundary conditions such that:

) (

0 0

∞

=

=

−

=

=

T T h q

T T

L L

x

x (7.18)

T₀

L

h T∞

x

T0 T_L

1 2 3 4

1 2 3

FIGURE 7.12 Boundary condition for a 1D steady state heat conduction For the first element:

1 2

x1

x

x₂

1 2

ξ = −1 ξ = +1

ξ

FIGURE 7.13 The first element in the finite element method

The transformation of the coordinate system to a local system is given by an isoparametric coordinateξ, i.e.

( )

2

1 1

2

−

− −

= x x

x

ξ x (7.19)

Thus, in order to calculate the temperature at the middle section, a linear interpolation function or a linear shape function can be used:

2 ) 1 (

2 ) 1 (

2 1

ξ ξ ξ ξ

= +

= −

N N

(7.20)

Hence, the temperature can be interpolated using the shape function as followed:

T e

N T N

T(ξ)= _{1 1} + _{2 2} =NT (7.21)

Ni

ξ = -1 ξ = +1

N1 N2

T

ξ = -1 ξ = +1

T1

T2

T = N1T1 + N2T2

FIGURE 7.14 Linear shape function Differentiation of Eq. (7.21) gives

x dx d x

1 2

2

= − ξ Using a chain rule:

[ ]

e

e

x x

dx d d dT dx

dT

BT

T

=

− −

=

=

1 , 1 1

1 2

ξ ξ

where

[

]

1

1 2

− −

= x x B

• In energy form, the heat conduction problem can be represented by

( )

2

2 Ω=

⎭⎬

⎫

⎩⎨

⎧ +

∂

∫

k T T

i

( )

∫

∫

=

∏_T ^L dx ^LTQ x dx h T_L T dx

k dT

0

2 2

1 0

2 2

1 ( ) (7.22)

Discretisation of Eq. (7.22) gives

( )

2 1 1

1 1

1 T 2

1

2 2

T

− ∞

− ⎥⎦⎤ + −

⎢⎣⎡

⎥⎦ −

⎢⎣ ⎤

= ⎡

∏

∑

∫

∑

∫

e

e e e e

e

e e e

T N T

B B

T ξ ξ (7.23)

Hence, the stiffness matrix for the element is

⎥⎦

⎢ ⎤

⎣

⎡

−

= −

=

∫

1 1 2

1 1

T

e e e

e

l d k l

k B B ξ

k (7.24)

and, the heat rate vector is

⎭⎬

⎫

⎩⎨

= ⎧

=

∫

1 2 2

1 1

e e e

e Ql

l d

Q N ξ

r (7.25)

The global stiffness matrix can be produced via an assembly of all stiffness matrices for all elements, i.e.

∑

←

e

ke

K (7.26)

Likewise, the global heat rate vector is an assembly by local heat rate vectors:

∑

←

e

re

R (7.27)

To combine with the boundary condition T1 =T0, a penalty method can be used:

( )

⎪⎪

⎭

⎪⎪

⎬

⎫

⎪⎪

⎩

⎪⎪

⎨

⎧

∞ + +

=

⎪⎪

⎭

⎪⎪

⎬

⎫

⎪⎪

⎩

⎪⎪

⎨

⎧

⎥⎥

⎥⎥

⎦

⎤

⎢⎢

⎢⎢

⎣

⎡

+ +

hT R

R

CT R

T T T

h K K

K

K K

K

K K

C K

L L

LL L

L

L L

M M

L M M

M

L L

2

0 1

2 1

2 1

2 22

12

1 12

11

(7.28)

or, in a matrix form

K⋅T = R (7.29) where the penalty parameter C can be estimated as:

104

max ×

= _ij

C K (7.30)

Example 7.5

Fig. 7.15 shows plate of three composites, where the temperature at the right-hand side is T0 = 20°C and the left-hand side is subjected to convection with T_∞ = 800°C and h = 25 W/m²°C. Obtain a temperature distribution across the plate using the finite element method.

k1 k2 k3

k1 = 20 W/m°C k₂ = 30 W/m°C k₃ = 50 W/m°C

T0 = 20°C

0.3m 0.15m 0.15m

1 2 3 4

1 2 3 T4 = 20°C T1

T∞ = 800°C

RAJAH 7.15 Plate of three composites used in Ex. 7.5 Solution

Divide the domain into three elements and construct local stiffness matrices for all elements:

⎥⎦

⎢ ⎤

⎣

⎡

−

= −

⎥⎦

⎢ ⎤

⎣

⎡

−

= −

⎥⎦

⎢ ⎤

⎣

⎡

−

= −

1 1

1 1 15 . 0

50

1 1

1 1 15 . 0

30

1 1

1 1 3 . 0

20

3 2 1

k k k

Combine all local stiffness matrices to form a global stiffness matrix:

⎥⎥

⎥⎥

⎦

⎤

⎢⎢

⎢⎢

⎣

⎡

−

−

−

−

−

−

=

5 5 0 0

5 8 3 0

0 3 4 1

0 0 1 1 7 . 66 K

whereas the heat rate vector is

[

]

= R

The penalty parameter C can be estimated by

4

4 66.7 8 10

10

max × = × ×

= _ij

C K

Hence, the finite element system can be solved as followed:

⎪⎪

⎭

⎪⎪

⎬

⎫

⎪⎪

⎩

⎪⎪

⎨

⎧

×

×

=

⎪⎪

⎭

⎪⎪

⎬

⎫

⎪⎪

⎩

⎪⎪

⎨

⎧

⎥⎥

⎥⎥

⎦

⎤

⎢⎢

⎢⎢

⎣

⎡

−

−

−

−

−

−

4 4

3 2 1

10 672 , 10

0 0

800 25

005 , 80 5 0 0

5 8

3 0

0 3

4 1

0 0

1 375 . 1 7 . 66

T T T T

[

]

= T

`

7.6 Boundary Element Method ⎯ An Introduction

• The Boundary Element Method (BEM) is a relatively new for PDE, where it consists only boundary elements, either line (2-D) or surface (3-D) elements.

• Consider a Laplace equation for a steady-state potential flow problems (ψ is the stream function):

2 = 0

∇ ψ

• For a multi-dimensional case, this equation leads to an analytical solution:

r ln1 2

1

= π

ψ∗ _(7.31)

FIGURE 7.15 Boundary elements

• Eq. (7.31) can be discretised to yield:

⎟⎠

⎜ ⎞

⎝⎛

⎟⎠

⎜ ⎞

⎝

⎛

∂

= ∂

⎟⎟⎠

⎜⎜ ⎞

⎝

⎛

∂

+

∑ ∫

∑ ∫

=

∗

= j

N

j j

j N

j j

i ds

ds n

n ψ ψ

ψ ψ ψ

1 1

2

1 (7.32)

FIGURE 7.16 Example of an analysis using BEM

Exercises

1. A quarter of disc having a radius of 8 cm as shown below has a variation of temperature at the boundary aligned with the principal axes, while the temperature is fixed at 100°C at its outer radius of 8 cm. If the temperature distribution follows the Laplace equation:

2 =0

∇ T

20°C T1 T2 T3

T4 T5 T6

T7 T8

100°C

0°C 40°C 65°C 85°C 100°C

50°C 75°C

100°C

100°C

By using the grid as shown:

a. Obtain a system of algebraic equations using the finite difference method, b. Solve the system to obtain T_i .

2. By using the implicit technique, solve the following heat conduction problem:

), , ( )

,

( ₂

2

t x x T t

x

tT ∂

= ∂

∂

∂ α 0<x<10, t > 0

using the following boundary conditions:

. 0 ) 0 , ( , 100 ) , 10 ( , 0 ) , 0

( t = T t = T x =

T

Use the constant α=10, the time step Δt=0.1, and a model of 10 grid including the grid at the boundary. Compar this solution with the solution obtained by using Δt=0.3.