Numerische Simulation mit finiten Elementen [-.2em...

Institut für Numerische Mathematik und Optimierung

Numerische Simulation mitfiniten Elementen

Antje Franke-Börner

Vorlesung im gleichnamigen ModulHörerkreis: 2. MNC, 2. MGPHY, 4. BGIP, 6. BEC-II, 2. MGINSommersemester 2013

Inhalt I

1. Einleitung1.1 Vorbemerkungen1.2 Rand- und Anfangswertaufgaben2. Variationsformulierung von Randwertaufgaben2.1 Vorbemerkungen2.2 Bezeichnungen und mathematische Hilfsmittel2.3 Referenzaufgaben in 2D2.4 Variationsformulierung der Poisson-Gleichung2.5 Galerkin-Approximation2.6 Variationsformulierung der Helmholtz-Gleichung3. Lagrange-Elemente3.1 Einleitung3.2 Konstruktion von Finite-Element-Räumen3.3 Assemblierung der Galerkin-Gleichungen3.4 Beispiel: 1D Poissongleichung

Antje Franke-Börner (TU Freiberg) Numerische Simulation mit finiten Elementen Sommersemester 2013 1

Inhalt II

3.5 Konvergenz3.6 Numerische Integration3.7 Beispiel: 2D Poissongleichung

Schema zur Assemblierung des Galerkin-GleichungssystemsBeispiel zur 2D Poissongleichung

4. Nédélec-Elemente4.1 Einleitung4.2 Neue Variationsräume

VariationsformulierungHelmholtz-Zerlegung

4.3 Divergenz- und Rotationskonforme ElementeGebietstransformationen - GradientoperatorRotationskonforme Finite Elemente


Inhalt

1. Einleitung1.1 Vorbemerkungen1.2 Rand- und Anfangswertaufgaben2. Variationsformulierung von Randwertaufgaben2.1 Vorbemerkungen2.2 Bezeichnungen und mathematische Hilfsmittel2.3 Referenzaufgaben in 2D2.4 Variationsformulierung der Poisson-Gleichung2.5 Galerkin-Approximation2.6 Variationsformulierung der Helmholtz-Gleichung3. Lagrange-Elemente3.1 Einleitung3.2 Konstruktion von Finite-Element-Räumen3.3 Assemblierung der Galerkin-Gleichungen3.4 Beispiel: 1D Poissongleichung3.5 Konvergenz3.6 Numerische Integration3.7 Beispiel: 2D Poissongleichung






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 sindfinite 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).


VorbemerkungenVorteile der FEM

Geometrische Flexibilität: Die Anpassung an komplizierte Geometrienwird in die Gittererzeugung verlagert, das grundlegendeVerfahren 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 (und weitereStrukturen) bleiben erhalten; hohe Ordnung möglich.

Weite Verbreitung: Es existieren inzwischen sehr viele Softwarepaketehoher Qualität, in denen die FEM realisiert ist, etwa MSCNastran, ANSYS, ABAQUS, STRESS CHECK, COMSOLneben vielen anderen.


VorbemerkungenHistorisches

Johann Bernoulli (1696) Brachistochronen-Problem.

Formulierung von Differentialgleichungen als Extremalprobleme.Beginn der Variationsrechnung.



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

932 SHORT COMMUNICATIONS

these, from the finite element point of view (in particular because it treats a problem in two independent variables, which Euler, Leibniz and others did not), was given by Karl Schellbach in 1851.13 The present author knows of no other instance wherein there is so early and pointed a demonstration of as many ingredients of the present-day finite element method. Schellbach’s presentation appears in a paper that is largely devoted to numerous examples in which differential equations are derived from extremum problems. One group of these examples

describes a finite element-like procedure applied to several problems in one and two independent variables, and the most interesting of these applications is Schellbach’s approach to Plateau’s problem.

He wishes to derive the partial differential equation describing the surface z(x, y ) of minimum area enclosed by a given closed curve in space. This is done through the approximation of the surface whose area is to be minimized by a piecewise linear trial function constructed with triangular elements. Instead of summing the entire collection of plane triangular surfaces in a particular approximation, Schellbach realizes that the sum of the areas for all the surfaces is

minimized when the sum of only half of them reaches a minimum. For example, in Figure 2, of

Figure 2. Schellbach’s triangular elements

the six triangles containing (xl, yl) , the areas of only the triangular surfaces associated with the

three shaded elements enter into the summation. That this is plausible is suggested immediately by the fact that the seven nodal values z(xi , yi) which determine the configurations of the three surfaces above the shaded triangles also completely determine the other three triangular surfaces which contain z(x1, y l ) . t If the sum begins with the surfaces above the three shaded triangles, and if zij = z (xi, yi), the sum is

s =4((x1-xo)2(Yl -Yo)2+(x1 -xo)2(z11 -zlo)2+(Y1 -Yo)2(z10-zoo~2~1’2

+&1 -xoY(y2 - Y d 2 + (Xl -x0)2(z12 - Z 1 d 2 + (Y2 - Yd2(Zll - zo1)2)1/2

+ 5 ~ ~ ~ 2 - x l ~ z ~ y z - y ~ ~ 2 + ~ ~ z - ~ l ~ 2 ~ ~ z z - ~ 2 1 ~ 2 + ~ Y 2 - Y 1 ~ 2 ~ z 2 1 - 2 1 1 ~ 2 ~ 1 ~ 2 + . . . .

S is differentiated with respect to 2 1 1 , or with respect to any other zii which lies at the centre of a similar hexagon, and the result is set to zero giving a difference equation which, as the mesh is

t Of course, this simplification is not applicable in most extremum problems.

K. Schellbach, ‘Probleme der Variationsrechnung’, J. Reine Angew.Math. 41,293-363 (1851).



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).



50er Jahre. FEM von Mechanikern neu entdeckt. Zerlegung von Festkörpernin endlich viele „Finite Elemente“, Berechnung derVerschiebungen unter gegebenen Lasten in den Knoten derFiniten Elemente.

60er Jahre. Theoretische Untermauerung der FEM seitens der Mathematik.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 Finite ElementMethod von Strang und Fix erscheint.


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. OxfordUniversity Press, 2003

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


Inhalt







Rand- und AnfangswertaufgabenTypen von partiellen Differentialgleichungen

Man unterscheidet zur Beschreibung der Ausbreitung unterschiedlichsterphysikalischer Felder im Gebiet Ω drei Typen partiellerDifferentialgleichungen (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.


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 von inhomogenenRandbedingungen.


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).


Rand- und AnfangswertaufgabenElliptische PDE - Anwendung in der Geophysik I

Divergenz-Gradient-Operator und Helmholtz-Gleichung−∇ · (k∇u) + bu = f

wobeiu : Ω→ R eine skalare (geo-)physikalische Größe,k, b : Ω→ R positive Koeffizientenfunktionen undf : Ω→ 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 = f


Rand- und AnfangswertaufgabenElliptische PDE - Anwendung in der Geophysik II

Geophysikalische Anwendungen zur elliptischen PDE:

Anwendung u k b fGeoelektrik U σ 0 I

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

Bemerkung: In vielen Anwendungen ist der Flussvektor ∇u von größeremInteresse 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 direktenumerische Berechnung des Flusses.


Rand- und AnfangswertaufgabenParabolische PDE - Anwendung in der Geophysik

Diffusionsgleichung−∆u+ a

∂u

∂t= f

wobei

u : Ω→ R eine skalare (geo-)physikalische Größe,a : Ω→ R eine positive Koeffizientenfunktion undf : Ω→ R den sog. Quellterm

darstellen.

Anwendung u a fTDEM E σ ∂j/∂t, 0 (AB)Wärmeleitung T κ Wärmequelle


Rand- und AnfangswertaufgabenHyperbolische PDE - Anwendung in der Geophysik

Wellengleichung

−∆u+1

c2

∂2u

∂t2= f

wobei

u : Ω→ R eine skalare (geo-)physikalische Größe,c : Ω→ R eine positive Koeffizientenfunktion undf : Ω→ R den sog. Quellterm

darstellen.

Anwendung u c fGeoradar (2D) E ε iωjSeismik ∆x v ∆x0


Inhalt







Variationsformulierung von RandwertaufgabenVariationsrechnung

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


Inhalt







Bezeichnungen und mathematische HilfsmittelGebiet, 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 Gebieten definiert:

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 Γ.


Bezeichnungen und mathematische HilfsmittelZulässige Gebiete

Die Beschaffenheit von ∂Ω beeinflusst Eigenschaften von Räumen auf Ω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ächeim R3) parametrisiert werden kann.1 Dies schließt im WesentlichenSchlitzgebiete oder solche mit Spitzen aus. Für die Praxis genügt es oft,zu wissen, dass etwa polygonal berandete Gebiete oder beschränktekonvexe Gebiete einen Lipschitz-stetigen Rand besitzen.

1Die genaue Definition findet der interessierte Leser etwa im Buch von Ciarlet.


Bezeichnungen und mathematische HilfsmittelAbleitungen

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α11 · · · ∂x

αdd

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


Bezeichnungen und mathematische HilfsmittelDifferentialoperatoren

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)∇·∇ub)∇×∇(w · v)

c)∇×∇× v!


Bezeichnungen und mathematische HilfsmittelNormalableitung

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 HilfsmittelFunktionenrä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 istebenfalls 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 .


Euklidischer VektorraumÜbung

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


Bezeichnungen und mathematische HilfsmittelDivergenzsatz

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 1Sei Ω ⊂ Rd ein zulässiges Gebiet sowie ui ∈ H1(Ω), i = 1, . . . , d, so gilt∫

Ω∇·u dx =

∫Ω

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

∫∂Ωn · u ds, (1)

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

über den Rand bezeichnen.


Bezeichnungen und mathematische HilfsmittelSchwache Ableitungen

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

Ωuφx dx = −

∫Ωuxφdx. (2)

(2) kann jedoch auch für nicht-differenzierbare Funktionen gelten: sindu und v integrierbare Funktionen mit der Eigenschaft∫

Ωuφx dx = −

∫Ωvφ dx, ∀φ, φ differenzierbar,

so bezeichnet man v als schwache Ableitung von u nach x.Antje Franke-Börner (TU Freiberg) Numerische Simulation mit finiten Elementen Sommersemester 2013 31

Schwache AbleitungenBeispiel

Die stückweise lineare Funktion

u(x) =

2x, 0 ≤ x ≤ 1

2 ,

2(1− x), 12 ≤ x ≤ 1

ist wegen des Knicks bei x = 12 auf 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, 12 ≤ 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 (1)∫

Ω

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 derklassischen Ableitung überein, was auf den Kanten geschieht ist unerheblich.Insbesondere können stückweise differenzierbare Funktionen auch stückweisedifferenziert werden.


Bezeichnungen und mathematische HilfsmittelSobolevsche 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 Funktion u ∈ Hm(Ω)eine stetige Funktion.


Bezeichnungen und mathematische HilfsmittelSobolevsche Einbettung II

Der Sobolevsche Einbettungssatz besagt, dass eine Hm-Funktion alsk-mal stetig differenzierbar angesehen werden kann.2Insbesondere sind im Eindimensionalen (d = 1) die H1-Funktionen(k = 0) stetig, d.h. deren punktweise (diskrete) Auswertung ist sinnvoll.

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


Bezeichnungen und mathematische HilfsmittelRandwerte

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 dieRandwerte von Hm-Funktionen.

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


Bezeichnungen und mathematische HilfsmittelSobolev-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:

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

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

Die homogenen Randbedingungen können sich auch nur über Teile desRandes erstrecken.


Inhalt







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 Platte dar.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π)

(sinh

kπ(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


Inhalt







Variationsformulierung der Poisson-GleichungElliptische partielle Differentialgleichung zweiter Ordnung

Viele physikalische Größen erfüllen eine elliptische Differentialgleichungzweiter Ordnung der Form

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

wobei

k, b : Ω→ R positive Koeffizientenfunktionen undf : Ω→ 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 und bzu bestimmen. (Hinzu kommen noch Randbedingungen auf ∂Ω.)


Variationsformulierung der Poisson-GleichungDie Poisson-Gleichung

Ist k konstant in (5), so kann dies in f zusammengefasst werden und esergibt sich die Poisson-Gleichung

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

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, (6b)∂u

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

mit zwei gegebenen, auf ΓD bzw. ΓN definierten Funktionen g und h.Die RWA (6) besitzt – unter geeigneten Voraussetzungen an das GebietΩ und die Daten f, gD, gN – eine eindeutig bestimmte klassische (d.h.in Ω zweimal stetig differenzierbare) Lösung.


Variationsformulierung der Poisson-GleichungL2-Innenprodukt und Divergenzsatz

Wir multiplizieren (6a) 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 (6c) ∂nu = gN . Insgesamt ergibtsich

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


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önnen alsou, v ∈ H1(Ω) wählen und setzen

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

Die Variationsformulierung von (6) lautet somit

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


Variationsformulierung der Poisson-GleichungBilinearformen I

Definition 3Sei V ein reeller normierter Vektorraum. Eine Bilinearform ist eineAbbildung

a : V × V → R,

welche linear in beiden Argumenten ist, d.h. es gilta(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 4Die 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 einer Sesquilinearform.Anstelle der Symmetrie fordert man hier a(u, v) = a(v, u) undspricht von einer Hermiteschen Sesquilinearform.

(b) Eine koerzive symmetrische (Hermitesche) Bilinearform(Sesquilinearform) definiert ein Innenprodukt auf dem reellen(komplexen) Vektorraum V . Oft wird es das zur Bilinearform a(·, ·)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 eine Bilinearformauf V sowie ` : V → R ein lineares Funktional auf V für die es KonstantenC,α 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 (8) 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 Ω, (9a)u = g auf Γ = ∂Ω, (9b)

betrachten wir die sog. verallgemeinerte Randwertaufgabe

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

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

(10)

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.

(11)


Variationsformulierung der Poisson-GleichungVariations- und Minimierungsaufgaben II

Es gilt nun:

Satz 7Seien g ∈ C(Γ) sowie f ∈ C(Ω) gegebene Funktionen. Sei ferneru ∈ C2(Ω). Dann gilt(a) Löst u die Minimierungsaufgabe (11), so löst u auch die

Randwertaufgabe (9).(b) Die Funktion u ist genau dann Lösung der RWA (9), wenn u Lösung

der verallgemeinerten RWA (10) ist.

Bemerkung 8Nach Satz 7 löst jede hinreichend glatte Lösung der Variationsaufgaben (10)oder (11) auch die RWA (9).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] (12a)und f(x) −→min x ∈ [a, b] ∩Q, (12b)

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

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


Variationsformulierung der Poisson-GleichungHomogenisierung wesentlicher Randbedingungen

Oft ist es praktischer, mit homogenen wesentlichen Randbedingungen zuarbeiten.

Insbesondere kann man dann in der Variationsformulierung (8)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 Temperatur undWärmefluss stetig, d.h. dort gilt

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

Zur Variationsformulierung wählen wir den AnsatzraumS = u ∈ H1(Ω) : u = g1 auf Γ1,

sowie den TestraumV = 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.


Inhalt







Galerkin-Approximation

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

a(u, v) = `(v) ∀v ∈ V . (13)Bei der Galerkin-Approximation der Lösung u von (13) erstetzt man denVariationsraum durch einen endlichdimensionalen Unterraum V h ⊂ Vund betrachtet die Lösung uh ∈ V h von (13) mit V h anstelle von V alsApproximation an u.Mit anderen Worten: wir bestimmen uh ∈ V h sodass

a(uh, v) = `(v) ∀v ∈ V h.

Bemerkung 9Im 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.


Céa-LemmaEine grundlegende Fehlerabschätzung

Satz 10 (Lemma von Céa)Gelten für die Variationsaufgabe (13) die Voraussetzungen desLax-Milgram-Lemmas, so gilt für den Fehler u− uh derGalerkin-Approximation

‖u− uh‖V ≤C

αinfv∈V h

‖u− v‖V . (14)

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


Céa-LemmaErgänzungen

(a) Ist die Bilinearform zusätzlich noch symmetrisch (Hermitesch), solässt sich (14) verbessern zu

‖u− uh‖ ≤√C

αinfv∈V h

‖u− v‖. (15)

(b) Für die Konvergenz einer Folge von Galerkin-Approximationen (imGrenzwert h→ 0) ist somit hinreichend, dass für die zugehörigeFamilie V hh>0 von Unterräumen gilt

limh→0

infv∈V h

‖u− v‖ = 0 ∀u ∈ V .

(c) In diesem Sinne wird durch das Céa-Lemma die Abschätzung desGalerkin-Fehlers zu einem Approximationsproblem.

(d) Die Tatsache, dass a(u− uh, v) = 0 für alle v ∈ V h wird auchGalerkin-Orthogonalität genannt.


Inhalt







Variationsformulierung der Helmholtz-GleichungMagnetotellurik

Die Magnetotellurik (MT) ist eine Methode im Bereich derGeo-Elektromagnetik, welche die natürlichen elektromagnetischen Felderals Quellen nutzt, um aus der Registrierung von induzierten elektrischenund magnetischen Feldern E bzw. H Aussagen über dieUntergrundstruktur hinsichtlich der Verteilung der elektrischenLeitfähigkeit (magnetischen Permeabilität) abzuleiten.

Anwendung

Frequenzbereich: f = 10−4 . . . 104 HzUntersuchung des Aufbaus der Erde von etwa 100 m bis 100 km(Kruste, oberer Mantel)Exploration von Rohstoffen (Erze), Kohlenwasserstoffen,Grundwasser, GeothermieInterpretationsgrößen: Übertragungsfunktionen, z.B. Z = EH−1 ...Impedanz, Hz = T(HxHy)

−1, T ... Tipper


Variationsformulierung der Helmholtz-GleichungMaxwell-Gleichungen

Ausgangspunkt sind die vier Maxwell-Gleichungen (mit den üblichenBezeichnungen)

∇×H = J + D (16a)∇×E = −B (16b)∇·D = ρ (16c)∇·B = 0. (16d)

Für die Geo-Elektromagnetik sind typischerweise folgende Annahmengerechtfertigt:

B = µH, µ konstant, (lineares Medium, konstante Permeabilität)D ≈ 0, (Vernachlässigung der Verschiebungsströme)J = σE + Je. (Ohmsches Gesetz)

Bei der MT ist Je = 0, die Anregung erfolgt durch die RBen.Antje Franke-Börner (TU Freiberg) Numerische Simulation mit finiten Elementen Sommersemester 2013 65

Variationsformulierung der Helmholtz-GleichungHarmonische Zeitabhängigkeit, Gleichungen zweiter Ordnung

Für harmonische Zeitabhängigkeit der Form eiωt gehen das Faradaysche undAmpèresche Gesetz (16b) bzw. (16a) über in

∇×E = −iωµH, (17a)∇×H = σE. (17b)

Oft – besonders im 2D-Fall – ist die Formulierung elektromagnetischerAufgaben in Bezug auf ein einziges Feld ökonomischer. Hierzu wird aus denbeiden Gleichungen (17) jeweils entweder E oder H eliminiert.So erhält man durch Anwendung von ∇× auf (17a) und Einsetzen von∇×H aus (17b)

∇×(∇×E) + iωµσE = 0 (18)

bzw. nach Division von (17b) durch σ > 0 und anschließender Anwendungvon ∇×

∇×(

1

σ∇×H

)+ iωµH = 0. (19)


Variationsformulierung der Helmholtz-Gleichung2D MT Gleichungen E-Polarisation

Annahme: alle Daten (speziell σ) sind in einer Koordinatenrichtung, dersog. Streichrichtung (hier x), uniform. Man unterscheidet zwei Spezialfälle:

E-Polarisation: Ist das elektrische Feld parallel zur Streichrichtung, so gilt

E =

Ex(y, z)00

und somit H =

0Hy(y, z)Hz(y, z)

.Mit

∇×E =

0∂zEx−∂yEx

bzw. ∇×(∇×E) =

−∂yyEx − ∂zzEx00

wird aus (18)

−∆Ex + iωµσEx = 0. (20)


Variationsformulierung der Helmholtz-Gleichung2D MT Gleichungen H-Polarisation

H-Polarisation: Zeigt H in die Streichrichtung, so erhält man mit

H =

Hx(y, z)00

, ∇×H =

0∂zHx

−∂yHx

,

∇×(

1

σ∇×H

)=

−∂y(σ−1∂yHx)− ∂z(σ−1∂zHx)00

aus (19) die skalare Gleichung

−∇·(

1

σ∇Hx

)+ iωµHx = 0. (21)

Die Komponenten des verbleibenden Feldes (E-Polarisation: Hy, Hz,H-Polarisation: Ey, Ez) werden mittels numerischer Differentiationberechnet.


Variationsformulierung der Helmholtz-GleichungMT Rechengebiet, Koeffizienten

Einfachste Situation: Rechteckgebiet Ω = [ymin, ymax]× [zmin, zmax]

Bei E-Polarisation gehört ein Teil des Lufthalbraums zu Ω(zmin < 0), bei H-Polarisation nicht (zmin = 0).Im Boden ist σ ortsabhängig, im Lufthalbraum konstant.

E-Polarisation

y

z

σ = σ(y, z)

σ ≡ 10−14S/mz=0

z=zmin

z=zmaxy=ymin y=ymax

H-Polarisation

y

z

σ = σ(y, z)

z=0

z=zmax

y=ymin y=ymax


Variationsformulierung der Helmholtz-GleichungMT Randbedingungen I

Grundannahme: im Bereich der lateralen Ränder ist σ = σ(z). Wirbetrachten hier den einfachen Fall, dass σ(z) sogar eine stückweisekonstante Funktion ist.

E-PolarisationIm Fall σ = σ(z) vereinfacht sich (20) zur gewöhnlichenDifferentialgleichung

E′′(z) + k2E(z) = 0, k2 = −iωµσ,

für die Feldkomponente E = Ex. Unter der Zusatzannahme σ ≡ consterhalten wir die allgemeine Lösung

E(z) = Aeikz +Be−ikz, A,B ∈ C.


Variationsformulierung der Helmholtz-GleichungMT Randbedingungen II

Wählt man für k =√−iωµσ denjenigen Zweig der Wurzelfunktion

mit Im(k) > 0 für ω > 0, so zieht die Abklingbedingung

E(z)→ 0 für z →∞

nach sich, dass A = E(0) und B = 0, und somit

E(z) = E0 eikz, z > 0, E0 = E(0).


Variationsformulierung der Helmholtz-GleichungMT Randbedingungen III

Zur Bestimmung von E0 beachte man, dass im Fall σ = σ(z) neben Hx

auch Hz verschwindet. Wie für E = Ex erhält man mit H0 := Hy(0) diez-Abhängigkeit

Hy(z) = H0 eikz

und durch Vergleich der ersten Komponenten der Gleichung (17b)

σEx = −∂zHy, d.h. Ex(z) =−ikσH0e

ikz, d.h. E0 =−ikσH0.

Wir erhalten somit die beiden lateralen Randwerte im Untergrund z > 0

Ex(ymin, z) =−ikσH0e

ikz, σ = σymin , k = kymin ,

Ex(ymax, z) =−ikσH0e

ikz, σ = σymin , k = kymin .

(22)


Variationsformulierung der Helmholtz-GleichungMT Randbedingungen IV, z < 0

Für den Bereich z < 0 beachten wir zunächst, dass aus den zweitenKomponenten von (17a) folgt

∂zEx = −iωµHy.

Hieraus und aus der Tatsache, dass im Lufthalbraum gilt Hy ≡ H0 folgtsomit für z < 0

Ex(z)− Ex(z = 0) =

∫ z

ζ=0∂zEx(ζ) dζ = −iωµ

∫ z

ζ=0H0 dζ = −iωµH0z.

Mit (22) erhalten wir somit an den lateralen Rändern im Lufthalbraum

Ex(z) = −iH0

(k

σ+ ωµz

), z < 0, (23)

wobei (die Konstante) H0 üblicherweise zu Eins normiert wird.


Variationsformulierung der Helmholtz-GleichungRandbedingungen V

An den horizontalen Rändern z = zmin und z = zmax werden die lateralenRandwerte durch ein kubisches Polynom φ in y interpoliert mit denNebenbedingungen, dass φ′(y) = 0 für y = ymin und y = ymax.Mit den Interpolationsbedingungen

φ(ymin) = φ0, φ′(ymin) = 0,

φ(ymax) = φ1, φ′(ymax) = 0

erhält man das Polynom

φ(y) = φ0 +(φ1−φ0)(3η2 − 2η3

), η =

y − ymin

ymax − ymin, ymin ≤ y ≤ ymax.


Variationsformulierung der Helmholtz-GleichungElliptische partielle Differentialgleichung zweiter Ordnung

Im Falle der Magnetotellurik (MT) erfüllt jeweils eine Komponente deselektrischen Feldes E und des Magnetfeldes H eine elliptischeDifferentialgleichung zweiter Ordnung der Form

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

wobei für

E-Polarisation : k = 1 und b = iωµσ und fürH-Polarisation : k = σ−1 und b = iωµ

gilt, wobei ω = 2πf die Kreisfrequenz und i die imaginäre Einheitbezeichnen.


Variationsformulierung der Helmholtz-GleichungRandwertaufgabe

Mit geeigneten Randbedingungen können wir folgende Randwertaufgabeformulieren

−∇·(k∇u(x)) + bu = 0, x ∈ Ω, (25a)

wobei Ω ⊂ Rd (d = 2, 3) ein zulässiges Gebiet sei. Der GebietsrandΓ := ∂Ω umfasst in diesem Fall nur ΓD. Wir stellen dieRandbedingungen

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

mit einer gegebenen, auf ΓD definierten Funktion g. Auch diese RWA(25) besitzt – unter geeigneten Voraussetzungen an das Gebiet Ω unddie Daten gD – eine eindeutig bestimmte klassische (d.h. in Ω zweimalstetig differenzierbare) Lösung.


Variationsformulierung der Helmholtz-GleichungL2-Innenprodukt und Divergenzsatz

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

(f, v) = 0 = −∫

Ω

(∇·(k∇u)v − buv)

)dx

= −∫

Ω

(∇·(vk∇u)− k∇u · ∇ v − buv

)dx

=

∫Ω

(k∇u · ∇ v + buv

)dx−

∫Γvk∂u

∂nds

= (k∇u,∇ v) + (bu, v),

wobei (·, ·) hier auch das L2-Innenprodukt vektorwertiger Funktionenauf Ω sowie (·, ·)Γ das L2-Innenprodukt auf Γ bezeichnen mögen. Wirwählen nun die Testfunktion v so, dass diese auf ΓD verschwindet.Insgesamt ergibt sich

(k∇u,∇ v) + (bu, v) = 0. (26)Antje Franke-Börner (TU Freiberg) Numerische Simulation mit finiten Elementen Sommersemester 2013 77

Variationsformulierung der Helmholtz-GleichungFunktionenräume

Wir setzen nun

a(u, v) := (k∇u,∇ v) + (bu, v) =

∫Ω

(k∇u · ∇ v + buv

)dx,

`(v) := 0.

Damit das Integral definiert ist, reicht es aus, dass die erstenAbleitungen von u und v quadratisch integrierbar sind, wir können alsou, v ∈ H1(Ω) wählen und setzen

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

Die Variationsformulierung von (25) lautet somit

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


Variationsformulierung der Helmholtz-Gleichung2D VLF FE-Lösung: σ = 0.01Sm−1, f = 20 kHz, E-Polarisation


Variationsformulierung der Helmholtz-Gleichung2D VLF FE-Lösung: σ = 0.01Sm−1, f = 20 kHz, H-Polarisation


Variationsformulierung der Helmholtz-Gleichung2D VLF FE-Lösung: σ1 = 0.01Sm−1, σ2 = 0.1Sm−1, f = 20 kHz, E-Polarisation


Variationsformulierung der Helmholtz-Gleichung2D VLF FE-Lösung: σ1 = 0.01Sm−1, σ2 = 0.1Sm−1, f = 20 kHz, H-Polarisation


Inhalt







Lagrange-ElementeLineare Dreiecks- (2D) und Tetraederelemente (3D)

1

3

2(0,0) (1,0)

(0,1)

1

4

2

3

(0,0,0) (1,0,0)

(0,1,0)

(0,0,1)

AnwendungSimulation skalarer Größen, z.B. 2D EM Feldkomponenten,elektrisches Potential2D MT, 2.5D CSEM, 2D und 3D Geoelektrik (DC)


Inhalt







Finite Elemente

Die Finite-Element-Methode (FEM) ist ein Galerkin-Verfahren, dessenzugehörige endlich-dimensionale Ansatz- und Testräume aus stückweisenPolynomen aufgebaut sind (d.h. bei Problemen in Rd sind diesPolynome in d Variablen).

Hierzu wird das zugrundeliegende Gebiet Ω in einfache Teilgebietezerlegt. Am häufigsten sind

Dreiecke und konvexe Vierecke (2D),Tetraeder, Hexaeder und Prismen (3D).

Bei krummlinig berandeten Gebieten können solche Zerlegungen dasGebiet lediglich durch polygonale bzw. polyedrische Gebieteapproximieren.

Das Zerlegen einer gegebenen Geometrie wird heute fast ausschließlichautomatisch durch sogenannte Netz- oder Gittergeneratoren (engl. meshgenerators) übernommen.


GittergeneratorDelaunay-Triangulierung

Delaunay-Triangulierung Erstellt aus einer Punktmenge einDreiecksnetz.

Umkreisbedingung Der Umkreis eines Dreiecks enthält keine weiterenPunkte der Punktmenge. Kleine Innenwinkel werdenvermieden.

Voronoi-Diagramm Ecken der Voronoi-Zellen sind Umkreismittelpunkteder Delaunay-Triangulierung.

Quelle: Wikipedia, „Delaunay-Triangulation“

3D Gebiete können mit Hilfe von Umkugeln in Tetraeder zerlegt werden.Antje Franke-Börner (TU Freiberg) Numerische Simulation mit finiten Elementen Sommersemester 2013 91

Zerlegungen

Unter einer Zerlegung Th des Gebiets Ω verstehen wir ein System vonTeilmengen K ⊂ Ω, welches folgende Bedingungen erfällt:

(Z1) Ω = ∪K∈ThK.

(Z2) Jedes K ∈ Th ist eine abgeschlossene Menge mitnichtleerem Inneren K.

(Z3) Für je zwei verschiedene K1,K2 ∈ Th gilt K1 ∩ K2 = ∅.(Z4) Jedes K ∈ Th besitzt einen Lipschitz-stetigen Rand ∂K.

Der Parameter h sei der maximale Durchmesser aller K ∈ Th.

Eine Zerlegung bezeichnen wir je nach Zusammenhang auch alsTriangulierung (nicht nur bei Dreiecken), Vernetzung, Netz oder Gitter.Es folgen einige grafische Beispiele von Zerlegungen.

Die einzelnen Teilgebiete werden Elemente genannt. (Später wird dieserBegriff jedoch erweitert!)


ZerlegungenBeispiele

Dreieckszerlegung des Äußeren eines Tragflächenquerschnitts.



Zerlegung aus Drei- und Vierecken.



Tetraedergitter einer 3D-Werkstücks.



3D-Gitter aus Quadern.



3D Tetraedervernetzung bei der FE-Analyse biologischen Gewebes.



Weitere Beispiele der Vernetzung komplexer Geometrien.


Konformität bei H1(Ω) I

Mit V h sei ein zunächst beliebiger endlichdimensionaler Raum von aufΩ definierten Funktionen bezeichnet. Entsprechend sei

PK := v|K : v ∈ V h

der durch sämtliche Einschränkungen von Funktionen aus V h auf Kaufgespannte Raum.Bei einer konformen FE-Diskretisierung ist V h ⊂ V erforderlich. BeiRandwertaufgaben zweiter Ordnung ist etwa V = H1(Ω) (bzw. einUnterraum hiervon). Eine Charakterisierung von Konformität liefertfolgender Satz.

Satz 11Sei Th eine Zerlegung des Gebietes Ω und V h ein endlichdimensionalerFunktionenraum. Gilt V h ⊂ C0(Ω) sowie PK ⊂ H1(K) für alleK ∈ Th, so gelten

V h ⊂ H1(Ω), sowie V h0 := v ∈ V h : v = 0 auf ∂Ω ⊂ H1

0 (Ω).


Konformität bei H1(Ω) II

Beweis: Aufgrund unserer Annahmen gilt bereits V h ⊂ L2(Ω). Damit auchV h ⊂ H1(Ω) müssen wir zeigen, dass jedes v ∈ V h schwache Ableitungen∂iv, i = 1, . . . , d besitzt, d.h. Funktionen vi ∈ L2(Ω) mit∫

Ω

vi φdx = −∫

Ω

v ∂iφ dx ∀φ, φ differenzierbar, φ|∂Ω = 0.

Elementweise gilt∫K

φ∂i(v|K) dx = −∫K

v|K ∂iφdx+

∫∂K

v|K φnK,i ds,

(nK,i die i-te Komponente der äußeren Einheitsnormalen längs ∂K).Summation über alle K ergibt (mit vi := ∂iv|K∀K)∫

Ω

φ vi dx = −∫

Ω

v ∂iφdx+∑

K∈Th

∫∂K

v|K φnK,i ds.

Die Summe verschwindet jedoch, da φ längs ∂Ω verschwindet und die(orientierten!) Randintegrale längs innerer Ränder je zweimal mitentgegengesetztem Umlaufsinn auftreten.


Konformität bei H1(Ω) III

Bemerkungen 12

1. Dass Stetigkeit auch notwendig ist, sieht man durch einenWiderspruchsbeweis.

2. Analog gilt V h ⊂ H2(Ω), falls PK ∈ H2(K) für alle K undV h ⊂ C1(Ω).

Fazit: Um konforme Diskretisierungen zu erhalten können wir also dieendlichdimensionalen Unterräume V h elementweise definieren undmüssen dabei nur die entsprechenden stetigen Übergänge derFunktionen bzw. deren Ableitungen gewährleisten.


Multivariate PolynomeVollständige Polynome

Für k ∈ N0 bezeichne Pk den Raum aller Polynome vom Grad ≤ k in dVariablen, d.h.

Pk :=

p : Rd → R : p(x) =

∑|α|≤k

bαxα

mit Koeffizienten bα, Multiindices α = (α1, . . . , αd) ∈ Nd0,|α| = α1 + · · ·+ αd sowie xα = xα1

1 · · ·xαdd .

Für M ⊂ Rd sei Pk(M) := p|M : p ∈Pk.

Es gilt dim Pk =(d+kk

). (dim Pk = dim Pk(M) sofern M 6= ∅)

Der Unterraum von Pk aus Polynomen in d Variablen vom exaktenGrad k (sog. homogene Polynome vom Grad k) besitzt die Dimension

dim Pk − dim Pk−1 =

(d+ k − 1

k

).


Simplexelemented-Simplices

Ein nichtentartetes Simplex K in Rd (kurz: d-Simplex) ist die konvexeHülle von d+ 1 nicht in einer Hyperebene gelegener Punkte ajd+1

j=1 ,den Ecken des d-Simplex:

K =

x =

d+1∑j=1

λjaj : 0 ≤ λj ≤ 1,

d+1∑j=1

λj = 1

.

Ein 2-Simplex ist ein Dreieck, ein 3-Simplex ein Tetraeder.

Ein Simplex ist genau dann nicht entartet, wenn die Spalten der Matrix

A :=

[a1 a2 · · · ad+1

1 1 . . . 1

]∈ R(d+1)×(d+1) (28)

linear unabhängig sind.


Baryzentrische KoordinatenEinführung

Für Punkte in einem Simplex K ⊂ Rd ist folgende Darstellung oftzweckmäßig: als baryzentrische Koordinaten λ = (λ1, . . . , λd+1)> einesPunktes x = (x1, . . . , xd)

> ∈ K bezeichnen wir den eindeutigenLösungsvektor λ des linearen Gleichungssystems

Aλ =

[x1

]mit Koeffizientenmatrix A aus (28).Sind a(−1)

i,j d+1i,j=1 die Einträge von A−1, so besitzen die baryzentrischen

Koordinaten die Darstellung

λi =

d∑j=1

a(−1)i,j xj + a

(−1)i,d+1, i = 1, . . . , d+ 1,

und sind damit affine Funktionen von x, d.h. λi = λi(x) ∈P1.


Baryzentrische KoordinatenInvarianz baryzentrische Koordinaten unter affinen Transformationen

Ist der Simplex K mit den Ecken a1, . . . ,ad+1 das Bild eines Simplex Kmit Ecken a1, . . . , ad+1 unter der Abbildung x 7→ x = Bx+ b mit einernichtsingulären Matrix B ∈ Rd×d und b ∈ Rd, so gilt für die zugehörigenbaryzentrischen Koordinaten[Ba1 . . . Bad+1

1 . . . 1

]λ =

[Bx1

]bzw.

[a1 . . . ad+1

1 . . . 1

]λ =

[x1

]also λ = λ.

Baryzentrische Koordinaten werden auch Dreiecks- bzw.Flächenkoordinaten im Fall d = 2 sowie Volumenkoordinaten im Falld = 3 genannt. Letztere Bezeichnungen rühren daher, dass λj dasFlächen- bzw. Volumenverhältnis zwischen K und dem durch(a1, . . . ,aj−1,x,aj+1, . . . ,ad+1) aufgespannten Simplex angibt.


Baryzentrische KoordinatenFlächenteilverhältnisse

a1 a2

a3

K

K = BK + b

a1

a2a3

x

x

λ1 = |K1||K|

= |K1||K|

λ2 = |K2||K|

= |K2||K|

λ3 = |K3||K|

= |K3||K|

Baryzentrische Koordinaten von x und x = Bx+ b in den affinen Dreiecken Kund K; die Flächen gleich eingefärbter Teildreiecke sind proportional.


Baryzentrische KoordinatenSchwerpunkt

a1

a2

a3

λ3 = 23

λ3 = 13

λ2 = 13

λ2 = 23

λ1 = 23

λ1 = 13

Linien konstanter baryzentrischer Koordinaten im Dreieck. Der Schwerpunktliegt bei λ1 = λ2 = λ3 = 1

3 .


Das lineare SimplexelementBasisfunktionen

Lemma 13Ein Polynom p ∈P1 in d Variablen ist durch dessen Funktionswertep(aj) an den d+ 1 Ecken eines d-Simplex im Rd eindeutig bestimmt.

Beweis: Zu zeigen ist, dass für jeden Vektor µ = (µ1, . . . , µd+1)> ∈ Rd+1 daslineare Gleichungssystem∑

|α|≤1

bαaαj = µj , j = 1, . . . , d+ 1,

eindeutig lösbar ist. Die zugehörige Koeffizientenmatrix ist quadratisch, somitist Existenz einer Lösung äquivalent mit deren Eindeutigkeit. Für diebaryzentrischen Koordinaten der Ecken gilt λi(aj) = δij (1 ≤ i, j ≤ d+ 1),weshalb das Polynom

p(x) :=

d+1∑i=1

µiλi(x)

das Gewünschte leistet.


Das lineare Simplexelement - Definition

Bemerkung 14Aus dem Beweis folgt insbesondere, dass jedes p ∈P1 bezüglich derbaryzentrischen Koordinaten eines Simplex die Darstellungp(x) =

∑d+1i=1 p(ai)λi(x) besitzt.

Das lineare Simplexelement in d Variablen ist nun wie folgt definiert:(i) Das Gebiet K ist ein nichtentartetes Simplex in Rd.(ii) Der endlichdimensionale Funktionenraum PK auf K ist P1(K).(iii) Die Freiheitsgrade ΨK , welche ein Polynom p ∈PK eindeutig

festlegen, sind dessen Funktionswerte an den Ecken, symbolischΨK = p(aj) : 1 ≤ j ≤ d+ 1.

Nach Bemerkung 14 bilden die Funktionenφj(x) = λj(x), j = 1, . . . , d+ 1,

eine Basis von PK mit der Eigenschaft φi(aj) = δij . Eine solche Basisheißt nodale Basis.


Das quadratische SimplexelementEckknoten

Seienaij := 1

2(ai + aj), 1 ≤ i < j ≤ d+ 1

die Seitenmittelpunkte des d-Simplex K. Es gilt

λk(aij) = 12(δki+δkj) =

12 ak und aij liegen auf gemeinsamer Kante0 sonst.

Damit gelten für die Funktionen φi ∈P2 definiert durch

φi(x) := λi(x)(2λi(x)− 1), i = 1, . . . , d+ 1,

die Beziehungen

φi(aj) = δij , φi(ajk) = 0, j < k.


Das quadratische SimplexelementMittelknoten

Für die d(d+ 1)/2 weiteren Funktionen φij ∈P2 definiert durch

φij(x) := 4λi(x)λj(x), 1 ≤ i < j ≤ d+ 1,

gilt

φij(aij) = 1,

φij(ak`) = 0, falls k 6= i oder ` 6= j,

φij(ak) = 0, für alle k.

Wegen dim P2 = (d+ 1)(d+ 2)/2 bilden damit die (d+ 1)(d+ 2)/2Funktionen

φ1, . . . , φd+1, φ1,2, . . . , φd,d+1

die nodale Basis von P2 bezüglich der Knoten

ai, 1 ≤ i ≤ d+ 1, und aij , 1 ≤ i < j ≤ d+ 1.


Das quadratische SimplexelementFreiheitsgrade

Insbesondere besitzt jedes Polynom p ∈P2 die Darstellung

p(x) =

d∑i=1

p(ai)λi(x)(2λi(x)− 1) +∑i<j

p(aij) · 4λi(x)λj(x).

Wir erhalten so das quadratische Simplexelement definiert durch dasSimplex K, den Funktionenraum P2(K) und als Freiheitsgrade dieFunktionswerte an den Ecken und Seitenmittelpunkten, d.h.

ΨK =p(ai), 1 ≤ i ≤ d+ 1;

p(aij), 1 ≤ i < j ≤ d+ 1.

a1

a2a3

a12a13

a23

Knoten im quadratischen DreieckAntje Franke-Börner (TU Freiberg) Numerische Simulation mit finiten Elementen Sommersemester 2013 112

Das quadratische SimplexelementNodale Basisfunktionen

Beispiel: Die nodalen Basisfunktionen im Dreieckelement vom Gradzwei.


Das kubische SimplexelementBasisfunktionen

Mit den Bezeichnungen

aiij := 13(2ai + aj), i 6= j,

aijk := 13(ai + aj + ak), i < j < k,

gilt analog für Polynome p ∈P3 die Beziehung

p =∑i

p(ai)λi(3λi − 1)(3λi − 2)

2+∑i 6=j

p(aiij)9λiλj(3λi − 1)

2

+∑i<j<k

p(aijk) · 27λiλjλk,

was auf das kubische Simplexelement mit Funktionenraum P3(K) undals Freiheitsgrade die Werte an den ai, aiij(i 6= j) und aijk(i < j < k)führt.


Das kubische SimplexelementKnoten

a1

a2a3

a221

a112a113

a331

a332 a223

a123

Knoten im kubischen Dreieck


SimplexelementeAllgemeiner Fall

Satz 15Sei K ein d-Simplex mit Ecken ajd+1

j=1 und p ∈Pk für ein k ∈ N.Dann ist p eindeutig bestimmt durch dessen Funktionswerte an denPunkten der Knotenmenge

Nk(K) :=

x =

d+1∑j=1

λjaj :

d+1∑j=1

λj = 1,

λj ∈ i/k, i = 0, . . . , k, j = 1, . . . , d+ 1

.

Beweis: Wegen |Nk(K)| = dim Pk ist durch die Forderungenφi ∈Pk, φi(xj) = δij , 1 ≤ i, j ≤ |Nk(K)|

die nodale Basis von Pk eindeutig bestimmt.


Finite-Element-Räume basierend auf SimplicesZerlegung Th

Zur Erinnerung: Wir betrachten eine Zerlegung Th eines polyedrischenGebiets Ω ⊂ Rd in disjunkte Simplices sodass die Eigenschaften(Z1)–(Z4) erfüllt sind:

(Z1) Ω = ∪K∈ThK.

(Z2) Jedes K ∈ Th ist eine abgeschlossene Menge mitnichtleerem Inneren K.

(Z3) Für je zwei verschiedene K1,K2 ∈ Th gilt K1 ∩ K2 = ∅.

(Z4) Jedes K ∈ Th besitzt einen Lipschitz-stetigen Rand ∂K.


Finite-Element-Räume basierend auf SimplicesRänder

Diesen Forderungen fügen wir eine weitere hinzu:(Z5) Jede Kante/Fläche eines Simplex K ∈ Th ist entweder

Teil des Randes ∂Ω oder gleichzeitig Kante/Fläche einesanderen Simplex von Th.

Erlaubt. Verboten.


Finite-Element-Räume basierend auf Simplices

Ein zu einer Zerlegung Th von Ω und einem finiten Element gehörenderFE-Raum V h besteht nun aus Funktionen auf Ω, deren Enschränkung aufjedes Element K ⊂ Th zu PK gehört.

Satz 16Sei V h der zu einer zulässigen Zerlegung Th von Ω und entweder demlinearen, quadratischen oder kubischen Simplex gehörende FE-Raum.Dann gilt

V h ⊂ C(Ω) ∩H1(Ω).

Beweis: Zu zeigen ist nur, dass die Funktionen v ∈ V h stetige Übergängezwischen den Elementen besitzen. Seien K1,K2 ∈ Th zwei benachbarte Simplicesmit gemeinsamer Fläche K ′ sowie v1 := v|K1

, v2 := v|K2für ein v ∈ V h. Die

Einschränkung von v1 bzw. v2 auf die Fläche K ′ ist jeweils ein Polynom vomGrad k in d− 1 Variablen (k = 1, 2, 3). Diese beiden Polynome stimmen an denauf K ′ liegenden Knoten überein. Die Anzahl Knoten auf (dem d− 1-Simplex)K ′ ist aber genau die erforderliche, um ein Polynom vom Grad k in d− 1

Variablen eindeutig festzulegen, also stimmt v1 auf K ′ mit v2 überein.


Affine Familien

Ziel des weiteren Vorgehens ist es nun, ein finites Element zunächst aufeinem Referenzelement K zu definieren und dann beliebige Elemente Keiner Triangulierung durch affine Abbildung von K nach K abzuleiten.

Vorteile dieses Ansatzes:Auf dem Referenzelement besitzt ein finites Element eine besonderseinfache Beschreibung.Berechnungen werden durch Transformation auf dasReferenzelement einheitlicher.Es genügt, theoretische Aussagen (z.B.Approximationseigenschaften) nur für das Referenzelement zubeweisen.

Wir bescheiben in diesem Abschnitt, wie dieser Ansatz im Allgemeinenangewandt werden kann.


Affine FamilienDefinition finites Element

Folgende abstrakte Definition eines finiten Elements geht auf Ciarlet zurück:

Definition 17Ein finites Element ist ein Tripel (K,P,Ψ); hierbei sind(i) K eine abgeschlossene Teilmenge des Rd mit nichtleerem Inneren und

Lipschitz-Rand,(ii) P ein linearer Raum auf K definierter Funktionen sowie(iii) Ψ eine endliche Menge linear unabhängiger linearer Funktionale

Φ1, . . . ,Φn auf P derart, dass jede Funktion φ ∈ P eindeutig durch dieWerte Φ1(φ), . . . ,Φn(φ) festgelegt ist. (Man sagt auch, Ψ seiunisolvent bezüglich P .)

Bemerkung 18Durch Ψ ist eine Basis φ1, . . . , φn von P ausgezeichnet durchΦi(φj) = δij , i, j = 1, . . . , n, die wir in Beispielen bereits als nodale Basisbezeichnet haben. Die Basen Φj und φj sind also zueinander dual.


Affine FamilienBedeutung der Freiheitsgrade

Dass es bei einem finiten Element auch auf die Freiheitsgrade ankommtveranschaulicht folgendes Beispiel: wir betrachten ein Dreieckelement mitPolynomraum P1 und als Freiheitsgrade die Funktionswerte an jeweils dreiKnoten, die wir wie folgt wählen.

Nur bei der ersten Wahl der Knoten (und damit der Freiheitsgrade) erhaltenwir einen FE-Raum mit stetigen Übergängen zwischen benachbartenDreiecken.

Insbesondere erhalten wir so drei verschiedene finite Elemente bei jeweilsgleichem K und P .


Affine FamilienDefinition affin äquivalent

Definition 19Zwei finite Elemente (K,P,Ψ) und (K, P , Ψ) heißen affin äquivalent,falls folgende Bedingungen erfüllt sind:(a) Es existiert eine affine Abbildung F : K → K, x 7→ Bx+ b, mit

B ∈ Rd×d nichtsingulär, b ∈ Rd sodass K = F (K),(b) Für die Funktionenräume gilt P = φ F−1 : φ ∈ P.(c) Für die Freiheitsgrade gilt Ψ = Φ · F−1 : Φ ∈ Ψ.

Proposition 20Affine Äquivalenz finiter Elemente ist eine Äquivalenzrelation.


Affine FamilienBemerkungen

Bemerkungen 21

1. Bedingung (c) besagt, dass die Funktionalmenge Ψ das Bild von Ψist unter der Abbildung

Φ 7→ Φ mit Φ(φ) = Φ(φ F−1) ∀φ ∈ P .

2. Die Zuordnung in (b) zwischen Funktionen auf K und solchen aufK über Verkettung mit F−1 wird auch „Pull-back“ genannt,geschrieben (F−1)∗(φ) := φ F−1.Damit lautet (b) P = (F−1)∗P oder P = F ∗P .

3. Analog nennt man die Zuordnung in (c) zwischen linearenFunktionalen von Funktionen auf K zu solchen von Funktionen aufK auch „Push-Forward“, geschrieben((F−1)∗Φ

)(φ) := Φ

((F−1)∗(φ)

). Eine entsprechende Formulierung

von (c) lautet also Ψ = (F−1)∗Ψ bzw. Ψ = F∗Ψ.Antje Franke-Börner (TU Freiberg) Numerische Simulation mit finiten Elementen Sommersemester 2013 124

Affine FamilienDefinition affine Familie

Beispiel: Lagrange-Elemente auf Simplices, deren Knotenmengen durchdie affine Gebietsabbildung F ineinander überführt werden, sind affinäquivalent.

Definition 22Eine Menge von finiten Elementen heißt affine Familie, falls alleElemente affin äquivalent zu einem Referenzelement (K, P , Ψ) sind.(Letzteres braucht selbst nicht zur Familie zu gehören.)

Bei d-Simplexelementen wählt man meist als Referenzelement (-Gebiet)den Einheitssimplex mit Ecken im Ursprung und an der Spitze der dKoordinaten-Einheitsvektoren.


Der InterpolationsoperatorDefinition

Der Interpolationsoperator ist ein wichtiges Hilfsmittel bei derKonstruktion von FE-Räumen aus einzelnen Elementen sowie bei derKonvergenzanalyse.

Definition 23Sei (K,P,Ψ) ein finites Element und φ1, . . . , φn die nodale Basis von Pzu Ψ. Ist u eine Funktion auf K, für die alle Freiheitsgrade Φj ∈ Ψdefiniert sind, so wird die lokale Interpolierende IKu ∈ P definiert durch

IKu =

n∑j=1

Φj(u)φj .


Der InterpolationsoperatorBeispiel

Beispiel: Sei K das lineare Dreieckelement bezüglich des Dreiecks Kmit den Ecken (0, 0), (1, 0) und (0, 1). Die nodale Basis ist gegebendurch

φ1(x, y) = 1− x− y, φ2(x, y) = x, φ3(x, y) = y.

Die lokale Interpolierende der Funktion u(x, y) = exy lautet damit(IKu)(x, y) = e0φ1(x, y)+e0φ2(x, y)+e0φ3(x, y) = (1−x−y)+x+y ≡ 1.

Proposition 24IK ist linear und es gilt Φj(IKu) = Φj(u), j = 1, . . . , n.

Mit anderen Worten: IKu ist diejenige Funktion in P mit denselbenFreiheitsgraden wie u.

Proposition 25Für alle p ∈ P gilt IKp = p, d.h. IK ist eine Projektion auf P .


Der globale Interpolationsoperator

Proposition 26Sind (K,P,Ψ) und (K, P , Ψ) zwei affin äquivalente finite Elementebezüglich der affinen Abbilding F : K → K und bezeichnen IK bzw. IK diezugehörigen lokalen Interpolationsoperatoren, so gilt IK F

∗ = F ∗ IK .

Definition 27Sei Th eine zulässige Zerlegung des polyedrischen Gebiets Ω ⊂ Rd und seiauf jedem Teilgebiet K ∈ Th ein finites Element (K,P,Ψ) definiert. Sindfür eine Funktion u : Ω→ R alle Freiheitsgrade ΨK ,K ∈ Th definiert, sodefinieren wir den globalen Interpolationsoperator Ih durch

(Ihu)|K = IKu, für alle K ∈ Th.

Man nennt ein finites Element Ck-Element, falls die zugehörigenglobalen Interpolierenden stets k mal stetig differenzierbar sind.Lagrange-Elemente sind (bei geeigneter Knotenwahl) C0-Elemente.


Inhalt







Das Galerkin-GleichungssystemDiskrete Variationsaufgabe

Das Galerkin-Verfahren zur Approximation der Lösung derPoissongleichung (6) besteht darin, Ansatz- und Testraum durchendlichdimensionale Räume S h bzw. V h (gleicher Dimension) zuersetzen.

Hierbei ist h > 0 der Diskretisierungsparameter.

Wir betrachten zunächst die homogenisierte Variationsaufgabe (8) undeinen n-dimensionalen Unterraum V h ⊂ V . (Oder, äquivalent, den FallgD ≡ 0.)

Die diskrete Variationsaufgabe lautet somit

Bestimme uh ∈ V h sodass a(uh, v) = `(v) ∀v ∈ V h. (29)


Das Galerkin-GleichungssystemMatrix-Vektor-Form

Ist φ1, φ2, . . . , φn eine Basis von V h sowie uh =∑n

j=1 ujφj , so ist(29) äquivalent mit

n∑j=1

uj a(φj , φi) = `(φi), i = 1, 2, . . . , n,

oder, mit A ∈ Rn×n gegeben durch [A]i,j = a(φj , φi), b ∈ Rn durch[b]i = `(φi) sowie u ∈ Rn durch [u]i = ui, zu dem Galerkin-System

Au = b. (30)

Beachte:Die diskrete Variationsaufgabe (29) bzw. (30) besitzt eineeindeutige Lösung.Die Galerkin-Matrix A aus (30) ist symmetrisch und positiv-definit.


Das Galerkin-GleichungssystemVerwendung der nodalen Basis

Wir beschäftigen uns nun mit der Aufstellung des Gleichungssystems(30) anhand des Beispiels linearer Dreieckelemente.

Funktionen aus V h sind eindeutig durch ihre Funktionswerte an denKnoten der Triangulierung (Ecken der Dreiecke) bestimmt. Bei linearenDreiecken stellen diese Funktionswerte sämtliche Freiheitsgrade dar.

Bei V h sind dies alle Knoten bis auf die, die auf ΓD liegen. Ihre Anzahlsei n.Eine besonders nützliche Basis φ1, . . . , φn, die sog. nodale Basis, istcharakterisiert durch

φj(xi) = δij i, j = 1, . . . , n.

Ist N h = x1, . . . , xn die Menge der Knoten mit xj 6∈ ΓD, so giltsuppφj = K ∈ T h : xj ∈ K.


Das Galerkin-GleichungssystemBeispiele zu Trägern nodaler Basisfunktionen

Triangulierung des L-förmigen Gebiets aus Abschnitt 1 (Beispiel 2) mit demTräger dreier nodaler Basisfunktionen.


Das Galerkin-Gleichungssystem

Konsequenz für das Galerkin-System (30):

[b]i = `(φi) =

∫suppφi

fφi dx+

∫suppφi∩ΓN

hφi ds

[A]i,j = a(φj , φi) =

∫suppφi∩suppφj

∇φi · ∇φj dx

Insbesondere: die Galerkin-Matrix A ist dünn besetzt.


Aufbau des Galerkin-SystemsVorteile

Großer Vorteil der FEM:Aufbau des FE Gleichungssystems – man nennt diesen VorgangAssemblieren – verläuft auch bei komplizierteren Problemen stetsnach dem gleichen Grundschema.Parametrisierung durch spezielle Eigenschaften der jeweiligenAufgabe.Grundbausteine der Assemblierung stets dieselben, kann bei derSoftwareumsetzung ausgenutzt werden (z.B. Objektorientierung).

Wir stellen nun die Grundschritte der Assemblierung für unserModellproblem zusammen. Das Vorgehen besteht darin, die Rechnungauf Operationen auf den einzelnen Elementen zurückzuführen, derenErgebnisse dann zusammengefügt (assembliert) werden.


Aufbau des Galerkin-SystemsVorgehensweise

Übliches Vorgehen:(1) Ignoriere zunächst die wesentlichen Randbedingungen, d.h. ein

gößerer Raum V h ⊃ V h wird zugrundegelegt, mit nodaler Basis

φ1, φ2, . . . , φn, φn+1, . . . , φn, wobei

n− n die Anzahl der Knoten auf ΓD ist.Liefert A ∈ Rn×n, b ∈ Rn.

(2) Eliminiere zum Schluß die Freiheitsgrade der wesentlichenRandbedingungen.Liefert A, b.

Zunächst naheliegend für (1) ist sukzessives Abarbeiten allerBasisfunktionen (= Einträge in A, b).

Aber: Lage und Anordnung der Träger variieren stark. Deshalb ist eseinfacher elementweise vorzugehen.


Aufbau des Galerkin-SystemsReduktion auf Elemente

Sei K ∈ T h: dann gilt für i, j = 1, 2, . . . , n:

a(φj , φi) =

∫Ω∇φj · ∇φi dx =

∑K∈T h

∫K∇φj · ∇φi dx =:

∑K∈T h

aK(φj , φi),

`(φi) =∑K∈T h

(∫Kfφi dx+

∫K∩ΓN

hφi ds

)=:

∑K∈T h

`K(φi).

Mit

[AK ]i,j := aK(φj , φi) i, j = 1, 2, . . . , n,

[bK ]i := `K(φi), i = 1, 2, . . . , n,

folgt alsoA =

∑K∈T h

AK , b =∑K∈T h

bK .


Aufbau des Galerkin-SystemsElementtabelle

Da jedes Element nur zum Träger dreier Basisfunktionen gehört, sindnur (höchstens) neun bzw. drei Einträge von AK bzw. bK von Nullverschieden.Welche Einträge dies sind kann durch Nachschlagen in einerElementtabelle ermittelt werden:

[ET (i, j)]i=1,2,3;j=1,...,nK :

Element K1 K2 . . . KnK

erster Knoten i(1)1 i

(2)1 . . . i

(nK)1

zweiter Knoten i(1)2 i

(2)2 . . . i

(nK)2

dritter Knoten i(1)3 i

(2)3 . . . i

(nK)3

Hierbei sei nK die Anzahl der Elemente in T h.Diese führt neben der bisherigen globalen Knotennumerierungx1, x2, . . . , xn in jedem Element zusätzlich eine lokale Nummerierungx

(K)1 , x(K)

2 , x(K)3 der zu K gehörenden Knoten ein.


Globale Nummerierung der Knoten (rot) und der Elemente (schwarz)einer Triangulierung des L-Gebiets.


Aufbau des Galerkin-SystemsElementmatrizen/-vektoren

Damit sind die von Null verschiedenen Teilmatrizen bzw. -vektoren von AK undbK gegeben durch

AK :=

aK(φ(K)1 , φ

(K)1 ) aK(φ

(K)2 , φ

(K)1 ) aK(φ

(K)3 , φ

(K)1 )

aK(φ(K)1 , φ

(K)2 ) aK(φ

(K)2 , φ

(K)2 ) aK(φ

(K)3 , φ

(K)2 )

aK(φ(K)1 , φ

(K)3 ) aK(φ

(K)2 , φ

(K)3 ) aK(φ

(K)3 , φ

(K)3 )

, bK :=

`K(φ(K)1 )

`K(φ(K)2 )

`K(φ(K)3 )

.

Trägt K in der Nummerierung der Elemente den Index k, so ist die Zuordnungder lokalen Nummerierung φ(K)

i i=1,2,3 der zu K gehörenden Basisfunktionenzur globalen Nummerierung φjnj=1 gegeben durch

φ(K)i = φj , j = ET (i, k), i = 1, 2, 3.

Man bezeichnet AK und bK auch als Elementsteifigkeitsmatrix bzw.Elementlastvektor.


Aufbau des Galerkin-SystemsAssemblierungsalgorithmus

Nach diesen Überlegungen erhalten wir folgenden Algorithmus4 zurAssemblierung:

Initialisiere A := O, b := 0foreach K ∈ Th

berechne AK und bKk := [Index des Elementes K]j1 := ET (1, k), j2 := ET (2, k), j3 := ET (3, k)

A([j1j2j3], [j1j2j3]) := A([j1j2j3], [j1j2j3]) +AK

b([j1j2j3]) := b([j1j2j3]) + bKend

4Hier wird folgende an MATLAB angelehnte Notation verwendet:

A([j1j2j3], [j1j2j3]) =

aj1,j1 aj1,j2 aj1,j3aj2,j1 aj2,j2 aj2,j3aj3,j1 aj3,j2 aj3,j3

, b([i1i2i3]) =

bi1bi2bi3

.Antje Franke-Börner (TU Freiberg) Numerische Simulation mit finiten Elementen Sommersemester 2013 141

Aufbau des Galerkin-SystemsReferenzelement

Hilfreich für Implementierung und Analyse: Bezug auf einReferenzelement K ⊂ R2. Für alle K ∈ T h gilt dann K = FK(K) mit

FK : K → K, K 3 ξ 7→ x ∈ K, x = FK(ξ) = BKξ + bK .

Bei Dreieckelementen üblich: Einheitssimplex

K = (ξ, η) ∈ R2 : 0 ≤ ξ ≤ 1, 0 ≤ η ≤ 1− ξ

Für jedes Dreieck K ∈ T h ist die affine Abbildung FK bestimmt durchdie Abbildungsvorschriften

(1, 0) 7→ (x1, y1),

(0, 1) 7→ (x2, y2),

(0, 0) 7→ (x3, y3), d.h.



K

ξ

η

(0, 0) (1, 0)

(0, 1)

x

y

FK

K

(x1, y1)

(x2, y2)

(x3, y3)

[xy

]=

[x1 − x3 x2 − x3

y1 − y3 y2 − y3

]︸︷︷︸

BK

[ξη

]+

[x3

y3

]︸︷︷︸bK



Lokale (nodale) Basisfunktionen in K:

φ1(ξ, η) = ξ, φ2(ξ, η) = η, φ3(ξ, η) = 1− ξ − η, (ξ, η) ∈ K.

Durch die Zuordnung

φ 7→ φ := φ F−1K , d.h. φ(x) := φ(ξ(x)) = φ(F−1

K (x))

wird jeder Funktion φ auf K eine Funktion φ auf K zugeordnet.Lokale Basisfunktionen auf K:

φj = φj F−1K : K → R, j = 1, 2, 3.


Aufbau des Galerkin-SystemsRückführung der Integration auf das Referenzelement

Die bei der Assemblierung der Galerkin-Matrix anfallenden Integralewerden ebenfalls auf das Referenzelement zurückgeführt. Dies ist auchfür die Verwendung von Quadraturformeln hilfreich.Aus φ(x) = φ(ξ(x)) folgt (Kettenregel5)

∇φ =

[φxφy

]=

[φξξx + φηηxφξξy + φηηy

]=

[ξx ηxξy ηy

] [φξφη

]= (DF−1

K )>∇φ.

Wegen x = FK(ξ) = BKξ + bK , d.h. DFK = BK ,

ξ = F−1K (x) = B−1

K (x− bK), d.h. DF−1K = B−1

K

folgt schließlich∇φ = B−>K ∇φ.

5∇ bedeutet, dass nach den Variablen ξ und η differenziert wird.



Hierbei sind (lineares Dreieckelement)

| detBK | = 2|K|,

B−>K =1

2|K|

[y2 − y3 y3 − y1

x3 − x2 x1 − x3

]=

1

2|K|

[y2 − y3 x3 − x2

y3 − y1 x1 − x3

]>,

[∇φ1 ∇φ2 ∇φ3

]=

[1 0 −10 1 −1

].


Aufbau des Galerkin-SystemsEinbringung wesentlicher Randbedingungen

Bei Lagrange-Elementen läßt sich die wesentliche Randbedingungu(x) = gD(x) ∀x ∈ ΓD

am einfachsten durch die Forderunguh(x) = gD(x) x ∈ N h \N h,

umsetzen, wobei mit N h := x1, . . . ,xn,xn+1, . . . ,xn die Mengealler Knoten des zu T h gehörenden FE-Raums bezeichnet sei.Zerlegt man N h in N h = NI ∪NW mitN hI := N h = x1, . . . ,xn, N h

W := N h \N h = xn+1, . . . ,xn

und nummeriert die Freiheitsgrade von uh entsprechend, so zerfällt dasbisher assemblierte Galerkin-System Au = b in[

AII AIW

AWI AWW

] [uIuW

]=

[bIbW

]. (33)


Aufbau des Galerkin-SystemsEinbringung wesentlicher Randbedingungen

Hierbei beschreiben die Teilmatrizen AII und AWW die Kopplung derden Knoten aus N h

I bzw. N hW zugeordneten Freiheitsgrade jeweils

untereinander, sowie AIW und AWI die gegenseitigen Kopplungen.

Bezeichnet gD ∈ Rn−n den Vektor mit den KomponentengD,i = gD(xi), so lauten unter Hinzunahme der Nebenbedingungen

uW = gD

die diskreten Gleichungen (33)[AII AIW

O I

] [uIuW

]=

[bIgD

].

Die nicht durch die wesentlichen Randbedingungen festgelegtenFreiheitsgrade uI ergeben sich somit aus dem reduzierten linearenGleichungssystem

AIIuI = bI −AIWgD.


Inhalt







Aufgabenstellung

Gesucht ist die Lösung u in Ω = (0, 1) für

−d2u

dx2= f in (0, 1) mit

u(0) = 0 und du

dx(1) = 0.

Die Variationsformulierung lautet: Gesucht ist u ∈ V , so dassa(u, v) = l(v) mit

a(u, v) =

∫ 1

0

du

dx

dv

dxdx und l(v) =

∫ 1

0fvdx.

Die Galerkin-Diskretisierung uh =∑N

j=1 ujφj liefert die diskreteVariationsaufgabe: Gesucht ist uh ∈ V h, so dass

N∑j=1

uj

∫ 1

0

dφidx

dφjdx

dx =

∫ 1

0fφidx ∀i = 1, ..., N.


Lineare Simplex-Elemente in 1D

1. Berechnen Sie die baryzentrischen Koordinatenfunktionen zum1-Simplex!

2. Wie ergeben sich daraus die zugehörigen linearen Basisfunktionen?3. Wie lässt sich die lokale Interpolierende Iku darstellen?

x1 x2

λ1(x) λ2(x)

e



Die baryzentrischen Koordinatenfunktionen

Aus Aλ =

[x1

]mit A =

[x1 x2

1 1

]ergibt sich

λ1(x) =x2 − xx2 − x1

und λ2(x) =x− x1

x2 − x1.

Basisfunktionen und InterpolierendeFür die Basisfunktionen gilt

φ1 = λ1(x) und φ2 = λ2(x).

Damit lässt sich uh = Iku auf e darstellen als

Iku = u1φ1(x) + u2φ2(x).


Zerlegung des Intervalls (0, 1) in drei Elemente

Wir zerlegen das Intervall (0, 1) in drei Elemente (Nk = 3) mit ∆x = 13

und f ≡ 1:

0 13

23

1x1

x2 x3 x4

φ1 φ2 φ3 φ4

e1 e2 e3

k = 1 k = 2 k = 3


Steifigkeitsmatrix und Lastvektor

Summation der Elementsteifigkeitsmatrizen akij und -lastvektoren bkiliefert die globale Steifigkeitsmatrix und den globalen Lastvektor:

N∑j=1

uj

Nk∑k=1

∫ek

dφidx

dφjdx

dx =

Nk∑k=1

∫ek

φidx ∀i = 1, ..., N mit

Ak = akij(φi, φj) =

∫ek

dφidx

dφjdx

dx und bk = bki (φi) =

∫ek

φidx sowie

aij =

Nk∑k=1

akij und bi =

Nk∑k=1

bki ,

wobei Nk = 3 die Anzahl der Elemente ist. Es gilt nun dasGleichungssystem Au = b zu lösen mit

A = aij(φi, φj), b = bi(φi),u = ui, i, j = 1, ..., N.


Elementsteifigkeitsmatrix und -lastvektor

Ak =

[ ∫ek

dφkdx

dφkdx dx

∫ek

dφkdx

dφk+1

dx dx∫ek

dφk+1

dxdφkdx dx

∫ek

dφk+1

dxdφk+1

dx dx

]

bk =

[ ∫ekφkdx∫

ekφk+1dx

]

Berechnen Sie φk, φk+1 sowie dφkdx ,

dφk+1

dx (affine Abbildung)!



Ak =

[ ∫ek

dφkdx

dφkdx dx

∫ek

dφkdx

dφk+1

dx dx∫ek

dφk+1

dxdφkdx dx

∫ek

dφk+1

dxdφk+1

dx dx

]

bk =

[ ∫ekφkdx∫

ekφk+1dx

]Berechnen Sie φk, φk+1 sowie dφk

dx ,dφk+1

dx !

φk(x) =xk2 − xk

xk2 − xk1=

xk+1 − xxk+1 − xk

=k∆x− x

∆x

φk+1 =xk − xk1xk2 − xk1

=x− xk

xk+1 − xk=x− (k − 1)∆x

∆x

dφkdx

= − 1

∆xdφk+1

dx=

1

∆xAntje Franke-Börner (TU Freiberg) Numerische Simulation mit finiten Elementen Sommersemester 2013 157


Ak =

[∫ek

(− 1∆x)(− 1

∆x)dx∫ek

(− 1∆x)( 1

∆x)dx∫ek

( 1∆x)(− 1

∆x)dx∫ek

( 1∆x)( 1

∆x)dx

]=

1

∆x

[1 −1−1 1

]

bk =

[ ∫ek

k∆x−x∆x dx∫

ek

x−(k−1)∆x∆x dx

]= ?

Berechnen Sie die Einträge von bk durch bestimmte Integration!



Ak =

[∫ek

(− 1∆x)(− 1

∆x)dx∫ek

(− 1∆x)( 1

∆x)dx∫ek

( 1∆x)(− 1

∆x)dx∫ek

( 1∆x)( 1

∆x)dx

]=

1

∆x

[1 −1−1 1

]

bk =

[ ∫ek

k∆x−x∆x dx∫

ek


]=

∆x

2

[11

]


Gleichungssystem

Elementtabelle und Koordinatenvektor

ET =

[1 2 32 3 4

], p = ∆x

[0 1 2 3

]=

1

3

[0 1 2 3

]

Zuordnung zwischen globaler und lokaler NummerierungA ∈ R4×4 und b ∈ R4

A1 →[a11 a12

a21 a22

], A2 →

[a22 a23

a32 a33

], A3 →

[a33 a34

a43 a44

]b1 →

[b1b2

], b2 →

[b2b3

], b3 →

[b3b4

]


Assemblierung des Gleichungssystems

1. Initialisierung von A und b als Null-Matrix bzw. -Vektor2. Addition der Elementmatrizen und -vektoren unter

Berücksichtigung der Elementtabelle

A =1

∆x

1 −1 0 0−1 1 0 00 0 0 00 0 0 0

+1

∆x

0 0 0 00 1 −1 00 −1 1 00 0 0 0

+1

∆x

0 0 0 00 0 0 00 0 1 −10 0 −1 1

=1

∆x

1 −1 0 0−1 2 −1 00 −1 2 −10 0 −1 1

b =∆x

2

1100

+∆x

2

0110

+∆x

2

0011

=∆x

2

1221


Dirichlet-Randbedingung

Gleichungssystem

Au = b→ 1

∆x

1 −1 0 0−1 2 −1 00 −1 2 −10 0 −1 1

u1

u2

u3

u4

=∆x

2

1221

, u1 = u(0) = 0

Reduziertes Gleichungssystem

1

∆x

1 O−1 2 −1 00 −1 2 −10 0 −1 1

u1

u2

u3

u4

=∆x

2

0221


Lösung des Gleichungssystems

Reduziertes GleichungssystemEs bleibt zu lösen:

1

∆x

2 −1 0−1 2 −10 −1 1

u2

u3

u4

=∆x

2

221

Lösung ?




1

∆x

2 −1 0−1 2 −10 −1 1

u2

u3

u4

=∆x

2

221

Lösung

u2 =5

2(∆x)2, u3 = 4(∆x)2, u4 =

9

2(∆x)2


Vergleich mit der exakten Lösung

Aus zweimaliger Integration von −∆u = 1 unter Berücksichtigung derRandbedingungen erhalten wir die exakte Lösung u(x) = x− 1

2x2.

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 10

0.05

0.1

0.15

0.2

0.25

0.3

0.35

0.4

0.45

0.5

x

u(x)

uh1, Nk = 3

uh2, Nk = 10

u(x) = x−1/2 x2


Inhalt







Konvergenz

Im Lemma von Céa (Satz 10) haben wir bereits gesehen, dass dieGalerkin-Approximation der Lösung einer Variationsaufgabe, welche denVoraussetzungen des Lax-Milgram-Lemmas genügt, bis auf eineKonstante die Bestapproximation der Lösung im verwendetenUnterraum liefert.

Somit konvergiert die FE-Approximation gegen die Lösung derentsprechenden RWA, sofern diese bei Vergrößerung des FE-Raumes(durch Verfeinerung des Gitters und/oder Erhöhung des Polynomgrades)immer besser durch Funktionen aus diesem Raum approximiert werdenkann.

Wir werden sehen, dass dies im Wesentlichen von der Regularität derLösung, d.h. der maximalen Ordnung ihrer (schwachen)Differenzierbarkeit, abhängt.


Konvergenz bei h-VerfeinerungFolgen von Triangulierungen

Wir betrachten den Fall, dass der FE-Raum ausschließlich durchVerfeinerung des Gitters vergrößert wird, d.h. der Polynomgrad derFE-Funktionen bleibt konstant.

Wir betrachten im Folgenden Triangulierungen aus Dreiecken (mitoffensichtlichen Verallgemeinerungen auf Tetraeder in 3D).

Die Feinheit h einer Triangulierung T h ist definiert als

h := max

diam(K) : K ∈ T h.

Bei der Konvergenzanalyse betrachtet man den Fehler ‖u− uh‖ imGrenzwert h→ 0. Letzterem liegt also eine Folge von FE-Räumen V hzugrunde, jeweils basierend auf einer Triangulierung T h, wobei in dieserFolge die Feinheit h beliebig klein wird.


Konvergenz bei h-VerfeinerungDie Form von Dreiecken

Die Approximierbarkeit einer Funktion in Räumen aus stückweisenPolynomen bezüglich einer Folge T h von Triangulierungen hängtwesentlich davon ab, dass die Dreiecke in dieser Folge nicht beliebig spitzwerden können.

Der Inkreisradius ρK eines Dreiecks K ist definiert durch

ρK := maxρ : K enthält einen Kreis mit Radius ρ

.

Das Verhältnis ρK/diam(K) gibt an, wie spitz das Dreieck K ist.

Definition 28Eine Folge von Triangulierungen T h heißt quasiuniform, falls es eineuntere Schranke σ > 0 gibt, mit

ρKdiamK

≥ σ für alle K ∈ T h für alle h.


Konvergenz bei h-VerfeinerungApproximation durch stückweise lineare Funktionen

Eine obere Schranke für den Fehler der Bestapproximation einerFunktion u ∈ H1(Ω) aus dem Raum V h stetiger stückweise linearerFunktionen bezüglich einer Triangulierung T h von Ω erhält manbeispielsweise durch den Approximationsfehler der InterpolierendenIhu ∈ V h.

Mit dem Céa-Lemma gilt also für den Fehler der Galerkin-Approximation

‖u− uh‖1 ≤C

αinfv∈V h

‖u− v‖1 ≤C

α‖u− Ihu‖1.


Konvergenz bei h-VerfeinerungApproximation durch stückweise lineare Funktionen

Satz 29Sei T h eine quasiuniforme Familie von Triangulierungen einespolygonalen Gebiets Ω ⊂ R2, V h der zugehörige Raum stetigerstückweise linearer Funktionen, Ih die globale Interpolierende sowieu ∈ H2(Ω). Dann existiert eine Konstante C = C(Ω, σ) sodass

‖u− Ihu‖1 ≤ Ch|u|2,‖u− Ihu‖0 ≤ Ch2|u|2.

Hierbei bezeichnet |u|2 die Halbnorm in H2(Ω) gegeben durch

|u|22 =

∫Ω

∑|α|=2

|Dαu|2 dx.


Konvergenz bei h-VerfeinerungApproximation durch stückweise Polynome höheren Grades

Satz 30Sei T h eine quasiuniforme Familie von Triangulierungen einespolygonalen Gebiets Ω ⊂ R2, V h der zugehörige Raum stetigerstückweiser Polynome vom Grad k, Ih die globale Interpolierende sowieu ∈ Hk+1(Ω). Dann existiert eine Konstante C = C(Ω, σ) sodass

‖u− Ihu‖1 ≤ Chk|u|k+1,

‖u− Ihu‖0 ≤ Chk+1|u|k+1.

Hierbei bezeichnet |u|k+1 die Halbnorm in Hk+1(Ω) gegeben durch

|u|2k+1 =

∫Ω

∑|α|=k+1

|Dαu|2 dx.


Konvergenz bei h-VerfeinerungKonvergenzsatz

Aus dem Céa-Lemma und den Interpolationsabschätzungen ergibt sichnun ein Konvergenzresultat.

Satz 31Für die Lösung der RWA (6) auf dem polygonalen Gebiet Ω ⊂ R2 gelteu ∈ Hk+1(Ω). Ist T h eine quasiuniforme Familie von Zerlegungen,V h der zugehörige Raum stetiger stückweiser Polynome vom Grad kund uh ∈ V h die zugehörige Galerkin-Approximation, so gilt

‖u− uh‖1 ≤ C hk |u|k+1.


Inhalt







Numerische Integration

Die Berechnung der bei der Assemblierung auftretenden Integrale ist oftzu aufwändig (bei komplizierten Elementen) oder unmöglich (etwa beinicht geschlossen integrierbaren Koeffizientenfunktionen). Man behilftsich deshalb zur Approximationen der Integrale mit Quadraturformelnder Bauart ∫

Kf(x) dx ≈

m∑i=1

γif(xi) (34)

mit Knoten xi = xKi ∈ K und Gewichten γi = γKi > 0 und erhält somitNäherungen∫

Ωf(x) dx =

∑K∈Th

∫Kf(x) dx ≈

∑K∈Th

m∑i=1

γKi f(xKi ).


Numerische IntegrationQuadraturformeln

Bei affinen Familien können alle Integrale auf solche über einReferenzelement K zurückgeführt werden. Daher genügt es in diesemFall, Quadraturformeln für Referenzelemente zu betrachten.

Quadraturformeln klassifiziert man nach deren Exaktheitsgrad, d.h. demhöchsten Polynomgrad der durch eine Formel noch exakt integriert wird.Eine Quadraturformel für zwei Raumdimensionen besitzt alsoExaktheitsgrad q ∈ N0, falls∫

Kξjηk dξdη =

m∑i=1

γiξji ηki ∀j, k : j + k ≤ q.

Beispiel: Die Gauß-Quadraturformeln für eindimensionale Integralebesitzen bei m Knoten und Gewichten den (maximalen) Exaktheitsgradq = 2m− 1.


Numerische IntegrationKonstruktion von Quadraturformeln

Bei Newton-Cotes-Quadraturformeln werden die Knoten xi vorgegebenund die Gewichte γi so gewählt, dass ein möglichst hoherExaktheitsgrad erreicht wird. Vorausgesetzt die Knoten sind so gewählt,dass sie ein eindeutiges Interpolationspolynom definieren (etwa wie inSatz 15) so sind für Exaktheitsgrad q höchstens n = (q + 1)(q + 2)/2Knoten erforderlich. Bei Integrationsgebieten mit Symmetrien reichenoft auch weniger Knoten aus.

Bei Gauß-Quadraturformeln wird neben den Gewichten auch die Lageder Knoten zur Maximierung des Exaktheitsgrades variiert. Dies führtoft zu wesentlich weniger Knoten als bei Newton-Cotes Formeln gleicherExaktheit. Diese stimmen aber meist nicht mit Knoten fürFreiheitsgrade überein.


Numerische IntegrationNewton-Cotes Formeln für das Referenzdreieck

Bei den folgenden Beispielen für Quadraturformeln bezeichnet

q den Exaktheitsgrad,m die Anzahl Knoten und|K| = 1/2 die Fläche des Referenzdreiecks.

Die Knoten werden sowohl in kartesischen als auch baryzentrischenKoordinaten angegeben. Letztere sind affin invariant, die Formel kann daherauch auf beliebige Dreiecke angewandt werden (nur |K| muss angepasstwerden).

(1) q = 1, m = 1 „Schwerpunktregel“ξ1 = (1

3 ,13) = (1

3 ,13 ,

13)

γ1 = |K|


Numerische IntegrationNewton-Cotes Formeln für das Referenzdreieck

(2) q = 2, m = 3

ξ1 = ( 12 , 0) = (1

2 ,12 , 0) γ1 = γ2 = γ3 = |K|/3

ξ2 = ( 12 ,

12 ) = (0, 1

2 ,12 )

ξ3 = (0, 12 ) = (1

2 , 0,12 )

(3) q = 3, m = 7

ξ1 = ( 13 ,

13 ) = (1

3 ,13 ,

13 ) γ1 = 27|K|/60

ξ2 = (0, 0) = (1, 0, 0) γ2 = γ3 = γ4 = 3|K|/60

ξ3 = (1, 0) = (0, 1, 0)

ξ4 = (0, 1) = (0, 0, 1)

ξ5 = ( 12 ,

12 ) = (0, 1

2 ,12 ) γ5 = γ6 = γ7 = 8|K|/60

ξ6 = (0, 12 ) = (1

2 , 0,12 )

ξ7 = ( 12 , 0) = ( 1

2 ,12 , 0)


Numerische IntegrationEinige Gauß-Formeln für das Referenzdreieck

Hier geben wir zum Vergleich nur d,m und die Lage der Knoten an.

q = 2,m = 3 q = 3,m = 4 q = 4,m = 6

q = 5,m = 7 q = 6,m = 12 q = 7,m = 13


Numerische IntegrationEinige Gauß-Formeln für das Referenzdreieck

Quelle: Die genauen Knoten und Gewichte sind zu finden inG. R. Cowper: Gaussian Quadrature Formulas for Triangles. Int. J. Num.Meth. Eng. 7 (1973) 405–408

oder im Buch von Hughes.

Zur Kontrolle der Exaktheitsgrade (Fehlersuche) ist folgende Formel fürbeliebige Dreiecke K mit Flächeninhalt A hilfreich:∫

Kλ1(x)α1λ2(x)α2λ3(x)α3 dx =

2Aα1!α2!α3!

(α1 + α2 + α3 + 2)!=

2Aα!

(|α|+ 2)!.

Dabei sind λi(x), i = 1, 2, 3, die baryzentrischen Koordinaten von x undαi ∈ N0.


Numerische IntegrationAuswirkungen numerischer Integration

Man kann recht allgemein zeigen, dass bei Verwendung vonQuadraturformeln hinreichend hohen Exaktheitsgrades (in Abhängigkeitvom Polynomgrades im FE-Raum) dieselbe Konvergenzrate in Bezug aufh→ 0 erreicht wird wie bei exakter Auswertung der Integrale.

Im Modellproblem (6) unter Verwendung eines FE-Raumes aus stetigenstückweisen Polynomen vom Grad k reicht hierfür ein Exaktheitsgradvon q = 2k − 2.


Numerische IntegrationBeispiel

Gesucht ist der Wert des Integrals∫ 1−1(2x+ 3)dx.

(1) exakte Integration:∫ 1

−1(2x+ 3)dx =

[x2 + 3x

]1−1

= 1 + 3− 1 + 3 = 6

(2) Gauß-Quadratur: ξ1 = −1/√

3, ξ2 = 1/√

3, γ1 = γ2 = 1/2∫ 1

−1(2x+ 3)dx = 2 ·

(1

2·(− 2√

3+ 3

)+

1

2·(

2√3

+ 3

))= − 2√

3+ 3 +

2√3

+ 3 = 6

(3) Newton-Cotes-Formeln: ξ1 = −1, ξ2 = 1, γ1 = γ2 = 1/2∫ 1

−1(2x+ 3)dx = 2 · (1/2 · (−2 + 3) + 1/2 · (2 + 3)) = 6


Inhalt







Vorgehensweise zur Assemblierung des Galerkin-Gleichungssystems

Ausgangspunkt klassisches Randwertproblem1. Schritt Aufstellen der Variationsformulierung2. Schritt Formulierung des diskreten Gleichungssystems mit Hilfe

der Galerkin-Diskretisierung uh =∑N

j=1 ujφj


Aufgabenstellung

Gesucht ist die Lösung u in Ω = (0, 1) für

−d2u

dx2= f in (0, 1) mit

u(0) = 0 und du

dx(1) = 0.

Die Variationsformulierung lautet: Gesucht ist u ∈ V , so dassa(u, v) = l(v) mit

a(u, v) =

∫ 1

0

du

dx

dv

dxdx und l(v) =

∫ 1

0fvdx.

Die Galerkin-Diskretisierung uh =∑N

j=1 ujφj liefert die diskreteVariationsaufgabe: Gesucht ist uh ∈ V h, so dass

N∑j=1

uj

∫ 1

0

dφidx

dφjdx

dx =

∫ 1

0fφidx ∀i = 1, ..., N.





j=1 ujφj

3. Schritt Bereitstellung der Basisfunktionen und ihrer Gradientenauf dem Referenzelement

Aλ =

[x1

]φi = λi(x) (abhängig von der Ordnung der FE)∇φ



Die baryzentrischen Koordinatenfunktionen

Aus Aλ =

[x1

]mit A =

[x1 x2

1 1

]ergibt sich

λ1(x) =x2 − xx2 − x1

und λ2(x) =x− x1

x2 − x1.

BasisfunktionenFür die Basisfunktionen gilt φ1 = λ1(x) und φ2 = λ2(x).

φk(x) =xk2 − xk

xk2 − xk1=

xk+1 − xxk+1 − xk

=k∆x− x

∆x,

dφkdx

= − 1

∆x

φk+1 =xk − xk1xk2 − xk1

=x− xk

xk+1 − xk=x− (k − 1)∆x

∆x,

dφk+1

dx=

1

∆x





j=1 ujφj


Aλ =

[x1


4. Schritt Berechnung der Elementsteifigkeitsmatrizen und-lastvektoren unter Berücksichtigung der affinen Abbildung



Ak =

[ ∫ek

dφkdx

dφkdx dx

∫ek

dφkdx

dφk+1

dx dx∫ek

dφk+1

dxdφkdx dx

∫ek

dφk+1

dxdφk+1

dx dx

]

bk =

[ ∫ekφkdx∫

ekφk+1dx

]

Ak =

[∫ek

(− 1∆x)(− 1

∆x)dx∫ek

(− 1∆x)( 1

∆x)dx∫ek

( 1∆x)(− 1

∆x)dx∫ek

( 1∆x)( 1

∆x)dx

]=

1

∆x

[1 −1−1 1

]

bk =

[ ∫ek

k∆x−x∆x dx∫

ek


]=

∆x

2

[11

]





j=1 ujφj


Aλ =

[x1



5. Schritt Assemblierung der Steifigkeitsmatrix und des Lastvektorsunter Berücksichtigung der Elementtabelle


Assemblierung des Gleichungssystems

1. Initialisierung von A und b als Null-Matrix bzw. -Vektor2. Addition der Elementmatrizen und -vektoren unter

Berücksichtigung der Elementtabelle

A =1

∆x

1 −1 0 0−1 1 0 00 0 0 00 0 0 0

+1

∆x

0 0 0 00 1 −1 00 −1 1 00 0 0 0

+1

∆x

0 0 0 00 0 0 00 0 1 −10 0 −1 1

=1

∆x

1 −1 0 0−1 2 −1 00 −1 2 −10 0 −1 1

b =∆x

2

1100

+∆x

2

0110

+∆x

2

0011

=∆x

2

1221





j=1 ujφj3. Schritt Bereitstellung der Basisfunktionen und ihrer Gradienten

auf dem ReferenzelementAλ =

[x1




6. Schritt Berücksichtigung der wesentlichen RB (reduziertesGleichungssystem)







[x1





Ergebnis LösungAntje Franke-Börner (TU Freiberg) Numerische Simulation mit finiten Elementen Sommersemester 2013 194



1

∆x

2 −1 0−1 2 −10 −1 1

u2

u3

u4

=∆x

2

221

Lösung

u2 =5

2(∆x)2, u3 = 4(∆x)2, u4 =

9

2(∆x)2







[x1

]φi = λi(x) (abhängig von der Ordnung der FE)∇φi




Ergebnis LösungAntje Franke-Börner (TU Freiberg) Numerische Simulation mit finiten Elementen Sommersemester 2013 196

Aufgabenstellung

Ausgangspunkt Gesucht ist die Lösung u in Ω = [−1, 1]2 ⊂ R2 für

−∆u = 0 in ([−1, 1]× [−1, 1]) mitu(x, y) = g auf ∂Ω = ΓD.

x

y

0

Ω

ΓD

−1 1

−1

1

vgl. Referenzaufgabe 3,Folie 41

Schritt 1 Wie lautet die Variationsformulierung?Schritt 2 Wie lautet die diskrete Variationsaufgabe?


Aufgabenstellung

Ausgangspunkt Gesucht ist die Lösung u in Ω = [−1, 1]2 ⊂ R2 für−∆u = 0 in ([−1, 1]× [−1, 1]) mit

u(x, y) = g auf ∂Ω = ΓD.

Schritt 1 Die Variationsformulierung lautet: Gesucht ist u ∈ V , sodass

a(u, v) = 0 mit

a(u, v) =

∫Ω∇u∇ vdx.

Schritt 2 Die Galerkin-Diskretisierung uh =∑N

j=1 ujφj liefert diediskrete Variationsaufgabe: Gesucht ist uh ∈ V h, so dass

N∑j=1

uj

∫Ω∇φi∇φjdx = 0 ∀i = 1, ..., N.


Referenzelement

Schritt 3 Basisfunktionen und ihre Gradienten auf demReferenzelement

φ1 = ξ, φ2 = η, φ3 = 1− η − ξ[∇φ1 ∇φ2 ∇φ3

]=

[1 0 −10 1 −1

]

ξ

η

1

1

01

2

3


Affine Abbildung

Schritt 4a Affine Abbildung

BK =

[x1 − x3 x2 − x3

y1 − y3 y2 − y3

], |detBK | = 2|K|

B−>K =1

2|K|

[y2 − y3 y3 − y1

x3 − x2 x1 − x3

]

∇φ = B−>K ∇φ

Schritt 4b Elementsteifigkeitsmatrizen

aK(φj , φi) =

∫K∇φj ∇φidx

=

∫KB−>K ∇φjB

−>K ∇φi|detBK |dx


Parametrisierung

Schritt 4 Erstellen Sie Element- und Koordinatenmatrix für folgendeParametrisierung!

x

y

0

Ω

ΓD

−1 1

−1

1

1 2

34

51

2

3

4


Parametrisierung

Schritt 4 Element- und Koordinatenmatrix

ET =

1 2 3 42 3 4 15 5 5 5

, p =

[−1 1 1 −1 0−1 −1 1 1 0

]

x

y

0

Ω

ΓD

−1 1

−1

1

1 2

34

51

2

3

4


Elementsteifigkeitsmatrizen

Schritt 4 Berechnen Sie die Elementsteifigkeitsmatrizen! 4 Elemente- 4 Gruppen.Überlegung: Die Berechnung wievieler Einträge ist nötigbei Berücksichtigung der Symmetrie?


Elementsteifigkeitsmatrizen

Schritt 4 Elementsteifigkeitsmatrizen

a1 = a2 = a3 = a4 = aK =1

2

1 0 −10 1 −1−1 −1 2

Wieso ist die Elementsteifigkeitsmatrix für alle Elementegleich?

Schritt 5 Assemblieren Sie die Steifigkeitsmatrix A!


Steifigkeitsmatrix

Schritt 5 Assemblierung der Steifigkeitsmatrix A ∈ R5×5

A =1

2

1 0 0 0 −10 1 0 0 −10 0 0 0 00 0 0 0 0−1 −1 0 0 2

+1

2

0 0 0 0 00 1 0 0 −10 0 1 0 −10 0 0 0 00 −1 −1 0 2

+1

2

0 0 0 0 00 0 0 0 00 0 1 0 −10 0 0 1 −10 0 −1 −1 2

+1

2

1 0 0 0 −10 0 0 0 00 0 0 0 00 0 0 1 −1−1 0 0 −1 2

=1

2

2 0 0 0 −20 2 0 0 −20 0 2 0 −20 0 0 2 −2−2 −2 −2 −2 8


Reduziertes Gleichungssystem

Schritt 6 Berücksichtigung der wesentlichen Randbedingungen:Referenzproblem 3 auf Folie 41 liefert Randwerte u1, u2, u3

und u4[AWW AWI

AIW AII

] [uwuI

]=

[g0

]mit uW = g

[I O

AIW AII

] [uwuI

]=

[g0

]AIWuW +AIIuI = AIWg +AIIuI = 0

AIIuI = −AIWg

Lösung Berechnen Sie u5 = uI ! Vergleichen Sie mit der exaktenLösung u(x, y) = 2(1+y)

(3+x)2+(1+y)2für (0, 0)!


Lösung

Ergebnis

g =

00

0.20.5

, AIW =[−1 −1 −1 −1

], AII = 4

4u5 = 0.2 + 0.5 = 0.7

u5 = 0.175

Vergleich mit exakter Lösung: u(0, 0) = 0.2


Inhalt







Nédélec-ElementeLineares Tetraederelement (3D)

1

2 34 5

6

(0,0,0) (1,0,0)

(0,1,0)

(0,0,1)

AnwendungSimulation vektorwertiger Größen, z.B. EM Felder3D MT, 3D TEM, 3D HEM


Inhalt







Nédélec-Elemente

Bisher: Differentialgleichungen mit Gradientenoperator imVariationsintegral. Diese führten auf Variationsraum H1(Ω).

In vielen Gleichungen der mathematischen Physik treten neben demGradienten auch Rotations- und Divergenzoperator auf.

Beispiel: Zeitharmonische Maxwell-Gleichung (16) nach Elimination desMagnetfeldes:

∇×(∇×E) + iωµσE = −iωµJe . (35)


Partielle Integration für ∇×

Zur FE-Diskretisierung von (35) gehen wir analog vor wie bei derPoisson-Gleichung und leiten zunächst eine schwache Formulierung her.Hierfür wird eine Formel zur partiellen Integration für Terme benötigt,die den Rotationsoperator beinhalten.

Für hinreichend glatte Vektorfelder u und v auf hinreichend regulärenbeschränkten Gebieten Ω ⊂ R3 gilt∫

Ωv · (∇×u) dx =

∫Ωu · (∇×v) dx+

∫∂Ω

(n× u) · vt ds . (36)

mitvt := v − (v · n)n .


Tangentialkomponenten

Bezeichnungen: vt : Tangentialkomponente von v,n× v : „rotierte“ Tangentialkomponente von v.

Es gilt:

(i) vt = n× (v × n)(ii) n× v = n× vt(iii) (n× u) · vt = (n× u) · v(iv) (n× u) · vt = −ut · (n× v)

Diese Beziehungen folgen unmittelbar aus folgenden beiden Identitätenfür das Vektorprodukt im Euklidischen Anschauungsraum: Für allea, b, c ∈ R3 gilt

a× (b× c) = (a · c)b− (a · b)c,a · (b× c) = b · (c× a) = c · (a× b).


FunktionenräumeQuadratisch integrierbare Rotation

Damit die Integrale in (36) definiert sind, müssen die Vektorfelder u undv sowie ∇×u und ∇×v quadratisch integrierbar sein. Dies führt aufdie Funktionenräume

L2(Ω) := u : Ω→ C3 :

∫Ω|u|2 dx <∞ = L2(Ω)3,

H(rot; Ω) := u ∈ L2(Ω) : ∇×u ∈ L2(Ω).Man nennt ein Vektorfeld v ∈ L2(Ω) schwache Rotation vonu ∈ L2(Ω), falls für jedes auf ∂Ω verschwindende differenzierbareVektorfeld φ gilt ∫

Ωu · (∇×φ) dx =

∫Ωv · φ dx.

(36) legt auch nahe, welche wesentliche Randbedingung in H(rot; Ω) zustellen ist: man definiert

H0(rot; Ω) := u ∈H(rot; Ω) : n× u = 0 längs ∂Ω.Antje Franke-Börner (TU Freiberg) Numerische Simulation mit finiten Elementen Sommersemester 2013 214

Konformität in H(rot; Ω)

Das Zusammenstückeln von Funktionen u ∈H(rot; Ω) aus Funktionenauf Teilgebieten von Ω erfordert an den Übergängen die Stetigkeit derTangentialkomponenten

n× u.

Satz 32Seien Ω1,Ω2 zwei nichtüberlappende Lipschitz-Gebiete mitgemeinsamem Randstück Γ von positivem Maß, d.h. Ω1 ∩Ω2 = Γ. Seienferner ui ∈ L2(Ωi), (i = 1, 2) und bezeichne

Ω := Ω1 ∪ Ω2 ∪ Γ, u(x) :=

u1(x), x ∈ Ω1,

u2(x), x ∈ Ω2.

Gilt ui ∈H(rot; Ωi) (i = 1, 2) und n× u1 = n× u2 längs Γ, so istu ∈H(rot; Ω).


Randbedingungen

Dirichlet RB:n×E = ED. (37)

Längs der Ränder von perfekten elektrischen Leitern gilt z.B.n×E = 0.Neumann RB: Vorgabe einer Flächenstromdichte n×H = JeF ,was nach der Faradayschen Gleichung auf

n× (∇×E) = −iωµJeF (38)

führt.Gemischte RB: Flächenstromdichte JeF bei Impedanz Z(Impedanz-RB):

n× (∇×E)− iωµZ−1Et = −iωµJeF . (39)


Variationsformulierung

Anwendung der partiellen Integrationsformel (36) auf dieMaxwell-Gleichung (35) zweiter Ordnung nach Multiplikation mit einerTestfunktion v führt auf

−iωµ∫

ΩJe · v dx =

∫Ω

[(∇×E) · (∇×v)− k2E · v

]dx

+

∫∂Ω

[n× (∇×E)] · vt ds.

Besteht der Rand ∂Ω aus drei Teilstücken ΓD, ΓN und ΓG, auf denenjeweils die Dirichlet-RB (37), die Neumann-RB (38) bzw. die gemischteRB (39) gestellt sind, so bieten sich als Ansatz- bzw. Testräume an:

S := E ∈H(rot; Ω) : n×E = ED längs ΓD (40a)V := v ∈H(rot; Ω) : n× v = 0 längs ΓD (40b)


Variationsformulierung

Die Variationsformulierung dieser RWA lautet somit:

Bestimme E ∈ S sodass

a(E,v) = `(v) ∀v ∈ V (41)

mit

a(E,v) =

∫Ω

[(∇×E) · (∇×v)− k2E · v

]dx+ iωµ

∫ΓG

Z−1Et · vt ds,

`(v) = −iωµ∫

ΩJe · v dx+ iωµ

∫ΓN∪ΓG

JeF · vt ds.

Die Bilinearform im zweiten Argument bzw. die Linearform sind hier (imKomplexen) antilinear.


FunktionenräumeQuadratisch integrierbare Divergenz

Man nennt eine Funktion v ∈ L2(Ω) schwache Divergenz vonu ∈ L2(Ω), falls für jede auf ∂Ω verschwindende differenzierbareFunktion φ gilt ∫

Ωu · ∇φdx = −

∫Ωvφ dx.

Der Raum aller Vektorfelder mit quadratisch integrierbarer Divergenzwird bezeichnet mit

H(div; Ω) := u ∈ L2(Ω) : ∇·u ∈ L2(Ω).


FunktionenräumeQuadratisch integrierbare Divergenz

Für hinreichend glatte Funktionen u und v auf hinreichend regulärenbeschränkten Gebieten Ω gilt∫

Ωv(∇·u) dx = −

∫Ωu · ∇ v dx+

∫∂Ωvun ds

mit un := n · u.Anhand dieser Greenschen Formel erkennt man auch, dass fürH(div; Ω)-Funktionen als Randwerte die Normalkomponente unvorzuschreiben ist. Daher setzt man

H0(div; Ω) := u ∈H(div; Ω) : un = 0 längs ∂Ω.


FunktionenräumeInnenprodukte

Die Räume H(rot; Ω) und H(div; Ω) sind Hilbert-Räume ausgestattetmit den Innenprodukten

(u,v)H(rot;Ω) :=

∫Ω

(u · v + (∇×u) · (∇×v)) dx

bzw.(u,v)H(div;Ω) :=

∫Ω

(u · v + (∇·u) · (∇·v)) dx.

Die zugehörigen Normen werden mit ‖ · ‖H(rot;Ω) bzw. ‖ · ‖H(div;Ω)

bezeichnet.


Skalarpotentiale

Es ist oft nützlich, Vektorfelder durch Potentiale in skalarer odervektorieller Form darzustellen. Dieser Abschnitt befasst sich damit, wanndies möglich ist.

Für glatte Funktionen ist bekannt, dass ein Vektorfeld durch ein skalaresPotential darstellbar ist, wenn seine Rotation verschwindet.

Satz 33Sei Ω ⊂ R3 ein einfach zusammenhängendes Lipschitz-Gebiet sowieu ∈ L2(Ω). Dann gilt ∇×u = 0 genau dann, wenn ein skalaresPotential φ ∈ H1(Ω) existiert mit u = ∇φ. Ferner ist φ bis auf eineadditive Konstante eindeutig bestimmt.


Vektorpotentiale

Um der Frage nachzugehen, wann ein Vektorpotential existiert, müssenwir das Grundgebiet Ω zunächst genauer beschreiben.

Wir nehmen an, Ω sei ein beschränktes, zusammenhängendesLipschitz-Gebiet und bezeichnen mit ∂ΩjJj=0 die zusammenhängendenKomponenten von ∂Ω. ∂Ω0 bezeichne dabei diejenige Komponente von∂Ω, welche den Rand der unbeschränkten Komponente von R3 \ Ωbildet.

Mit ΩjJj=0 bezeichnen wir das Teilgebiet von R3 \ Ω mit Rand ∂Ωj ,d.h. Ω0 ist unbeschränkt.


Vektorpotentiale

Satz 34Zu jeder Funktion u ∈H(div; Ω) mit

∇·u = 0 in Ω sowie∫∂Ωj

un ds = 0, j = 0, . . . , J,

existiert ein Vektorpotential A ∈ H1(Ω)3 mit

u = ∇×A und ∇·A = 0 in Ω.


Helmholtz-Zerlegung

Ebenfalls häufig von Nutzen ist es, ein Vektorfeld als Summe einesSkalar- und eines Vektorpotentials zu schreiben. Eine solche Darstellungheißt Helmholtz-Zerlegung.Wir benötigen für das Weitere noch folgende Notation bezüglich Ω,welches hier nicht als einfach zusammenhängend vorausgesetzt wird.Wir nehmen an, es gäbe L offene, zusammenhängende, in Ω enthalteneFlächen Σ`L`=1, genannt innere Schnitte, sodass für alle ` gilt:1. Jede Fläche Σ` ist der offene Teil einer glatten Fläche;2. ∂Σ` ⊂ ∂Ω;3. Σ` ∩ Σm = ∅ falls ` 6= m;4. die Menge Ω0 := Ω \ ∪L`=1Σ` ist einfach zusammenhängend und

pseudo-Lipschitz.


Helmholtz-Zerlegungde-Rham-Diagramm

Angefangen bei Unterräumen von H1(Ω) bilden die behandeltenDifferentialoperatoren ab in die Folge von Räumen

H1(Ω)/R ∇−→H(rot; Ω)∇×−→H(div; Ω)

∇·−→ L2(Ω) (42)

bzw. mit Nullrandbedingungen

H10 (Ω)

∇−→H0(rot; Ω)∇×−→H0(div; Ω)

∇·−→ L2(Ω)/R. (43)

Diese Darstellung wird als de-Rham-Diagramm bezeichnet.

Satz 35In den Diagrammen (42) und (43) ist der Bildraum des jeweilsvorangehenden Operators enthalten im Nullraum des nachfolgendenOperators. Der Bildraum ist jeweils ein abgeschlossener Unterraum desentsprechenden Nullraumes mit endlicher Kodimension6.

6Dimension des orthogonalen Komplements eines abgeschlossenen Unterraumes


Helmholtz-Zerlegung

Wir betrachten speziell die Folge (43). Da ∇H10 (Ω) ein abgeschlossener

Unterraum des Nullraums des Rotations-Operators ist, können wir nachdem Projektionssatz diesen Nullraum schreiben als

N (rot) = ∇H10 (Ω)⊕ (∇H1

0 (Ω))⊥,

wobei das orthogonale Komplement bezüglich desH(rot; Ω)-Innenprodukts zu vestehen ist.Sei nun u ∈ N (rot) ⊂H0(rot; Ω), und zwar u ∈ (∇H1

0 (Ω))⊥. Danngilt ∇×u = 0 in Ω, n× u = 0 längs ∂Ω, sowie

(u,∇φ) = 0 ∀φ ∈ H10 (Ω).

Partielle Integration liefert nun ∇·u = 0 in Ω. Diese Argumentationlässt sich umkehren und man erhält

(∇H10 (Ω))⊥ = u ∈H0(rot; Ω) : ∇×u = 0,∇·u = 0 in Ω.


Helmholtz-Zerlegung

Dieser Raum heißt normaler Kohomologieraum von Ω und wird mitKn(Ω) bezeichnet.

Analoge Betrachtung zeigt: das orthogonale Komplement von∇×H0(rot; Ω) im Nullraum des Divergenzoperators ist gegeben durchden tangentialen Kohomologieraum Kt(Ω) definiert durch

Kt(Ω) := u ∈H0(div; Ω) : ∇·u = 0,∇×u = 0 in Ω.

Satz 35 besagt, dass die Dimensionen von Kn(Ω) und Kt(Ω) endlichsind.


Helmholtz-ZerlegungBezug zur Elektrostatik

Diese Kohomologieräume treten auf natürliche Weise z.B. in derElektrostatik (und Magnetostatik) auf.

Gesucht sei eine statische/stationäre Lösung der Maxwell-Gleichungen ineinem beschränkten Gebiet mit metallischem Rand und ε = 1, σ = 0(Vakuum) sowie J = 0 (kein Quellstrom).

Es gelten also

∇×E = 0 in Ω (statisches Feld),∇·E = 0 in Ω (keine Quellen),n×E = 0 längs ∂Ω (metallischer Rand).

Damit liegt E in Kn(Ω).Analog liegt die magnetische Flussdichte B eines quellenfreienstatischen Feldes in Kt(Ω).


Helmholtz-ZerlegungZusammenfassung

Satz 36

1. Ist u ∈H0(rot; Ω) mit ∇×u = 0 in Ω, so existieren eineindeutiges Skalarpotential φ ∈ H1

0 (Ω) sowie eine Funktionfn ∈ Kn(Ω) mit

u = ∇φ+ fn.

2. Ist u ∈H0(div; Ω) mit ∇·u = 0 in Ω, so existiert einVektorpotential A ∈H0(rot; Ω) und eine Funktion f t ∈ Kt(Ω) mit

u = ∇×A+ f t.

Das Vektorpotential wird eindeutig, wenn wir verlangen, dass∇·A = 0 in Ω sowie

∫∂Ωj

n ·A ds = 0 für j = 0, 1, . . . , J .



Satz 37Jede Funktion in L2(Ω) besitzt die Zerlegung

u = ∇ p+ fn +∇×A

mit eindeutigem p ∈ H10 (Ω), fn ∈ Kn(Ω) sowie

A ∈ w ∈H(rot; Ω) :∇×w = 0 in Ω,n ·w = 0 längs ∂Ω,∫Σ`

n ·w ds = 0, ` = 1, . . . , L.



Bemerkungen 38

1. Eine analoge Zerlegung gilt mit p ∈ H1(Ω), f t ∈ Kt(Ω) sowieA ∈H0(rot; Ω).

2. Es ist dimKn(Ω) = J sowie dimKt(Ω) = L.3. Im Fall, dass Ω einfach zusammenhängend ist, ist L = 0.


Inhalt







Divergenz- und Rotationskonforme ElementeGebietstransformationen

Seien

K, K ⊂ R3 beschränkte Gebiete,FK : K → K, x 7→ x = FK(x) stetig differenzierbare Bijektion,detF ′K ohne Vorzeichenwechsel.

Einer skalaren Funktion u ∈ H1(K) wird durch

u = u FK (44)

eine skalare Funktion u ∈ H1(K) zugeordnet.

Seien ∇, ∇ die Gradienten bezüglich x, x, dann gilt nach derKettenregel

∇u = (F ′K)−> ∇u. (45)


Divergenz- und Rotationskonforme ElementeGebietstransformationen - Vektorfelder

Die Transformation von Vektorfeldern erfordert zusätzlich eineAnpassung der abhängigen Variablen:

Für u, u ∈ H1(Ω) wie oben gelten

∇u ∈H(rot; K) und ∇u ∈H(rot;K)

(vgl. (42)).

Damit die Gradienten durch die gesuchte Transformation ebenfallsaufeinander abgebildet werden, müssen zwei Vektorfelder

v ∈H(rot;K) und v ∈H(rot; K)

transformiert werden gemäß

v FK = (F ′K)−>v. (46)


Divergenz- und Rotationskonforme ElementeGebietstransformationen - Rotationsoperator

Lemma 39Für die Funktionen v ∈H(rot;K) und v ∈H(rot; K) gelte dieBeziehung (46) mit einer stetig differenzierbaren BijektionFK : K → K. Bezeichnet [∇×v] die 3× 3-Matrix

[∇×v]i,j =∂vi∂xj− ∂vj∂xi

, (i, j = 1, 2, 3),

so gilt[∇×v] FK = (F ′K)−> [∇ × v] (F ′K)−1.

Mit v ∈H(rot; K) ist somit auch v ∈H(rot;K).


Divergenz- und Rotationskonforme ElementeGebietstransformationen - Rotationsoperator

Korollar 40Unter den Voraussetzungen von Lemma 39 bestehe zwischenv ∈H(rot;K) und v ∈H(rot; K) die Beziehung (46). Dann istv ∈H(rot;K) und es gilt

[∇×v] FK =1

detF ′KF ′K [∇ × v].


Divergenz- und Rotationskonforme ElementeGebietstransformationen - Divergenzoperator

Analoge Transformation für Divergenz: nach (42) gilt mitv ∈H(rot; K) auch ∇ × v ∈H(div; K).

Eine Funktion w ∈H(div; K) muss also so transformiert werden, dass

w FK =1

detF ′KF ′K w. (47)

Lemma 41Seien w,w zwei differenzierbare Funktionen, welche (47) erfüllen miteiner stetig differenzierbaren Bijektion FK : K → K. Dann gilt

(∇·w) FK =1

detF ′K(∇ · w).

Insbesondere ist mit w ∈H(div; K) auch w ∈H(div;K).


Divergenz- und Rotationskonforme ElementeGebietstransformationen: Weitere nützliche Fakten

Ist n äußere Einheitsnormale an K und x ∈ ∂K, so ist n definiert durch

n(FK(x)) :=(F ′K)−>n

‖(F ′K)−>n‖(x),

äußere Einheitsnormale an K.

Ist t ein Einheitstangentenvektor an ∂K im Punkt x ∈ ∂K, so ist

t(FK(x)) :=F ′K t

‖F ′K t‖(x)

ein Einheitstangentenvektor an ∂K im Punkt FK(x).


Divergenz- und Rotationskonforme ElementeGebietstransformationen: Oberflächenintegrale

Gilt für u, u (44), so gilt mit Jn(x) := |detF ′K | ‖(F ′K)−>n‖∫∂K

u ds =

∫∂K

uJn ds. (48)

Gilt für w, w (47) und für φ, φ (44), so gilt∫∂K

(n ·w)φds = sgn(detF ′K)

∫∂K

(n · w)φ ds. (49)

Gilt für u, u (46) und ebenso für v, v, so gilt∫∂K

(n× u) · vt ds = sgn(detF ′K)

∫∂K

(n× u) · vt ds. (50)


Rotationskonforme Elemente

Whitney, 1957Nédélec, 1980„Kantenelemente“, (edge elements)


Rotationskonforme ElementePolynomraum

Sk :=p ∈ P3

k : x · p(x) = 0

P3k meint den Raum der homogenene Polynome p(x) (x ∈ R3) vom

exakten Grad k. Damit ist klar: p ∈ P3k ⇒ x · p ∈ Pk+1.

Umgekehrt kann jedes q ∈ Pk+1 geschrieben werden als q = x · p miteinem p ∈ P3

k .

Somit folgt:

dim Sk = 3 dim Pk − dim Pk+1 = k(k + 2).


Rotationskonforme Elemente. . . noch ein Polynomraum

Rk := P3k−1 ⊕Sk.

Es gilt dim Rk = 3 dim Pk−1 + dim Sk = 12k(k + 2)(k + 3).

Lemma 42P3k = Rk ⊕∇ Pk+1.

Lemma 43Gilt ∇×u = 0 für ein u ∈ Rk, so gilt u = ∇ p für ein p ∈Pk.


Rotationskonforme Elemente2D Fall

Bemerkung 44Analoger Polynomraum für zweidimensionale Gebiete bzw. Oberflächen:für ein Dreieck f in der (x1, x2)-Ebene:

Sk(f) :=p ∈ P2

k : x · p(x) = 0,

Rk(f) := Pk−1(f)2 ⊕Sk(f)

Lemma 43 gilt dann bezüglich des skalaren Rotationsoperators∇×u = ∂u2

∂x1− ∂u1

∂x2.



Definition 45Ein rotationskonformes finites Element (K, PK ,ΨK) der Ordnungk ∈ N auf dem Referenzelement wird definiert durch(a) K sei das Referenztetraeder,(b) der Funktionenraum sei gegeben durch PK = Rk,(c) die Freiheitsgrade ΨK = Me ∪Mf ∪MK seien gegeben durch die

Momente

Me(u) =

∫e

(u · t)q ds : q ∈Pk−1(e), e Kante von K, (51a)

Mf (u) =

1

|f |

∫f

u · q ds : q ∈Pk−2(f)3, n · q = 0, f Fläche von K,

(51b)

MK(u) =

∫K

u · q dx : q ∈P3k−3(K)

. (51c)


Rotationskonforme ElementeTransformation auf Gebietselement

Mit u ∈ Rk auf K assoziieren wir gemäß (46) im speziellen Fall einer affinenAbbildung FK : K → K

u FK = B−>K u. (52)

Gemäß Korollar 40 besteht zwischen den Rotationen der beiden Vektorfelder dieBeziehung

∇×u =1

detBKBK(∇ × u).

Ist t ein Einheitsvektor in Richtung der Kanten e von K, so ist

t :=BK t

‖BK t‖(53)

ein Einheitsvektor in Richtung des Bildes e der Kante e unter FK .

Lemma 46Rk ist invariant unter der Transformation (52).



Definition 47Das rotationskonforme finite Element (K,PK ,ΨK) auf einem beliebigenTetraeder K ⊂ R3 ist definiert durch(a) K sei ein beliebiges (nichtentartetes) Tetraeder.

(b) der Funktionenraum sei gegeben durch PK = Rk.

(c) Bezeichnen t und n Einheitsvektoren längs einer beliebigen Kante e bzw.normal zu einer beliebigen Fläche f von K, E (K) und F (K) die Mengeder Kanten bzw. Flächen von K, so seien die Freiheitsgrade gegebendurch

Me(u) =

∫e

(u · t)q ds : q ∈Pk−1(e), e ∈ E (K)

, (54a)

Mf (u) = 1

|f |

∫f

u · q ds : q ∈ BKPk−2(f)3,n · q = 0, f ∈ F (K),

(54b)

MK(u) =

∫K

u · q dx : q FK ∈ BK/(detBK)P3k−3(K)

. (54c)



Lemma 48Sei detBK > 0 und es mögen die Tangentenvektoren an die Kanten von Kund K in der Beziehung (53) stehen. Dann sind unter der Voraussetzung(52) die Freiheitsgrade (51) von u auf K identisch mit den Freiheitsgraden(54) von u auf K.

Konformität: Nach Lemma 32 zu zeigen: n× u stetig überElementgrenzen. Da die Kanten- und Flächenfreiheitsgrade an gemeinsamenFlächen übereinstimmen, sind die Freiheitsgrade der Differenz dort Null.

Lemma 49Verschwinden für eine Funktion u ∈ Rk alle zu einer Fläche f oder derenKanten e gehörenden Freiheitsgrade aus (54), so gilt n× u = 0 auf f .

Lemma 50Verschwinden für eine Funktion u ∈ Rk alle Freiheitsgrade (54), so istu = 0.


Hausaufgabe

1. Schritt Finden Sie die Variationsformulierung für folgendesklassisches Randwertproblem

∇×(µ−1∇×E) + (iωσ − ω2ε)E = 0 in Ω

n×E = 0 auf Γ⊥

n×H = 0 auf Γ||

µ−1∇×E = −iωHn(x, y, z) auf Γtop ∪ Γbottom!

2. Schritt Wie lautet die diskrete Variationsaufgabe unterBerücksichtigung von

uh(x) =

n∑i=1

lhi (uh)φhi (x)?


Hausaufgabe - Lösung

1. Schritt Die Variationsformulierung des Randwertproblems lautet:Finde E ∈ S , so dass

a(E,v) = l(v) ∀v ∈ V mit

a(E,v) =

∫Ω

(µ−1∇×E · ∇×v + (iωσ − ω2ε)E · v

)dx,

l(v) = iω

∫Γtop∪Γbottom

(n×Hn) · vtds

S := E ∈H(rot; Ω) : n×E = 0 auf Γ⊥ undV := v ∈H(rot; Ω) : n× v = 0 auf Γ⊥.


Hausaufgabe - Lösung

2. Schritt Das diskrete Gleichungssystem lautet:

Au = b mit

Ai,j = a(φi,φj), bi = l(φi), uj = lhj (Eh) und

a(φi,φj) =

∫Ωµ−1∇×φi · ∇×φjdx

+

∫Ω

(iωσ − ω2ε)φi · φjdx,

l(φi) = iω

∫Γtop∪Γbottom

n×Hn · φi,tds


Rotationskonforme ElementeDas rotationskonforme Element niedrigster Ordnung

Für den Fall k = 1 zeigt man leicht

R1 =u(x) = a+ b× x : a, b ∈ C3

.

Hier treten nur Kantenmomente auf, daher auch der Name„Kantenelemente“. Die 6 Freiheitsgrade (d.h. die Komponenten derVektoren a und b) werden durch die 6 Kantenmomente bestimmt.

Durch Nachrechnen bestätigt man: die nodale Basisfunktion mitKantenmoment Eins längs der Kante von Knoten i nach Knoten j istgegeben durch

φi,j = λi∇λj − λj ∇λi, i, j = 1, . . . 4, i < j.

wobei λi die baryzentrischen Koordinaten im Tetraeder bezeichnen.



Da die Mittelpunktregel exakt für Polynome vom Grad ≤ 1 ist, kannman durch Auswertung von φi,j an den Kantenmittelpunkten zeigen,dass ∫

ei,j

φi,j · t ds = 1.

Ferner gilt∇×φi,j = 2∇λi ×∇λj .


Rotationskonforme ElementeMasse- und Steifigkeitsmatrix auf dem Referenztetraeder

Massematrix M

Mi,j = (iωσ − ω2ε)

∫Kφi · φjdx

M =(iωσ − ω2ε)

120

4 1 0 −1 0 01 4 0 1 0 00 0 10 0 5 5−1 1 0 4 0 00 0 5 0 10 50 0 5 0 5 10


Rotationskonforme ElementeMasse- und Steifigkeitsmatrix auf dem Referenztetraeder

Steifigkeitsmatrix K

Ki,j = µ−1

∫K∇× φi · ∇× φjdx

K =2

3µ−1

1 0 −1 0 1 00 1 −1 0 0 1−1 −1 2 0 −1 −10 0 0 1 −1 11 0 −1 −1 2 −10 1 −1 1 −1 2


Rotationskonforme ElementeAffine Abbildung

Für jedes Tetraeder K ∈ T h ist die affine Abbildung FK bestimmtdurch die Abbildungsvorschriften

(1, 0, 0) 7→ (x1, y1, z1),

(0, 1, 0) 7→ (x2, y2, z2),

(0, 0, 1) 7→ (x3, y3, z3),

(0, 0, 0) 7→ (x4, y4, z4), d.h.xyz

=

x1 − x4 x2 − x4 x3 − x4

y1 − y4 y2 − y4 y3 − y4

z1 − z4 z2 − z4 z3 − z4

︸︷︷︸

BK

xyz

+

x4

y4

z4

︸︷︷︸bK

.



Theorem 51Das in Definition 47 definierte finite Element ist H(rot)-konform undunisolvent.

Ist Th ein zulässiges Tetraedergitter auf Ω, so definieren wir den globalenFE-Raum durch

V h := u ∈H(rot; Ω) : u|K ∈ Rk ∀K ∈ Th. (55)

Die lokale Interpolierende IK = I(R)K auf einem Element K ordnet einer

Funktion u ∈H(rot;K), für welche alle Freiheitsgrade (54) definiert sind,eine Funktion IKu ∈ Rk zu, welche charakterisiert ist durch

Me(u− IKu) = Mf (u− IKu) = MK(u− IKu) = 0.

Die globale Interpolierende Ih = I(R)h ist dann wie üblich definiert durch

Ihu = IKu für alle K ∈ Th.



Da dies für die Kantenmomente nicht der Fall ist, ist die Interpolierendenicht für alle Vektorfelder u ∈H(rot; Ω) definiert.Im Buch [Monk (2003)] werden folgende hinreichende Bedingungenangegeben:

u ∈ H12

+δ(K)3, ∇×u ∈ Lp(K)3 ∀K ∈ Th, (δ > 0, p > 2).

Von [Amrouche et al. (1998)] stammt die Abschwächung auf

u ∈ Lp(K)3, ∇×u ∈ Lp(K)3, n× u ∈ Lp(∂K)3, (p > 2).


Rotationskonforme ElementeEigenschaften der Interpolierenden

Lemma 52Sei W h der divergenzkonforme FE-Raum

Wh := u ∈H(div; Ω) : u|K ∈ Dk ∀K ∈ Th,

wobei Dk := P3k−1 ⊕ Pk−1x und V h der rotationskonforme FE-Raum

V h := u ∈H(rot; Ω) : u|K ∈ Rk ∀K ∈ Th.

Dann gilt ∇×V h ⊂ W h und, sofern für u sowohl I(R)h u als auch

I(D)h (∇×u) definiert sind, gilt auch

I(D)h (∇×u) = ∇×(I

(R)h u).


Rotationskonforme ElementeEigenschaften der Interpolierenden

Theorem 53Sei Th ein reguläres Tetraedergitter auf Ω. Dann gelten1. Ist u ∈ Hs(Ω)3 sowie ∇×u ∈ Hs(Ω) mit 1

2 + δ ≤ s ≤ k, δ > 0, sogilt

‖u− I(R)h u‖L2(Ω)3 + ‖∇×(u− I(R)

h u)‖L2(Ω)3

≤ Chs(‖u‖Hs(Ω)3 + ‖∇×u‖Hs(Ω)3).

2. Ist u ∈ H 12

+δ(K)3 mit 0 < δ ≤ 12 sowie ∇×u|K ∈ Dk, so gilt

‖u−I(R)h u‖L2(K)3 ≤ C

(h

12

+δ

K ‖u‖H

12+δ(K)3

+ hK‖∇×u‖L2(K)3

).


Rotationskonforme ElementeEinbeziehung wesentlicher Randbedingungen

Bei der Variationsformulierung (41) wurden auf einem Teil ΓD desGebietsrandes die wesentliche Randbedingung n× u = 0 (vgl. (40))gestellt. Ein weiterer Vorteil des Nédélecschen FE-Raumes ist, dass dieErfüllung dieser homogenen wesentlichen Randbedingung durchNullsetzen der Freiheitsgrade auf ΓD erreicht wird.

Lemma 54Liegt u im Variationsraum V für die in (41) definierte RWA, so gilt diesauch für I(R)

h u.

Demnach ist ein FE-Unterraum V h ⊂ V , wobei V wie in (40) definiertist, gegeben durch

V h := u ∈H(rot; Ω) : u|K ∈ Rk ∀K ∈ Th,n× u = 0 längs ΓD.


Rotationskonforme ElementeNédélec-Elemente zweiter Ordnung

Im Fall des (20-dimensionalen) Polynomraumes R2 kann man mit Hilfeder linearen Vektorpolynome

p1 =

0−x3

x2

, p2 =

x3

0−x1

, p3 =

−x2

x1

0

.und der Standardbasis ej3j=1 im kartesischen Koordinatensystemfolgende Basis angeben



v1 = e1, v2 = e2, v3 = e3,

v4 = x1e1, v5 = x1e2, v6 = x1e3,

v7 = x2e1, v8 = x2e2, v9 = x2e3,

v10 = x3e1, v11 = x3e2, v12 = x3e3,

v13 = x1p1, v14 = x2p1, v15 = x3p1,

v16 = x1p2, v17 = x2p2, v18 = x3p2,

v19 = x1p3, v20 = x2p3.

Jede Funktion u ∈ R2 kann somit dargestellt werden alsLinearkombination

u =

20∑j=1

uj vj (56)

mit eindeutig bestimmten Koeffizienten uj20j=1.



Zur Bestimmung der nodalen Basis von R2 bezüglich der Freiheitsgrade(51) auf dem Referenzelement bezeichnen wir die Ecken desReferenztetraeders mit

a1 =

100

, a2 =

010

, a3 =

001

, a4 =

000

.Jeder Kante e = [ai, aj ] wird der Tangenteneinheitsvektor

ti,j :=aj − ai‖aj − ai‖

, i, j = 1, . . . , 4, i < j.

zugeordnet. Somit sind die Kantenfreiheitsgrade des Referenzelementesgegeben durch ∫ aj

ai

(u · ti,j) q ds, q ∈P1(ei,j).



Da (u · ti,j)q ∈P3(ei,j) kann dieses Integral mit einer Gauß-LegendreQuadraturformel mit zwei Knoten exakt ausgewertet werden:∫ aj

ai

(u · ti,j)q ds = 12

(u(ξ

(1)

i,j )q(s(1)i,j ) + u(ξ

(2)

i,j )q(s(2)i,j ))· (aj − ai).

Dabei bezeichnen s(`)i,j (` = 1, 2) die Gauß-Knoten in der Bogenlänge auf

der Strecke von ai nach aj sowie ξ(`)i,j (` = 1, 2) die entsprechenden

dreidimensionalen Koordinaten der Gauß-Knoten auf [ai, aj ].

Wir wählen nun das (lineare) Polynom q so, dass

q(s(1)i,j ) = 1, q(s

(2)i,j ) = 0

bzw. umgekehrt und erhalten so auf jeder Kante die beidenFreiheitsgrade

u(ξ(1)

i,j ) · (aj − ai) and u(ξ(2)

i,j ) · (aj − ai).Antje Franke-Börner (TU Freiberg) Numerische Simulation mit finiten Elementen Sommersemester 2013 265


Für jede Seitenfläche f mit Ecken ai, aj , ak (i < j < k) wählen wirddie (nicht normierten) Vektoren

t(f)1 = aj − ai, t

(f)2 = ak − ai,

und erhalten somit die Freiheitsgrade1

|f |

∫fu · t(f)

1 ds und 1

|f |

∫fu · t(f)

2 ds

Da u höchstens vom Grad zwei ist, können diese Flächenintegrale durcheine Quadraturformel mit Knoten an den Seitenmittelpunkten berechnetwerden:

1

|f |

∫fu · t(f)

1 ds =1

3

(u

(ai + aj

2

)+ u

(ai + ak

2

)

+ u

(aj + ak

2

))· t(f)` , ` = 1, 2.



Mit diesen Freiheitsgraden können wir nun die nodale Basis von R2

konstruieren. Die Basisfunktion mit Freiheitsgrad Eins am erstenQuadraturknoten der Kante [a1, a2] erhalten wir beispielsweise, indemwir die Koeffizienten uj in der Darstellung (56) bestimmen aus denGleichungen

u(ξ(1)

1,2) · (a2 − a1) = 1,

u(ξ(2)

1,2) · (a2 − a1) = 0,

u(ξ(k)

i,j ) · (aj − ai) = 0, i < j, (i, j) 6= (1, 2), 1 ≤ i, j ≤ 4,

1

|fk|

∫fk

u · t(fk)` ds = 0, 1 ≤ k ≤ 4, 1 ≤ ` ≤ 2.

Hier bezeichnet fk die k-te Seitenfläche des Referenztetraeders K.


Rotationskonforme ElementeKonsistente Orientierung

Zwei mögliche Orientierungen bei Seitenflächen und Kanten, diesemüssen zwischen angrenzenden Tetraedern konsistent sein, damit jedeBasisfunktion in allen zu ihrem Träger gehörenden Tetraedern dasselbeVorzeichen besitzt.

Folgendes Vorgehen stellt dies sicher:Globale Nummerierung aiNhi=1 aller Tetraederecken im Gitter Th

In jedem Tetraeder K ∈ Th Nummerierung der Eckenai1 ,ai2 ,ai3 ,ai4 so, dass stets i1 < i2 < i3 < i4.Konstruktion der affinen Gebietsabbildung FK so, dass ak 7→ aikfür k = 1, 2, 3, 4, d.h.

BK =[ai1 − ai4 ,ai2 − ai4 ,ai3 − ai4

], bK = ai4 .


Numerische Simulation mit finiten Elementen [-.2em...

Documents

Transcript of Numerische Simulation mit finiten Elementen [-.2em...