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

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

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
计量
- 文章访问数: 0
- HTML全文浏览量: 0
- PDF下载量: 0
出版历程
- 收稿日期: 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]。虽然目前对于岩石脆性没有得到大家一致认可的定义，但大量学者基于不同的研究领域，对脆性表征行为给出了各自的描述。代表性学术观点有：从应变角度，Morley^[5]和Hetenyi等^[6]将脆性定义为材料塑性的缺失；从强度角度，Ramsey^[7]认为岩石黏聚力丧失时，材料即发生脆性破坏。另外，岩石脆性也与破裂特征有关，如，李庆辉等^[1]从岩石非均质性、应力、破裂特征的角度认为脆性是其综合特性。同时，岩石脆性与其细观结构的非均质性相关^[8]。

多年来，根据上述特征学者们提出了许多定量评价岩石脆性的方法，并对其影响因素进行了分析。如，史贵才等^[9]获得了大理岩的应力脆性跌落系数，并讨论了其与围压的关系。周辉等^[10]通过考虑峰后应力降的相对大小和绝对速率，建立了可评价简单及复杂应力条件下岩石脆性评价指数。夏英杰等^[11]从岩石峰后应力跌落速率及能量比的角度建立了新的脆性指数。侯鹏等^[12]结合红砂岩和花岗岩试验提出了基于声发射和能量演化规律评价岩石脆性的方法。侯振坤等^[13]基于室内页岩提出了包含多参数的综合脆性指标，并从扩容拐点、残余强度、全过程曲线等多角度对脆性指标进行了诠释。

随着岩石起裂的深入研究发现，脆性岩石破坏征兆和断裂能及起裂密切相关^[14]。如George^[15]指出岩石脆性是指岩石内部微裂纹产生并引起非永久形变而连续变形的能力。据此可知，岩石的脆性与微裂纹产生有关，这就包括起裂时的特征应力和应变。另外，岩石脆性也是指其在特定条件下经能量积累和峰后过程自我维持宏观破坏的能力^[16]。从这个定义看出，岩石的脆性和材料弹性能量积累有关。纵观岩石破坏全过程，岩石破坏可以认为是微裂纹压密、起裂、发展、汇合并贯通的过程，同时又是能量积累、消耗的过程。在微裂纹起裂前，或在应力水平低于起裂应力水平前，岩石内部微裂纹不扩展，是弹性能积累的过程；当应力水平高于起裂应力时，岩石内部微裂纹开始起裂和发展，在裂纹起裂初始阶段，外界荷载做功将转化为压缩的弹性能和维持裂纹扩展的消耗能，弹性能持续增加但增长速率降低；当应力水平达到屈服时，岩石内部裂纹快速扩展，外界荷载做功主要转化为消耗能；当岩石破坏时，外界做功总能量均转化为消耗能，亦不存在弹性势能。这就可以认为在微裂纹起裂前，微裂纹闭合，外界做功转化为弹性能；之后微裂纹起裂并扩展，外界做功部分转化为消耗能，直至最后破坏全部转化为消耗能。结合George^[15]和Tarasov等^[16]对脆性的定义，以及岩石破坏微裂纹发展和能量转化过程，认为可以从岩石起裂以及峰后特征的角度对岩石脆性进行评价。如王宇等^[17]从起裂应力的角度建立了新的岩石脆性指数。但从起裂应变的角度建立岩石峰前脆性特征评价指数的研究鲜有报道。本文拟从起裂应变的角度建立岩石峰前脆性评价分量，然后结合峰后应力跌落特征建立峰后评脆性价分量，并综合峰前和峰后脆性评价分量建立基于岩石峰前起裂及峰后应力跌落特征的脆性评价指数。

1. 现有脆性评价方法简述

大量学者从强度、应变、能量、硬度等角度对脆性指数进行了分类评述，详见文献[1，10]。本文简述如下：基于强度的脆性指标因参数获取方便，被广泛地运用于岩石工程脆性评价，但该类指数忽略了应变–强度间的内在联系，存在一定的内在矛盾。大部分基于应变的脆性指标没有综合考虑岩石应力–应变全过程曲线，存在一定局限性。基于应变能的脆性指标是目前的一个发展趋势。基于硬度的脆性指标对非匀质性、各向异性的岩石测定结果差异大。基于莫尔包络线的脆性指标但受试验情况影响。基于特殊试验的脆性指标鉴于各种原因并未被广泛采用。

2. 岩石峰前起裂特征分析方法

目前，关于岩石起裂特征的分析方法主要有观察法，裂纹体积应变法、和声发射方法。其中以裂纹体积应变法运用最为广泛。本文将结合裂纹体积应变法判定岩石起裂特征。

裂纹体积应变法最早由Martin^[18]中提出。该方法确定起裂应力的思想为，裂纹体积应变（ ${\varepsilon _{{\text{V }}}}$ ）由体积应变（ ${\varepsilon _{\text{V}}}$ ）减去弹性体积应变（ ${\varepsilon _{{\text{V elastic}}}}$ ）。其中，体积应变（ ${\varepsilon _{\text{V}}}$ ）通过试验中测量的轴向（ ${\varepsilon _{{\text{axial}}}}$ ）和横向（ ${\varepsilon _{{\text{lateral}}}}$ ）应变计算得出，弹性体积应变由试验相应阶段应力状态（ ${\sigma _{\text{1}}}$ ， ${\sigma _{\text{3}}}$ ）和弹性参数（E，μ）计算求得，各应变计算公式如下：

${\varepsilon _{{\text{V }}}} = {\varepsilon _{{\text{axial}}}} + {\text{2}}{\varepsilon _{{\text{lateral}}}} \text{，}$

(1)

${\varepsilon _{{\text{V elastic }}}} = {\text{(1}} - {\text{2}}\mu {\text{)(}}{\sigma _{\text{1}}} - {\sigma _{\text{3}}}{\text{)}}/E \text{，}$

(2)

${\varepsilon _{{\text{V crack }}}} = {\varepsilon _{\text{V}}} - {\varepsilon _{{\text{V ealstci}}}} 。$

(3)

图 1为裂纹体积应变特征。岩石压密过程中，岩石内部微裂纹、空洞和空隙体积被不断压缩。到线弹性变形阶段，裂纹体积被压缩至极限状态，岩石类似于各向同性材料，仅有弹性体积应变，没有裂纹体积应变。之后，随着应力增加，荷载应力达到岩石起裂应力时，岩石内部裂纹开始扩展，裂纹体积逐渐增大。

图 1 裂纹体积曲线特征^[18]

Figure 1. Characteristics of crack volumetric strain^[18]

下载: 全尺寸图片幻灯片

3. 基于岩石峰前起裂及峰后跌落特征的脆性评价方法构建

岩石压缩过程中，起裂应力是裂纹压密和扩展的分界点，若脆性表征的是岩石维持裂纹压密状态弹性能聚集的能力，则此时弹性能增长速率达到峰值，之后增长速率降低，微裂纹扩展。若将起裂应变 ${\varepsilon _{{\text{ci}}}}$ 表征岩石前期可压密状态的能力，而将峰值应变与起裂应变的差值（ $\Delta \varepsilon$ = ${\varepsilon _{\text{c}}}$ - ${\varepsilon _{{\text{ci}}}}$ ）表征裂纹扩展的能力，其应变比（ ${\varepsilon _{{\text{ci}}}}$ / $\Delta \varepsilon$ ）就可以表征岩石维持内部裂纹压密状态的能力。起裂应变 ${\varepsilon _{{\text{ci}}}}$ 越大，表明岩石压缩做功过程越长，弹性能积累越久，裂纹可压密状态能力越强； $\Delta \varepsilon$ 越小，表明岩石裂纹扩展迅速； ${\varepsilon _{{\text{ci}}}}$ / $\Delta \varepsilon$ 越大，表明岩石维持裂纹压密状态的能力越大。可以用 ${\varepsilon _{{\text{ci}}}}$ / $\Delta \varepsilon$ 来表征岩石峰前的脆性特征（图 2）。

图 2 岩石峰前脆性表征示意图

Figure 2. Schematic diagram of brittleness for pre-peak

下载: 全尺寸图片幻灯片

一般而言，岩石在经历压缩试验后，破坏形式有Ⅰ、Ⅱ两种情况，（图 3（a））。通常认为Ⅱ型破坏是非稳定的，脆性较强，脆而坚硬的岩石往往属于此种变形破坏型式。为此，结合岩石破坏后其应力–应变曲线特征，若岩石脆性特征越明显，则其峰后应力跌落越明显。如图 3（b）所示，峰后曲线跌落越陡，则峰后曲线与斜边MN所组成的区域面积S越大，可以用S/S_△（S_△为图 3（b）中红色和蓝色区域面积）来表征岩石峰后的脆性特征。其中，关于图 3（b）中斜边M、N的确定，M为应力应变曲线峰值点，N为应力应变曲线的终点。同时三角形MCN为直角三角形。

图 3 岩石峰后特征及峰后脆性表征示意图

Figure 3. Post-peak characteristics and schematic diagram of brittleness for post-peak

下载: 全尺寸图片幻灯片

基于岩石起裂及峰后特征的脆性评价方法如下：

（1）根据岩石应力–应变曲线所反映的特征，初步判别其脆性特征。

（2）结合试验结果和式（1）~（3），求解起裂应变 ${\varepsilon _{{\text{ci}}}}$ 。具体如下：

根据试验过程中量测的轴向、横向应变和式（1），计算体积应变 ${\varepsilon _{\text{V}}}$ ；

根据试验结果，计算弹性模量E和泊松比 $\nu$ ，并结合式（2），计算弹性应变；

最后结合式（3），裂纹体积应变为体积应变减去弹性应变。初始阶段体积被压缩，体积应变逐渐减小。当处于线弹性阶段时，该阶段岩石可视为弹性体，裂纹不再被压缩，体积应变为0，并持续到线弹性阶段结束。在线弹性结束时，开始起裂，裂纹体积应变偏离0，该点也就是起裂点，据该特征确定起裂应变 ${\varepsilon _{{\text{ci}}}}$ 。

（3）结合峰值应变 ${\varepsilon _{\text{c}}}$ ，求取峰前脆性特征评价值B_i= ${\varepsilon _{{\text{ci}}}}$ / $\Delta \varepsilon$ 。

（4）结合峰后特征，求取面积S和面积S_△，并计算峰后脆性特征评价值B_ii=S/S_△。具体过程见下式：

${S_{ACMNB}} = \int {{\sigma _1}{\text{d}}{\varepsilon _1}} \text{，}$

(4)

${S_{ACNB}} = {\sigma _N}({\varepsilon _N} - {\varepsilon _M}) \text{，}$

(5)

$S_{△} = ({\varepsilon _N} - {\varepsilon _M})({\sigma _M} - {\sigma _N})/2 \text{，}$

(6)

$S = S_{△} - ({S_{ACMNB}} - {S_{ACNB}}) \text{，}$

(7)

式中， ${S_{ACMNB}}$ 为的面积， ${S_{ACNB}}$ 为的面积，S_△为三角形的面积，S为的面积。

（5）基于岩石起裂及峰后特征的脆性评价方法其计算表达式如下：

$B_{I}=B_{i}·B_{ii} \text{，}$

(8)

式中，B_i和B_ii分别为岩石峰前和峰后脆性特征值。用岩石峰前和峰后特征表征岩石全过程脆性特征有多种组合方式，如B_i+B_ii、B_i·B_ii或者指数等其他表达形式。将峰前峰后特征评价值之和作为脆性评价指数的有李庆辉等^[1]、夏英杰等^[11]；采用积的形式作为脆性评价指数的有周辉等^[10]。本文采用积的形式表征其脆性评价指数主要是基于以下因素考虑的，B_i= ${\varepsilon _{{\text{ci}}}}$ / $\Delta \varepsilon$ 作为峰前脆性评价值是基于其维持裂纹处于压密状态的能力，一般而言，脆性岩石在经历屈服后，在较小应变内，应力快速跌落， $\Delta \varepsilon$ 较前期起裂应变 ${\varepsilon _{{\text{ci}}}}$ 小，其比值B_i通常要大于1；而B_ii作为峰后脆性评价值，其通常要小于1。若将两者之和（B_i+B_ii）作为脆性评价指数，则主要体现的是峰前脆性特征，因为峰前脆性特征值B_i将明显大于峰后脆性特征值B_ii，这显得有失合理。若将两者之积（B_i·B_ii）作为脆性评价指数，据前面分析，岩石脆性越大，B_i和B_ii值都将越大，两者乘积也越大，这样可以消除两者量值大小明显差异带来的不利影响。

4. 脆性指数验证

结合岩石变形及起裂特征资料，求取B_i和B_ii值及其乘积。然后从不同围压、加载方式及岩性的角度对基于岩石起裂及峰后特征脆性评价方法进行验证。

4.1 不同围压下岩石脆性评价方法验证

图 4给出了5~35 MPa（90℃）条件下大理岩三轴压缩试验曲线。据图 4，在90℃条件下大理岩随着围压的增加塑性增强，在5 MPa时，峰后有一小段应力的特征，之后呈现出软化性质。15~35 MPa时，随着围压的增大，塑性变形区间增大，塑性增强。从定性角度认为，大理岩90℃时，塑性随着围压的增大而增大，这与前人研究成果一致。表 1给出了对应条件下脆性指数B_I的计算参数及结果，脆性指数B_I的值在5~35 MPa下分别为0.152，0.110，0.082，0.078，脆性指数B_I值随着围压的增大而逐步降低，计算结果与试验结果吻合，有效地体现了脆性指数随围压的变化。

图 4 90℃各围压条件下大理岩三轴压缩试验曲线图

Figure 4. Curves of conventional triaxial compression tests on marble under various confining pressures at 90℃

下载: 全尺寸图片幻灯片

表 1 大理岩常规三轴压缩试验B_I计算结果表(90℃)

Table 1. Calculated results of B_I of marble under conventional triaxial compression test (90℃)

围压MPa	应变比B_i= ${\varepsilon _{{\text{ci}}}}$ / $\Delta \varepsilon$	面积比B_ii	脆性指数B_I=B_i·B_ii	试验曲线定性分析
5	0.46	0.33	0.152	弹-脆-塑性
15	0.43	0.25	0.110	弹-塑性
25	0.48	0.17	0.082	弹-塑性
35	0.30	0.26	0.078	弹-塑性

下载: 导出CSV

| 显示表格

4.2 不同加载方式下岩石脆性评价方法验证

图 5和表 2给出了5 MPa围压90℃、130℃条件下三轴压缩和卸荷试验曲线和B_I的计算结果。据图 5，在90℃加载条件下，黑色曲线峰值后首先表现出应变软化的特征，然后有一段应力跌落的阶段，最后进入残余阶段，总体表现为应变软化的特征；在90℃卸荷条件下，红色曲线峰后应力垂直跌落，之后进入残余阶段，表现出明显的脆性特征。在130℃加载条件下，绿色曲线总体表现出典型应变软化的特征，表现为应变软化；在130℃卸荷条件下，蓝色曲线峰后应力跌落明显，且曲线“左凸”，脆性特征最为明显。结合曲线变形特征认为这4个试样脆性由大到小分别为卸荷130℃＞卸荷90℃＞加载90℃＞加载130℃, 其B_I值分别对应为0.94，0.65，0.44，0.40，计算结果和试验曲线结果一致，有效反映了卸荷条件下岩石脆性比加载条件下岩石脆性强的事实。

图 5 不同加载方式件下花岗岩试验曲线图

Figure 5. Granite test curves under different loading modes

下载: 全尺寸图片幻灯片

表 2 不同加载方式下花岗岩脆性指数B_I计算结果表(5 MPa)

Table 2. Calculated results of B_I of granite under different loading modes (5 MPa)

围压	温度/℃	应变比B_i= ${\varepsilon _{{\text{ci}}}}$ / $\Delta \varepsilon$	面积比B_ii	脆性指数B_I=B_i·B_ii	试验曲线定性分析
加载 5 MPa	90	0.74	0.6	0.44	弹-软-塑性
加载 5 MPa	130	0.93	0.43	0.40	应变软化
卸荷 5 MPa	90	0.88	0.74	0.65	弹-脆-塑性
卸荷 5 MPa	130	1.29	0.73	0.94	脆性

下载: 导出CSV

| 显示表格

4.3 不同岩性条件下脆性评价方法验证

图 6和表 3给出了60℃、5 MPa条件下各岩性三轴压缩试验曲线和B_I的计算结果。据曲线图，蓝色花岗岩曲线峰后应力跌落较砂岩、大理岩明显，同时具有一定“左凸”现象，是脆性特征的表现；黑色砂岩曲线峰后表现出分级跌落的特征，但每次跌落范围较小，整体上显示的是应变软化的特征；红色大理岩曲线为明显的应变软化特征。总体可以认为3种岩性60℃、5 MPa条件下脆性由大到小为花岗岩＞砂岩＞大理岩。表 3给出了脆性指数B_I计算参数及结果，B_I值由大到小分别为0.684，0.336，0.186，分别对应花岗岩、砂岩、大理岩，计算结果和试验曲线结果一致。

图 6 不同岩性条件下三轴压缩试验曲线图

Figure 6. Curves of conventional triaxial compression tests under different lithologic conditions

下载: 全尺寸图片幻灯片

表 3 不同岩性三轴压缩试验脆性指数B_I计算结果表(5 MPa)

Table 3. Calculated results of B_I for different lithologies under conventional triaxial compression test conditions (5 MPa)

岩性	温度 /℃	应变比 B_i= ${\varepsilon _{{\text{ci}}}}$ / $\Delta \varepsilon$	面积比 B_ii	脆性指数 B_I=B_i·B_ii	试验曲线定性分析
花岗岩	60	0.94	0.73	0.684	弹-脆-塑性
砂岩	60	1.40	0.24	0.336	应变软化
大理岩	60	0.84	0.22	0.186	应变软化

下载: 导出CSV

| 显示表格

5. 讨论

目前，岩石脆性评价指数较多，本文将选取现有常用5种脆性评价指数与新建立的脆性评价指数（B_I）进行对比分析。鉴于岩石材料本身的非均质性和各项异性，为此选取基于同一试样或基于应力应变的脆性评价指数作为对比讨论分析的对象。选取的5种脆性指数表达式见下表 4。图 7给为脆性指数B₃₁，B₃₅和B₃₈的计算模型示意图。

表 4 用于对比分析的脆性指数列表

Table 4. Brittleness indices for comparative analysis

脆性指数	来源	公式	变量含义
B₃₁	刘恩龙等^[19]	B₃₁=1-exp(M/E)	M为软化模量，E为弹性模量
B₃₅	周辉等^[10]	${B_{35}} = \frac{{{\tau _{\text{p}}} - {\tau _{\text{r}}}}}{{{\tau _{\text{p}}}}}\frac{{\lg \left\| {{k_{ac(AC)}}} \right\|}}{{10}}$	${\tau _{\text{p}}}$ ， ${\tau _{\text{r}}}$ 分别为峰值强度和残余强度，k_{ac (AC)}为屈服起始点至残余起始点连线的斜率
B₃₆	Altindag 等^[20]	${B_{36}} = ({\tau _{\text{p}}} - {\tau _{\text{r}}})/{\tau _{\text{p}}}$	${\tau _{\text{p}}}$ ， ${\tau _{\text{r}}}$ 分别为峰值强度和残余强度
B₃₇	Altindag 等^[20]	${B_{37}} = ({\varepsilon _{\text{p}}} - {\varepsilon _{\text{r}}})/{\varepsilon _{\text{p}}}$	${\varepsilon _{\text{P}}}$ ， ${\varepsilon _{\text{r}}}$ 分别为峰值应变和残余应变
B₃₈	史贵才等^[9]	${B_{38}} = \frac{{({\varepsilon _{\text{B}}} - {\varepsilon _{\text{P}}})}}{{({\varepsilon _{\text{P}}} - {\varepsilon _{\text{M}}})}}$	${\varepsilon _{\text{P}}}$ ， ${\varepsilon _{\text{B}}}$ 分别为峰值应变和残余应变，ε_M为残余强度对于的峰前应变

下载: 导出CSV

| 显示表格

图 7 部分脆性评价指数计算模型示意图

Figure 7. Schematic diagram of some brittleness evaluation indices

下载: 全尺寸图片幻灯片

同时，刘恩龙等^[19]基于脆性指数B₃₁的研究结果，建立了相应的脆性评价等级分类建议值，结果分为6个等级，见表 5。

表 5 脆性指数B₃₁分类等级及描述^[19]

Table 5. Classification and descriptions of brittleness index B₃₁^[19]

等级	脆性指数B₃₁	脆性描述
1	1	M→-∞，理想脆性
2	0.6321~1	-∞＜M＜-E，脆性很大
3	0.6321	M=-E
4	0~0.6321	-E＜M＜0，脆性很小
5	0	M=0，理想塑性
6	＜0	M＞0，应变硬化

下载: 导出CSV

| 显示表格

下面将从不同围压、加载方式、岩性角度对比分析B_m，B₃₁，B₃₅~B₃₈的脆性评价结果。

图 8给出了90℃各围压条件下大理岩三轴压缩试验不同脆性评价方法的结果。B₃₁，B₃₅和B₃₆等3个脆性评价指数随围压增加，总体呈降低特征，能表达出大理岩随围压升高而降低的性质，结合表 5，可知B₃₁的评价和试验结果吻合较好。B₃₇随围压增大呈现增大趋势，据B₃₇表达式可知其值越小脆性越明显，也能在一定程度体现大理岩随围压升高脆性降低的特点。但B₃₈计算值对大理岩在90℃各围压条件下的脆性评价不是十分敏感。对比分析认为B_I，B₃₁，B₃₅~B₃₇能用于评价岩石脆性随围压的变化，B₃₈有待进一步分析。

图 8 大理岩常规三轴压缩试验各脆性指数计算结果图(90℃)

Figure 8. Calculated results of each brittleness index for marble under conventional triaxial compression test (90℃)

下载: 全尺寸图片幻灯片

图 9给出了围压5 MPa不同加载方式下花岗岩各脆性评价方法的结果。B₃₁在卸荷条件下其值明显大于压缩条件下计算值，最大最为90℃时的0.73，据表 5可判定其脆性很大。但针对卸荷5 MPa 130℃时的评价与试验实际不符（计算值0.57，据表 5应为脆性很小，实际试验脆性特征明显）。总体上B₃₁能在一定程度上体现不同加载方式对岩石脆性的影响。B₃₅~B₃₈在评价不同加载方式对岩石脆性影响时，计算结果规律性不强，对加卸载条件下岩石脆性评价有一定价值，但规律性不敏感。据图 9可认为，B_I和B₃₁可体现不同岩性其脆性的差别。

图 9 不同加载方式下花岗岩各脆性指数计算结果图

Figure 9. Calculated results of each brittleness index for granite under different loading modes

下载: 全尺寸图片幻灯片

图 10给出了5 MPa、60℃三轴压缩条件下不同岩性各脆性评价方法的结果。据试验结果可知，在5 MPa三轴压缩条件下各岩性脆性关系为花岗岩＞砂岩＞大理岩。B₃₁中砂岩计算值最大（0.88），若据表 5可以判断其为理想脆性，与实际试验不符，对花岗岩和大理岩评价准确。B₃₅~B₃₈计算结果规律性不强，计算结果与试验有一定差异，对不同岩性评价不敏感。据图 10可知，B_I可体现不同岩性脆性的差别。

图 10 不同岩性常规三轴压缩试验各脆性计算图(5 MPa)

Figure 10. Calculated results of each brittleness index for different lithologies (5 MPa)

下载: 全尺寸图片幻灯片

6. 结论

（1）将起裂应变 ${\varepsilon _{{\text{ci}}}}$ 表征岩石前期可压密状态的能力，将峰值应变与起裂应变的差值（ $\Delta \varepsilon$ = ${\varepsilon _{\text{c}}}$ - ${\varepsilon _{{\text{ci}}}}$ ）表征裂纹扩展的能力，提出了表征岩石峰前脆性特征的分量B_i（B_i= ${\varepsilon _{{\text{ci}}}}$ / $\Delta \varepsilon$ ）。若 ${\varepsilon _{{\text{ci}}}}$ / $\Delta \varepsilon$ 越大，表明岩石维持裂纹压密状态的能力越大，峰前脆性越大。

（2）基于岩石应力–应变曲线峰后应力跌落越明显其脆性越大的特征，提出了表征岩石峰后脆性特征的分量B_ii（B_ii=S/S_△）。若S/S_△越大，表明岩石峰后应力跌落越迅速，峰后脆性越大。

（3）建立了可表征峰前起裂及峰后应力跌落特征的脆性指数B_I（B_I=B_i·B_ii）。

（4）脆性指数B_I能有效体现不同加载方式、围压以及岩性对岩石脆性的影响。其他条件相同的情况下，三轴卸荷条件下岩石脆性比三轴压缩条件下大，低围压下岩石脆性比高围压下大。

图 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

施引文献(10)

期刊类型引用(7)

1.	陈贺，解纹，余相贵，尹淏，秦雨樵. 层状碎裂结构边坡灾变过程的离散元数值分析. 公路交通科技. 2025(03): 115-125+164 . 百度学术
2.	范超强，齐宏，刘家魁，田勇，孙宏泰，缪玉松. 青海油田某区块深层岩石可压性分析. 青岛理工大学学报. 2024(05): 43-51 . 百度学术
3.	冯龙飞，王双明，王晓东，解嘉豪，窦林名. 煤单轴峰后动态冲击破坏特征及差异机制模拟研究. 煤炭学报. 2024(S2): 714-730 . 百度学术
4.	蒋培林，丁书学，南华，乔静，刘江，田永超，李佳赛，李树杰. 中间主应力影响下加锚试样应力-应变曲线特征参数分析. 科学技术与工程. 2023(18): 7725-7733 . 百度学术
5.	辛子朋，柴肇云，孙浩程，李天宇，刘新雨，段碧英. 砂质泥岩峰后破裂承载特征与块体分布规律研究. 岩土力学. 2023(08): 2369-2380 . 百度学术
6.	宋昊，夏敏，任光明. 基于峰前能量演化及峰后侧向变形特征的岩石脆性评价方法. 成都理工大学学报(自然科学版). 2023(06): 723-733 . 百度学术
7.	牛子豪，朱珍德，阙相成，卢文斌. 基于三维多晶离散元柱状节理岩体压缩特性研究. 岩石力学与工程学报. 2022(09): 1874-1887 . 百度学术