CN113191066A - 基于无网格法的核反应堆燃料元件失效分析方法 - Google Patents

Abstract

一种基于无网格法的核反应堆燃料元件失效分析方法，步骤如下：1、建立粒子几何模型，初始化粒子参数；2、更新邻点粒子，计算权重函数和粒子数密度；3、能量守恒计算，更新粒子焓值、温度和相态；4、共晶反应计算，更新粒子物质含量和物性，更新粒子焓值、温度和相态；6、计算流体重力、粘性和表面张力，估算流体粒子速度和位置；7、固体运动计算，更新固体粒子速度和位置；8、粘性修正，更新流体粒子速度和位置；9、更新流体对固体的载荷；10、输出数据。本方法考虑核反应堆燃料元件失效过程中的所有现象；基于无网格法，能够精确捕捉界面变化，相比于网格法，避免大变形中存在网格畸变的问题；算法过程易于实现大规模并行计算。

Description

基于无网格法的核反应堆燃料元件失效分析方法

技术领域

本发明涉及核电厂严重事故堆芯燃料失效行为研究技术领域，具体涉及一种基于无网格法的核反应堆燃料元件失效分析方法。

背景技术

在核反应堆严重事故中，如果堆芯的热量无法正常导出，燃料元件和堆芯结构材料的温度会逐渐升高。随着温度的升高，核反应堆内会发生一系列反应，导致燃料元件失效。在严重事故早期阶段，温度的升高引起堆芯材料的热膨胀，局部变形产生局部应力，当应力超过对应材料的屈服应力时，堆芯材料可能发生断裂。并且，在堆芯高中子辐射条件下，堆芯材料发生脆性断裂的可能性大大增加。随着温度进一步升高，堆芯材料发生熔化、化学反应和共晶反应。当堆芯温度达到大约800℃时，堆芯内的Ag-In-Cd合金开始熔化；温度达到940℃，Fe-Zr和Ni-Zr共晶反应发生；温度达到1150℃，B₄C-Fe共晶反应发生；温度达到1200℃，Ni-Zr共晶反应，Ag-Zr反应，锆水氧化反应剧烈发生；温度达到1300℃，Fe-Zr，Al(Al₂O₃)-Zr共晶非反应发生；温度达到1450℃，不锈钢和铟科镍开始熔化；温度达到1760℃，锆-4合金开始熔化；温度达到1900℃，Al₂O₃-UO₂和Al₂O₃-ZrO₂共晶反应发生；温度达到1975℃，α-Zr(O)开始熔化；温度达到2050℃，Al₂O₃开始熔化；温度达到2350℃，B₄C开始熔化；温度达到2400℃，α-Zr(O)-UO₂和U-UO₂偏晶体产生；温度达到2600℃，U-Zr-O熔融物产生；温度达到2690℃，ZrO₂开始熔化；温度达到2850℃，UO₂开始熔化。熔融物产生后，会在堆芯通道内流动，与冷却剂和燃料元件、堆芯材料相互作用，改变堆芯状态。整个过程中存在巨大的不确定性，事故进程稍有不同，便可能存在不同的行为特性。通过以上一系列反应，严重事故进程中，核反应堆堆芯状态不断发生着变化，是一个复杂的多相态、多组分和多物理场问题。

针对如此复杂的问题，目前国内外常规的研究方法主要有实验和数值仿真。在实验方面，国外做了大量关于严重事故的实验，如QUENCH系列、Phébus系列、CORA系列等，针对核反应堆堆芯熔化和熔融物迁移过程展开了实验研究，实验中有的采用模拟材料，有的采用真实堆芯材料，获得了堆芯熔化过程中可能存在的一些现象和规律。但由于核反应堆燃料失效过程过于复杂，又由于实验无法全面考虑堆芯内可能发生的所有现象且测量难度较大，因此有必要展开数值仿真。在实验获得的现象、规律和数据的基础之上，大量学者展开了数值仿真方面的研究，开发了大量关于严重事故中燃料元件失效的数学物理模型并开发了相关程序，业内知名的严重事故分析程序有MELCOR、MAAP、ASTEC等。这类程序大多基于实验关系式进行建模，通过少量的网格，能够获得堆芯内大致的变化规律，对于工程中的核反应堆严重事故分析有重大意义。但由于严重事故分析程序不够精细，很难基于其结果展开机理分析，后续许多学者针对严重事故中燃料失效的现象展开计算流体力学(CFD)数值模拟分析，通过精细化几何建模和普适的控制方程，计算还原堆芯内可能出现的真实物理现象。传统成熟的CFD方法大多是基于网格进行几何建模，在欧拉坐标系下进行计算。但由于堆芯失效过程中，涉及到固体变形、流体运动、流体与固体之间的相互作用，采用网格法可能出现网格畸变的问题，且网格法很难考虑复杂的共晶传质和化学反应过程。无网格法由于其不依赖于网格进行建模、离散和计算，而是基于拉格朗日坐标系追踪粒子或节点运动从而还原宏观物体的运动，其在捕捉自由界面、相态变化、固体变形、流固耦合方面有着独特的优势。因此，本研究综合无网格法及核反应堆燃料失效过程的机理性分析提出一种核反应堆燃料失效分析方法。

发明内容

为了全面研究核反应堆燃料失效过程，揭示作用过程中可能存在的一些机理现象，本发明在对核反应堆燃料失效行为的机理性分析的基础，提出一种研究核反应堆燃料元件失效分析方法，该方法能够对核反应堆燃料失效过程中存在的传热相变、化学反应、共晶反应、固体变形、流体运动进行研究，获得堆芯燃料失效过程中燃料元件的变形特性、受力状态和温度场、熔融物在堆芯内的流动凝固换热过程、冷却剂的流动凝固换热过程，堆芯内复杂材料之间的化学反应和共晶反应程度，为核电厂反应堆严重事故安全特性研究提供重要依据。

为了实现上述目标，本发明采取了以下的技术方案予以实施：

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

步骤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⁰——参考粒子数密度；

T_j——粒子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——旋转矩阵；

m；

——粒子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根据堆芯内温度和化学、共晶反应过程，更新粒子的物性；通过以上步骤对核反应堆燃料失效过程中的传热相变、共晶反应、流体运动和固体变形过程展开机理性分析。

本发明方法为核反应堆燃料失效行为分析提供解决方案，为核电厂反应堆严重事故安全特性的研究提供重要依据。

和现有技术相比，本发明方法具备如下优点：

本发明的基于无网格法的核反应堆燃料元件失效分析方法，综合考虑了核反应堆燃料元件失效过程中可能存在的所有现象，包括传热相变、化学反应、共晶反应、流体运动、固体变形，能够对核反应堆燃料元件失效过程进行机理性分析。本方法能够很好地处理流体和固体计算过程中的耦合问题，在保证流体计算和固体计算精度的前提下，实现流体和固体界面的准确计算；压力计算采用全局隐式求解压力泊松方程，计算精度高，准确捕捉流体在压力作用下的运动趋势和流体对固体的外力载荷；本方法基于无网格法，能够方便地获得流体自由表面、流体-固体界面、不同组分物质界面、不同相态界面，为燃料元件失效机理分析提供重要依据；相比于网格法，无网格法有效规避了网格法中存在的网格畸变问题；本方法中仅压力采用隐式计算，其他均采用显式计算，容易实现大规模并行计算。综上，该方法能够更加全面、有效、高效地对核反应堆燃料失效过程进行分析。

附图说明

图1是本发明基于无网格法的核反应堆燃料失效分析方法的流程图。

具体实施方式

下面结合附图和具体实施例对本发明做进一步详细说明。

本发明基于无网格法的核反应堆燃料失效分析方法，如图1所示，以下以铅基核反应堆内单根燃料棒的失效过程分析为例展开介绍，步骤如下：

步骤1：对铅基核反应堆内单根燃料棒进行三维的粒子建模，所建模型长度2.2米，芯块简化为圆柱体，直径9毫米，包壳内径9.1毫米，外径10.2毫米，上下端采用不锈钢结构材料固定，整根燃料棒布置在铅铋冷却剂内；定义1、2、3号粒子表征燃料元件芯块的液相、固液混合相、固相，定义11、12、13号粒子表征燃料包壳的液相、固液混合相、固相，定义21、22、23号粒子表征不锈钢结构材料的液相、固液混合相、固相，定义31、32、33号粒子表征铅铋冷却剂的液相、固液混合相、固相；每种粒子根据具体表征的对象具有对应的密度、比热容、固液相线温度或熔点、沸点、弹性模量、泊松比等物性参数和力学参数；设置粒子直径l₀为0.05毫米；根据粒子类型定义Type属性，所有液相和固液混合相的Type为1，芯块的Type为2，包壳的Type为3，上下部分不锈钢结构的Type分别为4和5；设置燃料棒和结构材料的初始温度为500℃，铅铋温度为350℃；

通过步骤1的计算建立了铅基核反应堆内单根燃料棒的三维几何结构、关键参数、初始状态和边界条件；

步骤2：在步骤1所建立的计算域内建立虚拟链表，本算例中最大的粒子作用半径为0.15毫米，因此虚拟链表大小设置为0.15毫米；虚拟链表只起到定位粒子的作用，不参与实际的物理计算，通过检索粒子所在的虚拟链表能够快速高效地定位与其相邻的虚拟链表，并检索相邻虚拟链表内的所有粒子，从而实现快速检索与中心粒子发生作用的其他所有粒子，这些粒子被存储为中心粒子的邻点粒子，即中心粒子只和邻点粒子发生相互作用；本算例中，需要每个中心粒子需要检索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；

r_ij——粒子i和粒子j的距离，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⁰——参考粒子数密度；

T_j——粒子j的温度，K；

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

在铅基反应堆单根燃料元件失效算例中，铅铋冷却剂向环境的辐射换热对燃料元件失效过程影响较小，故此忽略公式(4)中的辐射换热项。在其他算例中，如果要考虑辐射换热，其计算如下：

式中

ε——发射率；

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

内热源是通过设定粒子的体积释热率实现的，体积释热率大小根据实际燃料功率大小确定，本算例中设置芯块的体积释热率为3×10⁸W/m³；

化学热的实现形式与内热源相同，也是设置粒子的体积释热率，体积释热率的大小是根据化学反应计算得到；在核反应燃料失效过程中，主要考虑锆与水的氧化还原反应，当锆包壳的温度达到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；

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

步骤5：在铅基反应堆单根燃料元件失效过程中，铅铋冷却剂、不锈钢结构材料、不锈钢包壳和UO₂芯块之间可能发生共晶反应，采用菲克第二定律对其进行描述，通过粒子离散得到公式(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的计算，模拟铅基反应堆单根燃料元件失效过程中铋冷却剂、不锈钢结构材料、不锈钢包壳和UO₂芯块之间的共晶反应；计算得到铅基反应堆单根燃料元件失效过程中铋冷却剂、不锈钢结构材料、不锈钢包壳和UO₂芯块的物质组成随时间的变化；

步骤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——旋转矩阵；

m；

——粒子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根据堆芯内温度和化学、共晶反应过程，更新粒子的物性；通过以上步骤对铅基反应堆单根燃料元件失效过程中的传热相变、共晶反应、流体运动和固体变形过程展开机理性分析。

Claims (1)

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根据堆芯内温度和化学、共晶反应过程，更新粒子的物性；

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

Images