Numerische Strömungsmechanik

7/23/2019 Numerische Strömungsmechanik

http://slidepdf.com/reader/full/numerische-stroemungsmechanik 1/71

Technische Universität Ilmenau

Fakultät für Maschinenbau

FG Thermo- und Magnetofluiddynamik

Skript zur Vorlesung

Numerische Strömungsmechanik

Lehrbeauftragter: Dr. Thomas Boeck

Raum: M 401

E-Mail [email protected]

Stand 6. August 2015



INHALTSVERZEICHNIS I

Inhaltsverzeichnis

1 Einleitung 1

1.1 Grundlagen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1

1.2 Ursprünge, Geschichte der CFD . . . . . . . . . . . . . . . . . . . . . . . . 2

1.3 Verwendung von CFD . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2

1.4 Schritte einer CFD-Analyse . . . . . . . . . . . . . . . . . . . . . . . . . . 2

1.5 Verwendete Methoden . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3

1.5.1 Finite-Differenzen-Methode . . . . . . . . . . . . . . . . . . . . . . 3

1.5.2 Finite-Volumen-Methode . . . . . . . . . . . . . . . . . . . . . . . . 3

1.5.3 Finite-Element-Methode . . . . . . . . . . . . . . . . . . . . . . . . 3

1.5.4 Spektralverfahren . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4

2 Grundgleichungen der Strömungsmechanik 5

2.1 Eulersche- und Lagrangesche Beschreibung . . . . . . . . . . . . . . . . . . 5

2.2 Divergenz des Geschwindigkeitfelds, materielles Volumen . . . . . . . . . . 6

2.3 Massenerhaltung für Fluidelemente . . . . . . . . . . . . . . . . . . . . . . 7

2.4 Impulserhaltung für Fluidelemente . . . . . . . . . . . . . . . . . . . . . . 7

2.5 Energieerhaltungsgleichung . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

2.6 Randbedingungen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12

3 Partielle Differentialgleichungen, Klassifizierung 13

3.1 korrekt gestelltes Problem . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

3.2 Skalare PDE 1. Ordnung mit zwei unabhängigen Variablen . . . . . . . . . 13

3.3 Skalare PDE 2. Ordnung mit zwei unabhängigen Variablen . . . . . . . . . 16

3.4 Hyperbolischer Fall . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

3.5 Parabolischer Fall . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21

3.6 Elliptischer Fall . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22

3.7 Beispiel für schlecht gestelltes Problem . . . . . . . . . . . . . . . . . . . . 23

Numerische Strömungsmechanik Dr. Thomas Boeck



INHALTSVERZEICHNIS II

3.8 Klassifikation bei mehr als zwei unabhängigen Variablen . . . . . . . . . . 25

4 Grundlegendes zu Finite-Differenzen-Verfahren 26

4.1 Rechengitter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26

4.2 Differenzquotienten und Abbruchfehler . . . . . . . . . . . . . . . . . . . . 27

4.3 Approximation von PDEs . . . . . . . . . . . . . . . . . . . . . . . . . . . 31

4.4 Modifizierte Gleichung, numerische Diffusion/Dispersion . . . . . . . . . . 34

4.5 Stabilität von FDA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36

4.6 Von-Neumann-Analyse in einer Raumdimension . . . . . . . . . . . . . . . 37

4.7 Crank-Nicolson-FDA der Wärmeleitungsgleichung . . . . . . . . . . . . . . 39

4.8 Methoden für die eindimensionale Advektionsgleichung . . . . . . . . . . . 40

4.8.1 Beispiel – Leapfrog-Verfahren . . . . . . . . . . . . . . . . . . . . . 41

4.8.2 Beispiel – Lax-Wendroff-Verfahren . . . . . . . . . . . . . . . . . . 42

4.8.3 Beispiel – MacCormack-Verfahren . . . . . . . . . . . . . . . . . . . 43

4.8.4 Stabilität für Burgers-Gleichung . . . . . . . . . . . . . . . . . . . . 44

4.9 Verfahren für zweidimensionale Diffusionsprobleme . . . . . . . . . . . . . 45

4.9.1 Alternative mit geringerem Aufwand = „approximate factorization“ 46

4.9.2 Weitere Alternative = ADI (alternating directions implicit) . . . . . 47

4.10 Methoden für elliptische Probleme . . . . . . . . . . . . . . . . . . . . . . . 48

5 Finite-Differenzen-Verfahren/Strömungsprobleme 55

5.1 Stromfunktions-Wirbelstärke-Formulierung . . . . . . . . . . . . . . . . . . 555.2 Rolle des Drucks in inkompressiblen Strömungen . . . . . . . . . . . . . . . 57

5.3 Projektionsverfahren für zeitabhängige inkompressible Strömungen . . . . . 58

5.4 Iterative Lösung von stationären Strömungen . . . . . . . . . . . . . . . . 61

6 Finite-Volumen-Methode 63

6.1 Integrale Form der Grundgleichungen . . . . . . . . . . . . . . . . . . . . . 63

6.2 globale Erhaltungseigenschaften . . . . . . . . . . . . . . . . . . . . . . . . 64




INHALTSVERZEICHNIS III

6.3 1D-Beispiel für Finite-Volumen-Diskretisierung . . . . . . . . . . . . . . . . 66

6.4 2D-Beispiel für Finite-Volumen-Diskretisierung . . . . . . . . . . . . . . . . 67




1 EINLEITUNG 1

1 Einleitung

1.1 Grundlagen

• Strömungsvorgänge werden durch Feldgrößen beschrieben, die von Ort und Zeit

abhängen → Geschwindigkeit, Druck, Temperatur, . . .

• diese Größen werden durch die Erhaltungsgleichungen der Kontinuumsmechanik

bestimmt → Erhaltung von Masse, Impuls und Energie

• diese Gleichungen sind nichtlineare partielle Differentialgleichungen, Ableitung nach

Ort und Zeit• analytische Lösungen nur in stark idealisierten Fällen möglich, z.B.: Rohrströmung

• für die meisten realen Geometrien müssen die Gleichungen numerisch gelöst wer-

den → Computational Fluid Dynamics (CFD) = Entwicklung und Anwendung von

numerischen Methoden zur Lösung der strömungsmechanischen Gleichungen am

Computer

• Computerberechnung erfordert „Umformulierung“ der partiellen Differentialgleichun-

gen in ein System algebraischer Gleichungen = Diskretisierung

• die kontinuierlichen Felder werden durch eine endliche Anzahl von Zahlenwerten

(z.B. auf einem Gitter) ersetzt

• man muss sicherstellen, dass die Lösung der algebraischen Gleichung die tatsächliche

Lösung der partiellen DGL approximiert

• für Simulation turbulenter Strömungen gibt es verschiedene Arten von Simulationen

• direkte numerische Simulation (DNS)

– alle Wirbelstrukturen der Strömung werden auf dem Rechengitter abgebildet

– keine zusätzlichen Annahmen/Terme in Erhaltungsgleichung notwendig

– saubere aber extrem aufwendige Lösung

• RANS-Simulation (Reynolds-averaged Navier–Stokes)

– Lösung der zeitgemittelten Erhaltungsgleichungen für das mittlere Strömungs-

feld




1 EINLEITUNG 2

– Erhaltungsgleichungen enthalten Reynoldsschen Spannungstensor, der model-

liert wird

– für praktische Anwendungen üblich, Ergebnisse mit Vorsicht zu genießen

• Large-Eddy-Simulation (LES) = Grobstruktursimulation Eddy = Wirbel

– Erhaltungsgleichungen werden räumlich gefiltert

– kleinskalige Strukturen müssen durch Schließungsterm dargstellt werden

– Entwicklung von Schließungsmodellen gehört auch zu CFD

• Aufwand an Speicher + Rechenzeit DNS » LES » RANS

1.2 Ursprünge, Geschichte der CFD

• erste Überlegungen zur numerischen Wettervorhersage von Richardson in den 1920er

• rasche Entwicklung nach 1945 durch Entwicklung von Computern und Bedarf im

Militärbereich

• 1950er – 1960er Jahre 2D-Berechnungen

• 1970er Jahre 3D-Berechnungen

• seit 1980er Jahre kommerzielle Programme

1.3 Verwendung von CFD

• Produktentwicklung in der Industrie

• Forschungswerkzeug für Strömungsmechanik, Klimaforschung, Astrophysik, . . .

1.4 Schritte einer CFD-Analyse

• Problemstellung/-definition (Geometrie, Gleichungen, Randbedingungen, . . . )

• Diskretisierung, Rechengitter

• Simulation

• Analyse der Resultate, Reduktion der Informationsmenge, Visualisierung („colors

for directors“)




1 EINLEITUNG 3

1.5 Verwendete Methoden

→ numerische Verfahren sind auf Basis verschiedener Diskretisierungen möglich

1.5.1 Finite-Differenzen-Methode

• räumliches Punktgitter

• ersetzen der partiellen Ableitungen durch Differentialquotienten

∂u

∂x →

ui+1 − ui

xi+1 − xi

i = 1, 2, 3,...

1.5.2 Finite-Volumen-Methode

• Zerlegung des Strömungsgebiets in Zellen

• Verwendung integraler Formulierung der Erhaltungsgleichungen in jeder Zelle

• auf komplexe Geometrien anwendbar, globale Erhaltungseigenschaften

1.5.3 Finite-Element-Methode

• Gebietszerlegung in Zellen, Verwendung von Basisfunktion Φi, global definiert, aber

nur auf einem Element (und dessen Nachbarn) = 0

• Verwendung einer integralen (schwächeren) Formulierung der Erhaltungsgleichun-

gen

u =n

i=1

ui, Φi(x) i = 1, 2, 3, ...




1 EINLEITUNG 4

1.5.4 Spektralverfahren

• Reihenentwicklung in Basisfunktionen, die nicht auf Teilgebiete konzentriert sind,

z.B.: Fourier-Reihe

u =n

k=1

uk eikx

• beste Konvergenzeigenschaften (sehr genau)

• nur auf einfache Geometrien anwendbar




2 GRUNDGLEICHUNGEN DER STRÖMUNGSMECHANIK 5

2 Grundgleichungen der Strömungsmechanik

→ kontinuumsmechanische Beschreibung durch Feldgrößen→ Konzept des Fluidelements, lokales thermodynamisches Gleichgewicht (Zustandsglei-

chung gilt lokal)

2.1 Eulersche- und Lagrangesche Beschreibung

Lagrangesche Beschreibung = Bewegung aller Fluidelemente erfasst

x(t, x0, t0) = Trajektoren

u = dx

dt = Geschwindigkeit

Eulersche Beschreibung = man verfolgt nicht die Bewegung der einzelnen Fluidelemen-

te, sondern betrachtet nur die Eigenschaften, die an einem gegebenen Ort r und Zeit t

vorliegen.

u(r, t)

substantielle Ableitung liefert zeitliche Änderung für ein materielles Fluidelement, z.B.

der Dichte.

(r, t) −−−−→dt

(r + u dt, t + dt) = Änderung der Dichte über Zeitintervall dt

→ D

Dt() =

1

dt [(r + u dt, t + dt) − (r, t)] = Dichtedifferenz für dt → 0

mit (r + u dt, t + dt) = (r, t) + ∇ · u dt + ∂∂t

dt = Taylorreihe

→ D

Dt () = ∂

∂t + (u · ∇)





2.2 Divergenz des Geschwindigkeitfelds, materielles Volumen

Materielles Volumen besteht immer aus denselben Fluidelementen, die mit der Strömungtransportiert werden → nicht ortsfest. Notation für materielles Volumen = υ.

Änderung von υ während eines infinitesimalen Zeitintervalls dt.

ds = Oberflächenelement mit Flächennomale n

υ = ds n · u dt n · u dt = Projektion auf Senkrechte über ds

Gesamte Änderung von υ im Zeitraum dt.

dυ =

sn · u dt ds →

dυ

dt =

s

u · n ds −−−→Gauß

dυ

dt =

υ

∇ · n dυ

Man schreibt auch DυDt

, da materielles Volumen betrachtet wird (substantielle Ableitung).

Betrachtung von Grenzfall υ → 0

→ limυ→0

υ

∇ · u dυ

= υ ∇ · u

Man bezeichnet jetzt υ als δυ, da υ → 0 geht.

1

δυ

D

Dt(δυ) = ∇ · u

= relative Volumenänderung eines infinitesimalen materiellen Volumens pro Zeit.





2.3 Massenerhaltung für Fluidelemente

Fluidelement hat Volumen δυ und Masse δm = δυ

D

Dt(δm) = 0 = Massenerhaltung

→ D

Dt( δυ) = δυ

D

Dt() +

D

Dt(δυ)

δυ ∇· u

= 0

→ D

Dt + (∇ · u) = 0

→ ∂

∂t + (u · ∇) + (∇ · u) = 0

→ ∂∂t

+ ∇ · ( u) = 0

Für = konst. gilt ∇ · u = 0 (= Massenerhaltung für inkompressible Strömung)

2.4 Impulserhaltung für Fluidelemente

Ziel = Anwendung des Newtonschen Gesetz F = m a auf δυ.

m = δυ und a = D

Dt (u)

=

∂u

∂t

→ F = m a = δυ

D

Dt (u)

Welche Kräfte wirken auf das Fluidelement?

Volumenkräfte werden durch Kraftdichte beschrieben

→ Schwerkraft = F = m g = δ g

→ Lorentzkraft = F = j × B δυ

Oberflächenkräfte werden durch Spannungstensor beschrieben.

δ F = t δA

t= Spannungsvektor

hängt von Orientierung der Fläche ab





Man kann durch Betrachtung der Kraftbilanz auf einem Tetraeder zeigen, dass die Kraft

δ F auf eine Fläche δA, die nicht in einer der Koordinatenebenen liegt, durch δ F = n↔t δ F

gegeben ist.

→↔t = Doppelvektor (Tensor)

Resultierende Kraft von der Oberfläche des würfelförmigen Fluidelements in x-Richtung.

Mittelpunkt bei x, y, z

Betrag der oberen Deckfläche = τ zx(x,y,z + δz2 ) δx δy

Betrag der unteren Deckfläche = τ zx(x,y,z − δz2 ) δx δy

→

τ zx(x , y , z + δz2 ) − τ zx(x,y,z − δz

2 )

δx δy = ∂τ zx

∂z δx δy δz

δυ

Analog für die anderen Ebenen

→ ∂τ

xx∂x δx δy δz

δυ

→ ∂τ

yx∂y δx δy δz

δυ





Kraftkomponenten in x-Richtung ist somit

F x = ∂τ xx

∂x + ∂τ yx

∂y + ∂τ zx

∂z δυ

Analog für y- und z-Richtung

Einsetzen ins Newtonsche Gesetz F = m a, erhält man für die x-Richtung

δυ D

Dt(ux) =

∂τ xx

∂x +

∂τ yx

∂y +

∂τ zx

∂z

δυ + δυ gx + weiter Volumenkräfte

In Vektorschreibweise erhält man DDt

(u) = g f

+ ∇ ·↔τ = Impulsbilanz

Problem mit ↔τ

• Bestimmung von ↔τ

•

↔

τ hat einen geschwindigkeitsunabhängigen Anteil

• Druck wirkt auch in ruhenden Flüssigkeiten

• Druck wirkt allseitig, hängt nicht von Orientierung ab (isotrop)

• ↔

τ (0)

= − p I I = Einheitsmatrix → τ ij(0) = − p δ ij

• geschwindigkeitsabhängiger Anteil von ↔τ hängt von Geschwindigkeitsgradienten ∇u

ab

Einfachster Ansatz, ↔τ

(I )hängt linear von ∇u ab. Der Spannungstensor erfüllt τ ij = τ ji

(Symmetrie) aufgrund der Drehimpulserhaltung, außerdem resultieren Spannungen nur

aus Scherung.

→ τ ij(I ) hängt nur vonS ij =

1

2

∂ui

∂x j+

∂u j

∂xi

ab (Scherratentensor)

→ τ ij(I ) =

u,l

Λijkl S kl





Aus Gründen der Isotropie verbleiben vom Λijkl nur ein Beitrag.

→ τ ij

(I )

= 2 µ S ij µ = dynamische Viskosität

Zusammengesetzt ergibt sich ↔τ

τ ij = τ ij(0) + τ ij

(I ) = − p δ ij + 2 µ S ij

= Materialgesetz für inkompressible Newtonsche Fluide

Einsetzen in Impulsbilanz ergibt

D

Dt(u) = f + ∇τ →

D

Dt(ui) = f i +

j

∂

∂x j

τ ij

mit

j

∂

∂x j

τ ij

=

j

∂

∂x j(− p δ ij)

+

j

∂

∂x j(2 µ S ij)

Getrennt auflösen ergibt

j

∂

∂x j

(− p δ ij)

=

j

−

δ ij

∂p

∂x j

= −

∂p

∂xi

= −(∇ p)i mit δ =

i = j → 1

i = j → 0

j

∂

∂x j

(2 µ S ij)

=

j

µ

∂

∂x j

(2 S ij)

=

j

µ

∂

∂x j

∂ui

∂x j

+ ∂u j

∂xi

= j

µ ∂ 2ui

∂x j 2 +

j

µ ∂ 2u j

∂x j ∂ xi

=µ ∂ ∂xi

∂uj∂xj

→

j

∂

∂x j

(2 µ S ij)

= µ ∇2 ui + µ

∂

∂xi

j

∂u j

∂x j

=0 ∇· u=0

mit ∇2 = ∂ 2

∂x2 +

∂ 2

∂y 2 +

∂ 2

∂z 2

→ es folgt die Navier-Stokes-Gleichung DDt

(u) = −∇ p + µ ∇2u + f





2.5 Energieerhaltungsgleichung

Für die Herleitung wendet man den ersten Hauptsatz der Thermodynamik auf ein Flui-delement an und bilanziert zu- bzw. abgeführte Wärme und Arbeit durch Oberflächen-

und Vol-menspannungen. Primär erhält man eine Gleichung für Gesamtenergie = innere

Energie + kinetische Energie (genauer für Energiedichten). Man zieht deshalb die reine

mechanische Energiegleichung für die kinetische Energiedichte ab.

D

Dt(e) = q + ∇ · (λ ∇T ) +

↔τ ∇ · u

=ijτ ij

∂ui

∂xj

ij

τ ij

∂ui

∂x j

= viskose Dissipation

q = volumetrische Heizung

λ ∇T = Wärmestromdichte

e = Dichte der inneren Energie

→ für inkompressible Fluide gilt e = cv T , wobei die viskose Dissipation vernachlässigbar

klein ist.

cv D

Dt(T ) = q + ∇ · (λ ∇T ) λ = Wärmeleitfähigkeit, T = Temperatur

→ D

Dt(T ) =

q

cv +

λ

cv ∇2T

λ

cv = κ = Temperaturleitfähigkeit

→ D

Dt(T ) =

∂T

∂t + (u · ∇) T = κ ∇2T ohne q





2.6 Randbedingungen

Betrachten nur inkompressible Strömungen mit ∇ · u = 0. Dadurch ist der Druck keineunabhängige dynamische Größe und es werden nur Randbedingungen für u benötigt.

• an festen Wänden gilt u = uWand (Haftbedingungen (no-Slip))

• an Grenzflächen zweier Fluide gilt u1 = u2, sowie lokales Kräftegleichgewicht, eine

Normalspannung und 2 Tangentialkräfte (Spannungen)

• kinematische Randbedingungen → Grenzfläche wird durch Strömung transportiert

• für Temperatur an Wänden gelten Randbedingungen 1., 2. oder 3. Art

– erster Art, Temperatur an Wand → T = T Wand

– zweiter Art, Wärmestromdichte an Wand → q Wand = λ n ∇T

– dritter Art, Beschreibung des konvektiven Wärmeübergangs von der Oberfläche

eines Körpers an ein vorbeiströmendes Fluid → q c = h (T W − T ∞)

• für Temperatur an Grenzflächen gilt T 1 = T 2 und Wärmestromdichte senkrecht zu

Grenzfläche stetig → λ1∂T

∂n 1 = λ2∂T

∂n 2




3 PARTIELLE DIFFERENTIALGLEICHUNGEN, KLASSIFIZIERUNG 13

3 Partielle Differentialgleichungen, Klassifizierung

3.1 korrekt gestelltes Problem

Gleichungen der Strömungsmechanik (Massen-, Impuls- und Energieerhaltung) sind parti-

elle Differentialgleichungen (PDE = partial differential equation). Man braucht zur Lösung

der PDE noch zusätzliche Bedingungen (Anfangs-/Randbedingungen). Dies ist analog zu

ge-wöhnlichen Differentialgleichungen.

z.B.: d2f

dx2 + a

df

dx + b = 0

Anfangswertproblem: f, df dx

an der Stelle x0 gegeben

Randwertproblem: Intervall (x1, x2) → f (x1) = A, f (x2) = B

→ f (x1) = A, df (x2)dx

= B

Frage: Wie sieht das für PDE aus?

Man möchte ein korrekt gestelltes Problem erhalten. Korrekt gestellt heißt, das eine Lö-

sung existiert, diese eindeutig ist und stetig von den Rand-/Anfangswerten abhängt und

ggf. noch von weiteren Parametern. Dies führt zur Betrachtung sogenannter charakteris-tischen Kurven/Flächen. Erläuterung anhand einfacher Beispiele.

3.2 Skalare PDE 1. Ordnung mit zwei unabhängigen Variablen

Beispiel für eine solche PDE → Betrachtung inkompressible 2D-Strömung in einem end-

lichen geschlossenen Gebiet.

u =

u(x, y)

v(x, y)

→ wir wollen Stromlinien als Isolinien einer Funktion ψ(x, y) darstellen





→ Eigenschaften einer Stromlinie

• dr muss parallel zu u seind dr = Länge der Stromlinie in x- und y-Richtung

dx = c u

dy = c v

→ dx

u =

dy

v → u dy = v dx

dψ = ψ(x + dx, y + dy) − ψ(x, y) = ∂ψ

∂x dx +

∂ψ

∂y dy

(u v)

→ dψ u v = u ∂ψ

∂x v dx + v

∂ψ

∂y u dy

• auf einer Stromline gilt u dy = v dx

dψ (u v) =

u

∂ψ

∂x + v

∂ψ

∂y

v dx

→ ψ = konst. ist eine Stromlinie, wenn u ∂ψ

∂x + v

∂ψ

∂y = 0 gilt → PDE

• damit ψ in dem Gebiet eindeutig bestimmt ist, muss man ψ auf einer Kurve g(x, y) =

0 vorgeben, die alle Stromlinien genau einmal kreuzt

• sie darf ein- und dieselbe Stromlinie nicht mehrfach schneiden oder tangieren

• sei nun x0(s), y0(s) einer Kurve mit ψ0(s) gegeben → diese Kurve ist eine charak-

teristische Kurve, wenn ψ nicht bestimmt ist

• anschaulich ist klar, dass die Kurve dann eine Stromlinie tangiert

• wir wollen zeigen, dass Ableitungen von ψ nicht berechnet werden können, wenn

x0(s), y0(s) lokal eine Stromlinie bilden

∂ψ0(s)

∂s =

∂ψ

∂x

∂x0

∂s +

∂ψ

∂y

∂y0

∂s gilt entlang der Kurve x0(s), ys(0)

→ u ∂ψ

∂x + v

∂ψ

∂y = 0





→ 2 Gleichungen für 2 Unbekannte ∂ψ∂x

, ∂ψ∂y

dx0ds

dy0ds

u v

∂ψ∂x∂ψ∂y

=

∂ψ0(s)∂s

0

→ keine Lösung, wenn Determinante

dx0ds

dy0ds

u v

= 0

→

dx0ds

dy0ds

||

u

v

parallel zueinander → auf Stromlinie

• Stromlinien sind also die charakteristischen Kurven

• solche Kurven beschreiben „Transport von Informationen“ durch die PDE

• Werte von ψ dürfen nicht auf charakteristischen Kurven gegeben werden

→ Bemerkung zur sogenannten Stromfunktion

Man kann für inkompressible 2D-Strömungen eine Stromfunktion ψ definieren. Die Isoli-

nien dieser Stromfunktion sind Stromlinien.

∂u

∂x +

∂v

∂y = 0 → u =

∂ψ

∂y und v = −

∂ψ

∂x

→ u ∂ψ

∂x + v

∂ψ

∂y = 0 → u v − v u = 0

Zur Berechnung der Stromfunktion ψ kann man z.B. so vorgehen

∂ 2ψ

∂x2 +

∂ 2ψ

∂y 2 =

∂u

∂y −

∂v

∂x





Überlegungen lassen sich auf quasi lineare PDE verallgemeinern.

a(x,y, Φ)

∂ Φ

∂x + b(x,y, Φ)

∂ Φ

∂y = c(x,y, Φ)

→ Φ(x, y) sei eine Lösungsfläche über x,y-Ebene

F (x,y, Φ) = Φ − Φ(x, y)

Die PDE kann auch so interpretiert werden

a

b

c

· ∇F = 0

Die Lösungsfläche wird vom Vektor [a b c]T tangiert. Wenn man eine Kurve x(s), y(s),Φ(s) hat, die tangential zu [a b c]T ist, dann liegt diese in der Lösungsfläche.

→ ∂ x

a =

∂ yb

= ∂ Φ

c = charakteristische Kurve

→ Startwerte dürfen nicht uf diesen Kurven gegeben werden.

3.3 Skalare PDE 2. Ordnung mit zwei unabhängigen Variablen

→ Charakteristische Kurve

a ∂ 2Φ

∂x2 + b

∂ 2Φ

∂x ∂y + c

∂ 2Φ

∂y 2 + d

∂ Φ

∂x + e

∂ Φ

∂y = f

PDE enthält 2. Ableitung, für Lösungen benötigt man Daten auf einer Kurve, die eindeu-

tige Berechnungen der zweiten Ableitung gestattet.

→ x0(s), y0(s), Φ0(s), ∂ Φ0(s)

∂x Φx(s)

, ∂ Φ0(s)

∂y Φy(s)

dΦx

ds =

∂ 2Φ

∂x2

dx0

ds +

∂ 2Φ

∂x ∂y

dy0

ds sowie

dΦy

ds =

∂ 2Φ

∂y 2

dy0

ds +

∂ 2Φ

∂x ∂y

dx0

ds





Zusammen mit der PDE sind das 3 Gleichungen für die zweite Ableitung von Φ

→ ∂ 2Φ

∂x2 ,

∂ 2Φ

∂y 2 ,

∂ 2Φ

∂x ∂y

→

a b c

dx0ds

dy0ds

0

0 dx0ds

dy0ds

∂ 2Φ∂x2

∂ 2Φ∂x ∂y∂ 2Φ∂y2

=

·

·

·

Das System ist nicht lösbar, wenn die Determinante = 0 ist

a

dy0ds

0dx0ds

dy0ds

− dx0

ds

b cdx0ds

dy0ds

= 0 f ür erste Spalte

→ a

dy0

ds

2

− b dx0

ds

dy0

ds + c

dx0

ds

2

= 0

:

dx0

ds

2

→ a

dy0

dx0

2

− b dy0

dx0+ c = 0 →

dy0

dx0=

b

2 a ±

b2

4 a2 −

c

a = λ1/2

Wenn gilt

• λ1/2 sind reell und verschieden → b2 − 4 a c > 0 (hyperbolisch)

• λ1 = λ2 → b2 = 4 a c (parabolisch)

• λ1/2 sind konjugiert komplex → b2 − 4 a c < 0 (elliptisch)

• Bezeichnungen kommen von den Kegelschnitten





Die Differentialgleichung für dy0dx0

stellt wegen λ1 = λ2 zwei verschiedene Lösungen dar.

• y(x, c1) aus λ1 und y(x, c2) aus λ2. c1 und c2 sind Parameter und können als neueKoordinaten verwendet werden.

• c1 = ξ und c2 = η

• neue Koordinaten sind ξ (x, y) und η(x, y). Dies ist nur für den hyperbolischen Fall

möglich, für elliptischen und parabolischen Fall muss ξ und η anders festgelegt.

→ Veranschaulichung

λ1, λ2 = konstant

λ1 = λ2

Man transformiert durch Verwendung von ξ (x, y) und η(x, y) die PDE in die sogenanntekanonische Form.

hyperbolisch: 1. Art ∂ 2Φ

∂ξ ∂η = h1

∂ Φ

∂ξ ,

∂ Φ

∂η , Φ, ξ , η

2. Art ∂ 2Φ

∂ξ 2 −

∂ 2Φ

∂η2 = h2

∂ Φ

∂ξ ,

∂ Φ

∂η , Φ, ξ , η

parabolisch:

∂ 2Φ

∂ξ 2 = h3 ∂ Φ

∂ξ ,

∂ Φ

∂η , Φ, ξ , η

elliptisch: ∂ 2Φ

∂ξ 2 +

∂ 2Φ

∂η2 = h4

∂ Φ

∂ξ ,

∂ Φ

∂η , Φ, ξ , η

Betrachten im Folgenden lineare PDEs in kanonischer Form mit konstanten Koeffizienten

zur Veranschaulichung.





3.4 Hyperbolischer Fall

Es gilt, dass a und c konstant sind → λ1, λ2 sind konstantdy

dx = λ1 → y − λ1 · x = ξ Charakteristiken zu λ1

dy

dx = λ2 → y − λ2 · x = η Charakteristiken zu λ2

→ ∂ 2Φ

∂ξ ∂η = ...

Wenn man ξ = 12

(ξ + η) und η = 12

(ξ − η) benutzt, erhält man die 2. Art ∂ 2Φ

∂ξ2 − ∂ 2Φ

∂η2 = ...

Beispiel Wellengleichung:

∂ 2u

∂t2 = −c0

2 ∂ 2u

∂x2 = 0 t → x, x → y zum Vergleich

a ∂ 2u

∂ x2 + b ∂ 2u

∂ x ∂ y + c ∂ 2u

∂ y2 = 0 = allgemeine Form der PDE

→ a = 1, b = 0, c = −c02 durch Vergleich mit allgemeiner Form einer PDE 2. Ordnung.

→ b2 − 4 a c > 0 (hyperbolisch) → dydx = b

2 a ±

b2

4 a2 − c

a = ± c0 = λ1/2

→ dxdt

= ± c0 = λ1/2, daraus ergeben sich y + c0 x = ξ und y − c0 x = η

Transformation

→ ∂u(ξ, η)

∂t =

∂u

∂ξ

∂ξ

∂t +

∂u

∂η

∂η

∂t

→ ∂

∂t

∂u

∂t

=

∂

∂t

∂u

∂ξ

∂ξ

∂t

+

∂

∂t

∂u

∂η

∂η

∂t

→ ∂

∂t ∂u

∂t = ∂ξ

∂t c0

∂ 2u

∂ξ 2

∂ξ

∂t c0

+ ∂ 2u

∂η ∂ξ

∂η

∂t −c0

+

∂η

∂t −c0

∂ 2u

∂η2

∂η

∂t −c0

+ ∂ 2u

∂η ∂ξ

∂ξ

∂t c0





Analog für x, es verleibt dann:

∂ 2u

∂η ∂ξ = 0 ξ−−→ ∂u

∂η = G(η) η−−−→ u = F 2(η) + F 1(ξ )

Allgemeine Lösung → u = F 1(ξ ) + F 2(η)

Anfangswertproblem bei t = 0 sei u(x,t=0) = f (x); ∂u(x,t=0)

∂t = g(x)

→ Ist verträglich mit der PDE, weil t = 0 keine charakteristische Kurve ist.

u = F 1(ξ ) + F 2(η) = f (x)

∂u∂t

= dF 1

dξ dξ dt

+ dF 2

dηdηdt

= c0dF 1

dξ − dF 2

dη

ξ=η=x

= g(x)

→ F 1(x) − F 2(x) = 1

c0

x

0g(τ ) dτ A

Auflösen nach F 1, F 2

F 1(x) = 1

2 f (x) +

1

c0 x

0g(τ ) dτ A

F 2(x) =

1

2

f (x) −

1

c0

x

0g(τ ) dτ A

Daraus ergibt sich die D’Alembertsche Lösung

u(x, t) = F 1(x + c0 t) + F 2(x − c0 t)

= 1

2 (f (x + c0 t) + f (x − c0 t)) +

1

2 c0 x+co t

x−c0 tg(τ ) dτ

u hängt nur von f und g im Intervall [x − c0 t; x + c0 t]

ξ = x + c0 t

η = x − c0 t

Diese Eigenschaft der charakteristischen Kurve (Begrenzung von Abhängigkeits- und Ein-

flussbereichs) trifft auch im allgemeinen Fall zu.





3.5 Parabolischer Fall

dy

dx =

b

2 a = λ1 = λ2

Charakteristiken können nicht zur Definition neuer Koordinaten verwendet werden. Man

braucht noch andere Kurvenscharen z.B.:

dy

dx = λ = λ2 → y − λ1 x = η

→ y − λ x = ξ

Ausgangsgleichung ∂ 2Φ∂ξ2

= g(...)

Typischerweise ist höchste Ableitung in η auf rechter Seite als linearer Term vorhanden.

∂ 2Φ

∂ξ 2 = α

∂ Φ

∂η + ... α = Konstante

Beispiel – Diffusionsgleichung

∂ Φ

∂t = κ

∂ 2Φ

∂x2 → a = κ, b = 0, y → t

dt

dx = 0 → t = konstant = charakteristische Kurve

→ können nicht Φ und ∂ Φ∂t

bei t = 0 vorgegeben, stattdessen Φt=0 + Randbedingungen

bei t > 0.





3.6 Elliptischer Fall

dy

dx = c1 ± i c2 → komplexe Charakteristiken machen keinen Sinn

formal: ξ = y − λ1 x = y − c1 x + i c2 x

η = y − λ2 x = y − c1 x − i c2 x

stattdessen: ξ = Re(ξ ) und η = I m(ξ )

ξ = y − c1 x η = y − c2 x

→ ∂

2

u∂ξ 2

+ ∂

2

u∂η 2 = ln(...) generische Form

Beispiel - Laplace-Gleichung

∂ 2u

∂x2 +

∂ 2u

∂y 2 = 0

Anfangswertproblem ist nicht wohlgestellt (u, ∂u∂y bei y = 0)

→ Vorgabe von Randwerten auf geschlossenem Rand

Vorgabe von u auf Rand → Dirichlet-Problem

Vorgabe von ∂u

∂n auf Rand → Neumann-Problem

→ Änderung der Randwerte wirkt Global





3.7 Beispiel für schlecht gestelltes Problem

∂ 2Φ

∂x2 +

∂ 2Φ

∂y 2 = 0

Versuchen Φ, ∂ Φ

∂y auf x-Achse vorzugeben

∂ Φ

∂y = 0 bei y = 0, Φ(x, 0) =

∞k=0

(ak cos(k x))

Versuchen Lösung mit Trennung der Variablen

Φ = f (x) g(y) ∂ 2Φ

∂x2 = g

d2f

dx2

∂ 2Φ

∂y 2 = f

d2g

dy2

Ansatz

∂ 2Φ

∂x2 +

∂ 2Φ

∂y 2 = 0 → g

d2f

dx2 + f

d2g

dy2 = 0

: (f g)

→ 1

f

d2f

dx2 +

1

g

d2g

dy2 = 0

→ d2f

dx2 = −k2 f → f = cos(k x)

→ d2

gdy2

= k 2 g → g = Ak cosh(k y) + Bk sinh(k y)

Φ =∞

k=0

cos(k x) Ak cosh(k y) + Bk sinh(k y)

→ ∂ Φ

∂y =

∞k=0

cos(k x) Ak k sinh(k y) + Bk k cosh(k y)

für y = 0 gilt

∂ Φ

∂y = 0 → 0 =

∞k=0

cos(k x)

Ak k sinh(k 0) =0

+ Bk k cosh(k 0) =1

→ 0 =

∞k=0

cos(k x) Bk k → Bk = 0





eingesetzt ergibt

→ Φ =

∞

k=0 Ak cos(k x) cosh(k y)

für y = 0 gilt

Φ(x, 0) =∞

k=0

ak cos(k x) → Ak = ak

→ Φ(x, y) =∞

k=0

ak cos(k x) cosh(k y)

Abschätzung Konvergenz, y » 1

→ cosh(k y) ∼ 1

2 ek y

→ Φ(x, y) ∼ 1

2

∞k=0

ak ek y cos(k x)

Konvergenz der Reihe erfordert, dass |ak| ek y cos(k x) → 0, mit k → ∞ Diese Forderungist für genügend Große Werte von y im Allgemeinen nicht erfüllt. → Problem schlecht

gestellt, keine kontinuierliche Abhängigkeit von den Anfangswerten.





3.8 Klassifikation bei mehr als zwei unabhängigen Variablen

Eine Klassifikation ist schwierig, Charakteristiken lassen sich formal als Hyperflächen ein-führen aber Auswertung/Interpretation kompliziert. Man kann sich behelfen, in dem man

z.B. eine der unabhängigen Variablen als homogen annimmt.

Beispiel – 2D-Wärmeleitung

∂T

∂t = κ

∂ 2T

∂x2 +

∂ 2T

∂y 2

→ elliptischer Fall ∂

∂t = 0

→ parabolischer Fall ∂ ∂x

= 0, bzw. ∂ ∂y

= 0

Weitere Komplikationen, wenn Vektorfelder gesucht sind. Generell ist Unterscheidung

zwischen hyperbolischen/parabolischen und elliptischen Problemen wesentlich. Hyperbo-

lische/parabolische Probleme sind durch Zeitintegration („time marching“) zu lösen, el-

liptische Probleme erfordern globales Lösungsverfahren. Es gibt auch Situationen, wo eine

Ortskoordinate als zeitartige Variable auftritt, das ist z.B. bei den Gleichungen einer sta-

tionären Grenzschicht der Fall.

u ∂

∂x + v

∂

∂y = −

dp

dx + +ν

∂ 2u

∂y 2

p(x) sei gegeben

Für inkompressible reibungsbehaftete Strömungen liegt ein parabolisches Problem vor.

→ Druck spielt eine Sonderrolle

Diskussion der Grundlagen mittels Modellgleichungen für die verschiedenen Typen.

lineare Advektionsgleichung ∂ Φ

∂t + c

∂ Φ

∂x = 0 hyperbolischer Fall

Wärmeleitungsgleichung ∂ Φ

∂t = κ

∂ 2Φ

∂x2 parabolischer Fall

Laplace-Gleichung ∂ 2Φ

∂x2

∂ 2Φ

∂y 2 = 0 elliptischer Fall




4 GRUNDLEGENDES ZU FINITE-DIFFERENZEN-VERFAHREN 26

4 Grundlegendes zu Finite-Differenzen-Verfahren

FDA – Finite-Differenzen-ApproximationWie erhält man FDA aus einer PDE?

Welcher Fehler tritt dabei auf? (Genauigkeit)

4.1 Rechengitter

Betrachten ein Problem, dass durch Zeitschnittverfahren zu lösen ist. Dazu erfolgt eine

Unterteilung von x und t in „Schichten“, die Lösung ist dann am Punkt (xi, tn) zu be-

stimmen.

tn = n ∆t konstanter Zeitabschnitt ∆t

äquidistantes Gitter in t

tn − tn−1 Variabel mit n

Variabler Zeitabschnitt nicht äquidistant

xi = i ∆x äquidistant in x

xi − xi−1 Variabel mit i, nicht gleichförmig

Man erhält auf diese Weise ein strukturiertes Gitter, bei mehreren Raumdimensionen sind

unstrukturierte/irreguläre räumliche Gitter möglich. Für FDV nicht üblich.

Diskussionen erfolgen anhand der klassischen FDA mit strukturierten Gittern.

Φni = Φ(xi, tn) = Lösungswerte





4.2 Differenzquotienten und Abbruchfehler

Umwandlung der PDE zu FDA in dem man Ableitungen durch Differenzquotienten dar-stellt.

∂u

∂x = lim

∆x→0

u(x + ∆x) − u(x)

∆x

→

∂u

∂x

xi

∼= u(xi+1) − u(xi)

xi+1 − xi

Dabei wird ∆x = xi+1 − xi als genügend klein vorausgesetzt.

Frage: Welcher Fehler tritt dabei auf?

→ Taylorentwicklung – Vorwärtsdifferenz

u(x + ∆x) = u(x) + ∂u

∂x

x

∆x + 1

2

∂ 2u

∂x2

x

(∆x)2 + ... + 1

n!

∂ nu

∂xn

ξ

(∆x)n

mit x < ξ < (x + ∆x)

→ u(x + ∆x) − u(x)

∆x =

∂u

∂x

x

+ 1

2

∂ 2u

∂x2

ξ

∆x

→ ∂u

∂x

xi

= u(xi+1) − u(xi)

xi+1 − xi−

1

2

∂ 2u

∂x2

ξ

∆x Abbruchfehler T.E.

T.E. steht für Truncation Error

→ Taylorentwicklung – Rückwärtsdifferenz

u(xi − ∆x) = u(xi) − ∂u

∂x

xi

∆x + 1

2

∂ 2u

∂x2

xi

(∆x)2 − ... + (−1)n 1

n!

∂ nu

∂xn

ξ

(∆x)n

→ ∂u

∂x

xi

= u(xi) − u(xi−1)

∆x +

1

2

∂ 2u

∂x2

ξ

∆x Abbruchfehler T.E.

Der Abbruchfehler ist der Betrag, der bei der Approximation der Ableitung durch finite

Differenzquotienten vernachlässigt wird. Größe des T.E. ist nicht bekannt, nur Verhalten

mit ∆x. Für Vorwärts- und Rückwärtsdifferenz ist der Abbruchfehler 0(∆x).

Das Symbol 0(h) bedeutet f (h) = 0(hn) ⇔ |f |hn

< k für h → 0, k = Koeffizient





Frage: Wie kann man den Abbruchfehler verkleinern?

Idee: Mehr Stützstellen für Approximation verwenden

→ Stützstellen xi−1 xi xi+1 mit ∆x = konst = xi − xi−1 unabhängig von i

Taylorentwicklung ergibt

für xi+1 → u(xi+1) = u(xi) + ∂u

∂x

xi

∆x + 1

2

∂ 2u

∂x2

xi

∆x2 + 0(∆x3)

für xi−1 → u(xi−1) = u(xi) − ∂u

∂x

xi

∆x + 1

2

∂ 2u

∂x2

xi

∆x2 + 0(∆x3)

für xi → u(xi) = u(xi)

es gilt → α u(xi+1) + β u(xi) + γ u(xi−1) = ∂u

∂x

xi

Einsetzen und Betrachtung der einzelnen Komponenten

0. Ableitung u(xi) α u(xi) + β u(xi) + γ u(xi) = 0→ α + β + γ = 0

1. Ableitung ∂u

∂x

xi

α

∂u

∂x

xi

∆x

+ β 0 + γ

− ∂u

∂x

xi

∆x

= ∂u

∂x

xi

→ α ∆x − γ ∆x = 1

2. Ableitung ∂

2

u∂x2 xi

α 12 ∂

2

u∂x2 xi

(∆x)2 + β 0 + γ 12 ∂

2

u∂x2 xi

(∆x)2 = 0

→ α + γ = 0 α = −γ

Zusammengefasst ergibt dies α = 12 ∆x

, β = 0 und γ = − 12 ∆x

.

Daraus ergibt sich 12 ∆x

(u(xi+1) − u(xi−1)) = ∂u∂x

xi+ 0(∆x2)

→ zentrale Differenz = Mittelwert von Vorwärts- und Rückwärtsdifferenz





Frage: Was macht man für höhere Ableitungen?

Beispiel

∂ 2u

∂x2 → v =

∂u

∂x →

∂ 2u

∂x2 =

∂v

∂x

→ ∂v

∂x

xi

= vi+1 − vi

∆x + 0(∆x)

→ vi+1 = ∂u

∂x

xi+1

= ui+2 − ui+1

∆x + 0(∆x)

→ vi = ∂u

∂x

xi

= ui+1 − ui

∆x + 0(∆x)

einsetzen ergibt

∂v

∂x

xi

= 1

∆x

ui+2 − ui+1

∆x + 0(∆x) −

ui+1 − ui

∆x + 0(∆x)

+ 0(∆x)

→ ∂v

∂x

xi

= 1

∆x2 (ui+2 − 2 ui+1 + ui) + 0(∆x0)

mit 0(∆x)∆x

= 0(∆x0), ∆x0 = 1, Abbruchfehler somit nicht verkleinerbar, zu pessimistisch

→ ui+2 − 2 ui+1 + ui

∆x2 kann jedoch auf Abbruchfehler untersucht werden

Ergebnis

∂ 2u

∂x2

xi

= ui+2 − 2 ui+1 + ui

∆x2 + 0(∆x)

→ für höhere Ableitungen sind mehr Zwischenschritte notwendig

Beispiel ∂u

∂x = ... + 0(∆x4) → 4. Ableitung

→ 5 Stützstellen xi+2 xi+1 xi xi−1 xi−2

→ Taylorentwicklung and diesen Stellen

u1+2 = ui + ...

→ α u(xi+2) + β u(xi+1) + γ u(xi) + δ u(xi−1) + u(xi−2) = ∂u

∂x

xi

+ 0(∆x4)

→ ∂u

∂x xi

= 1

12 ∆x (−u(xi+2) + 8 u(xi+1) − 8 u(xi−1) + u(xi−2))





Vorgehen für höhere Ableitungen im Prinzip analog, man benötigt aber mehr Stützstellen

z.B. 3 Stellen für 2. Ableitung usw.

→ zentrale Differenz 2. Ableitung ∂ 2u

∂x2

xi

= ui+1 − 2 ui + ui−1

∆x2 + 0(∆x2)

Frage: Wie verhält sich das bei gemischter Ableitung nach x und y?

∂ 2u

∂x ∂y

xi,yj

= ∂

∂x

∂u

∂y

=

1

∆x

∂u

∂y

xi+1

− ∂u

∂y

xi

Vorwärtsdifferenz

+ 0(∆x)

→ ∂u

∂y

yj

= u j+1 − u j

∆y (Vorwärtsdifferenz)

Vereinfachung u(xi, y j ) = ui,j

∂ 2u

∂x ∂y

xi,yj

= 1

∆x

ui+1,j+1 − ui+1,j

∆y −

ui,j+1 − ui,j

∆y

+ 0(∆x, ∆y)

→ Man kann auch Stützstellen verwenden, die nur auf einer Seite von xi liegen, z.B. für

Ableitungen am linken Gebietsrand für höhere Auflösung.

→ ∂u

∂x

xi

= 1

2 ∆x (−3 ui + 4 ui+1 − ui+2) + 0(∆x2)

→ ∂u

∂x

xi

= 1

2 ∆x (3 ui − 4 ui−1 + ui−2) + 0(∆x2)





4.3 Approximation von PDEs

Wärmeleitungsgleichung als Beispiel∂T

∂t = κ

∂ 2T

∂x2 + f (x, t)

mit äquidistanten Gitter → xi = i ∆x tn = n ∆t

T (xi, tn) → T ni

partielle Ableitung durch Differenzenquotienten ersetzen

∂T

∂t

xi,tn

= T n+1

i − T ni∆t

(Vorwärtsdifferenz in Zeit)

∂ 2T

∂t2

xi,tn

= T ni+1 − 2 T ni + T ni−1

∆x2 (zentrale Differenz 2. Ordnung im Ort)

f (xi, tn) = f ni

einsetzen in PDE ergibt

T n+1i − T ni

∆t = κ

T ni+1 − 2 T ni + T ni−1

∆x2 + f ni

= finite Differenzen Approximation (FDA)

wichtig: Ableitungen müssen immer an gleicher Stelle approximiert werden.

Abbruchfehler des Verfahrens kann wie gehabt mit Taylorentwicklung bestimmt werden,man muss aber sowohl in x als auch in t entwickeln.

T (x ± ∆x, t) = T (x, t) ± ∂T (x, t)

∂x ∆x +

1

2

∂ 2T (x, t)

∂x2 ∆x2 ±

1

6

∂ 3T (x, t)

∂x3 ∆x3 +

1

24 ...

T (x, t + ∆t) = T (x, t) +

∂T (x, t)

∂t ∆t +

1

2

∂ 2T (x, t)

∂t2 ∆t2

+

1

6

∂ 3T (x, t)

∂t3 ∆t3

+

1

24 ...





T.E. = T n+1

i − T ni∆t

− κ T ni+1 − 2 T ni + T ni−1

∆x2 − f ni

→ T.E. = ∂T (x, t)

∂t +

1

2

∂ 2T (x, t)

∂t2 ∆t − κ

∂ 2T (x, t)

∂x2 +

1

12

∂ 4T (x, t)

∂x4 ∆x2

− f ni

→ T.E. = ∂T (x, t)

∂t − κ

∂ 2T (x, t)

∂x2 − f ni

=0 exakte Lösung

+ 1

2

∂ 2T (x, t)

∂t2 ∆t − κ

1

12

∂ 4T (x, t)

∂x4 ∆x2

→ T.E. = 0(∆t) + 0(∆x2)

FDA ist konsistent, wenn Abbruchfehler T.E. → 0 geht, falls ∆x → 0, ∆t → 0

• Konsistenz ist notwendig, damit Lösung der PDE korrekt approximiert wird

• Konsistenz ist dafür jedoch im Allgemeinen nicht ausreichend

Die in der FDA verknüpften Nachbarpunkte von xi, tn bilden den sogenannten Rechens-

tern (Stencil, finite-difference-molecule).

Die FDA (Vorwärtsdifferenz für Zeit und zentrale Differenz für Ort) heißt auch FTCS

(forward in time, centered in space). FTCS ist ein sogenanntes explizites Verfahren, d.h.

man erhält die Lösung auf Zeitebene (n+1) von Zeitebene (n) ohne ein Gleichungssystem

zu lösen.





Gegenteil von explizit ist implizit, Gleichungssystem muss gelöst werden.

Beispiel – Rückwärtsdifferenz in Zeit t, Stelle xi, tn+1

∂T

∂t = κ

∂ 2T

∂x2 + f (x, t)

→ ∂T

∂t xi,tn

= T n+1

i − T ni

∆t

→ ∂ 2T

∂t2 xi,tn

= T n+1

i+1 − 2 T n+1i + T n+1

i−1

∆x2

→ T n+1i

1 +

2 κ ∆t

∆x2

α

− T n+1i+1

κ ∆t

∆x2 β

− T n+1i−1

κ ∆t

∆x2 β

= T ni + f n+1i ∆t

T sei gesucht zwischen x0 und xN , Randwerte von T bei x0, xN seien bekannt T n0 , T nN

1 0 0 · · ·

−β α −β 0 · ·

0 −β α −β 0 ·

0 0 · · · ·

· · · · · ·

· · · 0 0 1

T n+10

T n+11

T n+12

·

·

T n+1N

=

T (x0, tn+1)

T n1 + f n+11 ∆t

T n2 + f n+12 ∆t

·

·

T (xN , tN )

• Lösung durch Thomas-Algorithmus (spezielle Form der Gaußschen Elimination)

• implizites Verfahren ist aufwendiger, aber bessere Stabilitätseigenschaften (größeres

∆t möglich)





4.4 Modifizierte Gleichung, numerische Diffusion/Dispersion

Wegen T.E. liefert FDA nicht die exakte Lösung der PDE. Man kann sich fragen, welchePDE-Lösung an den Gitterpunkten exakt mit der FDA übereinstimmt. Diese PDE ist die

sogenannte modifizierte Gleichung.

Diskussion anhand der 1D-linearen Advektionsgleichung

∂u

∂t + c

∂u

∂x = 0 Lösung u(x − c t)

FDA mit “1st order upwind method“ (Aufwind-Methode 1. Ordnung)

∂u

∂t

xi,tn

= un+1

i − uni

∆t

∂u

∂x

xi,tn

= un

i − uni−1

∆x

→

un+1i − un

i

∆t + c

uni − un

i−1

∆x = 0

Frage: wie kommt man zur modifizierten Gleichung? (Taylorentwicklung und Differenzie-

ren über Kreuz)

un+1i = un

i + ∂u

∂t ∆t +

1

2

∂ 2u

∂t2 ∆t2 +

1

6

∂ 3u

∂t3 ∆t3 + ...

uni−1 = uni − ∂u

∂x ∆x + 1

2

∂ 2u

∂x2 ∆x2 − 1

6

∂ 3u

∂x3 ∆x3 + ...

einsetzen in FDA ergibt folgendes

I ) ∂u

∂t + c

∂u

∂x = −

1

2

∂ 2u

∂t2 ∆t −

1

6

∂ 3u

∂t3 ∆t2 + ... +

c

2

∂ 2u

∂x2 ∆x −

c

6

∂ 3u

∂x3 ∆x2 + ...

II ) ∂

∂t(I ) →

∂ 2u

∂t2 + c

∂ 2u

∂x ∂t = −

1

2

∂ 3u

∂t3 ∆t + ...





III ) c ∂

∂x(I ) → c

∂ 2u

∂t ∂x + c2 ∂ 2u

∂x2 = −

c

2

∂ 3u

∂t2 ∂x ∆t + ...

II - III ergibt IV

IV ) ∂ 2u

∂t2 = c2 ∂ 2u

∂x2 −

1

2

∂ 3u

∂t3 ∆t +

c

2

∂ 3u

∂t2 ∂x ∆t +

c

2

∂ 3u

∂t ∂x2 ∆x −

c2

2

∂ 3u

∂x3 ∆x + ...

durch einsetzen von IV in I entsteht

V )

∂u

∂t + c

∂u

∂x = −

∆t2

2 c

2 ∂ 2u

∂x2 −

∆t

2

∂ 3u

∂t3 +

c

2

∂ 3u

∂x2 ∂t ∆x + ...V I )

∂ 3u

∂t3 →

∂

∂t(IV ) →

∂ 3u

∂t3 = c2 ∂ 3u

∂x2 ∂t + ...

Durch einsetzen von Gleichung VI in V wird ∂ 3u∂t3

eliminiert.

Des Weiteren werden ∂ 3u∂t2 ∂x

und ∂ 3u∂t ∂x2

durch die Gleichungen II und III eliminiert.

Am Ende resultiert aus V

∂u

∂t + c

∂u

∂x =

c ∆x

2 (1 − ν )

numerische Diffusion

∂ 2u

∂x2 +

c ∆x2

6 (3 ν − 2 ν 2 − 1)

numerische Dispersion

∂ 3u

∂x3

+ 0(∆t3, ∆x ∆t2, ∆t ∆x2, ∆x3)

mit ν = c ∆t∆x





4.5 Stabilität von FDA

Stabilitätsbegriff im Sinne der numerischen Mathematik erfordert eine Norm zur Charak-terisierung des Abstands zwischen numerischen Lösungen.

→ Maximum-Norm ||u||∞ = maxo<j<∞

|u j|

→ diskrete L2-Norm ||u||2 =

N j=0

(∆x j u j 2 )

Stabilitätsbegriff betrifft Anfangswertproblem, d.h. Zeitschrittverfahren.

→ Wir betrachten im Folgenden nur lineare Probleme

∂ tu = L u L = linear mit räumlicher Ableitung

L (u1 + u2) = L u1 + L u2 L (α u) = α L u

Eine FDA der PDE ∂ tu = L u sei gegeben; u, v seien numerische Lösungen dieser FDA

für verschiedene Anfangswerte u∗, v∗ (Randwerte bleiben gleich).

FDA ist stabil für n ∆t < tmax, wenn es eine Konstante k gibt, so dass ||vn − un|| ≤

k ||u∗ − v+|| gilt. Wobei die Konstante k unabhängig von u∗, v∗ und vom "Verfeinerungs-

weg"(refinement path) → ∆t(∆x) ist.

→ Bemerkung

Aus praktischer Sicht kann man grob sagen, dass sich durch die FDA Rundungsfehler nicht„aufschaukeln“ dürfen. Der mathematische Stabilitätsbegriff ist aber unabhängig davon,

ob man die FDA nur mit einem Rundungsfehler lösen kann, oder nicht. Beweis der Stabi-

litätseigenschaft im Allgemeinen schwierig oder nicht durchführbar. Für einfache lineare

PDE mit konstanten Koeffizienten gibt es die sogenannte Von-Neumann-Stabilitätsanalye.





4.6 Von-Neumann-Analyse in einer Raumdimension

n

j = u

n

j − v

n

j = Abstand der Lösung u, v

un j , vn

j erfüllen FDA für verschiedene Anfangswerten mit gleichen Randwerten. erfüllt

FDA mit Randwerten = 0. kann per Fourier-Ansatz (*) dargestellt werden.

n j = λn ei k xj (∗) x j = j ∆x k = Wellenzahl

λ = zeitlicher Verstärkungsfaktor

λ ergibt sich durch Einsetzen in FDA

→ λ(k, ∆x, ∆t), Verfahren ist Stabil für |λ(k, ∆x, ∆t) ≤ 0|

→ Bemerkung

ist streng eine Superposition von Fourier-Moden (*) mit verschiedenen k. Wegen Linea-

rität können Moden getrennt betrachtet werden.

Beispiel – FTCS für Wärmeleitungsgleichung (explizit)

∂T

∂t = κ

∂ 2T

∂x2 →

T n+1 j − T n j

∆t = κ

T n j+1 − 2 T n j + T n j−1

∆x2

T n j → n j →

n+1 j − n

j

∆t =

κ

∆x2 (n

j+1 − 2 n j + n

j−1)

Ansatz für

n

j

λn ei k xjλ − 1

∆t =

κ

∆x2 λn ei k xj

ei k ∆x − 2 + e−i k ∆x

Eulerformel

Eulerformel = e±i ϕ = cos(ϕ) ± i sin(ϕ)

ei k ∆x − 2 + e−i k ∆x → cos(k ∆x) + i sin(k ∆x) − 2 + cos(k ∆x) − i sin(k ∆x)





durch Kürzen erhält man somit

λ = 1 +

2 κ ∆t

∆x2 (cos(k ∆x) − 1) ≤0

→ |λ| ≤ 1, diese Bedingung erfordert, dass λ ≥ −1

worst case cos(k ∆x) = −1 → λ = 1 − 4 κ ∆t

∆x2 ≥ −1

→ 4 κ ∆t

∆x2 ≤ 2 →

κ ∆t

∆x2 ≤

1

2 → Stabilitätsbedingung

Einschränkung der Zeitschrittweite ∆t in Abhängigkeit von ∆x, typisch für explizite FDA.

Beispiel – Rückwärtsdifferenz der Wärmeleitungsgleichung in der Zeit (implizit)

∂T

∂t = κ

∂ 2T

∂x2 →

T n+1 j − T n j

∆t = κ

T n+1 j+1 − 2 T n+1

j + T n+1 j−1

∆x2

T n j → n j →

λ − 1

∆t =

2 λ κ

∆x2 (cos(k ∆x) − 1)

→ λ =

1

1 + 2 κ ∆t∆x2

(1 − cos(k ∆x)) ≥0

≤ 1

→ Uneingeschränkt stabil für ∆t beliebig

Vergleich explizite mit implizitem Verfahren

r =

κ ∆t

∆x2 → λ = λ(r, cos(k ∆x))





Frage: Was wäre exaktes Resultat für λ aus der PDE?

= A eα t

ei k x

→ einsetzen in PDE →

∂

∂t = κ

∂ 2

∂x2 → α κ k2

(tn+1)

(tn) = λexakt = e−κ k2 ∆t → unabhängig von ∆x

→ Zunehmende Abweichung zwischen λexakt und λ aus FDA, wenn k ∆x nicht klein genug

ist (und r ebenfalls).

4.7 Crank-Nicolson-FDA der Wärmeleitungsgleichung

∂T

∂t = κ

∂ 2T

∂x2

→ T n+1

i − T ni∆t

= κ

2 ∆x2

(T n+1

i+1 + T ni+1) − 2 (T n+1i + T ni ) + (T n+1

i−1 + T ni−1)

Fehleranalyse zeigt, dass die FDA einen T.E. 0(∆x2

, ∆t2

) besitzt. Lösung mittels einfachenGleichungssystems.





4.8 Methoden für die eindimensionale Advektionsgleichung

∂u∂t

+ c ∂u∂x

= 0 u = u(x − c t) c > 0

tn = n ∆t

xi = i ∆x

explizites Verfahren un+1

i − uni

∆t = −c

uni −uni−1∆x

= upwind

uni+1−uni∆x

= downwind

uni+1−uni−12 ∆x

= zentral

implizites Verfahren un+1

i − uni

∆t = − c

un+1i+1 − un+1

i−1

2 ∆x = zentral

→ upwind, downwind ebenfalls möglich

Von-Neumann-Analyse n j = un

j − vn j mit n

j = λn ei k xj

nach einsetzen in die verschiedenen FDA stellt man fest, dass

explizit

upwind → stabil für |λ| ≤ 1, falls ν ≤ 1

downwind

zentral

immer stabil

implizites Verfahren ist uneingeschränkt stabil

ν ist die sogenannte Courant-Zahl bzw. auch CFL (Courant-Friedrich-Lewy)

→ Diese FDA sind auch mit hoher numerischer Dissipation behaftet (ungünstig), außer-

dem haben sie Phasenfehler.





exakte Lösung der Fouriermode ist ei k (x−c t)

λk,exakt

= ei k c ∆x

λk,FDA

λk,exakt= r ei ϕ

λk,FDA = λk,exakt

ϕ = Phasendifferenz bzw. -fehler

Verfahren zweiter Ordnung besser wegen geringeren Abbruchfehlers

4.8.1 Beispiel – Leapfrog-Verfahren

zentrale Differenzen für ∂u

∂t und

∂u

∂x an xi und tn

un+1i − un−1

i

2 ∆t = −c

uni+1 − un

i−1

2 ∆x

T.E. ist 0(∆x2, ∆t2)

Von-Neumann-Analyse liefert

|λ| =

± 1 − ν 2 sin(k ∆x)

12 − i ν sin(k ∆x)

→ instabil für |λ| > 1, wenn ν > 1, für ν < 1 immer stabil

→ |λk| ≡ 1, wenn ν ≤ 1 → keine Dämpfung, nur Dispersion

Man sieht dies auch an der modifizierten Gleichung

∂u

∂t + c

∂u

∂x =

c ∆x2

6 (ν 2 − 1)

∂ 3u

∂x3 −

c ∆x4

120 (9 ν 4 − 10 ν 2 + 1)

∂ 5u

∂x5 + ...

(nur ungerade Ableitungen rechts)





4.8.2 Beispiel – Lax-Wendroff-Verfahren

Idee = Taylor-Entwicklung in der Zeit

u(x, t + ∆t) = u(x, t) + ∂u

∂t =−c ∂u

∂x

∆t + 1

2

∂ 2u

∂t2 ∆t2

→ ∂ 2u

∂t2 =

∂

∂t

∂u

∂t

=

∂

∂t

−c

∂u

∂x

= −c

∂

∂x

∂u

∂t

= −c

∂

∂x

−c

∂u

∂x

= c2 ∂ 2u

∂x2

→ u(x, t + ∆t) = u(x, t) − c ∂u

∂x ∆t +

1

2 c2 ∂ 2u

∂x2 ∆t2

→ un+1

i − uni

∆t = −

c

2 ∆x (un

i+1 − uni−1) +

c2 ∆t

2 ∆x2 (un

i+1 − 2 uni + un

i−1)

T.E. ist 0(∆x2, ∆t2)

Lax-Wendroff ist „stabilisierte Version“ des instabilen Zentralen-Differenzen-Verfahren

1.Ordnung. Lax-Wendroff ist stabil für ν ≤ 1

modifizierte Gleichung für Lax-Wendroff

∂u∂t

+ c ∂u∂x

= −c ∆x2

6 (1 − ν 2) ∂

3

u∂x3

− c ∆x3

8 ν (1 − ν 2) ∂

4

u∂x4

→ Dominierend ist der Dispersionsterm mit ∂ 3u∂x3





4.8.3 Beispiel – MacCormack-Verfahren

basiert auf Prädiktor-Korrektor-Idee

un+1i − un

i

∆t =

1

2

∂u

∂t

t

+ ∂u

∂t

t+∆t

unbekannt

entspricht der Approximation bei t + 1

2 ∆t

Prädiktor für t + ∆t, u*

∂u

∂t t

= u∗i − un

i

∆t = −c

uni+1 − un

i

∆x → u∗i (rechtsseitige Differenz)

Korrektor

∂u

∂t

t+∆t

= −c u∗i − u∗i−1

∆x (linksseitige Differenz)

Zusammengefasst

→ un+1

i − uni

∆t =

1

2

u∗i − un

i

∆t −

c

∆x (u∗i − u∗i−1)

Prädiktorschritt

→ u∗i = uni − ν (un

i+1 − uni )

Korrektorschritt

→ un+1i =

1

2 (un

i − u∗i ) − ν

2 (u∗i − u∗i−1)

stabil für ν ≤ 1 und T.E. ist 0(∆x2, ∆t2)





4.8.4 Stabilität für Burgers-Gleichung

∂u

∂t + u

∂u

∂x = µ

∂ 2u

∂x2 (nichtlinear) → Von-Neumann-Analyse nicht anwendbar

FDA sei z.B.:

→ F DA sei z.B. : un+1

i − uni

∆t + un

i =C

uni+1 − un

i−1

2 ∆x = µ

uni+1 − 2 un

i + uni−1

Deltax2

Konstante C für Stabilitätsanalyse eingeführt

→ wenn nichtlinearer Term uni durch Konstante C ersetzt wird, dann ist Von-Neumann

Analyse möglich.

C 2 ∆t ≤ 2 µ → µ ∆t

∆x2 ≤

1

2 C = umax (heuristisches Vorgehen)

Für nichtlineare Gleichungen wird häufig heuristische Stabilität für einzelne Terme separatgefordert, z.B.: Advektionsterm und Diffusionsterm, man erwartet, dass das Verfahren

funktioniert.





4.9 Verfahren für zweidimensionale Diffusionsprobleme

∂T

∂t = κ

∂ 2T

∂x2 +

∂ 2T

∂y 2

xi = i ∆x, y j = j ∆y

T ni,j = T (x,y,t)

Stabilität kann mit von-Neumann-Analyse untersucht werden.

Beispiel – FTCS für Wärmeleitungsgleichung (explizit)

∂T

∂t = κ

∂ 2T

∂x2 +

∂ 2T

∂y 2

→ T n+1

i,j − T ni,j

∆t = κ

T ni+1,j − 2 T ni,j + T ni−1,j

∆x2 +

T ni,j+1 − 2 T ni,j + T ni,j−1

∆y2

ni,j = λn ei kx xi ei ky yj stabil für → κ ∆t

1

∆x2 +

1

∆y2

≤

1

2

Implizites Verfahren oder Crank-Nicolson-Verfahren führen auf lineares Gleichungssystemmit Nx · Ny Unbekannten, sehr aufwendig zu berechnen (Nx · Ny)3 Rechenschritte not-

wenig.





4.9.1 Alternative mit geringerem Aufwand = „approximate factorization“

vereinfachte Schreibweise

I ) Lx T ni,j = T ni+1,j − 2 T ni,j + T ni−1,j

∆x2

II ) Ly T ni,j = T ni,j+1 − 2 T ni,j + T ni−,j−1

∆y2

III ) ∆T i,j = T n+1i,j − T ni,j

implizite Methode 1∆t

∆T i,j = κ

Lx T n+1

i,j + Ly T n+1i,j

daraus ergibt sich 1

∆t ∆T i,j = κ

Lx T ni,j + Lx ∆T i,j + Ly T ni,j + Ly ∆T i,j

|∆t

→ (1 − κ ∆t Lx − κ ∆t Ly) ≈(1−κ ∆t Lx) (1−κ ∆t Ly)

∆T i,j = κ ∆t (Lx T ni,j + Ly T ni,j)

Fehler ist hierbei κ2 ∆t2 Lx Ly ∆T i,j d.h. T.E. ist 0(∆t2)

ersetzen von (1 − κ ∆t Ly) ∆T i,j = ∆ T i,j

durch Faktorisierung erhält man somit

→ (1 − κ ∆t Lx) ∆

T i,j = κ ∆t (Lx T ni,j + Ly T ni,j) 1D-Problem in i

→ (1 − κ ∆t Ly) ∆T i,j = ∆ T i,j 1D-Problem in j

Verfahren ist uneingeschränkt stabil, Aufwand ist deutlich geringer.





4.9.2 Weitere Alternative = ADI (alternating directions implicit)

Erklärung An 2D Diffusionsproblem auf Recheck∂T

∂t = κ

∂ 2T

∂x2 +

∂ 2T

∂y 2

1. Halbschritt T ni,j → T i,j

→ T i,j − T ni,j

∆t2

= κ

T i+1,j − 2

T i,j +

T i−1,j

∆x2 implizit

+ T ni,j+1 − 2 T ni,j + T ni,j−1

∆y2 explizit

2. Halbschritt T i,j → T n+1i,j

→ T n+1

i,j − T i,j

∆t2

= κ

T i+1,j − 2 T i,j + T i−1,j

∆x2

explizit

+ T n+1

i,j+1 − 2 T n+1i,j + T n+1

i,j−1

∆y2

implizit

daraus ergibt sich1 −

κ ∆t

2 Lx

T i,j =

1 +

κ ∆t

2 Ly

T ni,j 1D in x

1 −

κ ∆t

2 Ly

T n+1

i,j =

1 +

κ ∆t

2 Lx

T i,j 1D in y

Abbruchfehler ist 0(∆t2, ∆x2, ∆y2), ist uneingeschränkt stabil in 2D, gilt nicht in 3D,

dort nur approximate factorization ohne Einschränkung stabil.





4.10 Methoden für elliptische Probleme

elliptisch bedeutet Gleichgewicht, keine zeitartige Variable

→ einfachstes Problem = 1D stationäre Wärmeleitung

κ ∂ 2u

∂x2 = f

auf 0 < x < L

xi = i ∆x

∆x = L

N

zentrale Differenzen κ ui

+1−2 ui+ui−

1∆x2 = f i, lineares Gleichungssystem, tridiagonal

→ Lösung mit Thomas-Algorithmus, Aufwand dafür 0(N) Rechenoperationen

→ 2D Problem = ∇2 u = f (Poisson-Gleichung)

zentrale Differenzen

→

ui+1,j − 2 ui,j + ui−1,j

∆x2 +

ui,j+1 − 2 ui,j + ui,j−1

∆y2 = f i,j

Man erhält wieder ein lineares Gleichungssystem. Die Struktur ist aber komplizierter als

im eindimensionalen Fall.

→ Betrachten rechteckiges Gebiet als Beispiel

0 ≤ i ≤ N x 0 ≤ j ≤ N y

L = j (N x + 1) + i für Nummerierung

→ υL = ui,j, cL = f i,j

Fortlaufende Nummerierung der einzelnen Punkte, um lineares Gleichungssystem in üb-licher Weise zu erhalten.





Überführung der Gleichung für zentrale Differenzen in L-Schreibweise

i, j → j (N x

+ 1) + i = L

i + 1, j → j (N x + 1) + i + 1 = L + 1

i − 1, j → j (N x + 1) + i − 1 = L − 1

i, j + 1 → ( j + 1) (N x + 1) + i = L + N x + 1

i, j − 1 → ( j − 1) (N x + 1) + i = L − N x − 1

→ aL,L υL + aL,L+1 υL+1 + aL,L−1 υL−1 + aL,L+N x+1 υL+N x+1 + aL,L−N x−1 υL−N x−1 = cL

gilt nicht am Rand

Matrixschreibweise A υ = c A = (aL,m) υ = (υL) c = (cL)

Koeffizientenvergleich ergibt

aL,L = −2 1

∆x2

+ 1

∆y2 aL,L+1 = aL,L−1 = 1

∆x2

aL,L+N X

+1 = aL,L−N x−1 = 1

∆y2

aL,k = 0 = alle restlichen Werte (siehe Matrixschreibweise)

Struktur der Matrix (Bandmatrix, schwach besetzt (sparse)





Für solche Matrizen gibt es spezielle direkte Lösungsverfahren. Direkt heißt, man erhält

die exakte Lösung nach endlich vielen Rechenschritten, z.B.: Gaußsche Elimination.

Allgemeine direkte Verfahren, die die Bandstruktur nicht beachten, sind sehr aufwendigund deswegen im Allgemeinen ungeeignet. Der Aufwand der Lösung geht mit 0(N3)

N = Anzahl der Unbekannten

z.B.: N x = N y = 100 → N = 104 → 0(N 3) = 1012 Rechenoperationen

→ Alternative zu direkten Verfahren = iterative Verfahren

Iterativ heißt, keine exakte Lösung nach endlich vielen Rechenoperationen. Iterative Verfah-

ren sind häufig vorteilhaft (schneller + weniger Speicherintensiv).

Idee für Iteration

A υ = c A = A − B + B → B υ + (A − B) υ = c

→ B υn+1 = c + (B − A) υn n = Index der Itertion

Auflösung nach υn+1 soll möglichst einfach sein, rasche Konvergenz ist erwünscht. Ein-

fachste Wahl, B sei diagonaler Anteil von A, dies ist das Jacobi-Verfahren.

N j=0

aL,j υ j = cL = aL,L υn+1L +

L−1 j=0

aL,j υ n j +

N j=L+1

aL,j υ n j

→ υn+1L =

1

aL,L

cL −

L−1 j=0

aL,j υ n j −

N j=L+1

aL,j υ n j

Verbesserte Variante, man nutzt die schon bekannten neuen Werte υn+1

j auf der rechten

Seite. Dies ist das Gauß-Seidel-Verfahren.

→ υn+1L =

1

aL,L

cL −L−1 j=0

aL,j υ n+1 j −

N j=L+1

aL,j υ n j





In der Matrixnotation A = D + L + U

D = diagonaler Anteil

L = untere Dreiecksmatrix

U = obere Dreiecksmatrix

→ Jacob-Verfahren B = D

D υn+1 = c − (L + U ) υn

→ Gauß-Seidel-Verfahren

D υn+1 = c − L υn+1 − U υn

→ Gauß-Seidel-Verfahren mit Gewichtung (Diagonalanteil vom vorigen Schritt)

D υn+1 = (1 − ω) D υn + ω (c − L υn+1 − U υn)

ω =Relaxationsparameter

ω < 1 = Unterrelaxation

ω < 1 = Überrelaxation

Ziel ist es ω so zu wählen, dass Konvergenz möglich rasch erfolgt. (SOR successive over-

reaxation)

→ Konvergenzkriterien

δ n = ||υn+1 − υn|| =

N l=0

(υn+1l − υn

l )2

12

Falls ν n → ν konvergiert, dann geht δ n → 0. Abbruch der Iteration, wenn δ n < .

Frage: wie klein muss sein?





weitere Möglichkeit – Residuum

r

n

= c − A υ

n

c = A υ

→ rn = A υ − A υn = A (υ − υn)

→ ||rn|| = ||A|| ||(υ − υn)|| zur Definition der Matrixnorm, siehe lienare Algebra

Abbruch, wenn ||rn|| < , keine Garantie für ||υ − υn|| < Zielwert, aber trotzdem vorteil-

haft, weil Zusammenhang zwischen rn und υ − υn.

Frage: Wann konvergiert die iterative Verfahren Jacobi, Gauß-Seidel, SOR?

Hinreichend aber nicht zwingend notwendig ist Diagonaldominanz der Matrix A. Diago-

naldominanz liegt vor, wenn |aL,L| ≥N

j=1,j=L |aL,j| für alle L mit 1 ≤ L ≤ N .

Für mindestens einen Wert von L muss außerdem streng „>“ gelten.

Beispiel für Konvergenzraten des Jacobi-Verfahrens

∇2u = f auf 0 < x < 1, 0 < y < 1 u = 0 am Rand

→ ∆x = ∆y N ∆x = 1 (äquidistantes Gitter)

→ ui+1,j − 2 ui,j + ui−1,j

∆x2 +

ui,j+1 − 2 ui,j + ui,j−1

∆y2 = f i,j

für Jacobi gilt

− 4∆x2 u(n+1)i,j + 1∆x2 u(n)i+1,j + (u(n)i−1,j + (u(n)i,j+1 + (u(n)i,j−1 = f i,j ∆x = ∆y

Fehler = e(n)i,j = u

(n)i,j − u

(exakt)i,j exakte Lösung erfüllt die Gleichung

u(n)i,j = e

(n)i,j + u

(exakt)i,j in Jacobi einsetzen → e

(n+1)i,j = 1

4

e

(n)i+1,j + e

(n)i−1,j + e

(n)i,j+1 + e

(n)i,j−1





Fehler kann in Fourier-Moden entwickelt werden

e

(n)

i,j = kx,ky a

(n)

kx,ky sin(kx xi) sin(ky y j) Sinusansatz wegen Randbedingung an u

kx = π, 2π, ..., (N − 1)π = Wellenzahlen in x-Richtung

ky = π, 2π, ..., (N − 1)π = Wellenzahlen in y-Richtung

a(n)kx,ky

= Amplitude

Einsetzen der Fourierdarstellung ergibt Zusammenhang zwischen a(n)kx,ky und a

(n+1)kx,ky

e(n+1)i,j =

1

4

kx,ky

a(n)kx,ky

sin(kx xi − kx ∆x) sin(ky y j) + sin(kx xi + kx ∆x) sin(ky y j )

+sin(kx xi) sin(ky y j − ky ∆x) + sin(kx xi) sin(ky y j + ky ∆x)

mit sin(α + β ) = sin(α) cos(β ) + cos(α) sin(β )

und sin(α − β ) = sin(α) cos(β ) − cos(α) sin(β )

e(n+1)i,j =

1

4

kx,ky

a(n)kx,ky

[sin(kx xi) sin(ky y j) (2 cos(kx ∆x) + 2 cos(ky ∆x))]

→ e(n+1)i,j =

kx,ky

a(n+1)kx,ky

[sin(kx xi) sin(ky y j ) ]

Koeffizienten müssen gleich sein, weil Fourier-Moden linear unabhängig sind.

a(n+1)kx,ky

= 1

2 (cos(kx ∆x) + cos(ky ∆x))

µ

a(n)kx,ky

|µ| < 1, sonst Konvergenz





Betragsmäßig größte Werte von µ bestimmen das Abklingen des Fehlers mit Iterationsin-

dex n.

∆x = 1

N kx ∆x =

π

N , ...,

N − 1

N π ky = analog

→ betragsmäßig größte Werte von µ bei πN

und N −1N

π (siehe Kosinusfunktion)

max|µ| = cos

π

N

≈ 1 −

1

2

π

N

2

Konvergenzratee(n+1)

i,j

< max|µ|e(n)

i,j

asymptotisch (n → ∞)

e(n)i,j

∼ λ−n λ = max|µ| = 1 − 12

πN

2

Frage: Wie lange dauert es, den Fehler um einen bestimmten Faktor zu verringern?

λn = C C < 1

→ n ln(λ) = ln(C ) es gilt ln(1 + x) = x → ln(λ) ≈ −1

2

π

N

2

→ n = ln(λ)

ln(C ) = 2 N 2

ln(C )

π2

Mit steigenden N sehr rasche Zunahme der erforderlichen Iterationen

Problem – langwellige Anteile

Verbesserung der Konvergenz durch Mehrgitterverfahren (Multigrid)




5 FINITE-DIFFERENZEN-VERFAHREN/STRÖMUNGSPROBLEME 55

5 Finite-Differenzen-Verfahren/Strömungsprobleme

5.1 Stromfunktions-Wirbelstärke-Formulierung

inkompressible, reibungsbehaftete Strömung

Kontinuität ∇ · u = 0

Impulserhaltung ∂u

∂t + (u · ∇)u = −

1

∇ p + ν ∇2u

→ Problem = Druck und Inkompressibilitätsbedingung

Ausweg für 2D-Strömungen = Stromfunktion

→ ∂ux

∂x +

∂uy

∂y = 0

Ansatz

ux =

∂ψ

∂y und uy = −

∂ψ

∂x →

∂ 2ψ

∂x ∂y −

∂ 2ψ

∂x ∂y = 0

Wirbelstärke ω = ∇ × u =

0

0∂uy

∂x − ∂ux

∂y

→ ω = ω ez = (∂ x uy − ∂ y ux) ez = −∇2ψ

Rotation der Impulsbilanz um Druck zu eliminieren

x-Komponente

∂u

∂t + u

∂u

∂x + v

∂u

∂y = −

1

∂p

∂x + ν ∇2

u ∂

∂y ()

y-Komponente ∂v

∂t + u

∂v

∂x + v

∂v

∂y = −

1

∂p

∂y + ν ∇2v

∂

∂x()

→ ∂

∂y(x-Komp.) −

∂

∂x(y-Komp.) liefert

∂ω

∂t + u

∂ω

∂x + v

∂ω

∂y = ν ∇2ω

= Wirbeltransportgleichung mit ω = −∇2ψ, u = ∂ yψ und y = −∂ xψ

Für ψ werden Randbedingungen benötigt, ebenso für ω.





Beispiel - Lid-driven Cavity flow

ψ ist konstant auf Rand

→ geschlossene Strömung im Rechengebiet

v = 0 = −∂ xψ → ψ = konstant

u = 0 = ∂ yψ → ψ = konstant

Diskretisierung auf äquidistantem Gitter

1) ωn+1i,j − ωni,j

∆t = − un

i,jωni+1,j − ωni−1,j

2 ∆x − vn

i,jωni,j+1 − ωni,j−1

2 ∆y

+ ν

ωn

i+1,j − ωni,j + ωn

i−1,j

∆x2 +

ωni,j+1 − ωn

i,j + ωni,j−1

∆y2

2) ψn+1

i+1,j − 2 ψn+1i,j + ψn+1

i−1,j

∆x2 +

ψn+1i,j+1 − 2 ψn+1

i,j + ψn+1i,j−1

∆y2 = −ωn+1

i,j

3) un+1i,j =

ψn+1i,j+1 − ψn+1

i,j−1

2 ∆y vn+1

i,j = ψn+1

i+1,j − ψn+1i−1,j

2 ∆x

Man löst sukzessiv Gleichung 1, 2 und 3. Randwerte für Gleichung 2 sind bekannt. ωn+1i,j

am Rand für Gleichung 1 ergibt sich aus Gleichung 3.





5.2 Rolle des Drucks in inkompressiblen Strömungen

∂u

∂t + (u · ∇)u = −

1

∇ p + ν ∇2

u +

1

f

∇ · u

Druck erfüllt keine Evolutionsgleichung. Die Inkompressibilität ∇ · u stellt eine Nebenbe-

dingung dar, die durch geeignete Wahl eines Drucks erfüllt werden muss.

Bestimmungsgleichung für Druck aus Divergenz der Impulsbilanz:

∂

∂t(∇ · u) + ∇(u · ∇)u = −

1

∇2 p + ν ∇2(∇ · u) +

1

∇ · f

→ ∇2 p = ∇ · f − ∇(u · ∇)u = Poisson-Gleichung für p

Elliptisches Problem für p auch bei Verwendung expliziter Zeitschrittverfahren → hoher

Aufwand

Problem:

Wie kombiniert man die Bestimmung des Drucks mit Zeitschrittverfahren für die Ge-

schwindigkeitskomponenten?

Bemerkung:

Bei kompressiblen Strömungen stellt sich das Problem nicht. In diesem Fall existiert Evo-

lutionsgleichung für Dichte ρ und Temperatur T (Massen- und Energiegleichungen). DerDruck p wird über Zustandsgleichungen aus Dichte und Temperatur bestimmt.





5.3 Projektionsverfahren für zeitabhängige inkompressible Strö-

mungen

Idee basiert auf Zerlegung (splitting) der rechten Seite der Impulsbilanz. Analogie zu ADI

bzw. Prädiktor-Korrektor-Idee.

Zeitschritt un → un+1 in zwei Teilschritte zerlegt

1. Teilschritt u∗ aus Impulsgleichung ohne ∇ p

2. Teilschritt un+1 aus u∗ mittels ∇ p ohne restliche Terme der Impulsbilanz

∇ p im zweiten Teilschritt so wählen, dass ∇ · un+1 = 0 gilt.

∂u

∂t = −

∇ p

− N (u, u) + ν ∇2u

= F

N (u, u) = ∂

∂xi(ui ui)

→ 1. Teilschritt u∗ − un

∆t = F n

→ 2. Teilschritt un+1 − u∗

∆t = −

1

∇ pn+1 (∗)

Addition beider Teilschritte gibt konsistente Approximation der Impulsgleichung.

∇ · un+1 = 0 → ∇ · (∗) ergibt 1

∆t ∇ · un+1

=0

− 1

∆t ∇ · u∗ = −

1

∇2 pn+1

→ ∇2 pn+1 =

∆t ∇ · u∗ (∗∗)

Bestimmung von p erfordert Randbedingungen. Auf dem Strömungsrand gilt u · n = 0.





Anwenden auf (*) ergibt:

n ·un+1 − u∗

∆t Rand

= n 1

∇ pn+1

Rand

→ n · un+1 =0

− n · u∗ = ∆t

n ∇ pn+1

=∂pn+1

∂n

→ ∂pn+1

∂n =

∆t n · u∗ = Randbedingung für (**)

→ pn+1, aus pn+1 folgt un+1 aus (∗)

Räumliche Diskretisierung zweiter Ordnung erfordert sogenannte versetzte Gitter (stag-

gered grid). Ansonsten Problem mit Entkopplung von Untergittern.

→ Erläuterung des Entkopplungsproblems:

Verwenden zentrale Differenzen und u, v, p an ein und denselben Gitterpunkten.

un+1 = u∗ − ∆t

∂ x p vn+1 = v∗ −

∆t

∂ y p

→ un+1i,j = u∗i,j −

∆t

pn+1i+1,j − pn+1

i−1,j

2 ∆x → vn+1

i,j = v∗i,j − ∆t

pn+1i,j+1 − pn+1

i,j−1

2 ∆y

→ ∇ · un+1 = 0 = un+1

i+1,j − un+1i−1,j

2 ∆x +

vn+1i,j+1 − vn+1

i,j−1

2 ∆y

Diskrete Form der Poisson-Gleichung wird dann

pn+1i+2,j − 2 pn+1

i,j + pn+1i−2,j

(2 ∆x)2 +

pn+1i,j+2 − 2 pn+1

i,j + pn+1i,j−2

(2 ∆y)2 =

∆t

u∗i+1,j − u∗i−1,j

2 ∆x +

v∗i,j+1 − v∗i,j−1

2 ∆y

= (∗ ∗ ∗)

Druckwerte an i, i+2, i-2, j, j+2, j-2 sind verknüpft mit Geschwindigkeitswerten bei i+1,

i-1, j+1, j-1 → Entkopplung





Das ist ein Problem, weil man stark oszillierende Geschwindigkeitsfelder haben kann, die

bei lokaler Auswertung eigentlich nicht inkompressibel sind. Diese werden durch Berech-

nung von pn+1

mittels (***) nicht korrigiert.

u∗i+1,j − u∗i−1,j = 0

v∗i,j+1 − v∗i,j−1 = 0

Vermeidung des Problems bei versetzten Gittern.

→ zum Begriff Projektionsverfahren

Helmholtzsche Zerlegung eines Vektorfelds a auf Gebiet Ω mit Randwerten a|∂ Ω = 0

a = as + ∇φ (Zerlegung)

Wobei as divergenzfrei ist (solenoidal) und aS |∂ Ω = 0. Die Zerlegung in eine Summe aus

divergenzfreien Vektorfeld und einem Gradientenfeld ist eindeutig. Die Felder as bildeneinen linearen Raum (Vektorraum), die Felder ∇φ ebenfalls. Beide Räume sind orthogo-

nale Unterräume des linearen Raums aller Felder bezüglich der L2-Norm.

Das heißt:

Ω

aS ∇φ dΩ =

Ω∇ (aS · φ) dΩ −

Ω

φ ∇ · a =0

dΩ

Ω ∇ (aS · φ) dΩ = Ω φ n aS dA = 0 nach Gauß





Im Druckkorrekturschritt wird der divergenzbehaftete Anteil von u∗ entfernt. Das ent-

spricht genau einer Projektion auf den von den Feldern as aufgespannten Unterraum.

Formel un+1 = P (u∗)

P = Projektionsabbildung

P = P 2 definiert Projektionsabbildung

5.4 Iterative Lösung von stationären Strömungen

Ziel – Wollen stationäre Lösung der Navier-Stokes-Gleichung finden

∇ p

+ N (u, u) = ν ∇2u ∇ · u = 0 text(Konti)

Iteration = sukzessive Verbesserung einer Nährung um, pm

m → ∞ entspricht um → u und pm → p

Idee analog zum Projektionsverfahren, Verwendung der linearisierten Gleichung.

Schritt 1) um → u∗, pm bleibt fix

Impulsgleichung linearisiert → ∇ pm

+ N (u, u) = ν ∇2u∗ (∗)

→ kann nach u∗ aufgelöst werden

Schritt 2) um+1 = u∗ + u, pm+1 = pm + p

u und p sind Korrekturen, hängen über Impulsgleichung zusammen. u soll sichern, dass

∇ · um+1 = 0 ist.





→ setzen um+1 = u∗ + u und pm+1 = pm + p in linearisierte Impulsgleichung (*) ein.

∇ pm

+

∇ p

+ N (um, u∗) + N (um, u) = ν ∇2u∗ + ν ∇2u (∗∗)

∇ · (u∗ + u) = 0 (∗ ∗ ∗)

Wenn man Gleichung (**) diskretisiert vorstellt, dann gibt es Beträge von u am Git-

terpunkt und von Nachbarn (Rechenstern von ∇2). Vernachlässigen der Beiträge von

Nachbarpunkten.

→ ∇ p

= −C · u (heuristisches Vorgehen) C ist nicht dimensionslos

Einsetzen in Gleichung (***) liefert → ∇ ·

u∗ − 1C

∇ p

= 0, → ∇2 p = C

∇ · u∗ (∗ ∗ ∗ ∗)

Brauchen Randbedingungen für Gleichung (*) und (****)

Für (*) erfüllt u∗ die Randbedingungen von u

Für (****) gilt n · u = 0 = n ·

u∗ −

1

C

∇ p

→

∂p

∂n =

C

n · u∗

Für Druck macht man Unterrelaxation → pm+1 = pm + α p α << 1

Das ist (vereinfacht) die Idee des SIMPLE-Algorithmus (semi-implicit method for pressure-

linked equations)




6 FINITE-VOLUMEN-METHODE 63

6 Finite-Volumen-Methode

Finite-Volumen-Verfahren werden in kommerziellen Codes häufig verwendet, weil sie fürunstrukturierte Netze geeignet sind. Ein weiterer Vorteil ist, dass Finite-Volumen-Verfahren

bestimmte globale Erhaltungseigenschaften der PDE exakt sichern.

6.1 Integrale Form der Grundgleichungen

Bildet Ausgangspunkt für Finite-Volumen-Verfahren (z.B.: Advektions-Diffusionsgleichung

für Skalarfeld Φ).

∂ Φ

∂t + ∇ ( u Φ

Advektion

− β ∇Φ Diffusion

) = q Quellterm

(∗)

u Φ = Flussdichte aus Advektion von Φ

−β ∇Φ = Flussdichte aus Diffusion von Φ

Integration von (*) über ein raumfestes Volumen Ω liefert Ω

∂ Ω

∂t dV =

∂

∂t

Ω

Φ dV

Ω

∇ (u Φ − β ∇Φ) dV Gauß−−−→

∂ Ω

n (u Φ − β ∇Φ) dA

Ω

q dV

→ ∂

∂t

Ω

Φ dV = −

∂ Ωn u Φ dA

konvektiver Fluss

+

∂ Ωn β ∇Φ dA

diffusiver Fluss

+

Ωq dV

Quellen/Senken

(∗∗)





In dieser Form kann z.B. Massen- und Energiebilanz geschrieben werden

Φ = , β = 0, q = 0 für Massenbilanz

Φ = c p T , β = λ, q = 0 für Energiebilanz (λ = Wärmeleitfähigkeit)

Impulsbilanz ist etwas schwieriger. Diffusionsterm ist komplizierter und es gibt Druck.

x-Komponente ∂

∂t

ω

ux dV = −

∂ Ωn u ux dA −

∂ Ω

p dA +

∂ Ω(n

↔τ )x dA

viskoser Spannungstensor

Man wendet die integrale Form (**) der Erhaltungsgleichung auf Zellen des Gitters/Netzes

an. → Approximation von Integralen statt von Ableitungen wie bei Finite-Differenzen-

Verfahren.

6.2 globale Erhaltungseigenschaften

am eindimensionalen Beispiel ∂ Φ

∂t = ∇ F + q gilt auf 0 ≤ x ≤ L ∇ F = ∂F

∂x Integration liefert

→ L

0

∂ Φ

∂t dx =

L

0

∂F

∂x dx +

L

0q dx (∗)

mit L

0

∂ Φ

∂t dx =

∂

∂t L

0

Φ dx =M

und L

0

∂F

∂x dx = F (L) − F (0)

Die Größe M (Gesamtmenge von Φ) ist also zeitlich konstant. Eine Finite-Volumen-

Diskretisierung sichert, dass diese Eigenschaft auch für die numerische Approximation

gilt. Betrachten dazu Zerlegung in Intervalle (finites Strecken) Φi.





Finite-Volumen-Zerlegung (bzw. hier Intervall-Zerlegung, da 1D):

Mittelpunkt des des Intervalls Φi sei xi

Integrale Form →

Ωi dx anwenden auf Differentialgleichung

→ ∂

∂t

Ωi

Φ dx = ∂

∂t(Ωi hi) hi = xi+ 1

2− xi− 1

2= Länge von Ωi

→

Ωi

∂F

∂x dx = F

xi+12

xi−1

2

= F (xi+0,5) − F (xi−0,5)

F (xi+0,5) und F (xi−0,5) sind die numerischen Approximationen des Flusses F auf dem

rechten und linken Rand von Ωi. F (xi+0,5) und F (xi−0,5) müssen im Allgemeinen aus Φi

bestimmt werden, weil F zumeist von Φ abhängt.

Für Advektion z.B. F = −u Φ, für Diffusion F = D ∂ Φ∂x

.

→ ∂ ∂t

(Φi hi) = F (xi+0,5) − F (xi−0,5)

Bilden der Summe über alle Integrale

→

i

∂

∂t

Ωi

Φ dx =

i

F (xi+0,5) −

i

F (xi−0,5) +

i

Ωi

q dx

mit i

∂ ∂t

Ωi

Φ dx = ∂ ∂t

i

Ωi

Φ dx = ∂ ∂t

Ω

Φ dx

i

F (xi+0,5) −

i

F (xi−0,5) = F (L) − F (0)

i

Ωi

q dx =

Ωq dx

→ dies ist äquivalent zu (*)





Die Werte F(L) und F(0) sind gerade die vorgegebenen Randwerte. Weil F (xi+0,5) =

F (xi−0,5) für i = 1, . . . , L gilt, gilt somit F(L) - F(0) = 0.

→

i

∂

∂t

Ωi

Φ dx =

i

F (xi+0,5) −

i

F (xi−0,5) = 0

Die Eigenschaft ∂ ∂t

L0 Φ dx = 0 wird durch die räumliche Diskretisierung also erhalten.

Ursa-che ist, dass sich die Flussterme an den Rändern der Intervalle gegenseitig aufheben,

somit ist eine globale Erhaltung gewährleistet, z.B. für die Masse. Änderungen nur durch

äußeren Rand möglich, wo Randwerte vorgegeben sind. Das überträgt sich auch auf denzwei- und dreidimensionalen Fall.

6.3 1D-Beispiel für Finite-Volumen-Diskretisierung

∂T

∂t =

∂

∂x

κ

∂T

∂x

Ränder bei xi− 12 und xi+ 1

2 → xi+ 12 − xi− 1

2 = ∆x

∂

∂t

Ωi

T dx = κ ∂T

∂x

xi+1

2

xi− 1

2

T i sei Mittelwert Ωi

→ T n+1

i − T ni

∆t ∆x = κ

T ni+1 − T ni

∆x −

T ni − T ni−1

∆x

→ T n+1

i − T ni∆t

= κ

∆x2

T ni+1 − 2 T ni + T i−1

= FTCS





6.4 2D-Beispiel für Finite-Volumen-Diskretisierung

Ωi finite Fläche

N = Nord, S = Süd, W = West, O = Ost

∂ Φ

∂t = ∇ · F =

∂F x

∂x +

∂F y

∂y

→

Ωi

∂ Ω

∂t dA =

d

dt Ωi Ai Ai = Flächeninhalt von Ωi

→

Ωi

∇ · F dA =

∂ Ωi

n F ds ≈ F N LN + F S LS + F W LW + F O LO

L = Kantenlänge der Ränder, F = Mittelwerte des Flusses F n auf den Randstücken

Wenn man sich nun benachbarte Zellen vorstellt und summiert, heben sich die Flussterme

an den inneren Rändern wieder auf.d

dt

4i=1

Φi Ai = F N 1 LN

1 + F S 1 LS

1 + F W 1 LW

1 + F O1 LO1

+ F N 2 LN

2 + F S 2 LS

2 + F W 2 LW

2 + F O2 LO2

+ F N 3 LN

3 + F S 3 LS

3 + F W 3 LW

3 + F O3 LO3

+ F N 4 LN

4 + F S 4 LS

4 + F W 4 LW

4 + F O4 LO4

Es gilt: F S 2 LS

2 +F N 4 LN

4 = 0, F S 1 LS

1 +F N 3 LN

3 = 0, F O3 LO3 + F W

4 LW 4 = 0, F O1 LO

1 + F W 2 LW

2 = 0

Es bleiben als nur die Randterme außen.

d 4Φ A F N LN + F N LN + F O LO + F O LO + F S LS + F S LS + F W LW + F W LW

Numerische Strömungsmechanik

Documents

Transcript of Numerische Strömungsmechanik