姜悅寧，賈宏光，厲　明

JIANG Yuening1，2，3,JIAHongguang1，3,LI Ming1，3

1.中國科學(xué)院 長春光學(xué)精密機(jī)械與物理研究所，長春 130033

2.中國科學(xué)院大學(xué)，北京 100039

3.長光衛(wèi)星技術(shù)有限公司，長春 130033

1.Changchun Institute of Optics,Fine Mechanics and Physics,ChineseAcademy of Sciences,Changchun 130033,China

2.University of ChineseAcademy of Sciences,Beijing 100039,China

3.Chang Guang Satellite Technology Co.,Ltd.,Changchun 130033,China

1　引言

無人機(jī)以其結(jié)構(gòu)輕便、可操縱性強(qiáng)等獨特的技術(shù)優(yōu)勢，在航空測繪、攝影和森林植保等領(lǐng)域得到越來越廣泛的應(yīng)用[1]。氣動布局設(shè)計是無人機(jī)設(shè)計過程中的關(guān)鍵，通過對全機(jī)氣動特性進(jìn)行分析，可以用于進(jìn)一步分析飛機(jī)的飛行性能，為控制系統(tǒng)設(shè)計提供重要依據(jù)[2-4]。計算流體力學(xué)（Computational Fluid Dynamics，CFD）是指通過計算機(jī)利用數(shù)值方法求解流體動力學(xué)方程，獲得流動規(guī)律并解決流動問題。利用CFD技術(shù)可以模擬流體的真實狀況，快速獲得氣動特性數(shù)據(jù)[5-6]。

航空工程領(lǐng)域是CFD的最大推動者之一，隨著計算精度要求的提高，與此同時需要更大規(guī)模的計算網(wǎng)格，傳統(tǒng)串行的計算方法已經(jīng)難以滿足需求，必須采用并行計算的方法解決大規(guī)模流場數(shù)值計算耗時過多的問題。多核并行計算是指同時對多個任務(wù)、多條指令或多個數(shù)據(jù)項進(jìn)行處理[7]。文獻(xiàn)[8]對比了多核并行集群系統(tǒng)與單核集群的計算時間，驗證了多核集群的高效性和可擴(kuò)展性；文獻(xiàn)[9]基于以太網(wǎng)連接的工作站機(jī)群，采用FLUENT并行計算方法對某型號飛機(jī)的非結(jié)構(gòu)網(wǎng)格進(jìn)行了數(shù)值求解，并討論了計算節(jié)點數(shù)對計算時間、加速比和并行效率的影響；文獻(xiàn)[10]采用CFD并行計算的方法對乘波飛行器前體模型進(jìn)行了計算。然而基于多節(jié)點的計算集群受節(jié)點間通信的限制，加速比難以達(dá)到理想值[11]。

隨著硬件技術(shù)的快速發(fā)展，單機(jī)多核CPU應(yīng)運(yùn)而生，它將幾個甚至幾十個CPU集成到單個芯片上，每個CPU核作為一個單獨的處理器，從而進(jìn)一步提高并行計算的速度。

為了提高飛機(jī)總體設(shè)計中氣動布局設(shè)計的效率，需要獲得一種高效的數(shù)值模擬方法，在提高數(shù)值計算精確度的同時縮短計算時間。本文首先采用串行的計算方法分別對基于Euler方程和N-S（Navier-Stokes）方程的求解器進(jìn)行驗證，將求解結(jié)果與ONERA M6機(jī)翼[12]實驗結(jié)果進(jìn)行對比，以驗證主控方程對計算精度的影響。接著采用FLUENT并行的計算方法，討論了線程數(shù)對M6機(jī)翼流場數(shù)值計算時間和加速比的影響，尋找一種基于FLUENT的高效并行求解方法。最后得出一種準(zhǔn)確、高效的計算方案，并針對某V型尾翼、上單翼布局無人機(jī)的整機(jī)氣動外形進(jìn)行了計算。

2CFD并行計算技術(shù)

CFD技術(shù)與理論分析、實驗研究共同組成空氣動力研究的三種主要手段，相較于其他兩種方法具有研制周期短、損耗小的特點，同時不需考慮洞壁干擾和支架干擾，可以在計算機(jī)上實現(xiàn)各種條件下的流場計算。隨著計算機(jī)的飛速發(fā)展，使用計算機(jī)對工程現(xiàn)象進(jìn)行數(shù)值模擬的分析手段成為可能，并且計算精度和效率都大大依賴于計算機(jī)的硬件設(shè)備。

CFD任務(wù)能夠分解為可并行計算的多個子任務(wù)，在不同處理器上同時處理，這就為并行計算提供了可行性。采用并行計算一方面能夠解決串行計算時間長的問題[13-15]，另一方面可以在相同計算時間內(nèi)擴(kuò)大求解規(guī)模，提高計算精度[16-18]。

2.1　CFD數(shù)值模擬的控制方程及離散方法

Navier-Stokes方程組的求解是CFD的核心，該方程組是典型連續(xù)介質(zhì)流體力學(xué)的控制方程。Navier-Stokes方程組包括質(zhì)量守恒、動量守恒和能量守恒方程，直角坐標(biāo)系下，積分形式的三維守恒型N-S方程形式如下：

式中Q=[ρ,ρu,ρv,ρw,ρe]T表示守恒向量，n為邊界外法向量，?V為某一固定區(qū)域V的邊界。

分別表示無粘項和粘性項，ρ、e分別為密度和單位質(zhì)量氣動總能量，(u,v,w)為直角坐標(biāo)系下的速度分量。其中粘性應(yīng)力項：

無粘項空間離散采用迎風(fēng)通量差分分裂格式，對流項和壓力項可以通過一個網(wǎng)格單元的通量平衡來表示：

其中UL、UR是交界面的原始參數(shù)。

粘性項采用二階差分格式進(jìn)行離散：

時間離散采用近似因子分解法（Approximate Factorization Method），該方法是一種基于多維非線性控制方程的隱式方法，從而提高總體計算效率。

2.2　多核并行技術(shù)

不同于單核CPU多線程的交錯執(zhí)行，多核CPU的多線程在物理上是并行執(zhí)行的。要充分發(fā)揮多核CPU的硬件性能，必須采用多線程執(zhí)行，使得每個CPU核在同一時刻都有線程在執(zhí)行。

ANSYS FLUENT CFD（以下簡稱FLUENT）中的并行處理包括FLUENT、主機(jī)進(jìn)程和一套節(jié)點進(jìn)程之間的交互作用，通過CORTEX模塊與主機(jī)進(jìn)程及各個節(jié)點進(jìn)程進(jìn)行信息傳遞。串行求解器各模塊間邏輯關(guān)系如圖1所示。

圖1　串行求解器各模塊間的邏輯關(guān)系

FLUENT的并行計算可在大型并行計算機(jī)上運(yùn)行，也可在一個多CPU的工作站上運(yùn)行。CORTEX發(fā)出的指令通過主機(jī)進(jìn)行解釋，由主機(jī)進(jìn)程通過指定計算節(jié)點的數(shù)據(jù)通訊模塊傳遞給其他節(jié)點，被指定的計算節(jié)點稱為“計算節(jié)點1”，其他節(jié)點在接收主機(jī)進(jìn)程傳遞的命令后，在各自的數(shù)據(jù)集合上同時執(zhí)行相同的程序。各計算節(jié)點彼此相連，節(jié)點間數(shù)據(jù)傳輸、同步和執(zhí)行全局任務(wù)則需要通過數(shù)據(jù)通訊模塊實現(xiàn)，從而實現(xiàn)邏輯意義聯(lián)通。FLUENT并行求解器各模塊間的邏輯關(guān)系如圖2所示。

圖2　并行求解器各模塊間的邏輯關(guān)系

MPI（Message Passing Interface）是一種信息傳遞編程標(biāo)準(zhǔn)，為基于信息傳遞的并行程序提供一個高效統(tǒng)一的編程環(huán)境，是目前應(yīng)用最為廣泛的并行編程方式。FLUENT自帶的MPI并行機(jī)制不僅支持對稱多處理共享存儲并行計算，也支持網(wǎng)絡(luò)分布式并行計算，由于多核CPU屬于共享存儲并行機(jī)，因此，采用FLUENT內(nèi)置MPI并行機(jī)制能夠有效提高單機(jī)多核并行計算的效率。

3　計算應(yīng)用與分析

3.1　計算條件與并行環(huán)境

為了獲得高效準(zhǔn)確的計算方法，首先對算例機(jī)翼進(jìn)行數(shù)值求解，驗證求解器的準(zhǔn)確性。數(shù)值計算采用基于密度的耦合隱式算法，控制方程在計算域的離散方法采用有限體積方法，機(jī)體表面滿足無滑移邊界條件，遠(yuǎn)場為自由可壓縮流場，湍流模型為剪切滑移SST（Shear-Stress Transport）模型。

由于計算工況較多，網(wǎng)格數(shù)量較大，普通計算機(jī)串行計算將耗費(fèi)巨大的時間，因此需要并行計算進(jìn)行實驗。實驗采用的多核工作站配置兩個2.6 GHz的CPU，內(nèi)存為128 GB，48核數(shù)運(yùn)行處理器。

3.2　串行與并行求解器

本文選取用于驗證計算方法的ONERA M6實驗機(jī)翼作為驗證算例，計算條件如表1所示，采用串行的計算方法。

表1　ONERA M6機(jī)翼來流屬性及計算狀態(tài)

氣動網(wǎng)格為ANSYS ICEM CFD軟件生成的結(jié)構(gòu)網(wǎng)格，空間分區(qū)及機(jī)翼表面網(wǎng)格如圖3所示。分別采用N-S方程和Euler方程對流場進(jìn)行了數(shù)值求解，并與實驗數(shù)據(jù)進(jìn)行了對比。沿機(jī)翼展向截面的壓力系數(shù)如圖4所示，可以看出N-S解與實驗值吻合很好，這是由于該狀態(tài)下激波的作用導(dǎo)致流動分離，N-S解由于粘性的計入能夠更準(zhǔn)確地描述該狀態(tài)下截面的附體流動?？傊@種三維結(jié)構(gòu)網(wǎng)格劃分方法以及流場解算器效果是很好的，可以作為無人機(jī)流場數(shù)值求解的方法。

圖4　ONERA M6亞音速狀況下的壓強(qiáng)分布（z l=0.2）

加速比是衡量并行計算效果的核心指標(biāo)，加速比通常是指在相同規(guī)模下串行執(zhí)行時間與并行執(zhí)行時間的比值。若Tp表示并行系統(tǒng)下的計算時間，Ts表示串行系統(tǒng)計算時間，那么并行化加速比S可以用下述公式來表示：

分別采用串行和多核并行系統(tǒng)對模型進(jìn)行計算，設(shè)置一組FLUENT并行計算的線程數(shù)目，分別為1、6、12、24、36、48，表2、表3列出了1 000個時間推進(jìn)步的計算時間和計算加速比。由表2可以看出，并行計算可以大大縮短運(yùn)算時間，隨著線程數(shù)目的增加，計算效率得到進(jìn)一步提高，當(dāng)多線程所需的開銷增大到一定程度則會降低運(yùn)算速度。隨著網(wǎng)格數(shù)量的增加，多線程并行時加速比普遍略有所降低。由此可見，合理選擇并行計算的線程數(shù)，可以充分發(fā)揮并行計算的優(yōu)勢，并用于大規(guī)模的求解問題中。由圖5可以看出，不同線程的多核并行計算的壓強(qiáng)分布吻合一致，且與實驗值吻合較好。

表2　多核并行計算與串行計算執(zhí)行時間對比

表3　多線程數(shù)并行計算加速比

3.3　數(shù)值計算與分析

網(wǎng)格生成與數(shù)值模擬的結(jié)果密切相關(guān)，網(wǎng)格質(zhì)量的好壞直接關(guān)系到數(shù)值模擬結(jié)果的正確性。本次實驗前處理模型仍采用ANSYS ICEM CFD軟件生成計算域內(nèi)結(jié)構(gòu)網(wǎng)格。由于計算模型的對稱性，首先對半?yún)^(qū)幾何模型進(jìn)行分區(qū)結(jié)構(gòu)網(wǎng)格劃分，通過網(wǎng)格鏡像的功能生成全機(jī)網(wǎng)格。

數(shù)值計算內(nèi)容為某測繪無人機(jī)無動力狀態(tài)下的氣動力參數(shù)，采用上單翼和V型尾翼布局形式。計算條件為：Ma=0.065，Re=4.87×105，選擇三維可壓縮N-S粘性方程作為流場主控方程，采用FLUENT軟件進(jìn)行多核并行數(shù)值計算，通過設(shè)定線程數(shù)目，提高計算效率。求解時利用軟件自帶的腳本記錄功能，將不同飛行工況計算的步驟記錄下來，通過讀入腳本實現(xiàn)不同工況自動計算和保存，并在計算結(jié)束后自動切換下一工況。這種方法避免了計算資源的浪費(fèi)，同時減少因切換工況的時間花費(fèi)。圖6為整機(jī)升力系數(shù)和阻力系數(shù)隨迎角的變化曲線；圖7為機(jī)體表面的等值線分布和表面壓力系數(shù)分布云圖。

4　結(jié)論

（1）本文提出了一種用于求解無人機(jī)整機(jī)的精確CFD數(shù)值模擬方法。分別采用Euler法和N-S法對算例機(jī)翼進(jìn)行解算，驗證了以N-S方程作為流場主控方程的準(zhǔn)確性，并采用該解算器對某型號民用無人機(jī)流場的氣動特性進(jìn)行了數(shù)值模擬。

圖5　ONERA M6機(jī)翼不同線程數(shù)各截面的壓強(qiáng)分布

圖6　整機(jī)氣動特性曲線

圖7　機(jī)體表面等值線分布和壓強(qiáng)分布云圖

（2）提出了一種用于求解無人機(jī)整機(jī)的單機(jī)并行計算方法。本文基于多核系統(tǒng)對算例機(jī)翼的求解執(zhí)行時間進(jìn)行了探討，驗證了并行求解的準(zhǔn)確性和高效性，通過設(shè)置不同線程數(shù)目執(zhí)行求解程序，使CPU不同核心同時進(jìn)行計算，達(dá)到加速的目的。隨著線程數(shù)的增加，加速比的提升幅度先大后小，達(dá)到最大值后會有所減小。將這種思路用于大規(guī)模的無人機(jī)整機(jī)數(shù)值模擬問題中，使得計算效率得到大幅度提高。

（3）由于多核處理器多個CPU核的共享存儲模式，相對于集群計算平臺，節(jié)點間的通信時間幾乎可以忽略，可同時兼顧計算機(jī)硬件結(jié)構(gòu)和發(fā)揮CPU最大計算能力，具有很大工程應(yīng)用潛力。進(jìn)一步的研究內(nèi)容包括：流場網(wǎng)格分區(qū)數(shù)與計算節(jié)點數(shù)的關(guān)系，即如何對三維流場進(jìn)行分區(qū)，并選擇合適的計算節(jié)點，避免硬件資源的浪費(fèi)。

參考文獻(xiàn)：

[1]崔秀敏，王維軍，方振平.小型無人機(jī)發(fā)展現(xiàn)狀及其相關(guān)問題分析[J].飛行力學(xué)，2005，23（1）：14-18.

[2]祝小平，向錦武，張才文，等.無人機(jī)設(shè)計手冊[M].北京：國防工業(yè)出版社，2007.

[3]段鎮(zhèn)，高九州，賈宏光，等.無人機(jī)滑跑線性化建模與增益調(diào)節(jié)糾偏控制[J].光學(xué)精密工程，2014，22（6）：1507-1516.

[4]李迪，陳向堅，續(xù)志軍.增益自適應(yīng)滑模控制器在微型飛行器姿態(tài)控制中的應(yīng)用[J].光學(xué)精密工程，2013，21（5）：1183-1191.

[5]關(guān)鍵.低速無人機(jī)動態(tài)氣動特性數(shù)值模擬及布局研究[D].長沙：國防科技大學(xué)，2013.

[6]Baker A M，Ying S X.A fixed time performance evaluation of parallel CFD applications[J].Super Computing，2002，95：18-23.

[7]周偉明.多核計算與程序設(shè)計[M].武漢：華中科技大學(xué)出版社，2009.

[8]陳堅禎，陽平，李斌，等，多核并行計算下的流量傳感器流場模擬研究[J].衡陽師范學(xué)院學(xué)報，2011，32（6）：82-84.

[9]邵波，毛國勇，張武.基于Fluent的全機(jī)數(shù)值模擬及并行計算[J].計算機(jī)工程與設(shè)計，2006，27（17）：3178-3180.

[10]潘沙，李樺，夏智勛.高性能并行計算在航空航天CFD數(shù)值模擬中的應(yīng)用[J].計算機(jī)工程與科學(xué)，2012，34（8）：191-198.

[11]熊瑋，夏文龍，余曉鴻，等.多核并行計算技術(shù)在電力系統(tǒng)短路計算中的應(yīng)用[J].電力系統(tǒng)自動化，2011，35（8）：49-52.

[12]Schmitt V，Charpin F.Pressure distributions on the ONERAM6 wing at transonic mach numbers，AGARD-AR-138[R].1979.

[13]馮定華，潘沙，丁國昊，等.CFD并行計算方法在高超聲速數(shù)值模擬中的應(yīng)用[C]//第三屆高超聲速科技學(xué)術(shù)會議，2010.

[14]吳頎.GPU并行計算及其在飛行器設(shè)計中的應(yīng)用[D].北京：北京理工大學(xué)，2015.

[15]唐逸豪，高振勛，蔣崇文，等.基于LPT近似算法的CFD并行計算網(wǎng)格分配算法[J].工程力學(xué)，2015，32（5）：243-249.

[16]Anderson J D.Computational fluid dynamics[M].[S.l.]：Mc Graw Hill Education，2007.

[17]劉巍，張理論，鄧小剛.計算空氣動力學(xué)并行編程基礎(chǔ)[M].北京：國防工業(yè)出版社，2013.

[18]溫建華，朱自強(qiáng)，吳宗成.基于N-S方程串并行計算的機(jī)翼優(yōu)化設(shè)計[J].北京航空航天大學(xué)學(xué)報，2008，34（2）：127-130.