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

王存海 郑树 张欣欣

1) (北京科技大学能源与环境工程学院, 北京 100083)

2) (冶金工业节能减排北京市重点实验室, 北京 100083)

3) (华北电力大学能源动力与机械工程学院, 北京 102206)

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

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,Ωm) 在以下表达中缩写成Im.考虑介质发射和散射, (1)式中的源项 S (r,Ωm) 表示为

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

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

图1 (a)空间网格; (b)相邻单元间数值通量示意图Fig.1.(a) Spatial mesh; (b) sketch of numerical flux across the adjacent elements.

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

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

式中, nΓe表示图1(b)所示的离散单元边界外法向单位向量, Γe表示离散单元边界,表示相邻边界上的数值通量, 定义为

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

因此(4)式中的 -ImnΓe·Ωm可以写为

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

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

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

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

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

式中, k表示热导率; ∇ ·q(r) 为热辐射源项, 表示为

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

温度的数值通量表示为

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

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

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

3 模型验证

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

图2 方形介质对称线x/L = 0.5上无量纲温度T/Tb分布 (a)不同空间网格划分; (b)不同角度划分Fig.2.Dimensionless temperature T/Tb along the sym−metry line x/L = 0.5 for the cases：(a) Different spatial dis−cretization schemes; (b) different angular discretization schemes.

首先进行网格无关性检验.分别将方形区域的边界均匀划分为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.

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

图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结果更加精确.

图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 num−bers; (b) different scattering albedos.

图4 不同数值方法得到的方形介质对称线上的无量纲温度对比Fig.4.Comparison of dimensionless temperature along the square medium symmetry line x/L = 0.5 obtained by differ−ent 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.

在普朗克数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求解非规则形状介质内的辐射导热耦合问题更加精确.

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

图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.0Fig.7.Total flux distributions along the bottom boundary of the semicircle medium under the situations with different Plank num−bers：(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.

接下来考虑如图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.

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

在中心线上的同一位置处, 普朗克数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时, 显著的导热效应导致温度从高温圆环向外迅速降低, 因此外部边界附近较大范围内的介质温度处于较低水平.

图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 dis−tribution 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.

进一步考虑如图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.

在不同普朗克数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时介质温度分布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.

5 结 论

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