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

摘要

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

说明书

技术领域

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

背景技术

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

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

发明内容

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

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

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

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

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

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

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

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

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

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

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

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

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

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

式中：

j——第j个流道；

d——流道的上游；

r——流道的下游；

t——时间；

L

ΔP

K

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

式中：

ρ

ρ

g——重力加速度；

Z

第二部分为控制体k中水池面高程z

式中：

ρ

ρ

g——重力加速度；

Z

第三部分是各相

式中：

g——重力加速度；

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

式中：

将重力压头

式中：

Δ(ρgΔz)

Δm

g——重力加速度；

S

i——控制体i；

k——控制体k；

Z

Z

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

式中：

g——重力加速度；

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

S——控制体底面积；

Δt——时间步长；

σ

A

α

v

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

式中：

P

P

Δt——时间步长；

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

式中：

t——时间；

P——总压力；

a——等于

b——等于

R——理想气体常数；

m

M

ρ

m

ρ

a

c

T

V

u

m

ρ

c

T

u

β——水的膨胀系数；

σ

A

α

α

v

v

δm

δm

δm

δm

δU

δU

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

式中：

t——时间；

Δt——时间步长；

P——总压力；

a——等于

b——等于

(ρgΔz)

(ρgΔz)

Δ(ρgΔz)

S——控制体底面积；

ΔP

R——理想气体常数；

m

M

ρ

m

ρ

a

c

T

V

u

m

ρ

c

T

u

β——水的膨胀系数；

σ

A

L

K

α

α

v

v

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

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

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

式中：

Δt——时间步长；

v

v

v

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

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

F(v

F(v

Δt——时间步长；

v

v

v

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

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

其中：

式中：

J

I——单位矩阵；

λ

v

F

Δt——时间步长；

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

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

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

式中：

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

式中：

如果流道中没有双向流动，则流道的空泡数α

式中：

α

v

v

ρ

ρ

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

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

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

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

式中：

Δt——时间步长；

σ

A

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

附图说明

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

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

具体实施方式

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

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

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

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

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

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

流道参数：流道长度L，阻力系数K，面积A，上游、下游开口标高Z、连接处高程Z

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

式中：

j——第j个流道；

d——流道的上游；

r——流道的下游；

t——时间；

L

ΔP

K

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

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

式中：

ρ

ρ

g——重力加速度；

Z

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

式中：

ρ

ρ

g——重力加速度；

Z

第三部分是各相

式中：

g——重力加速度；

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

式中：

将重力压头

式中：

Δ(ρgΔz)

Δm

g——重力加速度；

S

i——控制体i；

k——控制体k；

Z

Z

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

式中：

g——重力加速度；

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

S——控制体底面积；

Δt——时间步长；

σ

A

α

v

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

压力包括两部分，一部分是蒸汽分压P

式中：

P——控制体总压力；

P

R——理想气体常数；

T

V

m

M

m

x

P

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

式中：

P——控制体总压力；

t——时间；

R——理想气体常数；

P

ρ

m

T

V

m

M

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

t——时间；

V

m

ρ

β——水的膨胀系数；

ρ

T

a

P——压力；

式中dρ

式中：

t——时间；

m

ρ

β——水的膨胀系数；

T

a

P——总压力；

记

式中：

a

T

R——理想气体常数；

M

ρ

P——总压力；

P

ρ

z

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

式中：

t——时间；

m

ρ

β——水的膨胀系数；

T

a

P——总压力；

T

V

R——理想气体常数；

m

M

P——总压力；

m

ρ

P

ρ

a

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

式中：

t——时间；

m

ρ

β——水的膨胀系数；

T

a

P——总压力；

T

V

R——理想气体常数；

m

M

P——总压力；

m

ρ

P

ρ

a

σ

A

α

α

v

v

进一步，记

式中：

m

T

V

R——理想气体常数；

M

m

M

m

z

P

ρ

ρ

a

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

t——时间；

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

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

m

ρ

β——水的膨胀系数；

T

a

P——总压力；

T

V

R——理想气体常数；

m

M

P——总压力；

m

ρ

P

ρ

a

σ

A

α

α

v

v

由于公式(17)中

式中：

t——时间；

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

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

m

ρ

β——水的膨胀系数；

T

a

P——总压力；

T

V

R——理想气体常数；

m

M

P——总压力；

m

ρ

P

ρ

a

σ

A

α

α

v

v

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

U

式中：

U

m

u

t——时间；

c

T

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

式中：

i——第i个控制体；

j——第j个流道；

d——流道的上游；

r——流道的下游；

σ

A

t——时间；

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

t——时间；

c

T

U

u

m

c

T

U

u

σ

A

α

α

v

v

δm

δm

δU

δU

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

式中：

t——时间；

P——总压力；

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

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

R——理想气体常数；

m

M

ρ

m

ρ

a

c

T

V

u

m

ρ

c

T

u

β——水的膨胀系数；

σ

A

α

α

v

v

δm

δm

δm

δm

δU

δU

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

式中：

P

P

Δt——时间步长；

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

式中：

t——时间；

Δt——时间步长；

P——总压力；

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

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

(ρgΔz)

(ρgΔz)

Δ(ρgΔz)

S——控制体底面积；

ΔP

R——理想气体常数；

m

M

ρ

m

ρ

a

c

T

V

u

m

ρ

c

T

u

β——水的膨胀系数；

σ

A

L

K

α

α

v

v

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

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

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

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

式中：

Δt——时间步长；

v

v

v

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

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

F(v

F(v

Δt——时间步长；

v

v

v

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

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

其中：

式中：

J

I——单位矩阵；

λ

v

F

Δt——时间步长；

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

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

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

式中：

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

式中：

如果流道中没有双向流动，则流道的空泡数α

式中：

α

v

v

ρ

ρ

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

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

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

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

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

式中：

Δt——时间步长；

σ

A

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

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

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