0 引　言

由温差导致的热对流现象广泛存在于自然界和各种工程应用中，比如热对流在海洋对流、大气环流、地幔对流、太阳对流等自然现象以及在电子设备的冷却和凝固过程、核反应堆设计、晶体生长等工程应用中扮演着重要角色。Rayleigh-Bénard（RB）对流模型是研究热对流问题的经典物理模型之一，RB对流指流体层中下壁面加热、上壁面冷却导致的流动现象。该系统的控制参数主要是瑞利数Ra和普朗特数Pr，响应参数为努塞尔数Nu和雷诺数Re。当Ra较小时，流动由羽流主导，大尺度环流并不明显^[1]；随着Ra的增加，系统逐渐进入湍流状态，此时上下壁面存在很薄的温度边界层，冷热羽流在温度边界层内生成并在浮力的作用下运动并自组织形成大尺度环流。随着流动复杂性的增加，对RB对流研究的重心也转为对湍流状态下系统各物理量统计特性的研究^[2-4]。

在RB对流的研究中通常使用侧壁绝热的温度边界条件，在侧壁绝热的充分发展的RB湍流中，体区的温度近似均匀，因此侧壁的温度在大部分区域近似为上下壁面温度的平均值。实际的对流换热装置中，侧壁的导热性往往不能忽略^[5-9]。极端条件下，当侧壁导热系数相比流体很大时，可以认为侧壁是理想导热壁面，此时侧壁的温度为线性分布^[10-22]。

前人在宽度∶深度∶高度之比为6∶2∶1^[23]和6∶4∶1^[24]的三维方腔内研究了低Ra下不同侧壁温度对系统流动结构的影响。Corcione^[25]在二维腔体研究了侧壁温度等于底壁温度对流动和传热的影响。Stevens等^[6]研究了圆筒中不同侧壁温度对Nu的影响。Moon等^[26]研究了侧壁温度等于上壁面温度的立方腔体内的流动。Vasu等^[27]研究了三维方腔中侧壁等温条件或等热流条件对流动模态的影响。也有研究讨论了壁面厚度、导热率对RB对流的影响，提出了侧壁修正模型并进行了分析^{[6-8,25,28-36]}。工程应用中热对流的温度边界条件各种各样，Pandey等^[37]对方腔温度边界条件方面的研究进行了详细的总结，指出侧壁温度边界条件主要分为绝热、加热、冷却、温度线性分布四种。目前，对侧壁等温且温度为上下壁面温度平均值的对流问题研究仍然不多。

本文使用直接数值模拟（DNS）研究了 $ Pr = 5.3 $ ，宽高比 $ \varGamma = 1 $ 的二维方腔和三维方腔内的RB对流。主要研究侧壁等温且温度为上下壁面平均温度的系统（ISO算例），并与侧壁绝热的经典RB流动（RB算例）进行对比。在二维方腔中，研究发现，Ra比较小时上下壁面Nu不相等，这与系统中出现的大尺度羽流在中心区域摆动的流动结构有关；而在三维方腔中，也发现了上下壁面Nu不相等的现象。文中也讨论了侧壁等温条件对系统的热量输运（Nu）和动量输运（Re）的影响。

1 控制方程和数值方法

腔体的高度为 $ \hat{H} $ ，宽度为 $ \hat{W} $ ，宽高比 ${\varGamma } = \hat{W}/\hat{H}$ 固定为1，三维腔体深度为 $ \hat{D} $ ，且宽高比 ${\varGamma } = \hat{W}/\hat{H} = \hat{D}/\hat{H}$ 同样固定为1。上壁面温度为 ${\hat{T}}_{{\rm{C}}}$ 、下壁面温度为 ${\hat{T}}_{{\rm{H}}}$ 。在Oberbeck-Boussinesq近似下，对流介质的物性参数认为是常数，仅在浮力项中考虑密度变化且 $ \; \hat{\rho } = {\hat{\rho }}_{0}-{\hat{\rho }}_{0}\hat{\alpha }(\hat{T}-{\hat{T}}_{0}) $ ，其中 ${\hat{T}}_{0} = ({\hat{T}}_{{\rm{C}}}+{\hat{T}}_{{\rm{H}}})/2$ ， $ \;{\hat{\rho }}_{0} $ 、 $ \hat{\alpha } $ 为 $ {\hat{T}}_{0} $ 下的密度和热膨胀系数。笛卡尔坐标系中无量纲化的不可压缩Navier-Stokes方程为：

$ \nabla \cdot {{\boldsymbol{u}}} = 0 $

$ \frac{\partial \boldsymbol{u}}{\partial t}+\left(\boldsymbol{u}\cdot \nabla \right)\boldsymbol{u} = -\nabla p+\frac{1}{\sqrt{Ra/Pr}}{\nabla }^{2}\boldsymbol{u}+T\boldsymbol{k} $

$ \frac{\partial T}{\partial t}+\left({\boldsymbol{u}}\cdot \nabla \right)T = \frac{1}{\sqrt{Ra\cdot Pr}}{\nabla }^{2}T $

(1)

其中， $ \boldsymbol{u} $ 为无量纲速度矢量； $ t $ 为时间； $ p $ 为压力， $ T $ 为无量纲温度； $ \boldsymbol{k} $ 为竖直方向单位矢量； $ Pr = \hat{\nu }/\hat{\kappa } $ 为普朗特数，表征了流体的动量扩散速度与热扩散速度间的比值；无量纲参数 $Ra = \hat{\alpha }\hat{g}\Delta \hat{T}{\hat{H}}^{3}/(\hat{\nu }\hat{\kappa )}$ 为瑞利数，刻画了温差驱动湍流流动的强度，其中 $ \hat{g} $ 为重力加速度， $\mathrm{\Delta }\hat{T} = {\hat{T}}_{{\rm{H}}}-{\hat{T}}_{{\rm{C}}}$ 为上下壁面的温度差值， $ \hat{\nu } $ 为运动学黏性系数， $ \hat{\kappa } $ 为热扩散系数。无量纲化选取自由落体速度 $ \hat{U} = \sqrt{\hat{\alpha }\hat{g}\mathrm{\Delta }\hat{T}\hat{H}} $ 作为参考速度，腔体高度 $ \hat{H} $ 作为参考长度，参考时间为 $ {\hat{t}}_{f} = \hat{H}/\hat{U} $ ，温度通过 $T = (\hat{T}-{\hat{T}}_{{\rm{C}}})/\mathrm{\Delta }\hat{T}$ 进行无量纲化。在各个固壁施加无滑移和无穿透的速度边界条件。

数值求解使用课题组自主开发的lMn2d和lMn3d程序。该程序采用交错网格，将矢量放置在网格界面上，将标量放在网格中心处，避免了速度和压力解耦。空间离散使用二阶精度的中心差分，对非线性项采用Adams-Bashforth格式，对黏性项和扩散项采用Crank-Nicolson格式。壁面附近采用拉伸网格以更好地分辨边界层。在二维方腔和三维方腔内分别在 $1\times {10}^{6}\leqslant Ra\leqslant 1\times{10}^{10}$ 、 $1\times{10}^{5}\leqslant Ra\leqslant 1\times{10}^{8}$ 范围内计算了34个、20个算例。表1和表2给出了二维和三维计算采用的网格数。

2 结果与讨论 2.1 方腔中侧壁等温对Nu数的影响

热量输运规律是热对流系统中关注的重要问题。无量纲数Nu表征了对流传热的效率，其计算公式为：

$ Nu = \frac{\hat{Q}}{\hat{k}\mathrm{\Delta }\hat{T}/\hat{H}} $

(2)

其中， $ \hat{Q} $ 是单位面积的平板传递给对流介质的实际热量， $ \hat{k}\mathrm{\Delta }\hat{T}/\hat{H} $ 是热传导状态下的热通量， $ \hat{k} $ 为热传导系数。 $ \hat{Q} $ 为对流导致的热通量 $ {\hat{Q}}_{\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{v}} $ 和传导导致的热通量 $ {\hat{Q}}_{\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{d}} $ 之和，即：

$ \hat{Q} = {\hat{Q}}_{\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{v}}+{\hat{Q}}_{\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{d}} = {\left\langle{\hat{\rho }{\hat{c}}_{p}{\hat{u}}_{z}\hat{T}}\right\rangle}_{{\rm{A}}}- {\left\langle{\hat{k}{\partial }_{z}\hat{T}}\right\rangle}_{{\rm{A}}} $

(3)

则公式（2）可写为：

$ Nu = \frac{\hat{Q}}{\hat{k}\mathrm{\Delta }\hat{T}/\hat{H}} = \frac{{\left\langle{\hat{\rho }{\hat{c}}_{p}{\hat{u}}_{z}\hat{T}}\right\rangle}_{{\rm{A}}}-{\left\langle{\hat{k}{\partial }_{z}\hat{T}}\right\rangle}_{{\rm{A}}}}{\hat{k}\mathrm{\Delta }\hat{T}/\hat{H}} $

(4)

其中 $ {\hat{c}}_{p} $ 是定压比热容， $ {\hat{u}}_{z} $ 为竖直速度分量，对三维构型 ${\left\langle{\cdot }\right\rangle}_{{\rm{A}}}$ 是对平面做平均，对二维构型则是对直线做平均。我们标记壁面向流体传热Nu为正，壁面从流体吸热Nu为负，在上、下壁面处得到 $Nu = \pm {\left\langle{{\partial }_{z}T{|}_{z = \mathrm{1,0}}}\right\rangle}_{{\rm{A}}}$ ，局部 $ {Nu}_{{\rm{L}}} = $ $ \pm {\partial }_{z}T{|}_{z = 1, 0} $ ；右、左侧壁Nu = $ \pm {\left\langle{{\partial }_{y}T{|}_{y = 1,0}}\right\rangle}_{{\rm{A}}}$ ，局部 ${Nu}_{{\rm{L}}} = \pm {\partial }_{y}T{|}_{y = 1,0}$ ；三维前、后侧壁Nu= $\pm {\left\langle{{\partial }_{x}T{|}_{x = 1,0}}\right\rangle}_{{\rm{A}}}$ ，局部 $ {Nu}_{{\rm{L}}} = \pm {\partial }_{x}T{|}_{x = 1,0} $ 。我们首先关注了侧壁等温边界条件对Nu的影响。使用 $ {Nu}_{{\rm{t}}} $ 、 $ {Nu}_{{\rm{b}}} $ 、 $ {Nu}_{{\rm{l}}} $ 、 $ {Nu}_{{\rm{r}}} $ 分别表示二维ISO算例上、下、左、右壁面的时均Nu；使用 $ {Nu}_{{\rm{t}}} $ 、 $ {Nu}_{{\rm{b}}} $ 、 $ {Nu}_{1} $ 、 $ {Nu}_{2} $ 、 $ {Nu}_{3} $ 、 $ {Nu}_{4} $ 分别表示三维ISO算例上、下、左、右、前、后壁面的时均Nu；用 $ {Nu}_{\mathrm{R}\mathrm{B}} $ 表示经典RB算例的时均Nu。上壁面的Nu始终小于0，为了方便比较， $ {Nu}_{{\rm{t}}} $ 特指上板Nu的绝对值。

图1给出了二维和三维方腔各个壁面Nu、 $|{Nu}_{{\rm{b}}}/{Nu}_{{\rm{t}}}-1|$ 以及各个侧壁Nu与 ${Nu}_{{\rm{t}}}$ 的比值随Ra的变化情况。从图1（a、e）中发现 ${Nu}_{{\rm{b}}} > {Nu}_{{\rm{RB}}}$ 而且 ${Nu}_{{\rm{t}}} > {Nu}_{{\rm{RB}}}$ ，说明在二维和三维方腔中侧壁等温边界条件增强了系统的传热效率。而且发现随着Ra的增加， ${Nu}_{{\rm{b}}}$ 、 ${Nu}_{{\rm{t}}}$ 和 ${Nu}_{{\rm{RB}}}$ 之间的差别逐渐减小。这是因为当侧壁温度为上下壁面平均温度时，一部分热量可以直接通过侧壁来传递，从而绕开温度边界层，这个过程增强了系统的传热效率。而随着Ra的增加，温度边界层的厚度急剧减小，导致侧壁与温度边界层之间的传热减弱，这使得两个系统的Nu逐渐相等。

当侧壁温度固定为上下壁面温度的平均值时，由于系统存在的上下对称性，直觉告诉我们上下壁面的Nu应该相等，然而有意思的是，从图1（a、e）可以看出二维和三维方腔中上下壁面的Nu在某些比较低的Ra区间都出现了明显的偏差，即 ${Nu}_{{\rm{b}}}\ne {Nu}_{{\rm{t}}}$ ，这说明系统的上下对称性被打破，系统的流动结构可能发生了变化。

图1（b、f）分别给出了二维方腔和三维方腔中 $ {Nu}_{b} $ 和 $ {Nu}_{t} $ 之间相对误差随Ra的变化。在二维方腔中发现 $|{Nu}_{{\rm{b}}}/{Nu}_{{\rm{t}}}-1|$ 在 $Ra = 1\times{10}^{6}$ 时最大，并且随着Ra的增加不断减小，当 $Ra\geqslant 1\times{10}^{7}$ 时 $\left|{Nu}_{{\rm{b}}}/{Nu}_{{\rm{t}}}-1\right| < $ $ 0.5\%$ ，误差可以忽略不计。在三维方腔中，上下壁面Nu不相等情形发生在有限的Ra区间内。从图1（c、d）可以看出，在二维方腔中，当上下Nu发生较大偏差时，左右侧壁Nu绝对值、侧壁Nu与 ${Nu}_{{\rm{t}}}$ 比值的绝对值较大，而且侧壁Nu正负相同。从图1（g、h）可知在三维方腔中，当Ra很小，比如 $Ra = $ $ 1\times {10}^{5}$ 时，侧壁Nu绝对值、侧壁Nu与 ${Nu}_{{\rm{t}}}$ 比值的绝对值较大，而且存在两对侧壁Nu互为相反数，导致侧壁Nu之和为0，最终 ${Nu}_{{\rm{b}}} = {Nu}_{{\rm{t}}}$ ；随着Ra的增加，当 $ Ra = 7\times {10}^{5} $ 和 $Ra = 1\times{10}^{6}$ 时，侧壁Nu绝对值仍然较大，但是不再互为相反数，侧壁Nu之和不为0，而且侧壁Nu与 ${Nu}_{{\rm{t}}}$ 比值的绝对值较大，从而导致 ${Nu}_{{\rm{b}}}\ne {Nu}_{{\rm{t}}}$ ；随着Ra进一步增加，侧壁Nu绝对值迅速减小，与 ${Nu}_{{\rm{t}}}$ 的比值的绝对值也迅速减小，使得上下壁面Nu恢复相等。

2.2 方腔中侧壁等温边界条件对时均流场的影响

时间平均的流场可以反映系统的整体流动特性。为了探究图1中上下壁面时均Nu不相等的原因，在二维方腔中选取 $Ra = 1\times{10}^{6}$ 和 $ Ra = 2\times {10}^{7} $ 两个典型算例讨论了时均温度场以及侧壁Nu局部分布，在三维方腔中选取 $Ra = 1\times{10}^{5}$ 、 $ Ra = 3\times {10}^{5} $ 和 $Ra = $ $ 1\times{10}^{6}$ 三个典型算例研究了水平截面上的竖直速度分布。

图2给出了二维方腔中 $Ra = 1\times{10}^{6}$ 时ISO系统和RB系统的时均温度场和速度场，以及ISO系统时均流场左右侧壁的局部Nu分布。通过图2可以观察到，ISO系统和RB系统的时均场的共同特点是，它们均呈现出四涡模态。然而ISO系统流场上下出现不对称，而经典RB系统的时均流场则呈现出完美的对称四涡结构。由图2（c）可知，ISO系统时均流场左侧壁局部Nu( ${Nu}_{{\rm{L l}}}$ )和右侧壁局部Nu( ${Nu}_{{\rm{L r}}}$ )分布完全重合，但是 ${Nu}_{{\rm{Ll}}}$ 和 ${Nu}_{{\rm{Lr}}}$ 上下不对称，进行积分后在左右侧壁处得到 ${Nu}_{{\rm{l}}} < 0$ 且 ${Nu}_{{\rm{r}}} < 0$ 。系统存在能量的守恒，即 ${Nu}_{{\rm{b}}}-{Nu}_{{\rm{t}}}+{Nu}_{{\rm{l}}}+{Nu}_{{\rm{r}}} = 0$ ，最终导致 ${Nu}_{{\rm{b}}}\ne {Nu}_{{\rm{t}}}$ 。

图3展示了 $ Ra = 2\times {10}^{7} $ 时ISO系统和经典RB系统的时均流场以及ISO系统时均流场的左右侧壁Nu分布。在计算过程中，两类系统都没有发生大尺度结构的流动反转。通过图3（a、b）可以发现两个系统均存在一个大的环流结构，且在对角线上存在一对角涡。图3（c）中左右侧壁局部Nu呈现出反对称分布，即 ${Nu}_{{\rm{Ll}}}\left(z\right) = -{Nu}_{{\rm{Lr}}}(1-z)$ ，最终积分得到 ${Nu}_{{\rm{l}}} = -{Nu}_{{\rm{r}}}$ ，侧壁整体传热为0。由于系统能量平衡，所以上、下壁面的Nu必然相等。

图4展示了三维RB算例和ISO算例在计算的Ra区间内定常、非定常的流动状态相图。可以发现上下壁面Nu的不相等情形主要发生在侧壁等温系统刚刚由定常转为非定常流动的区间内。

Pallares等^[14]指出可以根据水平截面上竖直速度分布判断系统内的流动结构。图5是 $Ra =1\times {10}^{5}$ 、 $ Ra = 3\times {10}^{5} $ 、 $Ra = 1\times{10}^{6}$ 时三维ISO算例的时均流场在z = 0.15、0.5、0.85处的竖直速度分布，可看出系统始终存在单涡模态。观察图5（a）中z = 0.5平面内的竖直速度分布，可以发现系统的定常上升流动和下降流动关于对角面对称。随着Ra增加，系统进入非定常流动状态，图5（b、c）中z = 0.5面上的竖直速度分布表明上升流和下降流之间失去了对称性，该非对称结构导致侧壁整体的Nu之和不为0，从而上下壁面的Nu不再相等，这与前面提到的二维情形比较类似。随着Ra进一步增加，侧壁Nu迅速减小，相对于上下Nu可以忽略不计，此时上下壁面的Nu再次相等。

2.3 腔侧壁等温对Re的影响

雷诺数是热对流系统另一个重要的响应参数，可以用于刻画湍流的强度，其计算公式为：

$ Re = \frac{\sqrt{{\left\langle \hat{\boldsymbol{u}}\cdot \hat{\boldsymbol{u}} \right\rangle } _{{\rm{V}},t}}\hat{H}}{\hat{\upsilon }} $

(5)

其中， $ \hat{\boldsymbol{u}} $ 为有量纲的速度矢量， ${\left\langle{\cdot }\right\rangle}_{{\rm{V}},t}$ 表示对整个计算区域和计算时间进行平均。

我们分别使用 ${Re}_{{\rm{ISO}}}$ 和 ${Re}_{{\rm{RB}}}$ 分别表示ISO系统和经典RB系统的Re，探究侧壁等温边界条件对Re的影响。图6给出了二维和三维方腔流动中 ${Re}_{{\rm{ISO}}}$ 、 ${Re}_{{\rm{RB}}}$ 和 ${Re}_{{\rm{ISO}}}/{Re}_{{\rm{RB}}}$ 随Ra的变化。通过图6（b、d）可以发现 ${Re}_{{\rm{ISO}}}/{Re}_{{\rm{RB}}} < 1$ 几乎在所有计算Ra值区间下成立，只在二维方腔中 $ Ra = 3\times {10}^{7} $ 时大于1，而且与1仅相差0.90%，说明侧壁等温一般会减小系统Re。

2.4 二维方腔中侧壁等温对瞬时流场结构的影响

本小节进一步讨论上下壁面Nu出现不相等时对应的瞬时流动状态。角动量的符号可以表征大尺度环流的瞬时流动方向，其计算公式为：

$ L\left(t \right) = {\left\langle{-\left(z-1/2\right)v\left(\boldsymbol{x},t\right)+\left(y-1/2\right)w(\boldsymbol{x},{t})}\right\rangle}_{{\rm{V}}} $

(6)

其中 $(v,w) $ 是速度矢量的分量， ${\left\langle{\cdot }\right\rangle}_{{\rm{V}}}$ 表示对整个计算区域进行平均。

图7展示了 $Ra = 1\times{10}^{6}$ 时二维ISO系统和RB系统的角动量变化以及典型的瞬时温度场。如图7（a）所示，可以看到经典RB对流系统的角动量随时间发生频繁的正负变化，说明系统在连续不断的发生流动反转。图7（b、c、d）展示了经典RB对流系统发生反转的连续过程，可以观察到此时流动由羽流主导，瞬时流场中存在大涡结构，随着角涡的增长并发生融合，大尺度环流发生反转。图7（e）给出了这个反转过程中角动量的变化。图7（f）展示了ISO系统的角动量变化，发现角动量同样频繁发生正负变化，但是在瞬时温度场中没有发现大尺度环流结构以及流动反转，此时流动也由羽流主导，没有形成明显的大尺度环流，而是出现如图7（g～i）所示的三种典型的瞬时温度场，分别为大尺度热羽流占据中心区域并进行摆动的流场、大尺度冷羽流占据腔体中心区域并进行摆动的流场和冷热羽流无序发射的流场。

当系统处于大尺度热羽流占据中心区域并进行摆动的流动结构时，下壁面Nu小于上壁面Nu；当系统处于大尺度冷羽流占据中心区域并进行摆动的流动结构时，上壁面Nu小于下壁面Nu；当系统处于冷热羽流无序发射的流场结构时，上下壁面Nu几乎相等。由于三种流动结构持续的时间不同，导致了时均Nu和时均流场的不对称。而随着Ra的逐渐增加，系统流动加强，瞬时流场呈现出稳定的大涡结构，上述三种结构消失，此时上下壁面的Nu相等。

3 结　论

本文结合二维和三维直接数值模拟方法研究了侧壁恒温且温度为上下壁面平均值的边界条件对RB对流的影响，并与传统的侧壁绝热条件下的RB对流进行了比较。宽高比 $\varGamma$ 固定为1，普朗特数固定为Pr = 5.3，二维和三维研究的Ra区间分别为 $1\times{10}^{6}\leqslant Ra\leqslant 1\times{10}^{10}$ 、 $1\times{10}^{5}\leqslant Ra\leqslant 1\times{10}^{8}$ 。主要结论归结如下：

1）Ra较小时，在特定Ra区间内存在上下壁面Nu不相等的现象。该现象出现时，流场结构存在明显的上下不对称，导致侧壁的Nu之和并不为0，从而上下壁面的Nu不相等。随着Ra的增加，侧壁的Nu减小且可以忽略，此时上下壁面的Nu恢复成相等。

2）侧壁恒温时系统的传热效率要高于侧壁绝热情况，这是由于侧壁恒温条件下，一部分热量可以直接通过侧壁来传递，从而绕开温度边界层。

3）当采用恒温边界条件时，系统的Re要小于侧壁绝热条件下对应的Re，从而表明侧壁恒温边界条件会降低系统的流动强度。