基于SPH-FEM算法的鋼板接觸爆破數(shù)值模擬研究

陳 誠，詹發(fā)民，周方毅，姜 濤

(海軍潛艇學(xué)院，山東 青島 266033)

1 引言

鋼結(jié)構(gòu)的接觸爆破在工程爆破領(lǐng)域有著非常廣泛的應(yīng)用。由于爆破過程中產(chǎn)生的大變形、物質(zhì)流動、界面接觸等問題較為復(fù)雜，使用傳統(tǒng)有限元法對此類問題進(jìn)行數(shù)值模擬存在一定困難。目前對于鋼結(jié)構(gòu)接觸爆破的數(shù)值模擬計算，通常采用ALE算法進(jìn)行數(shù)值模擬計算，其優(yōu)點是能夠克服爆破過程中帶來的網(wǎng)格畸變問題，計算準(zhǔn)確度較高，能夠較為準(zhǔn)確地模擬炸藥對鋼板的破壞過程，ALE方法的缺點在于需要建立較大的空氣網(wǎng)格，計算量較大，網(wǎng)格產(chǎn)生的大變形容易導(dǎo)致計算中止等。SPH算法是1977年由Lucy和Monaghan等人提出的一種無網(wǎng)格法，最初是用于解決天體運動問題，近年來被廣泛用于爆炸效應(yīng)的數(shù)值模擬計算。但由于SPH粒子的光滑長度一般取值較小，因此計算效率不高，且容易出現(xiàn)物質(zhì)穿透、計算崩潰等問題。

為了綜合FEM法與SPH算法的優(yōu)點，提高計算效率及準(zhǔn)確性，本文采用SPH-FEM耦合算法對鋼板接觸爆破過程進(jìn)行數(shù)值模擬計算。基本方法是在炸藥與鋼板接觸面附近變形較大的區(qū)域使用SPH粒子建模，鋼板外圍變形較小的區(qū)域采用FEM單元網(wǎng)格進(jìn)行建模。這樣既克服了FEM法導(dǎo)致的網(wǎng)格大變形問題，又能減少SPH粒子數(shù)量，提高計算的效率。

2 SPH算法基本理論

2.1 控制方程

SPH算法本質(zhì)上是一種Lagrange法，因此其控制方程采用Lagrange形式的Navier-Stockes方程，通過三大守恒定律得到的

(1)

其中，，=1，2，3，代表坐標(biāo)方向，遵循張量運算字母指標(biāo)法；、、、、、、、表示密度、內(nèi)能、速度分量、空間坐標(biāo)、總應(yīng)力張量和時間；表示作用在單位質(zhì)量流體上的體力。

引入SPH的核近似和粒子近似原理，對式(1)進(jìn)行粒子近似，得到用于爆炸數(shù)值模擬的離散化SPH控制方程

(2)

2.2 SPH粒子與FEM網(wǎng)格的耦合

對于SPH粒子與FEM網(wǎng)格的耦合，目前主流的方法有虛粒子法和罰函數(shù)法兩種。虛粒子法雖然具有較高的精度，但是于對于復(fù)雜邊界生成虛粒子困難，在某些情況下可能會發(fā)生粒子穿透邊界；而罰函數(shù)方法容易實現(xiàn)，且能夠處理復(fù)雜幾何邊界條件。

3 數(shù)值模擬計算模型

3.1 計算模型

數(shù)值模擬計算模型由炸藥、離散為有限元單元網(wǎng)格的鋼板和離散為SPH粒子的鋼板三部分組成。炸藥及鋼板接觸附近區(qū)域使用SPH粒子，外圍鋼板使用有限元網(wǎng)格。如圖1所示，鋼板尺寸為長40cm寬40cm厚度1cm；鋼板中心20cm×20cm的區(qū)域離散為SPH粒子，粒子間隔0.5cm；鋼板外圍區(qū)域離散為尺寸為1cm的六面體有限元單元；200gTNT炸藥尺寸為長10cm寬5cm厚度2.5cm，離散為SPH粒子，粒子間隔0.5cm。SPH粒子與有限元網(wǎng)格使用LS-DYNA中的*CONTACT-TIED-NODE-TO-SURFACE-OFFSET關(guān)鍵字進(jìn)行耦合。

圖1 SPH粒子與FEM網(wǎng)格的劃分

3.2 材料模型及狀態(tài)方程

炸藥選用LS-DYNA材料庫中的*MAT-HIGH-EXPLOSIVE-BURN材料進(jìn)行模擬，炸藥密度ρ=1.60g/cm，炸藥爆速D=6930m/s，C-J壓力P=27GPa。爆轟產(chǎn)物通過式(3)JWL狀態(tài)方程進(jìn)行描述，式(3)中、、、、為與炸藥性質(zhì)有關(guān)的常數(shù)。本例中取=374，=733，=415，=095，=030。

=(1-)-+(1-)-+

(3)

鋼板采用*MAT＿JOHNSON＿COOK材料模型，鋼板相關(guān)參數(shù)見表1，通過式(4)Gruneisen狀態(tài)方程進(jìn)行描述：

表1 鋼板JOHNSON＿COOK模型相關(guān)參數(shù)

(+)

(4)

式(4)中：μ-μ曲線的截距C=1.4569；μ-μ曲線斜率的系數(shù)S=1.49、S=-0、S=0；Gruneisen參數(shù)γ=2.17；一階體積修正量a=0.46；μ=(ρ/ρ)-1。

3.3 計算結(jié)果及分析

因炸藥20后爆轟完全，對計算影響較少，為避免炸藥的粒子出現(xiàn)速度超范圍的情況，設(shè)置炸藥粒子失效時間為20，設(shè)置總體計算時間為120。計算完成后使用后處理軟件查看計算結(jié)果。

圖2為20時的計算結(jié)果，從圖中能夠明顯看出炸藥爆轟過程和爆轟產(chǎn)物的飛散，以及鋼板產(chǎn)生形變的過程。數(shù)值模擬計算結(jié)果顯示爆炸發(fā)生后，爆轟產(chǎn)物與鋼板接觸后，鋼板首先發(fā)生凹陷，當(dāng)達(dá)到屈服極限后，鋼板背面發(fā)生脫落。圖3為120時的計算結(jié)果，可見鋼板中心15×8的區(qū)域范圍內(nèi)產(chǎn)生凹痕，12×6區(qū)域范圍內(nèi)的產(chǎn)生脫落穿孔。

圖2 20μs時炸藥爆轟及鋼板變形效果

圖3 120μs時鋼板的凹痕及穿孔

為了驗證粒子與網(wǎng)格的耦合效果，選取接觸面附近相鄰的一個粒子和單元(本文選取2512粒子和10864單元)，導(dǎo)出其軸方向上的速度時程曲線如圖4所示。對比壓力曲線發(fā)現(xiàn)其波形及峰值大致相同，說明該處粒子將壓力等物理量傳遞給了網(wǎng)格，耦合效果較好，進(jìn)一步說明了-耦合算法的有效性。

圖4 SPH粒子與FEM網(wǎng)格z軸方向速度對比

3.4 SPH-FEM耦合算法與FEM算法計算效果對比

為了比較-耦合算法與單純算法的優(yōu)劣性，使用純法建立有限元模型，通過算法進(jìn)行求解。有限元模型與31節(jié)中計算模型尺寸相同，在鋼板內(nèi)外建立30厚空氣層。炸藥及空氣使用網(wǎng)格，鋼板使用網(wǎng)格。材料模型及本構(gòu)方程與32節(jié)中相同。圖5為兩種算法計算效果的對比，可以看出，兩種算法都能模擬出炸藥對鋼板的破壞過程。在鋼板相同位置選取觀測點，導(dǎo)出其壓力時程曲線，本文選取鋼板上坐標(biāo)為(15，0，0)的點為觀測點，圖6為兩個觀測點的壓力時程曲線，從曲線可以得出-耦合算法得到的觀測點峰值壓力為0151，而純法得到的觀測點峰值壓力為0143，誤差在5左右，計算結(jié)果較為接近。從計算時間上看，建立單純網(wǎng)格，使用算法求解的計算時間為687秒，而-耦合算法計算時間僅為11秒，由此可見-耦合算法大大節(jié)省了算力，提高了計算效率。

圖5 120μs時兩種算法計算結(jié)果

圖6 鋼板相同位置單元的壓力時程曲線

4 試驗驗證

4.1 試驗工況布置

為了驗證-耦合算法的準(zhǔn)確性，設(shè)計全尺寸接觸爆破試驗。如圖5所示，鋼板選用235普通碳素結(jié)構(gòu)鋼，平放于地面上。200藥塊置于鋼板中心，用膠帶進(jìn)行固定，鋼板背面設(shè)置厚度為30的臨空面，使用8號雷管起爆，起爆點位于藥塊左端中心的雷管室。

圖7 裝藥設(shè)置

4.2 試驗效果及分析

爆破后效果如圖8所示，在爆炸作用下，鋼板產(chǎn)生穿孔，邊緣產(chǎn)生凹痕。經(jīng)測量，鋼板中心產(chǎn)生10×5的穿孔，12×8區(qū)域產(chǎn)生明顯的凹痕，與數(shù)值模擬計算結(jié)果吻合較好，驗證了該算法的準(zhǔn)確性。

圖8 爆破效果與數(shù)值模擬結(jié)果對比

5 結(jié)論

-耦合算法能夠基本模擬鋼板接觸爆破中的炸藥爆轟、壓力傳遞，鋼板形變等過程，計算準(zhǔn)確度較高，與單純法相比，-耦合算法因無需劃分復(fù)雜網(wǎng)格，大大減少了計算復(fù)雜度，提高了計算效率。-耦合算法實施的關(guān)鍵點是粒子與網(wǎng)格的耦合，本文采取罰函數(shù)法進(jìn)行處理，對粒子和單元接觸面進(jìn)行接觸定義。通過試驗驗證，-耦合算法能夠成功應(yīng)用在鋼板接觸爆破中。