SPH Integration of Continuum Equations

In order to keep an almost constant number of neighbors contributing at each particle, use smoothing length varying in time and in space.

Consider $d_{i}$ the smoothing length related to particle $i$ ;

(1)

$W_{j} (i) = \hat{W} (x_{i} - x_{j}, \frac{d_{i} + d_{j}}{2})$ and $\nabla W_{j} (i) = g r a d |_{x_{i}} [\hat{W} (x - x_{j}, \frac{d_{i} + d_{j}}{2})]$ if kernel correction

or (2)

$W_{j} (i) = W (x_{i} - x_{j}, \frac{d_{i} + d_{j}}{2})$ and $\nabla W_{j} (i) = g r a d |_{x_{i}} [W (x - x_{j}, \frac{d_{i} + d_{j}}{2})]$ without kernel correction

At each time step, density is updated for each particle

i

, according to:(3)

{\frac{d ρ}{d t} |}_{i} = - ρ_{i} \nabla \cdot {v |}_{i}

with (4)

\nabla \cdot {v |}_{i} = \sum \frac{m_{j}}{ρ_{j}} (v_{i} - v_{j}) \cdot \nabla W_{j} (i)

Where,

$m_{j}$: Mass of a particle $i$
$ρ_{i}$: Density
$v_{j}$: Velocity

Strain tensor is obtained by the same way when non pure hydrodynamic laws are used or in the other words when law uses deviatoric terms of the strain tensor:(5)

{\frac{d v^{α}}{d x^{β}} |}_{i} = \sum \frac{m_{j}}{ρ_{j}} (v_{i}^{α} - v_{j}^{α}) \frac{d W_{j}}{d x^{β}} (i), α = 1...3, β = 1...3.

Next the constitutive law is integrated for each particle. Then Forces are computed according to:(6)

{m_{i} \frac{d v}{d t} |}_{i} = - \sum_{j} V_{i} V_{j} [p_{i} \nabla W_{j} (i) - p_{j} \nabla W_{i} (j)] - \sum_{j} m_{i} m_{j} π_{i j} \frac{[\nabla W_{j} (i) - \nabla W_{i} (j)]}{2}

Where

ρ_{i}

and

p_{j}

are pressures at particles

i

and

j

, and

π_{i j}

is a term for artificial viscosity. The expression is more complex for non pure hydrodynamic laws.

Note: The previous equation reduces to the following one when there is no kernel correction:

(7)

{m_{i} \frac{d v}{d t} |}_{i} = - \sum_{j} V_{i} V_{j} [p_{i} + p_{j}] \nabla W_{j} (i) - \sum_{j} m_{i} m_{j} π_{i j} \nabla W_{j} (i),

since $\nabla W_{i} (j) = - \nabla W_{j} (i)$

Then, in order particles to keep almost a constant number of neighbors into their kernels (

ρ d^{3}

is kept constant), search distances are updated according to:(8)

\frac{d (d_{i})}{d t} = {d_{i} \frac{\nabla \cdot {v |}_{i}}{3}}_{}