Improvement of Radiolytic Gas Model for Criticality Accident Analysis of Nuclear Fuel Solution System
-
摘要: 对瞬态临界事故的准确模拟是核燃料溶液系统临界安全评估的关键因素。现有的辐解气体模型经验参数较多，导致功率特性预测存在较大偏差。为提高模拟精度和避免对模型中经验参数取值的依赖，需对辐解气体模型进行改进。基于对溶液中辐解气体行为的分析和简化假设，建立了包含辐解气体浓度、辐解气泡单位体积物质量和气泡数量密度的守恒模型，并将其与点堆中子动力学模型和二维导热模型相耦合，开发了溶液系统二维瞬态分析程序，通过日本TRACY实验进行了验证。结果表明，程序模拟值与实验数据符合较好，程序模型能够准确模拟溶液系统临界事故的功率变化。Abstract: The accurate simulation of transient criticality accident is a key factor in critical safety assessment of nuclear fuel solution system. However, the existing radiolysis gas models contain many empirical parameters, which result in significant deviations in the prediction of the power characteristics. To improve the simulation accuracy and avoid relying on the empirical parameters in the model, the radiolytic gas model needs to be improved. Based on the analysis of radiolysis gas behavior in solution and simplified assumptions, a conservation model including radiolytic gas concentration, mass per unit volume of radiolytic bubbles and number of density bubbles was established. This model was coupled with a point reactor kinetics model and a two-dimensional heat conduction model to develop a two-dimensional transient analysis code for solution systems. The code was verified with the Japanese TRACY experiment. The results show that the simulated values of the code are in good agreement with the experimental data, and the improved model can accurately simulate the power change of the solution system during critical accidents.
-
Key words:
- Nuclear fuel solution system /
- Critical safety /
- Radiolytic gas model
-
0. 引 言
在乏燃料后处理过程中，可裂变材料多以溶液的形式存在，对燃料溶液进行各种处理和操作时，可能会发生临界事故[1]。辐解气体是燃料溶液瞬发临界第一裂变脉冲后功率震荡的主要影响因素[2]，因此建立准确描述其行为的模型至关重要。现有的辐解气体模型主要分为两类：根据实验进行参数校准的模型（如FETCH[3]程序、TRACE[2]程序等）和采用简化假设来避免参数校准的模型（如TCCAHR[4]程序、Francisco J. Souto[5]等建立的模型）。前者预测能力受限于经验参数，而后者无法捕捉辐解气体的相关重要效应[6]。GeorgeE.Winter[6]等通过数值模拟开发了独立于系统特性的一维模型，然而该模型在模拟脉冲反应性引入瞬发临界工况时，二次功率峰值时刻过于延后。为减少对相关经验参数取值的依赖，并准确捕捉第一裂变脉冲峰后的功率变化特性，本研究在GeorgeE.Winter[6]模型的基础上，通过简化假设考虑了辐解气泡的形成，建立了包含溶解的辐解气体浓度、气体物质量和气泡数量的守恒模型，并将其与二维导热模型和点堆动力学模型相耦合，开发了溶液系统二维瞬态分析程序，并使用日本TRACY实验进行了验证以确保模型的准确性。
1. 理论模型
1.1 点堆中子动力学模型
点堆中子动力学方程为：
$$ \frac{{{\text{d}}P_{\mathrm{tol}}\left( t \right)}}{{{\text{d}}t}} = \frac{{R\left( t \right) - \beta }}{\varLambda }P_{\mathrm{tol}}\left( t \right) + \sum\limits_{i = 1}^6 {{\lambda _i}{c_i}\left( t \right)} $$ (1) $$ \frac{{{\text{d}}{c_i}\left( t \right)}}{{{\text{d}}t}} = \frac{{{\beta _i}}}{\varLambda }P_{\mathrm{tol}}\left( t \right) - {\lambda _i}{c_i}\left( t \right){\text{ }}\left( {i = 1,2, \cdots ,6} \right) $$ (2) 式中，$ P_{\mathrm{tol}}\left( t \right) $为t时刻系统总功率；$ R\left( t \right) $为堆内净反应性；$ \varLambda $为中子代时间；$ \beta $为缓发中子总份额；$ {\beta _i} $、$ {c_i}\left( t \right) $和$ {\lambda _i} $依次为第i组缓发中子份额、先驱核浓度和先驱核衰变常数。为了计算具有二维分布特点的空泡份额和温度，还需给出溶液区域内的功率密度分布函数$ f(r,{\textit{z}}) $：
$$ f(r,{\textit{z}}) = {J_0}\left( {2.405\frac{r}{{{R_{\text{F}}} + {\delta _r}}}} \right)\cos\left[\frac{\pi\left({\textit{z}}+\delta_{\textit{z}}\right)}{H_{\mathrm{F}}+2 \delta_{\textit{z}}}-\frac{\pi}{2}\right] $$ (3) 式中，δ为外推距离；下标r、z分别为径向和轴向坐标；J0为0阶贝塞尔函数；$ {R_{\text{F}}} $为溶液区域的半径；$ {H_{\text{F}}} $为溶液液位。反应性反馈计算时，各网格处的反应性反馈系数权重函数可取$ f(r,{\textit{z}}) $或$ {f^2}(r,{\textit{z}}) $，本文取为$ f(r,{\textit{z}}) $。考虑到空泡份额对功率分布的影响[6]，对二维功率密度分布函数进行了修正，获得如下功率密度P：
$$ P = P_{\mathrm{tol}}\left( t \right)\frac{{f\left( {r,{\textit{z}}} \right)\left( {1 - {\alpha _{\mathrm{g}}}} \right)}}{{\displaystyle \int\limits_V {f\left( {r,{\textit{z}}} \right)\left( {1 - {\alpha _{\mathrm{g}}}} \right)} {\mathrm{d}}V}} $$ (4) 式中，αg为空泡体积份额；V为溶液体积。
1.2 热工模型
溶液系统的结构通常为圆柱罐体，其几何结构如图1所示。
通过混合系数对温度分布的修正，溶液区域的自然对流可简化为二维导热，即：
$$ \alpha_{\text{l}}\frac{\partial T}{\partial t}=\left(a+D_{\text{mix}}\right)\nabla\cdot\left(\alpha_{\text{l}}\nabla T\right)+\frac{P}{\rho_1c_p} $$ (5) 式中，T为溶液温度；αl为溶液体积份额，其与αg满足关系式αl+αg=1；Dmix为混合系数；$ a $为溶液热扩散率；${\rho _{\text{l}}}$为溶液密度；cp为溶液比热容。溶液系统中气腔区域采用集总参数法处理，其他区域采用二维导热模型计算，与流体接触的固体壁面均考虑为自然对流换热边界条件，换热系数采用Churchill-Chu关系式[7]计算获得。此外本文还考虑了容器壁与外部环境的辐射换热。
1.3 辐解气体模型
在常见的核燃料溶液中，如硝酸铀酰溶液，水辐照分解的主要产物为H2和H2O2[8]，H2O2进一步分解形成H2O和O2，而硝酸则受辐照分解产生N2和O2[9]。因此，核燃料溶液中的辐解气体主要包含H2、O2和N2等。由于N2和O2的含量较少，本文辐解气体主要考虑H2。在过冷溶液中，水蒸气和H2主要沿着裂变碎片运动轨迹产生，并聚集成若干柱状或锥状气泡，水蒸气在10−8 s内凝结，形成一个半径约0.05 μm的微小H2气泡，在辐解气体浓度相对较低的溶液中，该H2气泡在10−6 s内溶解[10]。沿着裂变碎片运动轨迹形成的H2气泡半径可称为形核半径，其大小与燃料溶液浓度、温度、裂变产物粒子类型、溶液压力等因素有关，尺寸不一。
Henry定律给出了溶解于溶剂中气体摩尔浓度与其分压的关系：
$$ {p_{\text{g}}} = \frac{C}{H} $$ (6) 式中，pg为气体分压；H为Herry常数；C为气体浓度。当气泡内外压力平衡时，即：
$$ {p_{\text{g}}} = \left( {{p_{\text{l}}} + \frac{{2\sigma }}{{{R_{\text{b}}}}}} \right) = \frac{{{C_{{\text{eq}}}}}}{H} $$ (7) 式中，pl为溶液压力；σ为溶液的表面张力；Rb为气泡的半径；Ceq为气体平衡浓度。
在燃料溶液中，形核半径对应的平衡浓度又被称为临界浓度。辐解气体浓度达到临界浓度前，裂变碎片造成的微小气泡对功率的影响可忽略不计[11]。达到临界浓度后，部分小气泡直接生长或聚合形成大气泡，当大气泡流入浓度较低的局部区域时，还会发生气泡溶解。本研究假定达到临界浓度前，所有辐解气体均溶解于溶液中，则溶解气体的C、单位体积气体物质量M和气泡数量密度N满足以下守恒方程，即：
$$ {\alpha _{\text{l}}}\frac{{\partial C}}{{\partial t}} = GP - {S _{\text{M}}} - {Q_{\text{M}}} + {D_{{\text{mix}}}}\nabla \cdot \left( {{\alpha _{\text{l}}}\nabla C} \right) $$ (8) $$ \frac{\partial M}{\partial t}=S _{\text{M}}+Q_{\text{M}}-v_{\mathrm{b}}\frac{\partial M}{\partial {\textit{z}}}+D_{\text{mix}}\nabla^2M $$ (9) $$ \frac{{\partial N}}{{\partial t}} = \frac{{{S _{\text{M}}}}}{\xi } - {Q_{\text{N}}}- v_{\mathrm{b}}\frac{{\partial N}}{{\partial {\textit{z}}}} + {D_{{\text{mix}}}}{\nabla ^2}N $$ (10) 式中，G为辐解气体生成率；SM为形核气泡中生长气泡的气体物质量产生项；QM相间传质造成的单位体积气体物质量源项；QN为相间传质造成的气泡数量密度消失项； $ \xi $为单个形核气泡中的辐解气体物质量，其满足理想气体状态方程；vb为气泡流速。式（8）~式（10）右侧最后一项考虑了溶液自然对流的影响，其处理方式与式（5）类似。此外，式（9）、式（10）右侧第三项考虑了辐解气泡运动的影响，其中气泡流速可通过浮力与曳力达到平衡获得[12]，满足以下关系式：
$$ v_{\mathrm{b}}=\left[\frac{8R_{\text{b}}g\left(\rho_{\text{l}}-\rho_{\text{g}}\right)}{3\rho_{\text{l}}C_{\text{D}}}\right]^{0.5} $$ (11) 式中，${C_{\text{D}}}$为曳力系数；g为重力加速度；${\rho _{\text{g}}}$为气体密度。
大气泡的数量源项与裂变率成正比[13]，本研究假定在达到临界浓度后，作为大气泡形核位点的形核气泡产生量占总形核气泡产生量的比例为F1，则SM的表达式为：
$$ {S_{\text{M}}} = GP{F_1}\theta \left( {C - {C_{{\text{crit}}}}} \right) $$ (12) 式中，Ccrit溶解气体的临界浓度；θ(x)为Heviside函数，当x≤0时，θ=0，x>0时，θ=1。
由于气泡生长不会造成气泡数量变化，但气泡溶解会对气泡数量有影响。本文对气泡生长和溶解分别进行处理，式（9）和式（10）中的相间传质项如式（13）、式（14）所示。
$$ {Q_{\text{M}}} = {K_{\text{D}}}A\left( {C - {C_{{\text{eq}}}}} \right)\theta \left( {C - {C_{{\text{eq}}}}} \right) - \frac{M}{\tau }\theta \left( {{C_{{\text{eq}}}} - C} \right) $$ (13) $$ {Q_{\text{N}}} = \frac{N}{\tau }\theta \left( {{C_{{\text{eq}}}} - C} \right) $$ (14) 式中，KD为相间传质系数；A为相间界面面积浓度，满足A=4πNRb2；$ \tau $为气泡溶解所需时间。相间传质系数和气泡溶解时间可通过气体浓度扩散方程推导获得[14]：
$$ {K_{\text{D}}} = \kappa \left( {\frac{1}{{{R_{\text{b}}}}} + \frac{1}{{\sqrt {{\text{π}}\kappa {{{R_{\text{b}}}} \mathord{\left/ {\vphantom {{{R_{\text{b}}}} v}} \right. } v_{\mathrm{b}}}} }}} \right) $$ (15) $$ \tau = \frac{{{\rho _{\text{g}}}R_{\text{b}}^2}}{{2\kappa \left( {{C_{{\text{eq}}}} - C} \right){M_{\text{g}}}}} $$ (16) 式中，Mg为气体摩尔质量；$ \kappa $为辐解气体在溶液中的扩散速率。
气泡半径和单个形核半径中的气体物质量满足理想气体状态平衡方程，分别为：
$$ \left( {{p_{\text{l}}} + \frac{{2\sigma }}{{{R_{\text{b}}}}}} \right)N\left( {\frac{4}{3}{\text{π}}{R_{\text{b}}}^3} \right) = MRT $$ (17) $$ \left( {{p_{\text{l}}} + \frac{{2\sigma }}{{{r_0}}}} \right)\left( {\frac{4}{3}{\text{π}}{r_0}^3} \right) = \xi RT $$ (18) 式中，r0为形核半径；R为理想气体常数。基于式（17）迭代求解获得的气泡半径，可用于计算空泡体积份额：
$$ {\alpha _{\text{g}}} = \frac{4}{3}N{\text{π}}{R_{\text{b}}}^3 $$ (19) 微分方程式（8）~式（9）可采用有限容积法进行离散，并通过迭代求解获得各时刻的C、M、N和αg。
2. 计算结果分析与实验验证
基于上述的理论模型，本研究编制了二维柱坐标下核燃料溶液系统临界事故瞬态分析程序。由于辐解气体模型中未考虑溶液压力变化和沸腾等因素的影响，为避免这些因素对模型验证的干扰，本文选取日本TRACY实验装置的R100和R072工况进行模拟计算。
R100工况中，辐解气体浓度未达到临界浓度，溶液内不存在稳定的辐解气泡，系统净反应性的减小主要由温度的负反应性反馈所导致，该工况适用于点堆中子动力学和热工模块的验证。R072工况中，溶液内存在稳定的辐解气泡，功率波动受辐解气泡的影响较大，因此该工况适用于辐解气体模块的验证。
R100工况的计算值与实验值的对比见图2。尽管热工模型采用混合系数的方式来减小温度分布梯度，以减少扩散和自然对流两种换热方式导致的温度分布差异，但不可能完全消除这种偏差，且实验热电偶测量存在±1.5℃的误差，故温度计算值与实验值存在偏差是合理的。图2所示的功率和温度计算结果与实验值均符合较好，验证了点堆中子动力学模块和热工模块的正确性。
图 2 TRACY R100工况功率与温度模拟结果Figure 2. Power and Temperature Simulation Results for TRACY R100 ExperimentR072工况的计算结果如图3~图6所示。在图3所示阶段，辐解气体浓度未达到临界浓度，因此第一裂变功率峰值取决于系统的温度变化。
图 3 TRACY R072工况功率与反应性模拟结果Figure 3. Power and Reactivity Simulation Results for TRACY R072 Experiment
图 4 TRACY R072工况平均临界浓度、溶解氢气平均浓度、平均平衡浓度和气泡平均数量密度模拟结果Figure 4. Simulation Results of Average Critical Concentration, Average Dissolved Hydrogen Concentration, Average Equilibrium Concentration and Average Number Density of Hydrogen Bubbles for TRACY R072 Experiment
图 5 TRACY R072工况气泡平均半径和空泡份额模拟结果Figure 5. Simulation Results of Average Bubble Radius and Average Void Fraction for TRACY R072 Experiment
图 6 TRACY R072工况功率与温度模拟结果Figure 6. Power and Temperature Simulation Results for TRACY R072 Experiment在反应性引入后的约2 s时间内，当局部气体浓度达到临界浓度时，系统内会迅速形成大量可稳定存在的辐解气泡。如图4所示，气泡数量密度迅速增加，传质速率也相应增加。溶解的气体扩散进入大气泡后，如图5所示，大气泡快速生长至5×10−5 m，随后缓慢生长，最终尺寸约为1×10−4 m，平均空泡份额增加，平均平衡浓度相应减小（见图4）。气泡的生长会消耗溶解的辐解气体，当气体浓度减少至小于临界浓度时，稳定存在并生长的气泡数量密度源项为0，且受溶液对流和浮力的影响，气泡数量密度开始减少，传质速率也减小，平均空泡份额减小，导致功率在15 s左右开始回升（如图6所示）。
当气泡数量减少至溶液中几乎无大气泡时，气泡的平均半径减小为形核半径，同时平衡浓度由大气泡平衡浓度转变为临界浓度。此后，溶解的辐解气体重新累积，直至约70 s左右，局部浓度再次达到临界浓度，平均空泡份额再次增加，重复上述过程。由于气泡运动和溶液对流造成的浓度均匀化，尽管临界浓度逐渐降低，但局部浓度再次达到临界浓度时，溶解气体的平均浓度较前一次更高。
实验表明，二次功率峰后的功率下降速率加快，空泡份额再次增加的时刻大约在25 s左右，与计算结果偏差较大。其原因可能在于辐解产物中其他气体的浓度随着时间推移逐渐增加，导致H2临界浓度减小，使得本文所使用的临界浓度值偏高。R72工况的功率和温度计算结果如图6所示，温度计算值与实验值变化趋势一致，且功率的计算值与实验结果符合良好。虽然瞬态前期温度计算值与实验值有较大偏差，但此时空泡反馈占主导作用，温度偏差对功率影响不大。因此，在R72工况中，本研究开发的辐解气体模型能够再现气泡数量减少导致传质速率下降，空泡份额减少从而导致功率出现二次峰值的现象。
3. 结 论
本研究基于对辐解气体行为的分析和简化假设，建立了包括辐解气体浓度、单位体积气体物质量和气泡数量密度的辐解气体模型，通过与二维导热模型和点堆中子动力学模型的耦合，开发了用于溶液系统临界事故的二维瞬态分析程序，并基于TRACY实验装置的R100和R072工况对开发的模型进行了验证。结果表明，程序的模拟值与实验数据符合较好，模型能够准确再现辐解气泡在溶液中的瞬态行为。在后续研究中，还需考虑压力的影响，增加溶液注入及沸腾等相关模型，并与三维中子动力学模型相耦合，以尽可能展现溶液系统反应性引入后的功率变化趋势。
致谢：
本文研究内容受国防科工局乏燃料后处理科研专项中“后处理厂核临界安全事故研究” 项目资助，特此致谢。
-
图 2 TRACY R100工况功率与温度模拟结果
Figure 2. Power and Temperature Simulation Results for TRACY R100 Experiment
图 3 TRACY R072工况功率与反应性模拟结果
Figure 3. Power and Reactivity Simulation Results for TRACY R072 Experiment
图 4 TRACY R072工况平均临界浓度、溶解氢气平均浓度、平均平衡浓度和气泡平均数量密度模拟结果
Figure 4. Simulation Results of Average Critical Concentration, Average Dissolved Hydrogen Concentration, Average Equilibrium Concentration and Average Number Density of Hydrogen Bubbles for TRACY R072 Experiment
图 5 TRACY R072工况气泡平均半径和空泡份额模拟结果
Figure 5. Simulation Results of Average Bubble Radius and Average Void Fraction for TRACY R072 Experiment
-
[1] MCLAUGHLIN T P, MONAHAN S P, PRUVOST N L, et al. A review of criticality accidents 2000 revision: LA-13638[R]. Los Alamos, New Mexico, USA: Los Alamos National Lab., 2000. [2] BASOGLU B, YAMAMOTO T, OKUNO H, et al. Development of a new simulation code for evaluation of criticality transients involving fissile solution boiling: JAERI-Data/Code 98-011[R]. lbaraki-ken, Japan: Japan Atomic Energy Research Institute, 1998. [3] PAIN C C, DE OLIVEIRA C R E, GODDARD A J H, et al. Non-linear space-dependent kinetics for the criticality assessment of fissile solutions[J]. Progress in Nuclear Energy, 2001, 39(1): 53-114. doi: 10.1016/S0149-1970(01)00003-8 [4] 汪量子. 溶液堆的蒙特卡罗方法物理计算模型及特性研究[D]. 北京: 清华大学,2011. [5] SOUTO F J, KIMPLAND R H, HEGER A S. Analysis of the effects of radiolytic-gas bubbles on the operation of solution reactors for the production of medical isotopes[J]. Nuclear Science and Engineering, 2005, 150(3): 322-335. doi: 10.13182/NSE05-A2519 [6] WINTER G E, COOLING C M, EATON M D. A semi-empirical model of radiolytic gas bubble formation and evolution during criticality excursions in uranyl nitrate solutions for nuclear criticality safety assessment[J]. Annals of Nuclear Energy, 2022, 165: 108614. doi: 10.1016/j.anucene.2021.108614 [7] CHURCHILL S W, CHU H H S. Correlating equations for laminar and turbulent free convection from a vertical plate[J]. International Journal of Heat and Mass Transfer, 1975, 18(11): 1323-1329. doi: 10.1016/0017-9310(75)90243-4 [8] LANE J A, MACPHERSON H G, MASLAN F. Fluid fuel reactors[M]. Reading: Addison-Wesley, 1958: 102. [9] BIDWELL R M, KING L D P, WYKOFF W R. Radiolytic yields of nitrogen and hydrogen in water boilers[J]. Nuclear Science and Engineering, 1956, 1(6): 452-454. doi: 10.13182/NSE56-A18460 [10] SPIEGLER P, BUMPUS JR C F, NORMAN A. Production of void and pressure by fission track nucleation of radiolytic gas bubbles during power bursts in a solution reactor: NAA-SR-7086[R]. Canoga Park: Atomics International, 1962. [11] NAKAJIMA K, YAMANE Y, MIYOSHI Y. A kinetics code for criticality accident analysis of fissile solution systems: AGNES2: JAERI-Data/Code-2002-004[R]. Tokyo: Japan Atomic Energy Research Inst., 2002. [12] CELATA G P, D’ANNIBALE F, DI MARCO P, et al. Measurements of rising velocity of a small bubble in a stagnant fluid in one-and two-component systems[J]. Experimental Thermal and Fluid Science, 2007, 31(6): 609-623. doi: 10.1016/j.expthermflusci.2006.06.006 [13] KIMPLAND R H, KORNREICH D E. A two-dimensional multiregion computer model for predicting nuclear excursions in Aqueous Homogeneous Solution Assemblies[J]. Nuclear Science and Engineering, 1996, 122(2): 204-211. doi: 10.13182/NSE96-A24155 [14] EPSTEIN P S, PLESSet M S. Erratum: on the stability of gas bubbles in liquid‐gas solutions[J]. The Journal of Chemical Physics, 1951, 19(2): 256. -
下载:
下载:
计量
- 文章访问数: 535
- HTML全文浏览量: 58
- PDF下载量: 62
- 被引次数: 0
