明平洲 李治刚 潘俊杰 芦韡
摘 要
CHF（Critical Heat Flux）一体化开发系统根据临界热流密度的实验数据与子通道计算进行相互对比和统计学分析，最终对CHF经验关系式进行修正和优化。由于CHF一体化开发系统的运行效率偏低，且反复调用子通道计算模块，所以引入多线程并行编程技术来加速计算。通过OpenMP和POSIX线程两种编程方案取得不同的加速效果，区别在于后者编程难度更大。19组实验数据问题的CHF一体化开发系统运行效率被提升5倍左右，接近加速的理论上限。
关键词
CHF经验关系式;子通道计算;并行计算;多线程
中图分类号： TM75 文獻标识码： A
DOI：10.19694/j.cnki.issn2095-2457.2020.15.015
堆芯的热工水力设计准则之一是在正常、运行瞬态及中等频率事故工况下，堆芯不发生偏离泡核沸腾。当前计算偏离泡核沸腾比率（DNBR，Departure From Nucleate Boiling Ratio）的主流是采用子通道程序。目前采用子通道模型进行堆芯热工水力分析的程序很多，世界上主要的核电设计公司均拥有各自的子通道程序，如法国阿海玛公司开发的FLICA程序[1]，美国西北太平洋国家实验室开发的COBRA程序[2]是典型代表，以上程序都已在工程上得到应用。中国核动力研究设计院设计所近期也研制了具备全堆芯计算能力的子通道程序CORTH。
对于堆芯的稳态工况和预计事故工况都要分别计算出最小的DNBR值。子通道计算中要准确地得到DNBR值，首先必须准确地计算出临界热流密度值CHF，记为qCHF。预测堆芯的qCHF较为复杂，堆芯轴向和径向热流密度的分布、棒的间隙以及流体、定位格架等多种因素都会对临界热流密度产生一定影响。经验关系式在预测堆芯qCHF方面使用较多，例如著名的W-2关系式、W-3关系式和WRB-1公式等，它们对应于不同的工况和沸腾条件。查表方法近年来也是重要预测热流密度的重要方案，它应用范围更广，便于对特殊情况进行修正[3]。经验关系式涉及的大量实验数据和计算数据往往通过多个步骤，多种程序工具进行处理，耗时较长。CHF一体化开发系统是一种偏向于经验关系式方法的软件系统，它对堆芯的qCHF进行预测，确定和优化CHF的经验关系式。系统联合实验数据和子通道计算数据，缩短了数据处理周期。本文是使用多线程并行技术对研制的CHF一体化开发系统进行并行优化的研究，并归纳此类应用模式下的并行优化特点。
1 CHF一体化开发系统描述
1.1 计算原理
1.1.1 堆芯热工水力学参数
压水堆中燃料元件表面的最大热流密度应小于临界热流密度，DNBR实质上是指利用合适的热流量预测方法得到的冷却剂通道中燃料元件表面某一点的临界热流量与该位置的实际热流量比值[4]。
其中z为堆芯轴向的位置，子通道模型中径向则被划分为多个冷却剂通道。DNBR的值沿轴向是变化的，其最小值（DNBR）min称为该子通道的最小DNBR。在压水堆热工设计中，把整个堆芯（所有子通道）最小DNBR不小于某个值作为热工设计准则之一。
研制的CHF一体化开发系统构造了类似于已公开发表内容[5]的关系式。关系式所需的热工水力本地参数可以由子通道程序计算得出。现阶段关系式适用于压水堆的热工水力分析。
1.1.2 非线性回归
CHF经验关系式的准确性主要与关系式的各项系数相关。根据实验数据和子通道计算产生的中间结果，可以采用非线性回归方法进行系数的确定、修正或者优化。
常用的非线性回归方法有：共轭梯度法、Quasi-Newton法等。Quasi-Newton法是求解非线性优化问题最有效的方法之一，本文选用Quasi-Newton法拟合CHF经验关系式，该方法在系统使用的第三方工具Scilab[6]内的调用格式为：
[x，err]=datafit（g，z，x0，qn）（3）
其中，g表示拟合目标函数，z为存放测量数据的矩阵，err为误差返回值，x0为系数初始值，x为拟合所得系数。qn则制定使用的数据拟合方法，这里选定为Quasi-Newton法。
1.2 系统的程序结构
CHF一体化开发系统目前采用结构化程序设计，采用C++编程语言书写控制流程和数据存储，同时预留接口，对多个第三方程序进行集成。整个程序具备模块化编程的特点，各个重要功能可以进行合理替换。
上图是整个CHF一体化开发系统的控制流程框图，其中子通道程序计算和统计学检验依赖于第三方程序。图中子通道计算模块需要对多组数据调用外部的子通道程序进行计算处理，也是本文中需要进行并行编程应用的计算区域。
1.2.1 第三方程序调用
子通道程序的建模和计算采用中国核动力研究设计院自主研制的CORTH子通道程序。为了与CHF一体化开发系统较好的耦合，对CORTH程序进行代码级别的集成;CHF一体化开发系统调用解释型脚本来耦合SCILAB的接口，实现非线性回归等功能;调用第三方画图软件Gnuplot实现对关键参数分布、数据趋势进行图像化输出
1.2.2 模块功能
CHF实验数据预处理：对CHF实验获得的大量数据点进行预处理，如去除实验坏点等，形成开发数据库和验证数据库（对比验证）。
实验数据子通道建模和计算：批量调用子通道程序，对CHF实验数据进行计算和预测。
CHF经验关系式系数求解：根据CHF开发数据库对应的子通道计算结果确定并优化CHF关系式的形式，并求解得到CHF关系式的各项系数。
建立M/P数据库：将前一阶段计算得到的CHF关系式嵌入至子通道程序中，再次对CHF开发数据库调用子通道程序计算，提取临界位置的DNBR等参数，形成M/P数据库。
确定DNBR限值：对M/P数据库进行统计学分析和检验，计算DNBR限值。
CHF经验关系式的验证和评价：采用CHF验证数据库，调用已嵌入新CHF关系式的子通道程序进行建模和计算，获得验证数据库的M/P数据库，结合开发数据库的M/P数据库，评估新CHF关系式的预测率并对热工参数的变化趋势等进行评价。
2 并行编程方案描述
根据程序结构的描述，多组实验数据涉及的子通道程序建模和计算具备并行度。它们调用第三方程序CORTH完成子通道的计算相互独立，可以使用多线程并行编程技术来对计算加速。
CHF一体化开发系统占用的内存在2G以下，运行于工作站服务器，多线程并行编程属于共享式内存并行策略，可以按照如下两类算法形式来实施：
2.1 OpenMP并行
基于编译器提供的OpenMP并行编程[7]方法，对CHF一体化开发系统中的全局变量进行转换，降低数据相关性，然后添加OpenMP指导语句。
OpenMP并行编程使用的是编译器指导语句，开发人员一般不需要进行显式确定调度策略，由于实验数据CHFExpData与数据组数i相关，不均匀性会带来一定负载均衡问题。
2.2 POSIX线程并行
类UNIX计算服务器提供了POSIX标准线程库[8]，相比于OpenMP并行编程更贴近操作系统底层。CHF一体化开发系统在初始阶段生成多个线程，构成池式控制结构，通过互斥开关来使用这些线程完成并行计算。
POSIX多线程并行充分利用了线程细粒度运行任务的特点，以单次调用子通道计算为基本任务，避免了负载均衡等问题，但线程池方案存在资源互斥，编程难度大，且输出结果的次序也会影响后续的统计分析。
3 实验结果与分析
CHF一体化开发系统及其并行优化的初始版本目前已经研制完毕，这里选取HP Z800计算服务器作为运行环境。
整个实验涉及的CHF实验数据分为开发实验数据和验证实验数据两类，其中开发实验数据共19组，下表是每组实验数据对应的CHF实验样本点个数。
从输入来看，19组实验数据进行建模和计算存在着不均匀性，多线程并行的负载均衡需要进行考虑。按照整个CHF一体化开发系统的运行过程，完成整个计算需要对这19组实验数据前后迭代四次。实验结果先统计单次迭代运行效率，即多线程并行一次处理19组实验数据的性能。整个CHF一体化开发系统的运行时间（迭代四次）对比将在最后给出。
3.1 OpenMP并行编程
考虑到共有19组实验数据，而OpenMP并行仅仅针对外层循环实验组数进行，所以这里选取4个，8个和10个线程三种情况进行并行实验，观察OpenMP并行编程情况下的收益。
其中4个线程对应并行总共调用CORTH子通道次数为五次（=5），8个线程对应并行调用三次，10个线程对应并行调用两次。在OpenMP并行编程的情况下，计算时间的波动较为明显，这是由于CHF开发实验数据的不均匀性造成，OpenMP的默认调度方式无法做到负载均衡。
考虑10个线程运行情况，在OpenMP的编译指导语句中引入动态调度。图2是10个线程情况下OpenMP默认调度方式与动态调度方式的效率对比。由图中所示，动态调度方式比起默认调度方式有了一定改進，各个线程的负载相对均衡，但运行时间比起默认调度方式反倒有所恶化。这也从侧面说明动态调度方式在本计算例题情况下不是最优选择。
3.2 POSIX线程并行编程
利用POSIX标准线程库开发，对于每组实验数据内的单个实验样本点均从线程池中获取线程进行运行，当线程池中没有线程存在时则阻塞挂起，这种并行编程方式粒度较细。
为了兼顾通用性，OpenMP并行编程往往针对外层循环。CHF一体化开发系统中内层循环的任务数与外层循环的序号相关，具有一定随机性，所以直接采用POSIX线程的并行编程方式能够带来更多的收益。
对比表6和表7可以发现，OpenMP并行编程的编程难度较小，但从实验结果来看存在负载均衡问题;另一种POSIX线程的编程难度较大，但从实验结果来看并行收益更高。
此外，完整运行CHF一体化开发系统对19组实验数据进行处理。串行计算时除开调用子通道计算的内容，每次完成19组实验数据的处理耗时约320秒左右，整个计算时间为1687秒左右，则可并行区域占总计算的百分比为80%左右。根据阿姆达尔定律，并行优化之后带来的加速比上限约为5倍左右。实际工程实现之后接近这一理论值，证明可并行区域具有天然并行的特性，但各个任务之间存在着负载不均衡的缺点，工程编程实现时通过POSIX线程的策略能避免此缺点。
4 结束语
本文对中国核动力研究设计院的CHF一体化开发系统进行了描述。引入了两种并行编程方案，在调用子通道计算的过程中改善其计算效率。形成如下结论：
1）此问题域中，OpenMP并行编程实现简便，但存在负载均衡问题，POSIX线程池的方法实现困难，但负载较为平均。
2）CHF一体化开发系统的各个子任务与实验数据的规模相关，所使用的19组实验数据可并行区域带来的加速比上限约为5倍左右，目前的并行编程工程实现接近这一理论值，改善了CHF一体化开发系统的运行效率，为多次执行CHF实验数据的分析和处理奠定了一定基础。
后续该系统将应用于CHF实验数据的评价和筛选，且数据之间的相关性以及CHF关系式优化等内容将被进一步研究。
参考文献
[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]于平安，朱瑞安，喻真烷.核反应堆热工分析，上海交通大学出版社[M].2003.
[4]苏光辉，秋穗正，等.核动力系统热工水力计算方法，清华大学出版社[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]胡宝钢，赵星.科学计算自有软件—SCILAB教程[M].清华大学出版社，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.
