搜索

文章查询

x

留言板

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

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

非规则形状介质内辐射-导热耦合传热的间断有限元求解

王存海 郑树 张欣欣

非规则形状介质内辐射-导热耦合传热的间断有限元求解

王存海, 郑树, 张欣欣
PDF
HTML
导出引用
导出核心图
  • 采用间断有限元法(discontinuous finite element method, DFEM)求解非规则形状介质内的辐射导热耦合传热问题, 得到了典型非规则形状介质内辐射导热耦合传热问题的高精度数值结果. 和传统连续型有限元方法不同, DFEM将计算区域划分成相互独立的离散单元, 形函数的构造、未知量的加权近似以及控制方程的求解均在每一个离散单元上进行. 通过在单元之间施加迎风格式的数值通量, DFEM保证了整个计算区域的连续性, 因此这种方法兼具良好的几何灵活性和局部守恒性. 推导了辐射传输方程和能量扩散方程的DFEM离散格式, 验证了DFEM求解辐射导热耦合传热问题的正确性; 同时研究了不同几何形状介质内辐射导热耦合传热问题, 得到了典型非规则形状介质内辐射导热耦合传热的高精度数值结果.
      通信作者: 王存海, wangcunhai@ustb.edu.cn ; 郑树, shuzheng@ncepu.edu.cn
    • 基金项目: 国家自然科学基金(批准号: 51906014, 51890891)、中国博士后科学基金(批准号: 2018M641196)和中央高校基本科研业务费专项资金(批准号: FRF-TP-18-072A1)资助的课题
    [1]

    Viskanta R 1965 J. Heat Transf. 87 143

    [2]

    Viskanta R, Incropera F P 1985 J. Sol. Energy Eng. 107 29

    [3]

    Wang P Y, Tan H P, Liu L H, Tong T W 2000 J. Thermophys. Heat Transf. 14 512

    [4]

    Zhang J Q, Nie L R, Chen C Y, Zhang X Y 2016 AIP Adv. 6 075212

    [5]

    Chen R Y, Pan W L, Zhang J Q, Nie L R 2016 Chaos 26 093113

    [6]

    Chen R Y, Tong L M, Nie L R, Wang C J, Pan W L 2017 Physica A 468 532

    [7]

    李树, 李刚, 田东风, 邓力 2013 物理学报 62 249501

    Li S, Li G, Tian D F, Deng L 2013 Acta Phys. Sin. 62 249501

    [8]

    Sun B, Wang H, Sun X B, Hong J, Zhang Y J 2012 Chin. Phys. B 21 129501

    [9]

    Chen R Y, Nie L R, Chen C Y, Wang C J 2017 J. Stat. Mech.- Theory Exp. 013201

    [10]

    Chen R Y, Nie L R, Chen C Y 2018 Chaos 28 053115

    [11]

    梁子长, 金亚秋 2003 物理学报 52 1319

    Liang Z C, Jin Y Q 2003 Acta Phys. Sin. 52 1319

    [12]

    Ben X, Yi H L, Tan H P 2014 Chin. Phys. B 23 099501

    [13]

    赵军明 2007 博士学位论文 (哈尔滨: 哈尔滨工业大学)

    Zhao J M 2007 Ph. D. Dissertation (Harbin: Harbin Institute of Technology) (in Chinese)

    [14]

    赵军明, 刘林华 2007 化工学报 58 1110

    Zhao J M, Liu L H 2007 J. Chem. Ind. Eng. 58 1110

    [15]

    胡帅, 高太长, 刘磊, 易红亮, 贲勋 2015 物理学报 64 094201

    Hu S, Gao T C, Liu L, Yi H L, Ben X 2015 Acta Phys. Sin. 64 094201

    [16]

    Wang C H, Yi H L, Tan H P 2017 J. Quant. Spectrosc. Radiat. Transf. 189 383

    [17]

    高效伟, 王静, 崔苗 2011 中国科学: 物理学 力学 天文学 41 302

    Gao X W, Wang J, Cui M 2011 Sci. Sin. Phys. Mech. Astron. 41 302

    [18]

    孙杰 2016 博士学位论文 (哈尔滨: 哈尔滨工业大学)

    Sun J 2016 Ph. D. Dissertation (Harbin: Harbin Institute of Technology) (in Chinese)

    [19]

    Wang C H, Feng Y Y, Yue K, Zhang X X 2019 Int. Commun. Heat Mass Transf. 108 104287

    [20]

    Sun Y J, Zheng S Jiang B, Tang J C, Liu F S 2019 Int. J. Heat Mass Transf. 145 118777

    [21]

    Tan J Y, Liu L H, Li B X 2006 Numer. Heat Transf. Part B-Fundam. 49 179

    [22]

    Wang C H, Qu L, Zhang Y, Yi H L 2018 J. Quant. Spectrosc. Radiat. Transf. 208 108

    [23]

    Liu L H, Tan J Y, Li B X 2006 J. Quant. Spectrosc. Radiat. Transf. 101 237

    [24]

    Wang C H, Feng Y Y, Ben X, Yue K, Zhang X X 2019 Opt. Express 27 A981

    [25]

    Sun S C, Wang G J, Chen H, Zhang D Q 2019 Int. J. Heat Mass Transf. 134 574

    [26]

    Zheng S, Yang Y, Zhou H 2019 Int. J. Heat Mass Transf. 129 1232

    [27]

    张克瑾, 刘磊, 曾庆伟, 高太长, 胡帅, 陈鸣 2019 物理学报 68 194207

    Zhang K J, Liu L, Zeng Q W, Gao T C, Hu S, Chen M 2019 Acta Phys. Sin. 68 194207

    [28]

    Mishra S C, Krishna C H, Kim M Y 2011 Numer. Heat Transf. Part A-Appl. 60 254

    [29]

    Mishra S C, Roy H K 2007 J. Comput. Phys. 223 89

    [30]

    Howell J R, Menguc M P 2018 J. Quant. Spectrosc. Radiat. Transf. 221 253

    [31]

    Zabihi M, Lari K, Amiri H 2017 J. Braz. Soc. Mech. Sci. Eng. 39 2847

    [32]

    Hesthaven J S, Warburton T 2007 Nodal Discontinuous Galerkin Methods: Algorithms, Analysis, and Applications (New York: Springer Science & Business Media)

    [33]

    Cui X, Li B Q 2005 J. Quant. Spectrosc. Radiat. Transf. 96 383

    [34]

    Liu L H, Liu L J 2007 J. Quant. Spectrosc. Radiat. Transf. 105 377

    [35]

    王存海, 易红亮, 谈和平 2017 工程热物理学报 38 833

    Wang C H, Yi H L, Tan H P 2017 J. Eng. Thermophys. 38 833

    [36]

    Wang C H, Yi H L, Tan H P 2017 Appl. Opt. 56 1861

    [37]

    Wang C H, Yi H L, Tan H P 2017 Opt. Express 25 14621

    [38]

    Feng Y Y, Wang C H 2018 Int. J. Heat Mass Transf. 126 783

    [39]

    Mishra S C, Talukdar P, Trimis D, Durst F 2003 Int. J. Heat Mass Transf. 46 3083

    [40]

    Sun Y, Zhang X 2018 Int. J. Heat Mass Transf. 121 1039

  • 图 1  (a)空间网格; (b)相邻单元间数值通量示意图

    Fig. 1.  (a) Spatial mesh; (b) sketch of numerical flux across the adjacent elements.

    图 2  方形介质对称线x/L = 0.5上无量纲温度T/Tb分布 (a)不同空间网格划分; (b)不同角度划分

    Fig. 2.  Dimensionless temperature T/Tb along the symmetry line x/L = 0.5 for the cases: (a) Different spatial discretization schemes; (b) different angular discretization schemes.

    图 3  方形介质对称线x/L = 0.5上无量纲温度的DFEM结果和DTM结果对比 (a)不同普朗克数; (b)不同散射反照率

    Fig. 3.  Comparison of dimensionless temperature along the square medium symmetry line x/L = 0.5 obtained by DFEM and DTM for the cases: (a) Different Planck numbers; (b) different scattering albedos.

    图 4  不同数值方法得到的方形介质对称线上的无量纲温度对比

    Fig. 4.  Comparison of dimensionless temperature along the square medium symmetry line x/L = 0.5 obtained by different numerical methods.

    图 5  内含圆形热壁面的半圆形介质 (a)结构示意图; (b)网格划分

    Fig. 5.  Semicircle medium with an inner circle hot boundary: (a) Geometry sketch; (b) spatial discretization.

    图 6  不同数值方法得到的半圆介质对称线上温度分布 (a) y = [0.0, 0.2]; (b) y = [0.6, 1.0]

    Fig. 6.  Temperature distributions along the symmetric line of the semicircle medium obtained by different numerical algorithms: (a) y = [0.0, 0.2]; (b) y = [0.6, 1.0].

    图 7  不同普朗克数条件下半圆介质底边上总热流密度分布 (a) Npl = 0.1, (b) Npl = 1.0

    Fig. 7.  Total flux distributions along the bottom boundary of the semicircle medium under the situations with different Plank numbers: (a) Npl = 0.1; (b) Npl = 1.0.

    图 8  内含圆形热边界的非规则形状介质 (a)结构示意图; (b)网格划分

    Fig. 8.  Irregular medium with an inner hot boundary: (a) Geometry sketch; (b) spatial discretization.

    图 9  (a)普朗克数Npl = 0.1和1.0时内含圆形热边界的非规则形状介质中线上温度分布; (b) Npl = 0.1时介质温度分布; (c) Npl = 1.0时介质温度分布

    Fig. 9.  (a) Temperature distributions along the centerline of the irregular medium with an inner hot boundary; (b) temperature distribution within the computation domain for the case of Npl = 0.1; (c) temperature distribution within the computation domain for the case of Npl = 1.0.

    图 10  内含两个圆形热边界的矩形介质 (a)结构示意图; (b)网格划分

    Fig. 10.  The medium the square medium with two circular hot boundaries: (a) Geometry sketch; (b) spatial discr etization.

    图 11  (a)普朗克数Npl = 0.1和1.0时内含两个圆形热边界的矩形介质中线上温度分布; (b) Npl = 0.1时介质温度分布; (c) Npl = 1.0时介质温度分布

    Fig. 11.  (a) Temperature distributions along the centerline of the square medium with two circular hot boundaries; (b) temperature distribution within the computation domain for the case of Npl = 0.1; (c) temperature distribution within the computation domain for the case of Npl = 1.0

  • [1]

    Viskanta R 1965 J. Heat Transf. 87 143

    [2]

    Viskanta R, Incropera F P 1985 J. Sol. Energy Eng. 107 29

    [3]

    Wang P Y, Tan H P, Liu L H, Tong T W 2000 J. Thermophys. Heat Transf. 14 512

    [4]

    Zhang J Q, Nie L R, Chen C Y, Zhang X Y 2016 AIP Adv. 6 075212

    [5]

    Chen R Y, Pan W L, Zhang J Q, Nie L R 2016 Chaos 26 093113

    [6]

    Chen R Y, Tong L M, Nie L R, Wang C J, Pan W L 2017 Physica A 468 532

    [7]

    李树, 李刚, 田东风, 邓力 2013 物理学报 62 249501

    Li S, Li G, Tian D F, Deng L 2013 Acta Phys. Sin. 62 249501

    [8]

    Sun B, Wang H, Sun X B, Hong J, Zhang Y J 2012 Chin. Phys. B 21 129501

    [9]

    Chen R Y, Nie L R, Chen C Y, Wang C J 2017 J. Stat. Mech.- Theory Exp. 013201

    [10]

    Chen R Y, Nie L R, Chen C Y 2018 Chaos 28 053115

    [11]

    梁子长, 金亚秋 2003 物理学报 52 1319

    Liang Z C, Jin Y Q 2003 Acta Phys. Sin. 52 1319

    [12]

    Ben X, Yi H L, Tan H P 2014 Chin. Phys. B 23 099501

    [13]

    赵军明 2007 博士学位论文 (哈尔滨: 哈尔滨工业大学)

    Zhao J M 2007 Ph. D. Dissertation (Harbin: Harbin Institute of Technology) (in Chinese)

    [14]

    赵军明, 刘林华 2007 化工学报 58 1110

    Zhao J M, Liu L H 2007 J. Chem. Ind. Eng. 58 1110

    [15]

    胡帅, 高太长, 刘磊, 易红亮, 贲勋 2015 物理学报 64 094201

    Hu S, Gao T C, Liu L, Yi H L, Ben X 2015 Acta Phys. Sin. 64 094201

    [16]

    Wang C H, Yi H L, Tan H P 2017 J. Quant. Spectrosc. Radiat. Transf. 189 383

    [17]

    高效伟, 王静, 崔苗 2011 中国科学: 物理学 力学 天文学 41 302

    Gao X W, Wang J, Cui M 2011 Sci. Sin. Phys. Mech. Astron. 41 302

    [18]

    孙杰 2016 博士学位论文 (哈尔滨: 哈尔滨工业大学)

    Sun J 2016 Ph. D. Dissertation (Harbin: Harbin Institute of Technology) (in Chinese)

    [19]

    Wang C H, Feng Y Y, Yue K, Zhang X X 2019 Int. Commun. Heat Mass Transf. 108 104287

    [20]

    Sun Y J, Zheng S Jiang B, Tang J C, Liu F S 2019 Int. J. Heat Mass Transf. 145 118777

    [21]

    Tan J Y, Liu L H, Li B X 2006 Numer. Heat Transf. Part B-Fundam. 49 179

    [22]

    Wang C H, Qu L, Zhang Y, Yi H L 2018 J. Quant. Spectrosc. Radiat. Transf. 208 108

    [23]

    Liu L H, Tan J Y, Li B X 2006 J. Quant. Spectrosc. Radiat. Transf. 101 237

    [24]

    Wang C H, Feng Y Y, Ben X, Yue K, Zhang X X 2019 Opt. Express 27 A981

    [25]

    Sun S C, Wang G J, Chen H, Zhang D Q 2019 Int. J. Heat Mass Transf. 134 574

    [26]

    Zheng S, Yang Y, Zhou H 2019 Int. J. Heat Mass Transf. 129 1232

    [27]

    张克瑾, 刘磊, 曾庆伟, 高太长, 胡帅, 陈鸣 2019 物理学报 68 194207

    Zhang K J, Liu L, Zeng Q W, Gao T C, Hu S, Chen M 2019 Acta Phys. Sin. 68 194207

    [28]

    Mishra S C, Krishna C H, Kim M Y 2011 Numer. Heat Transf. Part A-Appl. 60 254

    [29]

    Mishra S C, Roy H K 2007 J. Comput. Phys. 223 89

    [30]

    Howell J R, Menguc M P 2018 J. Quant. Spectrosc. Radiat. Transf. 221 253

    [31]

    Zabihi M, Lari K, Amiri H 2017 J. Braz. Soc. Mech. Sci. Eng. 39 2847

    [32]

    Hesthaven J S, Warburton T 2007 Nodal Discontinuous Galerkin Methods: Algorithms, Analysis, and Applications (New York: Springer Science & Business Media)

    [33]

    Cui X, Li B Q 2005 J. Quant. Spectrosc. Radiat. Transf. 96 383

    [34]

    Liu L H, Liu L J 2007 J. Quant. Spectrosc. Radiat. Transf. 105 377

    [35]

    王存海, 易红亮, 谈和平 2017 工程热物理学报 38 833

    Wang C H, Yi H L, Tan H P 2017 J. Eng. Thermophys. 38 833

    [36]

    Wang C H, Yi H L, Tan H P 2017 Appl. Opt. 56 1861

    [37]

    Wang C H, Yi H L, Tan H P 2017 Opt. Express 25 14621

    [38]

    Feng Y Y, Wang C H 2018 Int. J. Heat Mass Transf. 126 783

    [39]

    Mishra S C, Talukdar P, Trimis D, Durst F 2003 Int. J. Heat Mass Transf. 46 3083

    [40]

    Sun Y, Zhang X 2018 Int. J. Heat Mass Transf. 121 1039

  • [1] 周瑜, 操礼阳, 马晓萍, 邓丽丽, 辛煜. 脉冲射频容性耦合氩等离子体的发射探针诊断. 物理学报, 2020, (): . doi: 10.7498/aps.69.20191864
    [2] 朱存远, 李朝刚, 方泉, 汪茂胜, 彭雪城, 黄万霞. 用久期微绕理论将弹簧振子模型退化为耦合模理论. 物理学报, 2020, (): . doi: 10.7498/aps.69.20191505
    [3] 王琳, 魏来, 王正汹. 垂直磁重联平面的驱动流对磁岛链影响的模拟. 物理学报, 2020, 69(5): 059401. doi: 10.7498/aps.69.20191612
    [4] 蒋涛, 任金莲, 蒋戎戎, 陆伟刚. 基于局部加密纯无网格法非线性Cahn-Hilliard方程的模拟. 物理学报, 2020, (): . doi: 10.7498/aps.69.20191829
    [5] 王培良. 蚁群元胞优化模型在路径规划中的应用. 物理学报, 2020, (): . doi: 10.7498/aps.69.20191774
    [6] 任县利, 张伟伟, 伍晓勇, 吴璐, 王月霞. 高熵合金短程有序现象的预测及其对结构的电子、磁性、力学性质的影响. 物理学报, 2020, 69(4): 046102. doi: 10.7498/aps.69.20191671
    [7] 张雅男, 詹楠, 邓玲玲, 陈淑芬. 利用银纳米立方增强效率的多层溶液加工白光有机发光二极管. 物理学报, 2020, 69(4): 047801. doi: 10.7498/aps.69.20191526
    [8] 王晓雷, 赵洁惠, 李淼, 姜光科, 胡晓雪, 张楠, 翟宏琛, 刘伟伟. 基于人工表面等离激元的厚度渐变镀银条带探针实现太赫兹波的紧聚焦和场增强. 物理学报, 2020, 69(5): 054201. doi: 10.7498/aps.69.20191531
  • 引用本文:
    Citation:
计量
  • 文章访问数:  58
  • PDF下载量:  2
  • 被引次数: 0
出版历程
  • 收稿日期:  2019-08-03
  • 修回日期:  2019-10-25
  • 刊出日期:  2020-02-01

非规则形状介质内辐射-导热耦合传热的间断有限元求解

  • 1. 北京科技大学能源与环境工程学院, 北京 100083
  • 2. 冶金工业节能减排北京市重点实验室, 北京 100083
  • 3. 华北电力大学能源动力与机械工程学院, 北京 102206
  • 通信作者: 王存海, wangcunhai@ustb.edu.cn ; 郑树, shuzheng@ncepu.edu.cn
    基金项目: 国家自然科学基金(批准号: 51906014, 51890891)、中国博士后科学基金(批准号: 2018M641196)和中央高校基本科研业务费专项资金(批准号: FRF-TP-18-072A1)资助的课题

摘要: 采用间断有限元法(discontinuous finite element method, DFEM)求解非规则形状介质内的辐射导热耦合传热问题, 得到了典型非规则形状介质内辐射导热耦合传热问题的高精度数值结果. 和传统连续型有限元方法不同, DFEM将计算区域划分成相互独立的离散单元, 形函数的构造、未知量的加权近似以及控制方程的求解均在每一个离散单元上进行. 通过在单元之间施加迎风格式的数值通量, DFEM保证了整个计算区域的连续性, 因此这种方法兼具良好的几何灵活性和局部守恒性. 推导了辐射传输方程和能量扩散方程的DFEM离散格式, 验证了DFEM求解辐射导热耦合传热问题的正确性; 同时研究了不同几何形状介质内辐射导热耦合传热问题, 得到了典型非规则形状介质内辐射导热耦合传热的高精度数值结果.

English Abstract

    • 参与性介质内辐射导热耦合传热过程广泛存在于众多工程实际应用中[1-3], 如高温熔岩融化、材料加工等, 因此辐射导热耦合问题的研究对相关传热分析和工程设计具有重要的意义. 辐射传输方程是积分-微分型方程, 这类方程多通过数值解法进行求解, 如蒙特卡罗法[4-8]和修正的欧拉算法[9,10]. 由于辐射传输问题常常涉及复杂边界条件和非规则计算区域[11,12], 国内外学者发展了多种数值方法[13-27]用以研究辐射传热问题. Mishra等[28,29]采用不同数值方法分别求解辐射传输方程和能量扩散方程, 发展了求解辐射导热耦合问题的混合方法.

      虽然目前多种数值方法被成功应用于辐射导热耦合问题的求解, 然而非规则形状介质内辐射导热耦合传热的高精度数值求解仍旧面临较大的挑战[30]. 对于含有弯曲壁面的非规则形状介质而言, 辐射强度在界面附近等局部区域的变化率较大, 而传统连续型数值方法则是基于全局形函数对未知量进行加权近似, 没有考虑局部区域的大梯度辐射强度分布, 计算结果往往存在较大误差. 如文献[31]中即使采用较为细密的网格划分, 采用有限体积法所得辐射热流和温度分布曲线仍旧存在非常明显的振荡, 因此, 扩展高精度数值算法求解辐射导热耦合传热方程, 获得非规则形状介质内辐射导热耦合问题的高精度数值解具有重要意义.

      间断有限元法(discontinuous finite element method, DFEM)[32]是在传统连续型有限元法(finite element method, FEM)基础上发展起来的一种高精度数值方法. 利用传统FEM的空间网格划分, DFEM利用相邻单元边界的数值通量从而释放了FEM强加的内部单元连续性. 由于DFEM形函数构造和控制方程求解均在每个离散单元上进行, 因此该方法兼具较好的几何灵活性和局部守恒性, 并被逐渐应用于辐射传输问题的求解[33-38]. 本文采用DFEM求解辐射传输方程和能量扩散方程, 将其扩展应用于辐射导热耦合传热问题的求解, 验证了数值结果的精确性并给出了典型非规则形状介质内辐射导热耦合传热问题的高精度数值解.

    • 辐射导热耦合传热过程的控制方程包括辐射传输方程和能量扩散方程, 本文的求解思路是先求解辐射传输方程得到辐射强度信息, 然后再把辐射源项代入到能量扩散方程以得到温度场分布.

      首先分析辐射传输方程. 参与性介质内辐射传输方程的离散坐标形式可写为

      ${{{\varOmega }}^m} \cdot \nabla I({{r}},{{{\varOmega }}^m}) + \beta I({{r}},{{{\varOmega }}^m}) = S({{r}},{{{\varOmega }}^m}),$

      式中, I表示辐射强度, r表示坐标, Ω表示离散方向, β为衰减系数, 上标m表示离散方向. 为了书写简洁, $I({{r}}, {{{\varOmega }}^m})$在以下表达中缩写成I m. 考虑介质发射和散射, (1)式中的源项$S({{r}}, {{{\varOmega }}^m})$表示为

      ${S^m} = {\kappa _{\rm{a}}}{I_{\rm{b}}} + \frac{{{\kappa _{\rm{s}}}}}{{4{\text{π}} }}\sum\limits_{m' = 1,m' \ne m}^M {{I^{m'}}{\varPhi ^{m',m}}{\nu ^{m'}}}, $

      式中, Ib表示黑体辐射强度, κa表示发射系数, κs表示散射系数, M = Nθ × Nφ表示角度离散个数, νm表示Ωm方向的权重, Φm′, m表示从方向Ωm到方向Ωm的散射相函数.

      在DFEM应用过程中, 计算区域首先被划分为紧密相邻的离散单元, 如图1(a)所示. 待求解的目标量(如辐射强度)在相邻单元边界上被认为是非连续的, 每个单元上的数值解是相互独立的, 如图1(b)所示. DFEM在每个离散单元上求解控制方程, 相邻单元之间通过数值通量相互连接, 从而保证了计算区域的连续性.

      图  1  (a)空间网格; (b)相邻单元间数值通量示意图

      Figure 1.  (a) Spatial mesh; (b) sketch of numerical flux across the adjacent elements.

      图1(b)中的单元e作为研究对象进行分析, 在该单元上对(1)式使用权函数$w\, ({{r}})$进行加权积分可得

      $\int_{{A_e}} {w\,({{r}}){{{\varOmega }}^m} \cdot \nabla {I^m}} {\rm{d}}A = \int_{{A_e}} {w\,({{r}})( - \beta {I^m} + {S^m})} {\rm{d}}A.$

      考虑辐射强度在相邻边界上的非连续性, 将高斯散度定理应用于(3)式可得

      $\begin{split} &\int_{{A_e}} {w\,({{r}})\,{{{\varOmega }}^m} \cdot \nabla {I^m}} {\rm{d}}A + \int_{{\varGamma _e}} {w({{r}})\,\rlap{–} I_{{\varGamma _e}}^m\,{{{n}}_{{\varGamma _e}}} \cdot {{{\varOmega }}^m}} {\rm{d}}\varGamma \\=\, & \int_{{A_e}} {( - \beta {I^m} + {S^m})w({{r}})} {\rm{d}}A, \\[-15pt]\end{split}$

      式中, ${{{n}}_{{\varGamma _e}}}$表示图1(b)所示的离散单元边界外法向单位向量, Γe表示离散单元边界, $\rlap{–} I_{{\varGamma _e}}^m$表示相邻边界上的数值通量, 定义为

      $\rlap{–} I_{{\varGamma _e}}^m = I_{{\varGamma _e},{\rm{nei}}}^m - I_{{\varGamma _e}}^m,$

      式中$I_{{\varGamma _e}, {\rm{nei}}}^m$表示图1(b)所示的与单元e相邻单元边界上的辐射强度, 定义为

      $I_{{\varGamma _e},{\rm{nei}}}^m = \mathop {\lim }\limits_{\delta \to 0} I_{{\varGamma _e}}^m({{r}} + \delta {{{n}}_{{\varGamma _e}}}).$

      本文数值通量的格式选取为迎风格式

      $\rlap{–} I_{{\varGamma _k}}^m = \left\{\!\!\begin{aligned} & \rlap{-} I_{{\varGamma _k}}^m = I_{{\varGamma _k},{\rm{nei}}}^m - I_{{\varGamma _k}}^m, \quad \,{{{n}}_{{\varGamma _k}}} \cdot {{{\varOmega }}^m} \! < \! 0, \\ & \qquad \qquad \qquad \quad\qquad ~~ k = 1,2,3, \\ & 0,\quad \quad {{{n}}_{{\varGamma _k}}} \cdot {{{\varOmega }}^m} \!> \!0, \;k \!=\! 1,2,3, \end{aligned} \right.$

      因此(4)式中的${\rlap{-} I^m}\, {{{n}}_{{\varGamma _e}}} \cdot {{{\varOmega }}^m}$可以写为

      $\begin{split} & \rlap{–} I_{{\varGamma _k}}^m{{{n}}_{{\varGamma _k}}} \cdot {{{\varOmega }}^m} = \max (0, - {{{n}}_{{\varGamma _k}}} \cdot {{{\varOmega }}^m})\left( {I_{{\varGamma _k},{\rm{nei}}}^m - I_{{\varGamma _k}}^m} \right),\\ & \qquad ~\qquad \qquad k = 1,2,3.\\[-10pt] \end{split}$

      在每一个单元上对未知量采用形函数近似, 单元e上的辐射强度分布可表示为

      ${I^m}({{r}}) \cong \sum\limits_{i = 1}^{{N_{\rm{n}}} = 3} {{\phi _i}({{r}})I_i^m} ,$

      式中, ${\phi _i}({{r}})$表示定义在单元e上的形函数, $I_i^m$表示第i个节点在第m个离散方向上的辐射强度, Nn = 3表示单元e的3个节点. 将(9)式代入(4)式中, 并采用形函数${\phi _i}({{r}})$作为权重函数$w\, ({{r}})$, 可得单元e上辐射传输方程的DFEM离散格式为

      ${{{K}}^m}{{{I}}^m} = {{{H}}^m},$

      式中Im = [I1, I2, I3]T表示离散节点上的辐射强度; 矩阵KmHm的元素表示为

      $\begin{split}{{K}}_{ji}^m =\; & {{{\varOmega }}^m} \cdot \int_{{A_e}} {\nabla {\phi _i}{\phi _j}{\rm{d}}A} + \int_{{A_e}} {\beta {\phi _i}{\phi _j}{\rm{d}}A} \\ & + \sum\limits_{k = 1}^{{N_{\rm{b}}} = 3} {\max (0, - {{{n}}_{{\varGamma _k}}} \cdot {{{\varOmega }}^m})\int_{{\varGamma _k}} {{\phi _i}{\phi _j}{\rm{d}}\varGamma } } , \end{split} \tag{11a}$

      $\begin{split} {{H}}_j^m =\; & {{{\varOmega }}^m} \cdot \int_{{A_e}} {{\phi _i}{\phi _j}S_i^m{\rm{d}}A} \\ &+ \sum\limits_{k = 1}^{{N_{\rm{b}}} = 3} {\max(0, - {{{n}}_{{\varGamma _k}}} \cdot {{{\varOmega }}^m})\int_{{\varGamma _k}} {{\phi _i}{\phi _j}I_{i,{\rm{nei}}}^m{\rm{d}}\varGamma } } , \end{split}\tag{11b}$

      其中Nb = 3表示单元e的3个边界.

      然后分析能量扩散方程. 能量扩散方程表示为

      $ \nabla \cdot [k\nabla T({{r}})] = \nabla \cdot {{q}}({{r}}),$

      式中, k表示热导率; $\nabla \cdot {{q}}({{r}})$为热辐射源项, 表示为

      $\nabla \cdot {{q}}({{r}}) = {\kappa _{\rm{a}}}\left[ {4{\text{π}}{I_{\rm{b}}}({{r}}) - \int_{4{\text{π}}} {I({{r}},{{\varOmega }}){\rm{d}}{ \varOmega}} } \right].$

      与辐射强度类似, 单元e上温度的形函数近似表示为

      $T({{r}}) \cong \sum\limits_{i = 1}^{{N_{\rm{n}}} = 3} {{T_i}{\phi _i}({{r}})} .$

      温度的数值通量表示为

      ${\rlap{-} T_{{\varGamma _k}}} = \left\{ \begin{aligned} & {T_{{\varGamma _k},{\rm{nei}}}} - {T_{{\varGamma _k}}}\quad \,\nabla \cdot {{{n}}_{{\varGamma _k}}} < 0, \; k = 1,2,3, \\ & 0\quad \quad \quad \quad\quad \;\;\;\;\nabla \cdot {{{n}}_{{\varGamma _k}}} > 0, \; k = 1,2,3. \end{aligned} \right.$

      最后可得能量扩散方程的DFEM离散格式为

      ${{MT}} = {{N}},$

      式中, T = [T1, T2, T3]T包含离散节点上的温度; 矩阵MN的元素表示为

      $\begin{split} {{{M}}_{ji}} =\; & k\int_{{A_e}} {\nabla {\phi _i}\nabla {\phi _j}{\rm{d}}A} \\ &+ k\sum\limits_{k = 1}^{{N_{\rm{b}}} = 3} {\max (0, - \nabla \cdot {{{n}}_{{\varGamma _k}}})\int_{{\varGamma _k}} {{\phi _i}\nabla {\phi _j}{\rm{d}}\varGamma } }, \end{split}\tag{17a} $

      $\begin{split} {{{N}}_j} =\, & \int_{{A_e}} {\left\{ {{\kappa _{\rm{a}}}\left[ {4{\text{π}}{I_{\rm{b}}} - \int_{4{\text{π}}} {I({{r}},{{\varOmega }}){\rm{d}}\varOmega } } \right]{\phi _j}} \right\}} {\rm{d}}A \\ &+ \sum\limits_{k = 1}^{{N_{\rm{b}}} = 3} {\max \left( {0, - \nabla \cdot {{{n}}_{{\varGamma _k}}}} \right)\int_{{\varGamma _k}} {{\phi _i}\nabla {\phi _j}{T_{i,{\rm{nei}}}}{\rm{d}}\varGamma } } .\end{split}\tag{17b}$

      至此, 辐射导热耦合传热控制方程的DFEM离散已经完成. 依次求解(10)式所示的辐射传输方程和(16)式所示的能量扩散方程即可得到离散节点的温度值.

    • 基于上述DFEM离散格式, 发展了MATLAB数值计算程序用以求解(11)和(17)式中的各项元素, 进而通过求解方程(10)和(16)得到辐射强度和温度分布. 本节将DFEM应用于求解二维方形介质内的辐射导热耦合传热以验证数值模型的正确性. 方形介质边长为L, 衰减系数为β, 光学厚度为τ = βL = 1.0. 所有壁面均为发射率ε = 1.0的黑壁面, 底边温度为Tb = 1000 K, 其他壁面温度为500 K, 普朗克数为Npl = /(4σTb3), 其中k表示热导率, σ表示Stefan-Boltzmann常数, σ = 5.67 × 10–8.

      首先进行网格无关性检验. 分别将方形区域的边界均匀划分为5, 10和15个单元, 二维计算区域离散为Ne = 62, 226和504个三角形单元. 图2(a)所示为角度划分为Nθ × Nφ = 10 × 20, 不同离散单元条件下介质对称线x/L = 0.5上的无量纲温度T/Tb的分布规律. 图2(b)所示为离散单元数Ne = 226, 不同角度划分Nθ × Nφ = 3 × 6, 5 × 10, 10 × 20和15 × 30条件下对称线上的无量纲温度分布. 由图2所示的结果可知, 离散单元Ne = 226, 离散角度Nθ × Nφ = 10 × 20条件下所得的结果满足网格无关性要求, 在该网格划分条件下采用DFEM求解本算例所需时间为25.38 s.

      图  2  方形介质对称线x/L = 0.5上无量纲温度T/Tb分布 (a)不同空间网格划分; (b)不同角度划分

      Figure 2.  Dimensionless temperature T/Tb along the symmetry line x/L = 0.5 for the cases: (a) Different spatial discretization schemes; (b) different angular discretization schemes.

      针对不同普朗克数Npl, DFEM所得无量纲温度分布如图3(a)所示, 并和文献[39]采用离散传递法(discrete transfer method, DTM)所得结果进行了比较. 针对不同的介质散射反照率ω = κa/β, DFEM和DTM所得无量纲温度的对比结果如图3(b)所示. 由图3(b)结果可知DFEM结果和文献DTM结果符合得很好, 验证了DFEM求解辐射导热耦合问题的正确性.

      图  3  方形介质对称线x/L = 0.5上无量纲温度的DFEM结果和DTM结果对比 (a)不同普朗克数; (b)不同散射反照率

      Figure 3.  Comparison of dimensionless temperature along the square medium symmetry line x/L = 0.5 obtained by DFEM and DTM for the cases: (a) Different Planck numbers; (b) different scattering albedos.

      图4所示为在Npl = 0.01, β = 1.0以及ω = 0.0条件下, 不同数值方法所得对称线x/L = 0.5上无量纲温度分布的对比结果. DTM是一种基于光线跟踪的数值方法[39], 避免了控制方程近似处理引起的误差, 因此采用该方法所得结果可视为该问题的基准解. 以DTM结果作为基准, 采用有限体积法(finite volume method, FVM)[29]所得结果最大误差为2.21%, 而本文DFEM结果的最大误差仅为0.68%, 说明DFEM结果更加精确.

      图  4  不同数值方法得到的方形介质对称线上的无量纲温度对比

      Figure 4.  Comparison of dimensionless temperature along the square medium symmetry line x/L = 0.5 obtained by different numerical methods.

    • 首先考虑如图5(a)所示的内含圆形壁面(半径Rc = 0.2 m)的半圆形(半径Rs = 1.0 m)介质内的辐射导热耦合问题. 介质衰减系数β = 1.0 m–1, 内部圆形壁面温度Tc = 400 K, 外部半圆形壁面和底边温度Ts = 300 K, 壁面均为黑壁面且发射率ε = 1.0. 采用DFEM方法求解沿着半径方向的无量纲温度分布. 方向离散为Nθ × Nφ = 20 × 40, 空间离散为如图5(b)所示的1768个三角形单元, 所得结果满足网格无关性要求, 采用DFEM求解本算例所需时间为182.69 s.

      图  5  内含圆形热壁面的半圆形介质 (a)结构示意图; (b)网格划分

      Figure 5.  Semicircle medium with an inner circle hot boundary: (a) Geometry sketch; (b) spatial discretization.

      在普朗克数Npl = /(4σTs3) = 0.1条件下, 图6所示为采用不同方法(Fluent[31], 不同网格处理方式的FVM[31], 混合格子Boltzmann-有限体积法(LBM-FVM)[40]以及本文DFEM)所得的半圆形介质对称线x = 0上的温度分布. 由图6可知采用Fluent和LBM-FVM以及本文DFEM所得结果符合得较好; 而文献[31]采用FVM所得结果则具有较大的偏差, 这是文献[21]采用四边形网格处理弯曲壁面引起的计算误差所致. 图6所示结果表明, 和FVM相比, 本文发展的DFEM求解非规则形状介质内的辐射导热耦合问题更加精确.

      图  6  不同数值方法得到的半圆介质对称线上温度分布 (a) y = [0.0, 0.2]; (b) y = [0.6, 1.0]

      Figure 6.  Temperature distributions along the symmetric line of the semicircle medium obtained by different numerical algorithms: (a) y = [0.0, 0.2]; (b) y = [0.6, 1.0].

      图7所示为普朗克数Npl = 0.1和1.0条件下半圆形介质底边y = 0上的总热流密度${q_{{\rm{total}}}}({{{r}}_{\rm{w}}}) = $$ - \left( {k\dfrac{{\partial T({{{r}}_{\rm{w}}})}}{{\partial n}} + \displaystyle\int_{{{{n}}_{\rm{w}}} \cdot {{\varOmega }} > 0} {I({{{r}}_{\rm{w}}}, {{\varOmega }})({{{n}}_{\rm{w}}} \cdot {{\varOmega }}){\rm{d}}{{\varOmega }}} } \right)$(负号表示方向)的分布规律. 由图7可知底边中点位置处的热流密度最强, 这是由于内部圆环温度较高, 底边中点距离高温圆环最近, 在辐射和导热共同作用条件下, 该位置处的换热强度达到最大. 由图7(a)结果可知, 当Npl = 0.1时, DFEM结果和Fluent结果吻合良好, 而LBM-FVM结果存在明显误差. 当普朗克数增大到Npl = 1.0时, 采用不同方法所得结果符合得比较好, 如图7(b)所示. 图7所示结果表明, 对于普朗克数较小即辐射占优的辐射导热耦合传热问题, DFEM所得结果更加精确.

      图  7  不同普朗克数条件下半圆介质底边上总热流密度分布 (a) Npl = 0.1, (b) Npl = 1.0

      Figure 7.  Total flux distributions along the bottom boundary of the semicircle medium under the situations with different Plank numbers: (a) Npl = 0.1; (b) Npl = 1.0.

      接下来考虑如图8(a)所示的内含圆形热边界的非规则形状介质, 该模型可用来研究圆形热管和环境之间的辐射导热耦合换热. 介质尺寸(单位为m)在图8(a)中给出, 介质衰减系数为β = 1.0 m–1, 内部圆形边界温度为400 K, 其他边界温度为自然环境温度300 K, 所有壁面均为黑壁面且发射率为ε = 1.0. 计算区域离散为图8(b)所示的1142个三角形单元、方向离散为Nθ × Nφ = 20 × 40条件下所得结果达到网格无关性要求, 采用DFEM求解本算例所需时间为118.38 s.

      图  8  内含圆形热边界的非规则形状介质 (a)结构示意图; (b)网格划分

      Figure 8.  Irregular medium with an inner hot boundary: (a) Geometry sketch; (b) spatial discretization.

      图9(a)所示为不同普朗克数Npl = 0.1和 1.0条件下非规则介质中线x = 0.5上的温度分布. 由图9(a)可知, 对于Npl = 0.1的辐射占优问题, 从壁面到圆环之间的温度梯度先减小后增加, 说明靠近圆环和弯曲壁面处的介质温度变化较为剧烈, 而中间区域介质的温度变化较为平缓; 对于Npl = 1.0的导热占优问题, 该区域内的温度梯度逐渐增加, 说明弯曲壁面附近处的介质温度变化较为平缓, 高温圆环附近处的介质温度变化剧烈; 圆环上方的介质温度变化剧烈程度也有所区别, 但由于该区间垂直高度较小, 因此该区间范围的介质温度沿y方向接近于线性减小.

      图  9  (a)普朗克数Npl = 0.1和1.0时内含圆形热边界的非规则形状介质中线上温度分布; (b) Npl = 0.1时介质温度分布; (c) Npl = 1.0时介质温度分布

      Figure 9.  (a) Temperature distributions along the centerline of the irregular medium with an inner hot boundary; (b) temperature distribution within the computation domain for the case of Npl = 0.1; (c) temperature distribution within the computation domain for the case of Npl = 1.0.

      在中心线上的同一位置处, 普朗克数Npl = 0.1对应的介质温度总是高于Npl = 1.0对应的介质温度. 这是由于当Npl = 0.1时, 辐射传热作用显著, 相同位置的辐射源项较强, 因此介质温度更高. 图9(b)图9(c)所示分别为同普朗克数Npl = 0.1和 1.0条件下的介质温度分布的等值线图, 图9中相邻两条等值线之间的温度差为10 K. 由图9(b)图9(c)的对比结果可知, 当普朗克数Npl = 0.1时, 强烈的辐射效应能够使高温圆环能量传播得更远, 因此介质内的能量分布更加均匀; 当普朗克数Npl = 1.0时, 显著的导热效应导致温度从高温圆环向外迅速降低, 因此外部边界附近较大范围内的介质温度处于较低水平.

      进一步考虑如图10(a)所示的内含两个高温圆环的矩形介质, 该模型可用来研究两个高温管道及其所处环境之间的辐射导热耦合换热. 矩形介质边长为Lx × Ly = 1.0 m × 1.0 m, 内部圆环半径为R = 0.1 m, 圆心位置分别为(x, y) = (0.3 m, 0.3 m)和(x, y) = (0.7 m, 0.7 m). 矩形介质边界和圆环边界之间充满衰减系数为β = 1.0 m–1的参与性介质, 内部两个圆形边界温度均为400 K, 其他边界温度为300 K, 所有壁面均为黑壁面且发射率为ε = 1.0. 计算区域离散为图10(b)所示的1262个三角形单元、方向离散为Nθ × Nφ = 20 × 40条件下所得到的结果达到了网格无关性的要求, 采用DFEM求解本算例所需时间为131.75 s.

      图  10  内含两个圆形热边界的矩形介质 (a)结构示意图; (b)网格划分

      Figure 10.  The medium the square medium with two circular hot boundaries: (a) Geometry sketch; (b) spatial discr etization.

      在不同普朗克数Npl = 0.1和 1.0条件下, 中线x = 0.5上的温度分布如图11(a)所示, 可知介质中线温度在y = [0, 0.3]逐渐升高且温度变化剧烈程度逐渐降低, 再y = [0.3, 0.7]的介质保持较高的温度, 在y = [0.7, 1.0]逐渐降低. 在同一位置处, Npl = 0.1对应的介质温度高于Npl = 1.0对应的介质温度, 说明辐射效应作用对该模型中线上的介质温度具有显著的提升作用. 图11(b)图11(c)所示分别为同普朗克数Npl = 0.1和 1.0条件下计算区域内温度分布的等值线图, 相邻两条等值线之间的温度差为5 K. 由图11(b)图11(c)所示结果可知, 介质中心点位置(x, y) = (0.5 m, 0.5 m)周围存在一个温度梯度较小的区域, 且普朗克数越小, 该区域的面积越大. 以上结果表明辐射效应会弱化中心点位置处的温度梯度, 而导热效应则会强化该位置处的温度梯度.

      图  11  (a)普朗克数Npl = 0.1和1.0时内含两个圆形热边界的矩形介质中线上温度分布; (b) Npl = 0.1时介质温度分布; (c) Npl = 1.0时介质温度分布

      Figure 11.  (a) Temperature distributions along the centerline of the square medium with two circular hot boundaries; (b) temperature distribution within the computation domain for the case of Npl = 0.1; (c) temperature distribution within the computation domain for the case of Npl = 1.0

    • 本文将DFEM扩展应用于求解辐射导热耦合传热问题, 采用非结构化三角形单元划分计算区域, 给出了辐射传输方程和能量扩散方程的DFEM离散格式, 验证了数值算法的正确性并研究了非规则形状介质内的辐射导热耦合传热问题. 研究结果表明: 本文发展的DFEM能够得到非规则形状介质内辐射导热耦合传热问题的高精度数值结果, 尤其是对于普朗克数较小的辐射占优问题; 非规则形状介质内的温度分布规律表明辐射效应会弱化热边界附近的温度梯度, 而导热效应则会强化热边界附近的温度梯度.

参考文献 (40)

目录

    /

    返回文章
    返回