一种严重事故下安全壳内控制体热力学响应综合分析方法

摘要

本发明涉及一种严重事故下安全壳内控制体热力学响应综合分析方法，包括：1、获取热工水力计算初始参数并初始化；2、根据计算时间控制参数判断是否结束计算；3、更新时间步长并求解边界条件；4、计算控制体内部总质量和总能量；5、判断控制体内部流体类型并分工况求解控制体内部状态参数；6、计算不可凝气体的状态参数及气相的平均温度；7、计算饱和水的闪蒸以及过饱和蒸汽的析出；8、计算气液两相与相界面之间的传热传质；9、迭代计算并输出各时刻状态参数，直至计算结束。本发明能准确模拟计算严重事故下安全壳内控制体在受到扰动后从不平衡状态向平衡状态过渡过程中的瞬态变化规律，如实再现实际相关的物理过程并反映系统的物理现象。

说明书

技术领域

本发明属于核安全分析技术领域，具体涉及一种严重事故下安全壳内控制体热力学响应综合分析方法。

背景技术

日本福岛核事故引发了社会上对反应堆严重事故的大量关注和研究，业界(包括核安全监管部门)对核电厂事故分析提出了新的要求，迫切需要加强机理模型和先进数值模拟技术的应用。安全壳热工水力响应分析作为核电厂反应堆严重事故进程分析的基础，精确计算其状态变化对于判断整体事故进程具有重要帮助。

当前常用的核电厂热工水力分析软件对于状态参数的求解多伴随控制方程的求解过程，同时受限于早期物性数据库的局限，特别是由于物性计算中偏导数计算数据的匮乏，常常采用半隐式的方式求解控制体中的温度、压力等状态参数，此种方式求解结果准确度偏低，且方程迭代过程容易不收敛造成方程求解失败。

发明内容

本发明的目的是针对现有技术中存在的缺陷，综合分析严重事故下安全壳内控制体热力学响应的变化规律，揭示实际相关的物理过程并反映系统的物理现象，提供一种严重事故下安全壳内控制体热力学响应综合分析方法，从而有效提高计算结果的准确度，避免因方程发散引起的求解失败问题。

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

一种严重事故下安全壳内控制体热力学响应综合分析方法，包括如下步骤：

步骤1、获取控制体热工水力计算初始参数以及计算所需的时间控制参数、边界条件参数，完成初始化；

步骤2、根据计算时间控制参数判断是否结束计算运行，当不满足结束条件时，进入步骤3继续执行计算，满足结束条件时，输出结果信息并结束计算；

步骤3、更新时间步长，根据边界条件参数计算出各时间内边界源项与指定控制体之间进行的能量与质量交换总额；

步骤4、根据质量和能量守恒方程，计算控制体内部总质量和总能量；

步骤5、判断控制体内部流体类型并分工况求解控制体内部状态参数；

步骤6、计算不可凝气体(如氮气、氧气、二氧化碳、氢气、一氧化碳等)状态参数及气相平均物性，相关物性数据可根据公开数据库和常用计算模型计算得到；

步骤7、计算饱和水的闪蒸以及过饱和蒸汽的析出；

步骤8、计算气液两相与相界面之间的传热传质；

步骤9、返回步骤2反复迭代计算。

进一步，如上所述的严重事故下安全壳内控制体热力学响应综合分析方法，步骤1中所述控制体热工水力的初始参数包括体积、标高、高度、水力直径、压力、气相初始温度、液相初始温度、各气体的体积份额、各组分的质量和气液两相的能量；

时间控制参数包括计算起始时间、终止时间、各时间段的最大、最小时间步长等；

边界条件参数包括外部源项的开启时间与关闭时间、作用控制体的编号、源项的状态、种类、组分、温度、压力、流量等。

进一步，如上所述的严重事故下安全壳内控制体热力学响应综合分析方法，步骤4中将微分、集总参数形式的质量守恒方程化为差分形式的方程，考虑质量源项后可得控制体内工质i的差分形式质量守恒方程：

式中：

f

d--流道的上游；

Δt--时间步长，s；

假设控制体内流体为理想流体，能量方程中不再考虑由于粘性作用引起的机械能转换为热能的部分，微分、集总参数形式的能量守恒方程化为差分形式的方程，可得：

式中：

f

d--流道的上游；

Δt--时间步长，s；

进一步，如上所述的严重事故下安全壳内控制体热力学响应综合分析方法，步骤5求解控制体状态参数时，首先判断控制体内部流体类型，根据不同工况选择选取不同的计算方程：

式中：

p--控制体压力，kPa；

p

m

R--摩尔气体常数，J/mol·K；

T

M

V--控制体体积，m

m

ρ

V--控制体体积，m

m

ρ

m

ρ

U

m

u

m

u

U

U

m

u

联立上述公式(3)-公式(6)并作转化，记方程为

分别对p，p

导数矩阵中第1行第1/3/4列元素分别为：

利用以上各式，根据牛顿迭代法，可得关于第k步的迭代方程如下：

上述方程中，除已知参数控制体总体积以及各工质的质量和气液两相的内能外，其它参数均通过当前状态下控制体内的温度和压力状态参数导出，不依赖于其它状态或外部参数，为全隐式形式的求解方程，通过对上述方程的求解，便可得到该条件下控制体内各相的温度和压力状态参数。

进一步，如上所述的严重事故下安全壳内控制体热力学响应综合分析方法，步骤7中，表示闪蒸质量时通过系统压力平衡与能量平衡两种方式进行计算，其中压力平衡使用的计算模型表示的闪蒸质量表示为：

式中：

m

m

m

M

V

p

p--控制体内部压力，kPa；

Z--压缩因子；

R--摩尔气体常数，J/mol·K；

T

采用能量平衡的方式计算时，认为闪蒸前后系统的能量保持平衡，根据能量守恒定律的表示形式：

U

式中：

U

m

m

h

h

转化后闪蒸质量表示为：

为保证计算结果的科学性与可靠性，采用上述两种计算结果的最小值作为最终的闪蒸质量；

当系统中蒸汽某时刻下对应的比内能小于其饱和比内能时，蒸汽将发生析出，根据能量平衡规律可以表示为：

m

式中：

m

u

m

u

u

转化后析出质量最终表示为：

进一步，如上所述的严重事故下安全壳内控制体热力学响应综合分析方法，步骤8中计算气液两相与相界面之间的传热传质时，首先考虑气相及液相到相界面的对流传热量，当相界面与气液两相存在温差时，对流传热量可根据牛顿冷却公式进行计算：

Q

Q

式中：

Q

ht

A--相界面面积，m

T

T

Q

ht

T

相界面的能量平衡方程写为：

Q

式中：

Q

Q

Γ--相变传质率，kg/s；

Δh--相变潜热，J/kg；

此外假定相变传质率是由气相和相界面之间的蒸汽质量扩散率决定的，因而根据热质比拟理论：

式中：

Γ--相变传质率，kg/s；

A--相界面面积，m

ht

λ--气相热导率，W/m·K；

D

M

p--气体总压，kPa；

R--摩尔气体常数，J/mol·K；

T

p

p

p

Sc--施密特数；

Pr--普朗特数；

在确定相界面温度后，即可得到气液相之间的传质率，最终获得气液两相之间的传热总量。

本发明的有益效果如下：本发明在核电厂热工水力计算软件开发过程中，根据集总参数模型中质量、能量方程构建本构关系式，提出了一种控制体中流体状态参数全隐式计算模型和方法。该方法将安全壳隔间按一定的规则等效划分成若干个控制体，从质量守恒方程和能量守恒方程出发，计算控制体内部总质量和总能量，经判断控制体内部流体类型后，分工况求解控制体内部状态参数，计算不可凝气体的状态参数及气相的平均温度，同时考虑饱和水的闪蒸、过饱和蒸汽的析出以及气液两相与相界面之间的传热传质，描述控制体在受到扰动后从不平衡状态向平衡状态过渡过程中的瞬态变化规律。本发明的提出能够有效提高计算结果的准确度，避免因方程发散引起的求解失败问题。

附图说明

图1为具体实施例中严重事故下安全壳内控制体热力学响应综合分析方法的流程图。

图2为具体实施例中控制体模型示意图。

具体实施方式

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

如图1所示，本发明提供了一种严重事故下安全壳内控制体热力学响应综合分析方法，步骤如下：

步骤1：获取控制体热工水力计算初始参数以及计算过程所需的控制参数、边界条件相关参数，并完成初始化，主要包括以下具体参数：

控制体参数：如图2中的控制体模型示意图所示，包括控制体的体积、标高、高度、水力直径、压力、气相初始温度、液相初始温度、各气体的体积份额、各组分的质量和气液两相的能量等；

时间控制参数：计算起始时间、终止时间、各时间段的最大、最小时间步长等。

边界条件参数包括外部源项的开启时间与关闭时间、作用控制体的编号、源项的状态、种类、组分、温度、压力、流量等。

步骤2：根据时间控制参数中的计算终止时间，判断是否结束计算运行，当计算时间小于给定的终止时间时，不满足结束条件，进入步骤3继续执行计算，直到达到计算终止时间，输出结果信息并结束计算。

步骤3：更新时间步长并求解边界条件，确定外部源项的状态、种类、组分组成及作用时间，例如，外部源项可以为过热蒸汽、不可凝气体、液态水或者三者之间的任意组合，对于给定源项，需进一步确定其作用时间、压力、气体温度、液相温度、气相流量、液相流量、各组分气体的份额，如相应参数缺省则作0处理，从而计算出各时间内边界源项与指定控制体之间进行的能量与质量交换总额。

步骤4：根据质量和能量守恒方程，计算控制体内部总质量和总能量，其中控制体内微分形式的质量守恒方程表示为：

式中：

ρ--控制体内流体密度，kg/m

t--时间，s；

Γ--质量源项的体密度，kg/m

微分形式的质量守恒方程在控制体内积分后，将微分、集总参数形式的质量守恒方程化为差分形式的方程，考虑质量源项后可得控制体内工质i的差分形式质量守恒方程：

式中：

f

Δt--时间步长，s；

上式右边第二项表示一个时间步长内通过所有与控制体相连的流道流入或流出控制体的工质i的质量，此项为正表示流体流入控制体，为负表示流体流出控制体；右边第三项表示控制体内非流动质量源项，此项包括控制体内部两相之间由于相变(闪蒸、析出与蒸发)引起的相间质量变化以及用户定义的外部源项；

控制体内微分形式的能量守恒方程为：

式中：

ρ--控制体内流体密度，kg/m

u--控制体内流体速度，m/s；

t--时间，s；

p--控制体压力，kPa

λ--流体的热导率，W/m·K；

T--控制体内流体温度，K；

S

上式右边第一项表示压强对流体微元体压缩所做的功，并假设控制体内流体为理想流体，能量方程中不再考虑由于粘性作用引起的机械能转换为热能的部分，将上述方程在控制体内积分后，将微分、集总参数形式的能量守恒方程化为差分形式的方程，可得：

式中：

f

Δt--时间步长，s；

上式右边第二项表示通过所有与控制体相连的流道流入或流出控制体的工质i的总焓，第三项表示控制体内非流动能量源项，此项包含因非流动质量源项引起的能量变化。

步骤5：求解控制体状态参数，首先判断控制体内部流体类型，如图1所示，控制体内部流体类别分为四种工况，分别是工况一：不可凝气体、水蒸气、水三者混合物；工况二：不可凝气体、水二者混合物；工况三：不可凝气体、水蒸气二者混合物；工况四：不可凝气体；以最复杂的工况一为例进行求解过程说明，其他工况下的求解过程将求解方程此进行简化即可；

控制体中的压力可以表示为：

式中：

p--控制体压力，kPa；

p

m

R--摩尔气体常数，J/mol·K；

T

M

V--控制体体积，m

m

ρ

控制体的体积可以表示为：

V--控制体体积，m

m

ρ

m

ρ

U

m

u

m

u

控制体内部液相总能量可以表示为：

U

U

m

u

联立上述5-8关系式并作转化，记方程为

分别对p，p

式中：c

上述公式中关于水和水蒸气的各导数项均为压力和温度的函数。

利用以上各式，根据牛顿迭代法，可得关于第k步的迭代方程如下：

上述方程中，除已知参数控制体总体积以及各工质的质量和气液两相的内能外，其它参数均通过当前状态下控制体内的温度和压力状态参数导出，不依赖于其它状态或外部参数，为全隐式形式的求解方程，通过对上述方程的求解，便可得到该条件下控制体内各相的各个时刻的温度和压力状态参数。

步骤6：计算不可凝气体(如氮气、氧气、二氧化碳、氢气、一氧化碳等)状态参数及气相平均物性，包括各组分气体的密度、分压、摩尔份额以及不可凝气体和水蒸气的混合物物性参数，相关物性数据可根据公开数据库和常用计算模型计算得到。

步骤7：过饱和水的闪蒸和过饱和蒸汽的析出，表示闪蒸质量时通过系统压力平衡与能量平衡两种方式进行计算，其中压力平衡使用的计算模型表示为：

发生闪蒸之前：

发生闪蒸以后：

式中：

p--控制体内部压力，kPa；

V

R--摩尔气体常数，J/mol·K；

T

m

M

p

m

Z--压缩因子，用来衡量蒸汽偏离理想气体的程度，Z的引入是由于上述推导过程是将蒸汽做理想气体的假设来处理的：

式中：

p

M

v

R--摩尔气体常数，J/mol·K；

T

闪蒸质量表示为：

式中：

m

m

m

M

V

p

p--控制体内部压力，kPa；

Z--压缩因子；

R--摩尔气体常数，J/mol·K；

T

采用能量平衡的方式计算时，认为闪蒸前后系统的能量保持平衡，根据能量守恒定律的表示形式：

U

式中：

U

m

m

h

h

转化后闪蒸质量表示为：

在本分析方法中，为保证计算结果的科学性与可靠性，采用上述两种计算结果的最小值作为最终的闪蒸质量；

当系统中蒸汽某时刻下对应的比内能小于其饱和比内能时，蒸汽将发生析出，根据能量平衡规律可以表示为：

m

式中：

m

u

m

u

u

转化后析出质量最终表示为：

步骤8：计算气液两相界面间传热传质，首先考虑气相及液相到相界面的对流传热量，当相界面与气液两相存在温差时，对流传热量可根据牛顿冷却公式进行计算：

Q

Q

式中：

Q

ht

A--相界面面积，m

T

T

Q

ht

T

其中，对流传热系数的计算方式可以根据水平面对流传热实验关联式，获得努塞尔数Nu后求得相应的对流传热系数：

强制对流：

Nu

自然对流：

Nu

Nu

式中：

Nu

Re--雷诺数；

Pr--普朗特数；

Nu

Gr--格拉晓夫数；

T

T

相界面的能量平衡方程写为：

Q

式中：

Q

Q

Γ--相变传质率，kg/s；

Δh--相变潜热，J/kg；

此处使用的相变潜热是相变潜热与相变前后气/液相达到饱和温度时消耗/放出的热量的叠加，而非常规定义的汽化潜热，此外还假定相变传质率是由气相和相界面之间的蒸汽质量扩散率决定的，因而根据热质比拟理论：

式中：

Γ--相变传质率，kg/s；

A--相界面面积，m

ht

λ--气相热导率，W/m·K；

D

M

p--气体总压，kPa；

R--摩尔气体常数，J/mol·K；

T

p

p

p

Sc--施密特数；

Pr--普朗特数；

在确定相界面温度后，即可得到气液相之间的传质率，最终获得气液两相之间的传热总量，实际在求解过程中，首先应假定一个相界面温度T

步骤9：返回步骤2反复迭代并判断是否达到计算结束条件，在满足结束条件后，输出计算结果，终止计算。

本方法从输入的安全壳计算相关参数着手，根据步骤1-9进行求解计算，最终得到各个时刻的安全壳控制体内的温度、压力、质量、能量等各项状态参数。这些求解的状态参数是进行安全壳安全分析的重要热工参数，是进行核电厂安全分析、严重事故计算基础条件。因此，作为求解热工参数的基础方法，本方法可应用于集总参数法的安全壳热工软件、严重事故计算分析软件中，用于安全壳内状态参数的求解。

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

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