CN115906699A - 超快速预测水或蒸汽管道破口处临界质量流速的方法 - Google Patents

Abstract

超快速预测水或蒸汽管道破口处临界质量流速的方法，主要步骤如下：1.输入管道、破口几何参数及管道边界参数用于初始化计算；2.对破口管道划分控制体，对每个控制体建立两相质量、动量、能量守恒微分方程组，根据破口上游控制体流量、压力、焓值计算当前时刻破口的上游滞止参数；3.根据上游滞止参数中的静态质量含气率选择相应的破口计算模块；4.计算当前时刻破口处的临界压力、临界质量流速；5.根据计算的破口流量更新边界条件，转至步骤2，进行下一时刻的计算，直至计算终止时刻。本发明的方法可以计算水或蒸汽管道破口处临界质量流速，并且在计算时可以得到更精确的计算结果。

Description

超快速预测水或蒸汽管道破口处临界质量流速的方法

技术领域

本发明涉及核反应堆破口事故分析技术领域，具体涉及超快速预测水或蒸汽管道破口处临界质量流速的方法。

背景技术

破口事故的瞬态性状和物理现象已经基本明确，关于研究集中在临界流动喷放实验以及模块化建模，而大量破口事故程序往往基于单相临界流动准则进行简化，其适用范围有限，难以准确计算水/蒸汽破口。对于破口事故的计算要求，亟需开发计算精度更高的破口数值计算模块，充分考虑喷放过程中各工况下临界流模型中热不平衡及声速出现第一类间断点的特征，准确计算临界流量，克服目前大量破口程序在低欠热度、低含气率条件下计算模型误差大的缺陷。

发明内容

为了克服上述现有技术存在的问题，本发明的目的在于提供超快速预测水或蒸汽管道破口处临界质量流速的方法，以破口上游的静态质量含气率为划分标准，对于不同质量含气率范围内的上游条件，采用不同计算模型求解破口处的临界质量流速，充分体现各个范围内的热工水力特性，提高了求解的准确性，有效增加了本方法对不同问题的适应性。

为了实现上述目的，本发明采取了以下技术方案：

超快速预测水或蒸汽管道破口处临界质量流速的方法，该方法包括以下步骤：

步骤1：输入破口处管道几何参数及管道边界条件上游流动参数初始化计算条件：破口处管道几何参数具体包括管道长度、面积、破口面积，初始化计算条件包括管道内初始流量，管道进出口初始压力；

步骤2：对破口处管道划分控制体，对每个控制体建立两相质量、动量、能量守恒微分方程组，根据破口上游控制体流量、压力、焓值计算当前时刻破口的上游滞止参数；

控制体质量、动量和能量守恒微分方程组具体形式如下：

式中，ρ——两相平均密度，kg/m³；

G——质量流速，kg/(m²·s)；

t——时间，s；

z——长度，m；

P——压降，Pa；

τ——剪切力，Pa；

U——湿周，m；

A——管道流动面积，m²；

θ——管道与竖直方向所呈夹角；

α——空泡份额；

Q——传热量，W；

ρ_g——气相密度，kg/m³；

ρ_f——液相密度，kg/m³；

h_g——气相比焓，J/kg；

h_f——液相比焓，J/kg；

pv_m——比推动功，J/kg；

求解上述方程组得管道各控制体参数即压力、速度、密度，破口处管道控制体参数即为破口上游参数，破口上游滞止压力的表达式如下：

式中，p₀——破口上游滞止压力，Pa；

p_up——破口上游压力，Pa；

v_up——破口上游流体速度，m/s；

ρ_l0——破口上游流体密度，kg/m³；

步骤3：根据破口上游的质量含气率选择相应的破口计算模型；

①质量含气率＜0.000001，选用Bernoulli方程，即单相过冷区计算模型；

②0.000001≤质量含气率＜0.2，选用Henry-Fauske模型，即两相低含气率区计算模型；

③0.2≤质量含气率＜0.998，选用Moody模型，即两相高含气率区计算模型；

④质量含气率≥0.998，选用连续性方程，即单相蒸汽区计算模型；

步骤4：计算当前时刻破口处的临界压力、临界质量流速；

(一)单相过冷区计算模型

式(5)-(10)为单相过冷区计算模型，当破口上游的流体处于过冷或饱和状态时，基于Bernoulli方程，临界质量流速的表达式为：

式中，C——喷放系数，对于孔口型破口，C＝0.61；

p_c——临界压力，Pa；

v_l0——破口上游流体比容，m³/kg；

x_c——出口质量含气率，m³/kg；

v_gc——临界压力下的饱和蒸汽比容，m³/kg；

G_c——临界质量流速，kg/(m²·s)；

为了确定临界压力，引入临界压力比η，其表达式为：

η＝p_c/p₀ (6)

在破口所在管道的长径比L/D＜12时，临界压力比随长径比增加而增大；在长径比L/D≥12时，临界压力比为一固定常数0.55；

为了求解出口质量含气率x_c，基于等熵假设，引入出口热平衡含气率x_e，其表达式为：

式中，s₀——由上游滞止压力与上游滞止温度确定的比熵，J/(kg·K)；

s_g——在出口临界压力下对应的饱和汽比熵，J/(kg·K)；

s_l——在出口临界压力下对应的饱和水比熵，J/(kg·K)；

由于热不平衡效应的影响，实际的出口质量含气率x_c应该小于出口热平衡含气率x_e，因此，需要对出口热平衡含气率x_e进行修正，令N＝x_c/x_e，用来表示热不平衡的程度，N的表达式为：

N＝(0.0376L/D-0.163)exp(-0.0322ΔT_sub) (9)

式中，L/D——破口所在管道的长径比；

ΔT_sub——过冷度，K；

当L/D＜5时，忽略闪蒸过程，整个管道内都是单相液体，临界质量流速直接用下式计算：

式中，ρ_l——破口上游流体密度，kg/m³；

p₀——破口上游滞止压力，Pa；

η——临界压力比；

当上游过冷度ΔT_sub＞50K时，用上游滞止温度对应的饱和压力p_sat(t₀)作为出口临界压力，此时，临界压力比表示为η＝p_sat(t₀)/p₀；

(二)两相低含气率区计算模型

式(11)-(13)为两相低含气率区计算模型，在两相低质量含气率下，基于Henry-Fauske模型计算临界质量流速，其表达式为：

式中，n——热力平衡多变指数；

x₀——静态质量含气率；

v_g——破口处临界压力下对应的饱和蒸汽比容，m³/kg；

v_f0——上游滞止压力下对应的比容，m³/kg；

s_ge——在热平衡时，临界压力下的饱和汽比熵，J/(kg·K)；

s_fe——在热平衡时，临界压力下的饱和水比熵，J/(kg·K)；

s_g0——在上游滞止压力下对应的饱和汽比熵，J/(kg·K)；

s_f0——在上游滞止压力下对应的饱和水比熵，J/(kg·K)；

γ——定熵指数；

N——热不平衡修正项；

c_pg——气相定压比热容，J/(kg·K)；

——热平衡时饱和水比熵关于压力的导数项；

热力平衡多变指数的表达式为：

式中，x——质量含气率；

c_f——液相比热容，J/(kg·K)；

c_pg——气相定压比热容，J/(kg·K)；

在Henry-Fauske模型中，N的计算公式为：

式中，x_e的计算方法与式(8)相同；值得注意的是，若计算出x_e＞0.14，则N＝1；

在低质量含气率下，由于热不平衡的影响，临界压力比由式(7)计算；

(三)两相高含气率区计算模型

式(14)-(15)为两相高含气率区计算模型，当破口上游质量含气率较高时，两相处于热力学平衡状态，基于Moody模型计算临界质量流速，其表达式为：

式中，h₀——破口上游滞止焓，J/kg；

s₀——破口上游流体比熵，J/(kg·K)；

s_f——破口处临界压力下对应的饱和水比熵，J/(kg·K)；

s_g——破口处临界压力下对应的饱和蒸汽比熵，J/(kg·K)；

v_f——破口处临界压力下对应的饱和水比容，m³/kg；

在破口上游滞止参数一定的情况下，滑速比S是一个独立变量，其表达式为：

由于两相处于热力学平衡状态，临界压力比取0.55；

(四)单相蒸汽区计算模型

式(16)为单相蒸汽区计算模型，通过破口的蒸汽流量根据连续性方程计算，临界质量流速的表达式：

式中，v₀——上游滞止压力下比容，m³/kg；

γ——与蒸汽性质相关的参数；对于过热蒸汽，γ＝1.3；对于干饱和蒸汽，γ＝1.135；

步骤5：根据计算的破口流量更新边界所有已知条件，转至步骤2，进行下一时刻的计算，直至计算终止时刻。

与现有技术相比，本发明有如下突出特点：

1、超快速预测水或蒸汽管道破口处临界质量流速的方法，以破口上游的静态质量含气率为划分标准，对于不同质量含气率范围内的上游条件，采用不同计算模型求解破口处的临界质量流速，充分体现各个范围内的热工水力特性，提高了求解的准确性，有效增加了本方法对不同问题的适应性。

2、超快速预测水或蒸汽管道破口处临界质量流速的方法，采用公认实验关系式确定临界压力，避免了常规的迭代求解，减小了计算量，有效地提高了计算效率，实现了超快速预测。

附图说明

图1为本发明方法流程图。

具体实施方式

下面结合附图和具体实施方式对本发明作进一步详细说明：

本发明为超快速预测水或蒸汽管道破口处临界质量流速的方法，以破口上游的静态质量含气率为划分标准，对于不同质量含气率范围内的上游条件，采用不同计算模型求解破口处的临界质量流速。如图1所示，该方法具体流程包括以下步骤：

步骤1：输入破口处管道几何参数及管道边界条件上游流动参数初始化计算条件：破口处管道几何参数具体包括管道长度、面积、破口面积，初始化计算条件包括管道内初始流量，管道进出口初始压力。例如，输入管道长度25m，管道面积3.14m²,破口面积0.00314m²，之后，进行初始化。

步骤2：对破口处管道划分控制体，对每个控制体建立两相质量、动量、能量守恒微分方程组，根据破口上游控制体流量、压力、焓值计算当前时刻破口的上游滞止参数。例如，计算后的上游滞止压力为16MPa，上游滞止温度为317.35℃，质量含气率为0。

控制体质量、动量和能量守恒微分方程组具体形式如下：

式中，ρ——两相平均密度，kg/m³；

G——质量流速，kg/(m²·s)；

t——时间，s；

z——长度，m；

P——压降，Pa；

τ——剪切力，Pa；

U——湿周，m；

A——管道流动面积，m²；

θ——管道与竖直方向所呈夹角；

α——空泡份额；

Q——传热量，W；

ρ_g——气相密度，kg/m³；

ρ_f——液相密度，kg/m³；

h_g——气相比焓，J/kg；

h_f——液相比焓，J/kg；

pv_m——比推动功，J/kg；

求解上述方程组得管道各控制体参数即压力、速度、密度，破口处管道控制体参数即为破口上游参数，破口上游滞止压力的表达式如下：

式中，p₀——破口上游滞止压力，Pa；

p_up——破口上游压力，Pa；

v_up——破口上游流体速度，m/s；

ρ_l0——破口上游流体密度，kg/m³；

步骤3：根据破口上游的质量含气率选择相应的破口计算模型。根据步骤2的质量含气率0，此时应该选择单相过冷区计算模型。

①质量含气率＜0.000001，选用Bernoulli方程，即单相过冷区计算模型；

②0.000001≤质量含气率＜0.2，选用Henry-Fauske模型，即两相低含气率区计算模型；

③0.2≤质量含气率＜0.998，选用Moody模型，即两相高含气率区计算模型；

④质量含气率≥0.998，选用连续性方程，即单相蒸汽区计算模型；

步骤4：计算当前时刻破口处的临界压力、临界质量流速。根据步骤3的选择进行计算临界质量流速。

(一)单相过冷区计算模型

式(5)-(10)为单相过冷区计算模型，当破口上游的流体处于过冷或饱和状态时，基于Bernoulli方程，临界质量流速的表达式为：

式中，C——喷放系数，对于孔口型破口，C＝0.61；

p_c——临界压力，Pa；

v_l0——破口上游流体比容，m³/kg；

x_c——出口质量含气率，m³/kg；

v_gc——临界压力下的饱和蒸汽比容，m³/kg；

G_c——临界质量流速，kg/(m²·s)；

为了确定临界压力，引入临界压力比η，其表达式为：

η＝p_c/p₀ (6)

在破口所在管道的长径比L/D＜12时，临界压力比随长径比增加而增大；在长径比L/D≥12时，临界压力比为一固定常数0.55；

为了求解出口质量含气率x_c，基于等熵假设，引入出口热平衡含气率x_e，其表达式为：

式中，s₀——由上游滞止压力与上游滞止温度确定的比熵，J/(kg·K)；

s_g——在出口临界压力下对应的饱和汽比熵，J/(kg·K)；

s_l——在出口临界压力下对应的饱和水比熵，J/(kg·K)；

由于热不平衡效应的影响，实际的出口质量含气率x_c应该小于出口热平衡含气率x_e，因此，需要对出口热平衡含气率x_e进行修正，令N＝x_c/x_e，用来表示热不平衡的程度，N的表达式为：

N＝(0.0376L/D-0.163)exp(-0.0322ΔT_sub) (9)

式中，L/D——破口所在管道的长径比；

ΔT_sub——过冷度，K；

当L/D＜5时，忽略闪蒸过程，整个管道内都是单相液体，临界质量流速直接用下式计算：

式中，ρ_l——破口上游流体密度，kg/m³；

p₀——破口上游滞止压力，Pa；

η——临界压力比；

当上游过冷度ΔT_sub＞50K时，用上游滞止温度对应的饱和压力p_sat(t₀)作为出口临界压力，此时，临界压力比表示为η＝p_sat(t₀)/p₀；

(二)两相低含气率区计算模型

式(11)-(13)为两相低含气率区计算模型，在两相低质量含气率下，基于Henry-Fauske模型计算临界质量流速，其表达式为：

式中，n——热力平衡多变指数；

x₀——静态质量含气率；

v_g——破口处临界压力下对应的饱和蒸汽比容，m³/kg；

v_f0——上游滞止压力下对应的比容，m³/kg；

s_ge——在热平衡时，临界压力下的饱和汽比熵，J/(kg·K)；

s_fe——在热平衡时，临界压力下的饱和水比熵，J/(kg·K)；

s_g0——在上游滞止压力下对应的饱和汽比熵，J/(kg·K)；

s_f0——在上游滞止压力下对应的饱和水比熵，J/(kg·K)；

γ——定熵指数；

N——热不平衡修正项；

c_pg——气相定压比热容，J/(kg·K)；

——热平衡时饱和水比熵关于压力的导数项；

热力平衡多变指数的表达式为：

式中，x——质量含气率；

c_f——液相比热容，J/(kg·K)；

c_pg——气相定压比热容，J/(kg·K)；

在Henry-Fauske模型中，N的计算公式为：

式中，x_e的计算方法与式(8)相同；值得注意的是，若计算出x_e＞0.14，则N＝1；

在低质量含气率下，由于热不平衡的影响，临界压力比由式(7)计算；

(三)两相高含气率区计算模型

式(14)-(15)为两相高含气率区计算模型，当破口上游质量含气率较高时，两相处于热力学平衡状态，基于Moody模型计算临界质量流速，其表达式为：

式中，h₀——破口上游滞止焓，J/kg；

s₀——破口上游流体比熵，J/(kg·K)；

s_f——破口处临界压力下对应的饱和水比熵，J/(kg·K)；

s_g——破口处临界压力下对应的饱和蒸汽比熵，J/(kg·K)；

v_f——破口处临界压力下对应的饱和水比容，m³/kg；

在破口上游滞止参数一定的情况下，滑速比S是一个独立变量，其表达式为：

由于两相处于热力学平衡状态，临界压力比取0.55；

(四)单相蒸汽区计算模型

式(16)为单相蒸汽区计算模型，通过破口的蒸汽流量根据连续性方程计算，临界质量流速的表达式：

式中，v₀——上游滞止压力下比容，m³/kg；

γ——与蒸汽性质相关的参数；对于过热蒸汽，γ＝1.3；对于干饱和蒸汽，γ＝1.135；

步骤5：根据计算的破口流量更新边界所有已知条件，转至步骤2，进行下一时刻的计算，直至计算终止时刻。

Claims (1)

1.超快速预测水或蒸汽管道破口处临界质量流速的方法，其特征在于：包括以下步骤：

步骤1：输入破口处管道几何参数及管道边界条件上游流动参数初始化计算条件：破口处管道几何参数具体包括管道长度、面积、破口面积，初始化计算条件包括管道内初始流量，管道进出口初始压力；

步骤2：对破口处管道划分控制体，对每个控制体建立两相质量、动量、能量守恒微分方程组，根据破口上游控制体流量、压力、焓值计算当前时刻破口的上游滞止参数；

控制体质量、动量和能量守恒微分方程组具体形式如下：

式中，ρ——两相平均密度，kg/m³；

G——质量流速，kg/(m²·s)；

t——时间，s；

z——长度，m；

P——压降，Pa；

τ——剪切力，Pa；

U——湿周，m；

A——管道流动面积，m²；

θ——管道与竖直方向所呈夹角；

α——空泡份额；

Q——传热量，W；

ρ_g——气相密度，kg/m³；

ρ_f——液相密度，kg/m³；

h_g——气相比焓，J/kg；

h_f——液相比焓，J/kg；

pv_m——比推动功，J/kg；

求解上述方程组得管道各控制体参数即压力、速度、密度，破口处管道控制体参数即为破口上游参数，破口上游滞止压力的表达式如下：

式中，p₀——破口上游滞止压力，Pa；

p_up——破口上游压力，Pa；

v_up——破口上游流体速度，m/s；

ρ_l0——破口上游流体密度，kg/m³；

步骤3：根据破口上游的质量含气率选择相应的破口计算模型；

①质量含气率＜0.000001，选用Bernoulli方程，即单相过冷区计算模型；

②0.000001≤质量含气率＜0.2，选用Henry-Fauske模型，即两相低含气率区计算模型；

③0.2≤质量含气率＜0.998，选用Moody模型，即两相高含气率区计算模型；

④质量含气率≥0.998，选用连续性方程，即单相蒸汽区计算模型；

步骤4：计算当前时刻破口处的临界压力、临界质量流速；

(一)单相过冷区计算模型

式(5)-(10)为单相过冷区计算模型，当破口上游的流体处于过冷或饱和状态时，基于Bernoulli方程，临界质量流速的表达式为：

式中，C——喷放系数，对于孔口型破口，C＝0.61；

p_c——临界压力，Pa；

v_l0——破口上游流体比容，m³/kg；

x_c——出口质量含气率，m³/kg；

v_gc——临界压力下的饱和蒸汽比容，m³/kg；

G_c——临界质量流速，kg/(m²·s)；

为了确定临界压力，引入临界压力比η，其表达式为：

η＝p_c/p₀ (6)

在破口所在管道的长径比L/D＜12时，临界压力比随长径比增加而增大；在长径比L/D≥12时，临界压力比为一固定常数0.55；

为了求解出口质量含气率x_c，基于等熵假设，引入出口热平衡含气率x_e，其表达式为：

式中，s₀——由上游滞止压力与上游滞止温度确定的比熵，J/(kg·K)；

s_g——在出口临界压力下对应的饱和汽比熵，J/(kg·K)；

s_l——在出口临界压力下对应的饱和水比熵，J/(kg·K)；

由于热不平衡效应的影响，实际的出口质量含气率x_c应该小于出口热平衡含气率x_e，因此，需要对出口热平衡含气率x_e进行修正，令N＝x_c/x_e，用来表示热不平衡的程度，N的表达式为：

N＝(0.0376L/D-0.163)exp(-0.0322ΔT_sub) (9)

式中，L/D——破口所在管道的长径比；

ΔT_sub——过冷度，K；

当L/D＜5时，忽略闪蒸过程，整个管道内都是单相液体，临界质量流速直接用下式计算：

式中，ρ_l——破口上游流体密度，kg/m³；

p₀——破口上游滞止压力，Pa；

η——临界压力比；

当上游过冷度ΔT_sub＞50K时，用上游滞止温度对应的饱和压力p_sat(t₀)作为出口临界压力，此时，临界压力比表示为η＝p_sat(t₀)/p₀；

(二)两相低含气率区计算模型

式(11)-(13)为两相低含气率区计算模型，在两相低质量含气率下，基于Henry-Fauske模型计算临界质量流速，其表达式为：

式中，n——热力平衡多变指数；

x₀——静态质量含气率；

v_g——破口处临界压力下对应的饱和蒸汽比容，m³/kg；

v_f0——上游滞止压力下对应的比容，m³/kg；

s_ge——在热平衡时，临界压力下的饱和汽比熵，J/(kg·K)；

s_fe——在热平衡时，临界压力下的饱和水比熵，J/(kg·K)；

s_g0——在上游滞止压力下对应的饱和汽比熵，J/(kg·K)；

s_f0——在上游滞止压力下对应的饱和水比熵，J/(kg·K)；

γ——定熵指数；

N——热不平衡修正项；

c_pg——气相定压比热容，J/(kg·K)；

——热平衡时饱和水比熵关于压力的导数项；

热力平衡多变指数的表达式为：

式中，x——质量含气率；

c_f——液相比热容，J/(kg·K)；

c_pg——气相定压比热容，J/(kg·K)；

在Henry-Fauske模型中，N的计算公式为：

式中，x_e的计算方法与式(8)相同；值得注意的是，若计算出x_e＞0.14，则N＝1；

在低质量含气率下，由于热不平衡的影响，临界压力比由式(7)计算；

(三)两相高含气率区计算模型

式(14)-(15)为两相高含气率区计算模型，当破口上游质量含气率较高时，两相处于热力学平衡状态，基于Moody模型计算临界质量流速，其表达式为：

式中，h₀——破口上游滞止焓，J/kg；

s₀——破口上游流体比熵，J/(kg·K)；

s_f——破口处临界压力下对应的饱和水比熵，J/(kg·K)；

s_g——破口处临界压力下对应的饱和蒸汽比熵，J/(kg·K)；

v_f——破口处临界压力下对应的饱和水比容，m³/kg；

在破口上游滞止参数一定的情况下，滑速比S是一个独立变量，其表达式为：

由于两相处于热力学平衡状态，临界压力比取0.55；

(四)单相蒸汽区计算模型

式(16)为单相蒸汽区计算模型，通过破口的蒸汽流量根据连续性方程计算，临界质量流速的表达式：

式中，v₀——上游滞止压力下比容，m³/kg；

γ——与蒸汽性质相关的参数；对于过热蒸汽，γ＝1.3；对于干饱和蒸汽，γ＝1.135；

步骤5：根据计算的破口流量更新边界所有已知条件，转至步骤2，进行下一时刻的计算，直至计算终止时刻。

Images