陳文博，房倩，張頂立，陳炫皓

(北京交通大學(xué) 城市地下工程教育部重點(diǎn)實(shí)驗(yàn)室，北京，100044)

非圓形隧道力學(xué)的解析方法能通過顯式解快速計(jì)算巖體任意位置處的力學(xué)響應(yīng)，克服數(shù)值軟件計(jì)算結(jié)果受單元尺寸、網(wǎng)格劃分、邊界距離和插值計(jì)算影響的缺點(diǎn)，在提取應(yīng)力跡線、分析隧道承載拱效應(yīng)等研究上具有應(yīng)用潛力和相對優(yōu)勢，是分析深部巖體內(nèi)部復(fù)雜力學(xué)行為的有效手段，多年來一直是研究的難點(diǎn)和熱點(diǎn)?；趶?fù)分析理論對非圓形隧道進(jìn)行力學(xué)解析的關(guān)鍵在于：1) 求解襯砌、圍巖等區(qū)域的保角映射函數(shù)；2) 求解各區(qū)域內(nèi)的應(yīng)力函數(shù)。

針對隧道力學(xué)解析中的保角映射，呂愛鐘等[1]提出復(fù)合形法、混合罰函數(shù)法等優(yōu)化方法；祝江鴻等[2]提出奇偶插值點(diǎn)相互迭代的計(jì)算方法；QI等[3]提出求解矩形隧道圍巖保角映射函數(shù)的螢火蟲尋優(yōu)算法；HE 等[4]基于黎曼映射定理和邊界對應(yīng)定理，提出求解不規(guī)則隧道保角映射函數(shù)的迭代算法。這些研究未對該保角映射方法的復(fù)分析理論基礎(chǔ)進(jìn)行嚴(yán)謹(jǐn)?shù)耐茖?dǎo)，認(rèn)為映射前后區(qū)域的邊界點(diǎn)相互對應(yīng)(無限接近)即為保角映射，所求映射函數(shù)僅具有將一個(gè)邊界映射為另一個(gè)邊界的功能，不能保證區(qū)域內(nèi)映射的保角性[5-6]。

在非圓形有支護(hù)隧道彈性力學(xué)解析的研究中，KARGAR 等[7-14]采用文獻(xiàn)[1-3]中保角映射函數(shù)的求解方法，基于MUSKHELISHVILI等[15]的復(fù)分析理論，探索有襯砌支護(hù)、注漿圈加固的非圓形隧道力學(xué)解析，利用復(fù)合閉路定理和解析域內(nèi)的冪級數(shù)性質(zhì)，提出了“同冪次消元”的應(yīng)力函數(shù)求解方法。文獻(xiàn)[7-14]對襯砌域、加固圈域和圍巖域采用同一個(gè)保角映射函數(shù)，使得相鄰邊界上不同區(qū)域的點(diǎn)在映射后仍重合，以便采用“消元法”求解應(yīng)力函數(shù)。由保角映射的唯一性定理[6]可知，僅具有一條共用邊界的有限雙連通域(襯砌)和無限雙連通區(qū)域(圍巖)的保角映射函數(shù)是不同的(嚴(yán)格意義上講，即不等價(jià)，2個(gè)保角映射函數(shù)不能通過旋轉(zhuǎn)、縮放和平移等線性保角變換相互轉(zhuǎn)化)。因此，文獻(xiàn)[7-14]中求應(yīng)力函數(shù)的方法是不適用的。

針對以上保角映射問題和應(yīng)力函數(shù)求解問題，首先，引進(jìn)了嚴(yán)格的基于邊界積分的保角映射方法，并編制了適用于隧道的保角映射程序；其次，根據(jù)隧道的邊界條件，提出了求解圍巖和襯砌應(yīng)力函數(shù)的多點(diǎn)逼近方法；再次，通過多心圓隧道的算例，驗(yàn)證了開發(fā)程序所求映射函數(shù)的保角性和保圓性，分析了圍巖與襯砌接觸面邊界上點(diǎn)在保角映射后的位置差異；最后，對比隧道力學(xué)解析和數(shù)值模擬計(jì)算結(jié)果，分析了側(cè)壓力系數(shù)為0.50時(shí)隧道的力學(xué)響應(yīng)規(guī)律，并驗(yàn)證了本文提出非圓形深埋隧道力學(xué)解析方法的正確性。

1 雙連通無限域、有限域的保角映射方法

圖1 所示為保角映射與一般映射差異說明圖。從圖1 可見：內(nèi)邊界為L1的雙連通無限域S1經(jīng)保角映射函數(shù)映射為內(nèi)邊界半徑為R的圓形雙連通無限域；外邊界為L1，內(nèi)邊界分別為L21、L22和L23的雙連通有限域S21、S22和S23分別經(jīng)保角映射函數(shù)和映射為外徑R和內(nèi)徑分別為r1、r2和r3的共心圓環(huán)。由單連通域及多連通域保角映射的唯一性定理[6,16]可知：互不等價(jià)，r1、r2和r3互不相等。映射函數(shù)和雖然同樣可以將無限域S1的邊界L1映射為半徑為R的圓，但均不是無限域S1的保角映射函數(shù)，僅有是其保角映射函數(shù)。

圖1 保角映射與一般映射差異說明圖Fig.1 Illustration of difference between conformal mapping and general mapping

文獻(xiàn)[1-3]所提保角映射方法存在問題的原因在于：

1) 曲解了黎曼存在性定理，將單連通域保角映射的邊界對應(yīng)定理錯(cuò)誤地應(yīng)用到雙連通無限域、有限域的保角映射上，單純地將邊界之間的映射作為兩區(qū)域之間的保角映射。

2) 其優(yōu)化方法、迭代計(jì)算方法僅可實(shí)現(xiàn)邊界之間的近似映射，不能保證區(qū)域映射的保角性和映射點(diǎn)之間的一一對應(yīng)關(guān)系。

因此，本文進(jìn)一步基于保角映射的復(fù)分析理論和數(shù)值計(jì)算方法，研究非圓形隧道襯砌域、圍巖域映射為共心圓環(huán)域、含圓無限域的保角映射。關(guān)于保角映射的復(fù)分析理論，CROWDY 等[17-20]基于Schwarz-Christoffel 公式的經(jīng)典函數(shù)論，引入Schottky-Klein 素函數(shù)，針對以圓為邊界的多連通域到多邊形域、狹縫域、多圓弧域的保角映射，建立了高階微分方程；NASSER[21]通過邊界積分的數(shù)值方法求解此微分方程，進(jìn)一步將保角映射問題轉(zhuǎn)化為邊界積分方程的求解問題。在數(shù)值計(jì)算方法方面，NASSER 等[22]提出了一種邊界積分方法解決Laplace 方程在多連通區(qū)域上的Dirichlet 邊界條件問題和Neumann 邊界條件問題，建立了第二、三、四、五類Koebe標(biāo)準(zhǔn)狹縫域保角映射的數(shù)值計(jì)算方法[23-24]，提出了含狹縫的多連通域的邊界積分方法[25]。NASSER 等[26-27]基于Nystr?m 方法對邊界積分方程進(jìn)行離散化，提出了一種求解含廣義Neumann 核函數(shù)和共軛廣義Neumann 核函數(shù)邊界積分方程的數(shù)值計(jì)算方法。

基于文獻(xiàn)[20]和[27]，本文采用MATLAB語言編寫程序，將以多段圓弧為邊界的多連通域映射為以圓為邊界的多連通域。對于如圖2所示的有限雙連通域，首先利用基于邊界積分方程的保角映射方法，將復(fù)平面Z上的有限雙連通域邊界上各點(diǎn)保角映射為復(fù)平面η上偏心圓環(huán)上的點(diǎn)；再利用M?bius 變換(分式線性變換)[16,28]，得到在ζ平面半徑分別為R和r的共心圓環(huán)上相對應(yīng)的各點(diǎn)，映射前后邊界上各點(diǎn)一一對應(yīng)。

圖2 不同復(fù)平面上襯砌邊界Fig.2 Lining boundary on different complex planes

為求得如式(1)Laurent 級數(shù)形式的保角映射函數(shù)，設(shè)待求映射函數(shù)的項(xiàng)數(shù)為2m1+1，在Z平面和ζ平面的邊界上選取n1對點(diǎn)，將點(diǎn)坐標(biāo)zt和ζt(t=1，2，…，n1)代入式(1)，建立求常系數(shù)x1的矩陣方程，見式(3)。調(diào)整m1和n1，即可獲得滿足精度要求的保角映射函數(shù)ω*(zt)。ω*(zt)的逆函數(shù)ω(ζt)也為保角映射函數(shù)(式(2))，將ζ平面圓環(huán)域上點(diǎn)映射到Z平面上襯砌區(qū)域上的點(diǎn)，其求法與ω*(zt)相同。

式中：ak和bk為復(fù)常系數(shù)。

圍巖區(qū)域的保角映射函數(shù)求解原理和方法與襯砌相同，但不需要圖2(b)的分式線性變換。

2 深埋隧道襯砌與圍巖的邊界條件

圖3 所示為襯砌邊界到同心圓環(huán)的映射示意圖。如圖3所示，L1和分別為襯砌內(nèi)邊界和外邊界，為圍巖內(nèi)邊界，圍巖外邊界為無限遠(yuǎn)處。保角映射函數(shù)ω1和ω2分別將ζ平面上的同心圓環(huán)域和含圓無限域映射為Z平面上的襯砌區(qū)域S1和圍巖區(qū)域S2。在邊界L1上，z11=ω1(ζ11)；在邊界上，z21=ω1(ζ21)；在邊界上，z22=ω2(ζ22)；其中，z11、z21和z22分別為Z平面上襯砌內(nèi)邊界、外邊界和圍巖邊界上的點(diǎn)，ζ11、ζ21和ζ22分別為ζ平面上襯砌內(nèi)邊界、外邊界和圍巖內(nèi)邊界上的點(diǎn)。

圖3 襯砌邊界到同心圓環(huán)的映射示意圖Fig.3 Schematic diagram of mapping of lining boundaries to concentric rings

在ζ平面極坐標(biāo)系下，應(yīng)力和位移可由應(yīng)力函數(shù)φ和ψ表示[29]。

式中：σθ、σρ和τθρ分別為環(huán)向應(yīng)力、軸向應(yīng)力和切應(yīng)力；Re 表示對復(fù)數(shù)取實(shí)部；uρ和uθ分別為軸向位移和環(huán)向位移；G為材料的切變模量，由材料彈性模量E和泊松比μ決定，；κ在平面應(yīng)變條件下取3-4μ。

對于深埋隧道的平面應(yīng)變彈性力學(xué)求解，在ζ平面上，襯砌區(qū)域和圍巖區(qū)域內(nèi)的應(yīng)力函數(shù)可表示為[7]

1) 區(qū)域S1內(nèi)

2) 區(qū)域S2內(nèi)

式中：ck、dk、ek、fk、hk、lk、mk和qk為應(yīng)力函數(shù)的常復(fù)系數(shù)；hk和mk由無窮遠(yuǎn)處地應(yīng)力確定。

在邊界L1上，對于任意ζ11應(yīng)滿足σρ=0和τθρ=0 的應(yīng)力條件，將ζ11代入式(4)和(5)，可得L1上邊界條件。同時(shí)，在圍巖內(nèi)趨于無窮遠(yuǎn)處應(yīng)力函數(shù)應(yīng)滿足式(11)和式(12)。

通過限定接觸面邊界上的應(yīng)力、位移連續(xù)條件及其范圍，可對襯砌與圍巖完全接觸、非完全解接觸或局部脫離接觸的隧道進(jìn)行彈性力學(xué)解析計(jì)算。本文假定圍巖和襯砌完全接觸，既不相互脫離又不相互滑動(dòng)[30]。在Z平面上，圍巖和襯砌接觸面上的法向應(yīng)力、切應(yīng)力和位移相同。由保角變換的性質(zhì)可知，在ζ平面極坐標(biāo)系下，邊界和上對任意滿足ω1(ζ21) =ω2(ζ22)的點(diǎn)ζ21和點(diǎn)ζ22，兩點(diǎn)上的軸向應(yīng)力、切應(yīng)力和位移也相等。將ζ21和ζ22分別代入式(5)和(6)，即可得接觸面邊界上的2個(gè)邊界條件。本研究基于復(fù)變理論進(jìn)一步拓展非圓性非靜水壓力下深埋隧道的彈性力學(xué)解析方法，暫不能考慮圍巖彈塑性、應(yīng)變軟化等問題。

3 應(yīng)力函數(shù)的求解方法

在Z平面邊界L1和L2上分別選取n2個(gè)點(diǎn)，經(jīng)映射得到ζ平面上的3 組點(diǎn)ζ11k=ω1(z1k)，ζ21k=ω1(z2k)，ζ22k=ω2(z2k)，k=1，…，n2。將ζ平面上3n2個(gè)點(diǎn)分別代入對應(yīng)的邊界條件，可提取得到應(yīng)力函數(shù)系數(shù)6組復(fù)數(shù)未知量ck、dk、ek、fk、hk和qk(k=1，…，m1)的3n2個(gè)方程。為能夠使用MATLAB的矩陣運(yùn)算進(jìn)行求解，將復(fù)數(shù)未知量拆分為以實(shí)數(shù)表示的形式，如ck=ck1+ick2(i 為虛數(shù)單位)，并代入方程組。

基于MATLAB 編寫程序，逐個(gè)搜索各等式方程中的常數(shù)項(xiàng)和各實(shí)未知量(ck1、ck2等)的系數(shù)，將各實(shí)未知量對應(yīng)的系數(shù)寫成6n2行、8+12m1列的矩陣Z2，對應(yīng)的實(shí)常數(shù)項(xiàng)依次序?qū)懗?n2行的列向量。

式中：x2=[c01，c02，…，q01，q02，c11，c12，…，q11，q12，…cm11，cm12，…，qm11，qm12]T；Φ2為各方程中常數(shù)項(xiàng)相反數(shù)的列向量。調(diào)整n2和m1，并選取合適方法求解矩陣，即可獲得滿足邊界條件的應(yīng)力函數(shù)系數(shù)。

4 隧道的保角映射及驗(yàn)證

選取某鐵路單線隧道襯砌斷面，其幾何尺寸見圖4，在以O(shè)為原點(diǎn)、x軸為實(shí)軸、y為虛軸的Z平面上，分別對以該斷面襯砌外輪廓為內(nèi)邊界的圍巖無限域、該斷面襯砌的雙連通有限域進(jìn)行保角映射。

圖4 多心圓隧道斷面圖Fig.4 Multi-center circle tunnel cross-section

4.1 圍巖區(qū)域保角映射

采用第3節(jié)保角映射方法，將圖4中圍巖無限域映射為含單圓孔無限域(圓孔半徑為1 m)，求得保角映射函數(shù)，其中z的1 次冪至-6 次冪的系數(shù)見表1。原求得保角映射函數(shù)的級數(shù)共17 項(xiàng)(小于10-4的系數(shù)未列出，表中以缺省符號表示，下同)。在Z平面上圍巖內(nèi)邊界上選取3 072 個(gè)點(diǎn)，經(jīng)過映射函數(shù)映射得到ζ平面上3 072 個(gè)點(diǎn)，各點(diǎn)復(fù)數(shù)模的均值為1.000 0，變異系數(shù)為2.10×10-6，說明該映射函數(shù)的保圓性較好。采用第3節(jié)中保角映射方法，求得的逆函數(shù)、保角映射函數(shù)ω2(ζ)，其中ζ的1次冪至-8次冪的系數(shù)見表1。

表1 圍巖域保角映射函數(shù)的系數(shù)Table 1 Coefficients of conformal mapping function in rock domain

由第3 節(jié)可知，解析計(jì)算中僅需要ω1(ζ)和ω2(ζ)，因此，本文僅檢驗(yàn)該映射函數(shù)的保角性。圖5所示為在ζ平面上取任意夾角的示意圖，為得到各角映射后的角度，令ds足夠小(本文取10-4m)，在夾角直線上各取2個(gè)點(diǎn)，將點(diǎn)a、點(diǎn)b和點(diǎn)c分別映射到Z平面，映射點(diǎn)連線的夾角可以近似認(rèn)為是映射后的曲線夾角。在ζ平面上任意選取8 個(gè)夾角，如表2所示，映射前后夾角的相對誤差量級等于或小于10-6，說明映射函數(shù)的保角性很好。

表2 圍巖域映射相對誤差Table 2 Relative error of the conformal mapping in rock domain

圖5 取三點(diǎn)逼近夾角的示意圖Fig.5 Schematic diagram of approximating angle by three points

4.2 襯砌區(qū)域保角映射

選取Z平面上坐標(biāo)為(4.8，0)的點(diǎn)作為積分參考點(diǎn)，采用第3節(jié)中保角映射方法，將圖4中襯砌有限域映射為外半徑為1 m、內(nèi)半徑為0.912 2 m的同心圓環(huán)，求得保角映射函數(shù)，其中z的各冪級的系數(shù)見表3。在Z平面上襯砌內(nèi)外邊界上分別選取1 048 個(gè)點(diǎn)，經(jīng)映射函數(shù)映射得到ζ平面上共心圓環(huán)上的2 組點(diǎn)，其復(fù)數(shù)模的均值分別為1.000 0 和0.912 2，變異系數(shù)分別為6.86×10-6和1.16×10-6，說明此映射函數(shù)的保圓性較好。采用第3 節(jié)中保角映射方法求得的逆函數(shù)、保角映射函數(shù)ω1(ζ)，其中ζ各冪的系數(shù)見表3。

表3 圍巖域保角映射函數(shù)的系數(shù)Table 3 Coefficients of conformal mapping function in lining domain

在ζ平面上任意選取8 個(gè)夾角，如表4 所示，映射前后夾角的相對誤差量級等于或小于10-6，說明此映射函數(shù)的保角性良好。

表4 襯砌域映射相對誤差Table 4 Relative error of conformal mapping in lining domain

4.3 襯砌圍巖接觸面上邊界點(diǎn)保角映射差異

文獻(xiàn)[7-14]中對圍巖和襯砌共用同一保角映射函數(shù)，在接觸面邊界條件的等式中忽略了接觸面邊界點(diǎn)保角映射后的點(diǎn)位差異。由4.1 節(jié)和4.2 節(jié)可知，圍巖和襯砌的保角映射函數(shù)是不同的，本節(jié)對接觸面邊界上2個(gè)區(qū)域的點(diǎn)經(jīng)保角映射后的差值進(jìn)行計(jì)算。圖6所示為接觸面邊界點(diǎn)保角映射點(diǎn)位差異示意圖。在Z平面上沿襯砌與圍巖接觸面邊界選取若干點(diǎn)z2k(k=1，…，n)，經(jīng)襯砌和圍巖的保角映射函數(shù)ω*1(z)和ω*2(z)分別映射為ζ復(fù)平面上單位圓上的點(diǎn)ζ21k=eiβ1k和ζ22k=eiβ2k(k=1，…，n)。圖7所示為Z平面圍巖與襯砌接觸面邊界上同一點(diǎn)分別映射后的點(diǎn)位差異。圖7縱坐標(biāo)為對應(yīng)點(diǎn)ζ21k與ζ22k在ζ 平面上的點(diǎn)位差，即角度差值β1k-β2k(k=1，…，n)。角度差值隨圍巖點(diǎn)位角度變化波動(dòng)較大，襯砌和圍巖接觸面邊界上點(diǎn)分別保角映射后的點(diǎn)位差異不能經(jīng)旋轉(zhuǎn)、縮放和平移等線性變換消除。因此，文獻(xiàn)[7-14]在忽視點(diǎn)位映射后差異的情況下所建立的圍巖與襯砌接觸面邊界上應(yīng)力、位移邊界條件等式是不合適的；圍巖與襯砌接觸面上邊界條件等式應(yīng)在左右兩側(cè)分別代入ζ21k與ζ22k。

圖6 接觸面邊界點(diǎn)保角映射點(diǎn)位差異示意圖Fig.6 Illustration of position difference of points of two domains on contact boundary after conformal mapping

圖7 接觸面邊界上相鄰區(qū)域點(diǎn)映射后的點(diǎn)位差Fig.7 Position difference of points of two domains on contact surface boundary after conformal mapping

5 隧道力學(xué)解析計(jì)算實(shí)例

選取圖4 中隧道開挖輪廓和襯砌的幾何參數(shù)，令豎向地應(yīng)力為80 MPa，側(cè)壓力系數(shù)為0.50，襯砌的泊松比與彈性模量分別為0.2 和30 000 MPa，圍巖的泊松比與彈性模量分別為0.3和25 000MPa。根據(jù)以上保角映射函數(shù)和應(yīng)力函數(shù)求解方法，基于MATLAB 編制求解隧道力學(xué)響應(yīng)的計(jì)算程序。在實(shí)際計(jì)算中，幾乎不可能實(shí)現(xiàn)n2和m1的取值，使矩陣Z1的秩恰好等于應(yīng)力函數(shù)系數(shù)的個(gè)數(shù)(有且僅有唯一解)。針對此大型稀疏線性方程組，本文采用求解算法LSQR[31]對矩陣求最小二乘解。

分別在邊界L1、和上選取2n2個(gè)點(diǎn)位，將其分別代入所求應(yīng)力函數(shù)和對應(yīng)的邊界條件式，左式減右式得fb21(k)、fb22(k)和fb23(k)(k=1，2，…2n2)這3 組數(shù)據(jù)，數(shù)據(jù)fb21的均值分別記為M1、M2和M3，數(shù)據(jù)的變異系數(shù)分別記V1、V2和V3。當(dāng)n2和m1取值不同時(shí)，解析解的精確性和計(jì)算時(shí)長也有較大差異。表5所示為隧道邊界條件誤差分析。由表5 可見：當(dāng)n2和m1分別取120 與60、160 與80、200 與100 的組合，邊界上誤差的量級隨應(yīng)力函數(shù)項(xiàng)數(shù)增加而大幅降低至10-6及以下，所求應(yīng)力函數(shù)能較好地滿足各邊界條件。

表5 隧道邊界條件誤差分析Table 5 Error data of boundary condition of lined tunnel

為進(jìn)一步驗(yàn)證本文解析計(jì)算方法，在有限差分軟件FLAC 3D 中對此含襯砌隧道進(jìn)行數(shù)值模擬計(jì)算。建立隧道開挖的平面應(yīng)變模型，上邊界作用豎向地應(yīng)力，左右邊界作用水平地應(yīng)力，下邊界固定；模擬全斷面開挖，開挖后立即支護(hù)，圍巖及襯砌均為實(shí)體單元的彈性模型[32-33]?？紤]到網(wǎng)格劃分和邊界條件對計(jì)算結(jié)果的影響，比較不同網(wǎng)格劃分及邊界方案的計(jì)算結(jié)果。上下左右邊界距離約為20倍洞徑時(shí)網(wǎng)格劃分見圖8，此時(shí)計(jì)算結(jié)果已十分穩(wěn)定。因此，本文選取圖8所示網(wǎng)格劃分方案和20 倍洞徑的邊界距離。隧道數(shù)值解的位移云圖見圖9。洞形非圓形不規(guī)則，隧洞拱頂和拱底變形存在差異，分別為2.40 cm和2.09 cm；受非靜水壓力的影響，拱腰變形遠(yuǎn)小于拱頂變形，拱腰變形僅為0.90 cm。

圖8 FLAC 3D計(jì)算模型的網(wǎng)格劃分Fig.8 Meshing drawing of FLAC 3D computation model

圖9 隧道數(shù)值解的位移云圖Fig.9 Displacement nephogram of tunnel numerical solution

在圖2 中的Z平面上從x軸正方向起始，以逆時(shí)針方向?yàn)檎?，以順時(shí)針方向?yàn)樨?fù)，沿襯砌內(nèi)邊界和外邊界選取若干點(diǎn)位，提取各點(diǎn)處微小單元體位移、應(yīng)力的解析解和數(shù)值解。圖10 所示為襯砌內(nèi)外邊界位移的解析解和數(shù)值解對比。數(shù)值解受網(wǎng)格劃分、插值提取應(yīng)力等的影響，常偏離解析解。內(nèi)外邊界位移在拱頂和拱底附近的相對誤差較小，在±5%以內(nèi)；拱腳處形狀較為復(fù)雜，相對誤差較大，在±20%以內(nèi)。整體上，從拱底到拱頂襯砌變形量的解析解和數(shù)值解的變化趨勢大致相同。

圖10 襯砌邊界上位移的解析解與數(shù)值解對比Fig.10 Comparison of analytical and numerical solutions of displacement on lining boundaries

圖11 對比了襯砌內(nèi)邊界上環(huán)向應(yīng)力、圍巖和襯砌接觸面邊界上法向應(yīng)力和切應(yīng)力的解析解和數(shù)值解。環(huán)向和切向均指在ζ平面上極坐標(biāo)系下的方向。在接觸面邊界上，環(huán)向應(yīng)力、法向應(yīng)力和切應(yīng)力的解析解和數(shù)值解的變化趨勢基本相同。環(huán)向應(yīng)力和法向應(yīng)力的相對誤差大多在±10%以內(nèi)；切向應(yīng)力的相對誤差大多在±20%以內(nèi)，離散元數(shù)值軟件采用網(wǎng)格劃分和插值計(jì)算方法，使得個(gè)別位置數(shù)值解的相對誤差較大。

圖11 襯砌邊界上應(yīng)力的解析解和數(shù)值解對比Fig.11 Comparison of analytical and numerical solutions of stress on lining boundaries

在襯砌臨空面上，拱頂和拱底處所受壓應(yīng)力最小，在拱腳處出現(xiàn)應(yīng)力集中，所受壓應(yīng)力達(dá)到拱底處的7倍左右，是襯砌結(jié)構(gòu)最易損傷區(qū)。在圍巖和襯砌接觸面邊界上，法向應(yīng)力均為負(fù)(壓力)，圍巖和襯砌未出現(xiàn)脫空；拱底處法向應(yīng)力最小，拱腳處附近法向應(yīng)力最大，從拱頂?shù)焦澳_附近的法向應(yīng)力在12 MPa 左右。接觸面邊界上切向應(yīng)力有正有負(fù)，不同位置的襯砌和圍巖發(fā)生相對滑動(dòng)或相對滑動(dòng)趨勢的方向不同；在拱腳附近靠近拱底處，切向應(yīng)力也最大，在拱頂、拱底處切應(yīng)力為零。

在圖2的Z平面上沿x軸正方向、y軸正方向分別選取若干點(diǎn)，提取各點(diǎn)處微小單元體應(yīng)力的解析解和數(shù)值解。圖12 所示為圍巖內(nèi)部各點(diǎn)主應(yīng)力的數(shù)值解和解析解。在x軸正方向上，圍巖最小主應(yīng)力逐漸增大、圍巖最大主應(yīng)力逐漸減小，直至等于原巖應(yīng)力；圍巖最大主應(yīng)力的最大值達(dá)到2倍以上的原巖應(yīng)力。在y軸正方向上，最大主應(yīng)力、最小主應(yīng)力均在豎向或水平方向上；在洞周y=0 m處，最大主應(yīng)力為水平應(yīng)力，水平應(yīng)力先略微增加后逐漸減?。恢羪=8.75 m 處，水平應(yīng)力與豎向應(yīng)力相等，均為52.0 MPa；y大于8.75 m后，水平應(yīng)力持續(xù)減小至40 MPa，且在此處變換為最小主應(yīng)力；豎向應(yīng)力則持續(xù)增加直至80 MPa。

圖12 圍巖內(nèi)部應(yīng)力解析解和數(shù)值解對比Fig.12 Comparison of analytical and numerical solutions of stress in rock mass

沿x軸正方向洞周位置處的最大主應(yīng)力的數(shù)值解和解析解分別為168.91 MPa 和170.57 MPa，相對誤差為9.7‰；最小主應(yīng)力的數(shù)值解和解析解分別為2.21 MPa和2.18 MPa，相對誤差為1.4%。沿y軸正方向洞周位置處的最大主應(yīng)力的數(shù)值解和解析解分別為56.41 MPa 和56.01 MPa，相對誤差為7.1‰；最小主應(yīng)力的數(shù)值解和解析解為1.52 MPa和1.61 MPa，相對誤差為5.6%。隨著距離增加，網(wǎng)格劃分逐漸規(guī)律，應(yīng)力變化漸小，主應(yīng)力的解析解和數(shù)值解趨于完全一致。

6 結(jié)論

1) 考慮隧道圍巖和襯砌保角映射的差異，采用基于邊界積分方程的保角映射數(shù)值計(jì)算方法，編制了以多段圓弧為邊界的雙連通無限域、雙連通有限區(qū)域的保角映射程序，提出多點(diǎn)逼近的應(yīng)力函數(shù)求解方法。

2) 以多心圓隧道為例，分別計(jì)算得到圍巖和襯砌的保角映射函數(shù)，并證明了所求映射函數(shù)具有良好的保角性；驗(yàn)證了圍巖和襯砌在接觸面邊界上的點(diǎn)在保角映射后并不會(huì)重合，證明了既有非圓形隧道力學(xué)解析研究中對圍巖和襯砌采用同一保角映射函數(shù)的做法是不合適的。

3) 采用多點(diǎn)逼近的應(yīng)力函數(shù)求解方法求得滿足襯砌臨空面、襯砌與圍巖接觸面邊界上位移和應(yīng)力條件的圍巖和襯砌應(yīng)力函數(shù)；隨著應(yīng)力函數(shù)項(xiàng)和邊界點(diǎn)位數(shù)量增加，解析精度隨之增加。隧道應(yīng)力、位移的解析解和數(shù)值解十分接近，驗(yàn)證了此非圓形隧道彈性力學(xué)解析方法的準(zhǔn)確性。

4) 在側(cè)壓力系數(shù)為0.50 的條件下，該隧道洞周收斂變形的差異顯著，拱頂和拱底變形最大，拱腳附近變形最小；襯砌拱腳位置附近所受圍巖壓力和切向力最大，并出現(xiàn)顯著的應(yīng)力集中現(xiàn)象；在水平方向上，圍巖最大、最小主應(yīng)力差異隨深度增加而減小至穩(wěn)定；在豎向上，圍巖最大、最小主應(yīng)力差值隨深度增加而呈現(xiàn)先減小再增加至穩(wěn)定的現(xiàn)象。