CN112613158B - 一种严重事故下安全壳内控制体热力学响应综合分析方法 - Google Patents
一种严重事故下安全壳内控制体热力学响应综合分析方法 Download PDFInfo
- Publication number
- CN112613158B CN112613158B CN202011351583.4A CN202011351583A CN112613158B CN 112613158 B CN112613158 B CN 112613158B CN 202011351583 A CN202011351583 A CN 202011351583A CN 112613158 B CN112613158 B CN 112613158B
- Authority
- CN
- China
- Prior art keywords
- mass
- phase
- gas
- energy
- controlling
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000004458 analytical method Methods 0.000 title claims abstract description 22
- 230000004044 response Effects 0.000 title claims abstract description 21
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 claims abstract description 89
- 239000012071 phase Substances 0.000 claims abstract description 82
- 238000004364 calculation method Methods 0.000 claims abstract description 65
- 238000012546 transfer Methods 0.000 claims abstract description 55
- 239000007791 liquid phase Substances 0.000 claims abstract description 40
- 238000000034 method Methods 0.000 claims abstract description 37
- 230000005501 phase interface Effects 0.000 claims abstract description 36
- 230000008859 change Effects 0.000 claims abstract description 27
- 239000012530 fluid Substances 0.000 claims abstract description 24
- 238000001704 evaporation Methods 0.000 claims abstract description 21
- 230000008020 evaporation Effects 0.000 claims abstract description 21
- 229920006395 saturated elastomer Polymers 0.000 claims abstract description 13
- 238000001556 precipitation Methods 0.000 claims abstract description 9
- 239000007789 gas Substances 0.000 claims description 115
- 239000007788 liquid Substances 0.000 claims description 26
- 238000011144 upstream manufacturing Methods 0.000 claims description 13
- 238000004134 energy conservation Methods 0.000 claims description 10
- 230000000704 physical effect Effects 0.000 claims description 9
- 238000009792 diffusion process Methods 0.000 claims description 8
- 239000000203 mixture Substances 0.000 claims description 8
- 238000006243 chemical reaction Methods 0.000 claims description 6
- 239000011159 matrix material Substances 0.000 claims description 6
- 230000009471 action Effects 0.000 claims description 5
- 238000007906 compression Methods 0.000 claims description 5
- 238000001816 cooling Methods 0.000 claims description 3
- 230000000694 effects Effects 0.000 claims description 3
- 239000011555 saturated liquid Substances 0.000 claims description 3
- 239000012808 vapor phase Substances 0.000 claims description 2
- 230000008569 process Effects 0.000 abstract description 14
- 230000007704 transition Effects 0.000 abstract description 6
- 230000001052 transient effect Effects 0.000 abstract description 2
- 239000000126 substance Substances 0.000 description 5
- IJGRMHOSHXDMSA-UHFFFAOYSA-N Atomic nitrogen Chemical compound N#N IJGRMHOSHXDMSA-UHFFFAOYSA-N 0.000 description 4
- 238000001727 in vivo Methods 0.000 description 4
- 238000010586 diagram Methods 0.000 description 3
- UGFAIRIUMAVXCW-UHFFFAOYSA-N Carbon monoxide Chemical compound [O+]#[C-] UGFAIRIUMAVXCW-UHFFFAOYSA-N 0.000 description 2
- UBAZGMLMVVQSCD-UHFFFAOYSA-N carbon dioxide;molecular oxygen Chemical compound O=O.O=C=O UBAZGMLMVVQSCD-UHFFFAOYSA-N 0.000 description 2
- 229910002091 carbon monoxide Inorganic materials 0.000 description 2
- 230000006835 compression Effects 0.000 description 2
- 239000001257 hydrogen Substances 0.000 description 2
- 229910052739 hydrogen Inorganic materials 0.000 description 2
- 125000004435 hydrogen atom Chemical class [H]* 0.000 description 2
- 229910052757 nitrogen Inorganic materials 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000009795 derivation Methods 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 238000002474 experimental method Methods 0.000 description 1
- 230000016507 interphase Effects 0.000 description 1
- 230000007246 mechanism Effects 0.000 description 1
- 238000012545 processing Methods 0.000 description 1
- 230000003014 reinforcing effect Effects 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 238000004088 simulation Methods 0.000 description 1
- 238000009834 vaporization Methods 0.000 description 1
- 230000008016 vaporization Effects 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
- G06F2119/08—Thermal analysis or thermal optimisation
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
- G06F2119/12—Timing analysis or timing optimisation
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
- G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces
-
- Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
- Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
- Y02E—REDUCTION OF GREENHOUSE GAS [GHG] EMISSIONS, RELATED TO ENERGY GENERATION, TRANSMISSION OR DISTRIBUTION
- Y02E30/00—Energy generation of nuclear origin
- Y02E30/30—Nuclear fission reactors
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Computer Hardware Design (AREA)
- Evolutionary Computation (AREA)
- Geometry (AREA)
- General Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Feedback Control In General (AREA)
Abstract
本发明涉及一种严重事故下安全壳内控制体热力学响应综合分析方法,包括:1、获取热工水力计算初始参数并初始化;2、根据计算时间控制参数判断是否结束计算;3、更新时间步长并求解边界条件;4、计算控制体内部总质量和总能量;5、判断控制体内部流体类型并分工况求解控制体内部状态参数;6、计算不可凝气体的状态参数及气相的平均温度;7、计算饱和水的闪蒸以及过饱和蒸汽的析出;8、计算气液两相与相界面之间的传热传质;9、迭代计算并输出各时刻状态参数,直至计算结束。本发明能准确模拟计算严重事故下安全壳内控制体在受到扰动后从不平衡状态向平衡状态过渡过程中的瞬态变化规律,如实再现实际相关的物理过程并反映系统的物理现象。
Description
技术领域
本发明属于核安全分析技术领域,具体涉及一种严重事故下安全壳内控制体热力学响应综合分析方法。
背景技术
日本福岛核事故引发了社会上对反应堆严重事故的大量关注和研究,业界(包括核安全监管部门)对核电厂事故分析提出了新的要求,迫切需要加强机理模型和先进数值模拟技术的应用。安全壳热工水力响应分析作为核电厂反应堆严重事故进程分析的基础,精确计算其状态变化对于判断整体事故进程具有重要帮助。
当前常用的核电厂热工水力分析软件对于状态参数的求解多伴随控制方程的求解过程,同时受限于早期物性数据库的局限,特别是由于物性计算中偏导数计算数据的匮乏,常常采用半隐式的方式求解控制体中的温度、压力等状态参数,此种方式求解结果准确度偏低,且方程迭代过程容易不收敛造成方程求解失败。
发明内容
本发明的目的是针对现有技术中存在的缺陷,综合分析严重事故下安全壳内控制体热力学响应的变化规律,揭示实际相关的物理过程并反映系统的物理现象,提供一种严重事故下安全壳内控制体热力学响应综合分析方法,从而有效提高计算结果的准确度,避免因方程发散引起的求解失败问题。
为了实现上述目的,本发明提供的技术方案如下:
一种严重事故下安全壳内控制体热力学响应综合分析方法,包括如下步骤:
步骤1、获取控制体热工水力计算初始参数以及计算所需的时间控制参数、边界条件参数,完成初始化;
步骤2、根据计算时间控制参数判断是否结束计算运行,当不满足结束条件时,进入步骤3继续执行计算,满足结束条件时,输出结果信息并结束计算;
步骤3、更新时间步长,根据边界条件参数计算出各时间内边界源项与指定控制体之间进行的能量与质量交换总额;
步骤4、根据质量和能量守恒方程,计算控制体内部总质量和总能量;
步骤5、判断控制体内部流体类型并分工况求解控制体内部状态参数;
步骤6、计算不可凝气体(如氮气、氧气、二氧化碳、氢气、一氧化碳等)状态参数及气相平均物性,相关物性数据可根据公开数据库和常用计算模型计算得到;
步骤7、计算饱和水的闪蒸以及过饱和蒸汽的析出;
步骤8、计算气液两相与相界面之间的传热传质;
步骤9、返回步骤2反复迭代计算。
进一步,如上所述的严重事故下安全壳内控制体热力学响应综合分析方法,步骤1中所述控制体热工水力的初始参数包括体积、标高、高度、水力直径、压力、气相初始温度、液相初始温度、各气体的体积份额、各组分的质量和气液两相的能量;
时间控制参数包括计算起始时间、终止时间、各时间段的最大、最小时间步长等;
边界条件参数包括外部源项的开启时间与关闭时间、作用控制体的编号、源项的状态、种类、组分、温度、压力、流量等。
进一步,如上所述的严重事故下安全壳内控制体热力学响应综合分析方法,步骤4中将微分、集总参数形式的质量守恒方程化为差分形式的方程,考虑质量源项后可得控制体内工质i的差分形式质量守恒方程:
式中:
--控制体内工质i在新时刻的质量,kg;
--控制体内工质i在旧时刻的质量,kg;
fj--经流道j流入到控制体的体积流量,m3/s;
--流道j上游工质i的密度,kg/m3;
d--流道的上游;
Δt--时间步长,s;
--控制体内工质i非流动质量源项变化率,kg/s;
假设控制体内流体为理想流体,能量方程中不再考虑由于粘性作用引起的机械能转换为热能的部分,微分、集总参数形式的能量守恒方程化为差分形式的方程,可得:
式中:
--控制体内工质i在新时刻的内能,J;
--控制体内工质i在旧时刻的内能,J;
fj--经流道j流入到控制体的体积流量,m3/s;
--流道j上游工质i的密度,kg/m3;
--流道j上游控制体内工质i的比焓,J/kg;
d--流道的上游;
Δt--时间步长,s;
--控制体内工质i非流动能量源项变化率,J/s。
进一步,如上所述的严重事故下安全壳内控制体热力学响应综合分析方法,步骤5求解控制体状态参数时,首先判断控制体内部流体类型,根据不同工况选择选取不同的计算方程:
式中:
p--控制体压力,kPa;
ps--水蒸气分压,kPa;
mi--第i种不可凝气体的质量,kg;
R--摩尔气体常数,J/mol·K;
Tg--气体温度,K;
Mi--第i种不可凝气体的摩尔质量,kg/mol;
V--控制体体积,m3;
mw--控制体内液态水的质量,kg;
ρw--水的密度,是压力p和液态水温度Tw的函数,kg/m3;
V--控制体体积,m3;
mw--控制体内液态水的质量,kg;
ρw--水的密度,是控制体压力p和液态水温度Tw的函数,kg/m3;
ms--控制体内水蒸气的质量,kg;
ρs--水蒸气的密度,是水蒸气分压力ps和气体温度Tg的函数,kg/m3;
Ug--控制体内所有气体的总内能,J;
ms--水蒸气的质量,kg;
us--水蒸气的比内能,是水蒸气分压力ps和气体温度Tg的函数,J/kg;
mi--控制体内第i种不可凝气体的质量,kg;
ui--第i种不可凝气体的比内能,是气体温度Tg的函数,J/kg;
Ul=mw·uw(p,Tw) (6)
Ul--控制体内水的总内能,J;
mw--水的质量,kg;
uw--水的比内能,是控制体压力p和液态水温度Tw的函数,J/kg;
联立上述公式(3)-公式(6)并作转化,记方程为
分别对p,ps,Tg,Tw进行偏导数计算,则其导数矩阵表示为:
导数矩阵中第1行第1/3/4列元素分别为:
利用以上各式,根据牛顿迭代法,可得关于第k步的迭代方程如下:
上述方程中,除已知参数控制体总体积以及各工质的质量和气液两相的内能外,其它参数均通过当前状态下控制体内的温度和压力状态参数导出,不依赖于其它状态或外部参数,为全隐式形式的求解方程,通过对上述方程的求解,便可得到该条件下控制体内各相的温度和压力状态参数。
进一步,如上所述的严重事故下安全壳内控制体热力学响应综合分析方法,步骤7中,表示闪蒸质量时通过系统压力平衡与能量平衡两种方式进行计算,其中压力平衡使用的计算模型表示的闪蒸质量表示为:
式中:
mfl--闪蒸质量,kg;
m2--闪蒸后控制体内部水蒸气的质量,kg;
m1--闪蒸前控制体内部水蒸气的质量,kg;
Mw--水的摩尔质量,kg/mol;
Vg--气相体积,m3;
psat--控制体内部饱和蒸汽压,kPa;
p--控制体内部压力,kPa;
Z--压缩因子;
R--摩尔气体常数,J/mol·K;
Tg--气体温度,K;
采用能量平衡的方式计算时,认为闪蒸前后系统的能量保持平衡,根据能量守恒定律的表示形式:
Usat=(mW-mfl)hW(p)+mflhs(p) (14)
式中:
Usat--系统中工质的总能量,J;
mW--系统中液相水的质量,kg;
mfl--闪蒸质量,kg;
hw--饱和水的比焓值,是压力p的函数,J/kg;
hs--饱和水蒸气的比焓值,是压力p的函数,J/kg;
转化后闪蒸质量表示为:
为保证计算结果的科学性与可靠性,采用上述两种计算结果的最小值作为最终的闪蒸质量;
当系统中蒸汽某时刻下对应的比内能小于其饱和比内能时,蒸汽将发生析出,根据能量平衡规律可以表示为:
msus=mrainoutuw(ps)+(ms-mrainout)us(ps) (16)
式中:
ms--过饱和蒸汽的质量,kg;
us--蒸汽的比内能,J/kg;
mrainout--析出的水的质量,kg;
uw(ps)--蒸汽分压下对应的饱和液相比内能,J/kg;
us(ps)--蒸汽分压下对应的蒸汽比内能,J/kg;
转化后析出质量最终表示为:
进一步,如上所述的严重事故下安全壳内控制体热力学响应综合分析方法,步骤8中计算气液两相与相界面之间的传热传质时,首先考虑气相及液相到相界面的对流传热量,当相界面与气液两相存在温差时,对流传热量可根据牛顿冷却公式进行计算:
Qg,i=htgA(Tg-Ti) (18)
Ql,i=htlA(Tw-Ti) (19)
式中:
Qg,i--气相对相界面的对流传热量,W;
htg--气相对流传热系数,W/m2·K;
A--相界面面积,m2;
Tg--气体温度,K;
Ti--相界面温度,K;
Ql,i--液相对相界面的对流传热量,W;
htl--液相对流传热系数,W/m2·K;
Tw--液相温度,K;
相界面的能量平衡方程写为:
Qg,i+Ql,i=ΓΔh (20)
式中:
Qg,i--气相对相界面的对流传热量,W;
Ql,i--液相对相界面的对流传热量,W;
Γ--相变传质率,kg/s;
Δh--相变潜热,J/kg;
此外假定相变传质率是由气相和相界面之间的蒸汽质量扩散率决定的,因而根据热质比拟理论:
式中:
Γ--相变传质率,kg/s;
--单位表面上质量变化率,kg/m2·s;
A--相界面面积,m2;
htg--气相对流传热系数,W/m2·K;
λ--气相热导率,W/m·K;
Ds--蒸汽的扩散系数,m2/s;
Ms--蒸汽摩尔质量,kg/mol;
p--气体总压,kPa;
R--摩尔气体常数,J/mol·K;
Tg--气体温度,K;
ps,i--蒸汽在相界面处的分压,kPa;
pstm,b--蒸汽在气相主流处的分压,kPa;
pNCG,avg--扩散边界层中不可凝气体平均分压,kPa;
Sc--施密特数;
Pr--普朗特数;
在确定相界面温度后,即可得到气液相之间的传质率,最终获得气液两相之间的传热总量。
本发明的有益效果如下:本发明在核电厂热工水力计算软件开发过程中,根据集总参数模型中质量、能量方程构建本构关系式,提出了一种控制体中流体状态参数全隐式计算模型和方法。该方法将安全壳隔间按一定的规则等效划分成若干个控制体,从质量守恒方程和能量守恒方程出发,计算控制体内部总质量和总能量,经判断控制体内部流体类型后,分工况求解控制体内部状态参数,计算不可凝气体的状态参数及气相的平均温度,同时考虑饱和水的闪蒸、过饱和蒸汽的析出以及气液两相与相界面之间的传热传质,描述控制体在受到扰动后从不平衡状态向平衡状态过渡过程中的瞬态变化规律。本发明的提出能够有效提高计算结果的准确度,避免因方程发散引起的求解失败问题。
附图说明
图1为具体实施例中严重事故下安全壳内控制体热力学响应综合分析方法的流程图。
图2为具体实施例中控制体模型示意图。
具体实施方式
以下结合附图对本发明方法的具体实施方式作进一步详细说明。在具体实施方式部分中,各公式的编号可能与发明内容部分相同公式的编号不同,相应编号只代表在该部分的公式顺序情况。
如图1所示,本发明提供了一种严重事故下安全壳内控制体热力学响应综合分析方法,步骤如下:
步骤1:获取控制体热工水力计算初始参数以及计算过程所需的控制参数、边界条件相关参数,并完成初始化,主要包括以下具体参数:
控制体参数:如图2中的控制体模型示意图所示,包括控制体的体积、标高、高度、水力直径、压力、气相初始温度、液相初始温度、各气体的体积份额、各组分的质量和气液两相的能量等;
时间控制参数:计算起始时间、终止时间、各时间段的最大、最小时间步长等。
边界条件参数包括外部源项的开启时间与关闭时间、作用控制体的编号、源项的状态、种类、组分、温度、压力、流量等。
步骤2:根据时间控制参数中的计算终止时间,判断是否结束计算运行,当计算时间小于给定的终止时间时,不满足结束条件,进入步骤3继续执行计算,直到达到计算终止时间,输出结果信息并结束计算。
步骤3:更新时间步长并求解边界条件,确定外部源项的状态、种类、组分组成及作用时间,例如,外部源项可以为过热蒸汽、不可凝气体、液态水或者三者之间的任意组合,对于给定源项,需进一步确定其作用时间、压力、气体温度、液相温度、气相流量、液相流量、各组分气体的份额,如相应参数缺省则作0处理,从而计算出各时间内边界源项与指定控制体之间进行的能量与质量交换总额。
步骤4:根据质量和能量守恒方程,计算控制体内部总质量和总能量,其中控制体内微分形式的质量守恒方程表示为:
式中:
ρ--控制体内流体密度,kg/m3;
t--时间,s;
--速度向量,m/s;
Γ--质量源项的体密度,kg/m3·s
微分形式的质量守恒方程在控制体内积分后,将微分、集总参数形式的质量守恒方程化为差分形式的方程,考虑质量源项后可得控制体内工质i的差分形式质量守恒方程:
式中:
--控制体内工质i在新时刻的质量,kg;
--控制体内工质i在旧时刻的质量,kg;
fj--经流道j流入到控制体的体积流量,m3/s;
--流道上游工质i的密度,kg/m3;
Δt--时间步长,s;
--控制体内工质i非流动质量源项变化率,kg/s;
上式右边第二项表示一个时间步长内通过所有与控制体相连的流道流入或流出控制体的工质i的质量,此项为正表示流体流入控制体,为负表示流体流出控制体;右边第三项表示控制体内非流动质量源项,此项包括控制体内部两相之间由于相变(闪蒸、析出与蒸发)引起的相间质量变化以及用户定义的外部源项;
控制体内微分形式的能量守恒方程为:
式中:
ρ--控制体内流体密度,kg/m3;
u--控制体内流体速度,m/s;
t--时间,s;
--速度向量,m/s;
p--控制体压力,kPa
λ--流体的热导率,W/m·K;
T--控制体内流体温度,K;
Sh--流体的内热源,kg/m2·s2
上式右边第一项表示压强对流体微元体压缩所做的功,并假设控制体内流体为理想流体,能量方程中不再考虑由于粘性作用引起的机械能转换为热能的部分,将上述方程在控制体内积分后,将微分、集总参数形式的能量守恒方程化为差分形式的方程,可得:
式中:
--控制体内工质i在新时刻的内能,J;
--控制体内工质i在旧时刻的内能,J;
fj--经流道j流入到控制体的体积流量,m3/s;
--流道上游工质i的密度,kg/m3;
--流道j上游控制体内工质i的比焓,J/kg;
Δt--时间步长,s;
--控制体内工质i非流动能量源项变化率,J/s;
上式右边第二项表示通过所有与控制体相连的流道流入或流出控制体的工质i的总焓,第三项表示控制体内非流动能量源项,此项包含因非流动质量源项引起的能量变化。
步骤5:求解控制体状态参数,首先判断控制体内部流体类型,如图1所示,控制体内部流体类别分为四种工况,分别是工况一:不可凝气体、水蒸气、水三者混合物;工况二:不可凝气体、水二者混合物;工况三:不可凝气体、水蒸气二者混合物;工况四:不可凝气体;以最复杂的工况一为例进行求解过程说明,其他工况下的求解过程将求解方程此进行简化即可;
控制体中的压力可以表示为:
式中:
p--控制体压力,kPa;
ps--水蒸气分压,kPa;
mi--第i种不可凝气体的质量,kg;
R--摩尔气体常数,J/mol·K;
Tg--气体温度,K;
Mi--第i种不可凝气体的摩尔质量,kg/mol;
V--控制体体积,m3;
mw--控制体内液态水的质量,kg;
ρw--水的密度,是压力p和液态水温度Tw的函数,kg/m3;
控制体的体积可以表示为:
V--控制体体积,m3;
mw--控制体内液态水的质量,kg;
ρw--水的密度,是压力p和液态水温度Tw的函数,kg/m3;
ms--控制体内水蒸气的质量,kg;
ρs--水蒸气的密度,是水蒸气分压力ps和气体温度Tg的函数,kg/m3;控制体内部气相总能量可以表示为:
Ug--控制体内所有气体的总内能,J;
ms--水蒸气的质量,kg;
us--水蒸气的比内能,是水蒸气分压力ps和气体温度Tg的函数,J/kg;
mi--控制体内第i种不可凝气体的质量,kg;
ui--第i种不可凝气体的比内能,是气体温度Tg的函数,J/kg;
控制体内部液相总能量可以表示为:
Ul=mw·uw(p,Tw) (8)
Ul--控制体内水的总内能,J;
mw--水的质量,kg;
uw--水的比内能,是控制体压力p和液态水温度Tw的函数,J/kg;
联立上述5-8关系式并作转化,记方程为
分别对p,ps,Tg,Tw进行偏导数计算,则其导数矩阵表示为:
式中:cvi为第i种不可凝气体的定容比热,导数矩阵中第1行第1/3/4列元素分别为:
上述公式中关于水和水蒸气的各导数项均为压力和温度的函数。
利用以上各式,根据牛顿迭代法,可得关于第k步的迭代方程如下:
上述方程中,除已知参数控制体总体积以及各工质的质量和气液两相的内能外,其它参数均通过当前状态下控制体内的温度和压力状态参数导出,不依赖于其它状态或外部参数,为全隐式形式的求解方程,通过对上述方程的求解,便可得到该条件下控制体内各相的各个时刻的温度和压力状态参数。
步骤6:计算不可凝气体(如氮气、氧气、二氧化碳、氢气、一氧化碳等)状态参数及气相平均物性,包括各组分气体的密度、分压、摩尔份额以及不可凝气体和水蒸气的混合物物性参数,相关物性数据可根据公开数据库和常用计算模型计算得到。
步骤7:过饱和水的闪蒸和过饱和蒸汽的析出,表示闪蒸质量时通过系统压力平衡与能量平衡两种方式进行计算,其中压力平衡使用的计算模型表示为:
发生闪蒸之前:
发生闪蒸以后:
式中:
p--控制体内部压力,kPa;
Vg--气相体积,m3;
R--摩尔气体常数,J/mol·K;
Tg--气体温度,K;
m1--闪蒸前控制体内部水蒸气的质量,kg;
Mw--水的摩尔质量,kg/mol;
psat--控制体内部饱和蒸汽压,kPa;
m2--闪蒸后控制体内部水蒸气的质量,kg;
Z--压缩因子,用来衡量蒸汽偏离理想气体的程度,Z的引入是由于上述推导过程是将蒸汽做理想气体的假设来处理的:
式中:
ps--水蒸气分压,kPa;
Mw--水的摩尔质量,kg/mol;
vs--水蒸汽的比容,m3/kg;
R--摩尔气体常数,J/mol·K;
Tg--气体温度,K;
闪蒸质量表示为:
式中:
mfl--闪蒸质量,kg;
m2--闪蒸后控制体内部水蒸气的质量,kg;
m1--闪蒸前控制体内部水蒸气的质量,kg;
Mw--水的摩尔质量,kg/mol;
Vg--气相体积,m3;
psat--控制体内部饱和蒸汽压,kPa;
p--控制体内部压力,kPa;
Z--压缩因子;
R--摩尔气体常数,J/mol·K;
Tg--气体温度,K;
采用能量平衡的方式计算时,认为闪蒸前后系统的能量保持平衡,根据能量守恒定律的表示形式:
Usat=(mW-mfl)hW(p)+mflhs(p) (19)
式中:
Usat--系统中工质的总能量,J;
mW--系统中液相水的质量,kg;
mfl--闪蒸质量,kg;
hw--饱和水的比焓值,是压力p的函数,J/kg;
hs--饱和水蒸气的比焓值,是压力p的函数,J/kg;
转化后闪蒸质量表示为:
在本分析方法中,为保证计算结果的科学性与可靠性,采用上述两种计算结果的最小值作为最终的闪蒸质量;
当系统中蒸汽某时刻下对应的比内能小于其饱和比内能时,蒸汽将发生析出,根据能量平衡规律可以表示为:
msus=mrainoutuw(ps)+(ms-mrainout)us(ps) (21)
式中:
ms--过饱和蒸汽的质量,kg;
us--蒸汽的比内能,J/kg;
mrainout--析出的水的质量,kg;
uw(ps)--蒸汽分压下对应的饱和液相比内能,J/kg;
us(ps)--蒸汽分压下对应的蒸汽比内能,J/kg;
转化后析出质量最终表示为:
步骤8:计算气液两相界面间传热传质,首先考虑气相及液相到相界面的对流传热量,当相界面与气液两相存在温差时,对流传热量可根据牛顿冷却公式进行计算:
Qg,i=htgA(Tg-Ti) (23)
Ql,i=htlA(Tw-Ti) (24)
式中:
Qg,i--气相对相界面的对流传热量,W;
htg--气相对流传热系数,W/m2·K;
A--相界面面积,m2;
Tg--气体温度,K;
Ti--相界面温度,K;
Ql,i--液相对相界面的对流传热量,W;
htl--液相对流传热系数,W/m2·K;
Tw--液相温度,K;
其中,对流传热系数的计算方式可以根据水平面对流传热实验关联式,获得努塞尔数Nu后求得相应的对流传热系数:
强制对流:
NuFC=0.037Re0.8Pr1/3 (25)
自然对流:
NuNC=0.54(Gr·Pr)1/4Tg>Tw (26)
NuNC=0.15(Gr·Pr)1/3Tg<Tw (27)
式中:
NuFC--强制对流努塞尔数;
Re--雷诺数;
Pr--普朗特数;
NuNC--自然对流努塞尔数;
Gr--格拉晓夫数;
Tg--气体温度,K;
Tw--液相温度,K;
相界面的能量平衡方程写为:
Qg,i+Ql,i=ΓΔh (28)
式中:
Qg,i--气相对相界面的对流传热量,W;
Ql,i--液相对相界面的对流传热量,W;
Γ--相变传质率,kg/s;
Δh--相变潜热,J/kg;
此处使用的相变潜热是相变潜热与相变前后气/液相达到饱和温度时消耗/放出的热量的叠加,而非常规定义的汽化潜热,此外还假定相变传质率是由气相和相界面之间的蒸汽质量扩散率决定的,因而根据热质比拟理论:
式中:
Γ--相变传质率,kg/s;
--单位表面上质量变化率,kg/m2·s;
A--相界面面积,m2;
htg--气相对流传热系数,W/m2·K;
λ--气相热导率,W/m·K;
Ds--蒸汽的扩散系数,m2/s;
Ms--蒸汽摩尔质量,kg/mol;
p--气体总压,kPa;
R--摩尔气体常数,J/mol·K;
Tg--气体温度,K;
ps,i--蒸汽在相界面处的分压,kPa;
pstm,b--蒸汽在气相主流处的分压,kPa;
pNCG,avg--扩散边界层中不可凝气体平均分压,kPa;
Sc--施密特数;
Pr--普朗特数;
在确定相界面温度后,即可得到气液相之间的传质率,最终获得气液两相之间的传热总量,实际在求解过程中,首先应假定一个相界面温度Ti并得到相应的传热量,然后根据上述公式计算出新的界面温度及传热量参数,当假定值的结果与计算值的结果在预先设置的偏差范围内时认为假定值是合理的,如果不满足条件则继续迭代更新,直至前后两次计算的结果满足偏差范围,最终得到界面温度的真实值;
步骤9:返回步骤2反复迭代并判断是否达到计算结束条件,在满足结束条件后,输出计算结果,终止计算。
本方法从输入的安全壳计算相关参数着手,根据步骤1-9进行求解计算,最终得到各个时刻的安全壳控制体内的温度、压力、质量、能量等各项状态参数。这些求解的状态参数是进行安全壳安全分析的重要热工参数,是进行核电厂安全分析、严重事故计算基础条件。因此,作为求解热工参数的基础方法,本方法可应用于集总参数法的安全壳热工软件、严重事故计算分析软件中,用于安全壳内状态参数的求解。
对于本领域技术人员而言,显然本发明方法不限于上述示范性实施例的细节,而且在不背离本发明的精神或基本特征的情况下,能够以其他的具体形式实现本发明方法。因此,无论从哪一点来看,均应将实施例看作是示范性的,而且是非限制性的,本发明方法的范围由所附权利要求而不是上述说明限定,因此旨在将落在权利要求的等同要件的含义和范围内的所有变化囊括在本发明方法内。不应将权利要求中的任何附图标记视为限制所涉及的权利要求。
此外,应当理解,虽然本说明书按照实施方式加以描述,但并非每个实施方式仅包含一个独立的技术方案,说明书的这种叙述方式仅仅是为清楚起见,本领域技术人员应当将说明书作为一个整体,各实施例中的技术方案也可以经适当组合,形成本领域技术人员可以理解的其他实施方式。
Claims (5)
1.一种严重事故下安全壳内控制体热力学响应综合分析方法,其特征在于,包括如下步骤:
步骤1、获取控制体热工水力计算初始参数以及计算所需的时间控制参数、边界条件参数,完成初始化;
步骤2、根据计算时间控制参数判断是否结束计算运行,当不满足结束条件时,进入步骤3继续执行计算,满足结束条件时,输出结果信息并结束计算;
步骤3、更新时间步长,根据边界条件参数计算出各时间内边界源项与指定控制体之间进行的能量与质量交换总额;
步骤4、根据质量和能量守恒方程,计算控制体内部总质量和总能量;
步骤5、判断控制体内部流体类型并分工况求解控制体内部状态参数;求解控制体内部状态参数时,首先判断控制体内部流体类型,根据不同工况选择选取不同的计算方程:
式中:
p--控制体压力;
ps--水蒸气分压;
mi--第i种不可凝气体的质量;
R--摩尔气体常数;
Tg--气体温度;
Mi--第i种不可凝气体的摩尔质量;
V--控制体体积;
mw--控制体内液态水的质量;
ρw--水的密度,是压力p和液态水温度Tw的函数;
V--控制体体积;
mw--控制体内液态水的质量;
ρw--水的密度,是控制体压力p和液态水温度Tw的函数;
ms--控制体内水蒸气的质量;
ρs--水蒸气的密度,是水蒸气分压力ps和气体温度Tg的函数;
Ug--控制体内所有气体的总内能;
ms--水蒸气的质量;
us--水蒸气的比内能,是水蒸气分压力ps和气体温度Tg的函数;
mi--控制体内第i种不可凝气体的质量;
ui--第i种不可凝气体的比内能,是气体温度Tg的函数;
Ul=mw·uw(p,Tw) (6)
Ul--控制体内水的总内能;
mw--水的质量;
uw--水的比内能,是控制体压力p和液态水温度Tw的函数;
联立上述公式(3)-公式(6)并作转化,记方程为
分别对p,ps,Tg,Tw进行偏导数计算,则其导数矩阵表示为:
cvi为控制体内第i种不可凝气体的定容比热;
导数矩阵中第1行第1/3/4列元素分别为:
利用以上各式,根据牛顿迭代法,可得关于第k步的迭代方程如下:
上述方程中,除已知参数控制体总体积以及各工质的质量和气液两相的内能外,其它参数均通过当前状态下控制体内的温度和压力状态参数导出,不依赖于其它状态或外部参数,为全隐式形式的求解方程,通过对上述方程的求解,便可得到该条件下控制体内各相的温度和压力状态参数;
步骤6、计算不可凝气体状态参数及气相平均物性;
步骤7、计算饱和水的闪蒸以及过饱和蒸汽的析出;
步骤8、计算气液两相与相界面之间的传热传质;
步骤9、返回步骤2反复迭代计算。
2.如权利要求1所述的严重事故下安全壳内控制体热力学响应综合分析方法,其特征在于,步骤1中所述控制体热工水力的初始参数包括体积、标高、高度、水力直径、压力、气相初始温度、液相初始温度、各气体的体积份额、各组分的质量和气液两相的能量;
时间控制参数包括计算起始时间、终止时间、各时间段的最大、最小时间步长;
边界条件参数包括外部源项的开启时间与关闭时间、作用控制体的编号、源项的状态、种类、组分、温度、压力、流量。
3.如权利要求1所述的严重事故下安全壳内控制体热力学响应综合分析方法,其特征在于,步骤4中将微分、集总参数形式的质量守恒方程化为差分形式的方程,考虑质量源项后可得控制体内工质i的差分形式质量守恒方程:
式中:
--控制体内工质i在新时刻的质量;
--控制体内工质i在旧时刻的质量;
fj--经流道j流入到控制体的体积流量;
--流道j上游工质i的密度;
d--流道的上游;
Δt--时间步长;
--控制体内工质i非流动质量源项变化率;
假设控制体内流体为理想流体,能量方程中不再考虑由于粘性作用引起的机械能转换为热能的部分,微分、集总参数形式的能量守恒方程化为差分形式的方程,可得:
式中:
--控制体内工质i在新时刻的内能;
--控制体内工质i在旧时刻的内能;
fj--经流道j流入到控制体的体积流量;
--流道j上游工质i的密度;
--流道j上游控制体内工质i的比焓;
d--流道的上游;
Δt--时间步长;
--控制体内工质i非流动能量源项变化率。
4.如权利要求1所述的严重事故下安全壳内控制体热力学响应综合分析方法,其特征在于,步骤7中,表示闪蒸质量时通过系统压力平衡与能量平衡两种方式进行计算,其中压力平衡使用的计算模型表示的闪蒸质量表示为:
式中:
mfl--闪蒸质量;
m2--闪蒸后控制体内部水蒸气的质量;
m1--闪蒸前控制体内部水蒸气的质量;
Mw--水的摩尔质量;
Vg--气相体积;
psat--控制体内部饱和蒸汽压;
p--控制体内部压力;
Z--压缩因子;
R--摩尔气体常数;
Tg--气体温度;
采用能量平衡的方式计算时,认为闪蒸前后系统的能量保持平衡,根据能量守恒定律的表示形式:
Usat=(mW-mfl)hW(p)+mflhs(p) (14)
式中:
Usat--系统中工质的总能量;
mW--系统中液相水的质量;
mfl--闪蒸质量;
hw--饱和水的比焓值,是压力p的函数;
hs--饱和水蒸气的比焓值,是压力p的函数;
转化后闪蒸质量表示为:
为保证计算结果的科学性与可靠性,采用上述两种计算结果的最小值作为最终的闪蒸质量;
当系统中蒸汽某时刻下对应的比内能小于其饱和比内能时,蒸汽将发生析出,根据能量平衡规律可以表示为:
msus=mrainoutuw(ps)+(ms-mrainout)us(ps) (16)
式中:
ms--过饱和蒸汽的质量;
us--蒸汽的比内能;
mrainout--析出的水的质量;
uw(ps)--蒸汽分压下对应的饱和液相比内能;
us(ps)--蒸汽分压下对应的蒸汽比内能;
转化后析出质量最终表示为:
5.如权利要求1所述的严重事故下安全壳内控制体热力学响应综合分析方法,其特征在于,步骤8中计算气液两相与相界面之间的传热传质时,首先考虑气相及液相到相界面的对流传热量,当相界面与气液两相存在温差时,对流传热量可根据牛顿冷却公式进行计算:
Qg,i=htgA(Tg-Ti) (18)
Ql,i=htlA(Tw-Ti) (19)
式中:
Qg,i--气相对相界面的对流传热量;
htg--气相对流传热系数;
A--相界面面积;
Tg--气体温度;
Ti--相界面温度;
Ql,i--液相对相界面的对流传热量;
htl--液相对流传热系数;
Tw--液相温度;
相界面的能量平衡方程写为:
Qg,i+Ql,i=ΓΔh (20)
式中:
Qg,i--气相对相界面的对流传热量;
Ql,i--液相对相界面的对流传热量;
Γ--相变传质率;
Δh--相变潜热;
此外假定相变传质率是由气相和相界面之间的蒸汽质量扩散率决定的,因而根据热质比拟理论:
式中:
Γ--相变传质率;
--单位表面上质量变化率;
A--相界面面积;
htg--气相对流传热系数;
λ--气相热导率;
Ds--蒸汽的扩散系数;
Ms--蒸汽摩尔质量;
p--气体总压;
R--摩尔气体常数;
Tg--气体温度;
ps,i--蒸汽在相界面处的分压;
pstm,b--蒸汽在气相主流处的分压;
pNCG,avg--扩散边界层中不可凝气体平均分压;
Sc--施密特数;
Pr--普朗特数;
在确定相界面温度后,即可得到气液相之间的传质率,最终获得气液两相之间的传热总量。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202011351583.4A CN112613158B (zh) | 2020-11-26 | 2020-11-26 | 一种严重事故下安全壳内控制体热力学响应综合分析方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202011351583.4A CN112613158B (zh) | 2020-11-26 | 2020-11-26 | 一种严重事故下安全壳内控制体热力学响应综合分析方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN112613158A CN112613158A (zh) | 2021-04-06 |
CN112613158B true CN112613158B (zh) | 2024-02-23 |
Family
ID=75225910
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202011351583.4A Active CN112613158B (zh) | 2020-11-26 | 2020-11-26 | 一种严重事故下安全壳内控制体热力学响应综合分析方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN112613158B (zh) |
Families Citing this family (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN113821998B (zh) * | 2021-08-31 | 2024-03-29 | 中国船舶重工集团公司第七0三研究所 | 利用牛顿迭代法进行凝汽器实时动态仿真模型壳侧压力求解的方法 |
Citations (6)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US5619433A (en) * | 1991-09-17 | 1997-04-08 | General Physics International Engineering Simulation Inc. | Real-time analysis of power plant thermohydraulic phenomena |
CN108536952A (zh) * | 2018-04-03 | 2018-09-14 | 东北大学 | 一种确定铁水包内气液两相流气含率的计算方法 |
CN110580375A (zh) * | 2019-07-29 | 2019-12-17 | 中广核工程有限公司 | 一种基于两相流模型的核电站安全壳仿真方法及系统 |
CN111651851A (zh) * | 2019-03-04 | 2020-09-11 | 国家电投集团科学技术研究院有限公司 | 安全壳求解方法及安全壳求解器 |
CN111680458A (zh) * | 2020-06-03 | 2020-09-18 | 西安交通大学 | 一种适用于钠水直流蒸汽发生器的热工水力瞬态计算方法 |
CN111723450A (zh) * | 2019-03-04 | 2020-09-29 | 国家电投集团科学技术研究院有限公司 | 核电厂安全分析方法及系统 |
-
2020
- 2020-11-26 CN CN202011351583.4A patent/CN112613158B/zh active Active
Patent Citations (6)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US5619433A (en) * | 1991-09-17 | 1997-04-08 | General Physics International Engineering Simulation Inc. | Real-time analysis of power plant thermohydraulic phenomena |
CN108536952A (zh) * | 2018-04-03 | 2018-09-14 | 东北大学 | 一种确定铁水包内气液两相流气含率的计算方法 |
CN111651851A (zh) * | 2019-03-04 | 2020-09-11 | 国家电投集团科学技术研究院有限公司 | 安全壳求解方法及安全壳求解器 |
CN111723450A (zh) * | 2019-03-04 | 2020-09-29 | 国家电投集团科学技术研究院有限公司 | 核电厂安全分析方法及系统 |
CN110580375A (zh) * | 2019-07-29 | 2019-12-17 | 中广核工程有限公司 | 一种基于两相流模型的核电站安全壳仿真方法及系统 |
CN111680458A (zh) * | 2020-06-03 | 2020-09-18 | 西安交通大学 | 一种适用于钠水直流蒸汽发生器的热工水力瞬态计算方法 |
Non-Patent Citations (5)
Title |
---|
AP1000典型事故工况瞬态热工水力特性研究;王伟伟;中国博士学位论文全文数据库 (工程科技Ⅱ辑)(第6期);全文 * |
华龙一号非能动安全壳冷却系统热工水力分析;丘锦萌等;原子能科学技术;第54卷(第1期);第72-80页 * |
基于MELCOR程序的AP1000核电厂安全壳瞬态事故分析;肖红;曹志伟;冯英杰;杨志义;朱建敏;;清华大学学报(自然科学版)(11);全文 * |
蒸汽排放系统蒸汽冷凝器动态特性仿真研究;孔夏明;王苇;孟海波;刘现星;陈保同;;原子能科学技术(12);全文 * |
非能动安全壳外部下降段的热工水力分析;李乐等;核动力工程;第37卷(第2期);第43-47页 * |
Also Published As
Publication number | Publication date |
---|---|
CN112613158A (zh) | 2021-04-06 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
Thomas | Simulation of industrial processes for control engineers | |
CN109918787A (zh) | 基于有限体积法的输水管道内水气两相均质流的模拟方法 | |
CN104834773A (zh) | 一种直管式直流蒸汽发生器换热性能的仿真方法 | |
CN108595752A (zh) | 一种面向核动力装置仿真应用的单相水流体网络仿真方法 | |
CN115859851B (zh) | 一种液态金属耦合超临界二氧化碳共轭传热的计算方法 | |
CN114444413B (zh) | 一种板状燃料堆芯亚通道级三维热工水力分析方法 | |
CN112613158B (zh) | 一种严重事故下安全壳内控制体热力学响应综合分析方法 | |
CN111859677A (zh) | 实验室尺度天然气水合物分解有效渗透率模型选择方法 | |
Bieder et al. | Large Eddy Simulation of the injection of cold ECC water into the cold leg of a pressurized water reactor | |
Zhong et al. | Numerical study on the local heat transfer characteristic of supercritical CO2 in semicircular/circular channel under cooling condition | |
Gango | Numerical boron mixing studies for Loviisa nuclear power plant | |
Walter | Dynamic simulation of natural circulation steam generators with the use of finite-volume-algorithms–A comparison of four algorithms | |
CN113027549B (zh) | 超临界二氧化碳发电系统动态模型的建模方法 | |
Tentner et al. | Modeling of two-phase flow in a BWR fuel assembly with a highly-scalable CFD code | |
Zhou et al. | Modeling and analysis of hydrodynamic instabilities in two-phase flow using two-fluid model | |
CN116341399A (zh) | 一种基于物理约束神经网络的热工水力换热系数预测方法 | |
Kurki | Simulation of thermal hydraulic at supercritical pressures with APROS | |
JPH11305646A (ja) | 初期値生成方法および装置 | |
Zhao et al. | Experimental study and analysis on the interfacial drag of two-phase flow in porous media | |
Shin et al. | Comparative study of constitutive relations implemented in RELAP5 and TRACE–Part I: Methodology & wall friction | |
Utberg Jr | Nitrogen concentration sensitivity study of the lock exchange phenomenon in the high temperature test facility | |
Blinkov et al. | Investigation on the interphase drag and wall friction in vertically oriented upward and downward two-phase flows under accident conditions in light water reactors | |
CN112417681B (zh) | 一种蒸汽发生器一二次侧对流换热系数分布的估计方法 | |
CN112613240B (zh) | 一种用于严重事故下安全壳内流动分析的计算方法 | |
Morel et al. | Modeling of multisize bubbly flow and application to the simulation of boiling flows with the Neptune_CFD Code |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
PB01 | Publication | ||
PB01 | Publication | ||
SE01 | Entry into force of request for substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant |