砂卵石地层盾构隧道掌子面稳定性理论分析

摘要

盾构隧道穿越砂卵石地层掘进施工时，由于砂卵石的地层特性及施工扰动，隧道掌子面容易失稳而引发事故。为研究砂卵石地层盾构隧道掌子面稳定所需的最小支护力，结合砂卵石地层的掌子面失稳模式，在三维Murayama模型的基础上，将掌子面前方土体的滑动面形式假设为对数螺旋滑移线，结合上部松动区域的破坏机制，引入椭球体放矿理论，对传统的棱柱体筒仓模型进行改进，进一步通过极限平衡法推导出砂卵石地层的掌子面极限支护力求解公式；最后通过离散元数值模拟验证所推公式，两者的误差为3.1 %。结果表明所得公式可对砂卵石地层的掌子面极限支护力进行较为精确的求解。研究成果以期为砂卵石地层掌子面极限支护力计算理论提供依据。

关键词

盾构隧道; 砂卵石地层; 掌子面稳定性; 理论分析; 极限平衡法

盾构隧道施工技术具有独特的优势，现阶段我国城市地铁常采用盾构法进行施工。然而在施工过程中，经常会出现因掌子面支护力不足引起的掌子面失稳，进而引发隧道掌子面涌水突泥、地层不均匀沉降甚至地表坍塌等伴生病害，给隧道施工带来了极大的威胁。合理的掌子面支护力（盾构机推力）将有利于维持掌子面的稳定，目前国内外已有大量学者对盾构隧道掌子面稳定性展开了研究，主要研究方法有模型试验、数值模拟以及理论分析。

在模型试验方面，Mair^［1］采用离心模型试验研究了砂土与黏土地层中的掌子面主动失稳机制；牛豪爽等^［2］研究了不同内摩擦角干土与不同水头压力差饱和土条件下的掌子面支护力变化规律和失稳机制；周鹏等^［3］使用泥水盾构模型装置，模拟泥水盾构开挖面的掘进失稳过程，分析地层水土应力变化规律，总结了不同水位下泥水盾构在砂土地层中开挖面主动失稳的破坏形式； Hu等^［4］学者基于成都地铁一号线开展了富含砂卵石地层的缩尺土压平衡（EPB）盾构隧道模型试验，充分考虑了卵石地层的离散性与土压平衡盾构机的机械驱动扰动影响，揭示了卵石土地层失稳运动特性。

在数值模拟方面，Vermeer等^［5］采用有限元数值模拟软件Plaxis3D，开展了盾构隧道掌子面稳定性影响因素分析；Weng等^［6］通过离心模型试验研究了渗流与隧道纵向倾角对软黏土地层盾构隧道掌子面稳定性的影响，并结合有限元数值仿真分析了掌子面的渐进破坏机理；Huang等^［7］通过有限元数值分析软件ABAQUS对隧道掌子面稳定性进行了相关计算；Li等^［8］通过FLAC3D的数值模拟结果提出了一种对数螺旋机制。

在理论分析方面，Mollon和Dias^［9］基于极限分析理论和可靠度分析对浅埋盾构隧道的掌子面稳定进行了研究；Mollon等^［10］基于极限分析理论提出了一种接近地层运动场轨迹的掌子面三维滑动破坏机制； Cheng等^［11］通过研究发现掌子面前方的破坏区分为破坏滑移区与上部松动区，其中破坏滑移区是由直线旋转和对数螺旋线旋转而成的，而上部松动区在破坏演化过程中被定义为不同阶段的锥体或圆柱体。

现阶段掌子面稳定性研究主要针对砂性土和黏性土，而砂卵石地层的研究较少，由于砂卵石地层的松散性、颗粒离散性、几乎没有黏聚力等特点，地层受到机械扰动破坏机理与其他地层有较大的差异性。目前国内外关于砂卵石地层隧道掌子面极限支护力的研究相对较少，但也有部分学者在这一领域取得了一些有价值的成果。最早Monnet和Chaffois^［12］首次提到了针对砂卵石地层盾构掌子面稳定性的研究；Anagnostou和Kovári^［13-14］基于楔型体模型，考虑土压平衡盾构施工中地下水渗流产生的渗透力作用于楔型体内的影响。国外地铁施工中遇到砂卵石地层的报道相对较少，后续关于砂卵石地层掌子面稳定研究以国内学者为主，王明年等^［15］基于颗粒离散元数值模拟方法，发现当掌子面支护力较小的工况下砂卵石地层的变形沿着曲面滑动；葛嘉诚^［16］、胡敏^［17］基于砂卵石土物理力学特性对盾构开挖面稳定性展开研究；范祚文^［18］、胡雄玉^［19］等通过室内模型试验对砂卵石地层盾构隧道掌子面破坏特性进行了研究；李伟平等^［20］基于室内模型试验研究浅埋砂卵石地层盾构隧道开挖力学响应和掌子面变形规律；王俊等^［21-22］通过室内模型试验结合PFC^3D（三维颗粒流数值模拟软件）离散元数值模拟等手段，对土压盾构隧道掌子面失稳诱发砂卵石地层中地层变形特征展开了研究；宋伟涛等^［23］对刀盘旋转扰动下的砂卵石地层掌子面失稳机理进行了探索。

现阶段研究砂卵石地层掌子面稳定性的方法主要有数值分析、模型试验。室内模型试验由于试验材料选取的偏差，边界条件的影响难以排除等缺点，导致试验数据的准确性存在一定问题，而单纯的数值模拟由于本构模型及求解程序等方面的限制往往导致其生成的结果难以让人信服，且建模及计算分析耗时较长。因此本文结合砂卵石地层盾构隧道掌子面失稳机理，开展掌子面稳定极限支护力的理论分析模型研究。本文是基于椭球体理论在砂卵石地层盾构隧道掌子面稳定性理论分析的应用，是对于盾构隧道穿越砂卵石地层掌子面极限支护力求解的理论补充。该方法能充分展现砂卵石地层在受到盾构隧道掘进过程中的扰动而变形及坍塌的特性。且该方法相较于其他研究手段更为便捷，后续工程应用只需根据实际工程调整理论公式中对应的地层参数及盾构掘进参数，即可得到相应的较为准确的极限支护力解析解。

1 椭球体理论

1938年前苏联学者米纳耶夫提出放矿放出的形状为椭球体；1952前苏联学者马拉霍夫的著作《崩落矿块放矿》中系统完整地论述了在松散介质的放出试验中，放出体的最终形态为椭球体；20世纪末期，我国学者对崩落放矿理论进行深入研究，使其在实际矿山生产中得到了较好应用，椭球体理论^［24-25］如图1所示，将松散介质颗粒放置于料斗中，当料斗底部打开一定宽度的放出口时，颗粒在重力作用下从放出口流出，并在一定时间后形成图中所示的流出椭球体，与此同时，料斗内还会形成一个对应的极限椭球体，两个区域之间的土体颗粒会产生相对松动，但并不会流出，极限椭球体范围以外的颗粒将保持静止不动。有关模型试验也证明，隧道开挖导致的掌子面前方松动区域在极限支护力的作用下近似为椭圆状^［26-27］。

图1  椭球体理论

Fig. 1  Ellipsoid theory

2 极限支护力计算

2.1　计算模型

结合椭球体理论与极限平衡法的对数螺旋线模型（Murayama model）建立一种新的掌子面稳定性分析模型，如图2所示，图中，h为上部塌落体高度；D为已知隧道直径；l_a为原点O到开挖面的水平距离；l_j为掌子面到椭球体中心的水平距离；l_s为原点O到开挖面中点的垂直距离；σ_t为维持掌子面稳定的极限支护力；h_s为椭球体最低点到隧道拱顶的高度；对数螺旋线初始半径为r₀；滑移面底端对应的对数螺旋线半径为r_a；a为椭球体的长半轴；b为椭球体的短半轴；r_a与r₀的夹角为θ_a；本模型将能够较好地模拟砂卵石地层盾构隧道开挖而导致的掌子面前方土体运动与失稳情况。其中，对数螺旋线方程表示为

$$r=r_{\mathrm{0}}\mathrm{e}\mathrm{x}\mathrm{p}\left(\theta \mathrm{t}\mathrm{a}\mathrm{n}\phi \right)$$

（1）

图2  掌子面稳定性分析模型

Fig. 2  Stability analysis model of excavation face

令滑移面顶端与水平面垂直，根据关联流动法则，r₀与水平面夹角为土体内摩擦角φ，且滑移面底端与水平面的夹角为π/4+φ/2，根据几何关系，可知θ_a=π/4-φ/2。在图2b中，σ_q为为上部松动荷载。由图2b中的几何关系，可得该模型各几何参数的关系如式（2）—（7）：

$$r_{\mathrm{0}}=\frac{D}{\mathrm{s}\mathrm{i}\mathrm{n}\left(\frac{\mathrm{\pi }}{\mathrm{4}}+\frac{\phi }{\mathrm{2}}\right)\cdot {\mathrm{e}}^{{\theta }_{\mathrm{a}}\mathrm{t}\mathrm{a}\mathrm{n}\phi }-\mathrm{s}\mathrm{i}\mathrm{n}\phi }$$

（2）

$$r_{\mathrm{a}}=r_{\mathrm{0}}{\mathrm{e}}^{{\theta }_{\mathrm{a}}\mathrm{t}\mathrm{a}\mathrm{n}\phi }$$

（3）

$$l_{\mathrm{a}}=r_{\mathrm{a}}\mathrm{c}\mathrm{o}\mathrm{s}({\theta }_{\mathrm{a}}+\phi )$$

（4）

$$l_{\mathrm{c}}=(r_{\mathrm{0}}\mathrm{c}\mathrm{o}\mathrm{s}\phi -l_{\mathrm{a}})/\mathrm{2}$$

（5）

$$h_{\mathrm{s}}=\frac{r_{\mathrm{0}}{l}_{\mathrm{j}}}{l_{\mathrm{a}}}\mathrm{s}\mathrm{i}\mathrm{n}\phi$$

（6）

$$h=\frac{\mathrm{1}}{\mathrm{1}-{\epsilon }^{\mathrm{2}}}\cdot \frac{L^{\mathrm{2}}}{\mathrm{4}{h}_{\mathrm{s}}}$$

（7）

式中：L=2l_j为开挖面前方滑动土体的宽度；ε为偏心率；假定椭球体的轴比等于侧压力系数K₀，则1-ε²=b²/a²。其中，$K_{\mathrm{0}}=\mathrm{1}-\mathrm{s}\mathrm{i}\mathrm{n}\phi$，或$K_{\mathrm{0}}=\frac{\nu }{\mathrm{1}-\nu }$。

由图2几何关系，进一步得椭球体的几何参数为：短轴长$b=a\sqrt[]{\mathrm{1}-{\epsilon }^{\mathrm{2}}}$，长轴长$a=\frac{\left(l_{\mathrm{a}}+\mathrm{1}\right)h_{\mathrm{s}}}{l_{\mathrm{j}}}$，其中$h+h_{\mathrm{s}}=\mathrm{2}a$。

2.2　上覆松动土压力计算

为简化分析，将上部松动截头椭球区域简化为半椭球体+椭圆台^［28］，其中椭圆台部分的受力分析如图3所示。

图3  椭圆台受力分析

Fig. 3  Force analysis of elliptic table

设椭圆台长轴方向侧面与竖直方向的夹角为β₁，椭圆台短轴方向侧面与竖直方向的夹角为β₂，β表示任意侧面与竖直方向夹角。

由图2b中的几何关系，可以求得：$\mathrm{t}\mathrm{a}\mathrm{n}{\beta }_{\mathrm{1}}=\frac{r_{\mathrm{0}}\mathrm{s}\mathrm{i}\mathrm{n}\phi }{a-l_{\mathrm{j}}}$，$\mathrm{t}\mathrm{a}\mathrm{n}{\beta }_{\mathrm{2}}=\frac{r_{\mathrm{0}}\mathrm{s}\mathrm{i}\mathrm{n}\phi }{b-l_{\mathrm{j}}}$。

令$\left\{\begin{array}{l}B\mathrm{\text{'}}\mathrm{\text{'}}=B+\mathrm{2}Z\mathrm{t}\mathrm{a}\mathrm{n}{\beta }_{\mathrm{1}}\\ L\mathrm{\text{'}}\mathrm{\text{'}}=L+\mathrm{2}Z\mathrm{t}\mathrm{a}\mathrm{n}{\beta }_{\mathrm{2}}\end{array}\right.$， 且B₁≥L₁。

式中：$B\mathrm{\text{'}}\mathrm{\text{'}}$与$L\mathrm{\text{'}}\mathrm{\text{'}}$分别为椭圆台顶面椭圆长轴和短轴长度；B与L分别为椭圆台底面椭圆长轴和短轴长度，且$B=\frac{\mathrm{\pi }D}{\mathrm{4}}$，$L=\mathrm{2}{l}_{\mathrm{j}}$； Z为椭圆台高。

在Z方向建立微分体的静力平衡方程：

$$\frac{\mathrm{\pi }}{\mathrm{4}}{\sigma }_{\mathrm{v}}{B}_{\mathrm{1}}{L}_{\mathrm{1}}+\frac{\mathrm{\pi }}{\mathrm{4}}\gamma B_{\mathrm{1}}{L}_{\mathrm{1}}\mathrm{d}z=\frac{\mathrm{\pi }}{\mathrm{4}}\left({\sigma }_{\mathrm{v}}+\mathrm{d}{\sigma }_{\mathrm{v}}\right)B_{\mathrm{1}}{L}_{\mathrm{1}}+\tau \left[\mathrm{\pi }{L}_{\mathrm{1}}+\mathrm{2}\left(B_{\mathrm{1}}-L_{\mathrm{1}}\right)\right]\mathrm{d}z$$

（8）

式中：σ_v椭圆台竖直方向正应力； $\tau$为侧面剪应力。

椭圆松动区域内滑动土体服从Mohr-Coulomb准则，故剪应力 $\tau$为

$$\tau =c+K_{\mathrm{0}}{\sigma }_{\mathrm{v}}\mathrm{c}\mathrm{o}\mathrm{s}\beta \mathrm{t}\mathrm{a}\mathrm{n}\phi$$

（9）

式中：c为内部黏聚力。

将式（9）代入式（8）中，得：

$${\sigma }_{\mathrm{v}}{B}_{\mathrm{1}}{L}_{\mathrm{1}}+\gamma B_{\mathrm{1}}{L}_{\mathrm{1}}\mathrm{d}z=\left({\sigma }_{\mathrm{v}}+\mathrm{d}{\sigma }_{\mathrm{v}}\right)B_{\mathrm{1}}{L}_{\mathrm{1}}+\left(c+K_{\mathrm{0}}{\sigma }_{\mathrm{v}}\mathrm{c}\mathrm{o}\mathrm{s}\beta \mathrm{t}\mathrm{a}\mathrm{n}\phi \right)\left[\mathrm{4}{L}_{\mathrm{1}}+\frac{\mathrm{8}}{\mathrm{\pi }}\left(B_{\mathrm{1}}-L_{\mathrm{1}}\right)\right]\mathrm{d}z$$

（10）

求其通解得：

$${\sigma }_{\mathrm{v}}={\left(B\mathrm{\text{'}}\mathrm{\text{'}}-\mathrm{2}z\mathrm{t}\mathrm{a}\mathrm{n}{\beta }_{\mathrm{1}}\right)}^{n_{\mathrm{1}}}{\left(L\mathrm{\text{'}}\mathrm{\text{'}}+\mathrm{2}z\mathrm{t}\mathrm{a}\mathrm{n}{\beta }_{\mathrm{2}}\right)}^{n_{\mathrm{2}}}[c+{\int }_{\mathrm{0}}^z\frac{\gamma -\frac{\mathrm{4}\left(\mathrm{\pi }-\mathrm{2}\right)c}{\pi \left(B\mathrm{\text{'}}\mathrm{\text{'}}-\mathrm{2}z\mathrm{t}\mathrm{a}\mathrm{n}{\beta }_{\mathrm{1}}\right)}-\frac{\mathrm{8}c}{\pi \left(L\mathrm{\text{'}}\mathrm{\text{'}}-\mathrm{2}z\mathrm{t}\mathrm{a}\mathrm{n}{\beta }_{\mathrm{2}}\right)}}{{\left(B\mathrm{\text{'}}\mathrm{\text{'}}-\mathrm{2}z\mathrm{t}\mathrm{a}\mathrm{n}{\beta }_{\mathrm{1}}\right)}^{n_{\mathrm{1}}}{\left(L\mathrm{\text{'}}\mathrm{\text{'}}-\mathrm{2}z\mathrm{t}\mathrm{a}\mathrm{n}{\beta }_{\mathrm{2}}\right)}^{n_{\mathrm{2}}}}\mathrm{d}z]$$

（11）

式中：$n_{\mathrm{1}}=\left(\mathrm{2}-\frac{\mathrm{4}}{\mathrm{\pi }}\right)\frac{K_{\mathrm{0}}\mathrm{c}\mathrm{o}\mathrm{s}{\beta }_{\mathrm{1}}\mathrm{t}\mathrm{a}\mathrm{n}\phi }{\mathrm{t}\mathrm{a}\mathrm{n}{\beta }_{\mathrm{1}}}$；$n_{\mathrm{2}}=\frac{\mathrm{4}}{\mathrm{\pi }}\frac{K_{\mathrm{0}}\mathrm{c}\mathrm{o}\mathrm{s}{\beta }_{\mathrm{2}}\mathrm{t}\mathrm{a}\mathrm{n}\phi }{\mathrm{t}\mathrm{a}\mathrm{n}{\beta }_{\mathrm{2}}}$由边界条件可知，当z=0时，则基于椭球体理论的掌子面上覆土压力理论计算式为公式（12）：

$${\sigma }_{\mathrm{v}}={\left(B\mathrm{\text{'}}\mathrm{\text{'}}-\mathrm{2}z\mathrm{t}\mathrm{a}\mathrm{n}{\beta }_{\mathrm{1}}\right)}^{n_{\mathrm{1}}}{\left(L\mathrm{\text{'}}\mathrm{\text{'}}+\mathrm{2}z\mathrm{t}\mathrm{a}\mathrm{n}{\beta }_{\mathrm{2}}\right)}^{n_{\mathrm{2}}}[{\int }_{\mathrm{0}}^z\frac{\gamma -\frac{\mathrm{4}\left(\mathrm{\pi }-\mathrm{2}\right)c}{\pi \left(B\mathrm{\text{'}}\mathrm{\text{'}}-\mathrm{2}z\mathrm{t}\mathrm{a}\mathrm{n}{\beta }_{\mathrm{1}}\right)}-\frac{\mathrm{8}c}{\pi \left(L\mathrm{\text{'}}\mathrm{\text{'}}-\mathrm{2}z\mathrm{t}\mathrm{a}\mathrm{n}{\beta }_{\mathrm{2}}\right)}}{{\left(B\mathrm{\text{'}}\mathrm{\text{'}}-\mathrm{2}z\mathrm{t}\mathrm{a}\mathrm{n}{\beta }_{\mathrm{1}}\right)}^{n_{\mathrm{1}}}{\left(L\mathrm{\text{'}}\mathrm{\text{'}}-\mathrm{2}z\mathrm{t}\mathrm{a}\mathrm{n}{\beta }_{\mathrm{2}}\right)}^{n_{\mathrm{2}}}}dz+\frac{q_{\mathrm{0}}}{B\mathrm{\text{'}}{\mathrm{\text{'}}}^{n_{\mathrm{1}}}L\mathrm{\text{'}}{\mathrm{\text{'}}}^{n_{\mathrm{2}}}}]$$

（12）

对于砂卵石地层而言，其内部黏聚力c=0，则应用于砂卵石地层的盾构隧道上覆土压力计算公式如下：

$${\sigma }_{\mathrm{v}}={\left(B\mathrm{\text{'}}\mathrm{\text{'}}-\mathrm{2}z\mathrm{t}\mathrm{a}\mathrm{n}{\beta }_{\mathrm{1}}\right)}^{n_{\mathrm{1}}}{\left(L\mathrm{\text{'}}\mathrm{\text{'}}+\mathrm{2}z\mathrm{t}\mathrm{a}\mathrm{n}{\beta }_{\mathrm{2}}\right)}^{n_{\mathrm{2}}}[{\int }_{\mathrm{0}}^Z\frac{\gamma }{{\left(B\mathrm{\text{'}}\mathrm{\text{'}}-\mathrm{2}z\mathrm{t}\mathrm{a}\mathrm{n}{\beta }_{\mathrm{1}}\right)}^{n_{\mathrm{1}}}{\left(L\mathrm{\text{'}}\mathrm{\text{'}}+\mathrm{2}z\mathrm{t}\mathrm{a}\mathrm{n}{\beta }_{\mathrm{2}}\right)}^{n_{\mathrm{2}}}}\mathrm{d}z+\frac{q_{\mathrm{0}}}{B\mathrm{\text{'}}{\mathrm{\text{'}}}^{n_{\mathrm{1}}}L\mathrm{\text{'}}{\mathrm{\text{'}}}^{n_{\mathrm{2}}}}]$$

（13）

式中：作用于椭圆台台顶面的土压力为$q_{\mathrm{0}}=\frac{\mathrm{2}}{\mathrm{3}}a\gamma$， 其中$\gamma$为土体的单位重度。

2.3　上覆围岩压力与隧道埋深的关系

在隧道埋深较浅时，掌子面前方的上部土体往往不能发展成完整椭球体，因此埋深对于椭球体的构建具有直接的关系，下面对几种不同埋深情况下的椭球发育情况与对应的掌子面上覆土压力进行介绍。图4展示了几种不同埋深情况下的椭球体发育程度，其中的隧道拱顶至坍落椭球体短轴的高度Z与上部塌落椭球体高度h可由2.1节所述公式求得，隧道上部松动土压力与椭球体发育程度直接相关，当盾构埋深Q小于椭球体的高度h时，视为将椭球体延地平线切割后剩余的部分作用在掌子面上方。

图4  隧道埋深与椭球体发育高度的关系（Q≥h）

Fig. 4  Relationship between tunnel buried depth and ellipsoid developing height(Q≥h)

（1） 隧道埋置深度（Q≥h）

如图4所示，埋深较大（Q≥h）时，塌落椭球体完全发育，当隧道埋深Q等于发育椭球体高度h时，地层内部的椭球体正好处于发育完成的临界状态，将此时临界发育状态的椭球体高度定义为砂卵石地层隧道椭球体发育的临界埋深。此时上覆松动土压力计算公式参照式（13），式中$z=Z, q_{\mathrm{0}}=\frac{\mathrm{2}}{\mathrm{3}}a\gamma$。

（2） 隧道埋置深度（Z<Q<h）

如图5所示，在一定埋深（Z<Q<h）时，椭球体部分发育。此时上覆松动土体分为两个部分，其中大于Z的部分，按全土柱理论进行计算；小于Z的部分，按公式（13）进行计算，式中$z=Z, q_{\mathrm{0}}=\gamma (Q-Z)+ q_{\mathrm{1}}$。

图5  隧道埋深与椭球体发育高度的关系（Z<Q<h）

Fig. 5  Relationship between tunnel buried depth and ellipsoid developing height(Z<Q<h)

（3） 隧道埋置深度（Q≤Z）

如图6所示，在埋深更浅（Q≤Z）时，上覆塌落椭球体较少部分发育。此时上覆松动土压力按公式（13）进行计算，式中$z=Q, q_{\mathrm{0}}=q_{\mathrm{1}}$。其中：q₁为地表荷载。

图6  隧道埋深与椭球体发育高度的关系（Q≤Z）

Fig. 6  Relationship between tunnel buried depth and ellipsoid developing height(Q≤Z)

2.4　掌子面稳定极限支护力求解

结合极限平衡理论建立掌子面前方滑移体整体力矩平衡，选取对数螺旋线滑移面的中心点O作为力矩平衡的中心点，根据滑移体上所受荷载关于中心点O的力矩之和为0，进行掌子面极限支护力的求解。由图2可知，掌子面前方滑移体作用力矩有：上覆椭球体松动区自重应力作用力矩M_v、滑动土体自重力矩M_g、两侧摩阻力矩M_ts、滑动面切向摩阻力矩M_nt、掌子面支护力作用力矩M_t。下面对各个力矩进行推导与计算。

（1） 上覆椭球体松动区自重应力作用力矩M_v

根据计算模型的几何条件，上覆荷载的作用力臂为

$$l_{\mathrm{g}}=l_{\mathrm{a}}+l_{\mathrm{j}}=\frac{r_{\mathrm{a}}\mathrm{c}\mathrm{o}\mathrm{s}({\theta }_{\mathrm{a}}+\phi )+r_{\mathrm{0}}\mathrm{c}\mathrm{o}\mathrm{s}\phi }{\mathrm{2}}$$

（14）

则顶部荷载产生的力矩为

$$M_{\mathrm{v}}=\frac{\mathrm{\pi }B{l}_{\mathrm{j}}}{\mathrm{2}}{\sigma }_{\mathrm{v}}{l}_{\mathrm{g}}={\sigma }_{\mathrm{v}}B\frac{(r_{\mathrm{a}}\mathrm{c}\mathrm{o}\mathrm{s}({\theta }_{\mathrm{a}}{+\phi ))}^{\mathrm{2}}-{\left(r_{\mathrm{0}}\mathrm{c}\mathrm{o}\mathrm{s}\phi \right)}^{\mathrm{2}}}{\mathrm{4}}$$

（15）

（2） 滑动土体自重力矩M_g

将对数螺旋滑动区域土体沿竖向进行微分，如图7所示，其中

$$\mathrm{d}g=\gamma \left[r\mathrm{c}\mathrm{o}\mathrm{s}\left(\phi +\theta \right)-l_{\mathrm{a}}\right]B\mathrm{d}z$$

（16）

$$\mathrm{d}z=r\mathrm{s}\mathrm{i}\mathrm{n}\alpha \mathrm{d}\theta$$

（17）

$$\alpha =\frac{\mathrm{\pi }}{\mathrm{2}}-\theta$$

（18）

式中：α为滑动土体切线方向与水平面夹角。

图7  滑动土体微分土条

Fig. 7  Differential soil slice of sliding soil

微分土条的单位弯矩可以表示为

$$\mathrm{d}{M}_{\mathrm{g}}=\frac{r\mathrm{c}\mathrm{o}\mathrm{s}\left(\phi +\theta \right)+l_{\mathrm{a}}}{\mathrm{2}}\mathrm{d}g$$

（19）

对式（19）积分，可得滑动区土体沿着圆心O的作用力矩为

$$\begin{array}{l}{M}_{\mathrm{g}}={\int }_{\mathrm{0}}^{{\theta }_{\mathrm{a}}}\gamma \left[r\mathrm{c}\mathrm{o}\mathrm{s}\left(\phi +\theta \right)-l_{\mathrm{a}}\right]Br\mathrm{s}\mathrm{i}\mathrm{n}(\frac{\mathrm{\pi }}{\mathrm{2}}-\theta )\frac{r\mathrm{c}\mathrm{o}\mathrm{s}\left(\phi +\theta \right)+l_{\mathrm{a}}}{\mathrm{2}}\mathrm{d}\theta \\ ={\int }_{\mathrm{0}}^{{\theta }_{\mathrm{a}}}\gamma \left[r\mathrm{c}\mathrm{o}\mathrm{s}\left(\phi +\theta \right)-r_{\mathrm{a}}\mathrm{c}\mathrm{o}\mathrm{s}({\theta }_{\mathrm{a}}+\phi )\right]Br\mathrm{s}\mathrm{i}\mathrm{n}(\frac{\mathrm{\pi }}{\mathrm{2}}-\theta )\frac{r\mathrm{c}\mathrm{o}\mathrm{s}\left(\phi +\theta \right)+r_{\mathrm{a}}\mathrm{c}\mathrm{o}\mathrm{s}({\theta }_{\mathrm{a}}+\phi )}{\mathrm{2}}\mathrm{d}\theta \end{array}$$

（20）

（3） 滑动面切向摩阻力矩M_nt

单位土条上沿着对数螺旋线滑移面上的切向力矩微分量dM_nt可由滑移面上的法向力微分量dT与剪切力微分量dN来表示：

$$\mathrm{d}{M}_{\mathrm{n}\mathrm{t}}=r\mathrm{c}\mathrm{o}\mathrm{s}\phi \mathrm{d}T-r\mathrm{s}\mathrm{i}\mathrm{n}\phi \mathrm{d}N$$

（21）

根据摩尔库伦强度准则：

$$\mathrm{d}T=\mathrm{t}\mathrm{a}\mathrm{n}\phi \mathrm{d}N+crB\mathrm{d}\theta$$

（22）

则求得对数螺旋滑移面上的切向力矩为

$$M_{\mathrm{n}\mathrm{t}}=\frac{c{r}^{\mathrm{2}}B({\mathrm{e}}^{\mathrm{2}{\theta }_{\mathrm{a}}\mathrm{t}\mathrm{a}\mathrm{n}\phi }-\mathrm{1})\mathrm{c}\mathrm{o}\mathrm{s}\phi }{\mathrm{2}\mathrm{t}\mathrm{a}\mathrm{n}\phi }$$

（23）

（4） 两侧摩阻力矩M_ts

掌子面前方滑动土体的水平土拱效应是根据两侧摩阻力矩来表征的，这里假设两侧竖向滑移面上的垂直土压力随深度呈线性变化，根据图7，微分土条上的剪切力微分量dT_s表示为

$$\mathrm{d}{T}_{\mathrm{s}}=\left[K_{\mathrm{a}}\left({\sigma }_{\mathrm{\nu }}+\gamma z\right)\mathrm{t}\mathrm{a}\mathrm{n}\phi +c\right]\left[r\mathrm{c}\mathrm{o}\mathrm{s}\left(\phi +\theta \right)-l_{\mathrm{a}}\right]\mathrm{d}z$$

（24）

其中：

$$z=r\mathrm{s}\mathrm{i}\mathrm{n}\left(\phi +\theta \right)-r_{\mathrm{0}}\mathrm{s}\mathrm{i}\mathrm{n}\phi$$

（25）

$$K_{\mathrm{a}}=\mathrm{t}\mathrm{a}{\mathrm{n}}^{\mathrm{2}}\left(\frac{\mathrm{\pi }}{\mathrm{4}}-\frac{\phi }{\mathrm{2}}\right)$$

（26）

式中：K_a为土体的主动土压力系数。

因此，微分土条上的竖向剪切力矩微分量dM_ts为

$$\mathrm{d}{M}_{\mathrm{t}\mathrm{s}}=\left[r\mathrm{c}\mathrm{o}\mathrm{s}\phi -\frac{r\mathrm{c}\mathrm{o}\mathrm{s}\left(\phi +\theta \right)-l_{\mathrm{a}}}{\mathrm{2}}\mathrm{c}\mathrm{o}\mathrm{s}\theta \right]\mathrm{d}T$$

（27）

对式（27）积分可得式（28）：

$$\begin{array}{l}{M}_{\mathrm{t}\mathrm{s}}={\int }_{\mathrm{0}}^{{\theta }_{\mathrm{a}}}\left\{K_{\mathrm{a}}\left[{\sigma }_{\mathrm{\nu }}+\gamma r\mathrm{s}\mathrm{i}\mathrm{n}(\phi +\theta )-\gamma r_{\mathrm{0}}\mathrm{s}\mathrm{i}\mathrm{n}\phi \right]\mathrm{t}\mathrm{a}\mathrm{n}\phi +c\right\}\cdot \left[r\mathrm{c}\mathrm{o}\mathrm{s}(\phi +\theta )-l_{\mathrm{a}}\right]\cdot r\mathrm{s}\mathrm{i}\mathrm{n}(\frac{\pi }{\mathrm{2}}-\theta )\cdot \\ \left[r\mathrm{c}\mathrm{o}\mathrm{s}\phi -\frac{r\mathrm{c}\mathrm{o}\mathrm{s}(\phi +\theta )-l_{\mathrm{a}}}{\mathrm{2}}\mathrm{c}\mathrm{o}\mathrm{s}\theta \right]\mathrm{d}\theta \end{array}$$

（28）

（5） 掌子面支护力作用力矩M_t

根据几何关系，掌子面水平支护力的作用力臂为

$$l_{\mathrm{s}}=r_{\mathrm{0}}\mathrm{c}\mathrm{o}\mathrm{s}\phi +D/\mathrm{2}$$

（29）

因此水平支护力在整个滑动土体上的作用力矩为

$$M_{\mathrm{t}}=\mathrm{\pi }{\sigma }_{\mathrm{t}}{D}^{\mathrm{2}}{l}_{\mathrm{s}}/\mathrm{4}=\mathrm{\pi }{\sigma }_{\mathrm{t}}{D}^{\mathrm{2}}\left(r_{\mathrm{0}}\mathrm{c}\mathrm{o}\mathrm{s}\phi +D/\mathrm{2}\right)/\mathrm{4}$$

（30）

（6） 水平极限支护力σ_t

根据式（14）—（30）求出各项力矩之后，根据滑移体所受荷载关于O点力矩为0，可得：

$$M_{\mathrm{v}}+M_{\mathrm{g}}-M_{\mathrm{n}\mathrm{t}}-\mathrm{2}{M}_{\mathrm{t}\mathrm{s}}-M_{\mathrm{t}}=\mathrm{0}$$

（31）

将式（15）、式（20）、式（23）、式（28）和式（30）各项代入式（31）并化简，得到掌子面水平极限支护力解析解如式（32）所示。根据前文推导过程，可以发现掌子面极限支护力是关于隧道直径D，土体内摩擦角φ以及破坏区土体所受外力的函数。

$${\sigma }_{\mathrm{t}}=\frac{\mathrm{4}{M}_{\mathrm{q}}+\mathrm{4}{M}_{\mathrm{g}}-\mathrm{4}{M}_{\mathrm{n}\mathrm{t}}-\mathrm{8}{M}_{\mathrm{t}\mathrm{s}}}{\mathrm{\pi }{D}^{\mathrm{2}}\left(r_{\mathrm{0}}\mathrm{c}\mathrm{o}\mathrm{s}\phi +D/\mathrm{2}\right)}$$

（32）

3 公式验证

3.1 数值模拟方法选取

常见数值模拟方法有有限元模拟、有限差分模拟以及离散元模拟。有限元法常用于处理连续介质地层的变形问题，对于砂卵石这一离散性介质并不适用；有限差分法根据其计算原理，可以解决部分离散介质的地层变形问题，但并不能反映地层中的力学响应关系；离散元法与前面两种方法有所不同，其建模计算体系是基于颗粒接触理论，能够对砂卵石地层的“点对点”接触模式进行较为准确的拟合，可以较好地解决非连续介质的变形问题。因此本节采用颗粒离散元程序（PFC3D）对所推公式进行验证。

由于砂卵石地层中，卵石作为主要的支撑骨架，砂土与卵石粒径级配差距过大，因此在模拟过程中，为提高计算效率，不考虑砂土颗粒的影响，仅对卵石颗粒粒径进行标定。

3.2 数值模型参数确定

本研究基于西南地区某盾构隧道，结合现场的实际工程概况，如图8所示，取定模型的几何尺寸为25 m×25 m×10 m（X×Y×Z），盾构隧道直径6.48 m，埋深15~22 m，这里取15 m，离散元颗粒直径参照实际卵石粒径范围：0.17 ~0.35 m，为更好地体现掌子面前方地层颗粒的流动状态，将隧道周围的颗粒尺寸粒径范围控制为0.17 ~0.27 m，卵石的内摩擦角φ=40.0°，模型颗粒总数为83 660。在模型的底面与四周施加wall单元约束，不允许发生位移；地表面作为自由面允许产生竖向变形；盾构机与管片均采用wall单元模拟，卵石颗粒采用ball单元，卵石和墙面以及卵石之间的接触本构采用线性接触，模型细观参数设置借鉴前人的经验^［20］。

图8  隧道开挖离散元数值模型

Fig. 8  Discrete element numerical model of tunnel excavation

3.3　结果对比

采用地表荷载为0 的工况，通过在掌子面处施加相应的梯形荷载来模拟支护力的施加（因为土压力呈梯形分布），并用掌子面支护应力比λ（掌子面中心位置处支护压力F与静止土压力σ₀的比值）来表征支护力F的大小，几个不同支护力情况下的掌子面土体变形情况如下：

由图9可知，当支护应力比λ=1.0时，掌子面支护力与原始地层静止土压力相等，土体变形最小，处于相对静止状态。

图9  支护应力比λ＝1.0

Fig. 9  Supporting stress ratio at λ＝1.0

由图10可知，当λ＜1.0时，随着支护力的减小，土体变形逐渐增大直至整体失稳，掌子面前方土体先在地层土压力的作用下挤出，并逐渐发展至隧道拱顶上方地层，形成图中的椭圆松动区域。

图10  支护应力比λ＝0.2

Fig. 10  Supporting stress ratio at λ＝0.2

如图11所示，当λ＞1.0时，掌子面支护力大于原始地层静止土压力，掌子面土体在支护力的作用下朝隧道开挖方向运动，此时前方土体的破坏模式近似于半圆形，随着支护力的增大，该区域会逐渐增大并扩展至地层深处引起地层被动失稳破坏。

图11  支护应力比λ＝1.4

Fig. 11  Supporting stress ratio at λ＝1.4

进一步，根据位移云图记录在不同支护应力比条件下的掌子面最大挤出位移，建立支护应力比与掌子面最大水平位移的关系曲线（图12，图中σ₀为初始压力，盾构隧道开挖直径6.48 m，土体黏聚力为0，土体重度为25 kg·m^-³），寻找掌子面支护力轻微变化引起水平位移陡增的点作为极限支护力点。

图12  应力比与最大水平位移关系曲线

Fig. 12  Relationship curve between stress ratio and maximum horizontal displacement

当掌子面支护应力比为0.17时，掌子面挤出位移变化最大，其对应的极限支护力为41.9 kPa，将相关参数代入式（32）中，算得掌子面极限支护力为40.6 kPa，两者的误差为3.1 %。故认为本文提出的对数螺旋椭球体模型可以对砂卵石地层的掌子面极限支护力进行较为精确的求解。

4 结论

为研究砂卵石地层盾构隧道掌子面稳定所需的支护力，结合椭球体理论与极限平衡法的对数螺旋线模型，推导了砂卵石地层盾构隧道的掌子面极限支护力理论解析解，主要得到以下结论：

（1）基于极限平衡理论，将掌子面前方失稳土体作为研究对象，考虑滑移体的整体力矩平衡，作用在滑移体上的力矩有上覆椭球体松动区自重应力作用力矩、滑动土体自重力矩、两侧摩阻力矩、滑动面切向摩阻力矩、掌子面支护力矩。

（2）掌子面前方土体在支护力的作用下发生变形。其中，支护应力比λ=1.0时土体变形最小；λ＜1.0时，随着支护力的减小，土体变形逐渐增大直至整体失稳；λ＞1.0时，掌子面土体在支护力的作用下朝隧道开挖方向运动。

（3）对比砂卵石地层盾构隧道掌子面极限支护力理论解析解与PFC3D软件数值模拟的结果，发现二者误差为3.1 %，较好地证明了推导公式的准确性及可行性。

砂卵石地层中进行隧道线形设计时，应特别考虑砂卵石地层综合地质情况，以此来尽量避免隧道穿越不良地层。在砂卵石地层中遇到盾构机上坡的情况，尤其要控制好盾构机土舱压力以及单环管片拼装过程中的出土量。

作者贡献声明

王立新：提出研究主题。

胡瑞青：论文框架搭设及背景调研。

任超：论文构思及撰写。

张才飞：论文图像信息处理。

张俊元：论文修订。

姜寅：论文修订。

刘畅：理论推导及数值模拟。