二维湍流热对流羽流运动路径对传热特性的影响*

包芸 何建超 高振源

(中山大学航空航天学院,广州 510275)

1 引 言

自然界和工业设计中存在广泛的热对流现象,热对流的传热特性研究有着重要的意义.Rayleigh-Bénard(RB)热对流是热对流研究的典型物理模型之一,即在一个封闭的空间内下底板加热上底板冷却产生热对流运动和热输运的系统[1].RB热对流系统存在丰富而复杂的流动和热输运现象,一直受到国内外学者的关注和研究.

国内外关于RB热对流的研究成果非常丰富,主要集中在传热特性变化规律以及温度边界层特性等问题.在理论研究方面,Grossman和Lohse[2,3]提出的GL理论最为成功,可以准确预测和描述大范围内传热Nusselt(Nu)数与Rayleigh(Ra)数和Prandtl(Pr)数的关系.实验研究中,He等[4]进行了Pr=0.8极高Ra数的系列实验,发现当Ra数接近1015时,反映系统参数变化关系的标度律有所增大,与GL理论预测的极高Ra数的结果一致.Zhou和Xia[5]对RB系统的温度边界层进行研究,发现温度边界层分布与Prandtl-Blasius理论一致.Stevens等[6]通过直接数值模拟(DNS)方法对三维圆柱RB热对流系统进行了模拟,Ra数达到2×1012.Zhu等[7]将二维湍流热对流的DNS模拟Ra数提高到Ra=1014,并发现羽流发射特性在极高Ra数发生变化,引起传热增强.Stevens等[8]研究了湍流热对流中热分布的超结构.黄茂静和包芸[9]计算了二维和三维的情况,发现二维热对流和三维展向平均热对流的时间平均场中都存在大尺度环流和角涡,两者温度边界层厚度关于Ra数的变化存在基本一致标度率关系.包芸等[10]讨论了Pr数对二维湍流热对流流动和传热特性的影响.Bao等[11]和Chen等[12]研究了隔板对流系统,仅在对流槽中加入一定数量顶端留有狭缝的竖直隔板,结果发现可以成倍增强传热效率.林泽鹏和包芸[13,14]对隔板对流系统的几何参数、压力特性等展开了系统的研究.Shishkina 等[15]考虑湍流脉动的影响,得到了湍流RB热对流温度边界层方程.何鹏等[16]利用湍流温度边界层方程解,与不同Ra数和Pr数下二维DNS数值解温度边界层拟合,得到了很好的效果,并对拟合参数特性进行了讨论.Gao等[17]在研究湍动能沿纵向分布时发现,二维和三维湍动能的差异仅发生在温度边界层内部.

Van der Poel等[18]对比了二维和三维热对流DNS计算结果,二维计算 R a≤1010,发现三维结果与GL理论的预测吻合良好.由于二维的Nu数值均小于三维的结果,二维Nu数在Ra < 1010的范围可以与向下平移的GL理论预测倍数线相符合,表明二维Nu数的变化规律与GL理论预测有一致性.但进一步增高Ra数时二维的传热结果偏离了GL理论预测倍数线.本文大范围Ra数二维湍流热对流的DNS计算结果研究发现,当Ra数进一步提高,传热Nu数随Ra数的变化规律出现新的转折,又回到了另一条GL理论预测倍数线上,并且这种传热特性的转折变化规律与大尺度环流羽流运动路径周长的变化特性相关联.

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

在Oberbeck-Boussinesq近似下,无量纲化的热对流方程为

其中V为无量纲速度矢量,θ为无量纲温度,p为压力,k为垂向单位矢量.无量纲参数Rayleigh(Ra)数为瑞利数,表征浮力驱动力与阻碍运动的力两者相对大小.Prandtl(Pr)数为普朗特数,表征流体黏性耗散和热耗散的关系.数值计算中边界条件为壁面速度均采用无滑移条件,温度为侧壁采用绝热条件,上下底板采用恒温条件,上底板冷却θ=—0.5,下底板加热θ=0.5.

二维热对流数值求解过程采用投影法.时间和空间均采用二阶格式,显式计算动量方程和温度对流扩散方程,压力泊松方程则需要全流场联立求解.利用FFT变换解耦泊松方程,而后求解三对角方程组所建立的直接解法,在大规模湍流热对流计算中必须通过并行计算进行.利用压力泊松方程的高效并行直接求解PDD算法,联合其他显示容易并行计算的动量方程等,建立了热对流DNS的并行直接求解方法(PDM-DNS,parallel direct method of DNS),并具有很好的并行效率[19].

3 高Ra数2D湍流热对流传热特性

热对流系统研究的核心问题是由热浮力引起的湍流流动传输热量的能力,用系统整体传热Nusselt数(Nu)表示.湍流热对流传热特性最有成效的理论研究成果是GL理论[2,3],GL理论很好地预测了传热Nu数随Ra数和Pr数的变化特性.Stevens等[20]通过近年大量新的实验和计算数据,重新对GL理论结果中的参数进行了修正,使得在很大的Ra数范围内Nu数对Ra数和Pr数的变化关系都能很好地满足理论预测结果.在Ra数的大范围中,传热Nu数随Ra数的变化标度律并不是常值,是随Ra数的变化而变化的.

当Ra≥1012以后,无论是实验还是DNS数值模拟研究都碰到过极大的困难,困扰研究的发展.因此一般称当Ra≥1012时为极高Ra情况.实验研究比DNS模拟较早突破并进行极高Ra数湍流热对流研究[4].Zhu等[7]在改进计算技术后将并行计算效率提高数十倍,从而在2018年实现了Ra=1014的二维湍流热对流DNS模拟,所花计算资源超过500万核时.

本文采用新的计算方法PDM-DNS,在“天河二号”超级计算机上进行了系列Ra数的二维湍流热对流DNS模拟,对应两组Pr=0.7和4.3.计算从Ra=107开始,最高Ra数达到Ra=1013.两组计算的Ra数如表1所列.

表1 计算Ra数Table 1. The Ra numbers.

大量的计算结果都发现,二维热对流的传热Nu数在同样的Ra数要小于三维的结果,但具有与三维结果随Ra数变化规律一致的特性.因此在探讨二维湍流热对流的传热特性时,用GL理论预测值乘以小数使其向下平移的倍数线,来反映二维传热特性的变化规律[18].

图1中给出了本文计算的两组Pr数情况下的二维湍流热对流的传热Nu随Ra数的变化情况,同时给出了四个三维计算结果.图中显示的是Nu/Ra0.3随Ra数的变化,其中黑线是GL理论预测线,虚线是GL理论预测倍数线.红色圆点是Pr=0.7的结果,蓝色三角是Pr=4.3的计算结果.图中还对应给出了Van der Poel等[18]和Zhang等[21]的二维计算结果以及Stevens等[6]的三维计算结果.

从图1中可以看到,对Pr=0.7的二维结果,原本随Ra数变化而增大的 Nu/Ra0.3在 Ra ＞109时开始减小,到Ra≈1010时达到最小,而后随着Ra数的增加,Nu/Ra0.3随Ra数的变化又开始相应地变大,回到另一条GL理论预测倍数线上,变化过程中出现两个转折.Pr=4.3的计算结果也出现同样的变化特征,只是比Pr=0.7时对应的Ra数后移一个量级.在Ra ≤ 1010的范围内Pr=4.3的Nu数与Ra数之间有较好的0.3标度律关系,当 Ra ＞ 1010后与 Van der Poel等[18]的二维结果一样,Nu/Ra0.3随Ra数的变化向下减小,但减小的幅度并没有随Ra数的增加而一直增大,当 Ra ≥ 2×1011时 Nu/Ra0.3的值又开始增加,并与GL理论预测倍数线变化走向一致.

荷里路德宫的一些建筑是三层的，一层是长长的连廊，二层和三层是房间。我惊奇地发现，宫殿的墙壁上居然有三种不同的装饰柱式，爸爸说，一层为多立克柱式，二层为爱奥尼柱式，三层则为科林斯柱式。在同一栋建筑上同时呈现出三种柱式，这也算是对古希腊建筑艺术风格的致敬吧！

图1 热对流传热Nu数随Ra数的变化情况Fig.1.The variation of the Nu number with Ra number for the turbulent convection.

图1中二维湍流热对流传热特性的这种在一定的Ra数范围,Nu/Ra0.3随Ra数变化的标度律相应减小而后增加的特殊情况,可能与二维湍流热对流的流动特性有关.为此,我们从分析二维湍流热对流的流场特性出发,探讨二维湍流热对流传热特性的特殊变化与二维湍流热对流流动特征之间的相关关系.

4 高Ra数2D湍流热对流流动特性

先以Pr=0.7为例,给出四个典型Ra数的二维湍流热对流的瞬时温度场以及温度和速度的平均场,讨论温度羽流在不同Ra数的变化情况.

瞬时温度场中的温度分布可以反映冷热羽流的状态.图2给出了Pr=0.7时四个典型Ra数瞬时温度场,可以看到在不同Ra数时羽流的形态变化.给出的四个典型Ra数都已不存在较稳定的大尺度环流和角涡.在Ra=109,角涡脱离角区随大尺度环流绕行,把羽流拉成条状并随机运动,并将部分羽流挤到方腔的中部,使得流动的最大速度不再沿侧壁而向方腔中部偏移.Ra=1010时仍以条状羽流为主,出现个别的涡状羽流.到Ra=1011和极高Ra数Ra=1012,羽流基本上全部变为小尺寸的涡状羽流,并相互缠绕,随大尺度环流做宏观的绕行.Ra数越高,涡状羽流尺寸越小,数量越多,在流动中绕行的时间越长.

平均场特性将滤掉羽流运动的瞬时影响,给出总体传热特性和运动特征.图3是上面四个典型Ra数的平均温度场和流线图.虽然四个不同Ra数的瞬时羽流形态和运动不一样,但它们的平均场分布较为相似.近底板的温度分布随Ra数的增高而越来越接近底板,表明温度边界层越来越薄.流线图形状基本一致,平均场流动都接近圆形,有四个小角涡存在.但平均场流线图并没有反映出平均速度场的大小分布和变化情况.

图4给出了对应的平均场的速度值分布.四个典型Ra数的平均场的速度值分布都基本为圆形,但速度值的最大值位置在不同Ra数时是不同的.Ra=1010速度值较大区域分布明显向方腔中部聚拢,速度的最大值位置随Ra数的变化是先向中部聚集再向外部扩展的过程.二维高Ra数的平均场速度分布的变化过程与瞬时羽流的运动形态相关联.高Ra数时的角涡脱离角区,将羽流挤入方腔的中部,使得由温度加速的流动最大速度分布区域也向方腔中部聚集.当Ra更高时羽流形态变化为小尺寸旋涡团状羽流,团状羽流随大尺度环流绕行不再挤压羽流运动,最大速度分布再次向方腔四周扩展.

图2 Pr=0.7典型Ra数的二维热对流瞬时温度场Fig.2.The instantaneous temperature fields with typical Ra number at Pr=0.7.

图3 Pr=0.7典型Ra数的二维热对流平均温度场和流线图Fig.3.The average temperature fields and stream lines with typical Ra number at Pr=0.7.

图4 Pr=0.7典型Ra数的平均速度场分布Fig.4.The average velocity field distributions with typical Ra number at Pr=0.7.

图5 Ra=1012时不同Pr数的瞬时温度场及羽流结构Fig.5.The instantaneous temperature field and plume structure with different Pr number at Ra=1012.

5 高Ra数2D湍流热对流大尺度环流路径与传热Nu数之间的相关性

定义大尺度环流路径,用于描述2D湍流热对流的平均场速度羽流运动路径特性.在平均场速度分布中找到速度值最大的那点,画出过这点的流线图,并以此流线反映大尺度环流路径[22].

图6给出了两个典型热对流的大尺度环流路径.可以看到,热对流大尺度环流路径在Ra=108是椭圆封闭曲线,路径较长,而在Ra=1011是一个近似圆形的封闭流线,路径较短.随不同Ra数大尺度环流的路径周长CLSC发生变化.

计算Pr=0.7时不同Ra数的大尺度环流路径周长CLSC,并给出它们随Ra数的变化,如图7中上半部分的红色圆点所示.图中在Ra数较低时椭圆形大尺度环流路径周长CLSC值较大.在Ra ≈109时有一个突然的减小,此时对应发生的是流态中角涡脱落羽流被挤向方腔的中心,大尺度环流变为圆形.在Ra ≈ 109至Ra ≈ 1010的区域大尺度环流路径周长CLSC值较小.随着Ra数的增加,羽流形态发生变化形成团状羽流,大尺度环流路径周长CLSC逐渐增大,到Ra=1013.

图6 Pr=0.7时大尺度环流路径图Fig.6.The large scale circulation path at Pr=0.7.

图7 大尺度环流路径周长与传热Nu数随Ra变化及其相关性,Pr=0.7Fig.7.The variation of large scale circulation path length and Nu number with Ra and its correlation at Pr=0.7.

讨论Nu数随Ra数的变化与大尺度环流路径变化的关联,将图2中Pr=0.7时二维湍流热对流的传热Nu/Ra0.3数随Ra数的变化局部放大并入图7中的下部.图中可以看到两者具有很好的相关性.Nu/Ra0.3数和大尺度环流路径周长CLSC随Ra数的变化都存在两个明显的转折点.Nu/Ra0.3数随Ra数的变化有个先增加再下降而后又增加的过程.当大尺度环流路径周长CLSC在Ra ≈109突然变小时,Nu/Ra0.3数的变化出现减小的转折.随着Ra数的增加到Ra ≈ 1010,大尺度环流路径周长CLSC开始增大,Nu/Ra0.3数的变化出现第二个转折并开始增加,而且与GL理论预测0.72的倍数线符合良好.这一结果表明当Ra数很高时,二维湍流热对流的传热特性变化规律仍可以与GL理论预测变化趋势一致.

图8中给出了Pr=4.3的二维湍流热对流大尺度环流路径周长CLSC随Ra数的变化,同样可以看到类似于图7中Pr=0.7的结果,存在明显的两个转折.

图8 大尺度环流路径周长与传热Nu数随Ra变化及其相关性,Pr=4.3Fig.8.The variation of large scale circulation path and Nu number with Ra and its correlation at Pr=4.3.

图8中同时给出了Van der Poel等[18]的二维计算结果,计算Pr=4.38,最大Ra=1011.Nu/Ra0.3数随Ra数的变化局部放大后可以明显的看到,当Ra数增高至Ra=1010时出现第一个转折,二维的结果偏离了GL理论预测倍数线,使得高Ra二维的计算结果不再与GL理论预测变化趋势一致,而且偏离情况随Ra数越来越大,因此他们认为二维的计算结果在高Ra数时可能没有意义.但他们没有继续提高Ra数的计算,错过第二个转折点的发现.

本文Pr=4.3的二维湍流热对流DNS结果表明,在一定的高Ra数范围内2D计算结果Nu/Ra0.3数随Ra数的变化会出现向下偏离GL理论预测倍数线的现象.但同时可以看到,这个范围的大尺度环流路径周长也是减小的.继续提高计算Ra数,当Ra ＞ 2X1011后大尺度环流路径周长开始增大,Nu/Ra0.3数的变化出现第二个转折,并且与GL理论预测0.7倍数线符合良好,又恢复到与GL理论预测变化基本一致.

Pr=4.3的两个转折点对应的Ra数与Pr=0.7时明显增大,分别为 Ra ≈ 5X109和Ra ≈2X1011.这个变化与流动的湍流特性有关,Pr数越小流动的湍流程度越强[9],因此对于较小Pr数,热对流流动更加容易发生角涡脱落和出现团状羽流.不同Pr数二维的结果在第二个转折点后相符合的GL理论预测倍数线略有不同,Pr=0.7时为0.72×GL理论预测值,Pr=4.3时为 0.70×GL理论预测值.Pr数较大GL理论预测倍数线略为下移.

以上研究结果可见,二维湍流热对流DNS计算结果Nu数随Ra数的变化,与三维的结果有些不同.当一定的Ra数范围由于羽流的不规则运动而向方腔中部聚拢,使得大尺度环流路径周长变小,传热Nu数随Ra数的变化会向下偏离GL理论预测倍数线.其原因可能是与GL理论的四壁边界层和中心区的假设不完全符合造成的.二维湍流热对流传热特性随Ra数变化出现两个转折现象的物理因素,还需要更多深入的研究.

综上所述,二维湍流热对流的计算结果在极高Ra数的情况下,传热Nu数随Ra数的变化依然可以与GL理论预测倍数线符合良好,即保持与GL理论预测相同的变化规律.这对进一步开展极高Ra数的二维湍流热对流DNS计算以及湍流热对流“终极态”的研究,有着重要的意义.

6 结 论

本文采用新的计算方法完成了系列二维高和极高Ra数的湍流热对流的DNS模拟.通过研究二维高和极高Ra数的传热变化特性,发现在一定的Ra数范围Nu数随Ra数的变化标度律减小,与GL理论预测倍数线偏离,但随着Ra数进一步提高偏离现象消失,而这一传热特性随Ra数的变化规律与大尺度环流路径的大小值相关.通过研究得到以下结论.

1)创建了高效二维热对流的并行直接求解方法(PDM-DNS),完成了系列不同Ra数的二维湍流热对流DNS模拟,包括Ra=1013的DNS模拟.

2)二维湍流热对流中,传热Nu/Ra0.3数随Ra数的变化与反映羽流运动的大尺度环流路径的变化具有很好的相关性,具有两个转折点.Pr=0.7时,Ra ≈ 109大尺度环流变为圆形,大尺度环流周长CLSC随Ra数变化突然减小,出现第一转折点,在Ra ≈ 1010时CLSC最小,出现第二转折点,而后随Ra增加变大.同样传热特性随Ra数的变化存在对应的两个转折.当Ra ＞ 109时Nu/Ra0.3随Ra变化的标度律减小,出现偏离GL理论预测倍数线的现象.在Ra数大于第二转折点Ra ≈1010,Nu/Ra0.3随Ra变化的标度律变为增加并再次与GL理论预测倍数线相符合,直到Ra=1013.Pr=4.3时的结果与上述变化规律一致,只是对应转折点的Ra数较高.

3)当Ra数大于第二转折点时,2D湍流热对流的传热Nu/Ra0.3随Ra数的变化与GL理论预测倍数线符合良好,表明极高Ra数情况下2D热对流的传热特性与GL理论预测变化趋势基本一致.