結(jié)冰條件下飛機(jī)氣動/運(yùn)動耦合特性

伍強(qiáng)，徐浩軍，魏揚(yáng)，裴彬彬, 3，薛源

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

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

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

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

1 控制方程及其離散格式

1.1 耦合控制方程組

建立描述飛機(jī)運(yùn)動與非定常流場耦合的控制方程組：

(1)

式中：為Navier-Stokes方程中流體的守恒變量;、和分別為直角坐標(biāo)系下、、這3個方向上的無黏通量；、、分別對應(yīng)3個方向黏性通量；為動力學(xué)方程狀態(tài)變量，包含(飛行速度)、(迎角)、(側(cè)滑角)、(滾轉(zhuǎn)角速度)、(俯仰角速度)和(偏航角速度)；為運(yùn)動學(xué)方程的狀態(tài)變量，包含(滾轉(zhuǎn)角)、(俯仰角)、(偏航角)和位移、、。

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

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

(2)

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

(3)

1.2 離散格式

1.2.1 非定常Navier-Stokes方程

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

1.2.2 飛行動力學(xué)方程

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

(4)

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

(5)

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

(6)

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

2 湍流模擬方法

針對飛行器整機(jī)級高雷諾數(shù)湍流流動的數(shù)值模擬，目前應(yīng)用最多的還是采用湍流模型的RANS方法，該方法?；怂谐叨鹊耐牧鹘Y(jié)構(gòu)，僅對流場平均量求解，計算量相對較小，且能夠保證一定的精度。本文采用工程上常用的湍流模型Spalart-Allmaras(SA)一方程湍流模型模擬飛行器的湍流效應(yīng)。

3 氣動/運(yùn)動一體化松耦合計算方法

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

1) 將第物理時刻的飛行器位移、姿態(tài)等參數(shù)(主要是當(dāng)前狀態(tài)下的高度、馬赫數(shù)、姿態(tài)角(如迎角、側(cè)滑角和滾轉(zhuǎn)角)等)傳遞給流體控制方程。

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

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

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

5) 按照上述步驟繼續(xù)進(jìn)行迭代，直到達(dá)到總的物理時間，仿真結(jié)束。

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

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

4 算例驗證

為了考核CFD計算程序模擬動姿態(tài)非定常流動的能力，以及剛體動力學(xué)/運(yùn)動學(xué)方程與流動控制方程耦合求解算法的有效性，采用二維NACA0012翼型強(qiáng)迫俯仰振蕩的標(biāo)準(zhǔn)算例來進(jìn)行驗證。

4.1 強(qiáng)迫俯仰振蕩動力學(xué)方程

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

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

(7)

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

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

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

(8)

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

4.2 計算條件設(shè)置

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

圖3 NACA0012翼型C型網(wǎng)格Fig.3 NACA0012 airfoil C-grid

來模擬這一非定常流動現(xiàn)象，通過與實驗數(shù)據(jù)的對比來驗證耦合算法和計算程序模擬非定常流動的準(zhǔn)確性。

4.3 計算結(jié)果分析

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

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

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

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

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

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

5 機(jī)翼前緣結(jié)冰條件下飛機(jī)的氣動/運(yùn)動耦合特性

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

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

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

3) 本文著重研究機(jī)翼前緣動態(tài)結(jié)冰條件下的氣動/運(yùn)動耦合特性，這里設(shè)定平尾前緣結(jié)冰冰形保持不變，不發(fā)生增長，發(fā)動機(jī)吊艙唇口處和垂尾前緣無結(jié)冰。

5.1 機(jī)翼帶冰構(gòu)型全機(jī)三維數(shù)模及結(jié)構(gòu)化網(wǎng)格生成方法

結(jié)合NASA Lewis中心IRT冰風(fēng)洞公開數(shù)據(jù)，形成了一組背景飛機(jī)在不同時刻的機(jī)翼二維結(jié)冰冰形的輸入如圖6所示。根據(jù)每一時刻的二維冰形，基于該數(shù)模能夠快速建立沿機(jī)翼展向的三維冰模如圖7所示。機(jī)翼帶冰構(gòu)型三維數(shù)模能夠反映不同時刻二維冰形的幾何細(xì)節(jié)特征。建立了適合背景飛機(jī)帶冰構(gòu)型氣動特性-飛行力學(xué)特性在線耦合分析的機(jī)翼帶冰構(gòu)型計算網(wǎng)格拓?fù)淙鐖D8所示。本文采用網(wǎng)格生成工具ICEM生成點對接多塊網(wǎng)格，針對不同時刻的機(jī)翼帶冰構(gòu)型，能夠在同一套網(wǎng)格拓?fù)湎律上鄳?yīng)的全機(jī)三維結(jié)構(gòu)化網(wǎng)格。

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

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

由于網(wǎng)格質(zhì)量直接影響數(shù)值模擬結(jié)果的好壞，這里充分考慮到冰增長帶來的網(wǎng)格重構(gòu)的復(fù)雜性，在生成三維結(jié)構(gòu)化網(wǎng)格時均進(jìn)行質(zhì)量檢查，對于幾何外形變化較大的區(qū)域進(jìn)行局部加密處理，以保證數(shù)值模擬計算的精度。

圖8 機(jī)翼結(jié)冰構(gòu)型三維結(jié)構(gòu)化計算網(wǎng)格拓?fù)渑c三維結(jié)構(gòu)網(wǎng)格的生成Fig.8 Grid topology and generation of three-dimensional structural grid for wing icing configuration

5.2 非定常氣動力求解器與飛行動力學(xué)耦合仿真平臺的開發(fā)

結(jié)合多塊結(jié)構(gòu)化網(wǎng)格高質(zhì)量計算網(wǎng)格生成技術(shù)，開發(fā)了求解非定常Navier-Stokes方程的CFD程序。該程序具備多核并行計算、計算條件輸入和流場后處理的功能，并能實現(xiàn)與飛行動力學(xué)仿真程序的數(shù)據(jù)交互。

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

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

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

5.3 算例分析

5.3.1 三維整機(jī)固定冰形條件下的氣動/運(yùn)動耦合特征

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

圖10給出了通過耦合計算得到的該狀態(tài)條件下飛行姿態(tài)的變化，并對迎角響應(yīng)過程選取了6個時刻點進(jìn)行流場分析?？梢钥吹街饕嵌讨芷趨?shù)迎角和俯仰角速度變化較為劇烈，飛行速度和高度雖然總體變化不大，但均呈現(xiàn)下降趨勢，且下降率越來越大，這是由于迎角超過失速迎角后升力系數(shù)陡降導(dǎo)致的。此時飛機(jī)處于失速狀態(tài)，如果駕駛員不能及時改出，可以預(yù)見，飛行高度下降率會急劇增大，飛機(jī)將在較短時間內(nèi)損失較多的高度，對飛行安全產(chǎn)生威脅。圖11為耦合計算過程中升力系數(shù)、阻力系數(shù)及俯仰力矩系數(shù)的變化曲線?？梢钥吹剑S著迎角不斷增大，氣動力及力矩系數(shù)的脈動越來越大，非定常效應(yīng)越來越顯著，導(dǎo)致俯仰角速度發(fā)生振蕩。當(dāng)迎角超過12°左右時，由于失速，升力系數(shù)逐漸降低。

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

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

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

圖12 固定冰形條件下耦合計算過程中不同時刻飛機(jī)翼面分離流動形態(tài)及空間流場結(jié)構(gòu)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給出了采用準(zhǔn)則計算得到的在=0 s迎角為6.26°及=1.98 s、迎角為16°時流動的渦結(jié)構(gòu)?？梢钥吹剑诖笥乔樾蜗聶C(jī)翼及平尾處的分離渦較為明顯，而小迎角條件下未發(fā)現(xiàn)明顯的分離渦。由于是采用非定常RANS方法和S-A模型對流場進(jìn)行求解，該方法對強(qiáng)非定常渦運(yùn)動的模擬無法達(dá)到較高精度。因此計算得到的分離渦只能作定性的分析，下一步需要精細(xì)化分析流動中渦脫落與非定常特征還需采用RANS/LES混合方法等。

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

5.3.2 三維整機(jī)動態(tài)結(jié)冰條件下的氣動/運(yùn)動耦合特征

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

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

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

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

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

圖16給出了不同時刻飛機(jī)翼面分離流動形態(tài)及空間流場結(jié)構(gòu)特征。由圖可以看出，在未結(jié)冰、小迎角情形下飛機(jī)翼面流動較為平穩(wěn)；隨著結(jié)冰的增長，2 min結(jié)冰構(gòu)型在迎角為5.62°時，翼面沿展向已出現(xiàn)分離特征，從空間流場結(jié)構(gòu)來看未見較大的分離流動；3 min結(jié)冰構(gòu)型的流動情況與之類似，在內(nèi)翼面后側(cè)出現(xiàn)一定的回流區(qū)域，同時外翼的展向橫流效應(yīng)均更加明顯；當(dāng)冰形增長到4 min結(jié)冰構(gòu)型時，在迎角為7.21°時內(nèi)翼面后側(cè)出現(xiàn)了顯著的分離，從空間流場分布來看，內(nèi)翼后緣出現(xiàn)再附效應(yīng)，并顯示了較高的回流強(qiáng)度，從表面壓力分布來看，分離區(qū)翼面壓力增加，形成了分離泡結(jié)構(gòu)。隨著迎角進(jìn)一步增加和冰形的增長，在5 min結(jié)冰構(gòu)型時，分離泡體積迅速膨脹，導(dǎo)致內(nèi)翼面形成了大范圍回流區(qū)域?？傮w來看，外翼處的分離流動趨勢要顯著弱于內(nèi)翼面，這是由于機(jī)翼/短艙等部件之間的氣動干擾削弱了分離泡的弦向擴(kuò)展過程，同時抑制了橫向流動，降低了當(dāng)?shù)亓鲃臃蛛x的趨勢。

圖16 動態(tài)結(jié)冰條件下耦合計算過程中不同時刻飛機(jī)翼面分離流動形態(tài)及空間流場結(jié)構(gòu)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 結(jié) 論

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

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