采用多边形壁面的MPS粒子分裂算法研究

张 静，熊进标

(上海交通大学 核科学与工程学院，上海 200240)

反应堆严重事故模拟中涉及大量含有自由界面或组分界面的多相流动，不具有拓扑结构的粒子法在这类流动的模拟中有其独特的优势。李勇霖等[1]基于移动粒子半隐式(MPS)法[2]模拟堆芯材料在高温熔化过程中的共晶反应现象。张明昊等[3]采用MPS-MAFL方法针对入口流量脉动条件下的单气泡垂直上升运动行为进行了特性分析。Zhu等[4-6]引入固体被动漂移模型、势能界面张力等模型验证MPS法对熔融现象计算的可靠性。但由于传统的MPS法中粒子尺寸相等，无法对局部区域(如相变界面附近)进行加密计算。因此提出多粒径模型，Tang等[7]针对不同分辨率界面上的粒子，对梯度算子模型和拉普拉斯算子模型添加修正系数来满足守恒性。Chen等[8]采用三次样条函数代替原先的核函数，添加额外的权重以及采取5步分裂方法来保证分裂计算的准确性和压力的稳定性。Shibata等[9]设置局部重叠的子区域，并添加区域边界条件来实现局部的精细化计算。

此外，为了改善壁面粒子对计算效率和准确性产生的影响，人们提出了多边形壁面模型。与传统的采用固定流体粒子作为壁面的方式相比，该模型不设置壁面粒子，因此在求解压力泊松方程时仅包含流体粒子，提高迭代计算效率；同时，壁面的压力采取合理的方式替代求解。对于壁面附近粒子数密度以及压力的计算方法，Harada等[10]基于粒子在壁面压力梯度的作用下回到平衡位置的假设提出多边形壁面形式。Zhang等[11]在Harada的方法上进行改进，提出一种新的虚拟壁面粒子生成方式，改进适用于多边形壁面的压力源项[12]。Mitsume等[13]通过计算镜面粒子的压力梯度代替实际壁面梯度。

在目前的研究进展中，粒子分裂一般需要设置较为复杂的边界条件，采取多边形壁面模型的算法中往往设置压力显示计算的方式，鲜有粒子分裂模型与多边形壁面模型相结合的研究。本文针对MPS进行算法改进，使其能适用于多边形壁面条件下基于压力隐式求解的粒子分裂计算，对比溃坝实验验证改进后模型的准确性和高效性，为进一步计算三维相变传热模型奠定基础。

1 移动粒子半隐式法

对于不涉及传热的流动过程，MPS的控制方程包括连续性方程和动量方程：

(1)

(2)

式中：ρ为密度；t为时间；u为速度；P为压力；ν为运动黏性；g为重力加速度；f为重力以外的其他外力。

MPS法将求解域离散为一系列的粒子，通过粒子间相互作用模型离散控制方程，基于追踪粒子得到计算域的信息。粒子相互作用模型包括核函数、梯度算子模型和拉普拉斯算子模型3个基本模型，常用的核函数模型为：

re=RDi

(3)

式中：|rij|为粒子i与j的间距；re为有效半径；R为有效半径系数，通常取2.1～4.1；Di为粒子直径。

对有效半径内粒子的核函数求和，得到目标粒子的粒子数密度ni为：

ni=∑wi(|rij|)

(4)

梯度算子模型为：

(5)

拉普拉斯算子模型为：

(6)

通过联立动量和连续性方程，可得压力泊松方程为：

(7)

在压力计算中，表面粒子的压力设置为0。通常采用粒子数密度阈值来判定表面粒子[14]，即：

〈n〉i<αn0

(8)

MPS算法流程如图1所示。在1个时间步长内，首先利用黏性力与重力显式更新速度和位置，然后通过式(7)隐式计算压力，并利用式(5)得到压力梯度，再次更新粒子的速度及位置信息，确认达到收敛标准后进入下一时间步的计算。

图1 MPS算法流程

2 模型开发

在多粒径计算中，直径不同的粒子质量存在差异，传统MPS法的粒子数密度计算方法无法准确反映粒子聚集程度(即密度变化)，从而导致压力计算不准确。同时大量的壁面粒子也增大了泊松方程的求解量。本文针对上述问题进行改进，开发适用于多边形壁面的粒子分裂模型。

2.1 有效半径

当计算域内的粒子尺寸不等时，采用传统的有效半径计算方法会影响计算的准确性。本研究中仅考虑两种粒径的粒子计算，当无外力作用下粒子均匀排布时，如图2a所示，粒径不同的i粒子与j粒子位于分辨率界面上，实线圈为i粒子的作用域，虚线圈为j粒子的作用域。理论上两者的粒子数密度应相等，粒子间保持相对静止；而由式(3)、(4)计算得到的i粒子的数密度远大于j粒子的数密度，两个粒子间会生成较大的作用力而运动，因此需要对有效半径进行改进。

本文参考Tanaka等[15]的处理方法，采用平均有效半径代替原固定的有效半径，即：

(9)

式中，R取值为3.1。

改进后的作用域如图2b所示。当大粒子i与小粒子相互作用时，作用域如小半圆实线圈所示，作用范围与改进前相比有所缩小；粒子i与其直径相等的大粒子作用时，作用域如大半圆实线圈所示，与原有效半径相等。同理，小粒子j的作用域也改进为大小半圆组合的不规则范围。改进后的有效半径计算方法减小了粒子i与粒子j间的数密度差值，同时也避免了由于两个粒子的有效半径的不同而可能造成的相互作用力不相等现象。

a——改进前；b——改进后

2.2 粒子数密度

粒子数密度ni包括流体粒子的密度nif及固体壁面的密度niw：

ni=nif+niw

(10)

nif的计算参考多粒径模型的处理方法[7]，根据牛顿第三定律对数密度添加系数项：

(11)

式中，m为粒子质量。

niw的求解与传统MPS法计算不同。首先在初始化过程中，对计算域建立大小均匀的背景网格，网格尺寸与小粒子尺寸一致。计算每个网格点到多边形壁面的最小距离，并在多边形壁面外侧生成虚拟壁面粒子。对位于壁面有效半径内的网格点，计算其曲率系数Cg[12]为：

(12)

式中：Zg为虚拟壁面粒子在g网格点处的数密度；Wg为同等距离下，固体粒子组成平板壁面时对应的g网格处的壁面数密度。当Cg=1时，在g网格点有效半径内壁面为平板；当Cg≠1时，该网格点有效半径内存在壁面转角。

在二维计算中，搜索流体粒子周围的4个网格点，并基于网格点的壁面距离插值得到流体粒子的壁面距离riw，进而计算该距离对应的平板壁面数密度Wiw；对网格的曲率系数进行插值得到粒子在该位置处的曲率系数Ciw，将平板壁面数密度Wiw与曲率系数Ciw相乘，得到流体粒子的实际壁面数密度niw。壁面数密度计算示意图如图3所示，红点为参与计算的网格点，绿线为多边形壁面，虚线粒子为虚拟壁面粒子。

图3 壁面数密度计算示意图

2.3 黏性拉普拉斯算子模型

(13)

(14)

(15)

式中，uw为壁面速度。

2.4 压力泊松方程

采用多边形壁面模型后，无需设置壁面粒子。因此在求解压力泊松方程时，仅考虑流体粒子间的相互作用。多粒径的压力拉普拉斯算子模型的计算采取式(16)代替式(6)：

(16)

参考Zhang等[12]对适用于多边形壁面的压力源项的处理，添加粒子流体数密度与粒子总数密度的比值这一系数，来反映壁面压力对流体的影响。本文将其改进为适用于分裂计算的表达形式：

(17)

式中：γ为修正系数，取0.015；k为当前时间步，k+1为下一时间步；u*为两个时间步间的临时速度。联立式(16)与(17)，隐式求解得到流体粒子的压力。

2.5 梯度算子模型

(18)

(19)

图4 镜面粒子压力梯度计算[16]

(20)

2.6 自由表面粒子识别

为了提高自由表面粒子识别的准确性，本文参考Shibata等[17]提出的基于弧度的表面识别方法进行模型改进。模型中，将目标粒子i假想为一个点光源，照亮半径为2.1di的圆周，计算有效半径内的相邻粒子由于遮挡光线而在圆周上形成的阴影弧长。定义表面率A为：

(21)

由于采用式(21)进行表面识别需要对所有粒子进行判断，计算量较大，本文采用式(8)与式(21)相结合的方式。首先对所有粒子进行基于粒子数密度的筛选，α取0.97，将符合条件的粒子再进行基于弧度的识别判断，从而在保证精度的同时提高效率。

2.7 分裂模型

粒子分裂计算的流程如图5所示。

图5 粒子分裂计算流程图

首先设置分裂条件，选取粒子位置为判断标准，在计算域内划分出需要精细计算的分裂区域。识别出符合分裂条件的大粒子，对分裂生成新的小粒子进行位置赋值。本文仅考虑1个大粒子均匀分裂成4个小粒子的情况，其排列分布如图6所示。为了保证分裂后粒子速度与周围粒子光滑过渡，本文参考Tanaka等[15]提出的梯度修正法，给定新粒子的速度：

图6 粒子分裂后位置及大小示意图

(22)

式中：ui,M为分裂后粒子速度；ui为分裂前速度；ri,M为分裂后粒子位置；ri为分裂前位置；M为分裂后粒子序号，M=1、2、3、4。

在对体积等常量进行赋值后，将发生分裂的大粒子的类型更新为Ghost，不参与后续的计算。分裂循环完成后，更新全计算域粒子的数密度，结束分裂程序的计算。

3 溃坝模拟

3.1 静压测试

静压计算中理论上粒子处于静止状态，监测点的压力具有精确的理论值，因此采用静压计算来评价改进后MPS模型的准确性以及稳定性。液体几何长度为0.36 m、高度为0.48 m，坐标原点位于左下角，监测(0.2 m，0.01 m)点处的压力。计算多边形壁面条件下，粒子直径分别为0.01 m及0.005 m工况(背景网格尺寸均为0.005 m)下的压力波动，如图7所示。计算达到稳定后，两个工况下压力均波动较小且略低于理论值，相对误差约为6.25%。因此，静压测试可验证多边形壁面模型的准确性和稳定性。

图7 不同工况在监测点的压力

3.2 无挡板溃坝

计算Martin等[18]的无挡板溃坝实验，验证粒子分裂模型的可靠性。无挡板溃坝模型如图8所示，水柱高0.44 m、长0.22 m。设置粒子的密度为1 000 kg/m3，运动黏性为1.0 mm2/s，重力加速度为9.81 m/s2，时间步长为1×10-4s，共计算了3种工况：1)粒子直径为0.01 m；2)粒子直径为0.005 m；3)粒子在波前距离大于0.4 m时分裂，分裂前直径为0.01 m、分裂后直径为0.005 m，分裂后粒子排布形式如图6所示。

图8 无挡板溃坝模型

图9 无量纲时间变化图

3.3 有挡板溃坝

计算有档板的溃坝模型[19]，以验证压力计算的准确性以及分裂模型的计算效率。有挡板溃坝模型如图10所示。水柱高0.12 m、长0.68 m，容器长1.18 m，压力测试点位于右侧壁面距离地面0.01 m处，即A点所示。

图10 有挡板溃坝模型

计算了3种工况：1)壁面设置为传统固定粒子方式，粒子直径为0.01 m；2)壁面设置为传统固定粒子方式，粒子直径为0.005 m；3)壁面设置为多边形壁面形式，当粒子的横坐标(长度方向)大于0.7 m时分裂，分裂前粒子直径为0.01 m，分裂后粒子直径为0.005 m，分裂后粒子排布形式如图6所示。设置时间步长为5×10-4s，选取0.5～1.0 s共6个时刻的压力计算结果进行比较，如图11所示。可看出，不同时刻的3种工况下的粒子流动具有相似性。对于自由表面的捕捉，采用均匀小粒子计算以及分裂计算的工况(图11b、c)能得到清晰的液面，而均匀大粒子的工况(图11a)得到的结果较为粗糙。

a——粒子壁面，均匀粒径0.01 m；b——粒子壁面，均匀粒径0.005 m；c——多边形壁面，初始粒径0.01 m，分裂后0.005 m

在计算过程中对图10中的A点进行压力检测，并与实验值进行比较(图12)。无分裂+粒子壁面工况的结果为蓝色线所示，撞击壁面时间、撞击后监测点的压力均与实验值吻合良好；分裂+粒子壁面工况的结果如红色线所示，除压力次高峰与实验值相比滞后约0.1 s外，其余均与实验值吻合良好；分裂+多边形壁面工况的结果为绿线所示，撞击壁面时间与实验值相近，但是检测点的压力波动幅值较大，且压力偏大。因此改进后模型的压力稳定性需要进一步的研究和改善。

图12 A点压力比较

对改进后模型的计算效率进行分析，计算均采用i9-10900处理器，内核频率为2.8 GHz。对数密度计算模块、压力求解模块以及计算总时长进行比较(表1)。其中，粒子个数为整个计算过程中的平均粒子个数，模块计算时长为1个时间步内、单次计算所需要的平均时间；总时长占比为改进后的模型与传统MPS模型(模型1)总时长的比值。由表1可看出，若仅采用分裂计算(模型1与2比较)，粒径增大、粒子个数减少使得数密度及压力求解模块的时长减小，总计算时长为原来的44.28%；若仅采用多边形壁面(模型1与3比较)，参与压力泊松求解的粒子数减少、迭代效率提高，总时长占比为37.42%；若采用多边形壁面与分裂模型相结合的形式(模型1与4比较)，所需粒子个数进一步减少，平均仅需1 486.3个，总时长占比为18.75%，略高于两个模型单独作用的时长占比的乘积(16.57%)。因此，采用多边形壁面的粒子分裂模型可大幅提高计算的效率。

表1 不同工况计算时长比较

4 结论

本文基于MPS法，开发适用于多边形壁面的粒子分裂模型，将改进后的模型进行无挡板以及含有挡板的溃坝实验计算，得出以下结论：1)改进后的模型流动结果与实验值吻合良好；2)改进后的模型压力波动较大且数值偏高，需要对压力稳定性进行改善；3)采用多边形壁面的粒子分裂模型可大幅提高MPS法的计算效率，改进后模型计算总时长与原模型的占比为18.75%。