基于神威超算平臺的三維方柱繞流大規(guī)模并行數(shù)值計算與分析

張亞英，吳乘勝，王建春，顧寒鋒

（1.中國船舶科學(xué)研究中心，江蘇無錫 214082；2.國家超級計算無錫中心，江蘇無錫 214072）

0 引 言

由于“主頻墻”、“通信墻”和“功耗墻”等的限制，近年來超級計算機(jī)架構(gòu)正在發(fā)生由同構(gòu)向異構(gòu)發(fā)展的變革。以2020年11月發(fā)布的超級計算機(jī)TOP500[1]榜單為例，其中約30%采用了加速卡或協(xié)處理器，且近十年來登頂?shù)某売嬎銠C(jī)均采用異構(gòu)處理器或異構(gòu)節(jié)點(diǎn)，這種趨勢表明異構(gòu)已經(jīng)成為超級計算機(jī)的發(fā)展趨勢，并在性能上有顯著的優(yōu)勢。

異構(gòu)與同構(gòu)的差異在于：同構(gòu)架構(gòu)中所有的計算核心都由CPU構(gòu)成，所有計算核心的邏輯處理能力和數(shù)據(jù)計算能力都很強(qiáng)，不受計算任務(wù)復(fù)雜度的影響；缺點(diǎn)是成本高、功耗大。異構(gòu)與同構(gòu)相比，計算核心的種類不同，一般包含CPU和協(xié)處理器，或者將兩種不同功能的計算核心集成在單個芯片上形成異構(gòu)眾核處理器。異構(gòu)系統(tǒng)通過CPU 調(diào)度加速硬件或加速核可以實(shí)現(xiàn)深層次的并行，使計算任務(wù)劃分更加細(xì)化。典型的異構(gòu)模式包括CPU+GPU、CPU+MIC，以及我國第一臺全部采用自主技術(shù)構(gòu)建的世界第一的超級計算機(jī)——“神威·太湖之光”所采用的SW26010處理器。

異構(gòu)體系使超級計算機(jī)的計算能力大幅度提升，但同時對CFD 高性能計算也是挑戰(zhàn)。要實(shí)現(xiàn)對計算資源的充分利用，需要從數(shù)值算法、數(shù)據(jù)結(jié)構(gòu)、計算流程等各個層面進(jìn)行重構(gòu)和優(yōu)化。在面向神威超算平臺的CFD 并行計算方面，國內(nèi)不少研究人員開展了研究工作。Ren等[2]將OpenFOAM 移植到“神威·太湖之光”上并進(jìn)行了優(yōu)化，優(yōu)化后主從核協(xié)同使用，相較于純主核應(yīng)用加速了3.55 倍，單個SW26010 處理器較Intel x86 處理器加速了1.18 倍；對于CFD 的核心計算過程，倪鴻等[3]基于SW26010處理器架構(gòu)對非結(jié)構(gòu)網(wǎng)格下稀疏下三角方程求解進(jìn)行了研究，提出了一種基于流水線串行-局部并行思想的通用眾核優(yōu)化方法，實(shí)現(xiàn)了單核組3 倍以上的加速效果；劉侃等[4]基于SW26010 處理器設(shè)計了針對三對角矩陣求解的并行算法swDCR，通過寄存器通信設(shè)計及雙緩沖等，相較于主核加速約2.07倍；在CFD 適應(yīng)性方面，李芳等[5]對CFD 中的隱式算法、多重網(wǎng)格等問題的眾核并行難點(diǎn)進(jìn)行了分析，給出了相應(yīng)的并行方案，并對OpenCFD、SWLBM 以及AHL3D 等流體力學(xué)軟件的眾核適應(yīng)性進(jìn)行了討論。

總體上看，國內(nèi)CFD領(lǐng)域尤其是船舶水動力學(xué)CFD領(lǐng)域，目前遠(yuǎn)未實(shí)現(xiàn)對國產(chǎn)異構(gòu)眾核超算能力的有效利用。究其原因，除了CFD應(yīng)用研究長期依賴國外商用軟件之外，還包括以下客觀因素：

（1）CFD 計算通常具有全局相關(guān)性的特點(diǎn)，并行規(guī)模的增大帶來了并行復(fù)雜度與通訊開銷的增加，導(dǎo)致并行效率下降；同時水動力學(xué)CFD常用的SIMPLE算法的流程相對復(fù)雜，增加了細(xì)粒度并行優(yōu)化的難度，眾核加速挑戰(zhàn)很大；

（2）CFD 軟件一般具有數(shù)據(jù)結(jié)構(gòu)復(fù)雜、計算流程復(fù)雜和代碼量龐大等特點(diǎn)，從程序的移植到優(yōu)化，都需要大量的重構(gòu)工作，難度和工作量巨大。

相較于上述對大型程序和局部迭代的移植和測試，本文更具針對性地對求解不可壓流動問題的SIMPLE 算法，進(jìn)行面向國產(chǎn)異構(gòu)超算平臺的并行計算研究。論文基于神威異構(gòu)超算平臺，采用自主開發(fā)代碼，開展三維非定常不可壓流動的大規(guī)模并行數(shù)值模擬研究。數(shù)值計算采用MPI+Athread多級并行方式，并在Athread眾核并行中，針對申威眾核處理器的特點(diǎn)對SIMPLE算法求解流程進(jìn)行優(yōu)化以提升加速效果。首先以三維方柱準(zhǔn)定常繞流問題為例開展并行計算測試，包括MPI 并行和MPI+Athread 多級并行，證實(shí)多級并行的有效性；隨后，針對三維方柱非定常繞流問題進(jìn)行直接數(shù)值模擬，開展網(wǎng)格數(shù)量從384 萬到2.46 億的大規(guī)模并行計算并對計算結(jié)果進(jìn)行分析，展現(xiàn)神威異構(gòu)超算平臺在非定常不可壓縮流動CFD大規(guī)模并行計算方面的應(yīng)用能力。

1 數(shù)值模擬并行計算方案

1.1 數(shù)值模擬方法簡述

本文的主要目的是研究神威異構(gòu)超算系統(tǒng)對三維不可壓粘性流動CFD 核心求解過程的適應(yīng)性，其積分形式控制方程組如下：

式中，V為控制體體積，S為控制體表面積，U→為速度矢量，其分量形式表示為（u，v，w），υ為流體運(yùn)動粘性系數(shù)，U∞為無窮遠(yuǎn)處自由流速度。

采用基于交錯網(wǎng)格的有限體積法離散控制方程，具體可參閱文獻(xiàn)[6]；控制方程組的求解采用SIMPLE算法[6]。

1.2 申威眾核處理器簡介

“神威·太湖之光”超級計算機(jī)，采用的是國產(chǎn)申威架構(gòu)SW26010處理器。處理器本身就包括控制核心和計算核心陣列，相當(dāng)于把CPU和加速處理器集成到一個芯片上，其內(nèi)部架構(gòu)如圖1所示。

圖1 SW26010處理器架構(gòu)示意圖Fig.1 Diagram of SW26010 architecture

SW26010 處理器包含四個核組（CG），各核組之間采用片上網(wǎng)絡(luò)（NOC）互聯(lián)，每個核組內(nèi)包含一個主控制核心（主核，MPE）、1 個從核（CPE）集群（由64 個從核組成）、1 個協(xié)議處理單元（PPU）和1 個內(nèi)存控制器（MC）。核組內(nèi)采用共享存儲架構(gòu)，內(nèi)存與主/從核之間可通過MC傳輸數(shù)據(jù)。

SW26010處理器的計算能力主要體現(xiàn)在從核上，但從核上的存儲空間和帶寬較小，使得數(shù)據(jù)傳輸往往成為程序運(yùn)行的瓶頸，這也是并行方案設(shè)計中的重點(diǎn)。

1.3 并行計算方案

并行方案采用MPI+Athread的多級并行方式，其中：主核（核組）采用基于區(qū)域分解的MPI并行；從核并行則采用神威專用加速線程庫Athread，方式為循環(huán)分解；SIMPLE 算法的迭代過程采用SOR 算法，并使用紅黑排序進(jìn)行分塊，避免數(shù)據(jù)相關(guān)。

由于迭代過程具有全局特性，要保證計算結(jié)果的準(zhǔn)確性需要實(shí)時更新影像區(qū)；而實(shí)時更新通信代價巨大，通?？梢栽诒ＷC結(jié)果精度的前提下適當(dāng)減少通信次數(shù)。在本文中，由于內(nèi)迭代次數(shù)較少，因而采用實(shí)時更新的方法保證結(jié)果的準(zhǔn)確性。由此，在并行方案設(shè)計上，主要面臨兩方面的問題：SIMPLE算法迭代過程中的通信；從核有限內(nèi)存空間的高效利用。

由于SW26010 處理器每個從核有64 kB的高速存儲由程序員自己控制，從核訪問這一局部存儲空間（LDM）的速度非常快，而從核訪問主存則會有很高的延遲（見表1）。所以在從核并行中，通常將主存中的數(shù)據(jù)以直接內(nèi)存訪問（DMA）的形式將數(shù)據(jù)傳輸?shù)絃DM空間（如圖2所示）。

表1 SW26010處理器從核訪存性能Tab.1 CPE performance of SW26010 processor

圖2 循環(huán)融合前后變量有效性及通信變化Fig.2 Change of variable validity and communication before and after loop consolidation

CFD程序中，通常都會使用結(jié)構(gòu)體。由于結(jié)構(gòu)體中一般不可避免地存在“非相關(guān)變量”，從而降低從核存儲空間的使用效率。如圖2 左圖所示，一個循環(huán)的眾核并行，其LDM 空間由相關(guān)變量（relevant variables）和非相關(guān)變量（irrelevant variables）構(gòu)成，結(jié)構(gòu)體內(nèi)“非相關(guān)變量”的傳輸既增加通信負(fù)擔(dān)又占用從核存儲空間；此外，過多的小循環(huán)也會增加通信頻次。這些都會導(dǎo)致眾核并行效率的降低。

為提高LDM 空間的使用和眾核并行效率，本文采用循環(huán)融合的方式，對SIMPLE 算法的求解流程進(jìn)行了優(yōu)化，以增加結(jié)構(gòu)體中變量的“有效性”，提高變量的復(fù)用度，避免同一變量的重復(fù)傳輸，如圖2 右圖所示。通過對二維方腔頂蓋驅(qū)動流的眾核并行測試，結(jié)果顯示優(yōu)化后的SIMPLE 算法求解流程，眾核并行加速比提高了約25%（如圖3所示）。

圖3 SIMPLE算法計算流程優(yōu)化前后眾核并行加速比Fig.3 Many-core parallel speedup for SIMPLE algorithm before and after loop consolidation

2 并行計算測試

2.1 算例與計算結(jié)果

三維不可壓流動問題的并行計算，以三維長方型截面柱體繞流為算例模型開展應(yīng)用測試（如圖4所示），表2 中為相關(guān)計算及模型參數(shù)，其中，h，l，m為柱體尺寸，并取d為特征長度。A，B，C表示計算域尺寸，D表示速度入口到柱體前壁面的距離，U為來流速度。設(shè)置Re=10，此時為準(zhǔn)定常流動狀態(tài)，便于對計算結(jié)果的準(zhǔn)確性進(jìn)行驗(yàn)證。

表2 三維方柱體繞流多級并行計算測試算例Tab.2 Test case of flow past 3D square cylinder for multilevel parallel computation

圖4 三維方柱體繞流算例示意圖Fig.4 Test case of flow past 3D square cylinder

圖5 為流動狀態(tài)達(dá)到穩(wěn)定后的流線分布以及x方向速度云圖，在中縱截面（y=1.5d）上柱體后方的渦心距離柱體后壁面約0.3d，高度約為0.25d，渦長度約為0.65d，與文獻(xiàn)[7]結(jié)果一致。

圖5 Re=10時三維矩形柱體繞流流線計算結(jié)果Fig.5 Computational result of streamline of flow past 3D square cylinder with Re=10

2.2 并行加速效果分析

分別對MPI 和MPI+Athread 的并行加速效果進(jìn)行測試分析。MPI 并行測試主要分析其并行效率，考察其擴(kuò)展性。不同網(wǎng)格數(shù)量算例下的MPI 并行加速比和并行效率分別見圖6 和表3。并行加速比和并行效率計算分析中，由于單核內(nèi)存無法處理大網(wǎng)格量算例，因而以2進(jìn)程（2核）作為對比基準(zhǔn)。

表3 不同網(wǎng)格數(shù)量算例MPI并行效率Tab.3 MPI parallel efficiency for cases with different grid numbers

圖6 不同網(wǎng)格數(shù)量算例MPI并行加速比Fig.6 MPI parallel speedup for cases with different grid numbers

從圖表中可以看出：MPI并行加速比和并行效率都隨著算例網(wǎng)格數(shù)量的增加而提高；相同網(wǎng)格規(guī)模下，MPI并行效率隨著進(jìn)程數(shù)（核數(shù)）增加而降低；對2025 萬網(wǎng)格算例，125 進(jìn)程并行計算加速約45倍（相對于2進(jìn)程），并行效率約為72.5%。

為了保證計算精度，在SIMPLE 算法迭代過程中，采用了實(shí)時通信的方式，一定程度上降低了MPI的并行效果。在實(shí)際的CFD 計算中，可以綜合考慮計算精度以及通信開銷，適當(dāng)減少通信頻次，能夠提高M(jìn)PI并行效率。

在MPI 并行計算測試的基礎(chǔ)上，開展了MPI+Athread 多級并行計算測試。不同網(wǎng)格數(shù)量算例下的多級并行加速比和眾核并行加速比分別見圖7和表4。

表4 不同網(wǎng)格數(shù)量算例眾核并行加速比Tab.4 Many-core parallel speedup for cases with different grid numbers

圖7 不同網(wǎng)格數(shù)量算例MPI+Athread多級并行加速比Fig.7 MPI+Athread parallelspeedup for cases with different grid numbers

從圖表中可以看出：MPI+Athread 多級并行和眾核并行加速比都隨著算例網(wǎng)格數(shù)量的增加而提高；對2025萬網(wǎng)格算例，125核組（主從核合計8125核）并行計算加速約174 倍（相對于2 主核），眾核（從核）并行加速7.9倍。

通過上述算例的測試與分析可見，對于三維不可壓流動的非定常數(shù)值模擬，無論是MPI并行還是MPI++Athread多級并行，都能夠獲得相當(dāng)可觀的加速效果，說明了本文的并行計算方案在神威超算平臺上的有效性。

3 三維方柱繞流大規(guī)模并行數(shù)值模擬與分析

鑒于前文中MPI+Athread的多級并行計算使并行加速比獲得顯著提升，將其用于大規(guī)模網(wǎng)格的非定常計算或能夠顯著提升計算速度。因此，采用直接數(shù)值模擬的方法，在神威異構(gòu)超算平臺上開展三維方柱非定常繞流數(shù)值模擬，柱體由圖4 中長方形截面改為方形截面，截面邊長同樣為d且為特征長度，計算域擴(kuò)大至A×B×C=24d×16d×10d。

Sakamoto 和Arie 等[8]研究發(fā)現(xiàn)，當(dāng)方柱高寬比h/d＞2 時，尾跡會形成交替脫落的渦即卡門渦街，反之則會形成對稱脫落的拱形渦。本文針對h/d=4的情況，開展Re=250三維方柱繞流的CFD模擬，對瞬時流場、渦系結(jié)構(gòu)、時均流場等進(jìn)行分析，并與相關(guān)文獻(xiàn)結(jié)果進(jìn)行對比驗(yàn)證。

數(shù)值計算中，計算域大小為24d×16d×10d，方柱位于距入口6d處；數(shù)值計算采用了三套網(wǎng)格，網(wǎng)格單元數(shù)量分別為384萬、3072萬和24 576萬（2.457 6億）；數(shù)值計算的網(wǎng)格劃分、時間步長、計算資源和計算耗時等相關(guān)參數(shù)列于表5中。

表5 大規(guī)模并行數(shù)值模擬相關(guān)參數(shù)Tab.5 Parameters of massive parallel CFD computation

3.1 網(wǎng)格影響分析

雖然網(wǎng)格單元數(shù)量384 萬已基本滿足Re=250 三維方柱繞流的直接數(shù)值模擬要求，本文仍然通過三維倍增細(xì)化的方式（見表5），從宏觀渦系結(jié)構(gòu)、典型位置的時均速度分布和橫向速度時歷頻譜分析等多個層面，進(jìn)行數(shù)值模擬結(jié)果的網(wǎng)格影響分析。

（1）宏觀渦系結(jié)構(gòu)

圖8 給出了根據(jù)Q判據(jù)（Q=0.01）得到的Re=250 三維方柱繞流渦系結(jié)構(gòu)，其中圖8（a）為384 萬網(wǎng)格模擬結(jié)果，圖8（b）為3072 萬網(wǎng)格模擬結(jié)果，圖8（c）為2.457 6 億網(wǎng)格模擬結(jié)果。從圖中可以看出三維方柱非定常繞流的典型特征：方柱與底面結(jié)合部前方形成馬蹄渦；方柱自由端有大包絡(luò)面的形成，表征自由端剪切層的分離；在方柱后方尾流區(qū)，形成了反對稱的卡門渦。從圖中同時可以看出，隨著網(wǎng)格規(guī)模增加（空間分辨率增強(qiáng)），數(shù)值模擬獲得的渦系結(jié)構(gòu)更為豐富、精細(xì)。

圖8 Re=250三維方柱繞流渦系結(jié)構(gòu)模擬結(jié)果Fig.8 Vortex system structure of flow past 3D square cylinder with Re=250

（2）時均速度的空間分布

為更好地分析網(wǎng)格空間分辨率對非定常三維方柱繞流時均流場的影響，取計算域中如圖9 所示的a、b、c三條直線，分別比較線上的流向、橫向、垂向時間平均速度分布。三條線的位置分別是：直線a平行于x軸，位于計算域的中縱截面上，距底面2.0d；直線b平行于y軸，位于距離方柱中心7.0d的垂直截面上，距底面3.5d；直線c則為距離方柱中心7.0d的垂直截面的垂直中心線。

圖9 時均速度對比分析位置示意圖Fig.9 Locations of time-averaged velocities for caparison

圖10～12 分別給出了a、b、c三條線上三套網(wǎng)格的流向、橫向、垂向時均速度分布計算結(jié)果，圖中同時給出了文獻(xiàn)[9]的結(jié)果。從圖中可以看出：對于三套網(wǎng)格，流向、橫向和垂向時均速度計算結(jié)果均體現(xiàn)出與文獻(xiàn)結(jié)果相同的空間分布特征；隨著網(wǎng)格數(shù)量增加，橫向和垂向時均速度隨空間變化的波動幅度越大；3072 萬網(wǎng)格和2.457 6 億網(wǎng)格的計算結(jié)果與文獻(xiàn)結(jié)果都較為接近，總體上3072萬網(wǎng)格的計算結(jié)果與文獻(xiàn)更為接近。

圖10 直線a上流向時均速度計算結(jié)果Fig.10 Computational results of time-averaged stream-wise velocity on Line a

由此可見，在滿足直接數(shù)值模擬要求的前提下，網(wǎng)格數(shù)量增加一方面會使數(shù)值模擬的宏觀渦系結(jié)構(gòu)更為精細(xì)，同時對時均速度的空間分布計算結(jié)果也有一定的影響。

圖11 直線b上橫向時均速度計算結(jié)果Fig.11 Computational results of time-averaged transvers velocity on Line b

圖12 直線c上垂向時均速度計算結(jié)果Fig.12 Computational results of time-averaged vertical velocity on Line c

（3）橫向速度時歷頻譜分析

前面提到，本文研究的三維方柱繞流工況，尾跡會形成交替脫落的渦即卡門渦街，必然會導(dǎo)致尾跡中橫向速度的振蕩。為分析網(wǎng)格數(shù)量對方柱尾跡卡門渦街這一典型非定常流動特征數(shù)值模擬結(jié)果的影響，取方柱正后方A、B兩個位置（如圖13所示），對橫向速度時間歷程使用快速傅里葉分析（FFT）進(jìn)行頻譜分析和比較。A、B兩點(diǎn)至方柱中心的距離分別為8.5d和15.5d，距底面皆為2.0d，都位于渦脫落影響區(qū)域內(nèi)。

圖13 橫向速度時歷頻譜分析位置示意圖Fig.13 Locations of time-history of transverse velocities for spectrum analysis

圖14 給出了不同網(wǎng)格數(shù)量下方柱后A、B兩點(diǎn)的橫向速度時間歷程數(shù)值模擬結(jié)果及其能量譜分析，其主頻與幅值列于表6中。

圖14 方柱后方不同位置橫向速度時歷曲線及能量譜分析Fig.14 Time history and spectrum of transverse velocity at different locations behind the cylinder

表6 方柱后不同位置橫向速度時歷頻譜分析結(jié)果Tab.6 Spectrum analysis of transverse velocity at different locations behind the cylinder

從圖表中可以看出：隨著網(wǎng)格數(shù)量增加，數(shù)值模擬得到的橫向速度振蕩主頻略有增大，且網(wǎng)格數(shù)量大于3072 萬后基本保持不變；而能量譜主頻處峰值則隨網(wǎng)格數(shù)量增加有所增大，且高頻成分明顯增強(qiáng)；此外，隨著流向距離的增加，能量譜的峰值有所降低，反映了因粘性導(dǎo)致的能量耗散和衰減。

結(jié)合前面宏觀渦系結(jié)構(gòu)對比分析結(jié)果，說明本文使用的三套數(shù)量不同的網(wǎng)格，一方面都能夠準(zhǔn)確地模擬出三維方柱非定常繞流的渦脫現(xiàn)象和過程，另一方面網(wǎng)格數(shù)量增加能夠更好地捕捉流場細(xì)節(jié)，特別是非定常流動中的高頻成分。

3.2 流場模擬結(jié)果分析

綜合以上網(wǎng)格數(shù)量對數(shù)值模擬結(jié)果影響的分析，接下來對三維方柱非定常繞流流場進(jìn)行分析，重點(diǎn)針對3072萬網(wǎng)格的數(shù)值模擬結(jié)果。以下從瞬時流場和時均流場兩個方面，對Re=250三維方柱非定常繞流流場的數(shù)值模擬結(jié)果進(jìn)行分析。

（1）瞬時流場分析

圖15 給出了t=450 s 時刻以流向渦分量表征的渦系結(jié)構(gòu)，與圖8 中根據(jù)Q判據(jù)得到的渦系結(jié)構(gòu)相比，能夠更為清晰地看出方柱頂端的梢渦。由于方柱根部渦系結(jié)構(gòu)復(fù)雜，在圖8 和圖15 中都難以觀察到底層渦結(jié)構(gòu)，不過在圖17 的方柱后壁面流線圖中，可以清楚地看到其根部由于底層渦形成的兩個對稱渦旋。

圖15 瞬時x方向渦分量等值面渦系結(jié)構(gòu)Fig.15 Instantaneous vortex system structure for x_vorticy

圖16 為t=450 s 時刻計算域中縱截面（y/d=8.0）和水平截面（z/d=1.0）上的流線計算結(jié)果。從圖中可以看出：水平截面上有明顯的非對稱渦的形成與脫落，與Sakamoto、Arie 等人的研究結(jié)果相符合；中縱截面上，在柱體前方，流動在此處分離，發(fā)生上沖和下洗；其中下洗流體在方柱與底面結(jié)合處回流，形成馬蹄渦；上沖流體經(jīng)方柱自由端發(fā)生分離，部分流體回流，在方柱后形成較大的回流區(qū)，延伸至方柱自由端處；同時在底部形成了反方向的回流區(qū)，且在回流下游靠近底面部分存在流線分離點(diǎn)即鞍點(diǎn)。

圖16 Re=250三維方柱繞流流線計算結(jié)果Fig.16 Computational results of streamline of flow past 3D square cylinder with Re=250

圖17為方柱前、后及側(cè)壁的極限流線。從圖中可以看出：在柱體前壁面有明顯的滯流現(xiàn)象，流體在此處流向發(fā)生偏轉(zhuǎn)并分離；在后壁面上，流動相較于前壁面更加復(fù)雜，可以看到存在一個奇點(diǎn)，流線均從該點(diǎn)發(fā)出，結(jié)合圖16 所示中縱截面上的流線分布，后壁面上下流線的分離則表示方柱和后方形成的兩個方向相反的回流區(qū)，并根據(jù)奇點(diǎn)位置，表明上升流仍然起主要作用，同時在底部可以看到一對底層渦；在側(cè)壁面上有兩組流線在此處相遇，且兩組流線均存在上下分離現(xiàn)象，向下流動的流線進(jìn)行匯聚，形成一個點(diǎn)，反映了底層渦的形成；向上流動的流體則匯聚成一條直線。

圖17 方柱前、后、側(cè)面流線計算結(jié)果Fig.17 Computational results of streamline on front,back and side faces of the square cylinder

圖18為不同高度水平截面上Q準(zhǔn)則渦量等值線，高度從0.2d到4.3d。圖18（a）為底部附近截面上的渦量等值線，由于非滑移壁面的存在，底部的渦系結(jié)構(gòu)非常復(fù)雜，但是可以明顯地觀察到柱體前方的馬蹄渦；圖18（b）中的渦量等值線則呈現(xiàn)出交替式的渦脫，但已不再是標(biāo)準(zhǔn)的卡門渦街分布，這表明自由端已經(jīng)對中部的渦旋產(chǎn)生了影響；圖18（c）-（e）則是方柱上部分以及頂面附近的渦量等值線圖，可以看到隨著高度向自由端靠近，卡門渦街越來越不明顯，但是仍能觀察到渦脫現(xiàn)象；同時在頂部附近可以觀察到成對的梢渦，并向上傾斜，與圖15中所示的梢渦形態(tài)相同。

圖18 水平截面上Q準(zhǔn)則渦量等值線圖Fig.18 Q criterion vorticity contour on horizontal section

（2）時均流場分析

相較于瞬時流場的隨機(jī)性，時均流場能夠相對定量地反映流場中的一些典型特征。圖19和圖20分別給出了中縱截面（y/d=8.0）和水平截面（z/d=1.0）上時均流場的流線分布。

圖19 Re=250三維方柱繞流中縱截面時均流線（y/d=8.0）Fig.19 Time-averaged streamline for flow past 3D square cylinder of Re=250 on y/d=8.0

圖20 Re=250三維方柱繞流中橫截面時均流線（z/d=1.0）Fig.20 Time-averaged streamline for flow past 3D square cylinder of Re=250 on z/d=1.0

表7 則給出了上部回流區(qū)和根部回流區(qū)特征信息及其與文獻(xiàn)結(jié)果的對比。由圖表可以看出：圖19 所示的水平截面上流線分布圖中，柱體后方流線均從距離柱體中心約6d的奇點(diǎn)流出，并呈現(xiàn)對稱狀態(tài)，該奇點(diǎn)對應(yīng)圖18中的方柱后上沖和下洗的分離點(diǎn)鞍點(diǎn)“+”，奇點(diǎn)位置以及流線分布與文獻(xiàn)[9]相當(dāng)接近。

表7 方柱后回流區(qū)大小及位置信息Tab.7 Size and position of reflux zone behind the square cylinder

4 結(jié) 論

本文針對不可壓縮流動常用的CFD 求解方法——SIMPLE 算法，開展了MPI+Athread 多級并行計算研究，并基于神威·太湖之光異構(gòu)超算平臺進(jìn)行了三維方柱繞流的大規(guī)模并行數(shù)值模擬，最大網(wǎng)格數(shù)量達(dá)到2.46億，并行規(guī)模達(dá)到13.3萬核。通過對數(shù)值模擬結(jié)果的分析，可以得出以下結(jié)論：

（1）基于計算域區(qū)域分解和子區(qū)域內(nèi)部循環(huán)分解的MPI+Athread 多級并行模式，能夠使申威超算平臺眾核處理器的從核發(fā)揮較為可觀的加速效果，加速比達(dá)到7.5～7.9 倍，且網(wǎng)格量越大加速效果越明顯；

（2）在三維有限長方柱繞流的直接數(shù)值模擬中，采用MPI+Athread 多級并行，網(wǎng)格數(shù)量增加到2.46 億，并行規(guī)模擴(kuò)展到13.3 萬核，可以將計算周期縮短至約160 小時，從而大幅度提升了CFD 大規(guī)模并行計算的能力；

（3）對于長徑比為4、Re=250 的三維方柱非定常繞流，采用網(wǎng)格數(shù)量從384 萬到2.46 億進(jìn)行直接數(shù)值模擬，對渦的脫落過程等宏觀流動現(xiàn)象的模擬結(jié)果沒有明顯區(qū)別，不過網(wǎng)格數(shù)量增加能夠更好地捕捉流場細(xì)節(jié)，模擬得到的渦系結(jié)構(gòu)更為精細(xì)。

需要說明的是，本文的研究工作還是初步的，無論是MPI并行、眾核并行還是MPI+Athread多級并行，都還存在優(yōu)化和提升的空間；同時，本文的多級并行技術(shù)要應(yīng)用到大型的CFD 計算程序/軟件并達(dá)到實(shí)用化的程度，還有很多的研究和開發(fā)工作需要去做。這些也是論文研究團(tuán)隊(duì)接下來的重點(diǎn)工作。