随着城市化建设高速发展，开发地下空间资源已成为构建资源节约型城市的重要手段之一^[1]，如依靠地铁建设缓解地上交通拥挤的情况。盾构法作为隧道施工的工法之一，以其施工速度快、周围环境不受干扰等优点而被广泛采用^[2]。盾构施工常常遇见复合地层的情况，盾构在复合地层掘进过程中若土仓压力控制不当，可能导致开挖面失稳，严重破坏周边建筑物，并严重威胁施工人员生命安全。由于地下岩土体特性及受力情况的复杂，关于盾构开挖面稳定性的研究主要集中于数值模拟与理论解析研究方面。数值模拟的优势在于能够模拟盾构施工中复杂的地层情况及各类荷载，可为实际工程评估提供借鉴。Senent等^[3]利用FLAC^3D对上软下硬复合地层的极限支护压力与开挖面的软硬比进行研究，得出上部地层较软时，支护压力不足会导致开挖面局部破坏的结论。何小辉等^[4]利用ABAQUS对上软下硬复合地层进行模拟，研究了隧道的变形特征。Zhang等^[5]基于网格自适应有限元法，研究了上软下硬地层中矩形隧道开挖面的稳定性。理论解析法大致可分为极限平衡法和极限分析法。已有研究多基于极限平衡法。陈强^[6]推导出上软下硬复合地层中极限支护压力的计算公式。李康^[7]引入了地层复合角的概念，研究了上软下硬复合地层中盾构隧道施工引起的地表沉降及开挖面支护力。同时，极限分析因其严密的理论框架也被部分学者采用^[8-9]。吕玺琳等^[10]采用极限分析法推导了2维村山氏模型盾构开挖面极限支护压力最小值的计算公式。冯利坡等^[11]采用极限分析上限法建立盾构隧道开挖面3维对数螺旋破面模型，并推导出极限支护压力计算公式。在极限分析的框架内，部分学者对盾构隧道上覆土体成层情况及土体黏聚力非均质的情况进行了研究^[12-14]。

可见，对于盾构穿越复合地层开挖面稳定性解析方法的研究仍然较少。解析方法可为参数研究提供极大便利，因此，有必要对盾构穿越复合地层开挖面稳定性理论开展进一步研究。

支护压力计算方面，俞超杰等^[15]运用FLAC^3D软件对不同泥浆液位下盾构隧道开挖面主被动极限支护力进行研究，得出盾构推进过程中支护力非均匀分布与均匀分布的计算结果十分接近，故本文将盾构开挖面上的支护力简化为均布荷载。同时，学者们给出不同的破坏模式。Jancsecz等^[16]采用楔形体极限平衡模型，计算最小极限支护力；Soubra^[17]采用多块体破坏模型并运用极限分析得到了极限支护力上限解；Takano等^[18]通过X射线断层扫描仪得出地层竖直破坏面类似于对数螺旋，故本文采用对数螺线破坏面作为复合地层盾构隧道开挖面破坏模式。

作者基于极限分析上限定理，构建了复合地层中盾构隧道开挖面对数螺旋破坏模式；将盾构开挖面上的支护力简化为均布荷载，计算了外力功率与内能耗损率，并利用非线性优化方法求得极限支护力的最优解，且与现有文献进行对比验证。着重讨论了地层参数、土体重度、隧道直径及深度系数对极限支护力及破坏模式的影响，并给出临界滑动面示意图。本文推导出的复合地层开挖面极限支护力解答简便而高效，为实际工程提供了一定参考。

1 计算原理和计算模型 1.1 计算原理

极限分析上限定理可表述为对于任意一个运动学许可的破坏机构，令其外力所做功率与内能耗散率相等，便能够得到破坏荷载或极限荷载的一个不安全的上限，所确定的荷载不小于真实极限荷载^[19]，即：

$\int_S {{{{T}}_i}} {{v}}_i^*{\rm{d}}S + \int_V {{{{F}}_i}} {{v}}_i^*{\rm{d}}V = \int_V {{{\sigma}} _{ij}^*} \dot {{\varepsilon}} _{ij}^*{\rm{d}}V$

(1)

式中，T_i和F_i分别为边界S上的面积分布力矢量（如支护力）和区域V内的体积力矢量（如土体自重）， $\dot{{ \varepsilon}} _{ij}^*$ 为运动许可速度场中的塑性应变率场， ${{\sigma}} _{ij}^*$ 为通过相关联流动法则与 $\dot {{\varepsilon }}_{ij}^*$ 关联的应力场， ${{v}}_i^*$ 为与 $\dot {{\varepsilon}} _{ij}^*$ 满足几何相容条件的速度场。

将土体自重与开挖面支护力做功功率作为外功率 $W$ ；同时，假设岩土体为理想弹塑性，服从摩尔库伦屈服准则，破坏时服从相关联流动法则。将土体内部能量耗散作为内功率D，开挖面上的最小极限支护力可由外功率等于内功率求得。同时，为简化支护力的计算，假设盾构推进过程中支护压力为均匀分布。

1.2 计算模型

如图1所示，根据相关联流动法则的要求^[20]，挖面前方塌落块体边界线由3条对数螺旋线AF、FE和EB组成，塌落块体AFEB以角速度 $\omega $ 绕点O作刚体旋转破坏，点O为极坐标原点。隧道埋深为C，隧道直径为D；盾构穿越上部地层高度为nD，下部地层为（1–n）D，其中n为深度系数。上下部分地层抗剪强度参数分别为 ${c_1}$ 、 ${\varphi _1}$ 和 ${c_{\rm{2}}}$ 、 ${\varphi _{\rm{2}}}$ 。

由此可得，对数螺旋线AF、FE和EB的表达式为：

$\left\{ \begin{array}{l} {r_1}\left( \theta \right) = {r_{OA}}{{\rm{e}}^{\left( {\theta - {\theta _A}} \right)\tan\; {\varphi _1}}}, \\ {r_{\rm{2}}}\left( \theta \right) = {r_{OE}}{{\rm{e}}^{\left( {{\theta _E} - \theta } \right)\tan \;{\varphi _1}}}, \\ {r_{\rm{3}}}\left( \theta \right) = {r_{OB}}{{\rm{e}}^{\left( {{\theta _B} - \theta } \right)\tan \;{\varphi _2}}} \end{array} \right.$

(2)

式中， ${r_{OA}}$ 、 ${r_{OB}}$ 、 ${r_{OE}}$ 分别为极径OA、OB、OE的长度， ${\theta _B}$ 、 ${\theta _E}$ 、 ${\theta _A}$ 分别为OB、OE、OA与竖直向的夹角。

2 外功率与内能耗散率 2.1 外功率

构建速度相容的破坏机构后，利用极限分析上限定理求解极限支护力的关键在于外功率和能耗的计算。本文所讨论的外力有土体重力和开挖面上的支护力。首先，计算土体重力功率，如图1所示，重力功率可由OEF、OEB区域功率之和减去OAF、OAB区域功率得到。考虑区域OEF，取一微元如图2所示。

微元面积 ${\rm{d}}A = {1 / 2}{r_2^2}{\rm{d}}\theta $ ，重心处的速度 ${v_G} = {2 / 3}\omega {r_2}$ ，因此，该微元土重做功功率为：

$\begin{aligned}[b] {\rm{d}}{W_{OEF}} =& \gamma {\rm{d}}A \cdot {v_G}\sin\; \theta = {\rm{ }}\gamma \cdot \frac{1}{2}{r_2^2}{\rm{d}}\theta \cdot \frac{2}{3}\omega {r_2}\sin\; \theta = \\ &\frac{{\rm{1}}}{3}\gamma \omega {r_2^3}\sin\; \theta {\rm{d}}\theta \end{aligned} $

(3)

在区域OEF上积分得：

$\begin{aligned}[b] {W_{OEF}} =& \int_{{\theta _E}}^{{\theta _F}} {\frac{2}{3}\omega {r_2}\gamma } \sin\; \theta \frac{1}{2}{r_2}^2{\rm{d}}\theta = \\ &\frac{1}{3}\gamma \omega \int_{{\theta _E}}^{{\theta _F}} {{r_2}^{\rm{3}}} \sin \;\theta {\rm{d}}\theta \end{aligned} $

(4)

相似地，可得区域OEB、OAF和OAB外力功率，叠加后得到总外力功率：

$\begin{aligned}[b] {W_\gamma } =& \frac{1}{3}\gamma \omega \int_{{\theta _E}}^{{\theta _F}} {{r_2}^{\rm{3}}} \sin \;\theta {\rm{d}}\theta {\rm{ + }}\frac{1}{3}\gamma \omega \int_{{\theta _B}}^{{\theta _E}} {{r_3^3}} \sin\; \theta {\rm{d}}\theta - \\ &\frac{1}{3}\gamma \omega \int_{{\theta _A}}^{{\theta _F}} {{r_{\rm{1}}}^{\rm{3}}} \sin\; \theta {\rm{d}}\theta - \gamma \omega \frac{2}{3}{r_{OA}}\sin\; {\theta _A}\frac{1}{2}D{r_{OA}} \cdot \\ &\sin\; {\theta _A} = \gamma \omega r_{OA}^3\left( {{f_1} + {f_2} - {f_3} - {f_4}} \right) \end{aligned} $

(5)

式中， $\gamma $ 为土体重度，f₁～f₄、 ${\theta _E}$ 及 ${\theta _F}$ 计算公式如下：

$\begin{aligned}[b] {f_1} =& \frac{{{{\sin }^3}{\theta _A}}}{{{{\sin }^3}{\theta _B}}}\frac{{{{\rm{e}}^{3\left( {{\theta _B} - {\theta _E}} \right)\tan \;{\varphi _2}}}}}{{{\rm{3}}\left( {{\rm{1 + 9ta}}{{\rm{n}}^2}{\varphi _1}} \right)}}\left[ {\cos \;{\theta _E} + 3\sin \;{\theta _E}\tan \;{\varphi _1}} \right. - \\ &{{\rm{e}}^{3\left( {{\theta _E} - {\theta _F}} \right)\tan\; {\varphi _1}}}\left. {\left. {\left( {\cos\; {\theta _F} + } \right.3\sin \;{\theta _F}\tan \;{\varphi _1}} \right)} \right]\\[-10pt] \end{aligned} $

(6)

$\begin{aligned}[b] {f_{\rm{2}}} =& \frac{{{{\sin }^3}{\theta _A}}}{{{{\sin }^3}{\theta _B}}}\frac{{\rm{1}}}{{{\rm{3}}\left( {{\rm{1 + 9ta}}{{\rm{n}}^2}{\varphi _2}} \right)}}\left[ {\left( {\cos\; {\theta _B} + 3\sin \;{\theta _B}\tan \;{\varphi _2}} \right) - } \right. \\ &\left. {{{\rm{e}}^{3\left( {{\theta _B} - {\theta _E}} \right)\tan \;{\varphi _2}}}\left( {\cos \;{\theta _E} + 3\sin\; {\theta _E}\tan \;{\varphi _2}} \right)} \right]\\[-12pt] \end{aligned} $

(7)

$\begin{aligned}[b] {f_{\rm{3}}} =& \frac{{\rm{1}}}{{{\rm{3}}\left( {{\rm{1 + 9ta}}{{\rm{n}}^2}{\varphi _1}} \right)}}\left[ {{{\rm{e}}^{3\left( {{\theta _F} - {\theta _A}} \right)\tan\; {\varphi _1}}}( {3\tan \;{\varphi _1}\sin \;{\theta _F} -} } \right. \\ &\;\;\;\;\;\;\;\cos\; {\theta _F} )\left. { - \left( {3\tan\; {\varphi _1}\sin \;{\theta _A} - \cos\; {\theta _A}} \right)} \right] \end{aligned} $

(8)

${f_{\rm{4}}} = \frac{{\sin \left( {{\theta _A} - {\theta _B}} \right){{\sin }^2}{\theta _A}}}{{{\rm{3}}\sin \;{\theta _B}}}$

(9)

由图1中的几何关系可得：

${r_{OA}}\sin \;{\theta _A} = {r_{OB}}\sin \;{\theta _B}$

(10)

$\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;D\sin \;{\theta _A} = {r_{OB}}\sin \left( {{\theta _A} - {\theta _B}} \right)$

(11)

$\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{r_{OB}}\cos\; {\theta _B} - {r_{OE}}\cos\; {\theta _E} = \left( {1 - n} \right)D$

(12)

在对数螺旋线BE上有：

${r_{OE}} = {r_{OB}}{{\rm{e}}^{\left( {{\theta _B} - {\theta _E}} \right)\tan \;{\varphi _{\rm{2}}}}}$

(13)

将式（13）代入式（12），并化简得：

$\cos \;{\theta _B} - \cos\; {\theta _E}{{\rm{e}}^{\left( {{\theta _B} - {\theta _E}} \right)\tan\; {\varphi _2}}} = \left( {1 - n} \right)\frac{{\sin \left( {{\theta _A} - {\theta _B}} \right)}}{{\sin\; {\theta _A}}}$

(14)

解此方程可求得 ${\theta _E}$ 。

此外：

$\left\{ \begin{array}{l} {r_{OF}} = {r_{OA}}{{\rm{e}}^{\left( {{\theta _F} - {\theta _A}} \right)\tan \;{\varphi _1}}}, \\ {r_{OF}} = {r_{OE}}{{\rm{e}}^{\left( {{\theta _E} - {\theta _F}} \right)\tan \;{\varphi _1}}} \end{array} \right.$

(15)

联立求解式（15）可得：

$\begin{aligned}[b] {\theta _F} =& \frac{{\rm{1}}}{{{\rm{2}}\tan\; {\varphi _1}}}\left[ {{\theta _A}\tan\; {\varphi _1} + {\theta _B}\tan\; {\varphi _2} + {\theta _E}} \right. \cdot \\ &\left. {\left. {\left( {\tan\; {\varphi _1} - } \right.\tan\; {\varphi _2}} \right) + \ln \left( {{{\sin\; {\theta _A}} / {\sin \;{\theta _B}}}} \right)} \right] \end{aligned} $

(16)

盾构推进过程中支护压力简化为均匀分布，如图3所示。

图3中，P为开挖面上一点，相应极角为 $\theta $ ，由几何关系得：

$OP = \frac{{{r_{OA}}\sin \;{\theta _A}}}{{\sin\; \theta }}$

(17)

P点处的速度大小为：

$v = \frac{{{r_{OA}}\sin\; {\theta _A}}}{{\sin \;\theta }}\omega $

(18)

则可得支护力做功功率为：

$\begin{aligned}[b] {W_T} =& - \int_{{\theta _B}}^{{\theta _A}} {{\sigma _T}{r_{OA}}\frac{{\sin\; {\theta _A}}}{{\sin\; \theta }}\cos\; \theta \omega r} \frac{{{\rm{d}}\theta }}{{\sin\; \theta }} = \\ & - \frac{{\rm{1}}}{{\rm{2}}}\omega {\sigma _T}r_{OA}^2\left( {\frac{{{{\sin }^2}{\theta _A}}}{{{{\sin }^{\rm{2}}}{\theta _B}}} - 1} \right) \end{aligned} $

(19)

则总的外力功率为：

$ {W}_{{\text{总}}}={W}_{\gamma }+{W}_{T}$

(20)

2.2 内能耗散率

由于塌落块体作刚体旋转破坏，土体内部的能量耗散率忽略不计。将各速度间断面上的微分面积与速度间断线切线方向速度及土体黏聚力相乘，并沿整个间断面积分，可得AF、FE、BE上的内能耗损率：

$ {D_{AF}} = \int_{{\theta _A}}^{{\theta _F}} {{c_1}} \cdot \omega {r_1}\cos\; {\varphi _1} \cdot \frac{{{r_1}{\rm{d}}\theta }}{{\cos \;{\varphi _1}}} = \frac{{{c_1}\omega r_{OA}^2}}{{2\tan\; {\varphi _1}}}{g_1} $

(21)

$ {D_{FE}} = \int_{{\theta _E}}^{{\theta _F}} {{c_1}} \cdot \omega {r_2}\cos\; {\varphi _1} \cdot \frac{{{r_2}{\rm{d}}\theta }}{{\cos\; {\varphi _1}}} = \frac{{{c_1}\omega r_{OA}^2}}{{2\tan\; {\varphi _1}}}{g_2} $

(22)

$ {D_{BE}} = \int_{{\theta _B}}^{{\theta _E}} {{c_2}} \cdot \omega {r_3}\cos\; {\varphi _2} \cdot \frac{{{r_3}{\rm{d}}\theta }}{{\cos \;{\varphi _2}}} = \frac{{{c_2}\omega r_{OA}^2}}{{2\tan \;{\varphi _2}}}{g_3} $

(23)

式中，g₁～g₃可表示为：

${g_1}{\rm{ = }}{{\rm{e}}^{2\left( {{\theta _F} - {\theta _A}} \right)\tan\; {\varphi _1}}} - 1$

(24)

$ {g_2}{\rm{ = }}\frac{{{{\sin }^2}{\theta _A}}}{{{{\sin }^2}{\theta _B}}}{{\rm{e}}^{2\left( {{\theta _B} - {\theta _E}} \right)\tan\; {\varphi _2}}}\left[ {1 - {{\rm{e}}^{2\left( {{\theta _E} - {\theta _F}} \right)\tan\; {\varphi _1}}}} \right] $

(25)

${g_3}{\rm{ = }}\frac{{{{\sin }^2}{\theta _A}}}{{{{\sin }^2}{\theta _B}}}\left[ {1 - {{\rm{e}}^{2\left( {{\theta _B} - {\theta _E}} \right)\tan\; {\varphi _2}}}} \right]$

(26)

则总的内能耗散率为：

$ {D}_{{\text{总}}}={D}_{AF}+{D}_{FE}+{D}_{BE}$

(27)

2.3 极限支护力的求解

令总外力功率与内能耗散率相等可得：

$\begin{aligned}[b] {\sigma _T}{\rm{ = }}&\frac{{{{\sin }^2}{\theta _B}}}{{{{\sin }^2}{\theta _A} - {{\sin }^2}{\theta _B}}}\left[ {{\rm{2}}\gamma D\frac{{\sin \;{\theta _B}}}{{\sin \left( {{\theta _A} - {\theta _B}} \right)}}} \right.\left( {{f_1} + {f_2}} \right. - \\ &\left. {\left. {{f_3} - {f_4}} \right) - \frac{{{c_1}}}{{\tan\; {\varphi _1}}}\left( {{g_1} + {g_2}} \right) - \frac{{{c_2}}}{{\tan\; {\varphi _2}}}{g_3}} \right] \end{aligned} $

(28)

由于支护力做负功，上限解求其最大值，直接利用求导的方式很难求得支护力最大值。因此，本文利用序列二次规划法来优化求其最大值，约束条件如下：

$\left\{ \begin{array}{l} 0 < {\theta _B} < {\theta _A} < {\text{π} / 2}, \\ {\theta _A} < {\theta _E} < {\text{π}} , \\ {\theta _B} < {\theta _E} < {\theta _F} \end{array} \right.$

(29)

3 验　证

Cheng等^[21]利用随机极限分析法研究了复合地层盾构开挖面稳定性。将本文计算结果与其数据对比，如表1所示。其他参数取值如下： ${\varphi _1} = {\varphi _2} = {30^ {\circ} }$ ， $\gamma = 20\;{\rm{kN/}}{{\rm{m}}^{\rm{3}}}$ ， $n = 0.5$ ， $D = 10$ m，c₁和c₂分别为上部地层和下部地层的黏聚力大小。

由表1可知，利用本文解答优化得到的结果非常接近文献[21]中的结果，最大误差不超过5%，由此验证了本文方法的可行性。

4 参数分析

为了解各因素对于盾构穿越复合地层稳定性的影响，计算不同参数下的开挖面极限支护力，并进行了相应的探讨。参数取值如下： ${\varphi _1} = {\varphi _2} = {20^ {\circ} }$ ， ${c_1} = {c_2} = 10\;{\rm{kPa}}$ ， $D = 8\;{\rm{m}}$ ， $C/D = 1.5$ ， $\gamma = 20\;{\rm{kN/}}{{\rm{m}}^{\rm{3}}}$ 。

4.1 地层强度参数

图4为不同深度系数n下，极限支护力随地层内摩擦角的变化。

由图4可知：当上部地层内摩擦角一定，随着下部地层内摩擦角的增大，极限支护力显著降低，这是由于土体自稳能力得到了改善；且深度系数越小，这一趋势越明显。深度系数为0.5时，上部地层内摩擦角为10°，下部地层内摩擦角由10°变化到30°，极限支护力减小56%，这比传统内摩擦角取均值时的计算结果减小了33%。同样地，当下部地层内摩擦角不变，随着上部地层内摩擦角增大，极限支护力逐渐降低。

图5为不同深度系数下，土体黏聚力对极限支护力的影响。

由图5可知：当上部地层黏聚力一定时，随着下部地层黏聚力的增大，极限支护力呈线性减小趋势；同样地，深度系数越小，下降幅度越大。另外，当下部地层黏聚力一定时，随着上部地层黏聚力的增大，极限支护力逐渐减小。

4.2 隧道直径及土体重度

图6为隧道直径对于极限支护力影响。

由图6可知：随着隧道直径的增大，极限支护力显著增大。另外，当土体内摩擦角较小时，极限支护力的变化幅度将更大；当地层内摩擦角为10°时，隧道直径由5 m变化到10 m，极限支护力增加了248%。

图7为土体重度对于极限支护力的影响。

由图7可知：随着土体重度的增加，极限支护力逐渐增大；此外，当土体内摩擦角较小时，极限支护力变化幅度增大。

4.3 深度系数

图8和9为深度系数对极限支护力的影响。

由图8、9可知：深度系数越小，极限支护力对于地层抗剪强度参数变化越敏感；反之，深度系数越大，极限支护力随抗剪强度参数变化而变化的幅度越小。当两地层抗剪强度参数相同，退化为单一地层，深度系数的变化对极限支护力没有影响。

4.4 临界滑动面示意图

图10为单一地层与复合地层中临界破坏面示意图。

由图10（a）可知，当上部地层内摩擦角不变时，随着下部地层内摩擦角的增大，破坏块体的体积减小，下部地层破坏面向隧道面移动，即较大的摩擦角产生较小的破坏范围。由图10（b）可知，当上部地层黏聚力不变，下部地层黏聚力增大，破坏面稍向隧道面移动，但变化并不明显。

5 结　论

在极限分析上限定理的框架内，构建了盾构穿越复合地层开挖面稳定性分析模型；根据内外功率相等推导出开挖面极限支护力表达式，利用优化方法求得极限支护力的上限解；讨论了一系列参数对极限支护力的影响。主要结论如下：

1）当上部地层内摩擦角一定，随着下部地层内摩擦角的增大，极限支护力显著降低，且深度系数越小，这一趋势越明显。当上部地层黏聚力一定时，随着下部地层黏聚力的增大，极限支护力呈线性减小趋势。同样地，深度系数越小，下降幅度越大。

2）随着隧道直径的增大，极限支护力显著增大；随着土体重度的增加，极限支护力逐渐增大。

3）当上部地层内摩擦角不变时，随着下部地层内摩擦角的增大，破坏块体的体积减小，下部地层破坏面向隧道面移动，即较大的摩擦角产生较小的破坏范围。当上部地层黏聚力不变，下部地层黏聚力增大，破坏面稍向隧道面移动，但变化并不明显。