侧壁效应对Rayleigh-Bénard对流影响的数值模拟研究

刘苗苗，王 启，万振华，孙德军

(中国科学技术大学 近代力学系，合肥 230027)

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，宽高比 Γ=1的二维方腔和三维方腔内的RB对流。主要研究侧壁等温且温度为上下壁面平均温度的系统（ISO算例），并与侧壁绝热的经典RB流动（RB算例）进行对比。在二维方腔中，研究发现，Ra比较小时上下壁面Nu不相等，这与系统中出现的大尺度羽流在中心区域摆动的流动结构有关；而在三维方腔中，也发现了上下壁面Nu不相等的现象。文中也讨论了侧壁等温条件对系统的热量输运（Nu）和动量输运（Re）的影响。

1 控制方程和数值方法

腔体的高度为， 宽度为，宽高比固定为1，三维腔体深度为， 且宽高比同样固定为1。上壁面温度为、 下壁面温度为。在Oberbeck-Boussinesq近似下，对流介质的物性参数认为是常数，仅在浮力项中考虑密度变化且其中为下的密度和热膨胀系数。笛卡尔坐标系中无量纲化的不可压缩Navier-Stokes方程为：

其中，u为 无量纲速度矢量；t为 时间；p为 压力，T为无量纲温度；k为竖直方向单位矢量；Pr=为普朗特数，表征了流体的动量扩散速度与热扩散速度间的比值；无量纲参数为瑞利数，刻画了温差驱动湍流流动的强度，其中为重力加速度，为上下壁面的温度差值，为运动学黏性系数，为热扩散系数。无量纲化选取自由落体速度作为参考速度，腔体高度作为参考长度，参考时间为，温度通过进行无量纲化。在各个固壁施加无滑移和无穿透的速度边界条件。

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

表1 不同Ra下二维数值模拟的网格数Table 1 Grid resolutions for two-dimensional simulations with different Ra

表2 不同Ra下三维数值模拟的网格数Table 2 Grid resolutions for three-dimensional simulations with different Ra

2 结果与讨论

2.1 方腔中侧壁等温对Nu数的影响

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

其中，是单位面积的平板传递给对流介质的实际热量，是热传导状态下的热通量，为热传导系数。为 对流导致的热通量和传导导致的热通量之和，即：

则公式（2）可写为：

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

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

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

图1 二维和三维方腔中各Nu随Ra的变化以及Nu之间关系随Ra的变化Fig. 1 Nu numbers and the relationships between them as functions of Ra in two- and three-dimensional square cavities

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

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

Ra=1×106

图2给出了二维方腔中 时ISO系统和RB系统的时均温度场和速度场，以及ISO系统时均流场左右侧壁的局部Nu分布。通过图2可以观察到，ISO系统和RB系统的时均场的共同特点是，它们均呈现出四涡模态。然而ISO系统流场上下出现不对称，而经典RB系统的时均流场则呈现出完美的对称四涡结构。由图2（c）可知，ISO系统时均流场左侧壁局部Nu(NuLl) 和右侧壁局部Nu(NuLr)分布完全重合，但是NuLl和NuLr上下不对称，进行积分后在左右侧壁处得到Nul<0且Nur<0。系统存在能量的守恒，即Nub−Nut+Nul+Nur=0 ，最终导致Nub≠Nut。

图2 二维数值模拟中R a=1×106时的时均温度场分析（图中黑色箭头代表速度矢量）Fig. 2 Analysis of time-averaged temperature field when R a=1×106 in two-dimensional simulation(Black arrows represent velocity vectors)

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

图3 二维数值模拟中R a=2×107 时的时均温度场分析Fig. 3 Analysis of time-averaged temperature field when R a=2×107 in two-dimensional simulation

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

图4 三维方腔内RB对流的流动状态相图Fig. 4 Phase diagram of flow states in three-dimensional RB convection

Pallares等[14]指出可以根据水平截面上竖直速度分布判断系统内的流动结构。图5是Ra=1×105、Ra=3×105、Ra=1×106时三维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再次相等。

图5 侧壁等温边界条件的三维RB对流时均竖直方向速度场在水平截面高度为z = 0.15，z = 0.5和z = 0.85处的分布Fig. 5 Time-avraged vertical velocity for the three-dimensional RB convection with isothermal sidewalls at z = 0.15, 0.5, 0.85

2.3 腔侧壁等温对Re的影响

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

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

图6 Re和R eISO/ReRB随着Ra的变化Fig. 6 Re and R eISO/ReRB as functions of Ra

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

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

其中 (v,w)是 速度矢量的分量，表示对整个计算区域进行平均。

Ra=1×106

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

图7 二维计算中R a=1×106时RB算例和ISO算例的角动量随时间的演化以及典型瞬时温度场，其中（b～d）为RB算例发生反转时不同时刻的瞬时流场，（g～i）为ISO算例三个不同时刻的典型瞬时温度场Fig. 7 Evolution of the angular momentum and typical instantaneous temperature fields for the adiabatic-sidewall case and the isothermal-sidewall case with R a=1×106 in 2D simulations. Figs. (b～d) instantaneous temperature fields during a flow reversal process for the adiabatic-sidewall case. Figs. (g～i) three typical instantaneous temperature fields for the isothermal-sidewall case

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

3 结 论

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

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

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

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