1. 引言
光滑粒子动力学SPH (Smoothed Particle Hydrodynamics)是一种纯Lagrange插值方法，最早用于解决天体物理学问题，和其他方法相比，此方法完全不需要网格，避免了网格扭曲和网格重构等问题，可以用于任何流体的模拟。SPH方法在Monaghan [3] 和Münzel [4] 的文献中有介绍。
在真实的流体模拟中，泡沫和喷雾效果是很重要的特征 [5] - [7] ，比如电影 [8] - [10] 和商用软件 [11] 。流体飞溅过程中会产生喷雾、气泡、泡沫等细节现象，而目前的方法忽略了这些细节现象的建模，致使动画效果的视觉逼真性缺失。所以，本文在传统SPH方法的基础上，采用了一种新的固壁边界处理方法和一种新的模型 [2] 。模型的主要作用是增加粒子扩散后的处理过程。
2. 基本方程
2.1. 流体控制方程
在二维Lagrange坐标系下，Naiver-Stokes质量守恒方程为：
(1)
动量守恒方程：
(2)
式中
为密度，
为重力加速度，
为速度矢量，
为压强，
为粘性系数。
2.2. 状态方程
在SPH方法中，常将不可压缩流体通过合适的压力状态方程看作微可压流体。本文采用文献 [12] 中的状态方程：
(3)
其中，
,
为声速，
。关于状态方程的讨论可参考文献 [13] 。
3. SPH离散
控制方程的SPH离散
在SPH方法中，对N-S方程有多种离散格式。对溃坝问题，边界处会出现大的伪振荡，所以本文采用文献 [14] 中比较光滑的离散形式：
(4)
(5)
式中
，
，
分别表示粒子
的速度、位置、压强、密度和粘性系数，
为核函数，
为光滑核函数距离。本文使用五次样条内核插值 [15] 。
文 [16] 介绍另一种内核函数，感兴趣可参考。
4. 固壁边界条件
为了保证固壁边界满足不可穿透和无滑移条件，本文采用一种新的处理方法。在固壁上，布置一层固壁粒子。对于固壁粒子
，其压力计算式为
(7)
式中
为
邻居流体粒子。
在固壁外再布置一层固定虚粒子，保持位置不变，压力和速度按照下列公式计算，见文献 [1]
(8)
表示向量
，为了精确求解
，再次用到(7)式。
可以看出，这种固壁处理方法利用周围流体粒子来计算固壁边界的压力，避免固壁边界的粒子出现伪震荡，并且不需要施加任何额外的排斥力就能够有效的防止粒子穿透，提高了边界处的精度。
5. 流体粒子扩散
流体粒子扩散通常发生在自由表面上，能够产生扩散现象的有三种情况 [17] 。一是空气滞留。二是在波峰处，流体粒子会飞溅，产生扩散粒子。三是流体颗粒自身有一定的能量，能量越大，产生扩散粒子的可能性越大。
5.1. 滞留气体
空气受到冲击，会滞留在浅水地区，高度紊乱的流体会把空气吸入水中，导致空气很有可能滞留在水中，形成气泡或泡沫粒子。
用控制函数来表示一个流体颗粒产生扩散粒子的概率：
(10)
流体粒子的相对运动，可以用式子
表示，其中
表示两个粒子单位相对速度，
表示单位距离向量。滞留空气的概率
，其中速度差分式子：
(11)
和
表示速度差分的最小最大值，由用户输入。
是对称权函数，定义为：
(12)
是粒子搜索的邻居半径。
5.2. 波峰
波峰不稳定，容易产生扩散粒子。在波峰处，表面曲率比较高，局部表面呈凸形。表面曲率
近似为
(13)
是单位表面法向。为区分凸凹区域，考虑
和
的向量积(如图1)。
所以波峰可以表示：
(14)
其中，
(15)
文 [17] 提出额外的条件对凸面进一步修正，
(16)
在波峰处，一个流体颗粒产生扩散粒子概率
，
,
由用户输入。
5.3. 能量
雾沫和液滴的计算需要精确知道表面张力的改变，这对于模型是一大挑战。本文未考虑表面张力的改变，文献 [18] 和 [19] 把表面张力设为常量。仅把内能
作为产生雾沫、液滴的衡量标准。所以内能产生扩散粒子的概率为
。
表示内能的最小最大值，由用户输入。
所以，一个流体粒子产生扩散粒子数为：
(17)
其中
表示每秒取样的波峰数，
每秒取样的滞留空气数，
为时间步长。
5.4. 扩散粒子的产生
图2.一个流体粒子产生许多扩散粒子，现在流体粒子位置是
5.5. 扩散粒子分类
根据产生出的扩散粒子位置，计算其周围邻居流体粒子的数量，把扩散粒子进行分类。扩散粒子的邻居流体粒子数小于4，记为喷雾粒子；粒子数大于12，记为气泡粒子；其它为泡沫粒子。
Figure 1. The normal direction of the particle does a point multiplication algorithm with its neighbour distance vector. The positive is concave surface, negative convex
图1. 粒子法向与粒子相对邻居粒子距离向量点积，为正则是凹面，为负是凸面
Figure 2. A fluid particle produces many diffusion particles. Now the fluid particle position is xf(t)
图2. 一个流体粒子产生许多扩散粒子，现在流体粒子位置是xf(t)
喷雾粒子仅受到重力和外力的作用，使用欧拉方法更新速度和位置，
，
。但是，泡沫和气泡粒子很大程度受到周围流体粒子的影响，泡沫位置的速度是通过局部流体粒子速度平均值体现的。所以扩散粒子的位置更新为
，速度不更新。其中
，
表示流体粒子，
表示扩散粒子。由于水和空气的密度差较大，气泡的运动还受到浮力的作用，所以气泡速度更新为
，
和
控制着浮力和重力的影响，由用户输入。位置更新为
。
6. 数值算例
本文以溃坝现象为例，说明此模型模拟流体细节的视觉效果，各参数见文献 [20] 。
6.1. 水柱倒塌模拟
图3给出了SPH方法得到的波前位置，并与文献 [21] [22] 得到的实验测量值进行比较。本文使用的时间步长是固定，容易计算。自适应步长运算比较繁琐，见 [23] 。可以看出，使用此种边界条件得到的波前位置更接近实验值。因此本文所采用的边界处理方法能够准确和稳定地模拟整个水柱倒塌过程。
图4给出了不同时刻的流场，并与实验 [20] 得到的实验现象(详见 [20] )进行比较，可以看出得到的流场和文献 [20] 结果一致。把图4(d)放大，可以从图4(g)中看出，增加流体粒子的扩散处理，使水柱倒塌的视觉效果更加逼真。
6.2. 有挡板的水柱倒塌模拟
为了研究SPH方法在实际问题中的应用，图5给出了有挡板的水柱倒塌过程。
与无挡板模型相比，此模型在水平中点上方，放置了一个挡板。当水柱流过挡板，挡板下方的水流类似实际问题中的泄洪过程；而在挡板上方，当水柱冲击挡板，类似实际问题中的流体碰撞固体，产生飞溅等现象。喷雾和泡沫粒子一般产生在流体表面，气泡粒子产生在流体粒子和墙壁相撞的角落处。数值结果表明，此模型在有挡板的水柱倒塌场景下，比无挡板的水柱倒塌产生的飞溅现象更明显，视觉效果更逼真。
7. 结论
本文通过使用了流体粒子扩散新模型，增加了流体粒子扩散的后处理过程，模拟了无挡板和有挡板
Figure 3. The motion of leading edge at different time
图3. 不同时刻对应的波前位置
(a) t = 0.1 s(b) t = 0.2 s 
(c) t = 0.35 s(d) t = 0.4 s 
(e) t = 0.45 s(f) t = 0.5 s
(g) t = 0.4 s
Figure 4. The collapsing simulation of the water column
图4. 水柱坍塌模拟
的溃坝问题。数值结果表明，此固壁边界能稳定准确的模拟溃坝问题，此模型能使溃坝的视觉效果更加逼真。新产生的喷雾、泡沫、气泡扩散粒子很大程度提高了流体在实际流动过程中的视觉效果，使流体动画更加逼真。