CN113192567B - 一种核反应堆板状燃料熔化流固耦合无网格分析方法 - Google Patents

一种核反应堆板状燃料熔化流固耦合无网格分析方法 Download PDF

Info

Publication number
CN113192567B
CN113192567B CN202110484617.5A CN202110484617A CN113192567B CN 113192567 B CN113192567 B CN 113192567B CN 202110484617 A CN202110484617 A CN 202110484617A CN 113192567 B CN113192567 B CN 113192567B
Authority
CN
China
Prior art keywords
particle
particles
fuel
solid
calculation
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
CN202110484617.5A
Other languages
English (en)
Other versions
CN113192567A (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 CN202110484617.5A priority Critical patent/CN113192567B/zh
Publication of CN113192567A publication Critical patent/CN113192567A/zh
Application granted granted Critical
Publication of CN113192567B publication Critical patent/CN113192567B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

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
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/25—Design optimisation, verification or simulation using particle-based methods
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/28—Design optimisation, verification or simulation using fluid dynamics, e.g. using Navier-Stokes equations or computational fluid dynamics [CFD]
    • 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
    • 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/90—Programming languages; Computing architectures; Database systems; Data warehousing
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2113/00—Details relating to the application field
    • G06F2113/08—Fluids
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces
    • 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)
  • Physics & Mathematics (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • General Physics & Mathematics (AREA)
  • Chemical & Material Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • General Engineering & Computer Science (AREA)
  • Crystallography & Structural Chemistry (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • Mathematical Analysis (AREA)
  • Pure & Applied Mathematics (AREA)
  • Algebra (AREA)
  • Fluid Mechanics (AREA)
  • Software Systems (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Physics (AREA)
  • Databases & Information Systems (AREA)
  • Health & Medical Sciences (AREA)
  • General Health & Medical Sciences (AREA)
  • Spectroscopy & Molecular Physics (AREA)
  • Analytical Chemistry (AREA)
  • Chemical Kinetics & Catalysis (AREA)
  • Monitoring And Testing Of Nuclear Reactors (AREA)

Abstract

一种核反应堆板状燃料熔化流固耦合无网格分析方法，具体步骤如下：1、在笛卡尔坐标系下对板状核燃料进行初始建模；2，采用链表法对邻居粒子进行检索，将邻居粒子信息保存在数组内；3，显式计算动量守恒方程中的重力项、粘性项和表面张力项，得到粒子的估算速度和位置；4，计算粒子间传热，使用焓值相变模型判断粒子的相态；5，检索是否存在固体粒子之间的碰撞，使用DEM模型进行计算；6，隐式计算压力梯度项，修正粒子的速度和位置；7，输出需要的结果，推进时间步长直到模拟结束。本方法考虑了核反应堆板状燃料熔化流固耦合的所有现象，粒子法能够对熔融物自由液面进行精确捕捉，固液耦合模型能够分析熔融物和燃料微粒之间的作用。

Description

一种核反应堆板状燃料熔化流固耦合无网格分析方法
技术领域
本发明涉及核电厂严重事故板状核燃料发生熔化时熔融物的迁徙行为的研究技术领域，具体涉及一种核反应堆板状燃料熔化流固耦合无网格分析方法。
背景技术
相较于一般核电站使用的棒状燃料元件，使用弥散颗粒压制的板状核燃料构成的核组件具有更加紧凑的结构，燃料芯体有更大的换热面积所以温度低，较高的换热效率和较深的燃耗等特点，因此被广泛的应用在一体化的反应堆、核动力反应堆和实验用研究堆中。在事故工况下，反应堆中出现如燃料元件在燃耗过深时发生辐照肿胀导致冷却剂流道变窄，而此时堆内材料碎片或异物随冷却剂的循环流入到堆芯内，就可能会引发冷却剂流道堵塞事故，造成堆芯燃料板无法被充分冷却，热量不能及时通过冷却剂带走，燃料板温度升高，进一步会导致局部冷却剂蒸干，完全裸露后的板燃料温度会急剧上升，极大的威胁了燃料包壳的完整性，甚至可能造成反射性外泄，引发严重的核事故后果。历史上，在1965年美国橡树岭研究堆发生过一起因橡胶垫片随冷却剂流入堆芯而引起的堵流事故，并最终造成一个板状燃料组件中的三块燃料板发生局部熔化。此外，由于目前为止国内外的研究主要集中在棒状核燃料的相关问题，因此对于棒状燃料元件在严重事故下的行为和机理的认识比较深入，结合已有的一些严重事故燃料和堆芯行为研究结果的机理性分析程序(如：MELCOR和SCDAP/RELAP5等)，得到了一些针对棒状燃料元件堆芯的较为准确的模拟结果，但国内外对板状核燃料元件的研究较少，已建立的模型存在相对简单、粗糙的问题，且现有的严重事故程序无法满足使用板状核燃料的反应堆严重事故分析的需求。综上所述，有必要针对板状燃料元件开展严重事故机理及关键技术的研究，为板状燃料元件核反应堆制定完善的严重事故缓解策略提供技术支持，进一步能够有针对性的制定一些防止严重事故发生的措施，减少严重事故发生的可能性和概率。
发明内容
为揭示板状核燃料发生熔化后熔融物的迁徙行为，从而能够给板状燃料元件核反应堆制定完善的严重事故缓解策略提供技术支持，制定一些防止严重事故发生的措施，减少严重事故发生的可能性和概率，本发明在粒子法的基础上，添加能够对固体之间作用进行分析的非连续介质模型的分析方法离散单元法并充分考虑固体和液体之间的耦合作用提出一种核反应堆板状燃料熔化流固耦合无网格分析方法，该方法能够对板状核燃料熔融物的迁徙运动进行研究，获得熔融物产生位置信息、开始迁徙后的运动信息等，板状燃料元件核反应堆制定完善的严重事故缓解策略提供技术支持。
为了实现上述目标，本发明采取了以下的技术方案予以实施：
一种核反应堆板状燃料熔化流固耦合无网格分析方法，步骤如下：
步骤1：对核反应堆板状燃料熔化流固耦合问题进行粒子建模；用不同类型的粒子代表核反应堆板状燃料的各个组成部分，0号粒子代表锆合金材质的燃料包壳和燃料基体，1号粒子代表弥散在燃料基体中的二氧化铀燃料微粒；对每一个粒子进行编号，所有的粒子均具有对应的物性参数，物性参数包括质量、密度、比热、熔沸点、温度、焓和初速度；粒子在笛卡尔坐标系内，依据核反应堆板状燃料的具体参数均匀布置，能够模拟核反应堆板状燃料的初始状态；
步骤2：使用链表法检索邻居粒子，具体方法是：将整个计算域均匀划分成正方形表格，正方形表格的边长与粒子最大作用半径相等；对于任意一个粒子而言，检索邻居粒子时只需检索包含粒子在内的表格以及包围此表格的共计9个表格，三维计算时，需要检索的就是由正方形表格围成的正方体，需要检索27个三维正方体；
步骤3：对时间步长进行稳定性分析，采用Courant-Friedrichs-Lewy条件简称CFL条件判断；CFL条件判断是差分方法保持计算稳定的必要条件，对于粒子法也适用，在移动粒子半隐式方法中采用下式进行判断
Figure BDA0003049813050000031
式中：
Δt——能够保持计算稳定的最小时间步长；
C——库朗常数；
l0——粒子初始布置时相邻两粒子间的距离；
umax——各时间步长中粒子速度最大值；
计算输入的初始条件中包含预设的时间步长大小，将输入的时间步长和使用CFL条件计算得到的能够保持计算稳定的最小时间步长进行比较，使用较大的时间步长作为当前时间步长；
步骤4：使用核函数表征粒子之间的相互作用程度，具体为在对动量守恒方程中的各项进行显式计算时，先使用核函数进行离散然后积分得到动量守恒方程中各项的具体值，使用的核函数如下所示：
Figure BDA0003049813050000041
式中
w(rij)——核函数值；
rij——粒子i与粒子j之间的距离；
re——粒子作用域半径；
步骤5：对核反应堆板状燃料熔化后的熔融物研究对象采用如下控制方程：
Figure BDA0003049813050000042
Figure BDA0003049813050000043
式中
ρ——粒子对应物质的密度kg/m3；
t——时间s；
P——压力Pa；
μ——动力粘度系数N·s/m2；
Figure BDA0003049813050000044
——速度矢量m/s；
Figure BDA0003049813050000045
——表面张力矢量N/kg；
Figure BDA0003049813050000046
——重力加速度m/s2；
公式(3)为连续性方程，针对的对象为不可压缩流体，公式(4)为动量守恒方程，从中看出分析对象只受到重力这一外力作用，对于任意一个粒子，通过核函数先将动量守恒方程即公式(4)中的重力项、粘性项和表面张力项离散，然后对计算域内全部粒子进行积分，显式得到粒子的估算速度和位置，其中粘性项需要使用采用如下多相粘度模型：
Figure BDA0003049813050000051
式中
μ——运动粘度系数m2/s
μij——粒子间动力粘度系数m2/s；
Figure BDA0003049813050000052
——速度m/s；
d——直径m；
re——粒子作用半径m；
nG——使用高斯核函数计算得到的粒子数密度；
Figure BDA0003049813050000053
——粒子j速度m/s；
Figure BDA0003049813050000054
——粒子i速度m/s；
G——高斯核函数；
rj——粒子j的位置矢量；
ri——粒子i的位置矢量；
步骤6：使用如下的焓值相变模型计算粒子的相态：
Figure BDA0003049813050000055
式中
T——需要计算的粒子温度K；
Ts——需要计算的粒子对应的熔点K；
h——当前粒子的焓值J/kg；
hs0——粒子开始熔化的焓值J/kg；
hs1——粒子结束熔化的焓值J/kg；
cp——粒子等压比热容J/(kg·K)；
通过步骤6的计算，能够模拟核反应堆板状燃料在最高温度超过燃料基体锆合金的熔点时的相变过程；计算得到每个粒子在不同时刻下的种类、焓值和温度，即得到锆合金基体和铀燃料微粒的相态、焓值和温度随时间的变化过程；
步骤7：得到粒子的相态后，通过非连续介质模型的分析方法离散单元法DEM对固体之间的碰撞作用进行如下的计算：
非连续介质模型的分析方法离散单元法DEM的计算对象为刚性固体，当计算对象中存在大固体时，需要先解决将离散颗粒聚集成刚性物体统一运动的问题；采用如下公式初始化大固体运动速度；
Figure BDA0003049813050000061
式中
Figure BDA0003049813050000062
——标记为ii的刚体的速度矢量m/s；
Figure BDA0003049813050000063
——粒子i当前的速度矢量m/s；
Figure BDA0003049813050000064
——标记为ii的刚体的角速度矢量rad/s；
mi——粒子i的质量kg；
Figure BDA0003049813050000065
——粒子i的位置矢量；
Figure BDA0003049813050000071
——标记为ii的刚体的位置矢量；
Iii——标记为ii的刚体的转动惯量；
对于任意固体单元mii，采用如下公式对单元受力作用后产生的运动进行分析；
Figure BDA0003049813050000072
Figure BDA0003049813050000073
式中
Figure BDA0003049813050000074
——标记为ii的刚体的加速度矢量；
Figure BDA0003049813050000075
——标记为ii的刚体的角加速度矢量；
Figure BDA0003049813050000076
——标记为ii的刚体的速度矢量；
Figure BDA0003049813050000077
——标记为ii的刚体的角速度矢量；
Figure BDA0003049813050000078
——标记为ii的刚体所受合力矢量；
Figure BDA0003049813050000079
——标记为ii的刚体所受合力矩矢量；
公式(8)和公式(9)中固体单元所受的合力与和力矩通过下式计算
Figure BDA00030498130500000710
Figure BDA00030498130500000711
式中
FX,jj→ii——标记为jj的刚体对标记为ii的刚体施加的法向合力；
en——法向弹性力；
dn——法向阻尼力；
kn——法向弹性系数；
cn——法向阻尼系数；
FY,jj→ii——标记为jj的刚体对标记为ii的刚体施加的切向合力；
es——切向弹性力；
ds——切向阻尼力；
ks——切向弹性系数；
cs——切向阻尼系数；
这样，固体单元mii就移动到了一个新的位置，并产生新的接触力和接触力矩，再次计算固体单元mii所受的合力
Figure BDA0003049813050000081
和合力矩
Figure BDA0003049813050000082
重复计算流程，即得到每个固体单元及整个散粒体的运动形态，此时固体粒子的速度和位置均为估算值；
步骤7：通过核反应堆板状燃料熔化后的熔融物的不可压缩特性和动量守恒方程中的压力项联立得到如下压力泊松方程，隐式求解如下的压力泊松方程
Figure BDA0003049813050000083
式中
ρ——粒子对应物质密度kg/m3；
P——压力Pa；
n0——初始粒子数密度常数；
n*——临时粒子数密度；
nk——当前时层的粒子数密度；
nk-1——上一时层的粒子数密度；
nk-1——下一时层的粒子数密度；
β——压力源项修正系数；
γ——压力源项修正系数；
通过步骤8的计算得到动量守恒方程中的压力项，使用该压力项所得值对步骤5和步骤7中得到的粒子估算速度和位置进行修正，得到当前时间步长全部粒子的真实速度和位置，这样就完成了一个时间步长的核反应堆板状燃料熔化流固耦合的模拟；
步骤9：为了准确分析模拟对象中固体和液体之间的耦合作用，还需要对步骤5中的粘性项计算和步骤8中的压力项计算进行处理，这是因为，固体和液体之间的作用通过下式表现
Figure BDA0003049813050000091
式中
Figure BDA0003049813050000092
——流体对固体的作用合力；
Figure BDA0003049813050000093
——流体对固体的压力作用；
Figure BDA0003049813050000094
——流体对固体的粘性作用；
在步骤5的粘性项计算中，通过多相粘度模型中使用的调和粘性系数来表现固体的物性，在计算时，这个粘性项只发生在固体和液体接触时，对固体内部不进行粘性项的计算；在步骤8的压力项计算中，流体粒子只和接触的固体粒子施加压力作用，后续压力作用将在非连续介质模型的分析方法离散单元法DEM的计算中扩展至全部固体粒子；
综上，通过步骤1对核反应堆板状燃料熔化流固耦合问题中，锆合金材质的燃料包壳和燃料基体和二氧化铀燃料微粒的位置、速度和初始物性参数进行设定；通过步骤6模拟板燃料的锆合金基体在温度超过其熔点后发生相变形成熔融物，计算得到每个粒子在不同时刻下的种类、焓值和温度，得到熔融物的相态、焓值和温度随时间的变化过程；通过步骤5至步骤9模拟核反应堆板状燃料熔化流固耦合过程，模拟熔融物和未熔化的燃料微粒之间的相互作用，模拟燃料微粒之间的相互作用过程，计算得到熔融物和燃料微粒粒子的速度、位置和压力，即得到熔融物迁徙过程中的运动和压力变化过程；综合以上步骤，分析了核反应堆板状燃料熔化流固耦合问题，得到熔化流固耦合过程中熔融物的位置、速度、压力、相态、温度、焓值随时间的变化，通过以上步骤对熔融物的迁徙行为展开机理性分析。
本发明方法能够给板状燃料元件核反应堆制定完善的严重事故缓解策略提供技术支持，制定一些防止严重事故发生的措施，减少严重事故发生的可能性和概率。
和现有技术相比，本发明方法具备如下优点：
本方法考虑了核反应堆板状燃料熔化流固耦合问题的所有现象，能够对板状核燃料熔融物的迁徙行为进行准确模拟；相比于传统的网格法，粒子法能够精确捕捉自由液面、方便对计算对象进行建模及根据焓值精确的计算相变；添加了固液耦合作用模型，能够分析熔融物和固体颗粒之间的作用；使用CPU+GPU异构并行对程序进行加速，有利于进行大规模粒子计算。
附图说明
图1是本发明对核反应堆板状燃料熔化流固耦合问题进行分析的流程图。
具体实施方式
下面结合附图和具体实施方式对本发明做进一步详细说明。
如图1所示，本发明一种核反应堆板状燃料熔化流固耦合无网格分析方法，步骤如下：
步骤1：对核反应堆板状燃料熔化流固耦合问题进行粒子建模；与网格方法不同，粒子法使用不同类型的粒子代表核反应堆板状燃料的各个组成部分，本方法中0号粒子代表锆合金材质的燃料包壳和燃料基体，1号粒子代表弥散在燃料基体中的二氧化铀燃料微粒；对每一个粒子进行编号，所有的粒子均具有对应的物性参数，物性参数包括质量、密度、比热、熔沸点、温度、焓和初速度；粒子在笛卡尔坐标系内，依据核反应堆板状燃料的具体参数均匀布置，能够模拟核反应堆板状燃料的初始状态；
步骤2：使用链表法检索邻居粒子，具体方法是：将整个计算域均匀划分成正方形表格，正方形表格的边长与粒子最大作用半径相等；对于任意一个粒子而言，检索邻居粒子时只需检索包含粒子在内的表格以及包围此表格的共计9个表格，三维计算时，需要检索的就是由正方形表格围成的正方体，需要检索27个三维正方体；
步骤3：对时间步长进行稳定性分析，采用Courant-Friedrichs-Lewy条件简称CFL条件判断；CFL条件判断是差分方法保持计算稳定的必要条件，对于粒子法也适用，在移动粒子半隐式方法中采用下式进行判断
Figure BDA0003049813050000111
式中
Δt——CFL条件认为能够保持计算稳定的最小时间步长；
C——库朗常数；
l0——粒子初始布置时相邻两粒子间的距离；
umax——各时间步长中粒子速度最大值；
计算输入的初始条件中包含预设的时间步长大小，将输入的时间步长和使用CFL条件计算得到的能够保持计算稳定的最小时间步长进行比较，使用较大的那个时间步长作为当前时间步长；
步骤4：使用核函数表征粒子之间的相互作用程度，具体表现为在对动量守恒方程中的各项进行显式计算时，先使用核函数进行离散然后积分得到动量守恒方程中各项的具体值，本方法使用的核函数如下所示：；
Figure BDA0003049813050000121
式中
w(rij)——核函数值；
rij——粒子i与粒子j之间的距离；
re——粒子作用域半径；
步骤5：对核反应堆板状燃料熔化后的熔融物采用如下控制方程：
Figure BDA0003049813050000122
Figure BDA0003049813050000123
式中
ρ——粒子对应物质的密度kg/m3；
t——时间s；
P——压力Pa；
μ——动力粘度系数N·s/m2；
Figure BDA0003049813050000131
——速度矢量m/s；
Figure BDA0003049813050000132
——表面张力矢量N/kg；
Figure BDA0003049813050000133
——重力加速度m/s2；
公式(3)为连续性方程，针对的对象为不可压缩流体，在粒子法中通过粒子数密度来代替真实的物质密度进行计算，通过保持粒子数密度一致来体现流体的不可压缩特性，公式(4)为动量守恒方程，从中可以看出分析对象只受到重力这一外力作用，对于任意一个粒子，通过核函数先将动量守恒方程即公式(4)中的重力项、粘性项和表面张力项离散，然后对计算域内全部粒子进行积分，显式得到粒子的估算速度和位置，其中粘性项需要使用采用如下多相粘度模型：
Figure BDA0003049813050000134
式中
μ——运动粘度系数m2/s
μij——粒子间调和运动粘度系数m2/s；
Figure BDA0003049813050000135
——速度m/s；
d——直径m；
re——粒子作用半径m；
nG——使用高斯核函数计算得到的粒子数密度；
Figure BDA0003049813050000136
——粒子j速度m/s；
Figure BDA0003049813050000137
——粒子i速度m/s；
G——高斯核函数；
rj——粒子j的位置矢量；
ri——粒子i的位置矢量；
步骤6：考虑到发生熔化的锆合金可以认为是一种纯净物，使用如下的焓值相变模型计算粒子的相态：
Figure BDA0003049813050000141
式中
T——需要计算的粒子温度K；
Ts——需要计算的粒子对应的熔点K；
h——当前粒子的焓值J/kg；
hs0——粒子开始熔化的焓值J/kg；
hs1——粒子结束熔化的焓值J/kg；
cp——粒子等压比热容J/(kg·K)；
通过步骤6的计算，能够模拟核反应堆板状燃料在最高温度超过燃料基体锆合金的熔点时的相变过程；计算得到每个粒子在不同时刻下的种类、焓值和温度，即得到锆合金基体和铀燃料微粒的相态、焓值和温度随时间的变化过程；
步骤7：得到粒子的相态后，对固体粒子之间是否发生碰撞进行判断，当固体粒子之间发生碰撞时，通过通过非连续介质模型的分析方法离散单元法(Discrete elementmethod简称DEM)对固体之间的碰撞作用进行如下的计算：
DEM模型的计算对象为刚性固体，当计算对象中存在大固体时，需要先解决将离散颗粒聚集成刚性物体统一运动的问题。采用如下公式初始化大固体运动速度
Figure BDA0003049813050000151
式中
Figure BDA0003049813050000152
——标记为ii的刚体的速度矢量m/s；
Figure BDA0003049813050000153
——i粒子当前的速度矢量m/s；
Figure BDA0003049813050000154
——标记为ii的刚体的角速度矢量rad/s；
mi——i粒子的质量kg；
Figure BDA0003049813050000155
——i粒子的位置矢量矢量；
Figure BDA0003049813050000156
——标记为ii的刚体的位置矢量；
Iii——标记为ii的刚体的转动惯量；
对于任意固体单元，采用如下公式对单元受力作用后产生的运动进行分析
Figure BDA0003049813050000157
Figure BDA0003049813050000158
式中
Figure BDA0003049813050000159
——标记为ii的刚体的加速度矢量；
Figure BDA00030498130500001510
——标记为ii的刚体的角加速度矢量；
Figure BDA00030498130500001511
——标记为ii的刚体的速度矢量；
Figure BDA00030498130500001512
——标记为ii的刚体的角速度矢量；
Figure BDA00030498130500001513
——标记为ii的刚体所受合力矢量；
Figure BDA0003049813050000161
——标记为ii的刚体所受合力矩矢量；
公式(8)和公式(9)中固体单元所受的合力与和力矩通过下式计算
Figure BDA0003049813050000162
Figure BDA0003049813050000163
式中
FX，jj→ii——标记为jj的刚体对标记为ii的刚体施加的法向合力；
en——法向弹性力；
dn——法向阻尼力；
kn——法向弹性系数；
cn——法向阻尼系数；
FY，jj→ii-—标记为jj的刚体对标记为ii的刚体施加的切向合力；
es——切向弹性力；
ds——切向阻尼力；
ks——切向弹性系数；
cs——切向阻尼系数；
进行一次DEM模型的计算后，固体单元mii在计算域中移动到了一个新的位置，和周围的单元之间产生新的接触力和接触力矩，再次计算固体单元mii所受的合力
Figure BDA0003049813050000171
和合力矩
Figure BDA0003049813050000172
重复计算流程，即得到每个固体单元及整个散粒体的运动形态，需要注意此时固体粒子的速度和位置均为估算值；
步骤8：通过核反应堆板状燃料熔化后的熔融物的不可压缩特性和动量守恒方程中的压力项联立得到如下压力泊松方程，隐式求解如下泊松方程
Figure BDA0003049813050000173
式中
ρ——粒子对应物质密度kg/m3；
P——压力Pa；
n0——初始粒子数密度常数；
n*——临时粒子数密度；
nk——当前时层的粒子数密度；
nk-1——上一时层的粒子数密度；
nk-1——下一时层的粒子数密度；
β——压力源项修正系数；
γ——压力源项修正系数；
通过步骤8的计算得到动量守恒方程中的压力项，使用该项所得值对步骤5和步骤7中得到的粒子估算速度和位置进行修正，得到当前时间步长全部粒子的真实速度和位置，这样就完成了一个时间步长的核反应堆板状燃料熔化流固耦合的模拟；
步骤9：为了准确分析模拟对象中固体和液体之间的耦合作用，还需要对步骤5中的粘性项计算和步骤8中的压力项计算进行处理，这是因为，固体和液体之间的作用通过下式表现
Figure BDA0003049813050000181
式中
Figure BDA0003049813050000182
——流体对固体的作用合力；
Figure BDA0003049813050000183
——流体对固体的压力作用；
Figure BDA0003049813050000184
——流体对固体的粘性作用；
在步骤5的粘性项计算中，通过多相粘度模型中使用的调和粘性系数来表现固体的物性，在计算时，这个粘性项只发生在固体和液体接触时，对固体内部不进行粘性项的计算；在步骤8的压力项计算中，流体粒子只和接触的固体粒子施加压力作用，后续压力作用将在DEM模型的计算中扩展至全部固体粒子。
综上，通过步骤1对核反应堆板状燃料熔化流固耦合问题中，锆合金材质的燃料包壳和燃料基体和二氧化铀燃料微粒的位置、速度和初始物性参数进行设定；通过步骤6模拟板燃料的锆合金基体在温度超过其熔点后发生相变形成熔融物，计算得到每个粒子在不同时刻下的种类、焓值和温度，得到熔融物的相态、焓值和温度随时间的变化过程；通过步骤5至步骤9模拟核反应堆板状燃料熔化流固耦合过程，模拟熔融物和未熔化的燃料微粒之间的相互作用，模拟燃料微粒之间的相互作用过程，计算得到熔融物和燃料微粒粒子的速度、位置和压力，即得到熔融物迁徙过程中的运动和压力变化过程；综合以上步骤，分析了核反应堆板状燃料熔化流固耦合问题，得到核熔化流固耦合过程中熔融物的位置、速度、压力、相态、温度、焓值随时间的变化，通过以上步骤对熔融物的迁徙行为展开机理性分析。

Claims (1)

1.一种核反应堆板状燃料熔化流固耦合无网格分析方法，其特征在于：步骤如下：
步骤1：对核反应堆板状燃料熔化流固耦合问题进行粒子建模；用不同类型的粒子代表核反应堆板状燃料的各个组成部分，0号粒子代表锆合金材质的燃料包壳和燃料基体，1号粒子代表弥散在燃料基体中的二氧化铀燃料微粒；对每一个粒子进行编号，所有的粒子均具有对应的物性参数，物性参数包括质量、密度、比热、熔沸点、温度、焓和初速度；粒子在笛卡尔坐标系内，依据核反应堆板状燃料的具体参数均匀布置，能够模拟核反应堆板状燃料的初始状态；
步骤2：使用链表法检索邻居粒子，具体方法是：将整个计算域均匀划分成正方形表格，正方形表格的边长与粒子最大作用半径相等；对于任意一个粒子而言，检索邻居粒子时只需检索包含粒子在内的表格以及包围此表格的共计9个表格，三维计算时，需要检索的就是由正方形表格围成的正方体，需要检索27个三维正方体；
步骤3：对时间步长进行稳定性分析，采用Courant-Friedrichs-Lewy条件简称CFL条件判断；CFL条件判断是差分方法保持计算稳定的必要条件，对于粒子法也适用，在移动粒子半隐式方法中采用下式进行判断
Figure FDA0003839039890000011
式中：
Δt——能够保持计算稳定的最小时间步长；
C——库朗常数；
l0——粒子初始布置时相邻两粒子间的距离；
umax——各时间步长中粒子速度最大值；
计算输入的初始条件中包含预设的时间步长大小，将输入的时间步长和使用CFL条件计算得到的能够保持计算稳定的最小时间步长进行比较，使用较大的时间步长作为当前时间步长；
步骤4：使用核函数表征粒子之间的相互作用程度，具体为在对动量守恒方程中的各项进行显式计算时，先使用核函数进行离散然后积分得到动量守恒方程中各项的具体值，使用的核函数如下所示：
Figure FDA0003839039890000021
式中
w(rij)——核函数值；
rij——粒子i与粒子j之间的距离；
re——粒子作用域半径；
步骤5：对核反应堆板状燃料熔化后的熔融物研究对象采用如下控制方程：
Figure FDA0003839039890000022
Figure FDA0003839039890000023
式中
ρ——粒子对应物质的密度kg/m3；
t——时间s；
P——压力Pa；
μ——动力粘度系数N·s/m2；
Figure FDA0003839039890000031
——速度矢量m/s；
Figure FDA0003839039890000032
——表面张力矢量N/kg；
Figure FDA0003839039890000033
——重力加速度m/s2；
公式(3)为连续性方程，针对的对象为不可压缩流体，公式(4)为动量守恒方程，从中看出分析对象只受到重力这一外力作用，对于任意一个粒子，通过核函数先将动量守恒方程即公式(4)中的重力项、粘性项和表面张力项离散，然后对计算域内全部粒子进行积分，显式得到粒子的估算速度和位置，其中粘性项需要使用采用如下多相粘度模型：
Figure FDA0003839039890000034
式中
μ——动力粘度系数N·s/m2；
μij——粒子间动力粘度系数m2/s；
Figure FDA0003839039890000035
——速度m/s；
d——直径m；
re——粒子作用半径m；
nG——使用高斯核函数计算得到的粒子数密度；
Figure FDA0003839039890000036
——粒子j速度m/s；
Figure FDA0003839039890000037
——粒子i速度m/s；
G——高斯核函数；
rj——粒子j的位置矢量；
ri——粒子i的位置矢量；
步骤6：使用如下的焓值相变模型计算粒子的相态：
Figure FDA0003839039890000041
式中
T——需要计算的粒子温度K；
Ts——需要计算的粒子对应的熔点K；
h——当前粒子的焓值J/kg；
hs0——粒子开始熔化的焓值J/kg；
hs1——粒子结束熔化的焓值J/kg；
cp——粒子等压比热容J/(kg·K)；
通过步骤6的计算，能够模拟核反应堆板状燃料在最高温度超过燃料基体锆合金的熔点时的相变过程；计算得到每个粒子在不同时刻下的种类、焓值和温度，即得到锆合金基体和铀燃料微粒的相态、焓值和温度随时间的变化过程；
步骤7：得到粒子的相态后，通过非连续介质模型的分析方法离散单元法DEM对固体之间的碰撞作用进行如下的计算：
非连续介质模型的分析方法离散单元法DEM的计算对象为刚性固体，当计算对象中存在大固体时，需要先解决将离散颗粒聚集成刚性物体统一运动的问题；采用如下公式初始化大固体运动速度；
Figure FDA0003839039890000042
式中
Figure FDA0003839039890000051
——标记为ii的刚体的速度矢量m/s；
Figure FDA0003839039890000052
——粒子i当前的速度矢量m/s；
Figure FDA0003839039890000053
——标记为ii的刚体的角速度矢量rad/s；
mi——粒子i的质量kg；
Figure FDA0003839039890000054
——粒子i的位置矢量；
Figure FDA0003839039890000055
——标记为ii的刚体的位置矢量；
Iii——标记为ii的刚体的转动惯量；
对于任意固体单元mii，采用如下公式对单元受力作用后产生的运动进行分析；
Figure FDA0003839039890000056
Figure FDA0003839039890000057
式中
Figure FDA0003839039890000058
——标记为ii的刚体的加速度矢量；
Figure FDA0003839039890000059
——标记为ii的刚体的角加速度矢量；
Figure FDA00038390398900000510
——标记为ii的刚体的速度矢量；
Figure FDA00038390398900000511
——标记为ii的刚体的角速度矢量；
Figure FDA00038390398900000512
——标记为ii的刚体所受合力矢量；
Figure FDA00038390398900000513
——标记为ii的刚体所受合力矩矢量；
公式(8)和公式(9)中固体单元所受的合力与和力矩通过下式计算
Figure FDA0003839039890000061
Figure FDA0003839039890000062
式中
FX,jj→ii——标记为jj的刚体对标记为ii的刚体施加的法向合力；
en——法向弹性力；
dn——法向阻尼力；
kn——法向弹性系数；
cn——法向阻尼系数；
FY,jj→ii——标记为jj的刚体对标记为ii的刚体施加的切向合力；
es——切向弹性力；
ds——切向阻尼力；
ks——切向弹性系数；
cs——切向阻尼系数；
这样，固体单元mii就移动到了一个新的位置，并产生新的接触力和接触力矩，再次计算固体单元mii所受的合力
Figure FDA0003839039890000063
和合力矩
Figure FDA0003839039890000064
重复计算流程，即得到每个固体单元及整个散粒体的运动形态，此时固体粒子的速度和位置均为估算值；
步骤7：通过核反应堆板状燃料熔化后的熔融物的不可压缩特性和动量守恒方程中的压力项联立得到如下压力泊松方程，隐式求解如下的压力泊松方程
Figure FDA0003839039890000071
式中
ρ——粒子对应物质密度kg/m3；
P——压力Pa；
n0——初始粒子数密度常数；
n*——临时粒子数密度；
nk——当前时层的粒子数密度；
nk-1——上一时层的粒子数密度；
nk-1——下一时层的粒子数密度；
β——压力源项修正系数；
γ——压力源项修正系数；
通过步骤8的计算得到动量守恒方程中的压力项，使用该压力项所得值对步骤5和步骤7中得到的粒子估算速度和位置进行修正，得到当前时间步长全部粒子的真实速度和位置，这样就完成了一个时间步长的核反应堆板状燃料熔化流固耦合的模拟；
步骤9：为了准确分析模拟对象中固体和液体之间的耦合作用，还需要对步骤5中的粘性项计算和步骤8中的压力项计算进行处理，这是因为，固体和液体之间的作用通过下式表现
Figure FDA0003839039890000072
式中
Figure FDA0003839039890000081
——流体对固体的作用合力；
Figure FDA0003839039890000082
——流体对固体的压力作用；
Figure FDA0003839039890000083
——流体对固体的粘性作用；
在步骤5的粘性项计算中，通过多相粘度模型中使用的调和粘性系数来表现固体的物性，在计算时，这个粘性项只发生在固体和液体接触时，对固体内部不进行粘性项的计算；在步骤8的压力项计算中，流体粒子只和接触的固体粒子施加压力作用，后续压力作用将在非连续介质模型的分析方法离散单元法DEM的计算中扩展至全部固体粒子；
综上，通过步骤1对核反应堆板状燃料熔化流固耦合问题中，锆合金材质的燃料包壳和燃料基体和二氧化铀燃料微粒的位置、速度和初始物性参数进行设定；通过步骤6模拟板燃料的锆合金基体在温度超过其熔点后发生相变形成熔融物，计算得到每个粒子在不同时刻下的种类、焓值和温度，得到熔融物的相态、焓值和温度随时间的变化过程；通过步骤5至步骤9模拟核反应堆板状燃料熔化流固耦合过程，模拟熔融物和未熔化的燃料微粒之间的相互作用，模拟燃料微粒之间的相互作用过程，计算得到熔融物和燃料微粒粒子的速度、位置和压力，即得到熔融物迁徙过程中的运动和压力变化过程；综合以上步骤，分析了核反应堆板状燃料熔化流固耦合问题，得到熔化流固耦合过程中熔融物的位置、速度、压力、相态、温度、焓值随时间的变化，通过以上步骤对熔融物的迁徙行为展开机理性分析。
CN202110484617.5A 2021-04-30 2021-04-30 一种核反应堆板状燃料熔化流固耦合无网格分析方法 Active CN113192567B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110484617.5A CN113192567B (zh) 2021-04-30 2021-04-30 一种核反应堆板状燃料熔化流固耦合无网格分析方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110484617.5A CN113192567B (zh) 2021-04-30 2021-04-30 一种核反应堆板状燃料熔化流固耦合无网格分析方法

Publications (2)

Publication Number Publication Date
CN113192567A CN113192567A (zh) 2021-07-30
CN113192567B true CN113192567B (zh) 2022-12-09

Family

ID=76983329

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110484617.5A Active CN113192567B (zh) 2021-04-30 2021-04-30 一种核反应堆板状燃料熔化流固耦合无网格分析方法

Country Status (1)

Country Link
CN (1) CN113192567B (zh)

Families Citing this family (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113836835B (zh) * 2021-08-28 2022-06-21 西安交通大学 一种核反应堆燃料熔化熔融物迁徙行为流固耦合分析方法
CN114757123B (zh) * 2022-04-20 2023-06-20 西安交通大学 一种用于板形核燃料堆芯的跨维度流固耦合分析方法
CN115099172B (zh) * 2022-07-08 2024-03-12 西安交通大学 一种用于熔融物碎片床形成过程特性分析的方法
CN115630559B (zh) * 2022-12-08 2023-03-10 中国人民解放军国防科技大学 一种基于粒子网格适配算法的流固耦合方法以及装置

Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2001039204A2 (en) * 1999-11-24 2001-05-31 Impulse Devices, Inc. Shaped core cavitation nuclear reactor
US6627445B1 (en) * 1999-12-15 2003-09-30 Uop Llc Process for simultaneously evaluating a plurality of catalysts
CN109670239A (zh) * 2018-12-18 2019-04-23 北京应用物理与计算数学研究所 基于pin-by-pin模型的压水堆生产同位素模拟方法及系统
CN110598324A (zh) * 2019-09-12 2019-12-20 西安交通大学 一种核反应堆弥散型板型燃料元件堆芯流固耦合计算方法
CN110867220A (zh) * 2019-11-07 2020-03-06 西安交通大学 利用粒子网格混合法研究堆芯内共晶反应及高温熔化行为的方法
CN111504588A (zh) * 2020-05-29 2020-08-07 上海交通大学 燃料棒束两相流动流固耦合试验回路
CN111832214A (zh) * 2020-06-29 2020-10-27 西安交通大学 基于粒子法的核反应堆严重事故碎片床熔化过程模拟方法
CN112102894A (zh) * 2020-09-04 2020-12-18 西安交通大学 基于粒子法的核反应堆堆芯材料熔池演变特性分析方法

Family Cites Families (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105023621B (zh) * 2015-06-12 2017-11-10 陈安海 快堆型耦合核反应的实施方法及其核反应堆
CN108491619B (zh) * 2018-03-19 2020-05-26 浙江大学 基于物理与非物理混合的复杂场景流固耦合高效模拟方法
CN110321641B (zh) * 2019-07-08 2020-08-04 西安交通大学 基于粒子法的熔融物与混凝土相互作用分析方法

Patent Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2001039204A2 (en) * 1999-11-24 2001-05-31 Impulse Devices, Inc. Shaped core cavitation nuclear reactor
US6627445B1 (en) * 1999-12-15 2003-09-30 Uop Llc Process for simultaneously evaluating a plurality of catalysts
CN109670239A (zh) * 2018-12-18 2019-04-23 北京应用物理与计算数学研究所 基于pin-by-pin模型的压水堆生产同位素模拟方法及系统
CN110598324A (zh) * 2019-09-12 2019-12-20 西安交通大学 一种核反应堆弥散型板型燃料元件堆芯流固耦合计算方法
CN110867220A (zh) * 2019-11-07 2020-03-06 西安交通大学 利用粒子网格混合法研究堆芯内共晶反应及高温熔化行为的方法
CN111504588A (zh) * 2020-05-29 2020-08-07 上海交通大学 燃料棒束两相流动流固耦合试验回路
CN111832214A (zh) * 2020-06-29 2020-10-27 西安交通大学 基于粒子法的核反应堆严重事故碎片床熔化过程模拟方法
CN112102894A (zh) * 2020-09-04 2020-12-18 西安交通大学 基于粒子法的核反应堆堆芯材料熔池演变特性分析方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
Analysis on the behavior of dispersed plate-type fuel based on fluid-solid coupling method;Chenxi Li等;《Progress in Nuclear Energy》;20200527;1-11 *
反应堆冷却剂泵多相流及流固耦合振动特性研究;陈向阳;《中国优秀硕士学位论文全文数据库 工程科技Ⅱ辑》;20111215;第2011年卷(第S1期);C029-128 *
矩形通道中单块柔性平板流固耦合实验研究;赵玉静;《中国优秀硕士学位论文全文数据库 工程科技Ⅱ辑》;20100915;第2010年卷(第9期);C040-7 *

Also Published As

Publication number Publication date
CN113192567A (zh) 2021-07-30

Similar Documents

Publication Publication Date Title
CN113192567B (zh) 一种核反应堆板状燃料熔化流固耦合无网格分析方法
Guo et al. Numerical simulation and parametric study on new type of high temperature latent heat thermal energy storage system
CN111832214B (zh) 基于粒子法的核反应堆严重事故碎片床熔化过程模拟方法
Eberhard et al. Simulation of cutting processes using mesh-free Lagrangian particle methods
Yamashita et al. A numerical simulation method for molten material behavior in nuclear reactors
Duan et al. A multiphase MPS method coupling fluid–solid interaction/phase-change models with application to debris remelting in reactor lower plenum
CN115587551B (zh) 一种高超声速飞行器防热结构烧蚀行为多尺度预测方法
Yeon et al. CFD analysis of core melt spreading on the reactor cavity floor using ANSYS CFX code
Kim et al. Structural assessment of reactor pressure vessel under multi-layered corium formation conditions
Takahashi et al. Analysis of hemispherical vessel ablation failure involving natural convection by MPS method with corrective matrix
Ye et al. Numerical investigation of the spreading and heat transfer characteristics of ex-vessel core melt
CN113836835B (zh) 一种核反应堆燃料熔化熔融物迁徙行为流固耦合分析方法
Zohrabi Chakaneh et al. Experimental and numerical characterization of drop impact on a hydrophobic cylinder
Mukherjee et al. Numerical study of an evaporating meniscus on a moving heated surface
Li et al. Pebble flow and coolant flow analysis based on a fully coupled multiphysics model
Rao et al. The lattice Boltzmann investigation for the melting process of phase change material in an inclined cavity
Alluhydan et al. On planar impacts of cylinders and balls
Park et al. Three-dimensional modeling of debris mixing and sedimentation in severe accidents using the moving particle semi-implicit method coupled with rigid body dynamics
Xiao et al. An improved MPS-DEM numerical model for fluid–solid coupling problem in nuclear reactor
Terrazas et al. The CFD Modeling of the Water Braking Phenomena for the Holloman High-Speed Test Track
Zhao et al. Simulation of melt spreading over dry substrates with the moving particle Semi-implicit method
Hao et al. Double-diffusive natural convection in a nuclear waste repository
Matsuo et al. Study on application of artificial neural network to debris bed coolability calculations for sodium-cooled fast reactors
Li et al. Estimation of debris relocation and structure interaction in the pedestal of Fukushima Daiichi Nuclear Power Plant Unit-3 with Moving Particle Semi-implicit (MPS) method
Chatzikyriakou et al. Three dimensional modeling of the hydrodynamics of oblique droplet-hot wall interactions during the reflood phase after a LOCA

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