楊文山，孟曉宇，王祖華

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

0 引言

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

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

1 SPH模擬多相流問題

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

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

2 SPH模擬三維問題

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

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

3 SPH處理邊界問題

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

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

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

4 結(jié)語

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

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

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

3)邊界條件難以施加、對邊界處的粒子近似時存在缺陷是影響SPH方法計算精度的主要因素。已有大量的文獻(xiàn)對邊界條件和邊界缺陷進(jìn)行了研究，提高了SPH方法的精度和適用范圍。然而，目前國際上仍沒有邊界處理的普適性方法，對于自由面、彈性邊界、無反射邊界等水下爆炸工程應(yīng)用必需的技術(shù)還極少有人研究。因此，邊界問題仍是SPH方法實現(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］.爆炸與沖擊，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方法在兩相流動問題中的典型應(yīng)用［J］.湖南大學(xué)學(xué)報，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］龔凱.基于光滑質(zhì)點水動力學(xué) (SPH)方法的自由表面流動數(shù)值模擬研究 ［M］.上海:上海交通大學(xué)，2009.86-130.

［24］楊剛.光滑粒子法的改進(jìn)及其若干典型應(yīng)用 ［M］.長沙:湖南大學(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方法在模擬線彈性波傳播中的運用 ［J］.計算力學(xué)學(xué)報，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］.國防科技大學(xué)學(xué)報，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固壁條件中的應(yīng)用 ［J］.應(yīng)用基礎(chǔ)與工程科學(xué)學(xué)報，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.