典型光學窗口流場的氣動光學效應數值模擬

譚小童，許和勇，*，田仁治

(1.西北工業大學 翼型、葉柵空氣動力學國家級重點實驗室，西安 710072；2.代爾夫特理工大學 航空航天工程學院，荷蘭 代爾夫特 999025)

0 引言

近些年，隨著激光武器和超聲速飛行器的發展，氣動光學受到越來越多的關注。根據光線的傳輸路徑，氣動光學系統可以分為主動光學系統（激光發射器）和被動光學系統（紅外尋的器）。在主動光學系統中，氣動光學效應導致高能激光發射器發出的激光束能量減弱，失去攻擊力或聚焦偏離攻擊目標；在被動光學系統中，氣動光學效應導致光學接收器上的圖像變形或者光學信號減弱，使跟蹤方向發生偏差。無論是主動光學系統還是被動光學系統都分為近場“氣動光學”區域和遠場“大氣光學”區域。

迄今為止，人們已經對湍流邊界層[1-2]、自由剪切層[3-4]、分離剪切層和湍流尾跡[5-6]等幾種典型流動進行了廣泛研究。截至目前，廣大研究人員已經設計出許多氣動光學系統裝置，其中圓柱半球轉塔是目前激光武器的首選發射裝置[7]。針對圓柱半球轉塔的氣動光學，科研工作者對具有平坦或共形的光學窗口的轉塔進行了許多風洞實驗和數值研究[8-10]。美國開展了ALL、ABL、ATL、SHIELD 等一系列項目來驗證各種轉塔系統的氣動光學性能，并做了一些流動控制研究[11]。Wang 等[12]、Jumper 和Gordeyev[13]也對近年來的計算和實驗研究工作進行了總結。流動控制也被應用于減輕三維轉塔的流動分離和相關氣動光學畸變[14-16]。

進入21 世紀后，我國的氣動光學研究開始起步，2003 年殷興良[17]出版的《氣動光學原理》一書奠定了我國氣動光學的基礎。李桂春[18]在其著作《氣動光學》一書中對氣動光學原理和實驗測量方法進行了詳細論述。史可天等[19]對新近的一些計算氣動光學方法進行了歸納總結。孫喜萬等[20]基于WCNS-E-5 對斜激波、超聲速混合剪切層、亞聲速混合流進行了廣泛的數值模擬驗證研究，并系統回顧了近年來氣動光學研究的進展[21]。丁浩林等[2]對超聲速湍流邊界層的氣動光學效應開展了風洞實驗研究；此外，他們還開展了帶超聲速氣膜的高超聲速光學頭罩的氣動光學效應抑制實驗[22-23]。董航等[24]對圓柱半球轉塔模型進行了氣動光學效應時空特性分析。

為了方便研究，通常將繞流流場分為平均流場分量和脈動流場分量，因此氣動光學效應也分為平均流場效應和脈動流場效應。平均流場效應決定平移、傾斜、離焦以及像散等低階項，容易通過自適應光學系統矯正，但是其絕對值較大；脈動流場效應影響氣動光學效應中的高階項，使用自適應光學系統矯正很困難[25-26]。

目前，采用Zernike 多項式擬合波前畸變并結合自適應光學思想矯正低階項的研究方法還較少，而且將近場氣動光學和遠場衍射結合的研究方法也不多。本文主要基于DDES 方法和光線追跡算法計算3D 共形窗口轉塔和3D 平面窗口轉塔的氣動光學效應；采用Zernike 多項式擬合波前畸變，并結合自適應光學方法初步分析了兩種光學窗口在不同轉角下的氣動光學效應；基于光線方程和點擴散函數，模擬光線經過近場氣動光學畸變后的遠場自由衍射光強。

1 物理模型和數值計算方法

1.1 流場計算方法

為了提高湍流模型在轉塔后部流動分離區域的預測能力，本文采用延遲脫體渦模擬（DDES），該方法是將雷諾平均Navier-Stokes（RANS）模擬和大渦模擬（LES）相結合。Strelets[27]在1997 年提出的SSTDES 模型，即在壁面附近采用SST 模型求解RANS方程模擬小尺度湍流結構；在分離區域轉化為LES求解大尺度湍流結構，并引入亞格子模型?；搮^域的小尺度湍流。而RANS/LES 的轉換則是通過當地的湍流尺度和網格尺度的相互對應關系來確定。單一的LES 方法求解邊界層區域需要耗費大量的計算資源，而DES 方法只需要在邊界層區域的網格尺度達到RANS 的精度即可，因此可以大大地減少計算量。

Strelets 提出的SST-DES 模型的原理是，當RANS模型預測的湍流長度Lt大于局部網格間距的區域時，從SST 模式切換到LES 模式。在這種情況下，湍動能方程中計算耗散項所用的長度尺度由局部網格間距替代，具體的轉換原理如下：

在DES 的公式中選擇當地網格最大邊長 δ是因為模型在邊界層中應使用RANS 模式，最大長度是符合要求的最安全估計。F1是用于區分k-ω和k-ε方程的混合函數。但是SST-DES 方法產生的實際問題便是當地網格的最大邊長 δ小于邊界層厚度時，RANS 模式提前過渡到LES 模式，導致?；瘧Σ蛔悖∕SD），容易引起網格誘導分離問題（GIS）。所以，本文在SSTDES 的基礎上對其耗散項進行修正，得到SSTDDES 模型[28]，將耗散項修正為如下的形式：

該修正推遲了RANS 模式轉換到LES 模式，有效地避免了GIS 問題。

1.2 物理模型及邊界條件

為了研究轉塔附近不同流態對氣動光學效應的影響，本文選取兩種不同的光學窗口（平面光學窗口和共形光學窗口），并對比研究了兩者在不同轉角下的氣動光學效應。轉塔模型的幾何模型如圖1 所示，底部圓柱半徑和頂部半球的半徑R=0.152 4 m，圓柱高H=0.114 3 m，光學窗口的平面半徑r=0.137 2 m，光線的發射仰角EL=30°，而光線的轉角AL=0°、90°、180°，其中共形光學窗口模型在三個轉角下的流場計算模型相同。

圖1 轉塔幾何模型示意圖及光學窗口視角的定義Fig.1 Schematic diagram of the turret geometry model and visual angle definition of the optical window

圖2 展示了轉角為 0?時，平面光學窗口轉塔的流場計算域和邊界條件設置，其余工況的流場計算域和邊界條件與該工況相同。為了穩健地計算流場，將上游距離圓心O7.5R處設置為計算域入口，采用速度入口邊界條件；而將下游距離圓心O20R處設置為計算域出口，采用靜壓出口邊界條件；上表面及兩側面采用絕熱自由滑移邊界條件，而下表面和轉塔表面則采用絕熱無滑移邊界條件[29]。

圖2 轉塔的流場計算域及其邊界條件示意圖Fig.2 Computational domain and boundary conditions for the flow simulation of a turret

數值計算選取的馬赫數為Ma∞=0.4，雷諾數為Re∞=1.43×106，其余計算條件為標準海平面條件。圖3 展示了兩個模型的表面結構化面網格分布，為了保證壁面第一層網格y+≤1，在無滑移邊界條件上的第一層網格高度為1 .67×10?5m。平面光學窗口轉塔的總網格量為1 .4×107，而共形光學窗口轉塔的網格量為1 .35×107。

圖3 轉塔表面面網格分布Fig.3 Turret surface grid distribution

1.3 光線追跡方法

通過GladStone-Dale 關系式可知，折射率和空氣中的密度直接相關，其表達式如下：

式中：r=xi+y j+zk，表示光線流場中的位置，而n(r) 和 ρ(r)分別為該位置處的折射率和密度；KGD為GladStone-Dale系數，在可見光條件下，KGD近似取值2.27×10?4m3/kg。在光線光學中，光線沿路徑對折射率的積分定義為光程（optical path length,OPL），其表達式如下：

式中，（x′,y′）為光學窗口平面中的光學坐標，z′為光線的發射方向。而光線在不均勻流場中的傳播導致了各條光線光程的不同，被定義為光程差（optical path difference,OPD），其表達式如下：

其中〈OPL(x′,y′)〉表示 OPL(x′,y′)的空間平均。

根據波動方程的射線近似，可以得到光線在不均勻介質中的傳播規律即光線方程，該方程的表達式如下：

由于該方程是典型的二階常微分方程，所以運用四級四階Runge-Kutta 方法[30]求解該方程。

為了驗證本文追跡算法的準確性，本文對螺旋線進行光線追跡[30]，該光線的表達式如下：

其中n0和α為常系數，表征了徑向變折射率介質中折射率分布的變化情況。上述徑向變折射率分布的介質是有解析解的，其光線軌跡的解析解表達式為：

其中，(x0,y0,z0)為光線初射點的坐標，l0、p0、q0分別為初始光點的折射率在對應的x、y、z軸的光學方向余弦的乘積。圖4 展示了光線追跡結果和解析解的對比情況，追跡得到的光線和理想光線的軌跡基本重合，所以本文將使用該追跡算法模擬光線在流場中的傳播。

圖4 光線追跡算法驗證Fig.4 Verification of the ray tracing algorithm

對于實際的光在流場中的傳播，初始平面光波穿過流場后不再是平面，這種實際波面與理想波面之間的偏差稱為波前像差。波前像差可以由一系列多項式的線性組合來表示。通常描述波前像差的多項式為Zernike 多項式[25]。光學系統像差、大氣湍流像差等靜態和動態像差都可以用Zernike多項式來表示。Zernike 多項式每一項都有明確的物理意義，并且在單位圓內正交。關于Zernike 多項式的描述都在極坐標下進行。

為了驗證Zernike 多項式擬合波前的合理性，圖5 給出了波前擬合Zernike 多項式前36 項的圖像和原始波前圖像的對比，其中Nx和Ny代表網格點數。兩者的誤差在 1×10?10量級，相比波前OPD 在1×10?7量級，誤差在千分之一，該誤差可以忽略不計。

圖5 原始波前和Zernike 擬合后的OPD 云圖對比Fig.5 Comparison of OPD contours between the original wavefront and the Zernike polynomial fitted wave-front

在后續的研究中將使用OPD的均方根 OPDrms評估氣動光學效應。在自適應光學中，通?？梢韵揭?、傾斜、離焦和像散等項來減弱氣動光學效應。所以本文采用自適應光學的方法，利用波前 OPD擬合36 項Zernike 多項式，并去掉平移、傾斜、離焦和像散項來消除低階項對氣動光學效應的影響，并將去掉低階項后的OPD 的均方根表示為。

1.4 光波空間自由傳輸理論

波陣面穿過湍流尾跡后，波前相位畸變通常表現為一階或更高，而量化遠場畸變的方法是直接計算遠場光強圖。在氣動光學區域，不均勻的傳播介質使得光線傳輸產生了光程差（OPD(r)），考慮光的波動性，光程差將使得波陣面上產生相位差，導致相位畸變，其表達式如下：

式中，λ表示光線波長，?φ表示光線在z平面的相位差。流場引起的光波相位畸變會嚴重影響光波的復振幅分布，光線的光波復振幅為E，則有：

假設光線經過的遠場大氣光學區域介質均勻，則光線的傳播遵循標量波動方程[18]。由波動光學可知，在強光或受到強擾動的條件下，光波的獨立傳播原理將不再滿足，光的衍射等波動現象將嚴重影響光波的光強輻照度（optical intensity）分布。其中光強為復振幅模的平方，即：

2 計算結果及其分析

2.1 流場計算結果

為了驗證網格拓撲方法和CFD 計算方法的準確性，本文對1 000 個時間步長內的流場進行了統計平均，其在物理時間上對應0.01 s。同時，每隔10 個時間步長統計一次瞬時流場，總共統計100 個瞬時流場用于后續氣動光學效應的計算。本文選取了文獻[31]中共形窗口轉塔的流場實驗結果與本文數值模擬結果進行對比驗證。數值模擬條件與實驗條件相同，馬赫數為Ma∞=0.4，雷諾數為Re∞=1.43×106，其余計算條件為標準海平面條件。圖6 展示了本文計算的共形窗口轉塔頂部半球中心截面壓力系數和實驗值的對比。本文的數值模擬結果在附著流區域和實驗值吻合較好；而在分離流區域，壓力系數與實驗值還存在微小差距。這是因為本文采用的是DDES 算法，在壁面位置處為了減少計算量采用RANS 求解，在分離區對流場壁面的解析能力不足，導致了分離區位置處的壓力系數與實驗值有所差距；但是DDES 相較于LES 可以大大減少計算量、計算時間，同時本文數值模擬結果與實驗值的誤差在可接受范圍內，所以仍然采用該種方法。

圖6 共形窗口轉塔頂部半球中心截面的時均Cp 分布Fig.6 Time-averaged pressure coefficient distribution along the central plane of the conformal window turret dome

2.2 近場氣動光學效應

為了研究共形光學窗口和平面光學窗口對氣動光學效應的影響，本文分析了兩種光學窗口分別在0°、90°和180°三個轉角下的氣動光學效應。圖7（a）展示了光學網格示意圖，光學傳輸的坐標軸分別為（x′,y′,z′），光學網格尺寸為 0.9R×0.9R×4R，光學網格數為 95×95×151，光學網格的解析依賴于流場的解析結果。圖7（b）展示了光線的傳輸過程，光線從轉塔內部發射，經過近場氣動光學區域，得到畸變的波前分布，然后遠場衍射至接收平面。

圖7 光線在流場中的傳輸過程Fig.7 Transmission of light in the flow field

基于Gladstone-Dale 關系式（式（6）），圖8 展示了兩種光學窗口轉塔在不同轉角下的折射率分布云圖。在本文的計算中，光波波長選取 λ=0.75 μm，則Gladstone-Dale 系數KGD=2.27×10?4m3/kg。為了更加真實地模擬光線發射過程，本文將光線追跡的起始位置設置在轉塔內部，如圖8 所示，其中紅色矩形框表示光學網格。在流場區域內的光學網格使用當地的折射率插值，而在流場區域外的光學網格使用自由來流的折射率插值。為了考慮湍流脈動效應的影響，本文將計算100 個瞬時流場的OPD，再將其平均，并計算平均后的OPD 的均方根 OPDrms。

圖8 不同轉角下的中心截面瞬時折射率分布Fig.8 Instantaneous refraction distributions in the central plane for different angles (the optical grids are indicated by red rectangles)

圖9 展示了共形光學窗口和平面光學窗口在轉角為0°時的 OPDrms和隨著追跡距離的變化趨勢。圖9（a）表明隨著追跡距離的增加，兩種光學窗口的 OPDrms都是先增加到1.5R處，然后再緩慢下降。但是，共形光學窗口的 OPDrms始終大于平面窗口的OPDrms；且隨著追跡距離的增加，兩者的差距還在不斷增大。這是因為共形窗口本身的曲面形狀會對氣動光學效應有一定的影響，而平面窗口的形狀則對氣動光學效應沒有影響。圖9（b）則表明去掉平移、傾斜、離焦和像散等低階項后，兩種光學窗口的隨著追跡距離的增加而增加至1.5R位置，后保持不變，且相差不大。結合圖8（a）和圖8（d），光線發射路徑都經過了轉塔后部的分離區，由于在該仰角下，共形轉塔和平面轉塔的分離位置相近，所以后部分離區也較為相似，從而使得兩者的折射率分布差異較小。因此，在去掉低階項后兩者的氣動光學效應相近。

圖9 兩種光學窗口在0°轉角下的O PDrms和 Fig.9 O PDrms and of the two optical windows at the angle of 0°

圖10 展示了共形光學窗口和平面光學窗口在轉角為90°時的 OPDrms和隨著追跡距離的變化趨勢。圖10（a）表明隨著追跡距離的增加，平面光學窗口的 OPDrms隨著追跡距離的增加先增加至1.5R處，之后緩慢減??；而共形光學窗口隨著追跡距離的增加先增加后保持不變，但是整體上平面光學窗口的OPDrms大于共形窗口的 OPDrms。圖10（b）則表明去掉平移、傾斜、離焦和像散后，仍舊是平面窗口氣動光學效應更強。對比圖8（b）和圖8（c），平面窗口轉塔附近的折射率梯度明顯大于共形窗口轉塔的折射率變化率，所以在此轉角下，共形窗口的氣動光學效應更小。

圖10 兩種光學窗口在90°轉角下的O PDrms和 Fig.10 O PDrms and of the two optical windows at the angle of 90°

圖11 展示了共形光學窗口和平面光學窗口在轉角為180°時的 OPDrms和。圖11（a）表明兩種光學窗口的 OPDrms隨著追跡距離的增加一直增大。結合圖8（a～f），此時光線經過非均勻折射率的流場區域遠遠大于0°和90°轉角下的非均勻流場區域；而且該轉角下折射率梯度也更大，導致了該轉角下OPDrms比前兩個轉角下的 OPDrms都更大。此時共形光學窗口的 OPDrms大于平面窗口的 OPDrms。對比圖8（c）和圖8（f），在光學窗口附近共形窗口轉塔的折射率梯度大于平面窗口轉塔；加之共形窗口形狀本身對光線傳輸的影響，兩者共同作用導致了共形窗口的 OPDrms更大。在圖11（b）中，去掉平移、傾斜、離焦和像散等低階項后，兩種光學窗口的都比0°和90°轉角下的小。這是因為在180°轉角下都是附著流，脈動效應較小，對氣動光學效應的低階項影響較大，而對高階項影響較小。同時，去掉低階項后在該角度下平面窗口的則會更大。

圖11 兩種光學窗口在180°轉角下的O PDrms和 Fig.11 O PDrms and of the two optical windows at the angle of 180°

表1 不同轉角下追跡末端的高階項OPD'rmsTable 1 OPD'rms at the end of the tracing for different angles

2.3 遠場衍射結果

圖12（c～h）展示了兩種光學窗口在不同轉角下無自適應光學矯正時，光線遠場衍射30 000R距離（真空）下的光強分布云圖；并將其與光學窗口出射處光線（圖12（a））以及未畸變光線傳播相同距離的理論光強分布（圖12（b））進行了比較。在0°轉角下，平面光學窗口和共形光學窗口發射的光線都出現了光斑的位置偏移，但是平面光學窗口還有明顯的光斑分散；在90°轉角下，平面光學窗口和共形光學窗口發射的光線也都出現了光斑的位置偏移，共形光學窗口還出現了光斑的抖動現象；而在180°轉角下，兩種光學窗口都出現了顯著的光斑偏移。

圖12 兩種光學窗口不同轉角下（0°、90°、180°）的遠場衍射光強分布云圖Fig.12 Far-field diffraction intensity contours at the angles of 0°,90° and 180° for the two optical windows

圖13 展示了在與圖12 相同的衍射距離下，不同光學窗口的遠場光強分布。在0°轉角下，兩種光學窗口的光強峰值都出現了一定的位置偏移，但共形光學窗口的光強峰值大于平面光學窗口；在90°轉角下，平面光學窗口的光強峰值發生了較嚴重的位置偏移且光強峰值有所增大，而共形光學窗口則是出現了光強峰值的減??；在180°轉角下，兩種光學窗口的光強峰值都發生了大幅的位置偏移和光強數值的增加。值得注意的是，在轉角為90°和180°時，畸變波面遠場衍射后的光強大于未畸變波面遠場衍射的光強。根據波動光學中的惠更斯原理，畸變波面可能使得子波向一方偏折和聚集，導致光強峰值增加和位置偏移。

圖13 不同轉角下（0°、90°、180°）的遠場衍射光強曲線圖Fig.13 Far-field diffraction intensity profiles at the angle of 0°,90° and 180°

3 結論

在本文的研究中，無自適應光學矯正下，共形光學窗口和平面光學窗口在各個角度下的氣動光學效應各有優劣。然而在利用自適應光學方法去掉平移、傾斜、離焦和像散等低階項后，流場的脈動效應決定了其高階項的大小，共形窗口和平面光學窗口在0°轉角下發射的光線經過了相似的分離區域，導致了此時兩種光學窗口的氣動光學效應相近。而在90°和180°轉角下，平面光學窗口的平面形狀加劇了附近流場的脈動效應，導致了平面光學窗口的氣動光學效應遠大于共形光學窗口。在3 種不同轉角下，光線傳輸區域位于分離區時（0°）的氣動光學效應大于位于附著流區域時（90°、180°）的氣動光學效應。綜上，共形光學窗口發射的光線的傳輸性能優于平面光學窗口。

對兩種光學窗口發射的光線進行遠場衍射計算，結果表明經過近場畸變的光線遠場自由衍射后可能比未畸變的光線遠場自由衍射的光強峰值更大，且此時峰值位置偏移也更為嚴重。