基于确定性抽样的过冷沸腾边界条件不确定性分析

张 翔，彭敏俊，丛腾龙，李孝佳，陈衣然

(哈尔滨工程大学 核安全与仿真技术国防重点学科实验室，黑龙江 哈尔滨 150001)

过冷沸腾具有换热系数高的特点，但随着热流密度的增加，可能达到沸腾临界，导致壁温飞升损坏设备[1]。近年来随着计算机技术的发展和对两相流机理的认识，计算流体力学(CFD)方法越来越多地应用于反应堆热工水力和流体流动传热的研究，并已得到认可[2]，因此，基于CFD方法的过冷沸腾流动研究引起广泛关注，包括基准实验模拟[3-5]、汽泡动力学研究[6-8]、模型对比验证[9-11]、棒束通道流动传热[12-13]等。

在上述的CFD模拟中，边界条件、模型经验系数等均作为确定值输入，计算结果无法反映在实际边界条件测量和模型系数简化过程中引入的不确定性带来的影响。然而Bestion等[14]指出，在沸腾流动的数值模拟中，边界条件不确定性和模型不确定性的影响不可忽略。因此，需在过冷沸腾的CFD模拟中引入不确定性分析，提高结果的可信度。Badillo等[15]和Prošek等[16]分别利用多项式混沌展开和随机抽样对单相GEMIX搅混实验的入口边界条件进行了不确定性分析。Dunn等[17]和Hedberg等[18]分别利用拉丁超立方抽样和确定性抽样对湍流模型在单相突扩流中的不确定性进行了分析。Zhang等[19]和Cong等[20]针对Bartolomei过冷沸腾实验进行了数值模拟，并分别对入口边界条件和沸腾模型进行了不确定性分析，但研究中采用的拉丁超立方抽样法(分别需要740、1 040个样本数据)产生的样本数较多，计算成本较高，不适合应用于复杂问题的计算，且二者实验数据均为截面平均值在轴向的分布，均未研究不确定性传递对沸腾流动参数径向分布的影响。

本文利用Fluent软件构建数学物理模型，对DEBORA实验[21]进行数值模拟，采用确定性抽样方法对边界条件输入参数进行抽样计算，获得主要径向参数分布的期望值和置信区间，并与实验数据对比验证，研究不确定性对过冷沸腾流动径向参数的影响。之后进一步分析得到每个不确定性源的影响权重，为模型改进方向提供参考。

1 数学物理模型

基于欧拉两流体模型，对过冷沸腾中的两相分别建立质量、动量和能量守恒方程，结合标准k-ε湍流模型、RPI壁面沸腾模型、相间作用力以及一些辅助模型对竖直圆管内过冷沸腾流动建立数值模型。

1.1 壁面沸腾模型

RPI壁面沸腾模型是由Kurul等[22]提出的，用于计算壁面向流体传递的热流。该模型将总热流密度qw分为单相对流热流密度qC、骤冷换热热流密度qQ和蒸发热流密度qE3个部分，即：

qw=qC+qQ+qE

(1)

qC=hC(Tw-Tl)(1-Ab)

(2)

(3)

qE=VdNwρghfgf

(4)

式中：hC为单相对流换热系数；Tw为壁面温度；Tl、kl、cp,l、ρl为液体温度、导热率、比热容、密度；ρg为气体密度；hfg为蒸发潜热；Vd为汽泡体积，可根据汽泡脱离直径Dw计算；Ab、f和Nw分别为汽泡影响面积、汽泡脱离频率和汽化核心密度。

壁面沸腾模型中选用的具体子模型如下：汽泡脱离直径Dw，Tolubinski-Kostanchuk模型[23]；汽泡影响面积Ab，Delvalle-Kenning模型[24]；汽泡脱离频率f，Cole模型[25]；汽化核心密度Nw，Lemmert-Chawla模型[26]。

1.2 相间模型

两相间由于速度差等因素存在相互作用力，主要包括曳力、升力和湍流耗散力，这些相间作用力通过在动量方程中引入源项来表示。

作用于汽泡上的曳力FD由两相相对运动引入，主要由局部流动参数(包括空泡份额αg、汽泡直径dg、汽泡速度vg和液体速度vl)计算得到：

(5)

式中，CD为曳力系数，可由Schiller-Naumann模型[27]计算得到。

在竖直圆管内过冷沸腾流动中，由于液相在径向的速度梯度不均匀，会对其中的汽泡产生升力FL作用，作用方向垂直于流动方向，表示为：

(6)

式中，CL为升力系数，可由Moraga模型[28]计算得到。

湍流耗散力用来描述液相湍流流动对汽泡的耗散作用。在湍流耗散力作用下，汽泡将从近壁面区域被带走，对径向空泡份额分布起到展平的作用。湍流耗散力FTD由Burns模型[29]计算：

(7)

式中：Kgl为相间动量交换系数；Dl为液体耗散标量；CTD为湍流耗散力系数；σgl为耗散普朗特数。

近壁面区域的相间质量交换率mE和主流区的相间质量交换率m分别由式(8)、(9)计算。

(8)

(9)

式中，ql和qg分别为两相界面到液体和汽泡的能量传递，可由式(10)、(11)计算。

ql=hgl(Tg-Tl)

(10)

qg=αgρgcp,g(Tsat-Tg)/δt

(11)

式中：hgl为体积传热系数，由Ranz-Marshall模型[30]计算；Tg和cp,g分别为气体温度和比热容；δt为时间尺度，取值为0.05[31]。

2 DEBORA实验工况及数值建模

Garnier等[21]对竖直圆管内R-12工质自下而上的过冷沸腾流动进行了研究，简称DEBORA实验。实验中加热管长3.5 m、圆管内径19.2 mm，实验测得了空泡份额、液体温度、汽泡直径等参数的径向分布，可用于数据验证。为排除单个工况的计算结果巧合，本文选取了4个实验工况[3,32-33](表1)。实验工况下，2.62 MPa的R-12具有与15.7 MPa高压水相近的沸腾相关无量纲数，可反映高压水的流动沸腾换热特性[3]。模拟中边界条件输入包括质量流速G、入口温度T、热流密度q和压力p，这些变量在测量过程中由于测量误差会引入不确定性，对模拟结果产生影响。此外，湍流强度I作为速度入口边界中的边界值，在充分发展段的数值计算中，若根据经验假定的充分发展段管长不足以使流体完全充分发展，可能会对出口结果产生影响，进而对过冷沸腾段的计算产生影响，因此湍流强度也作为不确定性源考虑。根据Fluent手册[31]推荐和有关数值模拟边界条件不确定性研究中湍流强度的选取[17,34]，模拟中假设其值为0.05。综上，本文主要分析这5个边界条件作为不确定性源对过冷沸腾径向参数计算结果的影响。假设5个边界条件相互独立且其概率密度函数服从正态分布，假设的输入变量在3σ区间内的误差列于表2，即5个输入变量低于此误差的概率为99.74%。

表2 输入变量在3σ区间内的误差Table 2 Error of input variable in 3σ region

考虑到实验圆管的对称性，为节省计算资源，截取圆管的二维轴对称平面作为计算域。采用压力-速度耦合算法Coupled进行稳态计算，使用伪瞬态方法帮助加速收敛，对动量、空泡份额、湍动能、湍流耗散率、能量采用一阶迎风的空间离散格式，对梯度采用基于单元体的最小二乘法插值离散格式。模拟实验中的迭代误差可通过计算获得充分收敛的解来忽略。此外，为降低网格离散对计算结果带来的影响，基于Case4工况针对不同网格(表3)进行了网格敏感性分析(图1，其中r表示径向位置，R表示圆管半径)。可看到，mesh3和mesh4、mesh6的结果基本一致，可认为达到网格无关解。综合考虑模拟结果的准确性和计算成本，选取mesh3对DEBORA实验进行数值模拟，并分析其边界条件不确定性对结果的影响。

表3 网格信息Table 3 Mesh information

图1 网格敏感性分析Fig.1 Grid sensitivity analysis

3 不确定性分析

3.1 确定性抽样方法

与蒙特卡罗随机抽样和拉丁超立方抽样通过抽取大量样本描述变量概率密度函数的方式不同，确定性抽样通过传递变量统计矩的方式抽取样本，使得样本数大幅减少，降低计算成本[35]。根据满足的统计矩阶数不同，确定性抽样具有不同抽样形式。按照顺序，前4阶统计矩分别为均值、方差、偏度和平直度。本文采用满足前4阶矩的抽样形式，并假设所研究的变量相互独立且其概率密度函数满足正态分布，可求出对于变量q，满足前4阶矩的抽样点分别为2个周边点qi1、qi3以及1个中心点qi2，并赋予相应的权重因子w[18,36]。

(12)

(13)

式中：μ和σ分别为平均值和标准差；下标i表示第i个不确定性源。

为合并所有中心点，中心点的权重因子需重新计算，计算公式如下：

(14)

式中，N为不确定性源的数目。

通过上述方法，针对表2中列举的带有不确定性的5个边界条件抽取样本，得到的样本及权重因子列于表4。

表4 抽样样本及权重Table 4 Sample and weight

3.2 边界条件不确定性对过冷沸腾模拟的影响

(15)

(16)

UB=3δφ

(17)

图2 边界条件不确定性对径向参数的影响Fig.2 Effect of boundary condition uncertainty on radial parameter

图2中空泡份额沿壁面向管道中心方向，呈逐渐下降趋势，这是由于流体在壁面受热蒸发产生汽泡，部分汽泡在相间力的作用下向管道中心运动，而管道主流区域处于过冷状态，汽泡运动过程中会受到冷凝作用。流道中心的空泡份额置信区间相比于壁面附近较小，说明和壁面附近相比，管道中心空泡份额的计算受边界条件不确定性的影响较低。相反，根据径向3σ区间的带宽分布来看，壁面附近气体速度和汽泡直径的不确定度带宽达到最小值，说明壁面附近气体速度和汽泡直径的计算对边界条件不确定性的敏感度较低。此外，边界条件不确定性对液体温度和液体速度的影响在径向上分布较为均匀，置信区间的带宽变化不大。

壁面轴向位置的最大空泡份额是判断临界热流现象发生的重要参数[39-40]，该参数的准确计算对预测临界热流密度具有重要意义。然而随着边界条件不确定性的引入，基于确定单一输入值的CFD方法得到壁面最大空泡份额，并将其用于预测临界热流密度的计算方式，不足以得到足够可靠的预测值。在考虑边界条件不确定性基础上，计算得到了壁面最大空泡份额期望值和3σ置信区间，如图3所示。Case1、Case2、Case3和Case4等4个工况的置信区间带宽分别为0.087、0.089、0.090和0.089，与各自期望值相比，所占的百分比在19%～24%之间。因此，从模拟结果来看，边界条件不确定性在壁面最大空泡份额计算中的影响不可忽略，也说明在今后临界热流密度模拟预测中引入不确定性分析具有一定必要性。

图3 边界条件不确定性对壁面最大空泡份额的影响Fig.3 Effect of boundary condition uncertainty on maximum void fraction at wall

3.3 边界条件不确定性源的影响权重

基于不确定性源相互独立假设，每个不确定性源对过冷沸腾径向参数的影响权重Ci可通过式(18)计算，计算结果如图4所示，图中不同径向位置用不同符号表示，为避免结果的偶然性，针对表1中的4个工况分别进行不确定性源的影响权重计算。

a——空泡份额；b——液体温度图4 不确定性源对径向参数的影响权重Fig.4 Uncertainty source contribution on radial parameter

(18)

UBi=3δφi

(19)

(20)

式中，δφi和UBi分别为第i个不确定性源的标准差和不确定度带宽。

对于径向空泡份额(图4a)，4个工况中不确定性源的影响权重计算结果基本相同，液体入口温度的不确定性起主要影响，在3个测量位置其影响权重分别约为30%、37%和42%。另外，湍流强度边界不确定性的影响权重几乎可忽略不计，这是因为本文考虑了流体的充分发展，将充分发展段出口参数作为加热沸腾段的入口边界，湍流强度边界影响足够小说明加热沸腾段的入口流体已充分发展。随着由圆管中心到管道壁面的径向位置的改变，可看到液体入口温度和运行压力不确定性的影响权重逐渐减小，而质量流速和壁面热流密度不确定性的影响逐渐升高。值得注意的是，对于壁面空泡份额(r/R=1)，液体入口温度和壁面热流密度不确定性的影响权重均达到约30%，因此，在依据壁面最大空泡份额预测临界热流密度的模拟实验中，液体入口温度和壁面热流密度须谨慎处理。

与空泡份额不同，在不确定性源对液体温度影响权重的计算中，不同工况的计算结果略有差异，但影响权重随径向位置变化的趋势相同(图4b)。这是由于4个工况实验条件不同，Case1、Case2、Case4、Case3的出口平衡含气率依次上升(表1)。从图4b可看出，随着由流道中心到四周的径向位置的改变，液体入口温度、壁面热流密度和质量流速不确定性的影响权重均逐渐减小，而运行压力不确定性的影响权重逐渐增加。这说明在圆管中心，液体温度受入口流体温度不确定性的影响较大，而靠近壁面区域的液体温度主要受运行压力的影响。此外，随着不同工况出口平衡含气率的上升，运行压力不确定性对壁面附近液体温度的影响也逐渐增大，对于4个工况中出口平衡含气率最高的Case3，运行压力的影响权重甚至达到约70%。

4 结论

本文采用确定性抽样方法，分析了边界条件不确定性对圆管过冷沸腾径向参数分布的影响。研究选用流动工质为R-12的DEBORA实验的4个工况作为模拟的基准实验。实验工况下，2.62 MPa的R-12具有与高压水相近的沸腾相关无量纲数，可反映高压水的流动沸腾换热特性。此外，计算得到了每个边界条件不确定性源的影响权重。得到如下结论：

1) 通过建立过冷沸腾模拟的数值模型，得到了边界条件不确定性对径向参数分布的影响趋势。与实验值相比，本文建立的模型可较好地模拟DEBORA实验中的空泡份额和流体温度径向分布，但对于气相速度和汽泡直径的计算存在一定偏差。

2) 通过计算边界条件不确定性影响下壁面最大空泡份额的期望值和置信区间，发现边界条件不确定性在壁面最大空泡份额计算中的影响不可忽略。

3) 通过对比分析边界条件不确定性源的影响权重，发现边界条件中流体入口温度的不确定性对径向空泡份额的计算影响较大，另外越靠近壁面位置，壁面热流密度的不确定性对空泡份额计算的影响越大。而运行压力和流体入口温度的不确定性是影响径向液体温度计算的主要因素。