摘要: 采用间断有限元法(discontinuous finite element method, DFEM)求解非规则形状介质内的辐射导热耦合传热问题, 得到了典型非规则形状介质内辐射导热耦合传热问题的高精度数值结果. 和传统连续型有限元方法不同, DFEM将计算区域划分成相互独立的离散单元, 形函数的构造、未知量的加权近似以及控制方程的求解均在每一个离散单元上进行. 通过在单元之间施加迎风格式的数值通量, DFEM保证了整个计算区域的连续性, 因此这种方法兼具良好的几何灵活性和局部守恒性. 推导了辐射传输方程和能量扩散方程的DFEM离散格式, 验证了DFEM求解辐射导热耦合传热问题的正确性; 同时研究了不同几何形状介质内辐射导热耦合传热问题, 得到了典型非规则形状介质内辐射导热耦合传热的高精度数值结果.
关键词: 辐射导热 /
耦合传热 /
数值模拟 /
间断有限元 English Abstract Discontinuous finite element solutions for coupled radiation-conduction heat transfer in irregular media Wang Cun-Hai 1,2 ,Zheng Shu 3 ,Zhang Xin-Xin 1,2 1.School of Energy and Environmental Engineering, University of Science and Technology Beijing, Beijing 100083, China 2.Beijing Key Laboratory of Energy Conservation and Emission Reduction for Metallurgical Industry, Beijing 100083, China 3.School of Energy, Power and Mechanical Engineering, North China Electric Power University, Beijing 102206, China Fund Project: the National Natural Science Foundation of China (Grant Nos. 51906014, 51890891), the China Postdoctoral Science Foundation (Grant No. 2018M641196), and the Fundamental Research Funds for the Central Universities (Grant No. FRF-TP-18-072A1) Received Date: 03 August 2019Accepted Date: 25 October 2019Published Online: 05 February 2020 Abstract: The discontinuous finite element method (DFEM) is used to investigate the coupled radiation-conduction heat transfer in an irregular medium, and the highly accurate solutions for several typical media are numerically obtained. Comparing with the traditional continuous finite element method, the computational domain in the DFEM application is discretized into unstructured meshes that are assumed to be separated from each other. The shape function construction, field variable approximation, and numerical solutions are obtained for every single element. The continuity of the computational domain is maintained by modeling a numerical flux with the up-winding scheme. Thus the DFEM has the salient feature of geometry flexibility and simultaneously supports locally conservative solutions. The DFEM discretization for the radiative transfer equation and the energy diffusion equation are first presented, and the accuracies of the DFEM for coupled radiation-conduction heat transfer problems are verified. Combined radiation-conduction heat transfer problems in several irregular media are afterward solved, and the highly accurate DFEM solutions are presented. Keywords: radiation-conduction /coupled heat transfer /numerical simulation /discontinuous finite element method 全文HTML --> --> --> 1.引 言 参与性介质内辐射导热耦合传热过程广泛存在于众多工程实际应用中[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求解辐射传输方程和能量扩散方程, 将其扩展应用于辐射导热耦合传热问题的求解, 验证了数值结果的精确性并给出了典型非规则形状介质内辐射导热耦合传热问题的高精度数值解.2.控制方程及DFEM离散 辐射导热耦合传热过程的控制方程包括辐射传输方程和能量扩散方程, 本文的求解思路是先求解辐射传输方程得到辐射强度信息, 然后再把辐射源项代入到能量扩散方程以得到温度场分布. 首先分析辐射传输方程. 参与性介质内辐射传输方程的离散坐标形式可写为 式中, I 表示辐射强度, r 表示坐标, Ω 表示离散方向, β 为衰减系数, 上标m 表示离散方向. 为了书写简洁, $I({{r}}, {{{\varOmega }}^m})$ 在以下表达中缩写成I m . 考虑介质发射和散射, (1 )式中的源项$S({{r}}, {{{\varOmega }}^m})$ 表示为 式中, I b 表示黑体辐射强度, κ a 表示发射系数, κ s 表示散射系数, M = Nθ × Nφ 表示角度离散个数, ν m ′ 表示Ω m ′ 方向的权重, Φ m ′, m 表示从方向Ω m ′ 到方向Ω m 的散射相函数. 在DFEM应用过程中, 计算区域首先被划分为紧密相邻的离散单元, 如图1(a) 所示. 待求解的目标量(如辐射强度)在相邻单元边界上被认为是非连续的, 每个单元上的数值解是相互独立的, 如图1(b) 所示. DFEM在每个离散单元上求解控制方程, 相邻单元之间通过数值通量相互连接, 从而保证了计算区域的连续性. 图 1 (a)空间网格; (b)相邻单元间数值通量示意图 Figure1. (a) Spatial mesh; (b) sketch of numerical flux across the adjacent elements. 以图1(b) 中的单元e 作为研究对象进行分析, 在该单元上对(1 )式使用权函数$w\, ({{r}})$ 进行加权积分可得 考虑辐射强度在相邻边界上的非连续性, 将高斯散度定理应用于(3 )式可得 式中, ${{{n}}_{{\varGamma _e}}}$ 表示图1(b) 所示的离散单元边界外法向单位向量, Γ e 表示离散单元边界, $\rlap{–} I_{{\varGamma _e}}^m$ 表示相邻边界上的数值通量, 定义为 式中$I_{{\varGamma _e}, {\rm{nei}}}^m$ 表示图1(b) 所示的与单元e 相邻单元边界上的辐射强度, 定义为 本文数值通量的格式选取为迎风格式 因此(4 )式中的${\rlap{-} I^m}\, {{{n}}_{{\varGamma _e}}} \cdot {{{\varOmega }}^m}$ 可以写为 在每一个单元上对未知量采用形函数近似, 单元e 上的辐射强度分布可表示为 式中, ${\phi _i}({{r}})$ 表示定义在单元e 上的形函数, $I_i^m$ 表示第i 个节点在第m 个离散方向上的辐射强度, N n = 3表示单元e 的3个节点. 将(9 )式代入(4 )式中, 并采用形函数${\phi _i}({{r}})$ 作为权重函数$w\, ({{r}})$ , 可得单元e 上辐射传输方程的DFEM离散格式为 式中I m = [I 1 , I 2 , I 3 ]T 表示离散节点上的辐射强度; 矩阵K m 和H m 的元素表示为 其中N b = 3表示单元e 的3个边界. 然后分析能量扩散方程. 能量扩散方程表示为 式中, k 表示热导率; $\nabla \cdot {{q}}({{r}})$ 为热辐射源项, 表示为 与辐射强度类似, 单元e 上温度的形函数近似表示为 温度的数值通量表示为 最后可得能量扩散方程的DFEM离散格式为 式中, T = [T 1 , T 2 , T 3 ]T 包含离散节点上的温度; 矩阵M 和N 的元素表示为 至此, 辐射导热耦合传热控制方程的DFEM离散已经完成. 依次求解(10 )式所示的辐射传输方程和(16 )式所示的能量扩散方程即可得到离散节点的温度值.3.模型验证 基于上述DFEM离散格式, 发展了MATLAB数值计算程序用以求解(11 )和(17 )式中的各项元素, 进而通过求解方程(10 )和(16 )得到辐射强度和温度分布. 本节将DFEM应用于求解二维方形介质内的辐射导热耦合传热以验证数值模型的正确性. 方形介质边长为L , 衰减系数为β , 光学厚度为τ = βL = 1.0. 所有壁面均为发射率ε = 1.0的黑壁面, 底边温度为T b = 1000 K, 其他壁面温度为500 K, 普朗克数为N pl = kβ /(4σ T b 3 ), 其中k 表示热导率, σ 表示Stefan-Boltzmann常数, σ = 5.67 × 10–8 . 首先进行网格无关性检验. 分别将方形区域的边界均匀划分为5, 10和15个单元, 二维计算区域离散为N e = 62, 226和504个三角形单元. 图2(a) 所示为角度划分为Nθ × Nφ = 10 × 20, 不同离散单元条件下介质对称线x /L = 0.5上的无量纲温度T /T b 的分布规律. 图2(b) 所示为离散单元数N e = 226, 不同角度划分Nθ × Nφ = 3 × 6, 5 × 10, 10 × 20和15 × 30条件下对称线上的无量纲温度分布. 由图2 所示的结果可知, 离散单元N e = 226, 离散角度Nθ × Nφ = 10 × 20条件下所得的结果满足网格无关性要求, 在该网格划分条件下采用DFEM求解本算例所需时间为25.38 s. 图 2 方形介质对称线x /L = 0.5上无量纲温度T /T b 分布 (a)不同空间网格划分; (b)不同角度划分 Figure2. Dimensionless temperature T /T b along the symmetry line x /L = 0.5 for the cases: (a) Different spatial discretization schemes; (b) different angular discretization schemes. 针对不同普朗克数N pl , 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)不同散射反照率 Figure3. 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 所示为在N pl = 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 不同数值方法得到的方形介质对称线上的无量纲温度对比 Figure4. Comparison of dimensionless temperature along the square medium symmetry line x /L = 0.5 obtained by different numerical methods. 4.结果和讨论 首先考虑如图5(a) 所示的内含圆形壁面(半径R c = 0.2 m)的半圆形(半径R s = 1.0 m)介质内的辐射导热耦合问题. 介质衰减系数β = 1.0 m–1 , 内部圆形壁面温度T c = 400 K, 外部半圆形壁面和底边温度T s = 300 K, 壁面均为黑壁面且发射率ε = 1.0. 采用DFEM方法求解沿着半径方向的无量纲温度分布. 方向离散为Nθ × Nφ = 20 × 40, 空间离散为如图5(b) 所示的1768个三角形单元, 所得结果满足网格无关性要求, 采用DFEM求解本算例所需时间为182.69 s. 图 5 内含圆形热壁面的半圆形介质 (a)结构示意图; (b)网格划分 Figure5. Semicircle medium with an inner circle hot boundary: (a) Geometry sketch; (b) spatial discretization. 在普朗克数N pl = kβ /(4σ T s 3 ) = 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] Figure6. 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 所示为普朗克数N pl = 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) 结果可知, 当N pl = 0.1时, DFEM结果和Fluent结果吻合良好, 而LBM-FVM结果存在明显误差. 当普朗克数增大到N pl = 1.0时, 采用不同方法所得结果符合得比较好, 如图7(b) 所示. 图7 所示结果表明, 对于普朗克数较小即辐射占优的辐射导热耦合传热问题, DFEM所得结果更加精确. 图 7 不同普朗克数条件下半圆介质底边上总热流密度分布 (a) N pl = 0.1, (b) N pl = 1.0 Figure7. Total flux distributions along the bottom boundary of the semicircle medium under the situations with different Plank numbers: (a) N pl = 0.1; (b) N pl = 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)网格划分 Figure8. Irregular medium with an inner hot boundary: (a) Geometry sketch; (b) spatial discretization. 图9(a) 所示为不同普朗克数N pl = 0.1和 1.0条件下非规则介质中线x = 0.5上的温度分布. 由图9(a) 可知, 对于N pl = 0.1的辐射占优问题, 从壁面到圆环之间的温度梯度先减小后增加, 说明靠近圆环和弯曲壁面处的介质温度变化较为剧烈, 而中间区域介质的温度变化较为平缓; 对于N pl = 1.0的导热占优问题, 该区域内的温度梯度逐渐增加, 说明弯曲壁面附近处的介质温度变化较为平缓, 高温圆环附近处的介质温度变化剧烈; 圆环上方的介质温度变化剧烈程度也有所区别, 但由于该区间垂直高度较小, 因此该区间范围的介质温度沿y 方向接近于线性减小. 图 9 (a)普朗克数N pl = 0.1和1.0时内含圆形热边界的非规则形状介质中线上温度分布; (b) N pl = 0.1时介质温度分布; (c) N pl = 1.0时介质温度分布 Figure9. (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 N pl = 0.1; (c) temperature distribution within the computation domain for the case of N pl = 1.0. 在中心线上的同一位置处, 普朗克数N pl = 0.1对应的介质温度总是高于N pl = 1.0对应的介质温度. 这是由于当N pl = 0.1时, 辐射传热作用显著, 相同位置的辐射源项较强, 因此介质温度更高. 图9(b) 和图9(c) 所示分别为同普朗克数N pl = 0.1和 1.0条件下的介质温度分布的等值线图, 图9 中相邻两条等值线之间的温度差为10 K. 由图9(b) 和图9(c) 的对比结果可知, 当普朗克数N pl = 0.1时, 强烈的辐射效应能够使高温圆环能量传播得更远, 因此介质内的能量分布更加均匀; 当普朗克数N pl = 1.0时, 显著的导热效应导致温度从高温圆环向外迅速降低, 因此外部边界附近较大范围内的介质温度处于较低水平. 进一步考虑如图10(a) 所示的内含两个高温圆环的矩形介质, 该模型可用来研究两个高温管道及其所处环境之间的辐射导热耦合换热. 矩形介质边长为L x × L y = 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)网格划分 Figure10. The medium the square medium with two circular hot boundaries: (a) Geometry sketch; (b) spatial discr etization. 在不同普朗克数N pl = 0.1和 1.0条件下, 中线x = 0.5上的温度分布如图11(a) 所示, 可知介质中线温度在y = [0, 0.3]逐渐升高且温度变化剧烈程度逐渐降低, 再y = [0.3, 0.7]的介质保持较高的温度, 在y = [0.7, 1.0]逐渐降低. 在同一位置处, N pl = 0.1对应的介质温度高于N pl = 1.0对应的介质温度, 说明辐射效应作用对该模型中线上的介质温度具有显著的提升作用. 图11(b) 和图11(c) 所示分别为同普朗克数N pl = 0.1和 1.0条件下计算区域内温度分布的等值线图, 相邻两条等值线之间的温度差为5 K. 由图11(b) 和图11(c) 所示结果可知, 介质中心点位置(x , y ) = (0.5 m, 0.5 m)周围存在一个温度梯度较小的区域, 且普朗克数越小, 该区域的面积越大. 以上结果表明辐射效应会弱化中心点位置处的温度梯度, 而导热效应则会强化该位置处的温度梯度. 图 11 (a)普朗克数N pl = 0.1和1.0时内含两个圆形热边界的矩形介质中线上温度分布; (b) N pl = 0.1时介质温度分布; (c) N pl = 1.0时介质温度分布 Figure11. (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 N pl = 0.1; (c) temperature distribution within the computation domain for the case of N pl = 1.0 5.结 论 本文将DFEM扩展应用于求解辐射导热耦合传热问题, 采用非结构化三角形单元划分计算区域, 给出了辐射传输方程和能量扩散方程的DFEM离散格式, 验证了数值算法的正确性并研究了非规则形状介质内的辐射导热耦合传热问题. 研究结果表明: 本文发展的DFEM能够得到非规则形状介质内辐射导热耦合传热问题的高精度数值结果, 尤其是对于普朗克数较小的辐射占优问题; 非规则形状介质内的温度分布规律表明辐射效应会弱化热边界附近的温度梯度, 而导热效应则会强化热边界附近的温度梯度.