Currently, Rocky uses the Weakly Compressible SPH (WCSPH) approach, which is suitable for compressible flows with Mach number < 0.1, thus, effectively simulating incompressible fluids.
The continuity equation in the WCSPH formulation used by Rocky is given by Eq. Equation 2–10, and describes how the density of an SPH element at
position evolves with time:
(2–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 2–11, the equation of motion is integrated explicitly to give the velocity variation over the time.
(2–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 particle 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. Rocky, like in most of the WCSPH formulations, uses an equation of state proposed for water to compute the relation between pressure and density:
(2–12) |
where:
is the pressure of the particle a.
is the background pressure.
is the local fluid density.
is initial fluid density, defined by Density, which is a method parameter.
is the material speed of sound defined by Sound Speed, which is a method parameter.
is a dimensionless value defining the material stiffness defined by Pressure Deg., which is a method parameter.
You must ensure that the material speed of sound, , is at least 10x larger than maximum simulated fluid velocity.
When the software is used with very small or high values of the Sound Speed , it can result in numerical issues. In the Rocky 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. One approach is to calculate the sound speed using
the Hydrostatic Pressure as a reference, employing the
formula .
The default value for is 7. The reference density is generally chosen as the initial density
associated to the SPH elements, which is based on the density of the fluid in the
fluid-solid mixture.
Particle acceleration due to pressure, , is calculated using the [1] formula:
(2–13) |
where:
and
are the pressure values corresponding to the positions
and
.
is the gradient of
evaluated with respect to the position
.
In Equation 2–13, is the tensile instability correction that prevents particles 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].
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:
(2–14) |
where:
is the distance between the SPH elements a and b.
denotes the average SPH element spacing in the neighborhood of particle a.
The tensile instability correction can then be written as:
(2–15) |
where n>0 is the Stab. Degree,
which is a method parameter, and the factor depends on the pressure and density. For the default value of Stab. 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:
(2–16) |
where is computed from:
(2–17) |
The rule for is obtained by replacing a with
b in Eq. Equation 2–17 and the
coefficient
is given as:
(2–18) |
In the SPH method,
Kpos
and
Kneg
are exposed as user inputs as Stab. Neg. Factor and
Stab. Pos. Factor, respectively. However, the tensile
instability correction was found unnecessary for Rocky and therefore the default values for
Stab. Neg. Factor and are zero, therefore making by default.
For the particle acceleration due to viscous forces, two approximations are available through the method's 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:
(2–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 and Ha (2011)[4], as shown in Eq. Equation 2–20:
(2–20) |
where:
where:
is given by Min. Dist. Factor, which is a method parameter with a default value of 0.01.
is the initial SPH elements spacing.