[发明专利]一种粒子法模拟MCCI中气泡夹带现象的边界处理方法

权利要求书

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。