考慮非共軸性的隧道開挖引起的地表沉降數值分析

陳洲泉, 陳湘生, 龐小朝, 蘇 棟, 林星濤

(1. 深圳大學土木與交通學院， 廣東 深圳 518060； 2. 深圳大學 濱海城市韌性基礎設施教育部重點實驗室，廣東 深圳 518060； 3. 深圳大學未來地下城市研究院， 廣東 深圳 518060； 4. 深圳市地鐵地下車站綠色高效智能建造重點實驗室， 廣東 深圳 518060； 5. 鐵科院(深圳)研究設計院有限公司， 廣東 深圳 518060)

0 引言

隨著我國城市化進程的迅猛發展，許多大城市人口激增，地上空間建(構)筑物愈發密集，地面交通擁堵問題越來越突出，已經嚴重影響市民生活質量以及社會經濟活動的正常運行。地下鐵路建設以其運量大、效率高、空間集約等特點正成為解決城市交通擁堵問題的重要手段[1]。然而，在地鐵開挖過程中勢必會對周圍土體造成擾動，超限的土體變形必然會影響到周圍建筑物及地下管線安全。因此，合理評估隧道開挖引起的地面沉降一直是工程安全風險控制的重點和難點[2]，尤其是在人口密集的城市商業圈和經濟圈等地區。

值得指出的是，在地下隧道工程建設中，普遍存在著隧道開挖誘發開挖面周圍應力場改變的情況。這種變化既存在應力大小的改變，也存在主應力方向的改變[3-5]。已有的模擬主應力軸旋轉作用的土單元試驗(如單剪[6-7]和空心圓柱環剪試驗[8-10])表明，土體的力學性質會受主應力軸旋轉的影響，顯現出變形非共軸特性，即材料當前的主應力方向與塑性主應變率方向不一致。在通用的工程數值分析中采用的材料本構模型一般基于傳統的共軸塑性理論建立，即假定塑性應變率的方向與當前應力狀態方向一致，不考慮非共軸性。在一定程度上，這種對非共軸塑性的忽略，會給基于數值模擬方法獲得的隧道工程設計和施工方案安全性評估結果帶來不確定性，這對于變形控制要求高的隧道工程而言是十分不利的。

為了在工程計算中考慮主應力軸旋轉引起的非共軸變形，首先，需要構建非共軸本構模型。為此，學者們已經開展了相當豐富的理論工作，當前非共軸塑性理論主要有以Rudnicki等[11]提出的所謂角點理論，Yu等[12]延伸的double shearing理論以及Lashkari等[13]在Mohr圓應力空間構建的圓切向流動理論。角點理論由于其表述形式的一般性，在后續研究中被進一步發展，形成了比較完備的體系。錢建固等[14]在考慮了第三應力不變量的影響后將該方法擴展到了更一般的三維應力狀態。Li等[15]、Tsutsumi等[16]以及陳洲泉等[17-19]將這類非共軸理論與各向異性本構模型相結合，很好地模擬主應力旋轉試驗中土體的非共軸力學行為。

其次，在采用非共軸模型對具體工程問題開展數值模擬方面，與本構理論研究相比，相關的數值研究在廣度和深度上顯得明顯不足。Yang等[20-22]較早地將角點型非共軸模型應用于具體工程問題的數值分析，比較了共軸與非共軸模型在分析淺基礎地基承載力和變形問題上的差異，指出不考慮非共軸塑性的計算結果偏于不安全；另外，還分析了非共軸塑性對倉筒問題中主應力軸旋轉和土壓力分布所產生的影響。Yuan等[23]則在double shearing理論框架下發展了二維應力狀態下各向異性非共軸模型，分析了地基承載問題。然而，將非共軸模型應用于隧道工程分析的案例屈指可數。只有袁冉等[24]將文獻[23]發展的各向異性非共軸模型應用于隧道開挖誘發地表沉降問題的分析。

考慮到這方面工作的缺乏，筆者針對隧道開挖問題進一步開展了探究。與袁冉等[24]的工作相比，本文的創新之處在于采用了截然不同的非共軸本構模型框架和模型數值積分算法。本文采用了基于一般應力狀態構建的角點型非共軸本構模型，既適用于平面應變情況，也適用于一般三維情況；其次，區別于文獻[24]中采用的基于誤差控制的自動劃分子增量步的顯式方法，本文采用了半隱式數值積分算法。通過將模型數值實現，針對隧道開挖問題，從主應力軸旋轉、塑性區發展以及地表沉降等方面，比較共軸與非共軸模型計算結果的差異，進一步提高對此類問題的認識。

1 角點型非共軸本構模型

(1)

(2)

(3)

或

(4)

(5)

式(3)—(5)中：Cijkl為非共軸張量；sij=σij-δijσkk/3為應力偏張量；hn為非共軸模量，表征應力增量引發非共軸塑性變形的難易程度。hn可以在共軸模型參數確定后，通過擬合單剪試驗應力應變曲線以及旋轉角變化曲線來獲得。

(6)

式中：rij為塑性流動方向；Q為塑性勢函數； dλ為塑性標量因子。在這里采用考慮雙曲硬化法則的Drucker-Prager模型來定義共軸塑性行為，該模型的屈服函數f和塑性勢函數形式Q為

f=q-η(p+ccotφc)；
Q=q-ptanψ。

(7)

(8)

參考HS模型[25]中彈性模量的定義方法，本文中剪切模量G同樣采用應力相關的形式

(9)

式中：Gref為參考剪切模量，即p=pref時的剪切模量;pref=100 kPa為參考壓力。

體積模量K為

(10)

式中ν為泊松比。

2 本構模型的數值實現

本文采用Moran等[26]提出的半隱式積分算法。這種算法在應力更新過程中，塑性標量因子為隱式表達，塑性流動方向則為顯式表達。這意味著在1個增量步中，塑性標量因子要在增量步中求解，而塑性流動方向在增量步初始即已給定。

2.1 應力更新算法

在ABAQUS/Standard計算平臺下，1個計算增量步的應力更新過程為根據增量步給定的應變增量dεij來計算應力增量dσij。根據式(1)和(2)，對于任意增量步時刻tn到tn+1的應力更新格式為

(11)

(12)

(13)

值得指出的是，方程組(13)包含了11個式子，式中與對應的未知量構成方程組解向量X，表示為

X=(n+1σ11,n+1σ22,n+1σ33,n+1σ12,n+1σ21,n+1σ13,

n+1σ31,n+1σ23,n+1σ32,n+1η,dλ)。

(14)

指定方程組(13)對應的殘差向量為R(X)，在采用Newton-Raphson迭代法求解該非線性方程組時，該殘差向量的Jacobi矩陣可表示為?R(X)/?X，這是一個11×11的方陣，其形式為

(15)

2.2 一致性剛度矩陣

由于ABAQUS/Standard分析模塊采用隱式平衡迭代方法求解平衡方程，在構造平衡迭代的剛度矩陣時需要提供與應力更新過程相適應的剛度矩陣，即一致性剛度矩陣。該剛度矩陣的定義為

(16)

結合式(13)可得：

(17)

(18)

(19)

(20)

2.3 數值算法驗證

在ABAQUS/Standard分析模塊中，基于用戶材料子程序開發接口UMAT，編寫文中介紹的非共軸模型的應力積分算法，并通過模擬單剪試驗的過程來驗證算法的有效性以及模型描述非共軸變形的能力。單剪有限元分析模型采用單個減縮積分平面應變單元CPE4R。本構模型參數的選取參考HS模型中黏土的取值，如表1所示。

表1 模型參數Table 1 Model parameters

試驗加載方式如圖1(a)所示，豎向壓力保持不變，σ22=500 kPa，初始靜止側壓力系數K0=0.5，在單元上邊緣施加水平位移實現剪切過程。

圖1 單剪試驗和張量主軸旋轉Fig. 1 Simple shear test and rotation of principal axis of tensor

以應力和應變張量為代表的張量Tkl主軸的旋轉過程可以用圖1(b)中坐標系的旋轉過程表示。求解張量主軸方向的過程在數學上就是對與張量對應的實對稱矩陣進行對角化處理的過程，其形式為

(21)

ABAQUS內置的命令能夠實現該對角化過程。這里規定主軸方向需要逆時針旋轉坐標系的情況為正，順時針為負，且旋轉角的取值范圍為[-180°，180°]。圖2和圖3展示了共軸和不同非共軸參數下模擬單剪試驗的應力應變關系和非共軸角變化關系。這些結果比較好地驗證了半隱式算法在程序實施時的有效性，并且保證了計算結果能夠符合非共軸模型反映的力學特性。

圖2 應力應變關系Fig. 2 Relationship between shear stress and shear strain

圖3 非共軸角的變化Fig. 3 Variation of non-coaxial angle

3 本構模型的數值實現

3.1 有限元分析模型

本文建立了平面應變條件下模擬半無限空間黏土地層中隧道開挖過程的有限元模型，如圖4所示。

圖4 隧道模型尺寸及網格劃分(單位： m)Fig. 4 Geometry and finite mesh of tunnel model (unit: m)

圖4展示了計算區域土體的有限元網格，考慮研究問題的對稱性，取模型的一半，即尺寸為60 m×60 m。土體的有限元網格劃分選用4節點平面應變完全積分單元(CPE4)。圓形隧道中心位于土體表面以下20 m位置，半徑為4 m；襯砌支護厚度為300 mm，采用平面應變不可壓縮完全積分單元(CPE4I)，襯砌與隧道開挖面用tie連接。襯砌采用彈性模型，彈性模量E=19 GPa，泊松比ν=0.2，重度γ=25 kN/m3。在隧道周圍選取6個土體單元作為監測點(A，B，C，D，E，F)。

模型分析過程分為3步： 1)地應力平衡。首先，約束模型左右兩側豎向邊界的水平位移和底部位移。整個土體施加重力，重度γ=20 kN/m3。然后，模型內部定義沿豎向線性分布的地應力為σv=γH(H為隧道中心距土體表面的豎向距離)，水平地應力為σh=K0σv(K0為靜止土壓力系數，取為0.5)。在這一步殺死襯砌單元，隧道范圍內的土體單元不殺死。2)應力釋放。這一步利用初始平衡態獲取開挖區域土體對于開挖邊界的支反力代替開挖土體的支護作用，通過釋放部分支反力的方法來模擬土體穩定后的地層損失。因此，這一步將開挖區域的土體單元殺死，并在隧道邊界節點處施加節點荷載，該節點荷載衰減為平衡態節點力的50%。 3)施加襯砌。激活襯砌單元，并將節點荷載衰減為0，開挖完成。

3.2 考慮非共軸性的主應力軸旋轉對隧道開挖的影響

隧道開挖必然會對周圍土體產生擾動。這種擾動包括土體塑性區發展的范圍和程度、應力應變場大小和方向的重分布。土體擾動引起的地表沉降和隧道收斂變形是實際工程中能夠直接而且準確監測到的響應，能在一定程度上反映開挖影響。下文將針對非共軸變形特性這一土體力學性質的差異對于隧道開挖引起土體擾動的影響進行討論。

3.2.1 塑性區發展

圖5展示了共軸和非共軸模型條件下隧道開挖后周圍土體塑性區的分布情況。圖中采用應力比表征土體的塑性發展程度。圖5彩圖展示了土體塑性發展程度的分布情況，應力比的取值為0.75～1。圖中紅灰2色圖展示了應力比大于0.75(紅色區域)的區域分布，即發生塑性屈服的區域。由圖可以看出： 開挖會引起隧道腰部土體塑性的顯著發展，沿斜上和斜下2個方向發展出2條塑性帶，形成交叉式分布；離隧道開挖區域越近，塑性發展的程度越深。土體考慮非共軸性后，隨著非共軸塑性模量的減小，2條塑性帶沿著各自的方向向更遠處延伸；同時，在2條塑性帶所夾的區域中也發展出更多的塑性區。

圖6展示了隧道周圍6個監測點在隧道開挖過程中的應力路徑變化。由于在不同非共軸參數條件下的計算得到的應力路徑變化規律大致相同，所以選取了具有代表性的塑性模量hn=0.5G的情況與共軸情況進行比較。圖中標注字母的位置為開挖的起始點，每條曲線上標注的點號代表開挖過程中應力釋放階段與施加襯砌階段的交接處。圖中斜直線的斜率為0.75，直線以下的區域為彈性區。由圖6可以看到所有應力路徑的起始點都位于直線附近，表明計算的初始狀態正確。整個開挖過程中，隧道頂部(A)和底部(F)都處于卸載狀態，即始終處于彈性區，這與圖5中塑性區分布規律是相符的。塑性區土體(B,C,D,E)的應力路徑主要在開挖的應力釋放階段發生塑性屈服，進入施加襯砌階段，B，E的應力路徑沒有顯著變化，但C，D的應力路徑使得土體的應力比減小，即發生彈性卸載。值得注意的是，發生塑性屈服的土體的應力路徑在共軸和非共軸條件下的差別較大，處于彈性狀態的則差別較小。

圖5 塑性區分布情況Fig. 5 Distribution of plastic parameter

圖6 監測點的應力路徑Fig. 6 Stress paths at monitoring points

3.2.2 應力主軸的旋轉

圖7和圖8展示了共軸和非共軸模型不同參數條件下隧道開挖后土體應力旋轉角和剪切應力的分布情況，圖中旋轉角的取值范圍為[-5°, 90°]。由圖可見： 應力旋轉發生在隧道頂部和底部，并且在隧道頂部的地表以下以及如圖5所示的2條塑性帶與彈性區交界的區域也都有顯著的應力旋轉發生；而隧道腰部土體的應力旋轉現象不顯著。同時，隨著非共軸塑性模量的降低，模型計算得到的主應力旋轉區域越廣，旋轉程度也越大。與應力旋轉分布相對應的是，在隧道周圍發生顯著應力旋轉的區域同時伴隨著顯著的剪切應力集中發展現象，并且隨著非共軸塑性模量的降低，剪切應力集中發生區域向外延伸，與應力旋轉區域大致重合。這與剪切應力引起主應力旋轉這一論述相符。

(a) 共軸情況 (b) 非共軸hn=2G (c) 非共軸hn=1G (d) 非共軸hn=0.5G (e) 非共軸hn=0.2G圖7 應力旋轉角分布情況(單位： kPa)Fig. 7 Distribution of rotation angle of stress (unit: kPa)

(a) 共軸情況 (b) 非共軸hn=2G (c) 非共軸hn=1G (d) 非共軸hn=0.5G (e) 非共軸hn=0.2G圖8 剪切應力分布情況(單位： kPa)Fig. 8 Distribution of shear (unit: kPa)

圖9展示了隧道開挖過程中監測點處應力旋轉在計算增量步中的變化情況。從圖中可以看到，在2個分析步結束位置，隧道周圍區域共軸模型與非共軸計算得到的應力旋轉情況相差不大。值得注意的是，土體的應力旋轉主要發生在開挖的應力釋放階段；而在襯砌施作階段，監測點A,B，C，D，E的應力旋轉角基本保持不變，點F則發生較大程度的應力旋轉。根據各監測點土體應力旋轉的方向和大小，各處的大主壓應力基本同隧道輪廓線平行。比如： 隧道正側面C，D點的應力在整個開挖階段都沒有發生顯著的應力旋轉，基本保持豎直向上的方向。這與圖7和圖8展示的分布規律相符。

圖9 開挖過程中監測點處的應力旋轉變化Fig. 9 Stress rotation at monitoring points during excavation

3.2.3 地表沉降和隧道收斂變形

3.2.3.1 地表沉降

圖10展示了在共軸和非共軸情況下隧道開挖的應力釋放階段和完成階段的地表沉降曲線。由圖10可以發現: 在模型條件相同的條件下，開挖完成的沉降曲線(圖10(b))比應力釋放階段的曲線(圖10(a))低；隨著非共軸模量的減小，地表沉降曲線的沉降槽更深，且各曲線都是在離對稱軸線足夠遠處才相交重合。

(a) 應力釋放階段

(b) 開挖完成階段圖10 開挖完成時的地表沉降曲線Fig. 10 Ground settlement curves at end of excavation

圖11進一步總結了圖10中各曲線斜率沿水平距離的變化規律。特別值得指出的是,各曲線斜率的最大值(曲線反彎點)在距離對稱軸7.7～10 m的區間取得，即沉降曲線的沉降槽寬度也處于該區間。表明非共軸塑性模量對地表沉降槽寬度的計算結果影響較小。根據O′Reilly等[27]歸納的沉降槽經驗公式

(a) 應力釋放階段

(b) 開挖完成階段圖11 沉降曲線斜率沿水平距離的變化Fig. 11 Slope of ground settlement curves along horizontal distance

i0=S0H。

(22)

式中：i0為沉降槽寬度；S0為經驗系數(黏土取0.4～0.7，砂土取0.2～0.3)；H為隧道中心到地表的距離。本文沉降槽寬度的經驗取值為8～14 m，與計算結果相當接近，表明了數值模擬結果的合理性。

3.2.3.2 隧道收斂變形

圖12比較了開挖2個階段中不同模型參數計算得到的隧道收斂變形差異。圖中隧道收斂變形是將隧道邊界點位移放大40倍后得到的結果。從圖12可以看出： 在開挖的2個階段中，模型非共軸參數主要影響隧道頂部和側面向內的收斂變形，不影響隧道底部的向上位移；同時，當模型非共軸模量小到一定程度(如0.2G)，土體對隧道側向的作用要大于對隧道頂部的作用，表現為隧道出現側向壓扁的傾向。對比圖12(a)與(b)可以看出，襯砌施作后，由于隧道管片的自重較大，隧道表現為豎向進一步受壓的傾向。因此，隧道豎向收斂變形的程度相對于應力釋放階段進一步加劇；而隧道側向被擠出，橫向收斂程度小于應力釋放階段。

(a) 應力釋放階段

(b) 開挖完成階段圖12 隧道收斂變形Fig. 12 Convergence of tunnel

4 結論與建議

采用半隱式應力更新算法，編寫了考慮雙曲剪切硬化的DP模型的非共軸模型用戶材料子程序UMAT。在ABAQUS/Standard的分析模塊中模擬半無限空間地層內襯砌隧道開挖過程，并針對非共軸性對隧道開挖的影響進行參數分析，得到以下結論。

1)隧道開挖過程會引起隧道腰部周圍土體沿斜上和斜下2個方向發展出2條交叉的塑性帶。考慮非共軸變形后，模型的非共軸塑性模量越小，這2條塑性帶延伸的越遠，同時2條塑性帶所夾區域發展出更多塑性區,而對于隧道頂部和底部的土體則一直處于彈性階段。

2)開挖過程中土體的應力主軸方向的改變是一個漸變的過程，開挖面周圍土體的大主壓應力方向由豎直方向逐漸向隧洞輪廓線形狀靠攏。同時在2條塑性帶上下外緣與彈性區也存在顯著的主應力旋轉，模型的非共軸變形越強，主應力旋轉的區域和程度也越大，這表明隧洞頂部和底部彈性區的土體沿塑性帶向隧道移動更多。

3)本構模型的非共軸效應越強，隧道開挖引起的最大地表沉降值越大，地表沉降曲線的沉降槽越深。但是，非共軸性并不顯著影響沉降槽的寬度。同時，更大的地表沉降也與土體內斜向上塑性帶向地表更廣地延伸，以及更大范圍的隧道上部應力旋轉區。

4)隧道收斂變形中，頂部和側面向隧道內收斂變形的程度會隨模型非共軸塑性模量的減小而增大，而非共軸性對于底部位移的影響較小。

總體來講，與主應力旋轉現象密切相關的非共軸變形特性會顯著影響隧道變形分析的結果，是變形預測研究中不可忽視的重要因素，值得開展更深入的研究。