Slope stability analysis: Barodesy vs linear elastic – perfectly plastic models

The results of slope stability analysis are not unique. Different factors of safety are obtained investigating the same slope. The differences result from different constitutive models including different failure surfaces. In this contribution, different strength reduction techniques for two different constitutive models (linear elastic perfectly plastic model using a Mohr-Coulomb failure criterion and barodesy) have been investigated on slope stability calculations for two different slope inclinations. The parameters for Mohr – Coulomb are calibrated on peak states of element tests simulated with barodesy for different void ratios. For both slopes the predictions of the factors of safety are higher with barodesy than with Mohr-Coulomb. The difference is to some extend explained by the different shapes of failure surfaces and thus different values for peak strength under plane strain conditions. The plane strain predictions of Mohr-Coulomb are conservative compared to barodesy, where the failure surface coincides with Matsuoka-Nakai.


Introduction
Strength reduction techniques are a widely used method to calculate factors of safety in Finite Element calculations. In this method, the strength parameters are reduced gradually until failure occurs. Mostly constitutive models including the Mohr -Coulomb failure criterion are used for strength reduction finite element analysis (SRFEA). The Mohr -Coulomb model (linear elasticperfectly plastic) does not take into account effects of changes in density and consequently effects like softening cannot be modelled.
In this work, results of slope stability analyses of conventional elasto -plastic models using a MC failure criterion are compared with results of strength reduction applied to barodesy. Barodesy is a material model, which has similarities to hypoplasticity. The model includes a failure criterion, which practically coincides with the Matsuoka -Nakai criterion [9]. This article investigates on the one hand the influence of the failure criterion and on the other hand the performance of the strength reduction procedures. Additionally, finite element limit analyses (FELA) are carried out. Since FELA provides rigorous solutions these results are used as reference solutions for comparisons.

Barodesy
Barodesy, introduced by Kolymbas [1], is a material model for granular materials which has similarities to hypoplasticity: As in hypoplasticity, the stress rate is formulated as a function of stress, stretching and void ratio: ̇= ( ,̇, ). Compared with elasto-plastic models, barodesy does not distinguish between elastic and plastic strain. Barodesy [2,3] includes concepts from Critical State Soil Mechanics: In order to define normally consolidated states, it includes an isotropic normal compression line (NCL) according to Butterfield [4]. Barodesy further includes a stress-dependent critical void ratio [5] (CSL), in order to distinguish between highly overconsolidated ( < ) and normal to slightly overconsolidated ( > ) clay.   1 shows simulations of drained triaxial tests with barodesy with a constant initial void ratio ini = constant and different initial mean stresses ini . Peak states of highly overconsolidated clay reach maximum mobilized friction angles that are larger than the critical friction angle, see Fig. 1 (b). It is thus clear that the peak strength envelope is dependent on the initial void ratio and initial mean stress. For a constant initial void ratio follows: The lower the initial mean stress is, the higher is the peak friction angle, e.g. [6,7]. In barodesy critical stress states almost coincide with predictions by Matsuoka-Nakai [9] and results from laboratory experiments [10]. In Fig. 2, the critical stress surface is shown for a critical friction angle of 26°. The failure surface according to Mohr-Coulomb is added for ´= 26° and ´= 0 and for ´= 22° and ´= 6 kPa. The discrepancies between the constitutive models, especially for plane strain conditions are apparent. The range of = 0° to = 15° marks plane strain states [10].

Calibration
In order to compare the peak strength envelope in barodesy with a material model, which includes the Mohr-Coulomb failure criterion, the peak envelope has to be linearized for the relevant/expected stress range. If the soil is highly overconsolidated, i.e. if the overconsolidation ratio (OCR) is larger than 2, strength is defined by a peak friction angle which is larger than the critical friction angle. For normally to slightly overconsolidated soil, the maximum mobilized friction angle is the critical friction angle. For the approximation of the peaks of the drained triaxial tests with the initial void ratio of ini = 0.65, it follows ´= 10 kPa and ´= 30°. These parameters are used for geometry 2. Table 1 summarizes all parameters used in the following studies. Table 1. Parameters and initial void ratio ini for barodesy and parameters for the Mohr-Coulomb model Note that barodesy has the same set of parameters for different densities. In the constitutive model the density is considered by the void ratio or overconsolidation ratio. For geometry 1 the initial void ratio is ini = 0.83, for geometry 2 the initial void ratio is ini = 0.65. For low stress levels the linearization with ´ and ´ predicts higher strength as the peak envelope of barodesy under axisymmetric compression. The Figures 1 and 3 only display axisymmetric conditions.

Calibration Results
Axisymmetric compression corresponds to a Lode angle of 30°, axisymmetric extension is defined by = −30°. In many geotechnical applications, plane strain conditions are more relevant than axisymmetric conditions. It is therefore interesting to investigate stress paths under plane strain compression, cf. Fig. 2. Stress paths under plane strain conditions are reported to occur within 0° < < 15° [10]. In [2,11] investigations under plane strain compression have been carried out with barodesy, which show that the resulting plane strain stress paths are within this range. For the Mohr-Coulomb model, the direction of the elastic stress paths depends also on , e.g. for = 0.3, it follows = 13°, cf. [12]. Note that in the range of 0° < < 15°, barodesy predicts higher maximum mobilized friction angles than Mohr-Coulomb due to the different failure surfaces, see Fig. 2 and Fig. 4. In Fig. 4 simulations of a conventional drained (CD) triaxial test and a plane strain biaxial test are shown. The initial stress states are chosen so that failure occurs at ′ = 72 kPa. Therefore, the models' predictions of the maximum mobilized deviatoric stress for axisymmetric compression coincide. However, under plane strain conditions, the maximum mobilized deviatoric stress is larger in the simulation with barodesy than in the simulations with the Mohr-Coulomb model.  Slope stability calculations of two geometries have been performed to evaluate the influence of both, the different failure criteria and strength reduction methods. One slope with an inclination of approx. 26° (Fig. 5) and one steep slope with an inclination of approx. 45° (Fig. 6).

Barodesy
The strength reduction finite element analyses have been performed using Plaxis 2D Version 2017 [13]. Identical meshes have been used for all SRFEA calculations. The final models (as shown in the figures 5 and 6) consist of 15 noded elements using a shape function of 4 th order. The FELA have been performed using Optum G2 [14]. The parameters were chosen according to the calibrations mentioned in Section 4, Tab. 1. For the barodetic material model, the parameters of the material were equal for both geometries, but different void ratios have been considered in the calculations. For the analyses performed with the linear elastic -perfectly plastic model using the failure criterion of Mohr -Coulomb two different parameter sets for geometry 1 and geometry 2 (to take into account the different peak states depending on the density) have been used. Additionally, one parameter set (Set 2) has been used for the calculations with the critical friction angle and a cohesion of 0 kPa, cf. Fig.4. The calibration details are explained in Section 3. Mohr -Coulomb failure criterion using an associated flow rule in combination with reduced strength parameters according to Davis [15]. This method is referred to as Davis A in the following. Finite Element Limit Analyses (FELA) with Mohr -Coulomb failure criterion using an associated flow rule in combination with reduced strength parameters according to Davis B [16]. Previous studies showed that Davis B yields more realistic results compared to Davis  on the effective strength parameters applied in the current iteration, thus the decrease if non-associativity is taken into account [16]. During the strength reduction procedure, the typical shear bands become visible in the slope cf. Fig. 7. For the displacement based Finite Element calculation with Mohr -Coulomb also the plastic points and the tension cut-off point can be plotted as calculation result, cf. Fig. 8.   Fig. 7. Typical shear band (incremental deviatoric strain)

Fig. 8. Plastic points (red dots) and tension cut off points (white squares)
The effect of tensile stresses on slope stability analyses has been investigated in [20]. Barodesy does not need an additional tension cut off formulation as the constitutive model excludes tensile stresses. Further research on the effect of tensile stresses on computed factors of safety is of interest as tensile stresses could affect the slope stability significantly, especially for steep slopes in combination with materials with high values of ′, since the critical sliding surface is short and the tension zone is relatively large (see Fig. 8).

Results of the Strength Reduction Analyses
The strength reduction method for barodesy has been introduced in [17] and barodesy has proofed to be well suited to simulate realistic shear bands [18,19]. In the SRFEA with barodesy in this article, the critical friction angle is gradually reduced until failure. Thus, the procedure is comparable to strength reduction applied to MC models. Due to the associated flow rule, strength reduction calculations using an associated flow rule compute higher factors of safety than SREFA using a non-associated flow rule (e.g. 1.28 vs. 1.24 for the flatter slope, cf. Table 2). As shown in previous studies, FELA with an associated flow rule are in very good agreement with SRFEA (assuming ′ = ′). These values are un-conservative and questionable, as soil does not follow associated flow. The results of the displacement based Finite Element calculations using a non-associated flow rule are in relatively good agreement with the Finite Element Limit analysis considering the non-associated behaviour by reduced parameters according to Davis [15]. The calculations with the method Davis B [16] perform even better (1.24 vs. 1.23 for the flatter slope). Those results apply for both slope inclinations.
By reducing strength in boundary value problems, plasticity and consequently displacements occur also before reaching failure. This results in load redistribution and may influence the factors of safety.
The strength reduction method of the barodetic model shows plausible results for the slope with the steeper inclination of 45°. The factor of safety is slightly higher than the factor of safety of the calculations with Mohr - Coulomb (1.18 vs. 1.10). This can be explained with the different shapes of the failure surfaces in the deviatoric plane, cf. Fig. 4. The parameters are fitted and coincide only for axisymmetric triaxial compression. In the range of interest for plane strain conditions, the cone of barodesy lies outside the failure surface according to Mohr -Coulomb, cf. Fig. 4. It has to be pointed out that experimental results confirm the shape of the barodetic cone [10].
For the slope with the lower inclination of 26°, the trend is the same, though more pronounced. For this slope ′ is 22°/26°, for the steeper slope ′ is 30°. The lower the friction angle is the larger is the difference between Mohr-Coulomb and Matsuoka-Nakai under plane strain conditions [12].
Future investigations will reveal how this effect might be superimposed by the effect of load redistribution.

Conclusion
In experiments, normally consolidated to slightly overconsolidated samples show contracting behaviour without any softening. Highly overconsolidated samples show dilatant behaviour and develop a peak with subsequent softening. As the void ratio is included as state parameter in barodesy, different relative densities can be considered by using different void ratios. For the Mohr -Coulomb parameters ′ and ′, we consider peak states of barodesy in the calibration. This leads to two parameter sets for Mohr -Coulomb and one parameter set for barodesy with two different initial void ratios to consider the influence of two different densities.
For both slopes the predictions of the FoS are higher with barodesy than with MC. The difference is to some extend explained by the different shapes of failure surfaces and thus different values for peak strength under plane strain conditions. The plane strain predictions of Mohr-Coulomb are conservative compared to Matsuoka-Nakai. The lower the friction angle is, the higher is the deviation between the two failure surfaces. Future research is necessary to quantify the different individual influences (effect of cohesion, failure criterion, and load redistribution) separately.