1 Introduction↩︎

Yield‑stress fluids are ubiquitous in natural and industrial contexts. Examples include food products such as ketchup, personal care products like toothpaste and gels, as well as natural materials like mud and clay. Such materials resist deformation up to a critical stress threshold, below which they behave like solids and above which they flow like viscous fluids. This behavior is generally attributed to microstructural rearrangements, with transitions between jammed and unjammed states [1], [2]. Recent reviews by [3]–[8] describe and discuss yield-stress fluids, and readers are referred to these excellent texts for further background on the topic.

The classical viscoplastic model attributed to Bingham [9] introduces a yield stress into a Newtonian framework. This is further extended by the Herschel–Bulkley model [10] which introduces an extra parameter controlling shear-thinning or shear-thickening. Another example of a yield-stress viscoplastic model is the \(\mu(I)\)-rheology, a now widely popular description of dense granular flows [11], [12], which has been implemented in various Navier-Stokes solvers (e.g., [13]–[15]). However, attributed in part to the non-smooth transition between flow and no flow by a sharp threshold, these well-known viscoplastic models should be regularized or augmented in numerical solvers [5], [16]–[19]. Moreover, such purely viscoplastic descriptions completely ignore elasticity, which not only governs the material response below yielding, but can also influence the flow dynamics (see, e.g., [20]–[22]). This recognition has motivated the development of elasto‑viscoplastic (EVP) constitutive models, which combine elastic deformation and rate‑dependent (viscous) plastic flow within a unified model.

In the non-Newtonian fluid mechanics community, Pierre Saramito [23]–[27] has made significant contributions to the development of modern EVP models. In his 2007 work, [23] introduced an EVP constitutive model that integrates elasticity with viscoplastic flow, providing a thermodynamically consistent description of yield‑stress fluids that behave as solids below yield and that flows above yield. Later, this was extended to incorporate Herschel–Bulkley viscoplasticity [24]. This model has been implemented in the finite element method [28] and more recently in the finite volume method [29]. It has been thoroughly benchmarked against alternative viscosity regularization methods and experimental data by [30], demonstrating great promise of this constitutive framework, although certain limitations have also been very recently identified [31], [32]. Subsequent developments explored brittle EVP behavior with a pressure‑sensitive yield surface [26], and very recent efforts, concurrent with the work presented here, have refined the thermodynamic environment of EVP models [27].

As strongly indicated by [33], yield-stress fluids are, at least in the slower limit, essentially plastic flow under large deformation, and thus the field of plasticity can benefit the constitutive and computational development of yield-stress fluids. This perspective suggests that concepts developed in finite strain plasticity can provide a natural framework for the modeling of such fluids. Although numerical modeling of materials with complicated yield criteria is well understood in solid mechanics, these established methods have not always been directly transferrable to flow problems. In granular flow modeling, some studies have introduced different yield surfaces into the \(\mu(I)\)-rheology in both a viscoplastic [15], [34], [35] and elasto-viscoplastic [36], [37] setting. Recognizing that the \(\mu(I)\)-rheology can be interpreted as a purely viscoplastic Drucker—Prager model with variable viscosity, [14] proposed a simplified viscoplastic Drucker—Prager formulation with constant viscosity, which also showed agreement with experimental observations of granular flows. In the context of landslide flow modeling, [38], [39] considered the combination of Perzyna viscoplasticity with different yield surfaces in an infinite EVP landslide model. This concept has been a source of motivation for the present work, in which we propose a general three-dimensional formulation and provide a verification and discussion of the resulting solutions, which has been outside the scope of the previous geophysical flows-focused studies.

Rather than focusing on a particular material rheology, the purpose of the present work is to compare and assess the challenges associated with three different yield surfaces, all well-known from solid mechanics, for the modeling of yield-stress flow problems. We address this by formulating a finite strain overstress‑type EVP model within a hybrid Eulerian–Lagrangian numerical model that can capture local transitions between solid- and fluid-like behavior, allowing two- and three‑dimensional simulations of free‑surface flows. In particular, the material point method (MPM) is used as the numerical scheme capable of handling large deformation and flow of EVP materials.

2 Theory↩︎

2.1 Finite strain elastoplasticity↩︎

We rely on a finite strain framework with a multiplicative decomposition of the deformation gradient \(\boldsymbol{F} = \boldsymbol{F}^E \boldsymbol{F}^P\) following [40]. From the deformation gradient, different strain measures can be defined, e.g., the Hencky strain tensor \(\boldsymbol{\varepsilon}\) defined as \[\boldsymbol{\varepsilon} = \frac{1}{2} \ln \boldsymbol{F} \boldsymbol{F}^T = \sum\limits_{\alpha=1}^{d} \ln \lambda_\alpha \text{ } \boldsymbol{n}_\alpha \otimes \boldsymbol{n}_\alpha \label{hencky}\tag{1}\] where \(d\) denotes the number of spatial dimensions, the unit vectors \(\boldsymbol{n}_\alpha\), \(\alpha = 1,...,d\), define the spatial principal directions and \(\lambda_\alpha\), \(\alpha = 1,...,d\), are the principal streches, i.e., the singular values of \(\boldsymbol{F}\). For simplicity, all models presented here will be based on Hencky’s elasticity model which relates the Kirchhoff stress \(\boldsymbol{\sigma} = \det(\boldsymbol{F}) \boldsymbol{\sigma}_*\), where \(\boldsymbol{\sigma}_*\) is the Cauchy stress, to the elastic Hencky strain \(\boldsymbol{\varepsilon}^E\) with a form familiar from isotropic linear elasticity, in particular, \[\boldsymbol{\sigma} = \mathcal{C} : \boldsymbol{\varepsilon}^E = \Lambda \mathop{\mathrm{tr}}\big(\boldsymbol{\varepsilon}^E\big) \boldsymbol{I} + 2 G \boldsymbol{\varepsilon}^E \text{,} \label{henckyelasticity}\tag{2}\] where \(\Lambda\) and \(G\) are the two Lamé parameters. These two parameters may be related to the perhaps more familiar Young’s modulus \(E = \frac{G(3\Lambda+2G)}{\Lambda+G}\) and Poisson’s ratio \(\nu = \frac{\Lambda}{2(\Lambda + G)}\), here assuming three dimensions. In addition, the bulk modulus can be defined as \(K=\Lambda + \frac{2}{d}G\). The properties of the finite strain hyperelastic model given by Eq. 2 and its relation to other models are given in, e.g., [41], [42].

A yield function \(y=y(\boldsymbol{\sigma})\) determines the onset of plastic deformations at \(y=0\), with elastic stress states satisfying \(y<0\). The velocity gradient can be additively decomposed as \[\nabla \boldsymbol{v} = \boldsymbol{l}^E + \boldsymbol{l}^P\] and a plastic flow rule imposes \[\boldsymbol{l^P} = \dot{\bar{\gamma}} \frac{\partial g}{\partial \boldsymbol{\sigma}} \label{plasticvelgrad}\tag{3}\] where \(g=g(\boldsymbol{\sigma})\) is known as a plastic potential, defining the direction of plastic flow, and the scalar \(\dot{\bar{\gamma}}\) is known as the plastic multiplier. The flow rule must induce a non-negative plastic rate of dissipation, i.e., \(\boldsymbol{\sigma} : \boldsymbol{l}^P \geq 0\). The choice \(g = y\) is called an associative flow rule and is a consequence of the principle of maximum plastic dissipation. In particular, it can be shown that maximizing \(\boldsymbol{\sigma} : \boldsymbol{l}^P\) subject to \(y \leq 0\) leads to \(\boldsymbol{l^P} = \dot{\bar{\gamma}} \frac{\partial y}{\partial \boldsymbol{\sigma}}\) where \(\dot{\bar{\gamma}}\) can be interpreted as a Lagrange multiplier [43], [44].

In this work, \(y\) and \(g\) will be chosen as functions of the principal stresses \[p = -\frac{1}{d} \mathop{\mathrm{tr}}(\boldsymbol{\sigma)} \text{ \;and \;} \tau = \frac{1}{\sqrt{2}} || \mathop{\mathrm{dev}}(\boldsymbol{\sigma}) || \label{p95and95q}\tag{4}\] where \(p\) is termed the isotropic pressure and \(\tau\) the equivalent shear stress. Here, we used the notation \(|| \boldsymbol{A} || = \sqrt{\boldsymbol{A} : \boldsymbol{A}} = \sqrt{\mathop{\mathrm{tr}}(\boldsymbol{AA})}\) and \(\mathop{\mathrm{dev}}(\boldsymbol{A}) = \boldsymbol{A} - \frac{1}{d}\mathop{\mathrm{tr}}(\boldsymbol{A})\boldsymbol{I}\) for a symmetric second-order tensor \(\boldsymbol{A}\). The particular form of Eq. 4 is chosen such that the two-dimensional case \(d=2\) provides plane strain conditions with the assumption \(\nu=0.5\) for yield and plastic flow, see, e.g., [45], [46] and references therein.

In the remainder of this section, three plastic models are outlined; von Mises, Drucker–Prager and modified Cam–clay. Their yield function can all be written on the general form \[y(p,\tau) = \tau - \tau_y(p) \label{yield}\tag{5}\] where \(\tau_y(p)\) is a yield stress that may depend on \(p\). The different yield surfaces \(y=0\) are sketched in Fig. 1.

In the von Mises model [47] the yield strength does not depend on pressure, i.e., \[\tau_y(p) = \tau_c \label{VM}\tag{6}\] where \(\tau_c>0\) is a constant (assuming no hardening).

The Drucker–Prager model [48] yield stress given by \[\tau_y(p) = \mu p + \tau_c \label{DP}\tag{7}\] where \(\mu>0\) is the internal friction and \(\tau_c \geq 0\) is a cohesive stress. With an associative flow rule, \(g=y\), this would, for any plastic deformation, imply a volume increase, i.e., \(\mathop{\mathrm{tr}}(\boldsymbol{l}^P) > 0\). Consequently, this model is often used with a von Mises plastic potential, i.e., \(g=\tau\).

The modified Cam–clay model [49] is characterized by an evolving yield surface of ellipsoidal shape in principal stress space, in particular, \[\tau_y(p) =\mu \sqrt{ (p-p_t) (p_c-p) } \text{,} \label{MCC}\tag{8}\] where \(\mu>0\) defines the so-called critical state line, \(p_c \geq 0\) is an isotropic compressive strength and \(p_t \leq 0\) is an isotropic tensile strength. Following [50], we define \(\beta \geq 0\) as a dimensionless measure of cohesion such that \(p_t = -\beta p_c\). Based on critical state soil mechanics [51], this model is conceived from the idea that there exists a critical state where the material can continue to shear without further changes to its volume or stress. This model is usually used with an associative flow rule, thus promoting the possibility of both dilation and compaction (until critical state is reached) and a hardening law based on the accumulated plastic volumetric deformation \(\varepsilon_V^P = \int \mathop{\mathrm{tr}}(\boldsymbol{l}^P) \;dt\). As such, the critical state must be the stress state at the apex of the ellipsoidal yield surface, which will always be along the critical state line defined by \(\mu\) (see Fig. 1c). The hardening law considered here is given by \[p_c(\varepsilon_V^P) = p_c^0 e^{-\xi \varepsilon_V^P} \label{hardening95exp}\tag{9}\] where \(\xi>0\) is a hardening parameter and \(p_c^0\) is an initial compressive strength [37].

2.2 The overstress approach↩︎

In rate-independent elastoplasticity, the Kuhn–Tucker conditions must be satisfied, i.e., \(\dot{\bar{\gamma}} \geq 0\), \(y \leq 0\) and \(\dot{\bar{\gamma}} y = 0\). This ensures that plastic states can only exist on the yield surface \(y=0\). In the overstress approach to viscoplasticity, plastic states may satisfy \(y > 0\) and these conditions no longer hold. Prominently, the overstress model of [52], [53] can be written in a general way as \[\boldsymbol{l}^P = \begin{cases} \frac{\Psi(y)}{t_v} \boldsymbol{n} , & \text{if } y > 0 \\ \boldsymbol{0}, & \text{otherwise} \end{cases} \label{Perzyna95general}\tag{10}\] where \(\Psi(y)\) is a scalar overstress function and \(\boldsymbol{n} = \frac{\partial g / \partial \boldsymbol{\sigma}}{|| \partial g / \partial \boldsymbol{\sigma} ||}\) determines the flow rule. The overstress approach stands in contrast to, e.g., the consistency approach [54], where, instead of relying on the rate-independent yield function \(y(\boldsymbol{\sigma})\), one formulates a rate-dependent yield function such that the Kuhn–Tucker conditions can still be enforced. Interestingly, this can be shown to be essentially equivalent to the Perzyna approach [55].

A popular special case of Eq. 10 is found by considering a von Mises plastic potential, i.e., \(g=\tau\) and yield surfaces described by a (potentially pressure-dependent) shear stress \(\tau_y(p)\) as outlined in Section 2.1, \[\dot{\gamma} = \begin{cases} \frac{1}{t_v} \Big( \frac{\tau-\tau_y}{\tau_y} \Big)^{1/s}, & \text{if } \tau > \tau_y \\ 0, & \text{otherwise} \end{cases} \label{Perzyna}\tag{11}\] where \(\dot{\gamma} = \sqrt{2} || \mathop{\mathrm{dev}}(\boldsymbol{l}^P)||\) is denoted the equivalent plastic shear strain rate, following the decomposition \(\boldsymbol{l}^P = \frac{1}{d} \mathop{\mathrm{tr}}(\boldsymbol{l}^P) \boldsymbol{I} + || \mathop{\mathrm{dev}}(\boldsymbol{l}^P)|| \frac{\mathop{\mathrm{dev}}(\boldsymbol{\sigma}) }{|| \mathop{\mathrm{dev}}(\boldsymbol{\sigma}) ||}\). In Eq. 11 , the two viscous parameters are the viscous time \(t_v\) and the viscous exponent \(s\).

[56] proposed the following slight modification, \[\dot{\gamma} = \begin{cases} \frac{1}{t_v} \Big( \Big(\frac{\tau}{\tau_y} \Big)^{1/s} -1 \Big), & \text{if } \tau > \tau_y \\ 0, & \text{otherwise} \end{cases} \label{Peric}\tag{12}\] which, unlike the previous formulation, recovers the rate-independent solution as \(s \rightarrow 0\) in addition to \(t_v \rightarrow 0\). This will be shown in Section 3.1

The overstress model typically attributed to [57] is often given by the general form \(\boldsymbol{\sigma} - \mathcal{P}\boldsymbol{\sigma} = t_v \;\mathcal{C} : \boldsymbol{l}^P\) where \(\mathcal{P}\boldsymbol{\sigma}\) is a projection of the stress onto the yield surface [58]. Typically, the latter is a closest point projection in stress space in the sense of an associative plastic flow rule. This Duvaut–Lions model can be rewritten as \[\boldsymbol{l}^P = \begin{cases} \frac{1}{t_v} \;\mathcal{C}^{-1} : (\boldsymbol{\sigma}-\mathcal{P}\boldsymbol{\sigma}), & \text{if } y > 0 \\ \boldsymbol{0}, & \text{otherwise} \end{cases} \label{DL}\tag{13}\] In the special case of a von Mises plastic potential and yield surfaces which can be described by a (pressure-dependent) shear stress \(\tau_y(p)\) as outlined in Section 2.1, this can be reduced to \[\dot{\gamma} = \begin{cases} \frac{1}{t_v} \frac{\tau-\tau_y}{G}, & \text{if } \tau > \tau_y \\ 0, & \text{otherwise} \end{cases} \label{DL95G}\tag{14}\] which, when compared to Eq. 11 , motivates a generalized expression \[\dot{\gamma} = \begin{cases} \frac{1}{t_v} \Big( \frac{\tau-\tau_y}{G} \Big)^{1/s}, & \text{if } \tau > \tau_y \\ 0, & \text{otherwise} \end{cases} \label{DL95G95s}\tag{15}\] or, even more generally, a generalized Duvaut–Lions formulation can be proposed as \[\boldsymbol{l}^P = \begin{cases} \frac{1}{t_v} \;\big( \mathcal{C}^{-1} : (\boldsymbol{\sigma}-\mathcal{P}\boldsymbol{\sigma}) \big)^{1/s}, & \text{if } y > 0 \\ \boldsymbol{0}, & \text{otherwise} \end{cases} \label{DL95general}\tag{16}\] This demonstrates the equivalence between the Perzyna and Duvaut–Lions formulations. From the point of view of Perzyna, if 1) the normalization constant (i.e., the denominator) in the overstress function \(\Psi\) is chosen appropriately (e.g., in the example above as \(G\) instead of \(\tau_y\)) and 2) the plastic potential \(g\) is chosen in line with the projection operator \(\mathcal{P}\), they are equivalent.

2.3 Steady-state velocity profiles↩︎

We consider here the flow on an inclined plane, as sketched in Fig. 2, in order to derive steady-state velocity profiles for Peric and Duvaut–Lions fluids with both von Mises and Drucker–Prager yield with a von Mises plastic potential. Assuming \(\sigma_{xx} = \sigma_{zz}\), we have \(\tau = \rho g (h-z)\sin \theta\) and \(p = \rho g (h-z)\cos \theta\). In the figure, \(h_y\) denotes the height until the start of the plug flow, and can be found as the elevation \(z=h_y\) where \(\tau = \tau_y\). In particular, \[h_y = h-\frac{\tau_c}{\rho g \omega}\] where we for convenience defined \[\omega = \begin{cases} \sin \theta, & \text{if von Mises} \\ \sin \theta - \mu \cos \theta , & \text{if Drucker--Prager} \end{cases}\] where we must assume \(\mu > \tan(\theta)\) for flow. Assuming \(\dot{\gamma} = \frac{\partial v_x(z)}{\partial z}\), the steady-state velocity profiles \(v_x(z)\) resulting from the von Mises and Drucker–Prager EVP models with a von Mises plastic potential can now be derived.

Considering first the Duvaut–Lions model, Eq. 15 , we obtain \[\frac{\partial v_x}{\partial z}(z) = \frac{1}{t_v} \Big( \frac{\rho g \omega}{G} \Big)^{1/s}(h_y-z)^{1/s}\] which can be integrated to give \[v_x(z)= \frac{1}{t_v} \frac{s}{1+s} \Big( \frac{\rho g \omega}{G} \Big)^{1/s} \Big( h_y^{\frac{1+s}{s}}-\big(h_y-z\big)^{\frac{1+s}{s}} \Big) \text{.} \label{ssDL}\tag{17}\] Remarkably, in the special case \(s=2\), the above result has a striking resemblance to the Bagnold velocity profile which can be derived from the \(\mu(I)\)-rheology, \[v_x(z) \propto \Big( h_y^{3/2}-\big(h_y-z\big)^{3/2} \Big) \label{bagnold}\tag{18}\] with a prefactor which depends on the material properties (density and grain size) and rheological parameters (that define the \(\mu(I)\) function) in addition to the angle \(\theta\). This is derived by, e.g., [59] using the expression for \(\mu(I)\) from [60].

With the Peric model, Eq. 12 , we obtain the following general expression for the shear strain rate \[\frac{\partial v_x}{\partial z}(z) = \frac{1}{t_v} \bigg( \Big(\frac{(h-z)\sin \theta}{\mu(h-z)\cos \theta + (h-h_y)\omega}\Big)^{1/s} - 1 \bigg) \label{velprofile95peric95generalderivative}\tag{19}\] where \(\mu=0\) must be replaced in the case of von Mises. In the case of a von Mises yield, this can be integrated to give, \[v_x(z)= \frac{1}{t_v} \bigg( \frac{s}{1+s} \frac{ h^{\frac{1+s}{s}} - (h-z)^{\frac{1+s}{s}} }{(h-h_y)^{1/s}} - z \bigg) \label{velprofile95peric95vm}\tag{20}\] where the first term in the parentheses shares a similar \(z\)-dependence as Eq. 17 . However, with a non-cohesive Drucker–Prager yield, we obtain a linear velocity profile, regardless of the exponent \(s\), \[v_x(z)= \frac{1}{t_v} \bigg( \Big( \frac{\tan \theta}{\mu} \Big)^{1/s} - 1\bigg)z \label{velprofile95peric95dp}\tag{21}\] Interestingly, the velocity profiles in Eqs. 20 and 21 do not depend on any material properties other than the viscous time and viscous exponent.

3 Numerical implementation↩︎

3.1 Elastic predictor – plastic corrector↩︎

Following the approach of [44] in elastoplasticity, an elastic predictor – plastic corrector scheme is employed. Denoting the Hencky strain \(\boldsymbol{\varepsilon}^{E, t}\) predicted in an elastic step from time \(T^n\), this is later corrected to \(\boldsymbol{\varepsilon}^{E, n+1}\) at time \(T^{n+1}\) as \[\boldsymbol{\varepsilon}^{E, n+1} = \boldsymbol{\varepsilon}^{E, t} - \boldsymbol{l}^P \Delta t \label{predictorcorrector}\tag{22}\] where \(\Delta t =T^{n+1}-T^n\) and \(\boldsymbol{l}^P\) is given in Eq. 3 . The derivation of Eq. 22 follows [61] and is given in detail in Appendix A of [37].

Equivalently, we can consider Eq. 22 in stress space as \(\boldsymbol{\sigma}^{n+1} = \boldsymbol{\sigma}^t - \mathcal{C} : \boldsymbol{l}^P \Delta t\). In the context of the Duvaut–Lions model, Eq. 13 , this can be used to immediately derive \[\boldsymbol{\varepsilon}^{E, n+1} = \boldsymbol{\varepsilon}^{E, t} - b \dot{\bar{\gamma}} \frac{\partial g}{\partial \boldsymbol{\sigma}} \Delta t \label{predictorcorrector95DL}\tag{23}\] where we defined \(b = \frac{1}{1+t_v/\Delta t}\) and \(\dot{\bar{\gamma}}\) is the plastic multiplier found in the equivalent rate-independent (inviscid) case. Note that when \(t_v = 0\), we get \(b=1\) and Eq. 23 reduces to the rate-independent solution as expected.

With a von Mises plastic potential, Eq. 22 can be reduced to \(\tau^{n+1} = \tau^t - G \dot{\gamma} \Delta t\) and \(p^{n+1} = p^t\). Combined with the Peric model, Eq. 12 , this gives \[(\tau^t - G \dot{\gamma} \Delta t)\bigg(\frac{1}{t_v \dot{\gamma} +1}\bigg)^s - \tau_y = 0\] which can be solved for \(\dot{\gamma}\) with iterative schemes, e.g., the Newton–Raphson method. Note that both \(t_v=0\) and \(s=0\) give \(\tau^{n+1} = \tau_y\) which is the rate-independent solution. In contrast, Eqs. 11 and 15 do not provide the rate-independent solution in the limit \(s\rightarrow 0\). For example, with Eq. 15 , the equation \[\tau^t - G \dot{\gamma} \Delta t - G(t_v \dot{\gamma})^s - \tau_y = 0\] must be solved for \(\dot{\gamma}\).

3.2 The volume increase in Drucker–Prager↩︎

As mentioned in Section 2.1, the Drucker–Prager model is often used with a von Mises plastic potential in order to avoid excessive volume increase. However, the non-smoothness of this yield surface at \(\tau=0\) necessitates a special plastic flow rule ensuring that, in the rate-independent case, the plastic state ends up on the yield surface, commonly, at the Drucker–Prager cone tip. This will induce some volumetric expansion, which may not be significant under small deformations, but may become excessive under large deformations (as will be shown in Section 4.3). To this end, a plastic volume correction technique similar to the one proposed by [62] can be utilized. This technique is based on tracking a temporary plastic volumetric strain \(\varepsilon_V^{P,*}\), which is initially zero, but can have both positive and negative values. The (rate-independent) return mapping can now be written \[\boldsymbol{\varepsilon}^{E, n+1} = \begin{cases} \frac{\tau_c}{dK\mu} \boldsymbol{I}, \text{ if } \mathop{\mathrm{tr}}(\boldsymbol{\varepsilon}^{E,t}) > \frac{\tau_c}{K\mu} - \varepsilon_V^{P,*,n} \\ \boldsymbol{\varepsilon}^{E, t} - \Delta t \dot{\bar{\gamma}} \boldsymbol{n} + \frac{1}{d} \varepsilon_{V}^{P, *,n} \boldsymbol{I}, \text{\;otherwise} \end{cases}\] where \[\begin{align} \varepsilon_{V}^{P,*, n+1} = \begin{cases} \varepsilon_{V}^{P,*, n} + \Delta \varepsilon_{V}^{P} , \text{ if } \mathop{\mathrm{tr}}(\boldsymbol{\varepsilon}^{E,t}) > \frac{\tau_c}{K\mu} - \varepsilon_V^{P,*,n} \\ 0, \text{\;otherwise} \end{cases} \raisetag{1.5\baselineskip} \end{align}\] i.e., \(\varepsilon_{V}^{P,*}\) is reset to zero once we are no longer in a plastic dilating state. This technique can be interpreted as letting some states undergo volume reduction by projecting them to the cone tip instead of to the yield surface through the given flow rule, consequently compensating for the otherwise strict volume expansion. This is sketched in Fig. 3, highlighting the plastic flow under the case \(\varepsilon_V^{P,*,n}>0\) and the case \(\varepsilon_V^{P,*,n}<0\).

3.3 The material point method↩︎

The constitutive laws have been implemented in the material point method (MPM), a numerical method for solving the momentum conservation equation \[\rho(\boldsymbol{x},t) \frac{D \boldsymbol{v}(\boldsymbol{x},t)}{D t} = \nabla \cdot \boldsymbol{\sigma}(\boldsymbol{x},t) + \rho(\boldsymbol{x},t) \boldsymbol{g} \text{,} \label{mom95cons95eulerian}\tag{24}\] where \(\rho\) is the mass density, \(\boldsymbol{v}\) is the velocity, \(\boldsymbol{\sigma}\) is the symmetric Cauchy stress tensor, \(D/Dt\) denotes the material time derivative and we have assumed that the only external body force present is the one resulting from the gravitational acceleration \(\boldsymbol{g}\). Typically attributed to [63], MPM is capable of handling large deformations and flow without mesh distortion issues and requiring no special treatment for the free surface. In this method, the material is discretized by Lagrangian (material) points associated with (at least) a mass \(m_p\), volume \(V_p\), velocity \(\boldsymbol{v}_p\) and strain \(\boldsymbol{\varepsilon}_p\). At every time step, the mass and velocities are interpolated to a background grid, and the momentum equation is solved on this grid.

As we are interested in flow dynamics, we focus here on an explicit MPM formulation. We adopt a regular grid with a grid cell size \(\Delta x\), initialized with eight material points per grid cell in three dimensions. Disregarding the boundary term, the discretized weak form of the momentum conservation equation can be written on the Eulerian grid nodes \(i\) for each time step \(n\) of size \(\Delta t\) as, assuming a forward Euler time discretization, \[\begin{align} \frac{ m^n_{i} \boldsymbol{v}^{n+1}_{i} - m^n_{i} \boldsymbol{v}^{n}_{i} }{\Delta t} &= - \sum_p V^0_{p} \boldsymbol{\sigma}(\boldsymbol{\varepsilon}_p^{E,n}) \nabla N_{i}(\boldsymbol{x}^n_p) \nonumber \\ &\quad + m_i^n \boldsymbol{g} \label{wf95discr} \end{align}\tag{25}\] where \(\boldsymbol{g}\) specifies gravity, \(N_{i}\) is an interpolation function and the mass on grid node \(i\) is \(m_i^n = \sum_p m_p N_i(\boldsymbol{x}_p^n)\). Since \(m_p\) remains constant in time, mass conservation is inherently satisfied. We let \(\Delta t\) be bounded by the elastic wave speed as well as the CFL condition, i.e., \[\Delta t = \min \Bigg( C_\text{cfl} \frac{\Delta x}{\max_p(|\boldsymbol{v}_p|)}, \text{ \;} C_\text{el} \frac{\Delta x}{\sqrt{E/\rho}} \Bigg)\] with the positive constants \(C_\text{cfl} \leq 1\) and \(C_\text{el} \leq 1\). Boundary conditions are imposed directly on the updated grid velocities, enforcing either no-slip or frictional slip governed by a prescribed Coulomb friction.

Initial critiques of MPM raised concerns about cell-crossing instabilities and numerical dissipation. In recent years, the MPM community has tackled these problems through the use of B-spline interpolation functions [64] as well as the adoption of affine transfer schemes such APIC [65] and subsequently AFLIP [66] to reduce dissipation. Further details about MPM and the implementation used in this work can be found in [67].

Although we have opted for MPM in this study, other implementations of the same constitutive framework can be accomplished in other methods with similar features, such as the Smoothed Particle Hydrodynamics (SPH), the Finite Element Method with Lagrangian Integration Points (FEMLIP) and the Particle Finite Element Method (PFEM), which all have been used for viscoplastic simulations of yield-stress fluids (e.g., [68]–[71]). Each method comes with its own advantages and disadvantages. For example, MPM demands extra memory due to its background grid, FEMLIP requires a treatment of the free surface, SPH is associated with expensive particle neighbor search and challenging handling of boundary condition and in PFEM frequent remeshing is necessary.

4 Examples↩︎

4.1 Steady-state flow↩︎

With the setup sketched in Fig. 2, using periodic boundary conditions in the \(x\)-direction, we investigate whether the numerical implementation can retrieve the steady-state velocity profiles found in Section 2.3. Using the Peric model, Eq. 12 , and the Duvaut–Lions model, Eq. 15 , this is shown in Figs. 4 and 5, respectively, using both a von Mises yield (a-c) and Drucker–Prager yield (d-e). In all cases the elastic parameters \(E=1\) MPa and \(\nu=0.3\) are used, the inclination angle is \(\theta=40^\circ\) and the velocity profile is measured at \(t=0.5\) s after starting from rest.

There is no general analytical solution for the steady-state velocity profile with a modified Cam–clay yield. However, with a sufficiently small initial value of the initial compressive strength \(p_c\) and a sufficiently large value of the hardening parameter \(\xi\), we propose that the plastic states will reside around the critical state line determined by \(\mu\), and, as such, we expect the solution to be equivalent to that of Drucker–Prager. Figs. 6a and d show that, indeed, the Drucker–Prager solution can be retrieved for various values of \(t_v\) and \(\mu\). Moreover, Figs. 6b and c show the departure from this solution when \(p_c\) and \(\beta\) are increased, providing a control on the emergence of plugged surface flow.

While the overstress approach ensures that steady-state flow can be induced, the purely inviscid model (\(t_v = 0\)) does not produce steady-state flow, rather it will accelerate indefinitely if \(\mu < \tan \theta\) or come to rest if \(\mu > \tan \theta\). As such, the overstress approach, similar to the \(\mu(I)\)-rheology, provides a means of regularizing this instability of the inviscid Drucker–Prager model.

4.2 Dam break↩︎

The preceding example demonstrated that modified Cam–clay can reproduce the solution with Drucker–Prager for infinite shear flow. To investigate this beyond this idealized setting, we consider the case of a dam break released on a no-slip surface, inclined an angle \(\theta = 10^\circ\) from the horizontal. Similar setups, with various inclinations, has been used by various researchers to study the behavior of yield-stress fluids, e.g., [72]–[74].

With the choice of an internal friction \(\mu > \tan(\theta)\), the dam break will experience a finite runout distance \(x_f\) along the inclined plane. This runout decreases as a function of viscosity, e.g., as illustrated in the example provided in Fig. 7 using the viscous modified Cam–clay model. In fact, with comparable parameters, there is no significant difference between the Drucker–Prager and modified Cam–clay models in the runout or collapse of the dam break, as shown in Fig. 8. This figure also includes the corrected Drucker–Prager model from Section 3.2.

Fig. 9 shows that the total energy in the system is approximately conserved by the numerical scheme (Fig. 9a) and that the small energy loss decreases with a finer spatial resolution (Fig. 9b). The runout distance also converges with a finer resolution, as shown in Fig. 9c. In Fig. 9a, the various contributions to the total energy are included, showing that the loss of potential energy is mainly converted into plastic dissipation, highlighting the minor contribution of the kinetic energy in the initial collapse.

The corrected Drucker–Prager model provides a similar result as the standard Drucker–Prager model, with only a very slightly reduced runout, which can be explained by the minimal volumetric expansion occurring in this example. Therefore, in the next example, we consider a case that displays a stronger volumetric deformation.

4.3 Stretching flows↩︎

In order to further demonstrate the differences between the corrected and uncorrected Drucker–Prager model, as well as the modified Cam–clay model, which all carry the notion of an internal friction parameter \(\mu\), we consider here a material released on an inclined plane as in the previous dam break simulations, here with an inclination of \(\theta = 30^\circ\). This flow is subsequently stretched from falling off the end of this plane, dropping vertically before impacting the ground, consequently inducing a folding of the material. The simulation with Drucker–Prager is illustrated in Fig. 10. We consider a frictional boundary condition, for simplicity with a Coulomb friction \(\mu_c=\mu=\tan(\theta)\), and we allow the material to separate from the boundary. Simply looking at the final rest state of the piled up material from the simulations using Drucker–Prager, the corrected Drucker–Prager and modified Cam–clay, we see in Fig. 11 that the Drucker–Prager has gained a significant amount of volume. In particular, Drucker–Prager induces an increase of the volume of the flowing material by around 50%. Surprisingly, the corrected Drucker–Prager and modified Cam–clay give almost identical volumetric behavior, increasing slightly during the stretching phase (until about \(t=0.9\) s after release), before compacting as the material impacts the ground and starts folding.

5 Discussion and conclusions↩︎

In this work, we have demonstrated, through representative examples, how the choice of yield surface (von Mises, Drucker–Prager and modified Cam–clay) influences the behavior of EVP fluids. Rate-dependent plasticity was considered through the overstress approach of Perzyna and Duvaut–Lions, between which we have highlighted their differences as well as their general equivalence under certain model choices. Moreover, we analytically derived velocity profiles for idealized shear flows and showed that the Bagnold velocity profile for granular flow can be considered a special case of the Duvaut–Lions model. We chose MPM as the numerical scheme in which we implemented the proposed finite strain framework, however, it must be emphasized that this framework could also be implemented in other schemes. For example, [69] provides a different EVP framework in FEMLIP, which, like the framework proposed here, does not require any regularization of the discontinuity due to the solid-fluid transition.

Drucker–Prager should generally not be accompanied by an associative flow rule due to excessive volume increase. As demonstrated, even with a von Mises plastic potential (i.e., a non-associative flow rule), Drucker–Prager can cause excessive volume increase during large deformations, due to the special treatment of the yield surface singularity (typically treated with a cone tip projection). While this volumetric expansion is typically negligible in small-strain solid mechanics, it must be carefully considered when modeling fluid flows. To this end, we have shown that the correction of [62], initially proposed in the computer graphics community, can help reduce excessive volume increase.

Alternatively, the (viscous) modified Cam–clay model can be used to retrieve the solution of the (viscous) Drucker–Prager, without any problems of excessive volume increase. This was shown for shear flow and dam break problems. However, it is important to highlight that Cam–clay models can also be used (with other parameter choices) to produce very different results than Drucker–Prager and, unlike Drucker–Prager, has great potential in modeling compressible flows. The compressible behavior can be specified through the hardening law, for which we here chose a relatively simple relationship given by Eq. 9 . We also restricted ourselves to relatively large values of the hardening parameter \(\xi\) in order to induce a hardening behavior that promoted stress states to rather quickly align with the critical state. This was a deliberate choice in order to critically compare modified Cam–clay with Drucker–Prager. The compressible properties of Cam–clay, especially for solid behavior, are beyond the scope of this contribution, and the reader is referred to classical geomechanics literature (e.g., [75], [76] and others) for the general characteristics of this model.

We have not studied the influence of elastic properties and we have restricted ourselves to one (hyper)elastic model. The chosen elastic model, Hencky’s elasticity model, has desirable features which are further discussed in [41], [42] and makes the numerical implementation particularly easy. Although the plastic properties are usually dominant during flow, the elastic properties may nevertheless be highly relevant and may also be rate-dependent. While we have considered yield surfaces that broadly represent the pressure-independent, frictional and fully capped families, our selection is not exhaustive. For example, the Tresca model, which is pressure-independent like von Mises, the Mohr–Coulomb model, which bears similarities to Drucker–Prager, the capped Drucker–Prager model [77], which extends Drucker–Prager to yield also under compression similar to modified Cam–clay, as well as more modern models such as Matsuoka—Nakai [78], [79] and Nor-sand [80] may also be considered. The (solid) properties of these models may be found in classical geomechanics literature.

Acknowledgments↩︎

The authors thank Johan Gaume and Guillaume Chambon for fruitful discussions. L.B. gratefully acknowledges financial support from the Swiss National Science Foundation (SNSF) through grant number P500PT_230265.

Data availability statement↩︎

All models presented here are implemented and freely available in the open-source MPM software Matter [67], [81].

CRediT authorship contribution statement↩︎

Lars Blatny: Conceptualization of this study, Methodology, Software, Validation, Formal analysis, Investigation, Writing - Original Draft, Visualization, Funding acquisition. Alexandre Pellet: Methodology, Software, Writing - Review & Editing.

References↩︎