一种可裂变水溶液系统二维辐解气体空泡份额分布计算方法

技术领域

本发明属于核与辐射安全技术领域，具体涉及一种可裂变水溶液系统二维辐解气体空泡份额分布计算方法。

背景技术

均匀铀盐水溶液常见于均质水溶液堆和乏燃料后处理过程中。当可裂变水溶液系统临界时，裂变产物与溶液中的物质碰撞，发生辐射分解，形成辐射分解产物。相比与其他裂变产物，裂变碎片造成的辐射分解产物约占96％，且其与水分子碰撞形成的初级产物主要是氢气和过氧化氢。此外沿裂变碎片运动轨迹快速释放的能量会造成水蒸气的产生，然而在溶液未达到饱和温度的情况下，水蒸气会快速溶解。因此沿裂变碎片运动轨迹形成的气体大多为氢气。

在未达到氢气饱和浓度时，氢气迅速产生并聚集形成的小气泡会发生溶解，直至达到其饱和浓度，这些沿裂变碎片运动轨迹产生的小气泡不再溶解，溶液中溶解的辐射裂解气体扩散进入气泡，造成气泡快速生长。

温度升高和气泡的快速生长会造成惯性压力，对容器壁造成安全性威胁。同时，这些稳定存在的气泡受浮力作用逸出溶液，会造成功率震荡。因此在均匀铀盐水溶液的瞬态分析中，对辐解气泡行为的建模分析十分重要。

发明内容

本发明基于辐解气泡在溶液中的行为简化建模，提供了一种可裂变水溶液系统二维辐解气体空泡份额分布计算方法。

本发明采用如下技术方案来实现的：

一种可裂变水溶液系统二维辐解气体空泡份额分布计算方法，包括如下步骤：

(1)确定溶液系统几何结构参数、材料物性、网格划分参数、单位能量氢气产生量G，单位为：mol/J，搅浑系数D

(2)溶液中的辐解气泡分为三种：形核气泡、溶解气泡和生长气泡；形核气泡为沿裂变碎片运动轨迹产生的气泡，依据其在溶液中的是否生长分为两组气泡，其中形核气泡和溶解气泡为第一组气泡，生长气泡为第二组气泡，气泡的以下信息均需要利用矩阵存储，以体现变量二维空间分布特性：

a.两组气泡的半径，r

b.两组气泡的流速，v

c.两组气泡的数量密度，N

d.两组气泡的单位体积气体物质量，N

e.控制体内两组气泡的体积，V

f.两组气泡的气体密度，ρ

g.形核气泡的半径r

h.生长气泡的平衡浓度C

(3)采用上标0和last标记相关参数的初始值和上一时刻数值；从0时刻开始计算，液位为h

(4)根据P

(5)根据溶液压力P

(6)根据P

(7)根据T与T

(8)根据步骤(4)获得的α和步骤(5)获得的ρ

(9)通过重复步骤(4)到(8)，直至迭代前后的控制体内溶液压力相对误差最大值ε

(10)根据NL

(11)计算编号NL以下网格区域内，步骤(2)中c、d以外的物理量；根据气泡半径r

(12)根据上述步骤获得的参数α、P

(13)通过重复步骤(4)到(12)，重复构建并求解步骤(12)中的线性微分方程组，以实现方程组之间的耦合计算，直至NL＝NL′，且两次迭代的N

(14)将当前时刻计算值作为上一时刻的值，更新刻物理量：

(15)从零时刻开始，推进时间步，重复输入不同时刻的功率P、T，并按照步骤(4)到(14)不断进行计算，直至达到瞬态计算所需时刻。

本发明进一步的改进在于，所述步骤(1)中，溶液系统几何结构参数为环形溶液区域R

本发明进一步的改进在于，所述步骤(3)中，根据以下关系式，计算初始时刻溶液总质量：

所述步骤(4)中，当各组气泡数量密度N

其中，r

根据以下关系式计算编号NL以下的网格区域内，各组气泡体积和空泡份额：

本发明进一步的改进在于，所述步骤(5)中，根据以下方程计算网格编号NL

其中，方程组右侧最后一项代表溶液与气泡相互作用导致的动量耗散项，其与α成正比，其中Λ代表动量耗散系数；

动量方程的边界条件为：

u

u

其中，P

根据以下方程计算网格编号NL

根据以下关系式，计算编号NL

本发明进一步的改进在于，所述步骤(6)中，根据以下关系式计算网格编号NL

上式中，κ

本发明进一步的改进在于，所述步骤(7)中，基于气液混合物状态方程，根据以下关系式计算溶液压力变化：

根据以下关系式，计算编号NL

本发明进一步的改进在于，所述步骤(8)中，根据以下关系式计算第J层网格以下的溶液总质量m：

其中，下标i、j分别代表径向和轴向网格位置；M为轴向总网格数量；

根据以下关系式确定液位高度：

其中，ΔH

根据以下关系式重新给出空泡份额分布：

本发明进一步的改进在于，所述步骤(10)中，根据以下关系式重新给出各物理量数值的分布：

其中，

根据以下关系式给出NL与NL′之间的气泡半径：

其中，若NL′＜NL

本发明进一步的改进在于，所述步骤(11)中，气泡内压与辐解气体浓度满足Heryy定律，当气泡内压与外压平衡时，辐解气体浓度称为平衡浓度，根据气泡半径，通过以下关系式计算平衡浓度：

式中，C

基于理想气体状态方程，根据以下关系式计算气泡内气体密度：

其中，ρ

采用气泡极限上升流速来描述气泡流速，根据曳力和浮力的平衡，气泡极限上升流速满足以下关系式：

式中，v

其中，Re为雷诺数，Eo为Etovos数，μ

根据以下关系式计算生长气泡与溶液的相间界面传质系数：

其中，D为气体在溶液中的扩散速率(m

根据以下关系式计算相间界面面积密度：

根据以下关系式计算第二组气泡溶解时间：

本发明进一步的改进在于，所述步骤(12)中，气体物质的量、气泡数量密度和辐解气体浓度守恒方程依次为：

上式中，ζ为单个形核气泡内的气体物质量，单位为：mol/个；V

根据理想气体状态方程，计算单个形核气泡内的气体物质量：

辐解气体浓度在各边界面上的边界条件为：

由于边界面为液面处，辐解气泡仅通过输运逸出溶液，扩散项在该边界面的通量为0，因此通过以下关系式描述液面边界条件：

其他边界面上的边界条件则均满足以下关系式：

其中X为相关物理量，代表N

和现有技术相比较，本发明至少具备如下有益的技术效果：

1、本发明所阐述的适用于可裂变水溶液系统的建模分析，该方法从可裂变水溶液系统的实际物理过程出发，考虑了溶解的辐解气体扩散，气泡形核和气液传质造成的气泡生长过程，结合物理实际和合理简化，对气泡产生、生长、迁移的过程进行了完整建模。

2、本发明能够用于针对可裂变水溶液系统的模拟程序之中，尤其对于瞬发临界引入的反应性较大而存在惯性压力脉冲的情况，本发明能够计算惯性压力大小，并模拟获得辐解气泡的分布。

3、由于常用的辐解气体模型未考虑气泡形核和生长的过程，依赖的经验参数较多，导致可裂变水溶液系统的模拟程序无法同时准确模拟第一裂变峰值和展宽，而采用CFD的两相流模拟建模复杂，计算耗时耗力，效率较低。本发明无需划分大量的网格也可以精确的计算，计算量小，计算速度快，在保证一定精度的同时确保计算效率。

附图说明

图1为本发明提出的一种可裂变水溶液系统二维辐解气体空泡份额分布计算方法的流程示意图。

图2为本发明实施方式中容器内部环形空腔区域网格划分试图。

图3为本发明实施方式中耦合程序计算的功率、压力与实验测量的结果对比图。

图4为本发明实施方式中耦合程序计算的功率、压力与实验测量的结果对比图。

图5为本发明实施方式中0s、0.04s、0.05s和0.55s时各网格内的空泡份额分布图。

具体实施方式

下面结合附图对本发明作进一步详细说明。

一种装载了铀浓度为C

(1)确定溶液系统几何结构参数、材料物性、网格划分参数、单位能量氢气产生量G(mol/J)、搅浑系数D

(2)将辐解气泡依据形核气泡在溶液中的是否生长分为两组气泡，其中形核气泡和溶解气泡为第一组气泡，生长气泡为第二组气泡，气泡的以下信息均需要利用二维数组存储，以体现变量二维空间分布特性：

a.两组气泡的半径，r

b.两组气泡的流速，v

c.两组气泡的数量密度，N

d.两组气泡的单位体积气体物质量，N

e.控制体内两组气泡的体积，V

f.两组气泡的气体密度，ρ

g.形核气泡的半径r

h.生长气泡的平衡浓度C

(3)从0时刻开始计算，液位为h

初始时刻溶液总质量

(4)根据P

a)当各组气泡数量密度N

其中，r

b)根据以下关系式计算编号NL以下的网格区域内，各组气泡体积和空泡份额：

(5)根据溶液压力P

a)忽略液位变化造成的压力变化，网格编号NL

其中，方程组右侧最后一项代表溶液与气泡相互作用导致的动量耗散项，其与α成正比，其中Λ代表动量耗散系数。动量方程的边界条件为：

u

u

其中，P

b)根据以下方程可计算网格编号NL

根据以下关系式，可计算编号NL

(6)根据P

网格编号NL

上式中，κ

(7)根据T与T

气液混合物状态方程：

根据以下关系式，可计算编号NL

(8)根据步骤(4)获得的α和步骤(5)获得的ρ

a)第J层网格以下的溶液总质量m：

其中，下标i、j分别代表径向和轴向网格位置；M为轴向总网格数量。

b)液位高度：

其中，ΔH

c)空泡份额：

(9)通过重复步骤(4)到(8)，直至迭代前后的控制体内溶液压力相对误差最大值ε

(10)根据NL

a)重新给出各物理量数值的分布：

其中，

b)给出NL与NL′之间的气泡半径：

其中，若NL′＜NL

(11)计算编号NL以下网格区域内，步骤(2)中c、d以外的物理量；根据气泡半径r

a)平衡浓度：

式中，C

b)气泡内气体密度：

其中，ρ

c)气泡极限上升流速满足以下关系式：

式中，v

/>

其中，Re为雷诺数，Eo为Etovos数，μ

d)生长气泡与溶液的相间界面传质系数：

其中，D为气体在溶液中的扩散速率(m

e)相间界面面积密度：

f)第二组气泡溶解时间：

(12)根据上述步骤获得的参数Δ、P

a)气体物质的量、气泡数量密度和辐解气体浓度守恒方程依次为：

/>

上式中，ζ为单个形核气泡内的气体物质量(mol/个)；V

b)边界条件：

辐解气体浓度在各边界面上的边界条件为：

由于边界面为液面处，辐解气泡仅通过输运逸出溶液，扩散项在该边界面的通量为0，因此可通过以下关系式描述液面边界条件：

其他边界面上的边界条件则均满足以下关系式：

其中X为相关物理量，可代表N

(13)通过重复步骤(4)到(12)，重复构建并求解步骤(12)中的线性微分方程组，以实现方程组之间的耦合计算，直至NL＝NL′，且两次迭代的N

(14)将当前时刻计算值作为上一时刻的值，更新刻物理量：

(15)从零时刻开始，推进时间步，重复输入不同时刻的功率P、T，并按照步骤(4)到(14)不断进行计算，直至达到瞬态计算所需时刻。

实施例

以法国实验装置SILENE进行实例展示，该装置环形容器内半径为0.038m，外径半为0.18m。在该系列实验编号为S2-173的试验中，装载体积为41.09L，铀浓度为71g/L(92.7wt％)、酸浓度为200ml/m

β

将本发明的计算方法与点堆动力学和热工水力学模型进行耦合计算。耦合程序计算的功率、压力与实验测量的结果对比图见图3-图4，0s、0.04s、0.05s和0.55s时各网格内的空泡份额分布如图5所示。

虽然，上文中已经用一般性说明及具体实施方案对本发明作了详尽的描述，但在本发明基础上，可以对之作一些修改或改进，这对本领域技术人员而言是显而易见的。因此，在不偏离本发明精神的基础上所做的这些修改或改进，均属于本发明要求保护的范围。