CN110044959A - 利用移动粒子有限容积法研究高瑞利数熔融池换热特性的方法 - Google Patents

利用移动粒子有限容积法研究高瑞利数熔融池换热特性的方法 Download PDF

Info

Publication number
CN110044959A
CN110044959A CN201910394189.XA CN201910394189A CN110044959A CN 110044959 A CN110044959 A CN 110044959A CN 201910394189 A CN201910394189 A CN 201910394189A CN 110044959 A CN110044959 A CN 110044959A
Authority
CN
China
Prior art keywords
particle
formula
triangle
tetrahedron
point
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN201910394189.XA
Other languages
English (en)
Other versions
CN110044959B (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 CN201910394189.XA priority Critical patent/CN110044959B/zh
Publication of CN110044959A publication Critical patent/CN110044959A/zh
Application granted granted Critical
Publication of CN110044959B publication Critical patent/CN110044959B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N25/00—Investigating or analyzing materials by the use of thermal means
    • G01N25/20—Investigating or analyzing materials by the use of thermal means by investigating the development of heat, i.e. calorimetry, e.g. by measuring specific heat, by measuring thermal conductivity

Abstract

一种利用移动粒子有限容积法研究高瑞利数熔融池换热特性的方法，主要步骤如下：1、针对具体问题进行粒子建模，设定粒子初始布置；2、根据各粒子位置划分网格；3、在时间步长Δt之后，使用有限容积法求解能量方程；4、粒子法处理显式计算动量方程中的粘性项和源项，估算粒子速度及位置；5、计算粒子数密度，求解压力Poisson方程；6、计算速度的修正值，并计算各粒子在下一时层的速度及位置；7、输出计算结果，并根据求解的具体问题输出壁面热流密度和熔融池内物质的相态。本方法具有粒子法精确捕捉自由液面及模拟相变问题的优点，又具有较高的计算精度，能够通过此方法来研究高瑞利数熔融池换热特性，为核电厂反应堆严重事故安全特性研究提供重要依据。

Description

利用移动粒子有限容积法研究高瑞利数熔融池换热特性的 方法
技术领域
本发明涉及核电厂严重事故下封头熔融池换热特性研究技术领域，具体涉及一种利用移动粒子有限容积法研究高瑞利数熔融池换热特性的方法。
背景技术
当核电厂压水堆发生严重事故时，堆芯如果不能得到有效冷却，会快速升温导致大量燃料组件和结构材料发生熔化并迁移至压力容器底部，进入到下封头内的高温熔融物逐渐积累形成碎片床或熔融池，不断向下封头壁面传递衰变热，当衰变热不能充分导出时，熔融池会将下封头壁面加热到很高温度，严重威胁压力容器的完整性。反应堆压力容器下封头内熔融池的流动还热特性对下封头壁面热负荷分布以及外部冷却能力有着重要的影响，熔融物堆内持留技术是缓解堆芯熔化事故后果的一项关键措施，运用不同的管理策略，可将熔融物滞留在压力容器下封头内，壁面后续事故的发生和大量放射性物质释放到外界环境中，从而可以有效终止反应堆严重事故。
关于熔融池换热特性的研究，国内外已经开展了一些研究，包括实验研究和数值模拟研究。对于实验研究，由于不同的研究侧重点和实验条件的限制，各实验采用不同比例的实验装置、不同的熔融物模拟物、不同的加热方式以及不同的边界条件，因此也具有不同的瑞利数范围和实验结果。瑞利数作为对流换热问题的关键参数，对熔融池内的对流换热特性具有重要的影响。对比不同实验得到的熔融池换热关系式可以发现，在瑞利数相对较小的范围内不同实验获得的熔融池换热特性参数差别较小。但是在高瑞利数范围内，实验结果的差别开始增大，特别是在反应堆量级下(1016～1017)，不确定性更大。所以国内外也开展了一些大型实验，但相应的实验花费也更大。除了实验研究，国内外还开展了许多理论研究及数值模拟研究，以求在低消耗低成本的基础上对高瑞利数熔融池换热特性进行比较深入的研究。目前使用的方法大多限于商业软件如Fluent等先对下封头熔融池进行建模，然后进行网格划分，其中使用较多的方法是有限体积法，因其具有很好的守恒性，对网格的适应性很好，可以很好的解决复杂的工程问题。但对于反应堆下封头熔融池换热特性而言，其内可能存在氧化层和金属层，还有因为冷却变化和再次受衰变热加热而导致的硬壳生成及再熔化的相变问题，网格法在处理这些问题上不及粒子法的优势大，而且粒子法建模简单，而且能够处理复杂形状的流体流动问题，在处理湍流问题上具有更好的鲁棒性。尽管粒子法在处理流体流动及存在相变的问题上的上有其独特的优势，但是粒子法在计算传热问题时，在计算域的边界因为支持域被截断而存在较大的误差。综上所述，网格法中的有限容积法在计算传热问题上精确度较高，粒子法在求解复杂流体流动问题及相变问题上优势较大。本发明综合粒子法和有限容积法提出了一种研究高瑞利数熔融池传热特性的移动粒子有限容积法。
发明内容
为了克服网格法与粒子法的缺陷，本发明的目的是为了精确并高效地对下封头熔融池换热特性进行研究，提供一种利用移动粒子有限容积法研究高瑞利数熔融池换热特性的方法，该方法能够对存在多种流体、存在相变问题、湍流流动的熔融池换热进行研究，获得熔融池内的温度场分布和壁面热流密度分布情况、熔融池内的硬壳分布特性，为核电厂反应堆严重事故安全特性研究提供重要依据。
为了实现上述目的，本发明采取了以下的技术方案予以实施：
一种利用移动粒子有限容积法研究高瑞利数熔融池换热特性的方法，步骤如下：
步骤1：对高瑞利数熔融池进行粒子建模，不同种类的粒子代表不同的物质，每种粒子根据所表示的物质具有其对应的质量、密度和熔化焓，使用粒子i来表示某一粒子，则其质量、密度、熔化焓对应标记分别为mi、ρi、设定粒子的初始布置，包括粒子的位置、速度、压力和温度，对应标记分别为ri、ui、pi和Ti，粒子初始布置时相邻两粒子间的距离记为l0；
步骤2：根据各粒子的位置划分网格：在本步骤中，设置一个粒子搜索半径，对于每一个粒子i，其坐标位置为(xi，yi，zi)，与其间距lij≤1.2l0的所有粒子作为搜索的对象，式中lij为j粒子与i粒子的距离，j粒子的坐标标记为(xj，yj，zj)，lij使用公式(1)进行计算，
并将搜索到的粒子使用j1，j2，……，jn符号进行标记，其中下标n表示搜索到的粒子总数；将j1，j2，……，jn粒子中的每三个粒子进行组合，将有种组合，由公式(2)进行计算，
对每个组合中的粒子进行考察，在这些组合中，每个组合内的三个坐标点必须满足构成面的要求，取其中的一个组合进行解析，其中一种组合粒子标记为j1，j2，j3，粒子中心分别标记为O1，O2，O3，空间坐标对应为(x1，y1，z1)，(x2，y2，z2)，(x3，y3，z3)，连接O1、O2为直线O1O2，连接O1、O3为直线O1O3，连接O2、O3为直线O2O3，若要判断O1、O2、O3能否构成一个平面，只需要判断直线O1O2与直线O1O3之间存在一个不为0的夹角，则只需满足公式(3)，
对所有的组合进行判断，将能够构成面的组合保留，不能构成面的组合删去，随后考虑i粒子与每个组合中三个粒子的关系，i粒子与组合中三个粒子必须能够构成四面体，因为组合中的三个粒子能够构成平面，所以要使这四个粒子能够构成四面体，只需满足i粒子不在组合中三个粒子所构成的平面上即可，添加i粒子，粒子中心标记为O，空间坐标定义为(xi，yi，zi)，若要使O、O1、O2、O3四个点构成四面体，因为O1、O2、O3三个点已经满足构成面的要求，则只需要满足i粒子到平面O1O2O3的距离不为0，使用公式(4)进行判定，
式中：
——i粒子到平面O1O2O3的距离m；
a、b、c、d——平面O1O2O3的平面方程的四个系数，通过公式(5)获得，
i粒子与j1、j2、j3粒子构成的平面O1O2O3的距离若满足公式(4)，则说明i粒子能够与三角形O1O2O3构成四面体OO1O2O3，四面体OO1O2O3的四个面分别为三角形OO1O2、三角形OO1O3、三角形OO2O3和三角形O1O2O3；使用海伦公式计算四个三角形的面积，四个面的面积对应标记为海伦公式如公式(6)所示，
式中，
S——三角形的面积；
L1、L2、L3——三角形的三边长；
p——三角形的半周长，即p＝(L1+L2+L3)/2；
四面体OO1O2O3的体积通过公式(7)进行计算，
式中，
V——四面体体积，在四面体OO1O2O3中，标记为
S——三角形O1O2O3的面积，标记为使用公式(6)计算；
h——O点与三角形所在平面O1O2O3的距离，记为使用公式(4)计算；
在使用粒子位置划分网格时，还有几个参量需要求解，其一是四面体OO1O2O3中四个顶点O、O1、O2、O3分别对四面体总体积的分属控制体积，对应标记为VO、这个分属控制体积由公式(8)进行计算，
实际上，四面体OO1O2O3的四个点的分属控制体积与其相对的三角形面积成反比，即有公式(9)所示的关系，
其二是对四面体OO1O2O3中四个面的三角形进行面积划分，在三角形O1O2O3中，三角形O1O2O3的面积通过公式(6)求出；点M12为线段O1O2的中点，点M13为线段O1O3的中点，点M23为线段O2O3的中点。此时需要在三角形中寻找一点P，连接P、O1，连接P、O2，连接P、O3，连接P、M12，连接P、M13，连接P、M23，此时三角形O1O2O3被分成了6个小三角形，分别标记为三角形O1PM12，三角形O1PM13，三角形O2PM12，三角形O2PM23，三角形O3PM13，三角形O3PM23；因为点M13为线段O1O3的中点，所以三角形O1PM13与三角形O3PM13面积相等，标记为同样的，三角形O1PM12与三角形O2PM12面积相等，标记为三角形O2PM23与三角形O3PM23面积相等，标记为对于各个三角形的面积，有如公式(10)和公式(11)所示的几何关系，
三角形的三个角的角度由公式(3)所示的余弦定理求出；由公式(10)和公式(11)求出三角形O1O2O3所分成的6个小三角形的面积；同理也能求得四面体OO1O2O3另外三个面所分成的共计18个小三角形的面积；
在四面体OO1O2O3中寻找一点Q，对于顶点O的控制体积而言，连接QP1，QP2，QP3，QM1，QM2，QM3，QO，至此点O的控制体积分成了六个四面体，分别为四面体OP1M1Q，四面体OP1M2Q，四面体OP2M1Q，四面体OP2M3Q，四面体OP3M2Q和四面体OP3M3Q，其体积对应标记为和对于四面体的体积，有公式(12)所示的关系，
设点Q距离平面O1O2O3、平面OO1O2，平面OO1O3，平面OO2O3的距离分别为H、H1、H2、H3，对于顶点O的控制体积列如公式(13)的方程，
同理对于顶点O1、O2、O3也列出相似的方程，联立这4个方程及公式(8)能够求出H、H1、H2、H3的值；
至此，在O、O1、O2、O3组成的四面体系统中，O1对O起影响的点O的控制体积为记为O2对O起影响的点O的控制体积为记为O3对O起影响的点O的控制体积为记为所以，在O、O1、O2、O3组成的四面体系统中，点O1对O的有效传热面积标记为点O2对O的有效传热面积标记为点O1对O的有效传热面积标记为通过公式(14)计算，
式中，
——O点到O1点的距离；
——O点到O2点的距离；
——O点到O3点的距离；
以上的划分网格的方法描述是针对其中一个四面体系统而言的，在计算中需要对i粒子的所有符合的四面体网格进行计算得到有效传热面积数据之后才能开始进行温度计算；
步骤3：在温度计算中，能量守恒方程如公式(15)所示，
式中，
T——温度K；
t——时间s
κ——导热系数W/(m2·K)；
ρ——密度kg/m3；
cp——比定压热容J/(kg·K)；
ST——温度的源项；
公式(16)给出能量方程的Δt时间步长的积分形式，
式中，
t——时间/s；
Δt——时间步长s；
T——温度K；
V——体积m3；
CV——控制体积m3；
κ——导热系数W/(m2·K)；
ρ——密度kg/m3；
cp——比定压热容J/(kg·K)；
ST——温度的源项；
对于粒子i而言，其已经通过步骤2划分为若干个四面体网格，所以公式(16)能够离散成公式(17)的形式，
公式(17)等式的左边对i粒子划分成的多个四面体的能量变化进行加和，等式的右边则是i粒子的能量变化率；
式中：
Sj→i——步骤2中网格划分完成后计算的j粒子对i粒子的有效传热面积；
κ——导热系数W/(m2·K)；
rj——j粒子的坐标矢量；
ri——i粒子的坐标矢量；
——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；
求解公式(17)即得到i粒子在时间步长Δt之后的温度Ti n+1；
步骤4：显式计算动量方程中的粘性项及源项，得到粒子速度及位置的估算值及动量方程如公式(18)所示
式中，
u——速度m/s；
t——时间s；
p——压力Pa
v——运动粘度m2/s；
β——热膨胀系数m-1；
g——重力加速度m2/s；
Tm——背景温度，定义为环境温度K；
τ——粒子应力项，表示不可解的亚粒子尺度的流动对可解尺度的贡献，这一项是用于求解高瑞利数熔融池内湍流的；
显式计算动量方程的粘性项及源项，即计算项，其中粘性项的拉普拉斯算子由公式(19)离散求解，
式中，
d——计算的维度；
n0——初始粒子数密度，通过计算，是核函数，使用公式(20)计算，只需将括号内的自变量修改为即可，为j粒子初始时刻的坐标矢量，为i粒子初始时刻的坐标矢量；
λ——是i粒子的邻居粒子和i粒子间距的方程，使用公式(21)计算，
式中，
re——粒子作用半径，位于半径内的粒子才会与i粒子发生相互作用，，re＝3.1l0；
速度和位置的估算值由公式(22)进行计算，
式中，
——i粒子在上一时间步的速度矢量m/s；
ri n——i粒子在上一时间步的位置矢量m；
Ti n——i粒子在上一时间步的温度K；
步骤5：根据估算得到的粒子位置，通过公式(23)计算此时的粒子数密度n* i，只需要将公式中的|rj-ri|使用位置的估算值替代即可，由公式(24)所示的压力Poisson方程求解得到此时各粒子的压力值pi；
公式(24)等式的左边使用如同公式(19)的形式离散，式中，
——i粒子在n+1时层的压力Pa；
α、γ——静态系数；
步骤6：由计算得到的压力值，根据公式(25)计算速度的修正值；
式中，
u′i——i粒子的速度修正值；
再由计算得到的速度修正值，根据公式(26)，公式(27)分别求解得到粒子在下一时层的速度及位置；
步骤7：输出计算所得到的高瑞利数熔融池中粒子的种类、速度位置压力及温度值并对熔融池内部粒子信息进行处理，按照不同的粒子种类表示不同的物质即得到熔融池内的硬壳分布，而对于壁面热流密度根据公式(28)求解；
式中，
q——穿过壁面的热流密度/W/m2；
κwall——壁面材料的导热系数/W/(m·K)；
——A、B粒子在n+1时层的温度K，A、B粒子分别为位于同一径向的内壁面粒子和外壁面粒子；
——A、B粒子在n+1时层的位置矢量m。
本发明方法为IVR熔融物堆内持留提供解决方案，为核电厂反应堆严重事故安全特性的研究提供重要依据。
和现有技术相比，本发明方法具备如下优点：
本发明的一种研究高瑞利数熔融池换热特性的移动粒子有限容积法，具有有限体积法计算高精度的优点，而且也具备粒子法中精确捕捉自由液面、方便建模不需人工划分网格、以及精确处理相变问题的优点，能够对核反应堆严重事故中下封头形成的高瑞利数熔融池换热特性进行准确的研究，从而为更加全面、有效地评估反应堆的安全性提供依据。
附图说明
图1是本发明利用移动粒子有限容积法研究高瑞利数熔融池换热特性的方法流程图。
图2是i粒子周围1.2倍半径内的随机三个粒子j1、j2、j3位置示意图。
图3是i粒子与周围随机三个构成平面的粒子j1、j2、j3所构成的四面体示意图。
图4是三角形O1O2O3面积划分示意图。
图5是三角形OO1O2面积划分示意图。
图6是三角形OO1O3面积划分示意图。
图7是三角形OO2O3面积划分示意图。
图8是四面体OO1O2O3中顶点O控制体积的体积划分示意图。
具体实施方式
本发明一种利用移动粒子有限容积法研究高瑞利数熔融池换热特性的方法，如图1所示，步骤如下：
步骤1：根据所研究的问题，对高瑞利数熔融池进行粒子建模，不同种类的粒子代表不同的物质，每种粒子根据所表示的物质具有其对应的质量、密度、熔化焓等信息，为方便描述，使用粒子i来表示某一粒子，则其质量、密度、熔化焓对应标记分别为mi、ρi、设定粒子的初始布置，包括粒子的位置、速度、压力和温度，对应标记分别为ri、ui、pi和Ti，粒子初始布置时相邻两粒子间的距离记为l0；
步骤2：根据各粒子的位置划分网格：在本步骤中，设置一个粒子搜索半径，对于每一个粒子i，其坐标位置为(xi，yi，zi)，与其间距lij≤1.2l0的所有粒子作为搜索的对象，式中lij为j粒子与i粒子的距离，j粒子的坐标标记为(xj，yj，zj)，lij可使用公式(1)进行计算，
并将搜索到的粒子使用j1，j2，……，jn符号进行标记，其中下标n表示搜索到的粒子总数；将j1，j2，……，jn粒子中的每三个粒子进行组合，将有种组合，由公式(2)进行计算，
对每个组合中的粒子进行考察，在这些组合中，每个组合内的三个坐标点必须满足构成面的要求，取其中的一个组合进行解析，如其中一种组合粒子标记为j1，j2，j3，粒子中心分别标记为O1，O2，O3，空间坐标对应为(x1，y1，z1)，(x2，y2，z2)，(x3，y3，z3)，如图2所示，连接O1、O2为直线O1O2，连接O1、O3为直线O1O3，连接O2、O3为直线O2O3，若要判断O1、O2、O3能否构成一个平面，只需要判断直线O1O2与直线O1O3之间存在一个不为0的夹角，只需满足公式(3)即可，
对所有的组合进行判断，将能够构成面的组合保留，不能构成面的组合删去，随后考虑i粒子与每个组合中三个粒子的关系，i粒子与组合中三个粒子必须能够构成四面体，因为组合中的三个粒子可以构成平面，所以要使这四个粒子能够构成四面体，只需满足i粒子不在组合中三个粒子所构成的平面上即可，为方便描述，如图3所示，在图2的基础上添加i粒子，粒子中心标记为O，空间坐标定义为(xi，yi，zi)。若要使O、O1、O2、O3四个点构成四面体，因为O1、O2、O3三个点已经满足构成面的要求，只需要满足i粒子到平面O1O2O3的距离不为0即可，可以使用公式(4)进行判定，
式中：
——i粒子到平面O1O2O3的距离m；
a、b、c、d——平面O1O2O3的平面方程的四个系数，可以通过公式(5)获得，
i粒子与j1、j2、j3粒子构成的平面O1O2O3的距离若满足公式(4)，则说明i粒子可以与三角形O1O2O3构成四面体OO1O2O3，四面体OO1O2O3的四个面分别为三角形OO1O2、三角形OO1O3、三角形OO2O3和三角形O1O2O3。使用海伦公式计算四个三角形的面积，四个面的面积对应标记为海伦公式如公式(6)所示，
式中，
S——三角形的面积；
L1、L2、L3——三角形的三边长；
p——三角形的半周长，即p＝(a+b+c)/2。
四面体OO1O2O3的体积可以通过公式(7)进行计算，
式中，
V——四面体体积，在四面体OO1O2O3中，标记为
S——三角形O1O2O3的面积，标记为使用公式(6)计算；
h——O点与三角形所在平面O1O2O3的距离，记为使用公式(4)计算；
在使用粒子位置划分网格的时候，还有几个参量需要求解，其一是四面体OO1O2O3中四个顶点O、O1、O2、O3分别对四面体总体积的分属控制体积，对应标记为VO、这个分属控制体积可以由公式(8)进行计算，
实际上，四面体OO1O2O3的四个点的分属控制体积与其相对的三角形面积成反比，即有公式(9)所示的关系，
其二是对四面体OO1O2O3中四个面的三角形进行面积划分，为方便表示，依然在图2的基础上添加线段进行描述，如图4所示，在三角形O1O2O3中，三角形O1O2O3的面积可以通过公式(6)求出。图4中的点M12为线段O1O2的中点，点M13为线段O1O3的中点，点M23为线段O2O3的中点。此时需要在三角形中寻找一点P，连接P、O1，连接P、O2，连接P、O3，连接P、M12，连接P、M13，连接P、M23，此时三角形O1O2O3被分成了6个小三角形，分别标记为三角形O1PM12，三角形O1PM13，三角形O2PM12，三角形O2PM23，三角形O3PM13，三角形O3PM23；因为点M13为线段O1O3的中点，所以三角形O1PM13与三角形O3PM13面积相等，标记为同样的，三角形O1PM12与三角形O2PM12面积相等，标记为三角形O2PM23与三角形O3PM23面积相等，标记为对于各个三角形的面积，有如公式(10)和公式(11)所示的几何关系，
三角形的三个角的角度可由如公式(3)所示的余弦定理求出。由公式(10)和公式(11)可以求出三角形O1O2O3所分成的6个小三角形的面积。同理也可求得四面体OO1O2O3另外三个面所分成的共计18个小三角形的面积。因为以下的方法描述需要用到每个面上的点，所以将其余三个面的三角形划分图进行展示，辅助点的添加及对应的标注与上述相似。图5为三角形OO1O2的划分图，图6为三角形OO1O3的划分图，图7为三角形OO2O3的划分图。
在四面体OO1O2O3中寻找一点Q，对于顶点O的控制体积而言，如图8所示为整合图5、图6、图7所得到的图形，为方便描述，已去掉不必要的线条，寻找点Q，连接QP1，QP2，QP3，QM1，QM2，QM3，QO，至此点O的控制体积分成了六个四面体，分别为四面体OP1M1Q，四面体OP1M2Q，四面体OP2M1Q，四面体OP2M3Q，四面体OP3M2Q，四面体OP3M3Q，其体积对应标记为 和对于四面体的体积，有公式(12)所示的关系，
设点Q距离平面O1O2O3、平面OO1O2，平面OO1O3，平面OO2O3的距离分别为H、H1、H2、H3，对于顶点O的控制体积可列如公式(13)的方程，
同理对于顶点O1、O2、O3也可列出相似的方程，联立这4个方程及公式(8)可求出H、H1、H2、H3的值。
至此，在O、O1、O2、O3组成的四面体系统中，O1对O起影响的点O的控制体积为记为O2对O起影响的点O的控制体积为记为O3对O起影响的点O的控制体积为记为所以，在O、O1、O2、O3组成的四面体系统中，点O1对O的有效传热面积标记为点O2对O的有效传热面积标记为点O1对O的有效传热面积标记为可以通过公式(14)计算，
式中，
——O点到O1点的距离；
——O点到O2点的距离；
——O点到O3点的距离；
以上的划分网格的方法描述是针对其中一个四面体系统而言的，在计算中需要对i粒子的所有符合的四面体网格进行计算得到有效传热面积数据之后才能开始进行温度计算。
步骤3：在温度计算中，能量守恒方程如公式(15)所示，
式中，
T——温度K；
t——时间s
κ——导热系数W/(m2·K)；
ρ——密度kg/3；
cp——比定压热容J/(kg·K)；
ST——温度的源项。
考虑到本发明方法使用有限体积法求解温度场，且前面已经对计算域进行了网格划分，公式(16)给出能量方程的Δt时间步长的积分形式，
式中，
t——时间s；
Δt——时间步长s；
T——温度K；
V——体积m3；
CV——控制体积m3；
κ——导热系数W/(m2·K)；
ρ——密度kg/m3；
cp——比定压热容J/(kg·K)；
ST——温度的源项。
对于粒子i而言，其已经通过步骤2划分为若干个四面体网格，所以公式(16)可以离散成公式(17)的形式，
公式(17)等式的左边对i粒子划分成的多个四面体的能量变化进行加和，等式的右边则是i粒子的能量变化率。
式中：
Sj→i——步骤2中网格划分完成后计算的j粒子对i粒子的有效传热面积；
κ——导热系数W/(m2·K)；
rj——j粒子的坐标矢量；
ri——i粒子的坐标矢量；
——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；
求解公式(17)即可得到i粒子在时间步长Δt之后的温度Ti n+1。
步骤4：显式计算动量方程中的粘性项及源项，得到粒子速度及位置的估算值及动量方程如公式(18)所示
式中，
u——速度m/s；
t——时间s；
p——压力Pa
v——运动粘度m2/s；
β——热膨胀系数m-1；
g——重力加速度m2/s；
Tm——背景温度，可定义为环境温度K；
τ——粒子应力项，表示不可解的亚粒子尺度的流动对可解尺度的贡献，这一项是用于求解高瑞利数熔融池内湍流的。
显式计算动量方程的粘性项及源项，即计算项，其中粘性项的拉普拉斯算子由公式(19)离散求解，
式中，d——计算的维度；
n0——初始粒子数密度，可通过计算，是核函数，使用公式(20)计算，只需将括号内的自变量修改为即可，为j粒子初始时刻的坐标矢量，为i粒子初始时刻的坐标矢量；
λ——是i粒子的邻居粒子和i粒子间距的方程，使用公式(21)计算，
式中，
re——粒子作用半径，位于半径内的粒子才会与i粒子发生相互作用，在本方法中，re＝3.1l0。
速度和位置的估算值ri *可由公式(22)进行计算，
式中，
——i粒子在上一时间步的速度矢量m/s；
ri n——i粒子在上一时间步的位置矢量m；
Ti n——i粒子在上一时间步的温度K；
步骤5：根据估算得到的粒子位置，通过公式(23)计算此时的粒子数密度n* i，只需要将公式中的|rj-ri|使用位置的估算值替代即可，由公式(24)所示的压力Poisson方程求解得到此时各粒子的压力值pi。
公式(24)等式的左边可以使用如同公式(19)的形式离散，式中，
——i粒子在n+1时层的压力/Pa；
α、γ——静态系数，对于不同的问题有不同的值。
步骤6：由计算得到的压力值，根据公式(25)计算速度的修正值。
式中，u′i——i粒子的速度修正值。
再由计算得到的速度修正值，根据公式(26)，公式(27)分别求解得到粒子在下一时层的速度及位置。
步骤7：输出计算所得到的高瑞利数熔融池中粒子的种类、速度位置压力及温度值Ti n+1。并根据所需要分析的问题对熔融池内部粒子信息进行处理，
按照不同的粒子种类表示不同的物质即得到熔融池内的硬壳分布，而对于壁面热流密度根据公式(28)求解；
式中，q——穿过壁面的热流密度/W/m2；
κwall——壁面材料的导热系数/W/(m·K)；
——A、B粒子在n+1时层的温度K，A、B粒子分别为位于同一径向的内壁面粒子和外壁面粒子；
——A、B粒子在n+1时层的位置矢量m。本发明方法为IVR熔融物堆内持留提供解决方案，为核电厂反应堆严重事故安全特性的研究提供重要依据。

Claims (1)

1.一种利用移动粒子有限容积法研究高瑞利数熔融池换热特性的方法，其特征在于：步骤如下：
步骤1：对高瑞利数熔融池进行粒子建模，不同种类的粒子代表不同的物质，每种粒子根据所表示的物质具有其对应的质量、密度和熔化焓，使用粒子i来表示某一粒子，则其质量、密度、熔化焓对应标记分别为mi、ρi、设定粒子的初始布置，包括粒子的位置、速度、压力和温度，对应标记分别为ri、ui、pi和Ti，粒子初始布置时相邻两粒子间的距离记为l0；
步骤2：根据各粒子的位置划分网格：在本步骤中，设置一个粒子搜索半径，对于每一个粒子i，其坐标位置为(xi，yi，zi)，与其间距lij≤1.2l0的所有粒子作为搜索的对象，式中lij为j粒子与i粒子的距离，j粒子的坐标标记为(xj，yj，zj)，lij使用公式(1)进行计算，
并将搜索到的粒子使用j1，j2，……，jn符号进行标记，其中下标n表示搜索到的粒子总数；将j1，j2，……，jn粒子中的每三个粒子进行组合，将有种组合，由公式(2)进行计算，
对每个组合中的粒子进行考察，在这些组合中，每个组合内的三个坐标点必须满足构成面的要求，取其中的一个组合进行解析，其中一种组合粒子标记为j1，j2，j3，粒子中心分别标记为O1，O2，O3，空间坐标对应为(x1，y1，z1)，(x2，y2，z2)，(x3，y3，z3)，连接O1、O2为直线O1O2，连接O1、O3为直线O1O3，连接O2、O3为直线O2O3，若要判断O1、O2、O3能否构成一个平面，只需要判断直线O1O2与直线O1O3之间存在一个不为0的夹角，则只需满足公式(3)，
对所有的组合进行判断，将能够构成面的组合保留，不能构成面的组合删去，随后考虑i粒子与每个组合中三个粒子的关系，i粒子与组合中三个粒子必须能够构成四面体，因为组合中的三个粒子能够构成平面，所以要使这四个粒子能够构成四面体，只需满足i粒子不在组合中三个粒子所构成的平面上即可，添加i粒子，粒子中心标记为O，空间坐标定义为(xi，yi，zi)，若要使O、O1、O2、O3四个点构成四面体，因为O1、O2、O3三个点已经满足构成面的要求，则只需要满足i粒子到平面O1O2O3的距离不为0，使用公式(4)进行判定，
式中：
——i粒子到平面O1O2O3的距离m；
a、b、c、d——平面O1O2O3的平面方程的四个系数，通过公式(5)获得，
i粒子与j1、j2、j3粒子构成的平面O1O2O3的距离若满足公式(4)，则说明i粒子能够与三角形O1O2O3构成四面体OO1O2O3，四面体OO1O2O3的四个面分别为三角形OO1O2、三角形OO1O3、三角形OO2O3和三角形O1O2O3；使用海伦公式计算四个三角形的面积，四个面的面积对应标记为海伦公式如公式(6)所示，
式中，
S——三角形的面积；
L1、L2、L3——三角形的三边长；
p——三角形的半周长，即p＝(L1+L2+L3)/2；
四面体OO1O2O3的体积通过公式(7)进行计算，
式中，
V——四面体体积，在四面体OO1O2O3中，标记为
S——三角形O1O2O3的面积，标记为使用公式(6)计算；
h——O点与三角形所在平面O1O2O3的距离，记为使用公式(4)计算；
在使用粒子位置划分网格时，还有几个参量需要求解，其一是四面体OO1O2O3中四个顶点O、O1、O2、O3分别对四面体总体积的分属控制体积，对应标记为VO、这个分属控制体积由公式(8)进行计算，
实际上，四面体OO1O2O3的四个点的分属控制体积与其相对的三角形面积成反比，即有公式(9)所示的关系，
其二是对四面体OO1O2O3中四个面的三角形进行面积划分，在三角形O1O2O3中，三角形O1O2O3的面积通过公式(6)求出；点M12为线段O1O2的中点，点M13为线段O1O3的中点，点M23为线段O2O3的中点。此时需要在三角形中寻找一点P，连接P、O1，连接P、O2，连接P、O3，连接P、M12，连接P、M13，连接P、M23，此时三角形O1O2O3被分成了6个小三角形，分别标记为三角形O1PM12，三角形O1PM13，三角形O2PM12，三角形O2PM23，三角形O3PM13，三角形O3PM23；因为点M13为线段O1O3的中点，所以三角形O1PM13与三角形O3PM13面积相等，标记为同样的，三角形O1PM12与三角形O2PM12面积相等，标记为三角形O2PM23与三角形O3PM23面积相等，标记为对于各个三角形的面积，有如公式(10)和公式(11)所示的几何关系，
三角形的三个角的角度由公式(3)所示的余弦定理求出；由公式(10)和公式(11)求出三角形O1O2O3所分成的6个小三角形的面积；同理也能求得四面体OO1O2O3另外三个面所分成的共计18个小三角形的面积；
在四面体OO1O2O3中寻找一点Q，对于顶点O的控制体积而言，连接QP1，QP2，QP3，QM1，QM2，QM3，QO，至此点O的控制体积分成了六个四面体，分别为四面体OP1M1Q，四面体OP1M2Q，四面体OP2M1Q，四面体OP2M3Q，四面体OP3M2Q和四面体OP3M3Q，其体积对应标记为 和对于四面体的体积，有公式(12)所示的关系，
设点Q距离平面O1O2O3、平面OO1O2，平面OO1O3，平面OO2O3的距离分别为H、H1、H2、H3，对于顶点O的控制体积列如公式(13)的方程，
同理对于顶点O1、O2、O3也列出相似的方程，联立这4个方程及公式(8)能够求出H、H1、H2、H3的值；
至此，在O、O1、O2、O3组成的四面体系统中，O1对O起影响的点O的控制体积为记为O2对O起影响的点O的控制体积为记为O3对O起影响的点O的控制体积为记为所以，在O、O1、O2、O3组成的四面体系统中，点O1对O的有效传热面积标记为点O2对O的有效传热面积标记为点O1对O的有效传热面积标记为通过公式(14)计算，
式中，
——O点到O1点的距离；
——O点到O2点的距离；
——O点到O3点的距离；
以上的划分网格的方法描述是针对其中一个四面体系统而言的，在计算中需要对i粒子的所有符合的四面体网格进行计算得到有效传热面积数据之后才能开始进行温度计算；
步骤3：在温度计算中，能量守恒方程如公式(15)所示，
式中，
T——温度K；
t——时间s
κ——导热系数W/(m2·K)；
ρ——密度kg/m3；
cp——比定压热容J/(kg·K)；
ST——温度的源项；
公式(16)给出能量方程的Δt时间步长的积分形式，
式中，
t——时间/s；
Δt——时间步长s；
T——温度K；
V——体积m3；
CV——控制体积m3；
κ——导热系数W/(m2·K)；
ρ——密度kg/m3；
cp——比定压热容J/(kg·K)；
ST——温度的源项；
对于粒子i而言，其已经通过步骤2划分为若干个四面体网格，所以公式(16)能够离散成公式(17)的形式，
公式(17)等式的左边对i粒子划分成的多个四面体的能量变化进行加和，等式的右边则是i粒子的能量变化率；
式中：
Sj→i——步骤2中网格划分完成后计算的j粒子对i粒子的有效传热面积；
κ——导热系数W/(m2·K)；
rj——j粒子的坐标矢量；
ri——i粒子的坐标矢量；
——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；
求解公式(17)即得到i粒子在时间步长Δt之后的温度Ti n+1；
步骤4：显式计算动量方程中的粘性项及源项，得到粒子速度及位置的估算值及ri *；动量方程如公式(18)所示
式中，
u——速度m/s；
t——时间s；
p——压力Pa
v——运动粘度m2/s；
β——热膨胀系数m-1；
g——重力加速度m2/s；
Tm——背景温度，定义为环境温度K；
τ——粒子应力项，表示不可解的亚粒子尺度的流动对可解尺度的贡献，这一项是用于求解高瑞利数熔融池内湍流的；
显式计算动量方程的粘性项及源项，即计算项，其中粘性项的拉普拉斯算子由公式(19)离散求解，
式中，
d——计算的维度；
n0——初始粒子数密度，通过计算，是核函数，使用公式(20)计算，只需将括号内的自变量修改为即可，为j粒子初始时刻的坐标矢量，为i粒子初始时刻的坐标矢量；
λ——是i粒子的邻居粒子和i粒子间距的方程，使用公式(21)计算，
式中，
re——粒子作用半径，位于半径内的粒子才会与i粒子发生相互作用，，re＝3.1l0；
速度和位置的估算值由公式(22)进行计算，
式中，
——i粒子在上一时间步的速度矢量m/s；
ri n——i粒子在上一时间步的位置矢量m；
Ti n——i粒子在上一时间步的温度K；
步骤5：根据估算得到的粒子位置，通过公式(23)计算此时的粒子数密度n* i，只需要将公式中的|rj-ri|使用位置的估算值替代即可，由公式(24)所示的压力Poisson方程求解得到此时各粒子的压力值pi；
公式(24)等式的左边使用如同公式(19)的形式离散，式中，
——i粒子在n+1时层的压力Pa；
α、γ——静态系数；
步骤6：由计算得到的压力值，根据公式(25)计算速度的修正值；
式中，
u′i——i粒子的速度修正值；
再由计算得到的速度修正值，根据公式(26)，公式(27)分别求解得到粒子在下一时层的速度及位置；
ri n+1＝ri *+u′iΔt 公式(27)
步骤7：输出计算所得到的高瑞利数熔融池中粒子的种类、速度位置压力及温度值Ti n+1；并对熔融池内部粒子信息进行处理，按照不同的粒子种类表示不同的物质即得到熔融池内的硬壳分布，而对于壁面热流密度根据公式(28)求解；
式中，
q——穿过壁面的热流密度/W/m2；
κwall——壁面材料的导热系数/W/(m·K)；
——A、B粒子在n+1时层的温度K，A、B粒子分别为位于同一径向的内壁面粒子和外壁面粒子；
——A、B粒子在n+1时层的位置矢量m。
CN201910394189.XA 2019-05-13 2019-05-13 利用移动粒子有限容积法研究熔融池换热特性的方法 Active CN110044959B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910394189.XA CN110044959B (zh) 2019-05-13 2019-05-13 利用移动粒子有限容积法研究熔融池换热特性的方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910394189.XA CN110044959B (zh) 2019-05-13 2019-05-13 利用移动粒子有限容积法研究熔融池换热特性的方法

Publications (2)

Publication Number Publication Date
CN110044959A true CN110044959A (zh) 2019-07-23
CN110044959B CN110044959B (zh) 2020-05-15

Family

ID=67281564

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910394189.XA Active CN110044959B (zh) 2019-05-13 2019-05-13 利用移动粒子有限容积法研究熔融池换热特性的方法

Country Status (1)

Country Link
CN (1) CN110044959B (zh)

Cited By (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110489712A (zh) * 2019-08-08 2019-11-22 中国核动力研究设计院 基于测量数据获取熔融池表面冷却热流密度的方法
CN110867220A (zh) * 2019-11-07 2020-03-06 西安交通大学 利用粒子网格混合法研究堆芯内共晶反应及高温熔化行为的方法
CN111832214A (zh) * 2020-06-29 2020-10-27 西安交通大学 基于粒子法的核反应堆严重事故碎片床熔化过程模拟方法
CN112102894A (zh) * 2020-09-04 2020-12-18 西安交通大学 基于粒子法的核反应堆堆芯材料熔池演变特性分析方法
CN113191066A (zh) * 2021-04-30 2021-07-30 西安交通大学 基于无网格法的核反应堆燃料元件失效分析方法
CN113435004A (zh) * 2021-05-25 2021-09-24 上海交通大学 堆芯熔融进程中熔融材料迁移质量的计算方法和装置

Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20070129919A1 (en) * 2005-10-26 2007-06-07 Akira Toda Fluid-structure coupled numerical simulation method and program for fluid-structure coupled numerical simulation storage device
CN102262689A (zh) * 2010-05-26 2011-11-30 利弗莫尔软件技术公司 能够实现实体/sph耦合作用的混合单元
US20120014491A1 (en) * 2009-07-13 2012-01-19 Mike Deeth Nuclear fusion power plant having a liquid reactor core of molten glass that is made laseractive and functions as a tritium breeding blanket which is capable of acousticly compressing/confining fuel so that it radiates and triggers outgoing laser cascades that will reflect from the blast chamber's spherical inside wall and return like photonic Tsunamis, crushing, heating, and causing thermonuclear ignition of the fuel so that heat engines and piezoelectric harvesters can convert the released energy into electricity
CN103594133A (zh) * 2013-10-21 2014-02-19 西安交通大学 模拟核反应堆堆芯熔化后堆内熔融物滞留的实验装置及方法
CN104657589A (zh) * 2015-01-16 2015-05-27 华北电力大学 一种管壳式换热器多孔介质系数计算方法
CN105806881A (zh) * 2016-05-11 2016-07-27 西安交通大学 一种高瑞利数熔融池换热特性测量实验装置
CN106055850A (zh) * 2016-07-18 2016-10-26 西安交通大学 一种获得偏离泡核沸腾型临界热流密度的方法
CN107563030A (zh) * 2017-08-22 2018-01-09 哈尔滨工程大学 一种针对两种流体传热交混破碎相变过程的无网格模拟方法

Patent Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20070129919A1 (en) * 2005-10-26 2007-06-07 Akira Toda Fluid-structure coupled numerical simulation method and program for fluid-structure coupled numerical simulation storage device
US20120014491A1 (en) * 2009-07-13 2012-01-19 Mike Deeth Nuclear fusion power plant having a liquid reactor core of molten glass that is made laseractive and functions as a tritium breeding blanket which is capable of acousticly compressing/confining fuel so that it radiates and triggers outgoing laser cascades that will reflect from the blast chamber's spherical inside wall and return like photonic Tsunamis, crushing, heating, and causing thermonuclear ignition of the fuel so that heat engines and piezoelectric harvesters can convert the released energy into electricity
CN102262689A (zh) * 2010-05-26 2011-11-30 利弗莫尔软件技术公司 能够实现实体/sph耦合作用的混合单元
CN103594133A (zh) * 2013-10-21 2014-02-19 西安交通大学 模拟核反应堆堆芯熔化后堆内熔融物滞留的实验装置及方法
CN104657589A (zh) * 2015-01-16 2015-05-27 华北电力大学 一种管壳式换热器多孔介质系数计算方法
CN105806881A (zh) * 2016-05-11 2016-07-27 西安交通大学 一种高瑞利数熔融池换热特性测量实验装置
CN106055850A (zh) * 2016-07-18 2016-10-26 西安交通大学 一种获得偏离泡核沸腾型临界热流密度的方法
CN107563030A (zh) * 2017-08-22 2018-01-09 哈尔滨工程大学 一种针对两种流体传热交混破碎相变过程的无网格模拟方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
S.HEO,ET AL.: "Numerical analysis of boiling on high heat-flux and high subcooling condition using MPS-MAFL", 《INTERNATIONAL JOURNAL OF HEAT AND MASS TRANSFER》 *
ZHONGYI LIU,ET AL.: "An Integrated Method of FVM and SPH for Treating Melting Process of Quartz ingot", 《INTERNATIONAL MECHANICAL ENGINEERING CONGRESS AND EXPOSITION》 *
李晨曦 等: "基于FOCUS程序的ADS颗粒流靶熔化事故分析", 《原子能科学技术》 *

Cited By (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110489712A (zh) * 2019-08-08 2019-11-22 中国核动力研究设计院 基于测量数据获取熔融池表面冷却热流密度的方法
CN110867220A (zh) * 2019-11-07 2020-03-06 西安交通大学 利用粒子网格混合法研究堆芯内共晶反应及高温熔化行为的方法
CN111832214A (zh) * 2020-06-29 2020-10-27 西安交通大学 基于粒子法的核反应堆严重事故碎片床熔化过程模拟方法
CN112102894A (zh) * 2020-09-04 2020-12-18 西安交通大学 基于粒子法的核反应堆堆芯材料熔池演变特性分析方法
CN112102894B (zh) * 2020-09-04 2021-11-30 西安交通大学 基于粒子法的核反应堆堆芯材料熔池演变特性分析方法
CN113191066A (zh) * 2021-04-30 2021-07-30 西安交通大学 基于无网格法的核反应堆燃料元件失效分析方法
CN113435004A (zh) * 2021-05-25 2021-09-24 上海交通大学 堆芯熔融进程中熔融材料迁移质量的计算方法和装置

Also Published As

Publication number Publication date
CN110044959B (zh) 2020-05-15

Similar Documents

Publication Publication Date Title
CN110044959A (zh) 利用移动粒子有限容积法研究高瑞利数熔融池换热特性的方法
Nithyadevi et al. Double diffusive natural convection in a partially heated enclosure with Soret and Dufour effects
Wu An approach combining body-fitted grid generation and conjugate gradient methods for shape design in heat conduction problems
Hossain et al. Buoyancy and thermocapillary driven convection flow of an electrically conducting fluid in an enclosure with heat generation
CN108416112B (zh) 压力容器下封头内多层熔融池传热分析方法
Zhang et al. The COPRA experiments on the in-vessel melt pool behavior in the RPV lower head
CN104834781A (zh) 一种基于平波电抗器多次换向失败时的暂态温度场仿真方法
CN113095004B (zh) 一种液态金属流动传热计算方法
CN112102894B (zh) 基于粒子法的核反应堆堆芯材料熔池演变特性分析方法
Sagara et al. First operation of the Flinak/LiPb twin loop OROSH2I-2 with a 3T SC magnet for R&D of liquid blanket for fusion reactor
CN103433448A (zh) 基于渣膜与气隙动态分布的连铸结晶器热流密度确定方法
CN106709137A (zh) 基于多项式极大似然估计的核电站堆芯温度场软测量方法
CN110867220A (zh) 利用粒子网格混合法研究堆芯内共晶反应及高温熔化行为的方法
CN111832214B (zh) 基于粒子法的核反应堆严重事故碎片床熔化过程模拟方法
Hlimi et al. Melting inside a horizontal cylindrical capsule
CN106897520A (zh) 一种含有模糊参数的传热系统可靠性分析方法
Capone et al. Reactor cavity cooling system (Rccs) experimental characterization
CN103983661A (zh) 强内热源大尺寸小球球床通道熔盐强迫对流换热实验装置
CN105808820B (zh) 一种求解区间热对流扩散问题的高精度数值方法
Shen et al. Soft sensor modeling of blast furnace wall temperature based on temporal–spatial dimensional finite-element extrapolation
Zhong et al. Experimental study of saturated pool boiling from downward facing structured surfaces
CN107084805A (zh) 一种确定连铸钢包辐射热流的方法
Chokri et al. Three-dimensional natural convection of molten Lithium in a differentially heated rotating cubic cavity about a vertical ridge
Raič et al. Coupled CFD simulation for the cooling process of steel tubes in a rake type cooling bed
Orlov et al. Mathematical modeling of process of obtaining the solid uranium hexafluoride

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