深埋马蹄形隧道开挖围岩应力与位移的复变函数解

施有志，高轩能

（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.