next up previous
Next: 8.1.3 Smoothed Particle Hydrodynamics Up: 8.1 Compressible Flow without Previous: 8.1.1 Explicit Eulerian Methods

8.1.2 Particle-in-Cell Method (PIC)

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 + \mbox{$\bf v$} \cdot \nabla c = 0$. Therefore
\begin{displaymath}
\frac{\partial }{\partial t} \left[ \rho c \right]+ \nabla \cdot
\left[ \rho c \mbox{$\bf v$} \right] = 0
\end{displaymath} (8.17)

(continuity equation for $c$).

$\Longrightarrow$Flow equations:
$\displaystyle \frac{\partial \rho }{\partial t}+ \nabla \cdot \left( \rho \mbox{$\bf v$} \right)$ $\textstyle =$ $\displaystyle 0$ (8.18)
$\displaystyle \frac{\partial \rho \mbox{$\bf v$}}{\partial t}+ \nabla \cdot
\left( \rho \mbox{$\bf v$} \mbox{$\bf v$} \right)$ $\textstyle =$ $\displaystyle - \nabla p$ (8.19)
$\displaystyle \frac{\partial }{\partial t} \left( \rho c \right) +
\nabla \cdot \left( \rho c \mbox{$\bf v$} \right)$ $\textstyle =$ $\displaystyle 0$ (8.20)

Strategy:
- Discretize the time: $t_{n+1} \equiv t_{n}+\Delta t$
- Assume a 2D (or 3D) Euler lattice with $\Delta x = \Delta y = \Delta l$
- Represent the (variable) local density by a number of particles in each cell
- Each particle $k$ represents a fluid element (not a molecule!) and carries a vector of properties,
\begin{displaymath}
\mbox{$\bf u$}_{k}^{n} \equiv \mbox{$\bf u$}_{k}(t_{n}) \equ...
...f v$}_{k}^{n}, c_{k}^{n} \right\}
\hspace{30pt} k=1, \dots N
\end{displaymath} (8.21)

The properties of the Eulerian (space fixed) cells $(i,j)$ are sums over the particles they contain:
$\displaystyle \rho_{i,j}^{n}$ $\textstyle =$ $\displaystyle \frac{m}{(\Delta l)^{2}} \sum_{k=1}^{N} \, \delta
\left[ \mbox{$\bf r$}_{k}^{n}(i,j) \right]$ (8.22)
$\displaystyle (\rho \mbox{$\bf v$})_{i,j}^{n}$ $\textstyle =$ $\displaystyle \frac{m}{(\Delta l)^{2}} \sum_{k=1}^{N} \,
\mbox{$\bf v$}_{k}^{n} \delta \left[ \mbox{$\bf r$}_{k}^{n}(i,j) \right]$ (8.23)
$\displaystyle (\rho c)_{i,j}^{n}$ $\textstyle =$ $\displaystyle \frac{m}{(\Delta l)^{2}} \sum_{k=1}^{N} \,
c_{k}^{n}\delta \left[ \mbox{$\bf r$}_{k}^{n}(i,j) \right]$ (8.24)

with
\begin{displaymath}
\delta \left[ \mbox{$\bf r$}(i,j) \right] \equiv
\delta \le...
...\right]
\delta \left[ \mbox{int}(\frac{y}{\Delta l})-j \right]
\end{displaymath} (8.25)

- Now rewrite 8.19 as
\begin{displaymath}
\rho \frac{\partial \mbox{$\bf v$}}{\partial t} = - \nabla p...
...abla \cdot
\left( \rho \mbox{$\bf v$} \mbox{$\bf v$} \right)
\end{displaymath} (8.26)

and first treat only the part $\rho \partial \mbox{$\bf v$} / \partial t = - \nabla p$:
$\displaystyle v_{x,i,j}^{n+1}$ $\textstyle =$ $\displaystyle v_{x,i,j}^{n}
- a \left( p_{i+1,j}^{n}- p_{i-1,j}^{n} \right)$ (8.27)
$\displaystyle v_{y,i,j}^{n+1}$ $\textstyle =$ $\displaystyle v_{y,i,j}^{n}
- a \left( p_{i,j+1}^{n}- p_{i,j-1}^{n} \right)$ (8.28)

with $a \equiv \Delta t \left/ 2 (\Delta l) \rho_{i,j}^{n} \right.$
- Update the particle properties $\mbox{$\bf v$}_{k}$ and $c_{k}$, thus:
\begin{displaymath}
\mbox{$\bf v$}_{k}^{n+1} = \mbox{$\bf v$}_{i,j}^{n+1} \quad ...
...and} \quad
c_{k}^{n+1}=p_{i,j}^{n}/(\rho_{i,j}^{n})^{\gamma}
\end{displaymath} (8.29)

- Now treat the Lagrangian part of equation 8.26, by letting the fluid particles move with appropriate velocities: Defining
\begin{displaymath}
\mbox{$\bf v$}_{i,j}^{n+1/2} = \frac{1}{2}
\left[ \mbox{$\bf v$}_{i,j}^{n+1} + \mbox{$\bf v$}_{i,j}^{n} \right]
\end{displaymath} (8.30)

compute the particle velocities as a weighted sum over the adjacent Eulerian cells:
\begin{displaymath}
\mbox{$\bf v$}_{k}^{n+1/2}= \frac{1}{(\Delta l)^{2}} \sum_{(ij)} a_{(ij)}
\mbox{$\bf v$}_{(ij)}^{n+1/2}
\end{displaymath} (8.31)

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 6.5.2; see Fig. 6.9.)
- Update the positions
\begin{displaymath}
\mbox{$\bf r$}_{k}^{n+1}=\mbox{$\bf r$}_{k}^{n} + \Delta t \,\mbox{$\bf v$}_{k}^{n+1/2}
\end{displaymath} (8.32)

to complete the time step.

\fbox{
\begin{minipage}{420pt}
{\bf PIC method (2-dimensional):}
At time $t_{n...
...{k}^{n+1}, c_{k}^{n+1} \right\}
\end{displaymath}\end{enumerate}\end{minipage}}
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.
next up previous
Next: 8.1.3 Smoothed Particle Hydrodynamics Up: 8.1 Compressible Flow without Previous: 8.1.1 Explicit Eulerian Methods
Franz J. Vesely Oct 2005
See also:
"Computational Physics - An Introduction," Kluwer-Plenum 2001