CN111832214A - 基于粒子法的核反应堆严重事故碎片床熔化过程模拟方法 - Google Patents

Abstract

基于粒子法的核反应堆严重事故碎片床熔化过程模拟方法，步骤如下：1、粒子建模，设定粒子初始布置和参数，用粒子簇表征固体碎片，且分别对粒子和碎片进行编号；2、划分背景网格，检索邻居粒子组；3、显式计算动量方程中粘性项、表面张力项和重力项，预估粒子速度和位移；4、采用泊松方程计算粒子压力，依据压力梯度模型修正粒子速度和位移；5、计算粒子焓值、温度和相态；6、计算碎片平动速度和旋转速度，校正碎片中固相粒子的速度和位移；7、输出计算结果。本发明用小尺度粒子代表大尺寸固体碎片，并对其进行编号，解决了粒子所属碎片的识别问题。此外，本发明依据两粒子间是否被其它粒子遮挡的判别方式，建立了固体碎片间辐射换热模型。

Description

基于粒子法的核反应堆严重事故碎片床熔化过程模拟方法

技术领域

本发明涉及核电厂严重事故碎片床熔化过程研究技术领域，具体涉及一种基于粒子法的核反应堆严重事故碎片床熔化过程模拟方法。

背景技术

核严重事故中，堆芯失去冷却，随着堆芯温度升高，燃料包壳、定位格架和控制棒等金属材料因共晶反应最先熔化，并且熔融物向下渗透，随着堆芯温度继续升高，燃料芯块等氧化物在低于各自组分熔点的温度下，逐渐被金属熔融物化学消熔。而燃料芯块、金属构件以及熔融物在冷却过程中凝结、粉碎、解体，以碎片形式堆积被称为碎片床。碎片床可存在于堆芯区域、压力容器下封头和安全壳内。其结构复杂，且空间分布不均匀，具有多组分、多孔隙和多相态等特征，而当碎片不能得到有效冷却的情况下会再次熔化，并且熔融物在孔隙内迁移、形成硬壳、堵塞流道、形成熔池。对于压力容器下封头内的碎片床而言，其熔化过程的形态演变直接影响内部传热特性和壁面热流密度，进而影响压力容器的破损模式和熔融物泄漏的状态，是反应堆严重事故中的重要现象。但是由于目前的研究对碎片床熔化现象的机理揭示不够深入，所以研究碎片床熔化的动态过程，对核电厂严重事故安全特性研究具有重要意义。

碎片床熔化是反应堆严重事故的重要现象，对堆芯冷却、压力容器下封头破损和熔融物泄漏等事故进程具有重要影响，但由于反应堆燃料组件损伤和碎片床熔化现象复杂且具有不确定性。所以对于碎片床熔化和熔融物迁移动态过程的模拟，采用传统基于网格划分的方法会显出诸多弊端，比如难以准确捕捉复杂条件下的多组分与相界面，难以体现碎片床的多孔隙特征。而移动粒子半隐式方法作为无网格粒子法的一种，在处理复杂边界以及流体大变形、破碎与融合等问题时会具有明显优势。所以本文发明了一种基于粒子法的碎片床熔化过程模拟方法。

发明内容

碎片床中的固体碎片的熔化、分层和迁移直接影响熔体的传热特性，进而影响包括压力容器故障模式在内的后续事故进程，这是核严重事故研究的关键问题。本发明为了研究碎片床熔化的过程，基于移动粒子半隐式方法，结合碎片之间的辐射传热模型，熔融相变模型和表面张力的模型，对碎片材料的流动、传热、传质、相变进行研究，进而分析碎片床中的熔融分层和迁移过程。本发明提出一种基于粒子法的核反应堆严重事故碎片床熔化过程模拟方法，该方法采用小尺度的粒子来代表较大尺寸的固体碎片，提高了计算精度，同时通过对每个固体碎片进行编号，解决了粒子所属固体碎片的识别问题。

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

基于粒子法的核反应堆严重事故碎片床熔化过程模拟方法，步骤如下：

步骤1：对碎片床熔化过程初始状态进行粒子建模；用不同类别的粒子代表不同的物质，采用0、1、2号粒子分别模拟熔融物粒子的液相、固液混合相、固相，采用3、4、5号粒子分别模拟混凝土的液相、固液混合相、固相，并且所有的粒子根据所代表的物质具有其对应的物性参数，物性参数包括质量、密度、比热、熔沸点、温度和焓；对每一个粒子依次进行编号，采用由粒子组成的粒子簇来表征固体碎片，并且对粒子簇所代表的固体碎片进行二次编号；用分别代表金属与氧化物固体碎片的粒子簇依次交错布置在压力容器下封头内，以此来模拟压力容器下封头内碎片床初始位置分布；

步骤2：对所建立的粒子模型区域划分成若干个均匀布置的正方形背景网格；背景网格的边长与粒子作用半径相等，背景网格不参与实际运算，只用于划分区域及检索粒子；

步骤3：利用核函数来表征粒子间相互作用的程度；

步骤4：计算粒子规则分布情况下的初始粒子数密度n₀，来替代流体的真实密度；

步骤5：估算临时速度场，通过显式求解动量方程中的粘性项、重力项和表面张力项，对每个粒子的临时速度进行估算；

在求解粘性项时，由于在碎片床熔化过程中粒子会发生相变，且粒子的粘性与固相率有关，所以对于粘度计算，采用粘度变化模型，表达式如公式(7)所示：

μ_app＝μ₀exp(2.5γC_R) 公式(7)

μ_app——相变过程中材料的近似动态粘度N·s/m²；

μ₀——液相线温度下流体的动态粘度N·s/m²；

C_R——经验系数，取值为4到8；

γ——粒子固相率；

显式求解表面张力项时，对于固-液、液-液的表面张力模型，通过微分粒子间势能来计算表面张力；

步骤6：通过步骤5计算出的粒子临时速度来估算粒子临时位置；

步骤7：计算粒子临时移动后的粒子数密度，然后建立并隐式求解压力泊松方程，更新压力场后，求解动量方程中压力梯度项；求解压力时，采用的是一个改进后的压力泊松方程，如公式(17)所示：

式中

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

P——粒子压力Pa；

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

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

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

n^*——移动粒子后的粒子数密度；

Δt——时间步长s；

β、γ——截断误差补偿系数，0≤r≤β≤1；

其中系数β、γ用于调整误差补偿部分，能有效抑制压力波动并减少颗粒数密度误差；

步骤8：利用求解出的压力梯度项对粒子速度和位置进行修正；

步骤9：基于能量守恒方程求解辐射换热模型，该能量守恒方程结合了传热模型、相变模型和辐射换热模型，如公式(23)所示：

式中

h——粒子焓值J/kg；

T——温度K；

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

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

在熔化过程中，只计算碎片表面粒子之间的辐射换热；首先应用公式(24)，计算目标粒子i的粒子数密度，在计算n_i时，仅计算与目标粒子i在同一碎片上的邻居粒子j，再应用公式(25)，检索碎片上的表面粒子；

碎片表面粒子之间是否进行辐射换热，取决于两粒子在空间上是否被其它粒子遮挡；首先目标粒子i与有效半径为3r_e的区域内的任一表面粒子j构成一条线段l_ij；然后计算有效区域内其它粒子k到该线段l_ij之间的距离

若

小于粒子半径，如公式(26)所示，并且∠_ijk和∠_jik都是锐角，则认为粒子i和j被粒子k所遮挡，则不计算粒子i与粒子j之间的辐射换热，反之则用斯忒藩-玻尔兹曼定律公式(27)计算粒子i与粒子j之间的辐射换热；

n_i＝∑_j≠iw(|r_j-r_i|) 公式(24)

n_i<β₀n₀ 公式(25)

式中

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

——k粒子到线段l_ij的距离m；

r_e——粒子作用范围m；

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

ε——辐射发射率；

σ₀——斯蒂芬-玻尔兹曼常数W·m^-2·K^-4；

d₀——粒子直径m；

T_i——目标i粒子的温度K；

T_j——j粒子的温度K；

β₀——自由表面判定系数，取0.9；

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

r_j——j粒子位置矢量；

r_i——i粒子位置矢量；

步骤10：对于传热过程，采用拉普拉斯模型公式对公式(23)导热微分方程中拉普拉斯算子进行离散，可以得到传热模型，如公式(28)所示：

式中

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

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

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

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

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

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

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

焓值为温度的函数，如公式(30)所示：

h(T)＝ρc_P(T-T₀)+h₀ 公式(30)

式中

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

T₀——基准温度，基准温度设置为0K；

h₀——基准焓，基准焓设置为0J/kg；

将公式(30)代入公式(28)，得到粒子温度计算方程式，如公式(31)所示：

式中

——下一时刻i粒子的温度K；

步骤11：由于在熔化过程中，粒子会发生相变，所以采用公式(32)求解粒子固态、液态以及相变阶段的温度和焓之间的关系：

式中

T_solidus——固相线温度K；

T_liquidus——液相线温度K；

c_P,s——固态粒子比热容J/(kg·K)；

c_P,l——液态粒子比热容J/(kg·K)；

c_P,sl——相变状态粒子比热容J/(kg·K)；

ρ_s——固态粒子密度kg/m³；

ρ_l——液态粒子密度kg/m³；

ρ_sl——相变状态粒子密度kg/m³；

h_s＝ρ_sc_P,sT_solidus，粒子固相线温度对应焓值J/kg；

h_l＝h_s+ρ_slc_P,sl(T_liquidus-T_solidus)，粒子液相线温度对应焓值J/kg；

ρ_sl＝(ρ_s+ρ_l)/2；

h_sl——相变状态粒子潜热J/kg；

而物质所处的相态则由焓值定义的固相率来表示，如公式(33)所示，用于跟踪粒子相态。

当粒子呈液态即γ＝0，其速度与位置由动量守恒方程控制，当粒子在相变过程中即0<γ<1，粒子像其他液态粒子一样通过动量守恒方程移动粒子；当粒子呈固态即γ＝1时，其速度与位置除了由动量守恒方程控制，还要考虑到碎片中其他固相粒子的偶联作用；

所以在以下的步骤中，仅对碎片中固相粒子执行附加过程，来校正它们的速度以维持碎片的完整构型；

步骤12：检索碎片上所有固相粒子，形成固相粒子组{i_i,i_i+1,i_i+2,i_i+3,……i_i+n-1}，采用公式(34)，计算碎片中间的粒子位置r_g,则单个碎片中任一固相粒子相对于中心粒子的坐标为q_i，如公式(35)所示，碎片的转动惯量为I_i，如公式(36)所示。

q_i＝r_i-r_g 公式(35)

式中

r_i——碎片中第i个粒子位置矢量m；

N——碎片中总的固相粒子数；

r_g——碎片中间粒子的位置矢量m；

q_i——碎片中第i个固相粒子相对于中心粒子的位置矢量m；

I_i——第i个碎片的转动惯量kg·m²；

步骤13：计算碎片的平动速度T_i和旋转速度R_i，如公式(37)和公式(38)所示：

式中

T_i——第i个碎片的平动速度m/s；

R_i——第i个碎片的旋转速度rad/s；

u_i——碎片中第i个固相粒子速度矢量m/s；

步骤14：通过计算碎片的平动速度T_i和旋转速度R_i，来重新校正碎片中各粒子的速度u_i′，如公式(39)所示：

u_i′＝T_i+q_i×R_i 公式(39)

式中

u_i′——碎片矫正后第i个粒子速度矢量m/s；

综上，通过步骤1中碎片床固体碎片的位置、速度和初始物性参数进行设定；通过步骤2-8计算得到碎片床中粒子的速度、位置和压力，即得到碎片床在熔融过程中的运动和压力变化过程；通过步骤9、10和11，模拟碎片床熔融过程中不同相粒子之间的传热、熔融物的相变过程，计算得到每个粒子在不同时刻下的种类、焓值和温度，即得到不同相粒子的相态、焓值和温度随时间的变化过程；通过步骤12、步骤13、步骤14，仅考虑固相粒子之间的偶联作用，来校正固相粒子的速度以维持固体碎片的完整构型；综合以上步骤，模拟了碎片床熔化过程，得到了熔化过程中熔融物的位置、速度、压力、焓值、温度和相态随时间的变化。

本发明方法为碎片床熔化过程碎片时间的传热问题提供了解决方案，解决了粒子所属固体碎片的识别问题，也解决了粒子法中固体碎片间辐射换热的计算问题，为核电厂反应堆严重事故的安全特性研究提供了重要依据。

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

1、本发明提出了采用若干个粒子来表征固体碎片的方法，并对每个固体碎片进行编号来辨别粒子归属问题，在计算过程中计算液态熔融物和固体碎片对粒子的作用力，并综合固体碎片所有固相粒子的受力来计算整个固体碎片扭矩和角动量，从而计算固体碎片的移动速度和转动速度。然后根据计算所得到的固体碎片的移动速度和转动速度来矫正固体碎片中粒子的运动状态，从而维持固体碎片的完整构型。该处理方式采用了小尺度的粒子来代表较大尺寸的固体碎片，提高了计算精度，同时通过对每个固体碎片进行编号，可以计算固体碎片的受力情况和移动速度，也解决了粒子所属固体碎片的识别问题。

2、本发明提出了一种计算碎片表面辐射换热的方法，在有效半径为3r_e的区域内，目标粒子i与任一表面粒子j构成一条线段l_ij，然后计算有效区域内其它粒子k到该线段l_ij之间的距离

若

小于粒子半径，并且∠_ijk和∠_jik都是锐角，则认为粒子i和j被粒子k所遮挡，则不计算粒子i与粒子j之间的辐射换热，反之则用斯忒藩-玻尔兹曼定律计算粒子i与粒子j之间的辐射换热，该处理方式解决了粒子法中固体碎片间辐射换热的计算问题。

附图说明

图1是本发明基于粒子法的碎片床熔化过程模拟方法的流程图。

图2是本发明中碎片床初始状态粒子建模图。

图3是背景网格划分及设定邻居粒子组示意图。

图4是辐射换热作用范围图。

图5是压力容器下封头碎片床熔融、迁移和分层图。

具体实施方式

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

本发明基于粒子法的碎片床熔化过程模拟方法，如图1所示，步骤如下：

步骤1：对碎片床熔化过程初始状态进行粒子建模，用不同类别的粒子代表不同的物质，采用0、1、2号粒子分别模拟熔融物粒子的液相、固液混合相、固相，采用3、4、5号粒子分别模拟混凝土的液相、固液混合相、固相，并且所有的粒子根据所代表的物质具有其对应的物性参数，物性参数包括质量、密度、比热、熔沸点、温度、焓。对每一个粒子依次进行编号。本发明采用由粒子组成的粒子簇来表征固体碎片，并且对粒子簇所代表的固体碎片进行二次编号。用分别代表金属与氧化物固体碎片的粒子簇依次交错布置在压力容器下封头内，以此来模拟压力容器下封头内碎片床初始位置分布，如图2所示。

步骤2：对所建立的粒子模型区域划分成若干个均匀布置的正方形背景网格。针对某一粒子i，其位置坐标为(x_i,y_i)，需计算粒子i所在网格及其邻近网格内的粒子与粒子i之间的距离，当粒子i与任一粒子j的距离小于粒子作用半径r_e时，则将粒子j设置为粒子i的邻居粒子，经计算可得任一粒子的邻居粒子组，邻居粒子组内的粒子数设为n_neigh，该检索粒子过程会得到每个粒子的邻居粒子组{j₁,j₂,j₃,j₄,……j_neigh}，本发明中r_e＝3.1l₀(l₀为粒子间最小距离)，如图3所示，粒子之间距离由公式(1)计算：

步骤3：在移动粒子半隐式方法中，每个粒子均会和其相邻的粒子发生一定的相互作用，因此定义核函数来表征粒子间相互作用的程度。核函数是粒子间距离的函数，如公式(2)所示：

式中

r_e——粒子作用半径m；

r——粒子间距m；

w——核函数；

粒子间距离越小，粒子间的相互作用就会越大；相反，粒子间距离越大，粒子间相互作用就会越小。当粒子间距离粒子大于粒子作用半径时，认为粒子间相互作用可以忽略。

步骤4：在计算过程中，对于不可压缩流体，由于粒子数密度与流体密度线性相关，所以采用初始粒子数密度替代流体的真实密度，初始时刻，计算规则的粒子分布情况的粒子数密度，即初始粒子数密度n₀，如公式(3)所示：

式中

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

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

r_i——碎片中第i个粒子位置矢量；

r_j——目标粒子i邻近粒子j的位置矢量。

步骤5：显式估算临时速度场，通过求解动量方程公式(4)中的粘性项、重力项和表面张力项，对每个粒子的临时速度进行估算：

式中

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

P——粒子压力Pa；

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

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

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

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

t——时间s；

D——随体导数；

显式求解动量方程公式(4)中的粘性项，采用拉普拉斯模型(公式(5))对公式(4)中的拉普拉斯算子进行离散，如公式(6)所示：

式中

d——维度；

φ_i——i粒子在位置r_i处的标量；

φ_j——j粒子在位置r_i处的标量；

式中

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

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

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

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

——目标粒子i邻近粒子j当前的位置矢量；

由于在碎片床熔化过程中粒子会发生相变，而且粒子的粘性与固相率有关，需要考虑熔融物粘度随温度的变化过程，尤其对于固液混合态，其粘度随固相率的变化较大。所以对于粘度计算，采用粘度变化模型，表达式如公式(7)所示：

μ_app＝μ₀exp(2.5γC_R) 公式(7)

μ_app——相变过程中材料的近似动态粘度N·s/m²；

μ₀——液相线温度下流体的动态粘度N·s/m²；

C_R——经验系数，取值为4到8；

γ——粒子固相率；

显式求解动量方程公式(4)中的表面张力项，对于固-液，液-液的表面张力模型，本发明采用的表面张力模型是通过粒子间势能进行计算，表面张力为势能的微分，表达式如公式(8)所示：

式中

E_p——粒子间势能；

n——界面表面法向量；

表面势能如公式(9)所示：

式中

C——根据不同物质计算出的拟合参数，对于液-液界面如公式(11)，对于固-液界面如公式(12)；

l₀——初始粒子最小间距m；

对公式(8)进行微分，得到表面张力计算公式(10)如下：

式中

r_ij——i粒子与其邻居粒子j之间的间距m；

C_ff——液-液界面的拟合参数；

C_fs——固-液界面的拟合参数；

θ——流体与固体的接触角；

σ——表面张力系数Nm^-1；

当粒子间距离小于l₀时粒子间作用力为斥力，当粒子间距离大于l₀为吸引力当粒子距离超过粒子间作用力的有效半径时，粒子间作用力非常小，可以忽略。

显式求解动量方程公式(4)中的重力项，如公式(13)所示：

G＝ρg 公式(13)

式中

G——重力矢量N/kg；

通过公式(6)、公式(10)和公式(13)，求解动量方程公式(4)中的粘性项、表面张力项和重力项，再通过公式(14)对粒子的临时速度进行估算：

式中

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

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

Δt——时间步长s；

步骤6：估算粒子位置，通过步骤5计算出的粒子临时速度来估算粒子临时位置，如公式(15)所示：

式中

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

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

步骤7：计算动量方程公式(4)中的压力项。对于不可压缩流体，它的密度应该保持为常数，由于粒子数密度与流体密度线性相关，所以粒子数密度应为常数。但是粒子临时移动后，粒子数密度n*不再等于初始粒子数密度n⁰，需要用质量守恒方程(公式(16))对其进行隐式修正。再联立质量守恒方程(公式16)、动量方程(公式4)以及拉普拉斯算子模型(公式5)，得到流体速度场及压力场的耦合——压力泊松方程，如公式(17)所示，然后隐式求解压力泊松方程：

式中

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

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

n^*——移动粒子后的粒子数密度；

β、γ——截断误差补偿系数，0≤r≤β≤1

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

——移动粒子后j粒子位置矢量；

——移动粒子后i粒子位置矢量；

公式(17)是一个改进的压力泊松方程，其中系数β、γ用于调整误差补偿部分，适当的β、γ对可有效抑制压力波动并减少颗粒数密度误差。

对公式(17)等式等号左侧的拉普拉斯算子通过拉普拉斯模型(式5)进行离散，如公式(19)所示：

式中

——下一时刻j粒子压力Pa；

——下一时刻目标i粒子的压力值Pa；

结合方程(17)和(19)，可以得出一个对称矩阵来计算压力。通过共轭梯度法求解矩阵，并将用于检验迭代收敛性的值设置为1.0×10^-9。

通过求解压力泊松方程更新压力场后，动量方程公式(4)中的压力梯度项可采用公式(20)进行计算：

式中

P_j——j粒子压力Pa；

——目标i粒子及其邻居粒子之间的最小压力值Pa；

采用所有邻居粒子中压力的最小值来替代目标粒子的压力，可以有效避免粒子的聚集行为，更好地保证计算的稳定性。

步骤8：修正粒子速度和位置，粒子速度和位置通过压力梯度项进行修正，如公式(21)和公式(22)所示：

式中

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

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

在每个时间层内，先通过步骤5显式计算动量守恒方程中的重力项、粘性项以及表面张力项，获得粒子速度的估算值；通过步骤6获得粒子位置的估算值；通过步骤7求解速度压力耦合关系式—

—泊松方程，得到压力场；通过步骤8利用压力梯度项修正得到粒子在下一时层的速度及位置。而且动量方程中压力项的计算是模拟流体不可压缩条件的重要内容，也是整个移动粒子半隐式方法求解过程中的关键步骤。

步骤9：碎片床熔化过程中存在传热相变过程，为了能够精确再现该过程，必须开发精确的能量守恒计算模型。本发明中的能量守恒方程结合了传热模型、相变模型和辐射换热模型，如公式(23)所示：

式中

h——粒子焓值J/kg；

T——温度K；

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

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

对于辐射换热模型。在熔化过程中，只计算碎片表面粒子之间的辐射换热。首先应用公式(24)，计算目标粒子i的粒子数密度，在计算n_i时，仅计算与目标粒子i在同一碎片上的邻居粒子j，再应用公式(25)，检索碎片上的表面粒子。

碎片表面粒子之间是否进行辐射换热，取决于两粒子在空间上是否被其它粒子遮挡。如图4所示，首先目标粒子i与有效区域内(有效半径为3r_e)的任一表面粒子j构成一条线段l_ij；然后计算有效区域内其它粒子k到该线段l_ij之间的距离

若

小于粒子半径(公式26)，并且∠_ijk和∠_jik都是锐角，则认为粒子i和j被粒子k所遮挡，则不计算粒子i与粒子j之间的辐射换热，反之则用斯忒藩-玻尔兹曼定律(公式27)计算粒子i与粒子j之间的辐射换热。

n_i＝∑_j≠iw(|r_j-r_i|) 公式(24)

n_i<β₀n₀ 公式(25)

式中

——粒子k到线段l_ij的距离m；

ε——辐射发射率；

σ₀——斯蒂芬-玻尔兹曼常数W·m^-2·K^-4；

d₀——粒子直径m；

T_i——目标i粒子的温度K；

T_j——j粒子的温度K；

β₀——自由表面判定系数，取0.9。

步骤10：对于传热过程，采用拉普拉斯模型公式(5)对公式(23)导热微分方程中拉普拉斯算子进行离散，可以得到传热模型，如公式(28)所示：

式中

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

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

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

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

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

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

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

焓值为温度的函数，如公式(30)所示：

h(T)＝ρc_P(T-T₀)+h₀ 公式(30)

式中

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

T₀——基准温度，基准温度设置为0K；

h₀——基准焓，基准焓设置为0J/kg；

将公式(30)代入公式(28)，可以得到粒子温度计算方程式，如公式(31)所示：

式中

——下一时刻i粒子的温度K；

该传热模型适用于固体粒子与固体粒子之间的导热、流体粒子和流体粒子之间的传热以及流体粒子和固体粒子之间的传热。

步骤11：相变模型。相对于气体和液体之间的密度差异，固体和液体密度差别非常小，所以，相变模型中仅采用两种拥有不同物性的粒子分别代表同一种物质的固态和液态，而不考虑相变过程中的体积变化。对于堆芯熔融物中金属物质和氧化物都为混合物，其发生相变时温度并不固定，而是在一定温度范围发生变化。所以固态、液态以及相变阶段的温度和焓之间的函数如公式(32)所示：

式中

T_solidus——固相线温度K；

T_liquidus——液相线温度K；

c_P,s——固态粒子比热容J/(kg·K)；

c_P,l——液态粒子比热容J/(kg·K)；

c_P,sl——相变状态粒子比热容J/(kg·K)；

ρ_s——固态粒子密度kg/m³；

ρ_l——液态粒子密度kg/m³；

ρ_sl——相变状态粒子密度kg/m³；

h_s＝ρ_sc_P,sT_solidus，粒子固相线温度对应焓值J/kg；

h_l＝h_s+ρ_slc_P,sl(T_liquidus-T_solidus)，粒子液相线温度对应焓值J/kg；

ρ_sl＝(ρ_s+ρ_l)/2；

h_sl——相变状态粒子潜热J/kg；

由焓值定义的固相率来表示物质所处的相态，如公式(33)所示，用于跟踪粒子状态。

当焓值低于h_s或高于h_l时，相态分别为固态(γ＝1)或液态(γ＝0)。否则，碎片处于相变状态。当粒子呈液态(γ＝0)，其速度与位置由动量守恒方程控制，当粒子在相变过程中(0<γ<1)，粒子像其他液态粒子一样通过动量守恒方程移动粒子。当粒子呈固态时，其速度与位置除了由动量守恒方程控制，还要考虑到碎片中其他固相粒子的偶联作用。

通过步骤9、10和11的计算，模拟碎片床熔融过程中粒子之间的传热及熔融物的相变过程。计算得到每个粒子在不同时刻下的种类、焓值和温度，即得到碎片床中熔融物的相态、焓值和温度随时间的变化过程。

上述步骤中由于没有考虑固体碎片中单个固相粒子之间的偶联作用，因此固体碎片不会立即移动。所以在接下来步骤中，仅对碎片中固相粒子执行附加过程，来校正它们的速度以维持碎片的完整构型。

步骤12：检索碎片上所有固相粒子，形成固相粒子组{i_i,i_i+1,i_i+2,i_i+3,……i_i+N-1}，采用公式(34)，计算碎片中间的粒子位置r_g,则单个碎片中任一固相粒子相对于中心粒子的坐标为q_i，如公式(35)所示，碎片的转动惯量为I_i，如公式(36)所示。

q_i＝r_i-r_g 公式(35)

式中

N——碎片中总的固相粒子数；

r_g——碎片中间粒子的位置矢量m；

q_i——碎片中第i个固相粒子相对于中心粒子的位置矢量m；

I_i——第i个碎片的转动惯量kg·m²；

步骤13：计算碎片的平动速度T_i和旋转速度R_i，如公式(37)和公式(38)所示：

式中

T_i——第i个碎片的平动速度m/s；

R_i——第i个碎片的旋转速度rad/s；

u_i——碎片中第i个固相粒子速度矢量m/s；

步骤14：通过计算碎片的平动速度T_i和旋转速度R_i，来重新校正碎片中各粒子的速度，如公式(39)所示：

u_i′＝T_i+q_i×R_i 公式(39)

式中

u_i′——校正后碎片中的第i个粒子速度矢量m/s；

通过步骤12、步骤13、步骤14的计算，可以模拟液态熔融物和固体碎片对粒子的作用力，并综合固体碎片所有粒子的受力计算整个固体碎片扭矩和角动量，从而计算固体碎片的移动速度和转动速度。通过计算所得的固体碎片的移动速度和转动速度可以重新矫正每个粒子的速度。

综上，通过步骤1中碎片床固体碎片的位置、速度和初始物性参数进行设定；通过步骤2-8计算得到碎片床中粒子的速度、位置和压力，即得到碎片床在熔融过程中的运动和压力变化过程。通过步骤9、10和11，模拟碎片床熔融过程中不同相粒子之间的传热、熔融物的相变过程，计算得到每个粒子在不同时刻下的种类、焓值和温度，即得到不同相粒子的相态、焓值和温度随时间的变化过程；通过步骤12、步骤13、步骤14，仅考虑固相粒子之间的偶联作用，来校正固相粒子的速度以维持固体碎片的完整构型。综合以上步骤，模拟了碎片床熔化过程，得到了熔化过程中熔融物的位置、速度、压力、焓值、温度和相态随时间的变化。

整合以上数据可以得到熔融物迁移、分层以及硬壳的分布。随着碎片温度的逐渐升高，中心部位的金属碎片最先熔化，并且熔融物向下渗透于碎片床的孔隙内，过程中不断释放热量，如图5(a)所示；熔融物向下渗透过程中，凝固形成了近似碗状的硬壳，如图5(b)所示；随后金属熔融物在硬壳内逐渐堆积，且熔融物的温度继续升高，硬壳也逐渐向下扩展，随后金属熔融物逐渐化学消熔氧化物，在金属熔融物与氧化物熔融物交界面形成一层氧化物硬壳，如图5(c)所示；最后除了紧挨下封头的边界颗粒未熔化，其他颗粒全部熔化，金属熔融物与氧化物分层，金属熔融物在上，氧化物融物在下，并且氧化物熔融物被由氧化物形成得硬壳所包围，其余金属熔融物也被由金属形成的硬壳包围，如图5(d)所示。

Claims (1)

1.基于粒子法的核反应堆严重事故碎片床熔化过程模拟方法，其特征在于：步骤如下：

步骤1：对碎片床熔化过程初始状态进行粒子建模；用不同类别的粒子代表不同的物质，采用0、1、2号粒子分别模拟熔融物粒子的液相、固液混合相、固相，采用3、4、5号粒子分别模拟混凝土的液相、固液混合相、固相，并且所有的粒子根据所代表的物质具有其对应的物性参数，物性参数包括质量、密度、比热、熔沸点、温度和焓；对每一个粒子依次进行编号，采用由粒子组成的粒子簇来表征固体碎片，并且对粒子簇所代表的固体碎片进行二次编号；用分别代表金属与氧化物固体碎片的粒子簇依次交错布置在压力容器下封头内，以此来模拟压力容器下封头内碎片床初始位置分布；

步骤2：对所建立的粒子模型区域划分成若干个均匀布置的正方形背景网格；背景网格的边长与粒子作用半径相等，背景网格不参与实际运算，只用于划分区域及检索粒子；

步骤3：利用核函数来表征粒子间相互作用的程度；

步骤4：计算粒子规则分布情况下的初始粒子数密度n₀，来替代流体的真实密度；

步骤5：估算临时速度场，通过显式求解动量方程中的粘性项、重力项和表面张力项，对每个粒子的临时速度进行估算；

在求解粘性项时，由于在碎片床熔化过程中粒子会发生相变，且粒子的粘性与固相率有关，所以对于粘度计算，采用粘度变化模型，表达式如公式(7)所示：

μ_app＝μ₀ exp(2.5γC_R) 公式(7)

μ_app——相变过程中材料的近似动态粘度N·s/m²；

μ₀——液相线温度下流体的动态粘度N·s/m²；

C_R——经验系数，取值为4到8；

γ——粒子固相率；

显式求解表面张力项时，对于固-液、液-液的表面张力模型，通过微分粒子间势能来计算表面张力；

步骤6：通过步骤5计算出的粒子临时速度来估算粒子临时位置；

步骤7：计算粒子临时移动后的粒子数密度，然后建立并隐式求解压力泊松方程，更新压力场后，求解动量方程中压力梯度项；求解压力时，采用的是一个改进后的压力泊松方程，如公式(17)所示：

式中

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

P——粒子压力Pa；

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

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

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

n^*——移动粒子后的粒子数密度；

Δt——时间步长s；

β、γ——截断误差补偿系数，0≤r≤β≤1；

其中系数β、γ用于调整误差补偿部分，能有效抑制压力波动并减少颗粒数密度误差；

步骤8：利用求解出的压力梯度项对粒子速度和位置进行修正；

步骤9：基于能量守恒方程求解辐射换热模型，该能量守恒方程结合了传热模型、相变模型和辐射换热模型，如公式(23)所示：

式中

h——粒子焓值J/kg；

T——温度K；

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

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

在熔化过程中，只计算碎片表面粒子之间的辐射换热；首先应用公式(24)，计算目标粒子i的粒子数密度，在计算n_i时，仅计算与目标粒子i在同一碎片上的邻居粒子j，再应用公式(25)，检索碎片上的表面粒子；

碎片表面粒子之间是否进行辐射换热，取决于两粒子在空间上是否被其它粒子遮挡；首先目标粒子i与有效半径为3r_e的区域内的任一表面粒子j构成一条线段l_ij；然后计算有效区域内其它粒子k到该线段l_ij之间的距离

若

小于粒子半径，如公式(26)所示，并且∠_ijk和∠_jik都是锐角，则认为粒子i和j被粒子k所遮挡，则不计算粒子i与粒子j之间的辐射换热，反之则用斯忒藩-玻尔兹曼定律公式(27)计算粒子i与粒子j之间的辐射换热；

n_i＝∑_j≠iw(|r_j-r_i|) 公式(24)

n_i<β₀n₀ 公式(25)

Q＝εσ₀(T_j ⁴-T_i ⁴)/d₀ 公式(27)

式中

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

——k粒子到线段l_ij的距离m；

r_e——粒子作用范围m；

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

ε——辐射发射率；

σ₀——斯蒂芬-玻尔兹曼常数W·m^-2·K^-4；

d₀——粒子直径m；

T_i——目标i粒子的温度K；

T_j——j粒子的温度K；

β₀——自由表面判定系数，取0.9；

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

r_j——j粒子位置矢量；

r_i——i粒子位置矢量；

步骤10：对于传热过程，采用拉普拉斯模型公式对公式(23)导热微分方程中拉普拉斯算子进行离散，可以得到传热模型，如公式(28)所示：

式中

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

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

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

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

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

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

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

焓值为温度的函数，如公式(30)所示：

h(T)＝ρc_P(T-T₀)+h₀ 公式(30)

式中

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

T₀——基准温度，基准温度设置为0K；

h₀——基准焓，基准焓设置为0J/kg；

将公式(30)代入公式(28)，得到粒子温度计算方程式，如公式(31)所示：

式中

——下一时刻i粒子的温度K；

步骤11：由于在熔化过程中，粒子会发生相变，所以采用公式(32)求解粒子固态、液态以及相变阶段的温度和焓之间的关系：

式中

T_solidus——固相线温度K；

T_liquidus——液相线温度K；

c_P,s——固态粒子比热容J/(kg·K)；

c_P,l——液态粒子比热容J/(kg·K)；

C_P,sl——相变状态粒子比热容J/(kg·K)；

ρ_s——固态粒子密度kg/m³；

ρ_l——液态粒子密度kg/m³；

ρ_sl——相变状态粒子密度kg/m³；

h_s＝ρ_sc_P,sT_solidus，粒子固相线温度对应焓值J/kg；

h_l＝h_s+ρ_slc_P,sl(T_liquidus-T_solidus)，粒子液相线温度对应焓值J/kg；

ρ_sl＝(ρ_s+ρ_l)/2；

h_sl——相变状态粒子潜热J/kg；

而物质所处的相态则由焓值定义的固相率来表示，如公式(33)所示，用于跟踪粒子相态。

当粒子呈液态即γ＝0，其速度与位置由动量守恒方程控制，当粒子在相变过程中即0<γ<1，粒子像其他液态粒子一样通过动量守恒方程移动粒子；当粒子呈固态即γ＝1时，其速度与位置除了由动量守恒方程控制，还要考虑到碎片中其他固相粒子的偶联作用；

所以在以下的步骤中，仅对碎片中固相粒子执行附加过程，来校正它们的速度以维持碎片的完整构型；

步骤12：检索碎片上所有固相粒子，形成固相粒子组{i_i,i_i+1,i_i+2,i_i+3,……i_i+n-1}，采用公式(34)，计算碎片中间的粒子位置r_g,则单个碎片中任一固相粒子相对于中心粒子的坐标为q_i，如公式(35)所示，碎片的转动惯量为I_i，如公式(36)所示。

q_i＝r_i-r_g 公式(35)

式中

r_i——碎片中第i个粒子位置矢量m；

N——碎片中总的固相粒子数；

r_g——碎片中间粒子的位置矢量m；

q_i——碎片中第i个固相粒子相对于中心粒子的位置矢量m；

I_i——第i个碎片的转动惯量kg·m²；

步骤13：计算碎片的平动速度T_i和旋转速度R_i，如公式(37)和公式(38)所示：

式中

T_i——第i个碎片的平动速度m/s；

R_i——第i个碎片的旋转速度rad/s；

u_i——碎片中第i个固相粒子速度矢量m/s；

步骤14：通过计算碎片的平动速度T_i和旋转速度R_i，来重新校正碎片中各粒子的速度u_i′，如公式(39)所示：

u_i′＝T_i+q_i×R_i 公式(39)

式中

u_i′——碎片矫正后第i个粒子速度矢量m/s；

综上，通过步骤1中碎片床固体碎片的位置、速度和初始物性参数进行设定；通过步骤2-8计算得到碎片床中粒子的速度、位置和压力，即得到碎片床在熔融过程中的运动和压力变化过程；通过步骤9、10和11，模拟碎片床熔融过程中不同相粒子之间的传热、熔融物的相变过程，计算得到每个粒子在不同时刻下的种类、焓值和温度，即得到不同相粒子的相态、焓值和温度随时间的变化过程；通过步骤12、步骤13、步骤14，仅考虑固相粒子之间的偶联作用，来校正固相粒子的速度以维持固体碎片的完整构型；综合以上步骤，模拟了碎片床熔化过程，得到了熔化过程中熔融物的位置、速度、压力、焓值、温度和相态随时间的变化。

Images