一种基于MPM和相场模型的多相流交互模拟算法

一种基于mpm和相场模型的多相流交互模拟算法
技术领域
1.本发明涉及计算机流体仿真领域，尤其是一种基于mpm和相场模型的多相流交互模拟算法。


背景技术：

2.基于物理的仿真方法是计算机图形学的研究热点方向，被广泛应用于影视动画、虚拟现实等领域。其中，多相流的交互是非常普遍的现象。针对不同的流体在交互过程中会产生不同的效果。如：油滴入水中，形成明显的分层效果，墨水滴入水中，会逐渐溶解扩散。更为复杂的交互，如萃取，部分溶解等。目前多相流模拟的方法大致可以分为两大类：拉格朗日法与欧拉网格法。
3.拉格朗日法侧重于捕捉粒子细节，大多扩展于粒子解算器中(如sph、pbf)。文献[ren b,li c,yan x,et al.multiple-fluid sph simulation using a mixture model[j].acm transactions on graphics(tog),2014,33(5):1-11.]使用体积分数表示不同相或组分的分布，引入漂移速度的概念处理组分之间的扩散。但是该文献方法未考虑能量，无法模拟萃取、部分溶解等较复杂的多相流交互现象。文献[yang t,chang j,ren b,et al.fast multiple-fluid simulation using helmholtz free energy[j].acm transactions on graphics(tog),2015,34(6):1-11.]首次将helmholtz自由能引入多相流交互系统中，扩展了多流体交互的模拟。但是基于拉格朗日的求解器引入了时间步限制，而且忽略了气液界面的模拟。文献[ren b,xu b,li c.unified particle system for multiple-fluid flow and porous material[j].acm transactions on graphics(tog),2021,40(4):1-14.]引入了一种新的虚相位概念，提出了一个统一的粒子框架模拟多种流体流动和多孔材料。
[0004]
欧拉网格侧重于对不同阶段的界面的处理，通过跟踪界面来获得合理的效果。文献[kang n,park j,noh j,et al.a hybrid approach to multiple fluid simulation using volume fractions[c]//computer graphics forum.oxford,uk:blackwell publishing ltd,2010,29(2):685-694.]将level set和体积分数结合起来，不仅捕捉了非混相流体之间的不连续界面，而且捕捉了混相流体之间的平滑过渡。文献[nielsen m b,o.a two-continua approach to eulerian simulation of water spray[j].acm transactions on graphics(tog),2013,32(4):1-10.]通过将欧拉空间的每个节点等效为气相与液相双向耦合的连续体，模拟出真实的水雾现象。
[0005]
作为混合欧拉-拉格朗日求解器的物质点法在多相流交互的模拟中也体现出了很好的可扩展性，文献[yan x,li c f,chen x s,et al.mpm simulation of interacting fluids and solids[c]//computer graphics forum.2018,37(8):183-193.]提出了一种模拟气相与固相的统一mpm(material point method，物质点法)框架。但是该文献并没有处理复杂的多相流交互，如萃取。
[0006]
综上，现阶段针对多相流交互的模拟依然存在以下问题：
[0007]
(1)引入能量之后，使得基于拉格朗日的统一求解器受到了时间步的限制，导致相扩散非常缓慢，甚至需要增加外力。
[0008]
(2)目前的混相、非混相统一处理的求解中，大多忽略了气液界面，降低了模拟的真实性。


技术实现要素：

[0009]
本发明需要解决的技术问题是提供一种基于mpm和相场模型的多相流交互模拟算法，能够解决上述背景技术中提到的技术问题。
[0010]
为解决上述技术问题，本发明所采用的技术方案是：一种基于mpm和相场模型的多相流交互模拟算法，包括如下步骤：
[0011]
步骤s1、将存储在拉格朗日粒子的质量参数、相场参数、动量参数分别转移到欧拉网格；
[0012]
步骤s2、根据化学势控制整个系统的能量，对欧拉网格相场进行更新；
[0013]
步骤s3：根据重力、弹性力和更新后的等效相场形式的表面张力对欧拉网格动量进行更新；
[0014]
步骤s4：更新完网格动量之后，在欧拉网格上对边界采用滑移边界条件进行处理；
[0015]
步骤s5、更新完相场与动量之后，将动量参数和相场参数从欧拉网格变量转换到拉格朗日变量；
[0016]
步骤s6、更新拉格朗日粒子状态：根据传回拉格朗日粒子的信息，对拉格朗日粒子的位置、质量状态进行更新。
[0017]
本发明技术方案的进一步改进在于：所述步骤s1的具体过程为：
[0018]
步骤s11、将存储在拉格朗日粒子的质量参数转移到欧拉网格：将拉格朗日质量转换为欧拉网格质量，使用apic作为传递方式，该过程满足以下公式：
[0019]
mi＝∑
pmp
ω
ip
，
[0020]
其中，mi表示欧拉节点i的质量，m
p
表示拉格朗日粒子p的质量,ω
ip
是拉格朗日粒子和欧拉网格映射的具有连续性的二次b-样条插值函数中的权重；
[0021]
步骤s12、将存储在拉格朗日粒子的动量参数转移到欧拉网格：转移的过程满足动量守恒，使用apic作为传递方式，该过程满足以下公式：
[0022]
(mv)i＝∑
pmp
ω
ip
(v
p
a
p
(x
i-x
p
))，
[0023]
其中，a
p
为拉格朗日粒子p的仿射速度，v
p
表示拉格朗日粒子p的体积，xi为欧拉节点i的位置，x
p
为拉格朗日粒子p的位置；
[0024]
步骤s13、将存储在拉格朗日粒子的相场参数转移到欧拉网格：使用质量分数作为相场参数反映多相流体的相容性，对于n相交互系统，只需要求解n-1相场方程，同时转移过程也只需要转移n-1个相场参量，选择mpm传递方式，相场参数转移满足以下公式：
[0025]ci
＝∑
pmpcp
ω
ip
，
[0026]
其中，ci表示欧拉网格的相场参数，c
p
表示拉格朗日粒子的相场参数。
[0027]
本发明技术方案的进一步改进在于：所述步骤s2中根据化学势更新欧拉网格相场的过程为：
[0028]
步骤s21、求解化学势：化学势表征相组分之间的转移趋势，控制整个系统由高化
学势向低化学势发展，将相场参数插值到欧拉网格之后，采用中心网格、中心差分对其进行离散化；
[0029]
相场的控制方程为：
[0030]
f(c)＝∫
ω
[b(c(x,t)) i(c(x,t))]dx
[0031][0032]
其中，ω表示系统的空间区域，k为常数，用于控制扩散速度，i表示界面能，用于减小质量分数场的梯度，使得表面变化更加平滑，b表示体能，控制流体如何进行交互，c(x,t)表示t时刻，x位置的相场取值；
[0033]
化学势的控制方程为：
[0034][0035]
μ表示化学势，b表示体能，k为常数，c表示相场参数，表示拉普拉斯算子；
[0036]
步骤s22、选择体能函数：用体能函数局部最小值控制流体的低能状态，对流体从高能状态演化为低能状态的过程进行模拟，设计不同的体能函数用以实现不用的交互效果，双相交互使用体能函数图像进行展示，对有效区间的极小值进行标记，三相交互采用吉布斯三角形作为能量的展示，设计体能函数遵循：能量的极小值位于有效区内，以确保最终可以达到能量最小化状态；化学势对于某一相的偏导数只应与该相的质量比有关；
[0037]
步骤s23、更新相场：针对不同的多相流交互设计对应的体能函数之后，使用化学势对相场进行更新，将简并迁移率设为常数，采用中心网格、中心差分对其进行离散化，相场更新满足以下公式：
[0038][0039][0040]
其中，表示物质导数，m表示简并迁移率，μ表示化学势。
[0041]
本发明技术方案的进一步改进在于：所述步骤s3中使用更新的相场，对不同组分之间的界面进行捕捉，施加表面张力，进一步对动量进行更新，表面张力控制方程满足：
[0042][0043][0044][0045]
其中，f
σ
是使用相场形式等效的表面张力，表示l1、l2相之间的表面张力系数，表示l1相的相场参数，表示l2相的相场参数，ε与扩散界面相关的系数。
[0046]
本发明技术方案的进一步改进在于：所述步骤s3中根据重力和弹性力对欧拉网格
动量进行更新遵循如下公式：
[0047][0048]
其中，fe表示弹性力，δx表示欧拉网格的网格间距，为拉格朗日粒子的初始体积，p表示皮奥拉基尔霍夫应力，表示t 1时刻拉格朗日粒子p的形变梯度，表示拉格朗日粒子p的形变梯度矩阵的转置，t表示转置，fg表示重力，g表示重力加速度，fi表示欧拉网格i受力的总和，t表示当前时刻，t 1表示下一时刻，δt表示一个时间步，表示欧拉网格i下一时刻的速度，表示欧拉网格i当前时刻的速度。
[0049]
本发明技术方案的进一步改进在于：所述步骤s4中更新完网格动量之后，在欧拉网格上对边界采用滑移边界条件进行处理，采用滑移边界条件：
[0050][0051]
其中，表示碰撞尚未结束的网格速度，在更新拉格朗日粒子位置之前，多次对拉格朗日粒子对应用碰撞处理，n表示曲面法线。
[0052]
本发明技术方案的进一步改进在于：所述步骤s6的具体过程为：
[0053]
步骤s51、在欧拉网格上更新完动量参数之后，将动量参数从欧拉网格传回拉格朗日粒子：采用apic作为传递方式，动量参数转移满足：
[0054]vp
＝∑iω
ip
vi[0055]
其中，v
p
表示拉格朗日粒子p的速度，vi表示欧拉网格i的速度；
[0056]
步骤s52、在欧拉网格上更新完相场参数之后，将相场参数从欧拉网格传回拉格朗日粒子，转移回拉格朗日粒子的相场信息为每种相组分的质量分数，该过程遵循以下公式：
[0057]cp
＝∑iω
ipci
。
[0058]
本发明技术方案的进一步改进在于：所述步骤s7的具体过程为：
[0059]
步骤s61、将动量信息传回拉格朗日粒子之后，根据动量信息更新拉格朗日粒子的位置，更新位置满足：
[0060][0061]
其中，表示当前时刻拉格朗日粒子p的位置，表示下一时刻拉格朗日粒子p的位置，δt表示一个时间步，表示下一时刻的拉格朗日粒子的速度；
[0062]
步骤s62、更新拉格朗日粒子质量：根据相场参数更新拉格朗日变量的质量，在转移过程不更新质量，一个时间步结束后，更新质量，更新质量满足以下公式：
[0063]
[0064]
其中，表示下一时刻拉格朗日粒子p的质量，c
l
表示l相的质量分数，m
l
表示l相的质量。
[0065]
由于采用了上述技术方案，本发明取得的技术进步是：
[0066]
1、本发明提出了一种对混相、非混相之间的扩散效应，化学反应统一处理的求解器。多相流交互会产生许多丰富的视觉效果，如果将每种效果都分开处理，将会是一个巨大的工作量，在系统的延展性、继承性方面也会有很大的限制。而本发明提出的系统高度统一了多相流交互的特点，该系统适用于几乎所有多相流交互的模拟，可以轻松进行扩展，而不需要做太多修改；
[0067]
2、本发明基于能量的求解器，允许用户自己定义体能函数，来控制整个系统的演化过程。在自然界中，系统朝着能量最小化的方向演化，基于这一原理，本发明提出一种基于能量的求解器，用来模拟部分溶解、萃取等现象，克服了传统多相流求解器的局限性。不同的体能函数驱动得到完全不同的视觉效果，允许用户自己定义体能函数，扩展多相流交互的模拟；
[0068]
3、本发明将mpm与pfm相结合，对整个mpm网格进行建模，可以快速捕捉不同相之间的界面，包括气相与液相的界面。mpm与相场pfm构成的复合模型可以有效结合拉格朗日法与欧拉法的优点。mpm使用欧拉网格和拉格朗日(材料)点来表示材质。在运动过程中，欧拉网格保持不变，拉格朗日点进行移动，可以保留运动的细节。引人pfm，对mpm的背景网格进行建模，不需要对界面进行跟踪，就可以快速捕捉到不同相的界面。包括流体与气体的界面。这种方法不仅在交界面上有物理意义，在体相上也有物理意义。
附图说明
[0069]
图1是本发明的流程框图；
[0070]
图2是本发明哈密顿算子的离散化示意图；
[0071]
图3(a)是本发明的双相流系统中双相不相溶的体能函数图像；
[0072]
图3(b)是本发明的双相流系统中双相部分溶解的体能函数图像；
[0073]
图4(a)是本发明的三相流系统中以吉布斯三角形作为展示的体能函数图像；
[0074]
图4(b)是本发明的三相流系统中以吉布斯三角形能量等高线作为展示的体能函数图像；
[0075]
图4(c)是本发明的吉布斯三角形能量大小与颜色的对比图；
[0076]
图5是本发明气液界面捕捉的效果图；
[0077]
图6是本发明双相流系统中双相不相溶的效果图；
[0078]
图7是本发明双相流系统中双相部分溶解的效果图；
[0079]
图8是本发明三相流系统中萃取的效果图。
具体实施方式
[0080]
下面结合实施例对本发明做进一步详细说明：
[0081]
如图1所示，一种基于mpm和相场模型的多相流交互模拟算法，包括如下步骤：
[0082]
步骤s1、将存储在拉格朗日粒子的质量参数、相场参数、动量参数分别转移到欧拉
网格；
[0083]
步骤s1的具体过程为：
[0084]
步骤s11、将存储在拉格朗日粒子的质量参数转移到欧拉网格：将拉格朗日质量转换为欧拉网格质量，使用apic作为传递方式，该过程满足以下公式：
[0085]
mi＝∑
pmp
ω
ip
，
[0086]
其中，mi表示欧拉节点i的质量，m
p
表示拉格朗日粒子p的质量,ω
ip
是拉格朗日粒子和欧拉网格映射的具有连续性的二次b-样条插值函数中的权重。
[0087]
在混合欧拉-拉格朗日法中，通常使用拉格朗日粒子存放所有的状态信息，欧拉网格用于辅助计算。该过程需要多次进行拉格朗日粒子与欧拉网格的信息传递。因此，选择具有良好数值属性的信息传递方法非常重要。
[0088]
apic是一种拉格朗日粒子与欧拉网格信息的传递方法，该方法可以消除粒子的耗散，防止信息的缺失，确保传输的角动量传递的精确守恒。除了位置和速度，apic方法为每个拉格朗日粒子添加了一个affine(仿射)矩阵，用于记录拉格朗日粒子的散度和角速度。
[0089]
步骤s12、将存储在拉格朗日粒子的动量参数转移到欧拉网格：转移的过程满足动量守恒，使用apic作为传递方式，该过程满足以下公式：
[0090]
(mv)i＝∑
pmp
ω
ip
(v
p
a
p
(x
i-x
p
))，
[0091]
其中，a
p
为拉格朗日粒子p的仿射速度，v
p
表示拉格朗日粒子p的体积，xi为欧拉节点i的位置，x
p
为拉格朗日粒子p的位置。
[0092]
步骤s13、将存储在拉格朗日粒子的相场参数转移到欧拉网格：在求解相场方程中需要多次计算哈密顿算子，为了更加快速的对其进行离散化，相场的更新将在欧拉网格上进行。为了更好的展示相组分的转移，使用质量分数作为相场参数反映多相流体的相容性，对于n相交互系统，只需要求解n-1相场方程，同时转移过程也只需要转移n-1个相场参量，选择mpm传递方式，相场参数转移满足以下公式：
[0093]ci
＝∑
pmpcp
ω
ip
，
[0094]
其中，ci表示欧拉网格的相场参数，c
p
表示拉格朗日粒子的相场参数。
[0095]
物质点法(mpm)是一种伽辽金方法，每一个mpm粒子代表一块连续的材料。使用控制方程的连续描述。物质点法使用欧拉网格进行碰撞检测与动量更新、拉格朗日粒子用于状态跟踪(如：平流和变形)。为了消除拉格朗日粒子的耗散，防止信息的缺失，确保传输的角动量传递的精确守恒，选择使用affine particle-in-cell(apic)进行欧拉网格与拉格朗日粒子之间的信息传递。为防止单元交叉不稳定性，使用具有连续性的二次b-样条插值函数。
[0096]
步骤s2、根据化学势控制整个系统的能量，对欧拉网格相场进行更新；
[0097]
相场法(pfm)是一种处理界面的数学方法。求解相场方程时存在四阶空间导数，由于高阶导数的敏感性，在粒子系统中会引入了时间步限制。在计算速率方面失去了相场可以快速求解的优势。同时，由于粒子系统使用核函数对周围粒子进行插值，因此很难对空气进行建模，对于捕捉气液界面等操作，需要进一步进行单独的处理。因此，选择在欧拉网格上更新相场。
[0098]
步骤s2更新欧拉网格相场具体过程为：
[0099]
步骤s21、求解化学势：化学势表征相组分之间的转移趋势，控制整个系统由高化
学势向低化学势发展，将相场参数插值到欧拉网格之后，采用中心网格、中心差分对其进行离散化，离散化示意图如附图2所示。
[0100]
相场法的精华思想之一就是采用序参量的梯度来追踪界面。多种流体的相态定义了一种称为自由能的势能。相态的演化使得整个系统的能量最小化。
[0101]
相场的控制方程为：
[0102]
f(c)＝∫
ω
[b(c(x,t)) i(c(x,t))]dx
[0103][0104]
其中，ω表示系统的空间区域，k为常数，用于控制扩散速度，i表示界面能，用于减小质量分数场的梯度，使得表面变化更加平滑，b表示体能，控制流体如何进行交互，c(x,t)表示t时刻，x位置的相场取值；
[0105]
化学势的控制方程为：
[0106][0107]
μ表示化学势，b表示体能，k为常数，c表示相场参数，表示拉普拉斯算子；
[0108]
步骤s22、选择体能函数：用体能函数局部最小值控制流体的低能状态，对流体从高能状态演化为低能状态的过程进行模拟，设计不同的体能函数用以实现不用的交互效果，双相交互使用体能函数图像进行展示，对有效区间的极小值进行标记，三相交互采用吉布斯三角形作为能量的展示，设计体能函数遵循：能量的极小值位于有效区内，以确保最终可以达到能量最小化状态；化学势对于某一相的偏导数只应与该相的质量比有关；
[0109]
不同的体能函数可以实现不用的交互效果，附图3(a)为双相流系统中双相不相溶的体能函数图像。该体能函数的极小值的分别在c＝0,c＝1处取得。系统向着能量最小化的方向演化，最终的状态将会达到有效区域的极小值处。也就是流体未发生溶解的纯状态。利用体能函数的区域极小值对该现象进行捕捉，附图6为双相流系统中双相不相溶的效果图。同时向容器中注入的蓝色液体与绿色液体，整个演化过程由体能函数所驱动，由附图3(a)可以看出，有效区域的极小值在两相未发生溶解的纯状态下取得，因此在该体能函数驱动下模拟出双相不相溶的效果。虽然双相不相溶，但是在动力学上依然存在交互，使用相场等效形式下的表面张力对两相之间的界面进行捕捉，从实验效果中可以看出，整个注入过程，蓝色液体与绿色液体呈现较好的分层效果。达到稳定后，可以观察出清晰的界面。由于选用了apic作为信息传递方式，二次b-样条函数作为插值函数，可以观察出，整个交互过程，流体非常“活跃”，不会出现流体越来越少的现象，能量损耗非常小。
[0110]
附图3(b)为双相流系统(l1，l2)部分溶解的体能函数。两个极小值，分别在处取得。假设流体初始为纯相。则对于纯l1相，到达极小值对于纯l2相，到达极小值自由能函数驱动系统达到能量最小化，利用局部极小值点捕捉到部分溶解现象。附图7为双相流系统中双相部分溶解的效果图。
[0111]
相同体积的红色液体与绿色液体，实现部分溶解。以红色液体作为纯l1相，绿色液体作为纯l2相。交互过程由附图3(b)给出的体能函数所驱动。对于红色液体而言，区域极小
值在处取得，因此达到稳态之后，呈现红色占比60％，绿色占比40％的偏红色的饱和物。对于绿色液体而言，区域极小值在处取得，因此达到稳态之后，呈现红色占比40％，绿色占比60％的偏绿色的饱和物，最终产生两种完全饱和的混合物。使用相场等效形式施加表面张力，可以观察出，两相之间的清晰界面。附图7与附图6实现了双相交互的两种不同状态，两者交互状态的不同完全由体能函数所驱动，因此，本发明提出的求解器允许用户设计不同的体能函数来对多相流丰富的交互效果进行模拟。对于混相、非混相不需要分开处理，只需要设计合适的体能函数即可，因此，本发明提出的求解器是一种混相、非混相统一处理的求解器。
[0112]
附图4(a)为三相流系统(a,b,c)萃取模拟中以吉布斯三角形作为展示的体能函数。选择使用吉布斯三角形进行能量展示。三角形顶点表示流体纯状态(100％质量的a,b,或c流体)平行线的质量分数相同。附图4(a)中“·”处表示ca＝0.4,cb＝0.5,cc＝0.1。附图4(b)三相流系统(a,b,c)萃取模拟中以吉布斯三角形能量等高线作为展示的体能函数。体能的大小由彩色绘制，白色虚线为等高线。每种颜色对应的能量大小由附图4(c)所示。由能量图可以观察出：相b溶解在相a中比其溶解在相c中具有更低的最小能量，因此，将相a添加到相b、相c的混合物中后，原本溶解在相c中的相b会重新溶解在相a中，实现化学萃取效果。附图8为三相流系统中萃取的效果图。初始状态下层液体为红色液体与绿色液体的混合物，上层蓝色液体为萃取剂。充分摇晃分液漏斗。整个交互过程由附图4(b)给出的体能函数所驱动。充分摇晃分液漏斗，使蓝色液体与混合物充分混合。由于红色液体在蓝色液体中的溶解度更大，蓝色液体将红色液体从混合物中提取出来。静止之后，上层液体是初始混合物红色液体被提取的剩余液体，下层液体为红色液体重新溶解于蓝色液体之后的混合物。实现化学萃取效果。
[0113]
步骤s23、更新相场：针对不同的多相流交互设计对应的体能函数之后，使用化学势对相场进行更新，将简并迁移率设为常数，采用中心网格、中心差分对其进行离散化，相场更新满足以下公式：
[0114][0115][0116]
其中，表示物质导数，m表示简并迁移率，μ表示化学势。
[0117]
步骤s3：对于mpm动量守恒的离散化基于材料的势能，在演化过程中考虑表面张力，用等效相场形式取代基于水平集的表面张力公式。假设材料的内能包括与表面张力，弹性力，重力相关的势能。
[0118]
根据重力和弹性力对欧拉网格动量进行更新遵循如下公式：
[0119][0120]
其中，fe表示弹性力，δx表示欧拉网格的网格间距，为拉格朗日粒子的初始体积，p表示皮奥拉基尔霍夫应力，表示t 1时刻拉格朗日粒子p的形变梯度，表示拉格朗日粒子p的形变梯度矩阵的转置，t表示转置，fg表示重力，g表示重力加速度，fi表示欧拉网格i受力的总和，t表示当前时刻，t 1表示下一时刻，δt表示一个时间步，表示欧拉网格i下一时刻的速度，表示欧拉网格i当前时刻的速度。
[0121]
等效相场形式的表面张力是指液体表面由于分子引力不均衡而产生的沿表面作用于任一界线上的张力。在相场方法中，两种不混相流体之间的经典尖锐界面被过渡区域取代。在过渡区域中流体的性质不断变化，采用等效相场替代表面张力。附图5是使用等效相场替代表面张力对气液界面进行捕捉的效果图。传统的基于拉格朗日粒子的多相流求解器，很难对气液界面进行捕捉。本发明提出的求解器通过对整个mpm网格进行建模，将气体视为双相流中的一种相，对气体进行处理。使用等效相场替代表面张力对气液界面进行捕捉。由效果图可以看出表面层粒子的移动使得表面积逐渐减小，符合表面张力的施加效果。根据重力、弹性力和更新后的等效相场形式的表面张力对欧拉网格动量进行更新。
[0122]
使用更新的相场，对不同组分之间的界面进行捕捉，施加表面张力，进一步对动量进行更新，表面张力控制方程满足：
[0123][0124][0125][0126]
其中，f
σ
是使用相场形式等效的表面张力，表示l1、l2相之间的表面张力系数，表示l1相的相场参数，表示l2相的相场参数，ε与扩散界面相关的系数。
[0127]
步骤s4：更新完网格动量之后，在欧拉网格上对边界采用滑移边界条件进行处理；
[0128]
由于基于粒子的碰撞会使形变梯度产生“不一致性”，导致模拟中出现“伪影”现象。为了避免“伪影”的现象，更新完网格动量之后，在欧拉网格上对边界采用滑移边界条件进行处理，采用滑移边界条件：
[0129][0130]
其中，表示碰撞尚未结束的网格速度，在更新拉格朗日粒子位置之前，多次对拉格朗日粒子对应用碰撞处理，n表示曲面法线。
[0131]
步骤s5、更新完相场与动量之后，将动量参数和相场参数从欧拉网格变量转换到拉格朗日变量，基本流程与拉格朗日粒子到欧拉网格的信息传递保持一致，并不对质量m进行转移，因为在一个时间步内，系统默认拉格朗日粒子质量保持不变。
[0132]
步骤s5的具体过程为：
[0133]
步骤s51、在欧拉网格上更新完动量参数之后，将动量参数从欧拉网格传回拉格朗日粒子：为了与拉格朗日粒子到欧拉网格的信息传递保持高度一致，采用apic作为传递方式，动量参数转移满足：
[0134]vp
＝∑iω
ip
vi[0135]
其中，v
p
表示拉格朗日粒子p的速度，vi表示欧拉网格i的速度；
[0136]
步骤s52、在欧拉网格上更新完相场参数之后，将相场参数从欧拉网格传回拉格朗日粒子，转移回拉格朗日粒子的相场信息为每种相组分的质量分数，该过程遵循以下公式：
[0137]cp
＝∑iω
ipci
。
[0138]
步骤s6、更新拉格朗日粒子状态：根据传回拉格朗日粒子的信息，对拉格朗日粒子的位置、质量状态进行更新。
[0139]
步骤s6的具体过程为：
[0140]
步骤s61、将动量信息传回拉格朗日粒子之后，根据动量信息更新拉格朗日粒子的位置，由于在欧拉网格上已经对边界进行了限制，因此在此步骤不需要对拉格朗日粒子的位置进行限制。更新位置满足：
[0141][0142]
其中，表示当前时刻拉格朗日粒子p的位置，表示下一时刻拉格朗日粒子p的位置，δt表示一个时间步，表示下一时刻的拉格朗日粒子的速度；
[0143]
步骤s62、更新拉格朗日粒子质量：根据相场参数更新拉格朗日变量的质量，由于每种流体物质的特性(如密度)不同，允许对不同相设置不同的参数。在本发明中，允许每种相具有不同的质量，根据相场参数更新拉格朗日变量的质量。为了与mpm解算器高度契合，在转移过程不更新质量，一个时间步结束后，更新质量，更新质量满足以下公式：
[0144][0145]
其中，表示下一时刻拉格朗日粒子p的质量，c
l
表示l相的质量分数，m
l
表示l相的质量。