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

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

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

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

1 SPH算法的基本理論

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

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

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

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

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

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

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

2 數(shù)值模擬

2.1 材料模型

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

高速碰撞過程中材料的應(yīng)變率可達105s－1甚至更高，在這樣高應(yīng)變率條件下，材料的動態(tài)力學性能通常會發(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ù)材料動力學實驗得到。

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

式中：pr（υ）和er（υ）分別為參考壓力和參考內(nèi)能，υ為比體積，Γ （υ）為 Grüneisen 系數(shù)。

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

表1 材料的力學性能Table1 Mechanical properties of materials

2.2 有限元模型

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

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

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

3 結(jié)果分析

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

下面主要以初速度為1500m／s的400g長桿彈為例，對整個侵徹過程進行詳細分析。

圖2 長桿彈侵徹3層靶板Fig.2 The three targets impacted by projectile

3.1 彈體對第1層靶板的侵徹

400g長桿彈以1500m／s初速度斜侵徹第1層靶時，鎢合金彈頭在侵徹初期碰撞速度最大，相應(yīng)的碰擊應(yīng)力也最大，遠遠超過了彈體和靶體的動態(tài)強度。鎢合金表現(xiàn)出良好的塑性流動特性，彈體頭部材料主要以反向塑性流動而發(fā)生損耗，不斷向抗力最小的方向飛濺排除，形成動態(tài)穩(wěn)定的侵徹頭部形狀。與此同時，靶板表面的彈坑不斷擴大，在靶內(nèi)建立起相對穩(wěn)定的高壓、高應(yīng)變和高應(yīng)變率狀態(tài)，形成有利于侵徹正常進行的條件，如圖3（a）所示。此時，由于在彈體和靶板內(nèi)各存在沖擊波向遠離接觸面的方向傳播，在彈體和靶板內(nèi)分別形成塑性和彈性等2個響應(yīng)區(qū)。在很短的時間內(nèi)，彈體與靶體接觸部分發(fā)生劇烈的破壞、變形，加之兩者之間強烈的相互摩擦，把彈體的一部分動能迅速轉(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é)束時轉(zhuǎn)而進入相對穩(wěn)定侵徹階段。彈、靶材料繼續(xù)破壞、飛濺。隨著彈體長度不斷縮短，彈坑也不斷擴大，出現(xiàn)新的接觸面，彈體不斷碰撞彈坑的新表面，使彈坑不斷加深。

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

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

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

4種工況下長桿彈貫穿第1層靶板后形成的彈孔尺寸如表2所示。可以看出，計算值與實際尺寸吻合得較好。

3.2 殘余彈體對第2、3層靶的侵徹

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

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

圖6 長桿彈速度與動能的時間歷程曲線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 剩余速度與剩余動能

4種工況的數(shù)值模擬與實驗結(jié)果基本吻合，如表3所示。可以看出，由于長桿彈具有較強的侵徹能力，在貫徹第1層靶板后，形成的破片仍然保持有非常高的動能，足以對后續(xù)靶板進行有效打擊。

3.4 彈靶的破片飛散

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

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

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

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

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

4 結(jié) 論

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

［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ì)點流體動力學（SPH）方法：綜述［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］.國防科技大學學報，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］毛益明，武文遠，陳廣林，等.高速碰撞問題的SPH方法模擬［J］.解放軍理工大學學報：自然科學版，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］.國防科技大學學報，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］隋樹元，王樹山.終點效應(yīng)學［M］.北京：國防工業(yè)出版社，2000.