CN112613158A - 一种严重事故下安全壳内控制体热力学响应综合分析方法 - Google Patents
一种严重事故下安全壳内控制体热力学响应综合分析方法 Download PDFInfo
- Publication number
- CN112613158A CN112613158A CN202011351583.4A CN202011351583A CN112613158A CN 112613158 A CN112613158 A CN 112613158A CN 202011351583 A CN202011351583 A CN 202011351583A CN 112613158 A CN112613158 A CN 112613158A
- Authority
- CN
- China
- Prior art keywords
- mass
- phase
- gas
- energy
- calculation
- 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.)
- Granted
Links
- 238000004458 analytical method Methods 0.000 title claims abstract description 22
- 230000004044 response Effects 0.000 title claims abstract description 22
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 claims abstract description 91
- 238000004364 calculation method Methods 0.000 claims abstract description 68
- 238000012546 transfer Methods 0.000 claims abstract description 54
- 230000005501 phase interface Effects 0.000 claims abstract description 33
- 238000000034 method Methods 0.000 claims abstract description 32
- 230000008859 change Effects 0.000 claims abstract description 30
- 239000012530 fluid Substances 0.000 claims abstract description 27
- 239000007788 liquid Substances 0.000 claims abstract description 26
- 238000001704 evaporation Methods 0.000 claims abstract description 20
- 230000008020 evaporation Effects 0.000 claims abstract description 20
- 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 119
- 239000012071 phase Substances 0.000 claims description 86
- 239000007791 liquid phase Substances 0.000 claims description 37
- 238000011144 upstream manufacturing Methods 0.000 claims description 13
- 238000004134 energy conservation Methods 0.000 claims description 11
- 230000009471 action Effects 0.000 claims description 8
- 230000000704 physical effect Effects 0.000 claims description 8
- 238000006243 chemical reaction Methods 0.000 claims description 6
- 238000009792 diffusion process Methods 0.000 claims description 6
- 239000011159 matrix material Substances 0.000 claims description 6
- 239000004576 sand Substances 0.000 claims description 6
- 238000007906 compression Methods 0.000 claims description 5
- 238000001816 cooling Methods 0.000 claims description 3
- 239000011555 saturated liquid Substances 0.000 claims description 3
- 239000000126 substance Substances 0.000 claims 1
- 230000008569 process Effects 0.000 abstract description 15
- 230000001052 transient effect Effects 0.000 abstract description 2
- 230000007704 transition Effects 0.000 abstract description 2
- 239000000203 mixture Substances 0.000 description 5
- IJGRMHOSHXDMSA-UHFFFAOYSA-N Atomic nitrogen Chemical compound N#N IJGRMHOSHXDMSA-UHFFFAOYSA-N 0.000 description 4
- 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
- 238000010586 diagram Methods 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
- 238000001727 in vivo Methods 0.000 description 2
- 229910052757 nitrogen Inorganic materials 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000006835 compression Effects 0.000 description 1
- 230000002950 deficient 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
- JEGUKCSWCFPDGT-UHFFFAOYSA-N h2o hydrate Chemical compound O.O JEGUKCSWCFPDGT-UHFFFAOYSA-N 0.000 description 1
- 230000007246 mechanism Effects 0.000 description 1
- 238000012545 processing Methods 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 238000004088 simulation Methods 0.000 description 1
- 238000005728 strengthening Methods 0.000 description 1
- 238000009834 vaporization Methods 0.000 description 1
- 230000008016 vaporization Effects 0.000 description 1
Images
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
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的差分形式质量守恒方程：
式中：
fj--经流道j流入到控制体的体积流量，m3/s；
d--流道的上游；
Δt--时间步长，s；
假设控制体内流体为理想流体，能量方程中不再考虑由于粘性作用引起的机械能转换为热能的部分，微分、集总参数形式的能量守恒方程化为差分形式的方程，可得：
式中：
fj--经流道j流入到控制体的体积流量，m3/s；
d--流道的上游；
Δt--时间步长，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；
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；
Γ--质量源项的体密度，kg/m3·s
微分形式的质量守恒方程在控制体内积分后，将微分、集总参数形式的质量守恒方程化为差分形式的方程，考虑质量源项后可得控制体内工质i的差分形式质量守恒方程：
式中：
fj--经流道j流入到控制体的体积流量，m3/s；
Δt--时间步长，s；
上式右边第二项表示一个时间步长内通过所有与控制体相连的流道流入或流出控制体的工质i的质量，此项为正表示流体流入控制体，为负表示流体流出控制体；右边第三项表示控制体内非流动质量源项，此项包括控制体内部两相之间由于相变(闪蒸、析出与蒸发)引起的相间质量变化以及用户定义的外部源项；
控制体内微分形式的能量守恒方程为：
式中：
ρ--控制体内流体密度，kg/m3；
u--控制体内流体速度，m/s；
t--时间，s；
p--控制体压力，kPa
λ--流体的热导率，W/m·K；
T--控制体内流体温度，K；
Sh--流体的内热源，kg/m2·s2
上式右边第一项表示压强对流体微元体压缩所做的功，并假设控制体内流体为理想流体，能量方程中不再考虑由于粘性作用引起的机械能转换为热能的部分，将上述方程在控制体内积分后，将微分、集总参数形式的能量守恒方程化为差分形式的方程，可得：
式中：
fj--经流道j流入到控制体的体积流量，m3/s；
Δt--时间步长，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；
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 (6)
1.一种严重事故下安全壳内控制体热力学响应综合分析方法，其特征在于，包括如下步骤：
步骤1、获取控制体热工水力计算初始参数以及计算所需的时间控制参数、边界条件参数，完成初始化；
步骤2、根据计算时间控制参数判断是否结束计算运行，当不满足结束条件时，进入步骤3继续执行计算，满足结束条件时，输出结果信息并结束计算；
步骤3、更新时间步长，根据边界条件参数计算出各时间内边界源项与指定控制体之间进行的能量与质量交换总额；
步骤4、根据质量和能量守恒方程，计算控制体内部总质量和总能量；
步骤5、判断控制体内部流体类型并分工况求解控制体内部状态参数；
步骤6、计算不可凝气体状态参数及气相平均物性；
步骤7、计算饱和水的闪蒸以及过饱和蒸汽的析出；
步骤8、计算气液两相与相界面之间的传热传质；
步骤9、返回步骤2反复迭代计算。
2.如权利要求1所述的严重事故下安全壳内控制体热力学响应综合分析方法，其特征在于，步骤1中所述控制体热工水力的初始参数包括体积、标高、高度、水力直径、压力、气相初始温度、液相初始温度、各气体的体积份额、各组分的质量和气液两相的能量；
时间控制参数包括计算起始时间、终止时间、各时间段的最大、最小时间步长；
边界条件参数包括外部源项的开启时间与关闭时间、作用控制体的编号、源项的状态、种类、组分、温度、压力、流量。
3.如权利要求1所述的严重事故下安全壳内控制体热力学响应综合分析方法，其特征在于，步骤4中将微分、集总参数形式的质量守恒方程化为差分形式的方程，考虑质量源项后可得控制体内工质i的差分形式质量守恒方程：
式中：
fj--经流道j流入到控制体的体积流量；
d--流道的上游；
Δt--时间步长；
假设控制体内流体为理想流体，能量方程中不再考虑由于粘性作用引起的机械能转换为热能的部分，微分、集总参数形式的能量守恒方程化为差分形式的方程，可得：
式中：
fj--经流道j流入到控制体的体积流量；
d--流道的上游；
Δt--时间步长；
4.如权利要求1所述的严重事故下安全壳内控制体热力学响应综合分析方法，其特征在于，步骤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进行偏导数计算，则其导数矩阵表示为：
导数矩阵中第1行第1/3/4列元素分别为：
利用以上各式，根据牛顿迭代法，可得关于第k步的迭代方程如下：
上述方程中，除已知参数控制体总体积以及各工质的质量和气液两相的内能外，其它参数均通过当前状态下控制体内的温度和压力状态参数导出，不依赖于其它状态或外部参数，为全隐式形式的求解方程，通过对上述方程的求解，便可得到该条件下控制体内各相的温度和压力状态参数。
5.如权利要求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)--蒸汽分压下对应的蒸汽比内能；
转化后析出质量最终表示为：
6.如权利要求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 true CN112613158A (zh) | 2021-04-06 |
| CN112613158B 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) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN113821998A (zh) * | 2021-08-31 | 2021-12-21 | 中国船舶重工集团公司第七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 |
|---|
| 丘锦萌等: "华龙一号非能动安全壳冷却系统热工水力分析", 原子能科学技术, vol. 54, no. 1, pages 72 - 80 * |
| 孔夏明;王苇;孟海波;刘现星;陈保同;: "蒸汽排放系统蒸汽冷凝器动态特性仿真研究", 原子能科学技术, no. 12 * |
| 李乐等: "非能动安全壳外部下降段的热工水力分析", 核动力工程, vol. 37, no. 2, pages 43 - 47 * |
| 王伟伟: "AP1000典型事故工况瞬态热工水力特性研究", 中国博士学位论文全文数据库 (工程科技Ⅱ辑), no. 6 * |
| 肖红;曹志伟;冯英杰;杨志义;朱建敏;: "基于MELCOR程序的AP1000核电厂安全壳瞬态事故分析", 清华大学学报(自然科学版), no. 11 * |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN113821998A (zh) * | 2021-08-31 | 2021-12-21 | 中国船舶重工集团公司第七0三研究所 | 利用牛顿迭代法进行凝汽器实时动态仿真模型壳侧压力求解的方法 |
| CN113821998B (zh) * | 2021-08-31 | 2024-03-29 | 中国船舶重工集团公司第七0三研究所 | 利用牛顿迭代法进行凝汽器实时动态仿真模型壳侧压力求解的方法 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN112613158B (zh) | 2024-02-23 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN105869685B (zh) | 一种基于模拟核反应堆中子反应性反馈过程的热工水力实验装置的模拟核反应堆中子反应性反馈过程的控制方法 | |
| Ambrosini | Assessment of flow stability boundaries in a heated channel with different fluids at supercritical pressure | |
| Roth et al. | Theory and implementation of nuclear safety system codes–Part I: Conservation equations, flow regimes, numerics and significant assumptions | |
| CN106682376A (zh) | 参数随工况变化实际特性的全过程汽轮机建模及辨识方法 | |
| CN108763670A (zh) | 一种求解超临界二氧化碳反应堆布雷顿循环瞬态过程方法 | |
| CN112613158A (zh) | 一种严重事故下安全壳内控制体热力学响应综合分析方法 | |
| Bae et al. | Comparison of gas system analysis code GAMMA+ to S-CO2 compressor test data | |
| CN113027549B (zh) | 超临界二氧化碳发电系统动态模型的建模方法 | |
| Tuunanen et al. | Assessment of passive safety injection systems of ALWRs. Final report of the European Commission 4th framework programme. Project FI4I-CT95-004 (APSI) | |
| CN116341399A (zh) | 一种基于物理约束神经网络的热工水力换热系数预测方法 | |
| 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 | |
| Roth et al. | Comprehensive analyses of nuclear safety system codes | |
| CN113792449B (zh) | 核反应堆监测方法及系统 | |
| CN112417681B (zh) | 一种蒸汽发生器一二次侧对流换热系数分布的估计方法 | |
| Jiang et al. | Extension, Verification and Validation of molten salt in contact with non-condensable gas problem on RELAP/SCDAPSIM/MOD4. 0 code | |
| Shin et al. | Comparative study of constitutive relations implemented in RELAP5 and TRACE–Part I: Methodology & wall friction | |
| Saini | Experimental and numerical investigation of supercritical flow instability in two vertical heated parallel channels using CO2 | |
| CN112417781B (zh) | 核电蒸汽发生器出口饱和蒸汽质量流量估计方法及系统 | |
| Bury et al. | Validation of the heat and mass transfer models within a pressurized water reactor containment using the International Standard Problem No. 47 data | |
| Bostanci et al. | Pyflow Wellbore Simulator | |
| Kang et al. | General Heat Exchanger Modeling for Quasi-Steady Small Modular Reactor Heat Source Simulation | |
| Winters | 70 MPa fast-fill modeling & validation experiments | |
| Hoover | SMR Full-Power Scaling Analysis and Numerical Simulation for Nuclear Hybrid Energy System Testing | |
| Oha et al. | KAIST Supercritical CO2 pressurizing Experiment modeling with Simulink | |
| Tahvola | Modelling of PKL test facility with TRACE 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 |