平板形、圆柱形、圆锥形石英灯阵列热流计算方法

张肖肖，张赐宝，赵旭升，丛琳华

中国飞机强度研究所，陕西 西安 710065

石英灯是地面热强度试验中常用的辐射加热元件，用于模拟高速飞行器飞行过程中的瞬态气动热效应[1]，进而开展相应热环境下的强度验证[2]。在加热器设计中需要根据试验件表面形状进行设计，如对于翼面、舵面等大面积平直表面的加热，一般设计为平板形加热器；对于等直舱段等圆柱形表面的加热，一般设计为圆柱形加热器；对于天线罩或变径非等直舱段等圆锥形或近似圆锥形表面的加热，一般设计为圆锥形加热器，如图1所示。

图1 加热阵列Fig.1 Heating array

试验件表面热流分布是评估加热器设计合理性的直接依据。杨晓宁等[3]建立了使用红外灯阵作为外热流模拟装置时到达被加热面热流密度的计算模型，对影响热流密度均匀性的因素进行了分析。刘守文等[4]采用蒙特卡罗法对红外灯能束辐射过程的物理模型进行了分析，建立了能束辐射随机位置、随机方向的数学模型和灯丝温度计算模型，给出了辐射波长、灯管反射率、吸收率、透射率的确定方法，进而提出了红外灯热流分布计算流程。同时，为了提高红外灯阵热流分布均匀度，应用单灯热流分布计算模型生成单灯热流分布数据库，解决了反光板反射热流问题，建立了红外灯阵热流分布仿真计算方法[5]。杨国巍等[6]利用蒙特卡罗法计算了红外单灯热流分布，并据此建立了数据库，综合考虑红外灯直射热流、挡板反射率以及间隙漏光量等影响因素，得到红外灯阵热流分布结果。王智勇等[7]基于蒙特卡罗法对单一石英灯加热器内部热流分布进行了数值模拟仿真，为更合理优化加热器内部热场分布提供了理论参考。夏吝时等[8]也采用蒙特卡罗法对平板形石英灯加热器中水冷反光板面积、水冷反光板与灯阵间距离、热源疏密程度、热源阵列与材料受热面间距等因素对辐射热场中典型隔热材料受热面温度分布均匀性和热流密度进行了模拟计算。朱言旦等[9]为了掌握石英灯单灯及灯阵热流分布规律，提高石英灯阵热流分布预测能力，对石英灯阵热流分布进行了计算分析与试验研究，发展了基于蒙特卡罗法的石英灯阵热流分布预测方法及计算程序。

以上研究对象多为平板形石英灯阵列，对圆柱、圆锥等热强度试验常见曲面形状的加热阵列热流计算方法研究报道较少。目前在加热器设计中受试验周期影响，往往依赖于经验选取加热阵列设计参数，缺少对表面热流分布状态的计算分析，而上述文献中广为采用的蒙特卡罗法在热流分布计算方面耗时较长，不能满足试验设计主管对于热流分布快速计算的需求。因此，需要建立一种能够满足平板形、圆柱形和圆锥形等常见形状阵列的辐射热流快速预测要求的计算方法，能对热强度试验中的加热器设计、试验件表面的热环境模拟快速给出直观的评估依据。

本文主要针对平板形、圆柱形和圆锥形的石英灯加热阵列，基于辐射传热学中的斯忒藩-玻耳兹曼定律和兰贝特定律，简化石英灯灯丝模型，忽略周围环境的折射、反射效应，将辐射加热元件离散成点光源，在试验件表面布置离散的网格节点，结合点光源与节点距离、加热元件和节点法线方向以及照射范围进行判断，分别计算各节点接收到各加热元件的辐射热流，之后将每个节点处的辐射热流进行累加，进而实现热试验常用石英灯阵列的热流分布快速计算。与蒙特卡罗法相比，在计算精度方面有所降低，在计算速度方面大幅提高。

1 理论基础

1.1 斯忒藩-玻耳兹曼定律

黑体的辐射力与热力学温度的关系由斯忒藩-玻耳兹曼定律所规定[10]

式中，σ为黑体辐射常数，其值为5.67×10-8W/(m2·K4)。

1.2 兰贝特定律

如图2 所示，用三维空间的立体角来表示某一方向的空间所占的大小[10]，定义为

图2 可见面积示意图Fig.2 Sketch map of visible area

设面积为dA的黑体微元面积向围绕空间纬度角θ方向的微元立体角dΩ内辐射出去的能量为dΦ(θ)，则试验测定表明

式中，I是从黑体单位可见面积发射出去的落到空间任意方向的单位立体角中的能量，称为定向辐射强度。其与黑体的辐射力之间的关系为

2 计算方法

考虑灯丝面向试验件表面的半圆柱段，将灯丝沿长度方向N等分，将每一半圆柱段视为平面微元处理，其向半球空间内辐射的总热流为

式中，L为灯管长度；r为灯丝半径；T为灯丝温度。

需要说明的是，第1节中所述公式虽针对黑体，在将其应用于灯丝时发射率的差异仅体现为灯丝温度的差异，因此为简化分析，本文中将灯丝发射率视为1。

2.1 平板形加热阵列

平板形加热阵列如图3所示，二维平面的法矢量a可写为[0，0，1]，对于平面上A点处面积为dA1的微元接收到某一微元dA（中心点为B）的辐射热流为

图3 平板形阵列点光源与节点位置关系示意图Fig.3 Sketch map of relation between light source point of plane array and receiving point

所以A点处受B点照射获得的热流密度为

将所有灯丝的微元对A点的热流密度累加，便得到了在此平板灯阵作用下A点的热流密度，按照此过程计算平面上所有布点获得的热流密度，便可以给出整个预测区域内的热流分布云图。

2.2 圆柱形外加热阵列

圆柱形外加热阵列如图4所示，对于圆柱面上接收点A点处的法矢量a可写为[Rcosα1，Rsinα1，0]，对于圆柱形阵列点光源B处的法矢量b可写为[-(R+d)cosα2，-(R+d)sinα2，0]，其中d为灯阵到圆柱面的距离，R为圆柱面半径。对于圆柱面上A点处面积为dA1的微元接收到面积为dA的微元B的辐射热流为

图4 圆柱阵列点光源与节点位置关系示意图Fig.4 Sketch map of relation between light source point of cylindrical array and receiving point

式中，θ1为点光源法线与AB连线的夹角；θ2为柱面上点的法线与AB连线的夹角。θ1与θ2二者的余弦值可通过矢量点积与矢量的模计算得到，如将A、B之间的矢量记为c，则

所以A点处受B点照射获得的热流密度为

需要注意的是，在圆柱灯阵中，每一根灯管的可照射范围是有限的，如在图4中，通过B点的切线BC、BD与圆的切点C、D之间的圆弧CD上的接收点可受到点光源B的照射。因此，在计算圆柱表面热流分布时，需要对表面节点A是否受到某个点光源B的照射进行判断，判断依据为灯管所在直线和接收点所在母线与圆柱圆心夹角小于相切时夹角（∠BOD），即

2.3 圆锥形外加热阵列

圆锥外加热阵列如图5所示，对于圆锥面上接收点A点处的法矢量a可写为[Rcosα1，Rsinα1，Rtanβ]，对于圆锥形阵列点光源B处的法矢量b可写为[-(R+d/cosβ)cosα2，-(R+d/cosβ)sinα2，-(R+d/cosβ)tanβ]，其中d为灯阵到圆锥面的距离，R为圆锥底面半径，β为半锥角，α1、α2分别为A点所在母线、B点所在阵列母线在xOy平面的投影与x轴的夹角。对于圆锥面上A点

图5 圆锥阵列点光源与节点位置关系示意图Fig.5 Sketch map of relation between light source point of conical array and receiving point

处面积为dA1的微元接收到面积为dA的微元B的辐射热流为

式中，θ1为点光源法线与AB连线的夹角；θ2为锥面上点的法线与AB连线的夹角。θ1与θ2二者的余弦值可通过矢量点积与矢量的模计算得到，如将A、B之间的矢量记为c，则

所以B点处受A点照射获得的热流密度为

在计算圆锥表面热流分布时，需要对表面节点是否受到某个点光源的照射进行判断。选定圆锥面上任一母线C'D'（C'为圆锥顶点，D'为圆锥下底面上一点），通过该母线做出圆锥面的切面，石英灯所在直线CD与C'D'平行，则该切面与除CD以外其他石英灯所在直线均产生一交点（如图6 中的E、F），可以看到各石英灯中位于该交点以上的部分可对母线C'D'上的点产生照射，而位于该交点以下的部分无法对母线C'D'上的点产生照射，而沿直线CD的石英灯各部分均可对母线C'D'上的点产生照射。如图7所示，圆锥面上母线C'D'处的法矢量d可写为（Rcosθ，Rsinθ，Rtanβ），θ为C'D'在xOy平面投影OD'与x轴的夹角。由于E点位于过母线C'D'的切面上，因此矢量d与由C'指向E的矢量e垂直，圆锥顶点C'坐标为（0，0，h），其中h为圆锥面高度。石英灯上点E坐标可写为（rcosθ3，rsinθ3，z1），其中z1为交点E的z向坐标，r为辐射加热元件所在圆锥体高度为z1处的切面圆半径，θ3为CG在xOy平面投影OE与x轴的夹角。矢量e可写为（rcosθ3，rsinθ3，z1-h），根据d⋅e= 0可得

图6 圆锥面辐射加热示意图Fig.6 Sketch map of radiation heating to conical surface

图7 圆锥面辐射加热几何模型Fig.7 Geometrical model of radiation heating to conical surface

对圆锥面有tanβ⋅h=R，代入式（17）整理可得

根据相似三角形特性，石英灯阵列所在圆锥面高度为z1处的截面圆半径r与z1的关系可表示为

式中，r1为石英灯阵列所在圆锥面的底面圆半径（如图7 中的OG），h1为石英灯阵列所在圆锥面的高度（如图7 中的OC）。根据解析几何知识，r1和h1分别可表示为

即对于投影与x轴夹角分别为θ的母线C'D'，各石英灯中位于z向坐标为z1的E点以上的部分可对母线C'D'上的点产生照射，而位于E点以下的部分无法对母线C'D'上的点产生照射，而与C'D'平行的直线CD的石英灯各部分均可对母线C'D'上的点产生照射。

2.4 内加热与外加热情况下的算法区别

圆柱形和圆锥形加热阵列如图8、图9 所示，与外加热相比，相当于是将灯管与试验件交换了方向，在照射范围判断上没有本质不同，都是判断光源与接收点能否在不受遮挡的情况下连出一条直线，所以只需在判断依据中将光源相关位置信息与接收点相关位置信息进行交换，就可以作为内加热情况的照射范围判断依据。需要注意的是，由于照射与接收的相对位置发生了变化，灯管点光源和试验件表面接收点的法矢量表达式也相应改变，见表1。

表1 内加热与外加热情况下法矢量表达式Table 1 Expression of normal vector of inner heating and outer heating

图8 圆柱内加热阵列Fig.8 Cylindrical inner array

图9 圆锥内加热阵列Fig.9 Conical inner array

2.5 具有解析表达式的曲面热流计算

对于形状复杂但具有解析表达式的空间三维曲面F(x,y,z)=0，其上一点(x,y,z)处的曲面法矢量为(Fx,Fy,Fz)，其中Fx、Fy、Fz分别为F(x,y,z)对x、y、z的偏导数。而石英灯阵列，受加热元件外形限制，所在平面或曲面的法矢量表达式不难确定，因此上述热流计算方法对于任意具有解析表达式的曲面也适用，均能得到类似式（11）的表达式，但需要根据具体曲面形状和阵列形状建立正确的照射范围判据。

3 热流计算结果

3.1 平板形加热阵列计算结果

取灯管长度为310mm，灯丝直径为3mm，灯丝温度为2500K，计算区域为350mm×300mm，灯管高度为50mm，灯管间距为30mm，阵列上下左右居中放置。图10 给出了1支和10支灯管组成的阵列在试验件表面的热流分布情况，图中虚线为灯管在计算区域平面上的投影。可以看到，单支灯管照射下试验件表面热流沿灯管中段向四周呈现由高到低的梯度分布；随着灯管排布密度的增加，加热阵列覆盖的试验件中部区域的热流密度接近同一水平，四周仍呈现出较为明显的热流梯度。

图10 平板形加热阵列热流分布结果Fig.10 Heat flux distribution of plane heating array

3.2 圆柱形外加热阵列计算结果

取灯管长度为310mm，灯丝直径为3mm，灯丝温度为2500K，圆柱试验件直径为400mm、高度为350mm，灯管距试验件表面高度为50mm，阵列上下居中放置。图11 给出了1支和20支灯管组成的阵列在试验件表面的热流分布情况。可以看到，单支灯管的照射范围较为有限，试验件表面热流呈现上下左右对称分布，并且具有显著的热流梯度；随着灯管排布密度的增加，加热阵列覆盖的试验件中部区域的热流密度趋于均匀，灯管正下方区域的热流偏高。

图11 圆柱形外加热阵列热流分布结果Fig.11 Heat flux distribution of cylindrical outer heating array

3.3 圆锥形外加热阵列计算结果

取灯管长度为310mm，灯丝直径为3mm，灯丝温度为2500K，试验件底面直径为150mm、高为500mm，灯管距圆锥试验件表面高度为30mm，阵列距试验件底面高度为50mm。图12 给出了1 支和9 支灯管组成的阵列在试验件表面的热流分布情况。可以看到，单支灯管照射下试验件表面热流呈现沿母线两侧左右对称分布，但受随高度渐变的截面圆曲率影响，沿母线上下并不对称，越接近圆锥底部，照射范围越大，但平均热流越小，这种趋势随着灯管排布密度的增加，可以观察得更为明显。

图12 圆锥形外加热阵列热流分布结果Fig.12 Heat flux distribution of conical outer heating array

4 结论

针对热强度试验中常见的平板形、圆柱形、圆锥形石英灯加热阵列，本文建立了三种内加热阵列的解析传热模型，提出了基于解析法的热流快速计算方法，分析了内加热与外加热情况下算法的区别，给出了典型灯管排布下热流分布结果，得出了以下结论:

（1）通过将平板形、圆柱形、圆锥形结构表面离散化以及石英灯管离散化，按照兰贝特定律逐点计算各接收点接收到灯管上各点光源的辐射热流，再计算各接收点接收到所有点光源的辐射热流之和，可以快速得到合理的结构表面热流分布结果。

（2）单支石英灯在圆柱面上照射范围有限，在热流分布计算过程中需要采用照射范围判断依据来判断圆柱面上母线与灯管所在直线的相对位置关系。

（3）不同石英灯对圆锥面上同一条母线的照射程度不同，在热流分布计算过程中需要采用照射范围判断依据来判断石英灯上光源点与圆锥面上母线之间的照射关系。