穿透对流是指热不稳定流体层中的对流穿透到相邻稳定层时发生的现象。在大多数实验室对流实验中，不稳定的流体层通常以刚性边界为界。然而，恒星对流区以稳定分层区域为界，了解对流穿过稳定层和不稳定层之间界面的穿透对流具有天体物理学意义^[1]。穿透对流通常在简化模型中研究，该模型通常是2D或3D简单几何体，其中充满了稳定(上方)和不稳定(下方)分层的流体。旋转Rayleigh-Bénard对流(Rayleigh-Bénard convection，RBC)是一个理想化的系统，通常用于描述地球和天体物理环境中的基本物理学，例如发生在恒星外层、海洋、大气层或卫星和行星的金属核心中的对流运动^[2-3]。旋转RBC系统通常简化为从下方加热、从顶部冷却，能够绕其垂直轴旋转的容器。然而，浮力和科氏力的联合作用可能会产生高度复杂的流动，流态及其性质在很大程度上取决于控制参数。当旋转RBC系统中的工作流体不是Oberbeck-Boussinesq(OB)流体(例如接近4 ℃的水)时，流体的密度随温度非线性非单调变化，在容器上方形成稳定层，并在浮力达到某一阈值时产生穿透对流。

前人的工作已经充分讨论了OB和non-Oberbeck-Boussinesq(NOB)流体从对流开始阶段到湍流的旋转RBC^[4-7]，以及2种不相溶流体或1种密度随温度非单调变化的流体的穿透对流^[8-11]。此前同时讨论旋转对流和穿透对流的大多数工作均以恒星对流为背景，并假设系统由2种互不相溶的流体填充，通过线性稳定性分析和实验研究，讨论在较低瑞利数(Ra)下，对流开始阶段的流动^[12-15]。截至目前，在中等和较高Ra下研究穿透性旋转RBC的工作尚属空白，本文采用直接数值模拟(direct numerical simulation, DNS)对该问题进行较为系统的研究，旨在探究施加轴向旋转对穿透RBC流态及传热的影响。

1 模型和数值方法

本文考察竖直圆筒内的旋转穿透RBC。圆筒尺寸固定，高$\hat{H}$ (本文中$\hat{·}$表示有量纲的物理量)，直径$\hat{D}$，宽高比$\mathit{Γ}=\hat{D} / \hat{H}=1$，下方恒温加热，上方恒温冷却，温度分别为$\hat{T}_{\text {heat }}$和$\hat{T}_{\text {cool }}$，侧壁绝热，旋转轴与圆筒对称轴重合。

在本文所研究的参数范围内，表征离心浮力作用效果的弗劳德数$F r \sim o\left(10^{-5}\right) \ll 1$，离心浮力效应可以忽略。使用的无量纲参数包括：普朗特数${Pr}=\hat{\nu} / \hat{\kappa}$，瑞利数$R a=\hat{g} \hat{\alpha} \hat{\Delta}^q \hat{H}^3 / \hat{\nu} \hat{\kappa}$，罗斯比数Ro= $\hat{U} / 2 \hat{\Omega} \hat{H}$，以及密度倒置参数θ_m= $\left(\hat{T}_m-\hat{T}_{\text {heat }}\right) /\left(\hat{T}_{\text {heat }}\right.$$\left.-\hat{T}_{\text {cool }}\right)$，其中$\hat{T}_m$为密度最大值对应的温度。

为简化模型，假设流体为不可压缩牛顿流体，忽略黏性耗散，仅在重力项中考虑密度对温度的依赖性，其余物性参数均设为常数。密度与温度之间的函数关系表示为$\hat{\rho}=\hat{\rho}_m(1-\hat{\alpha}(\hat{T}-$$\left.\left.\hat{T}_m\right)^q\right)$^[16]，其中$\hat{\rho}_m$=999.972 kg/m³为密度极大值，$\hat{\alpha}=9.297 \times 10^{-6}\left({ }^{\circ} \mathrm{C}\right)^{-q}$为热膨胀系数，$\hat{T}_m$=4.029 ℃为密度极值对应的温度，指数q=1.895。

以$\hat{H}$和$\hat{U}=\sqrt{\hat{g} \hat{\alpha} \hat{\Delta}^q \hat{H}}$分别作为长度、速度的特征尺度, 温度经由θ= $\left(\hat{T}-\hat{T}_{\text {cool }}\right) / \hat{\Delta}\left(\hat{\Delta}=\hat{T}_{\text {heat }}-\right.$$\left.\hat{T}_{\text {cool }}\right)$无量纲化，得到旋转参考系下的无量纲控制方程

$ \begin{gathered} \nabla \cdot \boldsymbol{u}=0, \\ \frac{\partial \boldsymbol{u}}{\partial t}+(\boldsymbol{u} \cdot \nabla) \boldsymbol{u}=-\nabla p+\sqrt{\frac{P r}{R a}} \nabla^2 \boldsymbol{u}-\frac{1}{R o} \boldsymbol{e}_z \times \boldsymbol{u}+ \\ \left(\theta-\theta_m\right)^q \boldsymbol{e}_z, \\ \frac{\partial \theta}{\partial t}+(\boldsymbol{u} \cdot \nabla) \theta=\frac{1}{\sqrt{P r \cdot R a}} \nabla^2 \theta . \end{gathered} $

(1)

其中: u为速度，t为时间，p为压力，θ为温度, e_z为轴向单位向量。求解域的热边界条件如上所述，速度场在边界上满足无滑移条件。本文固定冷水Pr为11.57，利用努塞尔数Nu评估对流系统的传热能力，定义为$N u=\sqrt{{Pr} \cdot R a}\left\langle u_z \theta\right\rangle_{V, t}-\left\langle\frac{\partial \theta}{\partial z}\right\rangle_{V, t}。$

本文使用有限体积方法进行数值求解，时间项采用二阶Crank-Nicholson格式进行离散，对流项使用二阶迎风格式进行离散，扩散项则采用二阶中心差分格式进行离散。该数值方案在冷水的旋转热对流研究中得到了较充分的验证^[16]。本文采用的初始条件为静态热传导(即速度为0，温度随高度线性变化)，并以较低参数的数值结果作为相邻高参数的初始条件。

2 结果分析与讨论 2.1 网格无关性检验

数值求解过程采用非均匀网格，沿水平方向采用O型网格，径向和轴向网格在壁面附近加密，网格示意图如图 1所示。对于无旋转的工况，网格分布满足Shishkina等^[17]提出的最小网格数要求。对于旋转穿透对流，涡柱的产生使得上壁面附近温度梯度增大，对上壁面附近的网格具有更高要求。因此以上、下壁面Nu一致作为网格充足的判断标准。表 1展示了Ra=1×10⁷, 1/Ro=6.0, θ_m=0.5时不同网格上下壁面平均努塞尔数($\overline{N_u}$)的计算值，其中全部网格均满足Shishkina等^[17]提出的最小网格数要求，综合考虑计算精度和效率，最终选用M2网格进行计算。

2.2 无旋转系统中冷水穿透对流的流动和传热特征

无旋转穿透RBC的流动和传热特性得到了较为广泛的实验^[18-19]和数值模拟研究^[20-23]，本节给出Ra∈[1×10⁴, 1×10⁸]，θ_m=0.0, 0.5时无旋转穿透对流的DNS结果，并将其作为在2.3节讨论旋转穿透RBC流动和传热特征的参考。

为了解不同Ra和θ_m下的流动结构，图 2展示了不同控制参数下的瞬时温度等值面。随着Ra的增加，流动从稳态流动转变为混沌态，并最终达到充分发展的湍流阶段。在稳态阶段，随着Ra缓慢增加，流动状态从静态热传导转变为对流辊，然后转变为六边形对流胞；当Ra进一步增加，流动进入非稳态区域，扰动的增长使六边形对流胞再次转变为对流辊，这一现象在20世纪通过稳定性分析得到了研究^[24-25]。如图 2所示，在非稳态对流区域，较强的浮力使流动具有更多的小尺度结构，在腔体中心区形成一个对流辊，表明存在大尺度环流(large-scale circulation，LSC)。LSC促进了腔体中心区域的流体混合，使得对流中心区域温度均匀分布，如图 3(a)所示。对于OB情况，中心区域的平均温度θ_c为(θ_top+θ_bottom)/2=0.5。但对于θ_m=0.0的冷水，由于密度随温度的非线性变化，θ_c≡〈θ|_z=H/2〉_{r, ϕ, t}略高于0.5。对于θ_m=0.5的冷水，由于密度随温度的非单调变化，上方稳定层抑制了冷边界与对流区的热交换，使得θ_c进一步上升。Wang等^[22]给出了二维方腔中冷水穿透对流中心温度随密度倒置参数的变化公式：θ_c=(1+θ_m²)/2。本文对于无旋转冷水穿透对流的中心温度在非稳态情况下与Wang等给出的结果基本一致。

如图 2所示，θ_m=0.0和0.5的流态随Ra的变化趋势具有一定的相似性，但θ_m=0.5的流态转变相对滞后。Sun等^[26]利用有效瑞利数Ra_E=Ra ·(1-θ_m)^q ·h_θm³解释流态转变的滞后现象，其中h_θm是密度倒置界面(即θ_m等值面)的平均高度。当θ_m≠0.0时，在给定θ_m的情况下，随着Ra的增大，越来越多的流体可以穿透进入上部稳定区域，稳定层越来越薄，密度倒置界面的平均高度h_θm随Ra单调增加，如图 4(a)所示。在高Ra下，h_θm接近1，稳定层对下方对流区流动的影响减弱，θ_m=0.0和0.5的流态更为接近。

这在一定程度上解释了穿透对流流态转变的滞后性，但并不全面。例如，当Ra=5×10⁶，θ_m=0.5时，其有效瑞利数约为Ra_E≈1.2×10⁶，但其径向速度及温度、速度脉动值都远小于Ra=1×10⁶，θ_m=0.0的情况，如图 3(b)~3(d)。我们认为上方稳定层造成了额外的动能耗散，是影响整体流动和传热特性的主要因素，在后续研究中值得进一步探究。

本文对流场的定量研究结果表明，不同于OB流体RBC的上下对称性，冷水的密度变化破坏了流场和温度场的上下对称性。当θ_m=0.0时，由于密度随温度的非线性变化使得上壁面附近密度变化相比下壁面更为剧烈，导致温度脉动在上壁面附近较强而在下壁面附近较弱，如图 3 (b)虚线所示。当θ_m=0.5时这种非对称性明显加剧，如图 3(b)实线所示。此外，稳定层额外的动能耗散使得θ_m=0.5的速度和速度脉动相比θ_m=0.0时明显减小，如图 3(b)~3(d)。

在无旋转热对流中，边界层厚度和Nu是描述流体运动和传热特性的重要参数。对于热边界层，通常根据斜率准则^[5]定义热边界层厚度，公式如下

$ \delta_{\mathrm{top}}^\theta=\frac{\theta_m-\theta_{\mathrm{top}}}{\left.\partial_z\langle\theta\rangle_{r, \phi}\right|_{z=1}}, \delta_{\mathrm{bottom}}^\theta=\frac{\theta_{\mathrm{bottom}}-\theta_m}{\left.\partial_z\langle\theta\rangle_{r, \phi}\right|_{z=0}} . $

(2)

对于速度边界层，Shishkina等^[17]认为利用水平速度或最大均方根速度都不能很好地反映速度边界层厚度，基于最大能量耗散的计算结果更为合理，将速度边界层厚度定义为最大能量耗散位置距上/下壁面距离的2倍，其中能量耗散率计算式为

$ \boldsymbol{\epsilon}_u^{\prime \prime}:=\left\langle\boldsymbol{u} \cdot \nabla^2 \boldsymbol{u}\right\rangle_{r, \phi, t} . $

(3)

图 4(b)~4(d)给出了Nu和边界层厚度$\delta_{\text {top/bottom }}^{\theta / u}$与Ra的标度关系。在本文研究的参数范围内, Nu~Ra^α和δ^θ~Ra^β的标度指数为α=0.307，β=-0.315，与Wang等^[22]得到的标度律基本一致。

2.3 旋转系统中冷水穿透对流的流动和传热特征

本小节以Ra=1×10⁷为例，利用DNS结果分析旋转系统中冷水穿透对流的流动和传热特征，并讨论密度倒置现象如何影响旋转对流。

2.3.1 流动结构和温度分布

图 5给出了不同Ro和θ_m时的瞬时温度场等值面。可以看出，随着1/Ro的增加，冷水呈现3种流态：大尺度环流(LSC)和羽流; 对流泰勒柱(convective Taylor columns，CTC)，也称为涡柱或Ekman涡; 以及涡柱崩溃, 即CTC崩溃。

正如2.2节关于无旋转冷水穿透对流研究得出的结论，θ_m=0.0的流动类型与OB流体相同，图 5(a)显示了低转速下的流态为LSC和羽流；随着1/Ro的增大，流态转变为CTC，在较高转速下CTC崩溃。当θ_m=0.5，即存在上方稳定层时，在中等转速下仍然存在涡柱结构，但是仅热羽流能够形成涡柱，冷羽流仍保持羽流的形态，这是由于冷流体(即θ < θ_m的流体)在下降过程中与热流体进行热交换并升温，当冷流体温度超过θ_m时，密度开始减小，阻碍冷流体继续下沉。

图 6(a)为沿垂向分布的时均温度剖面〈θ〉_{r, ϕ, t}，其形态在中低速旋转时几乎不发生变化(θ_m=0.0时，1/Ro≤8, θ_m=0.5时, 1/Ro≤2)，且中心温度与无旋转情况保持一致；而在快速旋转时，对流区存在非零温度梯度，中心温度明显升高。这是由于：当转速足够大时，流动几乎是二维的，在垂直方向上几乎不存在任何混合；混合只在水平方向发生涡合并时存在，因此不存在较为充分的三维混合。此外，浮升力和较大的科氏力在上方稳定区的作用效果均为抑制垂向流动，因此顶部冷壁面相比底部热壁面更难产生羽流，对流区热量积聚导致中心温度升高。图 6(b)为温度脉动剖面〈θ_rms〉_{r, ϕ, t}，当θ_m=0.0时，上下壁面附近均存在极大值，且2个极大值随1/Ro的增加同步变化，均为先增大后减小；当θ_m=0.5时，上下壁面的极大值变化趋势不同步，上壁面脉动极大值在1/Ro=2.0时达到最大，下壁面在1/Ro=6.0时达到最大；无论θ_m为何值，有旋转时的温度脉动始终大于无旋转情况。此外，1/Ro=0.0时，θ_m=0.5的温度脉动始终小于θ_m=0.0；而1/Ro≠0.0时，θ_m=0.5在下壁面附近的温度脉动始终更小，但当1≤1/Ro≤8时，上壁面附近温度脉动的极大值比同转速下θ_m=0.0时更大，其上壁面极大值位置约为上壁面与密度倒置界面平均高度h_θm的中间，即h(θ_{rms, max})≈1+h_θm/2。当1/Ro≤2时，在科氏力作用下密度倒置界面主要发生周向变形；当1/Ro>2时，涡柱结构的形成使得密度倒置界面发生较显著的垂向变形。因此，与θ_m=0.0的情况相比，θ_m=0.5时上壁面附近较强的温度脉动是由密度倒置界面的波动和变形引起的。图 6(c)、6(d)分别为径向速度和径向速度脉动的剖面。一方面，在θ_m=0.0和0.5时，速度脉动随转速的变化趋势与Horn和Shishkina^[5]关于OB和NOB流体的旋转对流结果一致；另一方面，在相同转速下，θ_m=0.5的速度和速度脉动始终小于θ_m=0.0的情况，与无旋转穿透对流的研究结果一致。

2.3.2 边界层

对于旋转热对流，我们利用上述斜率准则和耗散率方法分别定义热边界层和速度边界层厚度。如图 7所示，与旋转穿透对流情况相同，θ_m=0.5的热边界层厚度δ_{top/bottom, 0.5}^θ始终大于θ_m=0.0的热边界层厚度δ_{top/bottom, 0.0}^θ，并且随着1/Ro的增大具有不同的变化趋势。当θ_m=0.0时，上方热边界层厚度略大于下方热边界层厚度，当1/Ro增大时，δ_{top/bottom, 0.0}^θ先近似保持不变，随后突然增大；这种边界层厚度随转速的变化趋势也出现在OB和NOB流体中^[5]。当θ_m=0.5时，δ_{top/bottom, 0.5}^θ随转速的变化趋势与θ_m=0.0基本一致，但由于上方稳定区的存在，热边界层呈现明显的上下不对称性，即δ_{top, 0.5}^θ>δ_{bottom, 0.5}^θ。这与变物性(无密度倒置)流体旋转对流的结果(δ_top^θ < δ_bottom^θ)^[5]恰好相反。对于本文所研究的密度随温度非单调变化的NOB流体，始终存在δ_{top, 0.5}^θ>δ_{bottom, 0.5}^θ，并且(δ_{top, 0.5}^θ-δ_{bottom, 0.5}^θ)随转速提高而逐渐增加。在快速旋转时，δ^θ与1/Ro满足标度关系，其标度指数约为0.5。

随转速增加，上下速度边界层厚度先增大后减小，其增大阶段表明系统在为Ekman层的形成持续积累能量，而减小阶段表明Ekman层的激发。此时存在δ^u~1/Ro^-1/2的标度律关系，与OB流体的标度关系一致^[27]。穿透对流对流动上下对称性的破坏在速度边界层中也有体现，在θ_m=0.0时上下边界层厚度基本一致，而θ_m=0.5时的上下边界层厚度差|δ_{top, 0.5}^u-δ_{bottom, 0.5}^u|随转速先增大后减小。

2.4 系统整体传热

最后，讨论旋转穿透对流的归一化努塞尔数Nu/Nu₀随1/Ro的变化，其中Nu₀是无旋转情况下的Nu。如图 8所示，按照Nu/Nu₀的变化趋势，将1/Ro划分成3个区域，分别用罗马数字Ⅰ、Ⅱ和Ⅲ表示。

在区域Ⅰ(对应于低转速)，不同转速下的Nu/Nu₀几乎相同，表明传热由浮力主导。对于中等转速，即区域Ⅱ，Nu/Nu₀突然增大，达到最大值后，向区域Ⅲ过渡，热传输迅速下降。通常认为，区域Ⅱ中的传热增强归功于柱状涡结构的形成。柱状涡从边界层吸收额外的热量，这一过程被称为Ekman泵送。对于本文考虑的θ_m=0.0的冷水，Ekman泵送机制能够有效增强热传输，Nu的最大值约为无旋转时的118.5 %；但对于θ_m=0.5的情况，在1/Ro=2.0时达到最佳传热效果，此时Nu约为无旋转时的103.4 %。在本文所研究的参数范围内，旋转对穿透RBC传热的强化效果比较有限，这并不意味着Ekman泵送机制失效，而是由于仅在热边界附近能够产生较强的Ekman涡，该机制在冷边界的作用效果因密度倒置效应而明显减弱。在区域Ⅲ，随1/Ro的增大，传热性能骤降，表明Taylor-Proudman效应^[28-29]的增强削弱了Ekman泵送机制，对应于图 5所示高转速下涡柱的崩溃。

3 结论

本文采用DNS研究θ_m=0.0, 0.5，Ra∈[1×10⁴, 1×10⁸]范围内的冷水无旋转穿透对流，以及相同θ_m值，Ra=1×10⁷，1/Ro∈[0,1,2,4,6,8,10]范围内的旋转穿透对流。重点关注参数空间内的流动特征、边界层和传热特性，分析旋转和密度倒置对冷水热对流的影响。主要结论如下：

1) 在无旋转情况下，随着Ra的增大，流动由热传导过渡到混沌态，最终发展为湍流。在低Ra时，流动为稳态流动，一个规则且光滑的对流胞充满中心区域，θ_m=0.5与θ_m=0.0的流动具有显著差异；在中高Ra时，流动呈现非稳态特征且具有更多的小尺度结构，密度倒置效应对流动的影响减弱。此外，θ_m=0.5的穿透效应使得流态转变相比θ_m=0.0滞后，但文献中用于解释该滞后性的有效瑞利数Ra_E不能完全解释穿透对流中的动力学特性。我们认为稳定层是影响整体动力学特性的主要因素，其影响效果值得进一步探究和量化。Nu和热边界层厚度均与Ra满足标度关系，且该标度关系几乎不随θ_m发生变化，在本文中，标度关系为Nu~Ra^0.307，δ^θ~Ra^-0.315。

2) 在旋转穿透对流中，随转速增加，流态由LSC和羽流结构首先转变为CTC结构，然后CTC结构在高转速下崩溃。与无旋转时相比，中心温度在中低速旋转时保持不变，在高速旋转时有所升高。旋转的作用效果包括：加剧了θ_m=0.5时的上下非对称性，上下热边界层厚度差异随转速升高而增大，上下速度边界层厚度差异随转速升高先增大后减小。对于θ_m=0.0的无穿透工况，1/Ro=8.0的热传输性能最强，Nu高达无旋转时的118.5 %。对于θ_m=0.5的有穿透工况，当1/Ro=2.0时达到最佳传热效果，此时Nu约为无旋转时的103.4 %。