非连续变形法在重力坝深层抗滑稳定分析中的应用

王 静，覃克非，王 敏

(1.四川省水利水电勘测设计研究院，成都，610072；2.滕州市中等职业教育中心学校，山东 滕州，277500)

非连续变形法在重力坝深层抗滑稳定分析中的应用

王 静1，覃克非2，王 敏3

(1.四川省水利水电勘测设计研究院，成都，610072；2.滕州市中等职业教育中心学校，山东 滕州，277500)

重力坝的深层抗滑稳定问题是影响坝体安全、工程投资及工程建设工期的一项重大问题，需对其失稳模式及安全系数进行严格的、符合实际的计算。目前常用的刚体极限平衡法及有限元法各有其局限性，很多方法计算前，需首先根据经验判断滑动模式、选取危险的滑动组合，存在人为主观因素，可能遗漏更危险的破坏模式。非连续变形法可自动搜寻最危险的失稳模式，模拟破坏过程，使工程措施的制定更具针对性。非连续变形法严格满足平衡要求及能量守恒，符合运动学规律，计算结果较精确，可作为制定工程措施的依据。

重力坝 深层抗滑稳定 非连续变形法 破坏过程

1 前言

重力坝的深层抗滑稳定是关系到整个大坝安危的重大问题，在水利工程设计建设过程中经常遇到，既涉及安全又影响到工程投资及建设工期。过于保守的计算方法，既增加了工程量及工程投资，又耽误整体工期；若采取偏于冒进的计算方法，可能会危及工程安全。若使工程的经济性和安全性都达到最优，就需要有一个与实际相吻合的相对精确的计算方法。

现阶段我国工程中比较常用的计算方法有刚体极限平衡法、有限元法。基于刚体极限平衡的等安全系数法(等K法)，是目前现行规范《SL319-2005混凝土重力坝设计规范》所推荐的深层稳定分析方法。等K法概念明确，有相对成熟的评价标准，是规范推荐的计算方法，恰当运用可以得到一个近似的整体安全度，对地质条件不甚复杂的简单滑动模式便于操作应用。等K法局限性在于，假定滑动面上所有块体安全系数相等，对于局部渐进破坏、各个击破的问题，计算结果偏于危险，无法获得应力分布及变位结果，不便于制定针对性的处理措施，只适用于剪切滑动失稳的问题，很多复杂失稳模式无法计算；对于因局部岩体强度不足而被压碎、拉裂所引起的破坏，难以求出符合实际的成果；计算软弱滑动面需根据经验进行选取，存在人为主观因素，可能遗漏更危险的破坏模式。有限元计算适用于连续介质的小变形工况，对不连续介质，因其网格剖分及加入节理后数值分析方面的弱势，需根据经验选取少数可能的控制性软弱面建模，可能漏掉更危险的滑动组合。而且有限元计算成果受网格形式、剖分密度、计算方法、计算程序等的影响较大，因精度影响，滑裂面上的切向力和外力难以平衡，最大误差可达30%，因而安全系数离散性大，使有限元计算成果难以形成一个统一的安全衡量标准，所以具体应用中多用于辅助评价。刚体极限平衡法、有限元法等方法均有各自的适用范围及局限性，到目前为止，针对重力坝深层抗滑稳定问题，仍没有一个尽如人意的计算分析方法，实际工程中常用各种方法进行对比、分析、论证，耗费了大量的人力、物力。

非连续变形分析方法(DDA)[1]是石根华先生提出的一种适用于大变形的计算分析方法，对于岩体等不连续介质的计算有着显著的优点。非连续变形分析在块体边界可以是不连续的，形成网格的块体可以是任意几何形状，严格满足平衡要求，满足能量守恒，可进行大变形计算，完全符合运动学规律。DDA数值可靠性高，计算效率高，所作分析非常接近实际，力学现象的数学和数值描述与块体运动相一致。对深层抗滑稳定问题来说，块体间弹簧上的力能精确满足平衡条件，所求得的力是精确的，故可以精确求出滑裂面的安全系数。DDA可分析系统内各块体的应力、应变、滑动、块体接触力和位移，可计算块体界面的张开或闭合，可进行大变形计算，故可通过增加荷载或降低软弱面强度指标的方法进行破坏过程和破坏模式的模拟。

2 非连续变形法基本原理

非连续变形分析法将自然界中完整闭合的形体均视为独立单元，在每个单元内定义一个插值函数，并由该插值函数根据最小能量原理构造出每个单元的瞬时静力平衡方程。相接触的块体和块体间通过法向和切向弹簧连接成一个系统，进而构造整体平衡方程来求解。

设对任意一块体有常应力和常应变，块体任一点(x，y)的位移(u，ν)可用六个位移不变量表示：

(u0υ0r0εxεyγxy)

(1)

其中，(u0，ν0)是块体内特殊点(x0，y0)的刚体位移；角r0是块体绕转动中心(u0，ν0)的转动角，是用弧度给出的；εx，εy，γxy是块体的切向应变。

图1 平行位移 图2 转动

图3 法向应变 图4 剪切应变 考虑块体只包含平行移动(u0，ν0)时，块体任意点(x，y)的位移(u，ν)可表示为：

(2)

对小位移，当块体只包含绕点(x0，y0)转动角r0时，块体任一点(x，y)的位移(u，ν)可表示为：

(3)

当块体只有法向应变εx、εy时，块体任一点(x，y)的位移(u，ν)可表示为：

(4)

当块体只有剪应变γxy时，块体任一点(x，y)的位移(u，ν)可表示为：

(5)

点(x，y)的总位移(u，ν)是包括所有变量(u0，ν0，γ0，εx，εy，γxy)的位移累加，则块体变形矩阵写为：

(6)

对第i块体，块体变形矩阵为：

块体变形的全一阶近似公式为：

u=a1+a2x+a3y

ν=b1+b2x+b3y

(7)

式中，(u，ν)是点(x，y)的位移。

在点(x0，y0)上位移(u0，ν0)为：

u0=a1+a2x0+a3y0

图5 一般变形

ν0=b1+b2x0+b3y0

(8)

(9)

则

(10)

即有

(11)

任意点(x，y)的位移可表示为：

u=∑mj=1ajfj(x，y)ν=∑mj=1bjfj(x，y)

(12)

各块体互相接触、互相约束，形成一个系统，整个系统的平衡方程为：

(13)

对任意块体，根据最小能量原理有：

(14)

(15)

(16)

(17)

由此，可根据最小能量原理求解系统整体平衡方程。

3 非连续变形法应用实例

武都水库工程枢纽区位于龙门山主边界断裂之F5、F7断层之间，大坝建基岩体岩性为泥盆系中统白云岩、灰岩，拦河大坝为碾压混凝土重力坝，最大坝高120.34m。河床坝段坝基内断层和缓倾角结构面发育，存在深层抗滑稳定问题。因该工程利用三维地质力学模型试验对工程安全性进行了验证，为方便对比分析，利用DDA程序对同样边界条件的典型剖面(图6)进行分析，探究破坏模式、破坏过程及坝和地基系统的稳定安全系数。

图6 武都水库典型坝段地质剖面示意

断层F31主要由糜棱岩、断层泥、角砾岩和压碎岩组成，为岩屑夹泥型，抗剪断强度指标为f′=0.45、c′=0.06MPa。断层10f2顺河方向发育，高程535m以上为岩屑岩块夹泥型，抗剪断强度指标为f′=0.5、c′=0.06MPa；535m高程以下为岩块岩屑型，抗剪断强度指标为f′=0.55、c′=0.1MPa。断层f115在高程545m以上，为泥夹岩屑型，抗剪断强度指标为f′=0.35、c′=0.03MPa；545m高程以下为岩屑夹泥型，抗剪断强度指标为f′=0.45、c′=0.06MPa。断层f114在高程545m以上为岩屑夹泥型、以下为岩块岩屑型。

在坝体上下游及底部均取大于一倍坝高基础进行建模，在基础上游侧、下游侧、底部均施加固定约束(见图7)，可利用超载或强度折减法对破坏形式和破坏过程进行模拟。本次计算为与其它计算方法有对比性，采用强度折减法进行计算。利用DDA程序可清晰获得整个滑动过程，首先从上游面沿断层F31拉裂，进而断层10f2被拉裂，随后沿断层f115向下游滑动(破坏过程见图8、图9)。以位移为零时的临界值取为安全系数，当计算采取强度折减法时，计算所得安全系数为1.46。

图7 典型坝段未处理时的DDA计算模型 图8 破坏沿F31-10f2拉裂 图9 沿断层f115向下游滑动

因该坝段，断层10f2、f115埋深较浅，断层多含泥，抗剪断指标低，深层抗滑稳定性不满足规范要求。需对坝基进行加固处理，采用坝基上下游侧均开挖大齿槽，置换为混凝土的方案。处理后剖面图及相应DDA模型分别见图10、图11。

图10 典型坝段处理后剖面示意

利用强度折减法对采取工程处理后的破坏形式和破坏过程进行模拟。当强度折减系数为3.4时，对应破坏过程及滑面组合见图12、图13。首先沿图12面拉裂，继而滑动破坏。以位移为零时的临界值取为安全系数，最危险滑面组合的安全系数为3.39，加固处理后该坝段深层抗滑稳定性满足规范要求。

该工程同时利用以下方法进行深层抗滑稳定计算：(1)基于材料力学的等K法(将建基面上正应力简化为直线分布，剪应力简化为沿建基面均匀分布，移开坝体研究基础稳定)；(2)基于有限元的等K法(利用文献〔3〕中的方法，根据有限元计算获得建基面上的应力分布，继而移开坝体研究基础的稳定)。上述三种方法计算成果对比见表1。

图11 典型坝段处理后DDA计算模型 图12 处理后破坏形式(拉裂) 图13 处理后破坏形式(滑动破坏)

表1 三种方法计算成果

同时，该工程亦进行了地质力学模型试验，基本结论如下：地基未处理时采用超载破坏试验法，当超载系数为1.6时发生初裂，超载安全度为2.2；经置换处理后，采用以超载为主、降强为辅的方法进行破坏试验，断层的抗剪断强度降低15%～20%，再超载至3.0P0～4.0P0时，变位幅度明显增大，出现破坏失稳趋势。综合稳定安全度Kc=3.6～4.8。

综上可见，非连续变形法计算成果与地质力学模型试验结论、基于材料力学的等K法、基于有限元的等K法在工程定性上有良好的一致性，可作为对深层抗滑稳定问题进行定性和定量评价、制定处理方案的依据。同时，非连续变形法可自动判定可能的滑动组合，并模拟其破坏过程，使工程措施的制定更具备针对性。

4 结语

非连续变形法将每一个完整的块体作为一个单元，将节理、裂隙等软弱面及不同材料分区界面作为块体的边界，单元可以是任意形状的，允许不连续，适用性广，可用于任何地质情况；可精确得到任意一点的应力、变位；可自动搜寻出最危险的失稳模式、模拟破坏过程。

采用多种计算方法计算之前，首先需进行经验判断，选取可能的危险模式，这种情况下就有部分主观因素加入，因工程人员的经验不同而有可能漏判。非连续变形法可自行判定可能的滑动组合，并给出破坏、滑动过程，可发现因人力限制而通过主观判断想象不到的破坏模式，对工程安全起着重要作用。DDA在模拟破坏过程及应力精确性方面有优势，建议对可能存在深层稳定问题的工程，首先利用DDA程序进行破坏分析，观察破坏的渐进过程，在此基础上采用其他计算方法抽样验证。

〔1〕石根华.数值流形方法与非连续变形分析.北京：清华大学出版社，1997.

〔2〕张国新，刘 毅.坝基稳定分析的有限元直接反力法，水力发电，2007，(1).

〔3〕Sarma.S.K.Stabilityanalysisofembankmentsandslopes[J].Geotechnique.1973，(23)：423～433.

〔4〕张津生.浅析刚体极限平衡法——探讨重力坝深层抗滑稳定的安全判据.水力发电学报，2005，24(5).

〔5〕黄东军，聂广明.重力坝深层抗滑稳定安全评价若干问题的思考.水利发电学报，2005，24(2)：90～94.

王 静(1980-)，女，山东滕州人，高级工程师，主要从事水利水电工程设计工作；

覃克非(1967-)，男，成都人，教授级高工，主要从事水利水电工程设计及管理工作；

王 敏(1976-)，女，山东滕州人，讲师，主要从事计算分析及计算机教学工作。

■

TV642.3∶TV

B

2095-1809(2016)01-0040-05