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

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

0 引 言

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

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

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

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

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

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

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

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

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

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

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

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

1.2 申威眾核處理器簡(jiǎn)介

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

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

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

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

1.3 并行計(jì)算方案

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

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

由于SW26010 處理器每個(gè)從核有64 kB的高速存儲(chǔ)由程序員自己控制，從核訪問(wèn)這一局部存儲(chǔ)空間（LDM）的速度非?？?，而從核訪問(wèn)主存則會(huì)有很高的延遲（見(jiàn)表1）。所以在從核并行中，通常將主存中的數(shù)據(jù)以直接內(nèi)存訪問(wèn)（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程序中，通常都會(huì)使用結(jié)構(gòu)體。由于結(jié)構(gòu)體中一般不可避免地存在“非相關(guān)變量”，從而降低從核存儲(chǔ)空間的使用效率。如圖2 左圖所示，一個(gè)循環(huán)的眾核并行，其LDM 空間由相關(guān)變量（relevant variables）和非相關(guān)變量（irrelevant variables）構(gòu)成，結(jié)構(gòu)體內(nèi)“非相關(guān)變量”的傳輸既增加通信負(fù)擔(dān)又占用從核存儲(chǔ)空間；此外，過(guò)多的小循環(huán)也會(huì)增加通信頻次。這些都會(huì)導(dǎo)致眾核并行效率的降低。

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

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

2 并行計(jì)算測(cè)試

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

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

表2 三維方柱體繞流多級(jí)并行計(jì)算測(cè)試算例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 為流動(dòng)狀態(tài)達(dá)到穩(wěn)定后的流線分布以及x方向速度云圖，在中縱截面（y=1.5d）上柱體后方的渦心距離柱體后壁面約0.3d，高度約為0.25d，渦長(zhǎng)度約為0.65d，與文獻(xiàn)[7]結(jié)果一致。

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

2.2 并行加速效果分析

分別對(duì)MPI 和MPI+Athread 的并行加速效果進(jìn)行測(cè)試分析。MPI 并行測(cè)試主要分析其并行效率，考察其擴(kuò)展性。不同網(wǎng)格數(shù)量算例下的MPI 并行加速比和并行效率分別見(jiàn)圖6 和表3。并行加速比和并行效率計(jì)算分析中，由于單核內(nèi)存無(wú)法處理大網(wǎng)格量算例，因而以2進(jìn)程（2核）作為對(duì)比基準(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ù)）增加而降低；對(duì)2025 萬(wàn)網(wǎng)格算例，125 進(jìn)程并行計(jì)算加速約45倍（相對(duì)于2進(jìn)程），并行效率約為72.5%。

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

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

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

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

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

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

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

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

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

數(shù)值計(jì)算中，計(jì)算域大小為24d×16d×10d，方柱位于距入口6d處；數(shù)值計(jì)算采用了三套網(wǎng)格，網(wǎng)格單元數(shù)量分別為384萬(wàn)、3072萬(wàn)和24 576萬(wàn)（2.457 6億）；數(shù)值計(jì)算的網(wǎng)格劃分、時(shí)間步長(zhǎng)、計(jì)算資源和計(jì)算耗時(shí)等相關(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 萬(wàn)已基本滿足Re=250 三維方柱繞流的直接數(shù)值模擬要求，本文仍然通過(guò)三維倍增細(xì)化的方式（見(jiàn)表5），從宏觀渦系結(jié)構(gòu)、典型位置的時(shí)均速度分布和橫向速度時(shí)歷頻譜分析等多個(gè)層面，進(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àn)網(wǎng)格模擬結(jié)果，圖8（b）為3072 萬(wàn)網(wǎng)格模擬結(jié)果，圖8（c）為2.457 6 億網(wǎng)格模擬結(jié)果。從圖中可以看出三維方柱非定常繞流的典型特征：方柱與底面結(jié)合部前方形成馬蹄渦；方柱自由端有大包絡(luò)面的形成，表征自由端剪切層的分離；在方柱后方尾流區(qū)，形成了反對(duì)稱的卡門渦。從圖中同時(shí)可以看出，隨著網(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）時(shí)均速度的空間分布

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

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

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

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

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

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

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

（3）橫向速度時(shí)歷頻譜分析

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

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

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

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

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

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

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

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

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

（1）瞬時(shí)流場(chǎng)分析

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

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

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

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

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

圖17 方柱前、后、側(cè)面流線計(jì)算結(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)對(duì)中部的渦旋產(chǎn)生了影響；圖18（c）-（e）則是方柱上部分以及頂面附近的渦量等值線圖，可以看到隨著高度向自由端靠近，卡門渦街越來(lái)越不明顯，但是仍能觀察到渦脫現(xiàn)象；同時(shí)在頂部附近可以觀察到成對(duì)的梢渦，并向上傾斜，與圖15中所示的梢渦形態(tài)相同。

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

（2）時(shí)均流場(chǎng)分析

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

圖19 Re=250三維方柱繞流中縱截面時(shí)均流線（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三維方柱繞流中橫截面時(shí)均流線（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é)果的對(duì)比。由圖表可以看出：圖19 所示的水平截面上流線分布圖中，柱體后方流線均從距離柱體中心約6d的奇點(diǎn)流出，并呈現(xiàn)對(duì)稱狀態(tài)，該奇點(diǎn)對(duì)應(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é) 論

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

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

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

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

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