基于元胞自動機的氣動光學光線追跡算法*

雒亮夏輝劉俊圣費家樂謝文科?

1)(中南大學物理與電子學院,長沙410083)

2)(武漢大學高等研究院,武漢430072)

(2020年4月10日收到;2020年6月16日收到修改稿)

1 引 言

氣動光學是研究穩態或非穩態氣體流場與在其中傳輸的光束相互作用及傳輸規律的學科[1,2].流場中的激波、邊界層/剪切層和尾流等大梯度、非均勻折射率場分布會使在其中傳輸的光束的波前產生畸變,這就是氣動光學效應[3,4].

對于折射率脈動是微小量(10–6)的流場,流場特征長度遠大于傳輸光波長,因此光束在流場中的傳輸可認為滿足傍軸近似條件,進而可認為光線在流場中的傳輸路徑近似為直線并忽略流場對振幅的衰減,只需考慮非均勻流場對相位的畸變效應.因此在這種情況下,可以對直線光路徑上各點的折射率直接積分獲取光程(optical path length,OPL).但是對于大脈動量(由非定常流動和湍流大尺度結構產生)流場[5,6],例如超聲速光學頭罩繞流流場,由于流場激波的存在,簡單的對折射率場沿直線路徑進行積分計算OPL會存在較大誤差.因此就需要光線追跡得到光線在流場中的完整傳輸路徑,綜合多條光線的信息得到光波穿過流場后的OPL值,進而得到光程差(optical path difference,OPD)、斯特列爾比(Strehl ratio,SR)和光學傳遞函數(optical transmission function,OTF)等光學量并據此分析波前畸變.

Montagnino[7]在1968年通過泰勒級數展開法數值求解了光線方程,為光在非均勻介質中傳輸的數值計算提供了基礎.近年來,一些光線追跡的新方法不斷涌現.Chang等[8]基于三維Snell定律提出了一種在三維氣動光學流場中進行光線追跡的方法,該方法能夠有效地計算氣動光學的波前參數.Xu等[9]開發了一種反向光線追跡的方法,該方法以探測器為追跡初始位置,以遠距離的目標為追跡末位置, 可明顯簡化傳統的氣動光學光線追跡計算.

本文提出了一種基于元胞自動機(cellular automata,CA)的氣動光學流場光線追跡方法.CA是結構簡單但具有復雜自組織特性的離散動力學系統,是包含大量相同組分在局部相互作用的復雜自然系統的數學模型. 用CA來模擬一個物理過程的優點在于避免了用微分方程作為控制方程,而直接通過制定變換規則來模擬非線性物理現象[10].CA模型中,空間被離散成許多元胞,這些元胞只在有限的狀態集中取值,元胞的狀態更新遵循一定的變換法則.CA目前已被廣泛應用于如流體流動[11]、模式識別[12,13]、相位變化[14]和人口動力學[15,16]等多學科的研究中.對于氣動光學光線追跡,光線在流場中傳輸可以看成是光子在滿足變換規則網格間的傳輸,因此如果能通過制定合理的光子移動規則,利用CA可實現氣動光學流場中光線追跡計算.

本文基于納米粒子示蹤平面激光散射技術(nanotracer-based planar laser scattering,NPLS)[17,18]實驗得到的超聲速二維剪切層流場和脫體渦(detached eddy simulation,DES)[19]模擬得到的含激波的超聲速光學頭罩繞流流場的二維流場,計算、對比了CA光線追跡算法與數值求解光線方程光線追跡法(numerical solving the ray equation,NSRE)得到的光線路徑以及對折射率沿路徑積分得到的OPL數據,通過OPL得到了OPD的曲線.結果表明,CA算法對二維隨機流場光線追跡非常有效,并具有相對于光線方程數值求解方法更高的計算效率.

2 理論基礎

任意折射率分布介質中的光線傳輸路徑可用光線方程進行描述

其中,s為光線傳輸路徑上的弧長,r為光線傳輸路徑的位置矢量,n為折射率,?n為折射率梯度.該方程僅當流場折射率滿足特定分布時才有解析解.對于折射率隨機變化的氣動光學流場,光線傳輸路徑的計算通常采用NSRE算法,該方法的基本思想是:對過定點的一條光線矢量,基于光線方程,依次計算出下一個光線矢量[20,21],進而得到流場區域內光線的路徑.該方法能有效地處理激波處的折射率突變,具有較高的精度[22,23].然而由于該方法在處理非均勻網格時首先需要進行插值運算,且每一次迭代都要重置初始值并求解非線性微分方程,因此存在計算量大的不足.本文所采用的數值求解光線方程的差分方法為三階龍格庫塔算法[24],該算法主要通過求解光線位置和方向同時改變的復雜差分方程來獲得光線完整路徑.

CA算法中,當光線傳輸至某元胞(網格)時,該元胞的狀態為1,反之則為0.鄰居類型采用的是摩爾型鄰居[25],當前時刻的元胞狀態取決于上一時刻及其鄰居元胞的狀態.元胞的坐標和折射率可以分別用該元胞的節點坐標(j,k)和折射率nj,k來表示.光子在元胞間移動遵循一定位置變換規則和方向變換規則.位置變換規則用來獲得光線矢量矢端點的位置信息和判斷本次迭代后是否跨越網格,如果跨越網格,則需要根據方向變換規則計算光線偏轉角,反之則保持方向進入下一次位置變換.可見,CA追跡算法與差分數值求解光線方程追跡法最大的不同在于:將光線位置的改變與光線方向的改變拆分為兩個相互獨立的過程.避免了光線位置改變和方向改變的同時計算,因此編程更容易實現且具有較高的效率.CA光線追跡算法的具體步驟如下:

Yaks and sheep are busy grazing, rain or snow is coming then.

1)給定初始光線矢量r0=(0,0).V1表示初始元胞(j0,k0)處歸一化光速,| V1|=1.qx,1和qy,1分別表示V1和x,y軸之間的夾角.

2)位置變換規則:對于第i次迭代,光線矢量ri=ri?1+?ri(i=1,2,3...),其 中?ri=Vi?t(?t=1為迭代時間步長)表示每次迭代光線矢量ri的改變量;1/nj,k),θx,i和θy,i分別表示Vi和x,y軸之間的夾角.rx,i和ry,i為 ri在x,y軸方向的分量.定義[rx,i]和[ry,i] 分 別表示對rx,i和ry,i向整數取整,并且定義rcx=[rx,i]?[rx,i?1]和rcy=[ry,i]?[ry,i?1] 來 判 斷每一次迭代后是否跨越網格,從而判斷Vi+1的大小和方向.第i次迭代后:

如(2)式所示,當rcx和rcy中有一個為1時,表明元胞路徑向該方向推進一個網格;當rcx和rcy同時為1時,表明元胞路徑向傾斜方向推進一個網格.上述兩種情況需要根據(5)式來計算偏角θx,i+1和θy,i+1,否則, Vi+1的大小和方向保持不變.

3)方向變換規則:認為s傍近元胞路徑L,因此(1)式可以寫為

進而(3)式可寫成x和y方向的分量的形式:

其中,Lx,i和Ly,i分別表示L在x和y方向的分量.對(4)式兩邊同時積分

其 中,?θx,i=θx,i+1?θx,i和?θy,i=θy,i+1?θy,i分別表示計算后得出的在x和y方向的偏轉角的改變量.(5)式中的?nj,k/?x=(nj+1,k?nj?1,k)/2hx,i和?nj,k/?y=(nj,k+1?nj,k?1)/2hy,i采用中心差分來計算,其中hx,i和hy,i表示x和y方向相鄰元胞的間距.對于CA算法,dx(dy)和Lx(Ly)之間的位置關系只有平行、垂直和成45°角,因此(5)式中的dx/dLx,i和 dy/dLy,i利用三角關系來計算;dnj,k/dLx,i和 dnj,k/dLy,i分別由 dnj,k/dLi,x=(nj+1,k?nj,k)/hx,i和 dnj,k/dLi,y=(nj,k+1?nj,k)/hy,i計算得到.

3 實驗與結果

3.1 超聲速混合層流場

本文基于NPLS技術獲得高分辨率超聲速混合層和邊界層二維密度場分布.NPLS系統主要由光源、示蹤粒子發生器、CCD、風洞、同步控制器和數據采集模塊等組成[26].超聲速混合層流場來源于超聲速混合層風洞中隔板上下馬赫數不同的超聲速氣流混合[27];超聲速邊界層流場來源于超聲速風洞靠管壁的壁面薄層.利用NPLS技術獲得超聲速剪切層的密度場分布的關鍵在于使示蹤粒子的散射光準確地描述密度場.風洞及NPLS系統的主要參數、完整結構及實物圖易仕和等[28]和Tian等[29]已有詳細的論述,這里不再贅述.

實驗測量的超聲速混合層流場的對流馬赫數為0.5,其中兩股來流的馬赫數分別為3.509和1.400,對應的來流速度分別是654.7和421.1 m/s,屬于中等可壓縮流場.NPLS圖像的像素尺寸為1431×281, 單像素分辨率h=0.16 mm, 入射光波長l =1064 nm, 對應的KGD= 2.195 ×10–4m3/kg,跨幀間隔為15μs,樣本總數為50幀.流向為x正向,光線傳輸方向為y正向.實驗測量的超聲速邊界層流場的馬赫數為3.0,對應的來流速度為622.5 m/s.NPLS圖像的像素尺寸為2048×1436,單像素分辨率為h=0.013 mm,入射光波長l=0.5μm,對應的KGD=2.2×10–4m3/kg,跨幀間隔為15μs,樣本總數為30幀.流向為x正向,光線傳輸方向為y正向.NPLS技術拍攝的某時刻超聲速混合層和邊界層圖像如圖1所示.

圖1 NPLS獲得的超聲速剪切層圖像(a)混合層;(b)邊界層Fig.1.The flow visualization results of supersonic shear layers obtained by NPLS:(a)Supersonic mixing layer;(b)supersonic boundary layer.

3.2 超聲速光學頭罩繞流流場

在本文中,依據模型尺寸參數,采用基于SST兩方程湍流模型的DES模擬來得到超聲速光學頭罩繞流流場密度場分布[30].DES模擬在近物面采用雷諾平均方法(Reynolds averaged Navier-Stokes,RANS),在其他區域采用大渦模擬方法(large eddy simulation,LES),兼具前者計算量小的優點和后者能分離湍流流動的優勢,DES湍流模型方程和RANS控制方程采用全耦合求解[31].圖2所示為無冷卻噴流光學頭罩繞流流場歸一化密度分布,相關參數取值為海拔高度10 km,攻角為0°,馬赫數為6.0.對密度做了無量綱處理,假設r表示真實密度,則圖中的無量綱密度值可由r/0.41270得到.其中仿真參數為光波波長1μm,初始位置位于均勻來流區域.在仿真計算時,設定光線傳輸方向為z方向且垂直于光學窗口入射.計算時,沿著平行和垂直于光學窗口方向取二維截面,尺寸為600 mm×50 mm.

圖2光學頭罩繞流流場密度場Fig.2.The density field distribution of supersonic flow field surrounding the optical dome obtained by DES.

4 計算結果與分析

在混合層流場密度分布變化較大的區域(圖1(a)中紅框區域),使用CA和NSRE算法獲得的光線路徑如圖3所示.CA算法和NSRE算法計算得到的出射該區域后的偏角分別為31.1和30.7 μrad.

在氣動光學中,折射率n和密度r之間的關系可表示為

其中,KGD為Gladstone-Dale系數.對光線沿著路徑r上的折射率積分得到OPL

由表達式

可得到光在穿過流場后的光程差.(8)式中?OPL?代表空間平均.OPD的計算結果如圖4所示.

圖3 CA算法與NSER算法在混合層得的光線路徑Fig.3.Beam paths obtained by CA and NSRE in mixing layer.

圖4 CA算法與NSRE算法計算得到的OPD(a)混合層;(b)邊界層;(c)含激波的超聲速光學頭罩二維剖面流場Fig.4.The OPD results calculated by CA and NSRE:(a)Supersonic mixing layer;(b)supersonic boundary layer;(c)supersonic flow field surrounding the optical dome.

根據統計學中均方根誤差的定義

其中M表示取樣點總個數,τ表示取樣點,OPDNSRE和 O PDCA分別表示NSRE算法和CA算法計算的OPD.超聲速混合層、超聲速邊界層以及超聲速光學頭罩二維剖面流場的OPDrms分別為0.0086 l1,0.0014 l2和0.0146 l3,其中l1,l2和l3分別對應超聲速混合層、超聲速邊界層以及超聲速光學頭罩二維剖面流場的波長.結果表明:對于超聲速混合層、超聲速邊界層以及超聲速光學頭罩二維剖面流場,CA計算得到的OPD結果與NSRE算法計算結果具有較好的一致性.

上述結果顯示了CA算法在氣動光學流場計算中的高精度以及適用性.然而衡量一種算法的優劣除了其計算結果的正確性,計算效率也是一個重要的評價指標.tic和toc函數是MATLAB中的內置函數,tic用來保存當前時間,toc用來記錄程序完成時間.兩個函數的配套使用可以計算算法的實際運行時間,從而評價算法的效率.對于本文的光線追跡,計算光線傳輸一個網格的程序運行時間,混合層、邊界層和高速繞流流場的程序執行時間分別如表1、表2和表3所列.表中數據均保留三位小數.表中,t1?t5表示每種流場采用CA或NSRE算法的計算時間,ta表示對計算時間求平均值,K=5表示計算總次數,?表示計算次數,ta可表示為

表1—表3中的計算時間為兩種算法程序在配置為Inter(R)Core(TM)i9-9900k CPU@

3.6 GHz(內存32 GHz)的電腦上運行得到的.由表1—表3可見,NSRE算法的程序平均執行時間約為CA算法的4倍.從計算模擬時間來看,在相同的條件下,CA算法的計算效率要高于NSRE算法.從算法結構上來看,在每一次迭代周期,CA算法的加法和乘法次數分別為11次和7次,NSRE算法的加法和乘法次數分別為39次和32次,NSRE算法的加法和乘法次數約為CA算法的4倍,這與表1—表3中的時間之間的關系符合得很好.

表1 CA與NSRE算法計算混合層流場的程序執行時間Table 1.The program running time of CA and NSRE in mixing layer.

表2 CA與NSRE算法計算邊界層流場的程序執行時間Table 2.The program running time of CA and NSRE in boundary layer.

表3 CA與NSRE算法計算高速繞流流場的程序執行時間Table 3.The program running time of CA and NSRE in supersonic flow field surrounding the optical dome.

5 結論

本文研究并對比了CA光線追跡算法和NSRE算法計算二維超聲速氣動光學流場OPD的結果.計算結果顯示,對于二維超聲速NPLS流場和二維超聲速繞流流場,CA算法的OPD計算結果與NSRE算法計算結果符合較好,具有很好的計算精度;從算法效率角度來看,CA的執行時間約為NSRE算法的1/4,具有比NSRE更高的效率.研究結果表明,CA光線追跡算法對于二維流場的氣動光學計算是適用的,這為離散非均勻流場的氣動光學計算提供了新方案.

感謝國防科技大學空天科學學院易仕和教授、中國空氣動力研究與發展中心陳勇研究員在NPLS實驗數據及流場CFD數據等方面提供的支持與幫助.