Von Karman模型三維大氣紊流仿真理論與方法

高 靜 洪冠新 梁灶清

(北京航空航天大學(xué) 航空科學(xué)與工程學(xué)院，北京 100191)

Von Karman模型三維大氣紊流仿真理論與方法

高 靜 洪冠新 梁灶清

(北京航空航天大學(xué) 航空科學(xué)與工程學(xué)院，北京 100191)

基于三維紊流場相關(guān)函數(shù)矩陣，用蒙特卡羅法仿真生成Von Karman模型三維空間大氣紊流場，實(shí)現(xiàn)了對復(fù)雜的Von Karman模型的精確推導(dǎo)和應(yīng)用，仿真流場的相關(guān)特性在理論上可以任意逼近Von Karman理論模型.通過提高白噪聲序列的白化程度及優(yōu)化蒙特卡羅算法的方法提高了仿真精度和效率.數(shù)值仿真結(jié)果證明:所生成的空間大氣紊流場符合模型特性，具有較好的統(tǒng)計(jì)特性和仿真精度，滿足實(shí)時(shí)飛行仿真對紊流場數(shù)據(jù)的要求.

大氣紊流;蒙特卡羅法;Von Karman模型;數(shù)值仿真;飛行仿真

大氣紊流對在大氣層中飛行的各種飛行器，包括飛機(jī)、導(dǎo)彈、運(yùn)載火箭、空天飛機(jī)以及微小型飛行器、氣球、飛艇等的設(shè)計(jì)、使用安全等都有較大影響.大氣紊流模型中，與 Dryden模型相比，Von Karman模型的頻譜函數(shù)在高頻段的斜率不同，頻譜函數(shù)更為復(fù)雜，更符合大氣紊流的實(shí)際情況［1］.在研究涉及飛機(jī)結(jié)構(gòu)振動的飛行品質(zhì)、飛機(jī)的結(jié)構(gòu)疲勞等問題時(shí)，由于飛機(jī)結(jié)構(gòu)模態(tài)頻率通常恰好處在高頻范圍，因而高頻范圍內(nèi)紊流可能激發(fā)飛機(jī)的結(jié)構(gòu)振動，只要可行，最好使用Von Karman大氣紊流模型.

對大氣紊流進(jìn)行建模與仿真研究中，一維大氣紊流序列的仿真一般采用成型濾波器法，此方法要求模型的頻譜函數(shù)能夠進(jìn)行共軛分解，因此無法基于Von Karman模型進(jìn)行時(shí)域仿真，工程應(yīng)用較少.但在分析空中加油、編隊(duì)飛行等多機(jī)飛行仿真及飛行品質(zhì)等研究中，必須使用模型的二維或三維紊流場，且最好采用基于Von Karman模型.文獻(xiàn)［2－3］中生成二維紊流場的方法，在理論上有缺陷，各向同性不滿足.文獻(xiàn)［4］對 Von Karman模型進(jìn)行有理化逼近，簡化了 Von Karman模型，生成的并不是真正的空間三維大氣紊流場.文獻(xiàn)［5］中提出的空間三維大氣紊流場的理論模型和數(shù)值方法，擺脫了以上仿真模型的局限，該方法不用對模型進(jìn)行簡化，生成的三維大氣紊流場具有良好的均勻性和各向同性.但是文獻(xiàn)［5］只是基于 Dryden模型，沒有實(shí)現(xiàn)基于 Von Karman模型的三維紊流場仿真，且由于需要大樣本空間的計(jì)算以提高精度，此方法受計(jì)算機(jī)條件限制，仿真精度和效率都難以令人滿意，無法滿足工程計(jì)算需要.

本文應(yīng)用蒙特卡羅法基于相關(guān)函數(shù)矩陣，仿真生成了Von Karman模型三維空間大氣紊流場，驗(yàn)證了該方法的正確性，且在算法上進(jìn)一步完善，仿真生成了精度較高的可應(yīng)用飛行動力學(xué)仿真的三維大氣紊流場.

1 Von Karman模型

Von Karman模型的能量頻譜函數(shù)為

其中，a=1.33;L為紊流尺度;Ω為空間頻率;σ為紊流強(qiáng)度.

由此可推得大氣紊流Von Karman模型的縱向和橫向相關(guān)函數(shù)如下:

其中，ζ=ξ/(aL);Γ為 Gamma函數(shù);K為 Bessel函數(shù).由于Von Karman模型函數(shù)形式復(fù)雜，使得常規(guī)仿真方法受限.

根據(jù)Bachelor公式，得出基于Von Karman頻譜的空間相關(guān)函數(shù)如下:

基于式(4)可構(gòu)造三維空間大氣紊流場的空間相關(guān)函數(shù)矩陣.

2 三維大氣紊流場建模理論

2.1 產(chǎn)生三維高斯白噪聲場并將其離散化

用數(shù)學(xué)庫函數(shù)產(chǎn)生(0，1)上高斯分布白噪聲序列{W}，隨機(jī)變量為{W1，W2，…，Wmnl}.

白噪聲的相關(guān)函數(shù)等于δ函數(shù):

將隨機(jī)變量Wi按逐行排列方式排布在劃分為 m ×n×r的空間網(wǎng)格上，m，n，r分別為沿 X，Y，Z軸方向的網(wǎng)格點(diǎn)數(shù).由此生成了如圖1所示的高斯分布三維白噪聲場{Z(x，y，z)}.其中 x=1，2，…，m，y=1，2，…，n，z=1，2，…，r，且 Z(i，j，k)=Wp，其中p=(k－1)×m ×n+(j－1)×m+i.

2.2 離散化三維空間隨機(jī)場的相關(guān)函數(shù)矩陣

構(gòu)造一般相關(guān)系數(shù)矩陣MU={Cij}={R(ij)}，按照圖1所示，空間劃分為p=m×n×r個(gè)網(wǎng)格，則點(diǎn)點(diǎn)相關(guān)構(gòu)造相關(guān)系數(shù)矩陣的一般表達(dá)式為

圖1 隨機(jī)變量空間分布示意圖

矩陣中Cij為由式(4)得到的基于Von Karman模型的空間相關(guān)函數(shù).

2.3 產(chǎn)生滿足相關(guān)函數(shù)矩陣的紊流場

由于MU為對稱正定陣，采用Cholesky分解法，可將其分解為一個(gè)上三角陣和一個(gè)下三角陣的乘積，即MU=AAT.其中，A為下三角陣.設(shè)U是需要產(chǎn)生的隨機(jī)場，取U=AZ，其中，Z=［Z1，Z2，…，Zn］T為白噪聲隨機(jī)變量.

可以證明U是滿足相關(guān)系數(shù)矩陣要求的隨機(jī)場.

以上三維大氣紊流場建模方法，理論完備，可以精確生成任意需要的三維大氣紊流場.在實(shí)際數(shù)值仿真計(jì)算中，仿真精度取決于輸入的白噪聲精度以及紊流的樣本空間大小.可以證明，白噪聲的白化程度越高，紊流場的樣本空間越大，仿真精度就越高，表現(xiàn)在仿真結(jié)果上就是仿真曲線和理論曲線貼合得越好.

3 仿真算法

本文采用蒙特卡羅法對三維大氣紊流場進(jìn)行數(shù)值仿真.蒙特卡羅法可以對復(fù)雜的隨機(jī)性問題進(jìn)行求解，方法及程序結(jié)構(gòu)簡單，大量簡單重復(fù)抽樣，做出的解答能否反映實(shí)際，與原來所建立的數(shù)學(xué)模型是否正確及輸入信息的質(zhì)量有關(guān)，而與本方法無關(guān).由于選用的樣本容量有限，因此總會存在誤差，提高解的精度則需要成百倍的增長樣本的容量.隨著計(jì)算機(jī)的迅速發(fā)展，將蒙特卡羅隨機(jī)模擬的方法應(yīng)用于研究三維大氣紊流場數(shù)值仿真中，使得對空間三維大氣紊流速度場的高精度仿真得以實(shí)現(xiàn).

3.1 改進(jìn)的雙隨機(jī)交換法提高輸入白噪聲精度

由于仿真計(jì)算中采用數(shù)學(xué)庫函數(shù)生成的偽隨機(jī)數(shù)存在固有缺陷，導(dǎo)致對紊流場的仿真精度存在很大影響［6］.為了提高偽隨機(jī)數(shù)的白化程度，本文采用改進(jìn)的雙隨機(jī)交換法.主要步驟如下:

1)隨機(jī)抽取序列x(i－1)(n)的兩點(diǎn)數(shù)據(jù)并互換其位置，得到序列xi(n)，重復(fù)執(zhí)行N次;

2)計(jì)算xi(n)的自相關(guān)函數(shù) ri和平方和SSi;

3)若SSi＜ε或i達(dá)到預(yù)定執(zhí)行運(yùn)算的最大循環(huán)次數(shù)Nmax，則停止運(yùn)行;

4)若 SSi＜ SSi－1，則返回步驟 1)，繼續(xù)執(zhí)行運(yùn)算，反之，放棄在步驟1)中進(jìn)行的第i次隨機(jī)雙換，并返回步驟1)繼續(xù)執(zhí)行.

5)重復(fù)執(zhí)行步驟1)～步驟4)j次，其中Nj=N/10j，若 Nj＜1，則令 Nj=1.

該方法只是對隨機(jī)數(shù)序列進(jìn)行重新組合，因此不改變隨機(jī)數(shù)樣本的方差、均值以及概率密度特性.由于步驟2)、步驟3)的計(jì)算量很大，上述方法用較小的計(jì)算量獲得了較大的白化效果.例如處理一個(gè)6 000點(diǎn)序列耗時(shí)僅幾分鐘，隨機(jī)序列白化程度明顯提高.

3.2 相關(guān)函數(shù)矩陣的計(jì)算

Von Karman模型的三維空間大氣紊流場計(jì)算，是對計(jì)算機(jī)速度和內(nèi)存要求較高的海量數(shù)據(jù)計(jì)算.在具有足夠大樣本容量的計(jì)算條件下，計(jì)算結(jié)果的統(tǒng)計(jì)特性將更接近理論值，數(shù)據(jù)結(jié)果的精度也就更高.在計(jì)算相關(guān)函數(shù)矩陣，并Cholesky分解為下三角矩陣的環(huán)節(jié)，將占用大量計(jì)算機(jī)內(nèi)存并耗費(fèi)大量時(shí)間.需采用通過降低程序計(jì)算量和優(yōu)化內(nèi)存使用的優(yōu)化算法來提高計(jì)算效率，使得在相同的計(jì)算條件下可以計(jì)算盡可能大的樣本空間.

以10×10×10的空間隨機(jī)場為例，構(gòu)造的基于Von Karman頻譜的縱向相關(guān)系數(shù)矩陣MU為

其中

C1，1000表示網(wǎng)格空間中第1點(diǎn)與第1000個(gè)點(diǎn)間的相關(guān)函數(shù)，其余類推;Δh是網(wǎng)格間距步長，為簡化表達(dá)式，直接按照式(4)計(jì)算出了系數(shù);L是大氣紊流縱向尺度.

以上Von Karman模型相關(guān)函數(shù)表達(dá)形式雖然復(fù)雜，但是由于本建模方法無需對模型表達(dá)式進(jìn)行簡化，因此可以進(jìn)行數(shù)值計(jì)算.

3.3 仿真大氣紊流場的特性檢驗(yàn)

三維大氣紊流場仿真計(jì)算結(jié)果是否合理，精度是否滿足要求，需要對仿真結(jié)果進(jìn)行相關(guān)特性檢驗(yàn).檢驗(yàn)的準(zhǔn)則就是看仿真紊流序列的相關(guān)特性是否符合該紊流模型的相關(guān)函數(shù)理論表達(dá)式［7］.Von Karman模型的縱向和橫向相關(guān)函數(shù)的理論表達(dá)式見式(2)和式(3).仿真序列的相關(guān)函數(shù)如下式:

將仿真的相關(guān)函數(shù)曲線和理論相關(guān)函數(shù)曲線畫在同一幅圖上，來衡量它們的貼合程度.

3.4 抽樣-檢驗(yàn)方法

蒙特卡羅方法提高解的精度需要成百倍地增加樣本的容量，而增加樣本容量受內(nèi)存容量和CPU計(jì)算速度的限制，這是蒙特卡羅方法通過增加采樣點(diǎn)數(shù)提高精度的主要限制因素.另一方面，結(jié)果的誤差難以估計(jì)，更是蒙特卡羅法固有的特點(diǎn).注意到蒙特卡羅法是一個(gè)純粹的抽樣方法，本文提出了“抽樣-檢驗(yàn)方法”:抽樣方法采用蒙特卡羅法，檢驗(yàn)則是對抽樣結(jié)果的誤差進(jìn)行統(tǒng)計(jì)計(jì)算，進(jìn)而篩選誤差達(dá)到要求的仿真紊流場.本文對檢驗(yàn)長度分段逐段檢驗(yàn)來實(shí)現(xiàn)，對每段分別設(shè)定誤差，以實(shí)現(xiàn)精細(xì)的誤差控制，從而得到滿足相應(yīng)要求的仿真結(jié)果.按照抽樣-檢驗(yàn)方法的要求改進(jìn)后的算法，可以在短時(shí)間內(nèi)完成同一尺度的成千上萬個(gè)紊流場計(jì)算，按照數(shù)理統(tǒng)計(jì)理論的大數(shù)定律，可以從中挑選到相關(guān)特性與Von Karman模型吻合度很高的紊流場，這樣就不需要成百倍地增加樣本的容量而得到更高精度的仿真結(jié)果.圖2為仿真算法流程示意圖.

圖2 仿真算法流程示意

4 仿真算例和分析

選取兩個(gè)算例，用本文方法計(jì)算基于Von Karman模型的三維紊流場.這兩個(gè)算例的紊流場尺度，在常見的紊流場的大尺度范圍內(nèi)，具有典型的工程應(yīng)用意義.選取紊流強(qiáng)度σ=1m/s，縱向紊流尺度為L=530m.根據(jù)反復(fù)計(jì)算調(diào)試，設(shè)定比較合理的誤差率控制條件，以達(dá)到精度和計(jì)算效率的雙贏.

算例1 網(wǎng)格點(diǎn)數(shù):6000(1500×2×2).檢驗(yàn)長度:100點(diǎn)，分3段(40點(diǎn)，30點(diǎn)，30點(diǎn))逐段檢驗(yàn).U序列的控制誤差率為:第1段，0.05;第2段，0.05;第3 段，0.1.V，W 序列的控制誤差率為:第1段，0.1;第 2段，2;第 3段，1.計(jì)算時(shí)間:48min.相關(guān)性檢驗(yàn)結(jié)果見圖3.

圖3 三維大氣紊流場(1500×2×2)的相關(guān)性檢驗(yàn)

算例2 網(wǎng)格點(diǎn)數(shù):8000(2×2000×2).檢驗(yàn)長度:100點(diǎn)，分3段(40點(diǎn)，30點(diǎn)，30點(diǎn))逐段檢驗(yàn).U序列的控制誤差率為:第1段，0.05;第2段，0.05;第3 段，0.1.V，W 序列的控制誤差率為:第1 段，0.1;第 2 段，1.5;第 3 段，0.75.計(jì)算時(shí)間:0.25min.相關(guān)性檢驗(yàn)結(jié)果見圖4.

圖4 三維大氣紊流場(2×2000×2)的相關(guān)性檢驗(yàn)

圖5是算例1的Von Karman模型紊流速度序列仿真結(jié)果，圖6是仿真紊流場(30×30)剖面示意圖.

圖5 Von Karman模型紊流速度序列仿真結(jié)果

本方法未對模型進(jìn)行任何簡化，因此仿真生成的三維空間大氣紊流場完全符合所選用的Von Karman大氣紊流模型.由圖3～圖6的相關(guān)性檢驗(yàn)圖、仿真紊流速度圖及紊流場剖面圖等表明，仿真相關(guān)性曲線與理論曲線吻合程度高，與以往其他研究結(jié)果相比，仿真精度有了大幅度提高，非常令人滿意.另外本文通過具體算法的改進(jìn)，在計(jì)算條件、耗時(shí)和結(jié)果逼真度之間找到了一個(gè)比較理想的平衡方法，計(jì)算時(shí)間大大縮短，使仿真效率得到提高.

圖6 仿真三維空間大氣紊流場剖面示意圖

5 結(jié)論

通過大量計(jì)算分析表明，應(yīng)用蒙特卡羅法基于相關(guān)函數(shù)矩陣，仿真生成Von Karman模型三維空間大氣紊流場的模型方法正確，仿真結(jié)果合理，計(jì)算精度和效率高.

本文的仿真理論，實(shí)現(xiàn)了對復(fù)雜的Von Karman模型的精確推導(dǎo)和應(yīng)用，使得仿真結(jié)果在理論上可以無限逼近理論值.在仿真方法上，改進(jìn)的雙隨機(jī)交換法提高了白噪聲的精度，優(yōu)化了蒙特卡羅方法，采用抽樣-檢驗(yàn)方法，實(shí)現(xiàn)了精細(xì)的誤差控制，從而解決了蒙特卡羅方法提高解的精度需要海量樣本容量和誤差控制的困難.本文提出的仿真方法，大大提高了仿真精度和效率.

本方法可直接應(yīng)用于飛行仿真，使飛行仿真中的風(fēng)場模擬更加真實(shí)，對模擬器的研制具有實(shí)際意義，完全符合工程應(yīng)用需要.在實(shí)時(shí)飛行仿真中，可對白噪聲和確定紊流空間的相關(guān)函數(shù)矩陣進(jìn)行預(yù)計(jì)算，以滿足實(shí)時(shí)仿真需要.如計(jì)算條件有限，又需要較大紊流場，可采用基于對稱原理的紊流場擴(kuò)展方法，將生成的局部空間紊流場通過拼接的辦法有效擴(kuò)展為大范圍連續(xù)紊流場.

References)

［1］ Etkin B.Turbulent wind and its effect on flight［J］.Journal of Aircraft，1981，18(5):327 －345

［2］肖業(yè)倫.用于飛行仿真的二維紊流場的數(shù)字生成法［J］.航空學(xué)報(bào)，1990，11(4):124 －129 Xiao Yelun.Digital generation of two-dimensional field of turbulence for flight simulation［J］.Acta Aeronautica et Astronautica Sinica，1990，11(4):124 －129(in Chinese)

［3］陸宇平，胡亞海.基于空間相關(guān)函數(shù)的二維紊流場數(shù)值生成法［J］.南京航空航天大學(xué)學(xué)報(bào)，1999，31(2):139－145 Lu Yuping，Hu Yahai.Digital generation of two-dimensional field of turbulence based on spatial correlation function［J］.Journal of Nanjing University of Aeronautics＆ Astronautics，1999，31(2):139－145(in Chinese)

［4］張峰，汪沛，王沖，等.基于Von Karman模型的三維大氣紊流仿真［J］.計(jì)算機(jī)仿真，2007，24(1):35 －38 Zhang Feng，Wang Pei，Wang Chong，et al.Simulation of threedimensional atmospheric turbulence based on Von Karmanmodel［J］.Computer Simulation，2007，24(1):35 － 38(in Chinese)

［5］洪冠新，肖業(yè)倫.蒙特卡羅法仿真生成三維空間大氣紊流場［J］.航空學(xué)報(bào)，2001，22(6):542 －545 Hong Guanxin，Xiao Yelun.Monte Carlo simulation for3-D-field of atmospheric turbulence［J］.Acta Aeronautica et Astronautica Sinica，2001，22(6):542 －545(in Chinese)

［6］馬樹峰，岳曉奎.大氣紊流的數(shù)字仿真［J］.計(jì)算機(jī)輔助工程，2008，17(3):54 －56 Ma Shufeng，Yue Xiaokui.Digital simulation on atmospheric turbulence［J］.Computer Aided Engineering，2008，17(3):54 －56(in Chinese)

［7］肖業(yè)倫.大氣擾動中的飛行原理［M］.北京:國防工業(yè)出版社，1993:166－173 Xiao Yelun.Flight theory in atmosphere turbulence［M］.Beijing:National Defense Industry Press，1993:166 － 173(in Chinese)

(編 輯:李 晶)

Theory and method of numerical simulation for 3D atmospheric turbulence field based on Von Karman model

Gao Jing Hong Guanxin Liang Zaoqing
(School of Aeronautic Science and Engineering，Beijing University of Aeronautics and Astronautics，Beijing 100191，China)

atmospheric turbulence;Monte Carlo method;Von Karman model;numerical simulation;flight simulation

V 321.2

A

1001-5965(2012)06-0736-05

2011-03-28;網(wǎng)絡(luò)出版時(shí)間:2012-04-01 11:59

http://www.cnki.net/kcms/detail/11.2625.V.20120401.1159.003.html

凡舟青年科研基金資助項(xiàng)目(20080502)

高 靜(1981－)，女，陜西子洲人，助理研究員，gaojing@buaa.edu.cn.

Three-dimension-field atmospheric turbulence based on Von Karman model was generated with Monte Carlo simulation，in view of correlation function matrix.According to the derivation of the Von Karman model for mulae，the correlation curves of the simulated field can highly approach the theory results.The precision and efficiency of atmospheric turbulence simulation were improved by using double stochastic interchange minimization algorithm to improve Gauss white noise sequence，and also by the optimization of Monte Carlo method.The simulated atmospheric turbulence field，is identical with the Von Karman model by numerical test，and has perfect statistical characteristic and simulation precision.The simulation method could be used in real-time flight simulation.