[发明专利]基于粒子法的熔融物与混凝土相互作用分析方法

权利要求书

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的物质摩尔份额其中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^-9到10^-7；

σ——表面张力系数；

κ_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，得到了每个粒子的速度和位置，即得到了所有熔融物和混凝土的速度和位置，由此模拟得到了熔融物在和混凝土相互作用过程中的运动过程。