CN112699620A - 基于计算流体力学的反应堆堆芯热工水力特性分析方法 - Google Patents

Abstract

本发明公开了一种基于计算流体力学的反应堆堆芯热工水力特性分析方法，该方法主要内容如下：1、建立冷却剂流动换热特性分析模型；2、建立燃料棒热工特性分析模型；3、建立流固耦合换热模型；4、建立流固耦合求解方法。本方法通过对反应堆堆芯内复杂棒束结构实施简化处理，可以精确计算反应堆堆芯燃料棒温度分布和冷却剂流场和温度场，同时显著降低计算资源消耗，具有建模简单、计算速度快的特点，为在工程实际中运用计算流体力学数值模拟方法对反应堆全堆芯进行热工水力特性分析提供了一种新的计算方法。

Description

基于计算流体力学的反应堆堆芯热工水力特性分析方法

技术领域

本发明属于反应堆安全分析技术领域，具体涉及到正常运行及事故工况的反应堆功率变化条件下，反应堆堆芯热工水力特性分析方法。

背景技术

反应堆堆芯是核电厂一回路系统的核心部件，而堆芯燃料包壳是防止反应堆内的放射性物质向外泄露的第一道屏障，因此反应堆堆芯结构的完整性——尤其是燃料包壳的完整性对反应堆安全来说至关重要，在反应堆堆芯热工水力分析研究中，堆芯燃料棒及冷却剂温度分布一直是研究的重点。而在堆芯中核燃料及控制装置以棒束结构分布，造成堆芯内流动及传热特性的复杂化，对堆芯建模及热工水力数值模拟分析造成了较大的困难。

目前应用于反应堆热工水力分析程序中，对燃料棒温度分布计算的方法一般有三种：1.一维系统程序——如RELAP5程序、TRACE程序通过简单的热构件计算反应堆燃料的平均温度分布；2.子通道分析程序——如COBRA程序对燃料棒逐个建模，并通过数值求解得到燃料棒温度分布；3.计算流体动力学(CFD)模拟——如ANSYS FLUENT软件一般采用对燃料棒进行建模和网格划分，再逐个燃料棒求解温度分布。其中，一维系统分析程序无法获得堆芯燃料棒准确的温度分布，因此对功率分布不均匀的瞬态计算过程模拟效果很差；子通道程序建模复杂，而且在全堆模拟中，边界条件需要人为指定，对反应堆堆芯入口不规律瞬变过程模拟能力不足；CFD软件对燃料棒精确建模工作量大，网格量大，导致计算速度慢。

发明内容

本发明要解决的技术问题在于，提供一种低计算资源消耗、高计算精度的反应堆堆芯热工水力特性分析方法。

本发明解决该技术问题所采用技术方案为：提出一种基于计算流体力学的反应堆堆芯热工水力特性分析方法。针对反应堆堆芯内燃料组件复杂的棒束结构造成堆芯内部流动及传热特性复杂化的问题，本发明基于计算流体力学，建立冷却剂流动换热特性分析模型、燃料棒热工特性分析模型，建立流固耦合换热模型，通过流固耦合求解方法计算获得堆芯燃料棒温度分布和冷却剂流场和温度场，实现反应堆堆芯热工水力特性分析。

一种基于计算流体力学的反应堆堆芯热工水力特性分析方法，其特征在于：具体包括如下步骤：

步骤1：建立冷却剂流动换热特性分析模型

(1)建立反应堆堆芯的计算模型

建立反应堆堆芯的几何模型，将反应堆堆芯的几何模型导入网格生成软件中划分控制体，形成反应堆堆芯的网格模型，将反应堆堆芯的网格模型导入计算流体力学软件中，形成反应堆堆芯的计算模型；

(2)建立冷却剂流动换热特性分析模型

为了将反应堆堆芯内复杂棒束结构的影响简化为对冷却剂的流动阻力，提出如下假设：①堆芯的流动阻力使用包含粘性损失、惯性损失的压降描述：

式中：

是沿m坐标方向的动量方程中的压降，m和n是x,y,z坐标方向，D_mn和C_mn分别是粘性损失项和惯性损失项的系数矩阵，μ是冷却剂的动力粘性系数/kg·m^-1·s^-1，ρ是冷却剂的密度/kg·m^-3，|v|是速度的模/m·s^-1，v_n是沿n坐标方向的速度/m·s^-1；②计算域中燃料棒与定位格架固体部件体积不变，位置无移动；③孔隙率是各向同性的；基于上述假设，对流体控制方程中的动量方程和能量方程中添加源项：在动量方程中添加包含粘性损失、惯性损失的压降作为动量源项，在能量方程中添加燃料棒表面的对流换热量作为能量源项；从而建立了冷却剂流动换热特性分析模型；

步骤2：建立燃料棒热工特性分析模型

燃料棒热工特性分析模型中考虑燃料棒表面对流换热和燃料棒内部导热，包括表面对流换热模型和内部导热模型，分别建立表面对流换热模型和内部导热模型，步骤如下：

(1)建立表面对流换热模型

提出如下假设：①反应堆燃料棒及控制棒结构均匀地分布在堆芯中；②单个控制体内燃料棒几何结构不发生剧烈变化；③单个控制体内所有燃料棒具有相同的热工状态；④单个控制体内所有燃料棒具有相同的几何特性且已知；

基于上述假设，建立反应堆堆芯内的换热面积与反应堆堆芯体积的线性关系式：

式中：A_T,cell为控制体内的换热面积/m²，V_cell为控制体体积/m³，A_T,core为反应堆堆芯内的换热面积/m²，V_core为反应堆堆芯体积/m³；

将该线性关系式应用于步骤1中建立的反应堆堆芯的计算模型中，用控制体的体积表示控制体内的换热面积；最终将控制体内的换热面积代入到表面对流换热公式中，计算得到燃料棒表面的对流换热量：

Q_T＝h_TA_cell(T_r-T_f)

式中：Q_T为燃料棒表面的对流换热量/W，h_T为燃料棒表面的对流换热系数/W·m^-2·K^-1，A_cell为控制体内的换热面积/m²，T_f为流体主流温度/K，T_r为燃料棒表面温度/K；

通过上述步骤2的步骤(1)建立了表面对流换热模型；

(2)内部导热模型建立步骤

根据燃料棒的结构特点将燃料棒沿径向划分N个计算节点，其中外侧两个计算节点分别位于包壳内外两侧边界，其它计算节点沿燃料棒芯块中心向外侧依次布置，且有一个计算节点在芯块边界上，在包壳内外两侧各有一个计算节点；忽略燃料棒的轴向导热，沿燃料棒径向建立温度方程：

式中：ρ_i为计算节点i处材料密度/kg·m^-3，C_P,i为计算节点i处材料比热容/kJ·kg^-1·K^-1，V_i为计算节点i的等效控制体体积/m³，T_i为计算节点i处的温度/K，Q_i-1,i为从计算节点i-1传导到计算节点i的热量/kW，Q_i+1,i为计算节点i+1传导到计算节点i的热量/kW，Q_i为计算节点i处单位体积释热率/kW·m^-3；

燃料棒内部计算节点间传导的热量通过计算节点间温差得到：

Q_i-1,i＝K_i-1,i(T_i-1-T_i)

Q_i+1,i＝K_i+1,i(T_i+1-T_i)

式中：T_i-1为计算节点i-1处的温度/K，T_i+1为计算节点i+1处的温度/K，K_i-1,i为计算节点i-1和i之间的导热系数/kW·K^-1，等同于K_i,i-1，K_i+1,i为计算节点i+1和i之间的导热系数/kW·K^-1，等同于K_i,i+1，均为热导率和燃料棒几何参数的函数：

式中：k_i是计算节点i处的热导率/kW·m^-1·K^-1，k_i-1是计算节点i-1处的热导率/kW·m^-1·K^-1，k_i+1是计算节点i+1处的热导率/kW·m^-1·K^-1，r_i是计算节点i处的芯块直径/m，r_i-1是计算节点i处的芯块直径/m，ΔX_j为控制体轴向高度/m，

为计算节点i的体积平均半径/m，

为计算节点i-1的体积平均半径/m，

为计算节点i+1的体积平均半径/m，其表达式为：

根据燃料棒的导热特点，建立如下边界条件：①燃料棒中心为对称边界，因此对于计算节点1的控制方程有K_0,1＝K_i-1,i＝0；②在燃料芯块与包壳之间存在间隙导热，即在i＝N-1处Q_i,i-1＝H_gap(T_i-T_i-1)，其中H_gap为考虑换热面积的间隙导热率/kW·K^-1，Q_i,i-1为从计算节点i传导到计算节点i-1的热量/kW，T_i为计算节点i处的温度/K，T_i-1为计算节点i-1处的温度/K；③在燃料棒外表面，即在i＝N处，K_i+1,i＝H，T_i+1＝T_f，则有Q_i+1,i＝H(T_f-T_i)，其中，K_i+1,i为计算节点i+1和i之间的导热系数/kW·K^-1，H为考虑换热面积的燃料棒表面的对流换热系数/kW·K^-1，Q_i+1,i为计算节点i+1传导到计算节点i的热量/kW，T_f为流体主流温度/K，T_i+1为计算节点i+1处的温度/K，T_i为计算节点i处的温度/K；

将温度方程在各个计算节点上化为差分方程：

式中：T_i ⁿ为计算节点i处在第n个时间步下的温度/K，T_i ^n-1为计算节点i处在第n-1个时间步下的温度/K，T_i ⁿ⁺¹为计算节点i处在第n+1个时间步下的温度/K；

将差分方程变形为：

式中：系数A_i＝-K_i-1,i，系数

系数C_i＝-K_i+1,i，系数

联立各个计算节点的差分方程形成方程组，其中方程组的系数矩阵为三对角矩阵，通过高斯迭代法求解得到燃料棒内各计算节点温度；

通过上述步骤2的步骤(2)建立了内部导热模型；

步骤3：建立流固耦合换热模型

冷却剂与燃料棒之间的换热系数由Dittus-Bolter关系式计算：

Nu＝0.023Re^0.8Pr^0.4

式中：Nu为努塞尔数，Re为雷诺数，Pr为普朗特数；Nu和Re的表达式为：

式中：k为冷却剂热导率/W·m^-1·K^-1，L为特征长度/m，ρ为冷却剂密度/kg·m^-3，μ为冷却剂的动力粘性系数/kg·m^-1·s^-1，U为冷却剂流速/m·s^-1；

由于k、ρ、μ、U均为冷却剂热物性和流动状态的参数，因此冷却剂流场和温度场决定了燃料棒表面的对流换热系数，进而影响燃料棒表面的对流换热量；同时，燃料棒表面的对流换热量又影响冷却剂流场和温度场，二者存在耦合关系；因此在燃料棒表面的对流换热量的求解过程中，需对冷却剂流动换热特性分析模型和燃料棒热工特性分析模型之间进行流固耦合；

具体耦合过程为：通过燃料棒热工特性分析模型获得燃料棒温度分布的计算结果；根据冷却剂流场分布，通过上述Dittus-Bolter公式获得燃料棒表面的对流换热系数；将燃料棒表面的对流换热系数、冷却剂温度、燃料棒温度代入步骤2中的表面对流换热公式中，获得燃料棒表面的对流换热量；

通过上述步骤3建立了流固耦合换热模型；

步骤4：进行流固耦合求解

在计算开始时对步骤1中构建的反应堆堆芯的计算模型进行初始化，假定初始燃料棒温度分布、冷却剂流场和温度场；在计算过程中，燃料棒热工特性分析模型通过T_N-1时刻内的堆芯冷却剂流场和温度场获得燃料棒表面的对流换热系数，并以此作为边界条件完成T_N时刻的计算，得到T_N时刻燃料棒温度分布的计算结果；冷却剂流动换热特性分析模型根据燃料棒温度分布的计算结果，获得燃料棒表面的对流换热量，并将其添加为步骤1中的能量方程的源项进行求解，得到T_N时刻的冷却剂流场和温度场，自此完成T_N-1时刻至T_N时刻的求解；

步骤5：判断冷却剂流速、温度、燃料棒内部温度的残差是否满足收敛判据——即小于预设的残差值，若不满足则重复步骤4，直至收敛；最终获得在反应堆堆芯内燃料棒温度分布、冷却剂流场和温度场的计算结果，实现反应堆堆芯热工水力特性分析。

本发明具有以下优点和效果：

1.将堆芯内燃料棒和控制棒等复杂棒束结构简化，从而避免了逐个对燃料棒建立复杂几何模型；

2.克服了一维系统分析程序无法获得堆芯燃料棒准确的温度分布、对功率分布不均匀瞬态计算过程响应模拟较差的问题；

3.对堆芯堆芯内复杂流动状态进行简化，通过冷却剂的流动阻力反映堆芯区域中燃料棒、控制棒等固体结构的影响；

4.通过燃料棒热工特性分析模型与冷却剂流动换热特性分析模型之间的数据交换，形成流固耦合换热模型，能够在分析冷却剂流动换热特性的同时获得堆芯内燃料棒温度分布；

5.计算资源消耗少、计算速度快且计算结果精度高。

本发明已经通过实践证明，该方法能够准确获得堆芯内燃料棒温度分布、冷却剂流场和温度场的计算结果，本发明中提出的基于计算流体力学的反应堆堆芯热工水力特性分析方法可应用于三维计算流体力学数值模拟中，在显著降低计算资源消耗的同时实现反应堆堆芯燃料棒温度分布的精确计算，具有建模简单、计算速度快的特点。

附图说明

图1为反应堆堆芯燃料棒及燃料组件示意图。

图2为燃料棒热工特性分析模型计算节点划分示意图。

图3为流固耦合换热模型框图。

图4为反应堆堆芯热工水力特性分析方法框图。

具体实施方式

为了对本发明的技术特征、目的和效果有更加清楚的理解，现结合附图详细说明本发明的具体实施方式。

本发明提出了一种基于计算流体力学的反应堆堆芯热工水力特性分析方法，框架如图4所示。该方法具体实施方式如下：

步骤1：建立冷却剂流动换热特性分析模型

(1)建立反应堆堆芯计算模型

通过几何建模软件建立反应堆堆芯的几何模型，在建模时对反应堆堆芯整体建立几何模型，而不需要对单个燃料组件或单个燃料棒进行精细建模。将反应堆堆芯的几何模型导入网格生成软件——如ICEM CFD中，对几何模型沿径向和轴向划分控制体，形成反应堆堆芯的网格模型。将反应堆堆芯的网格模型导入计算流体力学软件——如Fluent中，建立反应堆堆芯的计算模型。

(2)建立冷却剂流动换热特性分析模型

为了将如图1所示的反应堆堆芯内复杂棒束结构的影响简化为对冷却剂的流动阻力，提出如下假设：①堆芯的流动阻力用包含粘性损失、惯性损失的压降描述：

式中：

是沿m坐标方向的动量方程中的压降，m和n是x,y,z坐标方向，D_mn和C_mn分别是粘性损失项和惯性损失项的系数矩阵，μ是冷却剂的动力粘性系数/kg·m^-1·s^-1，ρ是冷却剂的密度/kg·m^-3，|v|是速度的模/m·s^-1，v_n是沿n坐标方向的速度/m·s^-1；②计算域中燃料棒与定位格架固体部件体积不变，位置无移动；③孔隙率是各向同性的；基于上述假设，对流体控制方程中的动量方程和能量方程中添加源项：在动量方程中添加包含粘性损失、惯性损失的压降作为动量源项，在能量方程中添加燃料棒表面的对流换热量作为能量源项，对流换热量的具体值通过燃料棒热工特性分析模型计算；从而建立了冷却剂流动换热特性分析模型。

步骤2：建立燃料棒热工特性分析模型

燃料棒热工特性分析模型中考虑燃料棒表面对流换热和燃料棒内部导热，包括表面对流换热模型和内部导热模型，分别建立表面对流换热模型和内部导热模型，步骤如下：

(1)建立表面对流换热模型

提出如下假设：①反应堆燃料棒及控制棒结构均匀地分布在堆芯中；②单个控制体内燃料棒几何结构不发生剧烈变化；③单个控制体内所有燃料棒具有相同的热工状态；④单个控制体内所有燃料棒具有相同的几何特性且已知。

基于上述假设，反应堆堆芯内的换热面积与反应堆堆芯体积是线性相关的，建立反应堆堆芯内的换热面积与反应堆堆芯体积的线性关系式，即：

式中，A_T,cell为控制体内的换热面积/m²，V_cell为控制体体积/m³，A_T,core为反应堆堆芯内的换热面积/m²，V_core为反应堆堆芯体积/m³。

由于反应堆堆芯内的换热面积和反应堆堆芯体积是已知的，将该关系式应用于步骤1中建立的反应堆堆芯的计算模型中，用控制体的体积表示控制体内的换热面积，即：

最终将控制体内的换热面积代入到表面对流换热公式中，计算得到燃料棒表面的对流换热量：

Q_T＝h_TA_cell(T_r-T_f)

式中：Q_T为燃料棒表面的对流换热量/W，h_T为燃料棒表面的对流换热系数/W·m^-2·K^-1，A_cell为控制体内的换热面积/m²，T_f为流体主流温度/K，T_r为燃料棒表面温度/K。

通过上述步骤2的步骤(1)建立了表面对流换热模型。

(2)内部导热模型建立步骤

根据燃料棒的结构特点将燃料棒沿径向划分多个计算节点，其中外侧两个计算节点分别位于包壳内外两侧边界，其它计算节点沿燃料棒芯块中心向外侧依次布置，且有一个计算节点在芯块边界上，在包壳内外两侧各有一个计算节点；忽略燃料棒的轴向导热，沿燃料棒径向建立温度方程：

式中：ρ_i为计算节点i处材料密度/kg·m^-3，C_P,i为计算节点i处材料比热容/kJ·kg^-1·K^-1，V_i为计算节点i的等效控制体体积/m³，T_i为计算节点i处的温度/K，Q_i-1,i为从计算节点i-1传导到计算节点i的热量/kW，Q_i+1,i为计算节点i+1传导到计算节点i的热量/kW，Q_i为计算节点i处单位体积释热率/kW·m^-3。

燃料棒内部计算节点间传导的热量通过计算节点间温差得到：

Q_i-1,i＝K_i-1,i(T_i-1-T_i)

Q_i+1,i＝K_i+1,i(T_i+1-T_i)

式中：T_i-1为计算节点i-1处的温度/K，T_i+1为计算节点i+1处的温度/K，K_i-1,i为计算节点i-1和i之间的导热系数/kW·K^-1，等同于K_i,i-1，K_i+1,i为计算节点i+1和i之间的导热系数/kW·K^-1，等同于K_i,i+1，均为热导率和燃料棒几何参数的函数：

式中：k_i是计算节点i处的热导率/kW·m^-1·K^-1，k_i-1是计算节点i-1处的热导率/kW·m^-1·K^-1，k_i+1是计算节点i+1处的热导率/kW·m^-1·K^-1，r_i是计算节点i处的芯块直径/m，r_i-1是计算节点i处的芯块直径/m，ΔX_j为控制体轴向高度/m，

为计算节点i的体积平均半径/m，

为计算节点i-1的体积平均半径/m，

为计算节点i+1的体积平均半径/m，其表达式为：

根据燃料棒的导热特点，建立如下边界条件：①燃料棒中心为对称边界，因此对于计算节点1的控制方程有K_0,1＝K_i-1,i＝0；②在燃料芯块与包壳之间存在间隙导热，即在i＝N-1处Q_i,i-1＝H_gap(T_i-T_i-1)，其中H_gap为考虑换热面积的间隙导热率/kW·K^-1，Q_i,i-1为从计算节点i传导到计算节点i-1的热量/kW，T_i为计算节点i处的温度/K，T_i-1为计算节点i-1处的温度/K；③在燃料棒外表面，即在i＝N处，

K_i+1,i＝H，T_i+1＝T_f，则有Q_i+1,i＝H(T_f-T_i)，其中，K_i+1,i为计算节点i+1和i之间的导热系数/kW·K^-1，H为考虑换热面积的燃料棒表面的对流换热系数/kW·K^-1，Q_i+1,i为计算节点i+1传导到计算节点i的热量/kW，T_f为流体主流温度/K，T_i+1为计算节点i+1处的温度/K，T_i为计算节点i处的温度/K。

将温度方程在各个计算节点上化为差分方程：

式中：T_i ⁿ为计算节点i处在第n个时间步下的温度/K，T_i ^n-1为计算节点i处在第n-1个时间步下的温度/K，T_i ⁿ⁺¹为计算节点i处在第n+1个时间步下的温度/K。

将差分方程变形为：

式中：系数A_i＝-K_i-1,i，系数

系数C_i＝-K_i+1,i，系数

根据边界条件可知，i＝1时：

A₁＝0

对每个计算节点建立能量方程后联立求解，系数矩阵为：

该系数矩阵为三对角矩阵，可由高斯迭代法求解得到燃料棒内各计算节点温度。

步骤3：建立流固耦合换热模型

冷却剂与燃料棒之间的换热系数由Dittus-Bolter关系式计算：

Nu＝0.023Re^0.8Pr^0.4

式中：Nu为努塞尔数，Re为雷诺数，Pr为普朗特数；Nu和Re的表达式为：

式中：k为冷却剂热导率/W·m^-1·K^-1，L为特征长度/m，ρ为冷却剂密度/kg·m^-3，μ为冷却剂的动力粘性系数/kg·m^-1·s^-1，U为冷却剂流速/m·s^-1。

由于k、ρ、μ、U均为冷却剂热物性和流动状态的参数，因此冷却剂流场和温度场决定了燃料棒表面的对流换热系数，进而影响燃料棒表面的对流换热量；同时，燃料棒表面的对流换热量又影响冷却剂流场和温度场，二者存在耦合关系；因此在燃料棒表面的对流换热量的求解过程中，需对冷却剂流动换热特性分析模型和燃料棒热工特性分析模型之间进行流固耦合。

具体耦合过程为：通过燃料棒热工特性分析模型获得燃料棒温度分布的计算结果；根据冷却剂流场分布，通过上述Dittus-Bolter公式获得燃料棒表面的对流换热系数；将燃料棒表面的对流换热系数、冷却剂温度、燃料棒温度代入步骤2中的表面对流换热公式中，获得燃料棒表面的对流换热量。

通过上述步骤3建立了流固耦合换热模型。

步骤4：进行流固耦合求解

流固耦合求解方法如图3所示。在计算开始时对步骤1中构建的反应堆堆芯的计算模型进行初始化，假定初始燃料棒温度分布、冷却剂流场和温度场；在计算过程中，燃料棒热工特性分析模型通过T_K-1时刻内的堆芯冷却剂流场和温度场获得燃料棒表面的对流换热系数，并以此作为边界条件完成T_K时刻的计算，得到T_K时刻燃料棒温度分布的计算结果；冷却剂流动换热特性分析模型根据燃料棒温度分布的计算结果，获得燃料棒表面的对流换热量，并将其添加为步骤1中的能量方程的源项进行求解，得到T_K时刻的冷却剂流场和温度场，自此完成T_K-1时刻至T_K时刻的求解。

步骤5：判断冷却剂流速、温度、燃料棒内部温度的残差是否满足收敛判据——即小于预设的残差值，若不满足则重复步骤4，直至收敛；最终获得在反应堆堆芯内燃料棒温度分布、冷却剂流场和温度场的计算结果，实现反应堆堆芯热工水力特性分析。

Claims (1)

1.一种基于计算流体力学的反应堆堆芯热工水力特性分析方法，其特征在于：具体包括如下步骤：

步骤1：建立冷却剂流动换热特性分析模型

(1)建立反应堆堆芯的计算模型

建立反应堆堆芯的几何模型，将反应堆堆芯的几何模型导入网格生成软件中划分控制体，形成反应堆堆芯的网格模型，将反应堆堆芯的网格模型导入计算流体力学软件中，形成反应堆堆芯的计算模型；

(2)建立冷却剂流动换热特性分析模型

为了将反应堆堆芯内复杂棒束结构的影响简化为对冷却剂的流动阻力，提出如下假设：①堆芯的流动阻力使用包含粘性损失、惯性损失的压降描述：

式中：

是沿m坐标方向的动量方程中的压降，m和n是x,y,z坐标方向，D_mn和C_mn分别是粘性损失项和惯性损失项的系数矩阵，μ是冷却剂的动力粘性系数/kg·m^-1·s^-1，ρ是冷却剂的密度/kg·m^-3，|v|是速度的模/m·s^-1，v_n是沿n坐标方向的速度/m·s^-1；②计算域中燃料棒与定位格架固体部件体积不变，位置无移动；③孔隙率是各向同性的；基于上述假设，对流体控制方程中的动量方程和能量方程中添加源项：在动量方程中添加包含粘性损失、惯性损失的压降作为动量源项，在能量方程中添加燃料棒表面的对流换热量作为能量源项；从而建立了冷却剂流动换热特性分析模型；

步骤2：建立燃料棒热工特性分析模型

燃料棒热工特性分析模型中考虑燃料棒表面对流换热和燃料棒内部导热，包括表面对流换热模型和内部导热模型，分别建立表面对流换热模型和内部导热模型，步骤如下：

(1)建立表面对流换热模型

提出如下假设：①反应堆燃料棒及控制棒结构均匀地分布在堆芯中；②单个控制体内燃料棒几何结构不发生剧烈变化；③单个控制体内所有燃料棒具有相同的热工状态；④单个控制体内所有燃料棒具有相同的几何特性且已知；

基于上述假设，建立反应堆堆芯内的换热面积与反应堆堆芯体积的线性关系式：

式中：A_T,cell为控制体内的换热面积/m²，V_cell为控制体体积/m³，A_T,core为反应堆堆芯内的换热面积/m²，V_core为反应堆堆芯体积/m³；

将该线性关系式应用于步骤1中建立的反应堆堆芯的计算模型中，用控制体的体积表示控制体内的换热面积；最终将控制体内的换热面积代入到表面对流换热公式中，计算得到燃料棒表面的对流换热量：

Q_T＝h_TA_cell(T_r-T_f)

式中：Q_T为燃料棒表面的对流换热量/W，h_T为燃料棒表面的对流换热系数/W·m^-2·K^-1，A_cell为控制体内的换热面积/m²，T_f为流体主流温度/K，T_r为燃料棒表面温度/K；

通过上述步骤2的步骤(1)建立了表面对流换热模型；

(2)内部导热模型建立步骤

根据燃料棒的结构特点将燃料棒沿径向划分N个计算节点，其中外侧两个计算节点分别位于包壳内外两侧边界，其它计算节点沿燃料棒芯块中心向外侧依次布置，且有一个计算节点在芯块边界上，在包壳内外两侧各有一个计算节点；忽略燃料棒的轴向导热，沿燃料棒径向建立温度方程：

式中：ρ_i为计算节点i处材料密度/kg·m^-3，C_P,i为计算节点i处材料比热容/kJ·kg^-1·K^-1，V_i为计算节点i的等效控制体体积/m³，T_i为计算节点i处的温度/K，Q_i-1,i为从计算节点i-1传导到计算节点i的热量/kW，Q_i+1,i为计算节点i+1传导到计算节点i的热量/kW，Q_i为计算节点i处单位体积释热率/kW·m^-3；

燃料棒内部计算节点间传导的热量通过计算节点间温差得到：

Q_i-1,i＝K_i-1,i(T_i-1-T_i)

Q_i+1,i＝K_i+1,i(T_i+1-T_i)

式中：T_i-1为计算节点i-1处的温度/K，T_i+1为计算节点i+1处的温度/K，K_i-1,i为计算节点i-1和i之间的导热系数/kW·K^-1，等同于K_i,i-1，K_i+1,i为计算节点i+1和i之间的导热系数/kW·K^-1，等同于K_i,i+1，均为热导率和燃料棒几何参数的函数：

式中：k_i是计算节点i处的热导率/kW·m^-1·K^-1，k_i-1是计算节点i-1处的热导率/kW·m^-1·K^-1，k_i+1是计算节点i+1处的热导率/kW·m^-1·K^-1，r_i是计算节点i处的芯块直径/m，r_i-1是计算节点i处的芯块直径/m，ΔX_j为控制体轴向高度/m，

为计算节点i的体积平均半径/m，

为计算节点i-1的体积平均半径/m，

为计算节点i+1的体积平均半径/m，其表达式为：

根据燃料棒的导热特点，建立如下边界条件：①燃料棒中心为对称边界，因此对于计算节点1的控制方程有K_0,1＝K_i-1,i＝0；②在燃料芯块与包壳之间存在间隙导热，即在i＝N-1处Q_i,i-1＝H_gap(T_i-T_i-1)，其中H_gap为考虑换热面积的间隙导热率/kW·K^-1，Q_i,i-1为从计算节点i传导到计算节点i-1的热量/kW，T_i为计算节点i处的温度/K，T_i-1为计算节点i-1处的温度/K；③在燃料棒外表面，即在i＝N处，K_i+1,i＝H，T_i+1＝T_f，则有Q_i+1,i＝H(T_f-T_i)，其中，K_i+1,i为计算节点i+1和i之间的导热系数/kW·K^-1，H为考虑换热面积的燃料棒表面的对流换热系数/kW·K^-1，Q_i+1,i为计算节点i+1传导到计算节点i的热量/kW，T_f为流体主流温度/K，T_i+1为计算节点i+1处的温度/K，T_i为计算节点i处的温度/K；

将温度方程在各个计算节点上化为差分方程：

式中：T_i ⁿ为计算节点i处在第n个时间步下的温度/K，T_i ^n-1为计算节点i处在第n-1个时间步下的温度/K，T_i ⁿ⁺¹为计算节点i处在第n+1个时间步下的温度/K；

将差分方程变形为：

式中：系数A_i＝-K_i-1,i，系数

系数C_i＝-K_i+1,i，系数

联立各个计算节点的差分方程形成方程组，其中方程组的系数矩阵为三对角矩阵，通过高斯迭代法求解得到燃料棒内各计算节点温度；

通过上述步骤2的步骤(2)建立了内部导热模型；

步骤3：建立流固耦合换热模型

冷却剂与燃料棒之间的换热系数由Dittus-Bolter关系式计算：

Nu＝0.023Re^0.8Pr^0.4

式中：Nu为努塞尔数，Re为雷诺数，Pr为普朗特数；Nu和Re的表达式为：

式中：k为冷却剂热导率/W·m^-1·K^-1，L为特征长度/m，ρ为冷却剂密度/kg·m^-3，μ为冷却剂的动力粘性系数/kg·m^-1·s^-1，U为冷却剂流速/m·s^-1；

由于k、ρ、μ、U均为冷却剂热物性和流动状态的参数，因此冷却剂流场和温度场决定了燃料棒表面的对流换热系数，进而影响燃料棒表面的对流换热量；同时，燃料棒表面的对流换热量又影响冷却剂流场和温度场，二者存在耦合关系；因此在燃料棒表面的对流换热量的求解过程中，需对冷却剂流动换热特性分析模型和燃料棒热工特性分析模型之间进行流固耦合；

具体耦合过程为：通过燃料棒热工特性分析模型获得燃料棒温度分布的计算结果；根据冷却剂流场分布，通过上述Dittus-Bolter公式获得燃料棒表面的对流换热系数；将燃料棒表面的对流换热系数、冷却剂温度、燃料棒温度代入步骤2中的表面对流换热公式中，获得燃料棒表面的对流换热量；

通过上述步骤3建立了流固耦合换热模型；

步骤4：进行流固耦合求解

在计算开始时对步骤1中构建的反应堆堆芯的计算模型进行初始化，假定初始燃料棒温度分布、冷却剂流场和温度场；在计算过程中，燃料棒热工特性分析模型通过T_N-1时刻内的堆芯冷却剂流场和温度场获得燃料棒表面的对流换热系数，并以此作为边界条件完成T_N时刻的计算，得到T_N时刻燃料棒温度分布的计算结果；冷却剂流动换热特性分析模型根据燃料棒温度分布的计算结果，获得燃料棒表面的对流换热量，并将其添加为步骤1中的能量方程的源项进行求解，得到T_N时刻的冷却剂流场和温度场，自此完成T_N-1时刻至T_N时刻的求解；

步骤5：判断冷却剂流速、温度、燃料棒内部温度的残差是否满足收敛判据——即小于预设的残差值，若不满足则重复步骤4，直至收敛；最终获得在反应堆堆芯内燃料棒温度分布、冷却剂流场和温度场的计算结果，实现反应堆堆芯热工水力特性分析。

Images