Lecture 5 – Finite Element Method

Preliminary Math:

Integration by parts: | -

Integration theorem: ∀ (i.e., for arbitrary g(x))

If 0 => ≡ 0, ∈ ,

The Method of Weighted Residual:

Consider the problem of curve fitting where one might wish to fit the function sin ,
for ∈ 0, . Let ̅ be a polynomial approximation. We try to find a, b and
c.

Define the error or residual by ̅ sin

Consider weight function w(x) f(x)=sin(x)

1

0.9
We require: 0 0.8

0.7

0.6
f(x)

0.5
For several w(x): w1(x), w2(x) ⋯ wi(x) 0.4

0.3

If dx 0 (weak form) 0.2

0.1

We believe R(x) is approximately equal to zero. 0

0 0.5 1 1.5 2 2.5 3 3.5
x

Some common choices for the weight function are:

, is pre-selected points (collocation method)

0, ∉ subdomain
(subdomain method)
1, ∈ subdomain

(method of moments)

, are the constants in the assumed solution (method of least squares)

∅ , ̅ ∅ ∑ ∅ ,∅ is shape function (Galerkin’s method)

1
 
Where, the dirac delta function has the property

δ dx for ∈ ,

For the f(x) = sin(x) example: sin

Use collocation method:

0 0, /2 0, 0

What’s the essence of this method?

Answer: It is to force the residual equal to zero at fixed points. To force R(x) at x = 0, x =
π/2, x = π to be zero.

Therefore,

∙0 ∙0 sin 0 0

∙ ∙ sin 0

∙ ∙ sin 0

1 0 0 sin 0 0
=> 1 sin =>
1 sin

Then ̅

Approximation
1
exact function
0.9 collocation fit

0.8

0.7

0.6
f(x)

0.5

0.4

0.3

0.2

0.1

0
0 0.5 1 1.5 2 2.5 3 3.5
x

2
 
Or we use the mixed weight function (collocation + subdomain):

0 , ∙1 ,

Therefore,

∙0 ∙0 sin 0 0

cos | 0

∙ ∙ sin 0

1 0 0 0
sin 0
=> 1 cos 0 cos =>
1 sin

Then ̅

Approximation
1
exact function
0.9 collocation fit
mixed fit
0.8

0.7

0.6
f(x)

0.5

0.4

0.3

0.2

0.1

0
0 0.5 1 1.5 2 2.5 3 3.5
x

3
 
Use the moment weight function:

, ∙ ,

Therefore,

cos | 0

cos | 0

cos 2 sin 2 cos | 0

2
=>
4

Matlab code:

syms PI;

K = [PI, PI^2/2, PI^3/3; PI^2/2, PI^3/3, PI^4/4; PI^3/3, PI^4/4, PI^5/5];

f = [2; PI; PI^2-4]; 1.2

Approximation

exact function

u = K \ f;
moment fit
1

0.8

0.6
12 10 /
f(x)

=> 60 12 / 0.4

60 12 / 0.2

̅ -0.2
0 0.5 1 1.5 2 2.5 3 3.5
x

4
 
Apply this Weighted Residual technique to differential equations:

B.C. 0 0, 1 0, ∀ ∈ 0,1

The analytical solution is

Now let be an approximation to , then in general,

0

So we define the residual as

We try to minimize the error or residual by using the weighted residual method.

We choose approximation as:

∅ ∅

Some choices for the ∅ would be:

∅ 0, ∅ 1 ,∅ 1

If we use Galerkin method by choosing ∅ as weight function:

For the one term Galerkin solution

∅ ∅ 0 1

2 1 , 1

1 0

=>

2 1 1 0

=> 5/18

5
 
 5 1
18
Approximation
0.08
exact function
0.07 one term Galerkin

0.06

0.05

u(x)
0.04

0.03

0.02

0.01

0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
x

Two term Galerkin solution:

∅ ∅ ∅ 0 1 1

2 6

1 0

1 0

=>

So, 1 1
Approximation
0.08
exact function
0.07 one term Galerkin
two term Galerkin

0.06

0.05
u(x)

0.04

0.03

0.02

0.01

0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
x

6
 
The finite element by the Garlerkin approach:

Consider again the 1D cantilever problem, where

For the finite element method, as the finite difference method, the domain is firstly discretized
into subdomains (elements).

Linear Element

On an element x1 ≤ x ≤ x2, the primary variable is assumed to be of the form (for a linear
element)

To replace the constants and by nodal values of the function u, we evaluate at x1 and x2,

1
Or =>
1 1 1

Then

1
1 1 1
≡

, are shape functions or interpolation functions.

Then by Galerkin’s method of weighted residuals

7
 
 0

=>

And with integration by parts

So,

1, 0, 0, 1

|
|
|

=>

=>

∙ ∙

∙ ∙

Define ∆

8
 
 1 1
∆ ∆ 1 1
∙∆
1 1 ∆ 1 1
∆ ∆
∆ 1
2 1

The whole linear element equation for the bar element is

1 1 ∆ 1
∆ 1 1 2 1

1 1 0 0 0 1
1 1 0 0 0 ∆ 1
0 0 0 0 0 0 0
∆ 0 0 0 0 0 2 0 0
0 0 0 0 0 0 0

0 0 0 0 0 0 0
0 1 1 0 0 ∆ 1
0 1 1 0 0 1
∆ 0 0 0 0 0 2 0 0
0 0 0 0 0 0 0
Add up:

1 1 0 0 0 1
1 2 1 0 0 ∆ 2 0
0 1 2 1 0 2 0
∆ 0 0 1 2 1 2 2 0
0 0 0 1 1 1

By applying boundary condition: u0=0, σx4=0.

=>

9
 
 1 0 0 0 0 0
1 2 1 0 0 ∆ 2
0 1 2 1 0 2
∆ 0 0 1 2 1 2 2
0 0 0 1 1 1

By seeing the matrix form, can you recognize it?

Iso-parametric Elements

Consider again the derivation of the linear function on the typical element. If we map the x
coordinate to a ξ coordinate that it ranges from -1 ≤ ξ ≤ 1 on the element. We would lie to give a
standard element irrelevant to , .

Then the coordinate transformation is:

Where a, b are constants and meet

1 1
=>
1 1
1 1
=>
1 1

=>

1 1 1
Or
1 1
Therefore, we can rewrite shape function for the linear element:

2 2 1
∆ 2

10
 
 2 2 2 2 1
∆ ∆ 2

Or simply

1
2
1
2
After the coordinate transformation, the integral of a function f(x) will be transformed to

Where f(x) is an arbitrary function and ≡ / is the Jacobian function of the

transformation x(ξ).

The Jacobian function:

1 1 1 1 ∆
≡
2 2 2 2 2 2

Hence

∆
2

11
 
The Quadratic Element

In the intrinsic coordinate ξ,

1 ∙ 1 ∙ 1

0 ∙ 0 ∙ 0

1 ∙ 1 ∙ 1

1 1 1
1 0 0
1 1 1
=>

1 0 2 0
1 0 1
2
1 2 1
=>

1 0 2 0
1 ∙ 1 0 1
2
1 2 1

1 1
1 1 1 1
2 2

(1) (2) (3)

12
 
Term (1)

2
∆

2
∆

1 1 1 1
2
2 2 2 2
2 1 1
2 2 2
∆ 2 2
1 1 1 1
2
2 2 2 2
7 8 1
8 16 8
3∆
1 8 7
where

1 1 2
1 1 ∆
2 2

Term (2)

Term (3)

∆ 1
4
6
1

13
 
Put three terms together

7 8 1 ∆ 1
8 16 8 4 0
3∆ 6
1 8 7 1

7 8 1 0 0 0 0 0 0 1
8 16 8 0 0 0 0 0 0 4 0
1 8 7 0 0 0 0 0 0 1
0 0 0 0 0 0 0 0 0 0 0
∆
∆
0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0
0 0 7 8 1 0 0 0 0 1
0 0 8 16 8 0 0 0 0 4 0
∆
∆
0 0 1 8 7 0 0 0 0 1
0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0

Add up:

7 8 1 0 0 0 0 0 0 1
8 16 8 0 0 0 0 0 0 4 0
1 8 14 8 1 0 0 0 0 2 0
0 0 8 16 8 0 0 0 0 4 0
∆
0 0 1 8 14 8 1 0 0 2 0
∆
0 0 0 0 8 16 8 0 0 4 0
0 0 0 0 1 8 14 8 1 2 0
0 0 0 0 0 0 8 16 8 4 0
0 0 0 0 0 0 1 8 7 1

14
 
2D FEM

The governing equations – differential equations of equilibrium:

B.C.

Where ,

Geometric equation:

, ,

Hooke’s law:

For plane strain problem:

1 0
1 0
1 1 2 1 2
0 0
2
For the plane stress problem:

1 0
1 0
1 1
0 0
2
Substitute geometric equation and Hooke’s law to the differential equation of equilibrium:

15
 
For the u, v approximation u ( x, y )

u
̅ u3 u2
u1
Substitute , , , , , to y
3
2
1
x

=>

1
1
1

syms , syms , syms ;

syms , syms , syms ;

A = [1 ;1 ;1 ];

B = inv(A);

=>

1
2

A is the area of the element:

1
, 1
2

, , , , , ,

1
,
2

16
 
 1
,
2
1
,
2

, 1, , 0, , 0

, 0, , 1, , 0

, 0, , 0, , 1

What are the shape functions for , ?

0 0 0
0 0 0

Then, the strain:

0 0
0 0 0
0 0
0 0 0

0 0 0

0 0 0

=>

17
 
The weak form

As before, we integrate by parts by using Green’s 1st theorem and get:

Note:

Therefore

18
 
 0 0
0 0
0 0
0 0
0 0
0 0

1 0 0 0
0 0 0
2
1 1
0 0 0 0
∆ ∆
1 ∆ 0 ∆ 0 0 0 1 1
0 ∆ 0 0 0 ∆ 0 0 0 0
∆ ∆ ∆ ∆
∆ ∆ 0 ∆ ∆ 0
1 1 1 1
0 0
∆ ∆ ∆ ∆

1 1
0
∆ ∆
1 1
0 1 1
∆ ∆
0 0 0 0
1 1 0 ∆ ∆
∆ ∆ 0 0 1 1
∆ 1 0 0 0 0 0
2 1 1 1 2 1 2 ∆ ∆
0 0 0 0
∆ 2 1 1 1 1
1 0 0
∆ ∆ ∆ ∆
0 0
∆
1
0 0
∆

19
 
The order of assembly

Node order on elements:

1 2 3
(1) 1 2 3
(2) 3 2 4

20