考虑纵坡的隧道掌子面稳定性上限有限元分析

李少庭

(中铁十八局集团第五工程有限公司，天津市 300222)

1 引言

隧道掌子面稳定性是工程中非常关注的问题，文献调研发现，从理论分析或数值计算方面对此已开展不少工作。例如：杨峰、孙雁军、阳军生等分别以构建破坏模式、应用刚体平动运动单元和自适应上限有限元方法，研究了浅埋隧道掌子面稳定性问题；吕玺琳和黄茂松等采用极限分析上限法和极限平衡法，研究了盾构开挖面极限支护压力；王金麒等考虑掌子面旋转破坏模式，依据村山公式对数螺旋线假定及上限原理，分析了隧道掌子面极限支护压力；周森等构造了隧道双向开挖工作面刚体滑块体系破坏模式，获得隧道工作面临界支护反力变化规律。

上述研究将隧道掌子面简化为二维平面问题。按三维空间问题考虑，赵明华等开展了基于筒仓理论极限平衡法的上硬下软地层盾构隧道开挖面极限支护力分析；张明聚等按三维楔体模型和太沙基松动压力的解法进行了砂土地层深埋盾构隧道开挖面稳定性分析；刘文洁等应用“楔形体-棱柱体”极限平衡模型，分析了各向异性地基中盾构隧道开挖面稳定性和极限支护压力；凌同华等应用三维有限差分软件模拟隧道开挖过程，求解了开挖面极限支护压力。

现有的隧道掌子面稳定性分析均假定线路为平坡形式。实际上，工程中大量隧道设置不同程度的纵坡，有些纵坡坡度较大，如上海人民路越江隧道最大纵坡为5.9%，瑞士洛桑地铁2号线部分区段最大坡度达17.7%。在这种条件下，隧道纵向坡度对建设期间掌子面稳定性有不可忽视的影响，有必要深入考量。对此，雷明锋等基于筒仓理论，采用极限平衡分析法建立了迎坡施工条件下盾构隧道开挖面极限支护力的计算模型，认为极限支护力随线路坡度增大而呈线性显著增大，土体破裂角随线路坡度增大反而线性减小；周俊等采用极限分析上限法，建立土体破坏模型，推导上坡条件下的开挖面极限支护压力理论公式，得到与文献[14]类似的结论；赵智涛等采用上限法和数值分析法，研究了爬坡开挖与水平开挖的破坏区影响范围及水平位移变化规律，认为爬坡开挖影响范围较大，影响位置比水平开挖略低。

该文针对隧道掌子面稳定性和失稳破坏模式问题，考虑纵坡的影响，采用高阶单元上限有限元法开展多参数正交影响分析，探讨无量纲临界加载系数和耗散能密度表征的破坏模式及其演化规律，以期研究结果能为实际隧道工程建设提供一定的参考与借鉴。

2 考虑纵坡的隧道掌子面稳定性模型

2.1 问题描述与假定

考虑纵坡的隧道掌子面稳定性分析是三维空间问题。作为初步研究，将其简化为沿隧道纵向中心线剖切的二维平面应变问题(图1)，一方面能反映出掌子面稳定性和破坏演变基本规律；另一方面对于跨度较大的隧道而言较为适用。

图1 考虑纵坡的隧道掌子面稳定性分析模型

如图1所示，具体简化及相应参数如下：① 隧道高度D，埋深H，地表水平自由；② 围岩为均质体，破坏时服从相关联流动法则的摩尔-库仑材料，内摩擦角φ，黏聚力c；③ 平坡时，隧道掌子面竖直自由；下坡时，隧道后方整体绕掌子面中心顺时针旋转固定角度(如θ)；上坡时，变为逆时针旋转；④ 掌子面后方均为紧跟的刚性支护，掌子面自由无支护；⑤ 仅考虑围岩自重作用下掌子面失稳破坏，围岩干燥无水。

模型只承受围岩自重荷载，按重度增加法的思路开展掌子面稳定性分析，即为求解掌子面失稳临界重度γcr，转为无量纲形式：

Ncr=γcrH/c=f(φ,θ,H/D)

(1)

无量纲参数Ncr可称为临界加载系数，对比围岩实际重度对应的加载系数，即可评判掌子面稳定性。Ncr为内摩擦角φ、纵向坡度角θ和埋深比H/D的函数。其解析关系不易获得，该文采用上限有限元计算求解。

2.2 上限有限元分析模型

隧道掌子面稳定性上限有限元分析模型仍可参考图1。为清晰起见，未显示三角形网格。模型水平和竖向延伸范围取L2+L3和H+D+L1，延伸参数取值需消除边界效应的影响。

除掌子面和地表自由外，其余边界均同时约束x和y方向的速度，即u=0，v=0。模型的坐标原点取为掌子面中点，x和y正方向为向右和向上，速度正方向与此一致。

服从于：A11X1+A12X2=0

A21X1=B

(2)

X2≥0

2.3 围岩剪胀角的考虑

当围岩内摩擦角φ较大时，按相关联流动法则要求，将高估破坏发生时的剪胀效应。引入剪胀角φ并采用文献[17]提出的处理方法，设置等效参数φ*和c*，见公式(3)。上限有限元分析时取φ/φ=1和φ/φ=0两种情况。

(3)

2.4 计算参数确定

表1 隧道几何参数及围岩参数

3 临界加载系数Ncr分析

3.1 Ncr对比验证

为印证该文上限有限元结果的可靠性，将平坡条件下(θ=0°)隧道掌子面稳定性临界加载系数Ncr计算结果与文献[2]对比如表2所示。需说明的是，文献[2]应用刚体平动运动单元上限有限元，而该文采用基于六节点高阶三角形单元的塑性变形单元上限有限元方法。由表2可知：该文结果略大于文献[2]对应结果，按上限原理，精度略低。其中最大相对误差为6.29%，对应于H/D=3,φ=30°的情况；最小相对误差0.56%，对应于H/D=2,φ=5°。内摩擦角φ越大，该文结果相对误差越大，而与埋深比H/D关系不大。总体上看，Ncr相对误差在可接受范围内，适合进行后续考虑纵坡的隧道掌子面稳定性分析。

表2 掌子面临界加载系数Ncr结果对比(θ=0°)

3.2 临界加载系数Ncr变化规律

采用基于六节点高阶单元的塑性变形单元上限有限元，以表1各参数的正交组合开展系列计算，得到临界加载系数Ncr计算结果见表3，同时绘制Ncr散点图(图2)。工程实际应用时，可大致依据实际的N值进行判断，若N

表3 隧道掌子面稳定性临界加载系数Ncr上限有限元计算结果

注：括号内数据表示φ/φ=0对应的结果。

图2 隧道掌子面稳定性临界加载系数Ncr上限解

由表3、图2可知：随着内摩擦角φ增加，Ncr逐渐增大；而φ不变时，Ncr随着H/D的增加近似线性增加。单独考虑围岩临界重度γcr，φ较小时γcr随着H/D增加而逐渐减小，说明隧道掌子面上覆土层更倾向于荷载作用；而φ较大时，γcr随着H/D增加到特定值后不再变化，此时埋深影响基本消失。变化规律与文献[2]一致。

围岩剪胀角φ对掌子面稳定性亦有影响，考虑φ/φ=0这种特例情况，计算Ncr得到的相应结果均表示为表3括号内数据。总体上看，当内摩擦角较小时，φ/φ=0和φ/φ=1对应的Ncr数值差异不大。如θ=-10°、φ=5°、H/D=1时，φ/φ=0与φ/φ=1对应结果基本相同，而θ=-10°、φ=30°、H/D=4时，φ/φ=0比φ/φ=1计算结果小24.3%。说明围岩内摩擦角较大时，需要考虑剪胀效应对掌子面稳定性分析结果的影响。

结合图2和表3，分析纵向坡度对掌子面稳定性的影响。以纵向坡度角θ=0°为参照，发现θ=10°对应的Ncr值比相同条件下普遍小4.5%～23.1%，说明隧道上坡掘进掌子面稳定性对应降低；特别是内摩擦角φ越大，相对降低量越大，而埋深比H/D的变化对相对差异影响不明显。进一步对比θ=-10°的情况，发现下坡条件下Ncr值比平坡相同条件下普遍大4.4%～22.3%，说明隧道下坡掘进掌子面稳定性相对较好，而其影响规律与上坡条件基本近似。

由图2中Ncr散点图可看出：其值随围岩内摩擦角φ的变化呈非线性增长趋势。为便于寻找规律，下面进行Ncr与φ和H/D公式拟合。隧道掌子面稳定性可归结于土压力问题，针对不同的纵向坡度角θ，拟合公式为：

Ncr=a1H/D-a2tanφtan2(π/4-φ/2)-a3·tanφ·(H/D)2

(4)

式(4)中包含3个拟合系数a1、a2和a3，均是隧道埋深比H/D的函数，具体见表4。

表4 系数a1、a2和a3求解公式

图2中连续曲线即为拟合公式对应的计算结果。整体上看，该公式拟合效果较好，相对误差处于±7%范围；对于平坡、下坡和上坡情况，Ncr值拟合值最大相对误差分别为5.3%、5.7%和-7.0%。

4 破坏模式讨论

4.1 破坏模式及对比分析

考虑到破坏发生时将出现局部区域的塑性区和剪切带，为精确模拟，除采用高阶三角形单元外，在网格划分时，进行了细致的处理。

如图3所示，以θ=0°、φ=20°、H/D=2计算结果为例进行说明。从图3(a)可看出：在进行上限有限元网格划分时，对隧道掌子面前后、下方以及上方直至地表进行了网格加密；在此基础上，对掌子面前方局部范围又进行了整体加密，以提高计算精度。该例共划分六节点三角形单元总数为8 735个，按此规律，其他条件下网格划分单元总数最高约15 000个。

图3 平坡时掌子面稳定性分析网格划分及破坏模式(θ=0°,φ=20°,H/D=2)

图3(b)为上限有限元计算获得的围岩耗散能密度分布云图。耗散能密度可反映失稳临界状态围岩的塑性应变率。考虑到破坏发生时各区域耗散能密度差异很大，为清晰起见，将其数值归一化后进行了开方处理，再绘制图3(b)。图3(b)中白色区域(取值为0)代表围岩未发生塑性变形，而黑色区域(取值为1)即代表该处围岩塑性应变率最大。可看出，破坏时掌子面上下位置均形成两条主要的剪切带延伸至地表，掌子面本身形成三角形楔形区无耗散能，表明该范围为刚性运动，其前方形成耗散能密度逐渐减小的成片区域，表明其为破坏发生时，相互约束运动形成的塑性区。

将相同条件下，文献[2]应用运动单元上限有限元得到的滑移线网破坏模式示意如图3(c)所示。对比图3(b)、(c)可知：破坏模式表现形式不同，但大体范围和运动趋势一致。该文此例和文献[2]临界加载系数Ncr计算结果分别为11.34和11.09，相对误差2.25%。

4.2 耗散能密度表示的破坏模式演化规律

不同参数取值条件下以耗散能密度表示的隧道掌子面破坏模式示意如图4所示，讨论其演化规律。

从图4(a)可知，当θ=10°,φ=15°,H/D=2时，围岩耗散能密度与图3(b)的基本特征近似。由于隧道此时为上坡形式，掌子面附近的剪切带和塑性区形成逆时针的旋转，此时对应的临界加载系数Ncr=7.81。相同条件下，当θ=-10°，即转为下坡形式时，破坏模式示意如图4(b)所示。同样其特征与前者近似，不过剪切带和塑性区形成顺时针的旋转，其路径和范围有所增大，可以推断此时掌子面稳定性将变好一些。对应的临界加载系数Ncr=9.97，较之图4(a)的情况增加约21.7%。

当埋深比和内摩擦角均有所增加时，如H/D=4,φ=20°，掌子面破坏模式如图4(c)、(d)所示。从图4(c)、(d)可看出：此时破坏模式整体形式仍近似，与埋深比较小时差异体现在：① 掌子面附近剪切带和塑性区更加集中，说明局部范围破坏更容易发生；② 隧道上方直至接近地表的两条剪切带有收缩相交的趋势，推测埋深继续增加，破坏可能不至于达到地表；③ 地表下方破坏区域内部，存在大范围的整体刚性运动区，说明此范围为掌子面破坏的主要外荷载来源区，对应于掌子面附近为破坏抵抗区。

图4(c)、(d)为隧道上坡与下坡掌子面破坏的差异，与图4(a)、(b)类似，即下坡条件下掌子面前方剪切带和塑性区路径和范围增大，围岩抵抗破坏能力增加，相对稳定性变好。再次说明隧道纵坡坡度对施工期间隧道掌子面稳定性有不容忽视的影响。

5 结论

(1)自重作用下，隧道纵向坡度对掌子面稳定性存在不容忽视的影响。上坡掘进时(坡度角θ=10°)，Ncr值普遍小于平坡相同条件下4.5%～23.1%，且围岩内摩擦角φ越大，相对降低量越大。下坡(θ=-10°)条件下对应Ncr值大于平坡4.4%～22.3%。

图4 考虑纵坡时的隧道掌子面失稳破坏模式

(2)随着内摩擦φ增加，Ncr逐渐增大；当φ较小时，临界重度γcr随着H/D增加而逐渐减小；而φ较大情况需考虑围岩剪胀的影响。

(3)隧道自重作用下发生破坏时，掌子面形成两条主要剪切带、三角形刚性运动楔形区和前方耗散能密度逐渐减小的塑性区。上坡条件的掌子面附近剪切带和塑性区发生旋转，剪切带、塑性区路径和范围有所减小，围岩稳定性变差；下坡条件下规律与此相反，稳定性相对平坡时变好。