含氣地層不均勻沉降對盾構隧道結構影響的固-液-氣耦合分析

王暉皓, 葉 斌*, 涂新斌, 錢玉華, 程子睿

(1.同濟大學土木工程學院, 上海 200092； 2.國家電網有限公司, 北京 100031；3.江蘇省送變電有限公司, 南京 211100)

在中國沿海地區的地下工程施工過程中,經常發現富含高壓氣體的地層。這些氣體的主要成分為甲烷,且埋深較淺,被稱為淺層沼氣。例如,杭州地鐵1號線穿越了含淺層沼氣的地層,為此在施工過程中采取了專門的應對措施[1]。淺層沼氣在軟土中多以氣囊或氣層的形式存在。地下工程開挖將引起土體擾動,進而打破水-土-氣的三相平衡狀態,導致淺層沼氣的滲流或排出。在氣體滲流或排出過程中,土體骨架、地下水和淺層氣三相間將發生相互作用,引起土體變形甚至產生不均勻沉降,有可能對隧道結構的安全性造成威脅[2]。因此,為了保證隧道的安全運營,有必要研究含氣地層的不均勻沉降對盾構隧道結構的影響。

由于地層條件的復雜性,數值模擬方法通常被用于研究隧道的結構穩定性問題。近年來，考慮多相耦合的有限元方法逐漸興起[3]，該方法已被廣泛應用于油氣開采[4]、二氧化碳地下封存、降雨入滲[5]及核廢料處置庫存的安全性等領域。Lewis等[6]提出了熱-固-液-氣耦合模型模擬儲層變形和油氣開采過程中的兩相流動和地面沉降。Nagel等[7]基于BBM(barcelona basic model)模型建立了三相耦合有限元模型,并對壓氣法隧道施工過程進行了模擬。Rutqvist等[8]將FLAC3D與TOUGH軟件相結合,開發了FLAC-TOUGH耦合程序,實現了熱-流-固(THM)多相耦合的數值模擬。

針對淺層氣體滲流及排出引起的地層沉降問題,許多學者也開展了相關研究工作[9]。鐘方杰[10]針對杭州灣大橋工程地質勘探過程中發生的多次強烈井噴現象,應用FLAC2D兩相流模塊,模擬在含淺層氣的地層中進行鉆探引發淺層氣釋放對地層的影響效應。王勇[1]在此基礎上,以杭州地鐵典型含氣區段——錢塘江南岸某地層為背景,應用FLAC2D中的兩相流模塊,提出了放氣路徑下儲氣砂土的本構模型,并對工程放氣路徑(排氣井)下儲氣砂層變形進行了數值模擬。盧浩等[11]認為氣體泄漏引起的地層不均勻沉降將會導致管片受力變形，通過對沉降部分土體進行開挖來模擬氣體釋放對地層的影響,從而避開了兩相流的難點,模擬了氣體釋放對管片受力變形性能的影響。彭海波[12]在Rutqvist等[13]研究的基礎上開發了FLAC-TOUGH接口程序,控制TOUGH2和FLAC3D在數值模擬過程中有序運行并利用耦合關系式傳遞相關狀態變量,克服了FLAC3D不能進行兩相流分析計算的缺點,實現了對水-氣二相流-固耦合模擬。李曉敏[14]利用TOUGH2和FLAC3D建立水-氣-固三相耦合模型實現了對非飽和土邊坡的降雨入滲模擬。可以發現,目前的研究工作大多基于已有的商用軟件進行模擬,但這些軟件大多只考慮兩相的作用[15],無法模擬固-液-氣三相的共同耦合作用。此外,已有研究的土體本構模型多采用軟件自帶的簡單本構關系(如摩爾-庫侖模型等),無法模擬固-液-氣三相共存(非飽和狀態)土體的復雜力學行為特征。

以蘇通氣體絕緣金屬封閉輸電線路(gas-insulated metal enclosed transmission line,GIL)綜合管廊項目中盾構隧道穿越的長江江底淺層氣體分布區為研究對象,采用自主開發的水-土-氣三相滲流-變形耦合有限元程序,并基于Zhang等[16-17]提出的本構模型,模擬了含淺層沼氣地層的變形特性,研究了氣體滲流及排出引起的地層不均勻沉降對隧道結構受力的影響。

1 計算模型

1.1 工程背景

蘇通GIL綜合管廊工程位于G15沈海高速蘇通長江大橋上游1 km,是“淮南-南京-上海1 000千伏交流特高壓交流工程”下穿長江的重要節點,由兩岸引接站和GIL管廊組成。其中管廊部分采用盾構法施工穿越長江,隧道所在位置如圖1所示。

圖1 隧道所在位置Fig.1 The location of tunnel

在中國長江中下游及東南沿海地區,沿江、沿海的湖相、河流相及海相的沉積平原中廣泛分布有第四系淺層生物氣藏,其埋深一般小于100 m。根據《淮南-南京-上海1 000千伏交流特高壓輸變電工程蘇通GIL綜合管廊工程施工圖有害氣體勘察報告》,GIL綜合管廊工程的盾構隧道有部分區段位于淺層氣體分布區,其推定里程為DK1+0～DK1+780,如圖2所示。地層中存在這些淺層氣的成分甲烷(CH4)占比(85%～88%)、氮氣(N2)占比(8%～10%)、氧氣(O2)占比(2%～3%),為生物成因淺地層天然氣,成團塊狀、囊狀局部集聚分布,氣體壓力約為0.4～0.9 MPa,氣藏類型為原生超淺層微氣藏,有害氣體在土層中呈扁豆體狀、團塊狀、囊狀、蜂窩狀不連續分布,氣、水同層,在淤泥質、粉質黏土層中以團狀氣包的形式存在。隧道開挖將導致原有氣囊失衡,多個沼氣氣囊出現運移,有可能形成連續的含氣地層。根據地質勘察資料,分布區(DK1+0～DK1+780)地層由上至下依次分布有：③淤泥質粉質黏土層、④1粉質黏土混粉土層、④2粉砂層、⑤1粉細砂層、⑤2粉砂層、⑥1中粗砂層、⑦1粉細砂層、⑧1中粗砂層及⑧2粉細砂層。地層頂部距離自由海平面約7.0 m,隧道標高為-40.0～-51.6 m,中心標高為-48.5 m。

圖2 有害氣體分布區Fig.2 Harmful gas distribution aera

1.2 有限元模型

為了實現含氣地層在三維地層模型中固-液-氣三相耦合數值模擬計算,對上述地層進行概化并做如下假設簡化:①土層分布均勻且水平分布,氣壓分布均勻且各向同性；②地下水位取地表高程；③含氣地層位置緊貼隧道底面。

實際場地中所涉及的地層數較多,但部分土層物理性質相近,參數相差不大。因此,為簡化模擬計算步驟,提高模擬計算效率,將實際地層數劃分為2類,其中,③淤泥質粉質黏土層、④1粉質黏土混粉土和④2粉砂層劃分為第Ⅰ層粉土層,剩余下方地層統一劃分為第Ⅱ層細砂層。

簡化后的模型斷面如圖3所示,建模范圍取隧道半徑的10倍以上以消除邊界效應的影響。整個計算分析模型的尺寸的長度方向取131.6 m,寬度方向取100 m,深度方向取105.8 m。隧道外徑11.6 m,內徑10.5 m,襯砌厚度為0.55 m。由于現場無法獲取具體氣藏分布尺寸資料,故假設含氣地層分布長度取模型長度的1/2為60.7 m,分布寬度取隧道外徑長度為11.5 m,而氣藏厚度根據勘察報告取15 m。隧道上方布設觀察點(U1、U2、U3),隧道下方布設觀察點(D1、D2、D3、D4、D5、D6)和觀察單元(E1、E2、E3)。圖4為三維模型示意圖,圖5為建立的三維有限元網格模型。

圖3 基本地層二維模型及觀察點布置Fig.3 Basic soil layer 2D model and the observation points arrangement

圖4 三維模型示意圖Fig.4 3D model schematic diagram

圖5 三維有限元網格Fig.5 3D finite element method mesh

1.3 模型初始條件

滲流場的初始壓力條件主要包括初始孔隙水壓與初始氣壓兩部分。對于孔隙水壓的初始條件,根據實際場地地質鉆孔剖面圖,設定地層頂部為自由水平面,即自由水面的位置水頭為40 m。非含氣層初始水壓取自由水面下靜水壓力,則非含氣層中某點的初始孔隙水壓(pw0)表達式為

pw0=ρwg(40-z)

(1)

式(1)中：ρw為水的密度；g為重力加速度；z為該點的豎坐標。

對于含氣層,結合勘察報告靜力觸探CPTU試驗中的部分鉆孔孔隙水壓消散曲線(圖6),初始孔隙水壓不含毛細水壓力,則含氣層的初始孔隙水壓p′w0表達式為

圖6 孔隙水壓消散曲線Fig.6 Pore water pressure dissipation curve

p′w0=ρwg(10-z)

(2)

對于氣體壓力初始條件,假設氣體在含氣層以及滲流通道內各處氣壓連通。根據有害氣體勘察報告,含氣地層初始氣藏壓力取900 kPa,初始的水相飽和度設為0.8,非含氣層初始的水相飽和度設為最大飽和度,其中粉土層為1.0,細砂層為0.98。為了方便描述,文中飽和度均指水相飽和度。

對于土體應力場,假設隧道開挖前土體處于自重應力場作用。為簡化考慮隧道開挖的影響,數值模擬過程分3步進行:①計算完整土體在自重作用下的初始應力場；②去除隧道內的土體單元模擬隧道開挖,并在隧道邊界上設置薄板單元模擬襯砌添加,得到隧道開挖后應力場；③將②中得到的土體應力場作為滲流計算的初始應力場,計算沼氣滲漏模式下固-液-氣三相之間的相互作用。

1.4 模型邊界條件

模型上邊界設為自由透水透氣邊界,無位移約束；下邊界設為不透水不透氣邊界并約束z向位移；

左、右側邊界設為不透水不透氣邊界并約束x向位移；前、后側邊界也設為不透水不透氣邊界并約束y向位移；隧道邊界除滲氣口外皆設為不透水不透氣邊界；其余內部邊界均為滲透邊界。D1為滲氣邊界中心,如圖7所示。滲氣邊界為恒定氣壓邊界。

圖7 滲氣口附近的局部放大圖Fig.7 Partial enlargement of the vicinity of the gas permeation port

1.5 模型參數

數值計算采用自開發的有限元計算程序。該程序能夠對固-液-氣三相滲流-變形進行耦合分析,目前已經在多個巖土工程問題中得到了應用[18-20]。在三相耦合計算中,土層本構關系采用前章介紹的基于修正劍橋的非飽和彈塑性本構模型[16-17]。

模型參數包括基本地層和隧道襯砌材料的物理參數、水力學參數、土水特征曲線參數和劍橋模型基本參數,具體取值如表1～表3所示。地層物理及水力學參數取值主要依據項目的土工試驗報告,其中常規參數(如重度、泊松比等)可以直接采用,其他參數根據規范標準進行換算；土水特征曲線參數參考相似性質的土體參數進行取值；劍橋模型基本參數部分取自土工試驗報告,其他參數采用相同土類土體的經驗值；盾構隧道的襯砌材料按照項目工程結構設計文件選用的混凝土規格確定。

表1 基本地層及隧道襯砌材料的物理及水力學參數Table 1 Physical mechanics and hydraulic parameters of basic soil layer and tunnel lining

表2 土水特征曲線參數取值Table 2 Parameter value of soil water characteristic curve

表3 劍橋模型基本參數Table 3 Parameter value of Cambridge model

2 計算結果分析

盾構隧道在運營過程中容易發生過量的縱向沉降或不均勻沉降,引發隧道結構局部破壞或縱向扭曲變形,影響隧道的正常運營。從隧道豎向位移和截面附加彎矩兩方面評價隧道的穩定性,并且通過滲氣口飽和度與氣壓的變化觀察氣體的滲漏情況。

2.1 豎向位移分析

圖8為5年后地層最終豎向位移結果。圖8中位移正值表示隆起,位移負值表示沉降。從圖8可以看出,最大沉降發生在滲氣口(D1)位置兩側,其最大沉降值為26.0 mm。滲氣口處的最終沉降位移約為10.5 mm。其沉降位移比兩側位移小的原因是滲氣口下方由于氣體快速流出產生的向上滲流力作用比兩側土體更大所導致。該點處沉降隨時間變化的趨勢如圖9(a)所示。從圖9(a)可以看出,滲流開始后滲氣口的沉降位移迅速增大,而后產生少量回彈并逐漸穩定。圖9(b)為圖9(a)前期的局部放大。通過圖9(b)可以看出,沉降發生的時間主要在前50 h內,后續產生少量回彈直至達到基本穩定。回彈產生的原因是由于含氣地層的氣體快速滲漏引起整個地層的初期沉降,而后隨著地下水的涌入以及地層的卸荷效應產生回彈并逐漸達到平衡。

圖8 最終豎向位移Fig.8 Final vertical displacement

圖9 滲氣口豎向位移隨時間變化Fig.9 Distribution of the vertical displacement with time at the gas permeation port

圖10(a)、圖10(b)分別為隧道上部和下部觀察點的豎向位移隨時間的變化。其中,隧道上部的沉降主要發生在滲氣開始后的前50 h內,而后產生少量回彈直到穩定。在前50 h內,位于隧道頂部的觀察點(U3)沉降最大,而后3個觀察點均有所回彈；最終隧道上部所有觀察點的沉降穩定值相差較小,隧道頂部的觀察點(U3)沉降最大,最大值約為12.3 mm。而隧道下部的觀察點(D2～D5)的豎向位移開始有所波動,而后均有所回彈；其中位于含氣層底部的點回彈值最大,達到了15 mm。越接近含氣地層底部,回彈越明顯。引起回彈的原因可能是含氣地層下部的部分地層的卸荷效應。

圖10 觀察點豎向位移隨時間變化Fig.10 Distribution of the vertical displacement with time at the observation points

圖11為隧道中部(y=0)截面地表的沉降位移曲線。從圖11可以看出，地表沉降位移曲線呈漏斗狀,位于滲氣點正上方的地表點沉降最大,向兩側逐漸減小，其最大沉降為7.5 mm。

圖11 地表豎向位移曲線分布Fig.11 Distribution of the surface vertical displacement

2.2 截面附加彎矩分析

圖12為隧道中部(y=0)截面最終附加彎矩。截面附加彎矩是由排氣引起的彎矩。由圖12可知,滲氣口的位置兩側,截面附加彎矩最大,附加彎矩的最大值為41.7 kN·m/m；而越遠離含氣層的位置,截面附加彎矩越小。

根據蘇通GIL工程盾構隧道的結構設計文件,隧道截面彎矩的設計值為1 708.1 kN·m/m,因此可知附加彎矩的最大值僅占設計值的2.4%。因此,可以認為含氣地層的沉降對襯砌結構的受力影響較小。

2.3 飽和度與氣壓分析

最終飽和度分布如圖13所示。從圖13可以看出,整體含氣層飽和度最小值為0.97。滲氣口單元的飽和度隨時間變化如圖14所示,圖14(b)為圖14(a)前期的局部放大。從圖14可以看出,飽和度逐漸增大到0.98。前期飽和度快速增加的原因是由于氣體的迅速排出,土體飽和度增大,氣壓迅速下降,水壓隨之上升,導致基質吸力減小。

圖13 最終飽和度分布Fig.13 Distribution of final saturation

圖14 滲氣口飽和度隨時間的變化Fig.14 Distribution of saturation with time at the gas permeation port

滲氣口及附近單元氣壓隨時間變化如圖15所示,圖15(b)為圖15(a)前期的局部放大圖。從圖15可以看出,所有觀察單元的滲氣壓力都在初期迅速下降,而后緩慢下降直到穩定。其中,滲氣口單元(E1)的最終氣壓為436 kPa,與該單元靜水壓(450 kPa)相差較小；而滲氣口下方單元(E2、E3)的最終氣壓穩定在500 kPa左右，且深度越深,最終氣壓值也越大。

圖15 滲氣口及附近單元氣壓隨時間的變化Fig.15 Distribution of gas pressure with time at the gas permeation port and nearby elements

3 結論

以蘇通GIL綜合管廊項目中盾構隧道穿越的長江江底淺層氣體分布區為研究對象,采用自主開發的水-土-氣三相滲流-變形耦合有限元程序,并基于文獻[13-14]提出的本構模型模擬了含淺層沼氣地層的變形特性,研究了氣體滲流引起的地層不均勻沉降及其對隧道結構受力的影響,得出如下結論。

(1)根據位移計算結果,最大沉降發生在滲氣口位置兩側,約為26.0 mm。隧道上部土體的沉降值先增大后由于回彈而有所減小,且越靠近隧道結構的土體的最終沉降值越大；隧道下部土體的沉降值整體上先隨時間緩慢增大后產生回彈,且越接近含氣地層底部邊界的點回彈越大。而地表越靠近中心位置(滲氣口上方),對應的豎向沉降位移值越大,整體呈漏斗狀。

(2)根據截面附加彎矩的計算結果,滲氣口兩側的截面附加彎矩最大,而越遠離含氣層的位置,截面附加彎矩越小。附加彎矩的最大值為41.7 kN·m/m,僅占設計值的2.4%。因此,含氣地層的沉降對隧道結構的受力影響較小。

根據飽和度和氣壓變化的計算結果,氣體滲漏過程完成后,含氣層接近飽和。其中,滲氣口前期飽和度快速增加的原因是由于氣體的快速排出,土體飽和度上升,氣壓迅速下降,水壓隨之上升,導致基質吸力減小。滲氣口最終飽和度為0.98,最終氣壓與該處靜水壓相差較小。

受限于目前的勘查技術水平,僅針對氣體連續分布的典型氣藏模型條件下的隧道變形受力情況進行了分析,然而實際地層中的氣藏由于河流沉積作用或是松散堆積物的孔隙不同等原因并非單一均勻分布,其大小、位置以及形狀都有很大的不確定性。為了更準確地模擬滲流場的變化過程,一種方法是建立精細化的氣藏地質模型,這有賴于對土層和氣藏信息的精細勘探技術；另一種方法是采用隨機生成的地質模型描述分散分布的氣藏,這將是筆者進一步探討和研究的方向。