*Article* **Implementation of the Non-Associated Elastoplastic MSDPu Model in FLAC3D and Application for Stress Analysis of Backfilled Stopes**

**Feitao Zeng, Li Li \*, Michel Aubertin and Richard Simon**

Research Institute on Mines and Environment, Department of Civil, Geological and Mining Engineering, École Polytechnique de Montréal, C.P. 6079 Succursale Centre-Ville, Montréal, QC H3C 3A7, Canada; feitao.zeng@polymtl.ca (F.Z.); michel.aubertin@polymtl.ca (M.A.); richard.simon@polymtl.ca (R.S.) **\*** Correspondence: li.li@polymtl.ca

**Abstract:** The multiaxial Mises-Schleicher and Drucker-Prager unified (MSDPu) criterion has been shown to exhibit several specific features compared to other yield and failure criteria, including a nonlinear mean stress dependency, influence of the Lode angle, use of independent uniaxial compressive and tensile strength values and absence of an apex (singularity) on the envelope surface in the negative stress quadrant. However, MSDPu has been seldom used in practice to solve geotechnical and geomechanical engineering problems mainly because it had not yet been fully implemented into three-dimensional (3D) numerical codes. To fill this gap, a 3D elastoplastic MSDPu formulation is developed and implemented into FLAC3D. The proposed MSDPu elastic-perfectly plastic (EPP) constitutive model is then validated against existing analytical solutions developed for calculating the stress and displacement distributions around cylindrical openings. The FLAC3D MSDPu-EPP model is then applied to evaluate the vertical and horizontal stress distributions in a three-dimensional vertical backfilled stope. The numerical results obtained with the MSDPu-EPP model are compared with those obtained with the Mohr-Coulomb EPP model, to highlight key features of the new formulation.

**Keywords:** 3D nonlinear yield criterion; elastoplastic model; numerical modeling; circular opening; backfill; FLAC3D

#### **1. Introduction**

Elastoplastic constitutive models are widely used in geotechnical engineering to assess the mechanical response of geomaterials. The elastoplastic framework typically involves a yield criterion, a flow rule (with a plastic potential) and, in some cases, a hardening or softening function. Over the years, many elastoplastic constitutive models have been proposed and applied to analyze the complex mechanical responses of rocks and soils; the main ones are included in state-of-the-art review publications [1–9].

Elastic-perfectly plastic (EPP) models, with a fixed yield surface, are probably the most often used in practice to solve engineering problems involving geomaterials, due in a large part to their relative simplicity and ease of application. Several EPP models with different plastic (yield) criteria, such as Mohr-Coulomb (MC), Drucker-Prager (DP) and Hoek-Brown (HB), have already been built in commercial codes and applied in geotechnical engineering. In these models, the plastic criterion *F* is used to determine the limit of the stress state associated with plastic behavior. Criterion *F* also defines the plastic potential *Q* when an associated flow rule is used, while it can serve as a basis for *Q* (=*F*) in the case of a non-associated flow rule [9,10].

In addition to those mentioned above, a large number (>100) of failure and yield criteria have been proposed over the years [11–14], each having its advantages and limitations. The most commonly used criteria for soil and rock, MC and HB, are generally expressed with only two principal stresses and thus neglect the effect of the intermediate principal

**Citation:** Zeng, F.; Li, L.; Aubertin, M.; Simon, R. Implementation of the Non-Associated Elastoplastic MSDPu Model in FLAC3D and Application for Stress Analysis of Backfilled Stopes. *Processes* **2022**, *10*, 1130. https://doi.org/10.3390/pr10061130

Academic Editor: Blaž Likozar

Received: 15 May 2022 Accepted: 2 June 2022 Published: 5 June 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

stress. A few well-known 3D criteria (e.g., DP) simplify the effect of the stress geometry by neglecting the influence of the Lode angle (defined below). The open surface defined by function *F* in the principal stress space along the hydrostatic axis is another limitation of many existing criteria, which cannot describe the volumetric yield behavior of geomaterials under high mean stresses [14–16].

Initially proposed for intact rocks [17,18] and later modified and extended for various geomaterials including rock and rock mass, concrete and mine backfill [19–22], the multiaxial Mises-Schleicher and Drucker-Prager unified (MSDPu) criterion takes into account the effect of the three principal stresses, with a nonlinear (rounded) surface on the negative side of the mean stress axis and a cap on the positive (compressive) side. The MSDPu criterion exhibits four essential characteristics to define yielding or failure of cohesive/cemented geomaterials:


Additional features of MSDPu will be presented in the next section, after recalling the criterion basic formulation. Its applicability to describe the failure and yielding of a large variety of materials has been demonstrated by Aubertin et al. [20,21], Li et al. [14,22,23] and Aubertin and Li [12].

Despite its advantageous features, the practical use of the MSDPu criterion has been limited because it was only partially implemented in a numerical code. The implementation has been done within the 2D finite difference code FLAC through its user-defined model option with an external language, called "FISH" [24]. The ensuing FLAC2D MSDPu-EPP model can be used to analyze geotechnical problems under plane strain conditions only. Some of the key features and advantages of MSDPu thus cannot be exploited, particularly when facing three-dimensional problems. In addition, the implementation of the MSDPu-EPP model in FLAC2D (presented by Li et al. [24]) was based on an associated flow rule, which tends to overestimate the volumetric strains of geomaterial (and the related mean stresses). Moreover, the implementation of user-defined models in FLAC with the "FISH" language is no longer recommended by Itasca [25]. There is thus a need to implement a more complete 3D version of the MSDPu-EPP model, with a non-associated flow rule, in FLAC3D, a commercial code widely applied for the solution of three-dimensional problems in geotechnical engineering [26–30].

In this paper, the three-dimensional MSDPu-EPP model is formulated in terms of stress and strain increments, following the guidelines provided by Itasca [26] as well as Desai and Siwardane [1] and Chen and Baladi [2]. Its implementation in FLAC3D is done through its user-defined model option with the programing language C++. The model is then compiled into a DLL (dynamic link library) module that can be loaded and run as a built-in plug-in of FLAC3D. The FLAC3D MSDPu-EPP model is partially validated against some existing analytical solutions developed for evaluating the stresses and displacements distributions around cylindrical openings in plane strain. The new 3D model is applied to evaluate the vertical and horizontal stress distributions in three-dimensional backfilled stopes, which are compared to those obtained with the commonly used Mohr-Coulomb EPP (MC-EPP) model. The results comparison highlights some of the new model's distinctive features such as a better fit to the yield (failure) surface for a wide range of geomaterials and a more representative volumetric yield behavior at relatively high mean stress.

#### **2. The MSDPu Nonlinear Multiaxial Criterion**

The general form of a plastic (yield) criterion for isotropic materials can be expressed in terms of commonly used stress invariants [1,31,32]:

$$F(I\_1, I\_2, I\_3) = const \tag{1}$$

where *I*<sup>1</sup> is the is the first invariant of the stress tensor *σij*; *J*<sup>2</sup> and *J*<sup>3</sup> are respectively the second and third invariants of the deviator stress tensor *Sij* = *σij*–*pδij*; *δij* is the Kronecker delta (*δij* = 1 if *i* = *j* and *δij* = 0 if *i* = *j*); the mean stress *p* is defined as follows:

$$p = \frac{I\_1}{3} \tag{2}$$

The three invariants *I*1, *J*<sup>2</sup> and *J*<sup>3</sup> can be expressed using the general stress tensor components [1,8], or in terms of the major (*σ*1), intermediate (*σ*2) and minor (*σ*3) principal stresses:

$$I\_1 = \sigma\_1 + \sigma\_2 + \sigma\_3 \tag{3}$$

$$J\_2 = \frac{1}{6} \left[ \left( \sigma\_1 - \sigma\_2 \right)^2 + \left( \sigma\_2 - \sigma\_3 \right)^2 + \left( \sigma\_3 - \sigma\_1 \right)^2 \right] \tag{4}$$

$$J\_3 = (\sigma\_1 - p)(\sigma\_2 - p)(\sigma\_3 - p) \tag{5}$$

The MSDPu criterion can then be written as

$$F = f\_2 - F\_0^2 F\_\pi^2 = 0 \tag{6}$$

where *F*<sup>0</sup> is associated with the nonlinear surface for conventional triaxial compression condition (*σ*<sup>1</sup> > *σ*<sup>2</sup> = *σ*3), while *F<sup>π</sup>* defines the surface in the *π*-plane. These two functions are usually expressed as follows [21]:

$$F\_0 = \left[a^2 \left(I\_1^2 - 2a\_1 I\_1\right) + a\_2^2 - a\_3 \left(I\_1 - I\_c\right)^2\right]^{1/2} \tag{7}$$

$$F\_{\pi} = \frac{b}{\left[b^2 + (1 - b^2)\sin^2(45^\circ - 1.5\theta)\right]^{1/2}}\tag{8}$$

The Macaulay brackets 〈*x*〉 (= (*x* + |*x*|)/2, where *x* is a variable) are used in Equation (7) to avoid a negative term; *α* (taken from the DP criterion), *a*<sup>1</sup> and *a*<sup>2</sup> are material parameters related to shear failure or yielding:

$$\alpha = \frac{2\sin\phi}{\sqrt{3}(3-\sin\phi)}\tag{9}$$

$$a\_1 = \frac{\mathcal{C}\_0 - T\_0}{2} - \frac{\mathcal{C}\_0^2 - \left(\frac{T\_0}{\theta}\right)^2}{6a^2(\mathcal{C}\_0 + T\_0)}\tag{10}$$

$$a\_2 = \left\{ \left[ \frac{\mathbb{C}\_0 + \left( \frac{T\_0}{b^2} \right)}{\mathfrak{Z}(\mathbb{C}\_0 + T\_0)} - a^2 \right] \mathbb{C}\_0 T\_0 \right\}^{1/2} \tag{11}$$

where *φ* is internal friction angle, *C*<sup>0</sup> is uniaxial compressive strength (UCS), *T*<sup>0</sup> is uniaxial tensile strength (UTS, positive value) and *b* is a shape parameter defining the surface in the *π*-plane. Material parameter *a*<sup>3</sup> is related to the volumetric yield (cap) surface:

$$a\_3 = \frac{a^2 \left(I\_{1n}^2 - 2a\_1 I\_{1n}\right) + a\_2^2}{\left(I\_{1n} - I\_c\right)^2} \tag{12}$$

where *Ic* is the *I*<sup>1</sup> value where the cap starts and *I*1*<sup>n</sup>* is *I*<sup>1</sup> value where it intersects the hydrostatic axis (*σ*<sup>1</sup> = *σ*<sup>2</sup> = *σ*3) in the principal stress space.

The Lode angle *θ* is defined as

$$\theta = \frac{1}{3} \sin^{-1} \frac{3\sqrt{3}l\_3}{2\sqrt{l\_2^3}}, -30^\circ \le \theta \le 30^\circ \tag{13}$$

From this definition, one can deduce *θ* = 30◦ for conventional triaxial (axisymmetric) compression (CTC, with *σ*<sup>1</sup> > *σ*<sup>2</sup> = *σ*3) and *θ* = −30◦ for reduced triaxial (axisymmetric) extension (RTE, with *σ*<sup>1</sup> = *σ*<sup>2</sup> > *σ*3). Under CTC testing conditions, *θ* = 30◦ and the commonly defined deviatoric stress *q* = (3*J*2) 1/2 <sup>=</sup> *<sup>σ</sup>*<sup>1</sup> − *<sup>σ</sup>*3.

Figure 1 shows a typical representation of the MSDPu criterion in the *I*1–*J*<sup>2</sup> 1/2 plane (Figure 1a) and *π*-plane (Figure 1b). In the *I*1–*J*<sup>2</sup> 1/2 plane, the MSDPu envelope can be decomposed into two main parts. In the first part when *I*<sup>1</sup> ≤ *Ic*, *J*<sup>2</sup> 1/2 increases nonlinearly with *I*1. In the second part when *I*<sup>1</sup> > *Ic*, the cap controls the yield surface and *J*<sup>2</sup> 1/2 tends to decrease with an increasing *I*1. For dense geomaterials such as hard rocks, the values of *Ic* and *I*1*<sup>n</sup>* can be very high so the cap can be neglected for most applications; this is not the case however for porous materials such as soils, backfills and some weak rocks. In the *π*-plane, the usual shape of the MSDPu envelope defined above by Equation (8) depends on the value of parameter *b*. The typical value of *b* goes from 1 to about 0.7 and the corresponding envelope evolves from a circle (for *b* = 1) to a rounded triangle (*b* < 1) in the *π*-plane, as shown in Figure 1b.

**Figure 1.** MSDPu failure surface in the (**a**) *I*1–*J*<sup>2</sup> 1/2 plane and (**b**) *π*-plane (*σ*<sup>1</sup> \*, *σ*<sup>2</sup> \*, *σ*<sup>3</sup> \* are the projections of the principal stress axes in the *π*-plane).

The main characteristics of the MSDPu criterion can be summarized as follows:


Additional details can be found in Aubertin et al. [20,21], Aubertin and Li [12] and Li et al. [23].

#### **3. Implementation of the Multiaxial MSDPu-EPP Model in FLAC3D**

The constitutive model formulation in FLAC3D can be expressed in terms of the principal stresses *σ*1, *σ*2, *σ*<sup>3</sup> and the principal strains *ε*1, *ε*2, *ε*3, using procedures described by Itasca [25]. The total strain increment Δ*ε<sup>i</sup>* can be divided into elastic Δ*ε<sup>i</sup> <sup>e</sup>* and plastic Δ*ε<sup>i</sup> p* strain increments (with *i* = 1, 2, 3)

$$
\Delta \varepsilon\_i = \Delta \varepsilon\_i^{\varepsilon} + \Delta \varepsilon\_i^p \tag{14}
$$

The principal stress increments associated with the elastic strain increments can then be expressed as follows:

$$\begin{array}{l} \Delta \sigma\_1 = S\_1 \left( \Delta \boldsymbol{\varepsilon}\_1^{\varepsilon}, \Delta \boldsymbol{\varepsilon}\_2^{\varepsilon}, \Delta \boldsymbol{\varepsilon}\_3^{\varepsilon} \right) = \beta\_1 \Delta \boldsymbol{\varepsilon}\_1^{\varepsilon} + \beta\_2 \left( \Delta \boldsymbol{\varepsilon}\_2^{\varepsilon} + \Delta \boldsymbol{\varepsilon}\_3^{\varepsilon} \right) \\\ \Delta \sigma\_2 = S\_2 \left( \Delta \boldsymbol{\varepsilon}\_1^{\varepsilon}, \Delta \boldsymbol{\varepsilon}\_2^{\varepsilon}, \Delta \boldsymbol{\varepsilon}\_3^{\varepsilon} \right) = \beta\_1 \Delta \boldsymbol{\varepsilon}\_2^{\varepsilon} + \beta\_2 \left( \Delta \boldsymbol{\varepsilon}\_1^{\varepsilon} + \Delta \boldsymbol{\varepsilon}\_3^{\varepsilon} \right) \\\ \Delta \sigma\_3 = S\_3 \left( \Delta \boldsymbol{\varepsilon}\_1^{\varepsilon}, \Delta \boldsymbol{\varepsilon}\_2^{\varepsilon}, \Delta \boldsymbol{\varepsilon}\_3^{\varepsilon} \right) = \beta\_1 \Delta \boldsymbol{\varepsilon}\_3^{\varepsilon} + \beta\_2 \left( \Delta \boldsymbol{\varepsilon}\_1^{\varepsilon} + \Delta \boldsymbol{\varepsilon}\_2^{\varepsilon} \right) \end{array} \tag{15}$$

where *S*1, *S*2, *S*<sup>3</sup> are linear functions for the Hooke's law; *β*<sup>1</sup> and *β*<sup>2</sup> are material parameters (constants) defined in terms of the isotropic shear modulus *G* and bulk modulus *K*,

$$\begin{array}{l} \beta\_1 = K + \frac{4}{3}G\\ \beta\_2 = K - \frac{2}{3}G \end{array} \tag{16}$$

The plastic strain increment is given by the flow rule

$$
\Delta \boldsymbol{\varepsilon}\_i^p = \lambda \frac{\partial \mathcal{Q}}{\partial \sigma\_i} \tag{17}
$$

where *σ<sup>i</sup>* is the current stress state component (or initial stress state); *λ* is plastic coefficient; *Q* is the plastic potential (function), defined as follows

$$Q = l\_2 - \zeta F\_0^2 F\_\pi^2 \tag{18}$$

where coefficient *ξ* serves to control the plastic deviatoric and volumetric strain components. When *ξ* = 1, the plastic potential function is the same as the yield criterion (*Q* = *F*) so the flow rule is associated. When *ξ* 1 (very small value), the plastic potential leads to quasi-isovolumetric plastic strains, typical of critical state [15,36].

Combining Equations (14) and (17) leads to

$$
\Delta \varepsilon\_i^\varepsilon = \Delta \varepsilon\_i - \lambda \frac{\partial Q}{\partial \sigma\_i} \tag{19}
$$

Substituting Equation (19) into Equation (15) gives

$$
\Delta \sigma\_{\bar{i}} = S\_{\bar{i}}(\Delta \varepsilon\_1, \Delta \varepsilon\_2, \Delta \varepsilon\_3) - \lambda S\_{\bar{i}}(\frac{\partial Q}{\partial \sigma\_1}, \frac{\partial Q}{\partial \sigma\_2}, \frac{\partial Q}{\partial \sigma\_3}) \tag{20}
$$

Under varying stress conditions, the new stress state *σ<sup>i</sup> <sup>N</sup>* corresponding to the total strain increment Δ*ε<sup>i</sup>* is expressed by

$$
\sigma\_i^N = \sigma\_i + \Delta \sigma\_i \tag{21}
$$

In the user-defined model, the induced (postulated) elastic stresses *σ<sup>i</sup> <sup>I</sup>* are obtained by adding a "trial" stress increment to the current (initial) stress state *σi*. The "trial" stress increment is computed by using the incremental form of Hooke's law and the total strain increment Δ*εi*. The induced elastic stresses *σ<sup>i</sup> <sup>I</sup>* can be expressed as follows:

$$
\sigma\_i^l = \sigma\_i + S\_i(\Delta \varepsilon\_1, \Delta \varepsilon\_2, \Delta \varepsilon\_3) \tag{22}
$$

Substituting Equations (20) and (22) into Equation (21) leads to

$$
\sigma\_i^N = \sigma\_i^I - \lambda S\_i \left( \frac{\partial Q}{\partial \sigma\_1}, \frac{\partial Q}{\partial \sigma\_2}, \frac{\partial Q}{\partial \sigma\_3} \right) \tag{23}
$$

Equation (23) is called the plastic correction in FLAC3D.

The induced elastic stresses are initially taken as the new stress state and then adjusted if required. When the new stress state is within the elastic domain (i.e., *F* < 0), the new stress state is updated by directly using the incremental expression of Hooke's law. When the new stress state would lead to *F >* 0, the plastic coefficient *λ* is calculated from the MSDPu yield criterion to bring the new stress state *σ<sup>i</sup> <sup>N</sup>* determined by Equation (23) on the yield surface (*F* = 0), which results in

$$F\left(\sigma\_1^N, \sigma\_2^N, \sigma\_3^N\right) = f\_2^N - \left(F\_0^N\right)^2 \left(F\_\pi^N\right)^2 = 0\tag{24}$$

Substituting Equation (23) into Equation (24) gives the quadratic equation:

$$A\lambda^2 + B\lambda + \mathbb{C} = 0\tag{25}$$

The correct value of *λ* corresponds to the root with smaller absolute value of the two obtained after solving Equation (25).

Since *F*<sup>0</sup> includes two pieces and depends on the value of *I*1, the expressions of *A*, *B* and *C* in Equation (25) depend on whether (or not) *I*<sup>1</sup> ≤ *Ic*. Hereafter, their expressions are introduced with "*s*" denoting a shear response (for *I*<sup>1</sup> ≤ *Ic*) and "*v*" denoting the contribution of a volumetric response related to the cap (*I*<sup>1</sup> > *Ic*).

When *I*<sup>1</sup> ≤ *Ic*, *F*<sup>0</sup> becomes

$$F\_0^s = \left[a^2 \left(I\_1^2 - 2a\_1 I\_1\right) + a\_2^2\right]^{1/2} \tag{26}$$

The corresponding plastic potential is then expressed as

$$Q^s = \mathcal{J}\_2 - \tilde{\mathcal{J}} (F\_0^s)^2 F\_\pi^2 \tag{27}$$

The expressions for *A<sup>s</sup>* , *Bs* and *C<sup>s</sup>* are given in Appendix A. When *I*<sup>1</sup> > *Ic*, *F*<sup>0</sup> is written as

$$F\_0^v = \left[ a^2 \left( I\_1^2 - 2a\_1 I\_1 \right) + a\_2^2 - a\_3 \left( I\_1 - I\_3 \right)^2 \right]^{1/2} \tag{28}$$

and

$$Q^v = f\_2 - \zeta (F\_0^v)^2 F\_\pi^2 \tag{29}$$

The expressions for *Av*, *Bv* and *C<sup>v</sup>* are given in Appendix A.

To simplify the calculations of partial derivatives, variation of Lode angle *θ* is postulated to have a negligible effect on *F<sup>π</sup>* to obtain the new (updated) stress state (from the induced elastic stresses). This simplification leads to the following function in the *π*-plane:

$$F\_{\pi}^{N} \approx F\_{\pi}^{I} = \frac{b}{\left[b^{2} + (1 - b^{2})\sin^{2}(45^{\circ} - 1.5\theta^{I})\right]^{1/2}} \tag{30}$$

with

$$\theta^N \approx \theta^I = \frac{1}{3} \sin^{-1} \frac{3\sqrt{3}J^I\_3}{2\sqrt{\left(J^I\_2\right)^3}}\tag{31}$$

The computational procedure of the MSDPu-EPP model in FLAC3D v6.0 (Itasca, 2017) starts with adding the stress components, which are computed from the incremental

Hooke's law by using the total strain increments to the current (initial) stress state *σ<sup>i</sup>* and the induced elastic stresses *σ<sup>i</sup> <sup>I</sup>* are obtained (Equation (22)). Then, the *σ<sup>i</sup> <sup>I</sup>* values are substituted into the yield function *F* to determine if these are in the elastic domain (*F* < 0) or the plastic domain (*F* ≥ 0). If in the elastic domain, the new (updated) stress state *σ<sup>i</sup> <sup>N</sup>* equals to *σ<sup>i</sup> I* . If in the plastic domain, the new stress state *σ<sup>i</sup> <sup>N</sup>* is updated by applying the plastic correction (Equation (23)) to *σ<sup>i</sup> I* . The "corrected" new principal stresses *σ*<sup>1</sup> *<sup>N</sup>*, *σ*<sup>2</sup> *<sup>N</sup>* and *σ*<sup>3</sup> *<sup>N</sup>* can then be used to update the stress tensor *σij* in the system of reference axes, assuming that the principal directions have not been affected by the occurrence of a plastic correction.

Figure 2 shows the computational scheme for the implementation of the non-associated MSDPu-EPP model in FLAC3D 6.0 [25]. As indicated above, FLAC3D v6.0 provides the option to load and run user-written model in DLL (dynamic link library). The implementation of the MSDPu-EPP model was hence performed by creating a user-written DLL, which was created by compiling a program written with C++ language in Visual Studio 2017.

**Figure 2.** Computational scheme for the implementation of the non-associated MSDPu-EPP model in FLAC3D.

#### **4. Validation of the FLAC3D MSDPu-EPP Model**

After its implementation in FLAC3D, preliminary simulations with the user-defined MSDPu-EPP model were conducted to verify that the formulation and code programming were correctly done and to validate (in part) the numerical results. The main results from this assessment are summarized here.

The first validation is made against the analytical solution for the stresses around a cylindrical opening in a MSDPu-EPP material developed by Li et al. [23] for the open MSDPu surface (without cap) and the solution of Li and Aubertin [37] with the closed MSDPu yield surface (with cap). It should be noted that the analytical solutions were developed for a constant Lode angle *θ* = 0◦, even though the actual Lode angle tends to vary slightly in the plastic region (between about −19◦ and −25◦) along radial coordinate *r* around the cylindrical opening, while *θ* = 0◦ in the elastic region [23,37]. I

Figure 3a shows a cylindrical opening in an elastic-perfectly plastic material subjected to a hydrostatic far-field stress *P*0; an internal pressure *p*<sup>0</sup> is applied to the wall of the cylindrical opening (in plane strain condition). In the figure, *r*<sup>0</sup> is the radius of the opening, *R* is the radius of the interface between the plastic and elastic regions, *σ<sup>r</sup>* and *σψ* are the radial and tangential stresses respectively and *r* and *ψ* are the corresponding cylindrical coordinates. The calculations are made for *r*<sup>0</sup> = 1 m, *P*<sup>0</sup> = 30 MPa and *p*<sup>0</sup> = 2 MPa. The MSDPu-EPP model parameters are given in Table 1.

**Figure 3.** (**a**) A cylindrical opening under far-field isotropic stresses outside and an internal pressure applied to the inside wall, (**b**) a numerical model built with FLAC3D after taking advantage of the symmetry along two axes.

**Table 1.** Parameters used for comparing the results obtained with the FLAC3D MSDPu-EPP model and analytical solutions (taken from Li and Aubertin [37]).


Figure 3b shows a numerical model with domain size *d* × *d* and the boundary conditions. The model built with FLAC3D takes advantage of the quarter-symmetry geometry. The analytical solutions were developed for a plane-strain condition, so the numerical simulations are conducted in 2D even if the MSDPu-EPP model implemented in FLAC3D is three-dimensional. A commonly used method to do that is to isolate a thin domain in the direction perpendicular to the axis, with the displacements fixed in the direction parallel to the opening axis but allowed in the directions perpendicular to the axis. The thickness *t* of the modeling domain is taken as one-fifth of the cylinder radius, i.e., *t* = 0.2 m.

A sensitivity analysis was conducted to obtain an optimal configuration of the numerical model, which is based on an optimal domain size and optimal mesh. The optimal domain corresponds to the smallest size of the model to minimize the time of calculation and large enough to ensure stable and reliable numerical results. Similarly, an optimal mesh size *m* is associated with the coarsest elements (blocs) to minimize the time of calculation, with elements that are fine enough to ensure stable and reliable numerical results.

Figure 4 shows the variation of the radial displacements (Figure 4a) and stresses (Figure 4b) at a reference point *M* (shown in Figure 2) as a function of the mesh size (for domain size *d* = 10 m). In this figure, the mesh size is defined from the minimum size of the first layer of the grid around the opening, which increases at a constant ratio from the opening wall to the domain boundary (based on radial meshing; Itasca [25]). The results shown in Figure 4 indicate that, for both versions of the MSDPu-EPP model (with and without the cap), the numerical results tend to become stable once the mesh size is equal to or smaller than 0.02 m. The optimal mesh size corresponds to *m* = 0.02 m.

**Figure 4.** Variation of the (**a**) radial displacements and (**b**) stresses at point *M* (see Figure 3b) as a function of the mesh size (with *d* = 10 m), obtained from the MSDPu-EPP model with and without the cap.

Figure 5 shows the variation of the radial displacements (Figure 5a) and stresses (Figure 5b) as a function of the domain size (with *m* = 0.02 m). It can be seen that *d* = 10 m is the smallest value, which remains large enough to ensure stable and reliable numerical results. The optimal numerical model is then constructed with *m* = 0.02 m and *d* = 10 m, which are used to conduct the numerical simulations.

**Figure 5.** Variation of the (**a**) radial displacements and (**b**) stresses at point *M* (see Figure 3b) as a function of domain size (with *m* = 0.02 m) obtained from the MSDPu-EPP model with and without the cap.

Figure 6 shows the radial (*σr*) and tangential (*σψ*) stress distributions obtained from the numerical modeling and analytical solution for the cases of a closed MSDPu surface (Figure 6a) and an open MSDPu surface (Figure 6b). The agreements between the numerical and analytical results are excellent in all cases. Additional calculations were also made for other material parameters and loading conditions (results not shown here) to help assess the validity of the MSDPu-EPP model and its numerical implementation in FLAC3D. The good agreements obtained between the numerical and analytical results support the validation and indicate that the analytical solutions developed by considering a constant Lode angle provide a good approximation of the stress distribution around a cylindrical opening.

**Figure 6.** Distributions of the radial and tangential stresses, obtained from the numerical simulations and analytical solutions developed by Li and Aubertin [37] (**a**) for a closed MSDPu surface (with cap) and (**b**) for an open MSDPu surface (no cap).

The second component of this validation of the FLAC3D MSDPu-EPP model is performed against the analytical solution of Salençon [38,39], developed for evaluating the distribution of stresses and radial displacements around a cylindrical opening in a MC-EPP material (considering an associated flow rule). Parameter *ξ* = 1 (associated flow rule) and *a*<sup>3</sup> = 0 (without cap) are thus taken for the FLAC3D MSDPu-EPP model. As the MSDPu is nonlinear and MC is linear in the *I*1–*J*<sup>2</sup> 1/2 plane, the MC strength parameters *φ* and *c* were chosen so the two yield surfaces would correlate to each other (as well as possible) in the stress domain of interest around the cylindrical opening (for 50 MPa ≤ *I*<sup>1</sup> ≤ 100 MPa in this case). The resulting MC strength parameters are *φ* = 32◦ and *c* = 3.9 MPa for matching the open MSDPu surface with the parameters given in Table 1. The two yield surfaces in the *I*1–*J*<sup>2</sup> 1/2 plane (for *θ* = 0◦) are very close to each other for the stress states of interest, as shown in Figure 7. A constant Lode angle *θ* = 0◦ is again considered an acceptable approximation for the stresses around the cylindrical opening.

**Figure 7.** The open MSDPu yield surface (no cap) and corresponding Mohr-Coulomb surface in *I*1–*J*<sup>2</sup> 1/2 plane (for *θ* = 0◦).

Figure 8 illustrates the distributions of the stresses (Figure 8a) and radial displacements (Figure 8b), obtained from FLAC3D MSDPu-EPP model and calculated with the analytical solutions of Salençon [38,39]. The good agreement between the two types of solutions indicates that the FLAC3D MSDPu-EPP model correctly represents the problem at hand. This (partial) validation supports the use of the FLAC3D MSDPu-EPP model to analyze the elastoplastic response of geomaterials around underground openings.

**Figure 8.** Distributions of the (**a**) stresses and (**b**) radial displacements, obtained from the numerical simulation with FLAC3D MSDPu-EPP model and calculated with the analytical solutions of Salençon [38,39].

#### **5. Applications of the FLAC3D MSDPu-EPP Model**

In order to further assess the features of the new FLAC3D MSDPu-EPP model, simulations are conducted to analyze the mechanical response of cemented backfill in a 3D vertical mine stope.

Backfill is commonly used in underground mine excavations to ensure safe working conditions and improve ore recovery. Underground backfilling is also gaining momentum as a mine waste management approach [40]. Evaluating the stresses in backfilled stopes is of great importance and interest for underground mine stability. Until recently however, most of the numerical analyses were done under plane strain (2D) conditions. Only a few investigations have included 3D numerical analyses with the Mohr-Coulomb EPP model [27–29,41–43].

Figure 9 shows the conceptual model of a three-dimensional vertical mine stope in a semi-infinite rock mass before (Figure 9a) and after (Figure 9b) the addition of backfill. This stope is located at a depth of 500 m below the ground surface. It is excavated to a height of 45 m; the width is 6 m in the two horizontal directions (Figure 9a). After excavation, the stope is filled by a cemented (cohesive) backfill to a final height of 44.5 m, leaving a void space of 0.5 m at the top (Figure 9b). The stress distribution in the backfill and surround rock mass then depends on the fill settlement under its own weight and interaction with the rock walls. The stability of the rock walls of the empty stope is first evaluated using both MSDPu-EPP and MC-EPP models, respectively.

Figure 10 shows the corresponding numerical model of the stope, built by FLAC3D after excavation, before backfilling. Half of the stope is modelled to take advantage of the symmetric geometry. As indicated in the figure, the boundary conditions applied to the model domain are defined with a free surface on the top boundary and a fixed surface at the base. The conditions imposed along the four external vertical surfaces allow displacements within their respective plane, but not in the perpendicular direction. The rock mass is considered homogeneous, isotropic and perfectly elastoplastic; it is characterized by the following properties: *Er* = 30 GPa (Young's modulus), *ν<sup>r</sup>* = 0.3 (Poisson's ratio), *γ<sup>r</sup>* = 27 kN/m3 (unit weight). These parameters are used for both of the MC-EPP (with a tension cut-off at *T*0) and MSDPu-EPP models. The strength parameters are *φ* = 32◦, *c* = 3.9 MPa and *T*<sup>0</sup> = 0.2 MPa for the MC criterion and *φ* = 27◦, *C*<sup>0</sup> = 7 MPa, *T*<sup>0</sup> = 0.2 MPa and *a*<sup>3</sup> = 0 for the open MSDPu criterion. The rock mass is subjected to three-dimensional in-situ stresses, including the vertical in-situ stress *σ<sup>v</sup>* (*z*-direction) generated by the overburden (i.e., *σ<sup>v</sup>* = *γrh*; *h* is the depth below ground surface), maximum horizontal in-situ stress in the *x*-direction *σ<sup>H</sup>* taken as two times the vertical in-situ stress (i.e., *σ<sup>H</sup>* = 2*σz*) and minimum horizontal in-situ stress in the *y*-direction *σ<sup>h</sup>* (varying values in the different simulations).

**Figure 9.** Conceptual model of a three-dimensional vertical stope in rock mass (**a**) before and (**b**) after being backfilled.

**Figure 10.** A three-dimensional numerical model of the vertical stope in rock mass built by FLAC3D.

Sensitivity analyses are performed to ensure stable and reliable numerical results and obtain optimal numerical model for each case. For a numerical model such as the one shown in Figure 10, the sensitivity analysis gives the optimal mesh *m* and optimal domain size *d*.

The procedure is illustrated in Figure 11, which shows the variation of vertical displacement (*z*-direction) and horizontal stress *σxx* at a reference point *N* (at the center of the stope base, see Figs. 9a and 10) as a function of the mesh size (with *d* = 500 m) for simulations conducted with the MSDPu-EPP model, for an in-situ stress state *σ<sup>h</sup>* = *σ<sup>H</sup>* = 2*σz*. The results indicate that stable numerical results are achieved with *m* = 0.5 m. Figure 12 shows the variation of the vertical displacement (*z*-direction) and *σxx* as a function of the domain size (with *m* = 0.5 m). It is seen that the numerical results become stable when the

domain size *d* is equal to or larger than 500 m. The results thus indicate that the optimal numerical model has a mesh size of 0.5 m and a domain size of 500 m.

**Figure 11.** Variation of the (**a**) vertical displacement (*z*-direction) and (**b**) horizontal stress *σxx* at point *N* (see Figures 9a and 10) as a function of the mesh size *m* (for *d* = 500 m).

**Figure 12.** Variation of the (**a**) vertical displacement (*z*-direction) and (**b**) horizontal stress *σxx* at point *N* (see Figures 9a and 10) as a function of the domain size *d* (for *m* = 0.5 m).

Figure 13 shows the yielded areas around the stope along a vertical cross cut section in the *xz* symmetric plane, obtained from the numerical simulations conducted with the MSDPu-EPP and MC-EPP models, respectively. The results have been obtained for *σ<sup>h</sup>* equals to 1, 1.4 and 2 times *σv*, while *σ<sup>v</sup>* and *σ<sup>H</sup>* are kept unchanged. It can be seen that the size of yielded area obtained with the MSDPu-EPP model significantly increases as the intermediate in-situ stress *σ<sup>h</sup>* increases from *σ<sup>v</sup>* to 2*σv*, indicating that the intermediate in-situ stress plays an important role in the response and stability of openings. With the MC-EPP model however, the yielded area stays almost unchanged when *σ<sup>h</sup>* increases from *σ<sup>v</sup>* to 2*σ<sup>v</sup>* because the MC criterion is based on the 2D formulation and thus fails to consider explicitly the effect of the intermediate principal stress in the 3D extension used in FLAC3D.

**Figure 13.** Yielded areas around the stope along a vertical cross cut section in the *xz* symmetric plane, obtained from the numerical simulations with the MSDPu-EPP and MC-EPP models, respectively, considering (**a**) *σ<sup>h</sup>* = *σv*, (**b**) *σ<sup>h</sup>* = 1.4*σ<sup>v</sup>* and (**c**) *σ<sup>h</sup>* = 2*σv*.

Figure 14 shows the corresponding numerical model after the placement of backfill in the stope, built with FLAC3D. As was the case for the numerical model shown in Figure 10, half of the stope is modelled to take advantage of the symmetric geometry. The boundary conditions are defined by a free top boundary surface, a fixed bottom boundary surface and four vertical external surfaces whose displacements are allowed within their respective plane, but not allowed in the direction perpendicular to that plane. The stope backfilling is performed in multiple layers after the convergence of the rock walls has been completed. The thickness of each layer is 5 m (except for the top layer having a thickness of 4.5 m). A mesh size *m* of 0.5 m is used for the backfill; the domain size *d* = 500 m, determined for the rock wall stability analysis shown in Figure 10, is also used here. As the backfill is placed in the open stope after all the elastic and plastic strains have been released, the mechanical response of the backfill is almost independent on the in-situ stresses and rock model. The rock mass is thus considered homogeneous, isotropic and linearly elastic (without yield), characterized by *Er* = 30 GPa, *ν<sup>r</sup>* = 0.3 and *γ<sup>r</sup>* = 27 kN/m3. The vertical in-situ stress *σ<sup>v</sup>* is generated by the overburden and the horizontal in-situ stress is isotropic, with *σ<sup>h</sup>* = *σ<sup>H</sup>* = 2*σz*.

The stresses in the backfilled stope are evaluated using the MSDPu-EPP and MC-EPP models. The weakly cemented backfill is characterized by *Eb* = 300 MPa (Young's modulus), *ν<sup>b</sup>* = 0.327 (Poisson's ratio), ρ*<sup>b</sup>* = 1800 kg/m3 (density). The MSDPu criterion is defined by *φ* = 30◦, *C*<sup>0</sup> = 10 kPa, *T*<sup>0</sup> = 0.2 kPa, *b* = 0.75, *Ic* = 100 kPa and *a*<sup>3</sup> = 0.06. A value of *ξ* = 0.01 is used to express the plastic potential *Q* for the non-associated flow rule to represent the quasi-isovolumetric plastic strain of the backfill following settlement and mobilization of the frictional stress along the vertical walls. The MC yield parameters were selected so the surface would be close to the MSDPu yield surface for comparative purposes of the stresses obtained from numerical simulations; the corresponding MC strength parameters are *φ;* = 31◦, *c* = 3.8 kPa and *T*<sup>0</sup> = 0.2 kPa. A non-associated flow rule was imposed with *φ<sup>d</sup>* = 0◦ (dilation angle). It is noted that the values of the Poisson's ratio *ν<sup>b</sup>* and internal friction angle *φ* are interrelated through *ν<sup>b</sup>* = (1 − sin*φ*)/(2 − sin*φ*) to ensure a consistent at-rest earth pressure coefficient *K*<sup>0</sup> [28,29,44–46]. This results in a value of Poisson's ratio *ν<sup>b</sup>* = (1 − sin*φ*) / (2 − sin*φ*) = (1 − sin31◦) / (2 − sin31◦) = 0.327 and a Jaky's [47] earth pressure coefficient at rest *K*<sup>0</sup> = 1 − sin*φ* = 1 − sin31◦ = 0.485; the same value of *K*<sup>0</sup> is obtained from the equation based on Poisson's ratio. The theoretical horizontal stress due to the overburden can then be calculated as *σxx* = *K*0*σzz*.

**Figure 14.** The numerical FLAC3D model of the three-dimensional vertical backfilled stope in an elastic rock mass.

Figure 15 shows the vertical (*σzz*) and horizontal (*σxx*) stress distributions along the centerline (CL) of the backfilled stope obtained with the two models. It is seen that the general tendencies of the stress distributions obtained by numerical modeling with the MSDPu-EPP model and MC-EPP model are similar. The vertical and horizontal stresses follow the overburden stresses (marked as dashed lines in Figure 15) at shallow depth, but the overburden stresses significantly exceed the stresses in the backfill deeper in the opening due to the well-known arching effect associated with the frictional stress transfer along the rock walls.

**Figure 15.** Stress distributions along the centerline (CL) of the backfilled stope with the MSDPu-EPP model and the MC-EPP model obtained by FLAC3D for (**a**) vertical stress *σzz* and (**b**) horizontal stress *σxx*.

The trends obtained here are consistent with those given by previous numerical simulations and also by analytical solutions developed with the MC-EPP model [27,29,45,48,49]. The results shown here nonetheless indicate that the vertical and horizontal stresses obtained by the MSDPu-EPP model are smaller than those obtained by the MC-EPP model, due to the differences in the MSDPu and MC yield surfaces. Representative field experimental data are necessary to evaluate which one is more representative of the actual stress state.

One of the specific features of the MSDPu-EPP model is the introduction of a cap on the yield surface, which starts at the key parameter *Ic* (not included in other criteria and models). As indicated above, the value of *Ic* controls the cap position on the yield surface, which then departs from the (quasi) linear strength increase given by the most other models. Additional simulations were conducted to investigate the effect of *Ic* on the stress distributions in the backfilled stope. Three values (*Ic* = 20 kPa, 30 kPa, 100 kPa) were used with the MSDPu-EPP model, while the other parameters were kept constant for all cases (*Eb* = 300 MPa, *ν<sup>b</sup>* = 0.327, *ρ<sup>b</sup>* = 1800 kg/m3, *φ* = 30◦, *C*<sup>0</sup> = 10 kPa, *T*<sup>0</sup> = 0.2 kPa, *b* = 0.75, *a*<sup>3</sup> = 0.06 and *ξ* = 0.01).

Figure 16 shows the vertical (*σzz*) and horizontal (*σxx*) stress distributions along the CL obtained for different *Ic*. The results indicate that the vertical and horizontal stresses along the CL decrease with an increasing *Ic*, due to the higher shear strength (larger yield surface) of the backfill. The results also indicate that the arching effect tends to become more significant for stronger cemented (cohesive) backfill. These calculations highlight the importance of including the effect of the cap through *Ic*, which is a unique characteristic of the MSDPu-EPP model.

**Figure 16.** Stress distributions along the CL of the backfilled stope obtained from the numerical simulation conducted with the MSDPu-EPP model by considering different *Ic* for (**a**) vertical stress *σzz* and (**b**) horizontal stress *σxx*.

To further illustrate the features of the FLAC3D MSDPu-EPP model with associated and non-associated flow rules, the stresses in the backfilled stope shown in Figures 9b and 14 are analyzed by considering *ξ* = 1 for associated flow rule.

Figure 17 shows the vertical (*σzz*) and horizontal (*σxx*) stress distributions along the CL of the backfilled stope obtained by numerical modeling with the FLAC3D MSDPu-EPP model by considering *ξ* = 1 (associated flow) and *ξ* = 0.01 (non-associated flow), respectively. It can be seen that the vertical stress along the CL is much lower for the associated flow rule than that for the non-associated flow rule (Figure 17a), due to the significant volumetric strains (dilatancy) which tend to increase the arching effect. For the same reason, the horizontal stresses along the CL for the associated flow rule significantly exceed the overburden horizontal stress at shallow depth between around 0 to 5 m (Figure 17b). Similar results have been shown by Li and Aubertin [27] with the MC-EPP model for 2D backfilled stopes. Figure 18 shows the vertical stress distribution along the CL of the 3D backfilled stope (Figures 9b and 14) obtained by numerical modeling with the MC-EPP model by considering a non-associated (dilation angle *φ<sup>d</sup>* = 0◦) and associated (*φ<sup>d</sup>* = *φ* = 31◦) flow rules. The results suggest once again that the numerical modeling with an associated flow rule tends to lead an underestimation of the stresses in backfilled stopes.

**Figure 17.** Stress distributions along the CL of the backfilled stope obtained from the numerical simulation conducted with the MSDPu-EPP model by considering the non-associated (*ξ* = 0.01) and associated (*ξ* = 1) flow rules for (**a**) vertical stress *σzz* and (**b**) horizontal stress *σxx*.

**Figure 18.** Vertical stress *σzz* distributions along the CL of the backfilled stope obtained from numerical simulation conducted with the MC-EPP model by considering the non-associated (dilation angle *φ<sup>d</sup>* = 0◦) and associated (*φ<sup>d</sup>* = *φ* = 31◦) flow rules; other material parameters are *φ* = 31◦, *c* = 3.8 kPa and *T*<sup>0</sup> = 0.2 kPa.

#### **6. Discussion**

The non-associated MSDPu-EPP model recently implemented in FLAC3D is applied here for the analysis of underground backfilled stopes, to illustrate some of the features of this three-dimensional nonlinear model. The MSDPu criterion used here was previously shown to capture the essential characteristics of yielding and failure of a large variety of geomaterials, including soils, rocks and rockfills. The MSDPu-EPP model uses a closed yield surface with a controllable cap at high mean stress. It employs a non-associated flow rule to better describe the deviatoric and volumetric plastic strains. The availability of the FLAC3D MSDPu-EPP model will favor its application to help solve various geotechnical engineering problems such as those related to the behavior of slopes, tunnels, dams and barricades in underground mines.

Despite the various advantages of the model, there are a few limitations with the FLAC3D MSDPu-EPP model. For instance, the model has not yet been applied to cohesionless (granular) materials such as sand, rockfill, waste rock and uncemented backfill. This can be done by expressing the MSDPu criterion with *C*<sup>0</sup> = *T*<sup>0</sup> = 0 in Equations (10) and (11), which leads to *a*<sup>1</sup> = *a*<sup>2</sup> = 0. Equation (7) then becomes:

$$F\_0 = \left[a^2 I\_1^2 - a\_3 \langle I\_1 - I\_c \rangle^2\right]^{1/2} \tag{32}$$

This specific version of the criterion has not yet been validated in the context of the FLAC3D MSDPu-EPP model.

Another limitation of the model, common to all perfectly plastic models, is the absence of strain hardening and softening in the formulation. This is acceptable for many types of application, but not for some specific ones [3,4,8]. Additional work is considered for an evolving yield surface in the principal stress space with the MSDPu-EPP model.

Also, the model has not yet been validated for coupled problems involving pore water pressures (particularly for transient conditions).

Despite such limitations, the FLAC3D MSDPu-EPP model offers a powerful alternative to simulate the complex mechanical response of geomaterials as a built-in model in FLAC3D (and other commercial codes).

#### **7. Conclusions**

The non-associated MSDPu elastic-perfectly plastic constitutive model (MSDPu-EPP model) has been implemented in FLAC3D with a user-written DLL under the assistance of the C++ plug-in option. The resulting FLAC3D MSDPu-EPP model is validated, in part, by comparing simulation results with existing analytical solutions developed for evaluating the stress and displacement distributions around cylindrical openings. The FLAC3D MSDPu-EPP model is then used to analyze a specific three-dimensional geotechnical problem, taking into account the nonlinear stress-strain behavior, yielding and volumetric strains associated with the cap surface and considering an associated or non-associated flow rule. Application of the FLAC3D MSDPu-EPP model is illustrated with simulations conducted to analyze the stability of a 3D open stope and stress distribution in the backfilled stope. The results were compared with those obtained with the commonly used Mohr-Coulomb (MC-EPP) model. The numerical results obtained with the MSDPu-EPP model indicate that the size of yielded area around the empty stope significantly increases as the intermediate in-situ stress increases, indicating that the intermediate in-situ stress plays an important role in the response and stability of openings. This effect is not as clearly perceived with the MC-EPP model, however, as the simulated yielded area stays almost unchanged when the intermediate in-situ stresses are increased, hence showing that this model largely neglects this aspect. Other particular features of the FLAC3D MSDPu-EPP model related to the closed yield envelope (with a cap) are shown with the stress analysis of backfilled stopes. These results show that when a higher mean stress is necessary to reach backfill volumetric yield surface, defined by the closed (cap) envelope, the backfill becomes stronger, leading to more pronounced arching effect and smaller stresses in the backfilled stope. This is a unique characteristic obtained with numerical simulations conducted with the FLAC3D MSDPu-EPP model. It can thus be considered as a very useful numerical tool to analyze and solve geotechnical engineering problems. Nevertheless, improvements are still being considered for the model to minimize its limitations. For example, the model has not yet been applied to cohesionless (granular) materials; strain hardening and softening are not included in the formulation; the model has not been validated for couple problems involving pore water pressures. Work is underway to address these specific issues.

**Author Contributions:** F.Z.: formulation, programming, numerical modeling, literature and writing of the original draft. L.L.: project administration, supervision and editing of the original draft. M.A.: co-supervision and editing of the original draft. R.S.: co-supervision and editing of the original draft. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was financially supported by the Natural Sciences and Engineering Research Council of Canada (RGPIN89749-04, RGPIN-2018-06902, ALLRP-566888-21), Fonds de recherche du Québec-Nature et Technologies (2017-MI-202860), Mitacs (IT12569) and industrial partners of the Research Institute on Mines and the Environment (RIME UQAT-Polytechnique; http://rime-irme.ca/ accessed on 1 May 2022). The authors are grateful for their support.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Data are contained within the article.

**Acknowledgments:** The authors acknowledge the financial support from the Natural Sciences and Engineering Research Council of Canada, Fonds de recherche du Québec-Nature et Technologies, Mitacs and industrial partners of the Research Institute on Mines and the Environment.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A. Expressions for** *AS***,** *BS* **and** *C<sup>S</sup>*

When *I*<sup>1</sup> ≤ *Ic*, the expressions for *A*, *B* and *C* in Equation (25) are suffixed with a "*s*" to denote a shear response and expressed as follows:

$$\begin{split} A^{s} &= \left(\beta\_{1} - \beta\_{2}\right)^{2} \left[ \left(\frac{\partial Q^{s}}{\partial \sigma\_{1}} - \frac{\partial Q^{s}}{\partial \sigma\_{2}}\right)^{2} + \left(\frac{\partial Q^{s}}{\partial \sigma\_{2}} - \frac{\partial Q^{s}}{\partial \sigma\_{3}}\right)^{2} + \left(\frac{\partial Q^{s}}{\partial \sigma\_{3}} - \frac{\partial Q^{s}}{\partial \sigma\_{1}}\right)^{2} \right] \\ &- 6a^{2} \left(\beta\_{1} + 2\beta\_{2}\right)^{2} \left(F\_{\pi}^{I}\right)^{2} \left(\frac{\partial Q^{s}}{\partial \sigma\_{1}} + \frac{\partial Q^{s}}{\partial \sigma\_{2}} + \frac{\partial Q^{s}}{\partial \sigma\_{3}}\right)^{2} \end{split} \tag{A1}$$

$$\begin{cases} B^s = 12a^2 \left(\beta\_1 + 2\beta\_2\right) \left(I\_1^I - 2a\_1\right) \left(F\_{\pi}^I\right)^2 \left(\frac{\partial Q^s}{\partial \sigma\_1} + \frac{\partial Q^r}{\partial \sigma\_2} + \frac{\partial Q^s}{\partial \sigma\_3}\right) \\ - 2\left(\beta\_1 - \beta\_2\right) \left[\left(2\sigma\_1^I - \sigma\_2^I - \sigma\_3^I\right)\frac{\partial Q^r}{\partial \sigma\_1} + \left(2\sigma\_2^I - \sigma\_3^I - \sigma\_1^I\right)\frac{\partial Q^r}{\partial \sigma\_2} + \left(2\sigma\_3^I - \sigma\_1^I - \sigma\_2^I\right)\frac{\partial Q^s}{\partial \sigma\_3}\right] \end{cases} \tag{A2}$$

$$\mathbf{C}^{s} = J\_{2}^{I} - \left(F\_{0}^{\ast I}\right)^{2} \left(F\_{\pi}^{I}\right)^{2} \tag{A3}$$

where

$$\begin{array}{l} \frac{\partial Q^{\circ}}{\partial \sigma\_{1}} = \frac{1}{3} (2\sigma\_{1} - \sigma\_{2} - \sigma\_{3}) - 2\xi^{2}F\_{\pi}^{2} \left[ a^{2} (I\_{1} - a\_{1}) + \frac{3}{4} \left( \frac{1}{b^{2}} - 1 \right) \left( F\_{0}^{\circ} \right)^{2} F\_{\pi}^{2} \cos 3\theta \frac{\partial \theta}{\partial \sigma\_{1}} \right] \\\ \frac{\partial Q^{\circ}}{\partial \sigma\_{2}} = \frac{1}{3} (2\sigma\_{2} - \sigma\_{1} - \sigma\_{3}) - 2\xi^{2}F\_{\pi}^{2} \left[ a^{2} (I\_{1} - a\_{1}) + \frac{3}{4} \left( \frac{1}{b^{2}} - 1 \right) \left( F\_{0}^{\circ} \right)^{2} F\_{\pi}^{2} \cos 3\theta \frac{\partial \theta}{\partial \sigma\_{2}} \right] \\\ \frac{\partial Q^{\circ}}{\partial \sigma\_{3}} = \frac{1}{3} (2\sigma\_{3} - \sigma\_{1} - \sigma\_{2}) - 2\xi^{2}F\_{\pi}^{2} \left[ a^{2} (I\_{1} - a\_{1}) + \frac{3}{4} \left( \frac{1}{b^{2}} - 1 \right) \left( F\_{0}^{\circ} \right)^{2} F\_{\pi}^{2} \cos 3\theta \frac{\partial \theta}{\partial \sigma\_{3}} \right] \end{array} \tag{A4}$$

and

$$\begin{split} \frac{\partial\theta}{\partial\sigma\_{1}} &= \frac{\sqrt{3}(\sigma\_{1}-\sigma\_{2})(\sigma\_{1}-\sigma\_{3})(\sigma\_{2}-\sigma\_{3})^{2}}{6l\_{2}\sqrt{4l\_{2}^{3}-27l\_{3}^{2}}}\\ \frac{\partial\theta}{\partial\sigma\_{2}} &= \frac{\sqrt{3}(\sigma\_{2}-\sigma\_{3})(\sigma\_{2}-\sigma\_{1})(\sigma\_{1}-\sigma\_{3})^{2}}{6l\_{2}\sqrt{4l\_{2}^{3}-27l\_{3}^{2}}}\\ \frac{\partial\theta}{\partial\sigma\_{3}} &= \frac{\sqrt{3}(\sigma\_{3}-\sigma\_{1})(\sigma\_{3}-\sigma\_{2})(\sigma\_{1}-\sigma\_{2})^{2}}{6l\_{2}\sqrt{4l\_{2}^{3}-27l\_{3}^{2}}} \end{split} \tag{A.5}$$

#### **Appendix B. Expressions for** *Av, Bv* **and** *Cv*

When *I*<sup>1</sup> > *Ic*, the expressions for *A*, *B* and *C* in Equation (25) are suffixed with a "*v*" to denote volumetric response and expressed as follows:

$$\begin{split} A^{\overline{v}} &= \left(\beta\_{1} - \beta\_{2}\right)^{2} \Bigg[ \left(\frac{\partial Q^{v}}{\partial \sigma\_{1}} - \frac{\partial Q^{v}}{\partial \sigma\_{2}}\right)^{2} + \left(\frac{\partial Q^{v}}{\partial \sigma\_{2}} - \frac{\partial Q^{v}}{\partial \sigma\_{3}}\right)^{2} + \left(\frac{\partial Q^{v}}{\partial \sigma\_{3}} - \frac{\partial Q^{v}}{\partial \sigma\_{1}}\right)^{2} \Bigg] \\ &- 6\left(\alpha^{2} - a\_{3}\right) \left(\beta\_{1} + 2\beta\_{2}\right)^{2} \left(F^{I}\_{\pi}\right)^{2} \left(\frac{\partial Q^{v}}{\partial \sigma\_{1}} + \frac{\partial Q^{v}}{\partial \sigma\_{2}} + \frac{\partial Q^{v}}{\partial \sigma\_{3}}\right)^{2} \end{split} \tag{A6}$$

$$\begin{split} B^{\boldsymbol{\upsilon}} &= 12(\boldsymbol{\beta}\_{1} + 2\boldsymbol{\beta}\_{2}) \left[ a^{2} \left( \boldsymbol{I}\_{1}^{\boldsymbol{I}} - \boldsymbol{a}\_{1} \right) - a\_{3} \left( \boldsymbol{I}\_{1}^{\boldsymbol{I}} - \boldsymbol{I}\_{\boldsymbol{\varepsilon}} \right) \right] \left( \boldsymbol{F}\_{\pi}^{\boldsymbol{I}} \right)^{2} \left( \frac{\partial Q^{\boldsymbol{\upsilon}}}{\partial \boldsymbol{\sigma}\_{1}} + \frac{\partial Q^{\boldsymbol{\upsilon}}}{\partial \boldsymbol{\sigma}\_{2}} + \frac{\partial Q^{\boldsymbol{\upsilon}}}{\partial \boldsymbol{\sigma}\_{3}} \right) \\ &- 2(\beta\_{1} - \beta\_{2}) \left[ \left( 2\boldsymbol{\sigma}\_{1}^{\boldsymbol{I}} - \boldsymbol{\sigma}\_{2}^{\boldsymbol{I}} - \boldsymbol{\sigma}\_{3}^{\boldsymbol{I}} \right) \frac{\partial Q^{\boldsymbol{\upsilon}}}{\partial \boldsymbol{\sigma}\_{1}} + \left( 2\boldsymbol{\sigma}\_{2}^{\boldsymbol{I}} - \boldsymbol{\sigma}\_{3}^{\boldsymbol{I}} - \boldsymbol{\sigma}\_{1}^{\boldsymbol{I}} \right) \frac{\partial Q^{\boldsymbol{\upsilon}}}{\partial \boldsymbol{\sigma}\_{2}} + \left( 2\boldsymbol{\sigma}\_{3}^{\boldsymbol{I}} - \boldsymbol{\sigma}\_{1}^{\boldsymbol{I}} - \boldsymbol{\sigma}\_{2}^{\boldsymbol{I}} \right) \frac{\partial Q^{\boldsymbol{\upsilon}}}{\partial \boldsymbol{\sigma}\_{3}} \right] \end{split} \tag{A7}$$

$$\mathbf{C}^{v} = J\_2^I - \left(F\_0^{vI}\right)^2 \left(F\_\pi^I\right)^2 \tag{A8}$$

where

$$\begin{cases} \frac{\partial Q^{\mathcal{T}}}{\partial V\_{1}} = \frac{1}{3} (2\sigma\_{1} - \sigma\_{2} - \sigma\_{3}) - 2\frac{\nu}{5} F\_{\mathcal{R}}^{2} \left[ a^{2} (I\_{1} - a\_{1}) - a\_{3} (I\_{1} - I\_{\epsilon}) + \frac{3}{4} \left( \frac{1}{b^{2}} - 1 \right) \left( F\_{0}^{\mathcal{V}} \right)^{2} F\_{\mathcal{R}}^{2} \cos 3\theta \frac{\partial \theta}{\partial V\_{1}} \right] \\ \frac{\partial Q^{\mathcal{V}}}{\partial \nu\_{2}} = \frac{1}{3} (2\sigma\_{2} - \sigma\_{1} - \sigma\_{3}) - 2\frac{\nu}{5} F\_{\mathcal{R}}^{2} \left[ a^{2} (I\_{1} - a\_{1}) - a\_{3} (I\_{1} - I\_{\epsilon}) + \frac{3}{4} \left( \frac{1}{b^{2}} - 1 \right) \left( F\_{0}^{\mathcal{V}} \right)^{2} F\_{\mathcal{R}}^{2} \cos 3\theta \frac{\partial \theta}{\partial v\_{2}} \right] \\ \frac{\partial Q^{\mathcal{V}}}{\partial \nu\_{3}} = \frac{1}{3} (2\sigma\_{3} - \sigma\_{1} - \sigma\_{2}) - 2 \frac{\nu}{5} F\_{\mathcal{R}}^{2} \left[ a^{2} (I\_{1} - a\_{1}) - a\_{3} (I\_{1} - I\_{\epsilon}) + \frac{3}{4} \left( \frac{1}{b^{2}} - 1 \right) \left( F\_{0}^{\mathcal{V}} \right)^{2} F\_{\mathcal{R}}^{2} \cos 3\theta \frac{\partial \theta}{\partial v\_{3}} \right] \end{cases} \tag{A}$$

#### **Notations:**

The following symbols are used in the paper:



#### **References**

