核电站中气溶胶再悬浮的CFD研究

江 斌，黄 挺，陈 炼，常华健,2

(1.国核华清(北京)核电技术研发中心有限公司，北京 102209；2.清华大学 核能与新能源技术研究院，北京 100084)

压水堆核电站在严重事故晚期时，由于可能出现的氢气燃烧或爆燃等现象造成的瞬时气流冲击，会使得已沉积在安全壳壁面内的气溶胶再次悬浮到气空间中，导致安全壳气空间内裂变产物浓度再次上升，从而可能导致放射性物质的再次释放。由于气溶胶再悬浮会对事故晚期放射性源项产生影响，最终影响释放到环境中的放射性物质总量，因此有必要对该类现象进行深入研究。

已有用于反应堆严重事故评估的系统分析程序如ASTEC等[1]，根据不同原理开发了气溶胶再悬浮机理的计算模型。但该类程序是集总参数程序，不能精细获悉气溶胶不同时刻粒子的位置分布等物理量，这不利于了解事故细节。

GASFLOW是一款用于模拟流体动力、传热传质、化学反应、气溶胶输运与沉降及其他严重事故时反应堆内相关现象的CFD程序[3]。GASFLOW采用隐式连续欧拉-任意拉格朗日欧拉(ICE’d ALE)的方法求解三维可压缩多组分的纳维-斯托克斯方程，该方法属于有限差分法，最初由Hirt等[4]提出，可适用于求解所有流速的计算流体力学问题。每个版本的GASFLOW都经过了国际上多项标准实验的验证，其计算结果也被多个亚洲和欧洲国家的权威机构所认可。拉格朗日方法描述了物质颗粒的运动，它是时间和空间的函数；欧拉方法描述了空间固定体积中的流动情况，求解其控制方程相对于拉格朗日方法更容易。ICE’d ALE方法能求解包括压声速和超声速在内的所有流动问题，因为其具有与隐式连续欧拉法类似的隐式方程。其优势包括可采用比显式方法大得多的时间步长，以及能求解全拉格朗日方法无法求解的问题。控制方程可利用有限体积法在经典的交错网格上求解，其中，流体的状态量通过控制体中心网格求解，而速度矢量通过控制体面网格求解。由于采用了非贴体网格和亚网格技术，以及简化的工程建模方法，使得GASFLOW在计算大型安全壳内的流场时具有很高的效率。

本文以GASFLOW3.3为基础，对气溶胶粒子进行受力分析，利用力学平衡原理建立气溶胶再悬浮模型，并与国际标准例题ISP40的STORM试验SR11[5]结果进行对比，在此基础上对GASFLOW中再悬浮模型的改进提出合理化建议。

1 模型描述

根据经济合作与发展组织(OECD)能源署[6]的有关定义，再悬浮是指由于流体流动的拖曳力，结构表面已经沉积的气溶胶或凝聚物以细小碎片的形式悬浮起来。再悬浮也可由重新定位的熔体上注水产生的蒸汽尖峰或氢燃烧产生的流动所触发。

目前，有关气溶胶再悬浮的模型通常分为两类，一类为力学平衡(FB)模型，另外一类为能量累积(EA)模型[7]。GASFLOW中用于模拟气溶胶行为的模型包括拉格朗日粒子模型、随机湍流粒子扩散模型、粒子沉降模型、粒子再悬浮模型及粒子云模型[8]。GASFLOW3.3采用了FB模型，对单一气溶胶粒子进行受力分析，且为计算简便，将粒子近似为球体。实际两相流动中，经常采用基于体积的当量球半径来计算粒子拖曳力[9]。图1为稳定发展的湍流场中，沉积于壁面上的单一球形粒子的受力分析示意图。

对于单一球形粒子，其受到的各力表达式如下。

(1)

(2)

(3)

(4)

(5)

Ff=fsFn=fs(Fg+Fa-Fb-Fl)

(6)

目前，我国所开发利用的钼资源一般来自单一钼矿石类型、铜钼矿石类型和钨钼矿石类型等。从以铁为主的铁钼型矿石中回收钼的生产实践尚不多见。本技术为从铁钼型矿石中回收钼提供了一条切实可行的技术路线，对处理同类型性质矿石具有很大的借鉴价值，具有较好的推广前景。

图1 稳定且完全发展的湍流场中单一球形粒子的受力Fig.1 Force acting on single spherical particle with steady and fully developed turbulent flow

此外，流体中的粒子通常还会受到压力梯度力、虚拟质量力、Basset力、Magnus力、Saffman力和热泳力等作用。尽管作用在粒子上的力十分复杂，但在气溶胶悬浮体系中，由于气体的密度通常远小于粒子的密度，与粒子本身惯性相比，浮力、压力梯度力和虚拟质量力等均很小，可忽略不计[12]。在以上两相间的相互作用力中，拖曳力是十分重要的力，它是对两相间的相对流动、传热和传质等起着重要作用的力。根据Evett等[13]的公式可计算出拖曳力。

上述受力中，重力、浮力和黏性力仅取决于粒子和气体的物理性质，与气流速度无关。主要的黏性力是范德瓦尔斯力、静电力和其他黏性力。这些力主要受到下列因素的影响：粒子的材料、形状和尺寸，表面的材料、粗糙度和污染程度；环境气体的相对湿度和温度；粒子-表面接触的持续时间和初始接触速度。而拖曳力和升力取决于气流速度且随气流速度的增大而增大。当气流速度增大到一定程度，拖曳力等于静摩擦力时，气溶胶粒子开始运动并离开沉积表面，此时：

Fd=fs(Fg+Fa-Fb-Fl)

(7)

将式(1)～(5)代入式(7)中，可得：

(8)

式(8)是以Ugcp为未知数的方程，通过求解该方程可得出单一粒子的最小悬浮速度，将其用Ugcp0表示。当气流速度大于该速度时，气溶胶粒子就可悬浮起来。

2 几何建模及计算条件

2.1 几何建模

STORM试验台架位于意大利的Ispra，是由联合研究中心JRC运行维护的。STORM试验SR11分为两个独立阶段，第1阶段集中在气溶胶粒子由于热泳和湍流影响导致的沉积现象，第2阶段关注的是已沉积的气溶胶粒子在逐步增加的气流下导致的再悬浮现象。

GASFLOW是一个适合于大尺度、安全壳内的工程CFD程序，对于细长小管径是否适用需进一步研究。为便于对模型适用性进行分析，利用GASFLOW建立一长方形管道模型，其方形管口尺寸为0.5 m×0.5 m，长度为5.005 5 m，初始条件与SR11设置完全一致。设定初始状态下气溶胶的沉积区域如图2所示。

在保证气溶胶充分沉积后，在管口处设定一持续1 s的速度边界，使管内产生一水平方向的流场，该风流将使沉积的气溶胶产生再悬浮，再悬浮率η的计算公式为：

(9)

其中：N0为初始状态下沉积在表面的气溶胶颗粒数量；Nr为气流停止时沉积在表面的气溶胶颗粒数量。

图2 气溶胶初始沉积区域示意图Fig.2 Schematic diagram of aerosol initial deposition area

2.2 初始条件

表1列出计算分析的初始条件。其中边界层厚度、拖曳力系数粗糙度因子和颗粒与沉积表面的平衡间距粗糙度因子是GASFLOW3.3中与再悬浮模型相关的3个主要输入参数。核电站在严重事故条件下，安全壳内气溶胶粒子中包含了各种核素成分，粒径的典型值为1.0 μm[14]，因此本文选取该典型直径进行分析。

表1 计算的初始条件Table 1 Initial condition for calculation

由于气溶胶粒子运动具有随机性，因此数目越多的粒子数量越能更好地反映再悬浮的随机特性。但太多粒子数量会占用大量的计算资源，因此需对初始气溶胶粒子的数量进行合理选定。

图3示出在表1设定的初始条件下，气溶胶再悬浮率随气溶胶粒子数量的变化。由图3可看出，当气溶胶粒子数量大于等于1×105时，气溶胶再悬浮率的计算结果基本保持不变，因此选定气溶胶粒子数量为1×105个进行后续的分析计算。

图3 不同气溶胶粒子数量的再悬浮率Fig.3 Resuspension rate for different aerosol particle numbers

GASFLOW是利用拉格朗日粒子追踪法进行计算分析的，为计算的简便，同时在工程角度适用的情况下，不考虑粒子间的相互作用，且粒子与流场间的作用和传热微乎其微，因此也不予考虑。

3 计算结果及分析

3.1 再悬浮率

图4示出STORM试验SR11的结果及ASTEC和GASFLOW程序的计算结果。其中FB-1、FB-2、RNR-1和RNR-2分别代表ASTEC程序中SOPHAEROS模块采用的两种不同模型计算的结果。FB-1和FB-2均采用FB模型，区别在于FB-1的黏性力系数为1.0×10-6N/m，FB-2的黏性力系数为1.0×10-5N/m。RNR-1和RNR-2均采用EA模型，区别在于RNR-1的最大时间步长为1.0 s，RNR-2的最大时间步长为30.0 s。GASFLOW-1和GASFLOW-2采用的是FB模型，区别在于GASFLOW-1的杨氏模量为2.1×1012dynes/cm2, GASFLOW-2的杨氏模量为2.1×1013dynes/cm2。从图4可看出，GASFLOW与ASTEC程序的计算结果有所不同，这是由于两者采用的是不同模型，且图5示出两个不同时刻气溶胶粒子的位置分布。由图5可看出，距流场越近再悬浮粒子的数量越多，这是由于距流场不同，所受到的升力大小也不同。

图4 再悬浮率试验与计算结果的对比Fig.4 Comparison of resuspension rate between test and calculated results

ASTEC是集总参数程序，GASFLOW是CFD程序，因此计算结果有所差别。GASFLOW-1和GASFLOW-2的计算结果也有所差别，这是由于两者的杨氏模量不同，因而导致两者的黏性力不同，再悬浮率也不同，这一结论与FB-1和FB-2的计算结果不同是一致的，说明黏性力是影响再悬浮率的一个重要因素。同时还可看出，GASFLOW-1和GASFLOW-2与试验结果SR11的差距在20%以内，在工程应用范围内符合较好，GASFLOW能较好模拟气溶胶的行为，可在工程应用中提供参考。

3.2 位置分布

图5 不同时刻气溶胶的位置分布Fig.5 Position distribution of aerosol at different time

图6 气溶胶再悬浮率随时间的变化Fig.6 Change of aerosol resuspension rate with time

图6示出气溶胶再悬浮率随时间变化的定量分析。由图6可看出，气溶胶的再悬浮是一瞬态过程，再悬浮率随时间的变化是一阶梯状的过程。

4 结论

研究结果表明，GASFLOW3.3中的再悬浮模型基本可用来模拟核电站中气溶胶的再悬浮行为，根据计算结果可得到以下结论。

1) GASFLOW与ASTEC程序的计算结果有所区别，黏性力是影响再悬浮率的一个重要因素。

2) 气溶胶的再悬浮是一瞬态过程，再悬浮率随时间的变化是一阶梯状的过程。只要已沉积的气溶胶粒子最小悬浮速度小于气流速度，气溶胶粒子就可悬浮。

3) 由计算结果可看出，GASFLOW作为一个CFD程序，计算结果与SR11试验结果的差距在20%以内，在工程范围内是可接受的。GASFLOW能较好模拟气溶胶的行为，并展示不同时刻气溶胶粒子位置的分布，可为工程应用提供参考。

4) 气溶胶再悬浮现象是安全壳内严重事故晚期的一个重要现象，从试验结果和两个程序计算的结果可看出，均有超过50%的气溶胶粒子在严重事故晚期发生再悬浮，这一比例对于源项的影响很大，在源项分析中需多加关注。同时，更加合理精确的气溶胶再悬浮模型有待进一步开发。