Currently, the software uses the Weakly Compressible SPH (WCSPH), and Incompressible SPH (IISPH). The WCSPH is suitable for compressible flows with Mach number < 0.1, thus, effectively simulating incompressible fluids.
The continuity equation in the WCSPH formulation used by the software is given by Eq.
Equation 3–10, and describes how the density of an SPH element at
position evolves with time:
(3–10) |
where is the gradient of
evaluated with respect to the position
.
The equations of motion for Newtonian fluid elements are derived from the Navier-Stokes equations. As shown in Eq. Equation 3–11, the equation of motion is integrated explicitly to give the velocity variation over the time.
(3–11) |
where:
is the velocity of the SPH element.
is the SPH element acceleration due to pressure forces.
is the SPH element acceleration due to viscous forces.
is the SPH element acceleration due to body forces.
The pressure and viscous forces are defined by the SPH-discretized continuity equation, the Navier-Stokes equations, and the equation of state, which will be explained further in the sections that follow.
The continuity equation shows that the time rate of density change of an SPH element is closely related to the relative velocities between this element and all other SPH elements in the support domain. The gradient of the smoothing function determines the contribution of these relative velocities. Hence, while the mass of each SPH element is fixed, the density associated with that SPH element changes as the local volume associated with that mass changes. This assumption makes SPH simulations naturally compressible. Since the focus of this version involves incompressible fluid flow, an artificial equation of state can be introduced to approximate incompressibility by ensuring that the Mach number is sufficiently low.
Equations of state for WCSPH can be freely designed as long as it assures that the Mach number in the simulation is smaller than 0.1. the software, like in most of the WCSPH formulations, uses an equation of state proposed for water to compute the relation between pressure and density:
(3–12) |
where:
is the pressure of the SPH element a.
is the background pressure. It defines the absolute pressure level around which small pressure fluctuations occur. It is used internally in the equation of state to stabilize simulations by reducing large density variations.
is the local fluid density.
is initial fluid density or rest density, defined by Density, which is a property of the fluid material .
is the material speed of sound defined by Sound Speed, which is a parameter.
is a dimensionless value defining the material stiffness defined by Pressure Degree, which is a method parameter.
You must ensure that the material speed of sound, , is at least 10x larger than maximum simulated fluid velocity expected
during the simulation.
When the software is used with very small or high values of the Sound Speed , it can result in numerical issues. In the software, the current default value for the sound speed is set to 100 m/s, which often leads to unnecessarily small time steps in most cases. This becomes especially apparent when dealing with minimal fluid heights, as the density variation can be of similar magnitude as the numerical errors. To mitigate these issues, it is recommended to estimate the sound speed based on the specific problem at hand.
The SPH element acceleration due to pressure, , is calculated using the [1] formula:
(3–13) |
where:
and
are the pressure values corresponding to the positions
and
.
is the gradient of
evaluated with respect to the position
.
In Equation 3–13, is the tensile instability correction that prevents SPH elements from
clumping together for some flows. The tensile instability can be removed by adding certain
artificial pressure to the momentum equation, as proposed by Monaghan (2000) [1].
Important: The tensile correction term is considered only for WCSPH, not the incompressible solvers. And even for WCSPH it is hardly needed, as other alternatives, such as position correction, are usually preferred
The repulsive force must increase as the separation between the two SPH elements a and b decreases. That way, it is natural to express this repulsive force in terms of the kernel, as:
(3–14) |
where:
is the distance between the SPH elements a and b.
denotes the average SPH element spacing in the neighborhood of SPH element a.
The tensile instability correction can then be written as:
(3–15) |
where n>0 is the Stability
Degree, which is a method parameter, and the factor depends on the pressure and density. For the default value of Stability Degree, n=4, and with
h equal to
, the repulsive force increases by a factor of
as d decreases from
to zero. Additionally, as the ratio
decreases rapidly in the domain
, only the nearest neighbors are significantly affected by the artificial
pressure.
The factor is related to the pressure, and can be expressed as:
(3–16) |
where is computed from:
(3–17) |
The rule for is obtained by replacing a with
b in Eq. Equation 3–17 and the
coefficient
is given as:
(3–18) |
In the SPH method,
Kpos
and
Kneg
are exposed as user inputs as Stability Negative Factor and
Stability Positive Factor, respectively. However, the
tensile instability correction was found unnecessary for Rocky and therefore the default
values for Stability Negative Factor and Stability Positive Factor are zero, therefore making by default.
For the SPH element acceleration due to viscous forces, two approximations are available through the option: Morris and Cleary.
The Morris approximation is suitable for slow laminar flows. The viscous term is approximated according to the paper by Morris et al. (1997)[3] by means of the expression:
(3–19) |
where:
and
are the dynamic viscosity values of the fluid corresponding to positions
and
.
is the gradient of
evaluated with respect to the position
.
The Cleary approximation is suitable for fast turbulent flows. According to the work of Cleary (1998)[4], as shown in Eq. Equation 3–20:
(3–20) |
where:
where:
is given by Minimum Distance Factor, which is a method parameter with a default value of 0.01.
is the initial SPH elements spacing.
For viscous force solution, there are two methods available through Viscous Force Solution on the Advanced tab of SPH: Explicit and Implicit. For the Explicit solution, the viscosity is solved explicitly, so that the time step of the simulation is calculated by the procedures given on SPH time step calculation. However, for the Implicit solution, after the pressure forces are solved and the velocity is updated, the viscous force is calculated using an implicit discretization:
(3–22) |
where:
is the velocity calculated before the viscous force is
applied.