雒亮夏輝劉俊圣費(fèi)家樂謝文科?

1)(中南大學(xué)物理與電子學(xué)院,長(zhǎng)沙410083)

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

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

1 引 言

氣動(dòng)光學(xué)是研究穩(wěn)態(tài)或非穩(wěn)態(tài)氣體流場(chǎng)與在其中傳輸?shù)墓馐嗷プ饔眉皞鬏斠?guī)律的學(xué)科[1,2].流場(chǎng)中的激波、邊界層/剪切層和尾流等大梯度、非均勻折射率場(chǎng)分布會(huì)使在其中傳輸?shù)墓馐牟ㄇ爱a(chǎn)生畸變,這就是氣動(dòng)光學(xué)效應(yīng)[3,4].

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

Montagnino[7]在1968年通過泰勒級(jí)數(shù)展開法數(shù)值求解了光線方程,為光在非均勻介質(zhì)中傳輸?shù)臄?shù)值計(jì)算提供了基礎(chǔ).近年來,一些光線追跡的新方法不斷涌現(xiàn).Chang等[8]基于三維Snell定律提出了一種在三維氣動(dòng)光學(xué)流場(chǎng)中進(jìn)行光線追跡的方法,該方法能夠有效地計(jì)算氣動(dòng)光學(xué)的波前參數(shù).Xu等[9]開發(fā)了一種反向光線追跡的方法,該方法以探測(cè)器為追跡初始位置,以遠(yuǎn)距離的目標(biāo)為追跡末位置, 可明顯簡(jiǎn)化傳統(tǒng)的氣動(dòng)光學(xué)光線追跡計(jì)算.

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

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

2 理論基礎(chǔ)

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

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

CA算法中,當(dāng)光線傳輸至某元胞(網(wǎng)格)時(shí),該元胞的狀態(tài)為1,反之則為0.鄰居類型采用的是摩爾型鄰居[25],當(dāng)前時(shí)刻的元胞狀態(tài)取決于上一時(shí)刻及其鄰居元胞的狀態(tài).元胞的坐標(biāo)和折射率可以分別用該元胞的節(jié)點(diǎn)坐標(biāo)(j,k)和折射率nj,k來表示.光子在元胞間移動(dòng)遵循一定位置變換規(guī)則和方向變換規(guī)則.位置變換規(guī)則用來獲得光線矢量矢端點(diǎn)的位置信息和判斷本次迭代后是否跨越網(wǎng)格,如果跨越網(wǎng)格,則需要根據(jù)方向變換規(guī)則計(jì)算光線偏轉(zhuǎn)角,反之則保持方向進(jìn)入下一次位置變換.可見,CA追跡算法與差分?jǐn)?shù)值求解光線方程追跡法最大的不同在于:將光線位置的改變與光線方向的改變拆分為兩個(gè)相互獨(dú)立的過程.避免了光線位置改變和方向改變的同時(shí)計(jì)算,因此編程更容易實(shí)現(xiàn)且具有較高的效率.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)位置變換規(guī)則:對(duì)于第i次迭代,光線矢量ri=ri?1+?ri(i=1,2,3...),其 中?ri=Vi?t(?t=1為迭代時(shí)間步長(zhǎng))表示每次迭代光線矢量ri的改變量;1/nj,k),θx,i和θy,i分別表示Vi和x,y軸之間的夾角.rx,i和ry,i為 ri在x,y軸方向的分量.定義[rx,i]和[ry,i] 分 別表示對(duì)rx,i和ry,i向整數(shù)取整,并且定義rcx=[rx,i]?[rx,i?1]和rcy=[ry,i]?[ry,i?1] 來 判 斷每一次迭代后是否跨越網(wǎng)格,從而判斷Vi+1的大小和方向.第i次迭代后:

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

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

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

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

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

3 實(shí)驗(yàn)與結(jié)果

3.1 超聲速混合層流場(chǎng)

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

實(shí)驗(yàn)測(cè)量的超聲速混合層流場(chǎng)的對(duì)流馬赫數(shù)為0.5,其中兩股來流的馬赫數(shù)分別為3.509和1.400,對(duì)應(yīng)的來流速度分別是654.7和421.1 m/s,屬于中等可壓縮流場(chǎng).NPLS圖像的像素尺寸為1431×281, 單像素分辨率h=0.16 mm, 入射光波長(zhǎng)l =1064 nm, 對(duì)應(yīng)的KGD= 2.195 ×10–4m3/kg,跨幀間隔為15μs,樣本總數(shù)為50幀.流向?yàn)閤正向,光線傳輸方向?yàn)閥正向.實(shí)驗(yàn)測(cè)量的超聲速邊界層流場(chǎng)的馬赫數(shù)為3.0,對(duì)應(yīng)的來流速度為622.5 m/s.NPLS圖像的像素尺寸為2048×1436,單像素分辨率為h=0.013 mm,入射光波長(zhǎng)l=0.5μm,對(duì)應(yīng)的KGD=2.2×10–4m3/kg,跨幀間隔為15μs,樣本總數(shù)為30幀.流向?yàn)閤正向,光線傳輸方向?yàn)閥正向.NPLS技術(shù)拍攝的某時(shí)刻超聲速混合層和邊界層圖像如圖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 超聲速光學(xué)頭罩繞流流場(chǎng)

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

圖2光學(xué)頭罩繞流流場(chǎng)密度場(chǎng)Fig.2.The density field distribution of supersonic flow field surrounding the optical dome obtained by DES.

4 計(jì)算結(jié)果與分析

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

在氣動(dòng)光學(xué)中,折射率n和密度r之間的關(guān)系可表示為

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

由表達(dá)式

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

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

圖4 CA算法與NSRE算法計(jì)算得到的OPD(a)混合層;(b)邊界層;(c)含激波的超聲速光學(xué)頭罩二維剖面流場(chǎng)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.

根據(jù)統(tǒng)計(jì)學(xué)中均方根誤差的定義

其中M表示取樣點(diǎn)總個(gè)數(shù),τ表示取樣點(diǎn),OPDNSRE和 O PDCA分別表示NSRE算法和CA算法計(jì)算的OPD.超聲速混合層、超聲速邊界層以及超聲速光學(xué)頭罩二維剖面流場(chǎng)的OPDrms分別為0.0086 l1,0.0014 l2和0.0146 l3,其中l(wèi)1,l2和l3分別對(duì)應(yīng)超聲速混合層、超聲速邊界層以及超聲速光學(xué)頭罩二維剖面流場(chǎng)的波長(zhǎng).結(jié)果表明:對(duì)于超聲速混合層、超聲速邊界層以及超聲速光學(xué)頭罩二維剖面流場(chǎng),CA計(jì)算得到的OPD結(jié)果與NSRE算法計(jì)算結(jié)果具有較好的一致性.

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

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

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

表1 CA與NSRE算法計(jì)算混合層流場(chǎng)的程序執(zhí)行時(shí)間Table 1.The program running time of CA and NSRE in mixing layer.

表2 CA與NSRE算法計(jì)算邊界層流場(chǎng)的程序執(zhí)行時(shí)間Table 2.The program running time of CA and NSRE in boundary layer.

表3 CA與NSRE算法計(jì)算高速繞流流場(chǎng)的程序執(zhí)行時(shí)間Table 3.The program running time of CA and NSRE in supersonic flow field surrounding the optical dome.

5 結(jié)論

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

感謝國(guó)防科技大學(xué)空天科學(xué)學(xué)院易仕和教授、中國(guó)空氣動(dòng)力研究與發(fā)展中心陳勇研究員在NPLS實(shí)驗(yàn)數(shù)據(jù)及流場(chǎng)CFD數(shù)據(jù)等方面提供的支持與幫助.