淺埋海底隧道圍巖應力復勢函數顯式解1)

金 波 田俊彤 方棋洪

* (湖南大學土木工程學院,工程結構損傷診斷湖南省重點實驗室,長沙 410082)

? (湖南湖大建設監理有限公司,長沙 410082)

** (湖南大學機械與運載工程學院,長沙 410082)

引言

國內經濟持續發展、城市規模不斷擴大、城市之間聯系愈發緊密,客流量與貨物運輸量急劇增加,為了緩解城市交通壓力、節約土地資源和開發地下空間,建設一批市政管道和跨水域隧道成為最優方案.特別是在沿海城市,海底隧道具有抗震性能好、不影響貨船航運、隱蔽性高和不受氣候影響等優勢而被廣泛采用.這些隧道部分區段往往埋深較淺,建設中不可忽視地下水壓、巖體自重以及地表邊界上其他復雜工況的影響.準確得出隧道周邊應力分布與變形對指導隧道設計以及安全施工具有重要的理論意義和工程價值.

20 世紀初,俄國力學家Kolosov[1]采用復應力函數法解決二維彈性靜力學問題,Muskhelishvili[2]系統論述了彈性力學平面問題的復變函數法,極大促進該方法的應用與發展.保角映射能將復雜區域轉化為簡單規則區域,使得復變函數在研究平面含孔洞問題時具有明顯優勢,因此,平面問題復變函數方法在巖土工程領域得到了廣泛應用.

深埋單孔隧道屬于無限單連通域問題,研究中可忽略巖體重力梯度的影響,不同孔型的深埋隧道其映射函數均可采用Laurent 級數表示,對深埋隧道的理論研究已趨于成熟[3-17].采用復變函數方法研究淺埋隧道問題起步較晚,1997 年Verruijt[9-10]提出分式映射函數,將含圓孔的半平面巖石物理區域保角映射為像平面上同心圓環域,獲得了在不考慮巖體重力影響,僅有孔口均布徑向位移與均布徑向面力邊界條件時的解答;Verruijt 等[11]考慮周圍巖體自重應力對開挖淺埋單孔隧道的影響,將未被開挖的初始狀態、考慮巖體作用時開挖巖石等效合力狀態和孔口附加等效載荷3 種狀態的解疊加成為最終解;Verruijt 等[12]在原有復應力函數中增加了能夠反映彈性體在不平衡力系作用下的影響項,給出了二維平面問題在不平衡力系作用下復應力函數的一般形式,并對孔口位移邊值條件進行了研究.蔚立元等[13]、王志良等[14]和陸文超等[15]得到地面任意分布載荷下的圍巖應力場;韓凱航等[16]基于復變函數法并結合柯西-黎曼方程得出淺埋隧道應力及位移函數的級數顯式表達式;宋文杰等[17]求解了不排水黏土淺埋盾構隧道地層變位和地層應力;宋浩然等[18-19]與Fang 等[20]研究圍巖自重及水的作用,將問題分解為四部分求解,分析水對隧道圍巖應力的影響;楊公標等[21]研究在重力影響下的含空洞地層淺埋隧道圍巖應力及位移解析解.

最近,許多學者對隧洞不均勻邊界條件問題展開了研究.Lu 等[22-24]、Zeng 等[25-26]和Cai 等[27]研究了實際開挖過程引起的邊界不平衡力分布,提出不平衡力系由整個開挖邊界決定,并給出淺埋非圓孔隧道映射函數,對地表有水平初始應力作用的淺埋隧道、水工隧道以及馬蹄形非圓孔隧道等工程問題進行了研究.申航等[28]研究矩形孔口徑向位移邊界條件下的應力場和位移場;文明[29]提出在重力作用下淺埋非圓形隧道力學分析的解耦保角映射方法.Lu 等[30]和Kong 等[31-33]利用傅里葉級數來表示圓形淺埋隧道橫截面變形的統一位移函數,在地表水平以及邊坡邊界模型中考慮隧道“浮力效應”得到單孔隧道解析解,并在此基礎上采用Schwartz 交替法得到雙孔隧道解答;采用傅里葉級數表示出水下淺埋隧道周邊復雜應力狀態的統一應力函數,以此作為孔口應力邊界條件,給出了應力和位移的彈性解析解.

采用平面彈性復變函數方法求解淺埋隧道問題一直存在復勢函數系數求解困難的問題.Verruijt 等[12]在求解孔口應力邊界條件問題時借助復勢函數 φ(ζ)和 ψ (ζ)在像平面圓環域內的收斂性,假設1 階正冪項系數為0,利用迭代方程組經過多次迭代(比如1000 或10000 次)直至高冪次項系數為非零常數,以此得到其余系數.Lu 等[30]推導時取 2N個線性組方程求 2N個未知數,并假定其余正負高冪次項為零,其中N需要足夠大才能保證足夠的精度.

本文利用無窮遠點應力有界性對解析函數的Laurent 級數展開式冪次項進行了確定,從而得到1 階正負冪次項系數,再根據邊界條件所得迭代公式解出完整的解析函數.此方法區別于以上兩種方法,通過從低次冪系數迭代求解至高次冪系數,復勢函數每一項系數都能夠用公式明確表示,使淺埋隧道復勢函數求解更方便,更容易實現編程計算.以淺埋海底隧道孔口受不均勻應力邊界條件為例,得出孔口應力冪級數解,與有限元數值解對比驗證冪級數解的可靠性;研究了復勢函數 φ (ζ)和 ψ (ζ)級數展開后取不同冪次項和不同函數系數對求解結果的影響機理;分析了淺埋隧道埋深對環向壓應力的影響.此方法可求解半無限平面含孔洞問題,求解結果具有較高精度,將程序嵌入小型計算儀器中可方便工程應用.

1 問題描述和基本理論

根據文獻[11-32]提出的海底隧道力學模型可知,所研究的淺埋單孔海底隧道位于完全飽和且均勻、各向同性的彈性含水地層中,力學模型如圖1 所示.水平地表上作用有深度為hw,容重為 γw的海水,海水處于不可壓縮的穩定狀態.隧道半徑為r,隧道中心至地面距離為h,巖石容重為 γ,側壓力系數為k0.將該力學模型分解為兩部分(圖2 和圖3),其關系式如下

圖1 海底隧道力學模型Fig.1 Mechanical model of subsea tunnel

圖2 第1 部分模型Fig.2 Part I model

圖3 第2 部分模型Fig.3 Part II model

第1 部分解答為

運用平面彈性復變函數法研究第2 部分模型的圍巖應力.第2 部分為半無限平面含單圓孔模型,共包含地表和孔口兩個邊界條件,對于巖體任意點的應力和位移均可用復勢函數 φ (z) 和 ψ (z)表示.可求得第2 部分應力解

式中,G為剪切模量,G=E/[2(1+μ)],對于平面應變問題K=3-4μ,μ 為泊松比,k0=μ/(1-μ).

復勢函數 φ (z) 和 ψ (z)是區域內的單值解析函數,其形式如下

其中,φ0(z) 和 ψ0(z)為單值解析函數,可將其展開成Laurent 級數.在本文中考慮巖體自重和水壓力,則Fx=0,Fy=πr2(γ+γw).

2 復數形式的邊界條件

2.1 應力邊界條件和映射函數

采用平面彈性復變函數解法的應力邊界條件如下

其中,C為積分得到的復常數,Xn和Yn分別表示邊界上一點沿x軸與y軸平行的面力分量,l和m是所討論邊界點的單位外法向量的方向余弦,如圖4 所示.

圖4 隧道孔口任意點應力狀態Fig.4 Stress state at any point of tunnel orifice

第2 部分解中地表沒有施加載荷,則地表邊界L1上有Xn=0,Yn=0,此時不失計算結果的有效性,復常數C可以假設為0,地表邊界的應力邊界條件可以表示為

圖2 中的虛線表示未開挖的虛擬孔口,其與圖3兩模型疊加得到開挖后的模型如圖1,開挖后孔口面力分量為0.在隧道孔口邊界L2上施加面力分量Xn和Yn,該面力分量與圖2 中虛擬孔口邊界上面力分量大小相等,方向相反,孔口邊界L2的應力邊界條件表達式中復常數C不能再假設為0.

運用保角變換方法通過映射函數z=ω(ζ)將圖3中Z平面中巖土區域映射到復平面 ζ區域上的圓環中,圓環外徑為1,內徑為 α.α是與隧洞埋深h和隧洞半徑r有關的系數,其可通過如下確定.

將映射函數z=ω(ζ) 代入 φ0(z) 和 ψ0(z)中有

2.2 地表邊界條件求解

為方便求解和表示復平面圓環上的點,令ζ=ρσ, ρ為半徑,σ=eiθ′.圖3 中z平面地表邊界y=0 映射到 ζ平面的單位圓|ζ|=1,地表無窮遠點對應點 ζ=1,則地表邊界在復平面中用 ζ=σ表示.將φ(z) 和 ψ(z) 代入地表應力邊界條件得到z平面的表達式

將映射函數z=ω(ζ)代入式(20),其中

地表邊界 ζ=σ,可將式(20)左右各轉換為 ζ平面中有關 σ的表達式,分別如下所示

Ak可展開為

將式(23)與式(25)相應的k階冪指數對應有

為確定函數 φ0(ζ) 和 ψ0(ζ)中ak與bk的系數還需要通過孔口邊界條件得到其余關系式.

2.3 孔口邊界條件求解

隧道孔口x2+(y+h)2=r2經過映射成復平面ζ上半徑為 α的單位圓 ζ=ασ.在孔口邊界為

將復勢函數 φ (z),ψ (z)和映射函數式(16) 代入邊界條件(27)變換為有關 σ的表達式

將式(29)左側展開

式(29)右側用級數形式表示

將式(26)代入式(32)中可以得到迭代公式如下

當k=0和k=1時

可求得復常數C,其實部與虛部為

其中Ak已由地表邊界得出,Bk可根據隧道孔口邊界求得.將式(2)、式(13)和映射函數(16)代入孔口應力邊界條件(27)右側積分項中并將其展開為有關 σ的多項式得到式(39),則式(29)左側可以展開為式(40),式(40)通過同冪次項系數對應可得到級數系數,復常數C已知從而求得Bk

3 解析函數求解與說明

將復勢函數 φ0(ζ) 和 ψ0(ζ)Laurent 級數展開,其中未知系數ak和bk可由迭代式(33)和式(34)求出,但當k=0和k=1時得到相同的等式,則無法直接求解a-1和,其余ak和bk也無法根據迭代公式求得.為解決這一問題,利用無窮遠點應力有界性提出一種求解方法.在圖3 模型中僅有隧道孔口邊界L2施加載荷,則距離邊界L2無窮遠的巖體應力有界.為保證無窮遠點應力有界性這一條件,根據式(3)和式(4),則復勢函數可表達為

3.1 解析函數求解

將z=w(ζ)代入上式,并將映射函數中含 ζ分式級數展開為

從上式可知當k≤-2 時有ak=0和bk=0.

由式(35)得

式(26)中k=0時,有

將式(46)和式(47)實部與虛部分開計算,再采用復變函數方法能夠求出冪級數解工況中 Re(a-1)值為趨于0 的極小值,假定 Re(a-1)=0,則a1實部 Im(a-1)和 Im(a1)可由下式準確求得

根據迭代式(33),有

根據地表邊界條件(26)可求得系數bk

式(48)～式(53) 即為復勢函數 φ0(ζ)和ψ0(ζ)Laurent 級數展開中各未知系數的求解公式.

3.2 求解方法說明

從映射函數式中可知點 ζ=1為奇點,由于這一奇點的存在給解析計算帶來很大的困擾.從3 個不同的方向趨近于該奇點分別對應物理平面中巖體的3 個無窮遠處.在 ζ平面單位圓 ζ=eiθ′上,當 θ′→0和 θ′→2π其分別對應z平面地表邊界無窮遠點x=+∞和x=-∞處,ζ沿ξ軸ζ →1對應無窮遠點y=-∞處.同時表達式1/(1-ζ)在 |ζ|＜1區域內收斂,可以表示為冪級數,但其在奇點所處的單位圓上,不能將其展開為冪級數而無法采用冪級數解法.

為保證求解精準性和計算的可行性,我們重新定義地表邊界L*(如圖5 所示)為z平面中y→0的直線和 ζ平面中ζ=(1-t)σ,t∈(0,0.1)的近似單位圓.當t足夠小時,新定義的地表邊界L*將趨近于原地表邊界L1,這不僅能夠確保邊界條件與原地表的邊界條件相同,還能確保巖體映射后的討論域 ζ均在表達 式 1/(1-ζ)的收斂域內.雖然 ζ平面上邊界L*(ρ=1-t)和邊界L1(ρ=1)所表示的物理含義不同,但在解析計算時可以將 ρ=1-t近似為 ρ=1,這并不影響最后的計算結果.

圖5 L*邊界示意圖Fig.5 L* boundary diagram

4 應力求解

將映射函數代入式(3)和式(4)得到迪卡爾坐標系下的應力求解公式

在極坐標中各應力分量的關系式為

將復勢函數 φ (ζ) 和 ψ (ζ)代入式(54)～式(57)即可得到不同坐標系下的應力分量.

5 不同方法對比

當本文淺埋海底隧道模型所含參數hw=0 和γw=0 時即可簡化為文獻[22]中考慮重力和任意側向應力作用下的彈性半平面含圓孔模型.為了驗證本文方法的準確性,對側壓力系數k0=0.5下的工況求解,將其與Lu 等[22]得出的結果和有限元結果進行對比.復勢函數 φ0(ζ) 和 ψ0(ζ)分別計算至a6和b5.分析孔口應力時所取計算范圍與文獻[22]相同即隧道右半部分,α=0°為隧道底部,α=180°為隧道頂部.其余參數也與文獻[22] 中參數相同: μ=0.3,E=10.0 MPa,γ=25.0 kN/m3和h/r=2.0.極 坐標系中的計算結果對比見圖6.

圖6 不同方法計算結果對比Fig.6 Comparison of calculation results by different methods

從圖6 可知,本文解得的復勢函數所計算的孔邊環向應力與Lu 等[22]得出結果和有限元結果基本相同.與其他兩種方法相比,本文計算結果壓應力偏大.當k0=0.5,h/r=2.0時,Lu 等[22]計算復勢函數需求解32 元線性方程組.本文只需根據復勢函數系數的顯式求解式(52)和式(53)分別求出系數a6與b5,根據不同精度需求可確定k的取值,k取值越高精度越高,同時本計算過程也減少了依靠地表和孔口邊界條件所求系數Ak和Bk的需求.因此,在滿足工程精度需求下,運用復勢函數顯式計算公式求解,計算過程簡便,計算量小.

6 算例分析

為了驗證本文理論方法的可靠性,將本文在迪卡爾坐標系下的理論解答與數值計算結果進行對比.此例定義參數如下: 隧道孔洞半徑r=4 m,埋深h=10 m,海水深度hw=2 m,海水重度 γw=10 kN/m3,土體重度 γ=25 kN/m3,泊松比μ=0.3,彈性模量E=200 MPa.規定應力拉為正,壓為負.利用ANSYS 軟件建立第2 部分解答的有限元模型計算數值解,模型尺寸為100 m×60 m,采用PLANE183單元,共劃分單元3028 個,圖7 為網格劃分情況.模型上邊界為自由邊界,對模型左右邊界施加x向約束,下邊界施加y向約束,隧道圓形孔口施加與冪級數解邊界L2相同的面力.

圖7 有限元模型網格劃分結果Fig.7 Finite element model meshing results

本文主要是對第2 部分模型(圖3)采用復變函數方法求解,同樣利用ANSYS 對第2 部分解答模型進行模擬.通過有限元模擬得出應力分量,和的應力云圖,如圖8 所示.利用MATLAB 軟件編程計算得到冪級數解.分別將冪級數解和有限元結果與第1 部分解疊加得到隧道開挖后的圍巖應力.孔邊應力冪級數解和數值解對比見圖9.

(2) 環青藏高原東側第一階梯坡降區的復雜地形地貌，是川藏高速公路崩塌發育的根本原因之一，其地形陡、坡降大，岷江、大渡河流經的高山峽谷區，為崩塌危巖的發生奠定了良好的地形地貌基礎。

圖8 第2 部分模型有限元應力云圖Fig.8 Finite element stress nephogram of part II model

圖9 不同系數的冪級數解和數值解對比Fig.9 Comparison of power series and numerical solutions with different coefficients

6.1 數值解對比

從圖9 可以看出冪級數解(power series solutions 1)與數值解能較好吻合,驗證了本文方法的正確性.當 θ ∈[0°,180°]時為隧道上半部分,冪級數解與數值解得到的孔邊應力(σx,σy和 σxy)基本一致;但當 θ ∈(180°,360°)時為隧道下半部分,冪級數解與數值解具有一定的差異.冪級數解與數值解產生的計算差異是由所選映射函數造成的.與數值解相比,冪級數解得到的 σy略小于數值解,而孔口下半部分各點應力分量 σx與σxy的絕對值則比數值解的絕對值大.冪級數解最終計算結果比數值解相對保守.

6.2 解析函數系數影響分析

為了研究單值解析函數 φ0(ζ)與 ψ0(ζ)中系數ak與bk取不同項數對運算結果的影響,取3 組ak與bk展開不同項數時的冪級數解進行對比.系數展開項數見表1,冪級數解對比見圖9.

表1 項數展開表Table 1 Number of terms expansion table

從圖9 中可以看出power series solutions 1 與power series solutions 3 在3 個應力分量 σx,σy和 σxy的計算結果都基本保持一致,power series solutions 2與其結果相比在不同位置、不同應力分量上均有一定的差距.特別在隧道孔口上半部分 θ ∈[0°,180°],power series solutions 3 要優于power series solutions 2.

從中可以得出單值解析函數 φ0(ζ)對最終結果準確性的影響要大于 ψ0(ζ).求解時為了保證計算結果的準確性,φ0(ζ) 需要展開足夠多的項,ψ0(ζ)則至少要展開b-1～b1項.

6.3 隧道埋深影響分析

為了分析淺埋隧道埋深對環向應力的影響,隧道中心埋深h分別取8,10,12 和14 m,單值解析函數 φ0(ζ) 與 ψ0(ζ)中系數ak與bk展開項與power series solutions 1 相同,其余條件與參數不變,不同隧道埋深孔口應力情況見圖10.隧道底部(θ=270°)及兩側孔腰處(θ=200°或 θ=340°)環向壓應力隨著隧道埋深h增大而增大,而孔口頂部(θ=90°)環向壓應力無顯著變化.隧道埋深增加1.75 倍,腰部及底部環向壓應力增加3～4 倍.隧道腰部與底部環向應力的差值也隨著隧道埋深的增大而增大.由此可見隧道埋深對隧道底部及腰部環向壓應力的大小影響較大,且對腰部的影響最為顯著.

圖10 不同埋深孔口環向應力對比Fig.10 Comparison of annular stress at orifices with different burial depth

7 結論

本文運用平面彈性復變函數方法研究了淺埋圓孔海底隧道在圍巖重力和海水靜水壓力作用下的圍巖應力冪級數解.采用分式映射函數將半無限平面含圓孔模型映射為像平面上的圓環域,在圓環域內將復勢單值解析函數展開為Laurent 級數.利用無窮遠點應力有界性對級數展開式的冪次項進行確定,根據級數系數迭代公式得到復勢函數的顯式解,運用應力分量的復變函數表達式得到孔邊各點的應力分量,將冪級數解與有限元結果對比得出以下結論.

(1) 為保證無窮遠點應力有界,根據應力求解公式,將復勢單值解析函數 φ0(ζ) 和 ψ0(ζ)的Laurent 級數重新定義為負一項至正無窮項.在 R e(a-1)=0條件下,由地表與孔口邊界條件得出復勢解析函數系數顯式解,此解可迭代求得其余系數,從而實現級數系數從低次冪迭代至高次冪這一過程,使得淺埋隧道復勢函數求解更方便.

(2)冪級數解與有限元數值解能很好吻合.在隧道上半部分(θ ∈[0°,180°]) 冪級數解與數值解得到的孔邊應力(σx,σy,σxy)基本一致;在隧道下半部分(θ ∈(180°,360°))由于所選映射函數的影響,冪級數解得到的 σy略小于數值解,σx與 σxy的絕對值則比數值解的絕對值大.冪級數解最終計算結果較數值解偏保守.

(3) 單值解析函數 φ0(ζ)對最終結果準確性的影響要大于 ψ0(ζ).求解時為了保證計算結果的準確性,φ0(ζ)需 要展開足夠多項,ψ0(ζ)則 至少要展開b-1～b1項.

(4) 隨著海底隧道埋深增大,隧道底部及兩側孔腰處的環向壓應力隨之增大.隧道埋深增加1.75 倍,腰部及底部環向壓應力增加3～4 倍;而孔口頂部環向壓應力無顯著變化.隧道腰部與底部環向應力的差值也隨著隧道埋深的增大而增大.隧道埋深對隧道底部及腰部環向壓應力的大小影響較大,且對腰部的影響最為顯著.

數據可用性聲明

支撐本研究的科學數據已在中國科學院科學數據銀行(science data bank) ScienceDB 平臺公開發布,訪問地址為https://www.doi.org/10.57760/sciencedb.j00140.00011 或http://resolve.pid21.cn/31253.11.sciencedb.j00140.00011.