包 芸，何建超，方明卫

(中山大学 航空航天学院， 深圳 518107)

0 引 言

在自然界和生活中，对流是一种广泛存在的现象，而Rayleigh-Bénard（RB）热对流则是从中抽象出来的经典模型之一[1-3]。RB热对流系统是一个封闭的腔体，上导板恒温冷却，下导板恒温加热，腔体内的流体受浮力驱动而运动。RB热对流系统的控制参数有三个，分别是Rayleigh数（Ra）、Prandtl数（Pr)和宽高比(Γ)，其中，Ra数反映系统上下板的无量纲温差；Pr数仅与流体的性质有关，表示流体的黏性与热扩散的相对强度。系统的整体响应参数有两个，Nusselt数（Nu）和Reynolds数（Re），分别反映系统的传热效率和湍流强度。

多年来，Ra、Pr对Nu、Re的影响一直是RB热对流的核心问题。关于这个问题，前人提出了很多关于Nu(Ra,Pr)和Re(Ra,Pr)的模型[4-8]，其中最成功的是Grossmann和Lohse提出的GL理论[9-10]。多数关于RB热对流系统的研究都是基于实验或者三维数值模拟，而GL理论是基于二维的方程发展而来的。尽管二维RB热对流在小Pr数时与三维RB热对流的结果有所不同[11]，但是仍然体现很多三维RB热对流的重要特征[11-14]。因此，近年来有很多关于二维RB对流数值模拟的研究，包括传热效率[15-17]、边界层[18-19]、大尺度环流的翻转[20-21]、能量耗散[22-23]等。在GL理论中，有关Ra数的结果得到了大量二、三维数据的验证[12,15,24-25]，但是对于Nu数和Re数对Pr数的依赖性研究相对较少[4,26-27]。在前人的研究中，Nu数与Pr数的关系分为两个阶段，第一阶段是Pr数较小时(Pr≤ 1)，Nu数随Pr数增大快速上升，在Pr≈ 3时达到最大值；第二阶段是Pr数较大时，Nu数不随Pr数变化[12,28-29]，或部分结果认为Nu～Prβ，|β| ＜ 0.1[30-32]。而Re关于Pr的关系同样在小Pr数和大Pr数的时候有所差异[29-32]。在大Pr数下，有关RB热对流系统流动特征的研究相对较少。

在RB热对流中，常见的结构有冷热羽流、大尺度环流、角涡、温度边界层等。在特定Pr数下，在宽高比为1的方腔中会出现大尺度环流，并且冷热羽流会随着大尺度环流运动，对系统的传热有重要影响[33-34]。随着Pr数的变大，系统的流态会发生变化，从明显的大尺度环流转变成羽流主导的流态[17,31-32,35]，系统的流动会变得更加紊乱。研究二维大Pr数下的流态变化和结构转变，可以更清楚地认识对大Pr数下RB热对流系统的传热与流态之间的关系。

本文计算了Pr= 10、20，Ra= 1×108～1×1012范围内的多个二维RB热对流系统算例，并与小Pr数Pr=4.3的情况进行对比，讨论大Pr数的大尺度环流形态特征，定量刻画大尺度环流的形态以及系统Re数随Ra数的变化特征。

1 湍流热对流DNS的并行直接求解

在Oberbeck-Boussinesq近似下，引入特征长度H，特征温度 ΔT、特征速度以及特征时间τ=H/U对热对流方程进行无量纲化，得到:

其中，u为 无量纲速度矢量；θ为无量纲温度；k为单位垂向矢量；p为压力；无量纲参数Ra =(βgΔTH3)/(κν)为Rayleigh数，β为热膨胀系数，g为重力加速度，ΔT为上下壁面温差，H为系统装置的高度，W为宽度，Γ=W/H反映了对流系统的几何尺寸，κ为热扩散率；Pr = ν/κ为Prandtl数，ν为运动黏性系数。数值计算中，边界条件为：壁面速度均采用无滑移条件，温度为侧壁采用绝热条件，上下底板采用恒温条件，上底板恒定低温θ= -0.5，下底板恒定高温θ= 0.5。

本文采用Parallel Direct Method of DNS（PDMDNS）方法[36]对二维高Ra数湍流热对流进行数值模拟，在“天河二号”超级计算机上完成了多组Pr数和Ra数的2D湍流热对流DNS计算。在数值模拟中，计算网格与时间步长均满足Gao等[14]提出的标准，以便充分识别最小尺度的结构和边界层，并且计算网格与Zhu等[16]的计算网格量级一致（例如Ra=5×1010，网格为2048×2 304）。每个算例的统计时间至少为200个无量纲时间，Nu的相对误差在1%以内。本文中研究的Pr数为10、20，相应的Ra数1×108～1×1012，跨度为5个量级。

2 大Pr数大尺度环流形态特性研究

在RB热对流中，系统的传热效率主要受冷热羽流的运动影响，而多数情况下冷热羽流是随着大尺度环流运动的。因此，大尺度环流的形态是影响羽流运动及系统传热的重要因素之一。

在研究羽流运动路径对传热特性的影响时，曾讨论了Pr= 0.7和4.3的情况[37-38]。羽流运动的不稳定造成大尺度环流形态从椭圆形到圆形的突变，导致传热Nu数随Ra数的变化出现转折，只是不同Pr数对应的形态突变Ra数不一样[38]。这种差别可能与角涡处的羽流运动状态相关，在较低Ra数，椭圆型大尺度环流时角涡尺寸较大可形成独立的环流，而在较高Ra数时，圆形大尺度环流时角落中羽流运动不规则[22,37-38]。因此，小Pr数的大尺度环流形态有两种典型状态：椭圆形和圆形。

与较小Pr数对比，大Pr数（Pr≥10）的系统传热效率与Ra数的关系有所差异[17,32]。根据之前的工作[17]，大Pr数的流场、羽流形态与小Pr数的有明显不同，补偿Nu数(Nu/Ra0.3)的变化趋势在大Pr数时补偿Nu数会随Ra数增大而变小，与小Pr数的相反，因而造成低Ra数时大小Pr数的传热效率相差较大。

为了研究这种在低Ra数时不同Pr数对应的传热效率差异问题，对大Pr数时的流动形态特性进行研究，重点探讨羽流运动路径的变化以及造成的大尺度环流形态改变。在流场充分发展后，通过叠加多个瞬时场得到的平均速度场能体现系统整体的速度分布。为了方便比较，所有速度场在后处理时进行了统一流动方向（逆时针）的处理。由于本文的数值模拟计算是无量纲化的，因此后处理和分析的结果中所展示的速度为无量纲速度。

图1给出的是Pr= 10时不同Ra数的平均速度场，同样给出了最大速度点及对应的大尺度环流周长线。大尺度环流周长的定义为过最大速度点的闭环流线[33]。

图中明显地看到，在低Ra数时流场形态并不是倾斜的椭圆大尺度环流和两个角涡，大尺度环流的周长线是一个略为倾斜的圆角方形，最大速度点在侧壁上。随着Ra数的增加，图1（d）中Ra= 1×109时流场流态变为倾斜的椭圆形大尺度环流加两个角涡的形态，最大速度点出现在大尺度环流和角涡相剪切的位置。随着Ra数进一步增加，如图1（g、h），大尺度环流形态会从图1（f）中的椭圆形突变为圆形。由此可见，大Pr数时二维热对流的大尺度环流出现了三种形态：方形、椭圆形和圆形。

在图1（a、b）中方形的大尺度环流，较大的速度都分布在四个壁面附近，最大速度点在边壁附近，距离上下壁面有一定距离。之所以能形成这样的图形，是因为在大Pr数低Ra数时的热对流中，从底板边界层脱落的冷热羽流聚集后能直接沿侧壁运动到对板，推动羽流也沿底板向一个方向运动，从而形成了略为倾斜的圆角方形大尺度环流。由于羽流直接沿侧壁加速运动，导致最大速度点发生在侧壁附近。

随着Ra数的增加，逐渐有角涡出现，如图1（c）所示，并且较弱的角涡与大尺度环流的相剪切处的速度开始增大，但此时的最大速度点仍在侧壁。到图1（d）中显示的大尺度环流已变成了倾斜椭圆形加角涡的形态。由此，认为大Pr数的方形大尺度环流形态变化到倾斜椭圆形形态，存在一个过渡。也即，不像大尺度环流形态从椭圆形变为圆形是一个突变，大尺度环流形态方形变为椭圆形是一个渐变。

图1 Pr = 10, Ra = 1×108～1×1012的平均速度场与大尺度环流(每张图的色标范围是0到各个算例的最大无量纲速度值，蓝色表示低速，红色表示高速；白色点表示全场的最大速度点，黑色线表示最大速度点所在的闭环流线，即大尺度环流周长线)Fig. 1 Illustration of LSC in terms of averaged velocity coutour for Pr = 10, Ra = 1×108～1×1012（The color bar start from 0 to maximum value of dimensionless velocity of fluid, where blue represents low speed value, red represents high speed value, white point represents maximum velocity value, black line represents LSC）

在Pr= 20的算例中，我们同样发现了三种大尺度环流形态。

图2给出了Pr= 20的不同Ra数的速度场图。从前面Pr= 10时对方形大尺度环流形态的描述，大尺度环流周长线为圆角方形且最大速度点在侧壁上，可以看见，图2（a～d）均为方形的大尺度环流形态，且其中图2（d）中由于弱角涡的出现为过渡形态。图2（e、f）为椭圆形大尺度环流形态，图2（g、h）为圆形大尺度环流形态。Pr= 20的大尺度环流形态也为三种，方形、椭圆形和圆形，与Pr= 10的相同。只是Pr= 20的大尺度环流形态方形到椭圆形渐变所对应的Ra数相应增高。

图2 Pr = 20, Ra = 1×108～1×1012的平均速度场与大尺度环流Fig. 2 Illustration of LSC in terms of averaged velocity coutour for Pr = 20, Ra = 1×108～1×1012

综上所述，从不同Pr数的大尺度环流形态及周长线的变化特性可以明显的看出，小Pr数对应的大尺度环流形态只有两种，椭圆形和圆形，而大Pr数对应的大尺度环流形态有三种，方形、椭圆形和圆形。

3 大Pr数大尺度环流形态几何特性研究

从上节中的速度场图可以较直观的看到大尺度环流形态的不同类型。为了定量探讨大尺度环流形态变化特性，对大尺度环流周长线的变化特性进行研究。首先讨论不同Pr数大尺度环流的周长随Ra数的变化情况。

图3分别给出了Pr= 10和20的大尺度环流周长随Ra的变化情况，其中Ra数范围为1×108～1×1012，同时给出小Pr数Pr= 4.3相应的部分结果[38]。

在小Pr时仅存在两种大尺度环流形态，从图3（a）可以看出椭圆形的大尺度环流周长在2.5～3之间，在Ra≤ 5×109时总体趋势是变大。在椭圆形突变为圆形（Ra= 1×1010）时，大尺度环流周长会突然减小到2以下，随后有个缓慢减小后又缓慢变大的走向。可以明显地见到，大尺度环流形态从椭圆形突变为圆形会产生大尺度环流周长的间断性减小。

图3 大尺度环流周长与Ra数关系Fig. 3 Relation between LSC length and Ra

对于大Pr数的情况，如图3（b、c）所示，在大尺度环流形态为方形时大尺度环流周长的值较大在3左右，会随Ra增大而略微增大。在大尺度环流渐变为椭圆形时周长会略微减小。而大尺度环流形态从椭圆形突变为圆形时也存在间断性减小。从这里可以看出，大尺度环流从方形转变为椭圆形是一种渐变，从椭圆形转变为圆形是一种突变。

大尺度环流周长的变化不能完全描述大尺度环流形态的特征。椭圆形的典型特征是其长短轴之比。由于椭圆形大尺度环流是倾斜的，从系统方腔的对角线与大尺度环流周长线的交点，定义大尺度环流系统的长短轴之比，用来反映大尺度环流形态的几何特性。方形和圆形的对角线长短轴之比会比较接近1，而椭圆形的长短轴之比会显示出一个较大的值。

图4给出了对角线长短轴的定义，两个白点之间的实线为长轴，两个白点的距离为长轴长度，两个灰点之间的实线为短轴，两个灰点的距离为短轴长度。由于前文提到对大尺度环流进行统一方向，因此系统的流动方向不会影响长短轴的计算和变化趋势。需要特别说明的是，在大尺度环流为方形和圆形时，长轴和短轴并非对应真正的长短轴，只是用于刻画大尺度环流形态的特征。

图4 长短轴定义图(白色实线表示方腔的对角线，黑色线表示大尺度环流周长线，白色圆点和灰色圆点表示对角线与大尺度环流周长线的交点，两个白色圆点的距离定义为长轴长度，两个灰色圆点的距离定义为短轴长度)Fig. 4 Illustration of definition of the ratio between the major axis and the minor axis of LSC(White lines represent diagonal lines of the cell, black line represents the LSC, white and grey points represent the intersections of the LSC and diagonal lines. The distance of two white points represents the length of the major axis and the distance of two grey points represents the length of the minor axis)

图5给出了Pr数情况下大尺度环流周长的长短轴之比随Ra数的变化，图中虚线表示长短轴之比为1.2的位置。作为小Pr数情况的对比，对Pr= 4.3的结果进行了同样处理。

图5（a）中小Pr数Pr= 4.3的情况，低Ra数的区域长短轴之比均大于1.2，表明此时的大尺度环流形态均为椭圆形，随著Ra数的增高，长短轴之比突变到1左右，随后基本保持不变，此时大尺度环流形态为圆形。

图5（b、c）为大Pr数的情况。在较低的Ra数时，长短轴之比并不完全等于1，这是因为方形的大尺度环流形态略为倾斜且具有圆角，并不完全对称与方腔的对角线，但长短轴之比小于1.2。并且由方形到椭圆形是渐变过程，所以长短轴之比最后是逐渐增大的，图中小圆圈标出了大尺度环流从方形到椭圆形的过渡点。随着Ra数的增加，椭圆形的大尺度环流形态的长短轴之比都大于1.2，而后突然减小到1附近，反映出大尺度环流形态从椭圆形到圆形的突变过程。

图5 长短轴之比的变化特性(黑色虚线表示长短轴之比为1.2，黑色圆圈表示方形到圆形之间的过渡点)Fig. 5 Relation between aspect ratio and Ra（Dash line indicates that the ratio is 1.2, black circle represents the transition point from square LSC to circular LSC）

由以上结果可以看出，与小Pr数热对流具有典型的两种流动形态情况不同，大Pr数时热对流的大尺度环流形态分三种，其中方形到椭圆形是渐变，椭圆形到圆形是突变。方形大尺度环流形态与椭圆形是不同的两种形态，除了形状不一致，还存在流动中的最大速度位置完全不同的流动特征。

4 湍流热对流Re数特性研究

大尺度环流形态的变化可能会引起热对流系统湍流强度变化规律的改变。反映系统湍流强度的物理量为Re数。不同形态的最大速度的位置变化很大，以最大速度值为特征速度的Re数的变化规律中，会存在大尺度环流形态变化的影响。给出本文Re数的定义：

图6给出了三个Pr数下Re随Ra的变化情况。在图6（a～c）中可以看到，Re随Ra的变化都存在一个明显的间断，而间断处对应的是大尺度环流形态有椭圆形突变为圆形的位置。在较低Ra数时、在大尺度环流突变为圆形之前，三个Pr数下Re与Ra存在正幂律关系，依次为Re～Ra0.58、Re～Ra0.60、Re～Ra0.59。这些结果与他人研究的二维计算结果[12,15,39]以及Li等[32]在准二维的实验中得到的幂律关系基本一致，表明本文计算结果的合理性。

图6 Re数与Ra数的关系(黑色实线分别表示Re (Ra)的标度率)Fig. 6 Dimensionless relation between Re and Ra（Black line represents the best power-law fits of Re(Ra) to the corresponding data）

为了进一步明确表示Re数与Ra数的关系变化，使用转变为圆形之前的幂律对Re数进行补偿，如图（e、f）所示。从图中可以看出，大尺度环流形态由椭圆形突变为圆形前后的Re数变化明显的分为两段，突变处存在明显的间断式下降。补偿Re数在高Ra数大尺度环流形态转变为圆形后会逐渐变小，表明转变后的幂律要小于转变前的幂律。在图6（e、f）中大Pr数较低Ra数的区域，补偿Re数都会出现一个极小值。而这个极小值对应的正是大尺度环流形态从方形转变为椭圆形的过渡点。也即，在大Pr数时低Ra数区大尺度环流形态从方形渐变为椭圆形时，也会造成补偿Re数的下降。只是相对于椭圆形突变为圆形的间断式下降，渐变引起的补偿Re数变化要小很多。

由此，大尺度环流形态的变化会对热对流的湍流强度变化规律产生影响。

5 结 论

本文采用DNS模拟计算了大Pr数Pr= 10、20，Ra= 1×108～ 1×1012的系列二维湍流热对流算例，作为对比给出了小Pr数Pr= 4.3的部分结果，研究不同Pr数情况下大尺度环流形态特性。首次采用大尺度环流周长及形态特征量长短轴之比，定量描述大尺度环流形态变化，讨论大尺度环流形态特征以及随Ra数的变化规律，以及大尺度环流形态变化与系统Re数的关系。研究结论如下：

1） 相比之前发现的小Pr数时有椭圆形和圆形两种大尺度环流形态，大Pr的热对流有一种新的方形大尺度环流形态。随着Ra数增大，大尺度环流形态的转变依次为方形、椭圆形和圆形，并且Pr变大时，转变Ra数也相应会增高。方形大尺度环流与椭圆形的最大区别在于最大速度的位置不一样，方形的最大速度点在侧壁，而椭圆的在椭圆与角涡的剪切处。

2） 大Pr数时，大尺度环流的周长随Ra数增大椭圆形时略减小，从椭圆形突变为圆形时周长间断式大幅减小，与前人在小Pr数的研究结论一致，而在方形时随Ra数增大略有增大。大尺度环流形态特征量长短轴之比在不同形态时呈现不同规律，大尺度环流为方形时，比例在1.0～1.2之间，椭圆形在1.2以上，而圆形时基本为1。这一现象在Pr= 10和20均存在。

3） 大尺度环流形态的变化会对反映热对流湍流强度的Re数随Ra数变化规律产生影响。大Pr数时，补偿Re数变化规律可明显看到，大尺度环流形态由方形渐变为椭圆形时Re数出现小量下降，过渡点出现极小值，由椭圆形突变为圆形时Re数出现间断式大幅下降。这说明了不同大尺度环流形态时，Re数的标度律会发生变化，因此在进行Re数标度律研究时，需要考虑大尺度环流形态的影响。本文仅研究了大尺度环流为方形、椭圆形的Re数标度律，对于更高Ra数时圆形大尺度环流Re数的标度律仍需有进一步研究。