基于MRT-LBM与神经网络替代模型的三维多孔介质滑流传热降阶建模

《International Journal of Thermal Sciences》：Reduced-order modeling of slip-flow heat transfer in 3D porous media using MRT-LBM and neural network surrogates

　　摘要：本研究探讨了在连续流和滑移流状态下，由简单立方、体心立方和面心立方球形堆积构成的三维微观多孔介质中的流体流动与热传递现象。为描述流体与固体界面处的速度滑移和温度跃变，研究者开发了一种多松弛时间格子玻尔兹曼模型。在连续流状态下，固体边界附近会出现强烈的流动阻碍及较长的低速区域；而滑移流则能减少壁面引起的滞流现象，并抑制尾部涡旋的形成。热传导方面，连续流中温度上升较早，而滑移流由于近壁对流增强，其峰值温度较低。努塞尔数随雷诺数的增加而上升，但随着克努森数和孔隙率的增加而下降，这是因为有效表面积减小以及温度跃变效应加剧。对于所有晶格结构，摩擦因子均随雷诺数的增加而降低；不过，低孔隙率的体心立方结构由于入口流速较低，所受阻力最大。渗透率则随克努森数和孔隙率增加而上升，其中简单立方和面心立方结构因平均流速较高，其渗透率提升更为显著。为降低计算成本，研究者还训练了一种基于深度学习的替代模型，用以替代MRT-LB模型来预测速度场和温度场，该模型具有高精度且预测误差小，能够快速评估滑移流条件下的热传递情况，并实现微型多孔热系统的高效参数设计。

缩写说明：
下标说明：
BCC：体心立方
bb：边界
FCC：面心立方
f：最邻近流体节点
LBM：格子玻尔兹曼方法
ff：最邻近流体节点的相邻位置
MRT：多松弛时间
gg：气体
SC：简单立方
i：界面

术语说明：
inlet：入口
microscopic velocity：微观流速
jump：跃变
speed of sound：声速
outlet：出口
Cf：摩擦因子
solid：固体
distance ratio：距离比
slip：滑移
D：直径
wall：壁面
E：内能
τ：时间
tangent：切向
density distribution function：密度分布函数
eq：平衡态
F：外力
neq：非平衡态
g：温度分布函数

希腊字母说明：
h：对流热传递系数
α：热扩散率
H：高度
β：热膨胀系数
k：渗透率
γ：比热容比
K：热导率
ρ：密度
l：晶格单元长度
σv：切向动量调节系数
L：特征长度
σT：切向热调节系数
n：法向方向
τv：松弛时间
Nu：努塞尔数
τg：热松弛时间
P：压力
Ω：碰撞算子
Pr：普朗特数
Re：雷诺数
Ra：瑞利数
T：无量纲温度
u：无量纲流速
w：权重系数
x：位置

1. 引言
多孔介质在微观尺度上构成了众多高性能装置的基础，包括微反应器、催化床、MEMS热管理器件、多孔绝缘材料以及紧凑型换热器等。当孔隙尺度接近分子的平均自由程时，连续介质假设的有效性就会受到质疑，从而导致界面非平衡状态的出现。在这种情况下，界面处的速度滑移和/或温度跃变会影响粘性加热与对流在影响压降、壁面热通量乃至渗透率方面的作用，这一现象可通过对特定装置或材料的实验研究得到证实[[1], [2], [3], [4], [5], [6], [7]]。关于速度滑移和温度跃变的理论基础可追溯到动力学理论，大量微观流体装置的研究表明，当通道较窄或弯曲程度较大时，克努森数效应会使得速度滑移和温度跃变变得十分重要[8]——详见图13、图14。要在真实的三维几何环境中对这些现象进行建模，使用传统CFD模拟往往需要极高的计算成本，主要原因在于处理弯曲边界以及稀薄流体中的壁面效应时，需要极高的网格分辨率。格子玻尔兹曼方法则为这一难题提供了有效的解决方案，它不仅能够处理复杂的固体几何结构，还能针对速度滑移和温度跃变边界条件采用动力学边界描述方式。尤其是多重松弛时间方案，相比仅采用单一松弛时间的显式时间模拟方法，能够在模拟微观流体流动及耦合热问题时提高稳定性和精度，因为标准耗散机制可能会扭曲梯度值[[9], [10], [11], [12], [13], [14], [15], [16]]。在研究由球形颗粒构成的多孔介质时，简单立方、体心立方和面心立方排列结构被用作基准，用于研究孔隙率与特征长度对连续流/滑移流状态下的旁通流速、尾部涡旋结构、壁面热传递以及阻力的影响[8,[17], [18], [19], [20], [21], [22], [23]]。近年来，格子玻尔兹曼方法在理解多孔介质中的流动、热传递及多物理场传输现象方面取得了显著进展。张等人[24]通过耦合LBM-DEM框架研究表明，浆料堵塞现象是由颗粒尺寸以及与压力相关的堵塞与疏通转变机制决定的，而非单纯的 Pressure效应。在热应用领域，张等人[25]发现，微槽形吸液结构相较于烧结式和网状吸液结构，能够更有效地促进核态沸腾并抑制膜态沸腾，而周等人[26]则指出，螺旋形多孔支架的结构差异既能提升热传递效率，又能降低流动阻力。利用格子玻尔兹曼方法，人们还系统地研究了多相流动以及润湿性控制下的物质传输现象。张等人[27]发现，在双润湿性多孔介质中，强制吸液过程在低毛细数条件下以毛细作用主导的位移为主，而在较高毛细数条件下则转为由粘性作用控制的指形流动与滞留现象。在电化学系统中，徐等人[28]发现，质子交换膜水电解槽中多孔传输层内的氧渗透行为与孔隙率、孔径以及层厚度密切相关。格子玻尔兹曼方法还被应用于极端高温及多物理场环境中的研究。赵等人[29]建立了用于研究碳基多孔复合材料烧蚀-热解耦合过程的孔隙尺度格子玻尔兹曼模型，结果表明，热解过程中的热量吸收以及多孔炭的阻力效应能够显著提升材料的抗烧蚀性能。徐等人[30]提出了基于REV原理的多物理场格子玻尔兹曼框架，用于研究热化学储热问题，该框架能够将孔隙尺度上的传输现象与宏观反应器性能准确关联起来。在纳米尺度研究中，郑等人[31]结合FIB-SEM重建技术与格子玻尔兹曼方法，证明了在改性酚醛基质中，孔径大小引起的克努森效应是决定其热导率的关键因素。除了上述孔隙尺度格子玻尔兹曼研究之外，还有大量研究致力于建立不同工况下热传递、阻力以及渗透率预测的相关公式。惠特克[32]和桑原等人[33]率先建立了努塞尔数预测公式，随后吴等人[34]、夏等人[35]以及杜等人[36]将这些公式扩展到了更复杂的多孔结构以及类泡沫结构中。在微观层面，已有研究考虑了速度滑移和温度跃变对流体传输特性的影响，例如流体在球形障碍物周围以及微观多孔介质中的流动行为，这些研究结果都表明传统流体动力学理论存在一定偏差[[37], [38], [39]]。同样，埃贡[40]、塔尔马奇[41]以及奥扎希等人[42]提出的常用摩擦因子预测公式，尽管经过后续改进以考虑多孔泡沫及特殊结构的特性[43], [44], [45]，仍主要以连续介质理论为基础，无法很好地描述稀薄流体的特性。在预测多孔介质的渗透率方面，以克林肯伯格[7]的理论为基础，并通过孔隙尺度数值研究进一步完善的模型[46], [47], [48], [49], [50]，已经证明了克努森数和孔隙几何结构对此具有重要影响。而王等人[51]则指出，在燃料电池启动过程中，催化剂层的微观结构对冰的融化以及热管理效果有着决定性作用。在方法论方面，陈等人[52]引入了带有“虚拟流体”边界的分数阶多松弛时间格子玻尔兹曼模型，该模型有助于更好地模拟流体在多孔介质中的流动与热传递过程中的记忆效应。总体而言，这些研究都体现了格子玻尔兹曼方法在分析孔隙尺度传输现象方面的优势，同时也凸显出为了实现大规模参数分析，发展高效替代建模方法的必要性。近年来，孔隙尺度建模、格子玻尔兹曼方法以及机器学习技术的进步，极大地提升了人们对多孔介质和膜材料中热、质及流体传输现象的理解。唐等人[53], [54]通过对地热CO2-蒸汽分离过程的优化研究，指出逆流流动、大孔隙梯度以及高孔隙率是提升吸收效率的关键因素。李等人[55]发现，基于控制参数的MLP模型在预测各向异性多孔介质的方向性热导率方面，性能优于基于CNN的图像预测方法。苏艾等人[56]研究表明，在N型换热器中，混合纳米流体の性能会受到纳米颗粒占比、磁场强度、跨扩散效应以及孔隙率的影响。赖等人[57]发现，石墨烯基底的结构会直接影响复合材料中的相变热传递行为，其中二维有序结构具有更优异的垂直方向热导率。索达加尔-阿巴德和洛伊默[58]研究了存在冷凝/蒸发现象时的受限纳米流体流动，重点分析了液面曲率与壁面之间的相互作用。徐等人[59]在质子交换膜燃料电池的催化剂层中考虑了纳米尺度下的粘度衰减效应以及滑移现象，从而优化了氧气的传输效率以及反应动力学过程。洪等人[60]发现，经过PTFE改性的分层表面可以通过调控润湿性以及气泡动态来提升池沸腾效应。林等人[61]利用在格子玻尔兹曼数据基础上训练的3D CNN-U-Net模型，实现了对金属泡沫流动速度和温度的高精度预测，同时计算时间仅需传统方法的95%至99%。李等人[62]将U-net模型与格子玻尔兹曼方法相结合，用于模拟裂隙页岩中的油水两相流动，成功揭示了流体堵塞以及突然释放的现象。吴等人[63]则将CASSA-DELM方法与改进后的格子玻尔兹曼方法结合起来，用于研究三维食品挤出过程，建议对于单组分物料可使用活塞式挤出机，而对于多组分乳化物料则更适合使用螺杆式挤出机。第二个挑战在于，对于涉及固定范围内雷诺数、克努森数、孔隙率以及边界条件的参数研究而言，全规模的多松弛时间格子玻尔兹曼模拟在计算上几乎不可行，而这正是开展设计优化研究所必需的。因此，降阶建模以及替代建模已成为当前热流体研究中的重要手段。通过在高分辨率解决方案的基础上训练替代神经网络模型，就可以以较低的计算成本近似求解滑移流状态下的速度场和温度场，进而快速估算局部努塞尔数、摩擦因子以及渗透率修正值等相关参数[[64], [65], [66], [67], [68], [69], [70]]。基于物理原理的替代建模方法以及学习运算符，通过融入基于控制方程的结构或者直接学习解算子而非特定回归问题，进一步提升了模型的通用性[61,64,67]。基于以上考虑，将精确的多松弛时间格子玻尔兹曼模拟输入与经过充分验证的深度神经网络替代模型相结合，有望为高效模拟三维多孔介质中的滑移流传输现象提供可行路径。尽管已有大量基于格子玻尔兹曼方法的研究聚焦于孔隙尺度现象，但现有文献大多仅从连续介质视角或简化后的滑移流模型出发，来研究速度滑移和温度跃变现象，较少通过格子结构来分析热力学与流体力学现象之间的耦合关系。针对不同热力学与流体力学条件下，简单立方、体心立方以及面心立方晶格结构中球形颗粒的排列方式，目前进行的对比研究也非常少。此外，对于参数化研究而言，多次应用多松弛时间格子玻尔兹曼方法进行计算效率过低，而针对稀薄多孔材料中的孔隙尺度速度场和温度场进行预测的替代模型也尚不完善。表1汇总了2024年至2026年间关于孔隙尺度格子玻尔兹曼模拟、稀薄气体流动以及机器学习辅助模拟方面的最新研究情况。这份总结有助于明确当前研究的局限性，同时也凸显出本研究通过采用多松弛时间格子玻尔兹曼模型、考虑速度滑移与温度跃变效应、对比不同SC/BCC/FCC几何结构以及运用ANN替代模型所带来的独特价值。表1. 最近关于孔隙尺度热传递与稀薄流动的研究与本研究的对比。|参考文献编号|年份|数值方法|几何结构|速度滑移与温度跃变|热传递|人工智能/机器学习应用|局限性|Fan等人，[71]|2025|格子玻尔兹曼方法|随机多孔介质|?|?|?|侧重于反应传输现象，未涉及稀薄流体的滑移流现象。|Fan等人，[72]|2025|格子玻尔兹曼方法|随机多孔介质|?|?|?|未考虑稀薄效应，也未进行几何结构对比。|Ranjbarzadeh与Sappa|[73]|2025|综述|一般多孔介质|部分|?|部分|仅为综述文章，没有新的孔隙尺度模拟研究。|Zhao等人，[74]|2023|格子玻尔兹曼方法+概率神经网络|真实多孔介质|?|?|?|仅研究流动现象，未涉及热传递或滑移效应。|Wang等人，[75]|2024|孔隙尺度数值建模|复杂多孔介质|?|?|?|进行了传输建模，但未对比速度滑移与温度跃变现象。|李俊阳等人，[76]|2024|有限元方法|多孔介质|?|?|?|仅进行连续介质尺度分析。|Kirch, Alexsandro等人，[77]|2026|综述|一般多孔介质|?|?|?|指出了人工智能的应用前景，但未构建整合多松弛时间格子玻尔兹曼模型的框架。|余金波等人，[78]|2025|多松弛时间格子玻尔兹曼方法|多孔介质|部分|?|?|侧重于反应传输现象研究，未涉及滑移流热传递问题。|本研究|2026|多松弛时间格子玻尔兹曼方法+人工神经网络|SC/BCC/FCC球形孔隙几何结构|?|?|?|同时研究了滑移流、温度跃变、几何结构对比以及人工神经网络替代模型效应。|为克服上述缺陷，本研究提出了一种基于物理原理的降阶建模方法，该方法将多松弛时间格子玻尔兹曼方法与基于数据的代理模型相结合，用于求解滑移流问题。该建模方法既能够实现耦合流动现象的可靠模拟，又能在计算效率上有所提升，便于在不同稀薄条件下分析各种几何结构的热传递与压降特性。2. 用于气体流动与热传递的3D格子玻尔兹曼模型2.1. 用于气体流动的LB模型为模拟三维空间中的气体流动，本研究采用了D3Q19模型，即三维结构加上19个速度分量（详见图1）。气体粒子可以在立方晶格上移动，其微观流速用ci表示。密度分布函数fi(x,t)可通过以下公式确定：(1)fi′?fi=Ωi(f)其中fi′=fi(x+ciΔt,t+Δt)，fi=fi(x,t)，Δt表示时间增量。此外，还会应用多松弛时间碰撞算子Ωi(f)，该算子可根据松弛矩阵S和变换矩阵M表示为：(2)Ωi(f)=?M?1SM(fi?fieq)下载：下载高分辨率图像（339KB）下载：下载完整尺寸图像图1. 格子玻尔兹曼方法所用的D3Q19模型。变换矩阵通过公式=f=M将动量空间与密度分布函数联系起来。因此，方程(1)可改写为：(3)fi′?fi=?M?1S(m?meq)变换矩阵和松弛矩阵的具体表达式如下：(4)M=[1111111111111111111-30-11-11-11-11-11-1188888888888812-4-4-4-4-4-411111111111101-100001-11-11-11-100000-4400001-11-11-11-100000001-10011-1-100001-11-1000-440011-1-100001-11-1000001-1000011-1-111-1-100000-44000011-1-111-1-1022-1-1-1-111111111-2-2-2-20-4-4222211111111-2-2-2-200011-1-11111-1-1-1-10000000-2-2221111-1-1-1-1000000000001-1-11000000000000000000000001-1-110000000000011-1-1000000000001-11-1-11-1100000000000-1-11100001-11-10000000000011-1-1-1-111](5)S=diag(0,se,sε,0,sq,0,sq,0,sq,sυ,sπ,sυ,sπ,sυ,sυ,sτ,st,st,st)其中se=1/τe，sυ=1/τυ。在LB模型中，松弛时间τυ由微观尺度关系式给出：(6) τυ=1/2+λΔx/6π。剪切粘度和体积粘度与松弛率τυ和τe的关系为：(7) υ=cs2(τυ?1/2)δt，(8) ζ=cs2(τe?1/2)δt。此时，平衡分布关系可表示为：(9) fieq=wiρ(1+ci·ucs2+(ci·u)2/2cs??u·u2/cs2)，其中局部声速cs为1/3。D3Q19模型的权重wi和晶格速度ci如表2所示。

表2. D3Q19模型的权重函数与晶格速度。
i 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18
wi 13 118 118 118 118 118 118 118 136 136 136 136 136 136 136 136 136 136 136
cix 0 1 ?100 0 1 ?11 ?100 1 ?11 ?100
ciy 0 0 0 1 ?100 1 ?100 1 ?1 100 1 ?1
ciz 0 0 0 0 0 1 ?100 1 ?11 ?100 ?11 ?11

此时，宏观尺度下的密度和速度变量可通过以下公式计算：(10) ρ=∑ifi，(11) ρu=∑ifici2。

2.2 用于模拟速度滑移的边界方案
球形壁面通过非平衡外推边界条件进行模拟，该条件下固体边界处的分布函数被分解为平衡分布函数与非平衡分布函数，即：(12) fi,s=fi,seq+fi,sneq。固体壁面的平衡分布与方程(9)相同，其中靠近固体壁面的气体速度us是通过外推方法确定的，参见[12]。(13) us={us1d≥3/4dus1+(1?d)us2d<3/4}，其中距离比d=|xf?xint|/|xf?xs|，(14) us1=1duint+d?1du(xf)，us2=2d+1uint+d?1d+1u(xff)。类似地，非平衡分布函数fi,sneq的外推方法为：(15) fi,sneq={dfi,fneqd≥3/4dfi,fneq+(1?q)fi,ffneqd<3/4}。为在微观尺度上精确描述速度，需考虑球体气固界面处的滑移速度usl，其定义如下：(16) usl=ug?uw。滑移速度的定义方式与麦克斯韦一级滑移边界条件类似：(17) usl=2?σvσvλ(?uτ/?n?uτr)，其中n为沿固体壁面的法向矢量，τ为切向矢量，完全扩散壁面的动量适应系数σv取值为1。随后，在固体边界应用反外推方案来计算速度梯度。在垂直于壁面的矢量方向上考虑两个虚拟点n′和n″，这两个点均位于流体区域内，且满足δ=|n″?n′|=|n′?xint|，δ=1.5Δx。然后通过三线性插值法确定宏观气体速度u′和u″。方程(16)中的界面处气体速度可通过[79]进行转换。(18) ug=uw+[λ(2?σv)σv+λ(2?σv)(1.5δ+1Rw)(4u′?u″?3uw2δ?uwr)·τ]/τ，其中ug为壁面处的气体速度，uw为壁面速度，λ为分子平均自由程，σv为切向动量适应系数，δ为晶格间距，Rw为固体壁面（球面）的局部曲率半径，τ为壁面的单位切向矢量，u′和u″分别表示紧邻壁面的第一和第二流体节点处的速度。

2.3 用于热传递的LB模型
为利用LB模型模拟热传递，需使用热分布函数gi。分布函数的变化可表示为：(19) gi′?gi=ΩiT(g)，其中SRT碰撞算子用于定义热碰撞算子Ωi(g)：(20) ΩiT(g)=?1/τg[gi?gieq]。热平衡分布函数gieq和热松弛时间τg分别为：(21) gieq=wiT(1+ci·ucs2+(ci·u)2/2cs??u·u2/cs2)，(22) τg=1/2+υ/Prcs2Δt。此时，可通过将相邻节点处的所有热分布函数相加来计算晶格点的宏观温度：(23) T=∑igi。

2.4 用于模拟温度跃变的边界方案
球形边界处的热分布函数可与密度分布函数类似地表示为：(24) gi,s=gi,seq+gi,sneq，其中固体球形边界处的热平衡分布函数gi,seq是通过方程(21)确定的，而靠近球形节点处的速度则通过方程(13)确定。非平衡热分布函数则通过外推方法计算：(25) gi,sneq={dfi,fneqd≥3/4dfi,fneq+(1?d)gi,ffneqd<3/4}。对于该区域内的热传递，温度跃变的一级麦克斯韦型关联式可表示为：(26) Tj=Tg?Tw=2?σTσT2γγ+1λPr(?T/?n)，其中γ为比热容比，温度适应系数σT取值为1，普朗特数Pr为0.71。与速度计算类似，节点n′和n″处沿法向矢量n的温度也是通过三线性插值法确定的。此时界面处的气体温度可通过[79]表示为：(27) Tg=Tw+γλ(2?σT)δσTPr(γ+1)+3γλ(2?σT)(4T′?T″?3Tw)。

2.5 热流域的构建
利用改进后的D3Q19热玻尔兹曼模型，研究三维域内气体流动的速度滑移和温度跃变现象。构建了三种不同球形排列的三维孔隙尺度域，即简单立方结构、体心立方结构和面心立方结构。所研究的域中心为球形排列，入口和出口到最外层球心的距离始终等于球的直径D。所有轴方向上球心之间的距离为50lu。气体从温度较低的入口（Tin=0.0）沿正x轴方向流过该域，而球体所在区域的温度较高（Tw=1.0）。域的入口和出口表面采用非平衡外推边界条件进行建模，侧墙则应用对称边界条件。在球形曲面处，结合非平衡外推边界条件和复合外推方法，以捕捉球体处的速度滑移和温度跃变效应。通过改变Kn数，研究孔隙尺度域在无滑移和有滑移状态下的流动与热传递情况。通过改变球的直径，使所研究域的孔隙率在0.6到0.8之间变化。通过改变Re值从0.5到5.0以及相邻两个球体之间的最小距离，来评估不同流动参数对努塞尔数、摩擦因子系数和渗透率的影响，这些参数可通过以下公式确定：(28) Nu=D(?T/?n)w/Tw?Tbulk，(29) Cf=μ(?u/?n)wall/0.5ρuav2，(30) ki=uavυΔx/23NxΔρ。

3 结果与讨论
本节对简单立方、体心立方和面心立方孔隙尺度域内的流动与热传递行为进行了全面分析。利用MRT-LB模型，研究了在不同雷诺数和克努森数下速度分布、流线形态、温度演变、努塞尔数、摩擦因子以及渗透率修正因子的变化情况（见图2）。
下载：下载高分辨率图像（573KB）
下载：下载全尺寸图像
图2. 示意图，分别展示了(a)简单立方、(b)体心立方以及(c)面心立方球形排列结构。

3.1 速度分布
图3显示了在不同雷诺数和克努森数条件下三维孔隙尺度域内的速度场，反映了孔隙几何结构与稀薄效应对流场的影响。在连续流状态下，由于粘性作用带来的高阻力，无滑移边界条件会在球体表面附近形成较大的低速区域。在这三种几何结构中，由于特征长度较小，简单立方和面心立方结构具有较高的速度梯度和峰值速度，而体心立方结构由于特征长度较大，导致入口速度较低，但流速会随之加快。
下载：下载高分辨率图像（2MB）
下载：下载全尺寸图像
图3. 三维孔隙尺度域内的速度变化情况。
随着克努森数的增加，滑移流状态的重要性及其对速度场的影响日益显著。滑移边界条件降低了壁面摩擦力，使得气体在球体表面附近更易于流动。因此，低速区域变小，流动更加均匀，壁面附近的速度梯度增大，从而提升了流体在孔隙尺度域内的输送能力。在简单立方结构中，由于流道较窄，这种效应更为明显。

3.2 流线形态
图4展示了三种球形排列结构下的流线形态。这些流线形态有助于了解孔隙几何结构与滑移条件如何影响流动的循环现象。
下载：下载高分辨率图像（2MB）
下载：下载全尺寸图像
图4. 三维孔隙尺度域内的流线变化情况。
在连续流状态下，所有排列结构中，尾随球体后方都会形成涡旋结构。这些涡旋是由于流动分离和剪切层相互作用产生的，尤其是在流动遇到突然的几何收缩区域时更为明显。由于间距较小，简单立方和面心立方结构形成的循环区域更为紧凑，从而导致流动方向的突变更为剧烈。相比之下，体心立方结构形成的涡旋虽然规模较大，但强度较弱，其原因在于：①空隙空间更大；②入口速度较低；③球体周围的加速作用较为平缓。在滑移流状态下，涡旋的强度会显著减弱或完全消失。球体表面的滑移速度降低了剪切力，抑制了流动分离，使得流线更加平滑。尽管如此，由于多个尾流与转向路径之间的相互作用，整体三维流动依然较为复杂。

3.3 温度分布
图5展示了不同雷诺数和克努森数条件下的温度等值线，反映了孔隙的稀薄效应与几何形状对热传递过程的影响。当流体仍处于连续流状态时，无滑移边界条件的存在会增加球体表面的热阻，使得流体温度更快达到最大值。由于孔隙间距较小，简单立方结构会导致温度上升较快，而体心立方结构由于特征长度较大、流入域的流体速度较低，且对流作用较弱，因此其温度达到最大值的程度相对较低。
下载：下载高分辨率图像（2MB）
下载：下载全尺寸图像
图5. 三维孔隙尺度域内的归一化温度变化情况。
随着克努森数的增加，滑移与温度跃变效应对热场的影响愈发突出。滑移边界条件降低了热阻，而温度跃变效应则降低了壁面热通量。这使得峰值温度降低，温度场的形成时间有所延迟，且球体表面附近的温度梯度变得更加陡峭。在较高的克努森数条件下，最大温度的下降更为明显。虽然连续流状态能更快达到热平衡，但滑移流状态能使较冷的边界层在更长的距离上保持不变，这凸显了对流效应与温度跃变效应之间的相互作用。

3.4 努塞尔数分析
为准确理解努塞尔数随简单立方、体心立方和面心立方结构变化的趋势，需要认识到速度滑移与温度跃变对对流的影响方向是相反的。在克努森数较低时，热量传递主要通过球体表面的传导以及狭窄流道中的对流实现。随着克努森数的增加，速度滑移降低了壁面剪切力，有利于热传递，但温度跃变效应又会对热传递产生不利影响。不过，由于温度跃变效应占主导地位，整体热传递效率仍会出现下降。图6展示了三种球形排列结构下努塞尔数随雷诺数、克努森数以及孔隙率的变化情况。随着雷诺数的增加，由于流速升高，对流作用增强，努塞尔数也会随之上升，其中简单立方结构由于特征长度较小且速度梯度较高，因此努塞尔数达到最大值。然而，随着克努森数的增加，尽管滑移速度有所上升，但由于温度跃变效应降低了从固体壁面传递的热量，努塞尔数又会下降。此外，随着孔隙率的增加，努塞尔数也会持续下降。
下载：下载高分辨率图像（384KB）
下载：下载全尺寸图像
图6. 努塞尔数随(a)雷诺数以及(b)克努森数的变化情况。
简单立方和面心立方结构由于结构更为紧密，更依赖于表面积驱动的热交换，因此努塞尔数的下降更为明显。而体心立方结构由于流道较大，壁面间的相互作用较弱，因此孔隙率增加时努塞尔数的下降幅度相对较小。在滑移流状态下，温度跃变效应进一步降低了壁面热通量，使得图6中由孔隙率增加导致的努塞尔数下降更为显著。

3.5 与努塞尔数关联式的对比
为评估所研究模型的预测能力，将三维情况下数值计算得到的努塞尔数与方程(31)中所给出的通用关联式进行对比。如方程(31)所述，该方程给出了对流热传递过程与描述流体流动的物理参数——即雷诺数、克努森数和孔隙率——之间的函数关系。方程(31)可作为评估多孔材料中对流热传递过程的参考标准。如前文所述，方程(31)是基于理想化几何假设和量纲分析推导而来的，其中流动惯性与稀薄效应的作用是通过经验性比例律来描述的。然而，这样的推导并未考虑多孔材料内部三维空间中复杂的传输机制。本节将方程(31)与当前的数值计算结果进行对比。
(31) Nu=A+0.5ε+0.15εBReC?DReEKnε
方程(31)为多孔介质中的对流热传递提供了半经验性的比例关系，将努塞尔数与雷诺数、克努森数以及孔隙率联系起来。图7展示了在无滑移流和滑移流情况下，三维模拟得到的努塞尔数与方程(31)所描述的相关式之间的比较结果。从图中可以明显看出，三维情况下预测的努塞尔数要高得多，且更依赖于雷诺数，而使用相关式得到的结果则没有这种趋势。出现这种差异主要是因为相关式是基于二维分析得出的，未能考虑三维中存在的各种复杂性。例如对于球形颗粒而言，曲率、二次流以及三维涡流都会增强对流作用，从而导致更高的努塞尔数。表6. SC结构数据（计算时间：7244.61秒）。
数据R2RMSESEMSE
速度（m/s）0.98690.000010.00001
温度0.98960.02820.0008

表7. BCC结构数据（计算时间：7372.71秒）。
数据R2RMSESEMSE
速度（m/s）0.97490.000010.00001
温度0.98790.03030.0009

表8.对于FCC数据（经过时间：5749.63秒）。DATAR2RMSEMSE速度（米/秒）：0.9401、0.0001、0.00001；温度：0.9838、0.0337、0.0011。下载：高分辨率图像（1MB），下载：全尺寸图像。图15：对比MRT-LBM与深度学习模型对(a)简单立方、(b)体心立方以及(c)面心立方多孔结构中温度和速度的预测结果的差异图。下载：高分辨率图像（970KB），下载：全尺寸图像。图16：(a)简单立方、(b)体心立方以及(c)面心立方结构的训练集与验证集的损失收敛情况。需要指出的是，该替代模型是基于在固定流动条件（Kn=0、孔隙率=0.7、Re=0.5）下为三种多孔几何结构（简单立方、体心立方和面心立方）生成的MRT-LBM数据集训练而来的。因此，所构建的模型主要作为所研究参数空间内速度场和温度场的简化表示。尽管该模型对未见过的数据样本具有出色的预测能力，但其在外推至训练数据所涵盖范围之外的情况时的表现尚未得到评估。将这一框架扩展到包含多种雷诺数、克努森数、孔隙率以及几何特征，是未来开发能够全面覆盖设计空间的通用替代模型的重要方向。虽然所提出的前馈神经网络具有优异的预测性能，但其架构属于数据驱动的近似方法，并未明确体现守恒定律、滑移边界条件或能量方程。因此，模型预测结果受到来自MRT-LBM模拟的训练数据范围和质量的限制。更为先进的机器学习框架，如卷积神经网络、物理信息神经网络、深度算子网络以及傅里叶神经算子，或许能为高度非线性的传输现象提供更好的物理一致性和泛化能力。研究这类架构在微尺度多孔介质传热问题中的应用，也是未来工作的重要方向。4. 结论本研究利用多松弛时间格子玻尔兹曼模型，研究了由简单立方、体心立方和面心立方排列的微球体构成的三维微多孔介质中的气体动力学行为及热传导现象。本研究成功引入了基于非平衡状态的界面滑移边界方案以及反外推方法，用以模拟连续态和滑移态下的情况。研究结果表明，微多孔介质中的现象会受到稀薄效应的显著影响。而在连续态下，由于球体表层附近的粘性阻力，会出现低速区域，同时整个介质区域的温度会升高，且跟随在后的球体尾流中会形成涡旋。这些现象在密集排列的球体阵列中最为明显。随着滑移作用的重要性增加，靠近壁面的停滞区域会减小，跟随球体产生的涡旋会减弱或消失，此外，由于壁面对流的增强以及温度跃变带来的热阻，最高温度也会下降。在三种阵列几何结构中，于连续态下，体心立方阵列由于其特征长度更大、入口速度更低，因此会产生低速区域并产生拖曳效应。参数分析显示，随着雷诺数的增加，努塞尔数的增幅大于因克努森数和孔隙率增加而导致的降幅。虽然速度的滑移现象会随着对流传输而加剧，但在稀薄状态下，温度跃变的作用则会起相反的作用，从而降低高效传热的效率。对于各类晶格结构，摩擦因子都会随着雷诺数和相对稀薄程度的增加而减小，不过在孔隙率较低的情况下，体心立方结构的摩擦因子会更高。渗透率修正因子则随着克努森数和孔隙率的增加而上升，其中简单立方和面心立方结构由于平均速度较高而更具传热效率，而体心立方结构在相对稀薄程度较高时则更为敏感。为降低使用MRT-LB模型进行孔隙尺度模拟的高昂计算成本，本研究构建了一种深度学习替代模型，可通过笛卡尔坐标直接计算速度和温度分布图。该替代模型在匹配所有测试晶格结构的MRT-LB解方面表现出极高的精度，而且能够在大幅降低计算成本的同时，还原出壁面梯度、尾流以及温度分布图。本研究提出的物理模型与机器学习相结合的方法，为处于滑移流状态下的微多孔热系统提供了一种高效的简化模型。作者贡献说明：Ammar Tariq：概念构建、正式分析、研究设计、方法论、初稿撰写；Umar Farooq：数据整理、软件应用、验证、可视化、修改润色；Tauqir Muhammad：正式分析、可视化、初稿撰写、修改润色；Naveed Husnain：修改润色。Ammar Tariq|Umar Farooq|Tauqir Muhammad|Naveed Husnain