Smoothed Particle Hydrodynamics

Smoothed-particle hydrodynamics ( SPH ; German: smoothed particle hydrodynamics ) is a numerical method for solving hydrodynamic equations . It is used in astrophysics , ballistics and tsunami calculations , among other things . SPH is a Lagrange method, i.e. H. the coordinates used move with the fluid. SPH is a particularly easy to implement and robust method.

method

General

In Smoothed Particle Hydrodynamics, the fluid to be simulated is divided into elements. Similar to the Monte Carlo methods , the elements are randomly distributed over the liquid. This minimizes the error to be expected. The average distance between these elements is represented by the Smoothing Length . It is the most important parameter of the method. The kernel smooths the fluid between the particles, hence the name. Each size (e.g. density ) is calculated by summing over all particles. Each individual particle receives a share of this size in the form of a scalar. This turns the partial differential equations of hydrodynamics into ordinary differential equations, which greatly simplifies the calculations. SPH is a very empirical method. This means that many things are done because they work, not because there is a strict mathematical derivation. ${\displaystyle h}$ ${\displaystyle \rho }$

Derivation

The formal derivation runs either via a Lagrange function or via an integral interpolation. The integral interpolation for a quantity is based on an identity, where the Dirac delta distribution denotes: ${\displaystyle A_{1}}$ ${\displaystyle \delta ({\vec {r}})}$

${\displaystyle 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}}\,'}$

Then the distribution is approximated by a core , where is the smoothing length. So that the approximation remains valid in the limit case, normalization and identity with the distribution in the limit value for h  → 0 can be required: ${\displaystyle \delta }$ ${\displaystyle W({\vec {r}}-{\vec {r}}\,',h)}$ ${\displaystyle h}$ ${\displaystyle \delta }$

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

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

In fact, this is no longer the case with most of the cores used. In order to obtain the division into mass elements, one expands with the density and leaves it greater than 0. For the case of infinitely many, infinitely small particles, the sum goes over into the integral. Numerically, one will always have to be satisfied with a finite number of particles: ${\displaystyle \rho }$ ${\displaystyle h}$

${\displaystyle 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}}$

It is the mass of the particle b and the density at the location of the particle b: ${\displaystyle m_{b}}$ ${\displaystyle \rho _{b}}$

${\displaystyle \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)}$

With that we have derived the basic equation of Smoothed Particle Hydrodynamics (right part). The quantity A is calculated as a sum over all particles. You can see that the variable that depends on r has become a scalar multiplied by the kernel. This leads to a strong simplification of differential equations, since a derivative no longer affects the size, but only the kernel: ${\displaystyle A_{S}}$ ${\displaystyle A_{b}}$

${\displaystyle \nabla A({\vec {r}}\,)=\sum \limits _{b}m_{b}{\frac {A_{b}}{\rho _{b}}}\nabla W({\vec {r}}-{\vec {r}}\,',h)}$

Core and Smoothing Length

Smoothing length

Probably the most important parameter of the SPH is the smoothing length . It defines the resolution of the method and thus has a strong influence on the accuracy and computational effort in simulations. If the kernel is selected accordingly (see below), it also defines the number of neighbors to be included in the calculation. Up to a few tens of particles per size are common. For good results, use the mean density of the fluid as a guide: ${\displaystyle h}$

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

with particles, dimensions and ${\displaystyle n}$ ${\displaystyle \nu }$

${\displaystyle \langle \rho \rangle ={\frac {1}{n}}\sum \limits _{b}\rho _{b}}$

In modern codes, you choose time-dependent. With ${\displaystyle h=h(t)}$

${\displaystyle {\frac {\mathrm {d} h_{a}}{\mathrm {d} t}}=-\left({\frac {h_{a}}{\nu \rho _{a}}}\right){\frac {\mathrm {d} \rho _{a}}{\mathrm {d} t}}}$

a higher resolution is then used in areas of high densities, while the smoothing length increases in areas of low densities. This allows the computational effort to be reduced while maintaining the same level of accuracy.

core

The core is probably the most important structure of the SPH method. Different kernels correspond to different difference schemes in grid methods. To interpret SPH equations, it is advantageous to use a kernel in the form of a Gaussian curve:

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

In numerical terms, however, this approach is not very suitable, since in this case one often values ​​clear behavior with regard to the range of the core. I.e. one chooses a kernel that is zero above a certain level in order to be able to clearly determine the number of neighbors that are included in the calculation. This allows you to limit the computational effort required. As mentioned earlier, SPH is a very empirical method; H. very different cores are required for different applications. The exact choice is a matter of experience and is often made on the basis of trial and error . Since a core is often implemented in its own function, the effort to replace or change it is minimal. Cores based on splines are often used: ${\displaystyle {\vec {r}}}$

${\displaystyle 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}}}$

With , a normalization constant and the number of dimensions . Here only particles up to the next but one neighbor are included in the calculation. In addition, the 2nd derivative of this nucleus is not constant, which is why it does not depend on the disorder of the particles. ${\displaystyle q={\vec {r}}-{\vec {r}}\,'}$ ${\displaystyle \sigma }$ ${\displaystyle \nu }$

Error estimates

Two approximations were made for the derivation using integral interpolation functions. First, it was assumed that the summation only takes place over a finite number of particles. ${\displaystyle h>0}$

• For identity, i.e. H. with and any number of particles, a Taylor expansion gives an error of . ${\displaystyle h=0}$ ${\displaystyle O(h^{2})}$
• An error can also be calculated for the summation approximation with the help of the Shoenberg formula if the particles are distributed in an orderly manner in the fluid. ${\displaystyle O(h^{2})}$
• In the case of disordered particles, there is no traditional error estimation.

With simulations with SPH one is therefore always dependent on the comparison with other simulations, at least for an error assessment. Some publications mention that the errors are usually well below those of a Monte Carlo simulation ; this is also a matter of experience. In general, SPH tends to smear out discontinuities, so it is locally rather imprecise, especially in the case of simulations with few particles. The behavior becomes significantly better for large numbers of particles. However, the global behavior is very good even with low particle numbers, which corresponds to low computational effort. In other words, global variables such as energy are well represented. With SPH, a globally good simulation can often be programmed with little effort, which can be calculated on workstations in an acceptable time.

advantages and disadvantages

Advantages:

• SPH is a Lagrange method; the continuity equation is automatically satisfied.
• The code is very robust; H. almost always gives results
• The implementation of SPH is comparatively easy, as is the testing of different kernels.
• With the help of a Gaussian function as a kernel, theoretical results can be easily interpreted.
• In modern code there is a dependency of the computational effort on the number of particles. ${\displaystyle N\log N}$
• SPH shows good global results with low particle numbers.

Disadvantage:

• The code is often too robust; despite an incorrect model, SPH can deliver results that are then physically incorrect
• The error estimation is often problematic and can only be obtained in comparison with the results of other methods
• The method is highly dispersive
• High numbers of particles are required for good accuracy. The advantage of low computational effort is thus not applicable
• The treatment of discontinuities is often difficult as structures are smoothed on scales smaller than the smoothing length.

Hydrodynamic equations in SPH

Symmetrization

To formulate the hydrodynamics in SPH, the apparently simplest approach is the basic equation in the hydrodynamic equations such as B. to use the Navier-Stokes equation . However, the resulting equations are not symmetrical with respect to particle exchange. Therefore, in this case, many conservation laws for energy, angular momentum etc. no longer apply. However, it is often possible to save this by writing the density into the respective differential operator and using the product rule:

${\displaystyle \rho \nabla A=\nabla (\rho A)-A\nabla \rho }$

In this way, symmetrical equations can often be derived. All of this is not strictly formal, but only because it gives better results.

Movement of the fluid

The simplest way is to use the definition of speed:

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

The movement of one particle is not linked to that of the other, which can often lead to problems. That is why the XSPH method ("Extended SPH") was developed:

${\displaystyle {\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}}$

with an averaged density:

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

and a coupling parameter ε. The order of the particles is thus better preserved without the need to introduce additional viscosity.

Continuity equation in SPH

If we put the density in the basic equation, we get

${\displaystyle \rho _{a}=\sum \limits _{b}m_{b}W_{ab}}$

for a particle a. From this the SPH continuity equation can be calculated

${\displaystyle {\frac {\mathrm {d} \rho _{a}}{\mathrm {d} t}}=-\sum \limits _{b}m_{b}v_{ab}\nabla _{a}W_{ab}}$

Euler's equation in SPH

For the Euler equation we get:

${\displaystyle {\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}}$

This equation is not symmetrical about particle exchange: momentum and torque are not preserved. Therefore we use the trick indicated above for the pressure gradient:

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

From what we get the desired symmetric equation:

${\displaystyle {\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}}$

If we use a Gaussian function , the result is a central force that acts equally on both particles:

${\displaystyle 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}}$

viscosity

Like almost every numerical method, SPH also generates viscosity through calculation inaccuracies . However, this is often not sufficient for modeling. Therefore, similar to the transition from the Euler equation to the Navier-Stokes equation , a viscosity tensor is introduced. The exact choice of this tensor depends heavily on the model.

${\displaystyle {\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}}$

Applications

SPH is used in many different fields such as astrophysics . There are also relativistic and magnetic SPH methods:

• Gas dynamics
• Galaxy formation and merging
• Binary star systems , accretion disks and star collisions
• Moon formation
• Relativistic problems
• Magnetic problems

Further publications

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

Web links

Movies

• Visualizations of galaxies: evolution calculations (English)
• Mathias Steinmetz: Galaxies: Formation and Mergers. Leibniz Institute for Astrophysics Potsdam (English). 

code

• Website of the gadget code whose Gasdynamic uses SPH (English)
• SimPARTIX is a software package for SPH and DEM simulations from Fraunhofer IWM
• Pasimodo simulates various particle methods , including a. SPH
• ThreeParticle / CAE is a multiphysics software platform for DEM and SPH simulation from BECKER 3D GmbH
• SOFA - Simulation Open Framework Architecture ( Memento from April 8, 2008 in the Internet Archive ) simulates SPH with regard to real-time applications.
• OpenTissue ( Memento of August 24, 2011 in the Internet Archive ) also offers SPH with a focus on real-time application.