徐思博， 孟子飛， 劉文韜， 曹雪雁

(哈爾濱工程大學(xué) 船舶工程學(xué)院，哈爾濱 150001)

高速破片作為武器爆炸攻擊艦船的重要利器之一，其產(chǎn)生的高壓環(huán)境難以測量，初始速度通常在1 000 m/s左右，其形狀的隨機(jī)性以及強(qiáng)穿透性，使得如何防護(hù)破片毀傷成為研究的難點(diǎn)。如今國際公認(rèn)的防護(hù)手段是在艙壁外加設(shè)防護(hù)液艙用來衰減破片的沖擊，研究水下近場爆炸高速破片在液艙中的衰減特性有極其重要的軍事戰(zhàn)略意義[1]。針對高速破片衰減特性，目前可行的方法主要有基于理想模型的數(shù)值驗(yàn)證和一些實(shí)驗(yàn)驗(yàn)證。高速破片擊穿船舶液艙的整個(gè)過程大致分為沖擊、拉伸、空穴形成和潰滅以及穿透內(nèi)板四個(gè)階段[2]。沖擊階段，破片接觸液艙前壁并開始穿透，在艙壁后的液體內(nèi)產(chǎn)生擾動(dòng)；拉伸階段，高速破片撞擊液體并產(chǎn)生沖擊波，之后在液體中前進(jìn)；空穴的形成和潰滅階段，高速破片在液體中的速度被迅速衰減，形成空穴并演化直至潰滅；穿透內(nèi)板階段，高速破片抵達(dá)液艙內(nèi)壁并穿出。

穿透過程中有幾個(gè)關(guān)鍵性的難點(diǎn)，包括：材料的動(dòng)態(tài)力學(xué)性能和力學(xué)響應(yīng)，高速破片水中速度衰減問題，水中沖擊波形成的機(jī)理，空穴的形成和沖擊波的疊加等。20世紀(jì)初期， Worthington等[3]用高速攝像機(jī)記錄下圓球落水的動(dòng)態(tài)過程中的飛濺和空泡現(xiàn)象，并以此為基礎(chǔ)開啟水彈道學(xué)的發(fā)展。我國在這項(xiàng)領(lǐng)域的研究開始于20世紀(jì)末期，僅做過少量的入水沖擊、入水空泡實(shí)驗(yàn)[4-6]。陳先富[7]進(jìn)行了不同形狀彈丸入水的試驗(yàn)研究，記錄了空穴的產(chǎn)生潰滅以及脈動(dòng)，但由于試驗(yàn)工況較少，僅得出初始速度越高空穴增大越快這一結(jié)論。高速侵徹實(shí)驗(yàn)危險(xiǎn)性大，實(shí)驗(yàn)條件要求高，同時(shí)隨著計(jì)算機(jī)技術(shù)的飛速發(fā)展，對侵徹問題應(yīng)用數(shù)值仿真方式進(jìn)行模擬研究成為了解決該問題的有效且經(jīng)濟(jì)的途徑之一。在21世紀(jì)，徐雙喜等[8]依照破片穿甲運(yùn)動(dòng)方程以及德·瑪爾模型，分析得出圓柱形破片穿透背水板的剩余速度公式，并采用Autodyn中的數(shù)值耦合算法模擬該過程，孔祥韶等[9]同時(shí)對雙發(fā)破片的相互耦合作用模式進(jìn)行分析，其主要研究為穿透后的剩余速度，并未涉及空穴演化等現(xiàn)象。Varas等[10]開展了低速以及高速下的球型破片擊穿液艙的實(shí)驗(yàn)，并用高速成像儀記錄下了高速破片形成的沖擊波和空穴，并運(yùn)用ALE算法進(jìn)行了數(shù)值仿真，但彈片衰減速度、水中壓力變化等仍有較大偏差。張婧等[11-13]則關(guān)注到多層防護(hù)結(jié)構(gòu)及復(fù)合材料裝甲加持下破片對結(jié)構(gòu)的毀傷問題，與此同時(shí)沈曉樂等[14]在試驗(yàn)中對方形破片阻力系數(shù)進(jìn)行修正，并通過試驗(yàn)研究破片的侵徹能力以及剩余特性。近幾年，國內(nèi)針對破片形狀對侵徹的影響做了很多工作[15-17]，但并未引進(jìn)新的仿真方法。楊文山等[18-19]則應(yīng)用SPH方法對侵徹過程進(jìn)行仿真，這種無網(wǎng)格光滑粒子在處理破口界面處有極好的模擬效果。

本文根據(jù)Varas等開展的球型彈片擊穿液艙的實(shí)驗(yàn)數(shù)據(jù)，應(yīng)用有限元軟件ABAQUS中的CEL算法對該過程進(jìn)行數(shù)值仿真。試驗(yàn)針對不同初始速度和裝水量下的壓力測點(diǎn)、空穴徑寬、彈片速度衰減與實(shí)驗(yàn)數(shù)據(jù)進(jìn)行對比，驗(yàn)證了CEL在模擬高速侵徹過程的適用性問題。

1 CEL算法的基本方程

CEL方法同時(shí)具有拉格朗日方法和歐拉方法的優(yōu)勢，較好的處理大變形、碰撞及一些流固耦合問題[20]。CEL在處理耦合問題中，應(yīng)用拉格朗日算法計(jì)算結(jié)構(gòu)以保持其在響應(yīng)上的精度，用歐拉算法計(jì)算流體介質(zhì)以解決其流動(dòng)過程中的大變形問題。拉格朗日算法是一種離散化數(shù)值求解方法，其求解過程中，材料會(huì)依附在單元網(wǎng)格上，在網(wǎng)格的運(yùn)動(dòng)和變形的過程中伴隨著材料的流動(dòng)，計(jì)算過程中，單元會(huì)發(fā)生變形但單元的初始質(zhì)量不變。具體控制方程為

(1)

(2)

(3)

式中:ρ,u,e,σ分別為拉格朗日流體介質(zhì)的密度、速度、比內(nèi)能和應(yīng)力。

歐拉算法是一種固定網(wǎng)格而使材料在網(wǎng)格間運(yùn)輸?shù)那蠼夥椒?，?jì)算過程中一個(gè)歐拉單元按用戶定義的材料分配，有著不同的物質(zhì)，整個(gè)計(jì)算過程中空間網(wǎng)格將會(huì)固定，而材料網(wǎng)格將以拉格朗日單元的狀態(tài)變量實(shí)時(shí)傳輸?shù)娇臻g網(wǎng)格中，具體控制方程為

(4)

(5)

(6)

式中:ρ,u,e,σ,f分別為歐拉流體介質(zhì)的密度、速度、比內(nèi)能、應(yīng)力和質(zhì)量力。

在耦合算法之前，這兩種算法是相互獨(dú)立的，同時(shí)在耦合計(jì)算時(shí)通過耦合面的控制，歐拉的材料流動(dòng)也不會(huì)穿透拉格朗日單元。耦合作用則通過在界面交換壓力載荷和幾何約束來實(shí)現(xiàn)，如圖1所示。

圖1 耦合界面處理示意圖Fig.1 Coupling interface processing

2 高速破片穿透液艙的數(shù)值分析

高速破片穿透液艙的數(shù)值模型使用有限元仿真進(jìn)行建立運(yùn)算分析，符合工程應(yīng)用背景以及科學(xué)研究需求。

2.1 液艙結(jié)構(gòu)

鑒于有限元分析精準(zhǔn)度與網(wǎng)格尺寸密切相關(guān)，為提高運(yùn)算精度以及計(jì)算效率，仿真分析將采用實(shí)驗(yàn)整體模型的一半。根據(jù)Varas等實(shí)驗(yàn)中液艙的尺寸結(jié)構(gòu)，液艙用6063-T5鋁制成，長750 mm，寬150 mm，壁厚2.5 mm，液艙兩側(cè)用厚30 mm的PMMA有機(jī)玻璃封閉，在液艙內(nèi)部設(shè)有兩個(gè)壓力測點(diǎn)PTn與PTf，液艙外部前后壁板設(shè)有應(yīng)力及撓度測點(diǎn)G1與G2，測點(diǎn)位置如圖4所示。破片為球形，直徑12.5 mm，質(zhì)量8 g，實(shí)驗(yàn)中選取600 m/s和900 m/s了兩種初始速度。本次仿真針對該模型在壁厚方向(即彈片前進(jìn)方向)設(shè)置3層網(wǎng)格單元，接近彈片侵徹位置附近皆為正六面體網(wǎng)格，網(wǎng)格尺寸為0.83 mm，最終液艙模型網(wǎng)格單元數(shù)為49 540個(gè)，如圖2和圖3所示。

仿真試驗(yàn)采用J-C模型描述液艙的本構(gòu)關(guān)系。J-C模型的關(guān)系式為

(7)

(8)

式中：θ為當(dāng)前溫度；θmelt為熔化溫度；θtransition為轉(zhuǎn)化溫度。

在J-C失效模型中，斷裂應(yīng)變定義式為

(9)

表1 結(jié)構(gòu)材料表

圖2 液艙和破片數(shù)值模型Fig.2 Numerical modelling of fuel tanks and fragment

圖3 破片入射面網(wǎng)格圖Fig.3 Detail of the entry wall mesh

圖4 壓力測點(diǎn)位置示意圖Fig.4 Sketch of pressure transducers

2.2 歐拉域

鑒于高速?zèng)_擊問題帶來的大變形以及流體的不規(guī)則流動(dòng)，仿真試驗(yàn)將用歐拉網(wǎng)格建立水和空氣的模型。在彈片穿透過程中，結(jié)構(gòu)的大變形同時(shí)會(huì)伴隨著流體的運(yùn)動(dòng)，因此歐拉域的設(shè)置應(yīng)大于結(jié)構(gòu)的區(qū)域以保證流體在流動(dòng)過程中不會(huì)消失，保證能量守恒，所以歐拉域設(shè)置為190 mm寬。水和空氣的材料屬性通過體積分?jǐn)?shù)方式離散進(jìn)歐拉域內(nèi)，最終歐拉網(wǎng)格單元數(shù)為4 061 300個(gè)，如圖5所示。

圖5 歐拉網(wǎng)格展示圖Fig.5 Mesh of the fluids

水的狀態(tài)方程采用Mie-Grüneisen狀態(tài)方程，定義式為

(10)

式中：ρ0為流體密度；S，c0，Γ0為流體常數(shù)；η=1-ρ0/ρ?？諝獠捎美硐霘怏w模型。水和空氣的具體材料數(shù)據(jù)見表2。

表2 水和空氣的材料參數(shù)表

3 數(shù)值計(jì)算結(jié)果分析

近場水下爆炸產(chǎn)生的破片速度通常在1 000 m/s左右，因此在本節(jié)中，將給出裝水量75%和60%，彈片初始速度從600～1 200 m/s速度梯度下的工況數(shù)值仿真結(jié)果并與Varas等的實(shí)驗(yàn)數(shù)據(jù)進(jìn)行對比，并分析裝水量75%時(shí)不同入射速度下前后壁板應(yīng)力及撓度對比。

3.1 破片高速侵徹過程

由仿真的結(jié)果可知，破片侵徹過程分為四個(gè)階段，即沖擊、拉伸、空穴形成和潰滅以及穿透內(nèi)板，每個(gè)階段在數(shù)值仿真中得到了充足的展現(xiàn)。在彈片穿透的瞬間，結(jié)構(gòu)的應(yīng)力峰值達(dá)到257 MPa，并沿沖擊點(diǎn)向四周擴(kuò)散圓形波，如圖6(a)所示。在沖擊階段，當(dāng)彈片穿透壁板之后，彈片的動(dòng)能迅速在流體內(nèi)形成高壓的環(huán)境，并向四周擴(kuò)散沖擊波，如圖6(b)所示。在拉伸階段和空穴形成潰滅的階段，彈片穿透水箱在流體內(nèi)部留下一片空穴區(qū)，流體在慣性的作用下不斷擴(kuò)散使得空穴不斷變化，而空穴潰滅所需時(shí)間遠(yuǎn)遠(yuǎn)大于高速彈片穿透水箱所需時(shí)間，如圖6(c)所示。在穿透階段，由于流體內(nèi)沖擊波預(yù)先到達(dá)壁板，彈片在到達(dá)之前就在壁板處產(chǎn)生預(yù)應(yīng)力，圖6(d)給出穿透階段液艙的應(yīng)力云圖。

圖6 初始速度900 m/s裝水量75%的高速穿透過程Fig.6 HRAM phases in a tube impacted at 900 m/s and filled 75%

3.2 不同入射速度和裝載量對破片速度衰減的影響

圖7給出四組工況下速度衰減曲線與理論分析解和Varas等實(shí)驗(yàn)值得對比情況?？梢钥闯鰯?shù)值分析方法CEL的結(jié)果與理論分析結(jié)果非常接近，只是CEL算法在彈片穿透液艙時(shí)會(huì)有較大的減速。由圖7看出，CEL算法在彈片速度衰減趨勢上與理論分析和實(shí)驗(yàn)有較好的擬合；由于彈片在液艙內(nèi)滯留時(shí)間極短，液體阻滯作用有限，同時(shí)液面距彈道位置較遠(yuǎn)，液面效應(yīng)較小，液艙內(nèi)裝水量對速度衰減趨勢并沒有影響；彈片在穿透入射面時(shí)速度大致衰減15%左右。

圖7 速度衰減時(shí)歷曲線對比(實(shí)驗(yàn)值與理論值源自Varas等的研究)Fig.7 Comparison of velocity decay vs. time (data of experiment and analysis come from Varas et al)

3.3 不同液體裝載量對彈道空穴演化的影響

圖8和圖9給出由Varas等的實(shí)驗(yàn)數(shù)據(jù)中特定時(shí)刻空穴大小與CEL同一時(shí)刻空穴大小對比?？梢钥闯鯟EL仿真結(jié)果與實(shí)驗(yàn)數(shù)據(jù)空穴大小、水面隆起的波高等十分接近。同時(shí)也發(fā)現(xiàn)，入射面與射出面的箱體呈現(xiàn)向外隆起的變形，由于穿透瞬間與彈片接觸位置應(yīng)力高達(dá)257 MPa，高于壁板的屈服應(yīng)力使其瞬間形成穿透圓孔，而當(dāng)彈片穿透之后，在內(nèi)部流體的作用下，壁板所受應(yīng)力維持在100～200 MPa，使液艙壁板進(jìn)入塑性變形階段發(fā)生屈曲卻不至于形成斷裂破損，而此時(shí)液體量的多少則是影響壁板屈曲變形的主要因素，如圖8和圖9所示，圖中框線為液面初始位置。

液體裝載量的變化決定了彈道距離自由液面的相對位置，并影響及導(dǎo)致彈道出現(xiàn)明顯的非對稱空穴現(xiàn)象，液體裝載量從75%減小到60%，即彈道與自由液面的距離從6倍彈片半徑減小到2.4倍彈片半徑，非對稱差異可由20%增加到48%，這種非對稱差異導(dǎo)致空穴內(nèi)壓不平衡從而對彈道穩(wěn)定性的影響很大。

圖8 初始速度600 m/s裝水量75%下空穴演化0.26 ms,0.46 ms,0.60 ms時(shí)刻實(shí)驗(yàn)數(shù)據(jù)與CEL算法數(shù)據(jù)的對比Fig.8 Cavity evolution at 0.26 ms, 0.46 ms, 0.60 ms obtained from experiments and CEL simulation for a tube impacted at 600 m/s, filled 75%

圖9 初始速度600 m/s裝水量60%下空穴演化0.13 ms,0.24 ms,0.38 ms時(shí)刻實(shí)驗(yàn)數(shù)據(jù)與CEL算法數(shù)據(jù)的對比Fig.9 Cavity evolution at 0.13 ms, 0.24 ms, 0.38 ms obtained from experiments and CEL simulation for a tube impacted at 600 m/s, filled 60%

3.4 不同入射速度對液艙隔板壓力和變形的影響

根據(jù)實(shí)驗(yàn)中兩個(gè)壓力測點(diǎn)PTn，PTf與數(shù)值仿真結(jié)果進(jìn)行對比，對比壓力時(shí)歷曲線如圖10和圖11所示。由于實(shí)驗(yàn)與數(shù)值仿真中彈片速度存在偏差，因此反映到壓力時(shí)歷曲線中對應(yīng)出現(xiàn)峰值的時(shí)刻不盡相同，同時(shí)測點(diǎn)壓力的峰值有所偏差，尤其是PTf，但兩次壓力峰值的脈寬十分接近，結(jié)果基本吻合；在Ptn測點(diǎn)，幾次仿真結(jié)果均出現(xiàn)沖擊波的二次加載，其峰值大致為第一次波峰的15%，兩次波峰間脈寬小于0.05 ms，根據(jù)Ptn測點(diǎn)空間位置，推測二次加載為液艙底面反射稀疏波疊加形成(若是由射出面反射形成則脈寬應(yīng)在0.12 ms左右，假定沖擊波為球面波)。

圖10 裝水量75%初始速度900 m/s的實(shí)驗(yàn)值與數(shù)值的壓力時(shí)歷曲線對比Fig.10 Experimental and numerical pressure time historyPtn and Ptf in a tube 75% filled impacted at 900 m/s

圖11 裝水量75%初始速度600 m/s的實(shí)驗(yàn)值與數(shù)值的壓力時(shí)歷曲線對比Fig.11 Experimental and numerical pressure time historyPtn and Ptf in a tube 75% filled impacted at 600 m/s

結(jié)構(gòu)在達(dá)到塑性變形之前，分別經(jīng)歷線彈性階段以及非線性彈性階段，前置壁板在不同速度下應(yīng)力響應(yīng)如圖12所示，不同入射速度下結(jié)構(gòu)的應(yīng)力響應(yīng)僅因穿透時(shí)間不同而有一定的時(shí)域上的滯后性，而結(jié)構(gòu)的應(yīng)力響應(yīng)只與材料本構(gòu)方程有關(guān)，前置壁板在破片穿透過程中達(dá)到140 MPa左右應(yīng)力，并在反射稀疏波抵達(dá)前置壁板前有一定的應(yīng)力卸載，最終在稀疏波與水波的共同作用下達(dá)到動(dòng)態(tài)屈服強(qiáng)度時(shí)應(yīng)力達(dá)到最大值，此現(xiàn)象在后置壁板的應(yīng)力響應(yīng)中有更為明顯的表現(xiàn)，如圖13所示。圖14給出在前后壁板測點(diǎn)位置最終撓度，隨著初始速度的增加，前后隔板的變形撓度呈現(xiàn)出不同特征，即前置隔板的變形近乎呈線性關(guān)系，但幅度平緩，而后置隔板的變形呈非線性增長趨勢。

圖12 裝水量75%不同速度下G1應(yīng)力曲線Fig.12 Pressure time history G1 for a tube impacted at different velocity, filled 75%

圖13 裝水量75%不同速度下G2應(yīng)力曲線Fig.13 Pressure time history G2 for a tube impacted at different velocity, filled 75%

圖14 裝水量75%壁板測點(diǎn)最終撓度Fig.14 Deflection in the end on the entry and exit wall of a tube impacted at different velocity, filled 75%

4 結(jié) 論

本文基于CEL算法研究了破片的高速侵徹問題，數(shù)值結(jié)果呈現(xiàn)出明顯的四個(gè)階段，即沖擊、拉伸、空穴形成和潰滅以及穿透內(nèi)板。文中首先以Varas等開展的高速彈片穿透鋁制液艙實(shí)驗(yàn)為基礎(chǔ)模型，驗(yàn)證了數(shù)值模擬結(jié)果的可靠性，進(jìn)而模擬了不同破片速度、不同液艙裝載量時(shí)高速侵徹過程的空穴演化特征，得出了以下結(jié)論：

(1) 在當(dāng)前液艙參數(shù)下，彈片在穿透前置隔板時(shí)速度大致衰減15%左右；當(dāng)彈道遠(yuǎn)離自由液面時(shí)，液艙裝載量對破片速度衰減規(guī)律影響較小。

(2) 近自由液面入射彈道會(huì)出現(xiàn)明顯的非對稱空穴演化特征，隨著與自由液面距離變小，非對稱性差異可由20%增加到48%，這對彈道的穩(wěn)定性影響很大。

(3) 隨著破片初始速度增加，前后隔板的變形撓度呈現(xiàn)出不同特征，在當(dāng)前速度范圍內(nèi)，前置隔板的變形隨著破片入射速度增加而略有增加，但幅度平緩，而后置隔板的變形隨入射速度增加呈非線性增長趨勢。