一种基于最小二乘粒子法的堆芯熔化收集器设计优化方法技术

本发明专利技术公开了一种基于最小二乘粒子法的堆芯熔化收集器设计优化方法，该方法基于拉格朗日方法中具有更高模拟精度且没有空间一致性限制的最小二乘移动粒子法，在此基础上添加流固耦合模型(PMS模型)用于计算冷却剂与固体熔融物之间的相互作用，最后通过固体碰撞模型(DEM方法)计算流场中熔融物固体颗粒之间、熔融物固体颗粒与固体壁面之间的作用力。此外，采用优化粒子滑移模型缓解粒子粒子聚集效应；利用光滑化处理计算密度和动力粘度以提高计算稳定性；通过代数计算和几何计算的双检测方法来提高粒子类型判定的准确性。本发明专利技术方法是一种具备模拟堆芯熔化收集器中流固耦合问题和固体碰撞问题能力的方法，以用于堆芯熔化收集器的设计优化。集器的设计优化。集器的设计优化。

【技术实现步骤摘要】
一种基于最小二乘粒子法的堆芯熔化收集器设计优化方法


[0001]本专利技术属于核反应堆安全设施设计
，具体涉及一种基于最小二乘粒子法的堆芯熔化收集器设计优化方法。

技术介绍

[0002]钠冷快堆在发生严重事故时，由于其燃料富集度高以及冷却剂特性，堆芯解体事故尤其值得关注。此时，熔融物堆内滞留是保证压力容器完整性的重要缓解措施。然而，堆芯熔融物在下腔室的分布和堆积厚度严重影响其余热排出能力，设置一种具有合理的结构和布置的堆芯熔化收集器是一种有效的缓解方法。因此，堆芯熔化收集器的设计对钠冷快堆严重事故管理具有重要意义。
[0003]堆芯熔化收集器与堆芯熔融物相互作用的过程涉及强烈流固耦合和固体碰撞问题。传统的欧拉方法在使用网格处理流固耦合问题时，会面临相界面重构以及自由界面等难点，且收集器的设计中需要清晰地观察熔融物在堆芯熔化收集器上堆积的形状，因此欧拉方法不是最佳选择。而使用粒子方式计算流体的拉格朗日方法在处理上述问题时显示出特有的优势，其中MPS方法在发展过程中开发了PMS模型、表面张力模型、粘性力模型等一系列相关模型，在解决流固耦合或固体碰撞问题时已经取得了不错的成就，但传统的MPS方法无法摆脱离散算子在粒子不对称分布下误差过大的弊病，相关模型也无法同时对堆芯熔化收集器中固液耦合和固体碰撞问题进行较好的模拟。

技术实现思路

[0004]本专利技术所要解决的问题是在针对上述现有方法的不足，提供一种基于最小二乘粒子法的堆芯熔化收集器设计优化方法，是一种具备模拟堆芯熔化收集器中流固耦合问题和固体碰撞问题能力的方法，以用于堆芯熔化收集器的设计优化。该方法基于拉格朗日方法中具有更高模拟精度且没有空间一致性限制的最小二乘移动粒子法(LSMPS)，在此基础上添加流固耦合模型(PMS模型)、固体碰撞模型(DEM方法)形成最小二乘移动粒子
‑
离散要素法(LSMPS
‑
DEM)，对钠冷快堆堆芯熔融物在堆芯熔化收集器中的分布和堆积厚度进行模拟研究。
[0005]本专利技术采用以下技术方案：
[0006]一种基于最小二乘粒子法的堆芯熔化收集器设计优化方法，包括以下步骤：
[0007]S1、根据堆芯熔化收集器的模拟工况设置不同粒子的初始位置、速度、压力、相、粒子边界类型以及固体颗粒编号信息，并对设置信息进行预处理；
[0008]S2、每隔十个时间步长调用Link
‑
list算法，对步骤S1预处理后的每个粒子以及其截止半径内的粒子进行检索并存储粒子编号；
[0009]S3、采用几何计算和代数计算结合的双检测方法对粒子表面类型进行判定；
[0010]S4、根据S3步骤判定的粒子表面类型计算最小二乘移动粒子法所必需的微小修正位移量，保证计算精度的同时避免粒子位移扭曲；
[0011]S5、进入光滑化处理子程，对所有粒子的密度、动力粘度系数进行光滑化处理，降低压力计算不稳定的可能性；
[0012]S6、显式计算包括固体粒子在内的所有粒子的重力、粘性力、表面张力；忽略压力的影响，利用重力、粘性力以及表面张力求解粒子临时速度和临时位置；
[0013]S7、使用最小二乘移动粒子法的高精度离散算子对压力泊松方程中的项进行离散，再使用稳定双共轭梯度法隐式求解压力泊松方程得到下一时刻的压力值，根据压力值更新粒子的位置和速度；
[0014]S8、进入流固耦合子程对堆芯熔融物固体颗粒进行修正，计算堆芯熔融物固体颗粒的重心位置、速度以及角速度；根据计算的堆芯熔融物固体颗粒的重心位置、速度以及角速度对熔融物固体颗粒中所含固体粒子的速度和位置进行修正；
[0015]S9、进入固体碰撞子程循环，根据步骤S8计算出的熔融物固体颗粒的重心、速度以及角速度，计算固体碰撞时间步长内熔融物固体颗粒之间的相互作用力，由此更新熔融物固体颗粒的位置和速度，固体碰撞计算程序的时间步长要远小于最小二乘法离散时间步长，循环相同的时间退出固体碰撞计算程序进入下一最小二乘法离散时间步。
[0016]具体的，步骤S2中，用Link
‑
list算法来提高粒子检索效率，具体操作如下：
[0017]每个时步内对所有粒子包括墙边界粒子的区域划分网格，然后统计每个网格内的粒子以及粒子编号；如果是固体颗粒则统计固体颗粒的质心坐标；从而获得每个粒子或者固体颗粒质心所在的网格以及网格坐标；由于粒子的检索范围是一个有限大小的圆形域，在此基础上，每个粒子仅需搜索该粒子所在的网格以及周围的网格。
[0018]具体的，步骤S3中，采用几何计算和代数计算结合的双检测方法对粒子表面类型进行判定以提高效率和准确性，具体如下：
[0019]将粒子表面类型分为内部粒子、靠近表面的内部粒子、表面粒子和孤立粒子；首先通过几何计算，以判别粒子为原点生成N0个象限，计算该判别粒子所有邻域粒子的投影角θ
j
以及所属象限N
j
，计算公式如下：
[0020][0021][0022]其中，r
β
为粒子间距在y方向上的投影，r
α
为粒子间距在x方向上的投影，θ0为每个象限包含的角度大小；
[0023]设粒子i有邻域粒子存在的象限个数为Γ(i)，则没有邻域粒子分布的象限比例为
[0024][0025]当Ψ(i)大于0.2时，粒子i被初步判定为表面粒子；否则为内部粒子；
[0026]几何计算步骤结束后，进行代数计算；首先计算粒子的邻域粒子数N
sum
和粒子数密度n
′
，即
[0027][0028]其中w
ij
为核函数；当N
sum
＜6或n
′
＜0.6n0时，粒子被判定为孤立粒子；其中,n0表示初始粒子数密度；表面粒子和孤立粒子判断完成后，对靠近表面的内部粒子进行判定，所采用的判据为：
[0029][0030]ifsurf(i)为粒子i的粒子类型判定值，其值0、1、2、3分别表示粒子是内部粒子、靠近表面的内部粒子、表面粒子和孤立粒子。
[0031]具体的，步骤S4中，采用了优化粒子滑移模型计算最小二乘移动粒子法所必须的微小修正位移量：
[0032][0033]其中，δr
i
表示粒子i的修正位移，C
shift
取值小于0.5，h是扫描半径，C
i
为粒子i处的粒子浓度；表示粒子i处浓度梯度方向的单位法向量。
[0034]具体的，步骤S5中，通过对计算粒子i邻域粒子的密度和动力粘度系数进行加权平均，求得粒子i经过光滑处理后的粒子密度和动力粘度系数，计算公式如下：
[0035][0036][0037]其中和分别为粒子i经过光滑处理后的粒子密度和动力粘度系数，j表示邻域粒子，ρ
j
为粒子j的密度，为光滑化处理子程中使用的核函数，计算模型如下：
[0038][0039]其中，其中本文档来自技高网...

【技术保护点】

【技术特征摘要】
1.一种基于最小二乘粒子法的堆芯熔化收集器设计优化方法，其特征在于，包括以下步骤：S1、根据堆芯熔化收集器的模拟工况设置不同粒子的初始位置、速度、压力、相、粒子边界类型以及固体颗粒编号信息，并对设置信息进行预处理；S2、每隔十个时间步长调用Link
‑
list算法，对步骤S1预处理后的每个粒子以及其截止半径内的粒子进行检索并存储粒子编号；S3、采用几何计算和代数计算结合的双检测方法对粒子表面类型进行判定；S4、根据S3步骤判定的粒子表面类型计算最小二乘移动粒子法所必需的微小修正位移量，保证计算精度的同时避免粒子位移扭曲；S5、进入光滑化处理子程，对所有粒子的密度、动力粘度系数进行光滑化处理，降低压力计算不稳定的可能性；S6、显式计算包括固体粒子在内的所有粒子的重力、粘性力、表面张力；忽略压力的影响，利用重力、粘性力以及表面张力求解粒子临时速度和临时位置：响，利用重力、粘性力以及表面张力求解粒子临时速度和临时位置：其中，右侧括号内第一项为粘性力，表示表面张力和重力的矢量和，v为运动粘度，为粒子i的临时速度，为粒子i的临时位置，Δt为计算时间步长，为n时刻粒子i的速度，为n时刻粒子i的位置；S7、使用最小二乘移动粒子法的高精度离散算子对压力泊松方程中的项进行离散，再使用稳定双共轭梯度法隐式求解压力泊松方程得到下一时刻的压力值，根据压力值更新粒子的位置和速度：子的位置和速度：子的位置和速度：其中，表示由压力梯度引起的粒子i的速度变化量，ρ表示密度，P
in+1
、分别表示n+1时刻粒子i的压力、速度和位置；S8、进入流固耦合子程对堆芯熔融物固体颗粒进行修正，计算堆芯熔融物固体颗粒的重心位置、速度以及角速度：重心位置、速度以及角速度：
其中，为固体颗粒ii在n时刻的质心坐标，分别为n+1时刻的固体颗粒ii的速度、角速度；为n时刻固体颗粒ii的角速度；N为固体颗粒中包含的固体粒子的数量；m
i
为粒子i的质量；分别表示固体颗粒中固体粒子i在n时刻和n+1时刻的速度，在进入该步骤进行计算前，固体颗粒中的固体粒子被当作液体粒子进行计算；根据堆芯熔融物固体颗粒的重心位置、速度以及角速度对熔融物固体颗粒中粒子i的速度和位置进行修正，计算公式如下：正，计算公式如下：其中，坐标转换矩阵M的形式为：其中，为熔融物固体颗粒在n时刻角速度的大小；S9、根据步骤S8计算的熔融物固体颗粒的速度，熔融物固体颗粒将运动一个固体碰撞计算时间步长从而得到下一计算时刻的位置；根据这个位置，固体碰撞模块将计算每个熔融物固体颗粒与搜索域内的熔融物固体颗粒和收集器壁面的接触情况，并根据弹簧
‑
阻尼模型计算固体颗粒相互碰撞的作用力：模型计算固体颗粒相互碰撞的作用力：模型计算固体颗粒相互碰撞的作用力：其中，分别为熔融物固体颗粒i和熔融物固体颗粒j之间法向和切向的接触力；δ
n
、δ
s
分别是熔融物固体颗粒间接触点法向和切向的相对位移，k
n
、k
n
分别是法向和切向的刚度；c
n
、c
s
分别是法向和切向阻尼系数；μ为固体颗粒间的摩擦系数；t表示时间；产生重叠之后的熔融物固体颗粒会产生推力从而推开彼此，推开的加速度为：之后的熔融物固体颗粒会产生推力从而推开彼此，推开的加速度为：
其中是熔融物固体颗粒ii在k时刻的线加速度和角加速度，和分别是固体碰撞过程中k时刻作用在熔融物固体颗粒ii质心上的作用力和力矩，m
ii
是熔融物固体颗粒ii的质量；I
ii
为熔融物固体颗粒ii的转动惯量；利用计算出的加速度和角加速度，更新下一固体碰撞计算时刻的位置、角度、速度和角速度：速度：速度：速度：其中，分别表示固体颗粒ii在k+1时刻的速度、角速度、角度和位移；Δt
DEM
表...

【专利技术属性】
技术研发人员：张斌，胡铭豪，任前永，王文鹏，曹胜，单建强，
申请(专利权)人：西安交通大学，
类型：发明
国别省市：