基于物理的不可压缩流体模拟技术综述①

吴德阳 唐 勇 刘浩阳 刘培艺 汪国平**

(*燕山大学信息科学与工程学院 秦皇岛 066004)

(**北京大学信息科学技术学院 北京 100871)

(***河北省计算机虚拟技术与系统集成重点实验室 秦皇岛 066004)

(****北京市虚拟仿真与可视化工程研究中心 北京 100871)

0 引 言

近年来，随着数字信息、图形技术的快速发展，虚拟现实作为一种创建和体验虚拟世界的仿真技术，能够模拟出与自然景象极其逼真的效果，使人如同身临其境，该技术已经被广泛应用于消防、油画生成、游戏场景、影视特效等领域[1-3]。而流体模拟作为计算机图形学的一个研究分支，可以丰富虚拟场景，增强现实感、沉浸感，因此对于流体仿真的研究具有重要意义。随着计算机硬件的计算能力不断提升，流体仿真的计算速度得到了很大改善，但对于实时仿真还存在很大的提升空间。通常，许多研究者通过牺牲部分仿真效果来满足实时性的要求，但游戏、影视特效等领域对实时性和真实感两者都具有很高的要求，因此，如何快速、真实地模拟不可压缩流体仍是计算机图形学领域面临的挑战。

在单相的流体模拟中，求解不可压缩性一直是流体模拟的研究热点，虽然近年来随着对流体模拟算法的不断改进，取得了一些新颖的研究成果[4,5]，但由于受表面张力[6]、边界和空气压强等因素的影响，对于求解边界处的不可压缩性还存在一定的困难。而在真实感上，研究者为了达到更精细的效果，常常将网格进行细化或增加粒子数，虽然细节效果更明显，但与之带来的后果是增加了算法的计算耗时。为此，有些文献则利用GPU的并行计算能力来提高算法的实时性能，如文献[7,8]利用GPU对流体模拟进行加速的策略。在求解不可压缩流体的基础上，流体与固体的双向耦合成为了流固耦合的难点，相比于单种流体运动的模拟，流固耦合除了需要求解流体的不可压之外，还需要考虑流-固界面穿透问题。两相流或者多相流之间耦合的主要难点在于不同相之间的交界面特性，以及不同流体之间的混合、两种化学物质参与反应所产生的效果。多相流界面主要分为液-固界面、气-液界面、液-液界面以及气-液-固界面4种类型，不同相之间构成的界面也有所不同，如液相和固相构成的界面大多为尖锐性界面，而液相与液相构成的界面大多为弥散界面。目前处理尖锐性多相流界面的方法主要有有限体积法[9]和粒子水平集方法[10]，而描述弥散界面的方法主要包含相场法[11]和格子波尔兹曼法( lattice Boltzmann method，LBM )[12]等。目前国内外已有部分文献对流体仿真技术进行了总结[13-15]，为了与这些综述文献区分开，本文主要围绕2003年至2019年基于物理的不可压缩流体仿真算法进行梳理与总结，并通过比较分析，给出了当前流体仿真技术所面临的问题和挑战。

本文第1节介绍了不可压流体纳维-斯托克斯方程(Navier-Stokes equation，N-S)的求解过程。第2节介绍了基于物理的流体仿真方法，分别为Lagrange粒子的方法、基于Eular网格的方法以及基于Lagrange粒子-Eular网格的混合方法。第3节对第2节中的物理仿真方法进行对比分析和评价。第4节提出现有研究存在的不足并对未来值得关注的研究方向进行初步探讨。第5节总结全文。

1 不可压流体的N-S方程求解

不可压缩流体的运动通常由著名的N-S方程来描述，主要由质量守恒方程、动量守恒方程、能量守恒方程3大方程组成。在不可压流体仿真中，为了简化计算，通常将密度设为恒定值并忽略温度的影响，其基本形式如下[16]：

(1)

质量守恒：▽·u=0

(2)

其中，k▽2u为粘力项，(u·▽)u为平流项，▽p为压力项，f为外力项，式(2)为流体的不可压条件，N-S方程求解流程如图1所示。

图1 N-S方程求解流程图

具体求解步骤如下：(1)初始化速度场、压力场等。(2)施加外力，主要为重力。(3)对速度进行平流操作，主要有前向平流和后向平流。前向平流根据当前时刻的中心速度，向前计算下一时刻的位置，并根据下一时刻的物理量通过双线性插值操作更新周围4个网格的物理量，最后将下一时刻的物理量作为当前时刻的物理量。后向平流根据当前的中心速度，向后回溯Δt时刻的位置，然后根据Δt时刻位置的物理量进行双线性插值更新周围4个网格的物理量。(4) 进行压力投影，使流体满足不可压条件。(5) 添加粘力。(6) 添加步骤(4)得到的压强更新速度。

2 基于物理的不可压流体模拟方法

基于物理的不可压流体仿真方法主要由3大类构成，分别为Lagrange粒子法、Eular网格法和Lagrange-Eular混合法。Lagrange粒子法主要以光滑粒子动力学(smoothed particle hydrodynamics,SPH)法、基于位置的流体法(position based fluid,PBF)和LBM法为代表；Eular网格主要涉及高度场网格、水平集方法和相场法；Lagrange-Eular混合法主要以质点网格法(particle in cell，PIC)和流体隐式粒子法(fluid implicit particle，FLIP)为代表，不可压流体仿真方法框架如图2所示。

图2 流体仿真方法框架

2.1 Lagrange粒子法

Lagrange粒子法的基本思想是将连续的流体离散成一堆携带物理量(密度、速度、质量等)的粒子，然后通过计算所有粒子之间的相互作用力来描述流体质点随时间变化的规律。本节主要介绍光滑流体动力学和基于位置的流体模拟方法以及这2种方法的应用。

2.1.1 SPH法

在Lagrange粒子法中，SPH是流体模拟方法最为典型的方法之一。由于该方法实现简单、且能处理拓扑变化较大的边界效果，如喷溅现象等，因此被广泛应用于单种流体模拟[17]、流固耦合[18]、多相流[19]等，其基本模型如图3所示。

图3 SPH粒子模型

SPH的基本思想是通过在位置x处对光滑核半径r范围内的n个粒子累加求和，得到该处的物理属性的累积量[20]：

(3)

模拟的大致步骤如下：(1) 初始化粒子信息，如位置等。(2) 利用式(3)求出粒子的密度ρ。(3) 根据上一步骤得到的粒子密度ρ计算粒子的压强进而推导出由半径r范围内邻居粒子对该粒子施加的压力F。(4) 根据步骤(3)得到的压力F利用牛顿第二定律求出粒子的加速度a。(5) 利用粒子的加速度a求出粒子的速度。(6) 根据速度更新粒子的位置。

(1) 单种流体的仿真

最初的SPH由Stam[20]引入计算机图形学，用于模拟火焰等现象。随后Müller等人[21]将SPH用于流体模拟，模拟的流体效果相比于Eular网格法的真实感更强，其能模拟水珠飞溅的效果，如图4所示，但需要离线渲染，且该方法采用理想气体关联压力和密度，由于与实际的流体压力和密度不同，理想气体关联压强会导致模拟的流体存在较高的压缩性。

图4 SPH粒子法仿真效果[21]

微可压SPH法(weakly compressible SPH，WCSPH)[22]，针对文献[21]的高压缩性问题，Becker等人[22]引入一种新的表面张力模型，可以避免控制方程中导数项的计算，同时为了防止数值耗散，在模型中添加了人工粘度来提高数值的稳定性。WCSPH很好地改善了传统SPH的不可压问题，但需要用较小的时间步长来维持模拟的稳定性。

预测校正SPH法(predictive corrective incompressible SPH，PCISPH)[23]，为了克服WCSPH对时间步长的敏感性问题，Solenthaler等人[23]通过预测-矫正粒子的位置和密度，减少上一帧与下一帧之间的密度差，并关联密度与压强的关系，实现不可压缩性，相比WCSPH方案，其时间步长更长，数值求解更稳定。

局部泊松SPH(local Poisson SPH， LPSPH)[24]，He等人[24]为了弥补WCSPH大密度误差和PCISPH全局求解泊松方程耗时的不足，利用局部作用机理，将泊松方程转换成积分形式，然后通过离散化将积分方程变换为局部压力积分域，减少求解泊松方程的时间，同时为了减少密度误差，将局部压力集成到预测矫正器中，可以避免由于密度波动引起的不稳定问题。

隐式不可压SPH法(implicit incompressible SPH， IISPH)[25]，与WCSPH和PCISPH相比，LPSPH在时间步长和收敛速度上具有很大的改进，但LPSPH受模拟场景规模的限制，只能用于小规模流体的模拟。为此，Ihmsen等人[25]通过离散化压力泊松方程，并在最终的速度更新中使用压力项来表示，可以有效地改进大场景下的时间步长，且能使密度差降低至0.01%。He等人[26]提出一种基于SPH的自由表面流中小尺度薄片特征的鲁棒模拟方法，利用扩散模型将自由表面建模成非常薄的区域，即在界面区域内的粒子表示成1，界面外则表示为0，界面处则是从0到1光滑过渡，这种建模方法可以在界面粒子较少的情况下，仍然可以鲁棒地表示流体拉伸过程中的薄片特征。Bender等人[27]提出一种无散度的SPH方法(divergence-free SPH，DFSPH)，利用压强对密度进行 2 次调整， 使流体实现不可压缩性。Band等人[28]将偏微分方程离散化引入IISPH的压力求解方程中，用于求解流体的压力加速度，可以有效地解决IISPH的压力震荡问题。Li等人[29]利用Bender等人[27]提出的DFSPH对流体的不可压缩性进行求解N-S方程，通过求解密度约束来修正粒子的位置，以保证粒子密度在相对恒定状态。

(2) 流固耦合仿真

由于SPH在处理边界时具有较大的灵活性，该方法除了用于模拟单种流体的运动效果，也被广泛用于流固耦合的模拟。如Shao等人[30]提出一种流固耦合算法，在固体边界处给流体粒子施加一个相反的作用力，防止流体粒子穿透固体边界，解决流固界面的穿透问题，如图5所示，该方法能实现流-固双向耦合、湍流效果，但由于在边界施加这种排斥力并不属于物理力，因此在模拟时常常发生流体粒子堆积现象。Gao等人[31]使用基于位置的动力学求解的位置约束和考虑固体边界粒子对流体粒子的相对贡献，有效解决了粒子堆积的问题。Gissler等人[32]将SPH用于模拟流体与刚体之间的耦合，在该界面中，流体压力求解器每次迭代都会更新刚体颗粒的速度，与SPH相比可以支持更大的时间步长，减少了模拟流体所需的计算时间。

图5 SPH粒子法流固耦合效果[30]

通过对上述的SPH及其演变方法分析发现，SPH在流体模拟中具有如下优势：具有较大的灵活性，能够模拟不同类型流体，如粘性流体、不可压流体与流体之间的交互效果； SPH只需要计算邻域内局部粒子的信息，不需要通过求解全局的N-S方程即可得到流体的运动，因此在小规模流体中具有较高的计算效率。

2.1.2 PBF法

虽然SPH方法具有简单和易于实现的优势，但由于SPH对邻居粒子的密度较为敏感，通过直接求解压强来强制不可压缩性是比较耗时的，因此在粒子较多时无法达到实时的效果。然而在游戏和电影等场景中，往往对实时性要求较高，而且大多数算法通常需要设置小的时间步长来维持系统的稳定性，对于大规模的场景无法达到大时间步长的模拟，为此需要一种既能保证真实感又能达到实时的建模方法。

Müller等人[33]提出了一种基于牛顿第二运动定律的基于位置的动力学方法(position based dynamics，PBD)，基本思想是通过携带质量、位置、速度等物理量的N个粒子和M个约束来表示一个物体，并通过欧拉积分预测粒子的位置，同时将内力建模成一种约束，然后通过约束投影到最终的位置。由于采用显式积分求解约束方程，因此可以有效地消除模拟过程中的不稳定性，但该方法只用于布料的运动建模。随后，Macklin等人[34]基于PBD的思想，提出一种基于位置的流体方法PBF，为了强制不可压缩性，PBF通过使用状态方程建立当前粒子位置与邻居粒子位置之间的密度约束函数：

(4)

其中，ρ0为初始密度，ρi则通过SPH的密度估算器[22]求得，因此避免了SPH先计算力再由力计算粒子位置所带来的拉伸不稳定问题，如图6所示。Zhang等人[35]通过对输入的模型进行采样并自动生成控制粒子，用于控制流体生成的形状，该方法在传统的PBF上增加了流体形状的可控性，能够使工程师使用少量的控制参数去生成与目标形状的流体效果，而且能够满足实时性的要求。

图6 PBF粒子法仿真效果[34]

常规的PBF在模拟包含大量粒子的流体时，收敛速度慢的主要原因在于PBF为了求解流体的不可压缩性时，需要通过迭代不断地调整流体粒子的位置。为此，Köster等人[36]提出一种自适应PBF方法，根据每个粒子的水平细节信息自适应地调整求解器的迭代次数，以保证稳定性与运算效率的平衡。

PBF的优势在于其能够实时模拟流体运动，同时还能允许较大的时间步长，解决了传统SPH方法不能实时模拟和时间步长的限制。

2.1.3 LBM法

LBM法是一种用于求解N-S方程的新型数值方法，其基本思想是使用随机运动的微观粒子表示流体中的分子或宏观流体的一部分，并通过粒子流动和粒子碰撞完成动量的交换，LBM的基本模型如下[37,38]：

(5)

其中，式(5)的左边表示LBM的streaming操作，右边表示collision操作。虽然Eular网格法和Lagrange粒子在模拟流体的宏观运动上具有很好的灵活性和准确性，但对于流体的微观运动却很难捕获。相反，LBM方法则通过微观分子来描述流体的宏观运动状态，避免了传统数值方法的复杂性和精度问题，因此，近年来LBM方法在游戏、工程和物理仿真领域都有广泛的应用。在模拟大场景时，通常使用基于浅水方程的LBM，如Judice等人[39]利用LBM方法用于模拟游戏的流体动画，并将LBM应用于GPU，因此，该方法能够实时模拟大规模的流体动画场景，而LBM的简单高效特性在设计游戏场景中具有很好的优势。Bauza等人[40]提出一种基于LBM的实时流体交互方法，该方法结合浅水方程实现动画与人的实时交互。然而在工程应用中，传统的LBM在模拟多相流效果时，存在密度比低、雷诺数小的问题，为此许多研究者改进了粒子概率分布函数f，如Li等人[41]与Wu等人[42]利用LBM模拟多相流效果，在蒸气点附近利用平均插补的方法校正概率分布函数f，能够模拟密度比大于500的液滴运动和液滴震荡效果。

基于上述分析，LBM方法具有如下特点：

• 计算局部性好，易于扩展到GPU并行框架。

• LBM具有程序易实施、计算速度快的特点。

• 由于LBM基于微观理论，因此能够很好地描述多相流及复杂的边界效果。

2.2 Eular网格法

基于Lagrange粒子法的流体模拟方法虽然具有灵活、方便的优势，并且SPH经过近几年的不断改进，已经在效果和效率上有很大的提升。与Lagrange粒子法不同，基于Eular网格的流体模拟方法的核心思想是将模拟的流体固定在网格上，然后通过计算流过网格节点的物理量来描述流体的运动，在整个仿真过程中，网格不随流体的运动而运动。

2.2.1 高度场网格法

高度场网格法是基于浅水方程的一种流体模拟方法，其基本思想是忽略垂直方向的动力学问题，将流体产生的水波高度值存储于二维的网格中，进而降低求解N-S方程和数值解的复杂度，浅水方程的基本形式如下：

(6)

(7)

其中，u为速度，h为水体的深度，s为地形高度。从式(6)和(7)可以看出，水波的流动只与2个参数有关，一个是高度信息，另一个则是流体的速度，并且这2个公式只在x方向上进行求解，去除了水深方向的信息，简化了计算的复杂度，这也是基于高度场网格的主要优势。

Layton等人[43]提出一种用于产生水波的高效浅水方程方法，该方法通过非线性浅水方程来预测随时间变化的水波高度，实时地模拟了水波效果。由于该方法使用Semi-Lagrange方法[44]求解浅水方程，因此可以达到无条件稳定的大时间步长。随后，Wang等人[45]将该方法用于求解固体与水波之间的双向耦合，能够模拟表面波纹、水滴、流体与固体的耦合效果，但无法模拟波浪运动过程中的破碎效果。为此，Thurey等人[46]提出一种基于浅水方程的实时破碎波浪仿真方法，首先用浅水方程模拟高度场网格，其次，为了产生波浪的破碎效果，通过检测高度场中的陡峭波区域并用线段进行标记，然后在波峰处产生粒子，用于模拟波浪破碎的效果，如图7所示。

图7 高度场网格法仿真效果[46]

基于Thurey等人[46]的方法存在粒子伪影问题，Chentanez等人[47]直接在高度场网格产生喷雾粒子、飞溅粒子和泡沫粒子，并与网格进行信息传递，因此产生的破碎波现象更符合物理运动，但所描述的波长不得小于网格的大小，否则会导致生成的破碎波丢失部分细节。Nielsen等人[48]使用基于FFT的波浪模拟来生成高分辨率的纹理，用于替换渲染网格，产生的波浪细节效果比文献[47]的更逼真，且实时性更好。

为了进一步增强大规模场景中的流体细节效果，Chentanez等人[49]通过耦合3D Eular网格、粒子法和高度场，将高度场网格耦合在3D Eular网格中用于产生流体运动过程中的水波，而在3D网格表面区域添加流体在飞溅过程中所需要的细节粒子，该方法可以快速地模拟大规模场景下的波浪破碎效果。邵绪强等人[50]基于浅水方程提出一种大规模实时的、稳定的Eular网格法，对产生毛刺现象的高度场值进行平滑处理，而针对细节部分，其采用的是将高度场与SPH粒子进行结合，并通过隔点采样的方法降低计算复杂度，相比文献[49]的细节效果更加真实。

综上所述，基于高度场网格法的流体模拟方法相比粒子法在实时性能上具有较大的提升优势，由于高度场网格是基于浅水方程[51]，同时去除了垂直方向的信息，计算效率较高，因此特别适合模拟大规模场景。

2.2.2 界面模型法

界面捕获方法最常用的主要有流体体积法(volume of fluid，VOF)、水平集方法(level set，LS)以及相场法(phase field,PF)，VOF法通常使用流体体积函数F计算网格单元中流体的体积分数，然而在一个两相流系统中，每个单元格都有如下3种可能：

当F=0时，表示该网格单元被流体1所占据。

当F=1时，表示该网格单元被流体2所占据。

当0

并通过控制方程来跟踪两相界面变化：

(8)

其中u为速度，由于采用较为精确的数值方法计算网格单元的流体体积分数，因此可以保证质量守恒的特性，但由于VOF定义的界面是不连续的，因此很难精确地计算界面的曲率和法线等物理信息。

(1) LS法

另一种界面捕获方法LS法，将界面表示成连续光滑的零厚度的水平集函数φ，然后定义界面内侧的流体φ<0，在界面外侧的流体φ>0，在界面层上的φ=0，LS模型如下:

(9)

从LS模型可以看出，LS方法定义的界面是光滑的零水平集函数，这种光滑特性可以准确地计算界面的曲率和法线等物理信息。然而在水平集方法中最大的挑战是质量守恒问题，由于随着界面的演化，界面层的质量不断发生变化，为了保证界面在演化过程中质量守恒，研究者们提出了许多守恒的水平集算法用于减少界面在演化过程中的质量丢失。如Olsson等人[52]提出一种守恒的LS方法，该方法使用守恒的中间步对水平集函数进行平流，可以保证界面的轮廓和厚度是恒定的，不足的是该方法只用于2维的数值验证。Losasso等人[53]提出一种SPH和LS双向耦合的流体模拟方法，通过耦合SPH和LS法用于建模流体的流动，模拟的流体更符合实际的物理运动，但该方法在液体表面生成虚拟粒子，增加了系统开销。基于文献[52]的工作，文献[54,55]通过添加额外项修改重新初始化水平集函数来保持尖锐界面的质量守恒。文献[56]提出新的Heaviside函数求解水平集来保证质量守恒，虽然通过重新初始化水平集函数可以有效解决质量不守恒问题，但会降低算法的收敛速度。

在上述的改进算法中，大都关注解决LS方法的质量守恒问题，而在常见的流体现象中，由于多相流在运动过程中，界面的厚度和形状不断发生变化，因此如何控制这些界面特性也是LS方法所需要解决的问题。基于此，Santos等人[57]利用区域水平集法跟踪多相流的界面，并通过局部矫正来最小化总体积误差。相比于原有的LS方法能够处理更多不同流体间的界面，并且在界面层使用空气或者水进行填充用于控制界面的厚度，这种方式会造成一定的人工伪影现象。Balczar等人[58]针对流体变形问题提出一种新的多标记水平集方法，利用多个水平集函数φi表示混合流体中每个子区域Ωi:Ωd={Ω1, Ω2,…, Ωn}，其中Ωn表示被划分成区域个数，因此可以在相同的体积中控制不同的界面运动，同时将表面张力转换为一种体积力，避免了传统水平集方法的数值凝结问题。

(2) VOF-LS耦合方法

由于LS方法能够很好地捕获界面的光滑特性，VOF方法可以保证质量守恒，将两者耦合一起可以同时获得LS方法和VOF方法的优势，为此，近些年来，许多研究者将两者进行耦合用于处理更加复杂的界面演化现象。如文献[59,60]提出一种耦合水平集合流体体积法(coupled level set and volume of fluid，CLSVOF)方案，该类方法主要通过水平集函数计算界面的表面张力和曲率等，CLSVOF方法能够保证质量守恒的情况下准确地计算界面的曲率和法向等信息。文献[61,62]提出一种流体体积与水平集耦合法(volume-of-fluid-level-set，VOSET)，该方法直接利用几何方法计算水平集界面，可以有效解决界面因震荡引起的不连续问题。由于具有质量守恒和计算精确的优点，因此被广泛用于捕获不同多相流界面。

水平集方法及其改进的方案虽然能保证质量守恒和准确计算曲率等信息，但是该方法仍存在以下2个局限性：

• 由于LS方法定义多相流界面是光滑的零水平集函数，界面的厚度为零，通常只被用于两相流的模拟，且该方法难以表现出多相流混合过程中的扩散特性。

• 为了更加精确地表示界面的零厚度特性，LS方法通常通过网格细化来表示界面的细节特征度，因此增加了计算时间和复杂度。

(3) 相场法

为了解决水平集法界面扩散特性的局限性，Shen等人[63]提出一种用于模拟多相流扩散运动的相场模型，由于界面的厚度是一层薄而非零的过渡区域，因此在模拟过程中，这种扩散界面可以捕获更多的流体运动细节。描述弥散界面的方法主要有2种：一种是基于守恒场的相场方法(cahn-hilliadrd,CH)，另一种是基于非守恒场的相场法(allen-cahn,AC)。基于CH相场模型一般用扩散动力学表示相场浓度随空间和时间的变化，通常用于描述相与相之间演化过程中的守恒场变量，如两相流中的混合流体中的浓度场(也称为质量分数)等，起初用于材料领域的枝晶、共晶和包晶的生长[64]，后被广泛用于多相流领域。CH相场模型的控制方程为

(10)

(11)

式(10)左边第1项表示相场浓度随时间的变化，第2项表示多相流中平流对相场φ的影响，第3项表示扩散项。其中u为速度场；M是有关相场的迁移率，通常只与相场浓度相关，如文献[65]的M=1-φ2，或者为常数，即M=1，关于迁移的设置可以参考文献[66]；μ为化学势， 当达到平衡时，化学势必须满足μ=0。相场模型用于控制多相流中交界面的创建和演化，而不可压流体的流动则由N-S方程控制，将相场方程和N-S方程耦合一起得到著名的“模型H”[67]。

He等人[68]基于相场理论提出一种用于界面控制的流体混合方法，该方法主要包含2部分：第1部分是利用网格法求解N-S方程，获得流体平流的速度；第2部分在扩散方程中引入扩散因子和锐化因子，j=αjs+βjd，其中α为锐化因子，β为扩散因子，用于控制界面的弥漫效果。如图8所示，生成的流体界面为弥漫界面，与LS方法相比，这种方法的优点是多相流混合过程中界面的厚度是可控的，扩散界面的细节比较突出。Li等人[69]在原有的CH方程中添加矫正项使得两相之间的界面满足光滑、均匀分布的特性，该方法简化了复杂的网格拓扑，并且捕获的界面效果要优于原始的CH相场模型。文献[70,71]提出一种用于相场模型的全自适应能量稳定的数值模拟方法，文献[70]通过引用自适应网格对扩散界面区域进行细化，可以使得扩散界面变得非常薄，同时引入改进的欧拉方法用于提高二阶空间有限差分方法的精度， 相比于基于常规网格的相场法，精度更高，但由于界面处网格是自适应细化的，因此当界面厚度较小时，需要较大的计算开销。

图8 相场法仿真效果[68]

CH相场模型虽然在描述流体扩散的弥散界面具有较好的表达优势，但由于CH模型涉及4阶空间导数，因此需要昂贵的计算资源，相比CH相场模型，AC相场模型的计算只涉及2阶的空间导数，因此实现较为容易，AC相场模型通常用于描述非保守场变量的物理量， 则2维情况下非守恒的AC相场模型可表示为

(12)

由于传统的AC相场模型不能保持质量守恒，因此许多研究者通过引入矫正因子[72-74]来解决AC的不守恒性，守恒的AC方程可写成如下公式：

(13)

其中λ(t)为引入的矫正因子，虽然文献[72-74]引入的矫正因子可以保证相场在演化过程中保持质量守恒，但是引入的矫正因子会使得模型过于依赖参数的设置，同时由时间依赖的对流-扩散方程存在高梯度或者不连续的缺点，会造成界面的局部震荡问题。针对这一问题，Bazilevs等人[75]利用相场方程的自身残差对原有的模型进行修改，同时将不连续捕获YZβ-DC算子添加到有限元公式中，使标量保持在预期的物理范围内，解决了传统相场模型对参数的依赖性。

求解多相流动方程除了在不守恒的AC相场方程中引入矫正因子外，还需通过细化网格来提高弥散界面的精度，如文献[76,77]通过对均匀网格进一步细化用于求解多相流的弥散界面，相比非细化网格的精度要高。

通过对比LS模型和CH相场模型可以发现，2个模型的定义很相似，其最主要的区别在于LS法在界面处的φ=0，而CH模型-1<φ<1，同时CH模型的控制方程比LS法多了一项扩散因子。在描述锐化界面特性上，LS法相比于CH模型和AC模型使用更广泛，LS法在计算曲率和法向上具有较好的优势，同时LS法描述的界面是零厚度光滑函数，界面尖锐性可以直接通过求解LS模型即可。

2.3 基于Lagrange粒子-Eular网格的混合方法

基于网格法的流体模拟由于不需要考虑粒子与邻居粒子之间的信息传递，因此计算效率相对较高。然而由于网格法是将问题域固定在网格上，因此网格的大小决定着流体细节质量。对于粗网格模拟的效果往往较为粗糙，细节表现不够，但可以达到实时模拟，而细网格模拟的细节效果更为逼真，但需要较大的计算开销。相反，SPH则通过将流体域建模成离散的粒子，更擅长流体的细节模拟，尤其在边界处的细节，但当粒子数增多时，计算较为耗时。由于网格法和粒子法都各有不同的优势和不足，因此，将这2种方法结合起来，可以利用两者的优势，使模拟的流体效果更加符合物理的运动。Lagrange粒子-Eular网格混合法的基本思想是在网格中分布一定的携带物理量的流体粒子，首先将粒子的速度插值到网格边中心进行计算，在进行平流时再将网格速度插值到粒子上，完成速度的更新操作。

常见的Lagrange粒子/Eular网格混合法主要包含PIC法和FLIP法。基于标记网格(marker-and-cell， MAC)[81]的插值思想，两者都是将粒子信息插值到网格上，这种插值方法的过程如图9所示。具体步骤：(1) 初始化粒子信息；(2) 将速度插值到网格上，并在网格上对压力进行求解；(3) 将速度插回粒子，并更新粒子速度； (4) 更新粒子的位置。

图9 粒子/网格插值过程

2.3.1 PIC法

PIC方法是将粒子与网格结合的一种数值方法，原始的PIC方法[82,83]将流体表示成携带质量和位置信息的粒子，然后直接在网格上计算速度，将下一时间步的网格速度插值回当前粒子的速度，在稳定性上具有较好的优势，但是由于信息在粒子与网格之间传输会存在丢失信息问题，因此产生高耗散的现象。

Edwards等人[84]为了提高PIC的数值精度，提出一种高阶精确的PIC方法，先是利用高阶移动最小二乘法近似网格上的粒子，然后通过有限差分法[85]来估算网格上的物理量，最后将这些物理量从网格插回粒子，并通过在粒子上添加正则化项来消除PIC的噪声，该方法具有较低的数值耗散。虽然文献[84]能达到较高数值精度，但还存在一定的数值耗散问题，为此，Jiang等人[86]提出一种仿射式的PIC法(affine particle-in-cell method,APIC)，将速度和位移的增量从网格插值到粒子，并过滤从网格到粒子的传输信息，该方法不但能消除原始PIC的人工耗散问题，而且还能克服FLIP的不稳定性问题。Hammerquist等人[87]通过拓展原有的PIC方法，提出一种扩展的PIC方法(extended PIC，XPIC)，该方法在速度中添加平滑项来降低FLIP的噪声，并通过最小化新速度与平滑速度之间的误差来解决PIC的数值耗散问题。

2.3.2 FLIP法

利用PIC方法能够获得较好的稳定性能，但存在较高的数值耗散问题，为此，Brackbill等人[88]提出一种用于低耗散的FLIP法，用于解决PIC数值耗散问题。Zhu等人[89]将FLIP方法引入计算机图形学，与PIC方法不同，FLIP是将上一时刻和当前时刻的网格速度的差值叠加到当前粒子的速度上，因此FLIP相比于PIC模拟的浪花效果更加明显，同时能避免PIC方法的数值耗散问题。

Boyd等人[90]拓展了文献[89]的FLIP方法，提出一种用于模拟两相流的MultiFLIP方法，首先在常规网格上求解N-S方程获得平流速度，然后利用粒子水平集方法建模不同相之间的界面，最后关联相位与MultiFLIP粒子，由于FLIP需要在流体域上保持粒子，因此用粒子重建的液体表面更为自然。Yang等人[91]通过自适应FLIP方法计算每个流体粒子的密度、压强更新粒子的速度和位置，并采用线性混合PIC和FLIP的速度，最后通过深度图直接在网格周围寻找弥散粒子，增强了复杂流体间相互作用的细节。Ferstl等人[92]基于文献[91]的重采样思想提出一种窄边带的FLIP方法(narrow band-FLIP，NB-FLIP)，通过对表面固定窄带区域不断重采样粒子，以控制表面窄带的宽度，该方法不但减少了模拟域的FLIP粒子，还解决了PIC法的耗散与FLIP法的不稳定性问题。邹长军等人[93]利用NB-FLIP方法模拟多相流现象，模拟单相流时与文献[86]得到的仿真效果相似，但是在处理两种流体界面时，在不同相交界处的界面明显分离现象。

综上分析，由于PIC方法在稳定性能上优于FLIP方法，因此较适合用于处理大形变问题，但会存在较高的数值耗散问题，其主要通过2种方式来降低数值耗散的影响:(1)引入正则化项；(2)过滤粒子与网格之间传递的信息。相比之下，由于FLIP在粒子与网格之间传递的是2个时间步速度的变化值，采用叠加的方式更新速度，而不是用下一时间步的速度覆盖当前粒子的速度，因此不需要考虑PIC的数值问题。

3 对比分析

基于物理的流体模拟方法的分析与比较结果见表1，对比指标主要包含算法的分类、算法名称、发表年份、技术特点、效率以及最大的仿真场景等。

从表1可以看出，根据算法的仿真场景和技术特点，基于粒子的方法从传统的SPH方法演化至IISPH方法，不可压缩性都比较好，然而在效率上，由于受到粒子数的限制，该类型算法在模拟小规模流体流动场景上具有较好的仿真效果，且大都以离线渲染为主，无法达到实时模拟。而与Lagrange粒子法相比，Eular网格法由于不需要对粒子进行计算和存储，可用于解决流体仿真的实时性问题，并且大都用于大规模海浪场景的构建。Level-Set、CLSVOF、VOSET这3种类型方法主要用于建模两相或多相流尖锐性界面的数值建模，该类算法模拟的多相流界面具有足够的光滑特性，且能够准确计算表面张力；而CH相场法和AC相场法主要用于建模多相流的弥散界面，由于该类方法是基于菲克扩散定理的，因此在流体流动过程中能够准确捕获到精细的界面细节。

表1 流体仿真算法对比

4 存在的问题及未来展望

综上分析，虽然近年来基于物理的不可压缩流体模拟仿真领域取得了突出的成果，但在实时性、复杂流体耦合等方面仍存在一些问题，存在的不足及未来可能的研究方向如下。

4.1 大规模场景的实时绘制

在目前的流体模拟算法中，大多算法都能实时模拟小规模流体，然而对大规模场景的流体模拟却很难达到实时，虽然使用高度场网格法模拟在实时性上具有一定的优势，但其模拟的大都为光滑流体表面，在细节效果上仍很难达到逼真的程度，同时对于航海类的大规模场景，往往需要处理的是更多的流体粒子或者更大的网格。因此，如何将应用于小规模流体的仿真算法扩展于大规模场景，并能实时处理流体与固体的耦合关系，对于计算机仿真具有重要的研究意义。

4.2 增强边界细节模拟

本文介绍的流体模拟方法中，基于粒子法的流体模拟为了表现更加精细的细节效果，往往需要增加更多的流体粒子，但这会造成较大时间和计算的开销，同时，当边界处的粒子较为稀疏时，会导致拉伸不稳定。而对于网格法，很难处理拓扑变形严重的边界，同样，如果通过细化网格来增加边界的细节效果，会损失一定的实时性能，因此如何保证流体模拟过程中的实时性与细节效果之间的平衡是未来的面临的一项挑战。

4.3 多相流界面控制

多相流模拟可以为虚拟场景增加更丰富的细节特征，然而现有的多相流仿真算法为了简化模拟过程，大多忽略温度、气流等外界因素对多相流界面演化的影响，而这些外界因素往往对多相流的效果具有很大影响，同时常规的多相流算法难以控制界面的运动以及厚度等，造成界面细节丢失的现象，因此在模拟多相流时包含一定的外界因素和控制多相流界面的生成会给多相流模拟仿真领域带来新的研究方向。

4.4 多种类型方法的耦合

目前针对多种方法融合的模型较少，大多为粒子与网格之间的融合，如文献[46]在高度场网格的基础上添加粒子，文献[51]提出的仿射FLIP方法，文献[52]的扩展FLIP方法，在效果上具有很大的提升。但这些方法并没有解决不同方法耦合衔接问题，尤其在2维网格与3维网格的耦合，如果处理不好可能会导致如文献[46]的人工伪影问题，因此未来仍需要提出一些新的耦合模型来解决现有耦合模型衔接问题。

5 结 论

如今，基于物理的流体模拟已经成为增强现实领域的主流方法，尤其是基于粒子法和网格法的仿真方法，为实时、逼真的流体动画场景提供了技术支持。近年来这些方法在特定领域问题上已经取得了很好的研究成果，如SPH方法在小规模流体仿真领域上的应用，高度场网格法解决大规模实时性问题等。本文在分析基于物理的不可压缩流体模拟方法存在的问题的基础上，对现有的流体模拟方法进行分类和综述，将有助于理解和研究流体仿真方法中的关键技术。