CN110598304B - 一种空间核电推进系统球床反应堆物理热工耦合分析方法 - Google Patents

Abstract

一种空间核电推进系统球床反应堆物理热工耦合分析方法，主要步骤如下：1、划分通道和控制体，输入球床堆的结构参数和边界条件，设定计算时间和时间步长。2、堆芯初始化计算，得到零时刻的初始值。3、采用吉尔算法，计算当前t时刻的堆芯功率分布、流场、温度场和压力场。4、根据步骤3的计算结果，采用吉尔算法进行下一时刻的计算，直到设定的总计算时间，计算停止。5、若计算结果收敛，终止计算；否则读取上一次停止时刻的计算结果作为初值，重新启动计算，直至结果收敛。本发明可以分析球床堆堆芯的稳态热工水力特性及事故工况下的瞬态响应特性，为核电推进球床堆的设计及安全特性分析提供建议与指导。

Description

一种空间核电推进系统球床反应堆物理热工耦合分析方法

技术领域

本发明涉及空间核动力技术领域，具体涉及一种空间核电推进系统球床反应堆物理热工耦合分析方法。

背景技术

随着空间安全和探索任务的拓展，航天器将面临更为复杂、严峻的空间环境。传统的太阳能或蓄电池航天器功率范围小、环境依赖性强，已无法满足未来的任务需求。对于深空探测、载人火星飞行和一些对机动性要求高的军事任务等，环境适应性强、大功率长寿命、安全可靠的空间核动力系统成为最优甚至是唯一的选择。在20世纪60年代，美国和前苏联就成功在航天任务中使用了空间核电源，此后，美国、俄罗斯、法国、德国和日本等航天大国都积极开展了将核能应用于太空的探索和研究。为了满足我国未来军事需求与国家战略需要，研究空间核动力技术刻不容缓。

在空间核动力系统中，高温气冷球床反应堆具有功率密度高、结构紧凑、裂变产物完全包容等优点，可广泛应用于核热推进、核电推进及双模式推进，特别地，球床反应堆结合闭式布雷顿循环可满足空间兆瓦级的功率需求。针对用于核热推进的球床反应堆，新墨西哥大学的空间核动力研究所进行了一些堆芯物理计算和热工水力分析，而针对用于核电推进的球床反应堆，除提出相关概念设计外，并未进行深入研究。

发明内容

本发明公开一种空间核电推进系统球床反应堆物理热工耦合分析方法。该方法针对空间核动力系统中球床堆的结构特点，通过稳态热工水力计算，对堆芯的设计进行优化，对一些关键参数进行敏感性分析；通过瞬态计算，分析堆芯在一些典型事故工况下的响应特性，评估其安全性，为空间核电推进球床堆的设计及安全特性分析提供建议与指导。

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

一种空间核电推进系统球床反应堆物理热工耦合分析方法，步骤如下：

步骤1：将球床反应堆的堆芯根据需要沿轴向划分为若干个独立的通道，而后沿径向划分为若干控制体；输入球床反应堆的结构几何参数和边界条件，设定计算时间和时间步长；

步骤2：根据球床反应堆的结构几何参数和边界条件，进行堆芯的稳态初始化计算，得到堆芯的热功率以及每个控制体的冷却剂质量流量、温度和压力，作为零时刻的初始值；

堆芯的热功率：

P_t＝P_n (1)

第i个控制体的冷却剂质量流量：

M_i＝M_in (2)

第i个控制体的冷却剂压力：

P_i＝P_i+1+ΔP_fi+ΔP_ai (3)

第i个控制体的冷却剂温度：

第i个控制体的燃料球表面温度：

在堆芯的稳态初始化计算中，认为燃料球的各层材料温度等于燃料球的表面温度，即

T_z1,i＝T_gm,i＝T_z2,i＝T_hdc,i＝T_pyc,i＝T_u,i＝T_s,i (6)

式(1)-(6)中：

P_t——堆芯的热功率/W

P_n——堆芯的额定功率/W

M_i——第i个控制体的冷却剂质量流量/kg·s^-1

M_in——入口冷却剂质量流量/kg·s^-3

P_i——第i个控制体的冷却剂压力/Pa

P_i+1——第i+1个控制体的冷却剂压力/Pa

ΔP_fi——第i个控制体的摩擦压降/Pa

ΔP_ai——第i个控制体的加速压降/Pa

T_g,i——第i个控制体的冷却剂温度/T

T_g,i-1——第i-1个控制体的冷却剂温度/T

A_i——第i个控制体的流通面积/m²

A_i-1——第i-1个控制体的流通面积/m²

c_p,g——冷却剂比容/J·kg^-1·K^-1

ΔV_i——第i个控制体的体积/m³

q”'_c——堆芯的功率密度/W·m^-3

h_sg——对流换热系数/W·m^-2·K^-1

a_sg——对流换热面积与体积比，

T_s,i——第i个控制体的燃料球表面温度/T

T_z1,i——第i个控制体燃料球的ZrC涂层温度/T

T_gm,i——第i个控制体燃料球的活性区温度/T

T_z2,i——第i个控制体燃料球的TRISO颗粒ZrC涂层温度/T

T_hdc,i——第i个控制体燃料球的TRISO颗粒致密石墨层温度/T

T_pyc,i——第i个控制体燃料球的TRISO颗粒疏松热解碳层温度/T

T_u,i——第i个控制体燃料球的TRISO颗粒UC-ZrC核心温度/T；

步骤3：建立堆芯的中子物理模型、冷却剂热工水力模型、堆芯本体传热模型以及燃料球传热模型；

堆芯中子物理模型，采用点堆动力学方程求解堆芯的裂变功率，方程考虑六组缓发中子和反应性反馈；

式中：

P(t)——t时刻堆芯裂变功率/W

t——时间/s

ρ(t)——t时刻总反应性/$

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

Λ——中子代时间/s

λ_i——第i组缓发中子的衰变常数/s^-1

C_i(t)——t时刻第i组缓发中子先驱核的浓度/m^-3

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

冷却剂的热工水力模型：认为冷却剂在堆芯内的流动为不可压缩的径向一维流动，控制方程如下：

质量控制方程：

动量控制方程：

能量控制方程：

式中：

ρ_g——冷却剂密度/kg·m^-3

r——径向坐标/m

u——冷却剂径向速度/m·s^-1

P——冷却剂压力/Pa

S——动量源项

c_p,g——冷却剂比容/J·kg^-1·K^-1

T_g,——冷却剂温度/T

T_s,——燃料球表面温度/T

h_sg——对流换热系数/W·m^-2·K^-1

a_sg——对流换热面积与体积比，

d_p——燃料球直径/m

ε——球床孔隙率

球床堆的堆芯本体传热模型，不考虑球床的轴向导热，球床的能量方程如下：

式中：

ρ_s——燃料球密度/kg·m^-3

c_p,s——燃料球比容/J·kg^-1·K^-1

λ_eff——球床的有效热导率/W·m^-1·K^-1

q”'_c——堆芯的功率密度/W·m^-3

燃料球的传热模型，燃料球包括活性区和ZrC涂层，活性区由众多TRISO型燃料颗粒弥散在石墨基质中组成；

燃料球的ZrC涂层内没有内热源，其一维导热方程为：

式中：

ρ_z1——ZrC涂层密度/kg·m^-3

c_z1——ZrC涂层比容/J·kg^-1·K^-1

T_z1——ZrC涂层温度/T

λ_z1——ZrC涂层热导率/W·m^-1·K^-1

燃料球的活性区看作是一种TRISO型燃料颗粒和石墨基质的均匀混合物，导热方程为：

式中：

ρ_gm——活性区密度/kg·m^-3

c_gm——活性区比容/J·kg^-1·K^-1

T_gm——活性区温度/T

λ_gm——活性区的有效热导率/W·m^-1·K^-1

q”'_gm——活性区的功率密度/W·m^-3

TRISO型燃料颗粒以UC-ZrC燃料为核心，由内向外依次有三层包覆层：疏松热解碳层、致密石墨层和外ZrC涂层；

UC-ZrC核心的导热方程为：

式中：

ρ_u——UC-ZrC核心密度/kg·m^-3

c_u——UC-ZrC核心比容/J·kg^-1·K^-1

T_u——UC-ZrC核心温度/T

λ_u——UC-ZrC核心热导率/W·m^-1·K^-1

q”'_u——UC-ZrC核心的功率密度/W·m^-3

疏松热解碳层、致密石墨层和外ZrC涂层内没有内热源，具有相同的导热方程：

式中：

ρ——包覆层密度/kg·m^-3

c——包覆层比容/J·kg^-1·K^-1

T——包覆层温度/T

λ——包覆层热导率/W·m^-1·K^-1

采用控制容积积分法对控制方程(9)-(16)进行离散，消除空间微分项，仅保留时间微分项，离散后的控制方程为(17)-(24)：

M_i＝M_i-1 (17)

式中：

M_i-1——第i-1个控制体的冷却剂质量流量/kg·s^-1

M_i——第i个控制体的冷却剂质量流量/kg·s^-1

t——时间/s

P_i——第i个控制体的冷却剂压力/Pa

P_i+1——第i+1个控制体的冷却剂压力/Pa

ΔP_fi——第i个控制体的摩擦压降/Pa

ΔP_ai——第i个控制体的加速压降/Pa

Δr_i——第i个控制体的径向长度/m

A_i——第i个控制体的流通面积/m²

A_i-1——第i-1个控制体的流通面积/m²

c_g,i——第i个控制体的冷却剂比容/J·kg^-1·K^-1

c_g,i-1——第i-1个控制体的冷却剂比容/J·kg^-1·K^-1

ρ_g,i——第i个控制体的冷却剂密度/kg·m^-3

T_g,i——第i个控制体的冷却剂温度/T

T_g,i-1——第i-1个控制体的冷却剂温度/T

h_sg,i——第i个控制体的对流换热系数/W·m^-2·K^-1

a_sg,——对流换热面积与体积比，

T_s,i——第i个控制体的燃料球表面温度/T

T_s,i+1——第i+1个控制体的燃料球表面温度/T

T_s,i-1——第i-1个控制体的燃料球表面温度/T

c_s,i——第i个控制体的燃料球比容/J·kg^-1·K^-1

ρ_s,i——第i个控制体的燃料球密度/kg·m^-3

ΔV_i——第i个控制体的体积/m³

r_i——第i个控制体的径向坐标/m

r_i+1——第i+1个控制体的径向坐标/m

r_i-1——第i-1个控制体的径向坐标/m

L_i——第i个控制体的轴向高度/m

λ_eff,i——第i个控制体的球床有效热导率/W·m^-1·K^-1

λ_eff,i+1——第i+1个控制体的球床有效热导率/W·m^-1·K^-1

λ_eff,i-1——第i-1个控制体的球床有效热导率/W·m^-1·K^-1

q”'_c,i——第i个控制体的堆芯功率密度/W·m^-3

式中：

ρ_z1,i——ZrC涂层第i个控制体的密度/kg·m^-3

c_z1,i——ZrC涂层第i个控制体的比容/J·kg^-1·K^-1

ΔV_z1,i——ZrC涂层第i个控制体的体积/m³

T_z1,i——ZrC涂层第i个控制体的温度/T

T_z1,i+1——ZrC涂层第i+1个控制体的温度/T

T_z1,i-1——ZrC涂层第i-1个控制体的温度/T

r_z1,i——ZrC涂层第i个控制体的径向坐标/m

r_z1,i+1——ZrC涂层第i+1个控制体的径向坐标/m

r_z1,i-1——ZrC涂层第i-1个控制体的径向坐标/m

λ_z1——ZrC涂层热导率/W·m^-1·K^-1

式中：

ρ_gm,i——活性区第i个控制体的密度/kg·m^-3

c_gm,i——活性区第i个控制体的比容/J·kg^-1·K^-1

ΔV_gm,i——活性区第i个控制体的体积/m³

T_gm,i——活性区第i个控制体的温度/T

T_gm,i+1——活性区第i+1个控制体的温度/T

T_gm,i-1——活性区第i-1个控制体的温度/T

r_gm,i——活性区第i个控制体的径向坐标/m

r_gm,i+1——活性区第i+1个控制体的径向坐标/m

r_gm,i-1——活性区第i-1个控制体的径向坐标/m

λ_gm——活性区热导率/W·m^-1·K^-1

q”'_gm,i——活性区第i个控制体的功率密度/W·m^-3

式中：

ρ_u,i——UC-ZrC核心第i个控制体的密度/kg·m^-3

c_u,i——UC-ZrC核心第i个控制体的比容/J·kg^-1·K^-1

ΔV_u,i——UC-ZrC核心第i个控制体的体积/m³

T_u,i——UC-ZrC核心第i个控制体的温度/T

T_u,i+1——UC-ZrC核心第i+1个控制体的温度/T

T_u,i-1——UC-ZrC核心第i-1个控制体的温度/T

r_u,i——UC-ZrC核心第i个控制体的径向坐标/m

r_u,i+1——UC-ZrC核心第i+1个控制体的径向坐标/m

r_u,i-1——UC-ZrC核心第i-1个控制体的径向坐标/m

λ_u——UC-ZrC核心热导率/W·m^-1·K^-1

q”'_u,i——UC-ZrC核心第i个控制体的功率密度/W·m^-3

式中：

ρ_i——包覆层第i个控制体的密度/kg·m^-3

c_i——包覆层第i个控制体的比容/J·kg^-1·K^-1

ΔV_i——包覆层第i个控制体的体积/m³

T_i——包覆层第i个控制体的温度/T

T_i+1——包覆层第i+1个控制体的温度/T

T_i-1——包覆层第i-1个控制体的温度/T

r_i——包覆层第i个控制体的径向坐标/m

r_i+1——包覆层第i+1个控制体的径向坐标/m

r_i-1——包覆层第i-1个控制体的径向坐标/m

根据步骤2计算的初始值，采用吉尔算法对微分方程(7)、(8)、(17)-(24)进行求解，得到t时刻的堆芯功率、温度分布，冷却剂流场、温度场和压力场；

步骤4：根据步骤3的计算结果，利用吉尔算法进行下一时刻的计算，直到达到设定的总计算时间，计算停止；

步骤5：将输出的计算结果随时间的变化趋势绘制成曲线，若曲线已趋于平缓，则计算收敛，终止计算；否则读取上一次停止时刻的计算结果作为初值，再设定计算时间和时间步长，重新启动计算，直至结果收敛。

和现有技术相比较，本发明具备如下优点：

本发明采用的吉尔算法，克服了点堆中子物理方程中的刚性问题，实现了空间核电推进系统球床反应堆的中子物理与热工水力耦合计算。发明不仅考虑堆芯燃料球的表面温度，且采用等效热导率的方法计算燃料球中的最高温度，可以针对空间核电推进系统球床反应堆进行较为详细准确的设计优化和安全分析。本方法能计算通用的球床反应堆，如核热推进和双模式推进的球床反应堆，降低空间核动力的研究成本。

附图说明

图1为本发明的方法流程框图。

图2为空间核电推进系统球床堆的结构示意图。

图3为划分的堆芯通道和控制体。

具体实施方式

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

如图1所示，本发明一种空间核电推进系统球床反应堆物理热工耦合分析方法，步骤如下：

步骤1：空间核电推进系统球床堆的结构如图2所示，将球床反应堆的堆芯根据需要沿轴向划分为若干个独立的通道，而后沿径向划分为若干控制体，如图3所示；输入球床反应堆的结构几何参数和边界条件，设定计算时间和时间步长；

步骤2：根据球床反应堆的结构几何参数和边界条件，进行堆芯的稳态初始化计算，得到堆芯的热功率以及每个控制体的冷却剂质量流量、温度和压力，作为零时刻的初始值；

堆芯的热功率：

P_t＝P_n (1)

第i个控制体的冷却剂质量流量：

M_i＝M_in (2)

第i个控制体的冷却剂压力：

P_i＝P_i+1+ΔP_fi+ΔP_ai (3)

第i个控制体的冷却剂温度：

第i个控制体的燃料球表面温度：

在堆芯的稳态初始化计算中，认为燃料球的各层材料温度等于燃料球的表面温度，即

T_z1,i＝T_gm,i＝T_z2,i＝T_hdc,i＝T_pyc,i＝T_u,i＝T_s,i (6)

式(1)-(6)中：

P_t——堆芯的热功率/W

P_n——堆芯的额定功率/W

M_i——第i个控制体的冷却剂质量流量/kg·s^-1

M_in——入口冷却剂质量流量/kg·s^-3

P_i——第i个控制体的冷却剂压力/Pa

P_i+1——第i+1个控制体的冷却剂压力/Pa

ΔP_fi——第i个控制体的摩擦压降/Pa

ΔP_ai——第i个控制体的加速压降/Pa

T_g,i——第i个控制体的冷却剂温度/T

T_g,i-1——第i-1个控制体的冷却剂温度/T

A_i——第i个控制体的流通面积/m²

A_i-1——第i-1个控制体的流通面积/m²

c_p,g——冷却剂比容/J·kg^-1·K^-1

ΔV_i——第i个控制体的体积/m³

q”'_c——堆芯的功率密度/W·m^-3

h_sg——对流换热系数/W·m^-2·K^-1

a_sg——对流换热面积与体积比，

T_s,i——第i个控制体的燃料球表面温度/T

T_z1,i——第i个控制体燃料球的ZrC涂层温度/T

T_gm,i——第i个控制体燃料球的活性区温度/T

T_z2,i——第i个控制体燃料球的TRISO颗粒ZrC涂层温度/T

T_hdc,i——第i个控制体燃料球的TRISO颗粒致密石墨层温度/T

T_pyc,i——第i个控制体燃料球的TRISO颗粒疏松热解碳层温度/T

T_u,i——第i个控制体燃料球的TRISO颗粒UC-ZrC核心温度/T

步骤3：建立堆芯的中子物理模型、冷却剂热工水力模型、堆芯本体传热模型以及燃料球传热模型；

堆芯中子物理模型，采用点堆动力学方程求解堆芯的裂变功率，方程考虑六组缓发中子和反应性反馈；

式中：

P(t)——t时刻堆芯裂变功率/W

t——时间/s

ρ(t)——t时刻总反应性/$

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

Λ——中子代时间/s

λ_i——第i组缓发中子的衰变常数/s^-1

C_i(t)——t时刻第i组缓发中子先驱核的浓度/m^-3

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

冷却剂的热工水力模型：认为冷却剂在堆芯内的流动为不可压缩的径向一维流动，控制方程如下：

质量控制方程：

动量控制方程：

能量控制方程：

式中：

ρ_g——冷却剂密度/kg·m^-3

r——径向坐标/m

u——冷却剂径向速度/m·s^-1

P——冷却剂压力/Pa

S——动量源项

c_p,g——冷却剂比容/J·kg^-1·K^-1

T_g,——冷却剂温度/T

T_s,——燃料球表面温度/T

h_sg——对流换热系数/W·m^-2·K^-1

a_sg——对流换热面积与体积比，

d_p——燃料球直径/m

ε——球床孔隙率

球床堆的堆芯本体传热模型，不考虑球床的轴向导热，球床的能量方程如下：

式中：

ρ_s——燃料球密度/kg·m^-3

c_p,s——燃料球比容/J·kg^-1·K^-1

λ_eff——球床的有效热导率/W·m^-1·K^-1

q”'_c——堆芯的功率密度/W·m^-3

燃料球的传热模型，燃料球包括活性区和ZrC涂层，活性区由众多TRISO型燃料颗粒弥散在石墨基质中组成；

燃料球的ZrC涂层内没有内热源，其一维导热方程为：

式中：

ρ_z1——ZrC涂层密度/kg·m^-3

c_z1——ZrC涂层比容/J·kg^-1·K^-1

T_z1——ZrC涂层温度/T

λ_z1——ZrC涂层热导率/W·m^-1·K^-1

燃料球的活性区看作是一种TRISO型燃料颗粒和石墨基质的均匀混合物，导热方程为：

式中：

ρ_gm——活性区密度/kg·m^-3

c_gm——活性区比容/J·kg^-1·K^-1

T_gm——活性区温度/T

λ_gm——活性区的有效热导率/W·m^-1·K^-1

q”'_gm——活性区的功率密度/W·m^-3

TRISO型燃料颗粒以UC-ZrC燃料为核心，由内向外依次有三层包覆层：疏松热解碳层、致密石墨层和外ZrC涂层；

UC-ZrC核心的导热方程为：

式中：

ρ_u——UC-ZrC核心密度/kg·m^-3

c_u——UC-ZrC核心比容/J·kg^-1·K^-1

T_u——UC-ZrC核心温度/T

λ_u——UC-ZrC核心热导率/W·m^-1·K^-1

q”'_u——UC-ZrC核心的功率密度/W·m^-3

疏松热解碳层、致密石墨层和外ZrC涂层内没有内热源，具有相同的导热方程：

式中：

ρ——包覆层密度/kg·m^-3

c——包覆层比容/J·kg^-1·K^-1

T——包覆层温度/T

λ——包覆层热导率/W·m^-1·K^-1

采用控制容积积分法对控制方程(9)-(16)进行离散，消除空间微分项，仅保留时间微分项，离散后的控制方程为(17)-(24)：

M_i＝M_i-1 (17)

式中：

M_i-1——第i-1个控制体的冷却剂质量流量/kg·s^-1

M_i——第i个控制体的冷却剂质量流量/kg·s^-1

t——时间/s

P_i——第i个控制体的冷却剂压力/Pa

P_i+1——第i+1个控制体的冷却剂压力/Pa

ΔP_fi——第i个控制体的摩擦压降/Pa

ΔP_ai——第i个控制体的加速压降/Pa

Δr_i——第i个控制体的径向长度/m

A_i——第i个控制体的流通面积/m²

A_i-1——第i-1个控制体的流通面积/m²

c_g,i——第i个控制体的冷却剂比容/J·kg^-1·K^-1

c_g,i-1——第i-1个控制体的冷却剂比容/J·kg^-1·K^-1

ρ_g,i——第i个控制体的冷却剂密度/kg·m^-3

T_g,i——第i个控制体的冷却剂温度/T

T_g,i-1——第i-1个控制体的冷却剂温度/T

h_sg,i——第i个控制体的对流换热系数/W·m^-2·K^-1

a_sg,——对流换热面积与体积比，

T_s,i——第i个控制体的燃料球表面温度/T

T_s,i+1——第i+1个控制体的燃料球表面温度/T

T_s,i-1——第i-1个控制体的燃料球表面温度/T

c_s,i——第i个控制体的燃料球比容/J·kg^-1·K^-1

ρ_s,i——第i个控制体的燃料球密度/kg·m^-3

ΔV_i——第i个控制体的体积/m³

r_i——第i个控制体的径向坐标/m

r_i+1——第i+1个控制体的径向坐标/m

r_i-1——第i-1个控制体的径向坐标/m

L_i——第i个控制体的轴向高度/m

λ_eff,i——第i个控制体的球床有效热导率/W·m^-1·K^-1

λ_eff,i+1——第i+1个控制体的球床有效热导率/W·m^-1·K^-1

λ_eff,i-1——第i-1个控制体的球床有效热导率/W·m^-1·K^-1

q”'_c,i——第i个控制体的堆芯功率密度/W·m^-3

式中：

ρ_z1,i——ZrC涂层第i个控制体的密度/kg·m^-3

c_z1,i——ZrC涂层第i个控制体的比容/J·kg^-1·K^-1

ΔV_z1,i——ZrC涂层第i个控制体的体积/m³

T_z1,i——ZrC涂层第i个控制体的温度/T

T_z1,i+1——ZrC涂层第i+1个控制体的温度/T

T_z1,i-1——ZrC涂层第i-1个控制体的温度/T

r_z1,i——ZrC涂层第i个控制体的径向坐标/m

r_z1,i+1——ZrC涂层第i+1个控制体的径向坐标/m

r_z1,i-1——ZrC涂层第i-1个控制体的径向坐标/m

λ_z1——ZrC涂层热导率/W·m^-1·K^-1

式中：

ρ_gm,i——活性区第i个控制体的密度/kg·m^-3

c_gm,i——活性区第i个控制体的比容/J·kg^-1·K^-1

ΔV_gm,i——活性区第i个控制体的体积/m³

T_gm,i——活性区第i个控制体的温度/T

T_gm,i+1——活性区第i+1个控制体的温度/T

T_gm,i-1——活性区第i-1个控制体的温度/T

r_gm,i——活性区第i个控制体的径向坐标/m

r_gm,i+1——活性区第i+1个控制体的径向坐标/m

r_gm,i-1——活性区第i-1个控制体的径向坐标/m

λ_gm——活性区热导率/W·m^-1·K^-1

q”'_gm,i——活性区第i个控制体的功率密度/W·m^-3

式中：

ρ_u,i——UC-ZrC核心第i个控制体的密度/kg·m^-3

c_u,i——UC-ZrC核心第i个控制体的比容/J·kg^-1·K^-1

ΔV_u,i——UC-ZrC核心第i个控制体的体积/m³

T_u,i——UC-ZrC核心第i个控制体的温度/T

T_u,i+1——UC-ZrC核心第i+1个控制体的温度/T

T_u,i-1——UC-ZrC核心第i-1个控制体的温度/T

r_u,i——UC-ZrC核心第i个控制体的径向坐标/m

r_u,i+1——UC-ZrC核心第i+1个控制体的径向坐标/m

r_u,i-1——UC-ZrC核心第i-1个控制体的径向坐标/m

λ_u——UC-ZrC核心热导率/W·m^-1·K^-1

q”'_u,i——UC-ZrC核心第i个控制体的功率密度/W·m^-3

式中：

ρ_i——包覆层第i个控制体的密度/kg·m^-3

c_i——包覆层第i个控制体的比容/J·kg^-1·K^-1

ΔV_i——包覆层第i个控制体的体积/m³

T_i——包覆层第i个控制体的温度/T

T_i+1——包覆层第i+1个控制体的温度/T

T_i-1——包覆层第i-1个控制体的温度/T

r_i——包覆层第i个控制体的径向坐标/m

r_i+1——包覆层第i+1个控制体的径向坐标/m

r_i-1——包覆层第i-1个控制体的径向坐标/m

根据步骤2计算的初始值，采用吉尔算法对微分方程(7)、(8)、(17)-(24)进行求解，得到t时刻的堆芯功率、温度分布，冷却剂流场、温度场和压力场；

步骤4：根据步骤3的计算结果，利用吉尔算法进行下一时刻的计算，直到达到设定的总计算时间，计算停止；

步骤5：将输出的计算结果随时间的变化趋势绘制成曲线，若曲线已趋于平缓，则计算收敛，终止计算；否则读取上一次停止时刻的计算结果作为初值，再设定计算时间和时间步长，重新启动计算，直至结果收敛。

Claims (1)

1.一种空间核电推进系统球床反应堆物理热工耦合分析方法，其特征在于：步骤如下：

步骤1：将球床反应堆的堆芯根据需要沿轴向划分为若干个独立的通道，而后沿径向划分为若干控制体；输入球床反应堆的结构几何参数和边界条件，设定计算时间和时间步长；

步骤2：根据球床反应堆的结构几何参数和边界条件，进行堆芯的稳态初始化计算，得到堆芯的热功率以及每个控制体的冷却剂质量流量、温度和压力，作为零时刻的初始值；

堆芯的热功率：

P_t＝P_n (1)

第i个控制体的冷却剂质量流量：

M_i＝M_in (2)

第i个控制体的冷却剂压力：

P_i＝P_i+1+ΔP_fi+ΔP_ai (3)

第i个控制体的冷却剂温度：

第i个控制体的燃料球表面温度：

在堆芯的稳态初始化计算中，认为燃料球的各层材料温度等于燃料球的表面温度，即

T_z1,i＝T_gm,i＝T_z2,i＝T_hdc,i＝T_pyc,i＝T_u,i＝T_s,i (6)

式(1)-(6)中：

P_t——堆芯的热功率/W

P_n——堆芯的额定功率/W

M_i——第i个控制体的冷却剂质量流量/kg·s^-1

M_in——入口冷却剂质量流量/kg·s^-3

P_i——第i个控制体的冷却剂压力/Pa

P_i+1——第i+1个控制体的冷却剂压力/Pa

ΔP_fi——第i个控制体的摩擦压降/Pa

ΔP_ai——第i个控制体的加速压降/Pa

T_g,i——第i个控制体的冷却剂温度/T

T_g,i-1——第i-1个控制体的冷却剂温度/T

A_i——第i个控制体的流通面积/m²

A_i-1——第i-1个控制体的流通面积/m²

c_p,g——冷却剂比容/J·kg^-1·K^-1

ΔV_i——第i个控制体的体积/m³

q″′_c——堆芯的功率密度/W·m^-3

h_sg——对流换热系数/W·m^-2·K^-1

a_sg——对流换热面积与体积比，

T_s,i——第i个控制体的燃料球表面温度/T

T_z1,i——第i个控制体燃料球的ZrC涂层温度/T

T_gm,i——第i个控制体燃料球的活性区温度/T

T_z2,i——第i个控制体燃料球的TRISO颗粒ZrC涂层温度/T

T_hdc,i——第i个控制体燃料球的TRISO颗粒致密石墨层温度/T

T_pyc,i——第i个控制体燃料球的TRISO颗粒疏松热解碳层温度/T

T_u,i——第i个控制体燃料球的TRISO颗粒UC-ZrC核心温度/T

步骤3：建立堆芯的中子物理模型、冷却剂热工水力模型、堆芯本体传热模型以及燃料球传热模型；

堆芯中子物理模型，采用点堆动力学方程求解堆芯的裂变功率，方程考虑六组缓发中子和反应性反馈：

式中：

P(t)——t时刻堆芯裂变功率/W

t——时间/s

ρ(t)——t时刻总反应性/$

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

Λ——中子代时间/s

λ_i——第i组缓发中子的衰变常数/s^-1

C_i(t)——t时刻第i组缓发中子先驱核的浓度/m^-3

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

冷却剂的热工水力模型，认为冷却剂在堆芯内的流动为不可压缩的径向一维流动，控制方程如下：

质量控制方程：

动量控制方程：

能量控制方程：

式中：

ρ_g——冷却剂密度/kg·m^-3

r——径向坐标/m

u——冷却剂径向速度/m·s^-1

P——冷却剂压力/Pa

S——动量源项

c_p,g——冷却剂比容/J·kg^-1·K^-1

T_g,——冷却剂温度/T

T_s,——燃料球表面温度/T

h_sg——对流换热系数/W·m^-2·K^-1

a_sg——对流换热面积与体积比，

d_p——燃料球直径/m

ε——球床孔隙率

球床堆的堆芯本体传热模型，不考虑球床的轴向导热，球床的能量方程如下：

式中：

ρ_s——燃料球密度/kg·m^-3

c_p,s——燃料球比容/J·kg^-1·K^-1

λ_eff——球床的有效热导率/W·m^-1·K^-1

q″′_c——堆芯的功率密度/W·m^-3

燃料球的传热模型，燃料球包括活性区和ZrC涂层，活性区由众多TRISO型燃料颗粒弥散在石墨基质中组成；

燃料球的ZrC涂层内没有内热源，其一维导热方程为：

式中：

ρ_z1——ZrC涂层密度/kg·m^-3

c_z1——ZrC涂层比容/J·kg^-1·K^-1

T_z1——ZrC涂层温度/T

λ_z1——ZrC涂层热导率/W·m^-1·K^-1

燃料球的活性区看作是一种TRISO型燃料颗粒和石墨基质的均匀混合物，导热方程为：

式中：

ρ_gm——活性区密度/kg·m^-3

c_gm——活性区比容/J·kg^-1·K^-1

T_gm——活性区温度/T

λ_gm——活性区的有效热导率/W·m^-1·K^-1

q″′_gm——活性区的功率密度/W·m^-3

TRISO型燃料颗粒以UC-ZrC燃料为核心，由内向外依次有三层包覆层：疏松热解碳层、致密石墨层和外ZrC涂层；

UC-ZrC核心的导热方程为：

式中：

ρ_u——UC-ZrC核心密度/kg·m^-3

c_u——UC-ZrC核心比容/J·kg^-1·K^-1

T_u——UC-ZrC核心温度/T

λ_u——UC-ZrC核心热导率/W·m^-1·K^-1

q″′_u——UC-ZrC核心的功率密度/W·m^-3

疏松热解碳层、致密石墨层和外ZrC涂层内没有内热源，具有相同的导热方程：

式中：

ρ——包覆层密度/kg·m^-3

c——包覆层比容/J·kg^-1·K^-1

T——包覆层温度/T

λ——包覆层热导率/W·m^-1·K^-1

采用控制容积积分法对控制方程(9)-(16)进行离散，消除空间微分项，仅保留时间微分项，离散后的控制方程为(17)-(24)：

M_i＝M_i-1 (17)

式中：

M_i-1——第i-1个控制体的冷却剂质量流量/kg·s^-1

M_i——第i个控制体的冷却剂质量流量/kg·s^-1

t——时间/s

P_i——第i个控制体的冷却剂压力/Pa

P_i+1——第i+1个控制体的冷却剂压力/Pa

ΔP_fi——第i个控制体的摩擦压降/Pa

ΔP_ai——第i个控制体的加速压降/Pa

Δr_i——第i个控制体的径向长度/m

A_i——第i个控制体的流通面积/m²

A_i-1——第i-1个控制体的流通面积/m²

c_g,i——第i个控制体的冷却剂比容/J·kg^-1·K^-1

c_g,i-1——第i-1个控制体的冷却剂比容/J·kg^-1·K^-1

ρ_g,i——第i个控制体的冷却剂密度/kg·m^-3

T_g,i——第i个控制体的冷却剂温度/T

T_g,i-1——第i-1个控制体的冷却剂温度/T

h_sg,i——第i个控制体的对流换热系数/W·m^-2·K^-1

a_sg,——对流换热面积与体积比，

T_s,i——第i个控制体的燃料球表面温度/T

T_s,i+1——第i+1个控制体的燃料球表面温度/T

T_s,i-1——第i-1个控制体的燃料球表面温度/T

c_s,i——第i个控制体的燃料球比容/J·kg^-1·K^-1

ρ_s,i——第i个控制体的燃料球密度/kg·m^-3

ΔV_i——第i个控制体的体积/m³

r_i——第i个控制体的径向坐标/m

r_i+1——第i+1个控制体的径向坐标/m

r_i-1——第i-1个控制体的径向坐标/m

L_i——第i个控制体的轴向高度/m

λ_eff,i——第i个控制体的球床有效热导率/W·m^-1·K^-1

λ_eff,i+1——第i+1个控制体的球床有效热导率/W·m^-1·K^-1

λ_eff,i-1——第i-1个控制体的球床有效热导率/W·m^-1·K^-1

q″′_c,i——第i个控制体的堆芯功率密度/W·m^-3

式中：

ρ_z1,i——ZrC涂层第i个控制体的密度/kg·m^-3

c_z1,i——ZrC涂层第i个控制体的比容/J·kg^-1·K^-1

ΔV_z1,i——ZrC涂层第i个控制体的体积/m³

T_z1,i——ZrC涂层第i个控制体的温度/T

T_z1,i+1——ZrC涂层第i+1个控制体的温度/T

T_z1,i-1——ZrC涂层第i-1个控制体的温度/T

r_z1,i——ZrC涂层第i个控制体的径向坐标/m

r_z1,i+1——ZrC涂层第i+1个控制体的径向坐标/m

r_z1,i-1——ZrC涂层第i-1个控制体的径向坐标/m

λ_z1——ZrC涂层热导率/W·m^-1·K^-1

式中：

ρ_gm,i——活性区第i个控制体的密度/kg·m^-3

c_gm,i——活性区第i个控制体的比容/J·kg^-1·K^-1

ΔV_gm,i——活性区第i个控制体的体积/m³

T_gm,i——活性区第i个控制体的温度/T

T_gm,i+1——活性区第i+1个控制体的温度/T

T_gm,i-1——活性区第i-1个控制体的温度/T

r_gm,i——活性区第i个控制体的径向坐标/m

r_gm,i+1——活性区第i+1个控制体的径向坐标/m

r_gm,i-1——活性区第i-1个控制体的径向坐标/m

λ_gm——活性区热导率/W·m^-1·K^-1

q″′_gm,i——活性区第i个控制体的功率密度/W·m^-3

式中：

ρ_u,i——UC-ZrC核心第i个控制体的密度/kg·m^-3

c_u,i——UC-ZrC核心第i个控制体的比容/J·kg^-1·K^-1

ΔV_u,i——UC-ZrC核心第i个控制体的体积/m³

T_u,i——UC-ZrC核心第i个控制体的温度/T

T_u,i+1——UC-ZrC核心第i+1个控制体的温度/T

T_u,i-1——UC-ZrC核心第i-1个控制体的温度/T

r_u,i——UC-ZrC核心第i个控制体的径向坐标/m

r_u,i+1——UC-ZrC核心第i+1个控制体的径向坐标/m

r_u,i-1——UC-ZrC核心第i-1个控制体的径向坐标/m

λ_u——UC-ZrC核心热导率/W·m^-1·K^-1

q″′_u,i——UC-ZrC核心第i个控制体的功率密度/W·m^-3

式中：

ρ_i——包覆层第i个控制体的密度/kg·m^-3

c_i——包覆层第i个控制体的比容/J·kg^-1·K^-1

ΔV_i——包覆层第i个控制体的体积/m³

T_i——包覆层第i个控制体的温度/T

T_i+1——包覆层第i+1个控制体的温度/T

T_i-1——包覆层第i-1个控制体的温度/T

r_i——包覆层第i个控制体的径向坐标/m

r_i+1——包覆层第i+1个控制体的径向坐标/m

r_i-1——包覆层第i-1个控制体的径向坐标/m

根据步骤2计算的初始值，采用吉尔算法对微分方程(7)、(8)、(17)-(24)进行求解，得到t时刻的堆芯功率、温度分布，冷却剂流场、温度场和压力场；

步骤4：根据步骤3的计算结果，利用吉尔算法进行下一时刻的计算，直到达到设定的总计算时间，计算停止；

步骤5：将输出的计算结果随时间的变化趋势绘制成曲线，若曲线已趋于平缓，则计算收敛，终止计算；否则读取上一次停止时刻的计算结果作为初值，再设定计算时间和时间步长，重新启动计算，直至结果收敛。

Images