Smoothed Particle Hydrodynamics: Unterschied zwischen den Versionen

Smoothed Particle Hydrodynamics: Unterschied zwischen den Versionen

imported>Debenben
 
imported>Blaues-Monsterle
(Die 2 letzten Textänderungen von Bo Zhang tum wurden verworfen und die Version 205949006 von Orthographus wiederhergestellt. schaut mir etwas nach Werbung aus)
 
Zeile 1: Zeile 1:
'''Smoothed-particle hydrodynamics''' ('''SPH'''; deutsch: ''geglättete Teilchen-Hydrodynamik'') ist eine numerische Methode, um die [[ hydrodynamik|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.
'''Smoothed-particle hydrodynamics''' ('''SPH'''; deutsch: ''geglättete Teilchen-Hydrodynamik'') ist eine numerische Methode, um die [[Hydrodynamik|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 bewegen sich mit dem Fluid mit. SPH ist eine besonders einfach zu implementierende und robuste Methode.


== Methode ==
== Methode ==
=== Allgemeines ===
=== Allgemeines ===
In Smoothed Particle Hydrodynamics wird die zu simulierende Flüssigkeit in Elemente aufgeteilt. Dabei werden, ähnlich den [[Monte-Carlo-Simulation|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) <math> h </math> repräsentiert. Sie ist der wichtigste Parameter der Methode. Zwischen den Teilchen wird das Fluid durch den Kernel geglättet, daher der Name.
In Smoothed Particle Hydrodynamics wird die zu simulierende Flüssigkeit in Elemente aufgeteilt. Dabei werden, ähnlich den [[Monte-Carlo-Simulation|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) <math> h </math> 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.&nbsp;B. die Dichte <math> \rho </math>) 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.
Jede Größe (z.&nbsp;B. die Dichte <math> \rho </math>) 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 [[Empirie|empirische]] Methode. Das bedeutet, dass viele Dinge gemacht werden, weil sie funktionieren, nicht, weil es eine strenge mathematische Herleitung gibt.
SPH ist eine sehr [[Empirie|empirische]] Methode. Das bedeutet, dass viele Dinge gemacht werden, weil sie funktionieren, nicht, weil es eine strenge mathematische Herleitung gibt.
Zeile 9: Zeile 9:
=== Herleitung ===
=== Herleitung ===
Die formale Herleitung läuft entweder über eine Lagrange-Funktion oder über eine Integralinterpolation.
Die formale Herleitung läuft entweder über eine Lagrange-Funktion oder über eine Integralinterpolation.
Bei der Integralinterpolation für eine Größe <math> A_{1} </math>   geht man von einer Identität aus, wobei <math> \delta (\vec{r})</math> die [[Delta-Distribution|Diracsche Deltadistribution]] bezeichnet:
Bei der Integralinterpolation für eine Größe <math> A_{1} </math> geht man von einer Identität aus, wobei <math> \delta (\vec{r})</math> die [[Delta-Distribution|Diracsche Deltadistribution]] bezeichnet:


:<math>
:<math>
Zeile 25: Zeile 25:


Tatsächlich ist dies bei den meisten verwendeten Kernen nicht mehr der Fall.
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 <math>\rho</math> und belässt <math> h </math> 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:
Um daraus die Aufteilung in Massenelemente zu erhalten, erweitert man mit der Dichte <math>\rho</math> und belässt <math> h </math> 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:


:<math>
:<math>
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}
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}
</math>
</math>


Dabei ist <math> m_{b} </math> die Masse des Teilchens b und <math>\rho_{b} </math> die Dichte am Ort des Teilchens b:  
Dabei ist <math> m_{b} </math> die Masse des Teilchens b und <math>\rho_{b} </math> die Dichte am Ort des Teilchens b:
:<math>
: <math>
\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)
\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)
</math>
</math>
Zeile 38: Zeile 38:
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 <math> A_{S}</math> ein Skalar <math> A_{b}</math> 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:
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 <math> A_{S}</math> ein Skalar <math> A_{b}</math> 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:


:<math>
: <math>
\nabla A(\vec{r}\,) = \sum \limits_{b} m_{b} \frac{A_{b}}{\rho_{b}} \nabla W(\vec{r} - \vec{r}\,', h)
\nabla A(\vec{r}\,) = \sum \limits_{b} m_{b} \frac{A_{b}}{\rho_{b}} \nabla W(\vec{r} - \vec{r}\,', h)
</math>
</math>
Zeile 45: Zeile 45:
==== Glättungslänge ====
==== Glättungslänge ====
Der wohl wichtigste Parameter der SPH ist die Glättungslänge <math>h</math>. 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 sind bis zu einige zehn Teilchen pro Größe.
Der wohl wichtigste Parameter der SPH ist die Glättungslänge <math>h</math>. 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 sind bis zu einige zehn Teilchen pro Größe.
Für gute Ergebnisse orientiert man sich an der mittleren Dichte des Fluids :
Für gute Ergebnisse orientiert man sich an der mittleren Dichte des Fluids:


:<math>
: <math>
h \sim \frac{1}{\langle \rho \rangle^{\frac{1}{\nu}}}
h \sim \frac{1}{\langle \rho \rangle^{\frac{1}{\nu}}}
</math>
</math>
mit <math> n</math> Teilchen, <math> \nu</math> Dimensionen und
mit <math> n</math> Teilchen, <math> \nu</math> Dimensionen und


:<math>
: <math>
\langle \rho \rangle = \frac{1}{n} \sum \limits_{b} \rho_{b}
\langle \rho \rangle = \frac{1}{n} \sum \limits_{b} \rho_{b}
</math>
</math>


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


:<math>
: <math>
\frac{\mathrm d h_{a}}{\mathrm dt} = - \left( \frac{h_{a}}{\nu \rho_{a}} \right)\frac{\mathrm d \rho_{a}}{\mathrm dt}
\frac{\mathrm d h_{a}}{\mathrm dt} = - \left( \frac{h_{a}}{\nu \rho_{a}} \right)\frac{\mathrm d \rho_{a}}{\mathrm dt}
</math>
</math>
Zeile 68: Zeile 68:
Zur Interpretation von SPH-Gleichungen ist es vorteilhaft, einen Kern in Form einer gaußschen Kurve zu verwenden:
Zur Interpretation von SPH-Gleichungen ist es vorteilhaft, einen Kern in Form einer gaußschen Kurve zu verwenden:


:<math>
: <math>
W(\vec{r} - \vec{r}\,',h) \,\sim\, e^{-(\frac{\vec{r} - \vec{r}\,'}{h})^2}
W(\vec{r} - \vec{r}\,',h) \,\sim\, e^{-(\frac{\vec{r} - \vec{r}\,'}{h})^2}
</math>
</math>


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.&nbsp;h. man wählt einen Kern, der ab einem gewissen <math>\vec{r} </math> 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.
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.&nbsp;h. man wählt einen Kern, der ab einem gewissen <math>\vec{r} </math> 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.&nbsp;h. für unterschiedliche Anwendungen werden sehr unterschiedliche Kerne benötigt. Die genaue Wahl ist Erfahrungssache und erfolgt oft nach dem [[Versuch und Irrtum|Versuch und Irrtum]]-Prinzip. Da ein Kern oft in einer eigenen Funktion implementiert wird, ist der Aufwand ihn auszutauschen oder zu verändern minimal.
Wie bereits erwähnt, ist SPH eine sehr empirische Methode, d.&nbsp;h. für unterschiedliche Anwendungen werden sehr unterschiedliche Kerne benötigt. Die genaue Wahl ist Erfahrungssache und erfolgt oft nach dem [[Versuch und Irrtum|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 [[Spline|Splines]] verwendet:
Oft werden Kerne auf Basis von [[Spline]]s verwendet:


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


Mit <math> q = \vec{r} - \vec{r}\,'</math>, einer Normierungskonstante <math> \sigma</math> und der Anzahl der Dimensionen <math> \nu </math>.
Mit <math> q = \vec{r} - \vec{r}\,'</math>, einer Normierungskonstante <math> \sigma</math> und der Anzahl der Dimensionen <math> \nu </math>.
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.
Hier werden nur Teilchen bis zum übernächsten Nachbarn in die Berechnung mit einbezogen. Außerdem ist die 2.&nbsp;Ableitung dieses Kerns nicht konstant, weshalb er nicht von der Unordnung der Teilchen abhängt.


=== Fehlerabschätzungen ===
=== Fehlerabschätzungen ===
Zeile 102: Zeile 102:


* SPH ist eine Lagrange-Methode; die Kontinuitätsgleichung ist automatisch erfüllt.
* SPH ist eine Lagrange-Methode; die Kontinuitätsgleichung ist automatisch erfüllt.
* Der Code ist sehr robust, d. h. liefert fast immer Ergebnisse
* Der Code ist sehr robust, d.&nbsp;h. liefert fast immer Ergebnisse
* Die Implementation von SPH ist vergleichsweise einfach, ebenso das Testen verschiedener Kernels.
* 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.
* Mit Hilfe einer Gaußfunktion als Kernel lassen sich theoretische Ergebnisse leicht interpretieren.
* In modernem Code zeigt sich eine <math> N \log N</math> Abhängigkeit des Rechenaufwandes von der Teilchenzahl.
* In modernem Code zeigt sich eine <math> N \log N</math> Abhängigkeit des Rechenaufwandes von der Teilchenzahl.
* SPH zeigt gute globale Ergebnisse bei geringen Teilchenzahlen.
* SPH zeigt gute globale Ergebnisse bei geringen Teilchenzahlen.


Zeile 118: Zeile 118:
== Hydrodynamische Gleichungen in SPH ==
== Hydrodynamische Gleichungen in SPH ==
=== Symmetrisierung ===
=== Symmetrisierung ===
Um die Hydrodynamik in SPH zu formulieren, ist der scheinbar einfachste Ansatz die Grundgleichung in die hydrodynamischen Gleichungen wie z.&nbsp;B. die [[Navier-Stokes-Gleichungen|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.
Um die Hydrodynamik in SPH zu formulieren, ist der scheinbar einfachste Ansatz die Grundgleichung in die hydrodynamischen Gleichungen wie z.&nbsp;B. die [[Navier-Stokes-Gleichungen|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:
Oft ist es allerdings möglich, diese zu retten, indem man die Dichte in den jeweiligen Differentialoperator herein schreibt und die Produktregel nutzt:


:<math>
: <math>
\rho \nabla A = \nabla(\rho A ) - A \nabla \rho
\rho \nabla A = \nabla(\rho A ) - A \nabla \rho
</math>
</math>
Zeile 130: Zeile 130:
Die einfachste Möglichkeit ist die Verwendung der Definition der Geschwindigkeit:
Die einfachste Möglichkeit ist die Verwendung der Definition der Geschwindigkeit:


:<math>
: <math>
\frac{\mathrm d \vec{r}_{a}}{\mathrm dt} = \vec{v}_{a}
\frac{\mathrm d \vec{r}_{a}}{\mathrm dt} = \vec{v}_{a}
  </math>
  </math>
Zeile 136: Zeile 136:
Dabei ist die Bewegung eines Teilchens nicht an die der anderen gekoppelt, was oft zu Problemen führen kann. Deshalb hat man die XSPH-Methode ("Extended SPH") entwickelt:
Dabei ist die Bewegung eines Teilchens nicht an die der anderen gekoppelt, was oft zu Problemen führen kann. Deshalb hat man die XSPH-Methode ("Extended SPH") entwickelt:


:<math>
: <math>
\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}
\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}
</math>
</math>


mit einer gemittelten Dichte:
mit einer gemittelten Dichte:
:<math> \bar{\rho}_{ab} = \frac{1}{2}\left( \rho_{a} + \rho_{b} \right) </math>
:<math> \bar{\rho}_{ab} = \frac{1}{2}\left( \rho_{a} + \rho_{b} \right) </math>
und einem Kopplungsparameter&nbsp;ε. Damit wird die Ordnung der Teilchen besser erhalten ohne dass zusätzlich Viskosität eingeführt werden muss.
und einem Kopplungsparameter&nbsp;ε. Damit wird die Ordnung der Teilchen besser erhalten, ohne dass zusätzlich Viskosität eingeführt werden muss.


=== Kontinuitätsgleichung in SPH ===
=== Kontinuitätsgleichung in SPH ===
Setzen wir die Dichte in die Grundgleichung ein, so erhalten wir
Setzen wir die Dichte in die Grundgleichung ein, so erhalten wir


:<math>
: <math>
\rho_{a} = \sum \limits_{b} m_{b} W_{ab}
\rho_{a} = \sum \limits_{b} m_{b} W_{ab}
</math>
</math>
für ein Teilchen a. Daraus lässt sich die SPH-[[Kontinuitätsgleichung]] ausrechnen
für ein Teilchen a. Daraus lässt sich die SPH-[[Kontinuitätsgleichung]] ausrechnen


:<math>
: <math>
\frac{\mathrm d\rho_{a}}{\mathrm dt} = - \sum \limits_{b} m_{b} v_{ab} \nabla_{a} W_{ab}
\frac{\mathrm d\rho_{a}}{\mathrm dt} = - \sum \limits_{b} m_{b} v_{ab} \nabla_{a} W_{ab}
</math>
</math>


=== Euler-Gleichung in SPH ===
=== Euler-Gleichung in SPH ===
Für die [[Eulersche Gleichungen (Strömungsmechanik)|Euler-Gleichung]] ergibt sich:
Für die [[Euler-Gleichungen (Strömungsmechanik)|Euler-Gleichung]] ergibt sich:


:<math>
: <math>
\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}
\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}
</math>
</math>
Zeile 165: Zeile 165:
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:
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:


:<math>
: <math>
\frac{\nabla P}{\rho} = \nabla \left( \frac{P}{\rho} \right) + \frac{P}{\rho^{2}} \nabla \rho
\frac{\nabla P}{\rho} = \nabla \left( \frac{P}{\rho} \right) + \frac{P}{\rho^{2}} \nabla \rho
</math>
</math>
Zeile 171: Zeile 171:
Woraus wir die gewünschte symmetrische Gleichung erhalten:
Woraus wir die gewünschte symmetrische Gleichung erhalten:


:<math>
: <math>
\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}
\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}
</math>
</math>


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


:<math>
: <math>
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}
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}
</math>
</math>


=== Viskosität ===
=== 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 [[Eulersche Gleichungen (Strömungsmechanik)|Euler-Gleichung]] zur [[Navier-Stokes-Gleichung]], einen Viskositätstensor ein. Die genaue Wahl dieses Tensors hängt stark vom Modell ab.
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 [[Euler-Gleichungen (Strömungsmechanik)|Euler-Gleichung]] zur [[Navier-Stokes-Gleichung]], einen Viskositätstensor ein. Die genaue Wahl dieses Tensors hängt stark vom Modell ab.


:<math>
: <math>
\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}
\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}
</math>
</math>


== Anwendungen ==
== Anwendungen ==
SPH findet in vielen verschiedenen Bereichen wie der [[Astrophysik]] Anwendung. Es existieren auch relativistische und magnetische SPH-Methoden:
SPH wird in vielen verschiedenen Bereichen wie der [[Astrophysik]] angewendet. Es existieren auch relativistische und magnetische SPH-Methoden:


* [[Gasdynamik]]
* [[Gasdynamik]]
Zeile 205: Zeile 205:
== Weblinks ==
== Weblinks ==
=== Filme ===
=== Filme ===
* [http://www.mpa-garching.mpg.de/galform/data_vis/index.shtml Visualisierungen von Galaxien: Entstehungsrechnungen] (auf Englisch)
* [http://www.mpa-garching.mpg.de/galform/data_vis/index.shtml Visualisierungen von Galaxien: Entstehungsrechnungen] (englisch)
* [http://www.aip.de/People/MSteinmetz/Movies.html Galaxien: Entstehung und Verschmelzungen] (auf Englisch)
* {{Internetquelle
  |autor=Mathias Steinmetz
  |url=https://escience.aip.de/vis/category/movie/msteinmetz/
  |titel=Galaxien: Entstehung und Verschmelzungen
  |hrsg=Leibniz Institute for Astrophysics Potsdam
  |sprache=en
  |abruf=2018-02-26
  |abruf-verborgen=1}}


=== Code ===
=== Code ===
* [http://www.mpa-garching.mpg.de/galform/gadget/index.shtml Website des Gadget Codes], dessen Gasdynamic SPH nutzt (auf Englisch)
* [http://www.mpa-garching.mpg.de/galform/gadget/index.shtml Website des Gadget Codes], dessen Gasdynamic SPH nutzt (englisch)
* [http://www.simpartix.de SimPARTIX] ist ein Softwarepaket für SPH- und DEM-Simulationen vom Fraunhofer IWM
* [http://www.simpartix.de/ SimPARTIX] ist ein Softwarepaket für SPH- und DEM-Simulationen vom Fraunhofer IWM
*[http://www.itm.uni-stuttgart.de/research/pasimodo/pasimodo_de.php Pasimodo] simuliert verschiedene Partikelmethoden, u.a. SPH
* [http://www.itm.uni-stuttgart.de/research/pasimodo/pasimodo_de.php Pasimodo] simuliert verschiedene Partikelmethoden, u.&nbsp;a. SPH
*[http://www.sofa-framework.org/features SOFA - Simulation Open Framework Architecture] simuliert SPH im Hinblick auf Echtzeitanwendungen.
* [https://www.becker3d.com/ ThreeParticle/CAE] ist eine multiphysics Softwareplattform für DEM und SPH Simulation der BECKER 3D GmbH
*[http://www.opentissue.org/OpenTissue/media.php OpenTissue] bietet ebenfalls SPH mit Fokus auf Echtzeitanwendung an.
* {{Webarchiv |url=http://www.sofa-framework.org/features |wayback=20080408033019 |text=SOFA Simulation Open Framework Architecture}} simuliert SPH im Hinblick auf Echtzeitanwendungen.
* {{Webarchiv |url=http://www.opentissue.org/OpenTissue/media.php |wayback=20110824200803 |text=OpenTissue}} bietet ebenfalls SPH mit Fokus auf Echtzeitanwendung an.


[[Kategorie:Strömungsmechanik]]
[[Kategorie:Strömungsmechanik]]
[[Kategorie:Numerische Mathematik]]
[[Kategorie:Numerische Mathematik]]
[[Kategorie:Computerphysik]]
[[Kategorie:Computerphysik]]

Aktuelle Version vom 9. April 2021, 19:48 Uhr

Smoothed-particle hydrodynamics (SPH; deutsch: 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 bewegen sich mit dem Fluid mit. SPH ist eine besonders einfach zu implementierende und robuste Methode.

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 $ \rho $) 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_{1} $ 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 $ \delta $-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 $ \delta $-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 $ \rho $ 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 $ \rho _{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 sind 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, $ \nu $ 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 {d} t}}=-\left({\frac {h_{a}}{\nu \rho _{a}}}\right){\frac {\mathrm {d} \rho _{a}}{\mathrm {d} t}} $

nutzt man dann in Gebieten großer Dichten eine höhere Auflösung, während in Bereichen geringer Dichten die Smoothing Length größer 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 }{4h^{\nu }}}\cdot {\begin{cases}\left(4-6q^{2}+3q^{3}\right)&\quad \mathrm {f{\ddot {u}}r} \quad 0\leqslant {\frac {r}{h}}\leqslant 1\\\left(2-q\right)^{3}&\quad \mathrm {f{\ddot {u}}r} \quad 1\leqslant {\frac {r}{h}}\leqslant 2\\0&\quad \mathrm {sonst} \end{cases}} $

Mit $ q={\vec {r}}-{\vec {r}}\,' $, einer Normierungskonstante $ \sigma $ und der Anzahl der Dimensionen $ \nu $. 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^{2}) $.
  • Auch für die Summationsnäherung kann man mit Hilfe der Shoenberg-Formel einen Fehler $ O(h^{2}) $ 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 unter 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 $ N\log N $ 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 {d} t}}={\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 XSPH-Methode ("Extended SPH") entwickelt:

$ {\frac {\mathrm {d} {\vec {r}}_{a}}{\mathrm {d} t}}={\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 {d} t}}=-\sum \limits _{b}m_{b}v_{ab}\nabla _{a}W_{ab} $

Euler-Gleichung in SPH

Für die Euler-Gleichung ergibt sich:

$ {\frac {\mathrm {d} {\vec {v}}_{a}}{\mathrm {d} t}}=-{\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 {d} t}}=-\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 {d} t}}={\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 Euler-Gleichung 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 {d} t}}=-\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 wird in vielen verschiedenen Bereichen wie der Astrophysik angewendet. 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

Code

  • Website des Gadget Codes, dessen Gasdynamic SPH nutzt (englisch)
  • SimPARTIX ist ein Softwarepaket für SPH- und DEM-Simulationen vom Fraunhofer IWM
  • Pasimodo simuliert verschiedene Partikelmethoden, u. a. SPH
  • ThreeParticle/CAE ist eine multiphysics Softwareplattform für DEM und SPH Simulation der BECKER 3D GmbH
  • SOFA – Simulation Open Framework Architecture (Memento vom 8. April 2008 im Internet Archive) simuliert SPH im Hinblick auf Echtzeitanwendungen.
  • OpenTissue (Memento vom 24. August 2011 im Internet Archive) bietet ebenfalls SPH mit Fokus auf Echtzeitanwendung an.