旋转肥皂泡热对流能量耗散与边界层特性的数值模拟*

贺啸秋 熊永亮† 彭泽瑞 徐顺

1) (华中科技大学航空航天学院,武汉 430074)

2) (工程结构分析与安全评定湖北省重点实验室,武汉 430074)

将底部加热的半个肥皂泡作为一个新的热对流模型,结合了肥皂泡固有的球面与准二维特征,由此有助于理解行星大气流动中的复杂物理机制与热对流特性.本文使用直接数值模拟方法计算了旋转肥皂泡上的湍流热对流,研究了肥皂泡上的温度与黏性边界层以及拟热能和动能耗散规律.结合肥皂泡上温度场与速度场特征,分别根据温度脉动均方根最大值以及速度脉动边界处斜率延长线与最大值交点提出了肥皂泡上温度与黏性边界层的识别方法.研究发现,当肥皂泡从边界吸收能量时,拟热能耗散与动能耗散均集中在边界层中,肥皂泡上的温度边界层与黏性边界层厚度与瑞利数 Ra 存在明确的标度关系.相比经典Rayleigh-Bénrad对流(RB 对流)模型,温度标度指数具有较为接近的结果,但速度标度指数存在一定的差异.此外,在混合区,均方根温度(T*)随纬度(θ)具有近似 T*～θ0.5 的标度关系,这与RB 对流模型及其相应的理论预测一致.最后通过能量平衡方程发现,肥皂泡上拟热能内耗散率 和动能内耗散率 比拟热能外耗散率和动能外耗散率 大1 个量级,拟热能与动能的内部耗散率在边界层中具有支配地位,随着肥皂泡旋转速率的增加,热羽流难以输运到高纬度地区,进一步降低了拟热能与动能外耗散率的影响.

1 引言

海洋的平均深度约为数千米,但纵向长度可以达到数千公里;同样,地球大气圈的厚度只有数十公里,但是在经度与纬度方向的尺寸达到了数千公里,因而大尺度的海洋和大气流动具有一定的二维流动特征[1,2].先前的研究表明,二维湍流可以较好地近似描述宏观大尺度的大气与海洋流动规律[2].类似地,肥皂薄膜的宏观尺寸可以达到厘米量级,但是其厚度只在微米量级,因此流体在肥皂薄膜中的湍流运动非常好地近似了二维湍流,与行星大气与海洋的大尺度流动具有一定相似性[3,4],实验中广泛地将肥皂薄膜应用于二维湍流的研究[5].Wu等[6]使用肥皂薄膜研究了二维Couette 流动,并验证了被动标量在二维湍流中的扩散所遵守的标度率.Kellay等[7]则使用竖直肥皂薄膜验证了二维湍流的能量逆级串规律.近十年来,曲面肥皂薄膜[8]也被应用于实验研究中: 底部加热的肥皂泡作为一个较独特的热对流模型受到了很多学者的关注[3,4,9,10-14].Kellay[8]成功将半个肥皂泡长时间稳定地固定在一个水浴加热的铜制基座上.通过旋转铜质基座,可以实现肥皂泡绕对称轴进行旋转[13].肥皂泡被底座加热,肥皂水在浮力的驱动下进行对流输运;当加热足够强时,肥皂泡上就产生了湍流热对流[9].

在实验中,肥皂泡上出现了能长时间存在的岛涡,学者发现岛涡的运动轨迹与飓风轨迹均满足超扩散关系[9].实验结果表明肥皂泡上岛涡的运动轨迹特征以及强度与随机观察的一百多个飓风的轨迹以及强度之间具有非常类似的特征[9,12],学者依据这些特征建立了快速预报飓风轨迹锥的模型[3,12].肥皂泡模型刻画了流体在二维曲面中的湍流热对流运动,同时也有助于深入了解二维湍流热对流的物理机制[4,14].现阶段,对肥皂泡大尺度特征[3,9,12]和小尺度特征[4,10,13,14]都开展了一定研究.Bruneau等[4]使用直接数值模拟(DNS)方法研究了静止肥皂泡上的湍流热对流,发现了肥皂泡上能量串级符合Bolgiano-Oberbeck 的标度规律.而He等[14]使用DNS 方法研究了旋转对肥皂泡上流场与温度场的影响,获得了旋转对肥皂泡上小尺度脉动的影响规律.

在三维湍流系统中,能量由大尺度注入湍流,并在惯性区间内发生能量级串,最终在Kolmogrov尺度以下被消耗殆尽[15,16].而在二维情况下,湍流中存在双向能量级串: 湍流动能与拟涡能沿不同的方向传递[2].能量耗散对湍流场中的当地脉动具有重要的影响.惯性区间内的动能级串的过程中,二阶速度结构函数的标度直接受到动能耗散率的影响[17,18],其中为速度场,为运动黏度.在湍流热对流中,温度是主动标量,温度脉动在流体中的传播也具有级串形式[19],其标度规律与拟热能耗散率直接相关[20],其中为热扩散系数,而为温度场.同时,能量耗散的高低对湍流热对流系统的全局热通量具有调节作用[21].于是,研究能量耗散对揭示湍流热对流的物理机制具有重要的作用.例如,经典的Rayleigh-Bénrad 热对流(RB 对流)包含一层厚度为的流体和上下两个水平的恒温边界: 下边界温度高于上边界温度,温差为,热能被流体从下边界输运至上边界[22,23].RB 对流的控制参数为瑞利数Ra和普朗特数Pr:

其中是重力加速度的大小,是热膨胀系数.人们最关心Ra与Pr如何影响RB 对流的热输运效率[22,23].热输运效率用努塞尔特数Nu表征,即通过这一层流体的无量纲热通量:

式 中是选取的特征速 度.Nu与和有以下准确关系[21]:

式中,〈·〉V代表对流体体积求空间与时间平均.GL理论基于(5)式和(6)式推导出了Nu与Re随Ra和Pr变化的公式,并得到了大量研究的验证[24-28].GL 理论认为,在湍流热对流中能量耗散通常集中在边界层.事实上,Grossmann 与Lohse[24-26]在描述RB 对流中的能量耗散时,将边界层中的作用与湍流体中的作用分开讨论.人们在进一步发展GL 理论的过程中,将羽流的影响加入了这一物理图像,并将羽流归纳为边界层流体脱离边界而产生的运动[27,28].许多数值仿真[29,30]与实验研究[31-33]都验证了这一物理图像的正确性.最近,人们又利用区分研究边界层与湍流体中的平均能量耗散这一思想将GL 理论推广至水平热对流[34]、双标量RB 对流[35]、内加热对流[36]以及倾斜RB 对流[37].对于以上提及的多个热对流模型,GL 理论的推广理论成功得到了Nu和Re随Ra和Pr变化的规律,并被多个相互独立的实验和数值研究所验证[34-37].这说明了GL 理论具有较好的通用性,GL 理论所描述的物理图像在湍流热对流中广泛存在.然而,据作者所知,肥皂泡模型的研究中尚未有专门和系统的工作投入于研究能量耗散与边界层特性.所以本文希望借助GL 理论提出的物理图像,对肥皂泡上热对流的能量耗散与边界层特性进行深入的研究,并分析肥皂泡旋转所带来的影响,从而对肥皂泡模型的物理特点做进一步的揭示.

由于RB 对流简单的边界与几何特征中蕴含了热对流的本质规律,目前对热对流的研究主要集中在经典的RB 对流模型,人们对RB 对流中能量的输运与耗散做了深入细致的研究,取得了丰富的成果[22,23].在经典RB 对流模型的研究中,旋转所带来的影响也是热点方向之一,即所谓的旋转Rayleigh-Bénrad 对流(RRB 对流).旋转作用由罗斯比数的倒数 1/Ro表征,1/Ro越小,旋转作用越弱;反之 1/Ro越大,旋转作用越强.对于旋转速度为的RRB 对流,罗斯比数Ro的定义为

RRB 对流中,随着 1/Ro的逐渐增加,流场依次进入弱旋转态、中旋转态与强旋转态[38].在弱旋转态中,流场中旋转的作用基本可以忽略[39].而当流场处于中旋转态时,由于Ekman 泵吸效应,Nu随着1/Ro增大而有一定程度的上升[40-43].在强旋转态中,流体的运动被抑制,Nu随着 1/Ro的增大而剧烈下降[44].在靠近上下水平边界处,存在一层非常薄的温度边界层,RB 对流的主要热阻源于温度边界层[22].而RB 对流中能量的输运与耗散均与边界层有重要的关系[29,30,33].在RB 对流中,大尺度环流(largescale circulation,LSC)和角涡是两个独特的流动特征,对Nu与Re有着重要的影响.黄茂静和包芸[45]利用DNS 方法研究了二维与三维方腔在不同Ra下温度边界层的厚度变化.她们指出由于二维与三维情况下角涡数量的不同,导致了二者在温度分布变化上的差异,但二维与三维流动中边界层厚度都与Ra存在明确的标度关系[45].在二维较高Ra(Ra=1013)下,包芸等[46]发现Nu和LCS 受Ra的影响.他们发现,两者随Ra的变化规律有较强的关联性,而GL 理论可以部分准确地预测Nu随较高Ra的变化规律[46].另外,包芸等[47]研究了Ra较高的条件下(Ra=1010),温度边界层厚度和Nu随Pr变化的规律.其中温度边界层的厚度随Pr的变化并不明显[47].在较低Pr数下,Nu随Pr上升而增加;而当Pr数较高时,难以观察到Nu随Pr的变化[47].Wang等[48]通过实验测量了准二维系统中边界层的温度脉动的空间分布规律,并发现平均温度分布与边界层厚度存在一个普适的标度关系,并指出边界层方程对理解边界层脉动提供了一个理论框架.另外,Zhou 和Xia[49]对下边界的温度边界层进行了高精度测量,他们发现脉动温度均方根(RMS)的最大值距离下边界的距离,是一个可以非常好的描述温度脉动统计特征的空间尺度.

LSC 与黏性边界层有着密切关系,Sun等[50]利用粒子图像测速(particle image velocimetry,PIV)方法实验测量了在大尺度环流作用下的边界层特征,发现黏性边界层厚度随Re的标度指数是-0.5.Zhou等[51]也分析了温度与黏性边界层剖面曲线,他们发现在使用适当的物理量进行缩放后,RB 对流的温度与黏性边界层剖面与经典Prandtl-Blasius(PB)层流边界层符合得较好.最近,方明卫等[52]对温度边界层剖面曲线与PBP (Prandtl-Blasius-Pohlhausen)解进行了拟合,找到了拟合参数随Ra与Pr的变化规律.而在边界层外的混合区中,脉动温度的均方根(root mean square,RMS)剖面曲线与到边界的距离满足幂律函数的形式[50].Sun等[50]测量的脉动温度RMS 曲线与Adrian 的理论研究[53]所预测的结果一致.而最近的实验研究[54]也获得了满足幂律函数形式的脉动温度RMS 分布曲线.

对于RB 对流的能量耗散,学术界也取得了很多重要研究进展.在实验研究方面,He等[31,32]首先成功测量出RB 对流中空间不同位置的时间平均拟热能耗散率代表求时间平均.他们发现,在温度边界层中受平均温度的主导,而在边界层外的湍流体中受脉动温度的主导[31,32].Ni等[55]使用拉格朗日粒子追踪法直接测量了二阶速度结构函数,并利用结构函数间接获得了当地动能耗散率代表对空间某处局部体积 δV求空间与时间平均.他们发现,随Ra的变化可以用幂律函数进行标度[55].数值研究方面,Zhang等[29]研究了水(Pr=5.3)与空气(Pr=0.7)中,动能耗散率和拟热能耗散率的统计特性、空间分布规律以及与Ra的标度关系.而Xu等[30]研究了极低Pr的条件下(Pr=0.025)拟热能耗散率在湍流体与边界层中的统计特征.他们发现尽管边界层区域较小,但其对总的拟热能与动能耗散占到支配作用,并建立了和与Ra的标度关系[29,30].在理论研究方面,Petschel等[56]从能量耗散的角度定义了耗散边界层,并推广到了任意流动.

本文借助RB 对流中所建立起的经典研究方法与物理概念,利用DNS 可以直接获得肥皂泡上完整的速度场和温度场信息的优势[57],结合肥皂泡模型的特点,提出了肥皂泡温度与黏性边界层厚度的定义方法: 温度边界层的厚度为脉动温度RMS 最大值所处纬度距离赤道的大地距离;而黏性边界层厚度为纬度方向脉动速度RMS 曲线近赤道线性部分的延长线与其最大值的交叉纬度距离赤道的大地距离.同时本文也对肥皂泡模型中的能量耗散与边界层特征进行分析与研究,比较和RB 对流模型的异同.相比RB 对流,肥皂泡热对流模型中没有明显的LSC 与角涡的存在,且浮力强弱随着纬度而发生变化,其球面特征还使得不同纬度上空间的热通量强度不同.通过验证RB 对流中物理规律是否可以用于描述肥皂泡模型的特征,既可以促进充分认识肥皂泡模型的物理机制,同时也检验RB 对流模型中理论的通用性.本文的结构如下: 第2 节详细介绍DNS 所使用的数学模型;第3 节选取了多个不同数量级的Ra和罗斯比数开展数值研究,探索了肥皂泡上的能量输运和流场特征,Nu与Re随Ra以及 1/Ro的变化规律,同时还讨论了肥皂泡上边界层和能量耗散特征;第4 节进行总结分析.

2 DNS 数值模拟方法

2.1 数值模型

其中u和T分别为量纲一的速度场与温度场.肥皂泡控制方程的量纲一参数为Ra,Pr和Ro,它们的定义分别为

1/Ro的大小与转速成反比: 1/Ro越大,旋转的作用就越大.

与经典的Rayleigh-Bénard 方程相比,肥皂泡上流场的控制方程多了两项: 外冷却项-ST和外摩擦项-Fu.S和F分别为量纲一冷却系数与摩擦系数,外冷却项与外摩擦项表征了肥皂泡模型与外部环境中冷空气产生热交换与摩擦.当S=0且F ̸=0时,肥皂泡模型的能量方程和Oberbeck-Boussinesq 方程组中的能量方程一样.由于肥皂泡只有赤道一条恒温热边界,能量单向通过赤道流入肥皂泡,使肥皂泡上的温度随着时间的增加而不断增加,这一过程为非平衡过程.最终肥皂泡上的温度处处与赤道温度相同时,肥皂泡停止从赤道吸收能量而达到平衡态,此时流体失去浮力驱动形成热寂[4,14].冷却项使肥皂泡上的平均温度始终小于赤道温度,在适当的S取值范围下,可使极地平均温度在数量级上远低于赤道温度,近似为0.肥皂泡的赤道与顶部的温差近似保持恒定,使热对流在恒定的Ra下可以一直进行下去,并使肥皂泡进入统计稳平衡态[4,14].而如果S ̸=0 且F=0 时,肥皂泡的DNS 随着仿真时间步的推进而无法避免地产生爆炸[11,4,14].肥皂泡上的湍流热对流属于二维湍流,存在能量逆级串的现象[11,4,14].能量逆级串使流体从赤道吸收的热能以羽流相互作用的方式形成大尺度的经向流[2,11],最终不断快速汇聚的经向动能使得CFL 条件永远无法得到满足从而使DNS 仿真失败[11].向控制方程中添加外摩擦项,其意义在于抑制肥皂泡上大尺度流动结构的能量不受限制的积累,使肥皂泡进入统计平衡态.在以往的二维湍流DNS 的研究中,添加外摩擦项也是一种普遍采用的方法[2].以往的研究已经表明,S=F=0.06 及其所在的量级可以较好地模拟肥皂泡上的统计平衡热对流[4,14],此时北极部分的平均量纲一温度相比赤道量纲一温度T0=1 有多个数量级的差异,可近似为0.同时本文也会对两者带来的外耗散影响进行讨论,以验证外冷却项与外摩擦项对流动的作用不会影响系统的物理机制.

以肥皂泡的球心为原点建立量纲一三维笛卡尔坐标系 (x1,x2,x3),x1和x2坐标轴均在赤道平面.三个方向的基矢分别为e1,e2和e3.e3的方向为从原点指向肥皂泡的球顶.肥皂泡模型的边界为赤道,赤道上保持恒温无滑移边界条件(u=0与T=1).肥皂泡的曲面方程为

而描述赤道的曲线方程为

量纲一三维笛卡尔坐标系 (x1,x2,x3) 所对应的球坐标系为(r,ϕ,ψ) :r为径向长度,ϕ为方位角,ψ为极角.球坐标系三个方向的单位基矢分别为er,eϕ和eψ.为了方便描述数值计算结果,引入地理坐标系中的经度与纬度来描述肥皂泡上任意一点的空间位置: 经度坐标与方位角定义相同,也用ϕ表示;纬度坐标是极角ψ的余角,用θ=π/2-ψ表示.在赤道上纬度θ为0°,而极角ψ=90°;在肥皂泡的顶点纬度θ为9 0°,而极角ψ=0°.经度方向基矢与方位角方向基矢相同,为eϕ.纬度方向基矢与极角方向基矢方向相反,eθ=-eψ.

2.2 数值方法

控制方程(8)是向量形式,在求解时需要根据坐标系确定算符的具体形式.由于肥皂的几何外形为二维半球曲面,使用笛卡尔坐标系求解控制方程(8)需要额外的方法处理曲面,增加了求解难度.如果使用球坐标系求解,众所周知,球坐标系的NS 方程在极点处引入了奇异性,因而必须通过其他方式回避这一问题.球极投影坐标系 (x,y) 可以解决肥皂泡曲面几何所带来的难题.球极投影的几何关系如图1 所示.肥皂泡球面上任意A点与南极点B之间的连线在赤道平面上唯一相交于点C,通过该一一映射关系可以将球面上所有点映射在一个平面圆上.而球极投影坐标系的原点x=y=0也处于肥皂泡的球心.x方向与y方向分别与笛卡尔坐标系中x1方向和x2方向重合.

图1 不同坐标系下的计算域与球极投影平面投影示意图Fig.1.Computational domain under different coordinate system and the illustration of the stereographical projection.

通过简单的几何关系推导,可以得到笛卡尔坐标系 (x1,x2,x3) 与球极投影坐标系 (x,y) 的变换为

而逆变换为

为获得球极投影坐标系中的矢量算符的具体形式,首先需要计算球极投影坐标系的拉梅系数.x和y两个方向的拉梅系数可写为

由于两个方向的拉梅系数相同,下文定义L=Lx=Ly.球极投影坐标系在x与y方向上的单位基矢分别是hx与hy,分别为

速度矢量在球极投影坐标系中可以写为

其中u和v分别为速度矢量在x和y方向的投影.速度矢量在笛卡尔坐标系、球坐标系和地理坐标系中的各个分量为

由于肥皂泡为二维半球面,径向速度恒为0:ur ≡0.而根据基矢关系可以知道,uθ=-uψ.

根据曲线坐标系中梯度算符的通式,可以写出∇算符在球极投影坐标系中的形式:

而拉普拉斯算符 Δ 在球极投影坐标系中为

拉格朗日导数为

在获得以上单位基矢与算符的转换关系以后,可以将球极投影坐标系内的算符代入到矢量形式的控制方程中,从而获得球极投影坐标系下的控制方程.连续性方程可以写为

动量方程为

能量方程为

在球极投影坐标系 (x,y) 中,肥皂泡的曲面方程为x2+y2≤1,是一个半径为1 的圆.将计算域设定为边长为1.01 的正方形,使用交错结构化网格来离散计算域: 温度与压强结点在网格的中心,而速度结点位于网格边的中点.网格结点分布如图2 所示,其中i和j分别为x和y方向上的结点编号.

图2 网格结点分布图Fig.2.Node distribution inside a grid cell.

边界条件使用惩罚法来实现,时间离散格式为二阶Gear 格式,而非线性项使用三阶Murmanlike 格式来离散,更详细的数值方法与验证可参考文献[4,14].

2.3 算例设置

表1 列出了本文肥皂泡模型DNS 算例的全部参数信息.由于肥皂水中肥皂的浓度很低,其物性与水非常接近,算例的Pr固定为7 不变.算例的Ra的变化范围为3×106—3×109,小Ra下采用1024×1024 的数值解析度,在最高Ra下网格分辨率达到了 2 048×2048,以往的研究表明这一数值解析度可以充分获得网格独立解[4,11,14].

表1 算例详细信息Table 1.The detailed information of the cases.

2.4 统计稳态

在数值实验中,肥皂泡的初始状态为静止态,即u=0 且T=0.随着计算的进行,时间步不断向前迭代,肥皂泡上的流体开始受热并运动,产生湍流热对流逐渐达到统计稳态[4,14].此时肥皂泡上的湍流得到了充分发展,肥皂泡吸收的能量与其耗散的能量达到了动态平衡.在这样的条件下,流场量可以分解为平均值与脉动值之和,如温度T:

在下文,符号〈·〉表示时间和经度空间平均.肥皂泡在经度ϕ方向具有旋转对称性,流场变量在统计上并不依赖经度坐标ϕ.于是,在控制参数确定的情况下,所有流场变量的平均值是纬度θ的函数.

3 结果与讨论

3.1 瞬态流场

对所研究的不 同Ra都考虑了4 个 1/Ro: 0,0.1,1 和 10,以便分析旋转作用的影响.这里首先从最直观的流场瞬时云图来展示不同Ra与不同旋转角速度 1/Ro给流场带来的影响.图3 给出了不同Ra与 1/Ro条件下温度场T在肥皂泡上的瞬时分布云图.当Ra=3×106时,在没有旋转的情况下(1/Ro=0),可以从瞬时温度场中清晰地观察到羽流的形状特点: 羽流尺寸大,而数量相对少.随着Ra升高至Ra=3×109,肥皂泡表面的羽流的尺寸大幅减小,且更加卷曲,同时羽流数量大幅增加.在肥皂泡高速旋转的情况下(1/Ro=10),羽流的形状与尺度大小同未旋转的情况下相比并无明显变化,但是羽流造成的影响难以上升到肥皂泡顶部区域.

图3 不同 1/Ro 与 Ra 下,T 在肥皂泡上的分布 (a) 算例 Ra=3×106,1/Ro=0;(b) 算例 Ra=3×109,1/Ro=0;(c) 算例 Ra=3×109,1/Ro=10Fig.3.The instantanous temperature field with different 1/Ro and Ra : (a) The case of Ra=3×106,1/Ro=0;(b) the case of Ra=3×109,1/Ro=0;(c) the case of Ra=3×109,1/Ro=10.

同时,图4 给出了在不同Ra与 1/Ro设置下,肥皂泡上的瞬时动能的对数 l gEk=lg(1/2u2).在肥皂泡没有旋转的情况下,在Ra=3×106时,可通过区分动能集中的区域识别出流动结构的尺寸较大;而当Ra=3×109时,流动结构的尺寸明显变小,同时变得更加破碎.如果保持Ra大小不变,在高速旋转时 1/Ro=10,肥皂泡上靠近赤道的部分出现了大尺寸且动能集中的流动结构,这是旋转带来的稳定经向流.另一方面,在肥皂泡的顶部附近,流动结构的动能较小,但尺寸相比没有旋转的情况更大.

图4 不 同 1/Ro 与 Ra 下,log Ek 在肥皂泡上的分布 (a)算例 Ra=3×106,1/Ro=0;(b)算例 Ra=3×109,1/Ro=0;(c)算例 Ra=3×109,1/Ro=10Fig.4.The instantanous filed of log Ek with different 1/Ro and Ra : (a) The case of Ra=3×106,1/Ro=0;(b) the case of Ra=3×109,1/Ro=0;(c) the case of Ra=3×109,1/Ro=10.

3.2 Nu 与Re

在经典RB 热对流模型中,热能从下边界注入,从上边界流出,Nu被用来表征流体对流输运热能的效率.但在肥皂泡模型中,只有赤道边界:热量通过这条边界流入肥皂泡,然后通过外部和内部耗散在肥皂泡上消耗殆尽.本文使用Nu来表征肥皂泡赤道上的对流传热效率,Nu同时也代表了肥皂泡从外界吸收能量的效率,Nu越大意味着肥皂泡通过对流从外界吸收热能的能力越强.于是,Nu的定义为赤道处的无量纲热通量,其形式如下:

在热对流中,温度较高的流体受到浮力的作用从赤道向肥皂泡的顶部运动,将热能与动能输运至肥皂泡内部.但是,可以证明在经度方向上的平均速度为0.首先写出球坐标系下的平均连续性方程:

图6 给出了Nu与Re随Ra变化的规律.从图6可以看出,Nu随Ra变化近似满足标度率～Ra13,这与许多RB 对流中所观察到的结果非常接近.图6 中结果表明,1/Ro的变化对Nu的影响非常小,只有在Ra较低(Ra≤3×107),并且1/Ro=

图6 Nu(上)和 Re(下)随 Ra(左)和 1/Ro(右)的变化规律 (a) 在 1/Ro 恒定的条件下,Nu 随 Ra 的标度规律;(b) 在 Ra 恒定的条件下,Nu 随 1/Ro 的变 化规 律;(c) 在 1/Ro 恒定的条件下,Re 随 Ra 的标度规 律;(d) 在 Ra 恒定的 条件 下,Re 随 Ra 的变化规律Fig.6.The variation of Nu(up) and Re(down) with Ra(left) and 1/Ro(right): (a) The scaling behavior of Nu with Ra in condition of fixed 1/Ro;(b) the variation of Nu with 1/Ro in condition of fixed Ra;(c) the scaling behavior of Re with Ra in condition of fixed 1/Ro;(d) the variation of Re with 1/Ro in condition of fixed Ra.

10 的条件下,Nu有很微弱的上升.综合全部结果可以发现,1/Ro对Nu标度率的影响可以忽略.从定义公式(33)可以知道,Nu的大小主要取决于边界层中的平均温度分布曲线.图5 数据表明,越靠近赤道,平均温度随着旋转速度的变化越微弱[14].这暗示科里奥利力的作用随着纬度的减小而减小.如果将速度u=uϕeϕ+uθeθ代入动量方程的科里奥利力项,则可以得到:

图5 不同 1/Ro 和 Ra 下,平均温度 〈T〉 随纬度 θ 的变化规律 (a) 算例 Ra=3×106;(b) 算例 Ra=3×107;(c) 算例Ra=3×108;(d) 算例Ra=3×109Fig.5.The variation of mean temperature 〈T〉 with the latitude θ for the different 1/Ro and Ra : (a) The cases of Ra=3×106;(b)the cases of Ra=3×107;(c) the cases of Ra=3×108;(d) the cases of Ra=3×109.

当纬度为0 时,eθ=e3所以e3×eθ=0;另一方面,在赤道上有e3×eϕ=er;所以赤道上科里奥利力项可以写为.由于肥皂泡二维半球面的几何特征,径向上的力对流体的运动没有任何作用.在赤道上,科里奥利项的作用为0,而且在靠近赤道的区域内,科里奥利力的作用都很弱.这就揭示了旋转对Nu的影响非常微小的原因.肥皂泡上Nu随 1/Ro的变化规律与旋转RB 对流有较明显的不同.旋转RB 对流中Nu随着 1/Ro的逐渐增加呈现先升后降的变化规律: 在合理的参数范围内,1/Ro可以使Nu最大增加30%[41,40].在适当的旋转速度下,旋转RB 对流中的羽流变为竖直排列的涡,产生Ekman 抽吸效应,将边界层中的流体向湍流体中输运,使Nu有明显的上升[43,42].Ekman抽吸效应的效率随着Pr数的上升而得到加强,而随着Ra数的上升Nu被削弱[40].但是只有在三维流场中才会产生Ekman 抽吸效应,所以肥皂泡上Nu上升的机理与RB 对流具有较大不同.

另一方面,Re随着Ra的变化曲线也具有幂律函数的形式.当肥皂泡静止或者旋转速度较小时(1/Ro≤1),Re满足标度率Re～Ra0.48,而 且Re随1/Ro的变化可以忽略.但当 1/Ro=10 时,在较低Ra的条件下,Re出现了明显的下降,而在较高Ra的条件下,Re随 1/Ro的变化依然较小,这导致Re的标度系数从 0.48 上升到 0.51.Re同时受到Ra与1/Ro的影响,与随纬度分布的规律直接相关[14].在之后讨论黏性边界层的小节会展示,当1/Ro=10时,纬度方向速度的脉动受到极大的抑制,max()也 随之下降.但 随着Ra上升,m ax() 下降幅度逐渐缩小,使Re的降幅也变小.

3.3 均方根温度分布与温度边界层

RB 对流的温度边界层厚度δT可以根据平均温度〈T〉来定义,也可以根据脉动温度的均方根来定义[49].RB对流的实验发现,当非常靠近上下平板时,流体的平均温度〈T〉与到边界的距离d之间的变化曲线表现出较好的线性规律.对线性分布区域做拟合线并且延长,当拟合线达到湍流体区域的平均温度时,其与边界的距离可定义为边界层的厚度[50,49].但是,肥皂泡对流模型中只有一个热边界,所以在边界的法线方向上流场并不对称,且静止肥皂泡上的经向流非常微弱[14],这导致流场中不存在一个有明确物理定义的湍流体区域.RB 对流还可以用另外一种方式定义δT: 脉动温度的均方根达到最大值处到边界的距离[29,50,49].因此,本文根据肥皂泡的曲面外形特点,采取以下方式定义温度边界层厚度: 均方根脉动温度T*是θ的函数,当T*达到最大值时,θ=θT.纬度θT距离赤道的大地距离(即半球面上由赤道到纬度θT的最短距离: 圆周角为θT的大圆弧长)定义为温度边界层厚度δT.在无量纲坐标系下,肥皂泡的半径R ≡1,所以δT和θT之间的关系为θT=δT/R=δT.在肥皂泡上,纬度θT的物理含义与δT相同,代表了流体到边界(赤道)的距离,两者的差别在于量纲.

图7 的左边给出了Ra=3×106,Ra=3×107,Ra=3×108和Ra=3×109的静止肥皂泡上T*随纬度的分布曲线,并且用虚线标注了θT,即边界层的厚度.图7 的右边则给出不同罗斯比数下边界层厚度δT与Ra的关系.图7 中的数据显示,δT随Ra变化呈指数标度,而 1/Ro对温度边界层厚度的影响并不明显.δT与Ra之间的标度关系为δT～Ra-0.32.和3.2 节讨论的原因一样,动量方程中科里奥利的作用项随着纬度的降低而减小,在赤道处变为0.而边界层厚度很薄,而且随着Ra的上升而快速下降,所以边界层厚度受旋转作用的影响非常微弱.值得一提的是,在RB 对流模型的研究文献中,温度边界层厚度的标度规律δT～已经多次被研究者发现[45,49].在肥皂泡模型中,δT与Ra的标度关系与RB 对流非常相似.并且,肥皂泡的δT标度系数与Nu的标度系数完全相反,这也符合RB 对流中的特征.从几何角度观察,温度边界层靠近赤道且厚度很薄,流动尺度很小.之前的研究也表明,在边界层中受无滑移边界条件的影响,流体的运动速度较低且温度的脉动也很微弱[4,14],于是人工散热与人工摩擦项的作用在边界层中不明显.同时,边界层中浮力的投影与重力的夹角可以近似为1 80°,与RB 对流中的情况相同,这就说明了肥皂泡在边界层中热能输运的物理机制与RB 对流很相近,使得δT随着Ra的变化规律也与RB 对流接近.同时,由于温度边界层之外,温度梯度较低;因而对于这一系统,温度边界层的厚度也恰好与系统边界层的温度梯度成反比关系,因而其与Nu成反比关系.

图7 (a)均方根温度剖面分布与温度边界层定义示意图;(b)不同 1/Ro 下,温度边界层厚度随 Ra 变化规律Fig.7.(a) The variation of T* profile with Ra;(b) the thickness of the thermal boundary layers with Ra and 1/Ro.

图8为将不同算例的T*曲线用δT与max(T*)进行归一化后的结果.在保持 1/Ro恒定的条件下,图8 的每个子图中,Ra从 3×106增加至 3×109,提升了3 个数量级.结果显示,当旋转作用较小时(1/Ro≤1),即使Ra增加了3 个数量级,不同Ra的T*归一化曲线相互重合得较好,表现出较强的自相似性.当旋转作用达到最大时(1/Ro=10),在边界层中与靠近边界层的部分((Rθ/δT)≤2),不同Ra下的T*归一化曲线依然重合得较好.产生这种现象的原因在于科里奥利力在赤道附近作用微弱,如3.2 节所讨论.在离边界层较远的区域(2 ≤Rθ/δT≤4),除Ra=3×106外,其他Ra下的T*归一化曲线依然展现出较好的自相似性.在离边界层非常远的区域(Rθ/δT≥4),只有Ra=3×108与Ra=3×109的T*归一化曲线重合得较好.由此可以发现,在科里奥利力作用明显的条件下,随着Ra上升,T*的归一化曲线的重合度逐渐增加;而当纬度越靠近赤道,T*曲线的重合度越明显.在较强的旋转作用下,图中Ra=3×106在高纬度的差异可能是由于随着Ra的降低,其羽流尺度相对更大,当羽流延展到远离边界层区域时,科里奥利力作用逐渐增加,从而造成了羽流结合科里奥利力导致的温度分布变化.随着Ra上升,热羽流尺寸相对减小,这一耦合科里奥利力篡改温度分布的能力减弱,因而曲线的重合度更高.这些都表明,科里奥利力项的作用会改变T*的空间分布.另外,图8的结果也表明,m ax(T*) 是一个非常合适的归一化因子.m ax(T*) 直接描述了温度脉动的强弱,与Ra相关,于是可以较好地表征T*的自相似特性[49].

图8 不 同 1/Ro 和 Ra 条件下,使 用 δT 与 m ax(T*) 进行归一化之后的 T*曲线 (a) 1/Ro=0;(b) 1/Ro=0.1;(c) 1/Ro=1;(d)1/Ro=10Fig.8.The RMS temperature distribution normalized by δT and m ax(T*) for the different 1/Ro and Ra : (a) 1/Ro=0;(b) 1/Ro=0.1;(c) 1/Ro=1;(d) 1/Ro=10.

图9 采用对数坐标轴显示了不同Ra下T*在边界层外侧的随θ变化的曲线,每个子图的结果具有不同的 1/Ro并且使用幂律函数对T*进行了拟合,拟合结果用不同线型的黑线展示.对比拟合结果与T*曲线可以发现,在本文研究的罗斯比倒数范围内(0 ≤1/Ro≤10),边界层外的T*(θ) 都具有幂律函数的形式.在 1/Ro保持不变的条件下,尽管Ra变化范围达到了3 个数量级,但T*都在一定纬度范围内近似满足T*～θξ.当肥皂泡旋转作用较弱时(1/Ro≤1),近似满足T*～θξ的区间超过了一个数量级,拟合指数ξ与非常接近.但当肥皂泡高速旋转时(1/Ro=10),近似满足T*～θξ的区间长度缩短,小于一个数量级,并且ξ与的差距明显增加.对于所有算例,在θ≥1 时,流体非常靠近肥皂泡顶,脉动温度快速下降[14].于是T*(θ)不再具有幂律函数的形式.肥皂泡上的T*在边界层外的空间分布与RB 对流也具有一定的相似性.RB 对流中,边界层外可以分为对流区(或混合区)和对流核心区[58].Adrian[53]由惯性力与黏性力相互平衡的假定,推导得出对流区T*(d) 具有幂律函数的形式:T*～.多项RB 对流的实验研究也发现了T*的空间分布曲线具有幂律函数形式.然而Adrian 基于惯性力与浮力平衡的假设预测T*(d)同样可以具有对数函数的形式T*～lnd,并在实验中得到验证.最近,何宇昊与夏克青[54]使用高精度的实验方法测量了RB 对流槽中剪切主导区与羽流主导区的T*(d).剪切主导区中,大尺度环流带来的横向平均速度使流场处于强剪切状态中.而羽流主导区中,大量的羽流相互混合聚集并从热边界上升至冷边界.他们的研究表明,剪切主导区内T*～,而羽流主导区内T*～lnd.比较RB 对流和肥皂泡模型的结果可以发现,肥皂泡模型与RB 对流在边界层外具有较大的物理差异.由于大尺度环流非常微弱,于是肥皂泡上不存在与剪切主导区中类似的流动状态.另一方面,羽流在肥皂泡的赤道上随机生成并不断运动,羽流的混合随机发生在肥皂泡上不同的区域,羽流主导区中类似的流动状态也不存在于肥皂泡上.但是,肥皂泡上添加了外摩擦项,引入了较大的外摩擦黏性力;而随着θ的上升,在肥皂泡纬度方向上重力的投影越来越小.这些都暗示T*～θξ的分布规律可能取决于惯性力与黏性力相互平衡机制,但是详细物理图像还需要更加深入的研究.

图9 不同 1/Ro 和 Ra 条件下,T* 的剖 面曲线 (a) 1/Ro=0;(b) 1/Ro=0.1;(c) 1/Ro=1;(d)1/Ro=10Fig.9.The profile of T* for the different 1/Ro and Ra : (a) 1/Ro=0;(b) 1/Ro=0.1;(c) 1/Ro=1;(d) 1/Ro=10.

3.4 纬度脉动速度均方根的空间分布与黏性边界层

图10 不 同 Ra 与 1/Ro 条件下,均方根纬度速度 纬度剖面曲线 (a) Ra=3×106;(b) Ra=3×107;(c) Ra=3×108;(d)Ra=3×109Fig.10.The profile of the RMS velocity in the latitude direction for the differnet Ra and 1/Ro : (a) Ra=3×106;(b)Ra=3×107;(c) Ra=3×108;(d) Ra=3×109.

当肥皂泡上流体的湍流运动得到充分发展之后,流场的平均速度相比脉动速度非常微弱[14].之前的章节已经证明,纬度方向的平均速度〈uθ〉=0;静止肥皂泡上的大尺度环流非常微弱,其经度方向的平均速度同样非常微弱[14].因而,不便于通过平均速度场来定义黏性边界层.本文参考RB 对流的相关研究文献[50,51],基于纬度方向脉动速度的均方根来定义黏性边界层厚度δu.图11的左侧示例了通过曲线来定义黏性边界层厚度的方法.对于肥皂泡对流模型,随θ变化的曲线在靠近赤道的部分具有线性特征.将线性部分做拟合线并延长,当拟合线达到脉动纬度速度的均方根最大值 m ax() 时,纬度为θu.θu距离赤道的大地距离R×θu就是黏性边界层的厚度δu.

图11 (a)黏性边界层厚度 δu 的定义方法(示意算例Ra=3×107,1/Ro=0);(b)不 同 1/Ro 下 δu 随Ra变化规律Fig.11.(a) The definition of the viscous boundary layer thicknesses δu with the example case of Ra=3×107 and 1/Ro=0;(b) the variation of δu with Ra for the different 1/Ro.

图11 的右边给出了δu随Ra的变化规律.对于肥皂泡热对流模型,δu与Ra的关系呈现幂律函数的形式.而随着Ra增加,黏性边界层的厚度快速减小.对于静止肥皂泡,1/Ro=0,δu近似满足标度率δu～Ra-0.2.当 1/Ro≤1 时,肥皂泡上旋转作用较弱,δu随 1/Ro的变化可以忽略,此时δu随Ra变化的标度率保持恒定.当 1/Ro=10 时,肥皂泡高速旋转,Ra相对较低时(Ra≤3×108),δu明显下降.但随着Ra升高到 3×109,1/Ro对δu的影响可以忽略.结果,当 1/Ro=10 时,δu随Ra变化的标度指数从-0.2 变为-0.14.RB 对流中,由于大尺度环流的存在,使用水平脉动速度的均方根曲线来定义δu.在早期的文献中,Xin 和Xia[59]首先在圆柱形对流槽中测量得到δu～Ra-0.25,而Qiu 和Xia[60]测得δu～Ra-0.22.这与肥皂泡中的结果较为接近.最近,孙超等[50]在使用PIV 方法测量了矩形对流槽的黏性边界层,并得到δu～Ra-0.37.

3.5 拟热能耗散率和动能耗散率

令单位质量流体所包含的拟热能和动能分别为ET=T2/2与Ek=.将能量方程和动量方程两边分别乘以T或者u,就可以得到ET与动能Ek的控制方程:

图12 不同 Ra 与 1/Ro 条件下,拟热能内耗散率 在纬度方向的分布 (a) Ra=3×106;(b) Ra=3×107;(c)Ra=3×108;(d) Ra=3×109.子图为边界层附近的放大Fig.12.The distribution of the internal thermal energy dissipation rate in the latitude direction for the different 1/Ro and Ra :(a) Ra=3×106;(b) Ra=3×107;(c) Ra=3×108;(d) Ra=3×109.The insets are the zoom-in for the boundary layers.

随着Ra的上升,δT快速减小,温度边界层中和随θ下降得更快.肥皂泡在边界层中的平均温度梯度最大[14],于是拟热能内耗散集中在温度边界层中,而在肥皂泡的其他区域,拟热能内耗散较小.在边界层中,科里奥利力项的作用非常小,1/Ro的上升对拟热能耗散率的影响可以忽略.

图14 给出了不同Ra和 1/Ro参数组合下,和随纬度θ变化的曲线,并且使用黑色和红色虚线表示不同 1/Ro下黏性边界层的厚度δu.因为1/Ro≤1 的条件下,δu随 1/Ro的变化可以忽略,于是图14 使用黑色虚线标识δu(1/Ro=0),而使用红色虚线标识δu(1/Ro=10).可以看出,随着Ra数上升,δu(1/Ro=10) 与δu(1/Ro=0)之间的差别不断减小,最终在Ra=3×109时相互重合.从图14 还可以发现动能外耗散率比内耗散率要小一个量级,这暗示肥皂泡与外部的相互作用不会给湍流运动产生最本质的影响.此外,〈Tu3〉在边界层中随纬度快速增长,达到最大值后缓慢下降,在肥皂泡顶处为0.且随着Ra的上升,边界层变薄,〈Tu3〉的峰值越来越靠近赤道.在 1/Ro≤0.1 的条件 下,〈Tu3〉随 1/Ro并没有变化.而 1/Ro=1 时,〈Tu3〉在黏性边界层外随θ的上升有更明显的下降.在1/Ro=10 时,〈Tu3〉在黏性边界层外随θ的增加下降非常快,当θ≥45°时,〈Tu3〉的值约等于0.而在黏性边界层 内,〈Tu3〉的变化更加复 杂.在Ra=3×106或Ra=3×107的情况下,〈Tu3〉在边界层内的峰值有明显的上升.但随着Ra增加到3×108或 3×109时,〈Tu3〉峰值不再随 1/Ro变化.

图13 不同 Ra 与 1/Ro 条件下,拟热能内耗散率 在纬度方向的分布 (a) Ra=3×106;(b) Ra=3×107;(c)Ra=3×108;(d) Ra=3×109.内插图为边界层附近的放大Fig.13.The distribution of the internal thermal energy dissipation rate in the latitude direction for the different 1/Ro and Ra :(a) Ra=3×106;(b) Ra=3×107;(c) Ra=3×108;(d) Ra=3×109.The insets are the zoom-in for the boundary layers.

动能内耗散率的平均在边界层中有小幅的增长,之后随着θ的增加快速下降,在黏性边界层中形成一 个尖锐的峰值.当θ≥30°时,变 得非 常接 近于 0.随 着Ra的增 加,在黏性边界层中随θ上升下降更加剧烈,峰值变得越发尖锐,且不断接近赤道.图14 中可以观察到低速旋转(1/Ro≤ 1)对在黏性边界层内外的分布没有明显影响.但当肥皂泡高速旋转时(1/Ro=10),在黏性边界层中的峰值内有明显的下降,但在黏性边界层外又有轻微的增加.综合图14 中的结果可以发现,在黏性边界层中的值远大于边界层之外的值.

图14 不同 Ra 与 1/Ro 条件下,动能耗散率以及浮力项在 纬度 方向 的分 布 (a),(b),(c) Ra=3×106;(d),(e),(f)Ra=3×107;(g),(h),(i) Ra=3×108;(j),(k),(l) Ra=3×109.内插图为边界层附近的放大Fig.14.The distribution of the kinetic energy dissipation rate and the buoyancy term in the latitude direction for the different 1/Ro and Ra : (a),(b),(c) Ra=3×106;(d),(e),(f) Ra=3×107;(g),(h),(i) Ra=3×108;(j),(k),(l) Ra=3×109.The insets are the zoom-in for the boundary layers.

4 结论

通过简单加热肥皂泡可以获得一个新的热对流模型,该模型具有典型的准二维球面特征,相比RB 对流不存在由几何边界奇点形成的角涡,其浮力作用也随纬度不同而变化.本文提出了加热肥皂泡的数值模拟方法,并通过直接数值模拟研究了不同Ra与 1/Ro下旋转肥皂泡上的边界层特征和能量耗散规律.结合肥皂泡上的流动特征提出了相应的Nu,Re与温度和黏性边界层厚度的定义.其中δT定义为m ax(T*) 所在纬度距离赤道的大地距离,而δu是曲线近赤道线性部分延长至 m ax() 时所处纬度距离赤道的大地距离.DNS 数值结果表明加热肥皂泡上的湍流热对流具有以下特点.

1)当肥皂泡无旋转运动时,Nu与Re与Ra的标度关系为Nu～Ra0.30和Re～Ra0.48;尽管热对流模型与RB 对流不同,但其标度关系与标度指数相近,这进一步表明这一标度关系对热对流的普适性.旋转对Nu的影响可以忽略,但当旋转作用很强时(1/Ro=10),旋转在低Ra时对Re有一定的抑制,使得Re与Ra的标度关系变化为Re～Ra0.50.

2)温度边界层的厚度δT随Ra的标度关系为δT～Ra-0.32,且不随罗斯比数的变化而变化.在温度边界层外侧,均方根温度的剖面曲线与纬度θ的近似关系为而随着肥皂泡旋转作用的增加,满足标度关系的纬度区间逐渐缩短.

3)黏性边界层的厚度δu随Ra的标度关系为δu～Ra-0.20.只有当旋转作用足够强(1/Ro=10),并且Ra相对较小(Ra≤3×108)时,δu出现明显下降,使得与Ra标度关系变为δu～Ra-0.14.