非平直下壁面对辐射吸收引发的对流影响研究

Study on the Influence of Non-Uniform Lower Wall on Convection Induced by Radiation Absorption

摘要: 对非平直下壁面辐射驱动的自然对流进行了数值模拟，基于ANSYS FLUENT求解控制方程，分析了几何参数和Rayleigh数对辐射驱动自然对流流动和传热的影响。结果表明，随着Rayleigh数的增加，流函数的对称结构被破坏，上方流线图变得近似水平，小波长的工况会出现水平流动；几何参数对羽流穿透深度和下壁面温度分布有重要影响，随着振幅的增大，系统的温度均匀性变小，当周期个数为1时，羽流穿透深度最大；当波长长度为1.8时，下壁面温度分布出现过渡情况，由波峰温度最高变为波谷温度最高。

Abstract: A numerical simulation of radiation-driven convection on a non-flat lower wall is carried out. The influence of geometric parameters and Rayleigh number on radiation-driven natural convection flow and heat transfer is analyzed by solving the governing equation based on ANSYS FLUENT. The results show that with the increase of Rayleigh number, the symmetric structure of the stream function is destroyed, the upper flow plot becomes nearly horizontal, and horizontal flow occurs in the condition of small wavelength. The geometrical parameters have an important effect on the plume penetration depth and the temperature distribution on the lower wall. With the increase of the amplitude, the temperature uniformity of the system becomes smaller. When the number of periods is 1, the plume penetration depth reaches the maximum. When the wavelength length is 1.8, the temperature distribution of the lower wall appears a transition, from the highest peak temperature to the highest trough temperature.

1. 引言

自然对流作为流体力学领域的经典问题，广泛存在于自然界和工程中[1]，如湖泊、海洋、近岸水域、太阳能集热器中[2] [3]。RB对流是一个高度理想化的系统[4]，可以在控制良好的实验室中实验创建，但在自然界中很少遇到。实际系统中的边界板可能具有几何和加热边界的不规则性，并且它们对系统响应产生重大的影响[5]。

对流的性质取决于流动系统的几何形状和外部施加的热边界条件，模型几何形状的不均匀性会影响起始条件[6]。不均匀性的形式数不胜数，前人选择了以正弦波纹的简单几何形状作为参考，这种系统在自然界中经常出现，例如农村环境中的山丘和山谷以及微流体装置的几何图案[7] [8]。Abtahi等人[9]分析了两块波纹等温板形成的水平槽内的自然对流，研究发现波纹的存在影响传导热流并产生对流热流，并且无论加热强度如何，自然对流都会发生，对流热流只与波数有关，与瑞利数的三次方成正比。边界温度的不均匀性会影响RB对流的开始，从而导致各种形式的对称性破坏[10]。加热不均匀性的存在导致了系统响应的定性不同，不均匀加热产生水平密度变化，导致水平和垂直温度梯度，产生水平对流的运动，称为热漂移[11]。周期性加热和均匀加热的组合通常会导致二次对流开始的相关临界瑞利数的减少，研究表明，通过利用空间加热不均匀性可以增加热流，这在定性上改变了热传输[12]。Floryan [13]等人探讨了空间非均匀分布加热模式下两个水平板的相对运动，水平板的相对运动需要外力来维持，研究结果表明，足够强度的空间加热可以完全消除对驱动力的需要，减阻效果仅在小雷诺数流动中有效，并且在具有较小普朗特数的流体中更强。Inasawa等人[14]采用理论分析和实验方法分析了空间分布加热通道内的由压力驱动的二维稳定流动，通过改变整体流动分布，实现了减小层流剪切阻力的目的，研究结果表明，平均压力梯度的加入导致浮力驱动流和压力梯度驱动流之间的竞争。在许多工程应用中，有效地提高对流换热是至关重要的，壁面粗糙度的引入被认为是一种有效的手段。研究表明[15]，粗糙度并不总是可以增强传热，存在一个临界的粗糙元高度，只有当高度高于临界的粗糙元高度时，才可以增强传热。

上述对于辐射驱动对流的研究通常只考虑均匀边界，或者是非均匀边界的RB对流。由辐射吸收引发的自然对流更是吸引许多学者对自然对流现象的研究[16] [17]。在自然界和许多工业过程中，不规则表面的辐射驱动自然对流是经常遇到的情况，例如城市大气边界层中的城市表面通常不均匀、在深海中的海床和洋底也是具有非常不规则的地形结构；另一方面，在研究地幔对流时，大陆板块与海洋板块的传热效率也是存在差异的，此时需要考虑非均匀热边界条件下的辐射驱动对流[18]。目前，非均匀边界条件下的关于辐射驱动自然对流的研究极少。类似于研究非均匀的RB对流中采用正弦非均匀壁面结构，本文拟从采用下方余弦形状的非平直下壁面系统，假设入射辐射从系统上方照射，开展数值模拟研究，研究不同工况下的非平直下壁面边界的辐射驱动对流问题，揭示非均匀边界条件的辐射驱动自然对流规律，为工业应用中的辐射驱动自然对流现象提供初步的理论指导。

2. 物理模型和数学模型

2.1. 物理模型

图1示出本文所研究的非平直下壁面对辐射驱动对流影响研究的物理模型。模型的上方是平直壁面，下方是具有余弦函数分布函数形状的非平直结构，左右为周期性边界条件，上表面为自由液面，下表面为固定壁面，速度满足无滑移边界条件。到达下壁面的辐射被完全吸收，然后作为热流边界条件释放。

Figure 1. Physical model

图1. 物理模型

假设系统的几何形状为：

$y_L\left(x\right)=-h+A\cos \left(\frac{2\pi N}{L}x\right)$ (1)

$y_U\left(x\right)=0$ (2)

式中，h为高度，A为余弦函数的振幅，L为波长，N为周期个数。

2.2. 数学模型

根据比尔–朗伯定律，入射辐射随深度指数衰减，在不同深度的流体层接收到的辐射强度为：

$I=I_0\exp \left(\eta y\right)\left(y\le 0\right)$ (3)

式中，I₀为入射辐射强度，单位为W∙m⁻²，$\eta$ 为辐射衰减系数，单位为m⁻¹。

引入无量纲参数：

$\begin{array}{l}X=x\eta ,Y=y\eta ,U=\frac{u}{k\eta },V=\frac{v}{k\eta },P=\frac{p{c}_P\kappa {\eta }^2}{g\beta I_0},\tau =t\kappa {\eta }^2,\\ \theta =\frac{\left(T-T_{ref}\right)\rho c_p\kappa \eta }{I_0},Ra=\frac{g\beta I_0}{\rho c_P\nu {\kappa }^2{\eta }^4},Pr=\frac{\upsilon }{\kappa }\end{array}$ (4)

其中：x、y、u、v、p、t、T分别为有量纲的x方向长度、y方向长度、x方向的速度、y方向的速度、压力、时间、温度，X、Y、U、V、P、$\tau$ 、$\theta$ 为其对应的无量纲量。Ra为瑞利数，Pr为普朗特数；$\kappa$ 为流体热扩散系数，$c_P$ 为比热容，$\beta$ 为热膨胀系数，$T_{ref}$ 为流体的参考温度，$\upsilon$ 为运动粘度，g为重力加速度。

在Boussinesq假设下，描述辐射吸收引发的自然对流问题的无量纲控制方程为：

$\frac{\partial U}{\partial X}+\frac{\partial V}{\partial Y}=0$ (5)

$\frac{\partial U}{\partial \tau }+U\frac{\partial U}{\partial X}+V\frac{\partial U}{\partial Y}=-RaPr\frac{\partial P}{\partial X}+Pr{\nabla }^{2}U$ (6)

$\frac{\partial V}{\partial \tau }+U\frac{\partial V}{\partial X}+V\frac{\partial V}{\partial Y}=-RaPr\frac{\partial P}{\partial Y}+Pr{\nabla }^{2}V+RaPr\theta$ (7)

$\frac{\partial \theta }{\partial \tau }+U\frac{\partial \theta }{\partial X}+V\frac{\partial \theta }{\partial Y}={\nabla }^2\theta +\exp \left(Y\right)-\frac{1}{\eta h}$ (8)

无量纲初始条件及边界条件：

$\tau =0,U=V=0,\theta =\text{0}$ (9)

${\left(U,V,\theta \right)|}_{X=0}={\left(U,V,\theta \right)|}_{X=L_x\eta }$ (10)

$Y=0,\frac{\partial U}{\partial Y}=0,V=0,\frac{\partial \theta }{\partial Y}=0$ (11)

$Y=-L_y\eta ,U=V=0,\frac{\partial \theta }{\partial Y}=-\frac{{\text{e}}^{-L_y\eta }}{\sqrt{1+{\left[\frac{2\pi MA}{L}\sin \left(\frac{2\pi M}{L}x\right)\right]}^2}}$ (12)

3. 数值求解方法和数学模型验证

3.1. 求解方法简介

本文采用数值模拟方法，通过用户自定义函数(UDF)实现辐射源项和边界条件的输入，基于ANSYS FLUENT对无量纲控制方程(5)~(8)及其定解条件(9)~(12)进行数值求解。采用有限体积法进行辐射吸收引发的自然对流过程仿真。在仿真的过程中，采用控制体积法离散控制方程，QUICK格式离散对流项，中心差分格式离散扩散项，二阶隐式迎风格式离散时间项，使用耦合算法同时求解控制方程。在求解过程中，每个时间步内水平和垂直方向速度的残差低于10⁻⁴，能量方程残差低于10⁻⁶，视为计算收敛。

3.2. 程序有效性验证

使用程序与文献[19]算例进行了对比验证。验证算例工况为L_y = 2、Ra = 10⁹、Pr = 7、上表面绝热，分别对比了流动初始阶段(t = 0.00143)时的垂直方向温度分布、流动准稳态阶段(0.0144 ≤ t ≤ 0.0286)时的垂直方向时空平均温度分布。对比结果如图2所示，数字模拟所得温度分布情况与文献结果吻合良好，表明此程序运行结果的可靠性。

3.3. 网格无关和时间独立性验证

针对图1所示的辐射垂直射入的物理模型，分别选取了三组网格和三组时间步长进行网格无关性验证。图3给出了下壁面温度标准差随时间的变化，首先是边界层的扩散增长，随后流动失稳形成羽流，最后到达准稳态阶段。在准稳态阶段，下壁面的温度标准差趋于平均值附近，逐渐稳定。选取了500 × 250、800 × 408、1000 × 454的三组不同密度的网格，时间步长分别为2 × 10⁻⁶、8 × 10^-6、3 × 10⁻⁵可以看到，不同网格密度和不同时间步长的下壁面温度标准差相差不大，在误差范围内，结果合理。为了节约计算资源，本文选取网格密度为500 × 250，时间步长为3 × 10⁻⁵的网格进行后续的模拟。

(a) (b)

Figure 2. Spatiotemporal mean temperature distribution. (a) Quasi-steady state stage (0.0144~0.0286); (b) Initial stage

图2. 时空平均温度分布。(a) 准稳态阶段(0.0144~0.0286)；(b) 初始阶段

(a) (b)

Figure 3. Change of standard deviation of lower wall temperature with time. (a) Different mesh density; (b) Different time steps

图3. 下壁面温度标准差随时间的变化。(a) 不同网格密度；(b) 不同时间步长

4. 计算结果及分析

4.1. 羽流穿透深度

本节根据文献提出的羽流穿透深度为1的观点[19]，比较了非平直下壁面不同振幅和不同周期个数的羽流穿透深度，研究周期个数和振幅对羽流穿透深度的影响。

图4给出高度为1.5，长度为8，振幅为0.5的，周期个数为1的不同时刻的温度云图，初始阶段，边界层逐渐增长，羽流从下壁面最高点的波峰处开始产生，这是由于波峰离上表面更近，吸收到的辐射更多，温度更高，于是羽流先从此产生，类似于烟囱的抽吸效应。抽吸效应是指烟囱内外气体温度不同而引起气体密度差异，即烟囱抽力，它克服阻力推动烟气流动，在烟囱底部产生负压，形成抽吸效果。随后下壁面的其他地方也开始产生羽流，羽流逐渐上升，产生的羽流逐渐涌到上表面开始下沉，随着时间的发展，由于羽流的穿透高度为1，对于高度为1.5的模型，不能完全穿透，逐渐出现温度分层结构，上方是导热区，下方是混合区，到达准稳态阶段，分层明显。

(a) t = 0.01

(b) t = 0.015

(c) t = 0.00275

(d) t = 0.004

(e) t = 0.0085

(f) t = 0.03

Figure 4. Temperature fields of Ra = 10⁹, Pr = 0.7, h = 1.5 at different times

图4. Ra = 10⁹, Pr = 0.7, h = 1.5不同时刻的温度场

4.2. 系统温度均匀性

图5为高度0.5，相同振幅周期个数为1时不同波长时下壁面温度分布情况，图中可以看出当长度足够长时，由于波峰离上表面近，接收到的热流更多，因此波峰温度最高。随着长度逐渐减小，波谷区域变得更窄小，因此与波长长的情况相比，加热区域集中在波谷附近，因此接收的单位热流更多，存在过渡情况。通过模拟结果寻找过渡状态，可以看到长度L = 1.8时，波谷温度出现翻转现象，此时波峰波谷温度逐渐接近。随着长度进一步减小，由于下壁面波谷位置加热区域集中，波谷附近壁面温度最高。

Figure 5. Lower wall temperature distribution at different wavelengths

图5. 不同波长下壁面温度分布

图6给出高度1.5，周期个数均为1的不同振幅的系统温度标准差和温度最大偏差。可以看到，在流动达到失稳后，系统的温度标准差和最大偏差都会降低，这是因为下壁面热羽流产生，带走了大量热量，导致系统整体温度下降。当振幅为0时，达到失稳后，流动变得稳定，混合均匀，系统的温度标准差和最大偏差趋于稳定；当振幅增大，系统温度标准差和最大偏差随着时间发展而增大，这是由于振幅增大，导致下壁面最高温度与最低温度差别变大，混合变得越来越不均匀，导致系统的温度越来越不均匀，系统温度均匀性变差。

(a) (b)

Figure 6. Ra = 10⁹, Pr = 7, h = 1.5 with different amplitudes (a) Standard deviation of system temperature and (b) Maximum deviation of system temperature

图6. Ra = 10⁹, Pr = 7, h = 1.5不同振幅的(a) 系统温度标准差；(b) 系统温度最大偏差

4.3. 流函数分布

图7为高度0.5，振幅为0.06，长度0.5时不同瑞利数下的流函数分布，小瑞利数下流函数不对称，随着瑞利数增加，当瑞利数Ra = 10⁸时流函数出现混乱，出现类似水平的流动，随着之后瑞利数的进一步增加，上方出现更为水平的流线。

图8为高度0.5，振幅为0.06，长度1.8时不同瑞利数下的流函数分布，可以看出，瑞利数比较小的时候，流函数对称，形成对涡结构，随着瑞利数的增加，当瑞利数Ra = 10⁸时，流函数图的对称结构被破坏，上方流线从一边贯穿延伸到另外一边，随着瑞利数的进一步增加，上方流线越来越接近水平。

(a) Ra = 10⁵ (b) Ra = 10⁶

(c) Ra = 10⁷ (d) Ra = 10⁸

(e) Ra = 10⁹ (f) Ra = 10¹⁰

Figure 7. Stream functions for h = 0.5, A = 0.06, and L = 0.5 with different Rayleigh numbers

图7. h = 0.5, A = 0.06, L = 0.5时不同瑞利数下的流函数

(a) Ra = 10⁵

(b) Ra = 10⁶

(c) Ra = 10⁷

(d) Ra = 10⁸

(e) Ra = 10⁹

(f) Ra = 10¹⁰

Figure 8. Stream functions for h = 0.5, A = 0.06, and L = 1.8 with different Rayleigh numbers

图8. h = 0.5, A = 0.06, L = 1.8时不同瑞利数下的流函数

图9为高度0.5，振幅为0.06，长度3时不同瑞利数下的流函数分布，可以看出，瑞利数比较小的时候，流函数图对称，当瑞利数增大到10⁷时，流函数的对称结构开始被破坏，随后瑞利数进一步增大的时候，流函数的涡出现破裂，形成许多小涡，当瑞利数达到10⁹的时候，流线图从一边贯穿到另外一边，当瑞利数增加到10¹⁰时，可以观察到上方几乎水平的流动。

图10为高度0.5，振幅为0.06，长度6时不同瑞利数下的流函数分布，可以看出，小瑞利数下流函数对称，形成对涡结构，随着瑞利数增加，当瑞利数增大到10⁷的时候，流函数的对称结构被破坏，在两个大涡中间又形成许多小涡。当瑞利数达到10¹⁰的时候，流线图从一边贯穿到另外一边，没有观察到上方几乎水平的流动。

(a) Ra = 10⁵

(b) Ra = 10⁶

(c) Ra = 10⁷

(d) Ra = 10⁸

(e) Ra = 10⁹

(f) Ra = 10¹⁰

Figure 9. Stream functions for h = 0.5, A = 0.06, and L = 3 with different Rayleigh numbers

图9. h = 0.5, A = 0.06, L = 3时不同瑞利数下的流函数

(a) Ra = 10⁵

(b) Ra = 10⁶

(c) Ra = 10⁷

(d) Ra = 10⁸

(e) Ra = 10⁹

(f) Ra = 10¹⁰

Figure 10. Stream functions for h = 0.5, A = 0.06, and L = 6 with different Rayleigh numbers

图10. h = 0.5, A = 0.06, L = 6时不同瑞利数下的流函数

根据以上不同波长和不同瑞利数的流函数分布，可以看到波长小的更容易出现上方水平流动的情况，这是因为波长更小，相同振幅下的斜率就会更大，会有更多的水平速度分量推动流体向水平方向流动。小瑞利数时由于几乎没有流动或者流动比较弱，流函数对称。随着瑞利数增大，流函数图的对称结构逐渐被破坏，上方流线图越来越会往上变得近似水平。当波长长达6的时候，尽管瑞利数已经增加到10¹⁰，仍然没有出现上方的水平流动，而是出现从一边上方穿梭贯穿到另外一边的流动结构。

5. 结论

对垂直入射辐射吸收引发的自然对流进行数值模拟，考虑不同的瑞利数，探究了非平直下壁面的振幅、周期个数以及波长对辐射驱动对流的流动与传热的影响。研究羽流穿透深度、系统温度均匀性以及流函数的变化规律，得到以下结论：

(1) 高度h为1.5，振幅A为0.5时，羽流不能完全穿透，会出现温度分层结构。

(2) 当波长足够长时，波峰温度最高。随着波长逐渐减小，存在过渡情况，当波长L = 1.8时，波谷温度出现翻转现象，此时波峰波谷温度逐渐接近。随着长度进一步减小，波谷附近壁面温度最高。随着振幅的增大，系统的温度均匀性变小。

(3) 波长小的更容易出现上方水平流动的情况，小瑞利数时流函数对称，随着瑞利数增大，流函数图的对称结构逐渐被破坏，上方流线图越来越会往上变得近似水平。

NOTES

^*通讯作者。