李 立

(西安航空計算技術(shù)研究所 第七研究室, 陜西 西安 710068)

機動飛行數(shù)值模擬關(guān)鍵技術(shù)及初步驗證

李 立

(西安航空計算技術(shù)研究所 第七研究室, 陜西 西安 710068)

建立飛機機動飛行仿真及飛行性能評估能力是未來航空CFD發(fā)展的重要目標之一,為此,對機動飛行數(shù)值模擬的方法、原理、關(guān)鍵技術(shù)及潛在的應(yīng)用方向進行了探討。采用慣性系下任意拉格朗日-歐拉 (ALE)形式描述的Navier-Stokes方程與機體運動學(xué)方程相耦合的方法,結(jié)合動網(wǎng)格技術(shù)建立了通用的機動飛行數(shù)值模擬應(yīng)用框架。從微分方程形式的運動學(xué)方程出發(fā),建立了與流動控制方程進行緊耦合和松耦合計算的兩套計算方案,并提出了利用預(yù)報-校正技術(shù)提高耦合計算精度的方法。計算結(jié)果表明,所提方案合理、可行,為單自由度/多自由度機動飛行模擬、動導(dǎo)數(shù)分析等計算提供了可靠的新工具。

機動飛行; 動導(dǎo)數(shù); 慣性系; 非慣性系; 計算流體力學(xué)

0 引言

機動性是評價現(xiàn)代航空武器裝備性能優(yōu)劣的重要指標。如第四代戰(zhàn)機設(shè)計明確要求滿足高機動飛行要求,而第五代戰(zhàn)機設(shè)計要求進一步提升以滿足超機動飛行要求[1]。這表明,對飛行器機動飛行及其氣動特性的數(shù)值模擬將成為未來航空CFD應(yīng)用的重要方向。廣義的機動飛行可以定義為飛行狀態(tài)(包括飛行速度、高度和方向)隨時間變化的飛行,比如飛行表演中殲擊機做出的各種機動飛行動作(盤旋、橫滾、“眼睛蛇”等)?？梢?機動飛行數(shù)值模擬顯然是非定常的,必須采用基于非定常計算的手段才可能得到相對準確的計算結(jié)果。這也給算法的構(gòu)建及具體計算帶來實際困難,包括時間求解方面要有足夠的精度,并且由于采用非定常計算迭代,計算時間較長。

本文主要對機動飛行數(shù)值模擬的方法、原理、關(guān)鍵技術(shù)及具體解決途徑進行了探討,建立了一種機動飛行數(shù)值模擬方案,并進行了初步驗證,為進一步深入開展相關(guān)問題的數(shù)值分析奠定了基礎(chǔ)。

1 總體技術(shù)思路

按照機動飛行的定義,要對機動飛行進行數(shù)值模擬,必須首先解決飛行狀態(tài)的定義問題。具體包括3個要素:飛行速度、高度和方向。其中,飛行高度決定了計算中使用的當(dāng)?shù)貧饬鳁l件,與飛行器本身的受力狀態(tài)無關(guān),在一定條件下可以忽略;而飛行速度、飛行方向是模擬中需要密切注意的關(guān)鍵量,與飛行器的受力狀態(tài)直接相關(guān)。在飛行器飛行中,主要受力包括飛行器自身的重力、飛行器飛行中產(chǎn)生的氣動力,以及基于飛行控制律施加的外力。這些合力直接決定了飛行器的飛行狀態(tài)。由此不難分析,機動飛行模擬中,飛行狀態(tài)的定義問題最終可歸結(jié)為飛行器的機體運動學(xué)方程的建立問題。其次是機動飛行條件下流動控制方程如何建立的問題,即如何將傳統(tǒng)的、在靜態(tài)情況下建立的Euler或Navier-Stokes方程推廣到能夠處理機體運動的情況。此外,為了實現(xiàn)具體問題的分析求解,還必須解決上述兩類方程如何耦合進行數(shù)值求解的問題。

2 關(guān)鍵技術(shù)及方案

2.1 機體運動學(xué)方程的建立

飛行器飛行狀態(tài)可通過機體運動學(xué)方程來描述:描述飛行速度的方程對應(yīng)了機體質(zhì)心的運動學(xué)方程,即平動方程;描述飛行方向的方程對應(yīng)了機體姿態(tài)變化的運動學(xué)方程,即轉(zhuǎn)動方程。這兩類方程均可通過3種方式建立,即解析方法、離散數(shù)據(jù)方法和動力學(xué)方程方法。其中:前兩種方式本質(zhì)上屬于強迫運動形式,即飛行器飛行狀態(tài)隨時間的變化關(guān)系預(yù)先明確給定,計算中只關(guān)心機動飛行中飛行器氣動特性的變化,不再考慮飛行軌跡的模擬及如何具體實現(xiàn)的問題;第三種方式可看作自由運動形式,即飛行器飛行狀態(tài)隨時間的變化關(guān)系需要以數(shù)值求解的方式,通過與流動控制方程耦合求解得到,同時,計算中需要考慮飛行控制律及飛行中各操縱面的影響。

在自由運動形式下,飛行器的飛行狀態(tài)隨時間的變化關(guān)系可由六自由度方程給出[2]。描述質(zhì)心位移的平動方程和飛行姿態(tài)變化的轉(zhuǎn)動方程分別為:

(1)

(2)

其中:

式中:T為姿態(tài)角變化率與體軸系下角速度之間的轉(zhuǎn)換矩陣;Fouter,Mouter分別為合外力及合外力矩;m,I分別為飛行器機體的質(zhì)量和轉(zhuǎn)動慣量張量。

上述方程根據(jù)實際情況,還可只考慮平動或轉(zhuǎn)動,從而形成三自由度方程。

2.2 流動控制方程的建立

在機動飛行模擬中,流動控制方程的建立可以采用兩種思路。一種是在絕對坐標系(慣性系)下建立方程,這時求解的流體動力學(xué)方程為常規(guī)Euler/Navier-Stokes方程,同時考慮機體運動對流體的影響。其優(yōu)點是流動控制方程的公式表達簡單,可在現(xiàn)有求解器上直接開發(fā),只需要在通量計算和邊界處理中計及網(wǎng)格運動速度影響;缺點是必須顯式進行網(wǎng)格運動,因而必須引入動網(wǎng)格技術(shù)來實現(xiàn)完整的求解過程。另一種思路是在相對坐標系(非慣性系)下建立方程,這時求解的流體力學(xué)方程是關(guān)于流體相對運動速度的方程,機體運動的影響通過控制方程的源項來體現(xiàn)。其優(yōu)點是不需要顯式進行動網(wǎng)格,因而不需要使用動網(wǎng)格技術(shù),所有計算均在同一計算網(wǎng)格下進行;缺點是公式推導(dǎo)非常復(fù)雜,流場(在慣性系下)的動態(tài)變化過程必須通過對計算結(jié)果的后處理來完成。

本文選擇第一種思路。在這種情況下,計算中計算域與機體一起運動,機體運動對流體的影響體現(xiàn)為網(wǎng)格運動帶來的幾何守恒律問題?；趲缀问睾懵?不難推導(dǎo)得到采用任意歐拉-拉格朗日(ALE)形式統(tǒng)一描述的Euler/Navier-Stokes方程:

(4)

其中,守恒變量及對流通量項的形式為:

(5)

(6)

式中:U,at分別為流體的相對速度及網(wǎng)格運動的法向速度,定義U=(u-xt,v-yt,w-zt),at=xtnx+ytny+ztnz。

當(dāng)網(wǎng)格運動速度為0時,上述方程就是常規(guī)的Euler/Navier-Stokes方程。

2.3 運動學(xué)方程與流動控制方程的耦合

由ALE描述的Euler/Navier-Stokes方程(4),采用有限體積法可得其半離散形式為:

(7)

在慣性系下,網(wǎng)格運動速度根據(jù)

可得:

(8)

這樣,式(7)和式(8)連同運動學(xué)方程式(1)和式(2)就構(gòu)成了機動飛行模擬實際需要聯(lián)立求解的方程組。該方程組是一個關(guān)于時間微分的常微分方程組,可采用松耦合或緊耦合的策略進行求解。松耦合方法是在同一時間步內(nèi)獨立求解流體力學(xué)方程和運動學(xué)方程,但按同一物理時間步推進更新流動控制變量及運動學(xué)參數(shù)。緊耦合方法是在同一時間步上,聯(lián)立求解流動控制方程和運動學(xué)方程。理論上,緊耦合方法穩(wěn)定性優(yōu)于松耦合方法,但具體實踐表明,兩種方式的計算效果相當(dāng)。松耦合方法的優(yōu)勢是可獨立形成求解模塊,便于編程及程序移植。本文采用松耦合方法。

3 具體方案及驗證

整體采用松耦合求解思路,在同一物理時間步實現(xiàn)兩者的耦合,利用預(yù)報-校正的技巧提高整體耦合求解的計算精度。流動控制方程求解基于結(jié)構(gòu)網(wǎng)格有限體積法,采用時-空分開的半離散形式,空間離散采用3階迎風(fēng)偏置AUSM+格式(基于MUSCL插值和Albada限制器實現(xiàn)),時間推進采用τ-TS雙時間步方法實現(xiàn)非定常計算。其中,物理時間計算采用二階精度的后向差分法;虛擬時間計算采用LU-SGS格式。運動學(xué)方程求解通過獨立模塊實現(xiàn),基于統(tǒng)一策略采用4級Runge-Kutta方法求解,并采用預(yù)報-校正技巧及“小窗口”技術(shù)提高整體耦合的物理時間計算精度[3]。

3.1 運動學(xué)方程求解驗證

為了對運動學(xué)方程求解模塊進行驗證,采用兩個具有解析解的典型算例進行了計算。一個是自由落體算例,用于驗證平動方程求解的計算精度;另一個是自激振蕩算例,用于驗證轉(zhuǎn)動方程求解的計算精度[4]。算例參數(shù)均進行了無量化處理。

自激振蕩算例中,假設(shè)機體不受任何外力和外力矩的作用,初始速度為u0=(0,0,0)；初始轉(zhuǎn)速為ω0=(1.0,0.0,0.5),對應(yīng)的解析解為ωx=a× cos(λt),ωy=bsin(λt),ωz=c,其中,a=1,b=-1,c=0.5,λ=0.25。

圖1給出了不同時間步長得到的數(shù)值解與解析解的比較?？梢钥闯?本文運動學(xué)方程求解程序穩(wěn)定性非常好,對不同時間步長都有很好的適應(yīng)性。

圖1 運動學(xué)方程求解驗證結(jié)果Fig.1 Validation results for solving the kinetic equation

3.2 強迫運動及動導(dǎo)數(shù)算例計算

與常規(guī)強迫運動采用解析形式或離散數(shù)據(jù)形式給出不同,這里,強迫運動采用等價的微分方程形式給出,以便考核運動學(xué)方程與流動控制方程耦合求解的能力,采用NACA0012俯仰振蕩算例[5]和高超聲速動導(dǎo)數(shù)計算標模HBS算例。俯仰振蕩解析表達式θ=α+Δθsin(2πft)的等價微分形式可寫為:

(9)

這樣,通過聯(lián)立流動控制方程、式(9)和角位移方程即可實現(xiàn)俯仰振蕩計算。圖2給出了NACA0012翼型典型結(jié)果。選取的計算狀態(tài)為:Ma=0.6,α=2.89°,Re=4.8×106,并取式(9)中Δθ=1.41°,f=50.32 Hz,t無量綱。

圖2 NACA0012俯仰振蕩算例計算結(jié)果Fig.2 Typical results for pitching NACA0012 airfoil

由圖2可以看出,采用解析形式與本文采用微分方程形式得到的結(jié)果一致,表明本文運動學(xué)方程與流動控制方程耦合求解的程序?qū)崿F(xiàn)完全正確。

圖3給出了采用類似方法得到的高超聲速彈體外形(HBS)算例的靜導(dǎo)數(shù)和動導(dǎo)數(shù)計算結(jié)果,來流條件為:Ma=6.85,α=0°～15°。其中,動導(dǎo)數(shù)的后處理分析采用積分法、最小二乘法、T-F方法等多種方法作為對比?？梢钥闯?幾種處理方法得到的結(jié)果一致,并與試驗值和理論值符合良好,獲得了預(yù)期的計算效果。

圖3 HBS靜導(dǎo)數(shù)和動導(dǎo)數(shù)算例計算結(jié)果Fig.3 Results for static and dynamic derivatives of HBS

3.3 自由運動算例計算

為了對程序在自由運動情況下的機動飛行模擬能力進行驗證,選擇NACA0012翼型進行自由振蕩計算。這類計算的實際應(yīng)用包括通過計算尋找配平迎角的計算等。具體選擇Ma=0.8,α0=0°和Ma=0.8,α0=1.25°兩個狀態(tài)開展計算,t無量綱。模擬場景為翼型固定在0.25c位置,在氣動力矩作用下逐漸平衡。圖4給出了計算得到的迎角及升力系數(shù)隨時間的變化歷程?？梢钥闯?兩個計算條件下得到的配平狀態(tài)基本一致,結(jié)果合理。

圖4 NACA0012自由振蕩算例計算結(jié)果Fig.4 Calculation results for free oscillating NACA0012 airfoil

4 結(jié)束語

本文對飛行器(飛機、導(dǎo)彈)機動數(shù)值模擬的相關(guān)關(guān)鍵技術(shù)進行了研究,發(fā)展了一套開放的機動飛行數(shù)值模擬計算框架,并通過典型算例進行了程序驗證,得到了合理的計算結(jié)果,表明上述方案合理、可行,有望為單自由度/多自由度機動飛行模擬、動導(dǎo)數(shù)計算提供可靠的新工具。下一步有必要進一步開展上述程序的工程驗證及復(fù)雜工程應(yīng)用。同時,在該程序基礎(chǔ)上,研制在非慣性系下的機動飛行模擬程序,主要難點包括:推導(dǎo)非慣性系下基于相對運動速度描述的流動控制方程的一般形式;等效源項的隱式處理;計算結(jié)果后處理等。

[1] 高勁松,陳哨東.國外六代機發(fā)展情況研究[J].飛航導(dǎo)彈,2014(1):54-59.

[2] 肖業(yè)倫.飛行器運動方程[M].北京:航空工業(yè)出版社,1987:21-24.

[3] 鄧建中,葛仁杰,程正興.計算方法[M].西安:西安交通大學(xué)出版社,1985:269-270.

[4] Murman S M,Aftosmis M J.Simulations of 6-DOF motion with a cartesian method[J].AIAA-2003-1246,2003.

[5] Landon R H.NACA0012 oscillating and transient pitch,compendium of unsteady aerodynamic measurements[R].AGARD-R-702,1982.

(編輯：李怡)

Key technologies and preliminary validation for numerical simulation of maneuver flight

LI Li

(Seventh Department, ACTRI, Xi’an 710068, China)

Establishing the capability of flight simulation and flight performance estimation is one of the important goals for the future aviation CFD, some keynotes on numerical simulation of maneuvering flight with its methodologies, principles, key technologies and potential applications were discussed. Combined with moving grid techniques, an open framework for simulation of maneuvering flight was established based on the Navier-Stokes equation in arbitrary Lagrange-Euler (ALE) formulation coupled with the kinematic equation for aircraft. Based on kinematics equation in the form of differential equation, two numerical schemes to loose coupling and tight coupling with flow control equation were built, and a predictor-corrector method was further proposed to increase the prediction accuracy. Numerical results show that the schemes are reasonable and feasible, which can be used as a new robust tool for simulation of single or multiple freedom degree maneuver, prediction of dynamic derivatives.

maneuvering flight; dynamic derivative; inertial frame; non-inertial frame; CFD

2016-06-08;

2016-10-27;

時間:2016-11-10 09:10

國家863計劃資助項目(2012AA01A304)；航空科學(xué)基金資助(2015ZA31002)

李立(1977-)，四川成都人，高級工程師，碩士，研究方向為計算流體力學(xué)。

V211.3

A

1002-0853(2017)01-0089-04