非圓形隧洞解析應(yīng)力空間分布特征

杜小洲， 胡志平， 王振林， 安學旭*， 陶磊

(1.陜西省引漢濟渭工程建設(shè)有限公司， 西安 710024； 2.長安大學建筑工程學院， 西安 710061)

深埋高地應(yīng)力條件下隧洞開挖經(jīng)常伴隨有圍巖巖爆、剝落、大變形等一系列災(zāi)變現(xiàn)象，已給地下隧洞安全建設(shè)構(gòu)成了重大挑戰(zhàn)[1-2]，雖然引起這些災(zāi)變現(xiàn)象的因素較多，但現(xiàn)有理論和實踐證明這些災(zāi)害與隧洞掌子面推進過程圍巖重分布空間應(yīng)力分布密切相關(guān)[3]，因此，有必要對深埋復(fù)雜應(yīng)力場環(huán)境中隧洞圍巖空間應(yīng)力分布特征進行研究。

關(guān)于隧洞圍巖應(yīng)力分布特征，學者們采用了各種方法進行了大量研究。理論方面，已有諸多學者采用彈塑性理論對非均勻應(yīng)力場中圓形隧洞斷面應(yīng)力分布特征及塑性區(qū)形狀進行了各種研究[4-9]。對于復(fù)雜非圓形隧洞斷面，Lü等[10]、李巖松等[11]、李元輝等[12]通過復(fù)變函數(shù)理論對復(fù)雜隧洞斷面圍巖應(yīng)力分布特征進行了分析。目前雖然解析理論較多，但多集中在二維平面中。數(shù)值模擬方面，因其成本低、可重復(fù)以及較強的信息輸出能力已被廣泛地應(yīng)用于去分析隧洞開挖過程中出現(xiàn)的各種工程問題，通常采用邊界應(yīng)力釋放法、巖芯剛度折減法或巖芯置換技術(shù)等效手段來再現(xiàn)掌子面推進過程中的空間效應(yīng)作用[13-14]。如陶志剛等[15]、蔡武強等[16]等通過三維數(shù)值模型分析了深埋巖體隧洞開挖過程中圍巖應(yīng)力、位移以及掌子面的非線性擠出效應(yīng)。雖然數(shù)值模擬具有諸多優(yōu)點，由于隧洞開挖是一個三維空間問題，對三維模型進行過細網(wǎng)格劃分必然引起較高計算時間成本。試驗方面，一些學者在制作物理試驗?zāi)Ｐ蜁r，通過在預(yù)定位置預(yù)埋與開挖洞型一致的柱體，等試件干燥后將其拔出或頂出，然后再進行加載，或在預(yù)定位置通過手動或利用小型機械設(shè)備開挖出一定形狀的隧洞來模擬并分析掌子面推進過程中圍巖應(yīng)力空間分布特征[17-18]。三維模型試驗研究隧洞掌子面空間效應(yīng)雖效果較好，但其缺陷也較顯著，如造價高、可重復(fù)利用性差等。

由上述分析可知，對于隧洞圍巖應(yīng)力分布特征的研究雖然方法較多，但各有局限?，F(xiàn)基于收斂約束法，將非圓形隧洞圍巖應(yīng)力及位移解析理論與可視化軟件Tecplot通過程序相結(jié)合，提出一種高效、直觀的隧洞解析應(yīng)力空間分布方法，并采用該方法對引漢濟渭工程嶺北隧洞開挖過程中掌子面附近圍巖空間應(yīng)力分布特征以及不同初始主應(yīng)力場水平傾角下的隧洞圍巖空間應(yīng)力分布特征進行研究。

1 隧洞解析應(yīng)力空間分布方法

隧洞開挖過程中由于開挖面的卸荷作用使得圍巖中產(chǎn)生應(yīng)力重分布現(xiàn)象并伴隨位移釋放，且整個過程不會在開挖瞬間立即完成，而是在掌子面推進過程中逐漸完成，如果在無支護條件下，當掌子面推進一定距離后整個位移釋放才完成[19]。目前，這個過程可通過收斂約束法中縱向變形剖面(longitudinal deformation profile,LDP)與圍巖特征(surrounding rock characteristic，GRC)曲線來表示[20]，縱向變形曲線能反映與掌子面不同距離處隧洞斷面圍巖位移，一般可通過實測數(shù)據(jù)或數(shù)值模擬等方法獲取。圍巖特征曲線反映了隧洞斷面圍巖位移與虛擬支護力之間的關(guān)系[21]，本文中采用隧洞斷面位移釋放系數(shù)與虛擬支護力的之間的關(guān)系來反映，該曲線可通過考慮位移釋放系數(shù)的隧洞圍巖應(yīng)力及位移解析理論獲得。

解析應(yīng)力空間分布模型構(gòu)建：首先采用相同的變量位移釋放系數(shù)，建立隧洞LDP曲線與GRC曲線之間的關(guān)系如圖1所示，然后將平面空間位移釋放系數(shù)對應(yīng)至三維空間中支護與掌子面之間的距離，確定三維隧洞空間任意點處圍巖應(yīng)力及位移值，最后通過可視化軟件Tecplot輸出三維空間應(yīng)力分布圖，可供實際工程支護設(shè)計與施工提供參考。

為了使該方法通用化，具有可復(fù)制性，將解析理論通過編程在MATLAB軟件中實現(xiàn)快速計算，并將圍巖力學計算結(jié)果數(shù)據(jù)以及圍巖縱向變形曲線數(shù)據(jù)導(dǎo)入可視化軟件Tecplot，即可輸出三維隧洞空間應(yīng)力圖，整個計算過程只需要幾分鐘即可完成，計算流程圖如圖2所示。

η(x)=ur/urmax, ur為掌子面附近任意隧洞斷面圍巖位移，urmax為無支護條件下隧洞斷面最大圍巖位移；Pi為虛擬支護力，MPa; x(i)為距離掌子面的距離，m圖1 隧洞LDP與GRC曲線之間的關(guān)系Fig.1 Relationship between LDP and GRC curves

α為初始主應(yīng)力場水平傾角，(°)；σ1和σ3分別為最大和最小主應(yīng)力，MPa；E1和E2分別為圍巖與襯砌的楊氏模量，MPa；υ1和υ2分別為圍巖與襯砌的泊松比圖2 隧洞圍巖解析應(yīng)力空間分布方法流程圖Fig.2 Flow chart of analytical stress spatial distribution method for tunnel surrounding rock

2 考慮位移釋放系數(shù)的隧洞解析方程

2.1 隧洞圍巖位移解析表達

對于如圖3(a)所示的隧洞斷面豎向主軸與其主應(yīng)力場非平行的隧洞，可通過三角函數(shù)分解為如圖3(b)所示的力學模型。

σy為初始主應(yīng)力場水平應(yīng)力分量，MPa；σx為初始主應(yīng)力場豎向應(yīng)力分量；τxy為初始主應(yīng)力場剪切分量，MPa圖3 非圓形隧洞圍巖應(yīng)力場Fig.3 Surrounding rock stress field of a non-circular tunnel

L21表示z平面中支護外邊界，L11表示z平面中支護內(nèi)邊界；L22表示映射ζ平面中支護外邊界，L12表示映射ζ平面中支護內(nèi)邊界；θ為映射ζ平面中極坐標角度,(°)；ρ為映射ζ平面中極坐標半徑，mm；Ro為映射ζ平面中支護內(nèi)邊界半徑，mm；i為復(fù)數(shù)圖4 z平面上隧洞斷面形狀映射為ζ平面上單位圓周Fig.4 Tunnel section shape on the z plane is mapped to the unit circumference on the ζ plane

如圖4所示，對于單連通域內(nèi)復(fù)雜隧洞斷面形狀可通過保角映射轉(zhuǎn)換為單位圓外域，映射函數(shù)可以表示為

(1)

式(1)中：ω(ζ)為關(guān)于變量ζ的函數(shù)；R為反映隧洞大小的系數(shù)；ck為反映隧洞形狀的系數(shù)。

實際工程中，隧洞斷面雖然形態(tài)較多，但大多關(guān)于豎軸對稱，因此，式(1)中映射函數(shù)系數(shù)ck可用常實數(shù)表示。

對于無支護隧洞，圍巖應(yīng)力求解需通過兩個解析函數(shù)獲得，即

φ1(ζ)=Γ1ω(ζ)+φo(ζ)

(2)

ψ1(ζ)=(Γ2+iτxy)ω(ζ)+ψo(ζ)

(3)

式中：Г1=σx(1+λ)/4，Г2=σx(λ-1)/2；λ=σx/σy；φo(ζ)和ψo(ζ)分別為圓外解析函數(shù)，φo(ζ)、ψo(ζ)。

可以表示為

(4)

(5)

(6)

式中：ɑk(k=1,2,3,…)為解析函數(shù)的系數(shù)；j=2,3,…,n-2；系數(shù)lk求解見參考文獻[10]。

由于模型應(yīng)力邊界條件中存在非對稱初始主應(yīng)力場剪應(yīng)力分量，因此式(4)解析函數(shù)中系數(shù)ɑk一定是復(fù)數(shù)。ɑk可通過線性方程組求解，即

Xikak=Bi，i,k=1,2,…,n-2

(7)

式(7)中：當i=k時，Xik=kLi+k+1，若i+k+1>n，Xik=0；當i≠k時，Xik=kL2i+1-1，若i+k+1>n，Xik=-1。Bi列矩陣為

(8)

式(8)中：i=2,3,…,n-2。

由式(1)、式(4)及式(5)可知無支護隧洞圍巖中任意點處位移為

(9)

式(9)中：μ1和ν1分別為豎向和水平總位移，mm；G1=E1/[2(1+υ1)]；κ1=3-4υ1；E1和υ1分別為圍巖楊氏模量與泊松比。

隧洞開挖過程中隨著掌子面的推進，圍巖中會產(chǎn)生應(yīng)力重分布，圍巖中位移會逐漸釋放，在支護前假定隧洞圍巖產(chǎn)生的位移為

μ2+iv2=η(μ1+iv1)

(10)

式(10)中：η為位移釋放系數(shù)，0≤η≤1.0；μ2和ν2分別為支護前圍巖中產(chǎn)生的豎向與水平位移，mm。

支護施工后，由于圍巖位移的進一步釋放使圍巖與支護之間會產(chǎn)生相互作用，由于支護約束作用在圍巖中產(chǎn)生的位移可表示為

(11)

式(11)中：μ3和ν3分別為支護約束作用在圍巖中產(chǎn)生的豎向和水平位移，mm；Φ2(ζ)和ψ2(ζ)分別為支護對圍巖的約束作用的圓外解析函數(shù)，可表示為

(12)

(13)

式中：bo、bk、do、dk均為待求復(fù)系數(shù)。考慮到支護對隧洞無窮遠處的影響較小，因此，認為在無窮遠支護對圍巖的作用為0，由此可得k1bo=k1(bo1+ibo2)=conj(do)= (do1-ido2)，conj(i)表示對系數(shù)i共軛。

支護在圍巖作用下任意點處位移可表示為

(14)

式(14)中：G2=E2/[2(1+υ2)],κ2=3-4υ2，E2和υ2分別為支護結(jié)構(gòu)的楊氏模量和泊松比；μ4和ν4分別為支護中產(chǎn)生的豎向和水平位移，mm；Φ3(ζ)和ψ3(ζ)分別為支護在圍巖作用下的圓外解析函數(shù)，可表示為

(15)

(16)

式中：po、ek、fk、qo、gk、hk均為待求復(fù)系數(shù)。

2.2 邊界連續(xù)條件

為了使得解析函數(shù)滿足應(yīng)力及位移邊界連續(xù)條件，則解析函數(shù)[式(12)、式 (13)、式(15)、式(16)]在支護外邊界L21處應(yīng)符合邊界應(yīng)力及位移連續(xù)條件，在支護內(nèi)邊界L11處應(yīng)符合邊界應(yīng)力連續(xù)條件。

在支護外邊界L21處由隧洞徑向應(yīng)力連續(xù)條件可得

(17)

式(17)中：σ=eiθ。

在支護邊界L21處由隧洞徑向位移連續(xù)條件可知

(μ1+iv1)(1-η)+(μ3+iv3)=μ4+iv4

(18)

將式(9)、式(10)、式(11)、式(14)代入式(18)化解可得

(19)

將ζ平面中支護內(nèi)邊界用Roσ表示，由支護內(nèi)邊界處L11處應(yīng)力連續(xù)條件可知

(20)

將式(17)、式(19)及式(20)左右兩側(cè)ζ±j(j=0,1,…,N)系數(shù)相同的實部和虛部整理聯(lián)立方程可獲得12N+6個方程和12N+6個未知量，這些線性方程組可統(tǒng)一表示為

Aijxj=Fi

(21)

式(21)中:Aij為方程組左側(cè)項系數(shù)矩陣；Fi為方程組右側(cè)項列矩陣，整個方程組在MATLAB中通過程序進行求解。

2.3圍巖應(yīng)力及位移求解

2.3.1 支護

將式(21)求解的變量代入式(12)、式(13)、式(15)、式(16)解析函數(shù)，結(jié)合式(1)可得

σpL+σθL=Re[φ′|3(ζ)/ω′(ζ)]

(22)

(23)

(24)

式中：σθL、σρL、τρθL分別為支護環(huán)向、徑向及剪切應(yīng)力，MPa。

2.3.2 圍巖

同理將式(21)求解的變量代入式(13)和式(14)，可得圍巖中應(yīng)力計算公式為

σpR+σθR=4Re[φ′|4(ζ)/ω′(ζ)]

(25)

(26)

φ4(ζ)=φ1(ζ)+φ2(ζ)

(27)

ψ4(ζ)=ψ1(ζ)+ψ2(ζ)

(28)

(29)

式中：σθR、σρR、τρθR分別為圍巖中環(huán)向、徑向及剪切應(yīng)力，MPa。

上述解析理論中，N越大所求得的計算結(jié)果越精確，實際上N>50基本可以滿足多數(shù)工程計算精度要求，本文中N=120。

3 算列分析

3.1 工程概況

引漢濟渭工程是陜西境內(nèi)重大水利設(shè)施項目，隧洞全長81.779 m，其中秦嶺嶺北硬巖段隧洞地質(zhì)環(huán)境比較突出，埋深最深達2 012 m，巖石平均抗壓強度為170 MPa，隧洞開挖過程中發(fā)生巖爆、大變形、巖體剝落、卡機等諸多工程災(zāi)變問題。為了給隧洞設(shè)計及施工提供依據(jù)，在嶺北6號實驗洞附近進行了地應(yīng)力測試，測試結(jié)果如表1所示[22]。圍巖楊氏模量與泊松比參數(shù)分別采用E1=64.5 GPa，υ1=0.19，支護楊氏模量與泊松比參數(shù)分別采用E2=30 GPa，υ2=0.20。隧洞斷面形式如圖5所示，斷面映射函數(shù)通過復(fù)合形優(yōu)化算法獲得，具體公式為

z=ω(ζ)=3.983 77(ζ-0.051 9+

0.001 7ζ-1+0.039 1ζ-2-0.042 2ζ-3+

0.022 9ζ-4-0.001 1ζ-5-0.007 1ζ-6)

(30)

表1 6#試驗洞地應(yīng)力測試結(jié)果Table 1 Situ stress test results of No. 6 test tunnel

圖5 引漢濟渭工程隧洞斷面Fig.5 Tunnel section of Hanjiang to Weihe River valley water diversion project

3.2 隧洞斷面空間應(yīng)力分布特征

根據(jù)3.1節(jié)隧洞斷面形狀及基本力學參數(shù)，通過數(shù)值模擬法獲得了隧洞頂部位移釋放系數(shù)與掌子面不同距離之間的關(guān)系曲線如圖6所示。

將該曲線與解析理論相結(jié)合，應(yīng)用第1節(jié)提出的隧洞解析應(yīng)力空間分布法可獲得引漢濟渭工程隧洞斷面空間應(yīng)力分布云圖如圖7所示。

如圖7所示，隧洞掌子面推進過程中，隧洞左拱腳及右拱肩處環(huán)向應(yīng)力集中較大，尤其是左拱腳，最大可達5.0倍的初始主應(yīng)力場豎向應(yīng)力分量，而右拱腳及左拱肩處環(huán)向應(yīng)力集中相對較小。沿隧洞軸線縱向，距離掌子面0.5D范圍內(nèi)，隧洞左拱腳處環(huán)向應(yīng)力梯度變化較大，右拱肩處變化其次，而左拱肩處變化相對較小。距離掌子面0.5～2.0D范圍內(nèi)，隧洞左拱腳及右拱肩處環(huán)向應(yīng)力集中區(qū)域隨掌子面的推進在逐漸增大。距離掌子面2.0D范圍外，隧洞環(huán)向應(yīng)力分布區(qū)域基本趨于穩(wěn)定。由此可見，沿隧洞軸縱向，逐漸遠離掌子面時，隧洞環(huán)向集中應(yīng)力分布呈現(xiàn)出3階段特征，應(yīng)力集中程度快速增長區(qū)、應(yīng)力集中范圍擴展區(qū)、應(yīng)力集中穩(wěn)定區(qū)。

圖6 隧洞頂部位移釋放系數(shù)與掌子面之間的關(guān)系Fig.6 Relationship between displacement release coefficient of tunnel top and tunnel working face

D為隧洞斷面最大洞徑圖7 隧洞空間應(yīng)力分布云圖Fig.7 Spatial stress distribution of tunnel

4 初始主應(yīng)力場水平傾角的影響

為了探究隧洞初始主應(yīng)力場水平傾角α對隧洞解析應(yīng)力空間分布特征的影響，設(shè)置了如表2所示的3種應(yīng)力邊界工況。

基于3.1節(jié)隧洞斷面形狀及力學參數(shù)，通過數(shù)值模擬獲得了不同主應(yīng)力場水平傾角條件下的隧洞頂部位移釋放系數(shù)與掌子面不同距離之間關(guān)系曲線如圖8所示，相應(yīng)曲線方程如表3所示。

將表3中擬合曲線方程分別與隧洞解析理論相結(jié)合，通過應(yīng)用第1節(jié)提出的隧洞解析應(yīng)力空間分布方法可獲得不同初始主應(yīng)力場條件下隧洞圍巖解析應(yīng)力空間分布云圖如圖9所示。

表2 模型邊界不同應(yīng)力場參數(shù)Table 2 The different stress field parameters of the model boundary

圖8 隧洞頂部位移釋放系數(shù)與掌子面之間的關(guān)系Fig.8 Relationship between displacement release coefficient of tunnel top and tunnel working face

表3 擬合曲線方程Table 3 Fitting curve equations

圖9 不同初始主應(yīng)力場水平傾角下隧洞環(huán)向應(yīng)力集中空間分布云圖Fig.9 Spatial distribution nephogram of circumferential stress concentration of tunnel under different horizontal inclination angles of initial principal stress field

如圖9所示，在α=0° 時，隧洞空間應(yīng)力集中區(qū)域主要分布在隧洞左右拱腳及拱肩處；在α=30° 時，隧洞空間應(yīng)力集中區(qū)域分布主要集中隧洞左拱腳及右拱肩處，而隧洞左拱肩及右拱腳處環(huán)向應(yīng)力集中程度與α=0° 時相比，在降低；在α=60° 時，隧洞空間應(yīng)力集中區(qū)域分布主要表現(xiàn)在左拱腳及拱頂偏右側(cè)，與α=30° 時相比，隧洞頂部偏右側(cè)環(huán)向應(yīng)力集中程度在升高。由此可見，初始主應(yīng)力場水平傾角變化對隧洞環(huán)向應(yīng)力環(huán)向空間分布位置影響較大。沿隧洞軸線縱向逐漸遠離掌子面時，不同初始主應(yīng)力場水平傾角下隧洞環(huán)向應(yīng)力集中區(qū)域也呈現(xiàn)出應(yīng)力集中梯度快速升高、應(yīng)力集中區(qū)域逐漸擴展、應(yīng)力集中區(qū)域穩(wěn)定三階段特征，說明初始主應(yīng)力場水平傾角變化對隧洞環(huán)向應(yīng)力集中縱向變化影響較小。

5 結(jié)論

(1)基于收斂約束法，結(jié)合非圓形隧洞圍巖解析解與可視化軟件Tecplot，提出了一種隧洞解析應(yīng)力空間分布方法，通過該方法可以快速直觀地獲取隧洞空間應(yīng)力分布圖，可為實際隧洞支護設(shè)計及施工方案制定提供借鑒。

(2)隧洞開挖過程中，沿隧洞斷面縱軸逐漸遠離掌子面時，隧洞環(huán)向應(yīng)力集中區(qū)域分別呈現(xiàn)應(yīng)力集中梯度增長區(qū)、應(yīng)力集中區(qū)域擴展區(qū)、應(yīng)力集中區(qū)域穩(wěn)定區(qū)三階段特征，對于硬巖隧洞，三階段分布區(qū)域分別為0.5D范圍之內(nèi)，0.5D～2.0D范圍之間，2.0D范圍之外。

(3)初始主應(yīng)力場水平傾角變化對隧洞應(yīng)力集中區(qū)域環(huán)向分布方位影響較大，而對其沿隧洞軸縱向分布影響較小。