結冰條件下飛機全包線模態特性分析方法

伍 強, 徐浩軍, 裴彬彬,*, 朱和銓, 武 欣

(1. 空軍工程大學航空工程學院, 陜西 西安 710038;2. 中國人民解放軍93705部隊, 河北 唐山 064200)

0 引 言

波音公司統計結果表明,飛機結冰是誘發飛行失控嚴重事故的三大因素之一[1],飛機結冰問題長期以來一直受到人們的高度關注。據國際航空器擁有者及駕駛員協會的一項調查表明,1990—2000年間在所有天氣導致的飛行事故中,與結冰有關的飛行事故占了12%[2]。

結冰將導致飛機氣動特性和動力學特性惡化,進而影響飛行品質。早在1929年,Kopp就已指出,結冰導致的飛機空氣動力學特性惡化相較于飛機重量增加對飛行安全的影響更顯著[3]。Preston等人最早于1948年采用飛行試驗獲得了C-46飛機在自然積冰條件下阻力的變化特性[4]。美國宇航局Lewis研究中心與美國聯邦航空局技術中心聯合開展了尾翼積冰研究計劃[5],利用“雙水獺”(DHC-6)飛機開展了一系列的帶模擬冰型飛行試驗[6]和結冰風洞試驗[7],為積冰研究積累了大量的氣動數據。王良禹等人通過數值仿真對飛機結冰后的滾轉角變化進行了分析發現結冰后飛機飛行品質發生了降級[8]。研究結冰對飛行品質的影響可以為評估飛行安全提供依據[9]。但是現有的結冰后飛機飛行品質影響研究多集中在針對某些特定的狀態點進行分析,鮮有對全包線范圍內飛行品質變化的研究。

現代飛機均多采用高增穩系統,高階的飛機動力學模型使飛行品質的研究較為困難,通過等效系統擬配可將高階模型轉化為等效的低階模型[10]。侯世芳等人提出根據飛機在擾動作用下所固有的模態特性及其模態表達式來對飛機的時域響應基于最小二乘法進行擬合,從而計算出模態參數的方法[11]。最小二乘法憑借其原理簡單,編程容易且運算效率快的優點成為飛行品質分析中最常采用的低階等效系統參數擬配方法[12]。但其對初值的依賴性大的問題也十分明顯。傅慶慶等人在最小二乘法對初值依賴性研究中發現各參數的結果值受初值影響較大,要想獲得更好的擬配效果,需要經過多次試錯來得到最合適的初值[13]。

現有的文獻中對于等效擬配過程中的關鍵環節——擬配初值的選取問題鮮有專門的文獻進行敘述。在工程實際中,主要有兩種方法:一種是在合理的范圍內不斷地試錯出擬配初值,直至計算出合理的結果,這種方法無疑是費時費力的,而且一旦飛行狀態或操縱動作改變,初值就要重新設定;另一種方法是采用遺傳算法等智能優化算法獲取擬配初值,但從本質上來說,雖然遺傳算法不需要給出初值,但算法本身是通過隨機生成多組初值,再進行迭代進化計算的。從計算效率上來看,其計算效率并不高,而且所計算出來的初值僅是數學意義上的優化解,甚至可能出現不滿足實際情況情形。上述問題限制了等效擬配的使用,往往只能用于某個特定狀態點下進行品質分析,無法對整個包線范圍飛行品質進行評估。

本文針對上述問題,開發了自動數據預處理模塊,針對待評估的飛行數據的響應特征,初步計算出待擬配特征參數的初值,改善了現有的等效擬配中初值獲取方法無法涵蓋所有情形的問題,可適應不同響應曲線的等效擬配。在該方法的基礎上,開展了某型飛機全包線范圍內橫航向品質分析,驗證了該方法對于自適應計算飛機飛行品質具有廣泛的有效性。此外,分析了結冰條件下飛行包線和飛行品質的變化,獲取了結冰對整個飛行包線范圍內飛行品質的量化影響。

1 模型的建立

1.1 飛機本體模型的建立

飛機本體非線性動力學模型可表示為

(1)

式中:x為狀態向量,包含飛行速度、迎角、側滑角、四元數、俯仰角速率、滾轉角速率、偏航角速率和空間位置參數:

x=[V,α,β,q0,q1,q2,q3,p,q,r,xg,yg,zg]T

(2)

u為控制向量,包括油門偏度指令、升降舵偏度指令、副翼偏度指令和方向舵偏度指令:

u=[δth,δe,δa,δr]T

(3)

本文采用四元數法求解運動方程,能夠有效避免出現奇點,同時能夠減少三角函數的計算,提高運行效率[14],具體的方程表達形式參考文獻[14]。

1.2 電傳操縱系統模型的建立

電傳操縱系統采用大規模集成電路等數字技術,將駕駛員操縱指令經過運算處理(即飛行控制律)以電信號的形式傳遞給電液舵機,本文基于Simulink平臺,對某型飛機電傳操縱系統進行建模研究。

通常電傳操縱系統包括縱向操縱、橫向操縱與航向操縱3個主要通道,以及機動襟翼、極限狀態限制器等輔助操縱通道。某型飛機橫向通道結構圖如圖1所示。

圖1 橫向電傳操縱系統結構圖Fig.1 Structure diagram of lateral telex control system

其中,K1是位移信號按動壓修正的傳遞系數,K2是桿位移信號按動壓修正的傳遞系數,K3是p按迎角變化的傳遞系數,在迎角大于或等于某個數值后自動切斷傾斜阻尼反饋通道,K4是用動壓和靜壓修正的傾斜阻尼的傳遞系數,K5是按法向過載限幅的傳遞系數,K6是按動壓限幅的傳遞系數。num1/den1是對p信號中飛機結構彈性振動頻率信號濾波的濾波器,num2/den2是慣性濾波器,過濾掉駕駛桿傳來的噪聲信號。

根據橫向操縱系統結構圖,基于Simulink軟件平臺[15],建立橫向操縱系統的模型。縱、航向操縱系統和控制律的構建方法與橫向的相似,在此不再贅述。

1.3 非線性結冰影響模型的建立

結冰后動力學特性分析的前提首先要建立結冰影響模型,但由于結冰的復雜性和隨機性,直到目前仍然沒有能夠準確預測結冰后氣動參數的通用模型。目前只能通過冰風洞試驗、飛行試驗或者數值計算獲得給定氣象條件下的結冰冰形,然后通過風洞測力、氣動參數辨識或計算流體力學數值模擬等手段獲得該冰形影響下的氣動數據。最后在此基礎上構建不同冰形的氣動數據庫,根據飛行狀態進行插值計算出當前狀態下的氣動參數[16-20]。當前常用的是伊利諾伊州大學Bragg教授提出的結冰參量影響模型[16],模型可表示為

C(A)iced=(1+ηkCA)C(A)

(4)

式中：C(A)與C(A)iced分別表示任意的結冰前后飛機的性能、穩定性與控制參數或其導數；kCA表示結冰對飛機氣動參數的影響參數,與飛機自身的尺寸、飛行狀態或與飛機受結冰影響的敏感性相關；η為飛機結冰程度參數,只與氣象條件有關。該模型既考慮到了特定飛機的弦長、翼型、飛行速度等信息,又考慮到了大氣、結冰條件等信息,其作為分析結冰嚴重程度的一種估算模型,是一種相對而言比較科學合理的結冰后氣動參數計算方法,因此在國內外結冰后飛行動力學研究中,得到了廣泛的應用。但該模型在失速迎角前的迎角范圍內對結冰后氣動力描述并不準確,且無法描述飛機大迎角及過失速迎角階段的結冰后氣動特性[21]。由于涉及到全包線范圍內的品質分析,在飛行包線附近可能會出現非線性特性的影響。因此,需要建立非線性結冰影響模型。

傳統的線性結冰影響模型認為,結冰模型中的kCA值對于特定的氣動參數而言是一個常值,而參數η體現的是結冰嚴重程度,與迎角改變無關,因此不管迎角如何改變,該氣動參數對應的結冰后的值始終是等量地縮放。為此,對于飛機進入過失速區域之后,結冰前后氣動參數差異逐漸減小的現象,可以通過改變kCA在不同迎角區間上的取值來實現。文獻[21]具體分析了kCA在不同迎角情況下的取值。

傳統結冰影響模型中無法體現拐點的問題,這是由于式(4)本質上只是對氣動參數在原始干凈外形數值基礎上的縮放,也就無法體現一些拐點的移動。例如對于升力系數而言,結冰后最大升力系數的值降低了,但是其對應的失速迎角卻沒有改變。為此,可定義縮放因子ε,使得

α′=εαe,αe∈[α0,αp]

(5)

式中:αe為待縮放的迎角區間；α0為迎角插值表中的起始點；αp為大于拐點對應的迎角的設定值。對于升力系數CL而言,0<ε<1,通過計算可將失速迎角提前;對于俯仰力矩系數Cm而言,ε>1,通過計算可將俯仰力矩系數突變的迎角值推后。縮放因子的計算與結冰后失速迎角的估算結果直接相關。

2 飛機橫航向飛行品質評估方法

2.1 飛機橫航向擾動后運動參數表達式

按照美軍標(MIL-1797)中對于飛機橫航向運動模態中的分解,飛機的橫航向運動模態時域響應模型可以表示為

(6)

式中:等式右邊第一項為滾轉模態,滾轉角迅速衰減,滾轉模態特性可用一階時間常數τR[11]來表征:

φ=φR0e-t/τR

(7)

等式右邊第2項為螺旋模態,其時域響應模型的形式與滾轉模態相同,但螺旋模態與滾轉模態相比,其滾轉角變化較為緩慢,螺旋模態特性可用一階時間常數τS[11]來表示:

φ=φS0e-t/τS

(8)

(9)

在選擇根據哪一個橫航向的狀態參數進行時域等效擬配時,文獻[6]提出采用滾轉角對橫航向模態參數進行等效擬配,在實際工程應用中發現,由于滾轉角對外界擾動以及駕駛員輕微的動桿敏感性較強,在某些情況下擬配出來的結果誤差較大,甚至得出不合理的結果。而橫航向狀態參數中,飛機偏航角速率r對外界擾動的敏感性較低,可用來對飛機橫航向參數進行時域等效系統擬配。

2.2 基于數據特征的自適應橫航向擬配初值計算流程

圖2 傳統時域等效擬配方案Fig.2 Traditional time-domain equivalent allocation scheme

從典型的橫航向時域響應曲線可以看出,飛機的橫航向模態已經體現在曲線的特征當中,如曲線總體的發散收斂體現的是螺旋模態特性,曲線中波峰波谷極值點的位置及幅值體現的是荷蘭滾模態的周期及阻尼比。因此,如果能從飛機橫航向響應曲線中提取出飛機模態特征,來獲取初步的橫航向模態特征參數值,即可用于等效擬配時初值的選取。根據上述思想,本文提出了根據飛機時域響應曲線的特征,分步估算出滿足最小二乘的模態特征參數,達到根據試驗數據自動計算擬配初值的目的。

上述基于橫航向時域響應特征自適應計算橫航向等效系統擬配初值方法的具體工作流程如圖3所示。

圖3 基于數據特征的自適應橫航向擬配初值計算流程Fig.3 Flow chart of adaptively calculating the initial value of the lateral and directional equivalent system fitting based on data characteristics

在截取了方向舵倍脈沖試驗數據后,數據處理模塊均勻地在滾轉角響應曲線上選取n組待擬配的數據。以這n組數據為輸入,計算出滿足最小二乘的螺旋模態特征參數初值(如圖4(a)所示)。同時,找出響應曲線上所有的極值點,根據這些極值點對應的時間坐標估算出荷蘭滾模態的周期,進而計算荷蘭滾模態自振頻率ωn d。而后,再計算出每一個極值點到螺旋模態曲線的最小距離Δdi,根據Δdi計算出滿足最小二乘的荷蘭滾模態阻尼比,即ωn dζd,進而得到阻尼系數ζd，如圖4(b)所示。而荷蘭滾模態與螺旋模態幅值參數rR0和rS0的初值可按前兩個極值幅值的平均值進行估算。

圖4 橫航向關鍵模態參數初值計算原理圖Fig.4 Schematic diagram of initial value calculation of lateral and directional key modal parameters

飛機滾轉模態的時間常數通常為1～3 s,說明滾轉模態在橫航向響應中一般會很快衰減,因此在進行時域等效擬配時,可以從方向舵倍脈沖結束后的4～5 s以后開始截取數據,則可以認為此時只保留了荷蘭滾模態與螺旋模態項,滾轉模態項的影響忽略不計。

從兩種方式的流程圖可以看出,兩種方法最大的區別在于初值的選取和計算,傳統方法采用預設初值法或根據遺傳算法估算初值,往往只能用于某個特定狀態點下進行品質分析,均不能很好地解決飛機在所有情形下自動進行等效系統擬配并最終得到最優解的問題,本文擬采用的方案為根據數據估算初值,這對于適應各種飛行初始狀態具有不可比擬的優勢。

2.3 算例驗證

為驗證文中提出的通過基于自適應擬配初值計算的時域等效系統擬配的準確性,通過構造一典型的III型運輸類飛機橫航向傳遞函數模型,再利用文中提出的方法對其時域響應進行模態特性分析,來進行對比研究。典型的飛機橫航向低階等效系統傳遞函數形式為

(10)

根據GJB-18586中對飛機橫航向模態參數品質等級的劃分,設定15組模態參數如表1所示,其擬配結果如表2所示。

表1 參數設定Table 1 Parameters setting

表2 擬配結果Table 2 Fitted results

根據第1～3組,第5～7組和第9～11組的擬配結果可以發現,當螺旋模態參數發生變化時,失配度均處于10-5數量級,并且各模態參數計算結果與設定值基本一致,說明螺旋模態參數的變化對擬配結果影響不大。

根據第1、4、7組,第2、5、8組和第3、6、9組的擬配結果可以發現,當荷蘭滾模態參數發生變化時,失配度均處于10-5數量級,并且各模態參數計算結果與設定值基本一致,說明荷蘭滾模態參數的變化對擬配結果影響不大。

根據第10～15組的擬配結果可以發現,隨著滾轉模態時間常數的增加,失配度也迅速增大,當滾轉模態時間常數大于3 s后,失配度超過了0.001。說明滾轉模態參數增加會使得擬配效果變差,這是因為通常飛機滾轉模態的時間常數通常為1～3 s,而本文通過上文描述的時域等效擬配算法設計,使得飛機滾轉模態的時間常數小于3 s時的滾轉模態參數變化的影響可忽略不計。

考慮到飛機方向舵倍脈沖周期的不同,會對擬配結果產生一定影響。控制各模態參數的設定值不變,依次增大脈沖周期,擬配結果如表3所示。

表3 方向舵倍脈沖信號周期遞增時的擬配結果Table 3 Fitted results of rudder double pulse signal with increasing period

根據表3中的擬配結果可以發現當脈沖周期小于或等于4 s時,其變化對擬配結果總體精度的影響不大,失配度始終處于10-5～10-3,并且各模態參數計算結果與設定值基本一致。當脈沖周期大于4 s時,失配度迅速增大,這是因為本文所設計的飛機橫航向模態特性分析方法是基于飛機橫航向小擾動響應模型,當脈沖周期過大時,脈沖輸入可近似為階躍輸入,此時飛機橫航向小擾動響應模型不再適用。

通過以上驗證設計,可以得出以下結論:在某型飛機上應用本文所提出的基于等效系統擬配的飛機橫航向模態特性分析方法時,需要限制該機型的設計滾轉模態時間常數小于3 s,同時,倍脈沖輸入的周期小于4 s,就可以保證較高的精度。

根據文獻[11]中所述方法,本文通過對某次地面鐵鳥臺試飛數據進行分析來驗證文中提出的基于自適應擬配初值計算的時域等效系統擬配方法的準確性。擬配曲線與實際飛行曲線如圖5所示,失配度為0.000 148 9,滿足準確度要求。

圖5 計算結果對比Fig.5 Comparison of calculation result

綜上所述,控制脈沖周期小于或等于4 s且從方向舵倍脈沖結束后的4～5 s以后開始截取數據,文中提出的基于自適應擬配初值計算的時域等效系統擬配方法滿足準確度要求,其計算結果可以作為飛行品質評估的依據。

3 全飛行包線范圍內飛機橫航向模態特性分析及飛行品質評估

某型飛機的全飛行包線范圍內橫航向模態特性分析飛行品質評估具體工作流程如圖6所示。本文在已知背景飛機干凈構型的飛行包線的基礎上,以飛行包線左邊界點為采樣起點,飛行包線右邊界點為采樣終點,采樣高度從1 000 m至12 000 m,高度采樣間隔為1 000 m,在每一高度下均勻取15個采樣點。進行飛行動力學仿真之前需要進行配平,如果配平結果中發動機推力超過背景飛機實際可用最大推力或配平結果顯示飛機不能保持定直平飛狀態,則可以判斷該采樣點不在飛行包線內。通過飛行動力學仿真得到采樣點對于方向舵倍脈沖操縱的響應數據,便可通過上述基于數據特征的自適應橫航向等效系統擬配獲取背景飛機橫航向模態特性,進而進行背景飛機飛行品質評估。

圖6 全飛行包線范圍模態特性分析飛行品質評估流程圖Fig.6 Flow chart of mode characteristics analysis flight quality assessment for full flight envelope

3.1 并行數值仿真計算

在全包線范圍內進行飛機橫航向模態特性計算時,飛行包線即為計算區域,而飛行包線是由各個狀態點組成,可將計算區域離散成各個相互獨立計算單元進行并行計算,通過本文提出的基于自適應初值計算的飛機橫航向模態特性計算方法可得到該飛行狀態的各模態參數及失配度,進而根據飛機橫航向飛行品質評估標準可得到全包線范圍內的飛行品質評估結果。

3.2 某型飛機的全飛行包線范圍的橫航向模態特性分析及飛行品質評估

本文在建立典型的飛機本體模型和電傳操縱系統模型的基礎上,得到了某型飛機在干凈構型條件下其飛行包線內各個狀態點對于方向舵倍脈沖操縱的響應數據。

本文主要研究機翼對稱結冰對飛機的飛行性能的影響,在建立了某型飛機結冰線性氣動力模型的基礎上,得到了某型飛機分別在結冰嚴重程度參數η分別為0.1、0.15、0.2時,其飛行包線內各個狀態點對于方向舵倍脈沖操縱的響應數據。

運用上述提出的基于參數響應特征的時域等效系統擬配法,獲取了飛機在整個包線范圍內橫航向模態參數特性,每一個狀態點上失配度的計算結果如圖7所示。其中,綠色點代表該狀態點的失配度小于10-5,藍色點代表該狀態點的失配度的范圍在10-5～10-3,紅色點代表該狀態點的失配度的范圍在10-3～10-2。

圖7 各狀態點的失配度Fig.7 Degree of mismatch at each state point

由計算結果可知,在背景飛機干凈構型條件下絕大多數點的失配度都在10-5～10-3之間,個別的點小于10-5,只有少數點的失配度在10-3～10-2之間,說明所提出的等效擬配方法在全包線范圍內具有較好的適應性;在不同結冰嚴重程度條件下絕大多數點的失配度都在10-5～10-3之間,小部分點的失配度在10-3～10-2之間,少數點的失配度在小于10-5,飛機結冰后失配度的變化不大,說明所提出的等效擬配方法在結冰條件下全包線范圍內同樣具有較好的適應性。

背景飛機在干凈構型條件下,飛行包線內一共包含314個狀態點,當結冰程度參數達到0.1時,飛行包線所包含的狀態點減少為296個,背景飛機的飛行包線沒有明顯變化；當結冰程度參數達到0.15時,飛行包線所包含的狀態點減少為263個,飛行包線范圍縮小了6.69%；當結冰程度參數達到0.2時,飛行包線所包含的狀態點減少為211個,飛行包線范圍縮小了32.80%,計算結果驗證了本文所建立的結冰線性氣動力模型的準確性。

本文采用荷蘭滾模態阻尼作為飛行品質評估的依據,評估標準如表4所示。

表4 橫航向飛行品質評估標準Table 4 Evaluation standard for lateral and directional flight quality

以荷蘭滾阻尼為評估標準,某型飛機分別在干凈構型及不同結冰嚴重程度條件下的飛行品質評估結果如圖8所示,其中綠色點代表一級品質,藍色點代表二級品質,紅色點代表三級品質。

圖8 飛行品質評估結果Fig.8 Flight quality evaluation results

上述結果表明,在干凈構型條件下,以荷蘭滾模態阻尼比為飛行品質評估標準時,背景飛機的飛行品質大致以高度8 000 m為界限,低于8 000 m為一級飛行品質,高于8 000 m則基本為二級飛行品質,共有224個狀態點為一級飛行品質,90個狀態點為二級飛行品質,一級飛行品質區域占包線范圍的71.34%。當結冰程度參數達到0.1時,一級飛行品質的狀態點減少為71個,二級飛行品質的狀態點增加為225個,一級飛行品質區域占包線范圍的23.99%,一級飛行品質與二級飛行品質的分界線下降至3 000 m左右;當結冰程度參數達到0.15時,一級飛行品質的狀態點減少為30個,二級飛行品質的狀態點增加為233個,一級飛行品質區域占包線范圍的12.86%,一級飛行品質與二級飛行品質的分界線下降至2 000 m左右;當結冰程度參數達到0.2時,無一級飛行品質的狀態點,二級飛行品質的狀態點增加為211個。上述結果表明,以荷蘭滾阻尼為評估標準時,結冰嚴重影響某型飛機的飛行品質,隨著結冰嚴重程度的增加,某型飛機的飛行品質在高度由上至下逐漸由一級飛行品質轉變為二級飛行品質,直至所有的狀態點均變為二級品質。

限于篇幅,取高度為1 500 m的各狀態點的橫航向運動參數擬配結果分別如表5～表8所示。

表5 干凈構型Table 5 Clean configuration

表6 結冰嚴重程度參數η=0.1Table 6 Icing severity parameter η=0.1

續表6Continued Table 6

表7 結冰嚴重程度參數η=0.15Table 7 Icing severity parameter η=0.15

表8 結冰嚴重程度參數η=0.2Table 8 Icing severity parameter η=0.2

需要指出的是,飛機結冰一般發生在7 000 m以下,飛行速度較低的情形。文中對全包線范圍內結冰后動力學特性變化都進行了分析,主要是從理論上來驗證文中提出的算法的有效性,量化結冰后飛行包線的變化范圍。在實際運用中,分析具體的飛機結冰后飛行品質變化,需要結合具體的易結冰飛行范圍來進行分析。

4 結 論

本文首先建立了典型的人-機-環系統模型,通過數值仿真得到飛機在飛行包線內所有狀態點對于方向舵倍脈沖的響應數據。其次，開發了自動數據預處理模塊，改善了傳統算法因擬配初值設定無法涵蓋所有情形的問題，并通過具體算例驗證了改進后程序的有效性。在動力學仿真數據的基礎上,運用時域等效系統擬配法,分別獲取了飛機干凈構型及不同結冰嚴重程度下整個包線范圍內橫航向模態參數特性。每一個狀態點上失配度的計算結果表明,所提出的等效擬配方法在全包線范圍內具有較好的適應性，為快速分析飛行器縱橫向典型模態特性提供了新的思路,具有較好的工程應用前景。