王婕 張拴羊 徐洪濤 楊茉

摘要：基于格子Boltzmann方法，選取多弛豫時(shí)間（MRT）模型對具有局部冷熱壁面立體方腔中的三維自然對流換熱進(jìn)行模擬研究。分別采用三維D3Q19模型描述速度場，D3Q7模型描述溫度場，研究不同瑞利數(shù)Ra（103≤Ra≤105）和局部冷熱壁面位置變化對三維方腔內(nèi)自然對流的影響。結(jié)果表明：冷熱壁面的布置方式對流動(dòng)換熱有顯著影響，且瑞利數(shù)Ra越大，其影響效果越明顯。其中，冷熱壁面均處于中間位置（工況5）的自然對流換熱能力最強(qiáng)，并且隨著Ra逐漸增大，自然對流能力增強(qiáng)，流線逐漸復(fù)雜。在相同工況條件下，隨著Ra增大，平均努塞爾數(shù)Nu av也逐漸增加，換熱能力增強(qiáng)。在相同Ra條件下，工況5的Nu av最大，表明其自然對流換熱能力最強(qiáng)。

關(guān)鍵詞：Boltzmann模擬；多松弛時(shí)間模型；自然對流；三維方腔

中圖分類號(hào)：TK 124

文獻(xiàn)標(biāo)志碼：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

自然對流作為一種重要的物理現(xiàn)象，具有較強(qiáng)的應(yīng)用背景，涉及房間通風(fēng)、廢料冷卻及冶金工程等眾多領(lǐng)域[1]。對于自然對流的研究主要采用實(shí)驗(yàn)和數(shù)值模擬兩種方法。實(shí)驗(yàn)方法能夠觀察到宏觀的驅(qū)動(dòng)現(xiàn)象，但也存在實(shí)驗(yàn)時(shí)間長和費(fèi)用高的局限性。隨著計(jì)算機(jī)的快速發(fā)展，數(shù)值模擬方法在解決流體力學(xué)問題中發(fā)揮了重要作用，并已發(fā)展成為常用方法[2]。近年來興起的格子Boltzmann方法（lattice Boltzmann method，LBM）因其具有計(jì)算效率高、易并行和邊界處理簡單等優(yōu)點(diǎn)[3-4]被廣泛應(yīng)用于傳熱傳質(zhì)的研究中。因此，許多關(guān)于自然對流的研究均采用LBM。Xu等[5]研究了在具有Soret和Dufour效應(yīng)的方形二維封閉空間中，加熱圓柱體周圍的雙擴(kuò)散自然對流現(xiàn)象。研究表明，在高瑞利數(shù)Ra=105條件下，隨著Soret數(shù)Sr和Dufour數(shù)Df的同時(shí)增加，傳熱傳質(zhì)能力減弱，流動(dòng)的穩(wěn)定性增強(qiáng)。Dadvand等[6]探討了二維方腔內(nèi)不同熱源和熱沉位置、數(shù)量以及瑞利數(shù)Ra對自然對流的影響，結(jié)果表明，由于旋渦的形成，熱源和熱沉的布置會(huì)強(qiáng)烈影響方腔中的流場和溫度場，并且將熱源和熱沉交替放置在側(cè)壁面，可以實(shí)現(xiàn)最大的傳熱量和最低的熵產(chǎn)。朱建奇等[7]通過對不同瑞利數(shù)Ra和傾斜角θ條件下二維封閉方腔自然對流換熱進(jìn)行研究，發(fā)現(xiàn)隨著Ra的增大會(huì)提高自然對流的換熱效果，然而傾斜角θ的增加會(huì)導(dǎo)致自然對流換熱交替式增減。Farkach等[8]采用多松弛時(shí)間模型研究了具有離散隔熱表面的加熱圓柱與冷卻圓柱之間的自然對流換熱，研究發(fā)現(xiàn)，隨著Ra的增加，傳熱速率也隨之升高。另外，隨著半徑比的增大，平均努塞爾數(shù)Nu av也相應(yīng)增大，換熱速率增強(qiáng)。

上述文獻(xiàn)都是基于二維模型進(jìn)行LBM模擬研究，二維結(jié)果在一定程度上能夠近似反映真實(shí)流體的流動(dòng)，但基于三維空間的計(jì)算流體能更準(zhǔn)確地反映實(shí)際問題并驗(yàn)證實(shí)驗(yàn)結(jié)果。因此，近年來學(xué)者們已經(jīng)采用LBM研究三維條件下流體的流動(dòng)與換熱。Sheikholeslami等[9]利用LBM模擬了磁場對三維方腔中納米流體流動(dòng)和傳熱的影響，結(jié)果表明，Nu av隨著Ra的增加和哈密頓數(shù)Ha的減少而增加，并且隨著Ha的增加，傳熱效果得到進(jìn)一步增強(qiáng)。Karatas等[10]通過對三維矩形方腔內(nèi)的自然對流進(jìn)行研究，發(fā)現(xiàn)等溫線的密度會(huì)隨著高寬比的增大而增加，且高溫壁面附近的等溫線也相對密集。Zhao等[11]采用LBM方法揭示了加熱方向在不同Ra和普朗特?cái)?shù)Pr下，對三維方腔內(nèi)自然對流熔化的影響。數(shù)值結(jié)果表明，熔化效率隨著Ra的增加而提高，兩者的差異主要取決于Pr和加熱方向。

綜上所述，可以發(fā)現(xiàn)，利用雙MRT-LBM模型對具有局部冷熱壁面立方體方腔中三維自然對流的研究很少。因此，本文采用D3Q19-MRT和D3Q7-MRT雙分布函數(shù)模型，模擬了具有局部冷熱壁面的方腔內(nèi)自然對流，并且通過分析等溫面和流線的分布規(guī)律，研究了不同Ra和局部冷熱壁面位置變化對三維方腔內(nèi)自然對流的影響。

1計(jì)算模型

1.1物理模型

以邊長為H的三維正立方體方腔作為研究的物理模型。方腔內(nèi)介質(zhì)為空氣，方腔左壁面設(shè)置長度為H，高度為0.5 H且溫度為T h的熱壁面，右壁面設(shè)置與左壁面相同尺寸且溫度為T c的冷壁面，其余為絕熱壁面?；诖耍M(jìn)一步提出了5種具有不同冷熱壁面邊界條件的工況進(jìn)行模擬研究，如圖1所示。

1.2數(shù)學(xué)模型

假設(shè)方腔內(nèi)流體為不可壓縮的牛頓流體，基于Boussinesq假設(shè)，其無量綱方程描述為

式（1）中引入的無量綱參數(shù)分別為

式中：X，Y，Z為笛卡爾坐標(biāo)系中橫軸、縱軸、豎軸的無量綱數(shù)；U，V，W為X，Y，Z方向依次對應(yīng)的速度；β為熱膨脹系數(shù)；c s為聲速；Ma為馬赫數(shù)；P為壓力；τ為時(shí)間；g為重力加速度；Θ為無量綱溫度；T為溫度；下標(biāo)h表示高溫側(cè)，c表示低溫側(cè)。

高溫左壁面的局部努塞爾數(shù)Nu定義為

式中：h為對流換熱系數(shù)；k為導(dǎo)熱系數(shù)。

高溫?zé)岜诿娴钠骄麪枖?shù)Nu av定義為

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

壁面采用無滑移邊界條件，5種不同工況下的速度和溫度邊界條件設(shè)置如表1所示。

2 MRT-LBM模型

研究采用D3Q19-MRT和D3Q7-MRT雙分布函數(shù)模型，采用f和g這2個(gè)分布函數(shù)分別來描述流場和溫度場。速度場由D3Q19-MRT模型求解，其演化方程為[12]

式中：e i為離散速度；F i為源項(xiàng)；Ω i為碰撞項(xiàng)。

式中：G為有效重力項(xiàng)； feq;i（r; t）為平衡態(tài)分布函 數(shù)； V 為速度； T0 為參考溫度； ωi 為密度權(quán)重系 數(shù)； mi（r，t） 和meq;i（r; t）分別為分布函數(shù)和平衡態(tài)分 布函數(shù)； M 和 S 分別為變換矩陣和碰撞矩陣。

宏觀流動(dòng)密度、速度和矩陣為

式中，m為流場的宏觀矩陣。

溫度場由D3Q7-MRT模型求解，溫度場的分布g i可以表示為

式中： neq;i（r; t）為平衡態(tài)分布函數(shù)； N 和 Q 為能量 分布的變換矩陣和碰撞矩陣； ui 為離散速度。

能量權(quán)重系數(shù)

對應(yīng)于能量分布的宏觀參數(shù)為

式中，n為溫度場的宏觀矩陣。

使用該D3Q7-MRT模型求解溫度場，并根據(jù)參考文獻(xiàn)[13]中的設(shè)置求解熱邊界條件。

Li等[12]通過研究表明，當(dāng)網(wǎng)格數(shù)為50×50×50時(shí)，LBM在模擬三維方腔對流問題的研究中具有較高的準(zhǔn)確性。因此，本文同樣選取50×50×50的網(wǎng)格數(shù)來進(jìn)行模擬研究。

3程序驗(yàn)證

為了驗(yàn)證本文采用的雙MRT-LBM模型的準(zhǔn)確性，分別與文獻(xiàn)[14]和文獻(xiàn)[15]的立體方腔內(nèi)三維自然對流換熱進(jìn)行了對比驗(yàn)證。表2為Ra=103，104，105條件下，高溫左壁面的Nu av。結(jié)果表明，本文方法的計(jì)算結(jié)果與前人的研究結(jié)果吻合較好，證明本程序具有一定的準(zhǔn)確性。

4結(jié)果與分析

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

圖2分別為工況1和工況2在不同Ra條件下，方腔中等溫線的三維分布圖。工況1中的冷熱壁面呈左上和右下的對角分布，從圖2的工況1（a）可以看出，三維等溫面的分布與中截面X=0.5的溫度分布較為均勻。這是因?yàn)樵赗a=103時(shí)，流體換熱的方式以導(dǎo)熱為主，在有限空間內(nèi)自然對流受到了抑制。當(dāng)Ra=104時(shí)，傳熱能力進(jìn)一步得到增強(qiáng)，高溫面（低溫面）相對側(cè)的溫度逐漸升高（降低）。當(dāng)Ra=105時(shí)，高溫等溫面沿Y軸正方向水平方向移動(dòng)，其原因在于自然對流在Ra=105時(shí)變強(qiáng)，增強(qiáng)的熱浮升力使得流體向上流動(dòng)并和上壁面接觸時(shí)沿Y軸正方向發(fā)展，同時(shí)底部的流體沿Y軸負(fù)方向發(fā)展。從圖2的工況1（b）可以看出，頂部和底部壁面之間的溫差隨著Ra的增加而變大，此時(shí)方腔中部的等溫面幾乎和Y軸保持平行。同時(shí)，隨著Ra的增加，熱源面下端的熱邊界層與冷源面上端的冷邊界層厚度均變小，即溫度梯度增大，可見傳熱顯著增強(qiáng)。

工況2中的冷熱壁面均位于方腔兩側(cè)壁面的上半部。從圖2的工況2（c）可以看出，當(dāng)Ra=103時(shí)，靠近冷熱壁面附近的等溫面幾乎與X-Z平面平行，且在X-Z平面內(nèi)沿X軸方向的溫度差異幾乎可以忽略。當(dāng)Ra=104時(shí)，方腔上部的等溫面出現(xiàn)沿Y軸正向移動(dòng)的趨勢，而方腔底部的溫度梯度變小，傳熱能力較弱，這與冷熱壁面的位置密切相關(guān)。當(dāng)Ra=105時(shí)，冷熱壁面附近的溫度梯度變大，換熱能力增強(qiáng)，此時(shí)自然對流開始占主導(dǎo)地位。從圖2的工況2（d）可以看出，當(dāng)Ra=105時(shí)，對流現(xiàn)象主要集中在方腔的上部，底部的溫度分布則比較均勻，且溫差較小。隨著Ra的增加，熱源面下端的熱邊界層明顯變薄，且冷邊界層的厚度整體變薄，可見較大的Ra對換熱增強(qiáng)的效果較為顯著。

圖3分別為工況1和工況2在不同Ra條件下，方腔中流線的分布圖。從圖3（a）和（b）的工況1可以看出，隨著Ra的增大，三維流線的整體分布由整齊規(guī)律變得復(fù)雜。并且當(dāng)Ra=103時(shí)，在中截面X=0.5中心出現(xiàn)了一個(gè)渦流，且流線分布呈圓形；當(dāng)Ra=104時(shí)，由于方腔冷熱壁面的對角分布，中截面X=0.5處的渦流由圓形轉(zhuǎn)變?yōu)闄E圓形；當(dāng)Ra=105時(shí)，流動(dòng)變得更為復(fù)雜，在靠近高溫壁面和低溫壁面處各出現(xiàn)了一個(gè)渦流。從圖3（c）的工況1可以看出，在不同Ra下，在中截面Y=0.5內(nèi)均出現(xiàn)了一個(gè)交匯處，且?guī)缀醪淮嬖诓町悺膱D3（d）和（e）的工況2可以看出，工況2的三維流線的整體分布同樣隨著Ra增加變得較為復(fù)雜。當(dāng)Ra=103，104時(shí)，在中截面X=0.5的中心均出現(xiàn)了一個(gè)渦流；不同的是，Ra=104時(shí)的渦流出現(xiàn)向右平移的趨勢，且形狀呈扁平狀。當(dāng)Ra增大到105時(shí)，渦流繼續(xù)向右偏移，并且在中截面X=0.5處的左下角出現(xiàn)了一個(gè)較小的渦流。這種現(xiàn)象是由冷熱壁面的相對位置改變引起的。此外，從圖3（f）的工況2可以看出，當(dāng)Ra=103時(shí)，在中截面Y=0.5上形成了4個(gè)渦流；Ra=104時(shí)，底部的2個(gè)渦流消失且上部的渦流減?。浑S著Ra增大到105，所有渦流全部消失，且在平面的上方出現(xiàn)了一個(gè)流線交匯處。

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

圖4分別為工況3和工況4在不同Ra條件下，方腔中等溫線的三維分布圖。工況3中的冷熱壁面呈現(xiàn)與工況1相反的對角分布。從圖4（a）和（b）的工況3可以看出，隨著Ra的增大，靠近冷熱壁面附近的溫度梯度逐漸增大，流動(dòng)換熱增強(qiáng)。當(dāng)Ra=103時(shí)，大部分的等溫面幾乎與X-Z平面平行且均勻分布，因?yàn)榉角恢械膫鳠嶂饕问綖闊醾鲗?dǎo)。隨著Ra的逐漸增加，對流開始逐漸占據(jù)主導(dǎo)地位，中間區(qū)域的等溫面幾乎呈水平分布。從圖4（b）的工況3可以看出，隨著Ra增大到105，腔體中心區(qū)域的溫差逐漸變小，而方腔頂部和底部的溫差越來越大，這是對流換熱增強(qiáng)的結(jié)果。較工況1和工況2不同的是，隨著Ra的增加，熱邊界層與冷邊界層整體明顯變薄，溫度梯度增大，可見較大的Ra對工況3換熱增強(qiáng)的效果更加顯著。

工況4中的冷熱壁面均位于方腔的下半部分，與工況2相反。從圖4（c）的工況4可以看出，隨著Ra的增大，方腔內(nèi)換熱增強(qiáng)，方腔頂部與底部的溫差越來越大。這是由于熱壁面處于左壁面下側(cè)，其上方存在相對較大的空間，有利于傳熱的發(fā)展；而冷壁面位于右側(cè)壁面下側(cè)，底部壁面對流動(dòng)換熱起到了一定的阻礙作用。從圖4（d）的工況4可以看出，當(dāng)Ra=103時(shí)，此時(shí)熱傳導(dǎo)處于主導(dǎo)地位，其大部分等溫面幾乎垂直于X-Y平面，隨著Ra的增大，高溫等溫面逐漸向右側(cè)發(fā)生了彎曲。其中，在工況4方腔中，上部的等溫面非常稀疏，下部的等溫面相對密集，這是因?yàn)镽a=105時(shí)，方腔下部自然對流換熱的強(qiáng)化所致。隨著Ra的增大，熱邊界層與冷源面上端的冷邊界層厚度逐漸變小，溫度梯度增大，同樣可見傳熱顯著增強(qiáng)。

圖5分別為工況3和工況4在不同Ra條件下，方腔中流線的分布圖。從圖5（a）和（b）的工況3可以看出，當(dāng)Ra=103，104時(shí)，流體圍繞方腔的中心流動(dòng)，形成一個(gè)較大的渦流；當(dāng)Ra=105時(shí)，隨著自然對流的增強(qiáng)，方腔的內(nèi)部出現(xiàn)了2個(gè)渦流。從圖5（c）的工況3可知，在Ra=103時(shí)，流體流向平面Y=0.5的4個(gè)角，而Ra=104時(shí)，平面內(nèi)形成了4個(gè)渦流，這表明三維自然對流在不斷地增強(qiáng)。從圖5（d）和（e）的工況4可以看出，隨著Ra的增大，腔內(nèi)的三維效應(yīng)比較明顯，尤其是在Ra=105時(shí)，腔內(nèi)的流線形狀發(fā)生了明顯的變化。在中截面X=0.5上，不同Ra下均出現(xiàn)了一個(gè)渦流，在Ra=104時(shí)，渦流的形狀呈現(xiàn)出扁平的趨勢；隨著Ra增大到105，該渦流表現(xiàn)出向方腔的左上角移動(dòng)的趨勢。從圖5（f）的工況4可以看出，當(dāng)Ra=103時(shí)，在中截面Y=0.5上出現(xiàn)了4個(gè)渦流；隨著Ra增加到104，平面內(nèi)的渦流減少為2個(gè)，當(dāng)Ra=105時(shí)，所有的渦流全部消失，僅在下部出現(xiàn)了一個(gè)流線交匯處。這主要是由自然對流強(qiáng)度不斷加強(qiáng)導(dǎo)致的，這與工況2方腔的情況類似。

4.3工況5的特性分析

圖6為工況5在不同Ra條件下，方腔中等溫線的三維分布圖。工況5中的冷熱壁面分別位于方腔兩側(cè)壁面的中間部分。由圖6（a）可以看出，當(dāng)Ra=103時(shí)，由于冷熱壁面對稱布置且處于中部位置，等溫面從左到右的分布幾乎是垂直的。隨著Ra的增加，兩側(cè)壁面附近的等溫面變得密集，尤其是Ra=105時(shí)，方腔中心位置的等溫面分布與X-Y面幾乎平行。這表明自然對流換熱主要發(fā)生在方腔的中間區(qū)域。由圖6（b）可以看出，中截面X=0.5內(nèi)的溫度分布具有二維平面的特征。隨著Ra的增大，方腔頂部和底部的溫差越來越明顯，熱邊界層與冷邊界層整體明顯變薄，溫度梯度增大，可見較大的Ra能有效地增強(qiáng)換熱。

圖7為工況5在不同Ra條件下，方腔中流線的分布圖。從圖7（a）和（b）看出，當(dāng)Ra為103，104時(shí)，兩者的三維流場圖的差別很??；而當(dāng)Ra增加到105時(shí)，方腔中可以明顯看出2個(gè)渦流的存在。考慮到工況5方腔中的冷熱壁面的對稱布置，由圖7（c）可以看出，當(dāng)Ra＝103時(shí)，在中截面Y＝0.5處出現(xiàn)了4個(gè)對稱渦流，這同樣是因?yàn)閭鳠岱绞街饕詫?dǎo)熱為主引起的；當(dāng)Ra增加到104時(shí)，在中截面Y＝0.5內(nèi)出現(xiàn)了一個(gè)流線交匯處，這意味著流體流向X=0.5處；當(dāng)Ra=105時(shí)，中截面Y=0.5內(nèi)僅出現(xiàn)了一個(gè)流線交匯處，其原因是自然對流換熱的強(qiáng)化。中截面Y=0.5處的流線也證明了隨著Ra的增大，垂直于Y-Z平面的流動(dòng)越來越弱。

4.4不同冷熱壁面位置對平均努塞爾數(shù)的影響

表3為不同工況條件下，高溫壁面的平均努塞爾數(shù)。整體而言，平均努塞爾數(shù)Nu av隨著Ra的增大而增大，這是因?yàn)镽a表征著自然對流的強(qiáng)弱，Ra越大，換熱越強(qiáng)。當(dāng)Ra=103時(shí)，換熱方式主要是熱傳導(dǎo)，工況1?4方腔工況的Nu av差異很??；對于工況5方腔，其Nu av為1.712，相對于其他4種工況的Nu av更大，這是因?yàn)榉角坏纳舷卤砻婧屠錈岜诿嬷g的距離較大，有利于流體的流動(dòng)換熱。當(dāng)Ra=104時(shí)，不同工況下熱壁面的Nu av差異變得很明顯，因?yàn)槔錈岜诿嫖恢檬怯绊懽匀粚岬闹匾蛩亍９r1方腔的Nu av最小，其上壁面和下壁面阻礙了流體的流動(dòng)；工況5的Nu av最大。當(dāng)Ra增大到105時(shí)，不同工況的高溫壁面Nu av的差異更為顯著。工況5的Nu av仍為最大，表明其自然對流換熱能力最強(qiáng)，這是因?yàn)槔錈岜诿娴牟贾梅绞綔p小了上下壁面對流體流動(dòng)的抑制作用，使得浮升力的作用效果更為明顯。

5結(jié) 論

采用D3Q19-MRT和D3Q7-MRT雙分布函數(shù)模型，對具有局部冷熱壁面不同布置方式的三維方腔自然對流換熱進(jìn)行了數(shù)值模擬。通過分析三維方腔內(nèi)部溫度和流線的分布規(guī)律，研究了不同冷熱壁面位置和Ra（103≤Ra≤105）對三維方腔內(nèi)自然對流換熱的影響，并比較了不同工況下的平均努塞爾數(shù)Nu av，得到如下結(jié)論：

a.相同Ra條件下，冷熱壁面的相對布置方式對流動(dòng)換熱有重要影響。工況1的對流換熱能力最差，工況5的自然對流換熱能力最強(qiáng)。在相同工況條件下，隨著Ra從103逐漸增大105，自然對流能力增強(qiáng)，流線變得更為復(fù)雜。

b.在相同工況條件下，隨著Ra增大，Nu av也逐漸增加，換熱能力增強(qiáng)。在相同Ra條件下，工況5的Nu av均最大，表明其自然對流換熱能力最強(qiáng)，說明冷熱壁面的布置方式減小了上下壁面對流體流動(dòng)的抑制作用，使浮升力的作用效果更為明顯。

采用格子Boltzmann方法僅針對具有局部冷熱壁面不同布置方式的三維方腔自然對流進(jìn)行了數(shù)值模擬，后續(xù)作者將基于該模擬結(jié)果進(jìn)一步研究更為復(fù)雜條件下的自然對流換熱問題，以期為解決不同條件下的自然對流換熱問題提供一定參考。

參考文獻(xiàn)：

[1] 楊世銘， 陶文銓. 傳熱學(xué) [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效應(yīng)對方腔內(nèi)雙擴(kuò)散自然對流影響的格子 Boltzmann模擬 [J]. 空氣 動(dòng)力學(xué)學(xué)報(bào)， 2019， 37（2）： 207–215.

[4] 余端民 ， 趙明 . 圓內(nèi)開縫圓環(huán)形空間自然對流的格子 Boltzmann 模 擬 [J]. 上 海 理 工 大 學(xué) 學(xué)報(bào) ， 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].南京師范大學(xué)學(xué)報(bào)（工程技 術(shù)版）， 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.

（編輯：石 瑛）