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

文章编号：1007 6735(2023)05 0448 12DOI:10.13255/ki.jusst.20220531002不同冷热壁面对方腔内对流换热影响的格子 Boltzmann 模拟王婕1，张拴羊1，徐洪涛1，杨茉1,2（1.上海理工大学能源与动力工程学院，上海200093；2.上海建桥学院，上海201306）摘要：基于格子 Boltzmann 方法，选取多弛豫时间(MRT)模型对具有局部冷热壁面立体方腔中的三维自然对流换热进行模拟研究。分别采用三维 D3Q19 模型描述速度场，D3Q7 模型描述温度场，研究不同瑞利数 Ra（103Ra105）和局部冷热壁面位置变化对三维方腔内自然对流的影响。结果表明：冷热壁面的布置方式对流动换热有显著影响，且瑞利数 Ra 越大，其影响效果越明显。其中，冷热壁面均处于中间位置（工况 5）的自然对流换热能力最强，并且随着 Ra 逐渐增大，自然对流能力增强，流线逐渐复杂。在相同工况条件下，随着 Ra 增大，平均努塞尔数 Nuav也逐渐增加，换热能力增强。在相同 Ra 条件下，工况 5 的 Nuav最大，表明其自然对流换热能力最强。关键词：Boltzmann 模拟；多松弛时间模型；自然对流；三维方腔中图分类号：TK124文献标志码：ALattice Boltzmann simulation of convection heat transfer in a cubiccavity with different local cold and hot wall conditionsWANG Jie1，ZHANG Shuanyang1，XU Hongtao1，YANG Mo1,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 wasselectedtosimulatethethree-dimensionalnaturalconvectionheattransferinacubiccavitywithlocalhotandcoldwalls.TheD3Q19modelandtheD3Q7modelwereusedtodescribethevelocityfieldandtemperaturefield,respectively.TheeffectsofdifferentRayleighnumbersRa(103Ra105)andlocalhotandcoldwallpositiononthenaturalconvectioninthethree-dimensionalcubiccavitywerestudied.Theresultsrevealthatthearrangementofthehotandcoldwallshasasignificanteffectontheflowheattransfer,and the larger the Ra,the more obvious the effect.The strongest natural convection heattransfercapacityisobtainedbyCase5withthehotandcoldwallsinthemiddleposition.AstheRa上 海 理 工 大 学 学 报第45卷第5期J.UniversityofShanghaiforScienceandTechnologyVol.45No.52023收稿日期：20220531基金项目：上海市国际科技合作基金资助项目（18160743600）；上海市自然科学基金资助项目（20ZR1438700）第一作者：王婕（1987），女，博士研究生研究方向：传热传质数值模拟E-mail：wangjie_通信作者：徐洪涛（1976），男，教授研究方向：传热传质强化机理与多相流E-mail：increases gradually,the natural convection ability increases and the streamline becomes morecomplicated.Underthesameworkingcase,withtheincreaseofRa,theaverageNusseltnumberNuavalsoincreasesgradually,andtheheattransfercapacityisstronger.TheNuavvalueofCase5isthelargestatthesameRa,indicatingthatitsnaturalconvectionheattransfercapabilityisthestrongest.Keywords:Boltzmann simulation;multiple-relaxation-time model;natural convection;3D cubiccavity自然对流作为一种重要的物理现象，具有较强的应用背景，涉及房间通风、废料冷却及冶金工程等众多领域1。对于自然对流的研究主要采用实验和数值模拟两种方法。实验方法能够观察到宏观的驱动现象，但也存在实验时间长和费用高的局限性。随着计算机的快速发展，数值模拟方法在解决流体力学问题中发挥了重要作用，并已发展成为常用方法2。近年来兴起的格子 Boltzmann方法（latticeBoltzmannmethod，LBM）因其具有计算效率高、易并行和边界处理简单等优点3-4被广泛应用于传热传质的研究中。因此，许多关于自然对流的研究均采用 LBM。Xu 等5研究了在具有 Soret 和 Dufour 效应的方形二维封闭空间中，加热圆柱体周围的双扩散自然对流现象。研究表明，在高瑞利数 Ra=105条件下，随着 Soret 数Sr 和 Dufour 数 Df 的同时增加，传热传质能力减弱，流动的稳定性增强。Dadvand 等6探讨了二维方腔内不同热源和热沉位置、数量以及瑞利数Ra 对自然对流的影响，结果表明，由于旋涡的形成，热源和热沉的布置会强烈影响方腔中的流场和温度场，并且将热源和热沉交替放置在侧壁面，可以实现最大的传热量和最低的熵产。朱建奇等7通过对不同瑞利数 Ra 和倾斜角 条件下二维封闭方腔自然对流换热进行研究，发现随着Ra 的增大会提高自然对流的换热效果，然而倾斜角 的增加会导致自然对流换热交替式增减。Farkach 等8采用多松弛时间模型研究了具有离散隔热表面的加热圆柱与冷却圆柱之间的自然对流换热，研究发现，随着 Ra 的增加，传热速率也随之升高。另外，随着半径比的增大，平均努塞尔数 Nuav也相应增大，换热速率增强。上述文献都是基于二维模型进行 LBM 模拟研究，二维结果在一定程度上能够近似反映真实流体的流动，但基于三维空间的计算流体能更准确地反映实际问题并验证实验结果。因此，近年来学者们已经采用 LBM 研究三维条件下流体的流动与换热。Sheikholeslami 等9利用 LBM 模拟了磁场对三维方腔中纳米流体流动和传热的影响，结果表明，Nuav随着 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.5H 且温度为 Th的热壁面，右壁面设置与左壁面相同尺寸且温度为 Tc的冷壁面，其余为绝热壁面。基于此，进一步提出了5 种具有不同冷热壁面边界条件的工况进行模拟研究，如图 1 所示。1.2 数学模型假设方腔内流体为不可压缩的牛顿流体，基于 Boussinesq 假设，其无量纲方程描述为第5期王婕，等：不同冷热壁面对方腔内对流换热影响的格子 Boltzmann 模拟449UX+VY+WZ=0U+UUX+VUY+WUZ=PX+MaPr3Ra(2UX2+2UY2+2UZ2)V+UVX+VVY+WVZ=PY+MaPr3Ra(2VX2+2VY2+2VZ2)W+UWX+VWY+WWZ=PZ+MaPr3Ra(2WX2+2WY2+2WZ2)+Ma23+UX+VY+WZ=Ma13RaPr(2X2+2Y2+2Z2)（1）式（1）中引入的无量纲参数分别为X=xH，Y=yH，Z=zH，uc=g(ThTc)H，Ma=uccs，U=u3cs，V=v3cs,W=w3cs，Pr=a，P=p3c2s，=t3csH，Ra=g(ThTc)H3Pr2，=T TcThTc（2）式中：X，Y，Z 为笛卡尔坐标系中横轴、纵轴、竖轴的无量纲数；U，V，W 为 X，Y，Z 方向依次对应的速度；为热膨胀系数；cs为声速；Ma 为马赫数；P 为压力；为时间；g 为重力加速度；为无量纲温度；T 为温度；下标 h 表示高温侧，c 表示低温侧。高温左壁面的局部努塞尔数 Nu 定义为Nu=hk/H?Y=0,X0,1,Z0,1=Y?Y=0,X0,1,Z0,1（3）式中：h 为对流换热系数；k 为导热系数。高温热壁面的平均努塞尔数 Nuav定义为Nuav=wZ2Z1w10NudXdZ,Z1,Z2 0,1（4）式中，Z1，Z2分别为上、下壁面的取值。壁面采用无滑移边界条件，5 种不同工况下的速度和温度边界条件设置如表 1 所示。工况 1工况 2工况 3工况 4工况 5ZXYZXYZXYZXYZXY图 1 立方腔体中具有不同冷热壁面位置的物理模型Fig.1 Physical models with different hot and cold wall positions in a cubic cavity表 1 5 种不同工况的边界条件Tab.1 Boundary conditions for the five different cases工况边界条件15X=0,H;U=V=W=0Y=0,H;U=V=W=0Z=0,H;U=V=W=0Z=0,/Z=0;Z=H,/Z=0X=0,/X=0;X=H,/X=01Y=0,0.5 Z 1.0,=1;Y=H,0 Z 0.5,=0Y=0,0 Z 0.5,/Y=0;Y=H,0.5 Z 1.0,/Y=02Y=0,0.5 Z 1.0,=1;Y=H,0 Z 0.5,/Y=0Y=0,0 Z 0.5,/Y=0;Y=H,0.5 Z 1.0,=03Y=0,0.5 Z 1.0,/Y=0;Y=H,0 Z 0.5,/Y=0Y=0,0 Z 0.5,=1;Y=H,0.5 Z 1.0,=04Y=0,0.5 Z 1.0,/Y=0;Y=H,0 Z 0.5,=0Y=0,0 Z 0.5,=1;Y=H,0.5 Z 1.0,/Y=05Y=0,0 Z 0.25 0.75 Z 0.1,/Y=0;Y=H,0.25 Z 0.75,=0Y=H,0 Z 0.25 0.75 Z 0.1,/Y=0;Y=0,0.25 Z 0.75,=1450上海理工大学学报2023年第45卷 2 MRT-LBM 模型研究采用 D3Q19-MRT 和 D3Q7-MRT 双分布函数模型，采用 f 和 g 这 2 个分布函数分别来描述流场和温度场。速度场由 D3Q19-MRT 模型求解，其演化方程为12fi(r+eit,t+t)fi(r,t)=i+Fi,i=1,2,19（5）式中：ei为离散速度；Fi为源项；i为碰撞项。ei=c011000011111111000000011001111000011110000011000011111111（6）Fi=tG(eiV)pfeq,i(r,t)（7）G=(T T0)g（8）feq,i(r,t)=i1+eiVc2s+(eiV)22c2sVV2c2s,i=1,2,19（9）i=1/3,i=11/18,i=2,3,71/36,i=8,9,19（10）i=M1Smi(r,t)meq,i(r,t),i=1,2,19（11）feq,i(r,t)meq,i(r,t)式中：G 为有效重力项；为平衡态分布函数；V 为速度；T0为参考温度；i为密度权重系数；mi(r,t)和分别为分布函数和平衡态分布函数；M 和 S 分别为变换矩阵和碰撞矩阵。宏观流动密度、速度和矩阵为=19i=1fi(r,t),V=19i=1fi(r,t)m=Mf1,f2,f19T（12）式中，m 为流场的宏观矩阵。温度场由 D3Q7-MRT 模型求解，温度场的分布 gi可以表示为gi(r+uit,t+t)gi(r,t)=N1Qni(r,t)neq,i(r,t),i=1,2,7（13）neq,i(r,t)ui式中：为平衡态分布函数；N 和 Q 为能量分布的变换矩阵和碰撞矩阵；为离散速度。ui=c011000000011000000011（14）能量权重系数T,i=1/4,i=11/8,i=2,3,7（15）对应于能量分布的宏观参数为T=7i=1gi(r,t)n=Ng1,g2,g7T（16）式中，n 为温度场的宏观矩阵。使用该 D3Q7-MRT 模型求解温度场，并根据参考文献 13 中的设置求解热边界条件。Li 等12通过研究表明，当网格数为 505050时，LBM 在模拟三维方腔对流问题的研究中具有较高的准确性。因此，本文同样选取 505050 的网格数来进行模拟研究。3 程序验证为了验证本文采用的双 MRT-LBM 模型的准确性，分别与文献 14 和文献 15 的立体方腔内三维自然对流换热进行了对比验证。表 2 为 Ra=103,104，105条件下，高温左壁面的 Nuav。结果表明，本文方法的计算结果与前人的研究结果吻合较好，证明本程序具有一定的准确性。表 2 立体方腔内热壁面平均努塞尔数Nuav的比较Tab.2 Comparison of the average Nusselt number Nuav on thehot wall in a cubic cavityRa文献14文献15本文结果1031.0701.0751.0631042.0542.0852.0581054.3274.3784.370 4 结果与分析 4.1 工况 1 和工况 2 的特性分析图 2 分别为工况 1 和工况 2 在不同 Ra 条件下，方腔中等温线的三维分布图。工况 1 中的冷热壁面呈左上和右下的对角分布，从图 2 的工况1(a)可以看出，三维等温面的分布与中截面 X=0.5的温度分布较为均匀。这是因为在 Ra=103时，流体换热的方式以导热为主，在有限空间内自然对流受到了抑制。当 Ra=104时，传热能力进一步得到增强，高温面（低温面）相对侧的温度逐渐升高（降低）。当 Ra=105时，高温等温面沿 Y 轴正方向第5期王婕，等：不同冷热壁面对方腔内对流换热影响的格子 Boltzmann 模拟451水平方向移动，其原因在于自然对流在 Ra=105时变强，增强的热浮升力使得流体向上流动并和上壁面接触时沿 Y 轴正方向发展，同时底部的流体沿 Y 轴负方向发展。从图 2 的工况 1(b)可以看出，顶部和底部壁面之间的温差随着 Ra 的增加而变大，此时方腔中部的等温面几乎和 Y 轴保持平行。同时，随着 Ra 的增加，热源面下端的热边界层与冷源面上端的冷边界层厚度均变小，即温度梯度增大，可见传热显著增强。工况 2 中的冷热壁面均位于方腔两侧壁面的上半部。从图 2 的工况 2(c)可以看出，当 Ra=103时，靠近冷热壁面附近的等温面几乎与 X-Z 平面平行，且在 X-Z 平面内沿 X 轴方向的温度差异几乎可以忽略。当 Ra=104时，方腔上部的等温面出现沿 Y 轴正向移动的趋势，而方腔底部的温度梯度变小，传热能力较弱，这与冷热壁面的位置密0.950.900.850.800.750.700.650.600.550.500.450.400.350.300.250.200.150.050.100.950.900.850.800.750.700.650.600.550.500.450.400.350.300.250.200.150.050.100.950.900.850.800.750.700.650.600.550.500.450.400.350.300.250.200.150.050.100.950.900.850.800.750.700.650.600.550.500.450.400.350.300.250.200.150.050.10ZXYZXYZXYZXYZXYZXYZXYZXYZXYZXYZXYZXY(a)方腔内部温度分布云图(工况 1)(b)X=0.5 时温度分布云图(工况 1)(c)方腔内部温度分布云图(工况 2)(d)X=0.5 时温度分布云图(工况 2)Ra=103Ra=104Ra=105Ra=103Ra=104Ra=105Ra=103Ra=104Ra=105Ra=103Ra=104Ra=105图 2 工况 1 和工况 2 的温度分布云图Fig.2 Isotherm results of cases 1 and 2452上海理工大学学报2023年第45卷切相关。当 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 条件ZXYZXYZXYZXYZXYZXYZXYZXYZXY(a)3D 流线分布(工况 1)(b)X=0.5 时流线分布(工况 1)(c)Y=0.5 时流线分布(工况 1)Ra=103Ra=104Ra=105Ra=103Ra=104Ra=105Ra=103Ra=104Ra=105图 3 工况 1 和工况 2 的流线分布Fig.3 Streamline results of case 1 and 2第5期王婕，等：不同冷热壁面对方腔内对流换热影响的格子 Boltzmann 模拟453下，方腔中等温线的三维分布图。工况 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时，流体围绕方腔的ZXYZXYZXYZXYZXYZXYZXYZXYZXY(d)3D 流线分布(工况 2)(e)X=0.5 时流线分布(工况 2)(f)Y=0.5 时流线分布(工况 2)Ra=103Ra=104Ra=105Ra=103Ra=104Ra=105Ra=103Ra=104Ra=105图 3 工况 1 和工况 2 的流线分布（续）Fig.3 Streamline results of case 1 and 2454上海理工大学学报2023年第45卷中心流动，形成一个较大的涡流；当 Ra=105时，随着自然对流的增强，方腔的内部出现了 2 个涡流。从图 5(c)的工况 3 可知，在 Ra=103时，流体流向平面 Y=0.5 的 4 个角，而 Ra=104时，平面内形成了 4 个涡流，这表明三维自然对流在不断地增强。从图 5(d)和(e)的工况 4 可以看出，随着Ra 的增大，腔内的三维效应比较明显，尤其是在Ra=105时，腔内的流线形状发生了明显的变化。0.950.900.850.800.750.700.650.600.550.500.450.400.350.300.250.200.150.050.100.950.900.850.800.750.700.650.600.550.500.450.400.350.300.250.200.150.050.100.950.900.850.800.750.700.650.600.550.500.450.400.350.300.250.200.150.050.100.950.900.850.800.750.700.650.600.550.500.450.400.350.300.250.200.150.050.10ZXYZXYZXYZXYZXYZXYZXYZXYZXYZXYZXYZXY(b)X=0.5 时温度分布云图(工况 3)(d)X=0.5 时温度分布云图(工况 4)(a)方腔内部温度分布云图(工况 3)(c)方腔内部温度分布云图(工况 4)Ra=103Ra=104Ra=105Ra=103Ra=104Ra=105Ra=103Ra=104Ra=105Ra=103Ra=104Ra=105图 4 工况 3 和工况 4 的温度分布云图Fig.4 Isotherm results of cases 3 and 4第5期王婕，等：不同冷热壁面对方腔内对流换热影响的格子 Boltzmann 模拟455在中截面 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 中的冷热壁面分别位于ZXYZXYZXYZXYZXYZXYZXYZXYZXYZXYZXYZXY(a)方腔内部温度分布云图(工况 3)(b)X=0.5 时流线分布(工况 3)(c)Y=0.5 时流线分布(工况 3)(d)方腔内部温度分布云图(工况 4)Ra=103Ra=104Ra=105Ra=103Ra=104Ra=105Ra=103Ra=104Ra=105Ra=103Ra=104Ra=105图 5 工况 3 和工况 4 的流线分布Fig.5 Streamline results of case 3 and 4456上海理工大学学报2023年第45卷方腔两侧壁面的中间部分。由图 6(a)可以看出，当 Ra=103时，由于冷热壁面对称布置且处于中部位置，等温面从左到右的分布几乎是垂直的。随着 Ra 的增加，两侧壁面附近的等温面变得密集，尤其是 Ra=105时，方腔中心位置的等温面分布与X-Y 面几乎平行。这表明自然对流换热主要发生在方腔的中间区域。由图 6(b)可以看出，中截面X=0.5 内的温度分布具有二维平面的特征。随着ZXYZXYZXYZXYZXYZXY(e)X=0.5 时流线分布(工况 4)(f)Y=0.5 时流线分布(工况 4)Ra=103Ra=104Ra=105Ra=103Ra=104Ra=105图 5 工况 3 和工况 4 的流线分布（续）Fig.5 Streamline results of case 3 and 40.950.900.850.800.750.700.650.600.550.500.450.400.350.300.250.200.150.050.100.950.900.850.800.750.700.650.600.550.500.450.400.350.300.250.200.150.050.10ZXYZXYZXYZXYZXYZXY(a)方腔内部温度分布云图(b)X=0.5 时温度分布云图 Ra=103Ra=104Ra=105Ra=103Ra=104Ra=105图 6 工况 5 的温度分布云图Fig.6 Isotherm results of case 5第5期王婕，等：不同冷热壁面对方腔内对流换热影响的格子 Boltzmann 模拟457Ra 的增大，方腔顶部和底部的温差越来越明显，热边界层与冷边界层整体明显变薄，温度梯度增大，可见较大的 Ra 能有效地增强换热。图 7 为工况 5 在不同 Ra 条件下，方腔中流线的分布图。从图 7(a)和(b)看出，当 Ra 为 103，104时，两者的三维流场图的差别很小；而当Ra 增加到 105时，方腔中可以明显看出 2 个涡流的存在。考虑到工况 5 方腔中的冷热壁面的对称布置，由图 7(c)可以看出，当 Ra103时，在中截面 Y0.5 处出现了 4 个对称涡流，这同样是因为传热方式主要以导热为主引起的；当 Ra 增加到104时，在中截面 Y0.5 内出现了一个流线交汇处，这意味着流体流向 X=0.5 处；当 Ra=105时，中截面 Y=0.5 内仅出现了一个流线交汇处，其原因是自然对流换热的强化。中截面 Y=0.5 处的流线也证明了随着 Ra 的增大，垂直于 Y-Z 平面的流动越来越弱。4.4 不同冷热壁面位置对平均努塞尔数的影响表 3 为不同工况条件下，高温壁面的平均努塞尔数。整体而言，平均努塞尔数 Nuav随着 Ra 的增大而增大，这是因为 Ra 表征着自然对流的强弱，Ra 越大，换热越强。当 Ra=103时，换热方式主要是热传导，工况 14 方腔工况的 Nuav差异很小；对于工况 5 方腔，其 Nuav为 1.712，相对于其他 4 种工况的 Nuav更大，这是因为方腔的上下表面和冷热壁面之间的距离较大，有利于流体的流动 换 热。当 Ra=104时，不 同 工 况 下 热 壁 面 的Nuav差异变得很明显，因为冷热壁面位置是影响自然对热的重要因素。工况 1 方腔的 Nuav最小，其上壁面和下壁面阻碍了流体的流动；工况 5 的Nuav最大。当 Ra 增大到 105时，不同工况的高温壁面 Nuav的差异更为显著。工况 5 的 Nuav仍为最大，表明其自然对流换热能力最强，这是因为冷热壁面的布置方式减小了上下壁面对流体流动的ZXYZXYZXYZXYZXYZXYZXYZXYZXY(a)方腔内部温度分布云图(b)X=0.5 时流线分布(c)Y=0.5 时流线分布Ra=103Ra=104Ra=105Ra=103Ra=104Ra=105Ra=103Ra=104Ra=105图 7 工况 5 的流线分布Fig.7 Streamline results of case 5458上海理工大学学报2023年第45卷抑制作用，使得浮升力的作用效果更为明显。5 结论采用 D3Q19-MRT 和 D3Q7-MRT 双分布函数模型，对具有局部冷热壁面不同布置方式的三维方腔自然对流换热进行了数值模拟。通过分析三维方腔内部温度和流线的分布规律，研究了不同冷热壁面位置和 Ra（103Ra105）对三维方腔内自然对流换热的影响，并比较了不同工况下的平均努塞尔数 Nuav，得到如下结论：a.相同 Ra 条件下，冷热壁面的相对布置方式对流动换热有重要影响。工况 1 的对流换热能力最差，工况 5 的自然对流换热能力最强。在相同工况条件下，随着 Ra 从 103逐渐增大 105，自然对流能力增强，流线变得更为复杂。b.在相同工况条件下，随着 Ra 增大，Nuav也逐渐增加，换热能力增强。在相同 Ra 条件下，工况 5 的 Nuav均最大，表明其自然对流换热能力最强，说明冷热壁面的布置方式减小了上下壁面对流体流动的抑制作用，使浮升力的作用效果更为明显。采用格子 Boltzmann 方法仅针对具有局部冷热壁面不同布置方式的三维方腔自然对流进行了数值模拟，后续作者将基于该模拟结果进一步研究更为复杂条件下的自然对流换热问题，以期为解决不同条件下的自然对流换热问题提供一定参考。参考文献：杨世铭,陶文铨.传热学 M.4 版.北京:高等教育出版社,2006.1SHEN R Q,JIAO Z R,PARKER T,et al.Recentapplication of computational fluid dynamics(CFD)inprocesssafetyandlossprevention:areviewJ.JournalofLoss Prevention in the Process Industries,2020,67:104252.2张拴羊,徐洪涛,梁天生,等.Soret 和 Dufour 效应对方腔3内双扩散自然对流影响的格子 Boltzmann 模拟 J.空气动力学学报,2019,37(2):207215.余端民,赵明.圆内开缝圆环形空间自然对流的格子Boltzmann 模拟 J.上海理工大学学报,2019,41(2):108115.4XU H T,LUO Z Q,LOU Q,et al.Lattice Boltzmannsimulationsofthedouble-diffusivenaturalconvectionandoscillation characteristics in an enclosure with Soret andDufour effectsJ.International Journal of ThermalSciences,2019,136:159171.5DADVAND A,SARAEI S H,GHOREISHI S,et al.Lattice Boltzmann simulation of natural convection in asquareenclosurewithdiscreteheatingJ.MathematicsandComputersinSimulation,2021,179:265278.6朱建奇,侯佳煜,杲东彦,等.基于格子 Boltzmann 方法的倾斜方腔自然对流模拟 J.南京师范大学学报(工程技术版),2018,18(4):1926.7FARKACH Y,DERFOUFI S,AHACHAD M,et al.NumericalinvestigationofnaturalconvectioninconcentriccylinderpartiallyheatedbasedonMRT-latticeBoltzmannmethodJ.InternationalCommunicationsinHeatandMassTransfer,2022,132:105856.8SHEIKHOLESLAMI M,ELLAHI R.Three dimensionalmesoscopicsimulationofmagneticfieldeffectonnaturalconvection of nanofluidJ.International Journal of HeatandMassTransfer,2015,89:799808.9KARATAS H,DERBENTLI T.Natural convection inrectangular cavities with one active vertical wallJ.International Journal of Heat and Mass Transfer,2017,105:305315.10ZHAOY,WANGL,CHAIZH,etal.Comparativestudyofnaturalconvectionmeltinginsideacubiccavityusinganimprovedtwo-relaxation-timelatticeBoltzmannmodelJ.International Journal of Heat and Mass Transfer,2019,143:118449.11LIZ,YANGM,ZHANGYW.LatticeBoltzmannmethodsimulation of 3-D natural convection with double MRTmodelJ.InternationalJournalofHeatandMassTransfer,2016