深埋隧道前方承壓溶洞隔水巖體最小安全厚度研究

房忠棟，楊為民，王旌，石錦江，巴興之，王浩

(1.山東大學巖土工程中心，山東濟南，250061；2.山東大學齊魯交通學院，山東濟南，250002；3.濟南軌道交通集團有限公司，山東濟南，250101；4.山東高速集團有限公司，山東濟南，250098)

我國山區(qū)面積占國土總面積的2/3，并且是世界上巖溶分布最廣的國家，西南地區(qū)最發(fā)育。近年來，交通事業(yè)的快速發(fā)展促進了隧道工程的蓬勃發(fā)展。在巖溶山區(qū)將不可避免地修建大量深埋隧道，巖溶對于隧道的危害主要是誘發(fā)突水涌泥，突水涌泥誘發(fā)地面次生災害，如造成地面塌陷、井泉干涸等。巖溶突水已經(jīng)成為隧道、礦山巷道等地下工程建設中常見的地質(zhì)災害之一。深埋隧道導致圍巖自重應力處于較高水平，圍巖所處環(huán)境的流?固強耦合效應不可忽略。一般認為，巖溶區(qū)隧道與巖溶構(gòu)造間巖層防突安全厚度不夠，導致隔水巖體在開挖擾動作用下垮塌或者在高水壓作用下被擊穿是引發(fā)巖溶隧道突水涌泥等工程事故的主要因素[1?5]。因此，必須對深埋巖溶隧道最小防突巖層安全厚度進行研究。

巖溶隧道施工中經(jīng)常遇到不同發(fā)育程度的溶洞，隱伏溶洞的隱蔽性導致災害往往具有突發(fā)性，嚴重威脅工程和人員安全[6]。在隱伏溶洞安全性研究方面，國內(nèi)外學者針對溶洞的尺寸、分布位置、充填物性質(zhì)、水壓等因素對隧道穩(wěn)定性的影響進行了研究，取得了一些實質(zhì)成果。賴金星等[7]研究了不同位置的溶洞對隧道結(jié)構(gòu)的影響，發(fā)現(xiàn)溶洞位置對于收斂值的影響由大至小依次為隧底溶洞、隧頂溶洞和側(cè)壁溶洞，對于沉降值的影響由大至小依次為側(cè)壁溶洞、隧頂溶洞和隧底溶洞，該研究未考慮溶洞含水狀態(tài)下對隧道穩(wěn)定性的影響。一些研究者通過構(gòu)建耦合計算模型，借助數(shù)值模擬和模型試驗等方法，針對富水隧道圍巖穩(wěn)定性與防突巖層安全厚度開展了大量研究，如：李術(shù)才等[8]開展了多種工況下隱伏溶洞突水災變演化過程模擬試驗，將承壓隱伏溶洞突水過程劃分為群裂隙萌生擴展、優(yōu)勢導水通道形成、隔水巖體破裂失穩(wěn)3 個階段；潘東東等[9]綜合模型試驗和數(shù)值模擬手段，針對強充填滯后型溶洞突水模式，基于隧道圍巖穩(wěn)定性受不同充水條件隱伏溶洞的影響規(guī)律，揭示了溶洞突水致災機理；張橋[10]將圍巖的位移變化情況和塑性區(qū)分布情況作為評價隧道圍巖突水的判據(jù)，得出了隧道圍巖防突巖層的最小安全厚度；李濤等[11]通過分析掌子面位移演化規(guī)律獲得不同水壓下隔水巖體的最小安全厚度，認為位移突變是破裂到災變狀態(tài)的表征現(xiàn)象；SHAN等[12]分析了成—桂鐵路上隧道與周圍隱伏洞之間的安全厚度和巖石破碎機理，提出了基于玉景山隧道工程的綜合數(shù)值分析方法；王志杰等[13]通過研究隧道前方存在正交高壓溶洞時中間巖墻的承壓能力和破壞模式，建立了掌子面失穩(wěn)破壞的圓錐臺模型，并通過勢能判據(jù)的尖點突變理論得到掌子面失穩(wěn)時的溶洞臨界壓力。

現(xiàn)有研究主要集中在隧道洞身范圍外存在隱伏溶洞時圍巖結(jié)構(gòu)的穩(wěn)定性、防突巖層安全厚度及突水機理等方面。但當溶洞位于隧道前方時，隨著掌子面的推進，溶洞必然會被揭露，更易導致突水災害的發(fā)生。深埋巖溶隧道面臨高應力和強水壓的耦合作用，對隔水巖體最小安全厚度研究尚不夠深入，考慮隧道開挖擾動與流?固耦合綜合作用的分析較少。針對上述問題，本文以重慶歇馬隧道為工程背景，利用COMSOL Multiphysics軟件對流?固耦合計算和隧道開挖過程進行模擬分析，確定深埋隧道前方隱伏承壓溶洞隔水巖體最小安全厚度，并分析溶洞尺寸、水壓以及隧道埋深、跨度等因素對最小安全厚度的影響規(guī)律，最后通過多元線性回歸方法推導出巖溶隧道前方承壓水溶洞掌子面預留安全厚度公式，以期對類似工程施工安全提供技術(shù)支撐。

1 工程背景

歇馬隧道是重慶快速路工程中“一橫線”工程的控制性工程，設計為分離式雙洞雙向六車道行駛隧道，隧道開挖凈寬為14.5 m，凈高為9.6 m。隧道全長4 200 m，最大隧道埋深為392 m，屬于深埋特長隧道。隧道穿越以觀音峽背斜軸部隆起為主的中梁山，山頂部受碳酸鹽巖溶蝕作用的影響，東西方向形成了劉家槽和窯灣共2 條巖溶槽谷，呈“一山兩槽三嶺”形態(tài)。在2 條巖溶槽谷段，地表及地下水豐富，隧道存在巖溶突水涌水等地質(zhì)災害，地質(zhì)剖面圖如圖1所示。

圖1 歇馬隧道地質(zhì)剖面圖Fig.1 Geological profile of Xiema Tunnel

當隧道右線開挖至里程樁號YK9+268 時，圍巖體較破碎，節(jié)理裂隙較發(fā)育，經(jīng)超前鉆探得掌子面前方有1 個體積約(長×寬×高)為5 m×10 m×15 m的充水溶洞。根據(jù)地質(zhì)勘探報告顯示，YK9+266—YK9+296 段隧道平均埋深235 m，飽和抗壓強度為34.36 MPa，完整性系數(shù)為0.68，綜合判定為Ⅲ級圍巖。靜水壓力為1.8～2.0 MPa，地下水活動狀態(tài)多為大股涌水或噴射狀涌水。

2 突水計算模型的構(gòu)建

深埋隧道前方存在承壓水溶洞時，隨著開挖擾動作用的影響，圍巖內(nèi)部滲流場和應力場變化明顯，且相互作用，表現(xiàn)為2 個物理場的強耦合[14]：一方面，高水壓作用引起溶洞附近圍巖的應力場變化；另一方面，隧道的開挖導致圍巖的應力釋放，圍巖孔隙率和滲透率發(fā)生改變，從而改變滲流場的分布。在開挖擾動和滲流作用的雙重作用下，隧道圍巖強度逐漸降低，當開挖擾動損傷帶H1和水壓致裂帶H2連通時，即可認為隧道發(fā)生突水，此時的巖層厚度為最小安全厚度HS，最小安全厚度示意圖如圖2所示。

圖2 最小安全厚度示意圖Fig.2 Schematic diagram of minimum safe thickness

2.1 模型構(gòu)建與參數(shù)選取

本文使用COMSOL Multiphysics軟件中的固體力學模塊和達西定律模塊對流?固耦合作用時隧道開挖過程進行模擬計算分析。引入以下基本假設：1)圍巖材料服從Mohr?Coulomb屈服準則；2)地下水滲流發(fā)生時服從Darcy定律；3)溶洞經(jīng)過長期溶蝕作用已達到穩(wěn)定狀態(tài)，溶腔內(nèi)充滿水體。

以歇馬隧道突水段為工程背景，建立數(shù)值計算模型。為了消除邊界效應，整個模型長為120 m，寬為60 m，高為80 m；隧道開挖跨度為14.5 m，高度為9.6 m，高跨比為0.66。模型采用自由四面體進行網(wǎng)格剖分，總單元數(shù)10 218 個，具體計算模型如圖3所示。根據(jù)地質(zhì)勘察資料和巖石物理力學參數(shù)測定試驗，確定計算模型中圍巖物理力學參數(shù)如表1所示。

表1 圍巖的物理力學參數(shù)Table 1 Physical and mechanical parameters of surrounding rock

圖3 計算模型Fig.3 Calculation model

應力場邊界條件為：在巖體應力場變形計算中，模型的四周邊界設為輥支撐約束，底部邊界為固定約束，模型頂部和隧道開挖臨空面均設定為自由面；模型內(nèi)不考慮構(gòu)造應力，初始地應力按照自重應力計算，側(cè)壓力系數(shù)取1；隧道埋深為235 m，構(gòu)建的模型上表面至隧道拱頂為40 m，剩余195 m 巖體自重作為外部荷載作用于模型上表面，計算后為5 MPa。

滲流場邊界條件為：在多孔介質(zhì)滲流場計算中，模型的外邊界均設定為不透水邊界，質(zhì)量通量為0 kg/(m2·s)，溶洞內(nèi)部水壓設置為2 MPa，隧道開挖輪廓邊界設定為透水邊界，壓力水頭為0 m。

2.2 開挖過程實現(xiàn)

將隧道輪廓內(nèi)開挖巖體沿隧道軸線分割為多個柱體，設置多個計算步驟，每個計算步驟包含一組物理場(固體力學和達西定律)。每步計算開始時選擇計算范圍并設定參數(shù)。將待開挖的計算域鈍化，使其不參與計算，以實現(xiàn)隧道的開挖過程模擬。計算域分布見圖4。具體操作如下：

圖4 計算域分布圖Fig.4 Computational domain distribution map

1)求解固體力學和達西定律物理場信息，計算范圍為整個模型體，模擬隧道開挖前的平衡狀態(tài)，即計算模型體的初始應力狀態(tài)。

2)以步驟1)的計算結(jié)果作為該步驟的初始值，計算范圍為去掉計算域2的模型體，即計算域2不參與計算，模擬第一步開挖。

3)以步驟2)的計算結(jié)果作為該步驟的初始值，計算范圍為去掉計算域2和計算域3的模型體，即計算域2和計算域3 不參與計算，模擬第二步開挖。

4)依此類推，直至隧道開挖范圍內(nèi)的所有計算域全部刪除，完成隧道的開挖過程模擬。

3 突水計算結(jié)果分析

通過COMSOL Multiphysics計算模擬歇馬隧道突水段開挖過程，分析隧道開挖完成后圍巖位移變形分布規(guī)律，同時進一步分析考慮耦合作用下開挖后圍巖滲透壓力、滲流速度和塑性區(qū)的分布特征，總結(jié)隔水巖層最小安全厚度確定方法，并通過模型試驗結(jié)果對模擬結(jié)果的有效性進行驗證。

3.1 數(shù)值計算結(jié)果

綜合分析隧道開挖過程中圍巖滲流場與應力場的變化規(guī)律，結(jié)合開挖擾動和水壓作用下主應力差和塑性區(qū)分布特征，建立隔水巖層最小安全厚度確定方法。

3.1.1 滲流場變化規(guī)律

取隧道軸線所在豎直截面作為研究對象，圍巖滲流場分布如圖5所示(圖中箭頭的大小代表滲流速度的大小)。由圖5(a)可見：隧道開挖后導致掌子面臨空，此時隧道輪廓為透水邊界，故在掌子面附近出現(xiàn)較大的滲流速度集中區(qū)，但由于掌子面距離溶洞較遠，掌子面與溶洞間的巖體內(nèi)滲流速度很小，故在掌子面部位不會發(fā)生突水災害。由圖5(b)～(f)可見：隨著隧道逐步開挖，隔水巖體厚度逐漸降低，溶洞水流動所受阻力減小，圍巖內(nèi)的滲流速度逐漸增加，溶洞內(nèi)承壓水引起的滲流場逐漸擴展至掌子面，同時，溶洞除了靠近掌子面一側(cè)的邊界外，其他方向的滲流速度均逐漸減少，說明溶洞內(nèi)承壓水逐漸通過掌子面前方巖體向隧道內(nèi)流入。

圖5 開挖過程中圍巖滲流場分布圖Fig.5 Distribution of seepage field in surrounding rock during excavation

3.1.2 塑性區(qū)分布規(guī)律

溶洞水壓和隧道開挖擾動均導致溶洞周圍或隧道邊界附近形成塑性區(qū)，圖6所示為隧道軸線截面圍巖塑性區(qū)云圖。由圖6(a)～(c)可見：隨著掌子面逐漸靠近溶洞，原本開挖擾動控制的塑性區(qū)逐漸受溶洞的影響，該部分塑性區(qū)有增大的趨勢；溶洞承壓水控制的塑性區(qū)首先在溶洞邊界附近形成，在隧道的開挖過程中逐漸向掌子面方向發(fā)展。由圖6(d)可見：當隧道開挖至溶洞前方4 m 位置時，掌子面附近塑性區(qū)和溶洞附近塑性區(qū)發(fā)生貫通，認為此時突水通道已經(jīng)形成，即該施工工況對應的隔水巖體最小安全厚度在4～5 m之間。

圖6 開挖過程中圍巖塑性區(qū)分布圖Fig.6 Distribution of plastic zone in surrounding rock during excavation

3.1.3 隔水巖體最小安全厚度確定

圍巖主應力差與圍巖破壞的發(fā)生密切相關(guān)，相比于最大(小)主應力更能反映圍巖狀態(tài)[15]。通過模擬計算得到的隧道開挖過程中掌子面主應力差(最大主應力和最小主應力之差)與掌子面最大水平位移，其量綱和取值范圍存在較大差異。為了更方便分析這2個物理量變化產(chǎn)生的影響，將它們歸一化處理，并繪制于同一坐標系中，如圖7所示。由圖7可知：當掌子面與溶洞距離d>12 m 時，掌子面上主應力差和最大水平位移增長緩慢，即當掌子面距離溶洞較遠時，掌子面受溶洞影響較小，圍巖內(nèi)的流?固耦合作用不明顯，隧道開挖過程的擾動作用是引起應力和位移變化的主要因素；當5

圖7 開挖過程中掌子面信息曲線圖Fig.7 Curves of tunnel face information during excavation

模擬過程中靠近溶洞時的每循環(huán)開挖進尺為1 m，綜合以上模擬結(jié)果，認為隔水巖體最小安全厚度在4～5 m 之間。因此，改變最后一步開挖進尺，以0.1 m 為最小分度值，設計多組計算工況，重新模擬以上過程，以開挖擾動引起的塑性區(qū)和溶洞周圍塑性區(qū)貫通時掌子面與溶洞邊界距離作為隔水巖體的最小安全厚度，最終確定該工況最小安全厚度為4.7 m，此時的塑性區(qū)分布如圖8所示。

圖8 塑性區(qū)分布圖Fig.8 Distribution of plastic zone

3.2 溶洞突水模型試驗對比

文獻[16]依托歇馬隧道工程突水實例，采用1:50的模型相似比尺，開展了隧道前方充水溶洞突水災變演化模型試驗，總結(jié)了溶洞突水的3個階段即導水裂隙萌生、導水裂隙擴展和導水裂隙貫通，隧道開挖完成后，掌子面先后經(jīng)歷了滴狀滲水—股狀流水—突水等過程，見圖9。

圖9 試驗突水過程圖Fig.9 Diagrams of water inrush process during test

綜合試驗模型體位移場、應力場和滲流場信息變化規(guī)律，最終確定該工況下隔水巖體最小安全厚度為4.5 m。隔水巖體最小安全厚度數(shù)值模擬結(jié)果與該試驗結(jié)果相對誤差為4.4%，因此，本文所選數(shù)值模擬方法與計算模型可用于類似工程隔水巖體最小安全厚度的研究。

4 不同因素對隔水巖體最小安全厚度的影響

隔水巖體最小安全厚度的影響因素較多，并且需考慮具體的施工工況。為研究不同因素的影響效果，本文采用正交試驗進行研究。正交試驗設計是根據(jù)正交性從全面試驗中挑選部分具有代表性的點進行研究，是研究多因素多水平的試驗設計方法[17]，具有均勻分散、齊整可比的優(yōu)點，即通過正交試驗設計的各因素水平組合在全部組合中的分布是均勻的，并且每個因素的各水平之間具有可比性。

4.1 正交模擬試驗設計

基于巖溶隧道隔水巖層最小安全厚度的研究成果[18?20]，研究溶洞水壓P、溶洞直徑D、隧道跨度L和隧道埋深h這4 個因素對隔水巖層最小安全厚度的影響規(guī)律，各因素分別選擇4個水平，采用如表2所示的L16(44)正交試驗表設計模擬工況，具體參數(shù)設置見表3。

表2 L16(44)正交試驗表Table 2 L16(44)orthogonal test table

表3 因素水平設置Table 3 Factor level setting

4.2 模擬結(jié)果分析

4.2.1 不同工況最小安全厚度計算結(jié)果

基于4.1節(jié)中正交試驗表(表2)和各因素水平設置(表3)規(guī)則，確定各模擬工況的參數(shù)取值，開展16 組模擬試驗。根據(jù)隔水巖層最小安全厚度確定方法，確定各模擬工況隔水巖體最小安全厚度，如表4所示。

表4 最小安全厚度模擬結(jié)果Table 4 Simulation results of minimum safe thickness

4.2.2 各因素對最小安全厚度的影響規(guī)律

各個因素對最小安全厚度的影響程度各異，通過對上述模擬結(jié)果進行極差分析，確定其主次順序和重要程度，并分析各因素單獨作用效應。

1)極差分析。最小安全厚度的模擬結(jié)果極差分析如表5所示，其中，Ki為某個影響因素的第i個水平所對應模擬結(jié)果的總和，ki為某個影響因素的第i個水平模擬結(jié)果的平均值，Rj為第j列因素的各個水平所對應模擬結(jié)果的平均值的極差。Rj反映的是第j列因素的水平變動引起的觀測指標的變化幅度，Rj越大，則說明該因素對試驗指標的影響越大。通過對深埋隧道前方隱伏承壓溶洞的隔水巖體最小安全厚度模擬結(jié)果的極差分析，可得出如

表5 模擬結(jié)果極差分析Table 5 Range analysis of simulation results

下關(guān)系：R溶洞水壓>R隧道埋深>R溶洞直徑>R隧道跨度，即認為各因素中，對隔水巖體最小安全厚度影響最大的是溶洞水壓，其次是隧道埋深，影響最小的是隧道跨度。

2)各因素影響趨勢。根據(jù)隔水巖體最小安全厚度極差分析結(jié)果分別繪制各因素對最小安全厚度的影響規(guī)律曲線，如圖10所示。由圖10(a)可見：當其他因素確定時，最小安全厚度隨溶洞水壓增大而增大，并且水壓越大，最小安全厚度增大幅度越大。這是因為隧道突水通道是溶洞周圍巖體受水壓作用產(chǎn)生的破碎區(qū)和隧道開挖擾動產(chǎn)生的破碎區(qū)貫通而形成的，溶洞內(nèi)的水壓越大，則水壓致裂產(chǎn)生的破壞區(qū)域越大，從而需要更大的隔水巖體厚度來抵抗洞內(nèi)水壓的作用。

由圖10(b)可見：當其他因素確定時，最小安全厚度隨溶洞直徑增大而增大；隧道開挖前的地層各物理場均處于平衡狀態(tài)，每一步隧道開挖步都對隧道周圍巖體造成一定的擾動，圍巖通過一定變形形成新的平衡狀態(tài)，而溶洞的存在放大了該效應，溶洞越大，隧道開挖引起的溶洞周邊物理場變化范圍越大，對隔水巖體的穩(wěn)定性越不利，從而導致最小安全厚度越大。

圖10 各因素對最小安全厚度的影響Fig.10 Influence of various factors on the minimum safe thickness

由圖10(c)可見：當其他因素確定時，在隧道高跨比不變時，隨著隧道跨度增大，最小安全厚度也有增大趨勢，但相較于溶洞水壓和溶洞直徑而言，隧道跨度對最小安全厚度的影響效應并不突出。在通常情況下，由于巖體結(jié)構(gòu)面的存在，圍巖壓力隨著隧道跨度增大而增大。本文采用的有限元模擬方法將巖體假定為連續(xù)體，未考慮結(jié)構(gòu)面的影響，因此，最小安全厚度并未隨隧道跨度增大而大幅度增加。

由圖10(d)可見：當其他因素確定時，隔水巖體最小安全厚度隨著隧道埋深增大而減小。這是因為隧道埋深越大，開挖前隧道原平衡狀態(tài)中圍巖的地應力越大，對溶洞水壓致裂效應起到抑制作用。因此，模擬結(jié)果顯示，當溶洞內(nèi)水壓和溶洞直徑一定時，隧道埋深越大，溶洞周圍圍巖因水壓致裂導致的破壞區(qū)越小，從而導致最小安全厚度越小。

4.3 隔水巖體最小安全厚度預測模型

多元線性回歸分析方法是指對2個或2個以上的事件因素設置成自變量，而將事件結(jié)果設置為因變量，因變量與每一個自變量呈線性變化關(guān)系，通過相關(guān)性分析，建立因變量與若干自變量的線性關(guān)系式，由自變量變化進行回歸分析。通過采用多元回歸分析方法對正交數(shù)值模擬結(jié)果進行分析。

4.3.1 一元回歸分析

以影響因素溶洞水壓P、溶洞直徑D、隧道跨度L、隧道埋深h分別作為自變量，最小安全厚度HS作為因變量，選擇基本的回歸方程(線性函數(shù)、指數(shù)函數(shù)、對數(shù)函數(shù)、冪函數(shù)等)分別進行相關(guān)性分析，并取擬合優(yōu)度最大的函數(shù)作為其對應的回歸方程。

1)溶洞水壓P與最小安全厚度HS的關(guān)系：

擬合優(yōu)度R2=0.916 4。

2)溶洞直徑D與最小安全厚度HS的關(guān)系：

擬合優(yōu)度R2=0.930 6。

3)隧道跨度L與最小安全厚度HS的關(guān)系：

擬合優(yōu)度R2=0.934 9。

4)隧道埋深h與最小安全厚度HS的關(guān)系：

擬合優(yōu)度R2=0.981 6。

以上公式的擬合優(yōu)度均大于0.900 0，擬合效果較好，可認為所得經(jīng)驗公式有意義。

4.3.2 多元回歸分析

根據(jù)上述的單因素回歸關(guān)系，建立巖溶隧道掌子面前方隔水巖層最小安全厚度的預測模型。首先假設最小安全厚度HS與e0.0937P，ln(L)，D和h之間存在如下線性關(guān)系：

式中：ai(i=1，2，…，5)為待定系數(shù)。

通過多元線性回歸分析方法求出各待定系數(shù)，得最小安全厚度預測公式：

式中：1.5≤P≤6.0 MPa；9≤L≤18 m；4≤D≤13 m；300≤h≤900 m。

對式(6)進行顯著性檢驗，驗證回歸方程的線性相關(guān)性。式(6)對應的方差分析結(jié)果表明：F服從自由度為(4,11)的分布，由F分布檢驗表查得F0.05(4,11)=3.357，F(xiàn)0.01(4,11)=5.668；F>F0.01(4,11)，并且P<0.001，說明在顯著性水平α=0.01 下，HS與e0.0937P，lnL，D和h之間存在顯著的線性關(guān)系，即回歸方程顯著，F(xiàn)檢驗通過。

5 結(jié)論

1)當深埋隧道前方承壓溶洞發(fā)生突水時，掌子面圍巖主應力差和最大水平位移均發(fā)生突降，掌子面和溶洞間圍巖塑性區(qū)貫通，可以此確定隔水巖體最小安全厚度取值。

2)本文采用數(shù)值計算方法得到的隔水巖體最小安全厚度為4.7 m，與模型試驗結(jié)果相對誤差為4.4%。

3)對隔水巖體最小安全厚度HS影響重要性程度由大至小依次為溶洞水壓P、隧道埋深h、溶洞直徑D和隧道跨度L，并且最小安全厚度隨著溶洞水壓、溶洞直徑和隧道跨度增大而增大，隨著隧道埋深增大而減小。

4)最小安全厚度HS與e0.0937P，lnL，D和h之間存在極其顯著的線性關(guān)系。