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

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

中國船舶科學研究中心，江蘇 無錫 214082

0 引 言

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

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

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

1 并行方案

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

圖1 SW26010 處理器結構Fig. 1 The structure of SW26010 processor

圖2 并行計算任務分配示意圖Fig. 2 Schematics of data allocation in parallel computing

為檢驗上述并行方法的可行性，采用SIMPLE算法及直接數值模擬方法的并行效果，本文以三維頂板驅動方腔流為模型，將雷諾數設置為Re=1 000，對不同網格規模及不同進程數下的并行加速效果予以對比。

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

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

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

根據表1 所示結果，MPI 并行效率在128 進程時均高于50%；而MPI+Athread 基于核組的擴展效率偏低。結合文獻[4]對SW26010 處理器的從核并行加速影響因素的研究，經分析發現導致上述現象的原因包括：

1) 從核計算的通信過程使整個通信耗時增加，降低了并行效率。

2) 隨著進程數的增加，單個進程計算量減少。根據文獻[4] 所述結論，此時將造成計算耗時比例下降，通信耗時比例提升，導致從核并行計算加速比明顯降低。

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

2 計算模型與控制方程

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

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

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

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

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

采用基于交錯網格的有限體積法（finite volume method，FVM）離散控制方程。時間項離散采用一階顯式格式，對流項采用二階迎風插值格式（QUICK），擴散項采用中心差分格式，具體可見文獻[7]。控制方程組的求解采用SIMPLE算法。

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

3 網格規模效應分析

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

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

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

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

使用Q判據[8]顯示方柱周圍及尾流中的瞬時渦系結構，其形式如下：

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

圖6 t=450 s 時刻流場瞬時渦系結構（Q=0.01）Fig. 6 Instantaneous vortex structures at t=450 s when Q=0.01

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

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

圖7 不同網格規模下橫向速度時歷變化曲線及能譜分析（測點位置：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 流動特征分析

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

為研究流動現象，本文將進一步分析長徑比H/d=4 的方柱繞流流場，對比長徑比H/d=2，3，4的方柱繞流尾流渦結構差異，以探究自由端剪切流對尾流渦結構的影響。

4.1 尾流渦特征分析

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

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

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

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

圖9 不同高度水平截面上橫向速度時歷變化曲線及能譜分析(測點位置：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 時的瞬時（Q=0.01）渦系結構圖。由圖可見，在H/d=3時，尾流仍然為反對稱的卡門渦；在H/d=2 時，尾流中的反對稱渦脫現象完全消失，表現為細長的流向渦二次結構[10]。

圖10 不同長徑比方柱繞流瞬時渦系結構（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）上的橫向速度時歷變化曲線及相應的能譜分析結果，如圖11 所示。由圖可見，隨著方柱長徑比的減小，橫向速度振幅逐漸減小。在長徑比H/d=2 時，橫向速度時歷變化趨于一條直線。通過能譜分析可以明顯觀察到，隨著長徑比的減小，整個流場的交替渦脫特征逐漸減弱，并在H/d=2 時，上述渦脫特征幾乎完全消失。表4 給出了計算得到的不同長徑比方柱繞流流動平均阻力系數（Cd）和斯特勞哈爾數（St）。本文計算結果與文獻[10] 采用的局部網格加密計算結果十分接近，驗證了本文計算結果的正確性。

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

圖11 不同長徑比方柱繞流橫向速度時歷變化曲線及能譜分析（截面高度：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

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

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

5 結 論

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

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

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

3） 對比不同長徑比方柱繞流現象的結果表明，當長徑比為2 時，自由端剪切層分離的下洗流覆蓋了方柱后壁面附近區域，抑制了整個方柱的渦脫，導致尾流中的反對稱渦脫現象消失。

綜上所述，基于“神威·太湖之光”超級計算機的多級并行計算在船舶水動力學領域有較好的應用潛力，但在實際工程研究中采用的網格及幾何模型相較于本文更加復雜，若應用于實際工程領域，仍需要開展深入研究。