高速飞行器头罩气动热辐射传输效应计算

宋敏敏,王碧云,王 爽,吕 弢

(中国航天科技集团公司 红外探测技术研发中心,上海 201109)



高速飞行器头罩气动热辐射传输效应计算

宋敏敏,王碧云,王 爽,吕 弢

(中国航天科技集团公司 红外探测技术研发中心,上海 201109)

对头罩气动热辐射传输效应对红外探测系统的影响进行了研究。用四阶龙格-库塔法对头罩内传输的辐射光线进行追迹，根据理想光学系统成像特性对从头罩出射后在头罩后方理想光学系统中传输的空间辐射光线进行追迹。对头罩气动热辐射(温度场和光线)进行离散处理，建立了头罩内外表面辐射能量的计算模型，获得探测器各单元接收的头罩自身干扰辐射通量分布。仿真计算结果表明：头罩温度随飞行马赫数增大而上升，反之亦然；探测单元接收的头罩自身干扰辐射通量随飞行时间先增后减；气动热噪声随飞行时间先增后降，随头罩温度升高而增大。计算所得气动热噪声与产品实际采集的相符，方法正确。

高速飞行器； 光学头罩； 气动光学； 气动热； 光线追迹； 四阶龙格-库塔法； 头罩温度； 辐射通量； 气动热噪声

0 引言

随着精确制导技术的发展，精确打击武器的飞行速度已越来越快，同时对信息获取的精确性和时效性的要求也越来越高。飞行器的速度越快，气动热效应对头罩热应力和红外探测系统精度的影响就越大[1]。飞行器在大气层中高速飞行时，其光学头罩周围流场将产生真实气体效应、激波诱导边界层分离、无黏流与边界层的相互干扰等，从而引起气流密度、温度、组成成分等的变化，其中红外导弹头罩动热辐射效应对红外探测系统形成了辐射干扰，此即为气动热效应[2-4]。高速飞行器以超声速飞行时，头罩在气动热流的作用下温度迅速升高，高温使头罩自身热辐射增加，光学头罩表面的辐射主要集中于红外波段，这会对高速飞行器红外探测系统造成干扰，导致飞行器红外探测系统能力下降甚至探测系统饱和，使红外探测系统成像质量下降，制导精度也因此而降低。

目前关于气动热效应的研究多集中于气动光学头罩材料及形状等方面，研究气动热环境对头罩光传输影响的主要是相关高校，但其研究主要集中于光学头罩材料本身特性在气动热环境中红外辐射变化，未考虑气动热辐射在红外探测器中的响应[5-6]。同时传统气动热辐射计算方法不考虑头罩温度场随导弹不同飞行时刻的变化，一般只计算最高温和最低温两个状态，且计算时头罩表面只有一个温度，不细分整个头罩表面位置的温度场分布。本文采用的计算方法基于导弹头罩所用材料、具体结构设计等，充分考虑不同飞行弹道对头罩温度场分布的影响，通过建立热辐射传输模型可在飞行试验前分析计算不同飞行时刻气动热辐射噪声对红外探测的影响。在传输模型解算中，四阶龙格-库塔法具高精度，易收敛等优点，是仿真计算时常用的方法。针对本文研究的红外探测器接收的气动热辐射响应，先用CFD流体仿真软件计算不同飞行弹道的飞行器头罩在不同时刻的温度场分布，再用四阶龙格-库塔法对头罩热辐射传输的辐射光线进行追迹，建立传输模型，计算不同时刻因气动热而产生的头罩自身干扰辐射通量分布，同时综合探测器的响应特性(对应探测器的响应率、增益和偏置电压)，以获得探测器接受的气动热噪声。

1 头罩热辐射传输机理

气动热环境中高速飞行器光学头罩的温度场和热应变场分布并不均匀，由于热光效应和弹光效应的影响，光学头罩的折射率场分布亦不均匀。光学头罩可视作梯度折射率介质，而目前对梯度折射率介质内光线追迹的研究常采用欧拉法、泰勒级数展开法和龙格-库塔法。本文用精度最高的四阶龙格-库塔法对头罩内传输的空间辐射光线进行追迹。对从头罩出射的空间辐射光线传输，因头罩后方的光学系统为理想光学系统，故根据理想光学系统的成像特性对其进行追迹。

1.1 四阶龙格-库塔法

根据费马原理的拉格朗日描述，介质中的光线方程可表示为

(1)

式中：n为光线轨迹上某点的折射率；r为该点的位置矢量；ds为光线轨迹上的步长。式(1)在x、y、z轴方向分式为

对微分方程的初值问题，有

(2)

式中：[a,b]为取值闭区间，a为最小值，b为最大值；y0为x取a时y的初值。

由微分中值定理，取样点i处方程为

y(xi+1)-y(xi)=y′(ξ)(xi+1-xi).

(3)

式中：ξ∈[xi,xi+1]。令y′(ξ)=k，xi+1-xi=h，则式(3)可变为

y(xi+1)=y(xi)+k·h.

(4)

为获得较高精度的y(xi+1)，关键是确定合理的平均斜率值k。当k取值满足

k1=f(xi,yi)；

k4=f(xi+h,yi+hk3)

时，即为四阶龙格-库塔公式[7]。

1.2 头罩出射空间辐射光线追迹

从头罩出射的辐射光线为空间辐射光线。基于理想光学系统成像特性的空间辐射光线追迹的原理如图1所示。在头罩后方理想光学系统入瞳的中心取一根与出射空间辐射光线平行的辅助光线，根据理想光学系统成像特性，过入瞳中心的辅助光线在经理想光学系统传输后仍沿原传输方向传播，因此可根据计算直线与平面交点的公式得到辅助光线与理想像面的交点T。根据理想光学系统的成像特性，交点T同样也是入瞳处出射空间辐射光线在经过理想光学系统传输后与理想像面的交点。当头罩光学系统中心存在遮挡时，如头罩干扰辐射光线落在入瞳被遮挡的区域(图1中入瞳的ab区域)内，则不再对该光线进行光线追迹。

图1 基于理想光学系统成像特性的空间辐射光线追迹原理Fig.1 Ray tracing according to principle of ideal optical system imagine

2 头罩气动热辐射模型

2.1 头罩气动热辐射能离散

热辐射光谱是连续光谱，辐射源面上各点均向各方向发出连续波长范围的辐射，其辐射功率谱密度依赖于辐射源面的温度分布。辐射波面的法线可视作几何光学中的光线，光线携带辐射能向外传播。按光线的传播规律，能以有限能量的取样光线代表连续的辐射，以取样光线的能量因子作为其代表的一束辐射的能量。为准确地在探测器光敏面进行辐射噪声积分处理，先应完成一个辐射微分过程，即对头罩外内表面辐射线取样，包括辐射面元取样、辐射方向取样和波长间隔取样。

2.1.1 头罩温度场离散

为精确描述头罩内外表面温度场的非均匀分布情况，用不规则网格对头罩内外表面温度场进行划分(如图2所示)：圆周向分为Mφ份，高低向分为Mθ份，即对周向角和天顶角进行等分生成计算网格，有

(5)

式中：Mφ为周向角φ在[0，2π]内等分份数；Mθ为天顶角θ在[0，π/2]范围内等分份数。

图2 头罩温度场网格Fig.2 Temperature gridding of dome

设每个网格内温度为常量，可将光学头罩内外表面温度场离散成互不重叠的面元，看成由不同温度的面元组成，则头罩内外表面的热辐射可视作是各面元热辐射的叠加。

2.1.2 头罩热辐射能光线离散

将头罩的温度场按上述方法离散后，各离散面元辐射能的叠加即为头罩的干扰辐射[8]。对任一面辐射元来说，可将其视为朗伯辐射体，向2θ空间发出连续光谱的热辐射能，各方向的辐亮度相同。用几何光学分析，辐射波面的法线可视作光线，光线携带辐射能向外传播。因此，将任一面元发出的连续光谱热辐射能进行离散取样，当空间立体角及辐射波长的取样达到一定密度时，能以有限的离散光线替代面元向各方向发出的连续光谱热辐射能，取样后由该面元发出的具给定方向的一条光线，代表了此方向dΩm立体角内波长λ附近dλ范围内的一束辐射线簇。令辐射能量因子dW为该光线携带的能量，对每根离散的光线都有对应的方向矢量和能量因子表示其传输方向及携带的能量。

空间立体角的离散取样指将面元发出辐射的2π空间离散为互不重叠的立体角dΩm，如图3所示。

图3 立体角划分Fig.3 Solid angle partition

本文的离散取样方法为：φ(在如图3所示xoy平面内)在[0，2π]内等分为Mφ份；为使离散后每根光线代表的立体角大小均衡，θ(在如图3所示垂直于zos平面内)在[0，π/2]范围内按余弦等分为Mθ份，有

(6)

离散后每根光线对应的发射方向矢量F的方向余弦(Fx，Fy，Fz)可表示为

(7)

2.2 头罩外内表面辐射能量

设头罩外表面某一面元的温度为Ti，根据普朗克辐射定律，该面元光谱辐射亮度

(8)

式中：ε为头罩外表面的发射率；h为普朗克常数；c为光速；k0为波尔兹曼常数。

若该面元发出的辐射光线方向与该面元法线方向的夹角为θ，辐射面元面积为ds，λ处光谱辐射亮度为Lλ，则在立体角dΩ内dλ波段范围内该辐射光线的辐射功率

dW = Lλcos θdλdsdΩ.

(9)

用本文的模型对头罩温度场进行离散时，因面元曲率较小，可视为平面面元，则面元面积

ds = r2sin θΔθΔφ.

(10)

式中：r为面元所在球面半径。

根据本文的空间角离散方法，1根辐射光线代表的空间dΩ由其周围相邻4根光线确定，则

(φi + 1-φi)(cos θj + 1-cos θj).

(11)

当对辐射源面、热瞳和波长间隔以足够的密度取样时，可认为经取样由该面元发出的具给定方向的一条光线代表了此方向dΩ内λ附近dλ范围内的一束辐射线簇，这条光线携带的能量为dW，即为该光线的辐射能量因子。追迹各取样光线并在探测器光敏面上划分网格，叠加计算统计落在每一网格内各光线能量因子总和即为辐射噪声I，有

Lir2sin θjcos θjΔθΔφdλ.

(12)

3 气动热效应仿真

令探测器探元在光学系统轴向，设置飞行器飞行弹道，飞行高度10 km，整个弹道飞行时间36 s，初始时刻以0.7Ma飞行，18 s时马赫数达到最大(3.6Ma)，同时确定飞行器头罩和光学系统相关参数，如头罩材料、头罩内外半径等。

3.1 头罩光辐射通量计算

3.1.1 头罩温度场分布

基于CFD流体计算软件，根据弹道参数计算不同时刻头罩内外表面温度场分布[9-10]。其中：0～24 s飞行速度为升速，之后为降速。所得飞行时间18，36 s时内外表面温度场分别如图4、5所示。由图可知：头罩温度随飞行马赫数而变，马赫数升高时头罩温度变高，反之下降。

图4 内表面温度场Fig.4 Temperature field of inside surface

图5 外表面温度场Fig.5 Temperature field of outside surface

3.1.2 光辐射通量

根据各飞行时刻的温度场分布，由式(12)计算各时刻探测器探元接收的头罩自身干扰辐射通量，结果如图6所示。由图可知：随着飞行时间增加，各探测器单元接收的头罩自身干扰辐射通量先增后减。

图6 10 km弹道光辐射通量Fig.6 Radiation flux of 10 km trajectory

3.2 探测器探元响应

根据以上结果，由光线追迹可得探测器探元上的头罩干扰辐射通量φ，可得探测器探元的输出噪声电压

V=G·R·φ+VN.

(13)

式中：V为探元接收到的辐射通量经探测器响应后的输出噪声电压；G为前置放大器的增益；R为探元的响应率；φ为探元接收的头罩干扰噪声辐射通量；VN为探元的均方根噪声。

将探测器相应参数和各时刻的光辐射通量代入式(13)，所得不同飞行时间和头罩温度的气动热效应引起的热辐射噪声分别如图7、8所示。由图可知：气动热噪声随飞行时间先增后降，随头罩温度升高而增大。

图7 不同时间热辐射噪声Fig.7 Thermal radiation with time

图8 不同温度热辐射噪声Fig.8 Thermal radiation with temperature

4 结束语

本文对头罩气动热辐射传输效应计算方法进行了研究。基于CFD流体仿真软件计算得到相应弹道不同飞行时刻飞行器头罩的温度场分布，将飞行器头罩作为气动热辐射效应产生的辐射源，用四阶龙格-库塔法对头罩内传输的辐射光线进行追迹，根据理想光学系统成像特性对从头罩出射后在头罩后方理想光学系统中传输的空间辐射光线进行追迹，获得了探测器各单元接收的头罩自身干扰辐射通量分布。用仿真法计算得到了热辐射效应产生的热噪声。研究发现：头罩温度随飞行马赫数而变，当马赫数升高时头罩温度变高，反之下降；气动热噪声随飞行时间先升后降，随头罩温度升高而增加。理论计算的气动热噪声与相同弹道产品实际采集得到的热噪声基本相符，表明本文计算方法正确。同时数据处理发现：当头罩温度场分布均匀时，由此产生的热噪声经滤波后对红外探测并无大影响，但当头罩温度场分布不均匀时，热噪声就很难用滤波方法完全消除。用本文的方法能仿真计算不同飞行弹道导弹所受的气动热辐射噪声，可预测不同飞行时刻气动热噪声对红外探测的影响；可分析飞行器红外探测系统不同结构设计，头罩不同材料对气动热辐射效应的影响，可针对性地采用气动热效应校正方法，如头罩材料选取和结构设计改进等，以减弱气动热辐射效应。后续将基于本文方法研究临近空间飞行器因高速飞行产生的气动热效应，根据临近空间飞行器实际，分析本法的适用性和实用性。

[1] 刘锋， 唐晓峰. 超声速氟化镁红外头罩热强度失效判别准则研究[J]. 上海航天， 2014， 31(3)： 39-42+47.

[2] 殷兴良. 气动光学原理[M]. 北京: 中国宇航出版社, 2003: 194-196.

[3] FAY J A, RIDDLE F R. Theory of stagnation point heat transfer in dissociated air[J]. Journal of the Aeronautical Science, 1958, 25(2): 73-85.

[4] RAINE K W, PENFOLD A B. Instrument for measuring the thermal coefficient of refractive index in the infrared[J]. Journal of Physics E: Scientific Instruments, 1985, 18: 593-597.

[5] 吴琳, 房建成, 杨照华. 基于湍流涡模型的气动光学效应影响参数分析[J]. 红外与激光工程, 2007, 36(1): 97-101.

[6] 张义广, 杨军, 李小丹. 光学整流罩的热学性能评价工程试验[J]. 红外与激光工程, 2008, 37(S2): 556-559.

[7] 陈誌敏. 龙格-库塔法及其Mathematica实现[J]. 武汉工程职业技术学院学报, 2006, 18(2): 72-76.

[8] 范志刚, 肖昊苏, 李辉, 等. 气动热环境下高速飞行器光学头罩光传输分析[J]. 航空兵器, 2005(5): 59-64.

[9] 陈澄, 费锦东. 侧窗头罩高速层流流场光学传输效应数值模拟[J]. 红外与激光工程, 2005, 34(5): 548-552.

[10] 王发民, 沈月阳, 姚文秀. 高超声速升力体气动力气动热数值模拟[J]. 空气动力学学报, 2001, 19(4): 439-445.

Thermal Radiation Calculation of Optical Dome for High Speed Vehicle

SONG Min-min, WANG Bi-yun, WANG Shuang, LV Tao

(Infrared Detection Technology Research and Development Center, China Aerospace Science and Technology Corporation, Shanghai 201109, China)

The influence of aerodynamic thermal radiation effects on infrared detection system was studied in this paper. The fourth-order Runge-Kutta method was adopted for tracing the radiation ray transmitted in the dome. The space radiation ray transmitted in the ideal optical system on the rear of head cover after being exited from the dome was traced according to the characteristic of ideal optical system imaging. The aerodynamic thermal radiation of dome including temperature field and ray was discretized. The radiation computation models of the inside surface and outside surface of the dome were established. The distribution of the self-interference radiation flux of the dome received by each unit of detector was given out. The simulation results showed that the dome temperature would be higher as the Mach of the vehicle increasing, and vi'ce ver'sa; the self-interference flux of dome received by each unit of detector would increase firstly and then decrease during flying; the aerodynamic thermal noise would increase firstly and then decrease during flying and would increase as the dome temperature going higher. The calculated aerodynamic thermal noise was in accordance with data collected from the real product, which meant that the method proposed was correct.

High speed vehicle; Optical dome; Aerodynamic optics; Aerodynamic heating; Ray tracing; Fourth-order Runge-Kutta method; Dome temperature; Radiation flux; Aerodynamic thermal noise

2015-09-23；

2015-12-29

国家安全重大基础研究项目(973)资助(613271010204)

宋敏敏(1985—)，男，硕士，主要研究方向为红外探测技术。

1006-1630(2016)04-0050-06

TJ760.3； TK121

A

10.19328/j.cnki.1006-1630.2016.04.009