《1 前言》

1 前言

自20世纪20年代以来,岩土工程的极限分析方法(主要指滑移线场法、上下限分析法与极限平衡法)获得蓬勃发展,并广泛应用于工程实际。这些方法有的需要作一些人为假设,有的求解范围十分有限,限制了这种方法的发展与应用。而有限元法数值方法适应性强,应用范围宽,但无法求出工程设计中十分有用的稳定安全系数F与极限承载力,从而制约了有限元数值分析方法在岩土工程中的应用。

1975年,英国科学家 Zienkiewicz提出在有限元中采用增加荷载或降低岩土强度的方法来计算岩土工程的极限荷载和安全系数[1]。20世纪80年代、90年代曾用于边坡和地基的稳定分析[2],但是由于当时缺少严格可靠、功能强大的大型有限元程序以及强度准则的选用和具体操作技术掌握不够等原因,导致计算精度不足,而没有得到岩土工程界的广泛采纳。

20世纪末前后,国际上又发表了多篇文章[3~6],研究了有限元强度折减法求解均质土坡的稳定安全系数F,由于一些算例得到的结果与传统方法求解结果比较接近,逐渐得到学术界认可,有些国外学者认为有限元强度折减法使边坡稳定分析进入了一个新的时代。尤其是1999年美国科罗拉多矿业学院的D.V.Griffith等人用自编有限元程序对均质土坡进行稳定分析[3~5],与其他程序不同之处是该程序能够模拟水位和孔隙水压力的影响,还可进行库水下降情况下边坡的稳定分析。

1997年宋二祥介绍和研究了有限元强度折减法在土坡中的应用。21世纪初,国内学者开始致力于有限元强度折减法在边坡稳定分析中应用的研究[7~15],文献[8~10]是国内较早的研究文章。首先进行了该法基本理论和提高计算精度的研究,

随着计算精度的提高,这种方法受到国内岩土工程界和设计部门的广泛关注。一方面扩大了有限元极限分析法的应用范围,另一方面也开始被一些工程设计部门实际采用。目前,有限元极限分析法正进入方兴未艾的发展阶段。

国际上采用有限元强度折减法求解边坡的滑面与安全系数,用有限元增量加载(超载)法求地基的极限承载力,前者研究较多,并取得了可喜的成果,而后者还研究不多。将上述两种方法统称为有限元极限分析法,因为它们本质上都是采用数值分析手段求解的极限分析法。

国际上采用自编数值分析程序居多,其应用范围限于二维平面土坡与土基的分析。而国内趋向于采用国际大型通用程序,不仅计算方便,而且程序可靠,功能强大,计算精度高,表述清晰并便于工程应用。同时,将该方法的应用范围大为扩大,从均质的土坡、土基扩大到具有结构面的岩坡与岩基;从二维扩大到三维;还扩展到寻找边(滑)坡中多个潜在滑面;进行岩土与结构共同作用的支挡结构设计;用于计算机仿真地基承载板载荷试验,尤其是首次扩展到求隧道的稳定安全系数。更可喜的是有些工程设计部门已经采用这一方法进行边坡稳定分析与支挡结构设计,为有限元极限分析法开拓了灿烂的应用前景。

近来,我国在有限元极限分析法方面的发展极为迅速,国内许多学者作了有效的工作。目前,有限元极限分析法的发展刚处于起步阶段,离达到改革岩土工程设计方法的目的还有很大距离。所以,希望能与国内同仁们一道,为发展这一学科分支与革新岩土工程设计方法作出贡献。

《2 有限元极限分析法的原理》

2 有限元极限分析法的原理

《2.1 有限元极限分析法中安全系数的定义》

2.1 有限元极限分析法中安全系数的定义

有限元极限分析法中安全系数的定义依据岩土工程出现破坏状态的原因不同而不同。一类如边(滑)坡工程多数由于岩土受环境影响,岩土强度降低,导致边(滑)坡失稳破坏。这类工程宜采用强度贮备安全系数(也称强度安全系数),即可通过不断降低岩土强度使有限元计算最终达到破坏状态为止。强度降低的倍数就是强度贮备安全系数,因而这种有限元极限分析法称为有限元强度折减法。其实,无论何种岩土工程,凡是采用有限元强度折减法来求安全系数的都是强度贮备安全系数。

另一类,如地基工程由于地基上荷载不断增大而导致地基失稳破坏,这类工程采用荷载增大的倍数作为超载安全系数,称为有限元增量加载(超载)法。同样,凡是采用这一方法求解的安全系数的都是超载安全系数。显然,上述两种方法求得的安全系数是不同的,即不同的安全系数定义得到的安全系数是不同的,采用同一安全系数算出的支挡结构上的推力也是不同的。

《2.2 有限元极限分析法原理》

2.2 有限元极限分析法原理[11,12]

1)有限元强度折减法[13]    对于岩土中广泛采用的莫尔-库仑材料,强度折减安全系数ω可表示为

\(\begin{array}{c} \tau=(c+\sigma \tan \varphi) / \omega=c^{\prime}+\sigma \tan \varphi^{\prime} \\ c^{\prime}=c / \omega, \tan \varphi^{\prime}=(\tan \varphi) / \omega_{0} \end{array}\)       (1)

有限元计算中不断降低边坡中岩土抗剪强度直至达到破坏状态为止。程序根据有限元计算结果自动得到破坏滑动面,并获得强度贮备安全系数。
然而,岩土材料有2个强度指标c与tanφ,却采用一个强度贮备安全系数,这意味着两个指标按同一比例下降,而实际岩土并非这样。这是强度贮备安全系数的不足。

2)有限元增量加载法[14]     在实际工程中,岩土破坏往往是一个渐进的破坏过程,岩土体是由初始的线弹性状态逐渐过渡到塑性流动,直至达到极限破坏状态。采用增量加载的方式求解地基的极限承载力就是这一思路下的产物。随着荷载的逐步增加,岩土体由弹性逐渐过渡到塑性,最后达到极限破坏状态,对应的荷载就为所要求的极限荷载。这方法称为有限元增量加载法或有限元超载法。

《2.3 有限元极限分析法的优越性》

2.3 有限元极限分析法的优越性

有限元极限分析法具有数值方法与经典极限分析法两者的优点,既具有数值方法适应性广的优点,又具有极限分析法贴近岩土工程设计,实用性强的优点。

1)用有限元强度折减法求解边坡安全系数时,不需要假定滑面的形状和位置,也无需进行条分,而是由程序自动求出滑面与强度贮备安全系数。

2)用有限元超载法求解地基极限承载力时,不必假定破坏面位置并给出理论解答,而由程序自动给出破坏机构与极限承载力。

3)具有数值分析法的各种优点,能够对复杂地貌、地质条件的各种岩土工程进行计算,不受工程的几何形状、边界条件以及材料不均匀等的限制。

4)能考虑应力一应变关系,提供应力、应变、位移和塑性区等力和变形的全部信息。

5)能够考虑岩土体与支护结构的共同作用,模拟施工开挖过程和渐进破坏过程[15]

有限元极限分析法可以利用国际上通用程序的强大功能,把计算结果准确、清晰地表达出来,实用、方便,必将导致岩土工程设计方法的重大改革,因而是一种颇有前途的计算方法。

《3 基本理论》

3 基本理论

《3.1 关于有限元极限分析法中岩土工程(边坡、地基、隧道)整体失稳的判据》

3.1 关于有限元极限分析法中岩土工程(边坡、地基、隧道)整体失稳的判据

有限元极限分析法中,无论是采用强度折减法还是超载法都需知道岩土工程整体失稳的判据[16]

岩土体的整体失稳破坏是指岩土体沿滑面(破裂面)发生滑落或坍塌,整个滑面达到极限平衡状态,并且土坡整体不能继续承载;同时,滑面上的应变与位移发生突变,岩土体沿滑面快速滑动直至滑落,坍塌。然而,人们至今仍对岩土体的整体失稳破坏没有统一的认识,一般认为边坡整个滑面上都达到极限平衡状态,就是整体失稳破坏,因而建议把滑面上塑性区贯通作为整体失稳的判据。不过,也有一些人认为,即使滑面上每点都达到极限应力状态,但由于边界条件的约束,土体没有足够的位移时仍不会发生滑动破坏。按此观点,把滑面上每点都达到极限平衡作为整体破坏条件不够全面,滑面上塑性区贯通只是破坏的必要条件,而非充分条件,它表征着渐进破坏的开始。认为只有整个滑面上每点的应变也都达到极限应变才会发生滑动。显然,这一观点符合整体失稳破坏的实际情况。边坡失稳,滑体由稳定静止状态变为运动状态,滑面节点位移和塑性应变将产生突变,此时位移和塑性应变将以高速无限发展,直到滑体滑出。这一现象符合边坡破坏的概念,可把滑面上节点塑性应变或位移突变作为边坡整体失稳的标志。与此同时,笔者也发现在上述情况下,静力平衡有限元计算也正好表现出计算不收敛,因此也可将有限元静力计算是否收敛作为边坡失稳的判据。这也表明目前国际上惯用的以计算机不收敛作为破坏判据是合适的。当然,这一判据不适用于由于计算失误而引起的计算机不收敛。

图1a为有节理岩石边坡达到整体破坏状态后产生的直线滑动破坏形式,可见破坏后边坡由稳定状态转变为运动状态,滑体产生很大的位移,而且无限发展。图1b为边坡滑动面上单元节点水平位移(坡顶UX1、坡中UX2、坡脚UX3)随着荷载的逐步增加而逐渐增大的曲线走势图。由图1可见,随着荷载的逐渐增加,当达到破坏状态后,节点的水平位移产生了突变。如有限元程序继续迭代下去,该节点的水平位移和塑性应变还将继续无限发展下去。但有限元程序已无法从有限元方程组中找到一个既能满足静力平衡又能满足应力一应变关系和强度准则的解,此时不管是从力的收敛标准,还是从位移的收敛标准来判断,有限元计算都不收敛。

《图1》

图1 边坡失稳后的特征

Fig.1 The failure phenomenon of plane sliding rock slope

《3.2 本构关系与屈服准则(强度准则)的选取》

3.2 本构关系与屈服准则(强度准则)的选取

3.2.1 本构关系与莫尔-库仑准则[9,17~20]  

 有限元极限分析法一般选用理想弹塑性模型,因为岩土工程的稳定问题都是力和强度问题,而不是位移问题,因而对本构关系的选择不十分严格,可选用最简单的理想弹塑性模型,不考虑岩土的硬化与软化。但对强度准则的选取则有严格的要求,以前该法计算精度不高,多数是由于强度准则选取不当所致。岩土材料常用的准则是莫尔-库仑准则,其表达式为

\(\sigma_{1}(1+\sin \varphi)-\sigma_{3}(1-\sin \varphi)=2 c \cos \varphi \) 或

\(\quad\left(I_{1} \sin \varphi\right) / 3+\left(\cos \theta_{\sigma}-\right. \\ \left.3^{-1 / 2} \sin \theta_{\sigma} \sin \varphi_{2}\right) J_{2}^{1 / 2}=c \cos \varphi \)       (2)

莫尔-库仑准则在π平面上的图形为不等角六边形,如图2所示,存在尖顶给数值计算带来困难。所以计算程序中采用莫尔-库仑准则时常要作一些近似处理,或采用与莫尔-库仑准则相应的广义米赛斯准则。 

《图2》

图2 二维应力空间的莫尔-库仑屈服条件

Fig.2 The Mohr-Coulomb yield criterions in 2D stress space

3.2.2 广义米赛斯准则(Drucker-Prager)广义米赛斯准则是在米赛斯准则的基础上,考虑平均压应力而将米赛斯条件推广成为如下形式

\(aI_{1}+J_{2}^{1 / 2}=k \\ I_{1} = \sigma_{1}+\sigma_{2} + \sigma_{3}, \\ J_{2}=\left[\left(\sigma_{1}-\sigma_{2}\right)^{2}+\left(\sigma_{1}-\sigma_{3}\right)^{2}+\left(\sigma_{2}-\sigma_{3}\right)^{2}\right] / 6 \)            (3)

式中I1,J2分别为应力张量的第一不变量和应力偏张量的第二不变量,a,k是与岩土材料内摩擦角φ和粘聚力c有关的系数。不同的α,k在π平面上代表不同的圆(图3),各准则的a,k见表1。

《图3》

图3 各屈服准则在π平面上的曲线

Fig.3 The yield surface on the deviator plane

式(3)是1952年Drucker-Prager提出的,因此广义米赛斯准则也称为Drucker-Prager(D-P)准则。广义米赛斯准则在主应力空间的屈服面为一圆
锥面,在π平面上为圆形,不存在尖顶处的数值计算问题,因此目前国际上流行的大型有限元软件ANSYS以及美国MSC公司的MARC,NASTRAN等 均采用了广义米赛斯准则。

《表1》

表1 各准则α,k参数表

Table 1 The relationship of different yield criterions

编号

准则种类

a

k

DP1

外角点外接圆

(2sinφ)/

(6ccos φ)/

3(3-sinø)

3(3-sinφ)

DP2

内角点外接圆

(2sin φ)/

(6cco8 φ)/

3(3+sinφ)

3(3+sinφ)

DP3

莫尔-库仑等面

23 sin φ/(2·

63 ccoe φ/(2·

积圆

3x(9-si㎡φ))12

3π(9-si㎡φ))12

DP4

平面应变关联法

(sin φ)/

(3ccos φ)/

则下莫尔-库仑

 

匹配圆

3(3+sinφ)12

3(3+si㎡φ)12

 

平面应变非关联

 

 

DP5

法则下莫尔-库

(sin φ)/3

ccos φ

仑匹配圆

3.2.3 莫尔-库仑准则的几种特殊情况 当内摩擦角φ=0时,式(2)可写成

\(J_{2}^{1 / 2} \cos \theta_{\sigma}-C=0\)               (4)

此即屈瑞斯卡条件,它在π平面上的图形是外接米赛斯圆的六边形。

如式(4)\(\theta_{\sigma}\)再等于零,即得米赛斯条件

\(J_{2}^{1 / 2}-C=0\)   (5)

注意此C与米赛斯原式中C不同,差二次方,因含义不同。

\(\theta_{\sigma}\)为常数时,屈服函数不再与\(\theta_{\sigma}\)或应力偏量第三不变量J3有关。它在π平面上为一个圆,这时式(2)可写成

\(\alpha I_{1}+J_{2}^{1 / 2}-k=0\)     (6)

这就是广义米赛斯条件。

在式(2)中取不同的\(\theta_{\sigma}\)值,即有不同的α,k值,由此可得到大小不同的圆锥形屈服面。当取\(\theta_{\sigma}\)=π/6时可得

\(a=(2 \sin \varphi) / 3^{1 / 2}(3-\sin \varphi) \\ k=(6 c \cos \varphi) / 3^{1 / 2}(3-\sin \varphi)\)     (7)

它在π平面上的屈服曲线是通过莫尔-库仑不等角六角形外角点的外接圆。

当取\(\theta_{\sigma}\)=-π/6时可得

\(\alpha=(2 \sin \varphi) / 3^{1 / 2}(3+\sin \varphi) \\ k=(6 \cos \varphi) / 3^{1 / 2}(3+\sin \varphi) \),  (8)

它在π平面上的屈服曲线是通过莫尔-库仑不等角

六角形内角点的外接圆。

将式(2)对\(\theta_{\sigma}\)微分,并使之等于零,这时F取极小,可得

\(\operatorname{tg} \theta_{\sigma}=-(\sin \varphi) / 3^{1 / 2}\)               (9)

取此\(\theta_{\sigma}\)值时,可得:

\(\alpha=(\sin \varphi) / 3^{1 / 2}\left(3+\sin ^{2} \varphi\right)^{1 / 2} \\ k=(3 \cos \varphi) / 3^{1 / 2}\left(3+\sin ^{2} \varphi\right)^{1 / 2} \)    (10)

式(10)为平面应变条件下采用关联流动法则时的 Drucker-Prager准则,称为平面应变关联流动法则下莫尔-库仑匹配准则,在π平面上的屈服曲线是通过莫尔-库仑不等角六角形的内切圆。

当取\(\theta_{\sigma}\)=0时,可得

\(\alpha=(\sin \varphi) / 3, k=c \cos \varphi\)        (11)

此即为平面应变非关联法则下莫尔-库仑匹配准则,它就是平面应变下真实的莫尔-库仑准则。在π平面上的屈服曲线为一稍大于内切圆的圆。公式的推导过程参见文献[21]。在平面应变条件下,可得0。=0,而0。=0时,能得到式(11)。

徐干成、郑颖人(1990)提出的莫尔-库仑等面积圆DP3准则在π平面上的面积等于不等角六边形莫尔-库仑准则的面积,图3可见,莫尔-库仑准则构成的六角形面积可用正弦定理求得(17)

\(S_{\text {morl }}=3 r_{1} r_{2} \sin \pi / 3=3 \sqrt{3} r_{1} r_{2} / 2\),

式中r1,r2为外角圆与内角圆的半径。对于半径为r的圆锥面积为:S=πr2,令S=Smorl可得

\(\begin{array}{c} r=\left(3 \sqrt{3} r r_{1} r_{3} / 2 \pi\right)^{1 / 2}= \\ 6^{1 / 2}\left(6 c \cos \varphi-2 I_{1} \sin \varphi\right) /\left(2 \sqrt{3} \pi\left(9-\sin ^{2} \varphi\right)\right)^{1 / 2} \\ r / 2^{1 / 2}=J_{2}^{1 / 2}= \\ \quad 6 \sqrt{3} c \cos \varphi /\left(2 \sqrt{3} \pi\left(9-\sin ^{2} \varphi\right)\right)^{1 / 2}- \\ 2 \times 3^{1 / 2} I_{1} \sin \varphi /\left(2 \times 3^{1 / 2} \pi\left(9-\sin ^{2} \varphi\right)\right)^{1 / 2} \circ \end{array}\)

由此可得莫尔-库仑等面积圆屈服准则的表达式

\(\begin{aligned} \alpha &=2 \sqrt{3} \sin \varphi /\left(2 \sqrt{3} \pi\left(9-\sin ^{2} \varphi\right)\right)^{1 / 2} \\ k &=6 \sqrt{3} \cos \varphi /\left(2 \sqrt{3} \pi\left(9-\sin ^{2} \varphi\right)\right)^{1 / 2} \end{aligned}\)          (12)

3.2.4 有限元极限分析法中强度准则的选用     莫尔-库仑准则应用最为广泛,但也存在诸多缺点,例如莫尔-库仑准则没有考虑中间主应力的影响,它在三维应力空间中不是一个连续函数,它在主应力空间中由6个屈服函数构成。莫尔-库仑准则的屈服面为不规则的六角形截面的角锥体表面,在π平面上的图形为不等角六边形,存在尖顶和菱角,给数值计算带来困难。

广义米赛斯准则(Drucker-Prager准则)在主应力空间的屈服面为一圆锥,在π平面上为圆形,不存在尖顶处的数值计算问题,因此目前国际上许多流行的大型有限元软件均采用了广义米赛斯准则。美国大型有限元软件ANSYS采用的是莫尔-库仑不等角六边形外角点外接圆DP1屈服准则。研究表明,采用该准则与传统莫尔-库仑屈服准则的计算结果有较大误差,不管是评价边坡稳定性,还是地基极限承载力等等,在实际工程中如果采用该准则是偏于不安全的。依据理论分析和计算实例,对屈服准则的选用提出如下建议:

1)对于平面应变问题,可采用与传统Mohr-Coulomb准则相匹配的DP4与DP5准则,它们有很高计算精度,其计算误差一般在1%~2%。
当采用DP5准则时,应使用非关联法则,此时,理论上应取膨胀角ψ=φ/2,实际应用中也可在0~φ/2取值。当采用DP4准则时,应采用关联流动法则,膨胀角取ψ=φ。

2)对于三维空间问题,可采用莫尔-库仑等面积圆DP3准则,也可获得较好的计算结果。

3.2.5 不同D-P准则条件下安全系数的转换(22)求解岩土工程安全系数一般采用有限元强度折减法。因而,对于D-P准则也采用c/ω,tanφ/ω的安全系数定义。

D-P准则中a,k有多种表达形式,采用不同的屈服条件得到的边坡稳定安全系数是不同的,但这些屈服条件是可以互相转换的。目前国际上的通用程序最多也只有外角点外接圆、内角点外接圆、内切圆3种D-P准则,因而实施屈服条件的转换是十分必要的。

设c0,φ0为初始强度参数,在外角点外接圆屈服准则条件下的安全系数为ω1,在莫尔-库仑等面积圆屈服准则条件下的安全系数为ω2,经过变换可以得到

\(\begin{aligned} \omega_{2} &=\left\{\left[3 \sqrt { 3 } \left(3\left(\cos ^{2} \varphi_{0} \omega_{1}^{2}+\sin ^{2} \varphi_{0}\right)^{1 / 2}-\right.\right.\right.\\ &\left.\left.\left.\sin \varphi_{0}\right)^{2}-8 \sin \varphi_{0}\right] / 18 \pi \cos ^{2} \varphi_{0}\right\}^{1 / 2} \\ \end{aligned}\)    (13)

式(13)即为外角点外接圆DP1屈服准则和莫尔-库仑等面积圆DP3准则之间的安全系数转换关系式。只要求得了外角点外接圆屈服准则条件下的安全系数ω1,利用该表达式就可以直接计算出莫尔-库仑等面积圆准则条件下的安全系数ω2。表2为不同参数条件下两种准则之间安全系数的实际转换数据。

《表2》

表2 不同参数条件下两种D-P准则之间的安全系数转换数据示例

Table 2    Safety factor by different D-P yield criterions

等面积圆DP3准则的

 

 

 

外角

点外接圆DP

1准则安全

系数w

 

 

 

安全系数w2

1

1.1

1.2

1.3

1.4

1.5

1.6

1.7

1.8

1.9

0

0.909

1.000

1.091

1.182

1.273

1.364

1.455

1.546

1.637

1.728

10

0.854

0.945

1.036

1.127

1.218

1.310

1.401

1.492

1.583

1.674

内摩擦      15

0.822

0.914

1.006

1.097

1.188

1.280

1.371

1.462

1.553

1.644

角φo/(°)   20

0.786

0.879

0.971

1.063

1.155

1.247

1.339

1.430

1.521

1.613

25

0.742

0.837

0.931

1.024

1.117

1.210

1.302

1.394

1.486

1.578

30

0.685

0.784

0.881

0.977

1.072

1.166

1.259

1.352

1.445

1.537

采用同样的方法可以得到外角点外接圆DP1屈服准则(非关联流动法则)和平面应变莫尔-库仑匹配DP5准则(非关联流动法则)之间的安全系数转换关系式。设c0,φ0为初始强度参数,在外角点外接圆屈服准则(非关联流动法则)条件下的安全系数为ω1,在平面应变莫尔-库仑匹配准则DP5(非关联流动法则)条件下的安全系数为ω2,经过变换得到

\(\begin{aligned} \omega_{2} &=\left\{\left[\left(3\left(\cos ^{2} \varphi_{0} \omega_{1}{ }^{2}+\sin ^{2} \varphi_{0}\right)^{1 / 2}-\sin \varphi_{0}\right)^{2}-\right.\right.\\ \left.\left.12 \sin ^{2} \varphi_{0}\right] / 12 \cos ^{2} \varphi_{0}\right\}^{1 / 2} \end{aligned}\)           (14)

这样,只要求得外接圆屈服准则条件下的安全系数ω1,利用该表达式就可直接计算出平面应变莫尔-库仑匹配DP5准则条件下的安全系数ω2

3.3 提高有限元极限分析法计算精度的条件[23]

为了达到计算精度,一般需满足如下条件:

1)要有一个成熟可靠、功能强的有限元程序,尤其是选用国际上公认的通用程序,这些程序安全可靠、功能与通用性强。

2)有可供实用的岩土本构模型和强度准则。笔者已经较好地解决了这个问题。

3)计算范围、边界条件、网格划分等要满足有限元计算精度要求。这对有经验的计算人员不是难事,但一些缺少计算经验的计算人员常因这方面处理不当而导致计算精度不足。通过了边坡算例分析提出了下面的计范围与网格划分的处理。

边界范围的取值在有限元法中对计算结果有较大影响。当坡角到左端边界的距离为坡高的1.5倍,坡顶到右端边界的距离为坡高的2.5倍,且上下边界总高不低于2倍坡高时,计算精度较为理想。

计算时必须考虑适当的网格密度,如果网格划分太粗,将会造成很大的误差。究竟单元大小取多大为宜,一般根据具体的问题来解决。可以先执行一个认为合理的网格划分的初始分析,再在可能出现滑面的危险区域利用2倍的网格重新分析并比较两者的结果。如果这两者给出的结果几乎相同,则认为前次划分的网格密度是合适的。网格划分过程中,还可以对重要部位进行局部加密,不重要的地方,可以稀疏一些,需要注意的是从密集到稀疏最好要有一个平缓的过渡,单元大小不要突然急剧变化,如图4所示。

《图4 》

图4 网格加密

Fig.4 FEM model with local grid refinement

用有限元计算岩土工程稳定问题时,不仅需要有几何参数、土容重y与抗剪强度等参数,还需要填入泊松比v、弹性模量E等变形参数。研究表明v对边坡的塑性区分布范围有影响,v的取值越小,边坡的塑性区范围越大。但是计算表明,v的取值对安全系数计算结果的影响极小。E对边坡的变形和位移的大小有影响,但是对于稳定安全系数基本无影响。由此可见,只需按经验来选取E,v,即使选取有所不当,也不会影响稳定分析的结果。当然,如果计算变形与位移,必须尽量选准E,v。

《4 有限元强度折减法在均质边(滑)坡中应用》

4 有限元强度折减法在均质边(滑)坡中应用

《4.1 在二维边坡中的应用》

4.1 在二维边坡中的应用[9,24~26]

一个平面应变情况下的算例。均质土坡,坡高H=20m,土容重y=20 kN/m3,粘聚力c=42 kPa,内摩擦角φ=17°,求坡角β=30°,35°,40°,45°,50°时边坡的稳定安全系数以及对应的滑动面。

4.1.1 有限元模型的建立和计算 计算采用大型有限元 ANSYS 5.61软件。按照平面应变建立有限元模型,边界条件为左右两侧水平约束,下部固定,上部为自由边界,如图5所示。 

《图5》

图5 β=30°时的有限元模型

Fig.5 FEM model

为了和传统方法作比较,强度折减安全系数的计算统一采用c/ω,tan φ/ω的折减形式,力和位移的收敛标准系数均取0.00001,最大迭代次数为1000次。一次性施加全部重力荷载,即荷载增量步设置为1步。

4.1.2 安全系数计算结果及其分析     表3为各屈服准则采用非关联流动法则时的安全系数,表4为各屈服准则采用关联流动法则时的安全系数。平面应变莫尔-库仑匹配D-P准则在关联和非关联流动法则条件下分别采用不同的表达式DP5与DP4,而对于莫尔-库仑等面积圆DP3准则和外角点外接圆DP1准则均采用同一种表达形式,只是使用关联与非关联法则时,两者采用的膨胀角不同。传统极限平衡条分法计算采用加拿大的边坡稳定分析软件 SLOPE/W。
在平面应变条件下不管是采用非关联的莫尔-库仑匹配DP5准则还是采用关联的莫尔-库仑匹配DP4准则,求得的安全系数与传统极限平衡条分法中的Spencer法的计算结果十分接近,误差在2%以内,这是因为平面应变莫尔-库仑匹配D-P准则实际上就是在平面应变条件下的莫尔-库仑准则。
对于平面应变问题,莫尔-库仑等面积圆DP3

《表3》

表3 采用非关联法则时不同准则条件下的稳定安全系数

Table 3 Safety factor by different method with associated flow rule

坡角/(°)

30

35

40

45

50

DP1

1.91

1.74

1.62

1.50

1.41

DP3

1.64

1.49

1.38

1.27

1.19

DP5(非关联流动法则)

1.56

1.42

1.31

1.21

1.12

极限平衡 Spencer法(S)

1.55

1.41

1.30

1.20

1.12

(DP1-S)/S

0.23

0.23

0.25

0.25

0.26

(DP3-S)/S

0.05

0.06

0.06

0.06

0.06

(DP5-S)/S

0.01

0.01

0.01

0.01

0.00

《表4》

表4 采用关联流动法则时不同准则条件下的安全系数

Table 4 Safety factor by different method with deviatoric flow rule

坡角/(°)

30

35

40

45

50

DP1

1.93

1.77

1.65

1.54

1.44

DP2

1.66

1.51

1.40

1.30

1.21

DP4(关联流动法则)

1.56

1.42

1.32

1.22

1.13

极限平衡Spencer法(S)

1.55

1.41

1.30

1.20

1.12

(DP1-S)/S

0.25

0.26

0.27

0.28

0.29

(DP3-S)/S

0.07

0.07

0.08

0.08

0.08

(DP4-S)/S

0.01

0.01

0.01

0.02

0.01

屈服准则,当使用非关联流动法则时,计算结果与传统极限平衡方法中的Spencer法计算结果误差在6%左右;当使用关联流动法则时,误差在7%左右。而外角点外接圆DP1准则条件下的安全系数比传统的极限平衡条分法中的Spencer法大25%以上。

4.1.3 边坡临界滑动面的确定     根据边坡破坏的特征,边坡破坏时滑面上节点位移和塑性应变将产生突变,滑动面在水平位移和塑性应变突变的地方,因此可在ANSYS程序的后处理中通过绘制边坡水平位移或者等效塑性应变等值云图来确定滑动面。下面给出一个算例,除用上述两种方法确定滑动面外,并与传统确定滑面的方法进行比较,算例表明,上述3种方法确定的滑面是一致的。坡角β=30°时的滑动面形状和位置见图6至图8,其中边坡变形显示比例设为0。

《图6》

图6 用等效塑性应变等值云图表示的滑动面位置和形状

Fig.6 The failure surface by FEM

《图7》

图7 用水平位移等值云图表示的滑动面位置和形状

Fig.7 The failure surface by FEM

《图8 》

图8 用加拿大边坡稳定分析软件 SLOPE/W得到的滑动面形状

Fig.8 The failure surface by SLOPE/W

4.1.4 有地下水渗流作用时有限元强度折减法[27]

国内对渗流作用下边坡稳定性的分析的研究刚刚开始,一般采用国际通用软件,但是各种软件在分析方法以及显示功能等方面有所不同。

1)应用ADINA有限元程序求解渗流作用时边坡稳定安全系数。利用ADINA有限元程序进行渗流作用下的边坡稳定性分析,可以通过有限元强度折减法分析坡体在不同浸润面时的边坡稳定性来实现。

为分析坡体内含浸润面时的边坡稳定性,必须考虑坡体内孔隙水压力对土体受力与变形的影响。在ADINA程序中,岩土材料(骨架)可采用任何一种线性、非线性岩土材料或用户自定义的材料;而孔隙水的属性则由多孔介质属性进行定义,其主要参数为各个方向的渗透系数。因此,利用ADINA程序分析坡体内含浸润面时的边坡稳定性是可行的。

影响坡体内浸润面位置的因素很多,包括边界水头条件。当其他条件不变时,不同的边界水头条件对应坡体内不同的浸润面位置。因此,通过边界水头条件来控制坡体内的浸润面位置。但由于ADINA渗流程序中,无法显示出浸润面的位置,因而先利用渗流方程和温度方程相同的原理,在温度场中求解 seepage材料,采用热传导单元求解不同边界(水头边界)条件下的浸润面的形状,并把它显示出来。

在实际工程中,如果已知坡体内的浸润面位置,则可以通过假设不同的水头边界对应已知的浸润面置得到真正的水头边界,然后按孔隙水压力的形式施加在计算模型上进行稳定性分析;如果不知道坡体内的浸润面位置,则可以通过水头边界模拟浸润面的位置,然后按孔隙水压力的形式施加在计算模型上进行稳定性分析。这样就能通过有限元强度折减法求解边坡的安全系数和滑面位置。

算例:均质边坡,坡高H=20m,粘聚力c=42 kPa,土容重y=20 kN/m3,内摩擦角φ=17°,渗透系数k=4x10-3cm/s。在两侧边界上都施加水头荷载Ha=20m。该算例所对应的浸润面计算模型和稳定性分析计算模型分别如图9和图10所示。

《图9 》

图9 浸润面计算模型

Fig.9 The calculation model of phreatic surface

《图10 》

图10 稳定性分析计算模型

Fig.10 The calculation model of slope stability analysis

所对应的温度场荷载通过ADINA程序温度场模块的求解,可以得到坡体内浸润面的位置如图11所示。

《图11》

图11 浸润面位置示意图

Fig.11 The graph of phreatic surface

由于ADINA程序在考虑渗流作用时,其孔隙水压力是以节点力的形式进行计算的,因此它可以充分考虑孔隙水压力在各个方向上的影响,故在浸润面位置确定的情况下,可以把其所对应的水头荷载按孔隙水压力荷载的形式施加在稳定性分析的计算模型上。虽然计算模型中不能显示浸润面的位置,但由于计算模型上施加的约束与荷载条件和利用温度场计算时的约束与荷载条件相同,因此计算模型中浸润面的位置也应该和利用温度场计算出的浸润面位置相同。然后,采用有限元强度折减法进行安全系数的求解。

由于渗流作用下产生的孔隙水压力使计算模型的每个结点都产生了不同程度的塑性应变,因此在ADINA程序所绘制的等效塑性应变云图比较紊乱,无法显示滑面的位置。通过计算得到边坡的安全系数为1.550。

为了验证ADINA有限元程序分析的准确性,通过传统条分法GEO-SLOPE程序的 SLOPE/W和SEEP/W耦合以及直接在SLOPE/W程序中考虑浸润面的位置对该算例的安全系数分别进行了验算,该方法采用的是Spencer法,计算得到的安全系数分别是1.523和1.525。通过与ADINA程序安全系数的计算结果进行对比,其误差都在2%以内。

2)应用PLAXIS有限元程序求库水下降时边坡稳定安全系数[27]。PLAXIS程序是荷兰开发的岩土工程有限元软件,应用性非常强,能够模拟复杂的工程地质条件,尤其适合于变形和稳定性分析。

PLAXIS程序中孔隙水压力的渗流计算是以有限元原理为基础的,通过定义土体的渗透系数、划分网格、定义边界条件,建立水力计算模型进行计算。

PLAXIS程序假设孔隙水在多孔介质中的运动符合Darcy定律,并设有排水力学条件与不排水力学条件的计算功能,前者土体内将不产生超孔隙水压力,而后者超孔隙水压力不消散。

结合PLAXIS程序在渗流计算方面的强大功能,就目前工程上比较关心的库水位下降对边坡稳定性的影响进行了分析,并对水位下降速率(r)的影响进行了研究。

算例:均质边坡,坡高H=20m,c=24kPa,坡角 arctan 1/2,土容重Ydn=15 kN/m3,γwet=18 kN/m3, φ=20°,v=0.35,E=2000kPa。

在渗流计算模型中认为坡体后部地下水补给充足,坡体后部边界水头值保持Ha=30m不变,坡体前部Hb水位从初始水位40m开始下降。

有限元模型和渗流计算模型的网格划分示意图和渗流计算的模型示意图如图12、图13所示。

由于水位的下降是从初始水位40m开始的,因此,在初始水位时坡体已经经过长期的浸泡,坡体内的超孔隙水压力p.已经完全消散,所以土体在初始水位时是符合排水条件下的力学行为。

《图12 》

图12 有限元模型和渗流计算模型的网格划分示意图

Fig.12 The grid of finite element module and seepage calculation module

《图13》

图13 渗流计算模型示意图

Fig.13 The module of seepage calculation

为了分析水位下降速率对边坡稳定性的影响,是通过在程序中设置固结天数的方法实现的。即如果水位从30m处按1 m/d的速率下降,则程序设置不排水的力学条件:水位从30m下降到29m,

孔隙水压力不消散。然后稳定1d再下降,在1d的时间内程序进行固结计算,从而使产生的超孔隙水压力消散1d。同样,如果水位下降的速率为2 m/d,则水位从30m下降到28m后稳定2d再下降,在2d的时间内也进行固结计算,从而使产生的超孔隙水压力消散2d。

上述可见,在PLAXIS程序中如果单纯地考虑水位的变化是无法考虑时间因素的,所以采用结合固结计算的方法来考虑时间因素,即把水位下降到一定高度所经过的时间通过固结计算的时间来体现。因此,对于不同的水位下降速率,即水位下降到一定高度所经历的不同时间,就可以通过设置水位下降到一定高度后,再进行不同时间的固结计算来体现。

对水位分别按1 m/d,2 m/d和4 m/d不同速率下降时的边坡稳定性进行研究,计算结果见表5(渗透系数kx=ky=4x10-2 m/d)。

《表5》

表5 坡体前部不同水位下降速率对应的安全系数计算结果表

Table 5 The calculation results of the safety coefficient at different water declining rates in the former part

-1

安全

系数(当坡体前部水位

高度

3/m)

d

r/md

40

36            32            28

24

20

kN·m-2

1

2.42

1.91   1.59 1.39

1.28

1.23

1.36

2

2.42

1.90 1.58 1.35

1.23

1.20

1.36

4

2.42

1.87 1.53 1.34

1.20

1.15

1.36

表5可见,水位下降速率越快,边坡的稳定性越差。见图14至图18的浸润面和滑面位置示意图。

《图14》

图14 水位按1m/d下降到32m时浸润面和滑面位置示意图(安全系数为1.59)

Fig.14 The graph of phreatic and sliding surfaces when the water level drops to 32 m at 1m/d

《图15》

图15 水位按1m/d下降到24m时浸润面和滑面位置示意图(安全系数为1.28)

Fig.15 The graph of phreatic and sliding surfaces when the water level drops to 24 m at 1m/d

《图16》

图16 水位按4m/d下降到32m时浸润面和滑面位置示意图(安全系数为1.53)

Fig.16 The graph of phreatic and sliding surfaces ,when the water level drops to 32 m at 4m/d

《图17 》

图17 水位按4m/d下降到24m时浸润面和滑面位置示意图(安全系数为1.20)

Fig.17 The graph of phreatic and sliding surfaces when the water level drops to 24 m at 4 m/d

从图19中可以看出,水位下降速率越快,边坡安全系数下降的幅度也越大,说明在水位下降过程中,水位下降速率越快对边坡的稳定性越不利。
为了进一步说明水位下降过程中由于坡体内地下水位下降的滞后效应所产生的超孔隙水压力对边坡稳定性的影响,计算了该算例在完全排水条件下所对应的边坡安全系数,即假设水位下降到每一高度时,坡体内的超孔隙水压力都得到了充分消散,计算结果见表6。可以看出,水位下降过程中坡体

《图18》

图18 坡体内超孔隙水压力消散至最小值时的浸润面和滑面位置示意图

Fig.18 The graph of phreatic and sliding surfaces when the excess pore water pressure dissipates to the minimum
 

《图19》

图19 不同的水位下降速率所对应的水位和安全系数的关系曲线

Fig.19    The curve of water level and safety coefficient at different water declining rate

内产生的超孔隙水压力对边坡安全系数计算结果的影响十分明显,超孔隙水压力的产生对边坡稳定性十分不利。

《表6》

表6 是否考虑坡体内超孔隙水压力的安全系数计算结果表

Table 6 The calculation result of safety factor whether considering the excess pore water

Po的影响

 安全

系数(当

坡体前部水位

高度H

,/m)

Pepmin

rwd=1m/d

40

36

32          28

24

20

/kN·m-2

考虑

2.42

1.91

1.59 1.39

1.28

1.23

1.36

不考虑

2.42

2.08

1.77 1.52

1.45

1.37

1.37

《4.2 在三维边坡中的应用》

4.2 在三维边坡中的应用[28]

在边坡稳定分析领域,二维方法是常用的手段,但在岩土工程中很多边坡问题都属于三维边坡问题。有关边坡稳定三维极限平衡方法,已有众多
研究成果。Duncan 曾列表总结了20篇文献资料[29],列举了这些方法的特点和局限性。

4.2.1 简化为平面应变问题的空间模型 建立一个可以简化为平面应变问题的空间模型,计算模型如图20所示,坡高20m,坡角45°,坡角到左端边界的距离为坡高的1.5倍,坡顶到右端边界的距离为坡高的2.5倍,且总高为2倍坡高,在z方向取30m。计算采用ANSYS-5.61程序,有限元模型的边界条件为底面固定约束,坡体侧面约束相应的水平位移。采用SOLID45号实体单元和关联流动法则。土坡计算参数为y=25 kN/m3,c=42 kPa,φ为变量。

《图20》

图20 均质土坡计算模型

Fig.20 The model of soil slope

计算结果表明(见表7),在三维边坡计算中采用莫尔-库仑等面积圆屈服准则DP3是可行的,它所得到的计算结果与二维情况下得到的结果基本一致。

《表7》

表7 不同屈服准则得到的安全系数

Table 7    Safety factor by different method

H=20m β=45° c=42 kPa

φ/ ()

0.1

10

25

35

45

DP1

0.523

1.072

1.696

2.105

2.497

DP2

0.522

0.938

1.303

1.473

1.494

DP3(三维)

0.475

0.920

1.390

1.680

1.925

DP3(平面)

0.455

0.915

1.388

1.665

1.914

4.2.2 有限元强度折减法与三维极限平衡法比较通过典型算例对比三维极限平衡法与有限元强度折减法的计算结果,验证有限元强度折减法应用于三维边坡的可行性。

图21为Zhang Xing提供的椭球滑面三维极限平衡法算例[30],国内都选择该算例来检验各自的三维极限平衡法程序的合理性。计算参数见图22,椭球的长宽比为1:3,按原例要求,在对称轴平面

用一圆弧模拟滑裂面,在z方向,则以椭球面形成滑面,即滑面是给定的,在滑面四周约束土体的位移。对此算例采用不同的屈服准则计算,见表8。

《图21》

图21  椭球体滑面算例

Fig.21 An example with an ellipsoid failure surface

《表8》

表8 Zhang Xing 算例用不同屈服准则得到的稳定安全系数

Table 8 Safety factor by different method about The example of Zhang Xing  

屈服准则

DP1

DP2

DP3

Zhang Xing

安全系数

2.489

2.217

2.150

2.122

误差/%

17

5

1

 

可见,采用DP3准则计算所得的稳定安全系数与Zhang Xing求得的稳定安全系数非常接近,表明在分析三维边坡时采用莫尔-库仑等面积圆屈服准则是适宜的。封底图1为其计算不收敛时的x方向的位移云图。

在Zhang Xing原例中约束了滑面周围土体的位移,这种假定不符合实际情况。实际情况下滑面周边的土体并未约束,笔者对这种情况进行了计算。w/l=3椭球滑面计算所得安全系数为2.165,也与给定滑面情况下的安全系数2.15十分接近。

目前,虽然已有很多三极限平衡法的分析程序,但三维极限平衡法与二维相比作出了更多的假定,而且需要给定滑面。而有限元强度折减法不需要作任何假定,计算结果更可靠,它为三维边坡稳定性分析开辟了新的途径。

大型水电站岩石高边坡大都属于三维边坡。目前,有些设计部门正尝试将有限元强度折减法用于水电站的三维边坡的稳定分析。

《4.3 在确定滑坡多个滑面中的应用》

4.3 在确定滑坡多个滑面中的应用[3]

4.3.1 概述    确定滑坡滑动面位置和形状的传统方法主要是在现场钻探的基础上,通过技术人员的分析判断提出滑带位置。这种判断方法存在如下问题:a.当只有少量钻孔发现滑带特征时,依据少量滑带位置来判定整个滑带有时可能出现差错;b.当滑坡体处于蠕变阶段,滑面尚未形成,更无法通过勘查找出滑面;c.即使查明了滑带和剪出口,还可能存在一些次级滑面和次生剪出口。目前的判断方法容易造成滑面遗漏。因为次生滑面的遗漏而导致工程失败,已有了多次的教训。

4.3.2 寻找滑坡中的多个滑面    对一个复杂滑坡算例,通过依次约束已知滑面剪出口的方式,搜索出低于设定稳定安全系数的所有滑动面,可纠正错划、漏划现象,从而全面、准确地确定出复杂滑坡内的多个滑动面,为滑坡治理方案提供科学依据。
1)模型和计算参数 滑坡断面如图22所示,滑坡材料的物理力学参数见表9。

《图22》

图22 滑坡治理示意图

Fig.22 Managing landslide

《表9》

表9 材料物理力学参数

Table 9 Material properties

材料名称

y/kN·m-'

 E/MPa

V

c/kPa

(。)/

滑体

20.7

30

0.3

28

22.2

滑带

20

30

0.3

25

19.4

滑体下伏稳定岩层

23.7

1600

0.2

200

32.0

滑体、滑带和下伏稳定岩层均采用6结点二次三角形平面单元模拟。首先用有限元强度折减法计算,表明滑面位于滑带上,滑坡的稳定安全系数为1.00,而用极限平衡法(spencer)算得滑坡稳定安全系数为1.002。

2)通过约束剪出口寻找其他潜在滑面 滑坡治理的过程实际上是滑动面变化与稳定安全系数提高的动态过程。对于复杂滑坡,必须考虑多个次生滑动面的出现,只有所有次生滑动面的稳定安全系数都达到设计规定,该滑坡从工程意义上来说才是
安全的。算例中规定滑坡的稳定安全系数为1.20,为了寻求可能出现的多个次级滑动面的位置及滑动次序,在有限元计算中采用约束坡体剪出口附近某一地段的水平位移来表示对该部分的治理。依次对未达到设定稳定安全系数的所有剪出口进行约束,求出相应滑面与稳定安全系数,直至稳定安全系数达到设计规定值。表10列出了所有可能的滑面及滑面滑动先后顺序、剪出口位置与稳定安全系数,图22至图26列出了各滑面位置。

《表10》

表10 约束部位与滑坡稳定安全系数间的关系

Table 10 Parts restrained and safety factors of landslide  

序号

约束部位

滑面产

剪出口

稳定安

备注

生次序

位置

全系数

1

天然滑坡

1

A以上

1.000

 

2

ABC

2

C以上

1.023

滑坡设定稳定

3

+CDE

3

E以上

1.052

 

 

 

 

 

 

安全系数1.20

4

+EFG

4

M以上

1.135

 

5

NW+

5

G以上

1.203

 

《图23》

图23 约束ABC段滑坡极限状态的滑动面(F=1.023)

Fig.23 Landslide critical surface under restraining ABC part(F=1.023)

《图24》

图24 增加约束CDE段滑坡极限状态的滑动面(F=1.052)

Fig.24 landslide critical surface under restraining CDE part(F=1.052)

由此可见,复杂滑坡可能存在多个剪出口和滑

《图25》

图25 增加约束EFG段滑坡极限状态的滑动面(F=1.135)

Fig.25 Landslide critical surface under restraining EFG part(F=1.135)

《图26》

图26 增加约束MN段滑坡极限状态的滑动面(F=1.203)

Fig.26 Landslide critical surface under restraining MN part(F=1.203)

动面,而这些滑动面都未达到设计规定的稳定安全系数。显然,上述支护方案均不能达到工程治理的目的。因而合理的滑坡支护方案必须寻求出所有小于设计规定的稳定安全系数的滑动面与剪出口,抑制所有滑动面的贯通和剪出口的剪出,才能达到治理滑坡的目的。

《5 在岩质边坡中的应用》

5 在岩质边坡中的应用 [13,24,32,33]

《5.1 概述》

5.1 概述

岩质边坡的稳定分析历来是至为关注的重大课题,由于实际岩体中含有大量不同构造、产状和特性的不连续结构面(如层面、节理、裂隙、软弱夹层、岩脉和断层破碎带等),这就给岩质边坡的稳定分析带来了巨大的困难。岩质边坡的稳定性主要由岩体结构面控制,传统的用于土质边坡稳定分析的滑动面搜索方法很难用于岩质边坡。

岩体中的结构面,根据其贯通情况,可以将结构面分为贯通性、非贯通性两种类型。根据结构面的强弱和充填情况,可以将结构面分为硬性结构面和软弱结构面。由于岩体结构的复杂性,要十分准确地反映岩体结构的特征十分困难。基于这种考虑,对于一个实际工程来说,往往根据现场地质资料,根据结构面的长度、密度、贯通率,展布方向等,着重考虑2~3组对边坡稳定起主要控制作用的节理组或其他结构面。

《5.2 岩质边坡结构面模型》

5.2 岩质边坡结构面模型

1)软弱结构面有限元模拟 软弱结构面和岩体均采平面实体单元模拟,按照连续介质处理,只是结构面与岩体选用的参数不同而已。岩体以及结构面材料本构关系采用理想弹塑性模型,强度折减过程与均质土坡相同,即通过对岩体以及结构面强度参数同时进行折减使边坡达到极限破坏状态,此时可得到边坡的强度储备安全系数。

2)硬性结构面有限元模拟 无充填的硬性结构面可采用无厚度接触单元来模拟,程序通过覆盖在两个接触物体表面的接触单元来定义接触关系。两个接触面之间不抗拉,可以脱离,可以滑动。两个接触面的接触摩擦行为服从库仑定律:

\(\tau=c+\sigma \tan \varphi,  \sigma \geqslant 0 \) (规定压为正)。

《5.3 求两组贯通结构面边坡的滑面及安全系数》

5.3 求两组贯通结构面边坡的滑面及安全系数

如图27所示,两组方向不同的结构面,贯通率100%,第一组软弱结构面倾角30°,平均间距10m,第二组软弱结构面倾角75°,平均间距10m,岩体以及结构面计算物理力学参数见表11。

《图27》

图27 由有两组贯通的平行结构面控制的岩质边坡几何模型

Fig.27 Geometry model

岩体以及结构面均采用平面应变单元模拟,只是物理力学参数不同,计算步骤同上,通过有限元强度折减得到的破坏过程如图28所示,图28a是岩质边坡最先产生的破坏形式,接着出现第二条、第

《表11 》

表11 物理力学参数计算取值

Table 11    Material properties

材料名称

y/kN·m-3

E/MPa

V

c/ MPa

φ/(°)

岩体

25

10

0.2

1.0

38

第一组结构面

17

10

0.3

0.12

24

第二组结构面

17

10

0.3

0.12

24

三条次生滑动面(图28b),求得的稳定安全系数见表12,其中极限平衡方法计算结果是根据最先贯通的那一条滑动面求得的。计算表明,采用DP1准则计算结果偏大,采用DP3准则计算结果与极限平衡法算得的结果相近。

《图28》

图28 极限状态后产生的破坏形式

Fig.28 Failure progress of rock slope

《表12》

表12 岩质边坡稳定安全系数计算结果

Table 12 Safety factor by different method

计算方法

安全系数

有限元法(外角点外接圆DP1准则)

1.49

有限元法(莫尔-库仑等面积圆DP3准则)

1.21

极限平衡方法(spencer)

1.17

5.4 具有非贯通结构面岩质边坡的稳定分析

对于具有非贯通结构面岩质边坡,也可以采用有限元强度折减法来模拟其破坏机制。如图29a所示,结构面AB,CD,FG倾角均为45°,AB=21.21m,CD=14.14m,CE=35m,AD=10m, BAD=135°,AF=20m,FD=22.36m,计算采用的参数见表13。计算机自动使结构面从A与D之间贯通,如图29b所示,对应的安全系数为2.7。

《图29 》

图29 具有3条非贯通结构面岩质边坡破坏模式

Fig.29 Failure progress of rock slope with structural plane

《表13》

表13 物理力学参数计算取值

Table 13    Material properties

材料名称

y/kN·m-3

E/Pa

V

e/MPa

φ/()

岩体

25

1010

0.2

1.0

30

结构面

18

10

0.3

0.06

18

研究表明,在岩体及结构面参数相同的情况下,结构面之间的贯通破坏机制受结构面几何位置、倾角、结构面之间岩桥的倾角、岩桥长度等因素的影响。结构面及岩桥位于受剪力最大的地方容易贯通;在垂直边坡中,结构面倾角愈接近45°+φ/2就愈容易贯通;岩桥长度愈短愈容易贯通。无论何种状态的非贯通结构面,计算机都可自动获得首先贯通的结构面及其安全系数。计算中,岩体结构面的贯通率是一个重要参数。

《6 用有限元法进行边(滑)坡支挡结构的设计》

6 用有限元法进行边(滑)坡支挡结构的设计 [34~37]

《6.1 概述》

6.1 概述

用传统方法进行边(滑)坡支挡结构设计,通,先采用极限平衡法确定支护结构上的推力(岩土侧压力),但要求事先准确确定破坏滑动面的位置与形状,并采用合理的方法计算,这样才能准确算出岩土压力,还要明确岩土压力如何分布在支挡结构上。传统计算方法中,作用在支挡结构上的岩土压力分布是假定的,一般假设为矩形、三角形或梯形分布。假定不同的分布形式对支挡结构内力计算有很大差异,因而传统算法会有较大误差。然后,采用弹性地基梁法计算抗滑桩结构的内力,弹性地基梁法作了人为假设,而且地基系数也难以准确确定。采用有限元法可以严格按照弹塑性理论计算,并充分考虑岩土与支挡结构的共同作用,不需要对边(滑)坡推力的分布作任何假定,还可以直接计算出支挡结构内力。尤其当采用锚桩支护等复杂结构时,有限元法可以准确地进行结构优化计算,而传统方法是在假定桩上推力分布的基础上进行优化,因而优化计算的可信度很低。

用有限元法进行支护结构的设计计算,一般包括4个步骤:a.对抗滑桩上的推力进行验算;b.获得抗滑桩上的推力分布形式;c.计算抗滑桩的内力;d.进行结构优化设计。

《6.2 工程算例》

6.2 工程算例

6.2.1 工程概况 国道主干线重庆-湛江公路在贵州境内的崇溪河至遵义高速公路,高工天滑坡位于第五合同段K26+150-K26+260段,路基开挖时,下切滑体5~6m,即引起滑坡复活,而且还在不断发展,形成多级的滑面,发育在土层和强风化带内。如果按照设计开挖切脚,必将引起岩体滑动。根据设计,该滑坡的治理采用抗滑桩加预应力锚索的支挡措施,每根锚索设计锚固力800kN,每根桩上纵向布置两排锚索,每排3根,共6根,计算采用的典型断面如图30所示。

《图30》

图30 计算采用的典型断面

Fig.30 The analysis cross section of landslide

6.2.2 有限元模型的建立 计算采用的软件为美国ANSYS公司的大型有限元软件ANSYS(R)5.61商业版。有限单元网格划分见图31,计算按照平面应变问题建立模型,岩土体采用8节点平面单元PLANE183模拟,抗滑桩用梁单元 BEAM3单元模拟,桩的断面积、惯性矩等可以在其对应的实常数中定义,该单元可以输出轴力、弯矩、剪力等。

当抗滑桩与锚杆(索)联合使用或者单独使用作为边(滑)坡的支护结构时,采用有限元计算充分考虑了锚杆(索)、桩与岩土介质的共同作用。

《图31》

图31 滑坡有限元模型示意图

Fig.31 Sketch of Finite element mesh

一般锚杆不施加预应力,属于被动式支护,可采用杆单元模拟。而锚索一般是施加预应力的,属主动式支护,其施加的预应力,一般就是锚索的设计锚固力。传统的做法是在锚索的两锚固点,施加一对压力代表锚固力。这种情况下,锚索的作用力与岩土介质的变形无关。为了更好模拟锚索作用,也可采用杆单元来模拟锚索,锚索的预应力可以通过设置初应变来获得,初应变要根据设计锚固力来反算。施加预应力锚索后,随着滑体强度参数的降低,锚索的受力会逐渐增大,当锚索受力大于锚索设计抗拉强度时,锚索失效。桩与滑体之间的接触关系分别采用两种方案:方案一采用ANSYS程序提供的接触单元来模拟桩与土的接触行为;方案二为桩与土共节点但材料性质不同的连续介质模型。

根据设计,每根锚索设计锚固力800kN,每根桩上布置两排锚索,而该平面应变模型在纵向只布置了一排锚索,所以将锚索锚固力乘以2,即1600kN。锚索倾角10°,在有限元模型中,在锚索的外锚头节点的水平方向施加-1600cos 10°kN,在竖直方向施加-1600 sin 10°kN。抗滑桩截面为3mx4m。
锚索纵向间距为4m,而本次平面应变计算纵向只有1m,也就是说每根桩要承担4m宽的滑体的剩余下滑力,因此在有限元模型中可将土体重量乘以4,同时为了确保原有稳定安全系数不发生变化,将岩土体的粘聚力也乘以4,即保证y/c不发生变化。

6.2.3 材料本构模型及计算参数 岩土材料本构模型采用理想弹塑性模型,屈服准则采用平面应变摩尔-库仑匹配准则DP4,计算参数见表14。

6.2.4 开挖和支护过程的模拟 开挖和支护采用单元的“死活”来实现。所谓单元“杀死”,就是将单元刚度矩阵乘以一个很小的因子(10-6),死单元的荷载将为0,从而不对荷载向量生效,同

《表14》

表14 计算采用物理力学参数

Table 14    Material properties

材料名称

y/kN·m-3

E/MPa

V

c/kPa

(。)/6

滑体

21

30

0.30

25.5

24.5

滑床

24

10

0.25

200

30.0

桩(C25砼)

24

29x10

0.20

考虑为

弹性材料

样,死单元的质量也设置为0,单元的应变在“杀死”的同时也将设为0。与上面的过程相似,桩的施加采用单元的“出生”来模拟,并不是将单元增加到模型中,而是重新激活它们,其刚度、质量、单元荷载等将恢复其原始的数值,重新激活的单元没有应变记录,所有单元都要事先划分好。根据现场实际施工过程,有限元计算分4步:

1)计算未开挖前的初始应力场;

2)施工桩,激活桩元,同时施加锚固力;

3)开挖,杀死要开挖的土体单元;

4)取边坡稳定安全系数1.2,滑体强度参数折减1.2倍后计算。

6.2.5 滑坡推力计算与验算 滑坡推力安全系数采用强度贮备系数的定义,即将滑体强度参数折减1.2倍后计算滑坡水平推力。利用ANSYS软件提供的路径分析功能,沿桩从滑面到顶部设置路径,将水平应力映射到路径上,然后沿路径对水平应力进行积分,就可以得到总的滑坡水平推力,计算结果见表15。

《表15》

表15 不同方法计算得到的滑坡水平推力(kN)

Table 15 Landslide thrust force by different method  

 

极限平衡法

FEM

连续介

不平衡推力法

Spencer

6770

6 440

6 944

6 400

从表15看出,采用连续介质有限元模型的计算结果与接触单元模型中桩土粗造接触模型的计算结果比较接近,说明如果桩土之间没有明显的滑动时可以采用连续介质模型来模拟桩和土的接触关系,这样操作方便。总体看来,有限元法计算的推力与传统极限平衡方法计算结果比较接近,因而可采用有限元法中连续介质模型来计算支挡结构的内力。

6.2.6 滑坡推力分布 通过有限元法得到的滑坡推力分布如图32所示,推力呈弓形分布。当坡面倾角较大时,实测的推力分布一般也呈拱形分布。

《图32》

图32 用有限元法得到的滑坡推力分布

Fig.32 Landslide thrust pressure distribution by FEM

6.2.7 抗滑桩弯矩和剪力 采用上述方法计算得到施加锚固力后桩的最大弯矩为11900kN-m,最大剪力为2650kN,弯矩和剪力分布分别见封底图2、封底图3。

从表16看出,传统方法中采用不同的滑坡推力分布图式的计算结果有很大的差别,有限元计算结果与传统方法中滑坡推力分布假定为矩形时的计算结果相对接近;三角形分布计算结果偏于危险。另外通过锚索施加锚固力后,桩的弯矩和剪力都大大减小,可见锚索和抗滑桩联合使用显著地改变了桩的悬臂受力状态,可以使桩的截面积显著减小,大大地节约工程材料。

《表16》

表16 不同方法计算结果对比

Table 16 Internal force of pile by different method

传统方法

 

有法

无预应力锚索

剪力/kN

6276

8323

6 560

弯矩/kN·m

42 062

58 082

48 100

剪力/kN

875

1756

2650

有预应力锚索

 

 

 

弯矩/kN·m

5346

11 310

11 900

注①为传统抗滑桩计算中假定滑坡推力分布为三角形,②假定为矩形

6.2.8 支挡后安全系数 该滑坡采用预应力锚索加固后,如果锚索和桩不出问题,随着滑体强度参数的降低,就会出现如图33所表示的滑动面,滑动面出现在桩顶,滑体越过桩顶滑出。此时算出的稳定系数为1.39,大于安全系数1.2,表明滑体不会出现“越顶”破坏。传统方法中往往忽视这一验算,有可能出现“越顶”破坏。

6.2.9 锚固力优化 下面分别计算不同锚固力时

《图33》

图33 加固后的滑动面

Fig.33 The failure surface with anti-slide pile

桩的弯矩,以进行锚桩结构的优化,计算结果见表17。

《表17》

表17 采用不同锚固力时桩的弯矩

Table 17 Bending moments of pile with different cable force

 

每孔锚索

 

桩的弯矩/kN·m

 

 

锚固力/kN

有限元法

传统方法

传统方法

1

600

19 700

7 853

22 683

2

800

11 900

5346

11 310

3

900

4550

14 516

5 583

4

950

2650

17 249

2 967

5

1000

3 410

19 982

4532

6

1 100

7 300

25 447

8 110

7

1 200

11 700

30 913

13 575

注①为传统计算中假定滑坡推力分布为三角形,②假定为矩形

图34为不同锚固力时桩的弯矩变化曲线,从计算结果看出,锚固力并不是越大越好,存在一个极值。从桩的弯矩变化曲线的走势看,当锚固力变化时,有限元计算结果与传统方法中滑坡推力分布假定为矩形时计算结果的变化趋势接近,而三角形分布假定的计算结果则差异很大。当每根锚索的锚固力为900kN或950kN时,有限元计算结果为4 550 kN·m或2650kN·m,与传统方法计算结果5 583 kN·m或2967 kN·m比较接近,而且内力值小,因此经过比较分析认为每孔锚固力采用900kN或950kN为宜。

锚索的锚固力大小对桩的内力有较大影响,设计中可以通过不同方案对比进行优化设计,使结构更趋经济安全。由于锚固力会有衰减,很难准确确定锚固力,桩设计中必须考虑这点。

《图34》

图34 不同锚固力时桩的弯矩曲线

Fig.34 Relationship curve of bending moments with different anchor cable force

《7    有限元超载法在土基上应用》

7    有限元超载法在土基上应用[14,38,39]

《7.1 在光滑刚性条形地基上的极限承载力》

7.1 在光滑刚性条形地基上的极限承载力

光滑刚性条形地基的极限承载力问题的基本理论基础源于Prandtl解,对一个承受均匀垂直荷载的半无限、无重量地基,其极限承载力可以通过极限分析求得其精确解为

\(q_{\mathrm{u}}=c \cos \varphi\left[\exp (\pi \tan \varphi) \tan ^{2}(\pi / 4+\varphi / 2)-1\right]\)

\( \varphi=0  \)的情况下\(  q_{\mathrm{u}}=(\pi+2) c \)

下面通过算例用有限元超载法求解极限承载力,并验证其正确性。

有限元模型如图35所示。强度参数c=10,φ=0°~30°,采用平面应变摩尔-库仑匹配DP4准则(关联流动法则)进行求解。不同强度参数条件下的极限荷载的计算结果见表18。

《图35》

图35 有限元网格划分

Fig.35 FEM meshing

从表18看出,采用极限分析有限元法得到的极限承载力与Prandtl解非常接近。计算也可以采用非关联流动法则的平面应变摩尔-库仑匹配DP5准则,两者计算结果也十分接近。

《表18》

表18 极限承载力计算结果

Table 18    The results of bearing capacity  kPa

(o) /

Prandtl 理论解

有限元法

相对误差/%

0

51.42

52.19

1.50

5

64.89

65.96

1.66

10

83.45

84.98

1.83

15

109.77

111.90

1.94

20

148.35

151.75

2.29

25

207.21

212.08

2.35

30

301.40

310.00

2.85

图36为极限状态后的基础附近局部的位移矢量图(φ=0),图37为极限状态后的基础附近局部的等效塑性应变等值云图。图38为Prandtl解的破坏机构图(表19)。表20为Prandtl破坏机构的有限元解。由表19、表20可见,极限分析有限元法求解的破坏形式与Prandtl解的破坏形式非常一致。

《图36》

图36 位移矢量图

Fig.36 Displacement vector 

《图37 》

图37 极限状态时的破坏滑动面

Fig.37 The failure surface in limit status

图39为φ=0时地基中心点在增量加载全过程中的荷载一位移响应,可以看出,在地基的增量加载过程中,随着荷载的逐渐施加,地基顶部的位移也逐渐增大。当地基局部进入塑性状态后,位移增大得越来越快,当地基处于极限塑性状态时,位移将发生突变。

《表19》

表19 Prandtl 破坏机构有关参数(B0=1)

Table 19 some parameters about Prandtl failure mechanism(B0=1)

φ/()

0

5

10

15

20

25

30

m/'p

0.50

0.55

0.60

0.65

0.71

0.79

0.87

d2/m

0.71

0.79

0.89

1.01

1.16

1.35

1.59

h/m

1.00

1.25

1.57

1.99

2.53

3.27

4.29

《表20》

表20 Prandtl破坏机构的有限元解

Table 20 FEM results of Prandtl failure mechanism(B0=1)

φ/()

0

5

10

15

20

25

30

m/'p

0.49

0.53

0.60

0.65

0.70

0.75

0.89

d2/m

0.70

0.80

0.90

1.05

1.19

1.35

1.62

h/m

0.98

1.25

1.50

1.92

2.51

3.15

4.20

《图39》

图39 地基在增量加载全过程中荷载一位移响应

Fig.39 The load-displacement response of foundation in step-loading process

《7.2 求解有重地基的极限承载力》

7.2 求解有重地基的极限承载力

对土重的地基,其极限承载力在理论上无精确的解答。学者们提出了各种经验公式,主要有汉森、太沙基,梅耶霍夫,魏锡克等关于土重地基的极限承载力N,的经验公式,表21为增量加载有限元法与各经验公式对比,可见魏锡克经验式较为准确。

《7.3 节理岩石地基极限承载力求解》

7.3 节理岩石地基极限承载力求解

对复杂情况下的地基,如复杂几何边界、复杂边界荷载、含有结构面等非均质地基,传统诸多方法均不宜采用,而极限分析有限元法能方便应用。

《表21》

表21 Nγ 的有限元解与各经验公式对比

Table 21 Nγ got by FEM contrast to some classical empirical formulas

φ/()

5

10

15

20

25

汉森、太沙基

0.089

0.467

1.418

3.537

8.109

梅耶霍夫

0.069

0.366

1.129

2.870

6.765

魏锡克

0.449

1.224

2.647

5.386

10.876

FEM

0.502

1.557

3.422

6.249

12.325

地基岩块参数为c1=1.0MPa,φ1=40°。节理基本参数为c2=0.1MPa,φ2=10°。地基宽度为B,节理通过地基正下方2B深处,倾角为30°。节理用实体单元,屈服准则选用平面应变摩尔-库仑匹配DP4准则(关联流动法则)。有限元计算结果qu=18.14 MPa,若此地基中不存在节理,则其极限承载力的计算结果为77.75 MPa,可见节理的存在大大衰减了地基的极限承载力。图40为极限状态后的基础附近局部的等效塑性应变等值云图,图41为极限状态后的基础附近局部的位移矢量图。文献[39]列出了基础附近不同位置,不同强度,不同倾角结构面对地基极限承载力的影响。

《图40》

图40 极限状态时的破坏滑动面

Fig.40 The failure surface in limit status

《图41 》

图41 位移矢量图

Fig.41 Displacement vector

《7.4 地基承载板载荷试验有限元数值仿真模拟》

7.4 地基承载板载荷试验有限元数值仿真模拟

载荷试验是目前世界各国用以确定地基承载力的最主要方法,在地基处理效果检验中被广泛地采用,然而载荷试验在实际操作过程中也存在不少问题,如尺寸效应、费用高、工期长等。新兴的有限元等数值方法克服了以往耗时量大、操作复杂等缺点,有望逐渐成为岩土工程领域有效的实用方法之一。
土体的抗剪强度参数c,φ值是数值模拟中的关键参数,可通过现场直接剪切试验、室内直接剪切试验或室内三轴剪切试验来确定。图42及表22为某工程地基某一载荷试验点p-s试验结果。在压力380kPa前,p-s曲线近似直线;曲线出现第一拐点之后,p-s曲线呈逐渐增加到960kPa。在压力加到1.06MPa时,承压板周边土体出现明显隆起,表明地基土体已达到破坏,终止试验,此时,曲线呈下凹型。因此,该点地基土的极限承载力按规范取为960kPa,对应的沉降为29.16mm,比例界限取为380kPa。由于比例极限小于0.5倍的极限承载力,故取承载力特征值为380kPa。

《表22》

表22 试点压力一位移关系

Table 22 The relation about p-s  

压力/kPa

150

270

380

500

610

位移/mm

2.69

5.14

8.12

12.41

16.38

压力/kPa

730

840

960

1060

 

位移/mm

20.74

24.98

29.16

33.20

 

《图42》

图42 试点p-s曲线

Fig.42 p-s relation

该试点的相关参数:y=22kN/㎡,c=32.5kPa,φ=30.2°。泊松比v=0.27,土体的变形模量E=17 MPa。下面用增量加载有限元对载荷试验的过程进行模拟,有限元模型及网格剖分如图43所示,屈服准则选择平面应变摩尔-库仑匹配DP4准则(关联流动法则),应用理想弹塑性本构模型进行求解。
应用增量加载模拟载荷试验过程,共分12个载荷步,具体情况见表23所示。通过有限元求得其极限荷载值为1.14 MPa。

表24及图44为载荷试验p-s的有限元模拟结果,可以看出,计算结果与试验结果基本一致,

《图43》

图43 有限元模型及网格剖分

Fig.43 FEM model and meshing

《表23》

表23 增量加载过程

Table 23    The step-loading process

载荷步数

1

2

3

4

5

6

荷载增量/Mpa

0.15

0.12

 0.11

 0.12

0.11

0.12

载荷步数

7

8

9

10

11

12

荷载增量/MPa

0.11

0.12

 0.10

 0.07

0.01

0.01

在压力380kPa前两曲线吻合的很好,至960kPa这一段两曲线稍微有点出入,而960kPa后两曲线又吻合的很好。这说明只要E,v,c,φ等强度参数比较接近实际情况,载荷试验中沉降计算也与实测比较接近。

《表24 》

表24 试点压力-位移有限元计算结果

Table 24 The FEM results of p-s

压力/kPa

150

270

380

500

610

730

位移/mm

2.68

5.02

7.82

11.23

14.80

18.86

压力/kPa

840

960

1060

1 130

1140

1150

位移/mm

23.30

28.71

34.11

40.52

44.24

80.64

《图44》

图44 有限元模拟的p-s曲线与载荷试验的对比

Fig.44 The p-s relation got by FEM contrast to plate loading test

《8 有限元强度折减法在公路隧道中的应用》

8 有限元强度折减法在公路隧道中的应用 [40]

《8.1 引言》

8.1 引言

对隧道的稳定性评价一直缺乏一个合适的评判指标,传统的有限元法无法算出隧道工程的安全系数和围岩破坏面,仅凭应力、位移、拉应力区和塑性区大小很难确定隧道的安全度。当前工程上尚没有隧道稳定安全系数的概念,一般按照经验对隧道围岩的稳定性进行分级。极限分析有限元法在边坡稳定分析中取得了成功,笔者尝试将极限分析有限元法应用于求解隧道的稳定安全系数。从实际观察到的情况看,隧道受剪破坏的安全系数可分为两种:a.把围岩视作等强的均质体,引起隧道整体失稳,其相应的是整体安全系数;b.分别考虑围岩体的岩块强度与结构面强度,引起隧道周边局部失稳,一般发生在节理裂隙岩体中,其相应的是局部安全系数。笔者做一种探索性的尝试,只限于研究受剪破坏的整体安全系数,可采用有限元强度折减法求隧道安全系数与潜在破坏面。

《8.2 用有限元强度折减法求公路隧道的整体安全系数与潜在破坏面》

8.2 用有限元强度折减法求公路隧道的整体安全系数与潜在破坏面

8.2.1 工程概况 某半圆拱形公路隧道尺寸为9.4mx8.5m(宽x高),埋深50m,洞室所处位置岩体完整性较好,主要为花岗岩,根据国标《工程岩体分级标准》GB50218-94,分别属于II,III,IN类围岩。

计算准则采用摩尔-库仑等面积圆屈服准则DP3,按照平面应变问题来处理,边界范围取底部及左右两侧各4倍洞室跨度。岩体力学参数见表25,下标“上”、“下”分别表示围岩的上、下限值。

《表25》

表25 岩体物理力学参数

Table 25    Rock parameters with different rock mass

围岩类别

E/ GPa

V

y/kN·m-3

φ/()

c/ MPa

I

30

0.22

2700

60

2.0

II

20

0.25

2 700

50

1.5

I

10

0.30

2500

39

0.7

N

5

0.35

2400

27

0.35

8.2.2 计算结果与分析    隧道围岩破坏状态下塑性区分布如封底图4至封底图9所示,整体安全系数w1见表26。隧道稳定安全系数是指隧道整体安全系数,即把非等强度的真实岩体视为均质等强的岩体,据此求出安全系数。隧道中算出的塑性区往往是一大片,不像边坡岩土体中存在明显的剪切带,因而要找出围岩内的破坏面比较困难。但破坏时滑面上必将发生塑性应变和位移的突变。据此,只要找出围岩塑性应变发生突变时,将塑性区各断面中塑性应变值最大点的位置连成线,就可得到围岩的潜在破坏面。封底图4至封底图6给出了不同参数下围岩的塑性区及破坏面,破坏面为黑色点滑线。该线就在图中表示应变最大的黑色云图的周边。

从封底图4至封底图9塑性应变等值云图及其标尺可以看出,达到破坏状态时,II-类围岩的塑性区范围最大,但是破坏范围很小,安全系数最高;II-类围岩塑性区范围次之,破坏范围较小,安全系数较低;IV-类围岩塑性区范围最小,但是破坏范围最大,安全系数最小。这说明破坏状态下质量较好的岩体如II类围岩,塑性区即使出现一大片也可能保持整体稳定,而且破坏区也只是局部一小部分;相反,质量较差的岩石如N类围岩,塑性区范围并不很大就出现失稳,而且破坏区连成了一片,安全系数最低(见表26)。由此表明,单纯根据塑性区范围大小来评判隧道的安全性是值得商権的。上述参数中泊松比v与剪切强度c,φ值都会影响塑性区大小,岩质好的II类围岩破坏时的塑性区大于岩质差的II类、IV类围岩的塑性区很多,主要是因为II类围岩的泊松比小于皿类、IV类围岩的缘故。

《表26》

表26 不同围岩类别条件下的安全系数

Table 26    The safety factor with different rock mass

围岩类别

埋深/m

泊松比

安全系数

II

50

0.25

4.23

I

50

0.30

2.61

N

50

0.35

1.85

皿下

50

0.25

2.63

N

50

,0.25

1.87

II

150

0.25

2.05

Ш下

150

0.30

1.52

NF

150

0.35

1.19

不同围类别条件下的安全系数见表26。由表26可见,不同泊松比的条件下安全系数变化不大,但计算表明塑性区有较大变化,说明不能仅凭塑性区大小来评判围岩的稳定性。

将上述隧道的埋深变为150m,各类围岩的安全系数见表26,与上覆岩体厚度为50m相比,上覆岩体厚度对隧道的塑性区分布范围和安全系数有较大影响,同类围岩上覆岩体厚度越大,安全系数越小。这说明隧道的稳定性与埋深有很大关系,许多深层煤巷出现很大的地压就是例证。