随着国家对E级计算机的开发，高性能软件的重要性日趋显现，计算科学发展显著影响智能制造等国家战略实施。计算机建模与仿真已经被一些国家列为服务于国家利益的需要优先发展的核心技术。有学者指出，自主CAE软件的开发应被上升到国家战略发展的高度，国内外已有大量的专用和通用的CAE软件开发。颗粒流动是自然界中普遍存在的一种流动现象，如工业锅炉煤粉燃烧、内燃机燃油喷雾燃烧、大气环境中的悬浮颗粒污染物扩散等，都存在离散颗粒在连续介质中的运动现象。由于颗粒两相流在工业过程的广泛存在，高性能的颗粒两相流并行计算方法也是现代CAE软件的重要内容。OpenFOAM创始人Weller^[1]和思克莱德大学的Macpherson在OpenFoam的基础上，基于单元面的面心坐标及其法向量，开发了计算效率高且健壮性好的可用于任意多面体网格的并行颗粒追踪模块。Onishi等^[2]对均匀各向同性湍流中惯性粒子的颗粒追踪进行了并行化，其研究一个规则三维区域，通过MPI传递各个进程间的流动与颗粒影像信息。Buaria等^[3]在二维笛卡尔网格上应用一种改进的颗粒分配方法，基于PGAS通信模型对大尺度湍流中的颗粒追踪进行研究，该方法着眼于降低计算中的通信量，使得颗粒与流体单元的耦合只在本进程内发生，从而大大减少了通信时间和内存占用。自2003年起，哈尔滨工程大学启动CFD软件开发，于2007年取得软件著作权，并成功应用于传热与燃烧过程分析、燃烧室内的流动，船舶水动力学分析，柴油机缸内工作过程模拟等。2008年，进行计算声学和结构应力波模拟软件开发，并对结构振动与噪声问题进行数值预报。

1 GTEA并行模块开发

GTEA基于连续介质假设，对控制方程按时间、对流和扩散项以及源项进行离散，在此思想下，统一流体力学、结构动力学、声波动控制方程求解^[4]。软件采用非结构化网格离散控制方程^[5]，建立了通用的软件开发模式。本文在哈尔滨工程大学自主研发的GTEA通用数值计算平台中内燃机缸内喷雾燃烧模块，使用区域分解方法，开发适用于各种物理模型的通用并行计算模块。

1.1 GTEA的并行设计模式

在科学与工程计算中，有多种并行设计模式可选择，最常见的设计模式包括Fork-Join模式、区域分解模式等，如图 1所示。Fork-Join模式是最传统的并行计算方法，寻找程序的性能瓶颈，即代码中的“热点”^[6]是其主要步骤。针对“热点”设计并行算法，以减少计算资源的消耗。

针对独立处理的动网格、喷雾、流动、燃烧等模块，可采用Fork-Join模式^[7]。这种设计方法具有简单便捷的特点，无需考虑求解器的整体结构，只需使热点部分得到充分的并行化即可。热点部分的代码量通常较少，配合以编译制导语句为主的编程“语言”如OpenMP，可极大简化并行编程工作。但是由于这种设计模式需要首先识别程序中的“热点”，而“热点”依赖具体的模型复杂度，对通用CFD软件缺乏良好的自适应能力。

气相流动过程GTEA采用了欧拉方法，因此在程序开始求解之前，基于区域分解模式将网格根据并行进程数剖分，每个进程执行相同的求解器程序求解子网格的流动方程^[8]，并行程序需要在原串行程序基础上增加的程序模块如图 2所示。

1.2 GTEA的并行软件模块

GTEA采用了常用的MPI编程环境^[9]来实现并行，考虑到程序的可扩展性，对并行相关的数据结构和子程序进行封装和模块化，并形成独立的并行计算模块。

GTEA并行计算模块由并行公共数据结构，并行通用函数，前后处理的外部接口，面向求解器各模块的并行支持子模块等4个部分组成。

2 两相喷雾燃烧并行计算方法 2.1 欧拉相并行计算方法

流场计算是缸内喷雾燃烧过程模拟的基础，在欧拉相描述下，流场表现为求解区域中的物理量空间分布。GTEA采用非结构有限体积法求解连续介质物理场，物理量分布以多面体网格单元为基本单位，因此欧拉相并行计算需要对网格进行二次剖分。

GTEA实现了多种二次剖分策略，包括simple、metis和几何法。其中simple方法为简单平均方法，从非结构网格的单元序号出发，累加各单元的权重，通过目标子网格块数对总权重进行平均划分。几何法使得计算用户可以在二次剖分中考虑一定的几何特征，该方法从单元的形心坐标出发，考虑各单元所处的空间位置，根据剖分策略的几何特征来选择相应的切割面走向，比如沿笛卡尔坐标轴剖分的Cartesian X/Y/Z Coordinate策略，沿中心轴周向划分的Cylindrical Axes策略等。metis策略是GTEA默认的二次剖分方法，该方法基于多层次k路图划分算法(MultiLevel k-way Partitioning)，可以同时对负载平衡和通信量进行控制，首先将非结构化网格抽象为无向赋权图G(V, E)，其中V代表了每个单元，而E代表了单元与单元之间的通信关系，这样网格二次剖分问题就转化成了图剖分问题。同时这些节点和边还可以有各自的权值来代表每个单元的计算量和通信量。多层次k路图划分算法将图划分过程分为3个子过程：图的粗化、初始划分、图的细化，其过程如图 4所示。

计算网格经过二次剖分后形成较小规模的子网格，每个进程仅在其映射的子网格上进行求解任务，为了减少内存占用，在各进程上对子网格重新排序，形成本地单元索引。在计算任务完成之后，可以根据单元映射关系将各进程的计算结果进行合并，重建全局计算结果，从而保持和串行计算相一致的结果输出形式。

由于GTEA采用FVM求解，每个单元通过单元面与相邻单元产生数据关联。子网格的切割面两侧单元分别属于不同的进程，而MPI进程的内存空间又相互独立，因此必须设计相应的信息传递过程交换进程间的单元数据。GTEA的通信过程为相邻进程两两交换切割面两侧必要的数据，这种通信方式避免了全局通信操作收集全进程单元巨大的通信量，保证了仅有所有必要的数据得到了传递，因而是稀疏的。为了保证稀疏通信的进行，GTEA设计了进程间通信映射关系的SYSTEM对象和传递数据容器的ghost单元，并在梯度更新、方程组求解、压力修正等过程中调用通信程序。

2.2 拉格朗日相并行计算方法

GTEA喷雾基于离散液滴模型(discrete droplet model, DDM)，液相在拉格朗日体系下进行描述，每个样本都代表一定数目的液滴，这些具有相同大小和状态的液滴被看成一个颗粒，不考虑全部液滴，只处理其中若干具有代表性的统计样本。离散液滴在连续介质中运动的过程中，通过子模型模拟破碎、碰撞、蒸发、碰壁等一系列物理过程，从而与其相邻颗粒、容器单元产生数据依赖。

GTEA使用链表数据结构来描述离散颗粒，其长度可增减，且在内存中的分布不连续。有2种策略可以将颗粒映射给不同的MPI进程，如图 6所示，一种方法是根据链表长度将每一时刻新增的颗粒平均分配到不同的MPI进程；另一种方法是根据颗粒所属的空间位置，将其分配给对应容器单元所属的MPI进程，当颗粒运动路径穿过切割面时，将其移动到另一侧MPI进程，并还原其运动轨迹。

第1种分解方案是完全静态的，颗粒一旦产生就不会在进程之间产生传递，但是这种方案有2个缺陷：1)颗粒的相邻颗粒可能不在其同一进程上；2)颗粒所属的容器单元可能不在其同一进程上，为了保证碰撞、蒸发、动量交换等过程的计算，必须设计一个主进程对全单元全颗粒的物理信息进行收集，这样在欧拉相稀疏通信中所避免的全局通信在拉格朗日相中被再现了。颗粒计算在喷雾燃烧中所占的比例较小，为了保持其负载平衡导致全局通信量的增加是得不偿失的，因此GTEA采用了第2种分解策略，当颗粒新增进入计算域后，通过FFPT快速追踪算法定位其运动路径和容器单元，并映射给相应的MPI进程。

颗粒链表由于不具有索引，且长度不断变化，运动范围较广，其通信过程比欧拉单元更不易组织，为此通过建立通信缓冲区来实现颗粒的进程通信。由于欧拉相中建立了用于通信的ghost单元，这就为颗粒通信提供了判据，当FFPT算法计算颗粒进入ghost单元后，即可认为此颗粒已穿出进程，此时立刻将此颗粒移入缓冲区并从主链表中删除，同时标记该进程未完成颗粒的数量。此后，通过对缓冲区中的颗粒数据建立索引，动态更新映射，组织通信过程，还原运动路径，即可实现颗粒在进程间的传递。这种通信方式仅传递跨过切割面的颗粒，避免了全局通信，从而大大降低了通信量。

3 仿真算例

为了验证GTEA开发的并行计算模块的功能，搭建了HPC服务集群服务器(如图 7)。硬件上，采用分布式架构，机架式管理节点和高密度的刀片式计算节点混合组网，通过InfiniBand相连，同时，通过SAN网络与EMC存储节点相连。软件上，部署了RHEL6、OFED，ICC/IFORT、MVAPICH、PBS torque、xCAT、Ganglia等软件保障服务器的运行和管理。

并行计算准确性与加速比验证算例基于Sandia国家实验室的定容燃烧弹试验数据^[10]，该定容燃烧弹的计算网格由HyperMesh网格生成器完成，最后转化为CGNS标准格式，读入GTEA软件，网格如图 8所示。所生成的网格文件共有11 250单元，32 625个内部面，边界面1 500。湍流模型、雾化模型、碰撞模型、燃烧模型风别采用标准k-ε模型、KH-RT模型、O′Rourke模型和Shell-CTC模型^[11]。

定容燃烧弹内初始给定混合气体成分，其中22%氧气，77.8%氮气，0.1%二氧化碳，0.1%水蒸气。燃料从模型底部中心喷入，主要成分为C14H30，由预先充满的高温高压气体引燃。

计算中对时间步的喷雾模拟的结果都输出，因此详细地捕捉燃油喷雾的发展趋势，提取相关的计算结果进行对比验证。如图 9所示，为进程数分别为1、2、4、8相应时刻的液相雾化形态对比，可见其趋势基本一致。

定容燃烧弹内的温度和索特平均直径随时间变化见图 10，燃油约在0.25 ms被引燃，温度迅速升高，同时喷雾索特平均直径减小，温度与液滴直径变化相吻合，不同进程数计算结果吻合较好。计算进行到1 ms左右，定容弹内过程达到相对平衡，温度变化趋于稳定，同时索特平均直径平均值附近振荡变化。不同进程数得到的结果基本一致，证明了并行计算结果的可靠性。

定容燃烧弹内温度分布随时间的变化如图 11所示，温度分布可以表示火焰结构，从图中可以看出计算得到的火焰结构在相同的时刻是基本一致，与进程数无关。火焰浮起长度、火焰内的温度分布也类似，结果与进程数无关。结果表明，并行计算保持了着火与火焰扩散过程的计算准确性。

图 12给出了加速比和效率对比(具体定义见文献[6])，理想的线性加速比如图中的虚线所示。增加计算的进程数，相同网格数所需的总计算时间下降，与串行计算相比，并行计算大大缩短了计算时间。但是，受到了Amdahl′s定律的制约，加速比与理想加速比有一定偏差。

4 结论

1) 对于欧拉相，将网格二次剖分后映射给不同的进程，各进程在子网格上求解，并设计了稀疏通信保证各进程传递必要的数据。

2) 对于拉格朗日相，根据颗粒的空间位置决定其所在的进程，基于FFPT算法决定是否穿出进程，并通过缓冲区来组织进程间的颗粒通信。

3) 经过定容弹喷雾燃烧模拟算例，通过对比分析喷雾形态，温度变化，索特平均直径，火焰发展等计算结果，可知在不同的进程下计算结果都与串行结果相一致。

4) 并行计算的加速比以近似线性增长，有效地减少了计算时间，加速效率不低于60%。