CN110867220A - 利用粒子网格混合法研究堆芯内共晶反应及高温熔化行为的方法 - Google Patents

利用粒子网格混合法研究堆芯内共晶反应及高温熔化行为的方法 Download PDF

Info

Publication number
CN110867220A
CN110867220A CN201911083355.0A CN201911083355A CN110867220A CN 110867220 A CN110867220 A CN 110867220A CN 201911083355 A CN201911083355 A CN 201911083355A CN 110867220 A CN110867220 A CN 110867220A
Authority
CN
China
Prior art keywords
particle
particles
equation
formula
mass
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.)
Pending
Application number
CN201911083355.0A
Other languages
English (en)
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 CN201911083355.0A priority Critical patent/CN110867220A/zh
Publication of CN110867220A publication Critical patent/CN110867220A/zh
Pending legal-status Critical Current

Links

Images

Classifications

    • 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/70—Machine learning, data mining or chemometrics

Landscapes

  • Chemical & Material Sciences (AREA)
  • Engineering & Computer Science (AREA)
  • Computing Systems (AREA)
  • Theoretical Computer Science (AREA)
  • Crystallography & Structural Chemistry (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Health & Medical Sciences (AREA)
  • Chemical Kinetics & Catalysis (AREA)
  • Analytical Chemistry (AREA)
  • Artificial Intelligence (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Data Mining & Analysis (AREA)
  • Databases & Information Systems (AREA)
  • Evolutionary Computation (AREA)
  • General Health & Medical Sciences (AREA)
  • Medical Informatics (AREA)
  • Software Systems (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

一种利用粒子网格混合法研究堆芯内共晶反应及高温熔化行为的方法，主要步骤如下：1、针对具体问题设定粒子初始布置；2、对初始粒子位置进行网格划分；3、使用有限体积法求解质量扩散方程及能量方程，得到下一时层的质量和温度解；4、使用移动粒子半隐式法求解动量方程，得到粒子的位置及速度；5、迭代求解步骤2～4，直到求解得到所需时刻的结果；6、输出计算结果，后处理得到堆芯内的热量分布、熔融物分布及熔融物迁移情况。本发明方法具有移动粒子半隐式法精确处理自由表面及模拟相变的优点，同时具有有限体积法计算的高精度性，能够通过此方法来研究堆芯熔化的进程，为核电厂反应堆严重事故堆芯早期熔化安全特性的研究提供重要依据。

Description

利用粒子网格混合法研究堆芯内共晶反应及高温熔化行为的 方法
技术领域
本发明涉及核电厂严重事故堆芯内共晶反应及高温熔化行为特性研究技术领域，具体涉及一种利用粒子网格混合法研究堆芯内共晶反应及高温熔化行为的方法。
背景技术
近期我国核电发展迅速。为了改变我国现有的能源结构，缓解对环境的压力、实现减排目标，国家制定了庞大的核电中长期发展规划。日本福岛核电站事故对我国核电安全敲响了警钟。事故发生后，我国政府高度重视我国核电的安全问题。开展反应堆严重事故机理及关键技术研究可以了解和熟悉严重事故过程，为核电站制定完善的严重事故环节策略提供技术支持，制定一些防止严重事故发生的措施，减少严重事故发生的可能性和概率。
但严重事故本身是一个多相态、多组分的复杂物理和化学过程，严重事故的燃料熔化和烛化过程是严重事故进程的关键环节：它能提供堆内及堆外严重事故现象的初始条件，诱发对压力容器完整性的威胁，并能决定裂变产物及氢气源项。堆内燃料熔融物的迁移过程相关知识同样非常重要，它为堆芯冷却(堆芯再淹没)及压力容器失效分析提供帮助。堆芯熔化过程是一个非连贯的过程，它会导致堆芯材料在不同的温度下熔化或者液化。氧化锆包壳层能够明显地延缓熔融物(U，Zr，O)及其它熔化在内的物质的迁移。这些熔融物在不同的温度、不同的轴向位置受到冷却固化，并由陶瓷的包壳覆盖金属物质。如果事故不能得到缓解，这些包壳可能会失效并导致熔融物材料的再迁移，碎片进入下腔室。在预测堆芯熔化的这一系列过程中，存在很多的不确定性。
国际上关于严重事故燃料熔化和烛化的研究，可以追溯到上世纪60年代，核电工业刚刚兴起的时候。但是由于当时在使用核电的国家中，并不需要进行相关研究就可以取得核电执照，所以相关的研究并没有取得太大的进展。直到70年代末，美国发生三哩岛事故，给严重事故燃料熔化和烛化研究提供了直接的经验和丰富的素材。随后，德国进行了CORA实验的研究，法国CEA进行了PHEBUS实验的研究，积累了丰富经验。而在程序开发方面，包括MELCOR，MAAP，SCDAP/RELAP5等都吸收并包含了国际上堆芯熔融物材料的熔化和烛化研究的模型和成果。
移动粒子半隐式(Moving Particle Semi-implicit MPS)方法是日本东京大学Koshizuka教授上世纪90年代末提出的一种新型无网格方法，该方法是一种基于核近似的配点型纯Lagrange方法，其基本思想是以一系列离散的移动粒子来表示宏观的流体，流体运动控制方程中的各微分算子以粒子相互作用的形式表达。粒子相互作用模型包括一般标量的梯度离散模型、速度等矢量离散的散度模型及用于扩散项离散的拉普拉斯模型，粒子只与其相邻限定区域内的粒子相互作用，相互依赖强度通过核函数来实现。由于MPS方法采用移动粒子和Lagrange方法直接计算流体行为，使其在捕捉多相流体动力学相界面上具有独特的优势。目前该方法已在一些工程领域对大变形问题进行数值模拟方面得到了初步应用。在反应堆严重事故方面，MPS方法也得到了初步的应用。由于其计算得特性，其能详细描述严重事故中堆芯熔化和下封头扩展行为，并将其动态变化过程很形象得展示出来。日本学者Shibata对射流断裂进行了很详细的研究，Hirokazu对FCI(熔融物与冷却剂作用过程)的射流现象进行了模拟，早稻田大学Chen对管道中熔融物冷凝行为进行了数值模拟研究。西安交通大学的Chen和Li对燃料棒材料的高温共晶反应及消熔反应进行了纯粒子法的数值模拟研究。从上述的阐述中可以发现，堆芯材料在发生高温熔化时、熔融物流动时存在较多的流动自由面，使用移动粒子半隐式方法无疑是可以很精确地捕捉熔融物的流动形态，但是移动粒子半隐式在计算传热问题及质量扩散问题时，在计算域的边界及其附近的粒子会因为计算支持域被截断而使结果存在较大误差，这也必然会导致堆芯熔化行为研究结果的不准确性。因为堆芯熔化时，熔融物流动形态多变，在流动表面进行一系列的例如镜像粒子补偿等方法并不可取，有限体积法在计算传热问题及质量扩散问题时精度比移动粒子半隐式法的精度要好。基于此种考虑，本发明综合移动粒子半隐式方法和有限体积法提出了一种利用粒子网格混合法研究堆芯内共晶反应及高温熔化行为的方法
发明内容
为了克服移动粒子半隐式方法在计算堆芯内共晶反应及高温熔化行为时存在的质量扩散计算及传热计算的不精确问题，本发明提供一种利用粒子网格混合法研究堆芯内共晶反应及高温熔化行为的方法，该方法能够对堆芯内多相态、多自由面及存在质量扩散的堆芯内共晶反应及高温熔化行为进行研究，获得堆芯内部的热量分布、熔融物分布及熔融物迁移情况，揭示反应堆堆芯燃料熔化和熔融物扩展行为，为核电站制定完善的严重事故环节策略提供技术支持。
为了实现上述目的，本发明采取了以下的技术方案予以实施：
一种利用粒子网格混合法研究堆芯内共晶反应及高温熔化行为的方法，步骤如下：
步骤1：对反应堆堆芯进行粒子建模，不同的物质使用不同种类的粒子表示，粒子初始状态包括粒子的位置、速度、压力和温度，对应标记为ri、ui、pi和Ti，粒子初始布置时相邻两粒子间的距离记为l0，每个粒子内按所含元素标记对应元素的质量，如粒子i内包含A、B两种元素，则其质量对应标记为miA、miB，粒子i的质量mi使用公式(1)计算，
mi＝miA+miB 公式(1)
粒子i的密度ρi使用公式(2)计算，
Figure BDA0002264623960000031
对于熔化焓等则通过公式(3)求得，
Figure BDA0002264623960000032
步骤2：根据各粒子的位置划分网格：在本步骤中，对于每一个粒子i，设置粒子作用域半径lr，若粒子i的坐标位置为(Ix，Iy，Iz)，设定粒子i周围的某一粒子为粒子j，坐标位置为(Jx，Jy，Jz)，粒子i与粒子j的距离为lij，使用公式(4)计算，
Figure BDA0002264623960000033
若满足lij≤lr，则将该就j粒子标记为i粒子的有效作用粒子；将所有满足条件的m个j粒子进行汇总并标记为j1，j2，……，jm，然后对其中每3个粒子进行组合，按照排列组合法，将有
Figure BDA0002264623960000041
种组合，
Figure BDA0002264623960000042
由公式(5)进行计算，
Figure BDA0002264623960000043
对每个组合中的粒子进行考察，取其中一个组合的3个粒子，这3个粒子必须满足能且只能构成一个平面的要求，将取出的3个粒子标记为j1，j2，j3，粒子中心分别标记为J1，J2，J3，空间坐标对应为(Jx1，Jy1，Jz1)，(Jx2，Jy2，Jz2)，(Jx3，Jy3，Jz3)，连接J1、J2，连接J1、J3，连接J2、J3，J1、J2、J3三点构成平面的其中一个充要条件是直线J1J2与直线J1J3之间存在一个不为0的夹角，则只需满足公式(6)，
Figure BDA0002264623960000044
式中，||J1J2||、||J1J3||、||J2J3||分别为线段J1J2、J1J3、J2J3的长度，使用公式(7)进行计算，
Figure BDA0002264623960000045
对每一个i粒子，针对其
Figure BDA0002264623960000046
种组合进行判断，若组合内的3个粒子能够满足构成唯一平面的要求则留下该组合，随后针对留下的组合，考察i粒子与组合内3个粒子的关系，i粒子与组合中三个粒子必须能够构成四面体；对于留下的每个组合而言，这3个粒子已经满足构成唯一平面，那么要使这i粒子与组合中的3个粒子能够构成四面体，只需满足i粒子不在组合中的3个粒子所构成的平面上即可，添加i粒子，粒子中心标记为I，空间坐标定义为(Ix，Iy，Iz)，若要使I、J1、J2、J3四个点构成四面体，则只需要满足i粒子到平面J1J2J3的距离为正数，使用公式(8)进行判定，
Figure BDA0002264623960000051
式中：
Figure BDA0002264623960000052
——i粒子到平面J1J2J3的距离m；
A、B、C、D——平面J1J2J3的平面方程的四个系数，通过公式(9)获得，
Figure BDA0002264623960000053
i粒子与j1、j2、j3粒子构成的平面J1J2J3的距离若满足公式(8)，则说明i粒子能够与组合内3个粒子构成的三角形J1J2J3构成四面体IJ1J2J3，其四个面分别为三角形IJ1J2、三角形IJ1J3、三角形IJ2J3和三角形J1J2J3，面积对应标记为
Figure BDA0002264623960000054
三角形的面积使用如公式(10)所示的海伦公式计算
Figure BDA0002264623960000055
式中，
S——三角形的面积；
la、lb、lc——三角形的三边长；
p——三角形的半周长，即p＝(la+lb+lc)/2；
四面体IJ1J2J3的体积通过公式(11)进行计算，
Figure BDA0002264623960000056
四面体IJ1J2J3中四个顶点I、J1、J2、J3分别对四面体总体积的分属控制体积，对应标记为VI、
Figure BDA0002264623960000057
这4个体积有公式(12)所示的关系，
Figure BDA0002264623960000058
各控制体积由公式(13)进行计算，
Figure BDA0002264623960000061
除此之外，四面体IJ1J2J3的四个点的分属控制体积与其相对的三角形面积是成反比的，即有公式(14)所示的关系，
Figure BDA0002264623960000062
由以上公式联立求解即求得四面体IJ1J2J3的四个点的分属控制体积；
然后将对四面体IJ1J2J3中四个面的三角形进行面积划分，以三角形J1J2J3为例，在三角形J1J2J3中，三角形J1J2J3的面积
Figure BDA0002264623960000063
通过公式(10)海伦公式求出；点N12为线段J1J2的中点，点N13为线段J1J3的中点，点N23为线段J2J3的中点；此时需要在三角形中寻找P点，对于三角形J1J2J3而言，将其标记为Px，连接Px、J1，连接Px、J2，连接Px、J3，连接Px、N12，连接Px、N13，连接Px、N23，可以将三角形J1J2J3划分成6个小三角形，分别标记为三角形J1PxN12，三角形J1PxN13，三角形J2PxN12，三角形J2PxN23，三角形J3PxN13，三角形J3PxN23；因为点N13为线段J1J3的中点，所以三角形J1PxN13与三角形J3PxN13面积相等，标记为
Figure BDA0002264623960000064
同样的，有
Figure BDA0002264623960000065
对于各个三角形的面积，有如公式(15)和公式(16)所示的几何关系，
Figure BDA0002264623960000066
Figure BDA0002264623960000067
三角形的三个角的角度由类似于公式(6)所示的余弦定理求出；由公式(15)和公式(16)求出三角形J1J2J3所分成的6个小三角形的面积；同理也能求得四面体IJ1J2J3另外三个面所分成的共计18个小三角形的面积；
在四面体IJ1J2J3中寻找一点Qx，对于顶点I的控制体积而言，连接QxPx1，QxPx2，QxPx3，QxN1，QxN2，QxN3，QxI，其中Px1、Px2、Px3、分别对应为三角形IJ1J2、三角形IJ2J3、三角形IJ1J3所寻找的P点，N1、N2、N3则分别为线段IJ1、IJ2、IJ3的中点，至此点I的控制体积分成了六个四面体，分别为四面体IPx1N1Qx，四面体IPx1N2Qx，四面体IPx2N2Qx，四面体IPx2N3Qx，四面体IPx3N1Qx和四面体IPx3N3Qx，其体积对应标记为
Figure BDA0002264623960000071
Figure BDA0002264623960000072
和
Figure BDA0002264623960000073
对于四面体的体积，有公式(17)所示的关系，
Figure BDA0002264623960000074
设点Qx距离平面J1J2J3、平面IJ1J2，平面IJ2J3，平面IJ1J3的距离分别为H0、H1、H2、H3，对于顶点I的控制体积可列出如公式(18)的方程，
Figure BDA0002264623960000075
同理对于顶点J1、J2、J3也列出相似的方程，联立这4个方程及公式(13)能够求出H0、H1、H2、H3的值；
至此，在I、J1、J2、J3组成的四面体中，J1对I起影响的点I的控制体积为
Figure BDA0002264623960000076
记为
Figure BDA0002264623960000077
J2对I起影响的点I的控制体积为
Figure BDA0002264623960000078
记为
Figure BDA0002264623960000079
J3对I起影响的点I的控制体积为
Figure BDA00022646239600000710
记为
Figure BDA00022646239600000711
所以，在I、J1、J2、J3组成的四面体系统中，点J1对I的有效传输面积标记为
Figure BDA00022646239600000712
点J2对I的有效传输面积标记为
Figure BDA00022646239600000713
点J1对I的有效传输面积标记为
Figure BDA00022646239600000714
通过公式(19)计算，
Figure BDA00022646239600000715
式中，
Figure BDA00022646239600000716
——I点到J1点的距离；
Figure BDA00022646239600000717
——I点到J2点的距离；
Figure BDA00022646239600000718
——I点到J3点的距离；
以上的划分网格的方法描述是针对其中一个四面体系统而言的，在计算中需要对i粒子的所有符合的四面体网格进行计算得到有效传输面积数据之后才能开始进行质量扩散计算及温度计算；
步骤3：在质量扩散的计算中，使用的传质方程为菲克第二定律，方程如公式(20)所示，对于粒子i内的A、B元素分别列出方程，
Figure BDA0002264623960000081
式中，
m——粒子内的某元素质量kg；
t——时间s；
Ddiff——粒子内的某元素的质量扩散系数m2/s；
公式(21)给出质量扩散方程的Δt时间步长的积分形式，
Figure BDA0002264623960000082
式中，
m——粒子内的某元素质量kg；
t——时间s；
Δt——时间步长s；
V——体积m3；
CV——控制体积m3；
Ddiff——粒子内的某元素的质量扩散系数m2/s；
对于粒子i而言，其已经通过步骤2划分为若干个四面体网格，所以公式(21)能够离散成公式(22)的形式，
Figure BDA0002264623960000083
公式(22)等式的左边对i粒子划分成的多个四面体的某一元素的质量变化进行加和，等式的右边则是i粒子的该元素的质量变化率，计算中不考虑粒子的体积变化；式中，
Sj→i——步骤2中网格划分完成后计算的j粒子对i粒子的有效传输面积；
Ddiff——某元素的质量扩散系数m2/s；
rj——j粒子的坐标矢量；
ri——i粒子的坐标矢量；
Figure BDA0002264623960000091
——j粒子中x元素在n时层的质量kg；
Figure BDA0002264623960000092
——i粒子中x元素在n时层的质量kg；
Figure BDA0002264623960000093
——i粒子中x元素在n+1时层的质量，即在n时层的时间步长Δt后的x元素质量kg；
i粒子在n+1时层的质量
Figure BDA0002264623960000094
使用公式(23)计算，
Figure BDA0002264623960000095
在温度计算中，能量守恒方程如公式(24)所示，
Figure BDA0002264623960000096
式中，
T——温度K；
t——时间s
κ——导热系数W/(m2·K)；
ρ——密度kg/m3；
cp——比定压热容J/(kg·K)；
ST——温度的源项；
公式(25)给出能量方程的Δt时间步长的积分形式，
Figure BDA0002264623960000101
式中，
t——时间s；
Δt——时间步长s；
T——温度K；
V——体积m3；
CV——控制体积m3；
κ——导热系数W/(m2·K)；
ρ——密度kg/m3；
cp——比定压热容J/(kg·K)；
ST——温度的源项；
对于粒子i而言，其已经通过步骤2划分为若干个四面体网格，所以公式(25)能够离散成公式(26)的形式，
Figure BDA0002264623960000102
公式(26)等式的左边对i粒子划分成的多个四面体的能量变化进行加和，等式的右边则是i粒子的能量变化率；式中：
Sj→i——步骤2中网格划分完成后计算的j粒子对i粒子的有效传输面积；
κ——导热系数W/(m2·K)；
rj——j粒子的坐标矢量；
ri——i粒子的坐标矢量；
Figure BDA0002264623960000103
——j粒子在n时层的温度K；
Ti n——i粒子在n时层的温度K；
ρi——i粒子的密度kg/m3；
cpi——i粒子的比定压热容J/(kg·K)；
Ti n+1——i粒子在n+1时层的温度，即在n时层的时间步长Δt后的温度K；
求解公式(26)即得到i粒子在时间步长Δt之后的温度Ti n+1；
至此，i粒子在时间步长Δt之后的温度Ti n+1及其质量
Figure BDA0002264623960000113
已经求解得到，并且对于i粒子内的各元素进行的单独质量求解极易得到各元素在粒子内的质量百分比，若粒子内只含有A、B两种元素，i粒子中A元素的质量百分比由公式(27)计算，
Figure BDA0002264623960000111
式中，
CiA——i粒子中A元素的质量百分比；
miA——i粒子中A元素的质量kg；
miB——i粒子中B元素的质量kg；
此时使用A、B元素的二元相图进行粒子状态判定，在A-B二元相图中，横坐标为A元素的质量百分比，纵坐标为粒子的温度，曲线ApQSBp为固相线，曲线ApQLBp为液相线，粒子在多边形O1O2BpQSAp围成的区域内状态为固态，在多边形ApQLBpO3O4围成的区域内状态为液态，在多边形ApQSBpQL围成的区域内状态为固液两相混合状态；若在某一时刻，i粒子的温度为Tk，i粒子中A元素的质量百分比为Ck，在A-B二元相图中作温度Tk的等温线，做质量百分比Ck的等值线，交于点Qp，此时点Qp位于A-B二元相图的液相区域，则对应的，粒子相态为液相；对位于固液两相混合状态区域内的粒子，需要根据相图中的杠杆定律对其成分进行分析；同样的方法，对于温度为Tpx和粒子内A元素质量百分比为Cpx的粒子k，判断其对应的Qpx点位于固液两相区内，作线段TpxQpx的延长线交固相线和液相线于QS和QL，再过QS和QL作垂直线交横坐标于CS和CL点，则粒子k中的固体和液体质量之间的关系使用公式(28)计算得到，
Figure BDA0002264623960000112
式中，
mS——k粒子中处于固态的质量kg；
mL——k粒子中处于液态的质量kg；
Cpx——k粒子中A元素的质量百分比；
CS——点CS对应的横坐标数值；
CL——点CL对应的横坐标数值；
步骤4：使用移动粒子半隐式方法计算动量方程；动量方程如公式(29)所示
Figure BDA0002264623960000121
式中，
u——速度m/s；
t——时间s；
ρ——密度kg/m3；
p——压力Pa
ν——运动粘度m2/s；
f——外力加速度m2/s；
在此方程组中，动量方程中的粘性项及源项显式计算，其中粘性项的Laplace算子由公式(30)离散求解，
Figure BDA0002264623960000122
式中，
d——计算的维度；
n0——初始粒子数密度，通过
Figure BDA0002264623960000123
计算，
Figure BDA0002264623960000124
是核函数，使用公式(31)计算，只需将括号内的自变量修改为
Figure BDA0002264623960000125
即可，
Figure BDA0002264623960000126
为j粒子初始时刻的坐标矢量，
Figure BDA00022646239600001310
为i粒子初始时刻的坐标矢量；
λ——是i粒子的邻居粒子和i粒子间距的方程，使用公式(32)计算，
Figure BDA0002264623960000131
式中，
re——粒子作用半径，位于半径内的粒子才会与i粒子发生相互作用，，re＝3.1l0；
Figure BDA0002264623960000132
速度和位置的估算值
Figure BDA0002264623960000133
由公式(33)进行计算，
Figure BDA0002264623960000134
式中，
Figure BDA0002264623960000135
——i粒子在上一时间步的速度矢量m/s；
ri n——i粒子在上一时间步的位置矢量m；
Ti n——i粒子在上一时间步的温度K；
根据估算得到的粒子位置，通过公式(34)计算此时的粒子数密度n* i，只需要将公式中的|rj-ri|使用位置的估算值
Figure BDA0002264623960000136
替代即可，由公式(34)所示的压力Poisson方程求解得到此时各粒子的压力值pi；
Figure BDA0002264623960000137
Figure BDA0002264623960000138
公式(35)等式的左边使用如同公式(30)的形式离散，式中，
Figure BDA0002264623960000139
——i粒子在n+1时层的压力Pa；
α、γ——静态系数；
由计算得到的压力值，根据公式(36)计算速度的修正值；
Figure BDA0002264623960000141
式中，
ui'——i粒子的速度修正值；
再由计算得到的速度修正值，根据公式(37)，公式(38)分别求解得到粒子在下一时层的速度及位置；
Figure BDA0002264623960000142
Figure BDA0002264623960000143
步骤5：根据具体问题的求解需要，对步骤2～步骤4进行时间步长Δt推进进行计算；
步骤6：输出计算所得到的堆芯中粒子的种类、速度
Figure BDA0002264623960000144
位置
Figure BDA0002264623960000145
压力
Figure BDA0002264623960000146
及温度值Ti n+1，粒子质量
Figure BDA0002264623960000147
及粒子内各元素质量，按照要求进行数据后处理，得到所要求的结果进行输出分析堆芯熔化进程的形态变化，堆芯内部热量分布、堆芯熔化情况、熔融物流动及分布的情况。
本发明方法可研究堆芯内共晶反应及高温熔化行为，揭示反应堆堆芯燃料熔化和熔融物扩展行为，为核电站制定完善的严重事故环节策略提供技术支持。
和现有技术相比，本发明方法具备如下优点：
本发明的一种利用粒子网格混合法研究堆芯内共晶反应及高温熔化行为的方法，在计算中使用有限体积法精确计算了堆芯内材料在共晶反应过程中的质量扩散，及升温过程中的热量传递，使用相变准则判断状态变化，可以获得精确的堆芯内材料状态；再结合使用移动粒子半隐式方法计算堆芯内物质的流动状态，可精确处理自由表面的变化。整个过程不需人工划分网格，大大节省人力，而且相对于纯移动粒子半隐式方法计算堆芯熔化而言，其结果更加精确可信。本发明方法能够对反应堆严重事故中堆芯内共晶反应及高温熔化行为进行准确的研究，为后续事故现象的研究提供准确的条件，从而更全面、有效地评估反应堆的安全性提供依据。
附图说明
图1是本发明利用粒子网格混合法研究堆芯内共晶反应及高温熔化行为的方法流程图；
图2是i粒子周围1.2倍半径内的随机三个粒子j1、j2、j3对应标记为粒子中心点J1、J2、J3的位置示意图；
图3是i粒子与周围随机三个构成平面的粒子j1、j2、j3所构成的四面体IJ1J2J3示意图；
图4是三角形J1J2J3面积划分示意图；
图5是四面体IJ1J2J3体积划分示意图(顶点I控制部分)；
图6是某种类型粒子的A-B二元相图；
具体实施方式
本发明一种利用粒子网格混合法研究堆芯内共晶反应及高温熔化行为的方法，如图1所示，步骤如下：
步骤1：对反应堆堆芯进行粒子建模，不同的物质使用不同种类的粒子表示，粒子初始状态包括粒子的位置、速度、压力和温度，对应标记为ri、ui、pi和Ti，粒子初始布置时相邻两粒子间的距离记为l0，每个粒子内按所含元素标记对应元素的质量，如粒子i内包含A、B两种元素，则其质量对应标记为miA、miB，粒子i的质量mi使用公式(1)计算，
mi＝miA+miB 公式(1)
粒子i的密度ρi使用公式(2)计算，
Figure BDA0002264623960000151
对于熔化焓等则通过公式(3)求得，
Figure BDA0002264623960000152
步骤2：根据各粒子的位置划分网格：在本步骤中，对于每一个粒子i，设置粒子作用域半径lr，若粒子i的坐标位置为(Ix，Iy，Iz)，设定粒子i周围的某一粒子为粒子j，坐标位置为(Jx，Jy，Jz)，粒子i与粒子j的距离为lij，使用公式(4)计算，
Figure BDA0002264623960000161
若满足lij≤lr，则将该就j粒子标记为i粒子的有效作用粒子。将所有满足条件的m个j粒子进行汇总并标记为j1，j2，……，jm，然后对其中每3个粒子进行组合，按照排列组合法，将有
Figure BDA0002264623960000162
种组合，
Figure BDA0002264623960000163
由公式(5)进行计算，
Figure BDA0002264623960000164
对每个组合中的粒子进行考察，取其中一个组合的3个粒子为例，这3个粒子必须满足能且只能构成一个平面的要求，如图2所示，将取出的3个粒子标记为j1，j2，j3，粒子中心分别标记为J1，J2，J3，空间坐标对应为(Jx1，Jy1，Jz1)，(Jx2，Jy2，Jz2)，(Jx3，Jy3，Jz3)，连接J1、J2，连接J1、J3，连接J2、J3，J1、J2、J3 三点构成平面的其中一个充要条件是直线J1J2与直线J1J3之间存在一个不为0的夹角，则只需满足公式(6)，
Figure BDA0002264623960000165
式中，||J1J2||、||J1J3||、||J2J3||分别为线段J1J2、J1J3、J2J3的长度，可以使用公式(7)进行计算，
Figure BDA0002264623960000166
对每一个i粒子，针对其
Figure BDA0002264623960000167
种组合进行判断，若组合内的3个粒子能够满足构成唯一平面的要求则留下该组合，随后针对留下的组合，考察i粒子与组合内3个粒子的关系，i粒子与组合中三个粒子必须能够构成四面体。对于留下的每个组合而言，这3个粒子已经满足构成唯一平面，那么要使这i粒子与组合中的3个粒子能够构成四面体，只需满足i粒子不在组合中的3个粒子所构成的平面上即可，如图3所示，添加i粒子，粒子中心标记为I，空间坐标定义为(Ix，Iy，Iz)，若要使I、J1、J2、J3四个点构成四面体，则只需要满足i粒子到平面J1J2J3的距离为正数，使用公式(8)进行判定，
Figure BDA0002264623960000171
式中：
Figure BDA0002264623960000172
——i粒子到平面J1J2J3的距离m；
A、B、C、D——平面J1J2J3的平面方程的四个系数，通过公式(9)获得，
Figure BDA0002264623960000173
i粒子与j1、j2、j3粒子构成的平面J1J2J3的距离若满足公式(8)，则说明i粒子能够与组合内3个粒子构成的三角形J1J2J3构成四面体IJ1J2J3，其四个面分别为三角形IJ1J2、三角形IJ1J3、三角形IJ2J3和三角形J1J2J3，面积对应标记为
Figure BDA0002264623960000174
三角形的面积可以使用如公式(10)所示的海伦公式计算
Figure BDA0002264623960000175
式中，
S——三角形的面积；
la、lb、lc——三角形的三边长；
p——三角形的半周长，即p＝(la+lb+lc)/2；
四面体IJ1J2J3的体积通过公式(11)进行计算，
Figure BDA0002264623960000176
四面体IJ1J2J3中四个顶点I、J1、J2、J3分别对四面体总体积的分属控制体积，对应标记为VI、
Figure BDA0002264623960000181
这4个体积由公式(12)所示的关系，
Figure BDA0002264623960000182
各控制体积由公式(13)进行计算，
Figure BDA0002264623960000183
除此之外，四面体IJ1J2J3的四个点的分属控制体积与其相对的三角形面积是成反比的，即有公式(14)所示的关系，
Figure BDA0002264623960000184
由以上公式联立求解即可求得四面体IJ1J2J3的四个点的分属控制体积。
然后将对四面体IJ1J2J3中四个面的三角形进行面积划分，如图4所示，以三角形J1J2J3为例，在三角形J1J2J3中，三角形J1J2J3的面积
Figure BDA0002264623960000185
通过公式(10)海伦公式求出；点N12为线段J1J2的中点，点N13为线段J1J3的中点，点N23为线段J2J3的中点。此时需要在三角形中寻找P点，对于三角形J1J2J3而言，将其标记为Px，连接Px、J1，连接Px、J2，连接Px、J3，连接Px、N12，连接Px、N13 ，连接Px、N23 ，可以将三角形J1J2J3划分成6个小三角形，分别标记为三角形J1PxN12，三角形J1PxN13，三角形J2PxN12，三角形J2PxN23，三角形J3PxN13，三角形J3PxN23；因为点N13为线段J1J3的中点，所以三角形J1PxN13与三角形J3PxN13面积相等，标记为
Figure BDA0002264623960000186
同样的，有
Figure BDA0002264623960000187
对于各个三角形的面积，有如公式(15)和公式(16)所示的几何关系，
Figure BDA0002264623960000188
Figure BDA0002264623960000191
三角形的三个角的角度由类似于公式(6)所示的余弦定理求出；由公式(15)和公式(16)求出三角形J1J2J3所分成的6个小三角形的面积；同理也能求得四面体IJ1J2J3另外三个面所分成的共计18个小三角形的面积；
在四面体IJ1J2J3中寻找一点Qx，如图5所示，对于顶点I的控制体积而言，连接QxPx1，QxPx2，QxPx3，QxN1，QxN2，QxN3，QxI(Px1、Px2、Px3、分别对应为三角形IJ1J2、三角形IJ2J3、三角形IJ1J3所寻找的P点，N1、N2、N3则分别为线段IJ1、IJ2、IJ3的中点)，至此点I的控制体积分成了六个四面体，分别为四面体IPx1N1Qx，四面体IPx1N2Qx，四面体IPx2N2Qx，四面体IPx2N3Qx，四面体IPx3N1Qx和四面体IPx3N3Qx，其体积对应标记为
Figure BDA0002264623960000192
Figure BDA0002264623960000193
和
Figure BDA0002264623960000194
对于四面体的体积，有公式(17)所示的关系，
Figure BDA0002264623960000195
设点Qx距离平面J1J2J3、平面IJ1J2，平面IJ2J3，平面IJ1J3的距离分别为H0、H1、H2、H3，对于顶点I的控制体积可列出如公式(18)的方程，
Figure BDA0002264623960000196
同理对于顶点J1、J2、J3也列出相似的方程，联立这4个方程及公式(13)能够求出H0、H1、H2、H3的值；
至此，在I、J1、J2、J3组成的四面体中，J1对I起影响的点I的控制体积为
Figure BDA0002264623960000197
记为
Figure BDA0002264623960000198
J2对I起影响的点I的控制体积为
Figure BDA0002264623960000199
记为
Figure BDA00022646239600001910
J3对I起影响的点I的控制体积为
Figure BDA00022646239600001911
记为
Figure BDA00022646239600001912
所以，在I、J1、J2、J3组成的四面体系统中，点J1对I的有效传输面积标记为
Figure BDA00022646239600001913
点J2对I的有效传输面积标记为
Figure BDA00022646239600001914
点J1对I的有效传输面积标记为
Figure BDA00022646239600001915
通过公式(19)计算，
Figure BDA0002264623960000201
式中，
Figure BDA0002264623960000202
——I点到J1点的距离；
Figure BDA0002264623960000203
——I点到J2点的距离；
Figure BDA0002264623960000204
——I点到J3点的距离；
以上的划分网格的方法描述是针对其中一个四面体系统而言的，在计算中需要对i粒子的所有符合的四面体网格进行计算得到有效传输面积数据之后才能开始进行质量扩散计算及温度计算；
步骤3：在质量扩散的计算中，使用的传质方程为菲克第二定律，方程如公式(20)所示，对于粒子i内的A、B元素可以分别列出方程，
Figure BDA0002264623960000205
式中，
m——粒子内的某元素质量kg；
t——时间s；
Ddiff——粒子内的某元素的质量扩散系数m2/s；
公式(21)给出质量扩散方程的Δt时间步长的积分形式，
Figure BDA0002264623960000206
式中，
m——粒子内的某元素质量kg；
t——时间s；
Δt——时间步长s；
V——体积m3；
CV——控制体积m3；
Ddiff——粒子内的某元素的质量扩散系数m2/s；
对于粒子i而言，其已经通过步骤2划分为若干个四面体网格，所以公式(21)能够离散成公式(22)的形式，
Figure BDA0002264623960000211
公式(22)等式的左边对i粒子划分成的多个四面体的某一元素的质量变化进行加和，等式的右边则是i粒子的该元素的质量变化率，计算中不考虑粒子的体积变化；式中，
Sj→i——步骤2中网格划分完成后计算的j粒子对i粒子的有效传输面积；
Ddiff——某元素的质量扩散系数m2/s；
rj——j粒子的坐标矢量；
ri——i粒子的坐标矢量；
——j粒子中x元素在n时层的质量kg；
Figure BDA0002264623960000213
——i粒子中x元素在n时层的质量kg；
Figure BDA0002264623960000214
——i粒子中x元素在n+1时层的质量，即在n时层的时间步长Δt后的x元素质量kg；
i粒子在n+1时层的质量
Figure BDA0002264623960000215
即可使用公式(23)计算，
Figure BDA0002264623960000216
在温度计算中，能量守恒方程如公式(24)所示，
Figure BDA0002264623960000217
式中，
T——温度K；
t——时间s
κ——导热系数W/(m2·K)；
ρ——密度kg/m3；
cp——比定压热容J/(kg·K)；
ST——温度的源项；
公式(25)给出能量方程的Δt时间步长的积分形式，
Figure BDA0002264623960000221
式中，
t——时间s；
Δt——时间步长s；
T——温度K；
V——体积m3；
CV——控制体积m3；
κ——导热系数W/(m2·K)；
ρ——密度kg/m3；
cp——比定压热容J/(kg·K)；
ST——温度的源项；
对于粒子i而言，其已经通过步骤2划分为若干个四面体网格，所以公式(25)能够离散成公式(26)的形式，
Figure BDA0002264623960000222
公式(26)等式的左边对i粒子划分成的多个四面体的能量变化进行加和，等式的右边则是i粒子的能量变化率；
式中：
Sj→i——步骤2中网格划分完成后计算的j粒子对i粒子的有效传输面积；
κ——导热系数W/(m2·K)；
rj——j粒子的坐标矢量；
ri——i粒子的坐标矢量；
Figure BDA0002264623960000231
——j粒子在n时层的温度K；
Ti n——i粒子在n时层的温度K；
ρi——i粒子的密度kg/m3；
cpi——i粒子的比定压热容J/(kg·K)；
Ti n+1——i粒子在n+1时层的温度，即在n时层的时间步长Δt后的温度K；
求解公式(26)即得到i粒子在时间步长Δt之后的温度Ti n+1；
至此，i粒子在时间步长Δt之后的温度Ti n+1及其质量
Figure BDA0002264623960000232
已经求解得到，并且对于i粒子内的各元素进行的单独质量求解极易得到各元素在粒子内的质量百分比，以粒子内只含有A、B两种元素为例，i粒子中A元素的质量百分比可由公式(27)计算，
Figure BDA0002264623960000233
式中，
CiA——i粒子中A元素的质量百分比；
miA——i粒子中A元素的质量kg；
miB——i粒子中B元素的质量kg；
此时使用A、B元素的二元相图进行粒子状态判定，如图6所示，在A-B二元相图中，横坐标为A元素的质量百分比，纵坐标为粒子的温度，曲线ApQSBp为固相线，曲线ApQLBp为液相线，粒子在多边形O1O2BpQSAp围成的区域内状态为固态，在多边形ApQLBpO3O4围成的区域内状态为液态，在多边形ApQSBpQL围成的区域内状态为固液两相混合状态。若在某一时刻，i粒子的温度为Tk，i粒子中A元素的质量百分比为Ck，在A-B二元相图中作温度Tk的等温线，做质量百分比Ck的等值线，交于点Qp，此时可以看到点Qp位于A-B二元相图的液相区域，则对应的，粒子相态为液相。对位于固液两相混合状态区域内的粒子，需要根据相图中常用到的杠杆定律对其成分进行分析。同样的方法，对于温度为Tpx和粒子内A元素质量百分比为Cpx的粒子k，可判断其对应的Qpx点位于固液两相区内，作线段TpxQpx的延长线交固相线和液相线于QS和QL，再过QS和QL作垂直线交横坐标于CS和CL点，则粒子k中的固体和液体质量之间的关系可使用公式(28)计算得到，
Figure BDA0002264623960000241
式中，
mS——k粒子中处于固态的质量kg；
mL——k粒子中处于液态的质量kg；
Cpx——k粒子中A元素的质量百分比；
CS——点CS对应的横坐标数值；
CL——点CL对应的横坐标数值；
步骤4：使用移动粒子半隐式方法计算动量方程；动量方程如公式(29)所示
Figure BDA0002264623960000242
式中，
u——速度m/s；
t——时间s；
ρ——密度kg/m3；
p——压力Pa
ν——运动粘度m2/s；
f——外力加速度m2/s；
在此方程组中，动量方程中的粘性项及源项显式计算，其中粘性项的Laplace算子由公式(30)离散求解，
Figure BDA0002264623960000251
式中，
d——计算的维度；
n0——初始粒子数密度，通过
Figure BDA0002264623960000252
计算，
Figure BDA0002264623960000253
是核函数，使用公式(31)计算，只需将括号内的自变量修改为
Figure BDA0002264623960000254
即可，
Figure BDA0002264623960000255
为j粒子初始时刻的坐标矢量，r0为i粒子初始时刻的坐标矢量；
i
λ——是i粒子的邻居粒子和i粒子间距的方程，使用公式(32)计算，
Figure BDA0002264623960000256
式中，
re——粒子作用半径，位于半径内的粒子才会与i粒子发生相互作用，，re＝3.1l0；
Figure BDA0002264623960000257
速度和位置的估算值
Figure BDA0002264623960000258
由公式(33)进行计算，
Figure BDA0002264623960000259
式中，
Figure BDA00022646239600002510
——i粒子在上一时间步的速度矢量m/s；
ri n——i粒子在上一时间步的位置矢量m；
Ti n——i粒子在上一时间步的温度K；
根据估算得到的粒子位置，通过公式(34)计算此时的粒子数密度n* i，只需要将公式中的|rj-ri|使用位置的估算值
Figure BDA0002264623960000261
替代即可，由公式(34)所示的压力Poisson方程求解得到此时各粒子的压力值pi；
Figure BDA0002264623960000262
Figure BDA0002264623960000263
公式(35)等式的左边使用如同公式(30)的形式离散，式中，
Figure BDA0002264623960000264
——i粒子在n+1时层的压力Pa；
α、γ——静态系数；
由计算得到的压力值，根据公式(36)计算速度的修正值；
Figure BDA0002264623960000265
式中，
u'i——i粒子的速度修正值；
再由计算得到的速度修正值，根据公式(37)，公式(38)分别求解得到粒子在下一时层的速度及位置；
Figure BDA0002264623960000266
Figure BDA0002264623960000267
步骤5：根据具体问题的求解需要，对步骤2～步骤4进行时间步长Δt推进进行计算；
步骤6：输出计算所得到的堆芯中粒子的种类、速度
Figure BDA0002264623960000268
位置
Figure BDA0002264623960000269
压力
Figure BDA00022646239600002610
及温度值Ti n+1，粒子质量
Figure BDA00022646239600002611
及粒子内各元素质量，按照要求进行数据后处理，得到所要求的结果进行输出分析堆芯熔化进程的形态变化，堆芯内部热量分布、堆芯熔化情况、熔融物流动及分布等的情况。本发明方法能够对反应堆严重事故中堆芯内共晶反应及高温熔化行为进行准确的研究，为后续事故现象的研究提供准确的条件，从而更全面、有效地评估反应堆的安全性提供依据。

Claims (1)

1.一种利用粒子网格混合法研究堆芯内共晶反应及高温熔化行为的方法，其特征在于：步骤如下：
步骤1：对反应堆堆芯进行粒子建模，不同的物质使用不同种类的粒子表示，粒子初始状态包括粒子的位置、速度、压力和温度，对应标记为ri、ui、pi和Ti，粒子初始布置时相邻两粒子间的距离记为l0，每个粒子内按所含元素标记对应元素的质量，如粒子i内包含A、B两种元素，则其质量对应标记为miA、miB，粒子i的质量mi使用公式(1)计算，
mi＝miA+miB 公式(1)
粒子i的密度ρi使用公式(2)计算，
Figure FDA0002264623950000011
对于熔化焓等则通过公式(3)求得，
Figure FDA0002264623950000012
步骤2：根据各粒子的位置划分网格：在本步骤中，对于每一个粒子i，设置粒子作用域半径lr，若粒子i的坐标位置为(Ix，Iy，Iz)，设定粒子i周围的某一粒子为粒子j，坐标位置为(Jx，Jy，Jz)，粒子i与粒子j的距离为lij，使用公式(4)计算，
Figure FDA0002264623950000013
若满足lij≤lr，则将该就j粒子标记为i粒子的有效作用粒子；将所有满足条件的m个j粒子进行汇总并标记为j1，j2，……，jm，然后对其中每3个粒子进行组合，按照排列组合法，将有
Figure FDA0002264623950000014
种组合，
Figure FDA0002264623950000015
由公式(5)进行计算，
Figure FDA0002264623950000021
对每个组合中的粒子进行考察，取其中一个组合的3个粒子，这3个粒子必须满足能且只能构成一个平面的要求，将取出的3个粒子标记为j1，j2，j3，粒子中心分别标记为J1，J2，J3，空间坐标对应为(Jx1，Jy1，Jz1)，(Jx2，Jy2，Jz2)，(Jx3，Jy3，Jz3)，连接J1、J2，连接J1、J3，连接J2、J3，J1、J2、J3三点构成平面的其中一个充要条件是直线J1J2与直线J1J3之间存在一个不为0的夹角，则只需满足公式(6)，
Figure FDA0002264623950000022
式中，||J1J2||、||J1J3||、||J2J3||分别为线段J1J2、J1J3、J2J3的长度，使用公式(7)进行计算，
Figure FDA0002264623950000023
对每一个i粒子，针对其
Figure FDA0002264623950000024
种组合进行判断，若组合内的3个粒子能够满足构成唯一平面的要求则留下该组合，随后针对留下的组合，考察i粒子与组合内3个粒子的关系，i粒子与组合中三个粒子必须能够构成四面体；对于留下的每个组合而言，这3个粒子已经满足构成唯一平面，那么要使这i粒子与组合中的3个粒子能够构成四面体，只需满足i粒子不在组合中的3个粒子所构成的平面上即可，添加i粒子，粒子中心标记为I，空间坐标定义为(Ix，Iy，Iz)，若要使I、J1、J2、J3四个点构成四面体，则只需要满足i粒子到平面J1J2J3的距离为正数，使用公式(8)进行判定，
Figure FDA0002264623950000031
式中：
Figure FDA0002264623950000032
——i粒子到平面J1J2J3的距离m；
A、B、C、D——平面J1J2J3的平面方程的四个系数，通过公式(9)获得，
Figure FDA0002264623950000033
i粒子与j1、j2、j3粒子构成的平面J1J2J3的距离若满足公式(8)，则说明i粒子能够与组合内3个粒子构成的三角形J1J2J3构成四面体IJ1J2J3，其四个面分别为三角形IJ1J2、三角形IJ1J3、三角形IJ2J3和三角形J1J2J3，面积对应标记为
Figure FDA0002264623950000034
Figure FDA0002264623950000035
三角形的面积使用如公式(10)所示的海伦公式计算
Figure FDA0002264623950000036
式中，
S——三角形的面积；
la、lb、lc——三角形的三边长；
p——三角形的半周长，即p＝(la+lb+lc)/2；
四面体IJ1J2J3的体积通过公式(11)进行计算，
Figure FDA0002264623950000037
四面体IJ1J2J3中四个顶点I、J1、J2、J3分别对四面体总体积的分属控制体积，对应标记为VI、
Figure FDA0002264623950000038
这4个体积有公式(12)所示的关系，
Figure FDA0002264623950000039
各控制体积由公式(13)进行计算，
Figure FDA0002264623950000041
除此之外，四面体IJ1J2J3的四个点的分属控制体积与其相对的三角形面积是成反比的，即有公式(14)所示的关系，
Figure FDA0002264623950000042
由以上公式联立求解即求得四面体IJ1J2J3的四个点的分属控制体积；
然后将对四面体IJ1J2J3中四个面的三角形进行面积划分，以三角形J1J2J3为例，在三角形J1J2J3中，三角形J1J2J3的面积
Figure FDA0002264623950000043
通过公式(10)海伦公式求出；点N12为线段J1J2的中点，点N13为线段J1J3的中点，点N23为线段J2J3的中点；此时需要在三角形中寻找P点，对于三角形J1J2J3而言，将其标记为Px，连接Px、J1，连接Px、J2，连接Px、J3，连接Px、N12，连接Px、N13，连接Px、N23，可以将三角形J1J2J3划分成6个小三角形，分别标记为三角形J1PxN12，三角形J1PxN13，三角形J2PxN12，三角形J2PxN23，三角形J3PxN13，三角形J3PxN23；因为点N13为线段J1J3的中点，所以三角形J1PxN13与三角形J3PxN13面积相等，标记为
Figure FDA0002264623950000044
同样的，有
Figure FDA0002264623950000045
对于各个三角形的面积，有如公式(15)和公式(16)所示的几何关系，
Figure FDA0002264623950000051
Figure FDA0002264623950000052
三角形的三个角的角度由类似于公式(6)所示的余弦定理求出；由公式(15)和公式(16)求出三角形J1J2J3所分成的6个小三角形的面积；同理也能求得四面体IJ1J2J3另外三个面所分成的共计18个小三角形的面积；
在四面体IJ1J2J3中寻找一点Qx，对于顶点I的控制体积而言，连接QxPx1，QxPx2，QxPx3，QxN1，QxN2，QxN3，QxI，其中Px1、Px2、Px3、分别对应为三角形IJ1J2、三角形IJ2J3、三角形IJ1J3所寻找的P点，N1、N2、N3则分别为线段IJ1、IJ2、IJ3的中点，至此点I的控制体积分成了六个四面体，分别为四面体IPx1N1Qx，四面体IPx1N2Qx，四面体IPx2N2Qx，四面体IPx2N3Qx，四面体IPx3N1Qx和四面体IPx3N3Qx，其体积对应标记为
Figure FDA0002264623950000053
和
Figure FDA0002264623950000054
对于四面体的体积，有公式(17)所示的关系，
Figure FDA0002264623950000055
设点Qx距离平面J1J2J3、平面IJ1J2，平面IJ2J3，平面IJ1J3的距离分别为H0、H1、H2、H3，对于顶点I的控制体积可列出如公式(18)的方程，
Figure FDA0002264623950000056
同理对于顶点J1、J2、J3也列出相似的方程，联立这4个方程及公式(13)能够求出H0、H1、H2、H3的值；
至此，在I、J1、J2、J3组成的四面体中，J1对I起影响的点I的控制体积为
Figure FDA0002264623950000057
记为
Figure FDA0002264623950000058
J2对I起影响的点I的控制体积为
Figure FDA0002264623950000059
记为
Figure FDA00022646239500000510
J3对I起影响的点I的控制体积为
Figure FDA00022646239500000511
记为
Figure FDA00022646239500000512
所以，在I、J1、J2、J3组成的四面体系统中，点J1对I的有效传输面积标记为
Figure FDA0002264623950000061
点J2对I的有效传输面积标记为
Figure FDA0002264623950000062
点J1对I的有效传输面积标记为
Figure FDA0002264623950000063
通过公式(19)计算，
Figure FDA0002264623950000064
式中，
Figure FDA0002264623950000065
——I点到J1点的距离；
Figure FDA0002264623950000066
——I点到J2点的距离；
Figure FDA0002264623950000067
——I点到J3点的距离；
以上的划分网格的方法描述是针对其中一个四面体系统而言的，在计算中需要对i粒子的所有符合的四面体网格进行计算得到有效传输面积数据之后才能开始进行质量扩散计算及温度计算；
步骤3：在质量扩散的计算中，使用的传质方程为菲克第二定律，方程如公式(20)所示，对于粒子i内的A、B元素分别列出方程，
Figure FDA0002264623950000068
式中，
m——粒子内的某元素质量kg；
t——时间s；
Ddiff——粒子内的某元素的质量扩散系数m2/s；
公式(21)给出质量扩散方程的Δt时间步长的积分形式，
Figure FDA0002264623950000069
式中，
m——粒子内的某元素质量kg；
t——时间s；
Δt——时间步长s；
V——体积m3；
CV——控制体积m3；
Ddiff——粒子内的某元素的质量扩散系数m2/s；
对于粒子i而言，其已经通过步骤2划分为若干个四面体网格，所以公式(21)能够离散成公式(22)的形式，
Figure FDA0002264623950000071
公式(22)等式的左边对i粒子划分成的多个四面体的某一元素的质量变化进行加和，等式的右边则是i粒子的该元素的质量变化率，计算中不考虑粒子的体积变化；式中，
Sj→i——步骤2中网格划分完成后计算的j粒子对i粒子的有效传输面积；
Ddiff——某元素的质量扩散系数m2/s；
rj——j粒子的坐标矢量；
ri——i粒子的坐标矢量；
Figure FDA0002264623950000072
——j粒子中x元素在n时层的质量kg；
Figure FDA0002264623950000073
——i粒子中x元素在n时层的质量kg；
Figure FDA0002264623950000074
——i粒子中x元素在n+1时层的质量，即在n时层的时间步长Δt后的x元素质量kg；
i粒子在n+1时层的质量
Figure FDA0002264623950000075
使用公式(23)计算，
Figure FDA0002264623950000076
在温度计算中，能量守恒方程如公式(24)所示，
Figure FDA0002264623950000081
式中，
T——温度K；
t——时间s
κ——导热系数W/(m2·K)；
ρ——密度kg/m3；
cp——比定压热容J/(kg·K)；
ST——温度的源项；
公式(25)给出能量方程的Δt时间步长的积分形式，
Figure FDA0002264623950000082
式中，
t——时间s；
Δt——时间步长s；
T——温度K；
V——体积m3；
CV——控制体积m3；
κ——导热系数W/(m2·K)；
ρ——密度kg/m3；
cp——比定压热容J/(kg·K)；
ST——温度的源项；
对于粒子i而言，其已经通过步骤2划分为若干个四面体网格，所以公式(25)能够离散成公式(26)的形式，
Figure FDA0002264623950000091
公式(26)等式的左边对i粒子划分成的多个四面体的能量变化进行加和，等式的右边则是i粒子的能量变化率；式中：
Sj→i——步骤2中网格划分完成后计算的j粒子对i粒子的有效传输面积；
κ——导热系数W/(m2·K)；
rj——j粒子的坐标矢量；
ri——i粒子的坐标矢量；
Figure FDA0002264623950000092
——j粒子在n时层的温度K；
Ti n——i粒子在n时层的温度K；
ρi——i粒子的密度kg/m3；
cpi——i粒子的比定压热容J/(kg·K)；
Ti n+1——i粒子在n+1时层的温度，即在n时层的时间步长Δt后的温度K；
求解公式(26)即得到i粒子在时间步长Δt之后的温度Ti n+1；
至此，i粒子在时间步长Δt之后的温度Ti n+1及其质量
Figure FDA0002264623950000093
已经求解得到，并且对于i粒子内的各元素进行的单独质量求解极易得到各元素在粒子内的质量百分比，若粒子内只含有A、B两种元素，i粒子中A元素的质量百分比由公式(27)计算，
Figure FDA0002264623950000094
式中，
CiA——i粒子中A元素的质量百分比；
miA——i粒子中A元素的质量kg；
miB——i粒子中B元素的质量kg；
此时使用A、B元素的二元相图进行粒子状态判定，在A-B二元相图中，横坐标为A元素的质量百分比，纵坐标为粒子的温度，曲线ApQSBp为固相线，曲线ApQLBp为液相线，粒子在多边形O1O2BpQSAp围成的区域内状态为固态，在多边形ApQLBpO3O4围成的区域内状态为液态，在多边形ApQSBpQL围成的区域内状态为固液两相混合状态；若在某一时刻，i粒子的温度为Tk，i粒子中A元素的质量百分比为Ck，在A-B二元相图中作温度Tk的等温线，做质量百分比Ck的等值线，交于点Qp，此时点Qp位于A-B二元相图的液相区域，则对应的，粒子相态为液相；对位于固液两相混合状态区域内的粒子，需要根据相图中的杠杆定律对其成分进行分析；同样的方法，对于温度为Tpx和粒子内A元素质量百分比为Cpx的粒子k，判断其对应的Qpx点位于固液两相区内，作线段TpxQpx的延长线交固相线和液相线于QS和QL，再过QS和QL作垂直线交横坐标于CS和CL点，则粒子k中的固体和液体质量之间的关系使用公式(28)计算得到，
Figure FDA0002264623950000101
式中，
mS——k粒子中处于固态的质量kg；
mL——k粒子中处于液态的质量kg；
Cpx——k粒子中A元素的质量百分比；
CS——点CS对应的横坐标数值；
CL——点CL对应的横坐标数值；
步骤4：使用移动粒子半隐式方法计算动量方程；动量方程如公式(29)所示
Figure FDA0002264623950000111
式中，
u——速度m/s；
t——时间s；
ρ——密度kg/m3；
p——压力Pa
ν——运动粘度m2/s；
f——外力加速度m2/s；
在此方程组中，动量方程中的粘性项及源项显式计算，其中粘性项的Laplace算子由公式(30)离散求解，
Figure FDA0002264623950000112
式中，
d——计算的维度；
n0——初始粒子数密度，通过
Figure FDA0002264623950000113
计算，
Figure FDA0002264623950000114
是核函数，使用公式(31)计算，只需将括号内的自变量修改为
Figure FDA0002264623950000115
即可，
Figure FDA0002264623950000116
为j粒子初始时刻的坐标矢量，
Figure FDA0002264623950000117
为i粒子初始时刻的坐标矢量；
λ——是i粒子的邻居粒子和i粒子间距的方程，使用公式(32)计算，
Figure FDA0002264623950000118
式中，
re——粒子作用半径，位于半径内的粒子才会与i粒子发生相互作用，，re＝3.1l0；
Figure FDA0002264623950000121
速度和位置的估算值
Figure FDA0002264623950000122
由公式(33)进行计算，
Figure FDA0002264623950000123
式中，
Figure FDA0002264623950000124
——i粒子在上一时间步的速度矢量m/s；
ri n——i粒子在上一时间步的位置矢量m；
Ti n——i粒子在上一时间步的温度K；
根据估算得到的粒子位置，通过公式(34)计算此时的粒子数密度n* i，只需要将公式中的|rj-ri|使用位置的估算值
Figure FDA0002264623950000125
替代即可，由公式(34)所示的压力Poisson方程求解得到此时各粒子的压力值pi；
Figure FDA0002264623950000126
Figure FDA0002264623950000127
公式(35)等式的左边使用如同公式(30)的形式离散，式中，
Figure FDA0002264623950000128
——i粒子在n+1时层的压力Pa；
α、γ——静态系数；
由计算得到的压力值，根据公式(36)计算速度的修正值；
Figure FDA0002264623950000129
式中，
u'i——i粒子的速度修正值；
再由计算得到的速度修正值，根据公式(37)，公式(38)分别求解得到粒子在下一时层的速度及位置；
Figure FDA0002264623950000131
ri n+1＝ri *+u'iΔt 公式(38)
步骤5：根据具体问题的求解需要，对步骤2～步骤4进行时间步长Δt推进进行计算；
步骤6：输出计算所得到的堆芯中粒子的种类、速度
Figure FDA0002264623950000132
位置
Figure FDA0002264623950000133
压力
Figure FDA0002264623950000134
及温度值Ti n+1，粒子质量
Figure FDA0002264623950000135
及粒子内各元素质量，按照要求进行数据后处理，得到所要求的结果进行输出分析堆芯熔化进程的形态变化，堆芯内部热量分布、堆芯熔化情况、熔融物流动及分布的情况。
CN201911083355.0A 2019-11-07 2019-11-07 利用粒子网格混合法研究堆芯内共晶反应及高温熔化行为的方法 Pending CN110867220A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201911083355.0A CN110867220A (zh) 2019-11-07 2019-11-07 利用粒子网格混合法研究堆芯内共晶反应及高温熔化行为的方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201911083355.0A CN110867220A (zh) 2019-11-07 2019-11-07 利用粒子网格混合法研究堆芯内共晶反应及高温熔化行为的方法

Publications (1)

Publication Number Publication Date
CN110867220A true CN110867220A (zh) 2020-03-06

Family

ID=69653920

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201911083355.0A Pending CN110867220A (zh) 2019-11-07 2019-11-07 利用粒子网格混合法研究堆芯内共晶反应及高温熔化行为的方法

Country Status (1)

Country Link
CN (1) CN110867220A (zh)

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111832214A (zh) * 2020-06-29 2020-10-27 西安交通大学 基于粒子法的核反应堆严重事故碎片床熔化过程模拟方法
CN112102894A (zh) * 2020-09-04 2020-12-18 西安交通大学 基于粒子法的核反应堆堆芯材料熔池演变特性分析方法
CN113192567A (zh) * 2021-04-30 2021-07-30 西安交通大学 一种核反应堆板状燃料熔化流固耦合无网格分析方法
CN113435004A (zh) * 2021-05-25 2021-09-24 上海交通大学 堆芯熔融进程中熔融材料迁移质量的计算方法和装置
CN114757122A (zh) * 2022-04-19 2022-07-15 西安交通大学 建立钠冷快堆堆芯解体事故精细热工水力计算模型的方法

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102110188A (zh) * 2009-12-25 2011-06-29 鞍钢股份有限公司 一种连铸铸流温度及固相率分布计算方法
CN105081295A (zh) * 2014-05-21 2015-11-25 宝钢特钢有限公司 一种防止钢锭出现针孔缺陷的低碳结构钢冶炼方法
CN105808503A (zh) * 2016-03-07 2016-07-27 西安交通大学 针对反应堆逐棒计算中解析求解栅元不连续因子的方法
JP2018048352A (ja) * 2016-09-20 2018-03-29 新日鐵住金株式会社 高炉内の溶銑の流速推定方法および高炉の操業方法
US20180190395A1 (en) * 2016-12-30 2018-07-05 Nuscale Power, Llc Nuclear Reactor Protection Systems and Methods
CN110044959A (zh) * 2019-05-13 2019-07-23 西安交通大学 利用移动粒子有限容积法研究高瑞利数熔融池换热特性的方法

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102110188A (zh) * 2009-12-25 2011-06-29 鞍钢股份有限公司 一种连铸铸流温度及固相率分布计算方法
CN105081295A (zh) * 2014-05-21 2015-11-25 宝钢特钢有限公司 一种防止钢锭出现针孔缺陷的低碳结构钢冶炼方法
CN105808503A (zh) * 2016-03-07 2016-07-27 西安交通大学 针对反应堆逐棒计算中解析求解栅元不连续因子的方法
JP2018048352A (ja) * 2016-09-20 2018-03-29 新日鐵住金株式会社 高炉内の溶銑の流速推定方法および高炉の操業方法
US20180190395A1 (en) * 2016-12-30 2018-07-05 Nuscale Power, Llc Nuclear Reactor Protection Systems and Methods
CN110044959A (zh) * 2019-05-13 2019-07-23 西安交通大学 利用移动粒子有限容积法研究高瑞利数熔融池换热特性的方法

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
YONGLIN LI 等: "Numerical analysis of the dissolution of uranium dioxide by molten zircaloy using MPS method", 《PROGRESS IN NUCLEAR ENERGY》 *
徐长英、杨作文: "《电工电子实践基础教程》", 30 September 2013 *
李勇霖 等: "基于MPS方法的共晶反应研究", 《第十五届全国反应堆热工流体学术会议暨中核核反应堆热工水力技术重点实验室学术年会论文集》 *
谢文玲 等: "杠杆法在相图分析中的应用", 《广州化工》 *

Cited By (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111832214A (zh) * 2020-06-29 2020-10-27 西安交通大学 基于粒子法的核反应堆严重事故碎片床熔化过程模拟方法
CN112102894A (zh) * 2020-09-04 2020-12-18 西安交通大学 基于粒子法的核反应堆堆芯材料熔池演变特性分析方法
CN112102894B (zh) * 2020-09-04 2021-11-30 西安交通大学 基于粒子法的核反应堆堆芯材料熔池演变特性分析方法
CN113192567A (zh) * 2021-04-30 2021-07-30 西安交通大学 一种核反应堆板状燃料熔化流固耦合无网格分析方法
CN113192567B (zh) * 2021-04-30 2022-12-09 西安交通大学 一种核反应堆板状燃料熔化流固耦合无网格分析方法
CN113435004A (zh) * 2021-05-25 2021-09-24 上海交通大学 堆芯熔融进程中熔融材料迁移质量的计算方法和装置
CN114757122A (zh) * 2022-04-19 2022-07-15 西安交通大学 建立钠冷快堆堆芯解体事故精细热工水力计算模型的方法
CN114757122B (zh) * 2022-04-19 2024-02-23 西安交通大学 建立钠冷快堆堆芯解体事故精细热工水力计算模型的方法

Similar Documents

Publication Publication Date Title
CN110867220A (zh) 利用粒子网格混合法研究堆芯内共晶反应及高温熔化行为的方法
CN110044959B (zh) 利用移动粒子有限容积法研究熔融池换热特性的方法
CN110472355B (zh) 一种基于多场耦合建模与仿真求解的3d打印预览方法
Gartling et al. Simulation of coupled viscous and porous flow problems
Fan et al. Numerical modeling of the additive manufacturing (AM) processes of titanium alloy
Cordonnier et al. Benchmarking lava-flow models
Tran et al. The effective convectivity model for simulation of melt pool heat transfer in a light water reactor pressure vessel lower head. Part I: Physical processes, modeling and model implementation
Shams et al. Large eddy simulation of a randomly stacked nuclear pebble bed
Xu et al. Physics-informed neural networks for studying heat transfer in porous media
JP2009193110A (ja) グリッドフリー手法を用いた固気二相流シミュレーションプログラム及びそれを記憶した記憶媒体並びに固気二相流シミュレーション装置
Natsume et al. Evaluation of permeability for columnar dendritic structures by three-dimensional numerical flow analysis
TKADLEČKOVÁ et al. Testing of numerical model settings for simulation of steel ingot casting and solidification
Volkov et al. A coolant flow simulation in fast reactor wire-wrapped assembly
Kawakami et al. Improvement of solidification model and analysis of 3D channel blockage with MPS method
Wang et al. Analysis of focusing effect of light metallic layer in stratified molten pool under IVR-ERVC condition
Bakker Lecture 1-Introduction to CFD Applied Computational Fluid Dynamics
Koshizuka et al. R&D of the Next Generation Safety Analysis Methods for Fast Reactors With New Computational Science and Technology: 1—Introduction of the Project and Development of Structural Mechanics Module
Schiano et al. Modeling spreading with degassing using anisotherm viscoplastic multiphase shallow water approximation
Heidler et al. Eulerian–Lagrangian and Eulerian–Eulerian approaches for the simulation of particle-laden free surface flows using the lattice Boltzmann method
Zimbrod et al. Efficient Simulation of Complex Capillary Effects in Advanced Manufacturing Processes using the Finite Volume Method
Mistrangelo et al. Magnetoconvection in HCLL blankets
Li et al. Preliminary evaluation on the relocation phase of ex-vessel debris of fukushima daiichi nuclear power plant unit-3
Vispute et al. Utilizing flow simulation in the design phase of a casting die to optimize design parameters and defect analysis
Ohira et al. Numerical simulations of upper plenum thermal-hydraulics of Monju reactor vessel using high resolution mesh models
Kucala et al. A Control Volume Finite Element Approach for Modeling Molten Corium Spreading and Solidification.

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
WD01 Invention patent application deemed withdrawn after publication
WD01 Invention patent application deemed withdrawn after publication

Application publication date: 20200306