降落伞拉直阶段轨迹及拉直力计算

1、1合理假设为简化计算，假设：1）拉伞过程中，引导伞与回收物运动轨迹为一条直线，且物伞系统作平面运动2）不考虑风的影响，物伞系统没有升力。3）不考虑伞绳的弹性。4）引导伞、物体和拉直中的伞系统微元质量 dm作为三个质点处理。2动力学建模对于回收物及已拉出的主伞系统在回收物的Xq坐标上的运动方程dVwdtgsinDw De Dlmwme(1)其中，Dw， De分别为回收物和主伞系统已拉出部分的气动阻力; 绳内的张力，最大值即为拉直力。对于拉直中的主伞系统微元质量dm可写成如下的动量方程Fl为伞udm Fl Fsh dt其中，U为伞系统微元质量拉出速度。dLu v Vys dtv为拉出的微元质量运动

2、速度，在不考虑伞绳弹性情况下，方程（2）中忽略了主伞系统微元的自身重力。L为伞系统拉出长度速度。因此，（2）式可改写为当然在Vys为引导伞Fldm dL2Fshm uFshdt dt其中，m貉为拉动的伞系统单位长度质量。 对于引导伞及伞袋在回收物Xq坐标上的运动方程dVysg sindtDysDd - Fshmysmv其中，mv为未拉出的主伞系统质量，Dys，Dd分别为引导伞及伞袋的阻力。 为了消除式中Fl，将式代入式可得ys?dVwdtgsin2Dw De muFshmwme由于不考虑伞绳弹性，上式可改写为dvwdtgsinDw De m VwVysmwme又知ddtdXddtdyddtgc

3、os ；VwVw cosVwSi n(8)拉直阶段的轨迹计算可以根据式(3),(7)和(8)六个一阶微分方程组，已 时间t为计算步长，采用龙格-库塔法算出各个时间t内的个参数Vw， Vys， Xd， yd， L 。对于先拉伞绳法mshLm emsh Lsh mshmsh - mshb4Lxt-L14D。mv msh msy0 L LshLsh L Lxt(9)其中，msh msy为伞绳及伞衣总质量，msy为伞衣底边沿伞绳方向的单位长 度的质量，b为伞衣底边宽度，Lsh为伞绳长度，Lxt为伞系统全长，Do为伞衣名 义直径。3计算求解rho=1.293;%空气密度g=9.8;%重力加速度CAys=

4、1.5;%引导伞阻力系数CAd=0.2;%伞衣套阻力系数CAw=1;%物体阻力系数CAe=0.6;%已拉出物体阻力系数MATLAB中内置求解微分方程的ode函数，现选取ode45进行求解计算 计算输入参数如下：Qys=0.5*rho*x(4F2*CAys;%引导伞气动阻力Qd=0.5*rho*x(4)A2*CAd;%伞衣套气动阻力Qw=0.5*rho*x(5F2*CAw;%物体的气动阻力Qe=0.5*rho*x(5F2*CAe;%已拉岀伞系统的气动阻力Mw=60;%物体质量Mys=5;%引导伞质量（包括伞衣、套伞包）Msh1=0.6;%伞绳的质量密度Msy 仁5;%伞边的质量密度b=0.05;

5、%伞边的宽度Msh=6;%伞绳总质量Msy=10;%伞衣总质量Lsh=10;%伞绳全长Lxt=14.05;%伞系统全长Do=8;%伞衣名义直径Fsh=50;%拉出阻力4结果图1 Xd-时间历程图“时闾圈3*5图2 Yd-时间历程图15.1615.1415.1215.0815,Q615.CM1S.O2图3角度-时间历程图图4 Vys-时间历程图图5 Vw-时间历程图L-时间阳2015E2110t/s图6 L-时间历程图图7 FL-时间历程图(13)图8轨迹曲线图5结果优化在前面的讨论计算中，假设伞绳为非弹性体，无伸长。显然，伞绳假设为非 弹性体与实际情况出入较大。然而，由于在拉直过程中，除了伞绳

6、拉直瞬间拉直 力较大外，在拉动伞绳时伞绳内张力较小，因此，伞绳假设为非弹性体对计算拉 直阶段轨迹来说影响很小，但对于拉直力的计算影响较大，因而在拉直力的计算 中必须要考虑伞绳的弹性。假设伞绳为弹性的无质量的弹簧，则弹性纵波沿伞绳的传播速度iE 2c (10)m其中，E为伞绳的理论弹性模量，m为单根伞绳的单位长度质量。 现假设传播速度c为常数，则伞绳张力传播是一个典型的一维波动方程2(11)2C 2x其中，为伞绳任意横剖面的位移，X为伞绳剖面的坐标。为了求解式(11)，必须确定伞绳拉直瞬间的初始条件和边界条件。初始条件 可以写成(x,0)0-(x,0)0边界条件(0,t) (u VR)t(0,t

7、) u Vr假设伞绳的弹性变形符合虎克定律，则FlnT nE x x 0 由初始条件和边界条件可以解出式(11),得Vr U (x,t) (x ct)c将式(15)对x微分，并代入式(14)，可得FlnEc将式(4)与式(16)联立求解，可以求得拉出速度u其中，kB VrFsyu 2B 1 4k 1cnEmsyBmsh将式(17)代入式(16)，可以求得伞绳拉直力Vr1 4k 1Fl nEc 2B(14)(15)(16)(17)(18)6拉直力计算选取伞绳材料为凯夫拉29,弹性模量为62GPa,伞绳数量为40,其余参数 参考不考虑伞绳弹性时的计算输入输出。利用 MATLAB可以简便快捷的求解计

8、 算结果，则Fl 149990。求解程序见附录。附录A主函数clearclct,x=ode45(aero,0:0.001:1,0;0;15;100;100;0);% 输入初值%空气密度%重力加速度%引导伞阻力系数%伞衣套阻力系数%物体阻力系数%已拉出物体阻力系数%引导伞气动阻力%伞衣套气动阻力%物体的气动阻力%已拉出伞系统的气动阻力%物体质量%引导伞质量 （包括伞衣、套伞包 ）%伞绳的质量密度%伞边的质量密度%伞边的宽度%伞绳总质量%伞衣总质量%伞绳全长%伞系统全长%伞衣名义直径%拉出阻力Me=Msh1*Lsh+Msy1*b+(Msy-Msy1*b)*(1-4*(Lxt-x(6)/Do)A2)

9、;%未拉出主伞的质量% 参数rho=1.293;g=9.8;CAys=1.5;CAd=0.2;CAw=1;CAe=0.6;Qys=0.5*rho*x(4F2*CAys;Qd=0.5*rho*x(4F2*CAd;Qw=0.5*rho*x(5)A2*CAw;Qe=0.5*rho*x(5F2*CAe;Mw=60;Mys=5;Msh1=0.6;Msy1=5;b=0.05;Msh=6;Msy=10;Lsh=10;Lxt=14.05;Do=8;Fsh=50;m=length(t);Me=zeros(m,1);Dml=zeros(m,1);FL=zeros(m,1);if x(6)=LshMe=Msh1*x

10、(6);elseif x(6)=LxtendMv=Msh+Msy+Msy1*b-Me; for i=1:m;if x(i,6)=LshDml(i)=Msh1;elseif x(i,6)=Lsh+bDml(i)=Msy1;elseif x(i,6)=LxtDml(i)=-1.25*Msh1*x(i,6)+17.5*Msh1;endend k=menu(拉直阶段各变量时间历程曲线：,Xd-时间历程图,Yd-时间历程图,角度-时间历程图,Vys-时间历 程图:,Vw-时间历程图,L-时间历程图,FL-时间历程图,轨迹曲线图);switch kcase 1figure(1);plot(t,x(:,1)

11、;xlabel(t/s);ylabel(Xd/m);title (Xd- 时间图 );case 2figure(2);plot(t,x(:,2);xlabel(t/s);ylabel(Yd/m);title(Yd- 时间图 );case 3figure(3);plot(t,x(:,3);xlabel(t/s);ylabel(角度);title(角度-时间图)case 4figure(4);plot(t,x(:,4);xlabel(t/s);ylabel(Vys/m);title(Vys-时间图)case 5figure(5)plot(t,x(:,5);xlabel(t/s);ylabel(m/

12、s);title(Vw- 时间图 )case 6figure(6)plot(t,x(:,6);xlabel(t/s);ylabel(L/m);title(L- 时间图 )case 7for i=1:m;FL(i)=Dml(i)*(x(i,5)-x(i,4)F2+Fsh;endfigure(7);plot(t,FL);xlabel(t/s); ylabel(FL/N);title (拉直力-时间图 ); otherwisefigure(8); plot(x(:,1),-x(:,2); xlabel(Xd/m); ylabel(Yd/m); title ( 轨迹曲线图 );end附录 B求解函数f

13、unction X=aero(t,x)% rho=1.293;g=9.8;CAys=1.5;CAd=0.2;CAw=1;CAe=0.6;Qys=0.5*rho*x(4F2*CAys;Qd=0.5*rho*x(4F2*CAd;Qw=0.5*rho*x(5)A2*CAw;Qe=0.5*rho*x(5F2*CAe;Mw=60;Mys=5;Msh1=0.6;Msy1=5;b=0.05;Msh=6;Msy=10;Lsh=10;Lxt=14.05;Do=8;Fsh=50;参数%空气密度%重力加速度%引导伞阻力系数%伞衣套阻力系数%物体阻力系数%已拉出物体阻力系数%引导伞气动阻力%伞衣套气动阻力%物体的气动

14、阻力 %已拉出伞系统的气动阻力%物体质量%引导伞质量 （包括伞衣、套伞包 ）%伞绳的质量密度%伞边的质量密度%伞边的宽度%伞绳总质量%伞衣总质量%伞绳全长%伞系统全长%伞衣名义直径%拉出阻力m=length(t); Me=zeros(m,1); if x(6)=LshMe=Msh1*x(6);elseif x(6)=LxtMe=Msh1*Lsh+Msy1*b+(Msy-Msy1*b)*(1-4*(Lxt-x(6)/DoF2);endMv=Msh+Msy-Me;%未拉出主伞的质量X=zeros(6,1);X(1)=-x(5)*cos(x(3);X(2)=x(5)*sin(x(3);X(3)=-g*cos(x(3)/x(5);X(4)=-g*sin(x(3)-(Qys+Qd-Fsh)/(Mys+Mv);X(5)=-g*sin(x(3)-(Qw+Qe+Msh1*(x(5)-x(4)F2+Fsh)/(Mw+Me); X(6)=x(5)-x(4);附录 C%x(1) 表示 Xd%x(2) 表示 Yd%x(3) 表示角度%x(4) 表示 Vys%x(5) 表示 Vw%x(6)