基于改进对流粒子域插值物质点法的隧道大变形分析

王曼灵; 李树忱; 周慧颖; 王修伟; 彭科峰; 袁超

doi:10.11779/CJGE20230676

基于改进对流粒子域插值物质点法的隧道大变形分析 English Version

王曼灵^{1, 2,},
李树忱^{1, 2, 3, ,},
周慧颖^{1, 2},
王修伟^{1, 2},
彭科峰^{1, 2},
袁超^{2, 4}

1.
山东大学齐鲁交通学院, 山东济南 250002
2.
山东大学岩土与地下工程研究院, 山东济南 250061
3.
中国矿业大学力学与土木工程学院, 江苏徐州 221116
4.
山东大学土建与水利学院, 山东济南 250061

基金项目:

国家重点研发计划项目 2021YFC2902103

国家自然科学基金项目 52108373

山东省自然科学基金项目 ZR2021ZD36

山东省自然科学基金项目 ZR2021QE127

详细信息

作者简介:
王曼灵(1994—)，女，博士研究生，主要从事隧道大变形和数值计算方法等方面的研究工作。E-mail: manling_wang@163.com

通讯作者:
李树忱, E-mail: shuchenli@sdu.edu.cn

中图分类号: TU434
计量
- 文章访问数: 357
- HTML全文浏览量: 32
- PDF下载量: 69
出版历程
- 收稿日期: 2023-07-17
- 网络出版日期: 2023-12-19
- 刊出日期: 2024-07-31

Improved convective particle domain interpolation material point method for large deformation analysis of tunnels

WANG Manling^{1, 2,},
LI Shuchen^{1, 2, 3, ,},
ZHOU Huiying^{1, 2},
WANG Xiuwei^{1, 2},
PENG Kefeng^{1, 2},
YUAN Chao^{2, 4}

1.
School of Qilu Transportation, Shandong University, Jinan 250002, China
2.
Institue of Geotechnical and Underground Engineering, Shandong University, Jinan 250061, China
3.
School of Mechanics and Civil Engineering, China University of Mining and Technology, Xuzhou 221116, China
4.
School of Civil Engineering, Shandong University, Jinan 250061, China

摘要

摘要: 物质点法（MPM）在模拟大变形问题时具有很好的效果，然而传统的MPM在粒子穿越网格边界时存在单元穿越误差，导致精度降低。为克服传统MPM的单元穿越误差，基于对流粒子域插值物质点法（CPDI）理论框架，采用自适应正交改进插值移动最小二乘法（AOIIMLS），提出了改进CPDI方法。AOIIMLS通过构造加权正交基函数，并且忽略了新对角矩阵中的零元素或极小元素的贡献，以避免求解逆矩阵，增强了鲁棒性。改进CPDI采用速度梯度计算粒子域的速度场，粒子速度和粒子域角点速度用于重构背景网格速度函数。通过一维柱在自重作用下的压缩、砂柱坍塌和隧道坍塌离心机试验验证了改进CPDI方法的准确性和适用性，结果表明改进CPDI降低了单元穿越误差，得到了更高的精度。最后，采用改进CPDI方法模拟了青岛地铁4号线静沙区间地面塌陷全过程，验证了改进CPDI方法在岩土工程大变形领域的适用性及优势。
- 物质点法 /
- 对流粒子域插值 /
- 自适应正交改进移动最小二乘法 /
- 隧道大变形
Abstract: The material point method (MPM) has good effects in simulating large deformation problems. However, the conventional MPM suffers from cell-crossing errors when particles cross grid boundaries, resulting in reduced accuracy. In order to overcome the cell-crossing errors of the conventional MPM, an improved convective particle domain interpolation material point method (CPDI) is proposed based on the conventional CPDI framework and the adaptive orthogonal improved interpolation moving least squares method (AOIIMLS). By constructing weighted orthogonal basis functions and disregarding the minimal or zero elements in the new diagonal matrix, the inverse matrix computation is avoided, and the robustness is enhanced. In the improved CPDI method, the particle domain velocity field is calculated using the velocity gradients, and the AOIIMLS shape functions are employed to reconstruct the background grid velocity function using the particle velocity and particle domain corner point velocity. The accuracy and applicability of the improved CPDI method are verified through simulations of various scenarios such as the compaction of a one-dimensional column under self-weight, the collapse of a sand column and the centrifuge tests on tunnel collapse. The results show that the improved CPDI method reduces the cell-crossing errors caused by the particles cross grid boundaries and achieves higher accuracy. Finally, the improved CPDI method is employed to simulate the whole process of ground collapse in the Jinggang Road Station–Shazikou Station tunnel section of Qingdao Metro Line 4, effectively confirming the applicability and advantages of the method in addressing large deformation problems in geotechnical engineering.

HTML全文

0. 引言

维持开挖面稳定性是确保盾构施工安全的关键，在盾构施工中若开挖面支护力等施工参数设置不当，将会诱发开挖面失稳，造成地表的沉降^[1]。近年来，中国城市地下工程建设迅速发展，出现了越来越多的新建隧道近接既有隧道施工案例^[2]。盾构隧道近接施工的力学行为更加复杂，施工风险更加突出。与常规单一隧道不同，近接施工中的盾构开挖面稳定性分析受既有结构、近接形式和近接距离等因素的影响，因而呈现出多体、多形式、多过程的特点。当前针对近接施工盾构开挖面稳定性分析的研究成果匮乏，理论依据尚不完善，给盾构近接施工带来较大风险。

盾构隧道开挖面稳定性问题已成为地下工程施工安全研究的重点。国内外学者从模型试验、理论分析和数值模拟等方面开展了大量的研究。模型试验是研究盾构开挖面失稳破坏的主要手段之一，基于离心模型试验、大尺寸模型试验和小尺寸透明土模型试验等方法，学者们得到了砂土地层的烟囱状破坏模式^[3]和黏土地层的盆状破坏模式^[4]。在理论分析方面，学者们基于极限状态分析方法(极限平衡法和极限分析法)提出了二维及三维的失稳破坏模型，包括二维锥形模型^[5]、二维对数螺旋线模型^[6]、三维对数螺旋体模型^[7]、经典三维楔形体模型^[8]、各类改进的棱柱体-楔形体模型^[9-11]、Torus模型^[12]和三维水平球缺形模型^[13]等。二维模型和楔形体模型由于过于简化，因而计算结果误差较大，而三维旋转破坏理论考虑了破坏体三维空间形态，并且可用于砂土和黏土地层，因而具有明确的物理意义，准确性高，适用性广。近年来已有研究考虑不同影响因素对三维旋转破坏理论进行了改进，如Zou等^[14]考虑土拱效应影响，将对数螺旋线旋转模型在隧道拱顶以上部分修正为半椭球体；Ibrahim等^[15]考虑开挖面位于成层土地基，采用分段的对数螺旋形进行修正；Li等^[16]考虑渗流的影响，将水头方程引入三维旋转模型。然而，考虑近接既有隧道条件下的盾构开挖面稳定性问题及其相关计算方法还鲜见报道，因而本文基于近接既有隧道条件下开挖面破坏模式，对三维旋转模型的截面形状进行修正，可用于深入分析近接施工的盾构开挖面稳定性。

数值模拟是确定盾构隧道开挖面破坏模式的主要手段之一，主要包括有限元法(FEM)、有限差分法(FDM)和离散元法(DEM)等。有限元法在计算大变形时可能会出现网格畸变导致计算不收敛。离散元法虽然对大变形问题处理较好，但计算效率较低，并且对于足尺模型无法还原土体的真实粒径。有限差分法作为隧道围岩变形计算的常用方法，具有求解速度快，计算精度高和可处理大变形等优势而被广泛应用。

本文基于FLAC^3D有限差分软件揭示了近接既有隧道条件下的盾构开挖面破坏模式，构建了考虑水平方向破坏扩展不均匀性的对数螺旋体+椭圆柱的开挖面稳定性三维理论分析模型，给出了考虑土体黏聚力和重度影响系数的开挖面极限支护力表达式。分析了近接距离、隧道埋深和地层参数等因素对近接既有隧道条件下的盾构开挖面极限支护力的影响。

1. FLAC^3D数值模拟

1.1 数值模型建立

对于单隧道条件，可以建立半模型以节约计算资源，但本文主要考虑近接施工条件，若采用半模型会截断既有隧道，不利于设置边界条件，因此采用全模型进行数值模拟分析。为尽可能地削减边界效应，模型尺寸取为10倍盾构直径(D = 6 m)，如图 1所示。

图 1 FLAC^3D数值模型

Figure 1. Numerical model in FLAC^3D

下载: 全尺寸图片幻灯片

模型边界条件设置为：地表为自由边界，侧面为固定水平位移边界，底面为固定边界。土体选用基于Mohr–Coulomb准则的弹塑性本构模型。数值模型物理力学参数如表 1所示。不同地层条件下的土体参数根据Mollon等^[17]的建议取为4种，分别对应松砂、密砂、软黏和硬黏。土体的弹性模量根据不同土性分类进行取值^[18]。隧道结构(新建隧道与既有隧道)均采用衬砌结构单元(Liner)，厚度取为0.35 m。注浆层采用泥浆的相关参数。利用六面体单元划分网格，模型建立完成后，共有79360个单元，83681个节点。

表 1 数值模型参数(引自文献[17-18])

Table 1. Parameters of numerical model (data from References [17-18])

类别		重度γ/(kN·m^-3)	弹性模量E/MPa	泊松比 $\nu$	黏聚力c/kPa	内摩擦角φ/(°)
土体	松砂	17.0	17.5	0.397	0	20
	密砂	19.5	65.0	0.263	0	40
	软黏	16.0	3.5	0.414	7	17
	硬黏	18.0	12.5	0.366	10	25
衬砌	混凝土	25.0	33500	0.20	—	—
注浆层	泥浆	12.0	50	0.49	—	—

下载: 导出CSV

| 显示表格

模型生成后，首先施加重力加速度，使模型整体达到初始地应力平衡状态。其次，开挖既有隧道并施加衬砌，并再次进行地应力平衡。随后，开挖新建隧道至既有隧道附近并在开挖面处施加与初始水平地应力持平的支护力。本文针对盾构开挖面失稳破坏这一特定阶段进行研究，因此简化盾构施工过程，采用一次性开挖到位的形式。最后采用逐渐降低支护力的手段模拟开挖面失稳，这是由于在实际工程中，开挖面失稳往往是由于支护力不足导致的，即主动破坏模式。通过逐级降低开挖面支护力(共取20级)，待每级支护力下模型运算至平衡状态或出现塑性流动后方可调整为下一级支护力^[19]，直至支护力降为0。具体模拟过程如图 2所示。近接条件下由于设定既有隧道在新建隧道上方，限制了土体的变形发展，因而此时破坏模式的主要影响因素应为隧道间距，而非埋深，故在近接条件下暂不考虑埋深的影响。具体模拟方案如表 2所示。

图 2 数值模拟过程

Figure 2. Process of numerical simulations

下载: 全尺寸图片幻灯片

表 2 数值模拟方案

Table 2. Scheme of numerical simulations

土体类别	近接形式	新建隧道埋深比C/D	隧道垂直净距比d/D	开挖面水平间距比s/D
松砂	无	1.0，1.5，2.0	—	—
密砂	无	1.0，1.5，2.0	—	—
松砂	新建隧道垂直下穿既有隧道	2.0	0.25，0.50，0.75	-1，-0.5，0，0.5
密砂
软黏
硬黏
注：隧道垂直净距比d/D为新建隧道上表面与既有隧道下表面的垂直距离d与新建隧道直径D之比；开挖面水平间距比s/D为新建隧道开挖面与既有隧道中心轴线的水平距离s与新建隧道直径D之比，负值表示盾构尚未穿越既有隧道中心轴线。

下载: 导出CSV

| 显示表格

在模拟过程中，由于盾构开挖面支护力逐渐降低，开挖面土体产生朝向隧道内部的位移，但位移分布并不均匀，为了统一量化开挖面位移随支护力的变化情况，通过FISH函数获取开挖面各节点位移的平均值，以此作为盾构开挖面的代表性位移。

1.2 数值模型建立

为了验证数值模拟中采用的本构模型和参数的适用性和正确性，选取Kirsch等^[3]开展的离心模型试验结果进行对比分析。图 3所示为单隧道条件下，埋深比C/D = 1.5时的变形区域对比。可以看出，离心模型试验中沿开挖面前方和上方的主要变形范围与本文数值模拟结果基本一致，结果相差均在5%以内，因此可以认为数值模型可以较好地还原真实情况。

图 3 失稳变形区域对比^[3]

Figure 3. Comparison of instability areas^[3]

下载: 全尺寸图片幻灯片

图 4为不同埋深条件下开挖面支护力-位移关系对比结果，其中，虚线为本文单隧道条件下的开挖面极限支护力，取值原则为：当支护力变化不大时，开挖面位移显著增加，认为开挖面发生失稳，此时的支护力即为极限支护力。通过对比可以发现，松砂条件下数值模拟结果和离心试验结果^[3]相差不大；而密砂条件下，数值模拟结果未显现出支护力曲线回弹的现象，这是由于离心试验^[3]中采用“位移法”(控制开挖面的支护面板以一定速度后退模拟失稳的方法)，初始阶段由于土拱效应使得支护力迅速下降，超过临界位移后，周围土体出现“应变软化”，导致所需支护力回升^[20]。而在实际盾构施工中，开挖面失稳的直接原因是支护力的不足，因而采用“力法”(即本文中控制支护力逐级降低的方法)模拟开挖面失稳更加符合实际。

图 4 无量纲化开挖面支护力–位移关系曲线

Figure 4. Relationship between normalized displacement and normalized support force of tunnel face

下载: 全尺寸图片幻灯片

1.3 数值模拟结果

近接条件下(松砂，d/D = 0.5)盾构开挖面失稳变形模式如图 5所示。由于黏性土层失稳滑移时没有明显的剪切带，因此多通过采用位移轮廓线来确定破坏模式的范围^[21]。本文参考Li等^[22]数值模拟结果，取位移值大于1%盾构直径的区域为破坏模式包络面。由模拟结果可见，随着盾构开挖面接近→穿越→远离既有隧道(s/D逐渐增加)，失稳区呈现出沿开挖方向不对称的特点，这是由于既有隧道与土体刚度差异很大，限制了土体的位移发展。除s/D =-0.5外，其余3种条件下开挖面变形模式均与单隧道情况类似，即由开挖面前方对数螺旋体和上方筒仓状滑动体组合而成。当s/D =-0.5时，开挖面变形呈现出沿既有隧道方向扩展的形态。

图 5 近接既有隧道条件下开挖面变形模式(不同s/D)

Figure 5. Deformation characteristics of shield tunnel face adjacent to an existing tunnel under different s/D

下载: 全尺寸图片幻灯片

不同隧道垂直净距下(松砂，s/D = -0.5)的盾构开挖面失稳模式如图 6所示。垂直净距对破坏模式的影响较小，随着净距的增加，上方破坏影响范围有所增加。不同土性条件下的开挖面失稳模式(d/D = 0.5，s/D = -0.5)如图 7所示，可以发现不同土性中的失稳模式基本一致，只是沿开挖面前方和两侧的扩展范围有所不同：密砂相比松砂向两侧扩展范围更小，黏土相比砂土向前方扩展更大。

图 6 近接既有隧道条件下开挖面变形模式(不同d/D)

Figure 6. Deformation characteristics of shield tunnel face adjacent to an existing tunnel under different d/D

下载: 全尺寸图片幻灯片

图 7 近接既有隧道条件下开挖面变形模式(不同土性)

Figure 7. Deformation characteristics of shield tunnel face adjacent to an existing tunnel under different soil properties

下载: 全尺寸图片幻灯片

近接条件下的开挖面极限支护力如图 8所示。图 8(a)中虚线为相应地层中的单隧道开挖面极限支护力，可以看出既有隧道的存在使开挖面极限支护力有所降低，这是由于既有隧道的刚度远大于土体，在开挖面失稳过程中起到了“支撑”的作用。对于同种地层，随着隧道垂直净距的增加，开挖面极限支护力略有增加，这是由于隧道净距越大，隧道间松动土柱越高，导致松动土压力越大，也就需要更大的支护力。图 8(b)为盾构穿越既有隧道过程中的开挖面极限支护力，可以看出随着开挖面水平间距比s/D逐渐增加，开挖面极限支护力呈现先降低后升高的趋势，在s/D = -0.5时最低，这是由于此时失稳破坏区域与既有隧道的重合区域最大，既有隧道对开挖面前方土体的“支撑”作用最显著，因此所需的极限支护力越低。对于同一工况，软黏土地层的极限支护力最大，这是由于软黏土强度低且破坏范围较大。

图 8 不同近接距离下开挖面极限支护力

Figure 8. Limit support pressures on tunnel face under different adjacent distances

下载: 全尺寸图片幻灯片

2. 三维旋转破坏模型

2.1 理论模型描述

基于数值模拟结果，构建近接既有隧道条件下的开挖面三维旋转破坏模型。对于近接工况，当s/D = -0.5时，开挖面破坏模式受既有隧道的影响最大，因而先依据此工况建立理论模型进行推导。影响开挖面破坏模式的主要因素为隧道净距，但若垂直净距过大，则不符合近接施工条件。因此，为体现近接施工特点，同时参考相关规范中对于近接施工的定义^[23]，本文假定隧道垂直净距不大于1倍盾构直径。

根据数值模拟结果，将开挖面破坏模式简化为对数螺旋体+筒仓体，如图 9所示。对于筒仓体部分，为体现沿开挖面前方和两侧扩展范围的差异性，假设其为半长轴、半短轴和高度分别设为a，b和h的椭圆柱；对数螺旋体由对数螺旋线旋转得到，对数螺旋线又称等角螺旋线，曲线上每个点的切向量与对数螺旋线极点至该点所成的向量均成一固定角度^[21]，满足

。

$r = {r_0} \cdot \exp \left[ {(\theta - {\theta _0})\tan \varphi } \right] 。$

(1)

图 9 改进的三维旋转破坏模型

Figure 9. Improved three-dimensional rotational failure mode

下载: 全尺寸图片幻灯片

目前已有一些学者基于对数螺旋体模型开展了盾构开挖面稳定性的相关研究^{[14, 17]}，但其截面被假定为圆形，这在多数情况下与实际不符，且无法与本模型中的上部椭圆体对应，因此本文采用截面为椭圆形的对数螺旋体作为破坏体模型。对数螺旋体具体生成方式为：极角为θ的射线OG与对数螺旋线AH、BE分别相交，以两条对数螺旋线之间线段长度的一半为椭圆半短轴b。假定同一地层中椭圆离心率保持不变，即椭圆长短轴之比不变，设a = kb。则

。

$a = k\frac{{{r_2} - {r_1}}}{2}, {\text{ }}b = \frac{{{r_2} - {r_1}}}{2} 。$

(2)

其中，

，

${r_1} = D\tan {\beta _B} \cdot \exp \left[ {(\theta - {\text{π }}/2)\tan \varphi } \right] \text{，}$

(3a)

。

${r_2} = D\sec {\beta _B} \cdot \exp \left[ {{\text{(}}\theta - {\beta _B}{\text{)}}\tan \varphi } \right] 。$

(3b)

令射线OG与对数螺旋线BE相交于G点，与开挖面中心线AB相交于F点，过F点作平行于x轴的直线F₁F₂与开挖面圆周分别相交于F₁、F₂点，截去椭圆F₁GF₂外侧部分得到椭圆缺形截面，该截面关于FG两侧对称。令极角θ由 ${\beta _B}$ 增加至π/2，则生成的所有截面叠加即形成对数螺旋体。图 10为近接既有隧道条件下盾构开挖面失稳变形区域的数值模拟结果与理论模型对比，可以看出，所提出的对数螺旋体+筒仓体模型在各个方向上均可较好地拟合实际位移包络面。

图 10 数值模拟结果与理论模型对比

Figure 10. Comparison between FDM simulation and theoretical model

下载: 全尺寸图片幻灯片

2.2 理论模型求解

首先，求解椭圆柱对下部对数螺旋体的土压力。根据图 9所示条分法，并假设土体服从Mohr–Coulomb破坏准则，建立垂直方向受力平衡方程：

，

${\text{π }}{a_0}{b_0}\gamma {\text{d}}z = {\text{π }}{a_0}{b_0}{\text{d}}{\sigma _{\text{v}}} +\\ \left[ {2{\text{π }}{b_0} + 4({a_0} - {b_0})} \right]{\text{(}}c + K\tan \varphi {\sigma _{\text{v}}}{\text{)d}}z \text{，}$

(4)

。

$\gamma = \sum\limits_{i = 1}^n {{\gamma _i}{h_i}} /({C_0} + {D_0}) 。$

(5)

式中：a₀，b₀分别为椭圆柱的半长轴和半短轴；γ为椭圆柱体上方考虑既有隧道空心影响的土体加权平均重度；γ_i，h_i分别为既有隧道上方各层土体的重度和层厚；c，φ分别为土体黏聚力和内摩擦角； $K = 1 - \sin \varphi$ 为土体侧压力系数； ${\sigma _{\text{v}}}$ 为竖向土压力。

椭圆柱体上表面压力为既有隧道底部处的土体自重应力，即 ${\sigma _{\text{v}}}$ (0) = γ(C₀+D₀)，其中，C₀，D₀分别为既有隧道埋深和直径。由式(4)及a = kb可得

。

${\sigma }_{\text{v}}(z)=\frac{\left[\text{π}k{b}_{0}\gamma -2(\text{π}+2k-2)c\right]\left(1-{\text{e}}^{-\frac{2\left(\text{π}+2k-2\right)K\mathrm{tan}\varphi }{\text{π}k{b}_{0}}z}\right)}{2(\text{π}+2k-2)K\mathrm{tan}\varphi }+\gamma ({C}_{0}+{D}_{0})\text{ }。$

(6a)

进而可得椭圆柱体作用在对数螺旋体顶面的竖向力 ${\sigma _{{\text{v1}}}}$ (其中柱体高度h取为隧道垂直净距d)：

。

${\sigma }_{\text{v1}}=\frac{\left[\text{π}k{b}_{0}\gamma -2(\text{π}+2k-2)c\right]\left(1-{\text{e}}^{-\frac{2\left(\text{π}+2k-2\right)K\mathrm{tan}\varphi }{\text{π}k{b}_{0}}d}\right)}{2(\text{π}+2k-2)K\mathrm{tan}\varphi }+\gamma ({C}_{0}+{D}_{0})\text{ }。$

(6b)

式中，当 $c \geqslant {\text{π }}k{b_0}\gamma /(2{\text{π }} + 4k - 4)$ 时， ${\sigma _{{\text{v1}}}}$ 取0，b₀可由式(2)令θ = π/2得到。

其次，建立对数螺旋体对原点O的力矩平衡方程：

${M_{\text{w}}} + {M_{\text{v}}} - {M_{\text{T}}} - {M_{\text{c}}} = 0 \text{，}$

(7)

式中： ${M_{\text{w}}}$ 为重力力矩； ${M_{\text{v}}}$ 为上表面压力 ${\sigma _{{\text{v1}}}}$ 产生的力矩； ${M_{\text{T}}}$ 为开挖面支护力力矩； ${M_{\text{c}}}$ 为土体抗力矩。

如所示，由于直接计算 ${M_{\text{w}}}$ 较为复杂，因此采用完整对数螺旋体AHBE重力力矩减去不规则体AHB重力力矩进行计算：

${M_{\text{w}}} = {M_{AHBE}} - {M_{AHB}} 。$

(8)

图 11 M_w计算图

Figure 11. Calculation diagram for M_w

下载: 全尺寸图片幻灯片

完整对数螺旋体AHBE微元体的重力力矩为

${\text{d}}{M_{AHBE}} = \gamma \cdot {\text{π }}ab\left[ {\frac{1}{2}({r_{\text{1}}} + {r_{\text{2}}}){\text{d}}\theta } \right]({r_2} - b)\sin \theta \text{，}$

(9)

对式(9)积分得

${M_{AHBE}} = \int_{{\beta _B}}^{\frac{{\text{π }}}{2}} {\gamma \cdot {\text{π }}ab\left[ {\frac{1}{2}({r_{\text{1}}} + {r_{\text{2}}})} \right]({r_2} - b)\sin \theta {\text{d}}\theta } 。$

(10)

不规则体AHB可以近似看作底面为正方形的楔形体，根据面积等效原则，其边长为 $\sqrt {\text{π }} D/2$ ，则

${M_{AHB}} = \gamma \cdot \frac{1}{2}{\left( {\frac{{\sqrt {\text{π }} D}}{2}} \right)^2}{h_{\text{H}}} \cdot \left( {D\tan {\beta _B} - \frac{{{h_{\text{H}}}}}{3}} \right) \\ {h_{\text{H}}} = D\tan {\beta _B}\left\{ {1 - \exp \left[ {({\beta _B} - {\text{π }}/2)\tan \varphi } \right]} \right\} \text{，}$

(11)

式中， ${h_{\text{H}}}$ 为三角形AHB中底边AB的高。

上表面压力 ${\sigma _{{\text{v1}}}}$ 产生的力矩 ${M_{\text{V}}}$ 为

${M_{\text{V}}} = {\text{π }}kb_0^2{\sigma _{{\text{v}}1}}(OE - {b_0}) \text{，}$

(12)

其中， $OE = \frac{D}{{\cos {\beta _B}}} \cdot \exp \left[ {\left( {\frac{{\text{π }}}{2} - {\beta _B}} \right)\tan \varphi } \right]$ 。

开挖面支护力力矩 ${M_{\text{T}}}$ 为

${M_{\text{T}}} = {\sigma _{\text{T}}}{\text{π }}({D^2}/4)(D/2) 。$

(13)

采用条分法计算土体抗力矩 ${M_{\text{c}}}$ ，如图 12所示。

图 12 条分单元体受力平衡

Figure 12. Force equilibrium of soil slices

下载: 全尺寸图片幻灯片

假设滑动面上土体强度同时得到发挥，且服从Mohr–Coulomb破坏准则：

${\tau _{\text{f}}} = c + \sigma \tan \varphi 。$

(14)

将条分单元体近似看作圆缺形，则根据第i个土条的竖直方向受力平衡，有

${\sigma _{{\text{v1}}}}{l_i} \cdot {P_1}{P_2} + {W_i} = {T_i}\sin ({\theta _i} - \varphi ) + {N_i}\cos {\text{(}}{\theta _i} - \varphi {\text{)}} \text{，}$

(15)

${T}_{i}=c{l}_{i}\stackrel{⌢}{{P}_{1}Q{P}_{2}}+{N}_{i}\mathrm{tan}\varphi 。$

(16)

其中： ${l_i} = {r_i}{\text{d}}\theta$ ； ${W_i} = \gamma \cdot \Delta {S_i}{l_i}\cos ({\theta _i} - \varphi )$ 为第i个土条的重力；T_i，N_i分别为第i个土条的滑动面摩擦力合力和支持力合力。

ΔS_i，P₁P₂和 $\stackrel{⌢}{{P}_{1}Q{P}_{2}}$ 可根据几何性质进行求解：

$\begin{array}{l}\Delta {S}_{i}=\frac{\mathrm{arccos}\left(\frac{{P}_{1}{P}_{2}{}^{2}-4P{Q}^{2}}{{P}_{1}{P}_{2}{}^{2}+4P{Q}^{2}}\right)\cdot \left(PQ+\frac{{P}_{1}{P}_{2}{}^{2}}{4PQ}\right)+{P}_{1}{P}_{2}\left(PQ-\frac{{P}_{1}{P}_{2}{}^{2}}{4PQ}\right)}{4}\text{，}\\ {P}_{1}{P}_{2}=2{a}_{0}\sqrt{1-{\left(\frac{{r}_{i}\mathrm{sin}\theta -D\mathrm{tan}{\beta }_{B}-{b}_{0}}{{b}_{0}}\right)}^{2}}\text{，}\\ PQ={r}_{i}\mathrm{cos}\theta \end{array} \\ \stackrel{⌢}{{P}_{1}Q{P}_{2}}=\mathrm{arccos}\left(\frac{{P}_{1}{P}_{2}{}^{2}-4P{Q}^{2}}{{P}_{1}{P}_{2}{}^{2}+4P{Q}^{2}}\right)\cdot \left(PQ+\frac{{P}_{1}{P}_{2}{}^{2}}{4PQ}\right) 。$

(17)

由式(15)~(17)可得

$\begin{array}{l}{N}_{i}=\frac{{\sigma }_{\text{v}1}{r}_{i}\mathrm{cos}({\theta }_{i}-\varphi ){P}_{1}{P}_{2}\text{d}\theta +{W}_{i}-c{r}_{i}\stackrel{⌢}{{P}_{1}Q{P}_{2}}\mathrm{sin}({\theta }_{i}-\varphi )\text{d}\theta }{\mathrm{tan}\varphi \mathrm{sin}({\theta }_{i}-\varphi )+\mathrm{cos}({\theta }_{i}-\varphi )}\text{ }\text{，}\\ {T}_{i}=c{l}_{i}+\frac{{\sigma }_{\text{v}1}{r}_{i}\mathrm{cos}({\theta }_{i}-\varphi ){P}_{1}{P}_{2}\text{d}\theta +{W}_{i}-c{r}_{i}\stackrel{⌢}{{P}_{1}Q{P}_{2}}\mathrm{sin}({\theta }_{i}-\varphi )\text{d}\theta }{\mathrm{tan}\varphi \mathrm{sin}({\theta }_{i}-\varphi )+\mathrm{cos}({\theta }_{i}-\varphi )}\mathrm{tan}\varphi 。\end{array}$

(18)

滑动面上土体抗力的微力矩为

$\left. \begin{array}{l} {\rm{d}}{M_{{\rm{cT}}}} = {T_i}{r_i}{\rm{cos}}\varphi ,\\ {M_{{\rm{cN}}}} = {N_i}{r_i}{\rm{sin}}\varphi 。 \end{array} \right\}$

(19)

则滑动面上土体抗力的总力矩为

${M_{\text{c}}} = {M_{{\text{cT}}}} + {M_{{\text{cN}}}} = \int_{{\beta _B}}^{\frac{{\text{π }}}{2}} {({T_i}{r_i}\cos \varphi + {N_i}{r_i}\sin \varphi )}\\ ={\displaystyle {\int }_{{\beta }_{B}}^{\frac{\text{π}}{2}}\left(c{r}_{i}^{2}\stackrel{⌢}{{P}_{1}Q{P}_{2}}\mathrm{cos}\varphi +\frac{2{N}_{i}{r}_{i}\mathrm{sin}\varphi }{\text{d}\theta }\right)}\text{d}\theta 。$

(20)

最后，联立式(7)，(8)，(10)~(18)和(20)，可得

${\sigma _{\text{T}}} = c{N_c} + \gamma D{N_\gamma } 。$

(21)

式中： ${\sigma _{\text{T}}}$ 为开挖面支护力； ${N_c}$ ， ${N_\gamma }$ 分别为土体黏聚力和重度对开挖面支护力的影响系数：

$\begin{array}{l}{N}_{c}=\frac{8k{b}_{0}^{2}(OE-{b}_{0})\left[\mathrm{exp}\left(-\frac{2(\text{π}+2k-2)K\mathrm{tan}\varphi }{\text{π}k{b}_{0}}d\right)-1\right]}{{D}^{3}K\mathrm{tan}\varphi }-\\ \frac{8}{\text{π}{D}^{3}}{\displaystyle {\int }_{{\beta }_{B}}^{\frac{\text{π}}{2}}\{{r}_{2}^{2}\stackrel{⌢}{{P}_{1}Q{P}_{2}}\mathrm{cos}\varphi +\frac{2\mathrm{sin}\varphi \cdot {r}_{\text{2}}{}^{2}\mathrm{cos}(\theta -\varphi ){P}_{1}{P}_{2}}{K\mathrm{tan}\varphi \left[\mathrm{tan}\varphi \mathrm{sin}(\theta -\varphi )+\mathrm{cos}(\theta -\varphi )\right]}\cdot }\\ \left[\mathrm{exp}\left(-\frac{2(\text{π}+2k-2)K\mathrm{tan}\varphi }{\text{π}k{b}_{0}}d\right)-1\right]-\end{array} \\ \frac{2\mathrm{sin}\varphi \cdot {r}_{\text{2}}{}^{2}\stackrel{⌢}{{P}_{1}Q{P}_{2}}\mathrm{sin}\text{(}\theta -\varphi \text{)}}{\mathrm{tan}\varphi \mathrm{sin}\text{(}\theta -\varphi \text{)}+\mathrm{cos}\text{(}\theta -\varphi \text{)}}\}\text{d}\theta \text{，}$

(22)

${N_\gamma } = \frac{{4\int_{{\beta _B}}^{\frac{{\text{π }}}{2}} {ab({r_{\text{1}}} + {r_{\text{2}}})({r_2} - b)\sin \theta {\text{d}}\theta } }}{{{D^4}}} - \frac{{{h_{\text{H}}} \cdot \left( {D\tan {\beta _{\text{B}}} - \frac{{{h_{\text{H}}}}}{3}} \right)}}{{{D^2}}} +\\ \begin{array}{l} \frac{{8k{b_0}^2(OE - {b_0})({C_0} + {D_0})}}{{{D^4}}} + \hfill \\ \frac{{4{\text{π }}{k^2}{b_0}^3(OE - {b_0})\left[ {1 - \exp \left( { - \frac{{2({\text{π }} + 2k - 2)K\tan \varphi }}{{{\text{π }}k{b_0}}}d} \right)} \right]}}{{{D^4}({\text{π }} + 2k - 2)K\tan \varphi }} - \hfill \\ \frac{8}{{{\text{π }}{D^4}}}\int_{{\beta _B}}^{\frac{{\text{π }}}{2}} {\left\{ \begin{array}{l} \hfill \\ \hfill \\ \hfill \\ \end{array} \right.} \frac{{2\sin \varphi \cdot {r_{\text{2}}}^2\cos (\theta - \varphi )\Delta S}}{{\tan \varphi \sin (\theta - \varphi ) + \cos (\theta - \varphi )}} + \hfill \\ \frac{{2\sin \varphi \cdot {r_{\text{2}}}^2\cos (\theta - \varphi ){P_1}{P_2}}}{{\tan \varphi \sin (\theta - \varphi ) + \cos (\theta - \varphi )}} \cdot \hfill \\ \end{array} \\ \left. {\left( {\frac{{{\text{π }}k{b_0}\left[ {1 - \exp \left( { - \frac{{2({\text{π }} + 2k - 2)K\tan \varphi }}{{{\text{π }}k{b_0}}}d} \right)} \right]}}{{2({\text{π }} + 2k - 2)K\tan \varphi }} + C + {D_0}} \right)} \right\} \text{d}\theta 。$

(23)

当s/D为其他值时，可采用线性插值对开挖面极限支护力进行求解，具体为

${\sigma }_{\text{T}}^{*}=\{\begin{array}{l}(0.5-s/D){\sigma }_{\text{Tm}}+(s/D+0.5){\sigma }_{\text{T0}}\text{ (}-0.5 < s/D < 0.5)\\ 2(s/D+1){\sigma }_{\text{Tm}}-2(s/D+0.5){\sigma }_{\text{T0}}\text{ (}-1 < s/D < -0.5)\\ {\sigma }_{\text{T0}}\text{ (}s/D\le -1或s/D\ge 0.5)\end{array}。$

(24)

式中： $\sigma _{\text{T}}^*$ 为相应s/D处的极限支护力； ${\sigma _{{\text{Tm}}}}$ 为近接条件下s/D =-0.5时的极限支护力； ${\sigma _{{\text{T0}}}}$ 为相同土性和埋深的单隧道开挖面极限支护力，可令k = 1进行计算。

3. 理论模型验证及参数分析

3.1 理论模型验证

理论模型中自定义参数为两个：OB与垂直面的夹角 ${\beta _B}$ 及椭圆长短轴之比k。根据村山公式原理^[6]，对数螺旋线在盾构底面与水平方向夹角为 ${\text{π }}/4 + \varphi /2$ ，则 ${\beta _B} = {\text{π /}}4 - \varphi /2$ 。近接条件下k值主要受地层条件和既有隧道影响，可由数值模拟结果得到(松砂：k = 1.8；密砂：k = 1.3；软黏：k = 1.9；硬黏：k = 1.6)。由于式(22)和(23)较为复杂，可代入已知参数至Matlab中进行求解。当隧道埋深较大(通常C/D > 3.0)，则盾构施工引起的土拱效应、渐进性破坏特性与常规埋深隧道完全不同^[13]，因此本文暂不考虑大埋深情况。

图 13为本文理论分析和数值模拟结果与其他学者的结果对比。单隧道条件下(图 13(a))的理论分析结果与已有试验结果和本文数值模拟结果均符合较好。图 13(b)为松砂地层近接条件下(s/D = 0)的结果对比，可以看出，本文基于三维旋转破坏模型得到的极限支护力与数值模拟结果的偏差在10%以内，表明计算结果较好，且相比于数值模拟繁琐的建模过程可以大大提高不同工况下的计算分析效率。

图 13 开挖面极限支护力结果对比

Figure 13. Comparison of limit support pressure of tunnel face

下载: 全尺寸图片幻灯片

3.2 理论模型参数分析

在本文理论模型中，参数k体现了近接工况下破坏体横截面的椭圆程度，这一现象是由于土体变形向上方的扩展被既有隧道限制而只能沿水平方向扩展，并且沿既有隧道方向的土体存在与既有隧道“剥离”的趋势，更加容易产生变形，因而与垂直方向的土体变形扩展产生差异。土体c，φ值越高，变形越不易扩展，导致土体更趋于各向同性，因此两个方向的变形范围差异(k值)随土性抗剪强度增加而降低。

如图 14所示，开挖面极限支护力随黏聚力的增加而线性降低，随土体内摩擦角φ增加呈现出先增加后降低的趋势，随埋深增加近似线性增加。极限支护力影响系数N_c、N_γ与土体内摩擦角的关系如图 15所示。由于土体黏聚力起抗力作用，因而N_c为负数。N_c随土体φ值增加而先增加后降低，随近接距离增加而降低。N_γ随φ值增加先基本不变而后逐渐降低，随近接距离增加而增加。此外，当土体内摩擦角较大(φ > 30°)时，近接距离对N_c的影响很小，这是由于内摩擦角越大，破坏范围越小，邻近的既有隧道的影响越低。

图 14 土性参数对极限支护力的影响

Figure 14. Influences of soil parameters on limit support pressure

下载: 全尺寸图片幻灯片

图 15 极限支护力系数与土体内摩擦角关系

Figure 15. Relationship between coefficient of limit support pressure and internal friction angle of soils

下载: 全尺寸图片幻灯片

4. 结论

本文揭示了近接既有隧道条件下的盾构开挖面破坏模式，构建了对数螺旋体+椭圆柱的开挖面失稳三维破坏模型并对主要参数进行了分析，得到以下3点结论。

(1) 近接既有隧道条件下盾构开挖面失稳三维变形区呈现沿隧道纵向和横向不对称的特点。当盾构开挖面与既有隧道轴线距离(s/D)在-1~0.5时，破坏模式呈现螺旋体+筒仓体形态。隧道垂直净距和土性对近接条件下开挖面破坏模式影响不大，但会改变三维方向上破坏扩展范围。

(2) 改进的对数螺旋体+椭圆柱体模型可以反映开挖面失稳扩展沿隧道纵向和横向的非均匀性，因此更加符合近接既有隧道条件下的实际开挖面破坏模式，且可以通过设置参数k使其适用于单隧道条件。基于所提出的破坏模型，推导了可以反映土体黏聚力和重度影响系数的开挖面极限支护力计算公式。

(3) 理论模型中参数k可体现开挖面失稳变形沿隧道纵向和横向扩展的非均匀性。土体强度越高(黏聚力c和内摩擦角φ越大)、隧道埋深越浅、近接距离越小，则开挖面极限支护力越低。开挖面极限支护力影响系数N_c、N_γ随土体内摩擦角φ增加呈现先增加后降低的变化趋势。内摩擦角临界值约为11°~13°。

图 1 MPM计算流程

Figure 1. Simulation cycle of MPM

下载: 全尺寸图片幻灯片

图 2 CPDI的初始和更新的粒子域

Figure 2. Initial and updated particle domains in CPDI method

下载: 全尺寸图片幻灯片

图 3 粒子域速度场

Figure 3. Particle domain velocity field

下载: 全尺寸图片幻灯片

图 4 柱的几何模型和边界条件

Figure 4. Geometry and boundary conditions of column

下载: 全尺寸图片幻灯片

图 5 柱的垂直应力数值解和解析解(h_{x, y}=1/64 m, np=4)

Figure 5. Numerical and analytical solutions for vertical stresses of column (h_{x, y}=1/64 m, np=4)

下载: 全尺寸图片幻灯片

图 6 不同粒子密度下垂直应力解析解和数值解的误差(h_{x, y}=1.0 m)

Figure 6. Errors of analytical and numerical solutions for vertical stresses of column at different particle densities (h_{x, y}=1.0 m)

下载: 全尺寸图片幻灯片

图 7 不同网格尺寸下垂直应力解析解和数值解的误差(np=4)

Figure 7. Errors of analytical and numerical solutions for vertical stresses of column at different grid cell sizes (np=4)

下载: 全尺寸图片幻灯片

图 8 砂柱坍塌试验示意图、模型的网格和初始粒子设置

Figure 8. Experimental schematics, grid and initial particle setup of sand column

下载: 全尺寸图片幻灯片

图 9 砂柱坍塌改进CPDI模拟结果与试验结果^[30-31]对比

Figure 9. Comparison between improved CPDI simulation and experimental results^[30-31] for sand column collapse

下载: 全尺寸图片幻灯片

图 10 砂柱坍塌自由面对比

Figure 10. Comparison of free surface of collapsing sand column

下载: 全尺寸图片幻灯片

图 11 隧道几何模型和边界条件

Figure 11. Geometry and boundary conditions of tunnel

下载: 全尺寸图片幻灯片

图 12 离心机试验和改进CPDI计算的隧道坍塌

Figure 12. Centrifuge test and improved CPDI calculation for tunnel collapse

下载: 全尺寸图片幻灯片

图 13 地面塌陷

Figure 13. Ground collapse

下载: 全尺寸图片幻灯片

图 14 静沙区间地质剖面图

Figure 14. Geological profile of Jinggang Road Station-Shazikou Station section

下载: 全尺寸图片幻灯片

图 15 隧道几何模型和边界条件

Figure 15. Geometry and boundary conditions of tunnel

下载: 全尺寸图片幻灯片

图 16 隧道坍塌过程中的地面沉降

Figure 16. Ground subsidence during tunnel collapse

下载: 全尺寸图片幻灯片

图 17 隧道坍塌过程中的剪应变

Figure 17. Shear strains during tunnel collapse

下载: 全尺寸图片幻灯片

图 18 隧道坍塌后地面沉降曲线

Figure 18. Curve of ground subsidence after tunnel collapse

下载: 全尺寸图片幻灯片

图 19 监测点竖向位移和水平位移时程曲线

Figure 19. Time-history curves of vertical and horizontal displacements at monitoring points

下载: 全尺寸图片幻灯片

表 1 地层物理力学参数

Table 1 Physical and mechanical parameters of strata

地层	弹性模量E/MPa	泊松比 $\nu$	黏聚力c/kPa	内摩擦角φ/(°)	重度γ/(kN·m^-3)
杂填土	8.0	0.20	0	15	17.5
中粗砂	6.07	0.33	13.9	12.5	18.5
粉质黏土	5.671	0.30	8.2	12	19.7
强风化凝灰岩	20	0.30	3.0	30	22.5
中风化凝灰岩	50	0.25	3000	45	26.0
微风化凝灰岩	5000	0.22	11500	55	26.7

下载: 导出CSV

参考文献(35)

[1]	张成平, 张顶立, 王梦恕, 等. 城市隧道施工诱发的地面塌陷灾变机制及其控制[J]. 岩土力学, 2010, 31(增刊1): 303-309. https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX2010S1050.htm ZHANG Chengping, ZHANG Dingli, WANG Mengshu, et al. Catastrophe mechanism and control technology of ground collapse induced by urban tunneling[J]. Rock and Soil Mechanics, 2010, 31(S1): 303-309. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX2010S1050.htm
[2]	SULSKY D, CHEN Z, SCHREYER H L. A particle method for history-dependent materials[J]. Computer Methods in Applied Mechanics and Engineering, 1994, 18(1/2): 179-196.
[3]	SULSKY D, ZHOU S J, SCHREYER H L. Application of a particle-in-cell method to solid mechanics[J]. Computer Physics Communications, 1995, 87(1/2): 236-252.
[4]	HARLOW F H. The Particle-in-cell Method for Numerical Solution of Problems in Fluid Dynamics[R]. Los Alamos: Los Alamos National Lab, 1962.
[5]	BRACKBILL J U, KOTHE D B, RUPPEL H M. Flip: a low-dissipation, particle-in-cell method for fluid flow[J]. Computer Physics Communications, 1988, 48(1): 25-38. doi: 10.1016/0010-4655(88)90020-3
[6]	BRACKBILL J U, RUPPEL H M. FLIP: a method for adaptively zoned, particle-in-cell calculations of fluid flows in two dimensions[J]. Journal of Computational Physics, 1986, 65(2): 314-343. doi: 10.1016/0021-9991(86)90211-1
[7]	BANDARA S, SOGA K. Coupling of soil deformation and pore fluid flow using material point method[J]. Computers and Geotechnics, 2015, 63: 199-214. doi: 10.1016/j.compgeo.2014.09.009
[8]	CUOMO S, PERNA A D, MARTINELLI M. Material point method (MPM) hydro-mechanical modelling of flows impacting rigid walls[J]. Canadian Geotechnical Journal, 2021, 58: 1730-1743. doi: 10.1139/cgj-2020-0344
[9]	王兆南, 王刚. 饱和孔隙介质的耦合物质点-特征有限元方法[J]. 岩土工程学报, 2023, 45(5): 1094-1102. doi: 10.11779/CJGE20220332 WANG Zhaonan, WANG Gang. Coupled material point method and characteristic finite element method for saturated porous media[J]. Chinese Journal of Geotechnical Engineering, 2023, 45(5): 1094-1102. (in Chinese) doi: 10.11779/CJGE20220332
[10]	孙玉进, 宋二祥. 大位移滑坡形态的物质点法模拟[J]. 岩土工程学报, 2015, 37(7): 1218-1225. doi: 10.11779/CJGE201507007 SUN Yujin, SONG Erxiang. Simulation of large-displacement landslide by material point method[J]. Chinese Journal of Geotechnical Engineering, 2015, 37(7): 1218-1225. (in Chinese) doi: 10.11779/CJGE201507007
[11]	钟祖良, 贺凯源, 宋宜祥, 等. 基于仿射速度矩阵改进物质点法的大位移滑坡研究[J]. 岩土工程学报, 2022, 44(9): 1626-1634. doi: 10.11779/CJGE202209007 ZHONG Zuliang, HE Kaiyuan, SONG Yixiang, et al. Large-displacement landslides based on affine velocity matrix-improved material point method[J]. Chinese Journal of Geotechnical Engineering, 2022, 44(9): 1626-1634. (in Chinese) doi: 10.11779/CJGE202209007
[12]	CUOMO S, PERNA A D, MARTINELLI M. Modelling the spatio-temporal evolution of a rainfall-induced retrogressive landslide in an unsaturated slope[J]. Engineering Geology, 2021, 294: 106371. doi: 10.1016/j.enggeo.2021.106371
[13]	CORTIS M, COOMBS W, AUGARDE C, BROWN M, BRENNAN A, ROBINSON S. Imposition of essential boundary conditions in the material point method[J]. International Journal for Numerical Methods in Engineering, 2018, 113: 130-152. doi: 10.1002/nme.5606
[14]	BARDENHAGEN S G, KOBER E M. The Generalized interpolation material point method[J]. Computer Modeling in Engineering and Sciences, 2004, 5(6): 477-495.
[15]	SADEGHIRAD A, BRANNON R M, BURGHARDT J. A convected particle domain interpolation technique to extend applicability of the material point method for problems involving massive deformations[J]. International Journal for Numerical Methods in Engineering, 2011, 86(12): 1435-1456. doi: 10.1002/nme.3110
[16]	SADEGHIRAD A, BRANNON R M, GUILKEY J E. Second-order convected particle domain interpolation (CPDI2) with enrichment for weak discontinuities at material interfaces[J]. International Journal for Numerical Methods in Engineering, 2013, 95(11): 928-952. doi: 10.1002/nme.4526
[17]	PRUIJN N S. The improvement of the Material Point Method by Increasing Efficiency and Accuracy[D]. Delft: Delft University of Technology, 2016.
[18]	STEFFEN M, KIRBY R M, BERZINS M. Analysis and reduction of quadrature errors in the material point method (MPM)[J]. International Journal for Numerical Methods in Engineering, 2008, 76(6): 922-948. doi: 10.1002/nme.2360
[19]	STEFFEN M, WALLSTEDT P C, GUILKEY J E, et al. Examination and analysis of implementation choices within the material point method (MPM)[J]. Computer Modeling in Engineering and Sciences, 2008, 31(2): 107-127.
[20]	HU Y M, FANG Y, GE Z H, et al. A moving least squares material point method with displacement discontinuity and two-way rigid body coupling[J]. ACM Transactions on Graphics, 2018, 37(4): 1-14.
[21]	SONG J U, KIM H G. An improved material point method using moving least square shape functions[J]. Computational Particle Mechanics, 2021, 8(4): 751-766. doi: 10.1007/s40571-020-00368-9
[22]	WANG J F, SUN F X, CHENG Y M. An improved interpolating element-free Galerkin method with a nonsingular weight function for two-dimensional potential problems[J]. Chinese Physics B, 2012, 21(9): 090204. doi: 10.1088/1674-1056/21/9/090204
[23]	WANG J F, WANG J F, SUN F, et al. An interpolating boundary element-free method with nonsingular weight function for two-dimensional potential problems[J]. International Journal of Computational Methods, 2013, 10(6): 1350043. doi: 10.1142/S0219876213500436
[24]	WANG Q, ZHOU W, FENG Y T, et al. An adaptive orthogonal improved interpolating moving least-square method and a new boundary element-free method[J]. Applied Mathematics and Computation, 2019, 353: 347-370. doi: 10.1016/j.amc.2019.02.013
[25]	MIRZAEI D. Analysis of moving least squares approximation revisited[J]. Journal of Computational and Applied Mathematics, 2015, 282: 237-250. doi: 10.1016/j.cam.2015.01.007
[26]	TRAN Q A, SOŁOWSKI W, BERZINS M, et al. A convected particle least square interpolation material point method[J]. International Journal for Numerical Methods in Engineering, 2019, 121(6): 1068-1100.
[27]	WALLSTEDT P C, GUILKEY J E. An evaluation of explicit time integration schemes for use with the generalized interpolation material point method[J]. Journal of Computational Physics, 2008, 227(22): 9628-9642. doi: 10.1016/j.jcp.2008.07.019
[28]	WYSER E, ALKHIMENKOV Y, JABOYEDOFF M, et al. A fast and efficient MATLAB-based MPM solver: fMPMM-solver v1.1[J]. Geoscientific Model Development, 2020, 13(12): 6265-6284. doi: 10.5194/gmd-13-6265-2020
[29]	COOMBS W M, AUGARDE C. AMPLE: a material point learning environment[J]. Advances in Engineering Software, 2020, 139: 102748. doi: 10.1016/j.advengsoft.2019.102748
[30]	LUBE G, HUPPERT H E, SPARKS R S J, et al. Collapses of two-dimensional granular columns[J]. Physical Review E Statistical, Nonlinear, Biological and Soft Matter Physics, 2005, 72(4): 041301. doi: 10.1103/PhysRevE.72.041301
[31]	LUBE G, HUPPERT H E, SPARKS R S J, et al. Static and flowing regions in granular collapses down channels: insights from a sedimenting shallow water model[J]. Physics of Fluids, 2007, 19(10): 106601. doi: 10.1063/1.2773738
[32]	KAMATA H, MASHIMO H. Centrifuge model test of tunnel face reinforcement by bolting[J]. Tunnelling and Underground Space Technology, 2003, 18(2/3): 205-212.
[33]	CHENG X S, ZHENG G, SOGA K, et al. Post-failure behavior of tunnel heading collapse by MPM simulation[J]. Science China Technological Sciences, 2015, 58(12): 2139-2153. doi: 10.1007/s11431-015-5874-4
[34]	ZHANG Y J, ZHANG W G, XIA H S, et al. Case study and risk assessment of water inrush disaster in Qingdao Metro Line 4[J]. Applied Sciences, 2023, 13: 3384. doi: 10.3390/app13063384
[35]	YAN F Y, QIU W G, SUN K G, et al. Investigation of a large ground collapse, water inrush and mud outburst, and countermeasures during subway excavation in Qingdao: a case study[J]. Tunnelling and Underground Space Technology, 2021, 117: 104127. doi: 10.1016/j.tust.2021.104127