非线性侵彻动力过程的再生核质点法*

赵光明,宋顺成

(1.安徽理工大学煤矿安全高效开采省部共建教育部重点实验室,安徽 淮南232001;2.西南交通大学力学与工程学院,四川 成都610031)

1 引 言

动力侵彻过程是力学研究中的一个重要课题,它与军事、科技、国民经济的发展密切相关。因为侵彻实验速度有限、费用比较高,数值模拟成为侵彻过程研究的有效方法。由于侵彻过程涉及高温、高压和高速度,压力梯度、温度梯度、位移速度梯度很大,侵彻过程中材料破坏涉及冲塞、侵彻、层裂、崩落、溅飞等各种不同复杂形式,因此侵彻过程研究及其数值模拟是动力学中的难点[1-2]。目前常用的方法有拉格朗日(Lagrange)、欧拉(Euler)、拉格朗日-欧拉耦合(A LE)以及光滑质点流体动力学(SPH)法等。欧拉法难以跟踪界面位置,计算时间长,强度(失效)状态和位移历程关系计算精度差等;拉氏法当网格变形严重时需要网格重分,导致计算精度下降,甚至计算失败。而SPH 法更适合超高速撞击等大变形、高应变率现象的数值模拟,是一种无网格的拉格朗日法,比其他方法具有更大的实用性[3-5]。W.K.Liu等[6-7]的研究表明,SPH 法不能完全满足一致性条件,而且采用此法偏微分方程时,边界点的解不稳定,他们通入修正函数,对SPH 方法进行改进,提出了新型的再生核质点法(reproducing kernel particle method,RKPM)。

侵彻过程的数值模拟,在很大程度上取决于材料的本构模型。本文中,通过引入Johnson-Cook 本构模型及损伤模型,采用RKPM 无网格法实现高速侵彻过程数值模拟,同理论解的结果进行比较,验证数值模拟方法的有效性。

2 RKPM 形函数

由于SPH 对于边界点以及不规则离散点不能精确求解,W.K.Liu 等[6-7]对SPH 方法作了改进,提出了再生核质点法。在SPH 方法[3]中,任意函数可以近似为

式中:x 是空间点,s 是领域内的点。在Dirac Delta 条件下,很难解决离散化问题,为此,进行修改[6]

式中:uα(x)称为u(x)的再生函数,α是膨胀系数,修正核函数¯φα可以表示为

式中:C(x-s)是修正函数。通过修正,RKPM 方法可以避免SPH 方法在区域的适用性

式中

根据Taylor 级数展开式(2),可以得到

所以

式中

将再生方程离散为

式中:xi 表示质点,n 是计算域内质点数量,ΔV i 表示质点x i 的体积,Ψi 称为再生核质点法的形函数。

3 侵彻材料的本构关系

根据Johnson-Cook 本构方程[8],材料的Mises 等效应力

式中:ε为等效塑性应变,﹒ε*=﹒ε/﹒ε0,为量纲一等效塑性应变,通常取﹒ε0=1 s-1;T*=T/T m,为量纲一温度,T m 为材料熔化温度;A、B、C、n 和m 为材料常数。

弹丸材料的状态方程一般可以通过Grüneisen 状态方程表示

式中:K 1、K 2 和K 3 是材料常数,μ=ρ/ρ0-1。

根据Johnson-Cook 模型,单元内材料的损伤演化定义为

式中:D ≤1,当D=1 时材料出现断裂;Δε为每个积分步长中等效塑性应变的增量,εf为等效断裂应变,该应变与应变率、温度及压力有关

4 滑移面算法

通常滑移面处理时,滑移面条件(不可侵入条件、法向接触力为压力的条件及切向摩擦力的条件)都是不等式约束,需要通过拉格朗日乘子法或罚参数法将约束条件引入变分原理。具体处理时,事先通过试探的方法将不等式约束改为等式约束引入方程并进行方程求解,计算后利用不等式条件对解的结果进行检查,这样反复迭代求解,直至满足条件。传统的滑移面处理过程异常复杂,计算量非常大,精度低,对于RKPM 不适合,也不易实现,为此引入了新的滑移面处理法[9]。

图1 给出了二维问题滑移面处理方法。计算过程中,首先检查从滑移面质点A 在下一时间步时处于主滑移面的哪两个质点之间,并检查从滑移面质点与主滑移面的关系。在图2 中,ED和F N 是两条铅垂线,E F 与DF、与F H 相互垂直。

从滑移面质点A 移动后可能位于主滑移面的3 种位置,即A 可能在Ⅰ、Ⅱ或Ⅲ区域内点′ 或。如果在Ⅰ内位置,说明从质点移动后在滑移面之外,不发生接触,此时质点间的各种物理量不需要发生交换;如果在Ⅱ、Ⅲ内,说明从质点已侵入主滑移面内,按不侵入条件,则立即将质点A垂直拉到滑移面上。对于区域Ⅱ而言,)应垂直拉到EF 面的处,并且重新确定从质点A(A′)和主质点E、F 的法向速度分量,并保持切向速度分量不变

图1 滑移面Fig.1 Sliding surface

图2 二维滑移面分析Fig.2 Analysis on 2D sliding surface

式中:m、l 分别表示质点质量和线段长度,v+、v-分别是拉动后、前的质点的法向速度。可以证明,式(18)～(20)既保证了滑移面质点的动量守恒,也保证了动量矩守恒。如果质点A 位移后在区域Ⅲ内A′3位置,此时还须检查相邻主滑移面质点线段F G 的位置情况。如果FG 在F H 以下,如FG 在FG2 位置,则将垂直拉到FG2上的,速度可按上述原理发生交换;如果FG 在F H 以上(如F G在FG1位置)或与F H 在同一位置,质点A 只与质点F 在垂向发生速度调整,即A 和F 调整后的垂向速度

切向方向速度不发生改变。

对于三维问题计算,首先以质点为顶点,将主滑移面划分成相连的三角形区域。当主、从滑移面发生接触时,判断从滑移面上的某一质点处于主滑移面哪一个三角形区域,并按照二维问题类似的标准判断,将该从接触质点与三角形顶点的3 个质点发生速度调整。

5 速度配点法

对于给定的速度边界条件及滑移面质点速度的调整,本文中提出了速度配点法。在计算域Ω内的N P个质点可以分为NB个给定速度质点和NN个非给定速度质点,。RKPM 法中,质点速度可类似给出

给定质点速度

将式(23)代入式(22),可以得到

简记为

则实现方法为

本质边界条件及加速度边界条件也可采用与本配点法相同的方法实现。

6 弹丸侵彻靶板算例

根据理论编制了非线性侵彻动力过程数值分析的RKPM 法,对具体算例进行计算,验证上述算法。

钢质圆柱型弹丸直径15 mm,长度50 mm,周边固定的钢板厚度300 mm,弹丸以1 km/s 速度、与钢板表面法线方向倾斜角θ=10°入射钢板。状态方程参数分别为;本构和损伤参数分别为:A=790 M Pa,B=510 M Pa,C=0.015,n=0.27,m=1.05,=0.47×10-5,D2=-9.0,D3=3.0,D4=0,D5=0.78,溶化温度Tm=1.8 kK,初始温度T0=293 K。

由于具有对称性,取弹丸和靶板的一半作为分析模型,并划分成质点,作为计算基点。图3 给出弹丸侵彻前后的质点动态变化图,可以看出靶板在弹丸侵彻过程中破碎和弯曲变形。

图3 弹丸与靶板侵彻过程的质点动态图Fig.3 The dynamic particles' photo of pellet into steel plate

图4 ～5 为计算给出的弹丸和靶板不同时刻应力分布云图。可以看出,t=30 μs 时弹丸端头部出现明显的破坏断裂,局部应力达到1.2 GPa 以上。靶板应力成近乎圆形的分布规律向四周扩散。

图4 弹丸应力分布Fig.4 Nephogram of stress of pellet

图5 靶板应力分布Fig.5 Nephogram of stress of steel plate

图6 弹丸尾部中心点水平速度和加速度Fig.6 Velocity and acceleration at the center point of end part of pellet

图6 分别给出了弹丸尾部中心点水平方向的速度和加速度变化图,弹丸经过侵彻后水平方向的速度基本减至0。图7 计算给出了弹丸端头部中心点的速度变化和加速度变化图。

选择初始速度分别为0.6、1.0、1.2 和1.5 km/s 进行计算,并将弹体入射后的剩余速度与理论解[10]进行比较,结果如图8,两者完全一致。

图7 弹丸端头部中心点垂直速度和加速度Fig.7 Vertical velocity and acceleration at the center point of head part of pellet

图8 弹体剩余速度Fig.8 Residual velocity of pellet

7 结 论

利用再生核质点法计算模拟了侵彻过程。在侵彻阶段的计算中,将弹丸和靶板划分核质点,不需要划分成单元,利用RKPM 方法对侵彻过程进行大应变、高应变率计算,避免了有限单元法网格重分问题,减少了计算的难度,提高了计算的精度和速度。此外,为了方便解决RKPM 方法中的滑移面接触和边界条件问题,提出了新的滑移面算法和配点法。数值模拟结果表明,本方法易于编程,计算精度较高。

[1] 宋顺成,才鸿年.弹丸侵彻混凝土的SPH 算法[J].爆炸与冲击,2003,23(1):56-60.SONG Shun-cheng, CAI H ong-nian.SPH algorithm for projectile penetrating into concrete[J].Explosion and Shock Waves,2003,23(1):56-60.

[2] 宋顺成,李国斌,才鸿年,等.战斗部对混凝土先侵彻后爆轰的数值模拟[J].兵工学报,2006,27(2):230-234.SONG Shun-cheng,LI Guo-bin,CAI Hong-nian,et al.Numerical simulation of penetration-then-detonation of concrete target with projectile[J].Acta Armamentarii,2006,27(2):230-234.

[3] Lucy L B.A numerical approach to the testing of the fission hypothesis[J].The Astronomical Journal, 1977,82:1013-1024.

[4] H ayhurst C J,Clegg R A.Cylindrically symmetric SPH simulations of hypervelocity impacts on thin plates[J].International Journal of Impact Engineering, 1997,20:337-348.

[5] Dancygier A N, Yankelevsky D Z.High strength concrete response to hard projectile impact[J].International Journal of Impact Engineering,1996,18(6):583-599.

[6] Liu W K,Jun S,Zhang Y F.Reproducing kernel particle method[J].International Journal for Numerial Methods in Fluids, 1995,20(8-9):1081-1106.

[7] Liu W K, Jun S,Li S F,et al.Reproducing kernel particle methods for structural dynamics[J].International Journal for Numerial Methods in Engineering,1995,38(10):1655-1679.

[8] Holmquist T J, Johnson G R, Cook W H.A computational constitutive model for concrete subjective to large strain,high strain rates, and high pressure[C]∥14th International Symposium on Ballistics.USA:American Defense Preparedness Association, 1993:591-600.

[9] 赵光明,宋顺成.高速冲击过程数值分析的再生核质点法[J].力学学报,2007,39(1):63-69.ZH AO Guang-ming, SONG Shun-cheng.The reproducing kernel particle method for numerical analysis of highspeed impact process[J].Chinese Journal of Theoretical and Applied Mechanics, 2007,39(1):63-69.

[10] 马晓青.冲击动力学[M].北京:北京理工大学出版社,1991.