楊文山，孟曉宇，王祖華

(武漢第二船舶設(shè)計(jì)研究所，湖北武漢430064)

0 引言

隨著SPH(Smoothed Particle Hydrodynamics)方法在具有材料強(qiáng)度的動(dòng)態(tài)響應(yīng)問(wèn)題和具有大變形流體動(dòng)力學(xué)問(wèn)題上的成功應(yīng)用，許多學(xué)者開(kāi)始考慮用其模擬水下爆炸的可行性。Swegle［1］最早對(duì)SPH方法模擬水下爆炸的可行性進(jìn)行了論證，發(fā)現(xiàn)SPH的無(wú)網(wǎng)格性質(zhì)適合于解決水下爆炸產(chǎn)生的大變形問(wèn)題，而其拉格朗日性質(zhì)極易捕捉爆炸運(yùn)動(dòng)物質(zhì)交界面和自由面，在水下爆炸的模擬上具有先天優(yōu)勢(shì)，非常適合水下爆炸沖擊波階段的模擬。Liu［2］在SPH方法模擬水下爆炸領(lǐng)域作出了突出貢獻(xiàn)，對(duì)一維炸藥爆轟、二維爆炸氣體膨脹、錐形炸藥爆炸、二維水下爆炸、水介質(zhì)緩沖等問(wèn)題進(jìn)行了模擬。然而其模擬均是針對(duì)自由場(chǎng)爆炸等比較理想的模型，僅能對(duì)算法的穩(wěn)定性、精度、收斂性等進(jìn)行驗(yàn)證，且僅局限于二維水下爆炸的模擬，無(wú)法實(shí)現(xiàn)工程應(yīng)用。

國(guó)內(nèi)楊剛［3］、姚熊亮［4］分別對(duì)近自由面水下爆炸和沉底水下爆炸進(jìn)行了模擬，Xu［5］采用SPH和有限元結(jié)合的方法模擬了結(jié)構(gòu)在水下爆炸作用下的響應(yīng)。以上水下爆炸的模擬均是在一維或二維自由場(chǎng)條件下進(jìn)行的，且僅涉及到炸藥和水等密度差別較小的物質(zhì)，而工程中的水下爆炸均屬于三維問(wèn)題，存在著自由面、彈性邊界及無(wú)反射邊界，且多為水、空氣、爆炸氣體、鋼等多種密度差異較大物質(zhì)混合的多相流問(wèn)題。以上原因使得SPH方法難以解決許多工程問(wèn)題，例如:艦船接觸爆炸、艦船防雷艙爆炸的模擬對(duì)艦船防護(hù)結(jié)構(gòu)設(shè)計(jì)具有重要意義，但其涉及TNT裝藥、空氣、水、鋼等性質(zhì)差異較大物質(zhì)的相互作用，且鋼板厚度遠(yuǎn)小于水域和TNT裝藥的尺度，使其粒子難以均勻分布，材料和粒子的高度非均勻性使SPH方法模擬接觸爆炸存在較大困難;工程中常用的柱形和方形TNT裝藥爆炸均屬于三維問(wèn)題，且需考慮自由面和無(wú)反射邊界的影響，而當(dāng)前的研究?jī)H局限于二維自由場(chǎng)水下爆炸的模擬，難以解決實(shí)際工程問(wèn)題。因此，大密度比多相流模擬、三維或基于對(duì)稱算法的三維模擬、無(wú)反射邊界等邊界的處理是SPH方法實(shí)現(xiàn)水下爆炸工程應(yīng)用的主要技術(shù)瓶頸，需要進(jìn)行深入研究。

1 SPH模擬多相流問(wèn)題

由于SPH方法是拉格朗日性質(zhì)和粒子性質(zhì)的和諧結(jié)合，容易追蹤界面，因此特別適合模擬多相流等異種介質(zhì)交界面不斷運(yùn)動(dòng)的問(wèn)題。Monaghan［6］，Cleary［7］和 Ritchie［8］分別對(duì)火山噴發(fā)氣塵兩相運(yùn)動(dòng)、多物質(zhì)熱傳導(dǎo)、多物質(zhì)星云形成等多相流問(wèn)題進(jìn)行了模擬，并發(fā)現(xiàn)密度差別較小的多相流問(wèn)題可認(rèn)為其密度梯度小于光滑函數(shù)梯度，采用標(biāo)準(zhǔn)SPH方法進(jìn)行模擬便可達(dá)到精度要求。然而，當(dāng)多相流中的物質(zhì)密度差異較大時(shí)，交界面的計(jì)算精度則得不到保證。這是由于當(dāng)高密度介質(zhì)和低密度介質(zhì)互為交界面時(shí)，低密度介質(zhì)會(huì)因較大的密度梯度使近似結(jié)果偏大，而高密度介質(zhì)的近似值比實(shí)際值偏小，這種誤差的累積會(huì)使計(jì)算值失真，嚴(yán)重時(shí)會(huì)導(dǎo)致計(jì)算程序的發(fā)散。

許多學(xué)者對(duì)大密度比的多相流問(wèn)題進(jìn)行了研究，Valizadeh［9］對(duì)交界面處的粒子質(zhì)量進(jìn)行修正，模擬大密度比的多相流問(wèn)題，但其修正系數(shù)的取值需隨密度比的變化而變化，因此該方法通用性較差。Colagrossi［10］通過(guò)改進(jìn)SPH的近似方式，并周期性地初始化密度獲得了良好的模擬效果。Hu［11］將體積近似代替密度近似應(yīng)用于連續(xù)密度方程和動(dòng)量方程，來(lái)求解大密度比的多相流問(wèn)題。Solenthaler［12］采用“粒子密度”近似密度、壓力和粘性力亦得到了良好的效果。韓旭［13］對(duì)密度近似進(jìn)行了正則化，并對(duì)氣體狀態(tài)方程加入了范德華修正項(xiàng)成功模擬了氣泡在水中的上浮運(yùn)動(dòng)。以上研究均是對(duì)氣泡上浮、潰壩等低速流問(wèn)題進(jìn)行模擬，而水下爆炸為大密度比高速多相流問(wèn)題，物理模型中空氣和鋼的密度比達(dá)1/7800，且爆炸產(chǎn)生的物質(zhì)運(yùn)動(dòng)速度可達(dá)km/s，和低速流有較大差異，因此大密度比高速多相流的模擬是實(shí)現(xiàn)水下爆炸工程應(yīng)用的重要技術(shù)。

2 SPH模擬三維問(wèn)題

SPH方法的一大特點(diǎn)是其程序結(jié)構(gòu)簡(jiǎn)單，可自然擴(kuò)展到三維，但用SPH方法模擬三維問(wèn)題需要克服其計(jì)算量較大的缺點(diǎn)。影響SPH方法計(jì)算量的主要是粒子搜索，這是由于SPH方法是采用i粒子周邊支持域內(nèi)的粒子對(duì)i粒子進(jìn)行近似，且整個(gè)問(wèn)題域中的粒子是不斷運(yùn)動(dòng)的，因此在每個(gè)時(shí)間步內(nèi)均需要對(duì)粒子進(jìn)行搜索以確定i粒子對(duì)應(yīng)的近似粒子，使得粒子搜索成為影響SPH方法計(jì)算量的決定因素。在SPH研究的初期采用全配對(duì)搜索法搜索粒子，該算法的形式最為簡(jiǎn)單，但其計(jì)算量為O(N2)(N為問(wèn)題域的總粒子數(shù))，在計(jì)算粒子數(shù)較多的問(wèn)題時(shí)速度較慢。為提高粒子的搜索效率，Monaghan［6］提出了鏈表搜索法，該方法將粒子分布在胞元內(nèi)，通過(guò)胞元的存儲(chǔ)記錄搜索粒子，理想狀態(tài)下該方法的計(jì)算量為O(N)，但該方法對(duì)于問(wèn)題域尺度在計(jì)算過(guò)程中大幅擴(kuò)張、光滑長(zhǎng)度發(fā)生巨大變化的情況搜索效率較低。Hernquist［14］提出了適合求解變光滑長(zhǎng)度的樹(shù)形搜索法，該方法將問(wèn)題域分成多個(gè)卦限，直至每個(gè)卦限僅包含1個(gè)粒子，通過(guò)有序的樹(shù)形卦限對(duì)粒子進(jìn)行搜索，其計(jì)算量為O(N1gN)，樹(shù)形搜索算法的應(yīng)用在保證精度的前提下使計(jì)算效率大幅提高，采用該算法已可對(duì)數(shù)百萬(wàn)粒子的問(wèn)題進(jìn)行模擬。粒子搜索算法計(jì)算效率的提高使得SPH模擬三維問(wèn)題成為可能。Cleary［15］采用三維SPH方法模擬了輕金屬的高壓壓鑄過(guò)程;Dalrymple［16］模擬了潰壩水流對(duì)結(jié)構(gòu)的作用;L?nner［17］對(duì) 三 維 自 由 表 面 流 問(wèn) 題 進(jìn) 行 了 模 擬;Vuyst［18］采用有限元和SPH結(jié)合的方法模擬了三維高速?zèng)_擊問(wèn)題;國(guó)內(nèi)楊秀峰［19］等人也對(duì)三維潰壩問(wèn)題進(jìn)行了模擬。

模擬三維問(wèn)題需要的計(jì)算量較大，當(dāng)三維問(wèn)題具有對(duì)稱特征時(shí)，可采用一維球?qū)ΨQ、二維軸對(duì)稱或三維平面對(duì)稱算法進(jìn)行模擬，從而使計(jì)算量大幅降低。對(duì)稱SPH算法中用1個(gè)粒子代表周向所有粒子，故在對(duì)稱軸處粒子質(zhì)量分布極不均勻，導(dǎo)致對(duì)稱軸處的計(jì)算精度較低，因此其研究遠(yuǎn)沒(méi)有標(biāo)準(zhǔn)SPH算法成熟。Petschek 等對(duì)三維笛卡兒坐標(biāo)系中的角坐標(biāo)積分，將三維坐標(biāo)系轉(zhuǎn)化為二維軸對(duì)稱柱坐標(biāo)系，但其計(jì)算結(jié)果在對(duì)稱軸處存在問(wèn)題。Brookshaw［21］直接應(yīng)用二維軸對(duì)稱柱坐標(biāo)系的控制方程將三維問(wèn)題轉(zhuǎn)化為二維問(wèn)題，但該方法的精度在對(duì)稱軸附近得不到保證。Omang［22］在不同的坐標(biāo)系下推導(dǎo)了對(duì)稱問(wèn)題的核函數(shù)，然后采用拉格朗日公式離散了對(duì)稱問(wèn)題的控制方程，其計(jì)算結(jié)果在對(duì)稱軸處仍具有較高的精度。國(guó)內(nèi)龔凱［23］、楊剛［24］等人也開(kāi)發(fā)了對(duì)稱SPH方法計(jì)算程序，但均沒(méi)有討論對(duì)稱軸處計(jì)算結(jié)果的精度。從以上研究可知，基于對(duì)稱算法的三維SPH方法仍不成熟，且三維及基于對(duì)稱算法的三維水下爆炸模擬的相關(guān)研究極少。

3 SPH處理邊界問(wèn)題

邊界條件難以施加、對(duì)邊界處的粒子近似時(shí)存在缺陷是影響SPH方法計(jì)算精度的主要因素，這使得學(xué)術(shù)界一直致力于施加邊界條件和修正邊界缺陷的研究。Monaghan［6］處理固壁邊界時(shí)在固壁處施加一排虛粒子，通過(guò)虛粒子的排斥力防止粒子穿透固壁。Libersky以自由面為對(duì)稱面反射對(duì)稱虛粒子，通過(guò)虛粒子的近似修正自由面邊界處因缺少近似粒子導(dǎo)致的缺陷。Randles將虛粒子的所有參變量賦予和邊界處粒子等量的值，通過(guò)虛粒子和邊界內(nèi)粒子插值得到內(nèi)部粒子的變量值。Liu［2］在Monaghan和Libersky的基礎(chǔ)上采用2種虛粒子處理固壁，第一種虛粒子布置在固壁邊界上對(duì)內(nèi)部粒子施加排斥力，第二種虛粒子分布在邊界外參與粒子近似，該方法防止了粒子穿透邊界，且提高了邊界處的近似精度。以上方法均采用了虛粒子，可以解決近似邊界附近粒子時(shí)由于粒子不充足而引起的邊界缺陷問(wèn)題。

導(dǎo)致SPH方法邊界缺陷的另一個(gè)原因是核近似導(dǎo)函數(shù)時(shí)在邊界處會(huì)發(fā)生截?cái)?，使分部積分的面積分項(xiàng)不為0。為此，Campbell［25］將邊界余項(xiàng)引入到分部積分中，修正了邊界的截?cái)鄦?wèn)題。此外，一些SPH方法的改進(jìn)形式對(duì)修正邊界缺陷、提高計(jì)算精度起到了良好效果，如Chen［26］根據(jù)泰勒展開(kāi)構(gòu)造的正則化形式的修正光滑粒子法 (CSPH)，Dilts［27］基于伽遼金近似構(gòu)造的移動(dòng)最小二乘光滑粒子法(MLSPH)，Liu［2］在非連續(xù)區(qū)域兩端分段泰勒展開(kāi)構(gòu)造的非連續(xù)光滑粒子法 (DSPH)，Liu 采用校正函數(shù)對(duì)核函數(shù)進(jìn)行修正構(gòu)造的再生核粒子方法(RKPM)。國(guó)內(nèi)學(xué)者也開(kāi)展了邊界施加和邊界缺陷的相關(guān)研究，丁樺［28］通過(guò)統(tǒng)計(jì)體積和邊界粒子核函數(shù)的修正，提高了邊界近似的精度，且構(gòu)造了透射邊界條件，湯文輝［29］、張嘉鐘［30］構(gòu)造了滑移和無(wú)滑移固壁條件。

以上對(duì)邊界條件和邊界缺陷的研究使得SPH方法的精度和適用范圍大幅提高。然而，目前國(guó)際上仍沒(méi)有邊界處理的普適性方法，對(duì)于自由面、彈性邊界、無(wú)反射邊界等水下爆炸工程應(yīng)用必需的技術(shù)還極少有人研究。因此，邊界問(wèn)題仍是SPH方法實(shí)現(xiàn)水下爆炸工程應(yīng)用的主要技術(shù)瓶頸。

4 結(jié)語(yǔ)

SPH方法作為一種新興的數(shù)值方法，在模擬大變形問(wèn)題時(shí)較網(wǎng)格方法有重大優(yōu)勢(shì)，已成功應(yīng)用于水下爆炸的模擬，但目前SPH方法僅局限于二維自由場(chǎng)水下爆炸的模擬，無(wú)法工程應(yīng)用。SPH方法實(shí)現(xiàn)水下爆炸工程應(yīng)用仍存在以下主要技術(shù)問(wèn)題:

1)標(biāo)準(zhǔn)SPH方法模擬密度差異較小的多相流問(wèn)題時(shí)可獲得較高的精度，而大密度比多相流的模擬卻極為困難，成為水下爆炸實(shí)現(xiàn)工程應(yīng)用的主要技術(shù)瓶頸。國(guó)內(nèi)外許多學(xué)者提出了多種方法模擬大密度比多相流問(wèn)題，但其研究主要集中在氣泡上浮、潰壩等低速流問(wèn)題。而水下接觸爆炸等具有大密度比、強(qiáng)沖擊特征，狀態(tài)方程對(duì)參數(shù)變化極為敏感，用以上方法進(jìn)行模擬存在較大困難，需要對(duì)SPH算法加以改進(jìn)。

2)用SPH方法模擬水下爆炸還僅局限于二維情況，采用SPH方法模擬三維水下爆炸的研究十分罕見(jiàn)。由于SPH方法較大的計(jì)算量，根據(jù)具體問(wèn)題的對(duì)稱性，采用一維球?qū)ΨQ、二維軸對(duì)稱或三維平面對(duì)稱SPH算法模擬三維問(wèn)題，是SPH方法實(shí)現(xiàn)工程應(yīng)用的重要途徑。國(guó)內(nèi)外許多學(xué)者已開(kāi)展了對(duì)稱SPH算法的相關(guān)研究，但其研究主要為沖擊管等簡(jiǎn)單問(wèn)題的驗(yàn)證，且已形成的對(duì)稱SPH算法并不成熟，其對(duì)稱軸處的精度難以保證。因此，三維和基于對(duì)稱算法的三維SPH方法模擬水下爆炸的研究需進(jìn)一步開(kāi)展。

3)邊界條件難以施加、對(duì)邊界處的粒子近似時(shí)存在缺陷是影響SPH方法計(jì)算精度的主要因素。已有大量的文獻(xiàn)對(duì)邊界條件和邊界缺陷進(jìn)行了研究，提高了SPH方法的精度和適用范圍。然而，目前國(guó)際上仍沒(méi)有邊界處理的普適性方法，對(duì)于自由面、彈性邊界、無(wú)反射邊界等水下爆炸工程應(yīng)用必需的技術(shù)還極少有人研究。因此，邊界問(wèn)題仍是SPH方法實(shí)現(xiàn)水下爆炸工程應(yīng)用的主要技術(shù)瓶頸。

［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］楊剛，韓旭，龍述堯.應(yīng)用 SPH方法模擬近水面爆炸［J］.工程力學(xué)，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］姚熊亮，楊文山，陳娟，等.沉底水雷爆炸威力的數(shù)值計(jì)算［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方法在兩相流動(dòng)問(wèn)題中的典型應(yīng)用［J］.湖南大學(xué)學(xué)報(bào)，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］.計(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］龔凱.基于光滑質(zhì)點(diǎn)水動(dòng)力學(xué) (SPH)方法的自由表面流動(dòng)數(shù)值模擬研究 ［M］.上海:上海交通大學(xué)，2009.86-130.

［24］楊剛.光滑粒子法的改進(jìn)及其若干典型應(yīng)用 ［M］.長(zhǎng)沙:湖南大學(xué)，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方法在模擬線彈性波傳播中的運(yùn)用 ［J］.計(jì)算力學(xué)學(xué)報(bào)，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數(shù)值模擬中固壁邊界的一種處理方法［J］.國(guó)防科技大學(xué)學(xué)報(bào)，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］張嘉鐘，鄭俊，于開(kāi)平，等.碰撞-反射模型在SPH固壁條件中的應(yīng)用 ［J］.應(yīng)用基礎(chǔ)與工程科學(xué)學(xué)報(bào)，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.