深埋馬蹄形隧道開挖圍巖應力與位移的復變函數解

施有志，高軒能

（1.廈門理工學院 土木工程與建筑學院，福建 廈門 361021；2.華僑大學 土木工程學院，福建 泉州 362021）

地下工程穩定性分析和設計中，經常遇到可假定為平面應變問題的單洞及多洞問題。1920年，Jeffery［1］利用雙極坐標法給出了圓形洞室的平面應力和應變的解答，但其解是建立在一些限制條件基礎上的，過程繁瑣且未能給出位移邊界條件下該問題的解答。Mindlin［2－3］同樣采用雙極坐標解法，在應力勢函數的推導過程中加入了介質所受的重力場和初始應力場，從滿足隧道周邊的無法向和切向應力條件開始，提取出包括擴開挖體重量的應力勢函數，然后為了滿足位移單值和水平向無法向和切向應力邊界條件采用試湊的辦法開展一系列工作，由于雙極坐標法本身的局限性最終仍然未能給出位移場的公式；雙極坐標法一般只求得一些簡單孔洞的應力場，不能給出位移場。且當孔洞距離地表較近時，雙極坐標法往往給出錯誤的數值計算結果，從而迫使專家學者們尋求到解決固定邊界條件下半空間含括一個孔洞問題的另外一種更重要的方法，即復變函數法。自1777年Euler首創虛數單位符號“i”，并系統地建立復變函數理論以來［4］，到20世紀復變函數已經形成了非常系統的理論，并已經被廣泛地應用到彈性力學領域，如徐芝綸［5］、Timoshenko等［6］、Muskhelishvili［7］、路見可［8］等的力學專著中都有相關論述。復變函數求解單孔洞問題可以歸結為2種方法：柯西積分法和解析延拓法。Verruijt［9］采用柯西積分法對半平面含括單個孔洞的問題作了一系列研究。1997年他給出了滿足水平面無法向和切向應力的第一邊界條件和孔洞周邊的位移可展開成映射后域內單位復數的冪級數形式的第二邊界條件的解的遞推公式，給出了洞室周邊均布位移條件下解的遞推公式。1998年Verruijt［10］給出了滿足水平面無法向和切向應力的第一邊界條件、孔洞周邊的應力可展成單位復數的冪級數形式的第一邊界條件的解的遞推公式，并給出了洞室周邊在均布法向應力條件下的精確解。但Verruijt的工作僅局限于簡單形狀（圓形、橢圓形）的隧道。陳子蔭［11］，呂愛鐘［12］用柯西積分法求解出了任意形狀孔洞在半空間的應力解和位移解，但他們的模型邊界條件沒有考慮隧道實際施工情況，沒有考慮隧道開挖過程中襯砌力對圍巖的影響，且因其解的表達式復雜，僅給出了解的隱式表達，對于工程問題中常用的均為復雜的孔口問題（如馬蹄形斷面孔口）并沒有給出確定的顯式表達。王志良［13］、張頂鋒［14］、晏莉等［15］的推導，則直接將隧道簡化為圓形。童磊［16］推導了圓形隧道基于任意襯砌變形邊界條件的復變函數彈性解，來預測軟土中隧道開挖時短期地表豎向沉降與側向位移。筆者在呂愛鐘等［12］的研究基礎上，利用復變函數中的柯西積分法求解工程問題中常用的單心圓仰拱馬蹄形隧道在二維平面彈性半空間內任意一點處的值和位移值解析解；結合馬蹄形隧道的典型斷面，采用有限元數值分析來檢驗解析解精度，驗證了結果的精確性。

1 馬蹄形隧道開挖的柯西積分法求解

地下洞室埋置深度與孔徑比較大時，可不考慮重力梯度的影響，把重力作用化為無限遠處作用有P1、P2的外載來求解，圖1給出了z平面半無限空間單心圓仰拱隧道馬蹄形的構形，圖2為通過復變函數中保角映射后把z平面單心圓仰拱馬蹄形外隧道域轉化為ξ平面單位圓外域后的情況，其中P1、P2表示無窮遠處水平原巖和豎向原巖的應力值。

圖1 分析模型示意圖

圖2 映射函數與逆映射函數

式中：系數a0，a1…an和b0，b1…bn為單值解析函數的常數項系數，由具體邊界條件和復平面點所在的具體位置求得。

如圖2所示，將半無限平面z上單心圓仰拱馬蹄形隧道孔口的外域映射到ξ平面的單位圓外域，其映射函數采用最普遍的形式，可用Laurent表示為［12］

式中：R反映孔洞的大小；C1，C2，……，Cn反映了孔洞的形狀；C0反映了孔洞所處坐標系中的位置；z和ζ分別為物理平面和映射平面上的復數坐標。

根據已經確定映射函數z＝w（ζ），把z物理平面的應力函數轉化到ζ平面的應力函數，為此引用如下記號：

考慮邊界條件為無窮遠點剛體轉動為0，孔洞周邊沒有作用面力，則根據邊界條件：

把式（9）、（10）帶入式（5）和（6），經過化簡可求得應力函數

2 平面內應力解、位移解表達式

在求得應力函數后，利用到正交曲線坐標系下求解應力分量的表達式

式中：σθ、σρ和τρθ分別為正交曲線坐標系下環向應力、徑向應力和剪切應力；ρ為映射平面上徑向坐標；接著根據的求的σθ、σρ和τρθ求出在物理坐標下的σx、σy和τxy。

由應力邊界條件和位移邊界條件確定的應力復變函數有所差別，利用應力邊界確定的兩應力復變函數來計算位移場時會使整個圍巖多出一個剛性位移。研究表明，此剛性位移可通過一定的邊界條件予以消除，即距離兩洞足夠遠處（理論上應該是無窮遠處）的位移為0。以應力函數φ1▽（z1）、ψ1▽（z1）表示消除了剛性位移后的應力函數表達式，在第1映射坐標系x1O1y1，φ1▽（z1）、ψ1▽（z1）按照下式確定：

式中：φ1（z1）、ψ1（z1）為應力邊界條件下求出的復應力函數；m、n、m′、n′均為實常數，其具體大小由帶入的邊界條件確定。

根據位移場邊界條件式就可確定下來圍巖內任一點的位移：

式中u、v分別表示兩坐標軸方向即水平方向位移和垂直方向位移。

3 馬蹄形隧道映射函數參數的確定

根據式（4）求映射函數最終化為求解映射函數中的參數R和Ck，求解步驟如下：

1）在z平面選定直角坐標系，使單心圓仰拱隧道對稱于x軸分布，在x軸左邊圓周上逆時針任取點i，使得點i的直角坐標和相應極坐標為（xi，yi）和 （r0，ai），令i＝0對應坐標點 （r0，0），ζ平面單位圓圓周上對應的點的直角坐標和極坐標均為 （1，0），則最終可推出：

2）用復合最優化技術［17］來確定Ck，設目標函數為：

以Ck、Bi為未知量，通過最優化擬合來確定參數Ck；根據驗證，在式（4）映射函數的級數中只需取很少幾項計算的結果就已經足夠精確。經驗證，本例的中的單心圓仰拱馬蹄形隧道斷面（如圖3所示）的映射函數中取k＝5就已經足夠精確。

圖3所示為某工程單心圓仰拱隧道的典型斷面圖，其中R1＝6m，R2＝15m，無窮遠處作用的豎向地表外荷載取P2＝20MPa，彈性模量取E＝4.3GPa，泊松比μ＝0.28。水平向為固定邊界條件（P1＝0）。其余尺寸如圖3所示，經過最優化計算得單心圓仰拱隧道的保形變換函數為：

圖3 單心圓仰拱隧道斷面尺寸

將保形變換函數的系數Ck值代入式（4）、然后代入式（5）、（6），最終求的應力函數 φ（ζ）和 φ（ζ）。把應力函數帶入式（18）、（21）最終求得單心圓仰拱隧道的 （σρ，σθ，τρθ）和位移 （u ，v）。

4 解析解與有限單元法計算結果比較

采用三維有限元分析軟件ANSYS建立二維平面應變模型，對理論推導單心圓仰拱隧道解析解公式進行驗證，有限元模型如圖4。為了與解析解計算結果進行對比分析，模型中取單一巖層進行分析，根據相關工程經驗，參數在合理的經驗范圍內假定取值，巖層重度取24kN／m3，巖層彈性模量取E＝4.3GPa，豎向地表外荷載取P2＝2MPa，泊松比μ＝0.28。水平向為固定邊界條件（P1＝0）。模型計算范圍水平方向左右兩邊各取3倍洞跨寬；垂直方向上下邊界也各取3倍洞跨寬。為了便于分析，把單心圓隧道用角度劃分成12等份，分別比較了解析解所得的隧道洞周附近區域各的σr、σθ和τrθ等應力值和u、v等位移值（取上半平面），見圖5～7，并對其誤差大小進行分析，見表1。由計算結果可以看各點的周邊的應力值與位移值變化趨勢基本一致，隧道環向應力值σθ解析解計算的最大值為5.682MPa，最大值出現在逆時針方向120°，環向應力值的最大誤差為僅為8.8%，水平位移最大誤差為11.3%，而豎向位移的最大誤差為7.3%。結果證明了公式的可靠性。

表1 解析解和有限元結果比較

圖4 有限元模型

圖5 隧道周邊σθ值結果比較

圖6 隧道周邊u值結果比較

圖7 隧道周邊v值結果比較

5 結 語

分析馬蹄形這種不規則形狀的隧道，一般需要試驗或有限元法分析，筆者直接用隱式解析解進行分析。將深埋單心圓仰拱馬蹄形隧道所處的地下空間轉換為平面應變情況下的彈性半空間，考慮隧洞埋深與孔徑相比比較大從而不考慮重力梯度影響，將重力作用化為作用于遠場邊界上的外荷載，利用復變函數柯西積分法結合最優化理論求解出了地下半空間典型斷面的單心圓仰拱馬蹄形隧道在其所在的彈性半空間內任意一點的應力解和位移解解析解的表達式。解析解結果與有限元計算結果的驗證表明公式的正確性。

筆者的求解方法拓展了柯西積分法的應用，但還有若干問題須進一步研究，如考慮巖土體為粘彈性情況。同時，本文解只針對深埋隧道，淺埋馬蹄形斷面隧道的解的表達式也有待于進一步研究。

［1］Jeffery G B.Plane stress and plane strain in Bipolar Coordinates ［J］.Transaction of the Royal Society，SeriesA，1920，221：265－293.

［2］Mindlin R D.Stress distribution around a tunnel［J］.Transaction of the ASCE，1940，190：1117－1153.

［3］Mindlin R D.Stress distribution around a hole near the edge of a plate under tension［J］.Proceedings of the Society Experimental Srress Analysis，1948（5）：56－67.

［4］Green A E.General biharmonic analysis for a plate containing circular holes［J］.Proceedings of the Royal Society of London Series A：Mathematical and Physical Engineering Sciences，1940，964（176）：121－139.

［5］徐芝綸.彈性力學［M］.北京：高等教育出版社，2004.

［6］Timoshenko S P，Goodier J N.Theory of Elasticity［M］.北京：清華大學出版社，2004.

［7］Muskhelishvili N I.Mathematical theroy of elasticity［M］.Leyden：International Publishing，1954.

［8］路見可.平面彈性復變方法［M］.武漢：武漢大學出版社，2002.

［9］Verruijt A.A complex variable solution for a deforming circular tunnel in an elastic half plane ［J］.Geotechnique，1997，21（4）：77－89.

［10］Verruijt A.Deformations of an elastic half plane with a circular cavity［J］.International Journal of Solids and Structures，1998，35（21）：2795－2804.

［11］陳子蔭.圍巖力學分析中的解析方法［M］.北京：煤炭工業出版社，1994：70－76.

［12］呂愛鐘，張路青.地下隧洞力學分析的復變函數方法［M］.北京：科學出版社，2007：51－73.

［13］王志良，申林方，姚激，等.淺埋隧道圍巖應力場的計算復變函數求解法［J］.巖土力學，2010，31（Sup1）：86－90.Wang Z L，Shen L F，Yao J，et al.Calculation of stress field in surrounding rocks of shallow tunnel using computational function of complex variable method［J］.Rock and Soil Mechanics，2010，31（Sup1）：86－90.

［14］張頂鋒，杜守繼.盾構隧道并行施工對既有隧道影響的理論分析［J］.力學季刊，2010，31（4）：582－588 Zhang D F，Du S J.Theoretical analysis of existing tunnel affect under shield tunnel parallel construction［J］.Chinese Quarterly of Mechanics，2010，31（4）：582－588.

［15］晏莉，陽軍生，劉寶琛.淺埋雙孔平行隧道開挖圍巖應力和位移分析［J］.巖土工程學報，2011，33（3）：413－419 Yan L，Yang J S，Liu B C.Stress and displacement of surrounding rock with shallow twin－parallel tunnels［J］.Chinese Journal of Geotechnical Engineering，2011，33（3）：413－419.

［16］童磊，謝康和，盧萌盟，等.盾構任意襯砌變形邊界條件下復變函數彈性解［J］.浙江大學學報：工學版，2010，44（9）：1825－1830.Tong L，Xie K H，Lu M M，et al.Elastic complex variables solution for general arbitrary ground deformation of tunnels in clays［J］.Journal of Zhejiang University：Engineering Science，2010，44（9）：1825－1830.

［17］劉偉.信機械最優化理論［M］.北京：清華大學出版社，1997：143－149.