機(jī)動(dòng)飛行數(shù)值模擬關(guān)鍵技術(shù)及初步驗(yàn)證

李 立

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

機(jī)動(dòng)飛行數(shù)值模擬關(guān)鍵技術(shù)及初步驗(yàn)證

李 立

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

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

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

0 引言

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

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

1 總體技術(shù)思路

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

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

2.1 機(jī)體運(yùn)動(dòng)學(xué)方程的建立

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

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

(1)

(2)

其中:

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

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

2.2 流動(dòng)控制方程的建立

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

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

(4)

其中,守恒變量及對流通量項(xiàng)的形式為:

(5)

(6)

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

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

2.3 運(yùn)動(dòng)學(xué)方程與流動(dòng)控制方程的耦合

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

(7)

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

可得:

(8)

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

3 具體方案及驗(yàn)證

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

3.1 運(yùn)動(dòng)學(xué)方程求解驗(yàn)證

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

自激振蕩算例中,假設(shè)機(jī)體不受任何外力和外力矩的作用,初始速度為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í)間步長得到的數(shù)值解與解析解的比較?？梢钥闯?本文運(yùn)動(dòng)學(xué)方程求解程序穩(wěn)定性非常好,對不同時(shí)間步長都有很好的適應(yīng)性。

圖1 運(yùn)動(dòng)學(xué)方程求解驗(yàn)證結(jié)果Fig.1 Validation results for solving the kinetic equation

3.2 強(qiáng)迫運(yùn)動(dòng)及動(dòng)導(dǎo)數(shù)算例計(jì)算

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

(9)

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

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

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

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

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

3.3 自由運(yùn)動(dòng)算例計(jì)算

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

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

4 結(jié)束語

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

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

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

[3] 鄧建中,葛仁杰,程正興.計(jì)算方法[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;

時(shí)間:2016-11-10 09:10

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

李立(1977-)，四川成都人，高級工程師，碩士，研究方向?yàn)橛?jì)算流體力學(xué)。

V211.3

A

1002-0853(2017)01-0089-04