CN111274744B - 一种粒子法模拟mcci中气泡夹带现象的边界处理方法 - Google Patents

Abstract

一种粒子法模拟MCCI中气泡夹带现象的边界处理方法，步骤如下：1、建立几何模型，并输入所需的初始参数；2、初始化边界粒子所对应的密度信息；3、建立每个粒子的邻居粒子域；4、根据边界粒子的邻居粒子域内的流体粒子计算边界粒子密度；5、根据边界粒子密度进行压力计算；6、求解粒子的压力梯度，并更新流体粒子的速度及位置信息；7、根据计算时长是否达到设置的计算时长判断是否结束计算，若不结束计算，则重复步骤3‑7。计算流程如图1所示。本发明方法充分考虑了粒子法在模拟MCCI过程中熔融物界面波动对控制域边界压力计算稳定性的影响。通过该方法，能够提高对MCCI过程中熔融物界面处质量与能量交换的认识。

Description

一种粒子法模拟MCCI中气泡夹带现象的边界处理方法

技术领域

本发明涉及利用粒子法对核反应堆严重事故后熔化堆芯与混凝土相互作用(Molten-Core Concrete Interaction,MCCI)过程中混凝土熔化所产生的气泡在反应堆堆芯熔融物分层界面上夹带现象的模拟，具体涉及一种粒子法模拟MCCI中气泡夹带现象的边界处理方法，以提高粒子法对MCCI过程中气泡在熔融物分层界面处夹带现象模拟的准确性与可靠性，提高对MCCI过程中熔融物分层界面上质量与能量交换现象的模拟分析能力。

背景技术

目前，粒子法在捕捉剧烈变化的自由表明问题上展现了强大的数值模拟能力。在核反应堆严重事故领域，核反应堆堆芯熔融物与混凝土的相互作用(Molten-CoreConcrete Interaction,MCCI)过程中，由于熔化堆芯内含有的熔融金属氧化物与熔融金属的密度不同且不互溶，在重力作用下，会导致熔融物分层现象的发生，使密度大的熔融金属氧化物下沉，密度小的熔融金属上升。在熔融物下层的金属氧化物会加热熔化混凝土，产生大量的气体。这些气体在堆芯熔融物内上升并穿过熔融物分层界面时，会将下层熔融金属氧化物夹带至上层熔融金属，导致分层界面的波动，促进两层熔融物之间的换热。粒子法可以精确捕捉气泡上升过程中熔融物界面处的气泡夹带及界面波动现象，从而对MCCI过程中熔融物界面处的物质与能量交换进行精确模拟与分析。然而，在利用粒子法模拟气泡在熔融物界面处的夹带现象时，由于在相界面处上层熔融金属与下层熔融金属氧化物之间存在着严重的物性不连续性，会导致界面处的压力震荡和压力计算不准确性。目前，在粒子法求解压力及压力梯度的过程中，为提高计算过程中的准确性与稳定性，会使用计算粒子与其邻居粒子域内粒子的算数平均密度取代该粒子本身的固定密度。在边界粒子的密度处理问题上，通常也采用上述方法。对于现在公布的利用粒子法进行的多相流模拟中，边界粒子的应用情况可以分为两种：一种为两相流体密度比不是很大，两相流体均与边界粒子接触；另一种情况为高密度比两相模拟，然而，与边界粒子接触的只是其中一相，另一相并不与之接触。在以上两种情况下，边界粒子的密度对多相流模拟的压力计算结果的影响并不大。然而，在用粒子法模拟MCCI过程中气泡夹带现象时，由于熔融金属与熔融金属氧化物之间存在着大密度比，且均会与边界粒子接触，另外由于气泡夹带现象会导致两相界面的波动，与同一边界粒子直接接触的堆芯熔融物流体粒子在不同时刻并不相同，从而导致固定边界粒子密度的处理方法加剧压力计算的不稳定性与不准确性。而采用与流体粒子相同的密度光滑方法，虽然会减少粒子间的密度比，但由于考虑了边界粒子本身的密度，其与流体粒子之间仍然存在着巨大密度差，仍会对堆芯熔融物流体压力的计算产生不可避免的影响，导致熔融物压力计算的不稳定性与不准确性。

发明内容

为了克服上述现有技术存在的问题，本发明的目的是为了利用粒子法对MCCI过程中气泡夹带现象进行模拟时，避免几何模型边界处边界粒子密度与熔融物流体粒子密度的不连续性问题所导致的压力震荡与压力求解不准确性问题。本发明提供了一种用于提高粒子法模拟MCCI过程中气泡夹带现象压力稳定性的边界处理方法，使粒子法在模拟MCCI过程中气泡夹带现象时能够减小流体粒子压力计算的不稳定性与不准确性，促进对MCCI过程中熔融物分层界面处质量与能量交换现象的数值模拟分析。

为了实现上述目的，本发明采取了以下技术方案予以实施：

一种粒子法模拟MCCI中气泡夹带现象的边界处理方法，步骤如下：

步骤1：针对MCCI中气泡夹带现象建立初始粒子几何模型并确定模拟计算时长；用不同种类的粒子代表不同的物质，粒子直径为l，采用1号粒子模拟熔融物中密度小的熔融物流体，采用2号粒子模拟熔融物中密度大的熔融物流体，采用3号粒子模拟气泡，以下将1号粒子和2号粒子统称为熔融物流体粒子；1号粒子、2号粒子和3号粒子统称为流体粒子；根据建立的初始粒子几何模型，进行边界粒子布置，用4号粒子模拟壁面边界粒子，用5号粒子模拟虚拟边界粒子；壁面边界粒子和流体粒子直接接触，参与压力计算，作用为限制流体粒子的运动范围，为计算对象提供固体边界条件；虚拟边界粒子布置在壁面边界粒子外侧，布置3层，虚拟粒子不参与压力计算，作用为为壁面边界粒子和流体粒子的粒子数密度计算提供补偿，保证壁面边界粒子和流体粒子的粒子数密度计算的准确性，以下将4号粒子和5号粒子统称为边界粒子；

步骤2：每种熔融物流体粒子根据所表示的物质具有对应的密度信息ρ_m，边界粒子所对应的初始密度信息由公式(1)进行计算：

式中：

——边界粒子b的初始密度/kg·m^-3；

b——边界粒子b；

ρ_m——熔融物流体种类m的密度/kg·m^-3；

m——熔融物流体种类m；

s——熔融物流体种类数；

步骤3：建立每个粒子的邻居粒子域，粒子的邻居粒子域范围为以该粒子所在位置为圆心，以r_e为半径的圆所覆盖的区域，覆盖区域内的所有粒子均为该粒子的邻居粒子域内粒子，其中r_e为粒子作用半径，由公式(2)计算：

r_e＝3.1·l 公式(2)

式中：

r_e——粒子作用半径/m；

l——粒子直径/m；

步骤4：为保证边界处密度的连续性，防止发生压力计算的不稳定性，根据边界粒子邻居粒子域内的流体粒子的密度计算边界粒子密度，由公式(3)计算：

式中：

ρ_b——边界粒子b的密度/kg·m^-3；

N_f——边界粒子b邻居粒子域内的流体粒子总数；

ρ_f——流体粒子f的密度/kg·m^-3；

f——流体粒子f；

w_bf——边界粒子b与流体粒子f之间的核函数；

其中任意两粒子i与j之间的核函数w_ij的具体计算如公式(4)：

式中：

w_ij——粒子i与粒子j之间的核函数；

r_ij——粒子i与粒子j间的距离/m，有

——粒子i的位置/m；

i——粒子i；

——粒子j的位置/m；

j——粒子j；

另外，若N_f＝0，则该边界粒子的密度不进行更新，仍为上一时间步该边界粒子的密度；

步骤5：根据计算所得的边界粒子密度，对任意粒子i根据公式(5)进行压力计算：

式中：

n₀——初始粒子数密度，为常数；

N_i——粒子i邻居粒子域内的粒子总数；

ρ_i——粒子i的密度/kg·m^-3；

ρ_j——粒子j的密度/kg·m^-3；

P_ij——粒子i和粒子j的压力差/Pa，有P_ij＝P_i-P_j；

P_i——粒子i的压力/Pa；

P_j——粒子j的压力/Pa；

γ——松弛因子；

Δt——计算时间步长/s；

——粒子i速度的估算值/m·s^-1；

n_i ^*——粒子i的临时粒子数密度；

α——可压缩系数；

——下一时间步粒子i的压力/Pa；

k+1——下一时间步；

若粒子i为虚拟边界粒子，则该粒子的压力P_i＝0；

步骤6：计算压力梯度，通过压力梯度进一步计算得到下一个时间步下流体粒子的速度；边界粒子的速度和位置始终不发生改变；粒子i的压力梯度由公式(6)计算：

式中：

——粒子i的压力梯度/Pa·m^-1；

d——数值模拟维度；

C_ij——无量纲修正矩阵；

P_i,min——粒子i邻居粒子域内与粒子i相同种类编号粒子的最小压力/Pa；

每个流体粒子下一时刻的速度

由公式(7)计算：

式中：

——下一时间步流体粒子f的速度/m·s^-1；

——流体粒子f速度的估算值/m·s^-1；

——流体粒子f速度的修正值/m·s^-1；

其中，每个流体粒子速度的修正值

由公式(8)计算：

式中：

——流体粒子f的压力梯度/m·s^-1；

更新的每个流体粒子下一时间步的位置

则通过公式(9)计算获得：

式中：

——下一时间步流体粒子f的位置/m；

——流体粒子f位置的估算值/m；

——流体粒子f速度的修正值/m·s^-1；

步骤7：判断计算总时长是否达到步骤一设定的模拟计算时长，若达到，则结束计算，否则进行下一时间步的计算，重复步骤3-7。

和现有技术相比较，本发明方法是在求解粒子压力之时，采用特殊的粒子光滑方法对边界粒子密度进行光滑求解，通过只利用边界粒子邻居粒子域内的流体粒子的密度计算求解边界粒子密度，同时不考虑边界粒子本身的密度影响。本发明充分考虑了几何模型边界上边界粒子密度对流体粒子压力计算的影响，从而本发明的方法在粒子法模拟MCCI过程中气泡夹带现象时能够减小流体粒子压力计算的不稳定性与不准确性，促进对MCCI过程中熔融物分层界面处质量与能量交换现象的数值模拟分析。

附图说明

图1是一种粒子法模拟MCCI中气泡夹带现象的边界处理方法的计算流程图；

图2是粒子法模拟MCCI中气泡夹带现象的几何模型示意图；

图3是边界粒子布置示意图；

图4是邻居粒子域示意图；

图5是不同密度双层液池结构示意图。

具体实施方式

下面结合附图具体实施方式对本发明方法作进一步详细说明：

如图1所示，本发明一种粒子法模拟MCCI中气泡夹带现象的边界处理方法，包括如下步骤：

步骤1：针对MCCI中气泡夹带现象建立初始粒子几何模型并确定模拟计算时长，几何模型如图2所示；用不同种类的粒子代表不同的物质，粒子直径为l，采用1号粒子模拟熔融物中密度小的熔融物流体，采用2号粒子模拟熔融物中密度大的熔融物流体，采用3号粒子模拟气泡，以下将1号粒子和2号粒子统称为熔融物流体粒子；1号粒子、2号粒子和3号粒子统称为流体粒子；根据建立的初始粒子几何模型，进行边界粒子布置，用4号粒子模拟壁面边界粒子，用5号粒子模拟虚拟边界粒子；壁面边界粒子和流体粒子直接接触，参与压力计算，作用为限制流体粒子的运动范围，为计算对象提供固体边界条件；虚拟边界粒子布置在壁面边界粒子外侧，布置3层，虚拟粒子不参与压力计算，作用为为壁面边界粒子和流体粒子的粒子数密度计算提供补偿，保证壁面边界粒子和流体粒子的粒子数密度计算的准确性，边界粒子布置如图3所示，以下将4和5号粒子统称为边界粒子；

步骤2：每种熔融物流体粒子根据所表示的物质具有对应的密度信息ρ_m，边界粒子所对应的初始密度信息由公式(1)进行计算：

式中：

——边界粒子b的初始密度/kg·m^-3；

b——边界粒子b；

ρ_m——熔融物流体种类m的密度/kg·m^-3；

m——熔融物流体种类m；

s——熔融物流体种类数；

步骤3：建立每个粒子的邻居粒子域，粒子的邻居粒子域范围为以该粒子所在位置为圆心，以r_e为半径的圆所覆盖的区域，覆盖区域内的所有粒子均为该粒子的邻居粒子域内粒子，如图4所示，其中r_e为粒子作用半径，由公式(2)计算：

r_e＝3.1·l 公式(2)

式中：

r_e——粒子作用半径/m；

l——粒子直径/m；

步骤4：为保证边界处密度的连续性，防止发生压力计算的不稳定性，根据边界粒子邻居粒子域内的流体粒子的密度计算边界粒子密度，由公式(3)计算：

式中：

ρ_b——边界粒子b的密度/kg·m^-3；

N_f——边界粒子b邻居粒子域内的流体粒子总数；

ρ_f——流体粒子f的密度/kg·m^-3；

f——流体粒子f；

w_bf——边界粒子b与流体粒子f之间的核函数；

其中任意两粒子i与j之间的核函数w_ij的具体计算如公式(4)：

式中：

w_ij——粒子i与粒子j之间的核函数；

r_ij——粒子i与粒子j间的距离/m，有

——粒子i的位置/m；

i——粒子i；

——粒子j的位置/m；

j——粒子j；

另外，若N_f＝0，则该边界粒子的密度不进行更新，仍为上一时间步该边界粒子的密度；

步骤5：根据计算所得的边界粒子密度，对任意粒子i根据公式(5)进行压力计算：

式中：

n₀——初始粒子数密度，为常数；

N_i——粒子i邻居粒子域内的粒子总数；

ρ_i——粒子i的密度/kg·m^-3；

ρ_j——粒子j的密度/kg·m^-3；

P_ij——粒子i和粒子j的压力差/Pa，有P_ij＝P_i-P_j；

P_i——粒子i的压力/Pa；

P_j——粒子j的压力/Pa；

γ——松弛因子；

Δt——计算时间步长/s；

——粒子i速度的估算值/m·s^-1；

n_i ^*——粒子i的临时粒子数密度；

α——可压缩系数；

——下一时间步粒子i的压力/Pa；

k+1——下一时间步；

若粒子i为虚拟边界粒子，则该粒子的压力P_i＝0；

步骤6算压力梯度，通过压力梯度进一步计算得到下一个时间步下流体粒子的速度；边界粒子的速度和位置始终不发生改变；粒子i的压力梯度由公式(6)计算：

式中：

——粒子i的压力梯度/Pa·m^-1；

d——数值模拟维度；

C_ij——无量纲修正矩阵；

P_i,min——粒子i邻居粒子域内与粒子i相同种类编号粒子的最小压力/Pa；

每个流体粒子下一时刻的速度

由公式(7)计算：

式中：

——下一时间步流体粒子f的速度/m·s^-1；

——流体粒子f速度的估算值/m·s^-1；

——流体粒子f速度的修正值/m·s^-1；

其中，每个流体粒子速度的修正值

由公式(8)计算：

式中：

——流体粒子f的压力梯度/m·s^-1；

更新的每个流体粒子下一时间步的位置

则通过公式(9)计算获得：

式中：

——下一时间步流体粒子f的位置/m；

——流体粒子f位置的估算值/m；

——流体粒子f速度的修正值/m·s^-1；

步骤7：判断计算总时长是否达到步骤一设定的模拟计算时长，若达到，则结束计算，否则进行下一时间步的计算，重复步骤3-7。

通过以上步骤使粒子法在模拟MCCI中气泡夹带现象的过程中对几何模型的边界粒子所对应的密度信息进行了时时更新，实现了用粒子法对MCCI过程中气泡在堆芯熔融物分层界面处夹带现象的精确模拟，解决了由于堆芯熔融物分层界面波动所导致的边界处压力不稳定性与不准确性问题，促进了粒子法对核反应堆严重事故后熔化堆芯与混凝土相互作用(Molten-Core Concrete Interaction,MCCI)过程中堆芯熔融物界面处的物质与能量交换的数值模拟与分析。

下面以不同密度双层液池压力计算结果为例说明本发明的效果：

图5是不同密度双层液池结构示意图。由于在MCCI过程中，反应堆堆芯熔融物内既含有密度大的熔融金属氧化物，又含有密度小的熔融金属，熔融金属氧化物逐渐下沉而熔融金属逐渐上升，最终形成类似于图5所示的分层现象。现以图5为例进行结果分析。具体结构如下：上层为0.2m高的气体，密度为1.18kg/m³；下层为0.2m深的水，密度为997kg/m³，其密度比接近1000，远远超过MCCI过程中熔融金属氧化物与熔融金属的密度比。重力加速度为9.81m/s²。监视下层水最右下角粒子的计算压力变化。在该点，理论压力值为：195.84Pa。

表1是两种方法下随机选取的20个相同时刻计算获得的压力的误差平均值。可以从计算结果中看出，未采用本专利中提出的计算方法获得的压力平均误差为12.76％；而采用本专利中提出的方法，将压力计算平均误差缩小到了7.33％，准确性较原方法有明显的改进。表2是两种方法下一段时间内计算压力最大误差情况。如表2所示，采用本发明计算获得的压力的最大误差为9.30％，较未采用本专利方法中的最大误差38.39％有了显著的改善，充分说明了该发明在计算过程中能够显著提高压力计算的稳定性。

表1不同密度双层液池压力平均误差计算结果

压力计算	平均误差/％
		未采用本发明	12.76
采用本方法	7.33

表2不同密度双层液池压力最大误差计算结果

压力计算	最大误差/％
		未采用本发明	13.65
采用本发明	9.30

针对数据验证结果，本发明一种粒子法模拟MCCI中气泡夹带现象的边界处理方法能够有效提高粒子法在模拟MCCI过程中融化混凝土产生的气泡在堆芯熔融物分层界面处夹带现象压力稳定性的计算，有效解决了由于界面波动所导致的几何模型边界粒子与熔融物不同密度流体粒子之间的密度不连续性问题，从而提高了模拟过程中流体粒子压力计算的稳定性与准确性，实现了对MCCI过程中堆芯熔融物分层界面上的质量与能量交换的分析。

Claims (1)

1.一种粒子法模拟MCCI中气泡夹带现象的边界处理方法，其特征在于：该方法包括以下步骤：

步骤1：针对MCCI中气泡夹带现象建立初始粒子几何模型并确定模拟计算时长；用不同种类的粒子代表不同的物质，粒子直径为l，采用1号粒子模拟熔融物中密度小的熔融物流体，采用2号粒子模拟熔融物中密度大的熔融物流体，采用3号粒子模拟气泡，以下将1号粒子和2号粒子统称为熔融物流体粒子；1号粒子、2号粒子和3号粒子统称为流体粒子；根据建立的初始粒子几何模型，进行边界粒子布置，用4号粒子模拟壁面边界粒子，用5号粒子模拟虚拟边界粒子；壁面边界粒子和流体粒子直接接触，参与压力计算，作用为限制流体粒子的运动范围，为计算对象提供固体边界条件；虚拟边界粒子布置在壁面边界粒子外侧，布置3层，虚拟边界粒子不参与压力计算，作用为为壁面边界粒子和流体粒子的粒子数密度计算提供补偿，保证壁面边界粒子和流体粒子的粒子数密度计算的准确性，以下将4号粒子和5号粒子统称为边界粒子；

步骤2：每种熔融物流体粒子根据所表示的物质具有对应的密度信息ρ_m，边界粒子所对应的初始密度信息由公式(1)进行计算：

式中：

——边界粒子b的初始密度，单位为：kg·m^-3；

b——边界粒子b；

ρ_m——熔融物流体种类m的密度，单位为：kg·m^-3；

m——熔融物流体种类；

s——熔融物流体种类数；

步骤3：建立每个粒子的邻居粒子域，粒子的邻居粒子域范围为以该粒子所在位置为圆心，以r_e为半径的圆所覆盖的区域，覆盖区域内的所有粒子均为该粒子的邻居粒子域内粒子，其中r_e为粒子作用半径，由公式(2)计算：

r_e＝3.1·l 公式(2)

式中：

r_e——粒子作用半径，单位为：m；

l——粒子直径，单位为：m；

步骤4：为保证边界处密度的连续性，防止发生压力计算的不稳定性，根据边界粒子邻居粒子域内的流体粒子的密度计算边界粒子密度，由公式(3)计算：

式中：

ρ_b——边界粒子b的密度，单位为：kg·m^-3；

N_f——边界粒子b邻居粒子域内的流体粒子总数；

ρ_f——流体粒子f的密度，单位为：kg·m^-3；

f——流体粒子f；

w_bf——边界粒子b与流体粒子f之间的核函数；

其中任意两粒子i与j之间的核函数w_ij的具体计算根据公式(4)进行：

式中：

w_ij——粒子i与粒子j之间的核函数；

r_ij——粒子i与粒子j间的距离，单位为：m，有

——粒子i的位置；

i——粒子i；

——粒子j的位置；

j——粒子j；

另外，若边界粒子对应的N_f＝0，则该边界粒子的密度不进行更新，仍为上一时间步该边界粒子的密度；

步骤5：根据计算所得的边界粒子密度，对任意粒子i根据公式(5)进行压力计算：

式中：

n₀——初始粒子数密度，为常数；

ρ_i——粒子i的密度，单位为：kg·m^-3；

ρ_j——粒子j的密度，单位为：kg·m^-3；

P_ij——粒子i和粒子j的压力差，单位为：Pa，有P_ij＝P_i-P_j；

P_i——粒子i的压力，单位为：Pa；

P_j——粒子j的压力，单位为：Pa；

γ——松弛因子；

Δt——计算时间步长，单位为：s；

——粒子i速度的估算值，单位为：m·s^-1；

——粒子i的临时粒子数密度；

α——可压缩系数；

——下一时间步粒子i的压力，单位为：Pa；

k+1——下一时间步；

若粒子i为虚拟边界粒子，则该粒子的压力P_i＝0；

步骤6：计算压力梯度，通过压力梯度进一步计算得到下一个时间步下流体粒子的速度；边界粒子的速度和位置始终不发生改变；粒子i的压力梯度由公式(6)计算：

式中：

——粒子i的压力梯度，单位为：Pa·m^-1；

d——数值模拟维度；

C_ij——无量纲修正矩阵；

P_i,min——粒子i邻居粒子域内与粒子i相同种类编号粒子的最小压力，单位为：Pa；

每个流体粒子下一时刻的速度

由公式(7)计算：

式中：

——下一时间步流体粒子f的速度，单位为：m·s^-1；

——流体粒子f速度的估算值，单位为：m·s^-1；

——流体粒子f速度的修正值，单位为：m·s^-1；

其中，每个流体粒子速度的修正值

由公式(8)计算：

式中：

——流体粒子f的压力梯度，单位为：m·s^-1；

更新的每个流体粒子下一时间步的位置

则通过公式(9)计算获得：

式中：

——下一时间步流体粒子f的位置，单位为：m；

——流体粒子f位置的估算值，单位为：m；

——流体粒子f速度的修正值，单位为：m·s^-1；

步骤7：判断计算总时长是否达到步骤1设定的模拟计算时长，若达到，则结束计算，否则进行下一时间步的计算，重复步骤3-7。

Images