0.   引言

随着交通强国战略的深入推进和城市的迅速发展，越来越多的城市开始修建地铁缓解交通压力。但复杂的地质环境给地下工程的建造带来了巨大的挑战，地铁隧道施工过程中遇到复杂软弱围岩时，由于岩石强度低、自身的稳定性较差，隧道掌子面开挖时围岩变形表现为一定的流变特征，开挖会导致其产生大变形，甚至出现地面塌陷和建（构）筑物的破坏倒塌等事故，给城市的安全运行带来极大的威胁^[1]。

目前，国内外针对隧道掌子面的稳定性问题的研究手段主要有理论分析、模型试验和数值模拟，并取得了较多的成果。随着计算机硬件的发展，数值计算手段成为岩土工程领域的首选方案。隧道开挖掌子面的稳定性和变形行为的数值模拟方面主要为有限元法（FEM）和有限差分法（FDM）。然而隧道掌子面破坏过程涉及到大变形及破坏后行为，基于网格的数值方法由于网格畸变无法很好地模拟大变形行为。

MPM方法是一种结合拉格朗日和欧拉描述的无网格方法^[2-3]，由粒子单元法（PIC）^[4]和流体隐式粒子法（FLIP）^[5-6]发展而来。MPM结合了拉格朗日和欧拉方法的优点，避免了网格畸变问题，同时可以模拟材料的历史行为^[2-3]。MPM是处理大变形问题的有效方法，已成功应用于流固耦合问题^[7-9]、边坡破坏等工程领域^[10-12]。MPM可以消除其他无网格方法的一些缺点，但当粒子穿越网格边界时，MPM会出现网格穿越误差并导致数值噪声，这是由于MPM中粒子和网格之间插值和映射的网格基函数缺乏平滑性，同时MPM难以处理任意边界条件问题^[13]。

为克服传统MPM存在的网格穿越误差，Bardenhagen等^[14]通过引入粒子特征函数，并采用Petrov-Galerkin离散化方法建立了广义插值材料点法（GIMP），GIMP减少了单元穿越误差，但由于粒子域的大变形，GIMP可能会出现计算不稳定问题。为了消除GIMP的计算不稳定性，Sadeghirad等^[15]提出了对流粒子域插值物质点法（CPDI），CPDI将GIMP中的矩形粒子域转化为平行四边形。随后，提出了二阶对流粒子域插值物质点法（CPDI2）^[16]，消除了粒子域之间的重叠或间隙。采用传统的有限元形函数求解应力时存在计算误差^[17]，通过提高插值函数的阶数，可以有效提高MPM的计算精度，Steffen等^[18-19]提出了采用二次和三次B样条函数以减少单元交叉误差。

研究表明移动最小二乘函数（MLS）可以在MPM大变形计算中实现更高的精度^[20-21]。然而，MLS方法需要求解逆矩阵，这既耗时同时在遇到奇异或病态矩阵时会导致计算终止。为了避免传统MLS方法的不稳定性，需对其进行修正，增强其鲁棒性。

本文在改进插值移动最小二乘（IIMLS）基础上，引入了加权正交基函数，并忽略了新对角矩阵中的零元素或极小元素的贡献，构造了自适应正交改进插值移动最小二乘（AOIIMLS）形函数。在CPDI计算框架中，采用AOIIMLS改进CPDI。AOIIMLS避免了求解逆矩阵，同时AOIIMLS满足Kronecker$\delta$函数性质，允许直接施加本质边界条件。AOIIMLS形函数用于网格和粒子之间的插值和映射以减少CPDI计算中的单元穿越误差，提高了计算精度。通过一维柱在自重下的压缩、砂柱坍塌和隧道坍塌离心机试验算例验证了改进CPDI方法的准确性和适用性。最后，采用改进CPDI模拟了青岛地铁4号线静沙区间隧道地面坍塌破坏的全过程。

1.   物质点法基本理论

MPM采用粒子和背景网格双重离散，粒子放置在背景网格中，如图 1（a）所示。粒子（也可称为材料点或质点）携带与历史相关的所有材料属性，如速度、应力、应变等，背景网格用于求解运动方程。粒子和网格节点之间的信息映射通过有限元形函数建立。本文中网格节点用大写下标I表示，粒子用小写下标p表示。

图  1  MPM计算流程

Figure  1.  Simulation cycle of MPM

下载: 全尺寸图片 幻灯片

1.1   控制方程

MPM采用更新的拉格朗日式建立控制方程，质量和动量守恒方程如下：

式中：$\rho$为密度；$\boldsymbol{v}$为速度矢量；$\boldsymbol{\sigma }$为Cauchy应力张量；$\boldsymbol{b}$为体力。

边界条件为

式中：${\mathit{\Gamma } _\tau }$和${\mathit{\Gamma } _u}$分别为给定的表面力边界和位移边界；n为边界${\mathit{\Gamma } _\tau }$的外法向单位矢量；$\bar{\boldsymbol{\tau} }$为边界${\mathit{\Gamma } _\tau }$上的表面力；$\bar{\boldsymbol{v}}$为位移边界${\mathit{\Gamma } _u}$的速度。

MPM以弱形式和连续介质力学为基础，式（2）在问题域Ω中可以采用其改进拉格朗日弱形式求解：

式中：a为加速度矢量；$\delta \boldsymbol{u}$为虚位移矢量；x为位置矢量。

1.2   粒子离散及变量更新

在MPM中，连续体Ω由离散的拉格朗日粒子表示，如图 1（a）所示。问题域的密度为

式中：n_p为粒子数；m_p为粒子p的质量；$\delta$为Dirac delta函数；x_p为粒子p的当前位置矢量。

在时间步t，MPM计算流程的第一步为将粒子的质量和速度映射到网格中（如图 1（a）所示），以计算网格节点的质量$m_I^t$和速度$\boldsymbol{v}_I^t$：

式中：$\boldsymbol{v}_p^t$为粒子p的速度，${N_I}$为网格节点I的有限元形函数。

如图 1（b）所示，计算背景网格节点力，施加边界条件，在背景网格上显式求解运动方程，得到网格节点的速度$\boldsymbol{v}_I^{t + \Delta t}$：

式中：$\boldsymbol{f}_I^{e,t}$和$\boldsymbol{f}_I^{i,t}$分别为节点I的外力和内力；$\boldsymbol{b}_p^t$和$\bar{\boldsymbol{\tau} } _p^t$分别为施加在粒子p上的体力和面力；h为虚拟的边界层厚度；$\boldsymbol{\sigma } _p^t$为粒子p的Cauchy应力张量；${\rho _p}$为粒子p的密度；$\Delta t$为时间步长。

将网格的信息映射回粒子，如图 1（d）所示。更新粒子的速度$\boldsymbol{v}_p^{t + \Delta t}$和位置$\boldsymbol{x}_p^{t + \Delta t}$，并丢弃已经变形的网格，如图 1（c）所示：

通常，在MPM计算中，新的背景网格只是简单地转移回旧网格的位置，所以网格生成过程非常简单。在MPM计算的整个过程中，粒子与背景网格固联，不需要处理欧拉方法中的对流项，因此便于施加边界条件。此外，MPM在每个时间步采用新的背景网格，有效地避免了网格畸变问题，因此MPM是求解大变形问题的一种有效的数值计算方法。

2.   改进对流粒子域插值物质点法

2.1   自适应正交改进插值移动最小二乘形函数

在IIMLS^[22-23]中，对于给定的点$\boldsymbol{x} = \left[ {x,y,z} \right] \in \mathit{\Omega }$，定义新的基函数${q_j}(\boldsymbol{x})$：

式中：${\boldsymbol{P}^{\text{T}}}(\boldsymbol{x}) = [{p_1}(\boldsymbol{x}),{p_2}(\boldsymbol{x}), \cdots ,{p_m}(\boldsymbol{x})]$为单项式基函数，${p_1}(\boldsymbol{x}) \equiv 1$；m为基函数的项数。在二维中，对于线性基函数${\boldsymbol{P}^{\text{T}}}(\boldsymbol{x}) = [1,x,y]$，对于二次基函数${\boldsymbol{P}^{\text{T}}}(\boldsymbol{x}) =$$[1,x,y,{x^2},xy,{y^2}]$，本文算例中均采用线性基函数。$\sum\nolimits_{i = 1}^n {s(\boldsymbol{x},{\boldsymbol{x}_i})} = 1,{\text{ }}\forall x \in \mathit{\Omega }$，${\boldsymbol{x}_i}(i = 1,2, \cdots ,n){\text{ }}$是以x为中心的影响域内的节点，n为节点数量。

式中，如果n=1，则$\zeta = 1$。

IIMLS中的矩阵可能是奇异或病态的，可以通过采用加权正交基函数克服。基于Gram-Schmidt正交化，正交化基函数可以定义为

式中：${({q_j},{\bar q_k})_{\boldsymbol{x}}} = \sum\nolimits_{i = 1}^n {w(\boldsymbol{x} - {\boldsymbol{x}_i})} {q_j}({\boldsymbol{x}_i}){\bar q_k}({\boldsymbol{x}_i})$，${({\bar q_k},{\bar q_k})_{\boldsymbol{x}}} =$$\sum\nolimits_{i = 1}^n {w(\boldsymbol{x} - {\boldsymbol{x}_i})} {\bar q_k}({\boldsymbol{x}_i}){\bar q_k}({\boldsymbol{x}_i})$。$w(\boldsymbol{x} - {\boldsymbol{x}_i})$为具有紧支特性的加权函数，本文采用三次样条加权函数：

式中：$r = {{{d_i}} \mathord{\left/ {\vphantom {{{d_i}} {{r_i}}}} \right. } {{r_i}}}$，${d_i} = \left\| {\boldsymbol{x} - {\boldsymbol{x}_i}} \right\|$，${r_i}$为影响域的半径。

对于给定的函数$u(\boldsymbol{x})$，定义一个新的函数$\tilde u(\boldsymbol{x})$和局部插值逼近函数${\tilde u_l}(\boldsymbol{x})$：

式中：${\bar a_j}(\boldsymbol{x})$为系数向量。${\bar L_2}$范数误差可表示为

式中：$\boldsymbol{S}(\boldsymbol{x})$为一个$n \times n$的矩阵，其中每一行为

通过${{\partial {{\bar L}_2}} \mathord{\left/ {\vphantom {{\partial {{\bar L}_2}} \partial }} \right. } \partial }\bar{\boldsymbol{a}} = \mathbf{0}$，可得到$\bar{\boldsymbol{a}}(\boldsymbol{x})$：

式中：

因此由式（18），（19），（25）可以得到全局插值逼近函数u(x)：

式中：${\bar{\boldsymbol{A}}^{ - 1}}(\boldsymbol{x})$为对角矩阵；$\boldsymbol{\varPhi} (\boldsymbol{x})$为正交改进插值移动最小二乘（OIIMLS）近似的形函数。其中，

式（31）中的每项可表示为

OIIMLS的形函数${\mathit{\Phi} _i}\left( x \right)$满足Kronecker δ函数的性质和再现性质：

OIIMLS不需要直接求解矩阵$\bar{\boldsymbol{A}}(\boldsymbol{x})$的逆，但是无法改变矩阵$\bar{\boldsymbol{A}}(\boldsymbol{x})$的奇异或者病态的性质。本文通过忽略矩阵$\bar{\boldsymbol{A}}(\boldsymbol{x})$中非常小的对角元素，自适应地解决这个问题^[24]。如果${({\bar q_k},{\bar q_k})_{\boldsymbol{x}}} = 0$或${({\bar q_k},{\bar q_k})_{\boldsymbol{x}}} < {\varepsilon _0}$，则忽略式（32）中的${\bar q_k}(\boldsymbol{x}){\bar q_k}({\boldsymbol{x}_l})/{({\bar q_k},{\bar q_k})_{\boldsymbol{x}}}$项。${\varepsilon _0}$为规定的非常小的公差，本文取${\varepsilon _0} = {10^{ - 20}}$。因此，可以自适应地避免了奇异或病态矩阵。

通过忽略非常小的对角线元素改变矩阵的奇异或者病态性质可能会导致精度的降低，通过引入Mirzaei^[25]提出的位移和缩放多项式基函数，可以提高自适应正交改进插值移动最小二乘法的精度和稳定性。在二维中，对于线性基${\boldsymbol{P}^{\text{T}}}(\boldsymbol{x}) = \left[ {1,{\text{ }}{x_{zh}},{\text{ }}{y_{zh}}} \right]$，对于二次基${\boldsymbol{P}^{\text{T}}}(\boldsymbol{x}) = \left[ {1,{\text{ }}{x_{zh}},{\text{ }}{y_{zh}},{\text{ }}x_{zh}^2,{\text{ }}{x_{zh}}{y_{zh}},{\text{ }}y_{zh}^2} \right]$。其中${x_{zh}} = (x - {z_1})/h$，${y_{zh}} = (y - {z_2})/h$，$\boldsymbol{z}({z_1},{z_2})$取决于评估点，h为节点间距^[25]。

2.2   改进CPDI的基本理论

与MPM不同，CPDI方法中有粒子域的概念，CPDI通过粒子域角点将粒子的信息映射到背景网格节点^[15-16]。CPDI将GIMP中的矩形粒子域转换为平行四边形，由向量$\boldsymbol{r}_1^{\text{o}}$和$\boldsymbol{r}_2^{\text{o}}$描述，如图 2（a）所示。在时间步t，更新后的粒子域由向量$\boldsymbol{r}_1^t$和$\boldsymbol{r}_2^t$描述，如图 2（b）所示。随着变形梯度$\boldsymbol{F}_p^t$的更新，粒子域变形为

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

Figure  2.  Initial and updated particle domains in CPDI method

下载: 全尺寸图片 幻灯片

在时间步t，粒子域角点坐标$\left\{ {\boldsymbol{x}_{c1}^t,{\text{ }}\boldsymbol{x}_{c2}^t,{\text{ }}\boldsymbol{x}_{c3}^t,{\text{ }}\boldsymbol{x}_{c4}^t} \right\}$可由粒子坐标$\boldsymbol{x}_p^t$和向量$\boldsymbol{r}_1^t$，$\boldsymbol{r}_2^t$表示：

CPDI的形函数${N_{Ip}}$和形函数的梯度$\nabla {N_{Ip}}$为

式中：${V_p}$为粒子p的体积；$(r_{1x}^t,{\text{ }}r_{1y}^t)$和$(r_{2x}^t,{\text{ }}r_{2y}^t)$分别为$\boldsymbol{r}_1^t$和$\boldsymbol{r}_2^t$的分量。

传统CPDI采用式（37）将粒子的动量映射到背景网格节点，粒子域内的速度局部为常数^[26]。由于粒子域内的速度场为零阶多项式，将导致空间计算精度的降低。本文提出的改进CPDI方法中采用了类似Wallstedt等^[27]的方法，利用速度梯度计算粒子域内的速度场，如图 3所示。改进CPDI中，每个粒子域的速度场由粒子的速度和4个角点的速度$({\boldsymbol{v}_p},{\text{ }}{\boldsymbol{v}_{c1}},{\text{ }}{\boldsymbol{v}_{c2}},{\text{ }}$ ${\boldsymbol{v}_{c3}},{\text{ }}{\boldsymbol{v}_{c4}})$组成，4个角点的速度见式（39）。利用粒子和4个角点的数据，可以采用2.1节中构建的自适应正交改进插值移动最小二乘形函数重构背景网格中的速度函数。

图  3  粒子域速度场

Figure  3.  Particle domain velocity field

下载: 全尺寸图片 幻灯片

首先将粒子信息映射到相关的背景网格节点，计算节点的质量$m_I^t$，外力$\boldsymbol{f}_I^{e,t}$和内力$\boldsymbol{f}_I^{i,t}$：

与传统CPDI不同，本文的改进CPDI方法中背景网格节点速度由粒子速度$\boldsymbol{v}_p^t$和4个相关粒子域角点速度$\boldsymbol{v}_{cj}^t(j = 1,{\text{ }}2,{\text{ }}3,{\text{ }}4)$采用AOIIMLS形函数更新，且每个粒子和粒子域角点有16个相关节点。在背景网格上显式求解运动方程，得到节点I的加速度$\boldsymbol{a}_I^t$和速度$\boldsymbol{v}_I^{t + \Delta t}$：

式中：t+Δt为下一个时间步；Δt为时间步长；${\phi _I}({\boldsymbol{x}_p})$和${\phi _I}({\boldsymbol{x}_{cj}})$为自适应正交改进插值移动最小二乘形函数；$\boldsymbol{v}_{cj}^t$为由式（39）计算得到的粒子域角点速度。

将计算得到的节点信息映射回粒子，可得粒子的速度$\boldsymbol{v}_p^{t + \Delta t}$和位置$\boldsymbol{x}_p^{t + \Delta t}$：

然后，更新粒子p的变形梯度$\boldsymbol{F}_p^{t + \Delta t}$，体积$V_p^{t + \Delta t}$和密度$\rho _p^{t + \Delta t}$：

式中：${\boldsymbol{E}_I}$为单位矩阵，$V_p^o$为粒子p的初始体积；$\nabla \boldsymbol{v}_p^{t + \Delta t}$为速度梯度，$\nabla \boldsymbol{v}_p^{t + \Delta t}{\text{ = }}\sum\nolimits_I {\nabla {N_I}({\boldsymbol{x}_p})} \boldsymbol{v}_I^{t + \Delta t}$。

最后根据式（35）和（36）更新粒子域的角点坐标$\boldsymbol{x}_{cj}^{t + \Delta t}(j = 1,{\text{ }}2,{\text{ }}3,{\text{ }}4)$和拓扑参数$\boldsymbol{r}_1^{t + \Delta t}$和$\boldsymbol{r}_2^{t + \Delta t}$。其中粒子域的角点速度$\boldsymbol{v}_{cj}^{t + \Delta t}(j = 1,{\text{ }}2,{\text{ }}3,{\text{ }}4)$由式（39）计算。

2.3   算例验证：一维柱在自重下的压缩

第一个验证算例为柱在自重下的压缩^[28-29]，如图 4所示，计算模型高${h_0} = 10.0{\text{ m}}$，宽${l_0} = 1.0{\text{ m}}$，柱的侧面约束x方向位移，底面约束x，y方向位移，采用线弹性本构模型。柱的弹性模量$E = {10^4}{\text{ Pa}}$，泊松比$\nu = 0$，密度$\rho = 80{\text{ kg}}{\text{/m}}{^3}$。背景网格为四节点正方形单元，单元尺寸h_{x, y}分别为1，1/2，1/4，1/8，1/16，1/32，1/64 m，每个背景网格中粒子数np分别为2²，3²，4²，5²，6²，7²，8²，9²，10²。模型的计算时间步长$\Delta t = 0.4{{{h_{x,y}}} \mathord{\left/ {\vphantom {{{h_{x,y}}} C}} \right. } C}$，其中$C = \sqrt {{E \mathord{\left/ {\vphantom {E \rho }} \right. } \rho }}$，总计算时间${T_t} = 3{\text{ s}}$，重力g从0增加到9.81 m/s²。

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

Figure  4.  Geometry and boundary conditions of column

下载: 全尺寸图片 幻灯片

在网格尺寸h_{x, y}=1/64 m，粒子数np=4条件下，MPM、CPDI和改进CPDI得到的柱的垂直应力数值解和解析解如图 5所示。在MPM计算结果中出现了明显的应力振荡，且计算结果与解析解有着较大的差别。CPDI和改进CPDI计算结果与解析解吻合较好，且改进CPDI得到的柱的垂直应力精度更高。

图  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)

下载: 全尺寸图片 幻灯片

柱的垂直应力解析解和数值解之间的误差为^[29]

式中：${(\sigma _{yy}^n)_p}$为粒子p的垂直应力数值解；$\sigma _{yy}^a(y_p^0)$为粒子p在初始位置$y_p^0$的垂直应力解析解；$V_p^0$为粒子p的初始体积；${V_0}$为柱的初始体积，${V_0} = \sum\nolimits_{p = 1}^{{n_p}} {V_p^0}$。

图 6，7为对数坐标系下MPM、CPDI和改进CPDI计算得到的不同粒子密度和不同网格尺寸下柱的垂直应力解析解和数值解的误差。考虑到计算效率，选取网格尺寸1.0 m分析粒子密度对垂直应力解析解和数值解的误差影响，如图 6所示，MPM计算误差随着每个背景网格中粒子数的增加先缓慢减小再迅速增大，表明通过适当增加粒子的密度可以降低数值噪声，提高计算精度，但当粒子密度过大时，MPM计算结果会发生更剧烈的数值噪声，导致精度大幅度降低。CPDI和改进CPDI计算误差随着粒子密度的增加而降低，同时改进CPDI得到的误差更低。在粒子密度为100时，改进CPDI得到的垂直应力解析解和数值解的误差比CPDI低约6%。如图 7所示，在粒子数np=4时，MPM计算得到的柱的垂直应力数值解和解析解的误差随着网格尺寸的减小而增大，分析原因为网格尺寸较小时，易出现单元穿越误差，导致严重的应力振荡，导致MPM计算精度降低。CPDI和改进CPDI得到的误差随着网格尺寸的减小而减小。网格尺寸为1.0 m时，CPDI和改进CPDI得到的误差基本相同，当网格尺寸进一步减小时，改进CPDI得到更小的误差。当网格尺寸为1/64 m时，改进CPDI得到的垂直应力解析解和数值解的误差比CPDI低约20%，表明改进CPDI方法具有更高的计算精度。

图  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)

下载: 全尺寸图片 幻灯片

2.4   算例验证：砂柱坍塌问题

第二个验证算例为Lube等^[30-31]研究的砂柱坍塌自由表面的演变，试验示意图如图 8（a）所示。砂柱的初始宽度和高度分别为0.0905，0.6335 m，试验开始时下落的重物迅速提起挡板，砂在重力作用下自由运动。砂柱的Mohr-Coulomb参数：剪切模量${G_{\text{c}}} = 0.323{\text{ MPa}}$，体积模量${K_{\text{c}}} = 0.7{\text{ MPa}}$，泊松比$\nu = 0.3$，黏聚力$c = 0{\text{ MPa}}$，内摩擦角$\varphi =$31°，膨胀角$\psi =$1°，初始密度$\rho = 2650{\text{ kg/m}}{^3}$。考虑到计算效率，本算例背景网格尺寸为0.02 m×0.02 m，每个背景网格中有25个粒子，试验模型共有3634个粒子，模型的网格和粒子初始位置如图 8（b）所示。

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

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

下载: 全尺寸图片 幻灯片

图 9为改进CPDI得到的砂柱坍塌过程中的总位移云图，图中的黑色实线为砂柱坍塌试验中砂柱自由面的位置^[30-31]。由图 9可知，在相同时刻，改进CPDI得到的模拟结果和试验结果基本一致，改进CPDI可以准确模拟砂柱坍塌的整个过程，最终砂柱坍塌的滑移距离为0.8 m。同时可以观察到，在计算过程中，一些粒子分离，这可以通过更精细的网格和更多的粒子改进。

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

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

下载: 全尺寸图片 幻灯片

图 10为CPDI和改进CPDI模拟得到的砂柱堆积体的自由面与试验所得自由面的对比。由图 10可知，在0~0.4 m范围内，CPDI模拟得到的砂柱自由面高度大于试验值，而改进CPDI模拟结果与试验值较一致。在0.4~0.8 m范围内，CPDI和改进CPDI模拟结果与试验值一致。结果表明，改进CPDI模拟结果比CPDI更符合试验结果，证明了改进CPDI的精度。

图  10  砂柱坍塌自由面对比

Figure  10.  Comparison of free surface of collapsing sand column

下载: 全尺寸图片 幻灯片

2.5   算例验证：隧道坍塌离心机试验

为了验证改进CPDI方法模拟隧道坍塌的可行性，将改进CPDI模拟结果与Kamata等^[32]的隧道坍塌离心机试验结果进行了对比。隧道的几何模型和边界条件如图 11所示，模型两侧边界约束水平位移，底部边界约束水平和垂直位移。隧道模型长度L为400 mm，高度H为240 mm，隧道覆盖层厚度C为80 mm，隧道直径D为80 mm，隧道开挖长度L_d为60 mm。改进CPDI的背景网格尺寸为5 mm×5 mm，每个背景网格中有4个粒子，试验模型共有4368个背景网格，14592个粒子。砂的Mohr-Coulomb模型参数根据参考文献[32，33]确定：弹性模量$E =$ $5.0{\text{ MPa}}$，泊松比$\nu = 0.3$，黏聚力$c = 0{\text{ kPa}}$，内摩擦角$\varphi =$34.5°，膨胀角$\psi =$0°，密度$\rho = 1510$ kg/m³。

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

Figure  11.  Geometry and boundary conditions of tunnel

下载: 全尺寸图片 幻灯片

图 12（a）为未加固隧道在重力加速度为30g时，铝板被拉出触发的坍塌模式，图 12（b）为改进CPDI计算得到的隧道坍塌后的总位移云图，图 12（b）中的红色实线表示离心机试验得到的隧道坍塌形态。由图 12可知，在改进CPDI模拟中，隧道开挖面附近的砂流入隧道，并在地表形成沉降槽。隧道破坏体的两个滑移面一个从开挖面顶部几乎垂直延伸到地面，另一个从开挖面底部斜向上延伸到地面，破坏区域位于开挖面前方约0.3D范围内。隧道坍塌总位移呈拱形向上延伸到地面，最大地面沉降量约为0.15D，改进CPDI模拟的隧道坍塌模式与离心机试验的坍塌模式接近。同时可以观察到，改进CPDI得到的沉降槽比离心机试验结果宽约0.25D，分析原因为改进CPDI采用平面应变模型导致的几何效应影响，可以将改进CPDI扩展到三维以更准确地模拟离心机试验，尽管如此，改进CPDI可以模拟离心机试验中的隧道坍塌破坏机制。

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

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

下载: 全尺寸图片 幻灯片

3.   工程应用

3.1   工程概况

青岛地铁4号线崂山区静港路站—沙子口站区间由两条平行的单洞隧道组成，隧道间距为13.8 m，左洞全长1123.531 m，右洞全长1143.346 m，采用盾构法+矿山法施工，其中矿山法隧道净宽7.4 m，隧顶埋深19.6 m，采用两台阶法施工，上台阶高约3.7 m，下台阶高约3.85 m，单台阶长度为5 m。隧道初期支护为50 cm间距的钢拱架和30 cm厚度的C25喷射混凝土，二次衬砌为30 cm厚度的C45钢筋混凝土。当隧道穿越围岩性质较差的区域时，在隧道拱部120°范围内施作长3.5 m，倾角15°的Φ42超前小导管进行注浆加固。

2019年5月27日，隧道施工到左洞ZDK25+343时发生了塌陷事故，最终形成如图 13所示的长30.6 m，宽25.5 m，深6.0 m的大型地面塌坑，这次灾难造成了严重的生命财产损失。左洞ZDK25+343塌陷处隧道采用矿山法施工，地层从上到下依次为6.4 m的杂填土层，0.8 m的粉质黏土层，7.1 m的中粗砂层，4.9 m的粉质黏土层和风化程度不同的凝灰岩层，地质剖面图如图 14所示，地层物理力学参数见表 1^[34-35]。

图  13  地面塌陷

Figure  13.  Ground collapse

下载: 全尺寸图片 幻灯片

图  14  静沙区间地质剖面图

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

下载: 全尺寸图片 幻灯片

表  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 

| 显示表格

3.2   隧道坍塌过程分析

选取图 14中的数值模拟段进行隧道坍塌全过程的分析，隧道模型和边界条件见图 15。模型尺寸为100 m×47 m，隧道埋深为19.6 m，模型两侧边界约束水平位移，底部边界约束水平和垂直位移，上部为自由面。在上台阶掌子面顶部及上方各地层分界处布设4个监测点监测围岩水平和竖向位移，监测点坐标分别为A（60，27.5），B（60，32.5），C（60，39.5），D（60，47）。本节模拟过程包括初始地应力平衡阶段和隧道开挖阶段。首先将重力线性加载，在10 s加载完成，同时对所有活动单元施加局部阻尼，得到满足准静态平衡条件的初始应力场。此时，开挖至ZDK25+343处，并将所有节点的速度重置为0，开挖后隧道顶板衬砌结构采用两层粒子实现，且约束水平和垂直位移，采用两台阶法开挖，两台阶的开挖面为自由面。采用Mohr-Coulomb模型，模型相关参数见表 1。背景网格单元尺寸为1 m×1 m，每个网格内有4个粒子。计算时间步长为10^-4 s，计算总时间为50 s。隧道坍塌过程中的地面沉降和剪应变如图 16，17所示。

图  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

下载: 全尺寸图片 幻灯片

隧道开挖以后，随着围岩变形逐渐增大，围岩将经历弹性、塑性和失稳的过程。如图 16所示，计算初始阶段，隧道上台阶掌子面前方地层向下挤压隧道产生沉降变形，而下台阶出现轻微隆起现象。随着计算的进行，沉降破坏区发展至地面，且呈现上宽下窄的“杯状”。最大竖向位移由上台阶掌子面前方岩体转移至地面。随着隧道的进一步坍塌，上台阶掌子面前方岩体流入隧道内部，下台阶在上方岩体的挤压作用下开始滑动。上台阶掌子面前方形成狭窄的流通道且竖向位移最大，分析原因为此处为粉质黏土地层，岩体强度较低，黏聚力较小，易流动发生破坏，最终形成深6 m的地面塌陷区。

如图 17所示，剪切带最先在上台阶掌子面的顶部和底部形成，同时下台阶发生剪切变形，距离上台阶掌子面前方15 m处地面出现竖向剪切带。随着隧道坍塌的进一步演化，上台阶掌子面处的剪切带向上延伸，最终形成“杯状”剪切带，导致了“杯状”沉降破坏区。图 18为隧道坍塌后地面沉降曲线，图中竖向虚线为上台阶掌子面位置，可以观察到地面存在塌陷区、变形区和稳定区。地面塌陷区位于上台阶掌子面前方15 m，后方10.3 m范围内，呈宽25.3 m，深6 m的“U”形。变形区位于上台阶掌子面前方15~30 m和后方10.3~25 m范围，且后方变形区的变形较大，与隧道坍塌现场结果基本一致。

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

Figure  18.  Curve of ground subsidence after tunnel collapse

下载: 全尺寸图片 幻灯片

监测点的竖向位移和水平位移时程曲线变化规律如图 19所示。由图 19可知，监测点的竖向位移和水平位移变化规律一致，呈缓慢变形—急剧变形—缓慢变形—变形稳定的变化趋势。由于上台阶掌子面为自由面，监测点A处围岩在上方岩体的挤压作用下向隧道内流动，围岩变形主要以水平位移为主，而监测点B，C，D的围岩变形主要以竖向位移为主。地面监测点D的竖向位移最大，为6 m，且几乎不产生水平位移。以上结果均与现场实测结果较吻合，该结果充分证明了改进CPDI方法在岩土工程大变形问题中的适用性及准确性。

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

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

下载: 全尺寸图片 幻灯片

4.   结论

本文在CPDI理论框架中引入AOIIMLS形函数，提出了一种改进CPDI方法用于岩土工程大变形数值模拟，降低了传统MPM和CPDI方法中的粒子穿越网格误差，并通过经典算例验证了改进CPDI方法的正确性以及模拟隧道坍塌的可行性，得到以下3点结论。

（1）AOIIMLS避免了求解逆矩阵，与MLS相比，提高了鲁棒性。同时满足Kronecker δ函数性质，允许直接施加本质边界条件。AOIIMLS形函数用于网格和粒子之间的插值和映射以减少单元交叉误差，提高了计算精度。

（2）改进CPDI采用速度梯度计算粒子域内的速度场，并采用粒子和粒子域角点速度构造背景网格的速度函数。算例验证表明，改进CPDI方法模拟结果与理论值和试验结果较一致，与MPM和CPDI相比，改进CPDI提高了计算精度。

（3）最后采用改进CPDI模拟了青岛地铁4号线静沙区间地面塌陷全过程，模拟结果与隧道坍塌现场基本一致，表明了改进CPDI方法在岩土工程大变形领域的适用性及优势，为岩土工程大变形计算领域提供一种新的数值计算方法，可为隧道施工及支护设计提供一定的指导作用。