毛洞隧道稳定性和破坏模式是工程界普遍关注的问题，其研究方法可借助岩土极限分析上限理论。目前，极限分析上限法多用于研究地基^[1]、边坡^[2]、基坑^[3]等工程稳定性问题，而对于隧道稳定性，该法的应用相对较少。

利用上限法分析隧道稳定性常需假定破坏模式。如Atkinson、Davis、Osman和Klar等^{[4, 5, 6, 7]}以此研究了隧道环向或开挖面稳定性。杨峰、黄茂松等^{[8, 9, 10, 11]}也构建出多种隧道破坏模式，用于研究隧道环向和纵向开挖面稳定性问题。

当上限法与有限元技术相结合，便可在获得隧道稳定性上限解的同时得到其对应的破坏模式。这方面，Yamamoto等^[12]利用极限分析上、下限有限元研究了地表超载作用下双圆隧道稳定性。Abbo等^[13]利用上、下限法分析了不排水条件下方形隧道的稳定性。

上限有限元采用刚体单元与间断线的模式能自动搜索获得特征清晰的刚性运动块体破坏模式。在这方面，Milani等^[14]引入具有旋转自由度的刚性曲边三角形单元，建立序列线性规划模型求解上限解。Hambleton等^[15]采用刚性平动单元，引入节点坐标摄动思想并求解一系列二阶锥规划问题获得上限解。杨峰等^[16]直接建立了刚体平动运动单元上限有限元非线性规划模型，提出初值的获取方法，其实现过程便于编程求解。

刚体平动运动单元上限有限元^[16]在研究岩土、特别是隧道破坏模式方面具有明显的优势。鉴于此，为了探讨毛洞隧道稳定性和破坏模式等相关问题，本文展开毛洞隧道稳定性刚体平动运动单元上限有限元分析，获取土体内摩擦角、黏聚力、隧道埋深等因素对隧道稳定性及破坏模式的影响规律和计算图表。

如图 1所示，速度间断线两侧刚体平动运动单元速度分别为(u_z,v_z)和(u_y,v_y)， 其中节点①、②和③、④重合。

图 1 刚体平动运动单元间的速度间断线 Fig. 1 Velocity discontinuity line between adjacent rigid translator moving elements

图选项

如图 1所示，假定土体服从摩尔-库伦屈服准则，则速度间断线法向和切向相对速度Δv和Δu应满足：

点击浏览公式

点击浏览公式

式中下标i为第i条间断线对应的变量，利用式(2)可将式(1)所示速度间断线约束条件变换为：

点击浏览公式

位于边界上的单元应满足速度边界条件，即

点击浏览公式

式中u_i和v_i为边界单元i的x向和y向速度分量；u-_i和v-_i为该单元所在边界上的切向和法向速度分量；θ_i为边界切向与x轴夹角。

运动单元i面积为非负值，即：

点击浏览公式

设y=f(x)为几何边界函数，则该边界上节点i需满足：

点击浏览公式

由上限理论建立刚体平动运动单元上限有限元的非线性规划模型如下： 目标函数(最小化)：

点击浏览公式

约束条件为：

点击浏览公式

式(7)和(8)中，决策变量为单元速度u_i和v_i(i=1,…,n_t)、 节点坐标x_i和y_i(i=1,…,n_n)， 其中n_t为单元总数，n_d为速度间断线总数，n_v为速度边界上的单元总数，n_b为几何边界上的节点总数。 非线性规划模型采用序列二次规划法求解，初始值通过求解节点固定时的刚体单元上限有限元线性规划模型获得。

图 2为毛洞隧道二维稳定性分析模型。其中隧道直径为D，埋深为H；土体重度为γ，内摩擦角为φ，黏聚力为c。

图 2 毛洞隧道稳定性分析模型 Fig. 2 Analysis model of unlined tunnel stability

图选项

本例中，土体自重是决定隧道稳定性的关键因素。为便于分析，将参数无量纲化，隧道稳定性分析即为求解临界系数γD/c及其相应的破坏模式，此时毛洞隧道恰处于塑流发生时的失稳临界状态，而临界系数γD/c为φ和H/D的函数。

为分析土体内摩擦角φ、隧道埋深与直径比H/D这两个因素对临界系数γD/c的影响规律，选取计算参数如表 1所示。

毛洞隧道稳定性上限有限元网格划分见图 3，以H/D=1.5为例说明。利用对称性只取模型右侧一半，隧道下方和水平延伸长度L₁和L₂取值见表 2。模型左边界约束x向速度，即u=0；下边界和左边界x和y两个方向速度均约束，即u=0,v=0；地表和隧道轮廓边界自由。令土体单位容重的自重功率约束为∫_Av dA=－1， 此时上限有限元非线性规划模型的目标函数为土体临界重度γ，获得γ上限解后，临界系数γD/c随之得到。

图 3 毛洞隧道稳定性上限有限元模型与网格离散 Fig. 3 Finite element upper bound model and meshing of unlined tunnel stability

图选项

如图 3所示，为了获得一个良好的初始解，模型采用结构化网格进行离散化处理。同时，借助节点固定条件下的初始解信息，对可能的破坏区域进行了网格加密。模型单元总数230，速度间断线总数327，决策变量总数728。当H/D取其他值时，模型相应的网格形式与图 3类似。

利用自编刚体平动运动单元上限有限元程序，对表 1参数对应的工况进行计算，得到毛洞隧道临界系数γD/c的曲线如图 4所示。

图 4 毛洞隧道失稳临界系数γD/c计算结果 Fig. 4 Computing result of unlined tunnel unstable critical coefficient γD/c

图选项

由图 4可知，随着土体内摩擦角φ的增大，临界系数γD/c相应增大，反映出φ值增大毛洞隧道稳定性随之变好的特点。

从图 4还可看出，随隧道埋深与直径比H/D的增大，临界系数γD/c值稍有减小，即隧道埋深的增加使得稳定性变差。由于本文分析的隧道埋深小于4倍隧道直径D，且土体自重为唯一荷载，与隧道埋深直接相关，因此隧道埋深增大增加了自重荷载，由此隧道稳定性变差。不过当土体内摩擦角φ较大时，如图 4中φ≥30°的情况，不同H/D对应的γD/c差异变小，此时埋深对隧道稳定性的影响变小。

(1)刚体平动运动单元上限有限元破坏模式说明当H/D=1.5,φ=10°时，刚体平动运动单元上限有限元法所得毛洞隧道破坏模式示意如图 5所示。其中图 5(a)为初始解的网格变形图，由节点固定条件下的刚体平动单元上限有限元计算得到。删除破坏区域外围未发挥作用的速度间断线，得到初始解运动块体图如图 5(b)所示。由图可知，节点固定时初始解对应的运动块体破坏模式能反映出大致的破坏范围，不过滑动面不平滑，其对应的临界系数γD/c上限解为2.61。

图 5 刚体平动运动单元上限有限元破坏模式(H/D=1.5,φ=10°) Fig. 5 Failure mode obtained by FE upper bound method with rigid translatory moving elements (H/D=1.5,φ=10°)

图选项

应用该初始解求解节点可动条件下的刚体平动运动单元上限有限元非线性规划模型，得到毛洞隧道最终网格变形图如图 5(c)所示，此时破坏区域内相互错动的刚体块显示出平滑的滑动面。同样删除未发挥作用的速度间断线，得到最终刚性运动块体图如图 5(d)所示。该图较清晰地反映出了隧道破坏模式，破坏区域显示出类似于滑移线的两组平滑曲线，近似于文献^[9]对隧道工作面破坏所假定并优化后得到的网状块体破坏模式。

最终刚性运动块体对应的临界系数γD/c上限解为2.35。依据上限定理可知，上限解越小精度越高，由此说明了最终上限解精度和破坏模式精细化程度均有提高。

下面从上限有限元所得最终刚体运动块体图方面，探讨毛洞隧道破坏模式的参数影响。

(2)土体内摩擦角φ的影响

图 6为H/D=2时，不同的土体内摩擦角φ对应的毛洞隧道刚性运动块体破坏模式。从图 6(a)可知，隧道上方破坏区域内呈相互交叉错动的两组滑动面，不过隧道上方靠近地表处为整体刚性体。对比图 6(a)~(g)，当内摩擦角φ由15°增加到35°时，隧道上方的横向破坏范围明显收缩，靠近地表处的整体下滑范围显著缩小，刚性块体之间的相互错动更加显著。这些特征均反映出了土体内摩擦角对隧道破坏模式的影响。

图 6 刚体平动运动单元上限有限元所得毛洞隧道破坏模式(H/D=2；单位:m) Fig. 6 Failure mode of unlined tunnel obtained by FE upper bound method with rigid translatory moving elements (H/D=2; unit:m)

图选项

(3)毛洞隧道埋深与直径比H/D的影响

图 7(a)~(d)分别为H/D=1,2,3,4且φ=25°时，刚体平动运动单元上限有限元所得毛洞隧道刚性运动块体破坏模式。可以看出，H/D不同时，隧道破坏模式的形态大致相似，而当H/D增大时，隧道横向的破坏范围将明显增加，隧道边墙处的破坏范围向下扩展。

图 7 刚体平动运动单元上限有限元所得毛洞隧道破坏模式(φ=25°;单位：m) Fig. 7 Failure mode of unlined tunnel obtained by FE upper bound method with rigid translatory moving elements (φ=25°;unit:m)

图选项

利用自编刚体平动运动单元上限有限元程序研究了毛洞隧道稳定性和破坏模式，结论如下：

(1) 毛洞隧道失稳临界系数γD/c随土体内摩擦角φ的增大而增大，即随着内摩擦角的增加，隧道稳定性变好；γD/c随隧道埋深与直径比H/D的增加而减小；但H/D对γD/c的影响并不显著，特别是对于内摩擦角较大的情况。

(2) 刚体平动运动单元上限有限元所得毛洞隧道破坏模式大致呈两组交叉的滑动面形态；隧道上方靠近地表范围存在整体下滑体；随着土体内摩擦角φ增大，隧道横向破坏范围明显收缩，地表处刚性体范围显著缩小，下方刚性块体间的相互错动更加显著；隧道埋深与直径比H/D对破坏模式的形态影响不明显，当H/D增大时，隧道横向破坏范围增加，边墙破坏范围向下扩展。