非均質黏土隧道開挖環向穩定性的剛體平動有限元上限分析

劉國良， 黃飛衡， 羅靜靜， 楊 峰

(1. 中鐵隧道勘察設計研究院有限公司， 廣東 廣州 511458； 2. 中南大學土木工程學院， 湖南 長沙 410075)

0 引言

黏土地層普遍存在非均質、各向異性等特點[1-3]。其中，非均質性體現在土體抗剪強度沿深度方向線性增加，這將影響隧道開挖過程中的地層穩定性。通常環向角度下的隧道開挖面穩定性，適用于分析盾構法隧道開挖工作段較長以及礦山法隧道剛性支護未施作的情況。在非均質黏土地層條件下，關于隧道環向開挖面穩定性的研究已有不少文獻報道。黃茂松等[4]應用多塊體極限分析上限法揭示了非均質黏土地層圓形隧道開挖面穩定性最優上限解及其規律。Wilson等[5]采用剛性塊體上限法和有限元極限分析法，揭示了抗剪強度隨深度線性增加時黏土地層方形隧道開挖面穩定性。Ukritchon等[6]以三維有限元極限分析方法研究了抗剪強度沿深度線性增加條件下的隧道穩定性，總結了埋深比、土體非均質強度、土體重度等因素的影響規律。Chen等[7]利用離散元技術改進的三維破壞模式，在極限分析理論框架下結合土拱效應評估了方形隧道開挖面穩定性。Zhang等[8]運用極限分析運動學方法，構建出馬蹄形隧道地層連續速度場，探討了不排水條件下黏土地層隧道開挖面失穩破壞機制。

對于抗剪強度沿深度線性增加的非均質純黏土地層隧道縱向開挖面穩定性問題，文獻[9]應用剛體平動運動單元上限有限元進行了多因素綜合影響分析，揭示了滑移線網破壞模式的形態特征。對于非均質黏土地層隧道而言，沿著橫斷面的環向開挖面穩定性同樣值得關注，特別是不同橫斷面形狀的影響規律。此外，地層潛在破壞特征和擴展范圍的定量分析，對于穩定性評價及工程措施的制定具有重要的理論意義，現階段這方面的探討仍然少見。剛體平動運動單元上限有限元采用非線性數學優化方法，具有搜索獲取精細滑移線網破壞模式的特色功能，可作為破壞區域定量分析的有力手段[9-10]。

本文沿用剛體平動運動單元上限有限元方法[9-10]，對不排水抗剪強度沿深度線性增加的非均質黏土地層隧道環向開挖面極限狀態下的穩定性進行研究，重點探討圓形、馬蹄形及橢圓形等4種代表性隧道斷面形狀下，不同的隧道埋深比、土體非均質強度和土體重度等參數與臨界荷載比上限解及失穩破壞模式的關聯規律，揭示滑移線網破壞模式的特征規律，給出破壞范圍的定量幾何尺度，以期為地層加固等工程措施提供參考與借鑒。

1 問題描述與計算模型

1.1 非均質黏土地層隧道環向開挖面穩定性模型

將非均質黏土地層隧道環向開挖面穩定性問題簡化為沿橫斷面剖開的二維平面應變模型，即沿著隧道縱向力學狀態一致，如圖1和圖2所示。如此簡化與三維開挖面空間問題有所差異，不過用來分析沿橫斷面的環向開挖面潛在破壞模式的形態和擴展范圍是合適的，從穩定性評價結論方面也是保守的。

(a) 圓形隧道

(b) 橢圓形隧道 (c) 橢圓形隧道 (d) 馬蹄形隧道

圖2 非均質黏土地層圓形隧道環向開挖面穩定性分析模型

工程現場常見圓形和馬蹄形隧道斷面形式，而矢跨比小的瘦高型隧道(如單線鐵路隧道、連拱隧道中導洞等)以及矢跨比大的大跨型隧道(如3車道、4車道公路隧道等)，可近似為不同跨高比的橢圓斷面，以此簡化研究不同斷面形式隧道環向開挖面穩定性和破壞形態的演化規律，可為隧道斷面設計、變截面隧道施工等問題提供借鑒。

如圖1所示，考慮4種斷面形狀的隧道環向開挖面： ①橢圓形斷面(B/D=0.5)； ②圓形斷面(B/D=1)； ③橢圓形斷面(B/D=1.5)； ④馬蹄形斷面(上部半圓，下部矩形)。

由于隧道形狀和計算模型受力均沿豎向中心線對稱，故選取右側一半開展計算分析(以圓形為例，如圖2所示)。

本文考慮的非均質黏土地層為不排水抗剪強度變化規律僅沿深度線性增長的情況，參考文獻[9]定義如下：

cu(y)=cu0+ρ(C+D/2-y)。

(1)

式中：cu(y)為豎向坐標y處對應的土體不排水抗剪強度；ρ為抗剪強度隨深度變化率。

1.2 變量分析與參數確定

如圖2所示，以圓形斷面為例，隧道環向開挖面穩定性分析涉及7個計算參數：{C,D,σs,σt,cu0,γ,ρ}。將這些參數無量綱化[9]： 埋深比C/D，臨界荷載比(σs-σt)/cu0，無量綱土體重度γD/cu0，無量綱土體非均質強度ρD/cu0。其中，臨界荷載比(σs-σt)/cu0為待求解量，可視為穩定性評價指標。為方便計算，將(σs-σt)/cu0中的σt置為0；對于σt不為0的情況，將σs或σt具體值代入(σs-σt)/cu0(即σt=0時的σs/cu0值)，即可獲得對應的σt或σs。計算參數取值見表1。

表1 隧道環向開挖面穩定性計算參數

1.3 運動單元上限有限元模型

以圓形隧道為例，非均質黏土地層隧道環向開挖面穩定性剛體平動運動單元上限有限元模型如圖2所示。其中，地層不排水抗剪強度沿深度線性增長。上限有限元模型初始三角形網格類似文獻[10]，本文未顯示。利用初始網格可獲取σs/cu0初始解，再進行一系列網格更新操作得到最終較優解答。求解流程與文獻[9]一致，即依據上限定理將問題轉換為非線性規劃，在一系列約束條件下搜索σs/cu0最優(最小)上限解。

上限有限元模型非線性規劃目標函數表達式為：

(2)

2 臨界荷載比σs/cu0上限解

2.1 臨界荷載比σs/cu0上限解結果

采用自編剛體平動運動單元上限有限元程序[9]，對表1無量綱參數的多種組合工況進行系列計算，獲得4種斷面形狀的隧道于不同參數影響下的環向開挖面穩定性臨界荷載比σs/cu0上限解，如表2所示。

2.2 臨界荷載比σs/cu0上限解對比

為檢驗程序計算結果的可靠性，選取非均質地層圓形隧道環向開挖面穩定性臨界荷載比σs/cu0上限解(ρD/cu0=0.5,γD/cu0=0～2)，與黃茂松等[4]、Wilson等[11]求解對比如圖3所示。

極限分析上限解精度與破壞模式合理性直接關聯[4]，UBRB(Wilson等[11])和多塊體上限法(黃茂松等[4])需預先構造破壞機制進行上限分析，通常應用于求解圓形隧道時臨界荷載比σs/cu0上限解的數值稍大。如圖3所示，本文σs/cu0上限解曲線處于FEM(Wilson等[11])方法上、下限解之間，與LBFEM方法下限解平均相對差值僅2%左右，說明程序可靠性良好。

2.3 臨界荷載比σs/cu0上限解的分析與討論

以圓形隧道為例，選取典型數據繪制σs/cu0上限解曲線如圖4所示，以此分析無量綱土體重度γD/cu0、無量綱土體非均質強度ρD/cu0以及埋深比C/D對于隧道環向開挖面穩定性的影響規律。

由圖4(a)可知，C/D=1時σs/cu0上限解與γD/cu0呈顯著的線性負相關，γD/cu0增加，地層穩定性隨之下降；ρD/cu0=0.05～1.0對應曲線由下至上依次等間距平行排列，說明σs/cu0上限解與ρD/cu0呈線性正相關。圖4(b)和4(c)分別對應于C/D=3、5的情況，此時σs/cu0上限解曲線規律與圖4(a)相同，而對應的σs/cu0上限解數值增加。

表2 隧道環向開挖面穩定性臨界荷載比σs/cu0上限解

圖3 圓形隧道環向開挖面穩定性臨界荷載比σs/cu0上限解對比 (ρD/cu0=0.5)

由圖4(d)可知，ρD/cu0=0.05時，σs/cu0上限解與C/D呈近似線性的增加或減小規律，其中γD/cu0較大時近似線性負相關，而γD/cu0較小時近似線性正相關。圖4(e)和4(f)分別對應于ρD/cu0=0.5、1.0的情況，此時σs/cu0上限解與C/D近似呈線性正相關規律。

由圖4(g)可知，γD/cu0=0時，隨著ρD/cu0增大，σs/cu0上限解表現為線性增加[12]。比較圖4(g)—(i)發現，無論γD/cu0取何值，C/D較大時σs/cu0上限解增長速度更加明顯。圖4(h)—(i)各曲線近似交于一點，說明存在一種特定狀態，σs/cu0上限解隨C/D的發展趨勢在此刻發生轉折。

(a) C/D=1 (b) C/D=3 (c) C/D=5

(d)ρD/cu0=0.05 (e)ρD/cu0=0.50 (f)ρD/cu0=1.00

(g)γD/cu0=0 (h)γD/cu0=2 (i)γD/cu0=3

其他條件相同時，對比表2中不同斷面形狀隧道的臨界荷載比σs/cu0上限解，發現隧道環向開挖面穩定性規律如下： 橢圓形隧道和圓形隧道臨界荷載比σs/cu0上限解由大到小依次為①橢圓形隧道(B/D=0.5)、②圓形隧道、③橢圓形隧道(B/D=1.5)； 其中，①橢圓形隧道較③橢圓形隧道臨界荷載比σs/cu0上限解高出15%～40%，且隨著埋深比C/D增加，差值進一步增大，反映出隧道跨度是影響穩定性的關鍵因素之一。此外，與圓形隧道輪廓接近的④馬蹄形隧道，其臨界荷載比σs/cu0上限解普遍稍小一些。

3 隧道環向開挖面破壞模式

3.1 破壞形態對比

圖5示出C/D=4、ρD/cu0=1、γD/cu0=0時2種極限分析上限法得到的圓形隧道環向開挖面潛在破壞模式。其中，圖5(a)示出UBFEM(Wilson等[11])方法對應的破壞形態，以速度矢量形式展示。圖5(b)示出本文剛體平動運動單元上限有限元計算得到的滑移線網破壞模式。需要說明的是，這里的滑移線即為模型里面活動的或者有效的速度間斷線，與主應力的跡線關聯。而模型中的大量滑移線通常表現為網狀的形式，可等效為破壞區域中的塑性區。

(a) UBFEM法 (b) 本文方法 (c) 速度矢量圖

由圖5可知，2種上限法展示的隧道環向開挖面破壞形態相似，表現為： 圖5(a)洞頂正上方速度矢量一致的規則區、拱頂左上方速度矢量過渡區，分別對應于圖5(b)洞頂正上方整體下沉區及右拱腰上方網狀滑移區。此外，圖5(c)為速度矢量閉合圖，以模型中最大數值V為參考，對速度矢量值進行了歸一化處理。速度矢量閉合圖中每條線段對應于圖5(b)中特定的有效間斷線，而圖中右上方交點與下方任意交叉點的連線，反映了圖5(b)中特定剛性單元的速度矢量。

由圖5(b)可知，本文圓形隧道環向開挖面的失穩破壞形態，與太沙基隧道土壓力理論假定的塑性剪切帶和扇形塑性變形區的破壞形態基本吻合(此時內摩擦角為0)。圖5(b)展示的潛在失穩破壞模式，可提供破壞區域界限及面積、滑移面角度等信息，為破壞機制的量化研究提供數據。

3.2 破壞模式演變規律分析

選取圓形隧道，繪制典型參數對應的環向開挖面潛在破壞模式如圖6所示，以分析埋深比C/D、無量綱土體重度γD/cu0和無量綱土體非均質強度ρD/cu0等參數對破壞模式演變規律的影響。

圖6(a)示出C/D=2、γD/cu0=3、ρD/cu0=0.5時的滑移線網破壞模式，其形態同圖5(b)一致。即上方為整體下沉規則區，外側滑動面豎直；右拱腰附近為網狀滑移區，其2簇滑動面交線近似垂直，而外側滑動面構成平滑曲線過渡。將相同條件下γD/cu0=0～3對應的主要滑動面疊加繪制如圖6(b)所示，發現γD/cu0增大時，破壞區域面積有所增加，且稍有下移趨勢。

圖6(c)示出C/D=2、γD/cu0=1、ρD/cu0=0.5時的滑移線網破壞模式。將相同條件下ρD/cu0=0.05～1.0對應的主要滑動面疊加繪制如圖6(d)所示，發現隨著ρD/cu0增大，破壞區域面積反而明顯減小，破壞區域有向上和向中心線靠攏的趨勢，這符合土體強度增加、破壞范圍隨之減小的規律。

圖6(e)示出C/D=3、γD/cu0=1、ρD/cu0=0.5時的滑移線網破壞模式。將相同條件下C/D=1～3對應的主要滑動面疊加繪制如圖6(f)所示，發現隨著C/D增大，破壞區域面積顯著增加，破壞區域有向下和向外擴張的趨勢，說明對于非均質純黏土地層(內摩擦角為0)，隧道埋深是影響破壞范圍的關鍵因素之一。

選取C/D=2、ρD/cu0=0.5、γD/cu0=1時，繪制不同斷面形狀對應的隧道環向開挖面潛在破壞模式如圖7所示。由圖7(a)可知，直墻馬蹄形隧道破壞范圍易于擴展至墻角處； 而圖7(b)顯示的跨度較大的橢圓形隧道(B/D=1.5)，其潛在破壞范圍反而局限于隧道拱部區域。在相同參數條件下，將4種不同斷面形狀隧道對應的主要滑動面疊加繪制如圖7(c)所示，可以看出水平方向破壞區域范圍由大到小排序為： ④直墻馬蹄形隧道>①橢圓形隧道(B/D=0.5)>②圓形隧道>③橢圓形隧道(B/D=1.5)。

實際上，破壞區域的范圍僅反映了非均質黏土地層隧道環向開挖面破壞特征的一個方面，其他方面如破壞狀態下的速度場也是關鍵要素。圖7(d)和7(e)分別示出相同參數條件下圓形隧道(B/D=1.0)和橢圓形隧道(B/D=1.5)滑移網破壞模式對應的速度矢量閉合圖。圖中均以最大矢量值進行歸一化處理，可以看出圓形隧道對應的速度矢量值普遍大一些，特別是水平方向的速度，由此帶來臨界狀態下耗散能的增加，反襯出跨度較大的橢圓形隧道(B/D=1.5)盡管破壞范圍小，但速度矢量值稍小，系統耗散能降低，反而不利于地層穩定。

(a) (b) (c) (d) (e) (f)

(a) (b) (c) (d) (e)

3.3 破壞區域面積定量分析

破壞范圍的定量分析，結合對應速度場的特征規律，可為地層加固范圍的確定提供數據支撐。為進一步說明隧道失穩極限狀態下破壞區域面積隨綜合因素的演變規律，定義無量綱參數Sf/Sc為相對破壞范圍，Sf為破壞區域面積(篩選模型速度不為0的剛體單元，累計單元面積之和)，Sc為隧道輪廓面積，Sf和Sc均考慮右側的半模型。以圓形隧道為例，不同參數條件下的相對破壞范圍Sf/Sc計算值見表3，選取C/D=2、4的數據繪制曲線如圖8所示。

由圖8和表3可以看出，圓形隧道相對破壞范圍Sf/Sc隨著C/D和γD/cu0的增加而增大，并隨著ρD/cu0的增加而減小，這直接以定量數據反映了圖6所示的演變規律。

此外，埋深比C/D變化時，Sf/Sc隨之變化非常明顯。如C/D=5較之C/D=1的破壞范圍普遍增長了10倍以上，說明隧道埋深較大時，受地表超載作用，破壞面水平向延伸范圍較大； 同時，對于不排水抗剪強度沿深度線性增長的非均質黏土地層(內摩擦角為0)，潛在破壞面易于擴展至地表，于是Sf/Sc呈現出與(C/D)2近似線性相關的特征。

Sf/Sc隨γD/cu0增加而線性增長，說明自重使得隧道附近土體更易于產生剪切滑移破壞；不過該增長比較緩慢，尤其是ρD/cu0較大時。由圖8可知，Sf/Sc隨著ρD/cu0增加而非線性減小，即土體抗剪能力提高，破壞范圍受到極大的限制從而縮小。

4 σs/cu0上限解擬合公式

如圖4所示，環向開挖面穩定性臨界荷載比σs/cu0上限解與γD/cu0、ρD/cu0呈線性關系，且線性斜率值和常數項僅與埋深關聯，故C/D固定時σs/cu0上限解可認為是無重度均質土上限解在γD/cu0、ρD/cu02變量作用下的線性變化。于是，σs/cu0上限解可寫為：

σs/cu0=Bγ·γD/cu0+Bρ·ρD/cu0+B0。

(3)

式中：Bγ為σs/cu0上限解與γD/cu0線性關系斜率；Bρ為σs/cu0上限解與ρD/cu0線性關系斜率；B0為無重度均質土σs/cu0上限解。

表3 圓形隧道相對破壞范圍Sf/Sc計算值

以圓形隧道為例，對不同埋深對應的σs/cu0上限解公式擬合如下：

(4)

可以看出Bγ、Bρ、B0隨C/D變化規律性較強，數據擬合可得：

Bγ=-1.034·C/D-0.141；

(5)

Bρ=4.526·C/D-3.154；

(6)

B0=-0.096·(C/D)2+1.229·C/D+1.372。

(7)

故圓形隧道σs/cu0上限解最終擬合公式為：

σs/cu0=(4.526·C/D-3.154)·ρD/cu0+(-1.034·C/D-0.141)·γD/cu0-0.096·(C/D)2+1.229·C/D+1.372。

(8)

同樣地，擬合其余3種斷面隧道σs/cu0上限解公式見表4。其確定系數R2均大于0.99，表明擬合精度較好。上限解擬合公式只需少量參數(C,D,γ,ρ,cu0)即可初步評估抗剪強度隨深度線性增長的純黏土地層隧道環向開挖面穩定性，方便工程應用。

(a) C/D=2

(b) C/D=4

5 試驗驗證

為進一步檢驗剛體平動運動單元上限有限元方法的適用性，選取基于離心機模型試驗[13]的均質軟黏土地層單孔隧道臨界荷載比σs/cu0結果進行對比驗證。其試驗采用直徑D=100 mm的圓形隧道模型，土體重度γ=18.1 kN/m3，均質地層中抗剪強度隨深度變化率ρ=0。試驗過程中控制土體孔隙水壓力變化不超過1%以保證不排水狀態。試驗考慮不同埋深C、不同抗剪強度cu0條件下隧道環向開挖面失穩破壞，臨界荷載比σs/cu0試驗結果與本文運動單元上限解及擬合公式上限解對比如表5所示。

表4 4種斷面隧道環向開挖面穩定性臨界荷載比σs/cu0上限解擬合公式

表5 臨界荷載比σs/cu0試驗結果與本文上限解對比

由表5臨界荷載比σs/cu0數據對比可知，盡管數值上存在不同程度的差異，剛體平動運動單元上限解與離心模型試驗結果總體上吻合良好，從而驗證了本文計算方法的可靠性和計算結果的精度； 同時，擬合公式對應的σs/cu0上限解也具有較好的一致性。

6 結論與討論

1)在抗剪強度沿深度線性增加的非均質黏土地層條件下，不同斷面形狀隧道環向開挖面穩定性差異顯著，橢圓形隧道(B/D=0.5)較之橢圓形隧道(B/D=1.5)臨界荷載比σs/cu0上限解高出15%～40%，且埋深增加將增大這一差值。σs/cu0上限解與ρD/cu0正相關而與γD/cu0負相關，ρD/cu0和γD/cu0均影響C/D對上限解的作用趨勢。

2)由隧道環向開挖面潛在破壞模式形態演變及范圍定量分析發現，破壞范圍主要取決于斷面形狀和埋深大小。較之圓形隧道，直墻馬蹄形隧道失穩破壞區域易擴展至邊墻底部。此外，埋深成為影響隧道失穩破壞范圍的關鍵因素。

3)基于剛體平動運動單元上限有限元數據，獲得了非均質黏土地層不同形狀隧道的臨界荷載比σs/cu0上限解擬合公式。對比已有離心機模型試驗結果，驗證了運動單元上限解及擬合公式的可靠性，可為工程應用提供參考。

4)考慮到現階段分析模型與工程實際比照存在不少簡化，下一步的研究應考慮荷載的非均勻分布以及地層成層條件等因素的影響。