CN113486566A - 一种使用虚拟粒子相进行多孔介质和多相流体混合仿真的方法 - Google Patents

Abstract

一种使用虚拟粒子相进行多孔介质和多相流体混合仿真的方法，使用虚拟粒子相，计算固体外部流体和内部流体的吸收相，进而计算出不同相之间的质量和动量传输，并基于SPH方法得到离散化实现。本方法中，内部和外部流体使用体积分数和混合模型统一建模，避免了恒定的粒子删除和插入，实现了真实的固体流体混合仿真效果，并提供了艺术化控制的方案，对多孔介质和流体混合仿真效率和效果与当前业界先进方法的效果可比。

Description

一种使用虚拟粒子相进行多孔介质和多相流体混合仿真的 方法

技术领域

本发明涉及计算机图形学流体动力学模拟领域，尤其涉及多孔固体介质与多组分流体的动力学仿真方法与系统。

背景技术

对于多孔介质的研究工作在计算机图形学的很多领域都在进行着，例如涉及布料、沙子或泥土和体积固体的计算机图形学工作。多孔介质不同于只在表面与流体发生作用的不可渗透固体，在于流体发生作用时，多孔介质表现出更强的耦合性，并且和流体的相互作用同时发生在介质的表面和内部。这些流动现象远比开放空间中的单相流体流动复杂得多，它们可能涉及许多同时发生的物理过程，如多孔固体中的流体传质、固体表面的吸收和发射、疏水和亲水行为、毛细管效应和多相流。

Lenaerts在“Porous flowin particle-based fluid simulations”中通过删除和插入流体粒子并跟踪固体粒子内部的流体流动来模拟流体的吸收和逸出。因此，模拟器必须不断删除和插入不均匀的SPH粒子。这种方法对多流体模拟的简单扩展将拆分流体粒子，因为不同的流体相可能会以不同的方式与固体相互作用。这将是一个消耗大量计算资源过程，并且会产生许多破碎的粒子。

在不分裂流体粒子的情况下，需要处理不同相之间的质量和动量传递，并保持流体粒子穿过固体边界的一致性。Tampubolon在“MultiSpecies Simulation of PorousSand and Water Mixtures”使用两层网格分别模拟沙子和流体运动。对于砂状多孔材料来说，可以使用一组类似的纳维-斯托克斯方程来描述纯流体和固液混合物，因此重点放在这两个网格层之间的相互作用上。但是对于海绵等其他多孔材料，多孔固体内外的流体运动可能会有很大差异(多孔介质中的流动遵循达西定律)，还可能涉及其他物理过程，如毛细效应。因此，这种方法难以处理多种流体和多孔介质发生相互作用的情况。

发明内容

本发明的目的是解决当前粒子方法仿真技术均无法一般性地处理多相流体与多孔材质共同仿真的问题，提出了一种新颖的基于虚拟粒子相的流体模拟算法，用于多孔介质材料和多相流体进行联合动力学仿真。我们的方法可以使用体积分数和混合模型实现固体内部和外部流体运动模型的统一建模，并且可以避免在多孔介质材料和多相流体联合仿真过程中的粒子删除和插入。新方法能够描述流体在多孔介质中运动时产生的多种现象，包括多孔弹性现象、毛细管效应和多流体环境中的可变吸收。该算法框架在实现了高真实感的流体模拟的同时，保证了流体运动的高效求解。

本发明提供的使用虚拟粒子相进行多孔介质和多相流体混合仿真的方法是，使用SPH方法对单相或多相流体与多孔介质的联合仿真进行处理；所述SPH方法使用固体粒子集合表示多孔介质固体，使用液体粒子集合表示多相流体液体；

在计算液体粒子的运动时对多相流体的每一真实流体相，视其位于固体内部和固体外部的部分为两种“虚拟粒子相”，从而对所有液体粒子采用统一的混合流模型物理方程计算其运动，并且不在流体进出固体时进行粒子删除、增加操作。

所述的使用SPH方法对单相或多相流体与多孔介质的联合仿真进行处理的仿真步骤包括：

第1：初始化流体场景的粒子属性，包括位置，速度，粒子类别为固体粒子还是流体粒子，若为流体粒子，其所包含的各相的体积分数；

第2：计算所有固体粒子的密度ρ_s，i和流体粒子的密度ρ_f，i；

第3：计算流体粒子的N_f，即粒子邻域内通用SPH核函数的和；

第4：计算外部流体的流体粒子的粘度造成的加速度a^visc和由压力造成的加速度a^press；

该步中在分别计算由粘度和压力造成的加速度时，使用类似“使用压强边界的IISPH方法”来计算固体外部区域流体粒子的流体动压强和a^press，提供了更好的外部多相流体不可压缩性。计算固体外部流体动压强和a^press的方法是：

(1)：在算法中统一使用粒子中的流体所占体积替代原方法中的粒子整体体积进行计算；

(2)：在算法中统一使用调和平均方法替代原方法中的直接平均方法进行计算；

第5：计算固体粒子的孔隙压力、达西通量、固体内部流体速度、固体内应力，进而计算液体粒子与固体相互作用产生的加速度a^pore；

第6：计算混合流模型物理方程中其他控制流体粒子运动的项，即外部流体漂移速度u_mk，α，内部流体漂移速度um_k，β，吸收通量D_k，因毛细作用力引起的加速度a^cap；

第7：对固体与流体粒子分别根据弹性固体模型和混合流模型物理方程进行时间积分，推进一个时间步；

第8：重复第2到第7步，进行下一个模拟时间步计算。

本发明的技术特色：

1)基于SPH算法的多孔介质中的流体动力学仿真方法

该方法提出了一种基于SPH算法的新颖的多孔介质与流体联合仿真技术。在本技术方案所提出的模型中，考虑存在三种耦合物理机制的一般情况，即多孔固体运动、内部流体传输和外部流体运动。并使用三组控制方程来描述多流体流动和多孔介质材料的耦合物理。首先，对于多孔介质材料内部的流体流动，使用达西定律来模拟流动运动。接下来，对于多孔介质材料外的流体流动，采用多流体混合模型进行描述。最后，对于多孔介质，采用弹性固体模型进行描述。

2)基于虚拟粒子相概念的灵活本构建模策略，支持多流体流动、多孔介质材料变形以及流体和多孔介质固体之间的耦合物理。

本技术方案提出了虚拟粒子相的概念，以统一的方式处理穿过界面的流体粒子的状态变化以及相关的质量和动量传输，从而避免显式状态跟踪和运行时粒子删除/插入。借助虚拟分数概念，本方案可以以通用方式解决固体内部和外部流体粒子动力学问题。然后通过求解多流体模型中的连续性方程来处理与多孔介质固体有关的不同相行为下复杂的内外状态转变。

本发明的优点和有益效果：

本发明从理论上，提出了一个统一的粒子框架来模拟多孔介质中的流体流动现象。引入了一种新颖的虚拟粒子相概念，以避免显式粒子状态跟踪和运行时粒子删除/插入。该技术方案所提出的统一模型灵活且稳定，可以确保整个模拟空间的质量和动量传输一致，并且可以处理多种流体和多孔介质的耦合物理现象求解。所提出的技术方案可以很好的描述多孔介质材料和流体相互作用时的物理现象，例如多孔弹性现象，毛细现象以及可变吸收现象等。

附图说明

图1展示了仿真过程中与多孔行为相关的参数的影响。

图2是多孔介质和单相流体的相互作用。

图3是两种吸收率不同的流体同时通过多孔介质的仿真结果。

图4是多种性质的多孔介质和不同吸收率的混合流体之间的相互作用。

图5是沙状多孔介质材料和流体之间的相互作用。

图6是多相流体混合物的选择性过滤。

具体实施方式

由于本发明提出了一种新颖理论基础上的仿真方法，在本部分中将首先介绍本发明所提出的理论，而后说明如何使用各理论公式实施仿真计算。最后将提供实施效果附图。

具体实施中所需要的理论公式和物理模型方程

1.多孔介质与流体联合仿真的物理模型

首先，我们使用达西定律来模拟多孔介质内部的流体运动。具体而言，均质可渗透材料中的瞬时流量，即达西通量由下式确定：

其中下标k表示第k个流体的相，q_k表示第k相的达西通量，k_wk表示第k相的渗透率张量，μ_k表示第k相的流体粘度，

表示空隙空间中流体的压力。达西定律中的孔隙压力由固体和局部吸收流体的状态决定，并用公式1表示，它将吸收流体从有孔隙压力的区域驱动到低孔隙压力区域。

接下来，对于多孔介质材料外的流体流动，采用多流体混合模型进行描述。

其中

是第k相的体积分数，u_m和

是混合速度和第k相的漂移速度，而g、a^press、a^other分别表示重力，压力引起的加速度和其他影响因素(粘度)引起的加速度。上面的等式中使用波浪符号来将它们与我们将要介绍的“虚拟粒子相”量区分开来。混合模型通过局部体积上的平均值来表达连续性和动量方程，对于SPH来说，局部体积是粒子。它可以捕获相位速度差，但只需要通过解析计算漂移速度来求解粒子速度。在混合模型求解器中，可以使用标准粒子推进方案，粒子内的相体积分数变化由连续性方程自动处理。

最后，对于多孔介质，采用弹性固体模型进行描述。应变张量∈表示为：

其中，F是变形梯度张量，I是身份矩阵。为了考虑孔隙压力对固体变形的影响，我们在他们的压力张量中添加了一个

并表达了压力张量P如下：

其中ν、λ是Lam′e参数，而η是实常数。最后一项代表孔隙压力的影响，即吸收流体将使固体膨胀。

上述控制方程可以用现有方法单独求解，但它们在多孔介质材料表面附近变得强耦合，大大增加了求解的复杂性。本方案提出了虚拟粒子相的概念来以统一的方式处理穿过界面的流体粒子的状态变化以及相关的质量和动量传输，从而避免显式状态跟踪和运行时粒子删除和插入。

2.引入虚拟粒子相后的本方法所采用的统一流体物理方程

我们为每个流体粒子分配两个虚拟粒子相α_k和β_k，分别代表第k相在外部和内部区域的体积分数，并且它们满足关系∑_k(α_k+β_k)＝1。因此,∑_kα_k＝1表示流体粒子位于多孔介质材料外部，∑_kβ_k＝1表示粒子位于多孔介质材料内部，非零的α_k和β_k值表示流体粒子靠近多孔介质材料的表面。

考虑一个示例，其中处于两相状态的流体粒子从固体外部区域移动到内部区域，但只有一种流体相可以被吸收到多孔介质材料中，而另一种则不能。该流体粒子的质量传递、动量和位置变化受所有物质相的影响，每个相依次可能具有独立的内部和外部部分，遵循内部和外部区域的不同规律。虚拟粒子相概念区分了这些分数，并且虚拟粒子相可以遵循单独的控制方程。虚拟粒子相具有α_k和β_k可以通过吸收或其他机制在多孔介质材料表面附近转移而不违反质量守恒定律的性质，但实相之间的体积分数不能简单转移。

借助虚拟粒子相概念，可以使用通用方式解决固体内部和外部流体粒子动力学问题。然后通过求解多流体模型中的连续性方程来处理与多孔介质固体有关的不同相行为下复杂的内外状态转变。所有流体粒子的通用方程为：

其中α_k和β_k分别为第k相的外部和内部体积分数，u_m，umk，α，u_mk，β分别为流体粒子速度以及外部流体和内部流体的漂移速度。D_k和

分别是吸收通量和吸收源。g、a^press、a^pore、a^cap、a^other分别表示重力、压力引起的加速度、孔隙压力引起的加速度、毛细管力引起的加速度和其他因素(例如粘度)引起的加速度。

这种扩展的混合模型可以通用方式模拟流体粒子，并处理控制方程之间的耦合。公式(6)(7)是连续性方程，其中外部流体由公式(6)计算，内部流体由公式(7)计算，因为它们遵循的不同质量传输定律导致u_mk，α，um_k，β不同。纯内部流体粒子只会受a^pore影响，它将达西定律与弹性固体模型联系起来。纯外部流体粒子的α^pore为零，但固体粒子可以提供α^press中的边界压力，用于求解外区流体运动。对于固体表面附近的混合粒子，这些项都可以是非零的。

3.外部流体处理

方案使用类似“压强边界的IISPH方法”(Pressure Boundaries for ImplicitIncompressible SPH,Band et Al.2018)来计算外部区域流体粒子的流体动压，这提供了更好的不可压缩性。随后计算压力并实现双向耦合，并进行了一系列修改以支持多流体的多孔介质与流体联合模拟。

首先通过下式计算流体粘度：

弹性固体粒子的粘度可表示为：

我们使用i、j来表明计算适用于所有相邻的流体和固体粒子，角标f、s表明该量为液体粒子、固体粒子所携带存储。其中

μ_mf、μ_ms是粒子的聚集粘度系数，W是SPH算法的核，ρ_i，ρ_s是粒子密度，

α_fk为虚拟粒子相体积分数，m_s为固体质量。我们设置流体和固体粒子的中间加速度为

对于外部区域计算，固体粒子充当流体压力计算的场景边界。

我们将原方法中的流体粒子的体积更新为外部流体的有效体积，以在计算中排除吸收的流体。在每个时间步的开始，我们计算每个流体粒子的外部流体体积分数之和α_f＝∑_kα_k，之后它的剩余体积近似为V_f0＝α_fh³，其中h是光滑核半径。由于这个固体附近流体粒子的体积变化很大，为了更好的稳定性，我们使用调和平均法，这样速度的发散度就变成了

其中u^*是中间速度。使用中间速度，可以推导出一组线性方程(Ap)_f＝S_f用于粒子并求解粒子压力p，其中A是矩阵，S_f是导出项。例如对于粒子i：

我们进一步修改S_f以反映剩余体积的修正：

使用上述修改，我们可以求解每个流体粒子和固体粒子(被视为外部流体的边界粒子)的压力值。然后我们可以计算公式8中的a^press项。并计算固体粒子从外部流体接收的边界力以进行固体运动积分。由于V_f0＝α_fh³，在原方法与本方法中V_i，V_j是它的插值值，公式13确保a^press即使没有明确的区域标签，内部流体粒子上的该项计算结果也会自动为0。

该方案比原来的多流体压力计算可以获得更好的视觉不可压缩性，并在流体撞击固体时提供实用的固液机械运动解决方案。

4.计算孔隙压力和达西通量

在固体粒子上，我们计算与达西定律相关的物理项，以便以后计算多流体在多孔介质中流动的物理模型，并计算固体运动积分所需的项。

首先通过公式4计算固体粒子上的固体应变，为计算固体粒子的空隙压力做准备。内部区域固流体动量平衡的关键因素是孔隙压力和达西通量。当吸收更多的流体或固体受到挤压时，孔隙压力升高，而达西通量与其梯度成正比。基于多孔弹性理论，固体中的孔隙压力由下式获得：

其中p₀是静止孔隙压力，tr(∈)是固体的体积应变，B是Biot比率常数，φ_sk表示第k相的吸收体积率。上述孔隙压力是在固体粒子上计算的，可直接用于固体运动计算。然后我们通过收集附近的固体粒子信息来计算公式1中达西通量的流体粒子梯度。

吸收体积率φ_sk表示在固体粒子位置上有多少流体相k被局部吸收，它的计算方式为

其中V_s0表示固体粒子的静止体积，ρ_k表示第k相的静止密度，N_f＝∑_sW_fs表示每个流体粒子与其邻域中固体粒子之间的SPH核的总和。对于任何具有非零β_fk的流体粒子，我们使用SPH核权重有效地将其吸收质量分配给附近的固体粒子，公式17中使用的吸收质量是固体粒子从其周围流体粒子中收集到的吸收质量。

随后通过公式16和下面的公式进行达西通量的计算。其中N_s由s邻域中的固体粒子组成。

孔隙压力影响固体的行为，如公式5所述，它也是固体内部流体流动的驱动力。孔隙压力被施加到内部流体，如果没有其他力存在，它的贡献在下一个模拟步骤结束，使得流体和固体的相对速度u_rk与公式1中的达西通量一致，可以通过u_rk＝q_k/e得到这个结论。其中e称为孔隙率，描述空隙体积占总体积的比例。因此，我们在公式8中添加a^pore来模拟内部和外部流体模拟中的这个力源。具体来说，我们使用

其中

是内部静止密度，u_βk＝u_rk+u_s是第k相的世界坐标相速度，u_s是固体速度，ρ_fm＝∑_k(α_k+β_k)ρ_k是流体粒子的聚集密度。

最后，通过公式5获得固体应力，并通过外部流体的物理模型得到固体粒子的边界力。

5.流体物理方程中其他流体相关项的计算

使用固体粒子上计算的物理量，我们可以进一步的计算出流体控制方程(公式6，7，8)所需的剩余项。首先使用上节方法在固体粒子上计算出相关物理量值。然后，使用下列公式使用SPH公式插值到流体粒子上

uβ_k，f＝∑_su_βk，sW_fs/N_f (公式20)

q_k，f＝∑_sq_k，sW_fs/N_f (公式21)

当N_f为0时，将方程的左侧也置为0。随后，计算流体粒子的漂移速度。对于外部流体虚拟粒子相α_k，其漂移速度可以通过下式来进行计算:

其中

表示外部质量分数，

表示外部混合物密度，

τ，σ，

是常数权重因子，

是流体s孔隙压力的梯度，它可以是特定相的，因为不同的相通常吸收不同。公式22与原始漂移速度的区别仅在于最后一项，它模拟了由于多孔固体表面附近的吸收或发射引起的吸入或推进效应。

内部流体虚拟粒子相β_k在多孔固体内遵循达西定律，并且它没有与外部区域相同的流体动压。由于内部区域的孔隙压力与外部区域流体的流体动压力起着相同的作用，我们替换原始漂移速度方程的第二项，得到我们的内部流体漂移速度模型：

其中c_βk表示内部质量分数，τ′，σ′是常数权重因子，撇号表示它们可能因为内部和外部区域环境设置不同而与公式22中的不同。在公式23中，q_k是在局部立体坐标中还是在世界坐标中计算并不重要。只要固体孔隙率是均质的，由于∑_kc_βk＝1，可以使用任一坐标下计算的q_k得到相同的结果。

下面进行吸收源的计算。类似于公式1中的达西通量计算，通过计算孔隙压力梯度以获得吸收通量：

其中K是与固体相关的恒定吸收比系数。无量纲项(α_k+φ_sk)相乘以避免产生假通量。然后，通过公式7和公式6确定内部和外部流体的交换。

通过上面的式子。我们可以得到

然后通过下式计算孔隙压力梯度以得到吸收通量。

当一部分流体质量被固体吸收时，流体粒子将在固体表面附近受到毛细作用力，阻止流体从多孔固体中逸出。我们还需要计算因毛细作用力引起的加速度a^cap，计算公式如下：

a^cap＝ψ∑_kβ_kn (公式27)

其中ψ是恒定强度因子，n是从流体粒子到多孔固体内部方向的法向量。在多孔介质表面附近，毛细管力产生向内拖拽效应。粒子流入固体后，会逐渐减小到零，从而使流体粒子的运动完全受达西定律的支配。

在流体粒子上使用SPH公式，可以将公式27离散化为：

6.对固体粒子和流体粒子进行时间积分

在完成上述计算后，我们已经获得了固体、流体物理方程(公式4-公式8)中所有所需的项的结果，我们对SPH固体粒子根据公式4，5进行时间积分，对SPH流体粒子根据公式6，7，8进行时间积分，将时间步向前推进一步。

实施例1

本发明提供了一种使用虚拟粒子相进行多孔介质和多相流体混合仿真的方法。本发明的方法根据上面介绍的理论，可按如下步骤组织和编写算法程序，进而运行该算法程序实施实际仿真。程序算法执行步骤如下：

第1：初始化流体场景的粒子属性，包括位置，速度，粒子类别为固体粒子还是流体粒子，若为流体粒子，其所包含的各相的体积分数；

在初始化后，在每一个时间步的计算中，

第2：使用标准SPH插值方法计算所有固体粒子的密度ρ_s，i和流体粒子的密度ρ_f，i；

第3：计算流体粒子的N_f，即粒子邻域内通用SPH核函数的和；

第4：计算外部流体的流体粒子的粘度造成的加速度

(公式9-10)和由压力造成的加速度a^press；

该步中在分别计算由粘度和压力造成的加速度时，使用类似“使用压强边界的IISPH方法”来计算固体外部区域流体粒子的流体动压强和a^press，提供了更好的外部多相流体不可压缩性。计算固体外部流体动压强和a^press的方法是(公式11-15)：

(1)：在算法中统一使用粒子中的流体所占体积替代原方法中的粒子整体体积进行计算；

(2)：在算法中统一使用调和平均方法替代原方法中的直接平均方法进行计算；

第5：计算固体粒子的孔隙压力、达西通量、固体内部流体速度、固体内应力，进而计算液体粒子与固体相互作用产生的加速度a^pore(公式16-19)；

第6：计算混合流模型物理方程中其他控制流体粒子运动的项u_mk，α，u_mk，β，D_k，a^cap(公式20-28)；

第7：对固体与流体粒子分别根据弹性固体模型(公式4-5)和混合流模型物理方程(公式6-8)进行时间积分，推进一个时间步；

第8：重复第2到第7步，进行下一个模拟时间步计算。

实施效果

本部分展示的内容为使用上述实施方案进行多相流体与多孔固体共同仿真的输出经过渲染后的结果。仿真主体置于由不同纹理(白色、棋盘格、风景等)贴图生成的墙面与地面之间。

图1，展示了多孔介质相关的参数对仿真效果的影响。从左到右：标准，更高的吸收率，更高的渗透率，更高的p_0k。更高的吸收率使流体更快地进入固体。更高的渗透性导致固体内更高的流体传输速度，并在开始时导致更浅的颜色。更高的p_0k增加了流体吸收和运输速度，最终导致更深的颜色。

图2，展示了吸收单一流体和单一多孔介质的仿真结果。当流体到达固体表面时，它会引起固体变形，同时一部分流体被反弹回来。与此同时，一些流体被多孔固体材料吸收，当大量流体回落时，会留下可视化的痕迹。在左下角，我们可以观察到固体下端的橙色部分发生体积膨胀，显示出多孔弹性效应。由于不需要删除或插入粒子，我们的方法保证了整个模拟的质量守恒。

图3，展示了两种流体通过多孔介质的仿真结果。被吸收较少的流体相可以更快地通过固体区域。两个流体相同时通过穿孔材料，由于两种流体相被固体吸收的方式不同，因此其中一相大部分被固体吸收，而吸收较少的相几乎可以自由地穿过多孔介质。

图4，展示了有多种多孔介质时的仿真结果。固体由三部分组成，对流体的三种相有着不同的吸收特性。仿真结束后不同区域的固体会呈现不同的颜色，在这个例子中，我们只调整每个固体区域内的物理多孔参数，而没有在模拟中明确标记区域。当不同的固体区域接触时，模拟结果显示了我们方法的正确性，并且还实现了流体在固体中的分离行为。

图5，展示了我们的方法也可以实现沙状多孔材料与流体的联合仿真。结果表明，我们的方法也能够处理沙状多孔材料。

图6，展示了三相流体混合物的选择性过滤的仿真结果。黑色流体混合物流过由三层泡沫组成的垂直过滤器，每一层都吸收一种流体相。当混合物被选择性过滤时，不同的多孔介质会变成不同的颜色。

Claims (5)

1.一种使用虚拟粒子相进行多孔介质和多相流体混合仿真的方法，其特征在于，包括：

使用SPH方法对单相或多相流体与多孔介质的联合仿真进行处理；

所述SPH方法使用固体粒子集合表示多孔介质固体，使用液体粒子集合表示多相流体液体；

还包括：

计算液体粒子的运动时对多相流体的每一真实流体相，视其位于固体内部和固体外部的部分为两种“虚拟粒子相”，从而对所有液体粒子采用统一的混合流模型物理方程计算其运动，并且不在流体进出固体时进行粒子删除、增加操作。

2.如权利要求1所述的使用虚拟粒子相进行多孔介质和多相流体混合仿真的方法，其特征在于，所述的使用SPH方法对单相或多相流体与多孔介质的联合仿真进行处理的仿真步骤包括：

第1：初始化流体场景的粒子属性，包括位置，速度，粒子类别为固体粒子还是流体粒子，若为流体粒子，其所包含的各相的体积分数；

第2：计算所有固体粒子的密度ρ_s，i和流体粒子的密度ρ_f，i；

第3：计算流体粒子的N_f，即粒子邻域内通用SPH核函数的和；

第4：计算外部流体的流体粒子的粘度造成的加速度a^visc和由压力造成的加速度a^press；

第5：计算固体粒子的孔隙压力、达西通量、固体内部流体速度、固体内应力，进而计算液体粒子与固体相互作用产生的加速度a^pore；

第6：计算混合流模型物理方程中其他控制流体粒子运动的项，即外部流体漂移速度u_mk，α，内部流体漂移速度u_mk，β，吸收通量D_k，因毛细作用力引起的加速度a^cap；

第7：对固体与流体粒子分别根据弹性固体模型和混合流模型物理方程进行时间积分，推进一个时间步；

第8：重复第2到第7步，进行下一个模拟时间步计算。

3.如权利要求2所述的使用虚拟粒子相进行多孔介质和多相流体混合仿真的方法，其特征在于，第4步中在分别计算由粘度和压力造成的加速度时，使用类似“使用压强边界的IISPH方法”来计算固体外部区域流体粒子的流体动压强和a^press，提供了更好的外部多相流体不可压缩性。

4.如权利要求3所述的使用虚拟粒子相进行多孔介质和多相流体混合仿真的方法，其特征在于，所述的类似“使用压强边界的IISPH方法”计算固体外部流体动压强和a^press的方法是：

(1)：在算法中统一使用粒子中的流体所占体积替代原方法中的粒子整体体积进行计算；

(2)：在算法中统一使用调和平均方法替代原方法中的直接平均方法进行计算。

5.如权利要求2所述的使用虚拟粒子相进行多孔介质和多相流体混合仿真的方法，其特征在于，第5步中计算的孔隙压力由固体应变与固体局部吸收液体量共同决定。

Images