胡 瑾，章盛祺，夏振华，*

(1. 浙江大学 流体动力与机电系统国家重点实验室，杭州 310027；2. 浙江大学 航空航天学院，杭州 310027；3. 南方科技大学 广东省湍流基础研究与应用重点实验室，深圳 518055；4. 北京大学 湍流与复杂系统国家重点实验室，北京 100871)

0 引 言

热对流是一种常见的物理现象，与人们的日常生活、生产息息相关[1]。其不仅存在于北极夏季融冰湖[2]、大气及海洋环流[3-4]、地球内部地幔形成[5]等自然现象中；还广泛应用于核反应堆冷却[6]、新型储能材料[7]、芯片散热[8]、磁流体传热[9]等工业工程问题中。而Rayleigh-Bénard对流（Rayleigh-Bénard convection，RBC）系统，则是从众多复杂的自然现象及工程问题中抽象简化出来的研究对流问题的经典模型。同时，由于RBC系统中可以同时存在多种热量传递方式，因此它也是传热领域非常经典的课题，研究RBC系统对认识实际情况下各种各样对流过程中的传热机制有着非常重要的意义。

狭义上的RBC通常指在一个封闭腔体内，上表面冷却，下表面加热，形成恒定温差，从而导致腔体内流体产生特殊流动的系统，如图1（来源van der Poel等[10]）所示。当系统充分发展后，上下壁面附近会形成温度边界层[11]，同时不断生成冷热羽流，并在腔体中形成稳定的大尺度环流（large-scale circulation，LSC）结构。除了对经典标度律[12-15]、拟序结构[16-17]及能量耗散[18-19]等的研究，科研工作者在近二十年对LSC的反转问题也进行了深入的研究[20-22]。最近Jiang等[23]通过设置上下壁面的非对称棘轮状粗糙度和引入较小倾角对热对流系统中的LSC和传热进行控制；Zhang等[24-25]通过增加侧壁局部恒温条件也实现了对LSC反转及强弱的有效控制。

图1 RBC系统示意图[10]Fig. 1 Sketch of RBC system[10]

由于经典的RBC系统假设上下边界温度恒定且分布均匀，这与一些实际系统的热边界条件并不相符。例如在地球物理研究中模拟地幔对流时，需要考虑热传导性差异较大的大陆板块与海洋板块，此时热边界条件就不再均匀。另一方面，工程应用中一些对流换热优化问题也需要考虑非均匀热边界条件。在对流换热的优化问题研究中，已有的增加对流传热效率的途径主要分为两种：一是改变换热板表面的粗糙度，该方法已被证明可在有限Rayleigh数（Ra）范围内显著增强对流热流[26-29]；二是采用等温和绝热相结合的混合热边界条件，已有的数值和实验研究结果表明，系统的传热效率会随绝热面积的增大而减小[30]，当绝热面积超过50%时不同的绝热位置设置会对热传递产生显著差异[31]。除此之外，也有科研人员发现，在非均匀热边界条件下，当单个绝热块特征尺度与温度边界层厚度相近时，其Nusselt数（Nu）会与经典RBC系统相近[32-33]。近来，Vasiliev等[34]研究了Prandtl数Pr=6.46、Ra=1×107时仅下壁面局部加热条件对方腔内热湍流系统的流动结构及传热规律的影响，其发现加热面积相同时，LSC及形态与加热块数有较强的关联，传热效率随加热块数的增加而增大；Nandukumar等[35]研究了Pr= 1时，上下壁面均为一半加热一半冷却条件下热对流系统与经典RBC系统的关系，其发现此非均匀热边界条件在1×105＜Ra＜1×108时可明显增加系统传热效率。由此可知，研究非均匀热边界条件下热对流的流场及传热规律是更切合实际和富有价值的。

但目前已有文献中非均匀热边界条件主要考虑简单对称分布情况，所以进一步研究非对称分布的复杂热边界条件对高Rayleigh数热湍流系统流动结构及传热规律的影响是必要的。本文参考Zhang等[24-25]算例设置，基于Pr=2、Ra=1×108的经典二维RBC系统，改变其下壁面加热条件，初步研究了加热长度为方腔边长一半时，不同局部加热位置对系统LSC结构和传热效率的影响。

1 控制方程及数值设置

1.1 控制方程

假设方腔高为H，满足宽高比 Γ=1；u=(u,v)，其中u、v分别为水平方向及竖直方向速度分量；θ为温度，θl和 θu分别为下壁面加热区域及上壁面的恒定温度； Δθ=θl−θu为 上下壁面温差；给定运动黏度υ，热扩散系数κ，导热系数λ，重力加速度g，热膨胀系数β；坐标原点定义在方腔左下角。引入特征速度（自由落体速度）特征时间T=H/U、特征长度H、特征温度 Δθ，对控制方程进行无量纲化。假设系统满足Boussinesq近似，则可得无量纲控制方程组：

其中Ra和Pr为RBC系统的无量纲控制参数，其定义如下：

后文中变量如无特别说明，均已无量纲化。图2为本文计算的三种局部加热系统及RBC系统示意图。左右壁面为绝热条件，即 ∂θ/∂x= 0；上壁面恒温冷却，即 θ=0；下壁面非加热区域满足绝热条件，即∂θ/∂y= 0，加热区域满足 θ=1。

图2 局部加热系统及RBC系统计算示意图Fig. 2 Sketches of local heating and RBC cases

1.2 数值方法

本文数值模拟均采用经过部分修改的二阶有限差分法代码AFiD[10]进行计算，其中离散Poisson方程通过水平方向的离散余弦变换解耦并采用三对角求解器进行求解，时间推进采用二阶显式Adams-Bashforth格式实现。与de Vahl Davis等[36]的方腔自然对流结果对比，文献在Pr= 0.71、Ra= 1×104、1×105、1×106条件下的壁面时均Nu数为2.243、4.519、8.800，而本文所用程序计算的结果为2.256、4.537、8.839；与Zhang等[25]的经典RBC系统结果对比，文献在Pr=2、Ra= 1×108条件下的壁面时均Nu数为24.90，而本文所用程序计算的结果为24.92。可以认为与已有文献结果基本一致，证明了现有代码的准确性。

1.3 算例设置

本文在Pr= 2、Ra= 1×108条件下，以加热左端点无量纲坐标X为位置参数，对全局加热（经典RBC）系统和局部加热系统（加热位置分别为X= 0、0.125、0.25，加热长度l=0.5）进行统一的分析和比较。

所有算例中，水平方向为均匀网格，竖直方向为非均匀网格（靠近壁面处细化），保证边界层内网格大小Δ＜0.6min[ηK,ηB][37]；时间步长足够小，足以确保Courant-Friedrichs-Lewy数小于0.2；达到统计稳态后计算足够长时间，RBC系统可观察到反转现象。文中 〈···〉表示对括号内变量求平均，下标表示如何平均。

2 计算结果及分析

2.1 局部加热条件对流场变化的影响

首先考虑不同位置局部加热条件对热湍流达到统计稳态后的流动影响，分析总动能的时间演化情况。本文计算总动能为：

其中下角标S表示对整个计算域面积求平均。通过计算可得不同加热位置及经典RBC的Ek(t)在达到统计稳态后的变化情况，如图3所示。

图3 不同加热位置及经典RBC的Ek(t)变化情况Fig. 3 Time series of Ek(t) of local heating cases and traditional RBC case

其中RBC系统的总动能均值 〈Ek〉t=0.0200，最大值和最小值之间相差至少0.01，明显比局部加热系统的均值及振幅要大。但不同位置的局部加热系统之间的均值差异并不明显，因此我们进一步计算X=0、X=0.125、X=0.25三 个系统的 〈Ek〉t，可得其值分别为0.012 5、0.013 9、0.014 5，约 为RBC系统的62.5% ～72.5%。由于局部加热系统三条曲线重合率较高，故由加热位置从左至右分别给出系统Ek(t)的方差，结果为7.56×10−4、6.67×10−4、5.82×10−4，依次减小。综上所述，我们可以认为，局部加热系统中加热位置越靠下壁面中心，系统总动能Ek(t)越大、其对应振幅越小。

根据前文所提，传统RBC具有的流场结构特性之一便是当流动充分发展后，流场内会形成较为稳定的LSC结构且会发生反转现象。为了分析局部加热是否会影响LSC的反转，本文采用角动量符号对此进行判断。计算平均角动量定义如下[38]：

图4给出局部加热系统及RBC的L(t) 在1×104≤t≤4×104范围内的演化情况，其中L<0表示LSC呈顺时针状态，L>0为逆时针。从图4中可以明显地看到经典RBC存在大涡反转现象，而局部加热系统均未出现反转情况且LSC呈顺时针状态。放大局部加热系统的L(t)演化情况，如图5所示。结合图5，进一步计算局部加热系统达到统计稳态后角动量均值，可得加热位置从左壁面至下壁面中心 〈L〉t分别为−0.049 9、−0.052 6、−0.054 0，方差Lrms分别为0.002 4、0.001 8、0.001 5，由此可知局部加热系统中加热位置越靠近下壁面中心，其L(t)绝对值越大、振幅越小，系统主体LSC稳定性越高，LSC反转可能性越低。

图4 不同加热位置及经典RBC的L(t)变化情况Fig. 4 Time series of L(t) of local heating cases and traditional RBC case

图5 不同加热位置的L(t)变化情况Fig. 5 Time series of L(t) of local heating cases

Sugiyama等[38]曾描述了角涡在LSC反转过程中的作用，即由于角落边界层分离羽流的能量补给，小的角涡会逐渐增大，直至达到腔体高度的一半，然后破坏主体LSC结构，并在相反方向上建立另一个新的LSC。由此我们知道反转现象的产生与角涡的生长密不可分，为了观察角涡及中心大涡的变化情况，绘制了流场云图及动画（在网刊资源附件中提供）。从动画及图6中可看出无论是否局部加热，其基本流场结构均主要由中心大涡及角涡组成。其中局部加热系统中下壁面角涡不稳定，且大小受限。虽然RBC系统中的角涡关于方腔中心的对称性更强，θ、u、v的数值分布均具有更高的对称性，但其中心大涡也被更剧烈地挤压，更容易诱导涡的破碎和反转。而在局部加热系统演化过程中，系统热羽流只能沿某一（距离加热区域更近的）侧壁上升，遇到下降的冷羽流时发生分离，在上壁面附近形成与RBC系统中相似的正常形态角涡，而在另一侧由于缺少上升的热羽流，只能在下壁面附近产生仅由冷羽流诱导的角涡。此时系统中下壁面角涡由于缺乏能量补给条件，显然无法正常生长，使得系统LSC无法发生反转。我们猜测，本文研究的局部加热系统里由于下壁面的角涡缺乏持续的冷热羽流相互竞争机制，导致该角涡大小无法突破LSC反转的阈值，即达到腔体高度的一半，从而无法发生LSC反转，这与Sugiyama等[38]描述的机制一致。此外，局部加热系统中角涡的大小虽然受限制，但仍然会随着时间演化不断波动。

图6 不同加热位置及经典RBC瞬时温度云图Fig. 6 Temperature contours of local heating cases and traditional RBC case

2.2 不同位置的局部加热条件对传热效率的影响

由于左右壁面均为绝热条件，故只考虑竖直方向传热情况，通过Nu数对其进行量化，本文定义有效Nu数为：

其中下边界绝热壁面 ∂θ/∂y=0，由于进出热流的平衡，满足下边界加热区域平均热流密度为上边界的两倍，故根据上述公式仍能得到可靠的系统平均Nu数。图7给出了系统达到统计稳态后，不同加热位置及经典RBC系统的Nu(t)随时间演化情况。为了更直观地展示和对比，我们将RBC系统的Nu(t)上移4，而将X=0.125和X=0.25的Nu(t)下移了8和16。

图7 不同加热位置及经典RBC的上下壁面Nu(t)变化情况Fig. 7 Time series of Nu(t) at the bottom wall (a) and top wall (b)from different cases

其 中RBC系 统 〈Nu〉t=24.92，如 上 文 所 提，与Zhang等[25]文献给出的24.90基本一致。从图中可以明显看出的是，RBC系统的 〈Nu〉t比局部加热系统振荡更为剧烈，而局部加热系统内部则是随着加热位置越靠近中心而越小。为了得到定量的数据，我们计算X=0、X=0.125、X=0.25三个局部加热系统在统计稳态的 〈Nu〉t，结果分别为16.97、17.46、18.23。可以发现，虽然局部加热系统 〈Nu〉t明显小于RBC系统，但随着加热位置越靠近下壁面中心， 〈Nu〉t的值也会随之越来越大，传热效率依次升高，这与前文总动能 〈Ek〉t、角动量绝对值 |〈L〉t|的相对大小保持一致。其中虽然局部加热系统的加热长度l=0.5，仅为RBC系统的一半，但其 〈Nu〉t可以达到传统RBC的68.1% ～ 73.2%，远高于设想的50%。而且通过比较分析，我们发现X=0.25系 统 〈Nu〉t比X=0.125系 统 增 长 了4.4%，比X=0系统增长了7.4%。也就是说，在局部加热区域长度确定的情况下，除了Bakhuis等[33]通过将导热及绝热区域分布为数量尽量多的条纹图案以提升系统传热效率外，我们也可以通过控制整体加热位置，来使系统传热效率最大化。

观察图7（b）并与图7（a）作对比，可以看出在达到统计稳态后，RBC系统上下壁面Nu(t)特征相似；在局部加热系统中，虽然上壁面Nu(t)振幅区别不大，但加热位置越靠近下壁面中心，下壁面Nu(t)幅值越小导致上下壁面之间差距越明显。这是因为上壁面附近以侧壁分离出的冷热羽流共同作用为主，对其影响相近。而另一侧由于没有热羽流的阻碍，冷羽流下降后会先到达下壁面而后横向运动，再与相遇的热羽流一同上升。结合云图及动画可知，加热位置越靠近中心，即加热区域离冷羽流到达下边界位置越近，受其冲击作用和横向运动影响的范围就越大，从而导致自由生长的热羽流数量越少，因此加热区域温度脉动就越小，Nu(t)振幅也越小。

图8展示了统计稳态时，下壁面加热区域的时均温度梯度分布，其证明加热位置越靠近冷羽流冲击点，其受冷羽流影响端温度梯度就越大，即冷羽流作用程度越强。而中心加热系统中，加热位置最靠近冷羽流发生变向运动区域，其受影响最大，热羽流生成数量始终维持在1～2个，因此热羽流可以以近似一簇的状态快速从左侧壁上升进而与冷羽流进行相互作用，又因为其对于方腔中心而言具有更稳定的LSC结构，因此具有更高效率的对流换热，故整体传热效率更高。

图8 不同加热位置系统的下壁面时均温度梯度Fig. 8 Temperature gradient of the bottom wall of different heating position cases

3 结 论

本文以二维方腔RBC系统为基础，通过直接数值模拟计算了Pr=2、Ra=1×108条件下，加热区域长度l=0.5，加热位置分别为X= 0、0.125、 0.25的三种局部加热系统。分别研究了不同位置的局部加热条件对系统总动能、角动量、LSC反转及传热效率的影响。对结果进行分析，得出以下结论：

1）局部加热系统中，加热位置越靠近下壁面中心，其总动能和角动量绝对值越大、振幅越小，系统主体LSC稳定性越高；

2）局部加热条件会对系统下壁面角涡的生长产生限制，进而抑制LSC反转；

3）虽然加热长度为RBC系统的一半，但局部加热系统总动能〈Ek〉t可达到后者的62.5%～72.5%，传热效率可达到68.1%～73.2%；

4）在一定的加热长度下，可以通过调整加热位置令传热效率最大化，其中X=0.25系 统就比X=0系统增长了7.4%；

5）在没有热羽流上升侧，冷羽流下降后到达壁面处改为横向运动，该位置与加热区域越近则对热羽流数量和对应温度脉动影响越明显，其对应端点温度梯度也越大。

由于本文计算选取局部加热位置数量有限，因此想要得到更为确切的定量关系，还需大量不同加热位置及加热长度的计算结果。本文初步证明了加热长度确定的情况下，可以通过加热位置来控制传热效率。关于局部加热条件对LSC反转的影响，后续还会深入研究。除此之外，本文作为对非对称局部加热问题的初步研究，目前仅开展了二维直接数值模拟研究，后续我们将进一步开展三维问题的直接数值模拟与分析，以及相应的实验研究。