計算氣動光學(xué)研究進展

2019-05-08 02:00:04史可天馬漢東

空氣動力學(xué)學(xué)報 2019年2期

史可天, 馬漢東

(中國航天空氣動力技術(shù)研究院, 北京 100074)

0 引言

帶有光學(xué)成像探測制導(dǎo)系統(tǒng)的飛行器在大氣層內(nèi)高速飛行時，其光學(xué)窗口與來流之間發(fā)生劇烈的相互作用，形成激波、邊界層、混合層等復(fù)雜流場，對光學(xué)成像探測系統(tǒng)造成圖像傳輸干擾，引起目標圖像偏移、抖動、模糊以及能量散布，這種效應(yīng)稱為氣動光學(xué)效應(yīng)[1-3]。氣動光學(xué)效應(yīng)的存在會嚴重降低光學(xué)成像探測系統(tǒng)的性能，已成為紅外成像制導(dǎo)等高速流動光學(xué)傳輸系統(tǒng)研制的關(guān)鍵技術(shù)之一。

計算氣動光學(xué)是采用數(shù)值計算的方法研究空氣動力流場對光波傳輸和光學(xué)成像影響及其校正的一門交叉學(xué)科[4-5]，包括流場數(shù)值計算和光學(xué)傳輸數(shù)值計算兩大方面。光學(xué)折射率n和流場密度ρ通過Gladstone-Dale關(guān)系式相聯(lián)系：

n=1+KGDρ(1)

其中KGD為Gladstone-Dale系數(shù)，其值與入射光波波長和流體特性相關(guān)。

計算氣動光學(xué)的發(fā)展主要經(jīng)歷了三個階段，一是早期的基于CFD簡化方法的統(tǒng)計計算，二是基于RANS流場的光學(xué)傳輸計算，三是基于LES/DNS瞬態(tài)流場的光學(xué)傳輸計算。三種方法相輔相成，前兩種方法可以快速地預(yù)測氣動光學(xué)效應(yīng)，第三種方法可以對前兩種方法進行校核和驗證，并可以細致分析氣動光學(xué)效應(yīng)，研究其形成機理。

基于CFD簡化方法的統(tǒng)計估算方法采用流場簡化模型計算方法獲得光學(xué)傳輸模擬所需的流場信息，然后用統(tǒng)計光學(xué)原理構(gòu)建脈動流場的光學(xué)畸變統(tǒng)計模型，以波面均方值、光學(xué)傳遞函數(shù)平均值作為氣動光學(xué)效應(yīng)參數(shù)，對光學(xué)畸變進行平均或統(tǒng)計描述。

基于RANS流場的光學(xué)傳輸計算方法針對流場平均密度分布，采用幾何光線追跡法和Fourier光學(xué)相結(jié)合的方法，計算平均流場引起的光學(xué)畸變。通過波面光程差分布、像面光強(點擴散函數(shù))、瞄視誤差、Strehl比、環(huán)圍能量曲線等多種氣動光學(xué)效應(yīng)參數(shù)對流場光學(xué)傳輸特性進行定量描述。由于RANS計算采用的Reynolds平均抹平了脈動運動的細節(jié)，所以難以直接反映擬序結(jié)構(gòu)運動對光學(xué)傳輸?shù)挠绊憽?/p>

基于LES/DNS流場的光學(xué)傳輸計算方法采用幾何光學(xué)/波動光學(xué)和Fourier光學(xué)相結(jié)合的方法，對流場瞬態(tài)擬序結(jié)構(gòu)的光學(xué)傳輸特性進行計算，能夠了解光學(xué)畸變和流場結(jié)構(gòu)特征變化的時間歷程，更有利于建立二者間的內(nèi)在關(guān)聯(lián)。

以上三種方法都是目前氣動光學(xué)研究正在應(yīng)用的計算方法，但后兩種方法更為多見。這三種方法涉及的計算流體力學(xué)和光學(xué)傳輸方法的關(guān)系，如表1所示。

表1 計算氣動光學(xué)涉及的流場和光傳輸方法Table 1 CFD and optical transmission method in the computational aero-optics

本文對這三種方法的內(nèi)涵、研究發(fā)展歷程以及新近的一些研究進展進行綜述，最后指出未來的研究重點和方向。

1 基于CFD簡化方法的光學(xué)統(tǒng)計估算

早在1952年Liepmann[6]就開展了氣動光學(xué)研究的開創(chuàng)性工作，他為了確定紋影系統(tǒng)應(yīng)用于高速流場時的極限靈敏度，應(yīng)用幾何光學(xué)原理推導(dǎo)了小孔徑光束通過湍流流場后，偏折角均方值〈θ2〉的計算公式：

〈(?v/?y)2〉Rv(|y-ζ|)dydζ(2)

其中：折射率脈動值v=KGDρ′/n0，Rv(y)表示折射率變化相關(guān)函數(shù)。這一計算公式首次建立了光學(xué)畸變參數(shù)與流場參數(shù)之間的統(tǒng)計關(guān)聯(lián)，并且還在計算中考慮了沿光線傳播方向變化的相關(guān)函數(shù)。Liepmann提出的光束偏折角均方值的計算公式獲得了實驗證實[7]，實驗測量結(jié)果與Liepmann的計算十分吻合。

Sutton[8]根據(jù)統(tǒng)計光學(xué)理論提出了近場光波相位均方值的計算模型，首次建立了大孔徑條件下波面畸變參數(shù)與流場參數(shù)之間的定量聯(lián)系：

在這一計算公式中同樣考慮了密度脈動相關(guān)尺度對光學(xué)畸變的影響。

在隨后的10多年中，大部分的氣動光學(xué)研究工作都圍繞Sutton統(tǒng)計模型開展，主要是基于這一統(tǒng)計模型，應(yīng)用直接或間接的方法計算近場光波相位均方值[9]。

這一時期的計算氣動光學(xué)研究，其所需的流場密度脈動信息主要是通過計算流體力學(xué)簡化模型方法獲得的，包括點渦法、無粘位勢流方法等[10-11]。

雖然Sutton統(tǒng)計模型在工程應(yīng)用中取得了一定的成功，但仍有學(xué)者對統(tǒng)計方法在高超聲速飛行器氣動光學(xué)效應(yīng)研究中的適應(yīng)性提出質(zhì)疑。Havaner[12]通過理論分析認為Sutton統(tǒng)計模型僅適用于完全氣體狀態(tài)下的二維亞聲速均勻各向同性湍流流動，對于具有多尺度渦結(jié)構(gòu)的高速復(fù)雜流場，如超聲速湍流混合層、兩相流等非均勻各向同性流動，并不適用。

殷興良[13]從光學(xué)傳遞函數(shù)基本定義出發(fā)，基于密度脈動的空間平穩(wěn)性假設(shè)，提出了三維超聲速流場均勻脈動的光學(xué)畸變統(tǒng)計計算模型：

其中密度脈動相關(guān)函數(shù)C(x,y,z)常采用von Karman模型、指數(shù)模型和Gauss模型等形式。Gauss模型：

史可天等[14]放棄了密度脈動的空間平穩(wěn)性假設(shè)，僅認為密度脈動隨機過程是Gauss隨機過程，發(fā)展了一種適用于三維超聲速流場非均勻脈動的光學(xué)畸變統(tǒng)計模型:

[Rρ(ξ+x,η+y,ζ+z,0,0-z)/2+

Rρ(ξ,η,ζ;0,0,z)/2-

Rρ(ξ,η,ζ;x,y,z)]dzdζ}dξdη(2)

這一計算模型可用于復(fù)雜外形光學(xué)窗口脈動流場的光學(xué)傳遞函數(shù)平均值計算。如果密度脈動分布具有空間平穩(wěn)性，式(6)中最外層積分號內(nèi)的表達式將與ξ和η無關(guān)，公式則可以退化為殷興良提出的統(tǒng)計計算模型公式(4)。

2 基于RANS流場的光學(xué)傳輸計算

RANS方法認為湍流的脈動運動是隨機和無規(guī)則的，通過湍流模型描述流場脈動運動對平均運動的作用[15-16]。RANS方法可以針對復(fù)雜外形飛行器的繞流流場進行計算，提供流場平均密度分布，并在一定精度上對脈動密度進行重構(gòu)。

基于RANS方法獲得流場平均密度后，可以采用幾何光線追跡法計算光瞳范圍內(nèi)每一條光線的傳輸路徑，由此得到入射平面光波通過流場后產(chǎn)生的光程差，獲得光瞳函數(shù)；然后采用Fourier光學(xué)方法計算畸變波面通過成像系統(tǒng)到達像平面的光強分布，即點擴散函數(shù)，再由點擴散函數(shù)計算瞄視誤差、Strehl比、環(huán)圍能量曲線等氣動光學(xué)效應(yīng)參數(shù)對平均流場引起的成像偏移、模糊等效應(yīng)進行定量描述[17]。

光線追跡法根據(jù)幾何光學(xué)原理、光線微分方程對光線在氣體介質(zhì)內(nèi)的傳輸路徑進行求解[18]：

首先給定光線傳輸路徑上一點處的位置向量和光線單位切向量，然后根據(jù)該點處的密度以及密度梯度，計算得到傳輸路徑下一點的位置向量和光線單位切向量。為了提高光線方程的求解精度，可將光線方程化為1階微分方程組，然后采用4步4階Runge-Kutta方法顯式推進[18]：

獲得每條光線傳播路徑后，即可求出光程和光瞳函數(shù)[19]:

Fourier光學(xué)方法對光瞳函數(shù)進行Fourier變換，得到像平面光波復(fù)振幅分布，進而計算出點擴散函數(shù)(即像平面光強分布)和光學(xué)傳遞函數(shù)[20-21]：

PSF(x′,y′)=|U(x′,y′)|2(11)

圖1給出光學(xué)傳輸效應(yīng)計算采用的坐標系。

圖1 光學(xué)傳輸效應(yīng)計算坐標系示意圖[17]Fig.1 The coordinate for the optical effects simulation[17]

RANS方法除了提供光學(xué)傳輸計算所需的平均密度，還可以根據(jù)平均密度分布應(yīng)用工程經(jīng)驗公式對脈動密度進行估算，以此作為計算輸入條件，應(yīng)用統(tǒng)計方法對脈動密度引起的光學(xué)畸變進行數(shù)值計算。

傳統(tǒng)的脈動密度預(yù)測模型基于混合長度理論[22]:

其中:C為加權(quán)常數(shù)，對亞聲速流動C=0.32，對超聲速流動C=0.57。

潘宏祿等[23-24]對預(yù)測模型中的加權(quán)常數(shù)C進行了修正。其采用高階精度的大渦模擬方法，對超聲速剪切湍流場進行計算，由流場瞬態(tài)密度分布統(tǒng)計得到脈動密度，給出了混合長度模型中系數(shù)C的最適取值，見圖2。在轉(zhuǎn)捩區(qū)內(nèi)，C在0.1～0.55間劇烈變化，而在充分發(fā)展湍流區(qū)內(nèi)C在0.45～0.42間變化，波動范圍較小，建議取值0.43。

圖2 平板邊界層模型系數(shù)分布[23]Fig.2 Coefficients for the plane boundary layer[23]

閆溟等[25]針對氣動光學(xué)效應(yīng)的RANS方法進行了分析，對常用的幾種湍流模型進行了評估，將RANS方程和折射率輸運方程聯(lián)合求解，對折射率脈動均方值g=〈n′2〉進行了預(yù)測：

圖3為閆溟等計算得到的湍流密度脈動重構(gòu)結(jié)果。

(a) Ma=1.72

(b) Ma=3.56

應(yīng)用幾何光學(xué)、Fourier光學(xué)方法，針對RANS平均流場數(shù)據(jù)可以計算得到平均流場的光學(xué)傳輸特性；應(yīng)用統(tǒng)計光學(xué)方法，針對密度脈動重構(gòu)數(shù)據(jù)，可以獲得脈動流場的光學(xué)傳輸特性。如果將平均流場和脈動流場分別看作空間不變的線性光學(xué)系統(tǒng)，就可以通過卷積方法將這二者的點擴散函數(shù)結(jié)合，獲得湍流流場氣動光學(xué)傳輸特性[26]：

PSF(x′,y′)=PSFM(x′,y′)*PSFT(x′,y′)(14)

其中:PSFM(x′,y′)為平均流場產(chǎn)生的點擴散函數(shù)，PSFT(x′,y′)為脈動流場產(chǎn)生的點擴散函數(shù)，*代表卷積。

在式(14)對光學(xué)畸變的計算中，根據(jù)紅外成像探測系統(tǒng)的成像積分時間，確定截止頻率，可以將密度脈動區(qū)分為低頻脈動和高頻脈動分別進行處理，低頻脈動引起成像抖動，高頻脈動引起成像模糊。截止頻率kFZ根據(jù)流場特征速度VT和內(nèi)時間尺度τin確定：

kFZ=VT·τin(15)

基于RANS方法的氣動光學(xué)效應(yīng)計算是氣動光學(xué)效應(yīng)工程預(yù)測中常用的方法，可以針對復(fù)雜外形光學(xué)窗口繞流流動進行光學(xué)傳輸效應(yīng)預(yù)測，定量描述平均流場和脈動流場引起的光學(xué)畸變，為光學(xué)系統(tǒng)設(shè)計提供依據(jù)。

陳澄等[27]、史可天等 [28-29]、陳勇等 [30]、張黎等 [31]、萬士正等 [32]針對不同形式的光學(xué)窗口(凹窗、凸窗及頭罩等)、不同飛行彈道進行了高速飛行器氣動光學(xué)效應(yīng)計算分析，對光學(xué)傳輸效應(yīng)的影響因素進行了分析，包括飛行高度、飛行Mach數(shù)、光線入射角、成像系統(tǒng)工作波長及成像積分時間等參數(shù)，在工程設(shè)計中取得了成功的應(yīng)用。

由于RANS計算采用的Reynolds平均方法抹平了脈動運動的細節(jié)，丟失了很多包含在脈動運動中的重要信息[33]，所以難以直接反映擬序結(jié)構(gòu)運動等湍流特性對光學(xué)傳輸?shù)挠绊?，在氣動光學(xué)效應(yīng)機理研究中有一定的局限性。

3 基于LES/DNS瞬態(tài)流場的光學(xué)傳輸計算

隨著計算流體力學(xué)的發(fā)展和計算能力的提高，大渦數(shù)值模擬(LES)和直接數(shù)值模擬(DNS)以其具備的氣動流場瞬態(tài)描述能力，為分析波前畸變和流場結(jié)構(gòu)變化的時間歷程，研究氣動光學(xué)效應(yīng)機理提供了重要的手段。

DNS雖然能直接求解流場各種尺度的渦結(jié)構(gòu)瞬時信息，但對計算機的運算速度和容量要求很高，在現(xiàn)有的計算條件下，大多是對低Reynolds數(shù)下的簡單外形湍流問題或局部湍流流動的模擬[34]。例如時間發(fā)展剪切層的觀察點隨流體運動，?；氖怯^察點附近流場的時間演化過程(圖4)，與真實的物理過程并不完全等價[35]。

圖4 時間發(fā)展剪切層渦結(jié)構(gòu)DNS結(jié)果[35]Fig.4 Turbulence vortex structure from temporal developing mixing layer DNS[35]

LES是介于RANS和DNS之間的一種數(shù)值模擬方法，它通過空間濾波將流場結(jié)構(gòu)分為大小兩種尺度，直接模擬尺度大于計算網(wǎng)格的流場結(jié)構(gòu)，以亞格子模型近似小尺度渦的耗散行為[36-37]。大渦模擬方法可以獲得流場擬序結(jié)構(gòu)分布(圖5)，為光學(xué)傳輸計算提供密度信息[38-39]。

圖5 可壓縮混合層LES獲得的湍流渦結(jié)構(gòu)[37]Fig.5 Turbulence vortex structure from compressible mixing layer LES[37]

獲得瞬態(tài)流場后，應(yīng)用幾何光學(xué)、Fourier光學(xué)方法即可獲得每一時刻流場的點擴散函數(shù)等光學(xué)畸變定量特性。由于DNS和LES方法獲得豐富的流場結(jié)構(gòu)信息，如果光學(xué)計算采用異于流場計算的獨立網(wǎng)格，需要確保流場與光場間數(shù)據(jù)轉(zhuǎn)換的密度信息精確度。

由于DNS和LES方法對流動細節(jié)的分辨能力遠強于RANS方法，可以反映一些較小尺度的流動結(jié)構(gòu)，光波在這些小尺度結(jié)構(gòu)的傳輸過程中，有可能表現(xiàn)出一定的波動效應(yīng)，需要采用波動光學(xué)的方法來進行高精度光學(xué)傳輸計算。在氣動流場的光學(xué)傳輸計算中，只要求解傍軸近似條件下的光波方程即可獲得足夠精確的結(jié)果[40]：

White等[41-42]針對三維可壓縮剪切層流場進行LES計算，然后比較了幾何光學(xué)和波動光學(xué)的計算精度，計算結(jié)果表明由兩種方法獲得的光程分布差別甚微，僅有0.03λ。White認為目前大多數(shù)光學(xué)傳輸計算中，高階精度的幾何光學(xué)方法，如4步4階Runge-Kutta推進的光線追跡法，可以滿足計算精度需求。

LES為氣動光學(xué)機理研究提供了有力的研究手段。基于LES的氣動光學(xué)機理研究目前主要集中在可壓縮混合層流動[43-44]、超聲速邊界層流動[45]、后臺階流動[46]等高速流動的光學(xué)傳輸效應(yīng)計算分析。這些流動雖然幾何外形簡單，卻具有可壓縮湍流的重要特征。對這些流動的光學(xué)傳輸計算分析，可以認識和了解大尺度擬序結(jié)構(gòu)的空間發(fā)展變化、激波/激波干擾、激波/邊界層干擾等現(xiàn)象對光學(xué)傳輸?shù)挠绊憽?/p>

通過對可壓縮湍流流場的氣動光學(xué)效應(yīng)計算分析，得到了一些光學(xué)畸變的影響規(guī)律。例如，混合層流動從層流到轉(zhuǎn)捩直至完全湍流全過程的光學(xué)傳輸計算表明，隨著轉(zhuǎn)捩的發(fā)生，流場引起嚴重的光學(xué)畸變，而轉(zhuǎn)捩完成以后，光學(xué)畸變程度有所減弱，但仍然嚴重于轉(zhuǎn)捩以前的情況(圖6)；亞聲速邊界層氣動光學(xué)效應(yīng)的計算發(fā)現(xiàn)，黏性底層和緩沖層對光學(xué)傳輸?shù)挠绊懞苄?，光學(xué)畸變主要是由對數(shù)區(qū)流動引起的。這些光學(xué)傳輸?shù)挠绊懸?guī)律為氣動光學(xué)效應(yīng)校正提供了理論依據(jù)。

圖6 剪切層波面畸變空間分布[44]Fig.6 Wavefront distortion in the mixing layer[44]

在開展典型湍流流動氣動光學(xué)效應(yīng)機理研究的同時，有學(xué)者采用LES方法對真實光學(xué)窗口外形進行流場數(shù)值模擬，對其瞬態(tài)光學(xué)傳輸效應(yīng)進行分析，獲得了更精確的高速飛行器光學(xué)系統(tǒng)氣動光學(xué)效應(yīng)數(shù)據(jù)。其中，最值得關(guān)注的是美國空軍氣動光學(xué)飛行實驗室(Airborne Aero-Optics Laboratory)凸臺激光系統(tǒng)的氣動光學(xué)效應(yīng)計算研究[47]。

在這一研究中，Mathews等[48]采用帶壁面模型的LES方法對凸臺繞流場進行模擬，計算了流場引起的激光束分散特性。他們首先通過與試驗數(shù)據(jù)的對比驗證，確認了波面畸變均方值計算的準確性(圖7)，然后對不同高度角和周向角開展了共計390個光線入射角的光學(xué)傳輸計算，分析了流場渦結(jié)構(gòu)對不同口徑激光束的影響規(guī)律，給出了激光束失效的入射角度范圍，圖8為其計算獲得的繞流流線分布。為了獲得精準的氣動光學(xué)效應(yīng)數(shù)據(jù)，他們甚至考慮了凸臺系統(tǒng)安裝縫隙對流場結(jié)構(gòu)和光學(xué)傳輸效應(yīng)的影響[49]。

圖7 光程差均方根值計算與試驗結(jié)果比較[48]Fig.7 Comparison of OPDrms between computation and wind tunnel measurement[48]

根據(jù)數(shù)值計算和試驗測量獲得的基礎(chǔ)數(shù)據(jù)，美國空軍研究實驗室采用流體控制和光學(xué)補償方法，研制了航空自適應(yīng)光束控制臺系統(tǒng)，使激光武器可以在高速移動平臺上從任意角度精確打擊敵方目標。如果沒有這一系統(tǒng)，由于機體周圍強烈的湍流會造成激光束的嚴重分散，激光武器只能對正前方的目標實施精確打擊。

圖8 繞流流線分布[48]Fig.8 Streamlines of the velocity field[48]

4 結(jié) 論

在高速流場氣動光學(xué)效應(yīng)數(shù)值模擬中，三種計算方法相互補充，基于CFD簡化方法的光學(xué)統(tǒng)計估算和基于RANS流場的光學(xué)傳輸計算具有較高的效率，可以快速提供流場的光學(xué)畸變數(shù)據(jù)，適用于復(fù)雜外形光學(xué)系統(tǒng)設(shè)計的氣動光學(xué)傳輸效應(yīng)分析，而基于LES/DNS瞬態(tài)流場的光學(xué)傳輸計算雖然對計算資源要求較高，卻可以提供更精細的光學(xué)畸變信息。通過對三種方法的總結(jié)，得到以下結(jié)論：

(1)基于CFD簡化方法的光學(xué)統(tǒng)計估算方法通過近幾年的研究，從細光束偏折估算發(fā)展到大孔徑平面光波相位均方值、傳遞函數(shù)均值計算，適用范圍和計算精度都獲得了大幅提升。但是，計算過程中采用的密度相關(guān)函數(shù)在各向異性復(fù)雜非均勻流動條件下還有待進一步發(fā)展和驗證。

(2)基于RANS流場的氣動光學(xué)效應(yīng)計算方法可以針對復(fù)雜外形飛行器的繞流流場進行氣動光學(xué)計算，提供平均流場的成像偏移和模糊特性；通過脈動流場重構(gòu)方法，還可以對脈動密度引起的成像模糊進行評估，是氣動光學(xué)效應(yīng)應(yīng)用研究中常用的一種方法。由于光學(xué)傳輸計算需要考慮光波通過流場全路徑上的密度及其脈動的影響，所以如何精準描述流場的空間分布是RANS方法應(yīng)用于氣動光學(xué)效應(yīng)計算時有待解決的問題。此外，密度脈動重構(gòu)方法的精度還需要通過風(fēng)洞試驗進行確認。

(3)基于LES/DNS流場的光學(xué)傳輸計算方法能夠提供光學(xué)畸變的時間變化歷程，更有利于分析其與流場特性間的關(guān)聯(lián)。目前DNS方法受到計算條件限制，適用于局部湍流流動的氣動光學(xué)效應(yīng)機理分析。而LES方法的計算量相對DNS較小，同時其獲得的大尺度結(jié)構(gòu)又正是造成光學(xué)相位畸變的主要原因，因此LES方法在當(dāng)前及未來相當(dāng)長時期內(nèi)，都將是計算氣動光學(xué)機理和應(yīng)用分析的有力工具。為了能提供更精確的流場信息，LES方法在高精度數(shù)值格式、亞格子模型適用性等方面還需進一步研究。