基于復變函數(shù)方法的水下隧道圍巖彈性分析

蔚立元，陳曉鵬，韓立軍，王迎超

(1.中國礦業(yè)大學 深部巖土力學與地下工程國家重點實驗室，江蘇 徐州 221008；2.中國礦業(yè)大學 力學與建筑工程學院，江蘇 徐州 221116；3.中國水電顧問集團 成都勘測設計研究院，成都 610072)

1 引 言

近現(xiàn)代以來交通需求大大增加，江河湖海等水體成了社會進一步發(fā)展的制約和瓶頸，人類便利用輪渡、橋梁和隧道等手段進行跨越。最近幾十年來，水下隧道由于其獨有的特點得到了迅猛發(fā)展，逐漸成為與橋梁并舉的跨越方式[1－2]。

一般而言，隧道圍巖是復雜多變的地質體，但是，通過固結灌漿等人工改良，使得其中的結構面閉合，巖土體完整性有保障后，圍巖可視為均質、連續(xù)、各向同性介質[3]。覆蓋層厚度是水下隧道設計施工中的重要參數(shù)和關鍵指標，許多學者進行了相關研究[4－6]，一般認為其合理值不超過2倍洞徑，應屬于有地面荷載(上覆水體)作用的淺埋隧道，這樣水下隧道圍巖應力場、位移場的求解就可以簡化為含圓孔半平面體在水平邊界上受任意分布荷載的彈性力學問題，并可用解析方法求解。

與深埋隧道相比，淺埋隧道由于受到地表邊界和地面荷載的影響，其力學分析在數(shù)學處理上歷來存在較大的困難[7]。Jeffery[8]和Mindlin[9－10]曾采用雙極坐標法進行求解，地面載荷對于應力場的影響在他們的解答中得到了考慮，但是這種方法沒有給出位移場分布，而且不適用于隧道非常接近地表的情況。Bobet[11]和Chow[12]在諸多嚴格假定條件下，基于彈性力學的連續(xù)和位移諧調條件，研究了飽和圍巖中襯砌隧道的應力場和位移場。Verruijt等[13－16]提出了合理的映射函數(shù)，采用復變函數(shù)方法，分別研究了第二類邊界條件、不考慮巖土體自重情形下，無襯砌淺埋隧道變形導致的地表變形問題，以及考慮巖土體自重情況下，由于隧道構造浮力效應所導致的地表變形問題。王立忠等[17]應用Verruijt提供的基本解法，求得隧道周邊給定位移條件下的位移場和應力場，分析了不同埋深、泊松比對位移的影響以及不同埋深對應力場的影響。

上述淺埋隧道的復變函數(shù)方法研究中，都是在洞周位移或應力邊界條件下進行分析的，由水下隧道抽象而來的“含圓孔半平面體在水平邊界上受任意分布荷載”問題尚無文獻報導。本文將借助Verruijt提供的共形映射函數(shù)，采用復變函數(shù)方法對該問題進行求解。

2 問題的描述

水下隧道，或推廣為地面荷載下淺埋隧道的圍巖分析，可視為任意分布荷載q(z)作用下，邊界附近含圓形孔洞的半平面體的彈性分析，在直角坐標系(z平面)下的構形如圖1所示。圖中，圓形隧道半徑為 r0；圓心與地表相距 h；覆蓋層厚度為 d；半徑與圓心深度之比r0/h為一個重要的幾何參數(shù)，圍巖 區(qū)域(即 x 軸以下的圓外區(qū)域)標記為R。

根據(jù)Muskhelishvili的復變函數(shù)方法[18]，應力和位移分量(平面應變條件)可用在區(qū)域R內處處解析的復位勢φ1(z)、ψ1(z)來表示，見式（1）。

式中：σx、σy、τxy、ux、uy分別為平面問題的應力分量和位移分量為剪切模量；G為剪切模量；γ為泊松比。借助應力邊界條件的復變函數(shù)式（2），本問題的兩個應力邊界可表示為式（5）、（6）。

3 共形映射

如圖1所示水下隧道的幾何構形可以概化為含圓孔半平面體。通過保角映射可把物理平面 z 內的區(qū)域 R 變換為像平面ζ內的圓環(huán)域[13]。

分式線性映射具有保角性、保形性、保圓性(直線看成是圓的特例)和保對稱性，定義式為

式中：a、b、c、d為不全為 0 的待定常數(shù)。3個獨立條件便能確定一個分式線性映射[19]。如果要求將z 平面上的點 A(0，0)，B(∞，0)，C(0，-d)映射為ζ平面上的對應點 A’(-1，0)，B’(1，0)，C’(-α，0)，則得共形映射函數(shù)見式（6），其中，α為由r0/h 確定的參數(shù)，見式（7）。z 平面上的R區(qū)域映射為ζ平面的單位圓環(huán)域r，見圖2。容易驗證，圓|ζ|=1對應于直線 y=0，而圓|ζ|= α 對應于圓 x2+(y+h)2=若α→0，則r0/h→0，即隧道半徑較小，對應于深埋隧道；若α→1，則隧道半徑接近于圓心埋深，對應于淺埋隧道。

由于復位勢φ1(z)、ψ1(z)在區(qū)域R內解析，共形映射函數(shù)ω(ζ)在圓環(huán)域r內解析，由復變函數(shù)理論可知ζ平面上的復位勢φ(ζ)、ψ(ζ)在該圓環(huán)域上也處處解析，這意味著它們可用Laurent級數(shù)展開：

式中：系數(shù)ak、bk、ck、dk由邊界條件確定。

圖2 保角映射后的區(qū)域Fig.2 The region after conformal mapping

4 邊界條件

本問題是圓周邊界，其表達式比較簡單，而對于復雜的孔洞問題，該因子可能異常復雜以致問題難以求解。

4.1 地表邊界條件

地表受分布荷載q(z)作用，借助式（3），像平面ζ上的地表邊界條件表示為

當?shù)乇矸植己奢d的函數(shù)形式確定后，Q(ζ)一般可展開成Fourier級數(shù)。水下隧道地表邊界y=0，即z=x上作用著均布荷載q(z)=qx+i qy，當作用范圍為[l1，l2]時，式（3）右端的積分可表示為

在像平面ζ中，該邊界ρ=1，故ζ=ρσ=σ=exp(i)θ，式（12）轉換為

可見，Q(ζ)為周期為2π的周期函數(shù)，且滿足Dirichlet收斂條件，故可展開成Fourier級數(shù)：

將式（8）、（11）、（14）代入式（10），然后，令等號兩邊σ的同次冪相等，整理可得式（16）。可見，系數(shù)ck、dk可用ak、bk顯式表示。現(xiàn)在問題的關鍵在于系數(shù)ak、bk的確定，而它們可以通過洞周邊界條件求得。

4.2 洞周邊界條件

由式（4），洞周的自由邊界在ζ平面上表示為

將式（8）、（18）代入式（17），得到一個關于系數(shù)ak、bk、ck、dk的相當復雜的方程組，再利用式（16），就能整理出系數(shù)ak、bk的迭代方程組：

如果已知ak，bk，通過求解一個二元一次方程組，就能得到系數(shù)ak+1，bk+1。初始解a1，b1可以借助對σ0、σ1項的整理得到：

式（20）中的系數(shù)a0不產生應力，對應于剛體位移，可視為0，這樣可求得

5 返回物理平面求解

為便于表示，可令

然后，由式（6）可得

對于物理平面中的任一點z=x+iy，都可以由式（24）求出對應的像平面中的點ζ，將其代入式（6）、（8），再利用像平面上的應力、位移表達式（25），就能得到該點的應力和位移分量。

6 程序流程和算例

6.1 程序流程

按照上述方法，能夠計算上覆水體荷載作用下水下隧道圍巖的應力和位移，但計算過程相當繁瑣。通過編寫程序，借助計算機求解是一個方便的途徑。編寫fortran程序的流程見圖3（判斷ak是否趨近于非0極限時，推薦從l=1 000開始，每增加1 000判斷一次）。

6.2 水下隧道算例

對于水下隧道而言，由于上覆水體的黏性對巖土體表面的影響可略去不計[20]，故荷載為其產生的均布法向荷載，且作用范圍非常廣泛，但在實際運算中選擇 10倍洞跨左右，就可以滿足工程精度要求。在本例中，隧道半徑為5 m，覆蓋層厚20 m，上覆水體深30 m，具體參數(shù)可見表1。

6.3 正確性驗證

圖3 計算程序流程圖Fig.3 The flow chart of the calculation program

表1 水下圓形隧道計算參數(shù)Table 1 Calculation preferences for an underwater tunnel

為了驗證本方法的正確性，應用權威有限元分析軟件Ansys對上述案例進行分析。由于近似解析方法的求解域是含圓孔半平面體，數(shù)值計算只能選擇盡量大的區(qū)域(210 m×100 m，見圖4)來逼近，計算采用平面應變模式。人工邊界條件如下：上表面是應力邊界，而側邊和下表面為法向約束的位移邊界條件。其他條件和參數(shù)與章節(jié)6.2中相同。

圖4 水下隧道計算模型Fig.4 Calculation model for the underwater tunnel

由計算得出得兩個應力分量，利用式（25）可求出主應力。

覆蓋層中和洞周邊界上的6個關鍵點的兩種方法計算結果的對比見表 2，一對數(shù)據(jù)中前面的是數(shù)值計算結果，后面的是近似解析方法結果，其中負值表示壓應力，ux、uy分別為水平位移和豎直位移。由表可見，近似解析法和有限元方法的結果基本相同，說明本方法是可靠的。

表2 關鍵點計算結果對比Table 2 The calculation results comparison of key points

6.4 應力分析

主應力1σ的極值分別是拱頂附近的-8.14 kPa和右拱腳附近的-0.20 MPa，主應力2σ的極值分別是拱頂附近的-0.08 MPa和左右邊墻中部的-0.80 MPa。水平方向隧道中線左右70 m，垂直方向地表以下80 m范圍內的圍巖主應力等值線圖如圖5所示。由圖可見，中線左右2倍洞跨、底板以下兩倍洞跨到地表范圍內圍巖應力集中現(xiàn)象明顯，覆蓋層尤其應給予重點關注。

6.5 位移結果

由于問題的幾何形狀、約束情況、以及所受荷載在水平方向上的對稱性，隧道中線的水平位移ux為 0，圍巖的水平位移都可以精確計算。然而，半平面體垂直方向上不受約束，所以垂直位移（即沉陷）uy不能確定，這時只能求得相對沉陷。取隧道中線上的點 J(0，-100)為基點，圍巖中任一點對于基點 J的相對沉陷等于該點的沉陷減去 J點的沉陷。

水下隧道覆蓋層中不同深度位置的沉降槽如圖6所示，圖7為拱頂以上中線各點沉陷值隨深度的變化曲線，變形前后的隧道輪廓對比見圖8(為便于表達，洞周點位移是實際值的1 000倍)，圍巖位移的等值線如圖9所示。

綜合各圖，水下隧道圍巖的位移有以下特點：(1) 隧道中線兩側 2倍洞跨范圍內，沉降槽比較陡峭，然后逐漸平滑，最后趨近于水平。(2) 從深度方向來看，地層沉陷結果顯示出明顯的分層傳遞性，即沉降變形由地表到拱頂傳遞，變化逐漸減弱，沉降值由拱頂?shù)降乇磉f增，在近原點處取得最大值2.37mm，但沉降槽隨著深度增加越來越陡峭。(3)隧道的水平收斂非常小，可以忽略不計，洞周點垂直位移隨距地表距離增加而逐漸減小，從 2.09 mm(拱頂)依次變?yōu)?1.72 mm(邊墻中部)，1.36 mm(底板)，隧道斷面面積縮小為原面積的 91.7%。(4)中線兩側2倍洞跨、底板以下一倍洞跨到地表范圍內，圍巖位移受隧道開挖影響明顯，邊墻外側巖土體水平位移由向洞內收斂急劇變成向遠方移動，設計施工中應特別注意。

圖5 圍巖主應力等值線圖（單位：102 kPa）Fig.5 Isoline maps of principal stresses(unit: 102 kPa)

圖8 隧道變形后輪廓Fig.8 Contour of a underwater tunnle after deformation

圖9 圍巖位移分量等值線圖Fig.9 Isoline maps of displacement components

7 結 語

把水下隧道圍巖的彈性分析從力學角度理想化為“含圓孔半平面體在地表邊界上受任意分布荷載”問題，利用分式線性映射的特點，確定含圓孔半平面體變換為像平面上單位圓環(huán)域的共形映射函數(shù)。基于 Verruijt提供的基本方法，推導該問題的一般解法，給出均布荷載邊界條件在像平面上的級數(shù)形式。采用Fortain語言，編寫了計算程序，給出特定水下隧道算例的圍巖應力、位移結果，并分析了其受力變形特點，并指出設計施工中應該注意的地方。

本文從理論上完善了半無限平面含括一個孔洞在地表分布荷載邊界條件下的求解過程，對水下隧道的具體設計施工具有一定的指導意義。

[1] 王夢恕.水下交通隧道發(fā)展現(xiàn)狀與技術難題[J].巖石力學與工程學報, 2008, 27(11): 2161－2172.WANG Meng-shu.Current developments and technical issues of underwater traffic tunnel[J].Chinese Journal of Rock Mechanics and Engineering, 2008, 27(11): 2161－2172.

[2] 王夢恕.臺灣海峽海底鐵路隧道建設方案[J].隧道建設, 2008, 28(5): 517－526.WANG Meng-shu.Construction scheme of Taiwan Strait subsea railway tunnel[J].Tunnel Construction, 2008,28(5): 517－526.

[3] 薛守義, 劉漢東.巖體工程學科性質透視[M].鄭州:黃河水利出版社, 2002.

[4] 李術才, 李樹忱, 徐幫樹, 等.海底隧道最小巖石覆蓋厚度確定方法研究[J].巖石力學與工程學報, 2007,26(11): 2289－2295.LI Shu-cai, LI Shu-chen, XU Bang-shu, et al.Study on determination method for minimum rock cover of subsea tunnel[J].Chinese Journal of Rock Mechanics and Engineering, 2007, 26(11): 2289－2295.

[5] 張頂立, 李兵, 房倩, 等.基于風險系數(shù)的海底隧道縱斷面確定方法[J].巖石力學與工程學報, 2009, 28(1): 9－19.ZHANG Ding-li, LI Bing, FANG Qian, et al.Method for determining longitudinal section of subsea tunnel based on risk coefficient[J].Chinese Journal of Rock Mechanics and Engineering, 2009, 28(1): 9－19.

[6] 蔚立元, 徐幫樹, 李術才, 等.確定水下隧道覆蓋層厚度的經(jīng)驗公式及其應用[J].地下空間與工程學報, 2011,7(3): 497－503.YU Li-yuan, XU Bang-shu, LI Shu-cai, et al.Empirical formulas for deciding overburden thickness of underwater tunnel and relevant applications[J].Chinese Journal ofUnderground Space and Engineering, 2011, 7(3): 497－503.

[7] SAGASETA C.Analysis of undrained soil deformation due to ground loss[J].Geotechnique 1987, 37(4): 301－326.

[8] JEFFERY G B.Plane stress and plane strain in bipolar coordinates[J].Transactions of the Royal Society of London,(Series A), 1920, 221: 265－293.

[9] MINDLIN R D.Stress distribution around a tunnel[J].Transactions of the ASCE, 1940, 1117－1153.

[10] MINDLIN R D.Stress distribution around a hole near the edge of a plate under tension[C]//Proceedings of the Society of Experimental Stress Analysis, [S.l.]: [s.n.],1948, 5: 56－57.

[11] BOBET A.Analytical solutions for shallow tunnels in saturated ground[J].ASCE Journal Engineering Mechanics, 2001, 127: 1258－1266.

[12] CHOU W I, BOBET A.Predictions of ground deformations in shallow tunnels in clay[J].Tunnelling and Underground Space Technology 2002; 17: 3－19.

[13] VERRUIJT A, BOOKER J R.Surface settlements due to deformation of a tunnel in an elastic half plane[J].Geotechnique, 1996, 46: 753－756.

[14] VERRUIJT A.Deformations of an elastic half plane with a circular cavity[J].International Journal of Solids and Structures 1998; 35(21): 2795－2804.

[15] VERRUIJT A.A complex variable solution for a deforming circular tunnel in an elastic half-plane[J].International Journal for Numerical and Analytical Methods in Geomechanics 1997, 21: 77－89.

[16] STRACK O E, VERRUIJT A.A complex variable solution for a deforming buoyant tunnel in a heavy elastic half-plane[J].International Journal for Numerical and Analytical Methods in Geomechanics, 2002; 26: 1235－1252.

[17] 王立忠, 呂學金.復變函數(shù)分析盾構隧道施工引起的地基變形[J].巖土工程學報, 2007, 29(3): 319－327.WANG Li-zhong, Lü Xue-jin.A complex variable solution for different kinds of oval deformation around circular runnel in an elastic half plane[J].Chinese Journal of Geotechnical Engineering, 2007, 29(3): 319－327.

[18] MUSKHELISHVILI N I.Mathematical theroy of elasticity[M].Leyden: International Publishing, 1954.

[19] 劉建亞.復變函數(shù)與積分變換[M].北京: 高等教育出版社, 2005.

[20] 朱鏡清, 周建.海底地震動估計的一個流體力學基礎[J]地震工程與工程振動, 1991, 11(3): 87－93.ZHU Jing-qing, ZHOU Jian.A fluid mechanics basis for estimating undersea ground motion[J] Earthquake Engineering and Engineering Vibration, 1991, 11(3):87－93.