基于Code V的气动热环境下头罩的光传输仿真

史要涛，王玉雷，解 放，薛文慧，范志刚

(哈尔滨工业大学 空间光学工程研究中心，哈尔滨 150001)

0 引 言

在高速气动热环境下，光学头罩受到高速来流的压力和热作用，在热力耦合下会产生无规则的温度分布和表面变形。由于热光效应和弹光效应，头罩晶体的折射率会随温度场和应力场产生相应的无规则非均匀的分布。这一无规则非均匀的折射率分布势必会改变光线的传播方向和光程差，使得后继光学成像系统的入瞳处的波面产生畸变，从而造成成像失真和像偏移。高速气动热环境下的头罩光传输如同外流场的气动光学效应一样，都是高速光学头罩成像系统要研究的问题。

近年来，对高速气动热环境下的头罩光传输和外流场的气动光学效应的仿真研究引起了人们的关注［1］。肖昊苏和范志刚利用Ansys 和Matlab仿真了高速热环境下头罩的热响应［2］，他们还建立了头罩的光传输模型［3］。史可天等对头罩绕流的气动光学效应进行了数值模拟［4］。Xu Liang 和Cai Yuanli 根据有限元分析得到的外流场，编写程序研究了像偏移与海拔的关系［5］。Hao Chenglong和Chen Shouqian 使用有限元分析和四阶龙格库塔法对头罩和外流场的综合成像效应进行了研究［6］。虽然，对该问题有了一定的研究成果，但以上研究使用的光学成像仿真方法都是自己编写的，很难进行广泛交流和对比，且不能将后继光学成像系统考虑在内。本文利用广泛推广使用的商业化有限元分析软件(Ansys)和光学仿真软件(Code V)进行高速气动热环境下光学头罩光学成像影响分析。该方法简单、灵活，使用范围广，易于交流对比。

1 光学头罩的有限元分析

基于有限元分析和Code V 的气动光学头罩的成像性能评估方案中，首先使用Ansys 对高速飞行器光学头罩进行有限元分析。有限元分析包括两部分内容:一是使用Ansys 中的计算流体力学分析软件Fluent 仿真出光学头罩的外流场，以期得到光学头罩外表面的压强分布、温度分布和热流密度分布;二是以流体力学分析的结果为输入条件对光学头罩进行热力耦合分析，以期得到头罩的温度分布、表面形变、应力分布。

1.1 光学头罩的有限元分析

对半球头罩的高速飞行器进行了案例分析。半球头罩的材料为蓝宝石晶体;外表面半径为55.5 mm;内表面半径为49.5 mm;头罩厚度为6 mm，如图1 所示。根据高速飞行器的几何模型，采用GAMBIT 网格划分软件建立外流场计算网格模型，网格采用非结构体网格划分方式，而且要保证头罩表面附近边界层内的网格密度足够大。

图1 飞行器与头罩示意图

将网格计算模型导入Fluent 中，选择密度基求解器;迭代格式选择隐式格式;迭代方法选择基于节点的高斯克林函数求梯度法;流量类型选AUSM;迭代松弛因子为0.6;湍流亚松弛因子为0.55;固体亚松弛因子为1。流场计算边界条件如表1 所示。

头罩绕流壁面热流密度和静压分布如图2 ～3所示。

表1 流场计算边界条件

图2 头罩外表面热流密度分布

图3 头罩外表面气压分布

1.2 球头罩热力耦合分析

建立头罩的有限元模型，将流场分析得到的结果作为热力耦合分析的载荷，头罩的初始温度和参考温度为300 K;头罩内表面考虑自然对流，对流系数为10 W/(m2·K);内表面的大气压为0.1 MPa;而头罩底端面为固定约束;头罩工作时间为15 s。选择同时具有温度和形变计算参量的非规则六面体单元，计算头罩的温度场、热应力场、热形变场和热应变场，如图4 ～7 所示。

图4 15 s 时头罩温度场

图5 15 s 时头罩等效应力场

图6 15 s 时头罩等效形变场

图7 15 s 时头罩等效应变场

1.3 球头罩非均匀折射率场

三维变折射率固体介质折射率的变化量与温度和热应力应变等多种因素相关，且可以表示为

式中:Δnt为三维变折射率固体介质由于温度变化而引起的折射率变化量，这种现象被称为热光效应;Δnε为三维变折射率固体介质由于热应变存在而引起的折射率变化量，这种现象被称为弹光效应。应用头罩的温度分布和应力应变分布分别得到热光效应和弹光效应所引起的折射率改变。

2 Code V 自定义面型和自定义梯度折射率

本文中提出的气动光学头罩的成像性能评估方案的一个关键点在于使用现行广泛使用的光学追迹软件Code V。在复杂的高速气动热环境下，光学头罩的折射率因热力耦合作用而无规则分布，一般很难用解析式拟合;同时，由于热力耦合作用，头罩的内外表面产生无规则的形变，这种形变同样很难用解析式表示。在光线追迹成像仿真中，合理精确地得到折射率分布和面型分布至关重要。

光学追迹软件Code V 中提供自定义面型(UDS)和自定义梯度折射率(UDG)的功能给解决上述问题带来了方便。UDS 和UDG 是Code V 中自定义特征(UDF)的两个应用。该功能实现了交互设计和二次开发。具体方法是根据Code V 提供的语法规则，使用C 语言或Fortran 语言编写子程序，编译成动态链接库文件(DLL)供Code V 调用。本文采用的是C 语言在VS2013 下编译［7］。Code V 中自定义特征使用流程框图如图8 所示。

(1)自定义面型(UDS)

在Code V 中的光学表面定义如下:

式中:(x，y，z)为光线追迹点的坐标。除了表面面型方程，光线追迹点的偏导也需提供。

图8 Code V 中自定义特征使用流程框图

(2)自定义梯度折射率(UDG)

在Code V 中，光线追迹点的折射率与坐标有关，定义如下:

式中:(x，y，z)为光线追迹点的坐标。除此之外，还需提供折射率沿坐标轴的偏导。

需要注意的是，Ansys 有限元分析得到的气动头罩的内部折射率和内外表面的数据是离散网格节点数据。这要求在UDS 和UDG 的应用中，需要读取离散的外部数据。仿真结果表明基于离散节点的UDS 和UDG 的结果与解析式的结果高度一致，证明了该方法的可行性［8］。

3 光学头罩的面型和折射率数据插值处理

如上文所述，有限元分析得到的气动头罩的内部折射率和内外表面的数据是离散网格节点数据，光线追迹点并不一定正好位于离散节点上，这就要求考虑离散数据插值问题:光线追迹到头罩表面是二维插值问题，求取光线追迹点处的折射率属于三维插值。

需要指出的是有限元分析得到的头罩面型和折射率插值具有以下特点和要求:(1)插值点数据量巨大，面型数据量接近十万数量级，而三维折射率插值达到数十万数量级，插值算法应是基于局部数据的;(2)头罩面型和折射率插值属于多维插值;(3)根据式(1)～(2)中面型和折射率定义的要求，所用的插值方法需能给出各变量方向的偏导;(4)由于有限元分析所采用的网格是非结构体网格，故所用的插值方法应不依赖网格划分。

基于上述插值特点和要求的考虑，本文采用空间插值中常用的改进的谢别德插值［9］，与其他空间插值法相比，其有很大的运算效率和精度优势。该方法是一种局部插值方法，可以很容易地向多维空间扩展，且其插值函数的偏导一阶连续。

4 计算仿真实例

根据有限元分析得到的数据，求出该球头罩的折射率分布和内外表面的变形，再结合Code V的自定义表面和梯度折射率的功能，仿真出气动热环境不同视场下的导引头光学系统的成像性能。

4.1 离散数据导入Code V

第15 s 时，弹光效应导致的最大折射率改变量只占热光效应导致的最小折射率改变量的0.4%。故Code V 自定义梯度折射率时只考虑热光效应带来的变化，晶体折射率随温度变化的关系为

式中:nt［λ，t(x，y，z)］为晶体在考察温度下的相对折射率;t0为参考温度;nt(λ，t0)为晶体在参考温度下的相对折射率;dnt(λ，t)/dt 为晶体的相对折射率温度系数，即热光系数;Δt(x，y，z)为温度变化量;λ 为相应色光的波长。

由有限元分析得到气动热环境下蓝宝石晶体头罩的温度分布，结合热光效应，得到折射率分布如图9 所示，折射率的最大值为1.701 1，最小值为1.699 7，差值为0.001 4，分布范围并不大。

图9 头罩折射率分布

内外表面的形变数据可以作为Code V 自定义面型的离散数据。内外表面如图10 ～11 所示，经过15 s 工作时间后的最大等效形变约为47.7 μm，这一变化与头罩的几何参数相比是微小的。头罩的等效形变场分布与头罩温度场分布一致，第15 s时头罩具有最大等效形变分布的区域与头罩具有最高温度分布的区域一致。

导引头由光学头罩和光电成像系统组成。为突出考虑气动热环境下头罩的影响，在Code V 中用理想镜头模型代替光电成像系统。光学系统中，头罩的外表面和内表面，使用UDS 功能将离散数据导入Code V;而头罩的折射率分布则是使用UDG将离散数据导入Code V;需要指出的是二维和三维改进的谢别德插值法将在离散数据的导入中起重要作用。为考虑不同观测角度［10］，理想镜头模型设为基本偏心，系统的光学结构图如图12 所示。

图10 头罩内表面

图11 头罩外表面

4.2 内外表面形变像质影响分析

为探究气动热环境下球头罩表面变形与折射率分布对导引头光学系统的影响，分别对二者进行仿真。只考虑内外表面形变的影响，而折射率分布采用300 K 时的均匀折射率时的MTF 曲线，如图13 所示。结果表明该气动热环境下微小的内外表面形变对导引头光学系统的影响很小。

4.3 非均匀折射率分布像质影响分析

不考虑球头罩表面的形变，仅考虑非均匀折射率分布的影响。球头罩折射率均匀分布时导引头光学系统的MTF 和气动热环境下不同观察角度导引头光学系统的MTF 曲线如图14 所示。

图14 非均匀折射率分布时系统MTF 曲线

图14 表明气动热环境下，球头罩非均匀折射率分布将对导引头光学系统的成像质量带来严重不利影响，即使是这种非均匀的变化量级并不大(最大最小差仅为0.001 4);对于光学系统入瞳位于球心的导引头光学系统来说，0°攻角气动热环境带来的成像失真与观测角无关。需要指出的是，以上实例中的导引头攻角为0°，光学头罩是半球形且光电成像系统的入瞳中心在球心上，具有中心对称性，这是气动热环境带来的成像失真与观测角无关的原因;如不满足这种中心对称性(如二次曲面头罩或入瞳中心在球心上)，则不能断定气动热环境带来的成像失真与观测角应无关。

5 总 结

为提高高速飞行器导引头的成像识别精度，探究气动热环境对头罩光学系统的影响，使用广泛推广的商业化有限元分析软件Ansys 和光学设计软件Code V 进行高速飞行器在气动热环境下的光学头罩对导引头光电成像系统的影响仿真分析。球头罩0°攻角仿真结果表明气动热环境带来的非均匀的头罩折射率分布是成像失真的主要因素，且成像失真与观测角度无关，表明设法减少热效应是关键。本文介绍的方案简单、方便、通用且利于交流，为以后的研究提供了参考。

［1］Wang Meng，Mani A，Gordeyev S. Physics and Computation of Aero-Optics［J］. Annual Review of Fluid Mechanics，2012，44(1):299 -321.

［2］范志刚，肖昊苏，高豫强.气动热环境下高速飞行器光学头罩特性分析［J］. 应用光学，2009，30(3):361 -365.

［3］范志刚，肖昊苏，李辉.气动光学头罩光传输数值仿真［J］.应用光学，2011，32(2):189 -194.

［4］史可天，马汉东. 可压缩混合层气动光学效应研究［J］.计算物理，2010，27(1)，65 -72.

［5］Xu Liang，Cai Yuanli. Influence of Altitude on Aero-Optic Imaging Deviation［J］. Applied Optics，2011，50(18):2949 -2957.

［6］Hao Chenglong，Chen Shouqian，Zhang Wang，et al.Comprehensive Analysis of Imaging Quality Degradation of an Airborne Optical System for Aerodynamic Flow Field around the Optical Window［J］.Applied Optics，2013，52(33):7889 -7898.

［7］Optical Research Associates.Code V 10.2 Reference Manual［M］.California:ORA，1999.

［8］Frumker E，Pade O.Generic Method for Aero-Optic Evaluations［J］.Applied Optics，2004，43(16):3224 -3228.

［9］Basso K，de Ávila Zingano P R，Dal Sasso Freitas C M.Interpolation of Scattered Data:Investigating Alternatives for the Modified Shepard Method［C］∥XII Brazilian Symposium on Computer Graphics and Image Processing，IEEE，1999:39 -47.

［10］Zhang Wang，Chen Shouqian，Hao Chenglong，et al.Conformal Dome Aberration Correction with Gradient Index Optical Elements［J］. Optics Express，2014，22(3):3514 -3525.