巖體結(jié)構(gòu)變形破壞的內(nèi)在驅(qū)動力-不平衡力*

楊 強 王守光 李超毅 劉耀儒

(清華大學(xué)水沙科學(xué)與水利水電工程國家重點實驗室， 北京 100084，中國)

0 引 言

巖體結(jié)構(gòu)的破壞是一個漸進的破壞過程(胡啟軍等， 2011； 李世貴等， 2018)。規(guī)范采用的強度與極限設(shè)計主要針對的是彈性極限狀態(tài)和極限狀態(tài)，這兩個特征狀態(tài)之間的變形破壞過程并未涉及。但實際上在變形破壞過程中，出現(xiàn)開裂和顯著的非線性變形也是工程設(shè)計極為關(guān)注的狀態(tài)。法國的Malpasset拱壩潰壩就是拱壩壩踵開裂后，壩基內(nèi)應(yīng)力滲流耦合效應(yīng)所致(Serafim， 1987)。非線性變形控制也是巖體工程重要控制指標(biāo)，洞室大變形多和圍巖內(nèi)開裂損傷有關(guān)，拱壩的顯著非線性變形與壩趾的損傷破壞有關(guān)(程立， 2017)。

一般的開裂損傷分析計算方法難以有效應(yīng)用于巖體結(jié)構(gòu)，原因在于巖體自身具有的節(jié)理裂隙帶來的復(fù)雜不連續(xù)性。常規(guī)的連續(xù)介質(zhì)計算方法(有限元法、有限差分法等)假設(shè)巖體結(jié)構(gòu)連續(xù)，而把裂紋當(dāng)作計算邊界，這樣做的缺點是每當(dāng)裂紋擴展就要重新劃分網(wǎng)格邊界。擴展有限元法(Belytschko et al.， 1999； 李錄賢等， 2005)雖然解決了重新剖分網(wǎng)格的問題，但也難以分析復(fù)雜不連續(xù)體的多裂紋擴展貫通。為了克服連續(xù)性理論基礎(chǔ)帶來的種種約束，Silling(2000)放棄了連續(xù)性假設(shè)轉(zhuǎn)而把固體離散為一系列物質(zhì)點，通過求解空間積分方程來描述材料力學(xué)行為，這樣就從根本上解決了不連續(xù)導(dǎo)致的奇異性問題。但這種方法目前僅能處理簡單三維結(jié)構(gòu)的開裂問題。更深入的開裂分析需要考慮結(jié)構(gòu)的損傷和動態(tài)演化過程(Yang et al.， 2005； 孫琪皓等， 2019)，甚至還需考慮“遠距離物質(zhì)點的運動歷史”影響，即非局部理論(王林娟等， 2019)。上述方法雖然在理論上更加完備，但對于巖體結(jié)構(gòu)這種三維復(fù)雜結(jié)構(gòu)，計算效率和實用性目前尚無法保證。

從巖體工程的實踐來看，要準(zhǔn)確獲取巖體的彈性模量E、泊松比v、 摩擦系數(shù)f、黏聚力c等參數(shù)已屬不易，而這只能滿足理想彈塑性分析的要求，遠遠無法滿足先進的開裂損傷分析的要求。從巖體工程設(shè)計的角度來說，一般也不追求精確的破壞過程和最終破壞模式，設(shè)計重點放在加固設(shè)計上，以加固措施抑制和控制開裂和變形破壞的發(fā)生和發(fā)展。有鑒于此，準(zhǔn)確把握巖體結(jié)構(gòu)開裂和變形破壞的內(nèi)在驅(qū)動力對指導(dǎo)巖體工程設(shè)計具有重大意義，可大大提高加固設(shè)計的有效性和針對性。無疑這對控制巖爆等特殊破壞模式也有重要意義。

事實上，尋求巖體結(jié)構(gòu)變形破壞的內(nèi)在驅(qū)動力一直是巖石力學(xué)界的長期關(guān)切。如王思敬(2002)認(rèn)為“內(nèi)外動力耦合作用是重大地質(zhì)災(zāi)害的成因”； 何滿潮(2016)認(rèn)為“牛頓力突變，災(zāi)害發(fā)生”。

楊強等(2004，2008)， Yang et al.(2008，2012)提出了一個全新的思路，即結(jié)構(gòu)的開裂破壞可轉(zhuǎn)換為結(jié)構(gòu)彈塑性有限元分析解的存在性問題?；谖灰品ǖ挠邢拊治鍪冀K要求位移場保持連續(xù)，由此排除了開裂出現(xiàn)的可能性，因為開裂意味著出現(xiàn)不連續(xù)位移場。在位移場保持連續(xù)的前提下，開裂破壞在一定程度上可以用材料損傷來近似描述，如果采用無損本構(gòu)模型，則排除了出現(xiàn)開裂破壞。理想彈塑性模型是廣泛采用的無損本構(gòu)模型。所以結(jié)構(gòu)理想彈塑性有限元分析的解如果存在，該結(jié)構(gòu)一定是連續(xù)且無損傷； 反之如果解不存在，結(jié)構(gòu)必然出現(xiàn)損傷開裂。

結(jié)構(gòu)彈塑性有限元分析解不存在表現(xiàn)為迭代計算不收斂，存在無法消除的殘余不平衡力。該不平衡力即可視為巖體結(jié)構(gòu)開裂破壞的內(nèi)在驅(qū)動力。為抑制結(jié)構(gòu)開裂破壞，結(jié)構(gòu)所需加固力必定與不平衡力大小相等方向相反，這就是所謂的變形加固理論。

這種基于不平衡力的開裂破壞分析方法的好處是回避了直接模擬開裂破壞過程所帶來計算上的巨大的復(fù)雜性，且理想彈塑性模型和巖體工程實踐基本相適應(yīng)，包括巖體參數(shù)取值與加固設(shè)計。

一般說來，計算不收斂或者說出現(xiàn)無法消除的殘余不平衡力，這違背了經(jīng)典彈塑性理論及其有限元分析，計算成果是無意義的。為此本文對彈塑性結(jié)構(gòu)分析給出更一般性的提法：在保持結(jié)構(gòu)連續(xù)性及材料無損的前提下，結(jié)構(gòu)作用力與抗力的差值必然趨向最小值，該差值就是殘余不平衡力。這種提法包容了出現(xiàn)不平衡力情況，而常規(guī)結(jié)構(gòu)彈塑性分析的解就是該差值為0的一個特例。

本文先簡要介紹不平衡力的概念和求解方法，然后從幾個數(shù)值和試驗算例驗證不平衡力分布與巖體結(jié)構(gòu)變形破壞的相關(guān)性，接著從蓄水誘發(fā)巖體結(jié)構(gòu)非平衡演化的角度推廣不平衡力的概念。最后結(jié)合彈塑性結(jié)構(gòu)分析一般性的提法，對不平衡力的理論基礎(chǔ)和實質(zhì)進行說明。

1 彈塑性有限元計算中的不平衡力

對一個結(jié)構(gòu)而言，與外荷載平衡的應(yīng)力場σ1，對應(yīng)于作用力，可稱為作用應(yīng)力場：

(1)

其中，B為應(yīng)變矩陣； 下標(biāo)e表示對所有單元求和；F為外荷載節(jié)點力向量。

滿足屈服條件的應(yīng)力場σ可視為結(jié)構(gòu)抗力，可稱為抗力應(yīng)力場：

f(σ)≤0

(2)

f為屈服函數(shù)。如圖 1所示，兩者之間的差值即為塑性應(yīng)力Δσp：

Δσp=σ1-σ

(3)

塑性應(yīng)力的等效節(jié)點力ΔU即為不平衡力：

(4)

σ1由彈塑性迭代過程中的應(yīng)變增量Δε確定：

σ1=σ0+D︰Δε

(5)

其中，D為彈性張量。如果在結(jié)構(gòu)某一點上，f(σ1)≤0，則有σ=σ1。如果f(σ1)>0，楊強等(2004)指出正交流動法則要求σ與σ1在塑性余能(PCE)意義下距離最近：

(6)

其中，C為柔度張量。式(6)反映了在一個材料點上作用力與抗力差距最小的要求。

圖 1 彈塑性應(yīng)力調(diào)整示意圖Fig. 1 Diagram of elastic-plastic stress adjustment

對于理想彈塑性材料，當(dāng)屈服條件f采用Drucker-Prager屈服準(zhǔn)則時，σ的解析解如式(7)(楊強等， 2002)：

(7)

式中，I1，J2為應(yīng)力偏量；α為材料常數(shù)；K，G為體積模量和剪切模量。

2 數(shù)值與試驗算例

2.1 不平衡力與開裂相關(guān)性

破裂將導(dǎo)致不連續(xù)的位移場，而不連續(xù)位移場會導(dǎo)致彈塑性邊值問題解不存在，在彈塑性有限元分析中，無解就是計算不收斂，即存在無法消除的不平衡力。所以，由式(1)～式(7)的推導(dǎo)，不平衡力必然應(yīng)該是開裂狀態(tài)的一種表征。

為了驗證不平衡力與開裂在理論上存在的相關(guān)性，我們對一些已有的開裂破壞試驗進行數(shù)值模擬分析。數(shù)值仿真計算在課題組自行編寫的三維非線性有限元程序TFINE軟件中進行。圖 2是陳新等(2014)做的含預(yù)置裂紋石膏試件單軸壓縮試驗。石膏試件的尺寸為： 15icm×15icm×5icm，預(yù)制裂紋傾角為30°，預(yù)置裂紋平行排開，間距為30imm，試件的最終破壞如圖 2所示。對圖 2的破壞過程進行數(shù)值仿真，計算結(jié)果如圖 3所示，圖 3中紅色箭頭為不平衡力矢量。由圖 3可以看出，不平衡力基本都出現(xiàn)在裂紋擴展的位置，顯示出了它們之間很好的相關(guān)性。

圖 2 含預(yù)置裂紋石膏試件破壞過程裂縫分布圖 (陳新等， 2014)Fig. 2 Cracking maps of gypsum specimens with preset cracks (Chen et al.，2014)a. 加載過程1； b. 加載過程2； c. 加載過程3； d. 加載過程4

圖 3 基于不平衡力的石膏開裂過程數(shù)值模擬Fig. 3 Numerical simulation of gypsum cracking based on unbalanced forcea. 加載過程1； b. 加載過程2； c. 加載過程3； d. 加載過程4

2.2 基于不平衡力的白鶴灘邊坡開挖卸荷分析

拱壩建基面邊坡開挖卸荷將會打破邊坡原有的平衡狀態(tài)，進入非平衡狀態(tài)。在一些具有復(fù)雜地應(yīng)力條件的壩址區(qū)，如何采取合理的施工方案以避免嚴(yán)重的松弛失穩(wěn)現(xiàn)象出現(xiàn)一直是水利工程施工關(guān)注的重點問題之一。程立等(2017)基于不平衡力對白鶴灘拱壩左岸建基面邊坡開挖卸荷進行了系統(tǒng)研究。白鶴灘拱壩左岸建基面在開挖到628im高程時出現(xiàn)了卸荷松弛現(xiàn)象，為了弄清楚壩基繼續(xù)開挖卸荷松弛的演變規(guī)律，程立等(2017)使用不平衡力作為壩基巖體卸荷松弛的定量判據(jù)(圖 4)。

圖為由590im高程向下開挖至538im高程過程中，拱壩建基面某典型剖面的不平衡力矢量圖。由圖 4可知白鶴灘左岸建基面陡坎成型后，陡坎的底部和建基面坡腳附近出現(xiàn)明顯的不平衡力集中(圖 4b)，在陡坎底部不平衡力尤為突出。因此，建議將左岸建基面570～590im高程之間的陡坎坡度放緩，并加強對陡坎底部斷層LS331出露段附近的監(jiān)測與研究。

圖 4 白鶴灘拱壩左岸建基面繼續(xù)開挖至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二級拱壩壩踵和壩趾開裂破壞分析

拱壩開裂分析在理論和應(yīng)用上還存在諸多困難。由上述分析表明，不平衡力與裂紋位置和形貌有很好的相關(guān)性。因此不平衡力可以提供一種間接預(yù)測拱壩開裂的方法，確定可能出現(xiàn)的開裂位置并找到所有裂縫中最危險的裂縫。本節(jié)基于不平衡力研究JH二級拱壩的開裂破壞，并與模型試驗結(jié)果對比研究。

圖 5 JH二級拱壩數(shù)值模型圖Fig. 5 Numerical model of JH-2iarch dam

圖 6 模擬的斷層位置分布圖Fig. 6 Distribution map of faults

圖 7 JH二級拱壩地質(zhì)力學(xué)模型試驗圖Fig. 7 Geomechanical model test of JH-2iarch dama. JH二級模型試驗上游壩面圖； b. JH二級模型試驗下游壩面圖； c. 模型試驗加載裝置圖

JH二級雙曲拱壩壩高167.5im?？傮w來說，河谷較窄，壩體較為厚實。JH二級拱壩數(shù)值模型如圖 5所示，模型網(wǎng)格節(jié)點總數(shù)132i583，單元總數(shù)為120i048； 對于6條主要斷層進行了模擬：斷層jef4、jef36、jef61、jef8、jef51和jef38，如圖 6所示。同時，為了驗證數(shù)值模擬的結(jié)果，對JH二級拱壩進行地質(zhì)力學(xué)模型試驗研究，模型試驗圖如圖 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是拱壩上游面的不平衡力計算結(jié)果圖，圖 8b是拱壩上游面模型試驗破壞圖，圖 9a是拱壩下游面的塑性區(qū)計算結(jié)果圖，圖 9b是拱壩下游面模型試驗破壞圖。由圖 8a和圖8b可知，當(dāng)加載到兩倍水載時，不平衡力最先出現(xiàn)在左右壩踵與河床相接的位置，而模型試驗結(jié)果顯示，最先起裂的地方正是在這個位置并向河床延伸。由圖 9a可知，當(dāng)加載到5倍水載時，下游壩面壩趾區(qū)塑性區(qū)與上游貫通，壩趾應(yīng)該在此時破壞。模型試驗的結(jié)果也恰恰印證了這一點(圖 9b)。

3 蓄水誘發(fā)的不可逆谷幅變形-不平衡力作用機制

拱壩開挖卸荷、超載破壞是外荷載超出屈服面產(chǎn)生不平衡力。而另一種形式的不平衡力則是外荷載不變，屈服面收縮產(chǎn)生。近年來，蓄水誘發(fā)的谷幅變形就是后者的一個體現(xiàn)(楊強等， 2015； 周志芳等， 2019)。高拱壩蓄水期，裂隙的滲透性遠高于完整巖塊，庫水首先由裂隙快速入滲進入山體，使地下水位在短期內(nèi)就達到與庫水位大致相同的水平。

裂隙水壓力升高后，水推力或揚壓力都使坡體的不平衡力增大。裂隙水壓力使巖體內(nèi)部壓應(yīng)力減小，拉應(yīng)力增大。對于彈塑性本構(gòu)的連續(xù)介質(zhì)體而言，效果是使材料進入塑性屈服。因此，楊強等(2015)認(rèn)為高拱壩蓄水后，裂隙水壓力將使巖土材料屈服面收縮，不平衡力增大，從而產(chǎn)生不可逆的塑性變形：

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

(8)

其中，σ′為有效應(yīng)力張量；p為水壓力；I為單位張量。如式(8)所示，在屈服函數(shù)中總應(yīng)力應(yīng)替換為有效應(yīng)力，即考慮了裂隙水壓力后，應(yīng)力空間里的屈服面收縮，在π平面上的許可應(yīng)力范圍半徑減小，如圖 10b所示。因此，裂隙水的介入使原先穩(wěn)定狀態(tài)處于臨界狀態(tài)(圖 10a)、原先處于屈服或臨界屈服狀態(tài)的巖體應(yīng)力狀態(tài)超出屈服面，發(fā)生進一步的塑性變形(圖 1)，直到達到新的平衡，或者發(fā)生結(jié)構(gòu)失穩(wěn)。

錦屏一級拱壩蓄水誘發(fā)的谷幅收縮變形特征與上述思想不謀而合：蓄水增大了邊坡裂隙水壓力，從而使邊坡朝臨空面產(chǎn)生不可逆的塑性變形，形成了谷幅收縮現(xiàn)象(圖 11)。事實上，楊強等(2015)也正是基于上述思想解釋了錦屏一級拱壩蓄水誘發(fā)谷幅收縮現(xiàn)象的。

圖 10 考慮裂隙水壓力的屈服準(zhǔn)則 (楊強等， 2015)Fig. 10 Yield criterion considering pore water pressure (Yang et al.，2015)a. Mohr-Coulomb準(zhǔn)則； b. Drucker-Prager準(zhǔn)則

圖 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 不平衡力實質(zhì)的討論

變形與破壞是結(jié)構(gòu)非平衡的外在表現(xiàn)和必然結(jié)果，本章從結(jié)構(gòu)非平衡出發(fā)，討論不平衡力的實質(zhì)。

4.1 單滑塊體的不平衡力

結(jié)構(gòu)非平衡在一維條件下對應(yīng)單滑塊失穩(wěn)模型，如圖 12a和圖12b所示。對于穩(wěn)定狀態(tài)：滑動力T=阻滑力R≤fN+cA； 對于失穩(wěn)狀態(tài)：滑動力T>阻滑力R=fN+cA。

圖 12 一維單滑塊體不平衡力模型Fig. 12 One dimensional unbalanced force model of single slidera. 穩(wěn)定狀態(tài)；b. 失穩(wěn)及考慮加固力的穩(wěn)定狀態(tài)

滑塊失穩(wěn)本質(zhì)就是滑動力超過了最大阻滑力，平衡條件T=R被破壞，導(dǎo)致不平衡力U產(chǎn)生對穩(wěn)定狀態(tài)U=0； 對于失穩(wěn)狀態(tài)，阻滑力發(fā)揮到最大值，使得不平衡力最小化。由此我們可以看出穩(wěn)定和失穩(wěn)狀態(tài)的共性之處，就是阻滑力在其許可范圍內(nèi)R≤fN+cA，其取值必使不平衡力最小化。加固的本質(zhì)是在加固力Q的幫助下，滑塊恢復(fù)平衡狀態(tài)，T=R+Q。所以不平衡力和加固力大小相等，方向相反。

U=T-R

(9)

如果將滑動力和阻滑力分別視為作用力和抗力，則前述彈塑性結(jié)構(gòu)分析一般性提法顯然也適用于穩(wěn)定或失穩(wěn)狀態(tài)下的單滑塊模型：結(jié)構(gòu)作用力與抗力的差值必然趨向最小值，該差值就是殘余不平衡力。

4.2 材料的不平衡力

由式(6)可知，對于一給定的作用應(yīng)力和屈服條件，抗力應(yīng)力必使式(10)所示的余能范數(shù)最小，即minE(σ)：

(10)

式(10)是一個條件極值問題，這個問題的解對應(yīng)著彈塑性增量型的正交流動法則和一致性條件。為解出minE(σ)，構(gòu)造Lagrange函數(shù)：

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

(11)

則應(yīng)有：

(12)

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

(13)

式(13)就是彈塑性增量型的正交流動法則和一致性條件。由此可知作用力與抗力差值趨向最小值是彈塑性本構(gòu)關(guān)系的要求。

4.3 結(jié)構(gòu)的不平衡力

受式(10)的啟發(fā)，可以推測結(jié)構(gòu)作用應(yīng)力場與抗力應(yīng)力場，在總塑性余能意義下差距最小，即最小塑性余能原理：

(14)

式(14)是在整個結(jié)構(gòu)的積分。塑性余能是不平衡力的范數(shù)，所以最小塑性余能原理意味著結(jié)構(gòu)趨于不平衡力最小的狀態(tài)。

事實上，對于式(14)的推測可以在理論上給予嚴(yán)格的證明。Yang et al.(2012)從彈-黏塑性結(jié)構(gòu)出發(fā)，采用Duvaut-Lions模型證明了三維結(jié)構(gòu)非平衡狀態(tài)最終演化趨勢滿足最小塑性余能原理。因此，不平衡力無論在材料層次還是結(jié)構(gòu)層次都有嚴(yán)格的理論基礎(chǔ)。

圖 13 外力場和自承力場動態(tài)調(diào)整圖Fig. 13 Dynamic adjustment diagram of external and self-supporting force fields

有了最小塑性余能原理這一指導(dǎo)性原則，不平衡力的求解就成了一個雙場優(yōu)化的迭代過程：通過變形調(diào)整達到σ1和σ最接近的結(jié)構(gòu)整體塑性余能最小狀態(tài)(圖 13)。

潘家錚(1980)最大最小原理可由最小塑性余能予以說明：(1)最小值原理：滑坡如能沿許多滑面滑動，則失穩(wěn)時，它將沿抵抗力最小的一個滑面破壞——加固力最小。(2)最大值原理：滑坡體的滑面肯定時，則滑面上的反力(以及滑坡體內(nèi)的內(nèi)力)能自行調(diào)整，以發(fā)揮最大的抗滑能力——自承力最大。最小塑性余能原理反映了多自由度超靜定體系的非線性內(nèi)力分配原則； 而潘家錚最大最小原理用于確定在極限狀態(tài)下多滑塊體系的內(nèi)力分配。

從計算策略上來講，非線性有限元本質(zhì)是時空離散化?？臻g離散化是指有限單元在空間上近似無限自由度變形體，等效基礎(chǔ)是最小勢能原理； 而時間離散化是指無限小微分加載過程的離散化或增量化，等效基礎(chǔ)是最小塑性余能原理。最小塑性余能原理突破了經(jīng)典彈塑性理論，允許應(yīng)力超出屈服面。

5 結(jié) 論

基于理想彈塑性有限元分析，巖體結(jié)構(gòu)無法消除的殘余不平衡力就是其變形破壞的內(nèi)在驅(qū)動力。研究案例表明不平衡力分布與結(jié)構(gòu)開裂破壞關(guān)系密切。

本文對彈塑性結(jié)構(gòu)分析給出更一般性的提法：在保持結(jié)構(gòu)連續(xù)性及材料無損的前提下，結(jié)構(gòu)作用力與抗力的差值必然趨向最小值(體現(xiàn)為最小塑性余能原理)，該差值就是殘余不平衡力。這種提法包容了出現(xiàn)不平衡力情況，而常規(guī)結(jié)構(gòu)彈塑性分析的解就是該差值為0的一個特例。最小塑性余能原理可視為彈塑性理論的增量化的要求。施加反向殘余不平衡力作為加固力，理論上可完全避免結(jié)構(gòu)開裂損傷，保持結(jié)構(gòu)連續(xù)性。

需要說明的是，和單滑塊模型不同，變形體結(jié)構(gòu)出現(xiàn)不平衡力并不意味著結(jié)構(gòu)要整體失穩(wěn)，更多是出現(xiàn)結(jié)構(gòu)局部損傷破壞。