基于微分流形理論的結冰條件下飛機魯棒動力學邊界研究

李穎暉, 桂業偉, 左仁偉, 董澤洪, 張久星

(1.空軍工程大學航空工程學院，西安，710038；2.中國空氣動力研究與發展中心空氣動力學國家重點實驗室，四川綿陽，621000；3.93756部隊，天津，300131)

飛機結冰是過冷水滴(溫度低于冰點的液態水)撞擊機體表面并發生相變的一種復雜現象[1-3]。結冰能夠破壞飛機的氣動特性，降低飛機的操穩性能，對飛行安全造成嚴重威脅[4-6]。盡管國內外學界針對飛機結冰的研究已經持續了數十年，期間發展出多種防/除冰技術[7-8]，但由于國內外航班數量顯著增長，結冰引發的飛行事故屢屢發生：據美國運輸安全委員會(NTSB)統計，1978—2005年間結冰共造成645起飛行事故[9]，同時期美國航空航天局(NASA)報道了299起結冰引發的飛行事故。隨后在2006—2010年間，NTSB和NASA報道了228起結冰相關的飛行事故[10]。

由于飛機承載能力和負荷的限制，部分型號尚未裝備防除冰系統。此外，由于最大載荷的限制，部分裝備防除冰系統的飛機也不能將機翼表面的積冰完全除盡，且積冰融成的液態水在機翼后緣易再次凝結。為實現飛機帶冰安全飛行，帶冰飛行已成為當前亟待解決的課題。要解決飛機帶冰飛行的安全問題，首先需要掌握結冰對飛機飛行性能的影響機理。鑒于此，國外學者主要采取飛行試驗[11]、風洞試驗[12]和數值模擬[13]3種方式進行研究，并且已經在飛機結冰試驗設備、研究軟件等方面取得了顯著進展，形成了大量試驗數據和經驗規律。2000年，NASA在關于飛機結冰的報告中詳細分析了結冰對飛機氣動特性的影響[14]，2001年，Bragg團隊提出了能夠反映不同結冰程度下氣動參數變化的結冰因子影響模型[15-16]。近年來，國內學者提出了飛機結冰的多重安全邊界的概念，并指出飛機結冰防護系統應依據結冰致災鏈路中各環節的邊界對應的防護需求進行設計[17]。在此基礎上，學者們進一步提出一系列結冰條件下飛機安全邊界保護方法，主要涵蓋結冰條件下飛機風險評估、參數辨識和控制律設計等問題[18-25]。然而，上述研究均依賴于精準的飛機結冰模型，未考慮到帶冰飛機飛行環境的多變性和復雜性，由于實際飛行中很難獲取精確的冰型信息和結冰程度，因此亟需開展針對動力學邊界的魯棒性研究。

1 流形理論及結冰條件下飛機動力學模型

1.1 流形理論

流形理論是一種計算非線性系統穩定域的有效方法,其利用積分的方式獲取邊界上不穩定平衡點的穩定流形并將其作為穩定邊界[20]。計算二維流形的方法主要包括連續邊值問題法、測地圓法、偏微分方程法和擴展軌線法[26]。本文采用測地圓法刻畫系統的流形，具體流程見圖1。

圖1 測地圓法刻畫二維流形

1.2 縱向動力學模型

本文采用如下的飛機縱向通道動力學模型：

(4)

式中：

(5)

GTM飛機采用放寬穩定性技術，為提高穩定性，設計增穩控制律見式(6)，增穩補償控制器結構框圖如圖2所示。

圖2 增穩補償控制器結構框圖

δe=kαΔα+kqΔq

(6)

式中：δe為升降舵偏角；α0、q0分別是配平點處的迎角和俯仰角速度；Δα=α-α0和Δq=q-q0為狀態變化量；控制系數取kα=1、kq=0.01。

以H=4 000 m、V=162.25 m/s為飛機狀態，證明以上方法的準確性。本文研究的配平點處飛行狀態信息具體如下：

(7)

式中：δth為推力系數。

1.3 結冰模型

本文采用Bragg團隊于2000年提出的飛機結冰模型[15-16]，該模型初步解決了結冰氣動數據的獲取問題，具體表達式如下：

CA,iced=(1+ηfice)CA

(8)

式中：CA和CA,iced分別是結冰前、后飛機氣動導數值；fice是結冰系數，反映CA對結冰的敏感性，對于給定的飛機為常值，不同氣動系數對應的fice不同；0≤η<1是結冰因子，結冰程度越嚴重，η越大。結冰因子體現在飛機的氣動參數模型中，對于縱向通道，主要影響飛機的升、阻力系數和俯仰力矩系數，結冰程度越嚴重，結冰因子的值越大，升力系數和俯仰力矩系數下降越多，阻力系數增加越多。

圖3 結冰前后升力系數曲線變化示意圖

2 動力學的確定及準確性驗證

2.1 動力學的確定

為了詳細說明基于微分流形理論計算飛機動力學邊界的過程，首先計算該狀態下所有的穩定平衡點(Stable Equilibrium Point，SEP)和不穩定平衡點(Unstable Equilibrium Point，UEP)，如表1所示。

飛機的動力學邊界時由SEP穩定邊界上所有UEPs的流形組合而成，在求出所有的平衡點后需要驗證這些平衡點是否處于穩定邊界上。簡單的方法是在UEP處施加一個指向SEP的小擾動，這樣以來，UEP點處不穩定的平衡狀態被破壞，如果狀態最終能夠收斂至SEP，則該UEP處于穩定邊界上，反之不然[29]。將表1中的所有Ⅰ型UEP代入上述流程，結果如圖4所示。

表1 所有平衡點

圖4 確定穩定邊界上的UEPs

UEP1和UEP2在受到擾動后收斂至平衡點SEP，而UEP3受擾動后呈發散趨勢，因此UEP1和UEP2在穩定邊界上，UEP3不在SEP的穩定邊界上。由微分流形理論確定系統穩定域的方法具體如下：首先要確定SEP邊界上所有的UEP，然后計算每個UEP的穩定流形，將這些流形組合起來即圍成非線性系統的穩定域。本例中，UEP1和UEP2的穩定流形包圍的區域即為GTM飛機在所給工作狀態的三維縱向穩定域，如圖5所示。

圖5 微分流形法得到的穩定域

圖5中紅點表示SEP，藍點表示UEP，紫色曲面表示動力學邊界，紫色曲面圍成的區域即為穩定域?？梢婏w機的穩定域并不是規則的幾何形狀，在計算維度超過二維后難于使用解析方法(如Lyapunov能量函數法)得到準確的穩定域。

2.2 動力學邊界的準確性驗證

由于飛機的縱向三維穩定域不是規則的幾何形狀且難于使用解析方法進行驗證，本文采用動力學仿真法進行準確性驗證。Monte Carlo法可在選取范圍內以一定的間隔取不同初始狀態，力圖遍歷范圍內所有狀態點，通過大量動力學仿真，最終收斂至SEP的即為動力學邊界內的狀態點(簡稱內點)，所有內點的集合即為穩定域內的狀態點，從而可以得到飛機的動力學邊界。因Monte Carlo法的理論具備基礎性和廣泛性，通常被用作復雜理論的驗證。本文為了驗證基于微分流形理論所得穩定域的準確性，采用Monte Carlo法對其進行驗證，網格的劃分以及范圍的選取如表2所示。

表2 基于Monte Carlo法的狀態范圍及網格劃分

值得注意的是，Monte Carlo法雖然理論簡單易懂，但其計算效率低，計算耗時長，而飛機動力學特性復雜，且動力學因素間耦合性強，加大了計算的復雜度，為了說明微分流形法在計算效率上的優越性，利用微分流形法和Monte Carlo法在Intel(R)Core(TM)i7-4690 CPU、主頻3.6 GHz、8 GB內存的臺式機上進行相同工況的動力學邊界仿真計算，使用微分流形法用時1.34 min，而使用Monte Carlo法用時124 min，可說明微分流形法具有較高的計算效率。基于Monte Carlo法得到的動力學邊界驗證結果如圖6所示。

圖6 微分流形法和Monte Carlo法的穩定域對比

圖6中黑點構成的區域代表利用Monte Carlo法得到的穩定域，紫色曲面圍成的區域代表利用微分流形法得到的穩定域。從圖中可以看出，2種方法得到的穩定域吻合，驗證了微分流形法所得穩定域的準確性。同時其具備較高的計算效率，且能處理動力學因素耦合的復雜情形，因此在計算飛機動力學邊界時具有獨特優勢。

3 飛機的魯棒動力學邊界

3.1 結冰條件下飛機動力學邊界

結冰會導致在相同迎角下飛機的升力下降，受到的阻力增加，舵效降低。為了維持飛機飛行所必須的升力，結冰條件下飛機必須通過增大迎角來獲得足夠的升力，而迎角的增加勢必造成阻力的增加，此時飛機容易進入失速狀態。為保證飛行安全，必須在飛機陷入失速前及時的糾正飛機的狀態，因此構建一個準確的結冰動力學邊界對于保障結冰條件下飛機的飛行安全具有關鍵作用。

以H=4 000 m、V=162.25 m/s、α=0.087 7 rad、θ=0.087 7 rad、q=0 rad/s為研究狀態，考慮結冰程度η=0.2和η=0.3的情形，結冰條件下飛機的動力學邊界如圖7所示。

圖7中綠色曲面為干凈飛機，反映未結冰情況；藍色曲面為結冰因子η=0.2時的動力學邊界，反映輕度結冰情況；橘色曲面為結冰因子η=0.3時的動力學邊界，反映中度結冰情況。綠點、藍點和橘點分別表示干凈飛機、結冰因子η=0.2和η=0.3的配平點?？梢婋S著結冰程度的加劇，由于氣動性能下降，飛機需要更大的配平迎角以提供足夠的升力。同時，由圖7可知，飛機的動力學邊界隨著結冰程度的加劇不斷收縮。

3.2 動力學邊界的魯棒性研究

在飛機飛行過程中，由于飛行環境持續變化，飛機的氣動參數存在諸多不確定性。且當飛機結冰條件下機體表面的冰在不斷累積，以目前的技術手段難以獲取準確的冰型信息和結冰程度，因此，對于結冰條件下飛機動力學邊界的魯棒性研究相當重要且迫切。

為了使求得的動力學邊界具有魯棒性，定義魯棒氣動參數如下：

CA,robust=(1±Δ)CA

(9)

式中：Δ是飛行環境不確定性引發的攝動系數；CA,robust為考慮攝動后的氣動導數值。

圖8是考慮魯棒性的結冰條件下飛機氣動參數變化范圍，藍線和紅線分別是Δ=0和Δ=0.5時的升力系數CL、阻力系數CD及俯仰力矩系數Cm；綠色范圍為結冰條件下飛機氣動參數的容許變化范圍。這樣，當飛機在惡劣環境下飛行時，只需保證氣動參數在容許變化范圍內，飛機就可以在具有魯棒性的動力學邊界內根據任務需要進行機動，所以利用圖8所示的氣動參數計算得到的動力學邊界具有很強的魯棒性。

圖8 考慮魯棒性的氣動參數范圍

下面舉例說明魯棒性動力學邊界在保障結冰條件下飛機安全方面的優勢，以“結冰探測系統辨識結冰程度η=0.2，但由于機械誤差的存在，實際的結冰程度更加嚴重”為背景。飛機超出動力學邊界時將開啟邊界保護系統使飛行狀態收斂至平飛狀態。若以η=0.2時的動力學邊界作為邊界保護系統開啟的依據，飛機做機動至α=0.4 rad，θ=0.05 rad，q=-1 rad/s仍被認為是安全狀態，但此時飛機已經超出了實際的動力學邊界，這時飛機的關鍵姿態參數很快會發散，飛行狀態變化軌線如圖9～10所示。

圖9 結冰條件下飛機飛行狀態

圖9中綠色曲面表示Δ=0時的動力學邊界，紅色曲面表示Δ=0.5時的魯棒動力學邊界，藍色曲面表示η=0.2時的動力學邊界，紅色星號表示配平點，黑點表示飛機機動后的初始點，黑色曲線表示由初始點出發的飛機狀態軌線。圖10為以η=0.2對應的動力學邊界作為安全邊界時的結冰條件下飛機姿態角變化曲線，可見姿態角很快發散，飛機變得失穩失控。為了避免事故的發生，必須在飛機超出實際安全邊界前采取措施，使用本文提出的魯棒性動力學邊界可有效解決這一問題，結冰條件下飛機魯棒動力學邊界保護系統框圖如圖11所示。

圖10 結冰條件下飛機姿態角

圖11 結冰條件下飛機魯棒動力學邊界保護系統

圖11中紅色虛線框內表示結冰條件下飛機魯棒動力學邊界保護系統，當結冰檢測系統或氣象部門探測到飛機遭遇結冰惡劣飛行環境時，為了保證飛行安全，飛機的關鍵參數須保持在結冰魯棒動力學邊界內，當飛行狀態到達邊界時，舵偏控制決策器切換至限幅模式，阻止飛機超出動力學邊界。因為采用的動力學邊界具有魯棒性，所以飛機可在超出實際動力學邊界之前回穩，從而達到保護結冰條件下飛機飛行安全的目的。

圖9中飛機姿態角發散是因為其超出了實際的動力學邊界，而迎角、俯仰角速度屬于快變量，一旦超出安全閾值很快就會發散，留下的反應時間很短。采用圖11所示的結冰條件下飛機魯棒動力學邊界保護系統，在結冰程度η=0.2的情形下做機動，當飛行狀態到達邊界處，邊界保護系統中的限幅器作用使俯仰角速度不能進一步減小，從而使飛行狀態保持在安全閾值內，在完成機動后飛機可快速回穩。

使用結冰條件下飛機魯棒動力學邊界保護系統后的飛機姿態變化曲線見圖12。圖12中實線表示迎角變化曲線，點線表示俯仰角變化曲線，虛線表示俯仰角速度變化曲線?？梢姴捎帽疚奶岢龅慕Y冰條件下飛機魯棒動力學邊界保護系統，在機動結束后飛機的狀態可自發收斂至平衡點。所提方法可保障飛機在結冰情形下的飛行安全，允許飛機在一定狀態緊集內做機動，通過限制升降舵的操縱閾值實現限制飛行狀態的效果。

圖12 結冰條件下飛機魯棒動力學邊界保護方法下的姿態角

4 結語

本文提出了一種結冰條件下飛機魯棒動力學邊界保護系統設計方法，當飛機遭遇結冰不利飛行環境時，不需獲取準確的冰型信息和結冰嚴重程度，考慮當前結冰檢測技術手段，本文提出的方法具有更加廣泛意義的工程應用意義。在飛機超出動力學邊界之前邊界保護系統作用使飛行狀態保持在安全閾值范圍內，從而保證飛機安全可控。