脆塑性迭代逼近算法的改進(jìn)

金俊超，景來紅，楊風(fēng)威，宋志宇，尚朋陽(yáng)

(1.黃河勘測(cè)規(guī)劃設(shè)計(jì)研究院有限公司，河南 鄭州 450003；2.水利部黃河流域水治理與水安全重點(diǎn)實(shí)驗(yàn)室（籌），河南 鄭州 450003；3.中水北方勘測(cè)設(shè)計(jì)研究有限責(zé)任公司，天津 300222)

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

當(dāng)巖石強(qiáng)度弱化過程分解為一系列的脆性與塑性交互發(fā)生的過程時(shí)，軟化求解問題可以轉(zhuǎn)化為一系列脆塑性求解問題，規(guī)避經(jīng)典塑性力學(xué)求解中對(duì)于軟化速率的限制[18].Wang等[18]假設(shè)應(yīng)力跌落過程中最小主應(yīng)力不變，應(yīng)力從峰值強(qiáng)度跌落至殘余強(qiáng)度時(shí)Lode參數(shù)也不變，提出脆塑性迭代逼近算法及相應(yīng)的有限元求解過程.雖然脆塑性迭代逼近算法概念明確、方法簡(jiǎn)單，具有很強(qiáng)的工程應(yīng)用價(jià)值，但是最小主應(yīng)力不變跌落方法能否正確描述巖石脆塑性變形破壞過程中的應(yīng)力變化，引入其他應(yīng)力跌落方法是否必要，都有待進(jìn)一步分析論證.沈新普等[19]假定應(yīng)力跌落過程中各應(yīng)力偏量分量的原有比例保持不變，即應(yīng)力跌落是從初始屈服面沿著徑向向后繼屈服面跌落的，提出偏應(yīng)力等比例跌落方法，并推導(dǎo)彈脆塑性本構(gòu)積分的數(shù)值格式.Zheng等[7]認(rèn)為跌落時(shí)塑性應(yīng)變?cè)隽康姆较蚍纤苄晕粍?shì)理論，提出塑性位勢(shì)跌落方法，推導(dǎo)給出巖塊D-P(Drucker-Prager)準(zhǔn)則和結(jié)構(gòu)面M-C準(zhǔn)則的塑性跌落因子的計(jì)算方法.舒芹等[20]基于帶拉伸截?cái)嗟腗-C準(zhǔn)則，假設(shè)應(yīng)力跌落過程中應(yīng)力球量不變，推導(dǎo)破壞粒子應(yīng)力的計(jì)算方法.但是上述最小主應(yīng)力不變跌落方法、偏應(yīng)力等比例跌落方法、塑性位勢(shì)跌落方法等能否正確描述巖石脆塑性應(yīng)力跌落過程中的應(yīng)力變化，不同方法各存在什么問題，仍缺少系統(tǒng)性的論證分析.

在脆塑性迭代逼近算法[18]中，采用正確的應(yīng)力跌落方法是應(yīng)變軟化計(jì)算模擬的關(guān)鍵.金俊超等[21]初步分析了不同應(yīng)力跌落方法的合理性，但理論推導(dǎo)過程中主應(yīng)力和應(yīng)力張量概念不清，未分析應(yīng)力球量不變跌落方法合理性.本研究1）對(duì)已有不同方法存在的問題進(jìn)行系統(tǒng)嚴(yán)謹(jǐn)?shù)恼撟C分析；2）考慮脆塑性變形破壞過程中的泊松效應(yīng)，改進(jìn)塑性位勢(shì)跌落方法，克服已有脆塑性迭代逼近算法中缺陷，在此基礎(chǔ)上實(shí)現(xiàn)彈塑性變形破壞全過程的模擬；3）通過數(shù)值算例驗(yàn)證研究結(jié)果的正確性.

1 不同應(yīng)力跌落方法的問題分析

1.1 已有方法的理論假定與表達(dá)式

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

式中：sp為峰值強(qiáng)度面偏應(yīng)力張量.跌落后的殘余強(qiáng)度面應(yīng)力張量

最小主應(yīng)力不變跌落方法[18]針對(duì)以壓應(yīng)力為正的情況，假設(shè)應(yīng)力跌落過程中最小主應(yīng)力不變，應(yīng)力從峰值強(qiáng)度跌落至殘余強(qiáng)度時(shí)Lode參數(shù)不變.對(duì)于以拉應(yīng)力為正的情況，有

式中：下標(biāo)1、2、3分別為最大主應(yīng)力方向、中間應(yīng)力方向和最小主應(yīng)力方向；上標(biāo)p、r分別對(duì)應(yīng)峰值強(qiáng)度面和殘余強(qiáng)度面.

Zheng等[7]通過理論推導(dǎo)，證實(shí)在脆塑性變形破壞過程中應(yīng)力的功不負(fù)，脆塑性材料仍滿足Il’yushin公設(shè).在此基礎(chǔ)上，假定跌落過程產(chǎn)生的塑性應(yīng)變?cè)隽喀う舙的方向滿足塑性位勢(shì)理論：

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

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

式中：λ為L(zhǎng)amé參數(shù)，G為切變模量.

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

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

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

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

表3 塑性位勢(shì)跌落方法存在的問題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 塑性位勢(shì)跌落方法的改進(jìn)及驗(yàn)證

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

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

原有塑性位勢(shì)跌落方法[7]假設(shè)脆塑性變形破壞過程中總應(yīng)變保持不變，這與如圖2所示的試驗(yàn)結(jié)果一致.Zheng等[7]未考慮泊松效應(yīng)，簡(jiǎn)單認(rèn)為彈性應(yīng)變的增加等于負(fù)的塑性應(yīng)變，由此得到式(6).事實(shí)上，巖石作為彈塑性材料，存在泊松效應(yīng).以主應(yīng)力、主應(yīng)變型式進(jìn)行說明，在脆塑性變形破壞過程中，某個(gè)主方向發(fā)生彈性變形（等于負(fù)的塑性應(yīng)變）由于泊松效應(yīng)，必然在其他2個(gè)正交主方向產(chǎn)生橫向彈性應(yīng)變，對(duì)應(yīng)的回彈應(yīng)變?cè)隽?/p>

圖2 紅砂巖單軸拉伸試驗(yàn)結(jié)果[22]Fig.2 Uniaxial tension test of red sandstone[22]

因此，脆塑性變形破壞過程中的彈性應(yīng)變?cè)隽喀う舉應(yīng)為該方向上負(fù)的塑性變形Δεp和泊松效應(yīng)引起的回彈應(yīng)變?cè)隽喀う興e的疊加：

式(6)與式(9)存在顯著差異，證明考慮脆塑性變形破壞過程中的泊松效應(yīng)后，彈性應(yīng)變?cè)隽坑?jì)算公式發(fā)生變化，應(yīng)力增量計(jì)算公式也相應(yīng)改變：

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

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

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

2.2 塑性位勢(shì)跌落方法的改進(jìn)

在三維應(yīng)力空間考慮脆塑性變形破壞過程中的泊松效應(yīng)，對(duì)原有塑性位勢(shì)跌落方法進(jìn)行改進(jìn).塑性應(yīng)變?cè)隽喀う舙、回彈應(yīng)變?cè)隽喀う興e及彈性應(yīng)變?cè)隽喀う舉的計(jì)算公式分別為

跌落過程應(yīng)力增量

式中：D為彈性剛度矩陣，由彈性模量E和泊松比v構(gòu)成.殘余應(yīng)力

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

2.3 室內(nèi)試驗(yàn)?zāi)M驗(yàn)證

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

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

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

表6 2種應(yīng)力跌落計(jì)算方法的殘余應(yīng)力模擬結(jié)果對(duì)比Tab.6 Simulated residual stress of two stress-drop calculation method MPa

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

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

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

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

3 脆塑性迭代逼近法的改進(jìn)及驗(yàn)證

3.1 脆塑性迭代逼近法的改進(jìn)

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

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

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

3.2 算例驗(yàn)證

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

1.3.3 潛在生態(tài)危害指數(shù)法。瑞典科學(xué)家Hakanson[10]提出的生態(tài)危害指數(shù)法是目前最為流行的一種對(duì)土壤或沉積物中土壤重金屬污染進(jìn)行評(píng)價(jià)的方法。該法將重金屬的生態(tài)效應(yīng)、環(huán)境效應(yīng)與毒理學(xué)聯(lián)系在一起，不僅反映了某一特定環(huán)境中各種污染物對(duì)環(huán)境的影響，及多種污染物的綜合效應(yīng)，而且用定量的方法劃分出了潛在生態(tài)風(fēng)險(xiǎn)的程度。其計(jì)算公式為：

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

圖7 Tennessee大理巖及三峽花崗巖的三軸壓縮試驗(yàn)?zāi)MFig.7 Experimental simulation of triaxial compression for Tennessee marble and Sanxia granite

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

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

圖9 圓隧力學(xué)響應(yīng)有限元模擬結(jié)果與解析解對(duì)比Fig.9 Comparisons of numerical and theoretical results of circular tunnel excavation

4 工程算例

4.1 Mine-by試驗(yàn)洞破壞區(qū)模擬

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

表7 Mine-by隧洞圍巖計(jì)算參數(shù)Tab.7 Model parameters of surrounding rock of Mine-by tunnel

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

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

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

圖12 監(jiān)測(cè)的隧洞塑性區(qū)和破壞區(qū)[29]Fig.12 Monitored plastic and failure zones[29]

4.2 某水電站引水隧洞輔助洞松動(dòng)圈模擬

為了了解圍巖在開挖擾動(dòng)下的損傷區(qū)范圍，在某水電站引水隧洞輔助洞A和輔助洞B的多個(gè)斷面進(jìn)行聲波測(cè)試.以輔助洞A中的AK10+900斷面為研究對(duì)象，根據(jù)聲波測(cè)試結(jié)果確定的開挖松動(dòng)圈[31]如圖13所示，通過本研究所提應(yīng)變軟化有限元計(jì)算程序模擬開挖塑性區(qū)，檢驗(yàn)本研究所提方法的正確性.輔助洞的有限元模型如圖14所示，水平方向和豎直方向各取70 m，隧洞軸向取3 m；共劃分單元55 170個(gè)，開挖面附近單元寬度為0.2 m；模型底部及四周表面采用法向約束.初始地應(yīng)力場(chǎng)為σ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].參考已有研究及室內(nèi)試驗(yàn)結(jié)果[31-32]，綜合確定圍巖計(jì)算參數(shù)如表8所示.

表8 輔助洞圍巖計(jì)算參數(shù)Tab.8 Model parameters of surrounding rock of auxiliary tunnel

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

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

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

圖15 開挖模擬結(jié)果Fig.15 Simulated results by proposed method

5 結(jié) 論

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

（2）考慮脆塑性變形破壞過程中的泊松效應(yīng)，提出改進(jìn)塑性位勢(shì)跌落方法，推導(dǎo)具體應(yīng)力更新過程，并嵌入軟件Abaqus.

（3）巖石室內(nèi)多種類型破壞試驗(yàn)?zāi)M結(jié)果證實(shí)：改進(jìn)塑性位勢(shì)跌落方法能夠正確模擬多種情況下的巖石脆塑性變形破壞過程.改進(jìn)塑性位勢(shì)跌落方法克服了原有塑性位勢(shì)跌落方法的缺陷，為脆塑性迭代逼近算法的改進(jìn)奠定了基礎(chǔ).

（4）采用所提改進(jìn)塑性位勢(shì)跌落方法替換原有脆塑性迭代逼近算法中的應(yīng)力跌落計(jì)算方法，能夠?qū)崿F(xiàn)彈塑性應(yīng)變軟化過程的數(shù)值模擬.引入塑性強(qiáng)化算法，實(shí)現(xiàn)了彈塑性變形破壞全過程的正確計(jì)算.

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

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