CN111143996B - 一种快堆六边形闭式组件盒间流动换热计算方法 - Google Patents

Abstract

本发明公开了一种快堆六边形闭式组件盒间流动换热计算方法，主要包括以下步骤：1、对闭式组件的盒间通道进行分层；2、建立盒间通道内冷却剂流动换热的控制方程；3、采用交错网格进行网格划分；4、采用有限容积积分法对方程进行离散；5、采用SIMPLE算法求解连续方程与动量方程，采用Gear算法求解能量方程；6、对压力和速度的修正值进行亚松弛处理加快迭代收敛。本发明解决了以往只考虑盒间通道内冷却剂的导热所带来的计算不准确问题，或者是采用三维CFD方法需要大量计算机内存和计算时间的问题。本发明的计算方法可以与快堆系统分析程序进行耦合求解，更加有效地分析快堆六边形闭式组件盒间的流动换热情况，从而为快堆的设计与安全分析提供依据。

Description

一种快堆六边形闭式组件盒间流动换热计算方法

技术领域

本发明涉及快中子反应堆堆芯热工水力设计领域，具体涉及一种快堆六边形闭式组件盒间流动换热计算方法。

背景技术

快中子反应堆的堆芯由数量众多的六边形闭式组件排布而成，堆芯内的流道数量多、功率分布情况复杂，涉及复杂的流动换热问题。六边形闭式组件内的燃料棒束以正三角形栅格的形式排列，燃料棒束外由六边形闭式组件盒包围，组件与组件之间存在一定的间隙并以蜂窝状排布成堆芯。六边形闭式组件盒间隙的冷却剂会通过热传导与对流换热影响堆芯的径向与轴向的温度分布，从而影响堆芯的热工水力特性。因此准确计算快堆六边形闭式组件盒间隙的流动换热特性对快堆的设计与安全分析具有重要意义。

目前国际上针对快堆开展了大量的实验，实验主要关注核反应堆系统的热工水力特性，而对于组件盒间的流动换热的研究较少。针对快堆系统，国际上开发了许多快堆的系统分析程序。这些分析程序可以对快堆的系统进行计算，但是对于组件盒间流动换热对堆芯温度分布的影响都无法很好地模拟。美国阿贡国家实验室开发的快堆系统分析程序SAS4A/SASSYS-1将盒间流动考虑为结构之间的导热，日本原子能机构开发的快堆系统分析程序NETFLOW将盒间的流体考虑为静止并用多层模型计算组件盒间的导热，这两个程序都无法模拟出组件盒间的流动特征。对于组件盒间流动，采用三维CFD程序的精细建模由于组件数量多、盒间尺寸小而存在困难；并且如果对组件盒间流动进行三维模拟，在与一维系统程序耦合计算过程中计算速度会非常慢。因此需要开发一种既可以模拟快堆六边形闭式组件盒间流动换热特性，又不会耗费大量计算时间与计算机内存的计算方法，为准确、有效地进行快堆设计与安全分析提供依据。

发明内容

为了克服上述现有技术存在的问题，本发明提供了一种快堆六边形闭式组件盒间流动换热的计算方法，解决了现有计算方法无法模拟出组件盒间流动特征的问题；且该方法在与系统程序耦合进行计算时可以保证程序的计算速度，解决了三维CFD程序与一维系统分析程序耦合计算时耗费大量计算时间与计算机内存的问题。

为了实现上述目的，本发明具体包括以下步骤：

步骤1：根据快堆六边形闭式组件的布置方式将快堆六边形闭式组件的盒间通道进行分层，自内向外相邻两层组件的组件盒外壁面所形成的流动通道为一层；

步骤2：为了计算各层之间的流动换热机理，进行如下假设：

(1)流体粘性应力引起的耗散远小于燃料组件与盒间流体之间的换热，忽略不计；

(2)流体间轴向热传导远小于燃料组件与盒间流体之间的换热，忽略不计；

(3)冷却剂是不可压缩流体，冷却剂的密度只与温度有关；

(4)忽略组件盒间流动的周向不均匀性，简化为二维模型；

基于以上假设，建立组件盒间流动的连续方程、能量方程、轴向动量方程以及横向动量方程：

1)连续方程

式中：

t——时间/s；

i——通道i；

j——通道j；

z——通道i轴向高度/m；

r——通道i与j径向距离/m；

A_i——通道i流通截面积/m²；

ρ_i——通道i流体密度/kg·m^-3；

m_i——通道i轴向质量流量/kg·s^-1；

w_ij——通道i与j间横向质量流量/kg·s^-1；

2)能量方程

式中：

h——控制体流体比焓/J·kg^-1；

m——控制体流体质量/kg；

A——控制体与组件盒换热面积/m²；

K——组件盒间流体与组件盒对流换热系数/W·m^-2·K^-1；

T_w——组件盒壁温度/K；

T——控制体流体温度/K；

h_i——控制体或相邻控制体流体比焓/J·kg^-1；

w_i——穿过控制体界面的质量流量/kg·s^-1；

3)轴向动量方程

式中：

m_i——通道i轴向质量流量/kg·s^-1；

u_i——控制体界面的轴向速度/m·s^-1；

w_ij——通道i与j间径向质量流量/kg·s^-1；

v_ij——通道i与j间径向速度/m·s^-1；

P_i——通道i的压力/Pa；

g——重力加速度/m·s^-2；

f_i——通道i摩擦阻力系数；

D_i——通道i的轴向等效水力直径/m；

K_s——局部阻力系数；

Δz——轴向控制体高度/m；

4)横向动量方程

式中：

A_ij——横向流通面积；

f_ij——横向流动摩擦阻力系数；

D_ij——横向流动等效水力直径；

K_G——通道i与j间横向流动局部阻力系数；

Δr——径向控制体长度/m；

步骤3：为了避免流场中出现不合理的压力场，同时为了保证计算的准确度，组件盒间流动通道网格划分采用交错网格；此处的交错网格是指把轴向速度u和径向速度v以及压力P分布储存于三套不同的网格上；压力P所在的控制体为主控制体，轴向速度u存在于主控制体的南北界面上，称为轴向动量控制体，径向速度v存在于主控制体的东西界面上，称为径向动量控制体；若主控制体划分为N个区域，则主控制体节点数目为N+2，速度控制体的节点编号从2开始到N+2；

步骤4：盒间流的流道结构复杂，不同于简单的规则形状的计算，其径向流通面积、轴向流通面积、主控制体容积以及组件与盒间的换热面积的计算的正确性影响方程的质量守恒计算和能量守恒计算，进而影响矩阵求解的收敛性，因此能够准确的计算几何参数至关重要；

组件数目求解：

反应堆堆芯组件按照六边形布置，每层组件的组件数目M比内层均多6个：

第1层：M＝1

第i层：M＝6(i-1)

组件盒间的间隙面积A_axi求解：

式中：

L——六边形闭式组件盒边长/m；

S——组件盒间隙宽度/m；

径向流通面积求解：

每个盒间间隙的径向流通面积A_radial等于组件盒间隙宽度S乘以控制体轴向长度d_cv：

A_radial＝S×d_cv 公式(6)

步骤5：由于盒间结构不规则，为保证所得差分方程具有守恒特性，在此选用有限容积积分法，离散格式选用隐式迎风格式，壁面边界节点采用无滑移边界条件，控制方程组中各个方程的离散方程如下：

1)连续方程离散：

ρ_j,k-1/2A_j,k-1/2u_j,k-ρ_j,k+1/2A_j,k+1/2u_j,k+1+ρ_j+1/2,kS_j+1/2,kv_j+1,k-ρ_j-1/2,kS_j-1/2,kv_j,k＝0

公式(7)

式中：

ρ_j,k-1/2——第j通道，第k-1/2个轴向控制体的密度/kg/m³；

A_j,k-1/2——第j通道，第k-1/2个轴向控制体的轴向流通面积/m²；

S_j+1/2,k——第j+1/2通道，第k个轴向控制体的径向流通面积/m²；

u_j,k——第j通道，第k个轴向控制体的轴向流速/m/s；

v_j,k——第j通道，第k个轴向控制体的横向流速/m/s；

2)轴向动量方程离散：

公式(8)中的系数

和轴向动量源项

的表达式如下：

公式(8)中

和

分别表示一个轴向动量控制体的上、下、外和内界面位置的质量流量，其表达式为：

式中：

ΔV——轴向动量控制体体积/m³；

Δt——时间步长/s；

P_j,k——第j通道、第k轴向控制体的压力/Pa；

f_j,k——第j通道、第k轴向控制体的摩擦阻力系数；

Δz——轴向动量控制体长度/m；

——轴向动量控制体的等效水力直径/m；

k_j,k——第j通道、第k轴向控制体的局部阻力系数；

——第j通道、第k轴向控制体的轴向初始速度/m/s；

u_j,k——第j通道、第k轴向控制体的轴向速度/m/s；

v_j+1,k——第j+1通道、第k轴向控制体的径向速度/m/s；

A_j,k——第j通道、第k轴向控制体的轴向流通面积/m²；

S_j+1/2,k——第j+1/2通道、第k轴向控制体的径向流通面积/m²；

α_z——第k个主控制体和第k个轴向动量控制体重叠部分所占的比例；

β_z——第k-1个主控制体和第k个轴向动量控制体重叠部分所占的比例；

3)径向动量方程离散：

公式(9)中的系数

和径向动量源项

的表达式如下：

公式(9)中

和

分别表示一个径向动量控制体的上、下、外和内界面位置的质量流量，其表达式为：

式中：

Δr——径向动量控制体长度/m；

——径向动量控制体等效水力直径/m；

——第j通道、第k轴向控制体的轴向初始速度/m/s；

u_j,k+1——第j通道、第k+1轴向控制体的轴向速度/m/s；

v_j+1,k——第j+1通道、第k轴向控制体的径向速度/m/s；

A_j,k+1/2——第j通道、第k+1/2轴向控制体的轴向流通面积/m²；

S_j,k——第j通道、第k轴向控制体的径向流通面积/m²；

α_r——第j通道主控制体和第j通道动量控制体重叠部分所占的比例；

β_r——第j-1通道主控制体和第j通道动量控制体重叠部分所占的比例；

4)能量方程离散：

六边形闭式组件盒间流动的能量方程的求解是将方程离散成焓值的时间导数的表达形式，在整个堆芯及系统内耦合全场求解，离散后的能量方程为：

式中：

m_j,k——第j通道、第k轴向控制体的质量/kg；

h_j.k——第j通道、第k轴向控制体流体的比焓/J·kg^-1；

Q_j.k——组件盒间流体与组件或者是钢套屏蔽之间的换热功率/J/s；

公式(10)中G_N、G_S、G_O和G_IN分别表示一个主控制体的上、下、外和内界面位置的质量流量，其表达式为：

G_N＝ρ_j,k+1/2u_j,k+1A_j,k+1；

G_S＝ρ_j,k-1/2u_j,kA_j,k；

G_O＝ρ_j+1,kv_j+1,kS_j+1,k；

G_IN＝ρ_j-1,kv_j,kS_j,k；

步骤6：组件盒间流动的连续方程与动量方程采用压力修正法求解，对于能量方程则采用吉尔算法进行求解；其中组件盒间流动的连续方程与动量方程的瞬态计算的基本步骤如下：

1.读入或者传入上一时刻的轴向初始速度u⁰、径向初始速度v⁰和初始压力P⁰；

2.传入边界点的压力值；

3.根据上一时层的压力场和速度场以及压力边界求解动量方程，求解轴向速度修正值u^*、径向速度修正值v^*、压力修正值P^*；

4.求解压力修正方程，获得修正后的压力P′；

5.根据修正后的压力P′修正速度，获得修正后的轴向速度u′、修正后的径向速度v′；

6.利用改进后的速度场求解与当前速度和压力相关的物性和源项，重复步骤3-6，直至速度场收敛；

7.重复步骤1-6进行下一个时层的计算；

其中轴向速度与径向速度的修正方程如下：

式中：

u′_j,k——第j通道、第k轴向控制体修正后的轴向速度/m/s；

v′_j,k——第j通道、第k轴向控制体修正后的径向速度/m/s；

P′_j,k——第j通道、第k轴向控制体修正后的压力/Pa；

根据压力修正值P′改进的速度场应该满足连续性方程，将公式(11)、公式(12)代入连续方程得：

式中：

——第j通道、第k轴向控制体轴向速度修正值/m/s；

——第j通道、第k轴向控制体径向速度修正值/m/s；

整理得压力修正值的计算方程：

B_Pp′_j,k＝B_sp′_j,k-1+B_np′_j,k+1+B_op′_j+1,k+B_inp′_j-1,k+S_P 公式(14)

公式(14)中的系数B_P、B_s、B_n、B_o、B_in和源项S_P的表达式如下：

B_P＝B_s+B_n+B_o+B_in

式中：

P′_j,k——第j通道、第k轴向控制体修正后的压力/Pa；

——第j通道、第k轴向控制体轴向速度修正值/m/s；

——第j通道、第k轴向控制体径向速度修正值/m/s；

其中在边界上的压力修正值均为0；

步骤7：对压力和速度的修正值进行亚松弛处理从而加快计算的收敛，其中对于压力修正值直接进行亚松弛处理：

p＝α_pp′+p^* 公式(15)

式中：

α_p——压力亚松弛因子；

P——亚松弛处理后的压力；

为了加快非线性问题的迭代收敛，对速度也进行亚松弛处理；对于速度的亚松弛处理是将亚松弛过程组织到代数方程的求解过程中，因此动量方程亚松弛后的方程式为：

式中：

α_ra——径向亚松弛因子；

α_axial——轴向亚松弛因子；

u_nb——相邻控制体轴向速度/m/s；

v_nb——相邻控制体径向速度/m/s；

——由轴向相邻控制体进入的质量流量/kg/s；

——由径向相邻控制体进入的质量流量/kg/s；

u′_j,k——上一时层计算得到的轴向速度/m/s；

v′_j,k——上一时层计算得到的径向速度/m/s；

——本时层轴向的初场；

——本时层径向的初场；

终止每个层次流场迭代求解收敛的判据为连续性方程的余量范数与出口流量的比值小于等于收敛误差ε，表示为：

式中：

S_P——各控制体上连续方程的源项；

q_m——出口质量流量/kg·s^-1；

ε——收敛误差；

其中收敛误差ε建议的取值范围为：1.0e-3≤ε≤1.0e-5。

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

1)本方法将相邻两层六边形闭式组件之间的流道作为一层，考虑了层与层之间的流动换热以及冷却剂与六边形闭式组件盒壁面的换热，可以分析快堆六边形闭式组件盒间的冷却剂的流动换热对反应堆堆芯的径向与轴向温度分布的影响。

2)本方法可以与快堆系统分析程序进行耦合，并对快堆系统进行瞬态计算。在计算过程中，该方法计算速度快，不需要大量的计算机内存，可以满足快堆的设计与安全分析的需求。

附图说明

图1是本计算方法的计算流程图。

图2是快堆六边形闭式组件盒间流动通道径向分层示意图。

图3是快堆六边形闭式组件盒间流动通道轴向控制体划分示意图。

图4是快堆六边形闭式组件盒间流动通道径向控制体划分示意图。

图5是六边形闭式组件盒间间隙面积求解示意图。

图6是盒间与组件间换热示意图。

具体实施方式

以下结合附图对本发明的计算方法作进一步详细说明：

本发明提供了一种快堆六边形闭式组件盒间流动换热的计算方法，该方法的计算流程如图1所示，其主要包含以下步骤：

步骤1：如图2所示，根据快堆六边形闭式组件的布置方式将快堆六边形闭式组件的盒间通道进行分层，自内向外相邻两层组件的组件盒外壁面所形成的流动通道为一层，图中标注的是第一层至第三层；

步骤2：为了计算各层之间的流动换热机理，进行如下假设：

(1)流体粘性应力引起的耗散远小于燃料组件与盒间流体之间的换热，忽略不计；

(2)流体间轴向热传导远小于燃料组件与盒间流体之间的换热，忽略不计；

(3)冷却剂是不可压缩流体，冷却剂的密度只与温度有关；

(4)忽略组件盒间流动的周向不均匀性，简化为二维模型；

基于以上假设，建立组件盒间流动的连续方程、能量方程、轴向动量方程以及横向动量方程：

1)连续方程

式中：

t——时间/s；

i——通道i；

j——通道j；

z——通道i轴向高度/m；

r——通道i与j径向距离/m；

A_i——通道i流通截面积/m²；

ρ_i——通道i流体密度/kg·m^-3；

m_i——通道i轴向质量流量/kg·s^-1；

w_ij——通道i与j间横向质量流量/kg·s^-1；

2)能量方程

式中：

h——控制体流体比焓/J·kg^-1；

m——控制体流体质量/kg；

A——控制体与组件盒换热面积/m²；

K——组件盒间流体与组件盒对流换热系数/W·m^-2·K^-1；

T_w——组件盒壁温度/K；

T——控制体流体温度/K；

h_i——控制体或相邻控制体流体比焓/J·kg^-1；

w_i——穿过控制体界面的质量流量/kg·s^-1；

3)轴向动量方程

式中：

m_i——通道i轴向质量流量/kg·s^-1；

u_i——控制体界面的轴向速度/m·s^-1；

w_ij——通道i与j间径向质量流量/kg·s^-1；

v_ij——通道i与j间径向速度/m·s^-1；

P_i——通道i的压力/Pa；

g——重力加速度/m·s^-2；

f_i——通道i摩擦阻力系数；

D_i——通道i的轴向等效水力直径/m；

K_s——局部阻力系数；

Δz——轴向控制体高度/m；

4)横向动量方程

式中：

A_ij——横向流通面积；

f_ij——横向流动摩擦阻力系数；

D_ij——横向流动等效水力直径；

K_G——通道i与j间横向流动局部阻力系数；

Δr——径向控制体长度/m；

步骤3：为了避免流场中出现不合理的压力场，同时为了保证计算的准确度，组件盒间流动通道网格划分采用交错网格；此处的交错网格是指把轴向速度u和径向速度v以及压力P分布储存于三套不同的网格上；压力P所在的控制体为主控制体，轴向速度u存在于主控制体的南北界面上，称为轴向动量控制体，径向速度v存在于主控制体的东西界面上，称为径向动量控制体；如图3所示，将轴向划分为N个主控制体，则主控制体节点数目为N+2，速度控制体的节点编号从2开始到N+2；对于径向方向，控制体划分如图4所示，两层组件之间的一层盒间流道为一个控制体。第一层围板与组件盒壁之间为一个控制体，共有N_baff层围板因此径向共划分为N+N_baff+1个控制体；

步骤4：盒间流的流道结构复杂，不同于简单的规则形状的计算，其径向流通面积、轴向流通面积、主控制体容积以及组件与盒间的换热面积的计算的正确性影响方程的质量守恒计算和能量守恒计算，进而影响矩阵求解的收敛性，因此能够准确的计算几何参数至关重要；

组件数目求解：

反应堆堆芯组件按照六边形布置，每层组件的组件数目M比内层均多6个：

第1层：M＝1

第i层：M＝6(i-1)

如图5所示，图中标注部分为组件盒间间隙的一部分，组件盒间间隙总面积与该部分面积呈倍数关系，该部分的面积A_axi可由以下公式求解：

式中：

L——六边形闭式组件盒边长/m；

S——组件盒间隙宽度/m；

径向流通面积求解：

每个盒间间隙的径向流通面积A_radial等于组件盒间隙宽度S乘以控制体轴向长度d_cv：

A_radial＝S×d_cv 公式(6)

步骤5：由于盒间结构不规则，为保证所得差分方程具有守恒特性，在此选用有限容积积分法，离散格式选用隐式迎风格式，壁面边界节点采用无滑移边界条件，控制方程组中各个方程的离散方程如下：

1)连续方程离散：

ρ_j,k-1/2A_j,k-1/2u_j,k-ρ_j,k+1/2A_j,k+1/2u_j,k+1+ρ_j+1/2,kS_j+1/2,kv_j+1,k-ρ_j-1/2,kS_j-1/2,kv_j,k＝0

公式(7)

式中：

ρ_j,k-1/2——第j通道，第k-1/2个轴向控制体的密度/kg/m³；

A_j,k-1/2——第j通道，第k-1/2个轴向控制体的轴向流通面积/m²；

S_j+1/2,k——第j+1/2通道，第k个轴向控制体的径向流通面积/m²；

u_j,k——第j通道，第k个轴向控制体的轴向流速/m/s；

v_j,k——第j通道，第k个轴向控制体的横向流速/m/s；

2)轴向动量方程离散：

公式(8)中的系数

和轴向动量源项

的表达式如下：

公式(8)中

和

分别表示一个轴向动量控制体的上、下、外和内界面位置的质量流量，其表达式为：

式中：

ΔV——轴向动量控制体体积/m³；

Δt——时间步长/s；

P_j,k——第j通道、第k轴向控制体的压力/Pa；

f_j,k——第j通道、第k轴向控制体的摩擦阻力系数；

Δz——轴向动量控制体长度/m；

——轴向动量控制体的等效水力直径/m；

k_j,k——第j通道、第k轴向控制体的局部阻力系数；

——第j通道、第k轴向控制体的轴向初始速度/m/s；

u_j,k——第j通道、第k轴向控制体的轴向速度/m/s；

v_j+1,k——第j+1通道、第k轴向控制体的径向速度/m/s；

A_j,k——第j通道、第k轴向控制体的轴向流通面积/m²；

S_j+1/2,k——第j+1/2通道、第k轴向控制体的径向流通面积/m²；

α_z——第k个主控制体和第k个轴向动量控制体重叠部分所占的比例；

β_z——第k-1个主控制体和第k个轴向动量控制体重叠部分所占的比例；

3)径向动量方程离散：

公式(9)中的系数

和径向动量源项

的表达式如下：

公式(9)中

和

分别表示一个径向动量控制体的上、下、外和内界面位置的质量流量，其表达式为：

式中：

Δr——径向动量控制体长度/m；

——径向动量控制体等效水力直径/m；

——第j通道、第k轴向控制体的轴向初始速度/m/s；

u_j,k+1——第j通道、第k+1轴向控制体的轴向速度/m/s；

v_j+1,k——第j+1通道、第k轴向控制体的径向速度/m/s；

A_j,k+1/2——第j通道、第k+1/2轴向控制体的轴向流通面积/m²；

S_j,k——第j通道、第k轴向控制体的径向流通面积/m²；

α_r——第j通道主控制体和第j通道动量控制体重叠部分所占的比例；

β_r——第j-1通道主控制体和第j通道动量控制体重叠部分所占的比例；

4)能量方程离散：

六边形闭式组件盒间流体与组件之间的换热形式如图6所示，箭头方向代表组件盒间与组件之间的换热。组件盒间流动的能量方程的求解是将方程离散成焓值的时间导数的表达形式，在整个堆芯及系统内耦合全场求解，离散后的能量方程为：

式中：

m_j,k——第j通道、第k轴向控制体的质量/kg；

h_j.k——第j通道、第k轴向控制体流体的比焓/J·kg^-1；

Q_j.k——组件盒间流体与组件或者是钢套屏蔽之间的换热功率/J/s；

公式(10)中G_N、G_S、G_O和G_IN分别表示一个主控制体的上、下、外和内界面位置的质量流量，其表达式为：

G_N＝ρ_j,k+1/2u_j,k+1A_j,k+1；

G_S＝ρ_j,k-1/2u_j,kA_j,k；

G_O＝ρ_j+1,kv_j+1,kS_j+1,k；

G_IN＝ρ_j-1,kv_j,kS_j,k；

步骤6：组件盒间流动的连续方程与动量方程采用压力修正法求解，对于能量方程则采用吉尔算法进行求解；其中组件盒间流动的连续方程与动量方程的瞬态计算的基本步骤如下：

1.读入或者传入上一时刻的轴向初始速度u⁰、径向初始速度v⁰和初始压力P⁰；

2.传入边界点的压力值；

3.根据上一时层的压力场和速度场以及压力边界求解动量方程，求解轴向速度修正值u^*、径向速度修正值v^*、压力修正值P^*；

4.求解压力修正方程，获得修正后的压力P′；

5.根据修正后的压力P′修正速度，获得修正后的轴向速度u′、修正后的径向速度v′；

6.利用改进后的速度场求解与当前速度和压力相关的物性和源项，重复步骤3-6，直至速度场收敛；

7.重复步骤1-6进行下一个时层的计算；

其中轴向速度与径向速度的修正方程如下：

式中：

u′_j,k——第j通道、第k轴向控制体修正后的轴向速度/m/s；

v′_j,k——第j通道、第k轴向控制体修正后的径向速度/m/s；

P′_j,k——第j通道、第k轴向控制体修正后的压力/Pa；

根据压力修正值P′改进的速度场应该满足连续性方程，将公式(11)、公式(12)代入连续方程得：

式中：

——第j通道、第k轴向控制体轴向速度修正值/m/s；

——第j通道、第k轴向控制体径向速度修正值/m/s；

整理得压力修正值的计算方程：

B_Pp′_j,k＝B_sp′_j,k-1+B_np′_j,k+1+B_op′_j+1,k+B_inp′_j-1,k+S_P 公式(14)

公式(14)中的系数B_P、B_s、B_n、B_o、B_in和源项S_P的表达式如下：

B_P＝B_s+B_n+B_o+B_in

式中：

P′_j,k——第j通道、第k轴向控制体修正后的压力/Pa；

——第j通道、第k轴向控制体轴向速度修正值/m/s；

——第j通道、第k轴向控制体径向速度修正值/m/s；

其中在边界上的压力修正值均为0；

步骤7：对压力和速度的修正值进行亚松弛处理从而加快计算的收敛，其中对于压力修正值直接进行亚松弛处理：

p＝α_pp′+p^* 公式(15)

式中：

α_p——压力亚松弛因子；

P——亚松弛处理后的压力；

为了加快非线性问题的迭代收敛，对速度也进行亚松弛处理；对于速度的亚松弛处理是将亚松弛过程组织到代数方程的求解过程中，因此动量方程亚松弛后的方程式为：

式中：

α_ra——径向亚松弛因子；

α_axial——轴向亚松弛因子；

u_nb——相邻控制体轴向速度/m/s；

v_nb——相邻控制体径向速度/m/s；

——由轴向相邻控制体进入的质量流量/kg/s；

——由径向相邻控制体进入的质量流量/kg/s；

u′_j,k——上一时层计算得到的轴向速度/m/s；

v′_j,k——上一时层计算得到的径向速度/m/s；

——本时层轴向的初场；

——本时层径向的初场；

终止每个层次流场迭代求解收敛的判据为连续性方程的余量范数与出口流量的比值小于等于收敛误差ε，表示为：

式中：

S_P——各控制体上连续方程的源项；

q_m——出口质量流量/kg·s^-1；

ε——收敛误差；

其中收敛误差ε建议的取值范围为：1.0e-3≤ε≤1.0e-5。

Claims (1)

1.一种快堆六边形闭式组件盒间流动换热计算方法，其特征在于：步骤如下：

步骤1：根据快堆六边形闭式组件的布置方式将快堆六边形闭式组件的盒间通道进行分层，自内向外相邻两层组件的组件盒外壁面所形成的流动通道为一层；

步骤2：为了计算各层之间的流动换热机理，进行如下假设：

(1)流体粘性应力引起的耗散远小于燃料组件与盒间流体之间的换热，忽略不计；

(2)流体间轴向热传导远小于燃料组件与盒间流体之间的换热，忽略不计；

(3)冷却剂是不可压缩流体，冷却剂的密度只与温度有关；

(4)忽略组件盒间流动的周向不均匀性，简化为二维模型；

基于以上假设，建立组件盒间流动的连续方程、能量方程、轴向动量方程以及横向动量方程：

1)连续方程

式中：

t——时间/s；

i——通道i；

j——通道j；

z——通道i轴向高度/m；

r——通道i与j径向距离/m；

A_i——通道i流通截面积/m²；

ρ_i——通道i流体密度/kg·m^-3；

m_i——通道i轴向质量流量/kg·s^-1；

w_ij——通道i与j间横向质量流量/kg·s^-1；

2)能量方程

式中：

h——控制体流体比焓/J·kg^-1；

m——控制体流体质量/kg；

A——控制体与组件盒换热面积/m²；

K——组件盒间流体与组件盒对流换热系数/W·m^-2·K^-1；

T_w——组件盒壁温度/K；

T——控制体流体温度/K；

h_i——控制体或相邻控制体流体比焓/J·kg^-1；

w_i——穿过控制体界面的质量流量/kg·s^-1；

3)轴向动量方程

式中：

m_i——通道i轴向质量流量/kg·s^-1；

u_i——控制体界面的轴向速度/m·s^-1；

w_ij——通道i与j间径向质量流量/kg·s^-1；

v_ij——通道i与j间径向速度/m·s^-1；

P_i——通道i的压力/Pa；

g——重力加速度/m·s^-2；

f_i——通道i摩擦阻力系数；

D_i——通道i的轴向等效水力直径/m；

K_s——局部阻力系数；

Δz——轴向控制体高度/m；

4)横向动量方程

式中：

A_ij——横向流通面积；

f_ij——横向流动摩擦阻力系数；

D_ij——横向流动等效水力直径；

K_G——通道i与j间横向流动局部阻力系数；

Δr——径向控制体长度/m；

步骤3：为了避免流场中出现不合理的压力场，同时为了保证计算的准确度，组件盒间流动通道网格划分采用交错网格；此处的交错网格是指把轴向速度u和径向速度v以及压力P分布储存于三套不同的网格上；压力P所在的控制体为主控制体，轴向速度u存在于主控制体的南北界面上，称为轴向动量控制体，径向速度v存在于主控制体的东西界面上，称为径向动量控制体；若主控制体划分为N个区域，则主控制体节点数目为N+2，速度控制体的节点编号从2开始到N+2；

步骤4：盒间流的流道结构复杂，不同于简单的规则形状的计算，其径向流通面积、轴向流通面积、主控制体容积以及组件与盒间的换热面积的计算的正确性影响方程的质量守恒计算和能量守恒计算，进而影响矩阵求解的收敛性，因此能够准确的计算几何参数至关重要；

组件数目求解：

反应堆堆芯组件按照六边形布置，每层组件的组件数目M比内层均多6个：

第1层：M＝1

第i层：M＝6(i-1)

组件盒间的间隙面积A_axi求解：

式中：

L——六边形闭式组件盒边长/m；

S——组件盒间隙宽度/m；

径向流通面积求解：

每个盒间间隙的径向流通面积A_radial等于组件盒间隙宽度S乘以控制体轴向长度d_cv：

A_radial＝S×d_cv 公式(6)

步骤5：由于盒间结构不规则，为保证所得差分方程具有守恒特性，在此选用有限容积积分法，离散格式选用隐式迎风格式，壁面边界节点采用无滑移边界条件，控制方程组中各个方程的离散方程如下：

1)连续方程离散：

ρ_j,k-1/2A_j,k-1/2u_j,k-ρ_j,k+1/2A_j,k+1/2u_j,k+1+ρ_j+1/2,kS_j+1/2,kv_j+1,k-ρ_j-1/2,kS_j-1/2,kv_j,k＝0 公式(7)

式中：

ρ_j,k-1/2——第j通道，第k-1/2个轴向控制体的密度/kg/m³；

A_j,k-1/2——第j通道，第k-1/2个轴向控制体的轴向流通面积/m²；

S_j+1/2,k——第j+1/2通道，第k个轴向控制体的径向流通面积/m²；

u_j,k——第j通道，第k个轴向控制体的轴向流速/m/s；

v_j,k——第j通道，第k个轴向控制体的横向流速/m/s；

2)轴向动量方程离散：

公式(8)中的系数

和轴向动量源项

的表达式如下：

公式(8)中

和

分别表示一个轴向动量控制体的上、下、外和内界面位置的质量流量，其表达式为：

式中：

ΔV——轴向动量控制体体积/m³；

Δt——时间步长/s；

P_j,k——第j通道、第k轴向控制体的压力/Pa；

f_j,k——第j通道、第k轴向控制体的摩擦阻力系数；

Δz——轴向动量控制体长度/m；

——轴向动量控制体的等效水力直径/m；

k_j,k——第j通道、第k轴向控制体的局部阻力系数；

——第j通道、第k轴向控制体的轴向初始速度/m/s；

u_j,k——第j通道、第k轴向控制体的轴向速度/m/s；

v_j+1,k——第j+1通道、第k轴向控制体的径向速度/m/s；

A_j,k——第j通道、第k轴向控制体的轴向流通面积/m²；

S_j+1/2,k——第j+1/2通道、第k轴向控制体的径向流通面积/m²；

α_z——第k个主控制体和第k个轴向动量控制体重叠部分所占的比例；

β_z——第k-1个主控制体和第k个轴向动量控制体重叠部分所占的比例；

3)径向动量方程离散：

公式(9)中G_N、G_S、G_O和G_IN分别表示一个主控制体的上、下、外和内界面位置的质量流量，其表达式为：

G_N＝ρ_j,k+1/2u_j,k+1A_j,k+1；

G_S＝ρ_j,k-1/2u_j,kA_j,k；

G_O＝ρ_j+1,kv_j+1,kS_j+1,k；

G_IN＝ρ_j-1,kv_j,kS_j,k；

公式(9)中的系数

和径向动量源项

的表达式如下：

公式(9)中

和

分别表示一个径向动量控制体的上、下、外和内界面位置的质量流量，其表达式为：

式中：

Δr——径向动量控制体长度/m；

——径向动量控制体等效水力直径/m；

——第j通道、第k轴向控制体的轴向初始速度/m/s；

u_j,k+1——第j通道、第k+1轴向控制体的轴向速度/m/s；

v_j+1,k——第j+1通道、第k轴向控制体的径向速度/m/s；

A_j,k+1/2——第j通道、第k+1/2轴向控制体的轴向流通面积/m²；

S_j,k——第j通道、第k轴向控制体的径向流通面积/m²；

α_r——第j通道主控制体和第j通道动量控制体重叠部分所占的比例；

β_r——第j-1通道主控制体和第j通道动量控制体重叠部分所占的比例；

4)能量方程离散：

六边形闭式组件盒间流动的能量方程的求解是将方程离散成焓值的时间导数的表达形式，在整个堆芯及系统内耦合全场求解，离散后的能量方程为：

式中：

m_j,k——第j通道、第k轴向控制体的质量/kg；

h_j.k——第j通道、第k轴向控制体流体的比焓/J·kg^-1；

Q_j.k——组件盒间流体与组件或者是钢套屏蔽之间的换热功率/J/s；

步骤6：组件盒间流动的连续方程与动量方程采用压力修正法求解，对于能量方程则采用吉尔算法进行求解；其中组件盒间流动的连续方程与动量方程的瞬态计算的基本步骤如下：

1.读入或者传入上一时刻的轴向初始速度u⁰、径向初始速度v⁰和初始压力P⁰；

2.传入边界点的压力值；

3.根据上一时层的压力场和速度场以及压力边界求解动量方程，求解轴向速度修正值u^*、径向速度修正值v^*、压力修正值P^*；

4.求解压力修正方程，获得修正后的压力P′；

5.根据修正后的压力P′修正速度，获得修正后的轴向速度u′、修正后的径向速度v′；

6.利用改进后的速度场求解与当前速度和压力相关的物性和源项，重复步骤3-6，直至速度场收敛；

7.重复步骤1-6进行下一个时层的计算；

其中轴向速度与径向速度的修正方程如下：

式中：

u′_j,k——第j通道、第k轴向控制体修正后的轴向速度/m/s；

v′_j,k——第j通道、第k轴向控制体修正后的径向速度/m/s；

P′_j,k——第j通道、第k轴向控制体修正后的压力/Pa；

根据压力修正值P′改进的速度场应该满足连续性方程，将公式(11)、公式(12)代入连续方程得：

式中：

——第j通道、第k轴向控制体轴向速度修正值/m/s；

——第j通道、第k轴向控制体径向速度修正值/m/s；

整理得压力修正值的计算方程：

B_PP′_j,k＝B_sP′_j,k-1+B_nP′_j,k+1+B_oP′_j+1,k+B_inP′_j-1,k+S_P 公式(14)

公式(14)中的系数B_P、B_s、B_n、B_o、B_in和源项S_P的表达式如下：

B_P＝B_s+B_n+B_o+B_in

式中：

P′_j,k——第j通道、第k轴向控制体修正后的压力/Pa；

——第j通道、第k轴向控制体轴向速度修正值/m/s；

——第j通道、第k轴向控制体径向速度修正值/m/s；

其中在边界上的压力修正值均为0；

步骤7：对压力和速度的修正值进行亚松弛处理从而加快计算的收敛，其中对于压力修正值直接进行亚松弛处理：

P＝α_pP′+P^* 公式(15)

式中：

α_p——压力亚松弛因子；

P——亚松弛处理后的压力；

为了加快非线性问题的迭代收敛，对速度也进行亚松弛处理；对于速度的亚松弛处理是将亚松弛过程组织到代数方程的求解过程中，因此动量方程亚松弛后的方程式为：

式中：

α_ra——径向亚松弛因子；

α_axial——轴向亚松弛因子；

u_nb——相邻控制体轴向速度/m/s；

v_nb——相邻控制体径向速度/m/s；

——由轴向相邻控制体进入的质量流量/kg/s；

——由径向相邻控制体进入的质量流量/kg/s；

u′_j,k——上一时层计算得到的轴向速度/m/s；

v′_j,k——上一时层计算得到的径向速度/m/s；

——本时层轴向的初场；

——本时层径向的初场；

终止每个层次流场迭代求解收敛的判据为连续性方程的余量范数与出口流量的比值小于等于收敛误差ε，表示为：

式中：

S_P——各控制体上连续方程的源项；

q_m——出口质量流量/kg·s^-1；

ε——收敛误差；

其中收敛误差ε取值范围为：1.0e-3≤ε≤1.0e-5。

Images