近自由面水下爆炸沖擊載荷特性三維數(shù)值模擬＊

明付仁，張阿漫，楊文山

（哈爾濱工程大學(xué)船舶工程學(xué)院，黑龍江 哈爾濱 150001）

近自由面水下爆炸研究對現(xiàn)代國防建設(shè)及海洋開發(fā)有著重要意義。近自由面水下爆炸會伴隨一系列復(fù)雜的物理現(xiàn)象（炸藥爆轟、水面隆起、水柱上升及翻卷和水滴飛濺等），這些都使傳統(tǒng)數(shù)值方法很難進(jìn)行模擬，現(xiàn)用多數(shù)數(shù)據(jù)均來自于實(shí)驗(yàn)，對近自由面水下爆炸的模擬尚不成熟。符松等［1］從拓?fù)鋵W(xué)角度出發(fā)，采用位標(biāo)函數(shù)處理二維近自由面水下爆炸模擬時(shí)混合網(wǎng)格及捕捉運(yùn)動界面問題，點(diǎn)明了水下爆炸的基本物理現(xiàn)象；師華強(qiáng)等［2］利用Level Set方法對二維近水面水下爆炸的多介質(zhì)運(yùn)動邊界進(jìn)行了捕捉，研究了沖擊波與剛性結(jié)構(gòu)的相互作用；楊剛等［3］利用SPH方法對二維近水面爆炸的基本物理現(xiàn)象進(jìn)行了研究，比較了空氣對近水面水下爆炸壓力場和密度場的影響；以上數(shù)值模擬方法均需克服由于大變形，存在運(yùn)動交界面等經(jīng)常導(dǎo)致計(jì)算崩潰的難題，而且多數(shù)僅僅限于二維數(shù)值模擬，與真實(shí)的三維情況存在差距。顧文彬等［4］通過對實(shí)驗(yàn)數(shù)據(jù)的統(tǒng)計(jì)分析，對淺層水中爆炸界面沖擊波的壓力峰值進(jìn)行了研究，但卻依賴于耗費(fèi)很大的實(shí)驗(yàn)研究。

本文中，采用無網(wǎng)格SPH（smoothed particle hydrodynamics）方法對近自由面水下爆炸的沖擊載荷進(jìn)行研究。SPH方法具有拉格朗日性質(zhì)，對大變形、多相物質(zhì)交界面以及運(yùn)動界面的捕捉等傳統(tǒng)網(wǎng)格方法處理具有很大困難的問題具有獨(dú)特的優(yōu)勢［5］。基于三維SPH方法，通過應(yīng)用改進(jìn)的變光滑長度鏈表搜索算法和多相物質(zhì)交界面處理方法，模擬無限域和近自由面水下爆炸的基本過程，通過對比分析，獲得研究近自由面水下爆炸的沖擊載荷特征，為工程提供參考。

1 基本理論

1.1 控制方程

SPH方法在模擬水下爆炸時(shí)，涉及多相物質(zhì)、阻抗不匹配等問題［6－7］，所以應(yīng)用適合處理密度比較大問題的SPH控制方程。

連續(xù)性方程、動量守恒方程、能量守恒方程分別為

式中：ρ、m、t、v、x、p、e、σ、ε和μ 分別代表密度、質(zhì)量、時(shí)間、速度矢量、坐標(biāo)、壓力、能量、總應(yīng)力張量分量、應(yīng)變率張量和動力粘性系數(shù)，W 表示光滑函數(shù)，Π采用Monaghan型人工粘度項(xiàng)［7］，下標(biāo)i、j表示相互作用的一對粒子，上標(biāo)α、β表示沿坐標(biāo)軸的方向。

1.2 變光滑長度鏈表搜索算法

SPH算法的最近相鄰粒子搜索法中，直接搜索算法、樹形搜索算法、鏈表搜索算法的復(fù)雜度階數(shù)分別為O（N2）、O（NlgN）、O（N），其中N為問題域內(nèi)的粒子總數(shù)。在模擬三維水下爆炸問題時(shí)，一般要求問題域比較大，粒子數(shù)N較大。

不同搜索算法耗時(shí)比較見表1，表中τ1為直接搜索算法平均每步耗時(shí)，τ2為變光滑長度鏈表算法平均每步耗時(shí)。由表可見，應(yīng)用鏈表搜索算法時(shí)優(yōu)勢明顯，計(jì)算效率會有很大提高。

表1 不同搜索算法不同粒子數(shù)平均每步耗時(shí)對比Table1 Average time－consuming every timestep of different search algorithms and particle numbers

然而，傳統(tǒng)鏈表搜索算法一般應(yīng)用在光滑長度為常量時(shí)，如果應(yīng)用到水下爆炸問題中，必須增加變光滑長度處理。

本文中提出變光滑長度的鏈表搜索算法，具體辦法為：在計(jì)算模型離散為粒子后，每個(gè)時(shí)間步進(jìn)行計(jì)算時(shí)，首先在問題域鋪設(shè)臨時(shí)網(wǎng)格并編號，網(wǎng)格單元的邊長根據(jù)總粒子數(shù)和平均每個(gè)網(wǎng)格單元內(nèi)的粒子數(shù)以及問題域的長寬高比例確定，然后遍歷每一個(gè)粒子，確定自身所處的網(wǎng)格坐標(biāo)。對于任意粒子i進(jìn)行相鄰粒子搜索時(shí)，根據(jù)自身的當(dāng)前光滑長度換算為格子光滑長度，格子光滑長度等于自身光滑長度與網(wǎng)格單元邊長的比值取整。然后根據(jù)自身所處網(wǎng)格坐標(biāo)和格子光滑長度，確定最近相鄰格子，再利用自身的光滑長度在最近相鄰格子內(nèi)確定用于自身物理量近似的粒子。所謂的變光滑長度主要體現(xiàn)在：一是按每個(gè)粒子的光滑長度確定自身的最近相鄰格子；二是由自身光滑長度確定最近相鄰格子內(nèi)的最近相鄰粒子。

1.3 交界面處理

水下爆炸涉及到爆轟產(chǎn)物與周圍水介質(zhì)的交界面，而粒子間所有信息的傳遞都是通過交界面，所以對于交界面的處理十分關(guān)鍵［6］。G.R.Liu等［7］、M.B.Liu等［8－9］針對水下爆炸問題提出了不同粒子交界面處理作用力，但它依賴光滑長度h，在處理水下爆炸問題時(shí)，由于粒子的光滑長度需要更新，所以在進(jìn)行優(yōu)化松弛的過程中光滑長度是變化的，會出現(xiàn)矯枉過正現(xiàn)象。所以利用上式進(jìn)行交界面處理時(shí)，由光滑長度項(xiàng)hi、hj確定力的大小會不穩(wěn)定，這樣很難保證交界面向外能量傳遞的穩(wěn)定。

本文中利用Lennard－Jones形式的分子力進(jìn)行交界面處理

式中：r0為截止半徑，一般與粒子的初始間距相近。

由式（4）可以發(fā)現(xiàn)，利用r0／rij確定交界面的力不存在由光滑長度引起的誤差，能夠保證較好的數(shù)值穩(wěn)定性。

2 數(shù)值研究

2.1 數(shù)值驗(yàn)證

為了將近自由面水下爆炸的沖擊載荷特征與無限域進(jìn)行對比分析，先對無限域水下爆炸的沖擊載荷進(jìn)行研究，數(shù)值模擬中采用6.7kg方形TNT裝藥在邊長a＝b＝2m的方形水域中心起爆。如圖1所示，TNT處于方形水域中心，l為邊長，R為測點(diǎn)距爆心的距離。典型爆距R的沖擊波壓力時(shí)程曲線如圖2所示，可以看出，在不同半徑處壓力峰值與經(jīng)驗(yàn)公式基本吻合，脈寬基本相同，沖擊波的壓力沖量大致相同，可以作為研究近自由面水下爆炸的參照對象。數(shù)值計(jì)算的結(jié)果由于間斷面的存在，在初始時(shí)刻壓力的上升階段會存在一定的坡度。同時(shí)，由于TNT爆轟后，產(chǎn)生的沖擊波向外傳播的同時(shí)，向藥心傳播著稀疏波，在稀疏波相遇后又形成沖擊波向外傳播，所以在壓力下降的時(shí)候又出現(xiàn)了小峰現(xiàn)象。數(shù)值計(jì)算結(jié)果存在一定的數(shù)值震蕩，主要原因是由于TNT與水的交界面處理時(shí)產(chǎn)生的數(shù)值不穩(wěn)定，如果通過調(diào)節(jié)人工粘度消除震蕩，會影響沖擊波的壓力峰值，所以本文通過調(diào)整相關(guān)參數(shù)，兼顧沖擊波的壓力峰值及數(shù)值震蕩，使數(shù)值模擬存在了一定的震蕩［10］。

圖1 無限域水下爆炸的粒子分布和測點(diǎn)分布Fig.1 Particle and gauging point distribution of underwater explosion in free field

圖2 不同爆距的壓力時(shí)程曲線Fig.2 Pressure－time curves at different distances

2.2 近自由面水下爆炸

2.2.1 數(shù)值模型

近自由面水下爆炸的數(shù)值模型如圖3所示，其中a＝2m，R為爆距，H為爆深，h為測深，b＝a／2＋H為變量，TNT藥量與無限域數(shù)值模擬相同。為了使不同測點(diǎn)的沖擊載荷特征的研究具有普遍意義，將所有長度對TNT邊長的一半l／2進(jìn)行處理，定義相對測深＝2h／l，相對爆深ˉH＝2 H／l，相對爆距ˉR＝2R／l，相對水平位置＝2x／l，相對垂向位置＝2z／l，其中x、z為距離藥包中心的水平坐標(biāo)和垂向坐標(biāo)，坐標(biāo)原點(diǎn)位于TNT中心。ps為近自由面水下爆炸壓力峰值，pf為無限域壓力峰值，＝ps／pf；Is為近自由面水下爆炸壓力沖量，If為無限域壓力沖量＝Is／If。

圖3 近自由面水下爆炸的粒子分布和測點(diǎn)分布Fig.3 Particle and gauging point distribution of underwater explosion near free surface

2.2.2 基本物理過程

近自由面水下爆炸中，當(dāng)TNT引爆后，產(chǎn)生爆轟波和沖擊波分別在TNT和水中傳播，隨著沖擊波的不斷向外傳播，產(chǎn)生的高溫高壓的爆轟產(chǎn)物推動著水面向上凸起，形成水柱，水柱的上升速度按指數(shù)衰減。從圖4可以看到，在此爆深很小時(shí)，產(chǎn)生的水柱較寬，而且很快破碎，噴濺范圍較大。而從圖5可知，在水面形成水冢的同時(shí)，沖擊波開始產(chǎn)生和傳播。沖擊波波頭具有較高的壓力，以球面波的形式向外傳播，當(dāng)沖擊波傳到水面時(shí)，立即反射稀疏波，向深水中傳播，追趕直達(dá)波，近自由面形成低壓區(qū)。整個(gè)過程中，沖擊波的壓力衰減迅速。

圖4 相對爆深ˉH＝1.25時(shí)水柱形成過程Fig.4 The water column formation when＝1.25

圖5 相對爆深＝3.75時(shí)沖擊波的傳播過程Fig.5 Shock wave propagation when＝3.75

2.2.3 沖擊波壓力峰值和沖量

根據(jù)不同反射區(qū)的劃分［11］，很容易捕捉到水面下的不同反射區(qū)的沖擊波，以＝3.75為例，如圖6所示。可以看出＝4.33＝1.25時(shí)沖擊波受水面影響產(chǎn)生直接截?cái)嘈?yīng)，并不影響峰值＝1，屬于規(guī)則反射區(qū)；當(dāng)＝7.50＝1.25時(shí)，水面的影響已經(jīng)超過卸載波脈寬的20%，而峰值幾乎未受影響≈1，屬于不規(guī)則反射區(qū)；當(dāng)ˉR＝9.19＝1.25時(shí)，水面反射的稀疏波追上直達(dá)波，削弱了沖擊波的波峰＜1，且減小了脈寬，屬于不規(guī)則反射區(qū)［12－13］。

圖6 不同反射區(qū)的沖擊波特征Fig.6 Shock wave characteristics of different reflection zones

為了充分認(rèn)識近自由面水下爆炸的沖擊波載荷特征，考慮爆深H的影響。近自由面與無限域水下爆炸的沖擊波壓力峰值比與壓力沖量比ˉI如圖7所示，TNT中心坐標(biāo)為（0，0）。TNT在近自由面引爆后，由于自由面存在截?cái)嘈?yīng)，當(dāng)反射波到達(dá)時(shí)間遲于直達(dá)波峰壓到達(dá)時(shí)間，不會削弱直達(dá)波峰值，就會形成規(guī)則反射區(qū)，此時(shí)ˉp＝1；若反射波到達(dá)時(shí)間先于直達(dá)波峰壓，那么會形成不規(guī)則反射區(qū)，在近自由面的部分區(qū)域會直接削弱沖擊波的壓力峰值，致使＜1，不規(guī)則反射區(qū)（削弱峰值區(qū)域）與規(guī)則反射區(qū)（未削弱峰值區(qū)域）將形成一個(gè)內(nèi)含裝藥的柱面。隨著爆深的增加，這個(gè)柱面越貼近水面，且曲率逐漸變小，削弱峰值作用逐漸變?nèi)酰衣拥礁h(yuǎn)距離。當(dāng)相對爆深為＝2.50時(shí)，峰值最大可削減為無限域的0.32倍，當(dāng)＝3.75和ˉH＝5.00時(shí)，將分別為0.44和0.52倍。

相比自由面對沖擊波壓力峰值的影響，自由面對沖擊波沖量的影響區(qū)域更大。從圖7中可以看出，以TNT為中心，壓力沖量比等值線將形成斗笠狀的柱面，隨著柱面的外張，越靠近水面，自由面對沖擊波沖量的削減作用越強(qiáng)。同一測深，沿水平方向的一定范圍內(nèi)，由于自由面的截?cái)啵瑳_量的削弱效果呈增強(qiáng)趨勢。同一水平距離不同測深時(shí)，向深水方向，沖量的衰減效果逐漸變?nèi)酢ｋS著爆深的增加，自由面的反射變?nèi)酰瑖?yán)重衰減區(qū)域越貼近自由面，同一沖量比形成的柱面越貼近自由面。當(dāng)相對爆深為＝2.50時(shí)，最大衰減為0.14倍，當(dāng)3.75和＝5.00時(shí)，將分別為0.17和0.21倍。

圖7 近自由面與無限域的壓力峰值比與壓力沖量比ˉIFig.7 Peak pressure ratioand impulse ratio of explosion near free surface and in free field

2.2.4 爆深對水柱產(chǎn)生初期的影響

在近自由面水下爆炸中，當(dāng)炸藥爆轟后，形成的爆轟產(chǎn)物和水混合在一起，在水面形成水柱。一般有兩種效應(yīng)：一是高而窄的柱狀噴柱，這類噴柱沿徑向飛濺較小，飛濺的角度很小；二是相對較寬的，當(dāng)水柱上升一段時(shí)間就開始破碎，沿徑向飛濺很大的噴柱［14］。

為了探索不同爆深對不同形式水柱形成的影響，圖8給出了600μs時(shí)典型爆深水柱形成初期現(xiàn)象。可以看出，TNT爆轟后產(chǎn)生的爆轟產(chǎn)物迅速膨脹，推動自由水面向上凸起。當(dāng)爆深較大時(shí)，這種推動作用會引起水面較大范圍的凸起，當(dāng)相對爆深＝5.00＝3.75時(shí)水面產(chǎn)生較小的隆起，而當(dāng)3.13、＝2.50時(shí)水面已經(jīng)產(chǎn)生了明顯的水冢現(xiàn)象和逐漸上升的水柱，此類水柱會形成高而窄的柱狀噴柱，徑向飛濺很小，向上運(yùn)動的趨勢很大；當(dāng)＝1.86時(shí)產(chǎn)生的水柱已經(jīng)破碎，到＝1.25產(chǎn)生的已經(jīng)是分布較廣、沿徑向飛濺水花的噴柱。

圖8 近自由面不同爆深的水柱初期現(xiàn)象Fig.8 Water column initial phenomena of different explosion depths near free surface

3 結(jié) 論

基于傳統(tǒng)的三維SPH方法，提出了變光滑長度鏈表搜索算法和多相物質(zhì)交界面處理技術(shù)，大大提高了計(jì)算效率，同時(shí)也保證了計(jì)算精度和數(shù)值的相對穩(wěn)定性；在此基礎(chǔ)上，模擬了無限域水下爆炸并與經(jīng)驗(yàn)公式進(jìn)行了對比分析，研究了近自由面水下爆炸的沖擊載荷特征，通過研究得出以下結(jié)論：

（1）自由面最大可將沖擊波壓力峰值削減為無限域的1／3左右，沖量的削弱可達(dá)無限域的1／7左右；自由面下不規(guī)則反射區(qū)沖擊波沖量的影響區(qū)域要遠(yuǎn)大于對壓力峰值的影響。

（2）隨著爆深的增加，自由面的截?cái)嘈?yīng)逐漸變?nèi)酰粔毫Ρ群蜎_量比形成的斗笠狀等值柱面逐漸向自由面外張；嚴(yán)重削弱區(qū)域，沿著水平方向向遠(yuǎn)離裝藥方向擴(kuò)散，沿垂直方向則向自由面貼近。

（3）隨著爆深的減小，近自由面水下爆炸形成的水柱由高而窄、沿徑向噴濺很小的水冢向破碎程度很大、沿徑向分布較廣的噴柱過渡。

［1］符松，王智平，張兆順，等.近水面水下爆炸的數(shù)值研究［J］.力學(xué)學(xué)報(bào)，1995，27（3）：267－275.FU Song，WANG Zhi－ping，ZHANG Zhao－shun，et al.Numerical study of underwater explosion near air－water surface［J］.Acta Mechanica Sinica，1995，27（3）：267－275.

［2］師華強(qiáng)，汪玉，宗智，等.近水面水下爆炸二維Level set數(shù)值模擬［J］.計(jì)算力學(xué)學(xué)報(bào)，2010，27（3）：409－414.SHI Hua－qiang，WANG Yu，ZONG Zhi，et al.2Dnumerical simulation of underwater explosion near free surface based on level set method［J］.Chinese Journal of Computational Mechanics，2010，27（3）：409－414.

［3］楊剛，韓旭，龍述堯.應(yīng)用SPH 方法模擬近水面爆炸［J］.工程力學(xué)，2008，25（4）：204－213.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－213.

［4］顧文彬，葉序雙，劉文華，等.界面對淺層水中爆炸沖擊波峰值壓力影響的研究［J］.解放軍理工大學(xué)學(xué)報(bào)，2001，2（5）：61－63.GU Wen－bin，YE Xu－shuang，LIU Wen－h(huán)ua，et al.Peak pressure investigation of exploding wave influenced by interfaces in shallow－layer water［J］.Journal of PLA University of Science and Technology，2001，2（5）：61－63.

［5］劉謀斌，宗智，常建忠.光滑粒子動力學(xué)方法的發(fā)展與應(yīng)用［J］.力學(xué)進(jìn)展，2011，41（2）：217－231.LIU Mou－bin，ZONG Zhi，CHANG Jian－zhong.Developments and applications of smoothed particle hydrodynamics［J］.Advances in Mechanics，2011，41（2）：217－231.

［6］師華強(qiáng)，宗智，賈敬蓓.水下爆炸沖擊波的近場特性［J］.爆炸與沖擊，2009，29（2）：125－129.SHI Hua－qiang，ZONG Zhi，JIA Jing－bei.Short－range characters of underwear blast waves［J］.Explosion and Shock Waves，2009，29（2）：125－129.

［7］Liu G R，Liu M B.Smoothed particle hydrodynamics：meshfree particle method［M］.Singapore：World Scientific，2003.

［8］Liu M B，Liu G R，Lam K Y.Investigations into water mitigation using a meshless particle method［J］.Shock Waves，2002，12：181－195.

［9］Liu M B，Liu G R，Zong Z，et al.Computer simulation of the high explosive explosion using smoothed particle hydrodynamics methodology［J］.Computer ＆ Fluids，2003，32（3）：305－322.

［10］劉謀斌，常建忠.光滑粒子動力學(xué)方法中粒子分布與數(shù)值穩(wěn)定性分析［J］.物理學(xué)報(bào)，2010，59（6）：3654－3659.LIU Mou－bin，CHANG Jian－zhong.Particle distribution and numerical stability in smoothed particle hydrodynamics method［J］.Acta Mechanica Sinica，2010，59（6）：3654－3659.

［11］Zamyshlyayev B V，Yakovlev Y S.Dynamic loads in underwater explosion［R］.AD0757183，1973.

［12］Cole R H.Underwater explosions［M］.New Jersy：Princeton University Press，1948.

［13］馮剛，朱錫，張振華.近水面水下爆炸作用下艦艇結(jié)構(gòu)損傷數(shù)值仿真［J］.計(jì)算機(jī)輔助工程，2006，15（4）：81－84.FENG Gang，ZHU Xi，ZHANG Zhen－h(huán)ua.Numerical analysis of hull structure damage due to underwater explosion near water surface［J］.Computer Aided Engineering，2006，15（4）：81－84.

［14］郅斌偉，張志江，李健，等.近水面水下爆炸水柱效應(yīng)研究［J］.北京理工大學(xué)學(xué)報(bào)，2009，29（1）：5－8.ZHI Bin－wei，ZHANG Zhi－jiang，LI Jian，et al.A study on water columns produced by near water surface explosion［J］.Transactions of Beijing Institute of Technology，2009，29（1）：5－8.