SPH方法模拟水下爆炸研究进展

杨文山，孟晓宇，王祖华

(武汉第二船舶设计研究所，湖北武汉430064)

0 引言

随着SPH(Smoothed Particle Hydrodynamics)方法在具有材料强度的动态响应问题和具有大变形流体动力学问题上的成功应用，许多学者开始考虑用其模拟水下爆炸的可行性。Swegle［1］最早对SPH方法模拟水下爆炸的可行性进行了论证，发现SPH的无网格性质适合于解决水下爆炸产生的大变形问题，而其拉格朗日性质极易捕捉爆炸运动物质交界面和自由面，在水下爆炸的模拟上具有先天优势，非常适合水下爆炸冲击波阶段的模拟。Liu［2］在SPH方法模拟水下爆炸领域作出了突出贡献，对一维炸药爆轰、二维爆炸气体膨胀、锥形炸药爆炸、二维水下爆炸、水介质缓冲等问题进行了模拟。然而其模拟均是针对自由场爆炸等比较理想的模型，仅能对算法的稳定性、精度、收敛性等进行验证，且仅局限于二维水下爆炸的模拟，无法实现工程应用。

国内杨刚［3］、姚熊亮［4］分别对近自由面水下爆炸和沉底水下爆炸进行了模拟，Xu［5］采用SPH和有限元结合的方法模拟了结构在水下爆炸作用下的响应。以上水下爆炸的模拟均是在一维或二维自由场条件下进行的，且仅涉及到炸药和水等密度差别较小的物质，而工程中的水下爆炸均属于三维问题，存在着自由面、弹性边界及无反射边界，且多为水、空气、爆炸气体、钢等多种密度差异较大物质混合的多相流问题。以上原因使得SPH方法难以解决许多工程问题，例如:舰船接触爆炸、舰船防雷舱爆炸的模拟对舰船防护结构设计具有重要意义，但其涉及TNT装药、空气、水、钢等性质差异较大物质的相互作用，且钢板厚度远小于水域和TNT装药的尺度，使其粒子难以均匀分布，材料和粒子的高度非均匀性使SPH方法模拟接触爆炸存在较大困难;工程中常用的柱形和方形TNT装药爆炸均属于三维问题，且需考虑自由面和无反射边界的影响，而当前的研究仅局限于二维自由场水下爆炸的模拟，难以解决实际工程问题。因此，大密度比多相流模拟、三维或基于对称算法的三维模拟、无反射边界等边界的处理是SPH方法实现水下爆炸工程应用的主要技术瓶颈，需要进行深入研究。

1 SPH模拟多相流问题

由于SPH方法是拉格朗日性质和粒子性质的和谐结合，容易追踪界面，因此特别适合模拟多相流等异种介质交界面不断运动的问题。Monaghan［6］，Cleary［7］和 Ritchie［8］分别对火山喷发气尘两相运动、多物质热传导、多物质星云形成等多相流问题进行了模拟，并发现密度差别较小的多相流问题可认为其密度梯度小于光滑函数梯度，采用标准SPH方法进行模拟便可达到精度要求。然而，当多相流中的物质密度差异较大时，交界面的计算精度则得不到保证。这是由于当高密度介质和低密度介质互为交界面时，低密度介质会因较大的密度梯度使近似结果偏大，而高密度介质的近似值比实际值偏小，这种误差的累积会使计算值失真，严重时会导致计算程序的发散。

许多学者对大密度比的多相流问题进行了研究，Valizadeh［9］对交界面处的粒子质量进行修正，模拟大密度比的多相流问题，但其修正系数的取值需随密度比的变化而变化，因此该方法通用性较差。Colagrossi［10］通过改进SPH的近似方式，并周期性地初始化密度获得了良好的模拟效果。Hu［11］将体积近似代替密度近似应用于连续密度方程和动量方程，来求解大密度比的多相流问题。Solenthaler［12］采用“粒子密度”近似密度、压力和粘性力亦得到了良好的效果。韩旭［13］对密度近似进行了正则化，并对气体状态方程加入了范德华修正项成功模拟了气泡在水中的上浮运动。以上研究均是对气泡上浮、溃坝等低速流问题进行模拟，而水下爆炸为大密度比高速多相流问题，物理模型中空气和钢的密度比达1/7800，且爆炸产生的物质运动速度可达km/s，和低速流有较大差异，因此大密度比高速多相流的模拟是实现水下爆炸工程应用的重要技术。

2 SPH模拟三维问题

SPH方法的一大特点是其程序结构简单，可自然扩展到三维，但用SPH方法模拟三维问题需要克服其计算量较大的缺点。影响SPH方法计算量的主要是粒子搜索，这是由于SPH方法是采用i粒子周边支持域内的粒子对i粒子进行近似，且整个问题域中的粒子是不断运动的，因此在每个时间步内均需要对粒子进行搜索以确定i粒子对应的近似粒子，使得粒子搜索成为影响SPH方法计算量的决定因素。在SPH研究的初期采用全配对搜索法搜索粒子，该算法的形式最为简单，但其计算量为O(N2)(N为问题域的总粒子数)，在计算粒子数较多的问题时速度较慢。为提高粒子的搜索效率，Monaghan［6］提出了链表搜索法，该方法将粒子分布在胞元内，通过胞元的存储记录搜索粒子，理想状态下该方法的计算量为O(N)，但该方法对于问题域尺度在计算过程中大幅扩张、光滑长度发生巨大变化的情况搜索效率较低。Hernquist［14］提出了适合求解变光滑长度的树形搜索法，该方法将问题域分成多个卦限，直至每个卦限仅包含1个粒子，通过有序的树形卦限对粒子进行搜索，其计算量为O(N1gN)，树形搜索算法的应用在保证精度的前提下使计算效率大幅提高，采用该算法已可对数百万粒子的问题进行模拟。粒子搜索算法计算效率的提高使得SPH模拟三维问题成为可能。Cleary［15］采用三维SPH方法模拟了轻金属的高压压铸过程;Dalrymple［16］模拟了溃坝水流对结构的作用;Lönner［17］对 三 维 自 由 表 面 流 问 题 进 行 了 模 拟;Vuyst［18］采用有限元和SPH结合的方法模拟了三维高速冲击问题;国内杨秀峰［19］等人也对三维溃坝问题进行了模拟。

模拟三维问题需要的计算量较大，当三维问题具有对称特征时，可采用一维球对称、二维轴对称或三维平面对称算法进行模拟，从而使计算量大幅降低。对称SPH算法中用1个粒子代表周向所有粒子，故在对称轴处粒子质量分布极不均匀，导致对称轴处的计算精度较低，因此其研究远没有标准SPH算法成熟。Petschek 等对三维笛卡儿坐标系中的角坐标积分，将三维坐标系转化为二维轴对称柱坐标系，但其计算结果在对称轴处存在问题。Brookshaw［21］直接应用二维轴对称柱坐标系的控制方程将三维问题转化为二维问题，但该方法的精度在对称轴附近得不到保证。Omang［22］在不同的坐标系下推导了对称问题的核函数，然后采用拉格朗日公式离散了对称问题的控制方程，其计算结果在对称轴处仍具有较高的精度。国内龚凯［23］、杨刚［24］等人也开发了对称SPH方法计算程序，但均没有讨论对称轴处计算结果的精度。从以上研究可知，基于对称算法的三维SPH方法仍不成熟，且三维及基于对称算法的三维水下爆炸模拟的相关研究极少。

3 SPH处理边界问题

边界条件难以施加、对边界处的粒子近似时存在缺陷是影响SPH方法计算精度的主要因素，这使得学术界一直致力于施加边界条件和修正边界缺陷的研究。Monaghan［6］处理固壁边界时在固壁处施加一排虚粒子，通过虚粒子的排斥力防止粒子穿透固壁。Libersky以自由面为对称面反射对称虚粒子，通过虚粒子的近似修正自由面边界处因缺少近似粒子导致的缺陷。Randles将虚粒子的所有参变量赋予和边界处粒子等量的值，通过虚粒子和边界内粒子插值得到内部粒子的变量值。Liu［2］在Monaghan和Libersky的基础上采用2种虚粒子处理固壁，第一种虚粒子布置在固壁边界上对内部粒子施加排斥力，第二种虚粒子分布在边界外参与粒子近似，该方法防止了粒子穿透边界，且提高了边界处的近似精度。以上方法均采用了虚粒子，可以解决近似边界附近粒子时由于粒子不充足而引起的边界缺陷问题。

导致SPH方法边界缺陷的另一个原因是核近似导函数时在边界处会发生截断，使分部积分的面积分项不为0。为此，Campbell［25］将边界余项引入到分部积分中，修正了边界的截断问题。此外，一些SPH方法的改进形式对修正边界缺陷、提高计算精度起到了良好效果，如Chen［26］根据泰勒展开构造的正则化形式的修正光滑粒子法 (CSPH)，Dilts［27］基于伽辽金近似构造的移动最小二乘光滑粒子法(MLSPH)，Liu［2］在非连续区域两端分段泰勒展开构造的非连续光滑粒子法 (DSPH)，Liu 采用校正函数对核函数进行修正构造的再生核粒子方法(RKPM)。国内学者也开展了边界施加和边界缺陷的相关研究，丁桦［28］通过统计体积和边界粒子核函数的修正，提高了边界近似的精度，且构造了透射边界条件，汤文辉［29］、张嘉钟［30］构造了滑移和无滑移固壁条件。

以上对边界条件和边界缺陷的研究使得SPH方法的精度和适用范围大幅提高。然而，目前国际上仍没有边界处理的普适性方法，对于自由面、弹性边界、无反射边界等水下爆炸工程应用必需的技术还极少有人研究。因此，边界问题仍是SPH方法实现水下爆炸工程应用的主要技术瓶颈。

4 结语

SPH方法作为一种新兴的数值方法，在模拟大变形问题时较网格方法有重大优势，已成功应用于水下爆炸的模拟，但目前SPH方法仅局限于二维自由场水下爆炸的模拟，无法工程应用。SPH方法实现水下爆炸工程应用仍存在以下主要技术问题:

1)标准SPH方法模拟密度差异较小的多相流问题时可获得较高的精度，而大密度比多相流的模拟却极为困难，成为水下爆炸实现工程应用的主要技术瓶颈。国内外许多学者提出了多种方法模拟大密度比多相流问题，但其研究主要集中在气泡上浮、溃坝等低速流问题。而水下接触爆炸等具有大密度比、强冲击特征，状态方程对参数变化极为敏感，用以上方法进行模拟存在较大困难，需要对SPH算法加以改进。

2)用SPH方法模拟水下爆炸还仅局限于二维情况，采用SPH方法模拟三维水下爆炸的研究十分罕见。由于SPH方法较大的计算量，根据具体问题的对称性，采用一维球对称、二维轴对称或三维平面对称SPH算法模拟三维问题，是SPH方法实现工程应用的重要途径。国内外许多学者已开展了对称SPH算法的相关研究，但其研究主要为冲击管等简单问题的验证，且已形成的对称SPH算法并不成熟，其对称轴处的精度难以保证。因此，三维和基于对称算法的三维SPH方法模拟水下爆炸的研究需进一步开展。

3)边界条件难以施加、对边界处的粒子近似时存在缺陷是影响SPH方法计算精度的主要因素。已有大量的文献对边界条件和边界缺陷进行了研究，提高了SPH方法的精度和适用范围。然而，目前国际上仍没有边界处理的普适性方法，对于自由面、弹性边界、无反射边界等水下爆炸工程应用必需的技术还极少有人研究。因此，边界问题仍是SPH方法实现水下爆炸工程应用的主要技术瓶颈。

［1］SWEGLE J W，ATTAWAY S W.On the feasibility of using smoothed particle hydrodynamics for underwater explosion calculations［J］.Computational Mechanics，1995，17:151-168.

［2］LIU M B，LIU G R，LAM K Y.Meshfree particle simulation of the explosion process for high explosive in shaped charge［J］.Shock Waves，2003，12(6):509-520.

［3］杨刚，韩旭，龙述尧.应用 SPH方法模拟近水面爆炸［J］.工程力学，2008，25(4):204-208.YANG Gang，HAN Xu，LONG Shu-yao.Simulation of underwater explosion near air-water surface by SPH method［J］.Engineering Mechanics，2008，25(4):204-208.

［4］姚熊亮，杨文山，陈娟，等.沉底水雷爆炸威力的数值计算［J］.爆炸与冲击，2011，31(2):173-178.YAO Xiong-liang，YANG Wen-shan，CHEN Juan，et al.Numerical calculation of explosion power of mines lying on seabed［J］.Explosion and Shock Waves，2011，31(2):173-178.

［5］XU J X，LIU X L.Analysis of structural response under blast loads using the coupled SPH-FEM approach［J］.Journal of Zhejiang University Science A，2008，9(9):1184-1192.

［6］MONAGHAN J J.Implicit SPH drag and dusty gas dynamics［J］.Journal of Computational Physics.1997，138:801-820.

［7］CLEARY P W.Modelling confined multi-material heat and mass flows using SPH［J］.Applied Math.Model，1998，22:981-993.

［8］RITCHIE B W，THOMAS P A.Multiphase smoothed-particle hydrodynamics［J］.Mon.Not.R.Astron.Soc，2001，323:743-758.

［9］VALIZADEH A，SHAFIEEFAR M，MONAGHAN J J，et al.Modeling two-phase flows using SPH method［J］.J.Applied Sci，2008，8(21):3817-3826.

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

［11］HU X Y，ADAMS N A.A multi-phase SPH method for macroscopic and mesoscopic flows［J］.Journalof Computational Physics，2006，213:844-861.

［12］SOLENTHALER B，PAJAROLA R.Density contrast SPH interfaces［J］.Eurographics/ACM SIGGRAPH Symposium on Computer Animation.2008，211-218.

［13］杨刚，韩旭，龙述尧.SPH方法在两相流动问题中的典型应用［J］.湖南大学学报，2007，34(1):28-32.YANG Gang，HAN Xu，LONG Shu-yao.Typical application of SPH method to two-phase flow problems［J］.Journal of Hunan University，2007，34(1):28-32.

［14］HERNQUIST L，KATZ N.Tree SPH-A unification of SPH with the hierarchical tree method［J］.The Astrophysical Journal Supplement Series，1989，70:419-446.

［15］CLEARY P W，HA J，PRAKASH M，et al.3D SPH flow predictions and validation for high pressure die casting of automotive components［J］. Applied Mathematical Modeling，2006，30:1406-1427.

［16］DALRYMPLE R A，ROGERS B D.Numerical modeling of water waves with the SPH method［J］.Coastal Engineering，2006，53:141-147.

［17］LÖNNER R，YANG C，ON··ATE E.On the simulation of flows with violent free surface motion［J］.Computer Methods in Applied Mechanics and Engineering，2006，195:5597-5620.

［18］VUYST T D，VIGNJEVIC R，CAMPBELL J C.Coupling between meshless and finite element methods［J］.International Journal of Impact Engineering，2005，31:1054-1064.

［19］杨秀峰，彭世镠.溃坝流的光滑粒子法模拟［J］.计算物理，2010，27(2):173-180.YANG Xiu-feng，PENG Shi-liu.Simulation of dam-break flow with SPH method ［J］.Chinese Journal of Computational physice，2010，27(2):173-180.

［20］PETSCHEK A G，LIBERSKY L D.Cylindrical smoothed particle hydrodynamics［J］.JournalofComputational Physics.1993，109:76-83.

［21］BROOKSHAW L.Smooth particle hydrodynamics in cylindrical coordinates［J］.ANZIAM J，2003，44:114-139.

［22］OMANG M，BøRVE S，TRULSEN J.SPH in spherical and cylindrical coordinates［J］.JournalofComputational Physics，2006，213:391-412.

［23］龚凯.基于光滑质点水动力学 (SPH)方法的自由表面流动数值模拟研究 ［M］.上海:上海交通大学，2009.86-130.

［24］杨刚.光滑粒子法的改进及其若干典型应用 ［M］.长沙:湖南大学，2010.81-87.

［25］CAMPBELL P M.Some new algorithms for boundary values problems in smoothed particle hydrodynamics ［J］.DNA Report，1989，DNA-88-286.

［26］CHEN J K，BERAUN J E，CARNEY T C.A corrective smoothed particle method for boundary value problems in heat conduction ［J］.ComputerMethodsin Applied Mechanics and Engineering，1999，46:231-252.

［27］DILTS G A.Moving-least-squares-particle hydrodynamics-I:consistency and stability ［J］.International Journal for Numerical Methods in Engineering，1999，44:1115-1155.

［28］丁桦，龙丽平，伍彦峰.SPH方法在模拟线弹性波传播中的运用 ［J］.计算力学学报，2005，22(3):320-325.DING Hua，LONG Li-ping，WU Yan-feng.Application of smoothed particle hydrodynamics method to the simulations of elastic wave propagation in solid ［J］.Chinese Journal of Computational Mechanics，2005，22(3):320-325.

［29］汤文辉，毛益明.SPH数值模拟中固壁边界的一种处理方法［J］.国防科技大学学报，2001，23(6):54-57.TANG Wen-hui，MAO Yi-ming.A method of processing rigid boundary condition in SPH simulation ［J］.Journal of National University of Defense Technology，2001，23(6):54-57.

［30］张嘉钟，郑俊，于开平，等.碰撞-反射模型在SPH固壁条件中的应用 ［J］.应用基础与工程科学学报，2010，18(3):517-522.ZHANG Jia-zhong，ZHENG Jun，YU Kai-ping，et al.Collision-reflection model for modeling the wall boundary condition in SPH ［J］.Journal of Basic Science and Engineering，2010，18(3):517-522.