孔隙尺度下水饱和多孔介质的热-水动力学演化与颗粒重组：一种PSM耦合DEM与DDF-LBM的方法

《International Communications in Heat and Mass Transfer》：Pore-scale thermo-hydrodynamic evolution and particle reorganization in water-saturated porous media: a PSM-coupled DEM and DDF-LBM approach

　　

• 基于PSM耦合的DEM与DDF-LBM模拟含水层储热中的颗粒重组与热-水动力学过程。• 实验与模拟揭示了参数对孔隙结构与热-流演化过程的影响。• 对正交设计数据进行回归分析，量化各参数对出口温度的贡献。

引言

在双碳目标（碳达峰与碳中和）不断推进、能源需求日益增长的背景下，地热资源的开发与利用研究已成为能源领域的优先议题[1]。作为一种清洁、可再生且稳定的能源，地热资源不仅具有能源供应潜力，还可作为热能存储的天然介质[2]。在水饱和多孔介质中，颗粒重组、渗流与传热过程可能影响局部流动阻力与热传输行为，这些现象与地下储热过程的解释密切相关[3]。

在浅层松散含水层的热能储存运行过程中，水力梯度与流体温度的变化以及由此引发的密度驱动对流[4]，加之离子强度[5]与pH值[6]的变化，可诱发孔隙空间内的颗粒迁移、聚集、沉积与再悬浮[7]、[8]。这些过程反过来重塑孔隙结构、降低渗流速度并削弱传热效率[9]。由于饱和多孔介质中固液基质成分的变化可有效调节含水层蓄热系统（ATES）性能[10]、[11]，因此诱发饱和含水层中颗粒迁移行为与孔隙结构演化的机制已成为当前研究的焦点。

柱式物理模拟系统被广泛用于研究多孔介质中的热质传递以及悬浮颗粒的迁移与沉积[12]、[13]、[14]、[15]、[16]。在柱式流动试验中，Agbangla等人[17]确定了颗粒拱桥与沉积的临界通量，超过该值后渗透率急剧下降。利用饱和砂柱，Wang等人[18]发现颗粒尺寸与流速控制微塑料的穿透与滞留。使用循环回注砂柱装置，Wu等人[19]表明较高的离子强度与温度会加速化学沉积并快速降低渗透率。Lee等人[20]在柱式传热实验中监测了流体-固相温度演化过程，揭示局部热非平衡随颗粒尺寸和流速增大而加剧。总体而言，柱式研究能够捕捉颗粒释放、流场与温度场的演化以及孔隙结构变化，但在再现含水层中复杂的迁移与局部堵塞现象方面存在局限[21]。

离散元方法（DEM）将颗粒视为具有独立物理和力学特性的离散体，特别适用于研究颗粒迁移、碰撞、堆积与堵塞问题[22]。DEM在岩土材料颗粒尺度现象中被广泛应用：土体与砂体力学[23]、裂隙岩体演化[24]、含水层细颗粒迁移[25]以及地热储层中的热-水-力（THM）过程[26]。DEM常与计算流体力学（CFD）耦合用于多相流建模，可评估颗粒属性与流体条件对颗粒运动的影响[27]、[28]、[29]。利用经实验验证的DEM-CFD方法，Yang等人[30]表明页岩孔隙中纳米颗粒的迁移与堵塞主要受颗粒尺寸、浓度和孔隙结构控制。Liu等人[31]报告纳米级零价铁-粘土胶体堵塞与穿透主要受流速和颗粒电荷控制。

方法上，DEM-CFD根据流体网格分辨率与颗粒尺寸的关系分为未解析或完全解析两类[32]。完全解析方案可在微观尺度捕捉颗粒尺度流动与局部传热特征，而未解析方案更适合大规模流体-固体耦合[33]。未解析方案的一个关键局限在于无法解析介观尺度下颗粒附近边界层的流动与传热。此外，DEM-CFD模拟的计算效率常受限于DEM与CFD求解器之间同步时间步的制约[34]。

格子玻尔兹曼方法（LBM）是一种介观模拟方法，通过在规则格子上演化分布函数来恢复流体动力学方程，具有空间分辨率高、复杂孔隙几何处理简便以及并行效率优异等优点[35]。DEM与LBM的耦合已发展用于解析复杂多孔结构中的流体-固体相互作用，广泛应用于多相颗粒-流体系统。在孔隙尺度，研究范围从单独使用LBM研究油水液滴运动等现象[36]，到采用DEM-LBM耦合模拟多相流体-固体相互作用[37]并揭示可变形多孔介质中应力依赖的绝对渗透率与相对渗透率[38]。

为更好地捕捉颗粒-壁面相互作用及任意形状颗粒的行为，较新的DEM-LBM框架常纳入浸没边界方法（IBM），通常通过部分饱和法（PSM）实现。PSM能平滑地表征移动固体边界，已成功应用于模拟致密颗粒悬浮液和复杂孔隙几何[39]、[40]。对于热问题，LBM的双分布函数（DDF）格式分别演化动量和温度分布函数，增强了再现耦合流场与温度场的数值稳定性与精度[41]。该格式使Shao等人采用的有限尺寸颗粒双分布函数格子玻尔兹曼方法（DDF-LBM）能够再现孔隙尺度流动与传热，量化了雷诺数、孔隙率和固体温度对速度和温度场的影响[42]。Wang等人提出了一种先进的基于DDF的热LBM，能稳定恢复任意构型下的热-水动力学方程，表现出高精度和广泛适用性[43]。

尽管耦合DEM与DDF-LBM方法可揭示孔隙尺度颗粒迁移与THM耦合机制[44]，但其应用受限于高计算成本、有限数值稳定性以及难以扩展到复杂大规模或多物理系统[45]。相比之下，统计工具如社会科学统计软件包（SPSS）采用方法来直接量化实验或模拟数据中变量之间的相关性和统计显著性。Qin等人利用SPSS相关分析将微CT获取的孔隙结构参数与低阶煤的绝对渗透率关联[46]，而Lin等人利用其曲线估计和非线性回归模块，基于碎屑储层多源数据建立并验证渗透率预测模型[47]。更先进的统计技术有助于刻画多变量关系并估计设计实验或数值数据集中特定响应的相对贡献[48]、[49]。

本研究开发了孔隙尺度耦合框架，将离散元方法（DEM）与双分布函数格子玻尔兹曼方法（DDF-LBM）集成，以部分饱和法（PSM）为基础，确保流体-固体界面动量与能量的离散守恒。利用该框架，研究考察了五个选定参数对模拟水介质中颗粒迁移、渗流行为与传热演化的影响：介质长度（L）、渗流速度（v）、介质孔隙率（ε）、颗粒尺寸（d）和离子强度（I）。采用正交设计组织模拟工况，随后应用多元线性回归（MLR）和方差分析（ANOVA）作为补充方法分析出口温度响应。

各节片段

数学模型

在离散元方法（DEM）框架内，非球形颗粒采用多球模型进行模拟，将其表征为球形子单元的集合体以精确捕捉其动力学行为[50]、[51]。为提高预测精度，纳入了Derjaguin-Landau-Verwey-Overbeek（DLVO）相互作用——包括范德华吸引和电双层排斥——其参数取决于颗粒形态与尺寸[52]。对于流体

实验系统设置

色谱柱渗流与传热实验系统示意图如图4所示。柱体由亚克力材料制成，内截面为30×30 mm，流动长度为150 mm。为约束填充并最小化端部效应，在进出口两端齐平安装了孔径0.125 mm、由0.01 mm厚箔片制成的丝网筛。

进口处温度控制使用GDH0510

计算分析与讨论

基于我们先前团队工作中对十二个固相基质与液相参数的初步筛选——包括七个固相基质参数（介质长度、介质宽度、颗粒尺寸、颗粒数量、孔隙率、颗粒比表面积和颗粒长宽比）和五个液相参数（渗流速度、边界温度、离子强度、pH值和离子价态）——选取了五个相对敏感且具有代表性的因素进行进一步研究

基于多元线性回归的出口温度统计分析

因此，基于数值模拟获得的出口温度数据，采用Pearson相关分析评估所选参数与出口温度Tout之间的线性相关性[78]。建立了多元线性回归（MLR）模型，对所研究参数范围内的出口温度响应提供简洁的统计描述，随后进行方差分析（ANOVA）以估计相对

结论

本研究将色谱柱实验与孔隙尺度数学模型相结合，通过PSM方法耦合DEM与DDF-LBM，研究饱和水介质中的颗粒重组、孔隙堵塞、渗流衰减与传热演化过程。在不同实验条件下，模拟成功再现了观测到的瞬态行为。出口温度和渗流速度的峰值相对误差分别为4.36%和7.93%

符号表

Bf 与固相百分比相关的权重系数
BT 与松弛时间函数相关的权重系数
ci 方向i上的速度矢量（m/s）
Cs 格子声速（m/s）
D 最小表面间距（m）
fi 流分布函数（kg/m3）
fi^eq 流体平衡函数（kg/m3）
Fi 离散力项（N）
Fcn 法向接触力（N）
Fct 切向接触力（N）
Fad,s 球形颗粒间的范德华吸引力（N）
Fad,w 与...相关的范德华吸引力（N）

CRediT作者贡献声明

Jiuchen Ma：原稿撰写、数据整理、概念构思。Yinan Wang：审阅与编辑、资源获取、资金获取、形式分析。Kun Ma：验证、形式分析。Zhuoqian Li：软件、概念构思。Dianxin Zhang：形式分析。Changfeng Wang：监督指导。Xiangcui Lv：监督指导。Rifan Wang：调查研究。

利益冲突声明

作者声明不存在已知的经济利益或个人关系可能影响本文所报告的工作。

致谢

本研究受国家自然科学基金（项目编号41402228）、天津市企业科技特派员项目（项目编号19JCTPJC48100）及中国铁路建设集团限定研究项目课题（项目编号LX24-33C）资助。

Jiuchen Ma | Yinan Wang | Kun Ma | Zhuoqian Li | Dianxin Zhang | Changfeng Wang | Xiangcui Lv | Rifan Wang