CN110717269B - 一种基于网格和粒子耦合的流体表面细节保护方法 - Google Patents

Abstract

本发明涉及一种基于网格和粒子耦合的流体表面细节保护方法：(1)求解纳维‑斯特克斯方程(Navier‑Stokes，N‑S)采用网格方法，构成主体流体的仿真模型；(2)改进LBM‑VOF方法，并利用改进的VOF‑LBM耦合算法追踪流体的表面；(3)在异常表面网格位置生成粒子，然后通过粒子方法对粒子进行演化；(4)设计网格和粒子的耦合算法，将网格流体和粒子流体集成到同一个场景中，以此保证整个流场的物理守恒以及网格和粒子之间进行合理的物理信息传递；(5)最后利用屏幕空间方法进行真实感渲染，依次通过绘制球体，并计算每个像素点的深度值，进行深度滤波，根据深度值和像素点的位置信息求法向量，以及光照渲染，在GPU上实现一个逼真和实时的流体渲染。

Description

一种基于网格和粒子耦合的流体表面细节保护方法

技术领域

本发明涉及一种基于网格和粒子耦合的流体表面细节保护方法，具体地说是一种基于物理的图形学动画模拟方法，其部分技术涉及到构建高精度的主体流体仿真模型、构建与LBM耦合的流体自由表面追踪模型、粒子的生成与演化、网格与粒子的耦合算法以及力学相关理论。主要应用于各种真实物体三维逼真模拟，特别是应用于娱乐游戏、电影特效等领域。

背景技术

近些年来，计算机图形学领域研究的热点与难点之一就是真实感的场景、动画和现象的模拟，它在制造工程、影视特效、三维动画、电脑游戏与虚拟现实等领域应用广泛。为了模拟出和现实世界相差无几的真实感图形效果，可以基于物理原理和实验数据构建能准确描述物理现象的物理模型。基于物理模型的自然现象的绘制可以模拟出更加逼真的效果，但计算过于复杂，实时性得不到保证。因此，在基于物理的流体仿真中既要减少计算代价，又需要尽可能的保持逼真的视觉效果。

当前计算流体力学中描述流体运动的基本方法有三种：基于粒子的拉格朗日方法、基于网格的欧拉方法以及LBM方法。拉格朗日方法以SPH为代表，是一种基于粒子的方法，着眼于流体质点的运动追踪上，考虑所有质点的运动，进而形成整个流体运动，其优点在于简单直观，容易理解，可模拟流体的运动细节，但是计算量大，难以选取合适的光滑核函数以及在特大变形的情况下容易发生畸变；欧拉方法基于网格，它不具体跟踪某一质点的具体运动过程，而是将流场作为对象。欧拉法有利于分析物理场属性，容易构造液体运动表面，能够处理变形与扭曲问题，但难以表现流体的表面细节；LBM方法是一种基于欧拉的“介观”的流体仿真方法，它忽略了大量分子运动的细节，只保留了分子运动速度的部分统计特征，该方法便于处理复杂边界，算法简单，具有很高的并行度，但依旧是基于网格的描述方法，同时在细节仿真上依然存在不足。

文献1-Koerner C,Thies M,Hofmann T,et al.Lattice Boltzmann Model forFree Surface Flow for Modeling Foaming[J].Journal of Statistical Physics,2005,121(1-2):179-196，针对在表面网格重构分布函数的问题，提出了一种基于动量交换法的表面网格重构方法，但是此方法数值精度低，平衡状态下表面网格在视觉效果上还会发生抖动现象，同时文献1还针对异常表面网格提出了一套特殊的处理方法，一定程度上缓解了该问题，但是难以保证物理量的守恒。文献2-郝爱民,李帅,高阳.一种基于欧拉-拉格朗日耦合方法的流体仿真方法，该方法虽然也将LBM与SPH耦合在一起，但是与本发明有着本质的不同，首先是目标不同，文献2着重于增加流体细节，从而达到丰富的视觉效果，而本发明着重于保护流体仿真的数值精度，从而达到保护流体细节的目的；然后是粒子的生成不同，文献2为了达到增加细节的目的，在表面外生成粒子，而本发明主要是保护仿真的精度，用粒子解决异常网格问题，因此只在异常网格处生成粒子；最后是耦合方式不同，文献2采用PLSM进行表面追踪，并以此为媒介耦合LBM和SPH，本发明采用改进的VOF方法进行表面追踪，数值精度更高。简而言之，文献2为了增加流体细节，可能会破坏原流体的演化规律，而本发明是为了保护原流体演化规律的正常演化。

为了获得更好的真实感模拟，本发明提出了一种基于网格和粒子耦合的流体表面细节保护方法，该方法主要由五部分组成：构建基于LBM的主体流体仿真模型、VOF追踪流体表面，并标记异常界面网格、生成粒子并基于SPH方法进行演化、设计耦合算法耦合网格流体与粒子流体以及真实感渲染，最后实现一个实时、多细节、逼真和稳定的流体仿真。

发明内容

本发明的技术解决问题：克服现有技术的不足，提供一种基于网格和粒子耦合的流体表面细节保护方法，既可以摆脱传统网格法存在的诸多缺陷，又可以很好地保护流体的细节，同时满足了流体仿真对实时性的要求。

本发明的技术解决方案：一种基于网格和粒子耦合的流体表面细节保护方法，步骤包括：

(1)求解纳维-斯特克斯方程(Navier-Stokes，N-S)采用格子莫尔兹曼方程(Lattice Boltzmann Method，LBM)，LBM的碰撞模型采用线性Bhatnagar-Gross-Krook(BGK)碰撞模型，同时采用n维离散空间的m个速度多维离散网格模型(DnQm模型)，构成主体流体的仿真模型，然后经过演化得到每个网格的物理信息；

(2)根据步骤(1)得到的每个网格的物理信息，利用改进的LBM-VOF耦合算法追踪流体的表面，追踪流体表面，重构表面网格分布函数，该步骤得到了主体流体的表面以及表面新的分布函数，并标记流体表面网格中异常表面网格的位置；

(3)在步骤(2)获得异常表面网格位置后，首先对异常表面网格进行处理，即在该位置生成粒子代替原异常表面网格，同时粒子具备原异常表面网格的物理信息，然后通过Smooth Particle Hydrodynamics(SPH)方法对粒子进行演化，得到粒子新的物理信息；

(4)根据步骤(1)所得网格的物理信息和步骤(3)所得粒子的物理信息设计耦合算法，将网格流体和粒子流体集成到同一个场景中，并利用耦合算法将粒子流体转换为网格流体，并重新计算网格的物理信息，包括质量、速度、位置，以此保证整个流场的物理守恒以及网格和粒子之间进行合理的物理信息传递，并得到网格新的物理信息；

(5)将步骤(3)和步骤(4)所得到的粒子和网格的物理信息，利用屏幕空间方法，依次通过绘制球体，并计算每个像素点的深度值，深度滤波，根据深度值和像素点的位置求法向量，以及光照渲染，在GPU上实现一个逼真和实时的流体渲染。

本发明的原理：

(1)传统的LBM方法便于处理复杂边界、算法简单易实现以及并行性高等多个优势构建对主体流体仿真的模型，LBM方法可以分为碰撞和迁移两部分，公式分别为：

f_i(x+c_iδ_t，t+δ_t)＝f′_i(x，t)

其中，i为离散方向的下标，δ_t为单位时间步长，c_i为离散速度的方向向量，c_iδ_t为单位时间的迁移步长，τ为无量纲松弛时间。

(2)传统的LBM-VOF耦合算法通过计算每个流体格子的质量变化来追踪流体的自由表面，将模拟区域大致划分为三类网格：流体网格(记作”F”)，表面网格(记作”I”)和没有流体的网格(记作”G”)，严格意义上I网格必须与F网格相邻，F网格必须被F网格或I网格包围(即F网格不与G网格相邻)，同时F网格与G网格之间的转换必须先转换为I网格。但在实际情况中，会出现I网格不与F网格相邻的情况(本发明将此类I网格称为异常表面网格)，传统解决方案是特殊处理I网格的质量交换以及降低网格类型的转换条件，但是会导致物理量不守恒和数值解的不稳定；针对LBM-VOF方法产生的异常表面网格，本发明不再人工干预这些异常表面网格的演化，本发明采用粒子的方法处理这些异常表面网格，利用SPH方法利于刻画流体细节的优势进行粒子的演化，保证了物理守恒和数值解的稳定性，同时丰富了流体表面的细节。

(3)在主体流体仿真过程中，在I网格处需要考虑G网格对I网格的作用力，而G网格的物理量是未知的，传统的LBM-VOF耦合算法通过动量交换法估计G网格对I网格的作用力，但是误差较大。本发明借鉴固体边界问题的压力非平衡态外推格式方法，将其耦合到LBM-VOF耦合算法中，此方法具有二阶精度，其具体公式为：

f_i(I)＝f_i(G)

其中，表面网格“I”在空网格“G”的i方向，i与互为反向，f_i(I)为表面网格重构后的分布函数在/>方向上的分量，f_i(G)为G网格对I网格作用力在i方向上的分量，/>为G网格i方向平衡态分量，/>为G网格i方向非平衡态分量，ρ_G为G网格密度，一般为常数，u为I网格速度，ρ_I为I网格密度，/>为I网格/>方向上的分量，其中/>的计算不同于压力非平衡态外推格式方法，其仿真结果视觉效果更逼真，流体细节更多。

(4)对异常表面网格进行处理，即在该异常表面网格位置生成粒子代替原异常表面网格，同时粒子具备原异常表面网格的物理信息，然后通过SPH方法对粒子进行演化，得到粒子新的物理信息，具体步骤如下：

(41)计算网格转换为粒子后，粒子的位置信息，具体公式如下：

pos_particle＝pos_grid-(1-ε)·n

其中，pos_grid为表面网格的位置，pos_particle为粒子的位置，ε为I网格的体积积分，n为表面网格的单位法向量。

(42)物理量的传递，具体公式如下：

mass_particle＝mass_grid

v_particle＝v_grid

ρ_particle＝ε·m_grid

其中，mass_grid和v_grid分别为表面网格的质量与速度，mass_particle和υ_particle为所生成粒子的质量与速度，ε为表面网格的体积积分，ρ_particle为粒子的密度。

(43)粒子位置信息和物理信息计算完成后，删除转化为粒子的网格，并通过SPH方法对粒子进行演化，得到粒子新的物理信息；

(5)网格流体与粒子流体耦合算法的设计，主要是为了保证模拟区域内流体在物理规律的约束下演化，具体步骤如下：

(51)计算粒子所在的网格位置，粒子所在的网格位置具体计算公式如下：

其中，pos_particle为粒子演化后的位置信息，pos_grid为粒子所在网格的位置。

(52)如果步骤(51)所求的pos_grid所在位置的类型为表面网格，继续执行步骤(53)。

(53)为保证流体仿真的质量守恒，需要对表面网格的质量进行重新计算，计算公式如下：

m′_grid＝m_particle+m_grid

其中，m′_grid为表面网格重新计算后的质量，m_particle为粒子的质量，m_grid为表面网格重新计算前的质量。

(54)为保证流体仿真的动量守恒，根据步骤(53)计算的m′_grid，重新计算表面网格的速度，新的速度v′_grid要满足以下公式：

m′_gridv′_grid＝m_praticleυ_particle+m_gridv_grid

其中，v′_grid为表面网格重新计算后速度，m′_grid为步骤(53)重新计算后的表面网格质量，m_praticle，v_particle分别为粒子的质量和速度，m_grid，v_grid分别为表面网格重新计算前的质量和速度。

(55)根据步骤(54)所计算的v′_grid，重新计算表面网格的分布函数，计算公式如下：

其中，i为离散方向的下标，τ为无量纲松弛时间，f′_i为表面网格计算后的在i方向上分布函数分量，f_i为表面网格计算前在i方向上分布函数分量，v′_grid为重新计算后速度，为根据v′_grid计算的在i方向上的平衡分布函数。

(56)删除所在位置在表面网格的粒子，并由步骤(51)、(52)、(53)、(54)、(55)更新了表面网格的物理信息，得到新的表面网格；

本发明与现有技术相比的优点在于：

(1)本发明提出一种新的重建表面网格分布函数的方法，从而改进了LBM-VOF方法，使LBM-VOF算法具有更高的鲁棒性以及提高了流体表面追踪的精度。

(2)本发明针对异常表面网格问题，设计了一种新的网格-粒子双向耦合方法，在保证流体仿真的物理守恒的前提下，解决了异常表面网格的问题，保护了流体表面的细节，还通过粒子方法一定程度的丰富了表面细节，最终有效提高了视觉效果的逼真度。

(3)本发明提出的网格-粒子双向耦合方法，其中LBM的空间网格可以作为SPH领域搜索的空间网格，保证了LBM方法和SPH方法的契合度，减少了不必要的计算；同时继承了LBM具有很高并行性的特点，可以通过GPU并行技术加速算法。

附图说明

图1为本发明方法实现流程图；

图2为I网格迁移示意图；

图3为正常网格分布示意图；

图4为异常表面网格处理示意图；

图5为粒子与网格耦合示意图；

图6为I网格绘制位置计算示意图；

图7为本发明方法效果对比图；

图8为最终真实感渲染图。

具体实施方式

下面结合附图实施例对本发明进行具体描述。

本发明一种基于网格和粒子耦合的流体表面细节保护方法，如图1所示，具体实现过程如下：

(1)构建基于LBM的主体流体仿真模型

LBM方法是一种基于欧拉的“介观”的流体仿真方法，它忽略了大量分子运动的细节，只保留了分子运动速度的部分统计特征，该方法便于处理复杂边界，算法简单，具有很高的并行度。Boltzmann方程可以表示为：

其中，f为分子的分布函数，ξ为分子速度，a为分子加速度，Q(f，f)为碰撞项。

在boltzmann方程的基础上用分子分布函数代替分子本身进行演化，其演化方程直接采用格子boltzmann方程，并根据分布函数直接计算流体的密度和速度，其方程被分为碰撞和迁移两个步骤，具体方程分别为：

f′_i(x，t)＝f_i(x，t)+Q(f，f)

f_i(x+c_iδ_t，t+δ_t)＝f′_i(x，t)

其中，i为离散方向的下标，δ_t为单位时间步长，c_i为离散速度的方向向量，c_iδ_t为单位时间的迁移步长。

同时通过BGK碰撞模型将碰撞项线性化，新的碰撞步计算方程如下：

其中，τ为无量纲松弛时间。

根据Boltzmann H定理和DnQm网格离散模型，的计算公式可以近似为：

其中，ω_i为i方向分布函数的权重，e_i为i方向的方向向量，ρ为宏观密度，u为宏观速度。

(2)构建与LBM耦合的流体自由表面追踪

表面追踪是指提取流体的表面信息，表面的精度对仿真的精度有着重要的影响。本发明的流体自由表面追踪可以分为以下三个步骤：

(21)在迁移过程中，表面网格需要估计G网格对I网格的作用力f_i(G)，如图2所示，由于G网格未参与演化，其分布函数是未知的，因此需要估计G网格对I网格的作用力，本发明借鉴固体边界问题的压力非平衡态外推格式方法，将其耦合到LBM-VOF耦合算法中，此方法具有二阶精度，其具体公式为：

f_i(I)＝f_i(G)

其中，表面网格“I”在空网格“G”的i方向，i与互为反向，f_i(I)为表面网格重构后的分布函数在/>方向上的分量，f_i(G)为G网格对I网格作用力在i方向上的分量，/>为G网格i方向平衡态分量，/>为G网格i方向非平衡态分量，ρ_G为G网格密度，一般为常数，u为I网格速度，ρ_I为I网格密度，/>为I网格/>方向上的分量，其中/>的计算不同于压力非平衡态外推格式方法，其仿真结果视觉效果更逼真，流体细节更多。

(22)然后在LBM方法迁移步计算I网格与I网格、F网格与F网格和I网格与F网格之间的质量迁移，其中F网格与F网格和I网格与F网格之间的i方向质量迁移计算公式为：Δm_i＝f_in-f_out，其中f_in为i方向进来的量，f_out为从i方向出去的量，I网格与I网格的质量迁移计算公式为：Δm_i＝ε·(f_in-f_out)，其中ε₂，ε₁为两网格的体积积分即质量与密度之比。

(23)根据迁移后的质量，重新计算体积积分，再根据体积积分更新表面网格的类型，具体判断如下：ε＞1，I网格被充满，变为F网格；ε＜0，I网格被置空，变为G网格；同时要从新分配那些被充满或被置空网格的质量，以及保证F网格不与G网格相邻，正常流体的网格分布如图3所示，F网格被F网格或I网格包围，并且不与G网格相邻，I网格与F网格相邻。然而实际中，会出现不与F网格相邻的I网格，本发明称为异常表面网格，如图4所示，I网格就为异常表面网格。本发明提出了一种全新的方法处理异常表面网格，从粒子的角度出发，用粒子代替异常表面网格进行演化，因此本步骤要标记这些异常表面网格，为下一步做准备。

(3)粒子的生成与演化

(31)粒子的生成：如图4所示，对异常表面网格进行处理，即在该异常表面网格位置生成粒子代替原异常表面网格，同时粒子具备原异常表面网格的物理信息，具体步骤如下：

1)计算网格转换为粒子后，粒子的位置信息，具体公式如下：

pos_particle＝pos_grid-(1-ε)·n

其中，pos_grid为表面网格的位置，pos_particle为粒子的位置，ε为I网格的体积积分，n为表面网格的法向量。

2)物理量的传递，具体公式如下：

mass_particle＝mass_grid

v_particle＝v_grid

ρ_particle＝ε·m_grid

其中，mass_grid和v_grid分别为表面网格的质量与速度，mass_particle和υ_particle为所生成粒子的质量与速度，ε为表面网格的体积积分，ρ_particle为粒子的密度。

3)粒子位置信息和物理信息计算完成后，在pos_particle位置生成新的粒子，并且粒子质量、速度和密度分别为mass_particle、υ_particle和ρ_particle。

(32)粒子的演化：本发明采用光滑粒子流体动力学方法(SPH)，SPH方法是一种利用粒子插值法计算流体力学变量的拉格朗日无网格方法。在SPH方法中，整个流场被离散成一系列粒子，每个粒子都携带了一些物理信息，如质量、速度、位置等。SPH方程的构造通常按照两个关键步进行。第一步为积分表示法，又称“核函数逼近”；第二步为粒子近似法，称“粒子逼近”。从而将积分形式转化成粒子求和的级数形式，从而场内的每个粒子的物理属性可以通过支持域内所有粒子共同插值得到。因此，得到如下表达形式：

其中，范围j为x_i粒子支持域内的粒子，m_j为粒子j的质量，ρ_i是指的所求粒子j的密度。f(x)代表了粒子的任意一种物理属性，W是光滑核函数，h是光滑核函数的影响半径。粒子运动由经典的Navier-Stokes方程(N-S方程)控制，SPH需要对不可压缩流体的N-S方程进行离散化：

其中，若μ是动力粘性系数，v是μ/ρ为运动粘性系数，ρ是密度，p是压强，是速度，F是重力等外力。方程等号右边依次为粘滞力、压力和外力，每个时间步进行各种力的求解，获得粒子的速度，更新位置，驱动粒子运动。

(4)耦合网格流体与粒子流体

网格流体与粒子流体耦合算法的设计，主要是为了保证模拟区域内流体在物理规律的约束下演化，如图5所示，粒子演化后所在网格类型是I网格，符合耦合条件，耦合算法具体步骤如下：

(41)计算粒子所在的网格位置，粒子所在的网格位置具体计算公式如下：

其中，pos_particle为粒子演化后的位置信息，pos_grid为粒子所在网格的位置。

(42)如果步骤(41)所求的pos_grid所在位置的类型为表面网格，继续执行步骤(43)，否则执行步骤(47)。

(43)为保证系统的质量守恒，需要对表面网格的质量进行重新计算，计算公式如下：

m′_grid＝m_particle+m_grid

其中，m′_grid为表面网格重新计算后的质量，m_particle为粒子的质量，m_grid为表面网格重新计算前的质量。

(44)为保证系统的动量守恒，根据步骤(43)计算的m′_grid，重新计算表面网格的速度，新的速度v′_grid要满足以下公式：

m′_gridv′_grid＝m_praticlev_particle+m_gridv_grid

其中，v′_grid为表面网格重新计算后速度，m′_grid为步骤(43)重新计算后的表面网格质量，m_praticle，v_particle分别为粒子的质量和速度，m_grid，v_grid分别为表面网格重新计算前的质量和速度。

(45)根据步骤(44)所计算的v′_grid，重新计算表面网格的分布函数，计算公式如下：

其中，i为离散方向的下标，τ为无量纲松弛时间，f′_i为表面网格计算后的在i方向上分布函数分量，f_i为表面网格计算前在i方向上分布函数分量，v′_grid为重新计算后速度，为根据v′_grid计算的在i方向上的平衡分布函数。

(46)删除所在位置在表面网格的粒子。

(47)剩下的网格和粒子继续演化，然后重复步骤(41)、(42)、(43)、(44)、(45)、(46)、(47)。

(5)真实感渲染

本发明采用屏幕空间法渲染流体。将步骤(4)所得到的流体信息通过绘制精灵点，深度滤波，根据深度值和UV值求法向量等一系列步骤在GPU上实现一个逼真和实时的流体渲染。其中精灵点的半径与体积积分相关联，即：其中k与空间的维度有关系，二维空间为/>三维空间为/>如图6所示，为了避免仿真结果在视觉上流体的不连续感，绘制I网格的精灵点时，需要重新计算位置，使流体在空间上和视觉上更连续，具体计算公式如下：

pos′_new＝pos+(1-ε)·e_i

其中，pos为I网格的原位置，pos′_new为，ε为I网格的体积积分，e_i表示F网格在I网格的i方向的方向向量。

图7展示了传统LBM-VOF方法和本发明的方法在经过相同时间演化后的效果对比图，左图为传统LBM-VOF方法，可以看出相较与本发明丢失了很多细节，右图为本发明的方法，其中画圈内(不仅限于圈内)即为使用本发明所产生的细节，可以看出本发明通过网格-粒子耦合方法不仅保证了网格的仿真精度，保护了流体表面的细节，还通过粒子方法增加了表面的细节；图8为利用本发明方法进行流体仿真，再利用屏幕空间法真实感绘制的效果图，由于算法核心计算均在GPU上完成，因此帧率也达到105FPS(网格大小：40X40X40，GPU：NVIDIA Quadro K5200)，完全满足流体仿真对实时性的要求。

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

Claims (1)

1.一种基于网格和粒子耦合的流体表面细节保护方法，其特征在于，步骤包括：

(1)采用格子莫尔兹曼方程LBM(Lattice Boltzmann Method，LBM)求解纳维-斯特克斯方程(Navier-Stokes，N-S)，LBM的碰撞模型采用线性Bhatnagar-Gross-Krook(BGK)碰撞模型，同时采用n维离散空间的m个速度多维离散网格模型(DnQm模型)，构成主体流体的仿真模型，然后经过演化得到每个网格的物理信息；

(2)根据步骤(1)得到的每个网格的物理信息，利用改进的LBM-VOF耦合算法追踪流体表面，并重构流体表面网格的分布函数，得到主体流体的表面以及表面新的分布函数，并标记流体表面网格中异常表面网格位置；

(3)在步骤(2)获得异常表面网格位置后，首先对异常表面网格进行处理，即在该异常表面网格位置生成粒子代替原异常表面网格，同时粒子具备原异常表面网格的物理信息，然后通过Smooth Particle Hydrodynamics(SPH)方法对粒子进行演化，得到粒子新的物理信息；

(4)根据步骤(1)所得网格的物理信息和步骤(3)所得粒子的物理信息设计耦合算法，将网格流体和粒子流体集成到同一个场景中，并利用耦合算法将粒子流体转换为网格流体，并重新计算网格的物理信息，物理信息包括质量、速度、位置，保证整个流场的物理守恒以及网格和粒子之间进行合理的物理信息传递，并得到网格新的物理信息；

(5)将步骤(3)得到的粒子新的物理信息和步骤(4)网格新的物理信息，利用屏幕空间方法，依次通过绘制球体，并计算每个像素点的深度值，进行深度滤波，根据深度值和像素点的位置信息求法向量，以及光照渲染，在GPU上实现一个逼真和实时的流体渲染；

所述步骤(2)中，重构流体表面网格的分布函数f_i(I)，具体公式如下：

f_i(I)＝f_i(G)

其中，表面网格“I”在空网格“G”的i方向，i与i互为反向，f_i(I)为表面网格重构后的分布函数在i方向上的分量，f_i(G)为G网格对I网格作用力在i方向上的分量，为G网格i方向平衡态分量，/>为G网格i方向非平衡态分量，ρ_G为G网格密度，为常数，u为I网格速度，ρ_I为I网格密度，f_i(I)为I网格i方向上的分量；

所述步骤(3)中异常表面网格的处理方法是用粒子代替异常表面网格，网格转换为粒子的具体公式如下：

(31)计算网格转换为粒子后，粒子的位置信息，具体公式如下：

pos_particle＝pos_grid-(1-ε)·n

其中，pos_grid为表面网格的位置，pos_particle为粒子的位置，ε为I网格的体积积分，n为表面网格的法向量；

(32)物理量的传递，具体公式如下：

mass_particle＝mass_grid

v_particle＝v_grid

ρ_particle＝ε·m_grid

其中，mass_grid和v_grid分别为表面网格的质量与速度，mass_particle和v_particle为所生成粒子的质量与速度，ε为表面网格的体积积分，ρ_particle为粒子的密度；

(33)粒子位置信息和物理信息计算完成后，删除转化为粒子的网格，并通过SPH方法对粒子进行演化，得到粒子新的物理信息；

所述步骤(4)中将粒子流体转换为网格流体的具体步骤如下：

(41)计算粒子所在的网格位置，粒子所在的网格位置具体计算公式如下：

其中，pos_particle为粒子演化后的位置信息，pos_grid为粒子所在网格的位置；

(42)如果步骤(41)所求的pos_grid所在位置的类型为表面网格，继续执行步骤(43)，否则执行步骤(46)；

(43)为保证流体仿真的质量守恒，对表面网格的质量进行重新计算，计算公式如下：

m′_grid＝m_particle+m_grid

其中，m′_grid为表面网格重新计算后的质量，m_particle为粒子的质量，m_grid为表面网格重新计算前的质量；

(44)为保证流体仿真的动量守恒，根据步骤(43)计算的m′_grid，重新计算表面网格的速度，新的速度v′_grid满足以下公式：

m′_gridv′_grid＝m_praticlev_particle+m_gridv_grid

其中，v′_grid为表面网格重新计算后速度，m′_grid为步骤(43)重新计算后的表面网格质量，m_praticle，v_particle分别为粒子的质量和速度，m_grid，v_grid分别为表面网格重新计算前的质量和速度；

(45)根据步骤(44)所计算的v′_grid，重新计算表面网格的分布函数，计算公式如下：

其中，i为离散方向的下标，τ为无量纲松弛时间，f′_i为表面网格计算后的在i方向上分布函数分量，f_i为表面网格计算前在i方向上分布函数分量，v′_grid为重新计算后速度，为根据v′_grid计算的在i方向上的平衡分布函数；

(46)删除所在位置在表面网格的粒子，并由步骤(41)、(42)、(43)、(44)、(45)更新了表面网格的物理信息，得到新的表面网格。