CN114743607B - 一种计算核燃料中裂变气体释放及辐照肿胀行为的方法 - Google Patents

一种计算核燃料中裂变气体释放及辐照肿胀行为的方法 Download PDF

Info

Publication number
CN114743607B
CN114743607B CN202210332746.7A CN202210332746A CN114743607B CN 114743607 B CN114743607 B CN 114743607B CN 202210332746 A CN202210332746 A CN 202210332746A CN 114743607 B CN114743607 B CN 114743607B
Authority
CN
China
Prior art keywords
bubble
grain boundary
nuclear fuel
gas
bubbles
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.)
Active
Application number
CN202210332746.7A
Other languages
English (en)
Other versions
CN114743607A (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 CN202210332746.7A priority Critical patent/CN114743607B/zh
Publication of CN114743607A publication Critical patent/CN114743607A/zh
Application granted granted Critical
Publication of CN114743607B publication Critical patent/CN114743607B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16C—COMPUTATIONAL CHEMISTRY; CHEMOINFORMATICS; COMPUTATIONAL MATERIALS SCIENCE
    • G16C10/00—Computational theoretical chemistry, i.e. ICT specially adapted for theoretical aspects of quantum chemistry, molecular mechanics, molecular dynamics or the like
    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16C—COMPUTATIONAL CHEMISTRY; CHEMOINFORMATICS; COMPUTATIONAL MATERIALS SCIENCE
    • G16C20/00—Chemoinformatics, i.e. ICT specially adapted for the handling of physicochemical or structural data of chemical particles, elements, compounds or mixtures
    • G16C20/10—Analysis or design of chemical reactions, syntheses or processes
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02E—REDUCTION OF GREENHOUSE GAS [GHG] EMISSIONS, RELATED TO ENERGY GENERATION, TRANSMISSION OR DISTRIBUTION
    • Y02E30/00—Energy generation of nuclear origin
    • Y02E30/30—Nuclear fission reactors

Landscapes

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

Abstract

本发明提出了一种计算核燃料中裂变气体释放及辐照肿胀行为的方法，主要包括以下步骤：第一性原理计算核燃料物性参数；基于速率理论获取核燃料晶界气泡初始分布；建立核燃料晶界气泡演化相场模型；计算渗流参数及燃料孔隙率；获得核燃料辐照肿胀规律；建立核燃料晶界网络结构；建立核燃料裂变气体释放行为相场‑渗流耦合框架；获得核燃料裂变气体释放随时间变化等信息；获取核燃料裂变气体在不同晶界分布下、不同温度下、不同裂变密度下的释放规律。本发明方法将相场理论和渗流理论结合，能够定量分析燃料微观结构对裂变气体释放阈值、释放份额等宏观行为的影响，为建立更加机理化的核燃料多尺度裂变气体释放模型奠定了技术基础。

Description

一种计算核燃料中裂变气体释放及辐照肿胀行为的方法
技术领域
本发明属于核燃料材料计算领域。具体涉及一种计算核燃料中裂变气体释放及辐照肿胀行为的方法。
背景技术
核燃料堪称核反应堆的“心脏”，是核反应堆系统的核心部件之一。我国在核燃料的运行方面仍然相对缺乏经验，而核燃料中的晶粒长大、气泡演化及肿胀行为尤为重要。裂变反应所释放的大量的裂变气体通过高温下原子的扩散运动在燃料中聚集并形成气泡。在陡峭的温度梯度的作用下，大量的气泡在晶界上不断聚集，最终在晶界上形成特有的通孔结构，进而造成燃料肿胀。同时，晶界上聚集的大量的气泡对高温下晶界的运动产生拖拽作用。因此，必须对核燃料中裂变气体释放引起的气泡演化行为及燃料肿胀进行深入研究。
建立在堆内实验研究结果基础上的描述燃料中气泡演化过程的经验模型具有一定的局限性。然而，作为实验研究的有力补充，理论模型和计算机模拟可以动态地再现核燃料中再结晶、气泡演化及燃料肿胀过程，从理论上阐明再结晶、气泡产生、聚集、长大的机理，并预测燃料肿胀性能的变化。目前仍然缺少对核燃料的这一裂变气体释放和气泡形貌演变行为进行研究和预测的模型。
核燃料在服役时表现出复杂的堆内行为，辐照裂变会使其产生裂变气体产物(Xe、Kr)和固体裂变产物(Cs、Rb等)，裂变气体扩散到晶界或孔隙中，在晶内、晶间形成气泡。随着燃耗的增加，气泡数量增加、不断迁移和长大，被晶界捕获，在晶界聚集、长大、联合、联网，形成释放通道，再通过裂纹释放出来，导致燃料性能下降，限制了燃料的寿命。因此，对核燃料面界中气泡生长研究，一方面有助于了解在役燃料的运行状态和使用寿命，及时的发现并解决问题；另一方面根据辐照特性，可采取适当的措施增强燃料的性能，进一步提高核电的经济效益。
发明内容
本发明目的在于提出一种计算核燃料中裂变气体释放及辐照肿胀行为的方法，该方法能够展现辐照条件下核燃料晶界气泡演变过程，能够定量分析燃料微观结构对裂变气体释放阈值、释放份额等宏观行为的影响，能够基于裂变气泡孔隙率等参数获得燃料辐照肿胀规律，为建立更加机理化的多尺度裂变气体释放模型奠定了技术基础。
基于上述目的，本发明采用以下技术方案：
S1：基于第一性原理，计算核燃料物性参数，所述核燃料物性参数包括：空位形成能、迁移能；气体原子形成能、迁移能；空位扩散系数以及气体原子扩散系数；
S2：基于步骤S1获取的核燃料物性参数，建立描述核燃料晶面气泡气体原子行为的速率理论模型和晶棱气泡气体原子行为的速率理论模型，获得核燃料晶界气泡的初始分布；
S3：通过建立核燃料多晶多相系统总自由能密度以及演化方程建立核燃料晶界气泡演化相场模型：定义相场浓度场变量、相场序参量；基于热力学理论获取核燃料基体相的自由能密度以及核燃料气泡相的自由能密度；引入插值函数，获得核燃料多晶多相体系的体自由能密度；结合界面梯度能、多晶相互作用能、晶界与气泡相互作用能得到核燃料多晶多相系统总自由能密度；建立演化方程；结合步骤S2获取的核燃料晶界气泡初始分布初始化核燃料晶界气泡演化相场模型；求解演化方程；
S4：根据步骤S3演化方程计算结果，获得晶界气泡形貌演化图，统计计算渗流参数及燃料孔隙率，所述渗流参数包括：晶界气泡密度、晶界气泡平均尺寸、晶界气泡覆盖率、气泡连通阈值；
S5：根据步骤S4计算所得渗流参数以及燃料孔隙率，获得核燃料辐照肿胀规律；
S6：建立核燃料晶界网络结构，所述核燃料晶界网络结构包括：沿核燃料径向分布多个晶界；晶面随机分布接触角为40°与80°之间的气泡；
S7：基于气体原子晶内扩散方程，并考虑重溶效应，建立核燃料裂变气体释放渗流模型，所述核燃料裂变气体释放渗流模型与步骤S3所述核燃料晶界气泡演化相场模型基于同一套输入参数，所述输入参数包括：温度、气泡接触角、裂变气体密度；根据步骤S4所述计算渗流参数的方法，获取不同温度、不同裂变密度、不同气泡接触角下的气泡连通阈值并保存；基于获取的气泡连通阈值并输入给核燃料裂变气体释放渗流模型，结合步骤S6所述核燃料晶界网络结构，建立核燃料裂变气体释放行为相场-渗流耦合框架；
S8：基于步骤S7建立的核燃料裂变气体释放行为相场-渗流耦合框架，获取核燃料裂变气体在晶界气孔网络的长程连通过程，获得裂变气体释放随时间的变化、不同时刻晶界气体浓度随半径的变化以及不同时刻晶界气体浓度分布信息；
S9：根据步骤S8获得的裂变气体释放随时间的变化、不同时刻晶界气体浓度随半径的变化以及不同时刻晶界气体浓度分布信息，结合晶界网络结构演化过程，获取不同温度下、不同晶界分布下、不同裂变密度下的核燃料裂变气体释放的规律。
所述步骤S1包括以下内容：
采用基于第一性原理平面波方法的Vasp软件；根据研究的核燃料，描述核燃料体系的原子组态，对完美单胞、包含空位的单胞以及包含Xe气体原子的单胞进行能量计算，从而得到空位、气体原子的形成能以及迁移能；基于能量计算结果，根据统计热力学内容，获取空位、气体原子的扩散系数；
所述空位形成能计算公式如下：
式中，ΔH(D,q,EF)表示含空位体系与无空位体系焓差即空位形成能；E(D,q)表示含空位体系总能量；E(perfect)表示无空位体系总能量；μi表示原子化学势，ni表示空位所涉及原子的数量，i表示空位所涉及原子的种类；q(EVBM+EF)表示空位电荷态，EVBM表示无空位体系价带顶，EF表示费米能级；
所述气体原子形成能计算公式如下：
式中，ΔH(Xe)表示含Xe气体原子体系与无空位体系焓差即气体原子形成能；E(Xe)表示含Xe气体原子体系总能量；μj表示气体原子化学势，nj表示气体原子的数量，j表示气体原子的种类；
计算空位、气体原子迁移过程中马鞍点能量与体系迁移过程稳态时的能量之差，即为迁移能；求得迁移能后，基于热力学理论，求得空位、气体原子扩散系数；所述扩散系数求解公式为：
空位扩散系数：
气体原子扩散系数：
其中，Dv、Dg分别表示空位、气体原子的扩散系数；D0表示扩散系数前因子，分别表示空位、气体原子的迁移能；kB是玻尔兹曼常数；T表示绝对温度。
所述步骤S2包括以下步骤：
S2.1：建立核燃料晶面气泡气体原子行为的速率理论模型，具体如下：
晶面气泡气体原子浓度的微分方程为：
其中，Cf表示晶面气泡内的气体原子浓度；Sgb表示单位体积晶界面积；Vg表示晶内气体原子速度；Cg表示晶内气体原子浓度；Vn表示含有n个原子的气泡速度；Cn表示晶内含有n个原子的气泡内气体原子浓度；Nn表示晶内含有n个原子的气泡内的平均原子数；dg表示晶粒直径；Dn表示晶内含有n个原子的气泡扩散系数；Vgb表示晶界移动速度；Nf表示晶面气泡内的平均原子数；Vf表示晶面气泡的速度；SNB表示节点的横截面积除以节点的体积；Ngf表示每个晶粒的晶面数目；PA表示晶面通道互连份额；δ表示晶界上的气泡重溶系数；b表示重溶常数；t表示时间；
晶面气泡气体原子浓度的微分方程左侧表示晶面气泡浓度变化率；右侧第一行两项分别表示晶界随机吸收气体原子和气泡、晶界定向吸收气体原子和气泡；第二行第一项表示晶界气泡脱离，第二项和第三项表示气泡定向迁移，第四项表示晶面气泡向晶棱迁移；第三行第一项表示晶界扫掠，第二项表示晶面气泡重溶；
当晶面上的气体浓度积累到一定条件后，晶面气泡会向晶棱迁移，根据晶面气泡互连的份额，计算由晶面到达晶棱的气体的量；
单位体积的晶面气泡覆盖率计算公式为：
Af＝πRf 2CfF(θ)
其中，Af表示单位体积晶面气泡覆盖晶面的比例；Rf表示晶面气泡的半径；F(θ)表示晶面凸镜形气泡的几何因子；
当单位体积晶面气泡覆盖晶面的比例与单位体积晶界面积之比超过判定阈值时，认为晶面通道形成，晶面气泡开始迁移至晶棱；该判定阈值与核燃料的属性、初始孔隙率以及晶粒尺寸相关；
晶面通孔形成后，裂变气体向晶棱迁移，根据晶面通道互连份额来确定晶面气泡到达晶棱的量，具体如下：
其中，PA表示晶面通道互连份额，Af *表示晶面通道连通判定阈值；σf表示尺度参数；A表示位置参数；当晶面气泡连通比率达到阈值时，认为晶面气体达到饱和，之后所有到达晶面的，裂变气体全部迁移至晶棱；
S2.2：建立核燃料晶棱气泡气体原子行为的速率理论模型，具体如下：
晶棱气泡气体原子浓度的微分方程为：
其中，Ce表示晶棱气泡内的原子浓度；PI表示晶棱通道互连份额；Ne表示晶棱含有n个原子的气泡内的平均原子数；Ngf表示每个晶粒的晶面数目；
晶棱气泡气体原子浓度的微分方程左侧表示晶棱气泡内气体原子浓度变化率；右侧第一行第一项表示晶面气泡通过晶界通道向晶棱迁移，第二项表示晶面向晶棱定向性迁移；第二行三项分别表示气泡脱离、气泡重溶和气体向自由空间释放；
晶棱通道互连份额表示如下：
其中，σe表示计算晶棱气泡的几何因子；Bvedge表示位置参数1；Bvpor表示位置参数2；Re表示晶棱气泡半径；
S2.3：基于步骤S2.1建立的核燃料晶面气泡气体原子行为的速率理论模型和步骤S2.2建立的核燃料晶棱气泡气体原子行为的速率理论模型，计算核燃料晶界气泡分布情况并将晶面、晶棱气泡分布信息存储，作为核燃料晶界气泡演化相场模型气泡分布初始化的依据。
所述步骤S3包括以下步骤：
S3.1：为研究核燃料晶界气泡演化过程，定义两个相场浓度场变量：空位浓度cv(r,t)、气体原子浓度cg(r,t)；定义一个相场序参量η(r,t)，用于区分气泡相与基体相；在气泡相中，相场序参量η(r,t)取值为1；在基体相中，相场序参量η(r,t)取值为0；为描述核燃料体系的多晶结构，定义了一系列相场序参量代表p个不同取向的晶粒；在第i个晶粒内有：/>
S3.2：基于热力学理论，推导核燃料基体相的自由能密度函数fm(cv,cg,T)有以下表达式：
式中，表示基体中空位浓度；/>表示基体中气体原子浓度；/>表示基体中空位平衡浓度；/>表示基体中气体原子平衡浓度；kB是玻尔兹曼常数；取值为1.3806505×10- 23J/K；T是绝对温度，单位为K；
其中，在推导自由能时，认为空位是去除材料粒子的晶格位置，而气体原子则占据取代晶格位；所以，以原子百分位为单位表示，气体原子浓度、空位浓度、完美晶格浓度相加为1；
基于热力学理论，推导核燃料气泡相的自由能密度函数fb(cv,cg,T)有以下表达式：
式中，表示气泡中空位浓度；/>表示气泡中气体原子浓度；/>表示气泡中空位平衡浓度；/>表示气泡中气体原子平衡浓度；/>表示气泡中气体原子最大浓度；
其中，在推导自由能时，认为气泡晶格位置只由材料粒子和空位占据，而气体原子只能占据空位晶格位置；
S3.3：引入插值函数，结合基于热力学推导的基体相和气泡相自由能，获得核燃料多晶多相体系的体自由能密度函数fbulk(cv,cg,η,T)：
fbulk(cv,cg,η,T)＝[1-h(η)]fm(cv,cg,T)+h(η)fb(cv,cg,T)
其中，fbulk(cv,cg,η,T)表示基体相自由能密度，fb(cv,cg,T)表示气泡相自由能密度；h(η)是构造的一个插值函数，其表达式为h(η)＝η3(6η2-15η+10)；插值函数取值满足：在基体相中，即η＝0.0时，h(η)＝0.0；在气泡相中，即η＝1.0时，h(η)＝1.0；
S3.4：根据步骤S3.3获得的核燃料多晶多相体系体自由能密度，结合界面梯度能、多晶相互作用能、晶界与气泡相互作用能，得到核燃料多晶多相系统总自由能密度：
式中，表示多晶相互作用自由能密度，其表达式为：
其中，A、B、aGB、as是唯象参数；
表示界面梯度能，其表达式为：
其中，κv、κg、κη、是梯度项系数；
S3.5：考虑辐照条件下空位、气体原子产生的演化方程为：
空位浓度场演化方程：
气体原子浓度场演化方程：
气泡相序参量演化方程：
多晶序参量演化方程：
其中，Mv、Mg分别表示空位和气体原子的迁移系数；L是自由界面的迁移率；ξv、ξg、ξη分别为空位、气体原子和气泡相的热涨落项；Pv(r,t)、Pg(r,t)分别表示辐照条件下空位以及气体原子的产生率；
S3.6：将步骤S2获取的核燃料晶界气泡初始分布导入核燃料晶界气泡演化相场模型，作为核燃料晶界气泡演化相场模型初始化的依据；
S3.7：采用空间上的有限差分法以及时间上的显示欧拉法对演化方程进行求解。
所述步骤S4包括以下步骤：
S4.1：根据演化方程计算结果，将计算结果中相场变量信息以vtk文件形式存储；vtk文件中写入模拟区域大小信息；使用paraview软件导入vtk文件进行结果可视化，获得晶界气泡形貌演化图；paraview软件使用流程为：点击“file”，选择需要导入的vtk文件；导入文件后，点击“Apply”，点击后，主界面依据vtk文件信息显示可视化结果，获得晶界气泡形貌图；
S4.2：根据演化方程计算结果，统计计算渗流参数方法如下：
晶界气泡密度、晶界气泡平均尺寸：晶界气泡密度定义为整个面晶界上气泡数量与模拟区域面晶界面积之比，单位为个/μm2；晶界气泡平均尺寸定义为：晶界上所有气泡面积之和与气泡数量之比，单位为μm2；采用Two-pass连通性分析的算法，对晶界气泡密度以及晶界气泡平均尺寸进行统计计算；Two-pass连通性分析的算法通过两遍扫描，将图像中存在的所有连通区域找出并标记；在第一次扫描时，从左至右、从上往下将所有区域内每个像素位置赋予一个label；扫描过程中同一个连通区域内的像素集合中可能被赋予一个或多个不同label；第二遍扫描时，需要将这些属于同一个连通区域但具有不同值的label合并，具体规则如下：
第一次扫描时：访问当前像素点，如果该label值等于1，则：
a.如果像素点领域中label值都为0，则赋予当前像素点一个新的label值：label+1；
b.如果像素点领域有label值>1的像素Neighbors，将领域中像素值最小的值赋予当前像素点；
c.记录Neighbors中各个label值之间的相等关系，即这些label值同属同一个连通区域；
第二次扫描时：访问当前像素点，如果label值大于1，找到与该label值同属相等关系的最小label值，把该值赋予当前像素点；
最后，统计区域内不同label值的个数，得到气泡数量；统计不同label值包含的区域面积，获得每个气泡的面积；
晶界气泡覆盖率：为获取气泡连通阈值，首先定义两个晶界气泡覆盖率的含义：晶界覆盖率：定义为晶界区域平面上所有气泡的投影面积除以晶界区域平面总面积；排放晶界覆盖率：定义为预设圆内接触圆周的气泡所占面积百分比；该预设圆假设为描绘晶界面的三重连接处，圆心位于晶界区域中心，半径为晶界区域边长的90％；采用以下算法实现参数计算：
a.扫描模拟区域中心晶界面，判断该区域上每个格点是否处于气泡中；
b.统计处于气泡中的格点个数，计算气泡总面积，气泡总面积比上晶界面面积即为晶界覆盖率；
c.根据不同label值得到的气泡分布，扫描整个区域，分别判断每个气泡是否与上述预设圆圆周接触；与预设圆圆周接触的所有气泡面积之和与气泡总面积之比即为排放晶界覆盖率；
气泡连通阈值：基于上述气泡晶界覆盖率的两个参数计算结果，画出排放晶界覆盖率随晶界覆盖率的变化曲线；排放晶界覆盖率随着晶界覆盖率增大而向1急速增长时，认为此时该晶界面气泡与晶界连通；取该急速增长过程中斜率最大点，该点处的晶界覆盖率即为气泡连通临界阈值；
S4.3：根据演化方程计算结果，燃料孔隙率计算方法如下：
孔隙率定义为：模拟中心晶界区域气泡所占体积分数；气泡体积统计方法为：
a.扫描整个模拟中心晶界区域，判断该区域中每个格点是否处于气泡中；
b.统计处于气泡中的格点个数，计算气泡总体积；气泡总体积比上模拟中心晶界区体积即为孔隙率。
所述步骤S5包括以下内容：
结合步骤S4中可视化结果，获得不同辐照强度下的燃料孔隙率随时间变化曲线、气泡平均尺寸变化曲线、气泡密度变化曲线；获得不同温度下的燃料孔隙率随时间变化曲线、气泡平均尺寸变化曲线、气泡密度变化曲线，得到核燃料辐照肿胀规律。
所述步骤S6包括以下步骤：
S6.1：基于核燃料芯块结构，晶粒结构采用四晶粒结构，晶界尺寸固定为10微米；
S6.2：基于轴对称假设，取燃料芯块r-z平面，在径向上设置500个晶界；晶界存在三种状态：闭合、开放、排气；每个晶界的初始状态都是闭合的；
S6.3：步骤S6.2中所述的500个晶界，每个晶界的气泡接触角在40°到80°之间随机分布。
所述步骤S7包括以下步骤：
S7.1：基于气体原子晶内扩散方程，并考虑重溶效应，建立核燃料裂变气体释放渗流模型；模型建立过程如下：进入晶面的气体流速用下式表示：
式中，N表示晶界气体浓度；f0表示晶内到达晶界的平均气体原子流速；δ表示晶界上的气泡重溶系数；b表示重溶常数；βg表示单位体积内核燃料内产生裂变气体原子的速率；
晶界网络的温度分布为：
其中，Tmax是燃料中心最高温度，Tmin是燃料外表面温度，r0是芯块外径；
单位体积内核燃料内产生裂变气体原子的速率：
βg＝-2.218×1018+3.854×1015(T)
气体原子扩散系数与温度的关系为：
其中，(a)适用于T<1381K，(b)适用于1381K<T<1650K，(c)适用于T>1650K；
S7.2：根据步骤S4所计算渗流参数，获取不同温度、不同裂变密度、不同气泡接触角下的气泡连通阈值，并按照输入变量顺序将计算结果排列并保存为txt文件；
S7.3：将步骤S6建立的核燃料晶界网络以及步骤S7.2获取的txt文件输入到步骤S7.1所述核燃料裂变气体释放渗流模型中，核燃料晶界气泡演化相场模型计算的渗流参数作为核燃料裂变气体释放渗流模型判断晶界是否开放的判据条件，建立了核燃料裂变气体释放行为相场-渗流耦合框架；所述判据条件为：当某个晶界的气体浓度达到连通阈值时，该晶界就被认为是连通状态；每个时间步，扫描晶界网络是否通过开放的晶界与燃料表面相连，若相连，则认为晶界内气体发生宏观排；在下一个时间步，排气的晶界被重置为关闭状态，晶界气体浓度重置为0。
所述步骤S8包括以下内容：
基于步骤S7建立的核燃料裂变气体释放行为相场-渗流耦合框架，获取核燃料裂变气体在晶界气孔网络的长程连通过程，获得裂变气体释放随时间的变化、不同时刻晶界气体浓度随半径的变化以及不同时刻晶界气体浓度分布信息。
所述步骤S9包括以下内容：
根据步骤S8获得的裂变气体释放随时间的变化、不同时刻晶界气体浓度随半径的变化以及不同时刻晶界气体浓度分布信息，结合晶界网络结构演化过程，获取不同温度下、不同晶界分布下、不同裂变密度下的核燃料裂变气体释放的规律。
和现有技术相比较，本发明具备如下优点：
(1)本发明提供的一种计算核燃料裂变气体释放及辐照肿胀行为的方法，将相场方法与速率理论模型有效地联系起来；传统的相场方法将初始气泡随机置于晶界处，而基于速率理论获取核燃料晶界气泡初始分布，有效提升了研究的科学性；本方法发挥了速率理论在研究气泡原子形核以及相场方法在处理非均匀场形貌演化方面的优势，获得了准确的、具有预测能力的核燃料裂变气体气泡形貌演变和辐照肿胀的介观尺度模型。
(2)本发明提供的一种计算核燃料裂变气体释放及辐照肿胀行为的方法，能够计算不同温度、不同裂变密度下核燃料裂变气体释放和辐照肿胀行为；能够计算不同晶粒结构核燃料核燃料裂变气体释放和辐照肿胀行为；能够动态再现晶间气泡生长、连通过程。
(3)本发明提供的一种计算核燃料裂变气体释放及辐照肿胀行为的方法，将相场方法和渗流理论有效地联系起来，建立了核燃料裂变气体释放相场-渗流耦合框架；传统的渗流方法在研究裂变气体释放行为时，采用经验值的连通阈值；而使用相场方法获取连通阈值并传递给渗流模型，有效提升了研究的准确性和精确性；本方法建立的核燃料裂变气体释放相场-渗流耦合框架，能够定量分析燃料微观结构对裂变气体释放阈值、释放份额等宏观行为的影响，为建立更加机理化的核燃料多尺度裂变气体释放模型奠定了技术基础。
附图说明
图1是本发明模拟流程图。
图2是接触角100°晶间气泡演化到稳态的形貌图。
图3是接触角120°晶间气泡演化到稳态的形貌图。
图4是接触角135°晶间气泡演化到稳态的形貌图。
图5是接触角160°晶间气泡演化到稳态的形貌图。
图6是气泡接触角理论值与计算值对比图。
图7是晶面气泡演化连通图。
图8是四晶结构晶间气泡演化连通图。
具体实施方式
下面结合附图和实例对本发明进行说明；以下实例用于说明本发明，但不能用来限制本发明的使用范围；
本实施方式提出一种核燃料中裂变气体释放及辐照肿胀行为的方法。该方法主要有：第一性原理计算核燃料物性参数；基于速率理论获取核燃料晶界气泡的初始分布；构建核燃料多晶多相系统总自由能密度；建立演化方程；输出演化方程结果并计算渗流参数；核燃料裂变气体释放行为相场-渗流耦合框架；输出核燃料裂变气体释放行为相场-渗流耦合框架计算结果并对结果进行分析处理，具体如下：
(1)第一性原理计算核燃料物性参数
采用基于第一性原理平面波方法的Vasp软件；根据研究的核燃料，描述核燃料体系的原子组态，对完美单胞、包含空位的单胞以及包含Xe气体原子的单胞进行能量计算，从而得到空位、气体原子的形成能以及迁移能；基于能量计算结果，根据统计热力学内容，获取空位、气体原子的扩散系数；
所述空位形成能计算公式如下：
式中，ΔH(D,q,EF)表示含空位体系与无空位体系焓差即空位形成能；E(D,q)表示含空位体系总能量；E(perfect)表示无空位体系总能量；μi表示原子化学势，ni表示空位所涉及原子的数量，i表示空位所涉及原子的种类；q(EVBM+EF)表示空位电荷态，EVBM表示无空位体系价带顶，EF表示费米能级；
所述气体原子形成能计算公式如下：
式中，ΔH(Xe)表示含Xe气体原子体系与无空位体系焓差即气体原子形成能；E(Xe)表示含Xe气体原子体系总能量；μj表示气体原子化学势，nj表示气体原子的数量，j表示气体原子的种类；
计算空位、气体原子迁移过程中马鞍点能量与体系迁移过程稳态时的能量之差，即为迁移能；求得迁移能后，基于热力学理论，求得空位、气体原子扩散系数；所述扩散系数求解公式为：
空位扩散系数：
气体原子扩散系数：
其中，Dv、Dg分别表示空位、气体原子的扩散系数；D0表示扩散系数前因子，分别表示空位、气体原子的迁移能；kB是玻尔兹曼常数；T表示绝对温度。
(2)基于速率理论获取核燃料晶界气泡的初始分布
晶面气泡气体原子浓度的微分方程为：
其中，Cf表示晶面气泡内的气体原子浓度；Sgb表示单位体积晶界面积；Vg表示晶内气体原子速度；Cg表示晶内气体原子浓度；Vn表示含有n个原子的气泡速度；Cn表示晶内含有n个原子的气泡内气体原子浓度；Nn表示晶内含有n个原子的气泡内的平均原子数；dg表示晶粒直径；Dn表示晶内含有n个原子的气泡扩散系数；Vgb表示晶界移动速度；Nf表示晶面气泡内的平均原子数；Vf表示晶面气泡的速度；SNB表示节点的横截面积除以节点的体积；Ngf表示每个晶粒的晶面数目；PA表示晶面通道互连份额；δ表示晶界上的气泡重溶系数；b表示重溶常数；t表示时间；
晶面气泡气体原子浓度的微分方程左侧表示晶面气泡浓度变化率；右侧第一行两项分别表示晶界随机吸收气体原子和气泡、晶界定向吸收气体原子和气泡；第二行第一项表示晶界气泡脱离，第二项和第三项表示气泡定向迁移，第四项表示晶面气泡向晶棱迁移；第三行第一项表示晶界扫掠，第二项表示晶面气泡重溶；
当晶面上的气体浓度积累到一定条件后，晶面气泡会向晶棱迁移，根据晶面气泡互连的份额，计算由晶面到达晶棱的气体的量；
单位体积的晶面气泡覆盖率计算公式为：
Af＝πRf 2CfF(θ)
其中，Af表示单位体积晶面气泡覆盖晶面的比例；Rf表示晶面气泡的半径；F(θ)表示晶面凸镜形气泡的几何因子；
当单位体积晶面气泡覆盖晶面的比例与单位体积晶界面积之比超过判定阈值时，认为晶面通道形成，晶面气泡开始迁移至晶棱；该判定阈值与核燃料的属性、初始孔隙率以及晶粒尺寸相关；
晶面通孔形成后，裂变气体向晶棱迁移，根据晶面通道互连份额来确定晶面气泡到达晶棱的量，具体如下：
其中，PA表示晶面通道互连份额，Af *表示晶面通道连通判定阈值；σf表示尺度参数；A表示位置参数；当晶面气泡连通比率达到阈值时，认为晶面气体达到饱和，之后所有到达晶面的，裂变气体全部迁移至晶棱；
晶棱气泡气体原子浓度的微分方程为：
其中，Ce表示晶棱气泡内的原子浓度；PI表示晶棱通道互连份额；Ne表示晶棱含有n个原子的气泡内的平均原子数；Ngf表示每个晶粒的晶面数目；
晶棱气泡气体原子浓度的微分方程左侧表示晶棱气泡内气体原子浓度变化率；右侧第一行第一项表示晶面气泡通过晶界通道向晶棱迁移，第二项表示晶面向晶棱定向性迁移；第二行三项分别表示气泡脱离、气泡重溶和气体向自由空间释放；
晶棱通道互连份额表示如下：
其中，σe表示计算晶棱气泡的几何因子；Bvedge表示位置参数1；Bvpor表示位置参数2；Re表示晶棱气泡半径；
基于上述模型，计算核燃料晶界气泡分布情况并将晶面、晶棱气泡分布信息存储，作为核燃料晶界气泡演化相场模型气泡分布初始化的依据。
(3)构建核燃料多晶多相系统总自由能密度
定义两个相场浓度场变量：空位浓度cv(r,t)、气体原子浓度cg(r,t)；定义一个相场序参量η(r,t)，用于区分气泡相与基体相；在气泡相中，相场序参量η(r,t)取值为1；在基体相中，相场序参量η(r,t)取值为0；为描述核燃料体系的多晶结构，定义了一系列相场序参量代表p个不同取向的晶粒；在第i个晶粒内有：/>
基于热力学理论，推导核燃料基体相的自由能密度函数fm(cv,cg,T)有以下表达式：
式中，表示基体中空位浓度；/>表示基体中气体原子浓度；/>表示基体中空位平衡浓度；/>表示基体中气体原子平衡浓度；kB是玻尔兹曼常数；取值为1.3806505×10- 23J/K；T是绝对温度，单位为K；
其中，在推导自由能时，认为空位是去除材料粒子的晶格位置，而气体原子则占据取代晶格位；所以，以原子百分位为单位表示，气体原子浓度、空位浓度、完美晶格浓度相加为1；
基于热力学理论，推导核燃料气泡相的自由能密度函数fb(cv,cg,T)有以下表达式：
式中，表示气泡中空位浓度；/>表示气泡中气体原子浓度；/>表示气泡中空位平衡浓度；/>表示气泡中气体原子平衡浓度；/>表示气泡中气体原子最大浓度；
其中，在推导自由能时，认为气泡晶格位置只由材料粒子和空位占据，而气体原子只能占据空位晶格位置；
引入插值函数，结合基于热力学推导的基体相和气泡相自由能，获得核燃料多晶多相体系的体自由能密度函数fbulk(cv,cg,η,T)：
fbulk(cv,cg,η,T)＝[1-h(η)]fm(cv,cg,T)+h(η)fb(cv,cg,T)
其中，fbulk(cv,cg,η,T)表示基体相自由能密度，fb(cv,cg,T)表示气泡相自由能密度；h(η)是构造的一个插值函数，其表达式为h(η)＝η3(6η2-15η+10)；插值函数取值满足：在基体相中，即η＝0.0时，h(η)＝0.0；在气泡相中，即η＝1.0时，h(η)＝1.0；
结合界面梯度能、多晶相互作用能、晶界与气泡相互作用能，得到核燃料多晶多相系统总自由能密度：
式中，表示多晶相互作用自由能密度，其表达式为：
其中，A、B、aGB、as是唯象参数；
表示界面梯度能，其表达式为：
其中，κv、κg、κη、是梯度项系数。
(4)建立演化方程
考虑辐照条件下空位、气体原子产生的演化方程为：
空位浓度场演化方程：
气体原子浓度场演化方程：
气泡相序参量演化方程：
多晶序参量演化方程：
其中，Mv、Mg分别表示空位和气体原子的迁移系数；L是自由界面的迁移率；ξv、ξg、ξη分别为空位、气体原子和气泡相的热涨落项；Pv(r,t)、Pg(r,t)分别表示辐照条件下空位以及气体原子的产生率。
(5)输出演化方程计算结果并计算渗流参数
根据演化方程计算结果，将计算结果中相场变量信息以vtk文件形式存储；vtk文件中写入模拟区域大小信息；使用paraview软件导入vtk文件进行结果可视化，获得晶界气泡形貌演化图；
根据演化方程计算结果，统计计算渗流参数方法如下：
晶界气泡密度、晶界气泡平均尺寸：晶界气泡密度定义为整个面晶界上气泡数量与模拟区域面晶界面积之比，单位为个/μm2；晶界气泡平均尺寸定义为：晶界上所有气泡面积之和与气泡数量之比，单位为μm2；采用Two-pass连通性分析的算法，对晶界气泡密度以及晶界气泡平均尺寸进行统计计算；Two-pass连通性分析的算法通过两遍扫描，将图像中存在的所有连通区域找出并标记；在第一次扫描时，从左至右、从上往下将所有区域内每个像素位置赋予一个label；扫描过程中同一个连通区域内的像素集合中可能被赋予一个或多个不同label；第二遍扫描时，需要将这些属于同一个连通区域但具有不同值的label合并，具体规则如下：
第一次扫描时：访问当前像素点，如果该label值等于1，则：
a.如果像素点领域中label值都为0，则赋予当前像素点一个新的label值：label+1；
b.如果像素点领域有label值>1的像素Neighbors，将领域中像素值最小的值赋予当前像素点；
c.记录Neighbors中各个label值之间的相等关系，即这些label值同属同一个连通区域；
第二次扫描时：访问当前像素点，如果label值大于1，找到与该label值同属相等关系的最小label值，把该值赋予当前像素点；
最后，统计区域内不同label值的个数，得到气泡数量；统计不同label值包含的区域面积，获得每个气泡的面积；
晶界气泡覆盖率：为获取气泡连通阈值，首先定义两个晶界气泡覆盖率的含义：晶界覆盖率：定义为晶界区域平面上所有气泡的投影面积除以晶界区域平面总面积；排放晶界覆盖率：定义为预设圆内接触圆周的气泡所占面积百分比；该预设圆假设为描绘晶界面的三重连接处，圆心位于晶界区域中心，半径为晶界区域边长的90％；采用以下算法实现参数计算：
a.扫描模拟区域中心晶界面，判断该区域上每个格点是否处于气泡中；
b.统计处于气泡中的格点个数，计算气泡总面积，气泡总面积比上晶界面面积即为晶界覆盖率；
c.根据不同label值得到的气泡分布，扫描整个区域，分别判断每个气泡是否与上述预设圆圆周接触；与预设圆圆周接触的所有气泡面积之和与气泡总面积之比即为排放晶界覆盖率；
气泡连通阈值：基于上述气泡晶界覆盖率的两个参数计算结果，画出排放晶界覆盖率随晶界覆盖率的变化曲线；排放晶界覆盖率随着晶界覆盖率增大而向1急速增长时，认为此时该晶界面气泡与晶界连通；取该急速增长过程中斜率最大点，该点处的晶界覆盖率即为气泡连通临界阈值；
根据演化方程计算结果，燃料孔隙率计算方法如下：
孔隙率定义为：模拟中心晶界区域气泡所占体积分数；气泡体积统计方法为：
a.扫描整个模拟中心晶界区域，判断该区域中每个格点是否处于气泡中；
b.统计处于气泡中的格点个数，计算气泡总体积；气泡总体积比上模拟中心晶界区体积即为孔隙率；
基于上述计算方法，获取不同温度、不同裂变密度、不同气泡接触角下的渗流参数和孔隙率。
(6)建立核燃料裂变气体释放行为相场-渗流耦合框架
基于气体原子晶内扩散方程，并考虑重溶效应，进入晶面的气体流速用下式表示：
式中，N表示晶界气体浓度；f0表示晶内到达晶界的平均气体原子流速；δ表示晶界上的气泡重溶系数；b表示重溶常数；βg表示单位体积内核燃料内产生裂变气体原子的速率；
晶界网络的温度分布为：
其中，Tmax是燃料中心最高温度，Tmin是燃料外表面温度，r0是芯块外径；
单位体积内核燃料内产生裂变气体原子的速率：
βg＝-2.218×1018+3.85×1015(T)
气体原子扩散系数与温度的关系为：
其中，(a)适用于T<1381K，(b)适用于1381K<T<1650K，(c)适用于T>1650K；
根据步骤S4所计算的渗流参数，获取不同温度、不同裂变密度、不同气泡接触角下的气泡连通阈值，并按照输入变量顺序将计算结果排列并保存为txt文件；
将步骤S6建立的核燃料晶界网络以及步骤S7.2获取的txt文件输入到步骤S7.1所述核燃料裂变气体释放渗流模型中，核燃料晶界气泡演化相场模型计算的渗流参数作为核燃料裂变气体释放渗流模型判断晶界是否开放的判据条件，建立了核燃料裂变气体释放行为相场-渗流耦合框架；所述判据条件为：当某个晶界的气体浓度达到连通阈值时，该晶界就被认为是连通状态；每个时间步，扫描晶界网络是否通过开放的晶界与燃料表面相连，若相连，则认为晶界内气体发生宏观排；在下一个时间步，排气的晶界被重置为关闭状态晶界气体浓度重置为0。
(7)输出核燃料裂变气体释放行为相场-渗流耦合框架计算结果并对结果进行分析处理
基于燃料裂变气体释放行为相场-渗流耦合框架计算结果，取核燃料裂变气体在晶界气孔网络的长程连通过程，获得裂变气体释放随时间的变化、不同时刻晶界气体浓度随半径的变化以及不同时刻晶界气体浓度分布信息；根据获得的裂变气体释放随时间的变化、不同时刻晶界气体浓度随半径的变化以及不同时刻晶界气体浓度分布信息，结合晶界网络结构演化过程，获取不同温度下、不同晶界分布下、不同裂变密度下的核燃料裂变气体释放的规律。
实施示例一
本实例按照具体实施方式中的流程，使用UO2燃料参数，对UO2燃料中单个晶界气泡的演化进行了气泡接触角角度验证；核燃料多晶多相系统的总自由能表达式：
其中，表示多晶相互作用自由能密度，其表达式为：
理论上，式中参数aGB与as与晶界能γGB、界面能γs的关系为：
而气泡接触角/>与能量的关系为：/>
核燃料多晶多相系统的演化方程为：
空位浓度场演化方程：
气体原子浓度场演化方程：
气泡相序参量演化方程：
多晶序参量演化方程：
其中，Mv、Mg分别表示空位和气体原子的迁移系数；L是自由界面的迁移率；ξv、ξg、ξη分别为空位、气体原子和气泡相的热涨落项；Pv(r,t)、Pg(r,t)分别表示辐照条件下空位以及气体原子的产生率。在本实例中，空位以及气体原子的产生率为0；
通过空间上有限差分法以及时间上显示欧拉法求解演化方程，并利用Paraview对数值结果进行可视化处理；模拟结果中单个晶界气泡演化到稳态的二维中心截面形貌分布如图2、图3、图4、图5所示，气泡接触角计算值与理论值对比结果如图6所示，计算结果与理论分析一致。
实施示例二
本实例在实例一的基础上对多晶结构辐照条件下UO2晶间气泡演化进行了相场模拟研究；
晶面气泡气体原子浓度的微分方程为：
当晶面上的气体浓度积累到一定条件后，晶面气泡会向晶棱迁移，根据晶面气泡互连的份额，计算由晶面到达晶棱的气体的量；
单位体积的晶面气泡覆盖率计算公式为：
Af＝πRf 2CfF(θ)
晶面通孔形成后，裂变气体向晶棱迁移，根据晶面通道互连的比率来确定晶面气泡到达晶棱的量，具体如下：
晶棱气泡气体原子浓度的微分方程为：
晶棱通道互连份额表示如下：
通过有限差分法求解微分方程，得到气泡、气体原子浓度分布，获取核燃料晶间气泡初始分布；
核燃料多晶多相系统的总自由能表达式：
式中，表示多晶相互作用自由能密度，其表达式为：
系统的演化方程为：
空位浓度场演化方程：
气体原子浓度场演化方程：
气泡相序参量演化方程：
多晶序参量演化方程：
通过有限差分法求解演化方程，并利用Paraview对数值结果进行可视化处理；模拟结果中晶面气泡连通的形貌分布如图7所示，四晶结构晶间气泡连通的形貌分布如图8所示，计算结果与实验观测到晶间气泡连通和晶界三重交叉处出现三角状气泡的现象一致。
实施示例三
本实例基于实例二计算结果，获取获取不同温度、不同裂变密度、不同气泡接触角下的渗流参数，利用核燃料裂变气体释放行为相场-渗流耦合框架对UO2燃料裂变气体释放行进行了研究；
基于晶内扩散方程，并考虑重溶效应，进入晶面的气体流速用下式表示：
晶界网络的温度分布为：
其中，Tmax是燃料中心最高温度，Tmin是燃料外表面温度，r0是芯块外径；
气体产生率为：
βg＝-2.2181018+3.854×1015 (T)
气体原子扩散系数与温度的关系为：
其中，(a)适用于T<1381K，(b)适用于1381K<T<1650K，(c)适用于T>1650K；
计算出的渗流参数满足：不同温度、不同接触角、不同裂变密度下晶间气泡连通阈值在0.5-0.6区间内；基于开源框架OPENPNM，对UO2燃料中裂变气体长程连通过程进行了模拟；计算结果显示：在排放团簇出现前，晶界气体不断积累，裂变气体在模拟前期没有释放过程，在开始释放裂变气体后，随时间呈线性变化；芯块中心高温区域积累最快，晶界气体浓度随半径减小而增大；随着晶界进入饱和状态，开始出现与自由空腔连通的排放团簇，芯块中心处晶界气体首先从与芯块上表面连通的孔道释放，晶界气体浓度断崖式下降，因此气体浓度沿径向呈现两端低、中间高的“山峰形”分布，且随着燃耗加深，“峰顶”逐渐向芯块外侧移动；大约在1000天，排放团簇出现，大约在4000天，连接到燃料表面的释放通道形成。

Claims (8)

1.一种计算核燃料中裂变气体释放及辐照肿胀行为的方法，其特征在于：包括以下步骤：
S1：基于第一性原理，计算核燃料物性参数，所述核燃料物性参数包括：空位形成能、迁移能；气体原子形成能、迁移能；空位扩散系数以及气体原子扩散系数；
S2：基于步骤S1获取的核燃料物性参数，建立描述核燃料晶面气泡气体原子行为的速率理论模型和晶棱气泡气体原子行为的速率理论模型，获得核燃料晶界气泡的初始分布；
S3：通过建立核燃料多晶多相系统总自由能密度以及演化方程建立核燃料晶界气泡演化相场模型：定义相场浓度场变量、相场序参量；基于热力学理论获取核燃料基体相的自由能密度以及核燃料气泡相的自由能密度；引入插值函数，获得核燃料多晶多相体系的体自由能密度；结合界面梯度能、多晶相互作用能、晶界与气泡相互作用能得到核燃料多晶多相系统总自由能密度；建立演化方程；结合步骤S2获取的核燃料晶界气泡初始分布初始化核燃料晶界气泡演化相场模型；求解演化方程；
S4：根据步骤S3演化方程计算结果，获得晶界气泡形貌演化图，统计计算渗流参数及燃料孔隙率，所述渗流参数包括：晶界气泡密度、晶界气泡平均尺寸、晶界气泡覆盖率、气泡连通阈值；
S5：根据步骤S4计算所得渗流参数以及燃料孔隙率，获得核燃料辐照肿胀规律；
S6：建立核燃料晶界网络结构，所述核燃料晶界网络结构包括：沿核燃料径向分布多个晶界；晶面随机分布接触角为40°与80°之间的气泡；
S7：基于气体原子晶内扩散方程，并考虑重溶效应，建立核燃料裂变气体释放渗流模型，所述核燃料裂变气体释放渗流模型与步骤S3所述核燃料晶界气泡演化相场模型基于同一套输入参数，所述输入参数包括：温度、气泡接触角、裂变气体密度；根据步骤S4所述计算渗流参数的方法，获取不同温度、不同裂变密度、不同气泡接触角下的气泡连通阈值并保存；基于获取的气泡连通阈值并输入给核燃料裂变气体释放渗流模型，结合步骤S6所述核燃料晶界网络结构，建立核燃料裂变气体释放行为相场-渗流耦合框架；
S8：基于步骤S7建立的核燃料裂变气体释放行为相场-渗流耦合框架，获取核燃料裂变气体在晶界气孔网络的长程连通过程，获得裂变气体释放随时间的变化、不同时刻晶界气体浓度随半径的变化以及不同时刻晶界气体浓度分布信息；
S9：根据步骤S8获得的裂变气体释放随时间的变化、不同时刻晶界气体浓度随半径的变化以及不同时刻晶界气体浓度分布信息，结合晶界网络结构演化过程，获取不同温度下、不同晶界分布下、不同裂变密度下的核燃料裂变气体释放的规律。
2.根据权利要求1所述的计算核燃料中裂变气体释放及辐照肿胀行为的方法，其特征在于：步骤S1中，所述第一性原理计算核燃料物性参数的方法为：
采用基于第一性原理平面波方法的Vasp软件；根据研究的核燃料，描述核燃料体系的原子组态，对完美单胞、包含空位的单胞以及包含Xe气体原子的单胞进行能量计算，从而得到空位、气体原子的形成能以及迁移能；基于能量计算结果，根据统计热力学内容，获取空位、气体原子的扩散系数；
所述空位形成能计算公式如下：
式中，ΔH(D,q,EF)表示含空位体系与无空位体系焓差即空位形成能；E(D,q)表示含空位体系总能量；E(perfect)表示无空位体系总能量；μi表示原子化学势，ni表示空位所涉及原子的数量，i表示空位所涉及原子的种类；q(EVBM+EF)表示空位电荷态，EVBM表示无空位体系价带顶，EF表示费米能级；
所述气体原子形成能计算公式如下：
式中，ΔH(Xe)表示含Xe气体原子体系与无空位体系焓差即气体原子形成能；E(Xe)表示含Xe气体原子体系总能量；μj表示气体原子化学势，nj表示气体原子的数量，j表示气体原子的种类；
计算空位、气体原子迁移过程中马鞍点能量与体系迁移过程稳态时的能量之差，即为迁移能；求得迁移能后，基于热力学理论，求得空位、气体原子扩散系数；所述扩散系数求解公式为：
空位扩散系数：
气体原子扩散系数：
其中，Dv、Dg分别表示空位、气体原子的扩散系数；D0表示扩散系数前因子，分别表示空位、气体原子的迁移能；kB是玻尔兹曼常数；T表示绝对温度。
3.根据权利要求1所述的计算核燃料中裂变气体释放及辐照肿胀行为的方法，其特征在于：步骤S2中，包括以下步骤：
S2.1：建立核燃料晶面气泡气体原子行为的速率理论模型，具体如下：
晶面气泡气体原子浓度的微分方程为：
其中，Cf表示晶面气泡内的气体原子浓度；Sgb表示单位体积晶界面积；Vg表示晶内气体原子速度；Cg表示晶内气体原子浓度；Vn表示含有n个原子的气泡速度；Cn表示晶内含有n个原子的气泡内气体原子浓度；Nn表示晶内含有n个原子的气泡内的平均原子数；dg表示晶粒直径；Dn表示晶内含有n个原子的气泡扩散系数；Vgb表示晶界移动速度；Nf表示晶面气泡内的平均原子数；Vf表示晶面气泡的速度；SNB表示节点的横截面积除以节点的体积；Ngf表示每个晶粒的晶面数目；PA表示晶面通道互连份额；δ表示晶界上的气泡重溶系数；b表示重溶常数；t表示时间；
晶面气泡气体原子浓度的微分方程左侧表示晶面气泡浓度变化率；右侧第一行两项分别表示晶界随机吸收气体原子和气泡、晶界定向吸收气体原子和气泡；第二行第一项表示晶界气泡脱离，第二项和第三项表示气泡定向迁移，第四项表示晶面气泡向晶棱迁移；第三行第一项表示晶界扫掠，第二项表示晶面气泡重溶；
当晶面上的气体浓度积累到一定条件后，晶面气泡会向晶棱迁移，根据晶面气泡互连的份额，计算由晶面到达晶棱的气体的量；
单位体积的晶面气泡覆盖率计算公式为：
Af＝πRf 2CfF(θ)
其中，Af表示单位体积晶面气泡覆盖晶面的比例；Rf表示晶面气泡的半径；F(θ)表示晶面凸镜形气泡的几何因子；
当单位体积晶面气泡覆盖晶面的比例与单位体积晶界面积之比超过判定阈值时，认为晶面通道形成，晶面气泡开始迁移至晶棱；该判定阈值与核燃料的属性、初始孔隙率以及晶粒尺寸相关；
晶面通孔形成后，裂变气体向晶棱迁移，根据晶面通道互连份额来确定晶面气泡到达晶棱的量，具体如下：
其中，PA表示晶面通道互连份额，Af *表示晶面通道连通判定阈值；σf表示尺度参数；A表示位置参数；当晶面气泡连通比率达到阈值时，认为晶面气体达到饱和，之后所有到达晶面的，裂变气体全部迁移至晶棱；
S2.2：建立核燃料晶棱气泡气体原子行为的速率理论模型，具体如下：
晶棱气泡气体原子浓度的微分方程为：
其中，Ce表示晶棱气泡内的原子浓度；PI表示晶棱通道互连份额；Ne表示晶棱含有n个原子的气泡内的平均原子数；Ngf表示每个晶粒的晶面数目；
晶棱气泡气体原子浓度的微分方程左侧表示晶棱气泡内气体原子浓度变化率；右侧第一行第一项表示晶面气泡通过晶界通道向晶棱迁移，第二项表示晶面向晶棱定向性迁移；第二行三项分别表示气泡脱离、气泡重溶和气体向自由空间释放；
晶棱通道互连份额表示如下：
其中，σe表示计算晶棱气泡的几何因子；Bvedge表示位置参数1；Bvpor表示位置参数2；Re表示晶棱气泡半径；
S2.3：基于步骤S2.1建立的核燃料晶面气泡气体原子行为的速率理论模型和步骤S2.2建立的核燃料晶棱气泡气体原子行为的速率理论模型，计算核燃料晶界气泡分布情况并将晶面、晶棱气泡分布信息存储，作为核燃料晶界气泡演化相场模型气泡分布初始化的依据。
4.根据权利要求1所述的计算核燃料中裂变气体释放及辐照肿胀行为的方法，其特征在于：步骤S3中，包括以下步骤：
S3.1：为研究核燃料晶界气泡演化过程，定义两个相场浓度场变量：空位浓度cv(r,t)、气体原子浓度cg(r,t)；定义一个相场序参量η(r,t)，用于区分气泡相与基体相；在气泡相中，相场序参量η(r,t)取值为1；在基体相中，相场序参量η(r,t)取值为0；为描述核燃料体系的多晶结构，定义了一系列相场序参量i＝1→p，代表p个不同取向的晶粒；在第i个晶粒内有：/> j≠i；
S3.2：基于热力学理论，推导核燃料基体相的自由能密度函数fm(cv,cg,T)有以下表达式：
式中，表示基体中空位浓度；/>表示基体中气体原子浓度；/>表示基体中空位平衡浓度；/>表示基体中气体原子平衡浓度；kB是玻尔兹曼常数；取值为1.3806505×10-23J/K；T是绝对温度，单位为K；
其中，在推导自由能时，认为空位是去除材料粒子的晶格位置，而气体原子则占据取代晶格位；所以，以原子百分位为单位表示，气体原子浓度、空位浓度、完美晶格浓度相加为1；
基于热力学理论，推导核燃料气泡相的自由能密度函数fb(cv,cg,T)有以下表达式：
式中，表示气泡中空位浓度；/>表示气泡中气体原子浓度；/>表示气泡中空位平衡浓度；/>表示气泡中气体原子平衡浓度；/>表示气泡中气体原子最大浓度；
其中，在推导自由能时，认为气泡晶格位置只由材料粒子和空位占据，而气体原子只能占据空位晶格位置；
S3.3：引入插值函数，结合基于热力学推导的基体相和气泡相自由能，获得核燃料多晶多相体系的体自由能密度函数fbulk(cv,cg,η,T)：
fbulk(cv,cg,η,T)＝[1-h(η)]fm(cv,cg,T)+h(η)fb(cv,cg,T)
其中，fbulk(cv,cg,η,T)表示基体相自由能密度，fb(cv,cg,T)表示气泡相自由能密度；h(η)是构造的一个插值函数，其表达式为h(η)＝η3(6η2-15η+10)；插值函数取值满足：在基体相中，即η＝0.0时，h(η)＝0.0；在气泡相中，即η＝1.0时，h(η)＝1.0；
S3.4：根据步骤S3.3获得的核燃料多晶多相体系体自由能密度，结合界面梯度能、多晶相互作用能、晶界与气泡相互作用能，得到核燃料多晶多相系统总自由能密度：
式中，表示多晶相互作用自由能密度，其表达式为：
其中，A、B、aGB、as是唯象参数；
表示界面梯度能，其表达式为：
其中，κv、κg、κη、是梯度项系数；
S3.5：考虑辐照条件下空位、气体原子产生的演化方程为：
空位浓度场演化方程：
气体原子浓度场演化方程：
气泡相序参量演化方程：
多晶序参量演化方程：
其中，Mv、Mg分别表示空位和气体原子的迁移系数；L是自由界面的迁移率；ξv、ξg、ξη分别为空位、气体原子和气泡相的热涨落项；Pv(r,t)、Pg(r,t)分别表示辐照条件下空位以及气体原子的产生率；
S3.6：将步骤S2获取的核燃料晶界气泡初始分布导入核燃料晶界气泡演化相场模型，作为核燃料晶界气泡演化相场模型初始化的依据；
S3.7：采用空间上的有限差分法以及时间上的显示欧拉法对演化方程进行求解。
5.根据权利要求1所述的计算核燃料中裂变气体释放及辐照肿胀行为的方法，其特征在于：步骤S4中，包括以下步骤：
S4.1：根据演化方程计算结果，将计算结果中相场变量信息以vtk文件形式存储；vtk文件中写入模拟区域大小信息；使用paraview软件导入vtk文件进行结果可视化，获得晶界气泡形貌演化图；
S4.2：根据演化方程计算结果，统计计算渗流参数方法如下：
晶界气泡密度、晶界气泡平均尺寸：晶界气泡密度定义为整个面晶界上气泡数量与模拟区域面晶界面积之比，单位为个/μm2；晶界气泡平均尺寸定义为：晶界上所有气泡面积之和与气泡数量之比，单位为μm2；采用Two-pass连通性分析的算法，对晶界气泡密度以及晶界气泡平均尺寸进行统计计算；Two-pass连通性分析的算法通过两遍扫描，将图像中存在的所有连通区域找出并标记；在第一次扫描时，从左至右、从上往下将所有区域内每个像素位置赋予一个label；扫描过程中同一个连通区域内的像素集合中可能被赋予一个或多个不同label；第二遍扫描时，需要将这些属于同一个连通区域但具有不同值的label合并，具体规则如下：
第一次扫描时：访问当前像素点，如果该label值等于1，则：
a.如果像素点领域中label值都为0，则赋予当前像素点一个新的label值：label+1；
b.如果像素点领域有label值>1的像素Neighbors，将领域中像素值最小的值赋予当前像素点；
c.记录Neighbors中各个label值之间的相等关系，即这些label值同属同一个连通区域；
第二次扫描时：访问当前像素点，如果label值大于1，找到与该label值同属相等关系的最小label值，把该值赋予当前像素点；
最后，统计区域内不同label值的个数，得到气泡数量；统计不同label值包含的区域面积，获得每个气泡的面积；
晶界气泡覆盖率：为获取气泡连通阈值，首先定义两个晶界气泡覆盖率的含义：晶界覆盖率：定义为晶界区域平面上所有气泡的投影面积除以晶界区域平面总面积；排放晶界覆盖率：定义为预设圆内接触圆周的气泡所占面积百分比；该预设圆假设为描绘晶界面的三重连接处，圆心位于晶界区域中心，半径为晶界区域边长的90％；采用以下算法实现参数计算：
a.扫描模拟区域中心晶界面，判断该区域上每个格点是否处于气泡中；
b.统计处于气泡中的格点个数，计算气泡总面积，气泡总面积比上晶界面面积即为晶界覆盖率；
c.根据不同label值得到的气泡分布，扫描整个区域，分别判断每个气泡是否与上述预设圆圆周接触；与预设圆圆周接触的所有气泡面积之和与气泡总面积之比即为排放晶界覆盖率；
气泡连通阈值：基于上述气泡晶界覆盖率的两个参数计算结果，画出排放晶界覆盖率随晶界覆盖率的变化曲线；排放晶界覆盖率随着晶界覆盖率增大而向1急速增长时，认为此时该晶界面气泡与晶界连通；取该急速增长过程中斜率最大点，该点处的晶界覆盖率即为气泡连通临界阈值；
S4.3：根据演化方程计算结果，燃料孔隙率计算方法如下：
孔隙率定义为：模拟中心晶界区域气泡所占体积分数；气泡体积统计方法为：
a.扫描整个模拟中心晶界区域，判断该区域中每个格点是否处于气泡中；
b.统计处于气泡中的格点个数，计算气泡总体积；气泡总体积比上模拟中心晶界区体积即为孔隙率。
6.根据权利要求1所述的计算核燃料中裂变气体释放及辐照肿胀行为的方法，其特征在于：步骤S5中，包括以下内容：
结合步骤S4中可视化结果，获得不同辐照强度下的燃料孔隙率随时间变化曲线、气泡平均尺寸变化曲线、气泡密度变化曲线；获得不同温度下的燃料孔隙率随时间变化曲线、气泡平均尺寸变化曲线、气泡密度变化曲线，得到核燃料辐照肿胀规律。
7.根据权利要求1所述的计算核燃料中裂变气体释放及辐照肿胀行为的方法，其特征在于：步骤S6中，包括以下步骤：
S6.1：基于核燃料芯块结构，晶粒结构采用四晶粒结构，晶界尺寸固定为10微米；
S6.2：基于轴对称假设，取燃料芯块r-z平面，在径向上设置500个晶界；晶界存在三种状态：闭合、开放、排气；每个晶界的初始状态都是闭合的；
S6.3：步骤S6.2中所述的500个晶界，每个晶界的气泡接触角在40°到80°之间随机分布。
8.根据权利要求1所述的计算核燃料中裂变气体释放及辐照肿胀行为的方法，其特征在于：步骤S7中，包括以下步骤：
S7.1：基于气体原子晶内扩散方程，并考虑重溶效应，建立核燃料裂变气体释放渗流模型；模型建立过程如下：进入晶面的气体流速用下式表示：
式中，N表示晶界气体浓度；f0表示晶内到达晶界的平均气体原子流速；δ表示晶界上的气泡重溶系数；b表示重溶常数；βg表示单位体积内核燃料内产生裂变气体原子的速率；
晶界网络的温度分布为：
其中，Tmax是燃料中心最高温度，Tmin是燃料外表面温度，r0是芯块外径；
单位体积内核燃料内产生裂变气体原子的速率：
βg＝-2.218×1018+3.854×1015(T)
气体原子扩散系数与温度的关系为：
其中，(a)适用于T<1381K，(b)适用于1381K<T<1650K，(c)适用于T>1650K；
S7.2：根据步骤S4所计算渗流参数，获取不同温度、不同裂变密度、不同气泡接触角下的气泡连通阈值，并按照输入变量顺序将计算结果排列并保存为txt文件；
S7.3：将步骤S6建立的核燃料晶界网络以及步骤S7.2获取的txt文件输入到步骤S7.1所述核燃料裂变气体释放渗流模型中，核燃料晶界气泡演化相场模型计算的渗流参数作为核燃料裂变气体释放渗流模型判断晶界是否开放的判据条件，建立了核燃料裂变气体释放行为相场-渗流耦合框架；所述判据条件为：当某个晶界的气体浓度达到连通阈值时，该晶界就被认为是连通状态；每个时间步，扫描晶界网络是否通过开放的晶界与燃料表面相连，若相连，则认为晶界内气体发生宏观排；在下一个时间步，排气的晶界被重置为关闭状态，晶界气体浓度重置为0。
CN202210332746.7A 2022-03-31 2022-03-31 一种计算核燃料中裂变气体释放及辐照肿胀行为的方法 Active CN114743607B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210332746.7A CN114743607B (zh) 2022-03-31 2022-03-31 一种计算核燃料中裂变气体释放及辐照肿胀行为的方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210332746.7A CN114743607B (zh) 2022-03-31 2022-03-31 一种计算核燃料中裂变气体释放及辐照肿胀行为的方法

Publications (2)

Publication Number Publication Date
CN114743607A CN114743607A (zh) 2022-07-12
CN114743607B true CN114743607B (zh) 2024-04-09

Family

ID=82279994

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210332746.7A Active CN114743607B (zh) 2022-03-31 2022-03-31 一种计算核燃料中裂变气体释放及辐照肿胀行为的方法

Country Status (1)

Country Link
CN (1) CN114743607B (zh)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115266795A (zh) * 2022-07-29 2022-11-01 中国核动力研究设计院 一种强放射性燃料元件裂变气体产物扩散行为表征方法
CN117935994B (zh) * 2024-03-20 2024-05-17 哈尔滨工业大学(深圳)(哈尔滨工业大学深圳科技创新研究院) 一种预测陶瓷核燃料辐照肿胀行为的方法

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2001351919A (ja) * 2000-06-05 2001-12-21 Nec Corp 配線故障解析方法
CN104156630A (zh) * 2014-09-05 2014-11-19 西南科技大学 三维核素扩散的计算方法
CN107766641A (zh) * 2017-10-16 2018-03-06 中国核动力研究设计院 一种计算uo2燃料裂变气体热释放率的方法
CN109448799A (zh) * 2018-09-03 2019-03-08 岭东核电有限公司 金属冷却快堆金属燃料多物理场模型耦合方法
CN111508573A (zh) * 2020-04-17 2020-08-07 西安交通大学 一种分析裂变气体导致铀硅化合物核燃料膨胀行为的方法及系统
WO2020237977A1 (zh) * 2019-05-27 2020-12-03 北京工业大学 一种多相复合材料力学行为的多尺度模拟方法

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2001351919A (ja) * 2000-06-05 2001-12-21 Nec Corp 配線故障解析方法
CN104156630A (zh) * 2014-09-05 2014-11-19 西南科技大学 三维核素扩散的计算方法
CN107766641A (zh) * 2017-10-16 2018-03-06 中国核动力研究设计院 一种计算uo2燃料裂变气体热释放率的方法
CN109448799A (zh) * 2018-09-03 2019-03-08 岭东核电有限公司 金属冷却快堆金属燃料多物理场模型耦合方法
WO2020237977A1 (zh) * 2019-05-27 2020-12-03 北京工业大学 一种多相复合材料力学行为的多尺度模拟方法
CN111508573A (zh) * 2020-04-17 2020-08-07 西安交通大学 一种分析裂变气体导致铀硅化合物核燃料膨胀行为的方法及系统

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
弥散型燃料的裂变气体行为研究;邢忠虎, 应诗浩;核动力工程;20001228(第06期);全文 *

Also Published As

Publication number Publication date
CN114743607A (zh) 2022-07-12

Similar Documents

Publication Publication Date Title
CN114743607B (zh) 一种计算核燃料中裂变气体释放及辐照肿胀行为的方法
CN105139444B (zh) 基于岩心二维颗粒图像的三维颗粒结构重建方法
Martínez-Freiría et al. Integrative phylogeographical and ecological analysis reveals multiple Pleistocene refugia for Mediterranean Daboia vipers in north-west Africa
CN107423462B (zh) 工件考虑三维粗糙表面形貌的疲劳寿命预测方法及系统
Wan et al. Optimal micro-siting of wind turbines by genetic algorithms based on improved wind and turbine models
Kruse et al. Treeline dynamics in Siberia under changing climates as inferred from an individual-based model for Larix
CN105354873B (zh) 用于多孔介质三维重构的模式密度函数模拟方法
Gravner et al. Modeling snow-crystal growth: A three-dimensional mesoscopic approach
CN110119590B (zh) 一种基于多源观测数据的水质模型粒子滤波同化方法
CN110263418B (zh) 一种体心立方合金微观偏析数值预测方法
CN101853485A (zh) 一种基于近邻传播聚类的非均匀点云简化处理方法
CN105896578B (zh) 一种用于风光储联合发电系统的随机生产模拟方法
CN210194532U (zh) 用于推移质输沙率研究的植被群落河道模型
Sintes et al. Modeling nonlinear seagrass clonal growth: assessing the efficiency of space occupation across the seagrass flora
CN109781044A (zh) 边坡失稳的综合逐渐逼近预警方法
Karafyllidis Design of a dedicated parallel processor for the prediction of forest fire spreading using cellular automata and genetic algorithms
CN110222368A (zh) 一种利用二维切片计算岩心三维孔隙度和渗透率的方法
CN109359431B (zh) 一种流动海水中材料表面点蚀的模拟方法
CN106485030A (zh) 一种用于sph算法的对称边界处理方法
CN116401888A (zh) 一种高镍正极材料晶粒生长过程的模拟方法及系统
WO2022242852A1 (en) Method for optimizing material properties of components of a battery, manufacturing a fiber network, an electrode and a battery
CN105956389A (zh) 环境因素数据的采集方法及金属大气腐蚀等级图绘制方法
CN115688439A (zh) 一种基于数字孪生的水库模型构建方法
CN105740513B (zh) 一种gh4169合金热变形动态再结晶模拟方法
Lohner et al. Improved adaptive refinement strategies for finite element aerodynamic computations

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