，，

(华北电力大学 水电与岩土工程研究所，北京 102206 )

1 研究背景

隧洞被广泛应用于水电、交通、采矿等工程。由于隧洞的开挖，围岩原始的平衡状态将会被打破，致使岩体内部的应力和位移重新分布，尤其是隧洞的开挖边界会产生较大的应力集中，如果原始地应力较大和(或)围岩的强度较低时，洞周的岩体将发生破坏，造成地下结构的失稳破坏。为了加强围岩的稳定性，保证隧洞的安全使用，经常在隧洞内设置钢筋混凝土衬砌支护，这样可以改善围岩的应力状态，从而可以限制围岩的失稳破坏。

地下工程的分析方法主要有有限元法、解析分析法等。解析分析法中，Muskhelishvili[1]提出的平面弹性复变函数方法是进行地下工程围岩应力分析的常用解析方法；在弹性力学中，我们已经求出了无限平板中圆形孔洞的基尔西解。通过复变函数方法可以求出许多几何形状复杂的无支护深埋隧洞问题的封闭解，如椭圆形隧洞、矩形隧洞、半圆形隧洞、直墙半圆拱形隧洞、缺口的圆形开口隧洞等[2-6]。

当隧洞有衬砌支护时，利用复变函数方法求解应力位移比无支护隧洞困难得多。不考虑隧洞内的静水压力，从最早的单个圆形隧洞[7]到非圆形隧洞[8-9]，都已经得到了围岩和衬砌内的应力位移解析解。对于带支护的弹塑性岩体也获得了很多成果[10]。但对于隧洞内具有一定水头的有压水工隧洞，目前只获得了圆形衬砌隧洞的应力和位移解析解，陈振民[11]最早开展了这方面的研究，并获得了一些初步成果。本文作者对圆形压力隧洞考虑支护滞后的情况开展了较深入的研究[12]。

而对于带有衬砌的非圆形水工隧洞，当隧洞内有一定水头的静水压力作用时，目前尚未见到利用解析方法来求解它的应力和位移，这正是本文将要解决的问题。

以上的研究大都基于平面应变假定，本文也将问题简化为平面应变问题。当隧洞没有支护时，平面应变问题容易获得应力场和位移场，而当在隧洞中设置衬砌支护时，必须要考虑隧洞开挖后支护的滞后过程。理论上，平面问题只能解决洞室全部开挖完成后立即进行衬砌支护的问题，即认为围岩还没有发生变形。实际上，在设置衬砌之前，随着工作面的推进，围岩已经产生部分变形，设置衬砌后，随着洞室的继续开挖，围岩将继续变形而对衬砌产生作用力。为考虑这一支护滞后问题，可以根据工作面和支护断面的距离来确定位移释放系数η[13-14]。在已知支护前围岩所完成位移的前提下，再利用平面应变的假定，来计算围岩与衬砌的应力、位移解析解[8,12]。

进行应力、位移求解时, 除了平面应变假定外，本文还有如下假定：

(1)围岩与衬砌相互作用时，假定围岩与衬砌之间为完全接触，即认为接触面十分粗糙，岩体与衬砌之间不会产生任何滑动。在接触面上，围岩和衬砌之间的接触应力(法向正应力和切向剪应力)和位移均是连续的[8,12]。

(2)围岩与衬砌在原始地应力和水压力的共同作用下，始终处于各向同性的线弹性状态。

(3)隧洞埋置深度足够大，可以简化为无限域问题。

2 求解的基本原理和方程

2.1 无支护时开挖引起的全部位移

本文利用平面弹性复变函数方法中的保角变换法求解该问题。保角变换方法是通过映射函数把物理平面上的一个给定的区域映射到另一个象平面上具有简单边界形状的区域，如单位圆内或单位圆外。当映射到单位圆外时，其映射函数的一般形式为

(1)

式中：R为正实数；ck为一般复数；ζ=ρeiθ=ρ(cosθ+isinθ)为ζ平面上任一点,ρ为象平面内半径，i为虚数单位。

如图1所示，x轴为隧洞的对称轴，则式(1)中的ck必为实常数。

图1 支护横断面Fig.1 Supported cross-section

对于平面应变问题，位移表达式为[1]：

；

(2)

(3)

将式(1)、式(3)代入式(2)，位移表达式可统一用式(4)表示,即

(4)

式中：κ=3-4μ；G=E/[2(1+μ)]，E为杨氏模量。

(5)

式中：κ1=3-4μ1；G1=E1/[2(1+μ1)]，E1和μ1分别是围岩的杨氏模量和泊松比；φ1(ζ)，ψ1(ζ)是围岩在无支护时对应的单位圆外的解析函数，这2个解析函数可用式(6)、式(7)表示。

因为φ1(ζ)是圆外解析函数，故可以写成

(6)

由文献[2]可求出ψ1(ζ)，即

(7)

式中：P为垂直方向地应力；λ为侧压力系数。由于隧洞的几何形状和外荷载 (地应力) 都关于x轴对称，所以ak，Sk，S′0均为实数，可以由文献[2]得到。

2.2 围岩及衬砌位移的解析函数表示

针对本文的平面弹性接触问题，封闭式支护可视为镶嵌于弹性无限体内的环状弹性体。利用复变函数法中的保角变换方法 ，通过映射函数z=ω(ζ)，将物理平面中形状复杂的支护横断面(图1)，变换到象平面上的圆环区域。圆环的外边界半径为1，内边界半径为R0，圆环内外边界分别为γ1，γ2。其中映射函数由文献[15]确定。

对支护滞后问题，位移释放系数用η(0≤η≤1)表示，如果隧洞在完成了η倍的全部位移后，再进行支护，那么支护前围岩产生的位移由式(8)给出，即

(8)

围岩与衬砌相互作用时，衬砌限制了围岩的部分位移，这部分支护限制围岩产生的位移根据式(4)可写为

(9)

；

(10)

(11)

式中b0，d0，bk，dk为待求常数。

支护后衬砌内任一点的位移根据式(4)，可表示为

(12)

式中：uL，vL分别为衬砌内x，y方向的位移分量；κ2=3-4μ2；G2=E2/[2(1+μ2)]，E2和μ2分别是衬砌材料的杨氏模量和泊松比；φ3(ζ)，ψ3(ζ)是围岩对衬砌作用后，衬砌对应圆环内的解析函数。这2个解析函数可用幂级数式(13)、式(14)表示，即：

；

(13)

(14)

式中p0，q0，ek，fk，gk，hk是待求常数。

因为支护结构也关于x轴对称，所以式中待求的都是实常数。

2.3 求解解析函数的边值条件

当孔边作用的静水压力p较大时，可以忽略隧洞内水的自重对应力的影响，所以可以认为有σx=σy=p，τxy=0(本文规定应力以压为正)，则由应力边界条件可求得

(15)

可以去掉式(15)中的常数，而不影响计算结果。

对于水工隧洞，衬砌的内边界L1上受到静水压力的作用，所以L1上的应力边界条件可以用式(16)表示，即

(16)

式中:σ1=R0σ；R0为象平面上圆环内边界的半径；Xn,Yn为圆环内边界上x,y方向的面力分量。

在衬砌与围岩接触面L2上，因为假定围岩与衬砌为完全接触，则接触面上的法向应力和剪应力连续。L2上的应力连续条件可以用式(17)表示，即

(17)

式中：σ=eiθ；θ为象平面上的极角。

围岩与衬砌之间为完全接触，根据接触面L2上的位移连续条件，则可得到式(18)，即

(18)

由式(5)、式(8)、式(9)、式(12)、式(18)可得到L2上位移连续条件表达式。

2.4 求解φ2(ζ)，ψ2(ζ)，φ3(ζ)，ψ3(ζ)的过程

由于衬砌对围岩的影响随着ζ的增大而减小，则当ζ→时，物理平面中无穷远处的则由可得

b0=d0/κ1。

(19)

由于求解φ2(ζ)，ψ2(ζ)，φ3(ζ)，ψ3(ζ)的方程组个数为无穷多，所以在求解过程中必须取有限项。这里取ζ的最高次幂都为N,即bk，dk，ek，fk，gk，hk中的k=1,2,…,N。

(20)

(21)

3 围岩与衬砌内的应力与位移求解

解析函数φ2(ζ)，ψ2(ζ)，φ3(ζ)，ψ3(ζ)求出后，可以由此求出围岩与衬砌中的应力、位移。

3.1 衬砌内应力

衬砌内任一点的应力分量可由联立式(22)、式(23)求得，即：

σρ+σθ=4Re[Φ(ζ)] ；

(22)

(23)

其中：

式中σρ,σθ,τρθ是正交曲线坐标系下物理平面内的3个应力分量。

3.2 围岩内应力

围岩内任一点的应力分量可由联立式(24)、式(25)求得，即：

σρ+σθ=4Re[φ′(ζ)/ω′(ζ)]；

(24)

(25)

式中φ(ζ)，ψ(ζ)由式(26)、式(27)表示，即：

φ(ζ)=Γω(ζ)+φ1(ζ)+φ2(ζ) ；

(26)

ψ(ζ)=Γ′ω(ζ)+ψ1(ζ)+ψ2(ζ) 。

(27)

式中：Γ=P(1+λ)/4；Γ′=P(λ-1)/2；第1项Γω(ζ),Γ′ω(ζ)表示没开挖隧洞之前围岩所对应的复势函数；第2项φ1(ζ),ψ1(ζ)解析函数指的是支护前开挖隧洞对围岩造成的影响；第3项φ2(ζ),ψ2(ζ)解析函数指的是衬砌对围岩的影响。

3.3 支护后衬砌和围岩的位移

衬砌的位移uL，vL可以根据式(12)求出。

支护后，围岩产生的位移uR，vR可以根据式(5)、式(8)、式(9)求出。

4 算例及参数影响分析

由本文的解析解法可以求出具体算例的应力和位移解，本章只对算例进行应力分析。

以马蹄形隧洞为例，取衬砌厚度为0.3 m，由文献[15]可得到映射函数为

z=ω(ζ)=1.929 7(ζ-0.059 80-0.045 67ζ-1+

0.043 133ζ-2-0.001 77ζ-3-0.010 54ζ-4+
0.002 594ζ-5+0.001 579ζ-6-0.001 67ζ-7+
0.000 711ζ-8+0.000 84ζ-9) 。

(28)

从式(28)中可以看出：R=1.929 7；c0=-0.059 80；c1,…,c9为ζ的负幂次项前系数，可求出R0=0.869 57。用此映射函数可表示隧洞的形状、尺寸。

本算例中岩体的其他计算参数为：围岩弹性模量E1=20 GPa；衬砌弹性模量E2=30 GPa；围岩和衬砌泊松比μ1=μ2=0.2；垂直地应力分量P=5 MPa；水平地应力分量为λP。

算例中地应力、静水压力与孔形都关于x轴对称，所以只取x轴左半部分来讨论，也就是象平面中θ∈[0°,180°]的范围，将此区间等分为180份。

4.1 解析解与数值解的对比分析

本文利用ANSYS软件建模，求出该算例的数值解，并与上文解析解法的结果进行比对。对该算例，取侧压力系数λ为0.5，位移释放系数为0.2，静水压力p为0.5 MPa。

由映射函数绘制的衬砌断面形状如图2(a)所示，取有限元网格划分区域为40 m×40 m，隧洞的最大净空尺寸D为5 m[16]，考虑到支护的滞后效应，ANSYS数值结果应该是模型1(图2(a))和模型2(图2(b))的应力叠加。

图2 ANSYS示意图Fig.2 ANSYS models

对于模型1，在围岩与衬砌接触面上定义接触单元，设置为粗糙接触行为。如图2(a)所示，施加边界条件为：在接触面L2的节点上施加等效释放荷载，衬砌内边界L1上施加静水压力，在A和B点施加x和y方向的约束。由模型1可求出衬砌内外边界的应力。

对于模型2，裸洞尺寸及边界条件如图2(b)：在模型外边界施加垂直地应力P，水平地应力λP,裸洞内边界施加等效释放荷载，在A和B点施加x和y方向的约束。对于围岩边界的应力，需要将模型1求出的围岩边界应力与在原始应力作用下释放η倍等效荷载时围岩边界应力相叠加(模型2)。因此围岩边界的应力值就等于模型1与模型2的围岩边界应力值的叠加。解析解与数值解比对结果如图3所示。

图3 数值解与解析解的对比Fig.3 Comparison between numerical and analytical solutions

由图3可以看出：

(1)用有限元方法求出的马蹄形隧洞的数值解与用复变函数的保角变换法求出的解析解基本吻合，只在某些地方存在较小的差异。可以说明利用复变函数方法求得的应力位移解析解的正确性。

(2)衬砌内外边界和围岩边界的切向正应力比围岩与衬砌接触面上的剪应力和径向应力大得多。接触面上剪应力和径向应力的最大值为2.5 MPa；衬砌外边界上的切向正应力最大值为7.7 MPa，是2.5 MPa的3倍；围岩边界的切向正应力最大值为12 MPa,是2.5 MPa的4.8倍；衬砌内边界上的切向正应力最大值为15.5 MPa,是2.5 MPa的6倍。

4.2 不同侧压力系数下的应力分析

取侧压力系数λ分别为0.5，0.7，1.0，位移释放系数为0.2，衬砌内边界L1上的静水压力p为0.5 MPa。在不同侧压力系数下，马蹄形隧洞围岩与衬砌中的应力计算结果如图4所示。

图4 不同侧压力系数下各边界的应力分布Fig.4 Stress distribution along boundaries in the presence of different values of side compression

由图4可知：

(1)衬砌内外边界上的切向正应力，围岩边界的切向正应力以及接触面上的径向应力在马蹄形隧洞顶部和底部变化较大，且随着侧压力系数的增大而增大。但在拐角(θ=110°)附近相对变化较小，尤其接触面上的径向应力基本没有变化。

(2)围岩边界和衬砌外边界上的切向正应力随着侧压力系数的增大逐渐变得平缓。

(3)衬砌内边界拐角附近切向正应力为所有应力中的最大值。

(4)接触面上的剪应力在顶部两侧附近的区域，随着侧压力系数的增大，剪应力由负值变为正值，说明剪应力的方向发生了变化。在拐角下侧附近的区域，剪应力达到最大值，但这个最大值随侧压力系数的增大而减小。

4.3 不同位移释放系数下的应力分析

取位移释放系数η分别为0.2，0.5，0.8，侧压力系数为0.5，衬砌内边界L1上的静水压力p为0.5 MPa。在不同位移释放系数下，马蹄形隧洞围岩与衬砌中的应力计算结果如图5所示。

图5 不同位移释放系数下各边界的应力分布Fig.5 Stress distribution along boundaries in the presence of different values of displacement release coefficient

由图5可知：

(1)所有应力都在隧洞拐角(θ=110°)附近产生最大的应力集中。

(2)衬砌内外边界上的切向正应力，接触面上的径向应力随着位移释放系数的增大而减小，且变化幅度较大，尤其在拐角附近产生大幅度降低。

(3)在水压为0.5 MPa下，当位移释放系数为0.8时，衬砌内外边界上的切向正应力在顶部和底部出现较小的拉应力。

(4)围岩边界上的切向正应力基本随位移释放系数的增大而增大，但在顶部和底部少部分随位移释放系数的增大而减小。

(5)围岩与衬砌接触面上剪应力的绝对值随位移释放系数的增大而减小。

(6)接触面上的径向应力和剪应力相对于围岩和衬砌上的切向正应力来说，值要小得多。

由衬砌内切向正应力随位移释放系数的增大而减小的幅度较大、围岩边界的切向正应力增大幅度较小可知，如果支护较晚，会使衬砌内切向正应力较小，同时围岩边界上的切向正应力增幅较小。

4.4 不同静水压力下的应力分析

取衬砌内边界上的静水压力p分别为0.5，1.0，2.0，3.0，4.0 MPa，侧压力系数为0.5，位移释放系数为0.2。在不同静水压力下，马蹄形隧洞围岩与衬砌中的应力分布如图6所示。

图6 不同静水压力下各边界的应力分布Fig.6 Stress distribution along boundaries in the presence of different values of static water pressure

由图6可知：

(1)围岩与衬砌上的切向正应力随水压力的增大而减小，且都在拐角处(θ=110°)产生切向正应力集中。衬砌内切向正应力降幅大于围岩。

(2)当水压力为2.0 MPa时，衬砌内外边界上的切向正应力在顶部和底部开始出现拉应力。

(3)围岩与衬砌接触面上的径向应力随着水压力的增大而增大，当水压力较小时，最大的径向应力值发生在拐角处(θ=110°)。当水压力为3.0 MPa时，曲线较为平滑，径向应力在洞周分布比较均匀。

(4)当水压力由0.5 MPa增大到2.0 MPa时，接触面上剪应力的绝对值随着水压力的增大基本减小。而由2.0 MPa增大到3.0 MPa时，剪应力绝对值大幅度增大，当水压力继续增大到4.0 MPa时，剪应力变化不大。

5 结 论

(1) 本文通过复变函数方法获得了非圆形水工隧洞作用有内水压力以及考虑支护滞后的应力位移解析解。并且通过计算马蹄形隧洞，将解析结果与ANSYS数值结果进行了对比验证，计算结果吻合良好，说明本文的理论推导没有错误。

(2) 衬砌和围岩边界的切向正应力随侧压力系数的变化趋势基本一致。当侧压力系数增大，在顶部和底部附近的切向正应力大幅度增大，而在拐角附近的变化不明显。

(3) 衬砌边界的切向正应力随着位移释放系数的增大而减小，当位移释放系数过大时，边界会产生拉应力。围岩边界的切向正应力在拐角附近随位移释放系数增大而增大，其他位置的切向正应力变化不显著。

(4) 由于水压力在内边界产生的是径向应力，当水压力增大时，接触面上的径向应力也增大，而衬砌边界和围岩边界上的切向正应力减小。又因为径向应力比切向正应力值要小得多，水压力过大，衬砌边界将产生拉应力，所以可以适度增大水压力。

(5) 为了使隧洞更安全可靠，要尽量减小衬砌边界的切向正应力，并且不能出现拉应力。当隧洞运行时其内水压力较小时，可以考虑在建设隧洞时推迟支护，从而降低衬砌切向正应力的大小；当正常工况下水压力较大时，隧洞应提前支护，防止边界出现拉应力。

[1] MUSKHELISHVILI N I. Some Basic Problems of the Mathematical Theory of Elasticity[M]. Groningen: P.Noordhoff Ltd.,1953.

[2] 吕爱钟,张路青.地下隧洞力学分析的复变函数方法[M]. 北京:科学出版社,2007.

[3] 施高萍,祝江鸿,李保海,等. 矩形巷道孔边应力的弹性分析[J].岩土力学,2014,35(9):2587-2593.

[4] 祝江鸿.隧洞围岩应力复变函数分析法中的解析函数求解[J]. 应用数学和力学,2013,34(4):345-354.

[5] EXADAKTYLOS G E,STAVROPOULOU M C. A Closed-form Elastic Solution for Stresses and Displacements Around Tunnels[J]. International Journal of Rock Mechanics & Mining Sciences,2002,39(7):905-916.

[6] EXADAKTYLOS G E,LIOLIOS P A,STAVROPOULOU M C. A Semi-analytical Elastic Stress-Displacement Solution for Notched Circular Openings in Rocks[J]. International Journal of Solids & Structures,2003,40(5): 1165-1187.

[7] WANG M B,LI S C. A Complex Variable Solution for Stress and Displacement Field Around a Lined Circular Tunnel at Great Depth[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 2009,33(7):939-951.

[8] LU A Z,ZHANG N,KUANG L. Analytic Solutions of Stress and Displacement for a Non-circular Tunnel at Great Depth Including Support Delay[J]. International Journal of Rock Mechanics & Mining Sciences, 2014,70(9):69-81.

[9] 吕爱钟,覃 媛,陈虹宇. 马蹄形隧洞考虑支护滞后过程的应力分析[J]. 岩土力学,2014, 35(增1):42-48.

[10] 张黎明, 李 鹏, 孙林娜,等. 考虑地下水渗流影响的衬砌隧洞弹塑性分析[J]. 长江科学院院报, 2008, 25(5):84-87.

[11] 陈振民. 地应力对圆形压力隧洞衬砌应力的影响[J]. 江西水专学报,1984:23-30.

[12] LU A Z,ZHANG L Q,ZHANG N. Analytic Stress Solutions for a Circular Pressure Tunnel at Pressure and Great Depth Including Support Delay[J]. International Journal of Rock Mechanics & Mining Sciences, 2011,48(3):514-519.

[13] 李云鹏,刘怀恒. 典型洞室工作面端部围岩变形特征的有限元分析[J]. 西安矿业学院学报,1988,8(2):1-8.

[14] ZHU W S,LIN S S, ZHU J Q,etal. Some Practical Cases of Back Analysis of Underground Opening Deform Ability Concerning Time and Space Effects[C]∥ Proceedings of the Sixth International Conference on Rock Mechanics. Montreal,August 30-September 10,1987:1325-1356.

[15] 吕爱钟. 非圆形硐室封闭整体式支护映射函数确定的新方法[J]. 岩土工程学报,1995,17(4):38-44.

[16] 吕爱钟,蒋斌松,尤春安. 位移反分析有限元网格划分范围的研究[J]. 土木工程学报,1999,32(1):26-30.