一种板状燃料堆芯亚通道级三维热工水力分析方法

1.本发明属于板状燃料反应堆数值模拟技术领域，具体涉及到一种板状燃料堆芯亚通道级三维热工水力分析方法。


背景技术：

2.板状燃料由于其功率密度高、燃耗较深、排布紧凑、换料周期长的特点，被广泛应用于小型堆、测试堆与一体化实验堆中。其核心结构单元为由燃料芯块、包壳与冷却剂矩形窄缝流道组成的闭式通道。在进行板状燃料反应堆正常或事故工况下的瞬态安全分析时，燃料板最高温度、包壳表面温度、冷却剂流场等热工水力参数极为重要。但由于燃料组件排布的复杂性与矩形窄缝通道数量的繁多，实现针对板状燃料反应堆全堆芯的数值模拟目前仍具有较大挑战性。
3.现阶段国内外针对板状燃料反应堆的热工水力分析方法主要包括以下三类：
4.1、一维系统分析程序——如relap5、mitars、trace，此类程序能够构建核电厂宏观回路与各类热力构件模型，同时满足在较短时间内计算核电厂系统各处稳态参数分布与瞬态响应的要求。但该类程序使用离散的一维节点描述堆内各类系统级构件，如蒸汽发生器、稳压器、阀门、主泵，在程序中均采用单个节点表示，无法刻画构件内物理过程与参数分布。对于复杂的堆芯结构，仅使用单节点的集总值显然无法满足精细尺度下堆芯安全分析的需求。
5.2、子通道分析程序——如cobra、cosine，此类程序将堆芯划分为相邻的并联流道，并采用各类描述热工水力参数间关系的经验表达式计算子通道内冷却剂流动特性值，相比于一维系统程序能够较精确地反映堆芯各流道内部的流动换热细节，但此类程序对于瞬变工况模拟能力仍具有较大不足，且目前现有的子通道程序大多基于压水堆棒束通道进行开发，在板状燃料反应堆中应用较少。
6.3、计算流体力学(cfd)仿真软件——如fluent、cfx、star-ccm ，此类软件要求针对堆芯几何模型进行微尺度的网格划分，通过求解控制体质量、动量、能量守恒方程获得流场参数与温度分布。cfd仿真精度较高，结果可视化能力高，在采用恰当计算模型与方法的条件下能够捕捉冷却剂精细流动与传热现象。但同 cfd方法要求在几何建模软件中构建精细的堆芯结构，所需求的建模工作量较大；且计算速度缓慢，计算资源耗费量大。


技术实现要素：

7.本发明的目的在于提供一种板状燃料堆芯亚通道级三维热工水力研究方法，该方法通过全堆芯亚通道级网格划分，固体域节点温度场差分计算与流固耦合模型的建立，能够在保证计算精度的条件下，显著降低计算资源与计算时间。
8.为了达到上述目的，本发明采用如下技术方案：
9.一种板状燃料堆芯亚通道级三维热工水力研究方法，包括如下步骤：
10.步骤1：建立由全堆芯通道级网格与流体域亚网格两部分组成的亚通道级网格
11.(1)：通过几何建模软件solidworks建立板状燃料堆芯三维几何模型，此堆芯三维几何模型由众多相同的燃料芯块板、燃料包壳与冷却剂通道三部分组成；
12.(2)：在网格划分软件中基于步骤1-(1)中建立的板状燃料堆芯三维几何模型，在水平剖面上沿各燃料板中心对称轴进行纵向网格划分，生成既包含燃料板固体域同时也包含冷却剂流体域的通道级网格。需要特别说明的是，考虑到计算流体力学(cfd)软件的求解性质，对于网格内部介质的有限体积法数值求解必须针对流体；因此在后续处理中通道级网格内部原燃料板固体域被简化为与冷却剂性质相同的流体域，而由该简化方式引入的人工修正在步骤3中做了详细介绍；
13.(3)：为了精确捕捉矩形窄缝通道内冷却剂流动与换热细节，针对通道级网格中原流体域沿纵向均匀剖分若干节点，划分的网格中左右两列网格与毗邻固体域一同被合并为新的通道级网格，而流体域其余网格即为生成的亚网格，最后沿堆芯高度方向均匀划分预设数量节点；至此，完成全堆芯亚通道级网格的建立，通过该方式生成的网格数量远小于常规cfd模拟所需要的数目，为建模与网格划分过程节省了可观的工作量；
14.步骤2：对建立的全堆芯通道级网格与流体域亚网格进行网格标记
15.将生成的全堆芯亚通道级网格导入计算流体力学软件中，通过用户自定义函数分别对通道级网格与流体域亚网格进行网格标记，具体步骤如下：
16.设堆芯高度方向为z轴，垂直于燃料板长度方向为x轴，沿燃料板长度方向为y轴，假设每个标准燃料组件中燃料板的个数为nop，则根据计算需要确定的通道级网格节点个数也为nop，使用n表示某处通道级网格节点沿x方向所处的相对位置，1≤n≤nop，n＝1时x＝th
c0
thw 0.5thu，其中x为节点与燃料组件边界横向距离，th
c0
、thw、thu分别为组件边界包壳厚度、矩形窄缝厚度与铀燃料板厚度；当1＜n≤nop时，deltan＝δx(n,n-1)＝thw 2thc thu， thc为燃料板包壳厚度，即相邻燃料板中心距为一定值；
17.对于流体域亚网格，采用补集筛选法对其进行网格标记，设x方向每两块燃料板之间对冷却剂通道均匀划分亚网格节点数为n，则第n个通道网格节点处网格坐标为：
[0018][0019]
另设一变量e，令满足上述坐标的节点处e＝1，则其余节点均为亚通道级网格节点，赋值e＝0；至此，对亚通道级网格中两套不同网格的标记完成，在后续分类过程中只需引入简单的判断语句遴选变量e的取值即可；
[0020]
步骤3：通过针对标记的两类网格在控制守恒方程中附加源项、修正几何参数与引入通道间湍流搅混以构建冷却剂流动换热分析模型
[0021]
由于通道级网格包络了原燃料板固体域，因此在对控制体进行有限体积法积分求解时作以下几何参数修正：
[0022]
控制体体积修正：
[0023]
控制体侧面积修正：
[0024]
其中：
[0025]
v'、v0分别为修正后与未修正的控制体体积；
[0026]
sy'、s
y0
分别为修正后与未修正的控制体沿y方向表面积；
[0027]
l、d分别为控制体沿垂直板排列方向的高度与宽度；
[0028]
经过修正后，通道级网格内进行守恒方程的数值求解时所使用的几何参数均映射为真实的流体域参数值；
[0029]
为代替全堆芯亚通道级网格划分方式对板状燃料堆芯几何固体域结构的简化带来的冷却剂流动特性的影响，需在守恒方程中引入额外的动量、能量源项，并满足以下假设：
[0030]
1)冷却剂流动阻力使用由摩擦压降转换得到的动量源项表示：
[0031][0032]
式中：
[0033]
ρ——冷却剂密度/kg
·
m-3
；
[0034]
——冷却剂速度/m
·
s-1
；
[0035]de
——冷却剂通道水力直径/m，对于矩形通道通过计算得出，δx 为控制体径向厚度；
[0036]
f——摩擦阻力因子；
[0037]
在计算因流动摩擦引入的动量损失时，特征长度取通道级网格的水力直径；针对层流、湍流分别采用不同的阻力摩擦因子计算关系式，引入无量纲数re，re
[0038]
为雷诺数，通过下式计算得出：其中μ为冷却剂的动力粘性系数/kg
·
m-1
·
s-1
；
[0039]
当re《2300时，冷却剂流动视为层流，摩擦阻力因子由修正后hagen-poiseuill 公式计算得出；而当re≥2300时，流动转为湍流，摩擦阻力因子通过kakac经验关系式计算，具体表达式如下：
[0040][0041]
2)计算域中燃料板与包壳不发生结构形变与相对位移；
[0042]
3)计算过程中冷却剂不发生过冷沸腾，始终保持单相流动；
[0043]
基于上述假设，对流体控制方程中的动量方程和能量方程中添加源项：在动量方程中添加包含粘性损失、惯性损失的压降作为动量源项在能量方程中添加燃料板表面
的对流换热量q
hf
作为能量源项se；
[0044]
分析通道级网格与流体域亚网格的划分特性后发现，其毗邻界面为纯流体域，因此在交界面处考虑相邻通道p、n间湍流搅混带来的动量与能量交换，具体形式为：
[0045][0046][0047]
为相邻通道之间湍流搅混带来的动量交换/kg
·
m-1
·
s-2
；
[0048]
为界面间平均质量流密度/kg
·
m-2
·
s-1
，通过简单加权平均得出：
[0049][0050]
其中：
[0051]gp
、gn分别为p、n通道内流体质量流密度/kg
·
m-2
·
s-1
；
[0052]dp
、dn分别为p、n通道宽度/m；
[0053]
为相邻通道之间湍流搅混带来的能量交换/kj
·
m-2
·
s-1
；
[0054]hp
、hn为相邻p、n控制体内流体比焓/kj
·
kg-1
；
[0055]
δx
p,n
为p、n两点间距离/m；
[0056]
β为相邻通道间湍流搅混系数，为雷诺数与通道几何参数的函数，β通过以下关系式计算得出：
[0057][0058]
最终，考虑由通道级网格内部固体域的简化带来的动量损失、燃料板释热引入的能量源项与通道间湍流搅混作用后，给出亚通道网格与通道网格内流体域的守恒方程如下形式：
[0059]
质量守恒方程：
[0060][0061]
动量守恒方程
[0062]
能量守恒方程：
[0063]
[0064]
至此，建立了完整的冷却剂流动换热特性分析模型；
[0065]
步骤4：针对原固体域设置离散节点与构建燃料板热工特性数值求解模型
[0066]
分析步骤1中通道级网格的划分特性发现，其内部同时包含由燃料板组成的原固体域与冷却剂流体域，虽然在步骤3中已经讨论了固体域的简化引入的模型修正，但简化后的原固体域仍不适用于流体守恒方程；因此在通道级网格内部引入独立的板状燃料元件热工特性分析模型对其温度场进行求解，具体如下：
[0067]
首先提出以下假设：
[0068]
1)反应堆燃料板结构均匀分布在堆芯中；
[0069]
2)单个控制体内所有燃料板具有相同的几何特性；
[0070]
3)单个控制体内所有燃料板具有相同的热工特性；
[0071]
根据燃料板与包壳组成的固体域几何特性沿径向划分n个节点，其中最外层两个节点分布在包壳内、外侧边界区域，其余节点从燃料板对称中心向外侧依次排布，忽略燃料板的轴向导热，即对于单个控制体，其中添加的固体域温度求解算法是一维的；对于第i个计算节点，其能量守恒方程应满足以下形式：
[0072][0073]
式中：
[0074]
ρi为计算节点i处控制体密度/kg
·
m-3
；
[0075]cpi
为定压比热容kj
·
kg-1
·
k-1
；
[0076]
δxi为控制体径向宽度/m；
[0077]
ti为节点i处的温度/k；
[0078]qi-1,i
与q
i 1,i
分别为i-1节点与i 1节点传递到i处的面热流密度/kw
·
m-2
； qi为节点i处控制体内体积释热率/kw
·
m-3
。
[0079]
其中，两侧相邻计算节点传递给中心节点的面热流密度通过傅立叶导热定律求得：
[0080]qi-1,i
＝k
i-1,i
(t
i-1-ti)
[0081]qi 1,i
＝k
i 1,i
(t
i 1-ti)
[0082]
式中：
[0083]
t
i-1
、t
i 1
分别为计算节点i-1、i 1处的温度/k；
[0084]ki-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
，均表示为控制体内燃料板几何参数与导热系数λ二者的函数：
[0085][0086][0087]
式中xi为节点i距芯块中心的径向距离，某两节点间的平均导热系数由积分平均
值得出；
[0088]
根据由燃料板及包壳组成的固体域以及与冷却剂交界面的传热特点，建立如下边界条件：
[0089]
1)燃料板中心为对称边界，因此对于计算节点1的控制方程为绝热形式: k
0,1
＝k
i-1,i
＝0；
[0090]
2)在燃料芯块与包壳之间存在接触导热，即在i＝n-1处q
i,i-1
＝hs(t
i-t
i-1
)，其中hs为燃料芯块与包壳间的接触导热系数/kw
·
m-2
·
k-1
，由下式计算得出：
[0091][0092]
式中：
[0093]
ψ为带有量纲的常数，取10/m
1/2
；
[0094]
λu、λc分别为燃料芯块、包壳的导热系数/w
·
m-1
·
k-1
；
[0095]
pi为表面接触压力/pa；
[0096]
θ为燃料芯块的meyer硬度数/n
·
m-2
；
[0097]
δu、δc分别为在交界面处燃料芯块、包壳的表面粗糙度/m；
[0098]
3)在燃料板包壳外表面，即i＝n-1处，k
i 1,i
＝h，t
i 1
＝tf，则有 q
i 1,i
＝h(t
f-ti)；
[0099]
其中：k
i 1,i
为计算节点i 1和i之间的传热系数kw
·
m-2
·
k-1
；
[0100]
h为燃料板包壳表面的对流换热系数kw
·
m-2
·
k-1
；
[0101]qi 1,i
为计算节点i 1传导到计算节点i的热量kw
·
m-2
；
[0102]
tf为矩形窄缝内冷却剂主流温度/k；
[0103]
t
i 1
为计算节点i 1处的温度/k；
[0104]
ti为计算节点i处的温度/k；
[0105]
为了实现温度场的程序求解，将能量守恒方程转化为差分方程形式：式中：t
in
为计算节点i处在第n个时间步下的温度/k，t
in-1
为计算节点i处在第n-1个时间步下的温度/k，t
in 1
为计算节点i处在第n 1个时间步下的温度/k；
[0106]
将差分方程变形为固体域温度场差分方程：
[0107][0108]
式中：系数ai＝-k
i-1,i
，系数系数ci＝-k
i 1,i
，系数对于每个通道级网格中固体域内n个节点都存在上述关于其自身节点温度与相邻两节点温度的三元差分方程，整理为一n 1元一次方程组，并在初始、最末方程中添加了附加的本构边界条件作为方程闭合的必要条件，对矩阵方程的化简采用高斯消去法，求解算法采用高斯迭代法，最终得到n时刻各个节点的温度值，至此，完
成了针对原固体域离散节点的设置并构建了燃料板热工特性数值求解模型。
[0109]
步骤5：在通道级网格内固体域与流体域交界面处通过换热关系建立流固耦合换热模型
[0110]
为了进一步求解冷却剂热工水力参数，将通道网格中固体域边界，即燃料包壳最外层控制体与亚通道级网格中流体域节点联立并构建流固耦合换热模型：
[0111]
对于流固边界，热流密度满足dittus-boleter公式：
[0112][0113]
式中：λ为流体导热系数/w
·
m-1
·
k-1
；pr为普朗特数；上式中包含受温度场影响而变化的流体物性参数；因此对于冷却剂流场的计算需要通过 dittus-boleter关系式与固体域温度场进行耦合求解，具体实现步骤如下：
[0114]
1)：假定亚通道级网格内初始冷却剂流场与温度场并赋值，通过检索得到对应温度下的物性参数值，通过考虑附加的动量、能量源项和湍流搅混项的冷却剂守恒方程求得该时刻冷却剂流场；
[0115]
2)：将通道网格内冷却剂流场参数代入dittus-boleter公式求解燃料包壳外表面对流换热系数h；
[0116]
3)：以燃料包壳外表面对流换热系数h为边界条件求解步骤4中固体域温度场差分方程，获得燃料板及包壳温度分布；
[0117]
4)：根据燃料包壳、矩形窄缝通道内冷却剂温度、换热面积与表面对流换热系数求得表面对流换热量；
[0118]
5)：将得到的表面对流换热量作为通道网格内冷却剂控制方程中附加能量源项，根据冷却剂流动换热特性分析模型求解其下一时间步亚通道级网格内流场与温度场；
[0119]
6)：获得该时刻下流体物性参数，返回步骤
①
进行通道级网格流场参数的迭代求解，最终得到所需时间内板状燃料反应堆全堆芯热工水力参数分布。
[0120]
和现有技术相比较，本发明具有以下特点：
[0121]
1、堆芯控制体划分方式简单，计算网格数量少；
[0122]
2、在通道级网格的基础上，采用的进一步亚通道级网格划分能够精确捕捉矩形窄缝通道内流动换热细节；
[0123]
3、克服了一维系统分析程序无法获得板状燃料堆芯固体域准确的温度分布、对功率分布不均匀瞬态计算过程响应模拟较差的问题；
[0124]
4、针对燃料板温度场求解设计的算法模型普适性强，可在简单参数调节后应用于其他堆型的固体域传热计算；
[0125]
5、耗费计算资源少，计算速度快且计算精度高；
[0126]
6、构建的流固耦合换热模型可在获得冷却剂热工水力参数的同时计算得到燃料板、包壳的温度场，对板状燃料反应堆安全分析具有重要意义。
[0127]
本发明方法已通过程序完成实践验证，在该方法下针对板状燃料全堆芯的数值仿真能够准确地获得燃料板、包壳温度场以及冷却剂流场、温度场等热工水力特性参数，具有精确捕捉冷却剂通道内流体域流动换热细节的能力，且同时具有计算资源耗费量少，并行
能力高，计算速度快的特点。
附图说明
[0128]
图1为相邻燃料板间通道级/亚通道级网格划分示意图。
[0129]
图2为燃料板与包壳组成固体域内温度场差分算法节点选取示意图。
[0130]
图3为流固耦合模型流程图。
[0131]
图4为本发明板状燃料堆芯亚通道级三维热工水力分析方法系统框图。
具体实施方式
[0132]
为了更深层次地介绍本发明技术路线与实施意义，以下根据图4所示系统框图为主要流程，结合附图对本发明作进一步的详细描述。
[0133]
本发明提出了一种板状燃料堆芯亚通道级三维热工水力研究方法，包括以下步骤：
[0134]
步骤1：建立由全堆芯通道级网格与流体域亚网格两部分组成的全堆芯亚通道级网格
[0135]
(1)：通过几何建模软件建立板状燃料堆芯三维几何模型，此堆芯三维几何模型由众多相同的燃料板、燃料包壳与冷却剂通道三部分组成；常见的板状燃料反应堆燃料组件分别为板状燃料反应堆标准组件盒与控制组件盒结构示意图，主要区别在于控制组件盒两侧布置了硼材料制成的中子吸收体，而标准组件盒仅由铝制包壳包裹的铀燃料芯块板组成；
[0136]
(2)：在网格划分软件中基于步骤1-(1)中建立的板状燃料堆芯三维几何模型，在水平剖面上沿各燃料板中心对称轴进行纵向网格划分，生成既包含燃料板固体域同时也包含冷却剂流体域的通道级网格；考虑到计算流体力学软件的求解性质，对于网格内部介质的有限体积法数值求解必须针对流体；因此在后续处理中通道级网格内部原燃料板固体域被简化为与冷却剂性质相同的流体域，而由该简化方式引入的人工修正在步骤3中做了详细介绍；
[0137]
(3)：为了精确捕捉矩形窄缝通道内冷却剂流动与换热细节，针对通道级网格中原流体域沿纵向均匀剖分若干节点，划分的网格中左右两列网格与毗邻固体域一同被合并为新的通道级网格，而流体域其余网格即为生成亚网格，最后沿堆芯高度方向均匀划分预设数量节点；至此，完成全堆芯亚通道级网格的建立，建立的网格示意如图1所示，通过该方式生成的网格数量远小于常规cfd模拟所需要的数目，为建模与网格划分过程节省了可观的工作量；
[0138]
步骤2：对建立的全堆芯通道级网格与流体域亚网格进行网格标记
[0139]
将生成的全堆芯亚通道级网格导入计算流体力学软件中，通过用户自定义函数分别对通道级网格与流体域亚网格进行网格标记，具体步骤如下：
[0140]
设堆芯高度方向为z轴，垂直于燃料板长度方向为x轴，沿燃料板长度方向为y轴，假设每个标准燃料组件中燃料板的个数为nop，则根据计算需要确定的通道级网格节点个数也为nop，使用n表示某处通道级网格节点沿x方向所处的相对位置，1≤n≤nop，n＝1时x＝th
c0
thw 0.5thu，其中x为节点与燃料组件边界横向距离，th
c0
、thw、thu分别为组件边界包
壳厚度、矩形窄缝厚度与铀燃料板厚度；当1＜n≤nop时，deltan＝δx(n,n-1)＝thw 2thc thu， thc为燃料板包壳厚度，即相邻燃料板中心距为一定值；
[0141]
对于流体域亚网格，采用补集筛选法对其进行网格标记，设x方向每两块燃料板之间对冷却剂通道均匀划分亚网格节点数为n，则第n个通道网格节点处网格坐标为：
[0142][0143]
另设一变量e，令满足上述坐标的节点处e＝1，则其余节点均为亚通道级网格节点，赋值e＝0；至此，对亚通道级网格中两套不同网格的标记完成，在后续分类过程中只需引入简单的判断语句遴选变量e的取值即可；
[0144]
步骤3：通过针对标记的两类网格在控制守恒方程中附加源项、修正几何参数与引入通道间湍流搅混以构建冷却剂流动换热分析模型
[0145]
由于通道级网格包络了原燃料板固体域，因此在对控制体进行有限体积法积分求解时作以下几何参数修正：
[0146]
控制体体积修正：
[0147]
控制体侧面积修正：
[0148]
其中：
[0149]
v'、v0分别为修正后与未修正的控制体体积；
[0150]
sy'、s
y0
分别为修正后与未修正的控制体沿y方向表面积；
[0151]
l、d分别为控制体沿垂直板排列方向的高度与宽度；
[0152]
经过修正后，通道级网格内进行守恒方程的数值求解时所使用的几何参数均映射为真实的流体域参数值；
[0153]
为代替全堆芯亚通道级网格划分方式对板状燃料堆芯几何固体域结构的简化带来的冷却剂流动特性的影响，需在守恒方程中引入额外的动量、能量源项，并满足以下假设：
[0154]
1)冷却剂流动阻力使用由摩擦压降转换得到的动量源项表示：
[0155][0156]
式中：
[0157]
ρ——冷却剂密度/kg
·
m-3
；
[0158]
——冷却剂速度/m
·
s-1
；
[0159]de
——冷却剂通道水力直径/m，对于矩形通道通过计算得出，δx 为控制体径向厚度；
[0160]
f——摩擦阻力因子；
[0161]
在计算因流动摩擦引入的动量损失时，特征长度取通道级网格的水力直径；针对层流、湍流分别采用不同的阻力摩擦因子计算关系式，引入无量纲数re，re
[0162]
为雷诺数，通过下式计算得出：其中μ为冷却剂的动力粘性系数/kg
·
m-1
·
s-1
；
[0163]
当re《2300时，冷却剂流动视为层流，摩擦阻力因子由修正后hagen-poiseuill 公式计算得出；而当re≥2300时，流动转为湍流，摩擦阻力因子通过kakac经验关系式计算，具体表达式如下：
[0164][0165]
2)计算域中燃料板与包壳不发生结构形变与相对位移；
[0166]
3)计算过程中冷却剂不发生过冷沸腾，始终保持单相流动；
[0167]
基于上述假设，对流体控制方程中的动量方程和能量方程中添加源项：在动量方程中添加包含粘性损失、惯性损失的压降作为动量源项在能量方程中添加燃料板表面的对流换热量q
hf
作为能量源项se；
[0168]
分析通道级网格与流体域亚网格的划分特性后发现，其毗邻界面为纯流体域，因此在交界面处考虑相邻通道p、n间湍流搅混带来的动量与能量交换，具体形式为：
[0169][0170][0171]
为相邻通道之间湍流搅混带来的动量交换/kg
·
m-1
·
s-2
；
[0172]
为界面间平均质量流密度/kg
·
m-2
·
s-1
，通过简单加权平均得出：
[0173][0174]
其中：
[0175]gp
、gn分别为p、n通道内流体质量流密度/kg
·
m-2
·
s-1
；
[0176]dp
、dn分别为p、n通道宽度/m；
[0177]
为相邻通道之间湍流搅混带来的能量交换/kj
·
m-2
·
s-1
；
[0178]hp
、hn为相邻p、n控制体内流体比焓/kj
·
kg-1
；
[0179]
δx
p,n
为p、n两点间距离/m；
[0180]
β为相邻通道间湍流搅混系数，为雷诺数与通道几何参数的函数，β通过以下关系式计算得出：
[0181][0182]
最终，考虑由通道级网格内部固体域的简化带来的动量损失、燃料板释热引入的能量源项与通道间湍流搅混作用后，给出亚通道网格与通道网格内流体域的守恒方程如下形式：
[0183]
质量守恒方程：
[0184][0185]
动量守恒方程
[0186]
能量守恒方程：
[0187][0188]
至此，建立了完整的冷却剂流动换热特性分析模型；
[0189]
步骤4：针对原固体域设置离散节点与构建燃料板热工特性数值求解模型
[0190]
分析步骤1中通道级网格的划分特性发现，其内部同时包含由燃料板组成的原固体域与冷却剂流体域，虽然在步骤3中已经讨论了固体域的简化引入的模型修正，但简化后的原固体域仍不适用于流体守恒方程；因此在通道级网格内部引入独立的板状燃料元件热工特性分析模型对其温度场进行求解，具体如下：
[0191]
首先提出以下假设：
[0192]
1)反应堆燃料板结构均匀分布在堆芯中；
[0193]
2)单个控制体内所有燃料板具有相同的几何特性；
[0194]
3)单个控制体内所有燃料板具有相同的热工特性；
[0195]
根据燃料板与包壳组成的固体域几何特性沿径向划分n个节点，其中最外层两个节点分布在包壳内、外侧边界区域，其余节点从燃料板对称中心向外侧依次排布,最终的节点划分示意图如图2；忽略燃料板的轴向导热，即对于单个控制体，其中添加的固体域温度求解算法是一维的；对于第i个计算节点，其能量守恒方程应满足以下形式：
[0196][0197]
式中：
[0198]
ρi为计算节点i处控制体密度/kg
·
m-3
；
[0199]cpi
为定压比热容kj
·
kg-1
·
k-1
；
[0200]
δxi为控制体径向宽度/m；
[0201]
ti为节点i处的温度/k；
[0202]qi-1,i
与q
i 1,i
分别为i-1节点与i 1节点传递到i处的面热流密度/kw
·
m-2
； qi为节点i处控制体内体积释热率/kw
·
m-3
。
[0203]
其中，两侧相邻计算节点传递给中心节点的面热流密度通过傅立叶导热定律求得：
[0204]qi-1,i
＝k
i-1,i
(t
i-1-ti)
[0205]qi 1,i
＝k
i 1,i
(t
i 1-ti)
[0206]
式中：
[0207]
t
i-1
、t
i 1
分别为计算节点i-1、i 1处的温度/k；
[0208]ki-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
，均表示为控制体内燃料板几何参数与导热系数λ二者的函数：
[0209][0210][0211]
式中xi为节点i距芯块中心的径向距离，某两节点间的平均导热系数由积分平均值得出；
[0212]
根据由燃料板及包壳组成的固体域以及与冷却剂交界面的传热特点，建立如下边界条件：
[0213]
1)燃料板中心为对称边界，因此对于计算节点1的控制方程为绝热形式: k
0,1
＝k
i-1,i
＝0；
[0214]
2)在燃料芯块与包壳之间存在接触导热，即在i＝n-1处q
i,i-1
＝hs(t
i-t
i-1
)，其中hs为燃料芯块与包壳间的接触导热系数/kw
·
m-2
·
k-1
，由下式计算得出：
[0215][0216]
式中：
[0217]
ψ为带有量纲的常数，取10/m
1/2
；
[0218]
λu、λc分别为燃料芯块、包壳的导热系数/w
·
m-1
·
k-1
；
[0219]
pi为表面接触压力/pa；
[0220]
θ为燃料芯块的meyer硬度数/n
·
m-2
；
[0221]
δu、δc分别为在交界面处燃料芯块、包壳的表面粗糙度/m；
[0222]
3)在燃料板包壳外表面，即i＝n-1处，k
i 1,i
＝h，t
i 1
＝tf，则有q
i 1,i
＝h(t
f-ti)；
[0223]
其中：k
i 1,i
为计算节点i 1和i之间的传热系数kw
·
m-2
·
k-1
；
[0224]
h为燃料板包壳表面的对流换热系数kw
·
m-2
·
k-1
；
[0225]qi 1,i
为计算节点i 1传导到计算节点i的热量kw
·
m-2
；
[0226]
tf为矩形窄缝内冷却剂主流温度/k；
[0227]
t
i 1
为计算节点i 1处的温度/k；
[0228]
ti为计算节点i处的温度/k；
[0229]
为了实现温度场的程序求解，将能量守恒方程转化为差分方程形式：式中：t
in
为计算节点i处在第n个时间步下的温度/k，t
in-1
为计算节点i处在第n-1个时间步下的温度/k，t
in 1
为计算节点i处在第n 1个时间步下的温度/k；
[0230]
将差分方程变形为固体域温度场差分方程：
[0231][0232]
式中：系数ai＝-k
i-1,i
，系数系数ci＝-k
i 1,i
，系数对于每个通道级网格中固体域内n个节点都存在上述关于其自身节点温度与相邻两节点温度的三元差分方程，整理为一n 1元一次方程组，并在初始、最末方程中添加了附加的本构边界条件作为方程闭合的必要条件。对于初始节点与最末节点，分别不考虑前、后计算节点，即：
[0233]
i＝1时，a1＝0
[0234]
i＝n时，cn＝0
[0235]
将边界条件代入差分方程中整理可得：
[0236][0237]
左侧n
×
n系数为中心三对角矩阵，对矩阵方程的化简采用高斯消去法得到对角矩阵：
[0238][0239]
化简后对角矩阵求解算法采用高斯迭代法，最终得到n时刻各个节点的温度值；至此，完成了针对原固体域离散节点的设置并构建了燃料板热工特性数值求解模型；
[0240]
步骤5：在通道级网格内固体域与流体域交界面处通过换热关系建立流固耦合换热模型
[0241]
为了进一步求解冷却剂热工水力参数，将通道网格中固体域边界，即燃料包壳最外层控制体与亚通道级网格中流体域节点联立并构建流固耦合换热模型：
[0242]
对于流固边界，热流密度满足dittus-boleter公式：
[0243][0244]
式中：λ为流体导热系数/w
·
m-1
·
k-1
；pr为普朗特数；上式中包含受温度场影响而变化的流体物性参数；因此对于冷却剂流场的计算需要通过 dittus-boleter关系式与固体域温度场进行耦合求解，求解流程如附图3所示，具体实现步骤如下：
[0245]
1)：假定亚通道级网格内初始冷却剂流场与温度场并赋值，通过检索得到对应温度下的物性参数值，通过考虑附加的动量、能量源项和湍流搅混项的冷却剂守恒方程求得该时刻冷却剂流场；
[0246]
2)：将通道网格内冷却剂流场参数代入dittus-boleter公式求解燃料包壳外表面对流换热系数h；
[0247]
3)：以燃料包壳外表面对流换热系数h为边界条件求解步骤4中固体域温度场差分方程，获得燃料板及包壳温度分布；
[0248]
4)：根据燃料包壳、矩形窄缝通道内冷却剂温度、换热面积与表面对流换热系数求得表面对流换热量；
[0249]
5)：将得到的表面对流换热量作为通道网格内冷却剂控制方程中附加能量源项，根据冷却剂流动换热特性分析模型求解其下一时间步亚通道级网格内流场与温度场；
[0250]
6)：获得该时刻下流体物性参数，返回步骤
①
进行通道级网格流场参数的迭代求解，最终得到所需时间内板状燃料反应堆全堆芯热工水力参数分布。