二维湍流热对流最大速度Re 数特性及流态突变特征Re 数*

何建超 方明卫 包芸†

1) (北京航空航天大学,航空科学与工程学院,北京 100191)

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

本文计算系列二维湍流热对流,Prandtl(Pr)数和Rayleigh(Ra)数范围分别为0.25—100 和1×107—1×1012,研究Reynolds(Re)数的变化规律.以最大速度计算的Re 数与Ra 数存在标度律关系,但中间出现间断.研究表明,大尺度环流形态由椭圆形到圆形的突变引起流动失稳,导致最大速度值间断下降,影响Re 数变化趋势的连续性.所有Pr 数对应的流态突变特征Re 数为常值,Rec 约为1.4×104,即当Re 数达到特征Rec 时,大尺度环流形态会发生从椭圆形到圆形的突变.间断点对应的Rac 与Pr 数之间存在标度关系Rac-Pr1.5.对Ra 数进行补偿平移,所有Pr 数的Re 与RaPr—1.5的变化曲线重合,不同Pr 数有相同的间断临界点位置,RacPr—1.5=109.

1 引言

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

在RB 热对流系统中,Rayleigh 数(Ra)、Prandtl 数(Pr)和宽高比(Γ)是控制系统的3 个重要无量纲参数,而Nusselt 数(Nu)和Reynold 数(Re)分别是反映系统传热效率以及湍流强度的无量纲参数.多年来,两个响应参数与控制参数之间的关系是RB 热对流的研究重点[1-3],而其中影响最广泛的是Grossmann 和Lohse 提出的GL 理论[4-8].GL 理论在较大范围内预测了Nu(Ra,Pr)和Re(Ra,Pr)的变化趋势,至今也得到了比较多数据的验证[8-10].在GL 理论中,Ra-Pr相图被分为多个区域[4,8],在不同区域中耗散率与Ra,Pr的关系有所差异.也即,在同一Pr数下,随着Ra数增大,耗散率的发展趋势会出现变化,其他相关的物理量的行为也会有所变化[11].早期关于RB 热对流的研究中,发现到达临界Ra数Rac(1707)之后,流体运动会从一种定常(time-independent)状态转换变为随时间变化的非定常状态[11-15],而这个转变的Ra数与Pr数存在依赖关系[12].“芝加哥对流实验”[16-18]发现随着Ra数的增大,系统的流态会从对流状态变为湍流状态,同时系统的Nu数和Ra数关系也发生了转变.其中,当RB 系统转变为湍流状态后,系统中的主要结构,包括羽流、大尺度环流等,会随着Ra增大而发生突变[19-20].当大尺度环流由椭圆形转变为圆形时,系统的Nu数也会偏离GL 理论预测曲线[10],但在更大的Ra后,Nu数的变化趋势再次与GL 理论预测的一致[19-20].

有数值模拟研究[21]表明不同Pr数下,Re与Pr的标度律会出现变化,而Re与Ra的标度律均为0.53,与之前[18,22]获得的0.46 不同.最近,在准二维的实验中Li 等[23]在Pr=11.7—145.7的范围内发现Re的标度律在0.53—0.60 之间,与Chen等[24]在Pr=4—7 之间发现的0.55 十分接近.在二维DNS 中,Xu 等[25]发现在Pr=0.025 时Re与Ra的标度律为0.50,而Werne 等[26]在Pr=7 时的标度律为0.54.这些研究表明,Re(Ra)的关系也会随Pr数变化而略有变化,这变化值得去进一步深入研究.

本文采用高效并行直接求解方法PDM-DNS[27],在“天河二号”超级计算机上进行了多组Pr数和Ra数的二维湍流热对流DNS 模拟,Pr数从0.25到100,Ra数范围为1×107—1×1012,跨度为5 个量级,总共133 个计算算例.本文对多组Pr数和Ra数的二维湍流热对流中反映湍流特性的Re数变化规律进行研究,并探讨热对流流态突变时对应的典型特征Ra数和Re数的特性及其变化规律.

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

在RB 热对流的数值模拟中,通常引入Oberbeck-Boussinesq(OB)近似,对方程进行简化.在OB 近似下,无量纲化的热对流方程为

其中,u为无量纲速度矢量,θ为无量纲温度,p为压力,k为单位垂向矢量,方向与重力方向相反.数值计算中边界条件为壁面速度均采用无滑移条件,温度为侧壁采用绝热条件,上下底板采用恒温条件,上底板恒温冷却θtop=−0.5,下底板恒温加热θbot=0.5.

大规模湍流热对流的DNS 计算由于计算量巨大必须通过并行计算进行,其中压力泊松方程的并行求解是实现高效并行计算的关键.利用高效并行直接求解压力泊松方程的PDD 算法[28],建立了热对流DNS的并行直接求解方法(parallel direct method of DNS,PDM-DNS),并展现了很好的并行效率[27].在数值求解过程中,采用投影法进行计算,时间和空间均是二阶精度;压力泊松方程采用快速傅立叶变换(FFT)进行解耦,然后直接求解三对角方程组,而动量方程和温度方程采用显式格式推进.

本文使用高效的PDM-DNS 方法,在“天河二号”超级计算机上进行了多组Pr数和Ra数的2D 湍流热对流DNS 模拟.图1 给出了计算算例的Ra-Pr相图.本文中研究的Pr数从0.25 到100,Ra数从1×107到1×1012,跨度为5 个量级,总共133 个计算算例.本文的二维湍流热对流计算数据丰富,相图中呈现的结果是目前二维RB 系统较完整的数据.

图1 二维热对流计算算例Ra-Pr 相图Fig.1.The Ra-Pr Phase diagram explored in the present study.

为了充分识别系统的流动,本文采用了上下边界加密的网格进行数值模拟,计算网格大小满足Kolmogorov 尺度和Batchelor 尺度,即ηK=并且,边界附近的网格满足Shishkina 等[29]提出的要求,即:

其中,a≈ 0.482,Nθ,BL和Nv,BL分别表示温度边界层和速度边界层内的最少网格数.时间步长Δt小于Kolmogorov 时间尺度τK的1/1000.本文所有算例均在流场充分发展后进行统计,统计时间均大于200 个无量纲时间,并且在不同统计时间段求出的Nu数误差约为1%[20].

在本文的算例中,有少数算例会出现大尺度环流翻转的情况,具体包括Pr=1.2,Ra=2×107,Pr=2.0,Ra=1×107—1×108,Pr=4.3,Ra=5×107—2×108,Pr=100,Ra=2×108,共9 个算例.当仅统计顺时针旋转或逆时针旋转的Re数和Nu数时,两者的结果是一致的.因此,可以将顺时针旋转的数据进行水平翻转,叠加到逆时针的旋转的数据上,得到时间更长的时间平均场.如果统计数据包括大尺度环流的翻转过程,并且不作方向统一的处理,那么将对数据统计造成一定的影响,比如Re数会减小约50%,Nu数波动约为3%[30].因此,关于翻转算例的统计,仅统计顺时针和逆时针的流动,并作方向统一处理,不考虑发生大尺度环流翻转的过程.

3 二维湍流热对流Re 数特性

Re数是反映热对流系统的湍流流动特性的特征参数.本文研究系列Pr数及Ra数下的二维RB 系统湍流Re数与Ra数和Pr数之间的变化规律及特性.

图2 给出了一个典型流态的平均速度场云图.图2 可见,该速度分布成圆环状,方腔四个角落中的速度都较小,速度分布在中心速度很小,随着半径的变化速度快速增大,到接近一半的位置速度达到最大,而后速度随半径的增大速度减小.

图2 平均速度分布及最大速度点Fig.2.The time-averaged velocity filed and the maximum velocity point.

经过计算全场的速度大小,找到绝对速度的最大值,即图2 中白点处的速度,作为计算湍流热对流Re数的特征速度,其无量纲的计算公式如下:

3.1 二维湍流热对流Re 数变化规律及其间断特征

首先,以Pr=0.7的系列Ra数热对流为典型计算结果,探讨热对流的Re数随Ra数的变化规律.

图3 给出了Pr=0.7的系列Ra数热对流的Re数随Ra数的变化情况,图中为双对数坐标.可以看到,Re数随Ra数的增加逐渐增大,表明系统的湍流强度随Ra数增加逐渐变强.特别地,Re数的变化在Ra=109处出现了间断平移,分成了较明显的两段.

图3 Re 数与Ra 数关系Fig.3.Re as a function of Ra.

对于Re数与Ra数的变化关系中出现的间断平移,是由二维湍流热对流的流态突变造成的,将在下节内容中详细讨论.

3.2 二维湍流热对流流态突变及最大速度变化特性

在低Ra数时,二维湍流热对流的流态是倾斜的椭圆大尺度环流加两个角涡的流态.随Ra数增高,角涡会发生脱离导致流动失稳,系统的流态突变成圆形的大尺度环流流态[19].

图4 给出了Pr=0.7 时流态从椭圆突变成圆形前后典型流态的平均速度场云图,图中明显可见两种完全不同的速度分布.由于不同Ra数时速度值得变化范围较大,为了清楚地反映速度分布情况,每个Ra数对应的速度云图色标范围是不一致的,均设置为各个算例的最小速度到最大速度.在Ra数较低的图4(a)和(b)中,速度分布形态基本一致,均为倾斜的椭圆形和两个角涡,图中白点为最大速度点,都出现在椭圆和角涡的相切处.大尺度环流与角涡相切的位置速度都会比较大,而接近另外两个角落的速度较小.随着Ra数增高,Ra=1×109时流态发生突变,大尺度环流形态突变为圆形,如图4(c)和(d)所示.突变后高速区域呈圆环状,整体速度相近,最大速度随机产生在圆环上,并且最大速度相较突变前略有减小,这一点在图4的色标范围上可以看出.

图4 平均速度场变化及流态突变(Pr=0.7) (a) Ra=2×108;(b) Ra=5×108;(c) Ra=1×109;(d) Ra=2×109Fig.4.The time-averaged velocity fieldsc and the sudden change of flow pattern (Pr=0.7): (a) Ra=2×108;(b) Ra=5×108;(c) Ra=1×109;(d) Ra=2×109.

计算Pr=0.7 时所有Ra数情况下的最大速度值,探讨最大速度值随Ra数的变化规律.

图5 给出了Pr=0.7 时不同Ra数的平均速度场中最大速度值Umax随Ra数的变化情况.可以看到,最大速度的变化明显分为两段.当Ra≤5×108,大尺度环流流态为椭圆,随着Ra数的提高平均速度场中的最大值Umax逐渐变大.在Ra=1×109时,大尺度环流流态突变为圆形,最大速度Umax间断式突然下降.之后随着Ra数增大,平均速度场中的最大值Umax再次逐渐变大.当Ra≥5×1010后,最大速度Umax出现波动,但速度值变化不大.

二维湍流热对流在较低Ra数时呈椭圆形大尺度环流形态,由于羽流基本沿椭圆路径运动,流动相对集中稳定,使平均场速度较大.随着Ra数的提高,速度增大,带动角涡脱落造成流态失稳混乱,导致流态从椭圆到圆形的突变[31].当流态失稳突变后,羽流运动呈现较为混乱的绕行,会使平均后的速度场整体减小,所以导致图5 中最大速度间断式减小.

图5 最大速度随Ra 数变化Fig.5.The maximum velocity as a function of Ra.

热对流在流态突变时平均场最大速度Umax的间断式减小,是造成图3 中Re数随Ra数的变化规律出现间断的原因.

4 不同Pr 数的Re 数规律及临界Re 数

本文计算了10 个不同的Pr数系列的算例,数据较多.为清晰展示结果,在讨论二维湍流热对流Re数特性的Pr数影响时,首先选取3 个典型的Pr数0.7,4.3 和10 进行分析.

图6 给出了3 个Pr数系列的Re数的变化特性,蓝色六边形、红色方形和绿色实心圆点分别为Pr=0.7,4.3 和10,黑色线表示 1.4×104.由图6 中可见,3 个Pr数系列的Re数都分别随Ra数的增大而增大,同时Re数随Ra数分布有两个阶段,阶段之间存在明显的间断.将间断处的Re数定义为流态突变特征Re数Rec,Ra数定义为流态突变特征Ra数Rac.有意思的是,不同Pr数时流态突变对应的特征Re数Rec的值基本一致(Rec≈1.4×104).也即,无论Pr数为多少,当Re数达到Rec时,平均场的大尺度环流形态一定会从椭圆形突变为圆形.

图6 Re 数与Ra的关系图Fig.6.Re as a function of Ra.

这一现象的原因与上文讨论的Pr=0.7的一样,是由于最大速度值在流态突变处出现骤减,导致Re数变化出现间断平移.另外,随着Pr增大,发生间断平移的位置会向右移,表明随Pr增大,Rac会逐步增大.

计算所有Pr数下二维湍流热对流的Re数,Pr数范围为0.25—100,讨论这个范围内不同Pr数情况下,热对流的Re数与Rac的变化规律.

图7 给出相应的Pr从0.25 到100的系列Ra数下Re数分布.图中发现不同Pr数下湍流Re数分布特征与之前3 个典型Pr数情况相似.Re数从一百左右到几十万,跨度4 个量级.除了Pr=50 和100 时Re数数值不够大,其他的结果都在流态突变特征Re数处存在间断,间断点对应的Ra数依次往后移动.然而间断处的Re数,也就是流态突变特征Re数的值基本相同,是一个常数,Rec≈1.4×104.也即,对于不同系列Pr数和Ra数,当由最大速度定义的Re数达到特征Rec时,流态就会发生由椭圆型大尺度环流到圆形大尺度环流的流态突变,流体突变特征Re数Rec与Pr数和Ra数变化无关.

图7 不同Pr 数下Re 数与Ra 数的关系Fig.7.Re as a function of Ra at different Pr.

在不同Pr数下,Rac的位置逐渐右移,表明Rac随Pr增大而增大.在二维RB 热对流中,流态突变Ra数与Pr数有Ra-Pr1.5的标度律关系[30].为了探讨不同Pr数下Re数的转变共同特征,对横坐标Ra数进行补偿平移,研究Re与RaPr—1.5的关系.

图8 给出Re数随RaPr—1.5的分布.图中可见,10 个Pr数的Re(RaPr—1.5)分布完全重合在一起,表明在不同Pr数下,Re数的变化规律基本一致,呈现很好的自似性.

图8 不同Pr 数下Re 数与RaPr—1.5的关系Fig.8.Re as a function of RaPr—1.5 at different Pr.

图8 中黑色虚线表示Rec≈1.4×104,红色虚线表示RacPr−1.5=109.两条虚线的交点对应不同Pr数下的流态突变.这表明不同Pr数都有基本相同的流态突变特征点位置,RacPr−1.5=109,同时Re=Rec≈1.4×104.这是两个常数,可以用于预测不同Pr数和Ra数下流态发生突变的时机以及系统Re数的变化趋势,对区分不同参数下的流态以及流态突变特性的探讨有重要意义.当RaPr−1.5＜109时,系统的流态为椭圆形;当RaPr−1.5=109时,系统的流态会发生突变,系统的Re数达到Rec≈1.4×104;随后RaPr−1.5＞109,系统的流态为圆形.

另外,对图8 中所有数据进行Re-(RaPr—1.5)γ的拟合,得出: 流态突变前,Re数的标度律为Re-Ra0.61Pr—0.92,流态突变后为Re-Ra0.55Pr—0.83.

流态突变的临界Re数为常数,反映出二维湍流热对流特殊的流动形态变化特征,也将为二维湍流热对流的流动特性研究以及对应的传热特性变化特征等问题的研究提供新的思路.更多的价值还需要进一步深入的研究.

5 结论

本文采用高效并行直接求解计算方法PDMDNS 完成了系列Pr数和Ra数的二维湍流热对流的DNS 模拟,Pr数从0.25 到100,Ra数从1×107到1×1012,跨度为5 个量级,总共133 个计算算例,所呈现的结果是目前二维湍流热对流系统相当完整的数据.本文研究以平均场最大速度为特征速度的Re数特性以及最大速度的变化规律,发现大尺度环流形态由椭圆形变为圆形的突变对Re数特性的影响.研究结论如下:

1)典型算例Pr=0.7 时的Re数特性结果表明,Re随Ra数的变化存在标度律关系,但中间出现明显的间断现象.从基本流态特征出发进行研究发现,大尺度环流形态由椭圆形突变为圆形会引起最大速度值的间断式突然下降,导致Re数特性的间断现象出现.

2)不同Pr数的Re数特性研究表明,Re数随Ra数的变化特性有明显的自相似性.并且发现,流态突变对应的特征Re数Rec为常数,与Pr数和Ra数变化无关,Rec≈1.4×104.当Re数达到特征Rec时,大尺度环流形态会发生突变,从椭圆形变为圆形.这为RB 热对流的流态区分提供新的方法,并且进一步根据Re数的自相似性进一步预测Re数的变化规律.

3)流态突变特征Ra数Rac与Pr数之间存在标度关系Rac-Pr1.5.对横坐标Ra数进行补偿平移,Re与RaPr—1.5的变化曲线完全重合,表明不同Pr数时Re数随Ra数的变化规律基本一致,流态突变前后分别有Re-Ra0.61Pr—0.92和Re-Ra0.55Pr—0.83的标度律关系.不同Pr数有基本相同的间断特征点位置,RacPr−1.5=109.