搜索

x

留言板

尊敬的读者、作者、审稿人, 关于本刊的投稿、审稿、编辑和出版的任何问题, 您可以本页添加留言。我们将尽快给您答复。谢谢您的支持!

姓名
邮箱
手机号码
标题
留言内容
验证码

基于格子Boltzmann方法的二维气泡群熟化过程模拟

陈效鹏 冯君鹏 胡海豹 杜鹏 王体康

引用本文:
Citation:

基于格子Boltzmann方法的二维气泡群熟化过程模拟

陈效鹏, 冯君鹏, 胡海豹, 杜鹏, 王体康

Lattice Boltzmann method based simulation of two dimensional bubble group ripening process

Chen Xiao-Peng, Feng Jun-Peng, Hu Hai-Bao, Du Peng, Wang Ti-Kang
PDF
HTML
导出引用
  • Ostwald熟化(ripening)是指局部热力学平衡状态下, 颗粒/液滴/气泡系统为了减小界面能而自发地进行颗粒群尺度分布调整的过程, 具有重要研究价值. 针对目前数值模拟研究不充分的现状, 本文采用格子Boltzmann方法, 对相变速率主控的二维蒸气泡系统演化开展了数值模拟研究. 模拟结果与本文推导的二维气泡群演化标度律符合较好, 证实了格子Boltzmann方法对复杂相变-物质输运过程捕捉的准确性. 研究同时表明, 蒸气泡系统演化过程中物质输运为液相压力不平衡所驱动, 并且在小气泡“溃灭”过程中水动力学作用会影响气泡群半径分布函数的局部细节; 气-液状态方程参数对熟化过程的影响效果分析显示, 气液两相比内能差是驱动相变的核心要素, 此差异越大相变速率越快, 该结论进一步诠释了化学势驱动熟化过程的物理图像.
    Ostwald ripening refers to a process of a particle/droplet/bubble system under local thermal equilibrium state adjusting the size distribution spontaneously to reduce the total surface energy. A lattice Boltzmann approach is used to simulate the ripening process of a two dimensional vapor bubble cluster dominated by phase transition kinetics. By comparing the numerical results with the theoretical prediction derived in two-dimensional space, it is shown that the lattice Boltzmann method is accurate in the simulations. The results also indicate that the mass transfer in liquid phase is driven by hydrodynamic pressure distribution and the hydrodynamic collapse of the bubbles influences the size distribution function in a small size region. The influence of the parameters in the equation of state of the material is studied further. A positive relation between phase transition speed and specific internal energy is proposed, which enhances the thermal fundamental of phase transition.
      通信作者: 陈效鹏, xchen76@nwpu.edu.cn ; 胡海豹, huhaibao@nwpu.edu.cn
    • 基金项目: 国家自然科学基金 (批准号: 11872315, 51879218, 52071272)、基础前沿项目(批准号: JCKY2018-18)和陕西省自然科学基础研究计划(批准号: 2020JC-18)资助的课题.
      Corresponding author: Chen Xiao-Peng, xchen76@nwpu.edu.cn ; Hu Hai-Bao, huhaibao@nwpu.edu.cn
    • Funds: Project supported by the National Natural Science Foundation of China (Grant Nos. 11872315, 51879218, 52071272), the Basic Frontier Project, China (Grant No. JCKY2018-18), and the Natural Science Basic Research Program of Shaanxi Province, China (Grant No. 2020JC-18)
    [1]

    Wagner Voorhees P 1985 J. Stat. Phys. 38 231

    [2]

    Bray A J 1994 Adv. Phys. 43 357Google Scholar

    [3]

    Bender H, Ratke L 1998 Acta Mater. 46 1125Google Scholar

    [4]

    Alkemper J, Snyder V A, Akaiwa N, et al. 1999 Phys. Rev. Lett. 82 2725Google Scholar

    [5]

    Diddens C, Tan H, Lv P, et al. 2017 J. Fluid Mech. 823 470Google Scholar

    [6]

    Li Y, Garing C, Benson S M 2020 J. Fluid Mech. 889 889

    [7]

    Ardell A J 1990 Phys. Rev. B 41 2554Google Scholar

    [8]

    Voorhees P W, Glicksman M E 1984 Acta Metall. 32 2001Google Scholar

    [9]

    Fan D, Chen L, Chen S, Voorhees P W 1998 Comput. Mater. Sci. 9 329Google Scholar

    [10]

    Li J, Guo C, Ma Y, Wang Z, Wang J 2015 Acta Mater. 90 10Google Scholar

    [11]

    Wang Y, Li J, Zhang L, Wang Z, Wang J 2020 Model. Simul. Mater. Sci. Eng. 28 075007Google Scholar

    [12]

    Moats K A, Asadi E, Laradji M 2019 Phys. Rev. E 99 012803Google Scholar

    [13]

    Watanabe H, Suzuki M, Inaoka H, Ito N 2014 J. Chem. Phys. 141 234703Google Scholar

    [14]

    Watanabe H, Inaoka H, Ito N 2016 J. Chem. Phys. 145 124707Google Scholar

    [15]

    Aidun C K, Clausen J R 2010 Annu. Rev. Fluid Mech. 42 437

    [16]

    Guo Z, Shu C 2013 Lattice Boltzmann Method and Its Applications in Engineering (Beijing: World Scientific) pp4117–4134

    [17]

    Shan X, Chen H 1993 Phys. Rev. E 47 1815Google Scholar

    [18]

    Yuan P, Schaefer L 2006 Phys. Fluids 18 042101Google Scholar

    [19]

    Huang H, Krafczyk M, Lu X 2011 Phys. Rev. E 84 046710Google Scholar

    [20]

    Chen X, Zhong C, Yuan X 2011 Comput. Math. Appl. 61 3577Google Scholar

    [21]

    Shi Y, Luo K, Chen X, Li D 2020 J. Hydro. Ser. B 32 845Google Scholar

    [22]

    Shan M, Zhu C, Yao C, Cheng Y, Jiang X 2016 Chin. Phys. B 25 104701Google Scholar

    [23]

    Yang Y, Shan M, Han Q, Kan X 2021 Chin. Phys. B 30 024701Google Scholar

    [24]

    Li Q, Kang Q J, Francois M M 2015 Int. J. Heat Mass Transfer 85 787Google Scholar

    [25]

    Shen L Y, Tang G H, Li Q, Shi Y 2019 Langmuir 35 9430Google Scholar

    [26]

    Chang X, Huang, H, Cheng Y, Lu X 2019 Int. J. Heat Mass Transfer 139 588Google Scholar

    [27]

    Liu M, Chen X 2017 Phys. Fluids 29 082102Google Scholar

    [28]

    Lifshitz I M, Slyozov V V 1961 J. Phys. Chem. Solids 19 35

    [29]

    Chai Z, Sun D, Wang H and Shi B 2018 Int. J. Heat Mass Transfer 122 631Google Scholar

    [30]

    Carter A H 2007 Classical and Statistical Thermodynamics (Beijing: Tshinghua University Press) pp19–34

    [31]

    波林 B E, 普劳斯尼茨 J M, 奥康奈尔 J P著 (赵红玲, 王凤坤, 陈圣坤 译) 2001 汽液物性估算手册 (北京: 化学工业出版社) 第493页

    Poling B E, Prausnitz J M, O’Connell (translated by Zhao H L, Wang F K, Chen S K) 2001 The Properties of Gases and Liquids (Beijing: Chemical Industry Publishing House) p493 (in Chinese)

    [32]

    Anderson J 1995 Computational Fluid Dynamics (Berlin: McGraw-Hill Education) pp87–104

  • 图 1  双气泡熟化过程气泡演化 (a) 双气泡熟化过程相分部演化; (b) 大、小气泡半径 ($ R_2, R_1 $)演化, $T=0.80 T_{\rm{c}}, d=500$; (c) 双气泡熟化过程中气相区总面积 ($ S $)变化. 图(b)和图(c)中的红色曲线为滤波结果. 对应相变速率系数 $k_{\rm{s}}=0.6845$

    Fig. 1.  Ostwald ripening for two bubbles: (a) Evolution of the two bubbles; (b) relation of $ R_1, R_2\sim t $ at $T=0.80 T_{\rm{c}},\; d=500$; (c) evolution of the total vapor phase area ($ S $). The red curves are filted results in pannels (a) and (b). Corresponding $k_{\rm{s}}=0.6845$.

    图 2  双气泡熟化过程中密度、压力分布变化 (a) 双气泡中心连线上的密度分布及演化, 箭头标记了压力测点位置; (b) 图(a)中测点压力随时间的演化过程, 空心圈表示数值模拟结果, 曲线表示压力数据的平滑结果. 图标“SB”, “LB”分别表示小气泡和大气泡

    Fig. 2.  The pressure distribution in bubble and liquid phase during ripening: (a) Pressure distribution at different time, with the arrows marking the detected point in pannel (b); (b) temporal pressure on the four marked points. The open symbols denote the simulated results and the curves the smoothed results. Labels “SB” and “LB” denote large and small bubble, respectively

    图 3  双气泡增长速率与相变驱动作用的关系, 实心点表示小气泡, 空心点表示大气泡演化过程. “C1, C2, C3” 数据分别对应计算条件($ R_1/R_2/d $) = $(57/60/500), (57/63/500), $$ (57/60/1000)$. 线性拟合数据点(虚线)得到$ k_{\rm{s}}=0.6845 $. 图中箭头表示熟化过程发展方向

    Fig. 3.  The relation of $ {{\rm{d}}} R/ {{\rm{d}}} t\sim(1/R_{\rm{c}}-1/R) $. The closed symbols represent small bubble evolution, and the open ones the large bubble. “C1, C2” and “C3” correspond to the cases ($ R_1/R_2/d $) = $(57/60/500), (57/63/500), (57/60/ $$ 1000)$, respectively. A linear fitting shows the slope of $ k_{\rm{s}}=0.6845 $, and the arrows show the directions of the ripening processes for large and small bubbles, respectively

    图 4  $ k_{\rm{s}} $对热力学参数的依赖关系 (a)不同温度和表面张力下的气泡增长速率, 其中实心点表示$ (a=1, b=4) $条件结果, 空心点表示表示其他$ (a, b) $对结果 ($a=1.0, 1.05, $$ b=3.75—4.0$); (b) $ W\text{-}k_{\rm{s}} $关系, 其中$ W $表示介质相变所做的机械功

    Fig. 4.  Variation of $ k_{\rm{s}} $ depending on the thermal parameters in ripening process of dual bubble: (a) The $T\text{-}k_{\rm{s}}$ and $\sigma\text{-}k_{\rm{s }}$ relations. The close symbols correspond to $(a=1, $$ b=4)$, and the opens for various $(a, b)=(1.0-1.05, 3.75- $$ 4.0)$. (b) The $W\text{-}k_{\rm{s}}$ relation, where $ W $ denotes mechanical work done in phase transition

    图 5  气泡群熟化过程计算结果 (a) $ T=0.80 T_{\rm{c}} $气泡群熟化过程; (b) 气相面积演化 ($ T=0.80 T_{\rm{c}} $)和不同温度条件下气泡群临界半径增长趋势, $ \psi $全场气相面积占比; (c) 不同温度条件下计算域内气泡数量演化过程. 图中斜线显示了$ N\sim t^{-1} $标度律

    Fig. 5.  The simulated ripening process for vapor bubble cluster: (a) Bubble distribution in the ripening process as $ T=0.8 T_{\rm{c }}$; (b) vapor area ratio ($ \psi $) evolution ($ T=0.8 T_{\rm{c}} $) and $R_{\rm{c}}^2\text{-} t$ relations for different temperatures; (c) bubble number evolutions for different temperatures. The dashed line indicates the $ N\sim t^{-1} $ scaling

    图 6  气泡半径分布函数演化结果

    Fig. 6.  The evolution of $ {\cal{F}}-R $ relation in ripening process

  • [1]

    Wagner Voorhees P 1985 J. Stat. Phys. 38 231

    [2]

    Bray A J 1994 Adv. Phys. 43 357Google Scholar

    [3]

    Bender H, Ratke L 1998 Acta Mater. 46 1125Google Scholar

    [4]

    Alkemper J, Snyder V A, Akaiwa N, et al. 1999 Phys. Rev. Lett. 82 2725Google Scholar

    [5]

    Diddens C, Tan H, Lv P, et al. 2017 J. Fluid Mech. 823 470Google Scholar

    [6]

    Li Y, Garing C, Benson S M 2020 J. Fluid Mech. 889 889

    [7]

    Ardell A J 1990 Phys. Rev. B 41 2554Google Scholar

    [8]

    Voorhees P W, Glicksman M E 1984 Acta Metall. 32 2001Google Scholar

    [9]

    Fan D, Chen L, Chen S, Voorhees P W 1998 Comput. Mater. Sci. 9 329Google Scholar

    [10]

    Li J, Guo C, Ma Y, Wang Z, Wang J 2015 Acta Mater. 90 10Google Scholar

    [11]

    Wang Y, Li J, Zhang L, Wang Z, Wang J 2020 Model. Simul. Mater. Sci. Eng. 28 075007Google Scholar

    [12]

    Moats K A, Asadi E, Laradji M 2019 Phys. Rev. E 99 012803Google Scholar

    [13]

    Watanabe H, Suzuki M, Inaoka H, Ito N 2014 J. Chem. Phys. 141 234703Google Scholar

    [14]

    Watanabe H, Inaoka H, Ito N 2016 J. Chem. Phys. 145 124707Google Scholar

    [15]

    Aidun C K, Clausen J R 2010 Annu. Rev. Fluid Mech. 42 437

    [16]

    Guo Z, Shu C 2013 Lattice Boltzmann Method and Its Applications in Engineering (Beijing: World Scientific) pp4117–4134

    [17]

    Shan X, Chen H 1993 Phys. Rev. E 47 1815Google Scholar

    [18]

    Yuan P, Schaefer L 2006 Phys. Fluids 18 042101Google Scholar

    [19]

    Huang H, Krafczyk M, Lu X 2011 Phys. Rev. E 84 046710Google Scholar

    [20]

    Chen X, Zhong C, Yuan X 2011 Comput. Math. Appl. 61 3577Google Scholar

    [21]

    Shi Y, Luo K, Chen X, Li D 2020 J. Hydro. Ser. B 32 845Google Scholar

    [22]

    Shan M, Zhu C, Yao C, Cheng Y, Jiang X 2016 Chin. Phys. B 25 104701Google Scholar

    [23]

    Yang Y, Shan M, Han Q, Kan X 2021 Chin. Phys. B 30 024701Google Scholar

    [24]

    Li Q, Kang Q J, Francois M M 2015 Int. J. Heat Mass Transfer 85 787Google Scholar

    [25]

    Shen L Y, Tang G H, Li Q, Shi Y 2019 Langmuir 35 9430Google Scholar

    [26]

    Chang X, Huang, H, Cheng Y, Lu X 2019 Int. J. Heat Mass Transfer 139 588Google Scholar

    [27]

    Liu M, Chen X 2017 Phys. Fluids 29 082102Google Scholar

    [28]

    Lifshitz I M, Slyozov V V 1961 J. Phys. Chem. Solids 19 35

    [29]

    Chai Z, Sun D, Wang H and Shi B 2018 Int. J. Heat Mass Transfer 122 631Google Scholar

    [30]

    Carter A H 2007 Classical and Statistical Thermodynamics (Beijing: Tshinghua University Press) pp19–34

    [31]

    波林 B E, 普劳斯尼茨 J M, 奥康奈尔 J P著 (赵红玲, 王凤坤, 陈圣坤 译) 2001 汽液物性估算手册 (北京: 化学工业出版社) 第493页

    Poling B E, Prausnitz J M, O’Connell (translated by Zhao H L, Wang F K, Chen S K) 2001 The Properties of Gases and Liquids (Beijing: Chemical Industry Publishing House) p493 (in Chinese)

    [32]

    Anderson J 1995 Computational Fluid Dynamics (Berlin: McGraw-Hill Education) pp87–104

  • [1] 陈效鹏, 冯君鹏, 胡海豹, 杜鹏, 王体康. 基于格子Boltzmann方法的二维汽泡群熟化过程模拟. 物理学报, 2022, (): . doi: 10.7498/aps.71.20212183
    [2] 王伟, 揭泉林. 基于机器学习J1-J2反铁磁海森伯自旋链相变点的识别方法. 物理学报, 2021, 70(23): 230701. doi: 10.7498/aps.70.20210711
    [3] 周光雨, 陈力, 张鸿雁, 崔海航. 基于格子Boltzmann方法的自驱动Janus颗粒扩散泳力. 物理学报, 2017, 66(8): 084703. doi: 10.7498/aps.66.084703
    [4] 王佐, 张家忠, 王恒. 非正交多松弛系数轴对称热格子Boltzmann方法. 物理学报, 2017, 66(4): 044701. doi: 10.7498/aps.66.044701
    [5] 张娅, 潘光, 黄桥高. 疏水表面减阻的格子Boltzmann方法数值模拟. 物理学报, 2015, 64(18): 184702. doi: 10.7498/aps.64.184702
    [6] 刘邱祖, 寇子明, 贾月梅, 吴娟, 韩振南, 张倩倩. 改性疏水固壁润湿性反转现象的格子Boltzmann方法模拟. 物理学报, 2014, 63(10): 104701. doi: 10.7498/aps.63.104701
    [7] 黄桥高, 潘光, 宋保维. 疏水表面滑移流动及减阻特性的格子Boltzmann方法模拟. 物理学报, 2014, 63(5): 054701. doi: 10.7498/aps.63.054701
    [8] 解文军, 滕鹏飞. 声悬浮过程的格子Boltzmann方法研究. 物理学报, 2014, 63(16): 164301. doi: 10.7498/aps.63.164301
    [9] 史冬岩, 王志凯, 张阿漫. 任意复杂流-固边界的格子Boltzmann处理方法. 物理学报, 2014, 63(7): 074703. doi: 10.7498/aps.63.074703
    [10] 任晟, 张家忠, 张亚苗, 卫丁. 零质量射流激励下诱发液体相变及其格子Boltzmann方法模拟. 物理学报, 2014, 63(2): 024702. doi: 10.7498/aps.63.024702
    [11] 刘邱祖, 寇子明, 韩振南, 高贵军. 基于格子Boltzmann方法的液滴沿固壁铺展动态过程模拟. 物理学报, 2013, 62(23): 234701. doi: 10.7498/aps.62.234701
    [12] 郭亚丽, 徐鹤函, 沈胜强, 魏兰. 利用格子Boltzmann方法模拟矩形腔内纳米流体Raleigh-Benard对流 . 物理学报, 2013, 62(14): 144704. doi: 10.7498/aps.62.144704
    [13] 曾建邦, 李隆键, 蒋方明. 气泡成核过程的格子Boltzmann方法模拟. 物理学报, 2013, 62(17): 176401. doi: 10.7498/aps.62.176401
    [14] 张晋鲁, 李玉强, 赵兴宇, 黄以能. 用Weiss分子场理论对有外电场时铁电体系相变特征的研究. 物理学报, 2012, 61(14): 140501. doi: 10.7498/aps.61.140501
    [15] 曾建邦, 李隆键, 廖全, 蒋方明. 池沸腾中气泡生长过程的格子Boltzmann方法模拟. 物理学报, 2011, 60(6): 066401. doi: 10.7498/aps.60.066401
    [16] 宋婷婷, 何捷, 林理彬, 陈军. 氧化钒晶体的半导体至金属相变的理论研究. 物理学报, 2010, 59(9): 6480-6486. doi: 10.7498/aps.59.6480
    [17] 曾建邦, 李隆键, 廖全, 陈清华, 崔文智, 潘良明. 格子Boltzmann方法在相变过程中的应用. 物理学报, 2010, 59(1): 178-185. doi: 10.7498/aps.59.178
    [18] 卢玉华, 詹杰民. 三维方腔温盐双扩散的格子Boltzmann方法数值模拟. 物理学报, 2006, 55(9): 4774-4782. doi: 10.7498/aps.55.4774
    [19] 刘 红, 王 慧. 双轴性向列相液晶的相变理论. 物理学报, 2005, 54(3): 1306-1312. doi: 10.7498/aps.54.1306
    [20] 李华兵, 黄乒花, 刘慕仁, 孔令江. 用格子Boltzmann方法模拟MKDV方程. 物理学报, 2001, 50(5): 837-840. doi: 10.7498/aps.50.837
计量
  • 文章访问数:  518
  • PDF下载量:  21
  • 被引次数: 0
出版历程
  • 收稿日期:  2021-11-26
  • 修回日期:  2022-02-17
  • 上网日期:  2022-05-24
  • 刊出日期:  2022-06-05

基于格子Boltzmann方法的二维气泡群熟化过程模拟

    基金项目: 国家自然科学基金 (批准号: 11872315, 51879218, 52071272)、基础前沿项目(批准号: JCKY2018-18)和陕西省自然科学基础研究计划(批准号: 2020JC-18)资助的课题.

摘要: Ostwald熟化(ripening)是指局部热力学平衡状态下, 颗粒/液滴/气泡系统为了减小界面能而自发地进行颗粒群尺度分布调整的过程, 具有重要研究价值. 针对目前数值模拟研究不充分的现状, 本文采用格子Boltzmann方法, 对相变速率主控的二维蒸气泡系统演化开展了数值模拟研究. 模拟结果与本文推导的二维气泡群演化标度律符合较好, 证实了格子Boltzmann方法对复杂相变-物质输运过程捕捉的准确性. 研究同时表明, 蒸气泡系统演化过程中物质输运为液相压力不平衡所驱动, 并且在小气泡“溃灭”过程中水动力学作用会影响气泡群半径分布函数的局部细节; 气-液状态方程参数对熟化过程的影响效果分析显示, 气液两相比内能差是驱动相变的核心要素, 此差异越大相变速率越快, 该结论进一步诠释了化学势驱动熟化过程的物理图像.

English Abstract

    • Ostwald熟化现象是自然和工业界广泛存在的现象之一, 主要指在热力学局部平衡状态下颗粒、液滴或者气泡群体系为了进一步降低系统的表面能, 而自发地对气泡尺度分布进行调整; 发生大颗粒不断增长, 小颗粒减小直至消失的现象[1,2]. 生活中长期在冷冻环境下储存的冰激凌的口感会变硬、长期静置的自来水中会出现小气泡, 这些过程均包含了熟化机理. 工业生产中也采用熟化机理解释和控制金属凝固过程中固态晶粒的长大[3,4]、Ouzo鸡尾酒的乳化过程[5], 以及地下油藏中$ {\rm{CO}}_2 $的聚集[6]等.

      对Ostwald熟化现象的理论研究始于Lifshitz和Slyozov提出的扩散机制主控的控制方程, 并进一步给出了渐近解. 几乎同时, Wagner[1]提出了相变主控的数学描述, 简称为LSW理论体系. 此后, 针对该理论体系中要求颗粒分散度大、溶质在连续相中分布均匀等假设的局限性, 众多学者提出更新的模型, 使理论成果能被进一步应用于实际问题. Ardell[7]提出特殊半径法, 对颗粒发展的相互影响展开了讨论; Voorhees和Glicksman[8]采用点源和点汇模型刻画颗粒对环境浓度的影响, 理论探讨了颗粒体积有限条件下的熟化过程. 这些研究是对LSW理论的重要促进. 与理论同步发展的是实验研究. Bender和Ratke[3]对熔融Cu-Co合金中富Co固体颗粒的熟化行为进行了测试, 其中Co颗粒体积分数在25%—70%之间; Alkemper等[4]则实验观测了熔融Pb-Sn合金中的富Sn固体颗粒的长大过程, 并将结果与多种理论模型进行了对比. 实验研究一方面证实了理论模型, 并为理论体系的完善提供了方向; 另一方面也有力地支撑了生产技术的发展. 近几十年, 数值模拟被逐渐应用于Ostwald熟化研究. Fan等[9]采用相场方法对熟化过程进行了模拟, 所获得的颗粒形貌与实验观测相仿, 并且颗粒群形貌特征参数演化的标度律关系与理论结果一致. Li等[10]和Wang等[11]数值研究了初始半径分布对圆形颗粒的熟化过程的影响, 获得了颗粒增长的标度律关系, 同时也探讨了数值结果与理论预测的偏差问题. Moats等[12]采用晶粒相场模型, 开展了液-固体系的熟化过程研究, 他们除了验证了经典标度律关系, 同时基于数值结果分析了体系熵的演化. Watanab等[13,14]采用分子动力学方法对蒸气泡群的熟化过程开展了研究, 发现随着温度的提升熟化主控因素由相变速率转变为扩散速度. 目前, 数值模拟研究正成为构建理论体系的重要一环.

      格子Boltzmann方法 (lattice Boltzmann method, LBM)是基于格子气自动机方法发展起来的一种流体动力数值方法, 其本质是对Boltzmann方程的一种离散求解方法[15,16]. 通过进一步与湍流模型、热传导方程、电磁学方程等耦合, LBM已经被成功应用于多种复杂流动现象的模拟上. LBM多相流模型及计算是该领域比较活跃的一支, 其中以Shan-Chen模型[17]最为引人注目, 且历久弥新. 通过引入分子间作用力模型, LBM-Shan-Chen模型成功实现了多相流体的相分离模拟, 且在密度比、表面张力等关键参数上, 模拟结果与理论预测一致, 关键是该模型反映了两相分离的物理本质. 通过与经典状态方程相耦合, Yuan和Schaefer[18]进一步推进了LBM与现有热力学体系的融合, 同时也扩大了其能捕捉的多相系统密度比, 提高了数值稳定性. Huang等[19]对LBM多相流计算过程中, 体积力的施加方法、热力学参数影响的细节开展了对比研究, 为后期的数值参数的选取提供了依据. 前期研究中, 人们致力于LBM模拟能力的提升及采用LBM对典型多相流现象开展机理研究[15].

      相比而言, 采用LBM进行相变机理以及流动现象的研究偏少. Chen跟合作者[20,21]采用LBM对空化气泡、气泡群的演化过程进行了模拟研究. 他们发现, LBM对空泡动力学的模拟与Rayleigh-Plesset方程预测一致. Shan及合作者[22,23]研究了壁面润湿性、热效应对空泡溃灭的影响. 对于相变因素更为重要的过程, Li等[24]利用LBM模拟了沸腾过程, 准确捕捉到了流动模态的转换过程; Shen等[25]则研究了壁面浸润性对凝结过程的影响; Chang等[26]利用LBM开展了壁面结构对沸腾影响以及强化传热效果的研究, 探讨了其中流-热耦合的力学机理. 前述研究和分析中尽管涉及相变过程, 但多采用了相变速率无限快这一假设. 对于一些特殊情况, 如Ostwald熟化、气泡高速振荡等, 该假设无疑会引起一定误差.

      本文采用LBM-Shan-Chen多相流模型, 对蒸气泡群相变速率主控的Ostwald熟化过程开展了模拟研究. 探讨了不同物理参数对典型气泡群体系演化过程及相变速率的影响. 结合LSW熟化理论, 验证了数值方法的准确度; 数值结果同时也表明, 蒸气泡系统的演化过程中水动力学因素有明显的影响, 这使得熟化过程具有一些有别于经典熟化过程的特点.

    • LBM是近年来迅速发展起来的一种流体动力学数值模拟方法. 因其具有明确的微观粒子动理学基础, 可以对较为复杂的多相流动现象进行模拟[27]. 本文采用Bhatnagar-Gross-Krook (BGK)单松弛Shan-Chen多相流模型开展模拟[17]. LBM控制方程为

      $ \begin{split} &f_\alpha({\boldsymbol{x}}+{\boldsymbol{e}}_\alpha \text{δ} t, t+\text{δ} t)-f_\alpha({\boldsymbol{x}},t)\\ =\;&-{1\over\tau}\left[f_\alpha({\boldsymbol{x}},t)-f^{{\rm{eq}}}_\alpha(\rho,{\boldsymbol{u}})\right]+S_\alpha({\boldsymbol{x}},t), \end{split} $

      式中, $ f_\alpha $表示时刻t、坐标x处、具有微观速度$ {\boldsymbol{e}}_\alpha $的粒子分数, 即密度分布函数. 流体的局部密度表示为 $\rho({\boldsymbol{x}}, t)=\displaystyle\sum\nolimits_\alpha f_\alpha$, 局部动量表示为 $\rho{\boldsymbol{u}}= \displaystyle\sum\nolimits_\alpha{\boldsymbol{e}}_\alpha f_\alpha$. 本文采用LBM-D2Q9模型[16], 认为二维空间中粒子具有9种可能的离散速度, $\alpha=0, \ 1, \ \cdots, \ 8$. (1)式左端表示颗粒的迁移运动: 在一个时间步$ \text{δ} t\; (=1) $中, 颗粒从x跳跃到$ {\boldsymbol{x}}+{\boldsymbol{e}}_\alpha \text{δ} t $. D2Q9模型中格点上的离散速度表示为

      $ {\boldsymbol{e}}_\alpha=\left\{ \begin{aligned} & {{0}},&(\alpha=0),\qquad\quad\\ &[\cos(\pi (\alpha-1)/2), \sin(\pi (\alpha-1)/2)],&(\alpha=1,2,3,4),\;\\ &\sqrt{2}[\cos(\pi (2\alpha-1)/4), \sin(\pi(2\alpha-1)/4)],&(\alpha=5,6,7,8). \end{aligned} \right. $

      (1)式中等号左边第一项表示碰撞作用. 在准平衡态条件下, 颗粒的碰撞导致其在速度相空间的分布趋于Maxwell分布, 该项中τ表示碰撞的弛豫时间. 其中$ f_\alpha^{{\rm{eq}}} $为Maxwell分布函数在相空间中的离散形式,

      $ \begin{equation} f_\alpha^{{\rm{eq}}}=\rho\omega_\alpha\left[1+{{\boldsymbol{e}}_\alpha\cdot{\boldsymbol{u}}\over c_{\rm{s}}^2}+{({\boldsymbol{e}}_\alpha\cdot{\boldsymbol{u}})^2\over 2c_{\rm{s}}^4}-{{\boldsymbol{u}}^2\over 2c_{\rm{s}}^2}\right], \end{equation} $

      式中$ c_{\rm{s}} $是声速 (D2Q9模型中的$ c_{\rm{s}}=1/\sqrt{3} $); $ \omega_\alpha $是加权系数: $ \omega_0=4/9, \ \omega_{1-4}=1/9, \ \omega_{5-8}=1/36 $.

      LBM控制方程中, $ S_\alpha $表征了外力作用对粒子分布函数的影响, 可以是长程的体积力作用, 也可以是短程的分子间作用. 本文采用Kupershtokh作用力模型[19]获得作用力在微观粒子速度方向上的“投影”,

      $ \begin{array}{l} S_\alpha=f_\alpha^{{\rm{eq}}}(\rho,{\boldsymbol{u}}+\Delta{\boldsymbol{u}})-f_\alpha^{{\rm{eq}}}(\rho,{\boldsymbol{u}}),\quad \Delta{\boldsymbol{u}}={\boldsymbol{F}}\text{δ} t/\rho, \end{array} $

      其中$ {\boldsymbol{F}} $表示外力矢量. 相应的流体局部宏观速度表示为${\boldsymbol{u}}= \Big(\displaystyle \sum f_\alpha{\boldsymbol{e}}_\alpha \Big)/\rho+{\boldsymbol{F}}\text{δ} /(2\rho)$. Huang等[19]对比了Guo作用力模型和Kupershtokh作用力模型, 发现相对于Guo等提出的作用力模拟, Kupershtokh作用力模型具有密度捕捉准确、数值稳定性好的特点. Shan和Chen通过与密度相关的势函数ψ来构造分子间的不平衡作用力:

      $ {\boldsymbol{F}}({\boldsymbol{x}},t)=-G\psi({\boldsymbol{x}},t)\sum\limits_\alpha\omega_\alpha\psi({\boldsymbol{x}}+{\boldsymbol{e}}_\alpha\text{δ} t,t){\boldsymbol{e}}_\alpha. $

      给定介质的状态方程, 当地势函数可以表达为

      $ \psi=\sqrt{2(p-\rho c_{\rm{s}}^2)\over Gc_{\rm{s}}^2}, $

      其中G为常数. 本文采用Carnahan-Starling状态方程确定$ p(\rho) $关系:

      $ p(\rho)=\rho{\cal{R}}T{1+b\rho/4+(b\rho/4)^2-(b\rho/4)^3\over(1-b\rho/4)^3}-a\rho^2, $

      其中, $ a=0.4963({\cal{R}}T_{\rm{c}})^2/p_{\rm{c}} $, $ b=0.18727{\cal{R}}T_{\rm{c}}/p_{\rm{c}} $. 本文重点取$ a=1, \ b=4, \ {\cal{R}}=1 $开展计算, 由此得到临界密度和临界温度: $\rho_{\rm{c}}=0.1136, \ T_{\rm{c}}= 0.0943$. 将$ S_\alpha $代入LBM控制方程, 可以实现非理想气体流动的模拟, 并且可以捕捉到介质的相变过程[20,22,24].

    • 1958年Lifshitz和Slyozov[28]对扩散控制的熟化进行了机理研究, 他们首先数学描述了析出颗粒群在溶液中的生长过程: 引入描述粒子半径分布的函数$ {\cal{F}}(R, t) $ (R为粒子的半径, t为时间), 并建立了演化方程. Wagner进一步对反应主控与扩散主控的熟化过程进行了研究, 从而建立了描述颗粒群熟化过程的LSW理论体系[1,2]. 该理论体系由动理学方程、连续性方程和质量守恒方程组成. 动理学方程刻画了单个颗粒的增长规律, 其包含“溶质”在连续相中的迁移与析出过程的数学描述为

      $ \begin{equation} \nabla^2C=0,\quad { {{\rm{d}}} R\over {{\rm{d}}} t}\equiv\dot{R}=k_{\rm{s}}\left({1\over R_{\rm{c}}}-{1\over R}\right), \end{equation} $

      其中, $ k_{\rm{s}} $表征了相变速率, $ R_{\rm{c}} $为颗粒临界半径, C为溶液中溶质浓度. (4)式显示Ostwald熟化过程中, 当析出颗粒尺度小于临界尺度时, 会逐渐消融, 反之则长大. $ R_{\rm{c}} $与析出相在环境流体中的溶解浓度、析出相表面张力能量有关. 临界尺度可以通过下面的公式计算获得:

      $ \begin{equation} R_{\rm{c}}(t)={\displaystyle \int_0^\infty R{\cal{F}} {{\rm{d}}} R\over \displaystyle \int_0^\infty {\cal{F}} {{\rm{d}}} R}. \end{equation} $

      连续性方程描述了粒子半径分布函数的演化规律:

      $ \begin{equation} {\partial{\cal{F}}\over\partial t}+{\partial\over\partial R}\left({\cal{F}}\dot{R}\right)=0, \end{equation} $

      即具有半径R的颗粒数($ {\cal{F}} {{\rm{d}}} R $)随时间的变化率, 取决于气泡的增长: $ \left[\left({\cal{F}} {{\rm{d}}} R/ {{\rm{d}}} t\right)_{R^+}-\left({\cal{F}} {{\rm{d}}} R/ {{\rm{d}}} t\right)_{R^-}\right] $. 该表达式假设颗粒在熟化过程中没有发生合并或分裂现象. 质量守恒方程规定了析出相的总体积不变, 即

      $ \begin{equation} \int_0^\infty R^2{\cal{F}}(R,t) {{\rm{d}}} R={\rm const}. \end{equation} $

      (4)式—(7)式中, 假设析出相颗粒是半径为R的圆形 (二维).

      本文针对反应控制熟化过程开展研究, 即“溶质”在连续相中扩散得足够快, (4)式中相变过程成为制约颗粒生长的主要因素, 扩散方程可以忽略. Lifshitz和Pitaevskii[28]在数学上推导了$ t\rightarrow \infty $时扩散主控熟化过程中, 颗粒尺寸等参数增长的标度律. 参考Watanabe等[14]采用的标度律推导过程, 可以得到二维条件下的蒸气泡群几何特征参数演化规律 (附录A):

      $ \begin{array}{l} R_{\rm{c}}\propto t^{1/2}, \end{array} $

      $ \begin{array}{l} N\propto t^{-1}. \end{array} $

      同时, 分布函数演化满足相似关系:

      $ \begin{equation} {\cal{F}}=h(t)g(z),\quad z=R/R_{\rm{c}}, \end{equation} $

      其中

      $ g(z)=\left\{ {\begin{aligned} &{4z\over 3}\left(2\over2-z\right)^5\exp{\left(-3z\over 2-z\right)},& 0 < z < 2,\;\\& 0,&z > 2. \qquad\end{aligned}} \right. $

      熟化控制方程的完整解则需要通过数值求解方式获得.

    • 一般而言, Ostwald熟化过程是一种近平衡态的热力学过程, 不考虑水动力学因素. 数值计算中初始场需要尽量满足这些条件. 在双气泡熟化过程中, 首先根据设计条件对单气泡情况进行计算, 测量获得计算稳定后的气泡的半径$ R_1, \; R_2 $及相应的气液密度 $ \rho_{1 {\rm{L}}}, \; \rho_{1 {\rm{V}}} $$ \rho_{2 {\rm{L}}}, \; \rho_{2{\rm{ V}}} $. 在进行双气泡熟化过程初始化时, 取初始气泡半径为$ R_1, \; R_2 $, 气、液两相密度分别为$(\rho_{1 {\rm{V}}}+\rho_{2 {\rm{V}}})/2, \; (\rho_{1 {\rm{L}}}+ \rho_{2{\rm{ L}}})/2$. 对双气泡初始条件的准确设定, 可以减小初始化误差, 从而使有效数据的时间区间延长. 对于多气泡测试, 本文参考(11)式对气泡半径分布进行设定. 此时由于气泡数目较多, 直接根据介质状态方程结果设定初始气、液密度. 为了避免在气泡熟化过程中发生融合、破坏理论分析设定的前提条件, 气泡间距须保证足够大.

      本文针对二维空间的多气泡熟化过程开展LBM模拟研究, 计算域取为矩形, 针对不同的气泡数和气泡大小, 取适当的计算域大小: 双气泡计算域大小取$ 1000\times 500 $, 多气泡取$ 3600\times 3600 $, 气泡群初始数量约为800. 计算域四周边界取周期边界. 在后处理过程中, 本文采用图像处理技术 (如二值化、边界提取等), 对计算结果图片进行颗粒几何特征提取, 可以获得气泡数、各个气泡的直径和位置等信息.

    • 本节针对双气泡、多气泡熟化过程开展模拟和分析. 前者在反映熟化特征的同时, 便于我们对相关细节开展测试, 进而对 (模拟的)熟化过程物质输运细节进行分析; 后者对标一般熟化过程, 期望探讨蒸气泡群熟化过程的基本规律.

    • 本部分取$ \tau=1.0 $ (该参数决定了流体的黏性, 本文暂不讨论其影响), $ R_1=57, \; R_2=60 $, 气泡间距$ d=500 $, $ T=0.8 T_{\rm{c}} $为例, 给出双气泡熟化过程的计算结果. 该温度下, 气、液两相初始密度分别为 $ \rho_{\rm{V}}=0.01853, \; \rho_L{\rm{}}=0.3064 $, 气-液表面张力系数$\sigma= 0.0075$. 前期测算表明, 上述结果准确地再现了介质的热力学状态[20].

      图1(a)给出了双气泡熟化中大小气泡的发展, 显示了小气泡消失、大气泡增长的过程. 进一步统计气泡半径可以发现, 两个气泡的大小在发展过程中是振荡的, 这导致了气相区域总面积 (S )的振荡 (图1(b),(c)). 原因可能是在进行双气泡计算时, 密度、速度场初始化仍然不够精确, 引起了不平衡 (水动力)压力波在计算域内来回传播, 从而激励起气泡的非物理振荡. 当滤除掉气相总面积的波动成分后, 可以看出气相总面积是稳定的, 且熟化过程最终形成的大气泡面积十分接近预设气相面积, 这满足LSW理论的质量守恒约束条件. 该守恒律与水动力学作用占优的空化气泡发展过程有明显差别[20]. 将气相总面积波动滤除以后, 可以获得更为光顺的气泡半径演化过程 (见图1(b)中的红线).

      图  1  双气泡熟化过程气泡演化 (a) 双气泡熟化过程相分部演化; (b) 大、小气泡半径 ($ R_2, R_1 $)演化, $T=0.80 T_{\rm{c}}, d=500$; (c) 双气泡熟化过程中气相区总面积 ($ S $)变化. 图(b)和图(c)中的红色曲线为滤波结果. 对应相变速率系数 $k_{\rm{s}}=0.6845$

      Figure 1.  Ostwald ripening for two bubbles: (a) Evolution of the two bubbles; (b) relation of $ R_1, R_2\sim t $ at $T=0.80 T_{\rm{c}},\; d=500$; (c) evolution of the total vapor phase area ($ S $). The red curves are filted results in pannels (a) and (b). Corresponding $k_{\rm{s}}=0.6845$.

      另一方面, 提取熟化过程中不同时刻, 过两个气泡中心的轴线上的密度分布演化结果(图2(a)). 可以看出, 在熟化过程中气泡内的密度基本保持不变, 而在液相区压力有明显波动. 因此可以通过Laplace关系获得液相压力分布: 大气泡周围的液相压力较高 ($ p_{\rm{L}}=p_{\rm{B}}-\sigma/R_2 $, $ p_{\rm{B}} $表示气泡内的压力), 小气泡附近的压力 ($ p_{\rm{S}}=p_{\rm{B}}-\sigma/R_1 $)较低. 图2(b)给出了小气泡、大气泡中心, 以及计算域轴线上液相区两个点 (图2(a)中以箭头标记)上的压力演化时序. 图2(b)中, $ t < 2000 $时, 液相压力降低, 反映了精确密度场的建立过程, 同时由于初始阶段$ R_1\approx R_2 $, 因而液相中两个测点的压力差距几乎分辨不出来; 随着气泡半径差距拉大, 液相压力梯度逐渐明显起来. 这样的压力差异在量级上与Laplace公式的预测也是符合的: $(p_{\rm{L}}-p_{\rm{S}})_{t=1400}\approx 3\times10^{-4}$. 正是液相中的压力梯度, 驱动了介质往小气泡一侧移动, 从而使原来的小气泡区域逐渐被液相物质所填充. 这样的物质迁移机制替代了经典的Ostwald熟化理论模型中由离散相浓度梯度驱动的物质迁移机制. 从本文结果看, Watanabe等[14]采用传统的颗粒扩散假设对蒸气泡群熟化过程进行描述是有待商榷的.

      图  2  双气泡熟化过程中密度、压力分布变化 (a) 双气泡中心连线上的密度分布及演化, 箭头标记了压力测点位置; (b) 图(a)中测点压力随时间的演化过程, 空心圈表示数值模拟结果, 曲线表示压力数据的平滑结果. 图标“SB”, “LB”分别表示小气泡和大气泡

      Figure 2.  The pressure distribution in bubble and liquid phase during ripening: (a) Pressure distribution at different time, with the arrows marking the detected point in pannel (b); (b) temporal pressure on the four marked points. The open symbols denote the simulated results and the curves the smoothed results. Labels “SB” and “LB” denote large and small bubble, respectively

      进而基于平滑过的气泡演化数据对熟化过程进行定量分析. 在双气泡条件下, (4)式和(5)式变得比较简单:

      $ \left\{\begin{aligned} &\dot{R}_i=k_{\rm{s}}\left({1/ R_{\rm{c}}}-{1/ R_i}\right),\; i=1,\; 2\\ &R_{\rm{c}}=(R_1+R_2)/2. \end{aligned} \right. $

      可以直接通过数值方法求解. 图3给出了大、小气泡增长速率与相变驱动作用的关系, 曲线斜率即为$ k_{\rm{s}} $. 结果表明当前模拟中气泡增长过程服从LSW理论预测的规律. 进一步改变$ R_1, \; R_2 $d, 发现上述动力学关系仍然得到满足, 且在状态方程、温度(T )一定的情况下$ k_{\rm{s}} $值不变. 由此可知, $ k_{\rm{s}} $可以由介质的状态方程确定. 另外, 数据显示, 在熟化过程末期, 数据点存在一定程度波动, 这可能会影响$ k_{\rm{s}} $的拟合结果 (对比图3插图). 这有可能是由于小气泡在熟化末期, 界面运动速度过快, 从而引起了水动力学效应 (与空泡溃灭相似). 考虑到图3中熟化末期数据点斜率基本与拟合直线斜率相近, 以及对于不同算例气泡半径演化波动区较难判定, 下文仍然采用该统计方法.

      图  3  双气泡增长速率与相变驱动作用的关系, 实心点表示小气泡, 空心点表示大气泡演化过程. “C1, C2, C3” 数据分别对应计算条件($ R_1/R_2/d $) = $(57/60/500), (57/63/500), $$ (57/60/1000)$. 线性拟合数据点(虚线)得到$ k_{\rm{s}}=0.6845 $. 图中箭头表示熟化过程发展方向

      Figure 3.  The relation of $ {{\rm{d}}} R/ {{\rm{d}}} t\sim(1/R_{\rm{c}}-1/R) $. The closed symbols represent small bubble evolution, and the open ones the large bubble. “C1, C2” and “C3” correspond to the cases ($ R_1/R_2/d $) = $(57/60/500), (57/63/500), (57/60/ $$ 1000)$, respectively. A linear fitting shows the slope of $ k_{\rm{s}}=0.6845 $, and the arrows show the directions of the ripening processes for large and small bubbles, respectively

      值得一提的是, Chai等[29]结合相场方法与LBM对双气泡演化过程进行了模拟, 并基于相场方程推导获得了气泡演化的解析解. 这对本问题具有重要的参考价值, 然而也可以看出数学解在实际工程中的应用仍然具有障碍.

    • 本节围绕$ k_{\rm{s }}$与介质状态方程的关系展开讨论. 文献 [19]在探讨状态方程参数对LBM多相流模拟稳定性影响规律时发现, $ a $会改变气液界面的厚度; b参数可以影响热平衡状态下的气液密度, 但不改变气液密度比; 而$ {\cal{R}} $对气液状态没有影响. 经典热力学理论通常只考虑平衡状态, 本节研究从数值角度给出相变速率 (非平衡过程参数)与状态方程参数之间的依赖关系, 同时也可以为面向实际问题的参数设定提供一定的参考.

      图 4给出了不同温度下的$ k_{\rm{s}} $值, 其中每个$ k_{\rm{s}} $点综合了3种不同初始构型下的双气泡演化过程 (参见 3.1节)获得. 结果显示, 尽管对固定$ (a, \; b) $对, $ k_{\rm{s}} $总体随温度T呈降低趋势 (见$(a, \; b)=(1.0, \; 4.0)$数据, 实心点), 但下降的非线性特征明显. 并且综合所有数据, $ k_{\rm{s}} $对温度T的依赖性并不统一, 其同时也决定于其他状态方程参数.

      图  4  $ k_{\rm{s}} $对热力学参数的依赖关系 (a)不同温度和表面张力下的气泡增长速率, 其中实心点表示$ (a=1, b=4) $条件结果, 空心点表示表示其他$ (a, b) $对结果 ($a=1.0, 1.05, $$ b=3.75—4.0$); (b) $ W\text{-}k_{\rm{s}} $关系, 其中$ W $表示介质相变所做的机械功

      Figure 4.  Variation of $ k_{\rm{s}} $ depending on the thermal parameters in ripening process of dual bubble: (a) The $T\text{-}k_{\rm{s}}$ and $\sigma\text{-}k_{\rm{s }}$ relations. The close symbols correspond to $(a=1, $$ b=4)$, and the opens for various $(a, b)=(1.0-1.05, 3.75- $$ 4.0)$. (b) The $W\text{-}k_{\rm{s}}$ relation, where $ W $ denotes mechanical work done in phase transition

      经典熟化模型基于两相界面张力 (σ)解释了Ostwald熟化的机理[1,14], 其可以进一步引申为化学势驱动的物质输运机制. 对于单组份两相等温过程而言, 化学势可以表达为压力的函数 $ \mu(p) $[30]. Watanabe等[14]给出了气泡界面上液体气化过程的化学势变化, 相变速率具有包含σ的前置因子. 然而他们在随后的标度律分析中, 忽略了前置因子的作用, 因而热力学参数对于相变速率的影响并未深入讨论. 本研究数值模拟结果(图4(a))表明, 单纯的相变速率对表面张力的依赖关系同样也不能统一地刻画相变速率的变化. 需要补充说明: 前期测算表明, 对于固定$ (a, \; b) $对, $ T\text{-}\sigma $关系具有很好的线性度, 这与实验观测和相关理论预测是相符的[31].

      从前述温度越高相变速率越低的现象以及化学势具有能量的本质出发, 可以设想: 随着温度趋近于临界温度、两相之间密度差的减小, 驱动相变所需的机械功 ($ p {{\rm{d}}} V $)也会相应减小. 我们基于状态方程理论解提取了单位质量介质相变的功值,

      $ W=p_{\rm{V}}\left({1\over\rho_{\rm{V}}}-{1\over\rho_{\rm{L}}}\right). $

      图5(b)给出了不同热力学参数条件下$ k_{\rm{s}}\text{-} W $关系, 具有较好归一化特性, 并且也与线性变化趋势较为接近 (其中的偏离可能是由于气、液密度同时受表面张力影响所致). 由此推断, 等温条件下共存两相的比内能差可以被看成是驱动相变的根本原因, 相变速率与其有正相关关系.

      图  5  气泡群熟化过程计算结果 (a) $ T=0.80 T_{\rm{c}} $气泡群熟化过程; (b) 气相面积演化 ($ T=0.80 T_{\rm{c}} $)和不同温度条件下气泡群临界半径增长趋势, $ \psi $全场气相面积占比; (c) 不同温度条件下计算域内气泡数量演化过程. 图中斜线显示了$ N\sim t^{-1} $标度律

      Figure 5.  The simulated ripening process for vapor bubble cluster: (a) Bubble distribution in the ripening process as $ T=0.8 T_{\rm{c }}$; (b) vapor area ratio ($ \psi $) evolution ($ T=0.8 T_{\rm{c}} $) and $R_{\rm{c}}^2\text{-} t$ relations for different temperatures; (c) bubble number evolutions for different temperatures. The dashed line indicates the $ N\sim t^{-1} $ scaling

    • 基于前面两节的讨论, 进而对多气泡熟化过程开展分析. 如前文所述, 本节对密度场做简单初始化, 重点参考(11)式对气泡初始半径分布进行设定: 在半径$ 0—32 $范围内取32个分度, 每个分度中按照半径分布函数 (11)式分配气泡数量 (总数800), 这些气泡随机在计算域内分布. 为了与理论预测结果相对比, 计算域需要设得足够大, 以保证气泡群在熟化过程中不发生碰撞/合并. 典型的多气泡熟化过程如图5(a)所示. 在气泡群熟化过程中, 气泡数量逐渐减少(小气泡消失)并伴随着留存气泡尺度的增加, 相应气泡的临界 (平均)半径亦会增大.

      图5(b)给出了相应气泡群$ R_{\rm{c}} $和计算域内气相所占面积的比例 (ψ)随时间的变化, 可见Shan-Chen-LBM多相流模型在相变速率主控的熟化模拟中可以很好地保持气相面积的守恒性 ((7)式), 同时由于多泡的平均效应, 此条件下压力波造成的气泡总面积波动也被抑制了. 气泡总面积的守恒保证了气泡群临界半径演化 (数值结果)的$ t^{1/2} $标度律关系 ((8)式), 以及气泡数量的$ t^{-1} $标度律 ((9))式、图5(c)). 图5(b)中, 各组$ R_{\rm{c}}^2\sim t $的数据点对应的斜率与双气泡条件下的气泡增长斜率相匹配, 这说明多泡条件下的气泡增长与双泡时的机理是一致的. 针对当前计算参数, 本文同时也调节了计算域面积、气泡间距、气泡排布等几何条件, 所获标度律结果与展示结果一致 (误差小于5%).

      如前文所述, 在LSW理论分析基于$ {\cal{F}} $的时间相似性 ((11)式), 图6(a)给出了采用有限差分方法对联立方程组 (4)式—(6)式进行数值求解所获得的结果. 其中, 空间微分项$ \partial(\dot{R}{\cal{F}})/\partial R $采用迎风格式进行离散, 时间推进采用显示格式[32]. 方程中$ k_{\rm{s}} $值采用了图5(b)中测量获得的$ R_{\rm{c}}^2 $数据点斜率值:

      图  6  气泡半径分布函数演化结果

      Figure 6.  The evolution of $ {\cal{F}}-R $ relation in ripening process

      $ \begin{split} &{{\cal{F}}_i^{n+1}-{\cal{F}}_i^{n}\over\Delta t}=-k_{\rm{s}}{{\cal{F}}_i^n\over R_i^2}-k_{\rm{s}}\left({1\over R_i^n}-{1\over R_{\rm{c}}^n}\right)\cdot{{\cal{F}}_i^n-{\cal{F}}_{i-1}^n\over\Delta R},\, R_i^n > R_{\rm{c}}^n,\\ &{{\cal{F}}_i^{n+1}-{\cal{F}}_i^{n}\over\Delta t}=-k_{\rm{s}}{{\cal{F}}_i^n\over R_i^2}-k_{\rm{s}}\left({1\over R_i^n}-{1\over R_{\rm{c}}^n}\right)\cdot{{\cal{F}}_{i+1}^n-{\cal{F}}_i^n\over\Delta R},\, R_i^n\leqslant R_{\rm{c}}^n, \end{split} $

      其中, 离散化变量$ {\cal{F}}_i^n\equiv {\cal{F}}(R_i, n\Delta t) $, $ \Delta t, \, \Delta R $表示时间和“空间”离散分度. 与LBM结果对比发现, 两者所预测的气泡半径概率密度分布函数在各个时刻均符合较好. 由此说明, LBM方法不仅对气泡熟化过程的总体趋势进行了很好的再现, 同时也能对气泡发展过程中的细节加以捕捉, 这为进一步开展复杂条件下的气泡熟化过程模拟奠定了技术基础.

      进一步仔细对比理论预测曲线与计算结果可以发现: 计算结果在小尺度区域存在$ {\cal{F}} $上翘的趋势, 并且在气泡极小区域存在数据点缺失的现象. 这主要是由于小气泡收缩的最后阶段表面张力作用极大 ($ \propto\sigma/R $), 导致了水动力学意义上的“空泡溃灭”现象, 使小气泡过快消失. 综合整个R数据段的$ {\cal{F}} $分布结果可以知道, 小气泡的水动力学溃灭作用对于大半径区气泡影响较小, 这可能是空泡溃灭过程带来的压力波动在液相中传播. 可以想象当压力波波长大于某个气泡的尺度时, 会明显影响相变过程; 对于气泡而言这样的影响会比较小. 事实上, 在计算结果中确实发现, 在熟化过程中, 一些小气泡在收缩过程中存在抖动现象. 气泡溃灭的水动力学效应特征是蒸气泡群熟化与颗粒或者乳状液熟化过程的不同之处.

    • 本文采用格子Boltzmann方法对二维条件下的蒸气泡群Ostwald熟化现象开展数值研究. Shan-Chen多相流模型以及Carnahan-Starling状态方程用以捕捉真实介质的相变和相分离过程. 本文针对二维条件重新推导了基于Lifshitz-Slyozov-Wagner理论体系的气泡系统特征参数的标度律关系, 并针对相变速率主控的双气泡和多气泡系统熟化过程开展了模拟和分析. 研究结果如下.

      理论分析和数值模拟同时证实, 二维气泡群的反应速率主控的熟化过程中, 气泡半径增长率受气泡间化学势差异驱动, 服从$\dot{R}=k_{\rm{s}}(1/R_{\rm{c}}-1/R)$规律, 由此进一步导致气泡群统计参数服从$R_{\rm{c}}\sim t^{1/2}, \; N\sim t^{-1}$标度律. 该结果证实了格子Boltzmann方法对复杂多相流-相变过程捕捉的适用性, 相对于传统流体力学数值方法, 具有先天的优越性.

      对于双气泡系统气、液两相密度场演化的仔细观测表明, 在熟化过程中液相密度分布具有时变波动特征, 并且大气泡附近的压力高于小气泡周围的压力. 这样的压力分布一方面折射了Laplace压力定律的内涵, 同时也显示液相内的物质输运由该压力不均匀性驱动, 而非传统熟化理论中所提的浓度梯度驱动的输运机理.

      区别于经典熟化理论中忽略了界面移动的惯性效应, 本文结果显示小气泡收缩的最后阶段, 水动力学意义的空化气泡“溃灭”作用会加速小气泡的消失, 从而造成气泡群半径分布函数的小半径分支偏离理论预测值, 并且在细小气泡区造成函数曲线的截断现象.

      本研究基于格子Boltzmann方法重点研究了相变速率主控的蒸气泡群的熟化过程. 相对于基于热力学平衡态假设发展出的理论、数值方法而言, 本方法具有更明确分子动力学基础, 可以作为微观粒子运动与宏观物理现象的桥梁. 本研究为进一步发展格子Boltzmann数值方法、揭示相变机理, 乃至解决实际的工程问题奠定了坚实基础. 在本研究基础上, 可以进一步开展受限空间的颗粒群熟化现象、温度影响、流动影响, 以及扩散主控熟化现象等方面的研究.

      感谢西北工业大学材料学院刘峰教授、中国科学技术大学近代力学系黄海波教授、南方科技大学航空航天工程系单肖文教授给予本工作的关注和讨论.

    • 针对连续性方$程 $

      $ \begin{equation} {\partial{\cal{F}}\over\partial t}+{\partial\over\partial R}\left({\cal{F}}{{\rm d}R\over {{\rm{d}}} t}\right)=0, \end{equation}\tag{A1} $

      基于颗粒熟化过程中$ {\cal{F}} $的相似性, 引入标度率关系:

      $ R_{\rm{c}}\sim t^\alpha,\ {\cal{F}}(R,t)\sim t^\beta\tilde{f}(\tilde{R}),\ \dot{R}/R_{\rm{c}}\sim t^\gamma\tilde{\dot{R}}(\tilde{R}), $

      其中$ \tilde{R}=R/R_{\rm{c}} $ ($ R_{\rm{c}} $为临界/平均半径值). 将(A1)式代入 (二维)析出相体积守恒方程,

      $ \int R^2{\cal{F}} {{\rm{d}}} R=t^{\beta+3\alpha}\int\tilde{R}^2\tilde{f} {{\rm{d}}}\tilde{R}={\rm{const}}, $

      可以得到$ \beta=-3\alpha $. 同时也可以获得体系中颗粒总数满足

      $ N\equiv\int{\cal{F}} {{\rm{d}}} R=t^{-2\alpha}\int\tilde{f} {{\rm{d}}}\tilde{R}\sim t^{-2\alpha}. $

      另外, 将前述标度律关系代回到连续性方程 (A1), 配平t的指数, 可以得到$ \gamma=-1 $, 因此只剩α未知.

      在反应主控熟化过程中, 乳状液体系中小液滴的长大/缩小的速率已经由前人讨论过. 针对蒸气泡问题, Watanabe等[14]根据Gibbs-Thomson公式, 对单位面积相变质量通量进行了估算:

      $ \begin{equation} \dot{R}\propto(\rho_R^{{\rm{eq}}}-\rho_R),\ \rho_R^{eq}=\rho_\infty\left(1-\lambda/R\right), \end{equation} \tag{A2}$

      其中, $ \rho_R $表示气泡内蒸气的密度, $ \rho_\infty $表示平直液面条件下的气体平衡密度, $ \rho_R^{{\rm{eq}}} $表示半径为R的气泡界面上的平衡气体密度, λ为毛细尺度, 它由表面张力、温度、气体分子摩尔体积决定. 由此可得

      $ \dot{R}\propto{1\over R}\left[{R\over R_{\rm{c}}}-1\right],\ R_{\rm{c}}=\left(\rho_\infty\lambda\over\rho_\infty-\rho_R\right), $

      由于$\dot{R}/R_{\rm{c}}\sim\tilde{\dot{R}}t^{-1}$, 故

      $ {\dot{R}\over R_{\rm{c}}}\;\;\sim\;\;{1\over RR_{\rm{c}}}\tilde{R}, \tag{A3}$

      $ \displaystyle {\dot{R}\over R_{\rm{c}}}\;\;\sim \;\; {1\over R_{\rm{c}}^2}{1\over\tilde{R}}\tilde{R}\; \left(\sim t^{-2\alpha}\right) , \tag{A4} $

      ${\dot{R}\over R_{\rm{c}}}\;\; \sim \;\; t^{-1}, \tag{A5}$

      $ \alpha=1/2 $. 更进一步$ R_{\rm{c}}\sim t^{1/2} $, $ N\sim t^{-1} $ (注意, 二维情况下的标度律与三维有差别).

参考文献 (32)

目录

    /

    返回文章
    返回