SPH Multiphase Flow Solver

The SPH flow solver solves the integral conservation equations of mass and momentum in a sequential manner. The solver uses a pressure-correction method for the pressure prediction.

The pressure-correction equation is constructed from the continuity equation and the momentum equations such that a predicted velocity field is sought that fulfills the continuity equation, which is achieved by correcting the pressure. This method is also called a predictor-corrector approach. Pressure as a variable is obtained from the pressure-correction equation.

Prediction Step

At the prediction step, an intermediate particle position

r_{i}^{*}

is computed from:

图 1. EQUATION_DISPLAY

\frac{r_{i}^{*} - r_{i}^{n}}{Δ t} = v_{i}^{n}

(1117)

where

Δ t

is the time-step,

r_{i}^{n}

is the particle position at the current time-step, and

v_{i}^{n}

is the velocity at the current time-step.

For the explicit time, the intermediate velocity $v_{i}^{*}$ is computed from:

图 2. EQUATION_DISPLAY

m_{i} \frac{v_{i}^{*} - v_{i}^{n}}{Δ t} + V_{i} \nabla (p^{n}) - V_{i} μ \nabla^{2} (v^{n}) = V_{i} f_{b, i}^{n}

(1118)

And for the implicit time integration:

图 3. EQUATION_DISPLAY

m_{i} \frac{v_{i}^{*} - v_{i}^{n}}{Δ t} + V_{i} \nabla (p^{n}) - V_{i} μ \nabla^{2} (v^{*}) = V_{i} f_{b, i}^{*}

(1119)

The linear system can be solved using the Biconjugate Gradient Stabilized (BICGSTAB) method, the Conjugate Gradient (CG) method, or the Generalized Minimal Residual (GMRES) method. In terms of efficiency, BICGSTAB is generally the preferred method.

Free Surface Stabilization: The solver provides free surface stabilization for particles near the free surface. This stabilization introduces an additional dissipation term to the moment equation. A user-specified free surface stabilization factor can be applied to modify the magnitude of the dissipation term. Small values of the free surface stabilization factor can result in a highly dispersed free surface. Therefore, it is essential to identify a stabilization factor value that maintains the stability of the free surface.

Projection Step

To compute the pressure correction

p^{'} = p^{n + 1} - p^{n}

, the pressure-correction equation is solved:

图 4. EQUATION_DISPLAY

\nabla \cdot (\nabla (p^{'})) = \frac{ρ_{0}}{Δ t} \nabla \cdot (v^{*})

(1120)

The linear system can be solved using the Biconjugate Gradient Stabilized (BICGSTAB) method, the Jacobi method, or the Generalized Minimal Residual (GMRES) method. In terms of efficiency, BICGSTAB is generally the preferred method.

The velocity $v^{n + 1}$ is obtained from the corrected intermediate velocity:

图 5. EQUATION_DISPLAY

\frac{{v_{i}^{n + 1} - v}_{i}^{*}}{Δ t} = - \frac{1}{ρ_{o}} \nabla (p^{'})

(1121)

The final position

r^{n + 1}

is then updated from:

图 6. EQUATION_DISPLAY

\frac{r_{i}^{n + 1} - r_{i}^{n}}{Δ t} = \frac{v_{i}^{n} + v_{i}^{n + 1}}{2}

Pressure Stabilization: The solver applies pressure stabilization in the pressure-correction equation to minimize non-physical pressure oscillations that arise from chequer-board effects. The user-specified pressure stabilization coefficient adjusts the magnitude of the stabilization term. The coefficient value ranges from 0 (no stabilization) to 1 (high stabilization). Large values of the pressure stabilization coefficient can lead to particles leaving through the walls.

(1122)