高速弹体冲击侵彻混凝土靶数值模拟分析

高速弹体冲击侵彻混凝土靶数值模拟分析摘 要： 采用有限元程序LS-DYNA 对高速弹体侵彻混凝土靶进行了数值仿真研究，根据Hamilton 原理和有限元理论建立了由弹塑性动力学控制方程。

结果表明：接触冲击算法能较为真实反映高速弹体对混凝土靶的侵彻过程；不同速度及弹头长径比不同的弹体对混凝土靶的侵彻深度有较大影响；通过与经验公式的验算结果对比进行误差分析，证明数值分析所得结果可靠。

基于HJC 材料模型的损伤分析可以提取出靶体单元的损伤情况。

研究结果对我国高速弹体对混凝土靶体的侵彻设计具有重要参考价值。

关键词：哈密尔顿原理;动力学控制方程;侵彻；经验公式；损伤Numerical simulation of High-velocity projectilepenetrating concrete targetsAbstract:LS-DYNA finite element software was applied to simulate the concrete target penetrated by High-velocity rigid projectile ，According to Hamilton principle and finite element theory, the dynamical governing equation of non-conservative system consist of elastic-plastic material was built. The result show that: the process of High-velocity projectile penetrating concrete targets can be well demonstrated by the contact impact algorithm . The velocity and the Hr/d ratio of projectile significant effect penetrating depth ;there were a contrast between the results from the formulas and the result from Numerical simulation , prove the Numerical simulation are reliable .based on the Johnson-Holmquist-concrete constitutive relations and damage model ,concrete element damage can obtained. The simulation result provide some reference for the design of projectiles and targets.Key words :Hamilton principle ; dynamical governing equation ; penetration ;emprical formula ;damage1 引言高速弹体在军事上有着广泛的应用，它能够钻入地下，贯穿防护工事，进入目标内部发生爆炸。

弹体侵彻靶体的过载特性及侵彻深度关系到提高侵彻武器打击性能的研究和防护材料，防护结构的优化设计。

侵彻问题属于高速冲击动力学的研究范畴，是高应变率，大变形，非线性碰撞现象，其过程具有高速，高温，高压等基本特征。

进行侵彻的原型实验需要耗费大量的人力，物力，财力。

数值模拟因其经济性和高效性日益成为侵彻问题的重要研究手段。

本文采用大型非线性LS-DYNA 有限元程序，建立三维模型动态模拟了高速弹体侵彻混凝土靶体的侵彻过程。

模拟结果可为侵彻武器设计者和防护工事设计者提供很好的理论参考。

2 基于Hamilton 原理的弹塑性动力学控制方程动态有限元程序是以动力学控制方程为基础的。

本文将耗散能引入到Lagrange 函数中，根据Hamilton 原理推导出了考虑能量耗散的弹塑性体动力学控制方程的有限元格式。

首先取一单元体，设单元体的动能为T ，应变能为U ，外力势能为W e ,对于保守系统，Lagrange 函数可取为：e W U T L --= （1）对于弹塑性介质的非保守系统，设塑性耗散功为W d ，建立拉格朗日函数为d e L T U W W =--- （2）设u 为单元体中任一点的位移矢量，eu 为单元体上各节点的位移矢量，它是时间t 的函数。

令单元体中任一点的位移矢量u 用单元体上各节点的位移矢量eu 表示为:e Nu u = （3）式中N 为形函数矩阵，它是坐标z y x ,,的函数。

e u N u= （4） 于是可求得单元体的动能为()dV u N N u dV u u T eT T e vT v ρρ⎰⎰⎰⎰⎰⎰==2121 （5） 式中ρ为单元体积的质量。

按照几何关系，应变与节点位移之间的关系为e Bu =ε （6）式中B 为应变位移矩阵，它是几何矩阵，与时间t 无关。

则应变率与节点速率的关系为e Bu ε= （7） 对于弹塑性介质单元体的应变能为vU UdV =⎰⎰⎰ （8）式中：U 表示单位体积的应变能实际结构均存在着阻尼，计算结构在动载荷作用下的动态响应时，阻尼对结构动态响应的幅值和相位都有很大影响，正确的描述结构的阻尼对结构的动态响应分析至关重要。

结构的阻尼耗散可分为2个部分，（1）与质点的运动速度有关的耗散，其耗散力与质点的运动速度成正比，（2）与应变率有关的耗散对于阻尼耗散的第一部分，设单元体振动时，阻尼系数为c ，则单元体积上所受的阻尼力v f 为uc f v -= （9） 对于阻尼耗散的第二部分，与应变率有关的广义阻尼力r f 为εβ D f r = （10） 则单元体上阻尼力所消耗的能量为()()dV Nu N u c udV u c W eT T e vT v v ⎰⎰⎰⎰⎰⎰-=-=2121 （11）()dV u DB B u dV D W e VT T e V T r ⎰⎰⎰⎰⎰⎰-=-=βεβε2121 （12） 单元体上所受的外力分为两部分，即体积力),,(z y x V F F F F =和表面力),,(z y x S F F F F =。

它们的势能分别为21,e e W W 。

()dV F N u dV F u W V T TVe vV T e ⎰⎰⎰⎰⎰⎰-=-=1 （13）()dS F N u dS F u W SS T Te SS T e ⎰⎰⎰⎰-=-=2 （14）于是拉格朗日函数为()()()()12()T T T e T e e T e e T e vTe T T e T V S vSL u N Nu U c u N Nu u B DBu dV u N F dV uN F dSρβ⎡⎤=-++⎢⎥⎣⎦++⎰⎰⎰⎰⎰⎰⎰⎰ （15）由哈密尔顿原理[1-3]，将其在时间（t 1，t 2）上对L 积分，并使其变分等于0。

02121==⎰⎰t t t t Ldt Ldt δδ （16）Lagrange 函数中各项的变分可表达如下：()dV u N N uT ve T Te ⎰⎰⎰= ρδδ （17） T vvU UdV dV δδδεσ==⎰⎰⎰⎰⎰⎰ （18）()()()() =Te Te p v TTe Te e Te e p vvU uB D dV uBD Bu dV uBD dVδδεεδδε⎛⎫=- ⎪⎝⎭-⎰⎰⎰⎰⎰⎰⎰⎰⎰ （19）()dV Nu cN u W e T Tve v ⎰⎰⎰-= δδ （20） ()dV Bu D B uW e e T Tve r βδδ⎰⎰⎰-= （21） ()1Te T e V VW u N F dV δδ⎛⎫=- ⎪⎝⎭⎰⎰⎰ （22）()2Te Te S S W uN F dS δδ⎛⎫=- ⎪⎝⎭⎰⎰ （23）将表达式（17）~（23）代入表达式（16）采用分部积分，并有()01=t u e δ，()02=t u e δ可得()()()()()()()0][2121=⎪⎪⎭⎫ ⎝⎛+⎪⎪⎭⎫ ⎝⎛+-⎪⎪⎭⎫ ⎝⎛-⎪⎪⎭⎫ ⎝⎛+⎪⎪⎭⎫⎝⎛-⎪⎪⎭⎫ ⎝⎛-=⎰⎰⎰⎰⎰⎰⎰⎰⎰⎰⎰⎰⎰⎰⎰⎰⎰⎰⎰⎰⎰⎰dt dS F N u dV F N u dV u B D B u u dV N cN u u BdV D B u u BdV D B u u dV N N u Ldt S S T T e V V T T e ee T T v e e v T T e e v p T T e ev e T Te t t e vT Te t t δδβδδδδρδδ （24） 令,,,,eq eq eq eq p K M C F F 分别表示单元体的刚度矩阵，质量矩阵，阻尼矩阵，外载荷矩阵和由塑性变形产生的广义荷载矩阵则有T eq e vT p e p vT eq vT T eq e vvT T eq V s vSK B D BdVF B D dVM N NdVC cN NdV BD BdVF N F dV N F dSερβ====+=+⎰⎰⎰⎰⎰⎰⎰⎰⎰⎰⎰⎰⎰⎰⎰⎰⎰⎰⎰⎰ （25）由此得到()[]02121=++---=⎰⎰t t eq p e eq e eq e eq Te t t dt F F u C u K uM u Ldt δδ （26） 因为积分区间任意，所以被积函数()[]0=++---eq p e eq e eq e eq Te F F u C u K uM u δ （27） 由于单元体的变分()Teu δ任意，所以0=++---eq p e eq e eq e eq F F u C u K uM （28） 即e e e eq eq eq eq p M uC u K u F F ++=+ （29） 将式（29）写为增量形式的动力学控制方程：e e e eq eq eq eq p M uC u K u F F ∆+∆+∆=∆+∆ （30） 式（30）中，eq F ∆和p F ∆分别表示外载荷矩阵增量和由塑性变形产生的广义荷载矩阵增量,其表达式为：T p e p vF B D dV ε∆=∆⎰⎰⎰T T eq V s vSF N F dV N F dS ∆=∆+∆⎰⎰⎰⎰⎰动态有限元程序，是通过逐步时间积分来求解动力学方程的，需要选定合适的加载时间步长控制解法的稳定性，加载时间步长Δt 选得过小，将增加计算时间，加载时间步长Δt 选得过大将影响计算精度，对于冲击动态载荷，本文参考文献[4]中的公式确定时间步长为18lt λ∆∆≤（31） 式中：Δl 为结构中最小单元的特征尺寸，λ为纵波波速。