CN112102894A - 基于粒子法的核反应堆堆芯材料熔池演变特性分析方法 - Google Patents

Abstract

一种基于粒子法的核反应堆堆芯材料熔池演变特性分析方法，主要步骤如下：1、几何模型建模；2、物性初始化；3、计算粒子间的核函数和粒子数密度；4、计算物质混合过程，更新粒子物质含量、焓值和物性；5、计算化学反应过程，更新粒子物质含量和焓值；6、计算传热过程，更新粒子焓值；7、根据焓值，计算粒子温度和相态；8、显式计算粘性项和重力项；9、显式计算湍流应力项，估算粒子速度和位置，计算湍流热流项，更新粒子焓值；10、隐式迭代计算粒子压力；11、根据压力值，修正粒子速度和位置；12、输出计算结果。本方法考虑核反应堆堆芯材料熔池内的所有现象；基于粒子法，能够精确捕捉熔池内的物质变化、相态变化和流动特性。

Description

基于粒子法的核反应堆堆芯材料熔池演变特性分析方法

技术领域

本发明涉及核反应堆堆芯材料熔池演变特性分析研究技术领域，具体涉及一种基于粒子法的核反应堆堆芯材料熔池演变特性分析方法。

背景技术

在核反应堆运行过程中，如果冷却剂无法完全带走堆芯产生的热量，堆芯则会偏离正常运行工况，严重情况下，堆芯温度会大幅上升甚至发生熔化。堆芯熔化过程中，堆芯内材料熔点较低或能够发生化学反应的材料会优先液化。在约800℃时，银铟镉合金发生熔化；在约940℃时，铁-锆和镍-锆共晶反应形成液态共晶产物；约1150℃时，碳化硼和铁发生共晶反应；约1200℃，锆合金被水迅速氧化；约1450℃，不锈钢和铬镍铁合金开始熔化；约1760℃，锆-4合金开始熔化，且液态锆合金会与二氧化铀发生共晶反应，形成铀-锆-氧共晶熔融物；约1975℃，α-Zr(O)熔化；约2350℃，碳化硼熔化；约2690℃，氧化锆熔化；约2520℃，二氧化铀熔化。堆芯材料熔化过程中涉及多种物质，除了以上的物质外，还包含各种裂变产物。因此堆芯熔融物的组分复杂，其形成的熔池特性很难完全精确捕捉。

堆芯熔化后，熔融物会在压力容器下封头聚集形成熔池。目前针对我国先进核电站针对核反应堆严重事故大多采用熔融物对内滞留(IVR)，将堆芯熔融物滞留在压力容器下封头，从而阻止严重事故的进一步进行。这就要求必须提供有效的冷却方案来保证压力容器的完整性，压力容器外部冷却是实现IVR的有效措施。熔池的换热特性直接决定了压力容器壁面的热负荷，冷却措施能否保持压力容器完整性的重要判断依据之一。针对熔池的换热特性，国内外展开了大量的实验研究，大多采用水或盐水作为模拟物，获得相对高瑞利数情况下的熔池换热特性关系式，但这些实验无法模拟出熔融物在压力容器下封头内壁面上的凝固相变过程。由于熔池内物质材料种类繁多，大多数值模拟将熔池成分简化，将互溶的物质视为整体，主流的熔池结构模型有多层熔池模型和混合层熔池模型。多层熔池模型将熔池成分分为几种互不相溶的成分，主流有双层熔池结构(金属层+氧化物层)和三层熔池结构(轻金属层+氧化物层+重金属层)；混合层熔池模型将整体熔池视为具有等效物性的整体，不考虑熔池分层结构的影响。目前的数值模拟中很难考虑层与层之间的互溶过程，如重金属层析出形成过程，熔池结构倒转等现象。

此外，当压力容器失效后，熔融物会向核反应堆堆坑迁徙，并在堆坑处再次形成熔池。此时，熔融物会与安全壳内的混凝土直接接触，发生长期相互作用。混凝土会在高温熔融物的作用下分解或熔化，其分解产物会混入熔池，改变熔池的结构，影响较大的成分包括二氧化硅和氧化钙。随后熔融物与混凝土的作用过程，熔池内的物质组成会发生一系列改变，从而导致熔池结构发生改变。在反应初期，氧化物熔融物主要为二氧化铀和氧化锆，密度较大，熔池结构为上层金属熔融物，下层氧化物熔融物。随着混凝土分解的氧化物与氧化物熔融物混合后，氧化物熔融物层的密度会逐渐下降，当其密度小于金属熔融物时，熔池结构会发生倒置，变为上层氧化物熔融物，下层金属熔融物。随着金属熔融物被氧化，熔池结构逐渐模糊，不再存在明显的分层结构。不同熔池结构的换热特性和流动特性存在较大差异。目前的关于堆芯熔融物与混凝土反应时熔池结构改变的研究较少，大量实验关注混凝土烧蚀速率和冷却缓解措施，采用氧化物熔融物和金属熔融物开展实验，简单考虑了双层熔池与单一熔融物对混凝土烧蚀行为的差异。关于该过程的数值模拟也存在较大难度，大多简化考虑多层熔池的影响，无法很好的再现熔融物组分变化引起的熔池结构的改变。

基于拉格朗日方法的粒子法，在处理物质变化、能量变化及流动特性等方面具有独特的优势，精确计算熔池内多物质组分间的混合过程、化学反应过程及湍流流动行为。因此，本发明综合粒子法及核反应堆堆芯材料熔池演变特性提出一种核反应堆堆芯材料熔池演变特性分析方法。

发明内容

为了全面研究核反应堆堆芯材料熔池演变特性，揭示熔池内可能存在的一些机理现象，本发明在对核反应堆堆芯材料熔池的机理性分析的基础上结合粒子法、基础控制方程及相关的机理性化学物理模型，提出一种基于粒子法的核反应堆堆芯材料熔池演变特性分析方法，该方法能够对核反应堆堆芯材料熔池中存在的多种物质流动、传热、传质、相变、化学反应进行研究，获得熔池内熔融物物质变化、能量变化及流动特性，精确计算熔池内多物质组分间的混合过程、化学反应过程及湍流流动行为，有助于熔池分层结构的机理特性研究，并为核电厂反应堆严重事故安全特性研究提供重要依据。

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

一种基于粒子法的核反应堆堆芯材料熔池演变特性分析方法，步骤如下：

步骤1：几何模型建模，采用离散的粒子构建核反应堆芯材料熔池的几何模型，每个粒子携带有位置和速度信息，当粒子直径小于实际几何的一百分之一时，粒子构型能够精确还原真实几何的结构和运动特性；采用0、1、2号粒子表征液相、固相、固液混合相，采用3号粒子表征虚拟壁面粒子，虚拟壁面粒子保证壁面边界处的粒子数密度守恒，防止流体粒子向壁面的穿透；每个粒子携带有物质信息，核反应堆堆芯材料物质种类繁多，无法精确考虑所有物质，因此只考虑含量较高的物质：二氧化铀、不锈钢、锆和二氧化锆；此外，当压力容器失效后，熔融物进入安全壳，落入堆坑，在堆坑处重新形成熔池，该过程涉及熔融物与混凝土的相互作用，熔池内会混入混凝土材料，主要包括：二氧化硅和氧化钙；每个粒子携带有焓值信息和温度信息；

步骤2：粒子物性计算，设置每种物质的物性参数信息，根据每个粒子的物质信息，计算每个粒子的物性参数，包括密度、比热、熔点或固相线温度和液相线温度、沸点；

步骤3：定义粒子核函数，表征粒子间的相互作用程度，采用三阶尖峰核函数，如公式(1)所示，

式中

w——粒子核函数；

r——粒子间距m；

r_e——粒子作用半径m；

当粒子间距大于粒子作用半径时，忽略粒子间的微小作用力，认为粒子间不存在作用力；定义粒子数密度，表征粒子的密集程度，

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

式中

n——粒子数密度；

i——粒子ID；

j——粒子ID；

w——粒子核函数；

步骤4：物质混合计算，在核反应堆堆芯材料熔池内，物质种类繁多，不同物质与物质之间会相互混合，有些物质之间是互不相溶的，物质的互不相溶会使得熔池出现分层现象，物质的混合过程，决定了熔池结构特性，熔池的不同结构，其换热及流动特性会存在较大差异；物质的混合基于扩散定律，考虑两种物质之间的混合过程，如公式(3)所示，

式中

——粒子i下一时刻的物质浓度kg/m³；

——粒子i当前时刻的物质浓度kg/m³；

D——物质混合系数；

Δt——时间步长s；

d——维度；

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

——粒子j当前时刻的物质浓度kg/m³；

w_ij——粒子i和粒子j的核函数值；

r_ij——粒子i和粒子j的粒子间距；

物质的传递过程伴随着能量的传递，物质混合过程中的能量守恒如公式(4)和公式(5)所示，此处的物质混合过程不考虑化学反应；

式中

——粒子i下一时刻的体积焓J/m³；

——粒子i当前时刻的体积焓J/m³；

——物质x的焓J/kg；

T_x——物质x的温度K；

c_p,x——物质x的定压比热容J/(kg·K)；

T_s——固相线温度K；

Δh_sl,x——物质x的熔化潜热J/kg；

T_l——液相线温度K；

粒子的物性参数由各物质的物质浓度决定，其中粒子密度如公式(6)所示，粒子的比热、熔化潜热、汽化潜热、热导率、粘性如公式(7)所示，固相线温度和液相线温度根据不同的混合物的物性模型计算得到；

式中

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

x——物质种类号；

m_x——粒子中物质x的物质浓度kg/m³；

Y——物性参数；

步骤5：在堆芯材料熔池内，存在大量的金属，极易被氧化，主要考虑铁和锆的氧化反应：

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

基于以上化学反应方程式，当两个粒子间距足够小或粒子满足化学反应条件时，根据化学反应速率进行物质的转变，如公式(8)和公式(9)所示；同时以内热源的形式改变粒子的焓值来模拟化学热的过程，如公式(10)所示；物质种类改变后，同时对粒子物性进行更新计算；

式中

——粒子i下一时刻物质x的物质浓度kg/m³；

——粒子i当前时刻物质x的物质浓度kg/m³；

C_x,chem——物质x的变化速率kg/(m³·s)；

——粒子i下一时刻物质y的物质浓度kg/m³；

——粒子i当前时刻物质y的物质浓度kg/m³；

X——化学反应中物质x与物质y的物质变化比例；

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

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

Q_chem——化学反应热W/kg；

通过步骤4和步骤5，更新每个粒子的物质变化和焓值变化，并根据物质含量和焓值更新粒子的物性；

步骤6：传热计算，除了物质混合过程会伴随着能量传递外，粒子会通过对流换热、导热和辐射换热进行能量交换；在粒子法中，通过粒子的运动和导热过程的结合，还原对流换热过程，对流换热和导热过程采用公式(11)计算，辐射换热过程采用公式(12)和公式(13)计算；

式中

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

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

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

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

Δt——时间步长s；

h_r——辐射换热系数W/(m²·K)；

σ_stef——斯特藩-玻尔兹曼常数，5.67×10^-8W/(m²·K⁴)；

ε_i——自由表面粒子i的辐射发射率；

ε_env——周围环境的辐射发射率；

T_i——自由表面粒子i的温度K；

T_env——堆腔周围环境温度K；

Q_r——辐射换热量W/m³；

A——辐射换热面积m²；

l——粒子直径m；

通过步骤6，获得每个粒子下一时刻的物质焓值和温度；

步骤7：相变计算，通过能量守恒计算得到每个时刻下每个粒子的焓值，根据焓值计算得到粒子对应温度，如公式(14)所示，

式中

T——温度K；

h——焓值J/kg；

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

c_p——定压比热容J/(kg·K)；

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

定义固相率ε判断粒子的相态变化，固相率为焓值的函数：

式中

ε——固相率；

当固相率为1时，粒子完全为固态，当固相率为0时，粒子完全为液态，当固相率介于0和1之间时，粒子为固液混合态；固液混合态较为特殊，固相和液相的占比对粒子的行为特性影响较大；定义临界固相率ε_crit来判断固液混合态粒子的行为特征，当ε＞ε_crit时，粒子内固相占比较高，粒子行为趋向于固态；当ε＜ε_crit时，粒子内液相占比较高，粒子行为趋向于液态；临界固相率设为0.55；

通过步骤7计算得到每个粒子的温度和相态；

步骤8：质量和动量守恒计算，针对不可压缩流体展开计算，在质量守恒方程中忽略密度变化，如公式(16)所示；动量守恒主要考虑压力、重力、粘性力和湍流应力，如公式(17)所示，

式中

t——时间s；

——速度矢量m/s；

P——压力Pa；

ν——运动粘度m²/s；

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

——湍流应力N；

步骤9：湍流计算，在核反应堆熔池内，由于温度梯度较大，物质组分间密度差较大，熔池内部会存在剧烈的搅浑和对流现象，为了能够准确模拟熔池内的流动特性，采用大涡模拟方法对熔池湍流特性展开计算；以粒子尺寸为区分，将流体的湍流运动分为大尺度运动和小尺度运动，对于涡旋尺度大于粒子尺寸的，采用纳维斯-托克斯方程直接求解，对于涡旋尺度小于粒子尺寸的，采用湍流应力项进行考虑，因此将运动分为时均值和脉动值，如公式(18)所示，

式中

φ——瞬时值；

——时均值；

φ′——脉动值；

对于动量守恒方程中，如公式(17)，脉动值即为湍流应力项，对于能量守恒方程中，脉动值为湍流热流项，如公式(19)所示，

式中

T——温度K；

α——热扩散率m²/s；

——湍流热流项；

湍流应力项形式如公式(20)所示，为三阶矩阵，

式中

τ——湍流应力N；

x,y,z——表示湍流应力的不同方向；

每个方向的应力如公式(21)、公式(22)和公式(23)所示：

式中

τ_ab——ab方向上的湍流应力N；

u——速度m/s；

ν_t——涡粘性；

τ_γγ——三阶矩阵对角线上的切应力N；

δ_ab——颜色函数；

a——方向a；

b——方向b；

其中涡粘性采用Smagoringsky亚格子模型求解，如公式(24)、公式(25)、公式(26)、公式(27)、公式(28)和公式(29)所示，

式中

ν_t——涡粘性；

C_D——动态湍流系数；

Δ——滤波宽度；

Δ₂——二次滤波宽度；

——S的滤波值；

——u的滤波值；

——某参数滤波值；

r_ij——粒子i和粒子j的粒子间距m；

G——高斯核函数；

公式(22)中速度的偏导项根据梯度模型展开：

式中

u——速度m/s；

a——方向a；

s——长度m；

r_b,j——粒子j在方向b上的距离m；

r_b,i——粒子i在方向b上的距离m；

计算获得每个方向上的湍流应力后，湍流应力的散度项根据散度模型展开：

式中

——方向a上的切应力向量；

——粒子j在方向a上的切应力向量；

——粒子i在方向a上的切应力向量；

——粒子j的位置矢量；

——粒子i的位置矢量；

公式(19)中的湍流热流项计算如下：

式中

θ_a——方向a上的湍流热流；

Pr_t——普朗特数；

通过以上计算获得了湍流应力项和湍流热流项，从而能够计算得到每个粒子在考虑湍流影响的速度、位置和温度；

步骤10：压力泊松方程求解，针对公式(17)中的压力梯度项，采用隐式迭代求解压力泊松方程的方式计算得到每个粒子的位置和速度，粒子格式离散后的压力泊松方程如公式(33)所示，

式中

β——调节系数，常取0.4；

——粒子i当前时刻的粒子数密度；

——粒子i上一时刻的粒子数密度；

——粒子i上两个时刻的粒子数密度；

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

γ——调节系数，常取0.5；

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

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

获得的方程组形式为AP＝b，其中A为线性方程组的系数矩阵，其为对称的稀疏矩阵，P为未知的压力向量，b为压力泊松方程右侧的源项，随后采用不同的线性方程求解器对压力泊松方程进行求解，采用共轭梯度求解器，计算得到每个粒子的压力值；

步骤11：粒子的速度和位置计算，根据动量守恒方程，如公式(17)，显式计算重力项、粘性项和湍流应力项，得到速度和位置的估算值，如公式(34)和公式(35)所示；

式中

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

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

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

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

——粒子i的湍流应力矢量N；

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

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

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

随后采用速度和位置的估算值显式计算压力，对速度和位置进行修正，如公式(36)、公式(37)和公式(38)所示；

式中

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

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

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

P_j——粒子j的压力Pa；

——粒子i周围粒子作用半径内所有粒子压力的最小值Pa；

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

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

通过步骤8至步骤11，计算得到每个粒子的速度和位置；

综上，通过步骤1和步骤2完成核反应堆芯材料熔池的几何建模和物性初始化；通过步骤4和步骤5，模拟熔池内的混合过程和化学反应过程，计算粒子内的物质变化、质量变化和能量变化；通过步骤6，模拟熔池内的换热过程，考虑对流换热、导热和辐射换热，再结合步骤9中的湍流热流项，计算获得粒子下一时刻的焓值；通过步骤7计算每个粒子的温度和固相率，根据固相率判断粒子的相态；通过步骤8至步骤11，计算每个粒子的速度和位置，其中主要考虑在重力、压力、粘性力和湍流应力的作用下，熔池内部的流动情况；通过以上步骤，模拟了核反应堆堆芯材料熔池的动态变化过程，得到了熔池内熔融物的位置、速度、物质种类、压力、温度、焓值、相态随时间的变化，通过以上数据对熔池内的物质混合、化学反应和湍流流动过程展开机理性分析。

本发明方法为核反应堆堆芯材料熔池演变特性分析提供解决方案，为核电厂反应堆严重事故安全特性的研究提供重要依据。

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

本发明的基于粒子法的核反应堆堆芯材料熔池演变特性分析方法，综合考虑了核反应堆堆芯材料熔池中可能出现的所有现象，包括物质混合、化学反应、传热相变、流体运动和湍流流动，能够对核反应堆堆芯材料熔池演变特性进行机理性分析。并且该方法基于粒子法，能够捕捉熔池内的物质变化、相态变化和流动特性，并且具备方便建模、展示性高和实时追踪动态粒子的优点。综上，该方法能够更加全面、有效、高效地对核反应堆堆芯材料熔池演变特性进行分析。

附图说明

图1是本发明基于粒子法的核反应堆堆芯材料熔池演变特性分析方法的流程图。

具体实施方式

本发明基于粒子法的核反应堆堆芯材料熔池演变特性分析方法，如图1所示，步骤如下：

步骤1：几何模型建模，采用离散的粒子构建几何模型，每个粒子携带有位置和速度信息，当粒子直径小于实际几何的一百分之一时，粒子构型能够精确还原真实几何的结构和运动特性；采用0、1、2号粒子表征液相、固相、固液混合相，采用3号粒子表征虚拟壁面粒子，虚拟壁面粒子保证壁面边界处的粒子数密度守恒，防止流体粒子向壁面的穿透；每个粒子携带有物质信息，核反应堆堆芯材料物质种类繁多，无法精确考虑所有物质，因此只考虑含量较高的物质：二氧化铀、不锈钢、锆和二氧化锆；此外，当压力容器失效后，熔融物进入安全壳，落入堆坑，在堆坑处重新形成熔池，该过程涉及熔融物与混凝土的相互作用，熔池内会混入混凝土材料，主要包括：二氧化硅和氧化钙；每个粒子携带有焓值信息和温度信息；

步骤2：粒子物性计算，设置每种物质的物性参数信息，根据每个粒子的物质信息，计算每个粒子的物性参数，包括密度、比热、熔点或固相线温度和液相线温度、沸点；

步骤3：定义粒子核函数，表征粒子间的相互作用程度，采用三阶尖峰核函数，如公式(1)所示，

式中

w——粒子核函数；

r——粒子间距m；

r_e——粒子作用半径m；

当粒子间距大于粒子作用半径时，忽略粒子间的微小作用力，认为粒子间不存在作用力；定义粒子数密度，表征粒子的密集程度，

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

式中

n——粒子数密度；

i——粒子ID；

j——粒子ID；

w——粒子核函数；

r——粒子作用半径m；

步骤4：物质混合计算，在核反应堆堆芯材料熔池内，物质种类繁多，不同物质与物质之间会相互混合，如氧化物之间的混合、金属之间的混合，有些物质之间是互不相溶的，如金属与氧化物之间，物质的互不相溶会使得熔池出现分层现象，物质的混合过程，决定了熔池结构特性，熔池的不同结构，其换热及流动特性会存在较大差异；物质的混合基于扩散定律，考虑两种物质之间的混合过程，如公式(3)所示，

式中

——粒子i下一时刻的物质浓度kg/m³；

——粒子i当前时刻的物质浓度kg/m³；

D——物质混合系数；

Δt——时间步长s；

d——维度；

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

——粒子j当前时刻的物质浓度kg/m³；

w_ij——粒子i和粒子j的核函数值；

r_ij——粒子i和粒子j的粒子间距；

物质的传递过程伴随着能量的传递，物质混合过程中的能量守恒如公式(4)和公式(5)所示，此处的物质混合过程不考虑化学反应；

式中

——粒子i下一时刻的体积焓J/m³；

——粒子i当前时刻的体积焓J/m³；

D——物质混合系数；

Δt——时间步长s；

d——维度；

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

w_ij——粒子i和粒子j的核函数值；

r_ij——粒子i和粒子j的粒子间距；

——粒子j当前时刻的物质浓度kg/m³；

——粒子i当前时刻的物质浓度kg/m³；

——物质x的焓J/kg；

T_x——物质x的温度K；

c_p,x——物质x的定压比热容J/(kg·K)；

T_s——固相线温度K；

Δh_sl,x——物质x的熔化潜热J/kg；

T_l——液相线温度K；

粒子的物性参数由各物质的物质浓度决定，其中粒子密度如公式(6)所示，粒子的比热、熔化潜热、汽化潜热、热导率、粘性如公式(7)所示，固相线温度和液相线温度根据不同的混合物的物性模型计算得到；

式中

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

x——物质种类号；

m_x——粒子中物质x的物质浓度kg/m³；

Y——物性参数；

步骤5：在堆芯材料熔池内，存在大量的金属，极易被氧化，主要考虑铁和锆的氧化反应：

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

基于以上化学反应方程式，当两个粒子间距足够小或粒子满足化学反应条件时，根据化学反应速率进行物质的转变，如公式(8)和公式(9)所示；同时以内热源的形式改变粒子的焓值来模拟化学热的过程，如公式(10)所示；物质种类改变后，同时对粒子物性进行更新计算；

式中

——粒子i下一时刻物质x的物质浓度kg/m³；

——粒子i当前时刻物质x的物质浓度kg/m³；

C_x,chem——物质x的变化速率kg/(m³·s)；

Δt——时间步长s；

——粒子i下一时刻物质y的物质浓度kg/m³；

——粒子i当前时刻物质y的物质浓度kg/m³；

X——化学反应中物质x与物质y的物质变化比例；

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

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

Q_chem——化学反应热W/kg；

通过步骤4和步骤5，更新每个粒子的物质变化和焓值变化，并根据物质含量和焓值更新粒子的物性；

步骤6：传热计算，除了物质混合过程会伴随着能量传递外，粒子会通过对流换热、导热和辐射换热进行能量交换；在粒子法中，通过粒子的运动和导热过程的结合，还原对流换热过程，对流换热和导热过程采用公式(11)计算，辐射换热过程采用公式(12)和公式(13)计算；

式中

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

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

d——维度；

w_ij——粒子i和粒子j的核函数值；

r_ij——粒子i和粒子j的粒子间距m；

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

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

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

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

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

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

Δt——时间步长s；

h_r——辐射换热系数W/(m²·K)；

σ_stef——斯特藩-玻尔兹曼常数，5.67×10^-8W/(m²·K⁴)；

ε_i——自由表面粒子i的辐射发射率；

ε_env——周围环境的辐射发射率；

T_i——自由表面粒子i的温度K；

T_env——堆腔周围环境温度K；

Q_r——辐射换热量W/m³；

A——辐射换热面积m²；

l——粒子直径m；

通过步骤6，获得每个粒子下一时刻的物质焓值和温度；

步骤7：相变计算，通过能量守恒计算得到每个时刻下每个粒子的焓值，根据焓值可以计算得到粒子对应温度，如公式(14)所示，

式中

T——温度K；

T_s——固相线温度K；

h——焓值J/kg；

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

c_p——定压比热容J/(kg·K)；

T_l——液相线温度K；

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

定义固相率ε判断粒子的相态变化，固相率为焓值的函数：

式中

ε——固相率；

h——焓值J/kg；

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

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

当固相率为1时，粒子完全为固态，当固相率为0时，粒子完全为液态，当固相率介于0和1之间时，粒子为固液混合态；固液混合态较为特殊，固相和液相的占比对粒子的行为特性影响较大；定义临界固相率ε_crit来判断固液混合态粒子的行为特征，当ε＞ε_crit时，粒子内固相占比较高，粒子行为趋向于固态；当ε＜ε_crit时，粒子内液相占比较高，粒子行为趋向于液态；在本方法中，临界固相率设为0.55；

通过步骤7计算得到每个粒子的温度和相态；

步骤8：质量和动量守恒计算，在本方法中，针对不可压缩流体展开计算，在质量守恒方程中忽略密度变化，如公式(16)所示；动量守恒主要考虑压力、重力、粘性力和湍流应力，如公式(17)所示，

式中

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

t——时间s；

——速度矢量m/s；

P——压力Pa；

ν——运动粘度m²/s；

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

——湍流应力N；

步骤9：湍流计算，在核反应堆熔池内，由于温度梯度较大，物质组分间密度差较大，熔池内部会存在剧烈的搅浑和对流现象，为了能够准确模拟熔池内的流动特性，采用大涡模拟方法对熔池湍流特性展开计算；以粒子尺寸为区分，将流体的湍流运动分为大尺度运动和小尺度运动，对于涡旋尺度大于粒子尺寸的，采用纳维斯-托克斯方程直接求解，对于涡旋尺度小于粒子尺寸的，采用湍流应力项进行考虑，因此可以将运动分为时均值和脉动值，如公式(18)所示，

式中

φ——瞬时值；

——时均值；

φ′——脉动值；

对于动量守恒方程中，如公式(17)，脉动值即为湍流应力项，对于能量守恒方程中，脉动值为湍流热流项，如公式(19)所示，

式中

T——温度K；

t——时间s；

α——热扩散率m²/s；

——湍流热流项；

湍流应力项形式如公式(20)所示，为三阶矩阵，

式中

τ——湍流应力N；

x,y,z——表示湍流应力的不同方向；

每个方向的应力如公式(21)、公式(22)和公式(23)所示：

式中

τ_ab——ab方向上的湍流应力N；

u——速度m/s；

ν_t——涡粘性；

τ_γγ——三阶矩阵对角线上的切应力N；

δ_ab——颜色函数；

a——方向a；

b——方向b；

其中涡粘性采用Smagoringsky亚格子模型求解，如公式(24)、公式(25)、公式(26)、公式(27)、公式(28)和公式(29)所示，

式中

ν_t——涡粘性；

C_D——动态湍流系数；

Δ——滤波宽度；

Δ₂——二次滤波宽度；

——S的滤波值；

——u的滤波值；

——某参数滤波值；

r_ij——粒子i和粒子j的粒子间距m；

G——高斯核函数；

公式(22)中速度的偏导项根据梯度模型展开：

式中

u——速度m/s；

a——方向a；

s——长度m；

b——方向b；

d——维度；

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

r_ij——粒子i和粒子j的粒子间距m；

r_b,j——粒子j在方向b上的距离m；

r_b,i——粒子i在方向b上的距离m；

w_ij——粒子i和粒子j的核函数值；

计算获得每个方向上的湍流应力后，湍流应力的散度项根据散度模型展开：

式中

——方向a上的切应力向量；

d——维度；

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

——粒子j在方向a上的切应力向量；

——粒子i在方向a上的切应力向量；

——粒子j的位置矢量；

——粒子i的位置矢量；

r_ij——粒子i和粒子j的粒子间距m；

w_ij——粒子i和粒子j的核函数值；

公式(19)中的湍流热流项计算如下：

式中

θ_a——方向a上的湍流热流；

u——速度m/s；

T——温度K；

ν_t——涡粘性；

Pr_t——普朗特数；

a——方向a；

通过以上计算获得了湍流应力项和湍流热流项，从而可以计算得到每个粒子在考虑湍流影响的速度、位置和温度；

步骤10：压力泊松方程求解，针对公式(17)中的压力梯度项，采用隐式迭代求解压力泊松方程的方式计算得到每个粒子的位置和速度，粒子格式离散后的压力泊松方程如公式(33)所示，

式中

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

P——压力Pa；

t——时间s；

β——调节系数，常取0.4；

——粒子i当前时刻的粒子数密度；

——粒子i上一时刻的粒子数密度；

——粒子i上两个时刻的粒子数密度；

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

γ——调节系数，常取0.5；

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

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

获得的方程组形式为AP＝b，其中A为线性方程组的系数矩阵，其为对称的稀疏矩阵，P为未知的压力向量，b为压力泊松方程右侧的源项，随后可采用不同的线性方程求解器对压力泊松方程进行求解，本方法中采用共轭梯度求解器，计算得到每个粒子的压力值；

步骤11：粒子的速度和位置计算，根据动量守恒方程，如公式(17)，显式计算重力项、粘性项和湍流应力项，得到速度和位置的估算值，如公式(34)和公式(35)所示；

式中

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

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

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

——速度矢量m/s；

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

——粒子i的湍流应力矢量N；

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

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

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

Δt——时间步长s；

随后采用速度和位置的估算值显式计算压力，对速度和位置进行修正，如公式(36)、公式(37)和公式(38)所示；

式中

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

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

Δt——时间步长s；

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

P——压力Pa；

d——维度；

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

P_j——粒子j的压力Pa；

——粒子i周围粒子作用半径内所有粒子压力的最小值Pa；

r_ij——粒子i和粒子j的粒子间距m；

——粒子j的位置矢量；

——粒子i的位置矢量；

w_ij——粒子i和粒子j的核函数值；

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

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

通过步骤8至步骤11，计算得到每个粒子的速度和位置；

综上，通过步骤1和步骤2完成核反应堆芯材料熔池的几何建模和物性初始化；通过步骤4和步骤5，模拟熔池内的混合过程和化学反应过程，计算粒子内的物质变化、质量变化和能量变化；通过步骤6，模拟熔池内的换热过程，考虑对流换热、导热和辐射换热，再结合步骤9中的湍流热流项，计算获得粒子下一时刻的焓值；通过步骤7计算每个粒子的温度和固相率，根据固相率判断粒子的相态；通过步骤8至步骤11，计算每个粒子的速度和位置，其中主要考虑在重力、压力、粘性力和湍流应力的作用下，熔池内部的流动情况；通过以上步骤，模拟了核反应堆堆芯材料熔池的动态变化过程，得到了熔池内熔融物的位置、速度、物质种类、压力、温度、焓值、相态随时间的变化，通过以上数据对熔池内的物质混合、化学反应和湍流流动过程展开机理性分析。

Claims (1)

1.一种基于粒子法的核反应堆堆芯材料熔池演变特性分析方法，其特征在于：步骤如下：

步骤1：几何模型建模，采用离散的粒子构建核反应堆芯材料熔池的几何模型，每个粒子携带有位置和速度信息，当粒子直径小于实际几何的一百分之一时，粒子构型能够精确还原真实几何的结构和运动特性；采用0、1、2号粒子表征液相、固相、固液混合相，采用3号粒子表征虚拟壁面粒子，虚拟壁面粒子保证壁面边界处的粒子数密度守恒，防止流体粒子向壁面的穿透；每个粒子携带有物质信息，核反应堆堆芯材料物质种类繁多，无法精确考虑所有物质，因此只考虑含量较高的物质：二氧化铀、不锈钢、锆和二氧化锆；此外，当压力容器失效后，熔融物进入安全壳，落入堆坑，在堆坑处重新形成熔池，该过程涉及熔融物与混凝土的相互作用，熔池内会混入混凝土材料，主要包括：二氧化硅和氧化钙；每个粒子携带有焓值信息和温度信息；

步骤2：粒子物性计算，设置每种物质的物性参数信息，根据每个粒子的物质信息，计算每个粒子的物性参数，包括密度、比热、熔点或固相线温度和液相线温度、沸点；

步骤3：定义粒子核函数，表征粒子间的相互作用程度，采用三阶尖峰核函数，如公式(1)所示，

式中

w——粒子核函数；

r——粒子间距m；

r_e——粒子作用半径m；

当粒子间距大于粒子作用半径时，忽略粒子间的微小作用力，认为粒子间不存在作用力；

定义粒子数密度，表征粒子的密集程度，

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

式中

n——粒子数密度；

i——粒子ID；

j——粒子ID；

w——粒子核函数；

步骤4：物质混合计算，在核反应堆堆芯材料熔池内，物质种类繁多，不同物质与物质之间会相互混合，有些物质之间是互不相溶的，物质的互不相溶会使得熔池出现分层现象，物质的混合过程，决定了熔池结构特性，熔池的不同结构，其换热及流动特性会存在较大差异；物质的混合基于扩散定律，考虑两种物质之间的混合过程，如公式(3)所示，

式中

——粒子i下一时刻的物质浓度kg/m³；

——粒子i当前时刻的物质浓度kg/m³；

D——物质混合系数；

Δt——时间步长s；

d——维度；

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

——粒子j当前时刻的物质浓度kg/m³；

w_ij——粒子i和粒子j的核函数值；

r_ij——粒子i和粒子j的粒子间距；

物质的传递过程伴随着能量的传递，物质混合过程中的能量守恒如公式(4)和公式(5)所示，此处的物质混合过程不考虑化学反应；

式中

——粒子i下一时刻的体积焓J/m³；

——粒子i当前时刻的体积焓J/m³；

——物质x的焓J/kg；

T_x——物质x的温度K；

c_p,x——物质x的定压比热容J/(kg·K)；

T_s——固相线温度K；

Δh_sl,x——物质x的熔化潜热J/kg；

T_l——液相线温度K；

粒子的物性参数由各物质的物质浓度决定，其中粒子密度如公式(6)所示，粒子的比热、熔化潜热、汽化潜热、热导率、粘性如公式(7)所示，固相线温度和液相线温度根据不同的混合物的物性模型计算得到；

式中

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

x——物质种类号；

m_x——粒子中物质x的物质浓度kg/m³；

Y——物性参数；

步骤5：在堆芯材料熔池内，存在大量的金属，极易被氧化，主要考虑铁和锆的氧化反应：

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

基于以上化学反应方程式，当两个粒子间距足够小或粒子满足化学反应条件时，根据化学反应速率进行物质的转变，如公式(8)和公式(9)所示；同时以内热源的形式改变粒子的焓值来模拟化学热的过程，如公式(10)所示；物质种类改变后，同时对粒子物性进行更新计算；

式中

——粒子i下一时刻物质x的物质浓度kg/m³；

——粒子i当前时刻物质x的物质浓度kg/m³；

C_x,chem——物质x的变化速率kg/(m³·s)；

——粒子i下一时刻物质y的物质浓度kg/m³；

——粒子i当前时刻物质y的物质浓度kg/m³；

X——化学反应中物质x与物质y的物质变化比例；

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

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

Q_chem——化学反应热W/kg；

通过步骤4和步骤5，更新每个粒子的物质变化和焓值变化，并根据物质含量和焓值更新粒子的物性；

步骤6：传热计算，除了物质混合过程会伴随着能量传递外，粒子会通过对流换热、导热和辐射换热进行能量交换；在粒子法中，通过粒子的运动和导热过程的结合，还原对流换热过程，对流换热和导热过程采用公式(11)计算，辐射换热过程采用公式(12)和公式(13)计算；

式中

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

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

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

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

Δt——时间步长s；

h_r——辐射换热系数W/(m²·K)；

σ_stef——斯特藩-玻尔兹曼常数，5.67×10^-8W/(m²·K⁴)；

ε_i——自由表面粒子i的辐射发射率；

ε_env——周围环境的辐射发射率；

T_i——自由表面粒子i的温度K；

T_env——堆腔周围环境温度K；

Q_r——辐射换热量W/m³；

A——辐射换热面积m²；

l——粒子直径m；

通过步骤6，获得每个粒子下一时刻的物质焓值和温度；

步骤7：相变计算，通过能量守恒计算得到每个时刻下每个粒子的焓值，根据焓值计算得到粒子对应温度，如公式(14)所示，

式中

T——温度K；

h——焓值J/kg；

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

c_p——定压比热容J/(kg·K)；

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

定义固相率ε判断粒子的相态变化，固相率为焓值的函数：

式中

ε——固相率；

当固相率为1时，粒子完全为固态，当固相率为0时，粒子完全为液态，当固相率介于0和1之间时，粒子为固液混合态；固液混合态较为特殊，固相和液相的占比对粒子的行为特性影响较大；定义临界固相率ε_crit来判断固液混合态粒子的行为特征，当ε＞ε_crit时，粒子内固相占比较高，粒子行为趋向于固态；当ε＜ε_crit时，粒子内液相占比较高，粒子行为趋向于液态；临界固相率设为0.55；

通过步骤7计算得到每个粒子的温度和相态；

步骤8：质量和动量守恒计算，针对不可压缩流体展开计算，在质量守恒方程中忽略密度变化，如公式(16)所示；动量守恒主要考虑压力、重力、粘性力和湍流应力，如公式(17)所示，

式中

t——时间s；

——速度矢量m/s；

P——压力Pa；

ν——运动粘度m²/s；

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

——湍流应力N；

步骤9：湍流计算，在核反应堆熔池内，由于温度梯度较大，物质组分间密度差较大，熔池内部会存在剧烈的搅浑和对流现象，为了能够准确模拟熔池内的流动特性，采用大涡模拟方法对熔池湍流特性展开计算；以粒子尺寸为区分，将流体的湍流运动分为大尺度运动和小尺度运动，对于涡旋尺度大于粒子尺寸的，采用纳维斯-托克斯方程直接求解，对于涡旋尺度小于粒子尺寸的，采用湍流应力项进行考虑，因此将运动分为时均值和脉动值，如公式(18)所示，

式中

φ——瞬时值；

——时均值；

φ′——脉动值；

对于动量守恒方程中，如公式(17)，脉动值即为湍流应力项，对于能量守恒方程中，脉动值为湍流热流项，如公式(19)所示，

式中

T——温度K；

α——热扩散率m²/s；

——湍流热流项；

湍流应力项形式如公式(20)所示，为三阶矩阵，

式中

τ——湍流应力N；

x,y,z——表示湍流应力的不同方向；

每个方向的应力如公式(21)、公式(22)和公式(23)所示：

式中

τ_ab——ab方向上的湍流应力N；

u——速度m/s；

ν_t——涡粘性；

τ_γγ——三阶矩阵对角线上的切应力N；

δ_ab——颜色函数；

a——方向a；

b——方向b；

其中涡粘性采用Smagoringsky亚格子模型求解，如公式(24)、公式(25)、公式(26)、公式(27)、公式(28)和公式(29)所示，

式中

ν_t——涡粘性；

C_D——动态湍流系数；

Δ——滤波宽度；

Δ₂——二次滤波宽度；

——S的滤波值；

——u的滤波值；

——某参数滤波值；

r_ij——粒子i和粒子j的粒子间距m；

G——高斯核函数；

公式(22)中速度的偏导项根据梯度模型展开：

式中

u——速度m/s；

a——方向a；

s——长度m；

r_b,j——粒子j在方向b上的距离m；

r_b,i——粒子i在方向b上的距离m；

计算获得每个方向上的湍流应力后，湍流应力的散度项根据散度模型展开：

式中

——方向a上的切应力向量；

——粒子j在方向a上的切应力向量；

——粒子i在方向a上的切应力向量；

——粒子j的位置矢量；

——粒子i的位置矢量；

公式(19)中的湍流热流项计算如下：

式中

θ_a——方向a上的湍流热流；

Pr_t——普朗特数；

通过以上计算获得了湍流应力项和湍流热流项，从而能够计算得到每个粒子在考虑湍流影响的速度、位置和温度；

步骤10：压力泊松方程求解，针对公式(17)中的压力梯度项，采用隐式迭代求解压力泊松方程的方式计算得到每个粒子的位置和速度，粒子格式离散后的压力泊松方程如公式(33)所示，

式中

β——调节系数，常取0.4；

——粒子i当前时刻的粒子数密度；

——粒子i上一时刻的粒子数密度；

——粒子i上两个时刻的粒子数密度；

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

γ——调节系数，常取0.5；

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

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

获得的方程组形式为AP＝b，其中A为线性方程组的系数矩阵，其为对称的稀疏矩阵，P为未知的压力向量，b为压力泊松方程右侧的源项，随后采用不同的线性方程求解器对压力泊松方程进行求解，采用共轭梯度求解器，计算得到每个粒子的压力值；

步骤11：粒子的速度和位置计算，根据动量守恒方程，如公式(17)，显式计算重力项、粘性项和湍流应力项，得到速度和位置的估算值，如公式(34)和公式(35)所示；

式中

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

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

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

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

——粒子i的湍流应力矢量N；

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

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

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

随后采用速度和位置的估算值显式计算压力，对速度和位置进行修正，如公式(36)、公式(37)和公式(38)所示；

式中

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

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

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

P_j——粒子j的压力Pa；

——粒子i周围粒子作用半径内所有粒子压力的最小值Pa；

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

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

通过步骤8至步骤11，计算得到每个粒子的速度和位置；

综上，通过步骤1和步骤2完成核反应堆芯材料熔池的几何建模和物性初始化；通过步骤4和步骤5，模拟熔池内的混合过程和化学反应过程，计算粒子内的物质变化、质量变化和能量变化；通过步骤6，模拟熔池内的换热过程，考虑对流换热、导热和辐射换热，再结合步骤9中的湍流热流项，计算获得粒子下一时刻的焓值；通过步骤7计算每个粒子的温度和固相率，根据固相率判断粒子的相态；通过步骤8至步骤11，计算每个粒子的速度和位置，其中主要考虑在重力、压力、粘性力和湍流应力的作用下，熔池内部的流动情况；通过以上步骤，模拟了核反应堆堆芯材料熔池的动态变化过程，得到了熔池内熔融物的位置、速度、物质种类、压力、温度、焓值、相态随时间的变化，通过以上数据对熔池内的物质混合、化学反应和湍流流动过程展开机理性分析。

Images