脆塑性迭代逼近算法的改進

金俊超，景來紅，楊風威，宋志宇，尚朋陽

(1.黃河勘測規劃設計研究院有限公司，河南 鄭州 450003；2.水利部黃河流域水治理與水安全重點實驗室（籌），河南 鄭州 450003；3.中水北方勘測設計研究有限責任公司，天津 300222)

數值算法是有限元計算的關鍵，對計算結果的準確性和穩定性有直接影響.研究者已在巖石彈塑性及黏彈塑性算法上取得很多成果，如Clausen等[1-4]針對理想彈塑性積分算法中M-C(Mohr-Coulomb)準則和H-B(Hoek-Brown)準則屈服面拐點問題開展研究，Kindrachuk等[5-6]研究了黏彈性和黏塑性積分算法.由于應力-應變關系曲線表現為負斜率，導致剛度矩陣非正定，使得應變軟化行為的有限元數值模擬較為困難[7].S?rensen等[8]將應力狀態轉換為主應力空間，提出基于H-B準則的彈塑性應變軟化隱式積分算法，但文獻[8]并未給出有限元求解過程的詳細說明，考慮到經典塑性力學對軟化速率有限制[7]，其提出的隱式積分算法可能無法求解峰后軟化速率較大的情況.Karaveli?等[9]以線性強度準則為例，提出彈塑性應變強化-軟化隱式積分算法，但該算法也可能無法求解峰后軟化速率較大的情況[7].雖然Jia等[10-17]提出各自的彈塑性損傷模型，并編寫UMAT子程序將模型嵌入軟件Abaqus，但是當損傷變量為0時，模型退化為彈塑性應變軟化模型，所提積分算法也將遇到可能無法求解峰后軟化速率較大情況的問題[7].

當巖石強度弱化過程分解為一系列的脆性與塑性交互發生的過程時，軟化求解問題可以轉化為一系列脆塑性求解問題，規避經典塑性力學求解中對于軟化速率的限制[18].Wang等[18]假設應力跌落過程中最小主應力不變，應力從峰值強度跌落至殘余強度時Lode參數也不變，提出脆塑性迭代逼近算法及相應的有限元求解過程.雖然脆塑性迭代逼近算法概念明確、方法簡單，具有很強的工程應用價值，但是最小主應力不變跌落方法能否正確描述巖石脆塑性變形破壞過程中的應力變化，引入其他應力跌落方法是否必要，都有待進一步分析論證.沈新普等[19]假定應力跌落過程中各應力偏量分量的原有比例保持不變，即應力跌落是從初始屈服面沿著徑向向后繼屈服面跌落的，提出偏應力等比例跌落方法，并推導彈脆塑性本構積分的數值格式.Zheng等[7]認為跌落時塑性應變增量的方向符合塑性位勢理論，提出塑性位勢跌落方法，推導給出巖塊D-P(Drucker-Prager)準則和結構面M-C準則的塑性跌落因子的計算方法.舒芹等[20]基于帶拉伸截斷的M-C準則，假設應力跌落過程中應力球量不變，推導破壞粒子應力的計算方法.但是上述最小主應力不變跌落方法、偏應力等比例跌落方法、塑性位勢跌落方法等能否正確描述巖石脆塑性應力跌落過程中的應力變化，不同方法各存在什么問題，仍缺少系統性的論證分析.

在脆塑性迭代逼近算法[18]中，采用正確的應力跌落方法是應變軟化計算模擬的關鍵.金俊超等[21]初步分析了不同應力跌落方法的合理性，但理論推導過程中主應力和應力張量概念不清，未分析應力球量不變跌落方法合理性.本研究1）對已有不同方法存在的問題進行系統嚴謹的論證分析；2）考慮脆塑性變形破壞過程中的泊松效應，改進塑性位勢跌落方法，克服已有脆塑性迭代逼近算法中缺陷，在此基礎上實現彈塑性變形破壞全過程的模擬；3）通過數值算例驗證研究結果的正確性.

1 不同應力跌落方法的問題分析

1.1 已有方法的理論假定與表達式

偏應力等比例跌落方法[19]假設應力由峰值屈服面徑向跌落至殘余屈服面上，有且僅有各向應力偏量按同一比例β衰減.設峰值強度面應力張量為σp，殘余強度系數為β，跌落過程的應力偏量變化量

式中：sp為峰值強度面偏應力張量.跌落后的殘余強度面應力張量

最小主應力不變跌落方法[18]針對以壓應力為正的情況，假設應力跌落過程中最小主應力不變，應力從峰值強度跌落至殘余強度時Lode參數不變.對于以拉應力為正的情況，有

式中：下標1、2、3分別為最大主應力方向、中間應力方向和最小主應力方向；上標p、r分別對應峰值強度面和殘余強度面.

Zheng等[7]通過理論推導，證實在脆塑性變形破壞過程中應力的功不負，脆塑性材料仍滿足Il’yushin公設.在此基礎上，假定跌落過程產生的塑性應變增量Δεp的方向滿足塑性位勢理論：

式中：Δa為塑性流動因子；考慮到通常采用剪脹角來衡量巖體的剪脹變形，采用M-C準則型式塑性勢函數.假定脆塑性應力跌落過程中總應變保持不變，則彈性應變的增加等于負的塑性應變，以主應變形式表示為

根據式(5)計算跌落過程中的應力增量Δσ，以主應力、主應變形式表示為

式中：λ為Lamé參數，G為切變模量.

應力球量不變跌落方法[20]假設應力由峰值屈服面徑向跌落至殘余屈服面上，應力跌落前后的應力圓半徑發生變化，圓心保持不變，且應力跌落過程中應力主軸不旋轉，對于以拉應力為正的情況，有

1.2 不同方法存在的問題分析

如圖1所示，巖石的變形破壞分為單軸拉伸破壞、拉剪破壞、純剪破壞及壓剪破壞.由于直接從三維應力空間對已有應力跌落方法存在的問題進行理論推導分析較為困難，本研究以經典的線性M-C準則為例，從單軸拉伸破壞、單軸壓縮破壞及二向純剪破壞特征點，分析現有方法的模擬正確性.如果單軸壓縮破壞、單軸拉伸破壞及二向純剪破壞不能正確計算，那么拉剪、壓剪等復合破壞類型也無法正確計算.從主應力空間推導不同方法的脆塑性應力跌落計算過程，理論公式如表1～4所示.可以發現，已有脆塑性應力跌落方法，雖然能描述某種脆塑性變形破壞過程，但無法正確模擬全部破壞方式.在實際工程中，圍巖變形破壞往往多種破壞方式并存，已有方法無法適用于不同破壞方式的計算，阻礙了計算模擬的正確性.

表2 最小主應力不變跌落方法存在的問題Tab.2 Defects of existing calculation method based on constant minor principal stress in brittle-plastic process

表3 塑性位勢跌落方法存在的問題Tab.3 Defect of existing calculation method based on plastic potential theory

表4 球量不變跌落方法存在的問題Tab.4 Defect of existing calculation method based on invariant spherical stress

圖1 巖石破壞分類Fig.1 Schematic of rock failure types

2 塑性位勢跌落方法的改進及驗證

偏應力等比例跌落方法[19]、最小主應力不變跌落方法[18]及應力球量不變跌落方法[20]均采用人為假定應力跌落路徑.塑性位勢跌落方法[7]滿足彈塑性理論中的Il’yushin公設，無需人為假定應力跌落路徑，理論更為嚴謹，本研究將針對塑性位勢跌落方法存在的不足，考慮泊松效應進行方法改進.

2.1 脆塑性變形破壞過程中的泊松效應分析說明

原有塑性位勢跌落方法[7]假設脆塑性變形破壞過程中總應變保持不變，這與如圖2所示的試驗結果一致.Zheng等[7]未考慮泊松效應，簡單認為彈性應變的增加等于負的塑性應變，由此得到式(6).事實上，巖石作為彈塑性材料，存在泊松效應.以主應力、主應變型式進行說明，在脆塑性變形破壞過程中，某個主方向發生彈性變形（等于負的塑性應變）由于泊松效應，必然在其他2個正交主方向產生橫向彈性應變，對應的回彈應變增量

圖2 紅砂巖單軸拉伸試驗結果[22]Fig.2 Uniaxial tension test of red sandstone[22]

因此，脆塑性變形破壞過程中的彈性應變增量Δεe應為該方向上負的塑性變形Δεp和泊松效應引起的回彈應變增量Δεde的疊加：

式(6)與式(9)存在顯著差異，證明考慮脆塑性變形破壞過程中的泊松效應后，彈性應變增量計算公式發生變化，應力增量計算公式也相應改變：

在巖石彈塑性變形破壞過程中泊松比存在一定變化，為了簡化，本研究以泊松比v為定值進行理論推導，對于考慮泊松比變化的情況，理論推導過程一樣，只要將v考慮為非定值.

如表5所示，考慮脆塑性變形破壞過程中的泊松效應，分別針對單軸拉伸破壞、單軸壓縮破壞和二向純剪破壞等破壞類型進行計算理論推導.可以發現，此時塑性位勢跌落方法可以正確計算多種情況下的巖石脆塑性變形破壞過程.

表5 基于泊松效應的改進塑性位勢跌落方法的合理性Tab.5 Rationality of improved calculation method based on plastic potential theory considering Poisson’s effect

2.2 塑性位勢跌落方法的改進

在三維應力空間考慮脆塑性變形破壞過程中的泊松效應，對原有塑性位勢跌落方法進行改進.塑性應變增量Δεp、回彈應變增量Δεde及彈性應變增量Δεe的計算公式分別為

跌落過程應力增量

式中：D為彈性剛度矩陣，由彈性模量E和泊松比v構成.殘余應力

式中：n+1、n表示第n+1和第n次迭代；p1、p3為以主應力表示的迭代過程應力增量.設Δa0=0，當不超過預設精度時，跳出迭代；將Δan+1代入式(11)，計算得到脆塑性應力跌落過程的塑性應變增量Δεp；根據式(12)，計算得到回彈應變增量Δεde；根據式(13)，計算得到彈性應變增量Δεe；根據式(14)、(15)，計算得到跌落過程應力增量Δσ及最終殘余應力σr.根據上述計算步驟編制UMAT子程序，將改進塑性位勢跌落方法嵌入軟件Abaqus，實現脆塑性變形破壞問題的有限元求解.

2.3 室內試驗模擬驗證

2.3.1 紅砂巖單軸拉伸試驗模擬驗證 根據紅砂巖試驗資料[22]，構建標準圓柱體有限元模型，試件直徑為50 mm，高為100 mm.在模型底部設置法向鏈桿約束，頂部施加位移荷載，進行單軸拉伸試驗模擬.如圖3所示為本研究所提方法模擬結果與試驗結果的對比.可以看到，雖然由于將峰前簡化為線彈性，導致模擬結果與試驗結果存在一定差異，但整體上，本研究所提方法可以正確模擬紅砂巖單軸拉伸脆塑性變形破壞過程，模擬的殘余階段軸向及橫向應力均為0（點B、B'），與試驗結果一致.將原有塑性位勢跌落方法模擬的跌落后應力點繪于圖中（點C、C' ）.對于脆塑性變形破壞過程，在峰值及殘余強度參數給定的情況下，不同脆塑性應力跌落方法的模擬結果主要的差異是跌落后的應力模擬結果，因此本研究將分析重點放在跌落后的應力模擬結果正確性上.原方法計算的跌落后軸向拉應力為-0.03 MPa，橫向應力變為-0.20 MPa，與試驗結果不符，證明原有塑性位勢跌落方法不能正確模擬巖石單軸拉伸脆塑性變形破壞過程.

圖3 紅砂巖單軸拉伸應力-應變模擬曲線對比Fig.3 Comparison of simulated stress-strain curves of red sandstone under uniaxial tension

2.3.2 花崗巖三軸壓縮試驗模擬驗證 根據花崗巖試驗資料[23]，構建標準圓柱體有限元模型，在模型底部設置法向鏈桿約束.在頂部及環向施加靜水壓力至額定圍壓，再在頂部施加位移荷載，進行三軸壓縮試驗模擬.如圖4所示為本研究所提方法模擬結果與試驗結果的對比.由于殘余強度參數擬合值與試驗結果的誤差，導致本研究模擬的殘余強度點與試驗結果存在差異（點B1-B3和B′1-B′3）.如表6所示為提取殘余強度模擬結果與試驗數據擬合值.由表可知，模擬的殘余階段軸向應力與擬合結果一致，模擬的橫向應力始終與圍壓保持一致，證明本研究所提方法能夠正確模擬巖石三軸壓縮脆塑性變形破壞過程.將原有塑性位勢跌落方法計算的跌落后應力繪于圖中（點C1-C3、C1′-C3′），相應數據見表6，原方法在不同圍壓的跌落后應力模擬結果，均與試驗數據擬合值存在很大差異，證明原有塑性位勢跌落方法不能正確模擬巖石三軸壓縮脆塑性變形破壞過程.

表6 2種應力跌落計算方法的殘余應力模擬結果對比Tab.6 Simulated residual stress of two stress-drop calculation method MPa

圖4 花崗巖三軸壓縮應力-應變模擬曲線對比Fig.4 Comparison of simulated stress-strain curves of granite under triaxial compression

2.3.3 砂巖壓剪試驗模擬驗證 根據砂巖壓剪試驗資料[24]，構建40 mm×40 mm×40 mm的立方體試件.試件的底面采用法向鏈桿約束，頂面施加壓應力1.5 MPa；再在試件上部施加法向鏈桿約束，試件下部施加水平位移荷載，進行壓剪試驗模擬.如圖5所示為改進方法與原有塑性位勢跌落方法的應力-應變模擬曲線對比，點A、A′為峰值和殘余剪應力試驗點，點B為本研究所提方法模擬的殘余應力點，點C為原方法模擬的殘余應力點.可以看到，本研究所提方法模擬的殘余階段剪應力與試驗結果一致，證明本研究所提方法可以正確模擬砂巖壓剪脆塑性變形破壞過程；原有方法模擬的殘余剪應力與試驗值相差33.5%，驗證了原方法不能正確模擬砂巖壓剪脆塑性變形破壞過程.其中σn為法向壓應力.

圖5 砂巖壓剪應力-應變模擬曲線對比Fig.5 Comparison of simulated stress-strain curves of sandstone under compression and shear

單軸拉伸試驗、三軸壓縮試驗及壓剪試驗模擬涵蓋主要的變形破壞類型，充分證實了改進塑性位勢跌落方法能夠正確模擬多種情況下的巖石脆塑性變形破壞過程；原方法存在缺陷，無法正確模擬多種情況下的巖石脆塑性變形破壞過程.

3 脆塑性迭代逼近法的改進及驗證

3.1 脆塑性迭代逼近法的改進

脆塑性迭代逼近法的核心是脆塑性應力跌落過程的模擬[18].塑性強化及理想塑性積分算法研究較為成熟，Clausen等[1-4,21]通過建立強度參數隨塑性應變的演化方程，在向后歐拉隱式積分算法的基礎上引入塑性硬化參與迭代計算，使得屈服面在迭代過程中與應力一同更新(屈服面的大小會隨著塑性損傷的更新而變化)直至應力回到屈服面上.本研究將對如何在所提改進塑性位勢跌落方法基礎上實現應變強化-軟化數值模擬展開說明.

如圖6所示，將改進塑性位勢跌落方法與理想塑性流動算法結合，通過多級脆塑性跌落-塑性流動，實現彈塑性應變軟化的數值模擬，克服已有脆塑性迭代逼近算法[18]中應力跌落方法不正確的缺陷.將改進的彈塑性應變軟化算法與塑性強化算法結合，實現巖石塑性強化-峰后軟化-殘余強度的變形破壞全過程模擬.通過編譯UMAT子程序，將計算過程在軟件Abaqus中實現，其中塑性強化和理想塑性流動過程采用向后歐拉隱式積分算法，具體參見文獻[21].

圖6 巖石彈塑性變形破壞全過程計算流程Fig.6 Numerical procedures of full elasto-plastic deformation and failure process

3.2 算例驗證

已有脆塑性應力跌落方法均存在不足，以此為基礎建立的脆塑性迭代逼近算法也存在缺陷；所提改進塑性位勢跌落方法具備正確性，以此為基礎建立的脆塑性迭代逼近算法具備正確性.本研究將通過室內三軸壓縮試驗和應變軟化圓隧圍巖力學響應規律模擬，檢驗所提方法的正確性，不再與原方法的模擬結果進行對比.

1.3.3 潛在生態危害指數法。瑞典科學家Hakanson[10]提出的生態危害指數法是目前最為流行的一種對土壤或沉積物中土壤重金屬污染進行評價的方法。該法將重金屬的生態效應、環境效應與毒理學聯系在一起，不僅反映了某一特定環境中各種污染物對環境的影響，及多種污染物的綜合效應，而且用定量的方法劃分出了潛在生態風險的程度。其計算公式為：

3.2.1 室內三軸壓縮試驗模擬驗證 韓建新等[25]基于M-C準則，假設黏聚力和內摩擦角為最大主應變的分段線性函數，提出巖石應變軟化模型，并模擬Tennessee大理巖三軸壓縮試驗；沈華章等[26]分別假定峰前和峰后階段強度參數隨軟化參數按照指數函數演化，塑性變形遵守非關聯流動法則，對三峽花崗巖的三軸壓縮試驗進行數值模擬.如圖7所示為Tennessee大理巖及三峽花崗巖三軸壓縮試驗結果和模擬結果的對比.可以看到，模擬結果與試驗數據具有良好的可比性，實現了Tennessee大理巖峰后塑性軟化-殘余強度變形破壞全過程的正確計算，也實現了三峽花崗巖塑性強化-峰后脆性跌落-塑性軟化-殘余強度變形破壞全過程的正確計算.對比結果驗證了所提彈塑性變形破壞全過程數值算法的正確性.

圖7 Tennessee大理巖及三峽花崗巖的三軸壓縮試驗模擬Fig.7 Experimental simulation of triaxial compression for Tennessee marble and Sanxia granite

3.2.2 應變軟化圓隧圍巖力學響應規律模擬 如圖8所示，假定強度參數和剪脹角隨塑性剪應變分段線性演化規律，Lee等[27]將潛在塑性區圍巖按等圍壓釋放劃分為若干同心圓，采用塑性區按相等的應力增量劃分、圍壓逐漸遞減的差分方法，給出了應變軟化圓隧圍巖力學響應規律解析解；Park[28]用一階常微分方程代替應力平衡、本構關系和一致性條件的偏微分方程，采用Runge-Kutta求解常微分方程，給出圓隧圍巖力學響應規律解析解.如圖9所示為圓隧收斂位移及應力分布解析解和本研究有限元模擬結果對比.可以看到，本研究所提方法計算的圍巖特征曲線與Lee等[27]和Park[28]的解析解十分接近，所提方法計算的支護應力pi=0時的圍巖徑向及切向應力分布與Park[28]的解析解也具有較好的可比性，證明所述的數值實現工作在整體模型層次正確有效.

圖8 圓隧圍巖挖掘及有限元模型Fig.8 Example of circular tunnel excavation and FEM

圖9 圓隧力學響應有限元模擬結果與解析解對比Fig.9 Comparisons of numerical and theoretical results of circular tunnel excavation

4 工程算例

4.1 Mine-by試驗洞破壞區模擬

Mine-by試驗圓洞長46 m，直徑3.5 m，埋深420 m，圍巖為Lac du Bonnet花崗巖.試驗洞采用非爆破的機械開挖方法，開挖過程中圍巖不斷發生脆性剝落破壞，形成典型的V形脆性破壞區[29].由于監測資料完整，已被包括Hajiabdolmajid等[29-30]在內的多位學者用以驗證所建力學模型的合理性.本研究將利用所提有限元計算程序引用已有本構模型和參數，對試驗洞進行開挖模擬，并與監測結果對比，驗證所提方法的正確性.Mine-by隧洞的有限元模型如圖10所示，水平方向和豎直方向各取17.5 m，沿隧洞軸線方向取1 m；共劃分單元18 732個，開挖面附近單元寬度為0.1 m，厚度方向劃分3層網格；在模型底部及四周表面采用法向約束；圍巖計算參數如表7所示.

表7 Mine-by隧洞圍巖計算參數Tab.7 Model parameters of surrounding rock of Mine-by tunnel

圖10 Mine-by隧洞有限元模型Fig.10 Numerical model of Mine-by tunnel

如圖11所示，圖中εep為等效塑性應變，c為黏聚力.開挖導致高地應力集中在巷道頂部和底板，這些部位的強度參數變化最為明顯.模擬的頂拱塑性區最大深度為2.56 m，如圖12（a）所示，微震數據確定的最大深度為2.57 m，數值與分布區域基本一致.模擬的頂拱破壞區最大深度為2.27 m，如圖12（b）所示，接近實測的破壞區最大深度2.3 m，形狀與實測情況基本吻合.

圖11 模擬的隧洞塑性區和破壞區Fig.11 Simulated plastic and failure zones by proposed method

圖12 監測的隧洞塑性區和破壞區[29]Fig.12 Monitored plastic and failure zones[29]

4.2 某水電站引水隧洞輔助洞松動圈模擬

為了了解圍巖在開挖擾動下的損傷區范圍，在某水電站引水隧洞輔助洞A和輔助洞B的多個斷面進行聲波測試.以輔助洞A中的AK10+900斷面為研究對象，根據聲波測試結果確定的開挖松動圈[31]如圖13所示，通過本研究所提應變軟化有限元計算程序模擬開挖塑性區，檢驗本研究所提方法的正確性.輔助洞的有限元模型如圖14所示，水平方向和豎直方向各取70 m，隧洞軸向取3 m；共劃分單元55 170個，開挖面附近單元寬度為0.2 m；模型底部及四周表面采用法向約束.初始地應力場為σx=-48.98 MPa，σy=-55.67 MPa，σz=-66.16 MPa，τxy=-2.52 MPa，τyz=-0.30 MPa，τxz=7.17 MPa[31].參考已有研究及室內試驗結果[31-32]，綜合確定圍巖計算參數如表8所示.

表8 輔助洞圍巖計算參數Tab.8 Model parameters of surrounding rock of auxiliary tunnel

圖13 AK10+900斷面開挖松動圈[31]Fig.13 Plastic zone at section AK10+900[31]

圖14 輔助洞有限元模型Fig.14 Numerical model of auxiliary tunnel

如圖15所示，開挖導致高地應力集中在兩側邊墻及拱角位置，這些部位的強度參數變化最為明顯.本研究所提方法模擬的松弛深度平均值為1.35 m，現場試驗監測平均值為1.66 m，相差0.31 m.與已有模型相比，文獻[31]模擬的測線塑性區深度與監測結果的總方差為8.30，文獻[32]模擬的總方差為10.96，本研究所提方法模擬的為7.11，較已有研究分別減小了16.59%和32.13%.對比結果表明，本研究所提方法能夠較為合理地描述輔助洞圍巖的彈塑性變形破壞現象.

圖15 開挖模擬結果Fig.15 Simulated results by proposed method

5 結 論

（1）在主應力空間，結合特征點變形破壞特征，論證偏應力等比例跌落方法、最小主應力不變跌落方法及塑性位勢跌落方法均存在不足：不能正確模擬多種情況下巖石脆塑性變形破壞過程.

（2）考慮脆塑性變形破壞過程中的泊松效應，提出改進塑性位勢跌落方法，推導具體應力更新過程，并嵌入軟件Abaqus.

（3）巖石室內多種類型破壞試驗模擬結果證實：改進塑性位勢跌落方法能夠正確模擬多種情況下的巖石脆塑性變形破壞過程.改進塑性位勢跌落方法克服了原有塑性位勢跌落方法的缺陷，為脆塑性迭代逼近算法的改進奠定了基礎.

（4）采用所提改進塑性位勢跌落方法替換原有脆塑性迭代逼近算法中的應力跌落計算方法，能夠實現彈塑性應變軟化過程的數值模擬.引入塑性強化算法，實現了彈塑性變形破壞全過程的正確計算.

（5）不同地質條件的Mine-by試驗洞和某水電站引水隧洞輔助洞開挖模擬結果與現場監測結果均吻合良好，證明本研究所提方法能夠合理模擬工程中圍巖彈塑性變形破壞現象.

（6）本研究基于線性M-C準則，提出了改進塑性位勢迭代逼近算法，實現了彈塑性應變軟化過程的數值求解.非線性Hoek-Brown準則、巖石彈塑性損傷變形破壞過程，有待在本研究基礎上深入拓展.