A new design equation for drained stability of conical slopes in cohesive-frictional soils

Boonchai Ukritchon,Suraparb Keawsawasvong

Geotechnical Research Unit,Department of Civil Engineering,Faculty of Engineering,Chulalongkorn University,Bangkok,10330,Thailand

1.Introduction

Various civil construction projects involve slope stability problem,such as natural,filled and cut slopes.Determination of the factor of safety for these slopes is the most important aspect to ensure an adequate safety of the slopes in both short and long terms.Due to the practical importance,slope stability has drawn much attention from many investigators,including two dimensional slope stability under a plane-strain condition(Griffiths and Koutsabeloulis,1985;Yu et al.,1998;Michalowski,2002;Li,2006;Martin,2011;Griffiths and Yu,2015),and three dimensional(3D)slope stability(Griffiths and Marquez,2007;Li et al.,2009,2010;Michalowski,2010;Xie et al.,2011;Lu,2015;Lim et al.,2016;Kelesoglu,2016).In many practical and theoretical considerations of the problem,the plane-strain condition is commonly assumed in a conventional analysis of slope stability(e.g.Zheng et al.,2009;Chen and Huang,2011;Yang et al.,2011;Li et al.,2014;Yao et al.,2014;Ghanbari and Hamidi,2015).However,an inclined slope in an axisymmetric condition,called conical slope,can be generally found in various practical applications.For example,unsupported vertical circular excavations are fairly common especially in the construction of cast in situ piles or piers.In addition,an unsupported square excavation with all-around side slope may be reasonably approximated as a conical slope.A cut slope with a rectangular shape may also be simpli fi ed as a conical slope with an equivalent radius that produces the same crosssectional area of the rectangular slope.Furthermore,conical slopes can be employed to construct underground conical pits for energy storage applications(Dincer,2010;Lee,2013),known as pit thermal energy storages(Ochs et al.,2009;Schmidt and Miedaner,2012),and underground structures,such as conical concrete storage tanks for molten salt solar power plants(Salomoni et al.,2008).

A number of studies for unsupported axisymmetric excavation were performed in the past to study the stability problems,especially unsupported vertical circular excavations in both cohesive and cohesive-frictional soils.These include analytical lower bound(LB)and upper bound(UB)analyses for homogeneous clays(Britto and Kusakabe,1982;Pastor and Turgeman,1982),and LB and UB finite element limit analysis(FELA)for linearly increasing cohesive(Khatri and Kumar,2010;Kumar et al.,2014)and cohesive frictional soils(Kumar and Chakraborty,2012;Kumar et al.,2014).Recently,Keawsawasvong and Ukritchon(2017a)investigated the unsupported conical excavation in both homogeneous and nonhomogeneous clays using FELA for a wide range of height ratios,dimensionless strength gradients and slope inclination angles,while a closed-form approximate expression was proposed for the LB solutions to predict the stability factor for an unsupported conical excavation in a generalized linearly increasing undrained strength profile.Till now,there have been very few studies of the conical slope stability in cohesive-frictional soils in the literature,except for the special case of unsupported vertical circular excavations,as mentioned earlier(Kumar and Chakraborty,2012;Kumar et al.,2014).

This paper extends the previous study of Keawsawasvong and Ukritchon(2017a)for a conical slope in an undrained material to cover the general case of the problem in a drained material.In this study,OptumG2(Krabbenhoft et al.,2015),the computational limit analysis using FELA for axisymmetric problems,was employed to investigate the drained stability of unsupported conical slopes in cohesive-frictional soils.The OptumG2 was chosen based on several successful investigations of the authors on a variety of stability problems in geotechnical engineering under both planestrain and axisymmetric conditions,including planar suction caissons(Keawsawasvong and Ukritchon,2016a),cylindrical suction caissons(Ukritchon and Keawsawasvong,2016),circular shallow foundations(Ukritchon and Keawsawasvong,2017a),laterally loadedpiles(Keawsawasvongand Ukritchon,2016b,2017b;Ukritchon and Keawsawasvong,2017b,2017d),limiting pressure ofsoilgapsbetween stabilizingpiles(Keawsawasvong and Ukritchon,2017d;Ukritchon and Keawsawasvong,2017c),active trapdoors(Keawsawasvong and Ukritchon,2017c),cantilever f l ood walls(Keawsawasvong and Ukritchon,2017e),and opening in underground walls(Ukritchon and Keawsawasvong,2017c).Results of the parametric study of drained conical slopes were summarized in the form of the dimensionless parameters,including stability factor,slope height ratio,soil friction angle and slope inclination angle.The influences of these input parameters on the predicted failure mechanism of conical slopes were discussed and compared.A new design equation developed from a nonlinear regression of the LB solutions of conical slope was proposed as a convenient tool for a stability analysis of this problem in practice.An application of the proposed equation was demonstrated through examples of a conical slope with both associated and non-associated f l ow rules.

2.Description of numerical models in FELA

Fig. 1 shows the problem definition of an unsupported conical slope with height(H),radius at the bottom of slope(b)and slope inclination angle(β).The soil profile corresponds to a homogeneous cohesive-frictional soil with a constant unit weight(γ),effective cohesion(c)and effective soil friction angle(φ).The cohesive-frictional soil corresponds to an intermediate soil,and is modeled by a perfectly plastic Mohr-Coulomb material with the associated f l ow rule.Based on the dimensionless technique,it is shown that the stability factor(N)of conical slopes can be represented by three input parameters as

Fig. 1.Problem notation of an unsupported conical slope.

whereγH/cis the stability factor andH/bis the slope height ratio.

Note that the stability factorN=γH/cwas adopted by various researchers to develop slope stability charts under plane-strain condition.This parameter was also employed to represent the stability criterion for other geometric shapes of slopes in cohesive frictional soils,i.e.unsupported circular excavations(axisymmetric problem)(Kumarand Chakraborty,2012),and 3D slopes(Michalowski,2010).Consequently,the factorNserves as the stability criterion for slope problems under both plane-strain and axisymmetric conditions as well as fully 3D geometry.For a conical slope,the factorNis a function of slope height ratioH/b,effective soil friction angleφ,and slope inclination angleβ,according to Eq.(1).However,the stability factor under the axisymmetric condition is significantly higher that under the plane-strain condition due to the 3D effect of slope shape among the two problems.Numerically derived stability factors presented later are useful for determination of the maximum slope height or maximum slope inclination angle of a conical slope with known strength parameters of soils(e.g.candφ).

For the present study,the problem is concerned with an analysis of conical slope under drained loading conditions.Thus,the shear strength parameters of soils,i.e.candφ,correspond to the effective(drained)stress parameters.Consequently,the solutions to stability factor of conical slope are obtained by the drained effective stress analysis,and can be applied to long-term stability of the problem after construction.Note that the effective stress analysis is generally employed to develop the practical stability charts of drained slopes under plane-strain,axisymmetric(Kumar and Chakraborty,2012),as well as 3D(Michalowski,2010)conditions.Thus,the stability factorNis introduced in this study to represent a stability criterion of conical slopes under drained loading conditions.FELA(e.g.Sloan,2013),the numerical method based on the plastic bound theorems and finite element concept,is adopted in this study,since it is a convenient and powerful tool in solving stability problems in geotechnical engineering,and has been advocated by many researchers to study stability of slopes under undrained and drained loading conditions(Yu et al.,1998;Li et al.,2009,2010;Kumar and Chakraborty,2012;Sloan,2013;Kumar et al.,2014;Lim et al.,2015a,2015b,2016;Keawsawasvong and Ukritchon,2017a).

In this paper,OptumG2(Krabbenhoft et al.,2015),the computational limit analysis using FELA,was employed to investigate the drained stability factor of conical slopes in cohesive-frictional soils.This numerical approach combines the powerful capabilities of finite element discretization for handling complicated soil strati fications,loadings and boundary conditions with the plastic bound theorems tobracket the exact limit load byUB and LB solutions.The underlying bound theorems assume a perfectly plastic material with an associated fl ow rule.Sloan(2013)clearly demonstrated that FELA is a powerful numerical technique for analyzing various complex stability problems in geotechnical engineering.Historical developments of FELA and its applications for various practical stability problems can be found in Sloan(2013).The following summarizes the method of computational limit analysis used in this study regarding the stability analysis of conical slopes.

In the UB analysis,six-noded triangular elements are employed to discretize a conical slope.Two unknown velocities are associated with each node,and they have a quadratic variation within triangular element.Like a conventional displacement based finite element method,jump of velocities across adjacent triangular element is not permitted,and hence the unknown velocities are continuous across adjacent elements.The UB calculation of conical slope is formulated as the optimization problem in the form of second-order conic programming(SOCP)(Krabbenhoft et al.,2007a,2007b;Makrodimopoulos and Martin,2007),in which constraints on kinematically admissible velocity field are enforced.They are defined as the compatibility equations with the associated fl ow rule on average for six-noded triangular elements,and under velocity boundary conditions.The objective function of the UB SOCP is to minimize the soil unit weight of conical slope,and is obtained through invoking the principle of virtual work by equating the rate of work done by external loads(i.e.soil unit weight)with the internal energy dissipation at triangular elements.

In theLB analysis,three-noded triangular elements are employed to discretize a conical slope.Four unknowns of stresses under axisymmetric condition are associated with each node,and they have a linear variation within triangular element.In contrast to the UB formulation,a stress field in the LB formulation is continuous within triangular element using a piecewise linear variation,but is discontinuous across shared edges of adjacent elements.The latter is modeled by assigning each element to employ its own unique nodes.The LB calculation of conical slope is formulated as the optimization problem in the form of SOCP(Makrodimopoulos and Martin,2006;Krabbenhoft et al.,2007a),in which constraints on kinematically admissible stress field are enforced.They are defined as equilibrium equations at the centroid of triangular elements and along stress discontinuities,under stress boundary conditions,and there is no violation of the yield criterion for all nodes.The objective function of the LB SOCP is to maximize the soil unit weight of conical slope,which is enforced to satisfy the vertical equilibrium equation of triangular elements.Note that the axisymmetric formulation of the LB SOCP in OptumG2 is comparable to that in Tang et al.(2014).

Both the UB and LB SOCP optimization problems in OptumG2 were solved to yield UB and LB solutions by using the general purpose optimization solver,SONIC(Krabbenhoft et al.,2006),which is based on the interior-point method(Wright,2005;Krabbenhoft et al.,2007b).Full details of the FELA formulation in this study can be found in Krabbenhoft et al.(2015).

Three parametric studies of conical slopes were performed to cover the practical ranges of this problem,i.e.H/b=0.5-10,φ =10°-40°and β =45°-90°.Because the geometry of the problem is a conical shape,an axisymmetric condition was used in FELA,where only one half of the domain was considered in the analysis,as shown in Fig. 1.The soil slope was discretized as triangular elements and modeled as Mohr-Coulomb material with the associated f l ow rule.Note that the influence of the nonassociated f l ow rule on the solution of conical slopes will be described in Section 4.The side of the conical slope and the top of ground surface correspond to a free surface,where movements of the soil along these boundaries were freely permitted in both the horizontal and vertical directions.The boundary conditions of this problem were such that the left(axial symmetry)and right boundaries of the problem were allowed to move only in the vertical direction,while the bottom boundary of the model was fixed in both directions.The size of the problem domains was chosen to be so large that the plastic yielding zone was contained within the domain and did not intersect the right and bottom boundaries.Therefore,the UB and LB limit loads are not altered by an extension of the domain size,and there is no influence of the domain size on the computed solutions.

An automatically adaptive mesh refinement,a powerful feature in OptumG2,was employed in both the UB and LB analyses to determine the tight UB and LB solutions.Five iterations of adaptive meshing with the number of elements increasing from 5000 to 10,000 were used for all analyses,since initial calibration tests revealed that this setting was adequate enough to obtain a suitable accurate solution.

It is worth noting that all UB and LB axisymmetric formulations in the literature(Khatri and Kumar,2010;Kumar and Chakraborty,2012,2014;Chakraborty and Kumar,2014;Kumar et al.,2014;Tang et al.,2014)as well as OptumG2(Krabbenhoft et al.,2015)provide only estimates of the rigorous solutions,since their requirements of the f l ow rule and equilibrium are only satisf i ed on average over triangular elements.However,those previous studies on the axisymmetric FELA showed that these approximations hardly introduce any significant error to the bound solutions.In addition,recent applications employing OptumG2 to the axisymmetric stability problems,i.e.undrained stability of conical slopes(Keawsawasvong and Ukritchon,2017a),pullout capacity of cylindrical suction caissons(Ukritchon and Keawsawasvong,2016),and circular shallow foundation(Ukritchon and Keawsawasvong,2017a)also confirm such finding.

3.Results and comparisons

Fig. 2 shows the relationship betweenγH/candH/bof the drained stability of conical slopes at different conical slope angles(β =45°-90°)and soil friction angles(φ =10°-40°).In all cases,the exact stability factors of conical slopes can be accurately bracketed within 1%of the average between the computed UB and LB solutions.Thus,it is very difficult to see the differences between bound solutions in Fig. 2.It should be noted that for a given same set of input dimensionless parameters(H/b,β,and φ),all UB solutions are generally larger than the corresponding LB solutions,due to the rigorous bounding property of the plastic bound theorems in which the exact solution to stability problem is bracketed from above and below using UB and LB solutions.Very accurate bound solutions are achieved,because the numerical feature of automatically adaptive mesh refinement provided in OptumG2 was employed,which could significantly reduce the differences between UB and LB solutions.In all simulations,five iterations of adaptive meshing with the number of elements increasing from 5000 to 10,000 were set up using the shear dissipation control.In general,for every step of adaptive meshing,the number of elements was automatically increased in the area requiring mesh refinement with high plastic shearing strain.Consequently,very accurate UB and LB solutions to the stability factors of conical slopes in cohesive-frictional soils were solved in the present study,where they bracket the exact solutions very closely,as shown in Fig. 2.A nonlinear relationship betweenγH/candH/bwas observed,while a linear relationship of these parameters was found forβ=90°and φ=30°-40°.The influence of bothβandφonγH/cwas manifested bya significant increase inγH/cin a nonlinear fashion.As expected,for a certain excavated height,a narrower conical slope(largerH/b)and a higher friction angle result in a higher stability factor.Note that comparisons in γH/cof vertical circular excavations(β =90°)between the computed bound solutions and existing solutions by FELA(Kumar and Chakraborty,2012;Kumar et al.,2014)are also shown in Fig. 2a.Excellent agreement inγH/cbetween the present study and the previous solutions(Kumar and Chakraborty,2012;Kumar et al.,2014)were observed in all cases.Note that the stability factorγH/cis much larger than the others when the slope inclination angle is 45°.For example,when φ =30°,N=44.76-88.92 forβ=45°,andN=8.88-43.82 forβ=90°,where the former is larger than the latter at around 202.93%-526.33%.This result is reasonable,since when a conical slope has moderate inclination angle(β =45°),its driving pressure caused by gravity becomes less significant,as compared to a steep one(β =90°).As a result,the conical slope with moderate inclination angle is more stable,or has larger stability factor than that with steep one.

Fig. 2.Relationships between the stability factor and slope height ratio at different slope angles(β)of(a)90°,(b)75°,(c)60° and(d)45°.

The influences ofβ,H/band φ on the predicted failure mechanism deduced from the final adaptive mesh are illustrated in Figs.3-5,respectively.Generally,two failure modes are observed,i.e.the toe failure and the face failure.The toe failure was commonly encountered for small soil friction angle(Figs.3-5a-c).On the other hand,a face failure generally occurred for a conical slope with a larger soil friction angle(Fig. 5d).The smaller slope angle(β),larger slope height ratio(H/b)and larger friction angle(φ)resulted in a decreased plastic yielding zone of conical slope.Comparisons indicate that the effect of friction angleφon the failure mechanism is significant,as compared to other parameters.Fora certain set ofβ andH/b,the highest curve of the stability factor corresponds to the one with φ =40°,as shown in Fig. 2.Thus,it is expected that the face failure mechanism associated with the high value ofφiscaused by a significant increase in the frictional component of shear resistance of soil slope,which constrains the weakest failure zone to be close to the top inclined surface of conical slope.

It is worth noting that a below-toe failure was not observed for all cases ofφ =10°-40°.The boundary size of the slope problem is also shown in Fig. 1,where the rigid boundary at the bottom is defined by 0.4Hbelow the toe of slope.All final adaptive meshes were examined after the analyses,where there was no intersection of localized shear band of the mesh at the bottom boundary.This confirms that the depth of slope foundation is not set to be too shallow,and is sufficient enough to simulate the rigid base of the slope.Hence,the above-toe failure mechanisms were simulated correctly for all cases ofφ =10°-40°.On the contrary,the belowtoe failure mechanisms occurred for the conical slopes in cohesive soils,as reported by Keawsawasvong and Ukritchon(2017a).Thus,these results suggest that a rotational failure mechanism,such as a log-spiral failure surface,passing through either the toe or the face of slope and taking into account a cylindrical shape of conical slope,would provide an accurate stability estimate of the problem with a cohesive-frictional material.For a purely cohesive soil,a failure mechanism passing below the toe should be adopted instead.In both types of materials,the critical rotational failure mechanism shall be obtained by optimizing the geometrical parameters in order to produce the least solution to the stability factor.

4.Design equation of drained stability of conical slopes and its applications

Fig. 3.Comparison of the final adaptive meshes for different conical slope angles(β),where H/b=8 and φ =20°.

Fig. 4.Comparison of the final adaptive meshes for different excavated height ratios(H/b),whereβ =75°and φ =20°.

Fig. 5.Comparison of the final adaptive meshes for different friction angles(φ),where H/b=6 and β =60°.

A trial-and-error method of curve fitting was performed to determine the mathematical equation that best matched with the computed data from FELA.The LB solutions were chosen in the curve fitting method,since they provide a safe estimate of the stability factor of a conical slope.In this study,it was found that the exponential function of tanφprovides a best fit of the computed LB data.Hence,based on a nonlinear regression to the computed LB solutions,a new design equation that provides an approximate LB solution to the stability factor of conical slope in cohesive-frictional soils is proposed:

whereNcis the stability factor of conical slope in purely cohesive soils(i.e.φ =0°)(Keawsawasvong and Ukritchon,2017a),andaandmare the constant coefficients.

Nc,aandmare the functions of the slope height ratio and slope inclination angle.The least square method(Walpole et al.,2002;Sauer,2014)was employed to determine the optimal solution to these coefficients,as summarized in Table 1.In this nonlinear regression,the LB solutions to a conical slope withφ=0°by Keawsawasvong and Ukritchon(2017a)were utilized and merged with those of the present study with φ=10°-40°to obtain accurate coefficients in the range ofφ=0°-40°.Note that the approximate LB solution ofNcwas numerically derived by Keawsawasvong and Ukritchon(2017a)for the conical slope in clay with constant and linearly increasing strengths with depth.The accuracy of the approximate solution ofNin Eq.(2)is shown in Fig. 6,where it was compared with the compute dLB solutions.A good agreement in the factorNbetween these two solutions was observed,where the coefficient of determination(R2=99.97%)was very high.

The application of the proposed design equation was demonstrated by solving a conical slope in a cohesive-frictional soil as:γ =20 kN/m3,ci=20 kPa,φi=25°,H=5 m andb=1 m.For this stability evaluation,the calculations were done to determine the factor of safety(FS)against slope failure by varying slope inclinationangles(β)of 60°-90°.The conventional definition ofFSin slope stability is employed in this calculation as

Table 1Constant coefficients for the approximate lower bound solution to stability factor of conical slopes in homogeneous cohesive-frictional soils.

Fig. 6.Comparison of the stability factors between the proposed design equation and the computed lower bound solutions.

Note that a distinction must be made between the input strength parameters(ci,φi)and mobilized strength parameters(c,φ).Since the proposed design solution in Eq.(2)corresponds to the approximate LB solution in the limit state condition,the shear strength parameters of the soil used in the equation correspond to the mobilized parameters.To determine theFS,an iterative calculation is necessary,becauseFSis not the direct solution obtained from Eq.(2).A small code was written in EXCEL to perform the iterativeFScalculation.First,aFSwas assumed(FSk)and the mobilized strength parameters were calculated from Eq.(3).Next,using these mobilized strength parameters,the proposed equation(Eq.(2))was employed to compute the mobilized stability factor(γH/c),i.e.k.The resulting factor of safety(FSr)was calculated asFSr=kci/(γH).This process is repeatedly performed until the assumed and resulting factors converge to the same one,i.e.FSr≈FSk.

Fig. 7 shows a comparison of theFSbetween the prediction and the FELA for the selected example.The calculations of the former employ the proposed equation(Eq.(2))while those of the latter also use the strength reduction method(Krabbenhoft and Lyamin,2015),which is comparable to the conventionalFSdefinition in Eq.(3).Good agreement in theFSvalues between Eq.(2)and FELA was observed.Consequently,this result confirmed that the proposed equation(Eq.(2))was very accurate for the safety analysis of conical slopes in cohesive-frictional soils in practice.

Fig. 7.Comparisons of FS between the proposed design equation and FELA.

It is generally recognized that the assumption of the associated fl ow rule is not realistic in practice for drained stability of cohesive frictional soils or effective stress analysis,since this assumption gives rise to an excessive volume change at the limit state for soils with high friction angles.The consequences of the associated fl ow rule with stability calculations by limit analysis have been extensively discussed(Sloan,2013).The concept of Davis(1968)was adopted in this paper by following a recommendation suggested by Sloan(2013),and this concept provides a practical means for dealing with the non-associated f l ow rule in limit analysis.In the method of Davis(1968),the “reduced”strength parameters are computed as the “input”strength parameters for the limit analysis as follows:

Note that the dilation angle of soil generally varies during the plastic deformation under shearing and converges to zero at the critical state.However,the typical range of dilation angles(ψ)from 0 toφ/3 may be assumed for a constant rate of dilation when experimental data are not available,as suggested by Sloan(2013).The iterative code for theFSdetermination of a conical slope using Eq.(2)can be easily adapted to take into account the nonassociated f l ow rule,as proposed by Davis(1968).Basically,for a given dilation angle of soil,the “reduced”strength parameters are calculated according to Eq.(4)and used as “input”parameters for the determination of theFSagainst slope failure.Then,the iterative calculation of theFSas described earlier can be straightforwardly performed on the “reduced”input strength parameters.Fig. 7 also shows a comparison of theFSvalues for the same conical slope example with the non-associated f l ow rule with ψ =0°.Note that the FELA solutions for the non-associated f l ow rule were determined using “reduced”strength parameters as input ones for the calculations of the strength reduction method(Krabbenhoft and Lyamin,2015).For ψ =0°,there is also a good agreement ofFSbetween Eq.(2)and FELA solutions.The influence of the dilation angle on the predictedFScan be clearly observed in Fig. 7.For this conical slope,the non-associated f l ow rule results in a decreasedFSagainst slope failure by 9.4%forψ =0°,as compared to the associated f l ow rule(ψ = φ).

5.Conclusions

Drained stability of conical slopes in homogeneous cohesive frictional soils was investigated by FELA in an axisymmetric condition.Three parameters were considered,i.e.conical height ratio(H/b),conical slope angle(β),and soil friction angle(φ).Two failure mechanisms,i.e.toe and face failures,were observed and associated with these parameters.A new design equation employing the exponential function of tanφwas proposed for an approximate LB solution of stability evaluation of conical slopes in a general soil with both undrained and drained strengths in practice.Applications of numerical examples demonstrated that the proposed design equation was accurately enough for a practical stability evaluation of conical slopes in both cohesive and cohesive-frictional soils with both associated and non-associated fl ow rules.

Conflict of interest

The authors wish to confirm that there are no known Conflicts of interest associated with this publication and there has been no significant financial support for this work that could have influenced its outcome.

Britto AM,Kusakabe O.Stability of unsupported axisymmetric excavations in soft clay.Géotechnique 1982;32(3):261-70.

Chakraborty M,Kumar J.Solving axisymmetric stability problems by using upper bound finite elements limit analysis and linear optimization.Journal of Engineering Mechanics 2014;140(6).https://doi.org/10.1061/(ASCE)EM.1943-7889.0000726.

Chen X,Huang J.Stability analysis of bank slope under conditions of reservoir impounding and rapid drawdown.Journal of Rock Mechanics and Geotechnical Engineering 2011;3(Supp.1):429-37.

Davis EH.Theories of plasticity,failure of soil masses.In:Lee IK,editor.Soil mechanics selected topics.New York,USA:Elsevier;1968.

Dincer I.Thermal energy storage:systems and applications.2nd ed.Wiley;2010.

Ghanbari E,Hamidi A.Improvement parameters in dynamic compaction adjacent to the slopes.Journal of Rock Mechanics and Geotechnical Engineering 2015;7(2):233-6.

Grif fi ths DV,Koutsabeloulis N.Finite element analysis of vertical excavations.Computers and Geotechnics 1985;1(3):221-35.

Grif fi ths D,Marquez R.Three-dimensional slope stability analysis by elasto-plastic finite elements.Géotechnique 2007;57(6):537-46.

Grif fi ths DV,Yu X.Another look at the stability of slopes with linearly increasing undrained strength.Géotechnique 2015;65(10):824-30.

Keawsawasvong S,Ukritchon B.Finite element limit analysis of pullout capacity of planar caissons in clay.Computers and Geotechnics 2016a;75:12-7.

Keawsawasvong S,Ukritchon B.Ultimate lateral capacity of two dimensional plane strain rectangular pile in clay.Geomechanics and Engineering 2016b;11(2):235-52.

Keawsawasvong S,Ukritchon B.Stability of unsupported conical excavations in non-homogeneous clays.Computers and Geotechnics 2017a;81:125-36.

Keawsawasvong S,Ukritchon B.Undrained lateral capacity of I-shaped concrete piles.Songklanakarin Journal of Science and Technology 2017b;39(6):751-8.

Keawsawasvong S,Ukritchon B.Undrained stability of an active planar trapdoor in non-homogeneous clays with a linear increase of strength with depth.Computers and Geotechnics 2017c;81:284-93.

Keawsawasvong S,Ukritchon B.Undrained limiting pressure behind soil gaps in contiguous pile walls.Computers and Geotechnics 2017d;83:152-8.

Keawsawasvong S,Ukritchon B.Finite element analysis of undrained stability of cantilever fl ood walls.International Journal of Geotechnical Engineering 2017e;11(4):355-67.

Kelesoglu MK.The evaluation of three-dimensional effects on slope stability by the strength reduction method.KSCE Journal of Civil Engineering 2016;20(1):229-42.

Khatri VN,Kumar J.Stability of an unsupported vertical circular excavation in clays under undrained condition.Computers and Geotechnics 2010;37(3):419-24.

Krabbenhoft K,Lyamin AV.Strength reduction finite-element limit analysis.Géotechnique Letters 2015;5(4):250-3.

Krabbenhoft K,Lyamin AV,Sloan SW.SONIC:a conic programming code for plasticity problems.Version 1.0.Newcastle,Australia:University of Newcastle;2006.

Krabbenhoft K,Lyamin AV,Sloan SW.Formulation and solution of some plasticity problems as conic programs.International Journal of Solids and Structures 2007a;44(5):1533-49.

Krabbenhoft K,Lyamin AV,Sloan SW,Wriggers P.An interior-point method for elastoplasticity.International Journal for Numerical Methods in Engineering 2007b;69:592-626.

Krabbenhoft K,Lyamin AV,Krabbenhoft J.Optum computational engineering(OptumG2).2015.http://www.optumce.com.

Kumar J,Chakraborty D.Stability numbers for an unsupported vertical circular excavation inc-φsoil.Computers and Geotechnics 2012;39:79-84.

Kumar J,Chakraborty M.Upper bound axisymmetric limit analysis using Mohr-Coulomb yield criterion finite elements and linear optimization.Journal of Engineering Mechanics 2014;140(12).https://doi.org/10.1061/(ASCE)EM.1943-7889.0000820.

Kumar J,Chakraborty M,Sahoo JP.Stability of unsupported vertical circular excavations.Journal of Geotechnical and Geoenvironmental Engineering 2014;140(7).https://doi.org/10.1061/(ASCE)GT.1943-5606.0001118.

Lee KS.Underground thermal energy storage:green energy and technology.Springer;2013.

Li X.Finite element analysis of slope stability using a nonlinear failure criterion.Computers and Geotechnics 2006;34(2):127-36.

Li AJ,Merifield RS,Lyamin AV.Limit analysis solutions for three dimensional undrained slopes.Computers and Geotechnics 2009;36(8):1330-51.

Li AJ,Merifield RS,Lyamin AV.Three-dimensional stability charts for slopes based on limit analysis methods.Canadian Geotechnical Journal 2010;47(12):1316-34.

Li E,Zhuang X,Zheng W,Cai Y.Effect of graph generation on slope stability analysis based on graph theory.Journal of Rock Mechanics and Geotechnical Engineering 2014;6(4):380-6.

Lim K,Li AJ,Lyamin AV.Three-dimensional slope stability assessment of twolayered undrained clay.Computers and Geotechnics 2015a;70:1-7.

Lim K,Lyamin AV,Cassidy MJ,Li AJ.Three-dimensional slope stability charts for frictional fill materials placed on purely cohesive clay.International Journal of Geomechanics 2015b;16(2).https://doi.org/10.1061/(ASCE)GM.1943-5622.0000 526.

Lim K,Li AJ,Schmid A,Lyamin AV.Slope-stability assessments using finite-element limit-analysis methods.International Journal of Geomechanics 2016;17(2).https://doi.org/10.1061/(ASCE)GM.1943-5622.0000715.

Lu Y.Deformation and failure mechanism of slope in three dimensions.Journal of Rock Mechanics and Geotechnical Engineering 2015;7(2):109-19.

Makrodimopoulos A,Martin CM.Lower bound limit analysis of cohesive-frictional materials using second-order cone programming.International Journal for Numerical Methods in Engineering 2006;66(4):604-34.

Makrodimopoulos A,Martin CM.Upper bound limit analysis using simplex strain elements and second-order cone programming.International Journal for Analytical Methods in Geomechanics 2007;31(6):835-65.

Martin CM.The use of adaptive finite element limit analysis to reveal slip-line fields.Géotechnique Letters 2011;1(2):23-9.

Michalowski RL.Stability charts for uniform slopes.Journal of Geotechnical and Geoenvironmental Engineering 2002;128(4):351-5.

Michalowski RL.Limit analysis and stability charts for 3D slope failures.Journal of Geotechnical and Geoenvironmental Engineering 2010;136(4):583-93.

Ochs F,Heidemann W,Müller-Steinhagen H.Performance of large-scale seasonal thermal energy stores.Journal of Solar Energy Engineering 2009;131(4).https://doi.org/10.1115/1.3197842.

Pastor J,Turgeman S.Limit analysis in axisymmetrical problems:numerical determination of complete statistical solutions.International Journal of Mechanical Sciences 1982;24(2):95-117.

Salomoni VA,Majorana CE,Giannuzzi GM,Miliozzi A.Thermal- fluid fl ow within innovative heat storage concrete systems for solar power plants.Journal of Engineering Mechanics 2008;18:969-99.

Sauer T.Numerical analysis.London,UK:Pearson Education Limited;2014.

Schmidt T,Miedaner O.Solar district heating guidelines.Intelligent Energy Europe 2012.

Sloan SW.Geotechnical stability analysis.Géotechnique 2013;63(7):531-72.

Tang C,Toh K,Phoon K.Axisymmetric lower-bound limit analysis using finite elements,second-order cone programming.Journal of Engineering Mechanics 2014;140:268-78.

Ukritchon B,Kea wsa wasvong S.Undrained pullout capacity of cylindrical suction caissons by finite element limit analysis.Computers and Geotechnics 2016;80:301-11.

Ukritchon B,Keawsawasvong S.Unsafe error in conventional shape factor for shallow circular foundations in normally consolidated clays.Journal of Geotechnical and Geoenvironmental Engineering 2017a;143(6).https://doi.org/10.1061/(ASCE)GT.1943-5606.0001670.

Ukritchon B,Keawsawasvong S.Error in Ito and Matsui’s limit equilibrium solution of lateral force on a row of stabilizing piles.Journal of Geotechnical and Geoenvironmental Engineering 2017b;143(9). https://doi.org/10.1061/(ASCE)GT.1943-5606.0001753.

Ukritchon B,Keawsawasvong S.Design equations for undrained stability of opening in underground walls.Tunnellingand Underground SpaceTechnology 2017c;70:214-20.

Ukritchon B,Keawsawasvong S.Undrained lateral capacity of rectangular piles under a general loading direction and full f l ow mechanism.KSCE Journal of Civil Engineering 2017.https://doi.org/10.1007/s12205-017-0062-7.

Walpole RE,Myers RH,Myers SL,Ye K.Probability and statistics for engineering and scientists.7th ed.Prentice Hall;2002.

Wright MH.The interior-point revolution in optimization:history,recent developments,and lasting consequences.Bulletin of the American Mathematical Society 2005;42:39-56.

Xie M,Wang Z,Liu X,Xu B.Three-dimensional critical slip surface locating and slope stability assessment for lava lobe of Unzen volcano.Journal of Rock Mechanics and Geotechnical Engineering 2011;3(1):82-9.

Yang S,Ren X,Zhang J.Study on embedded length of piles for slope reinforced with one row of piles.Journal of Rock Mechanics and Geotechnical Engineering 2011;3(2):167-78.

Yao A,Tian Z,Jin Y.Sliding surface searching method for slopes containing a potential weak structural surface.Journal of Rock Mechanics and Geotechnical Engineering 2014;6(3):275-9.

Yu HS,Salgado R,Sloan S,Kim J.Limit analysis versus limit equilibrium for slope stability.Journal of Geotechnical and Geoenviron mental Engineering 1998;124(1):1-11.

Zheng Y,Tang X,Zhao S,Deng C,Lei W.Strength reduction and step-loading finite element approaches in geotechnical engineering.Journal of Rock Mechanics and Geotechnical Engineering 2009;1(1):21-30.