Weakly Compressible SPH

WeaklyCompressibleSPHSystem(initial_condition,
                            density_calculator, state_equation,
                            smoothing_kernel, smoothing_length;
                            viscosity=nothing, density_diffusion=nothing,
                            acceleration=ntuple(_ -> 0.0, NDIMS),
                            buffer_size=nothing,
                            correction=nothing, source_terms=nothing)

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. See ContinuityDensity and SummationDensity.
• state_equation: Equation of state for the system. See StateEquationCole.
• 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

• viscosity: Viscosity model for this system (default: no viscosity). See ArtificialViscosityMonaghan or ViscosityAdami.
• density_diffusion: Density diffusion terms for this system. See DensityDiffusion.
• acceleration: Acceleration vector for the system. (default: zero vector)
• buffer_size: Number of buffer particles. This is needed when simulating with OpenBoundarySPHSystem.
• 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)

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.

ArtificialViscosityMonaghan(; alpha, beta=0.0, epsilon=0.01)

Artificial viscosity by Monaghan ([16], [17]), given by

\[\Pi_{ab} = \begin{cases} -(\alpha c \mu_{ab} + \beta \mu_{ab}^2) / \bar{\rho}_{ab} & \text{if } v_{ab} \cdot r_{ab} < 0, \\ 0 & \text{otherwise} \end{cases}\]

with

\[\mu_{ab} = \frac{h v_{ab} \cdot r_{ab}}{\Vert r_{ab} \Vert^2 + \epsilon h^2},\]

where $\alpha, \beta, \epsilon$ are parameters, $c$ is the speed of sound, $h$ is the smoothing length, $r_{ab} = r_a - r_b$ is the difference of the coordinates of particles $a$ and $b$, $v_{ab} = v_a - v_b$ is the difference of their velocities, and $\bar{\rho}_{ab}$ is the arithmetic mean of their densities.

Note that $\alpha$ needs to adjusted for different resolutions to maintain a specific Reynolds Number. To do so, Monaghan (2005) defined an equivalent effective physical kinematic viscosity $\nu$ by

\[ \nu = \frac{\alpha h c }{2d + 4},\]

where $d$ is the dimension.

Keywords

• alpha: A value of 0.02 is usually used for most simulations. For a relation with the kinematic viscosity, see description above.
• beta=0.0: A value of 0.0 works well for most fluid simulations and simulations with shocks of moderate strength. In simulations where the Mach number can be very high, eg. astrophysical calculation, good results can be obtained by choosing a value of beta=2.0 and alpha=1.0.
• epsilon=0.01: Parameter to prevent singularities.

ViscosityAdami(; nu, epsilon=0.01)

Viscosity by Adami (2012). The viscous interaction is calculated with the shear force for incompressible flows given by

\[f_{ab} = \sum_w \bar{\eta}_{ab} \left( V_a^2 + V_b^2 \right) \frac{v_{ab}}{||r_{ab}||^2+\epsilon h_{ab}^2} \nabla W_{ab} \cdot r_{ab},\]

where $r_{ab} = r_a - r_b$ is the difference of the coordinates of particles $a$ and $b$, $v_{ab} = v_a - v_b$ is the difference of their velocities, $h$ is the smoothing length and $V$ is the particle volume. The parameter $\epsilon$ prevents singularities (see Ramachandran (2019)). The inter-particle-averaged shear stress is

\[ \bar{\eta}_{ab} =\frac{2 \eta_a \eta_b}{\eta_a + \eta_b},\]

where $\eta_a = \rho_a \nu_a$ with $\nu$ as the kinematic viscosity.

Keywords

• nu: Kinematic viscosity
• epsilon=0.01: Parameter to prevent singularities

ViscosityMorris(; nu, epsilon=0.01)

Viscosity by Morris (1997) also used by Fourtakas (2019).

To the force $f_{ab}$ between two particles $a$ and $b$ due to pressure gradients, an additional force term $\tilde{f}_{ab}$ is added with

\[\tilde{f}_{ab} = m_a m_b \frac{(\mu_a + \mu_b) r_{ab} \cdot \nabla W_{ab}}{\rho_a \rho_b (\Vert r_{ab} \Vert^2 + \epsilon h^2)} v_{ab},\]

where $\mu_a = \rho_a \nu$ and $\mu_b = \rho_b \nu$ denote the dynamic viscosity of particle $a$ and $b$ respectively, and $\nu$ is the kinematic viscosity.

Keywords

• nu: Kinematic viscosity
• epsilon=0.01: Parameter to prevent singularities

DensityDiffusion

An abstract supertype of all density diffusion formulations.

Currently, the following formulations are available:

See Density Diffusion for a comparison and more details.

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 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_{r_a} W(\Vert r_{ab} \Vert, h)\]

with

\[L_a := \left( -\sum_{b} V_b r_{ab} \otimes \nabla_{r_a} W(\Vert r_{ab} \Vert, h) \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.

DensityDiffusionFerrari()

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}{2h} \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$ is the smoothing length.

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

DensityDiffusionMolteniColagrossi(; 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.

AkinciFreeSurfaceCorrection(rho0)

Free surface correction according to Akinci et al. (2013). At a free surface, the mean density is typically lower than the reference density, resulting in reduced surface tension and viscosity forces. The free surface correction adjusts the viscosity, pressure, and surface tension forces near free surfaces to counter this effect. It's important to note that this correlation is unphysical and serves as an approximation. The computation time added by this method is about 2–3%.

Mathematically the idea is quite simple. If we have an SPH particle in the middle of a volume at rest, its density will be identical to the rest density $\rho_0$. If we now consider an SPH particle at a free surface at rest, it will have neighbors missing in the direction normal to the surface, which will result in a lower density. If we calculate the correction factor

\[k = \rho_0/\rho_\text{mean},\]

this value will be about ~1.5 for particles at the free surface and can then be used to increase the pressure and viscosity accordingly.

Arguments

• rho0: Rest density.

BlendedGradientCorrection()

Calculate a blended gradient to reduce the stability issues of the GradientCorrection as explained by Bonet (1999).

This calculates the following,

\[\tilde\nabla A_i = (1-\lambda) \nabla A_i + \lambda L_i \nabla A_i\]

with $0 \leq \lambda \leq 1$ being the blending factor.

Arguments

• blending_factor: Blending factor between corrected and regular SPH gradient.

GradientCorrection()

Compute the corrected gradient of particle interactions based on their relative positions (see Bonet, 1999).

Mathematical Details

Given the standard SPH representation, the gradient of a field $A$ at particle $a$ is given by

\[\nabla A_a = \sum_b m_b \frac{A_b - A_a}{\rho_b} \nabla_{r_a} W(\Vert r_a - r_b \Vert, h),\]

where $m_b$ is the mass of particle $b$ and $\rho_b$ is the density of particle $b$.

The gradient correction, as commonly proposed, involves multiplying this gradient with a correction matrix $L$:

\[\tilde{\nabla} A_a = \bm{L}_a \nabla A_a\]

The correction matrix $\bm{L}_a$ is computed based on the provided particle configuration, aiming to make the corrected gradient more accurate, especially near domain boundaries.

To satisfy

\[\sum_b V_b r_{ba} \otimes \tilde{\nabla}W_b(r_a) = \left( \sum_b V_b r_{ba} \otimes \nabla W_b(r_a) \right) \bm{L}_a^T = \bm{I}\]

the correction matrix $\bm{L}_a$ is evaluated explicitly as

\[\bm{L}_a = \left( \sum_b V_b \nabla W_b(r_{a}) \otimes r_{ba} \right)^{-1}.\]

KernelCorrection()

Kernel correction, as explained by Bonet (1999), uses Shepard interpolation to obtain a 0-th order accurate result, which was first proposed by Li et al. This can be further extended to obtain a kernel corrected gradient as shown by Basa et al. (2008).

The kernel correction coefficient is determined by

\[c(x) = \sum_{b=1} V_b W_b(x)\]

The gradient of corrected kernel is determined by

\[\nabla \tilde{W}_{b}(r) =\frac{\nabla W_{b}(r) - W_b(r) \gamma(r)}{\sum_{b=1} V_b W_b(r)} , \quad \text{where} \quad \gamma(r) = \frac{\sum_{b=1} V_b \nabla W_b(r)}{\sum_{b=1} V_b W_b(r)}.\]

This correction can be applied with SummationDensity and ContinuityDensity, which leads to an improvement, especially at free surfaces.

MixedKernelGradientCorrection()

Combines GradientCorrection and KernelCorrection, which results in a 1st-order-accurate SPH method (see Bonet, 1999).

Notes:

• Stability issues, especially when particles separate into small clusters.
• Doubles the computational effort.

ShepardKernelCorrection()

Kernel correction, as explained by Bonet (1999), uses Shepard interpolation to obtain a 0-th order accurate result, which was first proposed by Li et al. (1996).

The kernel correction coefficient is determined by

\[c(x) = \sum_{b=1} V_b W_b(x),\]

where $V_b = m_b / \rho_b$ is the volume of particle $b$.

This correction is applied with SummationDensity to correct the density and leads to an improvement, especially at free surfaces.

CohesionForceAkinci(surface_tension_coefficient=1.0)

This model only implements the cohesion force of the [25] surface tension model.

Keywords

• surface_tension_coefficient=1.0: Modifies the intensity of the surface tension-induced force, enabling the tuning of the fluid's surface tension properties within the simulation.

SurfaceTensionAkinci(surface_tension_coefficient=1.0)

Implements a model for surface tension and adhesion effects drawing upon the principles outlined by [25]. This model is instrumental in capturing the nuanced behaviors of fluid surfaces, such as droplet formation and the dynamics of merging or separation, by utilizing intra-particle forces.

Keywords

• surface_tension_coefficient=1.0: A parameter to adjust the magnitude of surface tension forces, facilitating the fine-tuning of how surface tension phenomena are represented in the simulation.