基于粒子法的熔融物与混凝土相互作用分析方法

摘要

一种基于粒子法的熔融物与混凝土相互作用分析方法，主要步骤如下：1、粒子建模，设定粒子初始布置和参数；2、建立背景网格并按节点划分，检索邻居粒子；3、计算粒子焓值、温度和相态；4、计算共晶和化学反应，更新粒子物质组分和物性参数，同时计算化学热，更新粒子焓值、温度和相态；6、计算气体运动，更新粒子类型、速度和位置；7、显示计算动量方程中的粘性项、表面张力项和重力项，估算粒子速度和位置；8、使用估算的粒子位置显示计算动量方程中的压力项并修正粒子速度和位置；9、输出计算结果。本方法考虑熔融物与混凝土相互作用的所有现象；基于粒子法，能够精确捕捉自由液面、方便建模及精确处理相变；采用显示压力模型有利于进行大规模计算。

说明书

技术领域

本发明涉及核电厂严重事故堆芯熔融物与混凝土相互作用研究技术领域，具体涉及一种基于粒子法的熔融物与混凝土相互作用分析方法。

背景技术

当核电厂压水堆发生严重事故时，堆芯如果得不到充分的冷却可能发生熔化并向下迁移，落入压力容器下封头的堆芯熔融可能会将下封头熔穿，随后，堆芯熔融物会进入安全壳并与安全壳内的混凝土发生长期的相互作用，这个过程涉及众多的化学和物理变化。高温的熔融物会不断加热混凝土，使其温度不断升高并发生熔化和化学分解。混凝土分解会产生水蒸气和其他不可凝气体，这些气体会影响熔池的流动换热行为并可能导致安全壳的超压失效。混凝土分解还会在混凝土表面形成熔渣层，影响其与堆芯熔融物的换热过程。与此同时，堆芯熔融物的温度会不断降低，并在混凝土-熔融物界面首先形成壳层，阻碍熔融物和混凝土进一步的相互作用。此外，熔融物和混凝土的相互作用还受到熔融物的初始温度、质量和成分、熔融物下落速率、注水的时间、混凝土的组成成分等因素影响。因此，熔融物与混凝土相互作用过程存在大量的不确定性，是核反应堆严重事故研究领域的热点和难点问题，至今其机理仍未被完全研究透彻，对其的研究有助于严重事故源项及安全壳完整性的分析，对于核电厂严重事故安全分析具有重大意义。

关于熔融物与混凝土相互作用的研究，国内外已经展开了一些研究，包括实验研究和数值模拟研究。对于实验研究，由于不同的研究侧重点和实验条件的限制，各实验采用不同比例的实验装置、不同的熔融物模拟物、不同的混凝土、不同的注入方式以及不同的加热方式，并且有的实验还考虑了熔融物衰变热、水的注入等的影响，主要探究了熔融物对混凝土的烧蚀过程、不可凝气体的产生。由实验研究可知，熔融物和混凝土的种类会很大程度影响熔融物对混凝土的烧蚀过程，不同的熔融物和混凝土会引起不同的烧蚀模式，如硅酸混凝土呈现烧蚀的各向异性，而石灰石混凝土呈现烧蚀的各向同性，再如金属熔融物和氧化物熔融物会出现明显的分层现象，形成多层熔池的结构特征。对于数值模拟研究，国内外基本都以集总参数法进行自编程分析，很少使用基于网格法的CFD软件对具体的相互作用过程进行分析，这是因为熔融物与混凝土相互作用过程过于复杂，气泡产生、熔融物流动行为、混凝土烧蚀界面、物质相变过程、熔池相界面等均难使用基于欧拉方法的网格法进行模拟。而对于基于拉格朗日方法的粒子法，在处理自由表面、物质流动、物质相变、气泡行为上有着独特的优势，能够很好的重现熔融物流动、熔融物和混凝土相变、熔融物-混凝土界面变化的过程。目前，有少数学者采用移动粒子法对熔融物和混凝土相互作用过程展开了模拟，但并未全面考虑熔融物和混凝土相互作用过程中的全部机理性现象，并且模拟过程中作了大量假设，特别是忽略了熔融物和混凝土相互作用过程中的化学分解、不可凝气体的产生。因此，本研究综合粒子法及熔融物与混凝土相互作用过程的机理性分析提出一种熔融物与混凝土相互作用分析方法。

发明内容

为了全面研究熔融物与混凝土相互作用过程，揭示作用过程中可能存在的一些机理现象，本发明在对熔融物与混凝土相互作用的机理性分析的基础上结合粒子法、基础控制方程及相关的机理性化学物理模型，提出一种研究熔融物与混凝土相互作用分析方法，该方法能够对熔融物与混凝土相互作用过程中存在的多种物质流动、传热、传质、相变、化学反应、气体运动进行研究，获得熔融物与混凝土相互作用过程中熔池组成及形态的变化、混凝土-熔融物界面的变化、不可凝气体的运动、化学反应程度，为核电厂反应堆严重事故安全特性研究提供重要依据。

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

一种基于粒子法的熔融物与混凝土相互作用分析方法，步骤如下：

步骤1：对熔融池与混凝土的初始状态进行粒子建模，用不同种类的粒子代表不同的物质，采用1、2、3号粒子模拟熔融物粒子的液相、固液混合相、固相，采用4、5、6号粒子模拟混凝土的液相、固液混合相、固相，每种粒子根据所表示的物质具有对应的质量、密度、比热、熔沸点、温度、焓信息；由于熔融池与混凝土相互作用过程中物质组成复杂，存在众多互可溶或互不可溶的物质组成，对于互相可溶的物质，对于粒子i添加物质x物质摩尔份额属性f_i,x，以判别单个粒子的物质组成成分，而对于互不可溶的物质，两两组分不会存在于同一个粒子中；定义某个粒子i的相关参数为Parameter_i，Parameter为参数量，则粒子i的质量、密度、比热、熔点即固相线温度和液相线温度、沸点、温度、焓、压力、速度矢量、位置矢量分别为m_i、ρ_i、Cp_i、Ts_i即Ts0_i和Ts1_i、Tb_i、T_i、h_i、P_i、粒子直径定义为l₀；由此得到熔融物和混凝土的初始位置分布和相关参数；

步骤2：在建立的粒子模型区域内建立背景网格，背景网格均匀布置，网格为正方形，边长为最大的粒子作用距离r_e,max；每个粒子均存在于一个网格上或网格面组成的正方体内；对于每一个粒子i，其坐标位置为(x_i，y_i，z_i)，其只能与包括它所处的网格体在内的27个网格体内的粒子发生作用；针对每个粒子i，检索其周围27个网格体内的所有粒子，当粒子i与粒子j的距离l_ij小于最大的粒子作用距离r_e,max，将粒子j列为粒子i的邻居粒子组，邻居粒子组内的粒子数设为邻居粒子数n_nei，即该检索过程会得到每个粒子的邻居粒子组粒子i与粒子j的距离l_ij由公式(1)计算：

步骤3：对步骤2建立的背景网格按节点数进行划分，保证每个节点计算域内的粒子数基本相同；划分得到每个节点计算域内所占网格体的上下限(X_up，Y_up，Z_up)、(X_down，Y_down，Z_down)；定义处于边界处的网格体为边界网格体，边界网格体向相邻节点计算域传递粒子信息，实现进程间的并行计算；

步骤4：定义一个加权函数来衡量粒子受附近粒子的作用程度，采用指数多项式型核函数，如公式(2)所示，

式中r_e为粒子作用范围，r为粒子间距，w为核函数；

进一步定义粒子数密度，如公式(3)，用以衡量粒子疏密程度，

n_i＝∑_j1iw(r)公式(3)

式中：n_i为i粒子的粒子数密度，j为i粒子周围邻居粒子符号，i为i粒子符号；

步骤5：能量守恒方程如公式(4)所示，

式中

h——粒子焓值J/kg；

t——时间s；

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

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

T——粒子温度K；

Q_radiation——辐射热源W/m³；

Q_out——外热源W/m³；

Q_chem——化学热W/m³；

对于辐射换热，首先检索表面粒子，令粒子数密度小于0.97倍的n₀为表面粒子，其中n₀为初始的粒子数密度；只对表面粒子进行辐射换热计算，采用斯忒藩-玻尔兹曼定律，如公式(5)所示，

式中

Q_radiation——辐射热源W/m³；

ε——发射率；

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

T——粒子温度K；

T_env——环境温度K；

l₀——粒子直径m；

对于传热过程，采用导热微分方程的离散格式，如公式(6)所示，

式中

——下一个时刻的粒子i的焓值J/kg；

——当前时刻的粒子i的焓值J/kg；

d——维度；

n₀——初始的粒子数密度；

ρ_i——i粒子密度kg/m³；

——当前时刻的粒子j温度K；

T_i^k——当前时刻的粒子i温度K；

Δt——时间步长s；

——j粒子对i粒子的核函数值，表达形式如公式(2)；——j粒子位置矢量；

——i粒子位置矢量；

——粒子i和粒子j热导率的调和平均值W/(m·K)；

k_i——粒子i热导率W/(m·K)；

k_j——粒子j热导率W/(m·K)；

Q＝Q_out+Q_chem——热量源项W/m³；

Q_out——外热源W/m³；

Q_chem——化学热W/m³；

通过焓值确定粒子的温度，对于单质如公式(7)所示，对于混合物如公式(8)所示

式中

T——粒子温度K；

T_s——粒子熔点K；

h——粒子焓值J/kg；

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

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

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

式中

T——粒子温度K；

T_s——粒子固相线温度K；

T_l——粒子液相线温度K；

h——粒子焓值J/kg；

h_s——粒子固相线温度对应焓值J/kg；

h_l——粒子液相线温度对应焓值J/kg；

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

由焓值定义固相率来表示物质所处的相态，如公式(9)所示，

γ——粒子固相率；

h——粒子焓值J/kg；

h_s——粒子固相线温度对应焓值J/kg；

h_l——粒子液相线温度对应焓值J/kg；

对于单质的固相率计算，只要将h_s和h_l分别用h_s0和h_s1替代；

当γ＝0时，粒子为液态；当γ＝1时，粒子为固态；当0<γ<1时，粒子为固液混合态；

通过步骤5的计算，模拟熔融物和混凝土相互作用过程中熔融池内的液相粒子的传热、熔融池与固体混凝土的接触界面的传热、熔融物和混凝土的相变过程；计算得到每个粒子在不同时刻下的种类、焓值和温度，即得到熔融物和混凝土的相态、焓值和温度随时间的变化过程；

步骤6：共晶反应计算，反应堆堆芯材料中存在的UO₂、锆合金、不锈钢两两之间可能发生共晶反应，因此定义堆芯熔融物组分包含UO₂、Zr、不锈钢份额，这些物质之间存在物质传递，传质过程满足菲克第二定律，如公式(10)所示，

式中

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

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

D——扩散系数m²/s；

Δt——时间步长s；

d——维度；

n₀——初始的粒子数密度；

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

——j粒子对i粒子的核函数值，表达形式如公式(2)；

——j粒子位置矢量；

——i粒子位置矢量；

由此获得每个粒子中物质x的物质摩尔份额f_x＝m_x/M_x/n_total，其中m_x为粒子中物质x的质量，M_x为物质x的摩尔质量，n_total为粒子中总的物质的量；通过伪二元共晶相图或三元相图即能够判定粒子的物性参数变化；

通过步骤6的计算，得到熔融物粒子中UO₂、Zr、不锈钢物质摩尔份额的变化，即熔融物的物质分布情况；并通过物质分布，得到熔融物的物性参数变化；

步骤7：化学反应计算，在熔融物和混凝土相互作用过程中会存在大量的化学方应，主要包括混凝土的分解反应和熔融物的氧化反应；

混凝土的分解反应主要有：

400℃下氢氧化钙脱水：Ca(OH)₂+1340kJ/kg→CaO+H₂O(g)

600℃下碳酸钙分解：CaCO₃+1637kJ/kg→CaO+CO₂(g)

1462℃下Fe₂O₃转变：6Fe₂O₃+480kJ/kg→4Fe₃O₄+O₂(g)

熔融物的氧化反应主要有：

Zr+2H₂O→ZrO₂+2H₂+6.3MJ/kg

Zr+2CO₂→ZrO₂+2CO+5.7MJ/kg

Fe+H₂O+3.0kJ/kg→FeO+H₂

Fe+CO₂+480kJ/kg→FeO+CO

Zr+SiO₂→ZrO₂+Si+1.6MJ/kg温度<1870℃

Zr+2SiO₂+4.7MJ/kg→ZrO₂+2SiO(g)温度>1870℃

Si+2H₂O→SiO₂+2H₂+15MJ/kg

Si+2CO₂→SiO₂+2CO+14MJ/kg

基于如上化学方程式，当两个粒子间相互接触，且温度达到反应温度或具有足够的内热源进行化学方应，则两个粒子进行物质转换生成新的粒子物质原子份额，期间通过控制粒子焓值的形式保证前后物质的能量守恒；在1870℃以下，对于Zr和SiO₂，当两个粒子中含有Zr和SiO₂，粒子相接触后发生反应，两种粒子中Zr和SiO₂的物质份额会转化为ZrO₂和Si物质份额，并以内热源的形式释放1.6MJ/kg的热量，此处需引入一个假设的释热速率，设为Δt时间完成全部化学热的释放；转变后的粒子根据温度更新对应粒子的物性，按初始粒子的温度比例计算转换后的粒子的温度，如公式(11)和公式(12)所示，

式中T表示粒子温度，h表示粒子焓值，下标Zr表示锆，下标SiO₂表示二氧化硅，下标ZrO₂表示二氧化锆，下标Si表示硅；

结合温度和焓值之间的转换关系，如公式(7)或公式(8)，即能够计算得到反应后的粒子温度和焓值；通过这种转化形式，在保证粒子能量守恒的前提下，尽可能避免粒子物性变化引起的温度突变而导致的温度计算震荡；

以上对化学反应的处理方法的前提是整个粒子均完全发生化学方应，也就是要求粒子足够小，认为粒子直径小于等于0.1mm时能够满足该前提条件；

通过步骤7的计算，得到了每个粒子内的物质份额变化，即得到了混凝土的分解反应过程和熔融物的氧化反应过程中物质的变化情况；

步骤8：气泡生长过程计算，在化学反应过程中，可能会产生不可凝气体，不可凝气体的存在会起到搅浑融池和局部增压；化学反应产生气体的过程是瞬间完成的，但是气体的膨胀过程是持续的，气体的生长过程，主要为两个过程，一是气泡直径不断变大的过程，二是气泡直径基本不变，继续运动的过程；

对于第一个过程，在气体生成处生成或转变出气体粒子，该粒子温度、压力取周围粒子平均值，定义一个气体生长时间Δt_bubble，该生长时间小于移动粒子法计算的时间步长，该粒子的计算半径随气体生长时间不断增大，增大速率取决于气泡的生长速率，直到生成的气体作用半径所处的球空间的体积乘以对应温度下气体的密度等于生成的气体质量，此时，该粒子计算半径等于气泡大小，随后，在计算半径区域内填充气体粒子；该过程避免了气体体积的突变而引起的压力震荡的过程；

对于第二个过程，基于多相流模型，如公式(13)、公式(14)、公式(15)所示，

公式(13)到公式(15)中

——j粒子对i粒子的高斯核函数值，表示形式如公式(16)；

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

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

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

——j粒子位置矢量；

——i粒子位置矢量；

d——维度；

r_e——粒子作用半径m；

——基于高斯核函数的初始粒子数密度；

μ_i——i粒子的动力粘度N·s/m²；

μ_j——j粒子的动力粘度N·s/m²；

P_i——i粒子的压力Pa；

P_j——j粒子的压力Pa；

——j粒子对i粒子的核函数值，表达形式如公式(2)；ρ_i——i粒子的密度kg/m³；

ρ_j——j粒子的密度kg/m³；

——当前时刻i粒子的速度m/s；

n^k+1——下一时刻的粒子数密度；

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

P_i^k+1——下一时刻i粒子的压力；

Δt——时间步长s；

β——人工调节系数，取值在0.01到0.05；

α——人工可压缩系数，取值在10-⁹到10-⁷；

σ——表面张力系数；

k_i——中心粒子处的局部等高线曲率；

C——颜色函数，表示形式如公式(17)；

算子<>为光滑算子，计算表达式如公式(18)；

式中

r_e——粒子作用半径m；

式中

Parameter_i——i粒子的相关参数；

w——核函数，表达形式如公式(2)；

r——粒子间距m；

V——i粒子以粒子作用半径的球空间内部；

通过以上计算，得到气体在液相中的流动行为；

步骤9：连续性方程如公式(19)

式中

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

t——时间s；

对于液相，将其视为不可压缩流体，只在计算压力过程中添加相关的弱压缩系数；步骤10：动量方程如公式(20)

式中

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

t——时间s；

P——粒子压力Pa；

μ——粒子动力粘度N·s/m²；

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

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

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

对于压力计算，采用显示压力模型进行计算，如公式(21)

式中

P_i^k+1——i粒子下一个时刻的压力Pa；

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

c——人造音速m/s，一般取为最大粒子速度的3倍；

B——显示压力计算模型系数，一般为7；

n_i——i粒子的粒子数密度；

n₀——初始的粒子数密度；

对于最大作用半径范围内不含气体粒子的粒子，采用压力梯度模型如公式(22)计算；对于附近含气体粒子的粒子，采用公式(14)计算；

式中

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

P——粒子压力Pa；

d——维度；

n₀——初始的粒子数密度；

——j粒子位置矢量；

——i粒子位置矢量；

ρ_i——i粒子密度kg/m³；

ρ_j——j粒子密度kg/m³；

P_j——j粒子压力Pa；

P_i,min——i粒子所有邻居粒子中压力的最小值Pa；

P_i^k+1——下一时刻i粒子的压力Pa；

α——人工可压缩系数，取值在10^-9到10^-7；

Δt——时间步长s；

——j粒子对i粒子的核函数值，表达形式如公式(2)；

对于粘度计算，采用粘度变化模型，如公式(23)

μ＝μ₀exp(2.5Aγ)公式(23)

式中

μ——粒子动力粘度N·s/m²；

μ₀——初始动力粘度N·s/m²；

A——粘度变化系数，对于Zr和UO₂设为3.0；对于不锈钢和熔融混凝土设为2.0；

γ——粒子固相率；

粘度项计算采用公式(13)；

对于固-液，液-液界面之间的表面张力采用基于自由能模型的表面张力模型计算，如公式(24)

式中

F——自由能系数，对于液-液界面如公式(25)，对于固-液界面如公式(26)；

l_ij——i粒子和j粒子的距离；

l_min——i粒子与周围的粒子的最小距离，采用1.5l₀；

r_e——粒子作用半径；

式中

F_ff——液-液界面的自由能系数

F_fs——固-液界面的自由能系数；

σ——粒子表面张力系数；

θ——粒子接触角°；

步骤11：速度和位置估算，计算完动量方程公式(20)中的重力项、粘性项、表面张力项后，对速度和位置进行估算，如公式(27)和公式(28)；

式中

——i粒子的估算速度矢量m/s；

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

μ——粒子动力粘度N·s/m²；

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

ρ_i——i粒子密度kg/m³；

——i粒子的表面张力矢量N/kg；

——粒子重力加速度矢量m/s²；

——i粒子的估算位置矢量m；

——当前时刻i粒子的位置矢量m；

Δt——时间步长s；

步骤12：速度和位置修正，计算完动量方程公式(20)中的压力项，对速度和位置进行修正，如公式(29)和公式(30)；

式中

——i粒子的估算速度矢量m/s；

——下一时刻i粒子的速度矢量m/s；

——i粒子的估算位置矢量m；

——下一时刻i粒子的位置矢量m；

ρ_i——i粒子密度kg/m³；

Δt——时间步长s；

P——粒子压力Pa；

通过步骤9至步骤12，得到了每个粒子的速度和位置，即得到了所有熔融物和混凝土的速度和位置，由此模拟得到了熔融物在和混凝土相互作用过程中的运动过程；

综上，通过步骤1对熔融物和混凝土的相互作用过程中熔融物和混凝土的位置、速度和初始物性参数进行设定；通过步骤5模拟熔融物和混凝土相互作用过程中熔融池内的液相粒子的传热、熔融池与固体混凝土的接触界面的传热、熔融物和混凝土的相变过程，计算得到每个粒子在不同时刻下的种类、焓值和温度，即得到熔融物和混凝土的相态、焓值和温度随时间的变化过程；通过步骤6模拟熔融物中UO₂、Zr、不锈钢两两之间的共晶反应，得到熔融物粒子中UO₂、Zr、不锈钢的物质摩尔份额的变化，即熔融物内的物质分布情况，并通过物质分布，更新熔融物的物性参数；通过步骤7模拟混凝土的分解反应过程和熔融物的氧化反应过程，计算得到每个粒子内的物质份额变化，即得到混凝土的分解反应过程和熔融物的氧化反应过程中物质的消失和生成过程；通过步骤8模拟气体的生成过程，计算得到气体粒子的出现和位置、速度的变化，即熔融池内由于混凝土分解产生的气体的生成和生长过程；通过步骤9至步骤12模拟熔融物和混凝土的运动过程，计算得到熔融物和混凝土粒子的速度、位置和压力，即得到熔融物和混凝土在相互作用过程中的运动和压力变化过程；综合以上步骤，模拟了熔融物和混凝土的相互作用过程，得到了作用过程中熔融物、混凝土和不可凝气体的位置、速度、压力、相态、温度、焓值、物质原子份额随时间的变化，通过以上数据对熔融物和混凝土相互作用过程中的传热相变、共晶反应、化学方应、气体产生和生长过程展开机理性分析。

本发明方法为堆芯熔融物和混凝土相互作用过程分析提供解决方案，为核电厂反应堆严重事故安全特性的研究提供重要依据。

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

本发明的基于粒子法的熔融物与混凝土相互作用分析方法，综合考虑了熔融物与混凝土相互作用过程中可能出现的所有现象，包括化学反应、共晶反应、气体产生、流体运动、传热相变，能够对熔融物与混凝土相互作用过程进行机理性分析。并且该方法基于粒子法，具备精确捕捉自由液面、方便建模以及精确处理相变问题的优点。除此之外，在压力计算中采用了显式压力求解替代隐式压力求解过程，在大规模计算中能够有效提高该方法的计算速率。综上，该方法能够更加全面、有效、高效地对熔融物与混凝土相互作用进行分析。

附图说明

图1是本发明基于粒子法的熔融物与混凝土相互作用分析方法的流程图。

具体实施方式

本发明基于粒子法的熔融物与混凝土相互作用分析方法，如图1所示，步骤如下：

步骤1：对熔融池与混凝土的初始状态进行粒子建模，用不同种类的粒子代表不同的物质，采用1、2、3号粒子模拟熔融物粒子的液相、固液混合相、固相，采用4、5、6号粒子模拟混凝土的液相、固液混合相、固相，每种粒子根据所表示的物质具有对应的质量、密度、比热、熔沸点、温度、焓信息；由于熔融池与混凝土相互作用过程中物质组成复杂，存在众多互可溶或互不可溶的物质组成，对于互相可溶的物质，对于粒子i添加物质x物质摩尔份额属性f_i,x，以判别单个粒子的物质组成成分，而对于互不可溶的物质，两两组分不会存在于同一个粒子中；定义某个粒子i的相关参数为Parameter_i，Parameter为参数量，则粒子i的质量、密度、比热、熔点即固相线温度和液相线温度、沸点、温度、焓、压力、速度矢量、位置矢量分别为m_i、ρ_i、Cp_i、Ts_i即Ts0_i和Ts1_i、Tb_i、T_i、h_i、P_i、粒子直径定义为l₀；由此得到熔融物和混凝土的初始位置分布和相关参数；

步骤2：在建立的粒子模型区域内建立背景网格，背景网格均匀布置，网格为正方形，边长为最大的粒子作用距离r_e,max；每个粒子均存在于一个网格上或网格面组成的正方体内；对于每一个粒子i，其坐标位置为(x_i，y_i，z_i)，其只能与包括它所处的网格体在内的27个网格体内的粒子发生作用；针对每个粒子i，检索其周围27个网格体内的所有粒子，当粒子i与粒子j的距离l_ij小于最大的粒子作用距离r_e,max，将粒子j列为粒子i的邻居粒子组，邻居粒子组内的粒子数设为邻居粒子数n_nei，即该检索过程会得到每个粒子的邻居粒子组粒子i与粒子j的距离l_ij由公式(1)计算：

步骤3：对步骤2建立的背景网格按节点数进行划分，保证每个节点计算域内的粒子数基本相同；划分得到每个节点计算域内所占网格体的上下限(X_up，Y_up，Z_up)、(X_down，Y_down，Z_down)；定义处于边界处的网格体为边界网格体，边界网格体向相邻节点计算域传递粒子信息，实现进程间的并行计算；

步骤4：定义一个加权函数来衡量粒子受附近粒子的作用程度，采用指数多项式型核函数，如公式(2)所示，

式中r_e为粒子作用范围，r为粒子间距，w为核函数；

进一步定义粒子数密度，如公式(3)，用以衡量粒子疏密程度，

n_i＝∑_j≠iw(r)公式(3)

式中：n_i为i粒子的粒子数密度，j为i粒子周围邻居粒子符号，i为i粒子符号；

步骤5：能量守恒方程如公式(4)所示，

式中

h——粒子焓值J/kg；

t——时间s；

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

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

T——粒子温度K；

Q_radiation——辐射热源W/m³；

Q_out——外热源W/m³；

Q_chem——化学热W/m³；

对于辐射换热，首先检索表面粒子，令粒子数密度小于0.97倍的n₀为表面粒子，其中n₀为初始的粒子数密度；只对表面粒子进行辐射换热计算，采用斯忒藩-玻尔兹曼定律，如公式(5)所示，

式中

Q_radiation——辐射热源W/m³；

ε——发射率；

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

T——粒子温度K；

T_env——环境温度K；

l₀——粒子直径m；

对于传热过程，采用导热微分方程的离散格式，如公式(6)所示，

式中

——下一个时刻的粒子i的焓值J/kg；

——当前时刻的粒子i的焓值J/kg；

d——维度；

n₀——初始的粒子数密度；

ρ_i——i粒子密度kg/m³；

——当前时刻的粒子j温度K；

T_i^k——当前时刻的粒子i温度K；

Δt——时间步长s；

——j粒子对i粒子的核函数值，表达形式如公式(2)；

——j粒子位置矢量；

——i粒子位置矢量；

——粒子i和粒子j热导率的调和平均值W/(m·K)；

k_i——粒子i热导率W/(m·K)；

k_j——粒子j热导率W/(m·K)；

Q＝Q_out+Q_chem——热量源项W/m³；

Q_out——外热源W/m³；

Q_chem——化学热W/m³；

通过焓值确定粒子的温度，对于单质如公式(7)所示，对于混合物如公式(8)所示

式中

T——粒子温度K；

T_s——粒子熔点K；

h——粒子焓值J/kg；

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

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

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

式中

T——粒子温度K；

T_s——粒子固相线温度K；

T_l——粒子液相线温度K；

h——粒子焓值J/kg；

h_s——粒子固相线温度对应焓值J/kg；

h_l——粒子液相线温度对应焓值J/kg；

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

由焓值定义固相率来表示物质所处的相态，如公式(9)所示，

γ——粒子固相率；

h——粒子焓值J/kg；

h_s——粒子固相线温度对应焓值J/kg；

h_l——粒子液相线温度对应焓值J/kg；

对于单质的固相率计算，只要将h_s和h_l分别用h_s0和h_s1替代。

当γ＝0时，粒子为液态；当γ＝1时，粒子为固态；当0<γ<1时，粒子为固液混合态；

通过步骤5的计算，模拟熔融物和混凝土相互作用过程中熔融池内的液相粒子的传热、熔融池与固体混凝土的接触界面的传热、熔融物和混凝土的相变过程；计算得到每个粒子在不同时刻下的种类、焓值和温度，即得到熔融物和混凝土的相态、焓值和温度随时间的变化过程；

步骤6：共晶反应计算，反应堆堆芯材料中存在的UO₂、锆合金、不锈钢两两之间可能发生共晶反应，因此定义堆芯熔融物组分包含UO₂、Zr、不锈钢份额，这些物质之间存在物质传递，传质过程满足菲克第二定律，如公式(10)所示，

式中

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

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

D——扩散系数m²/s；

Δt——时间步长s；

d——维度；

n₀——初始的粒子数密度；

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

——j粒子对i粒子的核函数值，表达形式如公式(2)；

——j粒子位置矢量；

——i粒子位置矢量；

由此获得每个粒子中物质x的物质摩尔份额f_x＝m_x/M_x/n_total，其中m_x为粒子中物质x的质量，M_x为物质x的摩尔质量，n_total为粒子中总的物质的量；通过伪二元共晶相图或三元相图即能够判定粒子的物性参数变化；

通过步骤6的计算，得到熔融物粒子中UO₂、Zr、不锈钢物质摩尔份额的变化，即熔融物的物质分布情况；并通过物质分布，得到熔融物的物性参数变化；

步骤7：化学反应计算，在熔融物和混凝土相互作用过程中会存在大量的化学方应，主要包括混凝土的分解反应和熔融物的氧化反应；

混凝土的分解反应主要有：

400℃下氢氧化钙脱水：Ca(OH)₂+1340kJ/kg→CaO+H₂O(g)

600℃下碳酸钙分解：CaCO₃+1637kJ/kg→CaO+CO₂(g)

1462℃下Fe₂O₃转变：6Fe₂O₃+480kJ/kg→4Fe₃O₄+O₂(g)

熔融物的氧化反应主要有：

Zr+2H₂O→ZrO₂+2H₂+6.3MJ/kg

Zr+2CO₂→ZrO₂+2CO+5.7MJ/kg

Fe+H₂O+3.0kJ/kg→FeO+H₂

Fe+CO₂+480kJ/kg→FeO+CO

Zr+SiO₂→ZrO₂+Si+1.6MJ/kg温度<1870℃

Zr+2SiO₂+4.7MJ/kg→ZrO₂+2SiO(g)温度>1870℃

Si+2H₂O→SiO₂+2H₂+15MJ/kg

Si+2CO₂→SiO₂+2CO+14MJ/kg

基于如上化学方程式，当两个粒子间相互接触，且温度达到反应温度或具有足够的内热源进行化学方应，则两个粒子进行物质转换生成新的粒子物质原子份额，期间通过控制粒子焓值的形式保证前后物质的能量守恒；在1870℃以下，对于Zr和SiO₂，当两个粒子中含有Zr和SiO₂，粒子相接触后发生反应，两种粒子中Zr和SiO₂的物质份额会转化为ZrO₂和Si物质份额，并以内热源的形式释放1.6MJ/kg的热量，此处需引入一个假设的释热速率，设为Δt时间完成全部化学热的释放；转变后的粒子根据温度更新对应粒子的物性，按初始粒子的温度比例计算转换后的粒子的温度，如公式(11)和公式(12)所示，

式中T表示粒子温度，h表示粒子焓值，下标Zr表示锆，下标SiO₂表示二氧化硅，下标ZrO₂表示二氧化锆，下标Si表示硅；

结合温度和焓值之间的转换关系，如公式(7)或公式(8)，即能够计算得到反应后的粒子温度和焓值；通过这种转化形式，在保证粒子能量守恒的前提下，尽可能避免粒子物性变化引起的温度突变而导致的温度计算震荡；

以上对化学反应的处理方法的前提是整个粒子均完全发生化学方应，也就是要求粒子足够小，认为粒子直径小于等于0.1mm时能够满足该前提条件；

通过步骤7的计算，得到了每个粒子内的物质份额变化，即得到了混凝土的分解反应过程和熔融物的氧化反应过程中物质的变化情况；

步骤8：气泡生长过程计算，在化学反应过程中，可能会产生不可凝气体，不可凝气体的存在会起到搅浑融池和局部增压；化学反应产生气体的过程是瞬间完成的，但是气体的膨胀过程是持续的，气体的生长过程，主要为两个过程，一是气泡直径不断变大的过程，二是气泡直径基本不变，继续运动的过程；

对于第一个过程，在气体生成处生成或转变出气体粒子，该粒子温度、压力取周围粒子平均值，定义一个气体生长时间Δt_bubble，该生长时间小于移动粒子法计算的时间步长，该粒子的计算半径随气体生长时间不断增大，增大速率取决于气泡的生长速率，直到生成的气体作用半径所处的球空间的体积乘以对应温度下气体的密度等于生成的气体质量，此时，该粒子计算半径等于气泡大小，随后，在计算半径区域内填充气体粒子；该过程避免了气体体积的突变而引起的压力震荡的过程；

对于第二个过程，基于多相流模型，如公式(13)、公式(14)、公式(15)所示，

公式(13)到公式(15)中

——j粒子对i粒子的高斯核函数值，表示形式如公式(16)；

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

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

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

——j粒子位置矢量；

——i粒子位置矢量；

d——维度；

r_e——粒子作用半径m；

——基于高斯核函数的初始粒子数密度；

μ_i——i粒子的动力粘度N·s/m²；

μ_j——j粒子的动力粘度N·s/m²；

P_i——i粒子的压力Pa；

P_j——j粒子的压力Pa；

——j粒子对i粒子的核函数值，表达形式如公式(2)；

ρ_i——i粒子的密度kg/m³；

ρ_j——j粒子的密度kg/m³；

——当前时刻i粒子的速度m/s；

n^k+1——下一时刻的粒子数密度；

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

P_i^k+1——下一时刻i粒子的压力；

Δt——时间步长s；

β——人工调节系数，取值在0.01到0.05；

α——人工可压缩系数，取值在10-⁹到10-⁷；

σ——表面张力系数；

k_i——中心粒子处的局部等高线曲率；

C——颜色函数，表示形式如公式(17)；

算子<>为光滑算子，计算表达式如公式(18)；

式中

r_e——粒子作用半径m；

式中

Parameter_i——i粒子的相关参数；

w——核函数，表达形式如公式(2)；

r——粒子间距m；

V——i粒子以粒子作用半径的球空间内部；

通过以上计算，得到气体在液相中的流动行为；

步骤9：连续性方程如公式(19)

式中

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

t——时间s；

对于液相，将其视为不可压缩流体，只在计算压力过程中添加相关的弱压缩系数；

步骤10：动量方程如公式(20)

式中

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

t——时间s；

P——粒子压力Pa；

μ——粒子动力粘度N·s/m²；

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

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

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

对于压力计算，采用显示压力模型进行计算，如公式(21)

式中

P_i^k+1——i粒子下一个时刻的压力Pa；

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

c——人造音速m/s，一般取为最大粒子速度的3倍；

B——显示压力计算模型系数，一般为7；

n_i——i粒子的粒子数密度；

n₀——初始的粒子数密度；

对于最大作用半径范围内不含气体粒子的粒子，采用压力梯度模型如公式(22)计算；对于附近含气体粒子的粒子，采用公式(14)计算；

式中

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

P——粒子压力Pa；

d——维度；

n₀——初始的粒子数密度；

——j粒子位置矢量；

——i粒子位置矢量；

ρ_i——i粒子密度kg/m³；

ρ_j——j粒子密度kg/m³；

P_j——j粒子压力Pa；

P_i,min——i粒子所有邻居粒子中压力的最小值Pa；

P_i^k+1——下一时刻i粒子的压力Pa；

α——人工可压缩系数，取值在10^-9到10^-7；

Δt——时间步长s；

——j粒子对i粒子的核函数值，表达形式如公式(2)；

对于粘度计算，采用粘度变化模型，如公式(23)

μ＝μ₀exp(2.5Aγ)公式(23)

式中

μ——粒子动力粘度N·s/m²；

μ₀——初始动力粘度N·s/m²；

A——粘度变化系数，对于Zr和UO₂设为3.0；对于不锈钢和熔融混凝土设为2.0；

γ——粒子固相率；

粘度项计算采用公式(13)；

对于固-液，液-液界面之间的表面张力采用基于自由能模型的表面张力模型计算，如公式(24)

式中

F——自由能系数，对于液-液界面如公式(25)，对于固-液界面如公式(26)；

l_ij——i粒子和j粒子的距离；

l_min——i粒子与周围的粒子的最小距离，采用1.5l₀；

r_e——粒子作用半径；

式中

F_ff——液-液界面的自由能系数

F_fs——固-液界面的自由能系数；

σ——粒子表面张力系数；

θ——粒子接触角°；

步骤11：速度和位置估算，计算完动量方程公式(20)中的重力项、粘性项、表面张力项后，对速度和位置进行估算，如公式(27)和公式(28)；

式中

——i粒子的估算速度矢量m/s；

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

μ——粒子动力粘度N·s/m²；

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

ρ_i——i粒子密度kg/m³；

——i粒子的表面张力矢量N/kg；

——粒子重力加速度矢量m/s²；

——i粒子的估算位置矢量m；

——当前时刻i粒子的位置矢量m；

Δt——时间步长s；

步骤12：速度和位置修正，计算完动量方程公式(20)中的压力项，对速度和位置进行修正，如公式(29)和公式(30)；

式中

——i粒子的估算速度矢量m/s；

——下一时刻i粒子的速度矢量m/s；

——i粒子的估算位置矢量m；

——下一时刻i粒子的位置矢量m；

ρ_i——i粒子密度kg/m³；

Δt——时间步长s；

P——粒子压力Pa；

通过步骤9至步骤12，得到了每个粒子的速度和位置，即得到了所有熔融物和混凝土的速度和位置，由此模拟得到了熔融物在和混凝土相互作用过程中的运动过程；

综上，通过步骤1对熔融物和混凝土的相互作用过程中熔融物和混凝土的位置、速度和初始物性参数进行设定；通过步骤5模拟熔融物和混凝土相互作用过程中熔融池内的液相粒子的传热、熔融池与固体混凝土的接触界面的传热、熔融物和混凝土的相变过程，计算得到每个粒子在不同时刻下的种类、焓值和温度，即得到熔融物和混凝土的相态、焓值和温度随时间的变化过程；通过步骤6模拟熔融物中UO₂、Zr、不锈钢两两之间的共晶反应，得到熔融物粒子中UO₂、Zr、不锈钢的物质摩尔份额的变化，即熔融物内的物质分布情况，并通过物质分布，更新熔融物的物性参数；通过步骤7模拟混凝土的分解反应过程和熔融物的氧化反应过程，计算得到每个粒子内的物质份额变化，即得到混凝土的分解反应过程和熔融物的氧化反应过程中物质的消失和生成过程；通过步骤8模拟气体的生成过程，计算得到气体粒子的出现和位置、速度的变化，即熔融池内由于混凝土分解产生的气体的生成和生长过程；通过步骤9至步骤12模拟熔融物和混凝土的运动过程，计算得到熔融物和混凝土粒子的速度、位置和压力，即得到熔融物和混凝土在相互作用过程中的运动和压力变化过程；综合以上步骤，模拟了熔融物和混凝土的相互作用过程，得到了作用过程中熔融物、混凝土和不可凝气体的位置、速度、压力、相态、温度、焓值、物质原子份额随时间的变化，通过以上数据对熔融物和混凝土相互作用过程中的传热相变、共晶反应、化学方应、气体产生和生长过程展开机理性分析。