具有湍流效應(yīng)的不規(guī)則微下?lián)舯┝鳂?gòu)建方法研究

劉時杰, 徐浩軍, 薛源, 張久星

(空軍工程大學(xué) 航空航天工程學(xué)院, 陜西 西安 710038)

0 引言

微下?lián)舯┝鲗儆谙聸_氣流的一種,是以垂直風(fēng)切變?yōu)橹饕卣鞯膹?fù)雜風(fēng)切變區(qū)域[1-2]。研究微下?lián)舯┝鲗τ陲w機起飛、著陸時的飛行品質(zhì)分析與飛行安全保障具有重要的意義。微下?lián)舯┝髂Ｐ偷难芯磕壳坝卸喾N方法,例如利用基于風(fēng)洞實驗與CFD數(shù)值計算的方式研究其速度場分布[3-4],以及基于融合重疊的方法[5]和基于渦環(huán)原理的方法[6-8]。其中Woodfield和Woods[6]最早提出的渦環(huán)原理模型是目前應(yīng)用最廣泛的一種微下?lián)舯┝髂Ｐ汀?/p>

目前的方法所構(gòu)建的微下?lián)舯┝髂Ｐ投酁檠刂行妮S對稱的規(guī)則形狀,但實際的微下?lián)舯┝黠L(fēng)場非常復(fù)雜,其下沉氣流并不一定垂直于地面,其外流范圍也不是軸對稱的,且一個區(qū)域內(nèi)可能伴隨有多個風(fēng)暴核,甚至還有局部上升氣流。并且在其風(fēng)速的分布范圍內(nèi)亦可能伴隨有其他大氣流動狀況,如不規(guī)則突風(fēng)與飛機尾流產(chǎn)生的湍流。故本文擬在構(gòu)建不規(guī)則微下?lián)舯┝髂Ｐ偷幕A(chǔ)上考慮突風(fēng)與湍流的影響,為起飛、著陸與低空飛行時的安全性仿真提供更加準(zhǔn)確的風(fēng)場輸入條件。

1 不規(guī)則微下?lián)舯┝黠L(fēng)場建模

1.1 單個渦環(huán)誘導(dǎo)下的風(fēng)場構(gòu)建

首先利用渦環(huán)原理構(gòu)造微下?lián)舯┝黠L(fēng)場模型。依據(jù)美式坐標(biāo)系,建立一個如圖1所示的渦環(huán)模型。在地面處,風(fēng)速垂直方向的分量為0。通過在地面上方布置一個強度為Γ的主渦環(huán),同時在對稱下方布置一個強度為-Γ的鏡像渦環(huán),即可滿足地面垂直速度為0的邊界條件。圖1中,主渦環(huán)中心在P(xP,yP,zP)處,鏡像渦環(huán)中心在(xP,yP,-zP)處。,R為渦環(huán)半徑,A為空間中某點,坐標(biāo)為(xA,yA,zA)。

圖1 單渦環(huán)示意圖Fig.1 Single vortex ring

首先討論單渦環(huán)影響下的微下?lián)舯┝鲌鰳?gòu)建方法,當(dāng)渦環(huán)面與地面平行時,渦環(huán)的流線方程為:

ψ=-(Γ/2π)(r1+r2)[F1(k)-F2(k)]

(1)

式中,Γ為渦環(huán)強度;r1為A點距渦環(huán)距離最近的距離;r2為A點距渦環(huán)最遠點的距離;k=|(r2-r1)/(r2+r1)|;F1(k)與F2(k)為圓環(huán)的積分函數(shù)。若0≤k2≤1,則:

(2)

將式(2)代入式(1),獲得主渦環(huán)的流線方程為:

(3)

鏡像渦環(huán)的流線表達式為:

(4)

最終將主渦環(huán)流線表達式與鏡像渦環(huán)流線表達式線性疊加,得到空間點A的流線方程為:

(5)

而后得到點A的誘導(dǎo)速度為:

(6)

(7)

wz=(-1/rA)?ψ/?rA

(8)

式中,rA為點A距渦環(huán)對稱軸的距離。

1.2 中心軸處的誘導(dǎo)速度計算

在渦環(huán)的中心軸處rA=0,用式(6)～式(8)計算誘導(dǎo)速度會產(chǎn)生奇點。因此,通過引入渦環(huán)的位函數(shù),利用式(10)計算垂直方向的速度wz;由于渦環(huán)沿中心軸的對稱性,水平速度wx和wy為0。

wz(z)=(Γ/2R)[1/(1+(z/R)2)3/2]

(9)

式中,z為中心軸上各點離渦環(huán)中心的距離。

1.3 渦核內(nèi)部的誘導(dǎo)速度計算

根據(jù)式(6)～式(8),在渦環(huán)的中心渦絲處,流速趨于無窮大,這不符合微下?lián)舯┝黠L(fēng)場的實際情況。由于流體粘性的影響,實際渦環(huán)有一個渦核的存在,且渦核內(nèi)速率逐漸減小,至渦核的中心渦絲處減為0。本文運用復(fù)合渦原理,將渦核看作半徑為r的圓柱,渦核內(nèi)部的渦量均勻分布,保證渦絲處的流速為0,而渦核外部仍然服從流線方程。從渦核中心到渦核半徑處,流速呈線性分布,如圖2所示。

圖2 渦核內(nèi)部誘導(dǎo)速度示意圖Fig.2 Velocity induced inside the vortex core

設(shè)渦核半徑為r,首先應(yīng)該判斷點A(xA,yA,zA)是否在渦核內(nèi),若滿足下式,則點A在渦核內(nèi)。

(rA-R)2+(zP-zA)2≤r2

(10)

根據(jù)點A的坐標(biāo)、渦環(huán)中心的坐標(biāo)求出與點A、渦環(huán)中心點共面的渦核絲坐標(biāo)點O(xO,yO,zO)。再求出與點A、點O共線,與點A、渦核中心點共面的渦核外壁坐標(biāo)N,如圖2所示。根據(jù)式(6)～式(8)求出坐標(biāo)N的誘導(dǎo)速度,再根據(jù)點A、點N距離點O的距離按比例求出A點的誘導(dǎo)速度。

1.4 微下?lián)舯┝黠L(fēng)場的不規(guī)則變換

根據(jù)美式坐標(biāo)系確定旋轉(zhuǎn)軸心為渦環(huán)的中心點處(見圖1),引入俯仰變換矩陣Lθ、滾轉(zhuǎn)變換矩陣Lφ、偏航變換矩陣Lψ,對坐標(biāo)系及速度進行旋轉(zhuǎn)變換,從而構(gòu)建不規(guī)則速度場。對于主渦環(huán)坐標(biāo)系的變換公式為式(14),對于速度矢量的變換公式為式(15)。鏡像渦環(huán)的俯仰變換角度和滾轉(zhuǎn)變換角度與主渦環(huán)正負相反,偏航變換角度與主渦環(huán)相同。

(11)

(12)

(13)

(14)

(15)

經(jīng)過變換后的坐標(biāo)排列亦是不規(guī)則的,因為飛行實時仿真需要根據(jù)微下?lián)舯┝髂Ｐ蛣討B(tài)插值每個坐標(biāo)點的速度矢量,故這種混亂不規(guī)則的坐標(biāo)排列在飛行實時仿真中是無法使用的。故需將變換后坐標(biāo)點對應(yīng)的速度矢量擬合到標(biāo)準(zhǔn)的升序排列三維坐標(biāo)點上,如式(16)。文中采用的擬合方法為線性插值法,由于坐標(biāo)點較多,達到了100萬以上,故擬合計算的過程相對于其他步驟來說耗時較長。

(16)

將標(biāo)準(zhǔn)坐標(biāo)系內(nèi)的主渦環(huán)與鏡像渦環(huán)的誘導(dǎo)速度場相疊加,得到渦環(huán)1總的誘導(dǎo)速度場[wx1,wy1,wz1]。

1.5 多個渦環(huán)誘導(dǎo)風(fēng)速的疊加

由于單個渦環(huán)誘導(dǎo)的微下?lián)舯┝鲌龃怪憋L(fēng)速相對較小,為了更加逼真,通常使用多個渦環(huán)誘導(dǎo)微下?lián)舯┝鲌?。多個傾斜渦環(huán)疊加而成的微下?lián)舯┝魉俣葓鰹?

(17)

式中,[wxn,wyn,wzn]為第n個渦環(huán)的誘導(dǎo)速度場。

2 低空大氣湍流風(fēng)場結(jié)構(gòu)的構(gòu)建

湍流風(fēng)場采用Dryden模型[9]得到,其單側(cè)一維頻譜函數(shù)為:

(18)

式中,Lwx,Lwy,Lwz為湍流尺度;σwx,σwy,σwz為湍流強度。文中Lwx=2Lwy=2Lwz=380 m,σwx=σwy=σwz=3.6 m/s。

在飛行實時仿真中,關(guān)鍵問題在于如何產(chǎn)生高保真度的大氣湍流。產(chǎn)生的湍流速度不僅要具有隨機性,還必須符合其頻譜特性。本文利用白噪聲通過一定形式的濾波器來產(chǎn)生湍流速度。白噪聲X(t)是一種具有特殊性質(zhì)的隨機過程,它的頻譜函數(shù)等于常數(shù),其相關(guān)函數(shù)為脈沖函數(shù),如下式:

ΦXX(ω)=N0,RXX(τ)=N0δ(τ)

讓白噪聲X(t)通過一個傳遞函數(shù)為G(s)的環(huán)節(jié),設(shè)其輸出為Y(t)。輸出Y(t)的頻譜函數(shù)為:

ΦYY(ω) =G*(iω)ΦXX(ω)G(iω)

=N0G*(iω)G(iω)

由上式可見,若要求輸出Y(t)具有特定的頻譜ΦYY(ω),應(yīng)將ΦYY(ω)作因式分解為G*(iω)和G(iω)的乘積,從而找到濾波器應(yīng)具有的傳遞函數(shù)G(s)。

以Dryden模型(式(18))的第一式為例進行因式分解。將其按空間頻率和時間頻率的轉(zhuǎn)換關(guān)系Φ(ω)=(1/v*)Φ(ω/v*)轉(zhuǎn)化為:

對其進行因式分解:

Φwxwx(ω)=

所以其對應(yīng)的傳遞函數(shù)為:

同理可得Gwy(s)和Gwz(s)。

根據(jù)以上論述,三個方向的大氣湍流速度的產(chǎn)生如圖3所示。

圖3 湍流速度的產(chǎn)生Fig.3 Generation of turbulence velocity

將微下?lián)舯┝魉俣葓雠c尾流引起的湍流速度場在坐標(biāo)系內(nèi)進行線性疊加,從而得到受湍流效應(yīng)影響的微下?lián)舯┝魉俣葓觥?/p>

3 計算結(jié)果及分析

選擇三個渦環(huán)來構(gòu)成不規(guī)則微下?lián)舯┝黠L(fēng)場,坐標(biāo)范圍為x:0～6 km,y:0～6 km,z:0～0.8 km。第1個渦環(huán)的中心點坐標(biāo)為(3.0,3.0,0.6)km,R=0.9 km,r=0.455 km,Γ=18 000 m2/s,俯仰偏角為15°,滾轉(zhuǎn)偏角為5°;第2個渦環(huán)的中心點坐標(biāo)為(2.5,2.5,0.5) km,R=0.7 km,r=0.3 km,Γ=-10 000 m2/s,俯仰偏角為10°,滾轉(zhuǎn)偏角為-5°;第3個渦環(huán)的中心點坐標(biāo)為(3.0,4.0,0.6) km,R=0.8 km,r=0.4 km,Γ=11 000 m2/s,俯仰偏角為25°,滾轉(zhuǎn)偏角為15°。所得結(jié)果如圖4～圖9所示。

從圖4～圖8中可以看出,相比于規(guī)則風(fēng)場,文中對不規(guī)則微下?lián)舯┝鞯臉?gòu)建是成功有效的,構(gòu)建的風(fēng)場不再沿軸線對稱,比較符合實際情況;并且由于加入了湍流的影響,速度場存在明顯的不規(guī)則湍動,較好地反映了低空大氣湍流對風(fēng)場隨機性的影響。

圖4 在剖面y=3 km上的風(fēng)速矢量圖Fig.4 Wind speed vector on the profile y=3 km

圖5 在剖面z=0.4 km上的風(fēng)速矢量圖Fig.5 Wind speed vector on the profile z=0.4 km

圖6 在剖面y=3 km上的wz云圖Fig.6 Cloud map of wz on the profile y=3 km

圖7 在剖面z=0.4 km上的wx云圖Fig.7 Cloud map of wx on the profile z=0.4 km

圖8 在剖面x=3.5 km上的wz云圖Fig.8 Cloud map of wz on the profile x=3.5 km

圖9 三維空間內(nèi)wz=-5 m/s的拓撲結(jié)構(gòu)圖Fig.9 Topology structure in the 3D space with wz=-5 m/s

4 結(jié)束語

本文綜合考慮了微下?lián)舯┝鞯碾S機性與不確定性,通過多渦環(huán)建模法構(gòu)建了具有旋轉(zhuǎn)傾斜效應(yīng)的不規(guī)則微下?lián)舯┝魉俣葓?并將其與低空大氣湍流相互疊加以反映風(fēng)場在湍流作用下的隨機性。所得風(fēng)場能較好地描述多因素耦合復(fù)雜情況下的大氣流動狀況,相比于以往對微下?lián)舯┝魉俣葓龅难芯扛由钊?、系統(tǒng),也更接近復(fù)雜風(fēng)場的實際情況。對于起飛、著陸階段飛行安全性仿真分析、復(fù)雜風(fēng)場中的風(fēng)險規(guī)避以及控制技術(shù)研究具有積極的意義,亦可以進一步補充完善大氣流動數(shù)據(jù)庫。

參考文獻：

[1] Fujita T T.Tornadoes and downbursts in the context of generalized planetary scales [J].Journal of Atmospheric Science,1981,38(8):1511-1534.

[2] 朱上翔.微下?lián)舯┝鞯牧黧w動力學(xué)模型[J].飛行力學(xué),1984,2(2):59-72.

[3] Qu W,Ji B.Numerical study on formation and diffusion wind fields for thunderstorm microburst[C]//Proc.of the International Conference on Mechanic Automation and Control Engineering.U.S.:IEEE,2010:1389-1392.

[4] 瞿偉廉,王錦文.下?lián)舯┝黠L(fēng)荷載的數(shù)值模擬[J].武漢理工大學(xué)學(xué)報,2008,30(2):70-74.

[5] Nguyen H,Manuel L,Veers P.Simulation of inflow velocity fields and wind turbine loads during thunderstorm downbursts[R].AIAA-2010-2651,2010.

[6] Woodfield A A,Woods J F.Worldwide experice of wind shear during 1981-1982[R].ADP002705,1983.

[7] DoD.MIL-HDBK-1797 Military Handbook[S].U.S.:DoD,1997.

[8] DoD.MIL-F-8785C Military Specification[S].U.S.:DoD,1980.

[9] Yeager J.Implementation and testing of turbulence models for the F18-HARV simulation[R].NASA CR-1998-206937,1998.