不同冷热壁面对方腔内对流换热影响的格子Boltzmann模拟

王婕 张拴羊 徐洪涛 杨茉

摘要：基于格子Boltzmann方法，选取多弛豫时间（MRT）模型对具有局部冷热壁面立体方腔中的三维自然对流换热进行模拟研究。分别采用三维D3Q19模型描述速度场，D3Q7模型描述温度场，研究不同瑞利数Ra（103≤Ra≤105）和局部冷热壁面位置变化对三维方腔内自然对流的影响。结果表明：冷热壁面的布置方式对流动换热有显著影响，且瑞利数Ra越大，其影响效果越明显。其中，冷热壁面均处于中间位置（工况5）的自然对流换热能力最强，并且随着Ra逐渐增大，自然对流能力增强，流线逐渐复杂。在相同工况条件下，随着Ra增大，平均努塞尔数Nu av也逐渐增加，换热能力增强。在相同Ra条件下，工况5的Nu av最大，表明其自然对流换热能力最强。

关键词：Boltzmann模拟；多松弛时间模型；自然对流；三维方腔

中图分类号：TK 124

文献标志码：A

Lattice Boltzmann simulation of convection heat transfer in acubic cavity with different local cold and hot wall conditions

WANG Jie1，ZHANG Shuanyang 1，XU Hongtao 1，YANG Mo 1，2

（1.School of Energy and Power Engineering，University of Shanghai for Science and Technology，Shanghai 200093，China;2.Shanghai Jian Qiao University，Shanghai 201306，China）

Abstract：Based on the lattice Boltzmann method，the multiple-relaxation-time（MRT）model was selected to simulate the three-dimensional natural convection heat transfer in acubic cavity with local hot and cold walls.The D3Q19 model and the D3Q7 model were used to describe the velocity field and temperature field，respectively.The effects of different Rayleigh numbers Ra（103≤Ra≤105）and local hot and cold wall position on the natural convection in the three-dimensional cubic cavity were studied.The results reveal that the arrangement of the hot and cold walls has asignificant effect on the flow heat transfer，and the larger the Ra，the more obvious the effect.The strongest natural convection heat transfer capacity is obtained by Case 5 with the hot and cold walls in the middle position.As the Ra increases gradually，the natural convection ability increases and the streamline becomes more complicated.Under the same working case，with the increase of Ra，the average Nusselt number Nu av also increases gradually，and the heat transfer capacity is stronger.The Nu av value of Case 5 is the largest at the same Ra，indicating that its natural convection heat transfer capability is the strongest.

Keywords：Boltzmann simulation;multiple-relaxation-time model;natural convection;3D cubic cavity

自然對流作为一种重要的物理现象，具有较强的应用背景，涉及房间通风、废料冷却及冶金工程等众多领域[1]。对于自然对流的研究主要采用实验和数值模拟两种方法。实验方法能够观察到宏观的驱动现象，但也存在实验时间长和费用高的局限性。随着计算机的快速发展，数值模拟方法在解决流体力学问题中发挥了重要作用，并已发展成为常用方法[2]。近年来兴起的格子Boltzmann方法（lattice Boltzmann method，LBM）因其具有计算效率高、易并行和边界处理简单等优点[3-4]被广泛应用于传热传质的研究中。因此，许多关于自然对流的研究均采用LBM。Xu等[5]研究了在具有Soret和Dufour效应的方形二维封闭空间中，加热圆柱体周围的双扩散自然对流现象。研究表明，在高瑞利数Ra=105条件下，随着Soret数Sr和Dufour数Df的同时增加，传热传质能力减弱，流动的稳定性增强。Dadvand等[6]探讨了二维方腔内不同热源和热沉位置、数量以及瑞利数Ra对自然对流的影响，结果表明，由于旋涡的形成，热源和热沉的布置会强烈影响方腔中的流场和温度场，并且将热源和热沉交替放置在侧壁面，可以实现最大的传热量和最低的熵产。朱建奇等[7]通过对不同瑞利数Ra和倾斜角θ条件下二维封闭方腔自然对流换热进行研究，发现随着Ra的增大会提高自然对流的换热效果，然而倾斜角θ的增加会导致自然对流换热交替式增减。Farkach等[8]采用多松弛时间模型研究了具有离散隔热表面的加热圆柱与冷却圆柱之间的自然对流换热，研究发现，随着Ra的增加，传热速率也随之升高。另外，随着半径比的增大，平均努塞尔数Nu av也相应增大，换热速率增强。

上述文献都是基于二维模型进行LBM模拟研究，二维结果在一定程度上能够近似反映真实流体的流动，但基于三维空间的计算流体能更准确地反映实际问题并验证实验结果。因此，近年来学者们已经采用LBM研究三维条件下流体的流动与换热。Sheikholeslami等[9]利用LBM模拟了磁场对三维方腔中纳米流体流动和传热的影响，结果表明，Nu av随着Ra的增加和哈密顿数Ha的减少而增加，并且随着Ha的增加，传热效果得到进一步增强。Karatas等[10]通过对三维矩形方腔内的自然对流进行研究，发现等温线的密度会随着高宽比的增大而增加，且高温壁面附近的等温线也相对密集。Zhao等[11]采用LBM方法揭示了加热方向在不同Ra和普朗特数Pr下，对三维方腔内自然对流熔化的影响。数值结果表明，熔化效率随着Ra的增加而提高，两者的差异主要取决于Pr和加热方向。

综上所述，可以发现，利用双MRT-LBM模型对具有局部冷热壁面立方体方腔中三维自然对流的研究很少。因此，本文采用D3Q19-MRT和D3Q7-MRT双分布函数模型，模拟了具有局部冷热壁面的方腔内自然对流，并且通过分析等温面和流线的分布规律，研究了不同Ra和局部冷热壁面位置变化对三维方腔内自然对流的影响。

1计算模型

1.1物理模型

以边长为H的三维正立方体方腔作为研究的物理模型。方腔内介质为空气，方腔左壁面设置长度为H，高度为0.5 H且温度为T h的热壁面，右壁面设置与左壁面相同尺寸且温度为T c的冷壁面，其余为绝热壁面。基于此，进一步提出了5种具有不同冷热壁面边界条件的工况进行模拟研究，如图1所示。

1.2数学模型

假设方腔内流体为不可压缩的牛顿流体，基于Boussinesq假设，其无量纲方程描述为

式（1）中引入的无量纲参数分别为

式中：X，Y，Z为笛卡尔坐标系中横轴、纵轴、竖轴的无量纲数；U，V，W为X，Y，Z方向依次对应的速度；β为热膨胀系数；c s为声速；Ma为马赫数；P为压力；τ为时间；g为重力加速度；Θ为无量纲温度；T为温度；下标h表示高温侧，c表示低温侧。

高温左壁面的局部努塞尔数Nu定义为

式中：h为对流换热系数；k为导热系数。

高温热壁面的平均努塞尔数Nu av定义为

式中，Z 1，Z 2分別为上、下壁面的取值。

壁面采用无滑移边界条件，5种不同工况下的速度和温度边界条件设置如表1所示。

2 MRT-LBM模型

研究采用D3Q19-MRT和D3Q7-MRT双分布函数模型，采用f和g这2个分布函数分别来描述流场和温度场。速度场由D3Q19-MRT模型求解，其演化方程为[12]

式中：e i为离散速度；F i为源项；Ω i为碰撞项。

式中：G为有效重力项； feq;i（r; t）为平衡态分布函 数； V 为速度； T0 为参考温度； ωi 为密度权重系 数； mi（r，t） 和meq;i（r; t）分别为分布函数和平衡态分 布函数； M 和 S 分别为变换矩阵和碰撞矩阵。

宏观流动密度、速度和矩阵为

式中，m为流场的宏观矩阵。

温度场由D3Q7-MRT模型求解，温度场的分布g i可以表示为

式中： neq;i（r; t）为平衡态分布函数； N 和 Q 为能量 分布的变换矩阵和碰撞矩阵； ui 为离散速度。

能量权重系数

对应于能量分布的宏观参数为

式中，n为温度场的宏观矩阵。

使用该D3Q7-MRT模型求解温度场，并根据参考文献[13]中的设置求解热边界条件。

Li等[12]通过研究表明，当网格数为50×50×50时，LBM在模拟三维方腔对流问题的研究中具有较高的准确性。因此，本文同样选取50×50×50的网格数来进行模拟研究。

3程序验证

为了验证本文采用的双MRT-LBM模型的准确性，分别与文献[14]和文献[15]的立体方腔内三维自然对流换热进行了对比验证。表2为Ra=103，104，105条件下，高温左壁面的Nu av。结果表明，本文方法的计算结果与前人的研究结果吻合较好，证明本程序具有一定的准确性。

4结果与分析

4.1工况1和工况2的特性分析

图2分别为工况1和工况2在不同Ra条件下，方腔中等温线的三维分布图。工况1中的冷热壁面呈左上和右下的对角分布，从图2的工况1（a）可以看出，三维等温面的分布与中截面X=0.5的温度分布较为均匀。这是因为在Ra=103时，流体换热的方式以导热为主，在有限空间内自然对流受到了抑制。当Ra=104时，传热能力进一步得到增强，高温面（低温面）相对侧的温度逐渐升高（降低）。当Ra=105时，高温等温面沿Y轴正方向水平方向移动，其原因在于自然对流在Ra=105时变强，增强的热浮升力使得流体向上流动并和上壁面接触时沿Y轴正方向发展，同时底部的流体沿Y轴负方向发展。从图2的工况1（b）可以看出，顶部和底部壁面之间的温差随着Ra的增加而变大，此时方腔中部的等温面几乎和Y轴保持平行。同时，随着Ra的增加，热源面下端的热边界层与冷源面上端的冷边界层厚度均变小，即温度梯度增大，可见传热显著增强。

工况2中的冷热壁面均位于方腔两侧壁面的上半部。从图2的工况2（c）可以看出，当Ra=103时，靠近冷热壁面附近的等温面几乎与X-Z平面平行，且在X-Z平面内沿X轴方向的温度差异几乎可以忽略。当Ra=104时，方腔上部的等温面出现沿Y轴正向移动的趋势，而方腔底部的温度梯度变小，传热能力较弱，这与冷热壁面的位置密切相关。当Ra=105时，冷热壁面附近的温度梯度变大，换热能力增强，此时自然对流开始占主导地位。从图2的工况2（d）可以看出，当Ra=105时，对流现象主要集中在方腔的上部，底部的温度分布则比较均匀，且温差较小。随着Ra的增加，热源面下端的热边界层明显变薄，且冷边界层的厚度整体变薄，可见较大的Ra对换热增强的效果较为显著。

图3分别为工况1和工况2在不同Ra条件下，方腔中流线的分布图。从图3（a）和（b）的工况1可以看出，随着Ra的增大，三维流线的整体分布由整齐规律变得复杂。并且当Ra=103时，在中截面X=0.5中心出现了一个涡流，且流线分布呈圆形；当Ra=104时，由于方腔冷热壁面的对角分布，中截面X=0.5处的涡流由圆形转变为椭圆形；当Ra=105时，流动变得更为复杂，在靠近高温壁面和低温壁面处各出现了一个涡流。从图3（c）的工况1可以看出，在不同Ra下，在中截面Y=0.5内均出现了一个交汇处，且几乎不存在差异。从图3（d）和（e）的工况2可以看出，工况2的三维流线的整体分布同样随着Ra增加变得较为复杂。当Ra=103，104时，在中截面X=0.5的中心均出现了一个涡流；不同的是，Ra=104时的涡流出现向右平移的趋势，且形状呈扁平状。当Ra增大到105时，涡流继续向右偏移，并且在中截面X=0.5处的左下角出现了一个较小的涡流。这种现象是由冷热壁面的相对位置改变引起的。此外，从图3（f）的工况2可以看出，当Ra=103时，在中截面Y=0.5上形成了4个涡流；Ra=104时，底部的2个涡流消失且上部的涡流减小；随着Ra增大到105，所有涡流全部消失，且在平面的上方出现了一个流线交汇处。

4.2工况3和工况4的特性分析

图4分别为工况3和工况4在不同Ra条件下，方腔中等温线的三维分布图。工况3中的冷热壁面呈现与工况1相反的对角分布。从图4（a）和（b）的工况3可以看出，随着Ra的增大，靠近冷热壁面附近的温度梯度逐渐增大，流动换热增强。当Ra=103时，大部分的等温面几乎与X-Z平面平行且均匀分布，因为方腔中的传热主要形式为热传导。随着Ra的逐渐增加，对流开始逐渐占据主导地位，中间区域的等温面几乎呈水平分布。从图4（b）的工况3可以看出，随着Ra增大到105，腔体中心区域的温差逐渐变小，而方腔顶部和底部的温差越来越大，这是对流换热增强的结果。较工况1和工况2不同的是，随着Ra的增加，热边界层与冷边界层整体明显变薄，温度梯度增大，可见较大的Ra对工况3换热增强的效果更加显著。

工况4中的冷热壁面均位于方腔的下半部分，与工况2相反。从图4（c）的工况4可以看出，随着Ra的增大，方腔内换热增强，方腔顶部与底部的温差越来越大。这是由于热壁面处于左壁面下侧，其上方存在相对较大的空间，有利于传热的发展；而冷壁面位于右侧壁面下侧，底部壁面对流动换热起到了一定的阻碍作用。从图4（d）的工况4可以看出，當Ra=103时，此时热传导处于主导地位，其大部分等温面几乎垂直于X-Y平面，随着Ra的增大，高温等温面逐渐向右侧发生了弯曲。其中，在工况4方腔中，上部的等温面非常稀疏，下部的等温面相对密集，这是因为Ra=105时，方腔下部自然对流换热的强化所致。随着Ra的增大，热边界层与冷源面上端的冷边界层厚度逐渐变小，温度梯度增大，同样可见传热显著增强。

图5分别为工况3和工况4在不同Ra条件下，方腔中流线的分布图。从图5（a）和（b）的工况3可以看出，当Ra=103，104时，流体围绕方腔的中心流动，形成一个较大的涡流；当Ra=105时，随着自然对流的增强，方腔的内部出现了2个涡流。从图5（c）的工况3可知，在Ra=103时，流体流向平面Y=0.5的4个角，而Ra=104时，平面内形成了4个涡流，这表明三维自然对流在不断地增强。从图5（d）和（e）的工况4可以看出，随着Ra的增大，腔内的三维效应比较明显，尤其是在Ra=105时，腔内的流线形状发生了明显的变化。在中截面X=0.5上，不同Ra下均出现了一个涡流，在Ra=104时，涡流的形状呈现出扁平的趋势；随着Ra增大到105，该涡流表现出向方腔的左上角移动的趋势。从图5（f）的工况4可以看出，当Ra=103时，在中截面Y=0.5上出现了4个涡流；随着Ra增加到104，平面内的涡流减少为2个，当Ra=105时，所有的涡流全部消失，仅在下部出现了一个流线交汇处。这主要是由自然对流强度不断加强导致的，这与工况2方腔的情况类似。

4.3工况5的特性分析

图6为工况5在不同Ra条件下，方腔中等温线的三维分布图。工况5中的冷热壁面分别位于方腔两侧壁面的中间部分。由图6（a）可以看出，当Ra=103时，由于冷热壁面对称布置且处于中部位置，等温面从左到右的分布几乎是垂直的。随着Ra的增加，两侧壁面附近的等温面变得密集，尤其是Ra=105时，方腔中心位置的等温面分布与X-Y面几乎平行。这表明自然对流换热主要发生在方腔的中间区域。由图6（b）可以看出，中截面X=0.5内的温度分布具有二维平面的特征。随着Ra的增大，方腔顶部和底部的温差越来越明显，热边界层与冷边界层整体明显变薄，温度梯度增大，可见较大的Ra能有效地增强换热。

图7为工况5在不同Ra条件下，方腔中流线的分布图。从图7（a）和（b）看出，当Ra为103，104时，两者的三维流场图的差别很小；而当Ra增加到105时，方腔中可以明显看出2个涡流的存在。考虑到工况5方腔中的冷热壁面的对称布置，由图7（c）可以看出，当Ra＝103时，在中截面Y＝0.5处出现了4个对称涡流，这同样是因为传热方式主要以导热为主引起的；当Ra增加到104时，在中截面Y＝0.5内出现了一个流线交汇处，这意味着流体流向X=0.5处；当Ra=105时，中截面Y=0.5内仅出现了一个流线交汇处，其原因是自然对流换热的强化。中截面Y=0.5处的流线也证明了随着Ra的增大，垂直于Y-Z平面的流动越来越弱。

4.4不同冷热壁面位置对平均努塞尔数的影响

表3为不同工况条件下，高温壁面的平均努塞尔数。整体而言，平均努塞尔数Nu av随着Ra的增大而增大，这是因为Ra表征着自然对流的强弱，Ra越大，换热越强。当Ra=103时，换热方式主要是热传导，工况1?4方腔工况的Nu av差异很小；对于工况5方腔，其Nu av为1.712，相对于其他4种工况的Nu av更大，这是因为方腔的上下表面和冷热壁面之间的距离较大，有利于流体的流动换热。当Ra=104时，不同工况下热壁面的Nu av差异变得很明显，因为冷热壁面位置是影响自然对热的重要因素。工况1方腔的Nu av最小，其上壁面和下壁面阻碍了流体的流动；工况5的Nu av最大。当Ra增大到105时，不同工况的高温壁面Nu av的差异更为显著。工况5的Nu av仍为最大，表明其自然对流换热能力最强，这是因为冷热壁面的布置方式减小了上下壁面对流体流动的抑制作用，使得浮升力的作用效果更为明显。

5结 论

采用D3Q19-MRT和D3Q7-MRT双分布函数模型，对具有局部冷热壁面不同布置方式的三维方腔自然对流换热进行了数值模拟。通过分析三维方腔内部温度和流线的分布规律，研究了不同冷热壁面位置和Ra（103≤Ra≤105）对三维方腔内自然对流换热的影响，并比较了不同工况下的平均努塞尔数Nu av，得到如下结论：

a.相同Ra条件下，冷热壁面的相对布置方式对流动换热有重要影响。工况1的对流换热能力最差，工况5的自然对流换热能力最强。在相同工况条件下，随着Ra从103逐渐增大105，自然对流能力增强，流线变得更为复杂。

b.在相同工况条件下，随着Ra增大，Nu av也逐渐增加，换热能力增强。在相同Ra条件下，工况5的Nu av均最大，表明其自然对流换热能力最强，说明冷热壁面的布置方式减小了上下壁面对流体流动的抑制作用，使浮升力的作用效果更为明显。

采用格子Boltzmann方法仅针对具有局部冷热壁面不同布置方式的三维方腔自然对流进行了数值模拟，后续作者将基于该模拟结果进一步研究更为复杂条件下的自然对流换热问题，以期为解决不同条件下的自然对流换热问题提供一定参考。

参考文献：

[1] 杨世铭， 陶文铨. 传热学 [M]. 4 版. 北京 ： 高等教育出版 社， 2006.

[2] SHEN R Q， JIAO Z R， PARKER T， et al. Recent application of computational fluid dynamics （CFD） in process safety and loss prevention： a review[J]. Journal of Loss Prevention in the Process Industries， 2020， 67： 104252.

[3] 张拴羊， 徐洪涛， 梁天生， 等. Soret 和 Dufour效应对方腔内双扩散自然对流影响的格子 Boltzmann模拟 [J]. 空气 动力学学报， 2019， 37（2）： 207–215.

[4] 余端民 ， 赵明 . 圆内开缝圆环形空间自然对流的格子 Boltzmann 模 拟 [J]. 上 海 理 工 大 学 学报 ， 2019， 41（2）： 108–115.

[5] XU H T， LUO Z Q， LOU Q， et al. Lattice Boltzmann simulations of the double-diffusive natural convection and oscillation characteristics in an enclosure with Soret and Dufour effects[J]. International Journal of Thermal Sciences， 2019， 136： 159–171.

[6] DADVAND A， SARAEI S H， GHOREISHI S， et al. Lattice Boltzmann simulation of natural convection in a square enclosure with discrete heating[J]. Mathematics and Computers in Simulation， 2021， 179： 265–278.

[7] 朱建奇， 侯佳煜， 杲东彦， 等. 基于格子 Boltzmann方法的 倾斜方腔自然对流模拟 [J].南京师范大学学报（工程技 术版）， 2018， 18（4）： 19–26.

[8] FARKACH Y， DERFOUFI S， AHACHAD M， et al. Numerical investigation of natural convection in concentric cylinder partially heated based on MRT-lattice Boltzmann method[J]. International Communications in Heat and Mass Transfer， 2022， 132： 105856.

[9] SHEIKHOLESLAMI M， ELLAHI R. Three dimensional mesoscopic simulation of magnetic field effect on natural convection of nanofluid[J]. International Journal of Heat and Mass Transfer， 2015， 89： 799–808.

[10] KARATAS H， DERBENTLI T. Natural convection in rectangular cavities with one active vertical wall[J]. International Journal of Heat and Mass Transfer， 2017， 105： 305–315.

[11] ZHAO Y， WANG L， CHAI Z H， et al. Comparative study of natural convection melting inside a cubic cavity using an improved two-relaxation-time lattice Boltzmann model[J]. International Journal of Heat and Mass Transfer， 2019， 143： 118449.

[12] LI Z， YANG M， ZHANG Y W. Lattice Boltzmann method simulation of 3-D natural convection with double MRT model[J]. International Journal of Heat and Mass Transfer， 2016， 94： 222–238.

[13] LI L K， MEI R W， KLAUSNER J F. Boundary conditions for thermal lattice Boltzmann equation method[J]. Journal of Computational Physics， 2013， 237： 366–395.

[14] TRIC E， LABROSSE G， BETROUNI M. A first incursion into the 3D structure of natural convection of air in a differentially heated cubic cavity， from accurate numerical solutions[J]. International Journal of Heat and Mass Transfer， 2000， 43（21）： 4043–4056.

[15] PENG Y， SHU C， CHEW Y T. A 3D incompressible thermal lattice Boltzmann model and its application to simulate natural convection in a cubic cavity[J]. Journal of Computational Physics， 2004， 193（1）： 260–274.

（編辑：石 瑛）