SPH方法在聚能裝藥射流三維數(shù)值模擬中的應(yīng)用＊

李 磊，沈兆武，李學(xué)嶺，倪小軍

（1.中國科學(xué)技術(shù)大學(xué)近代力學(xué)系，安徽 合肥 230027；2.安徽方圓機(jī)電股份有限公司，安徽 蚌埠 233010）

聚能裝藥技術(shù)廣泛應(yīng)用于彈藥武器設(shè)計(jì)和工業(yè)生產(chǎn)行業(yè)，例如破甲彈、聚能切割、石油開采等領(lǐng)域。聚能裝藥起爆后形成聚能射流過程的數(shù)值模擬對聚能裝藥技術(shù)的研究及設(shè)計(jì)具有重要的意義。聚能射流形成過程是彈塑性流體動力學(xué)問題，對于整個(gè)過程進(jìn)行數(shù)值模擬是一個(gè)難點(diǎn)。當(dāng)采用Lagrange有限元法模擬射流形成過程時(shí)，必然遇到網(wǎng)格發(fā)生大畸變和滑移面處理等一系列問題，最終導(dǎo)致計(jì)算精度降低甚至計(jì)算終斷。Euler方法不存在網(wǎng)格發(fā)生大畸變問題，但難以準(zhǔn)確描述各類界面。光滑粒子動力學(xué)（Smoothed particle hydrodynamics，簡稱SPH）方法是一種Lagrange無網(wǎng)格數(shù)值方法，它采用帶質(zhì)量、動量、能量的粒子構(gòu)成離散計(jì)算域，不同材料的粒子自然地構(gòu)成界面，材料間的相互作用可以由粒子間的相互作用來自然地模擬，材料的變形不依賴于網(wǎng)格而通過粒子的運(yùn)動來描述，因此在理論上，SPH方法能夠“自然地”模擬侵徹貫穿、高速碰撞等物理現(xiàn)象。SPH方法最早由L.B.Lucy［1］和R.A.Gingold等［2］于1977年提出并應(yīng)用于天體物理學(xué)，后來被廣泛地應(yīng)用于力學(xué)研究的各個(gè)領(lǐng)域，J.W.Swegle等［3］首先用SPH方法模擬爆炸問題，M.B.Liu［4－6］等應(yīng)用SPH方法模擬了聚能裝藥的爆轟過程。

本文中，首先介紹SPH方法的基本理論，其次采用有限元程序LS－DYNA中的歸一化粒子近似公式對聚能裝藥射流形成過程進(jìn)行三維數(shù)值模擬，并對射流形成過程和機(jī)理進(jìn)行分析；然后將SPH數(shù)值模擬結(jié)果與二維有限元數(shù)值模擬（FEM方法）結(jié)果及理論計(jì)算結(jié)果相比較，討論SPH方法的精確性、有效性以及相對于FEM方法的優(yōu)越性；最后分析和比較不同初始相鄰粒子數(shù)對于計(jì)算結(jié)果的影響。

1 光滑粒子法

SPH方法是通過構(gòu)造一個(gè)核函數(shù)對離散質(zhì)點(diǎn)位置進(jìn)行核估計(jì)來計(jì)算梯度的相關(guān)項(xiàng)，因此不需要網(wǎng)格來求解偏微分的差分，只需將微分形式的守恒方程轉(zhuǎn)化為積分方程形式，計(jì)算出在任意一點(diǎn)上的各個(gè)場變量的核估計(jì)。函數(shù)f（x）在空間某一點(diǎn)x上的核估計(jì)可以通過函數(shù)f（x）在域Ω中的積分獲得

式中：W（x－x′，h）為核函數(shù)或權(quán)函數(shù)，x為估計(jì)點(diǎn)的空間坐標(biāo)，x′為對估計(jì)點(diǎn)有貢獻(xiàn)作用的空間點(diǎn)坐標(biāo)，h為緊支域的度量尺寸即光滑長度。

函數(shù)導(dǎo)數(shù)的核估計(jì)可通過將式（1）里的函數(shù)f（x）視為導(dǎo)數(shù)?f（x）／?x，并利用分步積分和核函數(shù)W在積分域Ω邊界上為零的條件求得

由此可見光滑粒子法的基本思想之一是將函數(shù)導(dǎo)數(shù)的核估計(jì)轉(zhuǎn)換成核函數(shù)的導(dǎo)數(shù)，核函數(shù)是預(yù)先設(shè)定的已知函數(shù)。將解域Ω劃分為M個(gè)子域粒子，每個(gè)子域粒子j的質(zhì)量和密度分別為mj＝m（xj），ρj＝ρ（xj），設(shè)f（x）在粒子i、j上的值分別為fi＝f（xi），fj＝f（xj），則f（x）及其導(dǎo)數(shù)在粒子i上的核估計(jì)式（1）～（2）的離散式為

式中：Wij＝W（xj－xi，h）為離散核函數(shù)，積分微體元dx′＝mj／ρj，α表示空間維序數(shù)，下標(biāo)j為編號為i粒子的臨近粒子編號。

考慮人工粘性Πij和人工熱流Hi的影響，具有材料強(qiáng)度的全應(yīng)力張量空間中的SPH插值公式可表示為［7－8］

式中：u、σ分別表示速度和應(yīng)力張量，上標(biāo)α、β表示空間方向；E為比內(nèi)能；Γ、ε分別為偏應(yīng)力和應(yīng)變張量為狄拉克函數(shù)，當(dāng)α＝β時(shí)為1，當(dāng)α≠β時(shí)為0。

核函數(shù)是SPH方法中核心組成部分，常用的有B－spline核函數(shù)、Gauss核函數(shù)、二次核函數(shù)及指數(shù)核函數(shù)等，其中B－spline核函數(shù)是目前應(yīng)用最廣泛的核函數(shù)，其形式為［8］

2 數(shù)值模擬與分析

2.1 計(jì)算模型

圖1 四分之一聚能裝藥模型（剖面圖）Fig.1The quarter model of the shaped charge（profile map）

首先利用ANSYS有限元軟件建立聚能裝藥三維模型并劃分為Lagrange單元網(wǎng)格，根據(jù)對稱性只建立四分之一模型。其次，生成PART以及最終求解文件（K文件）。然后，修改K文件：刪除Lagrange單元信息，添加SPH質(zhì)點(diǎn)及其相關(guān)屬性設(shè)置，圖1所示為修改K文件之后得到的光滑粒子計(jì)算模型。最后，利用LS＿DYNA 971求解器求解K文件。然而，采用這種方式建立SPH模型時(shí)需要注意實(shí)際長度和建模長度的區(qū)別［9］。

模型中炸藥配置了151 827個(gè)粒子，藥型罩配置了23 816個(gè)粒子。粒子間距為0.05～0.06cm，初始光滑長度為0.08cm。炸藥粒子的最小、最大光滑長度分別為0.5倍和2.0倍初始光滑長度。藥型罩粒子的最小、最大光滑長度分別為0.01倍和1.20倍初始光滑長度。每個(gè)粒子初始相鄰粒子數(shù)設(shè)為1 200，時(shí)間步長系數(shù)為0.015。

炸藥采用 MAT＿HIGH＿EXPLOSIVE＿BURN高能燃燒模型，狀態(tài)方程為JWL［9－10］

式中：p為爆轟壓力，E為炸藥的比內(nèi)能，V為相對比容，狀態(tài)方程參數(shù)A＝371.2GPa，B＝3.21GPa，R1＝4.15，R2＝0.95，ω＝0.3，炸藥密度ρ＝1.60kg／m3，C－J爆速vD＝7.74km／s。

藥型罩材料為紫銅，采用 MAT＿JOHNSON＿COOK材料模型，狀態(tài)方程為 Mie－Grüneisen［9－10］

式中：vs－vp曲線截距C＝0.394；vs－vp曲線斜率的相關(guān)因數(shù)S1＝1.49，S2＝0，S3＝0，vs為沖擊波速度，vp為質(zhì)點(diǎn)速度；Grüneisen因數(shù)γ0＝2.02，Grüneisen因數(shù)的一階體積修正因數(shù)α′＝0.47，E′為材料單位體積內(nèi)能，體積變化率分別為紫銅的瞬時(shí)和初始密度，＝8.93kg／m3，剪切模量γ＝47.7GPa，厚度d0＝0.18cm，外直徑D＝4.0cm，頂角ξ＝53°，熔化溫度T＝1 360K。

2.2 射流形成過程和機(jī)理分析

以藥型罩內(nèi)、外壁上編號分別為154 202和164 459的節(jié)點(diǎn)（粒子）為例，在不同時(shí)刻粒子速度大小的SPH數(shù)值模擬結(jié)果與FEM方法模擬結(jié)果的對比如圖2～3中曲線所示。由圖2～3可知，這2種方法的數(shù)值模擬結(jié)果十分接近。下面以SPH數(shù)值模擬結(jié)果為例簡單地分析射流形成過程和機(jī)理。

圖2 粒子（154 202）速度與FEM方法計(jì)算結(jié)果比較圖Fig.2The velocity of the particle（154 202）calculated by SPH （3D）and finite element method（2D）

圖3 粒子（164 459）速度與FEM方法計(jì)算結(jié)果比較圖Fig.3The velocity of the particle（164 459）calculated by SPH and finite element method（2D）

當(dāng)炸藥起爆后產(chǎn)生球形爆轟波，爆轟波首先到達(dá)藥型罩的頂端，當(dāng)爆轟波前沿的沖擊波傳入藥型罩后，驅(qū)使藥型罩向軸線方向閉合，并在軸線處產(chǎn)生高壓、高溫、高速射流以及速度相對較低的杵體。粒子154 202于t＝5μs時(shí)刻開始向軸線方向閉合，由于該粒子在向軸線閉合過程中，不可避免地與其他軸線方向運(yùn)動的粒子相互碰撞擠壓，所以該粒子加速度是變化的，碰撞后粒子與罩內(nèi)壁的其他粒子形成能量密度很高的射流。由圖2中粒子154 202的速度曲線可以看出，當(dāng)t＝8.8μs時(shí)刻，沖擊波將越過粒子，其加速度逐漸降為0，而速度達(dá)到最大值為5 448m／s，然后在應(yīng)力波的作用下粒子的速度發(fā)生小幅度震蕩，并以約5 210m／s的平均速度沿著軸線方向運(yùn)動。由圖3可知，粒子164 459在t＝6.6μs時(shí)刻向軸線方向閉合運(yùn)動，并與相鄰粒子相互碰撞和擠壓變形；當(dāng)t＝12.0μs時(shí)速度達(dá)到最大值為1 184m／s，之后粒子速度逐漸降低并繼續(xù)向軸線方向閉合運(yùn)動；當(dāng)t＝25.0μs時(shí)閉合過程完成，與藥型罩外壁的其它粒子形成能量密度低的杵體，并以相對穩(wěn)定的速度（約429m／s）跟隨射流向前運(yùn)動。由于射流頭部至尾部存在速度梯度，且射流頭部速度遠(yuǎn)大于杵體速度，所以射流在運(yùn)動過程中不斷拉伸，并在拉伸過程中產(chǎn)生頸縮與斷裂，同時(shí)與杵體逐漸分離。當(dāng)采用SPH法計(jì)算時(shí)，計(jì)算過程十分穩(wěn)定，不存在網(wǎng)格畸變問題。而當(dāng)采用Lagrange有限元方法計(jì)算時(shí)，在t＝100.0μs時(shí)刻網(wǎng)格發(fā)生了嚴(yán)重的扭曲纏繞變形，并導(dǎo)致計(jì)算中斷。此外，采用FEM方法計(jì)算多介質(zhì)問題時(shí)，常使用接觸算法，然而該算法中參數(shù)較多，很難把握，經(jīng)常出現(xiàn)穿透現(xiàn)象，如圖4（b）所示，進(jìn)而導(dǎo)致計(jì)算終斷，而采用SPH方法計(jì)算多介質(zhì)問題時(shí)，2種不同材料之間同樣計(jì)算質(zhì)點(diǎn)近似，交界面處不需要特別的處理，所以穿透現(xiàn)象很少發(fā)生，如圖4（a）所示。因此，SPH方法比FEM方法更適合模擬此類的多介質(zhì)、大變形問題。

2.3 數(shù)值模擬結(jié)果比較與分析

圖4 藥型罩頂部與炸藥交界面處數(shù)值模擬圖Fig.4The simulation images in the interface between the liner’s top and the main charge

2.3.1 射流頭部速度和長度的理論計(jì)算

根據(jù)PER理論可以計(jì)算藥型罩射流頭部速度［11］

式中：vj為射流頭部速度，θ為藥型罩半錐角，φ、φ＋分別為動態(tài)壓合角、穩(wěn)態(tài)壓合角，x為罩微元在軸線方向坐標(biāo)（罩頂點(diǎn)為原點(diǎn)），v為指數(shù)壓合速度，vx′為v對x的偏導(dǎo)數(shù)，t0為罩微元初始壓垮時(shí)刻，t為罩微元壓垮至軸線的時(shí)刻，v0為常數(shù)壓合速度，τ為指數(shù)加速模型中的特征時(shí)間常數(shù)［12］

式中：m為單位面積藥型罩的初始質(zhì)量，pCJ＝28GPa，為C－J爆轟壓力，c1、c2為經(jīng)驗(yàn)參數(shù)。

對v0采用格尼模型計(jì)算［12］

式中：Q＝1.5（3χ2＋4χ＋1）／（4χ2＋4χ＋1），χ＝R0／Rr，μ＝mL／me，κ≈C0pCJτ（R0－Rr）［（R0／Rr）－1］，Rr、R0分別為罩微元x處炸藥內(nèi)表面和外表面半徑，me、mL分別藥型罩微元x處的炸藥和罩的質(zhì)量，C0為經(jīng)驗(yàn)參數(shù)，E為炸藥比內(nèi)能。

射流長度的計(jì)算采用P.C.Chou等［11］提出的公式

式中：r（x，t）為坐標(biāo)為x的罩微元在t時(shí)刻的位置，ζ為拋射角，t′為x罩微元?jiǎng)偟竭_(dá)軸線上的時(shí)間，L為罩頂部炸藥高度，η（x，t）為x罩微元在t時(shí)刻的軸向應(yīng)變。

2.3.2 SPH數(shù)值模擬結(jié)果、FEM方法模擬結(jié)果及理論計(jì)算結(jié)果的比較和分析

［11－12］，將x＝1.076 9cm處，θ＝26.5°，R0＝1.345 1cm，Rr＝0.560 6，μ＝0.629，m＝9.26g／cm2，me＝1.196g／cm2，φ＋＝43.3°［13］，φ＝55°［13］，L＝3cm 等計(jì)算參數(shù)代入式（11）～（18），計(jì)算射流頭部速度vj（即粒子154 202速度）和射流長度lj，vj的理論計(jì)算結(jié)果為5 332m／s，F(xiàn)EM 方法計(jì)算結(jié)果為5 197m／s，SPH 方法的計(jì)算結(jié)果為5 210m／s，可知分別采用SPH方法和FEM方法計(jì)算的射流頭部速度的誤差分別為2.29%、2.53%；分別計(jì)算的射流長度的比較結(jié)果，如表1所示，2種方法計(jì)算的射流長度的平均誤差分別為0.143%、4.425%，與理論計(jì)算結(jié)果的誤差均小于5%，因此可認(rèn)為SPH方法模擬聚能裝藥問題是可行的且具有更好的精度。

表1 不同時(shí)刻射流長度計(jì)算值Table 1Length of the jet at different time calculated by different method

2.4 光滑長度和初始相鄰光滑粒子數(shù)對計(jì)算結(jié)果的影響分析

光滑長度h對于計(jì)算精度和效率至關(guān)重要。在聚能裝藥爆炸形成聚能射流過程中，炸藥及藥型罩物質(zhì)變形速度和變形量很大。若采用固定光滑長度，在炸藥爆轟膨脹時(shí)，核函數(shù)影響域內(nèi)的粒子數(shù)量會不斷減少，可能出現(xiàn)數(shù)值斷裂等問題，而當(dāng)藥型罩在爆轟波作用下產(chǎn)生壓縮變形時(shí)，核函數(shù)影響域內(nèi)的粒子數(shù)會急劇增加，這將顯著地降低計(jì)算效率，而且可能產(chǎn)生壓縮不穩(wěn)定性［14］。因此采用變光滑長度的時(shí)間積分格式［10］

式中：h（t）為與時(shí)間相關(guān)的光滑長度，d表示空間維度，vp為粒子速度。

圖5 初始相鄰粒子數(shù)為200和300時(shí)藥型罩塑性應(yīng)變云圖Fig.5Plastic－strain contour of the liner when Nis 200or 300

初始相鄰光滑粒子數(shù)量N是由初始光滑長度h決定的，由于初始時(shí)刻粒子間距的非均勻性，導(dǎo)致每個(gè)粒子的緊支域內(nèi)的初始相鄰粒子數(shù)不相同，即每個(gè)粒子的緊支域內(nèi)包含不同的質(zhì)量，這將造成非一致的核近似精度。因此在采用LS＿DYNA程序模擬時(shí)，需要設(shè)定一個(gè)一致的初始相鄰粒子數(shù)量N用來調(diào)整初始化階段的內(nèi)存分配。在計(jì)算過程中發(fā)現(xiàn)變量N的設(shè)置除了影響計(jì)算效率外，對于計(jì)算精度也會產(chǎn)生一定影響。從能量角度分析并比較不同初始相鄰粒子數(shù)量（N＝200～3 500）對計(jì)算結(jié)果的影響，計(jì)算模型與2.1節(jié)相同。當(dāng)N設(shè)為200和300時(shí)，藥型罩的邊界自由面出現(xiàn)了“粒子飛散”的現(xiàn)象，如圖5所示，而且在計(jì)算十幾微秒之后，時(shí)間步長逐漸減小，最終導(dǎo)致計(jì)算困難。這可能是由于初始相鄰粒子數(shù)量太少，藥型罩邊緣粒子的球形緊支域內(nèi)沒有足夠的其他近似粒子，即光滑函數(shù)在該粒子的緊支域內(nèi)呈現(xiàn)出嚴(yán)重的非正則化性質(zhì)，進(jìn)而影響了計(jì)算精度。

為了避免以上出現(xiàn)的問題，增加初始相鄰粒子數(shù)量，當(dāng)N≥400時(shí)，“粒子飛散”的現(xiàn)象不再發(fā)生，射流成型質(zhì)量較好，計(jì)算過程穩(wěn)定。藥型罩的內(nèi)能e和動能ek時(shí)程曲線可知如圖6～7所示，炸藥爆轟能量主要以動能形式傳遞給藥性罩，而藥型罩的內(nèi)能增加很小。當(dāng)N取不同值時(shí)，計(jì)算得到的藥型罩內(nèi)能和動能在任意時(shí)刻的最大誤差分別小于5%和1%，因此當(dāng)N＝400～3 500時(shí)，計(jì)算結(jié)果之間的誤差較小。將N取不同值時(shí)計(jì)算得到的射流頭部速度vj與理論計(jì)算結(jié)果相比較如圖8所示，從圖中可以看出，當(dāng)N≤300時(shí)誤差最大，大于7.5% ；當(dāng)400≤N≤3 500時(shí)誤差小于4.5%。值得注意的是隨N值的增加，計(jì)算量也會增加，當(dāng)N＞1 800時(shí)計(jì)算效率大大地降低，而且從圖8的曲線走勢可以看出當(dāng)N≥1 800時(shí)誤差呈上升趨勢，因此N的取值不宜太大。

圖6 藥型罩（1／4模型）內(nèi)能時(shí)程曲線Fig.6Internal energy curve of the liner

圖7 藥型罩（1／4模型）動能時(shí)程曲線Fig.7Kinetic energy curve of the liner

圖8 不同N值計(jì)算的射流頭部速度與理論結(jié)果的誤差Fig.8Error between jet tip’s velocity calculated with different Nand theoretical results

3 結(jié) 論

（1）SPH方法與FEM方法相比，避免了網(wǎng)格扭曲、纏繞和物質(zhì)穿透等問題，計(jì)算過程更穩(wěn)定，更適合諸如聚能裝藥等多介質(zhì)、大變形問題。

（2）采用SPH法和FEM方法計(jì)算的射流頭部速度的誤差分別為2.29%、2.53%，射流長度的平均誤差分別為0.143%、4.425%，SPH 方法的計(jì)算結(jié)果更接近理論計(jì)算值。

（3）初始相鄰粒子數(shù)N影響計(jì)算精度和計(jì)算效率：當(dāng)N≤300時(shí)，位于邊界自由面上的粒子出現(xiàn)無物理意義的飛散現(xiàn)象，導(dǎo)致時(shí)間步長逐漸減小，與理論結(jié)果的誤差大于7.5%；當(dāng)400≤N≤3 500時(shí)，與理論結(jié)果的誤差小于4.5%；當(dāng)N≥1 800時(shí)計(jì)算效率大大降低，且誤差呈上升趨勢，因此N的最佳取值范圍為N＝600～1 800，在該范圍內(nèi)誤差小于3%。

參考文獻(xiàn)：

［1］Lucy L B.A numerical approach to the testing of the fission hypothesis［J］.Astronomical Journal，1977，82（12）：1013－1024.

［2］Gingold R A，Monaghan J J.Smoothed particle hydrodynamics－theory and application to non－spherical stars［J］.Monthly Notices of the Royal Astronomical Society，1977（181）：375－389.

［3］Swegle J W，Attaway S W.On the feasibility of using smoothed particle hydrodynamics for underwater explosion calculations［J］.Computational Mechanics，1995，17（3）：151－168.

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

［5］Liu M B，Liu G R，Lam K Y，et al.Meshfree particle simulation of the detonation process for high explosives in shaped charge unlined cavity configurations［J］.Shock Waves，2003，12（6）：509－520.

［6］Liu G R，Liu M B.Smoothed particle hydrodynamics：A meshfree particle method［M］.Singapore：World Scientif－ic，2003.

［7］Libersky L D，Petschek A G，Carney T C，et al.High strain Lagrangian hydrodynamics a three－dimensional SPH code for dynamic material response［J］.Journal of Computational Physics，1993，109（1）：67－75.

［8］Petschek A G，Libersky L D.Cylindrical smoothed particle hydrodynamics［J］.Journal of Computational Physics，1993，109（1）：76－83.

［9］李裕春，時(shí)黨勇，趙遠(yuǎn).ANSYS11.0／LS－DYNA基礎(chǔ)理論與工程實(shí)踐［M］.北京：中國水利水電出版社，2008.

［10］Hallquist J O.LS－DYNA keyword user’s manual［M］.Livermore，CA：Livermore Software Technology Corporation，2007.

［11］隋樹元，王樹山.終點(diǎn)效應(yīng)學(xué)［M］.北京：國防工業(yè)出版社，2000.

［12］Chou P C，Carleone J，F(xiàn)lis W J，et al.Improved formulas for velocity，acceleration，and projection angle of explosively driven liners［J］.Propellants，Explosives，Pyrotechnics，1983，8（6）：175－183.

［13］譚多望，孫承緯.大錐角罩聚能裝藥射流理論計(jì)算方法［J］.高壓物理學(xué)報(bào)，2006，20（3）：270－275.

TAN Duo－wang，SUN Cheng－wei.Analytical model for jet formation in shaped charge with wide cone angle［J］.Chinese Journal of High Pressure Physics，2006，20（3）：270－275.

［14］Swegle J W，Hicks D L，Attaway S W.Smoothed particle hydrodynamics stability analysis［J］.Journal of Computational Physics，1995，116（1）：123－134.