弹性半无限空间中矩形孔收缩的复变函数解答

申航 周航 刘汉龙

摘 要：隨着城市建设的发展，矩形隧道的应用越来越多，但针对矩形隧道的理论研究却鲜有见闻。针对矩形隧道，建立了半无限空间矩形隧道的弹性理论计算模型，采用最小二乘迭代方法确定共形映射函数的各项系数，并将计算区域映射为复平面上的一个同心圆环;运用Muskhelishvili复变函数方法，将计算区域内的应力函数展开成为Laurant级数的形式，给定了地表零应力边界和矩形孔口径向位移边界，求得了半无限空间矩形隧道在给定位移条件下的应力场和位移场。分析了不同高宽比、不同泊松比、不同埋深对位移场和应力场的影响，总结了矩形隧道位移场和应力场的一般规律。结果表明：高宽比偏小、泊松比偏大、埋深偏小都会使得沉降槽不再是类高斯曲线的形状，这些参数的变化也会在不同程度上影响应力场和位移场的大小和分布。

关键词：孔收缩;半无限空间;共形映射;位移条件;应力场;位移场

中图分类号：TU431   文献标志码：A   文章编号：2096-6717（2020）04-0001-11

收稿日期：2020-11-27

基金项目：国家自然科学基金（51978105）;重庆市技术创新与应用发展专项（cstc2019jscx-msxmX0107）

作者简介：申航（1996- ），男，主要从事隧道开挖对桩基的影响研究，E-mail：457570774@qq.com。

周航（通信作者），男，副教授，博士生导师，E-mail：zh4412517@163.com。

Received：2020-11-27

Foundation items：National Natural Science Foundation of China （No. 51978105）; Chongqing Technical Innovation and Application Development Special Project （No. cstc2019jscx-msxmX0107）

Author brief：SHEN Hang （1996- ）， main research interest： influence of tunnel excavation on pile foundation， E-mail： 457570774@qq.com.

ZHOU Hang （corresponding author）， associate professor， doctorial supervisor， E-mail： zh4412517@163.com.

Solution of complex function of rectangular hole contraction in elastic semi-infinite space

SHEN Hang， ZHOU Hang， LIU Hanlong

（Key Laboratory of New Technology for Construction of Cities in Mountain Area， Ministry of Education; School of Civil Engineering， Chongqing University， Chongqing 400045， P. R. China）

Abstract： With the development of urban construction in China， more and more applications of rectangular tunnel are emerging， but there are few theoretical studies on rectangular tunnels. In this paper， the elastic theoretical calculation model of the semi-infinite space rectangular tunnel was established. The coefficients of conformal mapping function were determined by the least squares iterative method. And the calculated area was mapped to a concentric ring on the complex plane. Afterwards， The Muskhelishvili complex function method was used to expand the stress function in the calculation area into the form of Laurant series， which gives the zero stress boundary on the ground surface and the radial displacement boundary of the rectangular hole. The stress field and displacement field of the rectangular tunnel in the semi-infinite space under the given displacement condition were also obtained by the method. In this paper， the influence of aspect ratios， Poisson's ratios， and buried depths on the displacement field and stress field was analyzed， and the general rules of the displacement field and stress field of rectangular tunnels has been summarized. The results show that a smaller aspect ratio， a larger Poisson's ratio， and a smaller buried depth will make the settlement curve no longer similar to a Gaussian curve. The variation of these parameters will also affect the size and distribution of the stress field and displacement field to varying degrees.

Keywords：constraction of hole; semi-infinite space; conformal mapping; displacement condition; stress field; displacement field

近年來，中国城市隧道建设突飞猛进，不仅数量上上了新台阶，而且越来越多的隧道类型不断涌现。其中，矩形顶管隧道相较传统圆形隧道，有着埋深浅、断面面积利用率高、施工时对地面交通影响小、无污染、无噪音等优点。近年来，在城市地下过街通道、地下综合管廊、地铁隧道、中短距离的城市地下道路、地下空间的互联互通、地铁暗挖车站等工程中，矩形顶管隧道都得到了广泛运用。矩形顶管隧道的发展，代表着中国未来城市中短程隧道建设的新方向。2002年，日本京都地铁工程采用了矩形盾构机，首次成功建成了矩形单洞双线隧道[1];1995年，中国的矩形隧道开始起步，2015年，宁波轨道交通3号线采用类矩形盾构，并圆满完工，标志着中国在类矩形盾构技术方面取得重大突破并处于世界领先行列。2004年，波士顿中央大道矩形顶管隧道工程顺利完成，这项工程连接了美国两条主要洲际公路干线，工程进行中没有影响道路的正常运营。中国也有许多矩形顶管的应用，比如武汉地铁2号线王家墩东站4号出入口工程，通道总长62.4 m，仅用22 d便施工完成[2]。

虽然当前矩形隧道的工程实践较为丰富，但对应的浅埋矩形隧道理论研究工作却明显薄弱。目前，浅埋隧道的理论研究以圆形隧道为主。Jeffery[3]和Mindlin[4]采用双极坐标法研究了圆形隧道洞周的应力分布，但由于双极坐标法的局限性，无法求出复杂隧道的应力场，也不能得出位移场的分布;而且双极坐标法在隧道埋深较浅时会产生较大的计算误差[5]。Sagasete[6]采用镜像法，假设土体为不可压缩的弹性半无限体，得到了土体的位移场和应力场，随后Verruijt等[7]在Sagasete研究的基础上进一步分析了土体位移场随隧道椭圆化变形的影响，但Sagasete的方法基于土体不可压缩的假设，只能计算泊松比为0.5时的情况，有其局限性。曾彬等[8-9]、魏刚等[10]基于随机介质理论，分析了双圆盾构隧道和类矩形隧道的土体位移规律，然而受制于随机介质理论的特点，只能得到位移场的结果，难以明确围岩的应力分布规律。

利用弹性力学复变函数方法，可以通过保角映射将一些较为复杂的几何单连通域映射为复平面内的简单几何单连通域（单位圆或同心圆环等），然后，在这些域内，可以结合给定的边界条件，较为方便地解出复应力函数。许多学者基于复变函数理论的这一特征，对圆形隧道以及其他异型隧道进行了很多有意义的研究。Verruijt[11-12]通过保角映射将半无限空间上的圆形孔洞映射为同心圆环，利用弹性力学方法，分别得到了不考虑体力情况下的给定孔边均布径向位移和给定孔边均布径向压力的两种解析解。随后，Strack等[13]考虑了体力的影响，在应力解析函数中添加了可以计算由于土体自重产生的开挖引起的不平衡力系的项，得到了孔边位移边值问题的解。路文超[14]、蔚立元等[15]、张永兴等[16]在此基础上考虑了不同荷载作用的情况，并作出了解答。王立忠等[17]基于Verruijt解法分析了Pack[18]提出的4种不同的圆形隧洞边界位移边界条件，并将其解答与隧道实测数据进行了比对，指出第3和第4边界条件更加符合实际工程结果。此前基于复变函数方法的研究多针对圆形隧道，而对于非圆孔隧道进行的研究较少。曾癸森[19]提出了新的保角映射函数，可以将半无限空间上的任意异型孔映射成同心圆环，并给出了椭圆形隧道的解答，验证了其可靠性。

笔者基于曾癸森[19]给出的保角映射函数，将半无限空间上的矩形隧道映射成复平面上的同心圆环，假设问题为平面应变问题，忽略体力的影响，假设地表边界为零应力边界，孔口边界为均匀径向收缩的矩形位移边界，通过弹性力学复变函数方法得到问题的解析解。通过与有限元结果对比验证解析解的可靠性，并基于得出的解答进一步分析矩形浅埋隧道开挖的应力场和位移场，总结矩形浅埋隧道位移场和应力场的一般规律，为之后的矩形隧道研究打下基础。

1 半无限空间矩形浅埋隧洞问题描述

图1为z平面（直角坐标系）下，半无限空间中矩形孔洞的示意图。b为矩形隧道的宽，h为隧道的高，d为孔洞中心到地表水平线的距离，A（O）点为坐标原点，B点为地表水平线的无穷远点，C点和D点分别为矩形隧道顶部和底部轮廓线与隧道中轴线的交点，E点和F点为矩形孔洞在右半平面的两个角点。θ为矩形孔洞的角点F相对矩形中心点所在的角度。R域为半无限空间中矩形孔洞以外的区域。假设土体为均质各向同性弹性材料，将问题简化成平面应变问题进行讨论。

根据Muskhelishvili[20]给出的平面应变问题的复变函数方法，弹性解答可以由R域内的两个解析函数φ（Z）和ψ（Z）表示出来。其中，应力分量为

σx+σy=2φ′（Z）+φ′（Z） （1）

σy-σx+2iτxy=2z-φ″（Z）+ψ′（Z） （2）

位移分量为

2G（u+iv）=κφ（Z）-zφ′（Z）-ψ（Z） （3）

式中：G=E/（2+2μ）是剪切模量;E为土体的弹性模量;κ是与泊松比μ相关联的参数，在该平面应变问题中，取κ=3-4μ。

在z平面中，地表平面无应力作用，故采用应力边界条件控制;矩形孔周采用均匀收缩的位移边界条件，故采用位移边界条件控制。边界条件为

z=z-：φ（Z）+z·φ′（Z）+ψ（Z）=0 （4）

z∈l：κ·φ（Z）-z·φ′（Z）-ψ（Z）=2G（u+iv）（5）

式中：l为矩形孔周线上的点形成的集合。

2 复变函数解法

2.1 共形映射

采用曾癸森[19]提出的共形映射公式（6），将z平面上的R區域映射为ζ平面上的同心圆环，其中，圆环内外半径分别为α和1，圆环内边界对应矩形孔洞边界，圆环外边界对应地表边界。如图2所示，z平面上的A、B、C、D、E、F点分别对应ζ平面上的A′、B′、C′、D′、E′、F′点。1和2分别为E′点和F′点在ζ平面上对应的幅角。

z=ω（ζ）=ia1+ζ1-ζ+i∑nk=1βk（ζk-ζ-k） （6）

式中：

z=x+iy=reiθ （7）

ζ=ξ+iη=ρei （8）

式中：a和β为待定系数，以及圆环域的内半径α也为待定系数。曾癸森[19]只计算了椭圆形隧洞，由于孔洞形状比较简单，所取的映射函数项较少，很容易求解，但矩形隧洞需要的映射函数项较多，大大提升了需要求解的非线性方程组的未知数数量，直接求解存在困难。Zhou等[21]曾采用基于最小二乘法迭代计算方法来确定任意空腔保角映射方程的常系数，参考其方法进行迭代计算，并作出一些改进。

1）首先，沿着矩形孔洞边界从D点开始按逆时针方向等间距取m+1个取样点（其中，第一个和最后一个取样点是同一个点），这样，就将孔洞轮廓均等地分成了m段，得到了m+1个取样点在z平面上的坐标。同时，在ζ平面上将圆环域的内边界从D′点开始，沿逆时针方向，同样等间距地划分成m段，这样，也就得到了ζ平面上的m+1个取样点的坐标。将z平面和ζ平面上得到的取样点，按照从起点处逆时针依次经过的先后顺序进行编号，并一一对应起来。在圆环域内边界取点的过程中，由于圆环域的内半径α暂且未知，不妨先假设一个初值（α0=0.5），再在之后的计算中迭代逼近真实值。

2）将m个取样点的坐标代入共形映射公式（6）中，将会得到一个由m+1个方程，k+1个未知数组成的方程组，这个方程组可以写成矩阵形式

AX=B （9）

式中：

A=i·1+ζ01-ζ0ζ0-ζ0-1ζ02-ζ0-2…ζ0k-ζ0-k1+ζ11-ζ1ζ1-ζ1-1ζ12-ζ1-2…ζ1k-ζ1-k1+ζ21-ζ2ζ2-ζ2-1ζ22-ζ2-2…ζ2k-ζ2-k              1+ζm1-ζmζm-ζm-1ζm2-ζm-2…ζmk-ζm-k（10）

X=[aβ1β2…βk]T （11）

B=[z0z1z2…zm]T （12）

根据最小二乘法原理，式（9）可以写成

X=（ATA）-1ATB （13）

3）通过过程1）中得到的z平面和ζ平面中取样点的坐标，可以得到一个初始的A和B矩阵，代入式（13）中便可解得第1组系数矩阵X。但得到的第1组系数是不准确的，需要进行进一步的迭代，以提高精度。

4）将上一步过程中得到的X矩阵再代入式（9）中，得到一组新的z平面上的参考点坐标。这组坐标代表了采用当前系数矩阵X的情况下，ζ平面中圆环内边界上的点投射到z平面上的情况。在仅进行了第一次迭代的情况下，会发现新的参考点偏离所需要的轮廓线较多。

5）为了进行下一次迭代，需要将新得到的参考点坐标进行修正。首先，需要计算当前参考点轮廓线的周长，将m个参考点中相邻的两两参考点间的距离都叠加起来便可得到;随后，计算从起点处沿逆时针方向到各个参考点所走过的路程;最后，根据各个参考点对应路程在周长中的占比，找到在精确的矩形轮廓上同样占比的取样点，由此便可得到一组修正后的取样点坐标。

6）由上一步可得到修正后的B矩阵，随后重复3）、4）、5）的迭代过程，直到第t次迭代产生的系数矩阵Xt和第t-1次迭代产生的系数矩阵Xt -1之间的误差满足精度要求，则认为迭代收敛。

7）以上过程是针对α0=0.5进行的迭代，此时在5）中得到的迭代收敛后的轮廓线还不是需要的矩形孔洞边界，需要进一步对α值进行修正，才能得到满足需要的α值。经过试算发现，当α值偏大时，轮廓线将会位于矩形孔洞边界的外侧，且α值越大，外扩的现象越明显;反之，当α值偏小时，轮廓线将会位于矩形孔洞边界的内侧，且α值越小，内缩的现象越明显。因此，通过计算判断轮廓线周长与精确矩形孔洞边界周长的大小关系，便可确定当前设定的α值与精确值之间的大小关系。当轮廓线周长大于孔洞周长时，说明当前α值偏大，将α值适当缩小，反之，则将α值适当放大，然后重新进行1）至6）的迭代过程。当得到的新轮廓线周长与孔洞周长之间的误差满足精度要求时，停止迭代过程，此时得到的待定系数a、βk、α能够精度较高地完成共形映射。

2.2 洞周位移收敛条件

隧道在开挖的过程中不可避免地会产生土体损失，一般来说，隧道的开挖面相比最终的隧道断面要略大一些。在矩形隧道开挖过程中，由于土体损失的存在，将会引发矩形孔收缩的问题，通过给定孔收缩的位移条件，对此问题进行分析。对于圆形隧道的位移边界条件，Park[18]给出过简化的4种径向位移边界条件，但矩形隧道的位移边界条件更加复杂，目前鲜有研究。采用如图3所示的矩形径向位移收敛模式，设定边界条件位移参数u0，图3中孔口位移值ub=b·u0、uh=h·u0，位移模式表述为

u+iv=[-b-2i（y+d）]·u0，∈（1，2）

（-2x-ih）·u0，∈（2，3）

[b-2i（y+d）]·u0，∈（3，4）

（-2x+ih）·u0，∈（0，1）∪（4，2π）（14）

式中：1、2、3、4为矩形4个角点在ζ平面中对应点的幅角。1和2已在图2中标出，根据对称性有3=2π-2、4=2π-1。

将ζ=αei=α（cos +isin ）代入共形映射函数式（6），使其展开成实部和虚部的形式。

z=ω（ζ）=x+iy （15）

x=2a·αsin 1+α2-2αcos -∑nk=1βk（αk+α-k）sin k （16）

y=a（1-α2）1+α2-2αcos +∑nk=1βk（αk-α-k）cos k （17）

图2中，通过矩形边界角点F（x，y）在z平面中对应的角度θ，可以推出

tan θ=hb=y+dx （18）

将式（16）、式（17）代入式（18）中，可以计算出F′在ζ平面对应的幅角1，同理可得E′对应的幅角2。

至此，矩形位移收敛模式分段函数的分段区间已经完全确定，可将其用傅里叶级数表示，方便后续计算。

u+iv=f（ασ）=∑+SymboleB@k=-SymboleB@Akσk（19）

2.3 应力函数的求解

由于φ（Z）和ψ（Z）是R域上的解析函数，映射函数ω（ζ）也是解析函数，因此，φ（Z）和ψ（Z）可以用ζ来表示。

φ（Z）=φ（ω（ζ））=φ（ζ）（20）

ψ（Z）=ψ（ω（ζ））=ψ（ζ）（21）

根据复变函数理论可以得出，在ζ平面的环形域γ上，φ（Z）和ψ（Z）可以展开成Laurent级数的形式。

φ（ζ）=∑+SymboleB@k=-SymboleB@akζk（22）

ψ（ζ）=∑+SymboleB@k=-SymboleB@bkζk（23）

同样，边界条件在进行一些微分代换后，也可以写成ζ表示的形式。

ζ=1∶φ（ζ）+ω（ζ）ω′（ζ）·φ′（ζ）+ψ（ζ）=0（24）

ζ=α∶κ·φ（ζ）-ω（ζ）ω′（ζ）·φ′（ζ）-ψ（ζ）=

2G（u+iv）=f（ζ）（25）

式（25）中，f（ζ）是由矩形孔口位移边界条件确定的位移函数。

在ζ平面上，有ζ=ρσ=ρei，因此ζ-=ρσ-1，由此可得

ω（ζ）ω′（ζ）=

a（1+ρ σ）/（1-ρ σ）+∑nk=1βk（（ρ σ）k-（ρ σ）-k）-2a/（1-ρ σ-1）2-∑nk=1kβk（（ρ σ-1）k-1+（ρ σ-1）-k-1）（26）

式中：ω（ζ）/ω′（ζ）较为复杂，无法直接写成级数形式，但可以转换为傅里叶级数来计算。

ζ=1∶ω（σ）ω′（σ）=∑SymboleB@k=-SymboleB@δkσk （27）

ζ=α∶ω（ασ）ω′（ασ）=∑SymboleB@k=-SymboleB@γkσk（28）

將式（22）、式（23）、式（27）代入地表应力边界条件式（24）中，可以得到

∑+SymboleB@k=-SymboleB@ak+∑+SymboleB@v=-SymboleB@（vδk+v-1av）+b-k·σk=0（29）

由此可得ak和bk的关系式

b-k=-ak-∑+SymboleB@v=-SymboleB@（v δk+v-1av）（30）

令f（ζ）=∑+SymboleB@k=-SymboleB@Akσk，再把式（22）、式（23）、式（28）、式（30）代入孔口位移边界条件式（25）中，解得

Ak=ak（καk+α-k）+∑+SymboleB@v=-SymboleB@v（δk+v-1α-k-γk+v-1αv-1）·av（31）

Verruij[11]根据半无限空间圆形孔洞问题的对称性，假定所有系数均为纯虚数，矩形边界的问题也具有同样的对称性，因此，可以做出同样的假定，有ak=-ak、bk=-bk。

通过式（30）、式（31）可以求解所有的ak、bk，由此获得问题的弹性解答。

根据式（3）、式（20）、式（21），得到位移表达式

u=Re12Gκφ（ζ）-ω（ζ）ω′（ζ）·φ′（ζ）-ψ（ζ）

v=Im12Gκφ（ζ）-ω（ζ）ω′（ζ）·φ′（ζ）-ψ（ζ）（32）

同样，可以得到应力表达式

σx=Re2φ′（ζ）ω′（ζ）-ω（ζ）·φ″（ζ）ω′（ζ）2-

φ′（ζ）·ω″（ζ）ω′（ζ）3-ψ′（ζ）ω′（ζ）

σy=Re2φ′（ζ）ω′（ζ）+ω（ζ）·φ″（ζ）ω′（ζ）2-

φ′（ζ）·ω″（ζ）ω′（ζ）3+ψ′（ζ）ω′（ζ）

τxy=Imω（ζ）·φ″（ζ）ω′（ζ）2-φ′（ζ）·ω″（ζ）ω′（ζ）3+ψ′（ζ）ω′（ζ）（33）

利用式（32）、式（33）可以求得土体的应力场和位移场。其中，位移场包括了土体整体的刚体位移，假设地表无穷远点的位移为零，则各点的位移再减去地表无穷远点的位移后才得到最终的位移场结果。

3 结果验证

为了验证理论模型的准确性，利用ABAQUS建立了矩形孔收缩的有限元模型，将其得到的有限元解答与理论解进行对比验证。模型宽取100 m，高取50 m，隧道中心点深度d取3 m，矩形隧道宽b取1.5 m，高h取1 m，泊松比μ取0.3，u0取0.1，弹性模量E取10 MPa。限制模型左右两侧的水平位移，Verruij指出有限元模型需要将底部的位移释放，以消除刚体位移对模型的影响，因此，底部不设置边界条件。矩形孔洞周边区域采用自由网格划分技术，远离孔口的区域采用结构化网格划分技术，网格划分后总计9 057个单元。有限元模型如图4所示。

图5为有限元解和本文解产生的地表沉降曲线、地表水平位移曲线、矩形孔周第一主应力分布以及位移场和最大切应力τmax云图。图5（a）中，有限元解和理论解都产生明显的沉降曲线，地表最大沉降值分别为70.1、68.7 mm，误差仅为2.85%，且两种方法计算的地表沉降曲线从中轴线一直到地表远端都保持着较好的一致性。图5（b）中，地表的土体都产生了朝中轴线方向的位移，有限元解和本文解的最大位移分别为距中轴线3.7、3.8 m处，误差为2.7%;有限元解和本文解的最大位移值分别为48、54 mm，误差为11.1%，相比沉降曲线的误差来说要偏大;这是由于有限元的边界条件限制了边界处的水平位移，但没有限制竖向位移，因此，有限元解得到的水平位移结果的误差将会大于竖向位移结果，本文解得到的水平位移解将比有限元解的精度更高。图5（c）中，横坐标为沿着矩形孔周，以矩形底边中点为起点，按逆时针方向所走过的路径长度占总周长的比值;纵坐标为第一主应力值;本文解和有限元解在底边中点、侧边中点和顶边中点处的第一主应力误差分别为7.4%、3%和1.5%，具有较好精度;在矩形角点处，有较为明显的应力集中现象，本文解的应力集中现象比有限元解要更加明显，这是由于有限元解在处理应力集中问题时，应力集中的结果受网格影响较大，同时，本文解的映射函数精度也会对应力集中的结果产生影响，但总体而言，这部分误差可以接受。

图5（d）、（e）为有限元解和本文解的位移场云图和最大切应力云图，左半部分为有限元解的结果，右半部分为本文解的结果。很明显可以看出，有限元解和本文解得出的位移场和应力场基本吻合，从而进一步验证了本文解的可靠性。

如果假设土体为弹塑性土体，采用Tersca屈服准则来判断土体的屈服面，那么最大切应力的等值线可以作为土体弹性区和塑性区的分界线。观察到塑性区最早会出现在矩形孔洞的4个角点，随着矩形孔洞的位移逐渐增加，塑性区将会从4个角点向外拓展相连，形成类矩形状的塑形区。

经过上述对比，可以看出本文的理论解与有限元解之间具有很高的一致性，从而验证了理论方法的正确性和可靠性，可以采用该理论方法继续深入矩形顶管的研究工作。

4 参数分析

4.1 位移场分析

4.1.1 高宽比对矩形隧道沉降槽的影响

算例保证矩形隧道的截面面积为π m2，取高宽比h/b分别为1/2、2/3、1、3/2、2，隧道深度d取3 m，土体泊松比μ取0.3，弹性模量E取10 MPa，位移参数u0取0.1。同时，根据Verruij给出的方法，计算一组圆形隧道来对比，圆形隧道的中心点深度和矩形隧道相同，圆形半径为1 m，圆形孔周收敛模式为均匀径向收敛，采用与矩形隧道相同的土体损失率来计算。

图6给出了圆形隧道以及不同高宽比的矩形隧道所形成的沉降槽，结果表明：若保证相同隧道的横截面积和土体损失率，则圆形隧道产生的沉降槽与高宽比h/b为1时的矩形隧道（正方形）产生的沉降槽基本接近。对于矩形隧道，随着高宽比的增加，沉降槽的宽度逐渐减小，而沉降槽的深度逐渐增加。随着高宽比的改变，沉降槽的形状将会发生变化，当高宽比减小到某一个值以下时，沉降槽将不再是一个类高斯曲线的形状，最大沉降点不再位于隧道中轴线处，轴线处的沉降值将小于两侧。这种现象随着高宽比的减少，将会越来越明显。这种沉降槽形状受高宽比影响的现象，可能是因为沉降槽受到矩形隧道侧边和上下边位移的共同作用叠加，其中顶边向下的位移对沉降槽的影响较大，而侧边的横向位移也会对位移场产生影响。当高宽比较小时，侧边较短，侧边位移对位移场影响较小，影响区域集中在兩条侧边处，所以，此时沉降槽的最大沉降点位于沉降槽中心两侧（对应两条侧边所在位置）;而当高宽比较大时，侧边较长，侧边位移对位移场影响较大，两条侧边的影响区域相互叠加，使得侧边横向位移产生的地表沉降最大值依旧在沉降槽中心处。

4.1.2 不同泊松比对沉降槽和横向位移的影响

选用泊松比μ分别为0.1、0.25、0.35、0.5，弹性模量E取10 MPa，矩形隧道宽b取1.5 m，高h取1 m，隧道中心点深度d取3 m，边界条件位移参数u0取0.1。

图7给出了不同泊松比条件下形成的沉降槽，结果表明：泊松比对沉降槽深度的影响较大，随着泊松比的增大，沉降槽的深度将会明显加大，同时，沉降槽的宽度也会增加。沉降槽的形状也会跟随泊松比的变化而变化，当泊松比很小时（如图7中μ=0.1时的曲线），沉降槽的形状很接近类高斯曲线;随着泊松比的增大，沉降槽的形状将会逐渐发生改变（如图7中μ=0.35时的曲线，明显沉降槽底部趋于平缓）;而泊松比的值超过某一界限，沉降槽的最低点将会向隧道轴线两侧偏移，沉降槽轴线处的土体相较于其两侧的土体发生轻微隆起（如图7中μ=0.5曲线所示），这种现象伴随泊松比的增大愈发明显。这种现象在高宽比分析中，推测是由于侧边横向位移造成的，而泊松比是横向正应变与轴向正应变绝对值的比值，泊松比的取值大小将会很大程度地反映侧边横向位移对地表竖向位移（沉降槽）的影响。泊松比越大，侧边横向位移对沉降槽的影响越大，因此，沉降槽形状发生改变的现象会越明显。

图8为不同泊松比条件下，x=b位置处土体的横向位移曲线图，结果表明：从整体上看，泊松比越大，得到的横向位移越小。横向位移的最大值点集中在隧道中心点深度d附近，泊松比对于横向位移的最大值影响很小。泊松比对横向位移的影响主要体现在隧道以上的土体，泊松比不仅会极大地影响隧道以上土体的横向位移值，而且会改变横向位移曲线的形状。当泊松比偏大时，隧道以上土体的最小横向位移点出现在地表;当泊松比减小到某个值以下，最小横向位移点将会向下偏移。

4.1.3 不同埋深对沉降槽的影响

选用隧道埋深D分别为1、3、5、7 m，隧道埋深D为隧道顶部到地表的距离;泊松比μ取0.3，弹性模量E取10 MPa，矩形隧道宽b取1.5 m，高h取1 m，边界条件位移参数u0取0.1。

图9为不同埋深的矩形隧道产生的地表沉降曲线，结果表明：隧道埋深对沉降槽最大沉降值影响很大，埋深越浅，最大沉降值越大;同时埋深也会影响沉降槽的宽度，沉降槽的宽度随着埋深的增加逐渐减小。埋深较浅时，沉降槽的形状会发生变化，变化情况和低泊松比和较大的高宽比时相同，中轴线上的地表沉降将会略大于其两侧的沉降值。

4.2 应力场分析

4.2.1 不同泊松比对孔周应力的影响

泊松比μ分别取0.1、0.3、0.5，对高宽比h/分别b为2/3、1、1.5进行分析。u0取0.1，弹性模量E取10 MPa，埋深D取3 m，隧道截面面积取π m2。

图10为高宽比分别为2/3、1、1.5的隧道在泊松比0.1、0.3、0.5下的孔周σ1分布。横坐标为孔周点的相对位置，0代表底边中点，0.5为顶边中点。分析结果可知：不同的泊松比将会影响孔周σ1的大小。泊松比越大，孔周σ1整体数值越小，但应力集中的现象将会更加明显，反之亦然。取不同泊松比时，矩形孔洞各个边中点处的第一主应力值相互的比值不会发生改变，这说明泊松比不会影响孔周应力的分布情况，只会影响孔周应力的大小。

4.2.2 不同埋深对孔周应力的影响

矩形隧道埋深D取1、3、5 m，对高宽比h/b为2/3、1、1.5进行分析。位移参数u0取0.1，弹性模量E取10 MPa，泊松比μ取0.3，隧道截面面积取π m2。

图11为高宽比h/b分别为2/3、1、1.5的隧道在埋深1、3、5 m下的孔周σ1分布。结果表明：无论对于哪种工况，埋深越小，矩形隧道孔周在各个位置的第一主应力值都会越小，反之亦然。很明显，埋深1 m相对于埋深3 m时的应力变化非常大，而埋深3 m相对于埋深5 m时的应力变化相对小很多。这是由于埋深较浅时，隧道孔周应力受地表影响较大，对于埋深相对敏感;而埋深较深时，隧道的孔周应力受地表影响较小，随着埋深的加大，孔周应力的变化将会趋于收敛。矩形隧道的埋深和高宽比都会影响孔周应力分布情况。埋深越小，隧道顶边第一主应力的最大值相对于底边的比值将会越大，随着埋深的增大，这一比值将会趋向于1，这与隧道从浅埋到深埋的变化规律相一致。同时，埋深越小，侧边第一主应力的最大值相对于底边的比值将会越大，随着埋深增大，这一比值也会趋向于一个固定值。这说明，当隧道属于超浅埋时，最大第一主应力往往会出现在顶部。

5 结论

基于复变函数方法得出了可以适用于平面应变条件下半无限空间中矩形隧洞的位移边界条件解析解，且通過有限元的计算结果验证了解析解的可靠性;然后基于该方法对不同埋深、不同高宽比、不同泊松比对位移场和应力场的影响进行了敏感性分析，得到以下主要结论：

1）提出的方法是一种解决半无限空间矩形孔收缩问题的复变函数解法，该方法经过与有限元计算结果的对比验证，保证了其可靠性。解析解与有限元解相比，应力场和位移场都比较一致，仅在应力集中的结果上偏于保守，在实际工程设计运用中，采用解析解的结果会更安全;分析了半无限空间矩形隧洞问题的塑性区发展规律，发现塑性区最早会在矩形孔洞的4个角点产生，随着位移的增大，逐渐拓展连通，最后，塑性分界面将在矩形孔洞外围呈类矩形分布。

2）在本文的边界位移条件下，对沉降槽的参数分析结果表明：沉降槽的深度与土体泊松比和矩形隧洞埋深呈负相关，与矩形隧洞的高宽比呈正相关;而沉降槽的宽度与土体泊松比和矩形隧洞的高宽比呈负相关，与矩形隧洞埋深呈正相关;在泊松比较大、埋深较浅、高宽比较大时，矩形隧洞的沉降槽曲线的形状将不再是类高斯曲线，主要表现为沉降槽的中心处相较于其两侧将会产生轻微隆起现象，沉降槽的最大沉降点不再位于沉降槽中心。

3）土体泊松比不会影响孔周应力的分布，只会影响孔周应力的大小。土体泊松比越大，孔周第一主应力的值越大，反之，泊松比越小，孔周第一主应力的值越小。隧道埋深不仅会影响孔周应力的大小，还会影响孔周应力的分布。隧道属于浅埋时，埋深对于孔周应力的影响较大，当埋深较大时，埋深对孔周应力的影响较小。隧道埋深越小，孔周第一主应力的值越小，反之亦然。随着埋深的加大，顶边最大第一主应力的值相较于底边将会逐渐减小，最终两者比值趋向于1，顶边最大第一主应力与侧边的比值也会逐渐减小，最终趋向于一个常数。参考文献：

[1] NAKAMURA H， KUBOTA T， FURUKAWA M， et al. Unified construction of running track tunnel and crossover tunnel for subway by rectangular shape double track cross-section shield machine [J]. Tunnelling and Underground Space Technology， 2003， 18（2/3）： 253-262.

[2] 彭立敏， 王哲， 叶艺超， 等. 矩形顶管技术发展与研究现状[J]. 隧道建设， 2015， 35（1）： 1-8.

PENG L M， WANG Z， YE Y C， et al. Technological development and research status of rectangular pipe jacking method [J]. Tunnel Construction， 2015， 35（1）： 1-8. （in Chinese）

[3]  JEFFERY G B. Plane stress and plane strain in bipolar co-ordinates [J]. Philosophical Transactions of the Royal Society of London Series A， 1921： 265-293.

[4]  MINDLIN R D. Stress distribution around a hole near the edge of a plate under tension[C]//Proceedings of the Society of Experimental Stress Analysis 5， 1948： 56-67.

[5] 施有志， 高轩能. 半无限空间隧道应力与位移的解析延拓法求解[J]. 厦门大学学报（自然科学版）， 2013， 52（5）： 722-727.

SHI Y Z， GAO X N. A complex variable solution for a lined tunnel with arbitrarily shaped section [J]. Journal of Xiamen University （Natural Science）， 2013， 52（5）： 722-727. （in Chinese）

[6]  SAGASETA C. Analysis of undrained soil deformation due to ground loss [J]. Geotechnique， 1987， 37（3）： 301-320.

[7]  VERRUIJT A， BOOKER J R. Surface settlements due to deformation of a tunnel in an elastic half plane [J]. Géotechnique， 1996， 46（4）： 753-756.

[8] 曾彬， 黄达， 刘杰， 等. 双圆盾构隧道施工偏转角对地表变形影响研究[J]. 岩石力学与工程学报， 2015， 34（12）： 2509-2518.

ZENG B， HUANG D， LIU J， et al. Influence of tilting angle on surface deformation during double-o-tube shield tunnel construction [J]. Chinese Journal of Rock Mechanics and Engineering， 2015， 34（12）： 2509-2518. （in Chinese）

[9]  ZENG B， HUANG D. Soil deformation induced by Double-O-Tube shield tunneling with rolling based on stochastic medium theory [J]. Tunnelling and Underground Space Technology， 2016， 60： 165-177.

[10] 魏綱， 张鑫海， 徐银锋. 考虑多因素的类矩形盾构施工引起土体竖向位移研究[J]. 岩石力学与工程学报， 2018， 37（1）： 199-208.

WEIG， ZHANG X H， XU Y F. Deriving vertical displacement of ground due to quasi-rectangular shield tunneling considering multiple factors [J]. Chinese Journal of Rock Mechanics and Engineering， 2018， 37（1）： 199-208. （in Chinese）

[11]  VERRUIJT A. A complex variable solution for a deforming circular tunnel in an elastic half-plane [J]. International Journal for Numerical and Analytical Methods in Geomechanics， 1997， 21（2）： 77-89.

[12]  VERRUIJT A. Deformations of an elastic half plane with a circular cavity [J]. International Journal of Solids and Structures， 1998， 35（21）： 2795-2804.

[13] 吕学金， 刘学斌. 半无限空间重力隧洞洞室临塑状态分析[J]. 应用力学学报， 2006， 23（3）： 500-503，518.

LV X J， LIU X B. Approaching plasticity in half space gravity tunnel [J]. Chinese Journal of Applied Mechanics， 2006， 23（3）： 500-503， 518. （in Chinese）

[14] 陆文超. 地面荷载下浅埋隧道围岩应力的复变函数解法[J]. 江南大学学报， 2002（4）： 409-413.

LU W C. The complex variable method on analysising the stress of surrounding rocks for shallow tunnel under ground load [J]. Journal of Southern Yangtze University， 2002（4）： 409-413. （in Chinese）

[15] 蔚立元， 韩立军， 陈晓鹏， 等. 含圆孔半平面体的弹性分析及其工程应用[J]. 工程力学， 2013， 30（7）： 167-172，199.

YU L Y， HAN L J， CHEN X P， et al. Elastic analysis of a half plane with a circular cavity and its engineering application [J]. Engineering Mechanics， 2013， 30（7）： 167-172，199. （in Chinese）

[16] 张永兴， 王桂林. 高层建筑岩石洞室地基稳定性分析方法与应用[J]. 工程力学， 2007， 24（Sup2）： 110-120.

ZHANG Y X， WANG G L. Stability analysis method and its application in high-rise building with rock cavern foundation [J]. Engineering Mechanics， 2007， 24（Sup2）： 110-120. （in Chinese）

[17] 王立忠， 吕学金. 复变函数分析盾构隧道施工引起的地基变形[J]. 岩土工程学报， 2007， 29（3）： 319-327.

WANG L Z， L X J. A complex variable solution for different kinds of oval deformation around circular tunnel in an elastic half plane [J]. Chinese Journal of Geotechnical Engineering， 2007， 29（3）： 319-327. （in Chinese）

[18]  PARK K H. Elastic solution for tunneling-induced ground movements in clays [J]. International Journal of Geomechanics， 2004， 4（4）： 310-318.

[19] 曾癸森. 弹性半无限域中非圆形隧洞的复变函数解[D]. 北京： 华北电力大学（北京）， 2018.

ZENG G S. Complex variable solution for A non-circular tunnel in an elastic half-plane [D]. Beijing： North China Electric Power University， 2018. （in Chinese）

[20] MUSKHELISHVILI N I. Mathematical theroy ofelasticity [M]. Leyden： International Publishing， 1954.

[21]  ZHOU H， LIU H L， YUAN J R. A novel analytical approach for predicting the noncylindrical pile penetration-induced soil displacement in undrained soil by combining use of cavity expansion and strain path methods [J]. International Journal for Numerical and Analytical Methods in Geomechanics， 2018， 42（11）： 1270-1305.

（編辑 王秀玲）