Time Step Control Stability

The stability conditions of explicit scheme in SPH formulation can be written over cells or on nodes.

Cell Time Step

In case of cell stability computation (when no nodal time step is used), the stable time step is computed as:(1)

Δ t = Δ t_{s c a} \cdot \min_{i} (\frac{d_{i}}{c_{i} (α_{i} + \sqrt{α_{i}^{2} + 1})}), w i t h α_{i} = (q_{b} + \frac{q_{a} \cdot {\bar{μ}}_{i} \cdot d_{i}}{c_{i}}), a n d {\bar{μ}}_{i} = \max_{j} (μ_{i j})

$Δ t_{s c a}$ is the user-defined coefficient (Radioss option /DT or /DT/SPHCEL). The value of $Δ T_{s c a}$ =0.3 is recommended. ¹

In case of nodal time step, stability time step is computed in a more robust way:(2)

$Δ t_{i} = \sqrt{\frac{2 m_{i}}{K_{i}}}$ at particle $i$

Use the following notations, if kernel correction:(3)

W_{j} (i) = \hat{W} (x_{i} - x_{j^{'}} \frac{d_{i} + d_{j}}{2}) a n d \nabla W_{j} (i) = g r a d |_{x i} [\hat{W} (x - x_{j^{'}} \frac{d_{i} + d_{j}}{2})]

Or, if no kernel correction:(4)

W_{j} (i) = W (x_{i} - x_{j^{'}} \frac{d_{i} + d_{j}}{2}) a n d \nabla W_{j} (i) = g r a d |_{x i} [W (x - x_{j^{'}} \frac{d_{i} + d_{j}}{2})]

Recalling that apart from the artificial viscosity terms:(5)

F_{i} = \sum_{j} F_{i j}, F_{i j} = V_{i} V_{j} [p_{i} \nabla W_{j} (i) - p_{j} \nabla W_{j} (j)]

write (6)

| K_{i j} | = ‖ \frac{d F_{i j}}{d (u_{i} - u_{j})} ‖ \leq \frac{d}{d (u_{i} - u_{j})} (V_{i} V_{j} [p_{i} ‖ \nabla W_{j} (i) ‖ + p_{j} ‖ \nabla W_{i} (j) ‖])

Where,

u_{i} - u_{j}

is the relative displacement of particles

i

and

j

. Keeping the only first order terms leads to:(7)

| K_{i j} | \leq V_{i} V_{j} [\frac{d p_{i}}{d (u_{i} - u_{j})} ‖ \nabla W_{j} (i) ‖ + \frac{d p_{j}}{d (u_{i} - u_{j})} ‖ \nabla W_{i} (j) ‖]

Where, (8)

V_{i} V_{j} \frac{d p_{i}}{d (u_{i} - u_{j})} ‖ \nabla W_{j} (i) ‖ = V_{i} V_{j} \frac{d p_{i}}{d ρ_{i}} \cdot \frac{d ρ_{i}}{d (u_{i} - u_{j})} ‖ \nabla W_{j} (i) ‖ = V_{i} V_{j} c_{i}^{2} \frac{d ρ_{i}}{d (u_{i} - u_{j})} ‖ \nabla W_{j} (i) ‖

that is(9)

V_{i} V_{j} \frac{d p_{i}}{d (u_{i} - u_{j})} ‖ \nabla W_{j} (i) ‖ = m_{i} c_{i}^{2} {\dot{V}}_{j}^{2} {‖ \nabla W_{j} (i) ‖}^{2}

Same reasoning leads to:(10)

V_{i} V_{j} \frac{d p_{j}}{d (u_{i} - u_{j})} ‖ \nabla W_{i} (j) ‖ = m_{j} c_{j}^{2} {\dot{V}}_{i}^{2} {‖ \nabla W_{i} (j) ‖}^{2}

So that (11)

| K_{i j} | \leq m_{i} c_{i}^{2} {\dot{V}}_{j}^{2} {‖ \nabla W_{j} (i) ‖}^{2} + m_{j} c_{j}^{2} {\dot{V}}_{i}^{2} {‖ \nabla W_{i} (j) ‖}^{2}

Stiffness around node

i

is then estimated as:(12)

| K_{i} | \leq \sum_{j} | K_{i j} |

¹ Monaghan J.J., “Smoothed Particle Hydrodynamics”, Annu.Rev.Astron.Astro-phys; Vol. 30; pp. 543-574, 1992.