一种熔盐堆精细物理热工耦合计算的方法与流程

本发明涉及熔盐堆堆芯设计和安全技术领域，具体涉及一种熔盐堆精细物理热工耦合计算的方法，是一种获取反应堆三维功率、温度及其它物理热工水力参数精细分布的方法。

背景技术：

在压水堆、沸水堆、重水堆等传统反应堆中，含裂变元素的固态组件作为燃料而轻水和重水等液态流体作为冷却剂，但是在熔盐堆中液态熔盐既是燃料又是冷却剂，且溶解在熔盐中的缓发中子先驱核会流出堆外在堆外发生衰变，从而使得熔盐堆表现出强烈的物理热工耦合特性，导致熔盐堆的堆芯计算与传统反应堆具有显著区别。

技术实现要素：

为了精确获得熔盐堆堆芯内的功率、温度及其它物理热工水力参数的精细分布，本发明提出了一种能精细且快速进行熔盐堆物理热工水力耦合计算的方法，弥补了现有技术在计算过程中几何适应性差、计算精度低、计算时间长而无法满足工程需要的缺点。

为了实现上述目的，本发明采取了以下技术方案予以实施：

一种熔盐堆精细物理热工耦合计算的方法，包括如下步骤：

步骤1：对熔盐堆堆芯进行网格划分，将堆芯分为熔盐区域和反射层区域，由熔盐和反射层的初始温度场，根据公式(1)计算熔盐区域和反射层区域的中子宏观反应截面包括中子扩散系数、宏观吸收截面、宏观裂变截面、宏观散射截面和中子能谱，由缓发中子先驱核浓度的初始分布，根据公式(2)计算熔盐区域的缓发中子源项，由公式(1)和公式(2)计算得到的结果，利用公式(3)采用有限体积法在熔盐区域和反射层区域进行包含缓发中子源项的中子扩散计算，求解得到中子通量密度的分布和有效增殖因子；

式中：

Σ(T)——中子宏观反应截面

T——温度

T_av——拟合状态点平均温度

T_std——拟合状态点方差

a₀a₁a₂a₃——拟合系数

式中：

——第g群缓发中子源项

χ_d,i,g——第i组缓发中子先驱核能谱

λ_i——第i组缓发中子先驱核衰变常数

——第i组缓发中子先驱核浓度

式中：

D_g——第g群中子扩散系数

——第g群中子通量密度

Σ_t,g——第g能群宏观中子移出截面

Σ_g′→g——第g'群到第g群的宏观中子散射截面

k_eff——有效增殖因子，即问题的特征值

β₀——总的缓发中子份额

χ_p,g——第g群中子裂变份额

(υΣ_f)_g′——第g′群宏观中子产生截面

——第g群缓发中子源项；

步骤2：由步骤1中计算得到的中子通量密度的分布，根据公式(8)计算得到堆芯中的功率分布，然后利用公式(9)在熔盐区域进行热工水力计算，得到熔盐区域温度场和速度场的分布；

式中：

S_T——功率密度

(κΣ)_fg——第g群中子能量产生截面

φ_g——第g群中子通量密度

式中：

ρ——熔盐密度

——熔盐速度

P——压力

——粘性应力张量

——重力加速度

T——熔盐温度

k——熔盐热导率

C_p——定压比热容

S_T——功率密度

步骤3：首先根据步骤1中计算得到的中子通量密度和有效增殖因子，利用公式(10)计算缓发中子先驱核源项分布，然后根据步骤2计算得到的熔盐区域的速度场分布，在熔盐区域利用公式(11)采用有限体积法进行缓发中子先驱核方程的求解，得到熔盐区域中的缓发中子先驱核浓度的分布；

式中：

S_d,i——第i组缓发中子先驱核源项

β_i——第i组缓发中子份额

k_eff——有效增殖因子

(υΣ_f)_g——第g群中子产生截面

——第g群中子通量密度

式中：

——熔盐区域速度场分布

——第i组缓发中子先驱核浓度分布

S_d,i——第i组缓发中子先驱核源项分布，由公式(6)计算得到

λ_i——第i组缓发中子衰变常数

步骤4：由步骤2计算得到的温度分布，根据公式(1)更新步骤1中的中子宏观截面，由步骤3计算得到的缓发中子先驱核浓度分布，根据公式(2)更新步骤1中的缓发中子源项，对步骤1、步骤2和步骤3进行迭代计算，直到堆芯有效增殖因子、功率分布和熔盐温度分布收敛，其中收敛条件为：有效增殖因子收敛限为10^-5,功率分布和熔盐温度分布收敛限为10^-4。

与现有技术相比，本发明有如下突出优点：

1.熔盐堆具有几何结构比较复杂的特点，传统的基于节块方法或者有限差分方法的堆芯计算方法不再适用，本发明中中子扩散计算基于非结构网格的有限体积法，对熔盐堆设计计算中采用的复杂几何有较强的适应性；

2.考虑了中子扩散和热工水力及缓发中子先驱核计算区域不同的特点，在熔盐区域和反射层区域进行中子计算，而热工水力计算和缓发中子先驱核计算只在熔盐区域进行，且中子学计算和热工水力计算重叠的区域采用相同的网格，避免了采用网格映射技术增加计算负担，节省内存和计算时间；

3.改进了熔盐堆物理热工水力耦合迭代计算流程，，缓发中子先驱核的求解在中子扩散计算和热工水力计算后进行，采用此流程有效减少了求解缓发中子先驱核对流问题的次数，显著减少了物理热工耦合计算的时间。

附图说明

图1熔盐堆物理热工耦合计算流程图。

图2熔盐堆中缓发中子先驱核浓度分布。

具体实施方式

下面结合附图和具体实施方式对本发明作进一步详细说明：

步骤1：堆芯中子通量计算

中子扩散计算中，中子宏观截面包括中子扩散系数、宏观吸收截面、宏观裂变截面、宏观散射截面和中子能谱是通过由组件均匀化计算得到的群常数拟合而成的多项式计算得到的，多项式的表达形式如公式(1)所示：

式中：

Σ(T)——中子宏观反应截面

T——温度

T_av——拟合状态点平均温度

T_std——拟合状态点方差

a₀a₁a₂a₃——拟合系数

熔盐堆中，熔盐既是燃料又是冷却剂，在中子裂变中产生的缓发中子先驱核会溶解在熔盐中，随熔盐流出堆外，在堆外发生衰变，所以熔盐堆中缓发中子先驱核的分布受熔盐流场的影响，这是熔盐堆区别与传统的固体燃料反应堆的一个显著特点，所以在熔盐堆稳态中子学计算时必须考虑缓发中子先驱核浓度分布的影响。在熔盐堆中，，缓发中子先驱核的源项采用公式(2)计算得到：

熔盐堆中多群中子扩散方程可以表示为：

式中：

D_g——第g群中子扩散系数

——第g群中子通量密度

Σ_t,g——第g能群宏观中子移出截面

Σ_g′→g——第g'群到第g群的宏观中子散射截面

k_eff——有效增殖因子，即问题的特征值

β₀——总的缓发中子份额

χ_p,g——第g群中子裂变份额

(υΣ_f)_g′——第g′群中子宏观中子产生截面

——缓发中子源项

在熔盐堆中没有燃料组件且燃料以液体形态存在，故堆芯中燃料区域的形状由堆芯容器的形状决定，这是显著区别于传统的固体燃料反应堆的一个特点，固体燃料反应堆中燃料组件一般是规则的矩形或者六边形，几何建模一般比较简单，而在熔盐堆中，由于液体燃料的存在，几何建模一般比较复杂，故采用基于非结构网格的有限体积法进行中子扩散方程的空间离散。在任意的控制体，第g群中子平衡方程可以写为：

式中：

D_g——第g群中子扩散系数

——第g群中子通量密度

Σ_t,g——第g能群宏观中子移出截面

——中子源项包括散射源项、裂变源项和缓发中子源项

在任意控制体上做积分，可以得到：

式中：

D_g——第g群中子扩散系数

——第g群中子通量密度

Σ_t,g——第g能群宏观中子移出截面

——中子源项包括散射源项、裂变源项和缓发中子源项

应用高斯定理，将体积分变为面积分可以得到：

式中：

D_g——第g群中子扩散系数

——第g群中子通量密度

Σ_t,g——第g能群宏观中子移出截面

——中子源项包括散射源项、裂变源项和缓发中子源项

最终可以得到离散的方程：

式中：

D_f,g——第g群中子控制体表面中子扩散系数

——第g群中子通量密度

Σ_t,g——第g能群宏观中子移出截面

——中子源项包括散射源项、裂变源项和缓发中子源项

——控制体表面单位面积矢量

V——控制体体积

为保证通量连续和流连续，对于控制体表面的扩散系数，采用调合平均的方式进行计算。

步骤2：堆芯热工水力计算

与传统固态燃料核反应堆相比而言，液态燃料熔盐堆燃料熔盐即作为燃料，又是堆芯冷却剂，熔盐堆燃料熔盐在堆芯的流动换热是有内热源的传热换热，熔盐堆内热源可以根据步骤1计算得到的中子通量密度，采用公式(8)进行计算：

式中：

S_T——功率密度

(κΣ)_fg——第g群中子能量产生截面

φ_g——第g群中子通量密度

在熔盐堆中由于物理热工耦合强烈必须进行耦合计算，熔盐热工水力流动控制方程如公式(9)所示：

式中：

ρ——熔盐密度

——熔盐速度

P——压力

——粘性应力张量

——重力加速度

T——熔盐温度

k——熔盐热导率

C_p——定压比热容

S_T——功率密度

采用基于压力-预测修正的SIMPLE算法求解热工水力方程，通过不断修正计算结果，反复迭代，最后求解出压力和速度的收敛解，SIMPLE算法的基本思想为：基于假定的或者上次迭代计算所得到的压力场，求解离散形式的动量方程，获得速度场。计算过程中所采用的压力场是假定的，如此计算所得的速度场一般是不能满足连续性方程，因此，必须对这种给定的压力场进行修正。为了速度场能够满足连续性方程，把由动量方程的离散形式获得的压力与速度的关系代入连续性方程的离散形式，从而得到压力修正方程，由此压力修正方程求解压力修正值。利用修正后的压力场，求得新的速度场。然后检查速度场是否收敛，若不收敛，用修正后的压力值作为给定的压力场，开始下一层次的计算。如此反复，直到获得收敛的解为止。

步骤3：缓发中子先驱核计算

首先根据步骤1中计算得到的中子通量密度和有效增殖因子，利用公式(10)计算缓发中子先驱核源项分布

式中：

S_d,i——第i组缓发中子先驱核源项

β_i——第i组缓发中子份额

k_eff——有效增殖因子

(υΣ_f)_g——第g群中子产生截面

——第g群中子通量密度

然后进行缓发中子先驱核浓度的计算，熔盐堆中考虑缓发中子先驱核的对流作用的缓发中子先驱核平衡方程为：

式中：

——熔盐区域速度场分布

——第i组缓发中子先驱核浓度分布

S_d,i——第i组缓发中子先驱核源项分布

λ_i——第i组缓发中子衰变常数

和步骤1中的中子扩散方程相同，公式(11)也采用有限体积法进行空间离散，对公式(11)在控制体上进行积分，得：

式中：

——熔盐区域速度场分布

——第i组缓发中子先驱核浓度分布

S_d,i——第i组缓发中子先驱核源项分布

λ_i——第i组缓发中子衰变常数

对式(12)应用高斯公式，最终可以离散为：

式中：

——熔盐区域速度场分布

C_i,f——第i组缓发中子控制体面上的缓发中子先驱核浓度

S_d,i——第i组缓发中子先驱核源项分布

λ_i——第i组缓发中子衰变常数

V——控制体体积

通过控制体体心的值可以计算得到面心的值及梯度，从而得到不同离散格式下的代数方程组，即可通过求解得到缓发中子先驱核的浓度值。

步骤4：耦合迭代计算流程

与传统的采用固体燃料的反应堆不同，熔盐堆的物理热工耦合分析有其独特的特点。熔盐堆的反馈包括(1)热工水力参数对物理计算中群常数的反馈；(2)熔盐流场对缓发中子先驱核浓度分布的反馈。

由于缓发中子的份额较少，稳态情况下缓发中子对熔盐堆的keff及中子通量分布及功率的影响很小，故缓发中子先驱核的计算堆芯物理计算和热工水力计算的影响都很小，传统方法中子通量和缓发中子先驱核是一起求解的，考虑到以上因素，本发明提出了在迭代过程中首先求解中子通量方程然后求解热工水力方程最后求解缓发中子先驱核对流扩散方程的迭代方案，避免了在求解缓发中子先驱核方程时，流场更新不及时的问题，减少了缓发中子先驱核的计算次数，有效节省了计算时间，最终的计算流程为：

由步骤2计算得到的温度分布，根据公式(1)更新步骤1中的中子宏观截面，由步骤3计算得到的缓发中子先驱核浓度分布，根据公式(2)更新步骤1中的缓发中子源项，对步骤1、步骤2和步骤3进行迭代计算，直到堆芯有效增殖因子、功率分布和熔盐温度分布收敛，其中收敛条件为：有效增殖因子收敛限为10^-5,功率分布和熔盐温度分布收敛限为10^-4。。

基于以上的理论模型，该发明的计算流程如图1所示，计算得到的缓发中子先驱核浓度分布如图2所示。