[发明专利]基于粒子法的核反应堆燃料早期变形分析方法

权利要求书

1.一种基于粒子法的核反应堆燃料早期变形分析方法，其特征在于：步骤如下：

步骤1：根据实际的核反应堆燃料结构建立粒子几何模型，定义材料数M区分不同的材料，M＝0,1,2分别表示核反应堆燃料元件的芯块、包壳和结构材料，粒子根据对应的材料携带对应的物性和力学参数；定义ID数区分固体边界，采用1号粒子模拟固体边界，2号粒子模拟固体内连续介质区域；不同固体之间采用参数Object区分，对于2号粒子，只有当Object数相同时才发生相互作用；初始化所有粒子的速度、位置、温度、应力张量和应变张量，分别为u_i、r_i、T_i、σ_i、ε_i；通过以上参数设定，建立燃料芯块、包壳与结构材料三者之间的固体边界，还原真实的燃料结构，并确定了燃料的初始位置、速度、温度、应力应变状态和相关物性参数；

步骤2：根据建立的初始粒子几何模型，为每个粒子建立固体内粒子邻居域，即存储每个粒子以自身为圆心，作用范围r_e为半径的区域内所有Object数与该粒子相同的粒子ID，固体内邻居域不随计算时间的推移而发生改变，即认为固体内连续介质发生相互作用的粒子始终保持不变；除了固体内粒子邻居域，还需要为所有1号粒子建立固体间粒子邻居域，即存储每个1号粒子以自身为圆心，作用范围r_e为半径的区域内所有Object数与该粒子不同的1号粒子ID，固体间邻居域会随着计算时间的推移而发生改变，即两个固体的相对位置会不断发生变化，当两者表面距离小于作用范围r_e时，两个固体的边界会发生相互作用，否则不会；定义粒子数密度为n，其值为粒子邻居域内所有粒子对中心粒子的核函数值之和，如公式(1)所示；

式中

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

w(r_ij,r_e)——核函数，采用的核函数形式为

r_ij——粒子i和粒子j的距离，m；

r_e——粒子作用范围，m；

下标i——任意粒子i；

下标j——粒子i的粒子邻居域内的任意粒子j；

由此，进一步定义固体内粒子数密度n_inner为粒子i的固体内邻居域的粒子数密度、粒子间粒子数密度n_outer为粒子i的固体间邻居域的粒子数密度、n⁰为参考粒子数密度，是假设粒子均匀分布情况下计算得到的粒子数密度；

步骤3：基于边界条件和固体碰撞模型计算固体边界粒子外力载荷，根据具体实际物理问题，为1号粒子设置外力的作用方向和大小，由公式(2)计算边界条件外力载荷的影响，边界条件外力载荷始终视为作用于粒子中心，因此不考虑该力所引起的粒子转动，公式采用张量表示法和下标表示法；

式中

v_α——α方向上的速度分量，m/s；

——α方向上的边界条件外力载荷引起的加速度，m/s²；

t——时间，s；

对于固体碰撞过程，采用固体碰撞模型计算1号粒子的碰撞载荷，如公式(3)所示，

式中

P_i^Contact——粒子i的接触应力，N/m²；

λ——拉梅常数的第一个参数，

(ε_γγ)_i——粒子i的应变张量对角项之和，即ε_γγ＝ε₁₁+ε₂₂+ε₃₃；

n_outer——粒子i的固体间邻居域的粒子数密度；

n_inner——粒子i的固体内邻居域的粒子数密度；

E——杨氏模量，N/m²；

υ——泊松比；

由碰撞载荷，计算对应的速度梯度，如公式(4)所示，

式中

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

t——时间，s；

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

d——空间维度；

n⁰——参考粒子数密度；

r_ij——粒子i和粒子j的相对位置矢量，m；

P_i^Contact——粒子i的接触应力，N/m²；

——粒子j的接触应力，N/m²；

——粒子i和粒子j的相对初始位置矢量，m；

——核函数，等同于

通过步骤3的计算，模拟核反应堆燃料早期变形过程中燃料元件受到的外力载荷和碰撞力，全面考虑了核反应堆燃料早期变形过程中可能发生的外力作用、燃料元件之间的碰撞、燃料元件与结构材料之间的碰撞、燃料包壳与芯块之间的碰撞、芯块与芯块之间的碰撞过程；计算得到每个粒子由外力载荷和碰撞引起的加速度，即燃料元件由外力载荷和碰撞引起的加速度；

步骤4：能量计算如公式(5)所示：

式中

h——焓值，J/kg；

t——时间，s；

ρ——密度，kg/m³；

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

T——温度，K；

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

Q_out——外热源，W/m³；

辐射换热采用斯忒藩-玻尔兹曼定律，如公式(6)所示

式中

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

ε——发射率；

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

T——温度，K；

T_env——环境温度K；

l₀——粒子直径m；

温度的拉普拉斯项采用拉普拉斯模型离散得到：

式中

d——空间维度；

n⁰——参考粒子数密度；

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

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

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

r_ij——粒子i和粒子j的相对位置矢量，m；

w(|r_ij|)——核函数，等同于w(|r_ij|,r_e)；

由焓值计算粒子温度：

式中

T——温度K；

T_s——熔点K；

h——焓值J/kg；

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

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

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

通过步骤4的计算，模拟核反应堆燃料早期变形过程中燃料元件内的导热、燃料元件间的辐射换热、芯块和包壳间的换热过程；计算得到每个粒子在不同时刻下的焓值和温度，即得到核反应堆燃料的焓值和温度随时间的变化过程；

步骤5：热膨胀应变计算如公式(8)所示，

式中

[ds^T]_i——粒子i的热膨胀应力增量张量，N/m²；

κ_i——粒子i的热膨胀系数；

T_i——粒子i的温度，K；

T_i^ref——粒子i的参考温度，K；

通过步骤5的计算，模拟核反应堆燃料早期变形过程中燃料元件的热膨胀变形过程；计算得到每个粒子在不同时刻下的热膨胀应变张量，即得到了燃料元件由于温度变化而发生碰撞变形的行为规律；

步骤6：计算总应变，固体内的相对变形是由当前相对位置减去旋转后的相对位置得到的，如公式(9)所示，

式中

u_ij——粒子i和粒子j的相对总变形矢量，m；

r_ij——粒子i和粒子j的相对位置矢量，m；

R——粒子i的旋转矩阵，R＝R_xR_yR_z；

——粒子i和粒子j的相对初始位置矢量，m；

R_x——粒子i绕x轴旋转的旋转矩阵，

R_y——粒子i绕y轴旋转的旋转矩阵，

R_z——粒子i绕z轴旋转的旋转矩阵，

θ_α,ij——粒子i和粒子j绕α轴旋转的相对角度，°，α可以取x，y，z；

由总变形可以计算得到总应变：

ε_ij——粒子i和粒子j的相对总应变矢量；

通过步骤6的计算，计算得到了每个粒子在不同时刻下的总应变张量，即得到了核反应堆燃料元件的总变形状态；

步骤7：基于弹性模型计算弹性应力，弹性模型如公式(11)所示，

式中

——弹性应力张量的分量，N/m²；

λ——拉梅常数的第一个参数，

μ——拉梅常数的第二个参数，

——弹性应变张量对角项之和，即

δ_αβ——克罗内克尔函数，

——弹性应变张量的分量；

其中记各向同性力：

式中

P_i^E——粒子i的弹性各向同性力，N/m²；

λ——拉梅常数的第一个参数，

——弹性应变张量对角项之和，即

d——空间维度；

n⁰——参考粒子数密度；

u_ij——粒子i和粒子j的相对总变形矢量，m；

r_ij——粒子i和粒子j的相对位置矢量，m；

——粒子i和粒子j的相对初始位置矢量，m；

——核函数，等同于

此外，在碰撞计算中，如果碰撞载荷P_i^Contact小于弹性各向同性力P_i^E，则令P_i^Contact＝P_i^E；

记各向异性力：

式中

——粒子i和粒子j间的弹性各向异性力，N/m²；

μ——拉梅常数的第二个参数，

——粒子i和粒子j的相对弹性应变矢量；

由弹性各向同性力和弹性各向异性力分别求得对应的速度梯度，如公式(14)和(15)所示，

式中

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

t——时间，s；

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

d——空间维度；

n⁰——参考粒子数密度；

r_ij——粒子i和粒子j的相对位置矢量，m；

P_i^E——粒子i的弹性各向同性力，N/m²；

——粒子j的弹性各向同性力，N/m²；

——粒子i和粒子j间的弹性各向异性力，N/m²；

——粒子i和粒子j的相对初始位置矢量，m；

——核函数，等同于

通过步骤7，模拟了核反应堆燃料早期变形过程中的弹性变形；计算得到每个粒子在不同时刻下的弹性应变张量、弹性应力张量、弹性力所引起的加速度，即得到燃料元件的弹性应力应变随时间的变化过程；

步骤8：基于角动量定理计算粒子旋转角速度，粒子j对粒子i的切应力为：

式中

F_ij——粒子j对粒子i的切应力，N；

d——空间维度；

n⁰——参考粒子数密度；

m_i——粒子i的质量，kg；

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

——粒子i和粒子j间的弹性各向异性力在与r_ij垂直方向上的投影，N/m²；

r_ij——粒子i和粒子j的相对位置矢量，m；

——粒子i和粒子j的相对初始位置矢量，m；

——核函数，等同于

粒子j对粒子i的切应力产生的扭矩为：

T_ij＝-r_ij×F_ij 公式(17)

式中

T_ij——粒子j对粒子i的切应力产生的扭矩，N·m；

r_ij——粒子i和粒子j的相对位置矢量，m；

F_ij——粒子j对粒子i的切应力，N；

粒子i的角速度由公式(18)得：

式中

ω_i——粒子i的角速度，rad/s；

t——时间，s；

I_i——粒子i的转动惯量，kg·m²；

T_ij——粒子j对粒子i的切应力产生的扭矩，N·m；

通过步骤8的计算，计算得到每个粒子的角速度，由角速度计算得到旋转角度，从而计算步骤6中的旋转矩阵R；该步骤保证了粒子的运动满足角动量定律，即得到了燃料元件变形过程中的固体旋转或扭转随时间的变化过程；

步骤9：总应变为弹性应变加塑性应变加热膨胀应变，求解塑性应变基于应变增量理论；假设在计算塑性应变之前，先假设无塑性应变，由该假设得弹性应变计算等效应力，并将等效应力与屈服极限比较，如果等效应力大于屈服极限，则认为发生塑性变形，等效应力求解如公式(19)至公式(21)所示，

[ds]_n＝2μ[dε^e]_n 公式(19)

[s]_n＝[s]_n-1+[ds]_n 公式(20)

式中

[ds]_n——第n时间步下的应力增量张量，N/m²；

μ——拉梅常数的第二个参数，

[dε^e]_n——第n时间步下的弹性应变增量张量；

[s]_n——第n时间步下的应变张量，N/m²；

[s]_n-1——第n-1时间步下的应变张量，N/m²；

——粒子i的等效应力，N/m²；

s_ij——应力矢量，N/m²；

通过步骤9的计算，得到了每个粒子在不同时刻下是否发生塑性变形的判定，即判定燃料元件的不同位置是否发生塑性变形；

步骤10：塑性应力应变计算，塑性应变增量如公式(22)和(23)所示，

[dε]_n＝[dε^p]_n+[dε^e]_n+[dε^T]_n 公式(23)

式中

[ds^p]_n——第n时间步下的塑性应力增量张量，N/m²；

[ds]_n——第n时间步下的应力增量张量，N/m²；

[ds^T]_n——第n时间步下的热膨胀应力增量张量，N/m²；

[ds^e]_n——第n时间步下的弹性应力增量张量，N/m²；

[s]_n-1——第n-1时间步下的应变张量，N/m²；

μ——拉梅常数的第二个参数，

dλ_n——增量参数，由材料的力学特性决定；

式中dλ_n的取值影响塑性应变和弹性应变的关系，不同特性的材料，其dλ_n的求解不同，以下针对线性弹塑性，dλ_n的求解如公式(24)所示，

式中

——假设所有应变均为弹性应变条件下的等效应力，N/m²；

H′——塑性硬化系数；

Y——屈服极限，N/m²；

——第n-1时间步下的等效塑性应变；

μ——拉梅常数的第二个参数，

由塑性模型获得塑性应变后，由公式(23)计算得到弹性应变，再采用步骤7中的弹性模型计算弹性应力；

通过步骤10的计算，模拟了核反应堆燃料早期变形过程中的塑性变形过程；计算得到了每个粒子在不同时刻下的塑性应变和弹性应变，即得到了燃料元件的塑性变形和弹性变形随时间的变化过程；

步骤11：断裂判断，当粒子间距大于或小于断裂阈值时，则认为发生断裂，断裂后断裂部分的2号内部粒子转化为1号边界粒子，且断裂粒子之间不再发生固体连续介质内的相互作用，仅作为固体边界参与边界载荷和碰撞的计算；根据外力载荷、碰撞和弹性力更新粒子速度、位置、角速度和旋转角度；

通过步骤11的计算，模拟了核反应堆燃料早期变形过程中的断裂过程；计算得到了每个粒子在不同时刻下的是否发生断裂的判定，即获得了燃料元件发生断裂的位置；

通过步骤1对实际的核反应堆燃料结构建立粒子几何模型，设定粒子的位置、速度、温度、物性参数和力学参数；通过步骤2存储粒子邻居域，为后续计算提供基础；通过步骤3至步骤10计算固体的温度和受力，受力包括外力、碰撞、弹性应力应变、塑性应变和热膨胀应变，更新粒子的速度和位置；通过步骤11判断粒子是否发生断裂；由步骤3至步骤10计算得到的加速度，更新所有粒子的速度和位置；对计算是否结束进行判定，若是，输出结果，若否，更新固体间粒子邻居域，再返回步骤3；综合以上步骤，模拟了核反应堆燃料早期变形过程，得到了变形过程中燃料元件的位置、速度、温度、应力状态、应变状态，通过以上数据对核反应堆燃料早期变形过程中的传热、热膨胀、弹性变形、塑性变形和受力情况展开机理性分析。