Riemann Particle Interaction Scheme

The Riemann problem can be defined as a category of initial value problems that involve a conservation equation and a piecewise data set with a single discontinuity.

In terms of SPH, you can solve a Riemann problem for each pair of particles, yielding a more isotropic particle distribution and smoother gradients in the solution. The following depiction is based on the work of C. Zhang. Given the Riemann problem definition, you start off by constructing a unit vector

e_{i j}

$e_{i j}$ between two particles

i

$i$ and

j

$j$ . The two particles can then be assumed to be the left and right sides of the Riemann problem, such that the physical values of density (

ρ

$ρ$ ), velocity (

U

$U$ ) and pressure (

p

$p$ ) are expressed as:(1)

\begin{array}{l} (ρ_{L}, U_{L}, p_{L}) = (ρ_{i}, v_{i} \cdot e_{i j}, p_{i}) \\ (ρ_{R}, U_{R}, p_{R}) = (ρ_{j}, v_{j} \cdot e_{i j}, p_{j}) \end{array}

$\begin{array}{l} (ρ_{L}, U_{L}, p_{L}) = (ρ_{i}, v_{i} \cdot e_{i j}, p_{i}) \\ (ρ_{R}, U_{R}, p_{R}) = (ρ_{j}, v_{j} \cdot e_{i j}, p_{j}) \end{array}$

Where,

v

$v$ is the velocity vector of the respective particle, which is illustrated in Figure 1. C. Zhang

The discontinuity is assumed to be located at the interface (Figure 1), half way between the particles. Assuming that particles are moving with arbitrary velocity vectors, the consequent solution of their motion is either a rarefaction or compression wave. However, there is a third wave solution at the interface denoted with * exponent, where

U_{L}^{*} = U_{R}^{*} = U_{}^{*}

$U_{L}^{*} = U_{R}^{*} = U_{}^{*}$ and

p_{L}^{*} = p_{R}^{*} = p

$p_{L}^{*} = p_{R}^{*} = p$ , as indicated by the dashed link in Figure 2.

From these basic assumptions, you can derive the

U^{*}

$U^{*}$ and

p^{*}

$p^{*}$ values as:(2)

U^{*} = U_{a v g} + \frac{1}{2} \frac{(p_{L} - p_{R})}{p_{a v g} C_{0}}

$U^{*} = U_{a v g} + \frac{1}{2} \frac{(p_{L} - p_{R})}{p_{a v g} C_{0}}$

(3)

p^{*} = p_{a v g} + \frac{1}{2} p_{a v g} C_{0} (U_{L} - U_{R})

$p^{*} = p_{a v g} + \frac{1}{2} p_{a v g} C_{0} (U_{L} - U_{R})$

Where, the $a v g$ $a v g$ index denotes the mean value of the respective value between the two particles.

With these two solutions to the Riemann problem, you can modify the continuity equation and the pressure gradient term in the momentum equation as follows:(4)

\begin{array}{l} \frac{d p_{i}}{d t} = 2 p_{i} \sum_{j} \frac{m_{j}}{p_{j}} (v_{i} - v^{*}) \cdot \nabla_{i} W_{i j} \\ \frac{d v_{i}}{d t} = 2 \sum_{j} m_{j} (\frac{p^{*}}{p_{i} p_{j}}) \nabla_{i} W_{i j} \end{array}

$\begin{array}{l} \frac{d p_{i}}{d t} = 2 p_{i} \sum_{j} \frac{m_{j}}{p_{j}} (v_{i} - v^{*}) \cdot \nabla_{i} W_{i j} \\ \frac{d v_{i}}{d t} = 2 \sum_{j} m_{j} (\frac{p^{*}}{p_{i} p_{j}}) \nabla_{i} W_{i j} \end{array}$

Where, $v^{*} = U^{*} e_{i j} + (v_{a v g - i j} - U_{a v g} e_{i j})$ $v^{*} = U^{*} e_{i j} + (v_{a v g - i j} - U_{a v g} e_{i j})$ .

It is well-known that 1st order Riemann solvers such as the one illustrated here are dissipative. To that extent, nanoFluidX is using a limiter so that the Riemann solutions are provided only when the fluid is under compression. This effectively reduces the dissipation and has been shown to improve the results in various situations. C. Zhang

Note: Although the Riemann solver option is experimental, canonical test cases, simple geometries, and problems involving single phase runs such as tank sloshing or water management should still greatly benefit from the new formulation. However, industrial gearboxes and violent multiphase flows in complex moving geometries in general are likely to experience local instabilities. This will be addressed in the coming versions of nanoFluidX.

Density filtering is an integral part of the Riemann solver implementation. In SPH simulation, the consistency between mass, density, and occupied area cannot be strictly enforced. Several density filtering schemes have been developed to tackle this. In nanoFluidX, two versions are implemented: INCREMENTAL and INSTANT filtering.

INSTANT: The instant density filtering is basically a 0th-order Shepard filter. At each evaluation point (according to the number set by RM_freq_rho_init), the density is re-calculated by:(5)
$ρ_{i} = \frac{\sum_{j} m_{j} W_{i j}}{\sum_{j} V_{j} W_{i j}}$ $ρ_{i} = \frac{\sum_{j} m_{j} W_{i j}}{\sum_{j} V_{j} W_{i j}}$
Where,

$ρ_{i}$ $ρ_{i}$

Density of particle I

$m_{j}$ $m_{j}$

Mass of particle j

$V_{j}$ $V_{j}$

Volume of particle j

$W_{i j}$ $W_{i j}$

Weight calculated from the kernel function
INCREMENTAL: A modified 1st-order density filtering scheme is used as an additional option. The operator is calculated following:(6)
$ρ_{i} = \frac{\sum_{j} m_{j} W_{i j}}{\min (1, \sum_{j} \frac{m_{j}}{ρ^{*}} W_{i j})}$ $ρ_{i} = \frac{\sum_{j} m_{j} W_{i j}}{\min (1, \sum_{j} \frac{m_{j}}{ρ^{*}} W_{i j})}$
Where, $ρ_{i}$ $ρ_{i}$ is the density from solving the continuity equation before applying the filter.