張天宇
(大連海事大學,遼寧 大連 116026)
近年來,隨著我國軌道交通事業的快速發展,有關隧道滲流場的研究越來越受到關注。隧道滲流場會與隧道圍巖應力場和隧道圍巖溫度場發生耦合,從而影響隧道襯砌內部應力[1]。當滲流與應力單獨耦合時,甚至還會發生涌水等重大地質災害[2]。隧道滲流還會對軟土隧道產生影響,可能會使隧道地表沉降以及地層部分土壤組成物損失[3]。在隧道工程中,滲流分析受土壤滲透性質、隧道襯砌防滲性質[4]以及地表水頭取值等一系列因素影響,其中隧道襯砌滲漏對隧道滲流場的影響較大,所以部分襯砌滲漏情況下的隧道滲流場具有重要研究意義。
目前,對于隧道滲流場的研究方法主要有數值法、模型試驗法和解析法。科研人員利用數值法和模型試驗,已經對隧道滲流場有了較為深入的研究[5-8]。
而解析法則由嚴謹的公式導出,可使分析結果更直觀精確,有利于結果的對比與評估,便于深入研究。國內外學者對隧道滲流場在不同的工況條件下,用不同的方法進行了解析研究。鏡像法對于解決具有對稱性的滲流解析問題具有重要意義。李鵬飛等[9]利用疊加原理以及綜合鏡像法,得出了考慮襯砌時平行四孔隧道滲流場的解析解;張冶國等[10]把鏡像法與復變函數理論相結合,得出了鄰近補水斷層的水頭分布解析式。相較于鏡像法而言,保角變換可以解決復雜邊界條件下的滲流問題。杜朝偉等[11]根據無線含水層豎井理論以及保角變換法,得出了隧道的襯砌水壓力計算公式;ZHANG等[12]使用保角變換法,得出了在非達西滲流條件下的帶襯砌隧道滲流場解析解;PARK等[13]以保角變換法為基礎,把半無限地層映射為圓環,從而得出隧道滲流場解析解在不同邊界滲流條件下的情況。
此外,還有諸多解析方法,趙建平等[14]采用雙極坐標來描述等勢線,得出了富水地區隧道滲流場的解析解。謝寒松等[15]把地下水水力學理論以及復變函數反演變換作為基礎,得出了充水巖溶隧道的滲流場解析解。然而,這些解析研究大多將隧道襯砌模型理想化,很少考慮隧道襯砌特定角度滲漏水對淺埋隧道滲流場的影響。
針對上述研究現狀,本文對淺埋圓形隧道部分襯砌滲漏條件下的滲流場進行了解析研究。本文應用保角變換,把非齊次邊界條件下的二維穩態達西滲流方程轉換到圓環坐標系下,在MATLAB軟件中得到圓環坐標系下的水頭空間函數后再轉換回原直角坐標系,從而得出隧道襯砌滲漏條件下的隧道滲流場解析解,與ANSYS中進行的仿真模擬結果吻合較好。通過對本文解析解進行詳細的參數分析,發現解析解結果受隧道半徑、隧道埋深、地表常水頭以及隧道襯砌滲漏部位的圓心角度等因素影響。本文參數分析的結果對實際隧道襯砌滲漏問題的分析有一定意義。
對于大多數三維隧道的達西滲流問題,由于其長度尺寸遠大于平面尺寸,故近似按照平面應變[16]問題來求解其橫截面的穩態滲流問題,建立的二維淺埋圓形隧道橫截面滲流模型見圖1。

圖1中,H為地表常水頭的取值。隧道中心到地表的距離為h,隧道半徑為r,C點是對應隧道極角為?處的襯砌單元,隧道襯砌位于方位角δ處的對應圓心角為α處的劣弧AB漏水,其余角度處襯砌完好。
在模型土體(x,y)處,二維穩態達西滲流方程為:
(1)
其中,φ(x,y)為位置坐標(x,y)處土體單元水頭。
在地表處有邊界條件:
φ(x,y)∣y=0=H
(2)
在隧道襯砌滲漏水的AB劣弧部分:

(3)
在隧道襯砌不漏水的AB優弧部分:
φ(x,y)∣x2+(y+h)2=r2=-h+rsinθ
(4)
在以下計算中采取如下假定:不考慮地表水頭隨著隧道襯砌滲漏水過程的變化;不考慮土層各高度處滲流性質的差異;不考慮隧道的漏水襯砌部分與不漏水襯砌部分連接處的突變性;不考慮隧道襯砌與滲漏水體接觸后滲漏性質的變化;隧道滲漏水部分的水頭等于其位置水頭;隧道非漏水部分在其表面的水頭函數法向梯度為0。
將圖1中的直角坐標系通過如下保角變換[17]式轉換到圖2的圓環坐標系中:
(5)


圖2中,除去圖1中已有標識外,A′,B′和C′點分別是原直角坐標系中A,B和C點保角變換后在圓環坐標系中對應的點,而原直角坐標系中?,δ以及角α則分別對應保角變換后圓環坐標系中的θ,Ω以及ε角。
由坐標轉換式(2),得到方程(1)在圓環坐標系中的表達形式為:
(6)
由數學物理方法[18]中的分離變量法,可得方程(3)的通解為:

(7)
相應的邊界條件式(2)~式(4)變為:
(8)
將式(8)中的第一式代入式(7)后,由圓環坐標系中達西滲流方程的對稱性可得:
(9)
采用邊界配點法[19]求解圓環坐標系下的方程(6)。對于隧道襯砌漏水的AB劣弧以及不漏水的AB優弧部分,沿隧道襯砌逆時針的方向選取M個控制點,其中有m個控制點落在AB劣弧部分,有(M-m)個控制點落在優弧部分。
對落在AB劣弧部分的m個控制點,將式(8)中第二式代入式(7)后得:
(10)
對落在AB優弧部分的(M-m)個控制點,將式(8)中第三式代入式(7)后得:

(11)
將式(7)的無窮級數項數取前N項來進行后續計算,聯立式(10),式(11),可得廣義矩陣方程如下所示:
[D]∣(M,2N+1)[Φ]∣(2N+1,1)=[H]∣(M,1)
(12)
式(12)中級數的前N項系數矩陣[Φ]∣(2N+1,1)、相應的左乘矩陣[D]∣(M,2N+1)以及水頭矩陣[H]∣(M,1)分別滿足:
[Φ]∣(2N+1,1)=[B0,A1,B1,…,AN,BN]T
(13)

(14)
[H]∣(M,1)=[0,…,0,hr(θm+1)-
H,…,hr(θM)-H]
(15)
在計算過程當中選取隧道襯砌計算點總數M遠大于所取級數項系數2N+1,以保證級數系數求解的精度。由式(13)可得級數的系數矩陣[Φ]∣(2N+1,1)滿足:
[Φ]∣(2N+1,1)=[D]-1∣(M,2N+1)[H]∣(M,1)
(16)
將式(16)中相應的級數系數代入式(7)可以得到圓環坐標系下的穩態滲流方程解。通過式(5)中的保角變換便可得地面坐標系中的穩態滲流方程通解。
為驗證本文解析解的正確性,令隧道中心到地表的距離為h=20 m,隧道半徑r=4 m。選取隧道襯砌控制點個數M=3 600,并且沿著隧道襯砌逆時針均勻分布。隧道襯砌漏水部分對應的圓心角α=20°,模型地表水頭高度H=0 m。經過試算,當式(6)中級數項數取N=17時,解析解結果收斂,并且與下文ANSYS有限元軟件的模擬結果吻合良好。
在ANSYS有限元軟件模擬的過程當中,選取平面模型矩形尺寸為50 m×50 m,矩形除上方邊界設置為常水頭外H=0 m,其余三面均設置為不透水及水頭函數的導數在其余三面為0。
選取隧道滲漏水的襯砌中心分別位于圖1中δ=80°,35°,-10°,-55°處,分別繪制出相應襯砌滲漏水位置的水頭函數與ANSYS有限元軟件的計算結果對比如圖3所示。

由圖3可知,4種角度的隧道滲漏解析解水頭與ANSYS有限元軟件模擬的相應隧道襯砌工況的結果在遠離隧道襯砌的部分十分符合,在隧道襯砌漏水部分的附近有些許差距,但也吻合較好。隨著隧道襯砌漏水部分距離地表越深,其附近的水頭等值線數值越小,這是符合初始條件式(11)的物理規律的。
取重力加速度為g=10 m/s2,水體密度為ρ0=103kg/m3,忽略水體滲流時的動能,以大氣壓強為計算零點,由伯努利方程,可得水頭為φ0處的土體滲流水壓力p(φ0)滿足以下方程:
p(φ0)=ρ0gH-ρ0gφ
(17)
聯立式(7)和式(17)可得隧道橫截面穩態滲流時任一點水壓力。
選取與圖3中相同的角度δ=80°,35°,-10°,-55°使襯砌部分滲漏水,分別繪制出相應襯砌滲漏水位置的水壓力函數與ANSYS有限元軟件的計算結果對比如圖4所示。

由圖4可知,4種角度的隧道滲漏解析解水壓力與ANSYS有限元軟件模擬的相應隧道襯砌工況的結果在遠離隧道襯砌的部分十分符合,在隧道襯砌漏水部分的附近有些許差距,但也吻合較好。隨著隧道襯砌漏水部分距離地表越深,其附近的水壓力等值線數值越大,這是符合式(17)的物理規律的。
為進一步驗證本文的解析解對于隧道襯砌滲漏部分的敏感性,現選取隧道襯砌上每隔120°的處于正三角形三個頂點上的3個滲漏點,每個滲漏點對應的圓心角仍為α=20°。所繪制出的解析解水頭與水壓力等值線與ANSYS有限元軟件對比如圖5所示。

由圖5可知,在正三角形頂點處均滲漏水的情況下,本文水頭與水壓力的解析解均與ANSYS有限元軟件模擬的相應隧道襯砌工況的結果較為吻合。隨著隧道襯砌漏水部分距離地表越深,其附近的水頭等值線數值越小,水壓力等值線數值越大,并且仍然呈現出左右對稱的趨勢。由此可見,本文無論是水頭還是水壓力的解析結果對于隧道襯砌多處漏水的情況有較高的靈敏度,這也更好地驗證了本文解析解的正確性與可靠性。
為對本文解析解進行進一步的參數分析,采用控制變量法,分別改變模型條件中的隧道半徑、隧道埋深、地表水頭以及隧道襯砌滲漏水的角度,在此過程當中其他參數不變,選取隧道拱頂中心滲漏的情況分析如下。
選取隧道半徑分別為r=3 m,4 m,5 m,保持其余變量不變,繪制出3種隧道半徑條件下的隧道橫截面水頭等值線如圖6所示。
分析圖6可知,隨著隧道半徑的增加,相同數值的水頭線逐漸向隧道外圍正法線方向移動。在半徑增幅相同的情況下,半徑基值越大,水頭等值線向外移動的幅度越小,這為同一地質條件下的變截面隧道滲流分析提供了參考。
選取隧道中心到地表的距離分別為h=19 m,20 m,21 m,保持其余變量不變,繪制出三種隧道埋深條件下的隧道橫截面水頭等值線如圖7所示。

分析圖7可知,隨著隧道埋深的增加,在隧道拱頂一定范圍內,相同數值的水頭線逐漸向隧道外圍負法線方向移動。隨著與隧道拱頂豎直直線夾角的增加,相同數值的水頭線逐漸變化至相交,最后向隧道外圍正法線方向移動。在半徑增幅相同的情況下,半徑基值越大,水頭等值線向外移動的幅度越小,為不同海拔的隧道橫截面滲流分析提供了基礎。
選取地表水頭分別為H=0 m,1 m,2 m,保持其余變量不變,繪制出三種地表水頭條件下的隧道橫截面水頭等值線如圖8所示。
分析圖8可知,隨著地表水頭的增加,相同數值的水頭線向隧道外圍負法線方向大幅移動。這為縱向處于非水平地表水流的隧道橫截面滲流分析提供了基礎。
選取隧道拱頂中心滲漏的角度分別為α=20°,25°,30°,保持其余變量不變,繪制出3種隧道拱頂中心滲漏角度條件下的隧道橫截面水頭等值線如圖9所示。

分析圖9可知,隨著隧道拱頂中心滲漏角度的增加,相同數值的水頭線向隧道外圍正法線方向移動,在隧道拱頂正上方,3種隧道拱頂滲漏角度相同位置的水頭差別不大。隨著與隧道拱頂豎直直線夾角的增加,不同拱頂滲漏角度的水頭等值線差別逐漸增大,為同一隧道不同拱頂滲漏情況下的橫截面滲流分析提供了基礎。
針對淺埋圓形隧道部分襯砌滲漏條件下的穩態滲流場,本文采用保角變換,結合極坐標條件下拉普拉斯方程的級數解,利用MATLAB的廣義矩陣運算求解出級數解中相應的系數,得出了隧道襯砌滲漏條件下的滲流解析解,并得出了以下結論:
1)在不同隧道襯砌滲漏條件下,本文用MATLAB軟件實現的穩態滲流解析解的水頭函數和水壓力函數均與ANSYS有限元軟件的模擬結果較為吻合,解析精度較高。
2)本文的解析解相對于如ABAQUS,PALXIS以及ANSYS之類的有限元軟件而言,具有更高的運行效率,在普通辦公電腦上用MATLAB得出解析解計算結果只需4 s左右。并且只要輸入相應的模型計算參數,即可得出隧道穩態滲流水頭和水壓力函數在任意處的具體取值,相比于有限元軟件從建模到得出結果的一系列過程來說更為簡潔、通用。
3)對本文模型的隧道半徑、隧道埋深、地表常水頭取值以及隧道襯砌滲漏部分對圓心角進行了一系列詳細的參數分析,討論了在各自變量變化時相應的滲流水頭函數的變化趨勢,并且給出了各個參數分析的結果對相應工況條件下隧道滲流場分析的意義。
4)將本文模型的土體條件進行更改,便可以在保角變換作用下求出在多層土體中的隧道滲流場解析解。因此,本文模型的適用范圍較為廣泛。