SPH

SPH (ang. smoothed particle hydrodynamics) – metoda numeryczna służąca do przeprowadzania symulacji numerycznych zachowania się płynów. W metodzie tej ciecz jest reprezentowana przez zbiór dyskretnych cząstek (bądź „pseudocząstek”), które oddziałują ze sobą poruszając się wraz z przepływem i niosąc ze sobą wszystkie informacje o płynie^[1].

Metoda ta została zaproponowana w roku 1977 niezależnie przez R.A. Gingolda & J.J. Monaghana oraz przez L.B.Lucy do prowadzenia obliczeń z dziedziny astrofizyki. Początkowo używana była do symulacji ściśliwych cieczy nielepkich. Z czasem została rozwinięta na nieściśliwe ciecze lepkie znajdujące się w polu grawitacyjnym, a także na zagadnienia z dziedziny magnetohydrodynamiki.

Równania ruchu

W metodzie SPH do opisu stanu cieczy używa się opisu Lagrange’a, gdzie siatka obliczeniowa porusza się wraz z przepływem cieczy. W takim wypadku równanie Naviera-Stokesa dla ${\displaystyle i}$-tej cząstki przybiera postać

${\displaystyle {\frac {d{\vec {v}}}{dt}}=-{\frac {\nabla p}{\rho }}+{\vec {a}}^{visc}+{\vec {\Phi }},}$

gdzie:

${\displaystyle {\vec {v}}}$ – prędkość,

${\displaystyle p}$ – ciśnienie,

${\displaystyle \rho }$ – gęstość,

${\displaystyle {\vec {a}}^{visc}}$ – przyspieszenie wynikające z istnienia sił lepkości,

${\displaystyle {\vec {\Phi }}}$ – przyspieszenie wynikające z obecności sił masowych (np. pola grawitacyjnego).

Metoda SPH

Podstawy

Metoda ta opiera się na teorii interpolacji. Ciągłe rozkłady takich parametrów jak gęstość czy ciśnienie cieczy zastępuje się odpowiednimi estymatami przy założonym pewnym jądrze interpolacji. Obliczenia wykonujemy dla dyskretnego zbioru ${\displaystyle N}$ cząstek płynu.

Estymata pewnej wielkości ${\displaystyle A}$ w pozycji ${\displaystyle i}$-tej cząstki jest dana jako

${\displaystyle A_{i}=\sum _{j=1}^{N}m_{j}{\frac {A_{j}}{\rho _{j}}}W_{ij}.}$

Natomiast estymata gradientu wielkości ${\displaystyle A}$ jako

${\displaystyle \nabla _{i}A_{i}=\sum _{j=1}^{N}m_{j}{\frac {A_{j}}{\rho _{j}}}\nabla _{i}W_{ij},}$

gdzie:

${\displaystyle \rho _{j}}$ – gęstość ${\displaystyle j}$-tej cząstki,

${\displaystyle m_{j}}$ – masa ${\displaystyle j}$-tej cząstki,

${\displaystyle W_{ij}=W(r,h)}$ – jądro interpolacji.

Ilość sąsiadów i długość wygładzania

Parametr ${\displaystyle h}$ jest nazywany długością wygładzania (smoothing length). Jest to wielkość, która określa na jaką odległość cząstka może oddziaływać z innymi cząstkami. Najczęściej w symulacji pozostaje stała podczas trwania obliczeń. Należy tylko uwzględnić, aby w promieniu ${\displaystyle 2h}$ znajdowała się odpowiednia liczba sąsiadów. Liczba ta powinna się wahać w granicach od ${\displaystyle N_{N}/2}$ do ${\displaystyle 2N_{N}.}$

Tabela 1. Ilość sąsiadów w promieniu ${\displaystyle 2h}$ w zależności od liczby wymiarów symulacji

	Jeden wymiar	Dwa wymiary	Trzy wymiary
Liczba sąsiadów ${\displaystyle N_{N}}$	5	15	55

Jądro interpolacji

Funkcja jądra interpolacji powinna mieć postać

${\displaystyle W(r,h)={\frac {1}{h^{\nu }}}f(u),}$

gdzie:

${\displaystyle \nu }$ – liczba wymiarów,

${\displaystyle u=r/h,}$

${\displaystyle r}$ – odległość między cząstkami.

Dodatkowo funkcja ${\displaystyle f(u)}$ powinna spełniać warunki

${\displaystyle \int f(u)dV=1,}$

${\displaystyle \lim _{h\to 0}f(u)=\delta (r).}$

przy czym ${\displaystyle dV}$ jest elementem objętości, równym odpowiednio ${\displaystyle du,}$ ${\displaystyle 2\pi udu}$ lub ${\displaystyle 4\pi u^{2}du}$ w jednym, dwóch lub trzech wymiarach.

Najczęściej stosuje się jednak jądro interpolacyjne zaproponowane przez Monaghana

${\displaystyle W(r,h)={\frac {\sigma }{h^{\nu }}}\left\{{\begin{matrix}1-3/2u^{2}+3/4u^{3}&\quad 0\leqslant u\leqslant 1\\1/4(2-u)^{3}&\quad 1<u\leqslant 2\\0&\quad u>2\end{matrix}}\right.}$

A, jego gradient w postaci

${\displaystyle W(r,h)={\frac {\sigma }{h^{\nu +1}}}{\hat {r}}\left\{{\begin{matrix}-3u+9/4u^{2}&\quad 0\leqslant u\leqslant 1\\-3/4(2-u)^{2}&\quad 1<u\leqslant 2\\0&\quad u>2\end{matrix}}\right.}$

gdzie:

${\displaystyle {\hat {r}}={\frac {{\vec {r}}_{i}-{\vec {r}}_{j}}{|{\vec {r}}_{i}-{\vec {r}}_{j}|}}.}$

Tabela 2. Wartości parametrów ${\displaystyle \sigma }$ i ${\displaystyle \nu }$ w zależności od liczby wymiarów symulacji

	Jeden wymiar	Dwa wymiary	Trzy wymiary
${\displaystyle \nu }$	1	2	3
${\displaystyle \sigma }$	${\displaystyle {\frac {2}{3}}}$	${\displaystyle {\frac {10}{7\pi }}}$	${\displaystyle {\frac {1}{\pi }}}$

Sztuczna lepkość

Lepkość w metodzie SPH jest uwzględniana poprzez dodanie w równaniu Naviera-Stokesa przyspieszenia w postaci

${\displaystyle {\vec {a}}_{i}^{visc}=-\sum _{j=1}^{N}m_{j}\Pi _{ij}\nabla _{i}W_{ij},}$

gdzie:

${\displaystyle \Pi _{ij}={\frac {-\alpha c\mu _{ij}+\beta \mu _{ij}^{2}}{{\overline {\varrho }}_{ij}}},}$

${\displaystyle \alpha ,\beta }$ – stałe, często przyjmowane jako ${\displaystyle \alpha =0{,}5,}$ ${\displaystyle \beta =1{,}0,}$

${\displaystyle {\overline {\varrho }}_{ij}=(\varrho _{i}+\varrho _{j})/2,}$

${\displaystyle \mu _{ij}=\left\{{\begin{matrix}{\frac {h^{2}{\vec {v}}_{ij}\cdot {\vec {r}}_{ij}}{r_{ij}^{2}+\eta ^{2}h^{2}}}&{\vec {v}}_{ij}\cdot {\vec {r}}_{ij}<0\\0&{\vec {v}}_{ij}\cdot {\vec {r}}_{ij}\geqslant 0\end{matrix}}\right.,}$

${\displaystyle \eta }$ – stała, często przyjmowane jako ${\displaystyle \eta =0{,}1,}$

${\displaystyle {\vec {v}}_{ij}={\vec {v}}_{i}-{\vec {v}}_{j},}$

${\displaystyle {\vec {r}}_{ij}={\vec {r}}_{i}-{\vec {r}}_{j},}$

${\displaystyle r_{ij}=|{\vec {r}}_{i}-{\vec {r}}_{j}|.}$

Równanie stanu

Najczęściej stosuje się równanie stanu w postaci

${\displaystyle c={\sqrt {\frac {p\gamma }{\rho }}},}$

gdzie:

${\displaystyle \gamma =c_{p}/c_{V},}$

${\displaystyle c}$ – prędkość dźwięku w tym opisywanym płynie,

przy założeniu warunku

${\displaystyle v_{i}\ll c.}$

Innym równaniem stanu stosowanym dla cieczy jest

${\displaystyle p=p_{0}\left[\left({\frac {\rho }{\rho _{0}}}\right)^{\gamma }-1\right],}$

gdzie:

${\displaystyle \rho _{0}}$ i ${\displaystyle p_{0}}$ – maksymalne przyjęte wartości gęstości i ciśnienia,

${\displaystyle \gamma =7.}$

Ostateczne równanie

Po podstawieniu estymat do równania Naviera-Stokesa (i nie uwzględniając pola grawitacyjnego) otrzymujemy do rozwiązania równanie dla ${\displaystyle i}$-tej cząstki

${\displaystyle {\vec {a}}_{i}={\frac {d{\vec {v}}_{i}}{dt}}=-{\frac {c^{2}}{\gamma }}\sum _{j=1}^{N}m_{j}\left({\frac {1}{\rho _{i}}}+{\frac {1}{\rho _{j}}}+\Pi _{ij}\right)\nabla _{i}W_{ij},}$

które można wyznaczać dowolną metodą numerycznego całkowania równań różniczkowych, ale przy założeniu, że krok czasowy spełnia warunek Couranta

${\displaystyle \Delta t\leqslant 0{,}25{\frac {h}{c}}.}$

Przypisy

Bibliografia

• J.J. Monaghan, Notes Smoothed Particle Hydrodynamics. ipam.ucla.edu. [zarchiwizowane z tego adresu (2007-03-12)].
• B. Schlatter, A pedagogical tool using Smoothed Particle Hydrodynamics to model fluid flow past a system of cylinders. physics.orst.edu. [zarchiwizowane z tego adresu (2006-09-06)].
• J.P. Morris, Analysis of Smoothed Particle Hydrodynamics with Applications
• M. Ellero, Smoothed Particle Dynamics Methods for the Simulation of Viscoelastic Fluids. edocs.tu-berlin.de. [zarchiwizowane z tego adresu (2007-04-17)].