基于Code V的氣動(dòng)熱環(huán)境下頭罩的光傳輸仿真

史要濤，王玉雷，解 放，薛文慧，范志剛

(哈爾濱工業(yè)大學(xué) 空間光學(xué)工程研究中心，哈爾濱 150001)

0 引 言

在高速氣動(dòng)熱環(huán)境下，光學(xué)頭罩受到高速來(lái)流的壓力和熱作用，在熱力耦合下會(huì)產(chǎn)生無(wú)規(guī)則的溫度分布和表面變形。由于熱光效應(yīng)和彈光效應(yīng)，頭罩晶體的折射率會(huì)隨溫度場(chǎng)和應(yīng)力場(chǎng)產(chǎn)生相應(yīng)的無(wú)規(guī)則非均勻的分布。這一無(wú)規(guī)則非均勻的折射率分布勢(shì)必會(huì)改變光線的傳播方向和光程差，使得后繼光學(xué)成像系統(tǒng)的入瞳處的波面產(chǎn)生畸變，從而造成成像失真和像偏移。高速氣動(dòng)熱環(huán)境下的頭罩光傳輸如同外流場(chǎng)的氣動(dòng)光學(xué)效應(yīng)一樣，都是高速光學(xué)頭罩成像系統(tǒng)要研究的問(wèn)題。

近年來(lái)，對(duì)高速氣動(dòng)熱環(huán)境下的頭罩光傳輸和外流場(chǎng)的氣動(dòng)光學(xué)效應(yīng)的仿真研究引起了人們的關(guān)注［1］。肖昊蘇和范志剛利用Ansys 和Matlab仿真了高速熱環(huán)境下頭罩的熱響應(yīng)［2］，他們還建立了頭罩的光傳輸模型［3］。史可天等對(duì)頭罩繞流的氣動(dòng)光學(xué)效應(yīng)進(jìn)行了數(shù)值模擬［4］。Xu Liang 和Cai Yuanli 根據(jù)有限元分析得到的外流場(chǎng)，編寫程序研究了像偏移與海拔的關(guān)系［5］。Hao Chenglong和Chen Shouqian 使用有限元分析和四階龍格庫(kù)塔法對(duì)頭罩和外流場(chǎng)的綜合成像效應(yīng)進(jìn)行了研究［6］。雖然，對(duì)該問(wèn)題有了一定的研究成果，但以上研究使用的光學(xué)成像仿真方法都是自己編寫的，很難進(jìn)行廣泛交流和對(duì)比，且不能將后繼光學(xué)成像系統(tǒng)考慮在內(nèi)。本文利用廣泛推廣使用的商業(yè)化有限元分析軟件(Ansys)和光學(xué)仿真軟件(Code V)進(jìn)行高速氣動(dòng)熱環(huán)境下光學(xué)頭罩光學(xué)成像影響分析。該方法簡(jiǎn)單、靈活，使用范圍廣，易于交流對(duì)比。

1 光學(xué)頭罩的有限元分析

基于有限元分析和Code V 的氣動(dòng)光學(xué)頭罩的成像性能評(píng)估方案中，首先使用Ansys 對(duì)高速飛行器光學(xué)頭罩進(jìn)行有限元分析。有限元分析包括兩部分內(nèi)容:一是使用Ansys 中的計(jì)算流體力學(xué)分析軟件Fluent 仿真出光學(xué)頭罩的外流場(chǎng)，以期得到光學(xué)頭罩外表面的壓強(qiáng)分布、溫度分布和熱流密度分布;二是以流體力學(xué)分析的結(jié)果為輸入條件對(duì)光學(xué)頭罩進(jìn)行熱力耦合分析，以期得到頭罩的溫度分布、表面形變、應(yīng)力分布。

1.1 光學(xué)頭罩的有限元分析

對(duì)半球頭罩的高速飛行器進(jìn)行了案例分析。半球頭罩的材料為藍(lán)寶石晶體;外表面半徑為55.5 mm;內(nèi)表面半徑為49.5 mm;頭罩厚度為6 mm，如圖1 所示。根據(jù)高速飛行器的幾何模型，采用GAMBIT 網(wǎng)格劃分軟件建立外流場(chǎng)計(jì)算網(wǎng)格模型，網(wǎng)格采用非結(jié)構(gòu)體網(wǎng)格劃分方式，而且要保證頭罩表面附近邊界層內(nèi)的網(wǎng)格密度足夠大。

圖1 飛行器與頭罩示意圖

將網(wǎng)格計(jì)算模型導(dǎo)入Fluent 中，選擇密度基求解器;迭代格式選擇隱式格式;迭代方法選擇基于節(jié)點(diǎn)的高斯克林函數(shù)求梯度法;流量類型選AUSM;迭代松弛因子為0.6;湍流亞松弛因子為0.55;固體亞松弛因子為1。流場(chǎng)計(jì)算邊界條件如表1 所示。

頭罩繞流壁面熱流密度和靜壓分布如圖2 ～3所示。

表1 流場(chǎng)計(jì)算邊界條件

圖2 頭罩外表面熱流密度分布

圖3 頭罩外表面氣壓分布

1.2 球頭罩熱力耦合分析

建立頭罩的有限元模型，將流場(chǎng)分析得到的結(jié)果作為熱力耦合分析的載荷，頭罩的初始溫度和參考溫度為300 K;頭罩內(nèi)表面考慮自然對(duì)流，對(duì)流系數(shù)為10 W/(m2·K);內(nèi)表面的大氣壓為0.1 MPa;而頭罩底端面為固定約束;頭罩工作時(shí)間為15 s。選擇同時(shí)具有溫度和形變計(jì)算參量的非規(guī)則六面體單元，計(jì)算頭罩的溫度場(chǎng)、熱應(yīng)力場(chǎng)、熱形變場(chǎng)和熱應(yīng)變場(chǎng)，如圖4 ～7 所示。

圖4 15 s 時(shí)頭罩溫度場(chǎng)

圖5 15 s 時(shí)頭罩等效應(yīng)力場(chǎng)

圖6 15 s 時(shí)頭罩等效形變場(chǎng)

圖7 15 s 時(shí)頭罩等效應(yīng)變場(chǎng)

1.3 球頭罩非均勻折射率場(chǎng)

三維變折射率固體介質(zhì)折射率的變化量與溫度和熱應(yīng)力應(yīng)變等多種因素相關(guān)，且可以表示為

式中:Δnt為三維變折射率固體介質(zhì)由于溫度變化而引起的折射率變化量，這種現(xiàn)象被稱為熱光效應(yīng);Δnε為三維變折射率固體介質(zhì)由于熱應(yīng)變存在而引起的折射率變化量，這種現(xiàn)象被稱為彈光效應(yīng)。應(yīng)用頭罩的溫度分布和應(yīng)力應(yīng)變分布分別得到熱光效應(yīng)和彈光效應(yīng)所引起的折射率改變。

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

本文中提出的氣動(dòng)光學(xué)頭罩的成像性能評(píng)估方案的一個(gè)關(guān)鍵點(diǎn)在于使用現(xiàn)行廣泛使用的光學(xué)追跡軟件Code V。在復(fù)雜的高速氣動(dòng)熱環(huán)境下，光學(xué)頭罩的折射率因熱力耦合作用而無(wú)規(guī)則分布，一般很難用解析式擬合;同時(shí)，由于熱力耦合作用，頭罩的內(nèi)外表面產(chǎn)生無(wú)規(guī)則的形變，這種形變同樣很難用解析式表示。在光線追跡成像仿真中，合理精確地得到折射率分布和面型分布至關(guān)重要。

光學(xué)追跡軟件Code V 中提供自定義面型(UDS)和自定義梯度折射率(UDG)的功能給解決上述問(wèn)題帶來(lái)了方便。UDS 和UDG 是Code V 中自定義特征(UDF)的兩個(gè)應(yīng)用。該功能實(shí)現(xiàn)了交互設(shè)計(jì)和二次開(kāi)發(fā)。具體方法是根據(jù)Code V 提供的語(yǔ)法規(guī)則，使用C 語(yǔ)言或Fortran 語(yǔ)言編寫子程序，編譯成動(dòng)態(tài)鏈接庫(kù)文件(DLL)供Code V 調(diào)用。本文采用的是C 語(yǔ)言在VS2013 下編譯［7］。Code V 中自定義特征使用流程框圖如圖8 所示。

(1)自定義面型(UDS)

在Code V 中的光學(xué)表面定義如下:

式中:(x，y，z)為光線追跡點(diǎn)的坐標(biāo)。除了表面面型方程，光線追跡點(diǎn)的偏導(dǎo)也需提供。

圖8 Code V 中自定義特征使用流程框圖

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

在Code V 中，光線追跡點(diǎn)的折射率與坐標(biāo)有關(guān)，定義如下:

式中:(x，y，z)為光線追跡點(diǎn)的坐標(biāo)。除此之外，還需提供折射率沿坐標(biāo)軸的偏導(dǎo)。

需要注意的是，Ansys 有限元分析得到的氣動(dòng)頭罩的內(nèi)部折射率和內(nèi)外表面的數(shù)據(jù)是離散網(wǎng)格節(jié)點(diǎn)數(shù)據(jù)。這要求在UDS 和UDG 的應(yīng)用中，需要讀取離散的外部數(shù)據(jù)。仿真結(jié)果表明基于離散節(jié)點(diǎn)的UDS 和UDG 的結(jié)果與解析式的結(jié)果高度一致，證明了該方法的可行性［8］。

3 光學(xué)頭罩的面型和折射率數(shù)據(jù)插值處理

如上文所述，有限元分析得到的氣動(dòng)頭罩的內(nèi)部折射率和內(nèi)外表面的數(shù)據(jù)是離散網(wǎng)格節(jié)點(diǎn)數(shù)據(jù)，光線追跡點(diǎn)并不一定正好位于離散節(jié)點(diǎn)上，這就要求考慮離散數(shù)據(jù)插值問(wèn)題:光線追跡到頭罩表面是二維插值問(wèn)題，求取光線追跡點(diǎn)處的折射率屬于三維插值。

需要指出的是有限元分析得到的頭罩面型和折射率插值具有以下特點(diǎn)和要求:(1)插值點(diǎn)數(shù)據(jù)量巨大，面型數(shù)據(jù)量接近十萬(wàn)數(shù)量級(jí)，而三維折射率插值達(dá)到數(shù)十萬(wàn)數(shù)量級(jí)，插值算法應(yīng)是基于局部數(shù)據(jù)的;(2)頭罩面型和折射率插值屬于多維插值;(3)根據(jù)式(1)～(2)中面型和折射率定義的要求，所用的插值方法需能給出各變量方向的偏導(dǎo);(4)由于有限元分析所采用的網(wǎng)格是非結(jié)構(gòu)體網(wǎng)格，故所用的插值方法應(yīng)不依賴網(wǎng)格劃分。

基于上述插值特點(diǎn)和要求的考慮，本文采用空間插值中常用的改進(jìn)的謝別德插值［9］，與其他空間插值法相比，其有很大的運(yùn)算效率和精度優(yōu)勢(shì)。該方法是一種局部插值方法，可以很容易地向多維空間擴(kuò)展，且其插值函數(shù)的偏導(dǎo)一階連續(xù)。

4 計(jì)算仿真實(shí)例

根據(jù)有限元分析得到的數(shù)據(jù)，求出該球頭罩的折射率分布和內(nèi)外表面的變形，再結(jié)合Code V的自定義表面和梯度折射率的功能，仿真出氣動(dòng)熱環(huán)境不同視場(chǎng)下的導(dǎo)引頭光學(xué)系統(tǒng)的成像性能。

4.1 離散數(shù)據(jù)導(dǎo)入Code V

第15 s 時(shí)，彈光效應(yīng)導(dǎo)致的最大折射率改變量只占熱光效應(yīng)導(dǎo)致的最小折射率改變量的0.4%。故Code V 自定義梯度折射率時(shí)只考慮熱光效應(yīng)帶來(lái)的變化，晶體折射率隨溫度變化的關(guān)系為

式中:nt［λ，t(x，y，z)］為晶體在考察溫度下的相對(duì)折射率;t0為參考溫度;nt(λ，t0)為晶體在參考溫度下的相對(duì)折射率;dnt(λ，t)/dt 為晶體的相對(duì)折射率溫度系數(shù)，即熱光系數(shù);Δt(x，y，z)為溫度變化量;λ 為相應(yīng)色光的波長(zhǎng)。

由有限元分析得到氣動(dòng)熱環(huán)境下藍(lán)寶石晶體頭罩的溫度分布，結(jié)合熱光效應(yīng)，得到折射率分布如圖9 所示，折射率的最大值為1.701 1，最小值為1.699 7，差值為0.001 4，分布范圍并不大。

圖9 頭罩折射率分布

內(nèi)外表面的形變數(shù)據(jù)可以作為Code V 自定義面型的離散數(shù)據(jù)。內(nèi)外表面如圖10 ～11 所示，經(jīng)過(guò)15 s 工作時(shí)間后的最大等效形變約為47.7 μm，這一變化與頭罩的幾何參數(shù)相比是微小的。頭罩的等效形變場(chǎng)分布與頭罩溫度場(chǎng)分布一致，第15 s時(shí)頭罩具有最大等效形變分布的區(qū)域與頭罩具有最高溫度分布的區(qū)域一致。

導(dǎo)引頭由光學(xué)頭罩和光電成像系統(tǒng)組成。為突出考慮氣動(dòng)熱環(huán)境下頭罩的影響，在Code V 中用理想鏡頭模型代替光電成像系統(tǒng)。光學(xué)系統(tǒng)中，頭罩的外表面和內(nèi)表面，使用UDS 功能將離散數(shù)據(jù)導(dǎo)入Code V;而頭罩的折射率分布則是使用UDG將離散數(shù)據(jù)導(dǎo)入Code V;需要指出的是二維和三維改進(jìn)的謝別德插值法將在離散數(shù)據(jù)的導(dǎo)入中起重要作用。為考慮不同觀測(cè)角度［10］，理想鏡頭模型設(shè)為基本偏心，系統(tǒng)的光學(xué)結(jié)構(gòu)圖如圖12 所示。

圖10 頭罩內(nèi)表面

圖11 頭罩外表面

4.2 內(nèi)外表面形變像質(zhì)影響分析

為探究氣動(dòng)熱環(huán)境下球頭罩表面變形與折射率分布對(duì)導(dǎo)引頭光學(xué)系統(tǒng)的影響，分別對(duì)二者進(jìn)行仿真。只考慮內(nèi)外表面形變的影響，而折射率分布采用300 K 時(shí)的均勻折射率時(shí)的MTF 曲線，如圖13 所示。結(jié)果表明該氣動(dòng)熱環(huán)境下微小的內(nèi)外表面形變對(duì)導(dǎo)引頭光學(xué)系統(tǒng)的影響很小。

4.3 非均勻折射率分布像質(zhì)影響分析

不考慮球頭罩表面的形變，僅考慮非均勻折射率分布的影響。球頭罩折射率均勻分布時(shí)導(dǎo)引頭光學(xué)系統(tǒng)的MTF 和氣動(dòng)熱環(huán)境下不同觀察角度導(dǎo)引頭光學(xué)系統(tǒng)的MTF 曲線如圖14 所示。

圖14 非均勻折射率分布時(shí)系統(tǒng)MTF 曲線

圖14 表明氣動(dòng)熱環(huán)境下，球頭罩非均勻折射率分布將對(duì)導(dǎo)引頭光學(xué)系統(tǒng)的成像質(zhì)量帶來(lái)嚴(yán)重不利影響，即使是這種非均勻的變化量級(jí)并不大(最大最小差僅為0.001 4);對(duì)于光學(xué)系統(tǒng)入瞳位于球心的導(dǎo)引頭光學(xué)系統(tǒng)來(lái)說(shuō)，0°攻角氣動(dòng)熱環(huán)境帶來(lái)的成像失真與觀測(cè)角無(wú)關(guān)。需要指出的是，以上實(shí)例中的導(dǎo)引頭攻角為0°，光學(xué)頭罩是半球形且光電成像系統(tǒng)的入瞳中心在球心上，具有中心對(duì)稱性，這是氣動(dòng)熱環(huán)境帶來(lái)的成像失真與觀測(cè)角無(wú)關(guān)的原因;如不滿足這種中心對(duì)稱性(如二次曲面頭罩或入瞳中心在球心上)，則不能斷定氣動(dòng)熱環(huán)境帶來(lái)的成像失真與觀測(cè)角應(yīng)無(wú)關(guān)。

5 總 結(jié)

為提高高速飛行器導(dǎo)引頭的成像識(shí)別精度，探究氣動(dòng)熱環(huán)境對(duì)頭罩光學(xué)系統(tǒng)的影響，使用廣泛推廣的商業(yè)化有限元分析軟件Ansys 和光學(xué)設(shè)計(jì)軟件Code V 進(jìn)行高速飛行器在氣動(dòng)熱環(huán)境下的光學(xué)頭罩對(duì)導(dǎo)引頭光電成像系統(tǒng)的影響仿真分析。球頭罩0°攻角仿真結(jié)果表明氣動(dòng)熱環(huán)境帶來(lái)的非均勻的頭罩折射率分布是成像失真的主要因素，且成像失真與觀測(cè)角度無(wú)關(guān)，表明設(shè)法減少熱效應(yīng)是關(guān)鍵。本文介紹的方案簡(jiǎn)單、方便、通用且利于交流，為以后的研究提供了參考。

［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］范志剛，肖昊蘇，高豫強(qiáng).氣動(dòng)熱環(huán)境下高速飛行器光學(xué)頭罩特性分析［J］. 應(yīng)用光學(xué)，2009，30(3):361 -365.

［3］范志剛，肖昊蘇，李輝.氣動(dòng)光學(xué)頭罩光傳輸數(shù)值仿真［J］.應(yīng)用光學(xué)，2011，32(2):189 -194.

［4］史可天，馬漢東. 可壓縮混合層氣動(dòng)光學(xué)效應(yīng)研究［J］.計(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.