煤系地层边坡开挖后，风化作用加强，当雨水渗入使得煤系土快速崩解软化，力学强度急速降低，边坡极易发生坍塌滑坡现象。目前，对煤系土的研究较少，尚处于起步阶段，而且现有的研究主要集中在对煤系土不良地质边坡的加固处理方法^[1-3]、边坡稳定性分析^[4-6]、以及煤系土强度特性等宏观层面上^[7-9]，从细观角度对煤系土边坡失稳的细观机理进行分析与解释的文献未有报道。煤系土边坡的失稳破坏运动是一个存在岩土体滑动、平移、转动的复杂过程，具有宏观上的不连续性和单个块体运动的随机性，用传统的宏观连续介质理论已经不能合理地分析散体过程，上述煤系土边坡风化之后的散、动特征都与传统的均匀、连续等假定冲突，将导致理论与实际的偏离。因此，要想揭示煤系土边坡的失稳机理，有必要从微细观角度进行研究^[10]。

降雨入渗是边坡发生失稳破坏的主要诱发因素之一，目前考虑流固耦合的计算方法主要有3类，第1类是Euler–Euler方法，该方法是将流体和固体都视为连续介质；第2类是Lagrange–Lagrange方法，该方法是把流体和固体都假设为离散介质；第3类是Euler–Lagrange计算方法，该方法将流体视为连续介质，采用流体动力学方法（computational fluid dynamic，CFD）进行模拟计算，固体视为离散介质，采用离散单元方法（discrete element method，DEM）进行模拟计算。这种DEM–CFD流固耦合方法早期主要应用于工业领域，Tsuji等^[11]首次将CFD方法和DEM方法进行耦合，对固体颗粒中的气体流动行为进行模拟，仿真计算结果很好地展现了DEM–CFD耦合方法模拟流体–固体相互运动问题的能力。Xu等^[12]提出了动力碰撞模型，将CFD和DEM的时间与空间尺度设置为一致，采用牛顿第三定律来描述流体与颗粒的相互作用。在此基础上，Xu^[13]、Yu^[14]等进一步将他们的方法进行了改进，对介于固体和流体之间的运动进行了模拟。

Shamy等^[15-16]首次将DEM–CFD耦合方法引入到岩土工程问题分析中，将Navier–Stokes方程进行了简化处理，采用3维DEM–CFD耦合方法计算分析了边坡渗流问题，对土体边坡渗流的细观机理进行了分析。国内一些研究人员在土体渗流和土体液化研究中也采用了DEM–CFD流固耦合的计算方法，取得了较理想的研究效果^[17-21]。蒋明镜等^[22]在2维离散单元软件（PFC2D）的基础上，将Tait状态方程引入到CFD，建立了弱可压缩流体DEM–CFD流固耦合模型，并通过模拟计算单个颗粒水下自由沉降和单面固结排水试验来验证所建耦合模型的可行性。Khalili等^[23]将离散化的DEM–CFD控制方程嵌入PFC2D软件中，对双轴不排水压缩试验进行了数值模拟计算。王胤等^[24]利用开源程序CFDEM建立了DEM–CFD流固耦合数值计算模型，并将颗粒滚动阻力机制引入到CFD–DEM耦合模型中，计算分析了滚动阻力模型的有效性以及所建立耦合模型的可行性。从以上文献可以发现，DEM–CFD耦合方法是一种非常有发展前景的数值模拟方法，特别是在土体与流体相互耦合作用的微细观分析方面具有明显的优势。

作者以江西省万载至宜春高速公路项目的煤系土边坡为研究背景，将离散元方法（DEM）与计算流体动力学方法（CFD）进行耦合，建立了煤系土边坡3维DEM–CFD流固耦合细观作用计算模型，分析降雨作用下煤系土边坡失稳破坏的细观机理，计算结果与模型试验结果比较吻合，研究成果可以为该区域煤系土边坡的防护设计与施工提供理论依据。

1 DEM–CFD流固耦合计算基本原理 1.1 流体–颗粒相互作用

采用流–固耦合理论（DEM–CFD）分析流体对岩土颗粒的作用，是将流体视为连续体，将岩土体视为不连续的颗粒（Euler–Lagrange法）^[25]。在CFD中流体的速度 ${{{v}}}$ 是宏观的速度，即单位横截面积上的体积流量，假定流体在整个横截面上产生，而实际流体流速只发生在颗粒孔隙空间中。作者将流体与颗粒之间的相互作用力进行了简化，仅考虑流体对颗粒的拖曳力和流体压力梯度力，流体与颗粒间作用方程如下^[26]：

$\frac{{\partial {{{u}}}}}{{\partial t}} = \frac{{{{{{{f}}}}_{{\rm{mech}}}} + {{{{{f}}}}_{{\rm{fluid}}}}}}{m} + {{{g}}}$

(1)

$\frac{{\partial {{{\omega }}}}}{{\partial t}} = \frac{M}{I}$

(2)

式（1）～（2）中， ${{{u}}}$ 为颗粒的速度， ${{{{f}}}_{{\rm{mech}}}}$ 为作用于颗粒的外力之和， ${{{{f}}}_{{\rm{fluid}}}}$ 为流体对颗粒的作用力（包括流体对颗粒的拖曳力和流体压力梯度力）， $m$ 为颗粒的质量， ${{{g}}}$ 为重力加速度， ${{{\omega }}}$ 为颗粒旋转角速度， $M$ 为作用于颗粒上的力矩， $I$ 为惯性矩。

假设流体施加于颗粒的力总是作用于形心，不考虑弯矩，则拖曳力 ${{{{f}}}_{{\rm{drag}}}}$ 为：

${{{{f}}}_{{\rm{drag}}}} = {{{{f}}}_0}{\varepsilon ^{ - \chi }}$

(3)

式中： ${{{{f}}}_0}$ 为单个颗粒所受的拖曳力； $\varepsilon $ 为流体所在单元的孔隙度； ${\varepsilon ^{{\rm{ - }}\chi }}$ 项为考虑局部孔隙度的经验系数，它可以使拖拽力同时适用于高孔隙度和低孔隙度的情况。

式（3）中单个颗粒所受的拖曳力为：

${{{{f}}}_0} = \frac{1}{2}{C_{\rm{d}}}{\rho _{\rm{f}}}\text{π} {r^2}\left| {{{{u}}} - {{{v}}}} \right|({{{u}}} - {{{v}}})$

(4)

式中， ${\;\rho _{\rm{f}}}$ 为流体密度，r为颗粒半径， ${{{u}}}$ 为颗粒的速度， ${{{v}}}$ 为流体速度， ${\;\rho _{\rm{f}}}$ 为流体的动力黏滞系数， ${C_{\rm{d}}}$ 为拖拽力系数，其表示为：

$ {C_{\rm{d}}}{\rm{ = }}{\left( {0.63{\rm{ + }}\frac{{4.8}}{{\sqrt {{R_{{\rm{ep}}}}} }}} \right)^2} $

(5)

式（3）中的经验系数 $\;\chi $ 为：

${\;\;\;\;\;\;\;\;\chi} = 3.7 - 0.65\exp \left( - \frac{{{{(1.5 - \lg {R_{{\rm{ep}}}})}^2}}}{2}\right)$

(6)

式中， ${R_{{\rm{ep}}}}$ 为颗粒的雷诺数，定义为：

${R_{{\rm{ep}}}} = \frac{{2{\rho _{\rm{f}}}r\left| {{{{u}}} - {{{v}}}} \right|}}{{{\mu _{\rm{f}}}}}$

(7)

施加在流体单位体积上的体力为：

${\;\;\;\;\;{{{f}}}_{\rm{b}}} = \frac{{\displaystyle\sum_{{j}} {{{{f}}}_{{\rm{drag}}}^{{j}}} }}{V}{\rm{ = }}\frac{4}{3}\text{π} {r^3}(\nabla p - {\rho _{\rm{f}}}{{{g}}})$

(8)

式中， $V$ 为流体单元体积， $\nabla p$ 为流体梯度，分子上求和对象为流体单元叠加的颗粒。

因此，流体对颗粒的相互作用力为：

$ {\;\;\;\;\;\;\;\;\begin{aligned}[b] {{{{f}}}_{{\rm{fluid}}}} =& {{{{f}}}_{{\rm{drag}}}} + {{{{f}}}_{\rm{b}}} = \frac{1}{2}{\varepsilon ^{{\rm{ - }}\chi }}{C_{\rm{d}}}{\rho _{\rm{f}}}\text{π} {r^2}\left| {{{{u}}} - {{{v}}}} \right|({{{u}}} - {{{v}}})+\\ &\frac{4}{3}\text{π} {r^3}(\nabla p - {\rho _{\rm{f}}}{{{g}}}) \end{aligned}}\!\!\!\!\!\!\!$

(9)

1.2 流体–流体相互作用

1）连续方程

平均的Navier–Stokes方程可以定量地描述孔隙流体特性，可用于流体的细观分析。流体的连续方程是根据质量守恒定律进行推导得出，如式（10）所示，其表达的物理意义为：单位时间内流入体积元的质量与流出体积元的质量差值等于体积元内质量的增加或减少。

$\frac{{\partial (\varepsilon {\rho _{\rm{f}}})}}{{\partial t}} + \nabla \cdot (\varepsilon {\rho _{\rm{f}}}{{{v}}}) = 0$

(10)

式中， $\nabla $ 为拉普拉斯算子， $\varepsilon $ 为流体所在单元的孔隙度。

2）动量方程

牛顿第二定律流体流动模型所获得的控制方程推导得到动量方程。Anderson等^[27]将外力在体积元内取平均，推导出流体的动量方程为：

$ \begin{aligned}[b] &\frac{{\partial (\varepsilon {\rho _{\rm{f}}}{{{v}}})}}{{\partial t}} + \nabla \cdot (\varepsilon {\rho _{\rm{f}}}{{ v v}}) - \varepsilon \nabla \cdot ({\mu _{\rm{f}}}\nabla {{{v}}})=\\ &\;\;\;\;\;\;\;\;\;\;\;\;\;\; - \nabla p - {{{{f}}}_{{\rm{drag}}}} + \varepsilon {\rho _{\rm{f}}} {{g}} \end{aligned}$

(11)

式中， $p$ 为流体压力。

1.3 颗粒运动控制方程

颗粒的运动方程是在牛顿第二定律的基础上推导出来的，为了考虑流体对颗粒的作用，在控制方程中引入颗粒与流体间的相互作用力 ${{{{f}}}_{{\rm{fluid}}}}$ ，得到：

$ {m_{{i}}}\frac{{{\rm{d}}{{{{{u}}}}_{{i}}}}}{{{\rm{d}}t}} = \sum\limits_{{{j}} = 1}^{{n_{{i}}}} {{{{{{f}}}}_{{{ij}}}}} + {{{{f}}}_{{{{\rm{fluid}}(i)}}}} + {{{{f}}}_{{{{{g}}(i)}}}} $

(12)

$ {I_{{i}}}\frac{{{\rm{d}}{{{{{\omega }}}}_{{i}}}}}{{{\rm{d}}t}} = \sum\limits_{{{j}} = 1}^{{n_{{i}}}} {{M_{{{ij}}}}} $

(13)

式（12）～（13）中， ${{{{u}}}_{{i}}}$ 为第 ${{i}}$ 个颗粒的速度， ${{{{f}}}_{{{ij}}}}$ 为作用于颗粒 $i$ 的接触力， ${{{{f}}}_{{{{{g}}(i)}}}}$ 为第 ${{i}}$ 个颗粒的重力， ${I_{{i}}}$ 为颗粒的转动惯量， ${M_{{{ij}}}}$ 为第 ${{i}}$ 个颗粒受到第 ${{j}}$ 个颗粒的转动力矩。

1.4 CFD–DEM耦合流程与程序实现

本文以离散单元法商业软件PFC3D为基础，采用Fish语言编写相关程序，使PFC3D与其内嵌的CFD求解器进行数据交换，从而实现流–固耦合计算。具体计算过程为：首先，根据边界条件，利用有限体积法将CFD模块中的流体控制方程式（10）、（11）进行离散，并采用PISO应力速度耦合算法进行求解^[28]，与此同时，将网格数据发送至DEM模块，通过DEM模块计算孔隙度和拖曳力并将数据发回CFD模块；然后，再将每一时间步控制单元内的流体与颗粒之间作用力发送至DEM计算模块，再运用离散元方法将流体作用力施加在土颗粒上，同时进行颗粒之间的力学计算；最后，将流体作用力及孔隙度传回CFD模型，以此循环计算直至结束。DEM–CFD耦合计算流程如图1所示，当两者计算时间相等时触发数据交互，完成流体与颗粒相互作用的计算过程。

2 煤系土边坡失稳DEM–CFD耦合数值模拟 2.1 煤系土边坡室外模型试验

为了研究降雨入渗条件下煤系土边坡失稳破坏机理，项目组进行了人工降雨室外边坡模型试验。试验边坡模型箱尺寸为4.5 m×3.0 m×2.4 m（长×宽×高），试验土样取自江西省万载至宜春高速公路项目A5标段K30+120处煤系土边坡开挖出露的原状土，边坡模型每40 cm填筑一层土体，坡率为1.0∶1.5。本次人工降雨模拟的降雨强度为江西宜春地区7月和8月份最大降雨强度为（9.84×10^–7 m/s）。试验从2018年11月19日上午8点开始，一直持续至晚上18点结束，现场观察边坡的失稳破坏模式主要是水土流失和雨水冲蚀，形成了不同深度的冲蚀沟，如图2（a）所示。

本次试验在边坡土体埋设了4行29个侧面示踪点，通过模型箱有机玻璃侧面可以清楚地观察示踪点的运动情况，如图2（b）所示。通过分析量测结果可知，在持续降雨过程中，边坡土体的主要破坏形式为雨水冲刷，在坡面形成深浅不一的冲蚀沟，最大冲蚀沟达0.36 m，发生在靠近坡中偏上的位置。详细的边坡模型试验结果分析参见文献[29]。

2.2 煤系土边坡数值计算模型建立 2.2.1 模型参数标定

作者采用数值模拟室内三轴试验的方法进行颗粒细观参数的标定。利用PFC^3D软件作为计算平台，自编程序建立三轴试验数值计算模型，试验模型墙体尺寸为高4.0 m，直径2.0 m的三轴墙体模型。为了节约计算时间和提高分析效率，将实际的土颗粒尺寸进行放大和集中处理。颗粒粒径在50～90 mm间等概率随机分布，生成的计算模型如图3所示。

根据煤系土的性质，标定与边坡数值建模均选用接触黏结模型，该黏结作用于颗粒间很小的接触点，可承受法向和切向应力，但不能传递弯矩^[30]，这与煤系土体性质比较接近。模型内生成颗粒之后，通过改变颗粒的细观参数反复试算，使得数值模拟的偏应力–轴应变曲线逼近室内试验结果，如图4所示。从图4中可以看出，模拟曲线的变化趋势及大小与室内试验结果相比都较为接近，说明该模型参数取值可以较好地用于模拟煤系土细观分析。颗粒参数细观标定结果如表1所示，其结果将直接用于边坡的数值模拟试验。

2.2.2 边坡计算模型建立

以室外人工降雨边坡模型试验为研究对象，通过fish语言编程，建立了3维边坡流固耦合（DEM–CFD）相互作用细观计算模型，如图5所示。图5（a）为雨水作用下边坡3维流固耦合计算模型，通过颗粒分组可以看出，整个模型颗粒分成两组，上部蓝色颗粒为雨水作用的饱和区域，生成2 646个颗粒；下部浅绿色颗粒为边坡土体的非饱和区域，生成4 572个颗粒。本次模拟的状态为模型边坡在人工降雨作用8 h后的工况，根据现场试验结果，雨水在模型边坡土体内形成了暂态饱和区，其渗透深度约为0.5 m，为方便建模及计算，饱和与非饱和土体的界限以直线代替曲线，土体渗流区设置见图5（b）。考虑到雨水作用于坡体之后，径流带动颗粒会在地面流动，所以将流体网格沿坡体前方进行了拓展，以符合实际情况，如图5（b）所示。在非饱和区域，其颗粒的细观参数按照第2.2.1节三轴试验标定的参数进行设置，由于饱和区的土体强度较非饱和土体强度低，故对饱和土体的黏结强度进行了折减，其余参数保持一致。

根据模型边坡人工降雨的流体特征及降雨强度，对作用于边坡饱和区域的流体参数按照表2所示的参数进行设置。

为使流体流动方向与实际一致， $ {\textit{Z}}$ 轴的流速与X轴的流速之比设为1.0∶1.5，流体流动的方向与室外降雨试验坡体的径流一致。

2.3 计算结果分析 2.3.1 煤系土边坡颗粒运动的宏观分析

1）边坡土颗粒运动轨迹计算结果分析

图6为模拟试验结束后从计算模型侧面观察到的土颗粒移动情况。从图6（a）可以看出，颗粒的移动主要发生在饱和区，非饱和区局部有颗粒移动，边坡整体会出现一个滑动带，这一现象可以由颗粒的移动清晰地看出。通过获取计算边坡模型的中截面，运用CAD将边坡的滑动带进行绘制，如图6（b）所示。从图6可以发现，上部土体发生的位移约为0.9 m，下部土体发生的位移约为1.7 m，同时形成一条近似直线状的滑动面。坡顶颗粒被冲刷，冲刷比较严重区域位于斜坡上部，冲刷深度在0.3～0.4 m之间，这与室外边坡模型试验结果比较一致。

2）边坡土体应力变化计算结果分析

为了解边坡土体在降雨入渗条件下内部受力变化，以便更好地分析土体内部的破坏机理，模拟计算时在土体内部设置了17个测量圆。对边坡施加流体作用后，通过模拟计算得到边坡土体X、Y、 $ {\textit{Z}}$ 3个方向应力变化监测数据，然后经surfer软件后处理生成云图，计算结果如图7～9所示。

通过比较各时刻边坡应力演变过程发现：当模型计算到20 000时步后，如图7所示，此时斜坡和坡顶上的应力较坡体内部应力小，应力集中主要发生在非饱和区，说明饱和区颗粒受流体作用，整体应力变小，稳定性变差；当模型计算到60000时步后，坡顶和斜坡上部的应力下降较多，坡脚部分应力增加，说明在这段时间坡体内部相对比较稳定，而饱和区土体遭到持续破坏，见图8；当模型计算到100 000时步后，此时边坡各方向应力变化主要发生在坡顶和坡脚，在坡脚处土体应力增加，坡顶部分土体应力减少幅度较大（图9），说明在雨水作用下，坡顶土体被大量冲刷至坡脚堆积，导致坡脚处颗粒大规模堆积。

2.3.2 煤系土边坡失稳破坏的细观机理分析

1）颗粒力链分析

边坡模型在降雨过程中力链的演变过程如图10所示。从图10可以看出，随着降雨的持续进行，边坡饱和土体区域的力链分布变得比较稀疏，这说明此部分土体颗粒间的接触力较小，表现出极不稳定的状态。在上部流体和颗粒的作用下，力链在非饱和土体下部区域表现得比较密集和粗大，说明边坡内部颗粒间的接触力较大，此区域土体表现得比较稳定。以此同时，力链沿坡体向前发生较大的延伸，斜坡变缓，坡顶及斜坡上力链减少，变得稀疏，说明饱和区颗粒发生了较大移动。当模型计算至100 000时步后，非饱和区土体颗粒的力链基本没有大的变化，如图10（d）所示。但处于饱和区的坡面上颗粒进一步减少，变得稀疏，说明坡面上的颗粒继续向下移动，造成坡脚处颗粒堆积增多，致使坡脚处和地面上的力链变粗以及变得密集。

2）颗粒配位数变化分析

降雨过程中边坡土体颗粒的配位数变化云图如图11所示。从图11中可以看出，随着降雨的持续作用，边坡土体颗粒的配位数逐渐变小，其中位于饱和区域坡顶和坡面上土体颗粒的配位数降幅较大，而坡体内部非饱和区域土体颗粒的配位数降幅较小，说明在雨水作用的区域，颗粒间接触不良，结构稳定性变差，在非饱和区的土体内部，颗粒间接触良好，该区域较为稳定。当模型计算到60 000时步后（图11（c）、（d）），饱和区的配位数均变得较小，尤其是斜坡顶部，而坡体内部和斜坡下部配位数较大，说明上部颗粒的接触间较差，在流体作用下，颗粒逐渐向下移动至坡脚处堆积，坡脚处变得密实，稳定性变好。

3）颗粒孔隙率变化分析

降雨过程中边坡模型孔隙率变化计算结果见图12。从图12可以看出，随着降雨持续进行，整个边坡土体的孔隙率都逐渐变大，其中，边坡饱和区部分（主要是坡顶位置）的孔隙率增加较大，非饱和区部分的孔隙率增加幅度较小。从局部上看，坡脚处土体的孔隙率在降雨过程中增加的幅度较小，当计算时步为60 000步以后，其孔隙率基本不再增加（最后增加至0.44），如图12（c）所示，说明在雨水作用下，坡顶处颗粒滑动，致使颗粒在跛脚处堆积，导致坡脚处土体孔隙率变化不大。而斜坡顶上部土体的孔隙率增加尤为明显，当计算时步为100 000步时，其孔隙率增加至0.8，如图12（d）所示，这说明在雨水的作用下坡顶土体颗粒发生移动，导致该区域土体颗粒变得疏松。

2.4 讨论

煤系土边坡由于其易风化，遇水容易软化和崩解，在计算模拟时将其视为非连续介质，采用离散元理论方法进行计算分析是比较合理的。本文通过采用离散元方法模拟的计算结果与室外降雨试验结果进行比较发现，两种研究方法得出的煤系土边坡的破坏模式都是雨水引起的冲刷破坏，边坡滑动面是近似的直线段，而且模型边坡雨水冲刷的范围和采用DEM–CFD耦合方法计算预测的滑动面位置非常接近，这说明本文采用的DEM–CFD流固耦合方法模拟煤系土边坡的降雨失稳破坏是合理的，这也验证了采用DEM–CFD流固耦合方法模拟分析岩土工程中的离散介质运动规律是可行的^[31]。

力链是反映颗粒运动宏观力学性能的重要指标，力链的粗细、位置变化、密集程度是颗粒受力及运动变化直观的显现。在降雨过程中，力链沿坡体向前产生较大的延伸，斜坡变缓，坡顶及斜坡上力链减少，且变得稀疏，说明饱和区颗粒极不稳定，整体发生了较大移动。而在上部流体和颗粒的作用下，力链在下部区域表现得比较密集和粗大，说明边坡内部接触力较大，此区域表现得比较稳定。在此次模拟边坡降雨试验中，通过力链变化的分析能够很好地解释雨水作用下边坡的破坏演变过程，对煤系土边坡的破坏机理有了更深入的理解。

3 结　论

1）室外降雨试验中，煤系土边坡坡率为1.0∶1.5时，其破坏的主要模式是雨水冲刷，在坡面形成深浅不一的冲蚀沟，最大冲蚀沟达0.36 m，发生在靠近坡中偏上的位置。本文采用DEM–CFD流固耦合方法模拟的煤系土边坡破坏形式与室外降雨试验结果比较一致，边坡滑动面预测结果为近似的直线段，与实际模型边坡雨水冲刷的范围非常接近，这说明采用DEM–CFD流固耦合方法对煤系土边坡的稳定性进行分析是可行的。

2）边坡土体颗粒的细观参数，如力链、配位数以及孔隙率，在降雨过程中都会随之发生变化，这些细观参数与边坡土体的宏观力学表现直接关联。因此，通过颗粒的细观参数变化分析，可以很好地解释雨水作用下煤系土边坡的破坏演变规律，从微细观角度分析非连续介质煤系土边坡的破坏机理是十分必要的。

3）本文研究结果表明，采用离散单元法和计算流体动力学方法相结合模拟分析离散介质边坡的稳定性是可行的，研究成果不仅为该区域煤系土边坡的防护设计与施工提供理论依据，而且为从微细观角度更好地分析离散介质岩土工程的宏观力学规律提供了一种新的思路。