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
— TypeWeaklyCompressibleSPHSystem(initial_condition,
density_calculator, state_equation,
smoothing_kernel, smoothing_length;
acceleration=ntuple(_ -> 0.0, NDIMS),
viscosity=nothing, density_diffusion=nothing,
pressure_acceleration=nothing,
transport_velocity=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
initial_condition
:InitialCondition
representing the system's particles.density_calculator
: Density calculator for the system. SeeContinuityDensity
andSummationDensity
.state_equation
: Equation of state for the system. SeeStateEquationCole
.smoothing_kernel
: Smoothing kernel to be used for this system. See Smoothing Kernels.smoothing_length
: Smoothing length to be used for this system. See Smoothing Kernels.
Keyword Arguments
acceleration
: Acceleration vector for the system. (default: zero vector)viscosity
: Viscosity model for this system (default: no viscosity). SeeArtificialViscosityMonaghan
orViscosityAdami
.density_diffusion
: Density diffusion terms for this system. SeeDensityDiffusion
.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, passtensile_instability_control
here.transport_velocity
: Transport Velocity Formulation (TVF). Default is no TVF.buffer_size
: Number of buffer particles. This is needed when simulating withOpenBoundarySPHSystem
.correction
: Correction method used for this system. (default: no correction, see Corrections)source_terms
: Additional source terms for this system. Has to be eithernothing
(by default), or a function of(coords, velocity, density, pressure, t)
(which are the quantities of a single particle), returning aTuple
orSVector
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 withBoundaryModelDummyParticles
andAdamiPressureExtrapolation
. The keyword argumentacceleration
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 orColorfieldSurfaceNormal()
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.
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
— TypeStateEquationCole(; 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 of7
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 withSummationDensity
but allows unphysical rarefaction of the fluid.
TrixiParticles.StateEquationIdealGas
— TypeStateEquationIdealGas( ;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 ratiobackground_pressure=0.0
: Background pressure.clip_negative_pressure=false
: Negative pressure values are clipped to 0, which prevents spurious surface tension withSummationDensity
but allows unphysical rarefaction of the fluid.
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 DensityDiffusion
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.
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$):
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.DensityDiffusion
— TypeDensityDiffusion
An abstract supertype of all density diffusion formulations.
Currently, the following formulations are available:
Formulation | Suitable for Steady-State Simulations | Low Computational Cost |
---|---|---|
DensityDiffusionMolteniColagrossi | ❌ | ✅ |
DensityDiffusionFerrari | ❌ | ✅ |
DensityDiffusionAntuono | ✅ | ❌ |
See Density Diffusion for a comparison and more details.
TrixiParticles.DensityDiffusionAntuono
— TypeDensityDiffusionAntuono(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 DensityDiffusion
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 DensityDiffusion
for an overview and comparison of implemented density diffusion terms.
TrixiParticles.DensityDiffusionFerrari
— TypeDensityDiffusionFerrari()
A density diffusion term by Ferrari (2009).
The term $\psi_{ab}$ in the continuity equation in DensityDiffusion
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 DensityDiffusion
for an overview and comparison of implemented density diffusion terms.
TrixiParticles.DensityDiffusionMolteniColagrossi
— TypeDensityDiffusionMolteniColagrossi(; delta)
The commonly used density diffusion term by Molteni (2009).
The term $\psi_{ab}$ in the continuity equation in DensityDiffusion
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 DensityDiffusion
for an overview and comparison of implemented density diffusion terms.
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. 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.
The Particle Shifting Technique can be applied in form of the ParticleShiftingCallback
.
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.
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.
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
— Functiontensile_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.