波形板干燥器内液滴模拟的GPU加速研究

张 迪，黄 伟，罗 琦，王 侃

（1．清华大学工程物理系，北京 100084；

2．中国核动力研究设计院核反应堆系统设计技术重点实验室，四川成都 610041；

3．中国核动力研究设计院，四川成都 610041）

波形板干燥器内液滴模拟的GPU加速研究

张 迪1，2，黄 伟2，罗 琦3，王 侃1

（1．清华大学工程物理系，北京 100084；

2．中国核动力研究设计院核反应堆系统设计技术重点实验室，四川成都 610041；

3．中国核动力研究设计院，四川成都 610041）

基于CUDA C，使用GPU模拟了波形板内大量液滴的二维运动。根据GPU硬件特性，优化了模拟算法。在CPU上实现相同算法，对比了单一液滴直径、不同液滴数量下GPU与CPU的计算时间。结果表明：液滴数量越大，GPU并行加速效果越明显，在液滴数量较大时，得到约60倍的加速。用GPU分析了典型液滴直径分布的情况，得到约30倍的加速。

波形板干燥器；液滴运动模拟；GPU并行加速；CUDA C

波形板干燥器是蒸汽发生器内保证蒸汽品质的关键设备。分析其内部液滴行为、分离效率，对其设计具有重要意义。分析的方法主要有实验、机理模型分析、CFD方法。CFD方法的优势在于成本低，可了解波形板内流场分布、液滴运动轨迹等细节。目前，波形板内两相流的CFD分析主要采用欧拉－拉格朗日方法［1－4］：连续相流场采用欧拉方法，并简化为二维流动；离散相的模拟采用拉格朗日方法，对液滴的运动轨迹进行模拟。由于液滴直径分布范围广［5］，入射位置不同，需模拟大量的液滴，方能得到准确的结果，导致计算量很大。而液滴的运动遵从相同的物理准则，对大量液滴的模拟具有良好的并行性，适合于GPU的大规模并行运算。

基于CUDA C的GPU并行计算是各领域数值计算的研究热点。在核工程领域，文献［6］使用GPU并行加速MOC算法，文献［7］使用GPU加速大量中子的模拟；对于粒子、液滴等的模拟，文献［8］使用GPU模拟SPH方法中的粒子行为，文献［9］使用GPU模拟流场中的固体颗粒，文献［10］使用GPU模拟自由表面流的MPS方法中的粒子，均获得良好的加速效果。很少有文献用GPU加速模拟波形板内的液滴，本文探讨其可行性。

1 液滴运动数值模拟

图1示出典型无钩波形板以及液滴在其中的运动轨迹。在离心力、曳力等作用下，液滴触壁被分离，达到干燥的效果。

图1 典型波形板及液滴运动轨迹示意图Fig．1 Scheme of typical wave－type vane and droplets’motion trails

液滴的模拟基于拉格朗日方法：将液滴看作质点，分析液滴的受力，记录运动过程中的速度、位移，以确定液滴轨迹。本文忽略液滴对气相流场的影响，不考虑液滴撞壁后的二次携带，即认为液滴碰到壁面就被分离，这适用于液滴密度低、入口流速低于临界流速的情况。

1．1 液滴受力分析

液滴密度远大于气相流体的密度，因此可忽略附加质量力、压力梯度力等，模拟液滴的二维运动，也不考虑重力和浮升力。因此，仅考虑液滴所受曳力。液滴近似为球形，其粒子雷诺数Rep为：

曳力系数CD与Rep有关，采用Schiller－Naumann曳力模型［5］计算：

考虑当Rep足够大时，对应标准阻力曲线的牛顿区，CD下降为0．44，但CD不应小于0．44，故有：

1．2 液滴运动方程求解

确定液滴受力情况后，根据牛顿第二定律，分析液滴的运动：

化简式（4）可得：

其中，τp为粒子松弛时间。

τp衡量液滴受周围流体作用而改变运动状态的快慢［11］。τp越大，液滴受到气相的影响越慢，更易保持原有运动路径，撞到波形板侧壁被分离。从式（6）可看出，直径越大的液滴越易被分离。

液滴位移公式为：

对于式（5），显式的数值求解方法在时间步较大时会不稳定，而隐式的方法需进行迭代求解，效率较低，因此采用积分法直接求解。在较短的时间步内，τp与uf近似保持不变，则对式（5）求积分可得：

给定液滴初始位置、速度，根据式（8）和（9）即可得到液滴的运动轨迹。

计算液滴所受曳力，需要液滴所处位置的气相流速、气相黏性系数等参数。这些气相参数由CFD方法模拟板间二维流场得到。由于CFD计算结果局限于有限的网格点上，需通过插值获取液滴所处位置的气相参数。首先计算液滴所处位置附近的网格编号（I，J）；然后读取编号为（I，J）的网格处气相参数，并通过二阶精度的泰勒展开插值获取液滴所处位置的气相参数。对于气相的x轴方向速度，有：

2 GPU与CUDA概述

GPU（图形处理器）是一种专用图像处理设备。GPU可与CPU一样用作通用计算。相较于CPU，GPU将更多的晶体管用作计算器，只将很少的晶体管用作控制器和缓存。因此GPU更擅长于浮点数运算，擅长数据并行类的通用计算，即相同的程序在大量的数据上执行，代数计算的强度远高于流程控制、内存访问。

早期的GPU通用计算基于OpenGL、DirectX等专用图像处理语言，用户需将实际问题转换为图像处理问题。为解决此不便，NVIDIA公司推出一通用并行计算架构，实现了异构资源下的大规模并行计算，即CUDA。CUDA并行架构可概括为SIMT，即多个线程对不同数据并发执行相同指令，线程具有“节点－块－线程”的层级结构。1个GPU上有多个流多处理器，每个流多处理器由多个标量处理核构成。单个线程运行于标量处理核，多个线程共同构成1个块，运行于1个流多处理器，多个块构成1个节点运行于1个GPU。在CUDA架构下，用户可使用类似C语言的CUDA C，极大地方便了编程。

液滴的运动遵从相同的物理准则，计算时相同计算流程应用于大量液滴，且计算量主要在于求解液滴运动方程的浮点运算，流程控制操作相对很少，因此大量液滴并行计算为数据并行类，适合于GPU并行加速。

3 算法实现与优化

3．1 计算流程

图2示出液滴模拟的计算流程。变量I记录了程序计算到的液滴组的编号。

图2 液滴模拟计算流程Fig．2 Droplet simulation process

主机负责计算流程的控制以及板间流场的CFD模拟。由于不考虑液滴对气相流场的影响，气相流场只需计算1次，然后将气相流场的速度等数据传递到GPU设备内存。由于入射液滴直径具有分布，并非单一数值，将入射液滴分为Ng组，1组内的液滴直径相同。每次调用GPU只模拟1组液滴。GPU程序使用CUDA C实现，每个线程模拟1个液滴，图2右侧框内为1个线程的计算流程，多层框表示多个线程同时进行。计算时记录液滴的运动轨迹，便于画图分析。计算结束后，统计逃出波形板的液滴数，以计算分离效率。

3．2 算法优化

针对GPU并行计算的特点，对算法进行以下优化。

1）使用单精度浮点数。GPU对于单精度浮点数的运算效率远高于双精度浮点数。程序中的实数变量均定义为单精度浮点数，数值实验表明，单精度的计算结果满足精度要求。

2）指令优化。CUDA C对于数学运算指令做了优化，对于精度要求不是特别高的数学运算，推荐使用内部函数提高运算效率。例如式（3）的指数运算Re0．687p耗时较长，用CUDA C内部函数＿powf（）代替，数值实验证明内部函数的精度足够。

文献［12］指出，存储器带宽是GPU计算的瓶颈，应尽量使用读写速度极快的共享内存。计算液滴所处位置附近的网格编号时，需反复读取网格坐标数据，将网格坐标数据储存在共享内存上，提高反复读取的速度。

4 算法验证与数值实验

4．1 算法验证

本文使用的GPU型号为NVIDIA Quadro 600，它具有2个流多处理器，每个处理器上有48个CUDA核，1．28GHz，1GB设备显存，计算能力为2．1。分析文献［13］给出的典型的波形板干燥器，其几何尺寸如图3所示。

图3 波形板二维几何图Fig．3 2Dgeometry of wave－type vane

气相取为25℃时的空气，使用商业软件ANSYS CFX计算气相二维流场。边界条件为：入口流速为1m／s，出口压力设为大气压，壁面为无滑移固壁，湍流模型取标准k－ε湍流模型。本文着重于液滴的模拟，气相流场计算结果不在此赘述。离散相为液滴，液滴直径取为40μm，入射速度与气相流速均为1m／s，在干燥器入口截面上均匀入射。

计算网格数为310×15（mesh1）、458×22（mesh2）、700×34（mesh3）的3套网格下，中央入射液滴轨迹，以确保网格无关性。图4示出计算结果，三者轨迹接近，因此取网格数为458×22可确保网格无关性。

为验证GPU模拟结果的准确性，使用商用软件ANSYS CFX的计算结果进行对比。两者采用相同的曳力模型，模拟由波形板中央进入的液滴，记录液滴运动轨迹，结果如图5所示。两者运动轨迹几乎完全重合，说明GPU模拟结果准确。

图4 网格无关解检验Fig．4 Check of mesh independent solution

图5 GPU与ANSYS CFX计算结果对比Fig．5 Comparison between GPU results and ANSYS CFX results

进一步模拟50个液滴的运动轨迹，以考察GPU对多个液滴的模拟效果，结果如图6所示。大部分液滴在前两个转弯处碰到壁面被分离，这与文献［13］的结论相同，进一步说明了GPU模拟的正确性。第3个拐角后，液滴几乎不再触壁，大量液滴的运动轨迹逐渐聚拢，导致整体轨迹逐渐变细。

图6 大量液滴的运动轨迹Fig．6 Trajectory of droplets

4．2 数值实验

为验证GPU算法相对于CPU的加速效果，在CPU上实现完全相同的算法。CPU为Intel i5－2320单核，3．0GHz，计算机有4GB内存，CPU上的程序进行了编译器优化。

定义加速倍数Acc为：

需指出的是，本文仅考虑液滴模拟的耗时，流场计算耗时不计入CPU计算时间，GPU计算时间包含GPU与主机数据交换的耗时。

采用4．1节的计算条件，考虑大量液滴均匀入射，对比不同入射液滴数N时，两者的计算时间，结果如图7所示。

可看出，CPU计算时间与液滴数呈正比。而GPU计算时间介于5～18ms之间，随液滴数增大而缓慢增加。GPU加速倍数随液滴数增大而增大。这是由于液滴越多，并行的线程越多，GPU的计算资源越易发挥。在液滴数量为8 000左右时，GPU的计算速度约为CPU的60倍。

为进一步验证加速效果，取入射液滴数为8 192，研究不同液滴直径、不同气相入口流速下的加速效果。图8示出了计算时间和加速倍数，可看出，随着液滴直径增大、气相流速增大，CPU与GPU计算时间均减小，这是因为液滴直径增大、气相流速增大使液滴更快分离，模拟的总路径变短。而两者均得到约60倍的加速。

由式（6）可知，液滴直径不同，松弛时间不同，从而影响波形板的分离效率。因此，精确的模拟需要入射液滴的直径分布。文献［5］给出了空气、水冷态试验时，液滴直径的分布。将此分布近似分为8组，每组的直径、液滴数百分比如图9所示。

基于图3所示波形板，研究不同入口流速下的分离效率。模拟1万个液滴，其直径服从图9所示的分布。GPU和CPU计算时间的对比如图10所示。

由图10可见，计算时间随着流速增大而减小，这是因为流速增大，分离效率增大，逃出波形板的液滴数减少，模拟的总路径变短；加速效果不如单直径液滴模拟，这可能是因为某些组的液滴数量较少，造成硬件资源使用不充分，但仍有30倍左右的加速。

图9 不同直径的液滴数分布Fig．9 Droplet number distribution under different diameters

5 结论与展望

根据GPU并行计算的特点，基于CUDA C编写程序并优化，在GPU上实现了典型波形板内大量液滴二维运动的模拟，得到了液滴在波形板内的运动路径。相对于CPU上实现的相同算法，在单一液滴直径时有60倍左右加速，在考虑液滴直径分布时有30倍左右的加速。

图10 液滴具有直径分布时GPU与CPU计算时间对比Fig．10 Comparison between GPU and CPU computing time with Droplets’diameter distribution

图6显示，大量液滴在较短的路程内被分离，相应的线程停止计算，计算量较小；少量液滴运动到波形板出口，相应的线程计算量大，这导致线程负载不平衡，可能会减弱加速效果，后续工作应解决这一问题。此外，本文着重于算法层面的加速，应指出，如果要得到准确的分离效率，还需考虑曳力模型的选取、液滴的二次携带、液滴对气相流场的影响、三维计算等问题。这势必会增大计算量，延长计算时间，而GPU并行加速可解决这一问题，使压水堆蒸汽发生器内波形板干燥器的设计、性能分析更加方便、快捷。

［1］ LI J，HUANG S，WANG X．Numerical study of steam－water separators with wave－type vanes［J］．Chin J Chem Eng，2007，15（4）：492－498．

［2］ DOSTAL V，TAKAHASHI M．Modeling and optimization of steam dryer for removal of leadbismuth droplets［J］．Progress in Nuclear Energy， 2008，50：631－637．

［3］ KIM J I，KIM M Y，BAE H S，et al．The performance of moisture separation system of a steam generator using computational fluid dynamics［C］∥Proceedings of ICAP’05．Seoul，Korea：［s．n．］，2005．

［4］ 李嘉．波形板汽水分离器的理论与实验研究［D］．武汉：华中科技大学，2007．

［5］ KOLEV N I．Multiphase flow dynamics，4：Nuclear thermal hydraulics［M］．Germany：Springer，2006．

［6］ 张知竹，李庆，王侃．GPU加速三维特征线方法的研究［J］．核动力工程，2013，34（S1）：18－23．

ZHANG Zhizhu，LI Qing，WANG Kan．Study on acceleration of three－dimensional method of characteristics by GPU［J］．Nuclear Power Engineering，2013，34（S1）：18－23（in Chinese）．

［7］ XU Qi，YU Ganglin，WANG Kan．GPU－accelerated 3Dneutron diffusion code based on finite difference method［C］∥Proceedings of PHYSOR 2012Conference．Knoxville，Tennessee，USA：American Nuclear Society，2012．

［8］ HERAULT A，BILOTTA G，DALRYMPLE R A．SPH on GPU with CUDA［J］．Journal of Hydraulic Research，2010，48：74－79．

［9］ XU Ming，CHEN Feiguo，LIU Xinhua，et al．Discrete particle simulation of gas－solid twophase flows with multi－scale CPU－GPU hybrid computation［J］．Chemical Engineering Journal，2012，207－208：746－757．

［10］HORI C，GOTOH H，IKARI H，et al．GPU－acceleration for moving particle semi－implicit method［J］．Computers ＆Fluids，2011，51：174－183．

［11］郭烈棉．两相与多相动力学［M］．西安：西安交通大学出版社，2002．

［12］仇德元．GPGPU编程技术——从GLSL、CUDA到OpenCL［M］．北京：机械工业出版社，2011．

［13］《蒸汽发生器》编写组．蒸汽发生器［M］．北京：原子能出版社，1982．

GPU Accelerated Droplets Simulation in Wave－type Vane Dryer

ZHANG Di1，2，HUANG Wei2，LUO Qi3，WANG Kan1
（1．Department of Engineering Physics，Tsinghua University，Beijing100084，China；
2．Science and Technology on Reactor System Design Technology Laboratory，
Nuclear Power Institute of China，Chengdu610041，China；
3．Nuclear Power Institute of China，Chengdu610041，China）

Droplets’two－dimensional dynamic simulation was performed by GPU based on CUDA C．The simulation algorithm was optimized considering hardware characteristics of GPU parallel computation．The same algorithm was carried out by CPU．The computing time of GPU was compared to that of CPU for the cases with single droplet diameter but different droplet numbers．The results show that the larger the number of water droplets is，the more obvious the GPU parallel acceleration effect is．In the case with large number of water droplets，about 60times speedup is achieved．And the computing time of GPU for the case with a typical distribution of multiple droplet diameters obtains about 30times speedup．

wave－type vane dryer；droplet dynamic simulation；GPU paraller acceleration；CUDA C

TL124

：A

：1000－6931（2015）03－0491－06

10．7538／yzk．2015．49．03．0491

2013－12－13；

2014－03－08

张 迪（1988—），男，河北晋州人，博士研究生，从事反应堆系统与设备研究