基于GDEM的隱伏巖溶隧道隔水巖體水壓致裂安全厚度及破裂演化規律分析

李 濤, 張 麗, 蔣 慶, 馮 春, 趙 然

(1. 山東高速濟萊城際公路有限公司， 山東 濟南 250014； 2. 中國科學院力學研究所， 北京 100190)

0 引言

在富水巖溶地區，隧道與溶洞間的巖層由于開挖引起的擾動而極易失穩，引發突水、涌泥災害，是隧道安全施工的重大問題，也是目前工程建設研究的重點之一。目前研究巖溶隧道隔水巖層穩定性的方法主要有模型試驗、理論研究和數值模擬等。

模型試驗需根據原型確定試驗模型的幾何相似比，然后根據力學方程和邊界條件確定應力相似比，并采用相似材料模擬一定工況下的突水現象。魏星等[1]以北崗隧道DK462+527～+537段為原型，采用1∶100幾何相似比模型，模擬了圍巖出現位移劇增、突泥涌水和破壞的現象，但材料的力學指標未能達到相似性要求。趙瑜等[2]建立了應力相似比為40、幾何相似比為55的某在建隧道模型，并與FLAC3D數值模擬結果進行對比，發現在各種趨勢上有較好的一致性，且物理模型的最終破壞及影響范圍與數值模擬有較好的一致性，但局限于材料特性和測量上的誤差，未能體現數值模擬的全部現象。李浪等[3]研制出一套深長隧道突水地質災害三維模型試驗系統，該系統不僅能夠模擬地應力場、水壓等初始環境，還可精確模擬隧洞開挖，但該模型試驗是在隧道模型簡化條件下進行的。Yang等[4]研發了高地應力高滲隧道突水模型試驗系統,并以重慶西馬隧道為研究對象，建立了三維流固耦合模型，通過監測應力、位移、滲流壓力等的變化趨勢，較為真實地模擬了突水發生的過程，但采用的是相似比材料和簡化模型。鑒于地層結構的復雜性和研究條件限制，模型試驗只能一定程度地模擬突水現象，對工程的指導性有限。大多數理論研究是通過對現場地質模型進行一系列簡化而得到相應的理論模型，對于簡化后的理論模型可以用結構力學、彈性力學和斷裂力學中的相關公式加以計算，從而得到臨界安全厚度[5]。李利平等[6]以斷裂力學、彈塑性力學理論為基礎，推導出隧道最小巖石保護厚度的半解析解表達式，并認為巖溶隧道裂隙突水具有明顯的時空效應特征。郭佳奇等[7-8]針對側壁溶腔，以塑性區和高滲透帶貫通與否作為判定中間巖柱最小安全厚度的標準，建立最小防突厚度計算公式；基于彈性厚板理論，推導了邊界條件分別為固支和簡支2種模式的巖溶隧道掌子面巖墻安全防突厚度及臨界水壓計算公式。Yang等[9]基于Hoek-Brown非線性破壞準則，利用上界定理和變分原理推導獲得了隔水巖柱的安全厚度，并研究了各參數的影響，繪制出了突水破壞區域。地質工程的復雜性使各種理論研究和模型試驗缺乏有力的支持，相比之下，數值模擬方法具有較廣泛的適用性,成為解決巖土工程問題的有效工具。隨著計算機技術的迅猛發展, 各種數值計算方法越來越普遍地被應用到圍巖穩定性分析中[10-11]。

目前，隧道圍巖采用的數值計算方法主要為有限元法、有限差分法和離散元法等，ABAQUS、MIDASGTS等為有限元軟件；FLAC3D為有限差分軟件；UDEC、3DEC等為塊體離散元軟件。佘健等[12]針對筆架山隧道建立了三維有限元模型，采用彈塑性模型進行計算，通過監測地表、拱頂、拱腰等的位移變化，分析開挖的影響。謝海文等[13]利用二維離散元數值模擬軟件UDEC，分析了貴州省德江隧道圍巖中裂隙的響應特征、流固耦合效應和滲流場分布特征。黃明利等[14]利用巖石破裂過程分析程序(RFPA), 建立了二維平面應變模型, 對隧道施工誘發隱伏溶洞破裂突水過程中的應力場、位移場和聲發射等特征進行了系統研究。徐長金等[15]基于MIDAS/GTS有限元分析軟件，分析了溶洞半徑以及溶洞與隧道底板距離的變化對圍巖應力分布和塑性區范圍的影響。雷霆等[16]利用FLAC3D有限差分軟件進行三維隧道開挖模擬，分析了頂部溶洞與隧道掌子面附近圍巖塑性區的發展趨勢，以及不同溶洞水壓下隧道的安全厚度。李紅衛[17]采用數值分析軟件ABAQUS對貴陽軌道交通建設中的隧道底部溶洞圍巖穩定性進行三維數值模擬，分析了距隧道斷面不同間距下溶洞對隧道底部、頂部和地表位移的影響。Shan等[18]利用FLAC3D建立了三維隧道模型，對3種不同位置的溶洞開展了數值模擬，并采用經驗公式分析了各因素對安全厚度的影響程度，認為數值方法更具合理性。巖體的破裂包含小變形、損傷演化、裂紋萌生、裂紋擴展以及大位移、大轉動等階段，是一個連續到非連續的破壞過程。因此，隔水巖體的破壞是開挖擾動引起的裂隙擴展、貫通直至破裂的漸進過程[6]。有限元法及有限差分法能夠較好地模擬材料在連續狀態下的特性，但不能模擬材料從連續到非連續的過程及在非連續狀態下的運動特性；塊體離散元在模擬非連續體的運動特性方面具有一定的優勢，但較難模擬材料的連續變形過程。目前大多數研究以溶洞的空間二維分布形態為主，三維模型受軟件本身計算能力的影響，網格較大，精度較差，計算耗時多。因此，大多數值模擬軟件無法真實模擬三維隧道的開挖突水破壞過程。

從上述分析可見，模型試驗、理論研究對于描述突水這種復雜地質體的力學行為具有一定局限性，而傳統數值模擬方法在描述巖體漸進破壞過程和計算精度方面均有不足。本文選取GDEM系列中的DAS軟件，建立百萬單元三維溶洞隧道模型，分析在高地應力、高溶洞水壓下，開挖擾動引起的巖體破壞演化規律；為了分析隔水巖體的漸進破壞過程，引入無量綱指標破裂度來定量化描述巖體的破裂狀態；通過監測掌子面位移變化特征，獲得隔水巖體的安全厚度，并分析在這一過程中隔水巖柱的破裂發展趨勢與位移之間的對應關系。

1 基于GPU加速的CDEM算法

GDEM應力分析系統(GDEM-DAS)是基于連續-非連續介質力學的離散元方法CDEM(continuum-based discrete element method)[19-21]的高性能有限元-離散元計算軟件。該軟件以CDEM理論為基礎，利用GPU(顯卡)進行加速計算。CDEM將有限元與離散元進行耦合，在塊體內部進行有限元計算，在塊體邊界進行離散元計算，通過塊體內部及塊體邊界的斷裂，不僅可以模擬材料在連續狀態下及非連續狀態下的變形、運動特性，更可以實現材料由連續體到非連續體的漸進破壞過程，目前已廣泛應用于滑坡、爆破、地下工程等行業[22-24]。

1.1 CDEM算法的基本概念

CDEM中的數值模型由塊體及界面2部分構成，如圖1所示。塊體由1個或多個有限元單元組成，用于表征材料的彈性、塑性、損傷等連續特征；2個塊體間的公共邊界即為界面，用于表征材料的斷裂、滑移、碰撞等非連續特征。界面包含真實界面及虛擬界面2個概念，真實界面用于表征材料的交界面、斷層、節理等真實的不連續面，其強度參數與真實界面的參數一致；虛擬界面主要有2個作用，一是連接2個塊體，用于傳遞力學信息，二是為顯式裂紋的擴展提供潛在的通道(即裂紋可沿著任意一個虛擬界面進行擴展)。

(a) 數值模型 (b) 塊體 (c) 界面

1.2 離散特征

采用CDEM方法描述連續-非連續問題時，連續部分采用有限元離散，非連續部分添加界面接觸表征其力學特征，相鄰接觸面通過空間拓撲信息進行查找，在公共接觸面兩側塊體的相應節點間建立彈簧(包括法向彈簧和切向彈簧)，接觸面之間的相互作用力由彈簧力和彈簧的特征面積來表征。假設計算區域被結構面T1、T2和T3切割后形成3個區域，各個區域內需要進行有限元離散，接觸邊界也需要進行離散化并建立界面接觸來表征結構面，如圖2所示。

(a) 塊體單元 (b) 接觸彈簧

1.3 破裂判斷準則

CDEM采用基于時程的動態松弛技術進行顯式迭代計算，因此可獲知每時步的單元應力和節點彈簧力，并根據破裂準則進行強度判斷，將發生破裂的節點彈簧力置為0，并記錄破裂彈簧的特征面積。本文擬在塊體部分采用線彈性模型，界面接觸彈簧采用張拉-壓剪復合準則判斷其破壞狀態。張拉模式的彈脆性破壞采用最大張力準則：

(1)

式中：σn為接觸法向應力；kn為法向彈簧剛度；u為彈簧的法向位移；σt為接觸抗拉強度。

當界面發生張拉破壞時，界面上的剪切應力也相應地置為0。剪切模式的彈脆性破裂采用摩爾-庫侖準則 ：

(2)

式中：τ為接觸切向應力；kτ為切向彈簧剛度；ν為彈簧的切向位移；c為接觸黏聚力；φ為內摩擦角。

當界面接觸發生破裂之后，接觸所連接的兩側塊體發生張開和滑動。

2 計算模型及材料參數

2.1 模型尺寸

根據趙明階等[25]的研究成果和彈塑性理論，隧道周邊影響范圍為3～5倍的洞徑，因此，建立三維隧道計算模型為100 m×100 m×100 m，三心圓隧道位于模型中間，開挖跨度為19 m，高度為11 m，襯砌厚度為55 cm，如圖3(a)所示。溶洞簡化為橢球模型，長軸為30 m，短軸為14 m，小于隧道跨度。考慮突水最不利位置，將溶洞置于隧道掘進的正前方，溶腔垂直隧道軸線方向(z向)的尺寸約為隧道高度的2倍，如圖3(b)所示。采用前處理功能較強的ANSYS建立幾何模型后，共劃分105萬個四面體單元、420萬個節點，然后通過數據轉換程序將節點單元導入到GDEM-DAS中，再進一步進行巖層屬性設定、開挖襯砌模擬和計算分析等。

(a) 三維溶洞隧道網格

2.2 邊界條件

計算模型設定為埋深130 m的深埋巖溶隧道，符合平面應變模型，因此，計算模型左右側邊界、前后邊界和底部均采用法向位移約束；模型上邊界(z=100 m)施加上端巖體自重；根據京滬高速某公路段地應力資料，圍巖側壓力系數取1.0，轉換為線性分布荷載分別施加于計算模型z=0 m、z=100 m、x=0 m和x=100 m 4個平面上。

鑒于溶腔軸的尺寸相對較小，可忽略溶腔內部上下水壓差別，水壓簡化為施加于橢球內表面的均布壓力。李利平等[6]研究發現瞬時突水突泥型溶洞水壓一般為1～3 MPa，本文擬分別計算水壓1、2、3 MPa時的隔水巖層破壞狀態。

2.3 材料參數

因地質情況復雜，各地巖石參數差異較大，本文圍巖參數選自京滬高速某公路段現場測試結果。該路段采用初期支護和二次襯砌支護，初期支護包含鋼筋網、鋼拱架、錨桿和混凝土等材料，二次襯砌包含C30模筑混凝土和鋼筋，各材料屬性均不同，故采用等效參數進行計算。圍巖和支護參數取值見表1。

表1 圍巖與支護參數

2.4計算過程

計算過程為先采用彈性模型計算獲得巖體初始應力狀態，再采用強度準則模擬出圍巖的初始破裂狀態，然后進行開挖過程計算： 采用全斷面開挖模型，開挖進尺2 m，在開挖下一進尺的同時完成上一進尺的襯砌支護，每次開挖計算均采用破壞準則迭代至40 000步。

選取掌子面為監測面，分別在每個監測面的拱頂、拱底、拱腰和中心位置設置監測點，如圖4所示。通過監測開挖過程中監測點的位移變化規律，分析開挖擾動對隔水巖層的影響。

(a) 監測面 (b) 監測點

2.5 隔水巖體安全厚度

從數值模擬的角度來講，位移能夠直觀地表征巖體的漸進破壞過程。在同一計算工況下，隨著掌子面向溶腔方向的推進，開挖引起隔水巖體出現裂縫，隨后逐漸擴展，此時掌子面的位移緩慢增加；當完全貫通時，掌子面的位移會迅速增加，存在突變現象，說明此時隔水巖體的厚度已不具備足夠的安全儲備。因此，可將位移值突變發生的前一開挖步視為安全厚度，在此時采取工程措施進行防治。

2.6 隔水巖體的破裂度

2.6.1 破裂度的定義

破裂是表征巖體災變過程的重要現象，隔水巖體的破裂程度是影響突水發生的重要因素。現有的測量手段可直接測得暴露的裂縫長度，但無法由此推測出巖體內部的損傷程度以及漸進破壞過程。

為了建立隔水巖體漸進破裂程度與開挖進程之間的關系，本文借助于連續-非連續數值計算方法，模型中塊體之間的界面即為潛在的破壞面，利用強度準則判斷界面狀態，通過彈簧斷裂模擬裂縫擴展，進而通過統計彈簧力為0的特征面積獲得當前狀態下巖體的破裂面積。將當前狀態下的破裂面積與界面總面積之比定義為破裂度，即破裂度

D=S/St×100%。

(3)

式中：S為研究區域內當前斷裂彈簧特征總面積；St為研究區域內界面總面積。

2.6.2 單軸壓縮破裂度計算和網格依賴性驗證

下面以單軸壓縮為例說明破裂度的計算方法，并進行網格依賴性驗證。計算模型為直徑50 cm、長100 cm的圓柱，上下端同時施加相同的加載速率0.5×10-9m/s，計算模型參數見表2，界面彈簧參數通過材料參數計算獲得。

表2 計算模型參數

計算模型采用張拉-壓剪復合準則進行計算，由于破裂發生在界面，通常選用隨機性較好的四面體網格進行計算。不同網格數量的計算模型如圖5所示。為了驗證計算結果的可靠性，對四面體網格進行加密，計算模型網格數量見表3。

(a) 模型1

表3 計算模型網格、接觸面和彈簧組數量

研究區域為整個圓柱，單元之間均采用界面接觸模型進行計算，通過彈簧連接，接觸面數量和彈簧組數見表3。采用1.3節的破裂準則進行強度計算，通過統計各監測時步彈簧力為0的特征面獲得當前時步的破裂總面積，彈簧力不為0的特征面積為當前時步的界面總面積，二者之比即為當前時步的破裂度。破裂度隨時步變化曲線如圖6所示。計算結果顯示，在同一計算工況下，不同密度網格的破裂度隨加載時步的變化趨勢基本一致，都存在突變點，在突變點后破裂度逐漸趨于穩定；相鄰網格數量的破裂度誤差小于4.3%。由此可見，破裂度對網格的依賴性較小。

圖6 破裂度隨時步變化曲線

3 隔水巖層破壞特征分析

3.1 隔水巖體監測點位移變化分析

以2 MPa溶洞水壓為例，計算模型在重力和側壓力作用下，彈性穩定時的y向位移云圖如圖7所示，開挖過程中的z向位移云圖如圖8—10所示，開挖進程中各監測點位移見表4，監測點處位移隨開挖進程的變化曲線如圖11所示。

計算結果說明：

1)隨著掌子面的掘進，軸向位移擾動區呈漏斗形，掌子面中心處位移值最大，破壞模式比較符合彈性圓板模型，以張拉破壞和剪切破壞為主。

2)掌子面中心點處的位移值最大，說明中心處最容易發生破壞，符合隔水巖體一側受溶洞水壓、另一側臨空的受力狀態。

3)側邊位移較小且變化不大，說明溶洞水壓對拱腰影響很小。

4)拱頂、拱底和中心處監測點位移均在距溶洞7 m處存在突變點。突變點左邊位移值較低，最大值在中心處(約0.003 m)，為溶洞水壓弱影響區，隔水巖體處于穩定階段；突變點右邊，位移先緩慢增加后迅速增加，為溶洞水壓強影響區，隔水巖體失穩。

5)掌子面距溶洞5 m處，掌子面中心處位移突然大幅增加至0.01 m左右，且有快速增加的趨勢，說明此時隔水巖體已開始失穩，巖體強度已不具備足夠的安全儲備。因此，可以認為該工況下的安全厚度為7 m。

圖7 彈性穩定時y向位移云圖

圖8 掌子面距溶洞17 m時z向位移云圖

圖9 掌子面距溶洞13 m時z向位移云圖

圖10 掌子面距溶洞7 m時z向位移云圖

表4 監測點位移

圖11 監測點處位移隨開挖進程的變化曲線

3.2 破裂狀態分析

由監測點位移分析結果顯示，掌子面距溶洞7 m是位移的突變點，為隔水巖體的安全厚度，下面將重點研究該部分隔水巖體的破裂狀態。隔水巖體隨掌子面推進的破裂面狀態如圖12所示。未開挖時，隔水巖體相對完整，無明顯破裂面；在距溶洞17 m時，開挖引起的應力重分布逐漸對巖體造成損傷，出現少量分散型破裂面，且隨著掌子面的推進逐漸增多；至距溶洞7 m時，隔水巖體已產生大量破裂面，呈放射狀分布。

根據上述分析結果，建立破裂度統計研究區域30 m×20 m×7 m，如圖12(a)的紅框所示。為分析開挖對研究區域的位移擾動，設立監測點于距溶洞3 m的掌子面中心。研究區域破裂度和監測點的z向位移隨開挖進程的變化曲線如圖13所示。

(a) 未開挖時

圖13 破裂度和監測點z向位移隨開挖進程變化曲線

計算結果顯示：

1)破裂度和位移曲線均有突變點，分為2階段。突變點左邊為開挖影響較弱區，二者變化均很小；突變點右邊為開挖影響強烈區，破裂度和位移均迅速增加。

2)從突變點開始，破裂度和位移均呈先緩慢增加后迅速增加，說明開挖引起的應力重分布對隔水巖體的損傷是逐漸加劇的，災變過程發生前有孕育階段。

3)破裂度突變點左邊，破裂面占比小于1%，呈零星分布，此時隔水巖體處于微損傷階段，如圖12(b)和圖12(c)所示；突變點右邊，破裂面占比快速增加，呈放射狀，掌子面處破裂面數量大于溶洞，失穩始發于掌子面。

4)破裂度突變點早于位移突變點6 m，說明隔水巖體破裂度達到一定程度后才會引起災變，位移突變是災變的表征現象。

5)距溶洞7 m處的破裂度約為8%，5 m處的破裂度約為14%，結合位移監測曲線可知，破裂度為8%時巖體損傷已處于臨界狀態。

3.3 不同水壓下隔水巖層的破裂特征分析

當溶洞水壓分別為1、2、3 MPa時，不同水壓下監測點z向位移和破裂度隨開挖進程變化曲線如圖14和圖15所示。計算結果顯示：

1)不同水壓作用下的破裂度突變點都在距離溶洞13 m處。突變點左邊破裂度基本無變化；突變點右邊，水壓越大，破裂度增大，但增幅較小，說明此種工況下水壓對隔水巖層破裂狀態的影響不大。

2)不同水壓作用下監測點位移的突變點都在距溶洞7 m處，突變點左邊位移相對穩定，水壓對遠端位移的影響很小；突變點右邊，水壓越大，位移增幅越大，因此，可認為3種水壓下的安全厚度均為7 m。

3)不同水壓下破裂度突變點均早于位移突變點6 m；在突變點右邊，二者均先緩慢增加后迅速增加。

圖14 不同水壓下監測點z向位移隨開挖進程變化曲線

圖15 不同水壓下破裂度隨開挖進程變化曲線

4 結論與建議

傳統數值方法僅能模擬巖體的單一狀態，且受軟件本身計算能力的影響，以溶洞的空間二維分布形態為主，精度較差，計算耗時多。本文采用能夠模擬隔水巖體漸進破壞的數值分析軟件GDEM-DAS，輔助以GPU(顯卡)加速，建立了百萬網格的三維隧道數值模型，通過計算界面上彈簧的強度，模擬了隔水巖體的漸進破壞過程。通過監測掌子面的位移演化規律，獲得隔水巖體的安全厚度，并利用算法的優勢提出了定量描述巖體漸進破壞過程的無量綱指標——破裂度，分析了不同水壓下隔水巖體的破裂度演化規律，以及與位移之間的對應關系。主要結論和建議如下：

1)突水為隔水巖體從裂紋萌生到貫通直至碎裂的漸進破壞過程，位移突變是破裂到災變狀態的表征現象。

2)本文給定的計算工況下，溶洞水壓為1、2、3 MPa時隔水巖體的安全厚度均為7 m，由此可見溶洞水壓非影響隔水巖體穩定性的首要因素，需結合圍巖強度、地應力特征、溶洞大小、隧道跨度和埋深等其他因素進行分析。

3)由于地質體本身的非均勻性、非連續性，現場監測所得的地表位移和勘察數據大都是表面的、局部的，很難獲取巖體內部復雜的破裂狀態和真實的材料參數，由此開展的巖體穩定性計算的可靠性大大降低。因此，可考慮將地表位移監測結果、地表破裂狀態與數值模擬結合，通過數值計算結果修正關鍵參數，獲得與現場觀測最為接近的一組數據作為材料的當前參數，并進一步計算地應力場和隔水巖體破裂狀態，得到當前破裂面積，從而開展對隔水巖體危險程度的定量化評價。

本文僅針對均質圍巖開展了隧道全斷面的開挖模擬，未考慮巖體結構的復雜性和開挖方式的多樣性，以及流固耦合作用的影響，具有一定的局限性。后續將分別針對節理裂隙隔水巖體破壞模式和溶洞水入滲誘發的隔水巖體漸進破壞效應開展進一步研究。