多孔介质内嵌热源引起的对流换热格子Boltzmann模拟

刘子阳，陈姿颖，王存海

(北京科技大学能源与环境工程学院，北京 100083)

多孔介质是由固体骨架和由骨架分隔成大量密集的微小空隙所构成的物质，是地球生物圈的基石之一.土壤、动物的微细血管网络、植物的根茎枝叶等都是天然多孔体，陶瓷、砖瓦、玻璃纤维等人造多孔介质材料也参与构成了人类社会的方方面面.多孔介质内的对流换热过程普遍存在于先进材料制造[1]、电子器件冷却[2]，以及生物医学[3]等相关技术与工程应用中，因此多孔介质内对流换热的研究具有重要意义.然后，由于多孔介质结构的复杂性，多孔介质内对流换热的实验研究难度较大且费时费力，因此数值模拟方法在多孔介质内对流换热问题的研究取得了广泛关注：Bayta等[4]采用一种隐式的有限差分法(Finite Difference Method，FDM)求解了达西方程和能量方程，数值模拟了梯形多孔介质内的自然对流过程并研究了瑞利数对温度场和流场的影响；Chu等[5]采用直接数值模拟(Direct Numerical Simulation，DNS)方法研究了多孔介质中的对流换热问题并分析了雷诺数(Reynolds number，Re)的影响；Habbachi等[6]应用有限体积法(Finite Volume Method，FVM)分析了填充多孔介质的三维方腔中的自然对流，发现固体与流体的热导率是影响对流换热的重要因素.

近年来，从格子气自动机[7]发展而来的格子-Boltzmann方法(Lattice Boltzmann Method，LBM)[8-14]是一种利用介观动力学模型来模拟复杂宏观输运现象的数值方法.该方法具有清晰的物理意义并且易于并行运算，因此广受关注并被成功地应用于多孔介质对流换热问题[15-20]的研究.目前，针对多孔介质中对流换热的研究多为壁面加热的情况[21-23]，但是在电子元件冷却器[24]、化学反应器[25]和热交换器[26]等诸多工程设施往往同时涉及多孔介质和内部热源.因此，研究多孔介质内嵌热源引起的对流换热问题有助于对相关工程设备的优化设计.然而，由于底部加热带来的热流分岔等非线性问题[27]以及内部热壁面之间的相互影响，多孔介质内嵌热源引起的对流换热问题[28 ，29]的求解和分析仍需进一步研究.

本文针对底部内嵌热源的多孔介质模型，基于Guo和Zhao[30]发展的求解流场和温度场的分布函数，充分考虑内嵌热源和方腔壁面的边界条件，建立了多孔介质内嵌热源引起的对流换热的格子-Boltzmann模型.通过实验结果对比，验证了本文数值模型和代码的正确性，进一步分析了瑞利数(Rayleigh number，Ra)、达西数(Darcy number，Da)和介质孔隙率(Medium Porosity，ε)等特征参数对多孔介质内流场和温度场的影响.

1 物理模型与数值方法

1.1 物理模型

考虑如图1所示的底部存在内嵌热源的多孔介质.方腔的尺寸为L×H，左右壁面为T=Tc的定温壁面，上下壁面为绝热壁面.内嵌热源位于底边中心位置，尺寸为l×h，温度恒定为T=Th(Th>Tc).考虑到方腔内的流体流速较小，该问题可看作是不可压缩的层流问题.

图1 物理模型

1.2 控制方程

在表征单元体系(Representative Elementary Volume，REV)尺度下，多孔介质对流换热问题的宏观控制方程表示为[31-33]

∇·u=0，

(1)

(2)

(3)

公式中：u为流体速度；T为温度；p为压力；ε为多孔介质的孔隙率；ρ为流体密度；υe为有效运动粘度；αe为有效热扩散系数；σ为热容比，此处为1；F为多孔介质和其他外力产生的总作用力，可表示为[32]

(4)

G=gβ(T-T0)j，

(5)

公式中：G为浮力，υ为运动粘度，g为重力加速度，β表示热膨胀系数，T0为基准温度，j为垂直方向的单位矢量；Fε为几何形状因子，K为渗透率，根据Ergun的经验关系[34]可以表示为[35]

(6)

(7)

对于上述控制方程所描述的物理问题，可以由无量纲参数表征，即达西数Da、瑞利数Ra、普朗特数Pr、努赛尔数Nu与粘度比J为

(8)

(9)

同时，无量纲速度和温度可以表示为

(10)

1.3 格子-Boltzmann模型

本文采用多分布函数的热格子-Boltzmann模型，边界条件采用非平衡外推格式.流场部分使用含外力项的D2Q9格子的单松弛模型处理，温度场部分使用D2Q5格子的单松弛模型处理，其演化方程与对应的平衡分布函数为[30]

(11)

(12)

(13)

(14)

公式中：Fi为外部体力项，可表示为

(15)

宏观变量可表示为

(16)

(17)

公式中：引入的中间速度v与修正参数c0和c1为

(18)

(19)

1.4 边界条件

对于流场，采用非平衡外推边界，即将边界节点处的分解为平衡部分与非平衡部分，边界节点的离散分布函数fi(xs，t)表示为[36]

(20)

公式中：xf=xs+eiδt表示和固体边界节点xs相邻的流体节点xf.

温度场的边界条表示为[37]

左边界

(21)

右边界

(22)

公式中：对于外部方腔的左右边界，无量纲温度取θ=θc；对于内嵌热源边界取θ=θh.内嵌热源顶壁的边界条件为

(23)

封闭方腔的上下边界是绝热的，故可以根据下式求取未知的分布函数.

上边界

g4(x，ym)=g4(x，ym-1)，

(24)

下边界

g2(x，y0)=g2(x,y1)，

(25)

公式中：下标m和0分别表示方腔顶壁和底壁在y方向的节点索引.

2 模型验证

2.1 网格无关性验证

首先设置孔隙率和达西数分别为ε=0.999 9和Da=108，将本文发展的多孔介质模型应用于底边局部加热自然对流过程[38]的求解：加热尺寸为d=0.8 L，模拟参数设置为Pr=0.71、Ra=1.836×105.计算所得的热壁面平均努赛尔数Nuavg如图2所示，模拟选用的网格类型分别有Nx×Ny=40×40，60×60，80×80，100×100，120×120，140×140和160×160.从图中可以看出，当网格数量达到120×120时，努赛尔数的变化逐渐趋于平稳.为确保收敛，Nx×Ny=300×300时的工况也被计算，其数值结果为Nuavg=7.501，以该结果为基准，采用120×120网格所得结果的相对误差小于1.0%.综合考虑计算效率和精度，本文接下来的数值模拟所选取的网格划分均为120×120.

图2 网格无关性验证

2.2 正确性验证

首先将采用本文模型得到的模拟结果与文献[38]的实验图像进行比对.由图3可知，本文采用的LB模型精确地复现了该方腔内的温度分布特征，验证了本文模型的准确性.

图3 底部局部加热方腔内温度分布

为定量分析模型的数值可靠性，本文模拟了不同瑞利数时底边加热的自然对流，分析高温壁面上的平均努赛尔数随瑞利数的变化.表1中列出了在加热尺寸分别为d=0.2L和d=0.8L时Ra=103、104、105和106的数值结果，与文献[38]中的数据符合程度良好.以上对比结果表明本文采用的LB模型能够准确地模拟自然对流换热现象.

表1 本文LB得到的加热壁面上的Nuavg和实验结果[38]对比

3 结果与讨论

3.1 瑞利数变化对换热的影响

瑞利数是与浮升力直接相关的无量纲准则数.为了研究瑞利数对具有底部矩形热源的多孔介质方腔内部自然对流的影响，本节将模拟瑞利数为Ra=101～107时的温度场至稳态.相关模拟参数设置为Da=10-3，ε=0.4，Pr=0.71，l=0.2L，h=0.2H，部分计算结果如图4所示.

图4 不同瑞利数条件下多孔介质内的等温线(左)和流线(右)

由图4(左)所示的温度分布结果可知，热量从热源壁面逐渐扩散到多孔介质的内部空间，随着Ra的增加，方腔顶部温度梯度更大.观察图4(a)可以发现，当Ra较低时，等温线在热源壁面附近密集，此时的换热性质主要为热传导；随着Ra的增加，热传导的作用减弱，热量逐渐从热源壁面附近向多孔方腔的顶部扩散：在Ra=106时，对流换热效果十分显著.图4(右)为对应图4(左)工况下的流线图.根据流线分布可以发现，多孔方腔中形成了沿x=0.5L对称的反向旋涡：当Ra=104时，旋涡的中心位置较低，同时流线主要在方腔底部密集；当Ra增大时，旋涡中心位置随之升高，在方腔顶部能够观察到更为密集的流线分布：多孔介质内的整体流动情况得到明显改善.

不同Ra下的垂直速度分布曲线与Nuavg变化曲线如图5所示.根据图5(a)展示的多孔介质腔体在x=0～0.5L、y=0.5H处的垂直速度分布可知：在低Ra时，速度几乎为零，流动处于停滞状态；而Ra=105～107时，由于浮升驱动力的增强，速度曲线产生明显的变化：最大速度出现在热源上方位置，最小速度则靠近低温侧壁，流体的流动行为较为显著.由此可见，增大Ra能够有效促进多孔介质中流体的运动，但是需要达到临界值：在本文的研究工况和数值范围内，该临界值为Ra=105.从图5(b)中可以发现，热源各壁面的Nuavg与Ra呈现正相关性，即Ra的增大能够提升多孔介质的整体对流换热强度.但考虑到高温热源存在多个热壁面，且壁面间的热交换可能相互影响，故还需要对热源的各个壁面进行对比分析.

图5 不同瑞利数条件下：(a)介质中心线y/H=0.5上流体垂直速度分布；(b) 热源壁面上的Nuavg

当Ra较低，例如Ra=101～104时，沿各壁面处的Nuavg均低于1.0，说明多孔介质内的换热形式以热传导为主.同时可以发现，沿热源侧壁(含热源左侧壁面与右侧壁面)处的Nuavg始终低于沿顶壁处的Nuavg.随着Ra的提升，沿热源侧壁处的Nuavg呈现增长趋势，沿顶壁处的Nuavg则逐渐降低，直到Ra=104时两者几乎相等.这意味着，在换热性质主要为热传导时，顶壁的换热情况并未随着浮升力的增强得到改善，反而是随侧壁换热能力的提升而出现减弱.由此可见，热源顶壁处的换热实际上受到了侧壁换热的影响：随着浮升驱动力的增大，热源侧壁向上方传递的热量增多.但由于热传导的影响，热源顶壁附近热边界层较厚，热流体处于流动滞止区，顶壁处的温度梯度低，从而使Nuavg的数值下降.

当Ra继续增加，Nuavg出现明显提升，同时沿热源侧壁处的Nuavg开始高于顶壁，这在Ra=107时尤为显著.由于Nuavg远高于1.0，热对流取代热传导成为了多孔介质内换热的主要形式.观察图像还可以发现，在临界值Ra=105处，沿顶壁处的Nuavg继续降低，而此时沿侧壁处的换热性质开始向对流转变.同时在Ra>105时，沿侧壁处的Nuavg随Ra增加而大幅提升；热源顶壁处的换热情况虽然得到改善，但其Nuavg仍远低于侧壁.即在高Ra下，热源侧壁处的换热是影响多孔介质内对流换热强度的重要因素，顶壁处的对流换热仍受到侧壁换热的影响，导致Nuavg仍处于较低水平.综上，对于本文讨论的多孔介质方腔对流换热问题，热壁面间的相互作用是需要考虑的关键因素.该研究能够为针对性地改善多孔介质内部换热提供相关的理论指导.

3.2 达西数变化对换热的影响

达西数是表征多孔介质渗透能力的无量纲准则数.在本文的工作中，讨论了Da=10-7～10-1时的多孔介质对流换热，其它模拟参数设置为Ra=106，ε=0.4，Pr=0.71，l=0.2L，h=0.2H，模拟温度场至稳态的部分结果如图6所示.

图6 不同达西数条件下多孔介质内的等温线(左)和流线(右)：(a) Da=10-4；(b) Da=10-3；(c) Da=10-2

根据模拟结果可知，Da较低时，流体的流动强度受限于多孔介质的渗透能力，导致对流换热强度处于较低水平；随着Da的增加，多孔介质的渗透能力逐渐增强，流体的流动强度逐渐增大.在Da=10-2时，由于多孔介质对流体部分的影响较小，温度与流线分布中存在显著的对流效应，整体流动换热情况与纯流体类似.故可以发现Da对多孔介质对流换热的影响与Ra类似，是决定多孔介质内部换热形式的主要因素.

Da的变化对垂直速度分布和Nuavg的影响如图7所示.根据速度曲线的变化可以发现，Da与Ra类似，它的增加能够强化流体的运动.但在Da>10-4后，随着Da的继续增加，方腔中的流体运动得到了整体增强，而Ra的增加则仅会使方腔的x/L=0与x/L=0.5附近出现显著的速度梯度.由图7(b)可知，热源各壁面的Nuavg随Da的增加而出现提升：在低Da的情况下，热源顶壁的换热情况优于侧壁，此时的换热性质主要为热传导.随着Da的增加，热源顶壁处的Nuavg逐渐降低，同时在Da>10-5时侧壁处的Nuavg开始高于顶壁处的Nuavg；当换热形式以热对流为主时，顶壁的换热能力开始与Da的提升成正比.结合上节的分析可以发现，在换热形式为热对流时，热源顶壁处的换热主要取决于Da和Ra.而在热传导时，热源侧壁的热交换成为了影响顶壁换热情况的主要因素.

图7 不同达西数条件下：(a)介质中心线y/H=0.5上流体垂直速度分布；(b)热源壁面上的Nuavg

3.3 孔隙率变化对换热的影响

孔隙率ε是由于孔隙的存在而引申出的密实固体所不曾具备的、多孔介质所特有的属性.为了研究不同的孔隙率对多孔介质内对流换热的影响，选取孔隙率为ε=0.1 ～ 0.9，同时设置对流换热参数为Da=10-3，Ra=106.其余模拟参数设置为Pr=0.71，l=0.2L，h=0.2H，模拟温度场至稳态的部分结果如图8所示.

图8 不同孔隙率条件下多孔介质内的等温线(左)和流线(右)：(a)ε=0.2；(b)ε=0.4；(c)ε=0.6

从图8中可以发现，孔隙率ε的变化对温度场的影响相对较小，但在孔隙率较高时，能够观察到温度更为有力的扩散现象，方腔顶部具有面积更大的高温区域.从图8的流线分布来看，多孔介质中的流体运动得到增强，但方腔内对称反向旋涡的形状与旋涡的中心位置几乎一致，整体的流动性质未能产生明显变化.由此可见，孔隙率的变化对于流场不具有显著影响，这与Ra和Da的情况不同.

垂直速度分布和Nuavg曲线随孔隙率的变化如图9所示.从图9中可以发现，随着孔隙率的增加，多孔介质中的流体运动得到显著强化，同时沿热源各个壁面的Nuavg出现明显提升，即孔隙率的增加能够有效改善多孔介质内对流换热效果.特别的是，由于在孔隙率较低时仍能够观察到明显的对流效应，这说明孔隙率的变化未能使多孔介质内的换热性质发生转变.故孔隙率并不是影响多孔介质内部的流动换热情况的主导因素.

图9 不同孔隙率条件下：(a)介质中心线y/H=0.5上流体垂直速度分布；(b)热源壁面上的Nuavg

4 结 论

本文采用LBM模拟底边中心处存在矩形热源的多孔介质方腔内的自然对流现象，讨论了不同的Ra、Da以及孔隙率ε对多孔介质对流换热特性的影响.研究结果表明：

1.当Ra>105时，多孔介质内对流换热水平随Ra增加而得到显著提升；当Ra≤105时，Ra对换热的影响可以忽略.

2.当Da>10-5时，Da的增加能有效改善多孔介质内的流动换热情况；当Da≤10-5时，Nuavg不随Da而产生明显变化.

3.孔隙率的增加能够强化流动换热，但不会改变多孔介质内换热性质.

4.当换热性质以热传导为主时，随Ra或Da的增加，沿热源侧壁处的Nuavg逐渐提升，并使热源顶壁处的热边界层变厚，导致沿顶壁处的Nuavg小幅降低.