1 Introduction

Successful modelling of the interaction between a mass of granular sediment and in a highly transient open flow, a situation commonly encountered in many natural formations and industrial processes, requires consideration of many complex physical phenomena that are challenging to incorporate into any contemporary numerical scheme. However, it is this complexity, even when dealing with specially designed physical experiments, which makes numerical modelling an extremely useful tool to gain insight into the evolution of such flow patterns, the magnitude of local speeds or pressures, faultline growth, stagnation points, and several other clearly defined variables that are very difficult to measure and visualize in a physical set-up. The discharge of a channel or river bedload dynamics at transitional points are some general cases of interest, along with sand traps or spillways that can become silted up during normal operation. Another important instance for accurate prediction of this family of flow dynamics is in landslide risk assessment, where literature on multi-phase modelling has focused on wave generation, as in the studies by [1,2,3]. In this type of case, the added benefit of a complex model for the granular phase can enhance both numerical performance and better similarity to physical observations.

In applications like these, since a sufficiently significant mass of saturated sediment, silt, or slurries can behave like a fluid under certain conditions but alternatively keep a relatively solid shape if the local forces are not enough to have it yield, it is not uncommon to apply nonlinear rheological models to a clearly identified granular phase along with any other physical characteristics, most importantly density and internal friction qualities. In the types of problems such as those addressed here, this nonlinear phase can become temporarily suspended, but not fully mixed or dissolved, within a linear water phase [4,5,6].

Aside from the fact that the movement of granular masses depends on tribological phenomena among clusters of grains, which are already difficult to describe deterministically, these and any other physical nonlinear characteristics that may arise can affect the stability of numerical schemes that otherwise exhibit good evolution of data fields and even excellent accuracy when validated with physical experiments if properly set up [7, 8]. Care must be taken to avoid numerical oscillations, which can become very difficult to overcome when supersonic speeds are expected, like the modelling of explosions in soil [9].

Previous approaches to the sediment scouring and transport problem as in applying the Boussinesq equation [3] or the shallow water equations (SWEs) [10] are useful only for long channels and riverbeds, but not for two or three-dimensional flows with rapidly changing geometries at a local level. For low aspect ratio geometries, studies using finite volume methods (FVM) have followed approaches that include changing viscosity using different rheological models [6, 8, 11] with varying degrees of success, but accuracy is hindered by sensitivity to the remeshing algorithms that must adapt continuously to moving boundaries formed between phases and at the free surface. A form of remeshing in SPH, that of particle shifting [12] or particle refinement techniques [13], cannot yet guarantee stability unless it is supported with other stabilization techniques. Transportation models of solid particles in multiphase flow were outlined early in SPH development [14], but the incorporation of more complex viscosity models and yield criteria did not begin until over a decade later, as accounted by [15], and is still a developing field in particle methods, as outlined in [16] regarding multiphase granular flow.

Due to its meshless nature, the modelling of a nonlinear generic fluid with Smoothed Particle Hydrodynamics (SPH) was attempted by [17], and later on, a good amount of work was done on nonlinear flows using variants of SPH [18, 19], greatly widening the applicability of SPH schemes. Hybrid approaches mentioning the presence of nonlinear flows that combine SPH and FVM have also been tested but with the FVM component applied to interaction with elastic structures [3, 8, 20].

Previous studies that have focused on the nonlinear nature of two-dimensional sediment flow [21,22,23] approached an application with SPH of the well-known Herschel–Bulkley (HB) or Herschel–Bulkley–Papanastasiou (HBP) rheological models; farther along, [7] implemented, among other aspects, a varying yield stress model for the granular phase according to its local effective pressure. Modification of the kernel function and adding correction terms to the conservation equation, such as the approach by [24] for an Oldroyd-B fluid, is another way to overcome instabilities and improve precision. These cases are also two-dimensional, and the agreement with experimental data is good. Other studies, such as that of [25], explore alternatives for a yield criterion at the nonlinear phase or prefer to solve a Laplace equation [26] to calculate effective stress in the unyielded nonlinear phase, with the purpose of combining an elastic model with a multiphase, liquid-only SPH formulation. Along the line of using a pseudoelastic model for the yielding process in a viscoplastic fluid, studies like [27] approach the problem by using an artificial stress term to evolve particle kinematics. Very recently, work like that by [28] has concentrated on smoother stress fields to dampen noise detrimental to other data fields. Treatment of the boundaries in this type of multiphase flow is another matter, which has received plenty of attention, like the application of normalization factors at free surfaces as described by [29], application of surface tension equations with the SPH approximation as described in [30], or the usage of artificial repulsion forces at the interface as in the proposal by [31] to deal with clumping issues.

Later applications of this combination of techniques to specific two-dimensional cases of two-phase dam break offer different perspectives on flow dynamics exhibited by the numerical models such as an underwater dam break using SPH [32], a silted-up dam break using FVM for both phases [11] with some comparison to experimental data or even the presentation of three-dimensional cases with limited spatial resolution, like the single-phase debris flow SPH model by [33], among other similar studies.

To address several of the challenges these studies have encountered when dealing with multi-phase flow, this study aims to lay out a numerical scheme that can accurately emulate complex flow dynamics while maintaining stability. To properly replicate pressure evolution, the scheme should prevent unstable fields that can be caused both by nonlinear variations of viscosity and by the sharp changes in density due to dissimilar phases.

The remaining sections are organized as follows: first, a description of the mathematical model is provided along with the methodological aspects of the implementation to cope with a nonlinear constitutive model. The focus is on the nature of the nonzero yield stress of the granular material and the presence of two distinct phases. Next, the model is tested and compared with two experimental series: a single-phase nonlinear dam break, and a multiphase dam break based on a silted-up spillway configuration. In each two-dimensional case, validation is done with published physical experiment data sets. Finally, the concluding section presents the main takeaways from this study and outlines future research in this line of work in accordance with the results obtained.

2 Numerical model

The present work will employ a weakly compressible SPH (WCSPH) formulation for a general fluid. All phases considered will be modelled under the same physical laws, either in a single phase or a multiphase context. For a weakly compressible flow, the Lagrangian form of Navier–Stokes equations can be written as follows:

$$\begin{aligned} \left\{ \begin{array}{l} \displaystyle { \frac{D \rho }{D t}=-\rho \ \nabla \cdot \textbf{u} } \\[0.5cm] \displaystyle { \frac{D \textbf{u} }{D t}= \frac{1}{\rho } \nabla \cdot \sigma + \textbf{g} }\\[0.5cm] \displaystyle { \textbf{u}=\frac{d \textbf{r}}{d t}, \qquad p = p(\rho ) \, , } \end{array} \right. \end{aligned}$$
(1)

where \(\rho \) is the density, \(\textbf{u}\) the velocity, \(\textbf{r}\) the position, p the pressure, \(\textbf{g}\) the external volumetric force due to gravity, and \(\displaystyle {\frac{D}{D t}}\) stands for the material derivative. The term \(\sigma \) represents the Cauchy stress tensor, which can be split into two terms, namely the normal stresses p and the shear stresses\(\tau \) as in Eq. (2):

$$\begin{aligned} \displaystyle { \sigma = - p \, \mathbb {I} + \mathbf {\tau }} \, . \end{aligned}$$
(2)

In SPH, the system of Eqs. (1) becomes system (3), where \(\langle {} \rangle \) is the customary notation for the SPH approximation:

$$\begin{aligned} \left\{ \begin{array}{l} \displaystyle { \left\langle {\frac{D \rho }{D t}} \right\rangle _i = -\rho _i \ \sum _j V_j \left( \textbf{u}_j - \textbf{u}_i \right) \nabla _i W_{ij}}, \\[0.5cm] \displaystyle { \left\langle {\frac{D \textbf{u} }{D t}} \right\rangle _i = \frac{1}{\rho _i} \ \sum _j V_j \left( \mathbf {\sigma }_i + \mathbf {\sigma }_j \right) \nabla _i W_{ij} + \textbf{g}, }\\ \displaystyle { \textbf{u}_i =\frac{d \textbf{r}_i}{d t}, \qquad p = c_s^2 \left( \rho _i - \rho _0 \right) ; } \end{array} \right. \end{aligned}$$
(3)

where subscript i represents the particle being computed, and j is the set of neighbouring particles. Term \(V_i = m_i / \rho _i\) is the volume of particle i, and \(W_{ij}:= W(r_i - r_j)\) and \(\nabla _i W_{ij}\) the kernel and kernel derivative, respectively. The kernel, which acts as a smoothing function, is a Wendland C2 kernel [34] with \(h/dr = 2\); h is the smoothing length and dr inter-particle distance unless otherwise stated.

With this set of equations fully defined, hypotheses are required for the forms of \(\mathbf {\tau }\) and p. The former will be achieved through a rheological model, discussed further along, and the latter through an equation of state (EOS). The EOS is the stiff equation shown in the system of Eq. (3). This way compressibility is adjusted so that the speed of sound is artificially lowered according to the scheme by [35]. Usually, the linear term is sufficient to retain the properties of the EOS.

In the present analysis, an apparent viscosity approach is used to account both for potential variation of viscosity as a function of strain rates and the existence of plastic behaviour, where yield stress \(\tau _0\) must be overcome before flow can begin. Therefore, viscous stresses (\(\tau \)) are computed as:

$$\begin{aligned} \tau _i = 2\,\mu \, \dot{\varepsilon }_i \, , \end{aligned}$$
(4)

where \(\mu \) represents the viscosity and \(\dot{\varepsilon }_i\) is the deviatoric strain rate tensor, defined as:

$$\begin{aligned} \dot{\varepsilon }_i = \frac{1}{2} \left( \nabla \textbf{u} + \nabla \textbf{u}^T \right) - \frac{1}{3} \nabla \cdot \textbf{u} \, \delta _{ij} \, , \end{aligned}$$
(5)

with \(\delta _{ij}\) being the Kronecker delta. In SPH, the \(\nabla \textbf{u}\) term can be computed as:

$$\begin{aligned} \left\langle { \nabla \textbf{u} } \right\rangle = \sum _j V_j \left( \textbf{u}_j - \textbf{u}_i \right) \otimes \nabla _i W_{ij} . \end{aligned}$$
(6)

When dealing with wet silt or sand, if the sediment is fully suspended, viscosity is a function of solid concentration [4] but does not exhibit any yield stresses since contact between the solid particles is limited and does not allow for frictional interlocking. This has been well documented and applied as an interface model for SPH by [4] and followed by later publications similar to that of [36].

Another important factor when dealing with granular materials is local packing pressure, as it significantly changes local yield stress. This is a well-known geotechnics phenomenon and has been applied in many studies using SPH, among them [4, 23, 26, 28]. These studies use the Drucker–Prager yield model that considers two factors in determining local yielding stress \(\tau _0\), as in Eq. (7):

$$\begin{aligned} \tau _0 = -\alpha p_{\textrm{eff}} + \kappa _s \, . \end{aligned}$$
(7)

In Eq. (7), \(p_{\textrm{eff}}\) is the local effective pressure, \(\alpha \) relates variation according to this pressure, and \(\kappa _s\) is the contribution of cohesiveness. Both of these constants can be estimated by projecting this model onto the Mohr–Coulomb yield criterion. For a yield surface at \(\pi /6\), internal angle of friction \(\phi \) and cohesiveness c, \(\alpha \) and \(\kappa _s\) can be calculated as follows:

$$\begin{aligned} \alpha&= \frac{2 \sqrt{3} \sin (\phi )}{3 - \sin (\phi )}, \end{aligned}$$
(8)
$$\begin{aligned} \kappa _s&= \frac{2 \sqrt{3} c \sin (\phi )}{3 - \cos (\phi )}. \end{aligned}$$
(9)

Several approaches may be used to compute effective pressure in SPH. Free surface and interface tracking for all the phases are considered here. The effective pressure is then computed from total pressure (\(p_\textrm{tot}\)) and pore water pressure for saturated sediment (\(p_\textrm{pw}\)) as, for instance, [23]:

$$\begin{aligned} p_{\textrm{eff}} = p_\textrm{tot} - p_\textrm{pw} = \rho g \left( h_\textrm{tot} - h_\textrm{pw} \right) \end{aligned}$$
(10)

where \(h_{tot}\) and \(h_{pw}\) are the total depth and the local depth, respectively, of the fluid.

In nonlinear fluids such as those with viscoplastic behaviour, a variation of the viscosity can be described using any equation that properly fits the phenomenon, among which the Herschel–Bulkley–Papanastasiou (HBP) rheological model [37] (a generic form shown in Fig. 1) is well known for its versatility since it can consider yielding stress and either a shear thickening or shear thinning fluid through flow index n.

Fig. 1
figure 1

Conventional apparent viscosity approach

Full size image

According to this rheological model, once movement begins, viscous shear stress \(\tau _i\) will vary as a function of angular shear rate \(\dot{\varepsilon }_i\); therefore, an apparent viscosity \(\mu _{i,\textrm{eff}}\) can be defined as the proportion \(\tau _i/\dot{\varepsilon }_i\), which will become large for shear rates close to zero. Therefore, the model approximates this viscosity as:

$$\begin{aligned} \displaystyle \mu _{i,\textrm{eff}} = \kappa \dot{\varepsilon }_i^{n-1} + \frac{\tau _0}{\dot{\varepsilon }_i} \left( 1 - e^{-m \dot{\varepsilon }_i} \right) , \end{aligned}$$
(11)

where parameters \(\tau _0\), flow index n, consistency index \(\kappa \), and stiffening rate m must all be determined by fitting the equation to appropriate experimental data.

To circumvent numerical problems with the second term when \(\dot{\varepsilon }=0\), an equivalent function that does not become undefined at that point is used, as in [38]:

$$\begin{aligned} \displaystyle \mu _{i,\textrm{eff}} = \kappa \dot{\varepsilon }_i^{n-1} + \tau _0 T(m \dot{\varepsilon }_i ) , \end{aligned}$$
(12)

where T(x) is a function defined as \((1-e^{-x})/x\) for values of x above a small threshold, but as its MacLaurin series (truncated after 8 terms), for smaller values. The threshold value, which can be visualised as an amount \(m \dot{\varepsilon }_0\) (see Fig. 1), should be chosen in accordance with the nonlinear properties of the fluid, but it does not produce a sharp discontinuity as the threshold is crossed, and does not require yield surface detection algorithms.

Whatever the method applied to obtain an effective viscosity, even small oscillations in the pressure field can produce apparent yielding of the fluid when, in reality, movement at a macroscopic level has not yet begun. This means that effective viscosity will decrease and then yielding happens before it should. This is especially important where pile-up, buckling, or scouring needs to be modelled. To avoid evolving an additional elastic stress equation and simultaneously track the yielded surface to determine whether such an equation applies or not, one approach used by [23] and followed up on by [7] is to evaluate the yield criterion and set a null acceleration to any unyielded particle. Yielded particles will then follow the applicable momentum equation, an approximation used here for cases where pileup or scouring is expected during the simulation. Therefore,

$$\begin{aligned} \left\{ \begin{array}{lll} \dot{\varepsilon } = 0 &{} \text{ if } &{} \tau \le \tau _{\textrm{y}} \\ \tau = 2 \mu _{\textrm{eff}} \dot{\varepsilon } &{} \text{ if } &{} \tau >\tau _{\textrm{y}} \end{array}\right. \end{aligned}$$
(13)

as proposed by [39], where \(\tau \) is evaluated as:

$$\begin{aligned} \tau =\sqrt{\frac{1}{2} \sum _{i, j} \tau _{i j} \tau _{i j}}. \end{aligned}$$
(14)

When applying the yield criteria from Eq. (13) to the nonlinear phase numerically, forcing zero acceleration at the nonyielded area tends to destabilize pressure fields. For that reason, this model proposes to keep the evolution of the conservation equation instead of assigning zero acceleration. In place of this, the conservation equation is multiplied by an attenuation factor, \(K_R\), which allows for a much smoother transition at the yielded/non-yielded interface, where sharp transitions are normally found in erosion processes.

To remove high-frequency pressure oscillations, the \(\delta \)-SPH term is added to the continuity equation; this also helps improve the overall stability of the scheme. The \(\delta \)-SPH term is applied independently to each of the phases involved, based on its impact on stability according to the multiphase study by [40]. The term is added to the continuity equation, yielding to:

$$\begin{aligned} \displaystyle {\frac{\textrm{d} \rho _i}{\textrm{d} t}} \displaystyle {= - \rho _i \sum _{j} \textbf{u}_{ji}\cdot \nabla _{i} W_{ij} V_j \, + h \, c_s \, \delta \sum _{j} \mathbf {\psi }_{j i} \cdot \nabla _{i} W_{ij} V_{j}}, \end{aligned}$$
(15)

where the term \(\varvec{\psi }_{ji}\) is defined as:

$$\begin{aligned} \displaystyle {\mathbf {\psi }_{ji} = 2 \left[ (\rho _j - \rho _i) - \frac{1}{2} \left( \left\langle {\nabla \rho } \right\rangle _{j}^{L} + \left\langle {\nabla \rho } \right\rangle _{i}^{L} \right) \cdot \textbf{r}_{ji} \right] \frac{\textbf{r}_{ji} }{\Vert \textbf{r}_{ji} \Vert ^2}}, \end{aligned}$$
(16)

where \(\left\langle {\nabla \rho } \right\rangle _{i}^{L}\) represents the renormalized density gradient, as defined by [41]. The value for \(\delta \) is always set to 0.1.

To improve the spatial distribution of particles, a particle shifting algorithm is used, following the solution proposed by [42], based on a diffusion law. Particle positions are slightly shifted from high-concentration areas to low-concentration areas to conserve a regular distribution. According to [42], the expression for the shifting term for all particles in domain D becomes:

$$\begin{aligned} \delta {\textbf {u}}_i&= -0.5 \max _{j\in D} \left( \left\| ({\textbf {u}}_j-{\textbf {u}}_i) \cdot \frac{{\textbf {x}}_j- {\textbf {x}}_i}{\left\| {\textbf {x}}_j- {\textbf {x}}_i \right\| } \right\| \right) \nonumber \\ {}&\quad \left( \frac{R}{\Delta x}\right) ^3 R \hat{\nabla } C_i \end{aligned}$$
(17)
$$\begin{aligned}&if \quad \left\| \left( \frac{R}{\Delta x}\right) ^3 R \hat{\nabla } C_i \right\| < \frac{R}{2 \Delta x}; \end{aligned}$$

otherwise,

$$\begin{aligned} \delta {\textbf {u}}_i = -0.5 \max _{j\in D} \left( \left\| ({\textbf {u}}_j-{\textbf {u}}_i) \cdot \frac{{\textbf {x}}_j- {\textbf {x}}_i}{\left\| {\textbf {x}}_j- {\textbf {x}}_i \right\| } \right\| \right) \left( \frac{R}{2\Delta x}\right) \hat{\nabla } C_i, \end{aligned}$$
(18)

where \(\hat{\nabla } C_i\) is the concentration gradient approximation by [43].

For free-surface flows, shifting at the free surface must be controlled to prevent overall fluid volume to increase along the simulation due to particles shifting towards low-concentration areas beyond the free surface. In this work, we follow the approach initially proposed by [44], and further developed by [42, 45].

To this means, identification of particles near the surface is done using the work developed by [45, 46]. These use thresholds to designate different regions depending on the distance of the particle from the free surface, based on the minimum eigenvalue (\(\lambda _i\)) of the renormalization matrix:

$$\begin{aligned} \varvec{L}_i\,:=\,\left[ \sum _j \textbf{r}_{ji}\otimes \nabla _{i} W_{i j} V_{j}\right] ^{-1} \, . \end{aligned}$$
(19)

In the first step, particles with eigenvalues \(\lambda _i \le 0.2\) are identified as free surface particles. For particles within an intermediate region (\(0.2< \lambda _i < 0.7\)), a second step is performed. Then, free surface particles can be fully identified, and the normal vector is found according to the following expression:

$$\begin{aligned} \varvec{n}_{i}&=\frac{\langle \nabla \lambda \rangle _{i}}{\left\| \langle \nabla \lambda \rangle _{i}\right\| } \quad \text{ with } \quad \langle \nabla \lambda \rangle _{i}\nonumber \\ {}&= -\sum _{j}\left( \lambda _{j}-\lambda _{i}\right) \varvec{L}_{i} \nabla _{i} W_{ij} V_{j} \end{aligned}$$
(20)

Boundaries will be treated here by using numerical boundary integrals, according to the approach initially proposed by [47] and further extended by [29]. In that technique, the boundary is split into elements distinguished by its geometric attributes: area, normal vector, and tangent vector. This implies an assumed truncation error, but it prevents node connectivity and is more efficient in computational cost [26, 48]. The Shepard normalization factor is computed geometrically with a semi-analytical approach.

Time integration is carried out by means of an explicit predictor-corrector scheme. Timestep size is chosen at each integration step according to the following restrictions:

$$\begin{aligned} \begin{array}{ll} \displaystyle { \Delta t \le \text {CFL} \frac{h^{2}\rho _i}{\mu _{i,\textrm{eff}}} } \\ \\ \displaystyle { \Delta t \le \text {CFL} \min _{i}\left( \frac{h}{c_{0}+h \max _{j} \pi _{i j}}\right) } \\ \\ \displaystyle { \Delta t \le \text {CFL} \min _{i} \left( \frac{h}{\left\| \varvec{u}_{i}\right\| }\right) } \\ \\ \displaystyle { \Delta t \le \text {CFL} \min _{i} \sqrt{\frac{h}{\left\| \varvec{a}_{i}\right\| }} } \end{array} \end{aligned}$$
(21)

The current model and the simulations presented below have been carried out with AQUAgpusph, an open-source SPH software developed by CEHINAV-UPM [47].

Table 1 Physical constants for numerical models
Full size table

3 Results and discussion

The test cases discussed here aim to cover several types of physical interaction to illustrate the effect of viscoplastic nonlinearity of dam break dynamics on flow profiles and the numerical performance of the proposed method. First, a single-phase dam break case will provide good insight as to the impact of viscoplasticity on the rest of the SPH scheme and is expected to replicate key physical behaviour elements that are expected, like a certain degree of pile-up and strong viscous damping of angular strain in the nonlinear phase.

Then a two-phase case, consisting of a silted-up dam break, is targeted at reproducing simultaneous instances of dynamic interactions between a linear and a nonlinear phase. Dam break evolution for both phases will be different, of course, and consequently, differences in speeds will bring up complex local erosion and scouring phenomena. Recall that the model for the granular phase implies the additional variability of yield stress as a function of local pressure, a critical variable in the stable evolution of the conservation equation in weakly compressible formulation. All these complications also mean that pressure and strain rate fields will tend to be unstable. These conditions tax the numerical scheme and therefore will be a good test of its adaptability to rapidly evolving yielded/nonyielded boundaries with multidirectional interaction between the linear and nonlinear phases.

Calculation results are compared directly to recently published experimental data for two distinct cases: (1) single-phase viscoplastic clay, a two-dimensional dam break with relatively low yield stress, and (2) the silted-up dam break, a two-dimensional sand and water interaction experiment which uses uniform, saturated quartz sand as the nonlinear phase. Each case used known physical properties applicable to each phase; since the nonlinearity of such properties can have a significant impact on accuracy and stability, the parameters for viscoplasticity are shown in Table 1 as a general reference for the rest of this article.

The implementation of the numerical model described above can handle n different nonlinear phases and is limited to nongaseous fluids, but here only two-phase cases are tested. This n-phase scheme has been implemented as a plugin within the fluid solver of AQUAgpuSPH.

3.1 Single-phase non-Newtonian dam break

The viscoplastic dam break studied here replicates the physical experiment by [49], which used homogeneous, water-saturated clay with rheological constants fitted to the Herschel–Bulkley (HB) model. This case was worked on recently by [50] who proposed the application of the volume-of-fluid (VOF) method to replicate the flow dynamics. The geometric set-up for the experiment is shown in Fig. 2 and will be considered two-dimensional since the tank width is constant and the water-saturated condition of the clay guarantees good lubrication in relation to all solid boundaries.

Fig. 2
figure 2

Case geometry for 2D nonlinear dam break

Full size image
Fig. 3
figure 3

Convergence of nonlinear dam break at \(t = 0.5\) s (\(t_d\) = 4.3)

Full size image
Fig. 4
figure 4

Convergence of nonlinear dam break at \(t = 1.0\) s (\(t_d\) = 8.7)

Full size image

The clay is initially confined by a bottom solid wall and two vertical solid boundaries; when the gate is instantaneously removed, the clay flows out along the horizontal bed, wet with a thin film of water. The height of the clay column is \(H_f =\) 0.13 m, while the length is \(L_f=\) 0.26 m.

The density of the fluid has been reported by [50] as \(\rho =\) 1540 kg/m\(^3\) along with HBP rheological parameters \(\tau _0 =\) 30.0 Pa, \(n =\) 0.472 and \(\kappa =\) 4.472. Gravity acts along the vertical axis and was considered as \(g =\) 9.81 m/s\(^2\).

Numerical parameters for the weakly compressible scheme are chosen according to low expected flow speeds, so an artificial speed of sound is chosen as \(c_s =\) 40 m/s, which is well over the criterion provided by [51] \(c_s = 10\) V\(_{\max }\). Reference speed \(V_{\max }\) can be estimated as \(\sqrt{(}2gH)=1.60\) m/s, but will likely be lower since viscosity is high.

Spatial resolution varies according to the number of particles used, but for the model with \(N =\) 80,000, for which numerical convergence was reached (as noted after several test runs with N between 5000 and 124,000 fluid particles), the initial particle spacing for the converged model was \(\Delta x = 7.5 \times 10^{-4}\) m. All solid boundaries were modelled with boundary integrals as described above, and no viscous or cohesive forces were considered at the walls. Convergence is sampled at t = 0.5 s (\(t_d\) = 4.3) with profiles seen in Fig. 3 and at t = 1.0 s (\(t_d\) = 8.7) in Fig. 4. Considering how simple the geometry is, convergence occurred at a relatively high resolution; this was probably related to the fact that enough particles must be present at the advancement front to replicate pileup properly, and even the converged models show higher error in this area compared to the mass farther back. This can be due to the fact that the mass that is in the front of the flow and close to the solid boundary has been through much higher shear strain than other areas; close inspection of the discrete distribution reveals that the SPH particles are packed into a higher, albeit smoothed, local density. If the mass is to remain constant, mean volume must decrease; this is an inaccuracy observed across test cases with high shear strain and the relatively low speeds implied with a viscoplastic collapse. When looking at these results, recall that yield stress for this fine, water-saturated clay is quite low, so static pileup angles would be much smaller than those expected from sand or other coarser granular materials.

Fig. 5
figure 5

Collapsed profiles for 2D nonlinear dam break

Full size image
Fig. 6
figure 6

Viscosity and pressure fields at \(t = 0.250\) s (\(t_d = 2.2\)) for 2D nonlinear dam break

Full size image

The results for this calculation are shown in Fig. 5 at times \(t =\) 0.500 s, 0.625 s, and 1.000 s (approximate dimensionless times \(t_d =\) 4.3, 5.4 and 8.7) in correspondence with published experimental data.

Several indicators of a viscoplastic collapse are present, such as the conservation of the initial rectangular shape for a longer time than a linear, less viscous fluid, with no splashing, rollover, or wave formation. To showcase contrasting yielded and unyielded areas in the fluid, the effective viscosity field is plotted out in Fig. 6 for t = 0.250 s (\(t_d = 2.2\)); in areas where particles are at low or zero shear rate, local viscosity is higher in accordance with Eq. (12), but at the advancement front, viscosity is correspondingly lower. In some areas, numerical oscillations resulted in locally higher angular shear rates, which reduces effective viscosity in a manner that the fluid will yield prematurely, but here it did not seem to change accuracy significantly.

This viscosity field is smooth, but not directly linked to the pressure field since a constant value for yield stress was used. This does not mean that a noisy pressure field could not affect scheme stability, though, because the calculation of local apparent viscosity is a function of angular shear rate \(\dot{\varepsilon }\) associated directly with the conservation equation (1). In this set of numerical experiments, once the resolution was high enough, pressure fields were very smooth during the \(t = \) 2.0 s of physical time (\(t_d = 17\)) tested. As an example, Fig. 7 shows a time instant of a simulation at \(t=\) 1.0 s (\(t_d = 8.7\)) where this field can be appreciated.

Fig. 7
figure 7

Viscosity and pressure fields at \(t = 1.0\) s (\(t_d = 8.7\)) for 2D nonlinear dam break

Full size image

Although this is a viscoplastic mass that is collapsing slowly in comparison with a linear, low-viscosity fluid, pileup angles will be nonzero but low since the order of magnitude of yield stress \(\tau _0\) is below \(10^2\) Pa. Hence, a clearly recognizable segment mass has not completely collapsed (for example, the area in Fig. 6 where apparent viscosity is over 50 Pa s) but a flow front that is completely dynamic has flattened and is moving quickly towards the right boundary. At \(t=\) 1.0 s, although barely recognizable, a faultline has formed between the main collapsing mass and the completely dynamic part of the fluid.

Aside from excellent accuracy in relation to physical experiments, no unexplained waves, splashes, or flow rebounding appeared within the simulation timeframe. However, since no artificial viscosity stabilization scheme is used for this nonlinear phase, pressure field instabilities will form for higher values of \(c_s\). The application of the \(\delta \)-SPH diffusion term and the particle shifting techniques greatly reduced the occurrence of fringes and nonphysical faultlines; likewise, free surface and interface treatment were key in preventing numerical cracking of the surface.

While testing sensitivity and tuning the numerical parameters of the model, several potential sources of pressure instability were recognized, some of which may be harmful feedback phenomena between pressure and strain. One possible source of instability can be values of \(\tau _0\) in the order of \(10^3\) Pa, and even at \(10^2\) Pa if flow speeds become high enough. In such cases, nonphysical, multiple faultlines will form or even render the pressure field catastrophically unstable. In fluids where apparent viscosity is described by the HBP rheology model, values for consistency index n over the vicinity of 1.2 also accelerate the occurrence of instability in the pressure field, perhaps due to stronger viscous stress differentials that are induced.

During the developmental stages of the scheme, lack of accuracy usually was in the form of extensive premature yielding, causing the whole mass to collapse much like a fluid with no yield stress. Not all of the fluid yielded; it did so, at initial stages, in a manner of alternating strips of low and high viscosity roughly distributed perpendicularly to the lower boundary. This yielding of course meant movement that was quickly diffused into nonyielded areas, which then correspondingly entered a fully dynamic state with no unyielded areas. Sometimes this occurred after the end of the simulation of the timeframe of interest, but barely noticeable evidence of potential instabilities was appearing even then, and eventually, the fluid would yield completely instead of showing static pileup. In this single phase case, the particle shifting algorithm kept several datafields stable, even if a unitary value was used for the relaxation factor \(K_R\). However, during test runs it was observed that if the physical values of \(\tau _0\) is over \(10^2\) Pa, a smaller value of \(K_R\) needs to be used to keep the unyielded phase numerically active and help stabilize the angular shear rate \(\dot{\varepsilon }\) field since it is used to determine local yield condition.

Finally, the potential for inaccuracies due to shearing boundary forces was verified by applying no-slip boundaries to a test run identical in all other aspects to the final run used to report these results. Neither accuracy nor stability seemed to be affected significantly, and in other test runs where unnatural sliding of the fluid mass occurred, it was not related to boundary interactions, but to unphysical cracking of the fluid at the free surface. Once again, the use of particle shifting and appropriate treatment at the free surface both contributed to reducing this problem significantly.

This behaviour of a single, nonlinear collapsing block of fluid is just a first approximation to the full intention of the proposed model. With the presence of a second, quicker reacting phase, and higher yield stress differentials, it is reasonable to expect much stiffer difficulties.

Fig. 8
figure 8

Snapshots of the two-phase dam-break evolution at different time instants. Sand phase takes 80% of the column and water the remaining 20%. Experiments reproduce that by [11]

Full size image

3.2 Multi-phase non-Newtonian dam break

This case deals with a multi-phase formulation of a dam break in which one of the phases is modelled with a non-Newtonian rheological model. Several complexities can be expected if compared to the previous single-phase, nonlinear dam break. First, since the nonlinear phase is sand, yield stresses will be higher in accordance with known internal friction characteristics; furthermore, the presence of a second phase means that all types of interactions can follow, including hydrostatic action (normal to the interface boundary) and erosive action (shearing). Of course, no dissolution between water and sand or physical state transformation effects such as evaporation or condensation are applicable.

In this experiment, following the geometric two-dimensional set-up shown in Fig. 9, the dam section is a layer of saturated, homogeneous sand with a comparable mass of water completing the height of this silted-up dam. A water bed will receive these collapsing masses, with the dam water eroding the silted bottom, a configuration that guarantees a good variety of interaction phenomena. As outlined briefly above, this is relevant to many engineering applications, both at large scales as in river mouths, port construction and operation, water overtopping landslide assessment, and at smaller scales as in gated channels, dam discharge outlets, large drainage pipe discharges and other instances where rapidly changing shape of an area occupied by sand or silt affects its condition or performance. On the other hand, in terms of SPH performance, this scheme must handle not only the alterations a highly non-linear rheology model introduces, but also has to deal with simultaneous dam break processes occurring at different speeds that are interacting with each other during the whole timeframe of the simulation. This is very different from bedload transportation models that can only model non-violent scouring or schemes where viscous stresses present smaller variations throughout the nonlinear phase.

To follow up on validating this scheme with these considerations in mind, physical experimental datasets were used; those published in [52] include a variety of silt depths and different wet bed depths for up to 6.00 s of physical time after the instantaneous opening of the gate. For this study, results from applying the proposed scheme to all silt depths with a downstream 0.05 m waterbed were examined, but to present relevant insights without repetition, the analysis will be carried out with results for initial silt depths \(H_s\) = 0.03 m; 0.15 m; 0.20 m; and 0.24 m. The body of water completes the total depth H, including silt, to \(H_w + H_s\) = 0.30 m depth for the dam section, which is \(L_d\) = 1.52 m long. The rest of the tank, with a total length of \(L_b\) = 6.00 m, holds a water bed with a depth of \(H_b\) = 0.05 m. In the physical experiment, the tank is D = 0.65 m wide. Drucker–Prager criterion constants associated with saturated quartz sand as described by [11] were applied here to calculate local yield stress \(\tau _0\) as prescribed by Eq. (11). Snapshots for the time evolution of the 0.24m initial silt depth case are reported in Fig. 8.

Fig. 9
figure 9

Case geometry for 2D silted-up dam break

Full size image

The numerical parameters associated with the weakly compressible scheme were chosen similarly to the single-phase case, with an artificial speed of sound of \(c_s =\) 40 m/s. The spatial resolution of the converged models that used about \(N =\) 180,000 total particles was \(\Delta x =\) 1.5 \(\times \) 10\(^{-3}\) m. All solid boundaries were modelled with boundary integrals as described above and, as before, no viscous or cohesive forces were considered at the walls, a reasonable assumption for water-saturated sand.

This set of numerical experiments showed greater sensitivity to many numerical parameters in terms of stability, especially in the pressure field. This may be related to the need to assure a sufficient degree of stress continuity throughout the different phases. Unlike the single-phase dam break, both numerical stability and accuracy in relation to experimental data were easily affected by the value of viscosity threshold factor K (m in Eq. (12)) and relaxation factor \(K_R\) that multiplies the conservation equation if an unyielded condition is detected. When properly tuned, accuracy was considered very good for dimensionless times \(t_d < 10\) for both water and sand profiles, and the scheme maintained stability for dimensionless times over \(t_d = 37\) (approximately \(t =\) 6.00 s), which is the available timeframe the physical experiment dataset provided.

In this set of tests, numerical parameter tuning was required to obtain the best degree of accuracy, starting out with the apparent viscosity factor K. This factor affects how apparent viscosity \(\mu _\textrm{eff}\) changes as a function of angular shear rates when threshold \(K \dot{\varepsilon }\) is below unity. For higher values of K, the resulting larger viscosities meant CFL criteria drove down the values for timestep, making them computationally too expensive to be practical. Since an artificially high viscosity might be regarded as an inaccuracy of sorts, it actually is taken care of differently; at low angular shear rates (that is, unyielded or very recently yielded fluid), particle movement evolves with the conservation equation but is multiplied by a relaxation factor \(K_R\). A nonzero value for \(K_R\) is needed to guarantee that the unyielded particles do not freeze or become too easily suspended in the phase that it does not belong to, where the absence of other nonlinear neighbouring particles means they cannot continue replicating viscous stresses that would have influenced it. Too large a value of \(K_R\) (if over 0.5, as an example for these test cases) meant that premature yielding occurred for the rheological models where \(\tau _0\) might be over an order of magnitude of \(10^2\) Pa. Within the scope of these tests, \(K_R = 0.01\) was a value that produced good results consistently.

Numerical convergence for this group of tests was obtained at approximately N = 180,000 or less, with little, if any, influence on the proportion of nonlinear particles present. Models with too few particles not only suffered from lower accuracy, but premature yielding occurred more readily. Care was also taken to check the convergence of both the nonlinear and linear phases separately, although some guidelines for convergence from the previous single-phase test were useful to set expectations. As an example of convergence rates for this group of tests, water and sand depth profiles for varying resolutions used for the 50% silted-up dam case can be seen in Figs. 10, 11, 12 and 13.

Fig. 10
figure 10

Convergence of water profiles for 50% silted-up dam at \(t = 1.0\) s

Full size image
Fig. 11
figure 11

Convergence of water profiles for 50% silted-up dam at \(t = 2.0\) s

Full size image
Fig. 12
figure 12

Convergence of sand profiles for 50% silted-up dam at \(t = 1.0\) s

Full size image
Fig. 13
figure 13

Convergence of sand profiles for 50% silted-up dam at \(t = 2.0\) s

Full size image

As a measure of the effect of parameters K and \(K_R\) on accuracy and consideration of any impact they may have on stability, hydrograms showing the evolution of total depth at certain dimensionless horizontal positions x/L as a function of dimensionless time are shown in Figs. 14, 15 and 16 for the case with 0.20 m of silt out of the total 0.30 m initial dam depth. Parameter x is the distance measured from the left tank wall, so position x/L = 1 is the location of the dam gate, x/L = 0.5 is within the dammed fluid volume and x/L = 1.7 is downstream relative to the gate.

Fig. 14
figure 14

Water depth at \(x/L = 0.5\) for 66% silted-up dam

Full size image

This first set of results, which remains accurate for approximately the first 2.0 s of physical time (dimensionless time \(t_d\) = 12.5) according to the experimental data, also shows that tuning parameter K has a limited effect on accuracy and can even become the cause of numerical instability. Both of these effects are to be expected, since if the HBP equation is used too soon, say for values of \(K > 0.1\), the apparent viscosity \(\mu _\textrm{eff}\) becomes high at low values of angular shear rate \(\dot{\varepsilon }\), in accordance with Eq. (12). This large value of apparent viscosity causes a stronger discontinuity close to the critical stage where movement is imminent and can interact harmfully with the yield criteria or drive the timestep to an excessively small value that is computationally impractical.

Fig. 15
figure 15

Water depth at \(x/L = 1.0\) for 66% silted-up dam

Full size image

That is also why smaller values of K produced good results even when below 0.001, where Eq. (12) led to lower apparent viscosity values and therefore kept the evolution of viscous forces gradual, so the model could rely on separate yield criteria to properly model the effect of yield stress \(\tau _0\). After the threshold is met, no significant problems can be expected to arise because the HBP equation, applicable to calculate apparent viscosity, is smoother once movement has begun because local yield stress has been surpassed. Adjustment of threshold \(K \dot{\varepsilon } = 1\) in Eq. (12) to different values results in similar effects and interactions. Changing the value of K is not equivalent to varying the threshold for \(K \dot{\varepsilon }\), but the effects are similar and do not contribute differently to accuracy or stability, at least when realistic values for density and viscosity are used.

Fig. 16
figure 16

Water depth at \(x/L = 1.7\) for 66% silted-up dam

Full size image

In every silted-up dam test case, no matter the initial silt depth, similar behaviour in accuracy was observed, both in terms of evolution in time and relative to the position. For downstream positions like that of Fig. 16 and for later stages in flow dynamics, say for dimensionless times \(t_d\) > 10, there is an obvious delay of the waterfront arrival, mainly due to the effect of a light artificial viscosity used on the linear phase only to stabilize the scheme for longer simulation times. The use of artificial viscosity is required to keep the nonlinear phase stable and helps attenuate pressure fluctuations; this, along with the application of the \(\delta \)-SPH scheme and the particle shifting technique described above, results in a smooth, almost hydrostatic pressure field, like that sampled in Fig. 17. For other cases with different initial silt levels and for most of the simulation time, pressure fields are just as smooth with small local instabilities that quickly diminish.

Fig. 17
figure 17

Pressure field at \(t = 1.0\) s (\(t_d = 6.3\)) for 80% silted-up dam (\(K=0.1\); \(K_R\) = 0.01)

Full size image

Smooth pressure fields are a key accomplishment of the proposed scheme and are essential to the successful functioning of the apparent viscosity equations. This is also linked to a smooth angular shear rate field, preventing massive premature yielding of the nonlinear phase because of numerical oscillations distributed through large portions of it. In the case of a granular phase, this is even more important since yield stress is a function of local pressure and even small fluctuations in the pressure field will be transferred to the viscosity field. This happens, for example, in the highly active area at the gate location seen in Fig. 18 that provides a sample of the viscosity field for an 80% silted-up dam.

Fig. 18
figure 18

Viscosity field at \(t = 1.0\) s (\(t_d = 6.3\)) for 80% silted-up dam (\(K=0.1\); \(K_R\) = 0.01)

Full size image

Water depth hydrograms are an indirect indicator of the amount of sand under that level, so nonlinear phase accuracy cannot be really evaluated unless it is compared directly with the evolving silt depth, which was available with the experimental datasets. The accuracy of the model for sand depth profiles was similar among different initial silt depth cases. The most demanding cases for the proposed method were 50%, 66%, and 80% initial silt depths, and even so, the results are very good. Figures 19 and 20 plot sand profiles for different significant instants for the cases with the least and the most sand. The latter case is not so important to evaluate accuracy but rather to show that even when a relatively small amount of nonlinear particles is present in the model, they will still exhibit proper dynamics.

Fig. 19
figure 19

Sand profile evolution for 80% silted-up dam

Full size image
Fig. 20
figure 20

Sand profile evolution for 10% silted-up dam

Full size image

The above water depth plots include three groups of results, generated using different values for numerical parameter K, which will act in accordance with threshold \(K \dot{\varepsilon } = 1\) from Eq. (12) to obtain the value of the local apparent viscosity. For this factor, many combinations of the threshold of \(K \dot{\varepsilon }\) and the value of K were tested, and in general terms, it was found that the scheme was stable for values \(K<0.5\) when the threshold value for \(K \dot{\varepsilon }\) was unitary. In this set of tests, whenever K became too small (say, \(K<0.001\)), it was more difficult to keep the scheme from becoming unstable as variations of apparent viscosity became sharper close to small values of angular shear rate \(\dot{\varepsilon }\). It was also observed that the higher the relative amount of nonlinear mass present in the model, the better a stronger value of K was for accuracy.

These considerations imply that the other numerical parameter, the relaxation factor \(K_R\), is constant and set at 0.01. This value, which produced good overall results, is safe to use for the problem types studied here, but that does not mean it will not interact negatively with other parts of the scheme. A sensitivity analysis carried out with the case with an 80% silted-up dam shows that the resulting hydrograms for relevant test values of \(K_R\) are plotted in Figs. 21, 22, 23 and 24 using a constant value of 0.01 for the apparent viscosity constant K.

Fig. 21
figure 21

Water depth at \(x/L = 0.5\) for 80% silted-up dam

Full size image
Fig. 22
figure 22

Water depth at \(x/L = 1.0\) for 80% silted-up dam

Full size image
Fig. 23
figure 23

Water depth at \(x/L = 1.7\) for 80% silted-up dam

Full size image
Fig. 24
figure 24

Water depth at \(x/L = 0.5\) for 66% silted-up dam

Full size image

Just as it occurs with the effect of varying K, changing \(K_R\) has a greater impact for longer simulation times and in areas where there was no fluid, like position x/L = 1.7. Other test cases with different initial silt depths replicated this general behaviour, where it seems that even if the changing shape of the silt bed is correctly replicated during the whole simulation time (as with this scheme) the waterfront is delayed. This cannot be completely blamed on artificial viscosity, since this was also the tendency during tests without activating that particular stabilization scheme. Once the main dam break event takes place, angular shear rate fields grow in mean value through the water phase from the dam that fell onto the wet bed section of the model.

Aside from these matters, values for \(K_R\) outside of a range of 0.1 to 0.005 resulted in inaccuracies because of premature yielding or a need for nonphysically high values of artificial viscosity to keep the scheme stable. Too low a value, in the vicinity of \(10^{-3}\), would not let the nonlinear phase yield and instead would favour particle suspension and mixing. Despite the fact that granular material can become temporarily suspended in flow, certain conditions must be met first, and one of the most important is enough speed; in these cases, even with low speeds of the suspending phase, with exceedingly low values of \(K_R\) nonlinear fluid particles became detached from the main mass and immediately lost nonlinear behaviour, so it was actually more of an unnatural boundary diffusion process rather than actual suspension. Another negative consequence of low values for \(K_R\) was the formation of unnatural cracks since the nonlinear mass does not yield to shear at the interface because of the constraint to the conservation equation. However, since the hydrostatic load still is present, any small oscillation ruptures the continuity and the nonlinear phase slowly collapses, with many growing cracks perpendicular to the boundary with little or no erosion processes present.

On the other hand, too high of a value for \(K_R\) can yield the nonlinear phase too soon because the conservation equation lets initial fluid yielding occur at a fast enough rate so that the threshold for angular shear rate is overcome. Since this rate is what feeds the calculation of viscous shear stresses, the yield criterion then will also be overcome without enough of a chance to act upon particle dynamics. Other problems arising when using high values for \(K_R\) include interaction with apparent viscosity constant K. If movement is difficult to start, the field of angular shear rates will remain small, and therefore apparent viscosity will grow across most of the nonlinear phase, driving down the value of the timestep required for numerical stability. With the values of physical constants associated with these test cases, this sometimes drove computational expense to impractical levels.

These constants were tuned using the test cases, but as the above discussion and diagrams show, within certain ranges, this scheme is not too sensitive to varying K or \(K_R\), so it is reasonable to apply these recommended values to other simulations in accordance to fluid type and the range of yield stress \(\tau _0\) where applicable. In summary, although it is enough to use values for \(K_R\) and K that are positive, nonzero and under unity as to not violate conservation laws, the set of tests used in this study resulted in the recommended values shown in Table 2.

Table 2 Relaxation and viscosity factor value guide
Full size table

Good versatility is expected since the test cases include violent interactions that come with any dam break. There are erosive phenomena, and the nonlinear phase did not end up artificially constrained but rather exhibited its own dam break dynamics. This variety of simultaneous types of behaviour has not been previously solved with this type of scheme that relies on an apparent viscosity scheme without recurring to assigning elastic behaviour to the nonyielded phase. This in turn is a satisfying development, since sand and silt in particular, and many real nonlinear fluids cannot exhibit elastic behaviour at all, so including it in a numerical scheme can introduce an entirely different type of inaccuracy.

In these multiphase cases, all numerical scheme components contributed to a relevant extent to stability and to smoothing data fields that are essential to maintaining unyielded areas of the fluid stable, but the most important change was the addition of the particle shifting scheme. This part of the method alone did not guarantee stability, but it certainly was a key factor for approaching a sufficiently versatile composition of a method for a multiphase, selectively yielding flow.

Another important contribution to good numerical behaviour was interface boundary treatment. For cases where high apparent viscosity was present close to the interface in the nonlinear phase, clumping and detachment of particle clusters were eradicated by updating state variables for particle i using data only from neighbouring particles that belong to the same phase as this particle. When applied to updating values of density, this prevented particle diffusion among phases.

The variety of conditions that the scheme had to react to in a stable manner without retuning numerical parameters means that it is robust enough to tackle many other multiphase cases and expect similar accuracy and replication of erosive and scouring phenomena. Although other approaches, like those using an equivalent elastic solid for the nonyielded phase, or the formulation of a specially formulated interfacial layer of particles along the yielded/nonyielded boundary, need to overcome similar challenges regarding stability and accuracy, this scheme provides comparable performance with a relatively simple assembly of proven SPH methods.

3.3 Method application insights

The scheme avoids using an elastic solid approximation for the unyielded portion of the nonlinear fluid. When dealing with nonconservative pseudosolids like a mass of granular material that piles up due to friction forces between grains, a nonelastic model not only is conceptually a better approximation but also will elude inaccuracies due to potential energy storage and exchange that certainly do not occur in real granular masses with a high elastic modulus. This is the case for soils, silt, and sand in a saturated state. Additionally, implementation is simpler as long as proper provisions are made to deal with smoothing pressure fields and applying the yield criteria at the right calculation stage.

Many studies have attempted the apparent viscosity approach to model nonlinear fluids, among which are viscoplastic fluids; these pose stiff challenges, like the modelling of the yielding process, that frequently causes researchers to implement algorithms that vary according to the target physical phenomena involved. In SPH schemes, using this intuitively direct concept has led to the consideration of additional factors to modify the value of apparent viscosity, such as weighted viscosity in mixed fluids, consideration of cohesive forces, or the addition of artificially high viscosities at low strain rates. Unfortunately, this approach alone has little effect on controlling unnatural pressure field variations and cannot deal with boundary effects that are harmful to the numerical stability of relevant variables. Successful yield point models often mean dealing with some sort of pseudosolid model acting according to a yield criterion, like an elastic solid model (SPH based or otherwise), a moving solid boundary with a particle addition/subtraction algorithm, or a particle refinement algorithm. Both elastic models with a high apparent modulus (like those mentioned in [15]) or a low apparent modulus (like that used in a colloidal suspension model in [53]) assume that the elastic energy will return with continuum relaxation, something that does not occur in granular flow where the discrete structure means dominance of friction forces in the dynamics. Some of these solutions involve algorithms that are complicated to implement, and many others, while relatively simple, fail to handle noisy data fields directly or cannot deal with a wide enough range of values for realistic rheological models.

Like most computational schemes, in this model some numerical parameters had to be added and tuned both for accuracy and stability, but a sensitivity analysis revealed that, for a certain family of cases, a recommended range can be set with confidence. Among additional numerical parameters to those pertaining to the weakly compressible SPH model is the apparent viscosity factor K, that proportionally affects apparent viscosity when angular shear rates are under a threshold as required by [38] in his approach to dealing with zero angular shear rates in any rheological model. Among other minor effects on data fields, higher values of K resulted in impractically small values for a stable timestep according to CFL criteria. This does not interfere with the proper modelling of fluid yielding, because high apparent viscosity is not used here to model the yielding point. Yielding processes are modelled using the selective application of the relaxation factor \(K_R\) to the conservation equation.

An essential finding of this study was that particle shifting based on the Fickian approach proposed by [42] is extremely important to ensure a smooth pressure field, and by extension, other data fields that are key to an accurate calculation of apparent viscosity and application of the yield criteria. If the particle shifting scheme is not used, in all cases studied pressure fields started to fluctuate in a manner that, locally, angular shear rates were high enough to yield the nonlinear phase, although average strain in an area about five to ten times the smoothing length was clearly below the yield threshold. Therefore, if even a small section of the nonlinear mass yielded, numerical oscillations artificially yielded the whole mass, and using a higher apparent viscosity, tuning conventional WCSPH numerical parameters or applying other artificial viscosity stabilization techniques did not contribute to accuracy. When attempted, either the fluid became hyperstatic or slightly slowed premature collapse, but never in a manner to achieve enough accuracy according to the physical experimental data. This became worse for the multiphase cases with values of yield stress \(\tau _0\) greater than \(10^2\) Pa, which are very common in sand and soil. However, particle shifting interferes with the modelling of erosion because it tends to keep particles spaced close enough so that particles expected to move with little to no effort, like grains of sand on a vertical surface, actually remained coherent with the unyielded mass as if it were composed of plastic clay. This meant that proper modelling of the yielding process could not be left only to the calculation of apparent viscosity, but also needed to be reapplied elsewhere.

All work done with apparent viscosity must deal with the impossibility of modelling nonzero yielding stress of a nonlinear fluid using viscosity alone. Since at zero strain rate, viscous stresses will be zero no matter the value of viscosity. Some approaches, such as forcing a minimum nonzero strain rate, usually lead to other difficulties, including incompatibility with higher values of yield stress and limited accuracy. Therefore, in this scheme, a modification of the yielded/unyielded treatment by [23] was applied, with the additional benefit of being able to introduce the yield criteria that will act after the application of the particle shifting calculation step, so scouring and erosion will be modelled properly.

When applying the yield criteria to the nonlinear phase, the direct approach to assign zero acceleration to the unyielded fluid particles was not suited to replicate the gradual collapsing of a mass, and corners exposed to erosion tended to destabilize pressure fields. For that reason, in the proposed scheme the nonyielded fluid still evolves with the conservation equation but was attenuated by the relaxation factor \(K_R\), another numerical parameter that required initial tuning. This provides a much smoother transition, numerically speaking, at the interface between different fluids, where high differentials in shear stresses are commonly encountered when modelling scouring and other erosion processes. Just as with the apparent viscosity factor K, a sensitivity analysis paved the way to a recommended range for \(K_R\) applicable to the group of cases studied here.

In this scheme, another crucial element for good numerical behaviour to prevent unnatural surface cracking or particle clumping was the selective updating of variables using information only from neighbouring particles belonging to the same phase. Nonphysical cracks and clumps can occur where apparent viscosity was high and pressure fields fluctuated. The solution was to omit data from particles of other phases or virtual boundary particles explicitly during physical variable updating processes. This same process is used as well when updating density to prevent diffusion among nonmixing fluid phases.

High differentials in shear stress magnitudes, such as those appearing at the boundary between different phases, seem to be among the main causes of noise in pressure fields. This noise cannot be controlled using artificial viscosity schemes without interfering with the very dynamics the viscosity model attempts to replicate. Other, less pronounced physical discontinuities like differences in density among phases, can also affect stability, but to a less significant extent.

4 Concluding remarks

A weakly compressible SPH scheme was developed successfully to model multi-phase, dam-break physics that can deal simultaneously with erosion and plastic collapse events, including cases where yield stresses can change according to local pressure. This requires the implementation of several methodological extensions since the smoothing kernel alone cannot keep data fields stable. The set of extensions used in this scheme, therefore, work together to guarantee sufficient numerical stability to replicate physically sound flow behaviour that can be experimentally validated for a certain set of cases.

The main contributions of this study were: (a) the integration of an SPH numerical scheme specifically designed to model a yielding/nonyielding granular material without approximating the unyielded mass to an elastic solid, avoiding the implied storage of elastic energy; (b) applying yield criteria to modify the response of the conservation equation, separately from the Particle Shifting Technique (PST) algorithm, to ensure proper modelling of scouring; (c) identifying a PST technique that, under a variety of erosive and flow dynamics, stabilizes pressure and strain fields to prevent irregular apparent viscosity fields and to replicate physical pressure distributions; and (d) validation of this straightforward approach with data from multiphase physical experiments that exhibit both nonlinear dam break collapse processes and erosive phenomena.

Although a good amount of attention was given here to the value of apparent viscosity and rheological models are certainly important to simulation accuracy, the development of this scheme was not centred on this variable, but on guaranteeing a smooth angular shear rate field, which does not depend upon a single element of the model or even the tuning of numerical parameters but on how all elements have been integrated. Future developments, before looking for completeness, must concentrate on keeping new elements that help accuracy from introducing noise to the strain rate fields, which was the central matter during the evolution of the model described in this study.