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 2019
Accepted Date:25 October 2019
Published 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
针对不同普朗克数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)不同散射反照率 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所示为在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 不同数值方法得到的方形介质对称线上的无量纲温度对比 Figure4. Comparison of dimensionless temperature along the square medium symmetry line x/L = 0.5 obtained by different numerical methods.
4.结果和讨论首先考虑如图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)网格划分 Figure5. Semicircle medium with an inner circle hot boundary: (a) Geometry sketch; (b) spatial discretization.
在普朗克数Npl = kβ/(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] 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所示为普朗克数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 Figure7. 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)网格划分 Figure8. 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时介质温度分布 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 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)网格划分 Figure10. 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时介质温度分布 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 Npl = 0.1; (c) temperature distribution within the computation domain for the case of Npl = 1.0