CN112613240A - 一种用于严重事故下安全壳内流动分析的计算方法 - Google Patents

Abstract

本发明涉及一种用于严重事故下安全壳内流动分析的计算方法，包括如下步骤：1、获取安全壳流动分析计算所需参数并初始化；2、根据参数建立流动方程；3、计算重力压头；4、进行压力线性化；5、利用准Newton迭代法求解动量方程；6、根据预测的新时刻速度求解流道空泡数；7、根据空泡数修正预测速度以得到新时刻速度；8、根据新时刻速度确定流道上、下游；9、根据新时刻速度求解质量、能量增量。本发明为针对严重事故下安全壳内流动分析问题，提出了一种基于集总参数法的动量方程计算方法，可准确并且快速的求解安全壳内流动情况，进而得到安全壳内由流动引起的质量能量变化情况，为安全壳内其他严重事故现象的分析打下基础。

Description

一种用于严重事故下安全壳内流动分析的计算方法

技术领域

本发明涉及核安全分析领域，具体涉及一种用于严重事故下安全壳内流动分析的计算方法。

背景技术

日本福岛核电事故引发了对反应堆严重事故新的关注和研究，其中对严重事故现象的分析显得更加重要，对提高电厂安全性有着重要意义。在严重事故现象中，安全壳热工水力现象是所有其它事故现象的根据，是安全壳安全分析的基础。目前常用的核电厂热工水力分析程序主要为集总参数模型。在集总参数模型中，连续的空间分割成互不重叠的各个控制体，控制体的两相流体满足质量、能量守恒方程和状态方程，控制体之间通过流道连接起来，流道中的流动遵循动量守恒定律。质量、动量、能量守恒方程组成控制热工水力学行为的非线性的常微分方程组，用于描述两相流体的质量、动量和能量等状态变量的演化。其中的动量方程计算，涉及到众多状态量和参数的相互作用，是方程组求解的关键。运用恰当的数值计算技巧求解动量方程，是当前众多严重事故分析软件的一个重要任务和主要差别所在。

目前常用的安全壳动量方程求解算法，有些考虑了过多状态变量和热工参数的影响，有些未考虑动量方程中的非定常项，因而计算效率或者计算精度受到影响。

发明内容

本发明的目的在于针对现有技术的缺陷，提供一种基于集总参数法的动量方程计算分析方法，用于严重事故下安全壳内流动分析计算。本发明提出的动量方程计算方法可准确并且快速的求解流道中的动量方程，从而得到准确的安全壳内由流动引起的质量能量变化情况，为安全壳内其他严重事故现象的分析打下基础。

为了实现上述目的，本发明提供的技术方案如下：

一种用于严重事故下安全壳内流动分析的计算方法，包括如下步骤：

步骤1、获取安全壳流动分析计算所需参数并初始化，即获取待分析安全壳流动分析算例中控制体、流道的几何参数、状态参数和连接关系等相关参数；

步骤2、建立流动状态方程，流动现象经过模型化后最终可以归结为一组集总参数形式的、非线性的、耦合的一阶常微分方程，构成了一个动力学系统；

步骤3、计算重力压头，连接控制体i和k的流道j中作用于相

的压力差为

其中P_i和P_k是相应控制体积中的热力学压力，

包含控制体内和沿着流道的所有重力头的项；

在不考虑流动时，重力压头

包含三部分，第一部分是控制体i中水池面高程z_w,i与连接处各相高程

之间的水头差，第二部分为控制体k中水池面高程z_w,k与连接处各相高程

之间的水头差，第三部分是各相

沿着流道的重力压头，纯重力压头定义为以上三部分之和；

步骤4、进行压力线性化，将动量方程中压力项线性化，将作为新时刻动力学的非线性函数的压力在旧时刻的压力的基础上展开为新时刻动力学变量的线性函数；

步骤5、根据步骤3、4，使用压力线性化表示新时刻压力，并将重力压头代入，利用准Newton迭代法求解动量方程；

步骤6、根据预测的新时刻速度求解流道空泡数；

步骤7、根据空泡数修正预测速度以得到新时刻速度；

步骤8、根据新时刻速度确定流道上、下游；

步骤9、根据新时刻速度求解质量、能量增量。

进一步，如上所述的用于严重事故下安全壳内流动分析的计算方法，步骤2中所述动力学系统相应的常微分方程如下：

式中：

j——第j个流道；

——某一相流体，可为气相或液相；

d——流道的上游；

r——流道的下游；

——流道j上游控制体相

密度；

——流道j相

速度；

t——时间；

L_j——流道j中非定常流动的长度，一般取为流道长度；

——流道上下游控制体中的热力学压强；

——流道j两端液相和气相各自的重力头差；

ΔP_j——流体机械产生的泵头；

K_j——无量纲的流道阻力系数。

进一步，如上所述的用于严重事故下安全壳内流动分析的计算方法，步骤3中所述重力压头

包含的三部分中，第一部分是控制体i中水池面高程z_w,i与连接处各相高程

之间的水头差

计算公式如下：

式中：

——控制体i中水池面高程z_p,i与连接处各相高程

之间的水头差；

ρ_w,i——控制体i中液相密度；

ρ_g,i——控制体i中气相密度；

g——重力加速度；

Z_w,i——控制体i液面高程；

——控制体i连接处各相高程；

第二部分为控制体k中水池面高程z_w,k与连接处各相高程

之间的水头差，计算公式如下：

式中：

——控制体k中水池面高程z_w,k与连接处各相高程

之间的水头差；

ρ_w,k——控制体k中液相密度；

ρ_g,k——控制体k中气相密度；

g——重力加速度；

Z_w,k——控制体k液面高程；

——控制体k连接处各相高程；

第三部分是各相

沿着流道的重力压头，其计算如下：

式中：

——各相

沿着流道的重力压头；

——流道中相

的密度，取为流道上游与下游控制体内该相密度的最大值；

g——重力加速度；

——控制体i连接处各相高程；

——控制体k连接处各相高程；

综上，纯重力压头定义为这三部分之和，如下所示：

式中：

——控制体内和沿着流道的所有重力压头的项；

——控制体i中水池面高程z_w,i与连接处各相高程

之间的水头差；

——控制体k中水池面高程z_w,k与连接处各相高程

之间的水头差；

——各相

沿着流道的重力压头；

将重力压头

进行线性化，考虑一个时间步长内流动引起的流道两端的重力压头的变化；对于气体而言，流动引起的重力压头的变化相较于压力变化是可以忽略的；水的重力压头的变化则相对不是很小，忽略流道中因为水的流动引起的重力压头变化，则流道两端水的重力压头变化为流道两端控制体水的重力压头变化之和，计算公式如下：

式中：

Δ(ρgΔz)_j,w——一个时间步长内因水的流动引起的流道两端水的重力压头的变化；

Δm_i、Δm_k——控制体i、k内水的质量变化，以控制体i为例，其等于

其中Δt为时间步长，σ_ij表示控制体与流道j的关系，A_j为流道面积，α_j,w为流道中液相空泡数，

流道上游控制体液相密度，v_j,w流道液相速度；

g——重力加速度；

S_i、S_k——控制体i、k的底面积；

i——控制体i；

k——控制体k；

Z_w,i、Z_w,k——控制体i、k中水池面高程；

Z_J,w,i、Z_J,w,k——控制体i、k连接处液相高程；

进一步，假定控制体内水的密度在一个时间步长内不变，以控制体i为例，一个时间步长内控制体中水的重力压头变化可由下式计算：

式中：

g——重力加速度；

Δm——控制体内水的质量变化；

S——控制体底面积；

Δt——时间步长；

σ_ij——控制体与流道j的关系，可以取为1、-1、0，分别表示流进、流出控制体或者与控制体不相连；

A_j——流道面积；

α_j,w——流道中液相空泡数；

——流道上游控制体液相密度；

v_j,w——流道液相速度；

进一步，如上所述的用于严重事故下安全壳内流动分析的计算方法，步骤4中进行压力线性化，其中新时刻压力可由下式计算

式中：

P^new——新时刻压力；

P^old——旧时刻压力；

——压力对时间的微分；

Δt——时间步长；

公式(8)中的压力对时间的微分可使用公式(9)计算

式中：

t——时间；

P——总压力；

a——等于

其中R为气体常数，T_g为气相温度，V_g为气相体积，z_s为蒸汽压缩因子，M_w为水的摩尔质量，ρ_s为蒸汽的密度，ρ_g为气相不可凝气体密度，x_k为气相不可凝气体组分k的质量份额，M_k为气相不可凝气体组分k摩尔质量；

b——等于

其中R为气体常数，z_s为蒸汽压缩因子，M_w为水的摩尔质量，ρ_s为蒸汽的密度，ρ_g为气相不可凝气体密度，x_k为气相不可凝气体组分k的质量份额，M_k为气相不可凝气体组分k摩尔质量；

R——理想气体常数；

m_k——气相不可凝气体组分k质量；

M_k——气相不可凝气体组分k摩尔质量；

ρ_k——气相不可凝气体组分k密度；

m_s——蒸汽的质量；

ρ₀——水密度的参考状态，1000kg/m³；

a_s——蒸汽中的声速；

c_vg——气相气体比热；

T_g——气相温度；

V_g——气相体积；

u_g——气相比内能；

m_w——水的质量；

ρ_w——水的密度；

c_vw——水的比热；

T_w——水的温度；

u_w——水的比内能；

β——水的膨胀系数；

σ_ij——控制体与流道的关系，可以取为1、-1、0，分别表示流进、流出控制体或者与控制体不相连；

A_j——流道面积；

α_j,g——流道气相的体积分数；

α_j,w——流道液相的体积分数；

——流道上游控制体气相的密度；

——流道上游控制体液相的密度；

v_j,g——流道中气相速度；

v_j,w——流道中液相速度；

——流道上游控制体气相的比焓；

——流道上游控制体液相的比焓；

δm_g——总气体质量源变化率；

δm_w——液相水质量源变化率；

δm_s——水蒸汽质量源变化率；

δm_k——气相不可凝气体k的质量源变化率；

δU_g——气相能量源项变化率；

δU_w——液相水能量源项变化率。

进一步，如上所述的用于严重事故下安全壳内流动分析的计算方法，步骤5中将公式(7)和公式(8)代入步骤1中的公式(1)，分别可得到气相和液相的方程，见公式(10)和公式(11)

式中：

t——时间；

Δt——时间步长；

P——总压力；

a——等于

其中R为气体常数，T_g为气相温度，V_g为气相体积，z_s为蒸汽压缩因子，M_w为水的摩尔质量，ρ_s为蒸汽的密度，ρ_g为气相不可凝气体密度，x_k为气相不可凝气体组分k的质量份额，M_k为气相不可凝气体组分k摩尔质量；

b——等于

其中R为气体常数，z_s为蒸汽压缩因子，M_w为水的摩尔质量，ρ_s为蒸汽的密度，ρ_g为气相不可凝气体密度，x_k为气相不可凝气体组分k的质量份额，M_k为气相不可凝气体组分k摩尔质量；

(ρgΔz)_j,g——流道中气相产生的重力头；

(ρgΔz)_j,w——流道中液相产生的重力头；

Δ(ρgΔz)_j,w——由于流动而在一个时间步长内产生的液相重力头变化；

——流道上游控制体旧时刻压力；

——流道下游控制体旧时刻压力；

S——控制体底面积；

ΔP_j——流体机械产生的泵头；

R——理想气体常数；

m_k——气相不可凝气体组分k质量；

M_k——气相不可凝气体组分k摩尔质量；

ρ_k——气相不可凝气体组分k密度；

m_s——蒸汽的质量；

ρ₀——水密度的参考状态，1000kg/m³；

a_s——蒸汽中的声速；

c_vg——气相气体比热；

T_g——气相温度；

V_g——气相体积；

u_g——气相比内能；

m_w——水的质量；

ρ_w——水的密度；

c_vw——水的比热；

T_w——水的温度；

u_w——水的比内能；

β——水的膨胀系数；

σ_ij——控制体与流道的关系，可以取为1、-1、0，分别表示流进、流出控制体或者与控制体不相连；

A_j——流道面积；

L_j——流道中非定常流动的长度，一般取为流道长度；

K_j——流道阻力系数；

α_j,g——流道气相的体积分数；

α_j,w——流道液相的体积分数；

——流道上游控制体气相的密度；

——流道上游控制体水的密度；

v_j,g——流道中气相速度；

v_j,w——流道中液相速度；

——流道上游控制体气相的比焓；

——流道上游控制体液相的比焓；

方程(10)和(11)是非线性常微分方程，它们构成了一个常微分方程系统，将常微分方程系统简记为公式(12)

式中f(v)是v的非线性函数，v为待求解速度矢量；

进一步，采用半显半隐式时间离散动力系统，公式(12)可离散为公式 (13)

式中：

Δt——时间步长；

vⁿ⁺¹——新时刻速度矢量；

vⁿ——当前时刻速度矢量；

v^n-1——旧时刻速度矢量；

φ——比例系数，当φ≠0时，时间离散格式具有二阶精度；

公式(13)是一个非定常、非线性方程组，可采用迭代法求解，因而可写为公式(14)

F(vⁿ⁺¹)＝(1+φ)(vⁿ⁺¹-vⁿ)-φ(vⁿ-v^n-1)-Δtf(vⁿ⁺¹)＝0 公式(14)

F(vⁿ⁺¹)——速度vⁿ⁺¹的非线性函数；

Δt——时间步长；

vⁿ⁺¹——新时刻速度矢量；

vⁿ——当前时刻速度矢量；

v^n-1——旧时刻速度矢量；

φ——比例系数，当φ≠0时，时间离散格式具有二阶精度；

公式(14)是一个非线性方程组，可以采用Newton迭代法求解，然而虽然Newton迭代法有较快的局部收敛性，但它的全局收敛性比较差，因此这里采用一种具有全局收敛性的准Newton迭代法：即修正的Levenberg-Marquardt 法，因而可进一步写为公式(15)，通过迭代求解速度矢量

其中：

F_k＝F(v_k)

式中：

J_k——第k次Newton迭代的Jacobian矩阵；

——表示J_k的转置；

——准Newton迭代步长；

I——单位矩阵；

λ_k—小参数，可由λ_k＝μ_k||F_k||^δ计算，其中δ∈[1,2]，μ_k是引入的为保持矩阵

的正定性的参数；

v_k，v_k+1——待求解的速度矢量在k，k+1次迭代中的值；

F_k——非线性函数F在k步的取值；

Δt——时间步长；

φ——比例系数，当φ≠0时，时间离散格式具有二阶精度；

f——v的非线性函数，见公式(12)。

进一步，如上所述的用于严重事故下安全壳内流动分析的计算方法，步骤6所述流道空泡数定义为流道中气体占据的体积份额，它只能是一个正数；一个流道连接到两个控制体，每一个连接处都可定义一个空泡数；逻辑上游的空泡数由公式(16)定义

式中：

——流道逻辑上游的空泡数；

——流道上游控制体流道开口顶部标高；

——流道上游控制体流道开口底部标高；

——流道上游控制体水面标高；

流道逻辑下游的空泡数由公式(17)定义

式中：

——流道逻辑下游的空泡数；

——流道下游控制体流道开口顶部标高；

——流道下游控制体流道开口底部标高；

——流道下游控制体水面标高；

如果流道中没有双向流动，则流道的空泡数α_j完全由流动的上游所确定，即是，如果流动速度为正值，则

否则，

当流道中气体与液体共存且流动方向相反时，流道空泡数可取为上、下游空泡数的某种加权平均，计算公式见公式(18)

式中：

α_j——流道空泡数；

——流道逻辑上游的空泡数；

——流道逻辑下游的空泡数；

v_g——流道气相速度；

v_w——流道液相速度；

ρ_g——流道气相密度；

ρ_w——流道液相密度；

进一步，如上所述的用于严重事故下安全壳内流动分析的计算方法，步骤7中，如果流道空泡数接近0，则令其为0，且将气相速度置0；如果流道空泡数接近1，则令其为1；如果液相份额接近0，则令其为0，且将液相速度置0；如果液相份额接近1，则令其为1。

进一步，如上所述的用于严重事故下安全壳内流动分析的计算方法，步骤8所述的上、下游包含两种概念，一种是规定的流动的正方向，它代表了流道连接的始末端，可以称为逻辑上、下游；另一种是真实流动的流动方向，当流动速度为正值时代表流动方向与规定的流动方向一致，为负值时代表流动方向与规定的流动方向相反，可以称为流动上、下游；

所述流动上、下游是质量、动量、能量交换的依据，对气体而言，当流道中的速度为正值时，流道的起始端(逻辑上游)为流动上游；当速度为负值时，终止端(逻辑下游)为流动上游；当速度为0时，取气体密度或者压强更高的一端为流动上游；对水而言，当速度不为0时，流动上游的判定与气体的情况是类似的，当速度为0时，上游取密度较大的控制体为上游。

进一步，如上所述的用于严重事故下安全壳内流动分析的计算方法，步骤9中，每个控制体中，每一相在一个时间步中产生的增量可由公式(35)(36) 计算

式中：

——控制体内相

的质量变化量；

——控制体内质量源引起的相

质量变化率；

Δt——时间步长；

——控制体内相

的内能变化量；

——控制体能量源引起的相

的内能变化率；

σ_ij——控制体与流道的关系，可以取为1、-1、0，分别表示流进、流出控制体或者与控制体不相连；

A_j——流道面积；

——流道相

的体积分数；

——流道上游控制体相

的密度；

——流道中相

速度；

——流道上游控制体相

的比焓。

本发明的有益效果如下：本发明是针对严重事故下安全壳内流动分析问题，提出的一种基于集总参数法的动量方程计算分析方法。在使用集总参数法进行安全壳流动分析时，通常将安全壳空间按一定规则划分为若干个控制体，各个控制体按照实际连接关系使用流道进行连接，通过对流道内动量方程的求解来得到实际的流动情况。本发明通过对基本方程进行本构分析，推导出一种考虑了主要状态变量和热工参数对压力影响的压力线性化方法，并保留了动量方程中的非定常项，从而提出一种新型的用于严重事故下安全壳内流动分析的动量方程计算方法，在提高计算效率的同时又保证了计算精度。本发明提出的动量方程计算方法可准确并且快速的求解流道中的动量方程，从而得到准确的安全壳内由流动引起的质量能量变化情况，为安全壳内其他严重事故现象的分析打下基础。

附图说明

图1为具体实施例中用于严重事故下安全壳内流动分析的计算方法的流程图。

图2为具体实施例中安全壳内控制体流道连接示意图。

具体实施方式

以下结合附图对本发明的具体实施方式作进一步的说明。在具体实施方式部分中，各公式的编号可能与发明内容部分相同公式的编号不同，相应编号只代表在该部分的公式顺序情况。

如图1所示，本发明用于严重事故下安全壳内流动分析的计算方法，步骤如下：

步骤1：获取安全壳流动分析计算所需参数并初始化，即获取待分析安全壳流动分析算例中控制体、流道的相关参数；图2展示了安全壳内控制体流道连接示意图，即流道两端分别连接一个控制体，一端定义为流道上游控制体，一端定义为流道下游控制体，控制体内的气相和液相通过流道进行流动；

根据图2所示，计算所需的初始参数如下：

流道上、下游控制体几何参数：底面积S，水面高程Z_w等；

流道上、下游控制体上一时刻状态参数：温度T，压力P，气相各组分(不可凝气体及蒸汽)及液相的质量m、密度ρ、比热c_v、比焓h，比内能u等；

流道参数：流道长度L，阻力系数K，面积A，上游、下游开口标高Z、连接处高程Z_J以及控制体与流道关系σ等；

步骤2：建立流动状态方程，流动现象经过模型化后最终可以归结为一组集总参数形式的、非线性的、耦合的一阶常微分方程，构成了一个动力学系统，其相应的常微分方程如公式(1)；

式中：

j——第j个流道；

——某一相流体，可为气相或液相；

d——流道的上游；

r——流道的下游；

——流道j上游控制体相

密度；

——流道j相

速度；

t——时间；

L_j——流道j中非定常流动的长度，一般取为流道长度；

——流道上下游控制体中的热力学压强；

——流道j两端液相和气相各自的重力头差；

ΔP_j——流体机械产生的泵头；

K_j——无量纲的流道阻力系数。

步骤3：计算重力压头，连接控制体i和k的流道j中作用于相

的压力差为

其中P_i和P_k是相应控制体积中的热力学压力，

包含控制体内和沿着流道的所有重力头的项；

在不考虑流动时，重力压头

包含三部分，第一部分是控制体i中水池面高程z_w,i与连接处各相高程

之间的水头差

计算公式如公式 (2)所示：

式中：

——控制体i中水池面高程z_p,i与连接处各相高程

之间的水头差；

ρ_w,i——控制体i中液相密度；

ρ_g,i——控制体i中气相密度；

g——重力加速度；

Z_w,i——控制体i液面高程；

——控制体i连接处各相高程；

第二部分为控制体k中相应的压力差，它与第一部分完全类似，如公式 (3)所示：

式中：

——控制体k中水池面高程z_w,k与连接处各相高程

之间的水头差；

ρ_w,k——控制体k中液相密度；

ρ_g,k——控制体k中气相密度；

g——重力加速度；

Z_w,k——控制体k液面高程；

——控制体k连接处各相高程；

第三部分是各相

沿着流道的重力头，其计算如公式(4)所示：

式中：

——各相

沿着流道的重力压头；

——流道中相

的密度，取为流道上游与下游控制体内该相密度的最大值；

g——重力加速度；

——控制体i连接处各相高程；

——控制体k连接处各相高程；

综上，纯重力压头定义为这三部分之和，见公式(5)：

式中：

——控制体内和沿着流道的所有重力压头的项；

——控制体i中水池面高程z_w,i与连接处各相高程

之间的水头差；

——控制体k中水池面高程z_w,k与连接处各相高程

之间的水头差；

——各相

沿着流道的重力压头；

将重力压头

进行线性化，考虑一个时间步长内流动引起的流道两端的重力压头的变化；对于气体而言，流动引起的重力压头的变化相较于压力变化是可以忽略的；水的重力压头的变化则相对不是很小，忽略流道中因为水的流动引起的重力压头变化，则流道两端水的重力压头变化为流道两端控制体水的重力压头变化之和，计算公式见公式(6)

式中：

Δ(ρgΔz)_j,w——一个时间步长内因水的流动引起的流道两端水的重力压头的变化；

Δm_i、Δm_k——控制体i、k内水的质量变化，以控制体i为例，其等于

其中Δt为时间步长，σ_ij表示控制体与流道j的关系，A_j为流道面积，α_j,w为流道中液相空泡数，

流道上游控制体液相密度，v_j,w流道液相速度；

g——重力加速度；

S_i、S_k——控制体i、k的底面积；

i——控制体i；

k——控制体k；

Z_w,i、Z_w,k——控制体i、k中水池面高程；

Z_J,w,i、Z_J,w,k——控制体i、k连接处液相高程；

进一步，假定控制体内水的密度在一个时间步长内不变，以控制体i为例，一个时间步长内控制体中的水的重力压头变化可由公式(7)计算

式中：

g——重力加速度；

Δm——控制体内水的质量变化；

S——控制体底面积；

Δt——时间步长；

σ_ij——控制体与流道j的关系，可以取为1、-1、0，分别表示流进、流出控制体或者与控制体不相连；

A_j——流道面积；

α_j,w——流道中液相空泡数；

——流道上游控制体液相密度；

v_j,w——流道液相速度；

步骤4：压力及其线性化；动力学系统可以采用显式的时间推进方法求解，但是为了时间稳定性，需要很小的时间步长；采用隐式的时间推进方法，则使用较大的时间步长也能维持数值稳定性，而进行隐式求解时重要的是对压力线性化，即是将作为新时刻动力学的非线性函数的压力在旧时刻的压力的基础上展开为新时刻动力学变量的线性函数；

压力包括两部分，一部分是蒸汽分压P_s(T,ρ)，它是蒸汽温度T和蒸汽密度ρ_s的函数(安全壳中蒸汽没有达到饱和，它的温度T_s与气体温度T_g相同，但与水温T_w不同；稳压器中蒸汽达到饱和，则蒸汽温度T_s与水温T_w相同，饱和蒸汽压力P(T_s)只是温度T_s的函数)，另一部分是气体分压，它服从理想气体状态方程，从而压力可由公式(8)表述

式中：

P——控制体总压力；

P_s(T,ρ)——蒸汽分压，其为温度T和密度ρ的函数；

R——理想气体常数；

T_g——气相温度；

V_g——气相体积；

m_k——气相组分k的质量；

M_k——气相组分k的摩尔质量；

m_g——气相各种不可凝气体的质量总和；

x_k——气相组分k的质量份额；

P_g——气相各种不可凝气体压力；

进一步，对公式(8)求导，可得公式(9)

式中：

P——控制体总压力；

t——时间；

R——理想气体常数；

P_s——蒸汽分压；

ρ_s——蒸汽密度；

m_s——蒸汽的质量；

T_g——气相温度；

V_g——气相体积；

m_k——气相不可凝气体组分k的质量；

M_k——气相不可凝气体组分k的摩尔质量；

其中气体体积随时间的变化率

可以与水的质量、密度的变化联系起来，见公式(10)

——气体体积随时间的变化率；

t——时间；

V_w——水的体积；

m_w——水的质量；

ρ_w——水的密度；

β——水的膨胀系数；

ρ₀——水密度的参考状态，1000kg/m³；

T_w——水的温度；

a_w——水中的声速；

P——压力；

式中dρ_w/dt主要是由水温度的变化引起的，因而可得到公式(11)：

式中：

——气体体积随时间的变化率；

t——时间；

m_w——水的质量；

ρ_w——水的密度；

β——水的膨胀系数；

T_w——水的温度；

a_w——水中的声速；

P——总压力；

记

式中：

a_w,a_s——水和蒸汽中的声速；

T_g——气相温度；

R——理想气体常数；

M_w——水的摩尔质量；

ρ_w——水的密度；

P——总压力；

P_s——蒸汽的分压；

ρ_s——蒸汽的密度；

z_s——蒸汽压缩因子；

将公式(12)带入公式(11)(10)(9)，并将公式(10)带入公式(9)，可得

式中：

t——时间；

m_w——水的质量；

ρ_w——水的密度；

β——水的膨胀系数；

T_w——水的温度；

a_w——水中的声速；

P——总压力；

T_g——气相温度；

V_g——气相体积；

R——理想气体常数；

m_k——气相不可凝气体组分k质量；

M_k——气相不可凝气体组分k摩尔质量；

P——总压力；

m_s——蒸汽的质量；

ρ₀——水密度的参考状态，1000kg/m³；

P_s——蒸汽的分压；

ρ_s——蒸汽的密度；

a_s——蒸汽中的声速；

进一步，将公式(13)连续性方程代入公式(14)可化成公式(15)：

式中：

t——时间；

m_w——水的质量；

ρ_w——水的密度；

β——水的膨胀系数；

T_w——水的温度；

a_w——水中的声速；

P——总压力；

T_g——气相温度；

V_g——气相体积；

R——理想气体常数；

m_k——气相不可凝气体组分k质量；

M_k——气相不可凝气体组分k摩尔质量；

P——总压力；

m_s——蒸汽的质量；

ρ₀——水密度的参考状态，1000kg/m³；

P_s——蒸汽的分压；

ρ_s——蒸汽的密度；

a_s——蒸汽中的声速；

σ_ij——控制体与流道的关系，可以取为1、-1、0，分别表示流进、流出控制体或者与控制体不相连；

A_j——流道面积；

α_j,g——流道气相的体积分数；

α_j,w——流道液相的体积分数；

——流道上游控制体水蒸汽密度；

——流道上游控制体不可凝气体组分k的密度；

——流道上游控制体水的密度；

v_j,g——流道中气相速度；

v_j,w——流道中液相的速度；

进一步，记

式中：

m_g——气相不可凝气体总质量；

T_g——气相温度；

V_g——气相体积；

R——理想气体常数；

M_w——水的摩尔质量；

m_k——气相不可凝气体组分k质量；

M_k——气相不可凝气体组分k摩尔质量； P——总压力；

m_s——蒸汽的质量；

z_s——蒸汽压缩因子；

P_s——蒸汽的分压；

ρ_s——蒸汽的密度；

ρ_g——气相不可凝气体密度；

a_s——蒸汽中的声速；

则公式(15)可化为公式(17)

t——时间；

a——公式(16)表示的系数a；

b——公式(16)表示的系数b；

m_w——水的质量；

ρ_w——水的密度；

β——水的膨胀系数；

T_w——水的温度；

a_w——水中的声速；

P——总压力；

T_g——气相温度；

V_g——气相体积；

R——理想气体常数；

m_k——气相不可凝气体组分k质量；

M_k——气相不可凝气体组分k摩尔质量；

P——总压力；

m_s——蒸汽的质量；

ρ₀——水密度的参考状态，1000kg/m³；

P_s——蒸汽的分压；

ρ_s——蒸汽的密度；

a_s——蒸汽中的声速；

σ_ij——控制体与流道的关系，可以取为1、-1、0，分别表示流进、流出控制体或者与控制体不相连；

A_j——流道面积；

α_j,g——流道气相的体积分数；

α_j,w——流道液相的体积分数；

——流道上游控制体水的密度；

v_j,g——流道中气相速度；

v_j,w——流道中水的速度；

由于公式(17)中

一般是一个小量，因而公式(17)可进一步化为公式(18)

式中：

t——时间；

a——公式(16)表示的系数a；

b——公式(16)表示的系数b；

m_w——水的质量；

ρ_w——水的密度；

β——水的膨胀系数；

T_w——水的温度；

a_w——水中的声速；

P——总压力；

T_g——气相温度；

V_g——气相体积；

R——理想气体常数；

m_k——气相不可凝气体组分k质量；

M_k——气相不可凝气体组分k摩尔质量；

P——总压力；

m_s——蒸汽的质量；

ρ₀——水密度的参考状态，1000kg/m³；

P_s——蒸汽的分压；

ρ_s——蒸汽的密度；

a_s——蒸汽中的声速；

σ_ij——控制体与流道的关系，可以取为1、-1、0，分别表示流进、流出控制体或者与控制体不相连；

A_j——流道面积；

α_j,g——流道气相的体积分数；

α_j,w——流道液相的体积分数；

——流道上游控制体水的密度；

v_j,g——流道中气相速度；

v_j,w——流道中水的速度；

进一步，根据公式(19)对时间求导，可得公式(20)，并对公式(20) 变形，最终得到公式(21)

U_i＝m_iu_i 公式(19)

式中：

U_i——物质i的总内能；

m_i——物质i的质量；

u_i——物质i的比内能；

t——时间；

c_vi——物质i的比热；

T_i——物质i的温度；

而根据控制体质量常微分方程和能量常微分方程，可得质量和能量变化量计算公式(22)和公式(23)

式中：

i——第i个控制体；

j——第j个流道；

——某一相流体，可为气相或液相；

d——流道的上游；

r——流道的下游；

——控制体i中表示控制体中相

流体的总质量；

σ_ij——表示控制体i与流道j的关系，可以取为1，-1，0，分别表示流进，流出控制体，或者与控制体不相连；

A_j——流道j面积；

——流道j相

的体积分数；

——流道j上游控制体相

密度；

——流道j相

速度；

——控制体i中相

的质量源产生的速率；

——控制体i中相

的总能量；

t——时间；

——流道j上游控制体中相

的比焓；

——控制体i中相

能量源产生的速率；

分别对气相和液相使用公式(21)，并带入公式(22)(23)，可得公式(24)

t——时间；

c_vg——气相气体比热；

T_g——气相温度；

U_g——气相内能；

u_g——气相比内能；

m_w——水的质量；

c_vw——水的比热；

T_w——水的温度；

U_w——水的内能；

u_w——水的比内能；

σ_ij——表示控制体i与流道j的关系，可以取为1，-1，0，分别表示流进，流出控制体，或者与控制体不相连；

A_j——流道j面积；

α_j,g——流道气相的体积分数；

α_j,w——流道液相的体积分数；

——流道上游控制体气相的密度；

——流道上游控制体水的密度；

v_j,g——流道中气相速度；

v_j,w——流道中水的速度；

——流道上游控制体气相的比焓；

——流道上游控制体水的比焓；

δm_g——总气体质量源变化率；

δm_w——液相水质量源变化率；

δU_g——气相能量源项变化率；

δU_w——液相水能量源项变化率；

最终，将公式(24)带入公式(18)，可得

计算公式(25)

式中：

t——时间；

P——总压力；

a——公式(16)表示的系数a；

b——公式(16)表示的系数b；

R——理想气体常数；

m_k——气相不可凝气体组分k质量；

M_k——气相不可凝气体组分k摩尔质量；

ρ_k——气相不可凝气体组分k密度；

m_s——蒸汽的质量；

ρ₀——水密度的参考状态，1000kg/m³；

a_s——蒸汽中的声速；

c_vg——气相气体比热；

T_g——气相温度；

V_g——气相体积；

u_g——气相比内能；

m_w——水的质量；

ρ_w——水的密度；

c_vw——水的比热；

T_w——水的温度；

u_w——水的比内能；

β——水的膨胀系数；

σ_ij——控制体与流道的关系，可以取为1、-1、0，分别表示流进、流出控制体或者与控制体不相连；

A_j——流道面积；

α_j,g——流道气相的体积分数；

α_j,w——流道液相的体积分数；

——流道上游控制体气相的密度；

——流道上游控制体液相的密度；

v_j,g——流道中气相速度；

v_j,w——流道中液相速度；

——流道上游控制体气相的比焓；

——流道上游控制体液相的比焓；

δm_g——总气体质量源变化率；

δm_w——液相水质量源变化率；

δm_s——水蒸汽质量源变化率；

δm_k——气相不可凝气体k的质量源变化率；

δU_g——气相能量源项变化率；

δU_w——液相水能量源项变化率；

因而，新时刻压力可由公式(26)计算，其中压力对时间的微分

的可由公式(25)给出

式中：

P^new——新时刻压力；

P^old——旧时刻压力；

——压力对时间的微分；

Δt——时间步长；

步骤5：利用准Newton迭代法求解动量方程，得到新时刻速度；根据步骤3、4，使用压力线性化表示新时刻压力，并将重力压头代入，即将公式(7) 和公式(25)代入步骤1中的公式(1)，分别可得到气相和液相的方程，见公式(27)和公式(28)

式中：

t——时间；

Δt——时间步长；

P——总压力；

a——公式(16)表示的系数a；

b——公式(16)表示的系数b；

(ρgΔz)_j,g——流道中气相产生的重力头；

(ρgΔz)_j,w——流道中液相产生的重力头；

Δ(ρgΔz)_j,w——由于流动而在一个时间步长内产生的液相重力头变化；

——流道上游控制体旧时刻压力；

——流道下游控制体旧时刻压力；

S——控制体底面积；

ΔP_j——流体机械产生的泵头；

R——理想气体常数；

m_k——气相不可凝气体组分k质量；

M_k——气相不可凝气体组分k摩尔质量；

ρ_k——气相不可凝气体组分k密度；

m_s——蒸汽的质量；

ρ₀——水密度的参考状态，1000kg/m³；

a_s——蒸汽中的声速；

c_vg——气相气体比热；

T_g——气相温度；

V_g——气相体积；

u_g——气相比内能；

m_w——水的质量；

ρ_w——水的密度；

c_vw——水的比热；

T_w——水的温度；

u_w——水的比内能；

β——水的膨胀系数；

σ_ij——控制体与流道的关系，可以取为1、-1、0，分别表示流进、流出控制体或者与控制体不相连；

A_j——流道面积；

L_j——流道中非定常流动的长度，一般取为流道长度；

K_j——流道阻力系数；

α_j,g——流道气相的体积分数；

α_j,w——流道液相的体积分数；

——流道上游控制体气相的密度；

——流道上游控制体水的密度；

v_j,g——流道中气相速度；

v_j,w——流道中液相速度；

——流道上游控制体气相的比焓；

——流道上游控制体液相的比焓；

方程(26)和(27)的主要区别是前者不包含一个时间步长内静水头的变化，相同点是两者都包含了一个时间步长内由流动引起的热力学压力的变化，即由热力学压力线性化产生的速度的线性项是相同的；

方程(26)和(27)是非线性常微分方程，它们构成了一个常微分方程系统，将常微分方程系统简记为公式(28)

式中f(v)是v的非线性函数，v为待求解速度矢量；

进一步，采用半显半隐式时间离散动力系统，公式(28)可离散为公式 (29)

式中：

Δt——时间步长；

vⁿ⁺¹——新时刻速度矢量；

vⁿ——当前时刻速度矢量；

v^n-1——旧时刻速度矢量；

φ——比例系数，当φ≠0时，时间离散格式具有二阶精度；

公式(29)是一个非定常、非线性方程组，可采用迭代法求解，因而可写为公式(30)

F(vⁿ⁺¹)＝(1+φ)(vⁿ⁺¹-vⁿ)-φ(vⁿ-v^n-1)-Δtf(vⁿ⁺¹)＝0 公式(30)

F(vⁿ⁺¹)——速度vⁿ⁺¹的非线性函数；

Δt——时间步长；

vⁿ⁺¹——新时刻速度矢量；

vⁿ——当前时刻速度矢量；

v^n-1——旧时刻速度矢量；

φ——比例系数，当φ≠0时，时间离散格式具有二阶精度；

公式(30)是一个非线性方程组，可以采用Newton迭代法求解，然而虽然Newton迭代法有较快的局部收敛性，但它的全局收敛性比较差，因此这里采用一种具有全局收敛性的准Newton迭代法：即修正的Levenberg-Marquardt 法，因而可进一步写为公式(31)，通过迭代求解速度矢量

其中：

F_k＝F(v_k)

式中：

J_k——第k次Newton迭代的Jacobian矩阵；

——表示J_k的转置；

——准Newton迭代步长；

I——单位矩阵；

λ_k—小参数，可由λ_k＝μ_k||F_k||^δ计算，其中δ∈[1,2]，μ_k是引入的为保持矩阵

的正定性的参数；

v_k，v_k+1——待求解的速度矢量在k，k+1次迭代中的值；

F_k——非线性函数F在k步的取值；

Δt——时间步长；

φ——比例系数，当φ≠0时，时间离散格式具有二阶精度；

f——v的非线性函数，见公式(28)；

步骤6：根据预测的新时刻速度求解流道空泡数；流道空泡数定义为流道中气体占据的体积份额，它只能是一个正数；一个流道连接到两个控制体，每一个连接处都可定义一个空泡数；逻辑上游的空泡数由公式(32)定义

式中：

——流道逻辑上游的空泡数；

——流道上游控制体流道开口顶部标高；

——流道上游控制体流道开口底部标高；

——流道上游控制体水面标高；

流道逻辑下游的空泡数由公式(33)定义

式中：

——流道逻辑下游的空泡数；

——流道下游控制体流道开口顶部标高；

——流道下游控制体流道开口底部标高；

——流道下游控制体水面标高；

如果流道中没有双向流动，则流道的空泡数α_j完全由流动的上游所确定，即是，如果流动速度为正值，则

否则，

当流道中气体与液体共存且流动方向相反时，流道空泡数可取为上、下游空泡数的某种加权平均，计算公式见公式(34)

式中：

α_j——流道空泡数；

——流道逻辑上游的空泡数；

——流道逻辑下游的空泡数；

v_g——流道气相速度；

v_w——流道液相速度；

ρ_g——流道气相密度；

ρ_w——流道液相密度；

步骤7：根据空泡数修正预测速度以得到新时刻速度，如果流道空泡数接近0，则令其为0，且将气相速度置0；如果流道空泡数接近1，则令其为1；如果液相份额接近0，则令其为0，且将液相速度置0；如果液相份额接近1，则令其为1；

步骤8：根据新时刻速度确定流道上、下游，这里上、下游有两种概念，一种是规定的流动的正方向，它代表了流道连接的始末端，可以称为逻辑上、下游；另一种是真实流动的流动方向，当流动速度为正值时代表流动方向与规定的流动方向一致，为负值时代表流动方向与规定的流动方向相反，可以称为流动上、下游；

流动上、下游是质量、动量、能量交换的依据，物质总是从上游流向下游，判定流动上、下游是重要的；对气体而言，当流道中的速度为正值时，流道的起始端(逻辑上游)为流动上游；当速度为负值时，终止端(逻辑下游)为流动上游；当速度为0时，取气体密度或者压强更高的一端为流动上游；对水而言，当速度不为0时，流动上游的判定与气体的情况是类似的，当速度为0时，上游取密度较大的控制体为上游；

因而，根据上一步骤求解出新时刻的速度确定流动上、下游；

步骤9：根据新时刻速度求解质量、能量增量；求得速度后，质量、能量方程的求解是简单而直接的，每个控制体中，每一相在一个时间步中产生的增量可由公式(35)(36)计算

式中：

——控制体内相

的质量变化量；

——控制体内质量源引起的相

质量变化率；

Δt——时间步长；

——控制体内相

的内能变化量；

——控制体能量源引起的相

的内能变化率；

σ_ij——控制体与流道的关系，可以取为1、-1、0，分别表示流进、流出控制体或者与控制体不相连；

A_j——流道面积；

——流道相

的体积分数；

——流道上游控制体相

的密度；

——流道中相

速度；

——流道上游控制体相

的比焓；

本方法根据输入的安全壳模型参数，根据步骤1-9进行求解计算，最终得到当前时间步长流道内的速度，以及对应的质量增量和能量增量，从而得到安全壳内的流动速度分布以及流动引起的质量和能量变化。这些求解的流动参数是进行安全壳安全分析的重要热工水力参数，是进行核电厂安全分析、严重事故计算基础条件。因此，作为求解热工水力参数的基础方法，本方法可应用于集总参数法的安全壳热工水力软件、严重事故计算分析软件中，用于安全壳内流动参数的求解。

对于本领域技术人员而言，显然本发明方法不限于上述示范性实施例的细节，而且在不背离本发明的精神或基本特征的情况下，能够以其他的具体形式实现本发明方法。因此，无论从哪一点来看，均应将实施例看作是示范性的，而且是非限制性的，本发明方法的范围由所附权利要求而不是上述说明限定，因此旨在将落在权利要求的等同要件的含义和范围内的所有变化囊括在本发明方法内。不应将权利要求中的任何附图标记视为限制所涉及的权利要求。

此外，应当理解，虽然本说明书按照实施方式加以描述，但并非每个实施方式仅包含一个独立的技术方案，说明书的这种叙述方式仅仅是为清楚起见，本领域技术人员应当将说明书作为一个整体，各实施例中的技术方案也可以经适当组合，形成本领域技术人员可以理解的其他实施方式。

Claims (9)

1.一种用于严重事故下安全壳内流动分析的计算方法，其特征在于，包括如下步骤：

步骤1、获取安全壳流动分析计算所需参数并初始化；

步骤2、建立流动状态方程，流动现象经过模型化后最终可以归结为一组集总参数形式的、非线性的、耦合的一阶常微分方程，构成了一个动力学系统；

步骤3、计算重力压头，连接控制体i和k的流道j中作用于相

的压力差为

其中P_i和P_k是相应控制体积中的热力学压力，

包含控制体内和沿着流道的所有重力头的项；

在不考虑流动时，重力压头

包含三部分，第一部分是控制体i中水池面高程z_w,i与连接处各相高程

之间的水头差，第二部分为控制体k中水池面高程z_w,k与连接处各相高程

之间的水头差，第三部分是各相

沿着流道的重力压头，纯重力压头定义为以上三部分之和；

步骤4、进行压力线性化，将动量方程中压力项线性化，将作为新时刻动力学的非线性函数的压力在旧时刻的压力的基础上展开为新时刻动力学变量的线性函数；

步骤5、根据步骤3、4，使用压力线性化表示新时刻压力，并将重力压头代入，利用准Newton迭代法求解动量方程；

步骤6、根据预测的新时刻速度求解流道空泡数；

步骤7、根据空泡数修正预测速度以得到新时刻速度；

步骤8、根据新时刻速度确定流道上、下游；

步骤9、根据新时刻速度求解质量、能量增量。

2.如权利要求1所述的用于严重事故下安全壳内流动分析的计算方法，其特征在于，步骤2中所述动力学系统相应的常微分方程如下：

式中：

j——第j个流道；

——某一相流体，可为气相或液相；

d——流道的上游；

r——流道的下游；

——流道j上游控制体相

密度；

——流道j相

速度；

t——时间；

L_j——流道j中非定常流动的长度，一般取为流道长度；

——流道上下游控制体中的热力学压强；

——流道j两端液相和气相各自的重力头差；

ΔP_j——流体机械产生的泵头；

K_j——无量纲的流道阻力系数。

3.如权利要求1或2所述的用于严重事故下安全壳内流动分析的计算方法，其特征在于，步骤3中所述重力压头

包含的三部分中，第一部分是控制体i中水池面高程z_w,i与连接处各相高程

之间的水头差

计算公式如下：

式中：

——控制体i中水池面高程z_p,i与连接处各相高程

之间的水头差；

ρ_w,i——控制体i中液相密度；

ρ_g,i——控制体i中气相密度；

g——重力加速度；

Z_w,i——控制体i液面高程；

——控制体i连接处各相高程；

第二部分为控制体k中水池面高程z_wk与连接处各相高程

之间的水头差，计算公式如下：

式中：

——控制体k中水池面高程z_w,k与连接处各相高程

之间的水头差；

ρ_w，k——控制体k中液相密度；

ρ_g,k——控制体k中气相密度；

g——重力加速度；

Z_w,k——控制体k液面高程；

——控制体k连接处各相高程；

第三部分是各相

沿着流道的重力压头，其计算如下：

式中：

——各相

沿着流道的重力压头；

——流道中相

的密度，取为流道上游与下游控制体内该相密度的最大值；

g——重力加速度；

——控制体i连接处各相高程；

——控制体k连接处各相高程；

综上，纯重力压头定义为这三部分之和，如下所示：

式中：

——控制体内和沿着流道的所有重力压头的项；

——控制体i中水池面高程z_w,i与连接处各相高程

之间的水头差；

——控制体k中水池面高程z_w,k与连接处各相高程

之间的水头差；

——各相

沿着流道的重力压头；

将重力压头

进行线性化，考虑一个时间步长内流动引起的流道两端的重力压头的变化；对于气体而言，流动引起的重力压头的变化相较于压力变化是可以忽略的；水的重力压头的变化则相对不是很小，忽略流道中因为水的流动引起的重力压头变化，则流道两端水的重力压头变化为流道两端控制体水的重力压头变化之和，计算公式如下：

式中：

Δ(ρgΔz)_j,w——一个时间步长内因水的流动引起的流道两端水的重力压头的变化；

Δm_i、Δm_k——控制体i、k内水的质量变化，以控制体i为例，其等于

其中Δt为时间步长，σ_ij表示控制体与流道j的关系，A_j为流道面积，α_j,w为流道中液相空泡数，

流道上游控制体液相密度，v_j,w流道液相速度；

g——重力加速度；

S_i、S_k——控制体i、k的底面积；

i——控制体i；

k——控制体k；

Z_w,i、Z_w,k——控制体i、k中水池面高程；

Z_J,w,i、Z_J,w,k——控制体i、k连接处液相高程；

进一步，假定控制体内水的密度在一个时间步长内不变，以控制体i为例，一个时间步长内控制体中水的重力压头变化可由下式计算：

式中：

g——重力加速度；

Δm——控制体内水的质量变化；

S——控制体底面积；

Δt——时间步长；

σ_ij——控制体与流道j的关系，可以取为1、-1、0，分别表示流进、流出控制体或者与控制体不相连；

A_j——流道面积；

α_j,w——流道中液相空泡数；

——流道上游控制体液相密度；

v_j,w——流道液相速度。

4.如权利要求1所述的用于严重事故下安全壳内流动分析的计算方法，其特征在于，步骤4中进行压力线性化，其中新时刻压力可由下式计算

式中：

P^new——新时刻压力；

P^old——旧时刻压力；

——压力对时间的微分；

Δt——时间步长；

公式(8)中的压力对时间的微分可使用公式(9)计算

式中：

t——时间；

P——总压力；

a——等于

其中R为气体常数，T_g为气相温度，V_g为气相体积，z_s为蒸汽压缩因子，M_w为水的摩尔质量，ρ_s为蒸汽的密度，ρ_g为气相不可凝气体密度，x_k为气相不可凝气体组分k的质量份额，M_k为气相不可凝气体组分k摩尔质量；

b——等于

其中R为气体常数，z_s为蒸汽压缩因子，M_w为水的摩尔质量，ρ_s为蒸汽的密度，ρ_g为气相不可凝气体密度，x_k为气相不可凝气体组分k的质量份额，M_k为气相不可凝气体组分k摩尔质量；

R——理想气体常数；

m_k——气相不可凝气体组分k质量；

M_k——气相不可凝气体组分k摩尔质量；

ρ_k——气相不可凝气体组分k密度；

m_s——蒸汽的质量；

ρ₀——水密度的参考状态，1000kg/m³；

a_s——蒸汽中的声速；

c_vg——气相气体比热；

T_g——气相温度；

V_g——气相体积；

u_g——气相比内能；

m_w——水的质量；

ρ_w——水的密度；

c_vw——水的比热；

T_w——水的温度；

u_w——水的比内能；

β——水的膨胀系数；

σ_ij——控制体与流道的关系，可以取为1、-1、0，分别表示流进、流出控制体或者与控制体不相连；

A_j——流道面积；

α_j,g——流道气相的体积分数；

α_j,w——流道液相的体积分数；

——流道上游控制体气相的密度；

——流道上游控制体液相的密度；

v_j,g——流道中气相速度；

v_j,w——流道中液相速度；

——流道上游控制体气相的比焓；

——流道上游控制体液相的比焓；

δm_g——总气体质量源变化率；

δm_w——液相水质量源变化率；

δm_s——水蒸汽质量源变化率；

δm_k——气相不可凝气体k的质量源变化率；

δU_g——气相能量源项变化率；

δU_w——液相水能量源项变化率。

5.如权利要求1所述的用于严重事故下安全壳内流动分析的计算方法，其特征在于，步骤5中将公式(7)和公式(8)代入步骤1中的公式(1)，分别可得到气相和液相的动量方程，见公式(10)和公式(11)

式中：

t——时间；

Δt——时间步长；

P——总压力；

a——等于

其中R为气体常数，T_g为气相温度，V_g为气相体积，z_s为蒸汽压缩因子，M_w为水的摩尔质量，ρ_s为蒸汽的密度，ρ_g为气相不可凝气体密度，x_k为气相不可凝气体组分k的质量份额，M_k为气相不可凝气体组分k摩尔质量；

b——等于

其中R为气体常数，z_s为蒸汽压缩因子，M_w为水的摩尔质量，ρ_s为蒸汽的密度，ρ_g为气相不可凝气体密度，x_k为气相不可凝气体组分k的质量份额，M_k为气相不可凝气体组分k摩尔质量；

(ρgΔz)_j,g——流道中气相产生的重力头；

(ρgΔz)_j,w——流道中液相产生的重力头；

Δ(ρgΔz)_j,w——由于流动而在一个时间步长内产生的液相重力头变化；

——流道上游控制体旧时刻压力；

——流道下游控制体旧时刻压力；

S——控制体底面积；

ΔP_j——流体机械产生的泵头；

R——理想气体常数；

m_k——气相不可凝气体组分k质量；

M_k——气相不可凝气体组分k摩尔质量；

ρ_k——气相不可凝气体组分k密度；

m_s——蒸汽的质量；

ρ₀——水密度的参考状态，1000kg/m³；

a_s——蒸汽中的声速；

c_vg——气相气体比热；

T_g——气相温度；

V_g——气相体积；

u_g——气相比内能；

m_w——水的质量；

ρ_w——水的密度；

c_vw——水的比热；

T_w——水的温度；

u_w——水的比内能；

β——水的膨胀系数；

σ_ij——控制体与流道的关系，可以取为1、-1、0，分别表示流进、流出控制体或者与控制体不相连；

A_j——流道面积；

L_j——流道中非定常流动的长度，一般取为流道长度；

K_j——流道阻力系数；

α_j,g——流道气相的体积分数；

α_j,w——流道液相的体积分数；

——流道上游控制体气相的密度；

——流道上游控制体水的密度；

v_j,g——流道中气相速度；

v_j,w——流道中液相速度；

——流道上游控制体气相的比焓；

——流道上游控制体液相的比焓；

方程(10)和(11)是非线性常微分方程，它们构成了一个常微分方程系统，将常微分方程系统简记为公式(12)

式中f(v)是v的非线性函数，v为待求解速度矢量；

进一步，采用半显半隐式时间离散动力系统，公式(12)可离散为公式(13)

式中：

Δt——时间步长；

vⁿ⁺¹——新时刻速度矢量；

vⁿ——当前时刻速度矢量；

v^n-1——旧时刻速度矢量；

φ——比例系数，当φ≠0时，时间离散格式具有二阶精度；

公式(13)是一个非定常、非线性方程组，可采用迭代法求解，因而可写为公式(14)

F(vⁿ⁺¹)＝(1+φ)(vⁿ⁺¹-vⁿ)-φ(vⁿ-v^n-1)-Δtf(vⁿ⁺¹)＝0 公式(14)

F(vⁿ⁺¹)——速度vⁿ⁺¹的非线性函数；

Δt——时间步长；

vⁿ⁺¹——新时刻速度矢量；

vⁿ——当前时刻速度矢量；

v^n-1——旧时刻速度矢量；

φ——比例系数，当φ≠0时，时间离散格式具有二阶精度；

公式(14)是一个非线性方程组，采用一种具有全局收敛性的修正的Levenberg-Marquardt法求解，因而可进一步写为公式(15)，通过迭代求解速度矢量

其中：

F_k＝F(v_k)

式中：

J_k——第k次Newton迭代的Jacobian矩阵；

——表示J_k的转置；

——准Newton迭代步长；

I——单位矩阵；

λ_k—小参数，可由λ_k＝μ_k||F_k||^δ计算，其中δ∈[1,2]，μ_k是引入的为保持矩阵

的正定性的参数；

v_k，v_k+1——待求解的速度矢量在k，k+1次迭代中的值；

F_k——非线性函数F在k步的取值；

Δt——时间步长；

φ——比例系数，当φ≠0时，时间离散格式具有二阶精度；

f——v的非线性函数，见公式(12)。

6.如权利要求1所述的用于严重事故下安全壳内流动分析的计算方法，其特征在于，步骤6所述流道空泡数定义为流道中气体占据的体积份额，它只能是一个正数；一个流道连接到两个控制体，每一个连接处都可定义一个空泡数；逻辑上游的空泡数由公式(16)定义

式中：

——流道逻辑上游的空泡数；

——流道上游控制体流道开口顶部标高；

——流道上游控制体流道开口底部标高；

——流道上游控制体水面标高；

流道逻辑下游的空泡数由公式(17)定义

式中：

——流道逻辑下游的空泡数；

——流道下游控制体流道开口顶部标高；

——流道下游控制体流道开口底部标高；

——流道下游控制体水面标高；

如果流道中没有双向流动，则流道的空泡数α_j完全由流动的上游所确定，即是，如果流动速度为正值，则

否则，

当流道中气体与液体共存且流动方向相反时，流道空泡数可取为上、下游空泡数的某种加权平均，计算公式见公式(18)

式中：

α_j——流道空泡数；

——流道逻辑上游的空泡数；

——流道逻辑下游的空泡数；

v_g——流道气相速度；

v_w——流道液相速度；

ρ_g——流道气相密度；

ρ_w——流道液相密度。

7.如权利要求1所述的用于严重事故下安全壳内流动分析的计算方法，其特征在于，步骤7中，如果流道空泡数接近0，则令其为0，且将气相速度置0；如果流道空泡数接近1，则令其为1；如果液相份额接近0，则令其为0，且将液相速度置0；如果液相份额接近1，则令其为1。

8.如权利要求1所述的用于严重事故下安全壳内流动分析的计算方法，步骤8所述的上、下游包含两种概念，一种是规定的流动的正方向，它代表了流道连接的始末端，可以称为逻辑上、下游；另一种是真实流动的流动方向，当流动速度为正值时代表流动方向与规定的流动方向一致，为负值时代表流动方向与规定的流动方向相反，可以称为流动上、下游；

所述流动上、下游是质量、动量、能量交换的依据，对气体而言，当流道中的速度为正值时，流道的起始端(逻辑上游)为流动上游；当速度为负值时，终止端(逻辑下游)为流动上游；当速度为0时，取气体密度或者压强更高的一端为流动上游；对水而言，当速度不为0时，流动上游的判定与气体的情况是类似的，当速度为0时，上游取密度较大的控制体为上游。

9.如权利要求1所述的用于严重事故下安全壳内流动分析的计算方法，其特征在于，步骤9中，每个控制体中，每一相在一个时间步中产生的增量可由公式(35)(36)计算

式中：

——控制体内相

的质量变化量；

——控制体内质量源引起的相

质量变化率；

Δt——时间步长；

——控制体内相

的内能变化量；

——控制体能量源引起的相

的内能变化率；

σ_ij——控制体与流道的关系，可以取为1、-1、0，分别表示流进、流出控制体或者与控制体不相连；

A_j——流道面积；

——流道相

的体积分数；

——流道上游控制体相

的密度；

——流道中相

速度；

——流道上游控制体相

的比焓。

Images