Fluid Models

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

Artificial viscosity by Monaghan ([9], [10]), 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

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.

ColorfieldSurfaceNormal()

Color field based computation of the interface normals.

CohesionForceAkinci(surface_tension_coefficient=1.0)

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