基于区间不确定性参数辨识的某软后坐火炮发射过程建模

2023-05-23郭书奇侯保林

兵工学报 2023年4期

郭书奇, 侯保林

(南京理工大学机械工程学院, 江苏南京 210094)

0 引言

火炮软后坐发射是一种利用火炮前冲击发原理减小后坐力的发射方式,通过火炮前冲产生的动能消耗部分后坐能量,合理应用此项技术可使火炮后坐力比同规格常规发射火炮减小50%左右,在火炮威力不变的情况下达到火炮轻量化和小型化的目的,提高火炮综合性能。与常规火炮相比,软后坐火炮发射需要先通过机械装置后拉至一挂卡位置,此时前冲机储存部分能量;火炮架体上有卡锁装置,可对后坐部分进行固定,完成前冲准备;当释放卡锁时,前冲机赋予火炮后坐部分前冲的能量,在达到某一前冲位置或速度时,火炮击发;后坐部分需要先抵消前冲运动积累的动能,再进行后坐。当火炮后坐至最大后坐长时复进,直至复进至卡锁位置挂卡完成一次发射循环[1-2]。国内外对软后坐火炮的研究一直在继续,但由于该技术的特殊性,国内外公开的研究成果不多,近几年可获取的相关论文也较少。早在20世纪60年代,美国岩岛兵工厂(Rock Island)就研制了世界上第1门软后坐火炮XM204,验证了软后坐发射原理在减轻后坐力方面的有效性,并且在测试中积累了较全面的设计经验与试验数据,奠定了火炮软后坐技术的基础。到21世纪初,同样位于岩岛的曼都斯公司应用软后坐技术成功研发105 mm“鹰眼”(Hawkeye)与155 mm“布鲁图斯”(Brutus)软后坐火炮,并申请了相关专利[3],其中“鹰眼”于2017年在阿布扎比防务展亮相,其部分装备了美国国民警卫队,成为第1种制式化列装的中大口径软后坐火炮。

除美国外,韩国于2007年前后针对软后坐发射技术进行了研究,Shin等[4]对软后坐火炮发射的故障发射现象进行了研究,通过MATLAB软件数值仿真的方式获得了软后坐火炮故障发射的动力学特性,提出并验证了应用磁流变缓冲器控制其冲击的有效性;Yang等[5]通过仿真与试验对比,建立基于软后坐火炮反后坐装置参数的运动微分方程,研究了反后坐装置不同参数对软后坐发射特性的影响;Kang等[6]采用数值仿真方法,对比了软后坐火炮与常规火炮的发射特性和后坐力数值,进一步验证了软后坐发射降低后坐力的可行性与有效性。

国内在软后坐火炮的研究上也进行了一定的探索。张志坚等[7]对软后坐发射原理和反后坐装置进行了系统性的分析,建立了比较完善的软后坐火炮理论体系;杜中华等[8]基于反后坐装置的反面问题模型,得到了常规火炮直接应用软后坐原理射击的动态特性曲线,对比获得一般火炮不能直接进行软后坐发射的结论。

随着各类高响应传感器应用于工程试验中,多种测试数据成为研究动态特性与规律的重要切入点,通过智能算法的帮助,可以更加细致地揭示系统特性。在各类研究中,李丹丹等[9]应用麻雀搜索算法对J-A磁滞模型参数进行了辨识;徐岩等[10]应用遗传粒子群算法对光伏逆变器参数进行了辨识;刘柏希等[11]应用区间辨识方法对LuGre摩擦模型参数进行了辨识;Li等[12]应用区间理论对制退机相关参数进行了优化,取得了较好的效果。

本文以某软后坐火炮为研究对象,详细地对制退机与前冲机的运动特性进行分析,建立其发射过程的运动解析模型;基于此模型和3发发射试验获得的数据,对该系统内的多个重要参数进行了区间参数辨识,并采用辨识参数和试验数据对模型进行了验证,其结果显示辨识取得了较好的结果,揭示了软后坐火炮的发射特性。

1 软后坐火炮发射过程建模

本文的研究对象为某型软后坐坦克炮,采用一个制退机和两个前冲机的反后坐装置方案,并配有卡锁和后拉装置。其结构简图及受力示意图如图1所示。图1中,θ为射角,Fpt为炮膛合力,2Ff为身管轴线方向上受前冲机力,f为等效摩擦力,Fψ为制退机力,G为后坐部分所受重力,FN为后坐部分垂直于身管轴线的重力分量,FG为后坐部分沿身管轴线的重力分量。

图1 软后坐火炮结构简图Fig.1 Structure diagram of a soft recoil gun

由图1可知,在发射过程中,软后坐火炮后坐部分在身管轴线方向上受前冲机力2Ff、制退机力Fψ、后坐部分运动所受的等效摩擦力f、后坐部分沿身管轴线的重力分量FG、炮膛合力Fpt等力的作用。对于火炮的发射过程中,以后坐部分为研究对象,假设所有力的方向均作用于射击平面内,且忽略各处间隙、力矩与弹性的影响,建立简化的一维运动微分方程。以炮口方向为正方向,挂卡位置为零点建立坐标系。因此软后坐发射过程总运动微分方程可写为

(1)

式中:FG=mhgsinθ;f=μmhgcosθ,μ为等效摩擦系数,mh为后坐部分质量,g为重力加速度。制退机力Fψ和等效摩擦力f的方向由速度方向决定,使用sign函数处理其方向。等效摩擦力包含摇架处的摩擦力、前冲机与制退机内密封件摩擦力等各处摩擦力。

1.1 前冲(复进)过程前冲机力和制退机力建模

1.1.1 前冲(复进)过程前冲机力表达式

前冲机提供前冲动力,也作为复进机复位火炮状态。本文系统前冲机采用了液体气压式前冲机,结构简图如图2所示。图2中,pf为前冲机内气体的瞬时压力,可视为前冲机游动活塞处压强,Af为游动活塞工作面积,ps为活塞杆处压强,As为前冲机活塞杆工作面积,wf为运动过程中外筒液体流速,wb为运动过程中内筒液体流速,v为火炮后坐部分运动速度。

图2 前冲机结构简图Fig.2 Structure diagram of preshooting mechanism

该前冲机由高压气体腔内的气体作为储能介质向外推动或吸收活塞杆的能量,在气体腔内高压气体的变化可描述为气体多变过程,即

(2)

式中:pf0为前冲机内气体的初压力;V、V0分别为气体某瞬时体积和初体积;n为气体多变指数。设X0为游动活塞位移,X为活塞杆位移,数值上等于火炮后坐部分的位移,则任意位置的气腔体积为

V=V0-AfX0

(3)

将式(3)代入式(2),可得

(4)

采用伯努利方程计算反后坐装置特性,常引入以下假设:1)液体不可压缩,内腔始终充满液体;2)流动是一维定常流动;3)以地球为惯性参考系;4)液体在某一点处的流速和压力视为液体在该点所在截面上的平均流速和平均压力。在此基础上,制退机内液体的流动满足连续性方程(质量守恒)与伯努利方程(能量守恒)。此假设适用于后文各处伯努利方程的推导。

根据如图2所示液体流动方向,前冲过程前冲机的伯努利方程为

(5)

(6)

将式(4)代入式(6),则前冲过程(v≥0 m/s)前冲机力Ff为

(7)

1.1.2 前冲(复进)过程制退机力处理

本文研究的软后坐制退机为带有筒壁沟槽和中心液流孔的多通道制退机,其结构简图如图3所示。图3中e为沟槽宽度,γ为中心流口开启角度。

图3 制退机结构简图Fig.3 Structure diagram of recoil mechanism

该制退机设计与一般制退机不同,在活塞上设有小型液流孔A和大型液流孔B,并有滑环在活塞上随活塞杆运动,该滑环功能类似于单向阀。在前冲和复进过程中,滑环受液体作用打开大型液流孔B,使液流阻力大幅下降;而在后坐过程中,滑环在液体的作用下封闭此大型液流孔B,液体沿筒壁沟槽和中心小型液流孔A流入制退机后腔。

前冲和后坐过程制退机内部液体流动如图4所示。图4中,pzf、pzb为制退机的工作腔压力与非工作腔压力,Axa、Axb分别为筒壁沟槽面积和活塞中心流口面积,ufa、ufb分别为筒壁沟槽和中心流口入口处的液体流速,uba、ubb分别为筒壁沟槽和中心流口出口处的液体流速。从图4中可以看出,该制退机工作状况下液体流动情况复杂,直接建立数学模型困难。为求解制退机内的液流特性,本文将结合Fluent软件仿真对制退机进行研究。

图4 制退机液体流动示意图Fig.4 Diagram of liquid flows in recoil mechanism

对于前冲过程,以制退机内液体为网格划分对象,设置活塞非工作腔端面处的绝对压强为0 MPa,应用标准k-epsilon湍流模型求解,层铺方式更新动网格,其网格划分如图5所示。

图5 Fluent软件动网格划分Fig.5 Fluent dynamic meshing

本文流体仿真采用Profile方式赋予动网格的运动速度,速度数据来源于试验测得的前冲速度。仿真结果如图6所示。在3 m/s的前冲速度下,制退机阻力仅为2 kN,整个过程中制退机液压阻力Fψ相对其他各力小很多,对全过程影响有限,为简化计算,采用拟合方程描述该力,处理为关于活塞杆速度v的二次方程:

Fψ=127.09v2+404.07v

(8)

图6 制退机前冲阻力Fig.6 Set-forward resistance of recoil mechanism

1.2 后坐过程前冲机力和制退机力建模

1.2.1 后坐过程前冲机力表达式

前冲机处于后坐过程时,内部液流方向与图2所示相反,根据伯努利方程写出后坐过程中压强关系为

(9)

对式(9)采用与1.1.1节所述参数值与简化原则,此处不再赘述,则后坐过程(v<0 m/s)前冲机力Ff为

(10)

1.2.2 后坐过程制退机力表达式

制退机活塞上小型液流孔A为定面积液流孔,6条筒壁沟槽沿圆周方向均匀分布,沟槽深度h=h(x)为一分段函数,在制退机轴线方向上连续变化,各段斜率不同。理论上,沟槽深度的变化将影响制退机力的大小,合适的沟槽深度可获得“平缓”的后坐阻力。

对于图4(b)所示的液体流动过程,可写出其连续性方程为

AOv=QO=qa+qb=ubaAxa+ubbAxb

(11)

式中:AO为制退机活塞工作面积;QO为液体的总体积流量;qa、qb分别为筒壁沟槽和中心流口的体积流量。针对制退机的工作腔压力pzf与非工作腔压力pzb的关系,并考虑总流与分流的体积流量时[13],伯努利方程可写为

(12)

(13)

式中:za、zb、zc、zd分别为制退液各段的水头高度;Hr1和Hr2为液体由工作腔向非工作腔流动过程中的比能损失,其与比动能呈正比,即

(14)

ξa、ξb为对应流动过程的能量损失系数。在1.1节的假设条件下可得ufa=ufb=0 m/s,pzb=0 MPa;忽略液流中心的高度变化,za≈zb,zc≈zd。将式(12)、式(13)相加并消去相等项,再引入筒壁沟槽和中心流口的液压阻力系数Ka、Kb,Ka=1+ξa,Kb=1+ξb,化简得

(15)

将式(11)代入式(15),得

(16)

则制退机后坐过程的液压阻力Fψ为

(17)

另外,由于制退机筒壁沟槽在X=0.07 m处消失,液体通道由两类通道变为一类通道,即X<0.07 m时,Axa=0。在此点前后,制退机内流场发生变化,设制退机中心流道阻力系数在两段分别为Kb和Kc。在X>0.07 m时,制退机中心流道的阻力系数为Kb,在X<0.07 m时为Kc。同时,制退机力表达式也发生变化,由伯努利方程可得第2段制退机力表达式,推导与上文所述过程相似。则制退机全工作段力的表达式如式(18)所示

(18)

1.2.3 后坐过程流量比计算

根据伯努利方程的推导可以看出,计算并获取两组液体通道的流量是求解制退机力的必要途径。对此,采用Fluent软件对制退机的工作过程进行仿真,监测其活塞受力与液体通道流量的变化规律。与1.1.2节中网格的划分和设置相同,对制退机匀速后坐和变速后坐运动进行仿真,其流量比与流口面积比如图7所示,其中变速运动的速度数据来源于发射试验获取的后坐速度。

图7 液体通道流量比仿真结果Fig.7 Simulation results of liquid flow ratio

由图7可以看出,匀速运动和变速运动时液体通道内液体的流量比关系基本相同。结合图8后坐过程中的制退机力仿真结果,验证了网格划分与边界条件设置的正确性。

图8 制退机力流体仿真结果Fig.8 Simulation results of recoil mechanism force

根据图7的仿真结果,设制退机6条筒壁沟槽的体积流量与中心流口的体积流量的比值为qα,该比值与液体通道面积比有一定的关系,因此可建立式(19)所示的qα表达式:

(19)

式中:α为针对流口流量比的修正系数,其数值在后坐过程中不是定值,作为影响制退机阻力的重要参数,将α作为待辨识参数由参数辨识过程求出贴近实际情况的最适值。在求解出流量比规律后,式(18)内的uba、ubb与后坐速度v的关系也可解出,进一步可求出制退机力Fψ。

1.2.4 液压阻力系数参数化

对于制退机内的筒壁沟槽,其截面积沿活塞运动方向变化。为更加精确地描述制退机力,提高式(18)的描述精度,将筒壁沟槽液压阻力系数Ka定义为连续变化量。为满足该连续性要求,对4段连续且坡度不同的筒壁沟槽对应的阻力系数进行参数化。结合流体力学中对液压阻力系数的描述方式,设Ka的变化规律为

(20)

式中:C为待求解的阻力系数;γa、γb、γc、bx、by为使方程连续的待求解常数。

当X>0.25 m时,筒壁沟槽截面积为定值,因此Ka=C为一定值。在此方程组的限制下,Ka的取值范围为[1,C]的子区间。为满足连续性,将节点处的截面积数值代入,通过数学处理,可以消去式(20)中的bx、by,即液压阻力系数Ka可表示为关于C、γx、γy的函数关系式。

1.3 炮膛合力

炮膛合力是火药燃气作用于火炮身管上,沿身管轴线方向向后的力,使火炮后坐部分产生后坐运动。参照现有的内弹道理论,一般以式(21)计算炮膛合力Fpt:

(21)

式中:tf为炮膛合力出现时刻;tk为膛压消失时刻;φ为次要功计算系数;ω为火药质量;Ab为炮膛横断截面积;pt为内弹道计算的平均压力;tg为弹丸出膛时刻;χ为炮口制退器冲量特征量;Fg为弹丸出膛瞬间膛内平均压力;b为反应炮膛合力衰减快慢的时间常数。经内弹道计算软件计算,可获得炮膛合力的基本曲线如图9所示,其中释放软后坐火炮卡锁的瞬时记为横轴零点。

图9 炮膛合力曲线Fig.9 Bore resultant force curve

需要特别说明的是,考虑到装药的特性,在不同工况下,炮膛合力的产生存在一定的延迟,这一延迟由多种因素共同影响,一般在数十毫秒到数百毫秒之间,特殊情况下该延迟可能更大。这一特性对于常规发射火炮的后坐过程影响不大,但是,由于软后坐火炮击发过程发生于前冲过程中,数十毫秒的延迟会使得炮膛合力出现于更远的前冲位置,导致后坐起始位置与预期不同,影响软后坐发射的一致性。对于定位置击发的软后坐火炮,设xe为设定的前冲击发距离,xf为对应tf时刻的后坐部分位移。tf和xf是关系到动态发射的起始点的重要参数,对预测火炮的软后坐发射过程有着重要意义。

2 软后坐发射模型不确定参数辨识

参数辨识是通过数学方式将模型与试验数据结合起来,用于验证模型或预测试验的重要方法。由建模过程可以看出,本软后坐火炮的发射模型包含由多个参数组成、复杂的力表达式,每一个参数对模型都存在影响。对这种存在多参数的复杂模型,部分凭经验给出的参数会导致建模精度低,试验中不得不进行多次试验,各类成本和复杂性成倍提升。

本文将基于区间相似度理论,对建立的软后坐发射模型进行区间参数辨识,以位移与速度等时间序列状态量为目标数据,采用相似度量化仿真值与试验值的相似性,再通过参数区间优化,建立与实际系统更为相似的模型。一方面可以验证建模的准确性,另一方面进一步的明确各参数所在区间,了解软后坐火炮特性,更加准确地预测发射。参数辨识与建模结合对研究复杂系统特性有着重要意义。

2.1 辨识目标函数的建立

辨识问题一般可以转化求解特定参数使目标函数值最大或最小的优化问题。辨识的理论与计算过程基于文献[14-15]。相似度S是量化两组数据间相似度大小的重要指标。对于描述软后坐火炮后坐部分运动状态的位移与速度状态量,单一考虑位移或速度会使得某些特征或运动趋势表现的不明显,易忽视系统的特性造成的影响。为此,可以综合考虑二者对运动的影响,定义综合相似度Sg,兼顾位移与速度在数值上的相似程度。Sg可以表示为

Sg=ξSv+(1-ξ)Sx

(22)

式中:Sv和Sx为后坐部分速度与位移的时间序列状态量的相似度;ξ为权重系数,可以人为提升或降低某一参量在辨识中的重要性,可根据需要调整。

在本文的参数辨识中,参数数量较多,根据所建模型,除后坐部分速度v与位移x外,解析模型中各个力包含的待辨识参数有:前冲机力(n,Kf),制退机力(C,Kb,Kc,γa,γb,α),摩擦力μ,炮膛合力xf。这些参数可以分为两类,第1类参数仅有1个,为n,这类参数是可通过试验测量或根据经验确定参数取值的,组成集合U;第2类参数有9个,包括C、Kb、Kc、γa、γb、α、Kf、xf、μ,这类参数在试验中无法直接测量,且不确定性强,难以确定为某一值,此类参数组成区间集合VI。因此辨识的目标函数为

fI(U,VI)=ξSv(U,VI)+(1-ξ)Sx(U,VI)

(23)

因此,软后坐系统参数辨识问题即转化为寻求目标函数fI(U,VI)最大值的优化问题。

2.2 区间参数辨识问题描述

根据参数辨识方法对参数的分类,目标函数的最大值寻优过程可分为两步。第1步为对第1类参数的辨识,采用区间可能度的方法计算第1类参数的最优值。设定参数的初始参数取值空间为ΩU与ΩV,在VI∈ΩV初始区间的同时获取最优的参数确定值的集合U,使得辨识结果与测试结果的时间序列的区间可能度最大,即辨识得到集合U的各参数值。这一步在明确第1类参数数值的同时,可以验证建模与初始区间选择的准确性,对复杂系统有重要作用。第2步辨识则通过代入第1步得到的集合U值,优化区间VI的上下界,即使空间ΩV在满足约束条件的前提下尽可能小,尽可能提高辨识精度。

在此过程中,保持第2类参数集合VI不变,应用区间分析理论对第1类参数的计算初值Ua进行迭代优化,不断获得新的解集Ub,使速度与位移的时间序列状态量的相似度大于设定值。此过程中,Ub不断替换Ua并再计算,直到P(fI(Ub,VI))满足要求,即Pmax=P(fI(U,VI))≥λ。根据需要,一般λ的取值范围为[0.5,1],该值为U的参数选取系数,在本文中λ取0.9。

该过程采用了遗传算法作为基础寻优算法,设定种群大小为40,进化代数20代,交叉概率0.9,变异概率0.1,解出确定的集合U。

辨识的第2步是在第1步获得的集合U的基础上,对第2类参数区间VI进行辨识。同第1类参数的辨识过程类似,第2类参数辨识的最终目的是在区间不确定性变量的初始给定区间VI范围内找出更优区间。在这一部分中,区间参数辨识所获得的结果也为区间。区间的优劣辨别问题也通过区间可能度方法进行比较[16-17],其基于模糊集来构造区间可能度,是一种定量比较区间之间的优劣程度的方法。数学模型可表述为

(24)

式中:ε为速度与位移的时间序列状态量的相似度的设定目标值,ε=[ε1,ε2]T。基于上述两步辨识过程,可获得所求参数的辨识结果。

3 试验与辨识结果

本文试验以某软后坐火炮为试验对象,进行固定台架射击,射角θ为2°,后坐部分质量mh为 2 150 kg,卡锁完全解锁时上下前冲机压力pf0均为4.7 MPa。使用常温正装药模拟穿甲弹进行2发软后坐发射试验和1发常规发射试验。其中,2发软后坐发射的前冲击发位置不同,第1发击发点在前冲243 mm处,第2发击发点在前冲265 mm处,即所述xe不同。除xe外,3发射击中各参数均保持相同。

试验数据采用压力与激光位移传感器获取。制退机筒壁设压力传感器,监测制退机工作腔压力值;激光位移传感器监测炮尾的位移,二者共同通过高响应采集仪记录数据。图10为软后坐火炮射击试验场景。

图10 软后坐火炮射击试验Fig.10 Soft recoil gun test

3.1 试验结果

试验采集数据如图11所示,其中制退机工作腔压力测试值已通过计算转换为制退机力。

图11 发射试验测试数据Fig.11 Test data of firing

本文选取第1发软后坐试验采集的0.4 s内的数据进行参数辨识,其他两发数据进行验证。通过有限差分法获取第1发射击的后坐速度v,经多重滤波后,保留重要特征,并修正不良数据带来的影响[18],获得图12所示平滑速度曲线。

图12 第1发软后坐射击速度Fig.12 Firing speed curve of the first round

为避免综合相似度区间上界值低,各参数初始区间选取较大范围,根据经验设定参数辨识的初始搜索区间取值如表1所示。

表1 参数辨识初始区间

3.2 辨识结果

以第1发射击试验数据为辨识目标数据,对第1类参数进行辨识可得如表2所示的辨识结果。

表2 第1发射击第1类参数辨识结果

气体多变指数n的辨识结果为1.394,与经验值相近。综合相似度区间Sg的区间下界为0.582 6,区间上界可达到0.936 3,初始区间有较大的优化空间。表明本组辨识所得的第1类参数值可信,从侧面验证了模型可以正确描述实际发射过程,可进行第2类参数辨识。进一步地,取相似度的辨识目标值ε=[ε1,ε2]T=[0.95,0.95]T对第2类参数进行辨识,获得如表3的第2类参数辨识结果。

表3 第1发射击第2类参数辨识结果

根据辨识结果,将最适值代入式(20),计算得出制退机的液压阻力系数Ka值如图13所示。

图13 沟槽液压阻力系数KaFig.13 Hydraulic resistance coefficient Ka of grooves

进一步地,经辨识过程可获得后坐部分位移X、速度v、后坐过程制退机力Fψ等数据的仿真结果,与测试数据的对比情况如图14所示。从图14中可以看出,两组数据的相似性较好,辨识结果与综合相似度值匹配。

图14 第1发辨识最适值仿真与试验数据对比Fig.14 Comparison between test data and optimal identification value simulation results of the first round

3.3 辨识参数验证

为验证辨识结果的准确性,采用第2和第3发射击数据对辨识所得的参数与模型进行验证,在此仅对试验获得的后坐位移与后坐过程制退机力两个测试值进行对比验证。

对于第2发软后坐发射试验,改变前冲击发距离xe为265 mm,考虑与第1发射击时相同的炮膛合力延迟时间,第2发发射的xf值为302 mm。在不改变其他参数值的条件下,采用此数值进行软后坐仿真,获得了图15(a)、图15(b)的对比结果。

对于第2发发射过程,保持反后坐装置的各参数不变,改软后坐为常规后坐,以火炮坐标系为基础,即在0.35 m处静止发射,获得的仿真结果与试验数据对比如图15(c)、图15(d)所示。仿真与试验数据存在部分差别,但基本趋势相同,可以认为辨识参数与模型能够反映该常规发射过程。

图15 第2和第3发辨识最适值仿真与试验数据对比Fig.15 Comparison between test data and optimal identification value simulation results of the second and third round

综合图15内各子图信息可以看出,第1发发射参数的辨识最适值代入第2和第3发进行仿真与同发射条件下的试验值相似性较好,证明了本文中软后坐系统的解析建模和参数辨识所得参数可以较好的反应软后坐发射的实际工作过程,参数辨识建模在复杂系统的特性研究中具备重要的使用价值。但仍要注意的是,在多不确定性的影响下,参数辨识建模与实际系统仍存在一定差异。在本文研究的软后坐系统中,这些差异在发射全过程的后半段逐渐累积,在制退机部分,系统的不确定性带来的差异较为明显。