马 佳,陈光宋,吉 磊,钱林方,徐亚栋
(1.长沙理工大学 土木工程学院,湖南 长沙 410114;2.南京理工大学 机械工程学院,江苏 南京 210094)
火炮发射过程中,弹炮间隙的存在使得身管内膛与弹丸前定心部间不可避免地发生接触碰撞,从而对身管产生激励作用。同时,身管和火炮其他部件的运动反过来又会影响弹丸的膛内运动,造成身管振动与弹丸膛内运动的相互耦合。这种相互耦合作用是火炮系统初始扰动的主要因素之一,影响了弹丸出炮口时的运动状态,进而降低了火炮系统的射击精度。因而,研究身管内膛与弹丸前定心部间的接触碰撞响应特性,对于分析弹丸膛内运动规律,揭示弹炮运动耦合机理,进而预测弹丸炮口扰动,最终提高火炮系统的射击精度具有非常重要的意义。
不同于常规间隙转动副中轴与轴承光滑表面间的接触碰撞问题,膛线的几何复杂性使得身管内膛与弹丸前定心部间的接触碰撞问题研究变得更加棘手;此外,火炮发射过程内膛环境恶劣,空间狭窄,现阶段尚不具备开展内膛接触碰撞实测的条件。针对身管内膛与弹丸前定心部间接触碰撞问题的研究,早期的理论建模较少考虑膛线影响。其主要原因在于线膛火炮膛线条数众多,内膛截面呈现类似齿形的不规则形状,导致难以直接建立身管内膛与弹丸前定心部间的接触碰撞力模型[1]。自上世纪70年代以来,众多研究者投身于接触碰撞力模型的研究之中,并先后提出了多种经典的接触碰撞力模型以试图更准确、更高效地模拟2个物体之间的接触碰撞过程,其中Lankarani-Nikravesh(L-N)[2]接触碰撞力模型已被多次用于描述火炮系统中的接触碰撞过程[3-4]。然而,需要注意的是:L-N模型作为“Hertz类模型”的代表,是以准静态球体间的接触碰撞为前提,2个接触碰撞体需满足“非协调接触”条件[5-6]。所谓“非协调接触”,指的是2个物体最初是在一个点或者一条线上发生接触碰撞,即使在载荷的作用下,接触面的尺寸相比于物体本身尺寸而言也是很小的,这种接触形式仅能在“大间隙小力”的情况下得到满足[5-6]。然而,身管内膛与弹丸前定心部几何形状皆为圆柱体,且弹炮间隙较小,因此不完全满足L-N模型前提条件;此外,现存接触碰撞力模型中接触碰撞动力学参数(接触碰撞刚度系数、阻尼系数)的解析计算公式大多针对表面规则几何体,并不适用于复杂表面间的接触碰撞过程。为此,有必要探索一种能够更加准确地模拟身管内膛与弹丸前定心部间接触碰撞过程的模型。在此之前,需充分了解并掌握身管内膛与弹丸前定心部间的接触碰撞响应特性。
本文以某口径火炮为研究对象,建立了身管内膛与弹丸前定心部有限元弹塑性动力学模型。考虑到弹丸膛内复杂的受力状态,为避免其中任一因素的不当处理而影响最终结果,本文仅考虑身管内膛与弹丸前定心部间的接触碰撞过程:从模型收敛性和罚因子设置两方面对所建立的接触碰撞仿真模型进行验证,以保证后续所得数值模拟结果的有效性;基于已验证的有限元模型,分析了不同初始接触碰撞速度及不同内膛磨损状态下相关接触碰撞物理量的变化规律。
目前描述接触碰撞的力学模型主要有2种:拉格朗日乘子法和罚函数法。拉格朗日乘子法以相互接触的物体互不嵌入为约束条件,优点在于符合实际情况,无需设置任何参数,缺点为约束方程随时间变化,方程维度也随之变化,从而引起求解困难。罚函数法是采用罚因子引入接触附加约束条件的数值方法,相比于拉格朗日乘子法,罚函数法不增加系统自由度数,且易于编程实现,因此应用非常广泛。下面列出2种方法的具体推导过程,需要注意的是:下文公式中出现多个下角标,无特殊说明的情况下,小写的下标表示分量,如i,j;大写的下角标表示节点,如I,J,默认对重复的指标求和。针对小写指标,指的是对空间的维数求和;而对于大写指标,则表示对节点的编号进行求和。
令拉格朗日乘子试函数为λ(ξ,t),其变分表示为δλ(ξ,t)。t为时间,ξ为接触碰撞物体表面坐标,δ为数学变分符号。不考虑切向摩擦情况下的拉格朗日乘子法弱形式可表示为[7]
δpL=δp+δGNL≥0
(1)
式中:δpL和δp分别为总的虚功率和变形虚功率,δGNL为侵彻率变分,下标L表示拉格朗日乘子法,N表示法向。计算公式如下:
(2)
(3)
式中:γN为接触碰撞速度的法向分量。
采用Lagrangian网格对接触体进行有限元离散化,以材料坐标的形式表示速度场的有限元近似,则物体A和B的速度场为[7]
(4)
(5)
若物体A和B的节点编号不同,则2个速度场可统一表示为
vi(X,t)=NI(X)viI(t)
(6)
在接触表面上,拉格朗日乘子场λ(ξ,t)是由一个C-1场近似的,因而有:
(7)
式中:λ(ξ,t)≥0,ΛI(ξ)是C-1形函数,Sc为接触界面。
相应的变分函数分别为
δvi(X)=NI(X)δviI, δλ(ξ)=ΛI(ξ)δλi
(8)
式中:δλ(ξ)≤0。于是可以推出[7]:
(9)
以节点速度形式表示的相互侵彻率为
γN=ΦiI(ξ)viI(t)
(10)
(11)
式中:nA和nB分别表示物体A和B的法向。
则接触弱形式为
(12)
(13)
结合式(1)、式(9)和式(12)可得以矩阵形式表示的运动方程[7]:
δvT(Fint-Fext+Ma)+δ(vTGTλ)≥0
(14)
考虑到δv和δλ的任意性,得到运动方程和相互侵彻条件为[7]
(15)
式中:O为零矩阵。
利用拉格朗日乘子法引入接触界面约束条件虽然可使其得到精确满足,然而因增加了方程的自由度数,且系数矩阵包含零对角元素,从而导致求解不便。此外,因涉及未知接触力的联立方程求解,需采用迭代方法计算,增加了计算费用,且与显式算法不相容,所以只限于隐式算法中,常用于静态和低速的接触问题。
罚函数法是采用罚因子引入接触附加约束条件的数值方法,相比于拉格朗日乘子法,罚函数法不增加系统自由度数,且易于编程实现,因此应用非常广泛。对于本文的身管内膛与弹丸前定心部接触碰撞问题,采用的即是罚函数法。
罚函数的弱形式为[7]
δpP=δp+δGNP=0
(16)
式中:下标P表示罚函数法。
(17)
式中:βP为罚因子,HP(γN)为Heaviside函数。
(18)
(19)
将式(10)代入式(19)可得:
(20)
式中:Fc为法向接触碰撞力,可表示为
(21)
当忽略摩擦影响时,结合式(9)、式(16)和式(20),有:
δpP=δvT(Ma+Fint-Fext)+δvTFc
(22)
考虑到非约束虚节点速度的任意性,可得罚函数法的有限元离散方程[7]:
Ma+Fint-Fext+Fc=O
(23)
以某口径火炮为原型,分别建立身管和弹丸的三维八节点六面体有限元数值仿真计算模型。有限元网格需准确刻画弹丸前定心部、身管膛线等关键部位的几何细节,以确保精准搜索小间隙下的接触碰撞现象。
在进行接触碰撞分析时,作如下假设:
①弹带处非本文研究重点,网格划分较疏,而且不考虑与身管内膛的相互作用;
②考虑到计算效率以及实用性问题,仅建立包裹弹丸前定心部的一段身管;
③鉴于接触碰撞问题和高速摩擦问题皆属于非常复杂的力学问题,参考现存文献中的普遍做法[4,8],本文暂不考虑切向摩擦力和法向接触碰撞力间的耦合效应,也即不考虑弹丸的旋转以及高速轴向运动对法向接触碰撞过程的影响。
仿真模型的边界条件:身管外表面全约束,身管左右侧面自由。初始条件设置是相对于图 1所示的弹尾坐标系。为模拟弹丸与身管间的斜碰撞现象,给弹丸施加沿yD轴正向的平动速度以及绕zD轴的转角速度。
图1 弹尾坐标系示意图
此外,为模拟实际结构中由于制造工艺、磨损等原因造成的弹炮间隙,将身管内膛沿半径方向缩短某一长度。为保证计算精度,添加质量点对弹丸进行配重,确保有限元模型的质量、质心位置和转动惯量与标准情况一致,建立的离散模型如图2所示。
图2 离散模型
设置弹丸前定心部和身管内膛接触碰撞区域,由于弹丸前定心部网格相较于身管内膛网格稀疏,因而选择弹丸前定心部接触表面为主面,身管内膛接触表面为从面,并采用平衡主从接触搜索算法进行接触节点的搜索。
JOHNSON和COOK[9]于1983年提出了一种用于金属大变形、高应变率和高温情况下的本构模型,简称J-C模型。得益于其简便的表达形式和明确的物理意义,一经提出便在工程领域得到了极为广泛的应用[10]:
(24)
(25)
本文中身管材料塑性本构采用J-C模型描述[11],弹丸材料塑性本构采用双线性硬化模型描述,身管与弹丸材料的弹性部分皆采用与应变率无关的胡克定律描述。
弹丸和身管材料具体参数设置如下:弹丸密度为7.75 g/cm3,弹性模量为206 GPa,泊松比为0.297 2。塑性设置:当εp=0时,σs=1 086 MPa;当εp=0.26时,σs=1 560 MPa。身管密度为7.83 g/cm3,弹性模量为206 GPa,泊松比为0.3,塑性设置如表1所示。
表1 身管材料塑性参数设置
为了确保所获有限元模拟结果具有较高精度,本小节分别从模型收敛性和罚因子设置两方面对所建立的身管内膛与弹丸前定心部接触碰撞模型进行验证。
为了检验仿真模型的收敛性,首先对比分析接触碰撞区域不同网格尺寸下相关动力学响应变化;然后基于选定的网格尺寸,检验仿真过程稳定步长设置的合理性。
鉴于罚函数法计算所得的穿透量值受罚因子控制,因此,利用罚函数法仿真所得结果的精度也很大程度上依赖罚因子的选择。太小的罚因子会导致过大的穿透量,从而造成计算错误,难以满足计算精度要求;而过大的罚因子,虽然能在一定程度上提高计算结果精度,但是会引起接触区域刚度增加,显式计算步长减小,严重影响计算效率,且会造成计算过程的不稳定。协调计算精度和计算成本之间的矛盾,一直以来都是有限元仿真过程的难题[12]。鉴于本文的身管内膛与弹丸前定心部接触碰撞过程采用罚函数法描述。因而需对罚因子的设置进行验证。
3.1.1 收敛性分析
分别建立了接触碰撞区域最小网格尺寸lmin为0.6 mm,0.4 mm,0.2 mm和0.1 mm的有限元仿真模型,对比分析了上述不同网格尺寸下相关接触碰撞响应的变化趋势,如图3所示。图中,δpe为接触碰撞穿透量。观察图3可知,当最小网格尺寸为0.6 mm时,接触碰撞持续时间以及接触碰撞响应峰值与其他几种情况相比有较为明显的不同。随着网格尺寸减小,计算结果逐渐趋于稳定,然而计算成本不断上升。本文后续仿真计算皆基于0.1 mm的网格模型。
根据文献[10]可知,显示计算的稳定步长应满足Δt 图3 身管内膛与弹丸前定心部接触碰撞区域不同网格尺寸影响 3.1.2 罚因子设置验证 本小节仿真对比了不同罚因子取值情况下的计算结果[12]:0.01,0.1,0.5,1,5。相应的仿真结果如图4所示,观察结果易知:当罚因子为0.01时,仿真计算结果失真;随着罚因子增大,计算结果逐渐趋于稳定。当罚因子等于5时,可以认为所得稳定结果即为Abaqus的计算结果。为平衡有限元仿真计算成本与计算精度之间的矛盾,对于后续身管内膛与弹丸前定心部的接触碰撞过程,统一将罚因子设置为1。 图4 身管内膛与弹丸前定心部接触碰撞过程不同罚因子设置对比 文献[4]研究表明,身管内膛与弹丸前定心部间的最大接触碰撞速度大约处于3~5 m/s范围,考虑到弹丸偏心、身管磨损以及强装药等情况,针对本文的身管内膛与弹丸前定心部接触碰撞问题,初始接触碰撞速度γN0研究范围设置为1~7 m/s。每隔0.1 m/s进行一次仿真分析,其中转角速度皆设置为1 rad/s,内膛半径磨损为0.1 mm。 图5为初始平动速度5 m/s,转角速度1 rad/s情况下,最大接触碰撞力时刻身管内膛表面接触压力和Mises应力云图。对应的弹丸表面接触压力和Mises应力云图结果如图6所示。图5和图6分别给出了身管内膛和弹丸表面沿yD轴正向最大接触碰撞时刻的接触压力和Von Mises应力云图。通过统计发现,最大接触碰撞时刻,与弹丸前定心部发生接触碰撞的身管膛线条数大约为16根(总条数为36根),参照文献[5-6]相关结论可知:身管内膛与弹丸前定心部不满足非协调接触条件。此外,在一次碰撞过程中(包括接触碰撞压缩阶段和恢复阶段),皆为单向碰撞,未发生双边碰撞的情况。 图7列举了几种初始接触碰撞速度下的身管内膛与弹丸前定心部接触碰撞数值模拟结果。观察易知:初始接触碰撞速度越大,相应的接触碰撞力、接触碰撞穿透量越大;此外,初始接触碰撞速度越大,碰撞结束时刻速度与初始时刻速度的差异越大,说明接触碰撞过程能量损失越大。 图8给出了不同初始碰撞速度下接触碰撞力与接触碰撞穿透量,以及接触碰撞速度间的变化关系。 图5 最大接触碰撞时刻身管内膛表面数值模拟结果 图6 最大接触碰撞时刻弹丸表面数值模拟结果 图7 身管内膛与弹丸前定心部不同初始接触碰撞速度计算结果对比Ⅰ 图8 身管内膛与弹丸前定心部不同初始接触碰撞速度计算结果对比Ⅱ 观察结果易得:随着初始碰撞速度的增加,接触碰撞力与接触碰撞穿透量所围面积逐渐增大,表明接触碰撞过程能量损失逐渐增加。此外接触碰撞力抖动的原因,参考文献[5,13],推断是由膛线几何形状非连续性以及应力波传递等综合作用导致。 参照靶场实验数据和文献[14-15]可知:身管内膛磨损主要集中于膛线起始端及炮口,其余位置磨损量较小。然而,为了检验内膛磨损对接触碰撞过程相关物理参量的影响,针对本文所研究的火炮结构,分别建立了磨损量为0.1 mm,0.3 mm,0.4 mm和0.5 mm(以上数据皆为半径差)的有限元网格模型,并对相应计算结果进行分析比较,具体如图9所示,图中,Wdepth为内膛磨损量。 观察图9易知:不同磨损状态下,各仿真模型计算结果差异较小。这主要归因于本文仿真模型的初始条件和边界条件施加方式,磨损量的不同仅会造成弹炮初始接触碰撞姿态的少许差异。此外,通过统计不同磨损状态下身管内膛与弹丸表面接触碰撞膛线条数可知:最大接触碰撞时刻,与弹丸前定心部发生接触碰撞的身管膛线条数皆为15或16根,接近身管总膛线条数的一半。 图9 不同磨损程度下身管内膛与弹丸前定心部接触碰撞过程仿真结果对比 本文针对某口径火炮,建立了身管内膛与弹丸前定心部接触碰撞有限元弹塑性动力学模型,模拟实际接触碰撞过程。通过研究身管内膛与弹丸前定心部间的接触碰撞响应特性,从而为后续发展一种精度更高的弹炮接触碰撞力模型奠定基础。对比分析了不同初始接触碰撞速度及不同内膛磨损程度下相关接触碰撞物理量变化规律,观察仿真结果可知: ①接触碰撞最大时刻,与弹丸前定心部发生接触碰撞的身管膛线条数皆接近身管总膛线条数的一半。 ②随着初始碰撞速度增大,相应的接触碰撞力和接触碰撞穿透量逐渐增大。此外,初始碰撞速度越大,接触碰撞力与接触碰撞穿透量图像所包围的滞回环面积也越大,表明整个接触碰撞过程能量损失越大。3.2 不同初始接触碰撞速度下的仿真结果
3.3 考虑不同内膛磨损程度的仿真结果
4 结论