基于倾斜滑裂面假设的浅埋盾构隧道松动土压力研究

Research on Loosening Earth Pressure of Shallow Buried Shield Tunnel Based on the Assumption of Inclined Sliding Surface

摘要

随着地铁建设的快速发展，隧道上方松动土压力的计算成了隧道工程领域的研究热点，而隧道松动区的确定对松动土压力的计算至关重要。为简便而准确地确定松动区的范围，首先假设松动区的边界为倾斜直线，然后基于Mohr-Coulomb破坏准则并考虑隧道尺寸、埋深等因素的影响推导滑裂角计算式，进而确定隧道周围地层的松动区范围。最后结合松动区土体主应力轴旋转理论以及倾斜滑裂面的假设，推导出滑裂面上的侧压力系数以及隧道松动土压力的计算式。结果表明：滑裂角随土体内摩擦角的增大而增大，随隧道埋深比的增大先减小然后趋于稳定；大主应力轨迹线为圆弧时计算更为简便、更易推广，可采用圆弧形大主应力轨迹的假设计算隧道松动土压力；松动土压力计算式所得结果和实际的试验结果吻合较好，因此该计算式可应用于砂土地质条件下浅埋盾构隧道上方的松动土压力计算。

Abstract

With the rapid development of metro construction, the calculation of loosening earth pressure acting on tunnels has become a research hotspot in the field of tunnel engineering. Accurate prediction of the loosening zone around tunnels is crucial for determining the loosening earth pressure. To simplify and precisely estimate the extent of the loosening zone around tunnels, its boundaries were first postulate as two inclined straight lines. Then, a formula for the sliding angle was derived considering tunnel size, buried depth of the tunnel, and other conditions based on the Mohr-Coulomb failure criterion, and thus the loosening zone of strata around the tunnel could be determined. Finally, the lateral pressure coefficient on the sliding surface and the formula for calculating the loosening earth pressure above tunnel were deduced by integrating the theory of principal stress axis rotation in the loosening zone soil mass with the assumption of an inclined sliding surface. The study shows that the sliding angle increases with the increase of the internal friction angle of soil, and initially decreases with the increase in tunnel burial depth ratio, subsequently stabilizing. When the trajectory line of the major principal stress is assumed to be circular, the calculation is more straightforward and conducive to wider application. The results obtained from the analytic solution presented in this paper align well with the actual test results, indicating its applicability for the calculation of loosening earth pressure above the shield tunnel in sandy strata.

Graphical abstract

关键词

隧道工程 / 松动区 / 土压力 / 土拱效应 / 滑裂面 / 数值模拟

Key words

Tunnel engineering / Loosening zone / Earth pressure / Soil arching effect / Sliding surface / Numerical simulation

计算和确定隧道松动区范围和松动土压力的大小对于衬砌的设计和盾构隧道的开挖具有重要意义。对于超浅埋隧道，作用在衬砌结构上的土压力可以利用岩柱法计算^［1］，但是随着上覆土层厚度的增大，土拱效应的产生以及地层破坏模式的变化，岩柱法不再适用。因此，考虑土拱效应和松动区破坏模式对隧道松动土压力大小的影响成为众多学者关注和研究的热点。

Terzaghi^［2］利用活动门试验提出松动区边界为垂直型滑裂面的模型，并考虑土拱效应推导得出了松动土压力的计算公式。此后该解析解被学者和工程师推广应用于隧道和地下管涵上覆土压力的计算。但是，Terzaghi对于松动土压力计算过程中围岩破坏范围与侧压力系数的确定仍然有待论证。为此，许多专家学者针对目前存在的问题展开了研究和修正。Chen等^［3］基于松动区土体主应力轴旋转理论，从侧压力系数的角度对Terzaghi松动土压力计算公式进行了修正。Chen等^［4］基于有限元计算结果，讨论了Terzaghi公式中存在的隧道周围土体位移、土压力分布、坍塌土体宽度和滑动摩擦等问题，并在此基础上提出了一种可行的隧道衬砌竖向土压力计算方法。Liang等^［5］认为土拱区的应力分布实际上是不均匀的，然后从主应力轴旋转的角度提出大主应力轨迹线为圆弧形状，并对Terzaghi松动土压力公式进行修正。

但是，以上解析解都是在Terzaghi提出的理论破坏模型的基础上推导建立的，一些学者^［6-8］在后来的研究中发现，隧道上方的破坏区域边界并非垂直型，而是曲线型，并且受到隧道尺寸，埋深以及地质条件的影响。Huang等^［9］和Gong等^［10］基于椭球体理论提出隧道上部松动区域为椭圆形，进而推导出隧道松动土压力计算公式，但是2位学者假设隧道上方的应力分布为梯形分布并不符合实际^［11］，从而导致其计算结果和实测数据有一定的差异。王将等^［12］基于滑移线场理论，推导出隧道松动区边界为曲线的解析表达式，并考虑主应力轴旋转提出了适用于盾构隧道的松动土压力计算模型，但其对于曲线边界的求解过程过于繁琐且参数较多，暂时无法在实际工程中推广应用。

为了简化松动土压力的求解过程同时又尽可能地提高计算精度，部分学者将松动区曲线边界简化为斜直线型，并以此推导土压力计算式。于丽等^［13］结合黄土地区浅埋隧道围岩破坏规律，推导出其破坏模式的计算式，然后结合极限平衡理论推导计算黄土隧道的围岩压力。刘翔等^［14］基于变坡面浅埋偏压隧道的破坏受力情况，推导了倾斜滑裂面破坏模式下的隧道松动压力计算式。Shukla等^［15］假设滑裂面为斜直线并考虑土拱效应推导出松动土压力的解析解。李瑞林等^［16］基于倾斜滑裂面的破坏模式，考虑土拱效应以及破坏区域土体处于非极限状态，推导出隧道松动土压力的计算式。但是，于丽等和刘翔等学者均是基于极限平衡理论推导的松动土压力，并未考虑土拱效应。Shukla等考虑了土拱效应的影响，但没有对滑裂面倾斜角度给出合理的计算方法。李瑞林等虽然给出了倾斜角度的计算，但是其仅考虑土体参数对倾斜角的影响。

综上所述，目前针对浅埋盾构隧道松动土压力计算的研究多以Terzaghi提出的垂直型破坏模式为主，但是这种破坏模式随着隧道埋深的增加而不再适用。部分学者提出了倾斜和曲线滑裂面的破坏模式，但是在破坏模式的确定上并没有给出合理的理论依据以及在松动土压力的计算方面考虑参数太多导致计算过程复杂且没有考虑土拱效应对松动土压力的影响。

本文为使隧道松动土压力推导过程简单化同时又尽可能地提高计算求解的精度，首先假设隧道上方松动区边界为斜线型，然后基于Mohr-Coulomb破坏准则，考虑隧道的尺寸、隧道埋深以及地质条件的影响推导出倾斜角度的解析解，最后基于所提出的倾斜滑裂面破坏模式和土拱效应，分别考虑圆弧和悬链线2种主应力轨迹线的情况推导计算松动土压力。

1 隧道围岩破坏模式

在盾构隧道掘进过程中，隧道上方的土体会发生松动产生不均匀沉降，距离隧道越近，土体的沉降越大，随着距离隧道越远，土体沉降逐渐减小直至为零^［8，17］。Terzaghi通过活动门试验得出土体松动滑移区的边界为曲线，但是为简化松动土压力的求解过程，Terzaghi提出滑移面与水平面夹角为\mathrm{45}°+\phi /\mathrm{2}（\phi为土的内摩擦角）、然后垂直向上的松动破坏形式。众多学者从试验以及实际工程的角度对此展开研究得出：隧道施工过程中，其上方土体的滑移面为曲线，并且滑裂角的大小受到隧道的尺寸、埋深以及地质条件的影响。由Moussaei等^［8］、方焘等^［17］对砂土地层隧道破坏模式的研究可知，隧道在埋深较浅的情况下，地层的破坏模式近似垂直土柱，但是随着隧道埋深的增大以及土体密实度的变化，地层的最终破坏形式并非垂直土柱，而是类似“漏斗”的形状。在这种情况下，如果依然采用Terzaghi提出的垂直滑移面破坏形式，则会导致假设滑动面和实际松动区滑动面之间存在相当大的差距，从而导致对隧道结构上的设计荷载的预测偏低。因此，本文为了准确计算隧道松动土压力，假设地层的松动区边界为倾斜直线。

现有研究表明，对于圆形盾构隧道，其破坏从两边的起拱处开始向上延伸至地表^［17］。因此，假设圆形隧道破坏模式如图1所示。图中：a，b分别为隧道起拱处的2点；i，d分别为地表沉降零界点；p，e为点a，b在地表的投影点；g为隧道拱顶位置点；f为隧道中轴线与地表的交点；n，l为过隧道拱顶的切线与隧道起拱点切线的交点；曲线ad和bi为地层的实际破坏面，虚直线ad和bi为本文假设破坏面；\beta为假设破坏面和水平方向的夹角；H为隧道的埋深；D为隧道直径；o为隧道中心； L 为地表沉降零界点和隧道起拱处在水平方向的投影距离。此外，令\gamma为土的容重，{\sigma }_x和{\sigma }_y分别为水平和竖直方向原岩应力。

因图1中破坏模式是对称的，取右侧一半进行分析。对于矩形lefg区域的土体，隧道开挖导致其失去支撑，其重力通过摩擦力转移至三角形bie部分的土体上。相对于矩形lefg部分的土体，glb部分土体可忽略不计，故将矩形lefg部分的土体重力平均化，折算到三角形bie土体中并忽略矩形lefg部分土体对{\sigma }_x的影响。取隧道纵向单位长度土体为研究对象，则矩形lefg部分土重为\mathrm{0.5}\gamma HD，三角形bie部分土重为\mathrm{0.5}\gamma L\left(H+\mathrm{0.5}D\right)。由此得到三角形bie区域的地应力场为

{\sigma }_y^{\text{'}}={\sigma }_y\left[\mathrm{1}+\frac{HD\mathrm{t}\mathrm{a}\mathrm{n}\beta }{{\left(H+\mathrm{0.5}D\right)}^{\mathrm{2}}}\right]（1）

{\sigma }_x^{\text{'}}={\sigma }_x（2）

式中：{\sigma }_x^{\text{'}}和{\sigma }_y^{\text{'}}分别为考虑lefg部分附加应力作用下bie区域的水平和竖直方向应力；K_{\mathrm{0}}为静止侧压力系数。

{\sigma }_x可由下式求得。

{\sigma }_x=K_{\mathrm{0}}{\sigma }_y（3）

根据Mohr-Coulomb破坏准则，滑裂面最大的抗滑强度{\tau }_{\mathrm{m}\mathrm{a}\mathrm{x}}为

{\tau }_{\mathrm{m}\mathrm{a}\mathrm{x}}={\sigma }_{\mathrm{n}}\mathrm{t}\mathrm{a}\mathrm{n}\phi +c（4）

式中：{\sigma }_{\mathrm{n}}为破坏面上的正应力；c为土的黏聚力。

定义滑裂面安全度F为

F={\tau }_{\mathrm{m}\mathrm{a}\mathrm{x}}-{\tau }_{\mathrm{s}}（5）

式中：{\tau }_{\mathrm{s}}为滑裂面上的剪应力。

由式（4）和式（5）可得

F=({\sigma }_y^{\text{'}}\mathrm{c}\mathrm{o}{\mathrm{s}}^{\mathrm{2}}\beta +{\sigma }_x^{\text{'}}\mathrm{s}\mathrm{i}{\mathrm{n}}^{\mathrm{2}}\beta )\mathrm{t}\mathrm{a}\mathrm{n}\phi +c-({\sigma }_y^{\text{'}}-{\sigma }_x^{\text{'}}\mathrm{s}\mathrm{i}\mathrm{n}\beta )\mathrm{c}\mathrm{o}\mathrm{s}\beta（6）

显然，{\sigma }_y^{\text{'}}与\beta相关，当F对\beta的1阶导数为0时，即可得出最危险滑裂角（极值点）。令F^{\text{'}}=\mathrm{0}，则有

\mathrm{t}\mathrm{a}\mathrm{n}(\mathrm{2}\beta )=\frac{K\mathrm{t}\mathrm{a}\mathrm{n}\phi +K_{\mathrm{0}}-\mathrm{1}}{K+\left(\mathrm{1}-K_{\mathrm{0}}\right)\mathrm{t}\mathrm{a}\mathrm{n}\phi }（7）

其中，

K=\frac{HD}{{\left(H+\mathrm{0.5}D\right)}^{\mathrm{2}}}

式中：K为隧道的几何形状参数，其大小由隧道埋深比H/D决定。

由图1可知\beta \in (\mathrm{\pi }/\mathrm{4}+\phi /\mathrm{2},\mathrm{\pi }/\mathrm{2})，并根据正切函数的二倍角关系可得

\mathrm{t}\mathrm{a}\mathrm{n}\beta =-\frac{\mathrm{1}+\sqrt[]{\mathrm{1}+\mathrm{t}\mathrm{a}{\mathrm{n}}^{\mathrm{2}}(\mathrm{2}\beta )}}{\mathrm{t}\mathrm{a}\mathrm{n}(\mathrm{2}\beta )}（8）

根据式（8）即可计算不同埋深、不同直径隧道以及不同地质参数条件下的隧道上方地层破坏倾角。

2 隧道松动土压力公式

2.1　计算模型

Terzaghi通过活动门试验证实了土拱效应的存在，并提出当土体之间产生不均匀沉降或相对位移以及存在作为支撑的拱脚时，地层内部就会产生土拱效应。同理，当隧道开挖后，隧道上方地层下沉，土体产生松动，松动土体和两侧未松动土体形成1个滑裂面，受滑移面两侧土体挤压地层拱效应逐渐产生，松动土体的部分自重会通过滑移面上剪力传递到滑移面外侧土体上。

由于土拱效应的影响，滑移面内土体水平方向应力和竖直方向应力不再为定值，故侧压力系数也不再为常数，而在Terzaghi所提出的松动土压力计算公式中侧压力系数取值为1。Liang等^［5］通过数值模拟得出侧压力系数因破坏区位置的不同而有所差异，当活动门下移时，中轴线处的侧压力系数大于1.0，然后逐渐向两侧递减，其大小受到土体空间位置以及土拱效应发挥程度的影响。

基于以上分析，本文以砂土为研究对象，建立如下浅埋隧道松动土压力的计算模型，如图2所示。模型假设可见文献［1］。图中：z为地表以下某位置覆土厚度；B_z为该位置处的松动区计算宽度的一半；\theta为滑裂面与竖直方向的夹角；q为地面超载；{\sigma }_{\mathrm{v}x}为隧道松动区水平微分土条单元任一位置x处的竖向应力，x为微分土条上任意一点到隧道中轴线的距离；{\sigma }_{\mathrm{h}}为作用在隧道上的侧压力。

需注意的是，Terzaghi提出同一深度处竖向应力大小相同，而根据主应力轴旋转理论可知，实际的竖向应力属于凹型曲线分布。

2.2　大主应力偏转轨迹

Handy^［18］假设小主应力拱迹线为悬链线，研究了挡墙土压力的土拱效应。Paik等^［19］则假设小主应力轨迹为圆弧，分析了挡墙土压力的土拱效应。Harrop-Williams^［20］研究发现小主应力拱迹线为介于悬链线与圆弧之间的曲线。Xu等^［21］假定小主应力轨迹分别为圆弧、抛物线和悬链线3种情况，研究了矿场回填竖向应力。Vivancoi等^［22］通过试验及数值模拟表明，活动门上方土体中水平应力为大主应力并指出采用大主应力拱假设能够较好地描述活动门上方土体中的主应力偏转情况。因此，本文基于Harrop-Williams的方法，从大主应力拱角度来研究土拱效应下隧道上方的松动土压力。

图3为隧道松动区土体应力偏转轨迹分析示意图。图中：A为任一微分土条与滑裂面的交点；{\sigma }_{\mathrm{1}}，{\sigma }_{\mathrm{3}}分别为大、小主应力；s为大主应力轨迹以隧道中轴线为起点计算的弧长，\psi为该位置处主应力偏转角；\omega为滑裂面处大主应力的偏转角，\omega =\mathrm{\pi }/\mathrm{4}+\phi /\mathrm{2}+\theta，当x\in [\mathrm{0},B_z]时，\psi \in [\mathrm{0},\omega ]。

1）圆弧形轨迹线

鉴于Harrop-Williams提出主应力拱轨迹线应是介于悬链线与圆弧线之间的曲线，本文首先假设图3中的大主应力偏转轨迹为圆弧，则其半径R和偏转角\psi分别为

R=\frac{B_z}{\mathrm{s}\mathrm{i}\mathrm{n}\omega }（9）

\psi =\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{s}\mathrm{i}\mathrm{n}(\frac{x}{B_z}\mathrm{s}\mathrm{i}\mathrm{n}\omega )（10）

由图2可得（9）和（10）中B_z为

B_z=\mathrm{0.5}D+(H+\mathrm{0.5}D-z)\mathrm{t}\mathrm{a}\mathrm{n}\theta（11）

2）悬链线形轨迹线

当大主应力轨迹线为悬链线时，假设主应力偏转角的正切值和大主应力轨迹呈线性关系，则有

\mathrm{t}\mathrm{a}\mathrm{n}\psi =ks+b_{\mathrm{1}}（12）

式中：k，b₁为待求函数系数。

由图3可知当s=\mathrm{0}时，\mathrm{t}\mathrm{a}\mathrm{n}\psi =\mathrm{0}，故b_{\mathrm{1}}=\mathrm{0}。

根据图3所示，对偏转轨迹线微分可得

\frac{\mathrm{d}z}{\mathrm{d}x}=\mathrm{t}\mathrm{a}\mathrm{n}\psi（13）

{(\mathrm{d}x)}^{\mathrm{2}}+{(\mathrm{d}z)}^{\mathrm{2}}={(\mathrm{d}s)}^{\mathrm{2}}（14）

联立式（12）—式（14）并积分可得

x=\frac{\mathrm{1}}{k}\mathrm{a}\mathrm{r}\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{h}(ks)+b_{\mathrm{2}}（15）

式中：b₂为待求参数。由图3可知当x=\mathrm{0}时，s=\mathrm{0}时，故b_{\mathrm{2}}=\mathrm{0}。

联立式（12）和式（15）并将边界条件x=B_z，\psi =\omega代入，可得

\mathrm{t}\mathrm{a}\mathrm{n}\psi =\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{h}(\frac{x}{B_z}\mathrm{a}\mathrm{r}\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{h}(\mathrm{t}\mathrm{a}\mathrm{n}\omega ))（16）

2.3　松动土压力的求解

图4中的摩尔应力圆表示的是滑裂面A点（图3）处土的应力状态。图中：{\sigma }_{\mathrm{v}A}，{\sigma }_{\mathrm{h}A}，{\tau }_A分别为A点处的竖向应力、水平应力和剪应力。

根据图中的三角函数关系可得

{\sigma }_{\mathrm{v}A}={\sigma }_{\mathrm{1}}\mathrm{s}\mathrm{i}{\mathrm{n}}^{\mathrm{2}}\omega +{\sigma }_{\mathrm{3}}\mathrm{c}\mathrm{o}{\mathrm{s}}^{\mathrm{2}}\omega（17）

{\sigma }_{\mathrm{h}A}={\sigma }_{\mathrm{1}}\mathrm{c}\mathrm{o}{\mathrm{s}}^{\mathrm{2}}\omega +{\sigma }_{\mathrm{3}}\mathrm{s}\mathrm{i}{\mathrm{n}}^{\mathrm{2}}\omega（18）

同理，松动区内部任一位置x处的竖向应力为

{\sigma }_{\mathrm{v}x}={\sigma }_{\mathrm{3}}(\mathrm{c}\mathrm{o}{\mathrm{s}}^{\mathrm{2}}\psi +K_{\mathrm{p}}\mathrm{s}\mathrm{i}{\mathrm{n}}^{\mathrm{2}}\psi )（19）

其中，

K_{\mathrm{p}}=\frac{{\sigma }_{\mathrm{1}}}{{\sigma }_{\mathrm{3}}}

式中：K_{\mathrm{p}}为被动土压力系数。

根据式（19）可求得水平微分单元所受到的平均竖向应力{\overline{\sigma }}_{\mathrm{v}}为

{\overline{\sigma }}_{\mathrm{v}}=\frac{{\int }_{\mathrm{0}}^{\mathrm{2}{B}_z}{\sigma }_{\mathrm{v}x}\mathrm{d}x}{\mathrm{2}{B}_z}=\frac{{\sigma }_{\mathrm{1}}}{B_z}{\int }_{\mathrm{0}}^{B_z}(\mathrm{c}\mathrm{o}{\mathrm{s}}^{\mathrm{2}}\psi +K_{\mathrm{p}}\mathrm{s}\mathrm{i}{\mathrm{n}}^{\mathrm{2}}\psi )\mathrm{d}x（20）

1）大主应力轨迹线为圆弧

对式（20）积分可得平均竖向应力{\overline{\sigma }}_{\mathrm{v}}为

{\overline{\sigma }}_{\mathrm{v}}={\sigma }_{\mathrm{3}}\left[\frac{\mathrm{1}}{\mathrm{3}}(K_{\mathrm{p}}-\mathrm{1})\mathrm{s}\mathrm{i}{\mathrm{n}}^{\mathrm{2}}\omega +\mathrm{1}\right]（21）

土拱效应发挥后土中侧向土压力系数K_{\mathrm{1}}为

K_{\mathrm{1}}=\frac{{\sigma }_{\mathrm{h}A}}{{\overline{\sigma }}_{\mathrm{v}}}=\frac{\mathrm{3}(K_{\mathrm{p}}\mathrm{c}\mathrm{o}{\mathrm{s}}^{\mathrm{2}}\omega +\mathrm{s}\mathrm{i}{\mathrm{n}}^{\mathrm{2}}\omega )}{(K_{\mathrm{p}}-\mathrm{1})\mathrm{s}\mathrm{i}{\mathrm{n}}^{\mathrm{2}}\omega +\mathrm{3}}（22）

将{\sigma }_{\mathrm{v}x}和2种不同主应力旋转轨迹下的应力均值{\overline{\sigma }}_{\mathrm{v}}的比值m定义为应力分布系数，用其表征不同主应力旋转轨迹下上覆土体作用在隧道上的压力在水平方向上的变化规律。当大主应力轨迹线为圆弧时，应力分布系数m为

m=\frac{{\sigma }_{\mathrm{v}x}}{{\overline{\sigma }}_{\mathrm{v}}}=\frac{\mathrm{3}(K_{\mathrm{p}}\mathrm{s}\mathrm{i}{\mathrm{n}}^{\mathrm{2}}\psi +\mathrm{c}\mathrm{o}{\mathrm{s}}^{\mathrm{2}}\psi )}{(K_{\mathrm{p}}-\mathrm{1})\mathrm{s}\mathrm{i}{\mathrm{n}}^{\mathrm{2}}\omega +\mathrm{3}}（23）

同理可得大主应力轨迹线为圆弧时，{\sigma }_{\mathrm{v}\mathrm{A}}和{\overline{\sigma }}_{\mathrm{v}}的比值K_{\mathrm{2}}为

K_{\mathrm{2}}=\frac{{\sigma }_{\mathrm{v}A}}{{\overline{\sigma }}_{\mathrm{v}}}=\frac{\mathrm{3}(K_{\mathrm{p}}\mathrm{s}\mathrm{i}{\mathrm{n}}^{\mathrm{2}}\omega +\mathrm{c}\mathrm{o}{\mathrm{s}}^{\mathrm{2}}\omega )}{(K_{\mathrm{p}}-\mathrm{1})\mathrm{s}\mathrm{i}{\mathrm{n}}^{\mathrm{2}}\omega +\mathrm{3}}（24）

2）大主应力轨迹线为悬链线

对式（20）积分可得平均竖向应力{\overline{\sigma }}_{\mathrm{v}}为

\begin{array}{l}{\overline{\sigma }}_{\mathrm{v}}={\sigma }_{\mathrm{3}}\frac{(\mathrm{1}-K_{\mathrm{p}})\mathrm{t}\mathrm{a}\mathrm{n}\mathrm{h}(\mathrm{a}\mathrm{r}\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{h}(\mathrm{t}\mathrm{a}\mathrm{n}\omega ))+}{\mathrm{a}\mathrm{r}\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{h}(\mathrm{t}\mathrm{a}\mathrm{n}\omega )}\to \\ \leftarrow \frac{K_{\mathrm{p}}\mathrm{a}\mathrm{r}\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{h}(\mathrm{t}\mathrm{a}\mathrm{n}\omega )}{\mathrm{a}\mathrm{r}\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{h}(\mathrm{t}\mathrm{a}\mathrm{n}\omega )}\end{array}（25）

土拱效应发挥后土中侧向土压力系数K_{\mathrm{1}}为

\begin{array}{l}{K}_{\mathrm{1}}=\frac{{\sigma }_{\mathrm{h}A}}{{\overline{\sigma }}_{\mathrm{v}}}=\frac{(K_{\mathrm{p}}\mathrm{c}\mathrm{o}{\mathrm{s}}^{\mathrm{2}}\omega +}{(\mathrm{1}-K_{\mathrm{p}})\mathrm{t}\mathrm{a}\mathrm{n}\mathrm{h}(\mathrm{a}\mathrm{r}\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{h}(\mathrm{t}\mathrm{a}\mathrm{n}\omega ))+}\to \\ \leftarrow \frac{\mathrm{s}\mathrm{i}{\mathrm{n}}^{\mathrm{2}}\omega )\mathrm{a}\mathrm{r}\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{h}(\mathrm{t}\mathrm{a}\mathrm{n}\omega )}{K_{\mathrm{p}}\mathrm{a}\mathrm{r}\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{h}(\mathrm{t}\mathrm{a}\mathrm{n}\omega )}\end{array}（26）

应力分布系数m为

\begin{array}{l}m=\frac{{\sigma }_{\mathrm{v}x}}{{\overline{\sigma }}_{\mathrm{v}}}=\frac{(K_{\mathrm{p}}\mathrm{s}\mathrm{i}{\mathrm{n}}^{\mathrm{2}}\psi +}{(\mathrm{1}-K_{\mathrm{p}})\mathrm{t}\mathrm{a}\mathrm{n}\mathrm{h}(\mathrm{a}\mathrm{r}\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{h}(\mathrm{t}\mathrm{a}\mathrm{n}\omega ))+}\to \\ \leftarrow \frac{\mathrm{c}\mathrm{o}{\mathrm{s}}^{\mathrm{2}}\psi )\mathrm{a}\mathrm{r}\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{h}(\mathrm{t}\mathrm{a}\mathrm{n}\omega )}{K_{\mathrm{p}}\mathrm{a}\mathrm{r}\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{h}(\mathrm{t}\mathrm{a}\mathrm{n}\omega )}\end{array}（27）

同理可得大主应力轨迹线为悬链线时，{\sigma }_{\mathrm{v}A}和{\overline{\sigma }}_{\mathrm{v}}的比值K_{\mathrm{2}}为

\begin{array}{l}{K}_{\mathrm{2}}=\frac{{\sigma }_{\mathrm{v}A}}{{\overline{\sigma }}_{\mathrm{v}}}=\frac{(K_{\mathrm{p}}\mathrm{s}\mathrm{i}{\mathrm{n}}^{\mathrm{2}}\omega +}{(\mathrm{1}-K_{\mathrm{p}})\mathrm{t}\mathrm{a}\mathrm{n}\mathrm{h}(\mathrm{a}\mathrm{r}\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{h}(\mathrm{t}\mathrm{a}\mathrm{n}\omega ))+}\to \\ \leftarrow \frac{\mathrm{c}\mathrm{o}{\mathrm{s}}^{\mathrm{2}}\omega )\mathrm{a}\mathrm{r}\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{h}(\mathrm{t}\mathrm{a}\mathrm{n}\omega )}{K_{\mathrm{p}}\mathrm{a}\mathrm{r}\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{h}(\mathrm{t}\mathrm{a}\mathrm{n}\omega )}\end{array}（28）

根据图4可得

{\sigma }_{\mathrm{n}}={\sigma }_{\mathrm{h}A}\mathrm{c}\mathrm{o}{\mathrm{s}}^{\mathrm{2}}\theta +{\sigma }_{\mathrm{v}A}\mathrm{s}\mathrm{i}{\mathrm{n}}^{\mathrm{2}}\theta +{\tau }_A\mathrm{s}\mathrm{i}\mathrm{n}(\mathrm{2}\theta )（29）

式（29）可写为

{\sigma }_{\mathrm{n}}=\frac{K_{\mathrm{1}}\mathrm{c}\mathrm{o}{\mathrm{s}}^{\mathrm{2}}\theta +K_{\mathrm{2}}\mathrm{s}\mathrm{i}{\mathrm{n}}^{\mathrm{2}}\theta }{\mathrm{1}-\mathrm{t}\mathrm{a}\mathrm{n}\phi \mathrm{s}\mathrm{i}\mathrm{n}(\mathrm{2}\theta )}{\overline{\sigma }}_{\mathrm{v}}（30）

图5为土体微分单元受力分析示意图。取一水平微分单元dz为研究对象，在竖直方向上，根据受力平衡可得

\begin{array}{l}\gamma \left(\mathrm{2}{B}_z+\mathrm{d}z\mathrm{t}\mathrm{a}\mathrm{n}\theta \right)\mathrm{d}z+{\overline{\sigma }}_{\mathrm{v}}\left(\mathrm{2}{B}_z+\mathrm{2}\mathrm{d}z\mathrm{t}\mathrm{a}\mathrm{n}\theta \right)=\\ \left({\overline{\sigma }}_{\mathrm{v}}+\mathrm{d}{\overline{\sigma }}_{\mathrm{v}}\right)\mathrm{2}{B}_z+\mathrm{2}{\tau }_{\mathrm{s}}\mathrm{d}z+\mathrm{2}{\sigma }_{\mathrm{n}}\mathrm{t}\mathrm{a}\mathrm{n}\theta \mathrm{d}z\end{array}（31）

假设地表存在超载q，则把边界条件z=\mathrm{0}时{\overline{\sigma }}_{\mathrm{v}}=q代入式（31）得

{\overline{\sigma }}_{\mathrm{v}}=\frac{\gamma B_z}{A_{\mathrm{1}}-\mathrm{t}\mathrm{a}\mathrm{n}\theta }+C{B}_z^{{}^{\frac{A_{\mathrm{1}}}{\mathrm{t}\mathrm{a}\mathrm{n}\theta }}}（32）

其中，

A_{\mathrm{1}}=\frac{(K_{\mathrm{1}}\mathrm{c}\mathrm{o}{\mathrm{s}}^{\mathrm{2}}\theta +K_{\mathrm{2}}\mathrm{s}\mathrm{i}{\mathrm{n}}^{\mathrm{2}}\theta )(\mathrm{t}\mathrm{a}\mathrm{n}\phi +\mathrm{t}\mathrm{a}\mathrm{n}\theta )}{\mathrm{1}-\mathrm{t}\mathrm{a}\mathrm{n}\phi \mathrm{s}\mathrm{i}\mathrm{n}\theta }-\mathrm{t}\mathrm{a}\mathrm{n}\theta

C={\left[\mathrm{0.5}D+\left(H+\mathrm{0.5}D\right)\mathrm{t}\mathrm{a}\mathrm{n}\theta \right]}^{-\frac{A_{\mathrm{1}}}{\mathrm{t}\mathrm{a}\mathrm{n}\theta }}\left\{q+\frac{\gamma \left[\mathrm{0.5}D+\left(H+\mathrm{0.5}D\right)\mathrm{t}\mathrm{a}\mathrm{n}\theta \right]}{\mathrm{t}\mathrm{a}\mathrm{n}\theta -A_{\mathrm{1}}}\right\}

假设隧道的埋深为H时，地表没有外荷载，则式（32）可写为

{\overline{\sigma }}_{\mathrm{v}}=\frac{\gamma {[\mathrm{0.5}D+\left(H+\mathrm{0.5}D\right)\mathrm{t}\mathrm{a}\mathrm{n}\theta ]}^{\mathrm{1}-\frac{A_{\mathrm{1}}}{\mathrm{t}\mathrm{a}\mathrm{n}\theta }}}{\mathrm{t}\mathrm{a}\mathrm{n}\theta -A_{\mathrm{1}}}{B}_z^{{}^{\frac{A_{\mathrm{1}}}{\mathrm{t}\mathrm{a}\mathrm{n}\theta }}}+\frac{\gamma B_z}{A_{\mathrm{1}}-\mathrm{t}\mathrm{a}\mathrm{n}\theta }（33）

联立（23）、式（27）和式（32）得

{\sigma }_{\mathrm{v}x}=m{\overline{\sigma }}_{\mathrm{v}}（34）

通过式（34）即可得出任一地层深度处竖向应力在水平方向上的分布。

联立式（22）、式（26）和式（34）即可得出作用在隧道衬砌结构上侧向土压力为

{\sigma }_{\mathrm{h}}=K_{\mathrm{1}}m{\overline{\sigma }}_{\mathrm{v}}（35）

通过对上式进行分析可知，侧向土压力值自隧道拱顶向下逐渐增大，其分布模式为上凸的曲线（见图2）。

另一方面，对于松动土压力计算过程中关于砂土地层隧道临界埋深的确定，可以根据Terzaghi理论的隧道深浅埋分界深度计算式确定，为

H_L=\frac{\mathrm{1.609}{a}_{\mathrm{1}}}{K_{\mathrm{0}}\mathrm{t}\mathrm{a}\mathrm{n}\phi }（36）

式中：H_L为深浅埋隧道的临界埋置深度；a_{\mathrm{1}}为施工过程中隧道上方土体形成的拱跨度的一半。

3 结果分析

3.1　地层破坏模式

由式（7）可知，当不考虑隧道几何参数影响时，可得出滑裂角β=45°+φ/2，即为Terzaghi提出的破坏模式下的滑裂角大小。为了更详细的描述松动区滑裂角与地层埋深比的关系，本文取埋深比为0.5~5.0，绘制不同内摩擦角（K₀=0.5时）和不同静止侧压力系数（φ=30°时）下的地层滑裂角随埋深比变化曲线，如图6所示。

由图6可知：随着隧道埋深比的增加，滑裂角总体上先增加后减小，但是其减小的速率逐渐降低，慢慢趋于平缓；隧道上方的受扰动区域随着埋深比的增加逐渐增大，也就是 L 在逐渐增大，该结果和方焘等学者基于GeoPIV技术从不同埋深情况下隧道施工引起的地层变形规律类似^［17］；相比于静止侧压力系数，土体的内摩擦角对松动区滑裂角的影响更大，即内摩擦角越大、滑裂角也越大，隧道上方地层受扰动的范围也就越小，该结果和Moussaei等^［8］学者的试验结果也类似，从而证明了本文滑裂角理论推导的合理性。

3.2　侧压力系数

图7为侧压力系数随内摩擦角变化曲线。Terzaghi建议将侧压力系数取值为1，但后来的很多学者经过研究得出侧压力系数并非恒值，其值大小受土体内摩擦角的影响较大。Chen等^［3］认为侧向土压力系数随着土体内摩擦角的增大而增大。Handy^［18］则认为侧向土压力系数随土体有效内摩擦角增加而减小。本文考虑主应力轴旋转和倾斜滑裂面角度对侧向土压力系数的影响，计算得到侧向土压力系数随着土体内摩擦角的增大而减小，其大小介于Terzaghi和Handy提出的数值之间。当不考虑松动区土体的应力分布形式时，即不考虑主应力偏转，则隧道上方同一深度的水平应力分布为1条水平直线，则根据侧压力系数公式可得K₁值退化为1，故Terzaghi提出的建议值为本文公式的特殊情况。

图8为不同主应力轨迹线时的侧压力系数随内摩擦角变化曲线。由图8可知：2种应力分布下侧压力系数均随着倾角θ的增大逐渐减小，与之相对应的是侧压力系数随着倾斜滑裂角β的增大而缓慢增大；当φ<10°时，2种轨迹线下计算的侧压力系数约等于1；当φ>10°时，随着内摩擦角φ的增大，侧压力系数逐渐减小，但相比较圆弧轨迹线，悬链线轨迹线下侧压力系数减小的速率更快；当φ=35°时，悬链线形分布对应的侧压力系数降为0.46，圆弧则为0.60。

3.3　松动土压力

采用本文方法计算隧道拱顶上方松动土压力，同时采用Phase2进行数值模拟以验证本文方法可行性。数值模型中，盾构隧道的直径为10 m，为了消除边界效应的影响，隧道距离左右和下边界的距离为3倍直径。下边界为水平和竖直方向约束，左右边界为水平方向约束。计算过程中采用Mohr-Coulomb破坏准则，地应力为自重应力。土体力学参数和衬砌结构力学参数见表1和表2，模型尺寸如图9所示。

图10为不同埋深比条件下隧道上方的应力分布。由图10可知：隧道上方的应力分布近似为一条下凹的曲线，在隧道轴线位置处竖向应力值最小，自轴线向两边竖向应力逐渐增大；2种应力分布（圆弧，悬链线）的解析解基本和数值模拟结果类似，相对悬链线轨迹而言，圆弧形轨迹下的解析解和数值结果更加吻合。从整个计算式的推导过程以及计算结果的对比来看，圆弧形轨迹下的解析解相对简单且和数值结果更加吻合，因此在计算松动土压力时推荐采用圆弧形大主应力轨迹的假设。

为进一步验证所推导的倾斜滑裂面假设下松动土压力计算式的合理性，本文选择大主应力轨迹为圆弧假设下得出的松动土压力计算式，与Iglesia^［24］和Dewoolkar^［25］的试验值以及现有的松动土压力计算式进行对比。对于Iglesia的试验结果，本文在其所开展的试验中选择其中2种试验工况（\gamma =\mathrm{16.64} kN · m^-3，B_z =25.4 mm；\gamma =\mathrm{16.25} kN · m^-3，B_z =50.8 mm）下对应的3种埋深比数据，相同埋深比下取均值；与此同时，将所选择的2种试验条件代入现有理论方法和本文计算式，计算相应的松动土压力并取均值，然后将试验结果和理论计算结果进行对比。

图11为松动土压力随隧道埋深比的变化曲线。由图11可知：在不同埋深情况下，本文方法和试验数据更加吻合，从而证明了本文理论推导的合理性。Terzaghi公式因其假设滑裂面为垂直型且对于侧压力系数的大小为经验性取值，因此导致其对于松动土压力的预测误差随着埋深比的增大而增大。Gong基于曲线滑裂面的破坏模式虽然更加符合实际情况，但是其假设的隧道上方应力分布为梯形分布和实际情况不符，导致其计算结果和试验数据有一定的误差。虽然Shukla也是基于倾斜滑裂面情况得出的松动土压力的解析解，但是没有给出合理的倾斜角取值，因此计算误差也偏大。

4 结论

（1）本文针对浅埋盾构隧道围岩破坏推导滑裂角计算式，可知滑裂角随土体内摩擦角的增大而增大，随着隧道埋深比的增大先减小然后趋于稳定，总体上大于\mathrm{45}°+\phi /\mathrm{2}。与之相对应，隧道上方的松动区体积随埋深比的增大先增大而后趋于稳定。该结果和现有试验规律基本相同，从而也证明文中滑裂角计算式的合理性。

（2）通过将2种应力分布形式（圆弧、悬链线）下的解析解和数值结果进行对比可知，2种应力分布形式假设都能很好地反映隧道上方松动土压力在水平方向上的分布规律，但是为计算的简便以及实际推广应用，建议采用圆弧形大主应力轨迹的假设计算隧道松动土压力。

（3）本文在对隧道松动区的确定过程中充分考虑了隧道尺寸、埋深以及地质条件的影响，并从大主应力轨迹线为圆弧和悬链线的角度考虑土拱效应对松动土压力的影响，使得本文计算式和试验结果更加吻合，可应用于砂土地质条件下盾构隧道上方的松动土压力计算。

参考文献

[1] 方焘，梁连，陈其志.基于修正椭球体理论的隧道松动区及松动土压力研究［J］.岩土工程学报，2023，45（6）：1113-1122.

[2] FANG Tao， LIANG Lian， CHEN Qizhi. Loosening Zone and Earth Pressure around Tunnels Based on Modified Ellipsoid Theory ［J］. Chinese Journal of Geotechnical Engineering， 2023， 45 （6）： 1113-1122. in Chinese

[3] TERZAGHI K. Theoretical Soil Mechanics ［M］. Hoboken， NJ， USA： John Wiley and Sons， Inc.， 1943： 66-76.

[4] CHEN R X， ZHU B， CHEN Y M， et al. Modified Terzaghi Loosening Earth Pressure Based on Theory of Main Stress Axes Rotation ［J］. Rock and Soil Mechanics， 2010， 31 （5）： 1402-1406.

[5] CHEN K H， PENG F L. An Improved Method to Calculate the Vertical Earth Pressure for Deep Shield Tunnel in Shanghai Soil Layers ［J］. Tunnelling and Underground Space Technology， 2018， 75 （5）： 43-66.

[6] LIANG L J， XU C J. Numerical and Theoretical Research on Stress Distribution in the Loosening Zone of the Trapdoor Problem ［J］. International Journal for Numerical and Analytical Methods in Geomechanics， 2019， 43 （7）： 1426-1447.

[7] WU J， LIAO S M， LIU M B. An Analytical Solution for the Arching Effect Induced by Ground Loss of Tunneling in Sand ［J］. Tunnelling and Underground Space Technology， 2019， 83 （1）： 175-186.

[8] LIN Q T， TIAN Y， LU D C， et al. A Prediction Method of Ground Volume Loss Variation with Depth Induced by Tunnel Excavation ［J］. Acta Geotechnica， 2021， 16 （11）： 3689-3707.

[9] MOUSSAEI N， KHOSRAVI M H， HOSSAINI M F. Physical Modeling of Tunnel Induced Displacement in Sandy Grounds ［J］. Tunnelling and Underground Space Technology， 2019， 90 （8）： 19-27.

[10] HUANG J Z， ZHANG Y， OUYANG X W， et al. Lagged Settlement in Sandy Cobble Strata and Earth Pressure on Shield Tunnel ［J］. Mathematical Biosciences and Engineering： MBE， 2019， 16 （6）： 6209-6230.

[11] 宫全美，张润来，周顺华，等.基于颗粒椭球体理论的隧道松动土压力计算方法［J］.岩土工程学报，2017，39（1）：99-105.

[12] GONG Quanmei， ZHANG Runlai， ZHOU Shunhua， et al. Method for Calculating Loosening Earth Pressure around Tunnels Based on Ellipsoid Theory of Particle Flows ［J］. Chinese Journal of Geotechnical Engineering， 2017， 39 （1）： 99-105. in Chinese

[13] 徐长节，梁禄钜，陈其志，等.考虑松动区内应力分布形式的松动土压力研究［J］.岩土力学，2018，39（6）：1927-1934.

[14] XU Changjie， LIANG Luju， CHEN Qizhi， et al. Research on Loosening Earth Pressure Considering the Patterns of Stress Distribution in Loosening Zone ［J］. Rock and Soil Mechanics， 2018， 39 （6）： 1927-1934. in Chinese

[15] 王将，袁大军，金大龙，等.基于非线型滑动面假设的盾构隧道松动土压力计算模型研究［J］.铁道学报，2021，43（6）：165-172.

[16] WANG Jiang， YUAN Dajun， JIN Dalong， et al. Research on Calculation Model for Loosening Earth Pressure of Shield Tunnel Based on Assumption of Non-Liner Sliding Surface ［J］. Journal of the China Railway Society， 2021， 43 （6）： 165-172. in Chinese

[17] 于丽，吕城，段儒禹，等.浅埋黄土隧道围岩压力计算方法［J］.中国铁道科学，2019，40（4）：69-76.

[18] YU Li， Cheng LÜ， DUAN Ruyu， et al. Calculation Method for Surrounding Rock Pressure of Shallow Buried Loess Tunnel ［J］. China Railway Science， 2019， 40 （4）： 69-76. in Chinese

[19] 刘翔，房倩.变坡面浅埋偏压隧道松动围岩压力计算方法［J］.中国铁道科学，2017，38（2）：64-70.

[20] LIU Xiang， FANG Qian. Calculation Method for Loose Surrounding Rock Pressure of Shallow Tunnel under Asymmetrical Pressure on Variable Slope ［J］. China Railway Science， 2017， 38 （2）： 64-70. in Chinese

[21] SHUKLA S， SIVAKUGAN N， GAURAV G. A Simplified Extension of the Conventional Theory of Arching in Soils ［J］. International Journal of Geotechnical Engineering， 2009， 3 （3）： 353-359.

[22] 李瑞林，周国庆，林超，等.考虑土拱效应的滑移面间非极限状态土压力解答［J］.岩土力学，2017，38（11）：3145-3153.

[23] LI Ruilin， ZHOU Guoqing， LIN Chao， et al. Solution of Earth Pressure between Slip Surfaces under Non-Limit State Considering Soil Arching Effect ［J］. Rock and Soil Mechanics， 2017， 38 （11）： 3145-3153. in Chinese

[24] 方焘，梁连，颜建伟.不同埋深下盾构隧道施工引起的地层变形试验［J］.长江科学院院报，2023，40（3）：85-92.

[25] FANG Tao， LIANG Lian， YAN Jianwei. Experimental Study on Stratum Deformation Caused by Shield Tunnelling at Different Buried Depths ［J］. Journal of Yangtze River Scientific Research Institute， 2023， 40 （3）： 85-92. in Chinese

[26] HANDY R L. The Arch in Soil Arching ［J］. Journal of Geotechnical Engineering， 1985， 111 （3）： 302-318.

[27] PAIK K H， SALGADO R. Estimation of Active Earth Pressure against Rigid Retaining Walls Considering Arching Effects ［J］. Géotechnique， 2003， 53 （7）： 643-653.

[28] HARROP-WILLIAMS K. Arch in Soil Arching ［J］. Journal of Geotechnical Engineering， 1989， 115 （3）： 415-419.

[29] XU C J， CHEN Q Z， LUO W J， et al. Analytical Solution for Estimating the Stress State in Backfill Considering Patterns of Stress Distribution ［J］. International Journal of Geomechanics， 2019， 19 （1）： 04018189.

[30] VIVANCO F， RICA S， MELO F. Dynamical Arching in a Two Dimensional Granular Flow ［J］. Granular Matter， 2012， 14 （5）： 563-576.

[31] 汪大海，贺少辉，刘夏冰，等.基于主应力旋转特征的浅埋隧道上覆土压力计算及不完全拱效应分析［J］.岩石力学与工程学报，2019，38（6）：1284-1296.

[32] WANG Dahai， HE Shaohui， LIU Xiabing， et al. A Modified Method for Determining the Overburden Pressure above Shallow Tunnels Considering the Distribution of the Principal Stress Rotation and the Partially Mobilized Arching Effect ［J］. Chinese Journal of Rock Mechanics and Engineering， 2019， 38 （6）： 1284-1296. in Chinese

[33] IGLESIA G R， EINSTEIN H H， WHITMAN R V. Validation of Centrifuge Model Scaling for Soil Systems via Trapdoor Tests ［J］. Journal of Geotechnical and Geoenvironmental Engineering， 2011， 137 （11）： 1075-1089.

[34] DEWOOLKAR M M， SANTICHAIANANT K， KO H Y. Centrifuge Modeling of Granular Soil Response over Active Circular Trapdoors ［J］. Soils and Foundations， 2007， 47 （5）： 931-945.

基金资助

国家自然科学基金资助项目(52168048)