CN104318598B - 一种三维流固单向耦合的实现方法及系统 - Google Patents

Abstract

本发明公开了一种三维流固单向耦合的实现方法及系统，其中，该方法包括：将固体边界的三角形网格离散成粒子，并进行边界粒子均匀化采样；计算该固体中流体粒子的状态，包括：流体与固体边界处流体粒子的密度、流体粒子的压力及粘力，以及流体粒子表面张力；并利用一数值积分模式对流体粒子的速度和位置信息进行更新；当所述流体粒子穿透所述固体边界时，对所述流体粒子的速度和位置进行校正，从而实现三维流固单向耦合。本发明公开的方法及系统增强了边界处理时的稳定性，成功模拟出压力、粘力和阻止穿透效果的同时大大减小了计算的复杂度，并且有效模拟出表面张力的宏观效果，而且并不会出现直接模拟聚合力情况下的网状失真现象。

Description

一种三维流固单向耦合的实现方法及系统

技术领域

本发明涉及基于物理的计算机实时动画技术领域，尤其涉及一种三维流固单向耦合的实现方法及系统。

背景技术

拉格朗日流体模拟在计算机动画领域一直是个热门的研究议题，其基于粒子的流体表示天然地支持模拟一些微小尺度的流体现象，同时无条件地保证质量守恒。光滑粒子动力学SPH(Smoothed Particle Hydrodynamics)方法是一种典型的拉格朗日方法，由Lucy、Monaghan&Gingold等人于1977年提出用来解决天体物理学中的行星运动问题，后来被引入到计算流体力学领域。在模拟流体时它是将连续的流体用相互作用的质点组描述，各个质点上承载各种物理量，包括质量、速度、密度等，通过求解质点组的动力学偏微分方程组和跟踪每个质点的运动，求得整个系统的力学行为。SPH方法使用核函数对密度、压力、粘力项进行离散化，得到Navier-Stokes方程组的离散化计算形式，从而在每步迭代过程近似解出各个物理量，模拟流体运动。

关于SPH方法模拟流固耦合现象属于流固两相交界处边界处理的范畴，由于SPH方法的粒子属性导致其在处理多相流问题时需要花费的计算代价很大，目前不存在一个统一的模型和方法能完美的解决固液耦合问题，几种具有代表性的方法如下：

1)Muller等人于2004年提出的基于Lennard-Jones势能的固液交互模型，其为流体粒子和固体粒子间施加一个距离相关的力，该方法在宏观上可以近似模拟，但在微观尺度上不满足流体运动的物理规律，而且计算所需的时间步长较大。

2)Seungtaik Oh等人于2009年提出的基于冲量边界力(IBF)方法，其为固液交互所施加的边界力的计算基于粒子间碰撞过程中产生的冲量大小和方向，该方法需要计算虚拟的冲量边界力。

3)Becker等人于2009年引入direct force来为边界力建模，该方法采用预测校正模式对与边界碰撞后粒子的位置和速度进行校正，达到没有穿透现象的发生并模拟了多种边界条件，该方法对于流固双向耦合计算的复杂高，每个迭代中需要两次周围邻居粒子查询。

发明内容

本发明的目的是提供一种三维流固单向耦合的实现方法及系统，增强了边界处理时的稳定性，成功模拟出压力、粘力和阻止穿透效果的同时大大减小了计算的复杂度，并且有效模拟出表面张力的宏观效果，而且并不会出现直接模拟聚合力情况下的网状失真现象。

本发明的目的是通过以下技术方案实现的：

一种三维流固单向耦合的实现方法，该方法包括：

将固体边界的三角形网格离散成粒子，并进行边界粒子均匀化采样；

计算该固体中流体粒子的状态，包括：流体与固体边界处流体粒子的密度、流体粒子的压力及粘力，以及流体粒子表面张力；

并利用一数值积分模式对流体粒子的速度和位置信息进行更新；

当所述流体粒子穿透所述固体边界时，对所述流体粒子的速度和位置进行校正，从而实现三维流固单向耦合。

进一步的，所述将固体边界的三角形网格离散成粒子，并进行边界粒子均匀化采样包括：

读取固体模型的三角形网格，并建立固体模型的周围符号距离场，具体的：将固体模型中每个三角形网格建立一个包围盒，并为每个包围盒以统一的方式划分网格；计算每个网格点到固体模型的符号距离，符号距离场以一个哈希表来存储，形成一个存在于固体模型周围的窄带宽的符号距离场；

抽取符号距离场的等值面，并在等值面的每个三角形网格上初始化边界粒子，每个三角形网格上布置gA/πr²个粒子；其中，g表示控制粒子密度的参数，A表示每个三角形网格的面积，r表示粒子的半径；

根据粒子的分布来决定每个粒子在固体表面法向和切向的速度；其中，X表示粒子的位置，和分别表示位置X处的符号距离和法向，Q表示一个类高斯核函数，S是整个粒子集合，Xl是粒子X周围第l粒子的位置；通过所述固体表面法向保证了粒子只能在边界表面移动，通过所述切向保证了粒子由密集向稀疏的区域移动；

对粒子的速度进行积分来更新粒子位置；

通过给整体粒子集合S的总位移设置一个粒子数量相关的阈值，来判断是否收敛；当总位移低于阈值时，表明此时迭代过程已经趋于稳定，均匀化的粒子采样已经形成。

进一步的，计算该固体中流体粒子状态的公式包括：

流体与固体边界处流体粒子的密度的计算公式为：

ρ_i＝∑_jm_jW_ij+∑_kΨ_bi(ρ₀)W_ik；

Ψ_bi(ρ₀)＝ρ₀V_bi；

其中，ρ_i表示流体与固体边界处流体粒子i的密度；m_j表示流体粒子i周围第j个流体粒子的质量，W_ij是一个支持域有限的核函数，ρ₀表示流体的静止密度，V_bi表示固体边界粒子bi的体积，k表示固体粒子bi周围的固体粒子；

流体粒子的压力及粘力的计算公式为：

其中，表示固体边界粒子bi对流体粒子i所施加的压力，m_i表示流体粒子i的质量，P_i表示流体粒子i的压强，ρ_i表示流体粒子i的密度，表示核函数W_ibi的梯度向量函数；表示固体边界粒子bi对流体粒子i的粘力，μ表示流体的粘度系数，v_i表示流体粒子i的速度，表示核函数W_ibi的拉普拉斯Laplacian函数；

流体粒子表面张力F_i←j ^cur的计算公式为：

F_i←j ^cur＝-γm_i(n_i-n_j)；

其中，γ表示控制力大小的参数，m_i表示流体粒子i的质量，n_i和n_j分别表示流体粒子i与流体粒子j的法向。

进一步的，所述当所述流体粒子穿透所述固体边界时，对所述流体粒子的速度和位置进行校正包括：

对固体边界区域进行标记，当流体粒子进入标记过的区域时判断流体粒子位置处的符号距离从而得到流体粒子是否穿透边界；

若穿透，则按照下述公式分别对粒子的速度和位置进行校正：

v_i＝v_i-(1+C_r)(v_i·n)n；

X_i＝X_i+(1+D_r)|s|n；

其中，v_i表示流体粒子i的速度，C_r表示控制动量损失的参数，值介于0～1之间；X_i表示流体粒子i的位置，D_r表示调节回弹距离的一个随机数，s表示X_i处的符号距离，n表示边界的法向。

一种三维流固单向耦合的实现系统，该系统包括：

固体边界粒子采样模块，用于将固体边界的三角形网格离散成粒子，并进行边界粒子均匀化采样；

流体粒子状态计算模块，用于计算该固体中流体粒子的状态，包括：流体与固体边界处流体粒子的密度、流体粒子的压力及粘力，以及流体粒子表面张力；

流体粒子信息更新模块，用于并利用一数值积分模式对流体粒子的速度和位置信息进行更新；

穿透粒子的速度和位置校正模块，用于当所述流体粒子穿透所述固体边界时，对所述流体粒子的速度和位置进行校正，从而实现三维流固单向耦合。

进一步的，所述固体边界粒子采样模块包括：

符号距离场建立模块，用于读取固体模型的三角形网格，并建立固体模型的周围符号距离场，具体的：将固体模型中每个三角形网格建立一个包围盒，并为每个包围盒以统一的方式划分网格；计算每个网格点到固体模型的符号距离，符号距离场以一个哈希表来存储，形成一个存在于固体模型周围的窄带宽的符号距离场；

边界粒子初始化模块，用于抽取符号距离场的等值面，并在等值面的每个三角形网格上初始化边界粒子，每个三角形网格上布置gA/πr²个粒子；其中，g表示控制粒子密度的参数，A表示每个三角形网格的面积，r表示粒子的半径；

法向与切向速度决定模块，用于根据粒子的分布来决定每个粒子在固体表面法向和切向的速度；其中，X表示粒子的位置，和分别表示位置X处的符号距离和法向，Q表示一个类高斯核函数，S是整个粒子集合，Xl是粒子X周围第l粒子的位置；通过所述固体表面法向保证了粒子只能在边界表面移动，通过所述切向保证了粒子由密集向稀疏的区域移动；

位置更新模块，用于对粒子的速度进行积分来更新粒子位置；

均匀化采样模块，用于通过给整体粒子集合S的总位移设置一个粒子数量相关的阈值，来判断是否收敛；当总位移低于阈值时，表明此时迭代过程已经趋于稳定，均匀化的粒子采样已经形成。

进一步的，所述流体粒子状态计算模块包括：

密度计算模块，流体与固体边界处流体粒子的密度的计算公式为：

ρ_i＝∑_jm_jW_ij+∑_kΨ_bi(ρ₀)W_ik；

Ψ_bi(ρ₀)＝ρ₀V_bi；

其中，ρ_i表示流体与固体边界处流体粒子i的密度；m_j表示流体粒子i周围第j个粒子的质量，W_ij是一个支持域有限的核函数，ρ₀表示流体的静止密度，V_bi表示固体边界粒子bi的体积，k表示固体粒子bi周围的固体粒子；

压力及粘力计算模块，流体粒子的压力及粘力的计算公式为：

其中，表示固体边界粒子bi对流体粒子i所施加的压力，m_i表示流体粒子i的质量，P_i表示流体粒子i的压强，ρ_i表示流体粒子i的密度，表示核函数W_ibi的梯度向量函数；表示固体边界粒子bi对流体粒子i的粘力，μ表示流体的粘度系数，v_i表示流体粒子i的速度，表示核函数W_ibi的拉普拉斯Laplacian函数；

表面张力计算模块，流体粒子表面张力F_i←j ^cur的计算公式为：

F_i←j ^cur＝-γm_i(n_i-n_j)；

其中，γ表示控制力大小的参数，m_i表示流体粒子i的质量，n_i和n_j分别表示流体粒子i与流体粒子j的法向。

进一步的，所述穿透粒子的速度和位置校正模块包括：

穿透判断模块，用于对固体边界区域进行标记，当流体粒子进入标记过的区域时判断流体粒子位置处的符号距离从而得到流体粒子是否穿透边界；

校正模块，用于穿透时，则按照下述公式分别对粒子的速度和位置进行校正：

v_i＝v_i-(1+C_r)(v_i·n)n；

X_i＝X_i+(1+D_r)|s|n；

其中，v_i表示流体粒子i的速度，C_r表示控制动量损失的参数，值介于0～1之间；X_i表示流体粒子i的位置，D_r表示调节回弹距离的一个随机数，s表示X_i处的符号距离，n表示边界的法向。

由上述本发明提供的技术方案可以看出，对固体边界布置边界粒子，整个计算过程不在流体模拟过程中，相对于模拟过程中动态生成边界粒子的方法而言具有较高的性能提升，同时边界粒子再经过均匀化后可以减小SPH数值模拟计算过程中的误差，增强算法在处理边界时的稳定性；将固体粒子与流体粒子之间的作用力以一种与流体粒子间作用力类似的方式来建模，这样可以将体粒子与固体粒子统一化处理，在成功模拟出压力、粘力和阻止穿透效果的同时大大减小了计算的复杂度；另外，还可有效模拟出表面张力的宏观效果，且不会出现直接模拟聚合力情况下的网状失真现象。

附图说明

为了更清楚地说明本发明实施例的技术方案，下面将对实施例描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域的普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他附图。

图1为本发明实施例一提供的一种三维流固单向耦合的实现方法的流程图；

图2为本发明实施例一提供的一种抽取符号距离场的等值面的示意图；

图3为本发明实施例一提供的一种边界粒子采样均匀化的示意图；

图4为本发明实施例二提供的一种三维流固单向耦合的实现系统的示意图；

图5为本发明实施例提供的一种流体粒子下落后与水箱边界发生碰撞的示意图；

图6为本发明实施例提供的一种流体粒子接触到水箱底部开始的模拟过程的示意图；

图7为本发明实施例提供的一种立方体液体在不施加重力的情况下迅速形变为球体过程的示意图。

具体实施方式

下面结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明的保护范围。

实施例一

图1为本发明实施例一提供的一种三维流固单向耦合的实现方法的流程图。如图1所示，该方法主要包括如下步骤：

步骤11、将固体边界的三角形网格离散成粒子，并进行边界粒子均匀化采样。

具体的实施步骤如下：

1)读取固体模型的三角形网格(主要是将固体模型的不规则三角形边界网格读入内存)，并建立固体模型的周围符号距离场，具体的：将固体模型中每个三角形网格建立一个包围盒，并为每个包围盒以统一的方式划分网格；计算每个网格点到固体模型的符号距离，符号距离场以一个哈希表来存储，形成一个存在于固体模型周围的窄带宽的符号距离场；

2)抽取符号距离场的等值面，示例性的，如图2所示，可以采用Marching Cubes算法来抽取每个网格的等值面(图2中的a)，得到一个均匀化的三角形网格表示(图2中的b)；

3)在等值面的每个三角形网格上初始化边界粒子，每个三角形网格上布置gA/πr²个粒子；其中，g表示控制粒子密度的参数，A表示每个三角形网格的面积，r表示粒子的半径；

4)根据粒子的分布来决定每个粒子在固体表面法向和切向的速度；其中，X表示粒子的位置，和分别表示位置X处的符号距离和法向，S是整个粒子集合，Xl是粒子X周围第l粒子的位置；Q表示一个类高斯核函数，示例性的，可以采用一个6次多项式函数，其中h是核函数的核半径，公式如下所示：

通过所述固体表面法向保证了粒子只能在边界表面移动，通过所述切向保证了粒子由密集向稀疏的区域移动；

5)对粒子的速度进行积分来更新粒子位置；

6)通过给整体粒子集合S的总位移设置一个粒子数量相关的阈值，来判断是否收敛；当总位移低于阈值时，表明此时迭代过程已经趋于稳定，均匀化的粒子采样已经形成。

示例性的，请参见图3所示的边界粒子采样均匀化结果，图3中的a是原始粒子分布，图3中的b为算法收敛后粒子的位置，已经趋于均匀。

步骤12、计算该固体中流体粒子的状态。

所述流体粒子的状态包括：流体与固体边界处流体粒子的密度、流体粒子的压力及粘力，以及流体粒子表面张力。

本发明实施例在计算流体粒子状态时考虑固体边界粒子对流体粒子的作用力，具体的计算公式如下：

流体与固体边界处流体粒子的密度的计算公式为：

ρ_i＝∑_jm_jW_ij ⁺∑_kΨ_bi(ρ₀)W_ik；

Ψ_bi(ρ₀)＝ρ₀V_bi；

其中，ρ_i表示流体与固体边界处流体粒子i的密度；m_j表示流体粒子i周围第j个流体粒子的质量，W_ij是一个支持域有限的核函数，ρ₀表示流体的静止密度，V_bi表示固体边界粒子bi的体积，k表示固体粒子bi周围的固体粒子；

边界处流体粒子的压力及粘力的计算公式为：

其中，表示固体边界粒子bi对流体粒子i所施加的压力，m_i表示流体粒子i的质量，P_i表示流体粒子i的压强，ρ_i表示流体粒子i的密度，表示核函数W_ibi的梯度向量函数；表示固体边界粒子bi对流体粒子i的粘力，μ表示流体的粘度系数，v_i表示流体粒子i的速度，表示核函数W_ibi的拉普拉斯Laplacian函数；

流体粒子表面张力F_i←j ^cur的计算公式为：

F_i←j ^cur＝-γm_i(n_i-n_j)；

其中，γ表示控制力大小的参数，m_i表示流体粒子i的质量，n_i和n_j分别表示流体粒子i与流体粒子j的法向。

步骤13、利用一数值积分模式对流体粒子的速度和位置信息进行更新。

本发明实施例中，可选择一种数值积分模式对流体粒子的加速度和速度做时间积分，从而更新粒子的速度和位置信息，关于SPH计算过程和时间积分模式选择业界已经有成熟的做法，不在赘述。

步骤14、当所述流体粒子穿透所述固体边界时，对所述流体粒子的速度和位置进行校正，从而实现三维流固单向耦合。

具体实现过程如下：

对固体边界区域进行标记，当流体粒子进入标记过的区域时判断流体粒子位置处的符号距离从而得到流体粒子是否穿透边界；

若穿透，则按照下述公式分别对粒子的速度和位置进行校正：

v_i＝v_i-(1+C_r)(v_i·n)n；

X_i＝X_i+(1+D_r)|s|n；

其中，v_i表示流体粒子i的速度，C_r表示控制动量损失的参数，值介于0～1之间；X_i表示流体粒子i的位置，D_r表示调节回弹距离的一个随机数，s表示X_i处的符号距离，n表示边界的法向。

为了便于理解本发明，下面针对流体粒子的状态计算过程，以及流体粒子穿透时的校正过程进行详细介绍。

1、边界流体粒子密度。

在处理流体边界处时，由于流体粒子周围没有足够的粒子来支撑SPH计算而导致边界处粒子密度估计不准，从而导致压强计算不准，这种误差会传播到整个计算域，因此需要对边界处流体粒子的密度做校正。本发明实施例中，考虑流体粒子i周围的边界粒子的积分贡献，下面将具体阐述。

目前，标准SPH方法只考虑流体粒子之间的作用力，其计算密度方法如下式所示：

其中，m_j是流体粒子i周围第j个流体粒子的质量，W是一个支持域有限的核函数。

为了解决边界处流体粒子密度估计偏低问题，我们将边界上的固体粒子(即执行步骤11后获得的固体边界粒子)纳入周围流体粒子的SPH计算的核半径当中，即使得固体粒子对流体粒子的密度计算产生贡献。由于边界粒子的分布不能完全均匀(曲率大区域，边界粒子较集中，曲率小的区域边界粒子较稀疏)，因此，考虑边界粒子的体积V_bi：

其中，m_bi表示固体边界粒子bi的质量，ρ_bi表示固体边界粒子bi的密度，固体边界粒子bi的密度的估计采用标准的SPH计算形式(ρ_bi＝∑_km_bkW_ik)，k表示固体粒子bi周围的固体粒子。

假定边界粒子的质量都相同，体积V_bi的表达式可简化为：

其中，δ_bi＝∑_kW_ik，令Ψ_bi(ρ₀)＝ρ₀V_bi，ρ₀表示流体的静止密度，此时Ψ_bi(ρ₀)在边界粒子集中的区域值较小，而稀疏的区域值较大，综合以上，得到校正边界处流体粒子密度的计算方法如下：

ρ_i＝∑_jm_jW_ij+∑_kΨ_bi(ρ₀)W_ik；

其中，ρ_i表示流体与固体边界处流体粒子i的密度；m_j表示流体粒子i周围第j个流体粒子的质量，W_ij是一个支持域有限的核函数。

2、流体粒子的压力与粘力。

本发明实施例中，将固体边界离散为均匀的粒子，在计算流体粒子的压力与粘力时，考虑其对周围流体粒子影响，下面将分布阐述压力模型和粘力模型。

1)固液间压力模型。标准SPH方法中仅考虑流体粒子之间的作用力，对应的压力计算公式如下：

其中，m_i和m_j分别表示流体粒子i、j的质量，P_i和P_j分别表示流体粒子i、j处的压强，ρ_i和ρ_j分别表示流体粒子i、j的密度，表示核函数W_ij的梯度向量函数。

本发明实施例中，固体边界粒子与流体粒子之间的压力计算采用类似上式方法，由于固体边界采用算法的特性，固体粒子的质量不能单纯地用一个固定值来衡量，我们采用边界流体粒子密度中的Ψ_bi(ρ₀)替代，另外由于对于不可压缩流体而言，密度震荡很小，我们将式子中ρ_iρ_j替代为ρ_i ²，用P_i替代。综上得到固液粒子间的压力(固体粒子对流体粒子所施加的压力)计算公式：

表示固体边界粒子bi对流体粒子i所施加的压力，m_i表示流体粒子i的质量，P_i表示流体粒子i的压强，ρ_i表示流体粒子i的密度，表示核函数W_ibi的梯度向量函数。

2)固液间粘力模型。标准SPH方法中，流体间粘力项用相对速度矢量的Laplacian来衡量，流体粒子间粘力的计算有以下式子：

其中，μ表示流体的粘度系数，v_i和v_i表示流体粒子i、j的速度，表示核函数W_ij的拉普拉斯(Laplacian)函数。

类比流体间粘力为固体对流体的粘力建模，有以下式子：

其中，表示固体边界粒子bi对流体粒子i的粘力，μ表示流体的粘度系数，v_i表示流体粒子i的速度，表示核函数W_ibi的拉普拉斯Laplacian函数；由于只考虑固液单向耦合，固体边界固定不变，因此相对速度就是流体粒子i的速度。

3、流体粒子表面张力。

表面张力从微观层面上看，流体分子和周围的流体分子之间存在相互的吸引力，对于流体内部的分子周围存在均匀的邻居分子，总体上净吸引力为0，处于表面的分子受到的合力不为0而是沿着表面法向指向流体内部。根据Laplace’s law，表面张力的宏观效果是最小化表面区域，使得在没有其他外力作用下液滴呈规则的球状，即表面张力抚平表面曲率大的区域。

我们的表面张力模型基于以下几点：第一，表面粒子曲率与通过颜色域求得的法向的模相关(流体内部法向的模接近0，曲率也接近0，流体表面法向模长较大，曲率也会较大)。第二，表面张力应该使得流体表面趋于光滑，抚平曲率大的区域，即使得粒子法向尽量一致。综合以上两点，考虑给相邻的粒子间施加对称的力F_i←j ^cur：

F_i←j ^cur＝-γm_i(n_i-n_j)；

其中，γ表示控制力大小的参数，m_i表示流体粒子i的质量，n_i和n_j分别表示流体粒子i与流体粒子j的法向。粒子法向通过计算颜色域的梯度得到，颜色域是在计算域内有流体处值为1，其他处值为0的连续的标量，其梯度用SPH方法计算如下式所示：

4、穿透时流体粒子的速度和位置的校正。

为了提高效率，我们并不会对所有流体粒子进行检测，而只会针对那些靠近固体边界的粒子，这里我们采用检测校正的模式来阻止穿透现象的发生。首先，对固体边界区域进行标记，当流体粒子进入标记过的区域时判断流体粒子位置处的符号距离从而得到流体粒子是否穿透边界；若穿透，则按照下述公式分别对粒子的速度和位置进行校正：

v_i＝v_i-(1+C_r)(v_i·n)n；

X_i＝X_i+(1+D_r)|s|n；

其中，v_i表示流体粒子i的速度，C_r表示控制动量损失的参数，值介于0～1之间；X_i表示流体粒子i的位置，D_r表示调节回弹距离的一个随机数，s表示X_i处的符号距离，n表示边界的法向。

本发明实施例通过对固体边界布置边界粒子，整个计算过程不在流体模拟过程中，相对于模拟过程中动态生成边界粒子的方法而言具有较高的性能提升，同时边界粒子再经过均匀化后可以减小SPH数值模拟计算过程中的误差，增强算法在处理边界时的稳定性；将固体粒子与流体粒子之间的作用力以一种与流体粒子间作用力类似的方式来建模，这样可以将体粒子与固体粒子统一化处理，在成功模拟出压力、粘力和阻止穿透效果的同时大大减小了计算的复杂度；另外，还可有效模拟出表面张力的宏观效果，且不会出现直接模拟聚合力情况下的网状失真现象。

实施例二

图4为本发明实施例二提供的一种三维流固单向耦合的实现系统的示意图。如图4所示，该系统主要包括：

固体边界粒子采样模块41，用于将固体边界的三角形网格离散成粒子，并进行边界粒子均匀化采样；

流体粒子状态计算模块42，用于计算该固体中流体粒子的状态，包括：流体与固体边界处流体粒子的密度、流体粒子的压力及粘力，以及流体粒子表面张力；

流体粒子信息更新模块43，用于并利用一数值积分模式对流体粒子的速度和位置信息进行更新；

穿透粒子的速度和位置校正模块44，用于当所述流体粒子穿透所述固体边界时，对所述流体粒子的速度和位置进行校正，从而实现三维流固单向耦合。

进一步，所述固体边界粒子采样模块41包括：

符号距离场建立模块411，用于读取固体模型的三角形网格，并建立固体模型的周围符号距离场，具体的：将固体模型中每个三角形网格建立一个包围盒，并为每个包围盒以统一的方式划分网格；计算每个网格点到固体模型的符号距离，符号距离场以一个哈希表来存储，形成一个存在于固体模型周围的窄带宽的符号距离场；

边界粒子初始化模块412，用于抽取符号距离场的等值面，并在等值面的每个三角形网格上初始化边界粒子，每个三角形网格上布置gA/πr²个粒子；其中，g表示控制粒子密度的参数，A表示每个三角形网格的面积，r表示粒子的半径；

法向与切向速度决定模块413，用于根据粒子的分布来决定每个粒子在固体表面法向和切向的速度；其中，X表示粒子的位置，和分别表示位置X处的符号距离和法向，Q表示一个类高斯核函数，S是整个粒子集合，Xl是粒子X周围第l粒子的位置；通过所述固体表面法向保证了粒子只能在边界表面移动，通过所述切向保证了粒子由密集向稀疏的区域移动；

位置更新模块414，用于对粒子的速度进行积分来更新粒子位置；

均匀化采样模块415，用于通过给整体粒子集合S的总位移设置一个粒子数量相关的阈值，来判断是否收敛；当总位移低于阈值时，表明此时迭代过程已经趋于稳定，均匀化的粒子采样已经形成。

进一步，所述流体粒子状态计算模块42包括：

密度计算模块421，流体与固体边界处流体粒子的密度的计算公式为：

ρ_i＝∑_jm_jW_ij+∑_kΨ_bi(ρ0)W_ik；

Ψ_bi(ρ₀)＝ρ₀V_bi；

其中，ρ_i表示流体与固体边界处流体粒子i的密度；m_j表示流体粒子i周围第j个粒子的质量，W_ij是一个支持域有限的核函数，ρ₀表示流体的静止密度，V_bi表示固体边界粒子bi的体积，k表示固体粒子bi周围的固体粒子；

压力及粘力计算模块422，流体粒子的压力及粘力的计算公式为：

其中，表示固体边界粒子bi对流体粒子i所施加的压力，m_i表示流体粒子i的质量，P_i表示流体粒子i的压强，ρ_i表示流体粒子i的密度，表示核函数W_ibi的梯度向量函数；表示固体边界粒子bi对流体粒子i的粘力，μ表示流体的粘度系数，v_i表示流体粒子i的速度，表示核函数W_ibi的拉普拉斯Laplacian函数；

表面张力计算模块423，流体粒子表面张力F_i←j ^cur的计算公式为：

F_i←j ^cur＝-γm_i(n_i-n_j)；

其中，γ表示控制力大小的参数，m_i表示流体粒子i的质量，n_i和n_j分别表示流体粒子i与流体粒子j的法向。

进一步的，所述穿透粒子的速度和位置校正模块44包括：

穿透判断模块441，用于对固体边界区域进行标记，当流体粒子进入标记过的区域时判断流体粒子位置处的符号距离从而得到流体粒子是否穿透边界；

校正模块442，用于穿透时，则按照下述公式分别对粒子的速度和位置进行校正：

v_i＝v_i-(1+C_r)(v_i·n)n；

X_i＝X_i+(1+D_r)|s|n；

其中，v_i表示流体粒子i的速度，C_r表示控制动量损失的参数，值介于0～1之间；X_i表示流体粒子i的位置，D_r表示调节回弹距离的一个随机数，s表示X_i处的符号距离，n表示边界的法向。

需要说明的是，上述系统中包含的各个功能模块所实现的功能的具体实现方式在前面的各个实施例中已经有详细描述，故在这里不再赘述。

所属领域的技术人员可以清楚地了解到，为描述的方便和简洁，仅以上述各功能模块的划分进行举例说明，实际应用中，可以根据需要而将上述功能分配由不同的功能模块完成，即将系统的内部结构划分成不同的功能模块，以完成以上描述的全部或者部分功能。

另一方面，基于本发明的上述实施例所描述的方案进行了实验。实验环境和相关参数如下：

CPU：Intel Core(TM)i5-2300 2.8GHz

内存：3.49GB

操作系统：Win7Service Pack1 32位

主程序语言：C\C++

开发与调试环境：Microsoft Visual Studio 2010

图形接口：OpenGL

调试和运行：32位的DEBUG版本和32位RELEASE版本

实验名称：水箱液柱下落

实验模拟在一个正方体水箱(体积1*1*1)中，液柱从水箱中间高度处下落的过程。实验中涉及的主要参数如表1所示：

参数名称	流体粒子数	核半径	气体状态常数	粘度系数	静止密度	时间步长
							值	8192	0.045	10.0	7.0	1000	0.005

表1 实验所涉及的主要参数

模拟的平均帧率达到58FPS，可获得实时流畅的模拟效果，图5-图7为模拟过程的关键帧截图：

图5中的a是水箱中流体粒子的初始分布情况，图5中的b为下落后与水箱边界发生碰撞后的效果。

图6中的a-h为液柱接触到水箱底部开始的模拟过程中一系列关键帧截图。

图7中的a-d为立方体液体在不施加重力的情况下迅速形变为球体过程的关键帧截图，此时立方体液体仅受表面张力作用，使得表面趋于平坦，最终形成一个规则球体。

通过以上的实施方式的描述，本领域的技术人员可以清楚地了解到上述实施例可以通过软件实现，也可以借助软件加必要的通用硬件平台的方式来实现。基于这样的理解，上述实施例的技术方案可以以软件产品的形式体现出来，该软件产品可以存储在一个非易失性存储介质(可以是CD-ROM，U盘，移动硬盘等)中，包括若干指令用以使得一台计算机设备(可以是个人计算机，服务器，或者网络设备等)执行本发明各个实施例所述的方法。

以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明披露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应该以权利要求书的保护范围为准。

Claims (6)

1.一种三维流固单向耦合的实现方法，其特征在于，该方法包括：

将固体边界的三角形网格离散成粒子，并进行边界粒子均匀化采样；

计算该固体中流体粒子的状态，包括：流体与固体边界处流体粒子的密度、流体粒子的压力及粘力，以及流体粒子表面张力；

并利用一数值积分模式对流体粒子的速度和位置信息进行更新；

当所述流体粒子穿透所述固体边界时，对所述流体粒子的速度和位置进行校正，从而实现三维流固单向耦合；

其中，所述将固体边界的三角形网格离散成粒子，并进行边界粒子均匀化采样包括：

读取固体模型的三角形网格，并建立固体模型的周围符号距离场，具体的：将固体模型中每个三角形网格建立一个包围盒，并为每个包围盒以统一的方式划分网格；计算每个网格点到固体模型的符号距离，符号距离场以一个哈希表来存储，形成一个存在于固体模型周围的窄带宽的符号距离场；

抽取符号距离场的等值面，并在等值面的每个三角形网格上初始化边界粒子，每个三角形网格上布置gA/πr²个粒子；其中，g表示控制粒子密度的参数，A表示每个三角形网格的面积，r表示粒子的半径；

根据粒子的分布来决定每个粒子在固体表面法向和切向的速度；其中，X表示粒子的位置，和分别表示位置X处的符号距离和法向，Q表示一个类高斯核函数，S是整个粒子集合，Xl是粒子X周围第l粒子的位置；通过所述固体表面法向保证了粒子只能在边界表面移动，通过所述切向保证了粒子由密集向稀疏的区域移动；

对粒子的速度进行积分来更新粒子位置；

通过给整体粒子集合S的总位移设置一个粒子数量相关的阈值，来判断是否收敛；当总位移低于阈值时，表明此时迭代过程已经趋于稳定，均匀化的粒子采样已经形成。

2.根据权利要求1所述的方法，其特征在于，计算该固体中流体粒子状态的公式包括：

流体与固体边界处流体粒子的密度的计算公式为：

ρ_i＝∑_jm_jW_ij+∑_kΨ_bi(ρ₀)W_ik；

Ψ_bi(ρ₀)＝ρ₀V_bi；

其中，ρ_i表示流体与固体边界处流体粒子i的密度；m_j表示流体粒子i周围第j个流体粒子的质量，W_ij是一个支持域有限的核函数，ρ₀表示流体的静止密度，V_bi表示固体边界粒子bi的体积，k表示固体粒子bi周围的固体粒子；

流体粒子的压力及粘力的计算公式为：

其中，表示固体边界粒子bi对流体粒子i所施加的压力，m_i表示流体粒子i的质量，P_i表示流体粒子i的压强，ρ_i表示流体粒子i的密度，表示核函数W_ibi的梯度向量函数；表示固体边界粒子bi对流体粒子i的粘力，μ表示流体的粘度系数，v_i表示流体粒子i的速度，表示核函数W_ibi的拉普拉斯Laplacian函数；

流体粒子表面张力F_i←j ^cur的计算公式为：

F_i←j ^cur＝-γm_i(n_i-n_j)；

其中，γ表示控制力大小的参数，m_i表示流体粒子i的质量，n_i和n_j分别表示流体粒子i与流体粒子j的法向。

3.根据权利要求1所述的方法，其特征在于，所述当所述流体粒子穿透所述固体边界时，对所述流体粒子的速度和位置进行校正包括：

对固体边界区域进行标记，当流体粒子进入标记过的区域时判断流体粒子位置处的符号距离从而得到流体粒子是否穿透边界；

若穿透，则按照下述公式分别对粒子的速度和位置进行校正：

v_i＝v_i-(1+C_r)(v_i·n)n；

X_i＝X_i+(1+D_r)|s|n；

其中，v_i表示流体粒子i的速度，C_r表示控制动量损失的参数，值介于0～1之间；X_i表示流体粒子i的位置，D_r表示调节回弹距离的一个随机数，s表示X_i处的符号距离，n表示边界的法向。

4.一种三维流固单向耦合的实现系统，其特征在于，该系统包括：

固体边界粒子采样模块，用于将固体边界的三角形网格离散成粒子，并进行边界粒子均匀化采样；

流体粒子状态计算模块，用于计算该固体中流体粒子的状态，包括：流体与固体边界处流体粒子的密度、流体粒子的压力及粘力，以及流体粒子表面张力；

流体粒子信息更新模块，用于并利用一数值积分模式对流体粒子的速度和位置信息进行更新；

穿透粒子的速度和位置校正模块，用于当所述流体粒子穿透所述固体边界时，对所述流体粒子的速度和位置进行校正，从而实现三维流固单向耦合；

其中，所述固体边界粒子采样模块包括：

符号距离场建立模块，用于读取固体模型的三角形网格，并建立固体模型的周围符号距离场，具体的：将固体模型中每个三角形网格建立一个包围盒，并为每个包围盒以统一的方式划分网格；计算每个网格点到固体模型的符号距离，符号距离场以一个哈希表来存储，形成一个存在于固体模型周围的窄带宽的符号距离场；

边界粒子初始化模块，用于抽取符号距离场的等值面，并在等值面的每个三角形网格上初始化边界粒子，每个三角形网格上布置gA/πr²个粒子；其中，g表示控制粒子密度的参数，A表示每个三角形网格的面积，r表示粒子的半径；

法向与切向速度决定模块，用于根据粒子的分布来决定每个粒子在固体表面法向和切向的速度；其中，X表示粒子的位置，和分别表示位置X处的符号距离和法向，Q表示一个类高斯核函数，S是整个粒子集合，Xl是粒子X周围第l粒子的位置；通过所述固体表面法向保证了粒子只能在边界表面移动，通过所述切向保证了粒子由密集向稀疏的区域移动；

位置更新模块，用于对粒子的速度进行积分来更新粒子位置；

均匀化采样模块，用于通过给整体粒子集合S的总位移设置一个粒子数量相关的阈值，来判断是否收敛；当总位移低于阈值时，表明此时迭代过程已经趋于稳定，均匀化的粒子采样已经形成。

5.根据权利要求4所述的系统，其特征在于，所述流体粒子状态计算模块包括：

密度计算模块，流体与固体边界处流体粒子的密度的计算公式为：

ρ_i＝∑_jm_jW_ij+∑_kΨ_bi(ρ₀)W_ik；

Ψ_bi(ρ₀)＝ρ₀V_bi；

其中，ρ_i表示流体与固体边界处流体粒子i的密度；m_j表示流体粒子ⁱ周围第j个粒子的质量，W_ij是一个支持域有限的核函数，ρ₀表示流体的静止密度，V_bi表示固体边界粒 子bi的体积，k表示固体粒子bi周围的固体粒子；

压力及粘力计算模块，流体粒子的压力及粘力的计算公式为：

其中，表示固体边界粒子bi对流体粒子i所施加的压力，m_i表示流体粒子i的质量，P_i表示流体粒子i的压强，ρ_i表示流体粒子i的密度，表示核函数W_ibi的梯度向量函数；表示固体边界粒子bi对流体粒子i的粘力，μ表示流体的粘度系数，v_i表示流体粒子i的速度，表示核函数W_ibi的拉普拉斯Laplacian函数；

表面张力计算模块，流体粒子表面张力F_i←j ^cur的计算公式为：

F_i←j ^cur＝-γm_i(n_i-n_j)；

其中，γ表示控制力大小的参数，m_i表示流体粒子i的质量，n_i和n_j分别表示流体粒子i与流体粒子j的法向。

6.根据权利要求4所述的系统，其特征在于，所述穿透粒子的速度和位置校正模块包括：

穿透判断模块，用于对固体边界区域进行标记，当流体粒子进入标记过的区域时判断流体粒子位置处的符号距离从而得到流体粒子是否穿透边界；

校正模块，用于穿透时，则按照下述公式分别对粒子的速度和位置进行校正：

v_i＝v_i-(1+C_r)(v_i·n)n；

X_i＝X_i+(1+D_r)|s|n；

其中，v_i表示流体粒子i的速度，C_r表示控制动量损失的参数，值介于0～1之间；X_i表示流体粒子i的位置，D_r表示调节回弹距离的一个随机数，s表示X_i处的符号距离，n表示边界的法向。