光滑粒子水动力学在水利与海洋工程中的应用研究进展

林鹏智，刘 鑫

(1.四川大学水力学与山区河流开发保护国家重点实验室，四川成都 610065;2.河海大学港口海岸与近海工程学院，江苏南京 210098)

1 光滑粒子水动力学简介

近年来，随着计算机硬件技术的不断进步，各种新的数值模拟方法逐渐被科研和工程领域重视并得到了广泛应用，其中光滑粒子水动力学(smoothed particle hydrodynamics，SPH)是最具代表性的一种。不同于传统的基于欧拉网格的数值计算方法，SPH基于拉格朗日系统建立控制方程，在整个模拟过程可以完全不依赖计算网格［1］，这使得SPH具有一些网格法所不具有的优势，如适于模拟具有复杂或者大变形自由面的水流问题和处理具有复杂流固交界面甚至是运动结构物的流固耦合问题。自SPH在1977年被两位学者［2-3］独立提出开始，其应用范围已经从最初的天体物理学发展到目前的流体动力学问题。

1.1 拉格朗日描述的控制方程

在计算流体运动时，SPH通过求解 Navier-Stokes方程(N-S方程)或它的其他演化形式(如水深平均的浅水方程等)解决问题。拉格朗日描述下的一般形式的N-S方程如下:

式中:u为速度;p为压强;g为重力加速度;ν0为分子运动黏性系数;ρ为密度。相比传统网格法中采用欧拉描述的控制方程，方程(2)中不存在对流项，网格之间的对流效应信息由拉格朗日粒子的自身速度携带。

在对方程(1)(2)进行数值求解时，SPH数值模型可细分为弱可压缩模型(weakly compressible SPH，WCSPH)和不可压缩模型(incompressible SPH，ISPH)两大类。WCSPH模型中，方程(1)包含流体密度变化项，辅以虚拟的流体状态方程直接将密度变化和压强建立联系，可以在不求解矩阵的情况下显式获得压强。但为了使时间步长不至于太小，往往会采用较小的人工弹性模量，在计算高速流动时可能产生较大的误差。ISPH模型则采用不可压假设，将方程(1)简化为速度散度为零的形式，即

然后再通过求解压力Poisson方程(pressure Poisson equation，PPE)获得压强，进而获得速度场。

一般讲，ISPH模型能够采用更大的时间步长进行计算，得到比WCSPH模型更加光滑稳定的压力场，在模拟剧烈的自由面变化等问题时压力场更加稳定，但对边界条件的处理要求也更加严格［4］。WCSPH模型则更适合并行计算的要求［5］，为了获取更加稳定的压力场，一般可以通过引入人工黏性等数值处理方法。当前学术界对二者的研究成果都很丰富。

1.2 计算域的SPH离散

根据狄拉克函数的数学性质，计算域中任意点x0处的变量(如密度、速度、压强)可以精确表示成其周围点的积分形式:

式(4)被称为场变量的积分表示，而δ即δ函数，或称狄拉克函数，Ω为积分域。在SPH模拟中，使用离散的拉格朗日点离散上述连续积分，并用多项式函数W代替超越函数δ，将上述积分写成离散形式，得到积分表示的粒子离散形式(图1):

式中:i为当前参考点的变量值;j为积分域内离散的拉格朗日点的变量;W为核函数，也有学者称其为权函数或者光滑函数，是用以代替狄拉克函数的多项式函数;h为光滑长度，κh为搜索域半径，系数κ根据模型不同而不同，通常取2～3;V为积分自变量，在二维情况下是圆的面积，而在三维情况下则是球的体积。关于粒子离散的详细介绍可见文献［6］，此处不再赘述。具体到水动力学问题中，拉格朗日点成为具有一定质量、体积、速度、压强等属性，并携带位置信息的粒子，而j则表示搜索域Ω内i的所有相邻粒子。这时，根据式(5)，粒子i的SPH表达式为

图1 SPH方法中的粒子离散

式中m为粒子质量。

总体来说，与一般的网格法相比，SPH方法有以下优势和特性:

a．由于SPH方法是一种纯拉格朗日粒子法，使得其可以方便地描述水流的自由表面，而不需要传统欧拉方法中对自由面的捕捉和重建操作，这也是SPH非常适合处理具有复杂自由面问题的原因。

b．由于SPH的无网格特性，使其可以灵活地处理流固耦合问题而不需要采用传统欧拉方法中复杂的网格重建或者网格切割操作，大大简化了流固耦合求解流程，使其非常适合模拟流固耦合问题。

c．基于拉格朗日法特性，SPH避免了对非线性对流项的离散操作，简化了求解过程并且避免了相应的数值误差。

d．不同于欧拉网格法中的处理，单相SPH模型不需要对计算域内非流体区域划分网格，减少了对计算资源的额外消耗，对于非连续流体(如飞溅水流)问题的模拟更加高效。

本文将围绕SPH的上述特点和优势，着眼于水利与海洋工程中各类实际问题，在对其适用的问题进行分类的基础上，从工程应用的角度论述SPH在水动力学领域的研究进展，并对SPH未来的发展方向进行展望。

2 SPH在纯水流问题中的应用

SPH被大规模应用在水利和海洋工程研究领域是从Monaghan开始的，他在文献［7］中详细讨论了SPH的自由面边界处理方法、固体边界处理方法以及N-S方程的求解流程等问题，为以后SPH在水利工程中的应用研究奠定了基础。SPH在水动力学问题中最广泛的应用领域是纯水流问题，即河流、海洋中的各种水体流动问题。

2.1 波浪模拟

单纯的波浪问题是海洋工程中最基本的水动力学问题，也是迄今为止SPH研究的热点问题之一。波浪属于典型的自由面大变形问题，传统欧拉网格法中对其模拟需要依赖高质量的界面追踪算法如VOF(volume of fluid)模型［8］、Level set方法等，这些方法均要求对自由面附近网格进行特殊处理，重建自由面位置并保证自由面位置清晰，而SPH基于无网格离散，流体自由面被拉格朗日粒子自动描述，省去了自由面重建操作，简化了模拟［9-11］。SPH波浪问题的应用研究主要集中在粒子法中造浪的方法及波浪与海岸结构物相互作用的模拟方面。

2.1.1 造波方法

SPH作为一种粒子法，其造波技术也有其相应的特点和要求。Monaghan［7］最早进行了相关工作，他根据孤立波的理论波形设置了初始粒子分布，并给定粒子群孤立波的理论速度场从而简单地产生波浪。很显然，该方法很难输出需要反复输入动量的周期波。Monaghan在文献［7］中也提出了推波板造波法，即使用一个类似物理固体边界的粒子墙，按照实验造波理论给定其运动规律，进而产生波浪。该方法应用非常广泛，且涉及的波浪形式也多种多样，如Gotoh等［12］在其MPS模型中(一种和 SPH非常类似的模型)使用了造波板产生椭圆余弦波(cnoidal wave)并研究了波浪浅化破碎的特性。在后续研究中，周期波［13］、孤立波［14］甚至是随机波［15］都被成功模拟。另外，还曾有学者对运动结构物运动产生的波浪［16］以及滑坡体涌浪［17］进行了一些研究，这些内容已经涉及了一些流固耦合的问题。

如前面所述，推板造波取得了巨大的成功，但也存在着无法克服的缺陷，即存在二次反射问题，当计算域中存在结构物时，结构物的反射波重新与固体造波板作用形成二次反射而扰乱入射波参数。一个比较简单的回避方式是增加计算域延后二次反射到达的时间［18］，但仅能用于短时间波浪模拟情况。Frigaard等［19］提出了一种主动吸波造波板理论，造波板运动的一个额外的分量恰好吸收掉反射波而避免了二次反射，但这种方法非常依赖对反射波的监测。Hayashi等［20］在 MPS模型中使用主动吸波造波板实现了长时间波浪破碎和越顶模拟;而Didier等［21］则在SPH中整合了无反射和黏性吸波层，研究了波浪的相互作用问题。

另一种方法被称为“内造波法”(internal wave maker)，它通过源项代替造波板产生波浪。源项只作用于源项区域(source region)内部的流体粒子，通过周期变化源项驱动源项区流体发生运动，产生水面波动，形成向外传播的行进波。结构物的反射波会穿过源项区域然后被消波层耗散，避免了二次反射。Larsen等［22］是最早使用该方法进行造波的人，他们基于Boussinesq方程实现了质量源项造波。之后的学者把质量源项拓展到了缓坡方程［23］和N-S方程［24］，又在其基础上发展了动量源项造波［25-27］。而由于SPH方法的拉格朗日特性，使得其难于进行质量源项造波法中的质量吞吐操作，所以动量源项造波法相比之下更加适合。基于动量源，Liu等［28］改进了内造波算法并将其引入到ISPH模型中，系统地测试了SPH内造波算法对各种源项区域尺寸、波浪参数的敏感性，并将其应用到典型的波浪结构物相互作用问题中，展示了ISPH内造波方法的应用潜力。

2.1.2 波浪破碎模拟

波浪破碎、水体抨击等属于典型的自由面剧烈变形、水体严重不连续的物理现象，是传统网格法比较难于处理的方向，也是SPH的优势所在。由于拉格朗日特性，SPH不需要进行FDM、FVM中的追踪自由面操作，算法简洁有效。对于破碎波的模拟，技术问题主要包含两方面内容，一是自由面处的边界处理，二是有效的紊流模型。

SPH不进行自由面重建操作，但对自由面处的积分域截断问题仍需要进行额外的处理，最传统的简单方法是根据自由面处粒子数密度(number density)的变化判断自由面粒子身份，然后对自由面粒子通过镜像附近粒子的方式补全积分域，施加边界条件。该方法作为SPH边界处理的基本方法，被现今几乎所有的SPH模型采用或者改进［29-30］。一种比较典型的改进方式是在上述理论基础上增加粒子对称性条件，判断粒子周围水粒子分布对称性，辅助判断自由面身份。该方法由Khayyer等［31］提出，采用相邻粒子坐标积分来加强判定准确率，类似的改进方法则被Liu等［30］用于流固耦合问题模拟，采用临近粒子坐标的SPH积分代替直接求和，并发现该方法对结构物附近自由面粒子误判问题改善显著。

基于这些技术，有大量的SPH破碎波模拟研究成果发表。使用 WCSPH 模型，Colagrossi等［32-33］分别研究了复杂形状结构物附近的波浪破碎问题，而Shao［34］则进行了一系列的工作，探索ISPH在波浪破碎模拟上的表现。Khayyer等［35］通过改进上述自由面判别条件和处理，增加了粒子角动量守恒，较大地提高了自由面判定精度。

对破碎波模拟的另一重要方面是紊流模型的建立。SPH被提出至今已经近40年，即便是其在流体动力学中的大范围应用也已经有20多年了，但仍然没有很好地解决紊流问题。学术界普遍认为的原因是由于传统SPH方法离散精度比网格法稍低，存在一些无法消除的数值误差而掩盖了紊流的脉动速度和压强。这使得在更高阶精度的SPH格式成熟应用之前的紊流模拟研究都不会达到很高的精度，这在一定程度上限制了SPH紊流模型的发展。2002年，SPH的鼻祖Monaghan提出了一种求解N-S方程的垂向二维SPH紊流模型，虽然模拟结果相对较好，但随后就被指出其计算效率太低。Violeau等［36］提出了两种改进的代数紊流模型和一种基于雷诺应力模型的SPH紊流处理方法，在模拟类似溃坝的问题中确实能够表现出更高的精度，但仍然弱于同等情况下网格法中的精度。同时，Lo等［37］也提出了基于大涡模拟(LES)和k-ε模型的ISPH紊流模型，研究了破碎波复杂自由面附近紊动能分布。Shao等使用这些模型模拟了一系列的波浪、水流问题［38］，其中LES模型在后续的SPH中应用相对比较广泛［39］。

2.2 明渠流模拟

上述波浪模型有一个共同的特点，即在模拟区域并无流体进出，计算域整体的质量也是守恒的或者是方便进行拉格朗日描述的。由于SPH的拉格朗日特性，能够保证进口和出口区域粒子无进出且自动守恒，而对于需要在进出口设置入流和出流等欧拉描述的边界条件的情况，SPH则会遇到一些困难。在传统欧拉网格法中，进出口边界被固定地设置在一系列(通常是直线或者平面)网格的边界处，边界另一侧布置一层虚拟网格，用来给定入流或者出流速度、VOF等信息从而施加入流、出流边界条件，相对比较方便和直接。而在拉格朗日模型下施加这类欧拉边界条件则需要一些额外的处理。在水利工程中，常常需要进行一段河道或者部分实验水槽等实际问题的模拟，对这类入流出流边界处理的依赖较大，故有一些学者做了这方面的尝试。

为了解决边界处粒子进出问题，KAJTAR等［40］采用对进出口流体粒子的直接引入或删除来简易地实现进出流，但他们并未给出进出流处理的细节，而是主要着眼于流固耦合问题。一般认为，直接删除这些粒子导致的边界处的截断误差(truncation kernel error［41］)会非常大，且对于压力场的误算往往影响很大一块区域。Lastiwka等［42］稍后提出了一个更加细致的处理方法，但这种方法仅仅适用于气体动力学一类无自由表面的问题，显然对于明渠流这类自由面大变形问题并不适用。Federico等［43］于2012年给出了详细的进出流边界条件处理细节，在他们的模型中，进出流边界外布置了数层(据核函数而定)inflow/outflow粒子，进口边界粒子被施以指定的速度、压强并流入计算域，而出口边界粒子则保持原有速度和压力场不变。这些边界粒子只用于避免截断误差和施加边界条件，本身并不通过N-S方程求解，所以流出内部粒子支持域的边界粒子随即被删除而不会影响计算域内任何粒子。Federico使用改进的模型进行了明渠流的测试，随后应用于水跃问题，该方法能够较准确地实现恒定速度进口和自由出流边界条件。类似地，Bouscasse等［44］模拟了明渠流中圆柱后方流速场，而de Padova等［45］则模拟了三维的水跃问题。总的来讲，对于明渠流问题的模拟研究仍处于起步阶段，有很多的问题尚未完全解决。

2.3 浅水流动模拟

在模拟河流甚至是流域尺度问题时，基于N-S方程的数值模型将面临运算量过大的问题。对网格法而言，此时可以选择基于水深平均的方程求解，如浅水方程等。这些方程中流体的垂向参数被在水深方向积分，也就不存在自由面的捕捉问题，但一般认为SPH方法在浅水方程中的应用优势在于无需处理非线性对流项，简化了求解流程，且拉格朗日方法良好的守恒性在大尺度模拟时也更有优势。浅水SPH模型具有一些不同于N-S方程模型的独特性质，虽然拉格朗日水粒子体积在模拟过程中不会变化，但由于粒子显含水深维度的参数并随时间变化，故单个水粒子的水平方向参数(x或x、y)必然发生变化，导致会牵涉到一系列诸如粒子光滑长度的实时改变甚至是粒子分裂合并的问题，因而浅水SPH模型面临一些不同于N-S方程模型的困难。但由于大尺度问题对数值模拟技术的需求不断增加，所以相关浅水SPH模型的研究工作也得到了很大的发展。对于浅水SPH，较早的工作始于Bonet等［46］，他们首次提出了SPH-SWE(shallow water equation)模型，随后de Leffe等［47］也进行了类似研究。近年来，很多浅水SPH研究集中在对诸如开放边界［48］、固体结构物边界、源项［49］等问题上，逐步地实现了对真实河流地形和洪水的模拟［48］。Vacondio等［50］指出浅水SPH模型不同于传统N-S方程SPH模型的特性之一是在模拟中粒子的光滑长度往往要根据每个粒子的情况随位置发生变化。处理河流问题的浅水SPH与基于N-S方程的传统SPH具有很大区别，同时也有更多的问题有待解决。

3 SPH在水流-结构物相互作用问题中的应用研究

3.1 水流-固定结构物相互作用

3.1.1 固定刚性不透水结构物

除了对自由面的模拟外，另外一个SPH擅长的领域是对结构物的模拟。传统的欧拉网格法中对于拉格朗日描述的结构物的处理常常通过网格重建［51］或切割网格［52］等方法实现结构物、网格之间的信息传递;而SPH中对流体和结构物统一采用拉格朗日描述的特性避免了上述问题，拉格朗日流体粒子直接从结构物边界粒子处获得信息，施加边界条件。处理固体结构物的重点在于解决流体粒子在边界处的截断误差和施加边界条件。在早期模型中采用对边界附近粒子施加斥力［7］的办法阻止粒子溢出，但是边界斥力的值完全是根据粒子间距计算的，并无物理依据，这导致了边界处虽然不发生粒子溢出，但粒子分布与内部差异很大。现在比较流行的边界处理方法被称为虚拟边界粒子法，包括固定虚粒子和镜像粒子。前者采用几层固定于边界附近的虚粒子阻止流体溢出并施加边界条件。Randles等［53］是较早改进固定边界粒子用以施加边界条件的研究者之一，该技术随后被多位学者采用［29，54］。从理论上说，由于固定边界粒子携带了质量、压强、密度等比边界斥力更多的物理信息，所以比斥力法具有更大的灵活性，也能够获得更好的边界附近流体粒子分布，该方法至今还被广泛采用。而镜像粒子法则是从固定边界粒子法改进而来，镜像粒子位置根据边界另一侧的流体粒子实时更新获得，同时具有一定的速度，携带比固定边界粒子更多的信息。Morris等［55］是最早使用镜像粒子概念的研究人员，Cummins等［56］则在其成果中描述了详细的镜像规则，Yildiz等［57］将镜像粒子边界处理用于相对复杂的结构物表面的模拟，但在结构物角点处镜像粒子相互重叠，发生过度镜像的问题。Liu等［18］改进了这个方法，防止了角点过度镜像的发生，并发现其在模拟溃坝问题上能够获得更加稳定光滑的压力场。

被 SPH 模型研究最多的问题是溃坝问题［13，18，29，58］。实验中的溃坝水流撞击一侧挡墙后发生翻卷水舌并破碎，期间存在大量的不连续和液滴飞溅现象，这种情况在传统网格法中模拟起来相对困难，而SPH方法能够轻易地捕捉到水舌的翻卷和不连续的液滴飞溅。SPH压力场计算结果相对稳定光滑，对抨击点高压区的捕捉也比较准确。另一类SPH应用热点则集中在周期性的波浪和海工建筑物(防波堤、海墙等)相互作用问题，如Shao［59］模拟了波浪与幕墙式防波堤的相互作用，得到了波浪的形变和各种水动力学参数;而Gotoh等［60］则通过模拟波浪和非淹没式防波堤的相互作用，指出防波堤对波浪的主要耗散作用来自于附近漩涡的生成和脱落;相比之下，三维SPH模型也得到了一定程度的发展，且通常能够应用在更加符合物理实际情况的模拟，如Crespo等［61］开发了3D-SPH模型模拟更大尺度的长波与海岸结构物的作用，而Altomare等［62］则使用并行化的DualSPHysics模型模拟了armour block防波堤的波浪作用特性。在Altomare等的模拟中，由于在块体间会掺入大量的流体粒子，所以该模拟需要巨大的粒子数和极大的内存消耗，计算中，总计2146095个流体粒子被用于波浪模拟，有效的并行计算框架是保证这样规模的模拟能够顺利进行的必要条件。

3.1.2 多孔介质结构物

上述波浪与结构物相互作用模拟仅限于不透水刚体结构物，还有部分研究人员着眼于多孔介质与水流的相互作用，这些模拟技术在沙床、多孔防波堤以及植物区水流模拟等领域应用需求较大。广义上的多孔介质可以是错综排列的微小结构物(如植物区)，同网格法中的处理类似，SPH并不直接离散多孔结构物细节，而是通过求解空间平均后的N-S方程来反应多孔介质对水流平均流速的影响。这方面较早的工作始于日本京都大学的 Gotoh等［12］的MPS模型，他们通过在N-S方程中引入一个简单的摩擦项来描述多孔介质影响，模拟了孤立波在透水沙坡上的变形问题。通过改进上述方法，很多相关研究［63-65］对SPH多孔介质模型的建立进行了系统的验证。Shao［66］通过在多孔介质区域内外使用不同方程的方法提出了ISPH多孔介质模型，并模拟了周期波在浅滩破碎等问题。

3.2 水流-刚体结构物耦合运动

对于结构物运动和水流耦合计算的处理，分为运动学耦合和动力学耦合两大类，前者指运动结构物速度等参数事先给定，典型的如液箱晃荡问题等［67-69］;相比之下，动力学耦合模型除了需要精确的边界条件处理方法和运动边界的处理方法外，还需要计算运动结构物受力的方法。在传统的网格法中，对于拉格朗日结构物的描述需要一整套额外的自适应网格或者切割网格理论而显得复杂，尤其是当运动结构物处于自由面附近甚至穿透自由面时，处理起来会遇到更多的困难，实际工程中港口船舶的锚固、各种浮式结构物的波浪响应、水库滑坡体入水涌浪等问题却都涉及到运动结构物穿破自由面的情况，对数值模型提出了更高的要求。在现有的大部分SPH模型中，边界条件都是由虚拟边界粒子(固定粒子或镜像粒子)施加的，加上流体直接用拉格朗日粒子描述，所以它更加适合对上述流固耦合问题的模拟，通过引入简单的运动边界处理和加速度求解模块，即可实现流固耦合。流固耦合计算的主要难点是准确地求解物体受力及运动轨迹，对此，已有的SPH算法一般可分为动量传递法和积分表面压力法两大类。

3.2.1 动量传递法

早期的SPH模型(尤其是MPS模型)由于各种原因，对流体压强计算时会产生较大的数值压力震荡。在进行流固耦合模拟时需要积分结构物表面的压强从而获得压力的合力，这个合力也会发生震荡而影响后续加速度计算。MPS方法的创始人Koshizuka提出过一个巧妙的办法回避压强积分这个操作［70］，这个方法随后被 Shao 等［71-72］引入 SPH方法用于运动结构物的受力及运动模拟。在他们的模型中，结构物被离散成和固体密度一样的固体粒子群，其初始间距和流体粒子相同，在计算中，固体粒子与流体粒子一样参与求解N-S方程并获得速度，在N-S方程求解步骤之后，按实际固体密度积分所有固体粒子的动量，计入结构物运动的总动量，再根据结构物刚体假设重新分配这部分动量。在这样处理之后，计算中的总动量并没有丢失而保持守恒，结构物的刚体假设也能够保持。在下一时刻初始，经过刚体假设更新的固体粒子的位置会对周围流体施加作用。这个方法的应用并不广泛，能够查询到的文献只有Bøckmann等发表的，他们进行了一系列流固耦合问题的模拟研究［73］。这个方法整个流程中并未涉及对压强的积分，也就避免了因为压强震荡而导致的误差，它通过动量守恒特性，先使用固体粒子参与计算“收集”由流体传来的动量，随后根据刚体假设和获得的总动量等物理特性重新分配而更新结构物参数。但这个方法的一个缺点是结构物附近流体“不可压”特性会被破坏，且由于粒子质量的差异，结构物附近粒子分布也会更不规则［30］。

3.2.2 积分力法

相比之下，另一种方法应用更加广泛，通过积分结构物附近(表面)粒子的压强(或结构物边界粒子压强)求得压力，结合结构物受力情况依据牛顿第二定律计算结构物加速度和角加速度，时间积分平动速度和转动速度及位置。理论上讲，牛顿第二定律本质也是动量守恒，故和第一种方法相比具有相同的理论基础，不过对流体力直接积分的方法更加符合流固耦合的物理实际。考虑到计算流程的简洁，大部分SPH流固耦合模型建立在WC-SPH下，如Vandamme等［74］采用固定边界粒子(见 3.1.1节)方法的WCSPH模型模拟了楔形体入水和圆柱入水问题;而Oger等［75］则基于镜像粒子继续了流固耦合的研究，他们引入了镜像边界附近压强积分的算法，详细介绍了流固耦合的处理并给出了一系列运动结构物的模拟结果，展示了模型的应用。Liu等［30］综合前人成果改进了一套用于处理流固耦合的ISPH模型，模拟了液箱晃荡和方柱入水等问题。液箱晃荡问题中，他们使用了真实的箱体运动外激励速度晃荡箱体，代替在网格法中普遍使用的非惯性坐标系［76］，模拟结果和实验对比良好，随后他们将其应用于更加复杂的方柱入水问题中，该问题结构物运动是包含平动和转动的全自由运动，ISPH结果较好地还原了方柱各个自由度的运动和自由面的变化，模拟结果和FLUENT软件计算结果相差不大。在水利工程实际问题方面，笔者也将该模型拓展到对滑坡体涌浪的产生、传播和翻坝问题的模拟上，并通过实验验证了SPH结果的精度，展示了ISPH在水利工程问题中的应用潜力。

3.3 水流-柔性结构物模拟

上述运动结构物的处理算法都仅限于刚体结构物，对于可变形的弹性/柔性结构物来说，情况则更加复杂一些。如前面所述，基于网格的数值方法在处理柔性物体时会面临比刚体更严重的困难。实际计算中对弹性结构物的模拟通常使用FEM方法或者类似的拉格朗日法，这使得原本就基于拉格朗日描述的SPH具有更大的理论优势。最早使用SPH模拟上述问题的是Müller等［77］，在他们的方法中，基于Lennard-Jones势引入边界力来模拟弹性结构物，但他们的模拟结果在弹性边界处会发生严重的粒子分布不均匀现象，精度较差。Lenaerts等［78］则模拟了柔性外壳与水的相互作用，他们采用了边界力和位置修正联合作用的办法防止固体与流体粒子的互相穿透，但也指出粒子穿透现象在材料弹性模量较小时很难避免。Amini等［79］则系统研究了多种柔性/弹性结构物的流固耦合问题，这些类似问题在传统网格法中是很难模拟的。

4 SPH在泥沙输运问题中的应用

水流中的泥沙输运问题是水利和海洋工程中极其重要的一大类问题。早期的关于泥沙的数值模拟研究都是基于网格法的，如基于欧拉网格的有限差分法(FDM)［80］和有限体积法(FVM)［81］，以及基于拉格朗日网格的有限元法(FEM)［82］。另外还有一种特殊的网格法，即格子玻尔兹曼方法(LBM)［83］。由于SPH的粒子特性，能够基于泥沙颗粒或者粒子底床等概念来描述及模拟泥沙问题，具有一定的优势，目前SPH方法对泥沙的处理模型大致分为泥沙颗粒两相流模型和底床高度冲刷模型。

4.1 SPH泥沙颗粒两相流模型

利用拉格朗日粒子概念描述非连续的泥沙颗粒具有一定的优势，它可以实现泥沙-水流的“双向耦合”，即在泥沙被水流冲蚀的同时，水流也会受到泥沙的反作用，且这一切是由SPH水粒子和泥沙粒子混合物自动完成的，不需太多的额外处理。Monaghan等［84］最早将这样的处理加入到SPH方法中，而Shakibaeinia等［85］使用了类似SPH的方法改进MPS方法研究了溃坝水流导致的泥沙输运和底床变形。在这些研究中，底床被可冲蚀的泥沙粒子代替，当底床剪力达到设定的临界值时，泥沙粒子开始启动，启动后的泥沙颗粒作为大密度的流体粒子参与流体计算，泥沙颗粒的沉降也完全按照其真实受力特性计算。

4.2 SPH底床高度冲刷模型

基于泥沙和水流的两相流模型具有处理简单、边界条件易于施加等特点，且实践证明其在具有底床大变形和水沙混掺剧烈的问题中(如3.1.1节中的溃坝)模拟精度良好、可靠性高。但单个泥沙颗粒的启动模型存在一定的缺陷:首先，在实际模拟中SPH泥沙粒子体积远远大于真实的沙粒，故必须等待底床冲蚀积累到一定量之后才能启动，在模拟中造成了显著的非连续底床变动，对于精细且缓慢的动床变化计算误差较大;其次，游离在水中的单个泥沙粒子的水动力学特性和它所表征的一团细颗粒泥沙差别较大，而将其作为一个整体通过自重沉降回底床的做法可能带来很大误差。从现有研究结果来看，大多数SPH泥沙输运模拟均只关注了比较剧烈的诸如溃坝、抨击等问题的泥沙输运，可能也与上述缺陷有一定关系。因此，Krištof等［86］尝试在 SPH 中引入类似网格法的连续底床冲刷模型来实现连续的底床变形和推移质输运。在他们的方法中，底床仍然用泥沙粒子来表征形变和施加边界条件，但该泥沙颗粒并不启动，而是根据水流条件实时计算冲蚀并由推移质输沙率和悬起沉降率共同保证质量守恒。同时，该方法还依据类似网格法中的处理加入了基于拉格朗日描述的悬移质对流扩散方程，将悬移质离散到水粒子系统中并通过拉格朗日粒子悬移质浓度概念求解对流扩散方程，实现真实的泥沙输运和底床变形模拟。

5 SPH气液两相流模拟应用研究

SPH方法在流体模拟中的另一个应用领域是对于多相流、特别是气液两相流的模拟［87］。由于传统单相SPH认为空气压强为常数，造成了当水流一侧压强过大时会将自由面向空气一侧推动，这在描述水中气泡时显然是错误的。气液两相SPH模型的最初并未局限于气液两相而是指广义上的两种不同密度流体的耦合模拟，Monaghan等［88］最早将两相流概念引入SPH，他们使用了两种不同粒子描述不同流相，相间的相互作用根据N-S方程计算而两相界面处则增加一个非物理的作用力实现耦合，用来模拟dust gas流。后续研究将其应用在很多领域［32，36，84］。这些早期的 SPH 方法均基于微可压缩模型，在处理密度差很大的两相流交界面时会遇到一些困难。如Colagrossi等［32］曾指出WCSPH在模拟大密度比的两相流时，需要在交界面处做非常大的限制(一个密度额外修正处理、一个施加在低密度流体上的增大表面张力的处理、一个速度光滑处理等)才能保证精度，并且需要在两种流体间设置差异极大的声速来均衡WCSPH中声速带来的影响。这些做法除了严重增加算法的非物理限制外，还大大增加了运算量，WCSPH方法中过多的人工耗散和修正的弊端在两相流问题面前暴露无遗。

不可压缩模型因为不需要采用人工声速，可以较好地处理水气两相流，其中比较有代表性的是Hu等［89-90］的一系列研究成果，其ISPH模型被成功地用于模拟大密度差的两相流问题。Shao［91］则提出了类似的ISPH模型，将其同时用于大密度差和小密度差两相流的模拟，展示了ISPH模型良好的通用性和稳定性。在Shao的模型中，依然采用两种粒子来分别模拟两相(如水和气)，各相分别独自求解一套N-S方程描述流体粒子的运动，而交界面上则使用另一相粒子位置充当边界粒子、施加边界条件，这样的处理类似多孔介质［66］甚至是运动结构物中的处理。这种ISPH模型在计算过程中不发生密度变化，保证了交界面两侧各自流体密度的恒定。现有的大部分ISPH模型对于气体相都处理成不可压流体，这在某些情况下并不符合物理实际，需要改进模型。

两相流问题中，虽然可以忽略水的压缩性，但气体的压缩性是不能忽略的，如水中大气泡的运动、破碎波中卷积的气囊等。这就需要采用不同的模型模拟水体和气体，目前这方面的工作还很少，只在近两年，新加坡国立大学Luo等通过粒子法实现了对可压缩气体与不可压缩液体的模拟。另一个水动力学中常见的问题是当水面紊动很强时，会在水面发生自掺气，可掺入的气泡又会在紊流变弱时逸出水面。这个问题的模拟尚未见报道。Luo等①LUO M，KOH C G，GAO M，et al.A particle method for two-phase flows with large density difference［J］.International Journal for Numerical Methods in Engineering，2015(已录用，待刊).基于两相CPM(consistent particle method，基于 MPS改进)方法改进了模型，使之能够处理晃荡问题中的气腔形变，模拟结果和实验吻合良好。

6 研究展望

通过分析发现，SPH具有稳定简洁的特点，在交界面大变形问题的模拟上具有优势，而在未来则会有更多的关于泥沙输运和水流掺气的SPH研究成果出现。而复杂的多物理过程的研究也将成为热点，如需要考虑材料特性的柔性体变形模拟，需要考虑破坏土体特性的滑坡体入水问题，以及需要考虑固液气三相的泥石流运动问题等。

在对纯水流问题与流固耦合问题的模拟中，SPH已经取得了很大的进展，计算方法相对成熟稳定，已经可以与传统的网格计算法分庭抗礼，而在交界面大变形(自由面和结构物表面)问题中，SPH还略占优势。但在泥沙模拟和气液两相流模拟方面，尤其是动态掺气过程的模拟中，SPH方法才刚刚起步，尚有大量工作需要逐步推进。考虑到SPH方法在处理多相流方面具有自身的优势，可以预见在未来会有越来越多的关于泥沙输运和水流掺气的SPH研究成果出现。同时采用SPH方法模拟复杂的多物理过程问题的研究也将出现，这些问题包含需要考虑材料特性的柔性体变形模拟，需要考虑破坏土体特性的滑坡体入水问题，以及需要考虑固液气三相的泥石流运动问题等。

未来SPH在水利与海洋工程研究领域的发展将主要集中在以下几个方面:①SPH计算效率与精度的提高。目前大部分成熟的水动力学SPH算法均基于使用各向同性均匀分布粒子的核函数，在计算多尺度物理问题时，面临计算效率低的问题。在结合N-S方程求解方法的基础上，开发出高阶精度的具有椭圆支持域的非均匀各向异性核函数的SPH模型，可以有效提高模型性能，用于模拟狭长计算域问题并提高局部模拟精度(如长波与结构物作用问题)，是 SPH未来的重要的理论创新方向之一。②在高精度SPH模型开发基础上，发展适用于SPH特点的紊流模型，如大涡模拟，更好地模拟各类实际问题，可以大大拓宽SPH的工程应用地位，摆脱目前SPH多用于算法研究的困境。③为了克服SPH相比网格法更加耗时的弱点，对三维SPH模型实现并行计算是必要的发展趋势。在这些研究成功的基础上，SPH方法将获得更为广阔的发展空间，并展现出其独特的优势与迷人的前景。

［1］LIU M B，LIU G R.Smoothed particle hydrodynamics(SPH):an overview and recent developments［J］.Archives of Computational Methods in Engineering，2010，17:25-76.

［2］GINGOLD R A，MONAGHAN J J.Smoothed particle hydrodynamics:theory and application to non-spherical stars［J］.Montly Notices of the Royal Astronomical Society，1977，181(3):375-389.

［3］LUCY L B.A numerical approach to the testing of fission hypothesis［J］.The Astronomical Journal，1977，82(12):1013-1024.

［4］LEE E S，MOULINEC C，XU R，et al.Comparisons of weakly compressible and truly incompressible algorithms for the SPH mesh free particle method［J］.Journal of Computational Physics，2008，227:8417-8436.

［5］MORRIS J P，ZHU Y，FOX P J.Parallel simulation of pore-scale flow through porous media［J］.Computers and Geotechnics，1999，25:227-246.

［6］LIU G R，LIU M B.Smoothed particle hydrodynamics:a meshfree particle method［M］.Singapore:World Scientific，2003.

［7］MONAGHAN J J.Simulating free surface flows with SPH［J］.Journal of Computational Physics，1994，110(2):399-406.

［8］LIN P Z，XU W L.NEWFLUME:a numerical water flume for two-dimensional turbulent free surface flows［J］.Journal of Hydraulic Research，2006，44(1):79-93.

［9］高睿，任冰，王国玉，等.孤立波浅化过程的SPH数值模拟［J］.水动力学研究与进展:A 辑，2010，25(5):620-629.(GAO Rui，REN Bing，WANG Guoyu，et al.SPH model of solitary waves shoaling on a mild sloping beach［J］.Chinese Journal of Hydrodynamics，2010，25(5):620-629.(in Chinese))

［10］刘谋斌，宗智，常建忠.光滑粒子动力学方法的发展与应用［J］.力学进展，2011，41(2):217-234.(LIU Moubin，ZONG Zhi，CHANG Jianzhong.Development and application of smoothed particle hydrodynamics［J］.Advances in Mechanics，2011，41(2):217-234.(in Chinese))

［11］郑兴，段文洋.K2_SPH方法及其对二维非线性水波的模拟［J］.计算物理，2011，28(5):659-666.(ZHENG Xing，DUAN Wenyang.K2_SPH method and application for 2D nonlinear water wave simulation［J］.Chinese Journal of Computational Physics，2011，28(5):659-666.(in Chinese))

［12］GOTOH H，SAKAI T，Lagrangian simulation of breaking waves using particle method［J］.Coastal Engineering Journal，1999，41(3/4):303-326.

［13］KHAYYER A，GOTOH H.Modified moving particle semiimplicit methods for the prediction of 2D wave impact pressure［J］.Coastal Engineering，2009，56(4):419-440.

［14］DAO MH，XU H，CHAN ES，TKALICH P.Modelling of tsunami-like wave run-up，breaking and impact on a vertical wall by SPH method［J］.Natural Hazards and Earth System Sciences，2013，13:3457-3467.

［15］ZHENG J H，SOE M M，ZHANG C，et al.Numerical wave flume with improved smoothed particle hydrodynamics［J］.Journal of Hydrodynamics:Ser B，2010，22(6):773-781.

［16］GALLATI M，BRASCHI G，FALAPPI S.SPH simulations of the waves produced by a falling mass into a reservoir［J］.Nuovo Cimento Societa Intalian di Fisica C:Geophysics and Space Physics，2005，28:129.

［17］DU X，W U W，GONG K，et al.Two dimensional SPH simulation of waterwaves generated by underwater landslide［J］.Journal of Hydrodynamics:Ser A，2006，21(5):579-586.

［18］LIU X，XU H H，SHAO S D，et al.An improved incompressible SPH model for simulation of wave-structure interaction［J］.Computers and Fluids，2013，71:113-123.

［19］FRIGAARD P，BRORSEN M.A time domain method for separating incident and reflected irregular waves［R］.Aalborg，Danish:Hydraulics ＆ CoastalEngineering Laboratory，Department of CivilEngineering，Aalborg University，1993.

［20］HAYASHI M，GOTOH H，MEMITA T，et al.Gridless numerical analysis of wave breaking and overtopping at upright seawall［C］//Coastal Engineering Conference.San Francisco:Asce American Society of Civil Engineers，2001:2100-2113.

［21］DIDIER E，NEVES M G.A semi-infinite numerical wave flume using smoothed particle hydrodynamics［J］.International Journal of Offshore and Polar Engineering，2012，22(3):193-199.

［22］LARSEN J，DANCY H.Open boundaries in short wave simulations-a new approach［J］.Coastal Engineering，1983，7(3):285-297.

［23］LEE C，SUH K D.Internal generation of waves for timedependent mild slope equations［J］.Coastal Engineering，1998，34(1):35-57.

［24］LIN P，LIU P L F.Internal wave-maker for Navier-Stokes equations models［J］.Journal of Waterway，Port，Coastal and Ocean Engineering，1999，125(4):207-215.

［25］WEI G，KIRBY J T，SINHA A.Generation of waves in Boussinesq models using a source function method［J］.Coastal Engineering，1999，36(4):271-299.

［26］CHOI J，YOON S B. Numerical simulations using momentum source wave-maker applied to RANS equation model［J］. Coastal Engineering，2009，56 ( 10 ) : 1043- 1060.

［27］HA T，LIN PZ，CHO Y S.Generation of 3D regular and irregular waves using Navier-Stokes equations model with an internal wave maker［J］.Coastal Engineering，2013，76:55-67.

［28］LIU X，LIN PZ，SHAOSD.ISPH wave simulation by using an internal wave maker［J］.Coastal Engineering，2015，95:160-170.

［29］SHAO S D，LO E Y M.Incompressible SPH method for simulating Newtonian and non-Newtonian flows with a free surface［J］.Advances in Water Resources，2003，26:787-800.

［30］LIU X，LIN P，SHAO S D.An ISPH simulation of coupled structure interaction with free surface flows［J］.Journal of Fluids and Structures，2014，48:46-61.

［31］KHAYYERA，GOTOH H，SHAOSD.Enhanced predictions of wave impact pressure by improved incompressible SPH methods［J］.Applied Ocean Research，2009，31:111-131.

［32］COLAGROSSI A，LANDRINI M.Numerical simulation of interfacial flows by smoothed particle hydrodynamics［J］.Journal of Computational Physics，2003，191(2)，448-475.

［33］MONAGHAN J J，KOS A，ISSA N.Fluid motion generated by impact［J］.Journal of Waterway，Port，Coastal and Ocean Engineering，2004，129:250-259.

［34］SHAO S D.Incompressible SPH simulation of wave breaking and overtopping with turbulence modelling［J］.International Journal for Numerical Methods in Fluids，2006，50(5):597-621.

［35］KHAYYERA，GOTOH H，SHAO SD.Corrected Incompressible SPH method for accurate water-surface tracking in breaking waves［J］.Coastal Engineering，2008，55(3):236-250.

［36］VIOLEAU D，BUVAT C，ABED-MERAÏM K，et al.Numerical modelling of boom and oil spill with SPH［J］.Coastal Engineering，2007，54(12):895-913.

［37］LO E Y M，SHAO S D.Simulation of near-shore solitary wave mechanics by an incompressible SPH method［J］.Applied Ocean Research，2002，24(5):275-286.

［38］SHAO S D，JI C，GRAHAM D I，et al.Simulation of wave overtopping by an incompressible SPH model［J］.Coastal Engineering，2006，53(9):723-735.

［39］DALRYMPLE R A，ROGERS B D.Numerical modeling of waterwaves with the SPH method［J］.Coastal Engineering，2006，53(2):141-147.

［40］KAJTAR J，MONAGHAN JJ.SPH simulationsof swimming linked bodies［J］.Journal of Computational Physics，2007，227:8568-8587.

［41］ANDREA A，JEAN-CHRISTOPHE M，FRANCIS L，et al.SPH truncation error in estimating a 3D function［J］.Computers and Fluids，2011，44:279-296.

［42］LASTIWKA M，BASA M，QUINLAN N J.Permeable and non-reflecting boundary conditions in SPH ［J］.International Journal for Numerical Methods in Fluids，2009，61:709-724.

［43］FEDERICO I，MARRONE S，COLAGROSSI A，et al.Simulating 2D open-channel flows through an SPH model［J］.European Journal of Mechanics-B/Fluids，2012，34:35-46.

［44］BOUSCASSE B，COLAGROSSI A，MARRONE S，et al. Viscous flow past a circular cylinder below a free surface ［C］/ /ASME 33rd International Conference on Ocean，Offshore and Arctic Engineering. San Francisco: merican Society of Mechanical Engineers ( AMSE ) ，2014: V002T08A080-V002T08A080.

［45］de PADOVA D，MOSSA M，SIBILLA S，et al.3D SPH modelling of hydraulic jump in a very large channel［J］.Journal of Hydraulic Research，2013，51(2):158-173.

［46］BONET J，KULASEGARAM S，RODRIGUEZ-PAZ M X，Variational formulation for the smooth particle hydrodynamics(SPH)simulation of fluid and solid problems［J］.Computer Methods in Applied Mechanics and Engineering，2004，193(12):1245-1256.

［47］de LEFFE M，LE TOUZÉ D，ALESSANDRINI B.SPH modeling of shallow-water coastal flows［J］.Journal of Hydraulic Research，2010，48(Sup1):118-125.

［48］VACONDIO R，ROGERS B D，STANSBY P K，et al.SPH modeling of shallow flow with open boundaries for practical flood simulation［J］.Journal of Hydraulic Engineering，2011，138(6):530-541.

［49］VACONDIO R，ROGERS B D，STANSBY P K，et al.Shallow water SPH for flooding with dynamic particle coalescing and splitting［J］.Advances in Water Resources，2013，58:10-23.

［50］VACONDIO R，ROGERS B D，STANSBY P K.Accurate particle splitting for smoothed particle hydrodynamics in shallow water with shock capturing［J］.International Journal for Numerical Methods in Fluids，2012，69(8):1377-1410.

［51］HASSAN O，PROBERT E J，MORGAN K，et al.Unsteady flow simulation using unstructured meshes［J］.Computer Methods in Applied Mechanics and Engineering，2000，189:1247-1275.

［52］LAI M C，PESKIN C S.An immersed boundary method with formal second-order accuracy and reduced numerical viscosity［J］.Journal of Computational Physics，2000，160:705-719.

［53］RANDLES P W，LIBERSKY L D.Smoothed particle hydrodynamics:some recent improvements and applications［J］.Computer Methods in Applied Mechanics and Engineering，1996，139:375-408.

［54］LIU G R，GU Y T.A local radial point interpolation method(LRPIM)for free vibration analyses of 2D solids［J］.Journal of Sound and Vibration，2001，246(1):29-46.

［55］MORRIS J P，FOX P J，ZHU Y.Modeling low Reynolds number incompressible flows using SPH［J］.Journal of Computational Physics，1997，136:214-226.

［56］CUMMINS SJ，RUDMAN M，An SPH projection method.Journal of Computational Physics［J］.1999，152:584-607.

［57］YILDIZ M，ROOK R A，SULEMAN A.SPH with the multiple boundary tangentmethod［J］.International Journal for Numerical Methods in Engineering，2009，77:1416-1438.

［58］王亚萍，光滑粒子流体动力学(SPH)方法在溃坝问题中的应用研究［D］.天津:天津大学，2012.

［59］SHAO S D.SPH simulation of solitary wave interaction with a curtain-type breakwater［J］.Journal of Hydraulic Research，2005，43(4):366-375.

［60］GOTOH H，SHAO S D，MEMITA T.SPH-LES model for numerical investigation of wave interaction with partially immersed breakwater［J］.Coastal Engineering Journal，2004，46(1):39-63.

［61］CRESPO A J C，GÓMEZ-GESTEIRA M，DALRYMPLE R A.3D SPH simulation of large waves mitigation with a dike［J］.Journal of Hydraulic Research，2007，45(5):631-642.

［62］ALTOMARE C，CRESPO A J C，ROGERS B D，et al.Numerical modelling of armour block sea breakwater with smoothed particle hydrodynamics［J］.Computers ＆Structures，2014，130:34-45.

［63］TARTAKOVSKY A M，MEAKIN P.Pore scale modeling of immiscible and miscible fluid flowsusing smoothed particle hydrodynamics［J］.Advances in Water Resources，2006，29(10):1464-1478.

［64］TARTAKOVSKY A M，MEAKIN P，SCHEIBE T D，et al.Simulations of reactive transport and precipitation with smoothed particle hydrodynamics［J］.Journal of Computational Physics，2007，222(2):654-672.

［65］TARTAKOVSKY A M. Langevin model for reactive transport in porous media［J］. Physical Review E，2010， 82( 2) : 026302.

［66］SHAO S D.Incompressible SPH flow model for wave interactions with porous media［J］.Coastal Engineering，2010，57(3):304-316.

［67］陈臻.SPH算法改进及在晃荡与入水中的应用［D］.大连:大连理工大学，2014.

［68］李大鸣，陈海舟，张建伟，等.基于SPH法的二维矩形液舱晃荡研究［J］.计算力学学报，2010，27(2):369-374.(LI Daming，CHEN Haizhou，ZHANG Jianwei，et al.A study of 2D liquid sloshing in rectangle tanks based on SPH method［J］.ChineseJournal of Computational Mechanics，2010，27(2):369-374.(in Chinese))

［69］崔岩，吴卫，龚凯，等.二维矩形水槽晃荡过程的SPH方法模拟［J］.水动力学研究与进展:A 辑，2008，23(6):618-624.(CUI Yan，WU Wei，GONG Kai，et al.Numerical simulation of sloshing in two dimensional rectangular tanks with SPH ［J］.Chinese Journal of Hydrodynamics，2008，23(6):618-624.(in Chinese))

［70］KOSHIZUKA S，NOBE A，OKA Y.Numerical analysis of breaking waves using the moving particle semi-implicit method［J］.International Journal for Numerical Methods in Fluids，1998，26(7):751-769.

［71］SHAO S D，GOTOH H，Simulating coupled motion of progressive wave and floating curtain wall by SPH-LES model［J］.Coastal Engineering Journal，2004，46(2):171-202.

［72］SHAO S D.Incompressible SPH simulation of water entry of a free-falling object［J］.International Journal for numerical methods in fluids，2009，59(1):91-115.

［73］BØCKMANN A，SHIPILOVA O，SKEIE G.Incompressible SPH for free surface flows［J］.Computers ＆ Fluids，2012，67:138-151.

［74］VANDAMME J，ZOU Q P，REEVE D E.Modeling floating object entry and exit using smoothed particle hydrodynamics［J］.Journal of Waterway，Port，Coastal，and Ocean Engineering，2011，137(5):213-224.

［75］OGER G，DORING M，ALESSANDRINI B，et al.Twodimensional SPH simulations of wedge water entries［J］.Journal of Computational Physics，2006，213:803-822.

［76］XUE M A，LIN P Z.Numerical study of ring baffle effects on reducing violent liquid sloshing［J］.Computers ＆Fluids，2010，52:116-129.

［77］MÜLLER M，SCHIRM S，TESCHNER M，et al.Interaction of fluids with deformable solids［J］.Computer Animation and Virtual Worlds，2004，15(3/4):159-171.

［78］LENAERTS T，DUTRÉ P.Unified SPH model for fluidshell simulations［R］.Leuven，Belgium:Department Computer Science，Katholieke Universiteit Leuven，2008.

［79］AMINI Y，EMDAD H，FARID M，A new model to solve fluid-hypo-elastic solid interaction using the smoothed particle hydrodynamics(SPH)method［J］.European Journal of Mechanics-B/Fluids，2011，30(2):184-194.

［80］ROACHE P J. Computational fluid dynamics ［M］.Albuquerque，New Mexico: Hermosa Publishers， 1998.

［81］EYMARD R，GALLOUET T，HERBIN R.Finite volume methods［M］.Wroclaw，Poland:University of Wroclaw，2008.

［82］GLOWINSKI R.Finite element methods for incompressible viscous flow［M］.Amsterdam:Elsevier Science B.V.，2003.

［83］THORNE D T，MICHAEL C.Lattice Boltzmann modeling:an introduction for geoscientists and engineers［M］.Berlin:Springer，2006.

［84］MONAGHAN J J，CAS RAF，KOS A M，et al.Gravity currents descending a ramp in a stratified tank［J］.Journal of Fluid Mechanics，1999，379:39-69.

［85］SHAKIBAEINIA A，JIN Y C.A mesh-free particle model for simulation of mobile-bed dam break［J］.Advances in Water Resources，2011，34(6):794-807.

［86］KRIŠTOF P，BENEŠ B，KRˇIVÁNEK J，et al.Hydraulic erosion using smoothed particle hydrodynamics［J］.Computer Graphics Forum，2009，28(2):219-228.

［87］龚凯，刘桦，王本龙.两相 SPH 方法［C］//吴有生，刘桦，许维临，等.第九届全国水动力学学术会议暨第二十二届全国水动力学研讨会文集.北京:海洋出版社，2009:340-345.

［88］MONAGHAN J J，KOCHARYAN A.SPH simulation of multi-phase flow［J］.Computer Physics Communications，1995，87(1):225-235.

［89］HU X Y，ADAMS N A.An incompressible multi-phase SPH method［J］.Journal of Computational Physics，2007，227(1):264-278.

［90］HU X Y，ADAMS N A，A constant-density approach for incompressible multi-phase SPH ［J］.Journal of Computational Physics，2009，228(6):2082-2091.

［91］SHAO S D. Incompressible smoothed particle hydrodynamics simulation of multifluid flows［J］.International Journal for Numerical Methods in Fluids，2012，69(11):1715-1735.