YH-ACT：热工流体力学并行应用程序*

刘 杰，龚春叶，杨 博，郭晓威，甘新标，李胜国， 李 超，陈旭光，肖调杰，穆利安，宋 敏，赵冬勇，鞠羽中

(1.国防科技大学计算机学院，湖南 长沙 410073；2.复杂系统软件工程湖南省重点实验室，湖南 长沙 410073)

1 引言

核动力装置的安全研究是发展核电的基础，利用数值模拟技术进行核动力装置的研究是当前核科学与工程中的主要手段[1]。早期的核反应堆热工水力的数值研究大多基于大尺度的几种参数或一维系统分析方法，主要研究核动力系统中物理量的平均特性。近年来，随着计算机的快速发展，基于更为精细化的计算流体力学(CFD)方法和热工流体力学软件的核动力热工水力研究得到广泛研究和应用[2]。

目前我国核反应堆热工CFD软件仍主要采用国外商业软件。由于西方国家对CFD软件核心技术的禁运，导致无法根据应用需求进行物理模型的修改和扩展，也无法高效利用国产高性能计算机提供的计算资源。商业CFD程序不能完全满足我国反应堆应用的需求。国外CFD软件起步早，20世纪80年代英国CHAM公司推出第1个具备通用性的CFD商业软件PHOENICS[3]。此后，STAR-CD、Fluent和CFX也相继问世于英国并成为主流CFD分析软件。2003年和2006年CFX和Fluent先后被ANSYS收购集成为一整套CFD解决方案[4]。作为后起之秀的NUMECA[5]，由于其高速流动计算性能强，在跨、超音速区的高速压缩机、汽轮机等模拟中具有明显优势。法国的STAR-CCM+[6]堪称为“下一代CFD解决方案”，主要支持多面体网格，内存开销少且求解速度快，但其在核反应堆领域中的应用尚处初期阶段。此外，Exa公司研发了一款基于格子玻尔兹曼方法LBM(Lattice Boltzmann Method)[7]的CFD商业软件POWERFLOW[8]，但仍有许多技术问题尚待解决。

目前，CFD商业软件已广泛用于核反应堆的热工水力模拟。2006年重水反应堆设计时就较早使用CFX 10.0模拟燃料冷却剂通道问题[9]。同年，德国安全研究所利用CFX模拟下腔室的混合流动[10]。CFX也在反应堆安全领域被广泛用于与REACFLOW等程序的耦合分析[11 - 13]。意大利比萨大学利用Fluent程序模拟超临界水堆中的不稳定现象[14]。德国反应堆安全研究所等使用CFX进行了多尺度CFD模型模拟[15]。上述工作均关注反应堆中的单相流问题，而对于多相流问题仍处于系统评价和基准校验阶段，代表性的工作是美国能源部组织的轻水堆先进仿真联盟(CASL)[16]和核能先进仿真建模(NEAMS)[17]2个项目，代表着目前核能仿真及数值反应堆领域的最新进展。

不同于商业软件，CFD开源软件源代码公开且可免费获取，最具代表性的是发源于英国帝国理工大学的OpenFOAM[18 - 20]。其基于有限体积法，能够求解湍流、热传递、不可压及可压缩流体等流动问题，可以针对特定应用问题方便地编写自定义的求解器。此外，由法国电力集团于1997年开始研发的通用CFD开源软件Code_Saturne[21]，能够完成流体动力学中的湍流、传热等流体计算，已用于核能源相关研究。国防科技大学计算机学院与中国核动力院合作研制的基于OpenFOAM的CFD模块能够预测卸压箱高速喷气过程等简化的工程应用场景[22]。但总体来说，这类开源软件的计算精度和计算效率等无法完全满足核反应堆工程设计的要求。

相比而言，国内面向核反应堆的CFD技术及其软件研发起步较晚。中国工程物理研究院高性能数值模拟软件中心自主研制了支持结构网格的JASMIN框架[23]、支持非结构网格的JAUMIN框架[24]和支持实体组合几何的JCOGIN框架[25]。中国空气动力研究与发展中心研发的风雷软件平台[26]和混合网格CFD软件HyperFlow[27]，主要面向气动力学仿真。中国科学院数学与系统科学研究院开发的系统PFEPG (Parallel Finite Element Program Generator)[28]，采用组件化程序设计方法和软件自动生成技术，可自动生成基于区域分解算法的并行有限元程序。尽管国内在CFD程序研发方面已取得一定成果，但与国外相比，在物理模型全面性、计算精度、计算效率以及软件易用性等方面仍存在差距，不能完全满足核反应堆热工水力应用开发的要求。

目前，国内核反应堆热工水力CFD分析仍主要采用国外商业软件。上海核工程研究设计院利用CFX模拟了反应堆安全阀结构及蒸汽发生器二次侧局部流场[29]；西安交通大学使用Fluent模拟应急堆芯冷却系统冷水注入压力容器时的瞬态流动特性[30]；中国核动力研究设计院利用CFX分析棒束定位格架的单相流场特性[31]；深圳中广核工程设计有限公司采用CFX模拟反应堆进行了压缩空气隔离阀开启瞬间流速计算[32]；为真实模拟复杂热工水力现象，清华大学工程物理系先后通过有效耦合COBEA-IV与CFX以及RELAP5与CFX程序实现了多尺度模拟[33,34]；海军工程大学船舶与动力工程学院实现了时空中子动力学模型与CFX联合模拟反应堆局部三维流场[35]。还有许多研究工作开始关注不同类型反应堆的热工水力特性，如借助CFX模拟第4代核能系统中多孔介质传热现象等[36]。

综上所述，商业CFD程序已广泛应用于反应堆的热工水力模拟，但不能完全满足反应堆的应用需求；开源CFD程序有部分应用，但与商业CFD程序相比，在物理模型全面性、计算精度、计算效率和易用性等方面仍存在差距。因此，为更好地满足局部精细热工水力分析的需求，需要更全面的物理模型、较高的计算精度和较好的并行计算效率，有必要开发自主热工CFD软件。本文详细描述了热工流体力学并行应用YH-ACT(Parallel Analysis Code of Thermohy-draulics)的设计、实现方案和测试结果。

2 热工CFD物理模型

在核反应堆热工水力分析过程中，需要对冷却剂的流动传热过程进行建模。完整描述流体中的流动以及传热过程需要的理论模型包含3部分：质量守恒方程、动量守恒方程和能量守恒方程[37]。当温度较高时，还需考虑热辐射，需要对热量的辐射过程进行建模。反应堆中的流体(冷却剂)区域，堆芯燃料组件及管壁等固体区域的热量传递过程对整个运行安全和性能都至关重要，因此，还需要研究流-固耦合传热模型。

压力容器内结构复杂，流动尺度范围较大，且工程实际应用中更关注物理量的平均量，因此，一般对上述基本三大方程进行雷诺平均，计算流体流动的平均物理量。平均化后的动量方程中出现雷诺应力项，需要增加湍流模型才能使得平均化后的质量、动量、能量方程组封闭，从而进行堆芯冷却剂流动和传热过程的模拟。

YH-ACT目前主要考虑雷诺时均湍流模型的实现，当温度较高时，还需考虑热辐射，求解辐射输运方程的目的是获得能量方程中的能量源项，并获得与固体壁面辐射热通量相关的物理量。至于流-固耦合传热模型，一般通过流体域和固体域的交界面建立流体-固体域间的热量传递模型，计算固体域内能量方程时，将流体域交界面上的相关物理量作为边界条件；而计算流体域内能量方程时，将固体交界面上的热量(温度)分布作为流体能量方程的源项，最终实现流体-固体的耦合传热。

3 热工CFD方程数值求解方法

假定流体为稳态不可压缩，且忽略重力影响，动量守恒方程和能量守恒方程可以表示成以下通用形式：

(1)

要使封闭方程得到定解，还需要满足一些边界条件，YH-ACT程序中包含的主要边界类型有进口边界：速度进口、压力进口等；出口边界：outflow、压力出口等；壁面边界：滑移壁面、无滑移壁面、旋转壁面；对称边界；周期边界等。

对任意形状的控制体积，引入高斯散度定理，得积分形式如下：

(2)

其中，V为控制体体积；S为控制容积的表面积；S为控制容积界面的面积矢量，其方向与外法线单位矢量一致。

假设对流项采用一阶迎风格式，Laplacian项采用非正交处理，源项进行线性化处理，经过一系列分析，最后可得离散方程。对离散后的控制方程组的求解可采用耦合式解法及分离式解法。

4 并行求解算法

4.1 无冗余网格数据分解并行求解算法

一种简单直接的并行方法就是对网格数据进行分解。将网格数据划分成多份后分配给不同进程计算，过程中通过进程间通信来实现数据交换。

网格数据分解如图1所示。求解区域被划分为互不重叠的网格单元，每个进程仅读取分配给本进程的网格数据。这样不可避免地需要进行进程之间的数据通信。无冗余的网格数据分解方法中各个进程之间并没有共享网格数据。由于采用的离散格式最高不超过2阶，计算过程中仅需要使用相邻1～2个网格单元的值，这样在并行数据划分之后，实际仅需要获取与边界相邻的网格单元中的数据。因此，可以将需要获取的跨进程数据存储在边界上(如图1中虚线所示)，此时，求解可分解为如下3步：边界初始化(获取相邻处理器上的数据)、本地计算和边界更新(更新相邻处理器上的数据)。

Figure 1 Decomposition of non-redundant grid data图1 无冗余网格数据分解

4.2 通信模式

基于循环遍历所有网格块的方式预先建立通信关系，再根据通信关系实施通信。对于每个网格块来说，它只需要与其邻居网格块进行通信，建立通信关系就是找到所有的邻居网格块，并存放到一个本地数组中。对于任意一个网格块来说，它可见的网格块就只有它本身及其邻居这些“局部化”的网格块。当并行规模提高时，网格块在通信时不必遍历所有的网格块，而只需遍历“局部化”的网格块即可，大大提高了通信的效率。采用MPI异步通信可进一步提高通信效率。一个通信模式示例如图2所示。

Figure 2 Illustration of communication mode图2 通信模式示例

通信伪代码如下所示：

//通信过程的伪代码

voidCommunicate(/*args*/){

for(intizone= 0;izone

intsrcZoneID=GetZoneID(izone);

if(srcZoneID == currentZoneID){

for(intidstZone= 0;idstZone

intdstZoneID=GetZoneID(idstZone);

if(dstZoneID==currentZoneID) continue;

MPI_Send(dstZoneID,/*other args*/);

}

}

else{

MPI_Recv(srcZoneID,/*other args*/);

}

}

}

4.3 多级并行架构

基于上述网格数据分解，可以进一步进行两级或多级分区。采用MPI/OpenMP混合并行模式。在一级分区上用 MPI并行，然后将一级分区划分为若干个二级子分区，二级子分区的数量一般等于计算结点内部的核数，在二级分区上进行粗粒度的 OpenMP 并行。

这种多层级的并行架构适应于目前的超级计算机集群的架构，一级分区分配在各个计算结点上，而二级分区在同一个计算结点之上。其同样适用于异构体系，第1层仍然使用 MPI 通信，即 CPU 核间或计算结点间、基于网格一级分区交界面信息的MPI 通信；第2层是 OpenMP 并行，即计算核心间、基于网格二级子分区的粗粒度并行，CPU和协处理器间采用相应的通道进行数据传输。当没有协处理器时，上述异构并行模式自动退化为同构并行模式。为了更好地控制不同的分区层级，将整个计算域分为2层：第1层为Region分区，第2层为Zone子分区。在并行计算时，Region分区被分配到计算结点，对应一个MPI进程，Zone子分区被分配到计算核心线程，对应一个OpenMP线程。

典型的CFD求解器执行流程主要包括残差初始化、时间步长计算、 通量计算、时间推进、流场更新、交界面通信等过程。求解器Solver被应用到计算域Zone上，针对不同的物理模型，一个Zone上可能还会加载多个Solver。在计算过程中，Region依次遍历其中的Zone，调用相应的求解器执行计算任务；当使用OpenMP并行时，则可以对每个Zone并发执行计算任务。

4.4 多级并行通信优化

在多级并行架构中，每一个MPI进程上可能有多个子网格块，如果直接按照前文所述的循环遍历网格块的模式进行MPI/OpenMP混合并行通信，会存在进程间多次通信的问题。以图3为例，Region 0和Region 1有2组相邻的子网格块，每次交换数据都需要2个收发过程。因此，可以将进程间的网格分区数据进行打包，如图4所示，以减少通信的次数。

Figure 3 Illustration of partition hierarchy图3 分区层级示例

Figure 4 Example of process data packaging图4 进程数据打包示例

5 并行应用程序软件架构

5.1 YH-ACT程序架构

YH-ACT程序采取基于组件的架构，针对核反应堆热工水力CFD程序需求，以CFD求解中的数据加载、物理前处理、求解、结果输出等关键步骤构成的CFD程序流程为框架蓝图。以共性并行算法和并行支撑为基础，结合现有先进的软件工程方法，形成了基于国产超算系统，以物理模型模块、数值离散模块、线性方程组模块、并行通信模块等为主要部件的CFD程序总体框架，如图5所示。

Figure 5 Framework of YH-ACT program图5 YH-ACT程序框架

前处理工具箱具有网格控制、物理前处理等功能；物理模型和数值求解构成CFD求解器，封装高保真度物理数学模型和高效数学求解方法，并具备参数化配置能力。共性算法实现稀疏迭代法与预条件子等并行计算，通过与计算机体系结构的紧密耦合，实现大规模并行计算优化的支持，为领域解法器提供基础性计算能力；并行支撑可以封装高性能计算机系统结构特点，提供基础性的计算、通信、存储、调试等功能，从而支持大规模CFD计算应用高效运行。

基于上述分析，本文将YH-ACT程序垂直划分为4个逻辑层，即用户交互层、模型算法层、并行计算层和数据层，如图6所示。其中用户交互层负责接收用户操作，展现用户操作接口，这一层要求实现时为用户提供友好的界面，符合科研人员的操作习惯，提高用户体验。而用户操作的具体业务逻辑则由模型算法层实现，这样将业务逻辑和用户操作分开的一个好处是当业务逻辑改变后可以仅改变模型算法层而不影响用户操作层，且也不影响用户原来的体验。同样，模型算法层仅涉及业务逻辑，关于实体数据的一切操作都封装进并行计算层，数据层负责数据的持久化并向并行计算层提供数据访问接口。

Figure 6 Logical framework diagram图6 逻辑框架图

YH-ACT程序采用了典型的MVC(Modal View Controler)设计模式，具有耦合性低、重用性高、生命周期成本低、部署快、可维护性高等优点。其中M是指数据模型，V是指用户界面，C则是指控制器，在YH-ACT程序的架构中业务逻辑层和数据访问层担当的就是控制器的角色。使用MVC的目的是将M和V的实现代码分离，从而使同一个程序可以使用不同的表现形式。比如一批统计数据可以分别用柱状图、饼状图来表示。C存在的目的则是确保M和V的同步，一旦M改变，V应该同步更新。

5.2 GUI界面

Figure 7 GUI interface of YH-ACT program图7 YH-ACT程序GUI界面

YH-ACT程序的GUI界面如图7所示，界面大致分为菜单工具栏、控制树、参数设置页、视图浏览窗口、信息窗口5个区域。菜单工具栏有组织地展示前处理过程的全部功能。控制树管理全部模型内容，包括网格、求解域、边界条件、材料数据库、化学反应等。视图浏览窗口以图形的方式直观呈现模型计算结果。信息窗口用于反馈CFD命令的执行情况。因常用的应用程序GUI中菜单栏与工具栏同时存在且功能上没有区别，所以界面此处将两者结合为菜单工具栏，是CFD前处理应用程序里所有功能的大集成，且以文字与图形按钮结合的形式展示，最常用的功能以图形按钮体现。其摒弃了传统逐层嵌套的多级下拉菜单组件，将功能项按类划分到不同的选项卡分页中，每个分页内部又继续包含多个功能组，各种相关功能命令视为一组。界面的样式采用目前桌面软件主流的Ribbon纽带扁平风格，以简洁清晰的方式凸显CFD前处理应用功能主题。

GUI设置界面主要包含2部分：左侧控制树显示网格和数值计算所需的相关设置目录，以树形结构布置；中间参数设置页为具体的参数详情属性值编辑而设置；而右侧视图浏览窗口以图形界面显示相关的几何和网格信息，并随用户的设置显示相关的边界信息。首先，导入计算的网格，导入成功后会在Mesh菜单下显示导入的网格信息，包含outdown、mashup、kong21和core2一共4个部件的网格。然后，在Materials菜单下创建流体介质的相关物理性质。接着，在Simulation菜单下的Analysis Type下设置为瞬态计算模型，并设置瞬态计算时间等。双击Default Domain菜单，设置计算域中的流动介质和计算模型，并在Default Domain菜单创建边界条件(Inlet和Outlet等)，与此同时对瞬态计算初始化。最后，在Solver→Solver Control菜单下设置数值求解的离散格式和收敛准则，并在Solver→Output Control菜单下设置求解结果的输出。通过以上步骤，完成tran.mesh算例的设置，并将其保存为tran.sim文件后，通过求解器就能利用tran.sim进行算例的求解。

6 验证与测试

6.1 功能验证

本节选取3个典型案例模型进行功能验证，测试平台系统配置为：单结点2个12核CPU，共24核。主频2.3 GHz，系统内存为125 GB。

(1)案例模型1：均匀热通量的管道层流流动，如图8所示。

计算域为半径R= 0.0025 m，长度L=0.1 m的圆管道。Vx表示x方向的速度，τ表示剪切应力，PI表示进口压力，PO表示出口压力。网格划分为非结构网格，网格数量总计为5 080个。网格剖分如图9所示。

Figure 8 Fluid computing domain of the case model 1图8 案例模型1的流体计算域

Figure 9 Meshing of the case model 1图9 案例模型1的网格划分

几何尺寸、材料属性和物理模型如表1所示。

Table 1 Geometry,material properties and physical model of the case model 1表1 案例模型1的几何尺寸、材料属性和物理模型

将网格化模型分别用Fluent和YH-ACT进行数值求解，得到在x=0截面上的温度云图。二者在x=0截面上的温度分布结果几乎相同，取出口截面中心点温度和Fluent结果进行对比，如表2所示。由表2可知，在出口中心点温度对比中，YH-ACT的计算结果比目标值偏小，相对误差为0.3%。

(2)案例模型2：圆管湍流压降，如图10所示。

计算域为半径R= 0.00125 m，长度L=0.1 m的圆管道。网格划分为非结构网格，网格数量总计为4 500个。由于管道非常细长，显示时长度方向经过缩放后网格如图11所示。

Table 2 Temperature at the center of the exit section of the case model 2表2 案例模型2的出口截面中心点温度对比

Figure 10 Fluid computing domain of the case model 2图10 案例模型2的流体计算域

Figure 11 Meshing of the case model 2图11 案例模型2的网格划分

几何尺寸、材料属性和物理模型如表3所示。

Table 3 Geometry,material properties and physical model of the case model 2表3 案例模型2的几何尺寸、材料属性和物理模型

将网格化模型分别用YH-ACT和Fluent软件进行数值求解，得到在zy方向上x=0截面上的压力云图。二者在x=0截面上压力云图结果差距非常小，取进出口截面压差结果进行对比，如表4所示可知，在进出口截面压差上YH-ACT的计算结果相比Fluent计算结果相对误差为1.1%。

(3)案例模型3：T型管中的层流流动，如图12所示。

计算域为L=3.0 m，W=1.0 m，深度H=1.0 m的T型管。网格划分为结构网格，网格数量总计为72 000个。网格如图13所示。

Table 4 Pressures of inlet and outlet of the case model 2表4 案例模型2的进出口截面压差对比

Figure 12 Fluid computation domain of the case model 3图12 案例模型3的流体计算域

Figure 13 Meshing of the case model 3图13 案例模型3的网格划分

几何尺寸、材料属性和物理模型如表5所示。

Table 5 Geometry,material properties and physical model of the case model 3表5 案例模型3的几何尺寸、材料属性和物理模型

将网格化模型分别用YH-ACT和Fluent软件进行数值求解，得到在z=0截面上的y方向速度分量云图。YH-ACT与Fluent在z=0截面上的y方向速度分量结果非常接近，取左边直管中心线上y方向速度分量与Fluent计算结果相比，如图14所示，最大相对误差为1.4%。

Figure 14 The y-velocity at z=0图14 z=0处y分量速度随y坐标变化曲线

6.2 性能测试

本节选取案例3中的模型进行网格加密后对YH-ACT进行并行性能测试，计算时每个结点用24个进程，最大计算规模达到400个结点共9 600个进程。结点数为17 664 643，四面体单元个数为103 038 583。

模型稳态情况下测试结果如表6和图15所示，共100次迭代，最大计算规模达到9 600进程。考虑程序主要耗时部分CFD Solver的计算时间，相比单结点24个进程加速比为111.7，并行效率为27.9%。

Table 6 Speedup and parallel efficiency (steady state)表6 模型(稳态)的加速比及并行效率

Figure 15 Speedup and parallel efficiency (steady state)图15 模型(稳态)测试的加速比及并行效率

模型瞬态情况下测试结果如表7和图16所示，共100次迭代。从结果可以看出，最大计算规模达到9 600进程，相比单结点24个进程程序的加速比为37.2，并行效率为9.3%。

可以看到，YH-ACT有较好的可扩展性。就稳态模型而言，当测试规模达到9 600个进程时，并行效率依然有28%左右。而同等网格规模的瞬态模型的并行效率会相对低一些，这与单个时间步的迭代次数设置及先后时间步之间的数据依赖等因素有关。

Table 7 Speedup and parallel efficiency (transient state)表7 模型(瞬态)的加速比及并行效率

Figure 16 Speedup and parallel efficiency (transient state)图16 模型(瞬态)测试的加速比及并行效率

7 结束语

本文详细描述了一种热工流体力学并行应用程序YH-ACT的设计、实现方案和测试结果。通过对程序的验证和测试可以看出YH-ACT具备与商业CFD软件相当的计算能力和精度，并行性能较好，具备大型超算平台并行拓展能力，能够满足局部精细热工水力分析的需求。