2.3. Fluid dynamics formulation

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.

2.3.1. Equation of state

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.

2.3.2. Pressure acceleration

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.

2.3.3. Viscous acceleration

For the particle acceleration due to viscous forces, two approximations are available through the method's option: Morris and Cleary.

2.3.3.1. Morris approximation

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 .

2.3.3.2. Cleary approximation

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:

  • is the dynamic viscosity of the SPH element i.

  • is a factor given by Cleary Visc. Factor, which is a method parameter with a default value equal to 4.96333.

  • is a small value that prevents division by zero, which is computed as:

    (2–21)

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.