高级搜索

基于格子Boltzmann方法的多热源受限圆腔内热对流的非线性特征演化

, 

downloadPDF

上一篇

下一篇

谢淦昊, 赵明. 基于格子Boltzmann方法的多热源受限圆腔内热对流的非线性特征演化[J]. 物理学报, 2026, 75(16): 160001. doi: 10.7498/aps.75.20260278
引用本文: 谢淦昊, 赵明. 基于格子Boltzmann方法的多热源受限圆腔内热对流的非线性特征演化[J]. 物理学报, 2026, 75(16): 160001. doi: 10.7498/aps.75.20260278
shu
Ganhao XIE, Ming ZHAO. Evolution of nonlinear characteristics in thermal convection within a confined circular cavity with multiple heat sources based on lattice Boltzmann method[J]. Acta Physica Sinica, 2026, 75(16): 160001. doi: 10.7498/aps.75.20260278
Citation: Ganhao XIE, Ming ZHAO. Evolution of nonlinear characteristics in thermal convection within a confined circular cavity with multiple heat sources based on lattice Boltzmann method[J]. Acta Physica Sinica, 2026, 75(16): 160001. doi: 10.7498/aps.75.20260278
shu

基于格子Boltzmann方法的多热源受限圆腔内热对流的非线性特征演化

  • 上海理工大学能源与动力工程学院, 上海　200093

    • 通讯作者: E-mail: lightzm@usst.edu.cn.

Evolution of nonlinear characteristics in thermal convection within a confined circular cavity with multiple heat sources based on lattice Boltzmann method

  • School of Energy and Power Engineering, University of Shanghai for Science and Technology, Shanghai 200093, China

    • Corresponding author: E-mail: lightzm@usst.edu.cn.

  • 摘要: 受限空间热对流的流动与传热特性具有对空间几何构型的强烈依赖性和复杂性. 本文基于格子Boltzmann方法针对菱形分布的4个等温热源的圆形封闭腔内热对流进行数值研究, 并通过相空间重构、功率谱分析及最大李雅普诺夫指数进行非线性特征演化规律的研究, 刻画出所研究系统的混沌演化路径. 研究发现, 该系统从对称稳态出发, 依次经历超临界叉式分岔和超临界Hopf分岔, 进入非对称稳态和单周期振荡状态. 随着瑞利数进一步增大, 系统发生异宿分岔引发异宿爆炸而转变为混沌态, 并伴随热源附近流型的转变. 随后, 系统进入一个拟周期窗口, 并在该窗口内观测到锁频现象, 以及频率解锁与环面破裂, 最终再次进入混沌状态.
    Abstract: The flow and heat transfer characteristics of thermal convection in confined spaces exhibit strong dependence on and complexity due to the spatial geometric configuration. This study numerically investigates the nonlinear dynamics and bifurcation transitions of natural convection in a two-dimensional circular enclosure containing four isothermal circular heat sources arranged in a diamond configuration, with the enclosure wall maintained at a constant cold temperature. The lattice Boltzmann method with a double-distribution function model (D2G9 for the flow field and D2Q4 for the temperature field) is employed to solve the governing equations under the Boussinesq approximation. Numerical simulations are conducted over a range of Rayleigh numbers (Ra) from 8.5×104 to 9.5×105 at a fixed Prandtl number of Pr = 0.706. The flow regime transitions are identified through temperature contours and streamlines, while the system state transitions and bifurcation types are determined by analyzing time series, phase space trajectories, power spectral density (PSD), and the maximum Lyapunov exponent.The results reveal a rich and complex bifurcation cascade as Ra increases. At low Ra, the flow remains in a steady, symmetric state. A supercritical pitchfork bifurcation occurs at Ra ≈ 9.5×104, where the symmetric steady state loses stability, giving rise to a pair of stable, mirror-image asymmetric steady solutions. This symmetry breaking is manifested by the dominance of a single vortex near the top heat source. As Ra increases further to Ra ≈ 1.1×105, the system undergoes a supercritical Hopf bifurcation, transitioning from a steady to a stable periodic oscillatory state represented by a limit cycle in phase space. This single-frequency oscillation is driven by the periodic expansion and contraction of the dominant vortex near the top heat source. Between Ra ≈ 1.75×105 and 3×105, the system maintains single-periodic oscillations, with vortices periodically generated near the top heat source and moving toward one side. A fundamental change in dynamics occurs at Ra ≈ 3.15×105, where a global heteroclinic bifurcation is identified. The limit cycle destabilizes, and the system transitions directly into chaos, confirmed by a broadband PSD and a positive maximum Lyapunov exponent. Phase space reconstruction reveals chaotic trajectories intermittently visiting regions corresponding to the two unstable mirror-image states, forming a heteroclinic network. As Ra increases through the chaotic regime, the flow structure near the top heat source reorganizes from intermittent dual-vortex motion to a persistent central thermal plume, causing the collapse of the chaotic attractor's topology. Remarkably, at Ra ≈ 7.5×105, the chaotic state destabilizes via a torus bifurcation, leading to a quasi-periodic window characterized by two incommensurate frequencies in the PSD. This is followed by a frequency-locking phenomenon at Ra ≈ 8.75×105, where the system reverts to a periodic state. The locked state eventually breaks down as Ra increases, re-entering quasi-periodicity and finally transitioning back to a fully developed chaotic state through torus breakdown at Ra ≈ 9.5×105.In conclusion, this study elucidates a complex and non-standard route to chaos for multi-source natural convection in a circular enclosure, involving successive steady symmetry breaking, Hopf bifurcation, a direct transition to chaos via heteroclinic bifurcation, and a subsequent quasi-periodic window with frequency locking. The findings highlight the profound influence of the interplay between multiple discrete heat sources and curved boundaries on the system’s nonlinear dynamics.
  • 加载中

  • 图 1  物理模型示意图

    Figure 1.  Physical model diagram.

    下载: 全尺寸图片 幻灯片

    图 2  曲边边界示意图

    Figure 2.  Schematc diagram of curved edges model.

    下载: 全尺寸图片 幻灯片

    图 3  不同网格$ \overline{Nu} $相空间重构　(a) $ Ra=4\times {10}^{5} $; (b) $ Ra=8\times {10}^{5} $

    Figure 3.  Phase space reconstruction diagram of $ \overline{Nu} $ from different grids: (a) $ Ra=4\times {10}^{5} $; (b) $ Ra=8\times {10}^{5} $.

    下载: 全尺寸图片 幻灯片

    图 4  不同Ra的温度云图　(a) $ Ra=8.5\times {10}^{4} $; (b) $ Ra=9\times {10}^{4} $; (c) $ Ra=9.5\times {10}^{4} $; (d) $ Ra=1.75\times {10}^{5} $; (e) $ Ra=3.1\times $$ {10}^{5} $; (f) $ Ra=6\times {10}^{5} $; (g) $ Ra=8\times {10}^{5} $; (h) $ Ra=9\times {10}^{5} $

    Figure 4.  Temperature contours at different Ra: (a) $ Ra=8.5\times {10}^{4} $; (b) $ Ra=9\times {10}^{4} $; (c) $ Ra=9.5\times {10}^{4} $; (d) $ Ra=1.75\times {10}^{5} $; (e) $ Ra=3.1\times {10}^{5} $; (f) $ Ra=6\times {10}^{5} $; (g) $ Ra=8\times {10}^{5} $; (h) $ Ra=9\times {10}^{5} $.

    下载: 全尺寸图片 幻灯片

    图 5  不同Ra的流线图　(a) $ Ra=8.5\times {10}^{4} $; (b) $ Ra=9\times {10}^{4} $; (c) $ Ra=9.5\times {10}^{4} $; (d) $ Ra=1.75\times {10}^{5} $; (e) $ Ra=3.1\times $$ {10}^{5} $; (f) $ Ra=6\times {10}^{5} $; (g) $ Ra=8\times {10}^{5} $; (h) $ Ra=9\times {10}^{5} $

    Figure 5.  Streamlines at different Ra: (a) $ Ra=8.5\times {10}^{4} $; (b) $ Ra=9\times {10}^{4} $; (c) $ Ra=9.5\times {10}^{4} $; (d) $ Ra=1.75\times {10}^{5} $; (e) $ Ra= $$ 3.1\times {10}^{5} $; (f) $ Ra=6\times {10}^{5} $; (g) $ Ra=8\times {10}^{5} $; (h) $ Ra=9\times {10}^{5} $.

    下载: 全尺寸图片 幻灯片

    图 6  速度分量u的极值随瑞利数Ra的分岔图

    Figure 6.  Bifurcation diagram of the extrema of velocity component $ u $ versus Rayleigh number Ra.

    下载: 全尺寸图片 幻灯片

    图 7  不同Ra下$ \overline{Nu} $的时间演化

    Figure 7.  Time evolutions of $ \overline{Nu} $ at different Ra.

    下载: 全尺寸图片 幻灯片

    图 8  不同Ra下$ \overline{Nu} $的时间演化

    Figure 8.  Time evolutions of $ \overline{Nu} $ at different Ra.

    下载: 全尺寸图片 幻灯片

    图 9  (a) $ \overline{Nu} $在$ Ra\approx 1.1\times {10}^{5} $和$ Ra=1.75\times {10}^{5} $下的PSD; (b) $ \overline{Nu} $在$ Ra=2.75\times {10}^{5} $和$ Ra=3\times {10}^{5} $下的PSD

    Figure 9.  (a) PSD of $ \overline{Nu} $ at $ Ra\approx 1.1\times {10}^{5} $ and $ Ra=1.75\times $$ {10}^{5} $; (b) PSD of $ \overline{Nu} $ at $ Ra=2.75\times {10}^{5} $ and $ Ra= $$ 3\times {10}^{5} $.

    下载: 全尺寸图片 幻灯片

    图 10  不同Ra下$ \overline{Nu} $相空间重构

    Figure 10.  Phase space reconstruction diagram of $ \overline{Nu} $ at different Ra.

    下载: 全尺寸图片 幻灯片

    图 11  不同Ra下中心点$ u\text-v $相图

    Figure 11.  $ u\text-v $ phase diagram of center point at different Ra.

    下载: 全尺寸图片 幻灯片

    图 12  不同Ra中心点$ u\text-v $相图　(a) $ Ra=4\times {10}^{5} $; (b) $ Ra= 6\times {10}^{5} $

    Figure 12.  $ u\text-v $ phase diagram of center point at different Ra: (a) $ Ra=4\times {10}^{5} $; (b) $ Ra=6\times {10}^{5} $.

    下载: 全尺寸图片 幻灯片

    图 13  在$ Ra=4\times {10}^{5} $和$ Ra=6\times {10}^{5} $下$ \overline{Nu} $的PSD

    Figure 13.  PSD of $ \overline{Nu} $ at $ Ra=4\times {10}^{5} $ and $ Ra=6\times {10}^{5} $.

    下载: 全尺寸图片 幻灯片

    图 14  不同Ra的$ u $相空间重构　(a) $ Ra = 4 \times {10}^{5} $; (b) $ Ra= $$ 6\times {10}^{5} $

    Figure 14.  Phase space reconstruction diagram of $ u $ at different Ra: (a) $ Ra=4\times {10}^{5} $; (b) $ Ra=6\times {10}^{5} $.

    下载: 全尺寸图片 幻灯片

    图 15  不同Ra下$ \overline{Nu} $数的时间演化

    Figure 15.  Time evolutions of $ \overline{Nu} $ at different Ra.

    下载: 全尺寸图片 幻灯片

    图 16  不同Ra数下$ \overline{Nu} $相空间重构

    Figure 16.  Phase space reconstruction diagram of $ \overline{Nu} $ at different Ra.

    下载: 全尺寸图片 幻灯片

    图 17  不同Ra的$ \overline{Nu} $的PSD

    Figure 17.  PSD of $ \overline{Nu} $ for different Ra.

    下载: 全尺寸图片 幻灯片

    图 18  不同Ra下$ \overline{Nu} $的时间演化

    Figure 18.  Time evolutions of $ \overline{Nu} $ at different Ra.

    下载: 全尺寸图片 幻灯片

    图 19  $ Ra=9.5\times {10}^{5} $下$ \overline{Nu} $的PSD

    Figure 19.  PSD of $ \overline{Nu} $ at $ Ra=9.5\times {10}^{5} $.

    下载: 全尺寸图片 幻灯片

    图 20  $ Ra=8.75\times {10}^{5} $时各热源$ \overline{Nu} $的时间演化

    Figure 20.  $ \overline{Nu} $ evolution of individual heat sources at $ Ra=8.75\times {10}^{5} $.

    下载: 全尺寸图片 幻灯片

    图 21  $ Ra=8.75\times {10}^{5} $时各热源$ \overline{Nu} $的PSD

    Figure 21.  PSD of $ \overline{Nu} $ for individual heat sources at $ Ra=8.75\times {10}^{5} $.

    下载: 全尺寸图片 幻灯片

    表 1  $ \overline{{K}_{\text{eq}}} $与文献对比结果

    Table 1.  Comparison of $ \overline{{K}_{\text{eq}}} $ with literature results.

    $ Ra $ 本文 文献[17] 误差/%
    $ 6.19\times {10}^{4} $ 3.32 3.30 0.50
    $ 2.25\times {10}^{5} $ 4.39 4.55 3.57
    $ 1.06\times {10}^{6} $ 6.47 6.52 0.72
    $ 1.30\times {10}^{6} $ 6.69 6.80 1.58
    $ 1.70\times {10}^{6} $ 7.05 7.15 1.39
    $ 2.14\times {10}^{6} $ 7.44 7.44 0.001
    下载: 导出CSV

    表 2  $ \overline{Nu} $网格无关性验证结果

    Table 2.  Results of grid independence test for $ \overline{Nu} $.

    Ra $ 201\times 201 $ $ 401\times 401 $ $ 801\times 801 $
    $ \overline{Nu} $ 相对误差/% $ \overline{Nu} $ 相对误差/% $ \overline{Nu} $ 相对误差/%
    $ 8\times {10}^{4} $
    $ 4\times {10}^{5} $
    $ 8\times {10}^{5} $    		 3.52
    5.40
    6.36    		 0.643
    2.03
    3.89    		 3.51
    5.29
    6.11    		 0.084
    1.07
    0.70    		 3.50
    5.23
    6.07
    —
    下载: 导出CSV

    表 3  不同Ra下的LLE

    Table 3.  LLE at different Ra.

    Ra LLE Ra LLE
    $ 1.1\times {10}^{5} $ –0.0021 $ 7\times {10}^{5} $ –0.004
    $ 1.75\times {10}^{5} $ 0.0019 $ 8.25\times {10}^{5} $ 0.0038
    $ 2.75\times {10}^{5} $ 0.0099 $ 8.5\times {10}^{5} $ 0.0003
    $ 3\times {10}^{5} $ 0.0036 $ 8.75\times {10}^{5} $ 0.0019
    $ 4\times {10}^{5} $ 0.0642 $ 9.2\times {10}^{5} $ 0.0043
    $ 6\times {10}^{5} $ 0.4274 $ 9.5\times {10}^{5} $ 0.1574
    下载: 导出CSV
  • [1] 王文博 2024 博士学位论文(昆明: 昆明理工大学) Wang W B 2024 Ph. D. Dissertation (Kunming: Kunming University of Science and Technology
    [2] 王晓梦 2024供热工程建设与高效运行研讨会 济南, 2024年4月24—26日 第134—138页 Wang X M 2024 Seminar on Heating Engineering Construction and Efficient Operation Jinan, China, 2024 pp134–138
    [3] 杨浩, 张销杰, 黄荣宗 2023 物理学报 72 134401 doi: 10.7498/aps.72.20230151 Yang H, Zhang X J, Huang R Z 2023 Acta Phys. Sin. 72 134401 doi: 10.7498/aps.72.20230151
    [4] Pandey S, Park Y G, Ha M Y 2019 Int. J. Heat Mass Transfer 138 762 doi: 10.1016/j.ijheatmasstransfer.2019.04.097
    [5] Angeli D, Barozzi G S, Collins M W, Kamiyo O M 2010 Int. J. Therm. Sci. 49 2231 doi: 10.1016/j.ijthermalsci.2010.08.002
    [6] Mizushima J, Hayashi S, Adachi T 2001 Int. J. Heat Mass Transfer 44 1249 doi: 10.1016/S0017-9310(00)00188-5
    [7] Rashid M H, Xu F 2025 Int. J. Therm. Sci. 214 109913 doi: 10.1016/j.ijthermalsci.2025.109913
    [8] 尹慧, 赵秉新 2021 物理学报 70 11 doi: 10.7498/aps.70.20201513 Yin H, Zhao B X 2021 Acta Phys. Sin. 70 11 doi: 10.7498/aps.70.20201513
    [9] 余端民, 赵明 2019 上海理工大学学报 41 2 doi: 10.13255/j.cnki.jusst.2019.02.002 Yu D M, Zhao M 2019 J. Univ. Shanghai Sci. Technol. 41 2 doi: 10.13255/j.cnki.jusst.2019.02.002
    [10] 赵明, 王柯, 余端民 2020 计算物理 37 667 doi: 10.19596/j.cnki.1001-246x.8159 Zhao M, Wang K, Yu D M 2020 Chin. J. Comput. Phys. 37 667 doi: 10.19596/j.cnki.1001-246x.8159
    [11] Mun G S, Park Y G, Yoon H S, Kim M, Ha M Y 2017 Int. J. Heat Mass Transfer 114 1234 doi: 10.1016/j.ijheatmasstransfer.2017.04.004
    [12] Seo Y M, Mun G S, Park Y G, Ha M Y 2017 Int. J. Heat Mass Transfer 114 1249 doi: 10.1016/j.ijheatmasstransfer.2017.06.101
    [13] 何雅玲, 王勇, 李庆 2009 格子Boltzmann方法的理论及应用(第1版) (北京: 科学出版社)第1, 141页 He Y L, Wang Y, Li Q 2009 Theory and Application of Lattice Boltzmann Method (1st Ed.) (Beijing: Science Press) pp1, 141
    [14] Guo Z L, Shi B C, Zheng C G 2002 Int. J. Numer. Methods Fluids 39 325 doi: 10.1002/fld.337
    [15] Guo Z L, Zheng C, Shi B 2002 Phys. Fluids 14 2007 doi: 10.1063/1.1471914
    [16] Li Z, Yang M, Zhang Y 2016 Int. J. Heat Mass Transfer 94 222 doi: 10.1016/j.ijheatmasstransfer.2015.11.042
    [17] Kuehn T H, Goldstein R J 1976 J. Fluid Mech. 74 695 doi: 10.1017/S0022112076002012
    [18] Wolf A, Swift J B, Swinney H L, Vastano J A 1985 Physica D 16 285 doi: 10.1016/0167-2789(85)90011-9
    [19] Rao Y F, Miki Y, Fukuda K, Takata Y, Hasegawa S 1985 Int. J. Heat Mass Transfer 28 705 doi: 10.1016/0017-9310(85)90193-0
    [20] Angeli D, Corticelli M A, Fichera A, Pagano A 2017 Int. Commun. Heat Mass Transfer 88 139 doi: 10.1016/j.icheatmasstransfer.2017.08.018
  • 加载中
图( 21) 表( 3)
计量
  • 文章访问数:  36
  • HTML全文浏览数:  36
  • PDF下载数:  1
  • 施引文献:  0
出版历程
  • 收稿日期:  2026-02-26
  • 刊出日期:  2026-08-20

基于格子Boltzmann方法的多热源受限圆腔内热对流的非线性特征演化

    通讯作者: E-mail: lightzm@usst.edu.cn.
  • 上海理工大学能源与动力工程学院, 上海　200093
  • 收稿日期:  2026-02-26

摘要: 受限空间热对流的流动与传热特性具有对空间几何构型的强烈依赖性和复杂性. 本文基于格子Boltzmann方法针对菱形分布的4个等温热源的圆形封闭腔内热对流进行数值研究, 并通过相空间重构、功率谱分析及最大李雅普诺夫指数进行非线性特征演化规律的研究, 刻画出所研究系统的混沌演化路径. 研究发现, 该系统从对称稳态出发, 依次经历超临界叉式分岔和超临界Hopf分岔, 进入非对称稳态和单周期振荡状态. 随着瑞利数进一步增大, 系统发生异宿分岔引发异宿爆炸而转变为混沌态, 并伴随热源附近流型的转变. 随后, 系统进入一个拟周期窗口, 并在该窗口内观测到锁频现象, 以及频率解锁与环面破裂, 最终再次进入混沌状态.

English Abstract

    全文HTML

    1.   引　言

    • 不依赖外部动力传热机制中的热对流(或自然对流)广泛存在于电子冷却[1]、建筑节能[2]、相变储能[3]等紧凑空间工程场景中. 受限空间内热对流因腔体外形、热源布局及浮升力的强耦合作用, 其流动和传热特征呈现复杂性, 且具有强烈的构型依赖性[4]. 探究其非线性动力学特征的演化行为及流动模式的转变, 不仅成为揭示其深层物理机制的核心, 也可以帮助人们加深对复杂流动的理解.

      在如方腔、环形腔等简单构型中, 其随瑞利数(Rayleigh number, Ra)主导的流动失稳路径已被系统研究. 如Angeli等[5]综述了受限空间内热对流的非线性行为通常会通过Ruelle-Takens, Feigenbaum和Pomeaue-Manneville三条路径过渡到混沌. Mizushima等[6]基于线性稳定性理论研究了不同半径比对环腔的分岔类型转变. Rashid等[7]通过有限体积法进行数值模拟, 利用相图和功率谱密度(power spectral density, PSD)刻画了侧壁分别为冷源和热源的方腔随Ra变化的Feigenbaum路径. 尹慧和赵秉新[8]采用有限差分法对不同倾斜角方腔自然对流进行非线性分析, 发现在倾斜角为70°—90°时存在解的分岔. 上述研究初步揭示了具有简单构型的受限空间内热对流失稳的影响机制, 针对复杂多热源热对流系统的非线性动力学研究还有待补充.

      系统内部存在多个离散热源时, 其动力学复杂性显著增强. 余端民和赵明[9]利用格子Boltzmann方法(lattice Boltzmann method, LBM)研究了不同开缝度对双羽流生成、合并机制及系统整体换热能力的影响规律. 赵明等[10]通过LBM进行数值模拟, 并采用相图和PSD揭示了开缝圆环腔内热对流遵循Ruelle-Takens-Newhouse的混沌转捩路径. Mun等[11,12]采用浸没边界-有限体积法对多组菱形分布的四热源热对流系统进行了研究, 发现几何排列对冷、热壁面平均努塞尔数(average Nusselt number, $ \overline{Nu} $)的影响较小, 但会诱发流动模式的不对称性. 这些研究共同证实, 受限空间内多热源热对流系统的流动模式的转变与演化路径强烈依赖于系统的具体构型与控制参数.

      然而现有关于受限空间内热对流的非线性动力学的研究, 大多聚焦于几何构型单一或热边界简单的系统, 对于兼具多热源布局与强曲壁约束的复杂结构, 其非线性动力学行为与分岔特性的研究尚存在明显空白[4]. 基于此, 本文采用格子Boltzmann方法, 针对抽象出的二维菱形分布的4个等温热源的圆形封闭腔进行研究, 系统分析其热对流稳定性的非线性动力学演化路径. LBM因其计算效率高、易于并行化实现[13], 且在求解Boussinesq近似的热对流这类弱可压缩问题上具有良好效果, 成为本研究理想的数值工具. 通过固定几何参数, 本文系统刻画了随Ra变化时, 热力学系统从稳态、经分岔演化至混沌的完整路径, 揭示了流动模式转变的非线性动力学机制, 从而为描述该类复杂约束下的多热源热对流系统的传热机理提供理论基础.

    2.   问题描述

      2.1.   物理模型

    • 物理模型如图1所示, 封闭腔体为圆形, 半径为$ {R}_{\text{o}} $, 壁面为恒定温度$ {T}_{\text{c}} $ 的冷源, 壁面速度为0.4个圆形热源半径为$ {R}_{\text{i}} $, 以模型几何中心为对称点呈菱形分布且对角线长度d相等, 取值为$ 1.4({R}_{\text{o}}-{R}_{\text{i}}) $, 壁温恒定$ {T}_{\text{h}} $, 壁面速度为0. 圆腔与热源半径比值为$ {R}_{\text{o}}/{R}_{\text{i}}=5 $. 重力加速度为$ \boldsymbol{g} $, 方向竖直向下. 定义无量纲温度$ \theta =(T-{T}_{\text{h}})/({T}_{\text{h}}-{T}_{\text{c}}) $, 圆心O为监测点.

      • 2.2.   LBM模型

    • 本文采用Guo等[14]提出的基于Boussinesq假设的耦合的格子BGK模型进行求解. 速度场采用D2G9模型, 其离散速度$ \boldsymbol{e}_{\alpha} $和权重系数$ {\omega }_{\alpha } $为标准形式:

      式中, $ c={\delta }_{x}/{\delta }_{\text{t}}=1 $, $ {\delta }_{x} $和$ {\delta }_{\text{t}} $对应网格步长和时间步长. 其中时间步长在格子单位中为1, 后续图形中的时间数据采样间隔为500个时间步.

      基于Boussinesq假设的模型外力项$ {F}_{\alpha } $定义如下:

      式中, $ \delta $为克罗内克尔符号, 有效外力$ \boldsymbol{G}=-\beta (T- {T}_{0})\boldsymbol{g} $, $ \beta $为热胀系数. 故而$ {f}_{\alpha } $的演化方程如下:

      其中, $ {\tau }_{\text{f}} $为无量纲松弛时间. 平衡态分布函数$ f_{\alpha }^{\text{eq}} $为

      其中密度$ \rho $、速度$ \boldsymbol{u} $、压力$ p $可由分布函数的矩计算得到:

      温度场采用D2Q4离散模型求解, 其演化方程为

      平衡态分布函数$ T_{\alpha }^{\text{eq}} $为

      宏观温度由下式计算:

      通过Chapman-Enskog多尺度展开, 可恢复描述热对流的Boussinesq方程组:

      式中, ν为运动黏度系数, χ为热扩散系数.

      热对流问题常采用无量纲参数描述问题, 普朗特数(Prandtl number)定义为$ Pr =\nu/\chi $; Ra定义为$ Ra=g\beta \Delta T{H}^{3} Pr /{\nu }^{2} $; 马赫数(Mach number, Ma)定义为$ Ma={u}_{\text{c}}/{c}_{\text{s}} $. 其中, $ H $为特征长度, $ {u}_{\text{c}} $为特征速度, 定义为

      取Ma = 0.1满足不可压缩, 取Pr = 0.706. 由上述无量纲参数可定义无量纲松弛时间$ {\tau }_{\text{f}} $, $ {\tau }_{\text{T}} $为

      • 2.3.   边界处理及$ \overline{Nu} $计算

    • 采用Guo等[15]提出的结合非平衡外推格式和空间插值的曲面处理方法处理曲面边界, 该方法中的边界插值节点示意如图2所示. 以速度边界为例, 边界节点的分布函数可表示为

      式中, 平衡态分布边界宏观量构造, 非平衡态部分通过流体节点分布函数插值获得. 该格式在时空上均具有二阶精度, 且数值稳定性良好[13]. 具体方式为

      因模型关于坐标轴对称, 将第一象限对应方向插值系数映射至其余3个象限, 避免浮点舍入误差影响物理模型的对称性.

      局部努塞尔数(local Nusselt number, $ N{u}_{\text{L}} $)定义为[16]

      因模型边界为曲面边界, 本文沿壁面内法线方向均匀选取 5个采样点, 通过双线性插值获得温度值并对最小二乘法拟合的径向温度分布进行求导, 可得到壁面温度梯度:

      式中, 正负号依壁面法向方向定义(外壁取负, 内壁取正).

      沿整个圆周均匀选取多个角度计算$ N{u}_{\text{L}} $, 得到局部努塞尔数分布$ N{u}_{\text{L}}(\phi ) $, 再通过积分平均可得$ \overline{Nu} $为

      内壁面$ \overline{Nu} $本质上是各局部热源贡献的算术和, 外壁面$ \overline{Nu} $反映了系统的整体换热效能, 下文给出的努塞尔数均为外壁面$ \overline{Nu} $.

      • 2.4.   模型验证

    • 上述双分布模型的数值方案通过C++编程实现, 并与文献[17]实验结果进行对比验证. 验证算例的几何模型、边界条件及无量纲参数均参照实验参数进行设置. 表1列出了不同Ra下, 计算所得平均表面导热系数(average surface heat transfer coefficient, $ \overline{{K}_{\text{eq}}} $)与文献[17]实验所得$ \overline{{K}_{\text{eq}}} $的对比. 结果显示, 两者在所有对比工况下均吻合良好, 最大相对误差不超过4%, 程序验证可靠.

      本研究采用均匀网格进行计算. 表2网格无关性验证结果表明: 3套不同密度网格所得到的$ \overline{Nu} $最大相对误差均在5%以内, 且逐级加密网格之间的相对误差随网格加密呈现收敛趋势. 在$ Ra= 8\times {10}^{4} $的定常状态下, 各网格结果差异很小; 而在更高的Ra下, 如图3(a), (b)所示, 网格401×401已很好收敛, 网格201×201与其他网格吸引子拓扑结构仍然保持一致性. 综合来看网格201×201虽然演化速度明显快于加密的网格, 但并不改变非线性演化的定性特征. 考虑到需要大量参数扫描, 且本文主要讨论的是流态演化的定性规律, 因此目前的网格精度是可以接受的.

      • 2.5.   非线性动力学分析方法

    • 采用功率谱密度、相图和最大李雅普诺夫指数(largest Lyapunov exponent, LLE)作为系统随Ra变化的非线性行为的分析方法. 选取几何中心点作为监测点, 该点与4个热源等距, 可削弱单一热源的局部边界效应; 且位于竖直中心线及环流分隔线附近, 其水平速度涨落与方向切换对系统环流拓扑变化敏感, 是捕捉非线性演化的理想观测点.

      具体分析方法如下.

      1)相图分析. 取$ \overline{Nu} $进行相空间重构判断系统周期性. 对于混沌态下的非线性分析, 取监测点水平速度$ u $与垂直速度$ v $绘制二维相图. 对于可能发生的异宿分岔, 利用$ u $进行相空间重构, 以判断分岔类型.

      2)功率谱密度分析. 对系统$ \overline{Nu} $的时间序列进行PSD分析, 以识别系统的特征频率及其演化. 采用Cooley-Tukey算法及矩形窗.

      3)最大李雅普诺夫指数计算. 采用Wolf 算法计算LLE. 该算法通过追踪相空间中相邻轨道的指数发散率来量化系统的混沌特性[18]. 由于数值模拟的离散性及有限精度, LLE的计算值难以严格为零, 故以LLE显著大于零作为系统进入混沌的判据. 采用自相关法计算延迟时间; 使用Cao方法确定嵌入维度; 演化时间取5倍相邻周期或拟周期Ra的PSD幅值最大的峰对应的时间. 计算结果见表3.

    3.   结果与讨论

    • 图4和图5整体展示了Ra在$ 8.5\times {10}^{4} $和$ 9\times {10}^{5} $范围内, 温度云图(图4)和流线(图5)的变化过程. 图6为系统中心点u的极值随瑞利数Ra演化的分岔图. 图中清晰展示了系统从稳态到混沌的完整演化路径. 由于系统在超临界叉式分岔后产生两个完全对称的分支, 本文仅展示其中一个分支. 以下各节将结合相图、功率谱及LLE分析, 详细阐述各分岔区间内的动力学特征.

      • 3.1.   稳态解

    • 如图6所示, 在$ Ra=8.5\times {10}^{4} $时, 系统处于定常态. 由图4(a)和图5(a)可知, $ Ra=8.5\times {10}^{4} $时系统温度场与流场均以竖直中心线镜像对称, 形成稳定的双环流结构. 上热源及左右热源与壁面间的狭窄区域呈现出一维线性温度分布, 与文献[19]观察到的“伪扩散”现象物理本质一致, 表明在整体对流主导的系统中, 局部强几何约束仍可诱发扩散主导的传热机制. 由图7也可以看出, 系统$ \overline{Nu} $于稳定后不再发生变化.

      回到图5(b), 可以观测到, 在$ Ra=9\times {10}^{4} $时圆腔顶部壁面竖直中心线附近产生一对相向旋转的亚稳态涡胞, 流场与温度场仍保持对称状态. 随Ra增大, 一个涡发展, 另一个涡耗散湮灭, 最终形成单一涡胞(图5(c)), 等温线向一侧倾斜(图4(c)). 由于边界条件与插值系数均严格对称, 该过程对应一次无外部偏置的超临界叉式分岔, 对称性自发破缺, 最终收敛分支由初始场决定, 在随机扰动下以等概率选择不同分支, 具有路径依赖性.

      • 3.2.   周期解

    • 由图6可见, $ Ra=1.1\times {10}^{5} $时系统已由定常态演化为稳定的周期振荡态. 该振荡表现为上热源主导的涡胞周期性运动: 大涡胞先扩展并与主环流贯通, 随后脱离并回缩. 此时, 图8系统$ \overline{Nu} $周期性振荡, 图9(a)系统PSD图谱的峰构成为$ {f}_{0} $的基频峰以及其整数$ 2{f}_{0} $, $ 3{f}_{0} $等的倍谐波峰, 图10系统$ \overline{Nu} $重构轨迹及图11系统u-v相图轨迹均表现为一个极限环. 且以高Ra终场作为初始场获得的解无滞后现象. 这些特征共同确认, 系统在$ Ra\approx 1.1\times {10}^{5} $附近经历了一个超临界Hopf分岔(supercritical Hopf bifurcation, Hopf分岔)进入单周期振荡.

      观测图5(d)可知, 经Hopf分岔后系统随Ra增大由涡核伸缩逐渐发展为同时向对称破缺方向移动的单侧双涡运动模式. Ra继续增大, 观测图5(e)可见, 流场中左右热源羽流另一侧开始生成一个小涡, 系统流动复杂性增大, 相应的图11内$ u\text-v $相图轨迹的拓扑结构变得更复杂. 且系统振荡周期随Ra的变化呈现非单调性: 在Hopf分岔初期, 周期随Ra增大而缩短, 表现为图8中$ \overline{Nu} $振荡周期缩短及图9(a)中PSD图谱中主频频率增大; 而在形成稳定的单侧涡对运动后, 周期转而随Ra增大而延长, 表现为图8中$ \overline{Nu} $振荡周期延长及图9(b)中PSD图谱主频频率减小. 尽管$ u\text-v $相图拓扑随Ra增大趋于复杂, 但$ \overline{Nu} $重构轨迹拓扑结构仍为单一的极限环, 且PSD图谱始终由基频峰及其谐波峰主导. 这些特征共同表明, 在单侧双涡运动模式的Ra区间内, 系统维持单周期极限环振荡.

      因系统存在一对的镜像对称的稳定解, 即P+及P-, 取$ Ra=1.75\times {10}^{5} $的两个分支进行非线性分析. 可发现其对应的$ \overline{Nu} $时间序列(图8)、PSD图谱(图9(a))、相空间重构的吸引子形态(图10)完全一致, 图中对应两个分支的点和线完全重合. 仅在以$ u=0 $为对称轴的$ u-v $相图(图11)上呈现镜像对称. 这表明, 除空间指向相反外, 两个分支的解具有完全相同的非线性动力学行为. 因此, 在后续分析中仅选取其中一个分支进行深入研究.

      • 3.3.   异宿分岔与混沌转变

    • 当Ra数增至$ 3.15\times {10}^{5} $时, 系统动力学发生本质改变. 由图6得知, 系统对称性开始恢复, 表现为上热源双涡开始随机且间歇地朝另一侧运动, 且随Ra增加频次增加直至近乎相等. 取$ Ra=4\times {10}^{5} $进行分析, $ u\text-v $相图(图12)的轨迹呈现出复杂的交叉结构, 并通过狭窄的通道连接着两个对应失稳前两个分支方向的高密度区域, 这一拓扑结构符合异宿轨道的特征. PSD图谱(图13)由离散谱峰演化为宽频连续谱且LLE = 0.0642, 确认了系统进入混沌状态.

      为更直观揭示其动力学拓扑, 以水平速度 $ u $重构相空间, 如图14所示. 结果表明, 混沌轨迹聚集在4个主要轨迹丛. 轨迹只在相邻轨迹丛间穿梭, 形成一个方形的连接网络, 这直观地印证了混沌运动源于在两个不稳定状态之间的异宿游走机制.

      综上所述, 系统经历了一次异宿分岔(heteroclinic bifurcation). 由于离散热源引入的复杂局部激励与系统全局耗散难以平衡, 导致形成的异宿环失稳, 从而引发了通向混沌的转变. 而文献[20]报道了在单热源对流中, 系统经历对称破缺及Hopf分岔后, 呈现双涡单侧振荡现象, 并最终通过蓝天灾难(Blue-sky catastrophe)步入混沌. 这表明尽管具有相似的物理现象, 系统几何形状及热源分布的复杂性直接决定了通向混沌的具体分岔类型.

      当Ra进一步增至$ 6\times {10}^{5} $, 上热源附近流型由 “双涡运动”向“中心羽流”过渡. 上热源附近两侧运动双涡涡核生成的时间间隔显著增加, 中心热羽流在上热源与顶部壁面之间形成并占据主导, 由图15可见, 系统$ \overline{Nu} $振荡的剧烈期的时间间隔随之增加. 这导致依赖于间歇性双涡运动的、大尺度的竖直方向速度振荡在流场中的整体贡献急剧降低; 系统的动力学能量与复杂性转而高度集中于水平方向的运动.

      该Ra下系统在相空间(图12)中体现为混沌吸引子拓扑结构的“坍缩”: 轨迹在竖直方向$ v $的分布范围大幅紧缩, 整体结构沿水平方向$ u $集中演化. 在此重构后的吸引子中, 密集的核心区域对应于系统长时间处于由上热源中心羽流振荡; 而稀疏的外围轨迹则对应着偶发的、残余的侧向双涡运动. 在$ u $相空间重构图(图14)可见, 早期混沌态中基于异宿游走的动力学特征逐渐消退, 原有的吸引子结构失稳并重组为一个新的吸引子. 此时系统LLE = 0.4274, 仍处于混沌状态.

      • 3.4.   窗口期

    • 上热源流型在Ra约$ 7.5\times {10}^{5} $时发展为完全为中心羽流, 由图15可见$ \overline{Nu} $振荡幅度减小, 图16相空间轨迹所有轨迹收敛永不重复、不发散、也不会收敛于任何一条闭合轨道, 呈现一个二维环面. 此时混沌吸引子失稳, 同时系统经由一次环面分岔(Neimark-Sacker bifurcation)进入拟周期窗口, 此时PSD图谱(图17)有两个不可公约的频率$ {f}_{1} $和$ {f}_{2} $, 由这两个频率的线性组合$ f=m{f}_{1}+n{f}_{2} $组成其他峰的频率, m和n为整数.

      Ra约为$ 8.75\times {10}^{5} $时, 如图18所示, 系统的$ \overline{Nu} $数以高低峰交替振荡且严格重复, 图16相空间轨迹呈现双环状且该双环已持续一小段Ra区间. PSD图谱(图17)拟周期态下存在的两个不可公约的独立频率峰$ {f}_{1} $与$ {f}_{2} $消失, 取而代之的是一个新的基频$ {f}_{0} $及其清晰的$ {f}_{0}/2 $次谐波峰. 故而可以判断此时系统进入锁频周期振荡.

      当Ra增至$ 9.2\times {10}^{5} $时, 从图18可看出系统$ \overline{Nu} $从周期振荡状态失稳, 锁频消失, 再次进入拟周期状态. PSD图谱(图17)两个不可公约的频率$ {f}_{1} $和$ {f}_{2} $再次出现, 对应图16的$ \overline{Nu} $重构轨迹的二维环面结构. Ra进一步增至 $ 9.5\times {10}^{5} $, 系统$ \overline{Nu} $重构轨迹的拟周期环面破裂, 转变为一个具有复杂分形结构的奇怪吸引子. 图19中PSD图谱中离散的线谱被宽频、连续的噪声背景谱所取代且LLE = 0.1574, 系统进入混沌态.

      由图5(g), (h)可见, 下热源羽流两侧的大涡随羽流摆动而上升, 导致左右半腔流体输运的延迟. 由于物理模型水平对称, 在拟周期阶段, 左右热源$ \overline{Nu} $的振荡模式(图20)相似但存在稳定相位差, 对应PSD图谱(图21)主峰频率一致. 同样可从PSD图谱中得出, 上下热源$ \overline{Nu} $的振荡模式不同, 且下热源$ \overline{Nu} $与系统$ \overline{Nu} $的振荡模式接近. 相对于下热源, 上热源$ \overline{Nu} $振荡频率高但振幅低.

      由于上下热源的羽流均在环流分界线上, 两热源羽流振幅及频率的差异均会影响系统左右半腔流体输运的延迟, 从而影响系统的动力学行为. 上热源的快速调制改变了主环流在顶部区域的瞬时偏转, 而下热源的慢速摆动则通过两侧大涡的交替上升驱动左右半腔的周期性输运延迟. 两者频率不可公度, 使得系统整体处于拟周期状态; 而左右热源因受相同的周期性强迫, 其相位差仍保持稳定.

      随着Ra 增大, 当上下热源的频率达到一定比值范围, 系统进入锁频周期, 且经过4个热源的频率协调, 最终呈现为有理数比. 继续增大Ra, 各热源羽流振荡强度持续增强, 导致维持频率锁定的动态平衡被打破, 系统失锁并从周期态中脱离. 随后, 系统动力学迅速失稳, 经由拟周期或环面破裂途径, 最终演变为混沌态.

    4.   结　论

    • 采用格子Boltzmann方法并通过相空间重构、功率谱密度分析及最大李雅普诺夫指数, 研究了随Ra增大时, 具有菱形分布的4个等温内热源的圆形封闭腔内热对流的非线性特征, 刻画出了从稳态经多重分岔通向混沌的复杂动力学路径. 主要结论如下.

      1)低Ra下, 系统呈现稳定的对称双环流定常态.

      2)在$ Ra\approx 9.5\times {10}^{4} $时, 系统经历一次超临界叉式分岔, 表现为上热源单涡流的自发对称破缺, 形成一对稳定的镜像对称定常态解.

      3)在$ Ra\approx 1.1\times {10}^{5} $时, 系统经由超临界Hopf分岔, 从定常态演化为一个稳定的极限环, 对应上热源双涡往复运动以及Ra增大后的单侧运动引发的全局性的单周期振荡. 该极限环的性质随Ra增大而复杂化.

      4)在$ Ra\approx 3.15\times {10}^{5} $时, 系统流动模式转变为上热源双涡双侧运动, 发生全局异宿分岔, 导致极限环失稳, 直接进入混沌态. 此后, 混沌吸引子内部结构持续演化, 上热源附近流型从间歇性双涡运动逐渐重组为持续的中心羽流.

      5)在$ Ra\approx 7.5\times {10}^{5} $附近, 混沌吸引子失稳, 同时系统经由一次环面分岔, 进入一个拟周期窗口. 在$ Ra\approx 8.75\times {10}^{5} $时, 系统发生频率锁定, 退化为锁频周期态. 随着Ra进一步增大, 锁频条件破坏, 系统再次进入拟周期态, 最终$ Ra\approx 9.5\times {10}^{5} $该拟周期环面破裂, 系统再次进入混沌态.

      本研究旨在阐明特定紧凑空间几何构型下的热对流稳定性的非线性演化路径, 但各热源间的相互作用及其对分岔行为的独立影响尚未充分揭示. 后续工作可系统探究热源间距与温度非均匀性等因素, 以更全面地理解具有多热源的热对流系统非线性动力学的复杂特性.

    参考文献 (20)

目录

/

返回文章
返回