岩体结构变形破坏的内在驱动力-不平衡力*

杨 强 王守光 李超毅 刘耀儒

(清华大学水沙科学与水利水电工程国家重点实验室， 北京 100084，中国)

0 引 言

岩体结构的破坏是一个渐进的破坏过程(胡启军等， 2011； 李世贵等， 2018)。规范采用的强度与极限设计主要针对的是弹性极限状态和极限状态，这两个特征状态之间的变形破坏过程并未涉及。但实际上在变形破坏过程中，出现开裂和显著的非线性变形也是工程设计极为关注的状态。法国的Malpasset拱坝溃坝就是拱坝坝踵开裂后，坝基内应力渗流耦合效应所致(Serafim， 1987)。非线性变形控制也是岩体工程重要控制指标，洞室大变形多和围岩内开裂损伤有关，拱坝的显著非线性变形与坝趾的损伤破坏有关(程立， 2017)。

一般的开裂损伤分析计算方法难以有效应用于岩体结构，原因在于岩体自身具有的节理裂隙带来的复杂不连续性。常规的连续介质计算方法(有限元法、有限差分法等)假设岩体结构连续，而把裂纹当作计算边界，这样做的缺点是每当裂纹扩展就要重新划分网格边界。扩展有限元法(Belytschko et al.， 1999； 李录贤等， 2005)虽然解决了重新剖分网格的问题，但也难以分析复杂不连续体的多裂纹扩展贯通。为了克服连续性理论基础带来的种种约束，Silling(2000)放弃了连续性假设转而把固体离散为一系列物质点，通过求解空间积分方程来描述材料力学行为，这样就从根本上解决了不连续导致的奇异性问题。但这种方法目前仅能处理简单三维结构的开裂问题。更深入的开裂分析需要考虑结构的损伤和动态演化过程(Yang et al.， 2005； 孙琪皓等， 2019)，甚至还需考虑“远距离物质点的运动历史”影响，即非局部理论(王林娟等， 2019)。上述方法虽然在理论上更加完备，但对于岩体结构这种三维复杂结构，计算效率和实用性目前尚无法保证。

从岩体工程的实践来看，要准确获取岩体的弹性模量E、泊松比v、 摩擦系数f、黏聚力c等参数已属不易，而这只能满足理想弹塑性分析的要求，远远无法满足先进的开裂损伤分析的要求。从岩体工程设计的角度来说，一般也不追求精确的破坏过程和最终破坏模式，设计重点放在加固设计上，以加固措施抑制和控制开裂和变形破坏的发生和发展。有鉴于此，准确把握岩体结构开裂和变形破坏的内在驱动力对指导岩体工程设计具有重大意义，可大大提高加固设计的有效性和针对性。无疑这对控制岩爆等特殊破坏模式也有重要意义。

事实上，寻求岩体结构变形破坏的内在驱动力一直是岩石力学界的长期关切。如王思敬(2002)认为“内外动力耦合作用是重大地质灾害的成因”； 何满潮(2016)认为“牛顿力突变，灾害发生”。

杨强等(2004，2008)， Yang et al.(2008，2012)提出了一个全新的思路，即结构的开裂破坏可转换为结构弹塑性有限元分析解的存在性问题。基于位移法的有限元分析始终要求位移场保持连续，由此排除了开裂出现的可能性，因为开裂意味着出现不连续位移场。在位移场保持连续的前提下，开裂破坏在一定程度上可以用材料损伤来近似描述，如果采用无损本构模型，则排除了出现开裂破坏。理想弹塑性模型是广泛采用的无损本构模型。所以结构理想弹塑性有限元分析的解如果存在，该结构一定是连续且无损伤； 反之如果解不存在，结构必然出现损伤开裂。

结构弹塑性有限元分析解不存在表现为迭代计算不收敛，存在无法消除的残余不平衡力。该不平衡力即可视为岩体结构开裂破坏的内在驱动力。为抑制结构开裂破坏，结构所需加固力必定与不平衡力大小相等方向相反，这就是所谓的变形加固理论。

这种基于不平衡力的开裂破坏分析方法的好处是回避了直接模拟开裂破坏过程所带来计算上的巨大的复杂性，且理想弹塑性模型和岩体工程实践基本相适应，包括岩体参数取值与加固设计。

一般说来，计算不收敛或者说出现无法消除的残余不平衡力，这违背了经典弹塑性理论及其有限元分析，计算成果是无意义的。为此本文对弹塑性结构分析给出更一般性的提法：在保持结构连续性及材料无损的前提下，结构作用力与抗力的差值必然趋向最小值，该差值就是残余不平衡力。这种提法包容了出现不平衡力情况，而常规结构弹塑性分析的解就是该差值为0的一个特例。

本文先简要介绍不平衡力的概念和求解方法，然后从几个数值和试验算例验证不平衡力分布与岩体结构变形破坏的相关性，接着从蓄水诱发岩体结构非平衡演化的角度推广不平衡力的概念。最后结合弹塑性结构分析一般性的提法，对不平衡力的理论基础和实质进行说明。

1 弹塑性有限元计算中的不平衡力

对一个结构而言，与外荷载平衡的应力场σ1，对应于作用力，可称为作用应力场：

(1)

其中，B为应变矩阵； 下标e表示对所有单元求和；F为外荷载节点力向量。

满足屈服条件的应力场σ可视为结构抗力，可称为抗力应力场：

f(σ)≤0

(2)

f为屈服函数。如图 1所示，两者之间的差值即为塑性应力Δσp：

Δσp=σ1-σ

(3)

塑性应力的等效节点力ΔU即为不平衡力：

(4)

σ1由弹塑性迭代过程中的应变增量Δε确定：

σ1=σ0+D︰Δε

(5)

其中，D为弹性张量。如果在结构某一点上，f(σ1)≤0，则有σ=σ1。如果f(σ1)>0，杨强等(2004)指出正交流动法则要求σ与σ1在塑性余能(PCE)意义下距离最近：

(6)

其中，C为柔度张量。式(6)反映了在一个材料点上作用力与抗力差距最小的要求。

图 1 弹塑性应力调整示意图Fig. 1 Diagram of elastic-plastic stress adjustment

对于理想弹塑性材料，当屈服条件f采用Drucker-Prager屈服准则时，σ的解析解如式(7)(杨强等， 2002)：

(7)

式中，I1，J2为应力偏量；α为材料常数；K，G为体积模量和剪切模量。

2 数值与试验算例

2.1 不平衡力与开裂相关性

破裂将导致不连续的位移场，而不连续位移场会导致弹塑性边值问题解不存在，在弹塑性有限元分析中，无解就是计算不收敛，即存在无法消除的不平衡力。所以，由式(1)～式(7)的推导，不平衡力必然应该是开裂状态的一种表征。

为了验证不平衡力与开裂在理论上存在的相关性，我们对一些已有的开裂破坏试验进行数值模拟分析。数值仿真计算在课题组自行编写的三维非线性有限元程序TFINE软件中进行。图 2是陈新等(2014)做的含预置裂纹石膏试件单轴压缩试验。石膏试件的尺寸为： 15icm×15icm×5icm，预制裂纹倾角为30°，预置裂纹平行排开，间距为30imm，试件的最终破坏如图 2所示。对图 2的破坏过程进行数值仿真，计算结果如图 3所示，图 3中红色箭头为不平衡力矢量。由图 3可以看出，不平衡力基本都出现在裂纹扩展的位置，显示出了它们之间很好的相关性。

图 2 含预置裂纹石膏试件破坏过程裂缝分布图 (陈新等， 2014)Fig. 2 Cracking maps of gypsum specimens with preset cracks (Chen et al.，2014)a. 加载过程1； b. 加载过程2； c. 加载过程3； d. 加载过程4

图 3 基于不平衡力的石膏开裂过程数值模拟Fig. 3 Numerical simulation of gypsum cracking based on unbalanced forcea. 加载过程1； b. 加载过程2； c. 加载过程3； d. 加载过程4

2.2 基于不平衡力的白鹤滩边坡开挖卸荷分析

拱坝建基面边坡开挖卸荷将会打破边坡原有的平衡状态，进入非平衡状态。在一些具有复杂地应力条件的坝址区，如何采取合理的施工方案以避免严重的松弛失稳现象出现一直是水利工程施工关注的重点问题之一。程立等(2017)基于不平衡力对白鹤滩拱坝左岸建基面边坡开挖卸荷进行了系统研究。白鹤滩拱坝左岸建基面在开挖到628im高程时出现了卸荷松弛现象，为了弄清楚坝基继续开挖卸荷松弛的演变规律，程立等(2017)使用不平衡力作为坝基岩体卸荷松弛的定量判据(图 4)。

图为由590im高程向下开挖至538im高程过程中，拱坝建基面某典型剖面的不平衡力矢量图。由图 4可知白鹤滩左岸建基面陡坎成型后，陡坎的底部和建基面坡脚附近出现明显的不平衡力集中(图 4b)，在陡坎底部不平衡力尤为突出。因此，建议将左岸建基面570～590im高程之间的陡坎坡度放缓，并加强对陡坎底部断层LS331出露段附近的监测与研究。

图 4 白鹤滩拱坝左岸建基面继续开挖至538im高程过程中 某典型剖面的不平衡力矢量图(程立等， 2017)Fig. 4 Unbalanced force vector map of a typical section during the excavation to 538im elevation of left slope of Baihetan arch dam(Cheng et al.，2017)a. 开挖至590im 高程； b. 开挖至560im 高程； c. 开挖至550im 高程； d. 开挖至538im 高程

2.3 JH二级拱坝坝踵和坝趾开裂破坏分析

拱坝开裂分析在理论和应用上还存在诸多困难。由上述分析表明，不平衡力与裂纹位置和形貌有很好的相关性。因此不平衡力可以提供一种间接预测拱坝开裂的方法，确定可能出现的开裂位置并找到所有裂缝中最危险的裂缝。本节基于不平衡力研究JH二级拱坝的开裂破坏，并与模型试验结果对比研究。

图 5 JH二级拱坝数值模型图Fig. 5 Numerical model of JH-2iarch dam

图 6 模拟的断层位置分布图Fig. 6 Distribution map of faults

图 7 JH二级拱坝地质力学模型试验图Fig. 7 Geomechanical model test of JH-2iarch dama. JH二级模型试验上游坝面图； b. JH二级模型试验下游坝面图； c. 模型试验加载装置图

JH二级双曲拱坝坝高167.5im。总体来说，河谷较窄，坝体较为厚实。JH二级拱坝数值模型如图 5所示，模型网格节点总数132i583，单元总数为120i048； 对于6条主要断层进行了模拟：断层jef4、jef36、jef61、jef8、jef51和jef38，如图 6所示。同时，为了验证数值模拟的结果，对JH二级拱坝进行地质力学模型试验研究，模型试验图如图 7所示。

图 8 拱坝上游面坝踵开裂图Fig. 8 Cracking map of upstream face and heel of arch dam

图 9 拱坝下游面坝趾破坏图Fig. 9 Damage map of dam toe at downstream face of arch dam

图 8a是拱坝上游面的不平衡力计算结果图，图 8b是拱坝上游面模型试验破坏图，图 9a是拱坝下游面的塑性区计算结果图，图 9b是拱坝下游面模型试验破坏图。由图 8a和图8b可知，当加载到两倍水载时，不平衡力最先出现在左右坝踵与河床相接的位置，而模型试验结果显示，最先起裂的地方正是在这个位置并向河床延伸。由图 9a可知，当加载到5倍水载时，下游坝面坝趾区塑性区与上游贯通，坝趾应该在此时破坏。模型试验的结果也恰恰印证了这一点(图 9b)。

3 蓄水诱发的不可逆谷幅变形-不平衡力作用机制

拱坝开挖卸荷、超载破坏是外荷载超出屈服面产生不平衡力。而另一种形式的不平衡力则是外荷载不变，屈服面收缩产生。近年来，蓄水诱发的谷幅变形就是后者的一个体现(杨强等， 2015； 周志芳等， 2019)。高拱坝蓄水期，裂隙的渗透性远高于完整岩块，库水首先由裂隙快速入渗进入山体，使地下水位在短期内就达到与库水位大致相同的水平。

裂隙水压力升高后，水推力或扬压力都使坡体的不平衡力增大。裂隙水压力使岩体内部压应力减小，拉应力增大。对于弹塑性本构的连续介质体而言，效果是使材料进入塑性屈服。因此，杨强等(2015)认为高拱坝蓄水后，裂隙水压力将使岩土材料屈服面收缩，不平衡力增大，从而产生不可逆的塑性变形：

f(σ′)=f(σ-pI) ∀p

(8)

其中，σ′为有效应力张量；p为水压力；I为单位张量。如式(8)所示，在屈服函数中总应力应替换为有效应力，即考虑了裂隙水压力后，应力空间里的屈服面收缩，在π平面上的许可应力范围半径减小，如图 10b所示。因此，裂隙水的介入使原先稳定状态处于临界状态(图 10a)、原先处于屈服或临界屈服状态的岩体应力状态超出屈服面，发生进一步的塑性变形(图 1)，直到达到新的平衡，或者发生结构失稳。

锦屏一级拱坝蓄水诱发的谷幅收缩变形特征与上述思想不谋而合：蓄水增大了边坡裂隙水压力，从而使边坡朝临空面产生不可逆的塑性变形，形成了谷幅收缩现象(图 11)。事实上，杨强等(2015)也正是基于上述思想解释了锦屏一级拱坝蓄水诱发谷幅收缩现象的。

图 10 考虑裂隙水压力的屈服准则 (杨强等， 2015)Fig. 10 Yield criterion considering pore water pressure (Yang et al.，2015)a. Mohr-Coulomb准则； b. Drucker-Prager准则

图 11 锦屏一级谷幅变形随库水位变化曲线 变形负向增加表示谷幅收缩加大Fig. 11 Curve of valley width reduction with water level of Jinping Ⅰ arch dam. Negative increase of deformation indicates increase of valley width reduction

4 不平衡力实质的讨论

变形与破坏是结构非平衡的外在表现和必然结果，本章从结构非平衡出发，讨论不平衡力的实质。

4.1 单滑块体的不平衡力

结构非平衡在一维条件下对应单滑块失稳模型，如图 12a和图12b所示。对于稳定状态：滑动力T=阻滑力R≤fN+cA； 对于失稳状态：滑动力T>阻滑力R=fN+cA。

图 12 一维单滑块体不平衡力模型Fig. 12 One dimensional unbalanced force model of single slidera. 稳定状态；b. 失稳及考虑加固力的稳定状态

滑块失稳本质就是滑动力超过了最大阻滑力，平衡条件T=R被破坏，导致不平衡力U产生对稳定状态U=0； 对于失稳状态，阻滑力发挥到最大值，使得不平衡力最小化。由此我们可以看出稳定和失稳状态的共性之处，就是阻滑力在其许可范围内R≤fN+cA，其取值必使不平衡力最小化。加固的本质是在加固力Q的帮助下，滑块恢复平衡状态，T=R+Q。所以不平衡力和加固力大小相等，方向相反。

U=T-R

(9)

如果将滑动力和阻滑力分别视为作用力和抗力，则前述弹塑性结构分析一般性提法显然也适用于稳定或失稳状态下的单滑块模型：结构作用力与抗力的差值必然趋向最小值，该差值就是残余不平衡力。

4.2 材料的不平衡力

由式(6)可知，对于一给定的作用应力和屈服条件，抗力应力必使式(10)所示的余能范数最小，即minE(σ)：

(10)

式(10)是一个条件极值问题，这个问题的解对应着弹塑性增量型的正交流动法则和一致性条件。为解出minE(σ)，构造Lagrange函数：

L=E(σ)-Δλf(σ)

(11)

则应有：

(12)

把式(10)和式(11)代入式(12)解得：

(13)

式(13)就是弹塑性增量型的正交流动法则和一致性条件。由此可知作用力与抗力差值趋向最小值是弹塑性本构关系的要求。

4.3 结构的不平衡力

受式(10)的启发，可以推测结构作用应力场与抗力应力场，在总塑性余能意义下差距最小，即最小塑性余能原理：

(14)

式(14)是在整个结构的积分。塑性余能是不平衡力的范数，所以最小塑性余能原理意味着结构趋于不平衡力最小的状态。

事实上，对于式(14)的推测可以在理论上给予严格的证明。Yang et al.(2012)从弹-黏塑性结构出发，采用Duvaut-Lions模型证明了三维结构非平衡状态最终演化趋势满足最小塑性余能原理。因此，不平衡力无论在材料层次还是结构层次都有严格的理论基础。

图 13 外力场和自承力场动态调整图Fig. 13 Dynamic adjustment diagram of external and self-supporting force fields

有了最小塑性余能原理这一指导性原则，不平衡力的求解就成了一个双场优化的迭代过程：通过变形调整达到σ1和σ最接近的结构整体塑性余能最小状态(图 13)。

潘家铮(1980)最大最小原理可由最小塑性余能予以说明：(1)最小值原理：滑坡如能沿许多滑面滑动，则失稳时，它将沿抵抗力最小的一个滑面破坏——加固力最小。(2)最大值原理：滑坡体的滑面肯定时，则滑面上的反力(以及滑坡体内的内力)能自行调整，以发挥最大的抗滑能力——自承力最大。最小塑性余能原理反映了多自由度超静定体系的非线性内力分配原则； 而潘家铮最大最小原理用于确定在极限状态下多滑块体系的内力分配。

从计算策略上来讲，非线性有限元本质是时空离散化。空间离散化是指有限单元在空间上近似无限自由度变形体，等效基础是最小势能原理； 而时间离散化是指无限小微分加载过程的离散化或增量化，等效基础是最小塑性余能原理。最小塑性余能原理突破了经典弹塑性理论，允许应力超出屈服面。

5 结 论

基于理想弹塑性有限元分析，岩体结构无法消除的残余不平衡力就是其变形破坏的内在驱动力。研究案例表明不平衡力分布与结构开裂破坏关系密切。

本文对弹塑性结构分析给出更一般性的提法：在保持结构连续性及材料无损的前提下，结构作用力与抗力的差值必然趋向最小值(体现为最小塑性余能原理)，该差值就是残余不平衡力。这种提法包容了出现不平衡力情况，而常规结构弹塑性分析的解就是该差值为0的一个特例。最小塑性余能原理可视为弹塑性理论的增量化的要求。施加反向残余不平衡力作为加固力，理论上可完全避免结构开裂损伤，保持结构连续性。

需要说明的是，和单滑块模型不同，变形体结构出现不平衡力并不意味着结构要整体失稳，更多是出现结构局部损伤破坏。