GPU加速的球床式氟鹽冷卻高溫堆堆芯熱工水力程序研發(fā)

鄂彥志 鄒 楊 郭 威 彭 玉 徐洪杰

1（中國(guó)科學(xué)院上海應(yīng)用物理研究所 嘉定園區(qū) 上海 201800）2（中國(guó)科學(xué)院大學(xué) 北京 100049）

GPU加速的球床式氟鹽冷卻高溫堆堆芯熱工水力程序研發(fā)

鄂彥志1,2鄒 楊1郭 威1彭 玉1,2徐洪杰1

1（中國(guó)科學(xué)院上海應(yīng)用物理研究所 嘉定園區(qū) 上海 201800）2（中國(guó)科學(xué)院大學(xué) 北京 100049）

球床式氟鹽冷卻高溫堆(Pebble Bed Fluoride-salt Cooled High Temperature Reactor, PB-FHR)是一種先進(jìn)的第四代反應(yīng)堆。三維堆芯熱工水力程序能夠模擬具有復(fù)雜空間效應(yīng)的工況，但計(jì)算耗時(shí)較高。圖形處理器(Graphics Processing Unit, GPU)具有大量計(jì)算單元，可有效提高程序的計(jì)算速度。本文研發(fā)了GPU加速的PB-FHR堆芯熱工水力程序(GPU-accelerated Thermal Hydraulic Code, GATH)，采用非熱平衡多孔介質(zhì)模型建立堆芯物理模型，研究并實(shí)現(xiàn)了GPU高速求解算法。對(duì)PB-FHR的堆芯模型進(jìn)行了熱工水力分析，與商用計(jì)算流體力學(xué)軟件ANSYS CFX的計(jì)算結(jié)果進(jìn)行了對(duì)比，驗(yàn)證了程序的正確性。GPU加速性能分析的結(jié)果表明，程序整體的加速比率可達(dá)8.39倍，證明所研發(fā)的GPU求解算法能有效提升堆芯熱工水力分析的計(jì)算效率。關(guān)鍵詞 球床式氟鹽冷卻高溫堆，GPU加速，熱工水力，非熱平衡多孔介質(zhì)

球床式氟鹽冷卻高溫堆(Pebble Bed Fluoridesalt Cooled High Temperature Reactor, PB-FHR)采用球形燃料元件及液態(tài)氟化鹽冷卻劑，是一種先進(jìn)的第四代反應(yīng)堆，能夠在高溫、低壓的工況下運(yùn)行，具有良好的安全及經(jīng)濟(jì)特性[1]。目前世界上已有若干科研機(jī)構(gòu)對(duì)PB-FHR做出了相關(guān)的研究。德國(guó)代爾夫特理工大學(xué)于2005年提出了一種球床熔鹽堆的概念設(shè)計(jì)[2]；美國(guó)伯克利大學(xué)于2014年提出了一種小型模塊化PB-FHR的概念設(shè)計(jì)[3]；中國(guó)科學(xué)院上海應(yīng)用物理研究所于2011年提出了一種PB-FHR小型實(shí)驗(yàn)堆設(shè)計(jì)[4]。

PB-FHR堆芯是一種典型的多孔介質(zhì)，球形燃料元件隨機(jī)堆積在PB-FHR的堆芯內(nèi)，形成了充滿空隙的球床區(qū)域，冷卻劑通過(guò)球床的空隙流出堆芯。堆芯熱工水力程序是反應(yīng)堆設(shè)計(jì)及安全分析的重要工具，用于計(jì)算堆芯流場(chǎng)及溫度信息，為中子物理模擬提供溫度反饋。然而PB-FHR目前并沒(méi)有專(zhuān)用的堆芯熱工水力程序，若干研究機(jī)構(gòu)通過(guò)對(duì)現(xiàn)有的熱工水力程序進(jìn)行二次開(kāi)發(fā)，使其能夠用于PB-FHR。代爾夫特理工大學(xué)對(duì)高溫氣冷堆的二維熱工水力分析程序THERMIX進(jìn)行了二次開(kāi)發(fā)并應(yīng)用于PB-FHR[2]。伯克利大學(xué)將熔鹽物性植入系統(tǒng)分析程序RELAP中，用于PB-FHR的系統(tǒng)熱工水力分析[5]。中國(guó)科學(xué)院上海應(yīng)用物理研究所使用商用計(jì)算流體力學(xué)軟件FLUENT建立PB-FHR堆芯多孔介質(zhì)模型，利用軟件提供的二次開(kāi)發(fā)接口添加熔鹽物性及換熱公式，進(jìn)行了PB-FHR三維堆芯熱工水力分析[6]。西安交通大學(xué)采用商用計(jì)算流體力學(xué)軟件ANSYS對(duì)PB-FHR球床中的燃料球進(jìn)行精細(xì)建模及模擬，分析了熔鹽在球床孔隙及燃料球表面的熱工水力特性[7-8]。堆芯燃料球的精確建模及模擬能夠提供最詳細(xì)的熱工水力信息，但需使用大量的網(wǎng)格對(duì)模型進(jìn)行劃分，導(dǎo)致求解非常耗時(shí)。因此，球床堆的全堆芯三維熱工水力模擬通常采用多孔介質(zhì)模型，而球床精確建模方法更適用于堆芯的局部分析。與二維計(jì)算程序及系統(tǒng)分析程序相比，三維熱工水力程序能模擬具有復(fù)雜空間效應(yīng)的工況，如堆芯功率分布非對(duì)稱(chēng)等情況，但計(jì)算耗時(shí)仍然較高，尤其在耦合中子物理程序及瞬態(tài)計(jì)算的情況。

從2006年英偉達(dá)公司推出CUDA (Compute Unified Device Architecture)運(yùn)算架構(gòu)開(kāi)始，圖形處理器(Graphics Processing Unit, GPU)逐步成為了大規(guī)模并行數(shù)值計(jì)算的重要硬件工具[9]。GPU具有大量的計(jì)算單元，將可并行的計(jì)算任務(wù)加載到GPU，可有效提高程序的計(jì)算速度。CUDA提供了GPU的硬件接口，為開(kāi)發(fā)者提供了通用的GPU計(jì)算模型，有效推動(dòng)了GPU加速技術(shù)在流體力學(xué)、材料科學(xué)、核科學(xué)等科研領(lǐng)域的應(yīng)用。例如，Bailey等[10]使用GPU將模擬流場(chǎng)的格子玻爾茲曼模擬算法加速了28倍，Boyd等[11]使用GPU將中子輸運(yùn)程序OpenMOC加速了25-35倍。

為PB-FHR研發(fā)專(zhuān)用、高效的三維堆芯熱工水力程序能有效推動(dòng)PB-FHR的研究及設(shè)計(jì)工作。本文將以計(jì)算流體力學(xué)方法為基礎(chǔ)，針對(duì)PB-FHR的熱工水力特性，研發(fā)GPU加速的堆芯熱工水力程序(GPU-accelerated Thermal Hydraulic Code, GATH)。本文采用非熱平衡多孔介質(zhì)模型[12]建立堆芯固體結(jié)構(gòu)及冷卻劑的物理模型；采用SIMPLEC (Semi-implicit Method for Pressure Linked Equation-Consistent)算法[13]求解流場(chǎng)及壓力場(chǎng)；為更準(zhǔn)確地模擬球床內(nèi)流體的熱彌散現(xiàn)象，引入多孔介質(zhì)的湍流模型；為了在大規(guī)模網(wǎng)格情況下高速求解堆芯熱工水力問(wèn)題，基于krylov子空間迭代算法及方程預(yù)處理算法，研究并實(shí)現(xiàn)了GPU高速求解算法；通過(guò)對(duì)PB-FHR堆芯模型進(jìn)行熱工水力分析，并與商用計(jì)算流體力學(xué)軟件ANSYS CFX的計(jì)算結(jié)果進(jìn)行對(duì)比，驗(yàn)證程序的正確性，同時(shí)詳細(xì)分析GPU的加速性能。

1 數(shù)理模型

本文以非熱平衡多孔介質(zhì)模型為基礎(chǔ)，將堆芯球床復(fù)雜的幾何近似為均勻的固相和液相介質(zhì)。固相介質(zhì)代表堆芯中的燃料球、反射層等固體結(jié)構(gòu)，液相介質(zhì)代表熔鹽冷卻劑。

1.1 冷卻劑數(shù)理模型

冷卻劑數(shù)理模型包括宏觀質(zhì)量、動(dòng)量、湍流動(dòng)能、耗散率及能量守恒方程。

熔鹽冷卻劑被作為不可壓縮流體處理，其質(zhì)量守恒方程如下：

式中：φ是多孔介質(zhì)孔隙率（液相體積分?jǐn)?shù)）；ρf是熔鹽密度；u為熔鹽的物理速度。

流體的流速通過(guò)求解宏觀動(dòng)量守恒方程獲得。該方程由微觀雷諾平均方程進(jìn)行體積平均化而獲得[14]。方程中包含由固相引起的阻力源項(xiàng)，具體表達(dá)式如下：

其中：

式中：Sij為宏觀雷諾應(yīng)力；ui、Ri、gi為i方向的速度、阻力源項(xiàng)、重力加速度；p、μ、k、μt代表壓力、速度、宏觀湍流動(dòng)能、湍流粘度；μt的表達(dá)式和普通流體相同：

式中：cμ為與湍流模型有關(guān)的無(wú)量綱數(shù)，一般取0.09；ε為湍流耗散率。k、ε和μt只在湍流區(qū)域存在，通過(guò)求解宏觀湍流動(dòng)能及耗散率方程獲得。Ri表達(dá)式采用了填充床中廣泛使用的Ergun公式[15]。

采用局部非熱平衡的多孔介質(zhì)模型，固相和液相的能量方程被分別建立，通過(guò)附加的對(duì)流換熱項(xiàng)耦合[16]。液相的溫度（能量）方程表達(dá)式如下：

式中：cpf、Tf、λf,eff及h代表比熱容、溫度、有效導(dǎo)熱系數(shù)及對(duì)流換熱系數(shù)；a為燃料球表面積體積比。λf,eff包含流體固有熱導(dǎo)率和湍流熱彌散效應(yīng)產(chǎn)生的導(dǎo)熱，本文采用Kuwahara等[17]提出的等效熱導(dǎo)率公式。

對(duì)流換熱項(xiàng)ha(Ts-Tf)代表流體從固體獲得的能量，本文采用Wakao[18]提出的對(duì)流換熱系數(shù)的經(jīng)驗(yàn)公式，該經(jīng)驗(yàn)公式在填充床中廣泛地應(yīng)用。

流體的等效粘度及等效導(dǎo)熱系數(shù)需根據(jù)湍流動(dòng)能及耗散率獲得，因此采用了N-K多孔介質(zhì)湍流模型[14]，多位研究人員證實(shí)該模型的結(jié)果適用于球型多孔介質(zhì)[19]，其表達(dá)式如下：

式中：Sk及Sε為N-K模型中的湍流動(dòng)能及耗散率生成率；c1、c2、σk和σε為湍流模型的無(wú)量綱常數(shù)。

1.2 堆芯固體數(shù)理模型

堆芯固體（球床、反射層等）的能量守恒方程為宏觀導(dǎo)熱方程，燃料被當(dāng)做固定在堆芯內(nèi)近似處理。堆芯內(nèi)固體的溫度通過(guò)求解該方程獲得，其表達(dá)式如下：

式中：ρs、cps、Ts、λs,eff代表固體的密度、比熱容、溫度及有效等效熱導(dǎo)率；qn為裂變產(chǎn)生的熱功率密度。方程中的對(duì)流換熱項(xiàng)是冷卻劑與燃料球、反射層換熱產(chǎn)生的冷卻功率密度。

其中，球床的等效導(dǎo)熱系數(shù)由Zehner-Bauer-Schlünder公式[20]給出，該公式在球床反應(yīng)堆中廣泛應(yīng)用。該公式考慮到所有傳熱機(jī)制：

式中：λs,cont為燃料球的接觸導(dǎo)熱；λs,f為空隙中滯止的流體導(dǎo)熱；λs,rad為球之間的輻射換熱等效熱導(dǎo)。

2 數(shù)值計(jì)算方法

2.1 離散方法

堆芯熱工水力物理量的控制方程（流場(chǎng)、溫度場(chǎng)等）可以表示成廣義的偏微分方程形式：

式中：ψ為待求物理量（流場(chǎng)或溫度場(chǎng)）；Γ為廣義擴(kuò)散系數(shù)；S為廣義源項(xiàng)。通用控制方程由瞬時(shí)項(xiàng)、對(duì)流項(xiàng)、擴(kuò)散項(xiàng)及源項(xiàng)組成。

采用有限體積法將偏微分方程在三維圓柱坐標(biāo)系的規(guī)則網(wǎng)格下進(jìn)行空間離散，采用交錯(cuò)網(wǎng)格存儲(chǔ)網(wǎng)格面上的流速。對(duì)流項(xiàng)采用一階迎風(fēng)格式，瞬時(shí)項(xiàng)采用一階向后差分格式。對(duì)通用控制方程進(jìn)行空間及時(shí)間上的積分后，便獲得離散方程：

式中：aP、aE、aW、aN、aS、aF、aB分別為網(wǎng)格中心、東側(cè)、西側(cè)、北側(cè)、南側(cè)、前側(cè)、后側(cè)的待求變量的離散系數(shù)；b為與源項(xiàng)有關(guān)的系數(shù)。三維空間網(wǎng)格的所有離散方程形成了一個(gè)大型七對(duì)角矩陣的線性方程組，需使用線性方程組的迭代求解算法計(jì)算方程組的解。

2.2 SIMPLEC算法

宏觀動(dòng)量方程是壓力與速度耦合的方程，采用SIMPLEC算法對(duì)其求解。SIMPLEC算法首先假設(shè)初始的壓力分布，求解動(dòng)量方程并獲得速度分布，然后根據(jù)質(zhì)量守恒方程推導(dǎo)出壓力修正方程，求解壓力修正方程獲得壓力修正值，并對(duì)速度及壓力進(jìn)行修正。整個(gè)過(guò)程反復(fù)迭代，直到速度及壓力收斂。x方向速度和壓力修正值的關(guān)系式如下：

式中：ae和anb為系數(shù)；u'e為表面上的速度；Ae為表面積；P'P和P'E為控制體積中心及東側(cè)控制體中心的壓力修正值，其余方向的關(guān)系式與其類(lèi)似。

2.3 求解流程

圖1顯示了堆芯熱工水力計(jì)算的求解流程。速度場(chǎng)及溫度場(chǎng)需要不斷地迭代求解，每次迭代求解后根據(jù)新的結(jié)果更新堆芯材料的熱物性，直到所有物理量收斂，以上過(guò)程又稱(chēng)外迭代過(guò)程。而求解物理量離散產(chǎn)生的線性方程組過(guò)程稱(chēng)為內(nèi)迭代過(guò)程。

3 GPU求解程序研發(fā)

三維堆芯熱工水力計(jì)算需使用大量網(wǎng)格，巨大的網(wǎng)格數(shù)量、多次的外迭代過(guò)程導(dǎo)致計(jì)算耗時(shí)很長(zhǎng)，其中最耗時(shí)的部分為離散的線性方程組的內(nèi)迭代過(guò)程，因此本文研究了基于GPU的迭代求解算法。

圖1 GATH程序求解流程Fig.1 Flowchart of solving thermal hydraulic procedure in GATH code.

方程預(yù)處理方法能進(jìn)一步加快求解過(guò)程的收斂速度。預(yù)處理算法的本質(zhì)為尋找線性方程組矩陣A的近似矩陣M，M≈A，使方程組Mx=b的求解速度很快。結(jié)合預(yù)處理算法的CG[21]與BICGSTAB[22]算法描述如下：

Conjugate Gradient Method

Equation: Ax=b

Initial guess of x: x0

r0=b-Ax0, z0=M-1r0, p0=z0

for j=0, 1,…, until convergence

αj=(rj, zj)/(Apj, pj)

xj+1=xj+αjpj

rj+1=rj-αjApj

zj+1=M-1rj+1

βj=(rj+1, zj+1)/(rj, zj)

pj+1=zj+1+βjpj

end for

BICGSTAB Method

Equation: Ax=b

Initial guess of x: x0

r0=b-Ax0, p0=r0

for j=0, 1,…, until convergence

qj=M-1pj

αj=(rj, r0)/(Aqj, r0)

sj=rj-αjAqj

zj=M-1sj

ωj=(Azj, sj)/(Azj, Azj)

xj+1=xj+αjqj+ωjzj

rj+1=sj-ωjAzj

βj=(rj+1, r0)/(rj, r0)×(αj/ωj)

pj+1=rj+1+βj(pj-ωjAqj)

end for

GPU求解器的兩種預(yù)處理算法為Neumann多項(xiàng)式算法(Neumann Polynomial Preconditioning, POLYN)[23]、幾何代數(shù)多重網(wǎng)格法(Geometric Algebraic Multigrid Method, GAMG)。

POLYN算法根據(jù)多項(xiàng)式展開(kāi)原理，將矩陣A-1用Neumann多項(xiàng)式近似表示，其計(jì)算公式如式(13)，本文采用三階多項(xiàng)式展開(kāi)。

本文提出了適用于規(guī)則網(wǎng)格的GAMG預(yù)處理算法，該算法不需要尋找矩陣A的近似矩陣。在求解方程組Ax=b時(shí)，該算法結(jié)合網(wǎng)格的空間信息及矩陣的代數(shù)信息，將求解過(guò)程化為由細(xì)到粗網(wǎng)格的迭代過(guò)程，通過(guò)在不同密度的網(wǎng)格進(jìn)行迭代，不同頻率的誤差分量得到有效衰減，加快了收斂速度。

該算法以標(biāo)準(zhǔn)的光滑聚集代數(shù)多重網(wǎng)格法[24]為基礎(chǔ)，利用網(wǎng)格的分布規(guī)律生成粗網(wǎng)，即空間上相鄰最近的27個(gè)網(wǎng)格聚集為一個(gè)粗網(wǎng)，這種方式雖然只適用于規(guī)則網(wǎng)格，但相比于標(biāo)準(zhǔn)方法速度更快。

GAMG預(yù)處理法分為建立過(guò)程和求解過(guò)程。多重網(wǎng)格建立過(guò)程描述如下，依次為網(wǎng)格粗化、插值矩陣計(jì)算、方程組矩陣及右端矢量轉(zhuǎn)換。其中，粗網(wǎng)格由插值矩陣Pk、限制矩陣PkT、粗網(wǎng)格層級(jí)的矩陣Mk+1和矢量bk+1組成。粗網(wǎng)格層級(jí)m的值一般不超過(guò)100，ω的值通常根據(jù)矩陣的譜半徑計(jì)算。

GAMG setup

M0=A, b0=b

f

or k=0, …, m

coarsen kthhierarchy of grid and get aggregates vector: aggk

construct tentative interpolation matrix Tkbased on bkand aggk

bk+1=TkTbk

Interpolation matrix Pk=(I-ωD-1Mk)Tk

Mk+1=PkTMkPk

end for

GAMG法的求解過(guò)程是一個(gè)遞歸過(guò)程，其算法描述如下。首先，方程組在第一層網(wǎng)格進(jìn)行幾次權(quán)重Jacobi光滑迭代，然后計(jì)算出殘差r并傳遞到下一層較粗的網(wǎng)格，在下一層網(wǎng)格遞歸求解x的修正值e，進(jìn)行到最后一層網(wǎng)格上時(shí)直接求解出結(jié)果，然后將修正值e插值到上一層網(wǎng)格并修正x，再進(jìn)行幾次權(quán)重Jacobi迭代。重復(fù)以上遞歸計(jì)算直到結(jié)果收斂。

GAMG Solve

GAMG_solver (Mk, Pk, xk, bk)

If k＜m

xk=weight_jacobi_iter (Mk, bk)

rk=bk-Mkxk

rk+1=PkTrk

GAMG_solver (Mk+1, Pk+1, ek+1, rk+1)

ek=Pkek+1

xk=xk+ek

xk=weight_jacobi_iter (Mk, bk)

else

solve Mkxk=bkdirectly

End

GAMG算法由若干矢量及矩陣線性運(yùn)算組成，流程復(fù)雜，單次計(jì)算耗時(shí)長(zhǎng)，但是在網(wǎng)格規(guī)模很大且物理方程難以收斂的情況下能有效加快方程的收斂速度。

3.1 GPU加速方法及程序結(jié)構(gòu)

采用的算法涉及了大量線性代數(shù)運(yùn)算，如向量加法、向量點(diǎn)積、對(duì)角矩陣-向量乘法等，這些線性代數(shù)運(yùn)算具有天然的并行性，使用CUDA編程語(yǔ)言將這些基本運(yùn)算移植到GPU，從而加速求解過(guò)程。

向量加減法、乘法、向量-標(biāo)量運(yùn)算的GPU并行計(jì)算方法相似。以向量加法為例，其并行計(jì)算過(guò)程為若干GPU的線程同時(shí)計(jì)算結(jié)果向量的某個(gè)元素，如圖 2所示。

圖2 GPU并行計(jì)算向量加法Fig.2 Parallel vector addition in GPU.

向量的點(diǎn)積運(yùn)算由向量元素的乘法及元素間的規(guī)約求和操作完成。采用線程順序?qū)ぶ返姆绞竭M(jìn)行并行規(guī)約操作，如圖3所示。對(duì)角矩陣-向量乘法算法的過(guò)程等效于每個(gè)對(duì)角線的元素與向量相乘并求和的過(guò)程。

求解程序分為三個(gè)層級(jí)（圖4），底層為GPU線性代數(shù)運(yùn)算模塊，中間層級(jí)為基于底層模塊的GPU線性方程求解模塊，頂層層級(jí)為流場(chǎng)及溫度場(chǎng)求解模塊。兩個(gè)物理求解模塊執(zhí)行完整的外迭代計(jì)算流程，并在對(duì)方程離散后都調(diào)用下層的GPU線性方程組求解模塊進(jìn)行計(jì)算。

圖3 GPU并行規(guī)約求和Fig.3 Parallel summation reduction in GPU.

圖4 GATH程序結(jié)構(gòu)Fig.4 GATH code structure.

4 結(jié)果及討論

4.1 程序校核

根據(jù)中國(guó)科學(xué)院上海應(yīng)用物理研究所的PB-FHR實(shí)驗(yàn)堆設(shè)計(jì)參數(shù)建立了簡(jiǎn)化的PB-FHR堆芯幾何模型，使用GATH程序?qū)ζ溥M(jìn)行熱工水力分析，并與計(jì)算流體動(dòng)力學(xué)(Computational Fluid Dynamics, CFD)軟件的結(jié)果進(jìn)行了對(duì)比。圖5為簡(jiǎn)化的堆芯模型的縱截面示意圖，堆芯半徑0.675 m，高1.91 m。下部區(qū)域?yàn)榍虼?，球床區(qū)域高1.61 m，球床出口半徑15 cm，球床出口處的斜面水平夾角為30°，上部區(qū)域?yàn)閹舾扇埯}通道的上反射層，堆芯燃料球直徑為0.06 m，堆芯球床及上反射層均視為多孔介質(zhì)。參考伯克利大學(xué)的實(shí)驗(yàn)結(jié)果[25]，設(shè)PB-FHR球床孔隙為0.4，由于缺少上反射層的孔隙率測(cè)量值，假設(shè)上反射層孔隙率和球床孔隙率近似相同。

圖5 PB-FHR堆芯幾何模型Fig.5 Reactor core model of PB-FHR.

堆芯內(nèi)冷卻劑為液態(tài)FLiBe熔鹽，由于FLiBe的密度略大于燃料球的密度，燃料球在堆芯處于漂浮狀態(tài)。堆芯底部為熔鹽入口，入口溫度873.15 K，入口流速0.052 m·s-1，堆芯頂部為熔鹽出口，出口采用流出邊界條件，模型外側(cè)壁面設(shè)為絕熱。堆芯熱功率為10 MW，堆芯內(nèi)的功率密度分布數(shù)據(jù)由中子輸運(yùn)程序提供。

圖6 堆芯固體(a)和冷卻劑(b)溫度分布Fig.6 Temperature distribution of solid phase (a) and coolant (b) in reactor core.

圖6(a)為GATH程序計(jì)算出的堆芯內(nèi)的固相溫度分布，包括球床及上反射層溫度。堆芯最高溫度位于中心處、高度1 m附近，堆芯入口附近區(qū)域及上反射層區(qū)域溫度較低，由于外側(cè)壁面功率有所上升，導(dǎo)致壁面溫度也較高。圖6(b)為GATH程序計(jì)算出的堆芯冷卻劑溫度分布。冷卻劑溫度沿著堆芯縱向一直上升，并于堆芯出口達(dá)到最高。壁面附近的冷卻劑也略有上升。球床多孔介質(zhì)強(qiáng)烈的湍流熱彌散效應(yīng)使得冷卻劑的徑向溫度分布比堆芯固體溫度分布均勻。

圖7-10對(duì)比了GATH程序與CFX軟件的計(jì)算結(jié)果，其中堆芯徑向的流場(chǎng)及溫度分布來(lái)自堆芯中高度為1 m處的網(wǎng)格點(diǎn)。GATH計(jì)算出的流場(chǎng)及溫度分布與CFX的結(jié)果符合得很好，所有位置的溫度差異均在2 K以?xún)?nèi)。表1顯示，GATH計(jì)算結(jié)果的相對(duì)誤差非常小，其中固體溫度的相對(duì)誤差小于冷卻劑溫度的誤差。

表1 GATH程序與CFX結(jié)果對(duì)比Table 1 Comparison between GATH and CFX.

圖7 堆芯固體(a)和冷卻劑(b)徑向溫度分布Fig.7 Radial temperature distribution of solid phase (a) and coolant (b) in reactor core.

圖8 堆芯流速(a)和壓力(b)徑向分布Fig.8 Radial velocity (a) and pressure (b) distribution of coolant in reactor core.

圖9 堆芯固體(a)和冷卻劑(b)軸向溫度分布Fig.9 Axial temperature distribution of solid phase (a) and coolant phase (b) in reactor core.

圖10 堆芯流速(a)和壓力(b)軸向分布Fig.10 Axial velocity (a) and pressure (b) distribution of coolant in reactor core.

4.2 GPU加速性能分析

針對(duì)每一個(gè)迭代求解算法實(shí)現(xiàn)了CPU串行求解器和相應(yīng)的GPU并行求解器，各個(gè)版本的求解器均采用§4.1的模型進(jìn)行了校核。CPU版本和GPU版本的程序均在64位Centos 6.7系統(tǒng)下由CUDA8.0的nvcc編譯器編譯。

為了分析GPU的加速性能，分別使用GPU并行求解器和CPU串行求解器進(jìn)行PB-FHR堆芯熱工水力計(jì)算，并統(tǒng)計(jì)不同求解器的耗時(shí)情況。用于測(cè)試的GPU型號(hào)為NVIDIA Tesla M2090，包含512個(gè)計(jì)算核心，運(yùn)算頻率1.3 MHz，雙精度浮點(diǎn)運(yùn)算速度665 Gflops（每秒6650億次浮點(diǎn)運(yùn)算）；CPU型號(hào)為Intel Xeon X5675，運(yùn)算頻率3.06 GHz，單線程雙精度浮點(diǎn)數(shù)峰值運(yùn)算速度達(dá)12.24 Gflops。

將堆芯模型劃分為2457600個(gè)網(wǎng)格，堆芯模型在徑向、軸向和角度方向的劃分?jǐn)?shù)目分別為480、160和32。表2顯示了POLYN和GAMG預(yù)處理的迭代算法耗時(shí)情況。其中，壓力修正方程和冷卻劑溫度方程的計(jì)算迭代次數(shù)較多，說(shuō)明壓力修正方程和冷卻劑溫度方程與其他物理方程相比不易收斂。POLYN預(yù)處理算法的加速比率最高可達(dá)20倍。GAMG預(yù)處理算法求解壓力修正方程和冷卻劑溫度方程的計(jì)算迭代次數(shù)明顯少于POLYN預(yù)處理算法，對(duì)應(yīng)的求解時(shí)間也大大減少。但是GAMG預(yù)處理算法求解其他易于收斂的物理方程時(shí)比POLYN預(yù)處理算法慢。GAMG預(yù)處理算法的加速比率最高可達(dá)10倍。為了使整體求解速度最快，程序應(yīng)使用GAMG預(yù)處理算法求解壓力修正方程及冷卻劑溫度方程，而使用POLYN預(yù)處理算法求解其他物理方程。

表2 POLYN預(yù)處理算法耗時(shí)Table 2 Elapsed time of POLYN preconditioning method.

根據(jù)上文的分析，使用GAMG預(yù)處理算法求解壓力修正方程及冷卻劑溫度方程，使用POLYN預(yù)處理算法求解其他物理方程，并統(tǒng)計(jì)了不同網(wǎng)格數(shù)量時(shí)程序整體的加速比率，統(tǒng)計(jì)的時(shí)間包含了方程系數(shù)的計(jì)算、數(shù)據(jù)拷貝等所有操作的耗時(shí)。結(jié)果如圖11所示，加速比率最高達(dá)8.39倍。在網(wǎng)格數(shù)量較小時(shí)，可并行任務(wù)較少，數(shù)據(jù)讀寫(xiě)的時(shí)間占據(jù)較大比例，因此加速比率較低；隨著網(wǎng)格數(shù)量的增大，可并行任務(wù)增多，加速比率開(kāi)始增大；當(dāng)網(wǎng)格數(shù)量到達(dá)一定規(guī)模后，GPU的計(jì)算資源飽和，加速比率趨于不變。

圖11 不同網(wǎng)格數(shù)量時(shí)GPU加速比率Fig.11 Speedup ratio vs. mesh quantity.

為進(jìn)一步證明GPU加速的有效性，將GATH程序與CFX的穩(wěn)態(tài)熱工水力計(jì)算耗時(shí)進(jìn)行了初步的對(duì)比。GATH與CFX模擬的PB-FHR堆芯模型采用相同的網(wǎng)格數(shù)量(449455)，兩個(gè)程序采用相同的收斂準(zhǔn)則（能量殘差小于10-6，其余物理量殘差小于10-3），CFX采用串行計(jì)算方式，GATH采用GPU并行計(jì)算方式。GATH進(jìn)行35次外迭代后結(jié)果收斂，總耗時(shí)46.25 s；CFX進(jìn)行30次外迭代后結(jié)果收斂，總耗時(shí)139.88 s。因此，GPU加速比率為3.02倍。

5 結(jié)語(yǔ)

本文基于計(jì)算流體力學(xué)方法，研發(fā)了GPU加速的PB-FHR堆芯熱工水力程序，采用非熱平衡多孔介質(zhì)模型建立堆芯固體結(jié)構(gòu)及冷卻劑的物理模型，采用SIMPLEC算法求解流場(chǎng)及壓力場(chǎng)，引入多孔介質(zhì)的湍流模型模擬球床流體的熱彌散效應(yīng)，實(shí)現(xiàn)了可在GPU上高速運(yùn)算的CG、BICGSTAB求解算法與POLYN、GAMG預(yù)處理算法。采用GATH程序?qū)B-FHR進(jìn)行了堆芯熱工水力分析，并與商用計(jì)算流體力學(xué)軟件CFX的計(jì)算結(jié)果進(jìn)行了對(duì)比，同時(shí)詳細(xì)分析了GPU的加速性能。

結(jié)果表明：1) GATH程序與CFX的結(jié)果差異很小；2) POLYN預(yù)處理的CG及BICGSTAB算法可達(dá)20倍左右的GPU加速比率，GAMG預(yù)處理的CG及BICGSTAB算法可達(dá)10倍左右的GPU加速比率；3) 同一算法求解不同物理方程具有不同的計(jì)算效率，采用GAMG預(yù)處理算法求解壓力修正方程及冷卻劑溫度方程而采用POLYN預(yù)處理算法求解其余物理方程可使得GATH程序整體計(jì)算速度最高；4) GATH程序整體的加速性能隨著網(wǎng)格數(shù)量而增大，最后趨于不變，最高可達(dá)8.39倍的GPU加速比率。通過(guò)以上結(jié)果證明了GATH程序中物理模型和求解算法的正確性，并且所采用的GPU加速方法能有效提升堆芯熱工水力分析的計(jì)算效率，為將來(lái)PB-FHR堆芯設(shè)計(jì)及分析提供了有效手段。

1 Forsberg C W, Peterson P F, Pickard P S. Molten-salt-cooled advanced high-temperature reactor for production of hydrogen and electricity[J]. Nuclear Technology, 2003, 144(3): 289-302. DOI: 10.13182/ NT03-1.

2 De Zwaan S, Boer B, Lathouwers D, et al. Static design of a liquid-salt-cooled pebble bed reactor (LSPBR)[J]. Annals of Nuclear Energy, 2007, 34(1): 83-92. DOI: 10.1016/j.anucene.2006.11.008.

3 Andreades A T C, Choi J K, Chong A Y K, et al. Technical description of the ‘Mark 1’ pebble-bed fluoride-salt-cooled high-temperature reactor (PB-FHR) power plant[D]. Berkeley: University of California, 2014.

4 Serp J, Allibert M, Bene? O, et al. The molten salt reactor (MSR) in generation IV: overview and perspectives[J]. Progress in Nuclear Energy, 2014, 77: 308-319. DOI: 10.1016/j.pnucene.2014.02.014.

5 Griveau A. Modeling and transient analysis for the pebble bed advanced high temperature reactor (PB-AHTR)[D]. Berkeley: Department of Nuclear Engineering, University of California, 2007.

6 牛強(qiáng), 宋士雄, 魏泉, 等. 熔鹽冷卻球床堆熱通道熱工水力特性數(shù)值分析[J]. 核技術(shù), 2014, 37(7): 070602. DOI: 10.11889/j.0253-3219.2014.hjs.37.070602. NIU Qiang, SONG Shixiong, WEI Quan, et al. Thermal-hydraulics numerical analyses of pebble bed advanced high temperature reactor hot channel[J]. Nuclear Techniques, 2014, 37(7): 070602. DOI: 10.11889/j.0253-3219.2014.hjs.37.070602.

7 Wang C, Xiao Y, Zhou J, et al. Computational fluid dynamics analysis of a fluoride salt-cooled pebble-bed test reactor[J]. Nuclearence & Engineering, 2014, 178(1): 86-102. DOI: 10.13182/NSE13-60.

8 Ge J, Wang C, Xiao Y, et al. Thermal-hydraulic analysis of a fluoride-salt-cooled pebble-bed reactor with CFD methodology[J]. Progress in Nuclear Energy, 2016, 91: 83-96. DOI: 10.1016/j.pnucene.2016.01.011.

9 Kirk D, Hwu W. Programming massively parallel processors[M]. Boston: Morgan Kaufmann, 2010.

10 Bailey P, Myre J, Walsh S D C, et al. Accelerating lattice boltzmann fluid flow simulations using graphics processors[C]. Proceedings of the International Conference on Parallel Processing, International Conference on Cultural Policy, Vienna, Austria, 2009: 550-557.

11 Boyd W R, Smith K, Forget B, et al. A massively parallel method of characteristic neutral particle transport code for GPUs[M]. La Grange Park, United States: American Nuclear Society, 2013.

12 De Lemos M J, Pedras M H. Recent mathematical models for turbulent flow in saturated rigid porous media[J]. Journal of Fluids Engineering, 2001, 123(4): 935-940. DOI: 10.1115/1.1413243.

13 陶文銓. 數(shù)值傳熱學(xué)[M]. 西安: 西安交通大學(xué)出版社, 2003. TAO Wenquan. Numerical heat transfer[M]. Xi’an: Xi’an Jiaotong University Press, 2003.

14 Nakayama A, Kuwahara F. A macroscopic turbulence model for flow in a porous medium[J]. Journal of Fluids Engineering, 1999, 121(2): 427-433. DOI: 10.1115/ 1.2822227.

15 Ergun S. Fluid flow through packed columns[J]. Chemical Engineering Progress, 1952, 48(2): 89-94.

16 Quintard M, Kaviany M, Whitaker S. Two-medium treatment of heat transfer in porous media: numerical results for effective properties[J]. Advances in Water Resources, 1997, 20(2): 77-94. DOI: 10.1016/ S0309-1708(96)00024-3.

17 Kuwahara F, Nakayama A. Numerical modelling of non-Darcy convective flow in a porous medium[C]. Heat Transfer Conference, Kyongju, Korea, 1998: 411-416.

18 Wakao N, Kaguei S, Funazkri T. Effect of fluid dispersion coefficients on particle-to-fluid heat transfer coefficients in packed beds: correlation of nusselt numbers[J]. Chemical Engineering Science, 1979, 34(3): 325-336. DOI: 10.1016/0009-2509(79)85064-2.

19 Guo B, Yu A, Wright B, et al. Simulation of turbulent flow in a packed bed[J]. Chemical Engineering & Technology, 2006, 29(5): 596-603. DOI: 10.1002/ceat. 200500292.

20 Kuzavkov N. Heat transport and afterheat removal for gas cooled reactors under accident conditions[D]. Vienna: International Atomic Energy Agency, 2000.

21 Hestenes M R, Stiefel E. Methods of conjugate gradients for solving linear systems[J]. Journal of Research of the National Bureau of Standards, 1952, 49(6): 409-436. DOI: 10.6028/jres.049.044.

22 Vorst H A V D. Bi-CGSTAB: a fast and smoothly converging variant of Bi-CG for the solution of nonsymmetric linear systems, SIAM J[J]. Siam Journal on Scientific & Statistical Computing, 1992, 13(2): 631-644. DOI: 10.1137/0913035.

23 Zhang J F, Zhang L. Efficient CUDA polynomial preconditioned conjugate gradient solver for finite element computation of elasticity problems[J]. Mathematical Problems in Engineering, 2013, (6): 1-12. DOI: 10.1155/2013/398438.

24 Stüben K. A review of algebraic multigrid[J]. Journal of Computational & Applied Mathematics, 2001, 128(1-2): 281-309. DOI: 10.1016/S0377-0427(00)00516-1.

25 Bardet P, An J Y, Franklin J T, et al. The pebble recirculation experiment (PREX) for the AHTR[C]. Proceedings of the Advanced Nuclear Fuel Cycles and Systems (GLOBAL 2007), La Grange Park: American Nuclear Society, 2007: 845-851.

Development of a GPU-accelerated thermal hydraulic code for pebble-bed fluoride-salt cooled high temperature reactor core

E Yanzhi1,2ZOU Yang1GUO Wei1PENG Yu1,2XU Hongjie1
1(Shanghai Institute of Applied Physics, Chinese Academy of Sciences, Jiading Campus, Shanghai 201800, China) 2(University of Chinese Academy of Sciences, Beijing 100049, China)

Background： Pebble bed fluoride-salt cooled high temperature reactor (PB-FHR) is a kind of Generation IV reactor. Three-dimensional thermal hydraulic code of reactor core is expected to simulate cases with complex spatial effect for PB-FHR but takes very heavy time-consuming. Graphics processing unit (GPU) contains numerous computing cores, thus can be used to efficiently accelerate numerical calculation if applied properly. Purpose： This study aims to develop a GPU-accelerated thermal-hydraulic code (GATH) for PB-FHR core. Methods： Thermal non-equilibrium porous media theory is adopted to build the reactor core physical model. Efficient iterative algorithms are researched and implemented on GPU. A PB-FHR core model is built for thermal hydraulic analysis with GATH. Simulation results are compared with ANSYS CFX software to verify GATH code and the GPU acceleration performance is analyzed. Results： The results between GATH and CFX are in good agreement. The speedup ratio of GATH can reach 8.39 times. Conclusion： The physical model and calculation method adopted inGATH code are right. The GPU accelerating methods proposed in this paper can efficiently accelerate thermal hydraulic simulation.

E Yanzhi, male, born in 1989, graduated from Tsinghua University in 2012, doctoral student, focusing on reactor thermal-hydraulics

GUO Wei, E-mail: guowei@sinap.ac.cn

date: 2017-03-07, accepted date: 2017-04-20

PB-FHR, GPU-accelerated, Thermal hydraulic, Thermal non-equilibrium porous media

TL333

10.11889/j.0253-3219.2017.hjs.40.070603

中國(guó)科學(xué)院戰(zhàn)略性先導(dǎo)科技專(zhuān)項(xiàng)(No.XDA02001002)、中國(guó)科學(xué)院前沿科學(xué)重點(diǎn)研究項(xiàng)目(No.QYZDY-SSW-JSC016)資助

鄂彥志，男，1989年出生，2012年畢業(yè)于清華大學(xué)，現(xiàn)為博士研究生，研究領(lǐng)域?yàn)榉磻?yīng)堆熱工水力

郭威，E-mail: guowei@sinap.ac.cn

2017-03-07，

2017-04-20

Supported by Strategic Pilot Science and Technology Project of Chinese Academy of Sciences (No.XDA02001002), the Frontier Science Key Program

of Chinese Academy of Sciences (No.QYZDY-SSW-JSC016)