基于流-固-損傷全耦合模型的裂隙型巖溶突涌水數值模擬

劉家慶, 林志, 邸小勇, 鞏雯, 劉先林, 邵羽

(1.廣西新發展交通集團有限公司, 南寧 530029; 2.重慶交通大學省部共建山區橋梁及隧道工程國家重點實驗室, 重慶 400074;3.重慶高速工程顧問有限公司, 重慶 404100; 4.廣西交通設計集團有限公司, 南寧 530029)

中國巖溶地貌面積約2.55×104km2,大多分布在云南、貴州等西南地區,是世界上巖溶地貌面積最大的國家[1]。隨著中國建設重心逐漸西遷,隧道工程的建設數量日益增多。大量工程實例表明,不良水文地址條件下,當施工揭露或即將揭露巖溶含水構造時,對局部隔水巖體產生擾動,易導致突涌水事故的發生。如何有效對巖溶地區隧道工程建設常見突涌水災害進行預防和控制成為熱點[2],巖溶突涌水機理與臨界條件的研究勢在必行。

大量學者對巖溶突涌水致災機理及臨界條件展開研究。馬國民等[3]采用非連續變形分析(discontinuous deformation analysis, DDA)方法對高壓巖溶水作用下突涌水事故過程進行分析,提出溶腔水頭與隧道突涌水事故的發生高度相關。武世燕[4]將不同模擬工況下巖溶隧道最小隔水巖盤厚度計算結果進行回歸分析,提出隔水巖盤最優設計厚度理論公式。李利平等[5]結合巖溶地質學、工程水力學和斷裂力學等相關理論,推導巖溶隧道巖溶裂隙水突出的最小巖體防突厚度半解析表達式。曹林衛等[6]利用巖石破裂過程分析程序(real failure process analysis, RFPA),對隧道施工誘發復合圍巖破裂突涌水展開研究,揭示了復合圍巖巖盤厚度與抗水壓能力的關系。Liang等[7]通過物理模擬實驗,記錄隧道開挖全過程圍巖應力、位移、水壓等數據,分析了應力、位移、水壓的變化規律,并結合工程案例建立突水災害是否會發生的判斷標準。Wang等[8]將滲流速率、圍巖應力值、圍巖位移量作為控制變量,分別在拱腰、拱肩布置各物理量監測傳感器,形成一種基于模糊理論的多信息區間融合預警新方法,并在實際工程中充分驗證其有效性與實用性。現有研究成果普遍認為,隧址區特大體量富水溶腔附近圍巖在溶腔高壓裂隙水作用下,發生裂隙發育、貫通并在一定范圍內形成破壞區[9],巖體滲透系數增大,當富水溶腔與開挖面間巖體厚度儲備不足時,極易發生突涌水災害。在隧道開挖擾動下,富水溶腔高壓裂隙水將激活并推動巖體裂隙、節理進一步發育,導致開挖階段對安全厚度判斷失誤,引發突涌水災害,但少有基于此展開的研究。

綜上,已有研究多將巖體視為均一介質,通過對巖體強度折減來模擬天然裂隙的存在,鮮有同時考慮裂隙存在,且裂隙隨水壓發育而誘發的裂隙型突涌水機理與臨界條件研究成果。現基于“流-固-損傷”多場耦合數值模型,對裂隙型巖溶突涌水機理與臨界條件展開研究,希望推動隧道品質工程的建設。

1 “流-固-損傷”全耦合模型

隧址區有巖溶儲水構造,隧道的開挖對圍巖的擾動導致巖體裂隙進一步發育,儲層數量增加、防突巖層整體破壞。圖1為裂隙型巖溶示意圖。

圖1 裂隙型巖溶示意圖Fig.1 Diagram of fissure Karst

多孔巖體的機械平衡的控制公式[10]為

?σ+f=0

(1)

式(1)中:σ為應力;f為體力。

巖石材料的損傷模型由以下模型控制,即

σ=(1-ω)De:ε

(2)

式(2)中:De為彈性剛度張量;ε為應變張量;ω為損傷變量,數值從0(未損壞巖石)至1(損壞巖石);“:”表示張量內積或外積。

1.1 達西定律物理場

巖體介質中的孔隙一般較小,水在巖體中流動時受到黏滯阻力的影響,流速較慢,因此其流動狀態多為層流。破碎多孔介質中單向流的連續性方程為

(3)

(4)

式中:ρw為水的密度;n為孔隙度(假定天然裂隙中無填充,裂隙孔隙度設置為1);t為時間;U為速度;Qm為源項;α為biot系數;εv為固體骨架的體積應變;k為滲透率;μw為水的動力黏度;P為流體壓力。

裂隙中滲透率的計算[11]公式為

(5)

式(5)中:bn為法相孔徑。

巖體的滲透率計算公式[12]為

(6)

式(6)中:n0與k0分別為初始孔隙度和初始滲透率;αk為損傷-滲透率相關系數,取5.0;n為巖體孔隙度。

裂隙及巖體的儲水模型為

(7)

式(7)中:S為儲能系數。

巖體的儲能系數為

Sm=ncw+(1-α)(α-n)cm

(8)

式(8)中:cw為水的壓縮系數,取4.74×10-10Pa-1;cm為巖石的壓縮系數。

1.2 固體力學物理場

應力場作用下,固體骨架孔隙或閉合,或變小,從而導致滲透系數的改變。裂隙在高應力狀態下也會趨于閉合。因此,應力場主要對固體骨架的孔隙度以及裂隙孔徑的改變其關鍵作用。

在法相壓縮荷載下,預存裂隙的孔徑遵循指數函數[12]為

bn=br+(b0-br)exp(-ξσ′n)

(9)

式(9)中:bn為法相裂隙孔徑;b0為初始裂隙孔徑;br為殘余孔徑;σ′n為有效法相壓應力;ξ為孔隙壓力系數。

σ′n=σN-p

(10)

式(10)中:σN為總應力;p為孔隙壓力。

巖石斷裂的剪切行為遵循Mohr Coulom理論,即

(11)

式(11)中:τs為剪應力;us為剪切位移;Ks為剪切剛度;τp為峰值剪應力;up為開始斷裂滑動的峰值剪切位移。

up=τp/Ks

(12)

τp=c′+σ′tanφ′

(13)

式中:c′、φ′為有效應力強度指標;σ′為作用在剪切面上的有效應力。

應力影響下,巖體的孔隙度受以下方程控制[12],即

(14)

1.3 “流-固-損傷”全耦合模型機理

通過梳理應力場與滲流場控制方程,發現實現流-固-損傷全耦合的關鍵在偏微分方程的表達與各參數之間的偏微分關系。

在固體力學物理場中,巖體介質的孔隙度受應力影響。根據上述微分方程,巖體的孔隙度與裂隙的張開量受應力影響,應力場的分布受滲流場的影響。隨著應力數值的變化,固體相也會相應產生損傷(圖2)。將上述控制方程應用在數值模型中,即可得到流-固-損傷全耦合模型[10]。

圖2 流-固-損傷全耦合作用機理Fig.2 Mechanism of ‘fluid-solid-damage’ full coupling

2 裂隙型巖溶突水數值分析

2.1 工程概況

以安康至來鳳國家高速公路奉節至巫山(渝鄂界)段小營特長隧道為工程背景,隧道進洞口位于重慶市巫山縣廟宇鎮附近,出洞口位于湖北省建始縣沙壩村附近,隧道線位整體呈北西-南東向走勢。隧道左線起訖樁號為K37+855～K50+241,長12 386 m;右線起訖樁號為YK37+857～YK50+158,長12 304 m,為奉建高速在建最長隧道。

初勘期間,發現隧址區山體中部區域出現大規模可溶巖(圖3)。為此,展開巖溶地質水文專項調查。根據“小營隧道”專水調查報告,擬建隧道嘉陵江組一、三段以質純的灰巖為主,特別是三段為中厚層狀,直接受大氣降雨補給,各種巖溶現象較為發育。調查區可溶巖出露范圍大,在灰巖與頁巖交界處的灰巖一側出現巖溶水聚集現象,形成巖溶強烈發育且十分富水地帶。區內巖層傾角較陡,層面裂隙溶隙十分發育,易于地下水的下滲補給。地質構造格局表明隧道施工過程中極易發生裂隙型巖溶突涌水,受大氣降雨直接補給下,巖溶水壓基本保持不變。分別根據大氣降雨入滲系數法[13]和徑流模數法[14]對隧道的突涌水量進行了預測。兩種計算方法的計算結果差異不大,左洞為47 627 m3/d,右洞為47 187 m3/d,雙洞為94 814 m3/d;豐水期涌水量可溶巖乘以系數3,非可溶巖乘以系數2,得左洞為140 523 m3/d,右洞為139 213 m3/d,雙洞為279 736 m3/d。

圖3 不良地質段縱剖面圖與代表性溶洞照片Fig.3 Longitudinal section of adverse geological conditions and photos of representative Karst cave

為分析富水區段巖溶裂隙型突涌水對隧道施工的影響,揭示其致災臨界條件,基于Comsol Multiphysics強大的偏微分方程求解能力,開展基于流-固-損傷耦合模型的裂隙型巖溶突涌水機理與臨界條件的研究。

2.2 模型建立

為體現流-固-損傷全耦合理論下裂隙型巖溶突涌水發育過程,利用comsol livelink with MATLAB接口在120 m×100 m的矩形地層區域隨機生成若干傾角分別為25°和65°的裂隙,在地層結構中布置隧道斷面模型與溶洞模型并與隧道斷面四周布置控制測點。建立幾何模型如圖4所示。為準確模擬裂隙作用,利用Comsol中裂隙流和彈性薄層對裂隙進行網格建模,如圖5所示。

圖4 裂隙型巖溶突水計算模型Fig.4 Calculation model of fissure Karst water inrush

圖5 裂隙網格單元Fig.5 Grid cell of fracture

在此模型中,將裂隙視為無質量、無體積,但其邊界存在對流體流速影響的網格模型。為此,在固體力學場中,選中所有裂隙,將其單位面積的質量設置為0,將其視為無質量的“空腔”。達西定律場中,需對裂隙網格設置孔徑以及儲水模型參數,將其孔徑設置為初始裂隙開度0.1 mm,且模型孔徑隨物理場更新而更新;裂隙網格的孔隙率設置為1,以模擬其“無體積”,滲透率設置為與裂隙孔徑相關,如1.1節中巖石滲透率公式一致;裂隙邊界對流體的摩擦力通過設置動力黏度系數來實現。

固體物理場中,地層上部設置法相荷載5 MPa,以模擬上覆200 m填土。

物理場中材料參數按照第1節推導公式進行設置。材料物理力學基本參數如表1所示。

表1 材料物理力學參數

裂隙型巖溶突涌水的本質是高水壓作用下,巖體裂隙發生擴展并相互貫通再進一步張開發生劈裂[15-17]。為便于計算,模擬過程中,設定溶洞水壓力恒定,隧道斷面為出水口并設置為零壓力。巖體為均勻各向同性,孔隙度φ=0.3,圍巖滲透率km=1×10-11m2,裂隙開度為0.1 mm,流體密度ρ=1 000 kg/m3,流體動力黏度系數μ=1×10-3Pa·s。

3 結果分析

3.1 臨界條件

分析不同溶腔水壓下拱頂控制點位移(圖6)可知,當溶腔水壓≤3.1 MPa時,拱頂控制點在富水溶腔作用下,先快速產生遠離溶洞的位移,并在0.6 s左右達到峰值3.432 4 mm(3.1 MPa)、2.826 6 mm(3.0 MPa),拱頂下沉量也同步的在0.6 s左右達到峰值2.880 8 mm(3.1 MPa)、2.475 5 mm(3.0 MPa)。在0.6 s之后,控制點處位移降低,并逐漸趨于平穩。

圖6 不同溶腔水壓下拱頂位移Fig.6 Displacements at vault under cavity pressures

當溶腔水壓大于3.1 MPa時,控制點處位移值持續增長,且增長速率逐漸增加,呈現明顯的不收斂特征。

圖7為不同溶腔水壓下拱頂控制點位移值。當溶腔水壓小于等于3.0 MPa時,拱頂控制點處x方向與y方向位移值很小,最大值分別為0.083 3 mm(y方向)和0.002 3 mm(x方向),且位移值隨溶腔水壓變化的幅度小。當溶腔水壓從3.0 MPa提升到3.1 MPa時,拱頂控制點位移量逐漸增加,提升幅度較為明顯。而當溶腔水壓從3.1 MPa提升至3.2 MPa時,拱頂位移增速達38.109 mm/MPa(y方向)和57.542 mm/MPa(x方向),遠高于溶腔水壓從3.0 MPa提升至3.1 MPa時拱頂控制點位移增速,增長幅度分別為204%和472%,且當溶腔水壓超過3.1 MPa后,位移增長速率基本保持不變,呈現不收斂的趨勢。

圖7 不同溶腔水壓下拱頂位移Fig.7 Displacements at vault under cavity pressures

通過分析不同溶腔水壓下拱頂控制點位移的演化規律可知,拱頂控制點位移隨溶腔水壓的增加而增加。當溶腔水壓較小時,隧道拱頂處圍巖能維持自身穩定,位移值小;當溶腔水壓接近臨界條件時,拱頂處圍巖在水壓作用下逐漸失去自穩能力,位移值增大,但增速較緩;當溶腔水壓超過臨界條件時,隧道拱頂圍巖失去自穩能力,位移量增長較快,位移增速提升顯著。

對比不同溶腔水壓下拱頂控制點位移演變規律認為該工況下,3.1 MPa為安全臨界水壓。將數值模擬結果和隧道與側部溶洞間巖層防突厚度計算公式[18]計算的臨界溶腔水壓3.446 MPa相比,臨界溶洞水壓降低幅度達11.1%。

3.2 臨界水壓下圍巖特性

分析已知溶腔水壓為3.1 MPa為臨界突水條件,為了解臨界突水圍巖特性,在隧道拱頂、拱肩、拱腰、起拱線及拱腳處共設置9個觀測點,以揭示臨界溶腔水壓下,圍巖位移分布規律。

分析圖8可知,當溶腔水壓為3.1 MPa時,隧道周邊控制點水平位移規律基本一致,均呈現在固定溶腔水壓,控制點位移值先增大后減小的趨勢。這是因為富水溶腔位于隧道左上角,隧道洞周控制點在水壓的作用下均產生遠離溶腔的位移。其中,位于拱頂的2號控制點的x方向位移值隨時間波動幅度最大,在0.64 s時達到最大值3.432 4 mm。0.64 s以后,洞周控制點位移值逐漸減少。

圖8 溶腔水壓為3.1 MPa時的控制點位移Fig.8 Displacements of control points under 3.1 MPa cavity pressure

當溶腔水壓為3.1 MPa時,洞周控制點y方向位移值演變規律與水平位移演變規律差異較大。分析圖9可知,在溶腔水壓作用下,隧道拱頂與拱肩處控制點y方向位移值最高位;兩側拱腰與右側拱腳處位移值較低,但與拱底及靠近溶腔一側的拱腳處控制點位移基本沒變化相比,位移波動幅度較大。

圖9 溶腔水壓為3.1 MPa時的控制點位移Fig.9 Displacements of control points under 3.1 MPa cavity pressure

因此,當溶腔位于隧道左上方時,可通過監測拱頂、拱肩處(尤其是靠近溶腔一側拱肩)結構位移來進行施工監測,以充分保障施工安全。

圖10展示了3.1 MPa(臨界水壓)與3.2 MPa(超過臨界值)下,拱頂控制點位移隨時間變化曲線。分析可知,x方向的位移分為三階段。第一階段,拱頂控制點x方向位移增長速率高,且兩種溶腔水壓作用下,位移值基本一致。該現象表明在臨界水壓與超過臨界值的水壓作用下,圍巖初期位移值變化差異性不大。第二階段為監測時間0.40～0.64 s,3.1 MPa水壓作用下,控制點x方向位移增長速率降低,位移逐漸達到峰值3.432 4 mm;3.2 MPa水壓作用下,控制點位移持續增長,且增長速率保持不變。第三階段為監測時間0.64 s以后,3.1 MPa水壓作用下,控制點位移在達到峰值后緩慢降低;而3.2 MPa溶腔水壓作用下,控制點位移繼續增長,且增長速率進一步提高,呈現明顯的不收斂趨勢。

圖10 拱頂控制點位移Fig.10 Displacements at vault

y方向的位移也可以分為三階段,但與x方向位移演化規律不同。y方向位移的三個階段可總結如下:緩慢增長段、快速增長段、位移失控段。在緩慢增長段,兩種溶腔水壓下,控制點位移增速均很低,位移幾乎沒有增長。這表明,在溶腔水壓作用初期圍巖仍能保持自穩。在快速增長段,控制點位移均快速增長,且增長速率幾乎一致。這表明在發生突涌水災害前,不同水壓下圍巖位移演化規律基本一致。在位移失控段,3.1 MPa溶腔水壓作用下,控制位移達到峰值后出現明顯的下降段,下降速率平緩;而3.2 MPa溶腔水壓作用下,控制點豎向位移先飆升一段,在經歷一段較為平緩的增長后,位移再次快速增長,直至圍巖失穩,發生突涌水災害。

圖11為3.2 MPa溶腔水壓工況下,達西流速度云圖。在發生裂隙型巖溶突涌水時,由于地層裂隙的存在,孔隙流體在擴散時,會找到一條連接路徑將地層裂隙、富水溶腔與隧道開挖面。突水過程中,富水溶腔流體出口及隧道開挖面流體入口處流體流速較快,達到1.4 m/s。隧道開挖面突水點集中出現在拱肩處,驗證了Wang等[8]提出的預警方案的合理性。由于富水溶腔底部高程位于隧道拱頂上方,發生突涌水災害時,與拱腳相比,拱頂處峰值流速云圖區域更大,在隧道邊界形成一高流速區域。

相比初始裂隙開度0.1 mm,發生突涌水災害時,裂隙網格模型最大孔徑為0.55 mm,最小裂隙孔徑為0.05 mm,如圖12所示。峰值孔徑裂隙出現在隧道開挖面與富水溶腔間,說明受高流體壓力激活、沖刷的影響,位于隧道開挖面與富水溶腔間的孔隙孔徑不斷發育、擴張,為水體流動提供條件。遠離隧道開挖面與富水溶腔的裂隙則在高應力作用下,逐漸趨于閉合,從而抑制流體流動。

4 結論

裂隙型巖溶突涌水是一種非常復雜的多物理場耦合過程。本文研究基于流-固-損傷全耦合及Comsol Multiphysics求解偏微分方程對突涌水臨界條件展開討論,為實際工程風險控制提供支撐。具體結論與建議如下。

(1)建立了流-固-損傷全耦合數值模型并進行計算分析,得出臨界溶腔水壓小于理論計算值,計算結果偏安全,更符合實際工程需求。

(2)裂隙會削弱巖體物理力學性能。在地勘階段,應充分考慮到巖體裂隙或裂隙發育區域對工程的危害,完備記錄探測所得裂隙及其產狀特征,以便施工安全控制。

(3)應根據富水溶腔與隧道斷面的相對位置關系設置位移監測點進行突涌水監測預警。當溶腔分布在隧道左上方時,宜在拱頂與靠近溶腔一側拱肩處設置監測點。

(4)裂隙型巖溶突涌水災害發生時,拱頂圍巖位移經歷三階段:緩慢增長段、快速增長段、位移失控段。

(5)裂隙型巖溶突涌水發生的本質是巖體裂隙在高水頭壓力作用下發育、擴張,最后失穩突水。溶腔水壓≤3.1 MPa時,圍巖位移增長穩定,達到峰值后會回落至穩定值;溶腔水壓≥3.1 MPa時,圍巖位移量陡增,且位移速率持續增長,呈現明顯的不收斂特征。