Fluid Models

Currently available fluid methods are the weakly compressible SPH method and the entropically damped artificial compressibility for SPH. This page lists models and techniques that apply to both of these methods.

Viscosity

Viscosity is a critical physical property governing momentum diffusion within a fluid. In the context of SPH, viscosity determines how rapidly velocity gradients are smoothed out, influencing key flow characteristics such as boundary layer formation, vorticity diffusion, and dissipation of kinetic energy. It also helps determine whether a flow is laminar or turbulent under a given set of conditions.

Implementing viscosity correctly in SPH is essential for producing physically accurate results, and different methods exist to capture both numerical stabilization and true viscous effects.

Artificial (numerical) viscosity

Artificial (numerical) viscosity is a technique used to stabilize simulations of inviscid flows, which would otherwise show unphysical particle movement due to numerical instability. To achieve this, a dissipative term is added to the momentum equations in a way that it does not significantly alter the physical behavior of the flow. This approach is especially useful in simulations such as high-speed flows with strong shocks or astrophysical scenarios, where other approaches are insufficient to stabilize the simulation.

Physical (real) viscosity

Physical viscosity is essential for accurately modeling the true viscous stresses within a fluid. It ensures that simulations align with a target Reynolds number or adhere to experimentally measured fluid properties. This is achieved by incorporating forces that replicate the viscous stress term found in the Navier–Stokes equations. As a result, the method is particularly effective for simulating low-speed, incompressible, or weakly compressible flows, where it is crucial to capture the actual behavior of the fluid.

Model comparison

ArtificialViscosityMonaghan

ArtificialViscosityMonaghan by Monaghan ([9], [10]) should be mainly used for inviscid flows (Euler), artificial stabilization or shock-capturing, for which Monaghan [10] originally designed this term to provide smoothing across shocks, intentionally overestimating the physical viscosity. The implementation includes a dissipation term that becomes more significant as particles approach one another. This helps suppress tensile instabilities, which can lead to particle clumping and effectively smooths out high-frequency pressure fluctuations. This increase in dissipation is triggered by the relative motion between particles: as particles come closer and compress the local flow, the artificial viscosity term becomes stronger to damp out rapid changes and prevent unphysical clustering. This ensures that while the simulation remains stable in challenging flow regimes with large density or pressure variations, the physical behavior is not overly altered. Several extensions have been proposed to limit the dissipation effect for example by Balsara ([11]) or Morris ([12]).

Mathematical Formulation

The force exerted by particle b on particle a due to artificial viscosity is given by:

\[\mathbf{F}_{ab}^{\text{AV}} = - m_a m_b \Pi_{ab} \nabla W_{ab}\]

where:

• \Pi_{ab} is the artificial viscosity term defined as:

  ```math \Pi{ab} = \begin{cases} -\frac{\alpha c \mu{ab} + \beta \mu{ab}^2}{\bar{\rho}{ab}} & \text{if } \mathbf{v}{ab} \cdot \mathbf{r}{ab} < 0, \ 0 & \text{otherwise} \end{cases} ```

• $\alpha$ and $\beta$ are viscosity parameters,

• $c$ is the local speed of sound,

• $\bar{\rho}_{ab}$ is the arithmetic mean of the densities of particles $a$ and $b$.

The term $\mu_{ab}$ is defined as:

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

with:

• $h$ being the smoothing length,
• $\epsilon$ a small parameter to prevent singularities,
• $r_{ab} = r_a - r_b$ representing the difference of the coordinate vectors,
• $v_{ab} = v_a - v_b$ representing the relative velocity between particles.

Resolution Dependency and Effective Viscosity

To ensure that the simulation maintains a consistent Reynolds number when the resolution changes, the parameter $\alpha$ must be adjusted accordingly. Monaghan (2005) introduced an effective physical kinematic viscosity $\nu$ defined as:

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

where $d$ is the number of spatial dimensions. This relation allows the calibration of $\alpha$ to achieve the desired viscous behavior as the resolution or simulation conditions vary.

ViscosityMorris

ViscosityMorris is ideal for moderate to low Mach number flows where accurately modeling physical viscous behavior is essential. Developed by Morris (1997) and later applied by Fourtakas (2019), this method directly simulates the viscous stresses found in fluids rather than relying on artificial viscosity. By approximating momentum diffusion based on local fluid properties, the method captures the actual viscous forces without excessive damping. This results in a more realistic representation of flow dynamics in weakly compressible scenarios.

Mathematical Formulation

An additional force term $\tilde{f}_{ab}$ is introduced to the pressure gradient force $f_{ab}$ between particles $a$ and $b$:

\[\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$ represent the dynamic viscosities of particles $a$and $b$ (with $\nu$ being the kinematic viscosity),
• $r_{ab} = r_a - r_b$ represents the difference of the coordinate vectors,
• $v_{ab} = v_a - v_b$ represents the relative velocity between particles.
• $h$ is the smoothing length,
• $\nabla W_{ab}$ is the gradient of the smoothing kernel,
• $\epsilon$ is a small parameter to prevent singularities.

ViscosityAdami

ViscosityAdami, introduced by Adami (2012), is optimized for incompressible or weakly compressible flows where precise modeling of shear stress is critical. It enhances boundary layer representation by better resolving shear gradients, increasing dissipation in regions with steep velocity differences (e.g., near solid boundaries) while minimizing compressibility effects. This results in accurate laminar flow simulations and accurate physical shear stresses.

Mathematical Formulation

The viscous interaction is modeled through a shear force for incompressible flows:

\[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 coordinate vectors,
• $v_{ab} = v_a - v_b$ is their relative velocity,
• $V_a = m_a / \rho_a$ and $V_b = m_b / \rho_b$ are the particle volumes,
• $h_{ab}$ is the smoothing length,
• $\nabla W_{ab}$ is the gradient of the smoothing kernel,
• $\epsilon$ is a small parameter that prevents singularities (see Ramachandran (2019)).

The inter-particle-averaged shear stress is defined as:

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

with the dynamic viscosity of each particle given by $\eta_a = \rho_a \nu_a$, where $\nu_a$ is the kinematic viscosity.

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

ViscosityAdami(; nu, epsilon=0.01)

ViscosityMorris(; nu, epsilon=0.01)

Corrections

AkinciFreeSurfaceCorrection(rho0)

BlendedGradientCorrection()

GradientCorrection()

KernelCorrection()

MixedKernelGradientCorrection()

ShepardKernelCorrection()

Surface Normals

Overview of surface normal calculation in SPH

Surface normals are essential for modeling surface tension as they provide the directionality of forces acting at the fluid interface. They are calculated based on the particle properties and their spatial distribution.

Color field and gradient-based surface normals

The surface normal at a particle is derived from the color field, a scalar field assigned to particles to distinguish between different fluid phases or between fluid and air. The color field gradients point towards the interface, and the normalized gradient defines the surface normal direction.

The simplest SPH formulation for a surface normal, $n_a$ is given as

\[n_a = \sum_b m_b \frac{c_b}{\rho_b} \nabla_a W_{ab},\]

where:

• $c_b$ is the color field value for particle $b$,
• $m_b$ is the mass of particle $b$,
• $\rho_b$ is the density of particle $b$,
• $\nabla_a W_{ab}$ is the gradient of the smoothing kernel $W_{ab}$ with respect to particle $a$.

Normalization of surface normals

The calculated normals are normalized to unit vectors:

\[\hat{n}_a = \frac{n_a}{\Vert n_a \Vert}.\]

Normalization ensures that the magnitude of the normals does not bias the curvature calculations or the resulting surface tension forces.

Handling noise and errors in normal calculation

In regions distant from the interface, the calculated normals may be small or inaccurate due to the smoothing kernel's support radius. To mitigate this:

1. Normals below a threshold are excluded from further calculations.
2. Curvature calculations use a corrected formulation to reduce errors near interface fringes.

ColorfieldSurfaceNormal(; boundary_contact_threshold=0.1, interface_threshold=0.01,
                          ideal_density_threshold=0.0)

Surface Tension

Surface tension is a key phenomenon in fluid dynamics, influencing the behavior of droplets, bubbles, and fluid interfaces. In SPH, surface tension is modeled as forces arising due to surface curvature and relative particle movement, ensuring realistic simulation of capillary effects, droplet coalescence, and fragmentation.

The surface tension coefficient $\sigma$ is a physical parameter that quantifies the energy required to increase the surface area of a fluid by a unit amount. A higher value of $\sigma$ indicates that the fluid resists changes to its surface area more strongly, causing droplets or bubbles to assume shapes (often spherical) that minimize their surface. In practice, $\sigma$ can be measured experimentally through techniques such as the pendant drop method, the Wilhelmy plate method, or the du Noüy ring method, each of which relates a measurable force or change in shape to the fluid’s surface tension. For pure substances, tabulated reference values of $\sigma$ at given temperatures are commonly used, while for mixtures or complex fluids, direct experimental measurements or values can be estimated from empirical equation (see Poling or Lange). In the following table some values are shown for reference. The values marked with a '~' are complex mixtures that are estimated by an empirical equation (see Poling).

Akinci-based intra-particle force surface tension and wall adhesion model

The Akinci model divides surface tension into distinct force components:

Cohesion force

The cohesion force captures the attraction between particles at the fluid interface, creating the effect of surface tension. It is defined by the distance between particles and the support radius $h_c$, using a kernel-based formulation.

Key features:

• Particles within half the support radius experience a repulsive force to prevent clustering.
• Particles beyond half the radius but within the support radius experience an attractive force to simulate cohesion.

Mathematically:

\[F_{\text{cohesion}} = -\sigma m_b C(r) \frac{r}{\Vert r \Vert},\]

where $C(r)$, the cohesion kernel, is defined as:

\[C(r)=\frac{32}{\pi h_c^9} \begin{cases} (h_c-r)^3 r^3, & \text{if } 2r > h_c, \\ 2(h_c-r)^3 r^3 - \frac{h^6}{64}, & \text{if } r > 0 \text{ and } 2r \leq h_c, \\ 0, & \text{otherwise.} \end{cases}\]

Surface area minimization force

The surface area minimization force models the curvature reduction effects, aligning particle motion to reduce the interface's total area. It acts based on the difference in surface normals:

\[F_{\text{curvature}} = -\sigma (n_a - n_b),\]

where $n_a$ and $n_b$ are the surface normals of the interacting particles.

Wall adhesion force

This force models the interaction between fluid and solid boundaries, simulating adhesion effects at walls. It uses a custom kernel with a peak at 0.75 times the support radius:

\[F_{\text{adhesion}} = -\beta m_b A(r) \frac{r}{\Vert r \Vert},\]

where $A(r)$ is the adhesion kernel:

\[A(r) = \frac{0.007}{h_c^{3.25}} \begin{cases} \sqrt[4]{-\frac{4r^2}{h_c} + 6r - 2h_c}, & \text{if } 2r > h_c \text{ and } r \leq h_c, \\ 0, & \text{otherwise.} \end{cases}\]

Morris surface tension model

The method described by Morris estimates curvature by combining particle color gradients (see surface_normal) and smoothing functions to derive surface normals. The computed curvature is then used to determine forces acting perpendicular to the interface. While this method provides accurate surface tension forces, it does not explicitly conserve momentum.

In the Morris model, surface tension is computed based on local interface curvature $\kappa$ and the unit surface normal $\hat{n}.$ By estimating $\hat{n}$ and $\kappa$ at each particle near the interface, the surface tension force for particle a can be written as:

\[F_{\text{surface tension}} = - \sigma \frac{\kappa_a}{\rho_a}\hat{n}_a\]

This formulation focuses directly on geometric properties of the interface, making it relatively straightforward to implement when a reliable interface detection (e.g., a color function) is available. However, accurately estimating $\kappa$ and $n$ may require fine resolutions.

Morris-based momentum-conserving surface tension model

In addition to the simpler curvature-based formulation, Morris introduced a momentum-conserving approach. This method treats surface tension forces as arising from the divergence of a stress tensor, ensuring exact conservation of linear momentum and offering more robust behavior for high-resolution or long-duration simulations where accumulated numerical error can be significant.

Stress tensor formulation

The surface tension force can be seen as a divergence of a stress tensor $S$

\[F_{\text{surface tension}} = \nabla \cdot S,\]

with $S$ defined as

\[S = \sigma \delta_s (I - \hat{n} \otimes \hat{n}),\]

with:

• $\delta_s$: Surface delta function,
• $\hat{n}$: Unit normal vector,
• $I$: Identity matrix.

This divergence can be computed numerically in the SPH framework as

\[\sum_b \frac{m_b}{\rho_a \rho_b} (S_a + S_b) \nabla W_{ab} \]

Advantages and limitations

While momentum conservation makes this model attractive, it requires additional computational effort and stabilization techniques to address instabilities in high-density regions.

API

CohesionForceAkinci(surface_tension_coefficient=1.0)

SurfaceTensionAkinci(surface_tension_coefficient=1.0)

SurfaceTensionMomentumMorris(surface_tension_coefficient=1.0)

SurfaceTensionMorris(surface_tension_coefficient=1.0)