(Kapitel 1) Im Allgemeinen:

(Kapitel 7) Spezifisch auf die Vorlesungskapitel:

• Modellierung

• Diskretisierung

• Lösung

• Der Modellierungsfehler bezeichnet die Differenz zwischen dem exakten Variablenverlauf beim physikalischen Problem und der exakten Lösung der Differenzialgleichung, d.h. |Φphy(x, t)−Φmod (x, t)|.

• Der Diskretisierungsfehler bezeichnet die Differenz zwischen der exakten Lösung der Differenzialgleichung (inkl. Rand– und Anfangsbedingungen) und der exakten Lösung der diskretisierten Gleichung, d.h. |Φmod(x, t) − Φn |.

• Der Lösungsfehler bezeichnet die Differenz zwischen der exakten Lösung der diskretisierten Gleichung und der numerisch tatsächlich berechneten Lösung, d.h. |Φn − �Φn|

• Der gesamte Numerikfehler als Summe aus Diskretisierungs– und Lösungsfehler besteht demnach aus der Differenz zwischen der exakten Lösung der Differenzialgleichung und der numerisch berechneten Lösung, d.h. |Φmod (x, t) − �Φn |

Konsistenz, Stabilität und Konvergenz müssen erfüllt sein.

Darunter versteht man den Zusammenhang des Differenzialproblems mit seiner diskreten Formulierung. Konkreter spricht man von einem konsistenten Schema, wenn die diskretisierte Gleichung gegen die zugrundeliegende Differenzialgleichung strebt bei kleiner werdenden Schrittweiten, bzw. der Fehler strebt gegen 0.

Stabilität beschreibt das Verhalten bei externen Störungen. Klingen die Störungen ab, spricht man von einem stabilen Schema, schaukeln sie sich auf, von einem instabilen.

Zur Kontrolle:

• Diskrete Störtheorie: Einbringen einer kleinen Störung (psi = psi + epsilon) an einer Diskreten Stelle. Berechnen zu zwei Zeitschritten: Unter welcher Bedingung hat sich die Störung verkleinert. (Seite 114)

• Von Neumannsche Stabilitätsanalyse: Einbringen einer Störungsfunktion(x,t), dann analytische Analyse gegen unendlich möglich. Genauer Kapitel 7.2.2.2

Konsistenz + Stabilität (bei lin. Anfangswertproblemen)

—> Das Problem konvergiert zu einer Lösung

Die Verfahren beruhen auf der Lösung der Erhaltungsgleichungen, ein Diskretisierungsverfahren sollte diese daher zumindest näherungsweise erfüllen. Ist die Erhaltung sichergestellt, kann es nur eine ungleiche Verteilung über das Lösungsgebiet geben.

Diskretisierungsverfahren, die die Erhaltungseigenschaften definitionsgemäß erfüllen, heißen konservativ (z.B. FVM).

Manche nicht-konservative Verfahren können mit Aufwand umgewandelt werden, um die Erhaltung zu gewährleisten, bei anderen geht das nicht. Hier können künstliche Quellen und Senken entstehen, die die lokalen und globalen Bilanzen stören können.

Trotzdem können auch Verfahren, die die Erhaltungseigenschaften nicht vollständig gewährleisten, ausreichend genaue Lösungen liefern

• Diffusion: elliptischen Charakter, d.h. pflanzen sich in alle Richtungen fort. Bei der Diskretisierung daher ausschließlich zentrale Fomulierungen.

• Konvektion: hyperbolischer Charakter, d.h. der Transport erfolgt in Strömungsrichtung. Diesen Transport versucht man über Upwindverfahren darzustellen.

  
  

Numerische Verfahren können physikalisch unsinnige Lösungen präsentieren. Dies muss überprüft werden, ob ein Verfahren das macht.

Beispiele:

• Dichte muss> 0

• Konzentration zwischen 0 und 100 %

• Ohne Quellen/Senken kann sich die Verteilung einer Größe nur innerhalb der Randwerte bewegen (z.B. Temperaturverlauf zwischen zwei Platten)

Bedingung nicht immer erfüllt, teilweise erkennbar durch Oszillationen, vergleiche Seite 61.

Sie können numerische Diffusion und Dispersion erzeugen:

• Dispersive Fehler: bei ungeraden Ableitungen, erzeugt Phasenverschiebung, keine Amplitudenänderung

• Diffusive Fehler: bei 2. Ableitung; stark dämpfend/glättend; numerische Diffusion

• Diffusive Fehler: höhere gerade Ableitungen; Hochfrequenzdämpfend

Zur Lösung der Gleichungen gibt es 2 grundsätzliche Ansätze: sequentielle Löser und gekoppelte Löser. Die Unterscheidung gibt an, ob die Erhaltungsgleichungen sequentiell oder gekoppelt gelöst werden. Gekoppelt ergeben sich komplexere Gleichungssysteme mit Matrizen als Diagonaleinträge, jedoch ist die Konvergenz besser als bei sequentiellen Lösern, bei denen iterativ die Gleichungen gelöst werden.

Zur Lösung des Gleichungssystems gibt es dann folgende Lösungsmöglichkeiten:

Direkte Löser, lösen das Gleichungssystem analytisch mit Berechnungsvorschrift:

• Cramersche Regel

• Gaußverfahren

• Thomas Algorithmus

Iterative Löser, die das Gleichungssystem nummerisch iterativ annähern:

• Jakobi Verfahren

• Punkt- und Linien-Gauß-Seidel-Verfahren

• Linieniteration mit alternierenden Richtungen

Speziellere Löser:

• Unvollständige LU-Zerrlegung: versucht die LU-Zerlegung mit Jakobi zu verknüpfen

• Konjugierte Gradientenverfahren: Überführung in ein Minimierungsproblem, Lösung mit z.B. Gradientenverfahren

• Mehrgitterverfahren: Erste Berechnung auf einem gröberen Gitter, kein eigenes Lösungsverfahren, aber deutliche Beschleunigung bei der Berechnung

Matrix A wird zerlegt in A = M - N.

M ist einfach zu invertieren, N nicht. Zur Berechnung von einem Schritt wird M*x^v = N*x^(v-1) + b gelöst. Es wird also ein Gleichungssystem gelöst, bei dem die nicht einfach zu berechnenden Terme aus dem letzten Schritt genommen werden, um den nächsten Iterationsschritt zu berechnen.

M*(x^v-x^(v-1)) = -Residuum(x^(v-1)) wird genommen, um ein Abbruchkriterium zu definieren.

Beispiel: Possionsgleichung, daher folgt folgendes Gleichungssystem:

Die Pentadiagonalmatrix wird zerlegt in die Diagonale, die untere und obere Matrix.

M ist die Diagonale der Bandmatrix, N der Rest.

Da die Matrix nur Werte auf der Diagonalen enthält, ist die Lösung damit trivial:

Dies ist das einfachste iterative Verfahren, die Konvergenz kann deutlich verbessert werden, wenn für die Werte, für die bereits der neue Wert bereitsteht, auch der neue Wert genommen wird (Punkt-Gauß-Seidel). Alternativ kann auch eine ganze Reihe an Punkten gleichzeitig gelöst werden, dies ist aufwendiger, da nicht jeder Punkt explizit berechnet werden kann, aber nicht zu aufwendig, da die Matrix nicht zu groß ist und einfache Löser verwendet werden können (Linien-Gauß-Seidel). Bei wechselnder Linienrichtung nennt sich das ADI Verfahren (Alternating Direction Implicit).

Stabilität/Konvergenz kann ähnlich zu Zeitintegration analysiert werden. Wie die Zeitintegration schreitet der Löser einen Schritt nach dem anderen fort, daher kann die gleiche Analyse angewand werden.

Zum Berechnungsaufwand:

Inkompressible Strömungen sind eine Sonderform der kompressiblen Strömungen. Wenn man die inkompressible Strömung als solche berechnen will, fällt allerdings die Machzahl stark ab, im tatsächlichen Inkompressibilitätsfall sogar auf null. Dies fürt zu schlechter Konvergenz.

Wenn man die Stabilität analysiert, zeigt sich, dass dieser Ansatz nur für sehr kleine Zeitschritte stabil ist. Die Zeit ist abhängig von den räumlichen Schritten und den Eigenwerten des Systems. Bei Ma gegen null, fällt auch der zulässige Zeitschritt auf null.

tl;dr:

Wir eliminieren mithilfe der Stromfunktion und Rotation den Druck vollständig aus den N-S-Gleichungen und berechnen dann über die Possionsgleichungen aus der Stromfunktion und die Wirbeltransportgleichungen das Geschwindigkeitsfeld.

1. Ausgangspunkt: Navier–Stokes für inkompressible Fluide

Kontinuitätsgleichung:div(v) = 0

Impulsgleichung:∂v/∂t + (v · ∇) v = -(1/ρ) ∇p + ν ∇²v

Das Problem: Es gibt keine direkte Gleichung für den Druck, weil in der Kontinuität kein p vorkommt.

2. Druckeliminierung durch Rotation

Wir definieren die Wirbelstärke:ω = ∇ × v

und wenden die Rotation auf die Impulsgleichung an. Wichtig:∇ × ∇p = 0

Die Rotation der einzelnen Terme:

Zeitliche Ableitung:∇ × (∂v/∂t) = ∂ω/∂t

Konvektion (Vektoridentität):∇ × [(v · ∇) v] = (ω · ∇) v - (v · ∇) ω

Viskoser Term:∇ × (ν ∇²v) = ν ∇²ω

Damit ergibt sich die Wirbeltransportgleichung:∂ω/∂t + (v · ∇) ω = (ω · ∇) v + ν ∇²ω

3. Rekonstruktion der Geschwindigkeit

Da der Druck eliminiert wurde, muss die Geschwindigkeit aus einem Potential rekonstruiert werden.

4. Fall 2D: Stromfunktion

Definition der Stromfunktion:u = ∂Ψ/∂y v = -∂Ψ/∂x

Damit gilt automatisch:∂u/∂x + ∂v/∂y = 0

(also Inkompressibilität)

Wirbelstärke in 2D:ω_z = ∂v/∂x - ∂u/∂y

Einsetzen ergibt die Poisson-Gleichung:∇²Ψ = -ω_z

5. Fall 3D: Vektorpotential

Definition des Geschwindigkeitsfeldes über ein Vektorpotential A:v = ∇ × A

Rotation ergibt:ω = ∇ × v = ∇ × (∇ × A)

Mit der Eichbedingung:div(A) = 0

folgt direkt:∇²A = -ω

Das bedeutet: drei Poisson-Gleichungen, eine für jede Komponente von A.

6. Gesamtalgorithmus

Das nummerische Vorgehen lässt sich so zusammenfassen:

1. Wirbelstärke berechnen:ω = ∇ × v

2. Wirbeltransportgleichung lösen:∂ω/∂t + (v · ∇) ω = (ω · ∇) v + ν ∇²ω

3. Potential berechnen:

   in 2D:∇²Ψ = -ω_z

   in 3D:∇²A = -ω

4. Geschwindigkeitsfeld rekonstruieren:

   2D:u = ∂Ψ/∂y v = -∂Ψ/∂x

   3D:v = ∇ × A

Damit ist das Geschwindigkeitsfeld wieder divergenzfrei.

7. Nachteile des Verfahrens

• Man löst mehr Gleichungen als im druckbasierten Ansatz.

• Die Randbedingungen für ω und A sind mathematisch anspruchsvoll.

• Die Lösung der Poisson-Gleichungen ist der rechenintensivste Teil.

Darum werden in der Praxis meist druckbasierte Methoden (SIMPLE, Fractional-Step usw.) bevorzugt.

tl;dr:

Bei der künstlichen Kompressibilität ergänzt man die Kontinuitätsgleichung um eine künstliche Druckzeitableitung, sodass ein hyperbolisches Gleichungssystem entsteht (ein Gleichungssystem ähnlich den kompressiblen Strömungen).

Man integriert dieses in einer Pseudozeit, bis ∂p/∂τ = 0 gilt und damit div(v) = 0 erreicht ist – wodurch die Inkompressibilität erfüllt wird, ohne eine Poisson-Gleichung lösen zu müssen.

Ausgangslage: Problem bei inkompressiblen Strömungen

Für inkompressible Strömungen lautet die Kontinuitätsgleichung:div(v) = 0

Der Druck besitzt keine eigene Zeitableitung, sodass es keine direkte Druckgleichung gibt. Daraus entsteht das Hauptproblem aller druckbasierten Verfahren.

Grundidee der künstlichen Kompressibilität

Chorin führt in die Kontinuitätsgleichung eine künstliche Zeitableitung des Drucks ein:

Original:div(v) = 0

Ersetzt durch:∂p/∂t + β² div(v) = 0

Hier ist β ein Parameter (künstliche Schallgeschwindigkeit), der das „Maß der Kompressibilität“ bestimmt.

Damit erhält das Gleichungssystem hyperbolischen Charakter – ähnlich einer kompressiblen Strömung.

Modifiziertes Gleichungssystem

Kontinuität mit künstlicher Kompressibilität:∂p/∂t + β² div(v) = 0

Impulsgleichung bleibt unverändert:∂v/∂t + (v · ∇) v = -(1/ρ) ∇p + ν ∇²v

Diese beiden Gleichungen werden nun simultan zeitlich integriert, als ob es sich um eine kompressible Strömung handelte.

Physikalische Bedeutung

In einem völlig inkompressiblen Fluid breiten sich Druckstörungen „instantan“ aus (theoretisch unendlich schnell).

Durch künstliche Kompressibilität gilt stattdessen:c = sqrt(u² + β²)

Dies ist eine „künstliche Schallgeschwindigkeit“, welche endlich ist. Somit werden Druckinformationen zeitverzögert, aber kontrolliert übertragen.

Wenn die Lösung stationär wird, gilt:∂p/∂t → 0 ⇒ div(v) = 0

Damit ist die Inkompressibilität wieder erfüllt – ohne eine Poisson-Gleichung lösen zu müssen.

Erweiterung für instationäre Probleme (Dual-Time-Stepping)

Man fügt eine Pseudozeit τ hinzu, um innerhalb jedes realen Zeitschritts das Gleichungssystem zu iterieren.

Modifizierte Gleichungen:

Kontinuität:(1/β²) ∂p/∂τ + div(v) = 0

Impuls:∂v/∂τ + [ ∂v/∂t + (v · ∇) v ] = -(1/ρ) ∇p + ν ∇²v

Hierbei:

• ∂/∂t ist die physikalische Zeit

• ∂/∂τ ist die künstliche Pseudozeit

• Für τ → ∞ konvergiert die Lösung wieder zum inkompressiblen NS-System

Vorteil: Man kann kompressible CFD-Methoden (Riemannlöser, Finite-Volumes, Faktorisierung) direkt verwenden.

Algorithmische Zusammenfassung

1. Druck und Geschwindigkeit zu Pseudozeit τ=0 initialisieren

2. Innerhalb eines Zeitschritts die modifizierten Gleichungen in τ integrieren

3. Konvergenzbedingung:∂p/∂τ → 0 ⇒ div(v) → 0

4. Voranschreiten in physikalischer Zeit t

5. Nächster Zeitschritt

Es entfällt die Poisson-Gleichung, dafür gibt es eine zusätzliche Pseudozeit-Iteration.

tl;dr:

Druckiterationsverfahren berechnen zuerst Geschwindigkeitsfeld mit einem weggelassenen Druckfeld/angenommenen Druckfeld, das die Kontinuitätsgleichung nicht erfüllt. Dann berechnen sie die Druckkorrektur, sodass das Geschwindigkeitsfeld die Kontinuitätsgleichung erfüllen würde. Mit dem neuen Druckfeld berechnet man das richtige Geschwindigkeitsfeld.

ChatGPT:

Grundidee der Druckiterations-/Druckkorrekturverfahren

Die Verfahren basieren auf folgendem Prinzip:

1. Geschwindigkeiten zunächst ohne korrekten Druck berechnen

   Man berechnet ein vorläufiges Geschwindigkeitsfeld, das die Kontinuität nicht erfüllt:

   v* = vⁿ + Δt [ -(vⁿ · ∇) vⁿ + ν ∇²vⁿ -1/ρ*∇p]

   v* ist im Allgemeinen nicht divergenzfrei.

2. Druck bestimmen, der das v-Feld korrigiert

   Man fordert:div(vⁿ⁺¹) = 0

   und nutzt die Korrekturgleichung:vⁿ⁺¹ = v* - (Δt/ρ) ∇p'

   (p' ist die Druckkorrektur)

3. Divergenzbildung liefert Poisson-Gleichung für p' (Druckkorrekur)

   div( vⁿ⁺¹ ) = div( v* - (Δt/ρ) ∇p' ) = 0

   → umgestellt:∇² p' = (ρ/Δt) div(v*)

4. Geschwindigkeit und Druck korrigieren

   vⁿ⁺¹ = v* - (Δt/ρ) ∇p' pⁿ⁺¹ = p* + p'

   Damit ist das korrigierte vⁿ⁺¹ divergenzfrei.

tl;dr:

Fractional-Step-Verfahren sind die Grundidee, umgesetzt mit expliziten Zeitintegrationsverfahren. Der Druck wird je nach Verfahren weggelassen oder approximiert im ersten Schritt. Wird er approximiert, wird in der Possionsgleichung nur ein Korrekturterm berechnet („Druckkorrekturverfahren“).

ChatGPT:

Schritt 1: Prädiktor (Impuls ohne neuen Druck lösen)

Die Impulsgleichung wird zunächst ohne korrekten Druckgradienten integriert:( v* - vⁿ ) / Δt = - (vⁿ · ∇) vⁿ + ν ∇² vⁿ

Hier wird der Druckgradient ausgelassen oder mit einem geschätzten Wert p* ersetzt.

v* erfüllt nicht:div(v*) = 0

Schritt 2: Druckkorrektur durch Poisson-Gleichung

Um vⁿ⁺¹ diverganzfrei zu machen, setzen wir:vⁿ⁺¹ = v* - (Δt/ρ) ∇pⁿ⁺¹

Forderung der Inkompressibilität:div(vⁿ⁺¹) = 0

Einsetzen liefert:div( v* - (Δt/ρ) ∇pⁿ⁺¹ ) = 0

Umstellen ergibt die Poisson-Gleichung für den Druck:∇² pⁿ⁺¹ = (ρ/Δt) div(v*)

Schritt 3: Geschwindigkeitskorrektur

vⁿ⁺¹ = v* - (Δt/ρ) ∇pⁿ⁺¹

Damit sind die Navier–Stokes im neuen Zeitschritt erfüllt und vⁿ⁺¹ ist diverganzfrei.

tl;dr:

SIMPLE löst die Impulsgleichung mit geschätztem Druck, berechnet dann eine Druckkorrektur über eine Poisson-Gleichung und korrigiert anschließend Druck und Geschwindigkeit.

Es ist die Umsetzung der Grundidee mit impliziten Zeitverfahren. Hierbei muss ein Zeitschritt iterativ gelöst werden. Ein Zeitschritt wird iterativ gelöst und dann wird der nächste Schritt berechnet.

Varianten wie SIMPLEC, SIMPLER und PISO verbessern Stabilität und Konvergenz durch unterschiedliche Formen der Druckkorrektur.

Zur besseren Konvergenz wird Unterrelaxation genutzt.

ChatGPT:

Schritt 1: Impulsgleichungen mit geschätztem Druck lösen

A_P u_P* + Σ A_N u_N* = b - (∂p*/∂x)_P

u* erfüllt NICHT div(u*) = 0.

Schritt 2: Druckkorrekturgleichung herleiten

Man definiert die Geschwindigkeit- und Druckkorrekturen:u' = uⁿ⁺¹ - u* p' = pⁿ⁺¹ - p*

Die Geschwindigkeitskorrektur ergibt sich aus:u' = -(1 / A_P) ∇p'

Unter Verwendung der Kontinuität:div(uⁿ⁺¹) = 0 = div(u* + u')

erhält man die Poisson-Gleichung für die Druckkorrektur:∇ · ( (1 / A_P) ∇p' ) = div(u*)

Schritt 3: Druck- und Geschwindigkeitskorrektur

pⁿ⁺¹ = p* + α_p p' uⁿ⁺¹ = u* + u'

(α_p ist der Unterrelaxationsfaktor)

Dieser Prozess wird iteriert, bis div(u) ≈ 0.

Variationen

SIMPLEC

• Verbessert die Approximation der Geschwindigkeitskorrektur.

• Schnellere Konvergenz.

• Kaum Bedarf für Druck-Unterrelaxation.

SIMPLER

• Zusätzliche Druckvorhersagegleichung.

• Stabiler, aber rechenintensiver.

PISO

• Mehrere „Pressure-Implicit Steps“ pro Iteration.

• Besonders geeignet für instationäre Strömungen.

• Deutlich schneller pro Zeitschritt.

Versetzte Gitter umgehen das Interpolieren, um die Oberflächenintegrale zu lösen, indem sie die Variablen direkt dort speichern, wo sie benötigt werden. Die Kontrollvolumina sind damit leicht verschoben. Die einzelnen Gleichungen werden auf unterschiedlichen Kontrollvolumina diskretisiert.

Dadurch entsteht auch eine bessere Kopplung zwischen den Größen. Das passiert durch die fehlende Interpolation, die Oszillationen verursachen kann.

Vor und Nachteile:

Alternative bei nicht versetzten Gittern

• Impulsinterpolation (Rhie-Chow-Methode):

  • Berechnet den Massenfluss durch Zellwände nicht nur durch lineare Interpolation der Geschwindigkeiten, sondern berücksichtigt zusätzlich den Druckgradienten.

  • Effekt: Druckoszillationen erzeugen einen künstlichen Massenfluss → werden im nächsten Iterationsschritt korrigiert.

  • Emuliert die Kopplung wie beim versetzten Gitter.

Strukturierte Gitter: orthogonale Gitter oder Nicht-orthogonale/Konturangepasste Gitter

unstrukturierte Gitter

Bockstrukturierte Gitter

• global unstrukturiert, im Block strukturiert,

• kombiniert Vor- und Nachteile der beiden Ansätze

• Grenzen: matching, non-matching, überlappend

• Blöcke natürliche Struktur zum parallelisieren

• Orthogonalität

• Glattheit

  • Expansionsrate (zwischen 0,9 und 1,1 optimal; nicht weiter als 0,7 bis 1,3)

  • Seitenverhältnis (zwischen 0,1 und 10)

• Algebraisch: Punkte werden auf dem Rand vorgegeben und dann verbunden und unterteilt.

• krummlinig über Koordinatentransformation

• Elliptische Methode: Körper ist als Stromfunktion definiert. Durch Orthogonalität der Strom- und Potenziallinien ergibt sich ein orthogonales, aber krummliniges Gitter