施有志，高軒能

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

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

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

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

圖1 分析模型示意圖

圖2 映射函數(shù)與逆映射函數(shù)

式中：系數(shù)a0，a1…an和b0，b1…bn為單值解析函數(shù)的常數(shù)項(xiàng)系數(shù)，由具體邊界條件和復(fù)平面點(diǎn)所在的具體位置求得。

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

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

根據(jù)已經(jīng)確定映射函數(shù)z＝w（ζ），把z物理平面的應(yīng)力函數(shù)轉(zhuǎn)化到ζ平面的應(yīng)力函數(shù)，為此引用如下記號(hào)：

考慮邊界條件為無窮遠(yuǎn)點(diǎn)剛體轉(zhuǎn)動(dòng)為0，孔洞周邊沒有作用面力，則根據(jù)邊界條件：

把式（9）、（10）帶入式（5）和（6），經(jīng)過化簡(jiǎn)可求得應(yīng)力函數(shù)

2 平面內(nèi)應(yīng)力解、位移解表達(dá)式

在求得應(yīng)力函數(shù)后，利用到正交曲線坐標(biāo)系下求解應(yīng)力分量的表達(dá)式

式中：σθ、σρ和τρθ分別為正交曲線坐標(biāo)系下環(huán)向應(yīng)力、徑向應(yīng)力和剪切應(yīng)力；ρ為映射平面上徑向坐標(biāo)；接著根據(jù)的求的σθ、σρ和τρθ求出在物理坐標(biāo)下的σx、σy和τxy。

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

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

根據(jù)位移場(chǎng)邊界條件式就可確定下來圍巖內(nèi)任一點(diǎn)的位移：

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

3 馬蹄形隧道映射函數(shù)參數(shù)的確定

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

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

2）用復(fù)合最優(yōu)化技術(shù)［17］來確定Ck，設(shè)目標(biāo)函數(shù)為：

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

圖3所示為某工程單心圓仰拱隧道的典型斷面圖，其中R1＝6m，R2＝15m，無窮遠(yuǎn)處作用的豎向地表外荷載取P2＝20MPa，彈性模量取E＝4.3GPa，泊松比μ＝0.28。水平向?yàn)楣潭ㄟ吔鐥l件（P1＝0）。其余尺寸如圖3所示，經(jīng)過最優(yōu)化計(jì)算得單心圓仰拱隧道的保形變換函數(shù)為：

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

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

4 解析解與有限單元法計(jì)算結(jié)果比較

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

表1 解析解和有限元結(jié)果比較

圖4 有限元模型

圖5 隧道周邊σθ值結(jié)果比較

圖6 隧道周邊u值結(jié)果比較

圖7 隧道周邊v值結(jié)果比較

5 結(jié) 語

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

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

［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］徐芝綸.彈性力學(xué)［M］.北京：高等教育出版社，2004.

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

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

［8］路見可.平面彈性復(fù)變方法［M］.武漢：武漢大學(xué)出版社，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］陳子蔭.圍巖力學(xué)分析中的解析方法［M］.北京：煤炭工業(yè)出版社，1994：70－76.

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

［13］王志良，申林方，姚激，等.淺埋隧道圍巖應(yīng)力場(chǎng)的計(jì)算復(fù)變函數(shù)求解法［J］.巖土力學(xué)，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］張頂鋒，杜守繼.盾構(gòu)隧道并行施工對(duì)既有隧道影響的理論分析［J］.力學(xué)季刊，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］晏莉，陽軍生，劉寶琛.淺埋雙孔平行隧道開挖圍巖應(yīng)力和位移分析［J］.巖土工程學(xué)報(bào)，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］童磊，謝康和，盧萌盟，等.盾構(gòu)任意襯砌變形邊界條件下復(fù)變函數(shù)彈性解［J］.浙江大學(xué)學(xué)報(bào)：工學(xué)版，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］劉偉.信機(jī)械最優(yōu)化理論［M］.北京：清華大學(xué)出版社，1997：143－149.