波形板干燥器內液滴模擬的GPU加速研究

張 迪，黃 偉，羅 琦，王 侃

（1．清華大學工程物理系，北京 100084；

2．中國核動力研究設計院核反應堆系統(tǒng)設計技術重點實驗室，四川成都 610041；

3．中國核動力研究設計院，四川成都 610041）

波形板干燥器內液滴模擬的GPU加速研究

張 迪1，2，黃 偉2，羅 琦3，王 侃1

（1．清華大學工程物理系，北京 100084；

2．中國核動力研究設計院核反應堆系統(tǒng)設計技術重點實驗室，四川成都 610041；

3．中國核動力研究設計院，四川成都 610041）

基于CUDA C，使用GPU模擬了波形板內大量液滴的二維運動。根據GPU硬件特性，優(yōu)化了模擬算法。在CPU上實現相同算法，對比了單一液滴直徑、不同液滴數量下GPU與CPU的計算時間。結果表明：液滴數量越大，GPU并行加速效果越明顯，在液滴數量較大時，得到約60倍的加速。用GPU分析了典型液滴直徑分布的情況，得到約30倍的加速。

波形板干燥器；液滴運動模擬；GPU并行加速；CUDA C

波形板干燥器是蒸汽發(fā)生器內保證蒸汽品質的關鍵設備。分析其內部液滴行為、分離效率，對其設計具有重要意義。分析的方法主要有實驗、機理模型分析、CFD方法。CFD方法的優(yōu)勢在于成本低，可了解波形板內流場分布、液滴運動軌跡等細節(jié)。目前，波形板內兩相流的CFD分析主要采用歐拉－拉格朗日方法［1－4］：連續(xù)相流場采用歐拉方法，并簡化為二維流動；離散相的模擬采用拉格朗日方法，對液滴的運動軌跡進行模擬。由于液滴直徑分布范圍廣［5］，入射位置不同，需模擬大量的液滴，方能得到準確的結果，導致計算量很大。而液滴的運動遵從相同的物理準則，對大量液滴的模擬具有良好的并行性，適合于GPU的大規(guī)模并行運算。

基于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足夠大時，對應標準阻力曲線的牛頓區(qū)，CD下降為0．44，但CD不應小于0．44，故有：

1．2 液滴運動方程求解

確定液滴受力情況后，根據牛頓第二定律，分析液滴的運動：

化簡式（4）可得：

其中，τp為粒子松弛時間。

τp衡量液滴受周圍流體作用而改變運動狀態(tài)的快慢［11］。τp越大，液滴受到氣相的影響越慢，更易保持原有運動路徑，撞到波形板側壁被分離。從式（6）可看出，直徑越大的液滴越易被分離。

液滴位移公式為：

對于式（5），顯式的數值求解方法在時間步較大時會不穩(wěn)定，而隱式的方法需進行迭代求解，效率較低，因此采用積分法直接求解。在較短的時間步內，τp與uf近似保持不變，則對式（5）求積分可得：

給定液滴初始位置、速度，根據式（8）和（9）即可得到液滴的運動軌跡。

計算液滴所受曳力，需要液滴所處位置的氣相流速、氣相黏性系數等參數。這些氣相參數由CFD方法模擬板間二維流場得到。由于CFD計算結果局限于有限的網格點上，需通過插值獲取液滴所處位置的氣相參數。首先計算液滴所處位置附近的網格編號（I，J）；然后讀取編號為（I，J）的網格處氣相參數，并通過二階精度的泰勒展開插值獲取液滴所處位置的氣相參數。對于氣相的x軸方向速度，有：

2 GPU與CUDA概述

GPU（圖形處理器）是一種專用圖像處理設備。GPU可與CPU一樣用作通用計算。相較于CPU，GPU將更多的晶體管用作計算器，只將很少的晶體管用作控制器和緩存。因此GPU更擅長于浮點數運算，擅長數據并行類的通用計算，即相同的程序在大量的數據上執(zhí)行，代數計算的強度遠高于流程控制、內存訪問。

早期的GPU通用計算基于OpenGL、DirectX等專用圖像處理語言，用戶需將實際問題轉換為圖像處理問題。為解決此不便，NVIDIA公司推出一通用并行計算架構，實現了異構資源下的大規(guī)模并行計算，即CUDA。CUDA并行架構可概括為SIMT，即多個線程對不同數據并發(fā)執(zhí)行相同指令，線程具有“節(jié)點－塊－線程”的層級結構。1個GPU上有多個流多處理器，每個流多處理器由多個標量處理核構成。單個線程運行于標量處理核，多個線程共同構成1個塊，運行于1個流多處理器，多個塊構成1個節(jié)點運行于1個GPU。在CUDA架構下，用戶可使用類似C語言的CUDA C，極大地方便了編程。

液滴的運動遵從相同的物理準則，計算時相同計算流程應用于大量液滴，且計算量主要在于求解液滴運動方程的浮點運算，流程控制操作相對很少，因此大量液滴并行計算為數據并行類，適合于GPU并行加速。

3 算法實現與優(yōu)化

3．1 計算流程

圖2示出液滴模擬的計算流程。變量I記錄了程序計算到的液滴組的編號。

圖2 液滴模擬計算流程Fig．2 Droplet simulation process

主機負責計算流程的控制以及板間流場的CFD模擬。由于不考慮液滴對氣相流場的影響，氣相流場只需計算1次，然后將氣相流場的速度等數據傳遞到GPU設備內存。由于入射液滴直徑具有分布，并非單一數值，將入射液滴分為Ng組，1組內的液滴直徑相同。每次調用GPU只模擬1組液滴。GPU程序使用CUDA C實現，每個線程模擬1個液滴，圖2右側框內為1個線程的計算流程，多層框表示多個線程同時進行。計算時記錄液滴的運動軌跡，便于畫圖分析。計算結束后，統(tǒng)計逃出波形板的液滴數，以計算分離效率。

3．2 算法優(yōu)化

針對GPU并行計算的特點，對算法進行以下優(yōu)化。

1）使用單精度浮點數。GPU對于單精度浮點數的運算效率遠高于雙精度浮點數。程序中的實數變量均定義為單精度浮點數，數值實驗表明，單精度的計算結果滿足精度要求。

2）指令優(yōu)化。CUDA C對于數學運算指令做了優(yōu)化，對于精度要求不是特別高的數學運算，推薦使用內部函數提高運算效率。例如式（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℃時的空氣，使用商業(yè)軟件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上的程序進行了編譯器優(yōu)化。

定義加速倍數Acc為：

需指出的是，本文僅考慮液滴模擬的耗時，流場計算耗時不計入CPU計算時間，GPU計算時間包含GPU與主機數據交換的耗時。

采用4．1節(jié)的計算條件，考慮大量液滴均勻入射，對比不同入射液滴數N時，兩者的計算時間，結果如圖7所示。

可看出，CPU計算時間與液滴數呈正比。而GPU計算時間介于5～18ms之間，隨液滴數增大而緩慢增加。GPU加速倍數隨液滴數增大而增大。這是由于液滴越多，并行的線程越多，GPU的計算資源越易發(fā)揮。在液滴數量為8 000左右時，GPU的計算速度約為CPU的60倍。

為進一步驗證加速效果，取入射液滴數為8 192，研究不同液滴直徑、不同氣相入口流速下的加速效果。圖8示出了計算時間和加速倍數，可看出，隨著液滴直徑增大、氣相流速增大，CPU與GPU計算時間均減小，這是因為液滴直徑增大、氣相流速增大使液滴更快分離，模擬的總路徑變短。而兩者均得到約60倍的加速。

由式（6）可知，液滴直徑不同，松弛時間不同，從而影響波形板的分離效率。因此，精確的模擬需要入射液滴的直徑分布。文獻［5］給出了空氣、水冷態(tài)試驗時，液滴直徑的分布。將此分布近似分為8組，每組的直徑、液滴數百分比如圖9所示。

基于圖3所示波形板，研究不同入口流速下的分離效率。模擬1萬個液滴，其直徑服從圖9所示的分布。GPU和CPU計算時間的對比如圖10所示。

由圖10可見，計算時間隨著流速增大而減小，這是因為流速增大，分離效率增大，逃出波形板的液滴數減少，模擬的總路徑變短；加速效果不如單直徑液滴模擬，這可能是因為某些組的液滴數量較少，造成硬件資源使用不充分，但仍有30倍左右的加速。

圖9 不同直徑的液滴數分布Fig．9 Droplet number distribution under different diameters

5 結論與展望

根據GPU并行計算的特點，基于CUDA C編寫程序并優(yōu)化，在GPU上實現了典型波形板內大量液滴二維運動的模擬，得到了液滴在波形板內的運動路徑。相對于CPU上實現的相同算法，在單一液滴直徑時有60倍左右加速，在考慮液滴直徑分布時有30倍左右的加速。

圖10 液滴具有直徑分布時GPU與CPU計算時間對比Fig．10 Comparison between GPU and CPU computing time with Droplets’diameter distribution

圖6顯示，大量液滴在較短的路程內被分離，相應的線程停止計算，計算量較小；少量液滴運動到波形板出口，相應的線程計算量大，這導致線程負載不平衡，可能會減弱加速效果，后續(xù)工作應解決這一問題。此外，本文著重于算法層面的加速，應指出，如果要得到準確的分離效率，還需考慮曳力模型的選取、液滴的二次攜帶、液滴對氣相流場的影響、三維計算等問題。這勢必會增大計算量，延長計算時間，而GPU并行加速可解決這一問題，使壓水堆蒸汽發(fā)生器內波形板干燥器的設計、性能分析更加方便、快捷。

［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］．北京：機械工業(yè)出版社，2011．

［13］《蒸汽發(fā)生器》編寫組．蒸汽發(fā)生器［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—），男，河北晉州人，博士研究生，從事反應堆系統(tǒng)與設備研究