For the 2025 R2 release, the solution of IISPH (Implicit Incompressible SPH) and DFSPH (Divergence-Free SPH) were unified into one code of incompressible formulations implemented in Rocky software. DFSPH can be considered an extension of the IISPH formulation, because both enforce a constant density value through the solution of a Poisson pressure equation. However, DFSPH introduces an additional step in which the divergence-free condition is enforced on the velocity field. Although both constraints are equivalent at the continuum level (constant density and zero velocity divergence), enforcing just one of them in a discrete model does not guarantee that the other one will be automatically satisfied. DFSPH enforces the divergence-free condition by solving a second Poisson pressure equation. The overhead of this extra step is usually compensated by a reduction on the overall number of iterations during a simulation and an increased stability [16]

For the IISPH, the equations described in this section are valid until Equation 3–84. All the other ones are used for DFSPH calculations.

The pressure Poisson equation solved on both stages of the DFSPH algorithm can be written for a generic SPH element i as:

(3–74)

where is the pressure Laplacian, whereas is a source term with a form that will depend on the constraint being enforced (see Equation 3–75 and Equation 3–85 below). For the constant density constraint, the source term will be given by:

(3–75)

where is the fluid's rest density and is the predicted density value, obtained from the discretized form of the continuity equation:

(3–76)

The velocity is an intermediate velocity obtained after integrating the motion equation, including only the accelerations caused by all non-pressure forces (gravity, viscous forces, surface tension forces, etc.) acting on the element i:

(3–77)

The predicted density value will be, in general, different from the specified fluid density . Thus, the role of the pressure Poisson equation Equation 3–74 with the source term given by Equation 3–75 is to determine pressure values that introduce additional accelerations that cancel out the density deviation caused by the non-pressure forces.

The pressure Laplacian is approximated by considering the following expression based on the pressure accelerations [17]:

(3–78)

which, in turn, are expressed by using the symmetric approximation commonly used in SPH [18]:

(3–79)

As usual in the SPH literature, in the previous expressions j stands for the elements inside the kernel support of the generic element i, while is the gradient of the kernel function.

The substitution of Equation 3–78 and Equation 3–79 into Equation 3–74 produces a linear equation in which the only unknowns are the pressure values of all elements in the kernel support of element i. This equation can be written as:

(3–80)

When considering the previous equation for all the n SPH elements in a simulation, a system of n linear equations with n pressure unknowns is obtained. As suggested in [13], a relaxed Jacobi iterative method is used in Rocky for solving the pressure system of equations. At each step of the iterative process, the pressure value at an element i is updated according to the following expression:

(3–81)

where k denotes the iteration index, and is an under-relaxation factor introduced to ensure the stability of the process by reducing the amount of change between iterations. If necessary, the value of this factor can be altered on the UI at Pressure Under-relaxation Factor.

The convergence criterion of the pressure iterative process is based on the average value of the corrected density after k + 1 iterations. The normalized value of the deviation of that average density in relation to is estimated as:

(3–82)

The pressure iterations will continue until , where the value of the tolerance can be altered on the UI at Density Relative Tolerance. If for some reason the convergence criterion above is not satisfied after iterations, the pressure iterative process will be aborted. Whenever that happens, Rocky will emit a warning alert on the Simulation Log, however, the simulation will continue being executed. Users can alter the value of on the UI at Maximum Number of Pressure Iterations.

After the iterative process has converged, a new set of velocity values is obtained by applying the pressure accelerations calculated with the pressure values coming from the solution of the linear system that enforces constant density:

(3–83)

Following the DFSPH algorithm proposed in [17], at this point the positions of all elements are updated with those velocities:

(3–84)

However, the divergence of velocity field may not be zero, as required by the continuity equation. In order to enforce that condition, an additional system of pressure equations is constructed based on Equation 3–74, as stated earlier. The source term for the new system of pressure equations is given by:

(3–85)

By solving this second system of equations, a new set of pressure values is found. This produces additional accelerations that removes the divergence error on the velocity field . The solution process and the convergence criterion for the second system of equations is the same described above for the system that enforces constant density. In fact, the same parameter values that control the iterative process are used on both systems of equations.

When the second system of equations has converged, the pressure accelerations obtained with the new set of pressure values are used to determine the velocity values at the end of the time step:

(3–86)

Additionally, the pressure values obtained enforcing zero velocity divergence are summed to the ones enforcing constant density, and the sum is considered to be the pressure . With the pressure, velocity and position of the SPH elements updated, the simulation proceeds to the next discrete time level.

 

3.14.1. Relaxed Incompressibility Constraint (BETA)

---

Important:   Note that beta features have not been fully tested and validated. Ansys, Inc. makes no commitment to resolve defects reported against these prototype features. However, your feedback will help us improve the overall quality of the product. We will not guarantee that the projects using this beta feature will run successfully when the feature is finally released so you may, therefore, need to modify the projects.

---

The Smoothed Particle Hydrodynamics method was initially developed for compressible flows. Hence, special treatment for solving the incompressible flow equations is necessary to mitigate or eliminate a few issues. On SPH incompressible solvers, one possible issue that can arise is the unrealistic pressure fluctuations.

For IISPH (Implicit Incompressible SPH), there is no source term to deal with this pressure problem. Some examples of the cases that we can observe pressure fluctuations are on discharge of fluids and open channel with obstacles (dam break case). However, other physical problems can show those high pressure fluctuations. For the Divergence-Free SPH (DFSPH), a special condition is implemented that mitigates unphysical pressure fluctuations, see Equation 3–85 on IISPH - Implicit Incompressible SPH and DFSPH - Divergence Free SPH (Beta).

To address this limitation the Relaxed Incompressibility Constraint is implemented for IISPH. The IISPH solver currently solves the following PPE (pressure Poisson Equation) to get the pressure field (see Equation 3–74 and ???)

(3–87)

The right hand side of the equation is known as density invariance condition or constant density condition. This is a condition that preserves the volume/density of the bulk of the fluid, but results in a pressure field exhibiting oscillations of larger amplitude.

On the other hand there is a condition named divergence free condition shown in the right hand side of the following PPE:

(3–88)

The Relaxed Incompressibility Constraint combines both conditions, Equation 3–87 and Equation 3–88, in the following manner

(3–89)

The condition shown in Equation 3–88 has a well behaved pressure field, less noisy but does not preserve well the volume of the bulk of the fluid and tends to shrink it. One can observe this this effect from doing a simple test, setting the parameter equals to 0.

In this case, only one PPE is solved in which the source term is a combination of the divergence free condition and a fraction of the constant density condition. In this condition we have an extra parameter which is a relaxation factor that must be tuned according to SPH element size and fluid properties. After turning on Relaxed Incompressibility Constraint, the relaxation factor is the only factor which the user can define.