考虑围岩和衬砌保角映射差异的深埋非圆形隧道力学解析

陈文博，房倩，张顶立，陈炫皓

(北京交通大学 城市地下工程教育部重点实验室，北京，100044)

非圆形隧道力学的解析方法能通过显式解快速计算岩体任意位置处的力学响应，克服数值软件计算结果受单元尺寸、网格划分、边界距离和插值计算影响的缺点，在提取应力迹线、分析隧道承载拱效应等研究上具有应用潜力和相对优势，是分析深部岩体内部复杂力学行为的有效手段，多年来一直是研究的难点和热点。基于复分析理论对非圆形隧道进行力学解析的关键在于：1) 求解衬砌、围岩等区域的保角映射函数；2) 求解各区域内的应力函数。

针对隧道力学解析中的保角映射，吕爱钟等[1]提出复合形法、混合罚函数法等优化方法；祝江鸿等[2]提出奇偶插值点相互迭代的计算方法；QI等[3]提出求解矩形隧道围岩保角映射函数的萤火虫寻优算法；HE 等[4]基于黎曼映射定理和边界对应定理，提出求解不规则隧道保角映射函数的迭代算法。这些研究未对该保角映射方法的复分析理论基础进行严谨的推导，认为映射前后区域的边界点相互对应(无限接近)即为保角映射，所求映射函数仅具有将一个边界映射为另一个边界的功能，不能保证区域内映射的保角性[5-6]。

在非圆形有支护隧道弹性力学解析的研究中，KARGAR 等[7-14]采用文献[1-3]中保角映射函数的求解方法，基于MUSKHELISHVILI等[15]的复分析理论，探索有衬砌支护、注浆圈加固的非圆形隧道力学解析，利用复合闭路定理和解析域内的幂级数性质，提出了“同幂次消元”的应力函数求解方法。文献[7-14]对衬砌域、加固圈域和围岩域采用同一个保角映射函数，使得相邻边界上不同区域的点在映射后仍重合，以便采用“消元法”求解应力函数。由保角映射的唯一性定理[6]可知，仅具有一条共用边界的有限双连通域(衬砌)和无限双连通区域(围岩)的保角映射函数是不同的(严格意义上讲，即不等价，2个保角映射函数不能通过旋转、缩放和平移等线性保角变换相互转化)。因此，文献[7-14]中求应力函数的方法是不适用的。

针对以上保角映射问题和应力函数求解问题，首先，引进了严格的基于边界积分的保角映射方法，并编制了适用于隧道的保角映射程序；其次，根据隧道的边界条件，提出了求解围岩和衬砌应力函数的多点逼近方法；再次，通过多心圆隧道的算例，验证了开发程序所求映射函数的保角性和保圆性，分析了围岩与衬砌接触面边界上点在保角映射后的位置差异；最后，对比隧道力学解析和数值模拟计算结果，分析了侧压力系数为0.50时隧道的力学响应规律，并验证了本文提出非圆形深埋隧道力学解析方法的正确性。

1 双连通无限域、有限域的保角映射方法

图1 所示为保角映射与一般映射差异说明图。从图1 可见：内边界为L1的双连通无限域S1经保角映射函数映射为内边界半径为R的圆形双连通无限域；外边界为L1，内边界分别为L21、L22和L23的双连通有限域S21、S22和S23分别经保角映射函数和映射为外径R和内径分别为r1、r2和r3的共心圆环。由单连通域及多连通域保角映射的唯一性定理[6,16]可知：互不等价，r1、r2和r3互不相等。映射函数和虽然同样可以将无限域S1的边界L1映射为半径为R的圆，但均不是无限域S1的保角映射函数，仅有是其保角映射函数。

图1 保角映射与一般映射差异说明图Fig.1 Illustration of difference between conformal mapping and general mapping

文献[1-3]所提保角映射方法存在问题的原因在于：

1) 曲解了黎曼存在性定理，将单连通域保角映射的边界对应定理错误地应用到双连通无限域、有限域的保角映射上，单纯地将边界之间的映射作为两区域之间的保角映射。

2) 其优化方法、迭代计算方法仅可实现边界之间的近似映射，不能保证区域映射的保角性和映射点之间的一一对应关系。

因此，本文进一步基于保角映射的复分析理论和数值计算方法，研究非圆形隧道衬砌域、围岩域映射为共心圆环域、含圆无限域的保角映射。关于保角映射的复分析理论，CROWDY 等[17-20]基于Schwarz-Christoffel 公式的经典函数论，引入Schottky-Klein 素函数，针对以圆为边界的多连通域到多边形域、狭缝域、多圆弧域的保角映射，建立了高阶微分方程；NASSER[21]通过边界积分的数值方法求解此微分方程，进一步将保角映射问题转化为边界积分方程的求解问题。在数值计算方法方面，NASSER 等[22]提出了一种边界积分方法解决Laplace 方程在多连通区域上的Dirichlet 边界条件问题和Neumann 边界条件问题，建立了第二、三、四、五类Koebe标准狭缝域保角映射的数值计算方法[23-24]，提出了含狭缝的多连通域的边界积分方法[25]。NASSER 等[26-27]基于Nyström 方法对边界积分方程进行离散化，提出了一种求解含广义Neumann 核函数和共轭广义Neumann 核函数边界积分方程的数值计算方法。

基于文献[20]和[27]，本文采用MATLAB语言编写程序，将以多段圆弧为边界的多连通域映射为以圆为边界的多连通域。对于如图2所示的有限双连通域，首先利用基于边界积分方程的保角映射方法，将复平面Z上的有限双连通域边界上各点保角映射为复平面η上偏心圆环上的点；再利用Möbius 变换(分式线性变换)[16,28]，得到在ζ平面半径分别为R和r的共心圆环上相对应的各点，映射前后边界上各点一一对应。

图2 不同复平面上衬砌边界Fig.2 Lining boundary on different complex planes

为求得如式(1)Laurent 级数形式的保角映射函数，设待求映射函数的项数为2m1+1，在Z平面和ζ平面的边界上选取n1对点，将点坐标zt和ζt(t=1，2，…，n1)代入式(1)，建立求常系数x1的矩阵方程，见式(3)。调整m1和n1，即可获得满足精度要求的保角映射函数ω*(zt)。ω*(zt)的逆函数ω(ζt)也为保角映射函数(式(2))，将ζ平面圆环域上点映射到Z平面上衬砌区域上的点，其求法与ω*(zt)相同。

式中：ak和bk为复常系数。

围岩区域的保角映射函数求解原理和方法与衬砌相同，但不需要图2(b)的分式线性变换。

2 深埋隧道衬砌与围岩的边界条件

图3 所示为衬砌边界到同心圆环的映射示意图。如图3所示，L1和分别为衬砌内边界和外边界，为围岩内边界，围岩外边界为无限远处。保角映射函数ω1和ω2分别将ζ平面上的同心圆环域和含圆无限域映射为Z平面上的衬砌区域S1和围岩区域S2。在边界L1上，z11=ω1(ζ11)；在边界上，z21=ω1(ζ21)；在边界上，z22=ω2(ζ22)；其中，z11、z21和z22分别为Z平面上衬砌内边界、外边界和围岩边界上的点，ζ11、ζ21和ζ22分别为ζ平面上衬砌内边界、外边界和围岩内边界上的点。

图3 衬砌边界到同心圆环的映射示意图Fig.3 Schematic diagram of mapping of lining boundaries to concentric rings

在ζ平面极坐标系下，应力和位移可由应力函数φ和ψ表示[29]。

式中：σθ、σρ和τθρ分别为环向应力、轴向应力和切应力；Re 表示对复数取实部；uρ和uθ分别为轴向位移和环向位移；G为材料的切变模量，由材料弹性模量E和泊松比μ决定，；κ在平面应变条件下取3-4μ。

对于深埋隧道的平面应变弹性力学求解，在ζ平面上，衬砌区域和围岩区域内的应力函数可表示为[7]

1) 区域S1内

2) 区域S2内

式中：ck、dk、ek、fk、hk、lk、mk和qk为应力函数的常复系数；hk和mk由无穷远处地应力确定。

在边界L1上，对于任意ζ11应满足σρ=0和τθρ=0 的应力条件，将ζ11代入式(4)和(5)，可得L1上边界条件。同时，在围岩内趋于无穷远处应力函数应满足式(11)和式(12)。

通过限定接触面边界上的应力、位移连续条件及其范围，可对衬砌与围岩完全接触、非完全解接触或局部脱离接触的隧道进行弹性力学解析计算。本文假定围岩和衬砌完全接触，既不相互脱离又不相互滑动[30]。在Z平面上，围岩和衬砌接触面上的法向应力、切应力和位移相同。由保角变换的性质可知，在ζ平面极坐标系下，边界和上对任意满足ω1(ζ21) =ω2(ζ22)的点ζ21和点ζ22，两点上的轴向应力、切应力和位移也相等。将ζ21和ζ22分别代入式(5)和(6)，即可得接触面边界上的2个边界条件。本研究基于复变理论进一步拓展非圆性非静水压力下深埋隧道的弹性力学解析方法，暂不能考虑围岩弹塑性、应变软化等问题。

3 应力函数的求解方法

在Z平面边界L1和L2上分别选取n2个点，经映射得到ζ平面上的3 组点ζ11k=ω1(z1k)，ζ21k=ω1(z2k)，ζ22k=ω2(z2k)，k=1，…，n2。将ζ平面上3n2个点分别代入对应的边界条件，可提取得到应力函数系数6组复数未知量ck、dk、ek、fk、hk和qk(k=1，…，m1)的3n2个方程。为能够使用MATLAB的矩阵运算进行求解，将复数未知量拆分为以实数表示的形式，如ck=ck1+ick2(i 为虚数单位)，并代入方程组。

基于MATLAB 编写程序，逐个搜索各等式方程中的常数项和各实未知量(ck1、ck2等)的系数，将各实未知量对应的系数写成6n2行、8+12m1列的矩阵Z2，对应的实常数项依次序写成6n2行的列向量。

式中：x2=[c01，c02，…，q01，q02，c11，c12，…，q11，q12，…cm11，cm12，…，qm11，qm12]T；Φ2为各方程中常数项相反数的列向量。调整n2和m1，并选取合适方法求解矩阵，即可获得满足边界条件的应力函数系数。

4 隧道的保角映射及验证

选取某铁路单线隧道衬砌断面，其几何尺寸见图4，在以O为原点、x轴为实轴、y为虚轴的Z平面上，分别对以该断面衬砌外轮廓为内边界的围岩无限域、该断面衬砌的双连通有限域进行保角映射。

图4 多心圆隧道断面图Fig.4 Multi-center circle tunnel cross-section

4.1 围岩区域保角映射

采用第3节保角映射方法，将图4中围岩无限域映射为含单圆孔无限域(圆孔半径为1 m)，求得保角映射函数，其中z的1 次幂至-6 次幂的系数见表1。原求得保角映射函数的级数共17 项(小于10-4的系数未列出，表中以缺省符号表示，下同)。在Z平面上围岩内边界上选取3 072 个点，经过映射函数映射得到ζ平面上3 072 个点，各点复数模的均值为1.000 0，变异系数为2.10×10-6，说明该映射函数的保圆性较好。采用第3节中保角映射方法，求得的逆函数、保角映射函数ω2(ζ)，其中ζ的1次幂至-8次幂的系数见表1。

表1 围岩域保角映射函数的系数Table 1 Coefficients of conformal mapping function in rock domain

由第3 节可知，解析计算中仅需要ω1(ζ)和ω2(ζ)，因此，本文仅检验该映射函数的保角性。图5所示为在ζ平面上取任意夹角的示意图，为得到各角映射后的角度，令ds足够小(本文取10-4m)，在夹角直线上各取2个点，将点a、点b和点c分别映射到Z平面，映射点连线的夹角可以近似认为是映射后的曲线夹角。在ζ平面上任意选取8 个夹角，如表2所示，映射前后夹角的相对误差量级等于或小于10-6，说明映射函数的保角性很好。

表2 围岩域映射相对误差Table 2 Relative error of the conformal mapping in rock domain

图5 取三点逼近夹角的示意图Fig.5 Schematic diagram of approximating angle by three points

4.2 衬砌区域保角映射

选取Z平面上坐标为(4.8，0)的点作为积分参考点，采用第3节中保角映射方法，将图4中衬砌有限域映射为外半径为1 m、内半径为0.912 2 m的同心圆环，求得保角映射函数，其中z的各幂级的系数见表3。在Z平面上衬砌内外边界上分别选取1 048 个点，经映射函数映射得到ζ平面上共心圆环上的2 组点，其复数模的均值分别为1.000 0 和0.912 2，变异系数分别为6.86×10-6和1.16×10-6，说明此映射函数的保圆性较好。采用第3 节中保角映射方法求得的逆函数、保角映射函数ω1(ζ)，其中ζ各幂的系数见表3。

表3 围岩域保角映射函数的系数Table 3 Coefficients of conformal mapping function in lining domain

在ζ平面上任意选取8 个夹角，如表4 所示，映射前后夹角的相对误差量级等于或小于10-6，说明此映射函数的保角性良好。

表4 衬砌域映射相对误差Table 4 Relative error of conformal mapping in lining domain

4.3 衬砌围岩接触面上边界点保角映射差异

文献[7-14]中对围岩和衬砌共用同一保角映射函数，在接触面边界条件的等式中忽略了接触面边界点保角映射后的点位差异。由4.1 节和4.2 节可知，围岩和衬砌的保角映射函数是不同的，本节对接触面边界上2个区域的点经保角映射后的差值进行计算。图6所示为接触面边界点保角映射点位差异示意图。在Z平面上沿衬砌与围岩接触面边界选取若干点z2k(k=1，…，n)，经衬砌和围岩的保角映射函数ω*1(z)和ω*2(z)分别映射为ζ复平面上单位圆上的点ζ21k=eiβ1k和ζ22k=eiβ2k(k=1，…，n)。图7所示为Z平面围岩与衬砌接触面边界上同一点分别映射后的点位差异。图7纵坐标为对应点ζ21k与ζ22k在ζ 平面上的点位差，即角度差值β1k-β2k(k=1，…，n)。角度差值随围岩点位角度变化波动较大，衬砌和围岩接触面边界上点分别保角映射后的点位差异不能经旋转、缩放和平移等线性变换消除。因此，文献[7-14]在忽视点位映射后差异的情况下所建立的围岩与衬砌接触面边界上应力、位移边界条件等式是不合适的；围岩与衬砌接触面上边界条件等式应在左右两侧分别代入ζ21k与ζ22k。

图6 接触面边界点保角映射点位差异示意图Fig.6 Illustration of position difference of points of two domains on contact boundary after conformal mapping

图7 接触面边界上相邻区域点映射后的点位差Fig.7 Position difference of points of two domains on contact surface boundary after conformal mapping

5 隧道力学解析计算实例

选取图4 中隧道开挖轮廓和衬砌的几何参数，令竖向地应力为80 MPa，侧压力系数为0.50，衬砌的泊松比与弹性模量分别为0.2 和30 000 MPa，围岩的泊松比与弹性模量分别为0.3和25 000MPa。根据以上保角映射函数和应力函数求解方法，基于MATLAB 编制求解隧道力学响应的计算程序。在实际计算中，几乎不可能实现n2和m1的取值，使矩阵Z1的秩恰好等于应力函数系数的个数(有且仅有唯一解)。针对此大型稀疏线性方程组，本文采用求解算法LSQR[31]对矩阵求最小二乘解。

分别在边界L1、和上选取2n2个点位，将其分别代入所求应力函数和对应的边界条件式，左式减右式得fb21(k)、fb22(k)和fb23(k)(k=1，2，…2n2)这3 组数据，数据fb21的均值分别记为M1、M2和M3，数据的变异系数分别记V1、V2和V3。当n2和m1取值不同时，解析解的精确性和计算时长也有较大差异。表5所示为隧道边界条件误差分析。由表5 可见：当n2和m1分别取120 与60、160 与80、200 与100 的组合，边界上误差的量级随应力函数项数增加而大幅降低至10-6及以下，所求应力函数能较好地满足各边界条件。

表5 隧道边界条件误差分析Table 5 Error data of boundary condition of lined tunnel

为进一步验证本文解析计算方法，在有限差分软件FLAC 3D 中对此含衬砌隧道进行数值模拟计算。建立隧道开挖的平面应变模型，上边界作用竖向地应力，左右边界作用水平地应力，下边界固定；模拟全断面开挖，开挖后立即支护，围岩及衬砌均为实体单元的弹性模型[32-33]。考虑到网格划分和边界条件对计算结果的影响，比较不同网格划分及边界方案的计算结果。上下左右边界距离约为20倍洞径时网格划分见图8，此时计算结果已十分稳定。因此，本文选取图8所示网格划分方案和20 倍洞径的边界距离。隧道数值解的位移云图见图9。洞形非圆形不规则，隧洞拱顶和拱底变形存在差异，分别为2.40 cm和2.09 cm；受非静水压力的影响，拱腰变形远小于拱顶变形，拱腰变形仅为0.90 cm。

图8 FLAC 3D计算模型的网格划分Fig.8 Meshing drawing of FLAC 3D computation model

图9 隧道数值解的位移云图Fig.9 Displacement nephogram of tunnel numerical solution

在图2 中的Z平面上从x轴正方向起始，以逆时针方向为正，以顺时针方向为负，沿衬砌内边界和外边界选取若干点位，提取各点处微小单元体位移、应力的解析解和数值解。图10 所示为衬砌内外边界位移的解析解和数值解对比。数值解受网格划分、插值提取应力等的影响，常偏离解析解。内外边界位移在拱顶和拱底附近的相对误差较小，在±5%以内；拱脚处形状较为复杂，相对误差较大，在±20%以内。整体上，从拱底到拱顶衬砌变形量的解析解和数值解的变化趋势大致相同。

图10 衬砌边界上位移的解析解与数值解对比Fig.10 Comparison of analytical and numerical solutions of displacement on lining boundaries

图11 对比了衬砌内边界上环向应力、围岩和衬砌接触面边界上法向应力和切应力的解析解和数值解。环向和切向均指在ζ平面上极坐标系下的方向。在接触面边界上，环向应力、法向应力和切应力的解析解和数值解的变化趋势基本相同。环向应力和法向应力的相对误差大多在±10%以内；切向应力的相对误差大多在±20%以内，离散元数值软件采用网格划分和插值计算方法，使得个别位置数值解的相对误差较大。

图11 衬砌边界上应力的解析解和数值解对比Fig.11 Comparison of analytical and numerical solutions of stress on lining boundaries

在衬砌临空面上，拱顶和拱底处所受压应力最小，在拱脚处出现应力集中，所受压应力达到拱底处的7倍左右，是衬砌结构最易损伤区。在围岩和衬砌接触面边界上，法向应力均为负(压力)，围岩和衬砌未出现脱空；拱底处法向应力最小，拱脚处附近法向应力最大，从拱顶到拱脚附近的法向应力在12 MPa 左右。接触面边界上切向应力有正有负，不同位置的衬砌和围岩发生相对滑动或相对滑动趋势的方向不同；在拱脚附近靠近拱底处，切向应力也最大，在拱顶、拱底处切应力为零。

在图2的Z平面上沿x轴正方向、y轴正方向分别选取若干点，提取各点处微小单元体应力的解析解和数值解。图12 所示为围岩内部各点主应力的数值解和解析解。在x轴正方向上，围岩最小主应力逐渐增大、围岩最大主应力逐渐减小，直至等于原岩应力；围岩最大主应力的最大值达到2倍以上的原岩应力。在y轴正方向上，最大主应力、最小主应力均在竖向或水平方向上；在洞周y=0 m处，最大主应力为水平应力，水平应力先略微增加后逐渐减小；至y=8.75 m 处，水平应力与竖向应力相等，均为52.0 MPa；y大于8.75 m后，水平应力持续减小至40 MPa，且在此处变换为最小主应力；竖向应力则持续增加直至80 MPa。

图12 围岩内部应力解析解和数值解对比Fig.12 Comparison of analytical and numerical solutions of stress in rock mass

沿x轴正方向洞周位置处的最大主应力的数值解和解析解分别为168.91 MPa 和170.57 MPa，相对误差为9.7‰；最小主应力的数值解和解析解分别为2.21 MPa和2.18 MPa，相对误差为1.4%。沿y轴正方向洞周位置处的最大主应力的数值解和解析解分别为56.41 MPa 和56.01 MPa，相对误差为7.1‰；最小主应力的数值解和解析解为1.52 MPa和1.61 MPa，相对误差为5.6%。随着距离增加，网格划分逐渐规律，应力变化渐小，主应力的解析解和数值解趋于完全一致。

6 结论

1) 考虑隧道围岩和衬砌保角映射的差异，采用基于边界积分方程的保角映射数值计算方法，编制了以多段圆弧为边界的双连通无限域、双连通有限区域的保角映射程序，提出多点逼近的应力函数求解方法。

2) 以多心圆隧道为例，分别计算得到围岩和衬砌的保角映射函数，并证明了所求映射函数具有良好的保角性；验证了围岩和衬砌在接触面边界上的点在保角映射后并不会重合，证明了既有非圆形隧道力学解析研究中对围岩和衬砌采用同一保角映射函数的做法是不合适的。

3) 采用多点逼近的应力函数求解方法求得满足衬砌临空面、衬砌与围岩接触面边界上位移和应力条件的围岩和衬砌应力函数；随着应力函数项和边界点位数量增加，解析精度随之增加。隧道应力、位移的解析解和数值解十分接近，验证了此非圆形隧道弹性力学解析方法的准确性。

4) 在侧压力系数为0.50 的条件下，该隧道洞周收敛变形的差异显著，拱顶和拱底变形最大，拱脚附近变形最小；衬砌拱脚位置附近所受围岩压力和切向力最大，并出现显著的应力集中现象；在水平方向上，围岩最大、最小主应力差异随深度增加而减小至稳定；在竖向上，围岩最大、最小主应力差值随深度增加而呈现先减小再增加至稳定的现象。