基于CFD系統(tǒng)辨識的氣彈分析及GPU并行算法初探*

黃燦 趙永輝

(南京航空航天大學機械結構力學及控制國家重點實驗室，南京 210016)

引言

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

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

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

1 計算方法

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

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

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

其中x和y為笛卡爾坐標系坐標，W為守恒變量:

F，G表示通量:

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

本文采用Jameson中心格式的有限體積法進行空間離散，利用雙時間推進法進行時間離散.物理時間采用顯示格式，擬時間采用四步龍格-庫塔推進格式.網格使用的是非結構網格，動網格技術使用了彈簧法［14-18］.

氣動彈性控制方程為:

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

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

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

式中

在每個時間步內，將視為時間的單值函數，這樣應用改進的龍格-庫塔方法可將方程(5)展開為:

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

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

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

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

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

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

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

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

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

式中:

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

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

2 算例

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

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

圖4 NACA0012翼型非定常氣動力系數與實驗值對比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，它是二維跨聲速氣動彈性計算的標準算例，其翼型為NACA64A010翼型.二元翼段氣動彈性系統(tǒng)的物理參數為:

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

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

圖6 NACA64A010翼型的顫振臨界響應(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，算例計算了其顫振速度隨馬赫數變化的邊界.圖7中可以看出，與眾多文獻中的結果相比，基于非定常Euler方程的顫振速度邊界和辨識模型計算的顫振邊界在的狀態(tài)下吻合的非常好.

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

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

3 結論

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

1 張偉偉，葉正寅.基于非定常氣動力辨識技術的氣動彈性數值模擬.航空學報，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方程算法研究［碩士學位論文］.南京:航空航天大學，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上的并行實現［碩士學位論文］.上海:上海交通大學，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并行計算方法.航空學報，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 鄧楓，伍貽兆，劉學強.基于混合動網格的二維非定常粘性流動數值模擬.南京航空航天大學學報，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 王軍利，白俊強，詹浩.基于非結構動網格的非定常氣動力計算.飛機設計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