3. Neutrino Physics

This document is a rough guide for the internal physics implementation of the Neutrino SPH solver. More details of the implementation could be found in [SAMPATH2016], [AKINCI2012] and [IHMSEN2014] and [MKSIISPH].

3.1. Introduction

Smoothed Particle hydrodynamics (SPH) is a common mesh-free Lagrangian method. It was originally introduced for dealing with astrophysical, compressible fluid flows [MM1994], and later, its range of application domain extended to a wider variety of flows, such as non-Newtonian fluids , granular flows , and even solid deformation and fracture . SPH has been applied to incompressible, free-surface flows for more than two decades along with the closely related moving particle semi-implicit (MPS) method. Mainly, two approaches have been used in SPH for enforcing incompressibility: (1) state equation-based SPH (SESPH), in which a compressible approach is employed with an equation of state, relating density, pressure and speed of sound in the fluid, and (2) in incompressible SPH (ISPH), in which pressure forces are implicitly computed and iteratively refined by solving a pressure Poisson equation until the desired compressibility is reached.

SPH requires the computation of sums over dynamically changing sets of neighboring particles. Spatial data structures are employed to accelerate the neighborhood search, and efficient construction and query are thus essential. While hierarchical data structures are used in simulations with a variable particle interaction radius , uniform grids are more efficient when the radius is fixed. The cell size is usually set to the interaction radius. Index sort-based grids are often preferred, because the basic grid suffers from low-cache-hit rate and race conditions during parallel construction.

3.2. SPH Theory

SPH is a numerical method for both interpolating quantities and approximating spatial differential operators from a set of known quantities at sampled positions. Among the different ways to interpolate from scattered data, SPH relies on the use of radial basis functions. The value of a field at any point of space is estimated based on the value of this field at neighbor sampling points and its distance to them.

Usually a particle in physics is a minute quantity of matter and could be one of the following

  • Smallest constituents of matter (Standard Model)
  • Nanoparticles, colloidal particles
  • Dust, powder, ashes
  • Sediment grains, water droplets etc.
However the particles in SPH have the following properties
  • They are material points
  • They have volume, mass, pressure, density, etc.

3.3. SPH Approximation

Particle a has position ra, mass ma, volume Va, etc. * Particle interaction are computed using the ’kernel’ w(r) * Cubic Spline Kernel is used in Neutrino

3.4. SPH Kernels

The kernel function used for interpolation in SPH must satisfy a certain number of properties. * It must be defined on a compact support * It must also be radial, feature a Gaussian-like shape, be sufficiently smooth * It is even and normalized leads to a first-order consistent continuous SPH interpolation. * Note that the particle distribution has an impact on the accuracy. * A regular and isotropic distribution is desired.

In Neutrino, the cubic B-spline kernel like the following is employed.


Fig. 3.1 Cubic Spline Kernel

3.5. SPH Navier-Stokes Equations

The single-phase, isothermal, incompressible, Newtonian fluid flows are the considered physical model. In addition, the viscosity is assumed constant in space, and the surface tension forces are ignored. The Navier-Stokes equations in their Lagrangian, velocity-pressure formulation then reads:

(3.1)\[\begin{split}\frac{dv}{dt} & =-\frac{1}{\rho}\nabla{} p+\nu\boldsymbol{\nabla^2}\text{v}+\text{g}\\\end{split}\]
(3.2)\[\begin{split}\frac{d\rho}{dt} & =-\rho\boldsymbol{\nabla}.\boldsymbol{v}\\\end{split}\]
\[\mathbf{v} = velocity\]
\[\mathbf{p} = pressure\]
\[\mathbf{\nu} = viscosity\]


\[\mathbf{g} = gravity\]

Equation (3.1) comes from the momentum conservation, and involves pressure, viscosity, and gravity forces. Equation (3.2) is derived from the mass conservation. The incompressibility condition is satisfied either by

\[\frac{d\rho}{dt} = 0\]


\[\boldsymbol{\nabla}.\boldsymbol{v} = 0;\]

these two ways are equivalent in the continuous representation.

3.6. SPH Solvers

3.6.1. Implicit Incompressible SPH (IISPH)

Neutrino base solver is named as IISPH. IISPH solver computes pressure by iteratively solving a linear system to meet the incompressibility criterion. IISPH is based on semi-implicit form of the density prediction using the time rate of change of the density. By discretizing the continuity equation

\[\frac{{D\rho}}{{Dt}} = - \rho\nabla\]
\[\frac{\rho_{i}(t + t)}{t} = \sum_{j}^{}{m_{j}\mathbf{v}_{\mathbf{\text{ij}}}\left( t + t \right)\mathbf{\bullet \nabla}W_{\text{ij}}(t)}\]

The unknown velocities \(\mathbf{v}_{\mathbf{\text{ij}}}\left( t + t \right)\) depend on unknown pressure values at time , it’s firstly predicted using any known non-pressure forces such as gravity, surface tension and viscosity. Then the predicted velocity is used to determine a predicted density.

\[\mathbf{v}_{\mathbf{\text{ij}}}^{\text{adv}} = \mathbf{v}_{\mathbf{i}}\left( t \right) + t\frac{\mathbf{F}_{i}^{\text{adv}}(i)}{m_{i}}\]
\[\rho_{i}^{\text{adv}} = \rho_{i}\left( t \right) + t\sum_{j}^{}m_{j}\mathbf{v}_{\mathbf{\text{ij}}}^{\text{adv}}\mathbf{\bullet \nabla}W_{\text{ij}}(t)\]

Now the pressure forces to resolve the compression is searched based on:

\[t^{2}\sum_{j}^{}{m_{j}\left( \frac{\mathbf{F}_{i}^{p}(t)}{m_{i}} - \frac{\mathbf{F}_{j}^{p}(t)}{m_{j}} \right)\mathbf{\bullet \nabla}W_{\text{ij}}\left( t \right) = \rho_{0} - \rho_{i}^{\text{adv}}}\]

Where \(\mathbf{F}_{i}^{p}(t)\) is calculated by1:

\[\mathbf{F}_{i}^{p}\left( t \right) = - m_{i}\sum_{j}^{}m_{j}\left( \frac{p_{i}\left( t \right)}{\rho_{i}^{2}\left( t \right)} + \frac{p_{j}\left( t \right)}{\rho_{j}^{2}\left( t \right)} \right)\mathbf{\nabla}W_{\text{ij}}(t)\]

From this point, a linear system can be got with unknown pressure value for each particle with the form:

\[\sum_{j}^{}{a_{\text{ij}}p_{j} = b_{i} = \rho_{o} - \rho_{i}^{\text{adv}}}\]

With the calculated pressure on each particle, the unknown velocity \(\mathbf{v}_{\mathbf{\text{ij}}}\left( t + t \right)\) can be calculated as:

\[\mathbf{v}_{\mathbf{i}}\left( t + t \right) = \mathbf{v}_{\mathbf{\text{ij}}}^{\text{adv}} + t\mathbf{F}_{\mathbf{i}}^{p}(t)/m_{i}\]

Then we will need this velocity to calculate the position of each particle. In Neutrino, different time integration schemes are applied.

3.6.2. Time Integration a. Euler Cromer Integration

In Euler Cromer integration scheme, the new position of particle is calculated as:

\[\mathbf{x}_{\mathbf{i}}\left( t + t \right) = \mathbf{x}_{\mathbf{i}}(t) + t\mathbf{v}_{i}(t + t)\] Verlet Integration

Verlet scheme calculate the new position as:

\[\mathbf{x}_{\mathbf{i}}\left( t + t \right) = \mathbf{x}_{\mathbf{i}}\left( t \right) + t\mathbf{v}_{i}\left( t + t \right) + 0.5t^{2}\mathbf{F}_{i}\]

The numerical tests show Verlet scheme has less numerical dissipation than Euler Cromer, this is because Verlet has a higher order calculation of position.

3.7. Momentum Equation

The conservation of momentum equation can be described as:

\[\frac{D\mathbf{v}}{\text{Dt}} = - \frac{1}{\rho}\mathbf{\nabla}P + \mathbf{g} + v_{0}\nabla^{2}\mathbf{v}_{\text{ab}}\]

\(\tau\) represents the viscosity term. In Neutrino, two different viscosity models are applied: a. Artificial viscosity; b. Laminar Viscosity.

3.7.1. Artificial Viscosity

The particle acceleration due to artificial viscosity proposed by Monaghan[MM1994] is defined as:


\[\mu_{\text{ab}} = \frac{h\mathbf{v}_{\text{ab}}\mathbf{r}_{\text{ab}}}{\mathbf{r}_{\text{ab}}^{2} + \eta}\]

Distance between particles a and b is


Average sound speed of particle pair a & b

\[\overset{\overline{}}{c_{\text{ab}}} = \frac{c_{a} + c_{b}}{2}\]

Average density of particle pair a & b

\[\overset{\overline{}}{\rho_{\text{ab}}} = \frac{\rho_{a} + \rho_{b}}{2}\]

\(\eta = 0.01h^{2},\ \alpha\) is the artificial viscosity parameter and it can be understand in term of kinematic viscosity as[LIU03]:

\[v_{0} = \frac{\text{αhc}}{8}\]

3.7.2. Laminar Viscosity

The laminar viscosity proposed by Morris[MORRIS97] is defined as:

\[\frac{d{{v}}_{a}}{{dt}} = - \sum_{b}^{}{m_{b}\left( \frac{4v_{0}\mathbf{r}_{\text{ab}}\mathbf{\bullet}\mathbf{\nabla}W_{{ab}}}{{(\rho}_{a} + \rho_{b})\mathbf{r}_{\text{ab}}^{2}} \right)}\mathbf{v}_{{ab}}\]
\[v_{0} \text{ is the kinematic viscosity of fluid simulated}\]

3.8. Boundary Handling

The fluid-rigid boundaries require special attention. First, discontinuities of physical quantities that occur at boundaries are problematic for the usual forms of SPH. Proper boundary handling is necessary to avoid underestimated densities and non-physical pressure forces. Furthermore, pressure and friction forces between the fluid and the rigid bodies must be accounted for, and non-penetration must be ensured. Handling of thin objects or with complex geometries is particularly challenging, as well as when multiple dynamic objects are involved. Different strategies have been proposed through use of distance-based penalty forces, frozen particles , mirror particles [LIU03] or a wall renormalisation factor.

In Neutrino an efficient technique based on frozen particles and able to cope with complex boundaries, multiple bodies and the discontinuity issue is employed. Moreover, the technique uses fluid-rigid pressure and friction force models conserving linear and angular momentum.

3.8.1. Single Layer Boundary Mode

Appropriate adjustments are made to the boundary conditions to ensure the particle deficiency in using a single layer of particles.

The boundary coefficient is impacting the magnitude of the pressure forces, and it was empirically calibrated.

3.8.2. Multi Layer Boundary Mode

Currently Under Development

3.9. Fluid Structure Interactions

Using a standard SPH method in the sense that the time-integrated scheme works as classically done, with the following chronological steps (at each time step) :

  1. compute the accelerations using the positions and velocities;
  2. update the velocities using the accelerations;
  3. update the positions using the velocities;

Typically, the steps 2) and 3) consists in a variant of Euler, Runge-Kutta, or velocity Verlet method.

3.9.1. IISPH vs Position Based Dynamics (PBD)

Conversely, PBD reformulates the force computations as a set of constraints on the positions and updates the positions directly without using the velocities. 1) compute the constraints using the positions and velocities; 2) update the positions by solving for the constraints; 3) update the velocities using the positions. Typically, the Stormer-Verlet method, or a higher-order equivalent one, is used.

3.10. Sample Cases & Validations

3.10.1. Dam Break

Dam break is simulated by lots of SPH program or software to demonstrate the SPH’s power of dealing free-surface slamming phenomena. In this case, we measure the force exerted by fluid particles onto the dam structure and compare the results with experimental data. The simulation is one-to-one scale to real experiments and set up as in Figure 1 according to Cummins work[CUMMINS12]. In neutrino, we firstly put the gate in position and wait for 1 sec until all fluid particles are settled down. Then we open the gate and let the fluid flow under the effects of gravity.


Figure : Comparison of the forces (Simulated Vs Experiment)


Figure : schematic diagram of the dam geometry[CUMMINS12]

Dual-Sphysics, Neutrino and LAMMPS-SPH can all fulfill dam-break case, which make it possible to compare their accuracy and time consuming. Besides, the force exerted from fluid particles to rigid body is an important factor to our purpose, by comparing the simulation data to the real experimental data, the feature of each software can be clearly seen. Table 2 shows the comparison of key parameters for each program. A comparison of simulated force to experimental data is shown in Figure 2 and a good agreement is found.

Table : comparison of experimental measurement, neutrino, LAMMPS-SPH and Dual-SPHysics

  Force Measurement # Particles Time Consumption (sec/step) Avg. unit Time step (sec)
Experimental Data Yes      
Neutrino Yes 67053 0.057 (4 cores) 0.00102
LAMMPS-SPH Need post-processing 64906 0.136 (8 cores) 0.00025
Dual-SPHysics Need post-processing 116795 1.454 (4 cores) 0.00004

Figure : comparison of Neutrino output to experimental data


Figure: Comparison of forces using several particle sizes in Neutrino.

The first peak, representing the first slamming from fluid to the dam structure, is higher than the experimental data. This is mainly caused by the repulsive boundary treatment at fluid-solid interface. This also shows SPH is capable of determining fluid-solid interaction accurately.


Figure: Estimate on the error based on the particle resolution

3.10.2. Aureli Dam Break

The Aureli dam break experiment is based the experiment from [Aureli2015]. The impact forces are about 5x lesser than the previous damn break experiment and there’s a significant amount of air entrapment.

The measurement of the impact forces is done by smoothing out the highest frequency noise . It seems like the forces are overestimated by about 1N and the secondary shock not captured very well.


3.10.3. Poisuelle Flow

3.10.4. Faltisen - Wave Sloshing Experiment

3.10.5. Falling Body in Water

3.10.6. Solitary Wave Past Shore

3.11. References

[SAMPATH2016]Sampath, R. Montanari, N. Akinci, N. Prescott, S. Smith, C., Large-scale solitary wave simulation with implicit incompressible SPH, Journal of Ocean Engineering and Marine Energy, 2016, 1-17, 10.1007/s40722-016-0060-8, http://dx.doi.org/10.1007/s40722-016-0060-8
[AKINCI2012]N. Akinci, M. Ihmsen, G. Akinci, B. Solenthaler, M. Teschner, “Versatile Rigid-Fluid Coupling for Incompressible SPH”, ACM Transactions on Graphics (Proc. SIGGRAPH 2012), vol. 31, no. 4, pp. 62:1-62:8, July 2012
[IHMSEN2014]Ihmsen, M. and Cornelis, J. Solenthaler, B. Horvath, C. Teschner, M, Implicit Incompressible {SPH}, IEEE Transactions on Visualization and Computer Graphics, 2014
[MM1994](1, 2) Monaghan, J. J., “Smoothed particle hydrodynamics”, Annual Rev. Astron. Appl., 30: 543- 574., 1992
[MKSIISPH]Markus, I, etc., Implicit Incompressible SPH,
[MORRIS97]Morris, J. P, etc., “Modeling Low Reynolds Number Incompressible Flows Using SPH”, Journal of Computational Physics 136, 214-226., 1997
[LIU03](1, 2) Liu. G. R., etc., “Smoothed Particl Hydrodynamics: a meshfree particle method”, World Scientific, 2003.
[CUMMINS12](1, 2) S. J. Cummins., etc., “Three-dimensional wave impact on a rigid structure using smoothed particle hydrodynamics”, International Journal for Numerical Methods in Fluids (2012); 68: 1471-1496, http://onlinelibrary.wiley.com/doi/10.1002/fld.2539/epdf
[CHERN]Chern, M.J., Borthwick, A.G.L. and Eatock Taylor, R. “Pseudospectral element model for free surface viscous flows”, Int. J. Num. Meth. For Heat & Fluid Flow (2005), 15(6), 517 – 554, http://www.researchgate.net/publication/235263308_Pseudospectral_element_model_for_free_surface_viscous_flows
[GHIA82]Ghia, U., etc., “High-Re solutions for incompressible flow using the Navier-Stokes equations and a multigrid method”, Journal of Computational Physics, v. 5, n. 48, p.387-411, 1982
[Aureli2015]Aureli, Francesca and Dazzi, Susanna and Maranzoni, Andrea and Mignosa, Paolo and Vacondio, Renato, 10.1016/j.advwatres.2014.11.009}, Adv. Water Resources, Volume 76, January, 2015