子通道程序对PSBT空泡分布实验计算的不确定性量化分析
Uncertainty Quantitative Analysis of Subchannel Code Calculation of PSBT Void Distribution Benchmark
通讯作者: 刘晓晶,男,研究员,博士生导师,电话(Tel.):021-34207121;E-mail:xiaojingliu@sjtu.edu.cn.
责任编辑: 孙伟
收稿日期: 2021-03-3
基金资助: |
|
Received: 2021-03-3
作者简介 About authors
张俊涛(1996-),男,河南省洛阳市人,硕士生,从事核反应堆热工水力研究.
为了评估子通道程序的准确性与可靠性,需要定量给出计算结果的不确定性.采用统计学上基于输入参数不确定性传递的方法进行不确定性分析,可以定量得到程序计算结果的不确定范围.在假设模型参数不确定性服从正态分布的基础上,采用统计学方法确定模型参数不确定性的分布以取代传统的专家判断.通过对压水堆子通道和棒束实验(PSBT)基准题空泡分布实验进行计算,分析子通道程序COBRA-IV 对实验结果的预测能力,同时得到满足容忍限的计算结果不确定性上下限.计算结果表明:评估得到的不确定带能较好地包络实验值;同时利用统计均值对模型进行标定后,可以得到比原模型更接近实验值的计算结果.
关键词:
In order to evaluate the accuracy and reliability of the subchannel code, it is necessary to quantitatively give the uncertainty of the calculation results. The uncertainty analysis is conducted by using the statistical method based on propagation of input uncertainties, and the uncertainty range of the subchannel code calculation results can be obtained quantitatively. Based on the assumption that the uncertainty of model parameters obeys normal distribution, the statistical method is used to determine the distribution of the uncertainty of model parameters to replace the traditional expert judgment. Through the calculation of the pressurized water reactor sub-channel and bundle tests (PSBT) benchmark, the ability of the subchannel code COBRA-IV to predict the experimental results is analyzed, and the uncertainty interval satisfying the tolerance limit of the calculation results is obtained. The results demonstrate that the experiment data is well enveloped by the obtained uncertainty bands and the model calibrated by the statistical mean value presents a good improvement of calculations.
Keywords:
本文引用格式
张俊涛, 刘晓晶, 张滕飞, 柴翔.
ZHANG Juntao, LIU Xiaojing, ZHANG Tengfei, CHAI Xiang.
1989年美国发布管理导则RG 1.157[1],允许在应急堆芯冷却系统分析中使用最佳估算程序.出于安全考虑,最佳估算程序必须与不确定性分析相结合,以判断计算与实际情况之间的差异.截至目前,国际上已有多种不同的不确定性分析方法,例如程序比例、适用性和不确定性分析方法(CSAU)[2]、德国核安全中心方法(GRS)[3]、基于精度外推的不确定性方法(UMAE)[4]、先进统计学处理的不确定性方法(ASTRUM)[5]等.这些不确定性分析方法大体可以分为两大类,一类是基于输入参数不确定性的传递,另一类是基于输出结果不确定性的外推[6].基于输入参数不确定性传递的统计性方法的研究和应用相对更成熟和广泛,如在由经济合作与发展组织(OECD)发起的不确定性研究项目“最佳估算方法-不确定性和敏感性分析”阶段V中,14个参与者中有12个采用输入不确定性传递的方法[7].
基于输入参数不确定性传递的方法需要知道输入参数的不确定性.然而,大多数热工水力程序都没有提供关于输入参数不确定性的足够信息,尤其是对于某些无法通过实验直接测量的输入模型参数,通常通过专家判断来进行量化.目前一些取代专家判断的评估方法被陆续提出,如比萨大学开发的基于快速傅里叶变换的方法(FFTBM)[8],法国开发的输入参数经验确定方法(DIPE)[8]和模型参数不确定性计算方法(CIRCÉ)[9],韩国原子力研究院开发的通过数据同化进行模型标定方法(MCDA)[10].李冬[11]对CIRCÉ方法可评估参数有限和MCDA方法求解非线性问题耗时过长这两个缺点进行了改进;Xiong等[12]开发了一种优化的最佳估算加不确定性分析方法,该方法包含了一个结构化的模型不确定性量化方法.
子通道程序已被广泛应用于堆芯热工水力分析,相较于计算流体力学(CFD),子通道分析消耗的计算资源更少,能够快速提供分析所需的热工参数.由于子通道程序对于物理模型存在一定简化假设以及数值计算误差等因素的影响,程序直接计算出的结果存在一定的不确定性,而目前针对子通道程序的不确定性分析还较少.本研究使用基于输入参数不确定性传递的统计学方法,利用期望最大化的思想来反推输入模型参数不确定性的分布,对子通道程序COBRA-IV 的计算进行不确定性量化分析.
1 不确定性量化分析方法
不确定性分析方法包括对影响计算程序结果的所有不确定性因素的确认和描述,并用统计学方法将这些不确定性因素整合得到输出结果的总体不确定性,其一般过程如图1所示.首先从众多输入参数中选取对计算结果影响较大的输入参数;其次量化这些输入模型参数的范围或不确定性分布;再次通过蒙特卡洛抽样或拉丁超立方抽样对这些输入参数进行抽样,并将这些抽样的输入参数输入子通道程序进行计算;最后对大量目标参数值进行统计,以得到输出目标参数的不确定性.其中,输入参数不确定性的确定和对输入参数进行抽样是关键步骤.
图1
1.1 输入模型参数不确定性量化
在分析过程中, 认为输出不确定性y是真实值yr与程序计算值yc之间的误差.实际可以得到的数据是实验值ye与对应的程序计算值yc,两者之差可以表示为
式中:e=(ej), j=1, 2,…, J 为实验误差,其中J为样本点数.
一般来说,模型参数的不确定性分析通常要确定模型参数基准值的无量纲乘数,记为k=(ki), i=1, 2, …, p,其中p为模型参数的个数;将待评估的模型参数不确定性记为x=(xi), i=1, 2, …, p.一般情况下, ki与xi存在以下关系:
因此可以采用ki=
输入模型参数不确定性x和输出不确定性y之间的传递关系可简化为
式中:u为已知分布的输入参数不确定性.一般情况下,正态分布是最常见的分布方式,因此假设x和u均服从正态分布.同时认为样本之间实验误差互相独立,实验误差满足均值为0且方差已知的正态分布,即ej~N(0,
近似认为x与y之间为线性关系,得到
式中:Sx和Su为样本点对应待求参数的敏感性矩阵,敏感性矩阵中的偏导数可以通过有限差分法求得[11].
当存在数据缺失问题时,期望最大化(EM)算法是极大似然估计的一种常用迭代算法.但是实际计算过程中EM算法一般迭代次数过多,收敛较慢,因此使用EM算法的改进算法,即ECME算法[13]进行计算.
将x的均值和方差记为μx和Cx,即待求参数,统一记为θ={μx, Cx}.完整数据的似然函数表示为L(θ; x1, x2, …, xJ, yw),简化表示为[11]
式中:xj为在第j个样本点对应的输入模型参数不确定性x的取值.由于L与 ln L 同时达到最大,所以求 ln L 达到最大值的点即可.考虑第m次迭代,令μx=
再由μx=
式中:ywi为yw的第j个分量.最后最大化似然函数,可以得到
在实际计算过程中,给定初值
由于以上是基于正态分布假设的,所以为了证明结果的可靠性需要进行假设检验.对残差进行检验,对于每个实验样本点,残差为
如果残差服从标准正态分布,则说明正态分布假设成立[8].
为测试该方法的可靠性,建立简单的数学模型对该方法进行测试.在测试中给定70个样本点,令已知分布的输入参数不确定性u的维度为1,待评估的模型参数不确定性x的维度为2,敏感性系数由均匀分布U(5, 50)得到,已知分布参数服从 u1~N(1.5, 0.42),认为实验误差为ej~N(0, 0.052),用于生成yw的真实分布满足x1~N(-0.6, 0.32), x2~N(1.2, 0.22),然后根据这些数据反推x的均值和方差,其计算结果如表1所示.可知,计算结果与真实值比较接近,相对误差不超过11%.
表1 评估所得参数与真实值的比较
Tab.1
参数 | 真实值 | 计算值 | 相对误差绝对值/% |
---|---|---|---|
-0.6 | -0.647 8 | 7.96 | |
0.3 | 0.281 7 | 6.23 | |
1.2 | 1.137 1 | 5.24 | |
0.2 | 0.221 3 | 10.65 |
1.2 输入参数抽样
不确定性的正向传递可以通过蒙特卡洛抽样进行.对确定的所有重要输入参数同时抽样,根据次序统计理论,抽样数目只与输出结果的容忍区间以及置信水平有关.满足特定容忍区间的最小抽样数目由Wilks公式[14]确定,在双侧容忍区间的情况下,其形式为
式中:β为置信水平;γ为概率;Z为最小样本数目.对于统计方法,美国核管会(NRC)建议采用95%的置信水平,不确定性评估必须证明在β=95%时,95%的计算结果(容忍区间)在安全裕度内[15].因此,令β=γ=95%,对于双侧容忍区间可得Z=93,即抽样计算93次后即可满足“95/95准则”, 得到的最大值和最小值可以认为是双95%置信区间的上下界.
2 PSBT空泡份额棒束基准
图2
表2 组件参数
Tab.2
参数 | 取值 | ||
---|---|---|---|
组件 | B5(图示①) | B6(图示②) | B7(图示③) |
棒束排列 | 5×5 | 5×5 | 5×5 |
加热棒数 | 25 | 25 | 24 |
套管数 | 0 | 0 | 1 |
加热棒外径/mm | 9.50 | 9.50 | 9.50 |
套管外径/mm | - | - | 12.24 |
加热棒间距/mm | 12.60 | 12.60 | 12.60 |
轴向加热长度/mm | 3 658 | 3 658 | 3 658 |
流道内宽度/mm | 64.9 | 64.9 | 64.9 |
径向功率分布类型 | A | A | B |
轴向功率分布类型 | 均匀 | 余弦 | 余弦 |
搅混叶片格架数 | 7 | 7 | 7 |
非搅混叶片格架数 | 2 | 2 | 2 |
简单格架数 | 8 | 8 | 8 |
搅混叶片格架 位置/mm | 471, 925, 1 378, 1 832, 2 285, 2 739, 3 247 | ||
非搅混叶片格 架位置/mm | 2.5, 3 755 | ||
简单格架数 位置/mm | 237, 698, 1 151, 1 605, 2 059, 2 512, 2 993, 3 501 |
3为实验组件截面图和子通道的具体划分情况.实验数据为位于距离加热底端 2 216、2 669、3 177 mm共3个轴向位置的4个中心子通道的空泡份额平均值.在距离加热底端2 216 mm的底端测量点的空泡份额较低,空泡聚集在受热表面,因此实验确定的空泡份额将低于实际的空泡份额,具有较大误差.对此,仅选取 2 669、3 177 mm 两个轴向位置的测量值.根据计算,所选用稳态实验工况的数据范围为压力7.36~14.71 MPa、质量流量1 380.6~2 288.9 kg/(m2·s)、加热功率 1 920~3 536 kW、入口温度173.5~287.8 ℃.共选取30组工况,并从中选取25组工况进行输入模型参数的不确定性量化,利用另外5组对计算得到的不确定性分析结果进行评估,以验证计算结果的有效性.
图3
从边界条件和计算模型两方面分析影响空泡份额分布的因素.对于边界条件,PSBT基准报告中给出了边界条件的测量精度,参照文献[17]的处理方式:压力、入口温度、质量流量、热流密度的不确定性的分布均为均值为0的正态分布,其标准差依次为0.333%、0.133%、0.5%、0.333%.
表3 PSBT基准题计算模型及参数选择
Tab.3
模型或参数 | 内容 |
---|---|
过冷空泡份额模型 | Levy模型 |
空泡份额模型 | 滑移模型 |
两相摩擦倍增因子 | Armand模型 |
湍流交混模型 | EM模型 |
两相湍流交混倍增因子 | Beus模型 |
单相湍流摩擦系数 | 0.184Re-0.2 |
格架压降因子(MV/NMV/SS) | 1.0/0.7/0.4 |
横流压降因子 | 0.5 |
模型参数使用基准值进行计算.如图4所示,横轴为空泡份额实验值ve,纵轴为子通道程序空泡份额计算值vc,越接近黑色实线代表计算值与测量值越符合,上、下两条虚线代表实验测量不确定性2σ(=0.08).可知,较多的计算值位于实验值的两倍标准差(2σ)范围之外,且整体预测值偏高,证明计算选用的模型过于保守,仍需要进一步改进.
图4
在进行模型参数不确定性量化过程中,计算敏感性矩阵是重要步骤.在使用有限差分法求偏导数的过程中,步长的选择并不是越小越好.因为步长太小可能会导致程序响应变化不明显,也会受到数值不稳定性的影响,步长过大又会导致误差过大.所以为了减小偏导数计算过程中的误差,需要进行平滑处理[11].例如在求解滑速比和湍流交混系数对应的敏感性系数时,分别取步长为 ±0.01, ±0.03, ±0.1 来进行计算,对于每个样本点可以得到6个偏导数,然后舍掉最大值和最小值,求平均值作为最终的偏导数.
3 计算结果分析
将压力、入口温度、质量流量、热流密度当作已知分布参数,将滑速比和湍流交混系数作为待评估的模型参数,利用25组工况共50个实验样本点计算得到滑速比的均值和标准差分别为 0.440 5、0.259 9;湍流交混系数的均值和标准差分别为-0.103 4、0.048 8.可知,滑速比不确定性的标准差较大.根据最小熵增理论,滑速比的大小与压力密切相关,而在进行模型参数的不确定性分析时所选用的工况压力范围为7.36~14.71 MPa,因而滑速比也存在较大的不确定性.为了验证正态分布假设是否成立,使用卡方检验来验证残差是否满足标准正态分布,将残差划分为10个区间,有:
则有95%的置信度认为参数服从正态分布.
在得到了输入参数不确定性的具体分布之后,对另外5组工况进行不确定性的正向传播.将得到的模型参数及边界条件参数按照其分布进行93次抽样,并通过子通道程序进行计算,将得到的最大值和最小值当作双95%置信区间的上下界.图5为4个中心子通道的空泡份额平均值计算值不确定带对实验值的包络情况.可知,当空泡份额较高时,不确定带对实验值的包络性非常好;当空泡份额较低时,不确定带对实验数据的包络情况不够理想,有部分未包络,但不确定带下界也非常接近实验值.
图5
图5
不确定带上下界对实验值的包络情况
Fig.5
Envelope of experimental values with uncertain upper and lower bounds
利用所求出的统计均值对模型系数进行修正,图6为基准模型修正前后4个中心子通道的空泡份额平均值计算结果对比.可知,样本点模型参数经过均值修正后,计算结果与实验值更接近,计算精度有了明显改善.
图6
4 结论
本文在假设输入参数不确定性与输出参数不确定性近似为线性关系的基础上,将边界条件压力、入口温度、质量流量、热流密度当作已知分布输入参数,对子通道程序COBRA-IV输入模型参数中滑速比和湍流交混系数的不确定性进行评估;然后进行输入参数不确定性的正向传递,得到了空泡份额的不确定带,并利用求得的模型参数不确定性统计均值对模型进行修正,计算结果表明:
(1) 在所选取的空泡分布实验工况范围内,滑速比存在较大的不确定性,而湍流交混系数的变化范围较小.
(2) 计算结果的不确定带对实验值的包络情况较好,尤其是当空泡份额较高时.
(3) 利用统计均值对模型进行修正后,可以得到比原模型更接近实验值的计算结果,对于输入模型参数的不确定性评估结果可以作为计算模型改进的重要参考.
参考文献
Quantifying reactor safety margins part 1: An overview of the code scaling, applicability, and uncertainty evaluation methodology
[J].DOI:10.1016/0029-5493(90)90071-5 URL [本文引用: 1]
Uncertainty and sensitivity analysis of a post-experiment calculation in thermal hydraulics
[J].DOI:10.1016/0951-8320(94)90073-6 URL [本文引用: 1]
Outline of the uncertainty methodology based on accuracy extrapolation
[J].DOI:10.13182/NT109-21 URL [本文引用: 1]
Thermal-hydraulics system codes uncertainty assessment: A review of the methodologies
[J].
Bemuse phase Ⅵ report: Status report on the area, classification of the methods, conclusions and recommendations
[R].
Quantification of the uncertainty of the physical models in the system thermal-hydraulic codes-PREMIUM benchmark
[J].DOI:10.1016/j.nucengdes.2019.110199 URL [本文引用: 6]
Identifying intrinsic variability in multivariate systems through linearized inverse methods
[J].DOI:10.1080/17415971003624330 URL [本文引用: 1]
Optimization of thermal-hydraulic reactor system for SMRs via data assimilation and uncertainty quantification
[J].DOI:10.13182/NSE11-113 URL [本文引用: 1]
Investigation of uncertainty quantification methods for constitutive models and the application to LOFT LBLOCA
[J].DOI:10.1016/j.anucene.2019.04.028 URL [本文引用: 1]
The ECME algorithm: A simple extension of EM and ECM with faster monotone convergence
[J].DOI:10.1093/biomet/81.4.633 URL [本文引用: 1]
Determination of sample sizes for setting tolerance limits
[J].DOI:10.1214/aoms/1177731788 URL [本文引用: 1]
An overview of Westinghouse realistic large break LOCA evaluation model
[J].
Assessment of the uncertainties of COBRA sub-channel calculations by using a PWR type rod bundle and the OECD NEA UAM and the PSBT benchmarks data
[J].DOI:10.3139/124.110460 URL [本文引用: 1]
Estimation of steady-state steam void-fraction by means of the principle of minimum entropy production
[J].DOI:10.1115/1.3687113 URL [本文引用: 1]
Single-phase subchannel mixing in a simulated nuclear fuel assembly
[J].
/
〈 |
|
〉 |
