- Email: [email protected]
An iterative method for the linearization of nonlinear failure criteria
ARTICLE IN PRESS
JID: ADES
[m5G;May 3, 2017;5:52]
Advances in Engineering Software 0 0 0 (2017) 1–6
Contents lists available at ScienceDirect
Advances in Engineering Software journal homepage: www.elsevier.com/locate/advengsoft
An iterative method for the linearization of nonlinear failure criteria S.K. Sharan Bharti School of Engineering, Laurentian University, Sudbury, Canada
a r t i c l e
i n f o
Article history: Received 16 February 2017 Accepted 24 April 2017 Available online xxx Keywords: Elastoplastic analysis Finite element method Nonlinear failure criteria Linearization Stress analysis Hoek–Brown failure criterion
a b s t r a c t Linearization of a nonlinear failure criterion may be required if a software incorporating the nonlinear failure criterion is not readily available or if an additional computational difficulty is encountered due to the nonlinearity of the failure criterion. In this paper, a novel method is proposed for the linearization of nonlinear failure criteria. The proposed method is based on an iterative procedure and the leastsquares regression is used for the linearization. The material is assumed to be elastic-perfectly plastic. By conducting analytical and finite element analyses of stresses and displacements around underground openings in rock mass governed by the generalized Hoek–Brown failure criterion, it is shown that the proposed method is convergent, efficient and effective. © 2017 Elsevier Ltd. All rights reserved.
1. Introduction For the elastoplastic finite element analysis of solids encountered in structural or geotechnical engineering, the material properties may be defined more accurately by using a nonlinear failure criterion [1–3] as shown in Fig. 1. In this figure, σ 1 and σ 3 are, respectively, the major and minor principal stresses and f is a nonlinear function of σ 3. The problem of non-convergence of non-elastic problems in implicit analysis is quite common. The nonlinearity of the failure criterion and the yield function may lead to additional problem of non-convergence. Therefore, an otherwise converging solution may become non-convergent due to the nonlinearity of the failure criterion [4–6]. Furthermore, due to this and other reasons, a suitable software incorporating the nonlinear failure criterion may not be readily available. For example, the conventional linear Mohr–Coulomb (M-C) failure criterion is the most commonly used failure criterion in geotechnical engineering and is implemented in most of the general purpose finite element software. However, this failure criterion is not suitable for many types of rocks, particularly a jointed rock mass. For such rocks, the nonlinear generalized Hoek–Brown (GH-B) failure criterion [2] is more appropriate. To the author’s knowledge, this failure criterion has not been implemented in any general purpose finite element software. Some special purpose finite element software, for example, Phase2 [6] have implemented this failure criterion. However, numerical problems such as nonconvergence of results may be encountered if the failure criterion is changed from the linear M-C to the nonlinear GH-B [4–6].
E-mail address: [email protected]
In order to circumvent these computational difficulties, researchers [2,4–5,7–17] have proposed various methods of linearizing the nonlinear failure criterion as shown in Fig. 2 in which g is a linear function of σ 3. Most of these methods are based on the curve-fitting approach within a predetermined range of stress. In the earlier methods [7–8], the linearization was assumed to be independent of the actual state of stress. Later, it was realized that the linearization should be based on the range of stresses existing within the yielded zone and thus improved methods of linearization were proposed [9–17] and applied to different types of problems. In all such methods, the range of major or minor principal stress within the yielded zone is either assumed or obtained by using available analytical solutions for simple cases. Thus the methods are valid for such specific cases only and for other problems, the results are highly approximate. Most recently, a general method was proposed [18] for the computation of the equivalent parameters for the linear M-C failure criterion for rock mass governed by the nonlinear GH-B failure criterion. The objective of this paper is to present the earlier work [18] in a general form so that it could be applied to a wide variety of other problems encountered in structural and geotechnical engineering. The paper has been revised further in the light of the comments and discussions following the presentation of the paper at the conference [18]. In the proposed method, the material behaviour is assumed to be elastic-perfectly plastic (Fig. 3) and the effect of strainhardening or softening behaviour is not considered. In this figure, ε1 is the major principal strain. The range of minor or major principal stresses within the yielded zone is determined by using an
http://dx.doi.org/10.1016/j.advengsoft.2017.04.008 0965-9978/© 2017 Elsevier Ltd. All rights reserved.
Please cite this article as: S.K. Sharan, An iterative method for the linearization of nonlinear failure criteria, Advances in Engineering Software (2017), http://dx.doi.org/10.1016/j.advengsoft.2017.04.008
ARTICLE IN PRESS
JID: ADES 2
[m5G;May 3, 2017;5:52]
S.K. Sharan / Advances in Engineering Software 000 (2017) 1–6
Fig. 1. A nonlinear failure criterion. Fig. 4. Linearization of a nonlinear failure criterion.
2. Yield and plastic potential functions It is assumed that the nonlinear failure criterion is independent of the intermediate principal stress and as mentioned earlier, it may be expressed in the following form:
σ1 = f ( σ3 )
(1)
The nonlinear function f depends on the mechanical properties of the material. Compressive stresses are considered to be positive. The yield function F is thus expressed as
F = σ1 − f ( σ3 ) = 0 Fig. 2. A linear failure criterion.
(1a)
The plastic potential function Q is assumed to be linear and non-associated and is given by
Q = σ1 − Kd σ3
(2)
where Kd is the dilation parameter of the material. The value of Kd = 1 corresponds to a nondilatant material. 3. The proposed method
Fig. 3. Elastoplastic model used for the analysis.
iterative procedure. The initial values of the range of minor principal stresses are computed by using an elastic analysis or it may be based on a reasonable estimate as proposed in [17]. The leastsquares regression is used for the linearization of the nonlinear yield function within the range computed by using the analytical solution or by using the finite element method. The linear yield function thus obtained is then used to conduct the elastoplastic analysis and the new range of values of the minor or major principal stresses is determined. The iterative process is continued until a convergence criterion is satisfied. Results of several numerical tests are presented to demonstrate the effectiveness of the proposed method.
The nonlinear function f has to be linearized within a particular range of values of σ 1 or σ 3 . As mentioned earlier, in the existing methods, the range is either assumed or determined analytically for simplified cases. For higher accuracy, the range of values of σ 1 or σ 3 should be the actual values existing within the yielded zone. However, in general, this is initially unknown. As mentioned earlier, the proposed method is based on an iterative procedure to determine the actual range. Initial values of the range are based on the elastic analysis. Alternatively, as mentioned earlier, the initial values may be based on approximations proposed in [17]. The linearization based on the range of σ 1 requires the solution of a nonlinear equation. Therefore, in the proposed method, it is preferred to use the range of σ 3 . For any given range of σ 3 from σ 3 min to σ 3max , the function f may be approximated by a linear function g as shown in Fig. 4 and the linearized failure criterion may be expressed as
σ1 = g ( σ3 ) = K σ3 + C
(3)
where K and C are linearization parameters for the material. The corresponding equivalent or virtual yield function is thus given by
Feq = σ1 − g(σ3 ) = 0
(3a)
As mentioned earlier, the proposed linearization is based on the least-squares regression and thus K and C are obtained by solving the following system of two linear equations:
Please cite this article as: S.K. Sharan, An iterative method for the linearization of nonlinear failure criteria, Advances in Engineering Software (2017), http://dx.doi.org/10.1016/j.advengsoft.2017.04.008
ARTICLE IN PRESS
JID: ADES
[m5G;May 3, 2017;5:52]
S.K. Sharan / Advances in Engineering Software 000 (2017) 1–6
3
Table 1 Properties of rock mass. Rock mass
E (GPa)
ν
σ ci (MPa)
mb
s
a
σ 0 (MPa)
A B C D E
8.944 1.414 30.115 12.989 0.867
0.25 0.30 0.25 0.25 0.30
80 20 51 30 15
2.0121 0.6567 6.6746 4.4695 0.6625
3.87E-03 4.19E-04 6.22E-02 2.05E-02 2.15E-04
0.506 0.522 0.501 0.502 0.533
10 1 10 10 1
Table 2 Cases analysed and number of iterations. Rock mass example
Case
Support pressure p0 (Mpa)
Number of iterations n
A
A1 A2 B1 B2 C1 C2 D1 D2 E1 E2
0 3.29E−01 0 4.36E−02 0 1.15E−01 0 3.70E−01 0 5.53E−02
3 3 5 4 5 4 4 3 6 4
B C D E
1 3
1 2 C σ3 max − σ32min σ33max − σ33min +
K
2
σ3 max
σ 3 f ( σ 3 )d σ 3
=
(4)
σ3 min
1 2
K
σ32max − σ32min + C (σ3 max − σ3 min ) =
σ3 max
f ( σ 3 )d σ 3
(5)
σ3 min
be almost identical to those obtained by using the analytical solution. The linearized function g(σ 3 ) is expressed in terms of the M-C failure criterion for which K and C are given by
K=
1 + sin φ 1 − sin φ
(7)
C=
2c cos φ 1 − sin φ
(8)
The integrations in the above equations may be performed analytically or numerically depending upon the function f(σ 3 ). By using the “new” values of K and C, the elastoplastic analysis is conducted by using the linear failure criterion and the “new” values of σ 3 min and σ 3max are computed. The iterative procedure is continued until a specified convergence criterion is satisfied. It is quite feasible to implement the proposed iterative procedure in a finite element program. However, for the numerical test example presented in this paper, a commercial finite element software was run separately for each iteration after obtaining the results for σ 3 min and σ 3max from the finite element analysis in the previous iteration and computing the “new” values of linearization parameters K and C by solving Eqs. (4)–(5), including the integrations involved, with the use of a commercial mathematical software.
where φ is the angle of internal friction and c is the cohesion of the material. Eqs. (7)–(8) may be written in the following form to determine the values of φ and c using the values of K and C computed by solving Eqs. (4)–(5):
4. Test examples
4.1. Analytical solutions
In order to study the effectiveness of the proposed method, plane strain elastoplastic analyses were conducted to compute stresses and displacements around underground openings in rock mass governed by the GH-B failure criterion [2] for which the function f(σ 3 ) is given by
Analytical solutions are available for the plane strain analysis of stresses and displacements around circular openings in rock mass subject to a hydrostatic stress field and governed by M-C [20] and GH-B [3] failure criteria. The radius of the opening was considered to be 5.0 m. Five different types of rock masses A-E were considered and their material properties are given in Table 1. These were taken from [11]. Rock mass cases analysed in [11] producing excessively large displacements were not considered because the analytical solutions are based on a small-displacement theory. In Table 1, E = Young’s modulus, ν = Poisson’s ratio and σ 0 = hydrostatic in situ stress. For each type of rock mass, both unsupported and supported openings were considered. Thus, as shown in Table 2, 10 cases (A1-E1 and A2-E2) were analysed. As in [11], the support pressure p0 was considered to be equal to 20% of the radial stress at the elastic-plastic interface and the values are shown in Table 2.
f (σ3 ) = σ3 + σci mb
a σ3 +s σci
(6)
where σ ci is the uniaxial compressive strength of the intact rock and mb , s and a are GH-B constants that depend on the properties of the rock mass. For this failure criterion, integrations involved in Eqs. (4)–(5) may be performed analytically [11]. However, for generality, the integrations were performed numerically by using the commercial mathematical software Maple [19] and the results were found to
φ = sin−1 c=
K − 1 K+1
C (1 − sin φ ) 2 cos φ
(9)
(10)
Several example cases were analysed by using analytical solutions for a circular opening in a hydrostatic stress field and the finite element method for a non-circular opening in a nonhydrostatic stress field. The rock mass was assumed to be nondilatant.
Please cite this article as: S.K. Sharan, An iterative method for the linearization of nonlinear failure criteria, Advances in Engineering Software (2017), http://dx.doi.org/10.1016/j.advengsoft.2017.04.008
ARTICLE IN PRESS
JID: ADES 4
[m5G;May 3, 2017;5:52]
S.K. Sharan / Advances in Engineering Software 000 (2017) 1–6 Table 3 Results of convergence test for case E1. Iteration number i
σ 3 min (Mpa)
σ 3max (Mpa)
φ (degrees)
c (Mpa)
R (m)
umax (mm)
Percentage change in umax (%)
1 2 3 4 5 6
0 0 0 0 0 0
1.0 0 0 0 0.2989 0.2668 0.2611 0.2600 0.2598
32.481 41.801 42.654 42.815 42.846 42.852
0.1945 0.0894 0.0834 0.0823 0.0821 0.0821
6.710 6.921 6.905 6.902 6.901 6.901
10.253 11.745 11.793 11.800 11.802 11.802
– 12.70 0.41 0.06 0.01 0.00
Table 4 Comparison of results for EM-C parameters. Case
φ Present (degrees)
Exact (degrees)
Percentage error (%)
Present (Mpa)
Exact (Mpa)
Percentage error (%)
A1 A2 B1 B2 C1 C2 D1 D2 E1 E2
50.89 48.88 46.60 44.22 57.45 56.86 49.14 46.83 42.85 40.04
50.60 48.71 46.20 44.07 57.38 56.80 48.74 46.67 42.39 39.89
0.6 0.4 0.9 0.3 0.1 0.1 0.8 0.3 1.1 0.4
1.095 1.323 0.099 0.126 1.883 1.945 1.070 1.329 0.082 0.114
1.116 1.338 0.102 0.128 1.888 1.950 1.098 1.343 0.085 0.116
−1.9 −1.1 −2.8 −1.1 −0.3 −0.3 −2.6 −1.0 −3.7 −1.2
c
Table 5 Comparison of results for the radius of elastic-plastic interface and maximum displacement. Case
A1 A2 B1 B2 C1 C2 D1 D2 E1 E2
R
umax
Present (m)
Exact (m)
Percentage error (%)
Present (mm)
Exact (mm)
Percentage error (%)
5.78 5.53 6.22 5.79 5.18 5.14 5.93 5.62 6.90 6.16
5.83 5.55 6.30 5.82 5.19 5.14 5.99 5.64 7.12 6.19
−0.8 −0.3 −1.2 −0.4 −0.1 0.0 −1.0 −0.4 −3.1 −0.6
8.36 7.40 6.07 5.05 2.12 2.08 5.98 5.17 11.80 8.80
8.41 7.41 6.16 5.05 2.12 2.08 6.03 5.17 12.41 8.81
−0.6 −0.1 −1.5 −0.0 0.0 0.0 −0.9 −0.0 −4.9 −0.2
Table 6 Results of convergence test for the finite element analysis of a D-shaped opening. Iteration number i
σ 3 min (Mpa)
σ 3max (Mpa)
φ (degrees)
C (Mpa)
umax (mm)
Percentage change in umax (%)
Elastic 1 2 3
… 0 0 0
… 0.50 0.44 0.45
… 37.86 38.85 38.68
… 0.124 0.114 0.116
0.0325 0.0392 0.0397 0.0396
… 20.7 1.1 −0.1
The iterative procedure was continued until the following convergence criterion was satisfied for the absolute value of the maximum total displacement umax :
i −1 (umax − uimax )/uimax × 100 < τ
(11)
where the superscript i stands for the iteration number and τ is the tolerance limit. The results presented for this analysis are for τ = 0.01%. As shown in Table 2, the number of iterations required for this accuracy varied from 3 to 6. As expected, the number of iterations required for supported openings was relatively less than that for the unsupported opening. Detailed results of convergence test for a typical case E1 are presented in Table 3. The values of the equivalent Mohr-Coulomb (EM-C) parameters φ and c were computed by using Eqs. (9)–(10) and a comparison of results for these parameters obtained by using the proposed method and analytical solutions [11] is shown in Table 4. The maximum error in φ was 1.1% and that in c was 3.7%. A comparison of results for the radius of the yielded zone R and the maximum total displacement umax computed by using the
proposed method and analytical solutions [3] are shown in Table 5. The maximum error in R was found to be 3.1% and that in umax was 4.9%. Since the linearization of a nonlinear function is only approximate, these errors may be considered to be insignificant. 4.2. Finite element solutions The proposed method was applied to the analysis of a noncircular opening in a non- hydrostatic stress field by using the finite element method. The properties of rock mass was considered to be those for rock mass E shown in Table 1. The in situ stresses σ h0 and σ v0 in the horizontal and vertical directions, respectively, were considered to be 1 MPa and 0.5 MPa. The opening was considered to be d-shaped and unsupported. The commercially available finite element software Phase2 [6] was used for the analysis. Eight-node isoparametric planestrain elements were used for the discretization. By using symmetry, only one half of the domain was analysed. Since a high ac-
Please cite this article as: S.K. Sharan, An iterative method for the linearization of nonlinear failure criteria, Advances in Engineering Software (2017), http://dx.doi.org/10.1016/j.advengsoft.2017.04.008
JID: ADES
ARTICLE IN PRESS
[m5G;May 3, 2017;5:52]
S.K. Sharan / Advances in Engineering Software 000 (2017) 1–6
5
Fig. 5. The finite element model for one-half of rock mass around a d-shaped opening.
Fig. 7. Comparison of results for the yielded zone and strength factors.
performed for the finite element analysis are presented in Table 6. As shown in the table, three iterations were required to achieve an accuracy of τ = 0.1%. The results thus obtained were compared with those obtained directly by using the GH-B failure criterion. Fig. 7 shows a comparison of results for the yielded Gauss points of elements (marked as ‘X’ ). The figure also shows strength factors [6] at some sampling points. The maximum error was found to be less than 7%. The comparison of results for displacements in the near-field is shown in Fig. 8 in terms of contour lines and values at some sampling points. The maximum error in displacement was found to be less than 0.01 m. Considering the facts that the analysis is based on finite element discretization and representing a nonlinear failure criterion by a linear one involves approximations, the errors are insignificant. For all the comparisons, the agreement in results was thus very good. 5. Conclusions
Fig. 6. Dimensions of the d-shaped opening and the near-field discretization.
curacy was not needed, infinite elements were used at the truncation boundary of the finite element model as shown in Fig. 5. The dimensions of the opening and the near-field discretization are shown in Fig. 6. The analysis was conducted by using the proposed method using the M-C failure criterion. Detailed results of convergence test
An iterative method was proposed for the linearization of nonlinear failure criteria. The principal advantage of the proposed method is that it is quite general and a general purpose finite element software that does not have the nonlinear failure criterion as an option for material properties may be used to analyze stresses and displacements in materials governed by such failure criterion. The method is also suitable for situations when a computational difficulty, such as divergence or non-convergence, is encountered in the direct use of a nonlinear failure criterion. The effectiveness of the method was demonstrated by analyzing several example cases related to the analysis of underground openings in rock mass governed by the generalized Hoek–Brown failure criterion. The method may be used to analyze a wide variety of problems in structural and geotechnical engineering involving materials governed by novel nonlinear yield functions.
Please cite this article as: S.K. Sharan, An iterative method for the linearization of nonlinear failure criteria, Advances in Engineering Software (2017), http://dx.doi.org/10.1016/j.advengsoft.2017.04.008
ARTICLE IN PRESS
JID: ADES 6
[m5G;May 3, 2017;5:52]
S.K. Sharan / Advances in Engineering Software 000 (2017) 1–6
Fig. 8. Comparison of results for displacements (m).
Acknowledgement The financial support provided by Laurentian University (Grant No. LURF S16 41-1-6090547) is gratefully acknowledged. References [1] Hoek E, Brown ET. Empirical strength criterion for rock masses. ASCE J Geotech Eng 1980;106(GT9):1013–36. [2] Hoek E, Carranza-Torres C, Corkum B. Hoek–Brown failure criterion - 2002 edition. In: Proceedings of the North American rock mechanics symposium; 2002. p. 267–73. [3] Sharan SK. Analytical solutions for stresses and displacements around a circular opening in a generalized Hoek–Brown rock. Int J Rock Mech Min Sci 2008;45:78–85.
[4] Naznin R. Equivalent Mohr–Coulomb dilation parameter for Hoek–Brown rock M.A.Sc. thesis. Sudbury: Laurentian University; 2007. [5] Sharan SK, Naznin R. Equivalent angle of dilation for the analysis of underground openings in rock mass obeying Hoek–Brown plasticity. In: Proceedings of the twenty-third Canadian congress of applied mechanics; 2011. p. 140–3. [6] Rocscience Inc. Phase2: finite element analysis and support design for excavations. Toronto: Rockscience Inc.; 2015. [7] Hoek E. Estimating Mohr–Coulomb friction and cohesion values from the Hoek–Brown failure criterion. Int J Rock Mech Min Sci Geomech Abstr 1990;27(3):227–9. [8] Hoek E, Brown ET. Practical estimates of rock mass strength. Int J Rock Mech Min Sci Geomech Abstr 1997;34(8):1165–86. [9] Sofianos AI, Halakatevakis N. Equivalent tunneling Mohr–Coulomb strength parameters for given Hoek–Brown ones. Int J Rock Mech Mining Sci 2002;39:131–7. [10] Sofianos AI. Tunneling Mohr–Coulomb strength parameters for rock masses satisfying the generalized Hoek–Brown criterion. Int J Rock Mech Mining Sci 2003;40:435–40. [11] Sofianos AI, Nomikos PP. Equivalent Mohr–coulomb and generalized Hoek–Brown strength parameters for supported axisymmetric tunnels in plastic or brittle rock. Int J Rock Mech Mining Sci 2006;43:683–704. [12] Yang XL, Yin JH. Linear Mohr–Coulomb strength parameters from the non-linear Hoek–Brown rock masses. Int J Non-Linear Mech 20 06;41:10 0 0–5. [13] Jimenez R, Serrano A, Olalla C. Linearization of Hoek and Brown rock failure criterion for tunnelling in elasto-plastic rock masses. Int J Rock Mech Mining Sci 2008;45:1153–63. [14] Yang XL, Yin JH. Slope equivalent Mohr–Coulomb strength parameters for rock masses satisfying the Hoek–Brown criterion. Rock Mech Rock Eng 2010;43:505–11. [15] Shen J, Priest SD, Karakus M. Determination of Mohr-Coulomb shear strength parameters from generalized Hoek–Brown criterion for slope stability analysis. Rock Mech Rock Eng 2012;45:123–9. [16] Sharan SK, Naznin R. Linearization of the Hoek–Brown failure criterion for non-hydrostatic stress fields. In: Topping BHV, editor. Proceedings of the eighth international conference on engineering computational technology, Stirlingshire, UK. Civil-Comp Press; 2012 Paper 111. [17] Nomikos PP, Sofianos AI. Response of some lined tunnels within original HB and equivalent MC rock masses. In: Proceedings of the ISRM-sponsored international symposium on rock mechanics: rock characterisation, modelling and engineering design methods, SINOROCK; 2009. [18] Sharan SK. Computation of the equivalent Mohr–Coulomb parameters for Hoek–Brown rock. In: Kruis J, Tsompanakis Y, Topping BHV, editors. Proceedings of the fifteenth international conference on civil, structural and environmental engineering computing, Stirlingshire, UK. Civil-Comp Press; 2015. Paper 254. doi: 10.4203/ccp.108.254. [19] http://www.maplesoft.com. [20] Park KH, Kim YJ. Analytical solution for a circular opening in an elastic-brittle-plastic rock. Int J Rock Mech Min Sci 2006;43:616–22.
Please cite this article as: S.K. Sharan, An iterative method for the linearization of nonlinear failure criteria, Advances in Engineering Software (2017), http://dx.doi.org/10.1016/j.advengsoft.2017.04.008
JID: ADES
[m5G;May 3, 2017;5:52]
Advances in Engineering Software 0 0 0 (2017) 1–6
Contents lists available at ScienceDirect
Advances in Engineering Software journal homepage: www.elsevier.com/locate/advengsoft
An iterative method for the linearization of nonlinear failure criteria S.K. Sharan Bharti School of Engineering, Laurentian University, Sudbury, Canada
a r t i c l e
i n f o
Article history: Received 16 February 2017 Accepted 24 April 2017 Available online xxx Keywords: Elastoplastic analysis Finite element method Nonlinear failure criteria Linearization Stress analysis Hoek–Brown failure criterion
a b s t r a c t Linearization of a nonlinear failure criterion may be required if a software incorporating the nonlinear failure criterion is not readily available or if an additional computational difficulty is encountered due to the nonlinearity of the failure criterion. In this paper, a novel method is proposed for the linearization of nonlinear failure criteria. The proposed method is based on an iterative procedure and the leastsquares regression is used for the linearization. The material is assumed to be elastic-perfectly plastic. By conducting analytical and finite element analyses of stresses and displacements around underground openings in rock mass governed by the generalized Hoek–Brown failure criterion, it is shown that the proposed method is convergent, efficient and effective. © 2017 Elsevier Ltd. All rights reserved.
1. Introduction For the elastoplastic finite element analysis of solids encountered in structural or geotechnical engineering, the material properties may be defined more accurately by using a nonlinear failure criterion [1–3] as shown in Fig. 1. In this figure, σ 1 and σ 3 are, respectively, the major and minor principal stresses and f is a nonlinear function of σ 3. The problem of non-convergence of non-elastic problems in implicit analysis is quite common. The nonlinearity of the failure criterion and the yield function may lead to additional problem of non-convergence. Therefore, an otherwise converging solution may become non-convergent due to the nonlinearity of the failure criterion [4–6]. Furthermore, due to this and other reasons, a suitable software incorporating the nonlinear failure criterion may not be readily available. For example, the conventional linear Mohr–Coulomb (M-C) failure criterion is the most commonly used failure criterion in geotechnical engineering and is implemented in most of the general purpose finite element software. However, this failure criterion is not suitable for many types of rocks, particularly a jointed rock mass. For such rocks, the nonlinear generalized Hoek–Brown (GH-B) failure criterion [2] is more appropriate. To the author’s knowledge, this failure criterion has not been implemented in any general purpose finite element software. Some special purpose finite element software, for example, Phase2 [6] have implemented this failure criterion. However, numerical problems such as nonconvergence of results may be encountered if the failure criterion is changed from the linear M-C to the nonlinear GH-B [4–6].
E-mail address: [email protected]
In order to circumvent these computational difficulties, researchers [2,4–5,7–17] have proposed various methods of linearizing the nonlinear failure criterion as shown in Fig. 2 in which g is a linear function of σ 3. Most of these methods are based on the curve-fitting approach within a predetermined range of stress. In the earlier methods [7–8], the linearization was assumed to be independent of the actual state of stress. Later, it was realized that the linearization should be based on the range of stresses existing within the yielded zone and thus improved methods of linearization were proposed [9–17] and applied to different types of problems. In all such methods, the range of major or minor principal stress within the yielded zone is either assumed or obtained by using available analytical solutions for simple cases. Thus the methods are valid for such specific cases only and for other problems, the results are highly approximate. Most recently, a general method was proposed [18] for the computation of the equivalent parameters for the linear M-C failure criterion for rock mass governed by the nonlinear GH-B failure criterion. The objective of this paper is to present the earlier work [18] in a general form so that it could be applied to a wide variety of other problems encountered in structural and geotechnical engineering. The paper has been revised further in the light of the comments and discussions following the presentation of the paper at the conference [18]. In the proposed method, the material behaviour is assumed to be elastic-perfectly plastic (Fig. 3) and the effect of strainhardening or softening behaviour is not considered. In this figure, ε1 is the major principal strain. The range of minor or major principal stresses within the yielded zone is determined by using an
http://dx.doi.org/10.1016/j.advengsoft.2017.04.008 0965-9978/© 2017 Elsevier Ltd. All rights reserved.
Please cite this article as: S.K. Sharan, An iterative method for the linearization of nonlinear failure criteria, Advances in Engineering Software (2017), http://dx.doi.org/10.1016/j.advengsoft.2017.04.008
ARTICLE IN PRESS
JID: ADES 2
[m5G;May 3, 2017;5:52]
S.K. Sharan / Advances in Engineering Software 000 (2017) 1–6
Fig. 1. A nonlinear failure criterion. Fig. 4. Linearization of a nonlinear failure criterion.
2. Yield and plastic potential functions It is assumed that the nonlinear failure criterion is independent of the intermediate principal stress and as mentioned earlier, it may be expressed in the following form:
σ1 = f ( σ3 )
(1)
The nonlinear function f depends on the mechanical properties of the material. Compressive stresses are considered to be positive. The yield function F is thus expressed as
F = σ1 − f ( σ3 ) = 0 Fig. 2. A linear failure criterion.
(1a)
The plastic potential function Q is assumed to be linear and non-associated and is given by
Q = σ1 − Kd σ3
(2)
where Kd is the dilation parameter of the material. The value of Kd = 1 corresponds to a nondilatant material. 3. The proposed method
Fig. 3. Elastoplastic model used for the analysis.
iterative procedure. The initial values of the range of minor principal stresses are computed by using an elastic analysis or it may be based on a reasonable estimate as proposed in [17]. The leastsquares regression is used for the linearization of the nonlinear yield function within the range computed by using the analytical solution or by using the finite element method. The linear yield function thus obtained is then used to conduct the elastoplastic analysis and the new range of values of the minor or major principal stresses is determined. The iterative process is continued until a convergence criterion is satisfied. Results of several numerical tests are presented to demonstrate the effectiveness of the proposed method.
The nonlinear function f has to be linearized within a particular range of values of σ 1 or σ 3 . As mentioned earlier, in the existing methods, the range is either assumed or determined analytically for simplified cases. For higher accuracy, the range of values of σ 1 or σ 3 should be the actual values existing within the yielded zone. However, in general, this is initially unknown. As mentioned earlier, the proposed method is based on an iterative procedure to determine the actual range. Initial values of the range are based on the elastic analysis. Alternatively, as mentioned earlier, the initial values may be based on approximations proposed in [17]. The linearization based on the range of σ 1 requires the solution of a nonlinear equation. Therefore, in the proposed method, it is preferred to use the range of σ 3 . For any given range of σ 3 from σ 3 min to σ 3max , the function f may be approximated by a linear function g as shown in Fig. 4 and the linearized failure criterion may be expressed as
σ1 = g ( σ3 ) = K σ3 + C
(3)
where K and C are linearization parameters for the material. The corresponding equivalent or virtual yield function is thus given by
Feq = σ1 − g(σ3 ) = 0
(3a)
As mentioned earlier, the proposed linearization is based on the least-squares regression and thus K and C are obtained by solving the following system of two linear equations:
Please cite this article as: S.K. Sharan, An iterative method for the linearization of nonlinear failure criteria, Advances in Engineering Software (2017), http://dx.doi.org/10.1016/j.advengsoft.2017.04.008
ARTICLE IN PRESS
JID: ADES
[m5G;May 3, 2017;5:52]
S.K. Sharan / Advances in Engineering Software 000 (2017) 1–6
3
Table 1 Properties of rock mass. Rock mass
E (GPa)
ν
σ ci (MPa)
mb
s
a
σ 0 (MPa)
A B C D E
8.944 1.414 30.115 12.989 0.867
0.25 0.30 0.25 0.25 0.30
80 20 51 30 15
2.0121 0.6567 6.6746 4.4695 0.6625
3.87E-03 4.19E-04 6.22E-02 2.05E-02 2.15E-04
0.506 0.522 0.501 0.502 0.533
10 1 10 10 1
Table 2 Cases analysed and number of iterations. Rock mass example
Case
Support pressure p0 (Mpa)
Number of iterations n
A
A1 A2 B1 B2 C1 C2 D1 D2 E1 E2
0 3.29E−01 0 4.36E−02 0 1.15E−01 0 3.70E−01 0 5.53E−02
3 3 5 4 5 4 4 3 6 4
B C D E
1 3
1 2 C σ3 max − σ32min σ33max − σ33min +
K
2
σ3 max
σ 3 f ( σ 3 )d σ 3
=
(4)
σ3 min
1 2
K
σ32max − σ32min + C (σ3 max − σ3 min ) =
σ3 max
f ( σ 3 )d σ 3
(5)
σ3 min
be almost identical to those obtained by using the analytical solution. The linearized function g(σ 3 ) is expressed in terms of the M-C failure criterion for which K and C are given by
K=
1 + sin φ 1 − sin φ
(7)
C=
2c cos φ 1 − sin φ
(8)
The integrations in the above equations may be performed analytically or numerically depending upon the function f(σ 3 ). By using the “new” values of K and C, the elastoplastic analysis is conducted by using the linear failure criterion and the “new” values of σ 3 min and σ 3max are computed. The iterative procedure is continued until a specified convergence criterion is satisfied. It is quite feasible to implement the proposed iterative procedure in a finite element program. However, for the numerical test example presented in this paper, a commercial finite element software was run separately for each iteration after obtaining the results for σ 3 min and σ 3max from the finite element analysis in the previous iteration and computing the “new” values of linearization parameters K and C by solving Eqs. (4)–(5), including the integrations involved, with the use of a commercial mathematical software.
where φ is the angle of internal friction and c is the cohesion of the material. Eqs. (7)–(8) may be written in the following form to determine the values of φ and c using the values of K and C computed by solving Eqs. (4)–(5):
4. Test examples
4.1. Analytical solutions
In order to study the effectiveness of the proposed method, plane strain elastoplastic analyses were conducted to compute stresses and displacements around underground openings in rock mass governed by the GH-B failure criterion [2] for which the function f(σ 3 ) is given by
Analytical solutions are available for the plane strain analysis of stresses and displacements around circular openings in rock mass subject to a hydrostatic stress field and governed by M-C [20] and GH-B [3] failure criteria. The radius of the opening was considered to be 5.0 m. Five different types of rock masses A-E were considered and their material properties are given in Table 1. These were taken from [11]. Rock mass cases analysed in [11] producing excessively large displacements were not considered because the analytical solutions are based on a small-displacement theory. In Table 1, E = Young’s modulus, ν = Poisson’s ratio and σ 0 = hydrostatic in situ stress. For each type of rock mass, both unsupported and supported openings were considered. Thus, as shown in Table 2, 10 cases (A1-E1 and A2-E2) were analysed. As in [11], the support pressure p0 was considered to be equal to 20% of the radial stress at the elastic-plastic interface and the values are shown in Table 2.
f (σ3 ) = σ3 + σci mb
a σ3 +s σci
(6)
where σ ci is the uniaxial compressive strength of the intact rock and mb , s and a are GH-B constants that depend on the properties of the rock mass. For this failure criterion, integrations involved in Eqs. (4)–(5) may be performed analytically [11]. However, for generality, the integrations were performed numerically by using the commercial mathematical software Maple [19] and the results were found to
φ = sin−1 c=
K − 1 K+1
C (1 − sin φ ) 2 cos φ
(9)
(10)
Several example cases were analysed by using analytical solutions for a circular opening in a hydrostatic stress field and the finite element method for a non-circular opening in a nonhydrostatic stress field. The rock mass was assumed to be nondilatant.
Please cite this article as: S.K. Sharan, An iterative method for the linearization of nonlinear failure criteria, Advances in Engineering Software (2017), http://dx.doi.org/10.1016/j.advengsoft.2017.04.008
ARTICLE IN PRESS
JID: ADES 4
[m5G;May 3, 2017;5:52]
S.K. Sharan / Advances in Engineering Software 000 (2017) 1–6 Table 3 Results of convergence test for case E1. Iteration number i
σ 3 min (Mpa)
σ 3max (Mpa)
φ (degrees)
c (Mpa)
R (m)
umax (mm)
Percentage change in umax (%)
1 2 3 4 5 6
0 0 0 0 0 0
1.0 0 0 0 0.2989 0.2668 0.2611 0.2600 0.2598
32.481 41.801 42.654 42.815 42.846 42.852
0.1945 0.0894 0.0834 0.0823 0.0821 0.0821
6.710 6.921 6.905 6.902 6.901 6.901
10.253 11.745 11.793 11.800 11.802 11.802
– 12.70 0.41 0.06 0.01 0.00
Table 4 Comparison of results for EM-C parameters. Case
φ Present (degrees)
Exact (degrees)
Percentage error (%)
Present (Mpa)
Exact (Mpa)
Percentage error (%)
A1 A2 B1 B2 C1 C2 D1 D2 E1 E2
50.89 48.88 46.60 44.22 57.45 56.86 49.14 46.83 42.85 40.04
50.60 48.71 46.20 44.07 57.38 56.80 48.74 46.67 42.39 39.89
0.6 0.4 0.9 0.3 0.1 0.1 0.8 0.3 1.1 0.4
1.095 1.323 0.099 0.126 1.883 1.945 1.070 1.329 0.082 0.114
1.116 1.338 0.102 0.128 1.888 1.950 1.098 1.343 0.085 0.116
−1.9 −1.1 −2.8 −1.1 −0.3 −0.3 −2.6 −1.0 −3.7 −1.2
c
Table 5 Comparison of results for the radius of elastic-plastic interface and maximum displacement. Case
A1 A2 B1 B2 C1 C2 D1 D2 E1 E2
R
umax
Present (m)
Exact (m)
Percentage error (%)
Present (mm)
Exact (mm)
Percentage error (%)
5.78 5.53 6.22 5.79 5.18 5.14 5.93 5.62 6.90 6.16
5.83 5.55 6.30 5.82 5.19 5.14 5.99 5.64 7.12 6.19
−0.8 −0.3 −1.2 −0.4 −0.1 0.0 −1.0 −0.4 −3.1 −0.6
8.36 7.40 6.07 5.05 2.12 2.08 5.98 5.17 11.80 8.80
8.41 7.41 6.16 5.05 2.12 2.08 6.03 5.17 12.41 8.81
−0.6 −0.1 −1.5 −0.0 0.0 0.0 −0.9 −0.0 −4.9 −0.2
Table 6 Results of convergence test for the finite element analysis of a D-shaped opening. Iteration number i
σ 3 min (Mpa)
σ 3max (Mpa)
φ (degrees)
C (Mpa)
umax (mm)
Percentage change in umax (%)
Elastic 1 2 3
… 0 0 0
… 0.50 0.44 0.45
… 37.86 38.85 38.68
… 0.124 0.114 0.116
0.0325 0.0392 0.0397 0.0396
… 20.7 1.1 −0.1
The iterative procedure was continued until the following convergence criterion was satisfied for the absolute value of the maximum total displacement umax :
i −1 (umax − uimax )/uimax × 100 < τ
(11)
where the superscript i stands for the iteration number and τ is the tolerance limit. The results presented for this analysis are for τ = 0.01%. As shown in Table 2, the number of iterations required for this accuracy varied from 3 to 6. As expected, the number of iterations required for supported openings was relatively less than that for the unsupported opening. Detailed results of convergence test for a typical case E1 are presented in Table 3. The values of the equivalent Mohr-Coulomb (EM-C) parameters φ and c were computed by using Eqs. (9)–(10) and a comparison of results for these parameters obtained by using the proposed method and analytical solutions [11] is shown in Table 4. The maximum error in φ was 1.1% and that in c was 3.7%. A comparison of results for the radius of the yielded zone R and the maximum total displacement umax computed by using the
proposed method and analytical solutions [3] are shown in Table 5. The maximum error in R was found to be 3.1% and that in umax was 4.9%. Since the linearization of a nonlinear function is only approximate, these errors may be considered to be insignificant. 4.2. Finite element solutions The proposed method was applied to the analysis of a noncircular opening in a non- hydrostatic stress field by using the finite element method. The properties of rock mass was considered to be those for rock mass E shown in Table 1. The in situ stresses σ h0 and σ v0 in the horizontal and vertical directions, respectively, were considered to be 1 MPa and 0.5 MPa. The opening was considered to be d-shaped and unsupported. The commercially available finite element software Phase2 [6] was used for the analysis. Eight-node isoparametric planestrain elements were used for the discretization. By using symmetry, only one half of the domain was analysed. Since a high ac-
Please cite this article as: S.K. Sharan, An iterative method for the linearization of nonlinear failure criteria, Advances in Engineering Software (2017), http://dx.doi.org/10.1016/j.advengsoft.2017.04.008
JID: ADES
ARTICLE IN PRESS
[m5G;May 3, 2017;5:52]
S.K. Sharan / Advances in Engineering Software 000 (2017) 1–6
5
Fig. 5. The finite element model for one-half of rock mass around a d-shaped opening.
Fig. 7. Comparison of results for the yielded zone and strength factors.
performed for the finite element analysis are presented in Table 6. As shown in the table, three iterations were required to achieve an accuracy of τ = 0.1%. The results thus obtained were compared with those obtained directly by using the GH-B failure criterion. Fig. 7 shows a comparison of results for the yielded Gauss points of elements (marked as ‘X’ ). The figure also shows strength factors [6] at some sampling points. The maximum error was found to be less than 7%. The comparison of results for displacements in the near-field is shown in Fig. 8 in terms of contour lines and values at some sampling points. The maximum error in displacement was found to be less than 0.01 m. Considering the facts that the analysis is based on finite element discretization and representing a nonlinear failure criterion by a linear one involves approximations, the errors are insignificant. For all the comparisons, the agreement in results was thus very good. 5. Conclusions
Fig. 6. Dimensions of the d-shaped opening and the near-field discretization.
curacy was not needed, infinite elements were used at the truncation boundary of the finite element model as shown in Fig. 5. The dimensions of the opening and the near-field discretization are shown in Fig. 6. The analysis was conducted by using the proposed method using the M-C failure criterion. Detailed results of convergence test
An iterative method was proposed for the linearization of nonlinear failure criteria. The principal advantage of the proposed method is that it is quite general and a general purpose finite element software that does not have the nonlinear failure criterion as an option for material properties may be used to analyze stresses and displacements in materials governed by such failure criterion. The method is also suitable for situations when a computational difficulty, such as divergence or non-convergence, is encountered in the direct use of a nonlinear failure criterion. The effectiveness of the method was demonstrated by analyzing several example cases related to the analysis of underground openings in rock mass governed by the generalized Hoek–Brown failure criterion. The method may be used to analyze a wide variety of problems in structural and geotechnical engineering involving materials governed by novel nonlinear yield functions.
Please cite this article as: S.K. Sharan, An iterative method for the linearization of nonlinear failure criteria, Advances in Engineering Software (2017), http://dx.doi.org/10.1016/j.advengsoft.2017.04.008
ARTICLE IN PRESS
JID: ADES 6
[m5G;May 3, 2017;5:52]
S.K. Sharan / Advances in Engineering Software 000 (2017) 1–6
Fig. 8. Comparison of results for displacements (m).
Acknowledgement The financial support provided by Laurentian University (Grant No. LURF S16 41-1-6090547) is gratefully acknowledged. References [1] Hoek E, Brown ET. Empirical strength criterion for rock masses. ASCE J Geotech Eng 1980;106(GT9):1013–36. [2] Hoek E, Carranza-Torres C, Corkum B. Hoek–Brown failure criterion - 2002 edition. In: Proceedings of the North American rock mechanics symposium; 2002. p. 267–73. [3] Sharan SK. Analytical solutions for stresses and displacements around a circular opening in a generalized Hoek–Brown rock. Int J Rock Mech Min Sci 2008;45:78–85.
[4] Naznin R. Equivalent Mohr–Coulomb dilation parameter for Hoek–Brown rock M.A.Sc. thesis. Sudbury: Laurentian University; 2007. [5] Sharan SK, Naznin R. Equivalent angle of dilation for the analysis of underground openings in rock mass obeying Hoek–Brown plasticity. In: Proceedings of the twenty-third Canadian congress of applied mechanics; 2011. p. 140–3. [6] Rocscience Inc. Phase2: finite element analysis and support design for excavations. Toronto: Rockscience Inc.; 2015. [7] Hoek E. Estimating Mohr–Coulomb friction and cohesion values from the Hoek–Brown failure criterion. Int J Rock Mech Min Sci Geomech Abstr 1990;27(3):227–9. [8] Hoek E, Brown ET. Practical estimates of rock mass strength. Int J Rock Mech Min Sci Geomech Abstr 1997;34(8):1165–86. [9] Sofianos AI, Halakatevakis N. Equivalent tunneling Mohr–Coulomb strength parameters for given Hoek–Brown ones. Int J Rock Mech Mining Sci 2002;39:131–7. [10] Sofianos AI. Tunneling Mohr–Coulomb strength parameters for rock masses satisfying the generalized Hoek–Brown criterion. Int J Rock Mech Mining Sci 2003;40:435–40. [11] Sofianos AI, Nomikos PP. Equivalent Mohr–coulomb and generalized Hoek–Brown strength parameters for supported axisymmetric tunnels in plastic or brittle rock. Int J Rock Mech Mining Sci 2006;43:683–704. [12] Yang XL, Yin JH. Linear Mohr–Coulomb strength parameters from the non-linear Hoek–Brown rock masses. Int J Non-Linear Mech 20 06;41:10 0 0–5. [13] Jimenez R, Serrano A, Olalla C. Linearization of Hoek and Brown rock failure criterion for tunnelling in elasto-plastic rock masses. Int J Rock Mech Mining Sci 2008;45:1153–63. [14] Yang XL, Yin JH. Slope equivalent Mohr–Coulomb strength parameters for rock masses satisfying the Hoek–Brown criterion. Rock Mech Rock Eng 2010;43:505–11. [15] Shen J, Priest SD, Karakus M. Determination of Mohr-Coulomb shear strength parameters from generalized Hoek–Brown criterion for slope stability analysis. Rock Mech Rock Eng 2012;45:123–9. [16] Sharan SK, Naznin R. Linearization of the Hoek–Brown failure criterion for non-hydrostatic stress fields. In: Topping BHV, editor. Proceedings of the eighth international conference on engineering computational technology, Stirlingshire, UK. Civil-Comp Press; 2012 Paper 111. [17] Nomikos PP, Sofianos AI. Response of some lined tunnels within original HB and equivalent MC rock masses. In: Proceedings of the ISRM-sponsored international symposium on rock mechanics: rock characterisation, modelling and engineering design methods, SINOROCK; 2009. [18] Sharan SK. Computation of the equivalent Mohr–Coulomb parameters for Hoek–Brown rock. In: Kruis J, Tsompanakis Y, Topping BHV, editors. Proceedings of the fifteenth international conference on civil, structural and environmental engineering computing, Stirlingshire, UK. Civil-Comp Press; 2015. Paper 254. doi: 10.4203/ccp.108.254. [19] http://www.maplesoft.com. [20] Park KH, Kim YJ. Analytical solution for a circular opening in an elastic-brittle-plastic rock. Int J Rock Mech Min Sci 2006;43:616–22.
Please cite this article as: S.K. Sharan, An iterative method for the linearization of nonlinear failure criteria, Advances in Engineering Software (2017), http://dx.doi.org/10.1016/j.advengsoft.2017.04.008