SPH 理论和方法在高速水动力学中的研究进展

钟诗蕴，孙鹏楠*,2，吕鸿冠,2，彭玉祥，张阿漫

1 中山大学 海洋工程与技术学院，广东 珠海 519082

2 南方海洋科学与工程广东省实验室，广东 珠海 519082

3 哈尔滨工程大学 船舶工程学院，黑龙江 哈尔滨 150001

0 引 言

面向海洋强国的战略需要，船舶与海洋工程得到了大力发展，与之相关的结构力学和水动力学等领域均产出了丰硕成果[1-3]。其中，高速舰船兴波[4]、航行体（器）出水与入水[5]、水下爆炸与结构毁伤[6-7]等典型高速水动力学问题与现代舰船及其武器装备的综合性能密切相关，但这些问题都具有大变形、动边界、强对流、多介质等复杂特征，而运用传统理论、试验或数值模拟方法较难得到彻底解决，因此，近年来成为了研究热点和难点[8-9]。

高速水面航行器的兴波、飞机水上滑行等过程涉及了复杂曲面结构物与自由液面或水气界面的复杂耦合作用，剧烈的液固耦合运动易产生非线性翻卷和破碎波浪[10]，甚至会引起液滴强烈飞溅，对航行器的结构及内部仪器零件的安全性构成威胁，且诱导的噪声[11]和尾迹[12-13]对航行器隐蔽性极为不利，精准捕捉兴波水花飞溅和预报砰击载荷关系到航行器的隐身性、安全性设计。

在军民用领域中，抛落式救生艇入水、飞行器海上回收，鱼雷空投入水、潜射导弹出水等航行器或结构物高速出、入水的问题[14-15]涉及瞬态强冲击载荷、结构物大幅运动、气泡生长与溃灭等复杂过程[16]，是典型的多相流−固耦合问题[17]。这种高速跨介质过程所产生的巨大砰击载荷，一方面可能造成装备结构损坏、内部零件失灵，另一方面也可能引起航行器弹道失稳、投放失效等。因此，研究航行器的高速出、入水问题对于跨介质武器装备等的研发有着重要意义。

研究表明，水下爆炸冲击波峰值及冲量相比空中爆炸有明显提高[18]。舰船遭受近场水下爆炸时，冲击波、高速射流和破片等首先对舰船造成局部毁伤[19]，随后的气泡脉动和射流冲击可能造成舰船整体失效[7]。近场水下爆炸导致船体出现破口以后，剧烈涌流还可能会进一步造成船体结构损坏，最终导致船体沉没。考虑到水下爆炸过程的复杂性和试验研究的高成本，研发精准的数值计算方法，预报水下爆炸载荷、结构毁伤和舰船破损沉没过程成为了各国海军研究的重点[19]。

对于上述高速航行器的兴波、航行体出水和入水、水下爆炸与结构毁伤等工程问题，其背后的基础力学问题离不开对复杂运动结构物与多相流体之间的流−固耦合过程的求解。研究人员针对不同场景下的力学过程进行简化和凝练，提出了多种理论模型[20-21]，解释了气泡动力学或简单构型结构物与流体作用过程背后的力学机理。然而，理论模型对于复杂的工程实际问题适用性有限。而模型试验虽然是获得可靠结果的重要手段[4,22]，但开展试验通常需要较大的场地和昂贵的设备，尤其是在高速水动力问题的模型试验中，尺度效应明显[23]，相关物理参数测量困难，试验数据的可重复性较差，试错成本也高。

近年来，随着计算机技术的发展，涌现出了许多数值计算方法。目前，基于不同的离散和求解形式，数值模拟方法可分为有网格法和无网格法[24]。前者在解决大变形流−固耦合问题时容易受到网格畸变的困扰；后者得益于天然的拉格朗日特性及逐渐完备的粒子近似理论优势，在模拟大变形问题时不受边界变形量的限制[25]，故在高速水动力问题中得到了广泛应用[24]。典型的无网格法包括离散元法（discrete element method, DEM）[26]、光滑粒子流体动力学（smoothed particle hydrodynamics, SPH）方法[27-29]、移动粒子半隐式（moving particle semi-implicit）方法[30]等。MPS 方法和SPH 方法在求解水动力问题方面难分伯仲[31-34]，但MPS 方法基于完全不可压缩假设[35-36]，难以应用于压缩效应明显的问题，例如水下爆炸的气泡脉动过程等；而DEM 是离散元算法[37]，主要用于处理离散粒子（例如沙砾、冰块等颗粒流）大变形问题[38-40]，并不适用于高速水力学研究。相较而言，SPH 理论和方法在处理流体大变形和可压缩问题时的优势明显，便于开展并行计算，工程适用性强，更适合于求解大型高速水动力问题。

海洋结构物在高速航行或跨介质运动时，对自由液面以下的单相流动区域，准确模拟湍流边界层是关键[41]。针对这方面的应用，有网格法因其湍流模型较为成熟、便于采用多分辨率网格技术而受到研究人员的青睐[42]。但对于自由液面或多相界面附近区域，因存在流体高速砰击引起的液滴飞溅和水气掺混等过程，有网格法的精度容易受到网格尺度的极大限制[24]。例如，一旦液滴飞溅到网格稀疏区域，易引起流体质量不守恒造成的液滴飞溅消失。相比而言，SPH 无网格粒子法用于模拟液滴飞溅问题更具优势[4,22]。该方法基于核近似、粒子近似和虚粒子边界等基本原理，辅之以日益完善的高精度修正理论和技术，在模拟运动边界和大变形界面等问题时具有较强计算能力和广阔发展前景。

尽管SPH 方法在计算的全过程中能够精确满足质量守恒要求并擅长捕捉自由液面翻卷、破碎和喷溅等复杂现象，但运用传统SPH 方法计算的压力载荷一般伴随严重的数值噪声[43]，涉及复杂区域时还容易引起张力不稳定性和数值空洞等[44]。为解决这些问题，经过近半个世纪的不断发展，SPH 方法在理论和数值技术方面得到了不断优化，计算精度和稳定性获得巨大提升，在新一代GPU 并行技术的加持下，显著提高了计算效率[45]，具备了求解工程问题的能力。并且，随着嵌入SPH的多算法耦合技术的出现和发展[6,46-47]，SPH 理论和方法的优势得以充分发挥，使其为复杂多元的高速水动力学问题提供了新的高精度、高效率求解方案。

综上所述，为进一步阐述SPH 理论和方法在船海高速水动力学相关问题中的研究进展，本文将从SPH 的基本原理及数值技术研究进展出发，围绕高速水面航行器的水动力以及高速跨介质航行器水动力问题、水下爆炸与结构毁伤等方面展开综述，并对SPH 理论和方法在求解高速水动力问题中的未来发展趋势进行展望，以期为涉及上述问题的相关研究提供参考。

1 SPH 理论及技术进展

1.1 SPH 基本理论与算法

1.1.1 SPH 基本理论

SPH 方法最初是为研究天体动力学而提出的[48]，上世纪90 年代开始应用于水动力问题的模拟[49]。以物体入水问题为例，图1 给出了基于流体体积（volume of fluid，VOF）的欧拉网格法[50]及拉格朗日描述的SPH 方法[51-52]模拟圆柱体及楔形体入水以及液面飞溅的过程。在欧拉网格法的计算结果中，可见喷溅射流尖端存在非物理性流体消失，水气界面较模糊的情况，而SPH 方法中粒子可以精确地捕捉到流体飞溅特征。上述结果体现了SPH 方法在处理大变形问题方面的优势，因此被广泛应用于船海工程领域[5,53]。

从原理上看，SPH 方法基于核函数近似理论计算空间导数，可对特定物理场的梯度、散度等采用积分形式进行计算[54]；在数值实现方面，SPH方法一般将流体域离散成紧凑的拉格朗日粒子，每个粒子均具有相应的质量和体积，通过粒子求和代替核近似中的积分，也称“粒子近似”。如图2所示，基于核近似理论和粒子近似理论，SPH 方法可对控制方程进行离散求解。图中：f和f分别代表标量和矢量； ρ，u，p，r，V分别表示密度、速度、压力、位置矢量和体积，下标i和j表示粒子编号；g，T， µ分别表示重力加速度、应力张量和动力黏性系数；W为核函数；h为光滑长度；c0， ρ0分别为流体声速及其参考密度。为获得更高的压力场、速度场计算精度，可在离散控制方程中进一步添加耗散项[24]。随后，在特定的初始条件和边界条件（自由液面边界[55]、壁面边界[56]、流入流出边界[57]、无反射边界[58]等）下，将离散控制方程在时域内进行积分，可获得SPH 控制方程的数值解。

图2 SPH 方法求解流程[24]Fig. 2 The process of SPH simulation[24]

SPH 粒子携带各自的速度、压力、密度等物理量，按照流体质点本身的速度进行拉格朗日运动。这些粒子的运动代表了流体的流动，粒子所到之处也代表了所处位置的流场速度、压力、密度等物理量的变化。鉴于SPH 理论和方法的拉格朗日特性及无网格粒子特性，在模拟自由液面流动问题时，自由面的运动学条件不仅能够自动满足[55]，而且也能够自动追踪水气和流−固界面，所以十分适合处理船海工程中的大变形、动边界、多介质和强对流的流体动力学问题[8]。

1.1.2 SPH 算法的分类

经过近半个世纪的发展，基于不同的理论框架，SPH 方法已发展出了多种形式，主要包含两种类型：弱可压SPH （weakly compressible SPH,WCSPH）方法[43]和不可压SPH （incompressible SPH, ISPH）方法[59]，并已在各自领域得到广泛的应用。对于水介质，一般认为其近似于不可压缩，但为便于求解压力，Monaghan 等[49]最早引入了状态方程，如图2 所示控制方程中的第4 个等式。基于状态方程求解压力时，积分的时间步长受到克朗条件（Courant–Friedrichs–Lewy condition），也即流体声速的极大限制。此外，Monaghan 等[49]还提出了弱可压假设，以马赫数小于0.1 为标准确定人工声速，确保流体密度变化不大于1%，以此维持流体的近似不可压特性，所以该方法也被称为WCSPH 方法。相反，在ISPH 方法中，流场的密度保持不变，亦即速度散度处处为0。Cummins 和Rudman[59]在SPH 方法框架内，建立了每个时间步的泊松压力方程，用于求解压力，并提出了ISPH 方法来模拟不可压缩流动问题。但是，ISPH 方法涉及大型稀疏矩阵运算和自由液面粒子搜索，计算量相对较大[60-61]，在大规模粒子数量的三维计算中挑战较大。WCSPH 相关理论和方法经过近20 年的发展，计算精度已得到明显提高，相比于ISPH 方法，WCSPH 方法涉及的GPU 存储量更小，计算效率更高，所以更适合于开源程序的设计。例如DualSPHysics[62]、SPHinXsys[63]等开源程序包，它们均是基于WCSPH 理论框架，目前已被广泛应用于模拟三维复杂自由液面流动等问题。

1.2 SPH 数值技术的进展

1.2.1 SPH 压力计算精度

在船海流体力学中，准确计算压力对流体砰击载荷预报至关重要。传统WCSPH 模型计算压力场时，一般会出现难以避免的压力场噪声，其主要原因是：

1） SPH 核函数的近似误差和数值离散误差会导致密度场出现波动，微小的密度波动被状态方程放大，变成剧烈的压力场噪声；

2） 因流体的弱可压假设，声速一般会被降低。在出现流体砰击时，砰击位置辐射的压力波以较小的声速在流场中传播，遇到边界时会发生多次反射和叠加，同样也会在流场中表现为压力场噪声。

针对压力场噪声问题，Colagrossi 等[43]率先提出SPH 密度过滤技术，成功缓解了SPH 压力场脉动。Vila[64]通过Riemann 求解器来计算粒子间的相互作用，提出了Riemann SPH 方法，此方法能够将SPH 数值误差和弱可压假设造成的压力脉动进行快速耗散，从而提高压力场的光滑性。Marrone 等[56]提出δ-SPH 方法并应用于自由液面流动的模拟，该方法引入了密度耗散项Dρ，该耗散项起到与Riemann 求解器相似的耗散作用，也同样能够提高压力场的计算精度。图3 所示的δ-SPH控制方程中，Tv表示黏性力张量。相比Riemann SPH 方法，δ-SPH 方法的优势在于数值耗散量可通过连续方程中的密度耗散项[65]进行定量控制，从而实时监测流体能量的耗散量，这对于SPH 数值波浪水池的建立十分有益。

1.2.2 SPH 粒子均匀化技术

传统SPH 理论基础和算法精度很大程度上依赖于粒子分布的均匀性[66]。为此，Nestor 等[67]最早提出了粒子位移修正方法（particle shifting techniques，PST），而Lind 等[68]将PST 方法引入到自由液面计算中，初步解决了ISPH 模拟中的粒子不均匀及流场撕裂的问题。如图3 中第3 个框图所示，Sun 等[69-70]将PST 技术与δ-SPH 进行结合，提出了δ+-SPH 方法（其中， δu为修正速度项，m表示粒子质量），随后又发展出了张力不稳定性控制（tensile instability control, TIC）技术[44]，克服了高负压区域的流场数值空洞问题，并成功应用到中、高雷诺数的黏性边界层流动的模拟中。其后，又进一步建立了满足一致性的δ+-SPH 方法，将PST 技术中位移修正以修正速度的方式表示，并在控制方程予以考虑，从而解决了早前PST 技术造成的自由液面流动中流体体积的非物理性膨胀问题，极大提高了WCSPH 理论和方法的精度及稳定性[70]。最近，Antuono 和Sun 等[66]又建立了任意拉格朗日−欧拉(arbitrary Lagrangian-Eulerian,ALE)形式的δ-ALE-SPH 理论与方法，将质量方程和质量数值耗散项Dm引入到控制方程中，提高了ALE 算法的稳定性，进一步完善了含粒子位移修正的该SPH 理论框架。针对海洋工程应用问题，Lü和Sun 等[5,17]针对液舱的晃荡、物体出水和入水、黏性绕流等问题进行了深入分析，从理论角度指出了此类问题中张力不稳定性产生的原因，有针对性地提出了完备的控制技术，以克服张力不稳定性问题。

图3 SPH 理论和方法的演化Fig. 3 The evolution of SPH theory and method

1.2.3 高精度SPH 理论和离散格式

制约SPH 方法计算精度的另一个重要原因是，传统SPH 理论和方法离散格式本身的精度不够，误差主要来源于核函数近似误差和粒子近似误差，其中后者又被称为数值离散误差。因此，为提高SPH 方法的计算精度，研究人员先后提出了修正SPH 算法[71]、移动最小二乘算法[72]、有限粒子法[73]以及近两年发展的高阶SPH 方法，例如WENO-SPH 方法[74]、TENO-SPH 方法[75]等。高阶SPH 理论是在Riemann SPH 的基础上发展而来，其通过引入网格法的高阶计算框架来提高求解精度。如图4 所示，在Riemann SPH 控制方程中，e表示比内能[75]，p∗和u∗分别表示粒子对之间的压力及速度[76]，是Riemann 方程求解的中间变量，也是重构框架的作用对象。Meng 和Zhang 等[75]提出的TENO-SPH 理论和方法，有效提高了SPH 方法的计算精度和鲁棒性，验证了高阶SPH 理论和应用的可行性。此外，为减少离散误差，可通过增加光滑长度与粒子间距的比例系数，来增加核函数半径范围内的粒子数量，达到减小离散误差的目的，使之满足理论要求，提高SPH 方法的计算精度，具体详见文献[77]。尽管高精度SPH 方法在解决特定问题方面明显提升了计算精度，但在物理量守恒性及算法稳定性上仍有很大的改善空间。对于工程化应用，因涉及了巨大的粒子数量和较长模拟时间，目前仍普遍采用传统的SPH数值离散理论，辅之以δ-SPH 或Riemann SPH 等数值方法，能够严格保障质量、动量的守恒特性，实现对自由液面流动的长时间的稳定数值模拟。

图4 高阶SPH 理论和方法示意图Fig. 4 Schematics of the high-order SPH theory and method

1.2.4 SPH 方法的工程应用

随着SPH 理论和方法的不断发展，其已被广泛应用于海洋工程领域。事实上，SPH 工程化应用的重要前提是其具备了复杂结构物的粒子化建模和离散能力，以及复杂形状边界条件的处理能力。经过数十年的发展，慕尼黑工业大学Zhu 等[78]开发出了基于.STL 格式文件的粒子化填充算法，并已成功用于复杂构型结构物的粒子化离散。图5 所示为飞机模型从CAD 到粒子化的离散过程，这为后续的工程化应用研究奠定了良好的基础。基于离散后的粒子，将其设置为虚拟粒子（ghost particle），采用基于固定虚粒子（fixed ghost particle）的固壁边界理论和技术[56,79]，即可实现固定、运动甚至变形的壁面边界的高精度施加处理[51, 80]。

图5 复杂构型结构物的CAD 建模和粒子化离散Fig. 5 CAD modeling and particle discretization of a complex structure

1.2.5 SPH 硬件加速技术

在SPH 数值计算中，精确反映三维复杂结构物构型通常需要设置精细的粒子分辨率，故也必然导致高昂的计算成本。随着计算机技术的快速发展，硬件并行计算技术的出现极大地提高了SPH 方法的计算效率，推动了SPH 方法的工程化应用[45,60]。硬件加速技术一般可基于CPU 或GPU两种并行计算框架。绝大部分网格法（例如有限体积法（FVM）等）均采用CPU 并行计算，原因是该方法需要求解大型稀疏矩阵和隐式时间积分，并不适合或较难采用GPU 并行计算。而弱可压SPH 方法作为一种纯拉格朗日粒子方法，其计算原理基于点对点的方式，且时间积分为显式格式，这两个特性决定了SPH 方法非常适合用于GPU并行计算。目前，大量文献已表明，在同等模拟规模下SPH 方法的GPU 并行可比CPU 并行在效率上提高1～2 个数量级[45]，极大提升了SPH 方法的工程实用性，尤其是基于该方法开发的计算软件能够在装有GPU 的个人电脑上运行[62]，并不完全依赖于大型超级计算机或服务器，因此在未来工程领域有较大的应用潜力。

当前，实现SPH 并行计算首先可基于并行计算框架，自主编制SPH 程序，例如Marrone 等[22]搭建的MPI-OpenMP 混合编程模型，设计了一种新型的三维并行SPH 求解器，成功实现了对三维船体兴波的全局仿真，为自主构造SPH 并行计算模型提供了参考。此外，还可以采用目前主流的并行SPH 开源程序包完成计算，例如SPHinXsys[63]、DualSPHysics[62]等开源程序包，它们均在三维复杂问题应用方面有出色表现[81]，为高精度SPH 理论和方法在船海工程中的实际应用提供了强大的计算能力保障。

2 水面航行器高速水动力问题

随着上述数值技术的发展和应用，运用SPH理论和方法求解实际工程问题的能力得到了不断加强。SPH 方法利用其在处理自由液面大变形和喷溅问题方面的优势，可以十分方便地对高速水动力问题开展相关研究。例如，图6 展示的是三维SPH 方法模拟的“中山大学”号（SYSU）科考船高速航行兴波和飞机水上迫降过程。由图可见，SPH 方法准确再现了自由液面的翻卷、破碎、喷溅过程。

图6 三维SPH 方法模拟的“中山大学”号科考船高速航行兴波和飞机水上迫降Fig. 6 3D SPH simulations of the breaking wave of SYSU research vessel in high-speed and the aircraft ditching process

以下将围绕航行器的兴波与非线性波浪砰击等问题，对SPH 理论和方法在水面航行器高速水动力问题中的应用研究进展进行综述。

2.1 高速水面航行器兴波

目前，计算流−固耦合问题通常使用有限元法（FEM）、FVM 等工业化程度高、应用发展较为成熟的求解器，例如，常使用VOF[82]或Level-Set[83]等方法模拟自由面流动。在单相、小变形、低弗劳德数的情况下，使用以上的网格法模拟航行器的兴波和尾流可以得到精准结果，但对高速航行过程中产生的强非线性自由液面流动、流体喷溅等情况，网格法较难精准捕捉到波浪破碎、浪花飞溅等细节。图7 展示了FVM 和SPH 方法模拟高速舰船兴波的结果与Tagliafierro 等[84]的模型试验结果的对比。由图可见，FVM 只能够模拟出微幅的液面变形，而SPH 方法能够很好地捕捉液面的飞溅细节。可见，SPH 充分发挥了拉格朗日粒子方法的理论优势，能更准确模拟兴波特征。

图7 高速滑行艇兴波各方法模拟结果对比Fig. 7 Comparison of wave-making of high-speed planning boat

早期的计算机内存和效率低，三维仿真计算成本高，故一般采用切片法、二维半[85-86]理论研究航行器的兴波问题。Marrone 等[4]通过2D+tSPH方法对细长船体兴波的问题进行研究，对不同弗劳德数下的船艏兴波演化和尾迹进行模拟，较成功地模拟了细长船舶两侧兴波特征，并通过试验对比验证了结果的准确度。但是，2D+tSPH 方法的研究对象局限于纵向速度变化小的细长船体，对于横剖面变化大的复杂构型，只有进行三维仿真才可以尽可能全面地保留兴波特征。最近，Mintu等[87]成功实现了实尺度高速舰船船艏破波的三维模拟，验证了三维SPH 计算方法的数值精度和技术可行性。

与真实试验结果[84]相比，SPH 方法已可以初步模拟船体两侧喷射的水花，如图7 所示，但目前仍未能对微小尺度的水气混合多相流实现精准捕捉。高密度比混合多相流模拟是未来SPH 方法的发展方向，需选用更精细的粒子分辨率、开发新的理论和数值模型进行求解。Monaghan 和Kocharyan[88]首次提出了基于体积分数理论的SPH混合物模型，经过十余年的发展，Fonty等[89-90]提出的SPH 混合物模型在模拟振荡水柱波浪能转换器的水−气高速耦合作用方面表现优越，可进一步拓展应用于航行器（体）的兴波模拟中。SPH混合物模型在高速水动力学应用领域发展前景广阔，相信未来会应用于模拟兴波水花及其他大密度比多相流体扩散、混合过程。

2.2 非线性波浪砰击

2.2.1 SPH 数值造波与消波

水面航行器高速航行时会受到反复交替的非线性波浪砰击载荷作用，尤其是在恶劣海况下，波浪力所引起的砰击压力和力矩可能会造成海洋结构物的局部毁伤甚至整体倾覆，因此研究波浪与浮体耦合作用意义重大。近年来，数值波浪水池以及波浪与结构物的耦合作用模型迅速发展，其中数值造波与消波技术成为研究重点。得益于SPH 方法模拟动边界的天然优势，SPH 方法可通过固定虚粒子理论和技术构建数值造波板，对其运动予以控制，模拟物理水池中的推板式[91]、活塞式[92]、摇板式[93]等造波机，生成多种形式的波浪[94]。

SPH 数值水池消波技术可分为被动消波和主动消波。被动消波是指通过斜坡和多孔结构[95]或在数值水池中设置人工阻尼区进行消波[92,96]，而主动消波方法是通过改变数值造波板的运动机制达到消波目的。虚粒子组成的造波板会根据反射波的变化重新调整运动方程，使反射波与新生成的波叠加形成与原波形一致的波浪[91]，从而保证结构物受恒定的数值波浪作用。研究表明，SPH 方法可以产生包括不规则波[91]、畸形波[94]、孤立波[92]等各类波形，造波质量较高，为进一步研究浪固耦合问题提供了技术支撑和保障。

2.2.2 甲板上浪与畸形波砰击

若巨幅波浪砰击船体、海洋平台等结构物，一般会造成严重的波浪爬升，甚至可能翻越甲板，发生甲板上浪现象，导致甲板设备严重毁坏或人员伤亡等，因此，受到了研究人员的高度重视[97-98]。早期的SPH 方法在研究舰船甲板上浪问题时，为降低计算难度，通常是通过二维仿真波浪冲击翻越固定水平平板[10]过程来进行研究。而随着计算机计算能力的提高，现已能够对舰船航行的上浪情况进行三维模拟，进一步研究甲板上浪对上层建筑的破坏力[99]。为更精确预报砰击载荷，Areu-Rangel 等[100]通过对比SPH 方法的仿真结果和精密实验仪器的测量结果，优化了甲板上浪的SPH 数值模型，该方法用于预报甲板上浪压力载荷的精度得以验证，有效性和可行性也得到证明。

传统的规则波和不规则波理论难以解释深海中时常出现的极大波高，因此，近年来兴起了对极端海况下畸形波和畸形波砰击海洋结构物的研究。Sun 等[94]基于规则波和畸形波理论，采用推板和摇板式造波数值技术联合人工阻尼消波技术，建立了数值波浪水池的SPH 数值模型，以研究规则波、聚焦波与结构物的浪固耦合作用，结果如图8（a）所示。针对Zhao 等[101-102]提出的畸形波−结构物耦合运动算例，SPH 方法可预报浮体结构物的纵荡、升沉、纵摇运动和波浪砰击、甲板上浪、回流等现象，计算结果与试验结果吻合，如图8（b）所示。

图8 SPH 数值波浪水槽模拟的非线性波浪砰击固定或浮动结构物的结果对比Fig. 8 The results comparison of nonlinear wave slamming on a fixed or floating structure simulated in a SPH numerical wave flume

畸形波砰击过程中会卷入气体，波浪翻卷上涌甲板后，受气相的影响出现负压区，如图8（a）所示。Sun 等[94]通过研究砰击过程中的负压现象，指出一次砰击过程相当于进水与出水过程的耦合，非线性波浪砰击退去后，会在结构下方（例如外漂底部）产生负压，而负压在SPH 方法中会带来张力不稳定性，随之造成非物理性的流场空洞并导致流场演化的错误。文献[94] 模拟结果还表明，当波浪砰击甲板边缘时，由于砰击速度高，卷入的气穴会形成负压（相对于背景压力）而导致流体吸附。此时，流体并不会立即在甲板上分离，而是在甲板表面形成高速吸附流，这就要求SPH 模拟能精确地预报流场中负压形成的吸附力。对于此问题的研究，研究人员建议使用SPH 张力不稳定性控制技术[44]。此外，波浪翻越甲板的过程还需考虑空气的气垫效应[94]，否则容易导致压力峰值预报失真。总之，涉及甲板上浪的浪固耦合问题还需进一步完善SPH 的理论和方法，建立多相流SPH 模型[17]，在考虑空气影响的情况下进行数值模拟研究。

3 航行体跨介质水动力问题

航行体跨介质过程是指物体出水和入水，实质上也是一个复杂的多相流−固耦合过程。流体与结构物的相互强耦合可能导致相变、湍流和空化，力学机理极具复杂性和挑战性。目前，对出水和入水力学机制的研究方法主要采用试验方法[103]、数值模拟[51]或两者结合的方法。研究表明，SPH 方法能很好地模拟此类大变形、多边界和多相流问题[15]。

3.1 细长体入水

细长航行体入水是典型的瞬态强非线性动力学过程，入水瞬间的砰击作用会造成强烈的液面变形，形成液滴喷溅。面对此类问题，网格法受网格尺度限制，对喷溅的模拟略显逊色。图9 给出了采用SPH 方法、VOF 方法计算模拟细长圆柱体（细长体）以100 m/s 高速入水的结果对比。由图9可见，网格法（VOF）难以捕捉到喷溅细节，多相界面相对模糊，而SPH 方法在流场喷溅特性的捕捉方面表现更佳。图10 所示为该算例中运用不同算法预报的细长体端部砰击压力和运动速度曲线。由图10 可见，SPH 方法能够精确地捕捉到瞬态冲击载荷的峰值，获得的细长体运动速度与CEL 方法更吻合。可见，SPH 方法能够精确地对航行体高速入水冲击进行仿真，尤其是液面喷溅模拟方面的优势更明显。

图9 不同方法模拟的细长体高速入水结果对比 [50]Fig. 9 Comparison of high-speed water entry of a slender body simulated by SPH and VOF [50]

图10 细长体高速入水砰击压力与运动速度变化曲线Fig. 10 Time evolution of impact pressure and velocity of the slender body during high-speed water entry

高速入水问题主要关注细长体结构局部砰击区域所受载荷及其附近流场的演化过程，而远场流动变化随着距离的增大而减小。因此，为了在保证局部高精度的同时提高计算效率，Sun 等[44]采用自适应粒子细化技术（adaptive particle refinement, APR）对大范围的计算域分设多级粒子分辨率。该技术也能根据需求定义高分辨率区域，通过此区域动态跟踪细长体结构运动的过程，并可对整个弹体周围流场区域进行加密，精准预报流体对弹体和弹道的影响，如图11（a）所示。此外，APR 技术也可仅对弹体头部进行局部加密，捕捉弹头区域的流动特征，如图11（b）所示（ ∆x表示粒子间距）。总之，APR 技术为SPH 方法实现高效精准的仿真提供了技术支撑。

图11 多级粒子分辨率SPH 技术模拟的细长体高速入水过程Fig. 11 High-speed water entry process of a slender body simulated by SPH with APR technique

前人基于SPH 方法的许多入水问题的研究，例如旋转[104]、倾斜[105]入水的仿真以及空腔演化特征等，大部分忽略了空腔闭合后的内部气体对结构运动的影响。首先，这是因为空腔闭合后会产生高速射流[17]，冲击结构物表面并造成运动姿态变化，所以空腔的产生及气相的影响对结构载荷与运动预报有着重要意义[103,106]，需要采用多相流SPH 模型开展计算，详见文献[17]。其次，砰击入水的高速运动会导致细长体周边产生空化，且目前SPH 理论框架内尚未建立完备的空化理论和数值模型，如何做到空化、湍流[107]、多相SPH 模型的结合以实现精准数值预报仍有待深入研究。

3.2 细长体出水

SPH 方法为细长体高速入水过程的模拟提供了可靠数值方案，对于细长体出水过程其同样具有较精确的预报能力。此外，高速出水过程还涉及复杂气−液−固三相的耦合，包括瞬态强冲击载荷、空泡产生及演化和流−固耦合效应等，这也是工程前沿所要研究的重要课题。为此，前人开展了丰富的数值和试验研究，包括球体出水[109]和圆柱出水仿真并结合试验校核[110]等。

在细长体高速出水时，其肩部或底部会出现相对背景压力的负压现象，负压在SPH 数值计算中会导致强烈的张力不稳定。通过建立张力不稳定性控制理论（tensile instability control, TIC）可避免细长体出水过程的非物理性数值空洞[5]，而且，还可通过多级粒子分辨率技术对细长体附近区域的粒子进行加密，提高局部区域的精度。 图12 给出了细长体高速出水过程的初步模拟结果。由图可见，对于高速出水问题研究，SPH 方法已具备较强的预报能力。

图12 SPH 方法模拟的细长体高速出水过程Fig. 12 Numerical simulation of high-speed water exit process of a slender body by SPH method

细长体从空中入水的瞬态过程中砰击载荷占主导地位，计算过程一般忽略黏性力，但因出水过程涉及了空化、湍流、强可压过程，通常受黏性严重影响，所以出水与入水并非是完全逆向过程[8]。Zhang 和Sun 等[8]采用δ+-SPH 方法研究二维圆柱体出水过程，使用TIC 技术成功模拟了圆柱体底部产生的涡旋分离，所得结果比边界元法（BEM）更精准。

高速出水问题还需要考虑细长体出水过程肩部产生的空泡，这就需要开发SPH 理论框架下的空化模型，精确预报空泡的产生、演化和溃灭过程。空化问题具有强可压特性，Sun 等[111]提出的体积自适应粒子撕裂与融合理论和技术在应对强可压的问题时具有较强的适用性，能够为开发SPH 空化模型提供基础支撑。然而，目前传统的计算流体力学（CFD）算法模拟空化相变的相关技术[112]暂时难以直接拓展应用到SPH 算法中，亟需在未来研究中重点开展此项工作。

4 水下爆炸与结构毁伤

水下爆炸是一个涉及高温高压爆炸产物剧烈膨胀的强非线性过程[19]，爆炸过程分别与船体结构和海况相关，研究极具挑战性，难以通过单纯的理论和试验方法开展研究，而数值模拟是实现实尺度预报的最佳策略之一。

从目前来看，远场水下爆炸载荷的预报理论相对完善，但近场和超近场水下爆炸载荷预报仍十分困难。SPH 方法模拟水下爆炸问题时，可自动捕捉多相界面，模拟从起爆、冲击波传播到气泡膨胀的全过程。当SPH 算法与结构求解器进行耦合时，可计算船体从损毁到沉没的全过程[6,113]。本节将对SPH 理论和方法在研究水下爆炸冲击波和气泡脉动方面取得的进展进行综述。

4.1 水下爆炸冲击波模拟

炸药在水下引爆后会产生高温高压的爆炸产物，辐射的高热量使附近的水瞬时汽化，形成高温高压气泡，并向外辐射形成冲击波。水下爆炸中冲击波的威力占据了总能量的一半左右。爆炸后产生的瞬态冲击波在多相中的传播，具有强可压、强非线性等特征，这对SPH 算法的精确性和稳定性要求极高。而以往的SPH 理论和方法基于不可压缩[59]或弱可压假设[49]，难以直接应用于水下爆炸这类强压缩性问题，尤其是对水下爆炸气泡的周期性膨胀和收缩过程难以精准模拟。

Riemann SPH 方法的出现推动了粒子方法处理可压缩问题的应用发展，但早期在捕捉激波时的精度不如成熟的网格法，在理论方面还需进一步提高。因此，基于传统SPH 理论，借鉴和引入网格法中处理可压缩流动的先进数值技术及高精度计算框架，Meng 和Zhang 等[75]提出的TENOSPH 模型与方法，在冲击波捕捉和分辨方面具备了较高的数值精度和计算能力，增强了Riemann SPH 方法处理此类问题的鲁棒性。Sun 等[111]针对强可压问题的数值特征，提出一种新型的粒子体积自适应理论，实现了粒子在膨胀作用下撕裂和收缩时的融合模拟，较准确地预报了冲击波载荷。但是，对于气体膨胀压缩导致的核函数影响域内粒子数的变化问题，通常还需采用变光滑长度技术[114]，来提高SPH 方法的计算精度和稳定性。

在水下爆炸的具体计算中，SPH 相关处理技术和方法取得了可观的进展。Wang 和Zhang 等[58]提出基于特征线原理的无反射边界条件，有效消除了水下爆炸冲击波计算中的边界反射，为SPH方法用于高精度模拟水下爆炸提供了依据。彭玉祥[115]采用SPH 方法模拟文献[116]中相同工况下的二维近平面固壁柱状水下爆炸算例，研究了高压气泡产生的冲击波和稀疏波的传播，求解的压力结果与参考结果高度吻合，为SPH 方法用于模拟强可压瞬态流动的可行性和有效性提供了参考。同时，文献[115]针对三维自由场的水下爆炸问题进行模拟，将计算得到的冲击波压力时历曲线与经验公式进行了对比，两者吻合较好，体现了SPH 方法预报冲击波压力的良好精度。

在计算三维轴对称水下爆炸问题时，还可以通过轴对称SPH 数值模型[117-119]进行模拟，以减小计算量。图13 给出了采用轴对称SPH 方法计算得到的水下爆炸几个典型时刻的冲击波压力云图。轴对称SPH 方法的计算量小，相对于直接三维模拟，可采用更加精细的粒子分辨率，因而计算得到的冲击波波阵面十分锐利。由图可见，SPH 方法模拟水下爆炸冲击波阶段已发展得相当成熟可靠，具备了水下爆炸冲击波阶段载荷的精确预报能力。

图13 轴对称SPH 方法模拟的水下爆炸几个典型时刻的冲击波压力云图Fig. 13 Pressure contours of the shock wave induced by underwater explosion at typical time instant by axisymmetric SPH model

4.2 水下爆炸气泡模拟

May 和Monaghan[120]最早采用SPH 方法研究了气泡引起的舰船沉没过程。在此研究领域，文献[10]主要考虑的是巨大的上浮气泡，而事实上更不可忽视的是水下爆炸气泡及其在舰船附近脉动激起的鞭状运动和射流对结构的冲击效应。研究水下爆炸产生的大尺度脉动气泡，需考虑环境中的静水压力、大气压力、重力及边界等的耦合作用，以及气泡脉动产生的热交换和能量耗散，此机理十分复杂。在上述方面，前人提出了理论、试验和计算等多种方法进行研究，取得了丰硕成果[7,121]。

SPH 方法尽管在处理大变形和动边界问题时具有明显优势，但用于模拟水下爆炸气泡脉动较为少见，其主要原因是，爆炸气泡涉及了爆炸气体的强压缩和膨胀过程，给粒子法模拟带来了挑战。Joshi 等[122]对固体边界的空泡崩塌进行了模拟，然而忽略了气体粒子以及气泡的膨胀过程；Pineda 等[123]也进行了多相模型模拟空泡崩塌，但因缺乏粒子体积自适应技术，气体粒子体积被过度压缩，导致水−气界面出现了明显的锯齿状。总之，水下爆炸气泡SPH 模拟的难点主要体现在气泡脉动过程中。SPH 理论的拉格朗日粒子特性决定了气体粒子分布会随气泡体积大幅度变化。为简化问题的复杂性，目前，SPH 方法暂未计入气泡脉动过程中的热力学效应和能量损失，首要任务是需要提出、验证和优化适用于强可压问题的SPH 理论，以增强SPH 方法模拟气泡脉动的可行性，提高对气泡脉动和射流的仿真精度。

采用SPH 方法计算爆炸气泡时，粒子不仅会因气泡膨胀变得十分稀疏，也会因气泡收缩而变得十分稠密，使得计算的载荷结果不准确[111]。如图14 所示，当气泡体积收缩至最小时（见图14 左侧），气体粒子严重集聚，核函数半径范围内粒子数量急剧增多；相反，当气泡体积膨胀到最大时，气泡粒子的间距随之增大，核函数半径范围内粒子数量急剧减少（见图14 右侧）。以上气泡收缩和膨胀过程中，核函数内部粒子数量大幅变化，但气泡外围的水粒子间距几乎不变（水的压缩性较小所导致）。由此，造成水−气界面的两侧粒子分布极不均衡，使得水−气界面不光顺，计算精度剧烈下降甚至计算中断。为解决SPH 方法模拟空泡脉动时气体粒子的体积和间距剧烈变化所带来的计算精度及稳定性问题，Sun 等[111]建立了基于SPH 体积自适应理论（volume adaptive scheme,VAS）的粒子撕裂−融合技术，成功模拟了水下爆炸气泡在自由场、自由液面、壁面等边界下的脉动过程。

图14 气泡脉动时气体粒子体积和间距变化示意图Fig. 14 Variation of the volume and spacing of gas particles during the pulsation of a bubble

尽管SPH 理论和方法拥有精确捕捉大变形水−气界面的优势，但在处理三维气泡脉动问题时，因需要较大的计算域以避免边界影响而使得计算量较大、效率较低。为此，Sun 等[117]构建了轴对称强可压多相流SPH 模型，实现了SPH 方法对水下爆炸气泡脉动的模拟，可精确捕捉高压气泡的膨胀、收缩、崩塌、射流等过程。图15 给出了自由场溃灭气泡半径演化和自由液面附近气泡膨胀及收缩时最大宽度演化的SPH 模拟结果与参考数据[124-126]的对比验证[117]。由图可见，轴对称SPH方法能准确预报不同条件下气泡的膨胀和收缩动态过程。

图15 水下爆炸气泡脉动的模拟SPH 模拟与验证[117]Fig. 15 SPH simulation and validation of bubble pulsation in underwater explosion [117]

本文采用Sun 等[117]的轴对称SPH 方法，模拟了Tian 等[127]给出的自由液面以下10 m 处500 kg装药所产生的水下爆炸气泡脉动和射流过程，如图16 所示。结果表明，计算得到的气泡膨胀和射流过程与Tian 等[127]提供的数值结果十分吻合。从数值模拟结果可以看出，基于Sun 等[111]提出的VAS 理论，采用粒子撕裂−融合技术，SPH 算法可用于模拟水下爆炸冲击波和气泡脉动的全过程。最近，Fang 等[119]将VAS 理论与Riemann SPH 方法结合应用于激波及入水问题，也取得了较好的计算精度，验证了SPH 理论和方法应用于水下爆炸及气泡动力学的可行性及有效性，为进一步深入研究水下爆炸力学和载荷机理提供了更多元的求解方案。

图16 自由液面附近水下爆炸气泡脉动和射流的SPH 模拟Fig. 16 SPH simulations of bubble pulsation and jetting near a free surface

4.3 水下爆炸结构毁伤模拟

如上所述，SPH 方法因其不依赖网格，极为适合用于模拟大变形流−固耦合问题。早在1991 年，Libersky 等[128]就将SPH 方法应用到了固体结构的冲击碰撞断裂的研究中，将SPH 方法的应用领域拓宽到了固体力学领域。早期运用SPH 方法研究固体断裂时都采用三维实体对结构进行建模，但考虑到汽车、飞机和舰船等工程结构中广泛采用的是薄壳结构，采用三维实体建模会导致壳体厚度方向的粒子间距极小，从而限制时间积分的增量，使得计算效率降低。在2008 年，Maurel 等[129]提出了SPH 壳单元理论，解决了时间增量受壳体厚度限制的缺点，极大提高了SPH 方法在薄壳结构动响应方面的求解能力。此后，Caleyron 等[130]又将该理论应用到了对薄壳结构冲击断裂响应的模拟中。同时，Ming 等[131]应用并发展了该理论，将该理论命名为光滑粒子壳单元（smoothed particle shell，SPS），并应用到了水下爆炸结构毁伤的模拟中[132]。

然而，因SPH 理论和方法采用的是强形式的控制方程，经过使用移动最小二乘法（moving least squares， MLS）形函数修正其核函数后，使之不能自动满足结构自由边界条件，需要人工施加自由边界条件[133-134]，增加了断裂自由边模拟的复杂度。针对此问题，Peng 等[135-136]采用重构核粒子方法（reproducing kernel particle method，RKPM）建立了任意曲壳曲梁三维无网格法的计算模型，实现了加筋壳结构的无网格模拟。而且，在上述基础上，建立了水下爆炸结构毁伤的流−固耦合计算模型[113]，分别采用SPH 和RKPM 方法求解水下爆炸载荷和舰船结构的断裂响应，并与试验结果进行对比，验证了所提计算模型的有效性和计算精度。图17 给出了水下爆炸时舰船结构断裂毁伤的SPH-RKPM 数值仿真结果。由图可见，该数值计算模型具备了预报水下爆炸时的舰船毁伤能力。目前，张阿漫教授团队建立了SPH-RKPM 流−固耦合动力学模型及计算方法，开发了具有自主知识产权的水下爆炸舰船毁伤数值计算FSLAB软件[6]，该软件被应用到了舰船结构强冲击毁伤的模拟之中，广泛论证了SPH-RKPM 耦合计算方法在水下爆炸领域的应用前景，能够为舰船防护结构设计提供理论依据和基础性技术支撑。

图17 SPH-RKPM 算法模拟水下爆炸引起的舰船结构毁伤过程Fig. 17 Numerical simulations of structural damage resulted from underwater explosion simulated by SPH-RKPM method

4.4 破舱涌流与结构动态沉没模拟

SPH 理论和方法除了能够精确预报船体在冲击载荷作用下的结构破损外，对破舱涌流和结构沉没的模拟也具有强适用性。在船体结构破损情况下，涌流灌入破口后与波浪的相互耦合作用，会使船体产生剧烈的横摇、升沉等多自由度的复杂运动。Guo 等[137]采用三维SPH 方法对静水中不同破口位置的破舱动态沉没过程进行了仿真，通过相关试验验证了SPH 方法的精度和稳定性，为破舱动态沉没问题提供了基准算例，其舷侧破口涌流的仿真结果与试验结果对比如图18（a）所示。 Cheng 等[138]、Ming 等[139]和Cao 等[140]构造了SPH 数值水池，对多水密舱体在波浪耦合作用下的沉没过程进行研究，准确模拟了舱内涌流晃动和外部波浪的耦合过程，进一步完善了破损结构在波浪环境下的动态沉没力学机理，如图18（b）所示的沉没过程[138]。综上，建议未来考虑舱内气体的真实压缩性，基于VAS 框架的多级粒子分辨率技术和多相流模型[141]，进一步开发适用于破舱涌流的高精度SPH 理论和方法，以提高模拟精度和稳定性。SPH 方法用于水下爆炸的各环节及船体破口与动态沉没过程均具有较强的数值模拟能力，预期未来能够实现将SPH 方法应用到水下爆炸全过程的精准仿真。

图18 三维SPH 方法模拟破损舱段涌流和沉没过程的结果与试验结果对比Fig. 18 Comparison between experimental and 3D SPH simulation results of flowing and sinking process due to a damaged ship cabin

5 高速水动力学耦合算法

考虑到SPH 数值模拟在实际工程应用中难以规避计算成本的问题，为充分发挥SPH 模拟大变形问题的优势，提高其在船海力学领域的应用范围，建议SPH 方法在未来发展上应重视与其他算法的耦合。例如，开发有/无网格耦合算法、黏流−势流耦合算法等，构建如图19 所示的耦合计算框架，以解决船海高速水动力领域常见的出水和入水以及水下爆炸与冲击问题。图19 所示的耦合算法可在流−固耦合作用剧烈的区域采用SPH 粒子法进行局部范围的精准计算，在较小变形的外围区域使用黏流网格法进行高效求解，在远场区域可采用势流方法进行大范围海域的模拟。事实上，文献[142] 中的部分耦合算法已为该框架的搭建和实现奠定了坚实的技术基础，例如，SPH-FVM 算法[143]在大变形流动仿真上拥有较好的稳定性，已拓展应用至航行体的兴波、高速入水等问题。此外，高阶谱（higher-order spectral, HOS）方法与网格法求解器的耦合也得到了验证[144]，其适合用于对大范围海域内的船海水动力问题进行黏流−势流耦合计算。图19 所示的多算法耦合计算框架通过结合有网格、无网格和黏流、势流数值算法的优点进行耦合求解，有望为SPH 理论和方法应用于实际船海工程问题提供新的技术方案，为实现高精度、高效率的船海工程应用仿真奠定基础。

图19 有/无网格耦合算法和黏流−势流耦合算法求解典型高速水动力问题示意图Fig. 19 Schematics of coupled algorithm between mesh-based and particle-based methods and between viscous-flow and potential-flow solvers for typical high-speed hydrodynamic problems

6 结 语

本文综述了SPH 理论与方法在船海高速水动力学研究中的最新进展。可见，SPH 理论和方法日趋得到完善，相关数值方法在高速舰船兴波、细长体出水和入水、水下爆炸与冲击等高速水动力问题中已得到广泛应用，取得了较高的计算精度，在精确捕捉大变形的多相界面和流−固界面等方面具备明显优势。

在未来研究中，SPH 理论与方法在以下几个方面仍有提高的空间，主要包括：

1） 在计算精度和效率方面，进一步提高SPH 理论和方法的完备性及数值格式的收敛阶数，实现以较少粒子数量达到预期的计算精度，提高计算效率。

2） 在边界处理方面，仍需进一步发展三维复杂构型壁面边界的SPH 数值处理技术，减弱边界效应对物理问题模拟结果的影响。

3） 在应用范围和领域方面，开发SPH 理论框架下的湍流模型和空化模型，加强高雷诺数流动和相变问题的模拟能力。同时，发展基于体积分数的SPH 水气混合物的理论和数值模型，进一步拓展应用于高速航行体的兴波和水汽泡混合流模拟中。

4） 在解决实际工程问题方面，进一步提高SPH 方法和计算程序的高性能并行计算能力，与国产大型超级计算机的能力结合，实现亿级粒子数的大规模并行化模拟，增强SPH 理论和方法的工程化应用能力。

7 致 谢

感谢哈尔滨工程大学明付仁副教授、王平平博士和中国海洋大学程晗博士对本文水下爆炸冲击波和细长体出入水计算提供的帮助。