基于三維重建飛機氣動特性的飛行包線

蘇光旭， 張登成， 張久星

(1.空軍工程大學航空工程學院， 西安， 710038； 2.94916部隊， 南京， 211500; 3.93756部隊， 天津， 300131)

飛行包線是以飛行速度、飛行高度、載荷等為邊界用來表示飛機飛行范圍的封閉圖形[1]，是飛機飛行性能關鍵指標。飛行包線是確保飛行安全的基礎，現代戰斗機通過多種包線限制可以實現飛行員的“無憂慮”飛行。文獻[2]按照飛行任務將飛行包線分為基本飛行包線、小表速包線、大表速包線、發動機空中啟動包線、速度載荷包線、突發載荷包線以及安全彈射救生包線等。基本飛行包線是以飛機平飛狀態下的失速限制、理論升限，以及最大平飛速度限制所組成的飛機飛行邊界。

對于飛行包線的計算與應用，國內外學者做了大量研究。文獻[3]提出了多種評估軍用飛機飛行邊界的方法，可以通過評估未知飛機的飛行性能進而給出飛行邊界。文獻[4]和[5]基于飛行數據，通過可達平衡集的評價方法求取飛機所有狀態點的機動能力，給出了飛機的機動飛行包線。文獻[6]基于飛機的包線數據庫，提出了一種針對機翼故障的估算動態飛行包線以確保飛行安全的方法。文獻[7]使用現有的飛機幾何模型計算了飛機結冰數據庫，提出了一種基于神經網絡自適應動態逆的結冰飛機飛行安全邊界保護方法。綜合現有研究成果，飛機的氣動特性是飛行包線計算的關鍵，而幾何模型是氣動特性計算的基礎，很少有針對氣動參數未知的飛機建立精確飛行包線的工作先例。

綜上所述，本文是以氣動參數未知飛機的飛行安全和空戰仿真為研究背景，對飛機幾何模型進行三維重建，應用計算流體力學(computational fluid dynamics, CFD)分析了飛機的縱向氣動特性，基于所計算的氣動數據通過“簡單推力法”[8]確定平飛包線，增加最大平飛速度限制條件，給出了飛機的基本飛行包線。氣動特性分析與包線計算充實了該飛機的飛行性能，為飛機的飛行安全提供了技術支撐，為其動作庫構建和空戰仿真奠定了基礎。本文的研究方法，是一個完整的從飛機逆向建模到飛行包線計算的技術途徑。

1 飛機幾何外形三維重建

某幾何外形參數未知的飛機是一種翼身融合的常規氣動布局雙發重型戰斗機，本文綜合運用圖像明暗恢復形狀和工程圖重建法建立飛機的三維幾何模型，為數值模擬計算進行模型準備。其主要步驟介紹如下：

1)透視參數求解。由圖形學知識，通過照相機標定可以求解圖像的透視參數，包括相機焦距、視點距離、坐標系轉換矩陣等。如圖1所示，對于飛機照片一般使用“僅知兩個滅點的相機標點”[9]方法，通過優化迭代計算第三滅點的位置，從而求解各透視參數。

圖1 飛機照片的兩個滅點

2)圖像明暗恢復形狀建立局部三維曲面。使用圖像明暗(灰度)恢復形狀(shape from shading，SFS)可以初步求解觀察體曲面的初始高度值[10-11]。在已知反射特性的前提下，依據圖像灰度約束方程[12]求解模型的表面梯度，得到物體的表面高度，即可恢復曲面的外形。本文使用三次參數樣條曲線擬合[13]曲面的高度數據構建局部三維外形，如圖2飛機機頭錐的線框模型。

圖2 飛機機頭錐的重建過程

3)工程圖重建法整合飛機整機模型。飛機整機模型的三維重建是一個復雜繁瑣的過程，需要分成不同部件和曲面進行重建，本文根據平行投影原理通過工程圖重建法[14]進行整機模型的拼接。另外，對于氣動參數影響較大的部件如機翼等，需要使用翼型數據輔助建模。圖3給出了整合后飛機的線框模型以及模型的渲染效果。

圖3 整機線框模型與渲染效果

2 氣動特性計算

使用重建的飛機三維幾何模型，在不同巡航高度、迎角和馬赫數的飛行狀態下進行數值模擬計算，得到氣動數據并分析其縱向氣動特性。

2.1 網格劃分與計算條件

使用成熟的CFD計算軟件進行數值模擬計算。網格劃分(圖4)為非結構網格，網格總數922萬。邊界層內第1層網格高度控制在0.002 mm，以滿足機體表面黏性邊界層的計算要求[15]。

圖4 網格劃分

求解計算條件中，飛機表面使用無滑移壁面條件，外場壁面定義為壓力遠場條件，湍流模型使用對渦黏性系數修正的方法以適應不同區域流動的Spalart-Allmars模型。

2.2 網格無關性與算例驗證

為了進行網格無關性驗證，在Ma=0.8、α=5°、飛行高度h=10 km工況下，選用3套非結構網格計算整機的升力系數和阻力系數。網格數和計算結果見表1，網格劃分情況見圖5。通過對比可見網格量變化對計算結果影響微弱，為了兼顧計算精度和計算時間，選取中網格。

表1 不同網格量升力系數和阻力系數計算結果對比

圖5 不同網格量的網格劃分

采用AIAA阻力會議的標模DLR-F6翼身組合體作為驗證模型[16]來檢驗本文計算方法的正確性。選用相同的S-A湍流模型和計算條件，該翼身組合體網格劃分和計算結果如圖6所示，分別計算了機體的升力系數CL、阻力系數CD、機身俯仰力矩Cm，其中Exp為試驗數據，S-A為數值計算結果。可以看出，仿真計算結果與試驗數據吻合較好，滿足工程應用的要求。

圖6 模型算例驗證

2.3 氣動特性分析

在軟件計算精度驗證的基礎上，對圖3所示飛機模型計算不同飛行高度，飛機從亞聲速到超聲速不同迎角下的縱向氣動特性，下面以飛行高度h=10 km為例進行分析(見圖7)。

圖7 h=10 km縱向氣動特性

圖7給出飛機的升力系數CL、阻力系數CD和升阻比K(K=CL/CD)隨迎角α的變化擬合曲線。從不同飛行馬赫數下升力系數隨迎角變化規律可以看出，升力系數隨著迎角增大成S型增長，零升迎角在-2°左右，失速迎角在30°左右，說明該飛機的氣動布局具有較好的大迎角失速特性。就不同的飛行馬赫數來看，飛機在跨聲速飛行時升力系數變化率最大，而低馬赫數飛行時飛機的升力系數變化率最小。

從阻力系數隨迎角的變化規律可以看出，在小迎角飛行(α<10°)狀態下，阻力系數隨迎角增大增長的較為緩慢，而在較大迎角飛行狀態阻力系數隨迎角增大呈線性增長。就不同的飛行馬赫數來看，隨著飛機飛行馬赫數增大，整機的阻力系數先增大后減小，在Ma=1時阻力系數達到最大值，且增長率最大。

從升阻比隨迎角的變化規律可以看出，隨著迎角增大，升阻比呈現出先減小后增大再減小的趨勢，最大升阻比的峰值均在α=5°左右，而且馬赫數越小這種變化越明顯。最大升阻比7.726出現在Ma=0.6時，隨飛行馬赫數增大，最大升阻比逐漸減小，Ma=1.8時最大升阻比減小為4.04，說明飛機在小迎角亞聲速飛行狀態下具有較好的亞聲速巡航特性。

根據升力特性可知飛機具有較好的大迎角失速特性，為了分析其大迎角狀態下的流動機理，圖8給出了Ma=0.8時不同迎角下飛機縱向的渦強分布。可以看出，來流經過邊條翼時形成了明顯的脫體渦，這個最明顯的脫體渦的尺度表示了流經機翼上方湍流的宏觀尺寸，脫體渦的形狀不規則，大體環繞在機翼上方，沿機身縱軸尺度不斷增大。通過對比可以發現，隨著迎角不斷增大，脫體渦強度不斷增大，渦流動速度加快，加強了切洗效應[17]，機翼的低壓區和低壓程度不斷增大，這是升力增大的主要原因。當迎角為30°時，機翼上方脫體渦的尺度擴大，出現了破裂的趨勢，切洗效應減弱，此時的升力系數接近峰值，當迎角再增大時，對應圖7(a)升力系數變化規律，升力系數開始減小。

圖8 0.8 Ma時不同迎角渦強分布

3 飛行包線的計算

基于CFD計算的氣動參數，使用“簡單推力法”加平飛最大速度限制條件的方法分別求解飛機的最小平飛速度、靜升限以及最大平飛速度，給出飛機的基本飛行包線。

3.1 最小平飛速度

最小平飛速度是飛機在某一高度下能夠保持平飛的最小速度，是基本飛行包線的左邊界，即失速限制。根據平飛條件，求某一固定高度下的最小平飛速度，需聯立以下公式:

(1)

CL max=fMa

(2)

式中：CL max為飛機最大允許升力系數；v′為當地聲速；ρ為空氣密度，取決于飛行高度h；S為機翼面積。式(1)為飛機的平飛約束條件，即升力與重力相等。式(2)體現了飛機的氣動特性，可以根據CFD計算結果進行擬合得到。

選取一系列適當的Mai(i=1，2，…),根據平飛約束條件式(3)可解得選取馬赫數Mai所對應的最大升力系數CL max i，求CL max i～Mai擬合曲線與CFD計算結果擬合曲線的交點，可以得到在某一高度下的平飛最小速度，即Mamin=0.342。

(3)

圖9給出了h=10 km時2條曲線的擬合結果，圖中Polyfit表示多項式擬合結果，下同。

圖9 h=1 km CL max隨Ma變化關系

求得不同高度的平飛最小速度之后進行擬合即得飛行包線左邊界，見圖10。

圖10 不同高度下的最小飛行馬赫數

3.2 升限的確定

飛機的靜升限，又稱為理論升限，定義為飛機的最大上升率為0時的飛行高度，即所能保持平飛的最大高度。飛機上升率為飛行速度在豎直方向上的速度分量，根據定義有:

(4)

式中：θ為飛機的航跡俯仰角。

當飛機在鉛錘面內作定常運動時，假設飛機的迎角α不大且飛機的推重比較小，則飛機在鉛錘面內的質心動力學方程可以表述為:

(5)

式中：PAVL為發動機的可用推力；L為升力；D為飛行阻力。另根據平飛的限制條件，其所需推力Pr等于飛行阻力，經推導可得:

(6)

所以，當飛機的最大上升率為0時，有:

(7)

因此，計算固定馬赫數下飛機平飛時的升阻比K，當在某一飛行高度下K滿足上式時就可以得到對應飛行馬赫數下的理論升限。

下面以Ma=1.2為例說明計算過程：

1)根據飛機平飛的限制條件，求解當Ma=1.2時不同飛行高度下平飛所需升力系數CL,L=G，如表2所示。

表2 Ma=1.2不同高度下的平飛升力系數

2)根據CFD計算結果，將不同飛行高度下飛機升阻比K隨升力系數的變化規律進行多項式擬合，如圖11所示。根據擬合結果求解不同飛行高度下平飛所需升力系數對應的平飛升阻比KL=G，如表3所示。

圖11 不同高度下升阻比隨升力系數的變化

表3 Ma=1.2不同高度下的平飛升阻比

3)根據2)的計算結果，得出G/KL=G隨飛行高度的變化規律并進行擬合，見圖12。其實，這里的G/KL=G即為Ma=1.2時飛機的平飛所需推力Pr隨飛行高度的變化關系，將平飛所需推力與可用推力繪制在同一曲線圖上，曲線的右交點所對應的飛行高度即為飛機在Ma=1.2時的理論升限hmax=16 077 m。

圖12 可用推力與平飛所需推力變化

分別計算不同飛行馬赫數下的理論升限，即可得到飛行包線的上邊界，如圖13所示。

圖13 理論升限隨馬赫數的變化關系

3.3 最大平飛速度

根據平飛條件，由PAVL=Pr可以得到由推力限制的最大平飛速度，但是現代高性能戰斗機由于氣動外形的改進和發動機推力的提升，按照簡單推力法計算得到最大平飛速度一般會超過飛機所能承受的結構強度、剛度和溫度極限，因此，本文綜合運用動壓限制、溫度限制和推力限制的方法確定飛機的最大平飛速度。

3.3.1 動壓的限制

具有超聲速飛行能力的飛機在飛行高度h≤12 km時最大平飛速度主要受到所能承受的最大氣動載荷，即動壓的限制[3]。動壓限制通常以最大允許表速的形式對飛行員顯示，它體現了飛機機體結構的剛度和強度。文獻[18]提出了一種使用氣動彈性靜力學模型求解飛機最大氣動載荷的方法，較為繁瑣，本文結合與該飛機具有相似構造的同類飛機A以及外軍具有類似戰術指標飛機的飛行性能[19]估算該型飛機的最大允許表速。

飛機在海平面飛行時最大允許表速等于給定最大氣動載荷的飛行真速，根據換算關系可以計算出最大允許表速在不同飛行高度所對應的飛行馬赫數，如圖14的動壓限制所示。

3.3.2 溫度的限制

飛機在高速飛行時，氣動增溫使機體表面溫度急劇升高，當溫度過高時機體結構材料會產生機械性破壞，因此機體表面的溫度也限制了飛機的最大平飛速度。現代軍用飛機蒙皮所使用的材料主要有鋁、鎂合金，在一些特殊部位會使用少量的鈦合金或者以碳纖維為代表的高性能復合材料[20]，其中能夠承受溫度最低的是鋁合金類材料，以此作為機體表面溫度的限制值是合理的。

根據空氣動力學知識，氣動增溫主要取決于環境溫度Ti和Ma，有：

T0=Ti(1+0.2Ma)2

(8)

式中：T0為機體表面溫度；Ti取決于飛行高度。因此，可以通過最大溫度限制計算出不同飛行高度下的Mamax，溫度限制如圖14所示。

圖14 動壓和溫度限制的最大平飛速度

3.3.3 推力限制

推力的最大值也限制了最大平飛速度，根據平飛時P=D，當飛機的阻力系數取得最小值CD min時，飛行速度達到最大值，即有

(9)

經過計算可知，以推力為限制的飛行速度超過了動壓限制和溫度限制，對于飛行包線無意義。

3.4 對比驗證

將由圖10、圖13和圖14所計算的結果繪制到一張曲線圖上并求交集，即為飛機的基本飛行包線，如圖15所示。飛行包線的左邊界為失速限制，超過左邊界時飛機的升力將不足以維持其飛行高度；上邊界為靜升限，飛機到達靜升限時上升率為0，不能夠繼續爬升；右邊界為動壓限制和溫度限制，超過右邊界時飛機結構發生破壞而不能安全飛行。

圖15 基本飛行包線計算結果

將同類飛機A的飛行包線與計算結果繪制在一張曲線圖上進行對比(見圖16)。分析可知，本文計算結果與已知飛機A的飛行包線基本相符，說明了計算結果的可靠性；計算結果的左邊界偏保守，主要是由于數值模擬計算以及曲線擬合產生的誤差造成的；計算結果的上邊界及右邊界范圍更大，這主要是由于該飛機是在飛機A的基礎上改進升級而來，推力、結構強度都有相應的提升，計算得到的溫度極限值和最大表速較大，包線邊界向右上擴展，擁有更大的極限飛行速度。

圖16 飛行包線對比驗證

4 結論

本文通過三維逆向建模建立了飛機的三維幾何模型，通過計算流體力學的方法計算分析了其縱向氣動特性，并基于氣動數據計算了體現飛機基本飛行性能的飛行包線，為該飛機的飛行安全提供了技術支撐，為其空戰仿真進行了數據準備。具體結論如下：

1)氣動特性方面，采用翼身融合氣動布局的該型飛機在小迎角飛行狀態下具有較好的巡航性能，升阻比較高，并且飛機的大迎角飛行性能突出，失速迎角在30°左右；

2)飛行包線方面，經過對比驗證，計算得到的基本飛行包線能夠體現飛機的飛行性能，具有較高的可靠性；

3)本文所做工作探索出了一條從三維逆向建模到氣動特性計算、飛行邊界求解來研究外軍缺乏氣動參數飛機飛行性能的完整技術途徑，具有較高的通用參考價值，可以應用于其他先進飛行器。