基于多核并行計算的無人機CFD數值模擬

姜悅寧，賈宏光，厲　明

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

1.中國科學院 長春光學精密機械與物理研究所，長春 130033

2.中國科學院大學，北京 100039

3.長光衛星技術有限公司，長春 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　引言

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

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

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

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

2CFD并行計算技術

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

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

2.1　CFD數值模擬的控制方程及離散方法

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

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

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

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

其中UL、UR是交界面的原始參數。

粘性項采用二階差分格式進行離散：

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

2.2　多核并行技術

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

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

圖1　串行求解器各模塊間的邏輯關系

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

圖2　并行求解器各模塊間的邏輯關系

MPI（Message Passing Interface）是一種信息傳遞編程標準，為基于信息傳遞的并行程序提供一個高效統一的編程環境，是目前應用最為廣泛的并行編程方式。FLUENT自帶的MPI并行機制不僅支持對稱多處理共享存儲并行計算，也支持網絡分布式并行計算，由于多核CPU屬于共享存儲并行機，因此，采用FLUENT內置MPI并行機制能夠有效提高單機多核并行計算的效率。

3　計算應用與分析

3.1　計算條件與并行環境

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

由于計算工況較多，網格數量較大，普通計算機串行計算將耗費巨大的時間，因此需要并行計算進行實驗。實驗采用的多核工作站配置兩個2.6 GHz的CPU，內存為128 GB，48核數運行處理器。

3.2　串行與并行求解器

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

表1　ONERA M6機翼來流屬性及計算狀態

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

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

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

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

表2　多核并行計算與串行計算執行時間對比

表3　多線程數并行計算加速比

3.3　數值計算與分析

網格生成與數值模擬的結果密切相關，網格質量的好壞直接關系到數值模擬結果的正確性。本次實驗前處理模型仍采用ANSYS ICEM CFD軟件生成計算域內結構網格。由于計算模型的對稱性，首先對半區幾何模型進行分區結構網格劃分，通過網格鏡像的功能生成全機網格。

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

4　結論

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

圖5　ONERA M6機翼不同線程數各截面的壓強分布

圖6　整機氣動特性曲線

圖7　機體表面等值線分布和壓強分布云圖

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

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

參考文獻：

[1]崔秀敏，王維軍，方振平.小型無人機發展現狀及其相關問題分析[J].飛行力學，2005，23（1）：14-18.

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

[3]段鎮，高九州，賈宏光，等.無人機滑跑線性化建模與增益調節糾偏控制[J].光學精密工程，2014，22（6）：1507-1516.

[4]李迪，陳向堅，續志軍.增益自適應滑?？刂破髟谖⑿惋w行器姿態控制中的應用[J].光學精密工程，2013，21（5）：1183-1191.

[5]關鍵.低速無人機動態氣動特性數值模擬及布局研究[D].長沙：國防科技大學，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]周偉明.多核計算與程序設計[M].武漢：華中科技大學出版社，2009.

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

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

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

[11]熊瑋，夏文龍，余曉鴻，等.多核并行計算技術在電力系統短路計算中的應用[J].電力系統自動化，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并行計算方法在高超聲速數值模擬中的應用[C]//第三屆高超聲速科技學術會議，2010.

[14]吳頎.GPU并行計算及其在飛行器設計中的應用[D].北京：北京理工大學，2015.

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

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

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

[18]溫建華，朱自強，吳宗成.基于N-S方程串并行計算的機翼優化設計[J].北京航空航天大學學報，2008，34（2）：127-130.