0 引言
冲击反应材料又称活性材料,指在冲击作用下发生化学反应而放热的一类特殊材料。因为冲击反应材料的反应特性而使其对于各种目标都具有良好的毁伤能力,所以基于冲击反应材料制成的毁伤元的侵彻特性一直是新型战斗部研究中的热门领域。
冲击反应材料可以为单质(铝热剂),也可以由多种物质经过烧结等工艺制成(Al/PTFE、AlNi等)。根据反应类型的不同,冲击反应材料大致分为铝热型、金属氧化型、金属合金化型3种[1]。
自1973年HughE[2]提出反应破片(reactive fragment)这一概念以来,国外学者就冲击反应材料的反应机理和毁伤特性都进行了较为深入的研究[3-6]。
国内相关研究起步较晚,但对于冲击反应材料毁伤特性的研究已经具有了一定成果:梁君夫[7]研究了活性破片撞击屏蔽盖板装药的引爆机理,分析了活性破片及盖板参数和引爆规律之间的定量关系,建立了活性破片引爆盖板装药的分析模型;黎勤基[8]在AutoDyn平台下建立了活性破片撞击双层靶的数值模型并探究了活性破片爆燃率和破片参数之间的关系;葛超[9]设计了一种以活性材料作为内芯的弹丸,并揭示了其对柴油油箱的引燃机理,证明了活性材料相对于传统惰性材料在毁伤性能上的优越性。
数值模拟是冲击动力学中主要的研究手段,然而在实际仿真中,完整的多层靶碰撞数值模型体量较为庞大,网格数量在几十万到几百万之间,在破片或弹丸的高速碰撞下往往会出现计算效率低下的问题,本研究中在LS-DYNA平台下利用材料模型的二次开发功能对经典Johnson Cook本构模型进行了温度项的优化,并加入化学反应机制得到冲击反应材料的自定义本构子程序,而后将自定义本构模型运用到单层防护结构冲击的数值模拟中,得到了不同形状破片对撞击速度对单层防护结构抗毁伤特性的影响规律。
1 本构模型的开发
1.1 Johnson Cook本构模型
在冲击动力学中,对于金属应用最广泛的材料模型就是JohnsonCook模型[10],其本构方程考虑了应变率、温度和应变对材料状态的影响:
(1)
式中:
为材料的初始屈服应力(参考应变率
和室温Tr条件下);B、n分别为同等条件下的应变硬化模量和硬化指数;m是和材料热软化行为相关的系数。
材料在冲击变形过程中温度逐渐升高,所以在数值计算过程中每个时间步下都需要对温度项进行更新,增加了计算的复杂程度,延长了计算时间。而在高应变率作用下,材料的绝热温升可归结为塑性做功量[11],以此得到修正后的Johnson Cook模型温升项:
(2)
式中: ρ为材料密度;Cp为材料的定压比热容; β为塑性功转热系数,在材料符合变形的均匀性假设前提下一般取0.9[12]。
1.2 返回映射与应力更新算法
在数值计算过程中,单元的应力更新往往是通过本构积分算法实现的。常用的返回映射算法分为两步:在试探步中假设材料在变形过程中始终处于弹性状态,此时的应力则由应变增量Δεn+1得到,称为试探应力
而在实际过程中材料的变形一般会偏离屈服面,因此在试探步后还要加入塑性步对屈服面下的应力状态进行修正,使其回到n+1时刻的屈服面上。
在向后欧拉格式中,塑性应变率
和内变量变化率
均取tn+1时刻的值,且要求更新后的应力状态
满足tn+1时刻的屈服条件,即[13]:
(3)
对于满足VonMises屈服条件的材料,在其变形过程中可以认为塑性应变rij和内变量qα不发生变化。
由此得到其应力更新格式为
sn+1=s*(n+1)-3GΔεp(n+1)
(4)
式(4)是非线性方程,所以在实际计算过程中可以用牛顿迭代法进行求解,在对应力进行塑性修正的过程中,单元的总应变
可视作常数,求解实际上是对Δλ也即塑性应变增量Δεp进行迭代计算的,因此由方程组(3)中的第二、三式得到迭代的目标函数为
(5)
将Johnson Cook本构的流动应力形式代入式中的σy,再把
作为自变量求导从而得到迭代函数为
(6)
Johnson Cook材料模型本构的具体编写流程如下:
1) 依据弹性状态下的应力应变关系求出n+1时刻的应力试探值
(7)
2) 构造屈服函数检查是否屈服
(8)
如果φ≤0,则说明材料的变形依然处于弹性阶段,变形中不含有塑性部分,此时
即真实应力等于试探应力,如果φ≥0,则表明材料变形进入塑性阶段,此时应通过修正有效塑性应变进而对材料的应力状态进行更新以使其满足屈服条件,此时通过上述式(6)进行迭代,如果φ和 f(k)同时小于给定的容差值则说明增量计算结果已收敛,停止迭代并完成塑性应变增量Δεp的修正,若不收敛,则在此步继续迭代。
3) 构造比例因子对偏应力部分进行修正,使之回到新的屈服面上
(9)
4) 更新n+1时刻的应力分量
(10)
1.3 状态方程
在爆炸冲击加载条件下,材料近似于流动状态,并且材料所受压力远远大于材料本身的强度,所以在这种条件下需要用状态方程进行材料压力的更新,状态方程在形式上一般为材料相对体积V和能量密度E的函数
P=P(V,E)
(11)
在材料发生压缩或膨胀行为的过程中,定义压缩系数为
(12)
而在弹塑性力学中,体应变的定义为材料在变后单位体积的变化,定义为
(13)
由此可知体应变和压缩系数的关系
(14)
Grunesien状态方程的一般形式为
(15)
若材料处于膨胀状态,则式(15)变为
p=ρ0c2μ+(γ0+aμ)E
(16)
S1、S2、S3为波后质点速度-冲击波速度曲线的斜率系数,α为对γ0的一阶体积修正。
根据上述推导出的体应变εv和压缩系数μ的关系代入Gruneisen方程中即可求出静水压力P。
1.4 失效模型
在强冲击和高应变率作用下,结构会发生失效和破坏,材料的损伤演化是个非常复杂的过程,不同种类的材料会产生不同的破坏形式,如剥落、冲塞、层裂、花瓣形失效等。一旦发生失效,材料的承载能力会迅速下降。
在本构算法中,对于材料的失效有2种处理方式:第一种是使失效单元只能承受压应力而不能承受拉应力,第二种是直接将失效单元删除。
为了完善材料本构,向本文中自定义材料模型中加入与LS-DYNA原有Johnson Cook本构相同的失效模型,即:
(17)
为材料的破坏应变,σ*为应力三轴度,下面再对单元损伤度进行定义:
(18)
当D值达到1时,单元失效,当失效单元的压力小于零时,将单元偏应力Sij设为零,并重新计算其内能。
1.5 LS-DYNA用户子程序的嵌入
由美国LLAL(lawrence livermore national laboratory)公司开发的LS-DYNA商用仿真软件被广泛用于爆炸冲击及侵彻等瞬态分析,为了适应用户需求,LS-DYNA中专门提供了相关子程序和链接库方便用户开发自定义本构材料模型,其中umat41-50是对应自定义本构的材料模型子程序编号,用户可在dyn21.f文件中添加或修改对应本构的积分算法,然后通过intelFortran对静态库重新编译生成lsdyna.exe程序,在提交计算时选择编译好的求解器并在K文件调用相关自定义材料进行计算[14]。
依照前文所推导的迭代形式,本文中所用的自定义本构的编写流程如图1所示。
图1 自定义本构模型计算流程
Fig.1 Calculation flow of custom constitutive model
在生成lsdyna可执行程序时,需要在静态链接库打开dyna21.f文件修改对应编号下自定义材料本构的代码,然后对整个静态链接库使用intelParallelStudio 2013进行重新编译,将生成的求解器主程序复制到LS-DYNA主目录并在提交计算时选取运行。
1.6 自定义本构的验证
通过上文中的自定义本构材料模型编写流程编译得到可执行的LS-DYNA主程序,下面将用本文中自定义材料模型与实验结果进行对比以验证二次开发材料模型的有效性。
以文献[15]中的Weldox 460E钢靶侵彻实验为依据,模型采用cm-g-μs单位制,具体弹靶几何参数如图2所示,弹体为圆柱平头弹,直径20 mm,弹体全长为80 mm,靶板为圆形靶板,直径0.5 m,厚度为12 mm,结构参数示意图和数值模型如图3所示。
图2 模型参数
Fig.2 Model parameters
图3 数值模型示意图
Fig.3 Schematic diagram of numerical model
全模型用三维拉格朗日网格划分,网格尺寸为1mm,定义侵蚀接触并使用自定义本构,相关材料参数[16-17]如表1、表2所示。
表1 靶板材料参数
Table 1 Material parameters of target plate
ρ/(g·cm-3)G/GPaA/GPaB/GPa7.85750.4900.807nCmT/KTr/K0.730.0120.941 800293
表2 弹丸材料参数
Table 2 Parameters of projectile materials
E/GPaνρ/(g·cm-3)σ0/MPaEt/MPa2040.337.851 90015 000
在LS-PrePost后处理中取1/4模型,计算结果如图4所示。
图4 自定义Johnson Cook本构的计算结果与实验对比
Fig.4 Comparison between the calculated results of the custom Johnson Cook constitutive and the experimental results
从图4中可以看出,靶板在弹丸冲击下由于发生绝热剪切而出现了冲塞现象并形成了塞块,这与文献[18]中描述的一致,也与实验中高速摄像机观察到的结果完全符合。
在进一步的验证中,根据试验数据设置了六种不同的工况,让弹丸以200～400 m/s的速度对靶板进行侵彻并评估弹丸在穿靶后的剩余速度。在表3中,v0为弹丸撞击靶板的初速,ve为试验下测定的弹丸剩余速度,vc为数值模拟计算出的弹丸剩余速度,对比结果如表3。
表3 工况设置及计算结果
Table 3 Working condition setting and calculation results
工况编号v0/(m·s-1)ve/(m·s-1)vc/(m·s-1)误差值/%1200.471.459.816.22224.7113.7108.44.73244.2132.6113.214.64285.4181.1172.44.85303.5199.7186.46.76399.5291.1284.62.2
从图5中可以看出,使用本研究中自定义本构模型的计算结果和实验数据吻合较好,最大误差为16%,在可接受的范围内,经过验证,该自定义本构模型可以较好地模拟出材料在冲击作用下的力学行为,并具有较高的准确性,说明该自定义本构模型可以用于侵彻冲击相关的数值计算。
图5 实验与仿真结果对比曲线
Fig.5 Comparison curve between experimental and simulation results
1.7 化学反应机制
在冲击反应材料受到冲击载荷时,由于温度升高、塑性变形增加等方式在材料颗粒间产生了能量的积聚,当累积能量达到一定阈值时,进而发生化学反应。本文中假设冲击反应材料的化学反应进程完全由温度控制,在反应速率方面,引入Arrhenius理论模型,即反应速率表示为[19]:
(19)
式中:k为化学反应速率常数;y为反应程度;t为反应持续时间。
对于固体反应类型,反应速率可表示为:
(20)
式中:Ru是摩尔气体常数;T为绝对温度;A和Ru分别是指前因子和表观活化能。
若考虑Avrami-Erofeev反应模型, f (y)可写为[20]:
f(y)=n(1-y)[-ln(1-y)](n-1)/n
(21)
根据Ortega A[21]的研究,由此得到了反应程度y和绝对温度T的一阶微分形式:
(22)
由式(22)可以求出反应程度y。而由化学反应引起的温升和压力的升高可以通过反应放出的热量来表示:
(23)
(24)
式中:QJ为单位质量的物质完全反应放出的热量;CP为定压比热; γ为Gruneisen系数,通常被认为是比容的函数。
再结合冲击压缩状态,便可得到考虑化学反应后的压力和温度:
(25)
(26)
其中,VH、TH、PH代表冲击压缩下的比容、温度和压力。
2 AlNi合金弹丸对单层靶毁伤特性分析
2.1 冲击反应弹丸和普通弹丸侵彻性能对比
在Johnson Cook本构模型的基础上,进一步加入1.5节推导出的化学反应机制,编写出考虑化学反应的冲击反应材料模型,添加了化学反应机制的本构模型计算流程如图6所示。
图6 AlNi合金自定义本构计算流程
Fig.6 AlNi alloy customized intrinsic calculation procedure
图7 弹丸撞击单层靶板结构有限元模型示意图
Fig.7 Schematic diagram of finite element model of projectile impact single-layer target structure
为了研究冲击反应材料和普通材料在侵彻性能上的不同,另建立单层靶模型如图7所示。其中靶板厚0.7 cm,材质为普通工具钢;弹丸形状为立方体,边长2.1 cm,分别由密度相近的AlNi冲击反应材料和普通工具钢制成,相关材料参数[14,16]如表4—表6,弹丸初速为1 000 m/s。计算结果如图8所示。
表4 工具钢材料参数
Table 4 Tool steel material parameters
ρ/(g·cm-3)G/GPaA/GPaB/GPa7.83770.7920.51nCmT/KTr/K0.260.0141.031 793294
表5 AlNi合金材料参数
Table 5 AlNi alloy material parameters
ρ/(g·cm-3)G/GPaA/GPaB/GPa7.800.300.9671.247nCmT/KTr/K0.660.0360.851 490294
(注:表中参数由中北大学任凯等提供。)
表6 AlNi合金化学反应参数
Table 6 Chemical reaction parameters of AlNi alloy
QJ/(kJ·g-1)Ea(kJ·mol-1)n16.1620.881.60
图8 AlNi合金弹丸和普通钢弹丸的侵彻过程
Fig.8 Penetration process of AlNi alloy projectile and ordinary steel projectile
从图8中可以看出,10 μs时,钢和AlNi合金开始嵌入靶板,靶板受力产生弯曲,由于此时的温度和压力还未达到AlNi合金的反应阈值,所以钢和AlNi合金的侵彻行为在这一时刻区别并不大。20 μs时由于立方体弹丸边缘的应力集中作用,钢制弹丸相关区域的单元逐渐失效,靶板有冲塞行为趋势,而这时可以明显看到AlNi合金开始发生反应,与靶板接触部分出现燃烧膨胀产生高压对靶板形成更大的挤压区域。30 μs时钢制弹丸冲击靶板形成了完整的塞块,AlNi合金弹丸的反应程度达到最大,其产生了剧烈的爆燃反应,靶板穿孔进一步扩大;40 μs时钢制弹丸完全脱离靶板并带动塞块一起向前运动,此时AlNi合金弹丸几乎完全反应,可以看到在强烈的冲击作用下靶板后方产生了碎片云。50 μs时2种工况下的靶板穿孔已趋于稳定,计算结束。
图9为2种材料的弹丸侵彻靶板穿孔直径对比,经过测量,钢制弹丸侵彻靶板形成的穿孔直径为3.07 cm,AlNi合金弹丸工况下穿孔直径最大处为6.89 cm,相对于钢制弹丸的穿孔直径增大了1倍左右,从穿孔形状看,钢制弹丸的穿孔形状接近于圆角矩形,也即与破片形状相似,而AlNi合金弹丸的穿孔直径表现为不规则孔洞,这是冲击反应材料在高速撞击时发生剧烈反应导致弹体迅速膨胀破裂进而制造出分布不规则的高温高压区域,此外,从其中还可以看出AlNi合金弹丸侵彻过后的靶板在穿孔周围形成了若干凹陷,这是冲击反应材料弹体碎裂飞溅出的破片四散开来又撞击到靶板上形成的,这说明了相对于普通金属,冲击反应材料在侵彻单层靶板时有更好的毁伤性能。
图9 2种材料的弹丸侵彻靶板穿孔直径对比
Fig.9 Comparison of the perforation diameters of the two materials
2.2 弹丸形状对冲击反应材料侵彻性能的影响
为了研究弹丸形状对冲击反应材料侵彻性能的影响,在立方体弹丸的基础上,再设计球形和圆柱形2种弹丸,其中球形弹丸半径1.3 cm,圆柱形弹丸底面半径1.14 cm,高为2.28 cm,3种弹丸质量相同,其他计算条件与前文设置相同。
计算结果3种弹丸的穿孔情况如图10所示。
图10 3种形状弹丸穿孔结果对比
Fig.10 Comparison of perforation results of three shaped projectiles
从图9中可以看出,3种形状的冲击反应材料弹丸侵彻所形成的穿孔在形状上相近,都表现为不规则的圆形,穿孔周围都分布有若干大小不一的凹坑。在穿孔直径上,根据表7可知3种工况都区别不大,原因是冲击反应材料弹丸在撞击靶板时由于发生反应而破裂,其破裂形式并没有固定规律,所以弹丸形状和最终的破孔形式也没有联系。以上结果说明,对于冲击反应材料来说,弹丸形状对侵彻效果并没有明显影响。
表7 3种形状弹丸穿孔直径统计
Table 7 Statistics of perforation diameters of projectile with three shapes
弹丸形状穿孔直径/cm立方体6.89球形6.61圆柱形6.78
3 结论
本研究中采用LS-DYNA平台下的二次开发功能编写了冲击反应材料模型,并以AlNi合金为例进行了单层防护结构撞击特性研究,得到以下结果:
1) 在Johnson Cook本构的基础上修改了温度项并加入了冲击反应材料的化学反应机制,开发出了适用于冲击反应材料侵彻计算的材料模型。
2) 分别用普通材料和冲击反应材料作为弹丸开展了单层钢靶侵彻仿真计算,体现了冲击反应材料相对于普通惰性材料在侵彻性能上的优势。
3) 研究了不同形状冲击反应材料弹丸撞击靶板的穿孔特性,得出弹丸形状对于冲击反应材料在单层靶工况下的侵彻性能并没有明显影响。
[1] 张先锋,赵晓宁.多功能含能结构材料研究进展[J].含能材料,2009,17(6):731-739.ZHANG Xianfeng,ZHAO Xiaoning.Review on multifunctional energetic structural materials[J].Chinese Journal of Energetic Materials,2009,17(6):731-739.
[2] MONTGOMERY H E.Reactive fragment[P].USA Patent:19730375246,1976-06-08.
[3] THADHANA N,CHEN E.Shock synthesis of materials workshop held in atlanta.Georgia on May 24-26,1995[Z].1995.
[4] MERZHANOV A G.Combustion and explosion processes in physical chemistry and technology of inorganic materials[J].Cheminform,2003,34(40):323-345.
[5] FERRANTI L,THADHANI N N.Dynamic impact characterization of Al+Fe2O3+30% epoxy composites using time synchronized high-speed camera and visar measurements[J].Mrs Proceedings,2005,896:0896-H06-05.
[6] TOMAR V,ZHOU M.Molecular Dynamics simulations of shock-induced thermite reaction[J].Materials Science Forum,2004,465-466:157-162.
[7] 梁君夫.活性破片作用屏蔽装药引爆增强效应研究[D].北京:北京理工大学,2016.LIANG Junfu.Research on enhanced initiation behavior of reactive material fragment impacting covered explosive[D].Beijing:Beijing Institute of Technology,2016.
[8] 黎勤.活性破片碰撞双层靶毁伤效应研究[D].北京:北京理工大学,2016.LI Qin.Research on damage effect of double-spaced plates by reactive material fragment impact[D].Beijing:Beijing Institute of Technology,2016.
[9] 葛超,余庆波,卢冠成,等.活性芯体子弹对柴油油箱引燃效应及机理研究[J].北京理工大学学报,2020,40(10):1072-1080,1087.GE Chao,YU Qingbo,LU Guancheng,et al.Igniting effects and mechanism of diesel oil tank by projectile with reactive core[J].Transactions of Beijing Institute of Technology,2020,40(10):1072-1080,1087.
[10] JOHNSON G R,COOK W H.A constitutive model and data for metals subjected to large strains,high strain rates and high temperatures[J].Engineering Fracture Mechanics,1983,21:541-548.
[11] MEYERS M A,SUBHASH G,KAD B K,et al.Prasad.Evolution of microstructure and shear-band formation in α-hcp titanium[J].Mechanics of Materials,1994,17(2/3):175-193.
[12] MACDOUGALL D.Determination of the plastic work converted to heat using radiometry[J].Experimental mechanics,2000,40(3):298-306.
[13] 周旭,张雄.物质点法数值仿真(软件)系统及应用[M].北京:国防工业出版社,2015.
[14] 白金泽.LS-DYNA 3D理论基础与实例分析[M].北京:科学出版社,2005.
[15] BORVIK T,LANGSETH M,HOPPERSTAD O S.Perforation of 12 mm thick steel plates by 20 mm diameter projectiles with flat,hemispherical and conical noses[J].International Journal of Impact Engineering,2002,27(1):19-35.
[16] BORVIK T,LANGSETH M,HOPPERSTAD O S,et al.Ballistic penetration of steel plates[J].International Journal of Impact Engineering,1999,22(9/10):855-886.
[17] B∅RVIK,TORE,HOPPERSTAD O S,BERSTAD T,et al.Computational model of viscoplasticity and ductile damage for projectile impact[J].European Journal of Mechanics-A/Solids,2001,20(5):685-712.
[18] JIANG D,LI Y C.Numerical simulation of adiabatic shear in blunt nose projectile plugging target[J].Explosion and Shock Waves,2011,31(1):1-5.
[19] WANG H,PRING A,NGOTHAI Y,et al.A low-temperature kinetic study of the exsolution of pentlandite from the monosulfide solid solution using a refined Avrami method[J].Geochimica Et Cosmochimica Acta,2005,69(2):415-425.
[20] FALEIROS A C,RABELO T N,THIM G P,et al.Kinetics of phase change[J].Materials Research,2000,3(3):51-60.
[21] ORTEGA A,MAQUEDA L P,CRIADO J M.The problem of discerning Avrami-Erofeev kinetic models from the new controlled rate thermal analysis with constant acceleration of the transformation[J].Thermochimica Acta,1995,254(8):147-152.
[22] 任柯融.Al/Ni基含能结构材料冲击释能行为实验与数值模拟研究[D].长沙:国防科技大学,2018.