A unified transport-velocity formulation for SPH simulation of cohesive granular materials

In the Lagrangian framework, the governing equations of fluids include the conservation of mass and momentum. When materials are incompressible and viscosity is neglected, these equations are expressed as (1) (2)where is the density, is velocity, is the pressure, and is the body force. The continuity Eq. (1) is the same for other materials, i.e., elastic solids and granular materials.

To simulate incompressible flow under the weakly compressible approximation (Monaghan, 1994, Morris et al., 1997), an artificial isothermal equation of state is employed to complete the formulation of Eq. (2) (3)where is the initial density and is the sound speed. Generally, in fluid simulations, when an artificial sound speed of is applied (where represents the anticipated maximum flow velocity), the density fluctuates by approximately 1% (Morris et al., 1997).

The momentum equation for elastic solids is given by (4)where is the Cauchy stress tensor. The total stress tensor can be decomposed into the sum of hydrostatic pressure and shear stress, as illustrated below. (5)where is the identity tensor, is the pressure and can be solved by Eq. (3) under the weakly compressible assumption (Gray et al., 2001, Zhang et al., 2017a), and is the shear stress tensor. The sound speed for solids is defined as (6)where is the Young’s modulus, and is the Poisson’s ratio. Shear stress is derived from the integral of the shear stress rate . (7)The shear stress rate is defined by the constitutive relation of the material. For linear elastic materials, the constitutive relation is given by (8)where is the shear modulus, and is the shear strain rate, with being the dimension of the space. The strain rate tensor is defined as (9)where is the velocity gradient tensor, and denotes the transpose of a tensor.

The elastic-perfectly plastic Drucker–Prager constitutive model is employed to describe the elastoplastic behavior of granular materials (Drucker and Prager, 1952, Borja, 2013). Although some rate-dependent constitutive models, such as the -rheological constitutive model (Jop et al., 2006), are also commonly used to describe the motion of granular materials (Lemiale et al., 2011, Hurley and Andrade, 2017, Madraki et al., 2017), SPH studies often assume cohesionless materials (Yang et al., 2021, Zhu et al., 2022) when employing rheological constitutive models to facilitate comparison with experimental results. However, cohesionless granular materials do not exhibit tensile instability, and therefore the UTVF scheme proposed in this study is not required in such cases. Accordingly, this study adopts the Drucker–Prager constitutive model, which is widely used in SPH simulations (Bui et al., 2008, Nguyen et al., 2017, Bui and Nguyen, 2021, Feng et al., 2021, Zhang et al., 2024c). Comparisons are conducted with experimental results under cohesionless conditions and with previous numerical results, to validate the constitutive model is correctly implemented in the SPH code. When cohesion is considered, we demonstrate that tensile instability is completely eliminated by the UTVF scheme, and compare the results with those obtained using the artificial stress method (Bui et al., 2008) and the particle shifting scheme (Lallemand et al., 2025).

The momentum equation for granular materials can be represented by Eq. (4). Unlike elastic solids, where stress tensor is decomposed into pressure and the shear stress tensor , with pressure estimated with equation of state, the total stress tensor for plastic solids is calculated directly based on constitutive relation, as pressure is also linked to the yield behavior in the Drucker–Prager model.

For non-associate flow rule, the stress rate is given by Borja (2013) (10) with the following yield criterion and plastic potential function . (11) (12)where the subscript indicates the time step . is the deviatoric stress tensor. is the bulk modulus and is the spin tensor. and are the first and second invariant of , respectively. The Jaumann stress rate is a commonly adopted approach in SPH-related studies (Bui et al., 2008, Nguyen et al., 2017, Feng et al., 2021). However, it has been shown to exhibit oscillatory behavior under simple shear conditions (Dienes, 1979), as also mentioned by del Castillo et al. (2024). Interestingly, del Castillo et al. (2024)’s investigation further demonstrates that such oscillations were absent in SPH simulations involving both simple shear and large deformation scenarios. Alternatively, other objective stress rates, such as those derived from the Lie derivative (Borja and Tamagnini, 1998), may also be considered in this context.

in Eq. (10) is the plastic multiplier and can be determined by the consistency condition. The change rate of the plastic multiplier is given by (13) , and are material constants in Drucker–Prager model, which are defined as (Borja, 2013, Nguyen et al., 2017) (14) Here, , , and are cohesion, friction angle, and dilation angle, respectively. Then the trial stress at the new time step can be obtained by (15) A two-step elastic predictor-plastic corrector scheme, known as the return mapping algorithm (Simo and Hughes, 2006, Bui et al., 2008, Borja, 2013), is employed to return the trial stress to the yield surface and get the final stress .

A simple benchmark case, specifically the evolution of a circular water drop into an elliptical shape (Monaghan, 1994, Khayyer et al., 2017), is conducted to validate the proposed UTVF for free-surface flows. The results are compared with those obtained using the TVF, as well as other particle regularization techniques, i.e., dynamic stabilization (DS) (Tsuruta et al., 2013), particle shifting (PS) (Xu et al., 2009), and optimized particle shifting (OPS) (Khayyer et al., 2017). The initial radius of the circular drop is 1 m, and it begins to move under the influence of the initial velocity field () m/s (as shown in Fig. 3), with and being the coordinates of each particle. The fluid density is set to 1000 , and viscosity is not considered. Since there are no external forces throughout the process, theoretically, the linear and angular momentum should be fully conserved.

In this section, we examine a rotating square water patch with significant free-surface deformation to evaluate the performance of the proposed UTVF method. This benchmark was first proposed in Colagrossi (2005) and has been widely utilized as a standard numerical test for particle-based simulations (Khayyer et al., 2017, Le Touzé et al., 2013). The initial configuration of the square patch is shown in Fig. 7, with a side length of . The initial velocity profile, as shown in Fig. 7a, is given by Colagrossi (2005) (43)where denotes the angular velocity corresponding to the pure rigid rotation of the fluid patch. The initial pressure distribution is obtained by solving the Poisson equation under the incompressible assumption (Colagrossi, 2005), and is defined as (44) where and . In this case, the fluid density , m and rad/s.

Next, we use a simple elastic case with a theoretical solution, to further validate the convergence and accuracy of the proposed UTVF method. The previous two fluid cases demonstrated that, compared to TVF, UTVF results in a more uniform distribution of free-surface particles, although it has little impact on the macroscopic deformation. We will further show that such free-surface particle clustering in solid materials can lead to severe consequences. As illustrated in Fig. 10, we first validate the proposed method using an elastic beam with one edge fixed. The results are then compared to existing theoretical (Landau and Lifshitz, 2013) and numerical (Gray et al., 2001, Zhang et al., 2024a) solutions. The beam has a length and thickness , with the left edge fixed to create a cantilever configuration. An observation point is placed at the midpoint of the rear end to record the vertical displacement (deflection). The material and dimensional parameters are taken from the literature (Gray et al., 2001, Zhang et al., 2024a), i.e., density , Young’s modulus , Poisson’s ratio , m, and m. An initial velocity , applied perpendicular to the beam, is given by (45)where is a constant that controls the magnitude of initial velocity, is the sound speed, and is a function defined as (46)where is determined by the equation . The frequency of the oscillating beam is theoretically given by

(47)

The temporal evolution of the deflection at the observation point is shown in Fig. 12 for various resolutions using the present UTVF. It can be observed that as the resolution increases, the difference between adjacent resolutions decreases, indicating the convergence of the present method.

Table 1 shows the first oscillation period obtained from the present UTVF, TVF, and the analytical solution. For comparison, results from the artificial stress method (SPH-AS) (Gray et al., 2001), commonly used in solid dynamics, have also been included. Clearly, the use of UTVF reduces the error relative to the theoretical value compared to TVF. Additionally, the error of SPH-UTVF is at the same level as that of SPH-AS. While the use of artificial stress can yield relatively accurate periodic values, it introduces significant dissipation. Fig. 13 shows the long-term temporal evolution of elastic strain energy, kinetic energy, and total energy in SPH simulations using both artificial stress and UTVF. When with artificial stress, total energy rapidly decays, whereas the present UTVF method demonstrates improved energy conservation. Additionally, after a certain period in the simulation, the elastic strain energy with artificial stress does not return to zero even when the beam returns to its initial position (where kinetic energy reaches its maximum). This is due to the hourglass issues (Zhang et al., 2024a) in the later stages of the simulation when using artificial stress, which causes a zigzag distribution of particles and stress. As a result, even when the beam returns to its initial position, the stress of some particles remains non-zero, as illustrated in the literature (Zhang et al., 2024a). The study by Zhang et al. (2024a) also demonstrated that in elastic solids, such particle clustering and non-physical fractures are caused by hourglass modes. Nevertheless, this case is presented here to illustrate that the proposed UTVF effectively eliminates these numerical instabilities.

The collapse of a 2D cohesionless granular column is simulated. The model setup is shown in Fig. 14, with material parameters detailed in Table 2 (Nguyen et al., 2017). The numerical simulations are designed to replicate the experimental conditions (Nguyen et al., 2017), where the granular column is released under the influence of self-gravity (gravity acceleration ). 2D granular flows are examined for model lengths () of 0.2 m and 0.1 m separately, while the model height () is consistently maintained at 0.1 m. The initial particle spacing is set to m.

A rectangular soil mass with a length of m and a height of m fails under its own weight, following the setup in Bui et al. (2008). Due to the cohesive nature of the soil, severe tensile instability will occur in the tensile region if only the original SPH formulation is used. In this section, we apply the proposed UTVF to this problem and compare the results with those obtained using TVF and the unmodified approach. Literature results based on artificial stress (Bui et al., 2008) and particle shifting (Lallemand et al., 2025) are also included to demonstrate the validity of the findings. The material parameters, consistent with those in Bui et al. (2008), are listed in Table 2.

Fig. 18 shows the simulation results at various time points when no correction, TVF, and UTVF are applied, respectively. In each subplot, an enlarged inset in the lower-left corner clearly shows the local particle distribution. It can be observed that the original SPH, when simulating cohesive granular materials, exhibits severe particle clustering and non-physical fractures in the tensile region, indicating tensile instability. This instability significantly affects the computational results, leading to an over-prediction of the sliding distance in this case. When TVF is applied, the tensile instability of inner particles is mitigated; however, particle clustering still occurs for free-surface particles. With the present UTVF, tensile instability is completely eliminated for both inner and free-surface particles, resulting in a uniform particle distribution. Upon reaching a stable state ( s), a shear band extending from top to bottom can be observed in the right region, while the area to the left of the shear band remains undisturbed. This result is consistent with findings reported in the literature (Bui et al., 2008, Lallemand et al., 2025).

As mentioned earlier, the extension of the artificial stress method to three dimensions is complex (Bui and Nguyen, 2021, Zhang et al., 2024a). In this section, we will demonstrate that the present UTVF can effectively eliminate tensile instability in 3D scenarios, without introducing additional complexity compared to the 2D case.

The flow of cohesive granular materials on an inclined flume with a barrier is simulated to demonstrate the capability of the proposed UTVF in preventing tensile instability under significant deformation conditions. The model setup is shown in Fig. 25a, where a rectangular soil mass is placed on the top of the flume. The inclined flume forms an angle of 50 degrees with the horizontal plane. The material parameters are as listed in Table 2.

Fig. 25 presents snapshots of the numerical simulation at various time steps, colored by velocity magnitude. The granular material moves downward under the influence of gravity and eventually comes to rest in front of the barrier. The magnified inset images reveal that tensile instability does not occur throughout the entire process, with the particle distribution remaining uniformly distributed. This indicates that the proposed UTVF is effective in preventing tensile instability even in 3D granular flows.