CN113191065B - 基于粒子法的核反应堆燃料早期变形分析方法 - Google Patents

Abstract

一种基于粒子法的核反应堆燃料早期变形分析方法，主要步骤如下：1、建立几何模型，设置粒子物性和力学参数；2、建立粒子邻居域；3、基于边界条件和固体碰撞模型计算固体边界粒子外力载荷；4、计算温度并基于热膨胀模型计算热膨胀应变；5、基于弹性模型、弹塑性模型和塑性模型计算固体内连续介质粒子的塑性应变和弹性应力应变；6、判断是否发生断裂并更新粒子的速度和位置；7、输出计算结果。本方法能够分析核反应堆在严重事故中燃料早期变形行为；基于粒子法，能够有效避免了网格法中存在的网格畸变问题；本方法为纯拉格朗日核方法，能够有效避免应力应变分析过程中存在的拉伸不稳定性。

Description

基于粒子法的核反应堆燃料早期变形分析方法

技术领域

本发明涉及核反应堆严重事故下燃料早期变形研究技术领域，具体涉及一种基于粒子法的核反应堆燃料早期变形分析方法。

背景技术

当核电厂发生事故时，核反应堆堆芯如果得不到充分的冷却，燃料的状态会偏离正常运行工况，其温度上升导致堆芯结构材料和燃料发生热膨胀，从而使得燃料整体发生变形并承受额外的应力。燃料的变形使得冷却剂通道形状和流通面积发生改变，影响后续应急冷却措施的效果。并且应急水的注入等措施会对堆芯结构及燃料产生外力载荷，这种情况下，处于高温的燃料极易发生变形，存在很大不确定性。如果堆芯温度进一步上升，燃料开始发生失效，堆芯结构材料、燃料包壳甚至芯块可能发生熔化并向下迁移，而燃料的早期变形为后续严重事故中燃料熔化及熔融物的迁移提供了初始条件和边界条件，影响严重事故进程。特别地，在新型燃料元件的情况中，如板状燃料元件、弥散体微球等，其燃料元件变形行为对严重事故安全分析格外重要。在核反应堆严重事故过程中，燃料受到的环境情况复杂，其变形行为存在很大不确定，变形行为对严重事故的影响机理仍未被完全研究透彻，对其的研究有助于严重事故应急措施的分析，对于核电厂严重事故安全分析具有重大意义。

关于核反应堆严重事故下燃料早期变形的研究，国内外已经展开了一些研究，包括实验研究和数值模拟研究。实验研究主要为分离效应性研究，研究燃料元件在不同外力载荷、温度、辐照条件下的应力应变情况，对于整体性实验，这类的侧重点在于燃料熔化后的行为，早期变形行为的研究较少。对于数值模拟研究，大多采用基于网格法的CFD软件对具体的相互作用过程进行分析，但基于网格法的CFD软件在分析大变形过程时，极易存在网格畸变的问题，且重塑网格的计算耗费过大。而对于基于拉格朗日方法的粒子法，其在处理自由表面运动上具有独特的优势，且不存在网格畸变的问题。在海洋研究领域中，粒子法成功应用于弹性结构材料与海浪相互作用的研究当中，但在核反应堆内严重事故下复杂情况应用的研究较少。因此，本研究综合粒子法及核反应堆严重事故下燃料早期变形行为的机理性分析提出一种核反应堆燃料早期变形分析方法。

发明内容

为了全面研究核反应堆燃料早期变形过程，揭示变形过程中可能存在的一些机理现象，本发明在对核反应堆燃料早期变形的机理性分析的基础上结合粒子法、基础控制方程及相关的机力学结构模型，提出一种研究核反应堆燃料早期变形的分析方法，该方法能够对核反应堆燃料早期变形中存在的弹性变形、塑性变形、热膨胀变形、传热过程进行研究，获得核反应堆燃料早期变形中燃料变形规律及程度，为核电厂反应堆严重事故安全特性研究提供重要依据。

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

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

步骤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为固体内粒子邻居域内所有粒子对中心粒子的核函数值之和、粒子间粒子数密度n_outer为固体间粒子邻居域内所有粒子对中心粒子的核函数值之和、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；综合以上步骤，模拟了核反应堆燃料早期变形过程，得到了变形过程中燃料元件的位置、速度、温度、应力状态、应变状态，通过以上数据对核反应堆燃料早期变形过程中的传热、热膨胀、弹性变形、塑性变形和受力情况展开机理性分析。

本发明方法为核反应堆严重事故中燃料早期变形行为分析提供解决方案，为核电厂反应堆严重事故安全特性的研究提供重要依据。

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

本发明的基于粒子法的核反应堆燃料早期变形分析方法，综合考虑了核反应堆严重事故中燃料早期变形行为可能出现的大部分现象，包括弹性变形、塑性变形、热膨胀和传热过程，能够对核反应堆燃料早期变形行为进行机理性分析。该方法适用范围广，不仅能够针对传统的棒状燃料元件进行分析，还能对板状燃料元件、弥散体微球等新型燃料进行早期变形分析。并且该方法基于粒子法，能够有效避免了网格法中存在的网格畸变问题。本方法为纯拉格朗日核方法，能够有效避免应力应变分析过程中存在的拉伸不稳定性。此外，整个分析过程均可采用显式求解，计算资源耗费小，并行效率高。综上，该方法能够更加全面、有效、高效地对核反应堆燃料早期变形行为进行分析。

附图说明

图1是本发明基于粒子法的核反应堆燃料早期变形分析方法的流程图。

具体实施方式

本发明基于粒子法的核反应堆燃料早期变形分析方法，如图1所示。以下以单根燃料元件早期变形为例，其具体实施步骤如下：

步骤1：根据单根燃料元件的结构建立粒子几何模型，采用0.1毫米粒子大小，建立的燃料棒几何模型长3.658米，芯块(采用M＝0的粒子)直径8.2毫米，高度13.5毫米，为蝶形加倒角的形式，高度方向布置270个芯块，包壳(采用M＝1的粒子)内径8.4毫米，外径9.5毫米，燃料棒上下端采用结构材料(采用M＝2的粒子)简化为固定结构，其速度和位置在整个计算过程中保持不变；270个芯块的Object数分别0至269；包壳视为整体，其Object数设为270；简化的结构材料分为上下两部分，Object数设为271和272；该算例中有273个固体，每个固体的最外一层粒子设置为1号粒子，其余粒子均为2号粒子；设置燃料元件整体速度为零，无预应力，无应变，整体温度相同，为1000K，芯块设置体积释热率2.4×10⁸W/m³，其余材料无热源；芯块、包壳和结构材料分别采用UO₂、Zr-4合金和不锈钢的真实物性；

步骤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为固体内粒子邻居域内所有粒子对中心粒子的核函数值之和、粒子间粒子数密度n_outer为固体间粒子邻居域内所有粒子对中心粒子的核函数值之和、n⁰为参考粒子数密度，是假设粒子均匀分布情况下计算得到的粒子数密度；

步骤3：基于边界条件和固体碰撞模型计算固体边界粒子外力载荷，在本算例中，没有外力载荷，只设置上下端的结构材料固定，因此外力载荷引入的加速度始终为零；

式中

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，计算得到了所有1号粒子由于碰撞引起的加速度，即每个固体之间碰撞引发的变形和运动情况；

步骤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；

以1000K为参考温度，计算燃料元件由于温度升高而引起的热膨胀应力，得到每个粒子在不同时刻下的热膨胀应变张量，即得到了燃料元件由于温度变化而发生碰撞变形的行为规律；

步骤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：断裂判断，本算例中，当粒子间距大于1.2倍粒子直径或小于0.8倍粒子直径时，认为发生断裂；断裂后将断裂部分的2号内部粒子转化为1号边界粒子，且断裂粒子之间不再发生固体连续介质内的相互作用，仅作为固体边界参与边界载荷和碰撞的计算；根据外力载荷、碰撞和弹性力更新粒子速度、位置、角速度和旋转角度；

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

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

综上综合以上步骤，模拟了单根燃料元件的早期变形过程，得到了变形过程中燃料元件的位置、速度、温度、应力状态、应变状态，通过以上数据对单根燃料元件早期变形过程中的传热、热膨胀、弹性变形、塑性变形和受力情况展开机理性分析。

参考单根燃料元件早期变形的分析过程，类似地，可以实现多根燃料元件或一个燃料组件的早期变形行为，只要额外考虑燃料元件与燃料元件之间的受力；在计算资源充足的情况下，甚至可以建立整个堆芯的粒子几何模型，综合考虑整个堆芯的早期变形行为；此外，通过添加外力载荷，可以研究更加恶劣的应力情况下，燃料元件的变形行为；在实际反应堆燃料元件的设计中，往往有预应力或预冲压的措施，对于这类情况，本发明可以通过设置表面粒子压力的方式实现，类似于外力载荷的方式。

Claims (1)

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；综合以上步骤，模拟了核反应堆燃料早期变形过程，得到了变形过程中燃料元件的位置、速度、温度、应力状态、应变状态，通过以上数据对核反应堆燃料早期变形过程中的传热、热膨胀、弹性变形、塑性变形和受力情况展开机理性分析。

Images