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 (Monaghan 1992, Monaghan 1989), 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 (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.

References

• Joseph J. Monaghan. "Smoothed Particle Hydrodynamics". In: Annual Review of Astronomy and Astrophysics 30.1 (1992), pages 543-574. doi: 10.1146/ANNUREV.AA.30.090192.002551
• Joseph J. Monaghan. "Smoothed Particle Hydrodynamics". In: Reports on Progress in Physics (2005), pages 1703-1759. doi: 10.1088/0034-4885/68/8/r01
• Joseph J. Monaghan. "On the Problem of Penetration in Particle Methods". In: Journal of Computational Physics 82.1, pages 1–15. doi: 10.1016/0021-9991(89)90032-6

ViscosityAdami(; nu, epsilon=0.01)

Viscosity by Adami (Adami et al. 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 et al. 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

References

• S. Adami et al. "A generalized wall boundary condition for smoothed particle hydrodynamics". In: Journal of Computational Physics 231 (2012), pages 7057-7075. doi: 10.1016/j.jcp.2012.05.005
• P. Ramachandran et al. "Entropically damped artificial compressibility for SPH". In: Journal of Computers and Fluids 179 (2019), pages 579-594. doi: 10.1016/j.compfluid.2018.11.023

ViscosityMorris(; nu, epsilon=0.01)

Viscosity by Morris et al. (1997).

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

References

• Joseph P. Morris, Patrick J. Fox, Yi Zhu. "Modeling Low Reynolds Number Incompressible Flows Using SPH". In: Journal of Computational Physics, Volume 136, Issue 1 (1997), pages 214–226. doi: doi.org/10.1006/jcph.1997.5776
• Georgios Fourtakas, Jose M. Dominguez, Renato Vacondio, Benedict D. Rogers. "Local uniform stencil (LUST) boundary condition for arbitrary 3-D boundaries in parallel smoothed particle hydrodynamics (SPH) models". In: Computers & Fluids, Volume 190 (2019), pages 346–361. doi: 10.1016/j.compfluid.2019.06.009

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 et al. (2010), also referred to as δ-SPH. The density diffusion term by Molteni & Colagrossi (2009) is extended by a second term, which is nicely written down by Antuono et al. (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.

References

• M. Antuono, A. Colagrossi, S. Marrone, D. Molteni. "Free-Surface Flows Solved by Means of SPH Schemes with Numerical Diffusive Terms." In: Computer Physics Communications 181.3 (2010), pages 532–549. doi: 10.1016/j.cpc.2009.11.002
• M. Antuono, A. Colagrossi, S. Marrone. "Numerical Diffusive Terms in Weakly-Compressible SPH Schemes." In: Computer Physics Communications 183.12 (2012), pages 2570–2580. doi: 10.1016/j.cpc.2012.07.006
• Diego Molteni, Andrea Colagrossi. "A Simple Procedure to Improve the Pressure Evaluation in Hydrodynamic Context Using the SPH." In: Computer Physics Communications 180.6 (2009), pages 861–872. doi: 10.1016/j.cpc.2008.12.004

DensityDiffusionFerrari()

A density diffusion term by Ferrari et al. (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.

References

• Angela Ferrari, Michael Dumbser, Eleuterio F. Toro, Aronne Armanini. "A New 3D Parallel SPH Scheme for Free Surface Flows." In: Computers & Fluids 38.6 (2009), pages 1203–1217. doi: 10.1016/j.compfluid.2008.11.012.

DensityDiffusionMolteniColagrossi(; delta)

The commonly used density diffusion term by Molteni & Colagrossi (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.

References

• Diego Molteni, Andrea Colagrossi. "A Simple Procedure to Improve the Pressure Evaluation in Hydrodynamic Context Using the SPH." In: Computer Physics Communications 180.6 (2009), pages 861–872. doi: 10.1016/j.cpc.2008.12.004

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.

References

• Akinci, N., Akinci, G., & Teschner, M. (2013). "Versatile Surface Tension and Adhesion for SPH Fluids". ACM Transactions on Graphics (TOG), 32(6), 182. doi: 10.1145/2508363.2508405

BlendedGradientCorrection()

Calculate a blended gradient to reduce the stability issues of the GradientCorrection.

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.

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}.\]

• Better stability with smoother smoothing Kernels with larger support, e.g. SchoenbergQuinticSplineKernel or WendlandC6Kernel.
• Set dt_max =< 1e-3 for stability.

References

• J. Bonet, T.-S.L. Lok. "Variational and momentum preservation aspects of Smooth Particle Hydrodynamic formulations". In: Computer Methods in Applied Mechanics and Engineering 180 (1999), pages 97–115. doi: 10.1016/S0045-7825(99)00051-1
• Mihai Basa, Nathan Quinlan, Martin Lastiwka. "Robustness and accuracy of SPH formulations for viscous flow". In: International Journal for Numerical Methods in Fluids 60 (2009), pages 1127–1148. doi: 10.1002/fld.1927

KernelCorrection()

Kernel correction 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.

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.

References

• J. Bonet, T.-S.L. Lok. "Variational and momentum preservation aspects of Smooth Particle Hydrodynamic formulations". In: Computer Methods in Applied Mechanics and Engineering 180 (1999), pages 97-115. doi: 10.1016/S0045-7825(99)00051-1
• Mihai Basa, Nathan Quinlan, Martin Lastiwka. "Robustness and accuracy of SPH formulations for viscous flow". In: International Journal for Numerical Methods in Fluids 60 (2009), pages 1127–1148. doi: 10.1002/fld.1927
• Shaofan Li, Wing Kam Liu. "Moving least-square reproducing kernel method Part II: Fourier analysis". In: Computer Methods in Applied Mechanics and Engineering 139 (1996), pages 159-193. doi:10.1016/S0045-7825(96)01082-1

MixedKernelGradientCorrection()

Combines GradientCorrection and KernelCorrection, which results in a 1st-order-accurate SPH method.

Notes:

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

References

• J. Bonet, T.-S.L. Lok. "Variational and momentum preservation aspects of Smooth Particle Hydrodynamic formulations". In: Computer Methods in Applied Mechanics and Engineering 180 (1999), pages 97–115. doi: 10.1016/S0045-7825(99)00051-1
• Mihai Basa, Nathan Quinlan, Martin Lastiwka. "Robustness and accuracy of SPH formulations for viscous flow". In: International Journal for Numerical Methods in Fluids 60 (2009), pages 1127–1148. doi: 10.1002/fld.1927

ShepardKernelCorrection()

Kernel correction uses Shepard interpolation to obtain a 0-th order accurate result, which was first proposed by Li et al.

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.

References

• J. Bonet, T.-S.L. Lok. "Variational and momentum preservation aspects of Smooth Particle Hydrodynamic formulations". In: Computer Methods in Applied Mechanics and Engineering 180 (1999), pages 97-115. doi: 10.1016/S0045-7825(99)00051-1
• Mihai Basa, Nathan Quinlan, Martin Lastiwka. "Robustness and accuracy of SPH formulations for viscous flow". In: International Journal for Numerical Methods in Fluids 60 (2009), pages 1127–1148. doi: 10.1002/fld.1927
• Shaofan Li, Wing Kam Liu. "Moving least-square reproducing kernel method Part II: Fourier analysis". In: Computer Methods in Applied Mechanics and Engineering 139 (1996), pages 159–193. doi:10.1016/S0045-7825(96)01082-1

CohesionForceAkinci(surface_tension_coefficient=1.0)

This model only implements the cohesion force of the Akinci 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.

References

• Nadir Akinci, Gizem Akinci, Matthias Teschner. "Versatile Surface Tension and Adhesion for SPH Fluids". In: ACM Transactions on Graphics 32.6 (2013). doi: 10.1145/2508363.2508395

SurfaceTensionAkinci(surface_tension_coefficient=1.0)

Implements a model for surface tension and adhesion effects drawing upon the principles outlined by Akinci et al. 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.

References

• Nadir Akinci, Gizem Akinci, Matthias Teschner. "Versatile Surface Tension and Adhesion for SPH Fluids". In: ACM Transactions on Graphics 32.6 (2013). doi: 10.1145/2508363.2508395