基于彈塑性損傷-滲流耦合模型的隧道洞口段圍巖穩定性研究

萬友生

(南昌軌道交通集團有限公司，江西 南昌 330013)

0 引言

隧道洞口段一般具有埋深較淺、圍巖質量差、偏壓作用明顯等特點，因此在開挖過程中容易失穩，發生工程事故。對此，相關學者對隧道洞口段的開挖穩定性問題展開研究。張惠民[1]利用FLAC軟件分析不同埋深和支護方式對隧道洞口段圍巖和支護結構的穩定性；宋洋等[2]通過數值模擬和分析現場監測數據的方法指出施工過程中的重點鞏固部位；皇民等[3]利用數值模擬和模型實驗研究地震對洞口段工程穩定性的影響規律；楊佳奇等[4]總結寒區隧道洞口段工程穩定性的影響因素，并引入模糊評判的方法，對其進行風險評估；鄭明新等[5]利用強度折減法分析鋼花管注漿在洞口段隧道施工中的加固效果；郭小紅等[6]通過歸納工程資料，總結出洞口段隧道洞內塌方和邊坡失穩的10大致災因素；任洋等[7]通過離心臺實驗對強震作用下的隧道洞口段動力響應特征和規律展開研究。

巖土工程中，隧道工程開挖產生的應力重分布作用會使周圍巖體中的裂隙、孔洞等微缺陷發育貫通，形成損傷[8-9]。損傷的存在會使巖石宏觀力學性質如彈性模量、單軸抗壓強度等發生弱化，威脅施工安全。同時，力學場的變化使得巖石骨架結構發生變化，裂隙等的貫通進一步加劇地下水的滲流作用，改變孔隙水壓分布規律，通過有效應力原理反作用于巖石力學場。因此，在施工設計時還應充分考慮應力-滲流的耦合作用。

現有研究中，針對孔隙水滲流作用下的邊坡或隧道單一結構類型的穩定性研究較多[10-12]，但針對隧道洞口段這一邊坡-隧道復合結構的相關研究報道較為少見,且研究過程中均未考慮到巖石的荷載損傷問題。對此，本文建立基于Morh-Coulomb(M-C)準則的巖石彈塑性損傷-滲流耦合模型，基于完全隱式的向后歐拉積分算法編寫應力求解程序，并通過ABAQUS的子程序接口實現應力-損傷-滲流的完全耦合計算；以福建某在建隧道為工程依托，建立三維有限元計算模型，分析降雨滲流作用對洞口段邊坡和隧道明洞的影響規律，旨在為類似工程的安全設計，提供一定的理論依據。

1 基于M-C準則的巖石彈塑性損傷-滲流耦合模型

1.1 基于M-C準則的巖石彈塑性損傷模型

考慮損傷的M-C準則如式(1)所示：

(1)

式中：f為屈服函數；I1為應力張量第一不變量；J2為應力偏量第二不變量；θ為羅德角；φ為內摩擦角，°；c為黏聚力，Pa；ω為損傷變量。

根據相關研究給出損傷演化方程，如式(2)所示[13]：

(2)

(3)

式中：εp1，εp2，εp3分別為3個主應力。

損傷變量演化規律如圖1所示。

圖1 不同κ值下損傷變量演化規律

此外，由Lemaitre有效應力和應變等效假設理論可知損傷材料的應力-應變關系如式(4)所示：

σ=(1-ω)Dε

(4)

式中：D為材料的彈性剛度矩陣；σ為應力矩陣；ε為應變矩陣。

損傷作用下巖體的應力應變曲線如圖2所示。

圖2 損傷作用下巖石應力-應變曲線

1.2 巖石應力-滲流耦合方程

由有效應力原理可知，滲流場產生的孔隙水壓值對巖石應力場造成影響，如式(5)所示：

(5)

根據Kozeny-Carman公式，非破裂區巖石滲透系數可以表示為體積應變的函數，如式(6)所示：

(6)

式中：K為巖石的滲透系數，m/s。K0為初始滲透系數，m/s；n0為初始孔隙度；εv為體積應變。

對于破裂損傷區的巖體，其滲透系數如式(7)所示[14]：

K=K0·10ξ(A′e-ω/α+B′)

(7)

式中：ξ為跳突系數；A′,B′，α分別為經驗參數。

2 耦合模型的數值求解算法

2.1 M-C準則的數值積分算法

由于M-C準則在屈服面上存在棱線、尖點等數值不連續特征，造成應力積分困難。為保證計算精度，摒棄角點光滑法等近似解法，從主應力空間出發，建立M-C準則的完全隱式向后歐拉算法，避免數值“奇異點”的問題。M-C準則在π平面上的圖形如圖3所示。

圖3 M-C準則在π平面上的圖形

為減少應力維數，簡化計算過程，在主應力空間中對問題進行討論。由于空間應力的對稱性，只在σ1>σ2>σ3區域分析問題即可。

主應力空間中M-C準則的表達式如式(8)所示：

(8)

式中：f1～f6分別為不同應力區域內屈服函數的表達式；σ1為最大主應力；σ2為中間主應力；σ3為最小主應力。

用膨脹角ψ代替式(8)中的內摩擦角φ可得到與屈服函數具有相同形式的塑性勢函數gi。

整個應力計算過程分為彈性預測、塑性修正和損傷修正3個部分：

1)彈性預測

彈性預測公式如式(9)所示：

(9)

(10)

式中：σn+1為tn+1時刻的計算應力。

2)塑性修正

映射應力可能存在的區域為應力平面f，棱線l1或l2及尖點P處。為判斷應力所在區域，利用邊界面法對應力區域進行劃分，邊界面方程的定義如式(11)所示[15]：

(11)

(12)

(13)

(14)

(15)

(16)

式中：p為向量長度，取1；a2為f2對應力的導數。

σa的表達式如式(17)所示：

(17)

式中：k為與摩擦角有關的常數。

同理可以建立邊界面方程如式(18)所示：

(18)

建立邊界面后，將試算應力帶入，當pⅠ-Ⅱ≥0且pⅠ-Ⅲ≤0時，應力返回至屈服面；當pⅠ-Ⅱ<0且pⅠ-Ⅲ<0，應力返回至棱線l1;當pⅠ-Ⅱ>0且pⅠ-Ⅲ>0時，應力返回值棱線l2；否則應力返回至尖點。

判定返回區域后，根據塑性增量理論，分別建立N-R方程對更新應力進行求解，更新應力表達式如式(19)所示：

(19)

式中：Δσp為塑性應力增量；Δλ為塑性因子。

3)損傷修正

塑性應變求解方程如式(20)所示：

(20)

將求得的塑性應變代入式(2)～(3)中計算損傷變量。在新求得的損傷變量基礎上，對應力進行再次修正，修正方程如式(21)所示：

(21)

2.2 流固耦合計算流程

基于ABAQUS軟件中的UMAT和USDFLD子程序接口實現流固耦合的計算過程，計算流程如圖4所示。

納入標準：①患者被本院醫師診斷為輸卵管妊娠；②患者的隨訪依從性較高；③患者自愿參與本次試驗，且簽署知情同意書。排除標準：①患者合并其他較為嚴重的疾?。虎诨颊叻亲栽竻⑴c本次試驗，或者未簽署知情同意書。兩組患者一般資料比較,差異無統計學意義(P>0.05)，具有可比性。本研究獲得醫院倫理委員會批準。

圖4 流固耦合計算流程

3 洞口段計算模型

3.1 工程概況

蘇橋隧道位于福建省蘇橋村，為雙線小凈距隧道，長360 m左右，最大埋深60 m。根據現場調查和鉆探揭露，所選典型斷面主要地質組成由上而下主要為強風化石英砂巖和中風化石英砂礫巖。

3.2 有限元計算模型

根據蘇橋隧道施工設計圖紙建立有限元計算模型，如圖5所示。建立的邊坡數值模型中總節點數為38 622個，單元數為29 199個。設置靜水位線，水位線以下孔隙水壓成梯形分布。模型底部為不透水邊界。坡面設置為降雨邊界條件。從最不利角度出發，巖石計算參數參照Ⅴ級圍巖進行選取，損傷參數κ=20，α=0.3。隧道主要支護形式為襯砌。計算參數如表1所示。

表1 洞口段模型計算參數

圖5 有限元計算模型

依據相關文獻及當地水文資料，擬定3種降雨條件分別進行計算，具體如表2所示。

表2 3種降雨類型參數

4 計算結果分析

4.1 降雨強度對復合結構安全系數的影響

利用重力加載法對不同降雨強度下的邊坡安全系數進行計算，結果如圖6所示。

由圖6可知，在同一降雨強度下，安全系數隨著降雨時間的增加而不斷下降，降雨初期的安全系數變化較降雨后期劇烈，在第28 h左右安全系數逐漸趨于穩定。對于3種不同類型降雨來說，孔隙水壓力變化幅度越大，邊坡安全系數下降的幅度越大，大于降雨條件下的安全系數變化幅度最大為6.46%。在第48 h，3種類型降雨的安全系數從小到大排序為：大雨類型<中雨類型<小雨類型，最小安全系數為1.07。

圖6 不同降雨類型下洞口段安全系數變化規律

4.2 降雨強度對隧道圍巖變形規律的影響

設置降雨量分別為0.002，0.004，0.006，0.008，0.01 mm/h，降雨時長為72 h，對降雨過程中的隧道拱頂沉降、拱底隆起和拱腰收斂變化規律進行計算，結果如圖7所示。

圖7 不同降雨強度下隧洞變形曲線

由圖7可知，降雨初期隧道內各個監測點的位移均發生明顯變化。隨著時間發展，土體逐漸飽和，位移變化速率逐漸平滑，最后趨近于零。降雨強度越大，最終位移值越大。由于洞口段的偏壓作用，隧道左洞和右洞的位移值變化規律不同。對拱頂監測點而言，在相同降雨強度下，右洞的沉降值均要大于左洞沉降值。而左洞拱底隆起值大于右洞拱底隆起值。2洞的拱肩向不同方向收斂。

圖7(a)中，對某時間段強降雨過程中拱頂位移變化監測值與計算值進行對比。由對比結果可以看出，數值計算結果在數值上略有偏差，變化規律基本一致，證明本文模型的準確性與工程實用性。

對降雨過程中(降雨量0.01 mm/h)的損傷區進行計算，結果如圖8所示。

由圖8可知,當降雨持續60 h后，巖體內部開始產生損傷區，損傷區范圍與損傷值大小隨著降雨進程的逐漸發展。塑性區發生在左洞與坡面之間，即孔隙水壓最大值所在處，計算結果再次證明降雨對洞口段邊坡造成的安全隱患。

5 結論

1)降雨降低洞口段整體結構的安全系數，降雨強度越大，安全系數降低的越快。當降雨達到一定時長后，由于土體發生飽和，安全系數逐漸趨于穩定。

2)降雨會造成隧洞內部圍巖發生變形，拱頂沉降、拱底隆起及拱肩收斂值均與降雨強度有關；淺埋與深埋段隧道的拱肩收斂規律相反。

3)降雨會使洞口段局部發生破壞，尤其是淺埋隧道處，施工時應密切注意，以防發生工程事故。