LOCA泄压工况下氮气的迁移析出建模研究

2023-10-27刘碧帆黄广源尹俊连王德忠赵剑刚鄢梦琪

原子能科学技术 2023年10期

刘碧帆,黄广源,吕欣,尹俊连,*,王德忠,赵剑刚,鄢梦琪

(1.上海交通大学机械与动力工程学院,上海 200240;2.中广核研究院有限公司,广东深圳 518026)

氮气加压具有非能动调节的优点,广泛应用于反应堆中。非能动氮气安注箱是VVER[1]、AP1000[2]等第3代核电堆型的重要专设安全设施。失水事故(LOCA)中,当一回路压力降至安注箱预设压力时,氮气安注箱可以非能动地向反应堆内注入冷却水,确保反应堆安全。氮气稳压器在先进一体化小型压水堆系统中具有极大应用前景[3-4],其利用充入氮气的可压缩性非能动地调节一回路的压力,具有结构简单可靠、运行能耗低、可快速启动等优点。但氮气加压带来了新的问题与挑战。氮气加压的冷却剂中的溶解氮为饱和状态。LOCA中,随着压力的持续下降,溶解氮进入过饱和状态,以微气泡的形式大量析出积聚,导致传热恶化,阻滞安注流量,影响反应堆安全[5-6]。为了保证反应堆安全,防止溶解氮的析出,需要对溶解氮的迁移和析出过程进行建模。

本文针对氮气析出的具体过程,建立氮气壁面析出模型,解释析出源项的物理意义。对氮气的溶解平衡、对流扩散、壁面析出、气泡上浮进行数学建模,基于MATLAB编写一维集总参数氮气迁移析出求解程序,对Sarrette实验中SBLOCA泄压工况下溶解氮迁移析出特性进行研究,获得溶解氮浓度和含气率变化。本研究建立的预测方法可为其他含氮系统内氮气析出的预测提供参考。

1 数学物理建模

本研究对氮气迁移析出过程中的各物理过程进行了数学物理建模。LOCA泄压工况中,氮气迁移析出主要经历了4个物理过程:溶解平衡、对流扩散、壁面析出和气泡上浮(图1)。由于溶解平衡的影响,液面处的溶解氮浓度始终为对应压力的平衡浓度。随着压力逐渐降低,液面处的溶解氮浓度也逐渐降低。由于浓度梯度和液相脉动,溶解氮在液相中扩散。液相的运动对溶解氮起到对流输运作用。在壁面附近,溶解氮向气核内析出氮气。氮气泡逐渐生长,最终脱离壁面上浮。

图1 泄压实验中的物理过程Fig.1 Physical process in pressure relief experiment

1.1 氮气质量守恒方程

考虑如图2所示的控制体V,在t时刻长度为δx、横截面积为A的控制体V中氮气的质量为ρC(1-α)Aδx。图2中:ρ为控制体中水的密度;C为溶解氮的浓度;α为空泡份额;Aδx为控制体的体积。在一维情况下,溶解氮只能在左右边界上通过对流作用和扩散效应进出控制体,且控制体内氮气满足如下质量守恒:

图2 控制体中的氮气质量守恒Fig.2 Nitrogen mass conservation control volume

Aδx((1-αt+δt)ρt+δtCt+δt-(1-αt)ρtCt)=

(1)

式中:上标表示时间,下标L、R分别表示左、右边界;u为液相速度;D为分子扩散系数Dm与涡流扩散系数Deddy之和。等号左侧表示控制体内氮气质量随时间的变化率,等号右侧第1项、第2项表示单位时间内通过对流和扩散作用分别在左边界进入及在右边界离开控制体的氮气质量,第3项表示由于析出源项减少的氮气质量。

当时间微元δt与控制体尺寸δx足够小,且AL=AR时,得到一维情况下描述溶解氮质量守恒方程:

(2)

式中:等号左侧表示单位长度的控制体内溶解氮在单位时间内的质量变化;等号右侧第1项为对流引起的质量变化,其物理意义为溶解氮由于液相运动而跟随运动的溶解氮质量,第2项为扩散引起的质量变化,第3项为析出引起的质量变化。

同理可得描述氮气对流和析出的质量守恒方程为:

(3)

式中:ρg为气相密度;ug为气相速度;等号左侧表示单位长度的控制体内氮气在单位时间内的质量变化;等号右侧第1项为对流引起的质量变化,第2项为析出引起的质量变化。

1.2 溶解模型

氮气溶解的动态平衡主要存在于气液界面。Jirka[8]通过实验测量发现,在十分靠近界面的区域中溶解气体的浓度满足亨利定律:

pB=HBCB

(4)

式中:pB为平衡压力;HB为亨利系数;CB为饱和浓度。

溶解度数值由Baranenko等[9]的数据插值得到。

1.3 扩散模型

D为分子扩散系数Dm与涡流扩散系数Deddy之和:

D=Dm+Deddy

(5)

分子扩散系数表达式为:

(6)

式中:kB为玻尔兹曼常数;T为热力学温度;nSE为Stokes-Einstein常数;η为水的黏度;a为氮气分子的动力学半径。

涡流扩散系数的表达式为:

(7)

式中:νt为湍流黏度;Sct为湍流施密特数,Sct取值1.36[10]。

1.4 壁面析出模型

1) 析出的主要机理

在泄压过程中压力变化较大,溶解氮气达到过饱和状态,因此需要额外考虑氮气的析出及其对氮气迁移的影响。初步估算表明,泄压过程中最大过饱和度不会超过10,远不足以克服在连续的液相介质中形成稳定的气液界面所需的能量壁垒。因此本文仅考虑第Ⅳ类非均相核化过程[11],即认为系统固体壁面上已经稳定存在着大于临界尺寸的微气核,当系统处于过饱和状态后,液相中溶解的氮气通过扩散作用进入气核并使其长大直至脱离壁面。下面针对这一过程进行物理建模。

2) 气泡成核条件

固体壁面上稳定存在的微气核半径往往在微纳米量级,其形貌可以近似简化为球帽。下面推导微气核生长所需的过饱和度。

根据杨-拉普拉斯公式,气泡内外压差为:

(8)

式中:σ为对应温度压力下的表面张力系数;R为气泡半径;pg为气泡内气体压力;p0为液相压力(系统压力)。根据亨利定律,此时气液界面上的平衡浓度Cg为:

(9)

式中:Hg为对应温度压力下的亨利系数;C0=p0/Hg,为系统当地温度压力下对应的平衡氮气浓度,即该处的氮气溶解度。

此时过饱和度ζ为:

(10)

初始气核半径一般在微米量级[12-14]。将SBLOCA中压安注工况(p0=5.5 MPa,σ=0.070 78 N/m)代入式(10)可得:在初始气核半径为100 μm时,过饱和度仅为2.57×10-4,当初始气核半径为1 μm时,临界过饱和度为2.57×10-2。可见,在壁面微气泡尺度下,过饱和析出所需的过饱和度非常小,不同初始气核半径对临界过饱和度的绝对值影响可忽略,可近似认为当系统一旦达到过饱和状态后,壁面上就会发生核化析出。

3) 气泡生长模型

气泡生长主要受扩散作用控制,用无限大过饱和流体中气泡的自由长大过程近似地描述[12]:

(11)

(12)

根据理想气体状态方程,气泡内的气相密度ρg可表示为:

(13)

式中,ρg0为p0压力下对应的气相密度。

将式(8)、(9)、(13)代入式(12)得:

(14)

该式具有R′=f(t,R)的形式,可通过显式的龙格-库塔方法对R进行数值求解。

为了分析半径R(t)随时间变化的趋势,可考虑将(C∞-Cg)/ρg视为常数,获得形式简单的解析解。如图3所示,数值计算结果表明,气泡生长周期尺度(数十毫秒)远小于析出过程的时间尺度(数十秒)。过饱和度在1个气泡周期内几乎没有变化。气泡尺寸R在10-6量级,2σ/(p0R)≪1,ρg=ρg0(1+2σ/p0R)≈ρg0。因此可将(C∞-Cg)/ρg视为常数,从而获得式(12)的解析解:

图3 SBLOCA工况中不同深度的气泡生长曲线Fig.3 Bubble growth curves at different depths in SBLOCA

(15)

4) 气泡脱落模型

考虑Jones等[11]和Klausner等[14]提出的力平衡模型。由于事故工况下回路的剪切效应较弱,认为析出气泡的脱离主要是浮升力引起的。

如图4所示,在竖直固体壁面上,气泡脱离时浮升力与三相接触线上的表面张力在垂直方向上达到平衡:

图4 气泡在竖直壁面上的受力平衡Fig.4 Force balance of a bubble on vertical wall

Fb=Fσ

(16)

(17)

Fb=(ρl-ρg)gV

(18)

式中:θ1、θ2分别为气液固三相接触线的前进接触角和后退接触角;dw为三相接触线的直径。

当核化点上的气泡析出长大达到临界直径脱落后,将其重置为初始尺寸,重新开始新一轮的析出长大以及脱落。

5) 析出质量源项

析出氮气质量源项S的物理意义为单位时间内单位体积的控制体中由于气泡在壁面析出而变化的氮气质量,用单位体积内的气泡数目乘以气泡体积的变化率表示为:

(19)

若忽略表面张力系数σ随时间的变化,有:

(20)

将式(14)代入可得:

(21)

式中,R为该计算时刻的析出气泡半径。

1.5 气泡上浮模型

气泡析出之后在水中上浮,受到浮力、曳力、虚拟质量力、壁面润滑力等力的作用。由于阻力的影响,气泡进入主流后作加速度逐渐减小的变加速运动,最后以某一特定速度匀速上升。Kamp和Kraume等[15]根据力平衡推导出气泡运动方程:

(22)

式中:u为气泡运动速度;ρc为连续相密度;ρd为离散相(气泡)密度;d为气泡直径;CVM为虚拟质量系数,取值为0.785[16];CD为曳力系数[17],可表示为:

(23)

2 数值求解与程序编制

2.1 数值方法

质量守恒方程的对流项采用一阶迎风格式离散。质量守恒方程的扩散项采用两次中心差分格式进行离散。质量守恒方程的析出质量源项由式(19)、(21)和气泡半径R代数求解得到,R通过式(14)由显式的龙格-库塔方法数值求解得到。式(14)是非线性的常微分方程,具有R′=f(t,R)的形式,通过显式的龙格-库塔方法对R进行数值求解,累积误差为Δt4阶,在时间步长较小时计算准确度很高。

2.2 程序编制

根据氮气的数学物理模型和相应的数值方法,采用MATLAB编写了一维集总参数氮气迁移析出求解程序,计算流程如图5所示。该程序考虑了氮气的的溶解平衡、对流扩散、壁面析出和气泡上浮过程,可获得一维管道中每一时刻每一深度的溶解氮浓度和空泡份额。输入初始条件和边界条件后,程序开始计算。根据氮气壁面析出模型,求解当前时间步气泡半径,若气泡半径小于气泡脱离半径,则基于当前气泡半径和质量守恒方程求解溶解氮浓度和空泡份额。若气泡半径大于等于气泡脱离半径,则判定该气泡脱离,将气泡半径重置为初始半径之后再求解溶解氮浓度和空泡份额。到达设定的仿真时间后停止计算。

图5 计算流程Fig.5 Calculation process

3 验证工况

选择Sarrett等[7]的SBLOCA泄压实验进行验证。Sarrett的实验在1个竖直放置的密闭圆管中进行,圆管的顶部与蓄压箱、泄压阀相连。压力测点和热电偶测点布置如图6所示。圆管高8.8 m、内径139.7 mm、水深5 780 mm,圆管上部分充有初始压力为6 MPa、温度为27 ℃的氮气。水的初始温度为33 ℃。实验在静水工况下进行,可认为水的速度为0。

图6 泄压实验装置示意图[7] Fig.6 Pressure relief experimental section geometry[7]

Loviisa核电站中氮气蓄压安注箱的预设压力为5.5 MPa,注入上排管路时的水温为40 ℃。为模拟该堆型SBLOCA下安注管道的泄压过程,实验过程中水的温度维持在32 ℃,在700～787 s,气体压力从5.5 MPa快速降至2.0 MPa(图7)。

图7 泄压实验段压力变化曲线[7]Fig.7 Pressure profile of experimental section under pressure relief condition[7]

4 验证计算

4.1 参数选择

如前文所述,对氮气迁移析出计算影响较大且目前尚不确定的参数包括:气泡脱离直径、气相运动速度、计算的时间步长Δt、初始气核尺寸。

根据前述气泡脱离的力平衡方程,初步计算得析出气泡脱离直径在10-1mm量级,与SBLOCA泄压实验测得气泡脱离直径(0.1 mm)基本吻合。取0.1 mm作为气泡脱离直径。

根据前述颗粒运动方程,初步计算得气泡在主流运动的终速度在10-2m/s量级,与SBLOCA泄压实验测得气相平均速度(0.01 m/s)基本吻合。对于10-1mm量级的析出气泡,计算得到变加速运动阶段在0.1 s以内,相比于气泡从产生到到达液面所需的时间(数百秒)而言十分短暂。因此可忽略变加速阶段,近似地认为气相以终速度作匀速运动,取0.01 m/s作为气相运动速度。

如图8所示,在满足显式计算库朗特数的前提下,随着计算时间步长的减小,氮气析出变化过程逐渐收敛,可认为,在计算的时间步长小于0.01 s后氮气的析出过程计算是准确的,能够获得时间无关解。

图8 计算时间步长对氮气析出的影响Fig.8 Effect of time step on nitrogen release

在经过一段时间t之后,某核化点析出的氮气总质量m可表示为:

(24)

式中:T为气泡生长周期;ρg为这段时间内的气泡中的气体平均密度。

(25)

假设气泡半径均在微米量级,满足1 mm>Rd>R1>R2>1 μm,且Rd=k1R1,R1=k2R2,其中k1>1,k2>1。

将R1、R2代入上式得:

(26)

由k1>1,k2>1,且:

(27)

(28)

可得到不等式关系:

(29)

将工况参数(p0=5.5 MPa,σ=0.070 78 N/m,1 mm>Rd>R1>R2>1 μm)代入,当初始气核半径比脱离半径小1倍,满足k1>2、k2>1时,代入上式可知,1.026>m2/m1>0.731,析出氮气质量最多只相差26.4%。当初始气核半径比脱离半径小1个量级,满足k1>10、k2>1时,1.026>m2/m1>0.965,析出氮气质量最多只相差3.5%。可见,由于初始气核半径一般都在微米量级[13,18-19],其具体数值的选取并不会显著影响氮气析出质量。出于保守设计考虑,取初始气核半径R0=1 μm。

4.2 初始条件与边界条件

一维计算液面边界节点处给定Dirichlet边界条件以描述氮气的溶解平衡,即C=CB。

一维计算圆管底部边界节点的气相控制方程离散时,不考虑流入控制体的氮气质量,只考虑流出控制体的氮气质量和析出的氮气质量,即满足:

(30)

根据SBLOCA泄压实验的实验数据,初始时刻的空泡份额为0.3%。初始时刻的溶解氮浓度为饱和浓度。

4.3 计算结果

程序计算结果如下。溶解氮平均浓度,平均空泡份额均与Sarrett实验值符合良好(图9)。如图9所示,不考虑氮气析出,仅考虑扩散作用的情况下,溶解氮的浓度几乎没有变化。

图9 SBLOCA泄压工况溶解氮浓度和含气率计算结果Fig.9 Dissolved nitrogen concentration and void fraction under pressure relief condition during SBLOCA

如图10所示,对于圆管中各位置,随着压力的下降,浓度逐渐下降,浓度的下降速率逐渐增加。平均浓度从815 ppm下降至约656 ppm,下降了19%。各时刻,圆管中浓度随深度的增加而线性增加,但变化幅度较小,不超过1%。

图10 SBLOCA泄压工况浓度-时间、浓度-深度数值模拟Fig.10 Numerical simulation results of concentration-time and concentration-depth under pressure relief condition during SBLOCA

如图11所示,除去圆管底部,圆管中各位置随着压力的下降,含气率逐渐上升,含气率的上升速率逐渐增加。含气率从约0.003升至约0.013。管道中析出了1.152 L的氮气。圆管底部的含气率-时间曲线与其他位置表现出不同的变化趋势。这是因为圆管底部流入的氮气质量更少。在降压初期,由于流出的氮气质量大于析出的氮气质量,圆管底部的含气率逐渐降低;在降压后期,由于过饱和度逐渐增加,析出速率增加,流出的氮气质量小于析出的氮气质量,圆管底部的含气率逐渐上升。从图11还可看出,各时刻,圆管中含气率随深度的增加而减小。含气率-深度曲线可分为线性下降段和加速下降段。由于圆管底部控制体节点的含气率求解采用了特殊的边界条件(不考虑对流流入项),该边界条件的影响会随着时间逐渐从上游向下游传递。线性下降段是未受影响的区域,加速下降段则是已经受到影响的区域。