李天勝，何 川，周子寒，陳子全，包燁明，汪 波

(1.西南交通大學 交通隧道工程教育部重點實驗室，四川 成都 610031；2.中鐵十二局集團有限公司，山西 太原 030024)

四川省西部地區(qū)隧道賦存的地質(zhì)環(huán)境普遍存在高地應力、高地溫的特點[1]。在溫度應力的疊加效應下,隧道開挖后的環(huán)向應力將會顯著增大,巖爆烈度隨之增強[2-3]。然而,在以往設計階段的巖爆預測模型中,作為應力邊界條件來源的巖體初始地應力場反演結(jié)果由于未考慮豎向溫度應力的影響,存在隧道巖爆烈度預測偏低的情況。鑒于此,研究高地溫、高地應力環(huán)境下考慮豎向溫度應力的巖爆瞬態(tài)預測方法對指導類似隧道工程的修建具有重要意義。

基于隧址區(qū)地勘及地應力實測資料,采用反演手段可獲取巖體初始地應力場。針對這一領域,國內(nèi)外學者在以下方面開展了大量研究工作。汪波等[4]利用多元線性回歸方法結(jié)合最小二乘尋優(yōu)準則求解了隧址區(qū)巖體初始地應力場,該方法的特點是能夠確定反演唯一解,且運算速率較快。蒙偉等[5]同樣基于這一方法,進一步考慮側(cè)壓力系數(shù)對一次反演結(jié)果進行修正,提高了反演精度。另一方面,隨著智能算法的興起,其在地應力反演領域同樣應用廣泛。劉泉聲等[6]采用支持向量回歸優(yōu)化算法,減小了實測地應力隨機誤差對地應力反演的影響。裴書鋒等[7]利用進化神經(jīng)網(wǎng)絡-遺傳算法,提出考慮還原白鶴灘水電站巖體應力型破壞特征以修正反演結(jié)果的方法。運用智能算法可以較好解決實測值與計算值之間的非線性擬合問題,但決定擬合效果的數(shù)據(jù)訓練樣本需求很大,因此計算效率較低。

根據(jù)反演結(jié)果,分析總結(jié)初始地應力場的分布特征對指導地下工程的修建意義重大。陳世杰等[8]討論了金川地下廠房斷層帶內(nèi)外局部地應力場的分布規(guī)律,認為靠近斷層帶應力量值增大,而斷層內(nèi)部應力減小。周子寒等[9]分析了花崗質(zhì)侵入巖地層初始地應力的分布特征,發(fā)現(xiàn)接近侵入體前沿區(qū)域巖體應力場集中現(xiàn)象顯著。Li等[10]通過分析地應力反演結(jié)果得出隧道穿越斷層或背斜時,主應力量值會突然下降的結(jié)論。

高能地質(zhì)環(huán)境下硬巖隧道往往會發(fā)生巖體瞬時破壞(巖爆)現(xiàn)象,國內(nèi)外學者在巖爆預測方面進行了相關研究:嚴健等[11]通過反演隧址區(qū)高精度地應力場,對巖爆發(fā)生類型、破壞形態(tài)進行了機理性探索。Meng等[12]研究了高地溫高地應力環(huán)境下的巖體初始應力場,認為高地溫環(huán)境下巖體初始地應力更大,可能會加劇巖爆程度。嚴健等[13]根據(jù)地應力場反演結(jié)果,采用熱力耦合模型計算了拉林鐵路桑珠嶺隧道的巖爆烈度,但其應力邊界條件未考慮溫度應力的影響。

綜上所述,利用反演求解的隧址區(qū)初始地應力場作為巖爆模型的應力邊界條件求解巖爆烈度是巖爆預測的常規(guī)手段。然而,在以往的研究中,由于地應力反演過程未考慮高地溫環(huán)境下豎向地溫梯度產(chǎn)生的溫度應力,進而造成巖爆預測結(jié)果“失真”。另一方面,采用熱力耦合本構(gòu)計算巖爆烈度也僅考慮穩(wěn)態(tài)結(jié)果,開挖瞬間的瞬態(tài)巖爆烈度還未見報道。本文以康定某隧道為依托,在水壓致裂法的豎向應力分量中疊加地溫熱應力,結(jié)合多元線性回歸法計算效率較高的特點,反演了隧址區(qū)初始地應力場。在此基礎上,以反演結(jié)果作為巖爆預測模型的應力邊界條件,考慮隧道洞壁與洞內(nèi)環(huán)境的對流換熱過程,采用熱力耦合本構(gòu)模型求解了巖爆段落的瞬態(tài)巖爆烈度。研究形成了適用于高地溫高地應力隧道初始地應力反演到瞬態(tài)巖爆預測的方法。研究成果可指導類似隧道工程的設計和施工。

1 工程概況

1.1 地形地貌及地層巖性

康定某隧道位于四川盆地與青藏高原過渡的西南邊緣。受岷江、大渡河等水系的強烈切割,地形高差大,溝壑密布,山嶺縱橫。隧址區(qū)地表高程3 460～4 730 m,為典型高原地貌。隧道進口位于康定市折多塘村西面山坡,出口位于瓦澤鄉(xiāng)水橋村318國道附近,隧道全長約20.8 km,最大埋深約1 215 m,地層巖性主要分布有黑云母二長花崗巖、石英砂巖、板巖。

1.2 研究區(qū)域地應力特征

隧道處于松潘—甘孜褶皺帶,SN向川滇構(gòu)造帶和NW向鮮水河構(gòu)造帶在區(qū)域內(nèi)形成中國西南著名的Y字形構(gòu)造格局,隧址區(qū)大地構(gòu)造背景處于Y字形構(gòu)造三岔口交接部位的NW側(cè)。隧道縱斷面見圖1。由圖1可以看出,強烈的地質(zhì)構(gòu)造作用下,隧址區(qū)共計12處斷層,其中,折多塘斷層位于隧道進口附近,影響范圍為3～5 km,因而僅以實際巖性表示。隧址區(qū)現(xiàn)今主應力方向為NW-NWW向,隧道大致走向為NW向,主應力方向與隧道走向近乎平行。為探明隧址區(qū)地應力分布情況,在地勘階段采用水壓致裂法進行了地應力測量,共計鉆孔5處,地應力實測結(jié)果見表1。

表1 水壓致裂法地應力測量結(jié)果

圖1 隧道縱斷面

由表1可知,隧址區(qū)實測地應力中最大水平主應力σH、最小水平主應力σh和豎向應力σV的范圍分別為4.1～32.2、3.9～17.4、3.7～16.8 MPa。受強烈的地質(zhì)構(gòu)造運動影響,隧址區(qū)初始地應力以水平向構(gòu)造應力為主導,5處鉆孔均呈現(xiàn)σH>σh>σV的分布特征。同時,實測地應力量值隨著測點埋深的增大而增大。

最大水平主應力σH的方位角范圍為N11°W～N68°W。其中,DZ-ZDS-S-02、DZ-DSFA-01、DZ-ZDS-S-03、DZ-SK-ZDS-09、DZ-ZDS-S-04鉆孔σH的平均方位角分別為N67°W、N51°W、N50°W、N28°W、N52°W。地應力測量結(jié)果與大地構(gòu)造走向基本一致。

1.3 隧址區(qū)高地溫測量

隧道穿越玉龍希斷裂有溫泉出露,同時隧道進口穿越折多塘導熱斷裂,因此針對隧道賦存的高地溫環(huán)境,在地勘階段共計進行了3處測溫鉆孔測量地溫梯度,見圖2。

圖2 隧道地溫鉆孔分布圖

2 隧址區(qū)地應力反演原理

2.1 高地溫環(huán)境水壓致裂法測量原理

要探討高地溫環(huán)境下的地應力反演原理,首先需要分析大地溫度梯度分布及水壓致裂法的測量原理。地表往下不受大氣溫度影響的地層溫度會隨著深度的增加而增加,其增長幅度稱為地溫梯度,如圖2中隧道的實測平均地溫梯度為3.1 ℃/100 m。而巖體受到大地熱流的影響將會發(fā)生熱膨脹,在周圍巖體的約束下產(chǎn)生熱應力[14],即

σT=αβEZ

(1)

式中:σT為巖體在高地溫環(huán)境下產(chǎn)生的熱應力;α為地溫梯度;β為巖體體膨脹系數(shù);E為巖體彈性模量;Z=H-h,其中,H為埋深,h為變溫帶與恒溫帶高度之和,一般約為30 m。

假定高溫膨脹無特定方向,則按式(1)產(chǎn)生的高溫熱應力呈靜水壓力的形式。

按照水壓致裂法的原理[15],測點最大水平主應力σH和最小水平主應力σh的測量結(jié)果為初始地應力的真實結(jié)果,見圖3,因此測量數(shù)據(jù)包含了熱應力分量和構(gòu)造應力分量

(2)

圖3 高地溫環(huán)境水壓致裂法測量原理

(3)

水壓致裂法豎向應力卻并非來自實際的測量結(jié)果,其是根據(jù)埋深對應的上覆荷載估計所得

σV=γH

(4)

顯然豎向應力忽略了高地溫環(huán)境下溫度熱應力分量,實際的豎向應力σ′V應為

σ′V=γH+σT

(5)

那么,以往在地應力反演過程中并未考慮豎向溫度應力分量的影響,僅以σH、σh、σV作為實測值進行數(shù)據(jù)擬合反演顯然是不妥的。

2.2 高地溫環(huán)境地應力反演原理

首先將實測應力中包含溫度應力的最大水平主應力與最小水平主應力剔除溫度應力分量。以未含溫度應力的水平應力參與地應力反演,則剔除溫度應力后的兩項水平主應力分別為

(6)

(7)

σv′w′=αv′vαw′wσvw

(8)

式中:σv′w′為三維數(shù)值模型坐標系的應力分量;σvw為實測應力坐標系的應力分量;v,w=x,y,z;αv′v和αw′w為坐標轉(zhuǎn)化系數(shù),如αv′w=cos(v′,w)。

(9)

式中:b0為常數(shù)項;bi為第i種邊界條件工況下的多元線性回歸系數(shù);n為工況總數(shù)。

利用統(tǒng)計學中最小二乘法原理,可以求解式(9)中的回歸系數(shù)。假定有m個測點,每個測點有l(wèi)個應力分量,則最小二乘法的殘差平方和Qc為

(10)

式中:σkj為第k個測點的第j個應力分量實測值。

(11)

(12)

2.3 三維數(shù)值模型的建立

基于克里金插值法,提取地勘階段獲取的隧址區(qū)平面圖三維坐標信息,導入ANSYS軟件中建立三維數(shù)值模型的頂面。為避免計算結(jié)果受模型尺寸效應的影響,三維模型的高度方向尺寸應與寬度方向尺寸基本相當。鑒于平面圖寬度方向的有效尺寸為4 000 m,因此,在高度方向上,三維模型的底面取至隧道進口高程以下2 000 m。模型長度方向包含隧道全部線路,約21 000 m。由于模型頂面及斷層的不規(guī)則分布,對部分區(qū)域的網(wǎng)格劃分進行加密處理。巖體材料的物理力學參數(shù)見表2。建立的三維數(shù)值模型見圖4。

表2 巖體物理力學參數(shù)

圖4 三維數(shù)值模型(單位:m)

2.4 模型邊界條件

考慮x、y方向的構(gòu)造作用,xoy平面的剪切構(gòu)造作用。對三維數(shù)值模型(圖5)施加x方向的水平均布荷載Fx和y方向的水平均布荷載Fy,在xoy平面施加剪切荷載;按照張奇華等[16]的方法對三維數(shù)值模型的側(cè)面施加位移邊界條件(Sx、Sy)實現(xiàn),即Sx=4.0 mm,Sy=21 mm。自重通過對三維模型施加重力加速度模擬。

圖5 模型邊界條件示意

邊界位移約束條件:頂面為自由邊界,底面施加法向位移約束,側(cè)面除需施加均布荷載面外,其余側(cè)面均施加法向位移約束。

2.5 反演結(jié)果及精度驗證

基于反演原理,通過施加圖5所示的邊界條件后,利用式(11)最小二乘法尋優(yōu)方法可求解得到剔除地溫熱應力的回歸方程為

σhg=1.32+0.9σzz-4.0σx,jb+9.5σy,jb+4.0σxoy

(13)

式中:σzz為自重荷載工況的計算值;σx,jb為x向水平均布荷載的計算值;σy,jb為y向水平均布荷載的計算值;σxoy為xoy面剪切荷載的計算值。

按照回歸統(tǒng)計原理,對回歸結(jié)果進行顯著性檢驗,F顯著性統(tǒng)計量F=1.8×10-24小于顯著性水平0.01,復相關系數(shù)R=0.91≈1,這表明回歸值與剔除熱應力的實測值之間相關性較高,回歸效果較好。

為驗證本次反演精度,以埋深最大的DZ-DSFA-01鉆孔為例,對比實測值與計算值,見表3,由于切向應力τxy量值較小,因此未參與討論。巖體熱力學參數(shù)見表4,則地溫熱應力σT=0.031×9×10-6×26×103×(H-30)。從表3可知,DZ-DSFA-01鉆孔σx分量的平均相對誤差為17.8%,σy分量的平均相對誤差為19.8%,σz分量的平均相對誤差為7.9%,總體平均誤差為15.2%,反演精度較為理想。由于該隧道地應力測點較多,限于文章篇幅有限,表5列出了其余鉆孔的平均相對誤差。由此可知,所有測點的平均相對誤差為18.2%,能夠控制在20%以下,滿足工程建設的需要。

表3 DZ-DSFA-01鉆孔實測地應力值與回歸值對比

表4 巖體熱力學參數(shù)

表5 各鉆孔平均相對誤差

注:σx為x方向地應力值;σy為y方向地應力值;σz為z方向地應力值;左側(cè)應力值代表剔除地溫熱應力的結(jié)果;右側(cè)應力值代表包含地溫熱應力的結(jié)果,相對誤差計算公式為:[(|實測值-回歸值|)/實測值]×100%,其中應力值均采用表中右側(cè)結(jié)果。

DZ-DSFA-01鉆孔各測點的地溫應力范圍為0.9～4.6 MPa,其平均占水平向應力的20.1%,平均占豎向應力的20.5%,由此可見溫度應力的占比不容忽視,其將會對巖爆預測產(chǎn)生影響。

3 瞬態(tài)巖爆預測

3.1 巖爆模型及邊界條件

根據(jù)地勘信息,隧道在里程DK281+880—DK282+400,DK283+000—DK283+860,DK284+260—DK284+720,DK286+200—DK286+700四個區(qū)段(圖1)圍巖級別為Ⅲ級,在高地應力環(huán)境下存在發(fā)生中等巖爆的風險。

為探討考慮地溫熱應力下瞬態(tài)巖爆烈度問題,采用COMSOL Multiphysics軟件熱力耦合模塊建立巖爆預測數(shù)值模型,見圖6。數(shù)值模型范圍取至隧道洞徑的3倍,長度約為90 m,高度約為75 m。

圖6 巖爆預測數(shù)值模型網(wǎng)格

通過式(13)可以查詢DK281+880—DK286+700范圍內(nèi)的地應力回歸值,再代入式(12)中,即可求出考慮地溫熱應力的初始地應力場。將三維數(shù)值模型坐標系下的地應力場轉(zhuǎn)化為主應力的形式為

(14)

DK281+880—DK286+700范圍內(nèi)考慮地溫熱應力的巖體初始地應力場見圖7。由圖7可知,當埋深增大以后,隧道走線賦存的地應力環(huán)境轉(zhuǎn)為以豎向應力為主導,整體呈現(xiàn)σV>σH>σh的特征。

圖7 隧道初始地應力場

以隧道里程處對應的初始地應力場作為應力邊界條件見圖8(a),溫度邊界條件見圖8(c),包括模型初始溫度T0和模型邊界溫度Tb。若是未考慮豎向溫度應力的邊界條件如圖8(b)所示。圖8(a)和8(b)中的帶下標1和2的應力分量分別為考慮溫度應力反演方法的地應力值和未考慮豎向溫度應力反演方法的地應力值?？紤]隧道洞壁與洞內(nèi)環(huán)境的對流換熱過程,對洞壁邊界設定對流換熱系數(shù)qd,熱力學計算參數(shù)見表4。

圖8 模型熱力耦合邊界條件

圍巖初始溫度T0和模型邊界溫度Tb后設為相同值,即T0=Tb,并考慮不利工況的影響取為表4中的最大值。洞內(nèi)環(huán)境溫度取洞口全年平均溫度(0.3 ℃,地勘資料中洞口氣溫統(tǒng)計見表6)與洞內(nèi)施工允許的最高溫度[13](28 ℃)的平均值14.2 ℃?？紤]熱力耦合瞬態(tài)求解,得到的應力解帶入陶振宇判據(jù)[17]進行巖爆烈度預測,見表7。表7中,Rc為單軸抗壓強度,里程DK281+880—DK284+720取135 MPa,里程DK286+200—DK286+700取130 MPa,σ1為最大主應力。

表6 洞口近30年月平均氣溫表

表7 陶振宇巖爆判據(jù)[17]

隧道考慮豎向地溫應力的瞬態(tài)巖爆預測結(jié)果見圖9,由圖9可知,與地勘預測的巖爆烈度不同,在考慮豎向地溫應力后,巖爆烈度顯著增大,除DK281+880—DK281+900范圍處于中等巖爆區(qū)外,其余區(qū)域均處于強烈?guī)r爆區(qū)。

圖9 不同方法的巖爆烈度預測結(jié)果

4 不同巖爆預測方法的對比分析

以往傳統(tǒng)地應力反演方法未考慮地溫熱應力的影響,反演直接以計算值按照多元回歸方法逼近實測值σH、σh、σV(最大水平主應力已包含溫度應力,豎向應力未包含溫度應力),欠缺了豎向地溫應力分量[18-19]。同時,以往求解僅考慮熱力耦合的穩(wěn)態(tài)解,對于開挖瞬間的瞬態(tài)巖爆烈度還未討論[20-22]?；诖?共設計3種計算工況巖爆強度(表8),選取四個巖爆區(qū)段中最大埋深里程DK282+400、DK283+130、DK284+320、DK286+700為算例進行分析,見表9。

表8 不同工況邊界和求解條件

表9 不同工況應力邊界條件

按照多元線性回歸原理,按傳統(tǒng)反演方法再進行一次回歸,得到未考慮豎向地溫應力分量的回歸方程為

σhg=1.7+0.95σzz-2.1σx, jb+12.1σy, jb+6.0σxoy

(15)

工況3已在前文進行了求解,工況2與工況3唯一的不同是其按照穩(wěn)態(tài)求解。工況1按照式(15)可求出各里程對應的巖體初始地應力作為應力邊界條件,見圖8(b),溫度邊界條件同樣按照圖8(c)所示進行施加,考慮穩(wěn)態(tài)進行求解。

圖10為里程DK286+700在三種工況下的計算云圖。工況1由于未考慮豎向溫度熱應力的疊加效應,應力邊界條件的量值小于其余工況,熱力耦合得到的洞壁第一主應力最大值為57 MPa,工況2考慮豎向溫度熱應力的疊加效應,洞壁第一主應力最大值顯著增大至86 MPa,工況3考慮瞬態(tài)求解后洞壁第一主應力最大值繼續(xù)增大,達到102 MPa。這說明,開挖瞬間圍巖初始溫度與洞內(nèi)環(huán)境溫度差值處于峰值狀態(tài),此時由溫差導致的熱應力釋放將增大洞壁處的應力集中效應,巖爆烈度將會增強。

圖10 里程DK286+700熱力耦合第一主應力云圖

四個算例在三種工況下的計算結(jié)果見圖11。由圖11可知,工況1未考慮豎向地溫應力穩(wěn)態(tài)解的巖爆預測,除里程DK286+700巖爆烈度指標小于2.5屬強烈?guī)r爆外,其余里程巖爆烈度指標均大于2.5,屬中等巖爆。工況2考慮豎向地溫應力穩(wěn)態(tài)解的巖爆預測,巖爆烈度顯著增大,里程DK282+400處巖爆烈度指數(shù)由3.1減小至1.9(減幅38.7%),里程DK283+130處巖爆烈度指標由2.8減小至1.7(減幅39.3%),里程DK284+320處巖爆烈度指標由2.8減小至1.9(減幅32.1%),里程DK286+700處巖爆烈度指標由2.3減小至1.5(減幅34.8%),巖爆烈度指標的減小代表巖爆烈度的增加,因此四處里程的巖爆烈度平均增幅為36.2%。工況3考慮豎向地溫應力瞬態(tài)解的巖爆預測,巖爆烈度繼續(xù)增強,較工況2巖爆烈度平均增幅為14.2%,較工況1巖爆烈度平均增幅為45.3%。

圖11 四個算例的巖爆預測結(jié)果

此外,按照表8中的三種工況,還計算了隧道四段巖爆區(qū)域的烈度指數(shù),見圖9。由圖9可知,按照工況1即傳統(tǒng)的巖爆預測手段,隧道幾乎所有巖爆區(qū)段都將落在中等巖爆區(qū)之內(nèi),這與地勘報告的結(jié)論相似。工況2考慮豎向溫度應力疊加效應后,除DK281+880—DK282+100段落在中等巖爆區(qū)內(nèi),其余區(qū)段的巖爆預測結(jié)果均為強烈?guī)r爆。工況3考慮瞬態(tài)解后,中等巖爆區(qū)域更短,僅DK281+880—DK281+900,其余區(qū)段全為強烈?guī)r爆,且?guī)r爆烈度指數(shù)小于工況2的結(jié)果。

考慮溫度熱應力的巖爆瞬態(tài)預測方法見圖12。

圖12 考慮溫度熱應力的瞬態(tài)巖爆預測方法

Step1實測水平主應力剔除溫度熱應力后與豎向應力一起參與多元線性回歸得到反演回歸值,將回歸值疊加溫度熱應力后得到巖體初始地應力場。

Step2按照隧道里程提取巖體初始地應力場作為巖爆預測模型的應力邊界條件,考慮洞壁與洞內(nèi)的對流換熱過程,以瞬態(tài)求解隧道開挖后的洞壁應力值,將其代入巖爆判據(jù)中進行巖爆預測。

5 結(jié)論

1)水壓致裂法的豎向應力(σV=γH)是根據(jù)上覆埋深的巖體重力估算而來,未包含溫度熱應力,因此針對高地溫環(huán)境的地應力反演應考慮豎向溫度熱應力的疊加效應。

2)計算表明,康定某隧道由于豎向地溫梯度導致的熱應力約占水平構(gòu)造應力的20.1%,約占豎向應力的20.5%,未考慮溫度熱應力的巖爆預測將會弱化巖爆烈度。

3)未考慮地溫熱應力的穩(wěn)態(tài)巖爆預測結(jié)果與地勘結(jié)論吻合(均為中等巖爆)。通過算例的對比分析可知,考慮地溫熱應力的穩(wěn)態(tài)巖爆預測烈度較傳統(tǒng)巖爆預測方法平均增強了約36.2%,而考慮地溫熱應力的瞬態(tài)巖爆預測烈度繼續(xù)增強了約14.2%。

4)提出的考慮地溫應力的熱力耦合瞬態(tài)巖爆預測方法可修正巖爆預測烈度,對指導相似高地溫隧道的巖爆預測具有一定的參考和指導意義。

5)提出的反演回歸應力與溫度應力疊加的方式獲取考慮地溫熱應力的隧址區(qū)地應力場雖原理簡便,但過程煩瑣。采用熱力耦合模型直接參與地應力場反演更加直接,這方面的研究有待下一步開展。