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 between two particles i MathType@MTEF@5@5@+= feaahqart1ev3aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiVCI8FfYJH8YrFfeuY=Hhbbf9v8qqaqFr0xc9pk0xbb a9q8WqFfeaY=biLkVcLq=JHqpepeea0=as0Fb9pgeaYRXxe9vr0=vr 0=vqpWqaaeaabiGaciaacaqabeaadaqaaqaaaOqaaabaaaaaaaaape GaamyAaaaa@36F7@ and j MathType@MTEF@5@5@+= feaahqart1ev3aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiVCI8FfYJH8YrFfeuY=Hhbbf9v8qqaqFr0xc9pk0xbb a9q8WqFfeaY=biLkVcLq=JHqpepeea0=as0Fb9pgeaYRXxe9vr0=vr 0=vqpWqaaeaabiGaciaacaqabeaadaqaaqaaaOqaaabaaaaaaaaape GaamOAaaaa@36F8@ . 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 ) and pressure ( p ) are expressed as:(1)
( ρ L , U L , p L ) = ( ρ i , v i e i j , p i ) ( ρ R , U R , p R ) = ( ρ j , v j e i j , p j )
Where, v is the velocity vector of the respective particle, which is illustrated in Figure 1. C. Zhang
Figure 1. Construction of the Riemann problem in SPH

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 and p L = p R = p , as indicated by the dashed link in Figure 2.


Figure 2. Simplified Riemann Problem Solution in Time. This example shows three possible waves.C. Zhang
From these basic assumptions, you can derive the U * and p * values as:(2)
U * = U a v g + 1 2 ( p L p R ) p a v g C 0
(3)
p * = p avg + 1 2 p avg C 0 ( U L U R )

Where, the 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)
d p i d t = 2 p i j m j p j ( v i v * ) i W i j d v i d t = 2 j m j ( p * p i p j ) i W i j

Where, 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)
p i = j m j W i j j V j W i j
Where, p i is the density of particle I. m j and V j is the mass and volume of particle j respectively. W i j is the 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)
p i = j m j W i j min ( 1 , j m j p * W i j )
Where, p i is the density from solving the continuity equation before applying the filter.