基于多切线技术的盾构隧道开挖面稳定性分析

钟军豪; CHIANSiau Chen; 李永鑫; 刘明芳; 陈明辉; 吴振元; 龙桂华; 杨小礼; 周德

doi:10.11779/CJGE20231160

基于多切线技术的盾构隧道开挖面稳定性分析 English Version

1.
中南大学土木工程学院, 湖南长沙 410075
2.
新加坡国立大学土木与环境工程系, 新加坡 117576
3.
合肥工业大学汽车与交通工程学院, 安徽合肥 230009
4.
深圳市市政工程总公司, 广东深圳 518000

基金项目:

国家自然科学基金项目 5210081105

安徽省自然科学基金项目 2108085QE250

湖南省研究生科研创新项目 CX20220239

国家留学基金项目（国家建设高水平大学公派研究生项目）

详细信息

作者简介:
钟军豪(1996—)，男，博士研究生，主要从事隧道开挖面稳定性和地层建模方面的研究工作。E-mail: zhongjunhao@csu.edu.cn

通讯作者:
周德, E-mail: 210026@csu.edu.cn

中图分类号: TU432
计量
- 文章访问数: 341
- HTML全文浏览量: 47
- PDF下载量: 77
出版历程
- 收稿日期: 2023-11-26
- 网络出版日期: 2024-04-18
- 刊出日期: 2024-12-31

Face stability of shield-driven tunnels using multi-tangent technique

1.
School of Civil Engineering, Central South University, Changsha 410075, China
2.
Department of Civil and Environmental Engineering, National University of Singapore, Singapore 117576
3.
School of Automotive and Transportation Engineering, Hefei University of Technology, Hefei 230009, China
4.
Shenzhen Municipal Engineering Corporation, Shenzhen 518000, China

摘要

摘要: 针对非线性条件下深埋隧道开挖面稳定性问题，采用多切线技术分段线性近似Power-Law（P-L）准则，基于极限分析上限定理，构建1种新的开挖面多曲线圆锥破坏机构，提出1种考虑拉应力截断Tension cut-off（T-C）效应的开挖面稳定性理论计算方法。通过与既有文献结果对比，验证所提方法的可靠性。根据反算的应力分布和破坏机构几何特点，探究不同参数条件下T-C对开挖面稳定性的影响。结果表明：①在线性条件下，黏聚力较大且内摩擦角较小时，T-C对开挖面临界支护力的影响较大；在非线性条件下，m的增加并不会放大T-C对开挖面稳定性的影响。②在抗拉强度较大时，T-C对开挖面稳定性影响较为显著，随着抗拉强度的折减，T-C控制的破坏机构尺寸逐渐减小，且顶部区域明显钝化。
- 开挖面稳定性 /
- 多切线技术 /
- Power-Law准则 /
- 拉应力截断 /
- 极限分析
Abstract: To address the face stability problem of deep-buried tunnels under nonlinear conditions, using Power-Law (P-L) criterion the approximate piecewise and by multi-tangent technique, a novel multi-cone failure mechanism is developed for tunnel faces based on the kinematic limit analysis theorem. A theoretical framework considering the tension cut-off (T-C) is proposed to evaluate the stability of tunnel faces. The reliability of the proposed method is validated through comparison with that in the existing literatures. Based on the back-calculated stress distribution and geometric feature of the failure mechanism, the influences of T-C under different parameters on the stability of tunnel faces are investigated. The results indicate that: (1) under the linear condition, the effects of T-C on the critical support pressure are more significant with high cohesion and small internal friction angle. Under the nonlinear condition, an increase in the nonlinear coefficient will not exacerbate the impact of T-C on the face stability. (2) For higher tensile strength, the effects of T-C on the stability of tunnel faces are significant. As the tensile strength decreases, the failure mechanism governed by T-C gradually shrinks, and the top of the failure mechanism exhibits noticeable blunting and looks like a dome.
- stability of tunnel face /
- multi-tangent technique /
- Power-Law criterion /
- tension cut-off /
- limit analysis

HTML全文

0. 引言

近年来，全球区域性暴雨频发，降雨量远高于以往水平^[1-2]，对铁路线路的安全运营带来挑战。1965年，日本新干线就曾因连续降雨等原因导致线路沉降过大，不得不将列车运行速度由最初设计的210 km/h下降至110~180 km/h^[3]。1997年3月份，哥伦比亚Conrad地区的铁路线路，由于强降雨导致路基和地基土体含水率迅速增大，在列车仅以平均车速43.45 km/h通过时，线路突然坍塌，造成重大安全事故^[4]。对于土体，水的存在既减小了土颗粒之间的摩擦力，影响线路的运营。由此可见，实际铁路线路的稳定性受移动荷载与含水路基/地基形成的水-动力耦合作用的影响。由于自然界中土体的含水率处于变化中，其饱和度也随之改变；因此，将路基/地基视为土颗粒、水、气组成的三相非饱和土，研究不同饱和度下，移动荷载作用时的动力响应可以对运行线路的稳定性进行更准确的预判。

移动荷载作用下的地基动力响应问题最初是将地基视为单相弹性介质^[5-7]。近年来，国内外许多学者基于Biot饱和多孔介质波传播理论^[8]，将单相介质扩展到两相介质，建立了将地基视为饱和两相介质的动力求解方法和模型。Lu等^[9]建立了三维解析解模型，考虑了饱和多孔介质及移动点荷载，研究了动力响应与土体剪切波波速的关系。Cai等^[10]采用双重傅立叶变换和逆变换，考虑了流固耦合作用，研究了荷载加载方式和饱和地基土参数对动力响应的影响。Gao等^[11]、Bian等^[12]和胡静等^[13]基于2.5维有限元法研究了饱和土渗透系数和列车速度对地基动力响应的影响。考虑到实际土体的三相性，徐明江^[14]基于连续介质力学理论，结合土-水特征曲线（soil-water characteristic curve，简称SWCC）以及Mualem^[15]理论，推导了移动荷载作用下的非饱和土波动方程。Lu等^[16]、Fang等^[17]采用积分变换的方式求解了该方程。同样基于该方程，Gao等^[18]、李绍毅^[19]研究了列车荷载作用下的非饱和地基振动。

在非饱和土动力响应的研究中^{[14, 16-19]}，需要基于SWCC推导不同饱和度下的孔隙流体渗透性和土的强度，因而SWCC对土体的动力响应具有决定性作用。SWCC描述的是非饱和土土体含水率和基质吸力之间的关系，早期的SWCC研究中，通常假设基质吸力和含水率之间存在着唯一关系^[20-22]。实际上，除了含水率，变形也会使土体饱和度发生改变。非饱和土在荷载作用下，土体必然会发生相应的变形，因而，将变形效应考虑进SWCC模型中更为合理。遗憾的是，已有关于非饱和土动力响应的研究均采用早期的V-G模型^[21]来描述土-水特征曲线，未考虑变形对SWCC的影响，不能完整的描述荷载作用下的非饱和土持水特性。

综上所述，目前还未出现基于考虑变形效应的SWCC的非饱和土动力响应求解方法。为此，本文在现有研究的基础上，建立了考虑变形效应的SWCC模型，并据此推导了考虑SWCC变形效应的非饱和土波动方程，利用2.5有限元法对该方程进行求解；求解结果与现有文献中的解析解进行了验证；最后对该方法的计算效率进行了对比，并对引入变形效应的影响进行了初步分析。

1. 理论求解

1.1 考虑变形的土-水特征曲线

孙德安^[23]的研究表明，非饱和土的土-水特征曲线与应力历史和应力状态无直接关系，在土体受力变形过程中，可以用孔隙比来表示孔隙的变化。蔡国庆等^[24]、胡冉等^[25]及张雪东等^[26]同样也是采用孔隙比来表征变形。图 1整理了不同初始孔隙比的土在等吸力下的等向压缩试验结果；其中，孙德安^[23]对珍珠黏土进行了基质吸力为147 kPa的等向压缩试验；蔡国庆等^[24]给出了膨润土/高岭土混合土体在基质吸力为100 kPa的等向压缩试验结果。图 1中s为基质吸力，e为孔隙比，S_r为孔隙水饱和度，简称饱和度。

图 1 饱和度与孔隙比的关系图

Figure 1. Relationship between saturation and void ratio

下载: 全尺寸图片幻灯片

从中可以发现，吸力一定时，土体受外力导致的孔隙比减小会使饱和度增加。不同初始孔隙比条件下，等吸力时的孔隙比与饱和度近似呈线性关系且直线斜率相近。可见，考虑变形效应的SWCC可近似认为是将原有SWCC进行偏移得到。假定中的直线斜率为 ${\lambda _{{\text{se}}}}$ ，结合蔡国庆等^[24]关于孔隙比与SWCC的研究，不同孔隙比下的SWCC可总结为图 2。

图 2 不同孔隙比下的土-水特征曲线示意图

Figure 2. Schematic diagram of soil-water characteristic curve under different void ratios

下载: 全尺寸图片幻灯片

根据图 2，考虑变形效应的SWCC模型可以表示为

，

${S_{\text{r}}} = {f_{{S_{\text{r}}}}}(e) + {f_{{S_{\text{r}}}}}(s) \text{，}$

(1)

。

${f_{{S_{\text{r}}}}}(e) = {\lambda _{{\text{se}}}}(e - {e_0}) 。$

(2)

式中： ${e_0}$ 为初始孔隙比； ${f_{{S_{\text{r}}}}}(e)$ 表示饱和度与孔隙比之间的关系； ${f_{{S_{\text{r}}}}}(s)$ 为 ${e_0}$ 时的SWCC。

在众多SWCC模型中，V-G模型适用范围广，并且用于求解非饱和动力问题时，能够形成关于时间的一阶常微分方程，实现求解，因此现有关于非饱和土的动力求解均是基于V-G模型^{[14, 16-19]}。当 ${f_{{S_{\text{r}}}}}(s)$ 为V-G模型，即 $\frac{{S_{\text{r}}^0 - {S_{{\text{w0}}}}}}{{(1 - {S_{{\text{w0}}}})}} = {\text{ }}{\left\{ {1 + {{(\beta s)}^k}} \right\}^{ - m}}$ ， $S_{\text{r}}^0$ 为 ${e_0}$ 时V-G模型中的饱和度， ${S_{{\text{w0}}}}$ 为残余饱和度，s为基质吸力， $\beta$ ，m，k为V-G模型的参数，其中m=1-1/k。则，考虑变形效应的SWCC模型可以写为

${S_{\text{r}}} = {S_{{\text{w0}}}} + {\lambda _{{\text{se}}}}(e - {e_0}) + (1 - {S_{{\text{w0}}}}){\left\{ {1 + {{(\beta s)}^k}} \right\}^{ - m}} 。$

(3)

此时，有效饱和度 ${S_{\text{e}}}$ 为

${S_{\text{e}}} = \frac{{[S_{\text{r}}^{\text{0}} + {\lambda _{{\text{se}}}}(e - {e_0})] - {S_{{\text{w0}}}}}}{{(1 - {S_{{\text{w0}}}})}} 。$

(4)

根据Mualem理论^[15]，非饱和土中流体的渗透系数是关于有效饱和度的函数为

${k_{\text{w}}} = \frac{{{\rho _{\text{w}}}g\kappa }}{{{\eta _{\text{w}}}}}\sqrt {{S_{\text{e}}}} {\left\{ {1 - {{\left[ {1 - {{({S_\text{e}})}^{\frac{1}{m}}}} \right]}^m}} \right\}^2} \text{，}$

(5)

${k_{\text{a}}} = \frac{{{\rho _{\text{a}}}g\kappa }}{{{\eta _{\text{a}}}}}\sqrt {1 - {S_{\text{e}}}} {\left[ {1 - {{({S_{\text{e}}})}^{\frac{1}{m}}}} \right]^{2m}} 。$

(6)

式中： ${k_{\text{w}}}$ ， ${k_{\text{a}}}$ 分别为孔隙水和孔隙气的达西渗透系数； ${\rho _{\text{w}}}$ ， ${\rho _{\text{a}}}$ 为孔隙水和孔隙气的密度； ${\eta _{\text{w}}}$ ， ${\eta _{\text{a}}}$ 分别为孔隙水和孔隙气的黏滞系数（Pa·s）；g为重力加速度； $\kappa$ 为固有渗透系数（m²）。

当考虑变形的影响时， ${S_{\text{e}}}$ 的表达式发生变化，流体渗透系数也发生相应的改变。

与此同时，徐明江^[14]根据试验结果和经验公式，建立了非饱和土的动剪切模量与有效饱和度之间的联系：

$G = {G_{\text{s}}} + 2050\tan \varphi '\int\limits_0^s {{S_{\text{e}}}} {\text{d}}s 。$

(7)

式中： $G$ 为非饱和土的动剪切模量； ${G_{\text{s}}}$ 为土体饱和状态下的动剪切模量， $\varphi '$ 为土体的有效内摩擦角。

1.2 非饱和多孔介质动力控制方程

固、液、气三相的质量守恒方程分别为

$\frac{{(1 - n)\partial {\rho _{\text{s}}}}}{{\partial t}} - {\rho _{\text{s}}}\frac{{\partial n}}{{\partial t}} + {\rho _{\text{s}}}(1 - n)\nabla \dot u = 0 \text{，}$

(8)

$\frac{{{S_{\text{r}}}{\rho _{\text{w}}}\partial n}}{{\partial t}} + \frac{{{\rho _{\text{w}}}n\partial {S_{\text{r}}}}}{{\partial t}} + \frac{{n{S_{\text{r}}}\partial {\rho _{\text{w}}}}}{{\partial t}} + {\rho _{\text{w}}}n{S_{\text{r}}}\nabla {\dot u^{\text{w}}} = 0\text{，}$

(9)

$\frac{{{S_{\text{a}}}{\rho _{\text{a}}}\partial n}}{{\partial t}} + \frac{{{\rho _{\text{a}}}n\partial {S_{\text{a}}}}}{{\partial t}} + \frac{{n{S_{\text{a}}}\partial {\rho _{\text{a}}}}}{{\partial t}} + {\rho _{\text{a}}}n{S_{\text{a}}}\nabla {\dot u^{\text{a}}} = 0 。$

(10)

式中： $n$ 为孔隙率；t为时间； ${\rho _{\text{s}}}$ 为土骨粒密度； $u$ ， ${u^{\text{w}}}$ ， ${u^{\text{a}}}$ 分别为土骨架位移，孔隙水位移，孔隙气位移，上标 $\cdot$ 表示对时间的导数； $\nabla$ 为汉密尔顿算子； ${S_{\text{a}}}$ 为孔隙气饱和度， ${S_{\text{a}}} = 1 - {S_{\text{r}}}$ 。

固、液、气三相应力引起土单元压缩变形的本构方程可以表示为

$\frac{{\partial {\rho _{\text{s}}}}}{{{\rho _{\text{s}}}\partial t}} = - \frac{{\partial \sigma _{ij}^{\text{s}}}}{{3{K_{\text{s}}}{\text{d}}t}} \text{，}$

(11)

$\frac{{\partial {\rho _{\text{w}}}}}{{{\rho _{\text{w}}}\partial t}} = \frac{{\partial {P^{\text{w}}}}}{{{K_{\text{w}}}{\text{d}}t}} \text{，}$

(12)

$\frac{{\partial {\rho _{\text{a}}}}}{{{\rho _{\text{a}}}\partial t}} = \frac{{\partial {P^{\text{a}}}}}{{{K_{\text{a}}}{\text{d}}t}} 。$

(13)

式中： $\sigma _{ij}^{\text{s}}$ 为作用在土颗粒上的应力； ${K_{\text{s}}}$ ， ${K_{\text{w}}}$ ， ${K_{\text{a}}}$ 分别为土颗粒、孔隙水、孔隙气的压缩模量； ${P^{\text{w}}}$ 为孔隙水压， ${P^{\text{a}}}$ 为孔隙气压。

考虑土颗粒由于孔隙流体压力作用产生的变形时，由弹性理论，有效应力分量可表示为

${\sigma '_{ij}} = \lambda \Theta {\delta _{ij}} + 2\mu {\varepsilon _{ij}} + \frac{{{K_{\text{b}}}}}{{{K_{\text{s}}}}}P{\delta _{ij}} 。$

(14)

式中： $\lambda$ ， $\mu$ 为Lame常数； $\Theta$ 为土骨架的体积应变； ${\varepsilon _{ij}}$ 为土颗粒的应变； ${\delta _{ij}}$ 为克罗内克符号； ${K_{\text{b}}}$ 为土骨架的压缩模量且 ${K_{\text{b}}} = \lambda + \frac{{2\mu }}{3}$ ；P为孔隙流体压力，包括孔隙水压及孔隙气压， $P = {S_{\text{r}}}{P^{\text{w}}} + {S_{\text{a}}}{P^{\text{a}}}$ 。

根据应力空间平均化方法，土单元的总应力可以表示为

${\sigma _{ij}} = (1 - n)\sigma _{ij}^{\text{s}} - n{S_{\text{r}}}{P^{\text{w}}}{\delta _{ij}} - n{S_{\text{a}}}{P^{\text{a}}}{\delta _{ij}} 。$

(15)

根据Bishop单一变量非饱和土有效应力原理：

${\sigma '_{ij}} = {\sigma _{ij}} + {\delta _{ij}}P \text{，}$

(16)

将式（16）代入式（14），可以得到

${\sigma _{ij}} + {\delta _{ij}}P = \lambda \Theta {\delta _{ij}} + 2\mu {\varepsilon _{ij}} + P{\delta _{ij}} 。$

(17)

将式（15）代入式（17），移项后可以得到

$(1 - n)\sigma _{ij}^{\text{s}} = \lambda \Theta {\delta _{ij}} + 2\mu {\varepsilon _{ij}} - (\alpha - n)P{\delta _{ij}} 。$

(18)

式中， $\alpha = 1 - \frac{{{K_{\text{b}}}}}{{{K_{\text{s}}}}}$ 。

将式（11）代入式（18）并展开，可以得到

$\begin{aligned} &-\frac{\partial \sigma_{i j}^{\mathrm{s}}}{3 K_{\mathrm{s}} \mathrm{~d} t}=\frac{1}{3 K_{\mathrm{s}}} \frac{1}{n-1}\left[2 \mu \frac{\partial \varepsilon_{i j}}{\partial t}+\lambda \frac{\partial \Theta}{\partial t} \delta_{i j}-\right. \\ &\left.(\alpha-n) S_{\mathrm{r}} \frac{\partial P^{\mathrm{w}}}{\partial t}-(\alpha-n) S_{\mathrm{a}} \frac{\partial P^{\mathrm{a}}}{\partial t}\right] \end{aligned}$

(19)

将式（19）代入式（8）中，最终可以求解出孔隙率的偏导数：

$\frac{{\partial n}}{{\partial t}} = (\alpha - n)\nabla \dot u + \frac{{(\alpha - n){S_{\text{r}}}}}{{{K_{\text{s}}}}}\frac{{\partial {P^{\text{w}}}}}{{\partial t}} + \frac{{(\alpha - n){S_{\text{a}}}}}{{{K_{\text{s}}}}}\frac{{\partial {P^{\text{a}}}}}{{\partial t}}。$

(20)

结合式（9），（10），（20），可得

$\begin{aligned} & (\alpha-n) \nabla \dot{u}+n S_{\mathrm{r}} \nabla \dot{u}^{\mathrm{w}}+n S_{\mathrm{a}} \nabla \dot{u}^{\mathrm{a}}+ \\ & {\left[\frac{(\alpha-n) S_{\mathrm{r}}}{K_{\mathrm{s}}}+\frac{n S_{\mathrm{r}}}{K_{\mathrm{w}}}\right] \frac{\partial P^{\mathrm{w}}}{\partial t}+\left[\frac{(\alpha-n) S_{\mathrm{a}}}{K_{\mathrm{s}}}+\frac{n S_{\mathrm{a}}}{K_{\mathrm{a}}}\right] \frac{\partial P^{\mathrm{a}}}{\partial t}=0 。} \end{aligned}$

(21)

从式（3）推导出饱和度对时间的导数：

$\begin{aligned} \frac{\partial S_{\mathrm{r}}}{\partial t}= & \lambda_{\mathrm{se}} \frac{\partial e}{\partial t}-\beta k m\left(1-S_{\mathrm{w} 0}\right) \\ & \left(\frac{S_{\mathrm{r}}^0-S_{\mathrm{w} 0}}{1-S_{\mathrm{w} 0}}\right)^{\frac{m+1}{m}}\left[\left(\frac{S_{\mathrm{r}}^0-S_{\mathrm{w} 0}}{1-S_{\mathrm{w} 0}}\right)^{-\frac{1}{m}}-1\right]^{\frac{k-1}{k}} \frac{\partial s}{\partial t}。 \end{aligned}$

(22)

根据定义，孔隙比 $e = \frac{n}{{1 - n}}$ ，基质吸力 $s = {P^{\text{a}}} - {P^{\text{w}}}$ ，则式（22）可以改写为

$\begin{aligned} \frac{\partial S_{\mathrm{r}}}{\partial t}= & \lambda_{\mathrm{se}} \frac{1}{(1-n)^2} \frac{\partial n}{\partial t}-\beta k m\left(1-S_{\mathrm{w} 0}\right)\left(\frac{S_{\mathrm{r}}^0-S_{\mathrm{w} 0}}{1-S_{\mathrm{w} 0}}\right)^{\frac{m+1}{m}}\cdot \\ & {\left[\left(\frac{S_{\mathrm{r}}^0-S_{\mathrm{w} 0}}{1-S_{\mathrm{w} 0}}\right)^{-\frac{1}{m}}-1\right]^{\frac{k-1}{k}}\left(\frac{\partial P^{\mathrm{a}}}{\partial t}-\frac{\partial P^{\mathrm{w}}}{\partial t}\right) 。} \end{aligned}$

(23)

令 ${A}_{\text{e}}={\lambda }_{\text{se}}\frac{1}{{(1-n)}^{2}}，{A}_{\text{s}}=-\beta km(1-{S}_{\text{w0}}){\left(\frac{{S}_{\text{r}}^{\text{0}}-{S}_{\text{w0}}}{1-{S}_{\text{w0}}}\right)}^{\frac{m+1}{m}}\cdot$ ${\left[ {{{\left( {\frac{{S_{\text{r}}^{\text{0}} - {S_{{\text{w0}}}}}}{{1 - {S_{{\text{w0}}}}}}} \right)}^{ - \frac{1}{m}}} - 1} \right]^{\frac{{k - 1}}{k}}}$ ；结合式（9），（10），（20）和（23），可得

$\begin{aligned} & A_{\mathrm{e}}(\alpha-n) \nabla \dot{u}+S_{\mathrm{r}} S_{\mathrm{a}} \nabla \dot{u}^{\mathrm{w}}-S_{\mathrm{r}} S_{\mathrm{a}} \nabla \dot{u}^{\mathrm{a}}+ \\ & {\left[\frac{A_{\mathrm{e}}(\alpha-n) S_{\mathrm{r}}}{K_{\mathrm{s}}}-A_{\mathrm{s}}+\frac{S_{\mathrm{r}} S_{\mathrm{a}}}{K_{\mathrm{w}}}\right] \frac{\partial P^{\mathrm{w}}}{\partial t}+} \\ & \quad\left[\frac{A_{\mathrm{e}}(\alpha-n) S_{\mathrm{a}}}{K_{\mathrm{s}}}+A_{\mathrm{s}}-\frac{S_{\mathrm{r}} S_{\mathrm{a}}}{K_{\mathrm{a}}}\right] \frac{\partial P^{\mathrm{a}}}{\partial t}=0。 \end{aligned}$

(24)

式（21），（24）组成求解 $\frac{\partial {P}^{\text{w}}}{\partial t}和\frac{\partial {P}^{\text{a}}}{\partial t}$ 的方程组：

$\left.\begin{array}{l} B_{11} \frac{\partial P^{\mathrm{w}}}{\partial t}+B_{12} \frac{\partial P^{\mathrm{a}}}{\partial t}+B_{13} \nabla \dot{u}+B_{14} \nabla \dot{u}^{\mathrm{w}}+B_{15} \nabla \dot{u}^{\mathrm{a}}=0, \\ B_{21} \frac{\partial P^{\mathrm{w}}}{\partial t}+B_{22} \frac{\partial P^{\mathrm{a}}}{\partial t}+B_{23} \nabla \dot{u}+B_{24} \nabla \dot{u}^{\mathrm{w}}+B_{25} \nabla \dot{u}^{\mathrm{a}}=0。 \end{array}\right\}$

(25)

式中： ${B_{11}} = \frac{{(\alpha - n){S_{\text{r}}}}}{{{K_{\text{s}}}}} + \frac{{n{S_{\text{r}}}}}{{{K_{\text{w}}}}}$ ， ${B_{12}} = \frac{{(\alpha - n){S_{\text{a}}}}}{{{K_{\text{s}}}}} + \frac{{n{S_{\text{a}}}}}{{{K_{\text{a}}}}}$ ， ${B_{13}} = \alpha - n$ ， ${B_{14}} = n{S_{\text{r}}}$ ， ${B_{15}} = n{S_{\text{a}}}$ ； ${B_{21}} =$ $\frac{{A}_{\text{e}}(\alpha -n){S}_{\text{r}}}{{K}_{\text{s}}}-{A}_{\text{s}}+\frac{{S}_{\text{r}}{S}_{\text{a}}}{{K}_{\text{w}}}，{B}_{22}=\frac{{A}_{\text{e}}(\alpha -n){S}_{\text{a}}}{{K}_{\text{s}}}+{A}_{\text{s}}-\frac{{S}_{\text{r}}{S}_{\text{a}}}{{K}_{\text{a}}}，$ ${B_{23}} = {A_{\text{e}}}(\alpha - n)$ ， ${B_{24}} = - {B_{25}} = {S_{\text{r}}}{S_{\text{a}}}$ 。当不考虑变形对SWCC的影响，即A_e为0时，方程（25）退化后的结果与Lu等^[16]和徐明江^[14]的结果一致。

令 $M=\frac{\alpha -n}{{K}_{\text{s}}}，{M}_{\text{w}}=\frac{n}{{K}_{\text{w}}}，{M}_{\text{a}}=\frac{n}{{K}_{\text{a}}}，{C}_{1}=\frac{{A}_{\text{e}}n+{S}_{\text{r}}}{n}，$ ${C_2} = \frac{{{A_{\text{e}}}n - {S_{\text{a}}}}}{n}$ ；采用相对位移， ${w}_{i}=n{S}_{\text{r}}({u}_{i}^{\text{w}}-{u}_{i}) ，$ ${v_i} = n{S_{\text{a}}}(u_i^{\text{a}} - {u_i})$ ，可求得

$- \frac{{\partial {P^{\text{w}}}}}{{\partial t}} = {D_{11}}\nabla \dot u + {D_{12}}\nabla \dot w + {D_{13}}\nabla \dot v \text{，}$

(26)

$- \frac{{\partial {P^a}}}{{\partial t}} = {D_{21}}\nabla \dot u + {D_{22}}\nabla \dot w + {D_{23}}\nabla \dot v 。$

(27)

式中： ${w_i}$ 为液相与土骨架的相对位移； ${v_i}$ 为气相与土骨架的相对位移，D₁₁，D₁₂，D₁₃，D₂₁，D₂₂，D₂₃为方程系数，具体表达式见附录1。

根据动量守恒原理，土体的运动方程可以表示为

${\sigma _{ij,j}} = {\rho _{\text{s}}}{\ddot u_i} + {\rho _{\text{w}}}{\ddot u_i}^{\text{w}} + {\rho _{\text{a}}}{\ddot u_i}^{\text{a}} 。$

(28)

式中：上标 $\cdot \cdot$ 表示对时间的二阶导数。

根据广义达西定律，孔隙水和孔隙气的渗流运动方程可以表示为

$n{S_{\text{r}}}(\dot u_i^{\text{w}} - {\dot u_i}) = \frac{{{k_{\text{w}}}}}{{{\rho _{\text{w}}}g}}( - P_{,i}^{\text{w}} - {\rho _{\text{w}}}\ddot u_i^{\text{w}}) \text{，}$

(29)

$n{S_{\text{a}}}(\dot u_i^{\text{a}} - {\dot u_i}) = \frac{{{k_a}}}{{{\rho _{\text{a}}}g}}( - P_{,i}^{\text{a}} - {\rho _{\text{a}}}\ddot u_i^{\text{a}}) 。$

(30)

将式（17）代入式（28）~（30），采用式（26）中的相对位移，可以得到三相介质运动方程：

$(\lambda + \mu )\nabla (\nabla u) + \mu {\nabla ^2}u - {\delta _{ij}}\alpha P = \rho \ddot u + {\rho _{\text{w}}}\ddot w + {\rho _{\text{a}}}\ddot v \text{，}$

(31)

$- P,_i^{\text{w}} = {\rho _{\text{w}}}{\ddot u_i} + \frac{{{\rho _{\text{w}}}}}{{n{S_{\text{r}}}}}{\ddot w_i} + \frac{{{\rho _{\text{w}}}g}}{{{k_{\text{w}}}}}{\dot w_i} \text{，}$

(32)

$- P,_i^{\text{a}} = {\rho _{\text{a}}}{\ddot u_i} + \frac{{{\rho _{\text{a}}}}}{{n{S_{\text{a}}}}}{\ddot v_i} + \frac{{{\rho _{\text{a}}}g}}{{{k_{\text{a}}}}}{\dot v_i} 。$

(33)

式中： $\rho$ 为非饱和土的总密度， $\rho = (1 - n){\rho _{\text{s}}} + n{S_r}{\rho _{\text{w}}} +$ $n{S_{\text{a}}}{\rho _{\text{a}}}$ 。

将式（26），（27）代入式（31）~（33），得到最终考虑变形效应的非饱和土动力控制方程为

$\begin{aligned} & \mu u_{i, j j}+\left(\lambda+\mu+\alpha S_{\mathrm{r}} D_{11}+\alpha S_{\mathrm{a}} D_{21}\right) u_{j, j i}+ \\ & \left(\alpha S_{\mathrm{r}} D_{12}+\alpha S_a D_{22}\right) w_{j, j i}+ \\ & \left(\alpha S_{\mathrm{r}} D_{13}+\alpha S_{\mathrm{a}} D_{23}\right) v_{j, j i}=\rho \ddot{u}+\rho_{\mathrm{w}} \ddot{w}+\rho_{\mathrm{a}} \ddot{v}, \end{aligned}$

(34)

${D_{11}}{u_{j,ji}} + {D_{12}}{w_{j,ji}} + {D_{13}}{v_{j,ji}} = {\rho _{\text{w}}}{\ddot u_i} + \frac{{{\rho _{\text{w}}}}}{{n{S_{\text{r}}}}}{\ddot w_i} + \frac{{{\rho _{\text{w}}}g}}{{{k_{\text{w}}}}}{\dot w_i} \text{，}$

(35)

$\begin{array}{*{20}{c}} {{D_{21}}{u_{j,ji}} + {D_{22}}{w_{j,ji}} + {D_{23}}{v_{j,ji}} = {\rho _{\text{a}}}{{\ddot u}_i} + \frac{{{\rho _a}}}{{n{S_{\text{a}}}}}{{\ddot v}_i} + \frac{{{\rho _{\text{a}}}g}}{{{k_{\text{a}}}}}{{\dot v}_i}} \end{array}。$

(36)

1.3 2.5维有限元求解

2.5维有限元法是由Yang等^[27]用于求解移动荷载作用下完全弹性半空间的响应。该方法的特点在于需要假定研究对象在沿荷载移动方向的几何形状以及材料属性都是连续一致的，通过傅里叶变换将荷载移动方向的空间坐标（假定为x方向）变换成波数，只在频域-波数域内求解出y-z二维平面的响应，再通过傅里叶逆变换就可以得到x方向的响应，从而获得三维空间响应的解答。由于轨道结构大都具备沿列车行进方向连续一致的特性，因此，2.5维有限元法近年来在交通岩土领域得到了广泛的应用^[11-13]。

在应用时，首先定义关于荷载移动方向x以及时间t的傅里叶变换：

$\tilde{\bar{u}}\left(\xi_x, y, z, \omega\right)=\int\limits_{-\infty}^{+\infty} \int\limits_{-\infty}^{+\infty} u(x, y, z, t) \mathrm{e}^{\mathrm{i} \xi_x x} \mathrm{e}^{-\mathrm{i} \omega t} \mathrm{~d} x \mathrm{~d} t。$

(37)

相应的傅里叶逆变换为

$u(x, y, z, t)=\frac{1}{4 \pi^2} \int\limits_{-\infty}^{+\infty} \int\limits_{-\infty}^{+\infty} \tilde{\bar{u}}\left(\xi_x, y, z, \omega\right) \mathrm{e}^{\mathrm{i} \xi_{x^x}} \mathrm{e}^{\mathrm{i} \omega t} \mathrm{~d} \xi_x \mathrm{~d} \omega。$

(38)

式中： ${\xi _x}$ 为沿x方向的波数（m^-1）； $\omega$ 为圆角频率（rad/s）；上标‘–’和‘~’分别为波数域和频域内的量。

对式（34）~（36）进行傅里叶变换后，非饱和土的动力控制方程在频域内的表达式为

$\begin{aligned} & \mu u_{i, j j}+\left(\lambda+\mu+\alpha S_{\mathrm{r}} D_{11}+\alpha S_{\mathrm{a}} D_{21}\right) u_{j, j i}+\rho \omega^2 u_i+ \\ & \left(\alpha S_{\mathrm{r}} D_{12}+\alpha S_{\mathrm{a}} D_{22}\right) w_{j, j i}+\rho_{\mathrm{w}} \omega^2 w_i+ \\ & \quad\left(\alpha S_{\mathrm{r}} D_{13}+\alpha S_{\mathrm{a}} D_{23}\right) v_{j, j i}+\rho_{\mathrm{a}} \omega^2 v_i=0, \end{aligned}$

(39)

$\begin{gathered} D_{11} u_{j, j i}+\rho_{\mathrm{w}} \omega^2 u_i+D_{12} w_{j, j i}+\frac{\rho_{\mathrm{w}}}{n S_{\mathrm{r}}} \omega^2 w_i-\frac{\rho_{\mathrm{w}} g}{k_{\mathrm{w}}} \mathrm{i} \omega w_i+ \\ D_{13} v_{j, j i}=0 \quad, \end{gathered}$

(40)

$\begin{gathered} D_{21} u_{j, j i}+\rho_{\mathrm{a}} \omega^2 u_i+D_{22} w_{j, j i}+\frac{\rho_{\mathrm{a}}}{n S_{\mathrm{a}}} \omega^2 v_i- \\ \frac{\rho_{\mathrm{a}} g}{k_{\mathrm{a}}} \mathrm{i} \omega v_i+D_{23} v_{j, j i}=0 \end{gathered}$

(41)

经过波数变换后整理得到的矩阵形式：

$\left[\begin{array}{ccc} \boldsymbol{K}_\boldsymbol{1}+\boldsymbol{K}_\boldsymbol{2}-\boldsymbol{M}_\boldsymbol{1} & \boldsymbol{L}_\boldsymbol{1}-\boldsymbol{M}_\boldsymbol{2} & \boldsymbol{G}_\boldsymbol{1}-\boldsymbol{M}_\boldsymbol{3} \\ \boldsymbol{K}_\boldsymbol{3}-\boldsymbol{M}_\boldsymbol{2} & \boldsymbol{L}_\boldsymbol{2}-\boldsymbol{M}_\boldsymbol{4} & \boldsymbol{G}_\boldsymbol{2} \\ \boldsymbol{K}_\boldsymbol{4}-\boldsymbol{M}_\boldsymbol{3} & \boldsymbol{L}_\boldsymbol{3} & \boldsymbol{G}_\boldsymbol{3}-\boldsymbol{M}_\boldsymbol{5} \end{array}\right]\left[\begin{array}{l} \boldsymbol{U} \\ \boldsymbol{W} \\ \boldsymbol{V} \end{array}\right]=\left[\begin{array}{l} \boldsymbol{F}^{\mathrm{\mathbf{s}}} \\ \boldsymbol{F}^{\mathrm{\mathbf{f}}} \\ \boldsymbol{F}^{\mathrm{\mathbf{g}}} \end{array}\right] 。$

(42)

式中：K，L，G分别为固相、液相和气相的刚度矩阵；M为各相的质量矩阵；U，W，V分别为固相、液相和气相的位移矩阵； $\boldsymbol{F}^\mathbf{s}，\boldsymbol{F}^\mathbf{f}，\boldsymbol{F}^\mathbf{g}$ 分别为固相、液相和气相的外力矢量矩阵。各矩阵的具体表达式见附录1。

求解时，首先建立y-z剖面的二维有限元模型，采用4节点等参单元划分模型网格。参考王勖成《有限单元法》^[28]，采用拉格朗日插值函数和高斯积分，根据式（42）得到单元刚度矩阵，最后集成总体刚度矩阵后进行求解。整个求解过程是根据式（42）在Matlab中通过编程实现。

2. 求解精度及求解效率

2.1 模型介绍

Yang^[27]讨论了网格尺寸与求解精度以及计算经济性的关系，认为计算区域2.5维有限元网格尺寸的最大值L(m)必须要小于等于 ${\lambda '_{\text{s}}}/6$ ，才能满足计算精度要求，网格尺寸的最小值大于等于0.5 ${\lambda '_{\text{s}}}$ 能兼顾计算的经济性；其中， ${{\lambda }^{\prime }}_{\text{s}}=2\text{π}/{{k}^{\prime }}_{\text{s}}，$ ${k'_{\text{s}}} = \sqrt {{k_{\text{s}}}^{\text{2}} - {{\left( {\frac{{\omega - {\omega _0}}}{c}} \right)}^2}}$ ， ${k_{\text{s}}}$ 为弹性半空间土体沿x方向上的剪切波波数， ${k'_{\text{s}}}$ 为移动荷载作用下土体沿x方向的剪切波波数， ${\lambda '_{\text{s}}}$ 为移动荷载作用下剪切波的波长，单位为m， ${\omega _0}$ 为移动荷载初始自振频率，c为外部荷载的移动速度。

本节建立了一个宽100 m，深50 m的2.5维有限元模型用以模拟非饱和半无限空间，用四边形等参单元模拟地基中的土单元，其中每个节点包含空间3个方向的位移自由度，如图 3所示，其中，由于x方向也被波数代替，因此只需在y-z平面进行网格划分。文献[16]中非饱和土的最小剪切波速为38 m/s，如外加荷载最大运行速度为200 m/s，则2.5维有限元网格的最大值不能超过0.2 m。因此，在建模时，本模型计算核心区域的网格尺寸为0.08 m×0.08 m，0.2 m×0.2 m；随后逐渐过渡至0.5 m×0.5 m，边缘网格尺寸为1.0 m×1.0 m。Bian等^[12]和胡静等^[13]在求解饱和地基动力响应时采用多层阻尼边界，很好地模拟了饱和半无限空间的动力响应。因此，本文同样采用多层阻尼边界对模型边界的波动进行吸收。整个2.5维有限元模型共计9052个网格，9250个节点。

图 3 模型网格分布示意图

Figure 3. Diagram of mesh distribution

下载: 全尺寸图片幻灯片

2.2 退化验证

（1）弹性地基验证

首先，将参数 $\alpha$ ， ${\rho _{\text{w}}}$ ， ${\rho _{\text{a}}}$ ，S_r和S_a的值都取无限小，则非饱和动力控制方程退化为单相弹性介质的Naiver波动方程。采用文献[12]中的单相弹性地基的参数（表 1所示）赋至上述2.5维有限元模型中进行计算。

表 1 单相弹性地基计算参数^[12]

Table 1. Parameters of single-phase elastic medium^[12]

K_s/GPa	K_w/GPa	K_a/kPa	G/MPa	$\nu$	n	ρ_s/(kg·m^-3)	ρ_w/(kg·m^-3)	ρ_a/(kg·m^-3)	S_r
35	2.25	145	25	0.125	0.001	2500	0.001	0.001	0.001

下载: 导出CSV

| 显示表格

图 4为单位荷载（1 N）以时速75 m/s移动时，观察点(x, 1, 1)处的位移响应时程曲线，x为荷载移动方向。由图 4可知，基于退化至弹性地基参数的2.5维非饱土模型计算得到的结果与半解析解结果吻合。

图 4 单相弹性地基验证

Figure 4. Verification of single-phase elastic medium

下载: 全尺寸图片幻灯片

（2）饱和地基验证

令S_r等于1，参数A_e，A_s，ρ_a和S_a的值都取无限小，非饱和土的动力控制方程退化为Biot波动方程。将文献[9]中双相饱和地基的参数（见表 2）赋至2.5维非饱和半空间模型中进行计算。

表 2 双相饱和地基计算参数^[9]

Table 2. Parameters of two-phase saturated medium^[9]

K_b/MPa	K_s/GPa	K_w/GPa	K_a/kPa	G/MPa	$\nu$	n	ρ_s/(kg·m^-3)	ρ_w/(kg·m^-3)	ρ_a/(kg·m^-3)	S_r	k_D/(m·s^-1)
8.77	11	2.25	145	2997.6	0.125	0.3	2500	1000	0.001	1	1×10^-6

下载: 导出CSV

| 显示表格

图 5为单位荷载（1 N）以121 m/s移动时观察点(x, 1, 1)处的位移及超静孔压响应。中，位移乘以 $G{a_{\text{R}}}{\text{/}}F$ ，超静孔压乘以 $a_{\text{R}}^2/F$ 进行了归一化，其中F=1 N，a_R=1.0 m，G为表 2中的剪切模量。从图 5中可以发现，基于退化至饱和地基参数的2.5维非饱土模型计算得到的位移响应和超静孔压响应都与半解析解吻合。

图 5 饱和地基验证

Figure 5. Verification of two-phase saturated medium

下载: 全尺寸图片幻灯片

2.3 与非饱和解析解验证

Lu等^[16]采用解析法求解了矩形荷载作用下非饱和半空间的动力响应，但该研究并未考虑变形对土-水特征曲线的影响。此处选取该研究中的参数进行计算，并与其结果进行对比以验证本模型对非饱和动力问题求解的准确性。模型中，荷载为0.32 m×0.32 m的矩形荷载，大小为30 kN，频率为2 Hz；地基土参数见表 3，地基上部为完全透水边界。由于非饱和土的剪切模量与饱和度有关，因而土体的剪切波波速也随着饱和度的变化而变化。为此，定义参数v_s^sat为不同饱和度下的土体剪切波波速，其上标sat代表饱和度。

表 3 非饱和地基计算参数^[16]

Table 3. Parameters of unsaturated medium^[16]

K_b/MPa	K_s/GPa	K_w/GPa	K_a/kPa	G_s/MPa	$\nu$	n	ρ_s/(kg·m^-3)	ρ_w/(kg·m^-3)	ρ_a/(kg·m^-3)
8.33	35	2.25	145	3.85	0.35	0.45	2650	1000	1.28

S_r	S_a	η_w/(Ns·m^-2)	η_a/(Ns·m^-2)	S_w0	κ/m²	β/(Pa^-1)	m	k	φ'
0.5，0.8	0.5，0.2	0.001	18×10^-6	0.05	5.3×10^-13	1.0×10^-4	0.5	2	10^°

下载: 导出CSV

| 显示表格

图 6为观察点（0, 0, 0.5）在不同饱和度和不同荷载移动速度下的响应。从图 6中可以看出，在低速和高速条件下，2.5维有限元在不同饱和度下的计算结果均与Lu等^[16]的解析解吻合良好。

图 6 非饱和地基验证

Figure 6. Verification of three-phase unsaturated medium

下载: 全尺寸图片幻灯片

2.2和2.3节共进行了3种介质模型的验证，数值结果都与已有解析/半解析结果吻合一致，证明了2.5维有限元方法求解多孔介质动力问题的准确性。

2.4 计算效率

采用Intel Core i9-10900k 10核20线程，基本频率为3.7 GHz的电脑对上述2.5维有限元模型进行计算。移动荷载为8节编组的列车荷载，移动速度为75 m/s，移动距离为384 m，结果表明，仅考虑每节点3个自由度的单相弹性介质模型所需的计算时间为60.43 s，考虑每节点6个自由度的双相饱和介质模型所需的计算时间为447.54 s，考虑每节点9个自由度的三相非饱和介质模型所需的计算时间为1807.24 s。可见，随着每节点自由度数的增加，模型所需的计算时间显著增大。Hall^[29]曾用ABAQUS建立了一个65 m×30 m的三维完全弹性地基模型，计算8节编组列车荷载下的动力响应，耗时约23 h。由此可见，采用有限元法求解多孔介质的动力响应问题不仅面临着二次开发困难，还必将存在计算用时过长的问题。与此同时，非饱和土波动方程求解较为复杂，将变形对土-水特征曲线的影响引入波动方程后，需要对关键控制项进行重新推导，解析和半解析求解难度增大。如果再考虑到路基/地基中饱和度横向及纵向的不均匀分布，解析和半解析方式则较难实现求解。对比而言，2.5维有限元方法体现出了高效的计算效率，与此同时还具备良好的准确性，因此，2.5维有限元方法是现阶段求解多孔介质动力问题的一种优势方法。

3. 数值分析

为研究考虑变形效应的SWCC对非饱和土动力响应的影响，基于孙德安^[23]研究，将孔隙比-饱和度直线图的斜率 ${\lambda _{{\text{se}}}}$ 取为-0.35，代入2.3节中的模型进行计算。对比了观察点（0, 0, 0.5）处，基于考虑变形与未考虑变形两种SWCC模型计算出的位移响应结果，如图 7所示。

图 7 SWCC考虑变形对位移响应的影响

Figure 7. Displacement responses under SWCC considering effect of deformation

下载: 全尺寸图片幻灯片

从图 7中可以看出，采用考虑变形效应的土-水特征曲线计算得到的位移响应大于未考虑变形的结果。根据图 2和式（3），（7）可知，考虑变形的SWCC相较于未考虑变形的曲线向右上方移动，导致同一基质吸力所对应的饱和度增大，土体剪切模量减小；因而在相同的荷载作用下其响应值更大。

图 8给出了多个饱和度和速度下，考虑变形效应所产生的位移增幅。

图 8 不同速度下考虑变形效应所产生的位移增幅随饱和度的变化

Figure 8. Amplification of deformation and saturation considering different velocities

下载: 全尺寸图片幻灯片

由图 8可知，当荷载移动速度没有超过土体的剪切波波速时，考虑变形带来的位移响应增幅与饱和度和列车速度均成正比。采用传统的土-水特征曲线模型去求解非饱和土的动力问题将低估土体的振动强度。

4. 结论

（1）本文在传统V-G模型的基础上，根据饱和度与孔隙比的关系，建立了考虑变形效应的土-水特征曲线模型，并据此推导了新的非饱和土动力控制方程，该方程能够完整描述动力作用下非饱土的水-力耦合作用。

（2）采用2.5维有限元方法对考虑SWCC变形效应的非饱和土动力控制方程进行了求解，求解结果分别与单相弹性介质，双相饱和介质和三相非饱和介质的解析/半解析结果进行对比，均验证了该求解方法的准确性。

（3）随着多孔介质模型的复杂程度即考虑相数的增加，动力求解模型的计算耗时会显著增长；2.5维有限元方法具备良好的准确性和高效的计算效率，并且能够根据实际地层特征进行建模，是目前求解多孔介质动力问题的一种优势算法。

（4）采用考虑变形效应的土-水特征曲线计算得到的位移响应大于未考虑变形的结果，说明传统的计算方法会低估非饱和土动力响应的强度。

（5）由于非饱和三相介质动力计算需要涉及众多参数，全套参数的测定也较耗时，因此文中进行的分析是基于前人已发表的参数。后续可进行非饱和土的试验研究，测定SWCC、渗透系数及强度，获得非饱和三相介质动力计算所需的参数，更加贴近实际应用。

附录

$\begin{aligned} & D_{11}=\left[A_{\mathrm{e}} n S_{\mathrm{a}}\left(M+M_{\mathrm{a}}\right)-\alpha C_1 S_{\mathrm{a}} M_{\mathrm{a}}+\alpha A_{\mathrm{s}}\right] /\left\{A_{\mathrm{s}} M+A_{\mathrm{s}}\left(S_{\mathrm{r}} M_{\mathrm{w}}+\right.\right. \\ & \left.\left.S_{\mathrm{a}} M_{\mathrm{a}}\right)+M S_{\mathrm{r}} S_{\mathrm{a}}\left[C_2 M_{\mathrm{w}}-C_1 M_{\mathrm{a}}\right]-\frac{1}{n} S_{\mathrm{r}} S_{\mathrm{a}} M_{\mathrm{w}} M_{\mathrm{a}}\right\} ; \\ & D_{12}=\left[C_2 S_{\mathrm{a}} M+\left(C_2-C_{\mathrm{r}}\right) S_{\mathrm{a}} M_{\mathrm{a}}+A_{\mathrm{s}}\right] /\left\{A_{\mathrm{s}} M+A_{\mathrm{s}}\left(S_{\mathrm{r}} M_{\mathrm{w}}+\right.\right. \\ & \left.\left.S_{\mathrm{a}} M_{\mathrm{a}}\right)+M S_{\mathrm{r}} S_{\mathrm{a}}\left[C_2 M_{\mathrm{w}}-C_1 M_{\mathrm{a}}\right]-\frac{1}{n} S_{\mathrm{r}} S_{\mathrm{a}} M_{\mathrm{w}} M_{\mathrm{a}}\right\} ; \\ & D_{13}=\left[C_1 S_{\mathrm{a}} M+A_{\mathrm{s}}\right] /\left\{A_{\mathrm{s}} M+A_{\mathrm{s}}\left(S_{\mathrm{r}} M_{\mathrm{w}}+S_{\mathrm{a}} M_{\mathrm{a}}\right)+\right. \\ & \left.M S_{\mathrm{r}} S_{\mathrm{a}}\left[C_2 M_{\mathrm{w}}-C_1 M_{\mathrm{a}}\right]-\frac{1}{n} S_{\mathrm{r}} S_{\mathrm{a}} M_{\mathrm{w}} M_{\mathrm{a}}\right\} ; \\ & D_{21}=\left[-A_{\mathrm{e}} n S_{\mathrm{r}}\left(M+M_{\mathrm{w}}\right)+\alpha C_2 S_{\mathrm{r}} M_{\mathrm{w}}+\alpha A_{\mathrm{s}}\right] /\left\{A_{\mathrm{s}} M+\right. \\ & \left.A_{\mathrm{s}}\left(S_{\mathrm{r}} M_{\mathrm{w}}+S_{\mathrm{a}} M_{\mathrm{a}}\right)+M S_{\mathrm{r}} S_{\mathrm{a}}\left[C_2 M_{\mathrm{w}}-C_1 M_{\mathrm{a}}\right]-\frac{1}{n} S_{\mathrm{r}} S_{\mathrm{a}} M_{\mathrm{w}} M_{\mathrm{a}}\right\} ; \\ & D_{22}=\left[-C_2 S_{\mathrm{r}} M+A_{\mathrm{s}}\right] /\left\{A_{\mathrm{s}} M+A_{\mathrm{s}}\left(S_{\mathrm{r}} M_{\mathrm{w}}+S_{\mathrm{a}} M_{\mathrm{a}}\right)+\right. \\ & \left.M S_{\mathrm{r}} S_{\mathrm{a}}\left[C_2 M_{\mathrm{w}}-C_1 M_{\mathrm{a}}\right]-\frac{1}{n} S_{\mathrm{r}} S_{\mathrm{a}} M_{\mathrm{w}} M_{\mathrm{a}}\right\} ; \\ & D_{23}=\left[-C_1 S_{\mathrm{r}} M-\left(C_1-C_2\right) S_{\mathrm{r}} M_{\mathrm{w}}+A_{\mathrm{s}}\right] /\left\{A_{\mathrm{s}} M+A_{\mathrm{s}}\left(S_{\mathrm{r}} M_{\mathrm{w}}+\right.\right. \\ & \left.\left.S_{\mathrm{a}} M_{\mathrm{a}}\right)+M S_{\mathrm{r}} S_{\mathrm{a}}\left[C_2 M_{\mathrm{w}}-C_1 M_{\mathrm{a}}\right]-\frac{1}{n} S_{\mathrm{r}} S_{\mathrm{a}} M_{\mathrm{w}} M_{\mathrm{a}}\right\} 。 \end{aligned}$

$\begin{aligned} & \boldsymbol{K}_1=\sum\limits_{\text {element }} \int\limits_{-1}^1 \int\limits_{-1}^1\left(\boldsymbol{B}^* \boldsymbol{N}\right)^{\mathrm{T}} \boldsymbol{D} \boldsymbol{B} \boldsymbol{N}|\mathbf{J}| \mathrm{d} \eta \mathrm{~d} \varsigma, \\ & \boldsymbol{K}_2=\left(\alpha S_{\mathrm{r}} D_{11}+\alpha S_{\mathrm{a}} D_{21}\right) \sum\limits_{\text {element }} \int\limits_{-1}^1 \int\limits_{-1}^1\left(\boldsymbol{B}^* \boldsymbol{N}\right)^{\mathrm{T}} m m^{\mathrm{T}} \boldsymbol{B} \boldsymbol{N}|\mathbf{J}| \mathrm{d} \eta \mathrm{~d} \varsigma, \\ & \boldsymbol{K}_3=D_{11} \sum\limits_{\text {element }} \int\limits_{-1}^1 \int\limits_{-1}^1\left(\boldsymbol{B}^* \boldsymbol{N}\right)^{\mathrm{T}} m m^{\mathrm{T}} \boldsymbol{B} \boldsymbol{N}|\mathbf{J}| \mathrm{d} \eta \mathrm{~d} \varsigma, \\ & \boldsymbol{K}_4=D_{21} \sum\limits_{\text {element }} \int\limits_{-1}^1 \int\limits_{-1}^1\left(\boldsymbol{B}^* \boldsymbol{N}\right)^{\mathrm{T}} m m^{\mathrm{T}} \mathbf{B} \mathbf{N}|\mathbf{J}| \mathrm{d} \eta \mathrm{~d} \varsigma。 \end{aligned}$

$\begin{aligned} & \boldsymbol{L}_1=\left(\alpha S_{\mathrm{r}} D_{12}+\alpha S_{\mathrm{a}} D_{22}\right) \sum\limits_{\text {element }} \int\limits_{-1}^1 \int\limits_{-1}^1\left(\boldsymbol{B}^* \boldsymbol{N}\right)^{\mathrm{T}} m m^{\mathrm{T}} \boldsymbol{B} \boldsymbol{N}|\mathbf{J}| \mathrm{d} \eta \mathrm{~d} \varsigma, \\ & \boldsymbol{L}_2=D_{12} \sum\limits_{\text {element }} \int\limits_{-1}^1 \int\limits_{-1}^1\left(\boldsymbol{B}^* \boldsymbol{N}\right)^{\mathrm{T}} m m^{\mathrm{T}} \boldsymbol{B} \boldsymbol{N}|\mathbf{J}| \mathrm{d} \eta \mathrm{~d} \varsigma, \\ & \boldsymbol{L}_3=D_{22} \sum\limits_{\text {element }} \int\limits_{-1}^1 \int\limits_{-1}^1\left(\boldsymbol{B}^* \boldsymbol{N}\right)^{\mathrm{T}} m m^{\mathrm{T}} \boldsymbol{B} \boldsymbol{N}|\mathbf{J}| \mathrm{d} \eta \mathrm{~d} \varsigma 。 \end{aligned}$

$\begin{aligned} & \boldsymbol{G}_1=\left(\alpha S_{\mathrm{r}} D_{13}+\alpha S_{\mathrm{a}} D_{23}\right) \sum\limits_{\text {element }} \int\limits_{-1}^1 \int\limits_{-1}^1\left(\boldsymbol{B}^* \boldsymbol{N}\right)^{\mathrm{T}} m m^{\mathrm{T}} \boldsymbol{B} \boldsymbol{N}|\mathbf{J}| \mathrm{d} \eta \mathrm{~d} \varsigma, \\ & \boldsymbol{G}_2=D_{13} \sum\limits_{\text {element }} \int\limits_{-1}^1 \int\limits_{-1}^1\left(\boldsymbol{B}^* \boldsymbol{N}\right)^{\mathrm{T}} m m^{\mathrm{T}} \boldsymbol{B} \boldsymbol{N}|\mathbf{J}| \mathrm{d} \eta \mathrm{~d} \varsigma, \\ & \boldsymbol{G}_3=D_{23} \sum\limits_{\text {element }} \int\limits_{-1}^1 \int\limits_{-1}^1\left(\boldsymbol{B}^* \boldsymbol{N}\right)^{\mathrm{T}} m m^{\mathrm{T}} \boldsymbol{B} \boldsymbol{N}|\mathbf{J}| \mathrm{d} \eta \mathrm{~d} \varsigma 。 \end{aligned}$

$\begin{aligned} & \boldsymbol{M}_1=\omega^2 \rho \sum\limits_{\text {element }} \int\limits_{-1}^1 \int\limits_{-1}^1 \boldsymbol{N}^{\mathrm{T}} \boldsymbol{N}|\mathbf{J}| \mathrm{d} \eta \mathrm{~d} \varsigma, \\ & \boldsymbol{M}_2=\omega^2 \rho_{\mathrm{w}} \sum\limits_{\text {element }} \int\limits_{-1}^1 \int\limits_{-1}^1 \boldsymbol{N}^{\mathrm{T}} \boldsymbol{N}|\mathbf{J}| \mathrm{d} \eta \mathrm{~d} \varsigma, \\ & \boldsymbol{M}_3=\omega^2 \rho_{\mathrm{a}} \sum\limits_{\text {element }} \int\limits_{-1}^1 \int\limits_{-1}^1 \boldsymbol{N}^{\mathrm{T}} \boldsymbol{N}|\mathbf{J}| \mathrm{d} \eta \mathrm{~d} \varsigma, \\ & \boldsymbol{M}_4=\left(\omega^2 \frac{\rho_{\mathrm{w}}}{n S_{\mathrm{r}}}-\mathrm{i} \omega \frac{\rho_{\mathrm{w}} g}{k_{\mathrm{w}}}\right) \sum\limits_{\text {element }} \int\limits_{-1}^1 \int\limits_{-1}^1 \boldsymbol{N}^{\mathrm{T}} \boldsymbol{N}|\mathbf{J}| \mathrm{d} \eta \mathrm{~d} \zeta 。 \end{aligned}$

$\boldsymbol{D}=\left(\begin{array}{cccccc} \lambda+2 \mu & \lambda & \lambda & 0 & 0 & 0 \\ \lambda & \lambda+2 \mu & \lambda & 0 & 0 & 0 \\ \lambda & \lambda & \lambda+2 \mu & 0 & 0 & 0 \\ 0 & 0 & 0 & \mu & 0 & 0 \\ 0 & 0 & 0 & 0 & \mu & 0 \\ 0 & 0 & 0 & 0 & 0 & \mu \end{array}\right),$

$\boldsymbol{B}^{\mathbf{T}}=\left[\begin{array}{cccccc} -\mathrm{i} \xi_x & 0 & 0 & \frac{\partial}{\partial_y} & 0 & \frac{\partial}{\partial_z} \\ 0 & \frac{\partial}{\partial_y} & 0 & -\mathrm{i} \xi_x & \frac{\partial}{\partial_z} & 0 \\ 0 & 0 & \frac{\partial}{\partial_z} & 0 & \frac{\partial}{\partial_y} & -\mathrm{i} \xi_x \end{array}\right],$

$\boldsymbol{m}=\left(\begin{array}{l} 1 \\ 1 \\ 1 \\ 0 \\ 0 \\ 0 \end{array}\right)|\mathbf{J}|=\left|\begin{array}{ll} \sum\limits_{i=1}^4 \frac{\partial N_i}{\partial \eta} y_i^{\text {element }} & \sum\limits_{i=1}^4 \frac{\partial N_i}{\partial \xi} y_i^{\text {element }} \\ \sum\limits_{i=1}^4 \frac{\partial N_i}{\partial \eta} z_i^{\text {element }} & \sum\limits_{i=1}^4 \frac{\partial N_i}{\partial \xi} z_i^{\text {element }} \end{array}\right| 。$

式中：B为应变矩阵；D为弹性矩阵； $N$ 为插值函数矩阵； ${\mathbf{J}}$ 为Jocobi矩阵； $\eta ,\zeta$ 表示局部坐标系中的变量；y_i，z_i为总体坐标系中的坐标值。

图 1 考虑拉应力截断的非线性Power-Law强度准则

Figure 1. Nonlinear Power-Law criterion considering tension cut-off

下载: 全尺寸图片幻灯片

图 2 隧道开挖面多曲线圆锥破坏模式

Figure 2. Multi-cone failure mechanism of tunnel face

下载: 全尺寸图片幻灯片

图 3 考虑拉应力截断修正的多锥破坏机构

Figure 3. Modified multi-cone failure mechanism considering tension cut-off

下载: 全尺寸图片幻灯片

${c_0}$ 对临界支护力的影响

Effects of ${c_0}$ on critical support pressure

下载: 全尺寸图片幻灯片

图 5 隧道直径对临界支护力的影响

Figure 5. Effects of tunnel diameter on critical support pressure

下载: 全尺寸图片幻灯片

图 6 σ_t对临界支护力的影响

Figure 6. Effects of σ_t on critical support pressure

下载: 全尺寸图片幻灯片

图 7 m对临界支护力的影响

Figure 7. Effects of m on critical support pressure

下载: 全尺寸图片幻灯片

图 8 不同非线性系数反算的应力分量

Figure 8. Back-calculated stress components for different nonlinear coefficients

下载: 全尺寸图片幻灯片

图 9 不同单轴抗拉强度反算的应力分量

Figure 9. Back-calculated stress components for different uniaxial tensile strengths

下载: 全尺寸图片幻灯片

图 10 不同拉应力截断系数反算的应力分量

Figure 10. Back-calculated stress components for different ξ

下载: 全尺寸图片幻灯片

图 11 非线性系数对多锥破坏机构的影响

Figure 11. Effects of nonlinear coefficient on geometry of multi-cone failure mechanism

下载: 全尺寸图片幻灯片

图 12 单轴抗拉强度对多锥破坏机构的影响

Figure 12. Effects of uniaxial tensile strength on geometry of multi-cone failure mechanism

下载: 全尺寸图片幻灯片

图 13 拉应力截断程度对多锥破坏机构的影响

Figure 13. Effects of T-C on geometry of multi-cone failure.mechanism

下载: 全尺寸图片幻灯片

图 14 考虑拉应力截断情况下隧道开挖面多曲线圆锥破坏机构

Figure 14. Multi-cone failure mechanism of tunnel faces considering T-C

下载: 全尺寸图片幻灯片

图 15 不同非线性参数对开挖面安全系数的影响

Figure 15. Effects of different nonlinear parameters on safety factors of the tunnel faces

下载: 全尺寸图片幻灯片

表 1 临界支护力与理论解对比

Table 1 Comparison between critical support pressures with theoretical solutions 单位: kPa

工况	${c_0}$ /kPa	$\varphi$ /(°)	${\sigma _{\text{t}}}$ /kPa	破坏准则	Li等^[16]	Mollon等^[3]	本文结果	差异/%
情况1	7	17	22.9	T-C（ξ=0）	33.92	37.44	33.84	9.6
情况1	7	17	22.9	M-C	33.61	37.07	33.53	9.5
情况2	10	25	21.4	T-C（ξ=0）	11.34	14.23	11.32	20.4
情况2	10	25	21.4	M-C	10.63	13.71	10.74	21.7
情况3	15	15	55.0	T-C（ξ=0）	14.01	18.22	13.97	23.3
情况3	15	15	55.0	M-C	10.77	15.20	10.72	29.5

下载: 导出CSV

表 2 开挖面临界支护力对比

Table 2 Comparison of critical support pressures of tunnel faces for σ_t=20 kPa 单位: kPa

c₀	m = 1		m = 1.2		m = 1.4		m = 1.6
c₀	T-C	M-C	T-C	P-L	T-C	P-L	T-C	P-L
4	75.29	74.88	123.81	123.23	192.81	191.97	291.53	291.01
6	38.03	37.81	62.30	61.86	95.35	94.57	138.84	137.99
8	19.96	19.64	33.10	32.67	51.46	50.73	75.22	74.11
10	9.80	9.23	17.01	16.42	27.33	26.57	41.25	40.13
12	3.65	2.57	7.56	6.56	13.05	12.09	20.91	19.75
14	—	—	1.92	0.10	4.70	3.01	8.52	7.01
16	—	—	—	—	0.03	—	1.69	—
注：—表示不需要支护。

下载: 导出CSV

表 3 线性条件下开挖面安全系数对比

Table 3 Comparison of safety factors of tunnel faces under linear condition

c/kPa	φ/(°)	M-C	ξ = 0	Vermeer等^[18]	Pan等^[19]
15	15	0.84	0.75	0.87	0.81
15	20	0.92	0.82	0.91	0.87
20	25	1.33	1.20	1.21	1.15
20	15	1.12	0.99	1.12	1.02
25	20	1.53	1.36	1.41	1.30
25	25	1.67	1.52	1.46	1.36

下载: 导出CSV

参考文献(19)

[1]	PERAZZELLI P, LEONE T, ANAGNOSTOU G. Tunnel face stability under seepage flow conditions[J]. Tunnelling and Underground Space Technology, 2014, 43: 459-469. doi: 10.1016/j.tust.2014.03.001
[2]	SAADA Z, MAGHOUS S, GARNIER D. Stability analysis of rock slopes subjected to seepage forces using the modified Hoek-Brown criterion[J]. International Journal of Rock Mechanics and Mining Sciences, 2012, 55: 45-54. doi: 10.1016/j.ijrmms.2012.06.010
[3]	MOLLON G, DIAS D, SOUBRA A H. Rotational failure mechanisms for the face stability analysis of tunnels driven by a pressurized shield[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 2011, 35(12): 1363-1388. doi: 10.1002/nag.962
[4]	吕玺琳, 王浩然, 黄茂松. 盾构隧道开挖面稳定极限理论研究[J]. 岩土工程学报, 2011, 33(1): 57-62. LÜ Xilin, WANG Haoran, HUANG Maosong. Limit theoretical study on face stability of shield tunnels[J]. Chinese Journal of Geotechnical Engineering, 2011, 33(1): 57-62. (in Chinese)
[5]	徐前卫, 唐卓华, 朱合华, 等. 盾构隧道开挖面极限支护压力研究[J]. 岩土工程学报, 2017, 39(7): 1234-1240. doi: 10.11779/CJGE201707009 XU Qianwei, TANG Zhuohua, ZHU Hehua, et al. Limit support pressure at excavation face of shield tunnels[J]. Chinese Journal of Geotechnical Engineering, 2017, 39(7): 1234-1240. (in Chinese) doi: 10.11779/CJGE201707009
[6]	黄阜, 潘秋景, 凌同华. 地下水渗流作用下的盾构隧道开挖面安全系数上限分析[J]. 岩土工程学报, 2017, 39(8): 1461-1469. doi: 10.11779/CJGE201708013 HUANG Fu, PAN Qiujing, LING Tonghua. Upper bound analysis of factor of safety for shield tunnel face subjected to underground water seepage[J]. Chinese Journal of Geotechnical Engineering, 2017, 39(8): 1461-1469. (in Chinese) doi: 10.11779/CJGE201708013
[7]	黄阜, 王芬, 张芝齐, 等. 基于Hoek-Brown破坏准则的盾构隧道开挖面支护力上限研究[J]. 铁道科学与工程学报, 2018, 15(11): 2892-2900. HUANG Fu, WANG Fen, ZHANG Zhiqi, et al. Upper bound solution of supporting pressure for shield tunnel face subjected to Hoek-Brown failure criterion[J]. Journal of Railway Science and Engineering, 2018, 15(11): 2892-2900. (in Chinese)
[8]	金大龙, 袁大军, 郑浩田, 等. 高水压条件下泥水盾构开挖面稳定离心模型试验研究[J]. 岩土工程学报, 2019, 41(9): 1653-1660. doi: 10.11779/CJGE201909009 JIN Dalong, YUAN Dajun, ZHENG Haotian, et al. Centrifugal model tests on face stability of slurry shield tunnels under high water pressures[J]. Chinese Journal of Geotechnical Engineering, 2019, 41(9): 1653-1660. (in Chinese) doi: 10.11779/CJGE201909009
[9]	徐涛, 史庆锋, 章定文, 等. 泥水盾构开挖面泥膜渗透特性与压力传递机制[J]. 岩土工程学报, 2023, 45(9): 1878-1887. doi: 10.11779/CJGE20220866 XU Tao, SHI Qingfeng, ZHANG Dingwen, et al. Permeability characteristics of filter cake and pressure transfer on face during slurry shield tunnelling[J]. Chinese Journal of Geotechnical Engineering, 2023, 45(9): 1878-1887. (in Chinese) doi: 10.11779/CJGE20220866
[10]	陈仁朋, 尹鑫晟, 李育超, 等. 泥水盾构泥膜渗透性及其对开挖面稳定性影响[J]. 岩土工程学报, 2017, 39(11): 2102-2108. doi: 10.11779/CJGE201711018 CHEN Renpeng, YIN Xinsheng, LI Yuchao, et al. Permeability of filter cake and its influence on face stability of slurry shield-driven tunnels[J]. Chinese Journal of Geotechnical Engineering, 2017, 39(11): 2102-2108. (in Chinese) doi: 10.11779/CJGE201711018
[11]	DRUCKER D C, PRAGER W. Soil mechanics and plastic analysis or limit design[J]. Quarterly of Applied Mathematics, 1952, 10(2): 157-165. doi: 10.1090/qam/48291
[12]	PAUL B. A modification of the Coulomb-Mohr theory of fracture[J]. Journal of Applied Mechanics, 1961, 28(2): 259-268. doi: 10.1115/1.3641665
[13]	HOEK E, MARTIN C D. Fracture initiation and propagation in intact rock: a review[J]. Journal of Rock Mechanics and Geotechnical Engineering, 2014, 6(4): 287-300. doi: 10.1016/j.jrmge.2014.06.001
[14]	CHARLES J A, SOARES M M. The stability of slopes in soils with nonlinear failure envelopes[J]. Canadian Geotechnical Journal, 1984, 21(3): 397-406. doi: 10.1139/t84-044
[15]	CHEN W F. Limit Analysis and Soil Plasticity[M]. Amsterdam: Elsevier Scientific Pub Co, 1975.
[16]	LI T Z, YANG X L. Three-dimensional face stability of shallow-buried tunnels with tensile strength cut-off[J]. Computers and Geotechnics, 2019, 110: 82-93. doi: 10.1016/j.compgeo.2019.02.014
[17]	MICHALOWSKI R L. Failure potential of infinite slopes in bonded soils with tensile strength cut-off[J]. Canadian Geotechnical Journal, 2018, 55(4): 477-485. doi: 10.1139/cgj-2017-0041
[18]	VERMEER P A, RUSE N, MARCHER T. Tunnel heading stability in drained ground[J]. Felsbau, 2002, 20(6): 8-18.
[19]	PAN Q J, DIAS D. Upper-bound analysis on the face stability of a non-circular tunnel[J]. Tunnelling and Underground Space Technology, 2017, 62: 96-102. doi: 10.1016/j.tust.2016.11.010