適于直升機飛行仿真的高原大氣紊流模型

吉洪蕾，蘇俊杰，陳仁良，孔衛紅

1. 重慶大學 航空航天學院，重慶 400044 2. 南京航空航天大學 航空學院 直升機旋翼動力學國家級重點實驗室，南京 210016

中國西部高原地區具有重要的國防戰略地位，其海拔高、地勢險峻、大風天氣多，獨特的地理氣候環境導致大氣紊流出現的概率以及紊流的強度和頻率寬度都遠高于一般平原地區，直升機在該地區執行近地飛行任務時受到大氣紊流的干擾，飛行品質下降、甚至危及飛行安全，準確掌握該地區大氣紊流對直升機運動的動態干擾特性顯得十分重要且迫切。然而，目前仍然缺乏一種適于直升機飛行仿真的高原大氣紊流模型。

大氣紊流對航空器飛行品質的影響始終是航空界重點關注的飛行安全問題之一。經過多年發展，用于固定翼飛機飛行仿真的大氣紊流建模技術已經日趨成熟。與之相比，適于直升機飛行仿真的大氣紊流模擬技術則略有滯后，主要原因之一是直升機升力主要由旋翼的旋轉產生，大氣紊流對直升機的干擾機理迥異于固定升力體的飛機，適于固定翼飛機的大氣紊流建模方法對于直升機來說不再合適。為解決這一問題，Gaonkar等提出槳葉紊流建模方法，該方法考慮各片槳葉氣動中心的相關性生成大氣紊流速度，并將其疊加各槳葉相對來流速度上建立用于直升機飛行仿真的大氣紊流模型。該方法雖然符合大氣紊流對直升機旋翼的氣動干擾機理，但是其計算效率不適合實時仿真應用。Mcfarland和Duisenberg提出了旋翼葉素紊流建模方法。該方法在旋翼前方固定兩個大氣紊流成型濾波器，基于“凍結場”假設，隨著直升機前飛則大氣紊流速度向后擴展，形成覆蓋整個旋翼的大氣紊流速度場。吉洪蕾等拓展了旋翼葉素紊流方法，在直升機前方垂直平面內布置多個紊流成型濾波器，基于“凍結場”假設，隨直升機前飛形成三維空間紊流模型。該方法由于兼顧廣泛的適用性和較高的計算效率而更加適于直升機飛行仿真應用，但是其采用低精度的德萊頓紊流譜模型進行數值模擬，且沒有考慮到橫切面內大氣紊流速度的空間相關性，因而計算精度受到限制。

除了以上基于物理機理的模型，控制等效紊流輸入(Control Equivalent Turblence Input,CETI)模型提供了一種數據驅動的直升機大氣紊流建模方法。其以具有特定頻譜的信號疊加到操縱輸入產生與大氣紊流干擾等效的直升機頻域響應，模型參數由飛行試驗數據辨識。CETI模型受益于其基于飛行試驗數據的特點，具有高計算精度且滿足實時仿真要求，但是也受限于飛行試驗數據缺乏，僅適合懸停、小速度飛行狀態。

直升機由于其獨特的飛行特點，常常在低空執行機動任務，經歷的大氣紊流環境常常比固定翼飛機惡劣得多，影響大氣紊流統計特性的因素也多。我國西部高原地區地貌復雜多變、晝夜溫差大、大風天氣多的特點導致其大氣紊流活動更加活躍，低空大氣紊流受下墊面地形地貌、大氣層結及平均風速的影響更加顯著，固定翼飛機飛行品質規范給出的僅隨高度變化的紊流統計參數無法反映以上因素的影響，成為限制高原大氣紊流環境直升機飛行仿真及飛行品質分析的主要問題之一。

本文發展一種考慮我國西部高原地形地貌、大氣層結和平均風速廓線影響的高原大氣紊流模型，針對已有三維空間大氣紊流模型精度低、無法考慮橫切面內紊流速度空間相關性的問題，采用高精度的馮·卡門紊流譜函數推導大氣紊流速度的高階成型濾波器及其時間序列生成算法，提出空間相關的隨機大氣紊流速度時間序列生成方法解決橫切面內紊流速度的空間相關性問題，提升大氣紊流模擬的計算精度，并通過數值仿真和與試驗數據對比驗證模型。

1 適于飛行仿真的高原大氣紊流模型

1.1 高原大氣紊流模型

低空大氣紊流，尤其是600 m高度以下的大氣紊流受下墊面地形地貌、大氣層結及平均風速的影響顯著。因此，首先建立平均風速模型。平均風速是指在一定時距內風速變化的平均值。600 m高度以下大氣邊界層內的平均風速度受到地表摩擦力和大氣層結狀態的影響，風速隨高度增加而增大。對于均一地貌下墊面，采用地面粗糙度表征地面摩擦力對風速的影響，指數律描述的大氣邊界層內風速廓線為

(1)

式中：為高度的平均風速度;為參考高度(通常為10 m)的風速;為風速廓線指數。

在高于600 m的大氣邊界層外，可以認為風速隨高度不再變化，即

=≥600

(2)

式中：為600 m高度的平均風速度。

風速廓線指數受到下墊面粗糙度、大氣層結狀態等因素的影響，采用式(3)計算:

=+lg+(lg)

(3)

式中：、、為風速廓線指數系數，其取值隨大氣層結狀態變化。國外研究人員測量了多個地區不同大氣層結狀態的風速廓線指數，采用測試數據擬合的方式獲得如表1所示的風速廓線指數系數。

表1 風速廓線指數系數取值[18]Table 1 Coefficients for wind speed profile exponent[18]

大量實測數據表明，高原大氣紊流速度的無量綱頻譜在慣性子區滿足-2/3次方衰減的卡莫夫條件，在低頻段滿足趨于零的極限條件，與馮·卡門紊流譜模型的基本理論一致，并且我國多個地區大氣紊流速度實測頻譜也與馮·卡門譜函數吻合良好。因此，采用馮·卡門紊流譜模型描述高原大氣紊流速度隨機變化的頻域特征。基于大氣紊流速度的一維空間頻譜函數以及“凍結場”假設，由角頻率和空間頻率形式的紊流譜轉化關系，可以得到角頻率形式的大氣紊流譜函數為

(4)

式中：()、()、()為縱、橫和垂向紊流速度譜函數;,,為3個方向紊流強度;=(為速度;為積分尺度)為特征頻率。在紊流強度保持一定的條件下，特征頻率決定了紊流能量的頻域分布，隨著值的增大，大氣紊流速度呈現出更多高頻成分。

由式(4)的紊流譜函數可知，大氣紊流強度和積分尺度是大氣紊流速度隨機變化的主要統計特性參數，通過調整兩者的變化即可描述不同頻率分布和幅值的大氣紊流。在300 m高度以下的大氣邊界層內，紊流活動發展強盛，并受地面摩擦、大氣層結狀態等因素的影響，隨高度、地面粗糙度、風速等因素變化，紊流強度和積分尺度采用國外統計的近地層變化規律來描述:

(5)

(6)

在600 m高度以上的自由層，大氣紊流不再受到下墊面的影響，可以認為其強度和積分尺度不再隨高度變化，紊流強度和積分尺度根據飛行品質規范給出

===01

(7)

===280

(8)

在300～600 m高度，紊流強度和積分尺度由大氣邊界層內取值到自由大氣層取值逐漸過渡:

(9)

(10)

式中：300、300、300為300 m高度的大氣紊流強度，其表達式分別為

由式(3)、式(5)和式(6)可以看出，高原大氣紊流模型考慮到了下墊面地貌和大氣層結狀態引起的平均風速度、紊流強度和積分尺度變化及其對大氣紊流統計特性的影響規律。

1.2 大氣紊流速度生成算法

由于馮·卡門紊流譜函數的無理式形式，直接將其用于生成大氣紊流速度序列較為困難。固定翼飛機飛行品質規范針對馮·卡門譜給出了一種三階有理形式的擬合結果，該結果的有效無量綱頻率為50 rad，對于飛行速度較高的固定翼飛機是合適的。與之相比，直升機飛行速度遠低于固定翼飛機，三階有理形式的擬合結果已無法滿足直升機懸停/低速狀態飛行品質評估精度的要求。采用最小二乘法對式(4)的馮·卡門譜擬合，獲得四階有理形式的大氣紊流譜函數，共軛分解之后得到高階大氣紊流成型濾波器如式(11) 所示:

(11)

式中：()、()、()分別為縱、橫和垂向高階紊流成型濾波器。系數,,,分別為

針對式(11)的高階紊流成型濾波器，采用離散變換的方式獲得大氣紊流速度時間序列的離散遞推生成算法。由于零極點匹配法離散精度高，甚至在奈奎斯特頻率附近也能較好地保持頻域離散精度，采用該方法對式(11)的縱向紊流濾波器()離散后得到

Δ+3=2Δ+2+1Δ+1+0Δ+(+2+1+1+0)

(12)

式中：和Δ分別表示=Δ時刻的白噪聲輸入和縱向紊流速度輸出。式(12)各系數分別為

且縱向紊流成型濾波器的零點和極點為

式中：=。

同理，可以得到橫、垂向大氣紊流速度時間序列生成算法為

Δ+4=3Δ+3+2Δ+2+1Δ+1+0Δ+

,(+3+2+2+1+1+0)

(13)

式中：Δ表示在=Δ時刻的橫向或垂向紊流速度分量。式(13)各系數分別為

式中:橫、垂向紊流速度成型濾波器的零點和極點為

其中：,=,。

由式(12)和式(13)的離散遞推算法可知，該大氣紊流速度時間序列生成算法計算簡單、效率高并具備較高的紊流頻譜，更加適于低空、低速飛行的直升機飛行仿真研究。

1.3 適于飛行仿真的三維空間紊流模型

基于大氣紊流速度時間序列，生成覆蓋整個直升機的三維空間紊流場，如圖1中長方體所示。為了計算每一時刻直升機各氣動面處的大氣紊流速度，定義了兩個坐標系：一個是地軸系，其中軸指向地球正北方向為正，軸指向地球正東方向為正，軸指向地心為正；另一個是紊流坐標系，其中軸指向空速的反方向，軸垂直于與軸指向機身左側為正，軸垂直于平面指向機身上方為正。紊流場長度為，寬度為2，高度為，固定于機體并隨直升機以速度前飛。紊流場前表面與旋翼槳尖平面相切，且垂直于直升機前飛速度。基于“凍結場”假設，將=×組大氣紊流成型濾波器等間隔置于紊流場前表面，每套成型濾波器生成3個分量的大氣紊流速度時間序列。當直升機以速度前飛時，由于大氣紊流相對地面靜止，所以置于前表面的紊流濾波器生成的大氣紊流時間序列將以一定的延遲時間向直升機尾部拖出，布滿整個大氣紊流場，形成適于直升機飛行仿真的三維空間大氣紊流場。

圖1 三維空間紊流模型及坐標系定義Fig.1 Three-dimensional turbulence model and definition of coordinate axes

由于式(12)和式(13)的大氣紊流速度時間序列生成算法是基于馮·卡門一維紊流譜函數獲得的，因此在沿縱向的紊流速度時間序列滿足縱向相關性，但無法滿足各橫截面內的相關性要求。隨機統計理論早已證明，基于協方差矩陣Cholesky因子分解的線性變換方法對多維正態分布的相關變量隨機數序列生成效果最好。因此，采用協方差矩陣Cholesky因子分解變換法將表面上相互獨立的紊流成型濾波器相互關聯。

假設置于表面的紊流成型濾波器分為行和列均勻排列。任意兩個互相關的成型濾波器生成的縱向紊流速度序列之間的相關系數為

(14)

由式(14)可以得到所有組成型濾波器的相關系數矩陣=[(,)]，對其進行Cholesky分解可得下三角矩陣，即

=

(15)

Δ=Δ

(16)

式中：Δ為×組相互關聯的大氣紊流速度隨機序列向量;Δ為×組獨立紊流速度序列組成的向量。其表達式為

同理，可以生成×組互相關的橫、垂向紊流速度時間序列，進而形成×組互相關的大氣紊流成型濾波器。

將以上互相關的大氣紊流成型濾波器布于大氣紊流場前表面上，結合“凍結場”假設,建立的三維空間大氣紊流場沿空間3個方向的分布均滿足馮·卡門模型的空間相關性要求。

2 集成大氣紊流的直升機飛行動力學模型

將上述三維空間大氣紊流模型集成到已有高階非線性直升機飛行動力學模型以開展飛行仿真研究，詳細的飛行動力學模型描述及驗證參見文獻[12-13,25]。大氣紊流速度主要通過干擾直升機各氣動面的相對氣流速度影響直升機氣動載荷及運動。由地軸系到紊流坐標系的坐標變換矩陣為

(17)

(18)

(19)

式中：、、為空速在地軸系的縱、橫和垂向分量。

以機身為例計算直升機各氣動面在大氣紊流坐標系的坐標。考慮到大氣紊流場固定于直升機重心，機身在紊流坐標系中的坐標為

(20)

(21)

采用類似方法可計算出旋翼各葉素、平尾、垂尾和尾槳處大氣紊流速度，采用線性疊加方式將其計入到直升機各氣動面的相對氣流速度上，形成集成高原大氣紊流的直升機飛行動力學模型。以4片槳葉的黑鷹直升機為例，飛行動力學模型可以表示為

(22)

3 模型驗證

3.1 高原大氣紊流模型驗證

采用四川省理塘縣氣象觀測站實測的風速廓線和大氣紊流度驗證高原大氣紊流模型，氣象站所在區域海拔高度為3 932 m，地表為高原草甸，下墊面粗糙度約為0.05 m。

圖2為理塘觀測站實測風速廓線與模型預測結果的對比。觀測數據分別在春季某天的02時、12時、18時和20時測得，預測風速廓線的大氣層結狀態分別為弱不穩定和中性狀態，參考高度10 m 的平均風速分別為6.5 m/s、4.9 m/s、4.7 m/s 和3.1 m/s。從圖2中可以看出，上午由于日照的逐漸增加，大氣呈現出弱不穩定狀態；下午到晚上由于日照逐漸減弱，大氣層結狀態趨于穩定，處于中性狀態。按照給定的下墊面粗糙度和大氣層結狀態，風速廓線模型能夠捕捉到高原大氣層結狀態和離地高度對平均風速的影響，預測結果可以與實測風速吻合的很好。

圖2 理塘模型預測風速廓線與實測結果對比Fig.2 Comparison between predicted and measured wind speed profiles in Litang

大氣紊流速度雖然呈現出很強的非定常性，但是其統計結果具有特定的分布規律。圖3為理塘觀測站實測大氣紊流度與模型預測值的對比結果。從圖3中可以看出，大氣紊流度的觀測結果存在多變性，尤其是在10 m以下高度，由于受到局部地貌和大氣狀態的干擾，大氣紊流度的實測結果散布范圍較大；隨著高度增加，同一高度處的大氣紊流度實測值趨于一致，呈現出隨高度降低的規律。雖然觀測結果存在多變性，但模型預測結果位于各散點中間，尤其是對直升機氣動特性影響最大的垂向紊流速度分量，模型能夠較好地預測大氣紊流度的統計值及其變化趨勢。由于大氣紊流的隨機性，采用具有代表性的統計特性參數開展大氣紊流對直升機飛行仿真和飛行品質的影響研究具有重要參考價值。

圖3 理塘模型預測大氣紊流度與實測結果對比Fig.3 Comparison between predicted and measured turbulence intensities in Litang

生成時域大氣紊流速度時間序列，計算其功率譜密度并與馮·卡門理論譜對比驗證高原大氣紊流模型的紊流速度生成算法。驗證以四川理塘和西藏日喀則兩地為例，其中，日喀則平均海拔超過4 000 m，地表主要是凹凸不平的丘陵地帶，下墊面粗糙度約為0.3 m。

圖4為仿真紊流速度頻譜與馮·卡門理論頻譜的對比結果。其中，參考高度10 m平均風速度為11.5 m/s，大氣層結狀態為中性，離地高度為12 m。在同樣參考平均風速和大氣層結狀態的條件下，由于地形地貌的差別，日喀則地區粗糙度(=0.3 m)更高，低空大氣紊流的幅值和頻率寬度都遠大于下墊面粗糙度更低(=0.05 m)的理塘地區大氣紊流。仿真得到的大氣紊流速度頻譜與理論馮·卡門頻譜模型吻合良好，證明大氣紊流仿真算法能夠準確捕捉高原地形地貌變化對大氣紊流速度特性的影響。

圖4 大氣紊流速度仿真頻譜與馮·卡門理論頻譜對比Fig.4 Comparison between simulated and theoretical von Krmn atmospheric turbulence velocity spectra

3.2 三維空間大氣紊流模型驗證

以黑鷹直升機為例生成西藏日喀則地區上空的三維空間大氣紊流驗證仿真模型。黑鷹直升機旋翼半徑為8.18 m，長度為19.76 m，高度為4.05 m。采用80個空間相關的紊流成型濾波器布于圖1前表面生成三維空間大氣紊流場，其中橫向分布個數=40，垂向分布個數=2，參考高度10 m的平均風速度=11.5 m/s，日喀則地區的下墊面粗糙度=0.3 m，大氣紊流場所在高度=12 m。

圖5為某一時刻旋翼所在平面的垂向大氣紊流速度場Δ。從圖中可以看出，雖然每個空間位置的垂向大氣紊流速度呈現出很強的隨機性，但是整個大氣紊流速度場表現出一定的空間相關性。

圖5 垂向紊流速度場仿真結果Fig.5 Simulated vertical turbulence velocity field

圖6為垂向大氣紊流速度場不同濾波器生成的時域紊流速度功率譜密度與馮·卡門譜的對比。從圖6中可以看出，雖然空間相關紊流成型濾波器生成的時域紊流速度序列并不完全相同(如圖5所示)，但其頻譜都與理論值一致，證明各空間相關紊流成型濾波器全部滿足馮·卡門模型的縱向相關性要求。

圖6 空間相關垂向紊流速度仿真頻譜與理論譜對比Fig.6 Comparison between simulated and theoretical spectra of spatially collected vertical turblence velocity

圖7為垂向大氣紊流速度場沿橫向的相關系數與理論值的對比。從圖7中可以看出，由各紊流成型濾波器生成的大氣紊流速度時間序列滿足空間相關性的理論要求，證明了空間相關紊流成型濾波器生成方法的準確性。

圖7 垂向紊流速度相關系數仿真結果驗證Fig.7 Validation of simulated vertical turbulence velocity correlation coefficients

3.3 直升機對大氣紊流的頻域響應驗證

以UH-60A黑鷹直升機為例模擬直升機對大氣紊流的響應，計算其功率譜密度，并與飛行試驗數據對比驗證集成大氣紊流的直升機飛行動力學模型。其中，飛行試驗為UH-60A黑鷹直升機在機庫后方大氣紊流環境中保持12 m高度懸停，平均風速度為11.6 m/s，垂向大氣紊流強度為1.37 m/s，機庫高度為12 m。飛行仿真的參考高度10 m平均風速度=11.5 m/s，根據機庫高度，下墊面粗糙度設為0.3 m，大氣紊流場由80個空間相關的紊流成型濾波器生成，其中橫向分布個數=40，垂向分布個數=2。飛行仿真中采用駕駛員模型實現直升機在大氣紊流干擾環境的姿態穩定和位置保持，駕駛員模型的詳細設計參見文獻[25]，姿態穩定回路的穿越頻率為2.5 rad/s。

圖8為飛行仿真的直升機對大氣紊流頻域響應結果與飛行試驗數據的對比。由于僅有3個方向角速度和垂向速度的飛行試驗數據，因此與試驗的對比也僅限于直升機對大氣紊流響應的角速度和垂向速度。從圖8中可以看出，雖然由于大氣紊流的干擾直升機的響應特性隨時間隨機變化，但直升機對大氣紊流的頻域響應是固定分布的，模型預測的直升機角速度和垂向速度功率譜密度與飛行試驗數據吻合良好，集成大氣紊流的飛行動力學模型能夠有效捕捉直升機對大氣紊流的頻域響應特性。

圖8 直升機對大氣紊流的頻域響應驗證Fig.8 Validation of helicopter frequency response to turbulence

4 結 論

1) 高原大氣紊流模型能夠準確捕捉高原地形地貌和大氣層結對平均風速廓線的影響，大氣紊流速度生成算法能夠反映下墊面粗糙度和平均風速的影響生成符合理論馮·卡門頻譜分布的時域大氣紊流速度。

2) 空間相關的大氣紊流成型濾波器仍然保持馮·卡門頻譜分布特性，各個成型濾波器生成的時域大氣紊流速度之間的相關系數與理論值吻合良好，采用空間相關成型濾波器和“凍結場”假設擴展的三維空間隨機紊流速度場滿足空間相關性要求。

3) 與飛行試驗數據的對比表明，集成大氣紊流的直升機飛行動力學模型能夠準確模擬算例直升機懸停狀態對大氣紊流的頻域響應特性，模型在直升機飛行仿真及飛行品質分析方面具有潛在的應用前景。