[发明专利]利用粒子网格混合法研究堆芯内共晶反应及高温熔化行为的方法

权利要求书

1.一种利用粒子网格混合法研究堆芯内共晶反应及高温熔化行为的方法，其特征在于：步骤如下：

步骤1：对反应堆堆芯进行粒子建模，不同的物质使用不同种类的粒子表示，粒子初始状态包括粒子的位置、速度、压力和温度，对应标记为r_i、u_i、p_i和T_i，粒子初始布置时相邻两粒子间的距离记为l₀，每个粒子内按所含元素标记对应元素的质量，如粒子i内包含A、B两种元素，则其质量对应标记为m_iA、m_iB，粒子i的质量m_i使用公式(1)计算，

m_i＝m_iA+m_iB 公式(1)

粒子i的密度ρ_i使用公式(2)计算，

对于熔化焓等则通过公式(3)求得，

步骤2：根据各粒子的位置划分网格：在本步骤中，对于每一个粒子i，设置粒子作用域半径l_r，若粒子i的坐标位置为(I_x，I_y，I_z)，设定粒子i周围的某一粒子为粒子j，坐标位置为(J_x，J_y，J_z)，粒子i与粒子j的距离为l_ij，使用公式(4)计算，

若满足l_ij≤l_r，则将该就j粒子标记为i粒子的有效作用粒子；将所有满足条件的m个j粒子进行汇总并标记为j₁，j₂，……，j_m，然后对其中每3个粒子进行组合，按照排列组合法，将有种组合，由公式(5)进行计算，

对每个组合中的粒子进行考察，取其中一个组合的3个粒子，这3个粒子必须满足能且只能构成一个平面的要求，将取出的3个粒子标记为j₁，j₂，j₃，粒子中心分别标记为J₁，J₂，J₃，空间坐标对应为(J_x1，J_y1，J_z1)，(J_x2，J_y2，J_z2)，(J_x3，J_y3，J_z3)，连接J₁、J₂，连接J₁、J₃，连接J₂、J₃，J₁、J₂、J₃三点构成平面的其中一个充要条件是直线J₁J₂与直线J₁J₃之间存在一个不为0的夹角，则只需满足公式(6)，

式中，||J₁J₂||、||J₁J₃||、||J₂J₃||分别为线段J₁J₂、J₁J₃、J₂J₃的长度，使用公式(7)进行计算，

对每一个i粒子，针对其种组合进行判断，若组合内的3个粒子能够满足构成唯一平面的要求则留下该组合，随后针对留下的组合，考察i粒子与组合内3个粒子的关系，i粒子与组合中三个粒子必须能够构成四面体；对于留下的每个组合而言，这3个粒子已经满足构成唯一平面，那么要使这i粒子与组合中的3个粒子能够构成四面体，只需满足i粒子不在组合中的3个粒子所构成的平面上即可，添加i粒子，粒子中心标记为I，空间坐标定义为(I_x，I_y，I_z)，若要使I、J₁、J₂、J₃四个点构成四面体，则只需要满足i粒子到平面J₁J₂J₃的距离为正数，使用公式(8)进行判定，

式中：

——i粒子到平面J₁J₂J₃的距离m；

A、B、C、D——平面J₁J₂J₃的平面方程的四个系数，通过公式(9)获得，

i粒子与j₁、j₂、j₃粒子构成的平面J₁J₂J₃的距离若满足公式(8)，则说明i粒子能够与组合内3个粒子构成的三角形J₁J₂J₃构成四面体IJ₁J₂J₃，其四个面分别为三角形IJ₁J₂、三角形IJ₁J₃、三角形IJ₂J₃和三角形J₁J₂J₃，面积对应标记为三角形的面积使用如公式(10)所示的海伦公式计算

式中，

S——三角形的面积；

l_a、l_b、l_c——三角形的三边长；

p——三角形的半周长，即p＝(l_a+l_b+l_c)/2；

四面体IJ₁J₂J₃的体积通过公式(11)进行计算，

四面体IJ₁J₂J₃中四个顶点I、J₁、J₂、J₃分别对四面体总体积的分属控制体积，对应标记为V_I、这4个体积有公式(12)所示的关系，

各控制体积由公式(13)进行计算，

除此之外，四面体IJ₁J₂J₃的四个点的分属控制体积与其相对的三角形面积是成反比的，即有公式(14)所示的关系，

由以上公式联立求解即求得四面体IJ₁J₂J₃的四个点的分属控制体积；

然后将对四面体IJ₁J₂J₃中四个面的三角形进行面积划分，以三角形J₁J₂J₃为例，在三角形J₁J₂J₃中，三角形J₁J₂J₃的面积通过公式(10)海伦公式求出；点N₁₂为线段J₁J₂的中点，点N₁₃为线段J₁J₃的中点，点N₂₃为线段J₂J₃的中点；此时需要在三角形中寻找P点，对于三角形J₁J₂J₃而言，将其标记为P_x，连接P_x、J₁，连接P_x、J₂，连接P_x、J₃，连接P_x、N₁₂，连接P_x、N₁₃，连接P_x、N₂₃，可以将三角形J₁J₂J₃划分成6个小三角形，分别标记为三角形J₁P_xN₁₂，三角形J₁P_xN₁₃，三角形J₂P_xN₁₂，三角形J₂P_xN₂₃，三角形J₃P_xN₁₃，三角形J₃P_xN₂₃；因为点N₁₃为线段J₁J₃的中点，所以三角形J₁P_xN₁₃与三角形J₃P_xN₁₃面积相等，标记为同样的，有对于各个三角形的面积，有如公式(15)和公式(16)所示的几何关系，

三角形的三个角的角度由类似于公式(6)所示的余弦定理求出；由公式(15)和公式(16)求出三角形J₁J₂J₃所分成的6个小三角形的面积；同理也能求得四面体IJ₁J₂J₃另外三个面所分成的共计18个小三角形的面积；

在四面体IJ₁J₂J₃中寻找一点Q_x，对于顶点I的控制体积而言，连接Q_xP_x1，Q_xP_x2，Q_xP_x3，Q_xN₁，Q_xN₂，Q_xN₃，Q_xI，其中P_x1、P_x2、P_x3、分别对应为三角形IJ₁J₂、三角形IJ₂J₃、三角形IJ₁J₃所寻找的P点，N₁、N₂、N₃则分别为线段IJ₁、IJ₂、IJ₃的中点，至此点I的控制体积分成了六个四面体，分别为四面体IP_x1N₁Q_x，四面体IP_x1N₂Q_x，四面体IP_x2N₂Q_x，四面体IP_x2N₃Q_x，四面体IP_x3N₁Q_x和四面体IP_x3N₃Q_x，其体积对应标记为和对于四面体的体积，有公式(17)所示的关系，

设点Q_x距离平面J₁J₂J₃、平面IJ₁J₂，平面IJ₂J₃，平面IJ₁J₃的距离分别为H₀、H₁、H₂、H₃，对于顶点I的控制体积可列出如公式(18)的方程，

同理对于顶点J₁、J₂、J₃也列出相似的方程，联立这4个方程及公式(13)能够求出H₀、H₁、H₂、H₃的值；

至此，在I、J₁、J₂、J₃组成的四面体中，J₁对I起影响的点I的控制体积为记为J₂对I起影响的点I的控制体积为记为J₃对I起影响的点I的控制体积为记为所以，在I、J₁、J₂、J₃组成的四面体系统中，点J₁对I的有效传输面积标记为点J₂对I的有效传输面积标记为点J₁对I的有效传输面积标记为通过公式(19)计算，

式中，

——I点到J₁点的距离；

——I点到J₂点的距离；

——I点到J₃点的距离；

以上的划分网格的方法描述是针对其中一个四面体系统而言的，在计算中需要对i粒子的所有符合的四面体网格进行计算得到有效传输面积数据之后才能开始进行质量扩散计算及温度计算；

步骤3：在质量扩散的计算中，使用的传质方程为菲克第二定律，方程如公式(20)所示，对于粒子i内的A、B元素分别列出方程，

式中，

m——粒子内的某元素质量kg；

t——时间s；

D_diff——粒子内的某元素的质量扩散系数m²/s；

公式(21)给出质量扩散方程的Δt时间步长的积分形式，

式中，

m——粒子内的某元素质量kg；

t——时间s；

Δt——时间步长s；

V——体积m³；

CV——控制体积m³；

D_diff——粒子内的某元素的质量扩散系数m²/s；

对于粒子i而言，其已经通过步骤2划分为若干个四面体网格，所以公式(21)能够离散成公式(22)的形式，

公式(22)等式的左边对i粒子划分成的多个四面体的某一元素的质量变化进行加和，等式的右边则是i粒子的该元素的质量变化率，计算中不考虑粒子的体积变化；式中，

S_j→i——步骤2中网格划分完成后计算的j粒子对i粒子的有效传输面积；

D_diff——某元素的质量扩散系数m²/s；

r_j——j粒子的坐标矢量；

r_i——i粒子的坐标矢量；

——j粒子中x元素在n时层的质量kg；

——i粒子中x元素在n时层的质量kg；

——i粒子中x元素在n+1时层的质量，即在n时层的时间步长Δt后的x元素质量kg；

i粒子在n+1时层的质量使用公式(23)计算，

在温度计算中，能量守恒方程如公式(24)所示，

式中，

T——温度K；

t——时间s

κ——导热系数W/(m²·K)；

ρ——密度kg/m³；

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

S_T——温度的源项；

公式(25)给出能量方程的Δt时间步长的积分形式，

式中，

t——时间s；

Δt——时间步长s；

T——温度K；

V——体积m³；

CV——控制体积m³；

κ——导热系数W/(m²·K)；

ρ——密度kg/m³；

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

S_T——温度的源项；

对于粒子i而言，其已经通过步骤2划分为若干个四面体网格，所以公式(25)能够离散成公式(26)的形式，

公式(26)等式的左边对i粒子划分成的多个四面体的能量变化进行加和，等式的右边则是i粒子的能量变化率；式中：

S_j→i——步骤2中网格划分完成后计算的j粒子对i粒子的有效传输面积；

κ——导热系数W/(m²·K)；

r_j——j粒子的坐标矢量；

r_i——i粒子的坐标矢量；

——j粒子在n时层的温度K；

T_iⁿ——i粒子在n时层的温度K；

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

c_pi——i粒子的比定压热容J/(kg·K)；

T_iⁿ⁺¹——i粒子在n+1时层的温度，即在n时层的时间步长Δt后的温度K；

求解公式(26)即得到i粒子在时间步长Δt之后的温度T_iⁿ⁺¹；

至此，i粒子在时间步长Δt之后的温度T_iⁿ⁺¹及其质量已经求解得到，并且对于i粒子内的各元素进行的单独质量求解极易得到各元素在粒子内的质量百分比，若粒子内只含有A、B两种元素，i粒子中A元素的质量百分比由公式(27)计算，

式中，

C_iA——i粒子中A元素的质量百分比；

m_iA——i粒子中A元素的质量kg；

m_iB——i粒子中B元素的质量kg；

此时使用A、B元素的二元相图进行粒子状态判定，在A-B二元相图中，横坐标为A元素的质量百分比，纵坐标为粒子的温度，曲线A_pQ_SB_p为固相线，曲线A_pQ_LB_p为液相线，粒子在多边形O₁O₂B_pQ_SA_p围成的区域内状态为固态，在多边形A_pQ_LB_pO₃O₄围成的区域内状态为液态，在多边形A_pQ_SB_pQ_L围成的区域内状态为固液两相混合状态；若在某一时刻，i粒子的温度为T_k，i粒子中A元素的质量百分比为C_k，在A-B二元相图中作温度T_k的等温线，做质量百分比C_k的等值线，交于点Q_p，此时点Q_p位于A-B二元相图的液相区域，则对应的，粒子相态为液相；对位于固液两相混合状态区域内的粒子，需要根据相图中的杠杆定律对其成分进行分析；同样的方法，对于温度为T_px和粒子内A元素质量百分比为C_px的粒子k，判断其对应的Q_px点位于固液两相区内，作线段T_pxQ_px的延长线交固相线和液相线于Q_S和Q_L，再过Q_S和Q_L作垂直线交横坐标于C_S和C_L点，则粒子k中的固体和液体质量之间的关系使用公式(28)计算得到，

式中，

m_S——k粒子中处于固态的质量kg；

m_L——k粒子中处于液态的质量kg；

C_px——k粒子中A元素的质量百分比；

C_S——点C_S对应的横坐标数值；

C_L——点C_L对应的横坐标数值；

步骤4：使用移动粒子半隐式方法计算动量方程；动量方程如公式(29)所示

式中，

u——速度m/s；

t——时间s；

ρ——密度kg/m³；

p——压力Pa

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

f——外力加速度m²/s；

在此方程组中，动量方程中的粘性项及源项显式计算，其中粘性项的Laplace算子由公式(30)离散求解，

式中，

d——计算的维度；

n⁰——初始粒子数密度，通过计算，是核函数，使用公式(31)计算，只需将括号内的自变量修改为即可，为j粒子初始时刻的坐标矢量，为i粒子初始时刻的坐标矢量；

λ——是i粒子的邻居粒子和i粒子间距的方程，使用公式(32)计算，

式中，

r_e——粒子作用半径，位于半径内的粒子才会与i粒子发生相互作用，，r_e＝3.1l₀；

速度和位置的估算值由公式(33)进行计算，

式中，

——i粒子在上一时间步的速度矢量m/s；

r_iⁿ——i粒子在上一时间步的位置矢量m；

T_iⁿ——i粒子在上一时间步的温度K；

根据估算得到的粒子位置，通过公式(34)计算此时的粒子数密度n^*_i，只需要将公式中的|r_j-r_i|使用位置的估算值替代即可，由公式(34)所示的压力Poisson方程求解得到此时各粒子的压力值p_i；

公式(35)等式的左边使用如同公式(30)的形式离散，式中，

——i粒子在n+1时层的压力Pa；

α、γ——静态系数；

由计算得到的压力值，根据公式(36)计算速度的修正值；

式中，

u'_i——i粒子的速度修正值；

再由计算得到的速度修正值，根据公式(37)，公式(38)分别求解得到粒子在下一时层的速度及位置；

r_iⁿ⁺¹＝r_i^*+u'_iΔt 公式(38)

步骤5：根据具体问题的求解需要，对步骤2～步骤4进行时间步长Δt推进进行计算；

步骤6：输出计算所得到的堆芯中粒子的种类、速度位置压力及温度值T_iⁿ⁺¹，粒子质量及粒子内各元素质量，按照要求进行数据后处理，得到所要求的结果进行输出分析堆芯熔化进程的形态变化，堆芯内部热量分布、堆芯熔化情况、熔融物流动及分布的情况。