Numerische Simulation mit finiten Elementen [-.2em] - O ... · Variationsformulierung, „Methode...

Einleitung

INHALT1. Finite Differenzen in 2D2. Einleitung2.1 Vorbemerkungen2.2 Rand- und Anfangswertaufgaben3. Variationsformulierung von Randwertaufgaben3.1 Vorbemerkungen3.2 Bezeichnungen und mathematische Hilfsmittel3.3 Referenzaufgaben in 2D3.4 Variationsformulierung der Poisson-Gleichung3.5 Galerkin-Approximation3.6 Variationsformulierung der Helmholtz-Gleichung4. Lagrange-Elemente4.1 Einleitung4.2 Konstruktion von Finite-Element-Räumen4.3 Assemblierung der Galerkin-Gleichungen4.4 Beispiel: 1D Poissongleichung4.5 Konvergenz4.6 Numerische IntegrationTU Bergakademie Freiberg | INMO | O. Rheinbach/O. Ernst | Numerische Simulation mit finiten Elementen | SS 2015 29

Einleitung

VORBEMERKUNGENEINORDNUNG

Die Finite-Elemente-Methode (FEM) bezeichnet eine große Klassenumerischer Verfahren zur näherungsweisen Lösung partiellerDifferentialgleichungen (PDG, PDE).

Andere gebräuchliche Verfahrensklassen hierfür sind

finite Differenzen (Differenzenquotienten, Tensorproduktgitter),finite Volumen (Integration mittels Gauß-Quadratur),Spektralverfahren (Integraltransformation, Anpassung vonSpektralkoeffizienten),Kollokationsverfahren (Interpolations-/Ausgleichsrechnung,implizite Runge-Kutta-Verfahren),Randelementverfahren (Diskretisierung der Oberflächen, 3D –>2D).

TU Bergakademie Freiberg | INMO | O. Rheinbach/O. Ernst | Numerische Simulation mit finiten Elementen | SS 2015 30

Einleitung

VORBEMERKUNGENVORTEILE DER FEM

Geometrische Flexibilität: Die Anpassung an komplizierte Geometrienwird in die Gittererzeugung verlagert, dasgrundlegende Verfahren bleibt davon unabhängig.

Mathematisches Fundament: Es existiert eine umfassende undausgereifte mathematische Konvergenztheorie, mittelsderer Konvergenzrate, Fehlerschätzer etc. analysiertwerden können.

praktische Handhabung: Randbedingungen; Symmetrien (undweitere Strukturen) bleiben erhalten; hohe Ordnungmöglich.

Weite Verbreitung: Es existieren inzwischen sehr vieleSoftwarepakete hoher Qualität, in denen die FEMrealisiert ist, etwa MSC Nastran, ANSYS, ABAQUS,STRESS CHECK, COMSOL neben vielen anderen.


Einleitung

VORBEMERKUNGENHISTORISCHES

Johann Bernoulli (1696) Brachistochronen-Problem.

Formulierung von Differentialgleichungen alsExtremalprobleme. Beginn der Variationsrechnung.


Einleitung


Karl Schellbach (1851) Lösung eines Minimalflächenproblems in 2D mit fürFEM typischen Teilschritten.

K. Schellbach, ‘Probleme der Variationsrechnung’, J. Reine Angew.

Math. 41,293-363 (1851).


Einleitung


Walter Ritz (1908) Minimierung eines quadratischen Funktionals(Energiefunktional) in einem endlichdimensionalenFunktionenraum.

Boris Galerkin (1915) Lösung einer Randwertaufgabe aufendlichdimensionalem Funktionenraum mittelsVariationsformulierung, „Methode der gewichteten Residuen“.

Richard Courant (1943) Verwendete zum ersten Mal Ansatzfunktionen mitkleinem Träger (lineares Dreieck).


Einleitung


50er Jahre. FEM von Mechanikern neu entdeckt. Zerlegung vonFestkörpern in endlich viele „Finite Elemente“, Berechnungder Verschiebungen unter gegebenen Lasten in den Knotender Finiten Elemente.

60er Jahre. Theoretische Untermauerung der FEM seitens derMathematik. Computer-Programm NASTRAN wird von derMacNeal-Schwendler Corporation vermarktet.

1967. Ingenieur-Monographie The Finite Element Method inContinuum and Structural Mechanics von Zienkiewicz undCheung erscheint.

1973. Mathematische Monographie An Analysis of the FiniteElement Method von Strang und Fix erscheint.


Einleitung

VORBEMERKUNGENLITERATUR

O. C. Zienkiewicz und R. L. Taylor. The Finite Element Method,Volume 1 and 2. 4th ed., McGraw-Hill, New York, 1989.

T. J. R. Hughes. The Finite Element Method – Linear Static andDynamic Finite Element Analysis. Prentice-Hall,Englewood Cliffs, NJ, 1987.

M. Jung und U. Langer. Methode der finiten Elemente für Ingenieure.Teubner, Stuttgart, 2001.

P. Monk. Finite Element Methods for Maxwell’s Equations.Oxford University Press, 2003

J. Jin. The Finite Element Method in Electromagnetics. JohnWiley & Sons, 2002


Einleitung


Einleitung

RAND- UND ANFANGSWERTAUFGABENTYPEN VON PARTIELLEN DIFFERENTIALGLEICHUNGEN

Man unterscheidet zur Beschreibung der Ausbreitungunterschiedlichster physikalischer Felder im Gebiet Ω drei Typenpartieller Differentialgleichungen (PDG, PDE) 2. Ordnung:

elliptisch Potentialverfahren (z.B. Geoelektrik), Elektromagnetik(EM) im Frequenzbereich

−∇ · (k∇u) + bu = f, −∆u = f,

parabolisch Diffusion/Wärmeleitung, EM im Zeitbereich

−∆u+ a∂u

∂t= f,

hyperbolisch Wellenausbreitung (z.B. Georadar, Seismik)

−∆u+1

c2∂2u

∂t2= f.


Einleitung

RAND- UND ANFANGSWERTAUFGABENTYPEN VON RANDBEDINGUNGEN

Drei Arten von Randbedingungen (RB, BC) sind auf dem Rand δΩandwendbar:

Dirichletsche RBu = r,

Neumannsche RB∂u

∂ne~n = r,

Gemischte RB (Robin-Typ)

α∂u

∂ne~n+ βu = r.

Für r = 0 spricht man von homogenen und für r 6= 0 voninhomogenen Randbedingungen.


Einleitung

RAND- UND ANFANGSWERTAUFGABENANFANGSWERTE

Für zeitabhängige Probleme müssen ebenfalls Anfangswerte für t = 0festgelegt werden, z.B.

u0 = u(t = 0) = s,

oderu0 = u(t = 0) = sin(π · x).

Es sind auch zeitabhängige Randwerte denkbar, z.B.

u = u0 + du · sin(kπ · t).


Einleitung

RAND- UND ANFANGSWERTAUFGABENELLIPTISCHE PDE - ANWENDUNG IN DER GEOPHYSIK I

Divergenz-Gradient-Operator und Helmholtz-Gleichung

−∇ · (k∇u) + bu = f

wobei

u : Ω → R eine skalare (geo-)physikalische Größe,

k, b : Ω → R positive Koeffizientenfunktionen und

f : Ω → R den sog. Quellterm

darstellen.

Die Koeffizienten k und b können skalar oder auchpositiv-definite Tensoren (d× d Matrix) sein undbeschreiben Material- und/oder Modelleigenschaften.

Laplace-Operator und Poisson-Gleichung (k ≡ const., b ≡ 0)

−∆u = fTU Bergakademie Freiberg | INMO | O. Rheinbach/O. Ernst | Numerische Simulation mit finiten Elementen | SS 2015 41

Einleitung

RAND- UND ANFANGSWERTAUFGABENELLIPTISCHE PDE - ANWENDUNG IN DER GEOPHYSIK II

Geophysikalische Anwendungen zur elliptischen PDE:

Anwendung u k b f

Geoelektrik U σ 0 Istat. Wärmeleitung T κ 0 Q2D MT E/H µ/σ f , ε, σ/µ 0 (RB)Gravimetrie Φ 1 0 ρ

Bemerkung: In vielen Anwendungen ist der Flussvektor ∇u von größerem

Interesse als die skalare Potentialfunktion u. Ersterer wird oft durch

(numerische) Differentiation der (numerisch berechneten) Lösung u

gewonnen, ein instabiler Vorgang. Sog. gemischte FE-Formulierungen

gestatten die direkte numerische Berechnung des Flusses.


Einleitung

RAND- UND ANFANGSWERTAUFGABENPARABOLISCHE PDE - ANWENDUNG IN DER GEOPHYSIK

Diffusionsgleichung

−∆u+ a∂u

∂t= f

wobei


a : Ω → R eine positive Koeffizientenfunktion und


darstellen.

Anwendung u a f

TDEM E σ ∂j/∂t, 0 (AB)Wärmeleitung T κ Wärmequelle


Einleitung

RAND- UND ANFANGSWERTAUFGABENHYPERBOLISCHE PDE - ANWENDUNG IN DER GEOPHYSIK

Wellengleichung

−∆u+1

c2∂2u

∂t2= f

wobei


c : Ω → R eine positive Koeffizientenfunktion und


darstellen.

Anwendung u c f

Georadar (2D) E ε iωjSeismik ∆x v ∆x0


Variationsformulierung von Randwertaufgaben



VARIATIONSFORMULIERUNG VON RANDWERT-AUFGABENVARIATIONSRECHNUNG

Funktional reelle Funktion einer Funktion, z.B. Integral über eineunbekannte Funktion und deren Ableitungen

stationäre Funktionen Funktional nimmt ein Extremum an (Minimum,Maximum, Sattelpunkt)

physikalische Extremalprinzipien z.B. Lagrange-Formalismus derklassischen Mechanik, Energieansätze

Variation kleine Veränderungen um die gesuchte Lösung herum



BEZEICHNUNGEN UND MATHEMATISCHE HILFS-MITTELGEBIET, RAND

Wir betrachten nun Randwertprobleme, deren Lösung (eine odermehrere) Funktionen von d unabhängigen Variablen sind (d = 2 oderd = 3, d ... Dimensionalität). Als Gebiet Ω ⊂ Rd bezeichnen wir eineoffene und zusammenhängende Menge. Lösungen vonRandwertproblemen seien hier stets auf beschränkten Gebietendefiniert:

u : Ω → R, x 7→ u(x),

wobei wir Punkte x ∈ Ω mit

x =

[

x1

x2

]

oder x =

[

xy

]

bzw. x =

x1

x2

x3

oder x =

xyz

bezeichnen. Die Menge aller Punkte auf dem Rand von Ω bezeichnenwir mit ∂Ω oder auch mit Γ.TU Bergakademie Freiberg | INMO | O. Rheinbach/O. Ernst | Numerische Simulation mit finiten Elementen | SS 2015 48


BEZEICHNUNGEN UND MATHEMATISCHE HILFS-MITTELZULÄSSIGE GEBIETE

Die Beschaffenheit von ∂Ω beeinflusst Eigenschaften von Räumenauf Ω definierter Funktionen sowie die Regularität von Lösungen aufΩ definierter Randwertprobleme. Wir schränken daher die Menge derzulässigen Gebiete wie folgt ein: Eine Funktion f : D ⊂ Rn → R heißtLipschitz-stetig auf D, falls es eine Konstante L gibt, sodass

|f(x)− f(y)| ≤ L|x− y| ∀x,y ∈ D.

Man sagt, ein Gebiet Ω besitze einen Lipschitz-stetigen Rand, falls ∂Ωlokal durch eine Lipschitz-stetige Funktion (als Kurve im R2 bzw.Fläche im R3) parametrisiert werden kann.1 Dies schließt imWesentlichen Schlitzgebiete oder solche mit Spitzen aus. Für diePraxis genügt es oft, zu wissen, dass etwa polygonal berandeteGebiete oder beschränkte konvexe Gebiete einen Lipschitz-stetigenRand besitzen.

1Die genaue Definition findet der interessierte Leser etwa im Buch von Ciarlet.TU Bergakademie Freiberg | INMO | O. Rheinbach/O. Ernst | Numerische Simulation mit finiten Elementen | SS 2015 49


BEZEICHNUNGEN UND MATHEMATISCHE HILFS-MITTELABLEITUNGEN

Partielle Ableitungen einer Funktion u : Ω → R bezeichnen wir mit

uxi=

∂u

∂xi, uxixj

=∂2u

∂xi∂xj, i, j = 1, . . . , d, oder allgemein:

Dαu =∂|α|u

∂xα1

1· · · ∂xαd

d

, α = (α1, . . . , αd),∈ Nd0, |α| = α1 + · · ·+ αd.



BEZEICHNUNGEN UND MATHEMATISCHE HILFS-MITTELDIFFERENTIALOPERATOREN

Ferner benötigen wir die Differentialoperatoren

∇u =

ux1

...uxd

Gradient,

∇·u = (u1)x1+ · · ·+ (ud)xd

Divergenz (u : Ω → Rd),

∆u = ∇·(∇u) = ux1x1+ · · ·+ uxdxd

Laplace-Operator,

∇×u =

∣

∣

∣

∣

∣

∣

ex1. . . exd

∂∂x1

. . . ∂∂xd

u1 . . . ud

∣

∣

∣

∣

∣

∣

Rotation (u : Ω → Rd).



DIFFERENTIALOPERATORENÜBUNG

Aufgabe Gegeben seien die Funktion u : R3 → R mitu(x, y, z) = xy + xz + yz und die Vektorfelderw, v : R3 → R3 mit

v(x, y, z) =

y2

32yx

und w(x, y, z) =

y − zz − xx− y

.

Berechnen Sie

a)∇·∇u

b)∇×∇(w · v)

c)∇×∇× v!



BEZEICHNUNGEN UND MATHEMATISCHE HILFS-MITTELNORMALABLEITUNG

Ist ∂Ω Lipschitz-stetig, so ist für fast alle x ∈ ∂Ω ein Normalenvektordefiniert. Den äußeren Einheitsnormalenvektor (kurz:Normalenvektor) im Punkt x ∈ ∂Ω bezeichnen wir mit n = n(x).

Die Richtungsableitung von u längs des äußerenEinheitsnormalenvektors im Punkt x ∈ ∂Ω

∂nu :=∂u

∂n:= n · ∇u

heißt Normalableitung von u im Punkt x.



BEZEICHNUNGEN UND MATHEMATISCHE HILFS-MITTELFUNKTIONENRÄUME

Wie im Eindimensionalen sei

L2(Ω) := u : Ω → R : ‖u‖0 < ∞, ‖u‖0 =

(∫

Ω

u(x)2 dx

)1/2

.

L2(Ω) ist ein Hilbert-Raum mit Innenprodukt

(u, v) =

∫

Ω

u(x)v(x) dx, insbesondere ‖u‖0 = (u, u)1/2.

Der Raum H1(Ω) := u ∈ L2(Ω) : uxi∈ L2(Ω), i = 1, . . . , d ist

ebenfalls ein Hilbert-Raum bezüglich des Innenprodukts

(u, v)1 :=

∫

Ω

(uv +

d∑

i=1

uxivxi

) dx =

∫

Ω

(

∑

|α|≤1

DαuDαv

)

dx

mit zugehöriger Norm ‖u‖1 = (u, u)1/21

.TU Bergakademie Freiberg | INMO | O. Rheinbach/O. Ernst | Numerische Simulation mit finiten Elementen | SS 2015 54


EUKLIDISCHER VEKTORRAUMÜBUNG

Überlegung Skalarprodukt, Norm/Länge, Winkel, Orthogonalität,Orthonormalbasen in unserem Anschauungsraum,dem euklidischen Vektorraum Rn



BEZEICHNUNGEN UND MATHEMATISCHE HILFS-MITTELDIVERGENZSATZ

Das Pendant zu partieller Integration ist im Mehrdimensionalen derDivergenzsatz (Gaußscher Integralsatz, Green’s theorem). Wie imEindimensionalen ist dies das wichtigste Hilfsmittel bei der Herleitungvon Variationsformulierungen.

Satz 1

Sei Ω ⊂ Rd ein zulässiges Gebiet sowie ui ∈ H1(Ω), i = 1, . . . , d, sogilt

∫

Ω

∇·u dx =

∫

Ω

((u1)x1+ · · ·+ (ud)xd

) dx =

∫

∂Ω

n · u ds, (9)

wobei u = [u1, . . . , ud]⊤, n den Normalenvektor und ds Integration

über den Rand bezeichnen.TU Bergakademie Freiberg | INMO | O. Rheinbach/O. Ernst | Numerische Simulation mit finiten Elementen | SS 2015 56


BEZEICHNUNGEN UND MATHEMATISCHE HILFS-MITTELSCHWACHE ABLEITUNGEN

Die auftretenden Ableitungen sind im schwachen Sinne zu verstehen,im Allgemeinen sind H1-Funktionen nicht stetig-differenzierbar. Wirbeschreiben hier kurz das Prinzip schwacher Ableitungen.Besitzt die Funktion u eine stetige partielle Ableitung ux nach x, sogilt nach dem Divergenzsatz für jede differenzierbare Funktion φwelche auf ∂Ω verschwindet (setze u1 = uφ, u2 = · · · = ud = 0 in (9))

∫

Ω

uφx dx = −

∫

Ω

uxφ dx. (10)

(10) kann jedoch auch für nicht-differenzierbare Funktionen gelten:sind u und v integrierbare Funktionen mit der Eigenschaft

∫

Ω

uφx dx = −

∫

Ω

vφ dx, ∀φ, φ differenzierbar,

so bezeichnet man v als schwache Ableitung von u nach x.TU Bergakademie Freiberg | INMO | O. Rheinbach/O. Ernst | Numerische Simulation mit finiten Elementen | SS 2015 57


SCHWACHE ABLEITUNGENBEISPIEL

Die stückweise lineare Funktion

u(x) =

2x, 0 ≤ x ≤ 1

2,

2(1− x), 1

2≤ x ≤ 1

ist wegen des Knicks bei x = 1

2auf dem Intervall (0, 1) nicht differenzierbar.

Für jede differenzierbare Funktion φ mit φ(0) = φ(1) = 0 gilt jedoch∫

1

0

uφ′ dx =

∫

1/2

0

2xφ′ dx+

∫

1

1/2

2(1− x)φ′ dx

= −

(

∫

1/2

0

2φ dx+

∫

1

1/2

(−2)φ dx

)

= −

(∫

1

0

vφ dx

)

,

und damit im obigen Sinn u′ = du/dx = v mit

v(x) =

2, 0 ≤ x ≤ 1

2,

−2, 1

2≤ x ≤ 1.



SCHWACHE ABLEITUNGENWEITERES BEISPIEL

Das Gebiet Ω bestehe aus der Vereinigung dreier Dreiecke Ki, i = 1, 2, 3, mitzwei gemeinsamen Kanten. Die Funktion u : Ω → R sei auf jedem TeildreieckKi stetig differenzierbar, auf Ω jedoch nur stetig. Dann gilt nach (9)∫

Ω

uφx dx =

3∑

i=1

∫

Ki

uφx dx =

3∑

i=1

∫

∂Ki

n ·

(

uφ0

)

ds−3∑

i=1

∫

Ki

uxφ dx.

Da die Randintegrale sich im Inneren des Gebietes wegheben (warum?) folgt,falls φ auf ∂Ω verschwindet,

∫

Ω

uφx dx = −3∑

i=1

∫

Ki

uxφ dx.

Wie man sieht, stimmt die Ableitung von u im Inneren von Ki mit der

klassischen Ableitung überein, was auf den Kanten geschieht ist unerheblich.

Insbesondere können stückweise differenzierbare Funktionen auch

stückweise differenziert werden.



BEZEICHNUNGEN UND MATHEMATISCHE HILFS-MITTELSOBOLEVSCHE EINBETTUNG I

Eine Funktion mit schwachen Ableitungen genügend hoher Ordnung besitztauch klassische Stetigkeits- bzw. Differenzierbarkeitseigenschaften. Seien

Hm(Ω) = u ∈ L2(Ω) : ‖u‖m < ∞, ‖u‖m =

(∫

Ω

∑

|α|≤m

|Dαu|2 dx

)1/2

,

L∞(Ω) = u : ‖u‖∞ < ∞, ‖u‖∞ = supx∈Ω

|u(x)|.

Satz 2 (Sobolevscher Einbettungssatz)

Sei Ω ⊂ Rd ein zulässiges Gebiet. Für m > k + d/2 existiert eineKonstante C mit

‖Dαu‖∞ ≤ C‖u‖m für alle |α| ≤ k.

Ferner enthält die L∞(Ω)-Äquivalenzklasse jeder Funktionu ∈ Hm(Ω) eine stetige Funktion.TU Bergakademie Freiberg | INMO | O. Rheinbach/O. Ernst | Numerische Simulation mit finiten Elementen | SS 2015 60


BEZEICHNUNGEN UND MATHEMATISCHE HILFS-MITTELSOBOLEVSCHE EINBETTUNG II

Der Sobolevsche Einbettungssatz besagt, dass eine Hm-Funktion alsk-mal stetig differenzierbar angesehen werden kann.2

Insbesondere sind im Eindimensionalen (d = 1) die H1-Funktionen(k = 0) stetig, d.h. deren punktweise (diskrete) Auswertung istsinnvoll.

2Genauer: in der L∞-Äquivalenzklasse von u liegt eine k-mal stetig differenzierbareFunktion.



BEZEICHNUNGEN UND MATHEMATISCHE HILFS-MITTELRANDWERTE

Die punktweise Auswertung von Funktionen auf ∂Ω – notwendig für dieFormulierung von Randwertaufgaben – ist auch mit Hilfe der SobolevschenEinbettungen im Allgemeinen nicht möglich.

Es ist jedoch möglich, die Randwerte etwa von Hm-Funktionen und derenAbleitungen als Funktionen in L2(∂Ω) aufzufassen. Für Randwertaufgabenzweiter Ordnung benötigen wir die Randwerte u und ∂nu, diese liegen füru ∈ H2(Ω) in L2(∂Ω).3

Die Zuordnung von Funktionen in Hm(Ω) zu Randwerten in L2(∂Ω)bezeichnet man als Spuroperator, Aussagen über die Regularität vonRandwerten als Spursätze.

Im Sinne dieser Spursätze verstehen wir im Folgenden stillschweigend die

Randwerte von Hm-Funktionen.3Ist u nur in H1(Ω), so existiert ∂nu in einem schwächeren Sinne.



BEZEICHNUNGEN UND MATHEMATISCHE HILFS-MITTELSOBOLEV-RÄUME

Folgende Sobolev-Räume von Funktionen stellen den natürlichenAnsatzraum für die Variationsformulierung von Randwertaufgabenzweiter bzw. vierter Ordnung mit homogenen Randbedingungen dar:

H1

0 (Ω) := u ∈ H1(Ω) : u = 0 auf ∂Ω, (11)

H2

0 (Ω) := u ∈ H2(Ω) : u = ∂nu = 0 auf ∂Ω. (12)

Die homogenen Randbedingungen können sich auch nur über Teiledes Randes erstrecken.



REFERENZAUFGABEN IN 2DBEISPIEL 1

Beispiel 1: Ω = (−1, 1)2, ΓD = Γ, f ≡ 1, gD ≡ 0.Dies stellt ein Modell für Wärmeausbreitung auf einer quadratischen Plattedar. Durch Trennung der Veränderlichen bestimmt man die Reihenlösung

u(x, y) =1− x2

2−

16

π3

∑

k∈N,k ungerade

[

sin(kπ(1 + x)/2)

k3 sinh(kπ)

(

sinhkπ(1 + y)

2+ sinh

kπ(1− y)

2

)]

.

−1

−0.5

0

0.5

1

−1

−0.5

0

0.5

1−0.05

0

0.05

0.1

0.15

0.2

0.25

0.3

−1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1−1

−0.8

−0.6

−0.4

−0.2

0

0.2

0.4

0.6

0.8

1




Beispiel 2: Ω = (−1, 1)2 \ [−1, 0]2, ΓD = Γ, f ≡ 1, gD ≡ 0.

−1

−0.5

0

0.5

1

−1

−0.5

0

0.5

10

0.02

0.04

0.06

0.08

0.1

0.12

0.14

0.16

−1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1−1

−0.8

−0.6

−0.4

−0.2

0

0.2

0.4

0.6

0.8

1




Beispiel 3: Ω = (−1, 1)2, ΓD = Γ, f ≡ 0, gD = u|Γmit exakter Lösung

u(x, y) =2(1 + y)

(3 + x)2 + (1 + y)2.

−1

−0.5

0

0.5

1

−1

−0.5

0

0.5

10

0.1

0.2

0.3

0.4

0.5

−1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1−1

−0.8

−0.6

−0.4

−0.2

0

0.2

0.4

0.6

0.8

1




Beispiel 4: Ω = (−1, 1)2 \ [−1, 0]2, ΓD = Γ, f ≡ 1, gD = u|Γ mitexakter Lösung

u(r, θ) = r2/3 sin2θ + π

3, x = r cos θ, y = r sin θ.

−1

−0.5

0

0.5

1

−1

−0.5

0

0.5

10

0.2

0.4

0.6

0.8

1

1.2

1.4

−1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1−1

−0.8

−0.6

−0.4

−0.2

0

0.2

0.4

0.6

0.8

1



VARIATIONSFORMULIERUNG DER POISSON-GLEICHUNGELLIPTISCHE PARTIELLE DIFFERENTIALGLEICHUNG ZWEITER ORDNUNG

Viele physikalische Größen erfüllen eine elliptischeDifferentialgleichung zweiter Ordnung der Form

−∇·(k∇u) + bu = f, u : Ω → R, (13)

wobei

k, b : Ω → R positive Koeffizientenfunktionen und

f : Ω → R einen sog. Quellterm

darstellen. Die Koeffizienten k und b können skalar oder auchpositiv-definite Tensoren (d× d Matrix) sein und beschreiben meistMaterial- und Modelleigenschaften.

Die typische Problemstellung besteht darin, u zu gegebenen f , k undb zu bestimmen. (Hinzu kommen noch Randbedingungen auf ∂Ω.)TU Bergakademie Freiberg | INMO | O. Rheinbach/O. Ernst | Numerische Simulation mit finiten Elementen | SS 2015 70


VARIATIONSFORMULIERUNG DER POISSON-GLEICHUNGDIE POISSON-GLEICHUNG

Ist k konstant in (13), so kann dies in f zusammengefasst werden undes ergibt sich die Poisson-Gleichung

−∆u(x) = f(x), x ∈ Ω, (14a)

wobei Ω ⊂ Rd (d = 2, 3) ein zulässiges Gebiet sei. Den GebietsrandΓ := ∂Ω zerlegen wir in ΓD und ΓN , Γ = ΓD ∪ ΓN , ΓD ∩ ΓN = ∅ undstellen die Randbedingungen

u(x) = gD(x) ∀x ∈ ΓD, kurz: u|ΓD= gD, (14b)

∂u

∂n(x) = gN (x) ∀x ∈ ΓN , kurz: ∂nu|ΓN

= gN (14c)

mit zwei gegebenen, auf ΓD bzw. ΓN definierten Funktionen g und h.Die RWA (14) besitzt – unter geeigneten Voraussetzungen an dasGebiet Ω und die Daten f, gD, gN – eine eindeutig bestimmteklassische (d.h. in Ω zweimal stetig differenzierbare) Lösung.TU Bergakademie Freiberg | INMO | O. Rheinbach/O. Ernst | Numerische Simulation mit finiten Elementen | SS 2015 71


VARIATIONSFORMULIERUNG DER POISSON-GLEICHUNGL2-INNENPRODUKT UND DIVERGENZSATZ

Wir multiplizieren (14a) mit einer Testfunktion v und wenden denDivergenzsatz an:

(f, v) = −

∫

Ω

v∆u dx = −

∫

Ω

(

∇·(v∇u)−∇u · ∇ v)

dx

=

∫

Ω

∇u · ∇ v dx−

∫

Γ

v∂u

∂nds = (∇u,∇ v)− (∂nu, v)Γ,

wobei (·, ·) hier auch das L2-Innenprodukt vektorwertiger Funktionenauf Ω sowie (·, ·)Γ das L2-Innenprodukt auf Γ bezeichnen mögen.

Wir wählen nun die Testfunktion v so, dass diese auf ΓD

verschwindet. Auf dem Neumann-Rand ΓN gilt nach (14c) ∂nu = gN .Insgesamt ergibt sich

(∇u,∇ v) = (f, v) + (gN , v)ΓN. (15)



VARIATIONSFORMULIERUNG DER POISSON-GLEICHUNGFUNKTIONENRÄUME

Wir setzen nun

a(u, v) := (∇u,∇ v) =

∫

Ω

∇u · ∇ v dx,

ℓ(v) := (f, v) + (gN , v)ΓN=

∫

Ω

fv dx+

∫

ΓN

gNv ds.

Damit alle Integrale definiert sind, reicht es aus, dass die erstenAbleitungen von u und v quadratisch integrierbar sind, wir könnenalso u, v ∈ H1(Ω) wählen und setzen

S = u ∈ H1(Ω) : u|ΓD= gD, V = v ∈ H1(Ω) : v|ΓD

= 0.

Die Variationsformulierung von (14) lautet somit

Bestimme u ∈ S so, dass a(u, v) = ℓ(v) für alle v ∈ V . (16)



VARIATIONSFORMULIERUNG DER POISSON-GLEICHUNGBILINEARFORMEN I

Definition 3

Sei V ein reeller normierter Vektorraum. Eine Bilinearform ist eineAbbildung

a : V × V → R,

welche linear in beiden Argumenten ist, d.h. es gilt

a(u1 + u2, v) = a(u1, v) + a(u2, v)

a(u, v1 + v2) = a(u, v1) + a(u, v2)

a(λu, v) = λa(u, v)

a(u, vλ) = a(u, v)λ

∀u, u1, u2, v, v1, v2 ∈ V und λ ∈ R.



VARIATIONSFORMULIERUNG DER POISSON-GLEICHUNGBILINEARFORMEN II

Definition 4

Die Bilinearform a heißt stetig, falls es eine Konstante C gibt mit

|a(u, v)| ≤ C‖u‖‖v‖ ∀u, v ∈ V .

Die Bilinearform a heißt symmetrisch, falls

a(u, v) = a(v, u) ∀u, v ∈ V .

Die Bilinearform a heißt koerziv, falls es eine Konstante α > 0 gibt mit

a(u, u) ≥ α‖u‖2 ∀u ∈ V .



VARIATIONSFORMULIERUNG DER POISSON-GLEICHUNGBILINEARFORMEN III

Bemerkungen 5

(a) Im Falle eines komplexen Vektorraumes fordert man Antilinearitätim zweiten Argument und spricht dann von einerSesquilinearform. Anstelle der Symmetrie fordert man hiera(u, v) = a(v, u) und spricht von einer HermiteschenSesquilinearform.

(b) Eine koerzive symmetrische (Hermitesche) Bilinearform(Sesquilinearform) definiert ein Innenprodukt auf dem reellen(komplexen) Vektorraum V . Oft wird es das zur Bilinearforma(·, ·) gehörende Energie-Innenprodukt genannt.

(c) Die Eigenschaft der Stetigkeit zieht die Stetigkeit in beidenArgumenten nach sich.



VARIATIONSFORMULIERUNG DER POISSON-GLEICHUNGEIN ABSTRAKTER EXISTENZSATZ

Der folgende Satz sichert Existenz und Eindeutigkeit der Lösung einergroßen Klasse von Variationsproblemen:

Satz 6 (Lax-Milgram-Lemma, 1954)

Sei V ein Hilbert-Raum mit Norm ‖ · ‖V , a : V × V → R eineBilinearform auf V sowie ℓ : V → R ein lineares Funktional auf V für diees Konstanten C,α und L gibt mit

|a(u, v)| ≤ C‖u‖V ‖v‖V ∀u, v ∈ V , („ a ist stetig “)

a(v, v) ≥ α‖v‖2V ∀v ∈ V , („ a ist koerziv “)

|ℓ(v)| ≤ L‖v‖V ∀v ∈ V , („ ℓ ist stetig “).

Dann besitzt das Variationsproblem

Bestimme u ∈ V sodass a(u, v) = ℓ(v) ∀v ∈ V

genau eine Lösung.



VARIATIONSFORMULIERUNG DER POISSON-GLEICHUNGVARIATIONS- UND MINIMIERUNGSAUFGABEN I

Anstelle von (16) könnte zunächst folgende Variante einer schwachenFormulierung naheliegender erscheinen: Im einfachsten Fall der reinenDirichletsche Randwertaufgabe (RWA) für die Poisson-Gleichung(Γ = ΓD)

−∆u = f auf Ω, (17a)

u = g auf Γ = ∂Ω, (17b)

betrachten wir die sog. verallgemeinerte Randwertaufgabe

Bestimme u ∈ C2(Ω) mit u = g längs Γ und∫

Ω

∇u · ∇ v dx =

∫

Ω

fv dx ∀v ∈ C∞0 (Ω).

(18)

Schließlich betrachten wir noch die Minimierungsaufgabe

Minimiere unter allen Funktionen u ∈ C2(Ω) mit u = g längs Γ

das Funktional J(u) :=1

2

∫

Ω

| ∇u|2 dx−

∫

Ω

fu dx.(19)



VARIATIONSFORMULIERUNG DER POISSON-GLEICHUNGVARIATIONS- UND MINIMIERUNGSAUFGABEN II

Es gilt nun:

Satz 7

Seien g ∈ C(Γ) sowie f ∈ C(Ω) gegebene Funktionen. Sei ferneru ∈ C2(Ω). Dann gilt

(a) Löst u die Minimierungsaufgabe (19), so löst u auch dieRandwertaufgabe (17).

(b) Die Funktion u ist genau dann Lösung der RWA (17), wenn u Lösungder verallgemeinerten RWA (18) ist.

Bemerkung 8

Nach Satz 7 löst jede hinreichend glatte Lösung der Variationsaufgaben(18) oder (19) auch die RWA (17).Entscheidend: in vielen Anwendungen tritt der Fall auf, dass keine derartglatte Lösung existiert.



VARIATIONSFORMULIERUNG DER POISSON-GLEICHUNGVOLLSTÄNDIGKEIT DER FUNKTIONENRÄUME

Die Situation entspricht der bei den Minimierungsaufgaben

f(x) −→min x ∈ [a, b] (20a)

und f(x) −→min x ∈ [a, b] ∩Q, (20b)

wobei −∞ < a < b < ∞ und f ∈ C[a, b].Aufgabe (20a) besitzt stets Lösungen. Sind diese jedoch allesamtirrationale Zahlen, so besitzt Aufgabe (20b) keine Lösung.

Auf analoge Weise muss im Fall der Variationsaufgaben i.A. derFunktionenraum geeignet vervollständigt werden. Dies führt hier aufdie sog. Sobolev-Räume.



VARIATIONSFORMULIERUNG DER POISSON-GLEICHUNGHOMOGENISIERUNG WESENTLICHER RANDBEDINGUNGEN

Oft ist es praktischer, mit homogenen wesentlichen Randbedingungenzu arbeiten.

Insbesondere kann man dann in der Variationsformulierung (16)denselben Funtionenraum für die Ansatz- und Testfunktionen wählen,d.h. S = V .Sei hierzu ug ∈ H1(Ω) eine auf ganz Ω definierte Funktion mit derEigenschaft

ug = g auf ΓD.

(Für zulässige Gebiete existiert eine solche Fortsetzung von g nachinnen.) Dann liegt aber für jede Funktion u ∈ V die Summe ug + u inS und es gilt

a(u, v) = ℓ(v)− a(ug, v) ∀v ∈ V .



BEISPIELWÄRMELEITUNG IN EINEM KÖRPERQUERSCHNITT I

Wärmeleitungsgleichung in Ω1 ∪ Ω2 ⊂ R2: −∇·(k∇u) = 0,zwei verschiedene Materialien mit

k(x) =

k1 = 1W(mK)−1, x ∈ Ω1,

k2 = 371W(mK)−1, x ∈ Ω2.

Randbedingungen:

u = g1 auf Γ1,

∂nu = 0 auf Γ2,

k∂nu+ α(u− g3) = 0 auf Γ3,

g1 ≡ 500K,

g3 ≡ 300K, α = 5.6 W/(m2K).0 0.17 0.35 0.65 0.83 1

0

0.15

0.3

0.45

0.6

0.8

Ω1

Ω2

Γ1

Γ3

Γ2

Γsym



BEISPIELWÄRMELEITUNG IN EINEM KÖRPERQUERSCHNITT II

Entlang des gemeinsamen Randes von Ω1 und Ω2 sind Temperaturund Wärmefluss stetig, d.h. dort gilt

u1(x) = u2(x), ∂n[k1u1(x)] = ∂n[k2u2(x)].

Zur Variationsformulierung wählen wir den Ansatzraum

S = u ∈ H1(Ω) : u = g1 auf Γ1,

sowie den Testraum

V = u ∈ H1(Ω) : u = 0 auf Γ1,

und suchen u ∈ S sodass a(u, v) = ℓ(v) ∀v ∈ V , mit

a(u, v) =

∫

Ω

∇ v · (k∇u) dx+ α

∫

Γ3

uv ds,

ℓ(v) = α

∫

Γ3

vg3 ds.



BEISPIELSYMMETRIE-RANDBEDINGUNG

Da die RWA symmetrisch zur Achse x = 0.5 ist reicht es aus, dieLösung u nur auf einer Hälfte des Gebiets zu bestimmen und in derverbleibenden Hälfte die Lösung aus der Beziehungu(0.5 + x, y) = u(0.5− x, y) zu bestimmen.

Längs der Symmetrieachse erfüllt die Lösung eine homogeneNeumann-Randbedingung:

∂nu(x) = 0, x ∈ Γsym.



GALERKIN-APPROXIMATION

Gegeben sei die Variationsaufgabe (ohne Einschränkung: S = V )Bestimme u ∈ V sodass

a(u, v) = ℓ(v) ∀v ∈ V . (21)

Bei der Galerkin-Approximation der Lösung u von (21) erstetzt manden Variationsraum durch einen endlichdimensionalen UnterraumV h ⊂ V und betrachtet die Lösung uh ∈ V h von (21) mit V h anstellevon V als Approximation an u.Mit anderen Worten: wir bestimmen uh ∈ V h sodass

a(uh, v) = ℓ(v) ∀v ∈ Vh.

Bemerkung 9

Im o.g. Fall V h ⊂ V spricht man von einer konformenGalerkin-Diskretisierung. Es werden jedoch auch nichtkonformeDiskretisierungen mit endlich-dimensionalen Unterräumen V h mitV h 6⊂ V verwandt.TU Bergakademie Freiberg | INMO | O. Rheinbach/O. Ernst | Numerische Simulation mit finiten Elementen | SS 2015 86


CÉA-LEMMAEINE GRUNDLEGENDE FEHLERABSCHÄTZUNG

Satz 10 (Lemma von Céa)

Gelten für die Variationsaufgabe (21) die Voraussetzungen desLax-Milgram-Lemmas, so gilt für den Fehler u− uh derGalerkin-Approximation

‖u− uh‖V ≤C

αinf

v∈V h‖u− v‖V . (22)

Hierbei bezeichnen C und α die Stetigkeits- bzw.Koerzivitätskonstante aus dem Lax-Milgram-Lemma.


Numerische Simulation mit finiten Elementen [-.2em] - O ... · Variationsformulierung, „Methode...

Documents

Transcript of Numerische Simulation mit finiten Elementen [-.2em] - O ... · Variationsformulierung, „Methode...