強洪夫,范樹佳,陳福振,劉 虎

(火箭軍工程大學動力工程系,陜西 西安 710025)

基于SPH方法的聚能射流侵徹混凝土靶板數(shù)值模擬*

強洪夫,范樹佳,陳福振,劉 虎

(火箭軍工程大學動力工程系,陜西 西安 710025)

在完全變光滑長度SPH(smoothed particle hydrodynamics)方法的基礎(chǔ)上，利用F.Ott等提出的修正SPH方法處理在求解多介質(zhì)大密度問題時的數(shù)值不穩(wěn)定性問題，運用Holmquist-Johnson-Cook本構(gòu)模型處理混凝土在沖擊載荷下的變形和損傷問題，對聚能裝藥射流侵徹混凝土靶板的過程進行了數(shù)值模擬，同時利用LS-DYNA非線性有限元程序進行對比，分析了2種方法得到的混凝土von Mises應力變化、射流頭部特定節(jié)點處的速度變化及裂紋演變，驗證了SPH方法的準確性。分析了另外2種不同尺寸的靶板在射流侵徹作用下的破壞形式，結(jié)果符合射流侵徹物理規(guī)律，表明該方法適合模擬聚爆炸與沖擊等大變形破壞等問題。

爆炸力學;侵徹;光滑粒子動力學;混凝土;破壞

二級串聯(lián)隨進戰(zhàn)斗部[1]作為攻擊跑道目標的有效手段，在現(xiàn)在戰(zhàn)爭中發(fā)揮著重要作用。其第1級戰(zhàn)斗部即采用聚能裝藥爆破后形成高速射流，對混凝土目標侵徹形成孔道，以便于第2級戰(zhàn)斗部跟進，繼續(xù)侵徹一定深度后爆破，從而對目標造成不可恢復的破壞。聚能射流侵徹混凝土的主要特點是：炸藥能量高，聚能效應明顯，形成高速射流,混凝土在高速沖擊作用下由于其拉壓不等的特性其破壞的模式復雜等。因此，對于該過程的數(shù)值模擬對于串聯(lián)隨進戰(zhàn)斗部的研究與設(shè)計具有重要意義。

傳統(tǒng)用于計算該過程的方法主要是基于現(xiàn)有商業(yè)軟件中的成熟算法，如多物質(zhì)Euler法、Lagrange法以及ALE(arbitrary Lagrangian-Eulerian)法等，這些方法大都采用有限差分法(finite difference method,FDM)、有限元法(finite element method,FEM)或有限體積法(finite volume method,FVM)等基于網(wǎng)格的數(shù)值方法進行離散求解，計算中不可避免地出現(xiàn)網(wǎng)格的扭曲和纏繞(針對Lagrange網(wǎng)格離散)或界面追蹤復雜、精度較低(針對Euler網(wǎng)格離散)的問題。光滑粒子流體動力學(smoothed particle hydrodynamics, SPH)方法作為一種無網(wǎng)格粒子方法，在對射流進行彈塑性流體動力學計算，對混凝土進行大應變、高應變率的變形計算時,可避免網(wǎng)格重分及算法耦合，因此非常適合此類問題的求解。最早采用SPH方法模擬爆炸的可追溯到J.W.Swegle等[2]，而L.D.Libersky等[3-4]最先將SPH方法運用于高速沖擊領(lǐng)域，隨后，M.B.Liu等[5-6]采用SPH方法模擬了聚能裝藥的爆轟過程，Qiang Hongfu等[7]運用F.Ott等[8]提出的修正SPH方法對聚能裝藥射流過程進行模擬，分析了不同的起爆方式對射流的影響。本文中，擬在強洪夫等[9-10]提出的完全變光滑長度SPH方法的基礎(chǔ)上，結(jié)合修正SPH方法[8]，處理爆炸和沖擊過程中密度和光滑長度變化劇烈的問題以及多介質(zhì)界面問題。采用SPH方法計算混凝土損傷大多引入Johnson-Cook[11]、SCG[12]等本構(gòu)模型，與混凝土的特性不符，未能捕捉混凝土在高速沖擊下出現(xiàn)的裂紋擴展等細節(jié)，且采用SPH方法計算包括炸藥的爆轟、藥罩的擠壓、射流的形成及發(fā)展、射流高速侵徹混凝土、混凝土損傷破碎等過程在內(nèi)的復雜全過程問題未見報道。因此，本文中將Holmquist-Johnson-Cook本構(gòu)模型[13]引入SPH方法，處理混凝土在高速沖擊作用下的變形損傷問題，對聚能炸藥爆炸擠壓形成射流，進而侵徹混凝土靶板的整個過程進行數(shù)值模擬，并對侵徹過程中射流頭部特定點處的速度變化進行分析，采用傳統(tǒng)網(wǎng)格方法進行對比驗證。同時為更充分地描述混凝土靶板在高速射流侵徹作用下的破壞損傷效應，對另外2種不同厚度的靶板進行設(shè)計，并開展數(shù)值實驗，探討SPH方法用于模擬諸如聚能射流侵徹混凝土等涉及爆炸與沖擊大變形多介質(zhì)問題的適用性。

1 SPH基本方程

在SPH方法中，通過對問題域的粒子離散來求解系統(tǒng)的狀態(tài)，系統(tǒng)中物理量f(r)及其導數(shù)·f(r)的SPH形式為：

f(ri)

(1)

(2)

式中：m、ρ、r分別為粒子的質(zhì)量、密度和位置矢量;Wij=W(ri-rj,h)為核函數(shù)，h為光滑長度，通常選用三次樣條核函數(shù)[14]。

為更好地模擬爆炸與沖擊、大變形大扭曲等密度和光滑長度變化劇烈的問題，本文中采用完全變光滑長度SPH方法[7]，同時為了很好地解決藥型罩和爆轟氣體間密度梯度較大帶來的間斷面不穩(wěn)定的問題，采用修正SPH方法[10]對聚能裝藥射流過程進行模擬，結(jié)合后的方程組如下：

(3)

式中:Wij=W(xi-xj,h)為核函數(shù)，它的選取直接影響計算的誤差和穩(wěn)定性，通常選用三次樣條核函數(shù)；iWij為核函數(shù)對粒子i坐標的空間導數(shù)，vij=vi-vj；σ為總應力張量；Πij為人工黏度；h為插值核寬度的一種度量，稱為光滑長度，表示W(wǎng)不顯著為零的取值范圍，控制著SPH粒子的影響域，通常設(shè)定時，W=0,這里h為空間和時間的函數(shù)，在連續(xù)方程中將其對時間求導，即：

(4)

式中:d為空間維數(shù)。對于光滑長度變化率dhi/dt與密度變化率dρi/dt之間相互耦合而帶來的難以顯式求解的問題，本文中采用迭代法[9]求解密度方程和光滑長度。

在計算過程中，為使粒子分布更有序，消除由于分布不均勻引來的粒子非物理聚集的問題，采用J.J.Monaghan[14]提出的XSPH對速度場進行光滑:

(5)

式中:ε(0≤ε≤1 )為常數(shù)，通過施加鄰近粒子的影響使自身的運動速度與鄰近粒子的平均速度相近，一般ε取0.3。

2 炸藥及藥型罩材料模型

本文模擬中，炸藥選用TNT，其爆轟速度由實驗[15]測得為6 930 m/s。對于爆轟氣體,使用標準JWL狀態(tài)方程[16]得到氣體壓力:

(6)

式中:ρ0=1 630 kg/m3為TNT炸藥的初始密度,v=ρ0/ρ,ρ為爆轟氣體的密度,e為爆轟氣體的比內(nèi)能;A、B、R1、R2和w為由實驗結(jié)果擬合得到的系數(shù)[17],且A=371.2 GPa,B=3.231 GPa,R1=4.15,R2=0.95,w=0.3。

對于藥型罩材料，選用銅，狀態(tài)方程選用Mie-Grüneisen狀態(tài)方程[16]:

(7)

(8)

式中:ρ0=8 530 kg/m3為銅的初始密度，Γ=1.99，cs=3 940 m/s,Ss=1.489。

3 混凝土本構(gòu)模型及狀態(tài)方程

圖1 混凝土屈服強度模型曲線Fig.1 Curves for yield strength model of concrete

圖2 混凝土累計損傷模型曲線Fig.2 Curve for cumulative damage model of concrete

HJC本構(gòu)模型[13]分為強度模型、累計損傷模型和狀態(tài)方程3部分，模型曲線分別如圖1～3所示。

3.1 屈服強度模型

歸一化的屈服強度函數(shù)為：

(9)

3.2 累計損傷模型

損傷值為由等效塑性應變和塑性體積應變引起的損傷累積起來的形式：

(10)

(11)

式中:D1和D2為混凝土的損傷常數(shù)。累計損傷模型曲線如圖2所示。

3.3 狀態(tài)方程

圖3 混凝土狀態(tài)方程曲線Fig.3 Curves for equation of state of concrete

混凝土靶板在沖擊載荷下的壓縮變形分為彈性階段、過渡階段和壓實階段3個階段，狀態(tài)方程曲線如圖3所示，圖中的plock和pcrush分別為下式中的pl和pc,μlock、μcrush和μp,lock分別為下式中的μl、μc和μpl。

(1)彈性階段(0≤μ≤μc)

在這一階段，混凝土發(fā)生可以恢復的彈性變形，加載與卸載的方程相同，為：

p=Keμ

(12)

式中:彈性體積模量Ke=pc/μc，pc和μc分別為單軸強度抗壓實驗中得到的壓碎壓力和壓碎體積應變，μ=ρ/ρ0-1為單元體積應變，ρ和ρ0分別為單元的密度和初始密度。

(2)過渡階段(μc<μ≤μpl)

在這一階段中混凝土中的氣體逐漸被擠壓出去，混凝土產(chǎn)生破碎性斷裂，加載時:

p=pc+Kp(μ-μc)

(13)

式中:Kp=(pl-pc)/(μpl-μc)為塑性體積模量，pl為壓實壓力;μpl=ρg/ρ0-1為壓實體積應變，對應于pl處的體積應變,ρg為顆粒密度。

該區(qū)的卸載是通過相鄰區(qū)域間插值的一條路徑進行的，卸載的狀態(tài)方程為：

p=[(1-F)Ke+FKp]μ

(14)

式中:插值因子F=(μmax-μc)/(μpl-μc) ，μmax為卸載前達到的最大體積應變。

(3)壓實階段(μ≥μpl)

該階段混凝土完全被壓碎，加載的狀態(tài)方程為:

(15)

卸載的狀態(tài)方程為：

p=K1μ

(16)

混凝土所能承受的最大拉伸強度為T(1-D),狀態(tài)方程曲線如圖3所示。

4 聚能射流侵徹混凝土靶板的數(shù)值模擬

數(shù)值模擬中采用的聚能裝藥模型及混凝土模型如圖4(a)所示，藥柱寬度為40.00 mm，裝藥頭長度為20.75 mm，藥孔長度為42.00 mm，張角為60 °,藥型罩厚度為3.00 mm，炸高為40.00 mm，混凝土靶板尺寸為100 mm×30 mm。炸藥采用點起爆的方式，起爆點為(0,0.07),具體粒子配置如圖4(b)所示,粒子間距均為0.5 mm,其中TNT炸藥粒子數(shù)為26 766,藥型罩粒子數(shù)為4 480,混凝土粒子數(shù)為12 000。光滑長度取1.5倍粒子間距，時間積分采用蛙跳格式，時間步長為0.1 μs。

本文中利用LS-DYNA非線性有限元程序進行對比驗證。計算模型由炸藥、藥型罩、空氣和混凝土靶板4部分組成。其中，對炸藥、藥型罩和空氣采用Euler網(wǎng)格建模，對混凝土靶板采用Lagrange網(wǎng)格建模，在混凝土靶板與空氣和藥型罩材料間采用耦合算法。聚能裝藥是線性的，考慮到模型的對稱性，為了減少計算時間，將模型簡化為二維平面對稱問題并建立其1/2模型。在對稱面上施加對稱邊界約束，在混凝土靶板邊界施加固定約束，把靶板側(cè)面設(shè)定為無反射邊界。對炸藥、藥型罩、空氣和混凝土靶板均選用與SPH數(shù)值模擬中參數(shù)相同的本構(gòu)模型和狀態(tài)方程。

圖4 算例模型結(jié)構(gòu)Fig.4 Construction of the model

圖5 聚能射流形成過程對比Fig.5 Comparison of the formation of jet flows

圖5給出了采用SPH方法和LS-DYNA程序計算過程中炸藥爆轟及藥罩被擠壓、射流形成過程的對比。圖6給出了采用SPH方法和LS-DYNA程序計算過程中射流侵徹混凝土形成開孔、混凝土損傷及混凝土完全破壞的過程的對比。可以看出，炸藥爆轟波逐漸將藥型罩向中軸線方向擠壓形成射流和杵體，隨后混凝土在射流的侵徹下發(fā)生壓縮和剪切變形，產(chǎn)生剝落，靶板彈坑的上部呈漏斗狀。隨著射流的不斷侵徹，靶板背面首先由于反射波到達壁面后的反射拉伸作用，導致背面出現(xiàn)拉伸損傷，部分脫落，開孔周圍則逐漸出現(xiàn)不同程度的層裂，并且裂紋逐漸向外擴展，最終形成一個漏斗狀損傷區(qū)域，中間部分出現(xiàn)粗細均勻的孔道，裂紋形態(tài)與圖7中的實驗結(jié)果基本相符。通過對比驗證，發(fā)現(xiàn)利用2種數(shù)值模擬方法得到的結(jié)果均符合射流侵徹物理規(guī)律，但LS-DYNA程序通過定義失效刪除單元產(chǎn)生的裂紋非常規(guī)則，數(shù)量少，基本沒有在裂紋上產(chǎn)生小裂紋的狀況，且不能模擬侵徹后碎片的物理特性，裂紋模擬效果受到網(wǎng)格劃分大小的影響，SPH方法則恰恰克服了以上困難。

由圖6可知，在侵徹過程中，前2段射流發(fā)揮了主要作用。因此，分別選取SPH方法及LS-DYNA程序計算結(jié)果中前2段射流的頭部頂點粒子A、B進行分析。圖8為2種方法的計算結(jié)果中A點和B點粒子速度隨時間的變化曲線，由此可知，藥型罩在0.02 ms時受到炸藥爆轟波的擠壓開始逐漸加速，在約0.10 ms時，在SPH方法的計算結(jié)果中A點粒子速度達到最大值3 247.9 m/s，在LS-DYNA程序的計算結(jié)果中A點粒子速度達到最大值3 457.4 m/s。射流頭部在0.14 ms時開始與靶板接觸，速度大幅降低，造成混凝土的損傷，隨著侵徹的深入，速度逐漸減小為零。2種方法的不同之處在于，SPH方法的計算結(jié)果中第1段射流在0.20 ms時速度基本降為零的同時第2段射流開始侵徹，而在LS-DYNA程序的計算結(jié)果中，第1段射流在0.14 ms時與混凝土靶板接觸，速度逐步降低，在0.17 ms左右時，第2段追上第1段射流，2段射流合為一段繼續(xù)侵徹，直到速度基本降至零。

圖6 聚能射流侵徹過程中不同時刻，混凝土靶板中Von Mises應力的分布Fig.6 Von Mises stress distribution in concrete target penetrated by shaped charge jet at different times

圖7 聚能射流侵徹的混凝土靶Fig.7 Concrete target penetrated by shaped charge jet

圖8 射流頭部特定點處速度-時間曲線Fig.8 Velocity-time curves at special points of jet head

相對于LS-DYNA中通過刪除網(wǎng)格單元形成裂紋的處理方法，SPH不刪除粒子單元而是通過顯示混凝土靶板的損傷值D來表述材料失效形成的裂紋，損傷后的SPH粒子按無黏性的流體粒子進行計算，可以完整地再現(xiàn)裂紋的損傷演化過程，損傷狀態(tài)及裂紋擴展分布情況如圖9所示，計算得到的靶板最終損傷狀態(tài)如圖10所示。由圖10可知：采用LS-DYNA得到的計算結(jié)果中彈坑直徑與采用SPH方法得到的計算結(jié)果中彈坑直徑差別較大，其主要原因在于LS-DYNA中理想空氣與混凝土靶板的重疊區(qū)域?。徊捎肧PH方法得到的計算結(jié)果中靶板左右壁面出現(xiàn)裂紋，而采用LS-DYNA得到的計算結(jié)果中沒有出現(xiàn)，其主要原因在于SPH中靶板的左右邊界是自由邊界，沖擊波到達左右壁面時由于反射波的作用開始出現(xiàn)拉伸損傷，而LS-DYNA靶板定義的是無反射邊界，反射波不能到達左右壁面。

圖9 在聚能射流的侵徹下不同時刻混凝土靶板的損傷情況Fig.9 Damage in concrete target penetrated by shaped charge jet at different times

圖10 不同方法模擬得到的裂紋長度及分布情況Fig.10 Length and distribution of crack simulated by different methods

為了更充分地描述混凝土靶板在高速射流侵徹作用下的破壞損傷效應，采用上述SPH方法對另2種尺寸的靶板進行了數(shù)值實驗，粒子數(shù)分別為4 000和40 000，計算結(jié)果見圖11，圖12為選取的射流頭部特定點進行追蹤得到的速度-時間曲線。從圖11～12可以看出，當靶板尺寸較小時，在較短的時間內(nèi)只需射流頭部一段即可將靶板完全擊穿，射流孔處損傷均勻；靶板較厚時，需要3段或更多段射流才可將靶板擊穿，而且在射流孔處不會形成大面積的損傷，而是形成不同程度的損傷裂紋延伸至靶板內(nèi)。

圖11 另外2種不同尺寸的混凝土靶板在聚能射流侵徹作用下的損傷Fig.11 Damage in other two concrete targets with different sizes penetrated by shaped charge jets

圖12 另外2種不同尺寸的靶板聚能射流侵徹作用下特定節(jié)點的速度-時間曲線Fig.12 Velocity-time curves at special points of other two concrete targets with different sizespenetrated by shaped charge jets

5 結(jié) 論

運用無網(wǎng)格SPH方法計算得到了聚能裝藥射流形成、發(fā)展及高速沖擊混凝土整個過程,并與傳統(tǒng)網(wǎng)格進行了對比驗證，得知：

(1)采用SPH方法能有效解決在高速沖擊、大變形問題的計算中不可避免地出現(xiàn)Euler網(wǎng)格界面追蹤復雜或Lagrange網(wǎng)格扭曲和纏繞的問題?；谕耆児饣L度SPH方法與修正SPH方法相結(jié)合的方法，有效解決了爆炸與沖擊中光滑長度變化劇烈的問題以及多介質(zhì)界面上大密度間斷帶來的計算不穩(wěn)定問題。

(2)運用HJC本構(gòu)模型和狀態(tài)方程描述混凝土在高速沖擊下的變形損傷特性，不僅得到了破孔的形狀，而且對由于壁面反射波造成的靶板背面拉伸損傷及破孔周圍裂紋的擴展等變形損傷細節(jié)描述精確，與實際物理現(xiàn)象符合。采用2種方法得到的計算結(jié)果基本一致，驗證了SPH方法的準確性。同時，采用SPH方法計算得到的裂紋形態(tài)與實驗結(jié)果吻合更好，更合理。

(3)隨著靶板尺寸的增大，射流孔徑逐漸減小，裂紋數(shù)目逐漸增加，裂紋的細度及延伸的長度逐漸增加，破壞的能量逐漸向周圍擴散。

[1] 段建,楊黔龍,周剛,等.串聯(lián)隨進戰(zhàn)斗部侵徹混凝土靶實驗研究[J].爆炸與沖擊,2007,27(4):364-369. Duan Jian, Yang Qianlong, Zhou Gang, et al. Experimental studies of a tandem follow-through warhead penetrating concrete target[J]. Explosion and Shock Waves, 2007,27(4):364-369.

[2] 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.

[3] Libersky L D, Petscheck A G. Smoothed particle hydrodynamics with strength of materials[C]∥Trease H, Fritts J, Crowley W. Proceedings of the Next Free Lagrange Conference. NY: SpringerVerlag, 1991,395:248-257.

[4] Libersky L D, Petscheck 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.

[5] 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.

[6] 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.

[7] Qiang Hongfu, Wang Kunpeng, Gao Weiran. Numerical simulation of shaped charge jet using multi-phase SPH method[J]. Transactions of Tianjin University, 2008,14(1):495-499.

[8] Ott F, Schnetter E. A modified SPH approach for fluids with large density differences[J]. Arxiv Physics E-prints, 2003:3112.

[9] 強洪夫,高巍然.完全變光滑長度SPH法及其實現(xiàn)[J].計算物理,2008,25(5):569-575. Qiang Hongfu, Gao Weiran. SPH method with fully variable smoothing lengths and implementation[J]. Chinese Journal of Computational Physics, 2008,25(5):569-575.

[10] 強洪夫,王坤鵬,高巍然.基于完全變光滑長度SPH方法的HE爆轟過程的數(shù)值試驗[J].含能材料,2009,17(1):27-31. Qiang Hongfu, Wang Kunpeng, Gao Weiran. Numerical study of high explosive detonation process using SPH method with fully variable smoothing lengths[J]. Chinese Journal of Energetic Materials, 2009,17(1):27-31.

[11] Johnson G R, Cook W H. A constitutive model and data for metals subjected to large strains, high strain rates and high temperatures[C]∥Proceedings of the Seventh International Symposium on Ballistics. Hague, Netherlands, 1983:571-574.

[12] Steinberg D J, Cochran S G, Guinan M W. A constitutive model for metals applicable at high strain rate[J]. Journal of Applied Physics, 1980,51(3):1498-1504.

[13] Holmquist T J, Johnson G R, Cook W H. A computational constitutive model for concrete subjected to large strains, high strain rates, and high pressures[C]∥Proceedings of the 14th International Symposium on Ballistics. Quebec, Canada, 1993:591-600.

[14] Monaghan J J. Smoothed particle hydrodynamics[J]. Reports on Progress in Physics, 2005,68(8):1703-1759.

[15] Liu G R, Liu M B.光滑粒子流體動力學:一種無網(wǎng)格粒子法[M].韓旭,楊剛,強洪夫,譯.長沙:湖南大學出版社,2005:195-197.

[16] Livermore Software Technology Corporation. LS-DYNA keyword user’s manual[M]. Livermore: Livermore Software Technology Corporation, 2012:17-45.

[17] Liu M B, Liu G R, Zong Z, et al. Computer simulation of high explosive explosion using smoothed particle hydrodynamics methodology[J]. Computers and Fluids, 2003,32(3):305-322.

(責任編輯 張凌云)

Numerical simulation on penetration of concrete target by shaped charge jet with SPH method

Qiang Hongfu, Fan Shujia, Chen Fuzhen, Liu Hu

(PowerEngineeringDepartment,Xi’anHi-TechInstitute,Xi’an710025,Shaanxi,China)

On the basis of the smoothed particle hydrodynamics (SPH) method with fully variable smoothing lengths, the modified Ott-Schnetter SPH method was used to cope with the computational instability in the large density gradient problems. The Holmquist-Johnson-Cook constitutive model was used to cope with the deformation and damage of concrete under impact. By comparing with the corresponding ones simulated by the non-linear finite element program LS-DYNA, the changing courses of the von-Mises stress and cracks in the concrete targets, and the velocities at the some special points of the jet head were analyzed, which proved the feasibility and accuracy of the SPH method. Other two concrete plates with different sizes penetrated by shaped charge jets were also simulated. The results are in good agreement with the physical principle. So this method can be used to deal with the multi-material and large deformation problems such as detonation and impact.

mechanics of explosion; penetration; smoothed particle hydrodynamics; concrete; damage

10.11883/1001-1455(2016)04-0516-09

2014-12-03；

2015-03-24

國家自然科學基金項目(51276192);國家重點基礎(chǔ)研究發(fā)展計劃(973計劃)基金項目(61338); 火箭軍工程大學創(chuàng)新型基金項目(EPXY0806)

強洪夫(1963— )，男，博士，教授，博士生導師;

范樹佳，fan_shu_jia@163.com。

O389國標學科代碼：1303530

A