压力容器内超临界二氧化碳喷放特性分析及模型验证

明 杨,方华伟,刘 凯,易经纬,刘玖松,赵富龙,*,田瑞峰

(1.哈尔滨工程大学 核科学与技术学院,黑龙江 哈尔滨 150001;2.哈尔滨工程大学 黑龙江省核动力装置性能与设备重点实验室,黑龙江 哈尔滨 150001;3.中国核动力研究设计院,四川 成都 610213)

超临界二氧化碳能量转换系统具有循环效率高、系统紧凑、灵活性高、应用场景广泛等优点。近年来,与第4代反应堆结合的超临界二氧化碳布雷顿循环已成为能源行业的研究热点方向之一[1-2]。然而,超临界二氧化碳的独特物性不仅对循环设备的设计与制造带来了特殊要求,同时也为系统出现泄漏或破口事故下的安全分析带来了一系列的未知与挑战。

超临界二氧化碳流体喷放过程随着上游滞止状态的不同,可能出现单相、两相,甚至多相临界流,其喷放进程较为复杂[3]。此外,喷放过程会导致循环压力的快速降低,二氧化碳的密度、比热容、导热系数等参数将出现大幅变化,这将对系统的稳定安全运行带来较大的威胁。因此,研究超临界二氧化碳的破口喷放特性,明确循环内冷却剂的流失速率,了解压力容器内的工质状态极为重要。

在破口事故的预防与处理方面,诸多学者进行了大量的研究。早在1963年,Faltetti与Moulton[4]开展了蒸汽与水两相混合物的临界流实验研究,并建立了两相均匀平衡模型。Moody[5]分别从动量守恒与能量守恒方程出发,推导出临界流动的理论模型。之后,Henry与Fauske[6]同样基于能量守恒方程,获得了更能更适用于计算短管喷放的Henry-Fauske模型。在此之后,众多学者针对以亚临界水工质的临界流模型进行了相关的实验与模拟研究,并对已有模型进行修正与优化。

近年来,随着超临界二氧化能量转换系统的兴起,超临界流体的破口喷放现象与临界流特性研究逐步成为关注重点之一。在实验方面,德国汉堡科技大学的Gebbeken等[7]进行了超临界二氧化碳喷放的实验研究,观察到了容器内超临界二氧化碳的闪蒸现象以及喷放过程中压力曲线的剧烈弯曲,他们认为发生闪蒸的压力主要取决于初始条件。Hébrard等、Ahmad等与Guo等[8-10]分别进行了中等规模尺寸容器与大型管道内超临界二氧化碳喷放过程的相关研究,在实验进程中均观测到了饱和状态下二氧化碳的长时间稳定喷放。Fan等[11]研究了超临界二氧化碳的临界流动行为,结果表明临界流量随上游温度的增大而减小,随上游压力的增大而增大。Zhou等[12]在容积为50 L的压力容器上进行了超临界二氧化碳喷放实验研究,结果表明喷放持续的时间主要与破口直径和初始压力有关。在模拟方面,Wang等[13]建立了超临界二氧化碳压力容器的喷放模型,并采用均相模型和相分离模型分析了容器内的相分布。Min等[14]基于Henry-Fauske模型计算了超临界二氧化碳通过喷嘴的临界流动,并与实验装置得到的实验数据进行了对比。章静等[15]建立了简单容器喷放的数学物理模型,并将计算结果与实验数据进行了比较,验证了模型的合理性。李伟卿等[16]基于超临界二氧化碳临界流实验的测量结果,验证了临界流热平衡通用模型的准确性。

总体而言,超临界二氧化碳喷放过程的机理尚未明确,利用现有的实验数据得到的经验关系式适用的几何参数与工况范围较为有限,喷放过程的特性分析与流量预测等方面的工作亟待开展。因此本文基于Moody与Henry-Fauske提出的单组分混合物两相临界流分析方法,采用Modelica语言开发一套瞬态分析程序(SCRSAP-LOCA),针对超临界二氧化碳喷放过程中的临界流现象进行模拟,与现有实验数据进行对比验证,并分析容器内不同初始条件对喷放过程造成的影响。

1 计算模型

1.1 简单容器模型

本研究借鉴压水堆与超临界水堆LOCA事故的研究方法,采取简单容器模型研究超临界二氧化碳热力循环中压力容器发生破口时的喷放进程。简单容器模型基于等熵假设,认为喷放过程中流体处于热力平衡状态,忽略容器与外界的热量交换、压力梯度和流动阻力,并采用可压缩流体模型[13 , 15]。简化后的瞬态过程由质量守恒方程、能量守恒方程及状态方程共同描述[13]:

(1)

(2)

h=h(ρ,s)

(3)

式中:V为容器的体积,m3;p为容器内的压力,Pa;Gout为破口处的流量,kg/s;ρ、s、T、h分别为工质的密度、熵、温度与焓;Q为总换热量,J。根据质量守恒与能量守恒方程可以计算得到容器内工质的密度与熵,并由状态方程得到焓、温度等参数。

1.2 容器内工质的相模型

根据已有的实验研究,喷放初始阶段压力容器内的压力较高,二氧化碳处于超临界状态或液相。而伴随喷放的进行,容器内压力将快速降至二氧化碳的饱和点,气化的二氧化碳将因密度差迁移至容器顶部,容器内可以观测到明显的气液界面[7,17]。

如图1所示,在初始的单相喷放阶段,容器内的二氧化碳并未产生相变,即破口处与容器内部的二氧化碳均处于超临界态或单相液态,可采取均匀流模型:

图1 简单容器模型及喷放进程Fig.1 Simple vessel model and ejection process

hout=h0

(4)

式中,h0为容器内工质的平均焓,J/kg。

而随着喷放过程的进行,气液界面出现后的相分离效应不可忽略,破口处的工质性质应使用分相流模型表征,当气液界面高于破口高度时,认为喷放出的工质为当前压力下的饱和液体,当气液界面低于破口高度时,则认为是饱和气体[13],具体的方程如下所示:

hout=hf(p)Z>Ze

hout=hg(p)Z

(5)

式中:hf为饱和液焓;hg为饱和气焓;α为空泡份额;Z为自由液面高度,m;Ze为破口高度;Zall为容器的总高度。

1.3 破口喷放流量模型

1) 非临界流动模型

在喷放初期,对于二氧化碳超临界区的非临界喷放,由于流体的压缩性较差,该过程可视为不可压缩流体的等熵流动,一般选用Bernoulli方程计算破口处的质量流量[11]:

p0>pc,T0>Tc,x0≤0

(6)

式中:p0为上游压力;pb为背压;pc与Tc分别为二氧化碳的临界压力与温度;x0为含气率;Cd为流量系数,其定义为真实流量与理想流量的比值,取值与破口的几何形状及长径比L/d相关,在建模过程中需根据计算结果进行修正。

现有的超临界二氧化碳喷放实验中测定的Cd的取值范围为0.61～0.9[11,16,18]。本文在进行多组工况的计算与对比后,最终发现Cd选取保守值0.61时,获得的非临界流动流量计算结果与Wang等[13]的实验结果最为接近。

2) 单相过冷临界流动模型

当容器内工质的压力高于临界压力,但温度低于临界温度时,容器内工质由超临界区进入密相区,破口处为过冷临界喷放状态。本文采取Henry-Fauske模型计算该阶段的破口流量[19-20],临界流量表达式为:

p0>pt,x0>0

(7)

式中:x为含气率;v为比容,m3/kg ;cp为比定压热容,J/(kg·K);下标0代表容器内的滞止参数,下标g与f分别代表气相与液相参数;N为表征破口处发生部分相变的参数,当含气率小于0.14时,N=x/0.14,而含气率大于0.14时,N=1;γ为等熵膨胀系数,γ=cp/cv;n的定义为热力平衡多变指数,其表达式为:

(8)

临界压比η的表达式为:

(9)

式中,β、α0与αt的表达式可在文献[19-20]中查询。而当容器内的二氧化碳处于过冷状态时,临界流量表达式可进一步近似为:

p0>pt,x0≤0

(10)

在容器内滞止参数p0与x0已知的情况下,可将临界流量计算公式与临界压比进行迭代求解,具体过程如下:首先假设破口处压力等于上游滞止压力,即临界压比η=1。同时调用pt下对应的物性参数,代入式(9),可求得新的临界压比。若两次压力比数值不同,则将下游压力减掉一个小量(0.001p0),并求得新的临界压比,直到两次压力比数值相同。此时得到的压力比为临界压力比,求得的下游临界压力pt=ηp0,进而计算得到临界流量。

3) 两相临界流动模型

随着容器内工质含气率的增长,破口处逐步达到稳定的两相临界流动状态,现有较为常用的两相临界流动模型为Henry-Fauske模型与Moody模型,二者均基于等熵假设。

本文在建模时同时考虑了上述两种模型,以验证其对二氧化碳工质的适用性。其中,Henry-Fauske模型的主要方程如式(7)～(10)所述,Moody模型的原理如方程(11)所述,该模型由能量守恒出发,认为气液两相处于热力平衡状态并拥有各自的流动速度,其主要方程[21]为:

p0>pt,x0>0

(11)

式中:hfg=hg-hf;sfg=sg-sf;h0为容器内的滞止焓;S为滑速比,S=(vg/vf)1/3。上式意味着,当容器内的滞止参数已知时,质量流量Gc可由滑速比S和压力p计算。此外滑速比S表达式中的液体比容vf和气体比容vg由压力p决定。因此Moody模型最终的计算过程如下。

首先给定上游的工况条件,确定容器内两相混合物的滞止焓与滞止压力。之后给下游出口端一个初始的压力,根据动量方程原理,破口处压力应小于滞止压力,令pt略低于p0并调用此时下游压力下对应的热力参数,并代入式(11)求得对应的质量流量。之后逐步降低下游压力,求出一系列的质量流量。选取所有pt对应的流量并查询极大值,即为两相临界流流量,其对应的破口压力为临界压力。

4) 单相气临界流动模型

随着容器内工质含气率的进一步增长,当自由液面高度在破口高度以下时,喷放出的工质将由两相混合物转变为单相二氧化碳气体。认为此时的二氧化碳气体为理想气体,临界压力与临界流量由下式[15]计算:

(12)

(13)

2 瞬态代码的实现

本文的超临界二氧化碳喷放瞬态分析程序(SCRSAP-LOCA)采用Modelica语言编译,同时基于模块化思想,根据计算功能分为简单容器模块、相判断模块、破口喷放计算模块、物性查询模块等子模块。程序主体算法与计算流程如图2所示,采取有限差分方法求常微分方程,并将单相过冷临界流模型与两相临界流模型中涉及到的迭代求解编译成函数调用的方式,时间步长设置在0.01～0.1 s之间。物性数据来源于NIST RefProp软件[22],并对二氧化碳在临界点附近的物性数据采取局部加密方法,最终实现的物性查询最大相对误差在5%以内。

图2 程序主体算法与计算流程Fig.2 Program main algorithm and calculation flow

3 结果与讨论

3.1 计算结果与实验数据对比

首先将瞬态计算结果与Wang等[13]的实验数据进行了比对,实验装置主要由压力容器、加热带、保温层、管道、阀门、喷嘴、采集系统等构成。压力容器的总体积为0.05 m3,喷嘴长度为5 mm,直径为1.01 mm,喷嘴高度位于压力容器总高度的1/2处。实验的初始工况列于表1[13],程序设置的初始参数与实验保持一致。

表1 实验压力容器内的初始参数Table 1 Initial parameter in experimental pressure vessel

不同初始条件下压力容器破口流量与压力瞬变如图3所示,喷放过程可大致分为3个阶段,破口流量与容器压力呈阶梯式下降。

图3 不同初始条件下流量与压力计算结果与实验数据[13]Fig.3 Calculation results and experimental data of mass flow rate and pressure with different initial conditions[13]

第1阶段为高质量流量快速泄压阶段,流量与压力曲线的斜率最大。破口流量的范围为0.045～0.06 kg/s,容器内压力迅速降至6.5～7 MPa,同时质量流量与工质密度迅速减小。初始阶段容器中二氧化碳为超临界态或液态,其密度较高且可压缩性差,破口处于非临界流动或单相过冷临界流动状态。第1阶段持续时长为100 s左右,阶段结束的标志为压力与流量曲线斜率突然降低后出现了拐点。

第2阶段为中质量流量稳定泄压阶段,流量与压力的曲线斜率最小。破口流量的范围在0.025～0.028 kg/s之间,呈阶梯平台状缓慢下降,而容器内压力接近线性下降。此时容器内二氧化碳为两相态,破口处于两相临界流动状态,容器内工质含气率逐渐增大,但此时自由液面高度仍高于破口高度。第2阶段持续时长为350 s左右,阶段结束的标志为压力与流量曲线出现第2个拐点,并且斜率出现一定程度的增加。

第3阶段为低质量流量匀速泄压阶段,流量与压力曲线均呈凹状,斜率先增加后逐步减小。破口流量的范围为0.005～0.025 kg/s,容器内压力于1 600 s左右降至3 MPa左右。此时容器内二氧化碳仍为两相态,但随着容器内工质含气率的增大,自由液面高度低于破口高度,因此破口处工质变为单相气体,质量流量明显降低。

图4为初始条件为10.0 MPa、308 K时容器内二氧化碳的P-D(压力-密度)与P-T(压力-温度)曲线。由图4可知,喷放第1阶段开始时破口处为单相喷放,容器内二氧化碳处于超临界状态,其密度为719.51 kg/m3。当喷放第2阶段开始时,容器内二氧化碳的密度降低至609.73 kg/m3,P-D曲线与二氧化碳饱和液线相交,容器内发生相变并出现气液界面。而在喷放第3阶段内,P-D曲线逐步接近饱和气体密度曲线,P-T曲线也愈发向气相区靠近。当自由液面降至破口以下时,喷放出二氧化碳饱和气体,因此容器内工质密度的降低速率逐步降低。

图4 初始条件为10 MPa、308 K的P-D曲线与P-T曲线Fig.4 P-D and P-T curves with initial condition of 10 MPa and 308 K

总体而言,结果表明程序计算的结果与实验数据吻合较好,流量与压力的计算误差均不超过23%。在两相临界流动模型方面,Moody模型计算的破口质量流量与实验最为接近,而Henry-Fauske模型计算得到的流量相对较高,这导致流量曲线第2拐点的较早出现。程序的压力计算结果与实验数据在喷放的第1、2阶段较为接近,但在第3阶段略低于实验数据,这是由于程序忽略了容器内的压力梯度,计算所得的平均压力可能与实验中采集点得到的局部压力存在一定误差。

3.2 初始压力对喷放过程的影响

图5示出了温度为297 K,压力分别为10、15、20 MPa的不同初始条件下的主要参数计算结果。如图5a与图5c所示,在初始温度相同时,喷放第1阶段的初始流量与初始压力呈正相关,这是由于更高的初始压力使流体具有高的初始密度。同时,随着初始压力的增加,容器内压力与流体密度降低的速率也逐步增加,图5c中P-D曲线与饱和曲线的交点即为喷放第2阶段的起始点。如图5a与图5b所示,随着初始压力的增大,喷放第2阶段开始时容器内的压力变低,但两相临界喷放的质量流量几乎不随初始压力的变化而变化。此外,如图5d所示,更高的初始压力使自由液面高度下降的速率更慢,这导致喷放更晚进入第3阶段,因此两相临界喷放持续时间将更长。

图5 不同压力初始条件下主要参数计算结果Fig.5 Main parameter calculation results under different initial pressure conditions

3.3 初始温度对喷放过程的影响

图6示出了压力为15 MPa,温度分别为280、297、330、350 K的主要参数计算结果,可看到在不同的初始温度条件下,容器内二氧化碳的喷放行为将出现较大的差别。

图6 不同温度初始条件下主要参数计算结果Fig.6 Main parameter calculation results under different initial temperature conditions

当初始温度低于临界温度(304.13 K)时,容器内的二氧化碳为液态,更低的初始温度使流体具有更高的初始密度与初始喷放流量,对应喷放第2阶段的压力更低,两相临界喷放的持续时间更长。而当初始温度高于临界温度时,容器内二氧化碳处于超临界状态。当初始温度为330 K时,流体初始密度为635.51 kg/m3,喷放过程中自由液面下降的速度显著加快,这使得两相临界喷放的时间大幅缩短。值得注意的是,由图6c与6d所示,当初始温度进一步升高至350 K时,流体密度为449.20 kg/m3,超临界二氧化碳的物性更接近于气体。由图6c可知,P-D曲线由二氧化碳的超临界区进入到气相区,随后与饱和线相交,这表明容器内的二氧化碳发生了相变。液化的二氧化碳在密度差的作用下沉积至容器底部,导致了图6d中的相对液面高度曲线在200 s左右升高,这与Tian等[17]的实验测量结果相符。由于破口处将持续维持较长时间的单相气体临界喷放状态,导致流量与压力曲线变得平滑,拐点几乎消失。

4 结论

因此本文基于Modelica语言,开发了一套针对超临界二氧化碳压力容器破口喷放的瞬态分析程序(SCRSAP-LOCA),并对喷放过程中不同的流动现象进行了建模与分析,选取了多组工况的计算结果与现有实验数据进行了对比验证,分析了容器内不同初始条件对喷放过程造成的影响,得到的主要结论如下。

1) 程序的计算结果与实验数据吻合较好,流量与压力的计算误差均不超过23%。二氧化碳压力容器的喷放可大致分为3个阶段,破口流量与容器压力曲线呈阶梯式下降。

2) 程序采用的两相临界流动模型中,Moody模型计算的破口质量流量与实验值最为接近,而Henry-Fauske模型计算得到的流量与实验值相比高15%～25%左右。

3) 初始温度为297 K但初始压力增加时,第1阶段的初始喷放流量更大,第2阶段开始时容器内的压力更低,两相临界喷放持续时间更长,但破口处的两相临界喷放流量几乎不随初始压力的变化而变化。

4) 不同的初始温度使容器内二氧化碳的喷放行为具有很大的差别。当压力容器内的二氧化碳为超临界态时,初始温度的升高会使流体密度大幅降低,自由液面下降的速度显著加快,两相临界喷放的时间大幅缩短。当初始温度进一步升高,超临界二氧化碳愈发趋近于类气区域时,流量与压力曲线将变得平滑,此外喷放过程中容器内会发生相变,液化的二氧化碳会在密度差的作用下沉积至容器底部。

本文研究结果可为压力容器内超临界二氧化碳喷放特性分析提供参考,并对超临界二氧化碳动力转换系统的泄漏事故进程研究与安全分析评价提供支撑。