謝文科 劉俊圣 費(fèi)家樂 周全 夏輝? 陳欣 張盼 彭一鳴 于濤

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

2) (國防科技大學(xué)前沿交叉學(xué)科學(xué)院,長沙 410073)

氣動光學(xué)的研究中,關(guān)聯(lián)方程(linking equation)是關(guān)聯(lián)湍流力學(xué)量與光學(xué)量的一個重要的方程.但是,基于模型簡化的關(guān)聯(lián)方程在亞聲速低速流場的應(yīng)用中通常忽略權(quán)重函數(shù)對波前方差估計(jì)精度的影響.本文在納米粒子示蹤平面激光散射技術(shù)獲得超聲速混合層流場密度數(shù)據(jù)的基礎(chǔ)上,應(yīng)用關(guān)聯(lián)方程計(jì)算超聲速混合層的流向波前方差,并進(jìn)行誤差分析.結(jié)果表明：基于關(guān)聯(lián)方程估計(jì)的流向波前方差與直接對密度場的積分計(jì)算結(jié)果具有較好的一致性;在適當(dāng)?shù)囟x相干長度、密度脈動協(xié)方差高斯模型近似的基礎(chǔ)上,分析了權(quán)重函數(shù)對關(guān)聯(lián)方程計(jì)算精度的影響,指出了權(quán)重函數(shù)對關(guān)聯(lián)方程在超聲速流場密度高度相關(guān)區(qū)域中應(yīng)用的必要性.研究的開展對于拓展關(guān)聯(lián)方程在高速流場中的應(yīng)用具有一定指導(dǎo)意義.

1 引 言

氣動光學(xué)是研究光波與窗口附近、非均勻薄層湍流流場相互作用的學(xué)科.由于湍流的隨機(jī)性,常用出射光束波前的平均值、方差、相關(guān)函數(shù)和相干長度等統(tǒng)計(jì)量來表征氣動光學(xué)效應(yīng)[1,2].由于波前方差可用來對上述其他統(tǒng)計(jì)量進(jìn)行建模,同時(shí)波前方差值還可以用來計(jì)算光學(xué)傳遞函數(shù)、斯特爾比等,因此波前方差是氣動光學(xué)統(tǒng)計(jì)建模首先要測量的量之一[3].早期的波前方差測量方法有直接光學(xué)測量和間接流體測量兩種,前者包含陰影法、紋影法、光學(xué)干涉法和哈特曼波前測量等.顯然,直接光學(xué)測量方法存在測試環(huán)境要求高、成本高等劣勢;后者是指Sutton[4]于1985年提出的基于關(guān)聯(lián)方程的間接測量方法,間接方法大大降低了測量環(huán)境要求和成本,因而在氣動光學(xué)研究歷史上具有重要意義.

Hugo和Jumper[5]研究了關(guān)聯(lián)方程在噴管出口流速為7 m/s的二維熱射流流場中的適用性,指出了流場特征長度的合理定義對用關(guān)聯(lián)方程計(jì)算精度的影響.Tromuer等[6]基于大渦模擬的流場數(shù)據(jù),討論了在馬赫數(shù)為0.9的邊界層流場應(yīng)用關(guān)聯(lián)方程計(jì)算波前方差的適用性.Fitzgerald等[7]發(fā)展了基于低速流氣動光學(xué)波前方差實(shí)驗(yàn)數(shù)據(jù)外推到高速流波前方差的工程模型.Yin等[8]基于關(guān)聯(lián)方程對亞音速機(jī)載氣動光學(xué)效應(yīng)進(jìn)行了研究.由于實(shí)驗(yàn)技術(shù)的限制,基于超聲速流場和高超聲速流場氣動光學(xué)效應(yīng)的研究較少[9].因此,開展關(guān)聯(lián)方程在高速復(fù)雜流場中應(yīng)用的適應(yīng)性研究和誤差分析工作具有現(xiàn)實(shí)理論和工程意義.

本文基于納米粒子示蹤平面激光散射(nanoparticle-based planar laser scattering,NPLS)技術(shù)實(shí)驗(yàn)測量了超聲速混合層密度場分布,計(jì)算、對比了基于關(guān)聯(lián)方程與直接對密度場積分的方法估計(jì)波前相位方差和誤差.研究結(jié)果表明,基于關(guān)聯(lián)方程的波前方差估計(jì)對于超聲速混合層流場同樣適用,估計(jì)誤差主要來源于在使用簡化了的關(guān)聯(lián)方程計(jì)算時(shí),是否考慮了權(quán)重函數(shù)的影響.

2 理論基礎(chǔ)

由于氣動光學(xué)流場的折射率脈動是微小量(10-6),且流場特征長度遠(yuǎn)大于傳輸光波長[10,11],因此光束在流場中的傳輸可認(rèn)為滿足傍軸近似條件,進(jìn)而可認(rèn)為光線在流場中的傳輸軌跡近似為直線并忽略流場對振幅的衰減,只需考慮非均勻流場對相位的畸變效應(yīng)[12-14].因此,光束沿y方向傳輸通過折射率分布為n(t,x,y)、厚度為L的流場后的光程(optical path length,OPL)為

光束孔徑范圍內(nèi)任意x處光程差(optical path difference,OPD)的方差為

其中E(OPL(t,x,y))代表時(shí)間平均,σ為方差.對于理想氣體,折射率n可通過Gladstone-Dale常數(shù)KGD與密度脈動量ρ′ 關(guān)聯(lián),

相位與光程差之間滿足關(guān)系式

其中k= 2π/λ,為波數(shù);φ是相位.由(1)—(4)式可得到任意x處的波前方差為

其中β=k·KGD.進(jìn)一步,由(5)式有

統(tǒng)計(jì)理論中,協(xié)方差函數(shù)的標(biāo)準(zhǔn)定義為

由(5)和(6)式可以得到廣義的關(guān)聯(lián)方程(linking equation)

上述推導(dǎo)過程用到密度脈動均值為零的假設(shè).由(7)式可見,協(xié)方差函數(shù)的計(jì)算需要遍歷測量光路上所有點(diǎn)的密度脈動量,因此基于(8)式所示的廣義關(guān)聯(lián)方程進(jìn)行相位方差估計(jì)是一個特別耗時(shí)的過程,進(jìn)而影響其實(shí)用性.為此,在實(shí)際應(yīng)用中,通常需要進(jìn)一步對湍流做必要近似以簡化協(xié)方差函數(shù)的測量.一種通常的做法是,對湍流做均勻、各向同性的近似假設(shè).基于湍流均勻、各向同性假設(shè)理論,(8)式中的協(xié)方差函數(shù)可以擬合為指數(shù)函數(shù)形式或高斯函數(shù)形式

或

其中為密度脈動方差,lρ是湍流密度脈動特征長度.相應(yīng)地,關(guān)聯(lián)方程(8)可以分別寫成指數(shù)型

或高斯型

圖1 混合層流場圖像 (a)原始圖像;(b)噪聲和背景光預(yù)處理后圖像Fig.1.Image of mixing layer：(a) Original image;(b) image after pre-processing.

這里體現(xiàn)了均勻假設(shè)條件下光路方向的分層思想.本文采用的相干長度定義為

在此基礎(chǔ)上,Havener等[3]進(jìn)一步推導(dǎo)了包含權(quán)重函數(shù)的指數(shù)型關(guān)聯(lián)方程

或高斯型

其中,(14)和(15)式中的權(quán)重函數(shù)可以分別表示為

或

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

3.1 超聲速混合層流場

本文基于NPLS技術(shù)獲得高分辨率超聲速混合層二維密度場分布.Yi等[15]對NPLS和實(shí)驗(yàn)流場裝置已有詳細(xì)的論述,這里不再贅述.

實(shí)驗(yàn)測量的超聲速混合層流場的對流馬赫數(shù)為0.5,其中兩股來流的馬赫數(shù)分別為3.509和1.400,對應(yīng)的來流速度分別是654.7和421.1 m/s,屬于中等可壓縮流場.NPLS圖像的像素尺寸為1431 × 281,單像素分辨率為h= 0.16 mm,入射光波長λ= 1064 nm,對應(yīng)的KGD= 2.195 × 10-4m3/kg,跨幀間隔為15 μs,樣本總數(shù)為50幀.流向?yàn)閤正向,光線傳播方向?yàn)閥正向.

NPLS技術(shù)拍攝的某時(shí)刻超聲速混合層密度場圖像如圖1(a)所示.圖1(b)則是對圖1(a)進(jìn)行去噪、對照明光源固有的空間不均勻性進(jìn)行校正后得到的流場密度分布圖[16,17].對比發(fā)現(xiàn),經(jīng)過上述過程處理后,如圖1(a)圓圈內(nèi)所示的大粒子噪聲得以去除,照明光源的空間不均勻性也得到了較好的校正.

3.2 超聲速混合層流場波前相位方差結(jié)果

基于上述參數(shù),分別采用(8),(11)和(12)式計(jì)算流向波前相位方差分布,如圖2所示.

從圖2可見,密度脈動協(xié)方差函數(shù)基于指數(shù)模型和高斯模型近似后的關(guān)聯(lián)方程與廣義關(guān)聯(lián)方程計(jì)算的流向波前相位方差分布具有較好的一致性,且高斯模型近似的結(jié)果優(yōu)于指數(shù)模型結(jié)果.因此,后續(xù)的計(jì)算及誤差分析均只針對高斯模型近似下的關(guān)聯(lián)方程進(jìn)行.需要說明的是,基于(8)式的廣義關(guān)聯(lián)方程與基于(5)式的直接積分計(jì)算曲線完全重合,因此(5)式的計(jì)算曲線未在圖2中給出.

4 結(jié)果分析

如圖2所示,曲線在x/h= 800附近區(qū)域波前方差估計(jì)誤差較大.根據(jù)本文第2節(jié)的理論推導(dǎo)過程,這里主要考慮流場的權(quán)重函數(shù)對估計(jì)誤差的影響.為此,定義波前方差估計(jì)誤差表達(dá)式

基于高斯模型近似的關(guān)聯(lián)方程計(jì)算的流向波前相位方差誤差曲線如圖3所示.另外,根據(jù)統(tǒng)計(jì)學(xué)中的擬合優(yōu)度定義式

計(jì)算出圖2中(8)和(12)式計(jì)算結(jié)果對應(yīng)的兩條曲線的擬合優(yōu)度為0.8810.

圖3 高斯型關(guān)聯(lián)方程和廣義關(guān)聯(lián)方程計(jì)算波前方差的誤差Fig.3.Errors of using gaussian linking equation and generalized linking equation to calculate wave-front variance.

流體力學(xué)中常應(yīng)用脈動空間相關(guān)函數(shù)對流場的尺度特性進(jìn)行定量描述,密度脈動空間相關(guān)函數(shù)定義為[18-20]

其中,(x,y)為選取的中心點(diǎn)坐標(biāo);δ(x) ,δ(y) 為相對于中心點(diǎn)的偏移量.結(jié)合圖1和圖2,選取坐標(biāo)為(200,148),(480,140),(800,112)和(1000,100)中心點(diǎn)并計(jì)算了中心點(diǎn)附近區(qū)域的密度脈動相關(guān)函數(shù),相應(yīng)的等值線分布見圖4.密度脈動特征長度可以表征流場中的渦尺度大小.根據(jù)文獻(xiàn)[21],常將其定義為當(dāng)相關(guān)函數(shù)值為最大值的1/e時(shí)所對應(yīng)的流場尺度.基于此,計(jì)算出上述各中心點(diǎn)的密度脈動特征長度的歸一化值分別為12,44,96和112.

圖4(b)—(d)可見,在流場的中、下游中心點(diǎn)附近區(qū)域,相關(guān)函數(shù)等值線呈傾斜近似橢圓形狀,這些傾斜的橢圓形狀表示湍流大尺度結(jié)構(gòu)的存在[22-24].同時(shí)從圖4可以看到,隨著流場逐漸向下游發(fā)展,密度脈動高度相關(guān)區(qū)域在逐漸擴(kuò)大.這是由于混合層流場開始發(fā)展的初期產(chǎn)生渦量堆積,失穩(wěn)之后產(chǎn)生大尺度結(jié)構(gòu)卷起,表征了該區(qū)域流場的高度非均勻、非各向同性性.

由(17)式可知,lρ與權(quán)重函數(shù)的計(jì)算密切相關(guān).結(jié)合圖3選取了混合層流場中誤差最大的點(diǎn)x/h= 772、誤差最小的點(diǎn)x/h= 380和誤差居中的點(diǎn)x/h= 932,計(jì)算了(17)式所示權(quán)重函數(shù)的y方向分布,如圖5所示.

圖4 密度脈動相關(guān)函數(shù)分布 (a)中心點(diǎn)(200,148);(b)中心點(diǎn)(480,140);(c)中心點(diǎn)(800,112);(d)中心點(diǎn)(1000,100)Fig.4.Correlations of density fluctuations：(a) Central point (200,148);(b) central point (480,140);(c) central point (800,112);(d) central point (1000,100).

圖5 部分流向點(diǎn)處的權(quán)重函數(shù)分布Fig.5.Distribution of weighting functions at some streamwise locations.

由圖5可見,在x/h= 772處,權(quán)重函數(shù)偏離1的區(qū)域最多,因此導(dǎo)致由于忽略權(quán)重函數(shù)的(12)式計(jì)算的相位方差誤差也最大;相反,x/h= 380位置處權(quán)重函數(shù)偏離1的區(qū)域最少,因此在該位置處(12)式計(jì)算的波前方差誤差值最小.進(jìn)一步,給出了(12)和(15)式中不包含權(quán)重函數(shù)的積分核函數(shù)和包含權(quán)重函數(shù)的積分核函數(shù)分布,如圖6(a)和圖6(b)所示,以及二者的差值分布如圖6(c)所示.由圖6(c)可知,權(quán)重函數(shù)影響集中在流場自由邊界一側(cè)、且位于x/h= 772處附近.由此可見,權(quán)重函數(shù)對于關(guān)聯(lián)方程在超聲速流場密度脈動高度相關(guān)區(qū)域中應(yīng)用的必要性.

對權(quán)重函數(shù)(17)式中的Erf(·)函數(shù)進(jìn)行進(jìn)一步分析發(fā)現(xiàn),只有當(dāng)L?lρ,即L/lρ較大時(shí),W(y)值越接近1;反之,W(y)值偏離1越遠(yuǎn).由于流場的厚度L是一不變量,因此密度脈動特征長度lρ的大小決定著權(quán)重函數(shù)W(y)的值.由于流場失穩(wěn)導(dǎo)致大渦結(jié)構(gòu)的存在,在流場的中游附近區(qū)域L/lρ較小,所以權(quán)重函數(shù)的值會逐漸偏離1較遠(yuǎn),此時(shí)忽略權(quán)重函數(shù)的影響自然會導(dǎo)致波前方差計(jì)算誤差的增加.

最后,基于上述分析,加入權(quán)重函數(shù)重新計(jì)算波前方差,繪出基于(15)式計(jì)算波前相位方差的曲線,結(jié)果如圖7所示.可以看到,高斯形式的關(guān)聯(lián)方程(12),在加入權(quán)重函數(shù)(17)式后,計(jì)算得到的波前相位方差曲線的擬合效果要明顯變好.(8)和(15)式計(jì)算結(jié)果對應(yīng)的兩條曲線的擬合優(yōu)度為0.9127.

圖6 高斯型關(guān)聯(lián)方程加入權(quán)重函數(shù)前后的積分核分布(a)未加入權(quán)重函數(shù);(b)加入權(quán)重函數(shù);(c)積分核分布差Fig.6.Integral kernel distribution calculated by Gaussian linking equation before and after adding weighting function：(a) Before adding the weighting function;(b) after adding the weighting function;(c) the integral kernel distribution differences.

5 結(jié) 論

圖7 高斯型關(guān)聯(lián)方程加入權(quán)重函數(shù)前后計(jì)算波前方差對比Fig.7.Wave-front variance calculated by Gaussian linking equation before and after adding weighting function.

本文分析了關(guān)聯(lián)方程在超聲速混合層流場的應(yīng)用中,密度相關(guān)函數(shù)采用高斯模型簡化處理所產(chǎn)生的權(quán)重函數(shù)對波前方差估計(jì)精度的影響.研究結(jié)果表明：由于混合層大尺度結(jié)構(gòu)的卷起,導(dǎo)致密度脈動高度相關(guān);在密度脈動高度相關(guān)區(qū)域附近,密度脈動特征長度值與權(quán)重函數(shù)偏離1的程度密切相關(guān);當(dāng)密度脈動特征長度較大時(shí),權(quán)重函數(shù)偏離1的程度較嚴(yán)重,反之則相反.研究結(jié)果表明,對于在超聲速混合層流場密度脈動高度相關(guān)區(qū)域應(yīng)用關(guān)聯(lián)方程計(jì)算波前相位方差考慮權(quán)重函數(shù)的必要性;加入權(quán)重函數(shù)的關(guān)聯(lián)方程計(jì)算波前方差的精度優(yōu)于未加權(quán)重函數(shù)的關(guān)聯(lián)方程計(jì)算結(jié)果.

感謝國防科技大學(xué)空天科學(xué)學(xué)院易仕和教授在實(shí)驗(yàn)設(shè)備和流場產(chǎn)生裝置等方面提供的支持與幫助.