Transcript 06_finite_elements_basics

Finite Elements




Basic formulation
Basis functions
Stiffness matrix
Poisson‘s equation
 Regular grid
 Boundary conditions
 Irregular grid
 Numerical Examples
Scope: Understand the basic concept of the finite element
method with the simple-most equation.
Finite element method
1
Formulation
Let us start with a simple linear system of equations
Ax  b
|*y
and observe that we can generally multiply both sides of this
equation with y without changing its solution. Note that x,y
and b are vectors and A is a matrix.
 yAx  yb
y  n
We first look at Poisson’s equation
 u( x)  f ( x)
where u is a scalar field, f is a source term and in 1-D
2

  2  2
x
Finite element method
2
Poisson‘s equation
We now multiply this equation with an arbitrary function v(x),
(dropping the explicit space dependence)
 uv  fv
... and integrate this equation over the whole domain. For
reasons of simplicity we define our physical domain D in the
interval [0, 1].
  uv   fv
D
1
D
1
  uvdx   fvdx
0
Das Reh springt hoch,
das Reh springt weit,
warum auch nicht,
es hat ja Zeit.
0
... why are we doing this? ... be patient ...
Finite element method
3
Discretization
As we are aiming to find a numerical solution to our problem it is
clear we have to discretize the problem somehow. In FE problems
– similar to FD – the functional values are known at a discrete set
of points.
... regular grid ...
... irregular grid ...
Domain D
The key idea in FE analysis is to approximate all functions in
terms of basis functions j, so that
N
u  u~   cij i
i 1
Finite element method
4
Basis function
N
~
u u
 ciji
i 1
where N is the number nodes in our physical domain and ci are real
constants.
With an appropriate choice of basis functions ji, the coefficients ci
are equivalent to the actual function values at node point i. This – of
course – means, that ji=1 at node i and 0 at all other nodes ...
Doesn’t that ring a bell?
Before we look at the basis functions, let us ...
Finite element method
5
Partial Integration
... partially integrate the left-hand-side of our equation ...
1
1
0
0
  uvdx   fvdx
1
1
  (  u )vdx  uv   vudx
1
0
0
0
we assume for now that the derivatives of u at the boundaries vanish
so that for our particular problem
1
1
0
0
  (  u )vdx   vudx
Finite element method
6
…
... so that we arrive at ...
1
1
0
0
 uvdx   fvdx
... with u being the unknown. This is also true for our
approximate numerical system
1
1
0
0
~vdx  fvdx

u


... where ...
N
u~   cij i
i 1
was our choice of approximating u using basis functions.
Finite element method
7
Partial integration
1
1
0
0
~vdx  fvdx

u


... remember that v was an arbitrary real function ...
if this is true for an arbitrary function it is also true if
v jj
... so any of the basis functions previously defined ...
... now let’s put everything together ...
Finite element method
8
The discrete system
N
The ingredients:
v  jk
u~   cij i
i 1
1
1
0
0
~vdx  fvdx

u


 n

cij i j k dx   f j k dx
0  
i 1

0
1
1
... leading to ...
Finite element method
9
The discrete system
... the coefficients ck are constants so that for one particular
function jk this system looks like ...
1
n
1
 c  j j dx   f j dx
i 1
i
i
0
k
k
0
... probably not to your surprise this can be written in matrix form
bi Aik  gk
A b  gk
T
ik i
Finite element method
10
The solution
... with the even less surprising solution
 
T 1
ik
bi  A
gk
remember that while the bi’s are really the coefficients of the basis
functions these are the actual function values at node points i as well
because of our particular choice of basis functions.
This become clear further on ...
Finite element method
11
Basis functions
we are looking for functions ji
with the following property
... otherwise we are
free to choose any
function ...
10
The simplest choice
are of course linear
functions:
7
1
j i ( x)  
0
for x  xi
for x  x j , j  i
9
8
6
5
+ grid nodes
blue lines – basis
functions ji
4
3
2
1
Finite element method
12
Basis functions - gradient
To assemble the stiffness matrix we need the gradient (red) of the basis
functions (blue)
10
9
8
7
6
5
4
3
2
1
Finite element method
13
Stiffness matrix
Knowing the particular form of the basis functions we can now
calculate the elements of matrix Aij and vector gi
1
n
1
 c  j j dx   f j dx
i 1
i
i
k
0
k
0
bi Aik  gk
1
1
Aik   j i j k dx
0
g k   f j k dx
0
Note that ji are continuous functions defined in the interval [0,1], e.g.
 x  xi 1
x x
 i i 1
x  x
j i ( x)   i 1
 xi 1  xi
 0


Finite element method
for xi 1  x  xi
for xi  x  xi 1
Let us – for now – assume a
regular grid ... then
elsewhere
14
Stiffness matrix –regular grid
 x  xi 1
x x
 i i 1
x  x
j i ( x)   i 1
 xi 1  xi
 0


for xi 1  x  xi
for xi  x  xi 1
elsewhere

x
~
 dx  1
~

x
~
j i ( x )  1 
 dx
 0

for  dx  ~
x 0
for
0~
x  dx
elsewhere
... where we have used ...
ji
~
x  x  xi
dx  xi  xi 1
xi
dx
Finite element method
15
Regular grid - gradient
 1 / dx

~
j i ( x )   1 / dx
 0

for  dx  ~
x 0
for 0  ~
x  dx
elsewhere
~
x  x  xi
dx  xi  xi 1
ji
1/dx
xi
dx
-1/dx
Finite element method
16
Stifness matrix - elements
10
1
9
8
7
6
Aik   j i j k dx
5
4
3
0
2
1
... we have to distinguish various cases ... e.g. ...
1
A11   j1j1dx 
x1  dx
 j j dx  
1
0
x1
1
x2
A22   j 2j 2 dx 
0
1
 2
dx
Finite element method
0
x1  dx
1
x1
dx

0
1
dx 
dx
x2  dx
 j j dx   j j dx
2
2
x2  dx
1
dx

dx dx2
1 1
1
dx  2
dx dx
dx
0
2
x2
dx

2
dx 
2
dx
17
Stiffness matrix
10
1
9
8
7
6
Aik   j i j k dx
5
4
3
0
2
1
......and
and......
1
A12   j1j 2 dx 
0
1
 2
dx
x1  dx
 j j dx  
1
x1
dx

0
dx 
x1  dx
2
x1
1 1
dx
dx dx
1
dx
A21  A12
... so that finally the stiffness matrix looks like ...
Finite element method
18
Stiffness matrix
10
1
Aik   j i j k dx
0
9
8
7
6
5
4
3
2
1
 1 1



 1 2 1

1

Aij  


dx 
 1 2  1


 1 1 

... so far we have ignored sources and boundary conditions ...
Finite element method
19
Boundary conditions - sources
... let us start restating the problem ...
 u( x)  f ( x)
... which we turned into the following formulation ...
1
n
1
 c  j j dx   f j dx
i 1
i
i
k
0
k
0
... assuming ...
N 1
N
u~   cij i
i 1
with b.c.
u~   cij i  u(0)j1  u (1)j N
i 2
where u(0) and u(1) are the values at the boundaries of the domain [0,1].
How is this incorporated into the algorithm?
Finite element method
20
Boundary conditions
1
n
1
 u( x)  f ( x)
 c  j j dx   f j dx
i 1
i
i
k
k
0
0
... which we turned into the following formulation ...
n 1
1
1
1
1
 c  j j dx   f j dx  u(0) j j dx  u(1) j j dx
i 2
i
i
0
k
k
0
1
k
0
n
k
0
... in pictorial form ...
AT
b
= g
boundary condition
source heterogeneity (f)
=
boundary condition
... the system feels the boundary conditions through the (modified) source term
Finite element method
21
Numerical Example
 u( x)  f ( x)
Domain: [0,1]; nx=100;
dx=1/(nx-1);f(x)=d(1/2)
Boundary conditions:
u(0)=u(1)=0
Matlab FEM code
% source term
s=(1:nx)*0;s(nx/2)=1.;
% boundary left u_1 int{ nabla phi_1 nabla phij }
u1=0;
s(1) =0;
% boundary right u_nx int{ nabla phi_nx nabla phij }
unx=0; s(nx)=0;
% assemble matrix Aij
A=zeros(nx);
Matlab FD code
f(nx/2)=1/dx;
for it = 1:nit,
uold=u;
du=(csh(u,1)+csh(u,-1));
u=.5*( f*dx^2 + du );
u(1)=0;
u(nx)=0;
for i=2:nx-1,
for j=2:nx-1,
if i==j,
A(i,j)=2/dx;
elseif j==i+1
A(i,j)=-1/dx;
elseif j==i-1
A(i,j)=-1/dx;
else
A(i,j)=0;
end
end
end
fem(2:nx-1)=inv(A(2:nx-1,2:nx-1))*s(2:nx-1)';
fem(1)=u1;
fem(nx)=unx;
end
Finite element method
22
Regular grid
 u( x)  f ( x)
Domain: [0,1]; nx=100;
dx=1/(nx-1);f(x)=d(1/2)
Boundary conditions:
u(0)=u(1)=0
FD (red) - FEM (blue)
0.25
0.2
Matlab FD code (red)
u(x)
0.15
0.1
Matlab FEM code (blue)
0.05
0
Finite element method
0
0.1
0.2
0.3
0.4
0.5
x
0.6
0.7
0.8
0.9
1
23
Regular grid - non zero b.c.
 u( x)  f ( x)
FD (red) - FEM (blue) -> Regular grid
Domain: [0,1]; nx=100;
dx=1/(nx-1);f(x)=d(1/2)
Boundary conditions:
u(0)=0.15
u(1)=0.05
0.4
0.35
0.3
% Quelle
0.25
s=(1:nx)*0;s(nx/2)=1.;
% Randwert
links
0.2
phij }
u(x)
Matlab FD code (red)
u1=0.15;
0.15
int{ nabla phi_1
nabla
s(2) =u1/dx;
% Randwert links
phij0.1
}
Matlab FEM code (blue)
u_1
u_nx int{ nabla phi_nx nabla
unx=0.05; s(nx-1)=unx/dx;
0.05
0
Finite element method
0
0.1
0.2
0.3
0.4
0.5
x
0.6
0.7
0.8
0.9
1
24
Stiffness – irregular grid
10
9
1
Aik   j i j k dx
8
7
6
0
5
4
3
2
1
1
x1  h1
x1  h1
0
x1
x1
1 1
dx
h1 h1
2
+
3
+
A12   j1j 2 dx 
 j1j 2dx 
1 1
1
 2  dx 
 A21
h1 0
h1

h
Aii 
Finite element method
1
1

hi 1 hi
i=1
+
h1 h2
4
+
h3
5
6 7
+ + +
h4 h5 h6
25
Example
Stiffness matrix A
 u( x)  f ( x)
for i=2:nx-1,
for j=2:nx-1,
if i==j,
A(i,j)=1/h(i-1)+1/h(i);
Domain: [0,1]; nx=100;
dx=1/(nx-1);f(x)=d(1/2)
Boundary conditions:
u(0)=u0; u(1)=u1
elseif i==j+1
A(i,j)=-1/h(i-1);
elseif i+1==j
A(i,j)=-1/h(i);
else
A(i,j)=0;
i=1 2
3 4 5 6
7
+ +
+
+ +
+
+
h1 h2 h3 h4 h5 h6
Finite element method
end
end
end
26
Irregular grid – non zero b.c.
 u( x)  f ( x)
FEM on Chebyshev grid
Domain: [0,1]; nx=100;
dx=1/(nx-1);f(x)=d(1/2)
Boundary conditions:
u(0)=0.15
u(1)=0.05
FD (red) - FEM (blue)
0.4
0.35
0.3
0.25
u(x)
Matlab FD code (red)
0.2
0.15
Matlab FEM code (blue)
0.1
0.05
+ FEM grid points
Finite element method
0
0
0.1
0.2
0.3
0.4
0.5
x
0.6
0.7
0.8
0.9
1
27
Summary
In finite element analysis we approximate a function defined in
a Domain D with a set of orthogonal basis functions with
coefficients corresponding to the functional values at some
node points.
The solution for the values at the nodes for some partial
differential equations can be obtained by solving a linear
system of equations involving the inversion of (sometimes
sparse) matrices.
Boundary conditions are inherently satisfied with this
formulation which is one of the advantages compared to finite
differences.
Finite element method
28