, , ,*,

(1. 西北工業(yè)大學(xué) 航空學(xué)院, 陜西 西安 710072; 2. 中國空氣動力研究與發(fā)展中心 計算空氣動力研究所, 四川 綿陽 621000)

0 引 言

大規(guī)模分布式并行計算快速發(fā)展，使計算流體力學(xué)(Computational Fluid Dynamics, CFD)具備了解決復(fù)雜外形工程問題的能力。大規(guī)模CFD并行計算成為飛行器氣動性能計算[1]、多體分離安全性評估[2]和飛行器外形優(yōu)化[3-5]的有效手段。

CFD并行計算屬于數(shù)據(jù)通信密集型，采用的數(shù)值格式和軟件并行設(shè)計對并行效率影響很大[6]。Knight首先系統(tǒng)性地闡述了并行計算機系統(tǒng)架構(gòu)和CFD并行程序設(shè)計模式[7]。 國內(nèi)陸林生和董超群等[8]研究了并行程序概念設(shè)計方法，陳剛和王磊等[9]提出了面向大規(guī)模并行的軟件實現(xiàn)方法。近年來，CFD并行計算廣泛應(yīng)用于航空航天等領(lǐng)域[10-12]，但并行規(guī)模都在萬核以下。

對于分布式并行系統(tǒng)，消息傳遞接口(Message Passing Interface, MPI)是最成熟和有效的并行程序?qū)崿F(xiàn)規(guī)范，廣泛應(yīng)用于計算力學(xué)領(lǐng)域。OpenMPI和MPICH是兩種廣泛使用的MPI編程程序包，都成功用于CFD并行程序編碼[13-14]。

在使用并行解算器模擬流動前，通常需要采用合適的方法對生成的網(wǎng)格進行多核并行分區(qū)。分區(qū)算法極大地影響著并行計算時的負載平衡和并行效率。對于非結(jié)構(gòu)網(wǎng)格，采用多級不規(guī)則圖算法的Metis程序庫[15]，是目前應(yīng)用最廣、最快速和高質(zhì)量的并行分區(qū)工具。由于沒有拓撲結(jié)構(gòu)的限制，分區(qū)后的非結(jié)構(gòu)網(wǎng)格通常比結(jié)構(gòu)網(wǎng)格有更高的負載平衡性能。

非結(jié)構(gòu)網(wǎng)格的有限體積方法通常采用積分形式的雷諾平均Navier-Stokes方程和湍流模型方程[16]。二階MUSCL (Monotonic Upwind Scheme for Conservation Laws) 類數(shù)值格式具有類似緊致性質(zhì)，即某個網(wǎng)格單元通量的計算只需要臨近共面單元的流場信息，因而廣泛應(yīng)用于CFD計算軟件且發(fā)展較為成熟。并行計算時，僅需要通過MPI交換并行交界面兩側(cè)單元的流場數(shù)據(jù)，有利于并行計算的實現(xiàn)和減少并行傳遞的數(shù)據(jù)量。

本文通過基于非結(jié)構(gòu)混合網(wǎng)格和有限體積方法，采用緊致型二階數(shù)值格式，發(fā)展了適用于復(fù)雜工程問題的萬核級并行實現(xiàn)方法，測試和評估了并行效率。

1 數(shù)值方法

雷諾平均Navier-Stokes方程(RANS)守恒形式的積分方程為：

(1)

對流通量可以表示為：

(2)

式中，p是靜壓，ρ是密度，V是逆變速度:

V=nxu+nyv+nzw

(3)

H為焓值，定義為：

(4)

式中，e為流動內(nèi)能。

黏性通量相對比較復(fù)雜，本文略去，具體內(nèi)容可以參考文獻[17]。

在網(wǎng)格單元內(nèi)積分，得到離散形式的控制方程：

(5)

式中，V是網(wǎng)格單元體積，N是與該網(wǎng)格單元共面網(wǎng)格單元總數(shù)。

控制方程的離散采用格心格式的有限體積法，時間項采用LU-SGS迭代[18]，湍流模型為一方程的SA模型[19]。對流項采用Roe格式[20]，黏性項采用中心格式[21]。

格心格式的有限體積方法，流場未知變量存儲在網(wǎng)格單元的形心上，如圖1所示的三角形形心點J1、J2和J3。然而，對流通量和黏性通量計算卻需要網(wǎng)格單元對應(yīng)的網(wǎng)格面心上的流場值，如K1、K2和K3。只有采用一階格式時，面上的流場值等于形心值。當(dāng)采用緊致型的二階格式時，網(wǎng)格面心上左右兩側(cè)的流場值由同側(cè)的網(wǎng)格單元體心值構(gòu)造。計算黏性通量時，面心流場值采用中心格式，以K1面為例，表示為:

(6)

圖1 通量計算模板示意圖Fig.1 Stencil for flux calculation

K2和K3面心流場值可以類似得到。黏性通量可以表示為:

Hv,K1=f(QK1)

(7)

黏性通量計算的具體表達式參考文獻[21]。

采用Roe格式時，對流通量可以表示為

Hc,K1=f(QK1,L+QK1,R)

(8)

式中，QK1,L和QK1,R分別為面K1左右兩側(cè)的流場值，通過各自同側(cè)的體心值構(gòu)造

(9)

采用Roe格式計算對流通量的具體表達式參考文獻[20]。

需要指出的是，如果采用更高階的數(shù)值格式，面心流場值的構(gòu)造將會用到更多臨近網(wǎng)格單元上的體心值。對于本文所述的二階格式，只需要緊鄰的網(wǎng)格單元體心值，因而并行計算時，只有并行交界面兩側(cè)的網(wǎng)格單元需要進行流場數(shù)據(jù)的交換。

式(5)中的第二項通常稱為殘差項，記為:

(10)

則式(5)表示為:

(11)

為了增加計算的數(shù)值穩(wěn)定性和提高收斂速度，時間項通常采用隱式格式[17]，即

(12)

式中：

(13)

通過對殘差項的線化，可以得到

(14)

上式求解采用LU-SGS迭代方法，具體內(nèi)容可參考文獻[18]。定常流場由邊界條件決定，迭代過程中，邊界信息由邊界單元逐漸向流場內(nèi)部傳遞，直到收斂。通常情況下，可能需要迭代上千步甚至更多流場才能收斂，當(dāng)網(wǎng)格量較大時，計算量巨大，單核計算需要很長的時間。對于工程復(fù)雜外形，并行計算是減少計算時間的有效方法。

2 并行實現(xiàn)算法

采用LU-SGS迭代方法，當(dāng)前步流場計算依賴于前一步的流場，當(dāng)流場為非定常時依賴更多。但對于某一迭代步，采用緊致二階格式，網(wǎng)格單元的流場值只與相鄰一層網(wǎng)格有關(guān)。因此對于分布式并行計算系統(tǒng)，采用對空間進行分區(qū)更適合CFD并行計算，并采用MPI進行多個分區(qū)間流場數(shù)據(jù)的交換。并行實現(xiàn)算法主要包括網(wǎng)格分區(qū)算法和流場數(shù)據(jù)交換算法。

2.1 網(wǎng)格分區(qū)技術(shù)

非結(jié)構(gòu)網(wǎng)格分區(qū)方法采用基于圖論算法的Metis程序庫，該程序庫使用多級二分法或多級多支分區(qū)算法，如圖2所示。

圖2 Metis多級分區(qū)算法示意圖[15]Fig.2 Multi-level partition used by Meits[15]

分區(qū)算法以網(wǎng)格數(shù)據(jù)構(gòu)成的圖作為輸入，由網(wǎng)格面和網(wǎng)格單元組成。包含n個節(jié)點V和m個邊E的圖G定義為:

(15)

需要指出的是，對于格心格式的有限體積法，圖定義中的節(jié)點和邊分別指網(wǎng)格中的網(wǎng)格單元和網(wǎng)格面。輸入的圖信息通常采用壓縮存儲格式(Compressed Storage Format， CSR)，該格式廣泛應(yīng)用于稀疏圖的存儲。網(wǎng)格連接關(guān)系表示成圖的信息可以用兩個整數(shù)數(shù)組來指定：xadj和adjncy，數(shù)組長度分別為n+1和2m。xadj[i+1]-xadj[i]的值表示與第i個網(wǎng)格單元共面的網(wǎng)格單元個數(shù)，數(shù)組adjncy中第xadj[i] 到第xadj[i+1] 位置存儲了相應(yīng)的共面網(wǎng)格單元的編號。以圖1網(wǎng)格數(shù)據(jù)為例，數(shù)組xadj和adjncy存儲數(shù)據(jù)如圖3所示。

圖3 非結(jié)構(gòu)網(wǎng)格連接關(guān)系CSR表示Fig.3 Graph data in CSR format of unstructured grid

圖4給出了M6機翼網(wǎng)格并行分區(qū)表面網(wǎng)格的示意圖，共3個并行分區(qū)，其中相同顏色的網(wǎng)格單元屬于相同的計算進程。

2.2 數(shù)據(jù)交換并行實現(xiàn)

本文并行算法基于分布式大規(guī)模并行計算機系統(tǒng)，采用所有CPU都參與迭代計算的對等并行模式，MPI數(shù)據(jù)傳遞采用非阻塞通信，盡量加大計算和通信的重疊，只在必要的地方進行數(shù)據(jù)同步以保證計算的正確性。

圖4 M6機翼網(wǎng)格分區(qū)圖Fig.4 Partitioned grid of M6 wing with Metis

采用對等并行模式，首先需要對網(wǎng)格進行前處理，按指定的并行規(guī)模采用Metis劃分網(wǎng)格。為了減少每個計算進程內(nèi)存消耗以提高程序的運行效率，每個分區(qū)網(wǎng)格只保存必需的局部網(wǎng)格數(shù)據(jù)和并行分區(qū)邊界上網(wǎng)格單元及頂點的對應(yīng)關(guān)系。網(wǎng)格分區(qū)完成后，每個計算進程讀入相應(yīng)分區(qū)的網(wǎng)格，獨立地逐一計算各自分區(qū)中每個網(wǎng)格單元的通量。

格心格式的有限體積法，并行邊界間數(shù)據(jù)傳遞需要的對應(yīng)關(guān)系包括網(wǎng)格單元與網(wǎng)格單元和網(wǎng)格頂點與網(wǎng)格頂點間的對應(yīng)關(guān)系，其中網(wǎng)格單元間對應(yīng)關(guān)系主要用于流場變量的迭代，網(wǎng)格頂點間的對應(yīng)關(guān)系用于二階格式中梯度的計算。并行網(wǎng)格對應(yīng)關(guān)系建立耗時長，然而對于定常問題只需要計算一次，程序設(shè)計時宜采用存儲該對應(yīng)關(guān)系，以空間換時間，提高并行效率。

對于串行計算，封閉的邊界條件決定了當(dāng)前流場的特征，邊界信息由邊界網(wǎng)格單元向內(nèi)部網(wǎng)格單元傳播。并行計算時，對于單個計算進程，只有部分邊界網(wǎng)格存在相應(yīng)的分區(qū)網(wǎng)格中，邊界條件不再封閉。將并行交界面上的網(wǎng)格面視為特殊類型的邊界面時，則邊界面滿足封閉性質(zhì)。同時，與真實邊界面類似地采用虛擬邊界單元，真實邊界面對應(yīng)的虛擬單元用于施加邊界條件，并行邊界面對應(yīng)的虛擬單元則用來存儲并行交換的數(shù)據(jù)。在計算網(wǎng)格單元的通量時，并行邊界面與真實邊界面采用相同的計算方法，程序設(shè)計亦能與串行程序保持一致。并行邊界與虛擬單元設(shè)計如圖5所示，則對應(yīng)任意并行計算進程，計算及存儲模式和串行時保持一致。

MPI并行接口中并行數(shù)據(jù)必須以連續(xù)內(nèi)存的方式提供，程序設(shè)計必須考慮減少并行計算時額外數(shù)據(jù)拷貝的開銷以提高并行效率。格心格式有限體積法任何分區(qū)內(nèi)的某個單元只對應(yīng)唯一其它分區(qū)內(nèi)的唯一單元，可以采用先物理網(wǎng)格單元后虛擬網(wǎng)格單元的網(wǎng)格排序策略，使對應(yīng)的每一個并行分區(qū)需要接收的數(shù)據(jù)在內(nèi)存中集中連續(xù)存放，如圖6所示，消除了網(wǎng)格單元并行數(shù)據(jù)的拷貝過程。

圖5 并行邊界及虛擬單元Fig.5 Ghost elements for data exchange between different processors

圖6 并行物理網(wǎng)格單元及虛擬網(wǎng)格單元排序存儲示意圖Fig.6 Memory pattern for physical elements and ghost elements of grid data

3 并行測試結(jié)果

采用兩個算例分別對本文非結(jié)構(gòu)流動解算器的并行正確性和并行效率進行測試。算例一使用網(wǎng)格量相對較少的旋成體外形，迭代計算直到流場收斂，比較分析串行計算和并行計算結(jié)果的一致性；算例二使用運輸機構(gòu)型，其網(wǎng)格量大，用于測試大規(guī)模并行計算的效率，其最大的并行規(guī)模達18816核。

并行加速比S和效率E定義為：

(16)

式中：n為測試并行核數(shù)，T為每步計算平均時間；n0為并行測試基數(shù)核數(shù)，T0為基數(shù)核數(shù)對應(yīng)的每步計算平均時間。

并行計算使用中國空氣動力研究與發(fā)展中心計算空氣動力研究所分布式并行計算機系統(tǒng)。該系統(tǒng)包含300個計算節(jié)點，每個節(jié)點64核，CPU主頻約為2.8 GHz。

3.1 并行正確性及收斂性測試

正確性測試模型為旋成體外形繞流模擬，旨在分析并行和串行計算結(jié)果的一致性。來流條件為：速度V= 1200 m/s，高度H=24 km，側(cè)滑角β= 0°，迎角α=2.0°。由于流動對稱，宜采用半模計算，網(wǎng)格如圖7所示，網(wǎng)格單元總數(shù)約為200萬，單元類型包括四面體、棱柱和金字塔。計算核數(shù)分別為：1、 32、64、128、256、512、1024、2048、4096和8192。

圖7 旋成體外形計算網(wǎng)格Fig.7 Computing grid of revolving body

圖8～圖10分別給出了升力、阻力和俯仰力矩系數(shù)的收斂歷程。由于并行計算和串行計算迭代過程中，網(wǎng)格單元通量計算的順序并不完全相同，因此收斂過程有小量的差別，但收斂趨勢一致且收斂后的值在誤差范圍內(nèi)相同。不同核數(shù)對應(yīng)收斂后的數(shù)值列于表1，可以看出，并行計算收斂數(shù)值和串行計算偏差小于10-5。

圖8 升力系數(shù)收斂曲線Fig.8 Convergence history of lift

圖10 俯仰力矩系數(shù)收斂曲線Fig.10 Convergence history of pitch moment

CoresCoefficient of liftDeviationCoefficient of dragDeviationCoefficient of pitch momentDeviation10.0710874--0.137338---0.0433570--320.0710849-2.50×10-60.137318-2.00×10-5-0.04335691.00×10-7640.07108972.30×10-60.137325-1.30×10-5-0.0433607-3.70×10-61280.07108972.30×10-60.137325-1.30×10-5-0.0433607-3.70×10-62560.0710849-2.50×10-60.137318-2.00×10-5-0.04335691.00×10-75120.07109083.40×10-60.1373435.00×10-6-0.0433632-6.20×10-610240.07108972.30×10-60.137325-1.30×10-5-0.0433607-3.70×10-620480.0710868-6.00×10-70.137324-1.40×10-5-0.0433588-1.80×10-640960.07108982.40×10-60.137325-1.30×10-5-0.0433608-3.80×10-681920.07108982.40×10-60.137325-1.30×10-5-0.0433608-3.80×10-6Max.0.07109083.40×10-60.1373435.00×10-6-0.04335691.00×10-7Min.0.0710849-2.50×10-60.137318-2.00×10-5-0.0433632-6.20×10-6

圖11和圖12分別給出了加速比和并行效率隨計算核數(shù)增加的變化，可以看出，512核是加速比快速偏離理論值的臨界核數(shù)，大于該值后效率快速降低，主要原因為網(wǎng)格單元很少，并行通信占用較多的時間。

圖13和圖14分別給出了并行核數(shù)為64、256、1024和4096四種情況下密度和殘差的收斂曲線，隨著并行核數(shù)的增加，除256核外其它并行規(guī)模下殘差收斂速度基本一致。

圖11 加速比隨計算核數(shù)變化Fig.11 Speed up versus processors

圖12 并行效率隨計算核數(shù)變化Fig.12 Efficiency versus processors

圖13 密度殘差收斂曲線Fig.13 Convergence history of density residual

圖14 能量殘差收斂曲線Fig.14 Convergence history of energy residual

3.2 大規(guī)模并行效率測試

用于并行效率測試的模型為運輸機繞流模擬，網(wǎng)格單元總數(shù)約為1.2億。64個并行分區(qū)后的表面網(wǎng)格如圖15所示，由飛機頭部和翼梢局部放大圖可以看出，計算網(wǎng)格很密，其有利于準(zhǔn)確模擬飛行阻力。來流條件為：飛行速度V=300 m/s，飛行高度H=11 km，側(cè)滑角β= 0°，迎角α= 2.0°。

圖15 運輸機構(gòu)型并行分區(qū)后的表面網(wǎng)格Fig.15 Partitioned surface grid for transport aircraft

為了增加可靠性，效率測試每步計算時間是取迭代200步的平均時間。由于網(wǎng)格量大，單核計算時間太長，因此并行加速比S和效率E測試都是以64核為基準(zhǔn)。

并行加速比隨并行核數(shù)增加變化曲線如圖16所示，可以看出，該解算器并行性能很好，加速比為線性甚至超線性。出現(xiàn)超線性的原因主要有：a) 非結(jié)構(gòu)網(wǎng)格沒有拓撲限制，采用Metis并行分區(qū)后容易達到很高的負載平衡性能；b) 在分布式并行系統(tǒng)中，每個計算進程擁有其單獨的緩存。隨著并行規(guī)模增加，緩存總量增加，因此每個核擁有的緩存增加，提高了緩存命中率；c) 當(dāng)并行規(guī)模增加后，每個核網(wǎng)格量減小，需要的內(nèi)存等資源減少，有利于提高并行系統(tǒng)的性能。

圖16 加速比隨計算核數(shù)變化Fig.16 Speed up versus processors

并行效率隨并行規(guī)模增加變化曲線如圖17所示，可以看出，由于加速比性能優(yōu)越，在并行核數(shù)小于14336時，并行效率都在100%以上，直到18816核，并行效率都保持在80%以上。

圖17 并行效率隨計算核數(shù)變化Fig.17 Efficiency versus processors

圖18給出了并行效率與每核網(wǎng)格單元數(shù)的關(guān)系，當(dāng)網(wǎng)格單元數(shù)太少，并行數(shù)據(jù)交換耗費的時間相比迭代計算時間不再是小量時，并行效率隨網(wǎng)格單元數(shù)減少快速下降。不同的數(shù)值格式和并行實現(xiàn)下降的拐點不同，本文在每核網(wǎng)格單元數(shù)大于1萬時，都能保持很高的并行效率。

圖18 并行效率隨每核網(wǎng)格單元數(shù)變化Fig.18 Efficiency versus numbers of grid elements

4 結(jié) 論

本文發(fā)展了基于分布式并行系統(tǒng)的大規(guī)模流動模擬方法，包括離散控制方程的緊致型數(shù)值格式、基于Metis的網(wǎng)格分區(qū)技術(shù)和并行邊界處理及虛擬單元技術(shù)等。數(shù)值結(jié)果表明，軟件并行計算與串行計算結(jié)果一致，并行效率高。

與串行計算相比，并行計算流動收斂過程有小量的差別，但收斂數(shù)值在合理的誤差范圍內(nèi)一致。在一定的并行核數(shù)范圍內(nèi)，表現(xiàn)出超線性的加速比特性。由于加速比性能好，所有的測試算例都具有很高的并行效率。另外，要保證有較高的并行效率，每核網(wǎng)格單元數(shù)不能太少。

總之，本文方法能有效地計算復(fù)雜飛行器構(gòu)型的氣動性能，并可充分利用大規(guī)模并行計算減少流動模擬時間和提高計算效率。