摘要:采用高效并行的直接数值模拟方法计算了三个Pr数的系列Ra数的二维热对流. 对所有计算结果的平均场温度边界层特性进行研究, 采用考虑脉动作用的双参数温度边界层理论对温度边界层剖面进行拟合, 得到了拟合参数a和参数c的分布. 参数a决定了温度剖面的基本特性, 而参数c起到了对温度剖面外区的修正作用, 使得在5个边界层厚度温度边界层剖面的计算值与理论解符合良好. 参数c的变化特性与参数a的相反, a值增大则c值减小. 不同的Pr数拟合参数a随Ra数变化具有不同的分布特性, 但都存在突然减小的间断. 参数a的间断变化是由湍流热对流流动状态由椭圆形大尺度环流突变到圆形大尺度环流造成的, 且随着Pr数变大间断点的特征Ra数增高. 相同Ra数时Pr数越大温度剖面的拟合参数a值越小, 表明温度边界层中脉动的影响越小. 传热特性$Nu/R{a^{0.3}}$、表征羽流运动特性的大尺度环流路径周长, 及温度边界层拟合参数a三者随Ra数的变化都存在转折或间断, 且对应相同的特征Ra数, 显示三者具有很好的相关性并与流动形态变化直接关联.
关键词: 二维湍流热对流/
温度边界层剖面/
间断特性/
传热特性/
大尺度环流路径周长

English Abstract

全文HTML

--> --> -->

热对流现象通常指在空间中通过流动介质传递热量, 普遍存在于自然界和工业界中. 从众多热对流现象中抽取出来的Rayleigh-Bénard(RB)热对流系统是研究热对流问题的经典模型之一^[1]. RB热对流系统由一个充满流体介质的封闭腔体组成, 下底板恒温加热, 上顶板恒温冷却, 流体介质在热浮力驱动下形成对流运动. 在RB系统的研究中, 一个核心问题是湍流流动的热量运输机制. 2000年Grossmann和Lohse^[2,3]提出著名的GL理论, 对系统的传热特性随Ra数和Pr数的变化规律给出了很好的预测.
湍流热对流中温度边界层特性是研究的重点问题之一. Lui和Xia^[4]对三维圆柱RB系统的温度边界层空间结构进行研究, 发现不同Ra数下温度边界层沿横向呈现V形分布. Zhou和Xia^[5]在实测瞬时速度边界层时, 引入动态重构法, 发现平均速度剖面分布与Prandtl-Blasius层流边界层剖面一致. 文献[6-8]在进行极高Ra数的系列实验中发现经典区间和终极区间的温度剖面都存在对数律. Zhou等^[9,10]通过直接数值模拟(DNS)方法对二维RB系统进行模拟, 分别讨论了速度剖面和温度剖面特性. Stevens等^[11]对三维圆柱系统进行数值模拟, 指出温度边界层更容易受到速度边界层引起的体积波动影响. Shishkina等^[12]在层流假设下求解出热对流温度边界层方程的PBP (Prandtl-Blasius-Pohlhausen)解, 发现与实验和数值计算的温度边界层剖面都有较大的差别. Shishkina等^[13,14]考虑加入温度脉动对边界层的作用, 引入涡热扩散系数($ {\kappa }_{t} $), 推导出新的温度边界层理论剖面, 包括单参数和双参数解, 大大改善了理论解与实际温度剖面的符合程度. 何鹏等^[15]研究了二维RB系统中温度边界层拟合参数的分布特性. 黄茂静和包芸^[16]对二、三维RB系统进行了对比研究, 讨论了流态特征、近底板温度分布和温度边界层等特性. Wang等^[17]从温度边界层方程出发, 对温度脉动边界层方程及其剖面分布展开研究. 包芸等^[18]进行了Ra = 10¹⁰系列Pr数的计算, 发现温度边界层厚度和Re数随Pr数变化都存在一定的标度率关系.
本文采用并行直接求解方法(PDM-DNS)计算了三个典型Pr数系列Ra数的二维RB热对流. 采用Shishkina等^[13,14]提出的双参数温度边界层剖面理论, 对计算的热对流温度边界层剖面进行拟合, 得到拟合参数的分布特性, 对拟合参数a值随Ra数和Pr数变化规律进行讨论. 并结合已有研究成果, 讨论温度边界层拟合参数a值、传热Nu数及羽流运动特性三者随Ra数的典型变化规律.

2

-->

2.1.DNS数值模拟

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

$\left\{\begin{aligned} &\nabla \cdot {{u}}=0,\\ &\frac{\partial {{u}}}{\partial t}+\left({{u}}\cdot \nabla \right){{u}}=-\nabla p+\frac{1}{\sqrt{Ra/Pr}}{\nabla }^{2}u+\theta {{k}},\\ &\frac{\partial \theta }{\partial t}+\left({{u}}\cdot \nabla \right)\theta =\frac{1}{\sqrt{Ra \cdot Pr}}{\nabla }^{2}\theta,\end{aligned}\right.$

其中u为无量纲速度矢量, θ为无量纲温度, t为无量纲时间, $ {{k}} $为单位垂向矢量, $ p $为压力. 无量纲参数$ Ra=\left(\beta g\Delta \theta {H}^{3}\right)/\left(\kappa \upsilon \right) $为瑞利数, $ Pr=\upsilon /\kappa $为普朗特数, $ \varGamma =W/H $反映了对流系统的几何尺寸. $ \beta $为热扩散系数, $ g $为重力加速度, $ \Delta \theta $为上下壁面温差, $ H $为系统装置的高度, $ W $为宽度, $ \upsilon $为运动黏性系数, $ \kappa $为热扩散率. 计算的边界条件为四壁无滑移, 侧壁绝热, 上、下底板恒温, 分别为–0.5和0.5.
本文采用PDM-DNS方法^[19]对二维RB热对流系统进行直接数值模拟. 利用投影法求解热对流方程, 其中的压力泊松方程求解需要全流场联立, 利用FFT解耦泊松方程, 结合并行三对角占优 (PDD)算法并行求解三对角泊松方程. 使用“天河二号”超级计算机计算了Pr = 0.7, 4.3和20.0, $ 1\times {10}^{7}\leqslant Ra\leqslant 5\times {10}^{12} $的51个算例. 获取系列RB热对流时间平均场数据, 用于研究一系列湍流热对流的传热特性、温度边界层特性、流动特性以及其他物理特性.
2 -->

2.2.温度边界层双参数拟合原理

温度边界层特性是湍流热对流的主要研究课题之一, 它的研究成果对理解湍流热对流的传热特性, 以及增强和控制传热效率等研究都具有重要的价值.
针对热对流温度边界层的理论研究工作, Shishkina等^[13]在原有层流假设下, 热对流温度边界层PBP理论解^[12]的基础上, 推导湍流温度边界层时加入脉动的作用, 并引入涡热扩散系数$ {\kappa }_{\rm{t}} $和相似变量$ \xi \left(x, z\right)=z/\lambda \left(x\right) $, $ \lambda \left(x\right) $为局部边界层厚度, 得到了温度边界层剖面的单参数控制理论. 受到线性速度边界层假设的约束, 认为温度边界层的单参数拟合不适用于小Pr数的情况. 何鹏等^[15]应用这一理论对系列的二维湍流热对流温度边界层剖面进行了研究, 发现新的理论相较于PBP解在拟合温度边界层的结果上有了显著的改善, 并且拟合参数有明显的变化规律.
Shishkina等^[14]重新改进对流动的假设, 引用Prandtl的混合长度理论模型和新的涡热扩散系数分布关系

${({\kappa }_{\rm{t}}/\kappa )}_{\xi }\approx KB\psi, \;\;\; {\left(\frac{{\kappa }_{\rm{t}}}{\kappa }\right)}_{\xi }\approx \frac{3{a}^{3}{\xi }^{2}}{1+{b}^{2}{\xi }^{2}},$

推导出全Pr数下温度边界层方程

$\theta \left(\xi \right)=\frac{1}{b}{\int }_{0}^{b\xi }{\left[1+\frac{3{a}^{3}}{{b}^{3}}(\eta -{\rm{arctan}}\eta )\right]}^{-c}{\rm{d}}\eta ,$

其中a, b, c为常数, 但这三个常数并不完全独立.
由温度边界层分布的边界条件得知, 当$ \xi =0 $时$ \theta \left(0\right)=0 $, 当$ \xi \to {\rm{\infty}} $时$ \theta \left(\infty \right)=1 $. 将边界条件代入(3)式, 可以很容易导出三个常数a, b, c满足下列函数关系:

$b={\int }_{0}^{\infty }\left[1+\frac{3{a}^{3}}{{b}^{3}}\left(\eta -\arctan\eta \right)\right]^{-c}{\rm{d}}\eta,$

(3)式和(4)式中的温度边界层剖面方程包含了两个可调参数. 由此, 建立了温度边界层剖面的双参数控制理论. 与单参数理论相比, 双参数理论中的速度分布突破了线性假设, 适用于全Pr数, 因此可以用来分析Pr < 1的情况, 以及改善温度边界层的外区解.
双参数理论的温度边界层剖面拟合工作有两个参数可以调整, 加大了拟合工作的不确定性. 仔细分析单参数和双参数理论的假设发现, 在ξ → 0时, 两个理论的涡热扩散系数有同样的渐进解$ {\kappa }_{t}/\kappa \approx {a}^{3}{\xi }^{3} $. 而且即使是小Pr数情况, 在近壁面速度线性假设仍然成立. 由此在进行双参数理论剖面拟合的过程中, 先采用单参数拟合确定参数a的值, 然后再调整c值拟合整个温度边界层剖面.
考虑到参数a和参数c的不同作用, 首先讨论参数a对温度边界层剖面理论解的影响. 图1(a)给出了不同参数a时理论温度剖面分布和PBP温度剖面, 给定的空间范围为5倍温度边界层厚度, 相似变量$ \xi \left(z\right)=z/\lambda $. 从图1(a)可以看到, 不同a值的温度剖面与PBP理论曲线在近壁面处一致, 而后在不到0.5个边界层厚度处逐渐分开. 随着控制参数a逐渐减小, 理论温度剖面形状以单调趋势不断上升, 逐渐接近PBP曲线. 由此表明, 参数a的值越小, 温度边界层中脉动的影响越小, 理论解越接近层流的PBP解.
图 1 (a) 不同控制参数a, 双参数理论温度剖面分布, 黑色实线为PBP理论解; (b) Pr = 4.3, Ra = 10⁹, 单、双参数温度剖面拟合
Figure1. (a) Two-parameter theoretical temperature profile distribution for different parameters a, the black solid line is the Prandtl-Blasius predictions; (b) Pr = 4.3, Ra = 10⁹, one and two-parameter temperature profile fitting.

双参数温度边界层剖面拟合的具体步骤为, 在确定参数a时采用Shishkina等^[13]中的单参数方法, 目标是修改参数a的值, 使温度剖面在2倍边界层厚度内理论解和计算值的误差达到最小值. 确定参数a后采用本文讨论的双参数理论, 调整参数c使温度剖面在4倍边界层厚度内的误差最小.
图1(b)给出了Pr = 4.3, Ra = 10⁹时双参数拟合和单参数拟合的温度剖面, 可见两者都远好于层流PBP理论解, 说明温度边界层湍流脉动的影响不可忽略. 双参数拟合的效果在5倍边界层厚度范围内好于单参数的拟合效果, 尤其是在温度边界层的外区. 依据理论及本文拟合过程, 参数a直接反映近壁面温度的变化情况以及温度边界层剖面的总体走向和基本特性, 参数c则修正了温度边界层剖面的外区解. 因此, 双参数拟合理论在流动假设上的改进是有效的.
本文采用双参数拟合理论研究三个Pr数的系列Ra数二维湍流热对流温度边界层特性, 分析讨论温度边界层拟合参数变化特性与物理流动特性之间的关系.

在RB对流系统中, 温度边界层拟合控制参数大小及变化规律反映了脉动对温度边界层的影响程度. 下面将详细讨论不同Pr数和Ra数下拟合控制参数值变化特性, 以及与温度剖面分布和流态特性之间的关系.
图2给出了Pr = 0.7系列Ra数的热对流温度边界层拟合参数a和参数c的变化特征. 可以看到, 在$ Ra < 1\times {10}^{9} $时参数a ≈ 2.2附近, 参数a随Ra的变化呈现出略为下降的特征, 而后出现明显的间断, $ Ra\geqslant 1\times {10}^{9} $之后参数a ≈ 1.6附近出现波动, 反映出此时的温度边界层剖面形状基本不变化. 参数c的值小于参数a的值, 且变化特性相反, $ Ra < 1\times {10}^{9} $时参数c ≈ 0.6附近呈现略为上升的变化特征, 同样有分布的间断, $Ra\geqslant 1\times {10}^{9}$之后参数c ≈ 0.7附近出现波动.
图 2 控制参数随Ra数分布
Figure2. Distribution of control parameter with Ra number.

拟合参数值的变化直接反映温度边界层剖面的变化特性. 图3给出Pr = 0.7时温度边界层剖面拟合线的分布, 明显可以看出, 在$ Ra < 1\times {10}^{9} $时温度剖面偏低并基本重合, 而$ Ra\geqslant 1\times {10}^{9} $时温度剖面偏高基本重合. 拟合参数的间断分布导致归一化的温度剖面分布出现了两组明显的不同形态.
图 3 Pr = 0.7时不同Ra数近底板温度剖面分布
Figure3. Temperature distribution close to the bottom plate with different Ra number for Pr = 0.7.

温度边界层理论解中参数a的值越小越接近层流PBP解. 但在Pr = 0.7低Ra数时, 参数a的值较大而高Ra数时较小, 而且出现间断分布, 与Ra数高湍流强度大的特点不一致. 由此, 温度边界层剖面形态的变化可能与热对流的流动形态发生根本性变化相关.
给出这两组典型的温度边界层剖面形态对应的不同热对流流动形态, 如图4所示^[20].
图 4 两个典型热对流平均温度场和流线图, 颜色表示温度分布, 带有箭头的曲线为流线, 其中红色流线代表大尺度环流路径周长^[20]　(a) Ra = 1 × 10⁸; (b) Ra = 1 × 10⁹
Figure4. The time-averaged temperature fields and streamlines in two typical RB convection, the temperature is coded in color, and the arrow indicates the flow direction. The red streamline represents the length of the large-scale circulation path^[20]: (a) Ra = 1 × 10⁸; (b) Ra = 1 × 10⁹.

图4给出两组形态变化中对应的典型Ra数热对流的平均温度场和流线, 其中红色流线是最大速度所在流线, 定义为大尺度环流特征周长. 图4(a)中较低Ra数时大尺度环流呈现椭圆状, 底板左右两侧角涡大小差异较大, 最大速度流线带起的边界附近高温层明显堆积, 使得温度在边界附近变化减缓. 图4(b)中较高Ra数时大尺度环流呈现为圆形, 下底板两侧角涡大小基本一致, 大尺度环流周长远离边界附近并且明显减小, 温度在边界附近变化很快. 两个典型热对流的流动形状完全不同. 因此, 在本文研究的Pr = 0.7的Ra数范围内, 不同的流动形态导致热对流的温度边界层产生了两种不同的剖面形状, 并可由剖面的拟合参数变化特性反映.
由此可见, 温度边界层剖面形状以及拟合参数的变化, 与热对流的流动形态变化有着密切的关联. 温度边界层拟合参数的变化特性可直接反映出这些典型的变化特征.
图5给出的是Pr = 4.3和20.0时的温度剖面拟合参数随Ra数的变化情况. 图5(a)显示的是Pr = 4.3时参数a和参数c随Ra数的变化过程. 与Pr = 0.7相比, 两个参数之间的数值差减小. 当Ra ≤ 5 × 10⁸时, 参数a随Ra数变化很小, a ≈ 1.5, 当5×10⁸ < Ra ≤ 5×10⁹时, 参数a随Ra增大而逐渐略为增大到$ a \!\approx \! 1.7 $, 在$Ra \!\approx\! 1 \!\times\! {10}^{10}$时参数a突然减小, 出现间断, 当Ra ≥ 1 × 10¹⁰时, 参数a相对比较稳定, a ≈ 1.5, 随着Ra数增高无明显波动. 参数c的值小于参数a的值, 且变化特性相反, c ≈ 0.85. 总体上参数a随Ra的变化幅度较小, 表明此时的温度剖面形状变化不大.
图 5 拟合参数a随Ra数分布　 (a) $Pr \!=\! 4.3$; (b) $ Pr \!=\! 20.0 $
Figure5. Distribution of fitting parameter a with Ra number: (a) Pr = 4.3; (b) Pr = 20.0.

图5(b)给出了Pr = 20时参数a和参数c随Ra数的分布图. 可以看出, 参数a的值在低Ra数时相对较小并随着Ra数增高而增加, 从a ≈ 1.0连续变化到 a ≈ 1.35, 表明此时的温度边界层剖面形状也是连续变化的. 从图2的结果可知, 大Pr数时温度边界层剖面存在从弱湍流化向强湍流化形态连续过渡的变化规律. 这种变化直到$ Ra\approx 1\times {10}^{11} $时参数a突然减小出现间断, 此时流动形态同样发生了变化. 参数c的变化与参数a的变化特性相反, 只是由于参数a减小而参数c增大, 在图中出现了两者相交的分布.
大Pr数温度边界层拟合参数的这种不同分布的物理含义, 还有待更多的研究.
下面结合三个典型Pr数, 讨论温度边界层拟合参数值随Ra数变化规律及间断位置特征.
图6给出了三组Pr数下拟合参数a随Ra数的变化分布, 其中小图是c值的分布. 可以看出, 相同Ra数时Pr数越大, 温度边界层的拟合参数a值越小, 表明温度边界层中脉动的影响越小. 不同Pr数的参数a随Ra数分布规律有明显差异. 总体变化趋势上, Pr = 0.7时参数a值由大变小, Pr = 4.3时参数a的变化在中间有明显抬高但总体呈现水平, Pr = 20时参数a值由小变大, 三者都存在突变式间断. 间断之后三者的拟合参数a值都呈现小幅波动, 变化幅值不大. 小图中c值的变化特点与a值成相反趋势, Pr数越大c值也越大, 随Ra数的变化为a值增加c值减小.
图 6 控制参数a随Ra数的分布
Figure6. Distribution of control parameter a with Ra number.

总结以上三个典型Pr数的温度边界层剖面拟合特性发现, 拟合参数ɑ虽然随Ra数的变化不完全一致, 但都存在一个共同的特性, 即出现了间断的现象, 并与热对流的大尺度环流形态从椭圆变为圆形的流态特性变化相关联.
拟合参数及其变化规律在温度边界层特性中的作用, 还需要更深入的研究.

Nusselt(Nu)数反映RB系统的整体传热. 在较大范围内, GL理论预测了传热Nu数与Ra数和Pr数的关系. 包芸等^[21]发现二维热对流$ Nu/{Ra}^{0.3} $数随Ra数变化具有两个转折点, 并在极高Ra数时其变化与GL理论的预测倍数线走向一致, 而反映羽流运动特性的大尺度环流路径周长变化和$ Nu/{Ra}^{0.3} $数随Ra数的转折有很好的相关性.
进一步分析温度边界层拟合参数a值与大尺度环流路径周长和传热Nu数随Ra数变化规律, 发现三个典型Pr数的温度边界层拟合参数a值的间断变化对应的Ra数与其他两者的变化特征都有着很好的对应关系.
图7分别给出了Pr = 0.7, 4.3, 20.0时大尺度环流路径周长、传热Nu数及参数a值随Ra数的变化, 红色圆点、深蓝色三角和绿色方块分别是Pr = 0.7, 4.3, 20.0的计算结果. 可以看出, 不同Pr数的三个物理量, 大尺度环流路径周长、传热Nu数及温度边界层拟合参数a值在红线标出的特征Ra数都具有相同的变化特性. 在红线处, 大尺度环流路径周长随着Ra数变化出现突然的减小, 导致传热Nu数随Ra数变化出现转折. 包芸等^[21]在分析Pr = 0.7和4.3的结果时指出, 这是由于二维热对流的流动形态由椭圆变为圆形造成大尺度环流路径周长变化产生的. 同样本文Pr = 20的结果也有类似现象. 大尺度环流路径周长的突然减小由羽流运动特性的变化引起, 造成了传热Nu数随Ra数变化的标度关系改变. 而这个流动形态的突然变化在热对流的温度边界层中也引起相应变化, 温度剖面的拟合参数a随Ra数的变化特性出现间断, 即温度剖面的形状变化出现不连续性. 三个物理量随Ra数变化具有很好的相关性.
图 7 大尺度环流路径周长、传热Nu数及参数a随Ra数变化, 其中红色空心圆点和蓝色空心三角是包芸等^[21]的计算结果, 黑色虚线是GL理论预测倍数线, 红线是特征Ra数位置　(a) Pr = 0.7; (b) Pr = 4.3; (c) Pr = 20
Figure7. The variation of large scale circulation path length, Nu number, and parameter a with Ra number, among which the red hollow dots and blue hollow triangles are calculated by Bao Yun^[21], etc. The black dotted line is the multiple predicted by the GL theory, the red line is the characteristic Ra number position: (a) Pr = 0.7; (b) Pr = 4.3; (c) Pr = 20.

综上所述, 二维湍流热对流流动形态的变化导致热对流的诸多物理量都发生了改变. 不同Pr数下大尺度环流路径周长、传热Nu数及温度边界层拟合参数a值随Ra数变化过程中, 都存在相关的典型变化, 其对湍流热对流特性的研究价值将进行更多深入的探讨.

温度边界层特性是Rayleigh-Bénard热对流研究的重点问题之一, 与传热机制紧密相关. 本文采用高效并行的DNS方法, 计算了Pr = 0.7, 4.3, 20.0在系列Ra数下的二维方腔热对流, 并采用Shishkina提出的双参数温度边界层理论, 对所有Pr数和Ra数的温度边界层剖面进行拟合, 得到了拟合参数a和参数c的分布. 通过分析和研究不同Pr数系列Ra数的控制参数a值、温度剖面特性、流态特性及传热特性, 得到的结论如下.
1)不同的Pr数温度剖面拟合参数a随Ra数变化具有不同的分布特性. Pr = 0.7, 参数a值基本分布在2.2和1.6两个值处, 中间出现突变减小, 由此造成温度剖面形成两组不同的剖面形状, 并且与热对流的流动特性大尺度环流的形状相关联. Pr = 4.3, 参数a值在1.4附近微小波动然后逐渐增加到1.6, 同样出现突降至1.4附近并保持基本不变, 此时的温度剖面将重叠在一起. Pr = 20, 参数a值随Ra数增高而增大, 直到Ra数增加到1 × 10¹¹时, 参数a同样突变减小而后基本不变. 参数c的变化特性和参数a的变化特性相反.
2)三个Pr数的温度剖面拟合参数a随着Ra变化, 在总体趋势上, Pr = 0.7时参数a值由大变小, Pr = 4.3时参数a的变化在中间有明显抬高但总体呈现水平, Pr = 20时参数a值由小变大. 三者都存在突变减小间断, 且随着Pr数变大, 间断点的特征Ra数增高. 相同Ra数时, Pr数越大温度剖面的拟合参数a值越小, 表明温度边界层中脉动的影响越小.
3)在同一Pr数下, 传热特性的$ Nu/{Ra}^{0.3} $、反映羽流运动特性的大尺度环流路径周长, 以及温度边界层特性的参数a都存在共同的特征, 由于二维热对流的流动形态由椭圆变为圆形造成大尺度环流路径周长变化, 三者在随Ra数的变化中出现明显的转折或间断, 且对应相同的特征Ra数. 三个典型Pr数系列Ra数的计算结果表明, 都存在传热、大尺度环流路径周长和温度剖面拟合参数a在同一特征Ra数发生转折或间断的特性, 三个物理量具有很好的相关性.

1)不考虑脉动影响的温度边界层PBP理论解
Shishkina等^[12]从无限长平板的二维流动出发, 得到温度方程

$\small\frac{\partial \theta }{\partial t}+u\frac{\partial \theta }{\partial x}+v\frac{\partial \theta }{\partial z}={\rm{\kappa}}\left(\frac{{\partial }^{2}\theta }{\partial x\partial x}+\frac{{\partial }^{2}\theta }{\partial z\partial z}\right),\tag{A1}$

考虑定常问题消去时间项$ \partial \theta /\partial t $, 并且假设温度在x方向的二阶偏导远小于z方向的二阶偏导, 即$ \left|\right({\partial }^{2}\theta )/\partial x\partial x|\ll \left|\right({\partial }^{2}\theta )/\partial z\partial z| $, 方程(A1)可化为

$u\frac{\partial \theta }{\partial x}+v\frac{\partial \theta }{\partial z}=\kappa \frac{{\partial }^{2}\theta }{\partial z\partial z}.\tag{A2}$

通过雷诺分解, 将瞬时量分解为时间平均项和时间脉动项. 即$u = \left\langle u \right\rangle + u'$, $v = \left\langle v \right\rangle + v'$和$\theta = \left\langle \theta \right\rangle + \theta '$, 并引入脉动连续性方程. 代入(A2)式后沿时间平均, 可得到

$\left\langle {{u}} \right\rangle \frac{{\partial \left\langle {{\theta }} \right\rangle }}{{\partial {{x}}}} + \left\langle {{v}} \right\rangle \frac{{\partial \left\langle {{\theta }} \right\rangle }}{{\partial {{z}}}} + \frac{{\partial \left\langle {{{u'\theta '}}} \right\rangle }}{{\partial {{x}}}} + \frac{{\partial \left\langle {{{v'\theta '}}} \right\rangle }}{{\partial {{z}}}} = {{\kappa }}\frac{{{\partial ^2}\left\langle {{\theta }} \right\rangle }}{{\partial {{z}}\partial {{z}}}},\tag{A3}$

进一步假设$\left| {\partial \langle {{u'\theta '}}\rangle /\partial {{x}}} \right| \ll \left| {\partial \langle {{v'\theta '}}\rangle /\partial {{z}}} \right|$, 有

$\left\langle {{u}} \right\rangle \frac{{\partial \left\langle {{\theta }} \right\rangle }}{{\partial {{x}}}} + \left\langle {{v}} \right\rangle \frac{{\partial \left\langle {{\theta }} \right\rangle }}{{\partial {{z}}}} + \frac{{\partial \left\langle {{{v'\theta '}}} \right\rangle }}{{\partial {{z}}}} = {{\kappa }}\frac{{{\partial ^2}\left\langle {{\theta }} \right\rangle }}{{\partial {{z}}\partial {{z}}}},\tag{A4}$

忽略脉动的影响, 通过直接对方程(A4)进行求解, 可以推导得到经典的PBP理论解^[12], 解的形式如下:

$\theta \left(\xi \right)={\int }_{0}^{\xi }\exp[-{\Gamma }^{\omega }(1+{\omega }^{-1}\left){\eta }^{\omega }\right]{\rm{d}}\eta,\tag{A5}$

其中$\Gamma$为伽马函数, $ \omega $为常数. 当$ Pr\ll 1 $时, $ \omega =2 $; 当$ Pr\gg 1 $时, $ \omega =3 $. 当为任意Pr数时, $ 2\leqslant \omega \leqslant 3 $.
2)考虑脉动的单参数温度边界层理论
与PBP理论不同, Shishkina等^[13]考虑脉动的影响, 引入涡热扩散系数

${\left\langle {v'\theta '} \right\rangle _t} = - {\kappa _t}\frac{{\partial \left\langle \theta \right\rangle }}{{\partial z}},\tag{A6}$

代入方程(A4)整理后得

$\left\langle u \right\rangle \frac{{\partial \left\langle \theta \right\rangle }}{{\partial x}} + \left( {\left\langle v \right\rangle - \frac{{\partial {\kappa _t}}}{{\partial z}}} \right)\frac{{\partial \left\langle \theta \right\rangle }}{{\partial z}} = \left( {\kappa + {\kappa _t}} \right)\frac{{{\partial ^2}\left\langle \theta \right\rangle }}{{\partial z\partial z}}.\tag{A7}$

引入相似变量$ \xi \left(x, z\right)=z/\lambda \left(x\right) $, 其中$ \lambda \left(x\right) $为温度局部边界层厚度. 当$ z\to \infty $时, $\left\langle u \right\rangle \left( {x, z} \right) \to {\left\langle u \right\rangle _0}$. 设平均速度的流函数$\varPsi = {\left\langle u \right\rangle _0}\lambda \left( x \right)\psi \left( \xi \right)$(其中$ \psi \left({\rm{\xi}}\right) $为流函数)可得

$\left\{ {\begin{aligned}&{\left\langle u \right\rangle = \frac{{\partial \varPsi }}{{\partial z}},}\\&{\left\langle v \right\rangle = - \frac{{\partial \varPsi }}{{\partial x}},}\\&{\left\langle \theta \right\rangle = {\theta _{{\rm{bot}}}} - \left( {\frac{\varDelta }{2}} \right)\theta \left( \xi \right).}\end{aligned}} \right.\tag{A8}$

同时, 给出温度和流函数的边界条件

$\left\{\begin{aligned}&\theta \left(0\right)=0,~ \frac{\partial \theta }{\partial \xi }\left(0\right)=1,~ \theta \left({\rm{\infty}}\right)=1,\\ & \psi \left(0\right)=0,~ \frac{\partial \psi }{\partial \xi }\left(0\right)=0,~ \frac{\partial \psi }{\partial \xi }\left({\rm{\infty}}\right)=1.\end{aligned}\right.\tag{A9}$

将平均速度的流函数公式(A8)式代入方程(A7)中得

$\left(\kappa +{\kappa }_{t}\right)\frac{{\partial }^{2}\theta }{\partial {\xi }^{2}}+\left[{\frac{\partial {\kappa }_{t}}{\partial \xi }+U}_{0}\psi \lambda \left(x\right)\frac{\partial \lambda \left(x\right)}{\partial x}\right]\frac{\partial \theta }{\partial \xi }=0,\tag{A10}$

令$ A=(\partial {\kappa }_{t}/\partial \xi )/\kappa $, $ B={\langle u\rangle }_{0}\lambda \left(x\right)(\partial \lambda (x)/\partial x)/\kappa $. 由于GL理论标度律假设, 存在局部边界层厚度$\lambda (x)\propto {(vx/{U}_{0})}^{1/2}\to \lambda \left(x\right)\propto {x}^{1/2}$, 代入B的定义式可得B为常数. 此时, 方程(A10)可化为以下形式:

$\left(1+{\kappa }_{t}/\kappa \right)\frac{{\partial }^{2}\theta }{\partial {\xi }^{2}}+\left[A+B\psi \right]\frac{\partial \theta }{\partial \xi }=0.\tag{A11}$

大量三维DNS的计算结果拟合得出, 近底板处$ {\kappa }_{t} $与$ \xi $存在标度律关系, 即$ {\kappa }_{t}\propto {\xi }^{3} $. 由壁面处所有脉动为0的条件, 可以进一步得到

${\kappa }_{t}\approx {a}^{3}{\xi }^{3}\kappa . \tag{A12}$

(A12)式是该理论的另一个重要假设, 称为涡热扩散系数分布假设. 将其代入方程(A10)可以得到

$\left(1+{a}^{3}{\xi }^{3}\right)\frac{{\partial }^{2}\theta }{\partial {\xi }^{2}}+\left(3{a}^{3}{\xi }^{2}+B\psi \right)\frac{\partial \theta }{\partial \xi }=0.\tag{A13}$

上述方程中包含了速度和温度的影响.
对于大Pr数, 热对流的温度边界层嵌入在速度边界层内, 流动速度变化假设为线性, 即$\left\langle u \right\rangle \propto z$. 化简后可以得到

$\psi \approx b{\xi }^{2},\tag{A14}$

其中$ b $为常数, 将(A14)式代入方程(A13)中, 并令$c=bB/(3{a}^{3})+1$, 方程(A13)的解为

$\theta \left(\xi \right)={\int }_{0}^{\xi }{(1+{a}^{3}\eta )}^{-c}{\rm{d}}\eta. \tag{A15}$

将边界条件$ \theta \left({\rm{\infty}}\right)=1 $代入(A15)式中, 得到常数a与常数c的关系为

$a=\frac{\varGamma \left(1/3\right)\varGamma (c-1/3)}{3\varGamma \left(c\right)}.\tag{A16}$

(A15)式和(A16)式称为考虑脉动影响的单参数理论剖面公式.