辐射下加热腔内自然对流传热特性研究

Characterization of Natural Convection Heat Transfer in a Heated Cavity under Radiation

摘要: 辐射吸收引发的自然对流现象广泛存在于自然界以及工业领域，为了探究辐射引发的自然对流传热特性，采用理论分析与数值模拟相结合，对由下方入射辐射加热二维方腔内半透明流体所引发的自然对流进行研究。结果表明：同一瑞利数下，流体失稳时间随光学厚度τ 的增加而减小，当τ=20 时，流体失稳时间最短；通过尺度分析推导出当光学厚度τ → 0 时，热边界层厚度与瑞利数及光学厚度之间满足标度率： ，并将数值模拟结果与之对比， Ra_I数的标度指数基本吻合；同时发现，不同光学厚度下， Nu数随着Ra_I 数的增加而逐渐增大，随着光学厚度的增加， 中的标度指数β 反而减小；最后采用数值模拟方法，得到了Nu 数与光学厚度及Ra_I 数之间的标度关系： ，并利用尺度分析方法推导出光学厚度τ → 0 时， Nu数与光学厚度τ 满足： 。

Abstract: The phenomenon of natural convection induced by radiation absorption exists widely in nature as well as in industry. In order to investigate the heat transfer characteristics of radiation-induced natural convection, a combination of theoretical analyses and numerical simulations is used to in-vestigate the natural convection induced by incident radiation from below heating a two-di- men-sional square cavity with a semi-transparent fluid. The results show that at a constant Rayleigh number, the fluid destabilization time decreases with the increase of the optical thickness, when τ=20, with a large fluctuation of the temperature at the lower wall surface; through the scale analysis, it is deduced that when the optical thickness τ → 0 , the scale rate is satisfied between the temperature boundary layer thickness and the Rayleigh number and the optical thickness: , and the numerical simulation results are compared with it, and the scale exponent of the Rayleigh number Ra_I is basically coincident. It is also found that the Nu num-ber gradually increases with increasing Ra_I number for different optical thicknesses, and the scaling index β in decreases instead with increasing optical thickness. By using the numerical simulation method, the scaling relationship between the Nu  number and the optical thickness and Ra_I number is obtained: , and the scaling analysis method is used to deduce that the Nu number and the optical thickness τ satisfy: when the optical thickness τ → 0.

1. 引言

自然对流是自然界和各类工程问题中的常见物理现象，与人们的日常生活息息相关。其不仅存在于冰川湖泊的融化 [1] 、大气及海洋环流 [2] [3] 、地幔对流板块运动 [4] 等自然现象中；还广泛应用于核反应堆冷却 [5] 、芯片散热 [6] 、磁流体传热和新型储能材料 [7] [8] 等工业问题中。而Rayleigh-Bénard (RB)对流系统 [9] ，则是从众多复杂的自然对流现象以及工程问题中简化出来的研究对流问题的经典模型。由于RB系统中存在多种热量传递方式，因此它也是传热领域经典的课题，研究RB系统对实际情况下各类对流过程中的传热机理有重要的意义。

由于经典的RB对流模型其上下边界温度恒定，这与一些实际系统的热边界条件并不相符。如太阳辐射在地面上反射对大气进行加热现象 [10] 、水冷壁受辐射热升温 [11] 等由辐射内热源驱动的热对流问题。在辐射驱动的热对流问题中，众多学者对其开展了模拟、理论和实验等方面研究。在理论分析方面，Lei等 [12] 对三角封闭腔进行标度分析，根据瑞利数Ra、底部三角倾斜角 $\tan \alpha$ 和临界瑞利数 $R{a}_c$ 对流动状态进行了描述。Mao等 [13] 对小斜坡水库吸收辐射引发的自然对流进行了尺度分析研究，获得瑞利数相对于水平位置的两个临界函数，并根据底部坡度和最大水深对流体流动区域进行划分，建立位置相关标度对不同区域的流量进行说明。Hattori等 [14] 通过线性理论与直接稳定性分析，对辐射驱动的自然对流进行了探究，并结合水平通流和横向约束，研究了两种不同形式的流动不稳定性，并对每种不稳定性模式形成的本征问题进行数值求解，获得了本征函数与失稳条件。在模拟方面：Miquel等 [15] 研究了努塞尔数Nu与普朗特数Pr之间的依赖关系，对于无应力下表面，传热标度率满足 $Nu\sim P{r}^{1/2}$ ，而对于无滑移下表面，可能呈现两种传热标度率， $Nu\sim P{r}^{1/2}$ 与 $Nu\sim P{r}^{1/6}$ 。Hattori等 [16] 采用理论分析与模拟验证相结合的方法，对上方直接入射辐射所形成的稳定热分层以及由底部边界热通量所产生的羽流进行了相关研究，并将标度关系与数值模拟进行了对比。在实验方面：Bouillaut等 [17] 搭建了辐射对流实验装置，通过改变光学厚度，对腔体内流体流动机理进行研究，结果表明：随着光学厚度 $\tau$ 的变化， $Nu~R{a}_T^{\gamma }$ 的指数 $\gamma$ 可在1/3至1/2范围内变动，该标度指数范围与RB系统所得到的范围相一致，标度指数 $\gamma$ 由1/3转变为1/2的临界点为 $R{a}_T~P{r}^{-3}{\tau }^6$ 与 $Nu~P{r}^{-1}{\tau }^2$ 。但仍有部分学者认为该标度区间并不适用于含有内热源的物理问题。Klein等 [18] 通过ODT模型对文献 [17] 进行更大范围的数值模拟，发现其标度指数上界限为0.55，并非是Bouillaut等所得到的标度指数。

以上学者均对辐射所引发的自然对流进行了研究。但在其研究过程中，未能充分考虑光学厚度对传热特性的影响。因此，本文将采用理论分析与数值模拟相结合的方式，探讨不同光学厚度、不同瑞利数 $R{a}_I$ 下辐射引发的自然对流传热特性。接下来第4.1节主要探讨光学厚度对流体失稳时间的影响，采用模拟方法对失稳时间进行探究。第4.2节简要分析光学厚度对温度边界层厚度的影响。第4.3节讨论努塞尔数、光学厚度以及瑞利数三者之间的标度指数关系。最后就辐射对自然对流的影响进行总结。

2. 问题描述及控制方程

2.1. 物理模型

考虑如图1所示几何模型。方腔左、右边界为周期边界，下表面为透明、无滑移绝热表面，上表面为透明、无应力绝热表面，内部为半透明流体。辐射从下方垂直均匀射入方腔内，流体吸收辐射后温度升高，形成密度差，近而在浮力的作用下发生自然对流。

Figure 1. Geometry model

图1. 几何模型

2.2. 数学模型

对于二维方腔区域，辐射从下方垂直均匀射入，辐射能随着流体吸收的深度逐渐衰弱，满足比尔–朗伯定律：

$Q=I_0\exp \left(-\eta y\right),y>0$ (1)

其中， $\eta$ 为流体的辐射吸收系数， $I_0$ 为入射辐射强度。下方垂直辐射加热效果等效于一个内热源的作用，其引发的自然流动现象可由Navier-Stokes方程组和能量方程描述。其中，浮升力的影响采用Boussinesq近似描述：

$\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}=0$ (2)

$\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial v}{\partial y}=-\frac{1}{\rho }\frac{\partial p}{\partial x}+\nu \left(\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial y^2}\right)$ (3)

$\frac{\partial v}{\partial t}+u\frac{\partial v}{\partial x}+v\frac{\partial v}{\partial y}=-\frac{1}{\rho }\frac{\partial p}{\partial y}+\nu \left(\frac{{\partial }^{2}v}{\partial x^2}+\frac{{\partial }^{2}v}{\partial y^2}\right)+\beta g\left(T-\bar{T}\right)$ (4)

$\frac{\partial T}{\partial t}+u\frac{\partial T}{\partial x}+v\frac{\partial T}{\partial y}=a\left(\frac{{\partial }^{2}T}{\partial x^2}+\frac{{\partial }^{2}T}{\partial y^2}\right)+\frac{I_0}{\rho c_p}\eta {\text{e}}^{-\eta y}$ (5)

其中， $\rho$ 为流体密度， $\nu$ 为流体运动粘度， $\beta$ 为热膨胀系数， $c_p$ 为流体比热容， $\bar{T}$ 为流体平均温度， $p$ 为流体压力， $a$ 为流体热扩散系数。

对应边界条件为：

${u|}_{x=0}={u|}_{x=L},{v|}_{x=0}={v|}_{x=L},{T|}_{x=0}={T|}_{x=L}$ (6)

$u=v=0,\frac{\partial T}{\partial y}=0,y=0$ (7)

$\frac{\partial u}{\partial x}=\frac{\partial v}{\partial y}=0,\frac{\partial T}{\partial y}=0,y=H$ (8)

上述系统中，流动在一定时间后可进入准稳态阶段，而系统平均温度由于壁面绝热随时间线性增加：

$\bar{T}=T_0+\frac{{\displaystyle {\int }_0^L{\displaystyle {\int }_0^H\frac{I_0}{\rho c_p}\eta {\text{e}}^{-\eta y}\text{d}y\text{d}x}}}{HL}t=T_0+\frac{I_0}{\rho c_{p}H}\left(1-{\text{e}}^{-\eta H}\right)t$ (9)

其中： $T_0$ 为系统初始温度。

采用以下标尺对控制方程进行无量纲化：

$x^{*}=\frac{x}{H},y^{*}=\frac{y}{H},t^{\ast }=\frac{ta}{H^2},v^{*}=\frac{vH}{a}，u^{*}=\frac{uH}{a},{\theta }^{*}=\frac{\left(T-\bar{T}\right)\rho c_{p}a}{I_{0}H},p^{\ast }=\frac{pa{c}_p}{H^{2}g\beta I_0}$ (10)

以及无量纲控制参数，瑞利数 $R{a}_I$ 、普朗特数 $Pr$ 以及光学厚度 $\tau$

$R{a}_I\text{=}\frac{g\beta H^4{I}_0}{\upsilon \rho c_p{a}^2},Pr=\frac{\nu }{a},\tau =H\eta$ (11)

代入原方程，整理后略去变量上标星号，可得无量纲控制方程组为：

$\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}=0$ (12)

$\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}=-R{a}_{I}Pr\frac{\partial p}{\partial x}+Pr\left(\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial y^2}\right)$ (13)

$\frac{\partial v}{\partial t}+u\frac{\partial v}{\partial x}+v\frac{\partial v}{\partial y}=-R{a}_{I}Pr\frac{\partial p}{\partial y}+Pr\left(\frac{{\partial }^{2}v}{\partial x^2}+\frac{{\partial }^{2}v}{\partial y^2}\right)+R{a}_{I}Pr\theta$ (14)

$\frac{\partial \theta }{\partial t}+u\frac{\partial \theta }{\partial x}+v\frac{\partial \theta }{\partial y}=\frac{{\partial }^2\theta }{\partial x^2}+\frac{{\partial }^2\theta }{\partial y^2}+\tau {\text{e}}^{-\tau y}-\left(1-{\text{e}}^{-\tau }\right)$ (15)

记 $\Gamma =\frac{L}{H}$ ，则相应无量纲边界条件如下：

${u|}_{x=0}={u|}_{x=\Gamma },{v|}_{x=0}={v|}_{x=\Gamma },{\theta |}_{x=0}={\theta |}_{x=\Gamma }$ (16)

$u=v=0,\frac{\partial \theta }{\partial y}=0,y=0$ (17)

$\frac{\partial u}{\partial x}=0,v=0,\frac{\partial \theta }{\partial y}=0,y=1$ (18)

易知对于任意时刻均有 ${\displaystyle {\int }_0^1{\displaystyle {\int }_0^{\Gamma }\theta \text{d}x\text{d}y}}=0$ 。

3. 数值求解方法与无关性验证

3.1. 数值求解方法

采用有限体积法对控制方程(12)~(15)进行离散，基于SIMPLE算法处理压力与速度耦合，扩散项采用中心差分格式离散，对流项采用QUICK格式离散，时间项应用二阶隐式迎风的离散方式。采用的无量纲控制参数见表1。

Table 1. Dimensionless parameter

表1. 无量纲参数

3.2. 网格无关性与时间独立性验证

分别采用4种网格尺寸400 × 130、640 × 185、880 × 246、1120 × 310进行网格无关性验证，选取瑞利数 $R{a}_I={10}^{10}$ ，通过时间空间平均的努塞尔数Nu进行网格无关性与时间独立性检验。详细信息如表2所示。

$Nu=\frac{1}{\Delta \theta }$ (19)

其中， $\Delta \theta ={\langle \bar{\theta }\rangle |}_{y=0}-{\langle \bar{\theta }\rangle |}_{y=\frac{1}{2}}$ 为下壁面与中心位置的时均温差， $\langle \theta \rangle$ 表示空间平均温度， $\bar{\theta }$ 表示时间平均温差。x方向网格均匀布置，y方向网格在边界层处进行加密，确保边界层内部至少存在15个网格。流体光学厚度 $\tau$ 取为无穷大，即此时流体不透明，辐射加热效果等效于下壁面恒热流加热。

Table 2. Grid information

表2. 网格信息

(a) (b)

Figure 2. (a) Comparison plot of Nu number for different grid sizes (b) Comparison plot of Nu number for different time steps

图2. (a) 不同网格尺寸下的Nu数对比图 (b) 不同时间步长下的Nu数对比图

图2(a)为四组网格系统在准稳态阶段Nu数随时间的发展。其中，通过计算可以得到，网格1、网格2以及网格3与网格4的Nu误差相差分别为3.23%、0.98%、1.03%，因此为了保证精度以及节省计算成本，本文所有工况选取网格2进行模拟计算。图2(b)为三组时间步长在准稳态阶段Nu数随时间的发展。时间步长为5.75 × 10⁻⁶，与时间步长为1.44 × 10⁻⁶相比，Nu误差相差5.83%，时间步长为2.87 × 10⁻⁶与时间步长为1.44 × 10⁻⁶相比，Nu误差相差0.43%，本文选取时间步长为2.87 × 10⁻⁶进行计算可以更好的节约时间成本。

3.3. 程序有效性验证

文献 [17] 采用实验的方式，测得多组工况下时间空间平均努塞尔数Nu与瑞利数 $R{a}_T$ 的标度关系，本文对照文献 [17] 中标度关系对程序进行有效验证。选取其中 $\tau =20、\mathrm{Pr}=7$ 时的工况，文献中努塞尔数Nu与瑞利数 $R{a}_T$ 的标度指数关系近似满足 $Nu\sim R{a}_T^{0.5}$ ，图3给出了本文努塞尔数Nu与瑞利数 $R{a}_T$ 的标度指数与文献中的标度指数结果对比图。对比结果表明，当前研究所用到的程序准确可靠。

Figure 3. Scaling relationship between Nusselt number and Rayleigh number (τ = 20)

图3. (τ = 20)时努塞尔数Nu与瑞利数Ra_T的标度关系

4. 计算结果与分析

4.1. 光学厚度对流动失稳时间影响

流体流动可分为，初始阶段、过渡阶段、准稳态阶段三个部分。在初始阶段，流体速度较低，此时能量主要通过导热方式传递。由于腔体下方辐射的作用，下壁面温度随着时间逐渐升高，当温度边界层温差达到一定值时，下表面将产生向上的热羽流，此时流体热边界层被破坏，流体将从初始阶段转为过渡阶段。

流体失稳需要满足热边界层温差超过一定范围，即当地瑞利数 $R{a}_L$ 超过临界值 $R{a}_{Tc}$ 。

$R{a}_L=R{a}_I\Delta {\theta }_b{\delta }^3\ge R{a}_{Tc}$ (20)

其中， $\Delta {\theta }_b$ 为热边界层温差， $\delta$ 为热边界层厚度。对于小光学厚度引发的对流，失稳临界瑞利数影响较小，并且由公式(15)可得：

$\delta \sim \sqrt{t}$ (21)

由公式(20)与(21)得：

$R{a}_L=t\ge {\left(\frac{R{a}_{Tc}}{R{a}_I\Delta {\theta }_b}\right)}^{\frac{3}{2}}$ (22)

(a) (b) (c) (d)

Figure 4. Change of critical time under different Rayleigh number

图4. 不同瑞利数下临界时间变化

图4给出了不同瑞利数、光学厚度下，腔体下表面流体平均温度随时间的变化情况。可以看出，在同一瑞利数下，随着光学厚度的增加，临界失稳时间逐渐减小。其原因在于，当光学厚度较大时，大部分辐射能被边界层吸收，导致边界层温度升高较快，并以热羽流形式带走边界层能量，导致流体达到失稳状态。而当光学厚度较小时，大部分辐射能将透过边界层对腔体内部流体进行加热，此时边界层内部温升较慢，流体达到失稳时间较长。同时，该结论可由公式(22)得到验证。

表3给出了不同工况下的具体失稳时间。可以看出，瑞利数 $R{a}_I={10}^7$ 时的失稳时间相比于其他工况更长。这是由于低瑞利数下流体的粘性力较强，流动强度较弱，因此流动失稳时间较长。

Table 3. Change of critical instability time under different Rayleigh number and optical thickness

表3. 不同瑞利数、光学厚度下失稳临界时间变化

4.2. 光学厚度对热边界层厚度的影响

早期Chillá等 [19] 人通过实验的方式得到 ${10}^6\le R{a}_T\le {10}^8$ ， $Pr=4$ 范围内热流边界层的厚度，得出结论，温度边界层内部的温度基本满足理论的线性变化。因此在提取边界层方法中，通常采用斜率法来取得温度边界层的厚度，温度场取时间空间平均，在下壁面处可以得到近似直线的温度分布。同时，在腔体中心 $y=H/2$ 处与下壁面温度边界层之间，温度分布同样可以近似的看成线性变化，两条切线的交点到下壁面的距离可被看作温度边界层厚度。

图5给出了不同光学厚度 $\tau$ 下，瑞利数 $R{a}_I$ 与热边界层厚度 $\delta$ 之间的关系，可以看出随着瑞利数 $R{a}_I$ 的增加，边界层厚度逐渐减小。结论式(23)与分析推导结果式(27)可以验证这一结论。与此同时，发现不同光学厚度 $\tau$ 下，瑞利数 $R{a}_I$ 的标度指数几乎一致。为了验证其标度率的正确性，下文将采用理论分析方法对其进行验证。为方便模拟与理论分析结果对比，模拟结果将取平均标度指数进行标度探究，如公式(23)所示：

$\delta \sim R{a}_I^{-0.216}$ (23)

Figure 5. Relation between Rayleigh number and thickness scaling rate of thermal Boundary layer

图5. 瑞利数Ra_I与热边界层厚度δ标度率关系

流体在准稳态阶段，由动量方程(14)，热流体脱离热边界层时惯性力与粘性力及浮力之间达到平衡：

$\frac{v^2}{\delta }\sim -Pr\frac{v}{{\delta }^2}+R{a}_{I}Pr\Delta {\theta }_b$ (24)

其中： $\Delta {\theta }_b$ 为热边界层温差。

又由能量方程(15)，对流项、扩散项以及源项三者之间达到平衡：

$v\frac{\Delta {\theta }_b}{\delta }\sim \frac{\Delta {\theta }_b}{{\delta }^2}\sim \tau {\text{e}}^{-\tau \delta }-\left(1-{\text{e}}^{-\tau }\right)$ (25)

若 $\tau \delta \to 0$ 成立，则式(25)可进一步简化：

$v\frac{\Delta {\theta }_b}{\delta }\sim \frac{\Delta {\theta }_b}{{\delta }^2}\sim \tau$ (26)

由式(24)、(26)可得：

$\delta \sim {\tau }^{-0.2}R{a}_I{}^{-0.2}{\left(1+P{r}^{-1}\right)}^{0.2}$ (27)

由模拟方法得到的结论式(23)与尺度分析推导的结果式(27)相比可得，标度指数量级基本相当，进一步说明数值模拟结果的准确性，同时也展现了瑞利数 $R{a}_I$ 对流体流动有重要影响。

4.3. 光学厚度对标度指数的影响

在湍流流动中，热量如何运输是研究热对流的重点关注问题，即努塞尔数与瑞利数之间的变化规律。对于有辐射驱动的对流系统，热量运输也与光学厚度有紧密的关系。因此，本节将主要探讨光学厚度对努塞尔数的影响，试图找到小光学厚度下 $Nu\sim {\tau }^{\alpha }R{a}_I^{\beta }$ 的表达式。图6(a)中给出了由模拟得到不同光学厚度下，Nu数随 $R{a}_I$ 数的变化关系，随着 $R{a}_I$ 数的增加，Nu数逐渐增大。随着光学厚度的增加， $Nu\sim R{a}_I^{\beta }$ 中的标度指数 $\beta$ 反而减小。在 $0.5\le \tau \le 20$ 范围，有 $0.302\le \beta \le 0.341$ 。从图6(b)可看出，以 $\tau =20$ 时结果为基准做归一化，可以得到关系式：

$Nu\sim {\tau }^{-0.918}N{u}_{20}\sim {\tau }^{-0.918}R{a}_I{}^{0.302}$ (28)

(a) (b)

Figure 6. (a) Relationship between Rayleigh number and Nusselt number under different optical thicknesses (b) Relationship between Nusselt number and optical thickness

图6. (a) 不同光学厚度下瑞利数Ra_I与努塞尔数Nu的变化关系 (b) 努塞尔数Nu与光学厚度τ之间关系

实际上根据式(15)，可对努塞尔数Nu与光学厚度 $\tau$ 之间的关系 $Nu\sim {\tau }^{\alpha }$ 做进一步理论推导，由式(15)可得：

$v\frac{2\Delta \theta }{H}\sim \frac{4\Delta \theta }{H^2}\sim \tau {\text{e}}^{-\frac{\tau }{2}}-\left(1-{\text{e}}^{-\tau }\right)$ (29)

当 $\tau \to 0$ 时，式(29)进一步可以得到：

$\Delta \theta \sim \frac{\tau }{4}$ (30)

根据式(19)，有：

$Nu\sim 4{\tau }^{-1}$ (31)

由理论推导结果式(31)与模拟结果式(28)对比可以看出，光学厚度 $\tau$ 的指数 $\alpha$ 几乎相当，也进一步说明模拟结果的准确性。

5. 结果与讨论

本文对下方垂直入射辐射引发的自然对流问题进行了相关研究，主要探讨不同瑞利数 $R{a}_I$ ( ${10}^7\le R{a}_I\le {10}^{10}$ )与不同光学厚度 $\tau$ ( $0.5\le \tau \le 20$ )对流体传热特性的影响，得到如下结论：

1) 流体在初始阶段，恒定瑞利数下，流体失稳时间随着光学厚度的增加而减小，当 $\tau =20$ 时，下壁面温度波动最大，失稳时间最小。当瑞利数 $R{a}_I={10}^7$ 时，流体失稳时间相比与其他工况下时间更长。

2) 通过尺度分析方法，推导出了光学厚度 $\tau \to 0$ 时温度边界层厚度 $\delta$ 与瑞利数 $R{a}_I$ 之间的标度关系，两者之间满足： $\delta \sim {\tau }^{-0.2}R{a}_I^{-0.2}{\left(1+P{r}^{-1}\right)}^{0.2}$ 。同时将模拟结果( $\delta \sim R{a}_I^{-0.216}$ )与尺度分析结果进行对比，结果

吻合较好，验证了数值模拟的准确性。

3) 通过数值模拟方法得到了光学厚度 $\tau \to 0$ 时，努塞尔数Nu与瑞利数 $R{a}_I$ 以及光学厚度 $\tau$ 三者之间标度关系： $Nu\sim {\tau }^{-0.918}R{a}_I^{0.302}$ 。并采用尺度分析方法推导出努塞尔数Nu与光学厚度 $\tau$ 之间满足： $Nu\sim 4{\tau }^{-1}$ 。两种结果中光学厚度的标度指数几乎相当。同时，由模拟得到不同光学厚度下，随着Nu数的增加， $R{a}_I$ 数逐渐增大。随着光学厚度的增加， $Nu\sim R{a}_I^{\beta }$ 中的标度指数 $\beta$ 反而减小。

基金项目

国家自然科学基金青年基金(52206103)、上海市青年科技英才扬帆计划(21YF1430400)。

NOTES

^*通讯作者。