CN116644628B - 判定弥散型板燃料元件起泡与沸腾临界先后顺序的数值模拟方法 - Google Patents

判定弥散型板燃料元件起泡与沸腾临界先后顺序的数值模拟方法 Download PDF

Info

Publication number
CN116644628B
CN116644628B CN202310566449.3A CN202310566449A CN116644628B CN 116644628 B CN116644628 B CN 116644628B CN 202310566449 A CN202310566449 A CN 202310566449A CN 116644628 B CN116644628 B CN 116644628B
Authority
CN
China
Prior art keywords
fuel
boiling
calculating
critical
temperature
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
CN202310566449.3A
Other languages
English (en)
Other versions
CN116644628A (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 CN202310566449.3A priority Critical patent/CN116644628B/zh
Publication of CN116644628A publication Critical patent/CN116644628A/zh
Application granted granted Critical
Publication of CN116644628B publication Critical patent/CN116644628B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00Computer-aided design [CAD]
    • G06F30/20Design optimisation, verification or simulation
    • G06F30/23Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00Computer-aided design [CAD]
    • G06F30/20Design optimisation, verification or simulation
    • G06F30/28Design optimisation, verification or simulation using fluid dynamics, e.g. using Navier-Stokes equations or computational fluid dynamics [CFD]
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00Details relating to CAD techniques
    • G06F2111/10Numerical modelling
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06FELECTRIC DIGITAL DATA PROCESSING
    • G06F2113/00Details relating to the application field
    • G06F2113/08Fluids
    • 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
    • YGENERAL 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
    • Y02TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02EREDUCTION OF GREENHOUSE GAS [GHG] EMISSIONS, RELATED TO ENERGY GENERATION, TRANSMISSION OR DISTRIBUTION
    • Y02E30/00Energy generation of nuclear origin
    • Y02E30/30Nuclear fission reactors

Abstract

本发明公开了一种判定弥散型板燃料元件起泡与沸腾临界先后顺序的数值模拟方法。包括以下步骤:1、填写并读入输入卡;2、设定时间循环次数与计算网格结构;3、计算燃料元件径向及轴向功率分布;4、进行稳态热工参数计算,更新燃耗和芯块最高温度;5、更新包壳和燃料板厚度;6、更新流道尺寸,计算燃料组件温度与热流密度;7、计算流道出口含汽率,根据沸腾临界类型计算临界热流密度;8、进行热流密度收敛判断,计算沸腾临界阈值温度;9、进行燃料颗粒开裂判断,计算裂变气体释放量;10、进行起泡判断,计算起泡阈值温度;11、将步骤8与步骤10温度进行对比,判定弥散型板燃料元件起泡与沸腾临界发生先后顺序。

Description

判定弥散型板燃料元件起泡与沸腾临界先后顺序的数值模拟 方法
技术领域
本发明属于反应堆安全分析技术领域,具体涉及一种判定弥散型板燃料元件起泡与沸腾临界先后顺序的数值模拟方法。
背景技术
弥散型燃料元件(dispersion fuel element)是指采用弥散型燃料的元件。弥散型燃料是将细颗粒状的核燃料弥散在基体材料中形成的核燃料,由核裂变颗粒(如铀、钚等的化合物)分布在金属、陶瓷或石墨基体中构成,在结构上类似于颗粒增强复合材料,其基本思想是用基体材料把燃料颗粒相互隔离开。核燃料可以是金属铀或二氧化铀等;基体材料可以是铝、锆合金或石墨等。弥散型燃料的辐照损伤只限于弥散相本身及其邻近区域,基体受损较轻,所以燃耗可以很高。弥散型燃料的热学性能和机械性能与基体材料基本相似,通常它有较高的强度和塑性,良好的导热性和抗腐蚀性,能承受较大的热应力。
相较于传统棒状UO2陶瓷型燃料元件,板状燃料元件具有传热好,燃料芯体温度低等特性,可大幅度提高堆芯的功率体积比;同时,板状燃料元件加工简单,毋需复杂的表面加工和处理,在通道内高速流体的冲刷下,换热表面不易发生杂质的沉淀和污染,是一种较为先进的紧凑型堆芯。在板状燃料元件反应堆堆芯内,燃料元件之间形成若干平行的矩形流道,燃料板间隙一般为1~3mm,构成一种大宽高比的扁平矩形窄缝流道。冷却剂从入口到出口可能经历单相水、包括过冷沸腾和饱和沸腾的汽水两相流及单相蒸汽几种流型的转变,流动和传热工况极为复杂。
弥散型板状燃料元件矩形窄缝通道热流局部集中和起泡现象会使得热流分布和流道几何形状均发生畸变,是弥散型燃料元件热工水力设计和安全分析的关键瓶颈技术。起泡是指弥散型燃料元件内辐照产生的裂变气体在局部聚集,不断增大的气体压力导致燃料发生塑性变形,对矩形窄缝热工水力特性会产生不利的影响。对矩形通道内热流局部集中下的流动换热特性进行研究,需关注临界热流密度这一重要的安全准则。临界热流密度是核燃料元件表面发生传热恶化时的热通量,是冷却剂流动沸腾机理发生转变的结果。对热流局部集中下高温高压沸腾临界影响机理的研究,需要关注起泡发生和生长演变机理,从而获得沸腾临界与起泡阈值关系以及影响机制。
板状核燃料元件起泡发生过程与热工水力学密切相关,不同阶段传热模式会影响元件温度分布。因此,可通过探索沸腾临界发生时燃料元件温度以及起泡阈值温度之间的关系,建立沸腾临界和起泡阈值的安全准则,如果沸腾临界发生时元件温度低于起泡阈值温度,则沸腾临界限值准则可确保燃料元件安全性;如果高于起泡阈值温度,则需要梳理出这些运行工况并分析起泡后传热特性,为起泡后热工安全准则制定提供技术基础。获得弥散型燃料元件沸腾临界前壁温与起泡阈值温度的先后关系,对于理清弥散型燃料元件在各种工况范围下安全传热和可冷却性具有重要的意义。
发明内容
为了解决上述问题,本发明提供一种判定弥散型板燃料元件起泡与沸腾临界先后顺序的数值模拟方法,在给定燃料元件参数和热工水力参数之后,可以计算服役期间核反应堆板状燃料元件的主要参数,并对元件起泡与沸腾临界发生的先后顺序进行判定,考虑燃料颗粒、芯块和包壳的行为包括颗粒开裂、裂变气体释放、包壳膨胀等。
为了达到上述目的,本发明采用如下技术方案:
一种判定弥散型板燃料元件起泡与沸腾临界先后顺序的数值模拟方法,包括以下步骤:
步骤1:读入输入文件,填写输入卡,输入文件包括燃料元件的几何结构参数和材料参数,反应堆堆芯的热工水力参数以及初始计算工况;
步骤2:使用输入文件参数,根据输入时间步长的数目设定时间循环次数,根据输入轴向节点数目设定计算的网格结构;
步骤3:依据线平均功率和堆芯形状因子得到燃料元件长度、宽度和高度方向上不同位置处的功率分布;
步骤4:假设芯块最高温度位置,进行稳态热工参数计算,依次计算通道和堆芯热流密度分布、冷却剂流量,根据热平衡计算得到冷却剂温度、包壳温度、通道压降,更新燃耗分布并计算芯块温度,根据得到的芯块温度重新更新芯块最高温度位置调整;
步骤5:根据燃料元件的温度分布和燃耗计算出芯块热膨胀、燃料颗粒密实化和辐照肿胀,并计算出芯块的应变、燃料板应力,进行屈服判断并更新等效弹塑性增量和等效应力,从而更新包壳和燃料板厚度;
步骤6:根据步骤5更新后的包壳和燃料板厚度更新流道尺寸,判断是否为最后几何尺寸节点,是则进行下一个步骤计算,否则更新计算节点,重复步骤4,更新燃料组件温度以及热流密度;
步骤7:将步骤4-6获得的稳态热工参数作为输入条件进行瞬态计算,计算冷却剂通道出口含汽率,从而判断沸腾临界类型,进行沸腾临界热流密度计算,将含汽率与当前工况下发生环状流起始含汽率进行对比判断沸腾临界类型,若出现环状流则进行Dryout型沸腾临界计算得到沸腾临界热流密度,若未出现环状流则进行DNB型沸腾临界计算得到沸腾临界热流密度;
步骤8:计算步骤7得到的沸腾临界热流密度与步骤6得到的热流密度之间的误差,如未达到收敛标准,则更新热流密度,重复步骤7计算;如达到收敛标准,则判断发生沸腾临界,得到沸腾临界阈值温度Tchf
步骤9:将步骤4得到的温度和燃耗分布作为输入数据,根据边界条件求解燃料颗粒弹性力学基本方程,对燃料颗粒开裂进行判断,如燃料颗粒未开裂则进行反冲击出释放量计算;如燃料颗粒开裂则进行裂纹连通、气泡连通以及原子扩散释放方式释放量的计算,从而得到总的裂变气体释放量;
步骤10:根据步骤3-6得到的热工、机械参数和步骤9得到的总的裂变气体释放量,进行包壳屈服应力和弹性模量计算,计算出自由空间气体的释放量并计算燃料板缺陷处气体温度和燃料板的体积变化,并得到缺陷处裂变气体总压力与可容纳气体体积,从而计算出包壳最大应力和应变进行包壳起泡判断;如包壳所受应力小于屈服应力,则进行温度、压力和角度迭代重复计算;如包壳应力大于屈服应力,则计算对应起泡阈值温度Tblister
步骤11:将步骤8得到的沸腾临界阈值温度Tchf与步骤10得到的起泡阈值温度Tblister进行对比,从而判定弥散型板燃料元件起泡与沸腾临界发生的先后顺序。
步骤1中所述的输入文件中,反应堆堆芯的热工水力参数是分时间步长输入的热工水力参数;输入文件输入长度、宽度和高度方向上的节点数目,通过输入的节点数改变划分网格的疏密程度并适应实际的计算工况;热工水力参数包括冷却剂温度、冷却剂压力、冷却剂质量流速、堆芯平均线功率、反应堆功率形状因子、快中子注量率;分时间步长输入的热工水力参数储存在定义的全局变量数组中,在每次时间循环迭代时,读取新的时间步长数目和新的时间步长的长度和新的热工水力参数;计算核电厂参数有波动变化时,燃料元件性能随时间的变化情况。
步骤7所述的DNB型沸腾临界热流密度和Dryout型沸腾临界热流密度计算中,采用有限差分法对环状流三流体模型的方程组进行求解,对偏微分方程组时间项和空间项进行离散,将偏微分方程组转换成方便精确求解的常微分方程组;对空间项采用交错式差分格式离散,对时间项采用半隐式差分离散,除了质量守恒方程中对流项的速度,以及动量守恒方程中的压力梯度和界面相关速度项采用隐式格式以外,其他项均采用显式格式。
步骤7在计算冷却剂沸腾临界热流密度过程中,根据与环状流特征热流密度进行对比判断发生沸腾临界类型,再根据沸腾临界类型进行对应计算,计算DNB型沸腾临界或者Dryout型沸腾临界的热流密度与发生沸腾临界的温度,这样使得沸腾临界阈值温度的计算更加准确。
步骤9中对于弥散在板燃料内部的燃料颗粒状态,由颗粒内部产生裂变气体的压力和颗粒本身的开裂强度相对大小进行判断;使用位移法求解空间球对称应变,计算出燃料颗粒所受最大应力,与颗粒本身的机械开裂强度进行比较判断弥散燃料颗粒的开裂状态,这使得判断更加准确。
步骤9中计算总的裂变气体释放量时,从细观和宏观层面多角度分析,综合考虑燃料颗粒开裂前后裂变气体释放方式多样性;对于弥散型板状燃料元件而言,裂变气体释放在燃料颗粒开裂之前通过“反冲-击出”方式释放,燃料颗粒开裂后裂变气体通过气体原子扩散作用、裂纹与气泡连通、气泡与气泡连通三种方式进行扩散。
步骤10中对于包壳发生起泡的相关计算,通过对包壳屈服应力与芯块自由空间缺陷处裂变气体压力的计算,并且通过压力、温度迭代方式精确的判断包壳发生塑性形变失效状态,准确计算出起泡的高度和起泡面积份额,通过模拟计算得到起泡发生时的阈值温度。
步骤8所述的收敛标准为根据步骤7计算得到的沸腾临界热流密度与根据步骤6得到的热流密度之间的偏差小于设定的误差时,则认为步骤8计算达到收敛;否则更新步骤6热流密度值并重新迭代直到两热流密度间偏差小于设定误差。
本发明与现有技术相比,具有如下优点:
1、本发明创新性提出了一种判定弥散型板燃料元件起泡与沸腾临界先后顺序的数值模拟方法,通过计算可以同时获得弥散型板状燃料元件沸腾临界阈值温度以及燃料元件的起泡阈值温度,从而判断板状燃料元件发生起泡与发生沸腾临界的先后顺序,对于理清弥散型燃料元件在各种工况下换热机理和制定相关安全准则具有奠基作用。
2、步骤1中所述的输入文件中,热工水力参数是分时间步长输入的热工水力参数;输入文件可以输入长度、宽度和高度方向的节点数目,通过输入的节点数可以改变划分网格的疏密程度和适应实际的计算工况。分时间步长输入的热工水力参数储存在定义的全局变量数组之中,在每次时间循环迭代时,读取新的时间步长数目和新的时间步长的长度和新的热工水力参数。在每次计算时采用新的热工水力参数,这样和现代核电厂采用的燃料元件分析程序的输入文件是一致的,可以计算核电厂参数有波动变化时,燃料元件性能随时间的变化情况。
3、步骤7所述的对DNB型和Dryout型沸腾临界热流密度计算中,采用有限差分法对环状流三流体模型的方程组进行求解,对偏微分方程组时间项和空间项进行合理的离散,将偏微分方程组转换成可方便精确求解的常微分方程组。对空间项采用交错式差分格式离散,对时间项采用半隐式差分离散,除了质量守恒方程中对流项的速度,以及动量守恒方程中的压力梯度和界面相关速度项采用隐式格式以外,其他项均采用显式格式。采用这种离散方法相对于其他方法更加简便,精度完全符合工程要求,实现计算收敛速度和计算效率的提升。
4、步骤7在计算冷却剂沸腾临界热流密度过程中,根据与环状流特征热流密度进行对比判断发生沸腾临界类型,再根据沸腾临界类型判断进行对应计算,计算DNB型沸腾临界或Dryout型沸腾临界热流密度与发生沸腾临界的温度,这样使沸腾临界阈值温度的计算更加准确。
5、步骤9中对于弥散在板内部的燃料颗粒状态,由颗粒内部产生裂变气体的压力和颗粒本身的开裂强度相对大小进行判断。使用位移法求解空间球对称应变,计算出燃料颗粒所受最大应力,与颗粒本身的机械开裂强度进行比较判断弥散燃料颗粒的开裂状态,这使得判断更加准确。
6、步骤9对于弥散型板状燃料元件而言,裂变气体释放在燃料颗粒开裂之前通过“反冲-击出”方式释放,燃料颗粒开裂后裂变气体通过气体原子扩散作用、裂纹与气泡连通、气泡与气泡连通三种方式进行扩散。在计算裂变气体释放量时,从细观和宏观层面多角度分析,综合考虑了燃料颗粒开裂前后裂变气体释放方式多样性。
7、步骤10中对于包壳发生起泡的相关计算,通过对包壳屈服应力与芯块自由空间缺陷处裂变气体压力的计算,并且通过压力、温度迭代方式可以较为精确的判断包壳发生塑性形变失效状态,准确计算出起泡的高度和起泡面积份额,并且可以通过模拟计算得到起泡发生时的阈值温度,对比通常实验只能将燃料板在不断升高的温度下退火直至燃料板起泡来确定阈值温度的方法更直接简便。
附图说明
图1是本数值模拟方法的流程图。
图2是弥散型板状燃料元件几何模型图。
具体实施方式
下面结合附图和具体实施方式对本发明作进一步详细说明。
如图1所示,本发明涉及一种判定弥散型板燃料元件起泡与沸腾临界先后顺序的数值模拟方法,包括以下步骤:
步骤1:读入输入文件,输入文件包括燃料元件的几何结构参数,燃料元件的材料参数,反应堆堆芯的热工水力参数。在计算过程中,弥散型板状燃料元件几何模型如图2所示,燃料元件的几何结构参数包括燃料元件活性段长度(mm)、燃料元件活性段宽度(mm)、燃料元件活性段高度(mm)、包壳厚度(mm)、芯块厚度(mm)、冷却剂通道宽度(mm)。燃料元件的材料参数包括燃料富集度、燃料百分数、燃料颗粒直径(mm)。反应堆堆芯的热工水力参数是按照时间步长输入的参数,首先输入时间步长的数目和轴向节点数目,然后将每个时间步长对应的反应堆运行的时间长度(天)输入,对应每个时间步的热工水力参数包括:冷却剂温度(℃)、冷却剂压力(MPa)、冷却剂质量流速[kg/(s·m2)]、堆芯平均线功率(W/m)、反应堆功率形状因子、中子注量率。程序读取分时间步长输入的参数之后,将其储存在数组之中,并在后续的分时间步计算中,每次计算调用对应于新的时间步长的热工水力参数。
步骤2:使用输入文件参数,根据输入时间步长的数目设定时间步长循环次数,根据输入轴向节点数目设定计算的网格结构。
步骤3:依据线平均功率和堆芯形状因子得到燃料元件长度、宽度和高度方向上不同位置处的功率分布。
步骤4:假设芯块最高温度位置,进行稳态热工参数计算,依次计算通道和堆芯热流密度分布、冷却剂流量,根据热平衡计算得到冷却剂温度、包壳温度、通道压降,更新燃耗分布并计算芯块温度。根据得到的芯块温度重新更新芯块最高温度位置调整。
步骤5:进行机械力学应力应变计算,根据燃料元件的温度分布和燃耗计算出芯块热膨胀、燃料颗粒密实化和辐照肿胀,并计算出芯块的应变、燃料板应力,进行屈服判断并更新等效弹塑性增量和等效应力,从而更新包壳和燃料板厚度。
步骤6:根据步骤5更新后的包壳和燃料板厚度更新流道尺寸,判断是否为最后几何尺寸节点,是则进行下一个步骤计算,否则更新计算节点,重复步骤4更新燃料组件温度以及热流密度。
步骤7:将步骤4-6获得的稳态热工参数作为输入条件进行瞬态计算,计算冷却剂通道出口含汽率,从而判断沸腾临界类型,进行沸腾临界密度计算。将计算得到的含汽率与当前工况下发生环状流起始含汽率进行对比判断沸腾临界类型,若出现环状流则进行Dryout型沸腾临界计算得到沸腾临界热流密度,若未出现环状流则进行DNB型沸腾临界计算得到沸腾临界热流密度。
其中,对于DNB型沸腾临界计算沸腾临界热流密度qDNB,采用如下公式(1)计算:
式中:δ——微液层厚度/m;ρl,sat——饱和液相密度/kg·m-3;hfg——汽化潜热/J·kg-1;LB——汽块长度/m;UB——汽块移动速度/m·s-1
汽块长度LB为Helmholtz临界波长:
式中:σ——表面张力/N·m-1;ρl——液相的密度/kg·m-3;ρg——汽相的密度/kg·m-3;UB——汽块移动速度/m·s-1;Usb——微液层中液体流速/m·s-1
由于加热壁面附近的微液层非常薄且微液层中液体流速相对汽块移动速度非常小,可假设Usb等于零,汽块长度变为:
汽块移动速度UB通过轴向方向施加在汽块上的浮力FBa(下标a表示轴向方向)和拖拽力FD间的平衡确定:
FBa=FD 公式(4)
式中:DB——汽块当量直径/m;g——重力加速度/m s-2;θ——管道倾斜角度/°;aa——附加加速度场的轴向分量/m s-2;t——时间/s;CD——拖拽系数;UBL——汽块中心线所处的主流速度/m s-1。aa(t)是一个随时间变化的量,(UB-UBL)表示汽块相对主流的移动速度。
联立以上三个关系式,汽块速度可以表达成如下形式:
汽块当量直径DB可由Qi Sun提出的适合于自然循环工况下基于脱离加热表面的汽泡受力模型计算:
式中:C9——经验系数;σ——表面张力/N·m-1;ρf——饱和液相密度/kg·m-3;De——通道水力当量直径/m;G——质量流速/kg·m-2·s-1
微液层厚度δ采用如下公式计算:
求出δ后,代入公式(7)和公式(2)可计算得到UB和LB,再将新的UB和LB代入公式(9)计算得到新的δ,通过迭代计算直到以上几个参数值达到收敛。最后将迭代计算得到的δ,UB和LB代入公式(1)计算得到DNB型临界热流密度值。
对于Dryout型沸腾临界的计算,其含汽率计算采用广泛应用的环状流起始含汽率关系式是Wallis关系式:
式中:xann——环状流起始点含汽率;De——通道水力当量直径/m;g——重力加速度/kg·s-2;ρf——饱和液相密度/kg·m-3;ρg——汽相的密度/kg·m-3;G——质量流速/kg·m-2·s-1
假设环状流液膜完全覆盖流道壁面,则液滴和汽相与壁面的摩擦力可以忽略,但需要包括空泡份额、相间摩擦力、夹带率和沉积率、环状流起始点判定等本构方程的选取和建立。三流体方程共七个变量:压力,液膜、液滴和汽相的份额和速度,除了建立三流体的质量和动量守恒方程外,还需要结合三流体的份额关系式:
αfdg=1 公式(11)
式中,αf——液膜份额,αd——液滴份额,αg——汽相份额。
微液层厚度δ可用下式求得:
式中,w——矩形通道的宽边尺寸/m;s——矩形通道的窄边尺寸/m;αf——液相份额。
汽相速度Ug和液相速度Ul的计算如下:
其中,q——热流密度/J·m-2·s-1;L——矩形通道的长边尺寸/m;G——质量流速/kg·m-2·s-1;Cpl——液相定压比热容/J·kg-1·K-1;Tsat——饱和温度/K;Tin——通道入口温度/K;hfg——汽化潜热/J·kg-1;x——通道内含汽率。
根据公式(13)可以求出汽相热流密度qg和液相热流密度ql
发生Dryout型沸腾临界时,汽相热流密度与液相热流密度相等,且此时临界热流密度与两热流密度也相等,因此有:
qDryout=qg=ql 公式(15)
由于液滴的夹带和液膜的蒸发,液膜厚度沿着流道逐渐变薄,当液膜份额αf小于10-9时,认为沸腾临界发生,此时通道内含汽率x与由公式(10)求得的环状流起始含汽率xann相等,由公式(14)和公式(15)获得Dryout型沸腾临界热流密度。
步骤8:计算步骤7得到的沸腾临界热流密度与步骤6得到的热流密度之间的误差,如未达到收敛标准,则更新热流密度,重复步骤7计算;如达到收敛标准,则判断发生沸腾临界,得到沸腾临界阈值温度Tchf
步骤9:根据步骤4得到温度和燃耗分布输入数据,根据边界条件求解燃料颗粒弹性力学基本方程对燃料颗粒开裂进行判断,如燃料颗粒未开裂则进行“反冲击出”释放量计算;如燃料颗粒开裂则进行“裂纹连通”、“气泡连通”以及“原子扩散”释放方式释放量的计算,从而得到总的裂变气体释放量。首先是在燃料颗粒未开裂条件下,芯块表面附近来自裂变直接产生或者裂变碎片碰撞过程击出的气体原子才能释放。这种方式的释放份额可由经验公式获得:
其中,F——反冲-击出释放份额,R——芯块的表面体积比,plate和rod下标分别表示板状和棒状燃料元件;Bu——百分燃耗深度/%FIMA;C——经验系数。
对于燃料颗粒开裂状态的判断,可以通过求解燃料颗粒所受最大应力与燃料颗粒的开裂强度相对大小进行判断。对于弥散在燃料板中的燃料颗粒所受最大应力的计算,通常是对其基本力学方程和本构方程求解得到。本专利采用位移法求解空间中气孔-燃料相球壳模型和燃料相-裂变损伤层-金属基体球壳模型两个球对称应变,从而得到燃料颗粒在燃料基体内受到最大压力:
式中,γ——燃料物质的表面自由能,0.78J/m2;ρ——燃料颗粒中气孔率;ρfu——燃料颗粒在弥散燃料芯体中的体积分数;Rg——平均气孔半径,m;Rf——燃料颗粒半径,m;Pg——燃料颗粒中裂变气孔内压力,MPa;Pf——燃料颗粒约束压力,MPa;Pm——金属球壳约束约束压力,MPa;Fl——裂变碎片损伤层限制燃料球开裂的约束压力,MPa;Fm——金属基体球壳层限制燃料球开裂的约束压力,MPa。
燃料颗粒中气孔率ρ和平均气孔半径Rg主要与百分燃耗深度相关,表达式分别如下:
ρ=0.021+0.66Bu 公式(19)
Rg=0.1+2.35Bu 公式(20)
燃料颗粒中裂变气孔内压力Pg利用超高压实际气体的状态方程计算,具体如下:
其中,
式中:ng——裂变气体总量/mol;Vg——气体总体积/m3;R=8.31J/(mol·K)为气体普适常数;a=5.57×10-5m3、b=2.39×10-5m3分别实际气体状态方程参数;Df=10.96×106/(1+ρg)g/m3为燃料颗粒的密度;β为裂变气体的裂变产额;Mf=269g/mol为燃料相的摩尔质量。
裂变碎片损伤层和金属基体球壳层限制燃料球起裂的约束压力:
式中,σθ是球壳模型在θ方向上所受应力分量,下标2和3分别表示裂变损伤层球壳和金属基体球壳。
由应力分量可计算得到燃料颗粒约束压力Pf关系式:
由于辐照后燃料颗粒的力学性能相关数据较少并且较为离散,拟合经典实验数据得到燃料颗粒到燃料球相的弹性模量和开裂强度分别如下:
式中:T——摄氏温度,BU——百分燃耗。
在金属基弥散燃料的失效初期,裂纹主要存在于燃料颗粒中,未受裂变碎片损伤的金属基体会保持完整性。裂纹穿过裂变气体气泡时,与裂纹连通的裂变气体气泡中裂变气体释放进入裂纹,裂纹附近的裂变气体可能会通过扩散或气泡连通等作用进入裂纹。开裂后裂变气体释放量由裂纹连通释放、气泡连通释放以及原子扩散释放三种释放途径的释放量之和组成:
F=FC-B+FB-B+FD 公式(26)
F为裂变气体释放量,下标C-B,B-B和D分别表示裂纹连通、气泡连通和原子扩散三种释放方式。单个燃料颗粒完全贯穿开裂后,裂纹面连通的裂变气体气泡的数量NC-B为:
Scrack——裂纹面积;NBpS——裂纹扩展时单位面积裂纹面连通的气泡数量:
其中,df——燃料颗粒直径/m,db——裂变气泡平均直径/m,NB——裂变气体气泡数量为气泡总体积与单个气泡体积的比值:
燃料相中气体释考虑由Xe和Kr两种裂变气体释放组成,且由于在燃料颗粒边缘,裂变气体产物会通过反冲作用逸出燃料相进入金属基体,不会形成裂变气体气泡。根据Weber对于逸出份额与燃料颗粒尺寸以及裂变产物在燃料相中的反冲射程关系式,可以得到燃料相中气体摩尔总量nfuel
nfuel=nXe(1-FXe)+nKr(1-FKr) 公式(30)
式中,nXe——Xe气体释放量;nKr——Kr气体释放量;FXe——Xe气体逸出到裂变损伤层份额;FKr——Kr气体逸出到裂变损伤层份额。
由裂纹与气泡连通引起的裂变气体释放量为:
燃料颗粒开裂后由于气泡连通引起的裂变气体释放量为:
其中,Pp——燃料颗粒内裂变气体气泡的逾渗概率,可由Roncki给出的宏观参数与微观结构的关系式以及逾渗模型结合进行求解。
原子扩散方式释放主要考虑逸出至裂变损伤层的气体原子通过扩散方式进行释放,裂变损伤层中的裂变气体总量为nmetal
nmetal=nXeFXe+nKrFKr 公式(33)
扩散层中的裂变气体原子均可能扩散进入裂纹,根据裂变损伤层和扩散层的几何关系,可以计算从裂变损伤层迁移进入裂纹的气体量,即裂变气体原子扩散引起的释放量:
其中,λg为在一定退火温度T和时间t下,气体原子在金属基体损伤层中的扩散距离:
式中,Dg——裂变气体原子在未损伤金属中的扩散系数。
步骤10:根据步骤3-6得到的热工、机械参数和步骤9得到的总的裂变气体释放量,进行包壳屈服应力和弹性模量计算,计算出自由空间气体的释放量并计算燃料板缺陷处裂变气体温度和燃料板的体积变化,并得到缺陷处气体总压力与可容纳气体体积,从而计算出包壳最大应力和应变进行包壳起泡判断。如包壳所受应力小于屈服应力,则进行温度、压力和角度迭代重复计算;如包壳应力大于屈服应力,则计算对应起泡阈值温度Tblister
步骤11:将步骤8得到的沸腾临界阈值温度与步骤10得到的起泡阈值温度进行对比,从而判定弥散型板燃料元件起泡与沸腾临界发生的先后顺序。
以上内容是结合具体的优选实施方式对本发明所作的进一步详细说明,但不能认定本发明的具体实施方式仅限于此,对于本发明所属技术领域的普通技术人员来说,在不脱离本发明构思的前提下,还可以做出若干简单的推演或替换,都应当视为属于本发明由所提交的权利要求书确定专利保护范围。

Claims (8)

1.一种判定弥散型板燃料元件起泡与沸腾临界先后顺序的数值模拟方法,其特征在于:包括以下步骤:
步骤1:读入输入文件,填写输入卡,输入文件包括燃料元件的几何结构参数和材料参数,反应堆堆芯的热工水力参数以及初始计算工况;
步骤2:使用输入文件参数,根据输入时间步长的数目设定时间循环次数,根据输入轴向节点数目设定计算的网格结构;
步骤3:依据线平均功率和堆芯形状因子得到燃料元件长度、宽度和高度方向上不同位置处的功率分布;
步骤4:假设芯块最高温度位置,进行稳态热工参数计算,依次计算通道和堆芯热流密度分布、冷却剂流量,根据热平衡计算得到冷却剂温度、包壳温度、通道压降,更新燃耗分布并计算芯块温度,根据得到的芯块温度重新更新芯块最高温度位置调整;
步骤5:根据燃料元件的温度分布和燃耗计算出芯块热膨胀、燃料颗粒密实化和辐照肿胀,并计算出芯块的应变、燃料板应力,进行屈服判断并更新等效弹塑性增量和等效应力,从而更新包壳和燃料板厚度;
步骤6:根据步骤5更新后的包壳和燃料板厚度更新流道尺寸,判断是否为最后几何尺寸节点,是则进行下一个步骤计算,否则更新计算节点,重复步骤4,更新燃料组件温度以及热流密度;
步骤7:将步骤4-6获得的稳态热工参数作为输入条件进行瞬态计算,计算冷却剂通道出口含汽率,从而判断沸腾临界类型,进行沸腾临界热流密度计算,将含汽率与当前工况下发生环状流起始含汽率进行对比判断沸腾临界类型,若出现环状流则进行Dryout型沸腾临界计算得到沸腾临界热流密度,若未出现环状流则进行DNB型沸腾临界计算得到沸腾临界热流密度;
步骤8:计算步骤7得到的沸腾临界热流密度与步骤6得到的热流密度之间的误差,如未达到收敛标准,则更新热流密度,重复步骤7计算;如达到收敛标准,则判断发生沸腾临界,得到沸腾临界阈值温度Tchf
步骤9:将步骤4得到的温度和燃耗分布作为输入数据,根据边界条件求解燃料颗粒弹性力学基本方程,对燃料颗粒开裂进行判断,如燃料颗粒未开裂则进行反冲击出释放量计算;如燃料颗粒开裂则进行裂纹连通、气泡连通以及原子扩散释放方式释放量的计算,从而得到总的裂变气体释放量;
步骤10:根据步骤3-6得到的热工、机械参数和步骤9得到的总的裂变气体释放量,进行包壳屈服应力和弹性模量计算,计算出自由空间气体的释放量并计算燃料板缺陷处气体温度和燃料板的体积变化,并得到缺陷处裂变气体总压力与可容纳气体体积,从而计算出包壳最大应力和应变进行包壳起泡判断;如包壳所受应力小于屈服应力,则进行温度、压力和角度迭代重复计算;如包壳应力大于屈服应力,则计算对应起泡阈值温度Tblister
步骤11:将步骤8得到的沸腾临界阈值温度Tchf与步骤10得到的起泡阈值温度Tblister进行对比,从而判定弥散型板燃料元件起泡与沸腾临界发生的先后顺序。
2.根据权利要求1所述的一种判定弥散型板燃料元件起泡与沸腾临界先后顺序的数值模拟方法,其特征在于:
步骤1中所述的输入文件中,反应堆堆芯的热工水力参数是分时间步长输入的热工水力参数;输入文件输入长度、宽度和高度方向上的节点数目,通过输入的节点数改变划分网格的疏密程度并适应实际的计算工况;热工水力参数包括冷却剂温度、冷却剂压力、冷却剂质量流速、堆芯平均线功率、反应堆功率形状因子、快中子注量率;分时间步长输入的热工水力参数储存在定义的全局变量数组中,在每次时间循环迭代时,读取新的时间步长数目和新的时间步长的长度和新的热工水力参数;计算核电厂参数有波动变化时,燃料元件性能随时间的变化情况。
3.根据权利要求1所述的一种判定弥散型板燃料元件起泡与沸腾临界先后顺序的数值模拟方法,其特征在于:
步骤7所述的DNB型沸腾临界热流密度和Dryout型沸腾临界热流密度计算中,采用有限差分法对环状流三流体模型的方程组进行求解,对偏微分方程组时间项和空间项进行离散,将偏微分方程组转换成方便精确求解的常微分方程组;对空间项采用交错式差分格式离散,对时间项采用半隐式差分离散,除了质量守恒方程中对流项的速度,以及动量守恒方程中的压力梯度和界面相关速度项采用隐式格式以外,其他项均采用显式格式。
4.根据权利要求1所述的一种判定弥散型板燃料元件起泡与沸腾临界先后顺序的数值模拟方法,其特征在于:
步骤7在计算冷却剂沸腾临界热流密度过程中,根据与环状流特征热流密度进行对比判断发生沸腾临界类型,再根据沸腾临界类型进行对应计算,计算DNB型沸腾临界或者Dryout型沸腾临界的热流密度与发生沸腾临界的温度,这样使得沸腾临界阈值温度的计算更加准确。
5.根据权利要求1所述的一种判定弥散型板燃料元件起泡与沸腾临界先后顺序的数值模拟方法,其特征在于:
步骤9中对于弥散在板燃料内部的燃料颗粒状态,由颗粒内部产生裂变气体的压力和颗粒本身的开裂强度相对大小进行判断;使用位移法求解空间球对称应变,计算出燃料颗粒所受最大应力,与颗粒本身的机械开裂强度进行比较判断弥散燃料颗粒的开裂状态,这使得判断更加准确。
6.根据权利要求1所述的一种判定弥散型板燃料元件起泡与沸腾临界先后顺序的数值模拟方法,其特征在于:
步骤9中计算总的裂变气体释放量时,从细观和宏观层面多角度分析,综合考虑燃料颗粒开裂前后裂变气体释放方式多样性;对于弥散型板状燃料元件而言,裂变气体释放在燃料颗粒开裂之前通过“反冲-击出”方式释放,燃料颗粒开裂后裂变气体通过气体原子扩散作用、裂纹与气泡连通、气泡与气泡连通三种方式进行扩散。
7.根据权利要求1所述的一种判定弥散型板燃料元件起泡与沸腾临界先后顺序的数值模拟方法,其特征在于:
步骤10中对于包壳发生起泡的相关计算,通过对包壳屈服应力与芯块自由空间缺陷处裂变气体压力的计算,并且通过压力、温度迭代方式精确的判断包壳发生塑性形变失效状态,准确计算出起泡的高度和起泡面积份额,通过模拟计算得到起泡发生时的阈值温度。
8.根据权利要求1所述的一种判定弥散型板燃料元件起泡与沸腾临界先后顺序的数值模拟方法,其特征在于:步骤8所述的收敛标准为根据步骤7计算得到的沸腾临界热流密度与根据步骤6得到的热流密度之间的偏差小于设定的误差时,则认为步骤8计算达到收敛;否则更新步骤6热流密度值并重新迭代直到两热流密度间偏差小于设定误差。
CN202310566449.3A 2023-05-19 2023-05-19 判定弥散型板燃料元件起泡与沸腾临界先后顺序的数值模拟方法 Active CN116644628B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202310566449.3A CN116644628B (zh) 2023-05-19 2023-05-19 判定弥散型板燃料元件起泡与沸腾临界先后顺序的数值模拟方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202310566449.3A CN116644628B (zh) 2023-05-19 2023-05-19 判定弥散型板燃料元件起泡与沸腾临界先后顺序的数值模拟方法

Publications (2)

Publication Number Publication Date
CN116644628A CN116644628A (zh) 2023-08-25
CN116644628B true CN116644628B (zh) 2023-11-07

Family

ID=87622343

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202310566449.3A Active CN116644628B (zh) 2023-05-19 2023-05-19 判定弥散型板燃料元件起泡与沸腾临界先后顺序的数值模拟方法

Country Status (1)

Country Link
CN (1) CN116644628B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116825408B (zh) * 2023-08-31 2023-11-03 哈尔滨工程大学 一种模拟板状燃料包壳起泡的可视化实验装置及实验方法

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2008261693A (ja) * 2007-04-11 2008-10-30 Toshiba Corp 熱的限界出力相関式作成方法および燃料集合体設計方法
CN106055850A (zh) * 2016-07-18 2016-10-26 西安交通大学 一种获得偏离泡核沸腾型临界热流密度的方法
CN115586209A (zh) * 2022-10-10 2023-01-10 西安交通大学 运动条件下热流局部集中窄矩形通道沸腾临界试验段及试验方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP7277176B2 (ja) * 2019-02-28 2023-05-18 キヤノン株式会社 ウルトラファインバブル生成方法、およびウルトラファインバブル生成装置

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2008261693A (ja) * 2007-04-11 2008-10-30 Toshiba Corp 熱的限界出力相関式作成方法および燃料集合体設計方法
CN106055850A (zh) * 2016-07-18 2016-10-26 西安交通大学 一种获得偏离泡核沸腾型临界热流密度的方法
CN115586209A (zh) * 2022-10-10 2023-01-10 西安交通大学 运动条件下热流局部集中窄矩形通道沸腾临界试验段及试验方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
板状燃料元件流道堵塞事故预防与探测技术研究;丁丽;骆贝贝;花晓;宁波;乔雅馨;;核技术(第04期);1-7 *

Also Published As

Publication number Publication date
CN116644628A (zh) 2023-08-25

Similar Documents

Publication Publication Date Title
Chang et al. A study of a nuclear hydrogen production demonstration plant
CN116644628B (zh) 判定弥散型板燃料元件起泡与沸腾临界先后顺序的数值模拟方法
CN108140432A (zh) 具有多重有效密度燃料的燃料元件
CN114444413A (zh) 一种板状燃料堆芯亚通道级三维热工水力分析方法
Espinosa-Paredes et al. Fuel rod model based on Non-Fourier heat conduction equation
Song et al. High burnup fuel technology in Korea
Sun et al. Various bypass flow paths and bypass flow ratios in HTR-PM
Cai et al. Research on thermal behaviors of several accident tolerant fuels based on 5× 5 bundle subchannel model
Duan et al. Multi-physics coupling analysis on neutronics, thermal hydraulic and mechanics characteristics of a nuclear thermal propulsion reactor
Ade et al. Particle Packing Characteristics in Dense TRISO/SiC Fuel Elements and the Impact of Neutronics and Thermomechanics
Cong et al. Development and evaluation of fuel performance analysis code FuSPAC
Gui et al. Experimental investigation on quenching behavior during reflooding in a 3× 3 dual-cooled annular rod bundle
Gutowska Study on depressurized loss of coolant accident and its mitigation method framework at very high temperature gas cooled reactor
CN113065241A (zh) 一种预测超临界二氧化碳冷却堆燃料元件主要参数的方法
Oh et al. Implications of Air Ingress Induced by Density-Difference Driven Stratified Flow
Jun et al. The benchmark calculations of the GAMMA+ code with the HTR-10 safety demonstration experiments
Katscher Coated particle fuel element for pressurized water reactors
Zhang et al. Numerical study on frictional resistance characteristics in a rod bundle channel under inclined conditions
Wu et al. Simulation on the Multidimensional Behavior of TRISO Particle Fuel With UN/UO2 Kernel
Deng et al. Subchannel thermal hydraulic analysis of 5× 5 rod bundle with CRUD layer
Likhanskii et al. Adaptability of RTOP fuel performance code for thermo-mechanical 2D and 3D simulations of fuel rods with various burnable neutron absorbers
Oh et al. Computational Fluid Dynamics analyses on very high temperature reactor air ingress
van Wijk Computational Modeling of the Flow Field in a Molten Salt Reactor Core
Zhang et al. Three-dimensional modeling and thermal hydraulic analysis of high temperature gas-cooled reactor core
Zhao et al. Single Phase CFD Method Validation Based on PIV Experimental Data

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