明平洲 李治剛 潘俊杰 蘆韡

摘 要

CHF（Critical Heat Flux）一體化開發(fā)系統(tǒng)根據(jù)臨界熱流密度的實(shí)驗(yàn)數(shù)據(jù)與子通道計(jì)算進(jìn)行相互對(duì)比和統(tǒng)計(jì)學(xué)分析，最終對(duì)CHF經(jīng)驗(yàn)關(guān)系式進(jìn)行修正和優(yōu)化。由于CHF一體化開發(fā)系統(tǒng)的運(yùn)行效率偏低，且反復(fù)調(diào)用子通道計(jì)算模塊，所以引入多線程并行編程技術(shù)來加速計(jì)算。通過OpenMP和POSIX線程兩種編程方案取得不同的加速效果，區(qū)別在于后者編程難度更大。19組實(shí)驗(yàn)數(shù)據(jù)問題的CHF一體化開發(fā)系統(tǒng)運(yùn)行效率被提升5倍左右，接近加速的理論上限。

關(guān)鍵詞

CHF經(jīng)驗(yàn)關(guān)系式;子通道計(jì)算;并行計(jì)算;多線程

中圖分類號(hào)： TM75 ? ? ? ? ? ? ? ? ?文獻(xiàn)標(biāo)識(shí)碼： A

DOI：10.19694/j.cnki.issn2095-2457.2020.15.015

堆芯的熱工水力設(shè)計(jì)準(zhǔn)則之一是在正常、運(yùn)行瞬態(tài)及中等頻率事故工況下，堆芯不發(fā)生偏離泡核沸騰。當(dāng)前計(jì)算偏離泡核沸騰比率（DNBR，Departure From Nucleate Boiling Ratio）的主流是采用子通道程序。目前采用子通道模型進(jìn)行堆芯熱工水力分析的程序很多，世界上主要的核電設(shè)計(jì)公司均擁有各自的子通道程序，如法國(guó)阿?，敼鹃_發(fā)的FLICA程序[1]，美國(guó)西北太平洋國(guó)家實(shí)驗(yàn)室開發(fā)的COBRA程序[2]是典型代表，以上程序都已在工程上得到應(yīng)用。中國(guó)核動(dòng)力研究設(shè)計(jì)院設(shè)計(jì)所近期也研制了具備全堆芯計(jì)算能力的子通道程序CORTH。

對(duì)于堆芯的穩(wěn)態(tài)工況和預(yù)計(jì)事故工況都要分別計(jì)算出最小的DNBR值。子通道計(jì)算中要準(zhǔn)確地得到DNBR值，首先必須準(zhǔn)確地計(jì)算出臨界熱流密度值CHF，記為qCHF。預(yù)測(cè)堆芯的qCHF較為復(fù)雜，堆芯軸向和徑向熱流密度的分布、棒的間隙以及流體、定位格架等多種因素都會(huì)對(duì)臨界熱流密度產(chǎn)生一定影響。經(jīng)驗(yàn)關(guān)系式在預(yù)測(cè)堆芯qCHF方面使用較多，例如著名的W-2關(guān)系式、W-3關(guān)系式和WRB-1公式等，它們對(duì)應(yīng)于不同的工況和沸騰條件。查表方法近年來也是重要預(yù)測(cè)熱流密度的重要方案，它應(yīng)用范圍更廣，便于對(duì)特殊情況進(jìn)行修正[3]。經(jīng)驗(yàn)關(guān)系式涉及的大量實(shí)驗(yàn)數(shù)據(jù)和計(jì)算數(shù)據(jù)往往通過多個(gè)步驟，多種程序工具進(jìn)行處理，耗時(shí)較長(zhǎng)。CHF一體化開發(fā)系統(tǒng)是一種偏向于經(jīng)驗(yàn)關(guān)系式方法的軟件系統(tǒng)，它對(duì)堆芯的qCHF進(jìn)行預(yù)測(cè)，確定和優(yōu)化CHF的經(jīng)驗(yàn)關(guān)系式。系統(tǒng)聯(lián)合實(shí)驗(yàn)數(shù)據(jù)和子通道計(jì)算數(shù)據(jù)，縮短了數(shù)據(jù)處理周期。本文是使用多線程并行技術(shù)對(duì)研制的CHF一體化開發(fā)系統(tǒng)進(jìn)行并行優(yōu)化的研究，并歸納此類應(yīng)用模式下的并行優(yōu)化特點(diǎn)。

1 CHF一體化開發(fā)系統(tǒng)描述

1.1 計(jì)算原理

1.1.1 堆芯熱工水力學(xué)參數(shù)

壓水堆中燃料元件表面的最大熱流密度應(yīng)小于臨界熱流密度，DNBR實(shí)質(zhì)上是指利用合適的熱流量預(yù)測(cè)方法得到的冷卻劑通道中燃料元件表面某一點(diǎn)的臨界熱流量與該位置的實(shí)際熱流量比值[4]。

其中z為堆芯軸向的位置，子通道模型中徑向則被劃分為多個(gè)冷卻劑通道。DNBR的值沿軸向是變化的，其最小值（DNBR）min稱為該子通道的最小DNBR。在壓水堆熱工設(shè)計(jì)中，把整個(gè)堆芯（所有子通道）最小DNBR不小于某個(gè)值作為熱工設(shè)計(jì)準(zhǔn)則之一。

研制的CHF一體化開發(fā)系統(tǒng)構(gòu)造了類似于已公開發(fā)表內(nèi)容[5]的關(guān)系式。關(guān)系式所需的熱工水力本地參數(shù)可以由子通道程序計(jì)算得出?，F(xiàn)階段關(guān)系式適用于壓水堆的熱工水力分析。

1.1.2 非線性回歸

CHF經(jīng)驗(yàn)關(guān)系式的準(zhǔn)確性主要與關(guān)系式的各項(xiàng)系數(shù)相關(guān)。根據(jù)實(shí)驗(yàn)數(shù)據(jù)和子通道計(jì)算產(chǎn)生的中間結(jié)果，可以采用非線性回歸方法進(jìn)行系數(shù)的確定、修正或者優(yōu)化。

常用的非線性回歸方法有：共軛梯度法、Quasi-Newton法等。Quasi-Newton法是求解非線性優(yōu)化問題最有效的方法之一，本文選用Quasi-Newton法擬合CHF經(jīng)驗(yàn)關(guān)系式，該方法在系統(tǒng)使用的第三方工具Scilab[6]內(nèi)的調(diào)用格式為：

[x，err]=datafit（g，z，x0，qn）（3）

其中，g表示擬合目標(biāo)函數(shù)，z為存放測(cè)量數(shù)據(jù)的矩陣，err為誤差返回值，x0為系數(shù)初始值，x為擬合所得系數(shù)。qn則制定使用的數(shù)據(jù)擬合方法，這里選定為Quasi-Newton法。

1.2 系統(tǒng)的程序結(jié)構(gòu)

CHF一體化開發(fā)系統(tǒng)目前采用結(jié)構(gòu)化程序設(shè)計(jì)，采用C++編程語(yǔ)言書寫控制流程和數(shù)據(jù)存儲(chǔ)，同時(shí)預(yù)留接口，對(duì)多個(gè)第三方程序進(jìn)行集成。整個(gè)程序具備模塊化編程的特點(diǎn)，各個(gè)重要功能可以進(jìn)行合理替換。

上圖是整個(gè)CHF一體化開發(fā)系統(tǒng)的控制流程框圖，其中子通道程序計(jì)算和統(tǒng)計(jì)學(xué)檢驗(yàn)依賴于第三方程序。圖中子通道計(jì)算模塊需要對(duì)多組數(shù)據(jù)調(diào)用外部的子通道程序進(jìn)行計(jì)算處理，也是本文中需要進(jìn)行并行編程應(yīng)用的計(jì)算區(qū)域。

1.2.1 第三方程序調(diào)用

子通道程序的建模和計(jì)算采用中國(guó)核動(dòng)力研究設(shè)計(jì)院自主研制的CORTH子通道程序。為了與CHF一體化開發(fā)系統(tǒng)較好的耦合，對(duì)CORTH程序進(jìn)行代碼級(jí)別的集成;CHF一體化開發(fā)系統(tǒng)調(diào)用解釋型腳本來耦合SCILAB的接口，實(shí)現(xiàn)非線性回歸等功能;調(diào)用第三方畫圖軟件Gnuplot實(shí)現(xiàn)對(duì)關(guān)鍵參數(shù)分布、數(shù)據(jù)趨勢(shì)進(jìn)行圖像化輸出

1.2.2 模塊功能

CHF實(shí)驗(yàn)數(shù)據(jù)預(yù)處理：對(duì)CHF實(shí)驗(yàn)獲得的大量數(shù)據(jù)點(diǎn)進(jìn)行預(yù)處理，如去除實(shí)驗(yàn)壞點(diǎn)等，形成開發(fā)數(shù)據(jù)庫(kù)和驗(yàn)證數(shù)據(jù)庫(kù)（對(duì)比驗(yàn)證）。

實(shí)驗(yàn)數(shù)據(jù)子通道建模和計(jì)算：批量調(diào)用子通道程序，對(duì)CHF實(shí)驗(yàn)數(shù)據(jù)進(jìn)行計(jì)算和預(yù)測(cè)。

CHF經(jīng)驗(yàn)關(guān)系式系數(shù)求解：根據(jù)CHF開發(fā)數(shù)據(jù)庫(kù)對(duì)應(yīng)的子通道計(jì)算結(jié)果確定并優(yōu)化CHF關(guān)系式的形式，并求解得到CHF關(guān)系式的各項(xiàng)系數(shù)。

建立M/P數(shù)據(jù)庫(kù)：將前一階段計(jì)算得到的CHF關(guān)系式嵌入至子通道程序中，再次對(duì)CHF開發(fā)數(shù)據(jù)庫(kù)調(diào)用子通道程序計(jì)算，提取臨界位置的DNBR等參數(shù)，形成M/P數(shù)據(jù)庫(kù)。

確定DNBR限值：對(duì)M/P數(shù)據(jù)庫(kù)進(jìn)行統(tǒng)計(jì)學(xué)分析和檢驗(yàn)，計(jì)算DNBR限值。

CHF經(jīng)驗(yàn)關(guān)系式的驗(yàn)證和評(píng)價(jià)：采用CHF驗(yàn)證數(shù)據(jù)庫(kù)，調(diào)用已嵌入新CHF關(guān)系式的子通道程序進(jìn)行建模和計(jì)算，獲得驗(yàn)證數(shù)據(jù)庫(kù)的M/P數(shù)據(jù)庫(kù)，結(jié)合開發(fā)數(shù)據(jù)庫(kù)的M/P數(shù)據(jù)庫(kù)，評(píng)估新CHF關(guān)系式的預(yù)測(cè)率并對(duì)熱工參數(shù)的變化趨勢(shì)等進(jìn)行評(píng)價(jià)。

2 并行編程方案描述

根據(jù)程序結(jié)構(gòu)的描述，多組實(shí)驗(yàn)數(shù)據(jù)涉及的子通道程序建模和計(jì)算具備并行度。它們調(diào)用第三方程序CORTH完成子通道的計(jì)算相互獨(dú)立，可以使用多線程并行編程技術(shù)來對(duì)計(jì)算加速。

CHF一體化開發(fā)系統(tǒng)占用的內(nèi)存在2G以下，運(yùn)行于工作站服務(wù)器，多線程并行編程屬于共享式內(nèi)存并行策略，可以按照如下兩類算法形式來實(shí)施：

2.1 OpenMP并行

基于編譯器提供的OpenMP并行編程[7]方法，對(duì)CHF一體化開發(fā)系統(tǒng)中的全局變量進(jìn)行轉(zhuǎn)換，降低數(shù)據(jù)相關(guān)性，然后添加OpenMP指導(dǎo)語(yǔ)句。

OpenMP并行編程使用的是編譯器指導(dǎo)語(yǔ)句，開發(fā)人員一般不需要進(jìn)行顯式確定調(diào)度策略，由于實(shí)驗(yàn)數(shù)據(jù)CHFExpData與數(shù)據(jù)組數(shù)i相關(guān)，不均勻性會(huì)帶來一定負(fù)載均衡問題。

2.2 POSIX線程并行

類UNIX計(jì)算服務(wù)器提供了POSIX標(biāo)準(zhǔn)線程庫(kù)[8]，相比于OpenMP并行編程更貼近操作系統(tǒng)底層。CHF一體化開發(fā)系統(tǒng)在初始階段生成多個(gè)線程，構(gòu)成池式控制結(jié)構(gòu)，通過互斥開關(guān)來使用這些線程完成并行計(jì)算。

POSIX多線程并行充分利用了線程細(xì)粒度運(yùn)行任務(wù)的特點(diǎn)，以單次調(diào)用子通道計(jì)算為基本任務(wù)，避免了負(fù)載均衡等問題，但線程池方案存在資源互斥，編程難度大，且輸出結(jié)果的次序也會(huì)影響后續(xù)的統(tǒng)計(jì)分析。

3 實(shí)驗(yàn)結(jié)果與分析

CHF一體化開發(fā)系統(tǒng)及其并行優(yōu)化的初始版本目前已經(jīng)研制完畢，這里選取HP Z800計(jì)算服務(wù)器作為運(yùn)行環(huán)境。

整個(gè)實(shí)驗(yàn)涉及的CHF實(shí)驗(yàn)數(shù)據(jù)分為開發(fā)實(shí)驗(yàn)數(shù)據(jù)和驗(yàn)證實(shí)驗(yàn)數(shù)據(jù)兩類，其中開發(fā)實(shí)驗(yàn)數(shù)據(jù)共19組，下表是每組實(shí)驗(yàn)數(shù)據(jù)對(duì)應(yīng)的CHF實(shí)驗(yàn)樣本點(diǎn)個(gè)數(shù)。

從輸入來看，19組實(shí)驗(yàn)數(shù)據(jù)進(jìn)行建模和計(jì)算存在著不均勻性，多線程并行的負(fù)載均衡需要進(jìn)行考慮。按照整個(gè)CHF一體化開發(fā)系統(tǒng)的運(yùn)行過程，完成整個(gè)計(jì)算需要對(duì)這19組實(shí)驗(yàn)數(shù)據(jù)前后迭代四次。實(shí)驗(yàn)結(jié)果先統(tǒng)計(jì)單次迭代運(yùn)行效率，即多線程并行一次處理19組實(shí)驗(yàn)數(shù)據(jù)的性能。整個(gè)CHF一體化開發(fā)系統(tǒng)的運(yùn)行時(shí)間（迭代四次）對(duì)比將在最后給出。

3.1 OpenMP并行編程

考慮到共有19組實(shí)驗(yàn)數(shù)據(jù)，而OpenMP并行僅僅針對(duì)外層循環(huán)實(shí)驗(yàn)組數(shù)進(jìn)行，所以這里選取4個(gè)，8個(gè)和10個(gè)線程三種情況進(jìn)行并行實(shí)驗(yàn)，觀察OpenMP并行編程情況下的收益。

其中4個(gè)線程對(duì)應(yīng)并行總共調(diào)用CORTH子通道次數(shù)為五次（=5），8個(gè)線程對(duì)應(yīng)并行調(diào)用三次，10個(gè)線程對(duì)應(yīng)并行調(diào)用兩次。在OpenMP并行編程的情況下，計(jì)算時(shí)間的波動(dòng)較為明顯，這是由于CHF開發(fā)實(shí)驗(yàn)數(shù)據(jù)的不均勻性造成，OpenMP的默認(rèn)調(diào)度方式無(wú)法做到負(fù)載均衡。

考慮10個(gè)線程運(yùn)行情況，在OpenMP的編譯指導(dǎo)語(yǔ)句中引入動(dòng)態(tài)調(diào)度。圖2是10個(gè)線程情況下OpenMP默認(rèn)調(diào)度方式與動(dòng)態(tài)調(diào)度方式的效率對(duì)比。由圖中所示，動(dòng)態(tài)調(diào)度方式比起默認(rèn)調(diào)度方式有了一定改進(jìn)，各個(gè)線程的負(fù)載相對(duì)均衡，但運(yùn)行時(shí)間比起默認(rèn)調(diào)度方式反倒有所惡化。這也從側(cè)面說明動(dòng)態(tài)調(diào)度方式在本計(jì)算例題情況下不是最優(yōu)選擇。

3.2 POSIX線程并行編程

利用POSIX標(biāo)準(zhǔn)線程庫(kù)開發(fā)，對(duì)于每組實(shí)驗(yàn)數(shù)據(jù)內(nèi)的單個(gè)實(shí)驗(yàn)樣本點(diǎn)均從線程池中獲取線程進(jìn)行運(yùn)行，當(dāng)線程池中沒有線程存在時(shí)則阻塞掛起，這種并行編程方式粒度較細(xì)。

為了兼顧通用性，OpenMP并行編程往往針對(duì)外層循環(huán)。CHF一體化開發(fā)系統(tǒng)中內(nèi)層循環(huán)的任務(wù)數(shù)與外層循環(huán)的序號(hào)相關(guān)，具有一定隨機(jī)性，所以直接采用POSIX線程的并行編程方式能夠帶來更多的收益。

對(duì)比表6和表7可以發(fā)現(xiàn)，OpenMP并行編程的編程難度較小，但從實(shí)驗(yàn)結(jié)果來看存在負(fù)載均衡問題;另一種POSIX線程的編程難度較大，但從實(shí)驗(yàn)結(jié)果來看并行收益更高。

此外，完整運(yùn)行CHF一體化開發(fā)系統(tǒng)對(duì)19組實(shí)驗(yàn)數(shù)據(jù)進(jìn)行處理。串行計(jì)算時(shí)除開調(diào)用子通道計(jì)算的內(nèi)容，每次完成19組實(shí)驗(yàn)數(shù)據(jù)的處理耗時(shí)約320秒左右，整個(gè)計(jì)算時(shí)間為1687秒左右，則可并行區(qū)域占總計(jì)算的百分比為80%左右。根據(jù)阿姆達(dá)爾定律，并行優(yōu)化之后帶來的加速比上限約為5倍左右。實(shí)際工程實(shí)現(xiàn)之后接近這一理論值，證明可并行區(qū)域具有天然并行的特性，但各個(gè)任務(wù)之間存在著負(fù)載不均衡的缺點(diǎn)，工程編程實(shí)現(xiàn)時(shí)通過POSIX線程的策略能避免此缺點(diǎn)。

4 結(jié)束語(yǔ)

本文對(duì)中國(guó)核動(dòng)力研究設(shè)計(jì)院的CHF一體化開發(fā)系統(tǒng)進(jìn)行了描述。引入了兩種并行編程方案，在調(diào)用子通道計(jì)算的過程中改善其計(jì)算效率。形成如下結(jié)論：

1）此問題域中，OpenMP并行編程實(shí)現(xiàn)簡(jiǎn)便，但存在負(fù)載均衡問題，POSIX線程池的方法實(shí)現(xiàn)困難，但負(fù)載較為平均。

2）CHF一體化開發(fā)系統(tǒng)的各個(gè)子任務(wù)與實(shí)驗(yàn)數(shù)據(jù)的規(guī)模相關(guān)，所使用的19組實(shí)驗(yàn)數(shù)據(jù)可并行區(qū)域帶來的加速比上限約為5倍左右，目前的并行編程工程實(shí)現(xiàn)接近這一理論值，改善了CHF一體化開發(fā)系統(tǒng)的運(yùn)行效率，為多次執(zhí)行CHF實(shí)驗(yàn)數(shù)據(jù)的分析和處理奠定了一定基礎(chǔ)。

后續(xù)該系統(tǒng)將應(yīng)用于CHF實(shí)驗(yàn)數(shù)據(jù)的評(píng)價(jià)和篩選，且數(shù)據(jù)之間的相關(guān)性以及CHF關(guān)系式優(yōu)化等內(nèi)容將被進(jìn)一步研究。

參考文獻(xiàn)

[1]Fillion P，Chanoine A，Dellacherie S，et al.FLICA-OVAP：A new platform for core thermal–hydraulic studies[J].Nuclear Engineering and Design，2011，241（11）：4348- 4358.

[2]Kelly J E，Loomis J N，Wolf L，et al.LWR core thermal-hydraulic analysis： assessment and comparison of the range of applicability of the codes COBRA IIIC/MIT and COBRA IV-I[R].OSTI-5493728，1978.

[3]于平安，朱瑞安，喻真烷.核反應(yīng)堆熱工分析，上海交通大學(xué)出版社[M].2003.

[4]蘇光輝，秋穗正，等.核動(dòng)力系統(tǒng)熱工水力計(jì)算方法，清華大學(xué)出版社[M].2013.

[5]FRAMATOME.FC-FRAMATOME Critical Heat Flux Correlation for AFA-2G and AFA-3G Fuel Assemblies-Topical Report （EPD/DC-454，Rev.A），December， 2000.

[6]胡寶鋼，趙星.科學(xué)計(jì)算自有軟件—SCILAB教程[M].清華大學(xué)出版社，2003.

[7]Salvini S.Unlocking the Power of OpenMP[C].Proc. Fifth European Workshop OpenMP （EWOMP03）， invited lecture，2003.

[8]March V，Zhang R，See S，et al.Survey on Parallel Programming Model[C].Network and Parallel Computing，2008.