[发明专利]基于无网格法的核反应堆燃料元件失效分析方法

权利要求书

1.一种基于无网格法的核反应堆燃料元件失效分析方法，其特征在于：步骤如下：

步骤1：对核反应堆燃料元件的初始状态进行粒子建模，得到粒子几何模型，定义1、2、3号粒子表征燃料元件芯块的液相、固液混合相、固相，定义11、12、13号粒子表征燃料包壳的液相、固液混合相、固相，定义21、22、23号粒子表征核反应堆结构材料的液相、固液混合相、固相，定义31、32、33号粒子表征核反应堆冷却剂的液相、固液混合相、固相；每种粒子根据具体表征的对象具有对应的密度、比热容、固液相线温度或熔点、沸点、弹性模量、泊松比物性参数和力学参数；由于核反应堆内不同材料会发生共晶或化学反应，不同的物质组合会形成众多不同的材料，每种材料的物性各不相同，为了正确反映堆芯材料，定义物质组分，即将粒子视为一个容器，用数种物质将其填充，根据所填充的物质组成计算对应的物性参数和力学参数；核反应堆内常见的物质有：铀U、锆Zr、氧O、氢H和不锈钢中常见的金属元素；根据燃料元件的初始状态设定粒子的初场，即设置位置、速度、温度、焓值、应力张力和应变张量；根据燃料元件的边界条件设置粒子的边界类型，即设置热源、绝热边界、固定边界和进出口边界；根据粒子类型定义Type属性，所有液相和固液混合相的Type为1，固相按连续介质进行区分，即不同的完整固体的Type数不同；粒子直径定义为l₀，其值必须小于等于实际分析对象最小尺寸的1/3，对于燃料元件，最小尺寸为燃料包壳的厚度，在燃料破裂过程的分析中，最小尺寸应取破裂形成的最小碎片的尺寸；

通过步骤1的计算建立了核反应堆燃料元件的几何结构、关键参数、初始状态和边界条件，能够还原任意复杂情况下核反应堆内燃料的真实状态；

步骤2：在预定的计算域内建立虚拟链表，预定的计算域大于步骤1所建的粒子几何模型所占的空间，虚拟链表大小由最大的粒子作用半径r_e,max决定，粒子作用半径为一个粒子能够影响其他粒子的最远距离；虚拟链表只起到定位粒子的作用，不参与实际的物理计算，通过检索粒子所在的虚拟链表能够快速地定位与其相邻的虚拟链表，并检索相邻虚拟链表内的所有粒子，从而实现快速检索与中心粒子发生作用的其他所有粒子，这些粒子被存储为中心粒子的邻点粒子，即中心粒子只和邻点粒子发生相互作用；根据维度的不同，所需检索的虚拟链表也不同，一维、二维和三维分别需要检索3、9、27个虚拟链表；

步骤3：定义权重函数来表征粒子间相互作用的大小程度，权重函数由粒子间距离决定，根据不同的相态，采用不同的权重函数和粒子作用半径；在液相和固液两相的计算中，流体的粒子作用半径r_e,liquid＝3.1l₀，采用双曲权重函数：

式中

l₀——粒子直径，m；

ω_ij,liquid——粒子i和粒子j间的流体权重函数；

r_e,liquid——流体的粒子作用半径，m；

r_ij——粒子i和粒子j的距离，m；

在固相的计算中，固体的粒子作用半径r_e,solid＝2.1l₀，采用核函数形式如下：

ω_ij,solid——粒子i和粒子j间的固体权重函数；

r_e,solid——固体的粒子作用半径，m；

在固体计算中采用纯拉格朗日体系，为了方便，定义以初始粒子位置计算得到的权重函数为初始权重函数

根据权重函数，定义粒子数密度n表征粒子疏密的程度，如公式(3)所示，在流体计算中，流体视为不可压缩或弱可压缩流体，即粒子数密度表征流体的密度变化情况，因此在后续的公式均是在将流体密度用粒子数密度n替代的假设上推导得到的；

n——粒子数密度；

ω_ij——粒子i和粒子j间的权重函数；

其中ω_ij根据不同的计算需求选择不同的权重函数，在流体计算中，采用流体权重函数，由此得到流体计算的粒子数密度n_liqud，在固体计算中，采用固体权重函数，由此得到固体计算的粒子数密度n_solid；此外，在固体计算中，固体内部的作用情况受固体初始状态决定，因此根据固体的初始状态，存储初始的固体粒子数密度为了定义归一化体积，定义参考粒子数密度n⁰，参考粒子数密度是假设粒子排布完全均匀情况下计算得到的粒子数密度，归一化体积则为由于该参数只在流体计算中出现，因此计算采用的权重函数为流体权重函数；

通过步骤3建立了粒子间相互作用的关系体系和粒子数密度，即得到了实际问题中，流体与流体、流体与固体、固体与固体之间的相互作用程度，流体和固体的密度，其中流体包括液相和固液混合相；实际问题中流体包括：冷却剂、堆芯材料熔融物、结构材料熔融物，固体包括：固相的堆芯材料和结构材料；

步骤4：核反应燃料元件失效的原因是热量无法被顺利地带出堆芯，由此出现燃料元件的热变形、破裂、熔化失效行为；针对堆芯中关键的能量传递过程构建的能量守恒方程包括传热项即流动换热和导热、辐射换热、内热源和化学热：

式中

H——粒子焓，J/kg；

t——时间，s；

ρ——粒子密度，kg/m³；

Q_inner——内热源，W/m³；

Q_chem——化学热，W/m³；

——由传热导致的焓值对时间的偏导，W/kg；

——由辐射换热导致的焓值对时间的偏导，W/kg；

对于传热项，是基于导热微分方程离散得到的粒子格式，如公式(5)所示，对于固定不动的粒子，该式表征导热过程，对于运动的物体，该式表征流动换热过程；

式中

k——粒子热导率，W/(m·K)；

d——模拟的空间维度；

n⁰——参考粒子数密度；

——粒子j的温度，K；

T_i——粒子i的温度，K；

对于辐射换热，辐射换热仅发生在表面与表面之间，因此首先检索固体表面粒子，对所有固体和固液混合态粒子进行判定，如果粒子作用半径内存在液体粒子，且粒子间距离小于1.1l₀，则判定为与液体直接接触的固体表面粒子；如果粒子数密度小于0.97倍的n₀，则判定为与空气直接接触的固体表面粒子；表面粒子与表面粒子之间存在辐射换热，应用灯源法检索存在辐射换热的粒子，即以中心粒子为中心向外放出射线，射线所触及的粒子如果为表面粒子，则该粒子与中心粒子会发生相互作用；辐射换热计算如下：

式中

ε——发射率；

σ——斯忒藩-玻耳兹曼常数；

内热源是通过设定粒子的体积释热率实现的，体积释热率大小根据实际燃料功率大小确定；

化学热的实现形式与内热源相同，也是设置粒子的体积释热率，体积释热率的大小是根据化学反应计算得到；在核反应燃料失效过程中，主要考虑锆与水的氧化还原反应，当锆包壳的温度达到1200℃以上时，锆水剧烈反应，产生大量的热量，锆水反应化学方程式如下：

两个粒子满足发生化学反应的条件为：粒子内均存在上述的两种反应物质，粒子温度达到化学反应发生的要求温度，粒子间距离足够近r_ij≤1.1l₀；发生化学反应后，粒子内的物质总类开始发生变化，通过保证两个粒子的总焓不变，计算粒子对应的温度和体积释热率，物质变化速率由化学反应速率决定，化学反应速率与温度相关；

通过能量计算获得粒子新的焓值，由焓值及物性计算对应的粒子温度：

式中

T_s,i——粒子i的固相线温度K；

T_l,i——粒子i的液相线温度K；

H_i——粒子i的焓值J/kg；

H_s,i——粒子i的固相线温度对应焓值J/kg；

H_l,i——粒子i的液相线温度对应焓值J/kg；

c_p,i——粒子i的比热容J/(kg·K)；

为了获得粒子的相态变化，定义液相率表征粒子内液相占粒子总物质的比率，如公式(9)所示，

式中

γ_i——粒子i的液相率；

当液相率为1时，粒子为纯液相，当液相率为0时，粒子为纯固相，当液相率介于两者之间时，粒子为固液混合相，且定义临界液相率γ_crit，当γ＜γ_crit时，粒子趋向固相的行为，否在趋向液相的行为，取γ_crit＝0.45；

对于冷却剂粒子，即31号粒子，当冷却剂为水且发生沸腾时，仅考虑能量守恒，当水达到沸点且吸收足够的能量时，令其消失；

通过步骤4的计算，模拟核反应堆燃料失效过程中燃料元件、堆芯结构材料和冷却剂间的传热相变过程；计算得到每个时刻下粒子的物质组成、焓值、温度和液相率，即获得了核反应堆燃料失效过程中燃料元件、堆芯结构材料和冷却剂的物质组成、焓值、相态和温度随时间的变化过程；

步骤5：在核反应堆燃料失效过程中，当堆芯材料之间发生共晶反应时，共晶反应产生的共晶产物往往具有更低的熔化温度，由此导致堆芯材料的低温液化；核反应堆失效过程中，考虑B₄C、二氧化铀、锆合金和不锈钢之间的共晶反应；共晶反应实际视为传质过程，采用菲克第二定律对其进行描述，通过粒子离散得到公式(10)：

式中

——k+1时刻的粒子i中物质x的质量kg；

——k时刻的粒子i中物质x的质量kg；

D_x——物质x的扩散系数m²/s；

Δt——时间步长s；

——k时刻的粒子j中物质x的质量kg；

由此获得了粒子内物质x的质量，通过质量计算物质x的物质摩尔分额：

式中

f_x——物质x的物质摩尔分额；

m_x——物质x的质量，kg；

M_x——物质x的摩尔质量，kg/mol；

n_x,total——粒子中物质x总的物质的量，mol；

通过物质摩尔分额，参考共晶相图，计算得到粒子的物性变化；

通过步骤5的计算，模拟核反应堆燃料失效过程中B₄C、二氧化铀、锆合金和不锈钢之间的共晶反应；计算得到每个时刻粒子内不同物质的质量，即得到核反应堆燃料、结构材料及冷却剂的物质组成随时间的变化；

步骤6：在核反应堆燃料失效过程中，存在冷却剂的流动、燃料及堆芯材料熔融物的流动，为了正确模拟该过程，考虑流体运动过程中受到的压力、粘性力、表面张力和重力，如公式(12)所示；流体计算中，将流体视为不可压缩流体，因此满足的连续方程如公式(13)所示，

式中

u——速度矢量，m/s；

P——压强，Pa；

μ——动力粘度，Pa·s；

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

g——重力加速度矢量，m/s²；

本方法中参考预估-校正思想，首先显式求解粘性项、表面张力项和重力项，计算得到粒子速度和位置的估算值，再显式求解压力项，对速度和位置进行修正；

粘度项的求解如公式(14)所示，

式中

μ_i——粒子i的动力粘度，Pa·s；

μ_j——粒子j的动力粘度，Pa·s；

u_ij——u_ij＝u_j-u_i，m/s；

u_j——粒子j的速度矢量，m/s；

u_i——粒子i的速度矢量，m/s；

r_ij——r_ij＝r_j-r_i，m；

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

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

表面张力项的求解如公式(15)所示

对于流体之间，表面张力模型修正系数计算如下：

对于固体与流体之间，表面张力模型修正系数计算如下：

式(15)至式(17)中

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

C_tension——表面张力模型修正系数；

r_min——表面张力的最小作用范围；

r_max——表面张力的最大作用范围；

m——质量，kg；

C_{tension,fluid}——流体之间的表面张力模型修正系数；

C_tension,fs——固体与流体之间的表面张力模型修正系数；

σ——表面张力系数；

θ——液体与固体之间的接触角，°；

由重力、粘性力、表面张力计算粒子的加速度，更新粒子的速度和位置，该速度和位置为估算值；

通过步骤6的计算，模拟核反应堆燃料失效过程中冷却剂、堆芯熔融物、结构材料熔融物流动过程中受到的粘性力、表面张力和重力；计算得到每个时刻粒子在粘性力、表面张力和重力作用下的速度和位置的估算值，即得到核反应堆燃料失效过程中冷却剂、堆芯熔融物、结构材料熔融物在粘性力、表面张力和重力作用下运动的趋势随时间的变化；

步骤7：在核反应堆燃料失效过程中，燃料元件、堆芯结构材料和熔融物凝固形成的壳层会由于温度升高、冷却剂冲击、固体碰撞而发生变形，因此要计算固体由于外力载荷、碰撞和温度升高而引起的热膨胀、弹性变形和塑性变形；

由热膨胀导致的应变如下：

式中：

[dε^T]——由热膨胀导致的应变增量张量；

——热膨胀应变系数；

T^ref——参考温度，K；

弹性应变和塑性应变之和为总应变减热膨胀应变：

[dε^e]+[dε^p]＝[dε]-[dε^T] 公式(19)

[dε^e]——弹性应变增量张量；

[dε^p]——塑性应变增量张量；

[dε]——总应变增量张量；

[dε^T]——由热膨胀导致的应变增量张量；

总应变由粒子的相对变形计算得到

张量和矢量之间的转化：

ε_ij——粒子i和粒子j的应变矢量；

R——旋转矩阵；

——粒子j的初始位置矢量，m；

——粒子i的初始位置矢量，m；

ε——应变矢量；

[ε]——应变张量；

——矢量与x轴方向的夹角；

——矢量与y轴方向的夹角；

——矢量与z轴方向的夹角；

假设[dε^p]＝0，计算弹性应力：

式中

——弹性应力张量第α行第β列的分量；

λ——拉梅常数的第一个参数；

——弹性应变张量第γ行第γ列的分量；

δ_αβ——克罗内克尔函数，

μ——拉梅常数的第二个参数；

——弹性应变张量第α行第β列的分量；

得到的弹性力与屈服极限进行比较，当弹性力大于屈服极限时，则发生塑性变形，否则，不存在塑性变形；

发生塑性变形的情况下，弹性应变和塑性应变的计算如下：

式中：

[dε^e]_k——k时刻的弹性应变增量张量；

[dε]_k——k时刻的总应变增量张量；

[dε^T]_k——k时刻的由热膨胀导致的应变增量张量；

[s]_k-1——k-1时刻的应力张量；

[dε^p]_k——k时刻的塑性应变增量张量；

μ——拉梅常数的第二个参数；

dλ_n——塑性变形增量系数；

由此得到的弹性应变，再采用公式(22)进行弹性应力计算；

由外力载荷和计算得到的弹性应力，计算粒子的加速度，更新固体粒子的速度和位置；

通过步骤7的计算，模拟核反应堆燃料失效过程中燃料由于外力载荷和温度升高引起的热膨胀、弹性变形和塑性变形；计算得到每个时刻粒子在热膨胀应变、弹性应力应变和塑性应变及其作用下的速度和位置，即得到核反应堆燃料失效过程中在外力载荷和温度升高作用下的变形情况随时间的变化；

步骤8：采用估算的速度和位置求解公式(12)中的压力项，压力项通过隐式求解压力泊松方程得到，压力泊松方程形式如公式(25)，

将公式(25)通过离散得到公式(26)，

式中

β——调节系数1；

ξ——调节系数2；

α——人工可压缩系数，常取10^-7；

n^*——k时刻计算得到的临时粒子数密度；

n^k——k时刻的粒子数密度；

n^k-1——k-1时刻的粒子数密度；

P_i^k+1——k+1时刻的粒子i的压强，Pa；

P_j——粒子j的压强，Pa；

P_i——粒子i的压强，Pa；

公式(26)采用不同求解器对其进行隐式求解，得到每个粒子的压强值；由压强计算得到压强梯度项，如公式(27)所示，

式中

P_i,min——粒子i的邻点粒子中压强的最小值，Pa；

由于流体在压力作用下存在运动趋势，一旦存在运动趋势，流体间可能存在速度差，由此需要考虑相对运动过程中粘性的影响，因此采用估算的速度和位置计算公式(14)，完成粘性项的修正；

通过步骤8的计算，模拟核反应堆燃料失效过程中冷却剂、堆芯熔融物、结构材料熔融物流动过程中受到的压力；计算得到每个时刻粒子在压力和粘性作用下修正得到的速度和位置，即得到核反应堆燃料失效过程中冷却剂、堆芯熔融物、结构材料熔融物在压力和粘性作用下运动的趋势随时间的变化；

步骤9：核反应堆燃料失效过程涉及复杂的多相、多组分及多物理场，需要考虑流体与固体之间的相互作用，从而实现固液边界处的计算准确；该相互作用为步骤7的固体计算提供外力载荷，为步骤8的流体压力计算提供固壁边界；

固液耦合计算是将固体表面视为流体，并参与步骤8压力泊松方程的隐式迭代求解中，并由计算得到的压力，获得流体和固体间的压力梯度，该压力梯度即为流体对固体的作用力；

对于固体的压力梯度计算，采用公式(28)，保证流体和固体之间的作用力大小相同方向相反；

式中

P_i——粒子i的压强，Pa；

P_j,min——粒子j的邻点粒子中压强的最小值，Pa；

该作用力为下一时刻步骤7的计算提供固体外力载荷；

通过步骤9的计算，模拟核反应堆燃料失效过程中冷却剂和堆芯结构、燃料的相互作用，堆芯熔融物、结构材料熔融物与堆芯结构、燃料的相互作用；计算得到每个时刻固体粒子受到流体粒子的压力作用，即得到核反应堆燃料失效过程中堆芯结构和燃料在冷却剂和熔融物压力作用下运动的趋势随时间的变化；

步骤10：输出粒子的位置、速度、温度、压力、应力场关键数据，根据设定的计算时长判定计算是否结束，若是，则结束计算，若否，则根据粒子的温度和物质组成更新粒子的物性，并返回步骤3；

通过步骤10根据堆芯内温度和化学、共晶反应过程，更新粒子的物性；

通过以上步骤对核反应堆燃料失效过程中的传热相变、共晶反应、流体运动和固体变形过程展开机理性分析。