基于復雜動力學仿真的結冰情形下飛行安全窗構建方法

裴彬彬， 徐浩軍， 薛源,*， 李哲， 劉東亮

1.空軍工程大學 航空航天工程學院， 西安 710038 2.北京航空工程技術研究中心， 北京 100076

基于復雜動力學仿真的結冰情形下飛行安全窗構建方法

裴彬彬1， 徐浩軍1， 薛源1,*， 李哲1， 劉東亮2

1.空軍工程大學 航空航天工程學院， 西安 710038 2.北京航空工程技術研究中心， 北京 100076

目前針對結冰情形下增強駕駛員情景感知的研究比較有限，現(xiàn)有的手段一般為通過評估部分飛行安全關鍵參數(shù)是否超出其極限值來對風險事件是否發(fā)生進行預測。建立了駕駛員操縱-飛機本體-積冰影響的動力學模型，通過對單個飛行情形預測時間段內(nèi)飛行參數(shù)風險度的疊加得到該情形下的飛行安全譜，并在此基礎上得到該情形下的風險值?；诮⒌牟⑿酗w行仿真平臺，獲取飛機在整個操縱范圍內(nèi)的風險拓撲圖，即安全窗。分析了飛機在對稱結冰情形和非對稱結冰情形下飛行安全窗的變化，并對結冰的致災機理進行了分析。仿真結果表明結冰導致安全飛行范圍縮減，對于非對稱結冰還會出現(xiàn)安全窗不對稱的現(xiàn)象。安全譜的提出可以為事故的演化分析提供一種全面直觀的分析方法，安全窗的構建可為飛機遭遇各種不利情形下的駕駛員操縱提供指示，也可為飛機設計人員優(yōu)化飛機性能提供一定的參考。

飛機結冰； 情景感知； 計算飛行動力學； 非對稱結冰； 并行飛行仿真； 安全窗

現(xiàn)代事故致因理論認為，飛行事故的發(fā)生往往是人-機-環(huán)復雜系統(tǒng)內(nèi)部因素之間耦合作用的結果[1]。駕駛員在事故的演化過程中扮演了重要的角色，駕駛員在各種情形下的操縱正確與否直接關系到飛行安全。波音公司對全球商用噴氣式飛機事故統(tǒng)計結果表明，超過50%的飛行事故主要由機組人員造成[2]，空客公司曾指出大約85%的飛行事件報告中均提到了駕駛員情景感知(Situational Awareness)能力的喪失[3]。提高駕駛員的情景感知能力對于飛行安全至關重要，這也是當今國際航空領域研究的熱點問題之一[4-8]。

在眾多的研究中，針對結冰情形下增強駕駛員情景感知的研究并不多見。Bragg[9]和Deters[10]等開創(chuàng)性地提出了智能結冰系統(tǒng)(Smart Icing System，SIS)的概念，該系統(tǒng)能夠?qū)崟r地感知積冰的存在及其影響，并計算出結冰狀態(tài)下飛機的飛行安全邊界提供給駕駛員及飛控系統(tǒng)。Gingras等[11-12]聯(lián)合開發(fā)了積冰污染邊界保護(Icing Contamination Envelope Protection，ICEPro)系統(tǒng)，該系統(tǒng)為駕駛員提供積冰信息及飛機的邊界信息，具備提升駕駛員在結冰情形下情景感知的能力，其在地面模擬器試飛的結果得到了駕駛員的肯定[13]。

國內(nèi)針對飛機結冰后飛行動力學特性改變對飛行安全影響的研究起步較晚。徐忠達等[14]就冰形對氣動參數(shù)及操穩(wěn)特性的影響進行了研究；Dong和Ai[15]研究了結冰后的參數(shù)辨識；應思斌和艾劍良[16]主要就結冰條件下的容冰控制進行研究。劉東亮等[17]提出了積冰條件下基于復雜系統(tǒng)動力學仿真的風險評估方法。王明豐等[18]就結冰對飛行動力學的影響進行了分析，張智勇[19]主要對智能防冰系統(tǒng)中的包線保護控制律設計問題進行了研究。由于研究手段等方面的限制，在致災機理與防護方面，國內(nèi)尚無系統(tǒng)性的研究。

當前結冰邊界保護系統(tǒng)為駕駛員提供的都是實時的飛行狀態(tài)及其邊界信息，如迎角和舵面操縱量及其限制值等。NASA Langley研究中心近期(2013年)的一項調(diào)查結果表明[20]，駕駛員更傾向于得到可用滾轉角和俯仰角等指示信息，特別是在安全邊界改變的情形下，同時，提供這些信息并不會帶來駕駛員負荷的顯著增加。另一方面，可用滾轉角和俯仰角等信息相對于實時的飛行狀態(tài)參數(shù)及其極限值而言，更具預測性。而這正是當前飛機結冰研究所缺乏的。因此，研究結冰情形下飛機俯仰角和滾轉角的可用安全范圍對于提高駕駛員情景感知能力、保障飛行安全具有重要的作用。

本文提出了基于計算飛行動力學的結冰條件下飛行安全窗的構建方法，將飛機在機動空間中可能飛行軌跡的飛行風險以色彩化的形式呈現(xiàn)出來，揭示人機閉環(huán)系統(tǒng)在結冰情形下的安全邊界以及不同結冰情形下的致災機理，為駕駛員的安全操縱提供直觀的指示告警。

1 模型的建立

對于結冰情形下飛行安全區(qū)域的確定要涉及臨界飛行狀態(tài)，往往伴隨著強耦合、非線性等特點，因此采用6自由度全量非線性運動方程更能準確地描述飛機的運動規(guī)律。分別對飛機本體動力學模型、駕駛員操縱策略模型、積冰影響模型進行了研究，建立起不同結冰情形下的人-機-環(huán)閉環(huán)系統(tǒng)模型，為進行動力學仿真打下基礎。

1.1 飛機動力學模型

1.1.1 飛機本體模型

飛機本體非線性動力學模型可用式(1)所示的向量形式表示[21-24]

(1)

式中：x為狀態(tài)向量;u為控制向量；且

x=[Vαβq0q1q2q3pqr

xgygzg]T

(2)

其中:V、α和β分別為飛行速度、飛機迎角和側滑角；xg、yg和zg為飛機在地面坐標系下的位置。

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

(3)

其中：δth為飛機油門偏度；δe、δa和δr分別升降舵、副翼以及方向舵偏轉的大小。

為避免奇點的產(chǎn)生，采用四元數(shù)法來建立飛機動力學模型，得到向量f的表達式為

(4)

(5)

(6)

(7)

此外，

(8)

(9)

飛機的姿態(tài)角：滾轉角φ、俯仰角θ和偏航角ψ可以通過四元數(shù)變換求得。

1.1.2 舵機模型

舵機模型可簡化為由一階慣性環(huán)節(jié)、速率限制器和位置限制器構成[25]，如圖1所示。圖中：T為一階慣性環(huán)節(jié)的時間常數(shù)；1/(Ts+1)為舵機的傳遞函數(shù)。

圖1 舵機動力學模型示意圖
Fig.1 Schematic diagram of dynamics model for actuator

舵機速率限制器以及位置限制器的具體參數(shù)可參照文獻[26]中的有關參數(shù)進行設置。

1.2 駕駛員模型

典型的駕駛員模型結構如圖2所示[27]，輸入信號是簡單的誤差信號，即基準的輸入量與人機閉環(huán)系統(tǒng)實際響應量之間的偏差，虛線框中的4個模塊構成了駕駛員主要行為特性的建模，各模塊的說明如圖2所示。

1) 含觀測噪聲的顯示環(huán)節(jié)

駕駛員通過儀表盤等接收飛機當前飛行狀態(tài)時引入的誤差。

2) 延遲環(huán)節(jié)

駕駛員在讀取信號、信號傳輸至大腦以及大腦做出決策的過程存在著一個不可控的最小反應延遲時間，常用e-τs表示，τ的范圍一般為0.06～0.30 s。

3) 補償操縱環(huán)節(jié)

補償操縱環(huán)節(jié)相當于是反饋控制系統(tǒng)中的控制器，構成了駕駛員模型的核心。在指定飛機目標俯仰角、滾轉角的情況下，根據(jù)文獻[28-29]中采用的駕駛員模型，本文在各通道上建立了以下駕駛員補償操縱策略模型。

俯仰角控制模型為

(10)

式中：Ixx、Iyy、Izz和Ixz為飛機的轉動慣量；Kθ1和Kθ2為常數(shù)值；其他中間變量的計算參照文獻[28-29]中所示。滾轉角控制模型為

(11)

式中：φc為指令滾轉角;Kpa、KDa和KIa為常數(shù)值，但在不同的φc區(qū)間段取值不同。

發(fā)動機油門控制受飛行速度影響，駕駛員操縱使得飛行速度盡量保持不變。

補償操縱環(huán)節(jié)中駕駛員控制模型的參變量直接影響最終的安全窗范圍，本文在進行參數(shù)的選取時，以盡量能夠使飛機按照給定操縱指令飛行為準繩。

4) 神經(jīng)肌肉動力學模型

當肌肉接受大腦發(fā)出的指令時，由于自身的慣性、黏性和收縮而呈現(xiàn)神經(jīng)肌肉的時延特性，用一階慣性環(huán)節(jié)1/(1+TNs)表示；此外駕駛員根據(jù)經(jīng)驗可采取提前的操縱量來補償這種滯后而呈現(xiàn)導前特性，可用一階微分環(huán)節(jié)1+TLs來表示[30]。TN和TL的值均約為0.1～0.2 s。另外，駕駛員在操縱過程中會出現(xiàn)自發(fā)性的手抖現(xiàn)象，同樣會引入噪聲信號。

圖2 駕駛員模型結構示意圖
Fig.2 Schematic diagram of structure of pilot model

觀測噪聲及操縱過程中的噪聲信號對動力學仿真的結果影響并不大，為適當?shù)睾喕嬎?，文中予以忽略?/p>

1.3 結冰影響模型

飛行中根據(jù)氣象條件來對冰形進行預測并估算對飛機空氣動力特性的影響是對飛機飛行動力學特性改變進行計算的基礎，也是最終實現(xiàn)飛行風險操縱區(qū)域劃分的基礎，國內(nèi)外相關的研究也比較多。本文的研究重點在于在已經(jīng)預測出冰形對飛機空氣動力特性影響的前提下，計算飛機安全操縱范圍，從而為駕駛員操縱提供有價值的信息。

1.3.1 對稱結冰影響模型

利用6自由度(6-DOF)運動方程對飛機在積冰后的動力學特性進行研究，積冰前后，運動方程的形式并沒有變，只是由于積冰影響了飛機的氣動特性，使得方程中的力和力矩發(fā)生了變化。

文中采用Bragg等[31]提出的積冰影響模型對結冰后氣動力進行建模。按照其理論，積冰前后氣動參數(shù)的關系為

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

(12)

式中：C(A)與C(A)iced分別為任意的結冰前后飛機的性能、穩(wěn)定性與控制參數(shù)或其導數(shù)；η為飛機結冰程度參數(shù)，取決于飛機的固有參數(shù)和結冰條件，其值越大，積冰后氣動參數(shù)的變化也就越大，表明飛機積冰情況越嚴重，根據(jù)η隨云層參數(shù)的變化曲線，積冰程度參數(shù)的變化范圍大致在0～0.3之間[31]；kC(A)為結冰對飛機氣動參數(shù)的影響參數(shù)，對于某架特定的飛機來說為一常數(shù)，其值可以通過數(shù)值仿真計算或飛行試驗獲得。

1.3.2 非對稱結冰影響模型

在實際飛行過程中，考慮到冰形形成的隨機性及除冰系統(tǒng)工作后冰形脫落的隨機性，并不能保證左右機翼氣動特性的一致。特別是當機翼除冰系統(tǒng)一側發(fā)生故障，導致左右機翼出現(xiàn)升力、阻力差，并由此產(chǎn)生滾轉及偏航力矩時，極易引發(fā)飛行事故。這種結冰情形不同于均勻結冰，因此，有必要單獨構建出非對稱結冰時的模型。

Lampton和Valasek[32-33]指出左右機翼的非對稱結冰主要會導致左右機翼非對稱升力與阻力的產(chǎn)生，從而產(chǎn)生附加的偏航與滾轉力矩，基于此，可以建立非對稱結冰情形下的積冰影響模型。文中假設右側機翼的除冰系統(tǒng)出現(xiàn)故障不能正常除冰，則左右機翼的升力及阻力差值為

(13)

式中：CL和CLice分別為飛機結冰前后的升力系數(shù)；CD和CDice分別為結冰前后的阻力系數(shù)。此時，兩側機翼的升力差產(chǎn)生了正的滾轉力矩，即使飛機向右滾轉的力矩；阻力差產(chǎn)生了正的偏航力矩，即使飛機向右偏航的力矩。附加的滾轉和偏航力矩可表示為

(14)

式中：dmgc為平均幾何弦長位置到飛機中心線的距離；Q為動壓；Sw為機翼面積。將式(13)和式(14)代入飛機6-DOF運動方程中的力和力矩項，通過仿真計算，即可得到左右機翼非對稱結冰時的飛機狀態(tài)變化。

1.3.3 結冰對其他參數(shù)的影響

結冰不僅會導致飛機氣動參數(shù)的改變，同時還會改變一些飛行安全有關的參數(shù)的限制值。如結冰會引起失速迎角降低[34]、最小平飛速度增大[35]等。

這些與安全相關的參數(shù)改變可以事先通過計算以數(shù)據(jù)庫的形式存儲在計算機中，對駕駛員安全飛行邊界進行計算時再進行調(diào)用。

需要指出的是上文所述的結冰影響模型是屬于目前國內(nèi)外比較常用的結冰后氣動參數(shù)估算模型，由于其包含了一定的物理意義，既考慮了飛機特性又考慮了結冰條件的影響，具有一定的通用性，該方法適合作為結冰后氣動參數(shù)尚未獲取的情況下對飛行動力學特性進行初步的分析。若要較為準確地預測出結冰后的飛行動力學特性，則需要借助于高精度數(shù)值模型、風洞試驗或真實試飛等手段來對特定的飛機預先獲取結冰后準確的氣動力方可進行。

2 單個飛行情形下飛行安全譜及風險量化

2.1 飛行安全譜計算

飛行風險的發(fā)生往往伴隨著飛行參數(shù)的異常變化，通過對單個飛行情形下的飛行參數(shù)數(shù)據(jù)變化情況進行分析，可以預判出飛行風險事件的發(fā)生與否。本節(jié)將飛行參數(shù)數(shù)據(jù)在不同區(qū)間上的值以色彩化的形式呈現(xiàn)，通過對單個飛行情形下不同飛行參數(shù)的色譜信號疊加，來得到該飛行情形下的飛行安全譜。安全譜的獲得可以很直觀地看出單個飛行情形下飛行風險的演化過程。

當前對于飛行安全參數(shù)的限制的描述往往是確定性的，如某型飛機的手冊中規(guī)定飛機的最大允許迎角α在馬赫數(shù)Ma<0.55時，不得超過16°。其潛在的含義意味著迎角小于16° 即是安全的，大于16° 將會導致飛行事故的發(fā)生。Burdun[36]認為實際中人們對于這種限制的認識是帶有模糊性質(zhì)的，如迎角達到15° 時同樣也是非常危險的狀態(tài)。為此，他們提出了用色彩來對飛行風險進行表示的方法，但他們所提出的方法僅能用來表示飛行參數(shù)處于風險狀態(tài)，而不能表示此時飛行參數(shù)是處于正的風險區(qū)間還是負的風險區(qū)間。為此，本文提出了一種考慮飛行參數(shù)風險區(qū)間正負性的表示方法，該方法可以更加合理地表示飛行參數(shù)變化及飛行風險的變化情況。

表1 與安全相關的飛行參數(shù)色彩化區(qū)間Table 1 Colored interval of safety related flight data

表格中舵面操縱量淺灰/深灰色對應的值為舵面的極限偏轉角度。對于舵面操縱量δe、δa和δr而言，舵面偏轉至極限位置飛行事故不一定會發(fā)生，其風險度設定為與淺紅/深紅一致，這里只是為了區(qū)分舵面操縱是否飽和，這一點與其他狀態(tài)參數(shù)不同。

通過人-機-環(huán)系統(tǒng)仿真，得到每一個關鍵的飛行參數(shù)在預測時間段內(nèi)的變化情況，根據(jù)表1中飛行參數(shù)的色彩化區(qū)間分布得到這些關鍵飛行參數(shù)的色譜圖Ci(t)，圖3所示為根據(jù)迎角在預測時間段內(nèi)的變化情況得到的迎角安全譜。

飛行安全是典型的木桶理論，任何一個關鍵飛行參數(shù)的異常都可能會導致飛行風險事件的發(fā)生，故最終的飛行安全譜上每一時刻對應的風險色與所有飛行參數(shù)在該時刻處在最危險狀態(tài)的飛行參數(shù)的風險色相同。以某型運輸類飛機在平飛狀態(tài)下以指定的航跡俯仰角(μc=8°)與滾轉角(φc=35°)協(xié)調(diào)爬升轉彎為例，計算出該飛行情形下各安全相關的飛行參數(shù)及整個飛行情形的飛行安全譜如圖4所示。其中，每個飛行參數(shù)對應的安全譜為該飛行參數(shù)數(shù)據(jù)在預測時間段內(nèi)的色譜，最下面一行表示該飛行情形整個時間段內(nèi)的風險變化情況，定義為該飛行情形總的飛行安全譜，其在每一時刻的風險值與該時刻處在最危險狀態(tài)的與安全相關的飛行參數(shù)風險值相同。

圖3 迎角變化曲線及相應的風險色譜圖
Fig.3 Time history and safety spectra for angle of
attack (AOA)

圖4 單個飛行情形的飛行安全譜(μ=8°, φ=35°)
Fig.4 Safety spectra for single flight condition (μ=8°, φ=35°)

由于淺灰/深灰、淺紅/深紅、淺黃/深黃只是表示風險區(qū)間的正負，其代表的風險值是相同的，故在計算總的飛行安全譜時，淺灰/深灰同用黑色表示，淺紅/深紅同用紅色表示，淺黃/深黃同用黃色表示。

由圖4的計算結果可知，在駕駛員按照目標飛行指令操縱時，整個飛行階段操縱舵面δe、δa和δr均保持在一個較低的風險水平，飛行風險偏高的階段，主要體現(xiàn)在操縱初期的迎角及法向過載值略微偏高；操縱后期爬升率以及滾轉角度維持在一個比較高的水平，但總體來說飛行仍是安全的，但除非有必要駕駛員最好能夠?qū)L轉角與高度變化率維持在一個較低的值。

2.2 飛行風險的量化

圖4中由各種顏色構成的飛行安全譜，代表了整個時間段內(nèi)的風險度演變過程。然而，不同飛行情形下的飛行安全譜都是不同的，如何根據(jù)安全譜中所含的風險信息進行不同飛行情形風險度的對比，是本章的研究內(nèi)容。提出了單個飛行情形風險量化的方法：采用對各種風險色所占百分比賦權相加的方法來得到該情形下飛行風險值。

假定風險色黑、紅、黃、綠分別占整個預測時間段的百分比為Pk、Pr、Py和Pg，將每個風險色代表的風險值定義為Vk、Vr、Vy和Vg，那么在整個預測時間段內(nèi)的風險值可表示為

R=PkVk+PrVr+PyVy+PgVg

(15)

考慮到安全譜中一旦有黑色區(qū)域出現(xiàn)表示至少有一個飛行參數(shù)的值超出了安全飛行的最大允許值，即認為飛行事故發(fā)生，這是飛行中要極力避免的，可將黑色代表的風險值Vk取一個較大的值，這樣一旦出現(xiàn)飛行事故，得到的飛行風險明顯要大于其他飛行情形。文中各風險色代表的風險值分別為：Vk=30、Vr=4、Vy=2、Vg=1。以圖4所示的飛行情形為例，得到該飛行情形的風險值為1.982。

3 飛行安全窗的構建

安全譜的獲得可以很直觀地看出飛機以某個目標航向指令飛行的風險演變過程。而單個飛行情形下風險的量化則是構建飛機在整個操縱空間的飛行風險拓撲云圖的基礎。以指令滾轉角與指令航向角構成的操縱空間為計算區(qū)域，將計算區(qū)域離散成許多計算單元，每個單元對應一個飛行情形，可以通過前述的方法計算出該飛行情形的風險值，進而可得到整個操縱空間風險場的拓撲云圖。顯然，每個計算單元之間是獨立的，風險值的計算過程相互之間并不干擾，因此可以通過并行計算的方式來加快計算進程。

對整個飛行操縱范圍的風險度進行計算需耗費一定的時長，采用并行仿真手段可以大大縮減計算時間，提高預測的時效性。本文采用的是基于MATLAB/Simulink平臺的并行飛行仿真方法，來解決大量蒙特卡羅飛行仿真的耗時性問題，為計算飛行動力學的實施提供了便利。其實施過程如下：

1) 建立飛機本體運動學模型，并計算在設定初始飛行狀態(tài)下的配平參數(shù)。

2) 在Simulink環(huán)境下搭建人-機-環(huán)復雜系統(tǒng)飛行動力學模型仿真平臺，仿真初始條件設定為配平值。

3) 將Simulink的仿真模式設置為在Rapid Accelerator模式下運行，并指定可調(diào)變量，如φc和μc等。

4) 建立該模型的Rapid Accelerator目標對象，初始化可調(diào)變量構成的結構體。

5) 在并行計算命令parfor下運行仿真模型，得到在整個操縱空間上每個計算單元的安全譜。

6) 依據(jù)式(12)計算出飛機在整個操縱范圍上每個計算單元對應的風險值，并在此基礎上，做出飛機操縱范圍的風險拓撲云圖。

需要指出的是，飛行事故的風險值明顯要大于一般飛行情形。例如在干凈構型下，當指令滾轉角φc=45°，指令航跡俯仰角μc=16° 時，通過上述方法得到的該飛行情形的風險值為18.447。過高的風險值會導致在建立風險拓撲云圖時，風險值較低的飛行情形的風險色幾乎為同一顏色，在拓撲圖中區(qū)別不明顯。為此，將所有飛行風險值大于某個限定值的飛行情形的飛行風險設定為該限定值，這樣使得風險色的分離度更加明顯。文中將該限定值取為4.5。根據(jù)2.2節(jié)中風險值的設定可知，只要總的風險值R>4，即意味著預測時間段內(nèi)出現(xiàn)了參數(shù)超限的情形，將限定值定為4.5既保證了對事故狀態(tài)點的區(qū)分，也確保了風險拓撲圖中低風險狀態(tài)點的風險色的分離度。

以飛機在H0=2 000 m，V0=120 m/s水平勻速直線飛行的初始飛行狀態(tài)下來計算飛機的安全飛行范圍。計算在曙光A840r-G服務器上實施，硬件為64核AMD Opteron 6276 2.3 GHz處理器，系統(tǒng)運行內(nèi)存為128 GB，仿真運行環(huán)境為MATLAB 2014a，通過測試，得到并行和串行兩種飛行仿真的計算信息對比如表2所示。

安全窗的計算結果如圖5所示，圖中不同的顏色代表了不同的風險度，圖右邊的色例為不同的風險值R對應的風險色。

表2并行與串行方法的計算條件與計算時間

Table2Computationconditionandtimeofparallelandserialmethods

ComputationmodeCalculationrangeNodenumberComputa?tiontime/sSerialcomputationμc∈[-6∶2∶18]?c∈[-55∶5∶55]299293．74Serialcomputationμc∈[-6∶0．5∶18]?c∈[-55∶2∶55]27442866．64Parallelcomputa?tion32workersμc∈[-6∶2∶18]?c∈[-55∶5∶55]29941．49Parallelcomputa?tion32workersμc∈[-6∶0．5∶18]?c∈[-55∶2∶55]2744232．79

圖5 干凈外形下的飛行安全窗(初始狀態(tài)H=2 000 m，V=120 m/s)
Fig.5 Flight safety windows for clean configuration
(trimmed at H=2 000 m，V=120 m/s)

4 案例分析

結冰會導致飛機穩(wěn)定性和操縱性變差，在飛機設計初期，借助于風洞試驗或者流場仿真等手段，獲得在不同結冰情形下飛機空氣動力學特性改變的前提下，通過前文描述的方法便可獲取飛行安全操縱范圍，進而為結冰條件下的邊界保護系統(tǒng)、飛行控制律、駕駛員情景感知增強系統(tǒng)等的設計提供一定的指導意義。本節(jié)以飛機遭遇機翼對稱結冰和非對稱結冰兩種飛行情形為例進行研究。飛行初始條件均設定為H0=2 000 m，V0=120 m/s，飛機保持平飛。

4.1 對稱結冰情形下的安全窗

4.1.1 安全窗的計算

假定結冰嚴重程度參數(shù)η=0.1，計算范圍φc∈[-55∶2∶55]，μc∈[-6∶0.5∶18]，可得計算節(jié)點數(shù)為46×49=2 254個，機翼對稱結冰時飛機的安全窗如圖6所示。

圖6 對稱結冰情形下的飛行安全窗(初始狀態(tài)H=2 000 m，V=120 m/s)
Fig.6 Flight safety windows for symmetric icing condition (trimmed at H=2 000 m，V=120 m/s)

從仿真結果可以很明顯地看出，結冰導致了飛機安全操縱范圍縮減。當飛機處于干凈構型時，縱向上，飛機最大正向航跡俯仰角大于18°，最大滾轉角約為53.5°。當飛機表面積聚冰形時，飛機在縱向和橫向上的安全操縱范圍分別縮減至小于12°和小于42°。

4.1.2 致災機理分析

為研究對稱結冰情形下飛行事故發(fā)生的致災機理，在圖6中選取4個典型狀態(tài)點((a)～(d))來進行研究。這4個狀態(tài)點的參數(shù)及其對應飛行

安全譜如圖7所示。

從圖7所示的仿真結果來看，發(fā)生對稱結冰時，飛機不管是僅目標航跡俯仰角設定過高(如圖7(a) 所示)，還是僅目標滾轉角設定過高(如圖7(c)所示)，或者是目標航跡俯仰角與目標滾轉角同時過高(如圖7(b)和圖7(d)所示)，往往伴隨著：迎角超過其右邊界，即大于失速迎角；飛行速度超過其左邊界，即小于最小飛行速度；高度變化率超過其左邊界，即飛機下降速度過快。即使目標航跡俯仰角設定的是正值，飛機最終也會陷入高度下降率變化過快的情形。

從事故演化過程來看，往往是由于飛機迎角超過結冰后的失速迎角，同時由于結冰帶來的氣動性能的惡化，飛機的推力不足以維持安全飛行的需要，導致飛行速度減小到其極限值，隨后飛機開始急速下降，此時很難避免事故的發(fā)生。因此，為確保對稱結冰情形下的飛行安全，首先要確保飛機處在安全的飛行范圍，駕駛員或者飛控系統(tǒng)還應當重點關注飛行過程中迎角、速度和高度變化率這3個安全關鍵參數(shù)的變化情況。

圖7 對稱結冰情形下典型事故狀態(tài)點的飛行安全譜
Fig.7 Flight safety spectrum of typical accident state points under symmetric icing conditions

4.2 非對稱結冰情形下的安全窗

4.2.1 安全窗的計算及分析

假定飛機在上述初始飛行條件及結冰氣象條件下，右側機翼除冰系統(tǒng)出現(xiàn)故障時，根據(jù)1.3.2節(jié)中所建立的非對稱結冰模型，通過文中所提出的飛行安全窗的計算理論，得到飛機在該情形下飛行風險的安全窗如圖8所示。

圖8 非對稱結冰情形下的飛行安全窗——右側除冰系統(tǒng)故障(初始狀態(tài)H=2 000 m，V=120 m/s)
Fig.8 Flight safety windows for asymmetry icing conditions—right wing half fatigue (trimmed at H=2 000 m，V=120 m/s)

從計算的結果來看，非對稱結冰情形下，結冰不但導致安全范圍的縮減，同時當結冰的不對稱性超過一定程度時，飛機的安全飛行范圍還出現(xiàn)了明顯的不對稱現(xiàn)象：對于綠色安全飛行區(qū)域，當飛機爬升時，飛機偏向除冰系統(tǒng)故障那一側的安全飛行范圍相對大一些，而當飛機下降時，飛機偏向除冰系統(tǒng)正常工作的一側進行操縱安全性相對大一些；對于安全窗整體而言，飛機向除冰系統(tǒng)正常一側滾轉時的安全飛行范圍比飛機向除冰系統(tǒng)故障一側滾轉時的安全飛行范圍要大些。下面分析造成結冰安全操縱范圍不對稱現(xiàn)象的原因。

1) 綠色飛行區(qū)域不對稱原因分析。

為分析非對稱結冰情形下，綠色飛行區(qū)域不對稱現(xiàn)象的原因，在圖8中選擇4個狀態(tài)點進行分析，它們的位置如圖8中藍色三角點所示，其對應的安全譜如圖9所示。

從圖9中對比分析的結果來看，造成機翼非對稱結冰時的非對稱安全飛行范圍的主要因素是由于飛機向左或向右滾轉時產(chǎn)生的高度變化率不同導致的：當飛機滾轉爬升時，向左滾轉爬升時的高度變化率(正值)要大于向右滾轉爬升時的高度變化率，導致向左滾轉爬升時的高度變化率更容易超限，因此飛機向右滾轉爬升時的安全飛行范圍要大些；而當飛機滾轉下滑時，向左滾轉下滑時的高度變化率(負值)要慢于向右滾轉下滑時的高度變化率，導致向右滾轉下滑時的高度變化率更容易超限，因此飛機向左滾轉下滑時的安全性相對而言更高些。

2) 安全窗整體不對稱原因分析。

飛機發(fā)生非對稱結冰時，除冰系統(tǒng)工作正常一側的可用飛行范圍明顯大于除冰系統(tǒng)故障的一側，為分析原因，在圖8中選擇4個狀態(tài)點進行對比，它們的位置如圖8中紫色三角點所示。這4個狀態(tài)點對應的飛行安全譜如圖10所示。

圖9 非對稱結冰情形下典型狀態(tài)點對應的飛行安全譜(φ=±25°)
Fig.9 Flight safety spectrum of typical state point under asymmetry icing conditions (φ=±25°)

圖10(a)和圖10(c)中，仿真時長均為60 s，而圖10(b)和圖10(d)中仿真均停止于滾轉角大于150°(滾轉角達到不可逆轉的值)。對比分析這4種飛行狀態(tài)的飛行安全譜可知，對于非對稱結冰情形而言，飛機向著除冰系統(tǒng)故障的一側滾轉機動時，當目標滾轉角設定較大時飛機很難達到并保持該目標滾轉角，極易引發(fā)滾轉角迅速發(fā)散的情形，進而導致飛行事故的發(fā)生。當飛機向著除冰系統(tǒng)工作正常的一側滾轉機動時，對于同樣大小的目標滾轉角，飛機仍可達到并保持該飛行狀態(tài)進行飛行，各安全關鍵參數(shù)變化穩(wěn)定。因此，對于非對稱結冰情形而言，尤其要關注飛機向除冰系統(tǒng)故障一側滾轉角的變化情況。

圖10 非對稱結冰情形下典型狀態(tài)點對應的飛行安全譜(φ=±42°)
Fig.10 Flight safety spectrum of typical state point under asymmetry icing conditions (φ=±42°)

需要指出的是，雖然飛機向左滾轉的飛行范圍與干凈構型下的情況相差不大，但相比于干凈構型的飛機，駕駛員模型中增益值相對較小，也就是說駕駛員在實際操作中要采用比較柔和的操縱方式才能保證飛行安全。

4.2.2 致災機理分析

為研究非對稱結冰情形下飛行事故發(fā)生的致災機理，選取4個典型的飛行事故狀態(tài)點來進行研究，它們的位置如圖8中白點所示。這4個狀態(tài)點對應的飛行安全譜如圖11所示。

圖11 非對稱結冰情形下典型事故狀態(tài)點的飛行安全譜
Fig.11 Flight safety spectrum of typical accident state points under asymmetry icing conditions

圖11(a)和圖11(b)中，指令滾轉角設定為向除冰系統(tǒng)正常一側以較大角度偏轉，仿真結果表明滾轉角出現(xiàn)了大幅發(fā)散振蕩，直至超出邊界；副翼正偏(即右副翼下偏)的操縱效能不足，使得駕駛員很難使?jié)L轉角維持在一個較大的角度；同時還伴隨著迎角的超限。圖11(c)表明當目標航跡俯仰角設定過大時，飛機的橫向滾轉會向著除冰系統(tǒng)故障的一側而無法糾正，同時迎角也會一直保持在一個較危險的區(qū)間直至飛行事故的發(fā)生。對于指令滾轉角設定為向除冰系統(tǒng)故障一側以較大角度偏轉時，如圖10(b)、圖10(d)和圖11(d)中所示(對應不同的航跡俯仰角)，滾轉角很快便發(fā)散至不可挽回的地步。

仿真結果表明，飛機出現(xiàn)非對稱結冰時，引起飛行事故的原因主要是滾轉角超限，而且向著除冰系統(tǒng)故障的一側滾轉角很容易超限，即使是目標航跡滾轉角設定為向除冰系統(tǒng)工作正常的一側偏轉，飛機也會由于副翼正向偏轉操縱效能的不足而最終演化為向著除冰系統(tǒng)故障的一側異常偏轉。此外，迎角也很容易在事故演化過程中達到其極限值。而此時飛機飛行速度、高度變化率基本處在安全的范圍之內(nèi)，只有在仿真的末期，才會出現(xiàn)超限的可能，這與對稱結冰情形下的致災機理有所差異。

因此，在非對稱結冰情形下，駕駛員或者飛行控制系統(tǒng)還應當首要關注飛行過程中的滾轉角和迎角這兩個安全關鍵參數(shù)的變化情況。

5 結 論

基于人-機-環(huán)復雜系統(tǒng)動力學仿真，研究了飛行安全窗的構建方法。建立了蒙特卡羅并行飛行仿真平臺，計算了不同結冰情形下的飛行安全窗口，揭示了不同結冰情形下，由于駕駛員可能的操縱失誤帶來的飛行事故發(fā)生機理，并給出駕駛員應對策略建議。

結冰對于不同飛機的氣動性能影響必然是不同的，文中只是以某典型構型的運輸類飛機為例對安全窗的構建方法及結冰致災機理進行研究分析。對于不同的飛機其在結冰情形下安全窗的變化情況必然會受到結冰影響模型準確性以及飛機固有特性的影響。本文側重于對風險譜的建立及風險的量化以及安全窗的建立方法進行研究，所提出的理論方法及研究結論對于結冰條件下致災與防護方面的研究具有一定的借鑒意義和價值。

計算飛機在不同結冰情形下的風險安全窗，可以以數(shù)據(jù)庫的形式存儲在機載計算機中。當飛機處于不同的飛行狀態(tài)下，遭遇不同程度的結冰時，可以通過插值計算方式得到飛機在該狀態(tài)下的安全窗，從而為駕駛員提供指示，提升結冰情形下駕駛員的情景感知能力。

駕駛員的操縱具有隨機性，在安全的飛行范圍內(nèi)也存在著駕駛員的粗暴操縱導致飛行事故發(fā)生的可能。文中在計算時并沒有考慮這一點。在對駕駛員模型進行建模時，考慮的只是盡可能地使飛機按照指定的飛行指令來飛行。因此得到的安全飛行范圍只是一個參考范圍，駕駛員在進行操縱時根據(jù)提供的安全窗可以對自己的操縱裕度產(chǎn)生直觀的感受，從而保證飛行安全。

此外，文中所提出的方法還可以用在飛行控制系統(tǒng)設計(控制律參數(shù)優(yōu)化)、驗證控制律、飛機性能計算、飛機遭遇各種不利情形時(如單發(fā)失效、多不利因素情形)的操縱范圍。

[1] BURDUN I Y. A method for accident reconstruction and neighborhood analysis using an autonomous situational model of flight and flight recorder data[C]//Advances in Aviation Safety Conference and Exposition. Warrendale, PA: SAE International, 1999: 1-13.

[2] Boeing. Airplanes statistical summary of commercial jet airplane accidents: Worldwide operations 1959-2005[R]. Chicago, IL: The Boeing Company, 2006: 1-25.

[3] Airbus Customer Services. Human performance: Enhancing situational awareness: FLT_OPS-HUM-PERF-SEQ 06-REV 01[R]. Blagnac Cedex: Airbus, 2007.

[4] BURDUN I Y. The intelligent situation awareness and forecasting environment (The S.A.F.E. Concept) a case study[C]//Advances in Aviation Safety Conference and Exposition. Warrendale, PA: SAE International, 1998: 1-15.

[5] ZOLGHADRI A. Early warning and prediction of flight parameter abnormalities for improved system safety assessment[J]. Reliability Engineering and System Safety, 2002, 76(1): 19-27.

[6] ROUWHORST W F J A, MARSMAN A P L A. A piloted investigation of an integrated situation awareness system (ISAS): AIAA-2006-6723[R]. Reston: AIAA, 2006.

[7] BORST C, SJER F A, MULDER M, et al. Ecological approach to support pilot terrain awareness after total engine failure[J]. Journal of Aircraft, 2008, 45(1): 159-171.

[8] BORST C, GROOTENDORST F H, BROUWER D I K, et al. Design and evaluation of a safety augmentation system for aircraft[J]. Journal of Aircraft, 2014, 51(1): 12-22.

[9] BRAGG M B, BASAR T, PERKINS W R, et al. Smart icing systems for aircraft icing safety[C]//40th AIAA Aerospace Sciences Meeting & Exhibit. Reston: AIAA, 2002.

[10] DETERS R, DIMOCK G A, SELIG M S. Icing encounter flight simulator with an integrated smart icing system: AIAA-2002-4599[R]. Reston: AIAA, 2002.

[11] GINGRAS D R, BARNHART B, RANAUDO R, et al. Envelope protection for in-flight ice contamination[C]//47th AIAA Aerospace Sciences Meeting including The New Horizons Forum and Aerospace Exposition. Reston: AIAA, 2009.

[12] GINGRAS D R, BARNHART B, RANAUDO R, et al. Development and implementation of a model-driven envelope protection system for in-flight ice contamination[C]//AIAA Guidance, Navigation, and Control Conference. Reston: AIAA, 2010.

[13] RANAUDO R, MARTOS B, NORTON B, et al. Piloted simulation to evaluate the utility of a real time envelope protection system for mitigating in-flight icing hazards[C]//AIAA Atmospheric and Space Environments Conference. Reston: AIAA, 2010.

[14] 徐忠達, 蘇媛, 曹義華.平尾積冰對飛機縱向氣動參數(shù)的影響[J]. 航空學報, 2013, 34(7): 1563-1571.

XU Z D, SU Y, CAO Y H. Effects of tailplane icing on aircraft longitudinal aerodynamic parameters[J]. Acta Aeronautica et Astronautica Sinica, 2013, 34(7): 1563-1571 (in Chinese).

[15] DONG Y, AI J. Research on inflight parameter identification and icing location[J]. Aerospace Science and Technology, 2013, 29(1): 305-312.

[16] 應思斌， 艾劍良. 飛機結冰包線保護對開環(huán)飛行性能影響與仿真[J]. 系統(tǒng)仿真學報, 2010, 22(10): 2273-2301.

YING S B, AI J L. Simulation of aircraft flight envelope protect in icing encounters effects on open loop dynamic [J]. Journal of System Simulation, 2010, 22(10): 2273-2301 (in Chinese).

[17] 劉東亮， 徐浩軍， 李嘉林， 等. 飛行結冰后復雜系統(tǒng)動力學仿真與風險評估[J]. 系統(tǒng)仿真學報, 2011, 23(4): 643-647.

LIU D L, XU H J, LI J L, et al. Dynamic simulation study of stalling in wing icing conditions and risk evaluation[J]. Journal of System Simulation, 2011, 23(4): 643-647 (in Chinese).

[18] 王明豐, 王立新, 黃成濤. 積冰對飛機縱向操穩(wěn)特性的量化影響[J]. 北京航空航天大學學報, 2008, 34(5): 592-595.

WANG M F, WANG L X, HUANG C T. Computational effects of ice accretion on aircraft longitudinal stability and control[J]. Journal of Beijing University of Aeronautics and Astronautics, 2008, 34(5): 592-595 (in Chines).

[19] 張智勇. 結冰飛行動力學特性與包線保護控制律研究[D]. 南京： 南京航空航天大學, 2006.

ZHANG Z Y. Research on iced aircraft flight dynamics characteristics and envelope protection control law[D]. Nanjing: Nanjing University of Aeronautics and Astronautics， 2006 (in Chinese).

[20] TRUJILLO A, GREGORY I. Pilot preferences on displayed aircraft control variables[J]. Lecture Notes in Computer Science, 2013, 8020(1): 203-211.

[21] SIBILSKI K, LASEK M, LADYZYNSKA-KOZDRAS E, et al. Aircraft climbing flight dynamics with simulated ice accretion: AIAA-2004-4948[R]. Reston: AIAA, 2004.

[22] SONNEVELDT L. Nonlinear F-16 model description[R]. The Netherlands: Delft University of Technology, 2010.

[23] COOK M V. Flight dynamics principles[M]. Burlington, MA: Butterworth-Heinemann, 2007: 66-95.

[24] 吳森堂. 飛行控制系統(tǒng)[M]. 北京: 北京航空航天大學出版社, 2013: 54-64.

WU S T. Flight control system[M]. Beijing: Beihang University Press, 2013: 54-64 (in Chinese).

[25] HUESCHEN R M. Development of the transport class model (TCM) aircraft simulation from a sub-scale generic tra nsport model (GTM) simulation: NASA/TM-2011-217169[R]. Washington, D.C.: NASA, 2011.

[26] HANKE C R, NORDWALL D R. The simulation of a jumbo jet transport aircraft. Volume 2: Modeling data/detail: NASA-CR-114494[R]. Washington, D.C.: NASA, 1970.

[27] DOGAN A, KAEWCHAY K. Probabilistic human pilot approach: Application to microburst escape maneuver[J]. Journal of Guidance, Control, and Dynamics, 2007, 30(2): 357-369.

[28] GAO Z, GU H. Simulation of microburst escape with probabilistic pilot model[J]. Lecture Notes in Computer Science, 2011, 7030(1): 663-670.

[29] DOGAN A, KABAMBA P T. Escaping microburst with turbulence altitude, dive, and pitch guidance strategies[J]. Journal of Aircraft, 2000, 37(3): 417-426.

[30] 高金源， 李陸豫， 馮亞昌. 飛機飛行品質(zhì)[M]. 北京: 國防工業(yè)出版社, 2003: 135-140.

GAO J Y, LI L Y, FENG Y C. Aircraft handling quality[M]. Beijing: National Defence Industry Press, 2003: 135-140 (in Chinese).

[31] BRAGG M, HUTCHISON T, MERRET J . Effect of ice accretion on aircraft flight dynamics[C]//38th AIAA Aerospace Sciences Meeting & Exhibit. Reston: AIAA, 2000.

[32] LAMPTON A, VALASEK J. Prediction of icing effects on the lateral directional stability and control[J]. Aerospace Science and Technology, 2012, 23(1): 305-311.

[33] LAMPTON A, VALASEK J. Prediction of icing effects on the coupled dynamic response of light airplanes[C]//AIAA Atmospheric Flight Mechanics Conference and Exhibit. Reston: AIAA: 2007.

[34] SHARMA V, VOULGARIS P G, FRAZZOLI E. Aircraft autopilot analysis and envelope protection for operation under icing conditions[J]. Journal of Guidance, Control, and Dynamics, 2004, 27(3): 454-465.

[35] Federal Aviation Administration. Pilot guide flight in icing conditions: AC 91-74A[R]. Washington, D.C.: FAA, 2007.

[36] BURDUN I Y. Automated planning, exploration and mapping of complex operational domains of flight using multifactor situational trees[J]. SAE International Journal of Aerospace, 2011, 4(2): 1149-1175.

(責任編輯： 李明敏)

URL：www.cnki.net/kcms/detail/11.1929.V.20161103.1633.002.html

Developmentofflightsafetywindowinicingconditionsbasedon

complexdynamicssimulationPEIBinbin1,XUHaojun1,XUEYuan1,*,LIZhe1,LIUDongliang2

1.AeronauticsandAstronauticsEngineeringCollege,AirForceEngineeringUniversity,Xi’an710038,China2.BeijingAeronauticalEngineeringTechnologyResearchCenter,Beijing100076,China

Studiesonenhancingsituationalawarenessofpilotsinicingconditionsarerelativelyrare.Currentmethodsnormallypredictoccurrenceofaccidentsbyestimatingwhetherthesafetyrelatedparametersexceedtheirlimitations.Complexdynamicsofthepilot-vehicle-icingeffectmodelisestablished.Safetyspectrumofasingleflightconditionisobtainedbycomparingtheriskdegreeoftheflightdatathroughthepredicttimeinterval,andthecoloredriskvaluefortheflightconditionisfurtheracquired.Safetywindowforflightsafetyinthewholemanipulationcanbecalculatedusingparallelflightsimulationplatform.Thesafetywindowsofsymmetricandasymmetricicingconditionsareresearched,andthedisaster-causingmechanismisanalyzed.Simulationresultsshowthaticewillleadtoshrinkofthesafetyflightscope,andtheasymmetricalicewillleadtotheunevensafetywindow.Theproposedmethodcanprovidetheoreticalsupportforenhancementofsituationalawarenessindifferentadverseevents,andengineeringtoolforoptimizingtheaircraftperformanceforaircraftdesigners.

aircrafticing;situationalawareness;computationalflightdynamics;asymmetricicing;parallelflightsimulation;safetywindow

2016-08-23；Revised2016-09-20；Accepted2016-10-25；Publishedonline2016-11-031633

s：NationalBasicResearchProgramofChina(2015CB755800);NationalNaturalScienceFoundationofChina(61374145,U1333131)

.E-mailszxy1986@163.com

2016-08-23；退修日期2016-09-20；錄用日期2016-10-25； < class="emphasis_bold">網(wǎng)絡出版時間

時間：2016-11-031633

www.cnki.net/kcms/detail/11.1929.V.20161103.1633.002.html

國家“973”計劃 (2015CB755800)； 國家自然科學基金 (61374145，U1333131)

.E-mailszxy1986@163.com

裴彬彬， 徐浩軍， 薛源， 等． 基于復雜動力學仿真的結冰情形下飛行安全窗構建方法J． 航空學報,2017,38(2):520695.PEIBB,XUHJ,XUEY,etal.DevelopmentofflightsafetywindowinicingconditionsbasedoncomplexdynamicssimulationJ.ActaAeronauticaetAstronauticaSinica,2017,38(2):520695.

http://hkxb.buaa.edu.cnhkxb@buaa.edu.cn

10.7527/S1000-6893.2016.0274

V212

A

1000-6893(2017)02-520695-14