基于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.