摘要:倾斜封闭腔内对流换热问题是非线性非平衡系统中研究的热点问题之一. 本文采用高精度数值方法对倾斜方腔内流体热对流进行了直接数值模拟, 研究了腔体倾角在$0^\circ— 180^\circ$之间变化时, 倾角的不同变化过程对流场非线性演化、传热效率以及流动分岔的影响. 所考虑的Rayleigh数范围为$10^3— 10^6$. 结果表明: 表征传热效率的Nusselt 数对Rayleigh数、Prandtl数及倾斜角度均具有较强依赖性, 在较高Rayleigh数时, Nusselt数会在80°和100°附近产生较大幅度的变化; 高Rayleigh 数下流场及温度场的演变更为复杂, 腔体内存在1—3个对流强度不等的涡卷; 低Rayleigh数下腔体倾角接近90°时流动状态为热传导状态. 当腔体倾角介于$70^\circ— 110^\circ$之间时, 在Rayleigh数$Ra\in(4949,314721)$内存在解的两条稳定分支.
关键词: 热对流/
直接数值模拟/
分岔/
倾角/
高精度

English Abstract

全文HTML

--> --> -->

腔体内自然对流传热有着极其广泛的工程应用, 例如建筑保温、电子仪器的冷却装置、太阳能集热器以及核反应堆设计等. 自1954年Batchelor^[1]对二维封闭腔内自然对流开创性的研究以来, 对于水平放置腔体内自然对流问题, 前人早期进行了理论或实验研究. 如Eckert和Carlson^[2]实验研究了两个具有不同温度的等温垂直板之间的空气层中的温度场, 给出了局部和平均传热与Grashof 数的关系. Patterson和Armfield^[3]研究发现, 随着Rayleigh数$ Ra $的增大, 腔体内对流状态开始分层并且逐渐由定常转向非定常. Xin和Quéré^[4]研究了水平放置的宽高比为4 的长方形腔体中混沌自然对流问题. 对于水平放置的外形为三角形、梯形或平行四边形等特殊腔体内的自然对流换热也有了一定的研究^[5]. 然而在实际应用中, 即使打算将腔体水平放置, 对流系统在自然界中往往是相对于重力倾斜的. Arnold等^[6]研究了倾斜腔体热对流中倾斜角对矩形区域内传热的影响. 之后不久, John和Jorg^[7]的研究也说明了腔体内的流动情况由Rayleigh数($ Ra $)、Prandtl数($ Pr $)、高宽比、倾斜角等参数决定. Khezzar等^[8]分析了初始条件和滞后效应对流场的影响, 认为产生分岔的临界倾斜角度对初始条件有较强的依赖性. Didier等^[9]研究了宽高比为4的二维腔体内Rayleigh-Bénard对流(腔体倾斜角度在$ 0^\circ $—$ 180^\circ $范围内变化)从稳态到混沌的分岔现象, 并将关于倾斜腔内热对流的研究拓展到三维腔体. Torres等^[10]对Prandtl数$ Pr = 5.9 $, Rayleigh数$Ra\leqslant150000$的倾斜立方腔内自然对流进行了三维分岔分析, 讨论了流体从初始分岔到振荡状态过程中倾斜角对立方体封闭空间内三维对流稳定性的影响. 紧接着, 他们证明了三维腔体在$ 0^\circ—90^\circ $转变过程中出现单一稳定涡卷的临界角的存在性, 并将这一结论推广到了任意Prandtl数$ Pr > 0.71 $的Boussinesq流体中^[11]. Miroshnichenko和Sheremet^[12] 对矩形腔内自然对流的研究进行了综述, 并对封闭腔体内自然对流问题的研究进展做了一些总结.
近年来, 随着计算机硬件水平的长足进步以及计算精度的大幅提高, 利用数值模拟研究自然对流已被更多地使用^[13], 这使得****们对自然对流的研究也更加深入. Hamady等^[14]和Kuyper等^[15]模拟研究了某些角度下方腔内的对流情况, 他们将局部和平均Nusselt数($ Nu $)的测量值与数值预测值进行了比较, 进一步说明腔内传热对倾斜角和Rayleigh数有较强的依赖性. Rasoul和Prinos^[16]分别对倾角在$ 0^\circ—180^\circ $和$ 40^\circ—160^\circ $范围的腔体进行了数值模拟, 分析了倾斜角对方形腔体内自然对流的影响, 并给出了倾角对Nu-Ra相关关系的影响. 文献[17-20]对二维腔体内稳定状态的自然对流换热情况进行了研究, 分别模拟了不同高宽比或者倾斜角下腔体内部的流场状态. Armfield和Janssen^[21]描述了$ Pr = 7.5 $, $ Ra = 10^{8} $时腔体内的流场状况, 他们的结论说明了充分展开的热边界层具有与启动流大致相同的稳定性. Zhao和Tian^[22]基于涡量-流函数型Navier-Stokes方程, 给出了一种求解非定常自然对流问题的高精度方法, 模拟研究了方腔内自然对流状态的过渡与分岔. 为了进一步探究不同参数对二维腔体内热对流的影响, 不少****通过改变腔体内流体类型来研究热对流的演变情况. 如Sheremet等^[23]对二维倾斜方腔内氧化铝水纳米流体的流动和传热进行了模拟研究, 他们的研究成果表明边界温度振荡频率的增加会导致平均Nusselt数振荡幅度的增加和振荡周期的减少. Boudjeniba等^[24]模拟研究了铜水纳米流体在倾斜方腔内从层流向混沌自然对流的过渡, 表明流动从稳定层流状态经超临界Hopf分岔后, 依次经历周期状态和准周期状态, 最后进入混沌状态. 最近, Wang等^[25-29]直接数值模拟研究了不同参数对二维腔体内热对流演变过程的影响, 其中Rayleigh数的范围为$10^7\leqslant $$ Ra\leqslant 10^{10}$; 高宽比为$ 0.5 $—$128$, Prandtl数为$0.01\leqslant $$ Pr\leqslant100$, 腔体倾斜角度的范围为$ 0^\circ— 90^\circ $. 他们分析了腔内热对流从层流到湍流的过渡状态、腔内大尺度环流翻转的演变过程等, 发现倾斜引起的对称破坏极大地抑制了方腔内大尺度环流的翻转^[25], 给出了大宽高比二维矩形腔内湍流热对流状态下影响涡卷数量的因素(即涡卷的椭圆形状和黏性阻尼)^[29]. 对于二维腔体, 环流翻转对Ra和Pr有很强的依赖性^{[25,27,28,30]}, 但当$ Ra > 5\times 10^8 $时, 数值模拟结果显示不存在环流翻转现象^[30]. 前人的研究工作集中在对倾斜腔体内流动演化和传热模式的模拟与分析, 很少将腔体倾斜角的各种变化情况综合考虑, 如腔体顺时针旋转与逆时针旋转对流场的影响是否一致等. 事实上, 考虑倾斜角不同的变化过程对流场的影响, 有助于对方腔自然对流问题的对流和传热特性进行更细致全面的分析, 有利于深入揭示自然对流换热的参数依赖性以及流动向湍流过渡的复杂机理. 本文采用高精度数值模拟的方法, 考虑腔体倾斜角在$ 0^\circ $—$ 180^\circ $的范围内, 研究不同角度变化过程对流动和传热的影响, 探讨流场及温度场随倾斜角变化时的分布规律, 并分析不同倾角变化过程下解的分岔情况.

2

-->

2.1.问题描述

本文主要研究倾斜二维方腔内的自然对流问题, 对流示意图如图1, 其中腔体的高为$ H $, 倾角为$ \beta $. 均匀重力场$ {{g}} $为竖直向下的正值. 水平壁面都是绝热的, 垂直壁面均恒温, 且左壁面$ T_{\rm h} $(热壁)温度高于右壁面$ T_{\rm c} $(冷壁)的温度.
图 1 带边界条件的倾斜方腔示意图
Figure1. Schematic diagram of the inclined square cavity with boundary conditions.

2 -->

2.2.控制方程

在本研究中, 假设流体满足Boussinesq近似, 即其密度与温度$ T $呈线性关系, 具体表达式如下:

$ \rho = \rho _0 [1 - \gamma (T - T_0 )], $

其中ρ₀为参考温度$ T_0 $下的流体密度, γ为热膨胀系数. 选取$ H $为特征长度, $ H^{2}/\kappa $为特征时间(其中$ \kappa $为热扩散系数), 引入以下无量纲量:

$\begin{split}&u = \frac{U}{\kappa/H},\; v = \frac{V}{\kappa/H}, \;\theta = \frac{{T - T_{\rm c} }}{{T_{\rm h} - T_{\rm c} }} = \frac{{T - T_{\rm c} }}{{\Delta T}}, \\&p = \frac{{p^* }}{{\rho _0 \kappa^2 /H^2 }}. \end{split}$

在笛卡尔坐标系中, 描述该对流系统的无量纲连续性方程、动量和能量方程如下:

$ \frac{{\partial u}}{{\partial x}} + \frac{{\partial v}}{{\partial y}} = 0, \tag{2a} $

$ \frac{{\partial u}}{{\partial t}} + u\frac{{\partial u}}{{\partial x}} + v\frac{{\partial u}}{{\partial y}} = - \frac{{\partial p}}{{\partial x}} + Pr \nabla ^2 u + RaPr \theta \sin \beta, \tag{2b}$

$ \frac{{\partial v}}{{\partial t}} + u\frac{{\partial v}}{{\partial x}} + v\frac{{\partial v}}{{\partial y}} = - \frac{{\partial p}}{{\partial y}} + Pr \nabla ^2 v + RaPr \theta \cos \beta, \tag{2c}$

$ \frac{{\partial \theta }}{{\partial t}} + u\frac{{\partial \theta }}{{\partial x}} + v\frac{{\partial \theta }}{{\partial y}} = \nabla ^2 \theta, \tag{2d} $

其中$ {{v}} = (u, v) $代表速度场, $ p $ 是流体的压力, $t, \; $$ \theta, \; \rho$和 $ \nu $分别为流体的时间、温度、密度和运动黏度. $ Ra $和$ Pr $分别为Rayleigh数和Prandtl数, 定义为:

$ Ra = \frac{{g\gamma \Delta TH^3 }}{{k\nu }}; \; \; Pr = \frac{\nu }{\kappa}. $

对于速度场, 腔体的所有壁面是无滑移的; 垂直壁面是完全导热的, 且在$ x = 0 $和$ x = H $时温度恒定, 横向壁面绝热. 因此, 边界条件为:
$ u = v = 0 $, $T = T_{\rm h}$ (在$ x = 0 $处),
$ u = v = 0 $, $T = T_{\rm c}$ (对于$ \forall y $, 在$ x = H $处),
$ u = v = 0 $, $\dfrac{{\partial T}}{{\partial y}} = 0$ (在$ y = 0 $处),
$ u = v = 0 $, $\dfrac{{\partial T}}{{\partial y}} = 0$ (对于$ \forall x $, 在$ y = H $处).
引入涡量$ \omega $和流函数$ \psi $来消除压力$ p $, 具体表达式为:

$ \omega = \frac{{\partial v}}{{\partial x}} - \frac{{\partial u}}{{\partial y}}, $

$ u = \frac{{\partial \psi }}{{\partial y}},\; v = - \frac{{\partial \psi }}{{\partial x}}. $

于是, 方程(2a)—(2d)可写为

$ -\omega = \nabla ^2 \psi, \tag{5a}$

$ \begin{split}&\frac{{\partial \omega }}{{\partial t}} + u\frac{{\partial \omega }}{{\partial x}} + v\frac{{\partial \omega }}{{\partial y}} \\=\;& Pr {\nabla ^2}\omega + RaPr \left( {\frac{{\partial \theta }}{{\partial x}}\cos \beta - \frac{{\partial \theta }}{{\partial y}}\sin \beta } \right), \end{split}\tag{5b}$

$ \frac{{\partial \theta }}{{\partial t}} + u\frac{{\partial \theta }}{{\partial x}} + v\frac{{\partial \theta }}{{\partial y}} = \nabla ^2 \theta. \tag{5c}$

相应地, 边界条件也可写成:
$\psi = \dfrac{\partial \psi }{\partial x} = 0$, $\omega = - \dfrac{\partial ^2 \psi }{\partial x^2 }$, $ \theta = 1 $ (在$ x = 0 $处),
$\psi = \dfrac{\partial \psi }{\partial x} = 0$, $\omega = - \dfrac{\partial ^2 \psi }{\partial x^2 }$, $ \theta = 0 $ 　　　　(对$\forall~ y$, 在$ x = 1 $处),
$\psi = \dfrac{\partial \psi }{\partial y} = 0$, $\omega = - \dfrac{\partial ^2 \psi }{\partial y^2 }$, $\dfrac{\partial \theta }{\partial y} = 0$ 　　　　(对$\forall ~x$, 在$ y = 0, 1 $处).
用Nusselt数来表征对流传热效率, 沿与$ y $轴平行的直线上Nusselt数定义为

$ Nu_x = \int_{0}^{1}{\left(u\theta - \frac{{\partial \theta }}{{\partial x}}\right){\rm d}y.} $

令(6)式中$ x = 0 $, 可得热壁Nusselt 数$ Nu_0 $, 其最大值和最小值分别由$ Nu_{\rm {max}} $和$ Nu_{\rm {min}} $表示. 若沿水平方向对$ Nu_x $积分, 可得整个腔体内Nusselt数的平均值:

$ \overline {Nu} = \int_0^1 {Nu_x } {\rm d}x. $

2

-->

3.1.数值方法

本文采用文献[22]中的高精度高分辨率有限差分法来求解控制方程(方程(5a)—(5c))及相应边界条件, 对倾斜腔内热对流进行直接数值模拟. 该方法中, 用优化的三阶迎风紧致格式(Opt-UCD3)离散对流项, 采用四阶对称Padé紧致格式对黏性项进行离散, 利用九点四阶紧致格式^[31]离散流函数泊松方程; 采用TVD 保持的三阶Runge-Kutta方法进行时间推进.
2 -->

3.2.数值检验

为了验证数值方法的有效性和计算结果的正确性, 并说明本文计算中所选取网格尺寸的合理性, 针对不同工作流体($ Pr = 0.71 $, $ Pr = 7.01 $), 在不同倾角($ \beta = 0^\circ $, $ \beta = 45^\circ $)和本文讨论的最高Rayleigh数($ Ra = 10^{6} $)时, 对多个工况进行了网格无关性检验, 并将本文结果与已有文献结果进行了对比.
3 -->

3.2.1.与其他文献的对比

表1将$ Ra = 10^{5} $和$ Ra = 10^{6} $ 时本文计算结果与其他文献结果进行了对比, 其中本文结果的计算网格尺寸为$ {\rm{121}} \times {\rm{121}} $. 结果表明本文所用数值方法的计算结果与文献结果高度符合, 由此, 本文所用数值方法计算结果的正确性得到了验证.

文献	$ \left| \psi \right|_{\rm {max}} $	$ \left| {\psi _{\rm {mid}} } \right| $	$ Nu_{0} $	$ \overline{Nu} $	文献	$ \left| \psi \right|_{\rm {max}} $	$ \left| {\psi _{\rm {mid}} } \right| $	$ Nu_{0} $	$ \overline{Nu} $
$ Ra=10^{5} $					$ Ra=10^{6} $
本文	9.615	9.115	4.520	4.522	本文	16.807	16.383	8.815	8.827
[32]	9.612	9.111	4.509	4.519	[32]	16.750	16.320	8.817	8.800
[33]	—	9.123	4.512	4.522	[33]	—	16.420	8.763	8.829
[34]	9.6173	9.1161	4.5195	—	[34]	16.8107	16.3863	8.8216	—
[35]	9.6202	9.1194	4.5214	—	[35]	16.8411	16.4183	8.8091	—

表1与其他文献结果的对比($ Pr = 0.71 $, $ \beta = {\rm{0}}^\circ $)
Table1.Comparison of the results by different numerical methods for $ Pr = 0.71 $ and $ \beta = {\rm{0}}^\circ $.

3 -->

3.2.2.网格无关性检验

本文在高Rayleigh 数下进行数值计算时选取的网格尺寸为$ 121\times121 $, 此计算的网格分辨率在表2—表5中分别针对$ Pr = 0.71 $和$ Pr = 7.01 $以及$ \beta = 0^\circ $和$ \beta = 45^\circ $(取$ Ra = 10^{6} $), 用$ 241\times241 $的精细网格进行了检验. 结果表明, 在高Rayleigh数下, 本文所取网格尺寸与精细网格下特征参数(流函数$ \psi $和Nusselt数$ Nu $等)的最大相对误差均小于0.2%, 保证了计算结果的准确性; 对于低Rayleigh数($ Ra = 10^{3} $和$ Ra = 10^{4} $), 在保证计算准确的前提下, 为了减少计算量, 本文采用$ 91\times 91 $的网格进行计算. 网格无关性检验结果显示, 低Rayleigh数下$ 91\times91 $网格与$ 241\times241 $精细网格下特征参数的最大相对误差小于0.01%, 与$ 121\times121 $网格下特征参数的最大相对误差小于0.011%.

网格尺寸	$\left| \psi \right|_{\rm {max}}$	误差/%	$\left| {\psi _{\rm {mid}} } \right|$	误差/%	$ Nu_0 $	误差/%
$ 31\times31 $	16.460	2.086	16.118	1.631	9.293	5.301
$ 61\times61 $	16.830	0.119	16.410	0.148	8.798	0.315
$ 91\times91 $	16.802	0.051	16.385	0.002	8.786	0.445
$ 121\times121 $	16.807	0.017	16.383	0.014	8.815	0.119
$ 241\times241 $	16.810	—	16.386	—	8.825	—

表2$ Pr = 0.71 $, $ \beta = {\rm{0}}^\circ $, $ Ra = 10^{6} $下的网格检验结果
Table2.Grid test results for $ Pr = 0.71 $, $ \beta = {\rm{0}}^\circ $ and $ Ra = 10^{6} $.

网格尺寸	$\left| \psi \right|_{\rm {max}}$	误差/%	$\left| {\psi _{\rm {mid}} } \right|$	误差/%	$ Nu_0 $	误差/%
$ 31\times31 $	32.400	3.276	27.974	3.306	9.077	9.345
$ 61\times61 $	33.252	0.734	28.707	0.771	8.332	0.381
$ 91\times91 $	33.438	0.176	28.874	0.195	8.301	0.001
$ 121\times121 $	33.477	0.062	28.911	0.068	8.304	0.039
$ 241\times241 $	33.498	—	28.931	—	8.301	—

表3$ Pr = 0.71 $, $ \beta = {\rm{45}}^\circ $, $ Ra = 10^{6} $下的网格检验结果
Table3.Grid test results for $ Pr = 0.71 $, $ \beta = {\rm{45}}^\circ $ and $ Ra = 10^{6} $.

网格尺寸	$\left| \psi \right|_{\rm {max}}$	误差/%	$\left| {\psi _{\rm {mid}} } \right|$	误差/%	$ Nu_0 $	误差/%
$ 31\times31 $	18.625	5.075	17.873	5.021	9.548	3.514
$ 61\times61 $	19.634	0.067	18.838	0.110	9.195	0.310
$ 91\times91 $	19.609	0.059	18.814	0.020	9.206	0.193
$ 121\times121 $	19.612	0.044	18.812	0.029	9.221	0.037
$ 241\times241 $	19.621	—	18.818	—	9.224	—

表4$ Pr = 7.01 $, $ \beta = {\rm{0}}^\circ $, $ Ra = 10^{6} $下的网格检验结果
Table4.Grid test results for $ Pr = 7.01 $, $ \beta = {\rm{0}}^\circ $ and $ Ra = 10^{6} $.

网格尺寸	$\left| \psi \right|_{\rm {max}}$	误差	$\left| {\psi _{\rm {mid}} } \right|$	误差	$ Nu_0 $	误差
$ {\rm{31}} \times {\rm{31}} $	38.233	6.689%	34.649	6.739%	9.791	7.723%
$ {\rm{61}} \times {\rm{61}} $	40.665	0.752%	36.858	0.793%	9.114	0.271%
$ {\rm{91}} \times {\rm{91}} $	40.902	0.174%	37.090	0.167%	9.089	0.001%
$ {\rm{121}} \times {\rm{121}} $	40.950	0.057%	37.131	0.058%	9.092	0.025%
$ {\rm{241}} \times {\rm{241}} $	40.973	—	37.152	—	9.089	—

表5$ Pr = 7.01 $, $ \beta = 45^\circ $, $ Ra = 10^{6} $下的网格检验结果
Table5.The grid test results for $ Pr = 7.01 $, $ \beta = {\rm{45}}^\circ $ and $ Ra = 10^{6} $.

本文重点研究了腔体倾斜角度($ 0^\circ $—$ 180^\circ $)对流动和传热的影响, 在计算时, 考虑了三种过程.
过程A　模拟腔体倾斜角度$ \beta $在$ 0^\circ—180^\circ $ 范围内腔内流体自然对流换热, 温度场和流场的初始值均取零, 即对每个给定的倾斜角度, 都从零场开始计算.
过程B　模拟腔体倾斜角度逐渐增大的过程. 计算时, $ \beta $以小步长$ \Delta\beta = 10^\circ $逐渐增大, $ \beta $每增大一次, 都以前一次$ \beta $较小时的稳定结果为初值, 直到$ \beta $增大到$ 180^\circ $时, 计算停止.
过程C　与过程B相反, 模拟倾斜角度逐渐减小的过程. 计算时, $ \beta $ 以小步长$ \Delta\beta = 10^\circ $逐渐减小, 以$ \beta = 180^\circ $时的稳定结果为初值, 之后$ \beta $每减小一次, 都以上一次$ \beta $较大时的计算结果为初值, 直到$ \beta $减小到$ 0^\circ $时, 计算停止.
2

-->

4.1.对流传热效率的变化

3 -->

4.1.1.Rayleigh数和Prandtl数的影响

图2为热壁面Nusselt数$ Nu $随Rayleigh数$ Ra $的变化情况(事实上, 整个腔体平均Nusselt数$ \overline{Nu} $具有相同的结论, 故以下只讨论热壁面Nusselt数的情形), 其中红线代表腔内流体为空气时$ Nu $的变化趋势, 蓝线则代表腔内流体为水时$ Nu $的变化趋势. 由图2可以清晰地看出, 随着$ Ra $的增大, Nusselt数也在逐渐增大, 且这种增大的趋势在$ Ra = 10^{5} $附近尤为明显. 另外, 对于不同的流体介质, 腔内自然对流换热的强弱程度略有差异, 且Ra越大, 这种差异越明显. 例如, 在$ Ra = 10^{3} $时, $ Pr = 0.71 $和$ Pr = 7.01 $下Nu的计算结果均约为1.118, 而在$ Ra $增大至$ 10^{6} $时, $ Pr = 0.71 $下$ Nu $的值为8.834, $ Pr = 7.01 $下$ Nu $的值约为 9.216, 表明在较高Rayleigh数下水的传热效率要比空气高一些.
图 2 Nusselt数随Rayleigh数的变化($ \beta = 0^\circ $)
Figure2. Variation of Nusselt number as a function of Rayleigh number ($ \beta = 0^\circ $).

$ Ra $分别为$ 10^{4} $, $ 10^{5} $和$ 10^{6} $时, $ Nu $随普朗特数$ Pr $的变化情况如图3所示, 可以看出, 当$ Ra $较小时$ Nu $随$ Pr $的变化并不明显, 几乎呈水平直线; 在高Rayleigh数下(如$ Ra = 10^{6} $), $ Nu $在$ Pr = 0.2 $到$ Pr = 1.0 $的范围内随$ Pr $增大而增大, 变化较为明显, 但当$ Pr > 1.0 $时, $ Nu $值则稳定在9.25附近.
图 3 Nusselt数随Prandtl数的变化($ \beta = {\rm{0}}^\circ $)
Figure3. Variation of Nusselt number as a function of Prandtl number ($ \beta = {\rm{0}}^\circ $).

3 -->

4.1.2.腔体倾斜角度$ \beta $的影响

研究发现腔体的倾斜角度对腔内传热效率的大小有着明显的影响. 从系统的控制方程来看, 倾斜角度的变化会引起浮力的改变, 这就使得腔内的自然对流换热情况发生改变, 最终导致传热效率也发生变化.
图4给出的是工作流体为空气时($ Pr = 0.71 $)不同Rayleigh数下Nusselt数随腔体倾角变化的曲线. 可以看出, 随着倾斜角度的增大, A, B, C三种过程中$ Nu $均会先增大, 然后开始减小, 大约在$ {\rm{80}}^\circ $—$ {\rm{100}}^\circ $的范围内达到最小值, 随后会继续增大, 接近$ {\rm{180}}^\circ $时又开始减小, 且这种变化趋势关于$ {\rm{90}}^\circ $对称. 但是当$ Ra $增大后, $ Nu $的变化规律也会明显发生改变. 例如在$ Ra = 10^{3} $时, 过程A, B, C中$ Nu $的变化曲线完全重合, 说明低$ Ra $下倾斜角度改变的趋势对$ Nu $的影响并不明显. 但是当$ Ra = 10^{4} $和$ Ra = 10^{5} $时, 过程A, B, C中$ Nu $的变化曲线在$ {\rm{70}}^\circ $—$ {\rm{110}}^\circ $的范围内截然不同, 此时的突变也说明流体的对流状态发生了明显改变. 另外, 在$ {\rm{70}}^\circ $—$ {\rm{110}}^\circ $的范围内, 过程B与过程C中$ Nu $的不同也说明此范围内出现了解的分岔(在4.3节详细讨论), 控制方程本身所具有的非线性性质是产生这种差异的根本原因^[31]. 值得一提的是, 当$ Ra $增大到$ 10^{6} $时, 三种过程中$ Nu $的变化曲线也几乎重合, 这说明在$ Ra $增大到一定程度后, 解的分岔也会消失.
图 4 不同Rayleigh数下Nusselt数随倾斜角度的变化($ Pr = 0.71 $)
Figure4. Variation of Nusselt number with inclination angle for different Rayleigh numbers ($ Pr = 0.71 $).

腔体内工作流体为水($ Pr = 7.01 $)时, Nusselt数随倾斜角度的变化关系具有类似的结论, 但因Rayleigh数不同而略有差异. 当$ Ra = 10^{3}, 10^{4} $和$ 10^{6} $时, 两种工作流体的$ Nu $的变化规律相同(即$Ra = 10^{3}, $$ 10^{6}$时, A, B, C三种过程中$ Nu $的曲线完全重合, $ Ra = 10^{4} $时在$ 80^\circ $—$ 100^\circ $的范围内存在解的分岔); 当$ Ra = 10^{5} $时则与空气中不同, 此时三种过程中$ Nu $曲线完全重合, 不存在解的分岔. 对于给定Rayleigh数, 结合图2可知不同流体下$ Nu $的值不同, 且$ Nu $随$ \beta $的变化也不同(如图5所示). 显然$ Ra = 10^{4} $和$ 10^{5} $时两种流体下$ Nu $随$ \beta $的变化曲线有明显差异. $ Ra = 10^{4} $时(图5左)两种流体下的$ Nu $曲线均在$ 40^\circ $和$ 140^\circ $附近取得极大值, 在$ 100^\circ $附近取得极小值, 但二者的极值大小不同, 如$ Pr = 0.71 $时$ Nu $的最小值为2.3759, 而$ Pr = 7.01 $时最小值为2.4230. 当$ Ra = 10^{5} $时(图5右)两种流体下$ Nu $曲线极值点不同, $ Pr = 0.71 $时极值点为$ \beta = 100^\circ $; 而$ Pr = 7.01 $时极值点为$ \beta = 90^\circ $且曲线关于该点具有对称性. 在该角度下, 空气$ Pr = 0.71 $时$ Nu $值为3.9094, 略高于$ Pr = 7.01 $时的$ Nu $值3.8777, 这表明在某些倾斜角度下, 对于给定的(较大) Rayleigh数, 空气的传热效率可能反而要比水的传热效率高一些.
图 5 不同Prandtl数下Nusselt数随倾斜角度的变化(以过程B为例)
Figure5. Variation of Nusselt number with inclination angle for different Prandtl numbers (Process B).

2 -->

4.2.对流斑图结构

3 -->

4.2.1.低Rayleigh数下对流状态随倾角的变化

图6—图8分别给出了$ Ra = 10^{3},\; 10^{4},\; 10^{5} $时典型倾斜角度下腔体内流场和温度分布($ Pr = 0.71 $). 可以看到, 对于足够大的Rayleigh数, 近热壁处流体沿壁面向上流动, 近冷壁处流体在重力场作用下沿壁面向下流动, 形成对流. 流体微团在运动过程中, 不断将其从热壁吸收得到的热量向近冷壁处流体传递, 直至腔内的热量分布趋于稳定. 由图6可知, $ Ra = 10^{3} $时低Rayleigh数下对流较弱, 流场和温度场的演变相对较平稳; 随着倾角从$ 0^\circ $增大到$ 90^\circ $, 等温线逐渐从弯曲状态转变为垂直状态, 这说明流体传递热量的方式由以对流换热为主逐渐转换为以热传导为主. 当倾角进一步增大时, 流场线及等温线的分布发生翻转, 与$ 0^\circ $—$ 90^\circ $的分布正好相反(如$ \beta = 60^\circ $与$ \beta = 140^\circ $). 在此Rayleigh数下, 不论倾角如何, 流场中只有一个主涡卷, 流函数的值介于(–1.2, 1.1)之间.
图 6 $Ra=10^{3}$时流线与等温线图
Figure6. Streamlines and isotherms for $Ra=10^{3}$.

图 8 $ Ra = 10^{5} $时流线与等温线图
Figure8. Streamlines and isotherms for $ Ra = 10^{5} $.

当$ Ra = 10^{4} $时, 腔内流场和温度分布如图7所示, 由等温线的弯曲程度可以看出, 随着$ Ra $数的增加, 对流明显增强. 随着倾角$ \beta $从$ 0^\circ $增大到$ 90^\circ $, 对流强度先增后减, $ \beta = 40^\circ $时对流最强. 若腔体继续从$ 90^\circ $倾斜到$ 180^\circ $, 流场线以及等温线的分布与$ 0^\circ—90^\circ $的分布相反. 在此Rayleigh数下, 无论何种倾角, 流场中仍然只有一个主涡卷, 流函数的值介于(–7, 7)之间.
图 7 $Ra=10^{4}$时流线与等温线图
Figure7. Streamlines and isotherms for $Ra=10^{4}$.

当$ Ra = 10^{5} $时, 腔内流场和温度分布如图8所示, 流体的对流运动比$ Ra = 10^{4} $时更加剧烈. 在$ \beta = 0^\circ $时, 流场仍具有明显的中心对称性, 但中心主涡分裂成两个向垂直壁面方向移动的二次涡, 在腔体的左上角和右下角出现再循环区域, 对流增强. 当倾角增大到$ \beta = 10^\circ $时, 从主涡分裂出的两个同向旋转二次涡分别向腔体左上角和右下角移动, 此时对流最强. 当倾角继续增大到$ \beta = 20^\circ $时, 两个二次涡又逐渐合并最终融入主涡之中, 流动回到只有一个主涡的状态, 但由于重力作用在温度梯度垂直方向上作用较大, 涡卷被拉伸为扁平状, 这种扁平状状态大约持续到$ \beta = 80^\circ $. 随着倾角进一步增大到$ 90^\circ $, 重力场在温度梯度方向上的作用最大, 使涡卷从扁平状逐渐变为圆形, 对流强度减弱. 当倾角$ \beta $在$ 90^\circ $—$ 180^\circ $范围内进一步增大时, 流体的对流状态和等温线的弯曲情况与$ 0^\circ $—$ 90^\circ $的状态刚好相反(如$ \beta = 10^\circ $与$ \beta = 170^\circ $), 此Rayleigh数下流函数$ \psi $的值介于(–24, 24)之间.
3 -->

4.2.2.高Rayleigh数下对流状态随倾角的变化

随着Rayleigh数的增大, 对流越来越强, 对流随腔体倾斜角度的变化也更加复杂. 图9给出了$ Ra = 10^{6} $时流场随倾角$ \beta $从$ 0^\circ $变化到$ 180^\circ $时的详细演变过程. 此时对流状态的变化较低Rayleigh数时要复杂得多, 在$ \beta = 0^\circ $时, 流线图中虽然只有一组闭合流线, 但是靠近冷热壁的流场线弯曲的程度较大, 有分化为两组闭合流线的趋势, 到$ \beta = 10^\circ $时流线图中已经分裂出两组闭合流线, 主涡卷内包含两个同向旋转的二次涡卷. 随着倾角的进一步增大, 二次涡卷分别向左上和右下移动, 在$ \beta = 30^\circ $时两个同向旋转的二次涡卷之间出现了反向旋转的三次涡. 这种具有一个主涡、两个同向旋转二次涡和一个反向旋转三次涡的对流状态一直持续到$ \beta = 70^\circ $, 但随着倾角的增大, 由于重力在温度梯度方向上的作用增大, 将对流涡卷在此方向上进行了拉伸, 致使涡卷从扁平状逐渐变为圆形. 随着倾角进一步增大到$ \beta = 80^\circ $, 二次涡和三次涡被主涡吸收合并为一个大涡卷, 流场演变为只有一组闭合流线的状态, 此状态一直持续到$ \beta = 100^\circ $. 流场在倾角为$ 110^\circ $—$ 180^\circ $范围内的变化情况与在$ 0^\circ $—$ 70^\circ $范围内的变化情况刚好相反, 这是腔体旋转$ 180^\circ $之后冷热壁面互换导致的结果. 在此Rayleigh数下, 流函数$ \psi $的值介于(–70, 75)之间.
图 9 $ Ra = 10^{6} $时流场结构随倾角的变化
Figure9. Variation of flow field with inclination angle for $ Ra = 10^{6} $.

2 -->

4.3.解的分岔

从4.1.2节知道, 对于相同的Rayleigh数和倾角, 采用不同的过程(如过程B和过程C)得到的流场稳定状态不同, 对应的对流传热Nusselt数不同, 表明存在解的分岔现象.
3 -->

4.3.1.Nusselt 数的变化

进一步模拟发现, 当$ Ra = 10^{4} $和$ 10^{5} $, 且倾角介于$ 70^\circ $—$ 110^\circ $之间时, 过程B和过程C的$ Nu $曲线不重合, 出现了分岔; 但是当$ Ra = 10^{6} $ 或$ 10^{3} $时, 这种分岔现象又消失了. 图10(a)给出了$ Ra = 4949 $时, 过程B与过程C中Nusselt数$ Nu $随倾角的变化情况. 可以看出, 两个过程的Nusselt数曲线完全重合, 此时并没有出现解的分岔现象. 而当$ Ra = 4950 $时, 如图10(b)所示, 过程B与过程C的$ Nu $曲线在$ 70^\circ $—$ 110^\circ $之间出现了分岔, 在此倾角区间内存在两条解的分支, 且二者的$ Nu $曲线是对称的(对称轴为$ 90^\circ $). 由此得到分岔出现的临界$ Ra $数下限值为$ Ra = 4949 $. 当$Ra\leqslant4949$ 时, 腔内流场变化情况相对稳定, 初始条件对流场变化无明显影响, 过程B和过程C中传热效率的变化情况几乎完全相同, 流场的稳定状态是唯一的. 当$ Ra > 4949 $时, 流场最终稳定状态对初始条件有很强的依赖性, 流场的非线性演变以及对流传热效率也不同.
图 10 分岔点附近Nusselt数随倾角的变化　(a) Ra = 4949; (b) Ra = 4950
Figure10. Variation of Nusselt number with inclination angle near the bifurcation point: (a) Ra = 4949; (b) Ra = 4950.

持续增大Rayleigh数直到$ Ra = 314061 $, 对流越来越复杂, 过程B和过程C的$ Nu $曲线仍在$ 70^\circ $—$ 110^\circ $之间出现了分岔, 且二者关于$ 90^\circ $对称; 当$314062\leqslant Ra < 314721$时, 分岔存在的倾角区间缩短到[$ 70^\circ $, $ 90^\circ $]之间, 两种过程的$ Nu $曲线出现了一种不对称现象. 此时过程B中$ Nu $曲线关于$ 90^\circ $对称, 而过程C中$ Nu $曲线不具有对称性, 它在$ 80^\circ $处产生跳跃, 如图11(a)所示. 进一步增大Rayleigh数到$ Ra = 314721 $时, 如图11(b)所示, 过程B与过程C的$ Nu $曲线已完全重合, 解的分岔消失, 分岔区间的上限为$ Ra = 314721 $. 综上可知, 随着倾角变化过程的不同, 分岔现象存在的Rayleigh数范围为$ Ra\in(4949, 314721) $.
图 11 分岔区间上界附近Nusselt数随倾角的变化曲线　(a) $Ra=314720 $; (b) Ra = 314721
Figure11. Variation of Nusselt number with inclination angle near the upper bound of bifurcation interval: (a) $Ra=314720$; (b) Ra = 314721.

3 -->

4.3.2.对流斑图结构的比较

图12给出了腔体倾斜角为$ 80^\circ $时给定Rayleigh数下的流场斑图结构. $ Ra = 4950 $时, 过程B对应的对流斑图(图12(a))只有一个主涡, 对流涡卷顺时针旋转, 主涡中心位于腔体几何中心, 流函数值为–4.412; 而过程C的对流斑图(图12(b))显示流场中除了中心处的大尺度主涡外, 在腔体的左上角和右下角出现了流动强度较弱的二次涡, 主涡的中心依然位于几何中心, 流函数值为2.122. 尽管二次逆向旋转的涡卷强度很弱, 但减缓了主涡的旋转速度, 降低了过程C的对流传热效率. 主涡中心处流函数值的大小关系表明, 过程B比过程C中的对流更强, 对流传热能力也更强, 对应于图10(b)中$ \beta = 80^\circ $时过程B的Nusselt数明显大于过程C的Nusselt数. 当$ Ra $数上升到314720时, 过程B与过程C的对流斑图区别更加明显, 过程B对应的对流斑图(图12(c))显示腔体左下角和右上角出现了极微弱的二次涡, 强度仅约为主涡的1/750; 过程C中二次涡的强度明显更强(约为主涡强度的1/14), 二次涡的成长将中心主涡向内挤压, 使得中心主涡变成了椭圆结构, 如图12(d)所示. 虽然两种过程中, 流场均出现了二次逆向涡, 但是过程C中更强的逆向二次涡大大降低了对流传热的能力, 致使过程C的对流传热能力较过程B更弱, 如图11(b)中$ \beta = 80^\circ $处. 两种过程的主涡中心位置均位于腔体几何中心, 过程B和过程C中腔体中心处流函数值分别为–45.273和41.450. 无论流场结构的差异还是对流传热的能力差异等方面均表明, 在$ Ra\in(4949, 314721) $的层流阶段, 存在解的分岔现象. 有趣的是, 若继续增大$ Ra $ 数, 这种分岔现象会消失, 但是随着$ Ra $增大到一定程度, 到达湍流状态时, Wang等^[26,28,29]近期的研究表明这种解的分岔现象, 或者说多态现象将会再次出现. 关于湍流状态下多态现象的研究, 感兴趣的读者可参阅文献[26,28,29]. 综合本文结果和文献结果, 我们希望未来对层流状态下多涡卷的稳态结构进行深入研究, 并从每一种稳态结构出发逐渐发展至湍流状态, 这或许能更深入地揭示层流状态和湍流状态下影响涡卷数量变化和流动翻转等的原因.
图 12 $ \beta = 80^\circ $时分岔点处过程B与过程C的对流斑图
Figure12. Flow field of process B and process C at the bifurcation point for $ \beta = 80^\circ $.

本文采用高精度高分辨率数值方法对二维倾斜方腔内自然对流进行了数值模拟研究, 讨论了不同倾角变化过程对流动状态和传热效率的影响, 得到以下主要结论.
1) 在给定的腔体倾角$ \beta $下, Nusselt数$ Nu $总是随着Rayleigh数$ Ra $的增大而增大. 若给定Ra, 随着倾角从$ {\rm{0}}^\circ $增大到$ {\rm{180}}^\circ $(或从$ {\rm{180}}^\circ $减小到$ {\rm{0}}^\circ $), $ Nu $随之先增大后减小再增大, 最后再次减小; 对于低Rayleigh数, $ Nu $曲线关于$ \beta = {\rm{90}}^\circ $对称.
2) 当腔体水平放置时, 对于不同的流体介质, 腔内自然对流换热的强弱程度不同, 且$ Ra $越大, 这种差异越明显. 此外, 腔体水平放置时, 相同Rayleigh数下不同流体介质的对流传热效率不同, 且流体的Prandtl数越大则传热越强, 例如流体介质为水时的传热效率大于流体介质为空气时的情况. 但在某些特殊倾角下, 结论可能相反, 这为实际应用提供了一些思路.
3) 腔体倾角对流动状态有很大的影响, 特别是对于高Rayleigh数$ Ra = 10^{6} $, 倾角$ \beta $在$ {\rm{10}}^\circ—{\rm{70}}^\circ $内, 腔体内会出现二次涡甚至三次涡, 对流传热明显增强.
4) 倾角在$ {\rm{70}}^\circ $—$ {\rm{110}}^\circ $之间时, 该对流系统存在解的分岔现象, 分岔现象存在的Rayleigh数范围为$ Ra\in(4949, 314721) $. 在该范围的上限附近有一个过渡区域[314062, 314721), 此区域内解的分岔出现的倾角范围缩小为$ {\rm{70}}^\circ $—$ {\rm{90}}^\circ $.