基于蒙特卡羅法的四維大氣紊亂仿真

基于蒙特卡羅法的四維大氣紊亂仿真

0紊流場的生成空氣擾動是飛機在飛行過程中經(jīng)常發(fā)現(xiàn)的空氣擾動現(xiàn)象。這是重疊定常風(fēng)中風(fēng)速變量的連續(xù)隨機變量。飛機在紊流場中飛行時會發(fā)生顛簸,直接影響乘座舒適性和飛行安全,同時也是飛機結(jié)構(gòu)疲勞損壞的一大誘因。在飛行模擬器飛行仿真系統(tǒng)中,需要建立紊流場的數(shù)學(xué)模型,同時在實時仿真中高逼真度地模擬紊流對飛機六自由度運動的影響。由于紊流形成的機理復(fù)雜,為適應(yīng)實時仿真的需要,應(yīng)引入若干的簡化假設(shè),建立工程化的數(shù)學(xué)模型。國外從20世紀50年代就對大氣紊流現(xiàn)象進行理論和試驗研究,內(nèi)容包括大氣紊流的建模方法,飛機對紊流的響應(yīng)分析和以改善乘座品質(zhì)為主要內(nèi)容的紊流顛簸問題等,并提出了VonKarman模型和Dryden模型,前者是基于實測大氣紊流而得出的,后者是對前者的近似。針對紊流場生成國外提出了蒙特卡羅法和統(tǒng)計離散風(fēng)兩種建模方法,其中,蒙特卡羅法是用高斯隨機白噪聲通過成型濾波器,獲得符合特定模型的統(tǒng)計特性,進而變換為時域有限差分方程形式的遞推算法,成型濾波器函數(shù)是通過對模型頻譜函數(shù)進行因式分解獲得的,但這種方法應(yīng)用在二維或三維紊流場中會造成各向不同性,此外,由于VonKarman模型頻譜函數(shù)不具有有理因式分解結(jié)果,難以推導(dǎo)為時域遞推算法。國內(nèi)自20世紀80年代起,先后有學(xué)者進行了符合Dryden模型的一維、二維紊流場的生成研究,其中,二維紊流場同樣存在各向不同性問題。近年來,有學(xué)者基于紊流場相關(guān)函數(shù)矩陣生成了Dryden模型的三維紊流場,其具有良好的均勻性和各向同性,但算法復(fù)雜,不適于實時仿真中在線生成。在以往的飛行仿真研究中,常常將飛機簡化為質(zhì)點,僅模擬飛機縱向的紊流顛簸,因而僅需要生成一維紊流場,對紊流在大范圍空間內(nèi)的強度和尺度變化,以及紊流場的連續(xù)擴展問題研究得不夠。對大型飛機而言,采用一維紊流場則忽略了沿翼展和機身方向的紊流風(fēng)梯度變化,無法準確模擬飛機對紊流擾動的六自由度動力學(xué)響應(yīng),因此,需要建立三維空間大氣紊流場。Dryden模型與VonKarman模型相比,無法模擬出飛機對紊流的高頻響應(yīng),因而對于高頻結(jié)構(gòu)彈性振動明顯的大型飛機而言,使用Dryden模型會造成模擬逼真度不夠。本文研究了可用于大型飛機實時飛行仿真的三維空間大氣紊流場的生成及擴展方法,針對以往紊流場生成方法上存在的問題,通過在頻域空間直接生成符合VonKarman模型的三維紊流場,避免了各向不同性問題,再基于三維傅立葉逆變換轉(zhuǎn)化為時域空間紊流場;根據(jù)傅立葉變換的復(fù)共軛對稱特性設(shè)計算法,保證獲得實的紊流場,同時降低計算開銷;通過預(yù)先生成無因次的紊流場,飛行仿真過程中再進行有因次化,可使紊流場強度和尺度隨飛行高度動態(tài)變化。本文研究的基于對稱擴展的紊流場擴展方法,可保證生成的局部空間紊流場有效擴展為大范圍連續(xù)紊流場。1dryhen模型本文研究生成的大氣紊流場是基于各向同性、高斯分布和凍結(jié)場假設(shè)的,若在頻域空間直接生成紊流場,則可避免時域差分方程的推導(dǎo),因此,本文在頻域內(nèi)運用蒙特卡羅法進行隨機抽樣,將三路高斯白噪聲序列ni(r)(i=1,2,3,r為空間矢量)直接在頻域中與三維成型濾波函數(shù)Hi(f1,f2,f3)相乘,得結(jié)果wi(r),再通過三維傅立葉逆變換(3D-IDFT),獲得時域中各向同性三維紊流場(u,v,w)T,具體過程見圖1。該算法由于無需推導(dǎo)差分方程,不要求獲得有理的頻譜因式分解形式,使得Dryden模型和VonKarman模型在這種方法下進行計算時計算量是一樣的。為進行數(shù)值仿真的需要,須構(gòu)建離散化三維紊流場(圖2)。圖2中,紊流場在頻域空間劃分為M1×M2×M3個網(wǎng)格點。由于紊流場須通過頻域到時域的變換生成,根據(jù)傅立葉變換的基本原理,其時域空間采樣頻率與頻域采樣頻率滿足如下關(guān)系Δfi=fsiΜii=1,2,3(1)Δri=1fsi(2)Δfi=fsiMii=1,2,3(1)Δri=1fsi(2)式中:fsi為xi方向的時域空間采樣頻率;Mi為xi方向空間網(wǎng)格點數(shù);Δfi為xi方向頻域空間間隔;Δri為xi方向時域空間間隔。鑒于紊流尺度和強度在真實大氣環(huán)境中隨飛行高度而變化,本文首先將空間距離無因次化,即■SymbolYCpfsi=fsiL,■SymbolYCpri=riL(3)ri=ΔriΜi式中:L為紊流場尺度;■SymbolYCpfsi為無因次空間采樣頻率;■SymbolYCpri為無因次空間距離。2基于時頻變換方法的混亂流場的生成2.1無因次及功率譜函數(shù)VonKarman模型自功率譜函數(shù)為?ii(f)=440π39σ2a4L5(f2-f2i)[1+(2πaLf)2]17/6(4)f2=f21+f22+f23式中:?ii(f)(i=1,2,3)分別為縱向和2個橫向的功率譜;參數(shù)a為1.339;σ為紊流強度;fi為三軸空間頻率。根據(jù)文獻可推導(dǎo)出無因次譜與有因次譜之間的關(guān)系為σ2=∫∞-∞∫∞-∞∫∞-∞■SymbolYCp?ii(■SymbolYCpf1,■SymbolYCpf2,■SymbolYCpf3)d■SymbolYCpf1d■SymbolYCpf2d■SymbolYCpf3=∫∞-∞∫∞-∞∫∞-∞?ii(■SymbolYCpf1L,■SymbolYCpf2L,■SymbolYCpf3L)d■SymbolYCpf1d■SymbolYCpf2d■SymbolYCpf3L3(5)■SymbolYCp?ii(■SymbolYCpf)=?ii(f)/L3(6)綜合式(3)、(4)、(6),最終得到無因次VonKarman模型功率譜函數(shù)為■SymbolYCp?ii(■SymbolYCpf)=440π39σ2a4(■SymbolYCpf2-■SymbolYCpf2i)[1+(2πa■SymbolYCpf)2]17/6(7)式中:■SymbolYCp?ii(■SymbolYCpf)、■SymbolYCpfi、■SymbolYCpf為相應(yīng)無因次變量。2.2維濾波器函數(shù)本文采用高斯白噪聲序列沖擊數(shù)字濾波器,白噪聲的采樣頻率與相應(yīng)的三軸時域空間采樣頻率一致。一般數(shù)值生成的白噪聲均為帶限白噪聲,其功率譜滿足如下關(guān)系∫fs/2-fs/2?Ν(f)dfs=σ2Ν(8)本文用白噪聲信號激勵線性濾波器H(f),其輸出頻譜為?w(f)=|H(f)|2?N(f)=H(f)H*(f)?N(f)(9)結(jié)合式(8)、(9),得到輸出濾波器傳遞函數(shù)為Η(f)=√fs?w(f)σΝejα(f)(10)式中:相角α(f)可不予考慮,不妨選擇為0;fs為空間采樣頻率;?N(f)為帶限白噪聲信號功率譜函數(shù);σN為帶限白噪聲信號標準差;H*(f)為H(f)的共軛。將以上關(guān)系推廣到三維情況,則有帶限白噪聲信號滿足∫fs1/2-fs1/2∫fs2/2-fs2/2∫fs3/2-fs3/2?N(f)dfs1dfs2dfs3=σ2N(11)經(jīng)類似的推導(dǎo),得出相應(yīng)的三維濾波器函數(shù)為Η(f)=√fs1fs2fs3?w(f)σΝ(12)式(12)為三維濾波器函數(shù)的一般形式。對于VonKarman模型,可令?w(f)=■SymbolYCp?ii(■SymbolYCpf),從而獲得相應(yīng)的傳遞函數(shù)H(f)。2.3傅立葉變換時域及頻域數(shù)額序列三路白噪聲信號在頻域中完成對各自功率譜函數(shù)的相乘后,須使用離散傅立葉逆變換進行頻域到時域的轉(zhuǎn)換。三維離散傅立葉逆變換表達式為x(n1,n2,n3)=Μ1-1∑k1=0Μ2-1∑k2=0Μ3-1∑k3=0X(k1,k2,k3)?ej2π(n1k1/Μ1+n2k2/Μ2+n3k3/Μ3)(13)式中:x(n1,n2,n3)、X(k1,k2,k3)分別為三維傅立葉變換時域及頻域數(shù)值序列,ni=0,1,…,Mi-1。傅立葉變換是針對一組序列值進行的,序列長度應(yīng)是2的整數(shù)倍?？紤]到必須生成實數(shù)形式的時域空間紊流場,根據(jù)傅立葉變換的基本性質(zhì),須使頻域內(nèi)的紊流數(shù)值序列滿足復(fù)共軛對稱關(guān)系。對于一維情況,要求頻域序列中心點為實數(shù),其他數(shù)值相對中心點呈共軛對稱;對于二維情況,兩正交軸上的序列要求與各自中心點對稱,網(wǎng)格內(nèi)的點則呈中心對稱;對于三維情況,在三軸和3個平面滿足各自對稱性的前提下,網(wǎng)格內(nèi)的點呈中心對稱,從而使頻域內(nèi)的數(shù)值序列滿足表1所示的關(guān)系。2.4大氣紊流場實時仿真通過圖1所示的過程可生成無因次的空間紊流場,該紊流場可預(yù)先生成,存放在擾動風(fēng)場數(shù)據(jù)庫中,在實時飛行仿真中,根據(jù)飛機飛行高度獲得紊流場的尺度,按式(3)進行無因次到有因次的轉(zhuǎn)換,最終乘以紊流強度獲得三維空間大氣紊流場。可以看出,有因次化后,不同高度的紊流場的時域網(wǎng)格的間隔將發(fā)生變化。若要生成時域等間隔的空間紊流場,采用不同的采樣頻率即可。3該算法與混沌流場擴展方法的研究3.1縱向流場試驗本文設(shè)計了32×32×32的網(wǎng)格空間,設(shè)置時域空間采樣頻率fs1、fs2、fs3均為100rad·m-1,進行算法仿真。以縱向紊流場試驗為例,紊流場的某一剖面的縱向風(fēng)速值見圖3,其他剖面的風(fēng)速值與圖3類似,呈隨機性地劇烈變化。計算結(jié)果顯示,產(chǎn)生的紊流場數(shù)值虛部全部近似為0,從而證實了根據(jù)復(fù)共軛對稱關(guān)系進行傅立葉逆變換的正確性。傳統(tǒng)的時域數(shù)值仿真算法是在仿真步長點處實時計算風(fēng)速矢量,而本文的方法是在仿真前生成紊流場并存儲在風(fēng)場數(shù)據(jù)庫中,仿真中只需查詢獲取數(shù)值即可,從而有效減少實時計算的開銷。3.2特定函數(shù)為檢驗所產(chǎn)生的三維紊流場的有效性,本文采用了在時域中檢驗紊流場數(shù)值相關(guān)性的方法,重點檢查了紊流場的一維自相關(guān)性和二維互相關(guān)性。對各向同性的三維紊流場,有R11(ξ1,0,0)=R22(0,ξ2,0)=R33(0,0,ξ3)=σ2f(ξ)(14)R11(0,ξ2,0)=R11(0,0,ξ3)=R22(ξ1,0,0)=R22(0,0,ξ3)=R33(ξ1,0,0)=R33(0,ξ2,0)=σ2g(ξ)(15)式中:Rij(·)為相關(guān)函數(shù);f(ξ)、g(ξ)分別為縱向及橫向相關(guān)函數(shù);ξi、ξ為空間距離變量。本文對生成的紊流場進行大量相關(guān)性檢驗,獲得某一剖面x1方向的相關(guān)性檢驗結(jié)果(圖4),圖中以ζ=ξ/aL作為橫坐標自變量。通過檢驗證明,相關(guān)性變化趨勢合理,產(chǎn)生的紊流場滿足相關(guān)性要求。對于互相關(guān)(i≠j)的情況,相關(guān)函數(shù)為Rij(ξ)=σ2{[f(ξ)-g(ξ)]ξiξjξ2}(16)當(dāng)ξi=ξj時,式(16)取到最大值,因此,本文重點考察了二維紊流場的對角線方向(如圖1所示的OA、OB方向等)的互相關(guān)特性。通過大量檢驗證明,計算結(jié)果與理論值接近(圖5),由此可見,其二維互相關(guān)特性亦符合VonKarman模型。3.3特定損害空間在飛行仿真過程中,若要模擬飛機連續(xù)穿越紊流場的情況,必然存在紊流場的空間擴展問題,對此,有學(xué)者提出一種簡易的方法,將飛機在航跡坐標系下的位置(xE,yE,zE)進行重新定位,其表達式為x′E=xEmodr1(17)y′E=yEmodr2(18)z′E=zEmodr3(19)式中:mod為取余符號。通過此方法可使飛機在任意空間運動時,保證紊流場不間斷,實現(xiàn)首尾相接。本文研究發(fā)現(xiàn),使用該方法連接后的空間紊流場的相關(guān)性檢驗結(jié)果存在一定問題。如對于32×32×32紊流場依次相接,拼接成96×96×96的紊流場后,其相關(guān)性檢驗結(jié)果見圖6。可以看出,在紊流場相接處的相關(guān)性檢驗結(jié)果與理論值有較大出入。為改善這種情況,本文提出了基于對稱原理的紊流場拼接方法,對二維紊流場(圖7(a)),設(shè)第Ⅰ象限為已生成的M1×M2紊流場,其內(nèi)部網(wǎng)格點上風(fēng)速矢量為V(i,j)。若需擴展這片紊流場,可利用對稱關(guān)系,設(shè)定第Ⅱ象限內(nèi)部的風(fēng)速矢量為V(i,M2-j),第Ⅲ象限內(nèi)部的風(fēng)速矢量為V(M1-i,M2-j),第Ⅳ象限內(nèi)部的風(fēng)速矢量為V(M1-i,j)。周圍的紊流場可按此方法循環(huán)擴展。對三維紊流場而言(圖7(b)),以象限Ⅰ為基礎(chǔ)進行擴展,共有8種情形,設(shè)第Ⅰ象限處為已生成的M1×M2×M3的紊流場,其內(nèi)部網(wǎng)格點上的風(fēng)速矢量為V(i,j,k),按照類似的方法可得Ⅱ～Ⅷ象限的風(fēng)速矢量分別為V(M1-i,j,k)、V(M1-i,M2-j,k)、V(i,M1-j,k)、V(i,j,M3-k)、V(M1-i,j,M3-k)、V(M1-i,M2-j,M3-k)、V(i,M2-j,M3-k)。本文同樣設(shè)計了32×32×32的紊流場,以該紊流場為中心,按照這種方法擴展成27個空間區(qū)域,成為96×96×96的紊流場。其相關(guān)性檢驗結(jié)果見圖8,可以看出,對稱拼接的紊流場更好地符合相關(guān)性理論。紊流場可經(jīng)簡單的對稱拼接方法擴展到無限大,從而有效解決