黃燦 趙永輝

(南京航空航天大學(xué)機(jī)械結(jié)構(gòu)力學(xué)及控制國家重點實驗室，南京 210016)

引言

顫振是飛行器氣動彈性問題中的重要課題.傳統(tǒng)氣彈研究在亞聲速上使用渦格法和偶極子網(wǎng)格法等頻域方法可以得到較為精確的氣動力解.而在跨聲速時，以上方法不再適用.非定常CFD的方法可以進(jìn)行大迎角、跨聲速等非線性計算，在跨聲速領(lǐng)域目前只有CFD方法計算最為精確［1］.然而，CFD計算所耗時間非常長，計算效率十分低下.

針對CFD計算離散化的特點，利用GPU并行計算可以天然地將任務(wù)劃分成多份同時進(jìn)行，大大的提高計算效率.GPU并行計算是近些年在計算領(lǐng)域快速發(fā)展的一項技術(shù).GPU與CPU相比有著更強(qiáng)的單精度浮點計算能力和更寬的內(nèi)存帶寬［2-3］.在浮點運算上，GPU比CPU有著相當(dāng)大的優(yōu)勢，同時期的GPU單精度浮點計算能力可達(dá)CPU的十倍.CUDA語言是GPU上的語言，可以通過CUDA命令調(diào)用GPU資源進(jìn)行計算，實現(xiàn)計算的并行化［4-10］.

近些年系統(tǒng)識別也被應(yīng)用到CFD計算中來，作為提高計算效率的另一種途徑.目前有多種模型降階(ROM，reduced-order model)方法，包括頻域模型、時域自回歸滑動平均模型(ARMA，auto-regressive-moving-average)和離散時間狀態(tài)空間模型等［11-13］.CFD技術(shù)屬于時域方法，本文選擇ARMA模型來進(jìn)行氣動力建模.建立精確的非定常時域氣動力降階模型可以繞過繁雜的CFD流場求解器，直接根據(jù)輸入輸出關(guān)系得出想要的結(jié)果，這樣可以大大降低計算所耗的時間.本文將這兩種提高氣動彈性分析效率的方法結(jié)合到了一起.首先利用GPU進(jìn)行并行CFD計算，得到非定常氣動力.然后利用非定常氣動力識別出降階模型.最后對降階模型的精度進(jìn)行了考核.

1 計算方法

圖1 基于氣動力辨識技術(shù)的氣動彈性仿真流程圖Fig.1 Flow chart of aero-elastic simulation based on the identification of aerodynamic force

基于CFD系統(tǒng)辨識的氣動彈性分析需要的流程如圖1所示.其中非定常求解器求解所耗時間最多，本文將對這一部分進(jìn)行并行化處理以提高效率.另外，對此部分進(jìn)行系統(tǒng)辨識，得到降階的氣動力模型，再次提高計算效率.需要注意的是，不同的初始條件下辨識得到的非定常氣動力降階模型只能使用一次，初始條件改變則需要重新識別.例如，在廣義氣動力的計算中，不同馬赫數(shù)的情況需要對應(yīng)不同的降階模型.

CFD計算部分采用Euler方程，守恒形式的Euler方程如下:

其中x和y為笛卡爾坐標(biāo)系坐標(biāo)，W為守恒變量:

F，G表示通量:

其中，ρ、P、H和E分別表示密度、壓強(qiáng)、單元總焓和單元總能量.U，V分別表示笛卡爾坐標(biāo)系下向x和y向的速度矢量.這些量由理想氣體的單位體積的總能量E和總焓H互相聯(lián)系，式中γ為比熱比.

本文采用Jameson中心格式的有限體積法進(jìn)行空間離散，利用雙時間推進(jìn)法進(jìn)行時間離散.物理時間采用顯示格式，擬時間采用四步龍格-庫塔推進(jìn)格式.網(wǎng)格使用的是非結(jié)構(gòu)網(wǎng)格，動網(wǎng)格技術(shù)使用了彈簧法［14-18］.

氣動彈性控制方程為:

式中M為質(zhì)量矩陣，K為剛度矩陣，F(xiàn)為廣義氣動力向量，q為廣義位移向量，它們分別可以表示為:

其中μ為質(zhì)量比，Cl和Cm分別為升力系數(shù)和力矩系數(shù)，rα為機(jī)翼對剛心的無量綱回轉(zhuǎn)半徑，xα為重心與剛心的無量綱距離，ωα和ωh分別為俯仰和浮沉兩個自由度的固有頻率為無量綱顫振速度.

為了便于時域求解，引入狀態(tài)變量E=(q1，q2，…，qN，˙q1，˙q2，…，˙qN)T，則方程(4)可以寫成狀態(tài)空間的形式:

式中

在每個時間步內(nèi)，將視為時間的單值函數(shù)，這樣應(yīng)用改進(jìn)的龍格-庫塔方法可將方程(5)展開為:

上式后幾步中的F(t+Δ/2)和F(t+Δt)可用前幾個時刻(F(t)，F(xiàn)(t-Δt)，F(xiàn)(t-2Δt))的值插值得到.該方法計算所得的精度要高于傳統(tǒng)的凍結(jié)氣動力的方法，而且在效率上大大超過標(biāo)準(zhǔn)的龍格-庫塔方法.

圖2 GPU任務(wù)劃分原理圖Fig.2 The principle of division of tasks on GPU

GPU并行計算本文應(yīng)用到CUDA語言和NVIDIA的GPU.并行計算的原理是將程序中的循環(huán)部分拆開，分別交給k(循環(huán)次數(shù))個計算單元同時進(jìn)行計算，再將計算結(jié)果讀回主機(jī)端.GPU做并行計算有著天然的優(yōu)勢，其最小計算單元是線程，GPU中的線程數(shù)遠(yuǎn)遠(yuǎn)大于CPU中的核心數(shù).其工作原理如圖2:

并行計算要求各線程塊和線程之間的計算互不干擾，于是隨時間變化的迭代過程是不能夠并行化的，因此只有適合的離散模型才能用GPU進(jìn)行并行計算.非定常Euler方程中，多處進(jìn)行了按單元或者按邊循環(huán)的計算，這種空間離散模型在循環(huán)計算時同一時間步內(nèi)邊與邊的計算相互獨立，可以用GPU進(jìn)行并行計算.隨著網(wǎng)格單元數(shù)的上升，CPU計算時間大幅增加，而對于GPU的計算時間卻增加很小.理論上并行部分的計算時間幾乎不變，只是網(wǎng)格數(shù)的增加會增加一定的數(shù)據(jù)傳輸時間.于是對這些部分進(jìn)行并行是可行的.其中主要包括通量和殘值計算以及動網(wǎng)格計算的并行化.

圖3 CFD計算流程Fig.3 Process of CFD

圖3為CFD的計算流程，其中陰影部分均為可并行部分.本文將算例的CFD計算部分并行化處理，大大提高計算效率.本文計算的硬件環(huán)境為:Intel Core 2 Duo CPU，2.94GHz，4GB內(nèi)存;NVIDIA Geforce GTX550Ti，192個流處理器(SM)，1GB顯存.軟件環(huán)境為Microsoft Windows7 SP1操作系統(tǒng)，VS2008編譯器，CUDA版本為4.2.

系統(tǒng)識別采用ARMA模型.多輸入多輸出的ARMA模型實際上是系統(tǒng)的離散差分模型，其表達(dá)式為:

式中，y(k)為系統(tǒng)輸出量(這里指廣義氣動力系數(shù)向量)的第k次觀測值;y(k-1)為系統(tǒng)輸出量的第k-1次觀測值，以此類推;u(k)為系統(tǒng)第k個輸入量(這里指廣義位移向量);u(k-1)為系統(tǒng)第k-1個輸入量，以此類推;e(k)為零均值的隨機(jī)噪聲;Ai和Bi為待辨識的參數(shù)矩陣;na和nb分別為輸出和輸入的延遲階數(shù).

設(shè)系統(tǒng)共采集了L組數(shù)據(jù)，即k=1，2，…，L，將式(7)化為如下方程:

式中:

根據(jù)最小二乘法，在隨機(jī)噪聲最小時，可得待識別參數(shù)矩陣θ的估計為:

本文選擇了比較容易實現(xiàn)并且具有很寬頻帶寬度的“3211”速度輸入.通過調(diào)整訓(xùn)練信號的時間步長，可以將頻帶移到試驗所希望激發(fā)的頻帶上去.由于輸入輸出數(shù)據(jù)是一次性批處理的，故采用上文所述的一次性最小二乘估計，其精度較高.參數(shù)辨識需要對系統(tǒng)進(jìn)行穩(wěn)態(tài)清零，使系統(tǒng)滿足零輸入-零輸出特性.

2 算例

為了驗證流場求解器的正確性，首先選用NACA0012翼型進(jìn)行非定常氣動力計算，用計算結(jié)果與實驗結(jié)果對比.以翼型的俯仰簡諧振蕩為例:

其中，初始攻角α0=0.016°，角度幅值αm=2.51°，角頻率ω=k·V∞/b，縮減頻率k=0.0814，半弦長b=0.5m，Ma=0.755.計算結(jié)果如圖4所示:

圖4 NACA0012翼型非定常氣動力系數(shù)與實驗值對比Fig.4 NACA0012 airfoil’s unsteady aerodynamic force simulation compared with experiment

圖5 3211信號輸入下Euler解和降階模型解比較Fig.5 The Euler solution compared with the ROM solution on 3211 input signal

顫振計算中本文選擇了Isogai Wing，它是二維跨聲速氣動彈性計算的標(biāo)準(zhǔn)算例，其翼型為NACA64A010翼型.二元翼段氣動彈性系統(tǒng)的物理參數(shù)為:

針對NACA64A010翼型，本文首先使用“3211”輸入信號進(jìn)行激勵，通過調(diào)整輸入信號的帶寬和幅值將信號訓(xùn)練到合適的頻帶上來.訓(xùn)練信號包含了俯仰和浮沉兩個自由度的輸入，所得到的輸出信號包括和兩個系統(tǒng)輸出.

圖5中給出了“3211”信號輸入下的Euler解和辨識模型解的比較.在相同輸入下辨識模型給出的輸出結(jié)果與非定常Euler求解器的輸出結(jié)果吻合很好.可以看出，該辨識模型在很寬的頻帶范圍內(nèi)對原系統(tǒng)都有很好的近似.

圖6 NACA64A010翼型的顫振臨界響應(yīng)(Ma=0.825=0.54)Fig.6 The critical flutter response of NACA64A010 airfoil(Ma=0.825=0.54)

圖7 NACA64A010翼型的顫振速度邊界Fig.7 The flutter velocity boundary of NACA64A010 airfoil

針對Isogai Wing，算例計算了其顫振速度隨馬赫數(shù)變化的邊界.圖7中可以看出，與眾多文獻(xiàn)中的結(jié)果相比，基于非定常Euler方程的顫振速度邊界和辨識模型計算的顫振邊界在的狀態(tài)下吻合的非常好.

CUDA語言的并行程序和串行C++程序計算的加速比約為2.8倍;ROM降階模型與全階CFD模型計算顫振速度加速比約為2.3倍;將并行程序和ROM降階模型結(jié)合與全階串行程序相比加速比約為6.4倍.對于二維翼型的計算，參考文獻(xiàn)［10］中的結(jié)果來看，其應(yīng)用CUDA語言并行計算得到了1.04倍的加速比，本文結(jié)果加速比更優(yōu).

表1 各方法下計算與原始程序計算加速比Table 1 The speedup of the program on different methods

3 結(jié)論

本文結(jié)合了兩種非定常Euler方程的加速方法:基于CFD的非定常氣動力采用GPU并行計算實現(xiàn)加速;之后利用CFD計算結(jié)果，得到了非定常氣動力的ARMA模型.仿真結(jié)果表明ARMA降階模型對原非定常氣動模型還原較好，在保證了一定精度的情況下提升了計算效率.此外，CUDA并行程序受限于網(wǎng)格量并未發(fā)揮全部潛力，但優(yōu)于文獻(xiàn)［10］的二維結(jié)果，在三維問題上GPU并行計算對效率的提升將會非?？捎^.本文將兩者結(jié)合可大大減少基于CFD的氣動彈性計算時間，值得進(jìn)一步深入研究.

1 張偉偉，葉正寅.基于非定常氣動力辨識技術(shù)的氣動彈性數(shù)值模擬.航空學(xué)報，2006，27(4)，579～583(Zhang W W，Ye Z Y.Numerical simulation of aeroelisticity basing on identification technology of unsteady aerodynamic loads.Acta Aeronoautica ET Astronautica Sinica，2006，27(4)，579～583(in Chinese))

2 張加樂.基于GPU并行計算的非定常Euler方程算法研究［碩士學(xué)位論文］.南京:航空航天大學(xué)，2012(Zhang J L.Numerical studies of unsteady euler equations based on GPU parallel computing［Master Thesis］.Nanjing:Nanjing University of Aeronautics and Astronautics，2012(in Chinese))

3 苗樹明.NS方程在GPU上的并行實現(xiàn)［碩士學(xué)位論文］.上海:上海交通大學(xué)，2011(Miao S M.Implementation of NSequations in parallel on GPU［Master Thesis］.Shanghai:Shanghai Jiao Tong University，2011(in Chinese))

4 Corrigan A，Camelli F，Lohner R，Wallin J.Running unstructured grid based CFD solvers in modern graphics hardware.AIAA，2009:4001

5 Antoniou A S，Karantasis K I，Polychronopoulos E D.Acceleration of a finite difference WENO scheme for largescale simulations on many-core architectures.AIAA，2010:525

6 Brandvik T，Pullan G.Acceleration of a 3D euler solver using commodity graphics hardware.AIAA，2008:607

7 Shinn A F，Vanka SP，Hwu WW.Direct numerical simulation of turbulent flow in a square duct using a graphics processing unit(GPU).AIAA，2010:5029

8 Elsen E，LeGresley P，Darve E.Large calculation of the flow over a hypersonic vehicle using a GPU.Journal of Computational Physics，2008，227:10148～10161

9 Hargen T R，Lie K A，Natvig JR.Solving the euler equations on graphic processing units.Computational Science，2006，3994:220～227

10 張兵，韓景龍.基于GPU和隱式格式的CFD并行計算方法.航空學(xué)報，2010，2:249～256(Zhang B，Han J L.Parallel computing methods for CFD using a GPU and implicit scheme.Acta Aeronoautica ET Astronautica Sinica，2010，2:249～256(in Chinese))

11 Daniella E Raveh.Identification of computational-fluiddynamics based unsteady aerodynamic models for aeroelastic analysis.Journal of Aircraft，2004，41(3):620～632

12 Gupta K K，Bach C.Systems identification approach for a computational-fluid-dynamics based aeroelastic analysis.AIAA，2007，45(12):2820～2827

13 Lai K L，Won K S，Koh E P C，Tsai H M.Flutter simulation and prediction with CFD-based reduced-order model.AIAA，2012:2006～2026

14 Stolcis L，Johnston L J.Solution of the euler equations on unstructured grids for two-dimensional compressible flow.Aeronautical Journal，1990:181～195

15 John T Batina.Implicit flux-split euler schemes for unsteady aerodynamic analysis involving unstructured dynamic meshes.AIAA，1991，29(11):1836～1843

16 鄧楓，伍貽兆，劉學(xué)強(qiáng).基于混合動網(wǎng)格的二維非定常粘性流動數(shù)值模擬.南京航空航天大學(xué)學(xué)報，2007，39(4):444～448(Deng F，Wu Y Z，Liu X Q.Numerical simulation of two-dimensional unsteady viscous flow based on hybird dynamic grids.Journal of Nanjing University of Aeronautics＆Astronautics，2007，39(4):444～448(in Chinese))

17 王軍利，白俊強(qiáng)，詹浩.基于非結(jié)構(gòu)動網(wǎng)格的非定常氣動力計算.飛機(jī)設(shè)計2005，9:24～29(Wang J L，Bai J Q，Zhan H.Calculation of unsteady flow using dynamic unsteructured grids.Aircraft Design，2005，9:24～29(in Chinese))

18 Jameson A.Time dependent calculations using multigrid，with applications to unsteady flows past airfoils and wings.10thAIAAComputational Fluid Dynamics Conference，1991:1596