Optimization and variational principles for the shear strength reduction method
Stanislav Sysala, Eva Hrubešová, Zden?k Michalec, Franz Tschuchnigg
OOptimization and variational principles for the shearstrength reduction method
Stanislav Sysala ∗ , Eva Hrubeˇsov´a , , Zdenˇek Michalec Institute of Geonics of the Czech Academy of Sciences, Ostrava, Czech Republic VˇSB - Technical University of Ostrava, Faculty of Civil Engineering, Ostrava, Czech Republic
January 5, 2021
Abstract
This paper is focused on a rigorous variant of the shear strength reductionmethod (SSR) and the corresponding determination of safety factors. The SSR-based safety factor is defined as a solution of an optimization problem that is in-dependent of the plastic flow rule and the space discretization. In case of nonas-sociative plasticity, a modified Davis approach is used. The optimization problemis analyzed and the corresponding duality between the static and kinematic princi-ples is derived. For numerical solution, a regularization method is introduced anda relation between the original and regularized problems is derived. The regular-ization method is combined with the finite element method, mesh adaptivity anda damped Newton method. In-house codes in Matlab are used for implementationof this solution concept. Two slope stability problems are considered, one of whichfollows from analysis of a real slope. Softwares Plaxis and Comsol Multiphysics areused for comparison of the results.
Keywords: slope stability, shear strength reduction method, convex optimization, staticand kinematic principles, regularization, finite elements and mesh adaptivity ∗ corresponding author, email: [email protected] a r X i v : . [ c s . C E ] J a n Introduction
This paper deals with slope stability assessment, which includes determination of thefactor of safety (FoS) and estimation of failure zones (slip surfaces) for a critical stateof the slope. Limit equilibrium (LE), shear strength reduction (SSR) or limit analysis(LA) belong among basic methods that can be used for the determination of FoS. Themethods arise from elastic-perfectly plastic models containing mainly the Mohr-Coulombyield criterion.The LE method is based on predefined failure zones, see, for example, [Duncan (1996)],[Yu (2006)]. It does not require finite element computation and stresses equilibrium atevery point in a domain arround the slope. Due to these facts, this method is simple andwidely used in geotechnical practice. On the other hand, accuracy of the solution cannotbe easily verified, especially if anisotropic or inhomogeneous materials are considered orif a complex geometry is defined.The SSR method [Zienkiewicz et al. (1975), Brinkgreve and Bakker (1991)],[Dawson et al. (1999), Griffiths and Lane (1999)] has been suggested mainly for the slopestability assessment. It is a conventional method based on a displacement variant of thefinite element method (FEM) and on reduction of strength parameters defining the Mohr-Coulomb model. It has also been implemented within some commercial softwares likePlaxis [Brinkgreve (2011)] or Middas GTS NX. Strong dependence on finite element meshdensity or even nonunique determination of FoS can happen in some cases with the non-associative plastic flow rule, see [Tschuchnigg et al. (2015a), Tschuchnigg et al. (2015b)].LA is a universal method that can be used for various stability problems, not only for theones in geotechnical practice. FoS is derived from a critical (limit) value of the load factor.Originally, this method was purely analytical, see, for example, [Chen and Liu (1990)].Now, it is rather a numerical method based on optimization, duality between kine-matic and static principles and on FEM [Christiansen (1990), Sloan (2013), Yu (2006),Haslinger et al. (2019)]. LA is supported by mathematical and numerical analyses, see[Temam (1985), Christiansen (1990), Haslinger et al. (2016a), Repin et al. (2018)],[Haslinger et al. (2019)]. On the other hand, this method is not so conventional in slopestability and may predict different FoS than the SSR or LE methods. The theory of LAmethod is not also fully developed for the non-associative plastic flow.The usage of the non-associative plastic flow rule is supported by laboratory experimentsin order to catch volume dilatancy of materials with internal friction. On the other hand,mathematical theory of non-associative elastic-plastic problems is missing or at least in-2omplete. Especially, standard implicit discretizations of time (pseudo-time) variable leadto problematic numerical behavior. The drawbacks of the non-associative models can besuppressed by variational approaches based on theory of bipotentials [Hjiaj et al. (2003),Hamlaoui et al. (2017)] or semi-implicit time schemes [Krabbenhoft (2012)]. Within theLA method, Davis [Davis (1968)] suggested to modify the strength parameters and con-sequently approximate the nonassociative model by the associative one. In[Tschuchnigg et al. (2015a), Tschuchnigg et al. (2015b), Tschuchnigg et al. (2015c)], theDavis approach has been modified for purposes of the SSR method.Inspired by the recent papers from Tschuchnigg et al. mentioned above, we suggest anoptimization definition of the SSR problem that is independent of the plastic flow ruleand a space discretization. Basic properties of the optimization problem are shown andduality between the static and kinematic settings of this problem is derived. Within thenonassociative plasticity, two different optimization approaches are suggested.For numerical solution, we propose a regularization of the optimization problems inspiredby recent papers [Sysala et al. (2015), Cermak et al. (2015), Haslinger et al. (2016a)],[Haslinger et al. (2016b), Repin et al. (2018), Haslinger et al. (2019), Reddy and Sysala (2020)].Convergence with respect to the regularization parameter is analyzed in order to relatethe original and regularized problems. The regularization enables to solve the problemby the standard finite element methods and by the damped Newton method suggestedin [Sysala (2012)]. Mesh adaptivity is also used to achieve more accurate results. Theproblem is implemented within in-house Matlab codes. These codes based on elastic-plastic solvers have been systematically developed and described in [Sysala et al. (2016),Sysala et al. (2017), Cermak et al. (2019)]. Some of the codes are available for down-load [Cermak et al. (2018)]. For comparison of the results with standard approaches, thecommercial softwares Plaxis and Comsol Multiphysics are used.The rest of the paper is organized as follows. In Section 2, we introduce a scheme ofthe elastic-perfectly plastic problem with the Mohr-Coulomb yield criterion and presentmain ideas of the SSR and LA methods. In Section 3, the optimization definitions of thesemethods are introduced and compared for the associative plasticity. Section 4 is devoted tothe non-associative case where the Davis approach is used and modified for purposes of theSSR method. In Section 5, variational principles, duality and the kinematic approachesto the SSR method are presented. The regularization method is introduced and analyzedin Section 6. Numerical examples illustrating the efficiency of the suggested numericalmethods are presented in Section 7. Concluding remarks to the paper are introduced inSection 8. Appendix contains a closed form of the regularized dissipative funtion to the3ohr-Coulomb yield criterion.
We consider an elastic-perfectly plastic problem based on the Mohr-Coulomb yield crite-rion, a non-associative plastic flow rule and on monotone, proportional loading. We brieflyintroduce this problem and refer to [de Souza Neto et al. (2011), Sysala et al. (2017)] formore details.It is assumed that an investigated body occupies a bounded domain Ω ⊂ R with boundary ∂ Ω. This boundary is split into two disjoint parts, ∂ Ω u and ∂ Ω f , where the body is fixedand the surface force f is prescribed, respectively. Volume forces acting in Ω are denotedby F . Let x ∈ Ω and t ∈ (0 , t max ) denote the space and time (pseudo-time) variables,respectively. The assumption on the monotone and proportional loading means that thevolume and surface forces are in the forms F = t F in Ω , f = t f on ∂ Ω f , t ∈ (0 , t max ) , (2.1)where F : Ω → R and f : ∂ Ω f → R denote given forces independent of t . So thepseudo-time variable t is also the load factor.Deformation of the body is described by the displacement vector field u satisfying theboundary and initial conditions u = on ∂ Ω u × (0 , t max ) , u = in Ω × { } , (2.2)and by the infinitesimal strain ε where ε := ε ( u ) = ( ∇ u + ( ∇ u ) (cid:62) ) . (2.3)The strain ε is decomposed into elastic and plastic components ε e and ε p , respectively,according to ε = ε e + ε p . (2.4)The stress state of the body is represented by the symmetric Cauchy stress tensor σ − div σ = F , (2.5)and the Neumann boundary conditions σn = f on ∂ Ω f × (0 , T ) , (2.6)where n denotes the outward unit normal to the boundary ∂ Ω. The stress σ also satisfiesthe Mohr-Coulomb yield criterion( σ − σ ) + ( σ + σ ) sin φ − c cos φ ≤ , (2.7)where σ and σ denote the maximal and minimal principle stresses of σ , respectively, c is the cohesion and φ denotes the friction angle. (Notice that the standard mechanicalsign convection has been used.) Using the well-known formulascos φ = 1 (cid:112) φ , sin φ = tan φ (cid:112) φ , (2.8)we arive at the following equivalent form of the yield criterion:Φ( σ ) := ( σ − σ ) (cid:112) φ + ( σ + σ ) tan φ − c ≤ , (2.9)which is more appropriate for consequent analysis of the shear strength reduction (SSR)method.The elastic-perfectly plastic model is completed by the initial-value constitutive model: σ = C ( ε − ε p ) in Ω × (0 , t max ) , ˙ ε p ∈ ˙ γ∂ Ψ( σ ) in Ω × (0 , t max ) , ˙ γ ≥ , Φ( σ ) ≤ , ˙ γ Φ( σ ) = 0 , in Ω × (0 , t max ) ε p = in Ω × { } . (2.10)Here, Ψ( σ ) := ( σ − σ ) (cid:112) ψ + ( σ + σ ) tan ψ − c, (2.11)is the plastic potential, ψ if the dilatation angle, γ denotes the plastic multiplier, the dotsymbol means the pseudo-time derivative of a quantity and ∂ Ψ( σ ) is the subdifferential5f Ψ at σ . The fourth order tensor C represents the linear isotropic elastic law: σ = C ε e = 13 (3 K − G )( I : ε e ) I + 2 G ε e , (2.12)where I is the unit 2nd-order tensor, “:” denotes the biscalar product and K, G >
K > G , denote the bulk and shear moduli, respectively.The elastic-perfectly plastic problem is thus defined by the formulas (2.1)–(2.12). Thegiven material parameters φ and ψ satisfy 0 ≤ ψ ≤ φ ≤ π/
2. In the case ψ = φ , we arriveat the associative model which is supported by mathematical theory based on the principleof maximum plastic dissipation, see, for example, [Han and Reddy (2012)]. Numericalmethods for the associative plasticity are also more stable and simpler. For example, it isnot necessary to use the arc-length control of the loading unlike the nonassociative case.Now we sketch the main idea of the shear strength reduction (SSR) method and the limitanalysis (LA) method. More precise definitions are introduced and derived in the nextsections. The SSR method is based on the choice t max = 1 and on a strength reductionparameter λ >
0. This parameter is applied on the strength parameters c , φ and ψ asfollows: c λ := cλ , φ λ := arctan tan φλ , (2.13)and ψ λ := ψ until ψ < φ λ , then ψ λ := φ λ . (2.14)The corresponding factor of safety (FoS) is formally defined as a maximum of λ for whichthe elastic-perfectly plastic problem has a solution with respect to the parameters c λ , φ λ , ψ λ and t max = 1.The FoS for the LA method can be formally defined as the maximum of the load factor t (see (2.1)) for which the elastic-perfectly plastic problem has a solution with respect tothe parameters c , φ and ψ . This method is mostly defined for the associative problemwith φ = ψ . For the non-associative case, the strength parameters can be modified byDavis’ approach in order to stabilize the computation, see Section 4.FoS determined by the SSR and LA methods are different, in general. The SSR-FoSis usually closer to one than the LA-FoS. Therefore, if LA-FoS is close to one, one canexpect that both the safety factors are almost the same. In [Tschuchnigg et al. (2015a),Tschuchnigg et al. (2015b), Tschuchnigg et al. (2015c)], a parametrization of the LA methodhas been proposed to be the factors LA-FoS and SSR-FoS almost the same. The parametriza-tion will be described below in more details.6 Shear strength reduction for associative plasticityand its comparison with limit analysis
In this section, we introduce the static principle of the SSR method for the associativemodel with φ = ψ . In this case, it is well-known that the influence of the plastic flow ruleand the elasticity tensor C (see (2.10), (2.12)) is negligible on the factor of safety (FoS)for the SSR method and also for the LA method. In the following three subsections, weintroduce definitions of FoS for the SSR and LA methods, analyze and compare them. Let λ ∗ denote FoS for the SSR method. If we neglect the influence of the plastic flow ruleand the elasticity tensor on λ ∗ then we arrive at the following time-independent definition: λ ∗ = supremum of λ ≥ subject to − div σ = F in Ω , σn = f on ∂ Ω f , ( σ − σ ) (cid:113) tan φλ + ( σ + σ ) tan φλ ≤ cλ in Ω . (3.1)It means that we look for the maximum over λ for which there exists the stress σ satisfying(3.1). In particular, the first line in (3.1) is derived from (2.5), (2.6) and is satisfied for statically admissible stresses. The second line is derived from (2.9), (2.13) and saysthat σ is plastically admissible with respect to λ . According to the literature on convexanalysis [Ekeland and Temam (1974), Temam (1985), Christiansen (1990)], we rather usethe supremum than the maximum in this definition because for the critical value λ ∗ , theadmissible stress σ satisfying (3.1) need not exist on functional spaces. Although thedefinition admits the case λ ∗ = + ∞ , one can expect that λ ∗ is finite in geotechnicalpractice.Multiplying the inequality in the second line of (3.1) by λ , one can arrange the constraints(3.1) into the following form: λ ∗ = supremum of λ ≥ subject to − div σ = F in Ω , σn = f on ∂ Ω f , Φ( λ ; σ ) ≤ . (3.2)7here Φ( λ ; σ ) := ( σ − σ ) (cid:112) λ + tan φ + ( σ + σ ) tan φ − c. (3.3)The form (3.2) is more convenient for analysis of the SSR problem than (3.1) as weshall see below. By comparing (2.9) with (3.3) we have Φ(1; σ ) = Φ( σ ). The followingstatement is important to be the safety factor λ ∗ well defined. Lemma 3.1. If (3.2) is satisfied for some λ := λ > then (3.2) holds for any λ < λ .Proof. For any σ fixed, σ − σ ≥ λ (cid:55)→ Φ( λ ; σ ) is nondecreasing.Hence, if there exists σ such that (3.2) holds for some λ := λ > σ also satisfies(3.2) for any λ < λ .One can also expect in slope stability or similar problems that there exists σ satisfying(3.2) for λ = 0, that is, − div σ = F in Ω , σn = f on ∂ Ω f ,σ tan φ ≤ c in Ω . (3.4)It suffices to choose the stress field derived from the solution of elastic or elastic-plasticproblem. The corresponding maximal principle stress σ should be even nonpositive in Ω.Let us recall that the Mohr-Coulomb yield surface is the pyramid aligned with thehydrostatic axis, see, for example, [de Souza Neto et al. (2011)]. From the inequalityΦ( λ ; σ ) ≤ λ . By reducingthe strength parameters (i.e., by enlarging λ ) the slope of the Mohr-Coulomb pyramid isreduced. For λ →
0, the pyramid varies to a half-space.
Although this paper is focused on the SSR method, we also introduce the LA methodand compare them in order to generalize the Davis approach for the SSR method, see thenext section.Let ζ ∗ denote FoS for the LA method. This value is defined as follows:8 ∗ = supremum of ζ ≥ subject to − div σ = ζ F in Ω , σn = ζ f on ∂ Ω f , ( σ − σ ) (cid:112) φ + ( σ + σ ) tan φ ≤ c in Ω . (3.5)Let us note that the parameter ζ is applied on the external forces and ζ ∗ is also called thelimit load factor. The definition above is the static principle of LA which is also knownin literature as the lower bound theorem of LA because for any ζ < ζ ∗ there exists σ satisfying (3.5). Further, the second line in (3.5) can be written in the formΦ( σ ) = Φ(1; σ ) ≤ , (3.6)see (2.9) and (3.3). The LA problem is well-defined from the following reasons:1. σ = satisfies the constraints (3.5) for ζ = 0.2. If (3.5) is satisfied for some ζ := ζ > ζ < ζ . Indeed, if σ satisfies (3.5) for some ζ := ζ > ζ/ζ ) σ satisfies (3.5) for any ζ < ζ .The value ζ ∗ = + ∞ is expected for larger values of the friction angle, in particular for φ ≥ ◦ . It is also well-known that ζ ∗ depends linearly on the cohesion. Indeed, by thesubstitution σ := ζ σ , one can rewrite the LA problem as follows: ζ ∗ = supremum of ζ ≥ subject to − div σ = F in Ω , σn = f on ∂ Ω f , ( σ − σ ) (cid:112) φ + ( σ + σ ) tan φ ≤ cζ in Ω . (3.7)Hence, if ζ ∗ is the limit load factor for the cohesion c then αζ ∗ is the limit value corre-sponding to αc for any α >
0. From (3.7), we also see that the proportional enlarging ofthe external forces f and F is equivalent to the reduction of the cohesion. During this re-duction, the apex of the Mohr-Coulomb pyramid is moved to the origin of the coordinatesbut the shape of the pyramid remains the same unlike the SSR method. In order to compare the SRR and LA methods, we arise from (3.7) and parametrize theLA problem as follows: 9 iven λ > , find (cid:96) ( λ ) = supremum of ζ ≥ subject to − div σ = F in Ω , σn = f on ∂ Ω f , ( σ − σ ) (cid:112) λ + tan φ + ( σ + σ ) tan φ ≤ cζ in Ω . (3.8)We shall analyze properties of the function (cid:96) . Clearly, ζ ∗ = (cid:96) (1). Lemma 3.2.
The function (cid:96) is nonincreasing.Proof.
Let λ ≤ λ and ζ ≥ ζ < (cid:96) ( λ ). Then thereexists σ satisfying (3.8) for λ := λ . Since the function λ (cid:55)→ (cid:112) λ + tan φ is increasing, itis easy to see that σ satisfies (3.8) also for λ := λ . Therefore, ζ ≤ (cid:96) ( λ ) and consequently (cid:96) ( λ ) ≥ (cid:96) ( λ ).For reasonable settings of the problem, one can also expect that (cid:96) is continuous anddecreasing. By analysis of the constraints (3.8), we deduce that (cid:96) (0) = + ∞ and (cid:96) ( λ ) → λ → + ∞ . The mentioned properties of the function (cid:96) are sketched in Figure 1. (cid:45)(cid:54) (cid:114) λλ ∗ (cid:96) ( λ ) (cid:114) ζ ∗ (cid:96) enabling to compare the safety factors for the SSRand LA methods.By comparison of the constraints (3.2) and (3.8), we derive that the SSR-based safetyparameter λ ∗ defined in Section 3.1 is a unique solution of the following equation: (cid:96) ( λ ∗ ) = 1 . (3.9)It means that the safety factor λ ∗ can be found by the parametrization of the limit analysisproblem. From (3.9) and the properties of (cid:96) , we have:10 λ ∗ = 1 if and only if ζ ∗ = 1; • if λ ∗ > ζ ∗ > • if λ ∗ < ζ ∗ < λ ∗ is closer to one than ζ ∗ . One can expect that thedifference between λ ∗ and ζ ∗ is significant for φ ≈ ◦ when ζ ∗ is too large, see Section3.2.Finally of this section, we show that it suffices to parametrize only the friction angle withinthe LA method in order to compute the SSR factor λ ∗ . Indeed, consider a parameter χ ≥ (cid:101) (cid:96) ( χ ) = supremum of ζ ≥ subject to − div σ = F in Ω , σn = f on ∂ Ω f , ( σ − σ ) (cid:113) tan φχ + ( σ + σ ) tan φχ ≤ cζ in Ω . (3.10)Since the LA-based safety factor depends linearly on the cohesion, it is easy to prove that (cid:96) ( λ ) = 1 λ (cid:101) (cid:96) ( λ ) ∀ λ ≥ . (3.11)From (3.9) and (3.11), it follows that the SSR factor of safety λ ∗ can be found as a solutionof the fixed-point problem (cid:101) (cid:96) ( λ ∗ ) = λ ∗ . (3.12)One can suppose that the function (cid:101) (cid:96) is also decreasing in slope stability or similar geo-application. Under this assumption it is possible to prove λ ∗ is closer to one than ζ ∗ .Using the equations (cid:96) ( λ ∗ ) = 1 or (cid:101) (cid:96) ( λ ∗ ) = λ ∗ and a solver on the LA problem, one cancompute the safety factor λ ∗ of the SSR method. In fact, such a treatment was used in[Tschuchnigg et al. (2015a), Tschuchnigg et al. (2015b), Tschuchnigg et al. (2015c)]. InSection 6, we shall propose another numerical method for finding λ ∗ without the necessityto solve repeatedly the LA problem. 11 Modified Davis approach for the shear strength re-duction method in non-associative plasticity
In order to use the limit analysis (LA) method for the non-associative model with ψ < φ , itwas suggested [Davis (1968)] to reduce the strength parameters by the following empiricalformulas: c β := βc, φ β = ψ β := arctan( β tan φ ) , (4.1)where β := cos ψ cos φ − sin ψ sin φ , (4.2)and then solve the LA problem with c β and φ β = ψ β . Notice that β = 1 for ψ = φ . Inorder to apply a similar approach for the shear strength reduction (SSR) method, we usethe function (cid:96) introduced in Section 3.3 which parametrizes the LA problem with respectto the strength reduction factor λ . However, this parametrization causes that the factor β depends on λ and thus the modification of the Davis approach for purposes of the SSRmethod is not so straightforward.In [Tschuchnigg et al. (2015b)], three different modifications of the Davis approach havebeen proposed. The aim of this section is to analyze in more details two of these modifi-cations. To this end, we use the static principle presented in Section 3.1. In the first approach (called Davis A in [Tschuchnigg et al. (2015b)]), we assume that theparameter β defined by (4.2) remains fixed. We denote the corresponding FoS of the SSRmethod by λ ∗ I and define it as follows: λ ∗ I = supremum of λ ≥ subject to − div σ = F in Ω , σn = f on ∂ Ω f , ( σ − σ ) (cid:113) β tan φλ + ( σ + σ ) β tan φλ ≤ βcλ in Ω . (4.3)Using the function Φ (see (3.3)), one can also write12 ∗ I = supremum of λ ≥ subject to − div σ = F in Ω , σn = f on ∂ Ω f , Φ (cid:0) λβ ; σ (cid:1) ≤ . (4.4)If we compare the definitions of λ ∗ (see Section 3.1) and λ ∗ I , use the properties of Ψ andthe fact that β ≤
1, then we arrive at the expected inequality λ ∗ I ≤ λ ∗ (4.5)saying that FoS decreases with the decreasing dilatation angle ψ . Using the function (cid:96) defined in Section 3.3, one can find λ ∗ I as a solution of the following implicit equation: (cid:96) (cid:16) λ ∗ I β (cid:17) = 1 . (4.6) In the second approach (called Davis C in [Tschuchnigg et al. (2015b)]), we assume thatthe parameter β depends on the strength reduction factor λ ≥
0. This dependence canbe represented by the following function: b ( λ ) = cos ψ cos φ λ − sin ψ sin φ λ , if φ λ ≥ ψ, , if φ λ ≤ ψ, (4.7)where the function λ (cid:55)→ φ λ is defined by (2.13). Clearly, β = b (1). Further, it is easy toshow that the function b is continuous, nondecreasing and bounded by one from above.We denote the FoS of the SSR method for approach II by λ ∗ II and define it as follows: λ ∗ II = supremum of λ ≥ subject to − div σ = F in Ω , σn = f on ∂ Ω f , Φ (cid:0) λb ( λ ) ; σ (cid:1) ≤ . (4.8)The following statement is important to be the safety factor λ ∗ II well defined. Lemma 4.1. If (4.8) is satisfied for some λ := λ > then (4.8) holds for any λ < λ . roof. Using the formulas (2.8), we derive that λb ( λ ) = √ λ +tan φ − tan φ sin ψ cos ψ , if φ λ ≥ ψ,λ, if φ λ ≤ ψ, (4.9)implying that the function λ (cid:55)→ λ/b ( λ ) is continuous, increasing and nonnegative.Let ¯ λ > σ satisfying (4.8) for λ := λ >
0. For any λ < ¯ λ , wehave Φ (cid:0) λb ( λ ) ; σ (cid:1) ≤ Φ (cid:0) ¯ λb (¯ λ ) ; σ (cid:1) ≤ . Therefore, σ satisfies (4.8) for any λ < λ .If we compare the definitions of λ ∗ (see Section 3.1) and λ ∗ II , use the properties of Ψ andthe fact that b ( λ ) ≤ λ ≥
0, then we arrive at the expected inequality λ ∗ II ≤ λ ∗ . (4.10)Using the function (cid:96) defined in Section 3.3, one can find λ ∗ II as a solution of the followingimplicit equation: (cid:96) (cid:16) λ ∗ II b ( λ ∗ II ) (cid:17) = 1 . (4.11)One can also relate the safety factors λ ∗ I and λ ∗ II . We have: • λ ∗ I = 1 if and only if λ ∗ II = 1; • if λ ∗ I > λ ∗ I ≤ λ ∗ II ; • if λ ∗ I < λ ∗ I ≥ λ ∗ II . The numerical methods presented below are rather based on the kinematic principle to theSSR method. Therefore, we derive this dual formulation of the SSR problem. The deriva-tion is similar to the one in the limit analysis problem [Temam (1985), Christiansen (1990),14aslinger et al. (2019)]. Therefore, some steps of the derivation will be skipped, for thesake of simplicity.First, we unify the notation of the problems defining the safety factors λ ∗ , λ ∗ I and λ ∗ II inthe following way: ω ∗ = supremum of λ ≥ subject to − div σ = F in Ω , σn = f on ∂ Ω f , Φ (cid:0) q ( λ ); σ (cid:1) ≤ , (5.1)where q is an increasing, continuous and nonnegative function. Setting q ( λ ) = λ , q ( λ ) = λ/β and q ( λ ) = λ/b ( λ ), we obtain ω ∗ = λ ∗ , ω ∗ = λ ∗ I and ω ∗ = λ ∗ II , respectively.Next, we introduce the following functional spaces: V = { v ∈ [ H (Ω)] | v = 0 on ∂ Ω u } , (5.2)Σ = { σ ∈ [ L (Ω)] × | σ ij = σ ji in Ω } . (5.3)The space V represents velocity fields satisfying the Dirichlet boundary condition from(2.2) and Σ is used for symmetric stress fields. L (Ω) and H (Ω) denotes the Lebesgue andSobolev spaces, respectively, with respect to square integrable functions. More advancedfunctional spaces are considered in [Christiansen (1990)].Using the space V we arrive at the weak form of (5.1) : (cid:90) Ω σ : ε ( v ) dx = L ( v ) ∀ v ∈ V, (5.4)where L is the load functional defined by L ( v ) = (cid:90) Ω F · v dx + (cid:90) ∂ Ω f f · v ds. (5.5)To be the integrals in (5.5) meaningful, it suffices to assume that the components of F and f belong to the L space. Let Λ denote the set of stresses σ ∈ Σ satisfying (5.4) andlet P q ( λ ) := { σ ∈ Σ | Φ (cid:0) q ( λ ); σ (cid:1) ≤ } . (5.6)15e see that the set P q ( λ ) represents the constraint (5.1) and thus we can write ω ∗ = sup { λ ≥ | P q ( λ ) ∩ Λ (cid:54) = ∅} = sup λ ≥ sup σ ∈ P q ( λ ) ∩ Λ { λ } . (5.7)From (5.4), we haveinf v ∈ V (cid:20)(cid:90) Ω σ : ε ( v ) dx − L ( v ) (cid:21) = (cid:40) , if σ ∈ Λ , −∞ , otherwise . Hence, one ca rewrite (5.7) as follows: ω ∗ = sup λ ≥ sup σ ∈ P q ( λ ) inf v ∈ V (cid:20) λ + (cid:90) Ω σ : ε ( v ) dx − L ( v ) (cid:21) = sup λ ≥ inf v ∈ V sup σ ∈ P q ( λ ) (cid:20) λ + (cid:90) Ω σ : ε ( v ) dx − L ( v ) (cid:21) = sup λ ≥ inf v ∈ V (cid:20) λ + (cid:90) Ω D ( q ( λ ); ε ( v )) dx − L ( v ) (cid:21) , (5.8)where D ( q ( λ ); ε ) = sup σ ∈ R × sym Φ( q ( λ ); σ ) ≤ σ : ε (5.9)denotes the local dissipation function depending on λ ≥
0. The function D ( q ( λ ); · ) isfinite-valued only on a convex cone belonging to R × sym . Therefore, the inner problem in(5.8) can be classified as cone programming. The problem (5.8) can be interpreted as the kinematic principle to the SSR method.Let us note that the ordering of inf and sup has been interchanged during the derivationof (5.8). The corresponding equality is expected and partially supported by the resultspresented in [Christiansen (1990), Haslinger et al. (2019)]. In [Sysala et al. (2015), Cermak et al. (2015), Haslinger et al. (2016a), Haslinger et al. (2016b),Repin et al. (2018), Haslinger et al. (2019)], a regularization method has been systemat-ically developed for solution of the limit analysis (LA) problem. This methods has also16een used in strain-gradient plasticity [Reddy and Sysala (2020)]. The aim of this sectionis to use the regularization for solution of the shear strength reduction (SSR) problemand study the relation between the original and the regularized problem.We arise from (5.7) and regularize this problem with respect to a parameter α > ω ∗ α = sup λ ≥ sup σ ∈ P q ( λ ) ∩ Λ (cid:20) λ − α (cid:90) Ω C − σ : σ dx (cid:21) , (6.1)where C is a positive definite fourth order tensor, for example, the elastic tensor from(2.12). We have the following result. Lemma 6.1.
The sequence { ω ∗ α } α> defined by (6.1) is nondecreasing and satisfying ω ∗ α ≤ ω ∗ , lim α → + ∞ ω ∗ α = ω ∗ , (6.2) where ω ∗ is defined by (5.7) .Proof. It is readily seen that the inequalities ω ∗ α ≤ ω ∗ α ≤ ω ∗ hold for any 0 < α ≤ α .Next, for any λ < ω ∗ the intersection P q ( λ ) ∩ Λ is nonempty as follows from (5.7). Thenthe inner sup-problem in (6.1) has a unique solution σ λ ∈ P q ( λ ) ∩ Λ because it containsthe quadratic functional. Consequently,lim α → + ∞ ω ∗ α ≥ lim α → + ∞ (cid:20) λ − α (cid:90) Ω C − σ λ : σ λ dx (cid:21) = λ and thus ω ∗ ≥ lim α → + ∞ ω ∗ α ≥ sup { λ ≥ | P q ( λ ) ∩ Λ (cid:54) = ∅} = ω ∗ . This implies the limit in (6.2).One can also write ω ∗ α = max λ ≥ [ λ − G α ( λ )] = λ ∗ α − G α ( λ ∗ α ) , (6.3)where G α ( λ ) = inf σ ∈ P q ( λ ) ∩ Λ α (cid:90) Ω C − σ : σ dx (6.4)= (cid:40) α (cid:82) Ω C − σ λ : σ λ dx, if P q ( λ ) ∩ Λ (cid:54) = ∅ , + ∞ , otherwise , (6.5)17nd λ ∗ α maximizes the middle term in (6.3). Since the value G α ( λ ∗ α ) is finite, we have P q ( λ ∗ α ) ∩ Λ (cid:54) = ∅ . This fact, (5.7) and (6.3) imply the following result. Lemma 6.2.
The sequence { λ ∗ α } α> defined by (6.4) satisfies ω ∗ α ≤ λ ∗ α ≤ ω ∗ , lim α → + ∞ λ ∗ α = ω ∗ , (6.6) where ω ∗ and ω ∗ α are defined by (5.7) and (6.1) , respectively. In practice, we observe that the value λ ∗ α is sufficiently close to ω ∗ for much smaller α than ω ∗ α . Therefore, the sequence { λ ∗ α } α is more convenient for the approximation of ω ∗ than { ω ∗ α } α .Let us note that the 1D optimization problem in (6.3) can be solved, for example, bysequential enlarging of λ for fixed α . It is convenient to use a continuation over α tobe λ ∗ α sufficiently close to ω ∗ . The sufficiently large value of α can be found by thecontinuation on a coarse finite element mesh and then used on finer meshes.Next, for solution of (6.3), it is crucial to evaluate the function G α . To this end, we use asimilar duality approach as in Section 5. We arrive at the following kinematic definitionof G α : G α ( λ ) = − inf v ∈ V (cid:20)(cid:90) Ω D α ( q ( λ ); ε ( v )) dx − L ( v ) (cid:21) , (6.7)where D α ( q ( λ ); ε ) = sup σ ∈ R × sym Φ( q ( λ ); σ ) ≤ (cid:20) σ : ε − α C − σ : σ (cid:21) (6.8)is the regularized dissipative function. In particular, D α is finite-valued and differentiablewith respect to ε unlike the original dissipation D , see, for example, [Sysala (2014)].Moreover, the second derivative of D α exists almost everywhere. Let T α ( q ( λ ); ε ) ∈ R × sym denote the derivative of D α ( q ( λ ); ε ) with respect to ε . Then the problem (6.7) is equivalentto the following nonlinear variational equation:find v q ( λ ) ∈ V : (cid:90) Ω T α ( q ( λ ); ε ( v q ( λ ) )) : ε ( v ) dx = L ( v ) ∀ v ∈ V. (6.9)It is convenient to solve it by a non-smooth and damped version of the Newton methodsuggested in [Sysala (2012)] because this method also finds descent directions of the func-tional in (6.7) which need not be bounded from below for some λ .Let D := D α and T := T α for α = 1. Then the following formulas hold for any α > ≥ ε ∈ R × sym : D α ( q ( λ ); ε ) = 1 α D ( q ( λ ); α ε ) , T α ( q ( λ ); ε ) = T ( q ( λ ); α ε ) . (6.10)This formulas simplifies the construction of the operators D α and T α if the continuationover α is used. In addition, T is practically the same as the operator which arises fromthe implicit Euler discretization of the constitutive problem (2.10). Its construction canbe found e.g. in [de Souza Neto et al. (2011), Sysala et al. (2017)]. The closed form of D is introduced in Appendix of this paper, for the sake of completeness. In this section, we present two different numerical examples on slope stability problems.The first example is a model homogeneous slope presented in [Tschuchnigg et al. (2015b)].The aim is to illustrate our theoretical results and verified that the computed FoS arein accordance with the published ones. The second example arises from analysis of areal slope. A heterogeneous material and the influence of the water pore pressure areconsidered.
We use and compare the results from three different softwares: in-house codes in Matlab,Plaxis and Comsol Multiphysics.The in-house Matlab codes are based on elastic-plastic solvers, the finite element methodand on mesh adaptivity. They have been systematically developed and described in[Sysala et al. (2016), Sysala et al. (2017), Cermak et al. (2019)]. Some of the codes areavailable for download [Cermak et al. (2018)]. Within these codes, we have implementedthe regularization method from Section 6 to achieve safety factors for the associativeplasticity and for the Davis I and Davis II approaches to the nonassociative plasticity. Inparticular, P2 finite elements with the 7-point Gauss quadrature have been used and com-bined with the mesh adaptivity introduced in [Haslinger et al. (2019), Sysala et al. (2019)].The software Plaxis enables to solve the shear strength reduction method for both asso-ciative and nonassociative plasticity. Their standard solver is based on the implicit Eulertime discretization, P4 elements (in 2D) and on the arc-lenght method [Brinkgreve (2011)].19ne can also easily implement the Davis I approach there. Due to this software we areable to compare suggested and current approaches to the shear strength reduction (SSR)method.The software Comsol Multiphysics with its Geomechanical module is not specialized nei-ther on the shear strenght reduction method nor the arc-length method. On the otherhand, one can add a global equation to the elastic-plastic system of equations with re-spect to an unknown parameter and optimize it. Due to these facts, the SSR methodfor the associative plasticity and the Davis I, II modifications can be implemented there.Beside the regularization method, the standard incremental procedure for solution of theelastic-plastic problem has been used. P4 elements are considered.
Following [Tschuchnigg et al. (2015b)], we consider a homogeneous model slope depictedin Figure 2. Its inclination is 45 ◦ and sizes (in meters) are introduced in the figure. The(effective) friction angle φ is 45 ◦ , the (effective) cohesion c is 6.0 kPa and the unit weight γ is 20.0 kN/m . The dilatancy angle is either ψ = 0 ◦ (nonassociative case) or ψ = 45 ◦ (associative case). We have chosen such values of φ and ψ in order to highlight maximalpossible differences between the associative and nonassociative models or between thesuggested and current approaches to the SSR method.Next, we set the following values of the Young modulus and the Poisson ratio: E = 40MPa and ν = 0 .
3. For purposes of the regularization method, the value α = 1000 is used.This value is sufficiently large as follows from Figure 5. (cid:64)(cid:64)(cid:64)(cid:64)(cid:64)(cid:63)(cid:63)(cid:63)(cid:63) F = − γ e (cid:52) (cid:52) (cid:52) (cid:52) (cid:52) (cid:52) (cid:52) (cid:52) (cid:52) (cid:52) (cid:52) (cid:52) (cid:52) (cid:52) (cid:52) (cid:100)(cid:100)(cid:100)(cid:100)(cid:100)(cid:100)(cid:100)(cid:100) (cid:100)(cid:100)(cid:100)(cid:100)
15 10 151010 10
Figure 2: Geometry of the model slope with inclination 45 ◦ .Table 1 summarizes computed factors of safety by different softwares and different ap-proaches. These values are practically the same as in [Tschuchnigg et al. (2015b)]. The20ast line contains the value of the current approach to the SSR method. The correspond-ing safety factor cannot be uniquely determined for finer meshes due to the oscillation ofthe method, see [Tschuchnigg et al. (2015b)] for more details.Table 1: Safety factors for the homogeneous slope and different approachesMATLAB COMSOL Plaxis ψ = φ , assoc. model 1.53 1.52 1.51 ψ = 0 ◦ , Davis I 1.08 1.08 1.08 ψ = 0 ◦ , Davis II 1.15 1.16 – ψ = 0 ◦ , current approach – – 1.27–1.35On the other hand, the results for the associative plasticity and for the Davis approachesare practically insensitive if sufficiently fine meshes are used. It is illustrated in Figure3 where the mesh adaptivity within the in-house Matlab codes is used. In particular, 20levels of meshes is considered and the corresponding safety factors remain almost constant.The finest mesh and the corresponding slip zone for the Davis II approach are depictedin Figure 4. The slip zone is visualized using the deviatoric strain.The dependence of λ ∗ α and ω ∗ α (see Section 6) for the Davis II approach on the regular-ization parameter α is depicted in Figure 5. We see that these curves are increasing andtending to the limit value ω ∗ . It confirms the theoretical results of the regularizationmethod. We also see that λ ∗ α approximates the safety factor ω ∗ even for relatively small α . On the other hand, much larger values of α has to be used in order to approximate ω ∗ by ω ∗ α . The second example follows from analysis of a real slope in locality Doubrava-Kozinec(near Karvina in the North-East part of the Czech Republic) where a landslide has beenobserved. In Figure 6, the investigated slope profile is illustrated. We see that the slopeis heterogeneous. Particular materials and their parameters are specified in Table 2. Thelevel of groundwater is estimated by the blue curve in the figure. We distinguish thespecific weights γ unsat and γ sat for unsaturated and saturated materials, respectively.The values of the dilatancy angle have not been available for us, therefore, we set ψ = 0 ◦ for all materials. Nevertheless, we also consider the associative case with ψ = φ in orderto analyze the influence of the dilatancy. 21igure 3: Safety factors for the homogeneous slope depending on mesh adaptivity.Figure 4: The finest mesh and the corresponding slip zone for the homogeneous slope.Figure 5: Dependence of λ ∗ α and ω ∗ α on the regularization parameter α .22able 2: Material parameters for the heterogeneous slope. neogene clay gravel quaternary clay sand clayed sand φ [ ◦ ] 26 45 13 33 27 c [kPa] 9 1 3 2 5 E [MPa] 16 140 10 14 27 ν γ unsat [kN/m ] 20.3 20.5 20.0 19.0 19.4 γ sat [kN/m ] 20.7 20.6 20.5 20.5 21.4The landslide has been observed in the central part of the slope, mainly in quaternary clayand in its interface with sand clay and neogene clay. It is worth noticing that the valueof the friction angle of the quaternary clay is much less than for the remaining materials.Numerical results presented below confirm the location of the landslide.Figure 6: Geometry of the heterogeneous slope.Computed safety factors for different approaches are listed in Table 3. We see that thesafety factors are very close to 1 for all investigated approaches and thus the slope isunstable. The values received by the in-house Matlab codes are slightly lower than theremaining ones due to the usage of the local mesh adaptivity.The initial mesh for computation in Comsol Multiphysics is depicted in Figure 7. Thismesh reflects the material heterogeneity and has been also imported to Matlab. In ComsolMultiphysis, this mesh was then locally refined in the central part of the slope to achievemore accurate results. In Matlab, the original mesh was adaptively refined and 15 meshlevels were considered.A detail of the finest Matlab mesh is depicted in Figure 8 together with the corresponding23able 3: Safety factors for the heterogeneous slope and different approachesMATLAB COMSOL Plaxis ψ = φ , assoc. model 1.05 1.09 1.08 ψ = 0 ◦ , Davis I 1.02 1.06 1.05 ψ = 0 ◦ , Davis II 1.02 1.06 – ψ = 0 ◦ , current approach – – 1.06Figure 7: Initial mesh for the heterogeneous slope created in Comsol Multiphysics.slip zone. For visualization of this zone, a norm of the deviatoric strain was used. Theblack curves in the figure depicts the material interface. We see that the zone is not toodeep and its large part lies on the interface of the quaternary and neogene clays.Figure 8: Detail of the finest mesh in Matlab and the corresponding slip zone on thematerial interface.In Figure 9, we see the dependence of the safety factors on mesh adaptivity computed inMatlab for all three approaches. We observe that these curves are practically constantafter a sufficiently large numbers of mesh refinements and thus the received values arealmost independent of the mesh density. We also observe that the safety factor for theDavis I approach is slightly less that the one for Davis II. This is in accordance with our24heoretical results.Figure 9: Safety factors for the homogeneous slope depending on mesh adaptivity. This work has been inspired by recent ideas presented in [Tschuchnigg et al. (2015b)]where the shear strength reduction (SSR) method was interpreted as a parametrizedlimit analysis (LA) method. Based on these ideas, we have built in this paper a rigorous variant of the SSR method, analyzed and used it in geotechnics, independently of the LAmethod. The suggested variant is based on optimization problems, variational principlesand on modifications of the Davis approach.For numerical solution, a regularization method has been used and combined with thefinite element and the damped Newton method. This solution concept leads to similarsolvers which are standardly used in computational plasticity and thus it can be easilyimplemented within existing elastic-plastic codes. In particular, we have arised fromour in-house Matlab codes [Cermak et al. (2019)] and completed them with local meshadaptivity. Softwares Plaxis and Comsol Multiphysics have been used for comparison ofthe results. One of the presented numerical examples has followed from analysis of a realslope.
Acknowledgment:
The authors acknowledge support for their work from the CzechScience Foundation (GA ˇCR) through project No. 19-11441S. The authors also thank to25r. Alexej Kolcun for fruitful discussions on mesh adaptivity for regular and irregularmeshes. D D ( q ( λ ); ε ) = sup σ ∈ R × sym Φ( q ( λ ); σ ) ≤ (cid:20) σ : ε − C − σ : σ (cid:21) . (9.1)from Section 6. In literature (see, for example, [de Souza Neto et al. (2011), Sysala et al. (2017)]),one can find closed form of the derivative T of D representing the stress-strain relationbut not D . Therefore, we try to fill this gap.Beside λ ≥
0, this function also depends on the inelastic parameters c , φ and the elasticparameters K , G . We construct the function D ( q ( λ ); ε ) only for such λ satisfying q ( λ ) =1. For other choices of λ , it suffices to replace c and φ with c q ( λ ) := cq ( λ ) , φ q ( λ ) := arctan tan φq ( λ ) . (9.2)see also (2.13). To be in accordance with the derivation presented in [Sysala et al. (2017)],we shall write the yield criterion Φ(1; σ ) ≤ φ ) σ − (1 − sin φ ) σ − c cos φ ≤ . (9.3)Next, we assume that the strain tensor ε is given and its eigenvalues ε , ε , ε satisfy ε ≥ ε ≥ ε . Let tr ε = ε + ε + ε denote the trace of ε and Λ = (3 K − G ) denotethe first Lam´e coefficient. As in the book [de Souza Neto et al. (2011)], we distinguishfife possible cases: the elastic response, the return to the smooth portion of the Mohr-Coulomb pyramid, the return to the left edge, the return to the right edge, and the returnto the apex of the pyramid. The elastic response.
This case happen if the elastic stress C ε satisfies Φ(1; C ε ) ≤ ε ) sin φ + 2 G (1 + sin φ ) ε − G (1 − sin φ ) ε − c cos φ ≤ . (9.4)26hen D (1; ε ) = 12 C ε : ε = 12 Λ(tr ε ) + G ( ε + ε + ε ) . (9.5)If the criterion (9.4) does not hold then the plastic response occurs and we distinguishfour possible cases of the return to the Mohr-Coulomb pyramid. We use the followingauxiliary notation: γ s,l = ε − ε φ , γ s,r = ε − ε − sin φ ,γ l,a = ε + ε − ε − sin φ , γ r,a = 2 ε − ε − ε φ . The return to the smooth portion.
This case happen if (9.4) is not satisfied and q s ( ε ) < S min { γ s,l , γ s,r } , (9.6)where q s ( ε ) = 2Λ(tr ε ) sin φ + 2 G (1 + sin φ ) ε − G (1 − sin φ ) ε − c cos φ,S = 4Λ sin φ + 4 G (1 + sin φ ) . Then, D (1; ε ) = 12 Λ(tr ε ) + G ( ε + ε + ε ) − S q s ( ε ) . (9.7) The return to the left edge.
This case happen if (9.4) is not satisfied and γ s,l < γ l,a , Lγ s,l ≤ q l ( ε ) < Lγ l,a , (9.8)where q l ( ε ) = 2Λ(tr ε ) sin φ + G (1 + sin φ )( ε + ε ) − G (1 − sin φ ) ε − c cos φ,L = 4Λ sin φ + G (1 + sin φ ) + 2 G (1 − sin φ ) . Then, D (1; ε ) = 12 Λ(tr ε ) + G (cid:104)
12 ( ε + ε ) + ε (cid:105) − L q l ( ε ) . (9.9)27 he return to the right edge. This case happen if (9.4) is not satisfied and γ s,r < γ r,a , Rγ s,r ≤ q r ( ε ) < Rγ r,a , (9.10)where q r ( ε ) = 2Λ(tr ε ) sin φ + 2 G (1 + sin φ ) ε − G (1 − sin φ )( ε + ε ) − c cos φ,R = 4Λ sin φ + 2 G (1 + sin φ ) + G (1 − sin φ ) . Then, D (1; ε ) = 12 Λ(tr ε ) + G (cid:104) ε + 12 ( ε + ε ) (cid:105) − R q r ( ε ) . (9.11) The return to the apex.
This case happen if (9.4) is not satisfied and q a ( ε ) ≥ A max { γ l,a , γ r,a } , (9.12)where q a ( ε ) = 2 K (tr ε ) sin φ − c cos φ, A = 4 K sin φ. Then, D (1; ε ) = 12 K (tr ε ) − A q a ( ε ) = c tan φ (tr ε ) − c K tan φ . (9.13) References [Brinkgreve and Bakker (1991)] Brinkgreve, R. B. J. & Bakker, H. L. (1991). Non-linear finiteelement analysis of safety factors. Proceedings of the international conference on com-puter methods and advances in geomechanics, pp. 1117–1122. Rotterdam, the Netherlands:Balkema.[Brinkgreve (2011)] Brinkgreve, R. B. J., Swolfs, W. M. & Engin, E. (2011). Plaxis 2D 2011 –user manual. Delft, the Netherlands: Plaxis bv.[Cermak et al. (2015)] Cermak, M., Haslinger, J., Kozubek, T., Sysala, S. (2015). Discretiza-tion and numerical realization of contact problems for elastic-perfectly plastic bodies. PARTII–numerical realization, limit analysis.
ZAMM-Journal of Applied Mathematics and Me-chanics/Zeitschrift f¨ur Angewandte Mathematik und Mechanik , 95(12), 1348–1371.[Cermak et al. (2018)] ˇCerm´ak, M., Sysala, S., & Valdman, J. (2018). MATLAB FEM packagefor elastoplasticity, https://github.com/matlabfem/matlab fem elastoplasticity. Cermak et al. (2019)] ˇCerm´ak, M., Sysala, S., & Valdman, J. (2019). Efficient and flexibleMATLAB implementation of 2D and 3D elastoplastic problems. Applied Mathematics andComputation, 355, 595–614.[Chen and Liu (1990)] Chen, W. and Liu, X.L. (1990).
Limit Analysis in Soil Mechanics . Else-vier.[Christiansen (1990)] Christiansen, E. (1996).
Limit analysis of colapse states . In P. G. Ciarletand J. L. Lions, editors,
Handbook of Numerical Analysis , Vol IV, Part 2, North-Holland,195–312.[Davis (1968)] Davis, E. H. (1968). Theories of plasticity and failure of soil masses. In Soilmechanics: selected topics (ed. I. K. Lee), pp. 341–354. New York, NY, USA: Elsevier.[Dawson et al. (1999)] Dawson, E. M., Roth, W. H. & Drescher, A. A. (1999).Slope stability analysis by strength reduction. G´eotechnique 49, No. 6, 835–840,http://dx.doi.org/10.1680/geot.1999.49.6.835.[Duncan (1996)] Duncan, J. M. (1996). State of the art: limit equilibrium and finite-elementanalysis of slopes. Journal of Geotechnical engineering, 122(7), 577–596.[de Souza Neto et al. (2011)] de Souza Neto, E. A., Peric, D., & Owen, D. R. (2011). Compu-tational methods for plasticity: theory and applications. John Wiley & Sons.[Ekeland and Temam (1974)] Ekeland, I. and Temam, R. (1974).
Analyse Convexe et Probl`emesVariationnels . Dunod, Gauthier Villars, Paris.[Griffiths and Lane (1999)] Griffiths, D. V. & Lane, P. A. (1999). Slope stabil-ity analysis by finite elements. G´eotechnique 49, No. 3, 387–403, http://dx.doi.org/10.1680/geot.1999.49.3.387.[Hamlaoui et al. (2017)] Hamlaoui, M., Oueslati, A., & De Saxc´e, G. (2017). A bipotentialapproach for plastic limit loads of strip footings with non-associated materials. InternationalJournal of Non-Linear Mechanics, 90, 1–10.[Han and Reddy (2012)] Han, W., & Reddy, B. D. (2012). Plasticity: mathematical theory andnumerical analysis (Vol. 9). Springer Science & Business Media.[Haslinger et al. (2016a)] Haslinger, J., Repin, S., Sysala, S. (2016). A reliable incrementalmethod of computing the limit load in deformation plasticity based on compliance: Contin-uous and discrete setting.
Journal of Computational and Applied Mathematics
Applications of Mathematics
61, 527–564. Haslinger et al. (2019)] Haslinger, J., Repin, S., Sysala, S. (2019). Inf-sup conditions on convexcones and applications to limit load analysis.
Mathematics and Mechanics of Solids
Comp. & Math. with Appl. (2018) 75: 199–217.[Sloan (2013)] Sloan SW (2013). Geotechnical stability analysis,
G´eotechnique , 63, 531–572.[Sysala (2012)] Sysala, S. (2012). Application of a modified semismooth Newton method to someelasto-plastic problems. Mathematics and Computers in Simulation, 82(10), 2004–2021.[Sysala (2014)] Sysala, S. (2014). Properties and simplifications of constitutive time-discretizedelastoplastic operators. ZAMM-Journal of Applied Mathematics and Mechanics/Zeitschriftf¨ur Angewandte Mathematik und Mechanik, 94(3), 233–255.[Sysala et al. (2019)] Sysala, S., Blaheta, R., Kolcun, A., ˇSˇcuˇcka, J., Souˇcek, K., & Pan, P. Z.(2019). Computation of Composite Strengths by Limit Analysis. Key Engineering Materi-als, 810, 137–142.[Sysala et al. (2017)] Sysala, S., ˇCerm´ak, M., & Ligursk´y, T. (2017). Subdifferential-basedimplicit return-mapping operators in Mohr-Coulomb plasticity. ZAMM-Journal of Ap-plied Mathematics and Mechanics/Zeitschrift f¨ur Angewandte Mathematik und Mechanik,97(12), 1502–1523.[Sysala et al. (2015)] Sysala, S., Haslinger, J., Hlav´aˇcek, I., Cermak, M. (2015). Discretiza-tion and numerical realization of contact problems for elastic-perfectly plastic bodies.PART I–discretization, limit analysis.
ZAMM-Journal of Applied Mathematics and Me-chanics/Zeitschrift f¨ur Angewandte Mathematik und Mechanik , 95(4), 333–353. Sysala et al. (2016)] Sysala, S., Cermak, M., Koudelka, T., Kruis, J., Zeman, J., & Blaheta,R. (2016). Subdifferential-based implicit return-mapping operators in computational plas-ticity. ZAMM-Journal of Applied Mathematics and Mechanics/Zeitschrift f¨ur AngewandteMathematik und Mechanik, 96(11), 1318–1338.[Temam (1985)] Temam, R. (1985).
Mathematical Problems in Plasticity . Gauthier-Villars,Paris.[Tschuchnigg et al. (2015a)] Tschuchnigg, F., Schweiger, H.F., Sloan, S.W., Lyamin, A.V., &Raissakis, I. (2015). Comparison of finite-element limit analysis and strength reductiontechniques. G´eotechnique, 65(4), 249–257.[Tschuchnigg et al. (2015b)] Tschuchnigg, F., Schweiger, H.F., & Sloan, S.W. (2015). Slope sta-bility analysis by means of finite element limit analysis and finite element strength reductiontechniques. Part I: Numerical studies considering non-associated plasticity. Computers andgeotechnics, 70, 169–177.[Tschuchnigg et al. (2015c)] Tschuchnigg, F., Schweiger, H.F., & Sloan, S.W. (2015). Slope sta-bility analysis by means of finite element limit analysis and finite element strength reduc-tion techniques. Part II: Back analyses of a case history. Computers and geotechnics, 70,178–189.[Yu (2006)] Yu, H.-S. (2006). Plasticity and Geotechnics, Springer Science+Bussiness Media,New York.[Zienkiewicz et al. (1975)] Zienkiewicz, O.C., Humpheson, C., and Lewis, R.W. (1975). Asso-ciated and non-associated visco-plasticity and plasticity in soil mechanics. G´eotechnique,25(4), 671–689.. Gauthier-Villars,Paris.[Tschuchnigg et al. (2015a)] Tschuchnigg, F., Schweiger, H.F., Sloan, S.W., Lyamin, A.V., &Raissakis, I. (2015). Comparison of finite-element limit analysis and strength reductiontechniques. G´eotechnique, 65(4), 249–257.[Tschuchnigg et al. (2015b)] Tschuchnigg, F., Schweiger, H.F., & Sloan, S.W. (2015). Slope sta-bility analysis by means of finite element limit analysis and finite element strength reductiontechniques. Part I: Numerical studies considering non-associated plasticity. Computers andgeotechnics, 70, 169–177.[Tschuchnigg et al. (2015c)] Tschuchnigg, F., Schweiger, H.F., & Sloan, S.W. (2015). Slope sta-bility analysis by means of finite element limit analysis and finite element strength reduc-tion techniques. Part II: Back analyses of a case history. Computers and geotechnics, 70,178–189.[Yu (2006)] Yu, H.-S. (2006). Plasticity and Geotechnics, Springer Science+Bussiness Media,New York.[Zienkiewicz et al. (1975)] Zienkiewicz, O.C., Humpheson, C., and Lewis, R.W. (1975). Asso-ciated and non-associated visco-plasticity and plasticity in soil mechanics. G´eotechnique,25(4), 671–689.