Notesale: Turn your study into money

Document Preview

Extracts from the notes are below, to see the PDF you'll receive please use the links above

---

Skript

Finite{Elemente{Methode
Eine Einfuhrung fur Ingenieurstudenten
Michael Jung
Technische Universitat Chemnitz{Zwickau
Ulrich Langer
Johannes{Kepler Universitat Linz

Autorenadresse:
Dr
...
nat
...
Dr
...
nat
...
Ulrich Langer
(O
...
{Professor)
Johannes{Kepler Universitat Linz
Technisch{Naturwissenschaftliche Fakultat
Institut fur Mathematik
Ordinariat "Numerische Mathematik\
Altenberger Stra e 69
A { 4040 Linz / Osterreich

Vorwort
Das vorliegende Skript entstand auf der Grundlage einer Vorlesung zum Thema "Numerik partieller Di erentialgleichungen\, welche die Autoren in den letzten Jahren fur
Ingenieurstudenten gehalten haben
...
Ein derartiges Skript kann naturlich nicht eine umfassende Darlegung
aller bekannten Diskretisierungsmethoden beinhalten
...
Unser Ziel besteht darin, anhand
von Beispielen den gesamten Proze vom Aufstellen des mathematischen Modells einer
physikalischen Erscheinung bis zur Computerrealisierung zu erlautern
...
Mit der Auswahl
von Problemen aus der Thermodynamik, der Elektrostatik und der Festkorpermechanik soll dem Leser verdeutlicht werden, da in vielen Anwendungsgebieten partielle
Differentialgleichungen bei der Modellierung entstehen
...
Die hergeleiteten Di erentialgleichungen dienen in den folgenden Kapiteln als Modellbeispiele bei
der Beschreibung der Finite{Elemente{Methode
...

Dabei beschreiben wir ausgehend von der Warmeleitgleichung in di erentieller Form
alle Schritte, die bei einer Finite{Elemente{Diskretisierung erforderlich sind, d
...
den
Ubergang zur verallgemeinerten Formulierung, die Diskretisierung des Gebietes, die
Wahl der Ansatzfunktionen, den Aufbau des Finite{Elemente{Gleichungssystems, die
Losung dieser Gleichungssysteme und Fehlerabschatzungen
...
Zunachst beschreiben wir allgemein das Ritz{ und
das Galerkin{Verfahren als mogliche Diskretisierungsverfahren
...
Wir
beschranken uns im wesentlichen auf Diskretisierungen mit linearen Dreieckselementen
...

Da jede Finite{Elemente{Diskretisierung letztendlich auf ein im allgemeinen gro dimensioniertes lineares Gleichungssystem fuhrt, haben wir einen Schwerpunkt auf die

4
Diskussion verschiedener Au osungsmethoden gelegt
...
Wir beschreiben die Algorithmen klassischer direkter Verfahren, klassischer iterativer Verfahren sowie moderner iterativer Verfahren und zeigen ihre Vorteile
und Nachteile auf
...

Wir mochten uns bei unseren Kollegen Herrn Dr
...
Heise, Herrn Dr
...
Queck, Frau
Dr
...
Weber, Herrn Dipl
...
T
...
{Math
...
Vogel fur
die Unterstutzung bei der Aufbereitung von Beispielen und fur zahlreiche Hinweise
bezuglich der Gestaltung des Skripts bedanken
...
Art
Randstuck mit Randbedingungen 2
...
Art
abgeschlossenes Randstuck ;i, i = 1 2 3
Triangularisierung des Gebietes
ein Element (Dreieck) der Triangularisierung
Referenzdreieck
Anzahl der Knoten in der Vernetzung Th
Anzahl der Knoten, die in
;2 ;3 liegen
Indexmenge, welche die Nummern aller Knoten enthalt
Indexmenge, welche die Nummern der Knoten in
;2 ;3 enthalt
Indexmenge, welche die Nummern der Knoten auf ;1 enthalt
Raum der Vektoren mit n Komponenten
Menge der im Intervall a b] stetigen Funktionen
Menge der im Intervall (a b) stetigen Funktionen
Menge der im Intervall a b] k{mal stetig di erenzierbaren Funktionen
Menge der im Intervall (a b) k{mal stetig di erenzierbaren Funktionen
Menge der in stetigen Funktionen
Menge der in stetigen Funktionen
Menge der in k{mal stetig di erenzierbaren Funktionen
Menge der in k{mal stetig di erenzierbaren Funktionen
Menge der uber dem Intervall (a b) quadratisch integrierbaren Funktionen
Menge der uber quadratisch integrierbaren Funktionen
Menge der L2{Funktionen, deren erste verallgemeinerte Ableitung existiert und ebenfalls Element des Raumes L2(a b) ist
Menge der L2{Funktionen, deren erste verallgemeinerte Ableitungen existieren und ebenfalls Element des Raumes L2( ) sind

6

V
Vg

Menge der Grundfunktionen
Menge der Funktionen aus V , welche die Randbedingungen 1
...
Art (naherungsweise) erfullen
V0h
Menge der Funktionen aus Vh , welche auf dem Randstuck ;1 gleich Null
sind
a(: :) Bilinearform
hF :i Linearform (rechte Seite)
v
ein Element eines Funktionenraumes
ein Element des Rn
v
u
exakte Losung eines Randwertproblems
uh
FE{Naherungslosung
@v
@xi (verallgemeinerte) Ableitung der Funktion v nach xi
@v
@N Konormalenableitung von v
v
Laplace{Operator angewendet auf eine Funktion v
grad v Gradient von v
div ~ Divergenz einer Vektorfunktion ~
v
v
rot ~ Rotor einer Vektorfunktion ~
v
v
Kh
Stei gkeitsmatrix
fh
Lastvektor
Losung des FE{Gleichungssystems
uh

Inhaltsverzeichnis
Vorwort

3

Liste verwendeter Bezeichnungen

5

Inhaltsverzeichnis

7

1 Einfuhrung

9

1
...
2 Zur Geschichte der FEM : : : : : : : : : : : : : : : :
1
...
3
...
3
...
3
...
3
...
3
...
1 Stationare Warmeleitprobleme : : : : : : : : : : : : : : : : : : : : : :
2
...
1 Das stationare 1D{Warmeleitproblem : : : : : : : : : : : : : :
2
...
2 Das stationare 2D{Warmeleitproblem : : : : : : : : : : : : : :
2
...
3 Das stationare 3D{Warmeleitproblem und einige Spezialfalle :
2
...
2
...
2
...
1
3
...
3
3
...
5
3
...
7
3
...
9

H 1(a

Die Funktionenraume L2(a b) und
b) : : : : : : : : : : : : :
Variationsformulierung von Randwertaufgaben : : : : : : : : : : :
Die FEM zur naherungsweisen Losung des Variationsproblems : : :
Der elementweise Aufbau der Stei gkeitsmatrix und des Lastvektors
1D{LAGRANGE{Elemente hoherer Ordnung : : : : : : : : : : : :
Au osung des Finite{Elemente{Gleichungssystems : : : : : : : : :
Diskretisierungsfehlerabschatzungen : : : : : : : : : : : : : : : : :
Zwei Anwendungsbeispiele : : : : : : : : : : : : : : : : : : : : : : :
Das Programm FEM1D : : : : : : : : : : : : : : : : : : : : : : : :
7

11
12
15
15
19
21
24
29

33
33
33
38
42
46
46
49

51

: : 51
: : 53
: : 55
: 61
: : 67
: : 71
: : 74
: : 83
: : 94

8

INHALTSVERZEICHNIS

4 FEM fur mehrdimensionale Randwertprobleme 2
...
1 Modellprobleme : : : : : : : : : : : : : : : : : : : : : : : : : : : : : :
4
...
3 Galerkin{Ritz{FEM : : : : : : : : : : : : : : : : : : : : : : : : : : : :
4
...
1 Galerkin{Verfahren : : : : : : : : : : : : : : : : : : : : : : : :
4
...
2 Ritz{Verfahren : : : : : : : : : : : : : : : : : : : : : : : : : : :
4
...
3 FEM = Ritz{Galerkin{Verfahren mit speziellen Ansatzfunktionen : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : :
4
...
4
...
4
...
4
...
4
...
4
...
4
...
1 Direkte Verfahren : : : : : : : : : : : : : : : : : : : : : : : : : : : : :
5
...
2
...
2
...
2
...
2
...
2
...
3 Ein Vergleich der Au osungsverfahren : : : : : : : : : : : : : : : : : :

99

99
101
104
104
106

107
111
111
118
121
135
136
142

147
148
154
155
156
158
161
166
168

6 Galerkin{FEM fur parabolische Anfangs{Randwertaufgaben

177

Literaturverzeichnis

183

6
...
2 Konvergenz und Stabilitat : : : : : : : : : : : : : : : : : : : : : : : : : 180

Kapitel 1

Einfuhrung
Viele technische und physikalische Prozesse, wie z
...
die Warmeleitung in festen Korpern, die Deformation von Bauteilen unter vorgegebenen Belastungen sowie elektrische
und magnetische Felder konnen durch partielle Di erentialgleichungen beschrieben werden
...
h
...
Der Einsatz von Computern erfordert den Ubergang
vom kontinuierlichen Problem (Di erentialgleichung) zu einem endlichdimensionalen
Ersatzproblem
...
Diskretisierungsmethoden fur partielle Di erentialgleichungen sind die
FEM { Finite{Elemente{Methode (Finite Element Method) 27, 92],
FDM { Di erenzenverfahren (Finite Di erence Method) 30, 70],
BEM { Randelemente{Methode (Boundary Element Method) 6, 35]
...
Ein Vorteil dieser Methode besteht
darin, da sie bei der Diskretisierung von Problemen in beliebigen beschrankten zwei{
und dreidimensionalen Gebieten erfolgreich angewendet werden kann
...
B
...
Die Randelemente{Methode ist besonders geeignet zur Losung von Feldproblemen in unbeschrankten Gebieten wie z
...
zur Berechnung der Fernwirkung von elektrischen und magnetischen Feldern
...
Ordnung beschaftigen
...
Dabei bezeichnen

@2 @2 @2
= @x2 + @x2 + @x2
1

2

3

10

KAPITEL 1
...
Die Funktionen f (x1 x2 x3)
sowie g(x1 x2 x3) sind vorgegeben
...
B
...

Allgemein bedeutet die Losung einer Randwertaufgabe 2
...
Ordnung in einem bestimmten Gebiet,
wobei die gesuchte Losung am Rand des Gebietes noch gewissen Bedingungen, den
Randbedingungen, unterliegt
...
Ordnung hangt
die gesuchte Funktion u(x1 x2 : : : xn) von mehreren Veranderlichen x1 x2 : : : xn ab,
und die hochste Ordnung der auftretenden partiellen Ableitungen ist 2
...
Dabei bezeichnen
die Funktionen f (x1 x2 x3 t) eine zeitlich veranderliche Warmequelle, g(x1 x2 x3 t)
die vorgegebene Randtemperatur und g0(x1 x2 x3) das Temperaturfeld zum Zeitpunkt
t = t0
...
Ordnung bedeutet allgemein die Bestimmung der Losung einer partiellen Di erentialgleichung 2
...

Weitere Beispiele fur Randwertaufgaben und Anfangs{Randwertaufgaben beschreiben
wir im Abschnitt 1
...


11

1
...
VOM TECHNISCHEN PROZESS BIS ZUR COMPUTERSIMULATION

1
...
Das folgende Schema charakterisiert alle Schritte, die notwendig
sind, um ausgehend von einer physikalischen Erscheinung oder einem technischen Proze ein adaquates mathematisches Modell aufzustellen
...


Physikalische Erscheinung
technischer Proze

Techniker,
Ingenieure,
Physiker

?

Physikalisch{technisches Modell
(Erhaltungssatze, Sto gesetze u
...
)
?

Mathematisches Modell

Mathematiker

Di erentialgleichungen (RWA, ARWA,
...

?

?

Funktionalanalytische Untersuchung des
mathematischen Modells

Frage nach der Existenz und Eindeutigkeit der Losung,
Ermittlung weiterer Eigenschaften der Losung
?

Numerische Analysis

1
...
a
...
Klassische numerische Untersuchungen
Frage nach der Existenz und Eindeutigkeit der
Naherungslosung, Konvergenzuntersuchungen,
Fehlerabschatzungen u
...

3
...
Algorithmisierung
?

Mathematiker,
Informatiker

Erstellung und Wartung von Software
?

Numerische Tests

-

Praktiker

12

KAPITEL 1
...
2 Zur Geschichte der FEM
Vorgeschichte
In der Arbeit 73] beschreibt Schellbach die Losung eines Minimal achenproblems
...
Der von Schellbach beschriebene Losungsweg beinhaltet Teilschritte, wie sie fur die Finite{Elemente{Methode charakteristisch
sind
...
3
...

Weiterhin zahlen die Arbeiten von Ritz und Galerkin zur Vorgeschichte der Finite{
Elemente{Methode
...
Galerkin{Verfahren
(siehe auch die Abschnitte 4
...
3)
...
W
...
B
...
1)
0

minimiert
...
2)
i=1

mit beliebigen reellen Zahlen ui in das Funktional J (u) ein
...
Das Minimum des Funktionals uber
allen Funktionen, die sich gema der Beziehung (1
...
Die Gute der Naherung u(n) fur die Funktion u hangt von n ab
...
1) kann man
beispielsweise
'i(x) = (1 ; x)xi
(1
...

2
...
B
...
2
...
Zur Bestimmung einer Naherungslosung wird wieder ein Ansatz der Gestalt (1
...
Es wird die Funktion u(n) gesucht, welche die Beziehungen
Z1
;(u(n))00 (x) ; f (x)] 'i (x) dx = 0
0

bzw
...


Geschichte

1
...
Courant (1943, siehe auch 14]): "Variational methods for the solution of
problems of equilibrium and vibrations\
Erstmals wurden im Ritzschen Verfahren Ansatzfunktionen mit lokalem Trager (sogenannte Hutchen{Funktionen, siehe Abschnitt 4
...
3) verwendet und FEM{Rechnungen durchgefuhrt, d
...
es wurden im Unterschied zu den Ansatzfunktionen (1
...

Die Arbeit von Courant fand zunachst keine besondere Beachtung
...
Neuentdeckung der FEM durch Mechaniker in den 50er Jahren:
Im Jahr 1956 beschrieb Turner eine Methode, die charakteristische Grundgedanken der FEM beinhaltete
...
h
...

3
...

Die ersten mathematisch fundierten Untersuchungen stammen von K
...
Friedrichs (1962, 21]), L
...
Oganesjan (1966, 62]), V
...
Korneev (1967, 45]),
M
...
a
...
Im Jahr 1967 wurde die erste Ingenieur{Monographie (Mechanik) zur FEM publiziert 90]
O
...
Zienkiewicz

"The Finite Element Method in Structural and Continuum Mechanics\
...

5
...
Strang und G
...
Fix
"An Analysis of the Finite Element Method\
...
EINFUHRUNG

Gegenwart
Die Finite{Elemente{Methode ist die standarde Diskretisierungsmethode fur Feldprobleme
...
Es existiert eine Vielzahl von FEM{Programmsystemen zur Losung
von Problemen der Festkorpermechanik, der Stromungsmechanik, der Elektrotechnik
u
...


Zukunft
Auf der Basis von FEM{Modellen und unter Nutzung modernster Rechentechnik (Parallelrechner, Vektorrechner) wird die Simulation komplexer Prozesse moglich (z
...
inverse Aufgabenstellungen)

Literatur
Ingenieur{FEM{Literatur:
{ Dankert (1977, 16])
"Numerische Methoden der Mechanik\
{ H
...
{G
...
Tobiska (1985, 27])
"Finite{Element{Methode\
{ N
...
Strang und G
...
Fix (1973, 77])
"An Analysis of the Finite Element Method\
{ V
...
Korneev (1977, 46])
"Schemy metoda konecnych elementov vysokich porjadkov tocnosti\
{ Ph
...
Weitere wichtige Literaturquellen
sind z
...
die Bucher von Axelsson/Barker 4], Braess 9], Breitschuh/Jurisch
11], Brenner/Scott 12], Fischer 20], Gallagher 22], Girault/Raviart 24],
Griffiths 29], Grossmann/Roos 30], Hackbusch 34], Hinton/Owen 41],
Kr zek/Neitaanmakki 48], Mercier 50], Mitchell/Wait 51], Hlavacek/
Necas 52], Norrie/de Vries 53], Oden 54], Oden/Reddy 61],Oden/Becker/
Carey 55, 56, 57, 58, 60, 59], Oganesjan/Ruchovec 63], Schwarz 74, 75],
Shajdurov 82], Whiteman 85], Zen sek 84], Zienkiewicz 91, 92, 93], Zienkiewicz/Morgan 95], Zienkiewicz/Taylor 96] u
...

Die Anzahl der Zeitschriften Artikel und Proceedings ist fast unuberschaubar
...


1
...
BEISPIELE

15

1
...
3
...
Die Temperatur u(x t) genugt der Di erentialgleichung
@
@u
@
@u
@
@u
c% @u ; @x @x ; @x @x ; @x @x = f
(1
...

Dabei bezeichnen c(x t) die spezi sche Warmekapazitat, %(x t) die Dichte, (x t) die
Warmeleitzahl und f (x t) die Intensitat der Warmequelle
...
5)
1
2
3
sowie
@v @v @v T
grad v = @x @x @x
(1
...
4) in der kompakteren Form
(1
...

Die Temperatur zum Zeitpunkt t = ta wird durch die Vorgabe einer Anfangstemperaturverteilung g0(x), der Anfangsbedingung
u(x ta) = g0(x) fur alle x 2
im mathematischen Modell berucksichtigt
...
Wir unterscheiden die drei folgenden Moglichkeiten:
1
...
h
...

2
...
h
...
Der Vektor ~ (x) = (n1(x) n2(x) n3(x))T ist der Vektor der au eren
n
Einheitsnormalen im Punkt x auf dem Rand ;
...
Mit g2 (x t) 0 kann man eine Warmeisolation modellieren
...
EINFUHRUNG

3
...
h
...

Naturlich ist es auch moglich, bei der Formulierung des Warmeleitproblems eine Kombination der obenbeschriebenen Randbedingungen zu nutzen, d
...
da auf einem Randstuck ;1 eine Temperaturverteilung und auf einem Randstuck ;2 ein Warmestrom vorgegeben werden sowie auf einem Randstuck ;3 der Warmeaustausch mit der Umgebung
modelliert wird (; = ;1 ;2 ;3 ;i \ ;j = fur i 6= j )
...
4) herleiten (siehe auch 17, 37])
...
Zuerst
betrachten wir das folgende stationare Problem:
Gesucht ist die Temperaturverteilung in einem Kuhlkorper
...
1 dargestellt
...
1: Schnitt durch den Kuhlkorper parallel zur (x1 x2){Ebene
Am Boden (Rand ;00) ie t der Warmestrom q(x) in den Kuhlkorper hinein, und uber
2
seine gro e Ober ache wird die Warme an die Umgebung abgegeben
...
3
...
Da der Querschnitt und die Eingangsdaten
wie Umgebungstemperatur, Warmeubergangszahl, Warmeleitzahl und Warmestrom
auch symmetrisch bezuglich der x2{Achse sind, konnen wir schlie lich unsere Berechnungen auf die Bestimmung des Temperaturfeldes in dem in der Abbildung 1
...
Es ist somit die folgende Randwertaufgabe zu losen:
Gesucht ist das Temperaturfeld u(x1 x2), fur das
@ @u ; @ @u = 0
;
fur alle x 2
@x1 @x1 @x2 @x2
@u := @u n + @u n = 0
fur alle x 2 ;02
@N
@x1 1 @x2 2
@u = q(x)
fur alle x 2 ;00
2
@N
@u + (x)u(x) = (x)u (x) fur alle x 2 ;
A
3
@N
gilt
...

Die Randbedingung auf dem Rand ;02 ist eine "kunstlich\ eingefuhrte Randbedingung
...

2

Die Abbildung 1
...
h
...


Abbildung 1
...
EINFUHRUNG

Als Beispiel fur ein instationares Warmeleitproblem betrachten wir den Abkuhlungsproze in einer abgesetzten Welle (siehe auch 8])
...
Wir wollen das zeitlich veranderliche Temperaturfeld u(x t) berechnen, d
...
wir bestimmen zu verschiedenen Zeitpunkten t > 0 die Temperaturverteilung
in der Welle
...
Die Meridianebene, d
...
die Flache, welche
durch Rotation um die z{Achse den Korper (Welle) erzeugt, ist in der Abbildung 1
...

z6

l

-

ra r

c =
670
% = 7:68 103
=
29
=
40
uA = 293:15
l =
0:8
ra =
0:35

Ws (kg K );1
kg m;3
W (m K );1
W (m2 K );1
K
m
m

;2 = f(r z) 2 @ : z = 0g f(r z) 2 @ : r = 0g

;3 = @ n ;2

Abbildung 1
...
4) lautet in Zylinderkoordinaten:
Gesucht ist das Temperaturfeld u(r ' z t), das der Di erentialgleichung

@u ; 1 @ r @u + 1 @ 2u + @ 2u = f
(1
...
Da das Gebiet, d
...
die Welle, rotationssymmetrisch ist sowie alle Eingangsdaten
wie spezi sche Warmekapazitat, Dichte, Warmeleitzahl, Warmeubergangszahl und
Umgebungstemperatur konstant sind, ist das Temperaturfeld unabhangig vom Winkel '
...
8) alle Ableitungen nach ' identisch Null
...
3), dann ergibt sich die folgende Randwertaufgabe:

19

1
...
BEISPIELE

Gesucht ist das Temperaturfeld u(r z t), fur das

@u ; a2 1 @ r @u + @ 2u = 0
fur alle (r z) 2
@t
r @r @r
@z2
u(r z 0) = 800 fur alle (r z) 2
@u = 0
@N

und fur alle t 2 (0 te)

fur alle (r z) 2 ;2 und fur alle t 2 0 te]

@u + u = u fur alle (r z) 2 ; und fur alle t 2 0 t ]
A
3
e
@N

gilt, mit

@u = a2 @u n + a2 @u n :
@N
@r r
@z z
In der Abbildung 1
...
Die Temperatur am au eren Rand betragt zu diesem Zeitpunkt 460o C ,
und in der Mitte hat die Welle noch eine Temperatur von 739o C
...
3
...
9)

und

~
div D = %
(1
...
Der Di erentialoperator div ist
durch die Gleichung (1
...
11)
rot ~ = @x3 ; @x2 @x1 ; @x3 @x2 ; @x1
v @v @v @v @v @v @v
2
3
3
1
1
2
mit ~ = (v1 v2 v3)T de niert
...
9) existiert eine skalare Funktion ', fur die
~
E = ; grad '

(1
...
Diese Funktion wird auch als elektrisches Potential bezeichnet
...
10) und (1
...
13)
(1
...
EINFUHRUNG

so kann die Di erentialgleichung (1
...
15)
"
geschrieben werden
...
Dabei reduzieren wir die Berechnungen wieder auf die Losung
eines Randwertproblems in einem zweidimensionalen Gebiet
...
4 ist
ein Schnitt in der (x1 x2){Ebene dargestellt
...
Die durch markierten Gebiete sind Abschirmplatten
...
h
...
Da wir zur numerischen Bestimmung des Potentials
die Methode der niten Elemente nutzen wollen, mussen wir unsere Untersuchungen
auf ein beschranktes Gebiet reduzieren
...
Auf diesem au eren Rand wird das elektrische Potential
als identisch Null angenommen
...
4: Aquipotentiallinien des elektrischen Potentials
In unserem Beispiel ist % = 0
...
Gesucht wird das elektrische Potential im Gebiet
( )
...
3
...

Die Abbildung 1
...


1
...
3 Berechnung elektromagnetischer Felder

Den Ausgangspunkt zur Aufstellung des mathematischen Modells bildet die Maxwellsche Gleichung

~ ~
rot H = S
(1
...
Weiterhin fuhren wir das Vektorpotential A ein, das durch die Gleichungen
und

~ ~
rot A = B

(1
...
18)
~
de niert ist
...
Der Zusammenhang zwischen
~
~
der magnetischen Induktion B und der magnetischen Feldstarke H ist durch
~
~ ~
B = 0 r (H + H0)
(1
...
Fur die nicht permanentmagnetischen Materialien wird H0 = ~ voraus~ 0 die magnetische Feldstarke, bei der
gesetzt
...
Besteht fur den Permanentmagneten ein
~
~
linearer Zusammenhang zwischen B und H ,
so ist

~
B=

~
0 rH

~
+ J0

(1
...
21)
~
wobei J0 die Permanentmagnetisierung bezeichnet
...


r

~
= r (jBj)

(1
...
EINFUHRUNG

Damit ergibt sich aus den Beziehungen (1
...
22) (siehe auch 39]) die Gleichung
1
~ ~
~
rot
(1
...
Die Erregung erfolgt durch radial magnetisierte Permanentmagnete,
und die Stromdichten in den einzelnen Nuten sind jeweils vorgegeben (siehe Abbildung 1
...
Au erdem sind die Magnetisierungskennlinien fur die ferromagnetischen
~ ~
Materialien Dynamoblech und Walzstahl fur gewisse Wertepaare (jH j jB j) tabellarisch
gegeben (siehe Tabelle 1
...

Dynamoblech
Walzstahl
~
~
~
~
jH j Am;1] jB j V s m;2 ] jH j Am;1] jB j V s m;2 ]
200
0
...
110
1000
1
...
0
2000
1
...
618
2500
1
...
750
5000
1
...
892
10000
1
...
062
20000
1
...
139
Tabelle 1
...
Der Querschnitt der Gleichstrommaschine ist in der Abbildung 1
...

~
Wir setzen voraus, da die x1{ und x2{Komponenten des Potentials A und der Strom~ identisch Null sind, d
...

dichte S

~
A = (0 0 A3)T

~
S = (0 0 S3)T :

und

Somit erhalten wir aus dem Di erentialgleichungssystem (1
...
3
...
5: Querschnitt des Gleichstrommotors
Damit ergibt sich das Randwertproblem:
~
Gesucht ist die x3{Komponente u = A3 des Vektorpotentials A, fur die
1
grad u = S3 ; @H01 + @H02 fur alle x 2
@x2 @x1
0 r (jgraduj)
gilt
...
24)

24

KAPITEL 1
...
6: Feldlinienbild fur die Gleichstrommaschine

1
...
4 Berechnung mechanischer Felder

Das Verschiebungsfeld ~ (x) = (u1(x) u2(x) u3(x))T unter vorgegebenen Belastungen
u
wird in der linearen Elastizitatstheorie als Losung eines Di erentialgleichungssystems
bestimmt (siehe 25, 26, 47])
...
Bezeichnen wir
~
mit f = (f1 f2 f3)T den Vektor der Volumenkrafte, ~1 = (g11 g12 g13)T den Vektor
g
der vorgegebenen Randverschiebungen, ~2 = (g21 g22 g23)T den Vektor der vorgegebeg
nen Ober achenkrafte und ~ = (n1 n2 n3)T den Vektor der au eren Einheitsnormalen,
n
dann kann die Randwertaufgabe der linearen Elastizitatstheorie in kartesischen Koordinaten folgenderma en geschrieben werden:
Gesucht ist das Verschiebungsfeld ~ (x), fur das
u
;

e

;

e

;

e

@ 2 u1 @ 2 u1 @ 2 u1
+ @x2 + @x2
@x2
1
2
3
2u
2u
2u
@ 2 @ 2 @ 2
+ @x2 + @x2
@x2
1
2
3
@2u

3

@x2
1

@
; ( e + e) @x

1

@
; ( e + e) @x

2

@u1 @u2 @u3
+
+
@x1 @x2 @x3
@u1 @u2 @u3
+
+
@x1 @x2 @x3

= f1 (x)

fur alle x 2

= f2 (x)

fur alle x 2

@u
@
+ @x23 + @x23 ; ( e + e) @x @u1 + @x2 + @u3 = f3 (x)
@x1
@x3
3
2
2
3
u
~ (x) = ~1 (x)
g
i1 n1 + i2 n2 + i3 n3 = g2i (x)
@2u

@2u

(1
...
3
...
Dabei bezeichnen

die Lameschen Elastizitatskonstanten, die durch
E
E
und
e=
e=
(1 + )(1 ; 2 )
2(1 + )
mit dem Elastizitatsmodul E und der Poissonschen Querkontraktionszahl de niert
sind
...
Fur die
Komponenten ij des Spannungstensors gilt
i=1 2 3
ii = e ("11 + "22 + "33) + 2 e "ii
i j = 1 2 3 i 6= j
ij = ji = 2 e "ij
mit den Verzerrungen
1 @u
"ij = 2 @xi + @uj
i j = 1 2 3:
(1
...
5) bzw
...
6) sowie
den Laplace{Operator eines Vektorfeldes ~ = (v1 v2 v3)T , d
...

v
e

und

e

~ = ( v1 v2 v3)T
v
(1
...


;

e

~ ; ( e + e) grad div ~ = f
u
u ~

(1
...
7: Querschnitt des Tragers in der (x1 x2){Ebene und die Deformation
der Kontur (stark vergro ert dargestellt)
Zuerst betrachten wir das folgende Beispiel
...
In der Abbildung 1
...


26

KAPITEL 1
...
Au erdem wirkt die Kraft F in Ebenen
senkrecht zur x3{Achse, und sie ist unabhangig von der x3{Richtung
...
29)
@x3 @x3 3
gilt und folglich ist
"13 = "23 = "33 = 0 :
Unter Beachtung der Beziehung (1
...
25)
das folgende ebene lineare Elastizitatsproblem (ebener Verzerrungszustand):
Gesucht ist das Verschiebungsfeld ~ (x) = (u1(x) u2(x))T , fur das
u
;

e

;

e

@ 2u1 + @ 2u1 ; ( + ) @
e
e
@x2 @x2
@x1
1
2
@ 2u2 + @ 2u2 ; ( + ) @
e
e
@x2 @x2
@x2
1
2

@u1 + @u2
@x1 @x2
@u1 + @u2
@x1 @x2

= 0

fur alle x 2

= 0

fur alle x 2

~ (x) = ~
u
0
11n1 + 12n2 = 0
21n1 + 22n2 = g22

(1
...
In unserem Beispiel ist E = 1:96 1011 Nm;2 , = 0:3 und
8 F = ; 106 x ; 4 105 Nm;2 auf der Oberseite des Tragers
>
1
<
(x1 2 0:0 0:18] x2 = 0:24)
g22 = >
:
0 Nm;2 sonst:
Als zweites Beispiel betrachten wir ein dickwandiges Rohr unter Innendruck
...
8 dargestellt
...

In diesem Koordinatensystem sind die Komponenten des Verzerrungstensors durch die
Beziehungen
"zz = @uz
"rr = @ur "'' = 1 @u' + ur
@r
r @'
@z
(1
...
3
...
32)

Es ist die folgende Randwertaufgabe zu losen:
Gesucht ist das Verschiebungsfeld ~ (r ' z) = (ur(r ' z) u'(r ' z) uz (r ' z))T , fur
u
das
2
2
2
r
; e @@ru2r + 1 @ur + r12 @ u2r + @@zu2r ; u2 ; r22 @u'
r @r
@'
r
@'
@
;( e + e) @r @ur + ur + 1 @u' + @uz = fr fur alle (r ' z) 2
@r r r @' @z
;

e

;

@ 2u' + 1 @u' + 1 @ 2u' + @ 2u' ; u' + 2 @ur
@r2 r @r r2 @'2 @z 2 r2 r2 @'
@
;( e + e ) 1 @' @ur + ur + 1 @u' + @uz
r
@r r r @' @z

= f'

fur alle (r ' z ) 2

@ 2uz + 1 @uz + 1 @ 2uz + @ 2uz
@r2 r @r r2 @'2 @z2
@
;( e + e ) @z @ur + ur + 1 @u' + @uz
@r r r @' @z

= fz

fur alle (r ' z ) 2

e

(1
...
Es bezeichnen f = (fr f' fz )T den Vektor der Volumenkrafte, ~1 = (g1r g1' g1z )T
g
den Vektor der vorgegebenen Randverschiebungen, ~2 = (g2r g2' g2z )T den Vektor
g
der vorgegebenen Ober achenkrafte und ~ = (nr n' nz )T den Vektor der au eren
n
Einheitsnormalen
...
EINFUHRUNG

Zwischen den Verschiebungen (ur u' uz )T und dem Verschiebungsfeld (u1 u2 u3)T in
kartesischen Koordinaten besteht der Zusammenhang
u1 = ur cos ' ; u' sin '
u2 = ur sin ' + u' cos '
u3 = uz
sowie zwischen den Zylinderkoordinaten (r ' z) und den kartesischen Koordinaten die
Beziehung
x1 = r cos ' x2 = r sin ' x3 = z :

p

Innenradius rI : 0
...
2 m

Elastizitatsmodul E
Querkontraktionszahl
Druck p

: 6:87 107 Nm;2
:
0
...
8: Geometrie des Rohres und Darstellung der verformten Kontur
Auf Grund der Geometrie des Rohres und der Wirkung des Druckes gilt
u' = uz = @ur = @ur = 0
@' @z
so da sich das Di erentialgleichungssystem (1
...
h
...
34)
fur r = rI
rr = ;p
gilt
...
3
...
34) ist
ur = C1r + C2 :
r
Die unbekannten Konstanten C1 und C2 werden unter Nutzung der Randbedingungen
bestimmt (Beziehungen (1
...
32)), d
...

( e + 2 e ) @ur + e ur = ;p fur r = rI
@r
r
( e + 2 e ) @ur + e ur = 0 fur r = rA :
@r
r
Wir erhalten schlie lich die Losung
2
2 r2
1
ur = p 2 + 2 r2 rI r2 r + 21 r2rI;Ar2 1 :
e
e A; I
e A
I r
Das Verschiebungsfeld fur das Rohr unter Innendruck konnten wir analytisch berechnen
...
Fur eine naherungsweise Bestimmung des Verschiebungsfeldes mussen wir eine Diskretisierung des jeweiligen Problems vornehmen
...
3
...
Wir berechnen
zuerst ein stationares Temperaturfeld und bestimmen dann die durch den Temperaturgradienten verursachten Verschiebungen
...
35)
@T = gT (x)
fur alle x 2 ;T
2
2
@N
@T + (x)T (x) = (x)T (x) fur alle x 2 ;T
A
3
@N
gilt, sowie das Verschiebungsfeld ~ (x), welches das Di erentialgleichungssystem
u
~
; e ~ (x);( e + e ) grad div ~ (x) = f e (x);(3 e +2 e ) l (x) grad T (x) fur alle x 2
u
u
(1
...
37)
e
fur alle x 2 ;e
21n1 + 22n2 + 23n3 = g22 (x)
2
e
fur alle x 2 ;e
31n1 + 32n2 + 33n3 = g23 (x)
2

30

KAPITEL 1
...
Fur den Rand des Gebietes gilt @ = ;T ;T ;T = ;e ;e mit ;T \
1
2
3
1
2
i
;T = ;e \ ;e = , i j = 1 2 3
...
Die Spannungskomponenten bei thermomechanischen Feldern sind
durch
i=1 2 3
ii = e ("11 + "22 + "33) + 2 e "ii ; (3 e + 2 e ) l T
i j = 1 2 3 i 6= j
ij = ji = 2 e "ij
de niert, mit den in den Beziehungen (1
...

Wir berechnen im weiteren das Verschiebungsfeld im oberen Teil des Kolbens eines
Verbrennungsmotors
...
Wir nehmen an, da der obere Teil
des Kolbens rotationssymmetrisch ist
...
9 und 1
...
Weiterhin setzen wir voraus, da alle Eingangsdaten
wie z
...
die Warmeleitzahl, die Warmeubergangszahl, die Umgebungstemperatur, die
Lameschen Konstanten, die vorgegebenen Ober achenkrafte und der lineare Warmeausdehnungskoe zient vom Rotationswinkel unabhangig sind
...
9: Niveaulinien des Temperaturfeldes

Zur Bestimmung des Temperatur{ und Verschiebungsfeldes nutzen wir die Formulierung des Warmeleitproblems und des Elastizitatsproblems in Zylinderkoordinaten
...
Wir erhalten die

31

1
...
BEISPIELE

folgenden Randwertprobleme (siehe auch die Beziehungen (1
...
33)):
Gesucht ist das Temperaturfeld T (r z), fur das
; (r z )

1 @ r @T + @ 2T = 0
fur alle (r z) 2
r @r @r
@z2
T
T (r z) = g1 (r z) fur alle (r z) 2 ;T
1
@T = 0
fur alle (r z) 2 ;T
2
@N

gilt, mit ;T = f(r z) 2 @ : r = 0g f(r z) 2 @ : z = 0g und ;T = @ n ;T
...

1
T wird durch die Losung der folgenDas Verschiebungsfeld ~ (r z) = (ur(r z) uz (r z))
u
den Randwertaufgabe bestimmt
...
Die Komponenten des Spannungstensors sind durch die Beziehungen
rr = e ("rr + "zz ) + 2 e "rr ; (3 e + 2 e ) l T
zz = e ("rr + "zz ) + 2 e "zz ; (3 e + 2 e ) l T
rz = zr = 2 e "rz
mit den Verzerrungen (1
...


32

KAPITEL 1
...
Weiterhin gilt
e
g2r = 0 auf ;e r
21

und

e
g2z = 0 auf ;e z :
21

e
e
Auf den Teilrandern ;e r und ;e z sind die Funktionen g2r sowie g2z durch den Gasdruck
22
22
gegeben
...
10: Darstellung der Deformation der Kontur (stark vergro ert)

Kapitel 2

Modellierung am Beispiel von Temperaturfeldern
2
...
1
...
B
...


Aufheizung (z
...
durch Umwandlung
von elektrischer Energie in Warme)
Kreis j@Qj = 2 r CC
Warmeabgabe uber
C
den Mantel
@Q
C
xC
6 6
6
C
C

a

Q = const:
x;
2
Kreis jQj = r

ga
Temperatur im linken
Randpunkt
u(a) = ga

;

x
2

x

uA = uA(x)
Umgebungstemperatur

@

x+

x
2

b

? ?
?

gb
Temperatur im rechten
Randpunkt
u(b) = gb

Abbildung 2
...
MODELLIERUNG AM BEISPIEL VON TEMPERATURFELDERN

Warmemengebilanz an einem "kleinen\ Stabstuck der Lange x (spater x ! 0)

2
...
ein ie ende
Warmemenge

-

1
...
aus ie ende
Warmemenge
???????????

x;

x
2

-

x+

x

x
2

Abbildung 2
...
Warmemenge, die durch Aufheizung entsteht:
xZ 2x
+
WH = jQj
f( ) d
x;

x

2

wobei f ( ) die Funktion der Intensitat der Warmequelle bezeichnet
...
Transportierte Warmemenge beim Warmeaustausch mit der Umgebung uber den
Mantel:
x+ 2x
xZ 2x
+
Z
j@Qj q ( ) (u( ) ; u ( )) d :
WA = j@Qj
q( ) (u( ) ; uA( )) d = jQj
A
jQj{z }
x
x |
x; 2
x; 2
q( )
Der Warmeaustauschkoe zient q( ) ist materialabhangig
...
Warmemenge, die am linken Rand des Stabstuckes, d
...
bei (x ; 2x ), in " x\
hinein ie t bzw
...
h
...
Fuhren wir als Proportionalitatsfaktor die Warmeleitzahl (x) ein, dann gilt
(x) = ; (x) u0(x) :
Somit ergibt sich fur die ein ie ende Warmemenge

W (x ; 2x ) = ; (x ; 2x ) u0(x ; 2x ) jQj
und fur die aus ie ende Warmemenge
W (x + 2x ) = ; (x + 2x ) u0(x + 2x ) jQj :

35

2
...
STATIONARE WARMELEITPROBLEME

Aus der Warmemengebilanz
Warmemenge,
Warmemenge,
Warmemenge,
die am linken
die am rechten
Warmemenge,
die uber den + die durch Auf- = 0
Rand, d
...
bei ; Rand, d
...
bei ;
Mantel abge(x ; 2x ), in
(x + 2x ), aus
heizung entgeben wird
\
\
steht
" iext hinein" iext heraus-

W (x ; 2x )

;

W (x + 2x )

WA

;

(2
...

Die Warmeleitgleichung in integraler Form

Fur alle x 2 (a b) und fur alle x > 0 , fur die x ;
Gleichung
; (x ;

x 0
x
2 )u (x ; 2 )jQj

; jQj

xZ
+

x;

+ (x +

x
2

x+

x 0
x
2 )u (x + 2 )jQj

x

2

q( )(u( ) ; uA ( )) d = ; jQj

xZ
+

x;

x

2

(2
...

Die di erentielle Form der Warmeleitgleichung bei homogenem Material und stetiger
Warmequelle

Wir gehen von folgenden Voraussetzungen aus:
1
...
h
...
B
...
0 (q = 0 entspricht einer Isolation)
2
...
h
...
> 0, z
...
= const
...
f 2 C a b], d
...
f (x) ist eine stetige Funktion im Intervall a b]
4
...
h
...
3)

Aus dem Mittelwertsatz der Integralrechnung (siehe 19]) folgt, da fur eine beliebige
Funktion v(x) 2 C a b] und fur ein beliebiges x 2 (a b) die Beziehung
lim0 1x
x!
gilt
...
4)

36

KAPITEL 2
...
2) mit ;1 und dividieren sie durch x sowie jQj,
dann erhalten wir
x+ 2x
Z
x )u0 (x + x ) ; (x ; x )u0(x ; x )
(x + 2
2
2
2 + 1
;
x
x x q( )u( ) d
x;

2

x+ x
x+ x
1 Z 2 f ( ) d + 1 Z 2 q( )u ( ) d :
= x
A
x x
x
x;

2

x;

2

Mit x ! 0 ergibt sich auf Grund der Voraussetzungen (2
...
4)
die Di erentialgleichung
;( (x)u0(x))0 + q (x)u(x) = f (x) + q(x)uA (x):

Bemerkung 2
...
bis 4
...
h
...
h
...

Die klassische di erentielle Form der Warmeleitgleichung lautet
Gesucht ist u(x) 2 C 2(a b) \ C a b], so da
; ( (x)u0 (x))0 + q(x)u(x) = f (x) + q (x)uA(x) fur alle x 2 (a b)

(2
...

Die di erentielle Form der Warmeleitgleichung bei stuckweise homogenem Material
und stuckweise stetiger Warmequelle

Wir betrachten o
...
d
...
einen Stab, bestehend aus zwei Materialien
...
3: Stab, bestehend aus zwei verschiedenen Materialien
Dabei setzen wir voraus, da 1, q1, uA1 und f1 die Voraussetzungen (2
...


37

2
...
STATIONARE WARMELEITPROBLEME

Eine zur Herleitung der Warmeleitgleichung (2
...

Gesucht ist das Temperaturfeld u(x), fur das die folgenden Beziehungen gelten:
Randbedingung
u(a) = ga
fur x = a
Di erentialgleichung
;( 1u0)0 + q1u = f1 + q1uA1
fur alle x 2 (a )
Interfacebedingungen
u( ; 0) = u( + 0)
fur x =
; 1( ; 0)u0( ; 0) = ; 2( + 0)u0( + 0) fur x =
Di erentialgleichung
;( 2u0)0 + q2u = f2 + q2uA2
fur alle x 2 ( b)
Randbedingung
u(b) = gb
fur x = b
Die Interfacebedingung ; 1( ; 0)u0( ; 0) = ; 2( + 0)u0( + 0) folgt unmittelbar
aus der Gleichung (2
...


Bemerkung 2
...
Art:
Freier Warmeubergang
Der Warme u ist proportional zur Di erenz
zwischen dem Wert der Temperatur am Rand
und der Umgebungstemperatur
...
rechten Rand
{ Gemischte Randbedingungen: z
...
u(a) = ga, u0(b) = 0 oder
u0(a) = 0, ; (b)u0(b) = b (u(b) ; ub)

{ Randbedingungen 1
...
Art:
(Neumannsche RB):

Bemerkung 2
...
Aus der Warme-

)u0( p ;

2

p

)jQj + ( p +

2

P

p

q( )(u( ) ; uA( )) d = ; jQj

)u0( p +

Z

p+ 2

p; 2

p

2

p

)jQj

p

f ( ) d ; jQj

Z

p+ 2

p; 2

p

p

(2
...
MODELLIERUNG AM BEISPIEL VON TEMPERATURFELDERN

ergibt sich fur

p

! 0 die Beziehung

^
; ( p ; 0) u0 ( p ; 0) = ; ( p + 0) u0 ( p + 0) ; f:
Hier bezeichnet ( ; p) die Diracsche Deltafunktion (siehe auch 55, 88])
...

Gesucht ist das Temperaturfeld u(x), fur das die Beziehungen
Randbedingung
u(a) = ga
fur x = a
Di erentialgleichung
;( u0)0 + qu = f + quA
fur alle x 2 (a p)
Interfacebedingungen
u( ; 0) = u( + 0)
fur x = p
; ( ; 0)u0( ; 0) = ; ( + 0)u0( + 0) ; f^ fur x = p
Di erentialgleichung
;( u0)0 + qu = f + quA
fur alle x 2 ( p b)
Randbedingung
u(b) = gb
fur x = b
gelten
...
1
...


x3

6
;
;

x2

;
;
;
;
;
;;
;
;;
;;
;
;
x1
;;;
;;
H
;
;
;;;
;;
;
x2 ;; ; ; ;6
;;;
;;
;;
; h
?
;
;;;
;;
6
;;;
;;
h ;;;
;;;
2 ;
;;
;
;h ;

r

2

@
(Rand der Mittel ache )
-

x1

P (x1 x2 0)
Abbildung 2
...
1
...
h
...


Warmeaustausch
mit der Umgebung
innere Warmequellen
x
1 (x1 ; 2 1

2)

2( 1

6

;
;
PP
;
PP
;
P
P
P
;
q
P
;
; ;
;

;
;
6
;
;
;
;
;
;
?
;
;
;; ;
;;
;;
2
;;
;;
;;
; ;

h

rP

2( 1

;

x2 ;

x
1 (x1 + 2 1

2)

x

;

;
;
;

x1

;
;

x2
2 )

;
;

;

;

x2 +

-

Warmeaustausch
mit der Umgebung

?

x2
2 )

Abbildung 2
...
7)

2

x2

2

2) d 2

2

x1

2

x1

x1 +
Z

x2 Z
+

q( 1 2)(u ; uA) d 2 d 1 + h

x1 +
Z

x1 ;

x1

2

x1

2

x2 +
Z

x2 ;

x2

2

f ( 1 2) d 2 d 1 = 0:

x2

2

Zusatzlich sind noch die Randbedingungen auf @ vorgegeben
...
MODELLIERUNG AM BEISPIEL VON TEMPERATURFELDERN

Nach dem Fourierschen Erfahrungsgesetz der Warmeleitung gilt fur ein orthotropes
Material
@u
@u
und
1=; 1
2=; 2
@x1
@x2
mit den Warmeleitkoe zienten 1 und 2
...
1(x) 2(x) 2 C 1( ), d
...
1 und 2 sind einmal stetig di erenzierbare Funktionen
uber dem Gebiet
1 2
0 = const
...
f (x) 2 C ( ), d
...
f ist eine stetige Funktion uber
3
...
uA(x) 2 C ( )
...
7) mit ;1=(h x1 x2) und lassen x1 sowie x2
gegen Null streben, dann erhalten wir die Warmeleitgleichung in di erentieller Form:
Gesucht ist u(x) 2 C 2( ) \ C ( ) , so da
@u ; @
@u + q(x)u(x) = f (x) + q(x)u (x)
@
(2
...


Bemerkung 2
...
Art:

@u n ; @u n = 0 auf ; = @
@x1 1 2 @x2 2

(a) Warmeisolation:

;

(b) vorgegebener Warme u :

; 1 @u n1 ;

1

@x1

2

@u n = g auf ;
@x2 2 2

Dabei bezeichnen ~ (x) = (n1(x) n2(x))T (j~ j = 1) den Vektor der au eren Einn
n
heitsnormalen im Punkt x = (x1 x2) 2 ; und 1(x) sowie 2(x) die Warmeleitkoe zienten
...
Art: Freier Warmeaustausch mit der Umgebung
@u
@u
; 1 n1 ; 2 n2 = (u ; uA ) auf ;
@x1
@x2
Es bezeichnen (x) die Warmeubergangszahl und uA(x) die Umgebungstemperatur
...
1
...
6: Vektor der au eren Einheitsnormalen am Rand ;
{ Gemischte Randbedingungen
Auf ;1 sind Randbedingungen 1
...
Art und auf ;3
Randbedingungen 3
...


s

;1

s

;3
;2

s

; = ;1 ;2 ;3
;i \ ;j = fur i 6= j i j = 1 2 3

Abbildung 2
...
5 Ein Spezialfall

Fur thermisch isotropes, homogenes Material, d
...
1 = 2 = = const: > 0, sowie
q 0 erhalt man aus der Gleichung (2
...


42

KAPITEL 2
...
1
...


x3

n1 !
~ (x) = n2
n
n3
I
@

6

innere Warmequellen

@r

x2 ;
;

x1

r

;
;

6

;
;

;
;
;

;

;=@

;
;
6 x3
?

Umgebungstemperatur uA
@

x2

Randbedingungen:
z
...
u(x) = g1 (x) x 2 ;

P (x1 x2 x3)

-

x1

Abbildung 2
...
9: Warmemengebilanz an einem Quader " x1

x2

x3\

43

2
...
STATIONARE WARMELEITPROBLEME

Fur alle x = (x1 x2 x3) 2 und fur alle x1, x2, x3 > 0 mit f(y1 y2 y3) :
x
x
x ; 2 y x + 2 = 1 2 3g
gilt
x2 Z
+
x2 ;

+

+

x2

2

x1 +
Z

x1 ;

x3 ;

x1

2

x1

2

x1 Z
+

x1 ;

x3 +
Z

x2

2

x1

2

x1

2

x3

2 3) d 3 d 2

;

x2;

2( 1

x2 ;

x3

x2
2

3) d 3d 1

;

3( 1 2

x3 ;

x2

2

x3
2 ) d 2d 1

;

x1 ;

2

+

x1 +
Z

x1 ;

x1

2

x1

2

x2 +
Z

x2 ;

x2

2

x2

2

x3 ;

x1

2

x3 ;

x1

x2 ;

x3 +
Z

2

x3 ;

x3

x
1 (x1 + 2 1

2 3) d 3 d 2

2

x3

2

2( 1

x2 +

x3

x2
2

3) d 3d 1

2

x2 Z
+

x1

2

2

x3

2

x3 Z
+

x1

2

x1 +
Z

x2

2

x2

x1 ;

2

x3 +
Z

x2

2

x1 +
Z

x3

2

x2 Z
+

x2 ;

x
1 (x1 ; 2 1

2

x3 Z
+

x3 ;

x2 Z
+

x3

2

(2
...

Auf Grund des Fourierschen Erfahrungsgesetzes der Warmeleitung gelten fur ein orthotropes Material die Beziehungen
@u i = 1 2 3 :
i=; i
@xi
Unter den Voraussetzungen:
1
...
h
...
0 > 0 und
2
...
h
...
9) mit ;1=( x1 x2 x3) und
dem Grenzubergang x1 ! 0, x2 ! 0 sowie x3 ! 0 die Warmeleitgleichung in
di erentieller Form:
Gesucht ist u(x) 2 C 2( ) \ C ( ) , so da
3
X @
;
i=1 @xi

gilt
...
10)

44

KAPITEL 2
...
6 Moglichkeiten zur Formulierung anderer Randbedingungen
{ Vorgabe des Warme usses (Randbedingungen 2
...
Art)
@u
; = (x) (u(x) ; uA (x)) fur alle x 2 ;
@N
Es bezeichnen (x) die Warmeubergangszahl und uA(x) die Umgebungstemperatur
...


Bemerkung 2
...
Herleitung der Warmeleitgleichung mittels des Gau schen Integralsatzes:
(falls @G ;2 ;3
Wir betrachten ein beliebiges Teilgebiet G
mit G
gilt, sind Sonderbetrachtungen notwendig)
...
B
...


:

45

2
...
STATIONARE WARMELEITPROBLEME

2
...
h
...
10) in die Poisson-Gleichung
; u = f (x) fur alle x 2
u(x) = g(x) fur alle x 2 ; = @
@2 @2 @2
mit dem LAPLACE-Operator = @x2 + @x2 + @x2 und f (x) = f (x)= uber
...
Inhomogenes, anisotropes Material (Materialsprunge):
{ Die Warmeleitgleichung in integraler Form gilt nach wie vor!
{ Die Di erentialgleichung gilt nur in Teilgebieten mit "glatten\ Daten
...
h
...

Spezialfalle bei 3D{Warmeleitproblemen

1
...
10: Beispiel fur einen "zylinderartigen\ Korper
Somit mu die Losung von x3 unabhangig sein
...
10) ergibt sich
deshalb die folgende Aufgabe:
Gesucht ist u(x) 2 C 2( a) \ C ( a), so da
@
@u ; @
@u = f (x)
;
1 (x)
2 (x)
@x1
@x1 @x2
@x2
u(x) = g1(x)

fur alle x = (x1 x2) 2
fur alle x = (x1 x2) 2 @

a
a

gilt
...
Rotationssymmetrische Probleme, d
...
rotationssymmetrisches Gebiet und Unabhangigkeit der Eingangsdaten vom Rotationswinkel
Wir gehen von den kartesischen Koordinaten (x1 x2 x3) zu Zylinderkoordinaten
(r ' z) uber ( ! a 0 2 )) und erhalten aus (2
...
MODELLIERUNG AM BEISPIEL VON TEMPERATURFELDERN

Gesucht ist u(r z), so da
1 @ r (r z) @u ; @ (r z) @u = f (r z) fur alle (r z) 2
;
a
r @r
@r @z
@z
u(r z) = g1(r z) fur alle (r z) 2 ;1
;1 = @

@u
@r

; (r z ) nr (r z ) ; (r z )

gilt
...
11: Meridianebene a eines rotationssymmetrischen Korpers

2
...
2
...
h
...

t6
t + 2t
6
t
t
?
t ; 2t
Q
@
@

a

x;

x
2

r

x
x

-

x

x+

x
2

b

-

Abbildung 2
...
2
...
h
...


Warmemenge, die
am linken Rand,
d
...
bei x ; 2x , in
" x\ wahrend der
Zeitspanne t
hinein ie t
; jQj

tZ
+

t

2

t; 2t

@u
@x

Warmemenge, die
am rechten Rand,
x
; d
...
bei x + 2 aus ;
" x\ wahrend der
Zeitspanne t
heraus ie t

d + jQj

t+
Z

t; 2t

x ; 2x

t

2

@u
@x

d ; jQj

+ jQj

t;

t

2

t

2

xZ
+

x;

fd d

t

2

xZ
+

x

2

q(u ; uA ) d d

x

2

Warmemengendi erenz
in " x\ zwischen End{
und Anfangszeit

x

2

t+
Z

t; 2t x;

x + 2x

Warmemenge, die
+ durch Aufheizung =
wahrend der Zeitspanne t entsteht
t+
Z

Warmemenge, die
uber den Mantel
wahrend der Zeitspanne t abgegeben wird

= jQj

xZ
+

x;

x

2

x

2

x

c% (u( t+ 2t ) ; u( t; 2t )) d

2

Hier bezeichnen % die Dichte, c die spezi sche Warmekapazitat, den Warmeleitkoefzienten und q = j@Qjj q den spezi schen Warmeaustauschkoe zienten
...

Fur alle x 2 (a b) und fur alle x > 0 mit x ; 2x x + 2x ] (a b) sowie fur alle
t 2 (0 T ) und fur alle t > 0 mit t ; 2t t + 2t ] (0 T ) gilt
xZ
+
x;

x

2

x

c% (u( t + 2t ) ; u( t ; 2t )) d ;

2

+

tZ
+

t;

t

2

t

2

xZ
+

x;

t+
Z
t;

t

2

t

2

@u
@u
;
@x jx+ 2x @x jx; 2x d

x

2

x

2

q( ) (u( ) ; uA( )) d d =

t+
Z
t;

t

2

t

2

x+
Z

x;

x

2

(2
...
B
...
Art)
u(a t) = ga(t) und u(b t) = gb(t) fur alle t 2 0 T ]
sowie die Anfangsbedingung
u(x 0) = u0(x) fur alle x 2 (a b)
gegeben
...
MODELLIERUNG AM BEISPIEL VON TEMPERATURFELDERN

Ubergang zur di erentiellen Form
Multiplizieren wir die Gleichung (2
...
B
...
3)) die di erentielle
Form der instationaren Warmeleitgleichung
...
12)
@t
@x
fur alle (x t) 2 QT
gilt und die Randbedingungen (z
...
1
...


t6
T
QT

0

a

-

b x
Abbildung 2
...
8 Der Raum C 2 1(QT ) umfa t alle Funktionen, die zweimal stetig dif-

ferenzierbar bezuglich des Ortes x 2 (a b) und einmal stetig di erenzierbar bezuglich
der Zeit t 2 (0 T ) sind
...
9 Analog zu der stationaren Warmeleitgleichung konnen auch bei der
instationaren Warmeleitgleichung Randbegingungen 2
...
Art bzw
...

Bemerkung 2
...
h
...
12), ist die Kompatibilitat zwischen der Anfangsbedingung und den Randbedingungen notwendig, d
...

lim g (t) = u0(a)
lim g (t) = u0(b)
t!+0 a
t!+0 b

49

2
...
INSTATIONARE WARMELEITPROBLEME

Bemerkung 2
...
2
...
auch die Abschnitte 2
...
1 und 2
...
2 bzw
...
1
...


Gesucht ist u(x t) 2 C 2 1(QT ) \ C (QT ) , so da
m
X @
@u
c(x t) %(x t) @u ; @x i (x t) @x = f (x t)
@t i=1 i
i
fur alle (x t) 2 QT =

(2
...
B
...
Art)

u(x t) = g1(x t) fur alle (x t) 2 @

0 T]

sowie die Anfangsbedingung

u(x 0) = u0(x)

fur alle x 2

erfullt werden (x = (x1 x2) bzw
...


Bemerkung 2
...


und einmal stetig di erenzierbar bezuglich

Bemerkung 2
...
Art oder 3
...
gemischte Randbedingungen gestellt werden
...
MODELLIERUNG AM BEISPIEL VON TEMPERATURFELDERN

Instationare (stationare) Warmeleit{Warmetransportprobleme

Gesucht ist u(x t) 2 C 2 1(QT ) \ C (QT ) , so da
m
m
X @
X
@u
@u
c(x t)%(x t) @u ; @x i(x t) @x +
ai(x t) @x + q(x t)u(x t)
@t i=1 i
i
i
i
| {z }
|
|=1 {z
{z
}
}
Warmetransport nur im 1D{ bzw
...

im 2D{Fall

fur alle (x t) 2 QT =

(0 T )

gilt und die Randbedingung

u(x t) = g1(x t)

fur alle (x t) 2 @

0 T]

sowie die Anfangsbedingung

u(x 0) = u0(x)

fur alle x 2

erfullt werden
...
C den VektorGebiet ie t (siehe auch welcher das Material
a @
A durch das
55]),
am(x t)
q(x t)
den spezi schen Warmeaustauschkoe zienten,
f (x t)
die Warmeintensitatsfunktion,
g1(x t)
die vorgegebene Randtemperatur und
u0(x)
die Anfangstemperatur
...
14 Bei der stationaren Warmeleit{Warmetransportgleichung wird

die Losung u als eine Funktion u(x) gesucht
...
Weiterhin entfallen der Term c% @u in der Di erentialgleichung sowie die
@t
Anfangsbedingung
...
1 Die Funktionenraume 2( ) und
L

a b

H

1 (a

)

b

Eine Funktion u(x) ist Element des Funktionenraumes L2(a b), wenn das Integral
Zb
(u(x))2 dx
a

existiert und endlich ist
...

u 6

1

;1 0

r

1

-

x

8
> 1 fur 0 < x 2
>
<
u(x) = sgn(x) = > 0 fur x = 0
> ;1 fur ;2 x < 0
:

;1

Abbildung 3
...

Fur beliebige Funktionen u(x) und v(x) aus L2(a b) gilt die Cauchy{Schwarzsche Ungleichung
v
v
u
uZb
Zb
u
u
u (u(x))2 dx uZb (v(x))2 dx:
t
u(x)v(x) dx t
a

Mit

a

a

v
uZb
u
t
jjujj0 2 (a b) = u (u(x))2 dx
a

ist eine Norm in L2(a b) de niert und mit

52

KAPITEL 3
...

Im weiteren wollen wir den Sobolev{Raum H 1(a b) de nieren
...

Wiederholung: Besitzt eine Funktion u(x) eine stetige Ableitung u0(x), so gilt nach
der Formel der partiellen Integration fur jede stetig di erenzierbare Funktion '(x) mit
'(a) = '(b) = 0
Zb
Zb
0
u(x)' (x) dx = ; u0(x)'(x) dx:
a

a

a

a

Wir de nieren mit Hilfe dieser Formel Ableitungen fur Funktionen, die im ublichen
klassischen Sinn nicht di erenzierbar sind
...

u 6

1

8
< x + 1 fur ;1 x 0
u(x) = :
1 ; x fur 0 < x 1

;@
;
@
;
@
;
@
;
@
;
@ -

;1

0

(3
...
2: Ein Beispiel einer Funktion aus dem Raum H 1(a b)
Die durch die Beziehung (3
...
Fur di erenzierbare Funktionen '(x)
mit '(;1) = '(1) = 0 erhalten wir mittels der Formel der partiellen Integration

Z1
;1

u(x) ' (x) dx =
0

Z0
;1

(x + 1) ' (x) dx +
0

Z0

Z1
0

(1 ; x) '0(x) dx

= ; '(x) dx + (x + 1)'(x)
;1

0

;1

Z1

; (;1) '(x) dx + (1 ; x)'(x)
0

3
2 Z0
Z1
= ; 4 '(x) dx + (;1) '(x) dx 5 + '(0) {z '(0) :
| ; }
0
;1
=0

1
0

53

3
...
VARIATIONSFORMULIERUNG VON RANDWERTAUFGABEN

Damit ist im verallgemeinerten Sinn
8
< 1 fur ;1 x < 0
u0(x) = :
:
;1 fur 0 < x 1
Im Punkt x = 0 kann u0(0) beliebig festgesetzt werden
...
Solche Funktionen treten gerade bei der Methode der niten
Elemente auf (siehe auch Abschnitt 3
...

Eine Funktion u(x) aus dem Raum L2(a b) ist Element des Raumes H 1(a b), wenn
u0(x) existiert und zum Raum L2(a b) gehort
...

Fur Funktionen u(x), die wenigstens in einem der Randpunkte x = a bzw
...

u(x H 1 ) nimmt im
"Die FunktionWert)g 2bzw
...
gx
...
b) =
a

b

c

a

c

b

Eine ausfuhrliche Darlegung der Theorie der Sobolev{Raume kann der interessierte
Leser z
...
in 1] und 86] nden
...
2 Variationsformulierung von Randwertaufgaben
Die Grundidee der Finite{Elemente{Methode beschreiben wir im weiteren anhand eines
stationares 1D{Warmeleitproblems (siehe Abschnitt 2
...
1)
...
B
...
Die klassische Formulierung des Warmeleitproblems
lautet somit:

54

KAPITEL 3
...
2)
;u0 (b) = b (u(b) ; gb )
gilt, wobei ga und gb vorgegebene Werte sind und f (x) eine gegebene Funktion ist
...

Wir wahlen eine beliebige Funktion v(x) 2 H 1(a b) , wobei v(x) in den Randpunkten mit Randbedingungen 1
...
h
...
2) mit einer solchen Testfunktion v(x)
...
3)
a

a

Gema der Formel der partiellen Integration
Zb
Zb
b
00
; u (x)v (x) dx = u0(x)v 0(x) dx ; u0(x)v (x) a
a

a

(3
...
5)
a

a

Aus den Beziehungen (3
...
5) folgt die Variationsformulierung (= schwache Formulierung = verallgemeinerte Formulierung) des Warmeleitproblems, d
...
:
Gesucht ist u(x) 2 Vg , so da
Zb
Zb
0
0
u (x)v (x) dx + bu(b)v(b) = f (x)v(x) dx + bgbv(b) fur alle v(x) 2 V0 (3
...


Bemerkung 3
...
6) enthalt im Unterschied

zur klassischen Formulierung (3
...
Ableitungen mehr
...
6) ist der Ausgangspunkt fur die FEM (Finite{Elemente{Methode)
und die Aufgabe (3
...
B
...


Bemerkung 3
...
Art) gehen

in die De nition von Vg und V0 ein, naturliche Randbedingungen (Randbedingungen
2
...
3
...


55

3
...
DIE FEM ZUR NAHERUNGSWEISEN LOSUNG DES VARIATIONSPROBLEMS

3
...


u

6

u = '(x)
-

|

{z
Trager

}

x

supp '(x) =

] (lokaler Trager)

Abbildung 3
...

Die Idee der FEM: Wie wird die Naherungslosung konstruiert ?

1
...
a
...


a
k
x0

xi;1

xi xi+1

b
k
xn

Abbildung 3
...
h
...
Jedem Knoten Pi ordnen wir eine im Intervall a b] stetige, nite Funktion 'i(x)
(Ansatzfunktion) zu, die nur im Intervall (xi;1 xi+1) von Null verschieden ist
...
5
dargestellten stuckweise linearen Funktionen
...
GRUNDPRINZIPIEN DER FINITE{ELEMENTE{METHODE: EIN 1D{BEISPIEL

u
1

'0(x)

6

a
k
x0
u
1

x2

xi;1 xi xi+1

-

x

xn;1 b
k
xn

'i(x)

6

a
k
x0
u
1

x1

x1

x2

xi;1 xi xi+1

-

xn;1 b
k
xn

x

'n(x)

6

a
k
x0

x1

x2

xi;1 xi xi+1

xn;1 b
k
xn

-

x

Abbildung 3
...
'n(x)
8 x ; x1
80
>;
>
<
<
a < x x1
h
'0(x) = >
'n(x) = > x ; xn;1
:0
:
x1 < x b
h

(3
...
8)

57

3
...
DIE FEM ZUR NAHERUNGSWEISEN LOSUNG DES VARIATIONSPROBLEMS

Die Ansatzfunktionen 'i(x), i = 0 1 : : : n, erfullen o enbar die Bedingung
(
'i(xj ) = ij = 1 fur i = j
(3
...

3
...
6) als Linearkombination

uh(x) =

n
X
i=0

ui'i(x)

(3
...

Der Parameter u0 ist in unserer Aufgabe auf Grund der vorgegebenen Randbedingung u(a) = ga bekannt
...

Setzen wir den Losungsansatz (3
...
6) ein, dann
erhalten wir zunachst die Ersatzaufgabe:
Gesucht ist uh(x) 2 Vgh , so da
Zb
Zb
0
0
uh(x)vh(x) dx + buh(b)vh(b) = f (x)vh(x) dx + bgb vh(b)
a

a

(3
...

Die Losung der Aufgabe (3
...
Diese sind die Losung des Finite{Elemente{
Gleichungssystems, welches aus der Ersatzaufgabe (3
...

Bestimmung der Knotenparameter u1 u2 : : : un

Wir wahlen nacheinander in der Aufgabe (3
...
Damit ergibt sich das zur Aufgabe (3
...
GRUNDPRINZIPIEN DER FINITE{ELEMENTE{METHODE: EIN 1D{BEISPIEL

Zb
n
X Zb 0 0
vh = '1 : ui 'i(x)'1(x) dx + buh(b)'1(b) = f (x)'1(x) dx + bgb'1(b)
i=0

a

a

Zb
n
X Zb 0 0
vh = '2 : ui 'i(x)'2(x) dx + buh(b)'2(b) = f (x)'2(x) dx + bgb'2(b)
i=0

a


...


...
Beachten wir au erdem
'j (b) = 0 fur j = 1 2 : : : n ; 1, 'n (b) = 1 und uh(b) = un , dann folgt
n
X Zb 0 0
ui 'i (x)'1(x) dx
i=1

=

a

n
X Zb 0 0
ui 'i (x)'2(x) dx
i=1

=

a


...


...
12)

Zb
a

Zb
a

Zb

f (x)'n;1 (x) dx ; ga '00(x)'0n;1(x) dx
a

f (x)'n (x) dx

Zb

; ga '00 (x)'0n (x)dx
a

+

b gb :

Die Beziehungen (3
...
13)

zur Bestimmung des unbekannten Koe zientenvektors uh = (u1 u2 : : : un)T 2 Rn
bezeichnet
...
3
...
11) und (3
...
13) sind auf Grund der folgenden Uberlegungen
aquivalent
...
11) suchen wir eine Funktion uh 2 Vgh , so da die
Integralidentitat fur jede Funktion vh 2 V0h gilt, also auch fur jede Funktion 'i(x),
i = 1 2 : : : n
...
11) die
Gleichungen (3
...
Bestimmen wir einen Losungsvektor uh = (u1 u2 : : : un)T des
n
Gleichungssytems (3
...
13) und somit u = P u ' (x) , d
...
eine Funktion u , die
h

i=0

i i

h

den Gleichungen (3
...
11) fur alle vh 2 V0h , da sich jede Funktion vh 2 V0h als Linearkombination
n
v = P v ' (x) darstellen la t
...

1
...
h
...

2
...
h
...
7) und (3
...
6: Veranschaulichung des gemeinsamen Tragers der Funktionen 'i;1 (x),
'i(x) und 'i+1(x)

60

KAPITEL 3
...
K ist eine tridiagonale Matrix
0
B
B
B0
B
B
...

B
B
B0
B
B0
@
0

0

0
0

...


0
0
0

...


...
7) und (3
...
4
...


...


...
CB
B
CB
hB 0
B
0 ;1 2 ;1
0 CB
B
CB
B 0
0 0 ;1 2 ;1 C B
@
A@
0
0 0 0 ;1 1 + b h

1 0 f~1 + 1 ga 1
C B f~2h C
C
C B
C B f~3 C
B
C
C B
C
C B
C
C=B

...

C B

...

mit f
xi;1

u1
u2
u3

...


...
4 Der elementweise Aufbau der Stei gkeitsmatrix und des
Lastvektors
Wir beschreiben im folgenden wie die Stei gkeitsmatrix und der Lastvektor elementweise berechnet werden konnen
...
Die Beachtung der Randbedingungen erlautern wir spater
...

Berechnung der Elementbeziehungen
Auf Grund der im Abschnitt 3
...
10) sowie unter Beachtung der De nition der Funktionen 'i(x) (siehe (3
...
8))

Zb
a

uh(x) vh(x) dx =
0

0

=

n
X Zxi
i=1 xi;1
n
X
i=1

(ui;1'0i;1(x) + ui'0i(x))(vi;1'0i;1(x) + vi'0i(x)) dx

(vi;1 vi)K (i)

ui;1
ui

!

(3
...
GRUNDPRINZIPIEN DER FINITE{ELEMENTE{METHODE: EIN 1D{BEISPIEL

Die Matrizen K (i), i = 1 2 : : : n,
0 (i) (i) 1
K K
K (i) = @ 11) 12) A
(i
(i
K21 K22
0 Zxi
Zxi
B '0i;1(x) '0i;1(x) dx
'0i(x) '0i;1(x) dx
B
B xi;1
xi;1
B
= B Zxi
B
Zxi
B
B
0
0
@
'i;1(x) 'i(x) dx
'0i (x) '0i(x) dx
xi;1

xi;1

1
C
C
C
C
C:
C
C
C
A

(3
...

Auf analoge Weise erhalten wir fur die rechte Seite

Zb
a

f (x) vh(x) dx =
=
=

n
X Zxi
i=1 xi;1
n
X Zxi

f (x) vh(x) dx
f (x) (vi;1'i;1(x) + vi'i(x)) dx

i=1 xi;1
n
X
(vi;1 vi)f (i)
i=1

(3
...
17)

xi;1

de niert
...
h
...
zum
globalen Lastvektor
...
14) und (3
...
Es gilt

Zb
a

uh (x) vh(x) dx =
0

0

n
X
i=1

(vi;1 vi)K (i)

ui;1
ui

!

= vT Kh uh =
h

63

3
...
DER ELEMENTWEISE AUFBAU DER STEIFIGKEITSMATRIX UND DES LASTVEKTORS
(1)
K11 u0 +

(1)
K12 u1

] v0

(1)
(1)
(2)
+ K21 u0 + (K22 + K11 )u1 +

...


...


(n
(n
(n
+ (K22 ;1) + K11 ))un;1 + K12 )un ] vn;1
(n
K21 )un;1

+


...


...


(n
+ K22 )un ] vn

und

Zb
a

f (x) v(x) dx =

n
X

(vi;1 vi)f (i) = vT f h =
h

i=1
f1(1) v0

(2)
(1)
(2)
+ (f1(2) + f2(1)) v1 + : : : + (fn;1 + fn ) vn;1 + fn vn

mit vh = (v0 v1 : : : vn)T , uh = (u0 u1 : : : un)T , f h = (f0 f1 : : : fn)T und
Kh = Kij ]nj=0
...

Zur Berechnung der globalen Stei gkeitsmatrix und des globalen Lastvektors fuhren
wir den folgenden Algorithmus durch:
Zuerst stellen wir den Zusammenhang zwischen der globalen und der lokalen Knotennumerierung her (siehe Abbildung 3
...
Dieser Zusammenhang wird in einer entsprechenden Zuordnungstabelle gespeichert
...


...


n

globale Knotennummer des Knotens mit
der lokalen Knotennummer
1
2
0
1
1
2
2
3

...


...


...

n;1
n

Tabelle 3
...
GRUNDPRINZIPIEN DER FINITE{ELEMENTE{METHODE: EIN 1D{BEISPIEL

'i;1(x)
A
A

x0

x1

'i(x)
A
A

A

A
A
A

x2

A
A

A
A

xi;1 xi xi+1
i;1 i i+1

globale Knotennummern:
0 1 2
lokale Knotennummern:
1 2
1 2

xn;1 xn
n;1 n

1

2
1

2

1

2

Abbildung 3
...
Danach bauen wir elementweise, d
...
Element
fur Element, die Elementstei gkeitsmatrizen K (i) und die Elementlastvektoren f (i) in
Kh bzw
...


K (i)

!

(i)
(i)
= K11) K12)
(i
(i
K21 K22

1

2

1
2

f (i) =

f1(i)
f2(i)

! 1
2

lokale Knotennummern

Nach dem Einbau von K (1) und f (1) erhalten wir
0
1
2

...

n

0 (1) (1)
B K11 K12
B (1) (1)
B K21 K22
B
B
B 0
0
B
B
B
...

B
...

B

...
C

...
C

...
C
A
0
0
2

n

und

0 (1)
B f1
B (1)
B f2
B
B
B 0
B
B
B
...

B
@
0

1
C
C
C
C
C
C:
C
C
C
C
C
A

65

3
...
DER ELEMENTWEISE AUFBAU DER STEIFIGKEITSMATRIX UND DES LASTVEKTORS

Der Einbau von K (2) und f (2) liefert
0 (1)
1
(1)
K11
K12
0 0
0C
B (1) (1) (2) (2)
B K K +K K 0
B 21 22
C
0C
11
12
B
C
B
C
(2)
(2) 0
B 0
K21
K22
0C
B
C
B
C
B 0
B
C
0
0 0
0C
B
C
B
...


...


...


...


...
C
@
A
0
0
0 0
0

und

0
f1(1)
B (1) (2)
B f +f
B 2
1
B
B
B f2(2)
B
B
B
B
0
B
B

...

B

...
f (i), i = 1 2 : : : n, entstehen die
Stei gkeitsmatrix
0 (1)
K (1)
0
0
0
B K11 (1) 12 (2)
B K (1) K + K
(2)
B 21 22 11
K12
0
0
B
B 0
(2)
(2)
(3)
(3)
B
K21
K22 + K11 K12
0
B
B
...


...


...

Kh = B
...


...

B
B (n;2) (n;1)
B f2 + f1
B
B (n;1) (n)
B f +f
B 2
1
@
f2(n)

globale

1
C
C
C
C
C
C
C
C
C
C
C
C
0 C
C
C
(n) C
K12 C
A
0
0
0

...


...

Einbau der Randbedingungen
Im Beispiel aus dem Abschnitt 3
...
Art bei x = a :
u(a) = ga
(b) Randbedingungen 3
...
GRUNDPRINZIPIEN DER FINITE{ELEMENTE{METHODE: EIN 1D{BEISPIEL

Mit den Beziehungen (3
...
16) konnen wir die Ersatzaufgabe (3
...
18)

mit uh = (ga u1 : : : un)T und vh = (0 v1 : : : vn)T gilt
...
18) aquivalent zu
= vT f~h + bgb vn fur alle vh 2 Rn :
(3
...

i
Da die Identitat (3
...

C u2 C B 2
...
C
B
...


...


...

B CB

...
C=B
C

...
CB
B
B
B0
C C Bf2(n;2) +f1(n;1)C
C
(n
(n
(n
(n
B
0
K21 ;2) K22 ;2) + K11 ;1) K12 ;1)
0 C
B
B CB
B0
C un;2C Bf (n;1) +f (n) C
C
B C @2
(n
(n
(n
(n
@
A
0
0
K21 ;1)
K22 ;1) + K11 ) K12 ) A
Bu A
1
@ n;1
(n)
(n)
(n)
0
0
0
0
K
K +
f + g

vT Kh uh +
h ~

b un vn

21

b

22

un

2

b b

gelten
...
Damit erhalten wir schlie lich das Finite{Elemente{Gleichungssystem
0 K(1) +K(2) K(2)
10 u1 1 0 f (1) +f (2) ;K(1) ga1
0
0
0
22
11
12
2
1
21
B K21 K22 +K11 K12
(2)
(2)
(3)
(3)
(2)
(3)
C
0
0 CB u2 C B f2 + f1
B
CB C B
C
B
CB
...


...
C B
B 0
C

...


...

B
C

...
CB
...

CB C=B
B
...

(n
(n
(n
(n
0
K21 ;2) K22 ;2) + K11 ;1) K12 ;1)
0 CB n;2 C B 2
B
CB C B (n;1) 1(n) C
C
B 0
CB C B
C
@
+f
A
0
K (n;1) K (n;1) + K (n) K (n) A@ u A @ f
0

0

21

0

22

(n
K21 )

11

n;1

12

(n
K22 ) +

b

un

2

1

(
f2n) + b gb

in dem die Randbedingungen u(a) = ga und ;u0(b) = b (u(b) ; gb ) eingearbeitet
sind
...
Zerlegung des Intervalls (a b) in nite Elemente (xi;1 xi), i = 1 2 : : : n
...
Fur alle Elemente (xi;1 xi)
{ Berechnung der Elementstei gkeitsmatrizen K (i), i = 1 2 : : : n
{ Berechnung der Elementlastvektoren f (i), i = 1 2 : : : n
{ Einbau von K (i) bzw
...
den globalen
Lastvektor
3
...
5
...
5 1D{LAGRANGE{Elemente hoherer Ordnung
Die im weiteren eingefuhrten Elemente hei en LAGRANGE{Elemente, weil zu ihrer
De nition LAGRANGE{Interpolationspolynome genutzt werden
...
Zur De nition der LAGRANGE{Elemente fuhren wir die
folgenden Schritte durch:
1
...
20)

Abbildung 3
...

2
...
9: Unterteilung des Basiselements b
und de nieren fur jeden Knoten k , k = 1 2 : : : m + 1, das Polynom
m+1
Y
j=1
j 6= k

( ; j ):

(3
...
Aufstellen des LAGRANGE{Interpolationspolynoms
Durch eine Normierung der Polynome (3
...


(3
...
GRUNDPRINZIPIEN DER FINITE{ELEMENTE{METHODE: EIN 1D{BEISPIEL

Im Fall m = 1 ergeben sich aus den Beziehungen (3
...
23)
Die Ansatzfunktionen uber den Elementen xi;1 xi] werden mittels einer Koordinatentransformation aus den Formfunktionen (3
...
Die so erhaltenen Ansatzfunktionen sind die im Abschnitt 3
...
5)
...
10: Bestimmung der linearen Ansatzfunktionen

Bemerkung 3
...

{ Der Einsatz dieser Ansatzfunktionen bei der Diskretisierung des Randwertproblems (3
...
3)
...
22) die quadratischen Formfunktionen
b 1( ) = 2 2 ; 3 + 1 b 2( ) = ;4 2 + 4 und b 3( ) = 2 2 ; :
(3
...
20) aus den Formfunktionen (3
...


69

3
...
1D{LAGRANGE{ELEMENTE HOHERER ORDNUNG

6

1

b1

b2

0

6

B
M
B
B
B
N
B
B

'2i;1

'2i;2

r

r

x0
x1
x2
p0 p1 p2 p3 p4

xi
p2i

;

1
2

;

-

1

1
2

x = x( )
1

b3

r
p2i

;

1

= (x)

'2i

xi
p2i

xn
p2n

1

;

;

2

r
p2n

;

-

1

xn x
p2n

Abbildung 3
...
Ansatzfunktionen

Bemerkung 3
...

{ Die Funktionen 'j (x), j = 0 1 : : : 2n, sind linear unabhangig
...
2 gestellten Randbedingungen die Struktur

0
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
B
@

0

ppp
0

ppp

p

1
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
A

70

KAPITEL 3
...
25)
9
b 4( ) = 2 ( ; 1 )( ; 2 ) :
b 3( ) = ; 27 ( ; 1 )( ; 1)
2
3
3
3
6

1

b1

b2

b3

b4

-

0

1
3

1

2
3

Abbildung 3
...
5 Die Ansatzfunktionen uber den Elementen xi;1 xi] werden wiederum durch eine Koordinatentransformation gema der Vorschrift (3
...
25) bestimmt
...

{ Diese Ansatzfunktionen sind linear unabhangig
...

In der Praxis werden auch Kombinationen von linearen, quadratischen und kubischen
Elementen auf unregelma igen Gittern genutzt, z
...


'2(x)

x0
p0

x1
p1

x2
p2

r
p3

'8(x)

x3
p4

r
p5

r

x4
x5
p6 p7 p8

r

r

p9

p10

r r

-

x6
x7 x
p11 p12 p13 p14

Abbildung 3
...
6
...
6 Au osung des Finite{Elemente{Gleichungssystems
Wir betrachten ein tridiagonales Gleichungssystems, das z
...
bei der FE{Diskretisierung von 1D{Randwertaufgaben 2
...
3)
...
3
...
CB
...
C
B
...


...


...
C B
...
CB
...
C
B
CB
C B
C
B 0
0 ;an;2 cn;2 ;bn;2
0 C B un;2 C B fn;2 C
B
CB
C B
C
B 0
0
0 ;an;1 cn;1 ;bn;1 C B un;1 C B fn;1 C
@
A@
A @
A
0
0
0
0 ;an
cn
un
fn
oder in der aquivalenten Schreibweise
c1u1 ; b1u2 = f1
; ai ui;1 + ci ui
; bi ui+1 = fi
; an un;1 + cn un = fn

i = 2 3 ::: n ; 1

Eliminationsschritte (LU{Faktorisierung)

1
...


3
...


...

(n ; 1)
...
GRUNDPRINZIPIEN DER FINITE{ELEMENTE{METHODE: EIN 1D{BEISPIEL

Infolge der Eliminationsschritte haben wir das tridiagonale Gleichungssystem in ein
Gleichungsystem mit einer Systemmatrix in Dreiecksgestalt uberfuhrt, d
...

b1
f1
u1 ; 2u2
= 2
2=
2=
c1
c1

ui ;

i+1 ui+1

un

=

=

= c ;bia
i
i
fi + ai
i+1 =
ci ; ai

i+1

i+1

n+1

n+1

i
i
i

i = 2 3 ::: n ; 1

+a
= fn ; a n n
...


...
C = B
...
4
...
C B
...
CB
B
C B
C

...
26)

Ruckwartseinsetzen

Das Gleichungsystem (3
...


un =
ui =

n+1
i+1

+

i+1 ui+1

i = n ; 1 n ;2 ::: 1

(3
...


73

3
...
AUFLOSUNG DES FINITE{ELEMENTE{GLEICHUNGSSYSTEMS

Speicherplatzbedarf
Die Systemmatrix, die rechte Seite des Gleichungssystems sowie die bei der Elimination
entstehende Dreiecksmatrix konnen wie folgt abgespeichert werden
...


...


...


...


...


...


...


...


...


...


...

an;1 cn;1 ; an;1 n;1 n
an
cn ; an n
0

2
3
4


...


...


...

un;1
un

Ist die Systemmatrix symmetrisch, dann entfallt der Speicherplatzbedarf fur die ai,
i = 2 3 : : : n
...
3n Speicherplatze benotigt, d
...
auch der
Speicherplatzbedarf ist proportional zur Anzahl der Unbekannten des zu losenden FE{
Gleichungssystems
...
1 werden hinreichende Bedingungen fur die Durchfuhrbarkeit und Stabilitat
des beschriebenen Au osungsalgorithmus formuliert
...
1

Voraussetzungen: jc1j > 0 jcn j > 0
jai j > 0 jbi j > 0 fur alle i = 2 3 : : : n ; 1
jci j jaij + jbi j fur alle i = 2 3 : : : n ; 1
jc1 j jb1j jcn j janj
wobei fur wenigstens eine der Ungleichungen jcij jaij + jbij
jc1 j jb1j bzw
...
ci ; ai i 6= 0 fur alle i = 2 3 : : : n (Durchfuhrbarkeit)
2
...

Die Bedingung ci ; ai i 6= 0 sichert die Berechenbarkeit der Gro en i+1 und i+1 in
den Eliminationsschritten und aus j ij 1 folgt, da Rundungsfehler, die in einem

74

KAPITEL 3
...
Sei
zum Beispiel beim Ruckwartseinsetzen fur ein i = i0 anstelle von ui0 ein Wert ui0 =
~
ui0 + i0 mit dem Fehler i0 berechnet worden
...
h
...
Unter der Bedingung j i0 j 1 erhalten wir fur den Fehler i0;1
i0
j i0 ;1 j = jui0 ;1 ; ui0 ;1 j = j i0 i0 j j i0 j :
~

Bemerkung 3
...
In unserem Beispiel (siehe Abschnitt 3
...
1
erfullt:
2
1
jhj > jhj
(jc1j > jb1j jc1j > 0)
2j
1j + j1j
jh h
(jcij jaij + jbij jaij > 0 jbij > 0 i = 2 3 : : : n ; 1)
jh
1 + j > j1j
jh b
(jcn j > janj jcn j > 0):
h

2
...

Bemerkung 3
...

C U = B
...

B
C
B
...


3
...
Dabei bezeichnen u die exakte Losung der Variationsformulierung des Randwertproblems und uh die Finite{Elemente{Naherungslosung
...
2 lautet (siehe
auch die Beziehung (3
...
28)
mit
Zb
Zb
0
0
a(u v) = u (x)v (x) dx + bu(b)v(b) hF vi = f (x)v(x) dx + bgb v(b)
a

a

Vg = fu 2 H 1(a b) : u(a) = gag und V0 = fv 2 H 1(a b) : v(a) = 0g
gilt
...
7
...
8 Der Ausdruck a(: :) ist eine Bilinearform, d
...

a( 1u1 + 2u2 v) =

1a(u1

v) + 2a(u2 v)

a(u 1v1 + 2v2) =

fur alle u1 u2 2 Vg fur alle v 2 V0
und fur alle 1 2 2 R1
fur alle u 2 Vg fur alle v1 v2 2 V0
und fur alle 1 2 2 R1

1a(u v1) + 2 a(u v2)

Finite{Elemente{Naherungslosung uh

Die FE{Naherungslosung erhalten wir aus der Ersatzaufgabe:
Gesucht ist uh 2 Vgh , so da

a(uh vh) = hF vh i
mit

a(uh vh) =
hF vh i

=

Zb
a

Zb
a

fur alle vh 2 V0h

0
u0h(x)vh(x) dx +

f (x)vh(x) dx +

(3
...
Die Bestimmung der Losung dieser Aufgabe erfolgt durch das Losen des linearen
FE{Gleichungssystems
n
X
uj a('j 'i) = hF 'i i ; u0 a('0 'i)
i = 1 2 : : : n:
(3
...
3 und 3
...

Beurteilung des Fehlers e(x) = u(x) ; uh (x) 2 V0

Zur Beurteilung des Fehlers benotigen wir eine Norm k : k
...
kvk 0 , kvk = 0 genau dann, wenn v = 0
2
...
ku + vk kuk + kvk fur alle u v 2 V0

76

KAPITEL 3
...
C {Norm k : kC
ku ; uh kC := xmax] ju(x) ; uh (x)j
2ab

2
...
H 1{Norm k : k1 = k : k1 2 (a b)
v
uZb h
u
i
t
ku ; uh k1 := u (u(x) ; uh (x))2 + ((u(x) ; uh (x))0 )2 dx
a

Falls u 2 Vg quadratisch integrierbare Ableitungen bis zur Ordnung k + 1 (zumindest elementweise) besitzt, d
...
wenn u0, u00, : : :, u(k+1) 2 L2(xj;1 xj ) gilt, und nite
Elemente p{ter Ordnung (p = 1 linear, p = 2 quadratisch, p = 3 kubisch) verwendet
werden, dann erhalten wir folgende Diskretisierungsfehlerabschatzungen:
ku ; uh ks

cs hminfk pg+1;s

ku ; uh kC

hminfk pg+0:5

c1

s=0 1
:

(3
...


Bemerkung 3
...
h
...
32)
fur alle u v 2 V0:

Im weiteren wollen wir die Abschatzung (3
...


u(0) = g0
u0(1) = 0

in (0 1)
(3
...
7
...
33) lautet:
Gesucht ist u 2 Vg , so da

a(u v) = hF vi
mit

a(u v) =

Z1

u (x)v (x) dx
0

0
Vg = fu 2 H 1(0

fur alle v 2 V0

hF v i =

0

Z1
0

(3
...

Die Finite{Elemente{Naherungslosung ist somit Losung der Aufgabe:
Gesucht ist die Funktion uh 2 Vgh , die der Beziehung

a(uh vh) = hF vh i
mit

a(uh vh) =

n

Z1
0

fur alle vh 2 V0h

uh(x)vh(x) dx
0

Vgh = uh : uh (x) =

n
V0h = vh : vh(x) =
genugt
...
32)
...
h
...
36)

1 Z1(v0(x))2 dx + 1 Z1 (v0(x))2 dx
= 2
20
0
1 Z1(v0(x))2 dx + 1 Z1(v(x))2 dx
20
2c2 0
F
(
)
1 1 kvk2
min 2 2c2
1
F

o

(3
...
Mit der Schwarzschen Ungleichung

(3
...
GRUNDPRINZIPIEN DER FINITE{ELEMENTE{METHODE: EIN 1D{BEISPIEL

Z1
0

v
uZ1
u
u (u0(x))2 dx
t

u0(x)v0(x) dx

0

v1
uZ
u
u (v0(x))2 dx
t
0

(3
...
39)

0

= kuk1 kvk1
d
...
2 = 1
...
40)
Subtrahieren wir die Gleichung (3
...
40), dann ergibt sich
a(u ; uh vh) = 0 fur alle vh 2 V0h:
(3
...
41)
2
1 ku ; uh k1

=
=
=
=

a(u ; uh u ; uh )
a(u ; uh u) ; a(u ; uh uh)
a(u ; uh u) ; a(u ; uh wh ; vh)
a(u ; uh u) ; a(u ; uh wh )
a(u ; uh u ; wh )
2 ku ; uh k1 ku ; wh k1

(3
...
h
...
43)

Setzen wir in (3
...
44)

3
...
DISKRETISIERUNGSFEHLERABSCHATZUNGEN

ein und beachten wir, da in unserer Modellaufgabe
gilt, so ergibt sich
ku ; uh k1

min

n1
1

1
2 2c2
F

1

= min

n1

1
2 2c2
F

o

79
sowie

o ku ; Inth (u)k1 :

2

=1

(3
...

Dabei erhalten wir
k(u ; wh ) k2
0
0

=
=

Z1
0

(u (x) ; wh
0

0

n
X Zxj
j =1 xj;1

(x))2 dx

=

n
X Zxj
j =1 xj;1

0
(u0(x) ; wh (x))2 dx

(z0(x))2 dx

mit z(x) = u(x) ; wh(x)
...
44)), gilt
o enbar

z(xj ) = u(xj ) ; wh (xj ) = 0

fur alle j = 0 1 : : : n:

Weiter folgt

Zxj
xj;1

(z (x))2 dx
0

=

Zxj
xj;1

1
z (x) ; h
0

Zxj

2

z 0( ) d
xj;1
| {z }

dx

z (xj ) ; z (xj 1) = 0
;

Zxj 1 Zxj
0
=
; 0 }
|
h xj;1 (z (xx) {z z ( )) d
xj;1
R
= z ( )d

2

dx

00

Zxj 1 Zxj
Zx
=
z00( ) d d
h xj;1 1
x j ;1
=

Zxj
xj;1

x
1 Z j 1 Zx z00( ) d d
h2 xj;1

2

dx
2

dx :

Die zweimalige Anwendung der Schwarzschen Ungleichung liefert

Zxj
xj;1

(z (x))2 dx
0

Zxj
xj;1

x
x
1 Z j 12 d Z j Zx z00( ) d
h2 xj;1
xj;1

2

d

dx

80

KAPITEL 3
...
45) und der Friedrichsschen Ungleichung (3
...
7
...
10 Falls der FE{Raum V0h die stuckweise linearen Funktionen ent-

halt, dann ist bei Randwertproblemen der Gestalt (3
...
Dies wollen wir kurz begrunden
...
46)

f (x)v(x) dx + gbv(b)

Vg = fu(x) 2 H 1(a b) : u(a) = gag und V0 = fv(x) 2 H 1(a b) : v(a) = 0g
gilt
...
Die Elemente der Funktionenmengen Vgh und V0h seien stetige Funktionen, die uber den niten Elementen xj;1 xj ] Polynome p{ten Grades sind
...
41) erhalten wir die
Relation

a(u ; uh vh) = 0
fur alle vh 2 V0h , d
...
fur die Bilinearform in der Aufgabe (3
...
47)

Wir wahlen als Testfunktion vh 2 V0h zunachst eine Funktion, welche die folgenden
Eigenschaften besitzt:
0
1
...
h
...
GRUNDPRINZIPIEN DER FINITE{ELEMENTE{METHODE: EIN 1D{BEISPIEL

2
...
vh(x) = v1 fur alle x 2 (x1 xn], d
...
vh(x) = 0 fur alle x 2 (x1 xn]
...
47)
Zx1
Zx1
0
0
0
u (x)vh(x) dx = u0h(x)vh(x) dx:
x0

x0

0
Dies ist wegen vh(x) = const: in x0 x1] aquivalent zu
Zx1
Zx1
0
u (x) dx = u0h(x) dx

x0

x0

d
...


u(x1) ; u(x0) = uh(x1) ; uh(x0) :
Wegen der vorgegebenen Randbedingung u(x0) = uh(x0) = ga erhalten wir schlie lich
u(x1) = uh(x1) :
Unter Ausnutzung der Beziehungen u(x0) = uh(x0) sowie u(x1) = uh(x1) konnen
wir beweisen, da auch u(x2) = uh(x2) gilt
...
Die Funktion vh ist stuckweise linear, d
...
linear in den Intervallen x0 x1] sowie
x1 x2],
2
...
vh = 0 fur alle x 2 x2 xn]
...
47) die Form

Zx1
x0

Zx2

Zx1

Zx2

0
u (x)vh(x) dx + u (x)vh(x) dx = uh(x)vh(x) dx + u0h(x)vh(x) dx:
0

0

x1

0

0

x0

0

0

x1

(3
...
Fur alle weiteren i = 3 4 : : : n kann auf vollig analoge Weise die Gleichheit
u(xi) = uh(xi) bewiesen werden
...
8
...
8 Zwei Anwendungsbeispiele
Als erstes Beispiel betrachten wir die stationare Warmeleitung in einem Stab, bestehend aus zwei Materialien
...

Dieses physikalische Problem wird durch die folgende Randwertaufgabe beschrieben
(siehe auch Abschnitt 2
...
1):

x=0
Randbedingung
x 2 (0 0:5) Di erentialgleichung
x = 0:5
Interfacebedingung
x 2 (0:5 1) Di erentialgleichung
x=1
Randbedingung

u(0)
;(50 u0 )0
u(0:5 ; 0)
;50 u0 (0:5 ; 0)
;(371 u0 )0
u(1)

=
=
=
=
=
=

273:15
0
u(0:5 + 0)
;371 u0 (0:5 + 0)
0
323:15

(3
...
49) u(x) = u(x) + 273:15 , dann erhalten wir die aquivalente Aufgabe:

x=0
Randbedingung
x 2 (0 0:5) Di erentialgleichung
x = 0:5
Interfacebedingung
x 2 (0:5 1) Di erentialgleichung
x=1
Randbedingung

u(0)
;(50 u0 )0
u(0:5 ; 0)
;50 u0 (0:5 ; 0)
;(371 u0 )0
u(1)

=
=
=
=
=
=

0
0
u(0:5 + 0)
;371 u0 (0:5 + 0)
0
50

(3
...
Da die rechten Seiten in den
Di erentialgleichungen identisch Null sowie die Warmeleitkoe zienten konstant sind,
wird folglich eine Funktion gesucht, deren 2
...
Wir vermuten,
da das Temperaturfeld u(x) in den Teilgebieten 1 und 2 jeweils eine lineare Funktion
ist und verwenden deshalb den folgenden Losungsansatz:

u(x) = cx + d
u(x) = ex + f

in
in

1

2:

(3
...
GRUNDPRINZIPIEN DER FINITE{ELEMENTE{METHODE: EIN 1D{BEISPIEL

u(0) = 0c + d = 0
u(1) = 1e + f = 50
Interfacebedingungen: 0:5c + d = 0:5e + f
;50c = ;371e
Randbedingungen:

=)
=)
=)
=)

d=0
e + f = 50
0:5c = 0:5e + f
c = 371 e :
50

Damit ergibt sich

9
e + f = 50 =
0:5e + f = 3:71e

=)

e + f = 50
;3:21e + f = 0 :

Die Losungen dieses Gleichungssystems sind:

e = 5000
421

und

f = 16050
421

so da c = 37100 gilt
...
50) hat folglich die Losung:

u(0) = 0
u(x) = 37100 x
421
u(0:5) = 18550
421
u(x) = 5000 x + 16050
421
421
u(1) = 50 :

fur x 2 (0 0:5)
(3
...
Deshalb stellen wir zuerst die Variationsformulierung auf
...
2, d
...

0
0
Z :5
Z1
Z:5
Z1
0
0
0
0
;(50u (x)) v (x) dx + ;(371u (x)) v (x) dx = 0 v (x) dx + 0 v (x) dx :
0

0:5

0

0:5

Mit der Formel der partiellen Integration folgt
0
Z :5
Z1
0:5
1
0
0
0
50u (x)v (x) dx ; 50u (x)v(x) + 371u0 (x)v0(x) dx ; 371u0 (x)v(x) = 0
0

0

0:5

0:5

85

3
...
ZWEI ANWENDUNGSBEISPIELE
0
Z:5
Z1
0
0
50u (x)v (x) dx + 371u0(x)v0(x) dx ; 50u0 (0:5)v(0:5) + 50u0(0)v(0)
0:5

0

; 371u0 (1)v (1) + 371u0 (0:5)v (0:5) = 0 :

Beachten wir noch die Interfacebedingung ;50u0(0:5 ; 0) = ;371u0(0:5 + 0) sowie die
Bedingung v(0) = v(1) = 0 , so erhalten wir die Variationsformulierung:
Gesucht ist u(x) 2 Vg = fu(x) 2 H 1(0 1) : u(0) = 0 u(1) = 50g so da

Z1
0

mit

(x)u0(x)v0(x) dx = 0 fur alle v 2 V0 = fv(x) 2 H 1(0 1) : v(0) = v(1) = 0g
(3
...

Diese Variationsformulierung ist der Ausgangspunkt einer Finite{Elemente{Diskretisierung
...
Um auch den Fall einer nicht aquidistanten
Unterteilung erlautern zu konnen, wahlen wir die folgende Zerlegung:
0 1] =

4

i=1

xi;1 xi]

(3
...
h
...
Anfangspunkt zweier
Teilintervalle gewahlt wird
...
Deshalb wurde bereits bei einer FE{Diskretisierung
mit den zwei niten Elementen (0 0:5) und (0:5 1) sowie stuckweise linearen Ansatzfunktionen die Losung des Randwertproblems exakt approximiert
...

Die stuckweise linearen Ansatzfunktionen sind durch
8
>0
x xi;1
>
>
>
> x ; xi;1
>
> x ;x
xi;1 < x xi
< i i;1
'i(x) = >
(3
...
GRUNDPRINZIPIEN DER FINITE{ELEMENTE{METHODE: EIN 1D{BEISPIEL

de niert
...
1)
8
>0
x xi;1
>
>
>
>
>
> x ;1x
xi;1 < x xi
< i i;1
0
'i(x) = >
(3
...

Gesucht ist uh(x) 2 Vgh , so da
Z1
0
(x)u0h (x)vh(x) dx = 0 fur alle vh 2 V0h
(3
...

Analog zu der im Abschnitt 3
...
58)
B Zi
C
Zxi
B
C
B
@
(x)'0i;1(x)'0i(x) dx
(x)'0i (x)'0i(x) dx C
A
xi;1

xi;1

0 Ri
1
x
f (x)'i;1(x) dx C
B xi;1
B
C
B
C
C:
f (i) = B
B x
C
B Ri
C
@
f (x)'i(x) dx A
x
i;1

(3
...
8
...
55) und (3
...
2: Elementstei gkeitsmatrizen und Elementlastvektoren
Bevor wir die globale Stei gkeitsmatrix und den globalen Lastvektor assemblieren,
stellen wir noch die Zuordnungstabelle zwischen der lokalen und der globalen Knotennumerierung auf
...
3: Zuordnungstabelle zwischen globaler und lokaler Knotennumerierung
Die Assemblierung der Elementstei gkeitsmatrizen und der Elementlastvektoren ergibt
die globale Stei gkeitsmatrix Kh und den globalen Lastvektor f h ,
0 500 ; 500
001
0
0
0 1
3
3
B
C
B C
B 500 500 500
C
B C
B;
C
B0C
500
;2
0
0 C
B 3 3 + 2
B C
B
C
B C
B
C
B C
B
C
B C
Kh = B 0
0 C fh = B 0 C
; 500 500 + 3170 ; 3170
B
C
B C
2
2
2
2
B
C
B C
B
C
B C
B 0
B0C
3710
3170 + 3170 ; 3170 C
B
C
B C
0
; 2
2
3
3 C
B
B C
@
A
@ A
3710
3170
0
0
0
0
; 3
3

88

KAPITEL 3
...
Der Einbau der Randbedingungen u(0) = u0 = 0 und u(1) = u4 = 50 liefert das FE{Gleichungssystem

0 500 + 500
B 3 2
B 500
B ;
B
2
B
@
0

; 500
2
500
2

+ 3170
2

; 3710
2

0
; 3170
2

+ 3170
3

3170
2

10 u
CB 1
CB
C B u2
CB
CB
A@
u3

1 0
C B
C B
C=B
C B
C B
A @

0
0
50 3170
3

1
C
C
C
C
C
A

~
w
w
w
0 1250 ;250
0
B 3
B
B ;250 2105 ;1855
B
B
@
0 ;1855 9275
3

10 u
CB 1
CB
C B u2
CB
CB
A@
u3

1 0 0 1
C B
C
C B
C=B 0 C
C
C B
C
C B
C
A @
A
185500
3

Zur Losung dieses Gleichungssystems verwenden wir den im Abschnitt 3
...
Wir erhalten
2

=3
5

2=0

3
u1 ; 5 u2 = 0

3

1855
= 2105 ; 3 250 = 371
391
5

3=0

u2 ; 371 u3 = 0
391

4=

185500
3

9275 ; 1855 371
3
391

= 19550
421

und somit

u3 = 19550
421

u2 = 371 19550 = 18550
391 421
421

3
u1 = 5 18550 = 11130 :
421
421

Die FE{Losung ist folglich

uh(x) =

4
X

i=0

ui'i(x)

mit

u = ui]4=0
i

u = 0 11130 18550 19550 50
421 421 421

T

(0 26:44 44:06 46:44 50)T :

89

3
...
ZWEI ANWENDUNGSBEISPIELE

Zum Abschlu vergleichen wir die analytische Losung mit der FE{Losung
...
4: Analytische Losung und FE{Losung der betrachteten Randwertaufgabe
u 6
50

-

0
0:5
1 x
Abbildung 3
...
60)
;u0(2) = ;50 :
Die analytische Losung dieser Aufgabe ist u(x) = 1 + 2x ; 3x2 + x3 ; x4 + x5 ,
denn es gilt
u(0) = 1
u0(x) = 2 ; 6x + 3x2 ; 4x3 + 5x4 d
...
; u0(2) = ;50
u00(x) = ; 6 + 6x ; 12x2 + 20x3 :
Als Variationsformulierung erhalten wir (siehe auch Abschnitt 3
...


(3
...
GRUNDPRINZIPIEN DER FINITE{ELEMENTE{METHODE: EIN 1D{BEISPIEL

Bei der FE{Diskretisierung der Aufgabe (3
...
5 de nierten stuckweise quadratischen Ansatzfunktionen nutzen
...
h
...
15: Zerlegung des Intervalls (0 2) und die Zuordnung zwischen globaler
und lokaler Knotennumerierung
Die stuckweise quadratischen Ansatzfunktionen 'j (x), j = 0 1 : : : 6, de nieren wir,
indem wir die auf dem Basiselement b = (0 1) de nierten Formfunktionen (siehe
Abschnitt 3
...
Dazu verwenden wir die
Transformationsvorschrift
x
1
x( ) = hi + xi;1
(x) = x ; xi;1 = h (x ; xi;1)
i ; xi;1
i
und die Formfunktionen

b 1( ) = 2 2 ; 3 + 1

b 2( ) = ; 4 2 + 4

b 3( ) = 2 2 ; :

Wir erhalten

8b
> 1( (x)) = b 1( h11 (x))
<
'0(x) = >
:0

8
>0
>
>
> b ( 1 (x ; x ))
> 3
<
i;1
'2i(x) = > hi
> b 1( hi1+1 (x ; xi))
>
>
>0
:

fur 0 x < 2
3
2
fur 3 x
fur x xi;1
fur xi;1 x < xi
fur xi x < xi+1
fur xi+1 x

i=1 2

91

3
...
ZWEI ANWENDUNGSBEISPIELE

8
>0
<
= >
4
: b 3( h1 (x ; 3 ))

'6(x)

fur x
fur

3

8
>0
>
>
<
'2i;1(x) = > b 2( h1i (x ; xi;1))
>
>0
:

4
3

4
3

x<2

fur x xi;1

i = 1 2 3:

fur xi;1 x < xi
fur xi x

Auf die im Abschnitt 3
...
62)
0

fur alle vh 2 V0h mit
Vgh = fuh (x) : uh(x) =

0

6
X

i=1

ui 'i(x) + '0(x)g und V0h = fvh (x) : vh (x) =

6
X

i=1

vi'i (x)g

gilt
...
4 beschriebenen
Vorgehensweise stellen wir das FE{Gleichungssystem auf
...
h
...
6)
" Z1
#
d'k (x( )) d (x) d'l(x( )) d (x) h d 2i
(i) =
K
d
dx
d
dx i
k l = 2i;2

0

=

" Z1
0

db 1 db 1 h d
d hi d hi i

#3
=1

" Z1 b
d
1
= h
d
i0

db
d

d

#3
=1

92

KAPITEL 3
...
Fur die 3 niten Elemente ergeben sich die
folgenden Elementstei gkeitsmatrizen und Elementlastvektoren:
nites Element
(0 2 )
3

2
(3 4)
3

( 4 2)
3

1
2

1
2

1
2

0 7
B
B ;8
B
@
1
0 7
B
B ;8
B
@
1
0 7
B
B ;8
B
@
1

K (i)

f (i)

1
C
C
C
A

;8

1
C
C
C
A

;8

1
C
C
C
A

0 146
B 405
B ; 544
B 135
@

1
C
C
C
A

1
C
C
C
A

0 ; 3302
1215
B 35216
B;
B 1215
@

1
C
C
C
A

1
16 ; 8
;8 7
1
16 ; 8
;8 7
;8

1
16 ; 8
;8 7

0
B
B
B
@

818
1215
2384
1215
278
1215

; 1154
405

; 15362
1215

Tabelle 3
...
6
Elementnummer
1
2
3

globale Knotennummer des Knotens
mit der lokalen Knotennummer
1
2
3
0
1
2
2
3
4
4
5
6

Tabelle 3
...
8
...
5
erhalten wir die globale Stei gkeitsmatrix Kh und den globalen Lastvektor f h,
0
1
0
1
818
7 ;8 1
0
0
0 0 C
1215
B
B
C
B ; 8 16 ; 8 0
B
C
2384
0
0 0 C
B
C
B
B
C
B 278 1215 146 C
C
B 1 ;8 7+7 ;8 1
B
B
C
B 1215 + 405 C
C
0 0 C
C
B
C
B
C
1 B 0 0 ; 8 16 ; 8 0 0 C f = B
544
B
C
B
C
Kh = 2 B
; 135
C
h B
B
C
B 1154 3302 C
B 0 0
B ; 405 ; 1215 C
C
1 ;8 7+7 ;8 1 C
B
C
B
C
B
C
B
B 0 0
C
B ; 35216 C
C
B
B
C
0
0 ; 8 16 ; 8 C
1215
@
A
@
A
15362
0 0
0
0
1 ;8 7
; 1215
ohne Berucksichtigung der Randbedingungen
...


analytische Losung u(x) FE{Losung uh (x) Fehler ju(x) ; uh (x)j

x=0
x=1
3
x=2
3
x=1
4
x=3

x=5
3
x=2

1
331
243
299
243
1
427
243
1403
243
17

1
551
405
299
243
401
405
427
243
259
45
17

0
2
1215
0
4
405
0
22
1215
0

Tabelle 3
...
GRUNDPRINZIPIEN DER FINITE{ELEMENTE{METHODE: EIN 1D{BEISPIEL

In diesem Beispiel ist die FE{Losung eine sehr gute Naherung fur die analytische
Losung
...
7 ist ersichtlich, da in den Intervallenden der Elemente
xi;1 xi] der Fehler u(x) ; uh (x) gleich Null ist
...
10 bewiesen
...
Dies ist darin begrundet, da die analytischen Losungen eine besonders
einfache Struktur hatten, d
...
da sie Polynome niedrigen Grades waren
...
9 vorgestellten Beispiele zeigen aber, da es i
...
notwendig ist, die Diskretisierungsschrittweite bzw
...
Das ist auch der typische Fall in praktischen Anwendungen
...
9 Das Programm FEM1D
Das Programm FEM1D kann als Lehrprogramm zur Demonstration des algorithmischen Ablaufs der Finite{Elemente{Diskretisierung elliptischer Randwertprobleme in
eindimensionalen Gebieten eingesetzt werden (siehe auch 79])
...
Funf verschiedene im Programm implementierte Beispiele, deren analytische Losung bekannt
ist, erlauben es, den Ein u der FE{Diskretisierung auf die Approximationsgute zu
untersuchen
...

Mit dem Programm FEM1D sind stationare Warmeleitprobleme in Gebieten, die aus
maximal funf verschiedenen Materialien bestehen, losbar
...
1
...
h
...

In den Unstetigkeitspunkten zj der Koe zientenfunktion k(x) sind die Interfacebedingungen:
u(zj ; 0) = u(zj + 0)
;k (zj ; 0)u0 (zj ; 0) = ;k (zj + 0)u0 (zj + 0)

gegeben, und in den Punkten x = j , in denen eine Punktquelle bzw
...
3)
...


sowie

u(b) = gb oder ; k(b)u0(b) = gb

95

3
...
DAS PROGRAMM FEM1D

Im weiteren stellen wir zwei Beispiele vor, die im Programm FEM1D implementiert
sind
...
63)

Die exakte Losung dieser Aufgabe ist u(x) = 1 + 2x ; 3x2
...
63) wahlen wir die folgenden
drei Diskretisierungen:
1
...
Variante

9 nite Elemente mit stuckweise linearen Ansatzfunktionen
u

u

6

1

0

6

1

-

1
...
Variante

1x

Abbildung 3
...


96

KAPITEL 3
...
Variante
1 nites Element mit quadratischen Anastzfunktionen

Die dritte Variante liefert eine exakte Approximation der analytischen Losung, da
diese eine quadratische Funktion ist
...

; (k u0(x))0 + c % v u0(x) = 0

fur x 2 (0 1)

(3
...

Die Aufgabe (3
...
65)
fur x = 0
fur x = 1
px ;
Die analytische Losung der Aufgabe (3
...

Wir untersuchen fur verschiedene Werte des Parameters p, welche FE{Diskretisierungen durchgefuhrt werden konnen, um eine gute Naherungslosung zu erhalten
...
Variante: p = 10, 10 gleichgro e Elemente mit linearen Ansatzfunktionen
u6

u(x) = 0
u(x) = 1

1

0

-

1x

1
...
17: Analytische und FE{Losung fur p = 10

97

3
...
DAS PROGRAMM FEM1D

Es zeigt sich, da die exakte Losung bei dieser Diskretisierung gut approximiert wird
...
Variante: p = 70, 10 gleichgro e Elemente mit linearen Ansatzfunktionen
3
...
Variante

6

-

1x

0

3
...
18: Analytische und FE{Losung bei verschiedener Feinheit der Diskretisierung
Die Naherungslosung bei der zweiten Variante weist ein stark oszillierendes Verhalten
auf
...
Die dritte
Variante liefert eine FE{Losung, welche die analytische Losung nur noch in der Nahe
des rechten Gebietsrandes schlecht approximiert
...
Variante: p = 70, 70 gleichgro e Elemente mit linearen Ansatzfunktionen

Diese Variante liefert eine sehr gute Approximation der exakten Losung
...
Wir zeigen mit der 5
...
Der Aufwand zur Generierung des FE{Gleichungssystems und zu seiner
Au osung ist auf Grund der niedrigen Anzahl von FE{Knoten naturlich wesentlich
geringer als bei der Diskretisierung mit 70 gleichgro en Elementen
...
GRUNDPRINZIPIEN DER FINITE{ELEMENTE{METHODE: EIN 1D{BEISPIEL

5
...
Variante

-

0

5
...
19: FE{Diskretisierungen mit gleichma iger und lokaler Netzverfeinerung
Eine weitere Moglichkeit zur Anpassung der FE{Diskretisierung an das Losungsverhalten ist die Erhohung des Polynomgrades der Ansatzfunktionen
...

6
...
9,0
...
Variante

Abbildung 3
...
Ordnung
4
...
1
...
Wie
die Beispiele aus den Abschnitten 1
...
2 und 1
...
3 zeigen, konnen auch elektrische und
magnetische Felder durch eine derartige Gleichung beschrieben werden
...

Klassische Formulierung (siehe auch Abschnitt 2
...
2)

Gesucht ist u(x) 2 C 2( ) \ C 1( ;2 ;3) \ C ( ), so da
@u ; @
@u = f (x)
fur alle x = (x1 x2) 2
; @ 1 (x)
2 (x)
@x1
@x1 @x2
@x2

u(x) = g1(x)

fur alle x 2 ;1

2
X

@u =
@u n = g (x)
fur alle x 2 ;2
i
2
@N i=1 @xi i
@u + (x)u(x) = (x)g (x) fur alle x 2 ;
3
3
@N
gilt mit ; = @ = ;1 ;2 ;3 , ;i \ ;j = fur i 6= j
...
= 500 K , g2 = 0 (Isolation),
@u = @u = ; @u = (g ; u),
I
3
@N I @n
@x2
= 5:6 W (m2K );1, g3 = 300 K
...
1)

100

KAPITEL 4
...
ORDNUNG

x2 6
0:8

(

II

(

I

;1
;2
;3
;In

{ Material 2)

{ Material 1)

-

0
1 x1
Abbildung 4
...
1) nur in den Teilgebieten I und II formulieren
...
Wir erhalten
das folgende Randwertproblem:
Gesucht ist das Temperaturfeld u(x) mit u(x) = uI (x) in I , u(x) = uII (x) in II ,
fur das
fur alle x = (x1 x2) 2 I
; @ I @uI ; @ I @uI = 0
@x1 @x1 @x2 @x2
; @ II @uII ; @ II @uII = 0
fur alle x = (x1 x2) 2 II
@x1
@x1
@x2
@x2

u(x) = g1(x)
@u = 0
@N
;

I

fur alle x 2 ;1
fur alle x 2 ;2

@uI + (x)u (x) =
I
@x2

(x)g3(x) fur alle x 2 ;3

uI (x) = uII (x)
I

@uI =
@x2

II

@uII
@x2

fur alle x 2 ;In

gilt mit ; = @ = ;1 ;2 ;3 , ;i \ ;j = fur i 6= j
...
2
...
2 Herleitung der Variationsformulierung
Die Raume L2( ) und H 1( )
Eine Funktion u(x) = u(x1 x2) ist Element des Raumes L2( ), wenn das Integral

Z

(u(x))2 dx

existiert und endlich ist
...
Seien die Funktionen u und w uber integrierbar, und es gelte
Z
Z
@'
u(x) @x dx = ; w(x)'(x)dx
i
fur jede stetig di erenzierbare Funktion '(x), die auf @ verschwindet
...

Der Raum H 1( ) umfa t alle die Funktionen u(x), die zum Raum L2( ) gehoren und
@u
deren verallgemeinerten Ableitungen @xi , i = 1 2, existieren sowie ebenfalls Elemente
des Raumes L2( ) sind (siehe auch 1, 86])
...
1 Zum Beispiel gehoren die bei der FEM verwendeten stetigen, stuckweise polynomialen Ansatzfunktionen zum Raum H 1( ) (siehe Abschnitt 4
...
3)
...
Wir multiplizieren die Di erentialgleichung (4
...
h
...
Wir wenden im Hauptteil (bei den 2: Ableitungen) die Formel der partiellen Integration
Z @w
Z @v
Z
@xi v dx = ; w @xi dx + @ wvni ds
an und erhalten
Z
@u @v + (x) @u @v dx ; Z
1 (x)
@x1 @x1 2 @x2 @x2
@

1

@u vn + @u vn ds
@x1 1 2 @x2 2
Z
= f (x)v(x) dx:

@u
Unter Beachtung der De nition der Konormalenableitung @N folgt
Z
@u @v + (x) @u @v dx ; Z @u v ds = Z f (x)v(x) dx:
1 (x)
@x1 @x1 2 @x2 @x2
@N
@

102

KAPITEL 4
...
ORDNUNG

3
...
Wir legen unter Beachtung der wesentlichen Randbedingungen die Menge der
zulassigen Funktionen, d
...
die Menge, in der wir die Losung suchen, fest:
Vg1 = fu(x) 2 H 1( ) : u(x) = g1(x) auf ;1g :
Die Schritte 1
...
fuhren zur Variationsformulierung der Aufgabe (4
...
2)

@u @v + (x) @u @v dx + Z (x)u(x)v(x) ds
@x1 @x1 2 @x2 @x2
;3
Z
Z
Z
hF v i = f (x)v (x) dx + g2 (x)v (x) ds +
(x)g3(x)v(x) ds
a(u v) =

Z

1 (x)

;2

;3

Vg1 = fu(x) 2 H 1( ) : u(x) = g1(x) auf ;1g
V0 = fv(x) 2 H 1( ) : v(x) = 0 auf ;1 g
erfullt ist
...
2 Die linke Seite in der Gleichung (4
...
h
...
2
...
3 Die Variationsformulierung (4
...
1)
...
B
...
Von der Losung u(x) der Aufgabe (4
...
Es mussen die verallgemeinerten Ableitungen @xi ,
i = 1 2, existieren und Elemente des Raumes L2( ) sein
...
2) die klassischen Glattheitsforderungen erfullt, d
...
u(x) 2
Vg1 \ (C 2( ) \ C 1( ;2 ;3 ) \ C ( )) , dann ist sie auch Losung der klassischen
Aufgabe (4
...

Bemerkung 4
...
Art
...

~
Die zur Aufgabe (4
...
3)
mit
Z
@w @v + (x) @w @v dx + Z (x)w(x)v(x) ds
a(w v) =
2
1 (x)
@x @x
@x @x
1

1

2

2

e
hF v i = hF v i ; a(~1 v )
g

;3

V0 = fv(x) 2 H 1( ) : v(x) = 0 auf ;1g
erfullt ist
...
2) ergibt sich somit aus u(x) = w(x) + g1(x)
...
5 Die Variationsformulierung (4
...
(4
...

Aquivalenz der Variationsformulierung zu einem Minimumproblem
Falls die Bilinearform a(: :) symmetrisch und positiv ist, d
...

a(w v) = a(v w) fur alle w v 2 V0 und

a(w w) > 0
fur alle w 2 V0 w 6= 0 ,
dann ist das Variationsproblem (4
...


J (w) = min J (v)
v2V0
1
e
J (v) = 2 a(v v) ; hF vi

(4
...
FEM FUR MEHRDIMENSIONALE RANDWERTPROBLEME 2
...
5)
1
e
e
= 1 a(w w) + 1 t a(w v) + a(v w)] + 2 t2a(v v) ; hF wi ; thF vi
2
2
1
e
= J (w) + t a(w v) ; hF vi] + 2 t2a(v v) :

(4
...
4)
Sei w 2 V0 die Losung des Variationsproblems (4
...
Aus (4
...

(4
...
3)
d
Sei w 2 V0 Minimalpunkt von J (:)
...
h
...
3)
...
3 Galerkin{Ritz{FEM
4
...
1 Galerkin{Verfahren

Ausgangspunkt: Variationsformulierung (siehe z
...
Abschnitt 4
...
6)

K
A
A

Menge der
Testfunktionen

Die beiden Funktionenmengen Vg1 und V0 sind Teilmengen der Menge der Grundfunktionen V
...
3
...


i2

Dabei ist
X
i 2 !h

X
i2

vipi(x) = 0 fur alle x 2 ;1

u ipi(x) = g1(x) bzw
...
Die Anzahl der Indizes in der Menge !h bezeichnen wir mit Nh und die der
Menge !h mit Nh
...
6) erhalten wir die Naherungsaufgabe
Gesucht ist uh 2 Vg1 h : a(uh vh) = hF vhi fur alle vh 2 V0h
~
w
w
X
X
w
w
uh =
uipi +
u ipi , vh = pk 2 V0h k 2 !h
w
w
i 2 !h
i2 h
w
w
w
a(: :) bilinear, hF :i linear
w
w
w
Gesucht ist uh = ui]i 2 !h 2 RNh :
X
X
uia(pi pk ) = hF pk i ;
u ia(pi pk )
i 2 !h

i2

fur alle k 2 !h

(4
...
8)

h

d
...
gesucht ist uh = ui]i 2 !h 2 RNh :
mit

Khuh = f h

(Galerkin{Gleichungssystem)

Kh = a(pi pk )]i k 2 !h
f h = hF pk i ; P u ia(pi pk )]k 2 !h 2 RNh
i2

h

(Stei gkeitsmatrix)
(Lastvektor)

(4
...
FEM FUR MEHRDIMENSIONALE RANDWERTPROBLEME 2
...
3
...
10)

Im Ritzschen Energiefunktional
1
e
J (v) = 2 a(v v) ; hF vi
setzen wir voraus, da a(: :) bilinear, positiv und symmetrisch ist
...
3
...

Idee des Ritz{Verfahrens

Gesucht ist uh = wh + g1h 2 Vg1 h mit J (wh) = vhmin0h J (vh)
~
2V
~
w
w
X
X
w
uh = wh + g1h =
~
uipi +
u ipi
w
w
i 2 !h
i2 h
w
w
w
w
@J (wh) = 0 fur alle k 2 !
w
w
h
w
@uk
w
Gesucht ist uh = ui]i 2 !h 2 RNh :
X
X
uia(pi pk ) = hF pk i ;
u ia(pi pk )
i 2 !h

i2

fur alle k 2 !h

(4
...
12)

h

d
...
gesucht ist uh = ui]i 2 !h 2 RNh :

Kh uh = f h

(Ritzsches Gleichungssystem)

mit

Kh = a(pi pk )]i k 2 !h
f h = hF pk i ; P u ia(pi pk )]k 2 !h 2 RNh
i2

h

Bemerkung 4
...
3
...
3
...
Galerkin{Verfahren zu verwenden, stammt
von R
...
B
...

Bei der Finite{Elemente{Methode wird das Gebiet in "kleine\, nichtuberlappende
Teilgebiete (r) (Dreiecke oder Vierecke) zerlegt
...
Als Ansatzfunktionen wahlen wir Funktionen, die nur
uber sehr wenigen Dreiecken bzw
...
Eine derartige
Ansatzfunktion ist in der Abbildung 4
...


pi(x)

s

K
A
A

A
A

Trager
supp(pi) :

s

A
A

xi { Knoten der zur Ansatzfunktion pi (x) "gehort\
pi (xj ) = ij
(! Knotenbasis)

h-

h { Diskretisierungsparameter

Abbildung 4
...
Dabei ist uh(x) eine stetige, aber keine stetig di erenzierbare Funktion
...


108

KAPITEL 4
...
ORDNUNG

Vorteile der Wahl von Ansatzfunktionen mit lokalem Trager

1
...

2
...
Folglich
ist die Stei gkeitsmatrix schwach besetzt
...
Die Integrale in der De nition der Stei gkeitsmatrix und des Lastvektors konnen
elementweise berechnet werden
...
Vh ! V (Vg1h ! Vg1 V0h ! V0 ) fur h ! 0 , d
...
die Funktionen aus V
lassen sich fur feiner werdende Vernetzungen immer besser durch Funktionen aus
Vh approximieren
...
7 Approximation krummliniger Rander
(a) Approximation des Randes durch einen Polygonzug, d
...
Vernetzung des Gebietes
mit geradlinig begrenzten Elementen
...
3: Approximation krummliniger Rander durch einen Polygonzug

109

4
...
GALERKIN{RITZ{FEM

(b) Verwendung krummlinig berandeter Dreieckselemente (siehe z
...
27, 91])
=)

=

h

Abbildung 4
...
8 Ansatzfunktionen hoherer Ordnung auf Dreieckselementen

Wir werden uns im weiteren nur mit stuckweise linearen Ansatzfunktionen beschaftigen
...
2 dargestellt
...
Die Abbildung 4
...
Weitere
Beispiele kann der Leser in 27] und 91] nden
(a) quadratisches Element

p(x) = a0 + a1x1 + a2x2 + a3x2 + a4x1x2 + a5x2
1
2

s

s

s

J
J
J
J

s

s

s

J
J
J

Abbildung 4
...
6: Ein Dreieckselement mit kubischen Ansatzfunktionen

110

KAPITEL 4
...
ORDNUNG

Bemerkung 4
...
B
...
Hau g genutzte Elemente sind z
...
das Viereckselement mit
bilinearen Ansatzfunktionen, das 8{Knoten Element mit Ansatzfunktionen der Gestalt
p(x) = a0 +a1x1 +a2x2+a3x1x2 +a4x2 +a5x2 +a6x2x2 +a7x1x2 (Serendipity{Element)
1
2
1
2
und isoparametrische Viereckselemente zur Approximation krummliniger Rander (siehe
auch 91])
...
Ordnung)

Abbildung 4
...
10 3D{Elemente

Zur FE{Diskretisierung von Randwertproblemen in 3D{Gebieten werden Tetraederelemente und Quaderelemente genutzt (siehe z
...
91])
...
Ordnung

Abbildung 4
...
4
...
4 FEM mit linearen Dreieckselementen
4
...
1 Gebietsdiskretisierung (Triangularisierung)

Wir zerlegen das Gebiet , in dem die Losung des betrachteten Randwertproblems
gesucht wird, in nite Elemente (r), z
...
in "kleine\ Dreiecke
...
Fur die
Triangularisierung Th soll gelten:
1
...


S
= r2R

h

(r)

bzw
...
a
...
),
{ den Eingangsdaten des Randwertproblems (z
...
unstetige Koe zienten ! Interfacelinien, rechte Seite, gemischte Randbedingungen ! Punkte, in denen der Typ
der Randbedingungen wechselt,
...
des Polynomgrades der Ansatzfunktionen)
Generelle Hinweise fur die Generierung einer Vernetzung

{ unzulassige Vernetzungen

s s
s s s
L

A

A

A
A

@
@

L
L
L

s

@
@
@

s s
s s s
s
L
A
L
A
L
L
A
A
A
@
@
@
@
@

Abbildung 4
...
FEM FUR MEHRDIMENSIONALE RANDWERTPROBLEME 2
...
Art

RB 1
...
Art

RB 1
...
10: Beachtung des Wechsels des Typs der Randbedingung
{ Vernetzung am Interface
falsch

richtig

s

Interface

@
@
@
(
(
((((
@
(((
@
@
!
!!
!!
!!

s

s

s

Interface

@
@
@ (((
(
(s
s
(((( @
J
@
@
J
!
!
J !!
Js
!
!
!!

s

s

Abbildung 4
...
a
...
B
...


Abbildung 4
...
4
...
13: Dreieck mit einem sehr stumpfen Innenwinkel
{ Generierung einer regularen Triangularisierung, d
...

fur alle r 2 Rh und fur alle Diskretisierungsschrittweiten h existieren Konstanten 0 und 0 mit 0 > 0 und
0 > 0, so da
0<

0

(r )
1

(r )
2

(r)
3

h(r)
2
P1(r)

h(1r) h(2r) h(3r) h

0h

;

P3(r)
(r)
3

(r)
1

0

h(r)
3

h(r)
1
(r)
2

P2(r)

Abbildung 4
...

global:
Numerierung aller Knoten und Elemente
!h = f1 2 : : : Nhg Nh die Anzahl der Knoten
Rh = f1 2 : : : Rhg Rh die Anzahl der Elemente
Festlegung der Knotenkoordinaten xi = (x1 i x2 i)
i = 1 2 : : : Nh
lokal:
Numerierung der Knoten in jedem Dreieck

sP

(r )
2

P3(r)

s

Knotenkoordinaten der Knoten P (r) :

(r)

s

P1(r)

x(r) = (x(1r) x(2r) )
2 A = A(r) = f1 2 3g

Abbildung 4
...
FEM FUR MEHRDIMENSIONALE RANDWERTPROBLEME 2
...
16: Beispiel einer Vernetzung mit Numerierung der Knoten und Elemente

115

4
...
FEM MIT LINEAREN DREIECKSELEMENTEN

Alle FE{Programme benotigen deshalb mindestens die folgenden zwei Felder, deren
Aufbau wir am Beispiel der in der Abbildung 4
...
1) erlautern
...


...

86
87

...


...
B
...


...


...


...


...


...


...


...


...


...


...


...


...
1: Zuordnungstabelle zwischen lokaler und globaler Knotennumerierung

i : (x1 i x2 i)
i
x1 i
x2 i

1
0
...
0

2
1
...
0

3
1
...
6

Nh = 72
0
...
7

Tabelle 4
...
11 Weitere Felder zur Charakterisierung von Knoten sind moglich,
z
...
zur Randbeschreibung und zur Randbedingungskodierung
...
Erzeugung der Felder fur die Zuordnung zwischen der lokalen und der globalen
Knotennumerierung (r :
$ i) und der Knotenkoordinaten (i : (x1 i x2 i))
von "Hand\
...
Halbautomatische Erzeugung der Daten, d
...
zum Beispiel Zerlegung des Gebietes
in Teilgebiete
...
Die Vernetzungen
der Standardgebiete werden auf die konkreten Teilgebiete abgebildet
...


116

KAPITEL 4
...
ORDNUNG

=S
e

e

e

{ Rechteck

e

{ krummliniges
Rechteck, d
...

'(Rechteck)

s
s
s
s
s
s

s
s
s
s
s
s

s
s
s
s
s
s

Y
HH
H

'

s s s
s s s
s s s

Abbildung 4
...
18: Abbildung der Vernetzung eines Standarddreiecks auf ein beliebiges
Dreieck
3
...
Als Eingangsdaten werden
i
...
entweder eine Beschreibung des Gebietsrandes @ (z
...
67]) oder eine Zerlegung des Gebietes in Teilgebiete e mit = Se e und die Beschreibung der
Teilgebietsrander @ e ( 15]) gefordert
...
B
...

Es werden folgende Schritte ausgefuhrt:
CAD{System ! gra sche Benutzer- ! Netzgenerator ! Vernetzung
ober ache
"

Gebietsbeschreibung

"

Daten fur den
Netzgenerator

#

Netzdaten le

117

4
...
FEM MIT LINEAREN DREIECKSELEMENTEN

Wir beschreiben diese Vorgehensweise am Beispiel der Generierung der Vernetzung
des Gebietes aus dem Abschnitt 4
...
Das Gebiet wird in zwei Basisteilgebiete
(BTG) zerlegt
...
Die Basislinien P1P2, P2P3, P3 P4, P4P7 , P7P8, P8P9, P9 P10 sowie P10P1 bilden
den Rand des Teilgebietes 1 und der Rand des Teilgebietes 2 wird in die Basislinien
P4P5, P5 P6, P6P7 sowie P7P4 unterteilt
...


x2
0:8

6

t b b t5
b BTG 2 b
b
7 t
b b bt4b b b t 3
bb
b
8 t
bb
b
9 t BTG 1
b
bb
10 t
bb
b
1t
b b b t2
6

-

x1
0
1
Abbildung 4
...
20 dargestellte FE{Vernetzung automatisch erzeugt
...
20: Automatisch generierte FE{Vernetzung

118

KAPITEL 4
...
ORDNUNG

4
...
a posteriori Informationen
H
9

HH
H
j
H

vor der FE{Rechnung

nach einer FE{Rechnung

Analyse der Eingangsdaten
(Existenz uberstumpfer Ecken am Gebietsrand, Sprunge in den Koe zienten
der Di erentialgleichung u
...
)
;! Festlegung einer Netzgraduierung
(siehe z
...
2, 63])

Analyse der FE{Losung,
Schatzung des Fehlers
;! Festlegung der Dreiecke, in denen
der Fehler gro ist und Forderung einer weiteren Verfeinerung dieser Dreiecke (siehe z
...
5, 7, 83, 94, 97])

Die Netzverfeinerung erfolgt durch
Dreiecksviertelung
;!

Abbildung 4
...
22: Halbierung eines Dreiecks

4
...
2 De nition der Ansatz- und Testfunktionen
Prinzip

Lokale De nition der Basisfunktionen (Ansatz-, Testfunktionen) uber die Formfunktionen, die ihrerseits durch Abbildung der Formfunktionen des Referenzdreiecks (Masterelement, Basiselement, Einheitselement) erhalten werden
...
4
...
auch die Abbildung 4
...
Es gilt
jdet J (r) j = 2 meas

(r )

wobei meas (r) den Flacheninhalt des Elements (r) bezeichnet
...

x2

s

2 6

x(r) = xk
3

6

x(r) = xi
1

s

(r)

x(r)
2

s= x

beliebiges Dreieck

1

s

0

= (x)-

s

=3

x = x( )

j

-

x1

(r) 2 T
h

=1

s

1
=2

-

1

Referenzdreieck

Abbildung 4
...
h
...
13)

mit Bi = fr 2 Rh : xi 2 (r)g und p(r)(x) = ' ( (r) (x)) x 2 (r)
...


120

KAPITEL 4
...
ORDNUNG

Fur die im Abschnitt 4
...
1 vorgestellte Vernetzung gilt beispielsweise B58 = f57 58
59 66 68 69g
...
24: De nition der Ansatzfunktion p58(x)
Mit

pi(xj ) =
Vh
V0h
Vg1h

fur alle i j 2 ! h ;! Knotenbasis
n
o
X
= vh(x) : vh(x) = vipi (x) V = H 1( )
i2!h
n
o
X
= vh(x) : vh(x) = vipi (x) V0
i2!h
n
o
X
X
= vh(x) : vh(x) = vipi (x) + g1(xi)pi (x)
ij

i2!h

i2

h

folgt fur die Koe zienten vi in der Darstellung der Funktionen vh(x)

vh(xi) = vi

fur alle i 2 !h:

In dem Beispiel aus dem Abschnitt 4
...
1 haben wir

!h = f1 2 : : : 72g

und

!h = ! h n

h

h

= f8 9 10 25 26 27 28g:

1

2

121

4
...
FEM MIT LINEAREN DREIECKSELEMENTEN

4
...
3 Aufbau des FE{Gleichungssystems

Wir beschreiben im weiteren den elementweisen Aufbau des Finite{Elemente{Gleichungssystems, d
...
des Ritz{Galerkin{Systems (siehe auch Abschnitt 4
...
14)

Kh = a(pi pk )]i k2!h
X
f h = hF pk i ; u ia(pi pk )]k2!h
i2

h

gilt, wobei die Bilinearform a(: :) und die Linearform hF :i wie folgt de niert sind :
Z @pi @pk
Z
@p @p
a(pi pk ) =
+ 2 @xi @xk dx +
(x)pi(x)pk (x) ds
1
@x1 @x1
2 2
;3
Z
Z
Z
hF pk i = f (x)pk (x) dx + g2 (x)pk (x) ds +
(x)g3(x)pk (x) ds:
;2

;3

Zuerst generieren wir eine Stei gkeitsmatrix und einen Lastvektor, in denen die Randterme noch nicht berucksichtigt sind
...
Den Ausgangspunkt unserer weiteren Uberlegungen bilden die Beziehungen
(Kh uh vh ) = a(uh vh) fur alle uh $ uh vh $ vh uh vh 2 Vh uh vh 2 RNh
und
(f h vh) = hF vhi fur alle vh $ vh vh 2 Vh vh 2 RNh :

122

KAPITEL 4
...
ORDNUNG

Berechnung der Elementstei gkeitsmatrizen

Es gilt fur beliebige uh vh 2 Vh mit uh =

a(uh vh) = a
=

Z

X
i2!h
1

uipi(x)

X

i2!h

X

k2!h

uipi(x) und vh =

X
k2!h

vk pk (x) :

vk pk (x)

@ X u p (x) @ X v p (x)
@x1 i2!h i i @x1 k2!h k k

@ X u p (x) @ X v p (x) dx
+ 2 @x
i i
@x2 k2!h k k
2 i2!h
=

X X X
r2Rh i2!h k2!h

uivk

Z

1

(r )

@pi @pk + @pi @pk dx:
@x1 @x1 2 @x2 @x2

Fuhren wir die Indexmenge ! (hr) = fj : xj = (x1 j x2 j ) 2 (r)g ein, und beachten wir,
da
Z
@pi @pk + @pi @pk dx = 0
1
@x1 @x1 2 @x2 @x2
(r )
gilt, falls xi 2
=

oder xk 2 (r) , so folgt
=
Z
X X X
a(uh vh) =
uivk
(r)

r2Rh i2!(r) k2!(r)
h
h

(r )

1

@pi @pk + @pi @pk dx:
@x1 @x1 2 @x2 @x2

Mit der Zuordnungsvorschrift r : i $ , k $ ,
2 A(r) = f1 2 3g , erhalten
wir (siehe auch die Beziehung (4
...


4
...
FEM MIT LINEAREN DREIECKSELEMENTEN

123

Wir haben die Transformationsvorschrift x (r) = x (r) ( ) :
!
!
!
!
!
!
x1 = x1 j ; x1 i x1 k ; x1 i
1 + x1 i
1 + x1 i
x2
x2 j ; x2 i x2 k ; x2 i
x2 i = J (r) 2
x2 i
2
wobei Pi (x1 i x2 i) Pj (x1 j x2 j ) Pk (x1 k x2 k ) die Eckknoten des Dreiecks (r) sind,
und die Zuordnung i $ 1 j $ 2 k $ 3 zwischen der globalen und der lokalen
Knotennumerierung besteht
...
h
...
13)
K (r ) =

Z
(r )

=

Z

@p(r) @p(r) + @p(r) @p(r) dx
1 @x @x
2 @x @x
1
1
2
2
@p(r)(x (r) ( )) @ 1 + @p(r)(x (r) ( )) @ 2
1(x (r) ( ))
@1
@x1
@2
@x1
@p(r)(x (r) ( )) @ 1 @p(r)(x (r) ( )) @ 2
@1
@x1 +
@2
@x1

+

(r) x
(r) x
@
@
(x (r) ( )) @p (@ (r) ( )) @x1 + @p (@ (r) ( )) @x2
2
1
2
2
2

@p(r)(x (r) ( )) @ 1 @p(r)(x (r) ( )) @ 2
@1
@x2 +
@2
@x2 jdetJ (r) j d

124
=

KAPITEL 4
...
ORDNUNG

Z

1 (x

(r )

( )) @' ( ) @ 1 + @' ( ) @

@

@x1

1

@

2
@x1
2

@
@
+ 2(x (r) ( )) @' ( ) @x1 + @' ( ) @x2
@1 2
@2 2

Bei linearen Formfunktionen gilt

'1( ) = 1 ; 1 ;
'2( ) =

2

so da
(r
K11)

1 Z
= jdetJ (r) j

1 (x

(r )

@' ( ) @ 1 + @' ( ) @ 2 jdetJ j d
(r )
@ 1 @x2
@ 2 @x2

@'1 = ;1
@1
@'2 = 1
@1
@'3 = 0
@1

2

1

'3( ) =

@' ( ) @ 1 + @' ( ) @ 2
@ 1 @x1
@ 2 @x1

@'1 = ;1
@2
@'2 = 0
@2
@'3 = 1
@2

( ))(;(x2 k ; x2 i) + (x2 j ; x2 i))2 +

( ))((x1 k ; x1 i) ; (x1 j ; x1 i))2] d
Z
Z
= jdet1 j (x2 j ; x2 k)2 1(x (r) ( )) d + (x1 k ; x1 j )2 2(x (r) ( )) d
J (r)
2 (x

(r
K22)

1 Z
= jdetJ (r) j

1 (x

(r )

( ))(x2 k ; x2 i)2 + 2(x (r) ( ))(;(x1 k ; x1 i))2] d

Z
1
2
= jdetJ j (x2 k ; x2 i)
(r )
(r
K33)

1 Z
= jdetJ j
(r )

1 (x

(r )

1 (x

(r )

1 (x

(r )

( )) d + (x1 i ; x1 k

)2

Z

2 (x

(r )

( )) d

( ))(;(x2 j ; x2 i))2 + 2(x (r) ( ))(x1 j ; x1 i)2] d

Z
1
2
= jdetJ j (x2 i ; x2 j )
(r )
Z
(r
K12) = jdet1 j
J (r)

(r )

1 (x

(r )

( )) d + (x1 j ; x1 i

)2

Z

2 (x

(r )

( )) d

( ))(x2 k ; x2 i)(;(x2 k ; x2 i) + (x2 j ; x2 i)) +

2 (x

(r )

( ))(;(x1 k ; x1 i)((x1 k ; x1 i) ; (x1 j ; x1 i)))] d

125

4
...
FEM MIT LINEAREN DREIECKSELEMENTEN

Z
= jdet1 j (x2 k ; x2 i)(x2 j ; x2 k )
J (r)

1 (x

(r )

( )) d +

(x1 i ; x1 k)(x1 k ; x1 j )
(r
K13)

1 Z
= jdetJ j
(r )

1(x

(r )

1
jdetJ (r) j

2 (x

(r )

( )) d

( ))(;(x2 j ; x2 i)(;(x2 k ; x2 i) + (x2 j ; x2 i))) +

2 (x

=

Z

( ))(x1 j ; x1 i)((x1 k ; x1 i) ; (x1 j ; x1 i))] d
Z
(x2 i ; x2 j )(x2 j ; x2 k ) 1(x (r) ( )) d +
(r )

(x1 j ; x1 i)(x1 k ; x1 j )

Z
(r
K23) = jdet1 j
J (r)

1(x

(r )

Z

2 (x

(r )

( )) d

( ))(;(x2 j ; x2 i)(x2 k ; x2 i)) +

2 (x

( ))(x1 j ; x1 i)(;(x1 k ; x1 i))] d
Z
1
= jdetJ (r) j (x2 i ; x2 j )(x2 k ; x2 i) 1(x (r) ( )) d +
Z
(x1 j ; x1 i)(x1 i ; x1 k ) 2(x (r) ( )) d
(r )

R
R
Es mussen noch die Integrale 1 d und 2 d berechnet werden
...

2 konstant sind, ist dies trivial
...
h
...
FEM FUR MEHRDIMENSIONALE RANDWERTPROBLEME 2
...
k1 = 1 1(x (r) ( s )) und k2 = 1
2
2
Elementstei gkeitsmatrix die Gestalt

2

bzw
...

Bemerkung 4
...
Eine explizite Darstellung der Formfunktionen uber dem Element (r) wurde nicht benotigt
...
h
...

Berechnung der Elementlastvektoren
Die Elementlastvektoren f (r) werden auf analoge Weise wie die Elementstei gkeitsmatrixen K (r) berechnet, d
...

Z
Z
X
(f h vh) = hF vhi = f (x)vh(x) dx = f (x) vk pk (x) dx
=

X X
r2Rh k2!h

vk

k2!h

Z

f (x)pk (x) dx :

(r )

Fuhren wir wiederum die Indexmenge !(r) = fk : xk = (x1 k x2 k ) 2
h
beachten wir, da pk (x)j (r) = 0 , falls xk 2 (r) , dann erhalten wir
=
X X Z
(f h v h ) =
vk f (x)pk (x) dx :
r2Rh k2!(r)
h

(r)g

(r )

Mit der Zuordnungsvorschrift r : k $ folgt
X X (r) Z
X ( ( (
(f h vh) =
v
f (x)p(r)(x) dx = (v1r) v2r) v3r))f (r)
r2Rh

2A(r)

r2Rh

(r )

wobei der Elementlastvektor f (r) durch
hZ
i3
f (r ) =
f (x)p(r)(x) dx =1
(r )

ein, und

4
...
FEM MIT LINEAREN DREIECKSELEMENTEN

127

de niert ist
...
Es gilt

Z

f (x)p(r)(x) dx

=

(r )

=

Z
Z

f (x (r) ( ))p(r) (x (r) ( ))jdetJ (r) j d
f (x (r) ( ))' ( )jdetJ (r) j d :

Zur Berechnung der Integrale kann beispielsweise die Mittelpunktsformel eingesetzt
werden, d
...

Z
1
f (x (r) ( ))' ( )jdetJ (r) j d 2 f (x (r) ( s))' ( s)jdetJ (r) j mit s = 1 1 :
3 3

Bemerkung 4
...
h
...


Assemblierung der Elementstei gkeitsmatrizen und der Elementlastvektoren
Zur Assemblierung der Elementstei gkeitsmatrizen bzw
...
h
...
Wir schreiben diese fur
unser Beispiel aus dem Abschnitt 4
...
1 auf
...


...

10
11

...


...
B
...


...


...


...


...


...


...


...


...


...


...


...


...
3: Zuordnungstabelle
Es seien die Elementstei gkeitsmatrizen K (r), r = 1 2 : : : Rh , berechnet
...
Nach dem Einbau von K (1) erhalten wir

128

KAPITEL 4
...
ORDNUNG

1

2

29 30 31

36 37 38

72

0

0
0

...


...


...

0
0
0

...


...


...

0 0
(1)
0 K33
0 0

...


...


...

0 0

(1)
0 K12
0 0

...


...


...

0 0
(1)
0 K22
0 0

...


...


...

0
0
0

...


...


...

0

0 (1)
B K11
B 0
B
B
...

B
...

B
...

B
...

@

1
2

...


29
30
31

...

36
37
38

...

72


...



...


...

0
0
0

...


...


...

0


...



...


...

0
0
0

...


...


...

0


...



...


...

0
0
0

...


...


...

0
0
0

...


...


...

0 0
(2)
0 K22
0 0

...


...


...

0 0
(2)
0 K32
0 0

...


...


...

0 0
0 0
0 0

...


...


...

0 0
(1)
0 K23
0 0

...


...


...


...


...


...


...


...


...


...


...


...


...


...


...


...
B
B
B 0
10 B
B
B
(2)
11 B K21
B
B 0
12 B
B
B

...
B

...


...


...


...
B
B
B
B 0
36 B
B (1) (2)
B
37 BK21 + K31
B
B 0
38 B
B
B

...


...
B

...



...



...


...

0
0
0

...


...


...

0
0
0

...


...



...



...


...

0
0
0

...


...


...

0
0
0

...


...



...



...


...

0
0
0

...


...


...

0
0
0

...


...



...



...
C

...
C
C
0C
C
C
C
0C
C
C
0C
C

...
C

...
C

...
C
C
0C
C
C
C
0C
C
C
0C
C

...
C

...
4
...

Fur die Assemblierung der rechten Seite f h nutzen wir den gleichen Algorithmus,
d
...
zunachst setzen wir f h identisch Null
...
Fur das Element f72 gilt z
...


f72 = f1(100) + f1(101) + f1(102) + f1(103) + f1(104) + f1(105) + f1(106) :

Bemerkung 4
...
Dabei ist C (r) wie folgt de niert
8
> 1 i globale Knotennummer eines Eckknotens des Dreiecks (r)
<
(
und j die zugehorige lokale Knotennummer
Cijr) = >
: 0 sonst

i = 1 2 : : : Nh, j = 1 2 3
...
FEM FUR MEHRDIMENSIONALE RANDWERTPROBLEME 2
...
und 3
...
Art nicht berucksichtigt
...

Berucksichtigung der Randbedingungen 2
...
Art im FE{Gleichungssystem

Wir setzen voraus, da das vernetzte Gebiet polygonal berandet ist
...
4
...
;(e3 ) Dreiecksseiten sind, die auf ;2 bzw
...
Analog wie bei
2
3
den Dreiecken numerieren wir alle Dreiecksseiten, die(3) ;2 bzw
...

Au erdem de nieren wir eine Zuordnungsvorschrift e : i $ , 2 A(e) = f1 2g ,
zwischen der globalen und der lokalen Knotennumerierung, d
...
in unserem Beispiel
Nummer der Dreiecksseite,
die auf ;3 liegt
1
2
3
4

globale Knotennummern der lokalen Knoten
P1(e3 )
P2(e3 )
1
11
11
12
12
13
13
2

Tabelle 4
...
;3 werden wie folgt berechnet:

Z
;3

(x)uh(x)vh(x) ds =
=

Z

X
e3 2Eh ;
(3)

(x)uh (x)vh(x) ds

(e3 )
3

X X X
e3 2Eh i2!h k2!h
(3)

uivk

Z

(x)pi (x)pk (x) ds :

;(e3 )
3

Fuhren wir die Indexmenge ! (e3) = fj : xj = (x1 j x2 j ) 2 ;(e3)g ein, und beachten
h
3
(e3 )
wir, da pj (x)j (e3) = 0 fur xj 2 ;3 gilt, dann folgt
=
;3

131

4
...
FEM MIT LINEAREN DREIECKSELEMENTEN

Z
;3

(x)uh(x)vh(x) ds =

X X X
e3 2Eh i2!h k2!h
(e3 )

(3)

(e3 )

uivk

Z

(x)pi(x)pk (x) ds

;(e3 )
3

und mit der Zuordnungsvorschrift e : i $ , e : k $

Z

X

(x)uh (x)vh(x) ds =

;3

X

(3)
e3 2Eh

X

=

(3)
e3 2Eh

mit

K (e3) =

(e )
(e )
2A 3 2A 3

Z

u(e3)v(e3)

(
(
(v1e3) v2e3))K (e3)

"Z
;

X

(x)p(e3 )(x)p(e3)(x) ds

;(e3 )
3

u(1e3)
u(2e3)

(x)p(e3 )(x)p(e3 )(x) ds

!

#2
=1

(e3 )
3

:

Die Elemente der Matrizen K (e3 ) berechnen wir, indem wir in den Kurvenintegralen eine
Koordinatentransformation durchfuhren
...
Es gilt
!
!
!
x1 = x1 k ; x1 i
x1 i :
x = x( 1) :
x2
x2 k ; x2 i 1 + x2 i
x2

s

2 6

(
x2e3 ) = xk

6

s

s

1

;(e3 )
3

x(e3 ) = xi
1

= 1(x)

x = x( 1)
-

x1

beliebige Dreiecksseite ;(e3 )
3

0

s
=1

s

1
=2

-

Referenzintervall 0 1]

Abbildung 4
...
FEM FUR MEHRDIMENSIONALE RANDWERTPROBLEME 2
...
h
...
4
...

Abschlie end sind noch die Randbedingungen 1
...


Berucksichtigung der Randbedingungen 1
...
Variante : Homogenisieren im Diskreten
Wir setzen uj = g1(xj ) fur alle j 2 h = !h n !h , d
...
fur alle Indizes j , die zu Knoten
auf dem Rand ;1 gehoren, und korrigieren die rechte Seite durch
X
fi := fi ; Kij uj fur alle i 2 !h
j2

h

wobei Kij die Matrixelemente der globalen Stei gkeitsmatrix Kh nach Einbau der Elementstei gkeitsmatrizen K (r) und der Matrizen K (e3) sowie fi die Komponenten des
globalen Lastvektors nach Einbau der Elementlastvektoren f (r) und der Vektoren f (e2)
bzw
...
Anschlie end streichen wir die Spalten j , j 2 h, und die Zeilen i,
i 2 h, aus der Stei gkeitsmatrix sowie die Zeilen i, i 2 h , aus dem Lastvektor
...
Variante:
Wir gehen analog wie in der Variante 1 vor, nur anstelle des Streichens der Spalten j
und der Zeilen i mit i j 2 h setzen wir
(
Kij = ij = 1 fur i = j
sowie fi = g1(xi)
0 fur i 6= j
fur alle i j 2 h
...
Variante: Straftechnik (siehe auch 44])
Fur alle i 2 h setzen wir
^
^
Kii := Kii + K K = 10p max jKij j Nh p hinreichend gro
j =1 2 ::: Nh

^
fi := fi + Kg1(xi) :

Zusammenfassung zum Aufbau des FE{Gleichungssystems
1
...
h
...


h

=

(r)

r2Rh

134

KAPITEL 4
...
ORDNUNG

2
...
den globalen Lastvektor
(
3
...
Art)

{ Berechnung der Matrizen K (e3)
"Z
#2
(e3 ) =
(e3 ) (x)p(e3 ) (x) ds
K
(x)p
;(e3 )
3

sowie der Vektoren

f (e3 ) =

"Z
;

(e3 )

(x)g3(x)p (x) ds

=1

#2

(e3 )
3

=1

{ Einbau der Matrizen K (e3 ) und der Vektoren f (e3) in die globale Stei gkeitsmatrix und den globalen Lastvektor
(
4
...
Art)

{ Berechnung der Vektoren

f (e2)

=

"Z
(e2 )
2

;

g2

#2

(x)p(e2 )(x) ds

=1

{ Einbau der Vektoren f (e2 ) in den globalen Lastvektor
5
...
Art im FE{Gleichungssystem

Bemerkung 4
...
dem Referenzintervall zuruckgefuhrt
...
4
...
4
...
Wir setzen voraus, da das Gebiet polygonal berandet ist, so
da
=
gilt
...
> 0

es existiert ein
es existiert ein

(r )

:
0 = const
...
4
...
14)
...
h
...
> 0 : ja(u v)j
1 = const
...
15)

fur alle u v 2 V0
...
h
...
B
...


136

KAPITEL 4
...
ORDNUNG

4
...
5 Ein Beispiel

Wir betrachten folgendes Randwertproblem:
Gesucht ist das Temperaturfeld u(x) 2 C 2( ) \ C 1(
; u = 0

x 2 = (0 1) (0 1)
x 2 ;11 = fx 2 @ : x1 = 1 0 < x2 1g
x 2 ;12 = fx 2 @ : 0 x1 1 x2 = 0g

u(x1 x2) = x2
u(x1 x2) = 0

gilt
...
26: Darstellung des Gebietes
Die exakte Losung dieser Aufgabe ist u(x1 x2) = x1x2 , denn

@u = x
@x1 2

@ 2u = 0
@x2
1

@u = x
@x2 1

@ 2u = 0
@x2
2

@u = @u n + @u n = ; @u = ;x
2
@n @x1 1 @x2 2
@x1

auf ;21

@u = @u n + @u n = @u = x
@n @x1 1 @x2 2 @x2 1

auf ;22 :

(4
...
4
...
16)
...
16) mit einer Testfunktion v(x) 2 V0 , V0 = fv(x) 2
H 1( ) : v(x) = 0 auf ;11 ;12g und integrieren uber , so da wir
Z
Z
;
u(x) v(x) dx = 0 v(x) dx
erhalten
...
16) lautet somit:
Gesucht ist u(x) 2 Vg1 = fu(x) 2 H 1( ) : u(x) = x2 fur alle x 2 ;11 und u(x) = 0
fur alle x 2 ;12g , so da
Z @u @v @u @v
Z
Z
+ @x @x dx = ; x2v ds + x1v ds
(4
...

Im weiteren fuhren wir die FE{Diskretisierung durch
...


3r

1

2r

2;
;3
1

r
1

2;
;3

6r

4

;
;
;
;

2

3

;
;
;
;

1

r9

3; 1
;2

8

2
51 r;;
3

3; 1
;2

r
4

;

2
1 ;;
3

;
;

;
;

6

3;
;2

7

r8

1
3;
;2

;
;
;

5

1

r7

Abbildung 4
...
FEM FUR MEHRDIMENSIONALE RANDWERTPROBLEME 2
...
4
...

Zunachst generieren wir das FE{Gleichungssystem ohne Berucksichtigung der Randbedingungen
...

Elementnummer
1
2
3
4
5
6
7
8

globale Knotennummern der lokalen Knoten
P1(r)
P2(r)
P3(r)
4
5
1
2
1
5
5
6
2
3
2
6
7
8
4
5
4
8
8
9
5
6
5
9
Tabelle 4
...
4
...
In unserem Beispiel
gilt 1 = 2 = 1
...


1
C
C
C
C
C
A

4
...
FEM MIT LINEAREN DREIECKSELEMENTEN

139

Damit erhalten wir die globale Stei gkeitsmatrix Kh
0
1
B 2+2
-1
-2
0
0
0
B1 1
2
B
B
B
B
B -1 1 + 1 + 1 -1
B 2
-1 - 1
0
0
B
2 2
2
2 2
B
B
B
B
B 0
1
1
-2
-2
B
1
0
0
B
B
B
B
B 1
1
B -2
1+ 1 + 2 -1 - 1
0
0
0
B
2
2 2
B
B
B
B
B
1 + 1 +1+
B
B 0
1
1
-1 - 1
- 2 - 1 2 21 1 - 2 - 1
B
0
2 2
2 1+ +
2
B
B
2 2
B
B
B
B
B
1
B 0
-1
-1 - 1 2 + 1 +1
0
0
2
2 2
2
B
B
B
B
B
B 0
1
-2
B
0
0
0
0
B
B
B
B
B
1
B 0
-2 - 1
0
0
0
0
B
2
B
B
B
B
B
1
@ 0
-2
0
0
0
0

1
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
A

d
...


0

0

0

0

0

0

0

0

0

1
-2

0

0

0

1
-1 - 2
2

0

0

0

-1
2

1

-1
2

0

1
-2

1 + 1 +1
2 2

-1
2

0

-1
2

0 1 -1 0 -1 0 0 0 0 0 1
B 1 2 1 2
C
B - 2 2 - 2 0 ;1 0 0 0 0 C
C
B
B
C
B 0 -1 1 0 0 -1 0 0 0 C
C
B
2
C
B 1 2
B1
C
B 2 0 0 2 ;1 0 - 2 0 0 C
C
B
C
B
Kh = B 0 ;1 0 ;1 4 ;1 0 ;1 0 C :
C
B
B
C
B 0 0 - 1 0 ;1 2 0 0 - 1 C
B
2
2 C
B
C
B 0 0 0 -1 0 0 1 -1 0 C
C
B
2
2
B
C
B 0 0 0 0 ;1 0 - 1 2 - 1 C
B
2
2 C
A
@
1 0 -1 1
0 0 0 0 0 -2
2

1
2

+

1
2

140

KAPITEL 4
...
ORDNUNG

Der globale Lastvektor f h ist identisch Null, da die rechte Seite der Di erentialgleichung
identisch Null ist
...
Art
berechnet, d
...

Z
Z
; x2vh ds und
x1vh ds :
;21

;22

Der Rand ;21 wird durch die zwei Dreiecksseiten ;(1) = P1P2 und ;(2) = P2P3 sowie
2
2
der Rand ;22 durch die Dreiecksseiten ;(3) = P6P3 und ;(4) = P9 P6 gebildet
...

Nummer der Dreiecksseite,
die auf ;2 liegt
1
2
3
4

globale Knotennummern der lokalen Knoten
P1(e2)
P2(e2 )
1
2
2
3
6
3
9
6

Tabelle 4
...
Es gilt mit der Abbildungsvorschrift
!
!
!
x1 = x1 k ; x1 i
x1 i
x = x( 1) :
x1
x2 k ; x2 i 1 + x2 i

141

4
...
FEM MIT LINEAREN DREIECKSELEMENTEN

Z

;x2 p(1)(x) ds
1

=

;

(1)
2

=

Z
;

;x2 p(1)(x) ds
2

=

(1)
2

=

Z1
0

Z1
0

Z1
0

Z1
0

;(x2 k ; x2 i) 1 ; x2 i](1 ; 1 )

;

1 (1 ; ) 1 d
1
2 1
2

1

1
2

2
1

1 d
2

1

1

1
= ; 24

;(x2 k ; x2 i) 1 ; x2 i]

;

1 d
2

1

1 d
2

1

1
= ; 12

und somit
1
1
f (1) = ; 24 ; 12
Auf analoge Weise erhalten wir
5 T f (3) = 1
f (2) = ; 1 ; 24
6
12

1
24

T

T

:
5
und f (4) = 24

1
6

T

:

Damit ergibt sich der Lastvektor
1 ; 1 ;1 ; 5 + 1 0 0 1 +1 0 0 5
24
12 6
24 24
12 6
24
1 ; 1 0 0 1 0 0 5 T:
; 1 ;
24
4
6
4
24
;

T

=

Die Randbedingungen 1
...
4
...
So erhalten wir schlie lich das Gleichungssystem
0
1 0 1 0 -1 1
1
2 - 2 ;1 0
B 1
C B u2 C B 4 C
B - 1 0 - 1 C B u3 C B - 1 C
B 2
C B 6 C
2 CB
B
B ;1 0 4 ;1 C B u5 C = B 0 + 1 C :
CB C B
B
C
CB C B 2 C
@
A
A@ A @
1+1
1 ;1
2
0 -2
u6
4 2

142

KAPITEL 4
...
ORDNUNG

In der folgenden Tabelle vergleichen wir die analytische Losung der Aufgabe (4
...


u1
u2
u3
u4
u5
u6
u7
u8
u9

analytische Losung FE{Losung
0
0
1
0
24
1
0
8
0
0
1
4
1
2

13
48
13
24

0

0:0

1
2

1
2

1

1

Tabelle 4
...
Insbesondere ist der Fehler in
der Komponente u3 noch sehr gro
...
h
...


4
...
6 Das Programm FEM2D

Mittels des Programms FEM2D konnen Randwertprobleme 2
...
Es ist moglich, dieses Programm als
Demonstrationsprogramm zur Veranschaulichung des Prozesses der Finite{Elemente{
Diskretisierung einzusetzen bzw
...

Eine ausfuhrliche Dokumentation des Programms kann der Leser in 65] nden
...
h
...
Im Verlaufe des
Demonstrationsprozesses wird die Vernetzung des Gebietes angezeigt, und es konnen
die Ansatzfunktionen uber den Dreiecken dargestellt werden
...

Als Rechenprogramm kann das Programm FEM2D zur Losung von physikalischen Problemen eingesetzt werden, welche durch die folgende Randwertaufgabe modelliert werden
...
4
...

Dabei wird vorausgesetzt, da die Materialkoe zienten kxx, kyy , k0 sowie die Funktionen g2(x y) und g3(x y) stuckweise konstant sind
...
Die maximal
mogliche Anzahl von Materialbereichen ist 10
...
Bei der Beschreibung des Randwertproblems haben wir
hier nur die Bezeichnungen verwendet, die auch im Demonstrationsprogramm genutzt
werden
...


:
;3

@
@

;2

;
;

Abbildung 4
...
FEM FUR MEHRDIMENSIONALE RANDWERTPROBLEME 2
...
Dieses File
hat die folgende Struktur:
NN { Anzahl der Knoten, NE { Anzahl der Dreiecke
Feld der Knotenkoordinaten, d
...
Eingabe der Koordinaten (xi yi), fur alle i =
1 2 : : : NN (siehe auch Abschnitt 4
...
1)
Feld der Zuordnung zwischen der lokalen und der globalen Knotennumerierung,
d
...
Eingabe der globalen Nummern der drei Eckknoten jedes Dreiecks sowie einer
Materialbereichsnummer
...
Eckknotens, Nummer des 2
...
Eckknotens, Materialbereichsnummer
fur jedes Dreieck erfolgen
...

{ NB { Anzahl der geschlossenen Rander (im Beispiel: NB = 2)
{ Anzahl der Knoten auf jedem der NB Teilrander (im Beispiel: Eingabe von
16 und 56, denn 16 Knoten liegen auf dem 1
...
geschlossenen Rand (;3))
{ Eingabe der Knotennummern der Knoten, die auf dem jeweiligen geschlossenen
Randstuck liegen, d
...

Eingabe der Knoten, die auf dem ersten geschlossenen Randstuck liegen,
Eingabe der Knoten, die auf dem zweiten geschlossenen Randstuck liegen,

...


...

Dabei ist zu beachten, da die Eingabe der Knotennummern auf dem jeweiligen
Randstuck in fortlaufender Reihenfolge entgegen dem Uhrzeigersinn erfolgt
...
4
...
dat abgespeichert werden
(J/N) ?
Falls (N), Eingabe eines Filenamens
...
Art mittels Straftechnik { siehe Abschnitt 4
...
3)
{ Knotennummern, vorgegebener Wert

146

KAPITEL 4
...
ORDNUNG

Bei Vorhandensein von Randbedingungen 2: Art
{ Wieviele Randkanten mit Randbedingungen 2: Art ?
{ Fur jede Randkante Eingabe von Anfangsknoten, Endknoten und g2 aus der
Randbedingung
@u
kxx @x nx + kyy @u ny = g2
@y

Bei Vorhandensein von Randbedingungen 3: Art
{ Wieviele Randkanten mit Randbedingungen 3: Art ?
{ Fur jede Randkante Eingabe von Anfangsknoten, Endknoten, k0 und g3 aus
der Randbedingung
@u
kxx @x nx + kyy @u ny + k0u = k0g3
@y
Ausgabe der Ergebnisse uber Drucker (J/N) ?
Speichern der Ergebnisse auf Diskette (J/N) ?
Falls (J), Eingabe eines Filenamens
Ausgabe der Eingangsdaten (J/N) ?
Falls (J), Drucken der Eingangsdaten (J/N) ?

Bemerkung 4
...
3: Art wird nach einer Option fur die Eingabe gefragt
...

Sind die Randbedingungsdaten uber mehrere Randkanten gleich, dann ist die Eingabeform "B\ empfehlenswert, da nur der Anfangsknoten der ersten Randkante und der
Endknoten der letzten Randkante sowie die entsprechenden Daten eingegeben werden
mussen
...
1)

mit u 2 RN , f 2 RN und einer (N N ){Matrix erhalten
...

Zur Losung der FE{Gleichungssysteme wollen wir sowohl direkte als auch iterative
Verfahren einsetzen
...
Iterative Verfahren bestimmen den Vektor u
ausgehend von einer Startnaherung u(0) als Grenzwert einer Folge von Naherungslosungen u(k)
...
Derartige Gleichungssysteme entstehen bei der in den Kapiteln 3 und 4 beschriebenen FE{Diskretisierung von
Randwertproblemen mit symmetrischer und positiver Bilinearform
...
Die Matrix K ist positiv de nit, wenn
(Ku u) > 0

fur alle u 2 RN u 6= 0

erfullt ist
...
h
...
AUFLOSUNG VON FINITE{ELEMENTE{GLEICHUNGSSYSTEMEN

Die FE{Gleichungssysteme besitzen eine Reihe spezieller Eigenschaften
...
Ihre Dimension N wachst mit der Verfeinerung
des Diskretisierungsparameters (Schrittweite) h
...
Das bedeutet bei 2D{Randwertproblemen
ein Anwachsen wie h;2
...
3
...
h
...
Bei einer geeigneten Knotennumerierung besitzen die FE{Matrizen eine Bandstruktur
...
h
...
Dies
werden wir im Abschnitt 5
...
Man kann erreichen, da die Bandweite
in der Gro enordnung von h;(m;1) liegt, d
...
da sie kleiner oder gleich ch;(m;1) mit
einer von h unabhangigen Konstanten c ist
...
1 Direkte Verfahren
Als direktes Au osungsverfahren ist der Gau {Algorithmus allgemein bekannt (siehe
z
...
27])
...
Beim Cholesky{Verfahren wird die Matrix
K in das Produkt einer unteren und einer oberen Dreiecksmatrix zerlegt
...
Seien
0k k
k1N 1
11
12
B k21 k22
C
B
B
...
N C
C
K =B

...

A
kN 1 kN 2
kNN

0
B
B
S =B
B
@

s11
0

...


...


...


...

0
sNN

und

0r
11
B 0
B
R = B
...

@
0

r12
r1N
r22
r2N

...


...


...


1
C
C
C
C
A

149

5
...
DIREKTE VERFAHREN

Algorithmus I :

K = ST S

q
s11 = k11
s1j = k1j =s11

fur j = 2 3 : : : N

fur i = 2 3 : : : N
v
u
i;1
u
X
sii = tkii ; s2
li
l=1

sij = kij ;
Algorithmus II :

i;1
X
l=1

slislj

sii

fur j = i + 1 i + 2 : : : N

K = RRT

q
rNN = kNN
riN = kiN =rNN

fur i = N ; 1 N ; 2 : : : 1

fur j = N ; 1 N ; 2 : : : 1
v
u
N
u
X 2
t
rjl
rjj = ukjj ;
l=j +1

rij = kij ;

N
X
l=j +1

rilrjl

rjj

fur i = j ; 1 j ; 2 : : : 1

Die Losung des Gleichungssystems Ku = f ist somit der Losung des Gleichungssystems

S T Su = f

bzw
...
2)

aquivalent
...
2) losen wir in zwei Etappen durch die Prozesse
des Vor{ und Ruckwartseinsetzens, d
...
wir losen zunachst

ST y = f

bzw
...
3)

Su = y

bzw
...
4)

und dann das System
Die Gleichungssysteme (5
...
4) sind Gleichungssysteme, deren Systemmatrix
eine Dreiecksgestalt hat
...

Wir demonstrieren dies am Beispiel von S T Su = f
...
AUFLOSUNG VON FINITE{ELEMENTE{GLEICHUNGSSYSTEMEN

ST y = f

Vorwartseinsetzen :

0
B
B
B
B
@

s11 0
0
s12 s22
0

...


...


...


...

yN

1 0
C B
C B
C=B
C B
A @

f1
f2

...


...

B
...


...


...

0
sNN

10
CB
CB
CB
CB
A@

u1
u2

...


...


...

yN

1
C
C
C
C
A

uN = yN =sNN
ui =

yi ;

N
X
j =i+1

sij uj

sii

fur i = N ; 1 N ; 2 : : : 1

Bei der Losung des Gleichungssystems RRT u = f verlaufen die Schritte des Ruckwarts{ und Vorwartseinsetzens analog
...
Bei der Durchfuhrung der Faktorisierung bleiben namlich gewisse Besetztheitsstrukturen erhalten
...
Bei der Faktorisierung K = S T S
gilt

slj = 0 fur alle l l0 < j falls klj = 0 fur alle l l0 < j k(l0+1)j 6= 0
denn
s1j = k1j =s11 = 0
s2j = (k2j ; s12s1j )=s22 = 0
s3j = (k3j ; s13s1j ; s23s2j )=s33 = 0

...


...
1
...
h
...
An Positionen, an
denen Matrixelemente klj mit l0+1 < l < j identisch Null sind, entstehen in der Matrix
S im allgemeinen von Null verschiedene Elemente
...
Auf Grund der obigen Tatsache konnen wir
die Matrizen K und S auf den gleichen Speicherplatzen folgenderma en abspeichern
...
1: Schematische Darstellung der VBS{Speichertechnik
Es wird ein Vektor K aufgebaut, auf dem nacheinander die Spalten des oberen Dreiecks
der Matrix K stehen jeweils beginnend mit dem ersten Nicht{Null{Element der Spalte
...
Eine derartige Speichertechnik
bezeichnen wir als "variable Bandweite spaltenweise (VBS)\
...
AUFLOSUNG VON FINITE{ELEMENTE{GLEICHUNGSSYSTEMEN

riN

= kiN =rNN = 0

falls kiN = 0

ri(N ;1) = (ki(N ;1) ; riN r(N ;1)N )=r(N ;1)(N ;1) = 0
falls kil = 0 l = N N ; 1
ri(N ;2) = (ki(N ;2) ; ri(N ;1)r(N ;2)(N ;1) ; riN r(N ;2)N )=r(N ;2)(N ;2) = 0
falls kil = 0 l = N N ; 1 N ; 2

...


...
h
...
An Positionen, an denen
Matrixelemente kil mit i < l < l0 ; 1 identisch Null sind, entstehen in der Matrix R im
allgemeinen von Null verschiedene Elemente
...

0
1
0 0
0 0 C
B
B@
C
B
0 0 0
0 C
B @
C
B @
C
B
0 0 C
B
C
@
B
C
B
@
B
C
0 0 0 C
B
C
@
B
C
@
B
C
0
B
C
@
B
C
B
C
B symmetrisch @@
C
@
Y
HA
HH
@
H Hulle der Matrix K
@

LK

K

?

|

0 0

?

0 0 0

?

{z
Lange der Hulle

? ?

0

?

Abbildung 5
...
1
...
Gleichzeitig wird ein Hilfsvektor LK eingefuhrt, dessen Elemente auf die
Position der Hauptdiagonalelemente der Matrix im Vektor K zeigen
...

Der Aufwand an notwendigen arithmetischen Operationen bei der Durchfuhrung des
Cholesky{Verfahrens hangt o enbar von der Lange des Pro ls bzw
...
Er verhalt sich bei der Faktorisierung wie Nb2 und beim Vor{ bzw
...
Da die Cholesky{Faktoren S bzw
...

Wir konnen den Aufwand an arithmetischen Operationen und den Speicherplatzbedarf
reduzieren, wenn wir die Bandweite verringern
...
Betrachten wir zum Beispiel die folgende Vernetzung
...
3: Eine Vernetzung mit einer ungeeigneten Knotennumerierung
Bei dieser Knotennumerierung hat die FE{Matrix die folgende Besetztheitsstruktur
0
0 0 0 0 0
0 0 0 0 01
B
0 0 0 0 0
0 0 0 0C
B
C
B0
B
C
0 0 0 0 0
0 0 0C
B
C
B0 0
0 0 0 0 0
0 0C
B
C
B0 0 0
0 0 0 0 0
0C
B
C
B0 0 0 0
C
B
C
0 0 0 0 0
B
C
B0 0 0 0 0
0 0 0 0 0 0 C
B
C
B 0 0 0 0 0 0
B
C
0 0 0 0 0C
B
C
B
0 0 0 0 0
0 0 0 0C
B
C
B0
B
0 0 0 0 0
0 0 0C
C
B0 0
C
B
0 0 0 0 0
0 0C
B
C
B0 0 0
0 0 0 0 0
0C
B
C
B0 0 0 0
C
@
A
0 0 0 0 0
0 0 0 0 0
0 0 0 0 0
und somit die Bandweite b = 8
...
der Hulle ist 75
...
AUFLOSUNG VON FINITE{ELEMENTE{GLEICHUNGSSYSTEMEN

2

1

s

s

6

;
;

;
;
;

;
;
;
;

s

s

4

s

;

;
;

3

s

8

s

10

s

s 14

s

s 13

12

;
;
;
;
;
;
;
;
;
;
;
;
;
;
;
;
;
;
;
;
;
;
;
;

s

5

s

7

s

7

11

Abbildung 5
...
4 angegebenen Knotennumerierung erhalten wir die folgende
Besetztheitsstruktur
...
Fur die Lange des Pro ls
bzw
...

Es existieren Algorithmen, die fur eine gegebene Vernetzung eine Neunumerierung der
Knoten erzeugen, um die Lange des Pro ls bzw
...
B
...
Es kann erreicht werden,
da die Bandweite in der Gro enordnung von h;(m;1) liegt und sich somit der Aufwand an arithmetischen Operationen des Cholesky{Verfahren bei 2D{Problemen wie
N 2 verhalt
...
2 Iterative Verfahren
Iterationsverfahren erzeugen ausgehend von einer Startnaherung u(0) eine Folge von
Naherungslosungen u(k) fur die exakte Losung u des Gleichungssystems Ku = f
...
2
...
B
...
5)
K
Somit ist zu untersuchen, ob ke(k)kK ;! 0 fur k ;! 1 gilt, d
...
ob das Iterations-

verfahren konvergiert
...
6)
erfullt ist mit einer vorgegebenen relativen Genauigkeit "
...

Im weiteren stellen wir einige Iterationsverfahren vor
...
2
...
Ihr Nachteil ist die sehr langsame Konvergenz (siehe z
...
72])
...

Jacobi{Verfahren
Ausgehend von einer Startnaherung u(0) werden die Naherungen u(k) = u(k)]N ,
i i=1
k = 1 2 : : : wie folgt berechnet:

u(ik)

= ;

N
X
j =1
j 6=i

kij u(jk;1) + fi

!

kii fur i = 1 2 : : : N

(5
...

Gau {Seidel{Verfahren
Gegeben sei eine Startnaherung u(0)
...


...


u(ik)

=

;


...


...
8)

kii fur i = 2 3 : : : N ; 1

156

KAPITEL 5
...
h
...
Um allerdings eine Naherungslosung u(k)
mit einer relativen Genauigkeit " zu erhalten, d
...
eine Naherungslosung, fur die
ku ; u(k)kK

" ku ; u(0)kK
(5
...
Ordnung resultieren, sehr viele Iterationen durchgefuhrt werden
...
h
...
Damit liegt der Gesamtaufwand an arithmetischen Operationen, der notwendig ist, um eine Naherungslosung mit einer vorgegebenen relativen Genauigkeit "
zu erhalten, in der Gro enordnung von h;4 ln ";1
...


5
...
2 Die Methode der konjugierten Gradienten ohne Vorkonditionierung
Das Verfahren der konjugierten Gradienten wurde 1952 von Hestenes und Stiefel
40] entwickelt
...

Der Ausgangspunkt dieses Verfahrens ist die Tatsache, da die Losung des Gleichungssystems Ku = f der Minimierung eines quadratischen Funktionals aquivalent ist
...


Satz 5
...
10)

gilt
...

Die Grundidee des Verfahrens der konjugierten Gradienten besteht darin, das Minimum
des Funktionals (5
...
Dabei wird in jedem Iterationsschritt
von der Richtung des Gradienten

r := grad F (v) = Kv ; f
Gebrauch gemacht
...
Als Relaxationsrichtung s(k), in welcher das Minimum
von F (:) in jedem Iterationsschritt k, k = 2 3 : : : , zu suchen ist, wahlen wir eine
Linearkombination aus der negativen Gradientenrichtung f ; Ku(k;1) und der vorhergehenden Relaxationsrichtung s(k;1), so da (Ks(k;1) s(k)) = 0 gilt
...
Wir erhalten den folgenden Algorithmus zur Losung des Gleichungssystems Ku = f
...
2
...


Satz 5
...
Au erdem gilt
ku ; u(k) kK

mit

(k) ku ; u(0)k

K

q
q
= 2 k =(1 + 2k )
= (1 ; )=(1 + )
= 1= 2 :
Es bezeichnet j x]j die kleinste ganze Zahl, die gro er oder gleich x ist
...

(k)

Fuhren wir die Konditionszahl (K ) = 2= 1 ein, dann gilt
q
q
= ( (K ) ; 1)=( (K ) + 1) :
Bei FE{Matrizen, die aus der Diskretisierung von Randwertproblemen 2
...
Folglich wachst beim CG{
Verfahren die Iterationszahl, die notwendig ist, um das Gleichungssystem mit einer
vorgegebenen relativen Genauigkeit " zu losen, wie h;1 ln ";1
...
AUFLOSUNG VON FINITE{ELEMENTE{GLEICHUNGSSYSTEMEN

eine Matrix Vektor{Multiplikation, drei Vektoradditionen und zwei Skalarprodukte
zu berechnen sind, ist der Aufwand an arithmetischen Operationen pro Iterationsschritt
proportional zur Anzahl der Unbekannten N
...

Eine schnellere Konvergenz wurde man erhalten, wenn die Konditionszahl (K ) schwacher von h abhangen wurde
...


5
...
3 Die Methode der konjugierten Gradienten mit Vorkonditionierung
Wir betrachten anstelle des zu losenden Gleichungssystems Ku = f das aquivalente
Gleichungssystem
~~
K u = f~

(5
...
Dabei wird die Vorkonditionierungsmatrix C so gewahlt,
~
da (K )
(K ) gilt
...
11) losen wir mittels des im Abschnitt 5
...
2 beschriebenen CG{Verfahrens
...
Damit
ergibt sich der folgende Algorithmus der Methode der konjugierten Gradienten mit
Vorkonditionierung
...
2
...


Satz 5
...
Weiterhin gelte
~1(Cv v) (Kv v) ~2(Cv v) fur alle v 2 RN ~1 > 0 :
Dann sind nicht mehr als

I (") = j ln(";1 + (";2 + 1)0:5)= ln ~;1 ]j
Iterationen notwendig, um den Anfangsfehler ku; u(0)kK auf das "{fache zu reduzieren
(0 < " < 1)
...
Um dies zu erreichen, ware C = K die gunstigste Wahl, denn dann
ware ~ = 1
...
h
...

Wir mussen also Matrizen C nden, so da ~ = ~1=~2 nahe bei 1 liegt, und da die
Gleichungssysteme Cw(k) = r(k) sehr leicht losbar sind, moglichst mit einem Aufwand
an arithmetischen Operationen, der proportional zur Anzahl der Unbekannten N ist
...

(a) Vorkonditionierung mit der Hauptdiagonale der Systemmatrix
Als Vorkonditionierungsmatrix C wahlen wir die Diagonale diag(K ) der Systemmatrix K
...
h
...
Dies hat zur Folge, da
sich auch beim PCCG{Verfahren mit der Diagonalvorkonditionierung die Anzahl
der notwendigen Iterationen wie h;1 ln ";1 verhalt
...
3)
...
1 vorgestellte Warmeleitproblem, bei dem die Warmeleitzahlen zwischen 1 und 371 variieren
...
Die Matrixelemente von R werden in Analogie zum Cholesky{Verfahren nach
folgender Vorschrift berechnet (siehe auch 31, 32, 69])
...
AUFLOSUNG VON FINITE{ELEMENTE{GLEICHUNGSSYSTEMEN

p

rNN = kNN
riN = kiN =rNN

fur i = N ; 1 N ; 2 : : : 1

fur j = N ; 1 N ; 2 : : : 1

v
u
N
u
X 2
t
rjj = ukjj ;
rjl
l=j +1

8
>0
<
N
rij = >
P
: kij ; l=j+1 rilrjl

falls kij = 0

rjj

falls kij 6= 0

fur i = j ; 1 j ; 2 : : : 1

d
...
die Cholesky{Zerlegung wird nur auf den Nicht{Null{Elementen der Matrix
K durchgefuhrt
...
h
...
Es ist allerdings zu beobachten, da bei Anwendung
dieser Vorkonditionierung die zum Erreichen einer vorgegebenen relativen Genauigkeit notwendige Anzahl an Iterationen meist wesentlich niedriger ist als beim
CG{Verfahren (siehe auch 5
...

(c) MIC{Vorkonditionierung (modi ed incomplete Cholesky factorization)
Wir wahlen die Vorkonditionierungsmatrix C in der Form C = S T S mit einer
oberen Dreiecksmatrix S
...

1: sij = kij

i = 1 2 ::: N j = i i+1 ::: N

2: fur i = 1 2 : : : N ; 1
sii = psii
fur j = i + 1 i + 2 : : : N
sij = sij =sii

sjj = sjj ; s2
ij
fur l = i + 1 i + 2 : : : j ; 1
fur m = l + 1 l + 2 : : : N
t = sil sim
falls slm 6= 0 : slm = slm ; t
falls slm = 0 : sll = sll ; t und smm = smm ; t

3: sNN = psNN

161

5
...
ITERATIVE VERFAHREN

Der Algorithmus der MIC{Zerlegung unterscheidet sich von dem der IC{Zerlegung
im wesentlichen darin, da die Hauptdiagonalelemente modi ziert werden, falls
slm = 0 und silsim 6= 0 fur wenigstens ein i < l m gilt
...
Ordnung ist, dann verhalt sich (K ) wie h;0:5 im Unterschied zu h;1 fur
q
(K )
...
Damit liegt der Gesamtaufwand an arithmetischen Operationen bei 2D{Problemen in der Gro enordnung von h;2:5 ln ";1 (N 1:25 ln ";1)
...
2
...
h
...
B
...

=) : : :

=)
Vernetzung
Schrittweite
Knotenanzahl

Th1

=)

Th2

Thl

h1

>

h2

>

:::

>

hl

N1

<

N2

<

:::

<

Nl

Abbildung 5
...
4
...

Im weiteren verwenden wir anstelle der Indizes "hq\ stets nur die Indizes "q\, q =
1 2 : : : l
...
Ein Iterationsschritt des (l l ; 1){
Zweigitter{Algorithmus lauft nach folgendem Schema ab:

162

KAPITEL 5
...

l

1
...
Wahle als Iteral
tionsverfahren z
...
das Gau {Seidel{Verfahren
...

l
2
...
h
...
h
...
Nachglattung
Fuhre 2 Iterationsschritte zur naherungsweisen Losung des Gleichungssystems Kl ul = f l durch, wobei mit u(k 2) gestartet wird
...
B
...
Die neue
"
Naherung ist u(k 3)
...

l
l
Um eine gute Konvergenz des Zweigitter{ und Mehrgitteralgorithmus zu erreichen,
genugt es im allgemeinen, ein oder zwei Vor{ bzw
...
h
...

Zur Interpolation der Korrektur nutzen wir die lineare Interpolation
...
2
...


Die Einschrankung (Restriktion) im Schritt 2(a) wird nach der Beziehung
k
d(l;)1 = Ill;1d(l k)

T

Ill;1 = Ill;1

mit

durchgefuhrt
...
Da das Grobgittergleichungssystem i
...
auch noch
gro dimensioniert ist, losen wir es naherungsweise mittels Iterationsschritten eines
(l ; 1 l ; 2){Zweigitter{Algorithmus, der analog zum (l l ; 1){Zweigitter{Algorithmus
de niert ist
...
Wir setzen diese Idee
solange rekursiv fort, bis im Schritt 2(b) das Gleichungssystem K1w(k) = d(k) steht,
1
1
welches wir beispielsweise mit einem direkten Verfahren losen konnen
...

Fur = 1 bezeichnen wir einen Iterationsschritt des Mehrgitteralgorithmus als V{
Zyklus und fur = 2 als W{Zyklus
...

Seien mit die Durchfuhrung der Glattung, die Losung des Gleichungssystems auf
dem grobsten Gitter, CW die Restriktion und
die Interpolation bezeichnet, dann
erhalten wir fur ein 4{Gitter{Verfahren folgende Schemata
...
Bei vielen Anwendungen
ist kleiner 0:2
...
AUFLOSUNG VON FINITE{ELEMENTE{GLEICHUNGSSYSTEMEN

vorgegebenen relativen Genauigkeit " notwendig ist, wie ln ";1, unabhangig von der
Feinheit der Diskretisierung
...
Damit verhalt sich der
Gesamtaufwand an arithmetischen Operationen wie Nl ln ";1 (siehe auch 33, 42, 78])
...
1 Die Mehrgitterverfahren konnen auch zur Konstruktion von Vor-

konditionierungen fur den PCCG{Algorithmus eingesetzt werden
...
Dadurch wird
implizit eine Vorkonditionierungsmatrix C de niert, fur die bei geeigneter Wahl der
Glattungsverfahren, der Interpolation und der Restriktion die Beziehung
(1 ; m)(Cv v) (Kl v v) (Cv v) fur alle v 2 RNl
gilt
...
3 fur den PCCG{Algorithmus
erhalten wir somit ~1 = (1 ; m) und ~2 = 1 , wobei die Konvergenzrate des
Mehrgitterverfahrens ist (siehe 43])
...
Der Gesamtaufwand an arithmetischen Operationen wachst also wie Nl ln ";1
...
Ob das MG(m){PCCG{Verfahren oder das
Mehrgitterverfahren e ektiver hinsichtlich der benotigten CPU{Zeit ist, hangt von der
Konvergenzrate ab
...
2 Eine weitere Moglichkeit, eine Folge von Diskretisierungen in den

Losungsproze einzubeziehen, sind Algorithmen mit hierarchischen Basen
...
4
...
Die Mengen !q beinhalten die Knotennummern
der Knoten des Netzes Tq
...
Wird diese Basis zum Aufbau des FE{Gleichungssystems genutzt, dann erhalt man das Gleichungssystem
^^
Kl ul = f^l :

165

5
...
ITERATIVE VERFAHREN

^
Zwischen der Matrix Kl in der hierarchischen Basis und der Matrix Kl in der Knotenbasis sowie zwischen f^l und f l besteht der folgende Zusammenhang:
^
Kl = QT Kl Ql
l

f^l = QT f l
l

und

~~
~
mit Ql = QlQl;1 : : : Q2 ,
8 1 fur i = j i j = 1 2 : : : N
>
l
>
> 1=2 fur j = i1 und j = i2 Nq;1 < i Nq wobei die Knoten
>
<
~ q )ij =
(Q
Pi1 und Pi2 die Begrenzungsknoten jener Dreiecksseite sind,
>
auf welcher der Knoten Pi de niert wurde
>
>
: 0 sonst
...
Yserentant 87] hat bewiesen, da bei 2D{Problemen fur alle v l 2 RNl
^

^^ ^
^^ v
^^ v
c (l + 1);2 (Clvl vl) (Klvl ^l) c (Clvl ^l)
mit
^
Cl = K1 0
0 I

!

gilt, wobei c und c vom Diskretisierungsparameter unabhangige, positive Konstanten sind
...
Die Durchfuhrung dieses PCCG{Verfahrens ist mit folgenden Problemen verbunden:
^
{ Die hierarchische Stei gkeitsmatrix Kl hat wesentlich mehr Nicht{Null{Elemente
als die Stei gkeitsmatrix Kl in der Knotenbasis
...

{ Die Komponenten ul i des Losungsvektors in der hierarchischen Basis fallen fur i =
^
b
N1 +1 N1 +2 : : : Nl nicht mit den Werten der entsprechenden FE{Funktion ul 2 Vl
im Punkt Pi zusammen
...
Den so
l ^ l
entstehenden PCCG{Algorithmus bezeichnen wir als HB{PCCG{Algorithmus
...

l
l

166

KAPITEL 5
...
2
...
Deshalb ist es
ausreichend, nur die Nicht{Null{Elemente (NNE) abzuspeichern
...
Dabei schreiben wir die Nicht{Null{Elemente zeilenweise nacheinander
auf einen Vektor K
...

LKZ

ppppppp

p p

1 2

?

LKS
K

|

{z

}|

?

{z

NNE
NNE
der 1
...
Zeile

}

p Np
pppppppp
pppppppp

? ?

6

NNE der
N-ten Zeile

Abbildung 5
...

Wann bricht man die Iteration ab ?
Wie wahlt man die Startnaherung u(0) ?
Sinnvoll ist der Abbruch eines Iterationsverfahrens sicher, wenn der Fehler der Naherungslosung des Gleichungssystems in der Gro enordnung des Diskretisierungsfehlers
liegt
...
4
...
4
...
3) und somit auch eine analoge Beziehung fur
u(hk) aufgeschrieben werden kann, dann folgt

167

5
...
ITERATIVE VERFAHREN

ku ; u(k)k1 2
h

ku ; uh k1 2 + kuh ; u(k) k1 2
h

c(u) h + kuh ; u(hk)k1 2
c(u) h +

;0:5

(a(uh ; u(k) uh ; u(k)))0:5
h
h

= c(u) h +

;0:5

(Kh (uh ; u(k) uh ; u(k)))0:5
h
h

= c(u) h +

;0:5

kuh ; u(k)kKh
h

1
1

1

wobei u die exakte Losung des Randwertproblems, uh die FE{Naherungslosung und
1 die Konstante aus der Koerzitivitatsabschatzung fur die Bilinearform a(: :) (siehe
die Beziehung 4
...
4
...

(
Um den Fehler ku ; uhk)k1 2 in der Gro enordnung des Diskretisierungsfehlers zu
erhalten, mu folglich kuh ; u(k)kKh in der Gro enordnung von h liegen, d
...
wir
h
mussen das " beim Abbruchtest im Iterationsverfahren in der Gro enordnung von h
wahlen
...
Bei einer konkreten Anwendung der Iterationsverfahren ist dieser Abbruchtest
jedoch i
...
nicht auswertbar, da zum Uberprufen dieses Kriteriums die exakte Losung uh
bekannt sein mu
...
Da die Matrix C die Matrix Kh approximieren soll, ist k:kKhC;1 Kh auch eine
gute Approximation fur k:kKh
...

h

168

KAPITEL 5
...


1
...
Setze q = 1, kq = 1 und
(
uqkq ) = u1
...
Interpoliere die Losung u(kq ) auf das Netz Tq+1, d
...

q
u(0) = Iqq+1u(qkq ) :
q+1
Setze q = q + 1 und lose das Gleichungssystem

Kq uq = f q
mittels kq Schritten eines Iterationsverfahrens
...

q
Falls q < l, gehe zu Schritt 2, sonst STOP
...


5
...
Dabei werden wir demonstrieren, wie sich fur die jeweiligen
Verfahren der Aufwand an arithmetischen Operationen und der Speicherplatzbedarf bei
wachsender Dimension des FE{Gleichungssystems verhalt
...
1 durch
...
Alle Rechnungen
wurden mit dem Programmsystem FEMGP ( 42, 66, 76]) durchgefuhrt
...
Die Tests wurden auf einem PC 80486 (33 MHz) unter
Nutzung des LAHEY{FORTAN{Compilers ausgefuhrt
...
Zu jeder Vernetzung wird das entsprechende FE{Gleichungssystem

Kq uq = f q
generiert und gelost
...
7
dargestellte Triangulation
...
3
...
7: Darstellung der Vernetzung T1
Die feineren Vernetzungen erhalten wir, indem wir jeweils alle Dreiecke vierteln, d
...
wir
halbieren alle Dreiecksseiten
...


Abbildung 5
...
Die Steigkeitsmatrizen werden dafur in der Speicherform "variable Bandweite spaltenweise\
abgespeichert
...
1 sind der Speicherplatzbedarf und die benotigten
CPU{Zeiten fur die Cholesky{Zerlegung sowie fur das Vorwarts{ und Ruckwartseinsetzen zusammengestellt
...
66 kB 92
...
61 kB 5039
...
01 s
0
...
18 s
60
...
01 s
0
...
38 s
2
...
1: Rechenzeit und Speicherplatzbedarf beim Cholesky{Verfahren

170

KAPITEL 5
...
h
...
h
...
h
...
h
...

Somit bestatigen die erzielten Resultate die im Abschnitt 5
...

Im weiteren untersuchen wir verschiedene Iterationsverfahren
...
Die
4
Iteration wurde beendet, wenn das Abbruchkriterium kr(k)k " kr(0)k mit dem
Defekt r(k) = K4u(k) ; f 4 fur " = 10;4 erfullt war
...
h
...

Wir betrachten zunachst zwei klassische Iterationsverfahren, das Jacobi{Verfahren und
das Gau {Seidel-Verfahren
...
08 kB 26
...
90 s 3841
...
2: Rechenzeit und Speicherplatzbedarf beim Jacobi{Verfahren
Die Anwendung des Gau {Seidel{Verfahrens erforderte die in der Tabelle 5
...

Anzahl der Unbekannten
133
482
Speicherplatzbedarf
7
...
87 kB
Anzahl der durchgefuhrten Iterationen
29269
114319
CPU{Zeit fur die durchgefuhrten Iterationen 133
...
88 s
Tabelle 5
...
h
...
Da der Aufwand an arithmetischen Operationen pro Iterationsschritt sich wie h;2 (N ) verhalt, wachst somit
der Gesamtaufwand an arithmetischen Operationen wie h;4 ln ";1 (N 2 ln ";1), d
...
mit

5
...
EIN VERGLEICH DER AUFLOSUNGSVERFAHREN

171

dem Faktor 16
...
2 und 5
...
O enbar konvergieren sowohl das Jacobi{Verfahren als auch
das Gau {Seidel{Verfahren nur sehr langsam
...
Hatten wir
namlich beispielsweise das FE{Gleichungssystem mit 7112 Unbekannten mittels des
Gau {Seidel{Verfahrens gelost, dann ware eine CPU{Zeit von etwa 124 Stunden, also
von rund 5 Tagen, notwendig gewesen
...
Theoretisch mu te sich die Iterationszahl
bei einem vorgegebenem " wie h;1 ln ";1 verhalten, d
...
sich jeweils verdoppeln
...
4 ist ersichtlich, da sich zunachst die Iterationszahl um das 3
...
1fache und schlie lich um das 2
...
Man kann erwarten, da sich
bei weiteren Netzverfeinerungen tatsachlich das theoretisch ermittelte Wachstum der
Iterationszahl einstellt
...
h
...
Asymptotisch nimmt also der Gesamtaufwand an arithmetischen Operationen wie h;3 ln ";1 (N 1:5 ln ";1) zu
...

Anzahl der Unbekannten 133
482
1828
7112
Speicherplatzbedarf
9
...
40 kB 132
...
44 kB
Anzahl der durchgefuhrten
125
442
1377
3871
Iterationen
CPU{Zeit fur die durch0
...
98 s
75
...
23 s
gefuhrten Iterationen
Tabelle 5
...
5)
...
16 kB 34
...
95 kB 522
...
22 s
1
...
95 s
84
...
5: Rechenzeit und Speicherplatzbedarf beim DIAG{PCCG{Verfahren
Aus der Tabelle 5
...
Daraus
resultiert ein Anwachsen der Anzahl an notwendigen arithmetischen Operationen wie
h;3 ln ";1, was sich in der Zunahme der CPU{Zeit um das 8fache widerspiegelt
...
AUFLOSUNG VON FINITE{ELEMENTE{GLEICHUNGSSYSTEMEN

Ein deutlich besser konvergierendes Verfahren ist die Methode der konjugierten Gradienten mit dem Einsatz einer unvollstandigen Cholesky{Faktorisierung als Vorkonditionierer (IC{PCCG{Verfahren)
...
6
...
62 kB 48
...
36 kB 739
...
22 s
0
...
05 s
44
...
6: Rechenzeit und Speicherplatzbedarf beim IC{PCCG{Verfahren
Wie beim CG{Verfahren gilt die theoretische Aussage, da die Iterationszahl bei vorgegebenem " wie h;1 ln ";1 wachst
...
6 deutlich erkennbar
...
Zusatzlich zu dem
im CG{Verfahren benotigten Speicherplatz mu hier noch der Speicherplatz fur die
Vorkonditionierungsmatrix bereitgestellt werden
...
62 kB 48
...
36 kB 739
...
16 s
0
...
03 s
13
...
7: Rechenzeit und Speicherplatzbedarf beim MIC{PCCG{Verfahren
Bei diesem Verfahren ist theoretisch ein Anwachsen der Iterationszahl wie h;0:5 ln ";1
zu erwarten
...
7 angegebenen Iterationszahlen wachsen auch ungefahr so an
...
Der
Speicherplatzbedarf ist dergleiche wie beim IC{PCCG{Verfahren
...
In der Tabelle 5
...


5
...
EIN VERGLEICH DER AUFLOSUNGSVERFAHREN

173

Anzahl der Unbekannten
482
1828
7112
Speicherplatzbedarf
49
...
76 kB 693
...
49 s
2
...
66 s
gefuhrten Iterationen
Tabelle 5
...
h
...
Dies ist in der Tabelle 5
...
Der Aufwand an arithmetischen
Operationen wachst wie h;2 ln h;1 ln ";1 (N ln N 0:5 ln ";1)
...

Die e ektivsten Verfahren zur Losung unserer Aufgabe sind Mehrgitterverfahren
...

Anzahl der Unbekannten
482
1828
7112
Speicherplatzbedarf
32
...
30 kB 530
...
33 s
1
...
24 s
gefuhrten Iterationen
Tabelle 5
...
h
...
9 bestatigt wird
...
h
...
Gleiches gilt auch fur den Speicherplatzbedarf
...
Im Vorbereitungsschritt fur den Mehrgitteralgorithmus sind zusatzlich die FE{Gleichungssysteme fur
alle Netze zu generieren
...
72 s
...
Dies fuhrt dann auf die sogenannten MG(m){PCCG{Verfahren, wobei
m die Anzahl der durchzufuhrenden Mehrgitterschritte bei der Losung der Vorkonditionierungsgleichungssysteme ist
...


174

KAPITEL 5
...
74 kB 190
...
90 kB
Anzahl der durchgefuhrten
3
3
3
Iterationen
CPU{Zeit fur die durch0
...
70 s
8
...
10: Rechenzeit und Speicherplatzbedarf beim MG(1){PCCG{Verfahren
Auch beim MG(1){PCCG{Verfahren ist die Iterationszahl bei vorgegebenem " unabhangig vom Diskretisierungsparameter h
...
Der Speicherplatzbedarf ist proportional zur
Anzahl der Unbekannten
...
a
...

Zum Abschlu vergleichen wir nochmals alle Verfahren hinsichtlich der benotigten
CPU{Zeit und des Speicherplatzbedarfs
...
02 s

0
...
08 kB

26
...
66 kB
92
...
90 s 3841
...
08 s

CG{Verfahren
DIAG{PCCG{Verfahren
IC{PCCG{Verfahren
MIC{PCCG{Verfahren
HB{PCCG{Verfahren
Mehrgitter{Verfahren

1828
7112
4
...
28 s
661
...
52 kB

7
...
55 s
9
...
22 s
9
...
22 s
12
...
16 s
12
...
88 s

26
...
98
34
...
15
34
...
72
48
...
44
48
...
49
49
...
75 s 1016
...
95 kB 522
...
95 s
84
...
95 kB 522
...
05 s
44
...
36 kB 739
...
03 s
13
...
36 kB 739
...
63 s
14
...
76 kB 693
...
33 s

32
...
59 s

8
...
30 kB 530
...
33 s
1
...
40 s
MG(1)-PCCG{Verfahren
47
...
42 kB 752
...
11: Vergleich aller vorgestellten Au osungsverfahren

5
...
EIN VERGLEICH DER AUFLOSUNGSVERFAHREN

175

Die Tabelle 5
...
Zur Losung von
gro dimensionierten Gleichungsystemen sollten nur Mehrgitterverfahren oder PCCG{
Verfahren eingesetzt werden
...

In der Abbildung 5
...


Abbildung 5
...
AUFLOSUNG VON FINITE{ELEMENTE{GLEICHUNGSSYSTEMEN

Kapitel 6

Galerkin{FEM fur parabolische
Anfangs{Randwertaufgaben
6
...
2
...


fur alle (x t) 2 QT
QT =
(0 T )
fur alle x 2 ;1 t 2 0 T ]
(6
...
Um diese zu erhalten, multiplizieren wir
die Di erentialgleichung und die Anfangsbedingung aus der Aufgabe (6
...
Damit erhalten wir
Z @u X @
2
@u v(x) dx = Z f (x t)v(x) dx :
@t ; i=1 @xi i(x t) @xi
Die partielle Integration im elliptischen Anteil liefert
Z @u
2
X
@u @v dx ; Z @u v(x) ds = Z f (x t)v(x) dx :
@t v(x) + i=1 i (x t) @xi @xi
@N
@

178

KAPITEL 6
...
1) lautet somit:
Gesucht ist u(x t) 2 Vg1 mit u2 L2( ) fur fast alle t 2 (0 T ), so da
(u v)0 + a(t u v) = hF (t) vi
(u(x 0) v)0 = (u0 v)0
mit
(u v)0

=

Z

fur alle v 2 V0 und fur fast alle t 2 (0 T )
fur alle v 2 V0
(6
...

Das Problem (6
...
Um eine
Diskretisierung der Aufgabe (6
...
Dazu verwenden wir den folgenden Losungsansatz
X
X
uh = uh(x t) = ui(t)pi(x) + u i(t)pi(x)
i2!h

i2

h

wobei wir hier wieder die im Abschnitt 4
...
1 eingefuhrten Mengen !h und
Die Koe zienten u i(t) sind durch u i(t) = g1(xi t) de niert
...


179

6
...
DIE STETIGE, DIE SEMIDISKRETE UND DIE VOLLDISKRETE AUFGABE

Damit erhalten wir
Gesucht ist uh(x t) 2 Vg1 h mit uh 2 L2( ) fur fast alle t 2 (0 T ), so da
(uh vh)0 + a(t uh vh) = hF (t) vhi fur alle v 2 V0h und fur fast alle t 2 (0 T )
(uh(x 0) vh)0 = (u0 vh)0 fur alle v 2 V0h
(6
...

Beachten wir die konkrete Gestalt der Funktion uh(x t), und setzen wir fur vh(x) die
Funktionen pk (x), k 2 !h , nacheinander ein, dann folgt:
Gesucht ist uh (t) = ui(t)]i 2 !h , so da fur alle k 2 !h und fur fast alle t 2 (0 T )
X
X
X
ui (t)(pi pk )0 +
ui(t)a(t pi pk ) = hF (t) pk i ;
u i(t)a(t pi pk )
i 2 !h
i 2 !h
i2 h
X
u i (t)(pi pk )0
;
i2 h
(6
...
h
...


i2

= (u0 pk )]k 2 !h

h

;

i

i k 0 k 2 !h

X
i2

h

u i(0)(pi pk )0

(6
...
GALERKIN{FEM FUR PARABOLISCHE ANFANGS{RANDWERTAUFGABEN

Zur Bestimmung der Vektorfunktion uh (t) = ui(t)]i2!h haben wir das System (6
...
Die Aufgabe (6
...
Zur
naherungsweisen Losung dieses Systems gewohnlicher Di erentialgleichungen nutzen
wir ein {gewichtetes Di erenzenschema
...
h
...
Wir wahlen hier der Einfachheit
halber eine feste Zeitschrittweite
...
Als Ergebnis erhalten wir die volldiskrete
Ersatzaufgabe ( {gewichtetes Schema)
...
6)

gilt mit 'j = f h(tj ) + (1 ; )f h(tj;1) sowie die Anfangsbedingung
erfullt wird
...
7)

1
Fur = 0 hei t das Schema (6
...

Aus der Aufgabe (6
...
8)
h
h
zu losen ist
...


Bemerkung 6
...
Diese Vorgehensweise ist aus der Literatur als "mass{lumping\ bekannt
...
8) fur = 0 auf
Grund der Diagonalgestalt von Dh besonders einfach losbar
...
2 Konvergenz und Stabilitat
Eine wichtige Frage bei derartigen Schemata ist die Frage nach der Stabilitat
...
Fur = 0 erhalten

6
...
KONVERGENZ UND STABILITAT

181

wir nur ein bedingt stabiles Schema
...
Fur
2 ist das Schema (6
...
h
...

Zum Abschlu geben wir zwei Fehlerabschatzung an
...
B
...

Falls u(x t) hinreichend glatt ist und wir stuckweise lineare Ansatzfunktionen bei der
FE{Diskretisierung bezuglich des Ortes nutzen, dann gilt in der H 1{ bzw
...


182

KAPITEL 6
...
A
...
Sobolev Spaces
...

2] Th
...
Finite{Elemente{Methoden uber lokal verfeinerten Netzen fur elliptische Probleme in Gebieten mit Kanten
...

3] J
...
Argyris
...
Aircraft Engineering,
27:125{154, 1955
...
Axelsson und V
...
Barker
...
Academic Press, Orlando Fla
...

5] I
...
C
...
Error estimates for adaptive nite element
computations
...
Num
...
, 15(4):736{754, 1978
...
K
...
Butter eld
...
McGraw{Hill Book Company, 1981
...
E
...
Weiser
...
Mathematics of Computation, 44(170):283{301, 1985
...
Bernhardt
...
Jahresarbeit, Technische
Universitat Karl{Marx{Stadt, Sektion Mathematik, 1988
...
Braess
...
Springer Lehrbuch, 1992
...
H
...
E
...
Xu
...
Mathematics of Computation, 55(191):1{22, 1990
...
Breitschuh und R
...
Die Finite{Element{Methode
...

12] S
...
R
...
The Mathematical Theory of Finite Element Methods
...

13] P
...
The nite element method for elliptic problems
...

14] R
...
Variational methods for the solution of problems of equilibrium and
vibrations
...

183

184

LITERATURVERZEICHNIS

15] A
...
Cybenko, N
...
Vascenko, N
...
Kriscuk und Ju
...
Lavendel
...
Golovnoe izdatel'stvo izdatel'skogo obedinenija "visca skola", Kiev, 1986
...
Numerische Methoden der Mechanik
...

17] N
...
Grundlagen der technischen Thermodynamik
...

18] G
...
Fichtenholz
...
3
...

19] G
...
Fichtenholz
...
2
...

20] U
...
Finite{Elemente{Programme in der Festkorpermechanik
...

21] K
...
Friedrichs
...
Technical report, N
...
Univ
...

22] R
...
Gallagher
...
Prentice{Hall, Englewood Cli s, New Jersey, 1975
...
George und J
...
Lui
...
Prentice Hall, Englewood Cli s, New Jersey, 1981
...
Girault und P
...
Raviart
...
Springer{Verlag, Berlin, 1979
...
Goldner
...
Band 1
...

26] H
...
Lehrbuch Hohere Festigkeitslehre
...
Fachbuchverlag, Leipzig{
Koln, 1992
...
Goring, H
...
Roos und L
...
Finite{Element{Methode, Wissenschaftliche Taschenbucher, Bd
...
Akademie{Verlag, Berlin, 1985
...
Griebel
...

Teubner Skripten zur Numerik
...
G
...

29] D
...
Gri ths, editor
...
The
Clarendon Press, Oxford University Press, New York, 1984
...
Gro mann und H
...
Roos
...
Teubner Studienbucher Mathematik
...

31] I
...
A class of rst order factorization methods
...

32] I
...
On modi ed incomplete factorization methods
...
Springer{Verlag, Berlin, 1982
...
Hackbusch
...
Springer{Verlag, Berlin, 1985
...
Hackbusch
...
Teubner
Studienbucher Mathematik
...

35] W
...
Integralgleichungen: Theorie und Numerik
...
Teubner{Verlag, Stuttgart, 1989
...
Hackbusch
...

Teubner Studienbucher Mathematik
...

37] G
...
Kozik
...
BSB B
...
Teubner Verlagsgesellschaft, Leipzig, 1971
...
Heise
...
Diplomarbeit, Technische Universitat Karl{Marx{Stadt, Sektion Mathematik, 1987
...
Heise
...
Wissenschaftliche Schriftenreihe 4/1991, Technische Universitat
Chemnitz, 1991
...
R
...
Stiefel
...
J
...
Nat
...
Standards, 49:409{436, 1952
...
Hinton und D
...
J
...
An Introduction to Finite Element Computations
...

42] M
...
Eine Einfuhrung in die Theorie und Anwendung von Mehrgitterverfahren
...

43] M
...
Langer, A
...
Queck und M
...
Multigrid preconditioners and their applications
...
Telschow, editor, Third Multigrid Seminar, Biesenthal 1988, 11{52, Berlin, 1989
...
Report
R{MATH{03/89
...
Kikuchi
...
Cambridge University Press,
1986
...
G
...
Sopostavlenie metoda konecnych elementov c variacionno{
raznostnym metodom resenija zadac teorii uprugosti
...
VNIIG im
...
E
...

46] V
...
Korneev
...

Izdatel'stvo Leningradskogo Universiteta, Leningrad, 1977
...
G
...
Langer
...
69 Teubner{Texte zur Mathematik
...

48] M
...
Neittaanmaki
...
Pitman Monographs and Surveys in Pure and Applied
Mathematics 50
...
), New York, 1990
...
D
...
M
...
Lehrbuch der theoretischen Physik
...
2
...

50] B
...
Topics in Finite Element Solution of Elliptic Problems
...

51] A
...
Mitchell und R
...
The nite element method in partial di erential
equations
...

52] J
...
Hlavacek
...
Elsevier, SNTL, 1981
...
H
...
de Vries
...
Academic Press, New York San Francisco London, 1978
...
T
...
Finite elements of nonlinear continua
...

55] J
...
Oden, E
...
Becker und G
...
Carey
...
Prentice Hall, Englewood Cli s, N
...
, 1982
...
T
...
F
...
Finite Elements: A Second Course, The Texas Finite
Element Series II
...
J
...

57] J
...
Oden und G
...
Carey
...
Prentice Hall, Englewood Cli s, N
...
, 1982
...
T
...
F
...
Finite Elements: Fluid Mechanics, The Texas Finite
Element Series VI
...
J
...

59] J
...
Oden und G
...
Carey
...
Prentice Hall, Englewood Cli s, N
...
, 1982
...
T
...
F
...
Finite Elements: Special Problems in Solid Mechanics,
The Texas Finite Element Series V
...
J
...

61] J
...
Oden und J
...
Reddy
...
John Wiley & Sons, New York, 1976
...
A
...
Schodomost' raznostnych schem pri ulucsennoj approksimacii
granicy
...
AN SSSR, 170(1):41{44, 1966
...

63] L
...
Oganesjan und L
...
Ruchovec
...
Izd
...
Nauk Armjanskoj SSR, Erevan, 1979
...
Oswald
...

Teubner Skripten zur Numerik
...
G
...

65] N
...
Das Demonstrationsprogramm FEM2D
...

66] W
...
FEMGP { Finite-Element-Multi-Grid-Package { Programmdokumentation und Nutzerinformation
...


LITERATURVERZEICHNIS

187

67] K
...
Reid
...
J
...
Math
...
, 9:1{13, 1972
...
Ritz
...
J
...
, 135, 1908
...
Rjasanow
...
Preprint 15, Technische Universitat Karl-Marx-Stadt, 1986
...
A
...
Theorie der Di erenzenverfahren
...

71] A
...
Samarskij und A
...
Gulin
...
Nauka, Moskva, 1989
...
A
...
S
...
Numerical Methods for Grid Equations
...

II: Iterative Methods
...

73] Schellbach
...
Journal der Mathematik, 41:293{
363, 1852
...

74] H
...
Methode der niten Elemente
...

75] H
...
FORTRAN{Programme der Methode der niten Elemente
...

76] T
...
Jung
...
90)
...

77] G
...
Fix
...
Prentice{Hall
Inc
...

78] K
...
Trottenberg
...
In W
...
Trottenberg,
editors, Multigrid Methods, Proceedings of the Conference held at Koln{Porz, November 23{27, 1981, Lecture Notes in Mathematics 960, 1{176, Berlin{Heidelberg{
New York, 1982
...

79] M
...
Das Demonstrationsprogramm FEM1D
...

80] V
...
Galerkin Finite Element Methods for Parabolic Problems, Lecture Notes in Mathematics 1054
...

81] M
...
Turner, R
...
Clough, H
...
Martin und L
...
Topp
...
J
...
Sci
...

82] V
...
Sajdurov
...
Nauka, Moskva,
1989
...
Verfurth
...
Technical report, Institut fur Angewandte Mathematik, Universitat Zurich, 1993
...
Zenisek
...
Academic Press, London, 1990
...
R
...
The Mathematics of Finite Elements and Applications
IV
...

86] J
...
Partielle Di erentialgleichungen { Sobolevraume und Randwertaufgaben
...

87] H
...
On the multi{level splitting of nite element spaces
...

88] E
...
Vorlesungen uber nichtlineare Funktionalanalysis II: Monotone Operatoren
...
B
...
Teubner Verlagsgesellschaft, Leipzig,
1977
...
Zhang
...
Numerische Mathematik, 63:521{539,
1992
...
C
...
The Finite Element Method in Structural and Continuum
Mechanics
...
McGraw{Hill, New York, 1967
...
C
...
The Finite Element Method in Engineering Science
...

92] O
...
Zienkiewicz
...
Fachbuchverlag, Leipzig,
1975
...
C
...
The Finite Element Method, 3rd ed
...

94] O
...
Zienkiewicz, W
...
Kelly, J
...
Babuska
...
In J
...
Whiteman, editor,
The Mathematics of Finite Elements and Applications IV, 313{346
...

95] O
...
Zienkiewicz und K
...
Finite Elements and Approximation
...

96] O
...
Zienkiewicz und R
...
Taylor
...
Volume 1: Basic
formulation and linear problems
...
McGraw{Hill Book Company, New York, 4 edition, 1991
...
C
...
Z
...
A simple error estimator and adaptive procedure
for practical engineering analysis
...
J
...
Meth
...
, 24:337{357, 1987
...
Zlamal
...
Numerische Mathematik, 12(5):394{408,
1968

---