結冰條件下飛機氣動/運動耦合特性

伍強，徐浩軍，魏揚，裴彬彬, 3，薛源

1. 空軍工程大學 航空工程學院，西安 710038 2. 中國人民解放軍 94639部隊，南京 211500 3. 西安電子科技大學 電子工程學院, 西安 710071

波音公司統計結果表明，飛機結冰是誘發飛行失控嚴重事故的三大環境因素之一。結冰會破壞飛機原有的氣動外形，結冰后的流動表現出強烈的非定常、非線性和隨機性，如邊界層轉捩位置的改變、大范圍分離流動、不同尺度旋渦的產生及脫落等。這種復雜流動使得飛機部件表面的氣動載荷處于不規則的變化之中，進而產生異常的氣動力和力矩。異常氣動特性作用下的飛機，其操縱性、穩定性等飛行力學特性被破壞，飛機的姿態和飛行狀態往往會產生不可控的變化，這種變化又反作用于飛機的流場和氣動特性，從而改變冰的演變過程。同時，結冰的動態增長又將導致飛機縱橫向氣動特性的耦合程度加劇、飛行模態改變，甚至舵面操縱困難。現有研究工作通常是把結冰之后的空氣動力學和飛行力學解耦，即是基于氣動力數據庫的飛行動力學仿真，飛行過程中的氣動力及力矩系數是由靜動態氣動數據庫插值得到。該方法雖然在氣動模型中加入了動導數項，一定程度上考慮了非定常效應，但本質上仍是一種基于線性疊加原理的準定常方法。該方法無法適用于飛行器處于強非線性、非定常流場的情形，且未深入考慮非定常空氣動力學和飛行動力學之間的相互耦合作用。

通常研究飛行器氣動/運動耦合特性的方法有以下3種：① 虛擬飛行試驗；② 風洞自由飛行試驗；③ 基于CFD的氣動/運動等多學科耦合的一體化數值模擬，即“數值虛擬飛行”技術。前兩種方法通過試驗手段能較為真實地模擬飛行器機動飛行過程，得到更為準確的動態氣動特性，然而試驗技術十分復雜，難度較大。隨著CFD基礎理論和算法的發展，數值虛擬飛行已成為當前CFD研究領域的熱點。2007年，美國國防部投資3.6億美元開展了計算研究和工程采辦工具與環境項目(CREATE)。研究人員通過該項目開發了固定翼飛機虛擬模擬工具Kestrel軟件，實現了固定翼飛機空氣動力學、結構動力學、飛行動力學、發動機推進等多學科的一體化耦合模擬，基于該軟件對不考慮舵偏的F-22飛機做快速拉升及失速過程進行了模擬，通過與飛行數據的對比，發現模擬結果準確度較高，英國的Allan等對某標模在給定舵偏規律下的自由俯仰運動進行了模擬。在非定常氣動力計算和飛行動力學方程耦合求解的研究方面，國防科技大學的楊小亮對飛行器的多自由度耦合搖滾運動進行了模擬；北京航空航天大學的閻超等對返回艙、帶翼導彈的動態特性進行了數值模擬；中國航天空氣動力研究院的楊云軍等基于數值虛擬飛行技術研究了三角翼在滾轉和側滑兩自由度下的耦合運動特征。中國空氣動力研究與發展中心的張來平等系統論述了虛擬飛行中氣動、運動和控制耦合的數值模擬技術，指出動態網格的生成、高精度湍流模擬方法及耦合求解算法是該技術的難點，并給出了典型的氣動/運動/控制耦合計算流程。

目前國內外沒有專門針對結冰條件下飛機氣動/運動耦合特性展開研究。國內已有的研究主要是對飛行器機動飛行過程的模擬，且大多是對飛行器單自由度或雙自由度自激運動的簡單模擬。本文研究了機翼前緣結冰條件下飛機氣動/運動的耦合特性，介紹了松耦合氣動力與飛行動力學耦合求解策略。通過NACA0012翼型強迫俯仰振蕩標準算例驗證了開發的非定常CFD求解程序及耦合求解算法的有效性和準確性。在此基礎上開發了高精度非定常流場求解器與飛行動力學的雙向耦合仿真平臺，對流場控制方程、飛行動力學方程和運動學方程組在時間尺度上耦合推進求解，能夠刻畫結冰增長特性-分離流動非定常特性-飛行動力學特性三者的相互關系，對于研究結冰條件下的飛行安全問題奠定了基礎，具有重要意義。

1 控制方程及其離散格式

1.1 耦合控制方程組

建立描述飛機運動與非定常流場耦合的控制方程組：

(1)

式中：為Navier-Stokes方程中流體的守恒變量;、和分別為直角坐標系下、、這3個方向上的無黏通量；、、分別對應3個方向黏性通量；為動力學方程狀態變量，包含(飛行速度)、(迎角)、(側滑角)、(滾轉角速度)、(俯仰角速度)和(偏航角速度)；為運動學方程的狀態變量，包含(滾轉角)、(俯仰角)、(偏航角)和位移、、。

從式(1)可以看出，非定常Navier-Stokes方程與飛行動力學方程組互相耦合，相互影響，互為輸入。具體來說就是：由求解非定常Navier-Stokes方程計算出飛行器的氣動力及力矩系數，傳遞給飛行動力學方程，再由飛行動力學方程組解算得到姿態、位移等物理量傳遞給非定常Navier-Stokes方程組，作為流場初始狀態的輸入條件，然后更新網格進行解算。

在耦合計算的過程中，飛機所受氣動力及氣動力矩系數可由Navier-Stokes方程計算得到。根據有限體積法的思想，飛行器受到的氣動力/力矩為物面所有微元所受氣動力/力矩的總和，而每個微元所受氣動力為該處應力張量和外法向矢量的點積，所受氣動力矩為矢徑與氣動力的叉積，具體可表示為

(2)

式中：表示一個微元；為外法向單位向量；為表面切應力；為表面壓應力；(，，)為飛機的氣動中心位置。根據式(2)可計算得到飛行器的氣動力及力矩系數，如式(3)所示：

(3)

1.2 離散格式

1.2.1 非定常Navier-Stokes方程

采用Jameson有限體積法對三維可壓縮非定常Navier-Stokes方程進行空間離散，對無黏項(對流和壓力項)的空間離散采用Roe-FDS通量差分分裂格式，對網格界面通量采用三階MUSCL格式，以提高空間分辨率，對黏性項(剪切應力和熱傳導項)的空間離散采用二階中心格式。采用隱式LU-SGS雙時間步法進行時間推進，該方法在經典隱式LU-SGS格式中引入偽時間導數項，借助偽時間方向的子迭代技術(-TS 迭代)，使得時間離散精度達到二階。

1.2.2 飛行動力學方程

飛機六自由度動力學和運動學方程組是典型的非線性常微分方程組，最常用的是顯式的四階Runge-Kutta數值解法。其計算公式可描述為

(4)

式中：Δ為時間步長；，，，表示為

(5)

參考文獻[17]中的方法，對飛行器剛體六自由度動力學方程采用與Navier-Stokes方程相同的時間離散格式，根據線性多步法松弛法迭代求解：

(6)

式中:為六自由度動力學/運動學方程中的變量、；為對應的六自由度動力學/運動學變量的通量(導數)；為松弛因子；、、、和為系數，通過改變系數值可以選取時間推進精度，本文依次取值為1、0、0、1、0。

2 湍流模擬方法

針對飛行器整機級高雷諾數湍流流動的數值模擬，目前應用最多的還是采用湍流模型的RANS方法，該方法模化了所有尺度的湍流結構，僅對流場平均量求解，計算量相對較小，且能夠保證一定的精度。本文采用工程上常用的湍流模型Spalart-Allmaras(SA)一方程湍流模型模擬飛行器的湍流效應。

3 氣動/運動一體化松耦合計算方法

本文采用松耦合氣動力與飛行動力學耦合求解策略進行氣動/運動耦合計算。松耦合是將流體控制方程Navier-Stokes方程和飛機動力學方程組在各自的時間域上獨立求解，交錯時間推進獲得系統耦合解的算法。其中飛機剛體六自由度運動/動力學方程采用顯式推進，在物理時間推進的每一步內，飛行動力學方程與Navier-Stokes方程之間只交互一次數據。具體的算法步驟如下：

1) 將第物理時刻的飛行器位移、姿態等參數(主要是當前狀態下的高度、馬赫數、姿態角(如迎角、側滑角和滾轉角)等)傳遞給流體控制方程。

2) 根據第時刻的位移、姿態等參數更新流場網格和邊界條件設置，求解Navier-Stokes方程并獲得第物理時刻的氣動力和氣動力矩。

3) 根據第物理時刻的氣動力和氣動力矩，飛行動力學方程進行顯式推進求解，獲得第+1物理時刻飛行器位移、姿態等參數，并傳遞給Navier-Stokes方程。

4) 進行第+1物理時刻的Navier-Stokes方程雙時間步的隱式推進求解，獲得第+1物理時刻飛行器的氣動力和氣動力矩，并傳遞給六自由度全量飛行動力學方程。

5) 按照上述步驟繼續進行迭代，直到達到總的物理時間，仿真結束。

松耦合計算過程中的參數傳遞過程如圖1所示。可以看到，松耦合的求解策略很大程度上保留了CFD程序和飛行仿真程序的獨立性和模塊化，只需編寫數據交互的程序即可，避免了聯立方程面臨的求解難度。但由于兩種方程積分時間的不同步，使得兩者之間始終存在一個物理時間步長的延遲，需要選擇足夠小的時間步長才能保持整個耦合求解的穩定性。

圖1 松耦合的參數傳遞與計算流程Fig.1 Parameter transfer and calculation flow of loose coupling

4 算例驗證

為了考核CFD計算程序模擬動姿態非定常流動的能力，以及剛體動力學/運動學方程與流動控制方程耦合求解算法的有效性，采用二維NACA0012翼型強迫俯仰振蕩的標準算例來進行驗證。

4.1 強迫俯仰振蕩動力學方程

如圖2所示，在亞聲速來流情況下模擬二維NACA0012翼型的強迫俯仰振蕩運動。常規強迫俯仰振蕩運動采用解析形式，定義為翼型迎角隨時間變化的正弦函數：

()=+sin()=+sin(2)

(7)

圖2 NACA0012翼型強迫俯仰振蕩示意圖Fig.2 Forced pitching oscillation diagram of NACA0012 airfoil

在本算例中，有黏、可壓縮自由來流的馬赫數為=0.755，雷諾數為=5.5×10。設定翼型俯仰振蕩的軸心處于距翼型前緣點0.25倍弦長的弦線上，翼型在=0016°初始迎角自由來流作用下開始繞俯仰軸做強迫俯仰振蕩運動，俯仰振蕩振幅=251°。

為了考核剛體動力學/運動學方程與流動控制方程耦合求解方法的有效性，本文將強迫俯仰振蕩運動采用等價的微分方程形式給出：

(8)

通過聯立流動控制方程與式(6)即可實現俯仰振蕩計算。設置物理時間步長為0.001 s，非定常子迭代步數設為100步。

4.2 計算條件設置

如圖3所示，為NACA0012翼型的C型拓撲網格，規模為500×411(流向×法向)，物面第一層網格節點距壁面的無量綱距離滿足<1。遠場邊界距翼型為20倍的弦長，弦長為1 m。邊界條件設置：計算域采用遠場自由來流條件，翼型物面采用無滑移壁面邊界條件。采用數值模擬手段

圖3 NACA0012翼型C型網格Fig.3 NACA0012 airfoil C-grid

來模擬這一非定常流動現象，通過與實驗數據的對比來驗證耦合算法和計算程序模擬非定常流動的準確性。

4.3 計算結果分析

為了提高計算效率，這里先在初始迎角的自由來流條件下計算定常流動，并作為非定常計算的穩定初始流場。計算得到的穩定初場壓力等值線分布如圖4所示。從圖4(a)中可以看到，由于初始迎角較小，同時NACA0012翼型又是對稱翼型，壓力分布大致呈現上下對稱。這與圖4(b)中文獻[13]的結果較為吻合，說明了本文CFD程序具備精確計算定常流動的能力。

圖4 NACA0012翼型壓力等值線分布對比Fig.4 Comparison of pressure isoline distribution of NACA0012 airfoil

在得到穩定的初始流場后，開始進行非定常的計算。采用松耦合算法進行了耦合計算，將計算得到的升力系數及俯仰力矩系數隨迎角的變化曲線與AGARD的實驗數據及Batina的計算結果進行了比較，如圖5所示。圖中：為升力系數；為俯仰力矩系數。

從圖5中可以看到，升力系數的數值模擬結果與Batina的計算結果吻合度較高，與實驗值也較為吻合，只存在較小的偏差；俯仰力矩系數的計算結果總體上與實驗基本吻合，存在小量偏差的原因可能是：翼型繞0.25弦長點進行俯仰振蕩對流場模擬的精度較為敏感，力矩特性的數值模擬難度較大，因此計算結果吻合度不如升力系數吻合度高。

圖5 NACA0012翼型強迫俯仰振蕩過程中的升力系數和俯仰力矩系數Fig.5 Lift coefficient and pitch moment coefficient of NACA0012 airfoil during forced pitching oscillation

上述案例驗證了本文的CFD數值模擬程序能夠較為精確的模擬非定常流動，同時說明了本文提出的耦合算法是可行的。

5 機翼前緣結冰條件下飛機的氣動/運動耦合特性

由于結冰動態增長的時間尺度、Navier-Stokes方程解算的時間尺度和飛行動力學方程解算的時間步長各不相同，且差距較大，因此選取合適的物理時間步長至關重要。由于冰形增長網格重構的復雜性及背景飛機整機帶冰構型巨量的網格將對計算效率造成不可忽視的影響，為了有效解決該問題，做出以下假設：

1) 由于冰的增長速率與形狀會受到外部氣象環境參數、飛行速度、飛行姿態等時變因素的影響，需要在每一物理時間步長進行積冰冰形計算。這里簡化了該項步驟，根據NASA Lewis中心IRT(Icing Research Tunnel)冰風洞公開實驗數據直接給定不同時刻的冰形，同時不考慮冰形沿飛機機翼展向的變化特征。

2) 根據NASA Lewis中心IRT冰風洞公開實驗數據和本文結冰實驗數據表明：冰形的變化速率雖然受到不同的來流條件和飛行參數的影響而發生變化，但總體來看，冰形的增長是緩慢的，通常長達幾十分鐘才能有較為明顯的變化特征。因此，為了兼顧計算效率，在計算機翼積冰動態增長條件下氣動/運動耦合特性時，每20 s更新一次結冰冰形網格，在20 s周期內假定冰形不發生變化。

3) 本文著重研究機翼前緣動態結冰條件下的氣動/運動耦合特性，這里設定平尾前緣結冰冰形保持不變，不發生增長，發動機吊艙唇口處和垂尾前緣無結冰。

5.1 機翼帶冰構型全機三維數模及結構化網格生成方法

結合NASA Lewis中心IRT冰風洞公開數據，形成了一組背景飛機在不同時刻的機翼二維結冰冰形的輸入如圖6所示。根據每一時刻的二維冰形，基于該數模能夠快速建立沿機翼展向的三維冰模如圖7所示。機翼帶冰構型三維數模能夠反映不同時刻二維冰形的幾何細節特征。建立了適合背景飛機帶冰構型氣動特性-飛行力學特性在線耦合分析的機翼帶冰構型計算網格拓撲如圖8所示。本文采用網格生成工具ICEM生成點對接多塊網格，針對不同時刻的機翼帶冰構型，能夠在同一套網格拓撲下生成相應的全機三維結構化網格。

圖6 機翼帶冰構型數模及不同時刻二維冰形輸入Fig.6 Numerical model of wing icing configuration and two-dimensional ice shape input at different time

圖7 基于二維冰形的機翼帶冰三維數模Fig.7 Three-dimensional numerical simulation of wing with ice based on two-dimensional ice shape

由于網格質量直接影響數值模擬結果的好壞，這里充分考慮到冰增長帶來的網格重構的復雜性，在生成三維結構化網格時均進行質量檢查，對于幾何外形變化較大的區域進行局部加密處理，以保證數值模擬計算的精度。

圖8 機翼結冰構型三維結構化計算網格拓撲與三維結構網格的生成Fig.8 Grid topology and generation of three-dimensional structural grid for wing icing configuration

5.2 非定常氣動力求解器與飛行動力學耦合仿真平臺的開發

結合多塊結構化網格高質量計算網格生成技術，開發了求解非定常Navier-Stokes方程的CFD程序。該程序具備多核并行計算、計算條件輸入和流場后處理的功能，并能實現與飛行動力學仿真程序的數據交互。

圖9給出了進行耦合計算的流程。在非定常氣動力計算模塊中，基于上一時間段的飛行力學數據輸入，將ICEM提前生成好的當前時刻結冰冰形的多塊結構化網格輸入到NASA著名的開源計算流體力學軟件CFL3D進行流場計算，得到本時間段的氣動力數據傳遞給飛行動力學特性計算分析模塊進行飛行仿真計算。進行精確的耦合分析時，時間段為每一物理時間步長，通常設置在10量級。如果兼顧計算效率，在結冰緩慢增長帶來氣動力變化幅度較小的前提下，時間段可定義為結冰增長的時間尺度，本文設為20 s，即每20 s更新一次結冰冰形的網格，在此期間與飛行力學交互時冰形的網格不做更新，只改變來流條件。

圖9 結冰條件下背景飛機氣動/運動耦合計算流程Fig.9 Aerodynamic and flight dynamics coupling calculation process of background aircraft under icing condition

本文在進行流場計算時，首先基于RANS方法獲得充分發展的初始流場，在穩定初場的基礎上進行后續非定常計算至非定常流場基本穩定，在進行動態結冰條件下的流場計算時，流場在新結冰外形下從本時間段的穩定流場開始計算至非定常趨勢明顯時結束(具體根據氣動力系數在穩定值上下波動的幅值)，為該時間段的非定常氣動力計算全過程，得到較為精確的大飛機結冰復雜外形下流動的非定常效應。然后進行耦合交互計算，實現非定常氣動力與飛行動力學的耦合分析。

5.3 算例分析

5.3.1 三維整機固定冰形條件下的氣動/運動耦合特征

對10 min結冰構型(中度結冰)下的背景飛機氣動/運動耦合特性進行了研究。背景飛機10 min結冰構型下的初始飛行狀態設定為：以高度6 000 m， 飛行速度為150 m/s(=0.473 5)定直平飛，飛機空重72 000 kg。假設初始時刻飛機迎角和舵偏角為未結冰時該飛行狀態的配平量，配平迎角為6.26°，升降舵偏角為=-6.71°，油門偏度為=34.91%。駕駛員在=0.2 s時感知飛機結冰帶來的低頭力矩，操縱舵面使飛機抬頭，此時=-14.11°。同時本仿真條件中不考慮非對稱結冰情形，通過耦合計算分析在2 s內飛行器姿態及流場的非定常變化特征。設置物理時間步長為0.01 s，進行非定常計算時子迭代步數設為25步，采用松耦合計算策略，因此在2 s的總物理時間域內耦合交互次數為200次。由于在此期間冰形不發生變化，因此不考慮網格的重構，只需改變來流條件。

圖10給出了通過耦合計算得到的該狀態條件下飛行姿態的變化，并對迎角響應過程選取了6個時刻點進行流場分析。可以看到主要是短周期參數迎角和俯仰角速度變化較為劇烈，飛行速度和高度雖然總體變化不大，但均呈現下降趨勢，且下降率越來越大，這是由于迎角超過失速迎角后升力系數陡降導致的。此時飛機處于失速狀態，如果駕駛員不能及時改出，可以預見，飛行高度下降率會急劇增大，飛機將在較短時間內損失較多的高度，對飛行安全產生威脅。圖11為耦合計算過程中升力系數、阻力系數及俯仰力矩系數的變化曲線。可以看到，隨著迎角不斷增大，氣動力及力矩系數的脈動越來越大，非定常效應越來越顯著，導致俯仰角速度發生振蕩。當迎角超過12°左右時，由于失速，升力系數逐漸降低。

圖10 固定冰形條件下耦合計算過程中飛行參數 變化曲線Fig.10 Change curves of flight parameters in coupling calculation under fixed ice shape conditions

圖11 固定冰形條件下耦合計算過程中氣動力及力矩系數變化曲線Fig.11 Change curves of aerodynamic force and moment coefficients in coupling calculation under fixed ice shape conditions

圖12給出了背景飛機在6個時刻下的翼面分離流動形態(左圖)及空間流場結構(右圖)。在右圖中截取了位于機翼9個不同位置的橫截面，畫出了不同時刻流向速度的流線分布，而整機是壓力系數的分布云圖。由圖可以看出，在10 min結冰構型、中等迎角狀態(迎角為6.26°時)時，內外翼分離同時出現，始發位置均位于機翼前緣，平尾處流動無分離。表明結冰使得分離流動始發迎角大幅提前，在中等迎角狀態下就表現出一定的翼面分離流動趨勢。隨著迎角增加，在迎角為8°時就已經出現了明顯的翼面分離流動。在迎角為12°時，空間流動形態發生較大變化，流動分離效應較為顯著，內翼和外翼幾個截面處的分離泡完全破碎，再附效應消失，產生大范圍的強回流區域。平尾處開始發生分離。隨著迎角的繼續增大，機翼處流動分離效應則更為明顯，而平尾處分離流動形態未發生較大變化。

圖12 固定冰形條件下耦合計算過程中不同時刻飛機翼面分離流動形態及空間流場結構Fig.12 Separation flow pattern and spatial flow field structure of aircraft airfoil at different time in coupling calculation process under fixed ice shape conditions

圖13給出了采用準則計算得到的在=0 s迎角為6.26°及=1.98 s、迎角為16°時流動的渦結構。可以看到，在大迎角情形下機翼及平尾處的分離渦較為明顯，而小迎角條件下未發現明顯的分離渦。由于是采用非定常RANS方法和S-A模型對流場進行求解，該方法對強非定常渦運動的模擬無法達到較高精度。因此計算得到的分離渦只能作定性的分析，下一步需要精細化分析流動中渦脫落與非定常特征還需采用RANS/LES混合方法等。

圖13 固定冰形條件下耦合計算過程中不同迎角下的分離渦結構Fig.13 Separated vortex structure at different angles of attack in coupled calculation under fixed ice shape conditions

5.3.2 三維整機動態結冰條件下的氣動/運動耦合特征

采用相同的方法對整機動態結冰條件下的氣動/運動耦合特征進行分析。在該模擬情形下，冰形是不斷增長變化的。仿真模擬的條件設置為：初始飛行高度6 000 m，飛行速度150 m/s，飛機空重42 000 kg，初始航跡俯仰角為-2.5°，并做直線下滑。1 min后穿越含大量過冷水滴的云層，開始逐漸積冰。在整個飛行過程中舵偏量及油門開度保持為初始配平值。假定飛機持續結冰5 min，而后保持不變。總的物理時間設定為360 s。

首先根據設置的條件計算配平量，配平迎角為2.98°，配平升降舵偏角為-4.59°，配平油門偏角為12.41(大小范圍為0～100)，側滑角、副翼偏角及方向舵偏角均為0°。由于在飛行過程中飛機逐漸結冰，導致飛行姿態發生變化。且結冰增長的時間較長，為了減小耦合計算的工作量，采用變步長的思想：在前1 min未結冰期間內，時間步長設為0.1 s，且不進行非定常氣動力的解算，只解算定常的氣動力，這是由于迎角較小且飛機未結冰，非定常效應較為微弱可忽略不計；在5 min結冰增長期間，視結冰氣動力系數變化量的大小來改變時間步長，當飛行姿態改變較小導致計算得到的氣動力變化不大時，時間步長可取為0.1 s，變化較大時可取為0.02 s，以實現減小計算量的同時盡量不影響計算精度。

圖14為通過耦合計算得到的飛行參數的變化曲線。圖15為計算得到的升力系數、阻力系數及俯仰力矩系數變化曲線。由圖可以看出，在動態結冰的過程中，飛機的阻力系數不斷增大，俯仰力矩系數不斷減小。在配平狀態下，氣動力產生的俯仰力矩與發動機推力產生的俯仰力矩平衡；在動態結冰狀態時，由于結冰導致氣動力產生的俯仰力矩出現變化，與發動機推力產生的俯仰力矩不能平衡，導致迎角不斷增大。而升力系數隨著迎角的增大而增大，直到飛機迎角超過失速迎角后，升力系數逐漸降低。隨著迎角的增大，非定常效應越來越顯著，氣動力及力矩系數脈動現象逐漸加強。由于總的俯仰力矩系數存在上下振蕩，導致俯仰角速度的波動較為明顯，從而使飛機的俯仰姿態也出現了一定的振蕩。在整個飛行過程中控制舵面和油門開度沒有變化，保持在初始配平位置，但阻力的增大，高度和速度整體呈現下降趨勢。

圖14 動態結冰條件下耦合計算過程中 飛行參數變化曲線Fig.14 Change curves of flight parameters in coupling calculation under dynamic icing conditions

圖15 動態結冰條件下耦合計算過程中氣動力及力矩系數變化曲線Fig.15 Change curves of aerodynamic force and moment coefficient in coupling calculation under dynamic icing conditions

圖16給出了不同時刻飛機翼面分離流動形態及空間流場結構特征。由圖可以看出，在未結冰、小迎角情形下飛機翼面流動較為平穩；隨著結冰的增長，2 min結冰構型在迎角為5.62°時，翼面沿展向已出現分離特征，從空間流場結構來看未見較大的分離流動；3 min結冰構型的流動情況與之類似，在內翼面后側出現一定的回流區域，同時外翼的展向橫流效應均更加明顯；當冰形增長到4 min結冰構型時，在迎角為7.21°時內翼面后側出現了顯著的分離，從空間流場分布來看，內翼后緣出現再附效應，并顯示了較高的回流強度，從表面壓力分布來看，分離區翼面壓力增加，形成了分離泡結構。隨著迎角進一步增加和冰形的增長，在5 min結冰構型時，分離泡體積迅速膨脹，導致內翼面形成了大范圍回流區域。總體來看，外翼處的分離流動趨勢要顯著弱于內翼面，這是由于機翼/短艙等部件之間的氣動干擾削弱了分離泡的弦向擴展過程，同時抑制了橫向流動，降低了當地流動分離的趨勢。

圖16 動態結冰條件下耦合計算過程中不同時刻飛機翼面分離流動形態及空間流場結構Fig.16 Separation flow pattern and spatial flow field structure of aircraft wing at different time in coupling calculation process under dynamic icing conditions

6 結 論

1) 本文針對機翼前緣結冰條件下空氣動力學與飛行力學耦合特性進行了分析研究。通過對三維整機固定冰形、動態結冰條件下的氣動/運動耦合特性分析，算例計算結果表明：結冰條件下非定常效應的影響在中等迎角階段就開始顯著，隨著結冰的動態增長和迎角的增大，內翼分離流動趨勢進一步增強，表現為分離泡體積的迅速膨脹。由于中等迎角條件下顯著的非定常效應將導致飛機出現一定程度的振蕩，影響飛行穩定性，此時駕駛員應該及時介入，解除危險的飛行狀態。

2) 考慮到計算量的限制，本文沒有考慮冰形沿機翼展向的變化特征。同時在進行動態結冰的流場計算時，網格在新結冰外形下是重新生成的。在后續的研究工作中需要減小這兩處的模擬誤差。