Consider an ideal gas; assume adiabatic equation of state (fast flow or
slow conduction of heat): $p/\rho^{\gamma}= c$ constant in a flowing
element $\longrightarrow$ Lagrangian time derivative
$\partial c / \partial t + \vec{v} \cdot \nabla c = 0$.
Therefore
where the weights $a_{(ij)}$ are the overlap areas of a square of
side length $\Delta l$ centered around particle $k\,$ and the nearest
Euler cells $(ij)$. (See the particle-mesh method of
Section 6.5.2;
see Fig. 6.9.)
Update the positions
$
\vec{r}_{k}^{n+1}=\vec{r}_{k}^{n} + \Delta t \,\vec{v}_{k}^{n+1/2}
$
(8.32)
to complete the time step.
Here is the procedure in concise form:
PIC method (2-dimensional):
At time $t_{n}$ the state of the fluid is represented by $N$ particles
with the property vectors $\vec{u}_{k}^{n} \equiv
\left\{ \vec{r}_{k}^{n},\vec{v}_{k}^{n}, c_{k}^{n} \right\} \;$
($k=1,\dots N$).
In each Eulerian cell of side length $\Delta l$ there should be at least
$\approx 100$ particles.
Compute, for each Euler cell $(i,j)$, the cell properties
with $ a \equiv \Delta t \left/ 2 (\Delta l) \rho_{i,j}^{n} \right.$.
For each fluid particle $k$ we now have
$\vec{v}_{k}^{n+1}=\vec{v}_{i,j}^{n+1}$
and $c_{k}^{n+1}=p_{i,j}^{n}/(\rho_{i,j}^{n})^{\gamma}$.
From the time-centered cell velocities
$\vec{v}_{i,j}^{n+1/2} \equiv$
$\left[ \vec{v}_{i,j}^{n+1} + \vec{v}_{i,j}^{n}
\left. \right] \right/ 2 \;$, compute for each particle $k$ an
intermediate velocity
Particle-in-cell method. Note that pressure gradients are evaluated
using an Eulerian grid, while the transport of mass, momentum and energy
is treated in continuous space.
8.1.3 Smoothed Particle Hydrodynamics (SPH)
Motivation:
PIC technique uses both Eulerian and Lagrangian elements.
Average density within a cell = number of point particles
in that cell.
Can we represent the local fluid density without a grid?
Lucy and
Gingold
&
Monaghan:
load each particle with a spatially extended
interpolation kernel $\longrightarrow$ Average local density
= sum over the individual contributions.
Let $w(\vec{r}-\vec{r_{i}})$ denote the interpolation kernel centered
around $\vec{r_{i}}$. Then the local density at $\vec{r}$ is
with a width $d$ chosen such that the number of particles within $d$
is $N \approx 5$ in 2 dimensions and $N \approx 15$ for 3 dimensions.
Now rewrite the Lagrangian equations of motion
8.8-8.11
in smoothed particle form.
Note:
In the momentum equation
$d \vec{v}/dt=- \nabla p / \rho$,
interpolating $\rho$ and $\nabla p$ directly would not conserve
linear and angular momentum
[Monaghan].
$\longrightarrow$ Use the identity
with $w_{ik} \equiv w(\vec{r}_{ik}) \equiv w(\vec{r}_{k}-\vec{r}_{i})$.
If $w_{ik}$ is Gaussian, this equation describes the motion of particle $i$
under the influence of central pair forces
Note: The density equation need not be integrated; just update all
positions $\vec{r}_{i}$, then invoke
equ.
8.33 to find $\rho(\vec{r}_{i})$.
To update $\vec{r}_{i}$ the obvious relation
$
\frac{\textstyle d \vec{r}_{i}}{\textstyle dt} = \vec{v}_{i}
$
(8.41)
with $\bar {\rho}_{ik} \equiv (\rho_{i}+\rho_{k})/2$. This relation
maintains angular and linear momentum conservation,
with the additional advantage that nearby particles will have similar
velocities
[Monaghan 89].
To solve equs.
8.39,
8.37,
8.41,
and
8.40,
use any suitable algorithm (see Chapter
4.)
Popular schemes are the leapfrog algorithm, predictor-corrector
and Runge-Kutta methods.
Example: Variant of the half-step technique
[Monaghan 89]:
Given all particle positions at time $t_{n}$, the local density
at $\vec{r}_{i}$ is computed from
8.39.
Writing equs.
8.37
and
8.40
as
$(i=1, \dots N)$.
In the simple case of an ideal gas undergoing adiabatic flow, the specific
energy $\epsilon$ may be replaced by $c \equiv p/\rho^{\gamma}$
$= \epsilon (\gamma -1)/\rho^{\gamma - 1} \,$.
A suitable interpolation kernel is assumed, e.g.
$ w(\vec{s}) = (1/\pi^{3/2}d^{3}) exp\{- s^{2}/d^{2}\}$,
with the width $d$ chosen so as to span about 5 (in 2 dimensions) or
15 (3-d) neighbors.
At each particle position $\vec{r}_{i}$ the density $\rho_{i}$
is computed by interpolation:
A modification of this scheme is obtained by including the density
$\rho_{i}$ in the state vector of particle $i$ and integrating
the pertinent equation of motion,
8.39. The time step
integrations for $\vec{r}$, $\vec{v}$, $\rho$ and $\epsilon$ may then be
performed simultaneously, and the evaluation of the density according to
the SPH averaging prescription is omitted. This procedure works faster,
but exact mass conservation is not guaranteed.
Note: SPH is not restricted to compressible inviscid flow.
Incompressibility may be introduced by using an equation of state
that keeps compressibility below a few percent
[Monaghan],
and viscosity is added by an additional term in the
equations of motions for momentum and energy, equs.
8.37
and
8.40:
and $c_{ik} \equiv c_{k}- c_{i}$ and
$\bar{c}_{ik} \equiv (c_{i}+c_{k})/2$ etc.
This form of $\Pi$ introduces the effects of shear and bulk viscosity.
The parameters $\alpha$ and $\beta$ should be near
$\alpha = 1$ and $\beta = 2$ for best results
[Monaghan].
The quantity $\eta$ prevents singularities for $r_{ik} \approx 0$. It should be
chosen such that $\eta^{2}=0.01 d^{2}$.
Additional features:
Thermal conduction may be included. See
[Monaghan 89].
Interfaces: Introduce dummy particles on the far side of
the boundary. By picking the properties of these particles appropriately
one can mimick a free surface or a "sticky" solid boundary.
See Nugent and Posch
[Nugent]
for free surfaces, and
[Ivanov]:
for rough interfaces.
Sample application:
"Rayleigh-Bénard" convection
A fluid layer is heated carefully from below and cooled from above.
$\longrightarrow$ Formation of stable convective rolls transporting heat from
the bottom to the top. See the Figure for a match
between SPH and an Euler-type calculation
[Hoover 99]:
- Computing times comparable for both calculations
- Results are in good agreement
- Fluctuations in SPH (like in any particle-type calculation), none
in Euler
- SPH code is quite simple -- similar to an MD program; Euler code
very massive
Figure: Comparison of Smoothed Particle Hydrodynamics with an
Eulerian finite-difference calculation. The density (above) and
temperature (below) contours for a stationary Rayleigh-Bénard
flow are shown. Left: SPH; right: Euler.
(From
[Hoover 99],
with kind permission by the author)