高速破片穿透液艙的數值模擬研究

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

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

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

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

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

1 CEL算法的基本方程

CEL方法同時具有拉格朗日方法和歐拉方法的優勢，較好的處理大變形、碰撞及一些流固耦合問題[20]。CEL在處理耦合問題中，應用拉格朗日算法計算結構以保持其在響應上的精度，用歐拉算法計算流體介質以解決其流動過程中的大變形問題。拉格朗日算法是一種離散化數值求解方法，其求解過程中，材料會依附在單元網格上，在網格的運動和變形的過程中伴隨著材料的流動，計算過程中，單元會發生變形但單元的初始質量不變。具體控制方程為

(1)

(2)

(3)

式中:ρ,u,e,σ分別為拉格朗日流體介質的密度、速度、比內能和應力。

歐拉算法是一種固定網格而使材料在網格間運輸的求解方法，計算過程中一個歐拉單元按用戶定義的材料分配，有著不同的物質，整個計算過程中空間網格將會固定，而材料網格將以拉格朗日單元的狀態變量實時傳輸到空間網格中，具體控制方程為

(4)

(5)

(6)

式中:ρ,u,e,σ,f分別為歐拉流體介質的密度、速度、比內能、應力和質量力。

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

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

2 高速破片穿透液艙的數值分析

高速破片穿透液艙的數值模型使用有限元仿真進行建立運算分析，符合工程應用背景以及科學研究需求。

2.1 液艙結構

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

仿真試驗采用J-C模型描述液艙的本構關系。J-C模型的關系式為

(7)

(8)

式中：θ為當前溫度；θmelt為熔化溫度；θtransition為轉化溫度。

在J-C失效模型中，斷裂應變定義式為

(9)

表1 結構材料表

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

圖3 破片入射面網格圖Fig.3 Detail of the entry wall mesh

圖4 壓力測點位置示意圖Fig.4 Sketch of pressure transducers

2.2 歐拉域

鑒于高速沖擊問題帶來的大變形以及流體的不規則流動，仿真試驗將用歐拉網格建立水和空氣的模型。在彈片穿透過程中，結構的大變形同時會伴隨著流體的運動，因此歐拉域的設置應大于結構的區域以保證流體在流動過程中不會消失，保證能量守恒，所以歐拉域設置為190 mm寬。水和空氣的材料屬性通過體積分數方式離散進歐拉域內，最終歐拉網格單元數為4 061 300個，如圖5所示。

圖5 歐拉網格展示圖Fig.5 Mesh of the fluids

水的狀態方程采用Mie-Grüneisen狀態方程，定義式為

(10)

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

表2 水和空氣的材料參數表

3 數值計算結果分析

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

3.1 破片高速侵徹過程

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

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

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

圖7給出四組工況下速度衰減曲線與理論分析解和Varas等實驗值得對比情況。可以看出數值分析方法CEL的結果與理論分析結果非常接近，只是CEL算法在彈片穿透液艙時會有較大的減速。由圖7看出，CEL算法在彈片速度衰減趨勢上與理論分析和實驗有較好的擬合；由于彈片在液艙內滯留時間極短，液體阻滯作用有限，同時液面距彈道位置較遠，液面效應較小，液艙內裝水量對速度衰減趨勢并沒有影響；彈片在穿透入射面時速度大致衰減15%左右。

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

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

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

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

圖8 初始速度600 m/s裝水量75%下空穴演化0.26 ms,0.46 ms,0.60 ms時刻實驗數據與CEL算法數據的對比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時刻實驗數據與CEL算法數據的對比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 不同入射速度對液艙隔板壓力和變形的影響

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

圖10 裝水量75%初始速度900 m/s的實驗值與數值的壓力時歷曲線對比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的實驗值與數值的壓力時歷曲線對比Fig.11 Experimental and numerical pressure time historyPtn and Ptf in a tube 75% filled impacted at 600 m/s

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

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

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

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

4 結 論

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

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

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

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