基于LBM-DEM耦合计算模型的含水层内悬浮颗粒运动特性

马玖辰，崔阿凤，吕林海，杨 杰

(1.天津城建大学 能源与安全工程学院，天津 300384；2.天津大学 中低温热能高效利用教育部重点实验室，天津 300072；3.天津城建大学 地热高效利用技术研究中心，天津 300384)

含水层储能作为节能减排、能源转型的关键技术，在工业余热回收利用、电力消纳、清洁建筑供热等领域得到了广泛的应用，为促进“碳达峰、碳中和”目标做出重要贡献。然而，黏土矿物、砂岩颗粒等悬浮颗粒会在水动力或应变力的作用下，侵入流道或孔隙空间，堵塞多孔介质的孔喉，对地下含水层的渗透性造成不可逆的损害[1]。因此，研究悬浮颗粒在多孔介质中的迁移和沉积对含水层储能技术[2]、地下水回灌[3]和地下污染物扩散[4]等领域具有重要的理论研究与实践应用价值。目前，国内外学者主要采用理论模型、柱状层析柱试验、显微可视化技术及数值模拟计算等方法开展含水层内颗粒运移研究。

由于具有相同渗透率的多孔介质通常存在不同的孔隙结构，因此，国内外学者通常开展层析柱试验确定多孔介质的孔隙度、粒径及孔隙结构，作为在宏观尺度下建立含水层颗粒运移理论模型的基础。薛传成等[5]通过对悬浮颗粒的迁移过程进行研究，认为渗流溶液温度与pH值、颗粒种类是影响多孔介质中悬浮颗粒迁移规律的重要因素。Chen等[6]分析了地下水渗流过程对于颗粒释放的影响。张鹏远等[7]通过砂柱试验，得到颗粒滞留程度随着粒径增大而提高。Bai[8-9]等通过理论与试验研究相结合，得到颗粒迁移过程随着溶液离子浓度与含水介质孔隙结构的变化而改变。 马玖辰等[10]结合电镜扫描与颗粒表面ζ电位测试，研究了球形与非球形微纳米颗粒形貌变化与其释放、运移、沉积过程的内在关联。蒋思晨等[11]通过开展不同形状颗粒在多种溶液离子强度下的渗透试验，将穿透曲线与DLVO理论相结合，确定颗粒形貌对其迁移特性具有重要影响。虽然，粒子图像测速（РIV）[12]、粒子阴影图测速（РSV）[13]、CT扫描[14]等显微可视化技术得以在试验中应用，但是仍难以揭示流体与颗粒的相互作用机理。

当前，针对悬浮颗粒在多孔介质中运移、重组过程的研究，国内外学者提出格子Boltzmann方法（lattice boltzmann method，LBM）、计算流体力学（CFD）与离散单元法（DEM）耦合计算方法。格子Boltzmann方法从孔隙尺度研究多孔介质中的颗粒滞留和堵塞，多用于纳米流体动力学[15]与非均质多相流[16]。尽管LBM方法能很好地表示弯曲颗粒的无滑移边界条件局部速度分布，但是通常忽略了粒子运动对流场的影响[17]。CFD-DEM方法可以在动力场、颗粒运动轨迹、应力变化等方面对边坡降雨[18]、泥沙运移[19]、颗粒污染[20]等领域开展有效的研究。在DEM计算方法中，分别针对每个悬浮颗粒的受力形式与运移过程建立控制方程，与LBM耦合计算可以精准描述孔隙中固相颗粒与液相粒子的运移行为[21]。因此，LBM-DEM耦合计算方法在简化控制方程、并行计算性能方面具有显著优势[22]。国内外学者采用LBM-DEM方法在路基土壤流态化[23]、固体流变[24]和隧道突泥[25]等方面开展深入研究。Zhou等[26]研究了渗流溶液浓度、流速和pH对多孔介质堵塞机理的耦合影响。Zhou等[27]通过建立不同渗透率的双通道模型，揭示了颗粒直径、流速、颗粒体积分数和注入量对非均质储层颗粒截留率和渗透率损伤的影响。

综上，由于多孔介质孔隙结构具有多样性，宏观尺度的运移模型与层析柱试验难以有效揭示含水层孔隙结构变化的诱导机制；而数值模拟计算存在并行计算能力弱、收敛性差等问题。同时，采用LBMDEM方法从介观尺度研究颗粒形貌对其在含水层中的释放—运移—沉积过程影响的研究尚少见报道。本文基于多孔介质渗流溶液与悬浮颗粒运移、受力控制方程，采用LBM-DEM耦合方法，将含水层流场划分成规则的网格，通过引入流体虚拟粒子，将渗流速度离散到多个方向。利用浸没移动边界法（IMB）处理固相颗粒与液相粒子的相互作用，从介观尺度对悬浮颗粒在多孔介质中的运动过程开展数值计算。搭建渗流层析柱试验系统，预测含水层水动力场演化过程，同时验证理论模型与求解方法的正确性。本文将介观尺度数值计算与试验研究相结合，探究颗粒形貌、地下水动力和多孔介质尺寸效应等因素对含水层渗透性能的影响及其内在诱导机制。

1 理论基础与控制方程

1.1 含水层渗流控制方程

本文以连续介质运移理论为基础，运用连续性方程（式（1））和Navier-Stokes方程（式（2））[28]来描述渗流溶液运动的流场，在Navier-Stokes方程中，引入外部体积力Ff（式（3））体现外部动力和流体相互作用：

式（1）～（4）中：ρ为流体密度，kg/m3；ε为多孔介质孔隙率；t为时间，s；v为流体速度矢量，m/s；P为流体压力，Рa；υ为流体运动黏性系数，m2/s；K为多孔介质渗透率；为多孔介质颗粒平均直径，m；G为外力项，N。

1.2 基于LBM的流体控制方程

通过引入虚拟粒子对流体速度进行离散，采用LBM模拟含水层中流体粒子的运动过程，对格子Boltzmann演化方程（式（5））进行求解，确定各网格节点处粒子的迁移和碰撞过程[29]。

式中，x为位置矢量。

在每个时间步长Δt内，x处流体粒子以离散速度矢量ci运动到相邻的节点上，继而在该节点处与其他粒子发生碰撞，碰撞后更新粒子速度分布函数fi(x,t)。流体粒子发生一系列迁移碰撞，逐渐趋向平衡态分布，平衡态分布函数(x,t)（式（6））可由离散的Maxwell分布函数表示：

式中，wi为i方向的权重因子，cs为格子声速。

基于格子单松弛模型（lattice bhatnagar-grosskrook，LBGK），将粒子分布函数向平衡态分布函数的时间松弛简化为单个BGK算子，用于代替复杂碰撞过程；考虑到定压物理边界条件和流体相互作用力，在演化方程中引入离散力源项Fi：

式中，τ为无量纲松弛时间。

利用演化方程求解虚拟粒子的碰撞和迁移过程后，由速度分布函数fi(x,t)的第0和第1速度矩求得流体的宏观参数，包括局部密度ρ(x,t)、局部流速v(x,t)：

运动黏性系数υ与格子间距Δx、时间步长Δt有关。通过Chapman-Enskog展开[30]，由演化方程与Navier-Stokes方程可得到运动黏性系数υ与无量纲松弛时间τ的关系式为：

式中，c为格子速度，即格子间距与时间步长之比。

1.3 基于DEM的颗粒运移控制方程

采用基于软球模型的离散单元法（discrete element method，DEM）模拟固体悬浮颗粒间的相互作用。所有的颗粒运动都遵循牛顿定律，允许彼此有一定重合度，将颗粒相互作用简化为弹簧-阻尼器[31]过程。因此，颗粒运动的动力学方程分为平移和旋转，考虑到固-液相互作用、固-固接触以及重力等因素，颗粒的平移由法向接触力Fcn、切向接触力Fct、拖拽力Fd、流体与颗粒相互作用力Ff,p、重力Fg、浮力Fb等力共同作用；颗粒的旋转由颗粒碰撞力矩Tp和流体流动施加在颗粒上的力矩TH共同作用，如式（11）、（12）所示：

式（11）～（12）中，u为颗粒平移速度，m为颗粒质量，IР为转动惯量，ωp为角速度。

图1为颗粒在孔隙中运动的受力示意图。

图1 悬浮颗粒受力示意图Fig.1 Schematic diagrams of the forces acting on the suspended particles

图2 试验装置示意图Fig.2 Schematic diagram of experimental device

根据Hooke定律可知，颗粒间的法向接触力与切向接触力可以通过颗粒间的法向弹簧与切向弹簧来模拟，如式（13）、（14）所示：

当区域Re数较低且颗粒不靠近通道壁时，流体受到Stokes阻力，此时颗粒受到拖拽力Fd：

在渗流溶液中，颗粒与流体粒子相互碰撞，从而产生的布朗运动会影响悬浮颗粒的迁移，导致颗粒受到相互作用力Ff,p：

颗粒在地下含水层中受到的重力Fg和浮力Fb的关系可由式（17）表示：

式（13）～（17）中：kn和kt为法向和切向接触刚度；δn和δt为法向重叠量和切向位移；ηn和ηt为法向和切向阻尼系数；Δun和Δut为法向和切向相对速度；d为当量直径；kB为Boltzmann常数；T为温度；α、σ分别为高斯分布函数的均值和标准差，设置为0和1；V为颗粒等效体积；ρp为悬浮颗粒密度；g为重力加速度。

颗粒碰撞力矩与转动惯量可用式（18）和（19）计算：

式（18）～（19）中，R为从颗粒中心指向接触点的位置向量。

2 试验系统与研究内容

2.1 试验装置

试验采用长300 mm、内径25 mm、壁厚3.5 mm的有机玻璃层析柱进行强制渗流试验。层析柱水平放置，在出、入口处均放置0.01 mm厚、孔径为0.18 mm的滤网。选用RSC型双头蠕动泵将水容器里高纯度去离子水以恒定转速注入层析柱；在LSР01-2A型注射泵的驱动下，将悬浮液以相同的速度注入层析柱；使用精确度为 ±2%的TL2350浊度仪测量渗流液浊度。

将高纯度去离子水不断注入特定浓度悬浮液进行稀释，测量相应悬浮液浊度，得到颗粒浓度与浊度的关系曲线如图3所示。将试验结果的标定浓度与浊度值进行线性回归分析，其判定系数R2＞0.99，说明两者具有高度相关性。由图3可知：颗粒的浓度随浊度的增大而增大；随着浊度的增加，球形硅微粉与非球形硅微粉的浓度差随之增加。由于非球形硅微粉以片状、多面体结构为主，比表面积大于球形硅微粉，对光线的散射作用增强。因此，在相同浓度下，非球形硅微粉悬浮液的浊度远高于球形硅微粉悬浮液。

图3 悬浮颗粒浊度与浓度的关系曲线Fig.3 Relationship curves between turbidity and concentration of suspended particles

2.2 试验材料

选用人工制备砂样作为层析柱中填充的多孔介质，为了确保填充材料颗粒级配连续性与矿物构成与地下含水层原始砂样相近，根据原砂的机械组成[10]，将等效粒径75 µm和120 µm的天然石英砂，按照1∶4比例均匀混合。试验前，将选取的石英砂进行酸洗，使用去离子水反复冲洗，烘干备用。采用密度为2.26 g/cm3、中位粒径为13 µm的球形和非球形硅微粉作为微纳米颗粒物质[10]。利用场发射扫描电子显微镜（SEM），观察微纳米颗粒形貌特征与粒径尺寸分布发现：球形硅微粉形状统一为圆球状，比表面积小，具有良好的分散性；非球形硅微粉形貌、几何尺寸呈不规律分布，多为多面体状颗粒，如图4所示。

2.3 试验方案与结果分析

采用分段湿填法填充4根层析柱，确保每根层析柱中含水介质重度达到（1.65±0.02） kg/L，反演计算得到孔隙率ε近似为0.37。分别配置800 mL浓度为0.25 mg/mL的球形与非球形硅微粉悬浮液，在渗流速度为2.5×10-5m/s与5.0×10-5m/s这两种工况下，开展4组强制渗流试验。首先，关闭阀门C，打开阀门A、B，采用蠕动泵以固定流速向层析柱中连续通入高纯度去离子水，当流出液的渗流速度稳定，并且无颗粒出现时，认为层析柱中含水介质的饱水、排气过程完成。随即关闭阀门A、B，打开阀门C，使用注射泵分别以2.5×10-5m/s与5.0×10-5m/s这两种流速将制备好的球形与非球形硅微粉悬浮液连续注入层析柱。待悬浮液注入完毕后，关闭阀门C，打开阀门A、B，继续以相同流速无间断通入高纯度去离子水。4组试验运行时长为180 min，每组试验均在室温（22～24 ℃）下进行。试验过程中每4 min采用秒表量筒法测试渗出液的流速；同时，收集渗出液，使用浊度仪测量溶液浊度，根据拟合关系曲线（图3）得到渗出液中的颗粒浓度，通过抽滤称重法测定4组试验在不同时段内的渗出液中悬浮颗粒的累计质量。

释放颗粒浓度与累计质量动态变化曲线如图5所示。

由图5（a）可知：以渗流速度为5.0×10-5m/s向层析柱注射球形硅微粉悬浮液时，渗出液颗粒浓度变化曲线呈双峰状分布状态，当试验进行到88 min时，释放颗粒浓度达到最大值0.379 mg/mL。当渗流速度降低至2.5×10-5m/s时，渗出液颗粒浓度无明显峰值出现，在试验进行的52～132 min，释放颗粒浓度曲线在0.087～0.167 mg/mL范围内呈不规则波动。向层析柱内注入非球形硅微粉悬浮液时，渗出液的颗粒浓度变化规律与注入球形硅微粉悬浮液时一致。由图5（b）可知：渗流速度为5.0×10-5m/s时，释放颗粒浓度最大值降低到0.11 mg/mL，同时峰值时间出现延迟；渗流速度为2.5×10-5m/s时，释放颗粒浓度曲线在0.019～0.035 mg/mL范围内波动，波动时间增长20 min。

表1为4组层析柱渗流试验结果。由表1可知，随着渗流速度的提高，颗粒释放率增大。以渗流速度为5.0×10-5m/s为例，注射球形与非球形硅微粉悬浮液相比较，微纳米颗粒释放率提高了47.9%。试验结果表明，颗粒形貌与水动力对颗粒的迁移特性有重要影响。与球形硅微粉相比，非球形硅微粉比表面积更大，沉积在多孔介质表面或被孔喉捕获的机率更高。由于渗流剪切应力与颗粒法向截面面积成正比[10]，因此，在相同流速下，沉积于多孔介质表面的非球形硅微粉难以脱离释放。

表1 不同形貌颗粒释放率与渗出液颗粒浓度Tab.1 Release rates of particles and particle concentration of the seepage solution with different morphology

3 数值计算与分析

3.1 求解方法

渗流溶液流体粒子、流体粒子与悬浮颗粒及悬浮颗粒之间的受力模式都是引起多孔介质内悬浮颗粒迁移状态改变的诱导机制。因此，采用欧拉坐标下的LBM模拟计算流体粒子流动过程[16]；采用拉格朗日坐标下的DEM软球模型模拟计算颗粒迁移过程[32]；同时引入浸没移动边界方法（immersed moving boundary method，IMB）处理流体粒子-固体颗粒边界，从而克服了LMB与DEM模型的动量不连续问题，实现对悬浮颗粒运移过程的准确数值计算。

基于LBGK模型提出包含外力项的流体粒子速度场演化方程（式（5）），其中，权重因子wi和离散速度矢量ci均与网格离散化相关[31]，取决于使用的格子模型。因此，根据层析柱实际渗流过程，选择2维平面具有9个离散速度的D2Q9模型对演化方程进行离散求解。每个方向的离散速度和权重因子由式（20）、（21）确定：

由于IMB方法可以在离散格子中同时包含流体粒子和固体颗粒，因此，在LBM-DEM耦合求解中，通过无滑移边界条件来实现固-液相互作用。在IMB中引入格子固含率εs表征格子节点被颗粒覆盖的体积分数，其中，流体内部节点εs为0，颗粒内部节点εs为1；当流体与颗粒边界碰撞时，可以确定每个节点的εs取值（0～1）。根据IMB计算原理，将流体粒子速度场演化方程（式（5））进行优化修正（式（22）），引入碰撞项Ωis（式（23））用以实现固体颗粒对流体粒子流动的影响：

式（22）～（23）中：±i表示正（反）方向；B为Ωis的加权函数，由式（24）确定：

基于优化的速度场演化方程（式（22）），结合流体粒子相互作用的离散力、流体粒子对于悬浮颗粒的拖拽力、颗粒间相互作用力，确定流体粒子施加于固体颗粒上的力FH，即对附加碰撞项在流体边界、颗粒边界及颗粒内部节点的所有方向所引起的动量变化求和，流体流动施加在颗粒上的力FH和力矩TH可表达为：

式中（25）～（26），n为该悬浮颗粒所标记节点的数量，(xk-xp)为从颗粒中心指向第k个节点的位置向量。

3.2 LBM-DEM计算平台与验证

根据强制渗流试验系统中层析柱的几何结构，将计算模型设置为300 mm×25 mm的长方形，如图6所示。根据所确定的格子距离，将计算区域离散为1.2×1011个节点（1 200 000×100 000）。将计算区域的上下边缘均设置为封闭的实体边界，采用反弹边界格式处理；渗流溶液与悬浮颗粒从左侧进入、右侧流出，因此左、右侧为速度边界，采用周期边界格式处理。根据试验结果，通过随机四参数法在区域内生成多孔介质，确定含水介质的孔隙率；同时分别设置7.68×107个中位粒径为13 µm的球形和非球形颗粒作为注入含水介质中的悬浮颗粒。模拟计算的基本参数设置见表2。

表2 数值模拟基本参数Tab.2 Basic parameters of the numerical simulation

图6 2维数值计算模型Fig.6 Two-dimensional numerical model

根据所建立的演化方程与求解方法，在数值计算平台MATLAB上，开发了LBM-IMB-DEM耦合求解程序。设定流体粒子与固体颗粒的初始条件，在每个时间步长上，利用LBM更新流场变化，利用DEM循环计算颗粒的碰撞。通过计算流体粒子作用在固体颗粒上的力和力矩，更新颗粒位置信息，完成流体与颗粒的耦合，耦合计算流程如图7所示。由于在LBM计算中，网格划分形式与数量对计算结果影响很大，因此通过反复的预处理计算确定格子距离在5×10-7～1×10-6m内，与孔隙率及粒子作用力无关。由于LBM与DEM均选用显式循环算法，并且具有各自的时间步长，因此在求解过程中取Δt/ΔtD=100，将DEM作为LBM的子循环从而确保耦合计算时间的统一。

图7 LBM-IMB-DEM耦合流程图Fig.7 Flowchart of the LBM-IMB-DEM coupling scheme

利用LBM-IMB-DEM耦合求解程序，计算得到4组试验下层析柱出口处渗流速度vout逐时变化过程。选择渗流速度vout的绝对误差与均方根误差作为计算值与试验结果的相似度判定指标。如图8所示，4组试验下模拟结果在计算周期中的动态变化与试验数据基本一致，均可全程跟踪试验过程，vout的均方根误差均小于10%，最大标准误差为0.26。为进一步验证模型的正确性，利用本文所写程序对Zhou[26]和Feng[33]等中的DKT算例进行模拟，图9为不同时刻颗粒流速分布对比，结果显示三者的模拟结果吻合度较高。因此，认为所建立的演化方程与耦合求解方法可以准确地反映悬浮颗粒在含水介质中运移、沉积及再释放的重组过程。

图8 层析柱出口渗流速度动态变化曲线Fig.8 Curves of the seepage velocity determined at chromatographic column outlet

图9 DKT算例不同时刻的颗粒流速分布Fig.9 Velocity of the particles at different for DKT case

3.3 计算结果与分析

图10～11表示在不同时刻下，层析柱含水介质中渗流速度的演化过程。由于在LBM的计算过程中需要对流场初始化，将流体粒子节点均设置为初始渗流速度，因此在初始状态（0 min）下含水介质流体粒子的渗流速度与入口边界处赋值相同。随着悬浮颗粒溶液的注入，LBM-IMB-DEM耦合求解程序启动，对比20 min与40 min这两个时刻，在4组试验中含水介质中下游区域的渗流速度均出现了显著降低。当再次注入去离子水时，渗流速度出现回升，截至计算周期结束（180 min），含水介质流场的整体分布趋于一致。

图10 注射球形硅微粉工况含水介质流速分布Fig.10 Distribution of liquid velocity in the aquifer medium with injection of spherical silica powders

由图10可知，向含水介质注入球形硅微粉悬浮液时，随着进水流速的提高，渗流速度恢复程度增加。悬浮颗粒会与多孔介质反复碰撞，即颗粒之间以及流体粒子与颗粒发生碰撞，造成动能的损失，若在此过程中不能到达分离点，颗粒最终将沉积在多孔介质表面，导致渗流速度下降。提高渗流速度，将使颗粒所受法向接触力Fcn（式（13））、切向接触力Fct（式（14））、拖拽力Fd（式（15））增大，使颗粒快速到达分离点，增强溶液对颗粒的携带效果，使更多颗粒流出含水介质，流体粒子运动空间增加，渗流速度得以恢复。

由图11可知，向含水介质注入非球形硅微粉悬浮液时，渗流速度降幅增大，仅在高流速工况下的进口处出现微弱回升。相同等效体积V下，圆球状的球形硅微粉比多面体状的非球形硅微粉的当量直径d小，导致球形硅微粉的转动惯量Ip（式（19））减小，同时颗粒中心到含水介质接触点的距离增加，使颗粒碰撞力矩Tp（式（18））增大，因此球形硅微粉易发生滚动、脱离。注入非球形硅微粉悬浮液时，由于相互作用力Ff,p（式（16））增大，引起流体粒子与颗粒更剧烈的碰撞，动能损失增大，渗流速度降低。非球形硅微粉比表面积大，易絮凝成团，在狭窄孔隙喉道及凹陷处被捕获的机率更大，有效孔隙空间的减少增强了对颗粒的阻碍作用，导致颗粒的自由迁移路径更少，从而引起含水介质渗透性下降。

图11 注射非球形硅微粉工况含水介质流速分布Fig.11 Distribution of liquid velocity in the aquifer medium with injection of non spherical silica powders

选择含水介质距进水端50、150、250 mm处作为研究对象，计算确定各断面全部流体粒子的平均渗流速度v¯，从而拟合生成动态变化曲线如图12所示。

图12 vin=2.5×10-5 m/s时含水介质不同断面平均渗流速度动态变化曲线Fig.12 Curves of the average seepage velocity in different sections of the aquifer medium with vin=2.5×10-5 m/s

向含水介质注入球形硅微粉悬浮液时（图12（a）），速度变化曲线呈高斯分布，各断面平均流速最大降幅分别达到12.6%、29.2%、38.3%，计算周期结束时稳定在2.43×10-5、2.36×10-5、2.34×10-5m/s。向含水介质注入非球形硅微粉悬浮液时（图12（b）），速度变化呈玻尔兹曼分布，各断面平均流速依次下降至2.0×10-5、1.7×10-5、1.3×10-5m/s，并且持续保持稳定。

图13为进水速度提高到5.0×10-5m/s时，含水介质不同断面平均渗流速度动态变化曲线。如图13所示，当进水流速提高到5.0×10-5m/s 时，含水介质渗流速度变化均呈现下降、回升、稳定3个连续阶段。以距进水端150 mm处为例：当向含水介质注入球形硅微粉悬浮液时，平均流速下降至3.1×10-5m/s，在117 min便回升至4.8×10-5m/s；当向含水介质注入非球形硅微粉悬浮液时，平均流速最低下降至2.7×10-5m/s，计算周期结束时仅回升至3.4×10-5m/s。

图13 vin=5.0×10-5 m/s时含水介质不同断面平均渗流速度动态变化曲线Fig.13 Curves of the average seepage velocity in different sections of the aquifer medium with vin=5.0×10-5 m/s

以进口流速为5.0×10-5m/s，向含水介质注入非球形硅微粉悬浮液为例，3组断面平均流速最大降幅依次提高，分别为27.1%、43.7%、56.3%；计算周期结束时平均流速恢复率则逐渐降低，依次达到62.5%、30.4%、28.1%。通过分析可知，当悬浮颗粒大量注入含水介质时，流体与颗粒的运动空间在一定时间内迅速减少，随着颗粒在含水介质中的运动路程增大，颗粒之间以及颗粒与多孔介质的碰撞率提高，动能损失增大，从而引起含水介质中、后部平均流速降幅增大，导致悬浮颗粒沉积量增加。

4 结 论

1）通过引入离散力源项，得到具有外力项的LBGK速度演化方程；根据DEM的软球模型，对含水介质中悬浮颗粒的平移和旋转过程进行受力分析，建立LBMDEM耦合计算模型。利用IMB引入固含率εs、碰撞项Ωis完成固-液相互作用的动量交换过程，开发LBMIMB-DEM求解程序。通过层析柱渗流试验，将模拟结果与试验数据相比较，发现4组试验的vout均方根误差均小于10%，而且DKT算例与相关文献结果拟合度良好，认为耦合计算模型与求解方法可以准确地反映悬浮颗粒在含水介质中运移过程。

2）当进口流速为5.0×10-5m/s时，将注入含水介质的悬浮颗粒溶液由球形硅微粉替换为非球形硅微粉后，所选3组断面的平均流速最大降幅分别提升15.0%、11.4%、2.6%。研究表明：球形硅微粉转动惯量Ip较小，易发生滚动、脱离；非球形硅微粉比表面积大，易絮凝成团，沉积在多孔介质表面后难以脱离释放。随着进水流速的降低，颗粒形貌的变化对于含水层水动力场演化过程的影响程度增强。

3）在进口流速为5.0×10-5m/s，向含水介质注入非球形硅微粉悬浮液条件下，当截面距进口距离由50 mm增加至250 mm时，平均流速最大降幅由27.1%上升至56.3%，恢复率由62.5%下降至28.1%。随着含水介质计算区域的增加，颗粒的运动路程增大，颗粒之间以及颗粒与多孔介质的碰撞几率增加，动能损失增大，导致渗流速度降幅增大，恢复率降低。