超音速悬索火箭橇下拉幅度运动仿真分析

游培寒, 胡 瑜, 缪 昕, 祝逢春, 李 洋

(1 95856部队, 南京 210000; 2 解放军理工大学, 南京 210007)

0 引言

悬索火箭橇是一种约束弹道试验装置(图1)。它利用两根平行布设的悬空悬索作为滑轨,火箭橇悬置于两根悬索中间,发动机推动火箭在悬索约束方向高速运动。高速悬索火箭橇可用于测试飞行状态下弹上设备的工作状态。国外资料表明火箭橇在悬索上的运动速度不能超过悬索的特征波速,两者接近时会使波动在火箭橇头部积累,造成悬索断裂。在前期试验中发现,火箭橇在高速运动中会对悬索形成向下“拉扯”,文中重点就这种向下拉扯进行分析。

悬索火箭橇的运动是火箭橇飞行运动和悬索动态应力传播的结合,牵涉到空气动力学和绳索弹性力学两个领域。文中采用这两个系统分开建模的方式,将火箭橇运动对悬索计算单元位置和速度的扰动引入悬索的动态计算;同时将悬索动态应力引入火箭橇运动模型和俯仰角计算中,使两个力学模型在同一个系统中仿真计算。文中第一和第二部分首先介绍悬索和火箭橇的受力模型;第三部分给出悬索动态仿真与火箭橇运动关系;第四部分给出了仿真方法和参数;第五部分给出仿真结果。第六部分计算了不同初始应力对悬索下拉作用的影响。

1 悬索动力学模型

悬索动力学建模实质是绳索的动力学建模,与拦阻索模型有相似之处[1],都可以采用基于牛顿定律的绳索动力学模型[2-3]进行建模。这里只考虑悬索在x、y两个方向的运动,忽略绳子的横切向应力,同时假定绳子在拉伸时横截面积一定,则其运动方程为:

(1)

(2)

式中：s0、sd分别表示静态和动态条件下悬索纵向各质点沿悬索伸展方向坐标;x、y为悬索上各质点惯性空间横坐标和纵坐标;xd、yd表示质点的动态坐标;Td为悬索动态条件下的拉力;m为悬索单位长度质量;X、Y分别为作用在悬索单位长度的横向外力和纵向外力。根据虎克定律[1]:

(3)

式中:εd表示悬索单位长度的动态伸长量;δd为动态应力;δ0为静态应力;E为悬索纵向弹性模量。悬索拉力与应力存在如下关系:

T0=δ0A

(4)

A为悬索横截面积,通过式(3)、式(4)可以得到悬索动态拉力为:

(5)

设αd为悬索纵轴切线在动态坐标sd处与x轴夹角,则:

(6)

(7)

将式(5)～式(7)式代入式(1)、式(2)可以得到

(8)

(9)

差分计算时,用计算单元来替换上式中的各微元,设悬索计算单元的位置为pi=(xi,yi),初始应力条件下计算单元重心距离为L0,i=∂s0,文中仿真时设L0,i=1 m,动态计算中的计算单元重心距离Li≈∂sd,αd为计算单元重心之间的夹角为:

(10)

2 火箭橇动力学模型

火箭橇建模不但要考虑火箭橇飞行过程中所受的空气阻力,同时要考虑悬索对火箭橇运动的扰动。一方面,由于在空间布设的悬索不是直线,其方向与火箭橇轴线存在夹角,造成运动中的火箭橇对悬索的横向拉扯。由于发动机不断消耗燃料,火箭橇质量中心前移,转动惯量减小,尾部支撑座向上力臂增加等因素会增加火箭橇俯仰角的变化速度。另一方面,在拉扯过程中悬索的总体长度增加,悬索的内部应力增加,其对火箭橇存在一定的拉升作用。除此以外还要考虑悬索与火箭橇之间还存在摩擦。采用火箭橇外形结构如图2所示,前后支架的距离为Lr=1 m。火箭橇在滑行过程中受到推力P、摩擦力M、空气阻力Z、悬索前后支承力T1和T2以及自身重力的影响。各力的方向如图3所示,火箭橇后端悬索与横坐标夹角为α1,火箭橇前端悬索与横坐标夹角为α2,火箭橇的俯仰角为α3。

则x、y两个方向的外力应为:

2T1cos(α1-α2)]sinα3+(P-M-Z)cosα3

(11)

2T1cos(α1-α2)]cosα3-(P-M-Z)sinα3

(12)

式中:xr、yr是火箭橇当前坐标;mr是火箭橇当前质量,火箭橇在滑行过程中攻角很小,火箭橇空气阻力主要来自零攻角阻力,其为[3]:

(13)

式中:Cr为阻力系数,取值与速度有关如图4所示;Sr=3 848.45 mm为参考面积;vr为火箭橇速度;ρ为空气密度。

火箭橇推力变化与火箭橇发动机型号有关,不是文中重点,这里不做详细介绍。根据前期应力摩擦试验结论,火箭橇运动时平均摩擦力M=10 N。随着发动机的工作火箭橇的质量会下降,这里假设质量下降与工作时间存在线性关系:

mr=m0-Kmt/Tr

(14)

式中：m0=10.35 kg为火箭橇的初始质量,t为已经工作时间;Tr=0.65 s为发动机总工作时间;Km=3.6 kg为发动机装药消耗质量。以上是火箭橇的质点运动模型,除此以外还需要考虑火箭橇在飞行过程中俯仰角变化,旋转主要是由于火箭橇前后的拉力角度变化和火箭橇质心的前移造成的。质心与火箭橇后支架的距离为:

Lx=Lr(0.5+0.3t/Tr)

(15)

同时,火箭橇转动惯量为:

Jr=J0,z-KJt/Tr

(16)

式中:KJ=0.362,则火箭橇俯仰角α3为:

T2(Lr-Lx)sin(α2-α1)

(17)

3 火箭橇与悬索运动关系处理

火箭橇在运动过程中悬索与火箭橇之间实际存在的应力关系很难直接得到。在仿真中采用位置速度修正方法,也就是用火箭橇的瞬时速度和位置修正悬索对应计算单元,然后进行悬索动态计算,并将计算结果转换为前后支承力T1、T2引入火箭橇姿态运算,如图5所示。图中pr1=(xr1,yr1)、pr2=(xr2,yr2)为火箭橇前支架和后支架的位置,pk为火箭橇上悬索计算单元的重心位置,Lk为pk到pk+1计算单元重心距离,vr为火箭橇速度。设pk=(xk,yk),vpk=(vxk,vyk)为该点速度,则有:

yk=yr2-d1sinα3

(18)

式中：d1为pk到pr1的距离,由于火箭橇在悬索上滑动,所以xk的运动相对较小,同时vxk≈0。

(19)

vyk=vrsinα3

(20)

悬索拉力由计算单元重心距离的变化得到,从图4可以看出:

Lk=Lr-d1+Dpr1,pk+1

(21)

Lk-1=d1+Dpr2,pk-1

(22)

式中:Dpr1,pk+1为pr1到pk+1的距离;Dpr2,pk-1为pr2到pk-1的距离。

4 仿真计算方法

文中采用龙格库塔法进行仿真计算[4]。以上各节介绍的计算公式可以集成为一个差分运算环节:

(23)

式中:χ=(pr1,pr2,vr,p1…pK,vp1…vpK),t为时间;φ(χ,t)是一个非线性函数,其计算步骤如下:

第一步,根据火箭橇受力状态,按照式(11)～式(14)计算火箭橇位置pr1、pr2、速度vr;

第二步,根据悬索当前应力得到火箭橇前后支承力T1、T2,按照发动机工作时间得到火箭橇质量、质心、转动惯量等参数,按照式(17)计算火箭橇俯仰角α3;

第三步,根据式(18)～式(20)计算火箭橇上悬索计算单元的位置和速度,并根据式(21)～式(22)计算对应计算单元距离;

第四步,根据式(8)、式(9)计算悬索上各计算单元的位置pk和速度vpk。

这里设仿真时间步长Δt=10-6s,龙格库塔法计算公式如下:

(24)

5 仿真计算结果

悬索长度为270 m,初始点比悬索末端高出5.32 m,参照国内某高分子聚乙烯绳的产品参数,悬索单位长度质量m=0.019 kg/m,悬索横截面积A=5.027×10-5m2,悬索纵向弹性模量E=3.4 GPa。发射时火箭橇头部距离悬索初始点10 m,初始状态下,火箭橇后端悬索与横坐标夹角α1=1.163 2°,火箭橇前端悬索与横坐标夹角为α2=1.094 4°,火箭橇的俯仰角为α3=1.128 8°。火箭橇推力变化参照了某型70 mm火箭弹发动机参数。火箭橇其他参数参照以上各节内容。悬索初始应力T0=2 400×9.8 N,初始状态下悬索前后支承力T1=T2≈T0,火箭橇瞬时飞行状态如图6所示。图中圆点代表火箭橇头部位置。

火箭橇飞行0.65 s悬索内各计算单元所受应力如图7所示,图中可以看出在火箭橇所在位置出现了应力突变,这主要是由于悬索与火箭橇之间存在摩擦产生的。

6 悬索下拉幅度分析

火箭橇飞行过程中对悬索的下拉与火箭橇初始俯仰角,悬索初始应力、悬索弹性模量等因素有关,当火箭橇参数和初始位置固定时,初始俯仰角变化很小。

悬索所用绳索材料的弹性模量也相对稳定,所以在实际试验中造成悬索下拉的主要考虑因素是悬索初始应力。我们对初始应力为T0=1 800×9.8 N、2 000×9.8 N、2 200×9.8 N、2 400×9.8 N、2 600×9.8 N五种状态下8对悬索的下拉幅度(图8)进行了计算。可见悬索初始应力越小,则火箭橇对悬索的下拉幅度越大。

7 结论

文中利用Newton法对悬索进行动态建模,参考某型火箭弹发动机参数对火箭橇进行飞行力学建模,采用龙格库塔法仿真计算高速火箭橇飞行过程中系统的相关状态参数,重点考察了火箭橇对悬索的下拉作用,并得到不同初始应力条件下的下拉幅度曲线,为后续火箭橇系统的进一步改进设计和试验方法研究提供了有利的理论分析依据。

参考文献:

[1] 张新禹. 阻拦索的动力学特性及仿真研究 [D]. 哈尔滨: 哈尔滨工程大学, 2011: 12-21.

[2] 金栋平, 胡海岩. 行进绳索在横向流体激励下的运动 [J]. 力学学报, 2001, 33(4): 525-529.

[3] 周新年. 工程悬索与柔性吊桥理论、设计、案例 [M]. 北京: 人民交通出版社, 2008: 18-35.

[4] 钱杏芳, 林瑞雄, 越亚男. 导弹飞行力学 [M]. 北京: 北京理工大学出版社, 2000: 36-48.