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

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

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

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

1 理論研究

1.1 SPH控制方程

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

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

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

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

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

表2 數(shù)值模型參數(shù)Tab.2 Numerical model parameters

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

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

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

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

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

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

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

1.3 SPH方法的前、后期處理

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

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

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

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

2 理論應(yīng)用實踐

2.1 數(shù)值模型

圖3 數(shù)值計算模型Fig.3 Numerical model

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

2.2 沖擊波的產(chǎn)生與傳播

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

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

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

表4 數(shù)值模型參數(shù)Tab.4 Numerical model parameters

2.3 鋼板破口分析

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

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

SPH方法在進(jìn)行數(shù)值模擬時，為了保證模擬小尺度構(gòu)件和不規(guī)則物體的可行性和數(shù)值連續(xù)性及其穩(wěn)定性，從客觀上和主觀上都避免不了粒子的不均勻化分布。因此，粒子不均勻分布對計算結(jié)果的影響就十分值得研究。

數(shù)值模型中，圓形裝藥半徑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，從而實現(xiàn)粒子的不同間距比例。

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

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

3 結(jié)論

本文基于傳統(tǒng)SPH算法，針對水下爆炸過程中經(jīng)常遇到的問題，進(jìn)行如下改進(jìn)：

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

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

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

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

［1］ Cole P.水下爆炸［M］.北京：國防工業(yè)出版社，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］.計算機技術(shù)與發(fā)展，2007，17（6）：36-38.

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

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

［7］強洪夫，高巍然，傅學(xué)金.考慮拉伸不穩(wěn)定性的圓柱形鋼彈丸穿透中厚鋼板的SPH三維數(shù)值模擬［A］.第17屆全國結(jié)構(gòu)工程學(xué)術(shù)會議論文集［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］.船舶力學(xué)，2004，8（5）：113-119.