子通道程序对PSBT空泡分布实验计算的不确定性量化分析

张俊涛， 刘晓晶， 张滕飞， 柴 翔

(上海交通大学 核科学与工程学院，上海 200240)

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 输入模型参数不确定性量化

在分析过程中, 认为输出不确定性y是真实值yr与程序计算值yc之间的误差.实际可以得到的数据是实验值ye与对应的程序计算值yc，两者之差可以表示为

yw=ye-yc=(ye-yr)+(yr-yc)=e+y

(1)

式中：e=(ej),j=1, 2，…,J为实验误差，其中J为样本点数.

一般来说，模型参数的不确定性分析通常要确定模型参数基准值的无量纲乘数，记为k=(ki),i=1, 2, …,p，其中p为模型参数的个数；将待评估的模型参数不确定性记为x=(xi),i=1, 2, …,p.一般情况下，ki与xi存在以下关系：

(2)

因此可以采用ki=exi来表示模型参数系数和不确定性之间的转化关系[8].

输入模型参数不确定性x和输出不确定性y之间的传递关系可简化为

y=h(x,u)

(3)

近似认为x与y之间为线性关系，得到

yw=Sxx+Suu+e

(4)

式中：Sx和Su为样本点对应待求参数的敏感性矩阵，敏感性矩阵中的偏导数可以通过有限差分法求得[11].

当存在数据缺失问题时，期望最大化(EM)算法是极大似然估计的一种常用迭代算法.但是实际计算过程中EM算法一般迭代次数过多，收敛较慢，因此使用EM算法的改进算法，即ECME算法[13]进行计算.

将x的均值和方差记为μx和Cx，即待求参数，统一记为θ={μx,Cx}.完整数据的似然函数表示为L(θ;x1,x2, …,xJ,yw)，简化表示为[11]

L∝

(5)

(6)

(7)

式中：ywi为yw的第j个分量.最后最大化似然函数，可以得到

(8)

由于以上是基于正态分布假设的，所以为了证明结果的可靠性需要进行假设检验.对残差进行检验，对于每个实验样本点，残差为

(9)

如果残差服从标准正态分布，则说明正态分布假设成立[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 Comparison of evaluated parameters and real values

1.2 输入参数抽样

不确定性的正向传递可以通过蒙特卡洛抽样进行.对确定的所有重要输入参数同时抽样，根据次序统计理论，抽样数目只与输出结果的容忍区间以及置信水平有关.满足特定容忍区间的最小抽样数目由Wilks公式[14]确定，在双侧容忍区间的情况下，其形式为

β=1-γZ+Z(1-γ)γZ-1

(10)

式中：β为置信水平；γ为概率；Z为最小样本数目.对于统计方法，美国核管会(NRC)建议采用95%的置信水平，不确定性评估必须证明在β=95%时，95%的计算结果(容忍区间)在安全裕度内[15].因此，令β=γ=95%，对于双侧容忍区间可得Z=93，即抽样计算93次后即可满足“95/95准则”, 得到的最大值和最小值可以认为是双95%置信区间的上下界.

2 PSBT空泡份额棒束基准

压水堆子通道和棒束实验(PSBT)是2010年发布的一套基准，主要用于验证现有子通道和计算流体力学程序模拟流体在单管和棒束的空泡份额、含气率以及偏离泡核沸腾(DNB)等.本研究计算对象为PSBT基准一阶段稳态棒束空泡份额分布问题.表2为本文计算涉及的3种组件具体参数[16].其中，B5和B6组件的径向功率分布场为类型A，如图2(a)所示；B7组件的径向功率分布为类型B, 如图2(b)所示.

图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组对计算得到的不确定性分析结果进行评估，以验证计算结果的有效性.

表2 组件参数Tab.2 Parameters of assemblies

①②③

从边界条件和计算模型两方面分析影响空泡份额分布的因素.对于边界条件，PSBT基准报告中给出了边界条件的测量精度，参照文献[17]的处理方式：压力、入口温度、质量流量、热流密度的不确定性的分布均为均值为0的正态分布，其标准差依次为0.333%、0.133%、0.5%、0.333%.

在子通道程序的计算过程中所选的重要模型和参数如表3所示.空泡份额模型采用滑移模型，滑速比可以设置.根据最小熵增原理[18]，在14 MPa下的滑速比计算结果为1.86，但是一般情况下计算值比大多数的实验值偏高[19]，而最典型的单相流模型是均相模型，其把两相流体看成某种单相混合物流体，认为气相与液相之间不存在相对速度，即滑速比为1，因此滑速比基准值最终取1.2.湍流交混模型中需要给定湍流交混系数，湍流交混系数主要与几何尺寸、棒束排布方式、流体类型以及雷诺数(Re)范围有关，根据与计算工况相近的Castellana等[20]的实验所得湍流交混系数关系式，取湍流交混系数的基准值为0.008.

表3 PSBT基准题计算模型及参数选择

模型参数使用基准值进行计算.如图4所示，横轴为空泡份额实验值ve，纵轴为子通道程序空泡份额计算值vc，越接近黑色实线代表计算值与测量值越符合，上、下两条虚线代表实验测量不确定性2σ(=0.08).可知，较多的计算值位于实验值的两倍标准差(2σ)范围之外，且整体预测值偏高，证明计算选用的模型过于保守，仍需要进一步改进.

在进行模型参数不确定性量化过程中，计算敏感性矩阵是重要步骤.在使用有限差分法求偏导数的过程中，步长的选择并不是越小越好.因为步长太小可能会导致程序响应变化不明显，也会受到数值不稳定性的影响，步长过大又会导致误差过大.所以为了减小偏导数计算过程中的误差，需要进行平滑处理[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个区间，有：

(11)

则有95%的置信度认为参数服从正态分布.

在得到了输入参数不确定性的具体分布之后，对另外5组工况进行不确定性的正向传播.将得到的模型参数及边界条件参数按照其分布进行93次抽样，并通过子通道程序进行计算，将得到的最大值和最小值当作双95%置信区间的上下界.图5为4个中心子通道的空泡份额平均值计算值不确定带对实验值的包络情况.可知，当空泡份额较高时，不确定带对实验值的包络性非常好；当空泡份额较低时，不确定带对实验数据的包络情况不够理想，有部分未包络，但不确定带下界也非常接近实验值.

利用所求出的统计均值对模型系数进行修正，图6为基准模型修正前后4个中心子通道的空泡份额平均值计算结果对比.可知，样本点模型参数经过均值修正后，计算结果与实验值更接近，计算精度有了明显改善.

4 结论

本文在假设输入参数不确定性与输出参数不确定性近似为线性关系的基础上，将边界条件压力、入口温度、质量流量、热流密度当作已知分布输入参数，对子通道程序COBRA-IV输入模型参数中滑速比和湍流交混系数的不确定性进行评估；然后进行输入参数不确定性的正向传递，得到了空泡份额的不确定带，并利用求得的模型参数不确定性统计均值对模型进行修正，计算结果表明：

(1) 在所选取的空泡分布实验工况范围内，滑速比存在较大的不确定性，而湍流交混系数的变化范围较小.

(2) 计算结果的不确定带对实验值的包络情况较好，尤其是当空泡份额较高时.

(3) 利用统计均值对模型进行修正后，可以得到比原模型更接近实验值的计算结果，对于输入模型参数的不确定性评估结果可以作为计算模型改进的重要参考.