1. 引言
波涛汹涌的大海、倾盆而下的暴雨、瞬间就要吞噬战舰的巨大漩涡等流体影视特效给好莱坞著名奇幻电影《加勒比海盗3》的观众带来了震撼的视觉享受。如图1给出的电影中宏大海战画面完全用计算机模拟生成，由美国著名的影视特效公司工业光魔ILM (Industrial Light and Magic)制作完成，该公司从1977年至2012年总共获得奥斯卡最佳视觉奖16次。在欣赏这些精彩特效大片的同时，我们领略到了计算机流体模拟技术的高超魅力。
计算机流体模拟的实际应用不仅需要外观上的真实感，也需要其运动模拟的真实感。而描述流体的真实感运动只有借助于真实世界的物理规律才能得以实现，物理规律描述了流体如何运动，如何动作以及它们之间如何相互影响。随着具有实时处理能力的超级图形工作站的出现，上世纪八十年代后期发展起来一种新的计算机流体模拟技术——基于物理的流体模拟。在基于物理的流体模拟过程中，各幅离散图片的产生是由计算机对控制流体运动的高阶连续偏微分物理方程进行数值离散求解，自动生成每一时刻流体的运动状态。由于真实物理模型的支撑，基于物理的流体模拟技术能逼真地模拟各种复杂流体现象，给人以震撼的视觉享受，因此吸引了国内外大量图形学研究学者进行深入研究。Foster等 [1] 首次提出使用有限差分的欧拉网格方法求解流体运动控制方程Navier-Stokes方程来实现基于物理的流体模拟。为了实现大时间步长下稳定的流体模拟，Stam [2] 提出半拉格朗日方法迭代求解Navier-Stokes方程。自此，基于Navier-Stokes物理方程的流体模拟技术开始普及起来。在计算机图形学顶级国际会议SIGGRAPH上，自2001年以来每年都有流体模拟这个主题，发表的学术论文平均在4篇以上，其中绝大多数都是基于物理的流体模拟。经过十几年的发展，基于物理的流体模拟作为图形学领域中一种主流的流体模拟技术，已经被广泛地应用于影视特效、电子游戏、灾难营救、军事仿真、科普教育、体育竞技等领域。
基于物理的流体模拟方法主要分为欧拉网格方法和拉格朗日粒子方法。欧拉网格方法 [3] [4] 把流体占据的计算空间按一定规则剖分为网格，然后分析被流体所充满的空间中每一个网格位置上的流体的速度、压强、密度等参数随时间和空间的变化。该方法可以逼真地模拟多相流 [5] 、不连续流体 [6] 、粘弹性流体 [7] 、流体细节 [8] 等现象，但由于模拟过程中要同时把固体边界网格离散化 [9] [10] ，不易于进行固流交互模拟。而拉格朗日粒子法从分析流体各个微团的运动着手，即研究每个流体微团的速度、压强、密度等描述流体运动的参数随时间的变化状况，这是一种以“实体”形式为流体建模的方法，建模的基本元素就是携带流场变量的粒子。该方法直观，易于进行固体边界处理和模拟流体细节。论文以拉格朗日粒子流体中的SPH方法为理论基础对固流交互模拟技术进行研究，下面结合研究内容对基于SPH的流体模拟技术的国内外研究现状进行综述。
2. SPH方法原理
SPH方法 [11] [12] 是一种用于求解连续介质动力学方程的无网格方法，它通过引入空间场函数和核函数的概念，将流体控制方程离散化。该方法避免了网格方法中存在的网格缠绕和扭曲问题，在计算机物理动画中主要用来进行流体的运动模拟。
SPH方法的基本思想是用粒子来离散化流体的物质空间，每个粒子代表了一定的流体质量，其携带的各种数值量将光滑地分布于以该粒子中心为球心，半径为h的区域内，该区域称为粒子的紧支域。如图2所示。
由图2可知，SPH方法的核心是核函数
，它是一种在一定光滑长度h范围内其他邻近粒子对目标粒子影响程度的权重函数。在SPH中任一宏观变量
(如速度、压力等)及其梯度都能通过一组无序质点的核函数插值集合来表示，进而可以将流体力学方程表示成SPH和插值形式，如公式(1)所示。
(1)
其中，
是代表粒子
的邻居粒子索引，
表示质量，
表示位置，
表示光滑半径，
表示粒子
的
Figure 1. Fluid effects in movie Pirates of the Caribbean 3
图1. 来自电影《加勒比海盗3》的流体场景特效
密度。公式(1)对应的梯度及拉普拉斯算子如公式(2)、(3)所示。
(2)
(3)
其中，
和
分别代表梯度和拉普拉斯算子。
3. SPH模拟国内外研究现状
3.1. 不可压缩SPH流体模拟
Müller等 [13] 首先利用SPH在计算机图形学领域实现了交互级流体自由表面模拟，该方法考虑了流体表面张力的作用，并通过理想气体状态方程来保证流体的微可压缩性。实现SPH流体的不可压缩性是SPH流体模拟技术的关键内容和技术难点。对SPH流体方法来说，模拟其不可压缩性主要有两种不同的策略：一种是微可压缩SPH流体 [14] (Weakly Compressible SPH, WCSPH)，该方法通过Tait方程由密度计算压强，因此为了提高流体的不可压缩性，必须在Tait方程中使用大的刚度系数，而这反过来要求必须使用较小的时间步长，严重影响了计算效率；另一种是通过计算压强的泊松方程来确保不可压缩SPH流体(Incompressible SPH, ISPH)模拟 [15] [16] ，该方法可以使用较大时间步长，但在每个时间步长内由于迭代求解线性方程组而导致计算量庞大。为提高泊松方程的计算效率，Ihmsen等 [17] 提出了一种新的投影方法，该方法集合了对称的SPH压力和连续方程离散方法来对离散化压强泊松方程，并对压力进行实际计算，提高了求解的收敛速度(图3(a))。针对WCSPH方法的步长限制问题，Solenthaler等 [18] 采用一种预测修正方案计算流体粒子压强，实现了较大时间步长下的不可压缩SPH流体模拟。该方法结合了WCSPH和ISPH的优点，允许较大时间步长，并简化了时间步内的计算复杂度，其核心思路是通过对密度进行预测修正的方式迭代计算流体压强和压力(图3(b))。
在对SPH流体施加不可压缩约束的同时会降低计算效率。自适应采样 [19] 和双层采样 [20] 通过对视觉重要的局部区域进行细分粒子，在远离视点的地方合并粒子的方式提高模拟效率，但复杂的剖分和合并操作也降低了部分性能，并且过小的粒子限制了时间步长的大小。由于具有粒子之间良好的数据独立特性，GPU并行计算技术 [21] - [23] 也被成功引入到SPH流体模拟中，实现了一定粒子规模下的实时模拟。
3.2. 固流交互模拟
SPH流体的固体边界处理方法一般分为三种：Penalty Force方法、Ghost Particle方法和Frozen Particle方法。Penalty Force方法通过对流体粒子施加一个与固体边界的距离成反比例关系的反向作用力来实现固流交互。2004年，Müller等 [24] 利用Penalty Force方法模拟了微可压SPH流体和有限元变形固体的交互，该法对固体的三角形表面采样生成边界粒子，并计算SPH粒子到这些边界粒子的距离。然而，在这种方法中，必须对计算惩罚力的刚度参数进行精心调整以防止固流交界面处的压强分布噪声和穿透现象。因此惩罚力方法通常必须使用比较小的时间步长来保证固体边界处流体压强的均匀分布。为了采用较大的时间步长，Becker等 [25] 提出了Direct Forcing方法，该方法利用一个预测修正方案来计算交互力和速度，实现了流体和刚体的单向和双向交互。直接作用力方法避免了固体交界面处的流体粒子堆积现象，但是不能处理可变形固体边界。2012年，Yang等 [26] 提出了一个基于GPU的实时耦合方法来处理SPH流体和非线性有限元变形模型的交互，该方法结合了直接作用力方法和预测修正方案，可以处理不同的
固体边界并避免穿透，但会受到时间步长和压强分布噪声的限制。
Ghost Particle方法通过为固体边界处的流体粒子动态生成镜向虚拟粒子，并把这些粒子引入到流体粒子的密度计算中来实现稳定的固流耦合，该方法可以保证固流交界面处流体压强的均匀分布。Hu等 [27] 和Morris等 [28] 使用镜向粒子方法分别实现了SPH流体和直线和弯曲固体边界的交互。2012年，Schechter等 [29] 提出一种粒子采样方法为流体的自由表面和固体边界生成一个狭窄的镜向粒子层。如图4(a)所示，该方法可以解决假性数值表面张力效果，并可以保证质量守恒，但是对于拓扑复杂的可变形边界很难生成镜向粒子。
2007年，Solenthaler等 [30] 提出了一个统一粒子框架来模拟流体、刚体、和可变形物体以及它们之间的耦合。该方法通过对固体边界进行粒子采样，并考虑这些粒子对流体粒子密度以及压力的贡献，因此保证了固流交界面处流体密度的均匀分布，但必须采用小的时间步长来保证固体变形的稳定性以及固体表面的非穿透。Ihmsen等 [31] 通过把上述方法与直接作用力方法结合实现了一个新的单向固流耦合方案，该方法不仅可以保证固体边界处流体密度的光滑分布，还可以支持较大时间步长。Akinci等 [32] 通过在固体表面采样一层边界粒子，并考虑这些边界粒子对流体的相对贡献实现了SPH流体和刚体的动量守恒保持的双向交互。该方法解决了由于粒子采样导致的边界非均匀问题，但是不能解决可变形固体边界(图4(b))。Akinci等 [33] 通过自适应采样边界粒子方案把该方法扩展到了可变形固体边界的处理(图4(c))。然而在较大时间步长和固流速度差情况下，仅仅依靠一层边界粒子的自适应采样仍不能保证边界的非穿透。2015年，Shao等 [34] 为解决SPH流体在固流交界面处的粒子不连续，以及较大时间步长或速度差下的穿透问题，提出一种基于统一粒子模型的稳定快速固流交互方法。该方法对固体表面和内部进行边界粒子采样，结合一种动量守恒保持的速度–位置修正方案，实现稳定的不可穿透固流交互模拟，并给出GPU并行计算方法，达到了120k粒子规模下10帧/秒的交互级模拟速度。
通过对固流交互模拟的综述和分析可知，固流交互界面的不稳定性包括压强噪声以及穿透仍然是基于SPH的流体模拟亟待解决的问题。
(a) ISPH不可压缩流体模拟 [17] (b) PCISPH不可压缩流体模拟 [18]
Figure 3. Incompressible SPH fluid simulation
图3. 不可压缩SPH流体模拟效果
(a) Ghost Particle方法 [29] (b) Frozen Particle方法 [32] (c) 自适应边界粒子采样方法 [33]
Figure 4. Fluid-solid simulation of SPH fluids
图4. SPH流体的固流交互模拟
3.3. 流体的漩涡细节模拟
由于数值方法固有的数值耗散，漩涡细节的丢失对于基于物理的流体模拟是不可避免的。真实有效地恢复流体模拟的漩涡细节一直是基于物理流体模拟的研究内容和技术难点。对于欧拉网格流体模拟方法，Fedkiw等 [35] 提出了一种Vorticity Confinement方法来恢复本来存在的漩涡细节，但是如果一致网格的分辨率不足以捕获所需精度的细节特征时，该方法不能恢复该精度以及更细小的漩涡细节。为了降低欧拉网格方法的数值耗散，一些研究者采用了更高阶的对流项 [36] - [38] ，但是这些方法仍然受网格分辨率的限制，并且由于拉哥朗日粒子方法中不需要计算对流项，因此该方法不适合粒子流体。自适应网格剖分方法 [4] [39] 通过对感兴趣的局部区域进行细分以增加计算节点数量来恢复细节特征，此类方法也应用在了拉格朗日粒子流体方法SPH中，如自适应SPH [19] 和双分辨率SPH [20] 。此类方法的缺点是只能恢复和保持局部区域内的漩涡细节，当这些细节脱离局部区域后就变得不能维持。
漩涡合成方法通过对求解流体控制方程得到的平均速度场进行合理扰动来恢复流体的漩涡细节特征，该类方法同时适用于欧拉网格方法和拉格朗日粒子方法。Stam [40] 利用一个Kolmogorov spectrum来产生欧拉网格流体中漩涡细节，而Bridson等 [41] 提出利用向量噪声场卷积来产生漩涡速度场。Kim等 [42] 利用小波分解来检测流体中丢失细节的位置，然后利用一个不可压湍流函数来重新生成这些细节。Schechter等 [43] 利用一个简单的线性模型来跟踪流体中的湍能分布带，然后使用噪声函数创建湍流速度场。Narain等 [44] 提出了一个快速有效的流体漩涡模拟方法，该方法通过一个可以集成到欧拉网格流体中的子网格湍流模型来自动地跟踪漩涡的产生和演化。
目前，在流体中引入带有涡度信息的涡度粒子(Vortex Particle)成为一类主流的漩涡合成方法。在图形学中的流体动画领域，Selle等 [45] 在欧拉网格流体中引入涡度粒子来增加流体的漩涡细节，但是该方法需要专业美工人员来确定涡度粒子的插入位置以达到预期的模拟效果。Park等 [46] 通过在整个计算空间引入涡度粒子并利用纯粹的拉格朗日方法求解涡度传递方程模拟了欧拉网格烟雾的漩涡，而文献 [47] - [49] 中方法在欧拉网格流体中引入涡度粒子来模拟具有高真实感的漩涡细节。固流交互过程中产生的漩涡模拟如图5所示。
由图5可以知，流体中产生漩涡的其中一个重要原因是由于固体对它的扰动，也就是固流交互会产生漩涡细节。Praff等 [50] 在欧拉网格流体中通过近似边界层理论在固体周围预计算一个涡度场，在运行阶段计算涡度场在固体表面的脱落点并在脱落点处插入涡度粒子来模拟固流交互产生的漩涡(图5(a))，但该方法中存在的预计算涡度场步骤比较耗时。上海交大的朱博等 [51] 模拟了SPH烟雾和刚体以及弹性体交互产生的漩涡细节(图5(b))，该方法通过对物体绑定一个高分辨率局部网格来生成漩涡涡度场，这个局部网格随物体进行刚性的旋转和平移运动。Jiang等 [52] 把艾米插值算法(Hermite-Interpolation)引入到粒子的对流处理过程，并提出了大尺度核函数和小尺度涡度模型来有效地捕获流体中的多尺度漩涡细节，该方法同样可以成功模拟固流交互中产生的漩涡。通过在物体周围插入SIP (Swirling Incentive Particle)粒子，Yuan等 [53] 提出了随机漩涡可能性模型来模拟固流交互中的随机漩涡细节(图5(c))，但该方法会在固体的前面产生不合理的漩涡细节，影响模拟的真实感。针对这一问题，Shao等 [54] 在物理模型上，给出一种基于边界层理论模型的漩涡合成方法以恢复SPH流体中固流交互产生的涡度场，并对其进行稳定演化；在几何表示上，给出一种高效的基于八叉树的自适应流体表面重建方法以进行漩涡的真实感绘制。
综上所述，真实地模拟SPH粒子流体中固体边界引起的漩涡，对于增强SPH流体模拟的真实感具有重要意义。
3.4. 气泡现象模拟
气泡作为液体的一种细节特征是真实世界中普遍存在的现象，捕获气泡细节有助于增强基于物理的
(a) 预计算边界层方法 [50] (b) SPH-Flip方法 [51] (c) 随机漩涡模型 [53]
Figure 5. Turbulence in fluid-solid coupling
图5. 固流交互中的漩涡
流体模拟的真实感。在纯粹的欧拉流体模拟器中，文献 [55] [56] 中算法利用Regional Level Set方法模拟相对于水体来说较大的气泡，但对于尺度小于网格分辨率的气泡以及泡沫，纯粹的欧拉方法不能够进行有效模拟。
通常采用混合模型来模拟水中气泡：把粒子表示的气泡耦合到各种类型的流体模拟器中。Kuck等 [57] 用一定大小的球形粒子表示气泡以避免处理气泡的拓扑变化，并着重处理气泡之间的交互和融合。Hong等 [58] 结合流体体积方法(Volume of Fluid)和前向跟踪方法(Front-tracking Method)模拟水中气泡的上升，并利用球形粒子模拟不变形的小气泡。Greenwood等 [59] 利用粒子水平集算法(Particle Level Set Method) 的标记粒子模拟气泡，但没有考虑气泡形状变化。Kim等 [60] 把粒子系统和低分辨率的欧拉网格结合从而有效地模拟了水中的大量气泡。上述气泡模拟方法通过利用一个个孤立的小球表示每个气泡，忽略了气泡之间的交互，模拟的气泡大小和形状始终保持不变。
由于气泡本身也是一种流体，因此很多流体模拟器中利用SPH方法模拟水中的气泡运动。为了模拟欧拉网格流体中小于网格分辨率的泡沫，Foster等 [61] 在粒子水平集方法中引入SPH粒子来模拟溢出空气产生的气泡。而Hong等 [62] 通过一个速度场把SPH方法与欧拉网格方法耦合以模拟水中的气泡，然而由于有限的网格分辨率，该方法不能模拟气泡上升过程中的路径不稳定性。为了高效地模拟水中的气泡，Thurey等 [63] 把表示气泡的粒子系统和表示水体的二维浅水模型(Shallow Water Model)结合，并作出了大量简化假定来产生三维的效果，只模拟气泡周围的水体和在水中的气泡。该方法计算高效，但存在很多不足，例如不能模拟流体的惯性。
统一的粒子框架广泛用于模拟水以及其中的气泡细节。Müller等 [64] 利用SPH方法模拟了最大密度比为10:1的多相流，并引入浮力模型模拟水中气泡，但其密度模型会在流体和气体交界面处产生错误的密度估计，进而导致错误的流体压强计算。针对这一问题，Solenthaler等 [65] 提出的非连续密度模型在计算流体粒子当前密度时忽略粒子的质量，可以在不同流体交界面处产生大的密度变化。该方法可以支持最大密度比为100:1的多相流模拟，但仍不能很好地模拟水和气泡(密度比接近1000:1)的交互。为了能更支持更大密度比的多相流模拟，Cleary等 [66] 提出了一种混合粒子模型来模拟由于溶解气体产生的气泡，例如啤酒、可乐。该方法独立地模拟气体和水：利用固定形状的离散实体模拟气泡，而用SPH方法模拟液体，并且液体对气泡的作用通过一个拉力模型来计算。类似地，Markus等 [67] 利用单相SPH方法模拟了气泡和水的运动，然后通过一个拉力模型计算两者之间的相互作用。固流交互中会产生气泡，如图6所示。
由图6可知，气泡的模拟可以提高固流交互模拟的真实感。Mihalef等 [68] 提出了一个新的欧拉网格方法来模拟沸腾水中的气泡，并计算温度驱动的气泡和水之间的质量转移。在模拟过程中，该方法在预先随机选定的气泡成核点处引入气体粒子，并考虑气体和固体之间的相互作用。Hong等 [6] 对欧拉网格方
Figure 6. Simulation effects of bubbles in solid-fluid interactions
图6. 固流交互中气泡的模拟效果
法进行了扩展以模拟不可压缩粘性多相流体，并重点模拟了多相流中的表面张力作用、交界面处的粘度变化以及浮力效果，但该方法只考虑了气泡和静态固体的相互作用(图6(a))。Mihalef等 [69] 提出了一个拉格朗日粒子模型和欧拉网格模型相结合的混合模型来模拟固流交互中由于带入的气体产生的气泡，该方法没有考虑气泡之间的粘滞力作用(图6(b))。Shao等 [70] 提出一种基于统一粒子多相流模型的气泡模拟方法，该方法在SPH流体中对气泡成核理论模型进行近似，给出一个综合考虑水中气体溶解度、固体材质以及固流速度差等物理因素的气泡生成模型，并定义拉力、表面张力、浮力以及凝聚力的计算模型来计算气体粒子受到的作用力，真实地模拟了气泡的变形、上升、融合和泡沫等现象。
综上所述，在有限的计算资源下，建立合理的近似模型对固流交互过程中产生的气泡进行模拟是SPH流体模拟的重要研究内容。
4. 结论
通过对其国内外研究现状的回顾与分析，我们可以得到，基于物理的流体模拟技术，特别是SPH流体模拟技术已经能够真实地模拟很多复杂的流体现象。但受非线性偏微分方程数值求解方法的收敛条件约束、固有耗散性以及计算求解复杂等因素影响，基于SPH的固流交互模拟仍存在界面不稳定、流体细节丢失以及计算效率较低等问题。
因此，SPH流体模拟的后续工作包括：
(1) 基于SPH的复杂流体现象的模型建立。真实世界的复杂流体现象是多种物质，多种物理规律共同作用的结果，如固流交互同时存在着液滴飞溅、漩涡、泡沫以及固体的浸润和侵蚀，真实地模拟这些现象有待在图形学角度提出近似物理模型。
(2) 实时SPH流体模拟。由于物理方程数值求解计算复杂，在大规模粒子下不能实现实时或交互性模拟是制约SPH流体模拟技术实现广泛应用的重要原因，因此设计高效的并行求解算法，构建多CPU和多GPU耦合的集群对复杂物理求解过程进行大幅度加速是当务之急。
(3) SPH流体的粘弹性固体边界模拟。固体弹性变形的一个重要应用是虚拟医学中的软组织模拟，而真实的软组织并不只是弹性的，还存在粘性以及非均匀一致性。因此为了保证虚拟手术的精确性，建立符合人体生物学特性的软组织变形模型至关重要。