艦船水下接觸爆炸的SPH算法研究

明付仁，張阿漫，楊文山，劉翠丹

（1.哈爾濱工程大學 船舶工程學院，哈爾濱 150001；2.武漢第二船舶設計研究所，武漢 430064）

艦船水下接觸爆炸是高能炸藥爆轟、高溫高壓爆轟產物的膨脹、以及水中沖擊結構的復雜過程。整個過程存在特別的性質如大變形、運動物質交界面、高度不均勻、可變形邊界和自由面等［1］，這使得傳統有網格算法常常因為網格變形扭曲等因素導致計算崩潰，而SPH算法具有無網格性質和拉格朗日性質，將問題域離散為攜帶物理量的粒子進行研究，能克服計算中與大變形等相關的困難。

SPH算法在進行數值模擬時，粒子的各物理量是通過光滑函數對緊支域內粒子進行光滑近似得到的，故當問題域內多相物質密度（或體積）極大不均勻時，進行粒子近似時很容易出問題，選取不同的SPH控制方程很關鍵。Ott等［2］給出了解決大密度梯度的連續性方程，但對于密度（或體積）極大不均勻問題尚沒有提出完備的思想。傳統SPH算法的最近相鄰粒子搜索法（NNPS）中，直接搜索算法、樹形搜索算法、鏈表搜索算法的復雜度階數［3］分別為O（N2）、O（NlgN）、O（N），其中N為問題域內的粒子總數。可見當N越大時，鏈表搜索算法優勢越明顯，但是傳統鏈表搜索算法一般應用在光滑長度為常量時，若要將其應用到水下爆炸問題中，必須增加變光滑長度處理。SPH算法將問題域離散為攜帶質量、密度、壓力等物理量的粒子，當研究對象為不規則形狀時，粒子初始化時比較困難。張姝慧等［4］進行了初始粒子配置分析，提出粒子應盡量分布均勻或者連續分布，但并未給出具體實施辦法；倪國喜等［5］提出了任意區域的粒子分布方法，雖然為均勻分布粒子方法，但三角剖分Delaunay方法計算相對復雜、效率很低，故需要提出具有更強適應性的初始化粒子分布的方法。利用SPH方法來研究含結構的水下接觸爆炸問題具有深遠的意義。

1 理論研究

1.1 SPH控制方程

SPH方法在模擬水下接觸爆炸問題時，存在炸藥爆轟產物、水、結構等多相物質，考慮結構材料的SPH控制方程可以分為兩種形式：一是粒子的近似式含有緊支域內粒子的體積，即具有體積思想的控制方程；另一種形式是近似式含有緊支域內粒子的質量，即具有質量思想的控制方程。

表1 不同的SPH控制方程Tab.1 Different SPH equations

表1 中：ρ，m，t，v，x，p，e，S，ε分別代表密度、質量、時間、速度、坐標位置、同性壓力、能量、偏應力張量、應變率張量，σαβ=-pδab+Sab為總應力張量，W表示光滑函數，Π表示人工粘度項，下標i、j表示相互作用的一對粒子，上標a、b表示沿坐標軸的方向。

在SPH方法近似中，對于任意粒子i，它自身的物理量是與緊支域內粒子息息相關的。具有體積思想的SPH控制方程，其中mj/ρj項表示的是緊支域內用于近似的粒子的體積，當粒子為近似均勻分布時，體積近似相同，這樣對于處理密度極大不均勻問題，粒子的物理量近似與最近相鄰粒子的質量、密度無關，具有體積思想的SPH控制方程具有很大優勢；而具有質量思想的SPH控制方程，粒子i的物理量與鄰近粒子的質量mj相關，當彼此密度（質量）差異較大時，即使粒子均勻分布，也將導致粒子i的各物理量近似存在問題，特別是粒子自身密度的變化，很可能導致計算的崩潰。需要說明的是，表1中控制方程針對的是考慮結構材料強度的SPH控制方程，當計算流體時，σab=-pδab+τab，τab為粘性剪應力。相應的本文應用的針對流體的能量方程如下（μ為動力粘性系數）：

建立2.1節中的數值模型，分別利用表1中具有體積思想與具有質量思想的SPH控制方程進行模擬，圓形裝藥中心起爆。模型中各部分具體參數見表2，對比結果如圖1所示。兩者除了控制方程應用不同外，其他各量均相同。結果表明，應用具有質量思想的SPH控制方程的程序運算到88步時就會發散，而應用具有體積的思想的控制方程程序可以保證運算無限運行下去。究其原因，這主要是由于當炸藥爆轟產生高溫高壓氣體后，在整個問題域內，各相物質的密度比增大，在以后的逐步近似過程中，密度近似誤差的累積最終導致了計算的崩潰。

表2 數值模型參數Tab.2 Numerical model parameters

圖1 應用不同控制方程的粒子分布圖，時間t為88 μs（88步）Fig.1 Particle distribution of different equations，t=88 μs

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

在解決水下接觸爆炸問題時，光滑長度的更新是必要的。然而在鄰近粒子搜索算法中，相對于全配對搜索算法和樹形搜索算法，鏈表搜索算法的復雜度階數最低，但鏈表搜索算法在光滑長度為空間常量時才非常有效。為了在提高計算效率的同時，又保證計算精度，本文提出變光滑長度的鏈表搜索算法。

具體辦法如下：如圖2所示，在進行鏈表搜索算法時，要在問題域鋪設臨時網格，并在坐標軸方向進行編號，網格單元的邊長是相同的，尺寸根據問題域的長寬高比例確定的。對于任意粒子i，首先將當前自身的光滑長度換算為“格子光滑長度”，即格子光滑長度等于自身光滑長度/網格單元邊長，然后取整數。再根據格子光滑長度確定粒子i進行臨近粒子搜索的網格單元，最后根據自身的光滑長度在網格單元內劃定緊支域進行粒子光滑近似。變光滑長度鏈表搜索算法對于每個粒子都按自身的光滑長度進行近似，在每個時間步結束后可按Liu［3］推導的廣義光滑長度進行更新，從而進行下一時間步的搜索計算。

圖2 變光滑長度鏈表搜索算法Fig.2 Linked-list search algorithm with variable smoothing length

2.1節數值模型進行研究時，問題域離散為不同數目粒子，使用全配對搜索算法與變光滑長度的鏈表搜索算法耗時對比如表3所示，所用計算機為Pentium Dual Core E6500 CPU，主頻2.93 GHz。從表3可以看出，隨著粒子數的增多，全配對搜索算法的計算效率是近似為二次拋物線增長，變光滑長度的鏈表搜索算法近似為直線增長，當粒子數越多時，變光滑長度的鏈表搜索算法的效率越明顯。

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

1.3 SPH方法的前、后期處理

SPH方法在進行粒子初始化時，對于小尺度、不規則構件很難實現均勻的粒子分布，很少有文獻專門提出粒子的初始化問題。通常借助于有限元網格來生成粒子，但卻很難實現粒子的均勻化分布，從而影響粒子的初始物理量，特別是光滑長度的影響甚大，直接關系到計算精度。本文提出利用通用軟件實現粒子初始化的方法。

具體的思想是：首先利用通用有限元軟件建模，然后進行網格劃分，單元邊長最好與預定的粒子間距相當，通常生成四面體（或三角形）網格，然后提取出每個單元頂點坐標存儲到相應文件，在SPH程序中讀入各個頂點的坐標，然后計算出每個單元形心的坐標位置作為SPH粒子，而每個四面體的質量、密度就是每個粒子的質量、密度。下面就是粒子的初始光滑長度的選取。

Liu等［3］提出如下公式計算初始光滑長度：

SPH方法的問題域以粒子形式存在，每個粒子攜帶著諸多物理量，如何直觀明了的展示計算結果，后處理十分關鍵。通用軟件EnSight提供了SPH方法后處理功能，將SPH方法的計算結果輸出為特定格式文件組，由.case文件讀入，繼而可以進行可視化處理。

2 理論應用實踐

2.1 數值模型

圖3 數值計算模型Fig.3 Numerical model

本文所用的SPH數值模型如圖3所示，其中鋼板厚度為t，水域半徑為R，裝藥半徑為r。球形裝藥在半圓形水域引爆，鋼板應用固體力學的Mie-Gruneisen 方程［6］及 Johnson-cook 材 料屈 服 模型［7-8］，炸藥 TNT 采用 Jones.Wilkins-Lee（JWL）狀態方程，水采用 Mie-Gruneisen 狀態方程［3］。

2.2 沖擊波的產生與傳播

當水中炸藥中心起爆后，產生一系列的爆轟波和沖擊波分別在炸藥中和水中和鋼板中傳播，以表4中參數為例，來說明水下接觸爆炸沖擊波的產生與傳播。

由圖4可以看到，當炸藥中心起爆后，由爆轟產物開始膨脹壓縮周圍的物體，并且傳出爆轟波。在水中和鋼板中開始了球面沖擊波的傳播，此時壓力高達幾GPa。在105 μs時鋼板中沖擊波傳播速度明顯高于在水中的傳播速度，沖擊波的壓力峰值由于水和鋼板的作用開始進行衰減，此時鋼板開始出現凹陷。直到230 μs時，在水域中的壓力發生了明顯的變化，在鋼板的法線方向上壓力大于同一半徑處的壓力，這是由于炸藥產生的爆轟波傳向鋼板方向后，鋼板的阻抗較大，很快就反射了沖擊波，追趕由炸藥直接傳入水中的沖擊波，經過疊加作用，出現了較大的壓力。此時鋼板中沖擊波已經傳至邊界，且已經開始斷裂，出現圓盤化現象。在255 μs時，鋼板中反射的稀疏波與傳出的入射波疊加形成了低壓區。由于稀疏波傳播的較快，在330 μs時鋼板中幾乎不存在沖擊波了，而水中的沖擊波壓力也逐漸地降低，近鋼板邊界處由于在鋼板面反射了一系列的稀疏波，形成了低壓區。在455 μs時鋼板已經出現了明顯的沖塞破口。

圖4 典型時刻沖擊波壓力分布圖Fig.4 Shock wave pressure distribution of typical time

表4 數值模型參數Tab.4 Numerical model parameters

2.3 鋼板破口分析

圖5為鋼板中沖擊波的傳播及其破壞過程。鋼板的接觸爆炸破壞一般分為三個過程：沖塞凹陷，圓盤化破口和花瓣開裂翻轉［9-10］。在鋼板的沖塞凹陷階段，鋼板的隆起范圍與初始載荷的作用面相當［11］，圖中255 μs時凹陷半徑與藥包半徑0.05 m相當，這與圖6是一致的。在255 μs時鋼板開始出現裂口，進入圓盤化破口階段。此后破口與沖塞物開始橫向運動，在破口的環向拉伸和彎曲以及徑向拉伸，消耗初始沖塞后的動能，破口逐漸擴大，如圖280 μs至480 μs的過程。當圓盤化破口的環向拉伸應變時，破口邊緣就將發生斷裂開始發生花瓣開裂，進而翻轉［12］，本文不予研究。

2.4 不均勻粒子分布對數值模擬的影響

SPH方法在進行數值模擬時，為了保證模擬小尺度構件和不規則物體的可行性和數值連續性及其穩定性，從客觀上和主觀上都避免不了粒子的不均勻化分布。因此，粒子不均勻分布對計算結果的影響就十分值得研究。

數值模型中，圓形裝藥半徑r和半圓形水域的半徑R分別為0.05 m和1 m，粒子間距均為0.005 m；鋼板厚度t為0.03 m，粒子間距為變量，分別取為0.005 m，0.004 m，0.003 m，0.002 m，0.001 m，從而實現粒子的不同間距比例。

如圖7所示，選擇炸藥正下方0.3 m處沖擊波的壓力時程曲線。炸藥爆轟后，水中產生的沖擊波開始向外傳播，很快就達到了峰值，隨后以指數形式開始衰減，壓力時程曲線與Cole經驗公式［1］所得形式基本吻合，由于存在鋼板的反射作用，在85 μs時反射的沖擊波與直接由炸藥傳出的沖擊波疊加，出現了第二次峰值。隨后繼續衰減。改變鋼板的粒子間距，當水域與鋼板的粒子間距比為1∶1時，沖擊波的壓力時程曲線較為平穩。隨著粒子的間距比增加，由炸藥爆轟后直接傳出的沖擊波在達到壓力峰值前與等間距粒子分布時相同，但當由鋼板反射的沖擊波與炸藥直接傳出的沖擊波疊加后，不同粒子間距比彼此之間出現了較大的不同。當粒子的間距比為5∶4時，壓力時程曲線與1∶1時非常相近；當粒子間距比為5∶3時，開始出現明顯的數值波動，當粒子的間距比為5∶2時，開始出現持續的波動；當粒子間距比為5∶1時，計算崩潰。

圖7 不同粒子間距比的沖擊波壓力時程曲線Fig.7 The shock wave pressure-time curve of different particle spacing ratio

3 結論

本文基于傳統SPH算法，針對水下爆炸過程中經常遇到的問題，進行如下改進：

（1）對水下爆炸問題中經常出現的大密度比問題，總結了基于體積思想的SPH控制方程，通過實例驗證，證明了其解決水下接觸爆炸問題的有效性；

（2）SPH方法解決水下爆炸問題存在問題域大，計算效率低等問題，為了保證計算精度又提高計算效率，提出了變光滑長度的鏈表搜索算法，同時將其應用到水下接觸爆炸問題，極大地提高了計算效率。粒子數的越多，計算效率的提高越明顯；

（3）對于水下接觸爆炸問題中，經常出現小尺度、不規則構件，在保證計算精度的前提下，提出了粒子初始化方法和后處理辦法；

（4）通過對二維水下接觸爆炸問題的應用實踐，水下爆炸沖擊波的產生和傳播符合沖擊波理論，爆炸破口在沖塞凹陷和圓盤化階段形式和大小與既有文獻論述一致，驗證了所提出理論的有效性和可行性。最后給出了不同粒子間距比對數值模擬的影響，在小于2時數值結果相對穩定。

［1］ Cole P.水下爆炸［M］.北京：國防工業出版社，1960：1-250.

［2］Ott F，Schnetter E.A modified SPH approach for fluids with large density differences［J］.Physics，2003，1：1-11.

［3］ Liu G R，Liu M B.Smoothed particle hydrodynamics：a meshfree particle method［M］.German：Springer Berlin/Heidelberg，2004：1-491.

［4］張姝慧，汪繼文.SPH法中初始時粒子配置的分析［J］.計算機技術與發展，2007，17（6）：36-38.

［5］倪國喜，林 忠.無網格方法中粒子分布與自適應研究［J］.計算物理，2006，26（4）：419-423.

［6］韓 旭，伍乾坤，胡德安.SPH方法模擬高速碰撞問題的實例［A］.計算爆炸力學進展［C］.2006.

［7］強洪夫，高巍然，傅學金.考慮拉伸不穩定性的圓柱形鋼彈丸穿透中厚鋼板的SPH三維數值模擬［A］.第17屆全國結構工程學術會議論文集［C］.2008.

［8］ Libersky L D，Petscheck G，Carney T C，et al.High strain lagrangian hydrodynamicsa three dimensional SPH code for dynamic material response［J］.Journal of Computational Physics，1993，109：67-75.

［9］ Lee Y W，Wierzbicki T.Fracture prediction of thin plates under localized impulsive loading.Part I：dishing［J］.International Journal of Impact Engineering，2005，31：1253-1276.

［10］ Lee Y W，Wierzbicki T.Fracture prediction of thin plates under localized impulsive loading.PartⅡ：discing and petalling［J］.International Journal of Impact Engineering，2005，31：1277-1208.

［11］ Rajendran R，Narasimhan K.Damage prediction of clamped circular plates subjected to contact underwater explosion［J］.International Journal of Impact Engineering，2001，25：373-386.

［12］張振華.剛塑性板在柱狀炸藥接觸爆炸載荷作用下的花瓣開裂研究［J］.船舶力學，2004，8（5）：113-119.