機動飛行數值模擬關鍵技術及初步驗證

李 立

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

機動飛行數值模擬關鍵技術及初步驗證

李 立

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

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

機動飛行; 動導數; 慣性系; 非慣性系; 計算流體力學

0 引言

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

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

1 總體技術思路

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

2 關鍵技術及方案

2.1 機體運動學方程的建立

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

在自由運動形式下,飛行器的飛行狀態隨時間的變化關系可由六自由度方程給出[2]。描述質心位移的平動方程和飛行姿態變化的轉動方程分別為:

(1)

(2)

其中:

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

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

2.2 流動控制方程的建立

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

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

(4)

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

(5)

(6)

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

當網格運動速度為0時,上述方程就是常規的Euler/Navier-Stokes方程。

2.3 運動學方程與流動控制方程的耦合

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

(7)

在慣性系下,網格運動速度根據

可得:

(8)

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

3 具體方案及驗證

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

3.1 運動學方程求解驗證

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

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

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

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

3.2 強迫運動及動導數算例計算

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

(9)

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

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

由圖2可以看出,采用解析形式與本文采用微分方程形式得到的結果一致,表明本文運動學方程與流動控制方程耦合求解的程序實現完全正確。

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

圖3 HBS靜導數和動導數算例計算結果Fig.3 Results for static and dynamic derivatives of HBS

3.3 自由運動算例計算

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

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

4 結束語

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

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

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

[3] 鄧建中,葛仁杰,程正興.計算方法[M].西安:西安交通大學出版社,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)；航空科學基金資助(2015ZA31002)

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

V211.3

A

1002-0853(2017)01-0089-04