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

高 靜 洪冠新 梁灶清

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

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

高 靜 洪冠新 梁灶清

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

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

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

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

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

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

1 Von Karman模型

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

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

由此可推得大氣紊流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}，隨機變量為{W1，W2，…，Wmnl}.

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

將隨機變量Wi按逐行排列方式排布在劃分為 m ×n×r的空間網(wǎng)格上，m，n，r分別為沿 X，Y，Z軸方向的網(wǎng)格點數(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 離散化三維空間隨機場的相關(guān)函數(shù)矩陣

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

圖1 隨機變量空間分布示意圖

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

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

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

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

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

3 仿真算法

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

3.1 改進的雙隨機交換法提高輸入白噪聲精度

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

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

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

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

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

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

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

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

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

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

其中

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

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

3.3 仿真大氣紊流場的特性檢驗

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

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

3.4 抽樣-檢驗方法

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

圖2 仿真算法流程示意

4 仿真算例和分析

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

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

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

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

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

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

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

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

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

5 結(jié)論

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

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

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

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

凡舟青年科研基金資助項目(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.