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

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

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

0 引言

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

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

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

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

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

1 角點型非共軸本構(gòu)模型

(1)

(2)

(3)

或

(4)

(5)

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

(6)

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

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

(7)

(8)

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

(9)

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

體積模量K為

(10)

式中ν為泊松比。

2 本構(gòu)模型的數(shù)值實現(xiàn)

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

2.1 應(yīng)力更新算法

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

(11)

(12)

(13)

值得指出的是，方程組(13)包含了11個式子，式中與對應(yīng)的未知量構(gòu)成方程組解向量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)對應(yīng)的殘差向量為R(X)，在采用Newton-Raphson迭代法求解該非線性方程組時，該殘差向量的Jacobi矩陣可表示為?R(X)/?X，這是一個11×11的方陣，其形式為

(15)

2.2 一致性剛度矩陣

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

(16)

結(jié)合式(13)可得：

(17)

(18)

(19)

(20)

2.3 數(shù)值算法驗證

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

表1 模型參數(shù)Table 1 Model parameters

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

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

以應(yīng)力和應(yīng)變張量為代表的張量Tkl主軸的旋轉(zhuǎn)過程可以用圖1(b)中坐標(biāo)系的旋轉(zhuǎn)過程表示。求解張量主軸方向的過程在數(shù)學(xué)上就是對與張量對應(yīng)的實對稱矩陣進(jìn)行對角化處理的過程，其形式為

(21)

ABAQUS內(nèi)置的命令能夠?qū)崿F(xiàn)該對角化過程。這里規(guī)定主軸方向需要逆時針旋轉(zhuǎn)坐標(biāo)系的情況為正，順時針為負(fù)，且旋轉(zhuǎn)角的取值范圍為[-180°，180°]。圖2和圖3展示了共軸和不同非共軸參數(shù)下模擬單剪試驗的應(yīng)力應(yīng)變關(guān)系和非共軸角變化關(guān)系。這些結(jié)果比較好地驗證了半隱式算法在程序?qū)嵤r的有效性，并且保證了計算結(jié)果能夠符合非共軸模型反映的力學(xué)特性。

圖2 應(yīng)力應(yīng)變關(guān)系Fig. 2 Relationship between shear stress and shear strain

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

3 本構(gòu)模型的數(shù)值實現(xiàn)

3.1 有限元分析模型

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

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

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

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

3.2 考慮非共軸性的主應(yīng)力軸旋轉(zhuǎn)對隧道開挖的影響

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

3.2.1 塑性區(qū)發(fā)展

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

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

圖5 塑性區(qū)分布情況Fig. 5 Distribution of plastic parameter

圖6 監(jiān)測點的應(yīng)力路徑Fig. 6 Stress paths at monitoring points

3.2.2 應(yīng)力主軸的旋轉(zhuǎn)

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

(a) 共軸情況 (b) 非共軸hn=2G (c) 非共軸hn=1G (d) 非共軸hn=0.5G (e) 非共軸hn=0.2G圖7 應(yīng)力旋轉(zhuǎn)角分布情況(單位： 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 剪切應(yīng)力分布情況(單位： kPa)Fig. 8 Distribution of shear (unit: kPa)

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

圖9 開挖過程中監(jiān)測點處的應(yīng)力旋轉(zhuǎn)變化Fig. 9 Stress rotation at monitoring points during excavation

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

3.2.3.1 地表沉降

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

(a) 應(yīng)力釋放階段

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

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

(a) 應(yīng)力釋放階段

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

i0=S0H。

(22)

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

3.2.3.2 隧道收斂變形

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

(a) 應(yīng)力釋放階段

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

4 結(jié)論與建議

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

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

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

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

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

總體來講，與主應(yīng)力旋轉(zhuǎn)現(xiàn)象密切相關(guān)的非共軸變形特性會顯著影響隧道變形分析的結(jié)果，是變形預(yù)測研究中不可忽視的重要因素，值得開展更深入的研究。