基于“神威·太湖之光”的三維有限長方柱繞流直接數(shù)值模擬

張亞英，吳乘勝，王建春，金奕星

中國船舶科學(xué)研究中心，江蘇 無錫 214082

0 引 言

隨著船舶工業(yè)技術(shù)和船舶水動(dòng)力學(xué)的發(fā)展，船舶設(shè)計(jì)及理論研究對數(shù)值模擬的精度要求越來越高，幾何模型也越來越復(fù)雜，導(dǎo)致計(jì)算所用網(wǎng)格規(guī)模呈現(xiàn)幾何量級增長。此外，隨著摩爾定律的失效，單核處理器的性能趨于極限，已無法滿足科研、設(shè)計(jì)工作的需求，并行計(jì)算已然成為計(jì)算流體力學(xué)（CFD）領(lǐng)域的關(guān)鍵技術(shù)之一。

目前，得益于計(jì)算硬件快速發(fā)展，超級計(jì)算機(jī)的體系結(jié)構(gòu)正在發(fā)生從同構(gòu)向異構(gòu)方向的轉(zhuǎn)變。相比于同構(gòu)體系的分布式并行計(jì)算，異構(gòu)系統(tǒng)通過在單個(gè)進(jìn)程內(nèi)使用圖形處理器（graphics processing unit，GPU）、眾核集成處理器（many integrated core，MIC）以及從核陣列等加速設(shè)備，實(shí)現(xiàn)了更細(xì)粒度的并行，突破了以往單核主頻對硬件算力的限制，使超級計(jì)算機(jī)在計(jì)算性能方面有了顯著提升。與此同時(shí)，更細(xì)粒度的并行將導(dǎo)致多層次的數(shù)據(jù)分配，使通信過程更加復(fù)雜。尤其是在CFD 領(lǐng)域，一方面，由于方程和算法與高性能計(jì)算的適配程度有限，對于不可壓縮流動(dòng)，控制方程本身存在速度壓強(qiáng)耦合的現(xiàn)象，隱式和半隱式算法所導(dǎo)致的數(shù)據(jù)相關(guān)性給多級并行帶來了較大的阻礙；另一方面，因加速設(shè)備本身相較于CPU 在內(nèi)存或指令處理能力上的不足，需在數(shù)據(jù)結(jié)構(gòu)及求解流程層面依據(jù)加速設(shè)備做出調(diào)整，導(dǎo)致工作量變得龐大且難度較高。所以，異構(gòu)超算平臺(tái)的應(yīng)用對船舶CFD 領(lǐng)域的并行計(jì)算既是機(jī)遇也是挑戰(zhàn)。

近年來，隨著船舶水動(dòng)力學(xué)CAE 軟件的自主化要求越來越高，自主軟硬件的結(jié)合是實(shí)現(xiàn)全面自主可控的必經(jīng)之路。在上述背景下，本文將基于“神威·太湖之光”異構(gòu)超算平臺(tái)，運(yùn)用MPI+Athread 的混合編程方法對自主代碼進(jìn)行并行處理，并以三維有限長方柱繞流為算例，使用SIMPLE（semi-implicit method for pressure linked equations）算法對繞流流動(dòng)進(jìn)行直接數(shù)值模擬。

1 并行方案

“神威·太湖之光”是我國自主研發(fā)的世界首臺(tái)峰值運(yùn)算速度超過十億億次的超級計(jì)算機(jī)，在2016～2017 年連續(xù)位居計(jì)算機(jī)性能世界500 強(qiáng)榜單的第1 位[1]。該計(jì)算系統(tǒng)基于SW26010 異構(gòu)眾核處理器構(gòu)建，可通過主核、從核實(shí)現(xiàn)異構(gòu)形式的多級并行計(jì)算，目前已應(yīng)用于地球數(shù)值模擬、分子模擬、矩陣求解[2]等領(lǐng)域。“神威·太湖之光”超級計(jì)算機(jī)擁有輔助計(jì)算系統(tǒng)和高速計(jì)算系統(tǒng)。前者采用x86 架構(gòu)處理器，后者采用我國自主研發(fā)的SW26010 處理器。如圖1 所示，該處理器共包含4 個(gè)核組，每個(gè)核組由一個(gè)主核和64 個(gè)從核組成，從核的并行由主核使用加速線程庫Athread[3]進(jìn)行控制，以實(shí)現(xiàn)從核的開啟、關(guān)閉和數(shù)據(jù)傳輸?shù)炔僮鳌W26010 處理器具有特殊的結(jié)構(gòu)，可以在核組間和從核間分別開展進(jìn)程級及線程級的并行計(jì)算。本文所述進(jìn)程級并行采用的是網(wǎng)格分區(qū)方法，按照進(jìn)程數(shù)對網(wǎng)格進(jìn)行分割，并將其分配到不同核組上進(jìn)行計(jì)算，使用消息傳遞接口（message-passing interface，MPI）實(shí)現(xiàn)進(jìn)程間的數(shù)據(jù)交換。由于從核內(nèi)存?。?4 kB）且通信只能在同行或同列間進(jìn)行，因此，線程級采用對循環(huán)并行方法將數(shù)據(jù)分批傳輸至從核上進(jìn)行計(jì)算。圖2所示為并行過程中計(jì)算數(shù)據(jù)的分配模式。

圖1 SW26010 處理器結(jié)構(gòu)Fig. 1 The structure of SW26010 processor

圖2 并行計(jì)算任務(wù)分配示意圖Fig. 2 Schematics of data allocation in parallel computing

為檢驗(yàn)上述并行方法的可行性，采用SIMPLE算法及直接數(shù)值模擬方法的并行效果，本文以三維頂板驅(qū)動(dòng)方腔流為模型，將雷諾數(shù)設(shè)置為Re=1 000，對不同網(wǎng)格規(guī)模及不同進(jìn)程數(shù)下的并行加速效果予以對比。

為表現(xiàn)多級并行的特點(diǎn)，本文展示的并行效果包括了MPI 單層次并行及MPI+Athread 多級并行的加速比及并行效率。其中，加速比計(jì)算參照對象均為單個(gè)主核計(jì)算耗時(shí)，考慮到主、從核在功能及內(nèi)存上的差異，MPI 并行效率計(jì)算以單主核為參照，MPI+Athread 并行效率則以單核組為參照。圖3 分別給出了MPI 及MPI+Athread 兩種并行模式在不同網(wǎng)格規(guī)模下的并行加速比變化曲線，表1 所示為并行效率。

圖3 采用SIMPLE 算法的3D 直接數(shù)值模擬加速比變化曲線Fig. 3 The speedup curves of direct 3D numerical simulation using SIMPLE algorithm

表1 3D 頂板驅(qū)動(dòng)方腔流MPI / MPI+Athread 并行效率Table 1 The parallel efficiency of MPI/MPI+Athread for 3D cavity driven flow

根據(jù)表1 所示結(jié)果，MPI 并行效率在128 進(jìn)程時(shí)均高于50%；而MPI+Athread 基于核組的擴(kuò)展效率偏低。結(jié)合文獻(xiàn)[4]對SW26010 處理器的從核并行加速影響因素的研究，經(jīng)分析發(fā)現(xiàn)導(dǎo)致上述現(xiàn)象的原因包括：

1) 從核計(jì)算的通信過程使整個(gè)通信耗時(shí)增加，降低了并行效率。

2) 隨著進(jìn)程數(shù)的增加，單個(gè)進(jìn)程計(jì)算量減少。根據(jù)文獻(xiàn)[4] 所述結(jié)論，此時(shí)將造成計(jì)算耗時(shí)比例下降，通信耗時(shí)比例提升，導(dǎo)致從核并行計(jì)算加速比明顯降低。

對比加速比變化的趨勢，在15.625 百萬網(wǎng)格數(shù)時(shí)，128 進(jìn)程多級并行加速比達(dá)到483.84 倍，且呈現(xiàn)出隨網(wǎng)格規(guī)模的增加，加速效果愈發(fā)顯著的變化趨勢。因此，采用圖2 所示并行方案，基于SIMPLE 算法使用億級網(wǎng)格規(guī)模對三維有限長方柱繞流進(jìn)行直接數(shù)值模擬確實(shí)可行，且能夠在有限的硬件資源條件下實(shí)現(xiàn)快速計(jì)算。

2 計(jì)算模型與控制方程

本文主要研究對象為三維有限長方柱繞流流動(dòng)現(xiàn)象，該現(xiàn)象普遍存在于高層建筑、海洋工程和交通工程等領(lǐng)域。以往的研究主要針對二維方柱繞流現(xiàn)象展開。但在工程實(shí)際中，通常為一端固定、另一端自由的有限長方柱，如圖4 所示。圖中：H為方柱高；橫截面為正方形，d為其寬度；u∞為來流速度。研究表明，方柱的長徑比H/d會(huì)對尾流有顯著影響。Sakamoto 和 Arie[5]研究發(fā)現(xiàn)方柱長徑比存在臨界值Hc，當(dāng)H/dHc時(shí)的三維有限長方柱繞流渦結(jié)構(gòu)模型。該模型顯示此時(shí)尾流為反對稱的卡門渦。

圖4 有限長方柱繞流3D 計(jì)算模型Fig. 4 The 3D model of flow around a finite square cylinder

圖5 有限長方柱繞流尾流渦3D 結(jié)構(gòu)[6]Fig. 5 The wake vortex structure of flow around a 3D finite square cylinder[6]

在本算例中，流動(dòng)介質(zhì)為不可壓黏性流體。根據(jù)圖4 選取方柱橫截面寬度d為特征長度，無窮遠(yuǎn)處速度為特征速度，積分形式的控制方程組如下：

式中：U為速度矢量；V為網(wǎng)格體積;S是網(wǎng)格面積；n為面外法向矢量；u，v，w為速度分量；P為壓力與密度之比；t為時(shí)間；υ 為流體運(yùn)動(dòng)黏性系數(shù)。

采用基于交錯(cuò)網(wǎng)格的有限體積法（finite volume method，F(xiàn)VM）離散控制方程。時(shí)間項(xiàng)離散采用一階顯式格式，對流項(xiàng)采用二階迎風(fēng)插值格式（QUICK），擴(kuò)散項(xiàng)采用中心差分格式，具體可見文獻(xiàn)[7]?？刂品匠探M的求解采用SIMPLE算法。

工況設(shè)定：雷諾數(shù)Re=u∞d/υ=250，邊界條件滿足左側(cè)為速度入口，來流速度為u∞=1.0；右側(cè)為壓力出口P=0；底面為固體壁面，滿足無滑移邊界條件，其他各面滿足對稱邊界條件。

3 網(wǎng)格規(guī)模效應(yīng)分析

相較于其他方法，直接數(shù)值模擬的優(yōu)勢在于其能夠獲取更完整的流場信息，計(jì)算域和網(wǎng)格分辨率需分別滿足積分尺度及耗散尺度的要求：

式中：L為計(jì)算域尺寸；l0為積分尺度； ?為網(wǎng)格尺寸；η 為耗散尺度。其中，積分尺度與計(jì)算域尺度同量級。

考慮到不同網(wǎng)格分辨率對于流場細(xì)節(jié)的捕捉能力有所不同，本文按照三向加密方式，縮小網(wǎng)格尺寸，以提升網(wǎng)格分辨率。一方面，統(tǒng)計(jì)計(jì)算耗時(shí)來表現(xiàn)并行計(jì)算在直接數(shù)值模擬中的可行性與優(yōu)越性；另一方面，對比長徑比H/d=4 時(shí)的流場差異，根據(jù)渦系結(jié)構(gòu)的豐富度及能譜特征來體現(xiàn)不同網(wǎng)格尺寸對小尺度流動(dòng)的捕捉能力。

表2 給出了不同網(wǎng)格規(guī)模下的網(wǎng)格分辨率、流場參數(shù)、并行規(guī)模及計(jì)算耗時(shí)。SW26010 處理器包含主、從核，使用的最大進(jìn)程數(shù)為2 048，核數(shù)達(dá)到了133 120。由表2 所示耗時(shí)可見，在時(shí)間推進(jìn)600 s（時(shí)間迭代步數(shù)6×105）情況下，億級網(wǎng)格規(guī)模（245.76 百萬網(wǎng)格數(shù)）的計(jì)算可在一周內(nèi)完成，千萬網(wǎng)格規(guī)模(30.72 百萬網(wǎng)格數(shù))的計(jì)算耗時(shí)縮短至3 天左右，而百萬網(wǎng)格規(guī)模（3.84 百萬網(wǎng)格數(shù)）計(jì)算耗時(shí)不足1 天。這表明在當(dāng)前有限的科研周期內(nèi)，使用“神威·太湖之光”超級計(jì)算機(jī)進(jìn)行CFD 數(shù)值計(jì)算，至少可將網(wǎng)格規(guī)模提升1～2 個(gè)量級。

使用Q判據(jù)[8]顯示方柱周圍及尾流中的瞬時(shí)渦系結(jié)構(gòu)，其形式如下：

圖6 所示為t=450 s 時(shí)刻瞬時(shí)流場渦系結(jié)構(gòu)（Q=0.01）。由圖可明顯觀察到方柱前方馬蹄形渦、方柱后包絡(luò)面的形成以及尾流中的反對稱渦，這與圖5 所示的渦結(jié)構(gòu)特征一致。對于不同網(wǎng)格規(guī)模下的計(jì)算結(jié)果，通過對比可以明顯看到，隨著網(wǎng)格分辨率的提高，渦系結(jié)構(gòu)更加清晰、豐富。這表明高分辨率網(wǎng)格能夠捕捉到更多的流動(dòng)細(xì)節(jié)。

圖6 t=450 s 時(shí)刻流場瞬時(shí)渦系結(jié)構(gòu)（Q=0.01）Fig. 6 Instantaneous vortex structures at t=450 s when Q=0.01

選取方柱正后方（z=2d，y=8d）兩個(gè)測點(diǎn)（x=15d和x=22d），獲取這兩個(gè)測點(diǎn)橫向速度時(shí)歷變化曲線，并進(jìn)行能譜分析，如圖7 所示。根據(jù)橫向速度時(shí)歷變化曲線，可見同一空間位置下的3 種網(wǎng)格規(guī)模算例所得到的速度振幅相近。在能譜分析曲線中，可以明顯觀察到3 種網(wǎng)格規(guī)模下的流場特征頻率相近。表3 給出的則是主頻及相應(yīng)幅值。由表可見，不同網(wǎng)格規(guī)模下計(jì)算得到的流場特征頻率幾乎完全一致，但隨著網(wǎng)格規(guī)模的擴(kuò)大，高頻區(qū)域的能量幅值逐漸增加，表明在較高網(wǎng)格分辨率下模擬能捕捉到更多小尺度運(yùn)動(dòng)。

表3 不同網(wǎng)格規(guī)模下方柱繞流數(shù)值模擬的流場特征頻率Table 3 The characteristic frequency of different grid sizes for flow around a square cylinder

圖7 不同網(wǎng)格規(guī)模下橫向速度時(shí)歷變化曲線及能譜分析（測點(diǎn)位置：z=2d，y=8d）Fig. 7 Time histories and energy spectrum analysis of transverse velocity of different grid sizes and where the detection point position is z=2d , y=8d

4 流動(dòng)特征分析

通過對比上述不同網(wǎng)格規(guī)模下的計(jì)算結(jié)果，表明較高網(wǎng)格分辨率更有利于對流動(dòng)現(xiàn)象的捕捉。綜合考慮計(jì)算耗時(shí)和硬件成本，本文對于方柱繞流流動(dòng)現(xiàn)象的研究采用了30.72 百萬網(wǎng)格規(guī)模進(jìn)行計(jì)算。其中，計(jì)算域尺寸和網(wǎng)格尺寸均滿足直接數(shù)值模擬的要求。

為研究流動(dòng)現(xiàn)象，本文將進(jìn)一步分析長徑比H/d=4 的方柱繞流流場，對比長徑比H/d=2，3，4的方柱繞流尾流渦結(jié)構(gòu)差異，以探究自由端剪切流對尾流渦結(jié)構(gòu)的影響。

4.1 尾流渦特征分析

圖8 所示為方柱繞流瞬時(shí)z方向的渦量（Ωz）等值線圖。由圖8（a）可以明顯觀察到方柱后方附近區(qū)域存在3 個(gè)截面上的渦量等值線。對比圖8（b）中的等值線，可明顯觀察到3 個(gè)截面上的等值線渦旋脫落（以下稱渦脫）幾乎相同。結(jié)合圖8（c）所示的方柱各高度截面渦脫模型[6]，也與本文計(jì)算的結(jié)果一致。

為進(jìn)一步量化對比不同高度截面上的渦脫特征，選取了30.72 百萬網(wǎng)格規(guī)模算例，對不同高度（z=0.5d，2.0d，3.5d）的水平截面，取方柱正后方x=15d和x=22d兩個(gè)位置的橫向速度（y方向速度）進(jìn)行檢測，并進(jìn)行能譜分析。由圖9 所示的橫向速度時(shí)歷變化曲線可見，無論是在近場區(qū)域（x=15d）或遠(yuǎn)場區(qū)域（x=22d），不同高度截面上的橫向速度具有近似一致的相位，與圖8（a）展示的關(guān)于方柱各截面同步渦脫的結(jié)論一致，且與文獻(xiàn)[9]中通過阻力計(jì)算所得結(jié)果相同。

圖8 方柱截面渦分布特征Fig. 8 The characteristics of vortex on square cylinder cross-section

對比橫向速度振幅，可見靠近底面或靠近自由端的截面上，橫向速度振幅相對于中部截面均有較大幅度的減小，表明自由端及底部的固體壁面對其附近的渦脫具有一定抑制作用。此外，對比這兩處的橫向速度振幅，遠(yuǎn)場處振幅均有一定程度的減小，表明黏性力起到了耗散作用。從圖9 所示的能譜分析結(jié)果可見，不同高度截面上的橫向速度具有相同的特征頻率，且中部截面上的橫向速度能量更高，與橫向速度時(shí)歷變化特征相對應(yīng)。

圖9 不同高度水平截面上橫向速度時(shí)歷變化曲線及能譜分析(測點(diǎn)位置：z=2d，y=8d)Fig. 9 Time histories and energy spectrum analysis of transverse velocity on horizontal sections at different heights and where the position of detection point is z=2d, y=8d

4.2 長徑比的影響

圖10 所示為長徑比H/d=2，3，4 時(shí)的瞬時(shí)（Q=0.01）渦系結(jié)構(gòu)圖。由圖可見，在H/d=3時(shí)，尾流仍然為反對稱的卡門渦；在H/d=2 時(shí)，尾流中的反對稱渦脫現(xiàn)象完全消失，表現(xiàn)為細(xì)長的流向渦二次結(jié)構(gòu)[10]。

圖10 不同長徑比方柱繞流瞬時(shí)渦系結(jié)構(gòu)（Q=0.01）Fig. 10 Instantaneous vortex structures of flow around a 3D finite square cylinder with different slenderness ratios (Q=0.01)

為了更加清晰地展示不同長徑比方柱繞流流場特征，給出了不同長徑比方柱1/2 高度截面（z/H=0.5）上的橫向速度時(shí)歷變化曲線及相應(yīng)的能譜分析結(jié)果，如圖11 所示。由圖可見，隨著方柱長徑比的減小，橫向速度振幅逐漸減小。在長徑比H/d=2 時(shí)，橫向速度時(shí)歷變化趨于一條直線。通過能譜分析可以明顯觀察到，隨著長徑比的減小，整個(gè)流場的交替渦脫特征逐漸減弱，并在H/d=2 時(shí)，上述渦脫特征幾乎完全消失。表4 給出了計(jì)算得到的不同長徑比方柱繞流流動(dòng)平均阻力系數(shù)（Cd）和斯特勞哈爾數(shù)（St）。本文計(jì)算結(jié)果與文獻(xiàn)[10] 采用的局部網(wǎng)格加密計(jì)算結(jié)果十分接近，驗(yàn)證了本文計(jì)算結(jié)果的正確性。

表4 不同長徑比方柱繞流流動(dòng)平均阻力系數(shù)及斯特勞哈爾數(shù)Table 4 The average flow resistance coefficient and Strouhal number of the flow around a 3D finite square cylinder with different slenderness ratio

圖11 不同長徑比方柱繞流橫向速度時(shí)歷變化曲線及能譜分析（截面高度：z/H=0.5）Fig. 11 Time histories and energy spectrum analysis of transverse velocity of flow around a 3D square cylinder with different slenderness ratio at the cross-section height of z/H=0.5

由于來流速度是恒定的，因此采用時(shí)間平均方法對湍流運(yùn)動(dòng)進(jìn)行統(tǒng)計(jì)，可更加直觀地反映出不同長徑比方柱繞流流動(dòng)的特征。圖12 所示為不同長徑比方柱繞流中縱截面（y=8d）上的流線分布。由圖可見，在長徑比H/d=3，4 的方柱繞流中，方柱后方均存在上、下兩部分回流區(qū)，其中上部回流區(qū)由自由端剪切層分離的下洗流形成，下（根）部回流區(qū)回流方向與上部回流區(qū)回流方向相反，且在流場下游均形成了流線分離點(diǎn)（鞍點(diǎn)）“+”，這是因?yàn)樯?、下回流區(qū)在橫截面上的渦脫方向相反所致。相較而言，隨著方柱長徑比由4 減小到3 時(shí)，上部回流區(qū)的尺寸顯著減小，表明由自由端剪切層分離的下洗流強(qiáng)度減弱。與長徑比為3，4 的方柱繞流現(xiàn)象不同的是，當(dāng)長徑比持續(xù)減小至2 時(shí)，尾流中沒有明顯的鞍點(diǎn)存在，尾流完全由自由端剪切層分離的下洗流控制，流場特征與臨界長徑比值均與文獻(xiàn)[5]中所述結(jié)論一致。

圖12 不同長徑比方柱繞流時(shí)間平均流場中縱截面流線分布（截面位置：y=8d）Fig. 12 Streamline distribution of longitudinal section (y=8d) in time average flow field with different slenderness ratios

5 結(jié) 論

本文通過基于“神威·太湖之光”超級計(jì)算機(jī)對雷諾數(shù)Re=250 的三維有限長方柱繞流進(jìn)行了直接數(shù)值模擬，網(wǎng)格規(guī)模為3.84 百萬～245.76 百萬，通過對長徑比H/d=4 的方柱繞流流動(dòng)進(jìn)行數(shù)值模擬，對網(wǎng)格規(guī)模效應(yīng)及其對應(yīng)的流場特征進(jìn)行了總結(jié)與歸納。在上述基礎(chǔ)上，對H/d=2，3 的有限長方柱繞流進(jìn)行了計(jì)算及對比分析，得到以下結(jié)論：

1） 基于“神威·太湖之光”超級計(jì)算機(jī)的多級并行計(jì)算能夠大幅度提高計(jì)算效率，顯著提升大規(guī)模并行計(jì)算能力。

2） 對比不同網(wǎng)格規(guī)模的計(jì)算結(jié)果表明，高分辨率網(wǎng)格有助于捕捉更多的流場信息。對于三維有限長方柱繞流，不同高度截面上的渦脫是同步的。而不同高度截面上的橫向速度幅值及能譜分析結(jié)果表明，自由端剪切層分離的下洗流和底面的黏性力均會(huì)抑制附近截面上的渦脫。

3） 對比不同長徑比方柱繞流現(xiàn)象的結(jié)果表明，當(dāng)長徑比為2 時(shí)，自由端剪切層分離的下洗流覆蓋了方柱后壁面附近區(qū)域，抑制了整個(gè)方柱的渦脫，導(dǎo)致尾流中的反對稱渦脫現(xiàn)象消失。

綜上所述，基于“神威·太湖之光”超級計(jì)算機(jī)的多級并行計(jì)算在船舶水動(dòng)力學(xué)領(lǐng)域有較好的應(yīng)用潛力，但在實(shí)際工程研究中采用的網(wǎng)格及幾何模型相較于本文更加復(fù)雜，若應(yīng)用于實(shí)際工程領(lǐng)域，仍需要開展深入研究。