隧道突泥破坏的耦合格子Boltzmann-离散元法模拟

金磊; 曾亚武; 程涛; 李晶晶

doi:10.11779/CJGE202106003

隧道突泥破坏的耦合格子Boltzmann-离散元法模拟 English Version

1.
湖北理工学院土木建筑工程学院,湖北　黄石　435003
2.
武汉大学土木建筑工程学院,湖北　武汉　430072

基金项目:

湖北省自然科学基金项目 2019CFB199

国家自然科学基金项目 41272342

国家自然科学基金项目 12002121

湖北省高等学校优秀中青年科技创新团队计划项目 T201823

详细信息

作者简介:
金磊(1989—),男,博士,副教授,主要从事岩土细观流–固耦合等方面的研究。E-mail: whujinlei@whu.edu.cn

中图分类号: TU457
计量
- 文章访问数: 274
- HTML全文浏览量: 58
- PDF下载量: 401
出版历程
- 收稿日期: 2020-09-24
- 网络出版日期: 2022-12-02
- 刊出日期: 2021-05-31

Numerical simulation of mud inrush of tunnels with coupled LBM-DEM

1.
School of Civil Engineering and Architecture, Hubei Polytechnic University, Huangshi 435003, China
2.
School of Civil Engineering, Wuhan University, Wuhan 430072, China

摘要

摘要: 为了探究隧道突泥灾害的灾变机制,开发了耦合格子Boltzmann–离散元法（LBM-DEM）数值计算平台并对隧道突泥破坏的演化过程进行了模拟,分析了致灾介质颗粒黏结强度、水压力及突泥口尺寸等因素对隧道突泥破坏特征的影响。结果表明：基于LBM-DEM模拟能再现隧道突泥破坏“启动、加速、缓慢和稳定”等连续4个阶段的演化过程；无黏结的致灾介质突泥破坏形态近似为直线,有一定黏结强度的致灾介质突泥破坏形态总体呈圆弧或抛物线状,突泥破坏区扩展范围和稳定后的突泥量随着颗粒间黏结强度的增大而逐渐减小；水压力越大,突泥灾害发生后突泥量增长越快,最终的突泥量也越大,且颗粒间黏结强度较大时水压力的这种影响越显著；当致灾介质颗粒间无黏结时,不同突泥口尺寸的模型在稳定后突泥量和破坏区范围基本相同,而当颗粒间形成一定强度的黏结后,突泥口尺寸越大,突泥灾害发生后突泥量增长越快,稳定后的突泥量也越多；隧道突泥破坏是致灾岩土介质、水压和开挖三者综合作用的结果。
- 隧道工程 /
- 突泥 /
- 格子Boltzmann方法（LBM） /
- 离散元法（DEM） /
- 颗粒–流体耦合
Abstract: To better understand the catastrophic mechanism of mud inrush disasters of tunnels, a numerical computing platform based on the coupled lattice Boltzmann method-discrete element method (LBM-DEM) is developed and used to simulate the evolution process of mud inrush of tunnels. According to the simulated results, the effects of particle bonding strengths of disaster-causing media, groundwater pressures and sizes of mud inrush holes on the characteristics of mud inrush of tunnels are analyzed. The results show that coupled LBM-DEM simulation can well reproduce the evolution process of four successive stages of mud inrush of tunnels: "starting, accelerating, decelerating and stabilizing". The failure form of unbonded disaster-causing media after mud inrush is approximately straight, whereas the failure form of disaster-causing media with certain bond strength is generally arc or parabolic. The expanded range of failure zone and mud inrush mass both decrease with the increase of the inter-particle bond strength of disaster-causing media. The higher the water pressure is, the faster the mud inrush mass increases after the occurrence of mud inrush disasters, and the more the final mud inrush mass is, which is more remarkable when the inter-particle bond is much stronger. When there is no bond between particles of the disaster-causing media, the models with different sizes of mud inrush holes have basically the same failure zone and mud inrush mass after stabilization. However, when a certain strength of bond is formed between particles, the mud inrush mass increases faster and the final mud inrush mass is more with the increase of the sizes of mud inrush holes. The mud inrush disasters are the result of the combined action of disaster-causing geo-materials, groundwater pressure and tunnel excavation.
- tunnel engineering /
- mud inrush /
- lattice Boltzmann method /
- discrete element method /
- particle-fluid coupling

HTML全文

0. 引言

在建筑工程领域中,对于水平荷载作用下桩基的受力分析,常见有极限地基反力法、弹性地基反力法、复合地基反力法、弹性理论法和p-y曲线法等^[1-2]。弹性地基反力法又包括地基系数常数法、k法、c法、m法以及吴恒立^[3]的双参数法。张有龄给出了地基系数为常数时的桩身响应解析解,N.B.ypdh与众多学者给出了桩身内力与变形的幂级数解,更有采用^[4]纽玛克法、有限差分法与有限元法来求解桩身内力与变形。

上述常用方法中对于多层地基情况的处理略显粗糙,如目前建筑桩基^[6]与公路桥涵桩基领域^[7]最常采用各层地基按其地基系数以权重进行折算,得到一个地基系数的等效值。近年来,Pise^[5]对双层地基水平受荷桩进行了数值求解,赵明华等^[8-10]对成层地基中桩的受力与试验做了大量工作,并尝试用无网格法分析计算,戴自航等^[11]采用有限元与有限差分进行数值计算,竺明星等^[12]利用矩阵传递法依次求解多层地基中的桩身各点内力,詹红志等^[13]也采用类似矩阵传递方法对抗滑桩嵌固段多层岩层进行了计算。

本文不同于先前学者从桩身形函数利用幂级数角度出发,引入张氏法的解析解函数形式,利用节点内力变形连续条件,建立全桩全节点统一矩阵线性方程,引入边界条件后一次性求解所有节点的变位与内力,并将该计算方法应用于多层地基桩基的水平响应计算。

1. 线性方程解力学模型

1.1 力学模型建立

不同于竺明星等^{[12, 14]}建立三参数地基系数模型并利用Laplace变换求解桩身响应的方法,本文在理论推导过程中不特别假定地基系数的分布模式,但考虑到设计人员使用上的便利性,以单层地基m值与多层地基m值分别演示计算过程。

根据Winkler理论,假定地基是服从胡克定律的弹性体,且每层地基厚度为h_j,如图1所示。

图 1 线性方程力学模型

Figure 1. Linear equation mechanics model

下载: 全尺寸图片幻灯片

将桩身沿深度方向分成 $n$ 段,桩单元依次编号为1,2,···, $n$ ,桩结点编号为0,1,2,···, $n$ ,结点对应各自坐标值。记桩身水平位移为 $y (z)$ ,桩身转角为 $φ (z)$ ,桩身弯矩为 $M (z)$ ,桩身剪力为 $F (z)$ , $M_{0}, F_{0}$ 表示桩顶作用的弯矩与水平力, $M_{i, z_{j}}, F_{i, z_{j}}$ 表示第i桩单元的 $z_{j}$ 节点处的弯矩与剪力,对于第 $i$ 桩单元,第 $i$ 段内地基系数 $k_{i}$ 以该段内的积分中值定理为原则,即

$k_{i} = \frac{\int_{z_{i - 1}}^{z_{i}} k (z) d z}{z_{i} - z_{i - 1}}$ 。

(1)

约定弯矩以桩左侧受拉为正,剪力以使桩顺时针转动方向为正,水平位移以坐标正向为正,而截面转角以逆时针转动为正。

1.2 微分方程及解答

对于第 $i$ 段,满足如下微分方程：

$E I \frac{d^{4} y_{i}}{d z^{4}} + k_{i} b y_{i} = 0$ 。

(2)

令 $A_{i} = \sqrt[4]{\frac{k_{i} b}{4 E I}}$ ,则可得到第 $i$ 段 $z \in [z_{i}, z_{i - 1}]$ 的挠曲线方程 $y_{i} (z)$ 解析解与挠曲线各阶导数：

$\begin{matrix} y_{i} (z) = C_{i}_{1} e^{- A_{i} z} \sin (A_{i} z) + C_{i}_{2} e^{- A_{i} z} \cos (A_{i} z) + \\ C_{i}_{3} e^{A_{i} z} \sin (A_{i} z) + C_{i 4} e^{A_{i} z} \cos (A_{i} z) \end{matrix}$ 。

(3)

将以上各函数表达式整理成矩阵形式,如式（4）：

$\begin{matrix} y_{i} = [\begin{matrix} f_{i 1} (z) & f_{i 2} (z) & f_{i 3} (z) & f_{i 4} (z) \end{matrix}] \cdot \\ {[\begin{matrix} C_{i 1} & C_{i 2} & C_{i 3} & C_{i 4} \end{matrix}]}^{T}, \\ y^{(1)}_{i} = [\begin{matrix} g_{i 1} (z) & g_{i 2} (z) & g_{i 3} (z) & g_{i 4} (z) \end{matrix}] \cdot \\ {[\begin{matrix} C_{i 1} & C_{i 2} & C_{i 3} & C_{i 4} \end{matrix}]}^{T}, \\ y^{(2)}_{i} = [\begin{matrix} p_{i 1} (z) & p_{i 2} (z) & p_{i 3} (z) & p_{i 4} (z) \end{matrix}] \cdot \\ {[\begin{matrix} C_{i 1} & C_{i 2} & C_{i 3} & C_{i 4} \end{matrix}]}^{T}, \\ y^{(3)}_{i} = [\begin{matrix} q_{i 1} (z) & q_{i 2} (z) & q_{i 3} (z) & q_{i 4} (z) \end{matrix}] \cdot \\ {[\begin{matrix} C_{i 1} & C_{i 2} & C_{i 3} & C_{i 4} \end{matrix}]}^{T} 。 \end{matrix}}$

(4)

令各系数矩阵表达式如下：

$\begin{array}{l} [P^{*}_{i, z_{i - 1}}] = [\begin{matrix} E I p_{i 1, z_{i - 1}} & E I p_{i 2, z_{i - 1}} & E I p_{i 3, z_{i - 1}} & E I p_{i 4, z_{i - 1}} \end{matrix}], \\ [Q^{*}_{i, z_{i - 1}}] = [\begin{matrix} E I q_{i 1, z_{i - 1}} & E I q_{i 2, z_{i - 1}} & E I q_{i 3, z_{i - 1}} & E I q_{i 4, z_{i - 1}} \end{matrix}], \\ [f^{*}_{i, z_{i - 1}}] = [\begin{matrix} f_{i 1, z_{i - 1}} & f_{i 2, z_{i - 1}} & f_{i 3, z_{i - 1}} & f_{i 4, z_{i - 1}} \end{matrix}], \\ [g^{*}_{i, z_{i - 1}}] = [\begin{matrix} g_{i 1, z_{i - 1}} & g_{i 2, z_{i - 1}} & g_{i 3, z_{i - 1}} & g_{i 4, z_{i - 1}} \end{matrix}] 。 \end{array}}$

(5)

则桩身 $n$ 段各个节点处的弯矩、剪力、挠度与转角记成： $[B] = [V^{*}] \cdot [C^{*}]$ ,其中

$[B] = [\begin{matrix} M_{1, 0} & ... & M_{n, n - 1} \\ M_{1, 1} & ... & M_{n, n} \\ F_{1, 0} & ... & F_{n, n - 1} \\ F_{1, 1} & ... & F_{n, n} \\ y_{1, 0} & ... & y_{n, n - 1} \\ y_{1, 1} & ... & y_{n, n} \\ φ_{1, 0} & ... & φ_{n, n - 1} \\ φ_{1, 1} & ... & φ_{n, n} \end{matrix}]$ ,

(6)

(7)

$[C^{*}] = {[\begin{matrix} [C_{1, j}] & [0] & ... & [0] \\ [0] & [C_{2, j}] & ... & [0] \\ [0] & [0] & ... & [0] \\ [0] & [0] & ... & [C_{n, j}] \end{matrix}]}_{j = 1, 2, 3, 4}$ 。

(8)

1.3 全桩线性方程组的构建与解答

保证桩身每一结点处内力与位移是连续的,以此思路建立桩身全结点的线性方程组如下：

$\begin{array}{l} M_{i, z_{i}} = M_{i + 1, z_{i}}, \\ F_{i, z_{i}} = F_{i + 1, z_{i}}, \\ y_{i, z_{i}} = y_{i + 1, z_{i}}, \\ φ_{i, z_{i}} = φ_{i + 1, z_{i}} 。 \end{array}}$

(9)

最终记成如下矩阵形式： $[ξ^{*}] \cdot [C] = [H]$ ,其中

$[C] = {[\begin{array}{l} \begin{matrix} [C_{1, j}] \\ [C_{2, j}] \\ [C_{3, j}] \\ ... \\ [C_{n - 1}_{, j}] \end{matrix} \\ [C_{n}_{, j}] \end{array}]}_{j = 1, 2, 3, 4}, [H] = [\begin{matrix} M_{0} \\ F_{0} \\ 0 \\ ... \\ 0 \\ H_{n - 1, z_{n}} \\ H_{n, z_{n}} \end{matrix}],$

(10)

$[ξ^{*}] = [\begin{matrix} [P^{*}_{1, z_{0}}] & [0] & [0] \\ [Q^{*}_{1, z_{0}}] & [0] & [0] \\ [P^{*}_{1, z_{1}}] & - [P^{*}_{2, z_{1}}] & [0] \\ [Q^{*}_{1, z_{1}}] & - [Q^{*}_{2, z_{1}}] & [0] \\ [f^{*}_{1, z_{1}}] & - [f^{*}_{2, z_{1}}] & [0] \\ [g^{*}_{1, z_{1}}] & - [g^{*}_{2, z_{1}}] & [0] \\ [0] & ... & ... \\ [0] & [0] & [ξ^{*}_{4 n - 1, n}] \\ [0] & [0] & [ξ^{*}_{4 n, n}] \end{matrix}]$ 。

(11)

上式矩阵运算表示了全桩全结点内力与位移值需要满足该线性方程组,引入桩顶与桩端的边界条件后,等式右侧矩阵也为常数阵,这样可通过Gauss消元等多种方法求解线性方程组,解得桩身每一段的四组参数 $C_{i 1}, C_{i 2}, C_{i 3}, C_{i 4}$ ,再将系数C矩阵回代式（8）即可。

2. 基于m法的设计计算与结果验证

2.1 单层地基算例验证

某建筑物^[2]采用桩基基础,直径d=1.5 m,埋入并支持在非岩石类土中,入土深度h=15 m,桩头在地面处自由,作用有水平荷载H₀=60 kN和M₀=700 kN·m,C25级混凝土的弹性模量E_c=2.8×10⁴ MPa= 2.8×10⁷ kN/m²,地基的反力系数的比例系数m=9400 kN/m⁴,土的内摩擦角 $φ$ =22°,黏聚力 $c = 15 {kN/m}^{2}$ ,重度 $γ = 20 k N / m^{3}$ 。

$b = K_{φ} K_{0} d = 0.9 (1.5 + 1) = 2.25 m$ ,

$E I = \frac{0.85 \times 2.8 \times 10^{7} \times π \times {1.5}^{4}}{64} = 59.3 \times 10^{5} kN \cdot m^{2}$ 。

分别将n取5,10,15,20进行了桩身内力与位移计算,本文法计算结果与传统m法计算的桩身弯矩绘制成曲线图,如图2所示。按规范法计算桩身最大弯矩为766.9 kN·m,将桩等分20段后桩身弯矩最大值为762.8 kN·m,相比规范法误差0.53%。从上图2看出,桩身弯矩随深度增加总体呈现先上升后下降的过程,弯矩极值出现在距离桩顶(1～3)d范围之间（d为桩身直径）,弯矩零点位于距离桩顶6d位置左右。

图 2 不同分段数桩身弯矩

Figure 2. Moments of pile body with different numbers of sections

下载: 全尺寸图片幻灯片

图3,4中看出3种端部约束下的弯矩、位移曲线重合度较高,仅在桩端附近处弯矩曲线出现了分叉发展的趋势。位移零点出现在距离桩顶4d位置处,相比桩身弯矩的6d变化范围缩小了33.3%。本算例所得的桩身弯矩极值、弯矩与位移零点所在的桩身位置符合目前国内外学者的研究结果,如赵明华等^[15-16]曾建议桩影响范围取3～5d,冯忠居等^[17]建议取（2～8）d等。

图 3 不同桩端约束的桩身弯矩

Figure 3. Moments of pile body restrained by different pile ends

下载: 全尺寸图片幻灯片

图 4 不同桩端约束的桩身位移

Figure 4. Displacements of pile body with different pile end constraints

下载: 全尺寸图片幻灯片

2.2 多层地基算例验证

某圆形^{[12, 18]}截面灌注桩^[10]桩径d=1.0 m,地面处桩顶剪力Q =150 kN,弯矩M =0,桩的弹性模量E =2.1675×10 kN/m²。桩侧有两层地基土体：第一层为流塑状回填土,层厚为2.0 m,相应的地基反力系数m为3000 kN/m⁴；第二层为硬塑状黏性土,桩身在该层土体中的长度为10.0 m,相应的地基反力系数为20000 kN/m⁴。

表1为不同计算方法的计算结果,从中可知本文线性方程解与精确解之间的桩顶位移误差为2.3%,最大弯矩误差为0.285%。

表 1 不同计算方法结果对比

Table 1. Results of different calculation methods

计算方法	桩顶位移/mm	最大弯矩/(kN·m)	最大弯矩位置/m
精确解^[18]	4.3735	336.81	3
规范解^[6]	3.0156	238.95	3
挠度加权换算^[9]	4.2093	343.83	3
有限差分法^[11]	4.2989	337.60	3
杆系有限元法^[11]	4.2520	336.07	3
矩阵传递法^[12]	4.2990	337.60	3
本文法	4.2703	335.85	3

下载: 导出CSV

| 显示表格

笔者在本算例基础上,改变地层情况再次进行桩身弯矩与位移计算,分别将地层视为全为上层土的单一地层与全为下层土的单一地层（简称“上层土地基”与“下层土地基”）,将计算结果分别绘制成曲线图5,6用以对比分析。由图5可以看出,桩身弯矩极值出现在距离桩顶（1～5）d之间,本例中双层地基情况与“下层土地基”情况都在8d位置处达到了弯矩零点；“下层土地基”与“两层土地基”均在深度5d处为位移零点。

图 5 不同地层情况桩身弯矩

Figure 5. Moments of pile body in different strata

下载: 全尺寸图片幻灯片

图 6 不同地层情况桩身位移

Figure 6. Pile displacements in different strata

下载: 全尺寸图片幻灯片

3. 结论

假设地基为弹性材料,分段建立梁挠曲线微分方程,通过结点内力与位移的连续条件一次性建立桩身全结点的线性方程组,求解得到各点内力与位移。以两个算例验证了线性方程解法在单层地基与多层地基中桩身响应计算的正确性,并对桩底不同边界条件、桩身周围不同地层进行了计算与讨论,得出如下结论：

（1）在桩顶水平荷载的作用下,桩身弯矩最大值出现在距离桩顶（1～5）d范围内,桩顶附近土层抗力越差,最大弯矩所出现的位置将越深。

（2）在距离桩顶（6～8）d位置附近将出现弯矩函数零点,且下降段所处区间受桩中部土层的抗力大小控制。

（3）桩身位移最大值出现在桩顶,距桩顶5d位置处出现位移零点。桩端不同的边界条件对桩身的位移影响较小,而桩周土层的抗力大小对桩身的位移起到控制作用。

（4）同一种情况下的桩身水平响应,其弯矩零点所出现的位置将比位移零点所出现的位置滞后（2～3）d。

本文在LBM-DEM耦合计算程序的开发和测试过程中得到了美国Los Alamos国家实验室Wang Min博士的悉心指导,在此致以衷心的感谢！

图 1 耦合LBM-DEM计算流程

Figure 1. Flow chart of coupled LBM-DEM

下载: 全尺寸图片幻灯片

图 2 DKT算例不同时刻的颗粒位置与流场

Figure 2. Position of particle and contours of fluid velocity at different instants for DKT case

下载: 全尺寸图片幻灯片

图 3 DKT过程中两颗粒的位置变化

Figure 3. Variation of horizontal and vertical positions of particles during DKT process

下载: 全尺寸图片幻灯片

图 4 隧道突泥的数值模型

Figure 4. Numerical model for mud inrush of tunnels

下载: 全尺寸图片幻灯片

图 5 不同黏结强度下的突泥量时变曲线

Figure 5. Time-varying curves of mud inrush mass of models with different inter-particle bond strengths

下载: 全尺寸图片幻灯片

图 6 B0-P2-O80模型突泥不同时刻的颗粒分布与流场

Figure 6. Distribution of particles and contours of fluid velocity at different instants during mud inrush for B0-P2-O80 model

下载: 全尺寸图片幻灯片

图 7 B14-P2-O80模型突泥不同时刻的颗粒分布与流场

Figure 7. Distribution of particle and contours of fluid velocity at different instants during mud inrush for B14-P2-O80 model

下载: 全尺寸图片幻灯片

图 8 不同黏结强度下模型最终的破坏形态

Figure 8. Final failure modes of models with different inter-particle bond strengths

下载: 全尺寸图片幻灯片

图 9 不同入口水压力下的突泥量时变曲线

Figure 9. Time-varying curves of mud inrush mass of models with different inlet water pressures

下载: 全尺寸图片幻灯片

图 10 不同入口水压力下模型最终的破坏形态

Figure 10. Final failure modes of models with different inlet water pressures

下载: 全尺寸图片幻灯片

图 11 不同突泥口尺寸下的突泥量时变曲线

Figure 11. Time-varying curves of mud inrush mass of models with different sizes of mud inrush hole

下载: 全尺寸图片幻灯片

表 1 隧道突泥的LBM-DEM模拟参数

Table 1 LBM-DEM simulation parameters for mud inrush of tunnel

固体参数	取值	流体参数	取值
颗粒密度/(kg·m^-3)	2650	流体密度/(kg·m^-3)	1000
法向刚度/(N·m^-1)	1×10⁷	运动黏度/(m²·s^-1)	1×10^-6
切向刚度/(N·m^-1)	5×10⁶	弛豫时间	0.5001
摩擦系数	0.5	格子步长/m	0.001

下载: 导出CSV

参考文献(27)

[1]	钱七虎. 地下工程建设安全面临的挑战与对策[J]. 岩石力学与工程学报, 2012, 31(10): 1945-1956. doi: 10.3969/j.issn.1000-6915.2012.10.001 QIAN Qi-hu. Challenges faced by underground projects construction safety and countermeasures[J]. Chinese Journal of Rock Mechanics and Engineering, 2012, 31(10): 1945-1956. (in Chinese) doi: 10.3969/j.issn.1000-6915.2012.10.001
[2]	李术才, 许振浩, 黄鑫, 等. 隧道突水突泥致灾构造分类、地质判别、孕灾模式与典型案例分析[J]. 岩石力学与工程学报, 2018, 37(5): 1043-1069. https://www.cnki.com.cn/Article/CJFDTOTAL-YSLX201805001.htm LI Shu-cai, XU Zhen-hao, HUANG Xin, et al. Classification, geological identification, hazard mode and typical case studies of hazard-causing structures for water and mud inrush in tunnels[J]. Chinese Journal of Rock Mechanics and Engineering, 2018, 37(5): 1043-1069. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YSLX201805001.htm
[3]	张志强, 阚呈, 孙飞, 等. 碎屑流地层隧道发生灾变的模型试验研究[J]. 岩石力学与工程学报, 2014, 33(12): 2451-2457. https://www.cnki.com.cn/Article/CJFDTOTAL-YSLX201412010.htm ZHANG Zhi-qiang, KAN Cheng, SUN Fei, et al. Experimental study of catastrophic behavior for tunnel in debris flow strata[J]. Chinese Journal of Rock Mechanics and Engineering, 2014, 33(12): 2451-2457. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YSLX201412010.htm
[4]	张庆松, 王德明, 李术才, 等. 断层破碎带突水突泥模型试验系统研制与应用[J]. 岩土工程学报, 2017, 39(3): 417-426. https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC201703007.htm ZHANG Qing-song, WANG De-ming, LI Shu-cai, et al. Development and application of model test system for inrush of water and mud of tunnel in fault rupture zone[J]. Chinese Journal of Geotechnical Engineering, 2017, 39(3): 417-426. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC201703007.htm
[5]	LIU J Q, CHEN W Z, LIU T G, et al. Effects of initial porosity and water pressure on seepage-erosion properties of water inrush in completely weathered granite[J]. Geofluids, 2018: 1-11.
[6]	张家奇, 李术才, 张庆松, 等. 基于透明土的隧道突泥破坏特征试验研究[J]. 中国公路学报, 2018, 31(10): 177-189. doi: 10.3969/j.issn.1001-7372.2018.10.017 ZHANG Jia-qi, LI Shu-cai, ZHANG Qing-song, et al. Experimental research on destruction characteristics of tunnel mud inrush using transparent soils[J]. China Journal of Highway and Transport, 2018, 31(10): 177-189. (in Chinese) doi: 10.3969/j.issn.1001-7372.2018.10.017
[7]	张庆艳, 陈卫忠, 袁敬强, 等. 断层破碎带突水突泥演化特征试验研究[J]. 岩土力学, 2020, 41(6): 1911-1923. https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX202006015.htm ZHANG Qing-yan, CHEN Wei-zhong, YUAN Jing-qiang, et al. Experimental study on evolution characteristics of water and mud inrush in fault fractured zone[J]. Rock and Soil Mechanics, 2020, 41(6): 1911-1923. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX202006015.htm
[8]	王媛, 陆宇光, 倪小东, 等. 深埋隧洞开挖过程中突水与突泥的机理研究[J]. 水利学报, 2011, 42(5): 595-601. https://www.cnki.com.cn/Article/CJFDTOTAL-SLXB201105015.htm WANG Yuan, LU Yu-guang, NI Xiao-dong, et al. Study on mechanism of water burst and mud burst in deep tunnel excavation[J]. Journal of Hydraulic Engineering, 2011, 42(5): 595-601. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-SLXB201105015.htm
[9]	ZHAO J D, SHAN T. Coupled CFD-DEM simulation of fluid-particle interaction in geomechanics[J]. Powder Technology, 2013, 239(17): 248-258.
[10]	蒋明镜, 张望城. 一种考虑流体状态方程的土体CFD-DEM 耦合数值方法[J]. 岩土工程学报, 2014, 36(5): 793-801. https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC201405002.htm JIANG Ming-jing, ZHANG Wang-cheng. Coupled CFD-DEM method for soils incorporating equation of state for liquid[J]. Chinese Journal of Geotechnical Engineering, 2014, 36(5): 793-801. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTGC201405002.htm
[11]	FENG Y T, HAN K, OWEN D R J. Coupled lattice Boltzmann method and discrete element modelling of particle transport in turbulent fluid flows: computational issues[J]. International Journal for Numerical Methods in Engineering, 2007, 72(9): 1111-1134.
[12]	BOUTT D F, COOK B K, MCPHERSON B J O L, et al. Direct simulation of fluid-solid mechanics in porous media using the discrete element and lattice-Boltzmann methods[J]. Journal of Geophysical Research, 2007, 112: B10209.
[13]	WANG M, FENG Y T, OWEN D R J, et al. A novel algorithm of immersed moving boundary scheme for fluid-particle interactions in DEM-LBM[J]. Computer Methods in Applied Mechanics and Engineering, 2019, 346: 109-125.
[14]	WANG M, FENG Y T, PANDE G N, et al. A coupled 3-dimensional bonded discrete element and lattice Boltzmann method for fluid-solid coupling in cohesive geomaterials[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 2018, 42: 1405-1424.
[15]	GALINDO-TORRES S A, SCHEUERMANN A, MUHLHAUS H B, et al. A micro-mechanical approach for the study of contact erosion[J]. Acta Geotechnica, 2015, 10: 357-368.
[16]	WANG M, FENG Y T, WANG C. Numerical investigation of initiation and propagation of hydraulic fracture using the coupled bonded particle-lattice Boltzmann method[J]. Computers & Structures, 2017, 181: 32-40.
[17]	BHATNAGAR P L, GROSS E P, KROOK M K. A model for collision processes in gases: I small amplitude processes in charged and neutral one-component systems[J]. Physical Review, 1954, 94(3): 511-525.
[18]	QIAN Y H, D’HUMIÈRES D, LALLEMAND P. Lattice BGK models for Navier-Stokes equation[J]. Europhysics Letters, 1992, 17(6): 479-484.
[19]	WANG M, FENG Y T, WANG C Y. Coupled bonded particle and lattice Boltzmann method for modelling fluid-solid interaction[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 2016, 40(10): 1383-1401.
[20]	WU J, SCHOLTÈS L, TINET A J, et al. A comparative study of three classes of boundary treatment schemes for coupled LBM/DEM simulations[C]//Proceedings of the 7th International Conference on Discrete Element Methods, Springer Proceedings in Physics, 2017: 551-560.
[21]	NOBLE D R, TORCZYNSKI J R. A lattice-Boltzmann method for partially saturated computational cells[J]. International Journal of Modern Physics C, 1998, 9(8): 1189-1201.
[22]	CUNDALL P A, STRACK O D L. A discrete numerical model for granular assemblies[J]. Géotechnique, 1979, 29(1): 47-65.
[23]	FENG Z G, MICHAELIDES E E. The immersed boundary-lattice Boltzmann method for solving fluid-particles interaction problems[J]. Journal of Computational Physics, 2004, 195: 602-628.
[24]	NIU X D, SHU C, CHEW Y T, et al. A momentum exchange-based immersed boundary-lattice Boltzmann method for simulating incompressible viscous flows[J]. Physics Letters A, 2006, 354: 173-182.
[25]	ZOU Q, HE X. On pressure and velocity boundary conditions for the lattice Boltzmann BGK model[J]. Physics of Fluids, 1997, 9: 1591-1598.
[26]	TANG Y, ZHU D Z, CHAN D H. Experimental study on submerged sand erosion through a slot on a defective pipe[J]. Journal of Hydraulic Engineering, 2017, 143(9): 04017026.1-04017026.14.
[27]	刘成禹, 张翔, 程凯, 等. 地下工程涌水涌砂诱发的沉降试验研究[J]. 岩土力学, 2019, 40(3): 843-851. https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX201903003.htm LIU Cheng-yu, ZHANG Xiang, CHENG Kai, et al. Experimental study of settlement caused by water and sand inrush in underground engineering[J]. Rock and Soil Mechanics, 2019, 40(3): 843-851. (in Chinese) https://www.cnki.com.cn/Article/CJFDTOTAL-YTLX201903003.htm