基于“神威·太湖之光”的三维有限长方柱绕流直接数值模拟

张亚英，吴乘胜，王建春，金奕星

中国船舶科学研究中心，江苏 无锡 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 时，自由端剪切层分离的下洗流覆盖了方柱后壁面附近区域，抑制了整个方柱的涡脱，导致尾流中的反对称涡脱现象消失。

综上所述，基于“神威·太湖之光”超级计算机的多级并行计算在船舶水动力学领域有较好的应用潜力，但在实际工程研究中采用的网格及几何模型相较于本文更加复杂，若应用于实际工程领域，仍需要开展深入研究。