快谱小堆堆芯三维瞬态核-热-力耦合计算方法研究

胡田亮,姜夺玉,李 达,王立鹏,曹 璐,张信一,苏春磊

(西北核技术研究所 强脉冲辐射环境模拟与效应国家重点实验室,陕西 西安 710024)

快谱小堆堆芯小、能谱硬、温度高、中子泄漏强,堆芯核-热-力之间存在着复杂的耦合关系。不同于热谱反应堆,温度载荷作用下的力学效应导致的几何变形是快谱小堆堆芯反应性反馈的重要来源[1],与传统堆芯中围绕工质(如水)的临界沸腾而延伸出来的传统堆芯安全分析方法不同,快谱小堆中的高力学行为与耦合效应也是快谱小堆安全分析中最为关键的内容之一[2]。因此必须研究基于核-热-力耦合的计算方法,为开展快谱小堆的瞬态安全分析奠定基础。

快中子脉冲反应堆在爆发脉冲时,可在瞬发超临界状态运行,裂变速率会在几毫秒内增加数个数量级,同时堆芯相对体积发生变化,具有热膨胀自熄灭机制。脉冲过程涉及了反应堆动力学、非稳态传热、力学分析及耦合过程。因此,开展快中子脉冲堆的瞬态分析,不仅可以验证中子学计算的精度,而且可对核-热-力耦合计算及力学变形反馈模型进行有效的验证[3]。本文开展快谱小堆堆芯三维瞬态核-热-力耦合计算方法初步研究,计算快中子脉冲堆Godiva的瞬发超临界瞬态过程。

1 多物理耦合计算模型

1.1 传热计算模型

对于传热计算,可用热传导方程表示:

(1)

式中:ρ为密度;cp为比定压热容;T为温度;k为导热系数;qf为体积释热率。

其中热源由中子学模型计算得到:

(2)

式中:G为中子能群数;φg为第g群中子通量密度;κΣf,g为第g群中子的能量产生截面。

1.2 力学计算模型

(3)

式中:u为位移;σ为应力;f为外力。

在本文中采用小变形假设,几何方程可表示为:

(4)

式中,ε为应变。

温度变化也能引起变形,由于温度变化引起的膨胀或收缩受到约束产生的应力就是热应力,当考虑温度变化时,根据Duhamel-Neumann原理,实际应力需在弹性应力中加入温度载荷产生的附加应力。若初始状态下,结构内部初始应力为零,此时结构的本构关系可表示为:

σij=Cijklεkl+βTΔTδij

(5)

式中:βT为热膨胀系数;ΔT为温度的变化。

1.3 中子动力学计算模型

在快中子脉冲堆Godiva中,由于采用高浓度的235U金属核燃料,中子能谱很硬,同时金属系统的温升相对较小,相关研究表明多普勒效应对反应性的影响可忽略[3-5]。主要的反馈效应是由于体积变形导致的中子泄漏引起的,因此必须考虑变形对中子学的影响。在本文中动力学计算采用小变形假设,则计算网格的体积应变可表示为:

(6)

式中:V为网格变形后的体积;V0为初始网格体积。

与其他类型的龙虾相比，淡水龙虾出现病害的情况比较少见。但是，如果稻田的水质较差或者饲养管理措施不合理，那么出现病虫害的可能性就会大大提高，为了最大限度地降低龙虾病害发生率，一方面要加强日常管理力度，另一方面要紧跟时代发展潮流，对龙虾病害问题进行预防，在施肥过程中尽量采用少量多次的方式，提高龙虾的抵抗力，保障龙虾健康。

假设在瞬态过程中,中子的能谱不变,则中子宏观反应截面的变化主要是由计算网格体积改变导致的原子核密度改变引起的,瞬态过程中中子截面更新计算方法如下式所示:

(7)

式中:Σ为网格变形后的中子宏观反应截面;Σ0为初始的中子宏观反应截面。

对于中子学计算,本文采用中子扩散时空动力学进行模拟:

(8)

式中:vg为第g群中子平均速度;φg为第g群中子通量密度;χp,g为瞬发中子到第g群中子的概率;χd,i,g为第i组缓发中子先驱核衰变后的缓发中子出现在第g群中的概率;λi为第i组缓发中子先驱核的衰变常量;βi为第i组缓发中子所占的份额;I为缓发中子先驱核总的组数。

1.4 中子动力学计算模型耦合计算模型

方程(1)、(3)和(8)及相应的初始条件、边界条件即构成了核-热-力耦合的方程组,联立求解该方程组即可得到中子通量、功率、温度、应力及位移等物理量的分布。在本文中空间离散采用有限元方法,对于时间离散,中子扩散时空动力学方程和瞬态传热方程均采用全隐式向后差分进行离散,动力方程采用直接积分法中的Newmark隐式算法[6]。经过空间和时间离散,则可将描述核-热-力的偏微分方程组转换为非线性方程组:

(9)

式中:φ为中子通量密度的节点参量;T为温度的节点参量;u为位移的节点参量;Kij和P为与中子通量、温度和位移相关的耦合系数。

方程(9)可直接联立求解,采用Picard或者Newton迭代法处理方程的非线性问题,但直接联立求解形成的矩阵规模很大,对计算机硬件要求高。本文首先采用算子分裂法,将中子场的求解分离出来,则可将核-热-力耦合问题转换为顺序求解热弹性动力学方程组(式(10))和中子动力学方程(式(11))。

(10)

(11)

采用JFNK方法(Jacobian-Free Newton Krylov method)求解热弹性动力学方程组(10),对于非线性方程组(10)将右端项移到左边,并写为算子形式:

F(X)=0 (X=[T,u])

(12)

采用牛顿法求解方程(式(12)),可将非线性问题转换为迭代求解如下的线性问题:

J(Xk)ΔXk+1=-F(Xk)

Xk+1=Xk+ΔXk+1

(13)

JFNK方法在外层采用Newton-Raphson方法处理偏微分方程组的非线性问题,内层采用克雷洛夫子空间法求解线性方程[7]。JFNK方法区别与传统Newton方法的是由于内层采用克雷洛夫子空间法,因此不需要显式的构造Jacobi矩阵,只需要Jacobi矩阵和向量的乘积,可采用差分近似方法计算得到,如式(14)所示。JFNK方法的求解效率与预处理方法相关,Krylov迭代中的预处理方法很多,包括应用于一般矩阵的代数预处理子到利用特殊问题类之特征的问题相关预处理子。本文中非线性方程的数值求解都是基于petsc程序包,计算中采用的预处理方法为代数多重网格方法[8-9]。

(14)

式中,ε为差分步长,计算方式如文献[9]所示。

瞬态过程中每一时间步的计算流程如图1所示。图1中:1) 根据上一步中子学计算得到的中子通量密度进行裂变热源计算;2) 采用JFNK方法求解热弹性动力学方程组;3) 根据热弹性动力学方程组求得的节点位移,计算每个网格的体积应变,根据式(7)更新计算中采用的宏观截面,同时将计算得到的节点位移传递给中子动力学计算的网格进行变形;4) 根据更新的宏观截面和网格进行中子动力学计算,得到通量分布。

2 基于MOOSE的多物理耦合程序开发与验证

本文中,三维核-热-力耦合求解程序是基于MOOSE多物理耦合框架开发的。MOOSE框架以面向对象的编程思想进行开发,对物理场偏微分方程中的各数学算子进行抽象封装,通过平台内不同的程序模块组合实现对某个物理场具体偏微分方程的描述。具备以下特征:1) 与维度无关的用户代码,一套代码对一、二、三维问题都适用;2) 基于连续或非连续Galerkin有限元的网格离散;3) 全隐式耦合;4) 允许采用非结构化网格,包含多种形状函数;5) 很强的网格自适应性;6) 内置并行功能,用户不用开发额外的并行代码;7) 内置的后处理等[10]。

作者基于MOOSE开发了基于中子扩散理论和偶阶输运理论的中子学计算模块,该中子学求解模块经过一系列基准题的验证[11]。传热和动力学方程求解是基于MOOSE中内置的HeatConduction模块和Tensormechanics模块,在本文中将这两个模块耦合用以求解热弹性问题。为验证传热和力学模块的正确性,本文选取劳伦斯利弗莫尔国家实验室(LANL)构建的金属铀-钼球壳堆[12]进行热应力分析计算,其中球壳的内径为5.08 cm,外径为7.62 cm。文献中给出的热加载关系为:

(15)

式中:Tmax为最大温升;b为脉冲半高宽;tpp为峰功率时刻。

在绝热边界条件下,计算脉冲半高宽为41 μs,堆芯平均温升为420 ℃,堆芯内外表面位移计算结果如图2所示。本文计算得到的结果与文献中给出的解析解吻合良好,说明程序可准确模拟脉冲堆热-力耦合过程。

中子学求解和热弹性问题的耦合求解是通过MOOSE的MultiApp和Transfer功能实现的。在计算中热力学求解程序作为MaterApp,中子学求解程序作为MultiApp,通过Transfer功能实现功率和位移在MaterApp和MultiApp之间的传递。

3 快中子脉冲堆Godiva-Ⅰ三维核-热-力耦合瞬态计算研究

Godiva-Ⅰ是由LANL建造的世界上第1座快中子脉冲堆。Godiva-Ⅰ是LANL为验证热膨胀自熄灭裂变脉冲的机理,将临界装置Lady Godiva改建成快中子脉冲堆。它是一个无反射层的高富集铀金属球形快中子脉冲堆。该装置的堆芯是由235U富集度为93.5%、密度18.75 g/cm3的铸造铀金属构成,缓发临界质量约为52.04 kg,堆芯金属球直径约为17 cm[13]。本文中对Godiva-Ⅰ初始周期为16.2 μs的瞬态过程进行了模拟计算。

3.1 计算模型

在本文的计算中,采用三维几何进行建模,网格如图3所示。计算中采用的传热和力学参数列于表1。传热计算中,金属球体表面设置为绝热边界条件,温度场初始条件设置为均匀温度场,温度为293.6 K。动力学计算中,金属球体表面设置为自由膨胀边界条件,初始位移和速度都设置为0[14-15]。

表1 Godiva-Ⅰ基本参数Table 1 Basic parameter of Godiva-Ⅰ

图3 Godiva-Ⅰ网格划分模型Fig.3 Meshing model of Godiva-Ⅰ

多群中子截面和缓发中子先驱核参数是采用Serpent程序在293.6 K的温度下计算制作得到的[16]。中子学计算采用的能群结构列于表2[17]。中子学计算先求解稳态特征值问题,并将功率归一化为500 W,作为瞬态计算的输入条件。瞬态计算的持续时间为0.001 s,时间步长为1×10-7s。

表2 能群结构Table 2 Group structure

3.2 脉冲瞬态过程计算结果

使用编写的多物理耦合程序对Godiva-Ⅰ初始周期为16.2 μs的瞬发超临界瞬态过程进行了计算,计算结果如图4所示。从图4可看出,程序计算值和实验值吻合良好。

图4 瞬态过程中裂变率计算值和实验值的比较Fig.4 Comparison between numerical and experimental results of fission rate

瞬态过程中平均温升随时间的变化如图5所示。由图5可看出,引入正反应性后,在开始阶段,由于系统功率上升较慢,平均温度上升也较小,随后堆芯功率急剧上升,累积的热量导致金属铀的温度升高,热膨胀导致球体发生膨胀,产生位移;随后由于堆芯膨胀,金属密度减少,中子泄漏增强,引入负的反应性,堆芯处于次临界状态,导致堆芯功率迅速降低,堆芯平均温度变化逐渐减缓。图6为外表面位移和速度。由图6可发现,由于裂变功率的迅速升高和降低,堆芯震荡,导致堆芯相继处于压缩和膨胀的状态,产生的热惯性效应,导致堆芯表面位移、速度和加速度产生震荡,同时由于变形导致堆芯反应性的改变,堆芯功率在500 μs后也产生周期性的震荡,如图7所示。在瞬态计算结果中还观察到堆芯功率分布形状变化较小,因此传统的基于点堆动力学的中子学分析方法也是一种误差较小的近似方法。

图5 平均温升随时间的变化Fig.5 Time evolution of average temperature increase

图6 外表面位移和速度随时间的变化Fig.6 Time evolution of outer surface displacement and velocity

图7 裂变率随时间的变化Fig.7 Time evolution of fission rate

在图1中可看出,本文建立的多物理耦合模型,力学变形效应的反馈作用包括两个方面:1) 热膨胀导致堆芯体积增大,核子密度减小,引起宏观反应截面的变化,中子泄漏增强,这是一种负反馈效应;2) 堆芯体积增大,导致堆芯尺寸的变化,这是一种正的反馈作用。在传统的快堆分析中,一般认为第1种反馈是主导因素,第2种是二阶效应,在小变形条件下,一般忽略堆芯尺寸变化引起的反应性反馈效应。本文分析了Godiva-Ⅰ在瞬发超临界瞬态条件下,变形导致的密度变化和几何形状变化对瞬态计算结果的影响。图8中的两个算例都考虑了密度负反馈效应,区别是是否考虑堆芯尺寸变化的反馈效应。从图8可看出,采用动网格对几何进行变形,几何尺寸增加引入的是正的反应性,会使峰值功率增加,而只考虑密度变化反馈效应得到的计算值偏低。

总之,在Godiva瞬发超临界瞬态计算过程中,由于多普勒反馈效应是可忽略的,热载荷导致的力学变形效应是计算中唯一考虑的反应性反馈机制。从计算结果可看出,本文的方法可准确考虑由于热膨胀导致的反应性反馈效应。

4 结论

本文基于MOOSE多物理耦合框架,开展了三维核-热-力耦合计算研究,开发了多物理耦合计算程序,并应用该程序计算了快中子脉冲堆Godiva-Ⅰ的瞬发超临界瞬态,实现了瞬态过程中,裂变率、温度、位移等物理量的准确模拟,计算结果与实验值符合较好。计算结果表明,该计算模型可准确模拟热膨胀引起的反应性反馈效应,为进一步开展快谱小堆堆芯多物理耦合安全分析奠定了基础。