1、1合理假設為簡化計算，假設：1)拉傘過程中，引導傘與回收物運動軌跡為一條直線，且物傘系統作平面運動。2)不考慮風的影響，物傘系統沒有升力。3)不考慮傘繩的彈性。4)引導傘、物體和拉直中的傘系統微元質量作為三個質點處理。2動力學建模對于回收物及已拉出的主傘系統在回收物的坐標上的運動方程 (1)其中，分別為回收物和主傘系統已拉出部分的氣動阻力；為傘繩內的張力，最大值即為拉直力。對于拉直中的主傘系統微元質量可寫成如下的動量方程 (2)其中，為傘系統微元質量拉出速度。 (3)為拉出的微元質量運動速度，在不考慮傘繩彈性情況下，。當然在方程(2)中忽略了主傘系統微元的自身重力。為傘系統拉出長度。為引導傘速

2、度。因此，(2)式可改寫為 (4)其中，為拉動的傘系統單位長度質量。對于引導傘及傘袋在回收物坐標上的運動方程 (5)其中，為未拉出的主傘系統質量，分別為引導傘及傘袋的阻力。為了消除(1)式中，將(4)式代入(1)式可得 (6)由于不考慮傘繩彈性，上式可改寫為 (7)又知 (8)拉直階段的軌跡計算可以根據式(3)，(5)，(7)和(8)六個一階微分方程組，已時間為計算步長，采用龍格-庫塔法算出各個時間內的個參數。對于先拉傘繩法 (9)其中，為傘繩及傘衣總質量，為傘衣底邊沿傘繩方向的單位長度的質量，為傘衣底邊寬度，為傘繩長度，為傘系統全長，為傘衣名義直徑。3計算求解 MATLAB中內置求解微分方程

3、的ode函數，現選取ode45進行求解計算。計算輸入參數如下：rho=1.293; %空氣密度g=9.8; %重力加速度CAys=1.5; %引導傘阻力系數CAd=0.2; %傘衣套阻力系數CAw=1; %物體阻力系數CAe=0.6; %已拉出物體阻力系數Qys=0.5*rho*x(4)2*CAys; %引導傘氣動阻力Qd=0.5*rho*x(4)2*CAd; %傘衣套氣動阻力Qw=0.5*rho*x(5)2*CAw; %物體的氣動阻力Qe=0.5*rho*x(5)2*CAe; %已拉出傘系統的氣動阻力Mw=60; %物體質量Mys=5; %引導傘質量(包括傘衣、套傘包)Msh1=0.6; %

4、傘繩的質量密度Msy1=5; %傘邊的質量密度b=0.05; %傘邊的寬度Msh=6; %傘繩總質量Msy=10; %傘衣總質量Lsh=10; %傘繩全長Lxt=14.05; %傘系統全長Do=8; %傘衣名義直徑Fsh=50; %拉出阻力4結果圖1 Xd-時間歷程圖圖2 Yd-時間歷程圖圖3角度-時間歷程圖圖4 Vys-時間歷程圖圖5 Vw-時間歷程圖圖6 L-時間歷程圖圖7 FL-時間歷程圖圖8 軌跡曲線圖5結果優化在前面的討論計算中，假設傘繩為非彈性體，無伸長。顯然，傘繩假設為非彈性體與實際情況出入較大。然而，由于在拉直過程中，除了傘繩拉直瞬間拉直力較大外，在拉動傘繩時傘繩內張力較小，因

5、此，傘繩假設為非彈性體對計算拉直階段軌跡來說影響很小，但對于拉直力的計算影響較大，因而在拉直力的計算中必須要考慮傘繩的彈性。假設傘繩為彈性的無質量的彈簧，則彈性縱波沿傘繩的傳播速度 (10)其中，為傘繩的理論彈性模量，為單根傘繩的單位長度質量。現假設傳播速度為常數，則傘繩張力傳播是一個典型的一維波動方程 (11)其中，為傘繩任意橫剖面的位移，為傘繩剖面的坐標。為了求解式(11)，必須確定傘繩拉直瞬間的初始條件和邊界條件。初始條件可以寫成 (12)邊界條件 (13)假設傘繩的彈性變形符合虎克定律，則 (14)由初始條件和邊界條件可以解出式(11)，得 (15)將式(15)對微分，并代入式(14)

6、，可得 (16)將式(4)與式(16)聯立求解，可以求得拉出速度 (17) 其中，將式(17)代入式(16)，可以求得傘繩拉直力 (18)6拉直力計算 選取傘繩材料為凱夫拉29，彈性模量為62GPa，傘繩數量為40，其余參數參考不考慮傘繩彈性時的計算輸入輸出。利用MATLAB可以簡便快捷的求解計算結果，則。求解程序見附錄。附錄A主函數clearclct,x=ode45('aero',0:0.001:1,0;0;15;100;100;0);%輸入初值%參數rho=1.293; %空氣密度g=9.8; %重力加速度CAys=1.5; %引導傘阻力系數CAd=0.2; %傘衣套阻力系

7、數CAw=1; %物體阻力系數CAe=0.6; %已拉出物體阻力系數Qys=0.5*rho*x(4)2*CAys; %引導傘氣動阻力Qd=0.5*rho*x(4)2*CAd; %傘衣套氣動阻力Qw=0.5*rho*x(5)2*CAw; %物體的氣動阻力Qe=0.5*rho*x(5)2*CAe; %已拉出傘系統的氣動阻力Mw=60; %物體質量Mys=5; %引導傘質量(包括傘衣、套傘包)Msh1=0.6; %傘繩的質量密度Msy1=5; %傘邊的質量密度b=0.05; %傘邊的寬度Msh=6; %傘繩總質量Msy=10; %傘衣總質量Lsh=10; %傘繩全長Lxt=14.05; %傘系統全長

8、Do=8; %傘衣名義直徑Fsh=50; %拉出阻力m=length(t);Me=zeros(m,1);Dml=zeros(m,1);FL=zeros(m,1);if x(6)<=Lsh Me=Msh1*x(6);elseif x(6)<=Lxt Me=Msh1*Lsh+Msy1*b+(Msy-Msy1*b)*(1-4*(Lxt-x(6)/Do)2);endMv=Msh+Msy+Msy1*b-Me; %未拉出主傘的質量for i=1:m;if x(i,6)<=Lsh Dml(i)=Msh1;elseif x(i,6)<=Lsh+b Dml(i)=Msy1;elseif

9、x(i,6)<=Lxt Dml(i)=-1.25*Msh1*x(i,6)+17.5*Msh1; endendk=menu('拉直階段各變量時間歷程曲線：','Xd-時間歷程圖', 'Yd-時間歷程圖','角度-時間歷程圖','Vys-時間歷程圖','Vw-時間歷程圖','L-時間歷程圖','FL-時間歷程圖','軌跡曲線圖');switch k case 1 figure(1); plot(t,x(:,1); xlabel('t/s

10、9;); ylabel('Xd/m'); title ('Xd-時間圖'); case 2 figure(2); plot(t,x(:,2); xlabel('t/s'); ylabel('Yd/m'); title('Yd-時間圖'); case 3 figure(3); plot(t,x(:,3); xlabel('t/s'); ylabel('角度'); title('角度-時間圖') case 4 figure(4); plot(t,x(:,4); xlabe

11、l('t/s'); ylabel('Vys/m'); title('Vys-時間圖') case 5 figure(5) plot(t,x(:,5); xlabel('t/s'); ylabel('m/s'); title('Vw-時間圖') case 6 figure(6) plot(t,x(:,6); xlabel('t/s'); ylabel('L/m'); title('L-時間圖') case 7 for i=1:m; FL(i)=Dml(i

12、)*(x(i,5)-x(i,4)2+Fsh; end figure(7); plot(t,FL); xlabel('t/s'); ylabel('FL/N'); title ('拉直力-時間圖'); otherwise figure(8); plot(x(:,1),-x(:,2); xlabel('Xd/m'); ylabel('Yd/m'); title ('軌跡曲線圖'); end附錄B求解函數function X=aero(t,x)%參數rho=1.293; %空氣密度g=9.8; %重力加速

13、度CAys=1.5; %引導傘阻力系數CAd=0.2; %傘衣套阻力系數CAw=1; %物體阻力系數CAe=0.6; %已拉出物體阻力系數Qys=0.5*rho*x(4)2*CAys; %引導傘氣動阻力Qd=0.5*rho*x(4)2*CAd; %傘衣套氣動阻力Qw=0.5*rho*x(5)2*CAw; %物體的氣動阻力Qe=0.5*rho*x(5)2*CAe; %已拉出傘系統的氣動阻力Mw=60; %物體質量Mys=5; %引導傘質量(包括傘衣、套傘包)Msh1=0.6; %傘繩的質量密度Msy1=5; %傘邊的質量密度b=0.05; %傘邊的寬度Msh=6; %傘繩總質量Msy=10; %

14、傘衣總質量Lsh=10; %傘繩全長Lxt=14.05; %傘系統全長Do=8; %傘衣名義直徑Fsh=50; %拉出阻力m=length(t);Me=zeros(m,1);if x(6)<=Lsh Me=Msh1*x(6);elseif x(6)<=Lxt Me=Msh1*Lsh+Msy1*b+(Msy-Msy1*b)*(1-4*(Lxt-x(6)/Do)2);endMv=Msh+Msy-Me; %未拉出主傘的質量X=zeros(6,1);X(1)=-x(5)*cos(x(3); %x(1)表示XdX(2)=x(5)*sin(x(3); %x(2)表示YdX(3)=-g*cos(x(3)/x(5); %x(3)表示角度X(4)=-g*sin(x(3)-(Qys+Qd-Fsh)/(Mys+Mv); %x(4)表示VysX(5)=-g*sin(x(3)-(Qw+Qe+Msh1*(x(5)-x(4)2+Fsh)/(Mw+Me); %x(5)表示VwX(6)=x(5)-x(4); %x(6)表示傘系統