Smoothed Particle Hydrodynamics

Smoothed Particle Hydrodynamics (SPH, im Deutschen: geglättete Teilchen-Hydrodynamik) ist eine numerische Methode, um die Hydrodynamischen Gleichungen zu lösen. Sie wird unter anderem in der Astrophysik, der Ballistik und bei Tsunami-Berechnungen eingesetzt. SPH ist eine Lagrange-Methode, d. h. die benutzten Koordinaten schwimmen mit dem Fluid. SPH ist eine besonders einfach zu implementierende, robuste Methode, da zur Berechnung von Ableitungen kein Gitter benötigt wird.

Methode

Allgemeines

In Smoothed Particle Hydrodynamics wird die zu simulierende Flüssigkeit in Elemente aufgeteilt. Dabei werden, ähnlich den Monte-Carlo-Methoden, die Elemente zufällig über die Flüssigkeit verteilt. Dies minimiert den zu erwartenden Fehler. Der mittlere Abstand dieser Elemente wird durch die Smoothing Length (Glättlänge) h repräsentiert. Sie ist der wichtigste Parameter der Methode. Zwischen den Teilchen wird das Fluid durch den Kernel geglättet, daher der Name. Jede Größe (z. B. die Dichte ρ ) wird durch Summation über alle Teilchen berechnet. Jedes einzelne Teilchen erhält einen Anteil, in Form eines Skalars an dieser Größe. Dadurch werden aus den partiellen Differentialgleichungen der Hydrodynamik gewöhnliche Differentialgleichungen, was die Berechnungen sehr vereinfacht. SPH ist eine sehr empirische Methode. Das bedeutet, dass viele Dinge gemacht werden, weil sie funktionieren, nicht, weil es eine strenge mathematische Herleitung gibt.

Herleitung

Die formale Herleitung läuft entweder über eine Lagrange-Funktion oder über eine Integralinterpolation. Bei der Integralinterpolation für eine Größe A₁ geht man von einer Identität aus, wobei $\delta (\vec{r})$ die Diracsche Deltadistribution bezeichnet:

$A_{1}(\vec{r}\,) = (A * \delta)(\vec r\,) = \int A(\vec{r}\,') \delta (\vec{r} - \vec{r}\,') \mathrm{d} \vec {r}\,' \approx \int A(\vec{r}\,') W(\vec{r} - \vec{r}\,', h) \mathrm d\vec{r}\,'$

Dann wird die δ-Distribution durch einen Kern $W(\vec{r} - \vec{r}\,', h)$ angenähert, wobei h die Glättungslänge ist. Damit die Näherung im Grenzfall gültig bleibt, kann man Normierung und Identität mit der δ-Distribution im Grenzwert für h → 0 fordern:

$\int W(\vec{r} - \vec{r}\,', h)\mathrm d\vec{r}\,' = 1$

$\lim \limits_{h \to 0} W(\vec{r} - \vec{r}\,', h) = \delta (\vec{r} - \vec{r}\,')$

Tatsächlich ist dies bei den meisten verwendeten Kernen nicht mehr der Fall. Um daraus die Aufteilung in Massenelemente zu erhalten, erweitert man mit der Dichte ρ und belässt h größer als 0. Für den Fall unendlich vieler, unendlich kleiner Teilchen geht die Summe in das Integral über. Numerisch wird man sich immer mit endlich vielen Teilchen zufriedengeben müssen:

$A_{S}(\vec{r}\,) = \lim \limits_{h \rightarrow \infty} \int \frac{A(\vec{r}\,')}{\rho(\vec{r}\,')} W(\vec{r} - \vec{r}\,', h) \rho(\vec{r}\,') \mathrm d\vec{r}\,' \propto \sum \limits_{b= 1}^{N} m_{b} \frac{A_{b}}{\rho_{b}} W(\vec{r} - \vec{r}\,', h) = A_{b}$

Dabei ist m_b die Masse des Teilchens b und ρ_b die Dichte am Ort des Teilchens b:

$\rho_b = \rho(\mathbf{r}_b) = \sum_j m_j \frac{\rho_j}{\rho_j} W(| \mathbf{r}_b-\mathbf{r}_j |,h) = \sum_j m_j W(\mathbf{r}_b-\mathbf{r}_j,h)$

Damit haben wir die Grundgleichung der Smoothed Particle Hydrodynamics hergeleitet (rechter Teil). Die Größe A wird durch eine Summe über alle Teilchen berechnet. Man sieht, dass aus der von r abhängigen Größe A_S ein Skalar A_b multipliziert mit dem Kernel geworden ist. Dies führt zu einer starken Vereinfachung von Differentialgleichungen, da eine Ableitung nun nicht mehr auf die Größe, sondern nur noch auf den Kernel wirkt:

$\nabla A(\vec{r}\,) = \sum \limits_{b} m_{b} \frac{A_{b}}{\rho_{b}} \nabla W(\vec{r} - \vec{r}\,', h)$

Kern und Glättungslänge (Smoothing Length)

Glättungslänge

Der wohl wichtigste Parameter der SPH ist die Glättungslänge h. Sie legt die Auflösung der Methode fest und hat damit starken Einfluss auf Genauigkeit und Rechenaufwand bei Simulationen. Bei entsprechender Wahl des Kernels (siehe unten) legt sie auch die Anzahl der bei Berechnung mit einzubeziehenden Nachbarn fest. Üblich bis zu einige Zehn Teilchen pro Größe. Für gute Ergebnisse orientiert man sich an der mittleren Dichte des Fluids :

$h \sim \frac{1}{\langle \rho \rangle^{\frac{1}{\nu}}}$

mit n Teilchen, ν Dimensionen und

$\langle \rho \rangle = \frac{1}{n} \sum \limits_{b} \rho_{b}$

In modernen Codes wählt man h = h(t) zeitabhängig. Mit

$\frac{\mathrm d h_{a}}{\mathrm dt} = - \left( \frac{h_{a}}{\nu \rho_{a}} \right)\frac{\mathrm d \rho_{a}}{\mathrm dt}$

nutzt man dann in Gebieten großer Dichten eine höhere Auflösung, während in Bereichen geringer Dichten die Smoothing Length kleiner wird. Dadurch lässt sich der Rechenaufwand bei gleich bleibender Genauigkeit verringern.

Kern

Der Kern ist die wohl wichtigste Struktur der SPH-Methode. Verschiedene Kerne entsprechen verschiedenen Differenzenschemata in Gittermethoden. Zur Interpretation von SPH-Gleichungen ist es vorteilhaft, einen Kern in Form einer gaußschen Kurve zu verwenden:

$W(\vec{r} - \vec{r}\,',h) \,\sim\, e^{-(\frac{\vec{r} - \vec{r}\,'}{h})^2}$

Numerisch ist dieser Ansatz allerdings nicht sehr geeignet, da man in diesem Fall oft auf ein klares Verhalten bezüglich der Reichweite des Kerns Wert legt. D. h. man wählt einen Kern, der ab einem gewissen $\vec{r}$ null ist, um die Anzahl der Nachbarn, die bei der Berechnung mit einbezogen werden, klar festlegen zu können. Damit kann man den benötigten Rechenaufwand eingrenzen. Wie bereits erwähnt, ist SPH eine sehr empirische Methode, d. h. für unterschiedliche Anwendungen werden sehr unterschiedliche Kerne benötigt. Die genaue Wahl ist Erfahrungssache und erfolgt oft nach dem Versuch und Irrtum-Prinzip. Da ein Kern oft in einer eigenen Funktion implementiert wird, ist der Aufwand ihn auszutauschen oder zu verändern minimal. Oft werden Kerne auf Basis von Splines verwendet:

$W(\vec{r}, h) = \frac{\sigma}{4 h^{\nu}} \cdot \begin{cases} \left( 4-6 q^{2}+3 q^{3} \right) & \quad \mathrm{f\ddot ur} \quad 0 \leqslant \frac{r}{h} \leqslant 1 \\ \left( 2 - q \right)^{3} & \quad \mathrm{f\ddot ur} \quad 1 \leqslant \frac{r}{h} \leqslant 2 \\ 0 & \quad \mathrm{sonst} \end{cases}$

Mit $q = \vec{r} - \vec{r}\,'$, einer Normierungskonstante σ und der Anzahl der Dimensionen ν. Hier werden nur Teilchen bis zum übernächsten Nachbarn in die Berechnung mit einbezogen. Außerdem ist die 2. Ableitung dieses Kerns nicht konstant, weshalb er nicht von der Unordnung der Teilchen abhängt.

Fehlerabschätzungen

Bei der Herleitung über Integralinterpolationsfunktionen wurden zwei Näherungen gemacht. Erstens wurde h > 0 angenommen und die Summation erfolgt nur über eine endliche Zahl von Teilchen.

• Für die Identität, d. h. mit h = 0 und beliebig vielen Teilchen, gibt eine Taylorentwicklung einen Fehler von O(h²)
• Auch für die Summationsnäherung kann man mit Hilfe der Shoenberg-Formel einen Fehler O(h²) ausrechnen, falls die Teilchen geordnet im Fluid verteilt sind
• Im Falle von ungeordneten Teilchen existiert keine traditionelle Fehlerabschätzung

Damit ist man bei Simulationen mit SPH immer auf den Vergleich mit anderen Simulationen angewiesen, zumindest für eine Fehlereinschätzung. Einige Veröffentlichungen erwähnen, dass die Fehler meist deutlich unten denen einer Monte-Carlo-Simulation liegen, auch dies ist Erfahrungssache. Generell neigt SPH zur Ausschmierung von Diskontinuitäten, ist also gerade im Falle von Simulationen mit wenigen Teilchen lokal recht ungenau. Für große Teilchenzahlen wird das Verhalten aber deutlich besser. Allerdings ist das globale Verhalten schon bei geringen Teilchenzahlen, was geringem Rechenaufwand entspricht, sehr gut. D. h. globale Größen wie die Energie sind gut wiedergegeben. Oft lässt sich mit SPH eine global gute Simulation mit wenig Aufwand programmieren, die in akzeptabler Zeit auf Workstations gerechnet werden kann.

Vorteile und Nachteile

Vorteile:

• SPH ist eine Lagrange-Methode; die Kontinuitätsgleichung ist automatisch erfüllt.
• Der Code ist sehr robust, d. h. liefert fast immer Ergebnisse
• Die Implementation von SPH ist vergleichsweise einfach, ebenso das Testen verschiedener Kernels.
• Mit Hilfe einer Gaußfunktion als Kernel lassen sich theoretische Ergebnisse leicht interpretieren.
• In modernem Code zeigt sich eine NlogN Abhängigkeit des Rechenaufwandes von der Teilchenzahl.
• SPH zeigt gute globale Ergebnisse bei geringen Teilchenzahlen.

Nachteile:

• Der Code ist oft zu robust; trotz eines falschen Modells kann SPH Ergebnisse liefern, die dann aber physikalisch inkorrekt sind
• Die Fehlerabschätzung ist oft problematisch und nur im Vergleich mit den Ergebnissen anderer Methoden zu erhalten
• Die Methode ist hoch dispersiv
• Für gute Genauigkeiten werden hohe Teilchenzahlen benötigt. Damit ist der Vorteil geringen Rechenaufwandes nicht zutreffend
• Die Behandlung von Diskontinuitäten ist oft schwierig, da Strukturen auf Skalen, die kleiner als die Smoothing length sind, geglättet werden.

Hydrodynamische Gleichungen in SPH

Symmetrisierung

Um die Hydrodynamik in SPH zu formulieren, ist der scheinbar einfachste Ansatz die Grundgleichung in die hydrodynamischen Gleichungen wie z. B. die Navier-Stokes-Gleichung einzusetzen. Die daraus resultierenden Gleichungen sind allerdings nicht symmetrisch gegenüber Teilchenvertauschung. Deshalb gelten in diesem Fall viele Erhaltungssätze für Energie, Drehimpuls, etc., nicht mehr. Oft ist es allerdings möglich, diese zu retten, indem man die Dichte in den jeweiligen Differentialoperator herein schreibt und die Produktregel nutzt:

$\rho \nabla A = \nabla(\rho A ) - A \nabla \rho$

Oft lassen sich so symmetrische Gleichungen herleiten. All dies geschieht nicht streng formal, sondern nur, weil es bessere Ergebnisse liefert.

Bewegung des Fluids

Die einfachste Möglichkeit ist die Verwendung der Definition der Geschwindigkeit:

$\frac{\mathrm d \vec{r}_{a}}{\mathrm dt} = \vec{v}_{a}$

Dabei ist die Bewegung eines Teilchens nicht an die der anderen gekoppelt, was oft zu Problemen führen kann. Deshalb hat man die XPSH-Methode entwickelt:

$\frac{\mathrm d \vec{r}_{a}}{\mathrm dt} = \vec{v}_{a} + \varepsilon \sum \limits_{b} m_{b} \left( \frac{\vec{v}_{ba}}{\bar{\rho}_{ab}} \right) W_{ab}$

mit einer gemittelten Dichte:

$\bar{\rho}_{ab} = \frac{1}{2}\left( \rho_{a} + \rho_{b} \right)$

und einem Kopplungsparameter ε. Damit wird die Ordnung der Teilchen besser erhalten ohne dass zusätzlich Viskosität eingeführt werden muss.

Kontinuitätsgleichung in SPH

Setzen wir die Dichte in die Grundgleichung ein, so erhalten wir

$\rho_{a} = \sum \limits_{b} m_{b} W_{ab}$

für ein Teilchen a. Daraus lässt sich die SPH-Kontinuitätsgleichung ausrechnen

$\frac{\mathrm d\rho_{a}}{\mathrm dt} = - \sum \limits_{b} m_{b} v_{ab}$

Euler-Gleichung in SPH

Für die Euler-Gleichung ergibt sich:

$\frac{\mathrm d \vec{v}_{a}}{\mathrm dt} = - \frac{1}{\rho_{a}} \sum \limits_{b} m_{b} \frac{P_{b}}{\rho_{b}} \nabla_{a} W_{ab}$

Diese Gleichung ist nicht symmetrisch gegenüber Teilchenaustausch: Impuls und Drehmoment sind nicht erhalten. Deswegen verwenden wir den oben angedeuteten Trick für den Druckgradienten:

$\frac{\nabla P}{\rho} = \nabla \left( \frac{P}{\rho} \right) + \frac{P}{\rho^{2}} \nabla \rho$

Woraus wir die gewünschte symmetrische Gleichung erhalten:

$\frac{\mathrm d \vec{v}_{a}}{\mathrm dt} = - \sum \limits_{b} m_{b} \left( \frac{P_{b}}{\rho_{b}^{2}} + \frac{P_{a}}{\rho_{a}^{2}} \right) \nabla_{a} W_{ab}$

Setzen wir einen Gauß-Funktion ein ergibt sich eine Zentralkraft, die auf beide Teilchen gleich stark wirkt:

$m_{a} \frac{\mathrm d \vec{v}_{a}}{\mathrm dt} = \frac{2m_{a}m_{b}}{h^{2}} \left( \frac{P_{b}}{\rho_{b}^{2}} + \frac{P_{a}}{\rho_{a}^{2}} \right) \left( \vec{r}_{a} - \vec{r}_{b} \right) W_{ab}$

Viskosität

Wie fast jede numerische Methode erzeugt auch SPH durch Rechenungenauigkeiten Viskosität. Zur Modellierung ist diese oftmals aber nicht ausreichend. Deswegen führt man, ähnlich wie beim Übergang von der Eulergleichung zur Navier-Stokes-Gleichung, einen Viskositätstensor ein. Die genaue Wahl dieses Tensors hängt stark vom Modell ab.

$\frac{\mathrm d \vec{v}_{a}}{\mathrm dt} = - \sum \limits_{b} m_{b} \left( \frac{P_{b}}{\rho_{b}^{2}} + \frac{P_{a}}{\rho_{a}^{2}} + \Pi_{ab} \right) \nabla_{a} W_{ab}$

Anwendungen

SPH findet in vielen verschiedenen Bereichen wie der Astrophysik Anwendung. Es existieren auch relativistische und magnetische SPH-Methoden:

• Gasdynamik
• Galaxie-Entstehung und -Verschmelzung
• Binäre Sternsysteme, Akkretionsscheiben und Sternkollisionen
• Mondentstehung
• Relativistische Probleme
• Magnetische Probleme

Weiterführende Veröffentlichungen

• Monaghan: Smoothed Particle Hydrodynamics; Annu. Rev. Astrophys. 1992
• Steinmetz, Müller: On the capabilities and limits of s.p.h.; Astronomy and Astrophysics 1993
• Alimi, Courty: Thermodynamic evolution of the cosmological baryonic gas pt.2; Astronomy and Astrophysics 2005

Weblinks

Filme

• Visualisierungen von Galaxien: Entstehungsrechnungen (auf Englisch)
• Galaxien: Entstehung und Verschmelzungen (auf Englisch)

Code

• Website des Gadget Codes, dessen Gasdynamic SPH nutzt (auf Englisch)
• SimPARTIX ist ein Softwarepaket für SPH- und DEM-Simulationen vom Fraunhofer IWM
• Pasimodo simuliert verschiedene Partikelmethoden, u.a. SPH
• SOFA - Simulation Open Framework Architecture simuliert SPH im Hinblick auf Echtzeitanwendungen.
• OpenTissue bietet ebenfalls SPH mit Fokus auf Echtzeitanwendung an.