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

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

Info

Publication number
CN113191065A
CN113191065A CN202110486476.0A CN202110486476A CN113191065A CN 113191065 A CN113191065 A CN 113191065A CN 202110486476 A CN202110486476 A CN 202110486476A CN 113191065 A CN113191065 A CN 113191065A
Authority
CN
China
Prior art keywords
particle
particles
strain
elastic
stress
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202110486476.0A
Other languages
English (en)
Other versions
CN113191065B (zh
Inventor
陈荣华
蔡庆航
肖鑫坤
王金顺
李勇霖
田文喜
苏光辉
秋穗正
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Xian Jiaotong University
Original Assignee
Xian Jiaotong University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Xian Jiaotong University filed Critical Xian Jiaotong University
Priority to CN202110486476.0A priority Critical patent/CN113191065B/zh
Publication of CN113191065A publication Critical patent/CN113191065A/zh
Application granted granted Critical
Publication of CN113191065B publication Critical patent/CN113191065B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00Computer-aided design [CAD]
    • G06F30/20Design optimisation, verification or simulation
    • G06F30/25Design optimisation, verification or simulation using particle-based methods
    • GPHYSICS
    • G16INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16CCOMPUTATIONAL CHEMISTRY; CHEMOINFORMATICS; COMPUTATIONAL MATERIALS SCIENCE
    • G16C10/00Computational theoretical chemistry, i.e. ICT specially adapted for theoretical aspects of quantum chemistry, molecular mechanics, molecular dynamics or the like
    • GPHYSICS
    • G16INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16CCOMPUTATIONAL CHEMISTRY; CHEMOINFORMATICS; COMPUTATIONAL MATERIALS SCIENCE
    • G16C20/00Chemoinformatics, i.e. ICT specially adapted for the handling of physicochemical or structural data of chemical particles, elements, compounds or mixtures
    • G16C20/10Analysis or design of chemical reactions, syntheses or processes
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/08Thermal analysis or thermal optimisation
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/14Force analysis or force optimisation, e.g. static or dynamic forces

Landscapes

  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • Computing Systems (AREA)
  • Physics & Mathematics (AREA)
  • Chemical & Material Sciences (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • General Engineering & Computer Science (AREA)
  • Health & Medical Sciences (AREA)
  • General Health & Medical Sciences (AREA)
  • Spectroscopy & Molecular Physics (AREA)
  • General Physics & Mathematics (AREA)
  • Geometry (AREA)
  • Evolutionary Computation (AREA)
  • Computer Hardware Design (AREA)
  • Analytical Chemistry (AREA)
  • Chemical Kinetics & Catalysis (AREA)
  • Crystallography & Structural Chemistry (AREA)
  • Monitoring And Testing Of Nuclear Reactors (AREA)

Abstract

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

Description

基于粒子法的核反应堆燃料早期变形分析方法
技术领域
本发明涉及核反应堆严重事故下燃料早期变形研究技术领域,具体涉及一种基于粒子法的核反应堆燃料早期变形分析方法。
背景技术
当核电厂发生事故时,核反应堆堆芯如果得不到充分的冷却,燃料的状态会偏离正常运行工况,其温度上升导致堆芯结构材料和燃料发生热膨胀,从而使得燃料整体发生变形并承受额外的应力。燃料的变形使得冷却剂通道形状和流通面积发生改变,影响后续应急冷却措施的效果。并且应急水的注入等措施会对堆芯结构及燃料产生外力载荷,这种情况下,处于高温的燃料极易发生变形,存在很大不确定性。如果堆芯温度进一步上升,燃料开始发生失效,堆芯结构材料、燃料包壳甚至芯块可能发生熔化并向下迁移,而燃料的早期变形为后续严重事故中燃料熔化及熔融物的迁移提供了初始条件和边界条件,影响严重事故进程。特别地,在新型燃料元件的情况中,如板状燃料元件、弥散体微球等,其燃料元件变形行为对严重事故安全分析格外重要。在核反应堆严重事故过程中,燃料受到的环境情况复杂,其变形行为存在很大不确定,变形行为对严重事故的影响机理仍未被完全研究透彻,对其的研究有助于严重事故应急措施的分析,对于核电厂严重事故安全分析具有重大意义。
关于核反应堆严重事故下燃料早期变形的研究,国内外已经展开了一些研究,包括实验研究和数值模拟研究。实验研究主要为分离效应性研究,研究燃料元件在不同外力载荷、温度、辐照条件下的应力应变情况,对于整体性实验,这类的侧重点在于燃料熔化后的行为,早期变形行为的研究较少。对于数值模拟研究,大多采用基于网格法的CFD软件对具体的相互作用过程进行分析,但基于网格法的CFD软件在分析大变形过程时,极易存在网格畸变的问题,且重塑网格的计算耗费过大。而对于基于拉格朗日方法的粒子法,其在处理自由表面运动上具有独特的优势,且不存在网格畸变的问题。在海洋研究领域中,粒子法成功应用于弹性结构材料与海浪相互作用的研究当中,但在核反应堆内严重事故下复杂情况应用的研究较少。因此,本研究综合粒子法及核反应堆严重事故下燃料早期变形行为的机理性分析提出一种核反应堆燃料早期变形分析方法。
发明内容
为了全面研究核反应堆燃料早期变形过程,揭示变形过程中可能存在的一些机理现象,本发明在对核反应堆燃料早期变形的机理性分析的基础上结合粒子法、基础控制方程及相关的机力学结构模型,提出一种研究核反应堆燃料早期变形的分析方法,该方法能够对核反应堆燃料早期变形中存在的弹性变形、塑性变形、热膨胀变形、传热过程进行研究,获得核反应堆燃料早期变形中燃料变形规律及程度,为核电厂反应堆严重事故安全特性研究提供重要依据。
为了实现上述目标,本发明采取了以下的技术方案予以实施:
一种基于粒子法的核反应堆燃料早期变形分析方法,步骤如下:
步骤1:根据实际的核反应堆燃料结构建立粒子几何模型,定义材料数M区分不同的材料,M=0,1,2分别表示核反应堆燃料元件的芯块、包壳和结构材料,粒子根据对应的材料携带对应的物性和力学参数;定义ID数区分固体边界,采用1号粒子模拟固体边界,2号粒子模拟固体内连续介质区域;不同固体之间采用参数Object区分,对于2号粒子,只有当Object数相同时才发生相互作用;初始化所有粒子的速度、位置、温度、应力张量和应变张量,分别为ui、ri、Ti、σi、εi;通过以上参数设定,可以有效建立燃料芯块、包壳与结构材料三者之间的固体边界,还原真实的燃料结构,并确定了燃料的初始位置、速度、温度、应力应变状态和相关物性参数;
步骤2:根据建立的初始粒子几何模型,为每个粒子建立固体内粒子邻居域,即存储每个粒子以自身为圆心,作用范围re为半径的区域内所有Object数与该粒子相同的粒子ID,固体内邻居域不随计算时间的推移而发生改变,即认为固体内连续介质发生相互作用的粒子始终保持不变;除了固体内粒子邻居域,还需要为所有1号粒子建立固体间粒子邻居域,即存储每个1号粒子以自身为圆心,作用范围re为半径的区域内所有Object数与该粒子不同的1号粒子ID,固体间邻居域会随着计算时间的推移而发生改变,即两个固体的相对位置会不断发生变化,当两者表面距离小于作用范围re时,两个固体的边界会发生相互作用,否则不会;定义粒子数密度为n,其值为粒子邻居域内所有粒子对中心粒子的核函数值之和,如公式(1)所示;
Figure BDA0003050542990000021
式中
ni——粒子i的粒子数密度;
w(rij,re)——核函数,采用的核函数形式为
Figure BDA0003050542990000022
rij——粒子i和粒子j的距离,m;
re——粒子作用范围,m;
下标i——任意粒子i;
下标j——粒子i的粒子邻居域内的任意粒子j;
由此,可以进一步定义固体内粒子数密度ninner为固体内粒子邻居域内所有粒子对中心粒子的核函数值之和、粒子间粒子数密度nouter为固体间粒子邻居域内所有粒子对中心粒子的核函数值之和、n0为参考粒子数密度,是假设粒子均匀分布情况下计算得到的粒子数密度;
步骤3:基于边界条件和固体碰撞模型计算固体边界粒子外力载荷,根据具体实际物理问题,为1号粒子设置外力的作用方向和大小,由公式(2)计算边界条件外力载荷的影响,边界条件外力载荷始终视为作用于粒子中心,因此不考虑该力所引起的粒子转动,公式采用张量表示法和下标表示法;
Figure BDA0003050542990000023
式中
vα——α方向上的速度分量,m/s;
Figure BDA0003050542990000031
——α方向上的边界条件外力载荷引起的加速度,m/s2
t——时间,s;
对于固体碰撞过程,采用固体碰撞模型计算1号粒子的碰撞载荷,如公式(3)所示,
Figure BDA0003050542990000032
式中
Pi Contact——粒子i的接触应力,N/m2
λ——拉梅常数的第一个参数,
Figure BDA0003050542990000033
γγ)i——粒子i的应变张量对角项之和,即εγγ=ε112233
nouter——粒子i的固体间邻居域的粒子数密度;
ninner——粒子i的固体内邻居域的粒子数密度;
E——杨氏模量,N/m2
υ——泊松比;
由碰撞载荷,计算对应的速度梯度,如公式(4)所示,
Figure BDA0003050542990000034
式中
vi——粒子i的速度矢量,m/s;
t——时间,s;
ρi——粒子i的密度,kg/m3
d——空间维度;
n0——参考粒子数密度;
rij——粒子i和粒子j的相对位置矢量,m;
Pi Contact——粒子i的接触应力,N/m2
Figure BDA0003050542990000035
——粒子j的接触应力,N/m2
Figure BDA0003050542990000036
——粒子i和粒子j的相对初始位置矢量,m;
Figure BDA0003050542990000041
——核函数,等同于
Figure BDA0003050542990000042
通过步骤3的计算,模拟核反应堆燃料早期变形过程中燃料元件受到的外力载荷和碰撞力,全面考虑了核反应堆燃料早期变形过程中可能发生的外力作用、燃料元件之间的碰撞、燃料元件与结构材料之间的碰撞、燃料包壳与芯块之间的碰撞、芯块与芯块之间的碰撞过程;计算得到每个粒子由外力载荷和碰撞引起的加速度,即燃料元件由外力载荷和碰撞引起的加速度;
步骤4:能量计算如公式(5)所示:
Figure BDA0003050542990000043
式中
h——焓值,J/kg;
t——时间,s;
ρ——密度,kg/m3
k——热导率,W/(m·K);
T——温度,K;
Qradiation——辐射热源,W/m3
Qout——外热源,W/m3
辐射换热采用斯忒藩-玻尔兹曼定律,如公式(6)所示
Figure BDA0003050542990000044
式中
Qradiation——辐射热源W/m3
ε——发射率;
σ——斯忒藩-玻耳兹曼常数;
T——粒子温度K;
Tenv——环境温度K;
l0——粒子直径m;
温度的拉普拉斯项采用拉普拉斯模型离散得到:
Figure BDA0003050542990000045
式中
d——维度;
n0——参考粒子数密度;
ρi——粒子i的密度,kg/m3
Figure BDA0003050542990000046
——当前时刻的粒子j温度K;
Ti k——当前时刻的粒子i温度K;
rij——粒子i和粒子j的相对位置矢量,m;
w(|rij|)——核函数,等同于w(|rij|,re);
Figure BDA0003050542990000051
由焓值计算粒子温度:
Figure BDA0003050542990000052
式中
T——温度K;
Ts——熔点K;
h——焓值J/kg;
hs0——开始熔化的焓值J/kg;
hs1——结束熔化的焓值J/kg;
cp——比热容J/(kg·K);
通过步骤4的计算,模拟核反应堆燃料早期变形过程中燃料元件内的导热、燃料元件间的辐射换热、芯块和包壳间的换热过程;计算得到每个粒子在不同时刻下的焓值和温度,即得到核反应堆燃料的焓值和温度随时间的变化过程;
步骤5:热膨胀应变计算如公式(8)所示,
Figure BDA0003050542990000053
式中
[dsT]i——粒子i的热膨胀应力增量张量,N/m2
κi——粒子i的热膨胀系数;
Ti——粒子i的温度,K;
Ti ref——粒子i的参考温度,K;
通过步骤5的计算,模拟核反应堆燃料早期变形过程中燃料元件的热膨胀变形过程;计算得到每个粒子在不同时刻下的热膨胀应变张量,即得到了燃料元件由于温度变化而发生碰撞变形的行为规律;
步骤6:计算总应变,固体内的相对变形是由当前相对位置减去旋转后的相对位置得到的,如公式(9)所示,
Figure BDA0003050542990000061
式中
uij——粒子i和粒子j的相对总变形矢量,m;
rij——粒子i和粒子j的相对位置矢量,m;
R——粒子i的旋转矩阵,R=RxRyRz
Figure BDA0003050542990000062
——粒子i和粒子j的相对初始位置矢量,m;
Rx——粒子i绕x轴旋转的旋转矩阵,
Figure BDA0003050542990000063
Ry——粒子i绕y轴旋转的旋转矩阵,
Figure BDA0003050542990000064
Rz——粒子i绕z轴旋转的旋转矩阵,
Figure BDA0003050542990000065
θα,ij——粒子i和粒子j绕α轴旋转的相对角度,°,
Figure BDA0003050542990000066
α可以取x,y,z;
由总变形可以计算得到总应变:
Figure BDA0003050542990000067
εij——粒子i和粒子j的相对总应变矢量;
通过步骤6的计算,计算得到了每个粒子在不同时刻下的总应变张量,即得到了核反应堆燃料元件的总变形状态;
步骤7:基于弹性模型计算弹性应力,弹性模型如公式(11)所示,
Figure BDA0003050542990000068
式中
Figure BDA0003050542990000069
——弹性应力张量的分量,N/m2
λ——拉梅常数的第一个参数,
Figure BDA0003050542990000071
μ——拉梅常数的第二个参数,
Figure BDA0003050542990000072
Figure BDA0003050542990000073
——弹性应变张量对角项之和,即
Figure BDA0003050542990000074
δαβ——克罗内克尔函数,
Figure BDA0003050542990000075
Figure BDA0003050542990000076
——弹性应变张量的分量;
其中记各向同性力:
Figure BDA0003050542990000077
式中
Pi E——粒子i的弹性各向同性力,N/m2
λ——拉梅常数的第一个参数,
Figure BDA0003050542990000078
Figure BDA0003050542990000079
——弹性应变张量对角项之和,即
Figure BDA00030505429900000710
d——空间维度;
n0——参考粒子数密度;
uij——粒子i和粒子j的相对总变形矢量,m;
rij——粒子i和粒子j的相对位置矢量,m;
Figure BDA00030505429900000711
——粒子i和粒子j的相对初始位置矢量,m;
Figure BDA00030505429900000712
——核函数,等同于
Figure BDA00030505429900000713
此外,在碰撞计算中,如果碰撞载荷Pi Contact小于弹性各向同性力Pi E,则令Pi Contact=Pi E;记各向异性力:
Figure BDA00030505429900000714
式中
Figure BDA00030505429900000715
——粒子i和粒子j间的弹性各向异性力,N/m2
μ——拉梅常数的第二个参数,
Figure BDA0003050542990000081
Figure BDA0003050542990000082
——粒子i和粒子j的相对弹性应变矢量;
由弹性各向同性力和弹性各向异性力分别求得对应的速度梯度,如公式(14)和(15)所示,
Figure BDA0003050542990000083
Figure BDA0003050542990000084
式中
vi——粒子i的速度矢量,m/s;
t——时间,s;
ρi——粒子i的密度,kg/m3
d——空间维度;
n0——参考粒子数密度;
rij——粒子i和粒子j的相对位置矢量,m;
Pi E——粒子i的弹性各向同性力,N/m2
Figure BDA0003050542990000085
——粒子j的弹性各向同性力,N/m2
Figure BDA0003050542990000086
——粒子i和粒子j间的弹性各向异性力,N/m2
Figure BDA0003050542990000087
——粒子i和粒子j的相对初始位置矢量,m;
Figure BDA0003050542990000088
——核函数,等同于
Figure BDA0003050542990000089
通过步骤7,模拟了核反应堆燃料早期变形过程中的弹性变形;计算得到每个粒子在不同时刻下的弹性应变张量、弹性应力张量、弹性力所引起的加速度,即得到燃料元件的弹性应力应变随时间的变化过程;
步骤8:基于角动量定理计算粒子旋转角速度,粒子j对粒子i的切应力为:
Figure BDA00030505429900000810
式中
Fij——粒子j对粒子i的切应力,N;
d——空间维度;
n0——参考粒子数密度;
mi——粒子i的质量,kg;
ρi——粒子i的密度,kg/m3
Figure BDA0003050542990000091
——粒子i和粒子j间的弹性各向异性力在与rij垂直方向上的投影,N/m2
rij——粒子i和粒子j的相对位置矢量,m;
Figure BDA0003050542990000092
——粒子i和粒子j的相对初始位置矢量,m;
Figure BDA0003050542990000093
——核函数,等同于
Figure BDA0003050542990000094
粒子j对粒子i的切应力产生的扭矩为:
Tij=-rij×Fij 公式(17)
式中
Tij——粒子j对粒子i的切应力产生的扭矩,N·m;
rij——粒子i和粒子j的相对位置矢量,m;
Fij——粒子j对粒子i的切应力,N;
粒子i的角速度可由公式(18)得:
Figure BDA0003050542990000095
式中
ωi——粒子i的角速度,rad/s;
t——时间,s;
Ii——粒子i的转动惯量,kg·m2
Tij——粒子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)
Figure BDA0003050542990000101
式中
[ds]n——第n时间步下的应力增量张量,N/m2
μ——拉梅常数的第二个参数,
Figure BDA0003050542990000102
[dεe]n——第n时间步下的弹性应变增量张量;
[s]n——第n时间步下的应变张量,N/m2
[s]n-1——第n-1时间步下的应变张量,N/m2
Figure BDA0003050542990000103
——粒子i的等效应力,N/m2
sij——应力矢量,N/m2
通过步骤9的计算,得到了每个粒子在不同时刻下是否发生塑性变形的判定,即判定燃料元件的不同位置是否发生塑性变形;
步骤10:塑性应力应变计算,塑性应变增量如公式(22)和(23)所示,
Figure BDA0003050542990000104
[dε]n=[dεp]n+[dεe]n+[dεT]n 公式(23)
式中
[dsp]n——第n时间步下的塑性应力增量张量,N/m2
[ds]n——第n时间步下的应力增量张量,N/m2
[dsT]n——第n时间步下的热膨胀应力增量张量,N/m2
[dse]n——第n时间步下的弹性应力增量张量,N/m2
[s]n-1——第n-1时间步下的应变张量,N/m2
μ——拉梅常数的第二个参数,
Figure BDA0003050542990000105
n——增量参数,由材料的力学特性决定;
式中dλn的取值影响塑性应变和弹性应变的关系,不同特性的材料,其dλn的求解不同,以下针对线性弹塑性,dλn的求解如公式(24)所示,
Figure BDA0003050542990000111
式中
Figure BDA0003050542990000112
——假设所有应变均为弹性应变条件下的等效应力,N/m2
H′——塑性硬化系数;
Y——屈服极限,N/m2
Figure BDA0003050542990000113
——第n-1时间步下的等效塑性应变;
μ——拉梅常数的第二个参数,
Figure BDA0003050542990000114
由塑性模型获得塑性应变后,由公式(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×108W/m3,其余材料无热源;芯块、包壳和结构材料分别采用UO2、Zr-4合金和不锈钢的真实物性;
步骤2:根据建立的初始粒子几何模型,为每个粒子建立固体内粒子邻居域,即存储每个粒子以自身为圆心,作用范围re为半径的区域内所有Object数与该粒子相同的粒子ID,固体内邻居域不随计算时间的推移而发生改变,即认为固体内连续介质发生相互作用的粒子始终保持不变;除了固体内粒子邻居域,还需要为所有1号粒子建立固体间粒子邻居域,即存储每个1号粒子以自身为圆心,作用范围re为半径的区域内所有Object数与该粒子不同的1号粒子ID,固体间邻居域会随着计算时间的推移而发生改变,即两个固体的相对位置会不断发生变化,当两者表面距离小于作用范围re时,两个固体的边界会发生相互作用,否则不会;定义粒子数密度为n,其值为粒子邻居域内所有粒子对中心粒子的核函数值之和,如公式(1)所示;
Figure BDA0003050542990000121
式中
ni——粒子i的粒子数密度;
w(rij,re)——核函数,本发明中采用的核函数形式为
Figure BDA0003050542990000122
rij——粒子i和粒子j的距离,m;
re——粒子作用范围,m;
下标i——任意粒子i;
下标j——粒子i的粒子邻居域内的任意粒子j;
由此,可以进一步定义固体内粒子数密度ninner为固体内粒子邻居域内所有粒子对中心粒子的核函数值之和、粒子间粒子数密度nouter为固体间粒子邻居域内所有粒子对中心粒子的核函数值之和、n0为参考粒子数密度,是假设粒子均匀分布情况下计算得到的粒子数密度;
步骤3:基于边界条件和固体碰撞模型计算固体边界粒子外力载荷,在本算例中,没有外力载荷,只设置上下端的结构材料固定,因此外力载荷引入的加速度始终为零;
Figure BDA0003050542990000131
式中
vα——α方向上的速度分量,m/s;
Figure BDA0003050542990000132
——α方向上的边界条件外力载荷引起的加速度,m/s2
t——时间,s;
在本算例中,需要计算芯块与芯块之间、芯块与包壳之间、芯块与结构材料之间、包壳与结构材料之间的碰撞过程,采用固体碰撞模型计算1号粒子的碰撞载荷,如公式(3)所示,
Figure BDA0003050542990000133
式中
Pi Contact——粒子i的接触应力,N/m2
λ——拉梅常数的第一个参数,
Figure BDA0003050542990000134
γγ)i——粒子i的应变张量对角项之和,即εγγ=ε112233
nouter——粒子i的固体间邻居域的粒子数密度;
ninner——粒子i的固体内邻居域的粒子数密度;
E——杨氏模量,N/m2
υ——泊松比;
由碰撞载荷,计算对应的速度梯度,如公式(4)所示,
Figure BDA0003050542990000135
式中
vi——粒子i的速度矢量,m/s;
t——时间,s;
ρi——粒子i的密度,kg/m3
d——空间维度;
n0——参考粒子数密度;
rij——粒子i和粒子j的相对位置矢量,m;
Pi Contact——粒子i的接触应力,N/m2
Figure BDA0003050542990000141
——粒子j的接触应力,N/m2
Figure BDA0003050542990000142
——粒子i和粒子j的相对初始位置矢量,m;
Figure BDA0003050542990000143
——核函数,等同于
Figure BDA0003050542990000144
通过步骤3,计算得到了所有1号粒子由于碰撞引起的加速度,即每个固体之间碰撞引发的变形和运动情况;
步骤4:本算例中,对于结构材料粒子和包壳粒子需要计算导热项和辐射换热项,而芯块粒子需要计算导热项、辐射换热项和外热源项,能量计算如公式(5)所示:
Figure BDA0003050542990000145
式中
h——焓值,J/kg;
t——时间,s;
ρ——密度,kg/m3
k——热导率,W/(m·K);
T——温度,K;
Qradiation——辐射热源,W/m3
Qout——外热源,W/m3
辐射换热采用斯忒藩-玻尔兹曼定律,如公式(6)所示
Figure BDA0003050542990000146
式中
Qradiation——辐射热源W/m3
ε——发射率;
σ——斯忒藩-玻耳兹曼常数;
T——粒子温度K;
Tenv——环境温度K;
l0——粒子直径m;
温度的拉普拉斯项采用拉普拉斯模型离散得到:
Figure BDA0003050542990000151
式中
d——维度;
n0——参考粒子数密度;
ρi——粒子i的密度,kg/m3
Figure BDA0003050542990000152
——当前时刻的粒子j温度K;
Ti k——当前时刻的粒子i温度K;
rij——粒子i和粒子j的相对位置矢量,m;
w(|rij|)——核函数,等同于w(|rij|,re);
Figure BDA0003050542990000153
由焓值计算粒子温度:
Figure BDA0003050542990000154
式中
T——温度K;
Ts——熔点K;
h——焓值J/kg;
hs0——开始熔化的焓值J/kg;
hs1——结束熔化的焓值J/kg;
cp——比热容J/(kg·K);
通过步骤4的计算,模拟单根燃料元件内包壳、芯块和结构材料之间的导热和辐射换热,以及芯块的释热过程;计算得到每个粒子在不同时刻下的焓值和温度,即得到单个燃料元件的焓值和温度随时间的变化过程;
步骤5:热膨胀应变计算如公式(8)所示,
Figure BDA0003050542990000161
式中
[dsT]i——粒子i的热膨胀应力增量张量,N/m2
κi——粒子i的热膨胀系数;
Ti——粒子i的温度,K;
Ti ref——粒子i的参考温度,K;
以1000K为参考温度,计算燃料元件由于温度升高而引起的热膨胀应力,得到每个粒子在不同时刻下的热膨胀应变张量,即得到了燃料元件由于温度变化而发生碰撞变形的行为规律;
步骤6:计算总应变,固体内的相对变形是由当前相对位置减去旋转后的相对位置得到的,如公式(9)所示,
Figure BDA0003050542990000162
式中
uij——粒子i和粒子j的相对总变形矢量,m;
rij——粒子i和粒子j的相对位置矢量,m;
R——粒子i的旋转矩阵,R=RxRyRz
Figure BDA0003050542990000163
——粒子i和粒子j的相对初始位置矢量,m;
Rx——粒子i绕x轴旋转的旋转矩阵,
Figure BDA0003050542990000164
Ry——粒子i绕y轴旋转的旋转矩阵,
Figure BDA0003050542990000165
Rz——粒子i绕z轴旋转的旋转矩阵,
Figure BDA0003050542990000166
θα,ij——粒子i和粒子j绕α轴旋转的相对角度,°,
Figure BDA0003050542990000171
α可以取x,y,z;
由总变形可以计算得到总应变:
Figure BDA0003050542990000172
εij——粒子i和粒子j的相对总应变矢量;
通过步骤6的计算,计算得到了每个粒子在不同时刻下的总应变张量,即得到了单个燃料元件的总变形状态;
步骤7:本算例中,由于温度的升高引起热膨胀,从而导致燃料元件发生变形,变形过程中,燃料元件先发生微小弹性变形。基于弹性模型计算弹性应力,弹性模型如公式(11)所示,
Figure BDA0003050542990000173
式中
Figure BDA0003050542990000174
——弹性应力张量的分量,N/m2
λ——拉梅常数的第一个参数,
Figure BDA0003050542990000175
μ——拉梅常数的第二个参数,
Figure BDA0003050542990000176
Figure BDA0003050542990000177
——弹性应变张量对角项之和,即
Figure BDA0003050542990000178
δαβ——克罗内克尔函数,
Figure BDA0003050542990000179
Figure BDA00030505429900001710
——弹性应变张量的分量;
其中记各向同性力:
Figure BDA00030505429900001711
式中
Pi E——粒子i的弹性各向同性力,N/m2
λ——拉梅常数的第一个参数,
Figure BDA00030505429900001712
Figure BDA00030505429900001713
——弹性应变张量对角项之和,即
Figure BDA00030505429900001714
d——空间维度;
n0——参考粒子数密度;
uij——粒子i和粒子j的相对总变形矢量,m;
rij——粒子i和粒子j的相对位置矢量,m;
Figure BDA0003050542990000181
——粒子i和粒子j的相对初始位置矢量,m;
Figure BDA0003050542990000182
——核函数,等同于
Figure BDA0003050542990000183
此外,在碰撞计算中,如果碰撞载荷Pi Contact小于弹性各向同性力Pi E,则令Pi Contact=Pi E;记各向异性力:
Figure BDA0003050542990000184
式中
Figure BDA0003050542990000185
——粒子i和粒子j间的弹性各向异性力,N/m2
μ——拉梅常数的第二个参数,
Figure BDA0003050542990000186
Figure BDA0003050542990000187
——粒子i和粒子j的相对弹性应变矢量;
由弹性各向同性力和弹性各向异性力分别求得对应的速度梯度,如公式(14)和(15)所示,
Figure BDA0003050542990000188
Figure BDA0003050542990000189
式中
vi——粒子i的速度矢量,m/s;
t——时间,s;
ρi——粒子i的密度,kg/m3
d——空间维度;
n0——参考粒子数密度;
rij——粒子i和粒子j的相对位置矢量,m;
Pi E——粒子i的弹性各向同性力,N/m2
Figure BDA00030505429900001810
——粒子j的弹性各向同性力,N/m2
Figure BDA0003050542990000191
——粒子i和粒子j间的弹性各向异性力,N/m2
Figure BDA0003050542990000192
——粒子i和粒子j的相对初始位置矢量,m;
Figure BDA0003050542990000193
——核函数,等同于
Figure BDA0003050542990000194
通过步骤7,模拟了单根燃料元件早期变形过程中的弹性变形;计算得到每个粒子在不同时刻下的弹性应变张量、弹性应力张量、弹性力所引起的加速度,即得到燃料元件的弹性应力应变随时间的变化过程;
步骤8:由于结构材料两端固定,发生热膨胀后,燃料包壳会向外侧弯曲,由此出现了固体的旋转或扭转运动;基于角动量定理计算粒子旋转角速度,粒子j对粒子i的切应力为:
Figure BDA0003050542990000195
式中
Fij——粒子j对粒子i的切应力,N;
d——空间维度;
n0——参考粒子数密度;
mi——粒子i的质量,kg;
ρi——粒子i的密度,kg/m3
Figure BDA0003050542990000196
——粒子i和粒子j间的弹性各向异性力在与rij垂直方向上的投影,N/m2
rij——粒子i和粒子j的相对位置矢量,m;
Figure BDA0003050542990000197
——粒子i和粒子j的相对初始位置矢量,m;
Figure BDA0003050542990000198
——核函数,等同于
Figure BDA0003050542990000199
粒子j对粒子i的切应力产生的扭矩为:
Tij=-rij×Fij 公式(17)
式中
Tij——粒子j对粒子i的切应力产生的扭矩,N·m;
rij——粒子i和粒子j的相对位置矢量,m;
Fij——粒子j对粒子i的切应力,N;
粒子i的角速度可由公式(18)得:
Figure BDA00030505429900001910
式中
ωi——粒子i的角速度,rad/s;
t——时间,s;
Ii——粒子i的转动惯量,kg·m2
Tij——粒子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)
Figure BDA0003050542990000201
式中
[ds]n——第n时间步下的应力增量张量,N/m2
μ——拉梅常数的第二个参数,
Figure BDA0003050542990000202
[dεe]n——第n时间步下的弹性应变增量张量;
[s]n——第n时间步下的应变张量,N/m2
[s]n-1——第n-1时间步下的应变张量,N/m2
Figure BDA0003050542990000203
——粒子i的等效应力,N/m2
sij——应力矢量,N/m2
通过步骤9的计算,得到了每个粒子在不同时刻下是否发生塑性变形的判定,即判定燃料元件的不同位置是否发生塑性变形;
步骤10:当等效应力大于屈服极限时,对粒子进行塑性应力应变计算,塑性应变增量如公式(22)和(23)所示,
Figure BDA0003050542990000211
[dε]n=[dεp]n+[dεe]n+[dεT]n 公式(23)
式中
[dsp]n——第n时间步下的塑性应力增量张量,N/m2
[ds]n——第n时间步下的应力增量张量,N/m2
[dsT]n——第n时间步下的热膨胀应力增量张量,N/m2
[dse]n——第n时间步下的弹性应力增量张量,N/m2
[s]n-1——第n-1时间步下的应变张量,N/m2
μ——拉梅常数的第二个参数,
Figure BDA0003050542990000212
n——增量参数,由材料的力学特性决定;
式中dλn的取值影响塑性应变和弹性应变的关系,不同特性的材料,其dλn的求解不同,以下针对线性弹塑性,dλn的求解如公式(24)所示,
Figure BDA0003050542990000213
式中
Figure BDA0003050542990000214
——假设所有应变均为弹性应变条件下的等效应力,N/m2
H′——塑性硬化系数;
Y——屈服极限,N/m2
Figure BDA0003050542990000215
——第n-1时间步下的等效塑性应变;
μ——拉梅常数的第二个参数,
Figure BDA0003050542990000216
由塑性模型获得塑性应变后,由公式(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数相同时才发生相互作用;初始化所有粒子的速度、位置、温度、应力张量和应变张量,分别为ui、ri、Ti、σi、εi;通过以上参数设定,建立燃料芯块、包壳与结构材料三者之间的固体边界,还原真实的燃料结构,并确定了燃料的初始位置、速度、温度、应力应变状态和相关物性参数;
步骤2:根据建立的初始粒子几何模型,为每个粒子建立固体内粒子邻居域,即存储每个粒子以自身为圆心,作用范围re为半径的区域内所有Object数与该粒子相同的粒子ID,固体内邻居域不随计算时间的推移而发生改变,即认为固体内连续介质发生相互作用的粒子始终保持不变;除了固体内粒子邻居域,还需要为所有1号粒子建立固体间粒子邻居域,即存储每个1号粒子以自身为圆心,作用范围re为半径的区域内所有Object数与该粒子不同的1号粒子ID,固体间邻居域会随着计算时间的推移而发生改变,即两个固体的相对位置会不断发生变化,当两者表面距离小于作用范围re时,两个固体的边界会发生相互作用,否则不会;定义粒子数密度为n,其值为粒子邻居域内所有粒子对中心粒子的核函数值之和,如公式(1)所示;
Figure FDA0003050542980000011
式中
ni——粒子i的粒子数密度;
w(rij,re)——核函数,采用的核函数形式为
Figure FDA0003050542980000012
rij——粒子i和粒子j的距离,m;
re——粒子作用范围,m;
下标i——任意粒子i;
下标j——粒子i的粒子邻居域内的任意粒子j;
由此,进一步定义固体内粒子数密度ninner为固体内粒子邻居域内所有粒子对中心粒子的核函数值之和、粒子间粒子数密度nouter为固体间粒子邻居域内所有粒子对中心粒子的核函数值之和、n0为参考粒子数密度,是假设粒子均匀分布情况下计算得到的粒子数密度;
步骤3:基于边界条件和固体碰撞模型计算固体边界粒子外力载荷,根据具体实际物理问题,为1号粒子设置外力的作用方向和大小,由公式(2)计算边界条件外力载荷的影响,边界条件外力载荷始终视为作用于粒子中心,因此不考虑该力所引起的粒子转动,公式采用张量表示法和下标表示法;
Figure FDA0003050542980000021
式中
vα——α方向上的速度分量,m/s;
Figure FDA0003050542980000022
——α方向上的边界条件外力载荷引起的加速度,m/s2
t——时间,s;
对于固体碰撞过程,采用固体碰撞模型计算1号粒子的碰撞载荷,如公式(3)所示,
Figure FDA0003050542980000023
式中
Pi Contact——粒子i的接触应力,N/m2
λ——拉梅常数的第一个参数,
Figure FDA0003050542980000024
γγ)i——粒子i的应变张量对角项之和,即εγγ=ε112233
nouter——粒子i的固体间邻居域的粒子数密度;
ninner——粒子i的固体内邻居域的粒子数密度;
E——杨氏模量,N/m2
υ——泊松比;
由碰撞载荷,计算对应的速度梯度,如公式(4)所示,
Figure FDA0003050542980000025
式中
vi——粒子i的速度矢量,m/s;
t——时间,s;
ρi——粒子i的密度,kg/m3
d——空间维度;
n0——参考粒子数密度;
rij——粒子i和粒子j的相对位置矢量,m;
Pi Contact——粒子i的接触应力,N/m2
Figure FDA0003050542980000031
——粒子j的接触应力,N/m2
Figure FDA0003050542980000032
——粒子i和粒子j的相对初始位置矢量,m;
Figure FDA0003050542980000033
——核函数,等同于
Figure FDA0003050542980000034
通过步骤3的计算,模拟核反应堆燃料早期变形过程中燃料元件受到的外力载荷和碰撞力,全面考虑了核反应堆燃料早期变形过程中可能发生的外力作用、燃料元件之间的碰撞、燃料元件与结构材料之间的碰撞、燃料包壳与芯块之间的碰撞、芯块与芯块之间的碰撞过程;计算得到每个粒子由外力载荷和碰撞引起的加速度,即燃料元件由外力载荷和碰撞引起的加速度;
步骤4:能量计算如公式(5)所示:
Figure FDA0003050542980000035
式中
h——焓值,J/kg;
t——时间,s;
ρ——密度,kg/m3
k——热导率,W/(m·K);
T——温度,K;
Qradiation——辐射热源,W/m3
Qout——外热源,W/m3
辐射换热采用斯忒藩-玻尔兹曼定律,如公式(6)所示
Figure FDA0003050542980000036
式中
Qradiation——辐射热源W/m3
ε——发射率;
σ——斯忒藩-玻耳兹曼常数;
T——粒子温度K;
Tenv——环境温度K;
l0——粒子直径m;
温度的拉普拉斯项采用拉普拉斯模型离散得到:
Figure FDA0003050542980000037
式中
d——维度;
n0——参考粒子数密度;
ρi——粒子i的密度,kg/m3
Figure FDA0003050542980000041
——当前时刻的粒子j温度K;
Ti k——当前时刻的粒子i温度K;
rij——粒子i和粒子j的相对位置矢量,m;
w(|rij|)——核函数,等同于w(|rij|,re);
λi——
Figure FDA0003050542980000042
由焓值计算粒子温度:
Figure FDA0003050542980000043
式中
T——温度K;
Ts——熔点K;
h——焓值J/kg;
hs0——开始熔化的焓值J/kg;
hs1——结束熔化的焓值J/kg;
cp——比热容J/(kg·K);
通过步骤4的计算,模拟核反应堆燃料早期变形过程中燃料元件内的导热、燃料元件间的辐射换热、芯块和包壳间的换热过程;计算得到每个粒子在不同时刻下的焓值和温度,即得到核反应堆燃料的焓值和温度随时间的变化过程;
步骤5:热膨胀应变计算如公式(8)所示,
Figure FDA0003050542980000044
式中
[dsT]i——粒子i的热膨胀应力增量张量,N/m2
κi——粒子i的热膨胀系数;
Ti——粒子i的温度,K;
Ti ref——粒子i的参考温度,K;
通过步骤5的计算,模拟核反应堆燃料早期变形过程中燃料元件的热膨胀变形过程;计算得到每个粒子在不同时刻下的热膨胀应变张量,即得到了燃料元件由于温度变化而发生碰撞变形的行为规律;
步骤6:计算总应变,固体内的相对变形是由当前相对位置减去旋转后的相对位置得到的,如公式(9)所示,
Figure FDA0003050542980000051
式中
uij——粒子i和粒子j的相对总变形矢量,m;
rij——粒子i和粒子j的相对位置矢量,m;
R——粒子i的旋转矩阵,R=RxRyRz
Figure FDA0003050542980000052
——粒子i和粒子j的相对初始位置矢量,m;
Rx——粒子i绕x轴旋转的旋转矩阵,
Figure FDA0003050542980000053
Ry——粒子i绕y轴旋转的旋转矩阵,
Figure FDA0003050542980000054
Rz——粒子i绕z轴旋转的旋转矩阵,
Figure FDA0003050542980000055
θα,ij——粒子i和粒子j绕α轴旋转的相对角度,°,
Figure FDA0003050542980000056
α可以取x,y,z;
由总变形可以计算得到总应变:
Figure FDA0003050542980000057
εij——粒子i和粒子j的相对总应变矢量;
通过步骤6的计算,计算得到了每个粒子在不同时刻下的总应变张量,即得到了核反应堆燃料元件的总变形状态;
步骤7:基于弹性模型计算弹性应力,弹性模型如公式(11)所示,
Figure FDA0003050542980000061
式中
Figure FDA0003050542980000062
——弹性应力张量的分量,N/m2
λ——拉梅常数的第一个参数,
Figure FDA0003050542980000063
μ——拉梅常数的第二个参数,
Figure FDA0003050542980000064
Figure FDA0003050542980000065
——弹性应变张量对角项之和,即
Figure FDA0003050542980000066
δαβ——克罗内克尔函数,
Figure FDA0003050542980000067
Figure FDA0003050542980000068
——弹性应变张量的分量;
其中记各向同性力:
Figure FDA0003050542980000069
式中
Pi E——粒子i的弹性各向同性力,N/m2
λ——拉梅常数的第一个参数,
Figure FDA00030505429800000610
Figure FDA00030505429800000611
——弹性应变张量对角项之和,即
Figure FDA00030505429800000612
d——空间维度;
n0——参考粒子数密度;
uij——粒子i和粒子j的相对总变形矢量,m;
rij——粒子i和粒子j的相对位置矢量,m;
Figure FDA00030505429800000613
——粒子i和粒子j的相对初始位置矢量,m;
Figure FDA00030505429800000614
——核函数,等同于
Figure FDA00030505429800000615
此外,在碰撞计算中,如果碰撞载荷Pi Contact小于弹性各向同性力Pi E,则令Pi Contact=Pi E
记各向异性力:
Figure FDA0003050542980000071
式中
Figure FDA0003050542980000072
——粒子i和粒子j间的弹性各向异性力,N/m2
μ——拉梅常数的第二个参数,
Figure FDA0003050542980000073
Figure FDA0003050542980000074
——粒子i和粒子j的相对弹性应变矢量;
由弹性各向同性力和弹性各向异性力分别求得对应的速度梯度,如公式(14)和(15)所示,
Figure FDA0003050542980000075
Figure FDA0003050542980000076
式中
vi——粒子i的速度矢量,m/s;
t——时间,s;
ρi——粒子i的密度,kg/m3
d——空间维度;
n0——参考粒子数密度;
rij——粒子i和粒子j的相对位置矢量,m;
Pi E——粒子i的弹性各向同性力,N/m2
Figure FDA0003050542980000077
——粒子j的弹性各向同性力,N/m2
Figure FDA0003050542980000078
——粒子i和粒子j间的弹性各向异性力,N/m2
Figure FDA0003050542980000079
——粒子i和粒子j的相对初始位置矢量,m;
Figure FDA00030505429800000710
——核函数,等同于
Figure FDA00030505429800000711
通过步骤7,模拟了核反应堆燃料早期变形过程中的弹性变形;计算得到每个粒子在不同时刻下的弹性应变张量、弹性应力张量、弹性力所引起的加速度,即得到燃料元件的弹性应力应变随时间的变化过程;
步骤8:基于角动量定理计算粒子旋转角速度,粒子j对粒子i的切应力为:
Figure FDA0003050542980000081
式中
Fij——粒子j对粒子i的切应力,N;
d——空间维度;
n0——参考粒子数密度;
mi——粒子i的质量,kg;
ρi——粒子i的密度,kg/m3
Figure FDA0003050542980000082
——粒子i和粒子j间的弹性各向异性力在与rij垂直方向上的投影,N/m2
rij——粒子i和粒子j的相对位置矢量,m;
Figure FDA0003050542980000083
——粒子i和粒子j的相对初始位置矢量,m;
Figure FDA0003050542980000084
——核函数,等同于
Figure FDA0003050542980000085
粒子j对粒子i的切应力产生的扭矩为:
Tij=-rij×Fij 公式(17)
式中
Tij——粒子j对粒子i的切应力产生的扭矩,N·m;
rij——粒子i和粒子j的相对位置矢量,m;
Fij——粒子j对粒子i的切应力,N;
粒子i的角速度由公式(18)得:
Figure FDA0003050542980000086
式中
ωi——粒子i的角速度,rad/s;
t——时间,s;
Ii——粒子i的转动惯量,kg·m2
Tij——粒子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)
Figure FDA0003050542980000091
式中
[ds]n——第n时间步下的应力增量张量,N/m2
μ——拉梅常数的第二个参数,
Figure FDA0003050542980000092
[dεe]n——第n时间步下的弹性应变增量张量;
[s]n——第n时间步下的应变张量,N/m2
[s]n-1——第n-1时间步下的应变张量,N/m2
Figure FDA0003050542980000093
——粒子i的等效应力,N/m2
sij——应力矢量,N/m2
通过步骤9的计算,得到了每个粒子在不同时刻下是否发生塑性变形的判定,即判定燃料元件的不同位置是否发生塑性变形;
步骤10:塑性应力应变计算,塑性应变增量如公式(22)和(23)所示,
Figure FDA0003050542980000094
[dε]n=[dεp]n+[dεe]n+[dεT]n 公式(23)
式中
[dsp]n——第n时间步下的塑性应力增量张量,N/m2
[ds]n——第n时间步下的应力增量张量,N/m2
[dsT]n——第n时间步下的热膨胀应力增量张量,N/m2
[dse]n——第n时间步下的弹性应力增量张量,N/m2
[s]n-1——第n-1时间步下的应变张量,N/m2
μ——拉梅常数的第二个参数,
Figure FDA0003050542980000101
n——增量参数,由材料的力学特性决定;
式中dλn的取值影响塑性应变和弹性应变的关系,不同特性的材料,其dλn的求解不同,以下针对线性弹塑性,dλn的求解如公式(24)所示,
Figure FDA0003050542980000102
式中
Figure FDA0003050542980000103
——假设所有应变均为弹性应变条件下的等效应力,N/m2
H′——塑性硬化系数;
Y——屈服极限,N/m2
Figure FDA0003050542980000104
——第n-1时间步下的等效塑性应变;
μ——拉梅常数的第二个参数,
Figure FDA0003050542980000105
由塑性模型获得塑性应变后,由公式(23)计算得到弹性应变,再采用步骤7中的弹性模型计算弹性应力;
通过步骤10的计算,模拟了核反应堆燃料早期变形过程中的塑性变形过程;计算得到了每个粒子在不同时刻下的塑性应变和弹性应变,即得到了燃料元件的塑性变形和弹性变形随时间的变化过程;
步骤11:断裂判断,当粒子间距大于或小于断裂阈值时,则认为发生断裂,断裂后断裂部分的2号内部粒子转化为1号边界粒子,且断裂粒子之间不再发生固体连续介质内的相互作用,仅作为固体边界参与边界载荷和碰撞的计算;根据外力载荷、碰撞和弹性力更新粒子速度、位置、角速度和旋转角度;
通过步骤11的计算,模拟了核反应堆燃料早期变形过程中的断裂过程;计算得到了每个粒子在不同时刻下的是否发生断裂的判定,即获得了燃料元件发生断裂的位置;
通过步骤1对实际的核反应堆燃料结构建立粒子几何模型,设定粒子的位置、速度、温度、物性参数和力学参数;通过步骤2存储粒子邻居域,为后续计算提供基础;通过步骤3至步骤10计算固体的温度和受力,受力包括外力、碰撞、弹性应力应变、塑性应变和热膨胀应变,更新粒子的速度和位置;通过步骤11判断粒子是否发生断裂;由步骤3至步骤10计算得到的加速度,更新所有粒子的速度和位置;对计算是否结束进行判定,若是,输出结果,若否,更新固体间粒子邻居域,再返回步骤3;综合以上步骤,模拟了核反应堆燃料早期变形过程,得到了变形过程中燃料元件的位置、速度、温度、应力状态、应变状态,通过以上数据对核反应堆燃料早期变形过程中的传热、热膨胀、弹性变形、塑性变形和受力情况展开机理性分析。
CN202110486476.0A 2021-04-30 2021-04-30 基于粒子法的核反应堆燃料早期变形分析方法 Active CN113191065B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110486476.0A CN113191065B (zh) 2021-04-30 2021-04-30 基于粒子法的核反应堆燃料早期变形分析方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110486476.0A CN113191065B (zh) 2021-04-30 2021-04-30 基于粒子法的核反应堆燃料早期变形分析方法

Publications (2)

Publication Number Publication Date
CN113191065A true CN113191065A (zh) 2021-07-30
CN113191065B CN113191065B (zh) 2022-07-26

Family

ID=76983646

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110486476.0A Active CN113191065B (zh) 2021-04-30 2021-04-30 基于粒子法的核反应堆燃料早期变形分析方法

Country Status (1)

Country Link
CN (1) CN113191065B (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115062525A (zh) * 2022-07-01 2022-09-16 西安交通大学 基于先进粒子法的核反应堆严重事故分析方法

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US4093682A (en) * 1976-02-14 1978-06-06 Hobeg Hochtemperaturreaktor - Brennelement Gmbh Process for the production of block-shaped fuel elements for high temperature reactors
CN110321641A (zh) * 2019-07-08 2019-10-11 西安交通大学 基于粒子法的熔融物与混凝土相互作用分析方法
CN110598324A (zh) * 2019-09-12 2019-12-20 西安交通大学 一种核反应堆弥散型板型燃料元件堆芯流固耦合计算方法
CN112102894A (zh) * 2020-09-04 2020-12-18 西安交通大学 基于粒子法的核反应堆堆芯材料熔池演变特性分析方法

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US4093682A (en) * 1976-02-14 1978-06-06 Hobeg Hochtemperaturreaktor - Brennelement Gmbh Process for the production of block-shaped fuel elements for high temperature reactors
CN110321641A (zh) * 2019-07-08 2019-10-11 西安交通大学 基于粒子法的熔融物与混凝土相互作用分析方法
CN110598324A (zh) * 2019-09-12 2019-12-20 西安交通大学 一种核反应堆弥散型板型燃料元件堆芯流固耦合计算方法
CN112102894A (zh) * 2020-09-04 2020-12-18 西安交通大学 基于粒子法的核反应堆堆芯材料熔池演变特性分析方法

Non-Patent Citations (6)

* Cited by examiner, † Cited by third party
Title
RONGHUA CHEN等: ""Numerical analysis of the melt behavior in a fuel support piece"", 《ELSEVIER》 *
YU JI等: ""Thermal performance optimization of a fuel element in particle bed reactors for nuclear thermal propulsion"", 《NUCLEAR ENGINEERING AND DESIGN》 *
丁淑蓉 等: ""弥散型燃料板辐照肿胀行为的有限元分析"", 《核动力工程》 *
段日强 等: ""采用无网格移动粒子法数值模拟闪蒸射流"", 《原子能科学技术》 *
谭军 等: ""Zr-4合金应力松弛过程中的热激活变形与动态应变时效"", 《金属学报》 *
陆维等: ""反应堆压力容器接管段三维弹塑性应力分析"", 《压力容器》 *

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115062525A (zh) * 2022-07-01 2022-09-16 西安交通大学 基于先进粒子法的核反应堆严重事故分析方法

Also Published As

Publication number Publication date
CN113191065B (zh) 2022-07-26

Similar Documents

Publication Publication Date Title
US11798658B2 (en) Multi-scale method for simulating mechanical behaviors of multiphase composite materials
Wang et al. An improved coupled thermo-mechanic bond-based peridynamic model for cracking behaviors in brittle solids subjected to thermal shocks
CN113191066B (zh) 基于无网格法的核反应堆燃料元件失效分析方法
Meguid et al. Development and validation of novel FE models for 3D analysis of peening of strain-rate sensitive materials
Zhu et al. Energy absorption of diamond lattice cylindrical shells under axial compression loading
CN115062525B (zh) 基于先进粒子法的核反应堆严重事故分析方法
CN113191065B (zh) 基于粒子法的核反应堆燃料早期变形分析方法
Liu et al. Multiphysics modeling of novel UO2-BeO sandwich fuel performance in a light water reactor
Yoo et al. Finite element analysis of the mechanical behavior of a nuclear fuel assembly spacer grid
Lei et al. Experimental model design and preliminary numerical verification of fluid–thermal–structural coupling problem
Sawei et al. Research progress on simulation modeling of metal foams
Liu et al. Reliability analysis of autonomous underwater vehicle aft pressure shell for optimal design and strength
Ma et al. A Hermite interpolation element-free Galerkin method for elasticity problems
Kim et al. Study of the mechanical properties and effects of particles for oxide dispersion strengthened Zircaloy-4 via a 3D representative volume element model
Jiang et al. Simulation of irradiation hardening of Zircaloy within plate-type dispersion nuclear fuel elements
Gan Thermo-mechanics of pebble beds in fusion blankets
Wang et al. Mechanical behaviors of the dispersion nuclear fuel plates induced by fuel particle swelling and thermal effect I: Effects of variations of the fuel particle volume fractions
Ding et al. Mechanical behaviors of the dispersion nuclear fuel plates induced by fuel particle swelling and thermal effect II: Effects of variations of the fuel particle diameters
Yan et al. Effect of fuel meat thickness on the non-uniform irradiation-induced thermo-mechanical behavior in monolithic UMo/Al fuel plates
Li et al. Sub-cell radiation-conduction (SCRC) model in densely packed beds of mono-sized pebbles for HTGR under vacuum condition
Hasanyan et al. On the buckling of a two-dimensional micropolar strip
Saibaba et al. A numerical study of the effect of extrusion parameters on the temperature distribution in Zr–2.5 Nb
Duan et al. Effect of strain gradients and heterogeneity on flow strength of particle reinforced metal-matrix composites
Wu et al. Validation of updated RANNS with effect of oxygen-dissolved metallic Zircaloy-4 under LOCA quench condition
Shmakov et al. A theoretical study of the kinetics of hydride cracking in zirconium alloys

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant