Date post: | 06-Feb-2018 |
Category: |
Documents |
Upload: | jonny-schnapsglas |
View: | 222 times |
Download: | 0 times |
7/21/2019 EFEM Kapitel 3-Begleitmaterial
http://slidepdf.com/reader/full/efem-kapitel-3-begleitmaterial 1/18
Einführung in die FEM Prof. Zehn SS 2014
75
3. Techniken in FEM-Programmen
3.1. Aufbau und Speicherung der Systemsteifigkeitsmatrix
3.1.1. Strategien zur Speicherung der Systemsteifigkeitsmatrix
Zwei einfache Beispiele zum Aufbau der Systemsteifigkeitsmatrix
(1) Modell Koinzidenzmatrix
Zur Vereinfachung je Knoten nur einen Freiheitsgrad, jedes Element gleich, hier
Elementsteifigkeitsmatrix für das Element 1
aij = a ji (symm.)
2 1 4 5
1 a11 a12 a13 a14 2
2 a22 a23 a24 1
3 a33 a34 4
4 a44 5
1 2 3 4
(1) (2) (3)
(4) (5) (6)
1
2
3
4
5
6
7
8
9
10
11
12
Knotennr.Elem.
1 2 3 4
1 2 1 4 52 5 4 7 8
3 8 7 10 11
4 3 2 5 6
5 6 5 8 9
6 9 8 11 12
GlobaleNumerierung
ElementlokaleNummerierung
7/21/2019 EFEM Kapitel 3-Begleitmaterial
http://slidepdf.com/reader/full/efem-kapitel-3-begleitmaterial 2/18
Einführung in die FEM Prof. Zehn SS 2014
76
Systemsteifigkeitsmatrix
1 2 3 4 5 6 7 8 9 10 11 12
1a22 a23 a24
2a11+a22+
a13 a14+a23
a24
3a11 a13 a14
4 a33+a22+
a34 a23 a24
5 a44+a33+a22+a11
a34 a13 a14+a23
a24
6a44+a11
a13 a14
7a33+a22
a34 a23 a24
8 a44+a33+a22+a11
a34 a13 a14+a23
a24
9a44+a11
a13 a14
10a33 a34
11a44+a34
a34
12a44
max. halbe Bandbreite =5, ( ) Knoten} jeadeFreiheitgr{1erenz}Knotendiff {maximale ⋅+=
ibw
Akt. Dreieck
7/21/2019 EFEM Kapitel 3-Begleitmaterial
http://slidepdf.com/reader/full/efem-kapitel-3-begleitmaterial 3/18
Einführung in die FEM Prof. Zehn SS 2014
77
Speichertechniken
1) Band, zeilenweise
(auch als eindimensionales Feld möglich,
spaltenweise bezogen auf die dargestellte
Form)
2) Spaltenweise, Skyline
3) Hypermatrixkonzept
4) Sparsematrixkonzept
zu 1) Modell
ibw = 9 bei einem Freiheitsgrad pro
Knoten (durch die Krümmung der Ränder
mit linearen Elementen eigentlich nicht
realisierbar)
(1) (2)
(3)(4)
1 23
4
5 6
7
8
a22 0 0 a23 a24 1.Zeile
2.
3.
4.
0
ibw –1=4
aktives Dreieck
12. Zeile
7/21/2019 EFEM Kapitel 3-Begleitmaterial
http://slidepdf.com/reader/full/efem-kapitel-3-begleitmaterial 4/18
Einführung in die FEM Prof. Zehn SS 2014
78
3.1.2. Bandweitenminimierung
⇒ Notwendigkeit aus dem Speicherschema erläutern
3.1.2.1 Verfahren von Cuthill–McKee
s. M.K.Schwarz „Methode der finiten Elemente“ s.171ff
Zur Bestimmung einer unteren und oberen Schranke für die Abschätzung der Güte der
erreichten Bandbreite.
untere Schranke
a priori: D maximaler Grad (Anzahl der an einem Knoten möglichen Kanten)
( )[ ]43421
121 +≥ Dibw
⇒ für ein Band ohne Nullelemente
obere Schranke
a posteriori: Aus dem Cuthill–McKee Algorithmus, aus der Anzahl der Knoten inden Stufen der Neunumerierung:
k → k+1 Nk → Nk+1
Im schlimmsten Fall kleinstindizierter Knoten im Niveau k mit größtindiziertem
Knoten im Niveau k+1 benachbart.
Die maximal mögliche Indexdifferenz Nk+Nk+1-1 über alle Stufen ist eine sichere
obere Schranke. Bei ν Stufen mit N0=1 als Startknoten gilt:
( )1max 1,...,1 −+≤ +
= k k
vk N N ibw
[...] ganzzahliger Teil
7/21/2019 EFEM Kapitel 3-Begleitmaterial
http://slidepdf.com/reader/full/efem-kapitel-3-begleitmaterial 5/18
Einführung in die FEM Prof. Zehn SS 2014
79
Cuthill – McKee Algorithmus (basiert auf der Graphentheorie)
1. Schritt: Man suche im Graphen G(S) einen Knoten mit minimalen Grad D als
Startknoten dieser bekommt die Nummer 1
2. Schritt: Zum Startknoten N0 bestimmt man alle benachbarten Knoten. Diese
Knoten werden mit zunehmenden Grad fortlaufend numeriert - Knoten
haben die Distanz 1 vom Startknoten
3. Schritt: Zu den Knoten der ersten Stufe mit aufsteigenden (neuen) Nummern
bestimmt man sukzessive ihre benachbarten Knoten und numeriert sie
jeweils mit zunehmenden Grad – Knoten haben die Distanz 2 vom
Startknoten
i. Schritt: Wie im 3. Schritt
Beispiel
A
B
C
M E
F D
H
I
G
O
K
L
N
P
Q
Steifigkeitsmatrix
Die Ausgangsnummerierung ist mit
Buchstaben gekennzeichnet
mögliche Startknoten:A,B,C,H,I,M,N,P,Q
7/21/2019 EFEM Kapitel 3-Begleitmaterial
http://slidepdf.com/reader/full/efem-kapitel-3-begleitmaterial 6/18
Einführung in die FEM Prof. Zehn SS 2014
80
Graph G(S)
5
2
1
127
63
8
13
4
9
10
11
14
15
16
1.Stufe
2.Stufe
3.Stufe
1
2
5
63
74
12
15
9
11
8
13
10
14
16
1.Stufe 2.Stu fe
3.Stufe
4.Stufe
Nummerierung und Stufen fürden Starknoten A
Nummerierung und Stufen fürden Startknoten C
7/21/2019 EFEM Kapitel 3-Begleitmaterial
http://slidepdf.com/reader/full/efem-kapitel-3-begleitmaterial 7/18
Einführung in die FEM Prof. Zehn SS 2014
81
Zum Ablauf des Algorithmus von Cuthill-McKee für den Startknoten A
Schritt Knoten Noch nicht nummerierte
Nachbarknoten
Grad Nummerk N
Anzahl der Knoten im
Schnitt
1 1 B
D
E
3
6
5
2
4
3
3
2 2
3
4
C
F
K
M
G
3
4
6
3
7
5
7
8
6
9
5
3 5
6
7
8
9
-
N
-
O
H
L
3
5
3
6
10
11
12
13
4
4 10
11
12
13
-
P
I
Q
3
3
3
14
15
16
3
Untere Schranke für unser Beispiel
7 D = 1 ( 1) 42
ganzzahliger Teil
ibw D ≥ + =
Obere Schranke für unser Beispiel
1max ( 1) 8
k=1,...,
k k ibw N N
ν
+≤ + − =
Die erreichte halbe Bandweite ist 6ibw = , damit ist die obere und untere Grenze
bestätigt
4 8ibw< <
7/21/2019 EFEM Kapitel 3-Begleitmaterial
http://slidepdf.com/reader/full/efem-kapitel-3-begleitmaterial 8/18
Einführung in die FEM Prof. Zehn SS 2014
82
3.2. Gleichungslösung – zentraler Prozess im FE –System
3.2.1. Direkte Löser
Nach Aufbau der Systemsteifigkeitsmatrix (Assemblierung) erfolgt das Lösen eines
großen Gleichungssystems
f vK =⋅
zur Bestimmung der globalen Knotenverschiebungen.
(a) EliminationsverfahrenBekannt aus der Mathematik ist sicher der Gaußsche–Algorithmus
K v f K v f
Symmetrie der Matrix nicht erforderlich, eine Bandstruktur geht verloren, die rechte
Seite wird verändert.
Durch Eliminaton (Addieren von Zeilen und Spalten, Zeilen- und Spaltentausch)
entsteht eine Rechtsdreiecksmatrix K ; durch Rückwärtseinsetzen erfolgt die
Berechnung des Vektors v .
Vorraussetzung für die Lösung: K ist positiv definit, da KK =T symmetrisch, muß die
quadratische Form 0>KvvT (nur für 0 }0,...,0{ ≡= Kvvv T T )
Eine indefinite Matrix K ist z.B. möglich, wenn 0≠v für f = 0 auftreten kann, d.h linear
abhängige Gleichungen im Gleichungssystem enthalten sind
⇒ z.B. Starrkörperverschiebungen möglich sind.
= =⇒ 0
7/21/2019 EFEM Kapitel 3-Begleitmaterial
http://slidepdf.com/reader/full/efem-kapitel-3-begleitmaterial 9/18
Einführung in die FEM Prof. Zehn SS 2014
83
Für die FEM – Anwendung häufiger Fehler:
K – nicht positiv definit
keine Gleichungslösung möglich
Ursachen:
1. fehlende Bedingungen (Randbedingungen),daher Starrkörperbewegungen
möglich, die zu v ≠ 0 bei f = 0 führen können,
2. fehlender Zusammenhalt in der Struktur, Mechanismus
3. fehlerhafter Elementaufbau (stark verzerrte Elemente, falsche Materialwerte)
Der Gaußsche–Algorithmus läßt sich auch als verketteter Gaußscher–Algorithmus
anwenden, indem eine Zerlegung der K Matrix erfolgt
KRL =⋅
Rechtsdreicksmatrix
Linksdreicksmatrix der Koeffizienten → (Diagonale 1)
Für eine symmetrische Matrix ergibt sich damit das
Cholesky – Verfahren
Zerlegung der Matrix K in
RTR=K
(durch negative Wurzel bei der Zerlegung
wird Indefinitheit angezeigt)
Unter Umständen wird die Indefinitheit der Matrix K durch Rundungsfehler bei der
nummerischen Lösung nicht bemerkt – Vorsicht bei der Ergebnisinterpretation!
R
2
1
1,5
1,414
-0,3536 1,62
0
1
3
3
4 2
5symm.
0
2 1
1,414
1,5
-0,3536
1,62
R
K
7/21/2019 EFEM Kapitel 3-Begleitmaterial
http://slidepdf.com/reader/full/efem-kapitel-3-begleitmaterial 10/18
Einführung in die FEM Prof. Zehn SS 2014
84
Es ist auch eine Zerlegung in der Form LTDL = K möglich (mit D=diag(d i)) und lii=1)
Zerlegung kann als eigenständiger Prozeß ausgeführt werden, wobei k ij kann durch r ij
überschrieben werden kann.
Matrix K ist eindeutig indefinit, wenn bei der Zerlegung eine negative Wurzel auftritt!
Nebenbemerkung:
( )
( ) ( ) ( ) ( )
( ) ( ) ( ) ( )2
1
2)det(detdetdetdet
...detdetdet...det
det
==⋅=⋅=
⋅⋅⋅=⋅⋅⋅
∏=
n
i
ii
T T r
C B AC B A
RRRRRK
K
Cholesky-Verfahren:
Kv = f = RTRv = f Rv =y
1. Vorwärtseinsetzen
RTy = f
2. Rückwärtseinsetzen
R v = y
In der Regel wird f →→→→ y und y →→→→ v im Speicher überschrieben ( Platzersparnis)
Vorzüge :
1. K wird auf Definitheit geprüft
2. Bandstruktur von K überträgt sich auch auf R
3. Speicher optimal, da k ij durch r ij überschrieben werden kann
4. Zusätzliche „rechte Seite“ in f auch nachträglich möglich, da Zerlegung
die „rechte Seite“ nicht betrifft (anders als beim allgem. Gaußschen–
Algorithmus)
5. Numerisch sehr stabil
=0
RT
y f
=0
R
v y
7/21/2019 EFEM Kapitel 3-Begleitmaterial
http://slidepdf.com/reader/full/efem-kapitel-3-begleitmaterial 11/18
Einführung in die FEM Prof. Zehn SS 2014
85
Cholesky – Verfahren ist für:
• Bandspeicherung
• Skylinespeicherung (nicht nur das Band sondern auch die Hülle bleibt
unverändert)
• Hypermatrixspeicherung (Blockbandverfahren)
mit speziellen Algorithmen und Speichertechniken einsetzbar. Nicht sinnvoll für
Kompaktspeicherung, da durch die Zerlegung Null- und Nichtnullelemente verändert
werden
→ aufwendige neue Speicherbelegung wäre erforderlich
Alternative: Frontlösungsmethode (Irons, Melosch, Beauford)
Kompilation und Elimination auf der Basis der Elemente
3.2.2. Fehler bei der Gleichlösung – Konditionszahl
Der numerische Fehler bei Lösung des Gleichungssystems lässt sich durch dieKonditionszahl abschätzen.
f vK =⋅
Konditionszahl von K ist ( ) 1−⋅= KKKκ ( beliebige Matrixnorm )
Deutung von ( )Kκ : Wenn mit einer festen Stellenzahl (Dezimalstellen) gerechnet wird,
dann ist im Ergebnis eine Abweichung von ( )Kκ Einheiten in den
letzten Dezimalstellen möglich.
Für symmetrische, positiv definite Matrizen gilt
( )min
max
λ
λ κ =K
λmax größter und λmin kleinster Eigenwert von K.
Unsichere Dezimalstellensind somit: ( )( )Kκ lg=a
7/21/2019 EFEM Kapitel 3-Begleitmaterial
http://slidepdf.com/reader/full/efem-kapitel-3-begleitmaterial 12/18
Einführung in die FEM Prof. Zehn SS 2014
86
Die größtmögliche Zunahme des Fehlers in v wird mit ( )Kκ angegeben.
• Für den direkten Löser → Dekompositionsverfahren
f vK =⋅ Näherungslösung ⇒ f vK ˆˆ =⋅
f f Kvv 1 ˆˆ −=− −
f f Kvv 1 ˆˆ −⋅≤− −
und
f vK ≥⋅
Daraus erhält man
1K
f v
≥
und damit erhält man für den relativen Fehler
ˆv v
v
−
die Gleichung
• Für iterative Löser gibt die Konditionszahl die Konvergenzgüte der Iteration an
Die Konditionszahl wächst mit der Feinheit der Diskretisierung
( ) ( )nhO2−
=Kκ
h – charakteristische Elementabmessung
n – Anzahl der Freiheitsgrade
Die feinere Diskretisierung ergibt auf der einen Seite eine bessere Approximation der
Näherungslösung erhöht allerdings den numerischen Fehler. Dies wird abgefangen
durch entsprechende Stellenzahlen (Wortlängen) im Rechner 32-Bit, 64-Bit,…(?).
( )f
f f
KKv
vv
K
ˆˆ1
−
⋅⋅≤
−−
43421
κ
7/21/2019 EFEM Kapitel 3-Begleitmaterial
http://slidepdf.com/reader/full/efem-kapitel-3-begleitmaterial 13/18
Einführung in die FEM Prof. Zehn SS 2014
87
3.2.3. Iterative Löser für K·v=f
Matrix K muss während der Gleichungslösung nicht verändert werden, damit ist auch
eine effiziente Kompaktspeicherung einsetzbar. Die Matrix müsste nicht zwingend
aufgebaut werden, man könnte auch mit den Elementmatrizen arbeiten.
Grundgedanke der Iteration
K+⋅= −
ia)()( 1ii vv iterative Verbesserung der Lösung in jedem Schritt
Große Verfahrensvielfalt
Vorteile :
• Zerlegung von K nicht erforderlich
• Optimale Speicherung realisierbar
• Nur Vektoroperationen erforderlich → (Skalarprodukte) → massive
Parallelisierung
• Für bestimmte Aufgabenklassen auch bei großen Dimensionen extrem
schnelle Lösung
Nachteile :
• Probleme mit der Stabilität ( schlechte Konvergenz bei großen
Konditionszahlen) und der Steuerung der Lösung (Abbruch)
• Keine gleichzeitige Bearbeitung von mehreren rechten Seiten oder schnelle
Abarbeitung neuer rechter Seiten ( wofür beim Cholesky-Verfahren nur
zusätzliches Vorwärts- und Rückwärtseinsetzen erforderlich ist)
⇒ Iteration immer nur für einen Lastfall
7/21/2019 EFEM Kapitel 3-Begleitmaterial
http://slidepdf.com/reader/full/efem-kapitel-3-begleitmaterial 14/18
Einführung in die FEM Prof. Zehn SS 2014
88
Es wird im Weiteren nun eine wichtige und weitverbreitete Verfahrensgruppe
vorgestellt:
Die Methode der konjugierten Gradienten (CG – Verfahren)
Funktionsanalytisch ist die Lösung des Gleichungssystems äquivalent mit der
Minimierung des Funktionals
( )
( )
T T12
T T0
00
F
F δ δ δ
δ
= ⋅ ⋅ − ⋅
= ⋅ − =
≠⋅ − =
v v K v f v
v K v f v
vK v f
Dieses Minimum des Funktionals wird nun iterativ bestimmt, indem die Richtung des
Gradienten der zu bestimmenden Funktion berechnet wird,
( ) rf Kvv =−=F grad
der dem Residuenvektor r zum Vektor v entspricht.
Der Gradient weist in Richtung der lokal stärksten Zunahme für F , deshalb wird mit dem
negativen Gradient zur Festlegung der Relaxationsrichtung gearbeitet, um das
Minimum zu erreichen.
Ziel: ( )( ) ( )( )ii vv F F <+1
Damit begründet sich auch der Name des Verfahrens: hier wird der konjugierte
Gradient benutzt wird.
Schema:
• Startvektor v(0)
• Relaxationsrichtung
( ) ( ) ( )( ) ( )( )1 0 0 0 grad F = − = − = − ⋅ −p r v K v f
• neuer Vektor
( ) ( ) ( )1
1
01 pvv q+=
(entspricht der Minimierung deselastischen Potentials)
konjugierter Gradient
7/21/2019 EFEM Kapitel 3-Begleitmaterial
http://slidepdf.com/reader/full/efem-kapitel-3-begleitmaterial 15/18
Einführung in die FEM Prof. Zehn SS 2014
89
Das neue Funktional ist damit
( ) ( ) ( )
( ) ( )
(1) (0) (1) (0) (1) (0) (1)
1 1 1
(1) (0) (0) (0) 2 (1) (1) (0) (1)
1 1 1
1( )
2
1
·
· · ·2
2
T T
T T T T
F v v q p q p v q p
q p K v K q p K p f v q
K v f
v pv
−= + + + =
= + + − +
Nebenrechnung:
( )
( ) ( )
(1) 2 (1) (1) (0) (1) (0) (1)
1 1
(1) (0) (1) (1) (0) (1) (0)
1 1 1
(0)
1( ) ( )
2
·
·
T
T
T T
T T T
F v q p K F v q p Kv f p
q p Kv f p q p Kv f q p r
r
p= + + −
− = − =
Die notwendige Bedingung für ein Minimum von F
( )( ) ( ) ( )( ) ( ) ( )( ) ( ) ( )( )( ) ( ) ( ) ( ) ( )( )
T1 0 1 0 1 0 1T1
1 1 12
T1 1 1 0 021
1 12 Min.
F q q q
q q F
= + ⋅ ⋅ + − +
= ⋅ ⋅ + ⋅ + →
v v p K v p f v p
p K p p r v
ist nun (Variation nach q1)
( ) ( ) ( ) ( )1
0T11T1
10 qqF δ δ rppKp +⋅⋅== 01 ≠qδ (es war (1) (0) p r = − )
( ) ( )
( ) ( )
( ) ( )
( ) ( )1T1
0T0
1T1
0T1
1
pKp
rr
pKp
rp
⋅⋅
⋅=
⋅⋅
⋅−=q
im i-ten Schritt wird die Relaxationsrichtung p(i) aus einer Linearkombination von
( ) ( )( ) ( )( )f vKvr −⋅−=−=− −−− 111 iii F grad
und ( )1−ip gebildet.
7/21/2019 EFEM Kapitel 3-Begleitmaterial
http://slidepdf.com/reader/full/efem-kapitel-3-begleitmaterial 16/18
Einführung in die FEM Prof. Zehn SS 2014
90
Algorithmus:
Start: Wahl von v(0)
( ) ( ) ( ) ( )
( ) ( ) ( )
( ) ( ) ( ) ( )
( ) ( ) ( )
( ) ( ) ( ) ( )( )
( ) ( )
( ) ( ) ( )( )
0 0 1 0
1 0 1
1
1 T 1 2 T 21
1 1
1
1 T 1 T
1 ( )
1
;
2, 3,
/ ; 2
/
i i i ii
i i i
i
i i i i
i
i i i
i i i
i
q
i
e i
e
q
q
q
− − − −
−
− −
−
− −
−
−
= ⋅ − = −
= +
=
= ⋅ ⋅ ≥
= − + ⋅
= ⋅ ⋅
= + ⋅
= + ⋅ ⋅
r K v f p r
v v p
r r r r
p r p
r r p K p
v v p
r r K p
K
Erklärung für die Operationen:
Abbruch, wenn ( )ε ≤
ir
( )( ) ( 1) ( )·i i i
ir r q K p−
= +
( 1) ( )i i
iKv f Kq p−
= − +
( )( 1) ( )i i
iK v q p f −
= + −
( ) ( )· i ivK f r = − =
Wichtig ist die richtige Wahl vonv(0) (je näher an der richtigenLösung desto besser)
Nur einmal zu berechnen
7/21/2019 EFEM Kapitel 3-Begleitmaterial
http://slidepdf.com/reader/full/efem-kapitel-3-begleitmaterial 17/18
Einführung in die FEM Prof. Zehn SS 2014
91
Theoretisch wird nach höchstens n Schritten eine Lösung erhalten.
Diese Zahl kann durch hohe Konditionszahlen von K numerisch
(große Konditionszahlen → schlechte Konvergenz !) bei iterativen Methoden höherliegen (bedingt durch endliche Zahlendarstellungen).
Für den Fehler in k-ten Schritt gilt:
( ) ( )k k e vv −=
und mit der Energienorm (potentielle Energie)
Tv v v= ⋅ ⋅
KK
die Fehlerabschätzung
( ) ( )
KK
0
1
12 ee
k
k
+
−≤
κ
κ ,
d.h. der Fehler verringert sich in k Schritten in Abhängigkeit von der Konditionszahl.
Kleine Konditionszahlen besitzen eine große Verringerung und große eine geringe.
Der für den Iterationsschritt k bereits dargestellten Fehler
( )
( )
k k
e
e
+
−≤
1
12
0κ
κ
K
K
ist abhängig von der Konditionszahl, d.h. die Anzahl der notwendigen Iterationen, um
eine Toleranz ε
( ) ( )ε ≤
KK
0 / ee k
zu erreichen, kann durch
12
ln21 +
≤
ε κ k
abgeschätzt werden.
7/21/2019 EFEM Kapitel 3-Begleitmaterial
http://slidepdf.com/reader/full/efem-kapitel-3-begleitmaterial 18/18
Einführung in die FEM Prof. Zehn SS 2014
92
Durch eine Vorkonditionierung in der Form
f CvKC ⋅=⋅⋅ −− 11
kann eine Verbesserung erreicht werden. Der Idealfall wäre C=K, da dafür die
Konditionszahl ( ) 1= Eκ
ist. Dies würde aber für die Bildung der Vorkondi-tionierungsmatrix den gleichen Aufwand wie zur vollständigen Lösung des
Gleichungssystems erfordern und somit unsinnig sein.
Deshalb werden für die Vorkonditionierung Näherungen von C-1 konstruiert. Im
Algorithmus würde sich mit einer Vorkonditionierungsmatrix nur
( ) ( )
( ) ( )
( ) ( )1T1
01T0
1
011
pKp
rCr
rCp
⋅⋅
⋅⋅=
⋅=
−
−
q
und
( ) ( ) ( ) ( )
( ) ( ) ( )
( ) ( ) ( ) ( )( )iiii
i
i
i
ii
iiii
i
q
e
e
pKprCr
prCp
rCrrCr
⋅⋅⋅⋅=
+⋅−=
⋅⋅⋅⋅=
−−−
−
−
−−
−−−−−−
−
/
/
11T1
1
1
11
21T211T1
1
ändern.
Die praktisch genutzten Verfahren benutzen eine solche Vorkonditionierung
→ PCG-Verfahren
(Preconditioned Conjugate Gradient)