Introduction FEM, 1D-Example
Universität Stuttgart
Institut für Wasserbau, Lehrstuhl für Hydromechanik und Hydrosystemmodellierung
home/lehre/VL-MHS-1-E/FOLIEN/VORLESUNG/3_FEM_INTRO/cover_sheet.tex
page 1 of 25
. – p.1/25
Table of contents 1D Example - Finite Element Method 1. 1D Setup – Geometry 2. Governing equation 3. General Derivation of Finite Element Method 4. Shape functions 5. Assembling the matrix 6. Calculation of the matrix entries 7. Solving the set of equations
Universität Stuttgart
Institut für Wasserbau, Lehrstuhl für Hydromechanik und Hydrosystemmodellierung
home/lehre/VL-MHS-1-E/FOLIEN/VORLESUNG/3_FEM_INTRO/toc_fem_intro.tex
page 2 of 25
. – p.2/25
1D-Example: Finite Element Method Given : 1. Geometry ∆x = 1 m
Hlef t Q
acements
x 4m
2
1
Universität Stuttgart
3
4
5
1
2
3
4
∆x
∆x
∆x
∆x
2. Permeability kf = 10−5 m/s 3. Boundary conditions Q = 10−4 m3 /s Hlef t = 50 m 4. Notation 1 denotes node 1 1 denotes element 1
Institut für Wasserbau, Lehrstuhl für Hydromechanik und Hydrosystemmodellierung
home/lehre/VL-MHS-1-E/FOLIEN/VORLESUNG/3_FEM_INTRO/1d_setup_fem.tex
page 3 of 25
. – p.3/25
Governing Equations (FEM – 1D) Continuity equation ∇·v−q =0
Momentum equation −→ (Darcy equation) v = −kf ∇h −→
∇ · (−kf ∇h) − q = 0
∂ ∂h 1D −→ (−kf ) − q = 0 ∂x ∂x Assumption: No sinks or sources −→ Universität Stuttgart
∂2h −kf · 2 = 0 ∂x
(Laplace-Equation)
Institut für Wasserbau, Lehrstuhl für Hydromechanik und Hydrosystemmodellierung
home/lehre/VL-MHS-1-E/FOLIEN/VORLESUNG/3_FEM_INTRO/gov_eqn_fem1.tex
page 4 of 25
. – p.4/25
The weighting function (FEM – 1D) ˜ We now introduce an approximate function h(x) Inserting this approximation into the Laplace equation yields a residuum (error) ε: ˜ −q =ε ∇ · (−kf · ∇h) We introduce a weighting function W. This function is multiplied with our residual and integrated over the whole domain. We choose this weighting function such, that the integral is zero. Z Z n o ˜ − q dΩ = 0 W · ε dΩ = W · ∇ · (−kf · ∇h) Ω
Ω
=
Z Ω
Universität Stuttgart
n
o
˜ dΩ − W · ∇ · (−kf · ∇h)
Z
W · q dΩ = 0
Ω
Institut für Wasserbau, Lehrstuhl für Hydromechanik und Hydrosystemmodellierung
home/lehre/VL-MHS-1-E/FOLIEN/VORLESUNG/3_FEM_INTRO/starting_equ.tex
page 5 of 25
. – p.5/25
Starting equation II (FEM – 1D) Use the product rule of differentiation in the reverse way Z n o ˜ dΩ = W · ∇ · (−kf · ∇h) Ω
Z
˜ ∇[W · (−kf · ∇h)]dΩ −
Ω
Z
˜ [∇W ]T · (−kf · ∇h)dΩ,
Ω
and apply Gauss theorem on the first term of the right hand side Z
˜ ∇[W · (−kf · ∇h)]dΩ =
Ω
Universität Stuttgart
Z
˜ · n dΓ. W · (−kf · ∇h)
Γ
Institut für Wasserbau, Lehrstuhl für Hydromechanik und Hydrosystemmodellierung
home/lehre/VL-MHS-1-E/FOLIEN/VORLESUNG/3_FEM_INTRO/starting_equ2.tex
page 6 of 25
. – p.6/25
Starting equation III (FEM – 1D) We collect all the terms Z Z Z ˜ · n dΓ − [∇W ]T · (−kf · ∇h)dΩ ˜ W · (−kf · ∇h) − W · q dΩ = 0 Γ
Ω
Ω
and rearrange them. On the left hand side, we have our unknowns on the right hand side we have the boundary conditions and the source and sink terms Z Z Z ˜ ˜ · n dΓ + W · q dΩ. [∇W ]T · (kf · ∇h)dΩ = W · (kf · ∇h) Ω
Γ
Ω
Please note, that at this point, we are still looking at the domain as a whole. Universität Stuttgart
Institut für Wasserbau, Lehrstuhl für Hydromechanik und Hydrosystemmodellierung
home/lehre/VL-MHS-1-E/FOLIEN/VORLESUNG/3_FEM_INTRO/starting_equ3.tex
page 7 of 25
. – p.7/25
Starting equation IV (FEM – 1D) Finally we can rewrite our boundary term, by introducing the ve˜ and get the weak formulation: locity again [vo = (−kf · ∇h)] Z
˜ dΩ = − [∇W ]T · kf · [∇h]
Ω
Z
W · vo · n dΓ +
Γ
Z
W · q dΩ
Ω
Weak formulation (mathematical concept): We reduced as in the IFDM the order of the differentiation. Derivatives don’t have to be available everywhere, because every term is an integral term. Existence and uniqueness of the problems can be shown. Universität Stuttgart
Institut für Wasserbau, Lehrstuhl für Hydromechanik und Hydrosystemmodellierung
home/lehre/VL-MHS-1-E/FOLIEN/VORLESUNG/3_FEM_INTRO/starting_equ4.tex
page 8 of 25
. – p.8/25
Shape functions (FEM – 1D) We now have to choose how our approximate solution shall look like. We choose it such, that it can be constructed from the shape functions and node values ˜ h(x) =
n X j=1
hˆj · Nj (x).
At this point we also introduce our discretization. We choose our elements and with that our nodes. n is here the number of nodes we have chosen. In the setting which is either called standard Galerkin method or Bubnov-Galerkin method the shape functions Ni are also used as the weighting functions Wi Wi = N i . Universität Stuttgart
Institut für Wasserbau, Lehrstuhl für Hydromechanik und Hydrosystemmodellierung
home/lehre/VL-MHS-1-E/FOLIEN/VORLESUNG/3_FEM_INTRO/ansatz_funct.tex
page 9 of 25
. – p.9/25
Equation in matrix from (FEM – 1D) Z
"
[∇Ni ]T · kf · ∇
Ω
As
j=1
n X j=1
Pn
hˆj ·
n X j=1
#
hˆj · Nj (x) dΩ =
Z
Ni · vo · n dΓ −
Γ
Z
Ni · q dΩ
Ω
ˆ j is not a function, we can take it out of the integral. h
Z
[∇Ni ]T · kf · [∇Nj (x)] dΩ = {z } | Aij
Z
|Γ
Ni · vo · n dΓ − {z
Z
Ni · q dΩ
Ω
bi
}
The equation has a vector-matrix structure representing a system of equations. The solution is our approximate function. n X ˆ=b hˆj · Aij = bi ⇐⇒ A · h Universität Stuttgart
j=1
Institut für Wasserbau, Lehrstuhl für Hydromechanik und Hydrosystemmodellierung
home/lehre/VL-MHS-1-E/FOLIEN/VORLESUNG/3_FEM_INTRO/stiffnessmat.tex
page 10 of 25
. – p.10/25
Shape functions (FEM – 1D)
cementsConstruction of the local shape functions N1 and N2 . We can use linear functions, as we have only first order derivatives in the “weak approximation”. Element 1: (1)
N1 hˆ1 ∆x
1 x=0
hˆ2 1
2 x=1
y=1
(1)
y=1
N1
(1)
x
1 x=0
∆x=1m
Universität Stuttgart
(1)
N2
1
2 x=1
x
=1−x
N2 = x 1 : Element 1
1,2 : nodes
∆x=1m
Institut für Wasserbau, Lehrstuhl für Hydromechanik und Hydrosystemmodellierung
home/lehre/VL-MHS-1-E/FOLIEN/VORLESUNG/3_FEM_INTRO/shape_funct.tex
page 11 of 25
. – p.11/25
Shape functions II (FEM – 1D) Construction of the local shape functions N1 and N2 , which are also called trial functions, interpolation functions or basis functions. We can use linear functions, as we have only first order derivatives in the “weak approximation”. These shape functions need to be adjusted for the other elements. The slope remains the same but the y-intercept changes. Element 2 for example has the following shape functions (2)
N2 = 2 − x
Universität Stuttgart
(2)
N3 = x − 1.
Institut für Wasserbau, Lehrstuhl für Hydromechanik und Hydrosystemmodellierung
home/lehre/VL-MHS-1-E/FOLIEN/VORLESUNG/3_FEM_INTRO/shape_funct2.tex
page 12 of 25
. – p.12/25
Shape functions III (FEM – 1D) By adding the products of piezometric head at the nodes hˆj and the local ansatz functions Nj one gets an approximate solution ˜ h(x) between the nodes Element 1: hˆ1 n X (1) ˜ h i i PSfrag replacements ˜ (x) = hˆ · N (x) h j
j
j=1
˜ (1) (x) = h1 (1 − x) + h2 x h
˜ (1) (x) = h1 + (h2 − h1 )x h
In this case we have a linear in- ∆ x terpolation. Universität Stuttgart
1 x=0
1
hˆ2
2 x=1
x
∆x=1m
Institut für Wasserbau, Lehrstuhl für Hydromechanik und Hydrosystemmodellierung
home/lehre/VL-MHS-1-E/FOLIEN/VORLESUNG/3_FEM_INTRO/shape_funct3.tex
page 13 of 25
. – p.13/25
System of equations (FEM – 1D)
A11 A21 A31 A 41 A51
A12 A22 A32 A42 A52
A13 A23 A33 A43 A53
A14 A24 A34 A44 A54
(1) Q1 h1 A15 (2) (1) A25 h2 Q2 − Q2 (3) (2) h3 = Q − Q A35 3 3 h (4) (3) A45 Q − Q 4 4 4 (4) A55 h5 −Q5
h1 is known (Dirichlet boundary conditions). → Delete the corresponding row and column. (1) A22 A23 A24 A25 h2 −h1 K21 A h A A A 0 32 33 34 35 3 = A42 A43 A44 A45 h4 0 A52 A53 A54 A55 h5 −Q5 Universität Stuttgart
Institut für Wasserbau, Lehrstuhl für Hydromechanik und Hydrosystemmodellierung
home/lehre/VL-MHS-1-E/FOLIEN/VORLESUNG/3_FEM_INTRO/set_of_eqn.tex
page 14 of 25
. – p.14/25
Calculating matrix entries (FEM – 1D) Example 1: Calculation of the matrix entry A23 (weighting function of node 2, shape function of node 3). The stiffness integral can be divided into different parts. At each element of A23 , except for element 2, the shape function or the weighting function equals zero. (1)
(2)
(3)
A23 = K23 + K23 + K23 (1)
(3)
K23 = K23 = 0
(2)
A23 = K23 =
Z2m
=⇒ (2)
(2)
A23 = K23 (2)
∂N2 ∂N3 · kf · dx ∂x ∂x
1m Universität Stuttgart
Institut für Wasserbau, Lehrstuhl für Hydromechanik und Hydrosystemmodellierung
home/lehre/VL-MHS-1-E/FOLIEN/VORLESUNG/3_FEM_INTRO/matr_entr_23a.tex.tex
page 15 of 25
. – p.15/25
Calculating matrix entries II (FEM – 1D) A23 =
(2) K23
=
Z2m
∂(x − 1) ∂(2 − x) · kf · dx ∂x ∂x
1m
(2)
A23 = K23
Z2m = −1 · kf · 1 dx = kf · [−x]2m 1m 1m
(2)
A23 = K23 = kf [−(2m − 1m)] = −kf · 1m A23 has the units
Universität Stuttgart
m2 [ s ].
Institut für Wasserbau, Lehrstuhl für Hydromechanik und Hydrosystemmodellierung
home/lehre/VL-MHS-1-E/FOLIEN/VORLESUNG/3_FEM_INTRO/matr_entr_23b.tex
page 16 of 25
. – p.16/25
Calculating matrix entries III (FEM – 1D) Example 2: Calculation of the matrix entry A55 . The stiffness integral can be divided into different parts. (1)
(2)
(3)
(4)
A55 = K55 + K55 + K55 + K55 (1)
(2)
(3)
(4) K55
Z4
K55 = K55 = K55 = 0
A55 =
=
=⇒ (4)
(4)
A55 = K55 (4)
∂N5 ∂N5 · kf · dx ∂x ∂x
3
Universität Stuttgart
Institut für Wasserbau, Lehrstuhl für Hydromechanik und Hydrosystemmodellierung
home/lehre/VL-MHS-1-E/FOLIEN/VORLESUNG/3_FEM_INTRO/matr_entr_55a.tex
page 17 of 25
. – p.17/25
Calculating matrix entries IV (FEM – 1D) A55 =
(4) K55
=
Z4
∂(x − 3) ∂(x − 3) · kf · dx ∂x ∂x
3
(4)
A55 = K55 =
Z4
1 · kf · 1 dx = kf · [x]43
3
(4)
A55 = K55 = kf [(4 − 3)] = kf
Universität Stuttgart
Institut für Wasserbau, Lehrstuhl für Hydromechanik und Hydrosystemmodellierung
home/lehre/VL-MHS-1-E/FOLIEN/VORLESUNG/3_FEM_INTRO/matr_entr_55b.tex
page 18 of 25
. – p.18/25
Calculating matrix entries V (FEM – 1D) Example 3: Calculation of the matrix entry A22 . The stiffness integral can be divided into different parts. In this case, only the weighting function and the shape function in element 3 and 4 equal zero. (1)
(2)
(3)
(4)
A22 = K22 + K22 + K22 + K22 (3)
(4)
K22 = K22 = 0 (1)
K22 =
=⇒
(1)
Z1
∂N2 ∂N2 · kf · dx ∂x ∂x
=
Z1
(1)
(2)
A22 = K22 + K22 (1)
0
(1) K22
∂x ∂x · kf · dx ∂x ∂x
0
Universität Stuttgart
Institut für Wasserbau, Lehrstuhl für Hydromechanik und Hydrosystemmodellierung
home/lehre/VL-MHS-1-E/FOLIEN/VORLESUNG/3_FEM_INTRO/matr_entr_22a.tex
page 19 of 25
. – p.19/25
Calculating matrix entries VI (FEM – 1D) (1)
K22 =
Z1
1 · kf · 1 dx = kf · [x]10
0
(1) K22 (2)
K22 =
Z2 1
(2) K22
=
Z2 1
Universität Stuttgart
= kf [1 − 0] = kf (2)
(2)
∂N2 ∂N2 · kf · dx ∂x ∂x
∂(2 − x) ∂(2 − x) · kf · dx ∂x ∂x
Institut für Wasserbau, Lehrstuhl für Hydromechanik und Hydrosystemmodellierung
home/lehre/VL-MHS-1-E/FOLIEN/VORLESUNG/3_FEM_INTRO/matr_entr_22b.tex
page 20 of 25
. – p.20/25
Calculating matrix entries VII (FEM – 1D) (2)
K22 =
Z2
−1 · kf · −1 dx = kf · [x]21
1
(2)
K22 = kf [(2 − 1)] = kf (1)
(2)
A22 = K22 + K22 = 2kf All other matrix entries can be calculated analogously. A32 = A34 = A43 = A45 = A54 = −kf A33 = A44 = 2kf
A21 = A12 = −kf Universität Stuttgart
;
A11 = kf
Institut für Wasserbau, Lehrstuhl für Hydromechanik und Hydrosystemmodellierung
home/lehre/VL-MHS-1-E/FOLIEN/VORLESUNG/3_FEM_INTRO/matr_entr_22c.tex
page 21 of 25
. – p.21/25
Inserting values in matrix (FEM – 1D) Inserting these values one obtains the following system of equations: h2 h 1 kf 2 −1 0 0 −1 2 −1 0 h 0 3 kf · = 0 −1 2 −1 h4 0 h5 0 0 −1 1 −Q5
Q2 = Q3 = Q4 equal zero, as there are no sinks and sources in the domain. Q5 is the known Neumann boundary condition.
Universität Stuttgart
Institut für Wasserbau, Lehrstuhl für Hydromechanik und Hydrosystemmodellierung
home/lehre/VL-MHS-1-E/FOLIEN/VORLESUNG/3_FEM_INTRO/values_matrix.tex
page 22 of 25
. – p.22/25
Solving the set of equations (FEM – 1D) By transforming the vector representation one obtains a linear system of equations. Inserting the boundary conditions we have four equations for four unknowns. 2h2 − h3 = 50m
−h2 + 2h3 − h4 = 0 −h3 + 2h4 − h5 = 0
−h4 + h5 = kf has the unit
Universität Stuttgart
m , s
−4 m3 −10 s
kf · 1m
but there is also the unit m from the integration.
Institut für Wasserbau, Lehrstuhl für Hydromechanik und Hydrosystemmodellierung
home/lehre/VL-MHS-1-E/FOLIEN/VORLESUNG/3_FEM_INTRO/lgse.tex
page 23 of 25
. – p.23/25
Solving the set of equations (FEM – 1D) With kf = 10−5 and using the Thomas algorithm, we express our unknown in terms of one other unknown and the left boundary condition. 50 1 h2 = + h3 2 2 50 2 h3 = + h4 3 3 25 3 h4 = + h5 2 4 −4 10 25 h5 = 4 − −5 + . 10 2 Substituting back yields the following result h1 = 50m ; h2 = 40m ; h3 = 30m ; h4 = 20m ; h5 = 10m. Universität Stuttgart
Institut für Wasserbau, Lehrstuhl für Hydromechanik und Hydrosystemmodellierung
home/lehre/VL-MHS-1-E/FOLIEN/VORLESUNG/3_FEM_INTRO/lgse2.tex
page 24 of 25
. – p.24/25
Test of the solution (FEM – 1D) Test: Q1 = Q5 has to be fulfilled. h1 ·
(1) A11
+ h2 ·
(1) A12
50 · kf · 1m + 40 · (−kf ) · 1m = Q1 =⇒
Universität Stuttgart
= Q1 =⇒
Q1 = 10−4 m/s
3 3 m m = 10−4 10−4 s s
Institut für Wasserbau, Lehrstuhl für Hydromechanik und Hydrosystemmodellierung
home/lehre/VL-MHS-1-E/FOLIEN/VORLESUNG/3_FEM_INTRO/lgse2a.tex
page 25 of 25
. – p.25/25
