GitHub

Weakly Compressible SPH

Weakly compressible SPH as introduced by Monaghan (1994). This formulation relies on a stiff equation of state that generates large pressure changes for small density variations.

TrixiParticles.WeaklyCompressibleSPHSystem — Type
WeaklyCompressibleSPHSystem(initial_condition,
                            density_calculator, state_equation,
                            smoothing_kernel, smoothing_length;
                            acceleration=ntuple(_ -> 0.0, NDIMS),
                            viscosity=nothing, density_diffusion=nothing,
                            pressure_acceleration=nothing,
                            shifting_technique=nothing,
                            buffer_size=nothing,
                            correction=nothing, source_terms=nothing,
                            surface_tension=nothing, surface_normal_method=nothing,
                            reference_particle_spacing=0.0))

System for particles of a fluid. The weakly compressible SPH (WCSPH) scheme is used, wherein a stiff equation of state generates large pressure changes for small density variations. See Weakly Compressible SPH for more details on the method.

Arguments

Keyword Arguments

  • acceleration: Acceleration vector for the system. (default: zero vector)
  • viscosity: Viscosity model for this system (default: no viscosity). See ArtificialViscosityMonaghan or ViscosityAdami.
  • density_diffusion: Density diffusion terms for this system. See AbstractDensityDiffusion.
  • pressure_acceleration: Pressure acceleration formulation for this system. By default, the correct formulation is chosen based on the density calculator and the correction method. To use Tensile Instability Control, pass tensile_instability_control here.
  • shifting_technique: Shifting technique or transport velocity formulation to use with this system. Default is no shifting.
  • buffer_size: Number of buffer particles. This is needed when simulating with OpenBoundarySystem.
  • correction: Correction method used for this system. (default: no correction, see Corrections)
  • source_terms: Additional source terms for this system. Has to be either nothing (by default), or a function of (coords, velocity, density, pressure, t) (which are the quantities of a single particle), returning a Tuple or SVector that is to be added to the acceleration of that particle. See, for example, SourceTermDamping. Note that these source terms will not be used in the calculation of the boundary pressure when using a boundary with BoundaryModelDummyParticles and AdamiPressureExtrapolation. The keyword argument acceleration should be used instead for gravity-like source terms.
  • surface_tension: Surface tension model used for this SPH system. (default: no surface tension)
  • surface_normal_method: The surface normal method to be used for this SPH system. (default: no surface normal method or ColorfieldSurfaceNormal() if a surface_tension model is used)
  • reference_particle_spacing: The reference particle spacing used for weighting values at the boundary, which currently is only needed when using surface tension.
  • color_value: The value used to initialize the color of particles in the system.
source

Equation of State

The equation of state is used to relate fluid density to pressure and thus allow an explicit simulation of the WCSPH system. The equation in the following formulation was introduced by Cole (1948) (pp. 39 and 43). The pressure $p$ is calculated as

\[ p = B \left(\left(\frac{\rho}{\rho_0}\right)^\gamma - 1\right) + p_{\text{background}},\]

where $\rho$ denotes the density, $\rho_0$ the reference density, and $p_{\text{background}}$ the background pressure, which is set to zero when applied to free-surface flows (Adami et al., 2012).

The bulk modulus, $B = \frac{\rho_0 c^2}{\gamma}$, is calculated from the artificial speed of sound $c$ and the isentropic exponent $\gamma$.

An ideal gas equation of state with a linear relationship between pressure and density can be obtained by choosing exponent=1, i.e.

\[ p = B \left( \frac{\rho}{\rho_0} -1 \right) = c^2(\rho - \rho_0).\]

For higher Reynolds numbers, exponent=7 is recommended, whereas at lower Reynolds numbers exponent=1 yields more accurate pressure estimates since pressure and density are proportional (see Morris, 1997).

When using SummationDensity (or DensityReinitializationCallback) and free surfaces, initializing particles with equal spacing will cause underestimated density and therefore strong attractive forces between particles at the free surface. Setting clip_negative_pressure=true can avoid this.

TrixiParticles.StateEquationCole — Type
StateEquationCole(; sound_speed, reference_density, exponent,
                  background_pressure=0.0, clip_negative_pressure=false)

Equation of state to describe the relationship between pressure and density of water up to high pressures.

Keywords

  • sound_speed: Artificial speed of sound.
  • reference_density: Reference density of the fluid.
  • exponent: A value of 7 is usually used for most simulations.
  • background_pressure=0.0: Background pressure.
  • clip_negative_pressure=false: Negative pressure values are clipped to 0, which prevents spurious surface tension with SummationDensity but allows unphysical rarefaction of the fluid.
source
TrixiParticles.StateEquationIdealGas — Type
StateEquationIdealGas( ;sound_speed, reference_density, gamma, background_pressure=0.0,
                       clip_negative_pressure=false)

Equation of state to describe the relationship between pressure and density of a gas using the Ideal Gas Law.

Keywords

  • sound_speed : Artificial speed of sound.
  • reference_density : Reference density of the fluid.
  • gamma : Heat-capacity ratio
  • background_pressure=0.0 : Background pressure.
  • clip_negative_pressure=false: Negative pressure values are clipped to 0, which prevents spurious surface tension with SummationDensity but allows unphysical rarefaction of the fluid.
source

Density Diffusion

Density diffusion can be used with ContinuityDensity to remove the noise in the pressure field. It is highly recommended to use density diffusion when using WCSPH.

Formulation

All density diffusion terms extend the continuity equation (see ContinuityDensity) by an additional term

\[\frac{\mathrm{d}\rho_a}{\mathrm{d}t} = \sum_{b} m_b v_{ab} \cdot \nabla W_{ab} + \delta h c \sum_{b} V_b \psi_{ab} \cdot \nabla W_{ab},\]

where $V_b = m_b / \rho_b$ is the volume of particle $b$ and $\psi_{ab}$ depends on the density diffusion method (see AbstractDensityDiffusion for available terms). Also, $\rho_a$ denotes the density of particle $a$ and $r_{ab} = r_a - r_b$ is the difference of the coordinates, $v_{ab} = v_a - v_b$ of the velocities of particles $a$ and $b$.

Numerical Results

All density diffusion terms remove numerical noise in the pressure field and produce more accurate results than weakly commpressible SPH without density diffusion. This can be demonstrated with dam break examples in 2D and 3D. Here, $δ = 0.1$ has been used for all terms. Note that, due to added stability, the adaptive time integration method that was used here can choose higher time steps in the simulations with density diffusion. For the cheap DensityDiffusionMolteniColagrossi, this results in reduced runtime.

density_diffusion_2d
Dam break in 2D with different density diffusion terms
density_diffusion_3d
Dam break in 3D with different density diffusion terms

The simpler terms DensityDiffusionMolteniColagrossi and DensityDiffusionFerrari do not solve the hydrostatic problem and lead to incorrect solutions in long-running steady-state hydrostatic simulations with free surfaces (Antuono et al., 2012). This can be seen when running the simple rectangular tank example until $t = 40$ (again using $δ = 0.1$):

density_diffusion_tank
Tank in rest under gravity in 3D with different density diffusion terms

DensityDiffusionAntuono adds a correction term to solve this problem, but this term is very expensive and adds about 40–50% of computational cost.

API

TrixiParticles.DensityDiffusionAntuono — Type
DensityDiffusionAntuono(initial_condition; delta)

The commonly used density diffusion terms by Antuono (2010), also referred to as δ-SPH. The density diffusion term by Molteni (2009) is extended by a second term, which is nicely written down by Antuono (2012).

The term $\psi_{ab}$ in the continuity equation in AbstractDensityDiffusion is defined by

\[\psi_{ab} = 2\left(\rho_a - \rho_b - \frac{1}{2}\big(\nabla\rho^L_a + \nabla\rho^L_b\big) \cdot r_{ab}\right) \frac{r_{ab}}{\Vert r_{ab} \Vert^2},\]

where $\rho_a$ and $\rho_b$ denote the densities of particles $a$ and $b$ respectively and $r_{ab} = r_a - r_b$ is the difference of the coordinates of particles $a$ and $b$. The symbol $\nabla\rho^L_a$ denotes the renormalized density gradient defined as

\[\nabla\rho^L_a = -\sum_b (\rho_a - \rho_b) V_b L_a \nabla W_{ab}\]

with

\[L_a := \left( -\sum_{b} V_b r_{ab} \otimes \nabla W_{ab} \right)^{-1} \in \R^{d \times d},\]

where $d$ is the number of dimensions.

See AbstractDensityDiffusion for an overview and comparison of implemented density diffusion terms.

source
TrixiParticles.DensityDiffusionFerrari — Type
DensityDiffusionFerrari()

A density diffusion term by Ferrari (2009).

The term $\psi_{ab}$ in the continuity equation in AbstractDensityDiffusion is defined by

\[\psi_{ab} = \frac{\rho_a - \rho_b}{h_a + h_b} \frac{r_{ab}}{\Vert r_{ab} \Vert},\]

where $\rho_a$ and $\rho_b$ denote the densities of particles $a$ and $b$ respectively, $r_{ab} = r_a - r_b$ is the difference of the coordinates of particles $a$ and $b$ and $h_a$ and $h_b$ are the smoothing lengths of particles $a$ and $b$ respectively.

See AbstractDensityDiffusion for an overview and comparison of implemented density diffusion terms.

source
TrixiParticles.DensityDiffusionMolteniColagrossi — Type
DensityDiffusionMolteniColagrossi(; delta)

The commonly used density diffusion term by Molteni (2009).

The term $\psi_{ab}$ in the continuity equation in AbstractDensityDiffusion is defined by

\[\psi_{ab} = 2(\rho_a - \rho_b) \frac{r_{ab}}{\Vert r_{ab} \Vert^2},\]

where $\rho_a$ and $\rho_b$ denote the densities of particles $a$ and $b$ respectively and $r_{ab} = r_a - r_b$ is the difference of the coordinates of particles $a$ and $b$.

See AbstractDensityDiffusion for an overview and comparison of implemented density diffusion terms.

source

Particle Shifting Technique

The Particle Shifting Technique (PST) is used to correct tensile instability in regions of low pressure, as observed in viscous flow around an object. Without PST, tensile instability causes non-physical separation of the fluid from the object, introducing void regions behind the object.

At lower resolutions, PST alone can be effective to correct a viscous flow around a cylinder, as shown in this figure. particle_shifting At higher resolutions, PST alone is not effective anymore; see the figure in Tensile Instability Control. We recommend using PST and Tensile Instability Control together in such simulations.

Mathematical formulation

We use the following formulation by Sun et al. (2018). After each time step, a correction term $\delta r_a$ is added to the position $r_a$ of particle $a$, which is given by

\[\delta r_a = -4 \Delta t \, v_\text{max} h \sum_b \left( 1 + R \left( \frac{W_{ab}}{W(\Delta x_a)} \right)^n \right) \nabla W_{ab} \frac{m_b}{\rho_a + \rho_b},\]

where:

  • $\Delta t$ is the time step,
  • $v_\text{max}$ is the maximum velocity over all particles,
  • $h$ is the smoothing length,
  • $R$ and $n$ are constants, which are set to $0.2$ and $4$ respectively,
  • $W(\Delta x_a)$ is the smoothing kernel of the particle size of particle $a$, which can be interpreted as the target particle spacing that we want to achieve.
  • $\nabla W_{ab}$ is the gradient of the smoothing kernel,
  • $m_b$ is the mass of particle $b$,
  • $\rho_a, \rho_b$ is the density of particles $a$ and $b$, respectively.

Note that we replaced $\text{CFL} \cdot \text{Ma}$ by $\Delta t \cdot v_\text{max} / h$, as explained in Sun2018 on page 29, right above Equation 9.

The $\delta$-SPH method (WCSPH with density diffusion) together with this formulation of PST is commonly referred to as $\delta^+$-SPH.

To apply particle shifting, use the keyword argument shifting_technique in the constructor of a system that supports it.

Transport Velocity Formulation (TVF)

An alternative formulation is the so-called Transport Velocity Formulation (TVF) by Adami (2013). Ramachandran (2019) applied the TVF also for the EDAC scheme.

The transport velocity $\tilde{v}_a$ of particle $a$ is used to evolve the position of the particle $r_a$ from one time step to the next by

\[\frac{\mathrm{d} r_a}{\mathrm{d}t} = \tilde{v}_a\]

and is obtained at every time step $\Delta t$ from

\[\tilde{v}_a (t + \Delta t) = v_a (t) + \Delta t \left(\frac{\tilde{\mathrm{d}} v_a}{\mathrm{d}t} - \frac{1}{\rho_a} \nabla p_{\text{background}} \right),\]

where $\rho_a$ is the density of particle $a$ and $p_{\text{background}}$ is a constant background pressure field. The tilde in the second term of the right-hand side indicates that the material derivative has an advection part.

The discretized form of the last term is

\[ -\frac{1}{\rho_a} \nabla p_{\text{background}} \approx -\frac{p_{\text{background}}}{m_a} \sum_b \left(V_a^2 + V_b^2 \right) \nabla_a W_{ab},\]

where $V_a$, $V_b$ denote the volume of particles $a$ and $b$ respectively. Note that although in the continuous case $\nabla p_{\text{background}} = 0$, the discretization is not 0th-order consistent for non-uniform particle distribution, which means that there is a non-vanishing contribution only when particles are disordered. That also means that $p_{\text{background}}$ occurs as pre-factor to correct the trajectory of a particle resulting in uniform pressure distributions. Suggested is a background pressure which is in the order of the reference pressure, but it can be chosen arbitrarily large when the time-step criterion is adjusted.

The inviscid momentum equation with an additional convection term for a particle moving with $\tilde{v}$ is

\[\frac{\tilde{\mathrm{d}} \left( \rho v \right)}{\mathrm{d}t} = -\nabla p + \nabla \cdot \bm{A},\]

where the tensor $\bm{A} = \rho v\left(\tilde{v}-v\right)^T$ is a consequence of the modified advection velocity and can be interpreted as the convection of momentum with the relative velocity $\tilde{v}-v$.

The discretized form of the momentum equation for a particle $a$ reads as

\[\frac{\tilde{\mathrm{d}} v_a}{\mathrm{d}t} = \frac{1}{m_a} \sum_b \left(V_a^2 + V_b^2 \right) \left[ -\tilde{p}_{ab} \nabla_a W_{ab} + \frac{1}{2} \left(\bm{A}_a + \bm{A}_b \right) \cdot \nabla_a W_{ab} \right].\]

Here, $\tilde{p}_{ab}$ is the density-weighted pressure

\[\tilde{p}_{ab} = \frac{\rho_b p_a + \rho_a p_b}{\rho_a + \rho_b},\]

with the density $\rho_a$, $\rho_b$ and the pressure $p_a$, $p_b$ of particles $a$ and $b$, respectively. $\bm{A}_a$ and $\bm{A}_b$ are the convection tensors for particle $a$ and $b$, respectively, and are given, e.g., for particle $a$, as $\bm{A}_a = \rho v_a\left(\tilde{v}_a-v_a\right)^T$.

To apply the TVF, use the keyword argument shifting_technique in the constructor of a system that supports it.

TrixiParticles.ParticleShiftingTechnique — Type
ParticleShiftingTechnique(; integrate_shifting_velocity=true,
                          update_everystage=false,
                          modify_continuity_equation=true,
                          second_continuity_equation_term=ContinuityEquationTermSun2019(),
                          momentum_equation_term=MomentumEquationTermSun2019(),
                          v_max_factor=1, sound_speed_factor=0)

Particle Shifting Technique by Sun et al. (2017) and Sun et al. (2019). The keyword arguments allow to choose between the original methods from the two papers and variants in between.

The default values of the keyword arguments provide the version of shifting that we recommend based on our experiments. The default values are subject to change in future releases.

See Particle Shifting Technique for more information on the method.

We provide the following convenience constructors for common variants of the method:

  • ParticleShiftingTechniqueSun2017(): Particle Shifting Technique by Sun et al. (2017). Shifting is applied as a position correction in a callback after each time step. No additional terms are added to the momentum or continuity equations.
  • ConsistentShiftingSun2019(): Consistent Particle Shifting Technique by Sun et al. (2019). Shifting is applied with a shifting velocity in each stage of the time integration. Additional terms are added to the momentum and continuity equations, most importantly to guarantee conservation of volume in closed systems, which is not the case for the original method by Sun et al. (2017).

Keywords

  • integrate_shifting_velocity: If true, the shifting is applied in each stage of the time integration method as a shifting velocity that is added to the physical velocity in the time integration. If false, the shifting is applied as a position correction in a callback after each time step.
  • update_everystage: If true, the shifting velocity is updated in every stage of a multi-stage time integration method. This requires integrate_shifting_velocity=true. If false, the shifting velocity is only updated once per time step in a callback, and the same shifting velocity is used for all stages. update_everystage=false reduces the computational cost, but may reduce the stability of the scheme and require a smaller time step.
  • modify_continuity_equation: If true, the continuity equation is modified to be based on the transport velocity instead of the physical velocity. This guarantees conservation of volume in closed systems, but is unstable at solid wall boundaries, according to our experiments. This requires integrate_shifting_velocity=true.
  • second_continuity_equation_term: Additional term to be added to the continuity equation to solve the stability problems with the modified continuity equation at solid wall boundaries. This requires modify_continuity_equation=true. See ContinuityEquationTermSun2019.
  • momentum_equation_term: Additional term to be added to the momentum equation to account for the shifting velocity. This requires integrate_shifting_velocity=true. See MomentumEquationTermSun2019.
  • v_max_factor: Factor to scale the expected maximum velocity used in the shifting velocity. The maximum expected velocity is computed as v_max_factor * max(|v|), where v is the physical velocity. As opposed to sound_speed_factor, the computed expected maximum velocity depends on the current flow field and can change over time. Only one of v_max_factor and sound_speed_factor can be non-zero.
  • sound_speed_factor: Factor to compute the maximum expected velocity used in the shifting velocity from the speed of sound. The maximum expected velocity is computed as sound_speed_factor * c, where c is the speed of sound. Only one of v_max_factor and sound_speed_factor can be non-zero.
Warning

The Particle Shifting Technique needs to be disabled close to the free surface and therefore requires a free surface detection method. This is not yet implemented. This technique cannot be used in a free surface simulation.

source
TrixiParticles.TransportVelocityAdami — Type
TransportVelocityAdami(; background_pressure::Real, modify_continuity_equation=false)

Transport Velocity Formulation (TVF) by Adami et al. (2013) to suppress pairing and tensile instability. See TVF for more details of the method.

Keywords

  • background_pressure: Background pressure. Suggested is a background pressure which is on the order of the reference pressure.
  • modify_continuity_equation: If true, the continuity equation is modified to be based on the transport velocity instead of the physical velocity. This guarantees conservation of volume in closed systems, but is unstable at solid wall boundaries, according to our experiments.
Warning

The Transport Velocity Formulation needs to be disabled close to the free surface and therefore requires a free surface detection method. This is not yet implemented. This technique cannot be used in a free surface simulation.

source
TrixiParticles.ConsistentShiftingSun2019 — Method
ConsistentShiftingSun2019(; sound_speed_factor=0.1, kwargs...)

Consistent Particle Shifting Technique by Sun et al. (2019).

This is a convenience constructor for:

ParticleShiftingTechnique(integrate_shifting_velocity=true,
                          update_everystage=true,
                          modify_continuity_equation=true,
                          second_continuity_equation_term=ContinuityEquationTermSun2019(),
                          momentum_equation_term=MomentumEquationTermSun2019(),
                          v_max_factor=0, sound_speed_factor=0.1)

See ParticleShiftingTechnique for all available options.

Keywords

  • sound_speed_factor: Factor to compute the maximum expected velocity used in the shifting velocity from the speed of sound. The maximum expected velocity is computed as sound_speed_factor * c, where c is the speed of sound. Since the speed of sound is usually chosen as 10 times the maximum expected velocity in weakly compressible SPH, a value of 0.1 corresponds to the maximum expected velocity.
  • kwargs...: All keywords are passed to the main constructor.

Examples

shifting_technique = ConsistentShiftingSun2019()
Warning

The Particle Shifting Technique needs to be disabled close to the free surface and therefore requires a free surface detection method. This is not yet implemented. This technique cannot be used in a free surface simulation.

source
TrixiParticles.ParticleShiftingTechniqueSun2017 — Method
ParticleShiftingTechniqueSun2017(; kwargs...)

Particle Shifting Technique by Sun et al. (2017). Following the original paper, the callback is applied in every time step and not in every stage of a multi-stage time integration method to reduce the computational cost.

This is a convenience constructor for:

ParticleShiftingTechnique(integrate_shifting_velocity=false,
                          update_everystage=false,
                          modify_continuity_equation=false,
                          second_continuity_equation_term=nothing,
                          momentum_equation_term=nothing,
                          v_max_factor=1, sound_speed_factor=0)

See ParticleShiftingTechnique for all available options.

Keywords

  • kwargs...: All keywords are passed to the main constructor.

Examples

shifting_technique = ParticleShiftingTechniqueSun2017()
Warning

The Particle Shifting Technique needs to be disabled close to the free surface and therefore requires a free surface detection method. This is not yet implemented. This technique cannot be used in a free surface simulation.

source

Tensile Instability Control

Tensile Instability Control (TIC) is a pressure acceleration formulation to correct tensile instability in regions of low pressure, as observed in viscous flow around an object. The technique was introduced by Sun et al. (2018). The formulation is described in Section 2.1 of this paper. It can be used in combination with the Particle Shifting Technique (PST) to effectively prevent non-physical separation of the fluid from the object.

As can be seen in the following figure, TIC alone can cause instabilities and does not improve the simulation. PST alone can mostly prevent separation at lower resolutions. A small void region is still visible, but quickly disappears. The combination of PST and TIC prevents separation effectively. low_res

At higher resolutions, PST alone is not effective to prevent separation, as can be seen in the next figure. Only the combination of PST and TIC is able to produce physical results. high_res

Mathematical formulation

The force that particle $a$ experiences from particle $b$ due to pressure is given by

\[f_{ab} = -m_a m_b \frac{p_a + p_b}{\rho_a \rho_b} \nabla W_{ab}\]

for the WCSPH method with ContinuityDensity.

The TIC formulation changes this force to

\[f_{ab} = -m_a m_b \frac{|p_a| + p_b}{\rho_a \rho_b} \nabla W_{ab}.\]

Note that this formulation is asymmetric and sacrifices conservation of linear and angular momentum.

TrixiParticles.tensile_instability_control — Function
tensile_instability_control

Pressure acceleration formulation to prevent tensile instability by Sun et al. (2018).

This formulation can be passed as keyword argument pressure_acceleration to the WeaklyCompressibleSPHSystem constructor. See Tensile Instability Control for more information on this technique.

Warning

Tensile Instability Control needs to be disabled close to the free surface and therefore requires a free surface detection method. This is not yet implemented. This technique cannot be used in a free surface simulation.

source