王金濤，余文力，王 濤，羅永鋒，王少龍

（1.第二炮兵工程大學(xué)，陜西 西安 710025；2.第二炮兵裝備研究院二所，北京 100085）

SPH算法是一種無(wú)網(wǎng)格的自由Lagrange方法，在求解空間導(dǎo)數(shù)時(shí)無(wú)需使用任何網(wǎng)格單元，通過一個(gè)稱為“核函數(shù)”的積分核進(jìn)行“核函數(shù)估值”近似將流體動(dòng)力學(xué)方程轉(zhuǎn)換成數(shù)值計(jì)算的SPH方程。在模擬物體大變形時(shí)，可以避免傳統(tǒng)Lagrange方法產(chǎn)生的網(wǎng)格扭曲和無(wú)物理意義的侵蝕算法對(duì)網(wǎng)格的刪除，同時(shí)也克服了Euler方法難以清晰跟蹤物質(zhì)變形和不能識(shí)別材料界面位形的缺點(diǎn)。SPH算法在結(jié)構(gòu)動(dòng)力學(xué)應(yīng)用方面得到了快速的發(fā)展，L.D.Libersky等［1］將材料強(qiáng)度效應(yīng)引入SPH算法，并成功應(yīng)用于高速碰撞數(shù)值模擬領(lǐng)域；G.R.Johnson等［2］和J.Campbell等［3］分別將SPH算法應(yīng)用于侵徹貫穿問題的研究。A.N.Parshikov［4］基于黎曼解的思想來(lái)描述粒子之間的相互作用，改進(jìn)了界面的計(jì)算精度。

本文中基于Autodyn顯式有限元分析軟件，應(yīng)用SPH算法，對(duì)長(zhǎng)桿彈斜侵徹多層金屬間隔靶板進(jìn)行了數(shù)值模擬，得到了侵徹過程的物理圖像，數(shù)值模擬與實(shí)驗(yàn)現(xiàn)象吻合較好。

1 SPH算法的基本理論

SPH算法的核心是一種插值算法。設(shè)任意宏觀變量（如密度、壓力、溫度等）在空間中某一點(diǎn)r處的場(chǎng)f（r），通常稱〈f（r）〉為f（r）的一個(gè)核估計(jì)，〈f（r）〉可以通過函數(shù)f（r）在一組無(wú)序點(diǎn)上的值表示成積分插值并計(jì)算得到

核函數(shù)W 的意義在于圍繞場(chǎng)f（r）產(chǎn)生一個(gè)光滑器或過濾器的作用，將其中的局部統(tǒng)計(jì)漲落都過濾掉，從而產(chǎn)生一個(gè)場(chǎng)f（r）的估計(jì)。

設(shè)函數(shù)f（r）在粒子j上的值為fj，則式（1）可以改寫成粒子求和的形式

式中：N為求解區(qū)域內(nèi)的粒子總數(shù)，mj、fj、ρj分別為位于空間點(diǎn)rj處粒子j的質(zhì)量、待求參數(shù)和密度。同樣，函數(shù)f（r）在粒子i上一階導(dǎo)數(shù)的核估計(jì)為

核函數(shù)的形式有很多種，例如，采用B樣條函數(shù)作為核函數(shù)

在全應(yīng)力張量空間中，連續(xù)介質(zhì)力學(xué)中質(zhì)點(diǎn)的運(yùn)動(dòng)方程以及粒子的質(zhì)量、動(dòng)量、能量守恒方程用SPH插值公式可表示為

式（5～8）中：vα、vβ為粒子速度矢量分量，α、β表示空間坐標(biāo)軸方向。

2 數(shù)值模擬

2.1 材料模型

高速碰撞過程是一個(gè)極為復(fù)雜的瞬態(tài)物理過程。對(duì)高速碰撞現(xiàn)象的描述涉及連續(xù)介質(zhì)力學(xué)的多個(gè)方面，需要綜合考慮如濺射等多種復(fù)雜的物理現(xiàn)象，對(duì)于數(shù)值模擬來(lái)說(shuō)，材料模型的選取和材料參數(shù)的設(shè)定是有效完成數(shù)值模擬的重要環(huán)節(jié)。

高速碰撞過程中材料的應(yīng)變率可達(dá)105s－1甚至更高，在這樣高應(yīng)變率條件下，材料的動(dòng)態(tài)力學(xué)性能通常會(huì)發(fā)生巨大的變化。Johnson－Cook模型能夠在較大的應(yīng)變率范圍內(nèi)較好地描述材料的本構(gòu)關(guān)系，它定義屈服應(yīng)力為

式中：εp為等效塑性應(yīng)變?yōu)橐?guī)范等效塑性應(yīng)變率，T＊＝（T－Tr）／（Tm－Tr），Tr為室溫，Tm為熔化溫度。參數(shù)A、B、C、n、m根據(jù)材料動(dòng)力學(xué)實(shí)驗(yàn)得到。

為了描述材料在高壓下的流動(dòng)壓力與密度、內(nèi)能等的關(guān)系，選用Mie－Grüneisen狀態(tài)方程

式中：pr（υ）和er（υ）分別為參考?jí)毫蛥⒖純?nèi)能，υ為比體積，Γ （υ）為 Grüneisen 系數(shù)。

長(zhǎng)桿彈與3層靶板的材料分別采用鎢合金、603裝甲鋼、45鋼及LY12鋁，相關(guān)力學(xué)性能參數(shù)如表1所示，其中ρ為材料密度，E為材料的楊氏模量，γ為材料的剪切模量，σy為材料屈服極限，ν為泊松比。

表1 材料的力學(xué)性能Table1 Mechanical properties of materials

2.2 有限元模型

設(shè)計(jì)2種長(zhǎng)徑比均為10的長(zhǎng)桿彈。彈種1，質(zhì)量為200g，直徑為12mm，長(zhǎng)為120mm；彈種2，質(zhì)量為400g，彈直徑為16mm，長(zhǎng)為160mm。第1層靶板厚45mm，第2層靶板厚14mm，第3層靶板厚2mm，第1、2層靶間隔2.7m，第2、3層靶間隔1.5m。2種長(zhǎng)桿彈的初始條件相同，彈著角均為60°，初始速度均分別為1000、1500m／s。對(duì)應(yīng)4種工況，進(jìn)行實(shí)彈打靶實(shí)驗(yàn)。

對(duì)長(zhǎng)桿彈和第1層靶板采用SPH算法建模，第2、3層靶板使用Lagrange有限元網(wǎng)格，以降低計(jì)算時(shí)間。長(zhǎng)桿彈與第1層靶板的SPH模型如圖1（a）所示。第2、3層靶板如圖1（b）所示。

圖1 有限元模型Fig.1 Finite elements model

3 結(jié)果分析

利用Autodyn軟件分別對(duì)4種工況進(jìn)行仿真，以400g長(zhǎng)桿彈、1500m／s初速的工況為例，數(shù)值模擬結(jié)果如圖2所示。長(zhǎng)桿彈貫穿第1層靶板后，殘余彈體與大量高能破片對(duì)第2層靶板造成貫穿打擊，形成5個(gè)貫穿孔洞及大量小彈坑。貫穿第2層靶板后，彈體已完全破碎，但仍然有相當(dāng)數(shù)量的高動(dòng)能破片，且破片散布范圍更大，在第3層靶板上形成大量大小不等的孔洞。

下面主要以初速度為1500m／s的400g長(zhǎng)桿彈為例，對(duì)整個(gè)侵徹過程進(jìn)行詳細(xì)分析。

圖2 長(zhǎng)桿彈侵徹3層靶板Fig.2 The three targets impacted by projectile

3.1 彈體對(duì)第1層靶板的侵徹

400g長(zhǎng)桿彈以1500m／s初速度斜侵徹第1層靶時(shí)，鎢合金彈頭在侵徹初期碰撞速度最大，相應(yīng)的碰擊應(yīng)力也最大，遠(yuǎn)遠(yuǎn)超過了彈體和靶體的動(dòng)態(tài)強(qiáng)度。鎢合金表現(xiàn)出良好的塑性流動(dòng)特性，彈體頭部材料主要以反向塑性流動(dòng)而發(fā)生損耗，不斷向抗力最小的方向飛濺排除，形成動(dòng)態(tài)穩(wěn)定的侵徹頭部形狀。與此同時(shí)，靶板表面的彈坑不斷擴(kuò)大，在靶內(nèi)建立起相對(duì)穩(wěn)定的高壓、高應(yīng)變和高應(yīng)變率狀態(tài)，形成有利于侵徹正常進(jìn)行的條件，如圖3（a）所示。此時(shí)，由于在彈體和靶板內(nèi)各存在沖擊波向遠(yuǎn)離接觸面的方向傳播，在彈體和靶板內(nèi)分別形成塑性和彈性等2個(gè)響應(yīng)區(qū)。在很短的時(shí)間內(nèi)，彈體與靶體接觸部分發(fā)生劇烈的破壞、變形，加之兩者之間強(qiáng)烈的相互摩擦，把彈體的一部分動(dòng)能迅速轉(zhuǎn)換為熱能，碎片被加熱到熾熱的程度飛濺出去，從而在彈、靶接觸初期看到一片“火光”，如圖3（b）所示。

圖3 彈體侵徹第1層靶板初期Fig.3 Beginning period of the first target impacted by projectile

圖4 彈體侵徹第1層靶板中后期Fig.4 Medium and late period of the first target impacted by projectile

圖5 靶板正面與背面的彈孔Fig.5 The front and back holes

在開坑階段結(jié)束時(shí)轉(zhuǎn)而進(jìn)入相對(duì)穩(wěn)定侵徹階段。彈、靶材料繼續(xù)破壞、飛濺。隨著彈體長(zhǎng)度不斷縮短，彈坑也不斷擴(kuò)大，出現(xiàn)新的接觸面，彈體不斷碰撞彈坑的新表面，使彈坑不斷加深。

在斜侵徹過程中，由于彈體頭部受到的抗力方向并不均勻，大部分破碎的彈體和靶板碎片向著抗力最小的方向飛濺出彈坑。在高速碰撞中，彈體的跳飛趨勢(shì)并不明顯，彈坑沒有呈現(xiàn)明顯的下凹形狀，數(shù)值模擬和實(shí)驗(yàn)觀察都體現(xiàn)了這一點(diǎn)，如圖4（a）所示。隨著侵徹的不斷深入，當(dāng)侵徹深度超過靶厚的1／2后，彈體和靶板的碎片都不容易順利排出，碎片擠在彈體頭部周圍使得侵徹孔洞不斷擴(kuò)大，同時(shí)在靶板背表面形成鼓包，并出現(xiàn)裂紋。

在侵徹的最后階段，由于靶板背表面鼓包的增大和裂紋的持續(xù)發(fā)展，彈體受到的抗力迅速下降，鼓包的薄弱處在剪切力作用下發(fā)生斷裂，出現(xiàn)沖塞塊和大量破片，伴隨著殘余

彈體穿孔噴出，侵徹過程結(jié)束，如圖4（b）所示。所形成的實(shí)際彈孔形狀如圖5所示。

4種工況下長(zhǎng)桿彈貫穿第1層靶板后形成的彈孔尺寸如表2所示?？梢钥闯?，計(jì)算值與實(shí)際尺寸吻合得較好。

3.2 殘余彈體對(duì)第2、3層靶的侵徹

殘余彈體對(duì)第2層靶板的侵徹作用是考核長(zhǎng)桿彈殺傷效能的重要指標(biāo)，彈體貫穿第1層靶板后，殘余彈體和破片對(duì)第2、3層靶板仍具有較強(qiáng)的侵徹能力。以殘余彈體為例，即使對(duì)于1g的破片，穿過第1層靶板后，破片的剩余速度為1 334.6m／s，剩余動(dòng)能為890J，如圖6所示，遠(yuǎn)遠(yuǎn)大于對(duì)人員的殺傷標(biāo)準(zhǔn)78J［9］，整個(gè)彈體的總動(dòng)能在貫穿第1層靶板后仍保持在200kJ的數(shù)量級(jí)（計(jì)算時(shí)采用1／2模型）。第3層靶板上的彈孔數(shù)量更多，如圖7所示，說(shuō)明了長(zhǎng)桿彈的破碎更為嚴(yán)重，因此數(shù)值仿真在模擬材料破碎方面與實(shí)驗(yàn)結(jié)果還有較大出入。

表2 第1層靶板彈洞尺寸Table2 Hole sizes in the first target

圖6 長(zhǎng)桿彈速度與動(dòng)能的時(shí)間歷程曲線Fig.6 Velocity and kinetic curves of a 400g projectile with the initial velocity of 1500m／s after its penetration into the first target

圖7 殘余彈體侵徹第2、3層靶板Fig.7 The second and third targets impacted by the residual projectile

3.3 剩余速度與剩余動(dòng)能

4種工況的數(shù)值模擬與實(shí)驗(yàn)結(jié)果基本吻合，如表3所示?？梢钥闯觯捎陂L(zhǎng)桿彈具有較強(qiáng)的侵徹能力，在貫徹第1層靶板后，形成的破片仍然保持有非常高的動(dòng)能，足以對(duì)后續(xù)靶板進(jìn)行有效打擊。

3.4 彈靶的破片飛散

根據(jù)數(shù)值模擬結(jié)果，彈體貫穿第1層靶板后形成大量破片（彈體和第1層靶板形成的破片均包含在內(nèi)），統(tǒng)計(jì)數(shù)據(jù)如表4所示。可以看出，在初速度相同的情況下，400g長(zhǎng)桿彈所形成的質(zhì)量大于1g的破片和高動(dòng)能破片（Ek＞100J）的數(shù)量均多于200g長(zhǎng)桿彈；而對(duì)于同種長(zhǎng)桿彈在初速度不同的情況下，以1500m／s的初速度貫穿第1層靶板后，雖然形成的質(zhì)量大于1g的破片數(shù)量少于初速度為1000m／s時(shí)的情形，但是高動(dòng)能破片數(shù)量明顯多于低速情況。這也說(shuō)明了鎢合金在高速侵徹時(shí)具有良好的易碎性，容易形成大量破片，從而擴(kuò)大殺傷范圍。

表3 貫穿第1層靶板后的剩余速度與剩余動(dòng)能Table3 The residual velocity and kinetic energy of projectiles after their penetration into the first target

表4 破片數(shù)量統(tǒng)計(jì)Table4 The number of fragments

圖8為400g長(zhǎng)桿彈以1500m／s的初速度侵徹第1層靶板后破片的飛散情況，除了幾塊質(zhì)量和動(dòng)能都比較大的破片外，絕大多數(shù)都是質(zhì)量非常小的金屬微粒，這些微粒具有非常高的溫度，看上去如同大片“火光”，“火光”的范圍與數(shù)值模擬結(jié)果所展示的微粒飛散范圍大致相同。對(duì)于質(zhì)量和動(dòng)能都較大的破片，能夠?qū)⒌?層靶板貫穿，如圖7所示，仿真結(jié)果顯示在第2層靶板形成了5個(gè)彈洞，彈洞的數(shù)量和分布與實(shí)驗(yàn)現(xiàn)象基本吻合。

圖8 400g長(zhǎng)桿彈貫穿第1層靶板后的破片飛散情況Fig.8 Fragments after a 400g projectile penetrating into the first target

4 結(jié) 論

本文中應(yīng)用SPH算法模擬了長(zhǎng)桿彈對(duì)多層靶板的侵徹過程，得到了2種長(zhǎng)桿彈高速侵徹的數(shù)值模擬結(jié)果，詳細(xì)分析了長(zhǎng)桿彈貫穿第1層中厚靶的3個(gè)階段及貫穿后破片的形成和飛散過程，并統(tǒng)計(jì)了破片的質(zhì)量、動(dòng)能及數(shù)量。動(dòng)態(tài)實(shí)驗(yàn)結(jié)果表明：數(shù)值計(jì)算結(jié)果與實(shí)驗(yàn)現(xiàn)象吻合很好；采用SPH算法對(duì)高速侵徹過程進(jìn)行數(shù)值模擬能夠較好地描述長(zhǎng)桿彈對(duì)多層間隔靶板的高速侵徹過程。

［1］Libersky L D，Petschek A G.Smooth particle hydrodynamics with strength of materials［C］∥Advances in the Free Lagrange Method，Lecture Notes in Physics.Oxford：AWE Hunting Press，1991：248－257.

［2］Johnson G R，Beissel S R.Normalized smoothing functions for SPH impact computations［J］.International Journal for Numerical Methods in Engineering，1996，39（16）：2725－2741.

［3］Campbell J，Vignjevic R，Libersky L D.A contact algorithm for smoothed particle hydrodynamics［J］.Computer Methods in Applied Mechanics and Engineering，2000，184（1）：49－65.

［4］Parshikov A N，Medin S A，Loukashenko I I，et al.Improvements in SPH method by means of interparticle contact algorithm and analysis of perforation tests at moderate projectile velocities［J］.International Journal of Impact Engineering，2000，24（8）：779－796.

［5］張鎖春.光滑質(zhì)點(diǎn)流體動(dòng)力學(xué)（SPH）方法：綜述［J］.計(jì)算物理，1996，13（4）：385－397.ZHANG Suo－chun.Smoothed particle hydrodynamics（SPH）method：A Review［J］.Chinese Journal of Computational，1996，13（4）：385－397.

［6］徐志宏，湯文輝，羅永.SPH 算法在高速侵徹問題中的應(yīng)用［J］.國(guó)防科技大學(xué)學(xué)報(bào)，2005，27（4）：45－48.XU Zhi－h(huán)ong，TANG Wen－h(huán)ui，LUO Yong.Smoothed particle hydrodynamics algorithm applied in penetration problem［J］.Journal of National University of Defense Technology，2005，27（4）：45－48.

［7］毛益明，武文遠(yuǎn)，陳廣林，等.高速碰撞問題的SPH方法模擬［J］.解放軍理工大學(xué)學(xué)報(bào)：自然科學(xué)版，2003，4（5）：86－89.MAO Yi－ming，WU Wen－yuan，CHEN Guang－lin，et al.Numerical simulation of high velocity impact problems with SPH method［J］.Journal of PLA University of Science and Technology：Natural Science Edition，2003，4（5）：86－89.

［8］崔偉峰，曾新吾.SPH 算法在超高速碰撞數(shù)值模擬中的應(yīng)用［J］.國(guó)防科技大學(xué)學(xué)報(bào)，2007，29（2）：47－50.CUI Wei－feng，ZENG Xin－wu.Smoothed particle hydrodynamics algorithm applied in numerical simulation of hypervelocity impact［J］.Journal of National University of Defense Technology，2007，29（2）：47－50.

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