Computational Physics Praktikum: Numerische Hydrodynamik Wilhelm Kley Institut fur ¨ Astronomie & Astrophysik & Kepler Center for Astro and Particle Physics Tubingen ¨

Hydrodynamik:

Hydrodynamische Gleichungen

Die Euler-Gleichungen der Hydrodynamik lauten in Erhaltungsform ∂ρ + ∇·(ρ~u ) ∂t ∂(ρ~u ) + ∇·(ρ~u ⊗ ~u ) ∂t ∂(ρ) + ∇·(ρ~u ) ∂t

= 0

(1)

= −∇p + ρ~k

(2)

= −p∇·~u

(3)

~u : Geschwindigkeit, ~k: außere ¨ ¨ Krafte,  innere spezifische Energie Die Gleichungen beschreiben die Erhaltung der Masse, Impuls und Energie. ¨ Vervollstandigung durch Zustandsgleichung: p = (γ − 1) ρ

(4)

Forme damit die Energie-Gleichung (3) in eine Gl. fur ¨ den Druck um ∂p + ∇·(p~u ) = −(γ − 1)p ∇ · ~u ∂t W. Kley

Computational Physics Praktikum

(5) 2

Hydrodynamik:

Umschreiben

Entwickle die Divergenzen auf der linken Seite und benutze fur ¨ die Impuls¨ und Energiegleichung die Kontinuitatgleichung ∂ρ + (~u · ∇)ρ ∂t ∂~u + (~u · ∇)~u ∂t ∂p + (~u · ∇)p ∂t

= −ρ ∇ · ~u

(6)

1 = − ∇p + ~k ρ

(7)

= −γp∇·~u

(8)

¨ Da alle Großen Funktionen von Ort (~r ) und Zeit (t) sind, z.B. ρ(~r , t), kann fur ¨ die linke Seite die totale Zeitableitung geschrieben werden. z.B. fur ¨ die Dichte Dρ ∂ρ = + (~u · ∇)ρ = −ρ ∇ ·~u (9) Dt ∂t der Operator D ∂ = + ~u · ∇ (10) Dt ∂t heißt substantielle Zeitableitung (entspricht der totalen Zeitableitung, d/dt) W. Kley

Computational Physics Praktikum

3

Hydrodynamik:

Lagrange-Formulierung

Benutze die substantielle Ableitung Dρ Dt D~u Dt Dp Dt

= −ρ∇ · ~u

(11)

1 = − ∇p + ~k ρ

(12)

= −γp∇·~u

(13)

¨ ¨ ¨ Beschreibt zeitliche Anderung der Großen in einem mit der Stromung mitbewegten System = Lagrange-Formulierung Die Lagrangeformulierung kann z.B. gut bei radialen Stern-Oszillationen verwendet werden. Ist 1D Problem. Hier durch mitbewegte Masseschalen Bei Euler-Formulierung (auf Gitter): orstfest ! W. Kley

Computational Physics Praktikum

4

Numerische Hydrodynamik:

Problemstellung

¨ Betrachte die volle Entwicklung der zeitabhangigen hydrodynamischen Gleichungen. Die nicht-linearen partiellen Differentialgleichungen der ¨ Hydrodynamik werden numerisch gelost: Kontinuum ⇒ Diskretisierung

W. Kley

Computational Physics Praktikum

5

Numerische Hydrodynamik:

¨ Losungsverfahren

Gitter-Methoden (Euler)

festes Gitter ¨ - Stromung durch Gitter  ρ

∂~u + ~u · ∇~u ∂t

bewegtes Gitter/Teilchen ¨ - Stromung bewegt Gitter ρ

 = −∇p

Methoden: • Finite Differenzen keine Erhaltungseigenschaften • Control Volume Erhaltungseigenschaften ¨ • Riemann-Loser Welleneigenschaften ¨ • Problem: Diskontinuitaten W. Kley

Teilchen-Methoden (Lagrange)

Computational Physics Praktikum

d ~u = −∇p dt

Bekannte Methode: Smoothed Particle Hydrodynamics, SPH

’ausgeschmierte Teilchen’ ¨ Gut fur Eigengraviation ¨ freie Rander, 6

Numerische Hydrodynamik:

betrachte: 1D Eulergleichungen

Beschreiben Erhaltung von: Masse, Impuls und Energie ∂ρ ∂ρu + ∂t ∂x ∂ρu ∂ρuu + ∂t ∂x ∂ρ ∂ρu + ∂t ∂x

=

0

=

−

=

∂p ∂x ∂u −p ∂x

ρ: Dichte u: Geschwindigkeit p: Druck : innere spezifische Energie (Energie/Masse) Mit Zustandsgleichung p = (γ − 1)ρ

(14) (15) (16)

(17)

γ: Adiabatenexponent Partielle Dgl. in Raum und Zeit. → Brauche Diskretisierung in Raum und Zeit. W. Kley

Computational Physics Praktikum

7

Numerische Hydrodynamik:

Diskretisierung

ψ

Betrachte Funktion: ψ(x, t) Diskretisierung im Raum ¨ Uberdeckung mit Gitter ∆x =

j−1

ψjn

j

xmax − xmin N

X

j+1

Zellmittelwert von ψ(x, t) am Gitterpunkt xj zum Zeitpunkt t n ψjn = ψ(xj , t n ) ≈

1 ∆x

Z

(j+1/2)∆x

ψ(x, n∆t)dx (j−1/2)∆x

¨ ψjn ist stuckweise konstant. j raumlicher Index, n Zeitschritt. ¨ W. Kley

Computational Physics Praktikum

8

Numerische Hydrodynamik:

Zeitintegration

Betrachte allg. Gleichung ∂ψ = L(ψ(x, t)) (18) ∂t ¨ mit einem (raumlichen) Differentialoperator L. Typische Diskretisierung (1. Ordnung in der Zeit), z.Zt.: t = t n = n∆t ψ(t + ∆t) − ψ(t) ψ n+1 − ψ n ∂ψ ≈ = = L(ψ n ) ∂t ∆t ∆t

(19)

Jetzt am Ort, dem Gitterpunkt xj mit Umformung ψjn+1 = ψjn + ∆tL(ψkn )

(20)

L(ψkn ): diskretisierter Differentialoperator L (hier explizit) ¨ - k in L(ψk ): Satz von raumlichen Indizes: - typisch bei 2. Ordnung: k ∈ {j − 2, j − 1, j, j + 1, j + 2} (brauche Information von links und rechts, 5 Punkt ’Stencil’) W. Kley

Computational Physics Praktikum

9

Numerische Hydrodynamik:

Operator-Splitting

~ ∂A ~ + L2 (A) ~ = L1 (A) ∂t

(21)

~ i = 1, 2 einzelne (Differential-)Operatoren Li (A), ~ = (ρ, u, ). ¨ angewandt auf die Großen A Hier bei 1D idealer Hydrodynamik L1 : Advektion ¨ L2 : Druck, bzw. ext. Krafte ¨ Zur Losung in einzelne Unterschritte unterteilt

~ n+1 A

~1 A ~2 =A

~ n + ∆tL1 (A ~ n) = A ~ 1 + ∆tL2 (A ~ 1) = A

(22)

Li ist Differenzenoperator zu Li . W. Kley

Computational Physics Praktikum

10

Numerische Hydrodynamik:

Advektions-Schritt

∂ρ ∂t ∂(ρu) ∂t ∂(ρ) ∂t

∂ρu ∂x ∂(ρuu) = − ∂x ∂(ρu) = − ∂x = −

In expliziter Erhaltungsform ∂~u ∂~f (~u ) + =0 ∂t ∂x

(23)

Fur ¨ ~u = (u1 , u2 , u3 ) und ~f = (f1 , f2 , f3 ) gilt: ~u = (ρ, ρu, ρ) und ~f = (ρ u, ρuu, ρu). Dieser Schritt ergibt: W. Kley

ρn → ρ1 = ρn+1 ,

Computational Physics Praktikum

un → u1,

n → 1 11

Kraftterme

Numerische Hydrodynamik: Impulsgleichung

1 ∂p ∂u =− ∂t ρ ∂x ujn+1

= uj − ∆t

1 ρ¯jn+1

pj − pj−1 ∆x

(24)
fur ¨

j = 2, N

(25)

Energiegleichung ∂ p ∂u =− ∂t ρ ∂x n+1 j

= j − ∆t

pj ρjn+1

uj+1 − uj ∆x

(26)
fur ¨

j = 1, N

(27)

auf der rechten Seite werden jeweils die momentanen Werte fur ¨ u,  und p eingesetzt, also hier u 1 , p1 , 1 . Dieser Schritt ergibt: W. Kley

Computational Physics Praktikum

u 1 → u n+1 ,

1 → n+1 12

Numerische Hydrodynamik:

Modellgleichung fur ¨ Advektion

¨ Kontinuitatsgleichung lautete ∂ρ ∂ρu + =0 ∂t ∂x

(28)

Hier ist F m = ρu der Massenfluss Geht mit ρ → ψ und u → a = const. uber in die Lineare ¨ Advektionsgleichung ∂ψ ∂ψ +a = 0. ∂t ∂x ¨ Bei konstantem a ist die Losung eine nach rechts laufende Welle mit

(29)

ψ(x, t = 0) = f (x) wird ψ(x, t) = f (x − at)

Hierbei ist f (x) die Anfangsbedingung zur Zeit t = 0, welche durch Advektion mit konstanter Geschwindigkeit a verschoben wird. ¨ Die Numerik sollte diese Eigenschaft bestmoglichst erhalten. W. Kley

Computational Physics Praktikum

13

Numerische Hydrodynamik: FTCS:

Lineare Advektion

Forward Time Centered Space Algorithmus ∂ψ ∂ψ +a =0 ∂t ∂x

(30)

ψj−1

ψj

ψj+1

x j−1

xj

x j+1

Spezifiziere Gitter : und schreibe

n ψjn+1 − ψjn ∂ψ = ∂t j ∆t n n n ψj+1 − ψj−1 ∂ψ = ∂x j 2 ∆x

(31) (32)

es folgt ψjn+1 = ψjn −


a∆t n n ψj+1 − ψj−1 2∆x

(33)

Methode sieht gut motiviert aus: ist aber Instabil fur ¨ alle ∆t ! W. Kley

Computational Physics Praktikum

14

Numerische Hydrodynamik: ∂ψ ∂aψ + =0 ∂t ∂x

Upwind-Methode I (34)

oder

∂ψ ∂ψ +a =0 (35) ∂t ∂x a: Konstante (Geschwindigkeit) > 0 ¨ ψ(x, t) beliebige Transportgroße

ψ

α

a ∆t

¨ Anderung von ψ in Zelle j in out

ψjn+1 ∆x = ψjn ∆x +∆t(Fin −Fout ) (36) Der Fluss Fin ist fur ¨ konstante ψj Fin Fout

n = a ψj−1

=

a ψjn

j−1

(37) (38)

j

j+1

X

Violette Bereiche werden in Nachbarzelle transportiert.

Upwind-Methode ¨ Information kommt von Stromaufwarts W. Kley

Computational Physics Praktikum

15

Upwind-Methode II

Numerische Hydrodynamik: Erweiterung ¨ Zustande

fur ¨

nicht

konstante

(39)

ψI (x) Interpolationspolynom. Hier lineare Interpolation (Gerade) Damit   n 1 (40) Fin = a ψj−1 + (1 − σ)∆ψj−1 | {z } 2 1st Order

{z

2nd Order

mit σ ≡ a∆t/∆x ∆ψj ≈ W. Kley

∂ψ ∆x ∂x xj

Computational Physics Praktikum

a ∆t

ψ I (x)

  a∆t Fin = a ψI xj−1/2 − 2

|

α ψ

in out

j−1

ψI (x) = ψjn +

j

j+1

X

x − xj ∆ψj (41) ∆x

} ∆ψj undividierte Differenz, approximiert Ableitung Zweite Ordnung Upwind ψI (x) wird jeweils in der Mitte der violetten Bereiche ausgewertet. 16

Numerische Hydrodynamik:

Undividierte Differenz

Verschiedene Methoden: a) ∆ψj = 0 Upwind, 1st Order, stuckweise konstant ¨
1 b) ∆ψj = 2 ψj+1 − ψj−1 Fromm, zentrierte Ableitung c) ∆ψj = ψj − ψj−1 Beam-Warming, Upwind Steigung d) ∆ψj = ψj+1 − ψj Lax-Wendroff, Downwind Steigung Meist 2nd Order Upwind (van Leer Schema) ¨ die Monotonizitat) ¨ (erhalt

∆ψj =

 (ψj+1 −ψj )(ψj −ψj−1 )   2 (ψj+1 −ψj−1 )

falls

 

sonst.

0

Geometric Mean

(ψj+1 − ψj )(ψj − ψj−1 ) > 0

(42) Die Ableitungen werden jeweils zum entsprechenden Zeitschrittlevel bzw. Zwischenschritt berechnet. W. Kley

Computational Physics Praktikum

17

Numerische Hydrodynamik:

t t

n+1

n+1

ψj

n+1/2

tn

Lax-Wendroff-Methode

n ψj−1

n+1/2 ~ ψ j−1/2

¨ Schematischer Uberblick uber ¨ das Verfahren ¨ verwendet raumlich und zeitlich zentrierte Differenzen

n+1/2

~ ψ j+1/2

ψjn

n ψj+1

Dadurch zweite Ordnung in Raum und Zeit

Durch zwei Schritte: ¨ Pradiktor-Schritt (z.Zt. t n+1/2 )  σ  1 n n+1/2 n n ψ˜j+1/2 = ψj + ψj+1 − ψj+1 − ψjn 2 2 n+1 Dann Korrektor-Schritt (auf Zt. t )   n+1/2 n+1/2 ψjn+1 = ψjn − σ ψ˜j+1/2 − ψ˜j−1/2 W. Kley

Computational Physics Praktikum

(43)

(44) 18

Numerische Hydrodynamik:

Beispiel: Linare Advektion

ψ(x)

Upwind - (Diffusiv) Linear Advection (N=600,t=40): Upwind, sigma=0.8 1 Analytic Upwind

0.8

0.6 Phi-Achse

Rechteckfunktion: Breite 0.6 im Intervall [−1, 1] Geschwindigkeit a = 1, bis t = 40 periodische Randbedingung σ = a∆t/∆x = 0.8 - (Courantzahl)

0.4

0.2

0 -1

-0.5

0 X-Achse

0.5

1

Lax-Wendroff - (Dispersiv)

Van Leer

Linear Advection (N=600,t=40): Van Leer, sigma=0.8

Linear Advection (N=600,t=40): Lax Wendroff, sigma=0.8

1

1.4

Analytic Van Leer

Analytic Lax Wendroff 1.2

0.8 1 0.8 Phi-Achse

Phi-Achse

0.6

0.4

0.6 0.4 0.2 0

0.2

-0.2 0 -1

W. Kley

-0.5

0 X-Achse

Computational Physics Praktikum

0.5

1

-0.4 -1

-0.5

0 X-Achse

0.5

1

19

Numerische Hydrodynamik:

¨ Stabilitatsanalyse I

¨ Setze fur eine Fourier-Reihe ein (von Neumann 1940/50er) ¨ Losung Betrachte vereinfachend eine Komponente, untersuche deren Wachstum ψjn = V n eiθj

(45)

¨ ¨ mit der Definition von θ uber die Gittergroße ∆x und der Gesamtlange L ¨ θ=

2π∆x L

(46)

Betrachte jetzt Upwind Verfahren mit σ = a∆t/∆x n )=0 ψjn+1 − ψjn + σ(ψjn − ψj−1

(47)

Einsetzen von Gleichung (45) h i V n+1 eiθj = V n eiθj + σV n eiθ(j−1) − eiθj Teilen durch V n und eiθj liefert   V n+1 −iθ − 1 = 1 + σ e Vn W. Kley

Computational Physics Praktikum

(48) 20

Numerische Hydrodynamik:

¨ Stabilitatsanalyse II

¨ man Fur ¨ das Betragsquadrat erhalt n+1 2 h  i h  i V λ(θ) ≡ n = 1 + σ e−iθ − 1 1 + σ eiθ − 1 V     = 1 + σ e−iθ + eiθ − 2 − σ 2 e−iθ + eiθ − 2 = 1 + σ(1 − σ)(2 cos θ − 2)   θ 2 = 1 − 4σ(1 − σ) sin 2

(49)

¨ Ein Verfahren ist stabil, falls der Betrag des Verstarkungsfaktors λ(θ) kleiner als eins ist. Das Upwind-Verfahren ist also stabil fur ¨ 0 < σ < 1, dann ist |λ(θ)| < 1. Umgeschrieben ∆t < fCFL

∆x a

(50)

¨ ¨ mit dem Courant-Faktor fCFL < 1. Ublich ist fCFL = 0.5. (in 1D auch großer) Satz: Courant-Friedrich-Levy Es existiert kein explizites, konsistentes, stabiles finites Differenzenschema, welches bedingungslos stabil ist (d.h. fur ¨ alle ∆t). W. Kley

Computational Physics Praktikum

21

Numerische Hydrodynamik:

Modifizierte Gleichung I

Betrachte Upwind Verfahren mit σ = a∆t/∆x n ψjn+1 − ψjn + σ(ψjn − ψj−1 )=0

(51)

Ersetze Differenzen durch Ableitungen, d.h. Taylor-Entwicklung (bis 2. Ordnung)

∂ψ 1 ∂2ψ 2 ∆t+ ∆t +O(∆t 3 )+σ ∂t 2 ∂t 2



 ∂ψ 1 ∂2ψ 2 ∆x − ∆x +O(∆t∆x 2 ) = 0 ∂x 2 ∂x 2 (52)

Teile durch ∆t, ersetze σ   ∂2ψ ∂ψ ∂ψ 1 ∂ 2 ψ ∆t − a 2 ∆x + O(∆t 2 ) + O(∆x 2 ) = 0 (53) +a + ∂t ∂x 2 ∂t 2 ∂x Benutze Wellengleichung ψtt = a2 ψxx ⇒ modifizierte Gleichung (Index M) ∂ψM ∂ψM 1 ∂ 2 ψM +a = a∆x (1 − σ) ∂t ∂x 2 ∂x 2 D.h. die FDG addiert einen diffusiven Term zur ursprunglichen PDG ¨ W. Kley

Computational Physics Praktikum

(54)

22

Modifizierte Gleichung II

Numerische Hydrodynamik: mit Diffusions-Koeffizienten

1 a∆x (1 − σ) (55) 2 Bem: Nur fur ¨ D > 0 ist dies eine Diffusionsgleichung, woraus σ < 1 ¨ folgt (Hirt-Methode). Fur fur ¨ die Stabilitat ¨ Upwind-Methode ist D > 0 ⇒ Diffusion. D=

Lax-Wendroff liefert  ∂3ψ ∂ψM ∂ψ ∆t 2 a  2 M +a M = σ −1 ∂t ∂x σ ∂x 3

(56)

Die Gleichung hat also die Form ψt + aψx = µψxxx
∆t 2 a 2 σ −1 σ Dies verursacht Dispersion. Hier: Wellen zu langsam (µ < 0) ¨ (vgl. Rechteckfunktion) ⇒ Oszillationen hinter Diskontinuitat µ=

W. Kley

Computational Physics Praktikum

(57) (58)

23

Numerische Hydrodynamik:

¨ Zeitschrittgroße

¨ Aus obiger Analyse folgt, dass die Zeitschrittgroße ∆t limitiert sein muss fur ¨ eine stabile numerische Entwicklung. Bei der linearen Advektion (mit der Geschwindigkeit a) gilt ∆x (59) ∆t < a Im allgemeinen Fall geht die Schallgeschwindigkeit (cs ) ein, und es ergibt sich die Courant-Friedrich-Lewy-Bedingung ∆x ∆t < (60) cs + |~u | Anschaulich heißt dies, dass sich die Information in einem Zeitschritt nicht uber eine Gitterzelle hinweg ausbreiten darf. Man schreibt ¨ ublicherweise ¨ ∆x ∆t = fC (61) cs + |~u | mit dem Courant-Faktor fC . Typisch ist: fC ∼ 0.5. ¨ Nur bei impliziten Verfahren gibt es (theoretisch) keine Beschrankung von ∆t. W. Kley

Computational Physics Praktikum

24

Numerische Hydrodynamik:

¨ Zeitschrittgroße graphisch

¨ Der numerische Abhangigkeitsbereich (gestrichelter Bereich) muss ¨ großer als der physikalische (grauer Bereich) sein (∆x/∆t > a). Die gesamte Information des physikalischen ’Schallkegels’ muss berucksichtigt werden. ¨ W. Kley

Computational Physics Praktikum

25

Numerische Hydrodynamik:

Multi-dimensional

Gitterdefinition (in 2D, staggered): Skalare in Zellzentren (hier: ρ, , p, v3 , ψ) Vektoren an Interfaces (hier: v1 , v2 )

¨ Flusse uber Zellrander: ¨ ¨ Oben: Masseflusse ¨ Unten: X-Impuls (Gitter verschoben)

aus: ZEUS-2D: A radiation magnetohydrodynamics code for astrophysical flows in two space dimensions. I in The Astrophysical Journal Suppl., von Jim Stone und Mike Norman, 1992. Benutze Operator-Splitting und Directional Splitting: Die x und y Richtung werden nacheinander abgehandelt. Erst x-scans, dann y -scans. W. Kley

Computational Physics Praktikum

26

Numerische Hydrodynamik:

Zusammenfassung

Numerische Methoden sollten Erhaltungseigenschaften wiedergeben. Gleichungen in Erhaltungsform schreiben Numerische Methoden sollten Welleneigenschaften wiedergeben. ¨ Shock-Capturing Methoden, Riemann-Loser ¨ Numerische Methoden mussen Diskontinuitaten unter Kontrolle halten. ¨ ¨ Brauche dazu Diffusion (⇒ Stabilitat) ¨ oder intrinsisch (durch Verfahren) entweder explizit (kunstliche Viskositat) ¨ Numerische Methoden sollten genau sein mind. 2. Ordnung in Zeit und Raum Frei verfugbare Codes: ¨ ZEUS: http://www.astro.princeton.edu/~jstone/zeus.html klassischer Upwind-Code, zweiter Ordnung, staggered grid, RMHD ATHENA: https://trac.princeton.edu/Athena/ ¨ Nachfolger von ZEUS: Riemann Loser, zentriertes Gitter, MHD NIRVANA: http://nirvana-code.aip.de/ 3D, AMR, Finite Volume Code, MHD PLUTO: http://plutocode.ph.unito.it/ ¨ 3D, relativistisch, Riemann-Loser/Finite Volume, MHD W. Kley Computational Physics Praktikum GADGET: http://www.mpa-garching.mpg.de/galform/gadget/

27

Hydrodynamik:

Wellenstruktur

Betrachte eindimensionale Gleichungen (Bewegung in x-Richtung): Aus Eulergleichungen: Mit p = (γ − 1)ρ und Trennen der Ableitungen ∂ρ ∂ρu ∂t + ∂x ∂ρu ∂ρuu ∂t + ∂x ∂ρ ∂ρu ∂t + ∂x

 = 0  ∂p = − ∂x ∂u  = −p ∂x

=⇒

∂ρ ∂ρ ∂u ∂t + u ∂x + ρ ∂x ∂u ∂u 1 ∂p ∂t + u ∂x + ρ ∂x ∂p ∂p ∂u ∂t + u ∂x + γp ∂x

= 0 = 0 = 0

Als Vektorgleichung

mit



 ρ W= u  p

∂W ∂W +A =0 ∂t ∂x  und

u A= 0 0

(62) ρ u γp

 0 1/ρ  u

(63)

Gleichungen sind nichtlinear und gekoppelt. Versuche Entkopplung: ⇒ Diagonalisierung von A

W. Kley

Computational Physics Praktikum

28

Hydrodynamik:

Diagonalisierung

Eigenwerte (EW) = (u − λ) u − λ γp h i = (u − λ) (u − λ)2 − γp/ρ = 0

u−λ 0 det(A) = 0

ρ u−λ γp

0 1/ρ u−λ

1/ρ u−λ (64)

Es folgt: λ0 = u λ± = u ± cs

(65)

mit der Schallgeschwindigkeit

γp ρ Die Eigenwerte geben charakteristische Geschwindigkeiten, mit denen sich die Information ausbreitet. Setzt sich aus Fluid- (u) und Schallgeschwindigkeit (cs ) zusammen 3 reelle Eigenwerte ⇒ A diagonalisierbar cs2 =

Q−1 AQ = Λ

(66)

(67)

Q setzt sich aus Eigenvektoren zu EW (λi , i = 0, +, −) zusammen, Λ ist Diagonalmatrix. W. Kley

Computational Physics Praktikum

29

Hydrodynamik:

Charakteristische Variablen

Fur ¨ Q folgt



1 Q= 0 0

1 ρ 2 cs 1 2 1 2 ρcs

− 12 cρs

1 2 1 − 2 ρcs

Hatten



 1 0 − c12 s   1 = 0 1  ρcs 0 1 − ρc1 s 

 und

Q−1

∂W ∂W +A =0 ∂t ∂x

(68)

und Q−1 AQ = Λ Definiere: dv ≡ Q−1 dW Multipliziere Gl. (68) mit Q

also

dW = Qdv

(69)

−1

∂v ∂v +Λ =0 ∂t ∂x v = (v0 , v+ , v− ) sind charakteristische Variable: vi = const. auf Kurven

(70)

dx = λi dt W. Kley

Computational Physics Praktikum

30

Hydrodynamik:

Variable v0

aus den Definitionen dv0 = dρ −

1 dp cs2

∂v0 ∂v0 + λ0 =0 ∂t ∂x

mit

(71) λ0 = u

was ist dv0 ? ¨ Aus Thermodynamik (1. Hauptsatz fur ¨ spezifische Großen)     p 1 1 Tds = d + p d = d − 2 d ρ ρ ρ mit p = (γ − 1)ρ,  = cv T , γ = cp /cv folgt   cp cp dp ds = − dρ − 2 = − dv0 ρ ρ cs ∂s ∂s +u =0 ∂t ∂x d.h. s ist const. entlang Stromlinien, also =⇒

ds =0 dt W. Kley

Computational Physics Praktikum

(72)

(73)

(74) (75)

(76) 31

Hydrodynamik:

Riemann-Invarianten

Fur ¨ weitere Variablen

∂v± ∂v± + (u ± cs ) =0 ∂t ∂x

mit dv± = du ±

(77)

1 dp ρcs

(78)

dp . ρcs

(79)

folgt Z v± = u ± Sei Entropie uberall konstant (d.h. p = K ργ ) ¨ =⇒

v± = u ±

2cs γ−1

(80)

Riemann-Invarianten: v± = const. auf Kurven dx = u ± cs dt

W. Kley

Computational Physics Praktikum

32

Hydrodynamik:

Aufsteilen von Schallwellen

¨ Linearisierung der Euler-Gleichungen fuhrt auf Wellengleichung fur ¨ ¨ Storung: ∂ 2 ρ1 ∂ 2 ρ1 = cs2 2 ∂t ∂x 2 aber: cs nicht konstant ⇒ Aufsteilen

Bspl. fur Stoßwelle ¨ (rucklaufende) ¨ (81)

¨ ⇒ Diskontinuitaten ¨ ≡ Sprung: UberUnterschall W. Kley

Computational Physics Praktikum

33

Beispiele:

Stoßrohr - Shocktube

¨ in einem Rohr am Ort x0 (eindimensional) Anfangsdiskontiuitat Bereich 2

Riemann-Problem

Sprung in Druck (p) und Dichte (ρ) es entwickelt sich:

Bereich 1 Diskontinuität

P2

ρ

2

ρ

1

P1

x0

W. Kley

Computational Physics Praktikum

- eine Stoßwelle nach rechts (X4 ) ¨ (mit Uberschall ush > cs ) ¨ - eine Kontaktdiskontinuitat Dichtesprung (entlang X3 ) - eine Verdunnungswelle ¨ (zwischen X1 und X2 )

34

Sod-Shocktube

Beispiele:

Standard Testproblem fur ¨ numerische Hydrodynamik, x ∈ [0, 1] mit X0 = 0.5, γ = 1.4 ρ1 = 1.0, p1 = 1.0, 1 = 2.5, T1 = 1 und ρ2 = 0.1, p2 = 0.125, 2 = 2.0, T2 = 0.8 Shock-Tube: Sod; Mono: Geometric Mean; Nx=100, Nt=228, dt=0.001 1

1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0

0.9 0.8 0.7 Density

Velocity

¨ Hier Losung mit van Leer Verfahren (z.Zt. t = 0.228 nach 228 Zeitschritten:)

0.3 0.2 0.1 0.2

0.4

0.6

0.8

1

1.15

1

1.1

0.9

1.05

0.8

1

0.7

Pressure

Temperature

0.5 0.4

0

0.95 0.9 0.85 0.8

0

0.2

0.4

0.6

0.8

1

0

0.2

0.4 0.6 X-Axis

0.8

1

0.6 0.5

Rot: Exakt Grun: ¨ Numerik ¨ Losung ist: ¨ selbstahnlich gehen durch Streckung auseinander hervor.

0.4 0.3

0.75

0.2

0.7

0.1 0

W. Kley

0.6

0.2

0.4 0.6 X-Axis

0.8

Computational Physics Praktikum

1

35

Beispiele:

Sedov-Explosion

¨ Als Beispiel fur ¨ Bombenexplosionen (Sedov & Taylor, 1950er), analytische Losung (Sedov) ¨ ¨ Grundproblem fur z.B. Abschatzung der Remnant-Große ¨ Supernovae-Ausbruche, ¨ Standard Testproblem mehr-dimensionale Hydrodynamik fur ¨ x, y ∈ [0, 1] × [0, 1] Energie-Input am Ursprung, E = 1, in ρ = 1, γ = 1.4, 200 × 200 Gitterpunkte ¨ ¨ fur Hier Losung mit van Leer Verfahren (lose ¨ totale Energie-Variable). Dargestellt: Dichte

W. Kley

Computational Physics Praktikum

36

Beispiele:

Wassertropfen: SPH

Wasserkugel (R=30cm), Wanne (1x1 m, 60cm hoch) Incl. ¨ Oberflachenspannung, Zeit in Sekunden (TU-Munchen, 2002) ¨

(Website) W. Kley

Computational Physics Praktikum

37

Beispiele:

Sternentstehung: SPH Molekulare Wolke Masse: 50 M Durchmesser: 1.2 LJ = 76,000 AE Temperatur: 10 K

(M. Bate, 2002)

W. Kley

Computational Physics Praktikum

38

Beispiele:

¨ Prinzip Kelvin-Helmholtz Instabilitat:

¨ Sprung in Tangential-Geschwindigkeit. Die KHI ist eine Scherinstabilitat. W. Kley

Computational Physics Praktikum

39

Beispiele:

¨ In Atmosphare ¨ KH-Instabilitat:

(Boulder (NCAR), USA) W. Kley

Computational Physics Praktikum

40

Beispiele:

¨ Simulation KH-Instabilitat:

Direkter Vergleich: Bewegtes < − > Festes Gitter Links: Gitter bewegt (Voronoi-Tesselation) Rechts: festes quadratisches Gitter (Euler)

mit Gitterbewegung

(Kevin Schaal, Tubingen) ¨

Youtube channel

W. Kley

Computational Physics Praktikum

41

Beispiele:

¨ Rayleigh-Taylor Instabilitat

PPM Code, 128 Knoten, ASCI Blue-Pacific ID System at LLNL 5123 Gitterpunkte (LLNL, 1999)

(Web-Link) W. Kley

Computational Physics Praktikum

42

Beispiele:

Diesel Injektion

Finite Volumen Methode (FOAM) ¨ Geschwindigkeit, Temperatur, Teilchen (+Isoflache) (Nabla Ltd, 2004)

W. Kley

Computational Physics Praktikum

43

Beispiele:

Kataklysmische Variable: Gitter

.6

.4

y

.2

0

−.2

−.4

−.6

−1

−.8

−.6

−.4

−.2

0

.2

.4

.6

x W. Kley

Computational Physics Praktikum

44

Beispiele:

Kataklysmische Variable: Scheibenbildung

RH2D Code, Mono-Verfahren ¨ 5122 Gitterpunkte, Massenverhaltnis: q = m2 /m1 = 0.15

W. Kley

Computational Physics Praktikum

45