[发明专利]一种核反应堆板状燃料熔化流固耦合无网格分析方法

权利要求书

1.一种核反应堆板状燃料熔化流固耦合无网格分析方法，其特征在于：步骤如下：

步骤1：对核反应堆板状燃料熔化流固耦合问题进行粒子建模；用不同类型的粒子代表核反应堆板状燃料的各个组成部分，0号粒子代表锆合金材质的燃料包壳和燃料基体，1号粒子代表弥散在燃料基体中的二氧化铀燃料微粒；对每一个粒子进行编号，所有的粒子均具有对应的物性参数，物性参数包括质量、密度、比热、熔沸点、温度、焓和初速度；粒子在笛卡尔坐标系内，依据核反应堆板状燃料的具体参数均匀布置，能够模拟核反应堆板状燃料的初始状态；

步骤2：使用链表法检索邻居粒子，具体方法是：将整个计算域均匀划分成正方形表格，正方形表格的边长与粒子最大作用半径相等；对于任意一个粒子而言，检索邻居粒子时只需检索包含粒子在内的表格以及包围此表格的共计9个表格，三维计算时，需要检索的就是由正方形表格围成的正方体，需要检索27个三维正方体；

步骤3：对时间步长进行稳定性分析，采用Courant-Friedrichs-Lewy条件简称CFL条件判断；CFL条件判断是差分方法保持计算稳定的必要条件，对于粒子法也适用，在移动粒子半隐式方法中采用下式进行判断

式中：

Δt——能够保持计算稳定的最小时间步长；

C——库朗常数；

l₀——粒子初始布置时相邻两粒子间的距离；

u_max——各时间步长中粒子速度最大值；

计算输入的初始条件中包含预设的时间步长大小，将输入的时间步长和使用CFL条件计算得到的能够保持计算稳定的最小时间步长进行比较，使用较大的时间步长作为当前时间步长；

步骤4：使用核函数表征粒子之间的相互作用程度，具体为在对动量守恒方程中的各项进行显式计算时，先使用核函数进行离散然后积分得到动量守恒方程中各项的具体值，使用的核函数如下所示：

式中

w(r_ij)——核函数值；

r_ij——粒子i与粒子j之间的距离；

r_e——粒子作用域半径；

步骤5：对核反应堆板状燃料熔化后的熔融物研究对象采用如下控制方程：

式中

ρ——粒子对应物质的密度kg/m³；

t——时间s；

P——压力Pa；

μ——动力粘度系数N·s/m²；

——速度矢量m/s；

——表面张力矢量N/kg；

——重力加速度m/s²；

公式(3)为连续性方程，针对的对象为不可压缩流体，公式(4)为动量守恒方程，从中看出分析对象只受到重力这一外力作用，对于任意一个粒子，通过核函数先将动量守恒方程即公式(4)中的重力项、粘性项和表面张力项离散，然后对计算域内全部粒子进行积分，显式得到粒子的估算速度和位置，其中粘性项需要使用采用如下多相粘度模型：

式中

μ——动力粘度系数N·s/m²；

μ_ij——粒子间动力粘度系数m²/s；

——速度m/s；

d——直径m；

r_e——粒子作用半径m；

n_G——使用高斯核函数计算得到的粒子数密度；

——粒子j速度m/s；

——粒子i速度m/s；

G——高斯核函数；

r_j——粒子j的位置矢量；

r_i——粒子i的位置矢量；

步骤6：使用如下的焓值相变模型计算粒子的相态：

式中

T——需要计算的粒子温度K；

T_s——需要计算的粒子对应的熔点K；

h——当前粒子的焓值J/kg；

h_s0——粒子开始熔化的焓值J/kg；

h_s1——粒子结束熔化的焓值J/kg；

c_p——粒子等压比热容J/(kg·K)；

通过步骤6的计算，能够模拟核反应堆板状燃料在最高温度超过燃料基体锆合金的熔点时的相变过程；计算得到每个粒子在不同时刻下的种类、焓值和温度，即得到锆合金基体和铀燃料微粒的相态、焓值和温度随时间的变化过程；

步骤7：得到粒子的相态后，通过非连续介质模型的分析方法离散单元法DEM对固体之间的碰撞作用进行如下的计算：

非连续介质模型的分析方法离散单元法DEM的计算对象为刚性固体，当计算对象中存在大固体时，需要先解决将离散颗粒聚集成刚性物体统一运动的问题；采用如下公式初始化大固体运动速度；

式中

——标记为ii的刚体的速度矢量m/s；

——粒子i当前的速度矢量m/s；

——标记为ii的刚体的角速度矢量rad/s；

m_i——粒子i的质量kg；

——粒子i的位置矢量；

——标记为ii的刚体的位置矢量；

I_ii——标记为ii的刚体的转动惯量；

对于任意固体单元m_ii，采用如下公式对单元受力作用后产生的运动进行分析；

式中

——标记为ii的刚体的加速度矢量；

——标记为ii的刚体的角加速度矢量；

——标记为ii的刚体的速度矢量；

——标记为ii的刚体的角速度矢量；

——标记为ii的刚体所受合力矢量；

——标记为ii的刚体所受合力矩矢量；

公式(8)和公式(9)中固体单元所受的合力与和力矩通过下式计算

式中

F_X,jj→ii——标记为jj的刚体对标记为ii的刚体施加的法向合力；

e_n——法向弹性力；

d_n——法向阻尼力；

k_n——法向弹性系数；

c_n——法向阻尼系数；

F_Y,jj→ii——标记为jj的刚体对标记为ii的刚体施加的切向合力；

e_s——切向弹性力；

d_s——切向阻尼力；

k_s——切向弹性系数；

c_s——切向阻尼系数；

这样，固体单元m_ii就移动到了一个新的位置，并产生新的接触力和接触力矩，再次计算固体单元m_ii所受的合力和合力矩重复计算流程，即得到每个固体单元及整个散粒体的运动形态，此时固体粒子的速度和位置均为估算值；

步骤7：通过核反应堆板状燃料熔化后的熔融物的不可压缩特性和动量守恒方程中的压力项联立得到如下压力泊松方程，隐式求解如下的压力泊松方程

式中

ρ——粒子对应物质密度kg/m³；

P——压力Pa；

n⁰——初始粒子数密度常数；

n^*——临时粒子数密度；

n^k——当前时层的粒子数密度；

n^k-1——上一时层的粒子数密度；

n^k-1——下一时层的粒子数密度；

β——压力源项修正系数；

γ——压力源项修正系数；

通过步骤8的计算得到动量守恒方程中的压力项，使用该压力项所得值对步骤5和步骤7中得到的粒子估算速度和位置进行修正，得到当前时间步长全部粒子的真实速度和位置，这样就完成了一个时间步长的核反应堆板状燃料熔化流固耦合的模拟；

步骤9：为了准确分析模拟对象中固体和液体之间的耦合作用，还需要对步骤5中的粘性项计算和步骤8中的压力项计算进行处理，这是因为，固体和液体之间的作用通过下式表现

式中

——流体对固体的作用合力；

——流体对固体的压力作用；

——流体对固体的粘性作用；

在步骤5的粘性项计算中，通过多相粘度模型中使用的调和粘性系数来表现固体的物性，在计算时，这个粘性项只发生在固体和液体接触时，对固体内部不进行粘性项的计算；在步骤8的压力项计算中，流体粒子只和接触的固体粒子施加压力作用，后续压力作用将在非连续介质模型的分析方法离散单元法DEM的计算中扩展至全部固体粒子；

综上，通过步骤1对核反应堆板状燃料熔化流固耦合问题中，锆合金材质的燃料包壳和燃料基体和二氧化铀燃料微粒的位置、速度和初始物性参数进行设定；通过步骤6模拟板燃料的锆合金基体在温度超过其熔点后发生相变形成熔融物，计算得到每个粒子在不同时刻下的种类、焓值和温度，得到熔融物的相态、焓值和温度随时间的变化过程；通过步骤5至步骤9模拟核反应堆板状燃料熔化流固耦合过程，模拟熔融物和未熔化的燃料微粒之间的相互作用，模拟燃料微粒之间的相互作用过程，计算得到熔融物和燃料微粒粒子的速度、位置和压力，即得到熔融物迁徙过程中的运动和压力变化过程；综合以上步骤，分析了核反应堆板状燃料熔化流固耦合问题，得到熔化流固耦合过程中熔融物的位置、速度、压力、相态、温度、焓值随时间的变化，通过以上步骤对熔融物的迁徙行为展开机理性分析。