付丽荣，张文平，明平剑，罗跃生

(1.哈尔滨工程大学动力与能源工程学院，黑龙江哈尔滨150001;2.哈尔滨工程大学理学院，黑龙江哈尔滨150001)

高负荷柴油机的辐射传热占总传热量的10%～45%［1］，直接关系到发动机热效率以及因传热引起的各种热负荷、热强度问题;同时，辐射热流量对燃烧性能的影响很大，因此需对缸内辐射传热进行深入研究［2-3］。柴油机缸内辐射传热的研究主要基于实验［4-6］和数值模拟［7-9］;实验研究法，即用热辐射传感器直接测量气缸内辐射传热量或者用二色法测量气缸内火焰温度，再求得辐射传热量。数值模拟主要利用商用软件，如 Fluent［8］、AVL-Fire［9］。目前主要的辐射传热模拟方法有:区域法、蒙特卡洛法、离散传递法、有限体积法、无网格法［10］等。而有限体积法，可确保辐射能量整体守恒，对不规则边界适应性强，易处理各向异性散射，可以采用与流场计算相同的网格，在迭代计算中可以不对温度和热流进行插值计算。因此本文采用有限体积法进行编程，从而可以与通用流体力学数值模拟软件(general transportation equation analysis，GTEA)［11］的流场计算采用同一套网格，简化了计算程序，提高了计算的通用性。

本文在前期工作基础上，针对三维辐射传热进行数值模拟。本文借助有限体积法计算三维封闭空腔内的辐射传热，并用3个算例——六面体、圆柱体和TBD620型柴油机的ω燃烧室来验证此方法和程序的可靠性与准确性，同时考虑了网格离散个数和立体角离散个数对计算结果的影响。

1 控制方程

1.1 辐射传热方程

对于吸收、发射、散射性灰介质的辐射传递方程表达式如下［12］

式中:κ、σs分别为吸收系数和散射系数，β0=κ+σs为衰减系数;I为空间位置r、方向s处的辐射强度;Ib为黑体辐射强度，Ib=6T4/π;σ为斯蒂芬-玻尔兹曼常数，为 5.67×10-8W/(M2·K4)Φ(s，si)为辐射从入射方向s到散射方向si的散射相函数;Ωi是中心方向为si的立体角。

本文所研究的介质是非散射的，所以方程(1)可以简化为

1.2 辐射传热方程的边界条件

RTE的边界条件定义如下

式中:εw为壁面发射率，下标w表示壁面位置，nw为单位法向量。式(3)表示当辐射能离开一个光学反射面时有2部分的贡献，一是高温表面放出的辐射量，二是投射到此表面的反射量。

2 数值离散

2.1 有限体积法

为了得到离散化方程，将方程(2)对控制体ΔV和立体角ΔΩmn积分，可得ΔΩmn内辐射能量守恒方程的有限体积表达式:

应用Gauss定理把体积分转换成面积分，假设辐射强度I和Ib在控制体和立体角内为常数，同时，源项在之内为均匀的，则方程(4)的左端和右端分别可以表示为如下形式:

将式(5)、(6)代入到方程(4)，可得

式中:

将空间区域离散为互不重叠的四面体或三棱柱单元，所有变量信息存储在单元中心，如图1(a)的点1、2、3、4所示。将4π角域离散为互不重叠的控制立体角 ( Nθ×Nφ)，把方向矢量smn定义为立体角的中心方向，θ是0～π变化的极角，φ是0～2π变化的方位角，如图1(b)所示。控制角被平分成Δθm=θm+-θm-= π/Nθ，Δφn= φn+-φn-=2π/Nφ。当立体角是确定时，方向权的符号决定了控制体表面的辐射能是流进还是流出，它由下式表示:

本文采用一阶迎风格式将单元界面上的辐射强度与单元中心的辐射强度关联起来，则它们的关系表示如下

式中:

将式(13)代入式(7)，则式(7)可写成如下线性方程组:

式中:

图1 混合网格和立体角示意图Fig.1 Diagram of hybrid grids and the solid angle

2.2 边界条件离散

边界条件可以离散成如下形式

式中:

2.3 壁面辐射热流密度

壁面辐射热流密度表示为

3 数值计算结果及分析

本文的辐射换热数值模拟是应用GTEA软件，在文献［13］基础上，为了验证本文的UFVM方法，并且在三维辐射问题上检验所开发的计算程序的可靠性，本文选择2个基本的几何体，即六面体、圆柱体。

3.1 不规则六面体辐射模拟

验证算例1是对不规则六面体区域内的辐射进行仿真并与文献中结果进行对比。几何体结构的建立和参数的设定与文献［12］中的数据一致。如图2(a)所示，其中z0=1 m，所有边界为黑体面，温度为0，壁面发射率为εw=1，内部介质温度为100 K。空间离散为六面体单元，见图 2(b)，(Nx×Ny×Nz)= (27×27×27)，立体角划分为( Nθ×Nφ)= ( 8×16)，吸收系数分别采用10、1和0.1 m-1进行计算，不考虑散射，计算底面的无量纲热流量分布。本文将计算结果(如图3所示)与精确解和数值解［12］进行对比，可以看到本文的计算结果与精确解吻合很好。

图2 六面体模型与网格Fig.2 A hexahedral enclosure and mesh

由图3可见，当吸收系数是10 m-1时，无量纲热流在底面中心处的值接近于1.0，主要是由于底面此时只受到周围介质的影响;在两边处热流迅速减小，因为此时两边受到边界冷壁面的影响。当吸收系数为0.1 m-1时，热流急剧下降，这是因为介质被冷边界冷却的结果。

图3 底面上AB线段无量纲热流比较Fig.3 Comparison of radiative wall heat flux along line AB on the bottom wall of a quadrilateral enclosure

图4 网格的影响Fig.4 The effect of different discrete grids

图5 网格和立体角划分均不同时的误差Fig.5 The error of different both grids and solid angles

图4给出了吸收系数为1 m-1时，立体角离散个数均为 ( Nθ×Nφ)= ( 4×8)，空间区域分别离散为19 683个单元和12 167个单元时计算结果与精确解的相对误差分布。由图可知网格加密后的相对误差不但没有明显的变小，反而在某些位置还有所变大。当空间离散为12 167个单元，立体角离散个数为 ( Nθ×Nφ)= ( 4×8)时，与单元个数为 19 683，立体角离散个数为 ( Nθ×Nφ)= ( 8×16)时，计算结果与精确解的相对误差分析见图5。由图可看出当网格加密且立体角离散个数增多时，计算误差明显减小，计算结果更接近于精确解。因此，在数值计算中，空间离散的网格并不是越密越好，而是要和立体角离散的个数相对应，即如果空间网格离散的个数增加，那么相应的立体角离散的个数也要增加，才可以得到更高精度的数值解。

3.2 圆柱体辐射模拟

验证算例2是半径为1m，z轴方向长度为2m的圆柱体。边界为黑体，温度为0，壁面发射率为εw=1，内部气体温度为100 K，圆柱体和空间网格见图6。在吸收系数分别为0.1，1.0和5.0m-1情况下，计算圆柱体右侧边界的无量纲辐射热流量，空间域离散为74 203个四面体单元，立体角划分为 ( Nθ×Nφ)= ( 8×16)。本文的计算结果与文献［14］提供的精确解和文献［15-16］的数值解进行对比，结果如图7所示，可以看出它与精确解和数值解吻合很好。

图6 圆柱体和网格Fig.6 The cylinder enclosure and mesh

图7 无量纲辐射热流量的对比Fig.7 Comparison of dimensionless radiative heat flux along the lateral side of the cylinder

图7显示了在3种吸收系数0.1、1.0和5.0 m-1下右侧边界的无量纲辐射热流量，当吸收系数是5.0 m-1时，到达右侧壁面中间处的辐射热流密度接近于介质的黑体辐射强度，这是由于壁面的辐射强度只受封闭壁面周围的热气体的发射本领影响。然而，在两边角处，辐射热流迅速减小，因为受周围顶部和底部冷壁面影响。但当吸收系数是0.1 m-1时，介质的辐射本领是弱的，并且无量纲辐射热流量大幅度下降，这是因为受其他冷壁面和光学薄气体可忽略的自身衰减的深远影响。

4 ω型燃烧室辐射模拟

验证算例3是TBD620型柴油机ω型燃烧室。发动机的基本参数:气缸数为16，冲程为195 mm，缸径为170 mm，连杆长度为350 mm，压缩比为13.5。在燃烧结束至排气开始前这一阶段，此时活塞在下止点位置，缸盖、缸壁和活塞温度均为350 K，假设气体平均温度为800 K;利用GAMBIT软件划分网格，通过软件将网格文件转换为CGNS格式，并读入本文开发的求解器;空间域内采用混合网格，气缸内为三棱柱非结构化网格，ω型燃烧室内采用四面体网格，网格总数为63 227，图8(a)、(b)分别显示了计算网格的主视和俯视图;立体角划分为( Nθ×Nφ)= ( 8×16)。计算燃烧室纵向中心截面边界线左侧、下侧和右侧的辐射热流密度，图9中s为0～200、200～432.7、432.7～632.7mm 区间分别代表边界线左侧、下侧和右侧。

图8 计算网格Fig.8 Computational grid

图9 s的定义Fig.9 Definition of s

图10显示的是工质吸收系数对缸内纵向中心截面的左侧、下侧和右侧边界无量纲辐射热流密度的影响。从432.7～632.7 mm区间可以看出3种吸收系数的热流量图都是对称的，与算例2的圆柱体模型结果图7对比可知，算例3的底面的ω型燃烧室与算例2底面是圆形的效果一样，所得的侧面热流量图都是对称的，即底面的形状对侧边的辐射热流量没有影响。在0和200 mm周围，由于两边角处受顶部和底部冷壁面影响，所以热流量值变小。当吸收系数为0.1 m-1时，工质相当于是透明的，所以各位置的热流量几乎一样大。

图10 不同吸收系数下的无量纲辐射热流量的对比Fig.10 Comparison of dimension less radiative heat flux of different absorption coefficients

图11显示的是吸收系数为10.0 m-1时，工质温度分别为600、800、1 000 K时缸壁无量纲辐射热流量。介质温度对中心部位的热流量影响较大，对边角处的影响较小;边角处由于受冷壁面的影响，热流量值变小。

图11 介质不同温度下的无量纲辐射热流量的对比Fig.11 Comparison of dimensionless radiative heat flux of medium of different temperatures

5 结论

本文通过六面体、圆柱体及ω型燃烧室的数值模拟，得出如下结论:

1)本文方法的数值解与精确解和文献中的数值解吻合较好，可见，本文的非结构有限体积法和求解器GTEA适用于三维辐射换热的数值模拟。

2)不同的网格划分和立体角划分对计算结果有相应的影响，空间离散的网格并不是越密越好，而是要和立体角离散的个数相对应，即空间网格离散的个数要与立体角离散的个数同时增加，空间离散产生的假散射与角度离散产生的射线效应在一定程度上才会互相抵消，所以可以得到更高精度的数值解。

3)可将程序应用于柴油机燃烧室的辐射换热模拟，燃烧室缸壁的辐射热流量会受到缸内工质的吸收系数和温度的影响;无量纲辐射热流密度随着工质吸收系数的增大而增大，工质温度对中心部位的热流量影响较大，对边角处的影响较小。为后续考虑碳黑颗粒的辐射传热对柴油机工作性能的影响打下了基础。

［1］沈季胜，沈瑜铭，陈红岩，等.柴油机气缸内传热模型研究历程与展望［J］.内燃机工程，1999，20(1):72-77.SHEN Jisheng，SHEN Yuming，CHEN Hongyan，et al.Prospect and research process on radiated heat transfer models in cylinder of diesel engine［J］.Chinese Internal Combustion Engine Engineering，1999，20(1):72-77.

［2］严兆大，沈季胜，刘震涛，等.用蒙特卡洛法计算柴油机缸内多元热辐射［J］.内燃机学报，2000，18(2):133-136.YAN Zhaoda，SHEN Jisheng，LIU Zhentao，et al.Calculation of the radiation heat transfer in cylinder of diesel engine by the Monte-Carlo method［J］.Transactions of CSICE，2000，18(2):133-136.

［3］熊仕涛，陈国华.柴油机缸内传热计算［J］.内燃机学报，2001，19(3):215-218.XIONG Shitao，CHEN Guohua.Calculation of in-cylinder heat transfer of diesel engine［J］.Transactions of CSICE，2001，19(3):215-218.

［4］潘克煜，周龙保，杨中乐.直喷式柴油机气缸内辐射传热的试验研究［J］.内燃机工程，1996，17(4):59-67.PAN Keyu，ZHOU Longbao，YANG Zhonglu.Experimental investigation of in-cylinder heat radiation in diesel engine［J］.Chinese Internal Combustion Engine Engineering，1996，17(4):59-67.

［5］SAID A M，BUTTSWORTH D R，YUSAF T F.A review of radiation heat transfer measurement for diesel engines using the two-colour method［C］//Proceeding of 3rd International Conference on Energy and Environment.Malaysia，2009:202-207.

［6］高天，潘克煜，陈飞.直喷式柴油机缸内气体辐射传热的研究［J］.西安交通大学学报，2006，40(3):302-306.GAO Tian，PAN Keyu，Chen Fei.Investigation on in-cylinder radiation from gases in DI diesel engine［J］.Journal of Xi＇an Jiaotong University，2006，40(3):302-306.

［7］宋有生，聂宇宏.燃烧室型线对缸内燃烧和辐射传热影响的计算分析［J］.科学技术与工程，2008，8(3):760-763.SONG Yousheng，NIE Yuhong.Calculation investigation on influence of combustion chamber pr of ile to combustion process and radiative heat transfer［J］.Science Technology and Engineering，2008，8(3):760-763.

［8］姜晓光，聂宇宏，谢凯弘，等.基于Fluent的柴油机缸内辐射换热分析［J］.科学技术与工程，2008，8(15):4099-4103.JIANG Xiaoguang，NIE Yuhong，XIE Kaihong，et al.Diesel engine cylinder internal radiation heat transfer analysis based on Fluent［J］.Science Technology and Engineering，2008，8(15):4099-4103.

［9］刘健，吕继组，计时鸣.柴油机缸内辐射换热的多维数值模拟研究［J］.机械工程学报，2009，45(12):311-317.LIU Jian，LYU Jizu，JIShiming.Study onmulti-dimensional numerical simulation of in-cylinder radiation heat transfer of diesel engine［J］.Journal of Mechanical Engineering，2009，45(12):311-317.

［10］SADAT H，WANG C A，DEZ V L.Meshless method for solving coupled radiative and conductive heat transfer in complexmulti-dimensional geometries［J］.Applied Mathematics and Computation，2012，218:10211-10225.

［11］明平剑.基于非结构化网格气液两相流数值方法及并行计算研究与软件开发［D］.哈尔滨:哈尔滨工程大学，2008:16-72.MING Pingjian.Development of numerical modeling for gas-liquid two-phase flows based on unstructured grids and parallell computing［D］.Harbin:Harbin Engineering University，2008:16-72.

［12］BEAk SW，KIM M Y，KIM JS.Nonorthogonal finite-volume solutions of radiative heat transfer in a three-dimensional enclosure［J］.Numer Heat Transfer B，1998，34:419-437.

［13］MING Pingjian，ZHANG Wenping.Numerical simulation of natural convection and radiation heat transfer in 2D enclosure on hybrid grids［J］.Numerical Heat Transfer:Part B，2012，59:116-137.

［14］KIM C，KIM M Y，YU M J，et al.Unstructured polygonal finite volume solution of radiative heat transfer in a complex axisymmetric enclosure［J］.Numerical Heat Transfer:Part B.2010，57:227-239.

［15］MURTHY JY ，MATHUR SR.Radiative heat transfer in axisymmetric geometries using an unstructured finite-volume method［J］.Numer Heat Transfer B，1998，33:397-416.

［16］KIM M Y.Assessment of the axisymmetric radiative heat transfer in a cylindrical enclosure with the finite volume method［J］.Int JHeat Mass Transfer，2008，51:5144-5153.