自然对流作为一种重要的物理现象，具有较强的应用背景，涉及房间通风、废料冷却及冶金工程等众多领域^[1]。对于自然对流的研究主要采用实验和数值模拟两种方法。实验方法能够观察到宏观的驱动现象，但也存在实验时间长和费用高的局限性。随着计算机的快速发展，数值模拟方法在解决流体力学问题中发挥了重要作用，并已发展成为常用方法^[2]。近年来兴起的格子Boltzmann方法（lattice Boltzmann method，LBM）因其具有计算效率高、易并行和边界处理简单等优点^[3-4]被广泛应用于传热传质的研究中。因此，许多关于自然对流的研究均采用LBM。Xu等^[5]研究了在具有Soret和Dufour效应的方形二维封闭空间中，加热圆柱体周围的双扩散自然对流现象。研究表明，在高瑞利数Ra=10⁵条件下，随着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假设，其无量纲方程描述为

$ \left\{ \begin{gathered} \frac{{\partial U}}{{\partial X}} + \frac{{\partial V}}{{\partial Y}} + \frac{{\partial W}}{{\partial Z}} = 0 \\ \frac{{\partial U}}{{\partial \tau }} + U\frac{{\partial U}}{{\partial X}} + V\frac{{\partial U}}{{\partial Y}} + W\frac{{\partial U}}{{\partial Z}} = - \frac{{\partial P}}{{\partial X}} + Ma\sqrt {\frac{{Pr}}{{3Ra}}} \left( {\frac{{{\partial ^2}U}}{{\partial {X^2}}} + \frac{{{\partial ^2}U}}{{\partial {Y^2}}} + \frac{{{\partial ^2}U}}{{\partial {Z^2}}}} \right) \\ \frac{{\partial V}}{{\partial \tau }} + U\frac{{\partial V}}{{\partial X}} + V\frac{{\partial V}}{{\partial Y}} + W\frac{{\partial V}}{{\partial Z}} = - \frac{{\partial P}}{{\partial Y}} + Ma\sqrt {\frac{{Pr}}{{3Ra}}} \left( {\frac{{{\partial ^2}V}}{{\partial {X^2}}} + \frac{{{\partial ^2}V}}{{\partial {Y^2}}} + \frac{{{\partial ^2}V}}{{\partial {Z^2}}}} \right) \\ \frac{{\partial W}}{{\partial \tau }} + U\frac{{\partial W}}{{\partial X}} + V\frac{{\partial W}}{{\partial Y}} + W\frac{{\partial W}}{{\partial Z}} = - \frac{{\partial P}}{{\partial Z}} + Ma\sqrt {\frac{{Pr}}{{3Ra}}} \left( {\frac{{{\partial ^2}W}}{{\partial {X^2}}} + \frac{{{\partial ^2}W}}{{\partial {Y^2}}} + \frac{{{\partial ^2}W}}{{\partial {Z^2}}}} \right) + \frac{{M{a^2}}}{3}\varTheta \\ \frac{{\partial \varTheta }}{{\partial \tau }} + U\frac{{\partial \varTheta }}{{\partial X}} + V\frac{{\partial \varTheta }}{{\partial Y}} + W\frac{{\partial \varTheta }}{{\partial Z}} = Ma\sqrt {\frac{1}{{3RaPr}}} \left( {\frac{{{\partial ^2}\varTheta }}{{\partial {X^2}}} + \frac{{{\partial ^2}\varTheta }}{{\partial {Y^2}}} + \frac{{{\partial ^2}\varTheta }}{{\partial {Z^2}}}} \right) \\ \end{gathered} \right. $

(1)

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

$ \begin{split}& X=\frac{x}{H}\text{，}Y=\frac{y}{H}\text{，}Z=\frac{z}{H}\text{，}{u}_{{\rm{c}}}=\sqrt{g\beta \left({T}_{{\rm{h}}}-{T}_{{\rm{c}}}\right)H}\text{，}\\&Ma=\frac{{u}_{{\rm{c}}}}{{c}_{{\rm{s}}}}\text{，}U=\frac{u}{\sqrt{3}{c}_{{\rm{s}}}}\text{，}V=\frac{v}{\sqrt{3}{c}_{{\rm{s}}}},\; W=\frac{w}{\sqrt{3}{c}_{{\rm{s}}}}\text{，}\\&Pr=\frac{\upsilon }{a}\text{，}P=\frac{p}{3\rho {c}_{{\rm{s}}}^{2}}\text{，}\tau =\frac{t\sqrt{3}{c}_{{\rm{s}}}}{H}\text{，}\\&Ra=\frac{g\beta ({T}_{{\rm{h}}}-{T}_{{\rm{c}}}){H}^{3}Pr}{{\upsilon }^{2}}\text{，}\varTheta =\frac{T-{T}_{{\rm{c}}}}{{T}_{{\rm{h}}}-{T}_{{\rm{c}}}}\\[-12pt] \end{split}$

(2)

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

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

$ \begin{split} N u=&\left.\frac{h}{k / H}\right|_{Y=0, X \in[0,1], Z \in[0,1]}=\\&-\left.\frac{\partial \varTheta}{\partial Y}\right|_{Y=0, X \in[0,1], Z \in[0,1]} \end{split}$

(3)

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

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

$ N{u_{{\rm{av}}}} = \int_{{Z_1}}^{{Z_2}} {\int_0^1 {Nu{\rm{d}}X{\rm{d}}Z} },\;\; {\text{ }}{{{Z}}_{\text{1}}}{\text{,}}\;{{{Z}}_{\text{2}}} \in \left[ {0,1} \right] $

(4)

式中，Z₁，Z₂分别为上、下壁面的取值。

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

2 MRT-LBM模型

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

$ \begin{split}& {f_i}\left( {{\boldsymbol{r}} + {{\boldsymbol{e}}_i}\Delta t,t + \Delta t} \right) - {f_i}\left( {{\boldsymbol{r}},t} \right) = {{\boldsymbol{\varOmega}} _{{\text{ }}i}} + {F_i}{\text{ }},\\&\qquad \qquad i = 1,2, \cdots ,19 \end{split} $

(5)

式中：e_i为离散速度；F_i为源项；Ω_i为碰撞项。

$ \boldsymbol{e}_i=c\left[\begin{array}{ccccccccccccccccccc} 0 & 1 & -1 & 0 & 0 & 0 & 0 & 1 & -1 & 1 & -1 & 1 & -1 & 1 & -1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & -1 & 0 & 0 & 1 & 1 & -1 & -1 & 0 & 0 & 0 & 0 & 1 & -1 & 1 & -1 \\ 0 & 0 & 0 & 0 & 0 & 1 & -1 & 0 & 0 & 0 & 0 & 1 & 1 & -1 & -1 & 1 & 1 & -1 & -1 \end{array}\right] $

(6)

$ {F_{{i}}}{\text{ = }} \Delta t{\boldsymbol{G}}\frac{{({{\boldsymbol{e}}_i} - {\boldsymbol{V}})}}{p}f_{{\rm{eq}},i}({\boldsymbol{r}},t) $

(7)

$ {\boldsymbol{G}} = - \beta \left( {T - {T_0}} \right){\boldsymbol{g}}{\text{ }} $

(8)

$ \begin{split}& f_{{\rm{eq}},i}({\boldsymbol{r}},t) = \rho {\omega _i}\left[ {1 + \frac{{{{\boldsymbol{e}}_i}{\boldsymbol{V}}}}{{c_{\rm{s}}^2}} + \frac{{{{({{\boldsymbol{e}}_i}{\boldsymbol{V}})}^2}}}{{2c_{\rm{s}}^2}} - \frac{{{\boldsymbol{VV}}}}{{2c_{\rm{s}}^2}}} \right]{\text{ , }}\\&\qquad i = 1,2, \cdots ,19 \end{split} $

(9)

$ {\omega _i} = \left\{ \begin{gathered} 1/3,{\text{ }}i = 1 \\ 1/18,{\text{ }}i = 2,3, \cdots ,7 \\ 1/36,{\text{ }}i = 8,9, \cdots ,19 \\ \end{gathered} \right. $

(10)

$ \begin{split}& {{\boldsymbol{\varOmega}} _{{\text{ }}i}} = - {{\boldsymbol{M}}^{ - 1}}{\boldsymbol{S}}\left[ {{m_i}\left( {{\boldsymbol{r}},t} \right) - m_{{\rm{eq}},i}\left( {{\boldsymbol{r}},t} \right)} \right]{\text{ ,}}\\&\;\;\qquad \quad i = 1,2, \cdots ,19 \end{split}$

(11)

式中：G为有效重力项； $f_{{\rm{eq}},i}({\boldsymbol{r}},t)$ 为平衡态分布函数；V为速度；T₀为参考温度；ω_i为密度权重系数；m_i(r,t )和 $m_{{\rm{eq}},i}({\boldsymbol{r}},t)$ 分别为分布函数和平衡态分布函数；M和S分别为变换矩阵和碰撞矩阵。

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

$ \left\{ \begin{gathered} \rho {\text{ = }}\sum\limits_{i = 1}^{19} {{f_i}\left( {{\boldsymbol{r}},t} \right),{\text{ }}\quad\rho {\boldsymbol{V}} = \sum\limits_{i = 1}^{19} {{f_i}\left( {{\boldsymbol{r}},t} \right)} } \\ {\boldsymbol{m}} = {\boldsymbol{M}}{\left[ {{f_1},{\text{ }}{f_2},{\text{ }} \cdots, {f_{19}}} \right]^{\rm{T}}} \\ \end{gathered} \right. $

(12)

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

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

$ \begin{split}& {g_i}\left( {{\boldsymbol{r}} + {{\boldsymbol{u}}_i}\Delta t,t + \Delta t} \right) - {g_i}\left( {{\boldsymbol{r}},t} \right) =\\&\quad - {{\boldsymbol{N}}^{{{ - }}1}}{\boldsymbol{Q}}\left[ {{n_i}\left( {{\boldsymbol{r}},t} \right) - n_{{\rm{eq}},i}\left( {{\boldsymbol{r}},t} \right)} \right]{\text{, }}i = 1,2, \cdots ,7 \end{split} $

(13)

式中： $n_{{\rm{eq}},i}\left( {{\boldsymbol{r}},t} \right)$ 为平衡态分布函数；N和Q为能量分布的变换矩阵和碰撞矩阵； ${{\boldsymbol{u}}_i}$ 为离散速度。

$ {{\boldsymbol{u}}_i} = c\left[ \begin{gathered} 0\quad \;\; 1\quad \;\; - 1\quad \;\; 0\quad \;\; 0 \quad \;\; 0 \quad \;\; 0 \\ 0\quad \;\; 0\quad \;\; 0\quad \;\; 1\quad \;\; - 1\quad \;\; 0\quad \;\; 0 \\ 0\quad \;\; 0\quad \;\; 0\quad \;\; 0\quad \;\; 0\quad \;\;1 \quad \;\; - 1 \\ \end{gathered} \right] $

(14)

能量权重系数

$ {\omega _{T,i}} = \left\{ \begin{gathered} 1/4,{\text{ }}i = 1 \\ 1/8,{\text{ }}i = 2,3, \cdots ,7 \\ \end{gathered} \right. $

(15)

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

$ \left\{ \begin{gathered} T = \sum\limits_{i = 1}^7 {{g_i}\left( {r,t} \right)} \\ {\boldsymbol{n}} ={\boldsymbol{ N}}{\left[ {{g_1},{g_2}, \cdot \cdot \cdot {g_7}} \right]^{\rm{T}}} \\ \end{gathered} \right. $

(16)

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

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

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

3 程序验证

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

4 结果与分析 4.1 工况1和工况2的特性分析

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

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

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

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

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

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

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

4.3 工况5的特性分析

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

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

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

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

5 结　论

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

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

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

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