一种针对两种流体传热交混破碎相变过程的无网格模拟方法与流程

本发明涉及的是多相流传热相变模拟方法，具体地说是一种针对两种流体传热交混破碎相变过程的无网格模拟方法。

背景技术：

金属燃料堆芯破损事故及其导致的熔融金属材料与冷却剂接触的相互作用fci(fuel-coolantinteraction)，是国际上金属燃料钠冷快堆堆芯安全研究面临的一个重点和难点问题。在fci过程中，两种或两种以上的热熔融物质与冷却剂液态钠之间的巨大温差，强化了冷热流体间的多相传热传质，使得流体接触界面发生剧烈相变，并伴随能量的瞬间释放，以及熔融物质的变形、破碎和凝固相变等，使得fci过程极为复杂。

目前国内关于金属燃料与冷却剂钠相互作用的模拟程序比较缺乏，而国外的模拟程序主要关注fci过程中，温度变化以及产生的压力脉冲,对于熔融金属材料的凝固破碎尺寸并没有太多的模拟研究。传统网格面对此类存在自由界面的流体传热交混破碎相变问题时会出现网格大变形以及数值扩散等问题，而无网格方法可以很好的解决这一问题。

技术实现要素：

本发明的目的在于提供一种可实现两种流体传热交混破碎相变过程，尤其针对两相流动、凝固传热以及流固混合等复杂问题的模拟的针对两种流体传热交混破碎相变过程的无网格模拟方法。

本发明的目的是这样实现的：

步骤一，设定时间步长δt、粒子间初始距离d0、粒子的控制半径re以及总的模拟时间ttotal，设定tn＝n·δt，则在t0时刻输入初始化参数，包括熔融金属流体的初始温度th、速度v0、形态，冷却流体的初始温度tc；

步骤二，在tn时刻采用基于焓方法的凝固相变模型进行传热计算，n＞0，根据tn-1时刻的粒子温度t^n-1、焓值h^n-1、位置r^n-1、密度ρ^n-1以及热导率k^n-1，计算得到tn时刻粒子的温度tⁿ、焓值hⁿ以及液相分数αⁿ；

步骤三，采用物性参数变化模型对处于相变区域粒子的物性参数进行修正，根据粒子的液相分数αⁿ，更新得到tn时刻粒子的密度ρⁿ、热导率kⁿ以及粘性系数vⁿ；

步骤四，采用基于表面自由能的表面张力模型计算tn时刻粒子所受的表面张力

步骤五，根据步骤三和步骤四所得的粒子物性参数ρⁿ、vⁿ和表面张力采用原始mps方法计算出tn时刻所有粒子的位置rⁿ和速度uⁿ，分为以下几个小步骤：

1，不考虑压力梯度项，显式求解动量方程中的粘性项、表面张力项以及重力项，计算得出tn时刻粒子的中间速度u^*和中间位置r^*；

2，计算压力泊松方程，得出tn时刻粒子的压力值pⁿ；

3，利用粒子的压力值pⁿ，对粒子速度进行修正，得到tn时刻粒子的速度uⁿ和位置rⁿ；

步骤六，根据步骤五中计算所得粒子位置rⁿ，采用碎片结合判据确定破碎相变产生碎片的粒子组成；包括液相分数判据和粒子间距离判据，只有两个判据同时满足，才能判定两个粒子结合成碎片；液相分数判据为，两个熔融金属粒子的液相分数均需满足α1＜α＜α2，其中α1、α2是两个可变参数且0≤α1＜α2≤1，根据具体模拟情况进行调整；粒子间距离判据为，两个熔融金属粒子之间的距离需小于1.2l0；

步骤七，采用流固混合模型(pms方法)修正tn时刻碎片中组成粒子的位置rⁿ和速度uⁿ，分为以下几个步骤：

1，计算碎片的速度和位置

2，计算碎片的转动惯量

3，基于角动量守恒定律，计算碎片的角速度

4，计算tn时刻碎片中粒子的位置rⁿ和速度uⁿ；

步骤八，判断是否达到设定的总的模拟计算时间ttotal，若tn＜ttotal，则跳至步骤二，进行tn+1时刻的计算；若tn＝ttotal，则完成模拟。

本发明还可以包括：

1、步骤二中，焓方法的传热方程为：

2、步骤二中，粒子的液相分数定义如下：

其中下标s和l分别代表固态以及液态。

3、步骤三中，处于相变区域的粒子的密度、热导率根据液相分数加权计算，表示如下：

ρ＝ρs*(1-α)+ρl*α

k＝ks*(1-α)+kl*α。

4、步骤三中，处于相变区域的粒子的粘性系数根据液相分数以指数形式变化，表示如下：

v＝vl·e^5(1-α)。

5、步骤四中，表面张力模型是基于表面自由能，且以粒子间相互作用力的形式呈现，当两个粒子间的距离在(0,re)范围内，则粒子间相互作用力表示如下：

f＝c·(r-l0)·(r-re)/m。

6、步骤六中，碎片结合判据是用来确定破碎凝固产生碎片的粒子组成，包括液相分数判据，以及粒子间距离判据，只有两个判据同时满足才能判定两个粒子结合成碎片；液相分数判据为，两个熔融金属粒子的液相分数均需满足α1＜α＜α2，其中α1、α2是两个可变参数且0≤α1＜α2≤1，根据具体模拟情况进行调整；粒子间距离判据为，相邻的两个金属粒子之间的距离需小于ε·l0，其中l0是粒子间初始距离，ε是一个系数，根据初始布置中一个粒子i和其周围8个粒子之间距离的平均值，推算出ε≈1.2。

7、步骤七中，流固混合模型是基于角动量守恒定律，对尺寸不断变化的碎片的位置和速度进行计算。

本发明针对堆芯熔融金属破碎凝固的模拟方法的不足，提供了一种可实现两种流体传热交混破碎相变过程，尤其针对两相流动、凝固传热以及流固混合等复杂问题的模拟方法。本发明基于无网格移动粒子半隐式法提出了一种针对两种流体传热交混破碎相变过程的无网格模拟方法(multiphasemovingparticlesemi-implicitmethodforfragmentationandphasetransitionsimulation,mmps-fpt)。

采用基于焓方法的凝固相变模型，对相变问题进行了简化；物性参数变化模型对熔融物在凝固破碎过程中的传热及流动物性参数进行修正，使计算更加准确；采用粒子间相互作用的表面张力模型，能够更简单地与无网格法进行结合；采用合理的碎片结合判据确定破碎产生碎片的粒子组成；采用流固混合模型对凝固后碎片的运动进行模拟。

本方法基于移动粒子半隐式法(mps)，提出一种改进的无网格方法，综合各个模型可实现两种流体传热交混破碎相变过程，尤其针对两相流动、凝固传热以及流固混合等复杂问题的模拟计算求解。

附图说明

图1本发明两种流体传热交混破碎相变过程的无网格模拟方法流程图；

图2本发明金属粒子温度、焓值及液相分数之间对应关系图。

图3本发明中碎片结合判据中的粒子间距离判据原理。

具体实施方式

下面结合附图对本发明及其实施方式作进一步详细描述。

如图1所示,本发明实现如下：

步骤一，对于一个特定工况的模拟计算，设定计算中使用的时间步长δt、粒子间初始距离l0、粒子的控制半径re以及总的模拟计算时间ttotal。规定tn＝n·δt，则在t0时刻输入初始化参数，包括熔融金属流体的初始温度th、速度v0、形态，冷却流体的初始温度tc。

步骤二：在tn时刻(n＞0)采用基于焓方法的凝固相变模型进行传热计算。根据tn-1时刻的粒子温度t^n-1、焓值h^n-1、位置r^n-1、密度ρ^n-1以及热导率k^n-1等，计算得到tn时刻粒子的温度tⁿ、焓值hⁿ以及液相分数αⁿ。具体计算表示如下：

其中cps,cpl,l,tmp分别代表固态定压比热、液态定压比热、潜热以及熔点温度，下标s和l分别代表固态以及液态。

图2所示为本发明金属粒子温度、焓值及液相分数之间对应关系图。焓值和温度的对应关系表示如下：

步骤三，采用物性参数变化模型对处于相变区域粒子的物性参数进行修正，根据粒子的液相分数αⁿ，更新得到tn时刻粒子的密度ρⁿ、热导率kⁿ以及粘性系数vⁿ。表示如下：

ρⁿ＝ρs*(1-αⁿ)+ρl*αⁿ

kⁿ＝ks*(1-αⁿ)+kl*αⁿ

步骤四，采用基于表面自由能的表面张力模型计算tn时刻粒子所受的表面张力以粒子间相互作用的形式呈现。当两个粒子间的距离在(0,re)范围内，则粒子间相互作用了表示如下：

其中c为系数，r,m,σ分别为粒子间距离、粒子的质量和表面张力系数。

步骤五，根据步骤三和步骤四所得的粒子物性参数ρⁿ、vⁿ和表面张力采用原始mps方法计算出tn时刻所有粒子的位置rⁿ和速度uⁿ。可分为以下几个小步骤：1，不考虑压力梯度项，显式求解动量方程中的粘性项表面张力项以及重力项g，计算得出tn时刻粒子的中间速度u^*和中间位置r^*；2，计算压力泊松方程，得出tn时刻粒子的压力值pⁿ；3，利用粒子的压力值pⁿ，对粒子速度进行修正，得到tn时刻粒子的速度uⁿ和位置rⁿ。具体计算流程如下：

r^*＝r^n-1+δt·u^*

rⁿ＝r^n-1+δt·uⁿ

步骤六，根据步骤五中计算所得粒子位置rⁿ，采用碎片结合判据确定破碎相变产生碎片的粒子组成。具体包括液相分数判据和粒子间距离判据，只有两个判据同时满足才能判定两个粒子结合成碎片。液相分数判据为，两个熔融金属粒子的液相分数均需满足α1＜α＜α2，其中α1、α2是两个可变参数且0≤α1＜α2≤1，可根据具体模拟情况进行调整。图3所示为本发明中碎片结合判据中的粒子间距离判据原理，相邻的两个金属粒子之间的距离需小于ε·l0。其中l0是粒子间初始距离，ε是一个系数。根据初始布置中一个粒子i和其周围8个粒子之间距离的平均值，可以推算出

步骤七，采用流固混合模型(pms方法)修正tn时刻碎片中组成粒子的位置rⁿ和速度uⁿ。可分为以下几个步骤：1，计算碎片的速度和位置2，计算碎片的转动惯量3，基于角动量守恒定律，计算碎片的角速度4，计算tn时刻碎片中粒子的位置rⁿ和速度uⁿ。对于一个包含n个粒子的碎片，具体计算表示如下：

其中，下标f代表碎片，下标i代表粒子编号，m为旋转矩阵。

步骤八，判断是否达到设定的总的模拟计算时间ttotal。若tn＜ttotal，则跳至步骤二，进行tn+1时刻的计算；若tn＝ttotal，则完成计算。

需要说明的是，以上介绍的是本发明的实施方案而并非限制。任何对本发明技术法案的修改或者等同替代都不脱离本发明技术法案的精神和范围，均应该涵盖在本发明的权利要求范围之内。