Die Methode der finiten Elemente

Kapitel 15
Die Methode der finiten
Elemente
Die Darstellung der Methode der FE stützt sich im wesentlichen auf die
Kapitel II und III des Buches über FE von D. Braess.
15.1
Erläuterung der Grundideen anhand des MP’s
Bevor wir die Methode der finiten Elemente systematisch einführen, wollen
wir vorab schon die Grundideen dieses Ansatzes anhand des Modellproblems
−∆u = f
Ω ⊂ IR2
, u = 0 (∂Ω)
(15.1)
erläutern.
1.Grundidee: Die erste Grundidee besteht darin, der Randwertaufgabe
eine Variationsaufgabe zuzuordnen und zwar in dem Sinne, dass die (noch
zu spezifizierende) Lösung des Variationsproblems gleichzeitig auch Lösung
der Randwertaufgabe ist. Für das Modellproblem betrachten wir hierzu das
Variationsfunktional
Z
Z
J [u] :=
u2x + u2y dΩ − 2 uf dΩ
Ω
Ω
auf dem Raum D0 := u ∈ C2 (Ω) ∩ C Ω : u = 0 auf ∂Ω . Gesucht ist nun
das Minimum u ∈ D0 , so dass J [u∗ ] ≤ J [u] für alle u ∈ D0 .
Satz 15.1 (RWA ↔ Variationsproblem). Falls u∗ ∈ D0 obige Variationsaufgabe löst, dann ist u∗ auch eine klassische Lösung der Randwertaufgabe (15.1).
102
Die Methode der finiten Elemente
103
Beweis: u∗ löse die Variationsaufgabe. Betrachte u∗ +tv mit t ∈ IR, v ∈ D0 .
Dann gilt:
Z
Z
∗
∗
2
2
2
J [u ] ≤ J [u + tv] = t
vx + vy dΩ + t 2
u∗x vx + u∗y vy dΩ
Ω
Ω
Z Z
Z
∗2
∗2
−2 f vdΩ +
ux + uy dΩ − 2 f u∗ dΩ
Ω
Ω
Ω
Dabei sei v fest und t beliebig. Dann nimmt die Funktion Φv (t) = J [u∗ + tv]
für t = 0 ihr Minimum an. Damit ergibt sich
Z
Z
∂Φv (0)
∗
∗
ux vx + uy vy dΩ − f vdΩ = 0 =
.
(15.2)
∂t
|Ω
{z
} | Ω {z }
=:a(u∗ ,v)
=:hf,vi
Partielle Integration liefert dann
Z
Z
Z
∗
∗
−
v uxx + uyy dΩ +
v . . . − f vdΩ = 0
Ω
∂Ω
Ω
| {z }
=0
Z
⇒
v (−∆u − f ) dΩ = 0
∀v ∈ D0 .
Ω
Wäre jetzt −∆u − f 6= 0, also z.B. > 0 an einer Stelle (x∗ , y ∗ ), dann wäre
−∆u − f > 0 auch in einer Umgebung von (x∗ , y ∗ ) aufgrund von Stetigkeitsargumenten. Man kann dann v so wählen, dass der Integrand in dieser
Umgebung größer als Null und außerhalb der Umgebung gleich Null ist. Dies
führt zu einem Widerspruch. Folglich ist
−∆u(x, y) = f (x, y)
((x, y) ∈ Ω) .
Die Randbedingungen werden offensichtlich wegen u ∈ D0 automatisch
erfüllt.
2.Grundidee: Im zweiten Schritt drückt man die Variationsaufgabe mithilfe einer Bilinearform aus, siehe (15.2): Falls
a (u∗ , v) = hf, vi
∀v ∈ D0 ,
(15.3)
dann löst u∗ auch die RWA (15.1). Diese Formulierung der Aufgabe mithilfe einer geeigneten Bilinearform bezeichnet man als schwache Formulierung der RWA. Bezüglich der Formulierung (15.3) reicht offensichtlich, dass
Die Methode der finiten Elemente
104
u∗ ∈ C1 (Ω) im Gegensatz zu einer klassischen Lösung u∗ ∈ C2 (Ω) der RWA
(15.1). Eine derartige Abschwächung des Lösungsbegriffs wird systematisch
in Abschnitt 15.3 behandelt.
3.Grundidee: Als nächstes ersetzt man den unendlichdimensionalen Raum
D0 durch einen endlichdimensionalen Raum D0M = span {ψ1 , . . . , ψM }. Dies
bedeutet also nichts anderes als eine
P Diskretisierung der Variationsaufgabe,
indem wir u∗ ∈ D0 durch u∗M = M
i=1 ci ψi approximieren. Hierbei werden
die Koeffizienten ci mithilfe der M Gleichungen
a (u∗M , ψk ) = hf, ψk i
⇔
M
X
ci a (ψi , ψk ) = hf, ψk i
(k = 1, . . . , M )
i=1
bestimmt. Dies liefert ein lineares Gleichungssystem
 
c1
 .. 
A .  = b
cM
zur Berechnung der ci und zwar mit
A = (aik )i,k=1,...,M = (a (ψi , ψk ))i,k=1,...,M
und b = (hf, ψk i)k=1,...,M .
Bemerkung 15.1 (Galerkin-, Ritz-Verfahren). Die unter 3. beschriebene Vorgehensweise nennt man Galerkin-Verfahren. Bei inhomogenen Randbedingungen wählt man zusätzlich eine Funktion ψ0 , die die inhomogenen
Randbedingungen erfüllt:
u∗M = ψ0 =
M
X
ci ψi .
i=0
Dies führt zu einer Modifikation der rechten Seite b. Falls die Bilinearform
a spd ist, ergibt sich das sogenannte Ritz-Verfahren. Die Verwendung dieser
Begriffe ist jedoch nicht einheitlich in der Literatur. Häufig spricht man auch
allgemein vom Ritz-Galerkin-Verfahren unabhängig davon, ob a spd ist oder
nicht.
I
Klassische Ansatzfunktionen zur Konstruktion von D0M :
1. Polynome
2. Produkte von Orthogonalpolynomen
Die Methode der finiten Elemente
105
3. Produkte trigonometrischer Funktionen
4. Produkte von Spline-Funktionen
4.Grundidee. Bei der Methode der finiten Elemente wählt man die Basis/Ansatzfunktionen ψk derart, dass das resultierendes GS dünnbesetzt ist.
Hierzu
• unterteilt man Ω näherungsweise in Dreiecke und/oder Vierecke und
• verwendet Ansatzfunktionen, die nur lokal (d.h. nur auf wenigen Elementen der Zerlegung) ungleich Null sind, aber global (d.h. über Elementgrenzen hinweg) wohldefinierte Glattheitsbedingungenen erfüllen.
Insbesondere seien die ψk durch Polynome niedrigen Grades darstellbar.
Beispiel 15.1 (Lineare FE, reguläre Triangulierung, MP).
Als spezielles Beispiel betrachten wir das Modellproblem auf dem Einheitsquadrat mit einer gleichmäßigen regulären Triangulierung und linearen finiten Elementen:
−∆u = f
Ω = (0, 1)2
u=0
(∂Ω)
D0M := v ∈ C Ω , v ist in jedem Element linear, und v = 0 auf ∂Ω .
In jedem Dreieck hat v ∈ D0M die Form v (x, y) = a + bx + cy und ist durch
die Werte in den 3 Ecken der Dreiecke eindeutig bestimmt. v ist somit global
durch die Werte
an den M Gitterpunkten festgelegt und es ergibt sich somit:
M
dim D0 = M .
Eine geeignete Basis (d.h. eine Basis, die möglichst viele Nulleinträge
in der resultierenden Matrix A nachsichzieht) stellt die sogenannte nodale
Basis dar, für die gilt:
ψk (x` , y` ) = δk`
für k = 1, . . . , M .
Die Methode der finiten Elemente
11111
00000
00000
11111
00000
11111
00000
11111
N
NW
2.
W 3.
Z 1. O
4.
6.
5.
SO
S
106
Die Ansatzfunktion ψk ist offensichtlich nur im schraffierten Bereich der
nebenstehenden Zeichnung ungleich
Null. Hierbei sei Z = Zk = (xk , yk ) =
h(ik , jk ). Derartige Ansatzfunktionen
werden auch als Hutfunktionen bezeichnet.
h
Wir wollen nun die Darstellung von ψk in jedem der sechs schraffierten Dreiecke berechnen. Hierzu müssen wir die zugehörigen Koeffizienten a, b, c der
allgemeinen Darstellung v(x, y) = a + bx + cy der linearen finiten Elemente
in jedem der sechs Dreiecke bestimmen. Für die Berechnung der Matrixeinträge a(ψi , ψk ) benötigen wir die entsprechenden Gradienten von ψk , die wir
deshalb direkt mit angeben.
1. Dreieck:
a + bhik + chjk = 1 ,
a + bh (ik + 1) + chjk = 0 ,
a + bhik + ch (jk + 1) = 0 .
x y
ψk (x, y) = 1 + ik + jk − −
.
h h
1
1 T
⇒∇ψk (x, y) = − , −
.
h h
⇒
2. Dreieck:
a + bhik + chjk = 1 ,
a + bhik + ch (jk + 1) = 0 ,
a + bh (ik − 1) + ch (jk + 1) = 0 .
⇒
3. Dreieck: . . .
y
.
ψk (x, y) = 1 + jk −
h
1 T
⇒∇ψk (x, y) = 0, −
.
h
x
ψk (x, y) = 1 − ik +
h
⇒
∇ψk =
1
,0
h
T
.
Die Methode der finiten Elemente
4. Dreieck: . . .
107
x y
ψk (x, y) = 1 − ik − jk + +
h h
6. Dreieck: . . .
x
ψk (x, y) = 1 + ik −
h
∇ψk =
y
ψk (x, y) = 1 − jk +
h
5. Dreieck: . . .
⇒
⇒
∇ψk =
⇒
∇ψk =
1
0,
h
1 1
,
h h
T
T
.
T
1
− ,0
.
h
Schließlich bestimmen wir die Einträge von AM mit
Z
a (ψi , ψk ) =
∇ψi ∇ψk dx .
Ω
Z
Z
∇ψz ∇ψz dx = 2
∇ψz ∇ψz dx
Ω
1.+2.+3.
Z
Z
Z
2
1
1
=2
dx +
dx +
dx
2
2
2
1. h
2. h
3. h
2 h2
1 h2
1 h2
=2 2
+ 2
+ 2
=4
h 2
h 2
h 2
Z
Z
a(ψZ , ψN ) =
∇ψz ∇ψN dx =
∇ψz ∇ψN dx .
a(ψZ , ψZ ) =
Ω
1.+2.
Z 1
1 T
1 T
1 T 1 1 T
=
− ,−
0,
dx +
0, −
,
dx
h h
h
h
h h
1.
2.
Z
Z
1
1 h2
1 h2
1
− 2 dx = − 2
− 2
=
− 2 dx +
h
h
h 2
h 2
2.
1.
= −1 = a (ψZ , ψS ) = a (ψZ , ψW ) = a (ψZ , ψO )
Z
Z
a(ψZ , ψN W ) =
∇ψZ ∇ψN W dx =
∇ψz ∇ψN w dx .
Z Ω
Z 2.+3.
T 1
1
=
0, −
− ,0
h
h
2.
= 0 = a (ψZ , ψSO ) .
T
Z dx +
3.
1
,0
h
T 1
0,
h
T
dx
Alle anderen Einträge mit nicht benachbarten ψk sind offensichtlich ebenfalls
Null. In Sternschreibweise ergibt sich somit dasselbe lineare Gleichungssystem, wie wir es schon bei den FD und den FV kennen gelernt haben:


−1
AM =
b −1 4 −1 .
−1
h
.
Die Methode der finiten Elemente
108
Dieser Zusammenhang zwischen den verschiedenen Diskretisierungsansätzen
ist jedoch lediglich für dieses einfache Modellproblem charakteristisch und
läßt sich auf allgemeinere Fälle nicht übertragen.
J
15.2
Allgemeinere Aufgaben, weitere wichtige Beispiele
15.2.1
Inhomogene Dirichlet-Randbedingungen
Häufig lassen sich RWA’en mit inhomogenen Dirichlet-RB’en auf RW’en mit
homogenen Dirichlet-RB’en zurückführen. Man Betrachte den elliptischen
Differentialoperator 2. Ordnung
Lu = −
N
X
∂i (aik ∂k u) + a0 u
i,k=1
mit a0 (x) ≥ 0 für x ∈ Ω. Ziel ist es das Dirichletproblem
Lu = f
u=g
(x ∈ Ω)
(x ∈ ∂Ω)
auf ein solches mit homogenen RB’en zu transformieren. Hierzu sei eine
zulässige Funktion u0 bekannt mit u0 = g auf ∂Ω. Daraus ergibt sich die
modifizierte RWA
Lw = f1
(x ∈ Ω)
w=0
(x ∈ ∂Ω)
mit w = u − u0 , f1 = f − Lu.
Insbesondere für theoretische Zwecke wird daher oft angenommen, dass
die Aufgabe mit homogenen Randbedingungen gegeben ist.
15.2.2
Poisson-Gleichung mit Neumann-Randbedingungen
Betrachte die RWA
−∆u = f
un = 0
(Ω)
(∂Ω) .
Wir setzen voraus, dass ∂Ω hinreichend glatt ist (also z.B. ∈ C∞ ). J[u] sei
b = C2 (Ω) ∩
wie in Abschnitt 15.1 definiert, wobei anstelle von D0 jetzt D
Die Methode der finiten Elemente
109
b keine Randbedingungen
C1 (Ω) betrachtet wird. Man beachte, dass in D
auftauchen. Dafür wird jetzt einmal stetige Differenzierbarkeit bis zum Rand
gefordert.
Der Gedankengang ist wieder ähnlich wie in Abschnitt 15.1, mit
Φv (t) = J[u∗ + tv]
mit v ∈ C2 (Ω) .
Damit erhält man dann
Z
Z
∗
∗
f v d(x, y) = 0 ,
(ux vx + uy vy ) d(x, y) −
Ω
Ω
und hieraus folgt durch partielle Integration
Z
Z
∗
(−∆u − f )v d(x, y) +
u∗n v do(x, y) = 0
Ω
für alle v ∈ C2 (Ω) ,
∂Ω
0 (Ω). Mit dem gleichen Schluß wie oben folgt
insbesondere für alle v ∈ C∞
−∆u = f
(Ω) .
Ein ähnlicher Schluß für das Randintegral (die Annahme un 6= 0 führt man
durch geeignete Wahl von η zu einem Widerspruch) liefert
un = 0
(∂Ω) .
b erfüllt zu
Offenbar brauchen hier die Randbedingungen nicht im Raum D
sein; sie ergeben sich automatisch“. Es sind natürliche“ Randbedingungen.
”
”
Die allgemeine Definition von wesentlichen und natürlichen Randbedingungen erfolgt später.
15.2.3
Minimalflächengleichung
Betrachte die nichtlineare RWA
(1 + (ux )2 )uxx − 2ux uy uxy + (1 + (uy )2 )uyy = 0
u=g
(Ω)
(∂Ω)
und das zugehörige Variationsfunktional
Z q
J[u] =
1 + (ux )2 + (uy )2 d(x, y) .
Ω
Dieses Integral beschreibt den Flächeninhalt der durch u beschriebenen
Fläche; er wird minimal für die Lösung der Minimalflächengleichung. Physikalisch beschreibt die RWA eine elastische Membran, die in eine Raumkurve
Die Methode der finiten Elemente
110
g eingespannt wird. Sie stellt sich so ein, dass die eingespannte Fläche minimal wird.
Bei einer eingespannten elastischen Membran unter einer Belastung f ist
das Funktional
Z q
Z
2
2
J[u] = ( 1 + (ux ) + (uy ) − 1) d(x, y) −
f u d(x, y) .
Ω
Ω
Dabei beschreibt das zweite Integral die Arbeit der äußeren Kräfte, die −1“
”
im ersten Integral repräsentiert die Flächenänderungsenergie.
Der Zusammenhang zur Poissongleichung wird deutlich, wenn man (für
kleine Auslenkungen ux , uy ) eine Taylorentwicklung durchführt. In erster
Näherung ist
q
1 + (ux )2 + (uy )2 − 1 =
˙ (ux )2 + (uy )2 ,
d.h. die Poissongleichung stellt eine globale Linearisierung dar.
15.2.4
Biharmonische Gleichung
Die DGL
∆∆u = uxxxx + 2uxxyy + uyyyy = f
beschreibt die Auslenkung einer elastischen dünnen Platte unter einer Belastung f (bei kleinen Auslenkungen).
Die Randbedingungen
u = g1
(∂Ω)
un = g2
(∂Ω)
(15.4)
beschreiben eine eingespannte Platte (u beschreibt die Lage, un die Richtung
der Platte). Man benötigt hier zwei Randbedingungen, da die DGL von
vierter Ordnung ist.
Die Randbedingungen
u = g1
(∂Ω)
∆u = 0
(∂Ω)
(15.5)
beschreiben eine frei aufliegende Platte (∆u beschreibt das Drehmoment).
Das zugehörige Variationsfunktional lautet
Z
J[u] = ((∆u)2 − 2f u) d(x, y) .
Ω
Die Methode der finiten Elemente
111
Dabei sind im Falle von (15.4) beide Randbedingungen wesentlich und
müssen im Raum Dg berücksichtigt werden. Im Falle von (15.5) ist nur
die Bedingung u = g1 eine wesentliche Randbedingung, −∆u = 0 ist eine
natürliche Randbedingung (ergibt sich automatisch).
Allgemein haben bei einer DGL 2m-ter Ordnung die Randbedingungen
eine Ordnung ≤ 2m − 1. Dabei sind die Ableitungen mit einer Ordnung
≤ m − 1 wesentliche Randbedingungen.
15.3
Variationsformulierung elliptischer RWA’en
Eine systematische Verallgemeinerung der ersten beiden Grundideen aus
Abschnitt 15.1 führt auf den Begriff der schwachen Lösung (Lösung der Variationsaufgabe bzw. der Formulierung mit Bilinearform). Hierbei wird eine
Abschwächung der Glattheitsvoraussetzungen an die Lösung vorgenommen
(vgl. Beispiel 15.1). Insgesamt ergeben sich somit deutlich schwächere Voraussetzungen für die Existenz von Lösungen als bei der Methode der FD
(bzw. FV).
15.3.1
Sobolevräume
Grundlegend für den Begriff der schwachen Lösung sind die sogenannten
Sobolevräume. Sobolevräume basieren auf dem L2 (Ω), der mit dem Skalarprodukt
Z
hu, viL2 =: hu, vi0 =
uv dΩ
Ω
und der zugehörigen Norm
kuk0 =
q
hu, vi0
zu einem Hilbertraum wird.
Sobolevräume fassen Funktionen mit sogenannten schwachen Ableitungen zusammen. Der Begriff der schwachen Ableitung wird durch die partielle
Integration motiviert: Für u, φ ∈ C1 [a, b] gilt
Z b
Z b
0
uφ dx = −
u0 φ dx + uφ|ba .
(15.6)
a
a
0 [a, b] gilt, dass
Falls für u, v ∈ L2 [a, b] für alle φ ∈ C∞
Z b
Z b
0
uφ dx = −
vφ dx
a
a

Download Report