SPH.tex

\section{The SPH method}
\label{sec:sph}
Smoothed particle hydrodynamics starts with representing a function at a point
as an integral by taking a convolution with the Dirac delta function,
%
\begin{equation}
  \label{eq:func:approx}
  f(\ten{x}) = \int_{V} f(\ten{y}) \delta(\ten{x} - \ten{y}) d\ten{y}
\end{equation}
%
where $f$ is a function which takes a vector $\ten{x}$ as an input,
$\delta(\ten{x} - \ten{y})$ is the Dirac delta function, and the integral is
over the volume $V$, which contains $\ten{x}$.  The delta function is then
approximated with a weighting or ``smoothed'' function called a smoothing
kernel, $W(\ten{x} - \ten{y}, h)$, where $h$ is the smoothing length which
determines the fall of the smoothing kernel. The smoothing kernel should
satisfy few properties to accurately and effectively approximate the function
namely, $\int_{V} W(\ten{x} - \ten{y}, h) d\ten{y} = 1$ (i.e., partition of
unity), it should approximate the Dirac delta function in the limit $h
\rightarrow 0$, smoothing kernel should be symmetric with respect to $(\ten{x}
- \ten{y})$, compact to be computationally efficient. Now the above
representation can be approximated as,
%
\begin{equation}
  \label{eq:func:smoothed:approx}
  <f(\ten{x})> = \int_{V} f(\ten{y}) W(\ten{x} - \ten{y}, h) d\ten{y}
\end{equation}
%
one of the kernel which is widely used in SPH is the Cubic spline kernel which
belongs to the B-spline family is given by,
%
\begin{align}
  \label{sph:kernel:cs}
  W(\ten{x} - \ten{y}, h) =
  \begin{cases}
     \sigma [{(2-q)}^3 - 4{(1-q)}^3] &\text{if}~0 \le q \le 1, \\
     \sigma {(2-q)}^3 &\text{if}~1 < q \le 2, \\
     0 &\text{if}~q > 2, \\
  \end{cases}
\end{align}
%
\todoin{references for kernels needed.}
where $q = \frac{|\ten{x} - \ten{y}|}{h}$, $\sigma = \frac{1}{6h}$ in 1-D,
$\frac{15}{14\pi h^2}$ in 2-D, and $\frac{1}{4\pi h^3}$ in 3-D. There are many
other kernels like Quintic spline, Quartic spline, Wendland quintic, Gaussian,
Super Gaussian which are supported by PySPH. Now the function is discretized
and the integral in the function approximation is discretized by a summation,
where the summation is over all the points that falls under the support of
the kernel, $\ten{y} \in
N(\ten{x})$,
\todoin{Write about all the variables introduced N, \ten{y}}
\begin{equation}
  \label{eq:sph:approx}
  <f(\ten{x})>_{sph} = \sum_{\ten{y} \in N(\ten{x})} \Delta V f(\ten{y})
  W(|\ten{x} - \ten{y}|, h)
\end{equation}
where $\Delta V = m(\ten{y})/\rho(\ten{y})$ is the volume at $\ten{y}$.

Now we use the SPH function approximation to discretize the governing
equations. We consider here the incompressible fluid, but in general SPH
approximation can be used for any partial differential equation. The equations
which govern the incompressible fluid motion in Lagrangian formulation are, the
continuity equation,
\begin{equation}
  \label{eq:ns:mass}
   \frac{d\rho}{dt} = - \rho \nabla \cdot \ten{u}
\end{equation}
where $\ten{u}$ is the velocity of the fluid. The momentum equation is given by,
%
\begin{equation}
  \label{eq:ns:mom}
  \frac{d\ten{u}}{dt} = - \frac{1}{\rho}\nabla p + \nu \nabla \cdot
  (\nabla \ten{u}) + \ten{f}
\end{equation}
%
where $\frac{d(.)}{dt}$ is the material derivative, $\rho$ is the density,
$p$ is the pressure, $\nu$ is the dynamic viscosity, $\ten{f}$ is
the external body force.

We simulate the incompressible fluid flow using Weakly compressible
\todoin{It should be clear on what equations we are solving}
formulation, where a stiff equation of state describes the relation between
pressure and density, and completes the above equations. The equation of state
we use is,
\todoin{Should this be called Tait equation? See Monaghan 2012 annual review. If
not we need to change the code listings too!}
%
\begin{equation}
  \label{eq:eos}
  p = \frac{c^2_0 \rho_0}{\gamma}
  \left[\left(\frac{\rho}{\rho_0}\right)^{\Gamma} - 1\right]
\end{equation}
%
where $\Gamma = 7$.
We now use the SPH approximation to bring the governing equation to
a discretized form.

The density is computed by,
\begin{equation}
  \label{eq:sph:summation_density}
  \rho_i = \sum_{j \in N(x_i)} m_j W(\ten{x}_i - \ten{x}_j, h)
\end{equation}
where $\rho_i = \rho(\ten{x}_i)$, is the density of the fluid particle, $m_j =
m(\ten{x}_j)$ is the mass of the particle in the neighbourhood of $\ten{x}_i$,
$h = hdx \times \Delta x^{d}$ is the kernel radius, where $hdx$ is a constant
\todoin{Remove the d}
in the range of $1 - 2$, $\Delta x$ is the initial particle spacing, and $d$ is
the spatial dimension. Alternatively one could also use the continuity
equation~\ref{eq:ns:mass} to compute density,
%
\begin{equation}
  \label{eqn:continuity:sph}
  \frac{d\rho_i}{dt} = \sum_{j \in N(x_i)} m_j (\ten{v}_i - \ten{v}_j) \cdot
  \nabla W_{ij}
\end{equation}
%
The momentum equation in the Lagrangian form is given by,
\begin{equation}
  \label{eq:mom:sph}
  \frac{d\ten{u}_i}{dt} = \ten{f}_{\text{p}} + \ten{f}_{\text{visc}}
  + \ten{f}_{\text{body}}
\end{equation}
where $\ten{u}_i$ is the velocity of the $i^{\text{th}}$ particle,
$\ten{f}_{\text{p}}$ is the acceleration due to pressure gradient,
$\ten{f}_{\text{visc}}$ is the acceleration due to viscous forces, and
$\ten{f}_{\text{body}}$ is the acceleration due to external body forces, which
are discretized as,
%
\begin{equation}
  \label{eq:mom:sph:press}
  \ten{f}_{\text{p}} = - \sum_{j \in N(i)} m_j
  \left( \frac{p_i}{\rho_i^2} + \frac{p_j}{\rho_j^2} \right) \nabla_i W_{ij}
\end{equation}
%
where $p_i$ is pressure of $i^{\text{th}}$ particle, $p_j$ is pressure, $m_j$
is the mass, and $\rho_j$ is the density of the $j^{\text{th}}$ particle,
$\nabla_i W_{ij} = \frac{\partial W}{\partial \ten{x}_i}(\ten{x}_i - \ten{x}_j,
h)$ is the gradient of the smoothing kernel.
%
\begin{equation}
  \label{eq:mom:sph:visc}
  \ten{f}_{\text{visc}} = \sum_{j \in N(i)} m_j \frac{4 \nu}{(\rho_i + \rho_j)}
  \frac{\ten{r}_{ij} \cdot \nabla W_{ij}}{(|\ten{r}_{ij}|^{2} + \eta
      h^{2})} \ten{u}_{ij}
\end{equation}
%
where $\ten{r}_{ij} = \ten{x}_i - \ten{x}_j$, finally the equation of motion for the
particles is,
%
\begin{equation}
  \label{eq:eom}
  \frac{d \ten{x}_i}{dt} = \ten{u}_i.
\end{equation}
%
The problem we are going to consider in the next section is the simulation of
2-D Dam break~\cite{wcsph-state-of-the-art-2010} where a column of fluid is
placed in a tank and simulated under the influence of gravity. The boundary
conditions for the purpose of simplicity in explaining the framework are
implemented in a naive inaccurate way where the tank is represented by solid
particles on which we solve the continuity equation with the influence of only
fluid particles and then update the pressure on the solid by solving the
equation of state.

\subsection{Time Integration}
The equations of motion and momentum equations are then integrated in time to
obtain the positions and velocities respectively, A Predictor-Corrector
integrator~\cite{monaghan-review:2005} is used to integrate in time which is
described as,
%
\begin{equation}
  \ten{u}^{t+\frac{1}{2}}_i = \ten{u}^{t}_i + \frac{\Delta t}{2}{\left(\frac{d
              \ten{u}_i}{dt}\right)}^t
\end{equation}
%
\begin{equation}
  \ten{x}^{t+1}_i = \ten{x}^{t}_i + \Delta t\ten{u}^{t + \frac{1}{2}}_i
\end{equation}
%
\begin{equation}
  \ten{u}^{t+1}_i = \ten{u}^{t+\frac{1}{2}}_i + \frac{\Delta
    t}{2}{\left(\frac{d \ten{u}_i}{dt}\right)}^{t + 1}
\end{equation}
%
\subsection{Schemes}
Although WCSPH scheme works for many problems it has its own disadvantages, to
name a few, stiff equation of state, tensile instability i.e., particles start
to pair together due to negative pressures. Many new schemes are proposed to
improve the simulations in different areas using wide variety of methods.  PySPH
has many number of the schemes tested and implemented into the framework here we
will describe some of them.

The standard WCSPH as described above, WCSPH with hg
correction~\cite{hughes-graham:compare-wcsph:jhr:2010} which imposes a minimum
value of density on the solid particles to make sure the fluid particles are not
sticking to the boundary, WCSPH with tensile instability
correction~\cite{sph:tensile-instab:monaghan:jcp2000},
$\delta$-SPH~\cite{marrone-deltasph:cmame:2011} a variant of WCSPH which adds
artificial diffusive terms to the continuity equation to reduce undesired
high-frequency oscillations in the pressure field.

Transport velocity formulation (TVF)~\cite{Adami2013} which ensures uniform
particle distribution by adding a constant background pressure to the convective
velocity, Generalized TVF (GTVF)~\cite{zhang_hu_adams17} a variant of TVF with
variable background pressure and supports simulation of solid mechanics and
captures of free surface flows accurately.

Entropically Damped Artificial compressibility (EDAC)~\cite{edac-sph:cf:2019}
where a pressure evolution equation is solved instead of the equation of state
resulting in accurate pressure distribution.  Dual time SPH
(DTSPH)~\cite{pr:dtsph:2019} a variant of WCSPH which solves for pressure using
pseudo time iterations.

Incompressible SPH (ISPH)~\cite{sph:psph:cummins-rudman:jcp:1999} which solves
the coupled incompressible fluid equations by solving a pressure Poisson
equation leading to higher timesteps and accurate pressures, uniform particle
distribution is maintained by shifting
algorithms~\cite{acc_stab_xu:jcp:2009,fickian_smoothing_sph:skillen:cmame:2013},
a simple itrative ISPH (SISPH)~\cite{muta2019simple} is a variant of ISPH which
solves the pressure Poisson equation using matrix free iterative approach and
ensure particle uniformity using GTVF. Implicit ISPH
(IISPH)~\cite{iisph:ihmsen:tvcg-2014} widely used in the graphics community to
simulate fluid problems with larger timesteps.

Conservative reproducible kernel SPH~\cite{crksph:jcp:2017} which
uses a first-order consistent reproducing kernel to capture linear variation
exactly meanwhile ensuring conservation properties in reproducing kernels (which
are not spatially symmetric).  Kernel and gradient
corrections~\cite{bonet_lok:cmame:1999} which ensure exact interpolation of a
linear field. Density corrections using Shepard filters which ensures partition
of unity of smoothing kernels.

Godunov SPH (GSPH)~\cite{inutsuka2002} is artificial-viscosity free Riemann
solver with a consistent density estimate, Approximate Riemann solver with
Godunov SPH (AGSPH)~\cite{puri2014}. Adaptive density kernel estimate
(ADKE)~\cite{sigalotti2006} which is a variant of WCSPH with an additional
heating term to the energy equation and also uses variable smoothing lengths,
$h$.

\todoin{cite these and explain.}
- Parshikov
- Zhang Hu Adams (Riemann Solver based)

Schemes with governing equations solving solid mechanics are also a part of the
framework like, Gray-Monaghan~\cite{gray2001}. Equations modelling surface
tension~\cite{morris:surface:tension:2000, adami:surface:tension:2010} are also
included in the framework.

\todoin{Add references to DEM, FSI, Akinci, Rigid body if they are implemented
  in PySPH and write a short explanation.}

\todoin{Citation needed for IO bcs, and a short explanation if possible.}
Boundary conditions are also implemented in the framework which include, wall
boundary conditions~\cite{Adami2012}, inlet boundary conditions and various outlet
boundary conditions which includes, Do nothing, Hybrid, Method of
characteristics, Mirror~\cite{tafuni2018}

%%% Local Variables:
%%% mode: latex
%%% TeX-master: "paper"
%%% End: