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

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

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

0 引言

在富水巖溶地區(qū)，隧道與溶洞間的巖層由于開挖引起的擾動(dòng)而極易失穩(wěn)，引發(fā)突水、涌泥災(zāi)害，是隧道安全施工的重大問題，也是目前工程建設(shè)研究的重點(diǎn)之一。目前研究巖溶隧道隔水巖層穩(wěn)定性的方法主要有模型試驗(yàn)、理論研究和數(shù)值模擬等。

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

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

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

1 基于GPU加速的CDEM算法

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

1.1 CDEM算法的基本概念

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

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

1.2 離散特征

采用CDEM方法描述連續(xù)-非連續(xù)問題時(shí)，連續(xù)部分采用有限元離散，非連續(xù)部分添加界面接觸表征其力學(xué)特征，相鄰接觸面通過空間拓?fù)湫畔⑦M(jìn)行查找，在公共接觸面兩側(cè)塊體的相應(yīng)節(jié)點(diǎn)間建立彈簧(包括法向彈簧和切向彈簧)，接觸面之間的相互作用力由彈簧力和彈簧的特征面積來表征。假設(shè)計(jì)算區(qū)域被結(jié)構(gòu)面T1、T2和T3切割后形成3個(gè)區(qū)域，各個(gè)區(qū)域內(nèi)需要進(jìn)行有限元離散，接觸邊界也需要進(jìn)行離散化并建立界面接觸來表征結(jié)構(gòu)面，如圖2所示。

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

1.3 破裂判斷準(zhǔn)則

CDEM采用基于時(shí)程的動(dòng)態(tài)松弛技術(shù)進(jìn)行顯式迭代計(jì)算，因此可獲知每時(shí)步的單元應(yīng)力和節(jié)點(diǎn)彈簧力，并根據(jù)破裂準(zhǔn)則進(jìn)行強(qiáng)度判斷，將發(fā)生破裂的節(jié)點(diǎn)彈簧力置為0，并記錄破裂彈簧的特征面積。本文擬在塊體部分采用線彈性模型，界面接觸彈簧采用張拉-壓剪復(fù)合準(zhǔn)則判斷其破壞狀態(tài)。張拉模式的彈脆性破壞采用最大張力準(zhǔn)則：

(1)

式中：σn為接觸法向應(yīng)力；kn為法向彈簧剛度；u為彈簧的法向位移；σt為接觸抗拉強(qiáng)度。

當(dāng)界面發(fā)生張拉破壞時(shí)，界面上的剪切應(yīng)力也相應(yīng)地置為0。剪切模式的彈脆性破裂采用摩爾-庫侖準(zhǔn)則 ：

(2)

式中：τ為接觸切向應(yīng)力；kτ為切向彈簧剛度；ν為彈簧的切向位移；c為接觸黏聚力；φ為內(nèi)摩擦角。

當(dāng)界面接觸發(fā)生破裂之后，接觸所連接的兩側(cè)塊體發(fā)生張開和滑動(dòng)。

2 計(jì)算模型及材料參數(shù)

2.1 模型尺寸

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

(a) 三維溶洞隧道網(wǎng)格

2.2 邊界條件

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

鑒于溶腔軸的尺寸相對(duì)較小，可忽略溶腔內(nèi)部上下水壓差別，水壓簡(jiǎn)化為施加于橢球內(nèi)表面的均布?jí)毫?。李利平等[6]研究發(fā)現(xiàn)瞬時(shí)突水突泥型溶洞水壓一般為1～3 MPa，本文擬分別計(jì)算水壓1、2、3 MPa時(shí)的隔水巖層破壞狀態(tài)。

2.3 材料參數(shù)

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

表1 圍巖與支護(hù)參數(shù)

2.4計(jì)算過程

計(jì)算過程為先采用彈性模型計(jì)算獲得巖體初始應(yīng)力狀態(tài)，再采用強(qiáng)度準(zhǔn)則模擬出圍巖的初始破裂狀態(tài)，然后進(jìn)行開挖過程計(jì)算： 采用全斷面開挖模型，開挖進(jìn)尺2 m，在開挖下一進(jìn)尺的同時(shí)完成上一進(jìn)尺的襯砌支護(hù)，每次開挖計(jì)算均采用破壞準(zhǔn)則迭代至40 000步。

選取掌子面為監(jiān)測(cè)面，分別在每個(gè)監(jiān)測(cè)面的拱頂、拱底、拱腰和中心位置設(shè)置監(jiān)測(cè)點(diǎn)，如圖4所示。通過監(jiān)測(cè)開挖過程中監(jiān)測(cè)點(diǎn)的位移變化規(guī)律，分析開挖擾動(dòng)對(duì)隔水巖層的影響。

(a) 監(jiān)測(cè)面 (b) 監(jiān)測(cè)點(diǎn)

2.5 隔水巖體安全厚度

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

2.6 隔水巖體的破裂度

2.6.1 破裂度的定義

破裂是表征巖體災(zāi)變過程的重要現(xiàn)象，隔水巖體的破裂程度是影響突水發(fā)生的重要因素?，F(xiàn)有的測(cè)量手段可直接測(cè)得暴露的裂縫長(zhǎng)度，但無法由此推測(cè)出巖體內(nèi)部的損傷程度以及漸進(jìn)破壞過程。

為了建立隔水巖體漸進(jìn)破裂程度與開挖進(jìn)程之間的關(guān)系，本文借助于連續(xù)-非連續(xù)數(shù)值計(jì)算方法，模型中塊體之間的界面即為潛在的破壞面，利用強(qiáng)度準(zhǔn)則判斷界面狀態(tài)，通過彈簧斷裂模擬裂縫擴(kuò)展，進(jìn)而通過統(tǒng)計(jì)彈簧力為0的特征面積獲得當(dāng)前狀態(tài)下巖體的破裂面積。將當(dāng)前狀態(tài)下的破裂面積與界面總面積之比定義為破裂度，即破裂度

D=S/St×100%。

(3)

式中：S為研究區(qū)域內(nèi)當(dāng)前斷裂彈簧特征總面積；St為研究區(qū)域內(nèi)界面總面積。

2.6.2 單軸壓縮破裂度計(jì)算和網(wǎng)格依賴性驗(yàn)證

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

表2 計(jì)算模型參數(shù)

計(jì)算模型采用張拉-壓剪復(fù)合準(zhǔn)則進(jìn)行計(jì)算，由于破裂發(fā)生在界面，通常選用隨機(jī)性較好的四面體網(wǎng)格進(jìn)行計(jì)算。不同網(wǎng)格數(shù)量的計(jì)算模型如圖5所示。為了驗(yàn)證計(jì)算結(jié)果的可靠性，對(duì)四面體網(wǎng)格進(jìn)行加密，計(jì)算模型網(wǎng)格數(shù)量見表3。

(a) 模型1

表3 計(jì)算模型網(wǎng)格、接觸面和彈簧組數(shù)量

研究區(qū)域?yàn)檎麄€(gè)圓柱，單元之間均采用界面接觸模型進(jìn)行計(jì)算，通過彈簧連接，接觸面數(shù)量和彈簧組數(shù)見表3。采用1.3節(jié)的破裂準(zhǔn)則進(jìn)行強(qiáng)度計(jì)算，通過統(tǒng)計(jì)各監(jiān)測(cè)時(shí)步彈簧力為0的特征面獲得當(dāng)前時(shí)步的破裂總面積，彈簧力不為0的特征面積為當(dāng)前時(shí)步的界面總面積，二者之比即為當(dāng)前時(shí)步的破裂度。破裂度隨時(shí)步變化曲線如圖6所示。計(jì)算結(jié)果顯示，在同一計(jì)算工況下，不同密度網(wǎng)格的破裂度隨加載時(shí)步的變化趨勢(shì)基本一致，都存在突變點(diǎn)，在突變點(diǎn)后破裂度逐漸趨于穩(wěn)定；相鄰網(wǎng)格數(shù)量的破裂度誤差小于4.3%。由此可見，破裂度對(duì)網(wǎng)格的依賴性較小。

圖6 破裂度隨時(shí)步變化曲線

3 隔水巖層破壞特征分析

3.1 隔水巖體監(jiān)測(cè)點(diǎn)位移變化分析

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

計(jì)算結(jié)果說明：

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

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

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

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

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

圖7 彈性穩(wěn)定時(shí)y向位移云圖

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

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

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

表4 監(jiān)測(cè)點(diǎn)位移

圖11 監(jiān)測(cè)點(diǎn)處位移隨開挖進(jìn)程的變化曲線

3.2 破裂狀態(tài)分析

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

根據(jù)上述分析結(jié)果，建立破裂度統(tǒng)計(jì)研究區(qū)域30 m×20 m×7 m，如圖12(a)的紅框所示。為分析開挖對(duì)研究區(qū)域的位移擾動(dòng)，設(shè)立監(jiān)測(cè)點(diǎn)于距溶洞3 m的掌子面中心。研究區(qū)域破裂度和監(jiān)測(cè)點(diǎn)的z向位移隨開挖進(jìn)程的變化曲線如圖13所示。

(a) 未開挖時(shí)

圖13 破裂度和監(jiān)測(cè)點(diǎn)z向位移隨開挖進(jìn)程變化曲線

計(jì)算結(jié)果顯示：

1)破裂度和位移曲線均有突變點(diǎn)，分為2階段。突變點(diǎn)左邊為開挖影響較弱區(qū)，二者變化均很??；突變點(diǎn)右邊為開挖影響強(qiáng)烈區(qū)，破裂度和位移均迅速增加。

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

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

4)破裂度突變點(diǎn)早于位移突變點(diǎn)6 m，說明隔水巖體破裂度達(dá)到一定程度后才會(huì)引起災(zāi)變，位移突變是災(zāi)變的表征現(xiàn)象。

5)距溶洞7 m處的破裂度約為8%，5 m處的破裂度約為14%，結(jié)合位移監(jiān)測(cè)曲線可知，破裂度為8%時(shí)巖體損傷已處于臨界狀態(tài)。

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

當(dāng)溶洞水壓分別為1、2、3 MPa時(shí)，不同水壓下監(jiān)測(cè)點(diǎn)z向位移和破裂度隨開挖進(jìn)程變化曲線如圖14和圖15所示。計(jì)算結(jié)果顯示：

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

2)不同水壓作用下監(jiān)測(cè)點(diǎn)位移的突變點(diǎn)都在距溶洞7 m處，突變點(diǎn)左邊位移相對(duì)穩(wěn)定，水壓對(duì)遠(yuǎn)端位移的影響很??；突變點(diǎn)右邊，水壓越大，位移增幅越大，因此，可認(rèn)為3種水壓下的安全厚度均為7 m。

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

圖14 不同水壓下監(jiān)測(cè)點(diǎn)z向位移隨開挖進(jìn)程變化曲線

圖15 不同水壓下破裂度隨開挖進(jìn)程變化曲線

4 結(jié)論與建議

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

1)突水為隔水巖體從裂紋萌生到貫通直至碎裂的漸進(jìn)破壞過程，位移突變是破裂到災(zāi)變狀態(tài)的表征現(xiàn)象。

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

3)由于地質(zhì)體本身的非均勻性、非連續(xù)性，現(xiàn)場(chǎng)監(jiān)測(cè)所得的地表位移和勘察數(shù)據(jù)大都是表面的、局部的，很難獲取巖體內(nèi)部復(fù)雜的破裂狀態(tài)和真實(shí)的材料參數(shù)，由此開展的巖體穩(wěn)定性計(jì)算的可靠性大大降低。因此，可考慮將地表位移監(jiān)測(cè)結(jié)果、地表破裂狀態(tài)與數(shù)值模擬結(jié)合，通過數(shù)值計(jì)算結(jié)果修正關(guān)鍵參數(shù)，獲得與現(xiàn)場(chǎng)觀測(cè)最為接近的一組數(shù)據(jù)作為材料的當(dāng)前參數(shù)，并進(jìn)一步計(jì)算地應(yīng)力場(chǎng)和隔水巖體破裂狀態(tài)，得到當(dāng)前破裂面積，從而開展對(duì)隔水巖體危險(xiǎn)程度的定量化評(píng)價(jià)。

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