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

明付仁，張阿漫，楊文山

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

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

本文中，采用無網格SPH（smoothed particle hydrodynamics）方法對近自由面水下爆炸的沖擊載荷進行研究。SPH方法具有拉格朗日性質，對大變形、多相物質交界面以及運動界面的捕捉等傳統網格方法處理具有很大困難的問題具有獨特的優勢［5］?；谌SSPH方法，通過應用改進的變光滑長度鏈表搜索算法和多相物質交界面處理方法，模擬無限域和近自由面水下爆炸的基本過程，通過對比分析，獲得研究近自由面水下爆炸的沖擊載荷特征，為工程提供參考。

1 基本理論

1.1 控制方程

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

連續性方程、動量守恒方程、能量守恒方程分別為

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

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

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

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

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

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

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

1.3 交界面處理

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

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

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

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

2 數值研究

2.1 數值驗證

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

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

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

2.2 近自由面水下爆炸

2.2.1 數值模型

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

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

2.2.2 基本物理過程

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

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

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

2.2.3 沖擊波壓力峰值和沖量

根據不同反射區的劃分［11］，很容易捕捉到水面下的不同反射區的沖擊波，以＝3.75為例，如圖6所示?？梢钥闯觯?.33＝1.25時沖擊波受水面影響產生直接截斷效應，并不影響峰值＝1，屬于規則反射區；當＝7.50＝1.25時，水面的影響已經超過卸載波脈寬的20%，而峰值幾乎未受影響≈1，屬于不規則反射區；當ˉR＝9.19＝1.25時，水面反射的稀疏波追上直達波，削弱了沖擊波的波峰＜1，且減小了脈寬，屬于不規則反射區［12－13］。

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

為了充分認識近自由面水下爆炸的沖擊波載荷特征，考慮爆深H的影響。近自由面與無限域水下爆炸的沖擊波壓力峰值比與壓力沖量比ˉI如圖7所示，TNT中心坐標為（0，0）。TNT在近自由面引爆后，由于自由面存在截斷效應，當反射波到達時間遲于直達波峰壓到達時間，不會削弱直達波峰值，就會形成規則反射區，此時ˉp＝1；若反射波到達時間先于直達波峰壓，那么會形成不規則反射區，在近自由面的部分區域會直接削弱沖擊波的壓力峰值，致使＜1，不規則反射區（削弱峰值區域）與規則反射區（未削弱峰值區域）將形成一個內含裝藥的柱面。隨著爆深的增加，這個柱面越貼近水面，且曲率逐漸變小，削弱峰值作用逐漸變弱，且蔓延到更遠距離。當相對爆深為＝2.50時，峰值最大可削減為無限域的0.32倍，當＝3.75和ˉH＝5.00時，將分別為0.44和0.52倍。

相比自由面對沖擊波壓力峰值的影響，自由面對沖擊波沖量的影響區域更大。從圖7中可以看出，以TNT為中心，壓力沖量比等值線將形成斗笠狀的柱面，隨著柱面的外張，越靠近水面，自由面對沖擊波沖量的削減作用越強。同一測深，沿水平方向的一定范圍內，由于自由面的截斷，沖量的削弱效果呈增強趨勢。同一水平距離不同測深時，向深水方向，沖量的衰減效果逐漸變弱。隨著爆深的增加，自由面的反射變弱，嚴重衰減區域越貼近自由面，同一沖量比形成的柱面越貼近自由面。當相對爆深為＝2.50時，最大衰減為0.14倍，當3.75和＝5.00時，將分別為0.17和0.21倍。

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

2.2.4 爆深對水柱產生初期的影響

在近自由面水下爆炸中，當炸藥爆轟后，形成的爆轟產物和水混合在一起，在水面形成水柱。一般有兩種效應：一是高而窄的柱狀噴柱，這類噴柱沿徑向飛濺較小，飛濺的角度很?。欢窍鄬^寬的，當水柱上升一段時間就開始破碎，沿徑向飛濺很大的噴柱［14］。

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

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

3 結 論

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

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

（2）隨著爆深的增加，自由面的截斷效應逐漸變弱；壓力比和沖量比形成的斗笠狀等值柱面逐漸向自由面外張；嚴重削弱區域，沿著水平方向向遠離裝藥方向擴散，沿垂直方向則向自由面貼近。

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

［1］符松，王智平，張兆順，等.近水面水下爆炸的數值研究［J］.力學學報，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］師華強，汪玉，宗智，等.近水面水下爆炸二維Level set數值模擬［J］.計算力學學報，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］楊剛，韓旭，龍述堯.應用SPH 方法模擬近水面爆炸［J］.工程力學，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］.解放軍理工大學學報，2001，2（5）：61－63.GU Wen－bin，YE Xu－shuang，LIU Wen－hua，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］劉謀斌，宗智，常建忠.光滑粒子動力學方法的發展與應用［J］.力學進展，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］師華強，宗智，賈敬蓓.水下爆炸沖擊波的近場特性［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］劉謀斌，常建忠.光滑粒子動力學方法中粒子分布與數值穩定性分析［J］.物理學報，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］馮剛，朱錫，張振華.近水面水下爆炸作用下艦艇結構損傷數值仿真［J］.計算機輔助工程，2006，15（4）：81－84.FENG Gang，ZHU Xi，ZHANG Zhen－hua.Numerical analysis of hull structure damage due to underwater explosion near water surface［J］.Computer Aided Engineering，2006，15（4）：81－84.

［14］郅斌偉，張志江，李健，等.近水面水下爆炸水柱效應研究［J］.北京理工大學學報，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.