[发明专利]一种板状燃料堆芯亚通道级三维热工水力分析方法

权利要求书

1.一种板状燃料堆芯亚通道级三维热工水力分析方法，其特征在于：包括如下步骤：

步骤1：建立由全堆芯通道级网格与流体域亚网格两部分组成的全堆芯亚通道级网格

(1)：通过几何建模软件建立板状燃料堆芯三维几何模型，此堆芯三维几何模型由众多相同的燃料板、燃料包壳与冷却剂通道三部分组成；

(2)：在网格划分软件中基于步骤1-(1)中建立的板状燃料堆芯三维几何模型，在水平剖面上沿各燃料板中心对称轴进行纵向网格划分，生成既包含燃料板固体域同时也包含冷却剂流体域的通道级网格；考虑到计算流体力学软件的求解性质，对于网格内部介质的有限体积法数值求解必须针对流体；因此在后续处理中通道级网格内部原燃料板固体域被简化为与冷却剂性质相同的流体域，而由该简化方式引入的人工修正在步骤3中做了详细介绍；

(3)：为了精确捕捉矩形窄缝通道内冷却剂流动与换热细节，针对通道级网格中原流体域沿纵向均匀剖分若干节点，划分的网格中左右两列网格与毗邻固体域一同被合并为新的通道级网格，而流体域其余网格即为生成亚网格，最后沿堆芯高度方向均匀划分预设数量节点；至此，完成全堆芯亚通道级网格的建立；

步骤2：对建立的全堆芯通道级网格与流体域亚网格进行网格标记

将生成的全堆芯亚通道级网格导入计算流体力学软件中，通过用户自定义函数分别对通道级网格与流体域亚网格进行网格标记，具体步骤如下：

设堆芯高度方向为z轴，垂直于燃料板长度方向为x轴，沿燃料板长度方向为y轴，假设每个标准燃料组件中燃料板的个数为NOP，则根据计算需要确定的通道级网格节点个数也为NOP，使用N表示某处通道级网格节点沿x方向所处的相对位置，1≤N≤NOP，N＝1时x＝Th_c0+Th_w+0.5Th_U，其中x为节点与燃料组件边界横向距离，Th_c0、Th_w、Th_U分别为组件边界包壳厚度、矩形窄缝厚度与铀燃料板厚度；当1＜N≤NOP时，Th_c为燃料板包壳厚度，即相邻燃料板中心距为一定值；

对于流体域亚网格，采用补集筛选法对其进行网格标记，设x方向每两块燃料板之间对冷却剂通道均匀划分亚网格节点数为n，则第N个通道网格节点处网格坐标为：

另设一变量E，令满足上述坐标的节点处E＝1，则其余节点均为亚通道级网格节点，赋值E＝0；至此，对亚通道级网格中两套不同网格的标记完成，在后续分类过程中只需引入简单的判断语句遴选变量E的取值即可；

步骤3：通过针对标记的两类网格在控制守恒方程中附加源项、修正几何参数与引入通道间湍流搅混以构建冷却剂流动换热分析模型

由于通道级网格包络了原燃料板固体域，因此在对控制体进行有限体积法积分求解时作以下几何参数修正：

控制体体积修正：

控制体侧面积修正：

其中：

V'、V₀分别为修正后与未修正的控制体体积；

S_y'、S_y0分别为修正后与未修正的控制体沿y方向表面积；

L、D分别为控制体沿垂直板排列方向的高度与宽度；

经过修正后，通道级网格内进行守恒方程的数值求解时所使用的几何参数均映射为真实的流体域参数值；

为代替全堆芯亚通道级网格划分方式对板状燃料堆芯几何固体域结构的简化带来的冷却剂流动特性的影响，需在守恒方程中引入额外的动量、能量源项，并满足以下假设：

1)冷却剂流动阻力使用由摩擦压降转换得到的动量源项表示：

式中：

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

——冷却剂速度/m·s^-1；

D_e——冷却剂通道水力直径/m，对于矩形通道通过计算得出，Δx为控制体径向厚度；

f——摩擦阻力因子；

在计算因流动摩擦引入的动量损失时，特征长度取通道级网格的水力直径；针对层流、湍流分别采用不同的阻力摩擦因子计算关系式，引入无量纲数Re，Re为雷诺数，通过下式计算得出：

其中μ为冷却剂的动力粘性系数/kg·m^-1·s^-1；

当Re2300时，冷却剂流动视为层流，摩擦阻力因子由修正后Hagen-Poiseuill公式计算得出；而当Re≥2300时，流动转为湍流，摩擦阻力因子通过Kakac经验关系式计算，具体表达式如下：

2)计算域中燃料板与包壳不发生结构形变与相对位移；

3)计算过程中冷却剂不发生过冷沸腾，始终保持单相流动；

基于上述假设，对流体控制方程中的动量方程和能量方程中添加源项：在动量方程中添加包含粘性损失、惯性损失的压降作为动量源项在能量方程中添加燃料板表面的对流换热量Q_hf作为能量源项S_E；

分析通道级网格与流体域亚网格的划分特性后发现，其毗邻界面为纯流体域，因此在交界面处考虑相邻通道P、N间湍流搅混带来的动量与能量交换，具体形式为：

为相邻通道之间湍流搅混带来的动量交换/kg·m^-1·s^-2；

为界面间平均质量流密度/kg·m^-2·s^-1，通过简单加权平均得出：

其中：

G_P、G_N分别为P、N通道内流体质量流密度/kg·m^-2·s^-1；

D_P、D_N分别为P、N通道宽度/m；

为相邻通道之间湍流搅混带来的能量交换/kJ·m^-2·s^-1；

h_P、h_N为相邻P、N控制体内流体比焓/kJ·kg^-1；

Δx_P,N为P、N两点间距离/m；

β为相邻通道间湍流搅混系数，为雷诺数与通道几何参数的函数，β通过以下关系式计算得出：

最终，考虑由通道级网格内部固体域的简化带来的动量损失、燃料板释热引入的能量源项与通道间湍流搅混作用后，给出亚通道网格与通道网格内流体域的守恒方程如下形式：

质量守恒方程：

动量守恒方程

能量守恒方程：

至此，建立了完整的冷却剂流动换热特性分析模型；

步骤4：针对原固体域设置离散节点与构建燃料板热工特性数值求解模型

分析步骤1中通道级网格的划分特性发现，其内部同时包含由燃料板组成的原固体域与冷却剂流体域，虽然在步骤3中已经讨论了固体域的简化引入的模型修正，但简化后的原固体域仍不适用于流体守恒方程；因此在通道级网格内部引入独立的板状燃料元件热工特性分析模型对其温度场进行求解，具体如下：

首先提出以下假设：

1)反应堆燃料板结构均匀分布在堆芯中；

2)单个控制体内所有燃料板具有相同的几何特性；

3)单个控制体内所有燃料板具有相同的热工特性；

根据燃料板与包壳组成的固体域几何特性沿径向划分N个节点，其中最外层两个节点分布在包壳内、外侧边界区域，其余节点从燃料板对称中心向外侧依次排布，忽略燃料板的轴向导热，即对于单个控制体，其中添加的固体域温度求解算法是一维的；对于第i个计算节点，其能量守恒方程应满足以下形式：

式中：

ρ_i为计算节点i处控制体密度/kg·m^-3；

为定压比热容kJ·kg^-1·K^-1；

Δx_i为控制体径向宽度/m；

T_i为节点i处的温度/K；

q_i-1,i与q_i+1,i分别为i-1节点与i+1节点传递到i处的面热流密度/kW·m^-2；

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、T_i+1分别为计算节点i-1、i+1处的温度/K；

K_i-1,i为计算节点i-1和i之间的传热系数/kW·m^-2·K^-1，等同于K_i,i-1；K_i+1,i为计算节点i+1和i之间的传热系数/kW·m^-2·K^-1，等同于K_i,i+1，均表示为控制体内燃料板几何参数与导热系数λ二者的函数：

式中x_i为节点i距芯块中心的径向距离，某两节点间的平均导热系数由积分平均值得出；

根据由燃料板及包壳组成的固体域以及与冷却剂交界面的传热特点，建立如下边界条件：

1)燃料板中心为对称边界，因此对于计算节点1的控制方程为绝热形式:

K_0,1＝K_i-1,i＝0；

2)在燃料芯块与包壳之间存在接触导热，即在i＝N-1处q_i,i-1＝h_s(T_i-T_i-1)，其中h_s为燃料芯块与包壳间的接触导热系数/kW·m^-2·K^-1，由下式计算得出：

式中：

ψ为带有量纲的常数，取10/m^1/2；

λ_U、λ_c分别为燃料芯块、包壳的导热系数/W·m^-1·K^-1；

P_i为表面接触压力/Pa；

Θ为燃料芯块的Meyer硬度数/N·m^-2；

δ_U、δ_c分别为在交界面处燃料芯块、包壳的表面粗糙度/m；

3)在燃料板包壳外表面，即i＝N-1处，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·m^-2·K^-1；

H为燃料板包壳表面的对流换热系数kW·m^-2·K^-1；

q_i+1,i为计算节点i+1传导到计算节点i的热量kW·m^-2；

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，系数对于每个通道级网格中固体域内N个节点都存在上述关于其自身节点温度与相邻两节点温度的三元差分方程，整理为一N+1元一次方程组，并在初始、最末方程中添加了附加的本构边界条件作为方程闭合的必要条件，对矩阵方程的化简采用高斯消去法，求解算法采用高斯迭代法，最终得到n时刻各个节点的温度值，至此，完成了针对原固体域离散节点的设置并构建了燃料板热工特性数值求解模型；

步骤5：在通道级网格内固体域与流体域交界面处通过换热关系建立流固耦合换热模型

为了进一步求解冷却剂热工水力参数，将通道网格中固体域边界，即燃料包壳最外层控制体与亚通道级网格中流体域节点联立并构建流固耦合换热模型：

对于流固边界，热流密度满足Dittus-Boleter公式：

式中：λ为流体导热系数/W·m^-1·K^-1；Pr为普朗特数；上式中包含受温度场影响而变化的流体物性参数；因此对于冷却剂流场的计算需要通过Dittus-Boleter关系式与固体域温度场进行耦合求解，具体实现步骤如下：

1)：假定亚通道级网格内初始冷却剂流场与温度场并赋值，通过检索得到对应温度下的物性参数值，通过考虑附加的动量、能量源项和湍流搅混项的冷却剂守恒方程求得该时刻冷却剂流场；

2)：将通道网格内冷却剂流场参数代入Dittus-Boleter公式求解燃料包壳外表面对流换热系数H；

3)：以燃料包壳外表面对流换热系数H为边界条件求解步骤4中固体域温度场差分方程，获得燃料板及包壳温度分布；

4)：根据燃料包壳、矩形窄缝通道内冷却剂温度、换热面积与表面对流换热系数求得表面对流换热量；

5)：将得到的表面对流换热量作为通道网格内冷却剂控制方程中附加能量源项，根据冷却剂流动换热特性分析模型求解其下一时间步亚通道级网格内流场与温度场；

6)：获得该时刻下流体物性参数，返回步骤①进行通道级网格流场参数的迭代求解，最终得到所需时间内板状燃料反应堆全堆芯热工水力参数分布。