飞行器建模与仿真

1、飞行器建模与仿真作业班级：姓名：飞设12潘周周学号：2110702015作业一、题目已知微分方程组厂为=X2X2=X3X3=X4X4=X5XX6X6皿a6-XXi-a X2Wo2X3_ 2Woa3X4a4X53_4WoWoa5X6_5Wo用四阶龙哥库塔法求解该方程组。程序1.编写rk4f.m函数fun ctio nr=rk4f(x1,x2,x3,x4,x5,x6,xs,ts,tf,N)%四阶龙格库塔法计算微分方程组程序%x1,x2,x3,x4,x5,x6分别是初值%xs是方程组中的X%ts，tf代表代表左端点和右端点%N是迭代步数w0=50;a1=3.86;a2=7.46;a3=9.13;a4

2、=7.46;a5=3.86;a6=1;aa=zeros(1,N+1);ab=zeros(1,N+1);ac=zeros(1,N+1);ad=zeros(1,N+1);ae=zeros(1,N+1);af=zeros(1,N+1);aa(1)=x1;ab(1)=x2;ac(1)=x3;xd(1)=x4;xe(1)=x5;af(1)=x6;h=(tf-ts)/N;t=ts:h:tf;for i=1:1:N;%ki1(i=1,2,.6)k11=x2;k21=x3;k31=x4;k41=x5;k51=x6;k61=w0A6/a6*(xs-x1-a1*x2/w0-a2*x3/w0A2-a3*x4/w0A

3、3-a4*x5/w0A4-a5*x6/w0A5)；%ki2(i=1,2,.6)k12=x2+h/2*k21;k22=x3+h/2*k31;k32=x4+h/2*k41;k42=x5+h/2*k51;k52=x6+h/2*k61;k62=w0A6/a6*(xs-(x1+h/2*k11)-a1*k12/w0-a2*k22/w0A2-a3*k32/w0A3-a4*k42 /w0A4-a5*k52/w0A5);%ki3(i=1,2,.6)k13=x2+h/2*k22;k23=x3+h/2*k32;k33=x4+h/2*k42;k43=x5+h/2*k52;k53=x6+h/2*k62;k63=w0A6

4、/a6*(xs-(x1+h/2*k12)-a1*k13/w0-a2*k23/w0A2-a3*k33/w0A3-a4*k43 /w0A4-a5*k53/w0A5);%ki4(i=1,2,.6)k14=x2+h*k23;k24=x3+h*k33;k34=x4+h*k43;k44=x5+h*k53;k54=x6+h*k63;k64=w0A6/a6*(xs-(x1+h*k13)-a1*k14/w0-a2*k24/w0A2-a3*k34/w0A3-a4*k44/w 0A4-a5*k54/w0A5);%计算下一步的值x1=x1+h*(k11+2*k12+2*k13+k14)/6;x2=x2+h*(k21+

5、2*k22+2*k23+k24)/6;x3=x3+h*(k31+2*k32+2*k33+k34)/6;x4=x4+h*(k41+2*k42+2*k43+k44)/6;x5=x5+h*(k51+2*k52+2*k53+k54)/6;x6=x6+h*(k61+2*k62+2*k63+k64)/6;%形成向量aa(i+1)=x1;ab(i+1)=x2;ac(i+1)=x3;xd(i+1)=x4;xe(i+1)=x5;af(i+1)=x6;end%画出图型subplot(2,3,1) TOC o 1-5 h z plot(t,aa, -)subplot(2,3,2)plot(t,ab,-)subplo

6、t(2,3,3)plot(t,ac,-)subplot(2,3,4)plot(t,ad,-)subplot(2,3,5)plot(t,ae,r)subplot(2,3,6)plot(t,af,b)2编写主函数调用该函数rk4f(0,0,0,0,0,0,1,0,0.5,100)三、运行结果x1与t的关系0 6010.5、 题目0.50.5莅 o0tLilttn作业2弹道式再入轨迹仿真。已知太空舱质量为 350Kg,从近地点高度为200Km勺轨道返回。偏心率为0.2，轨道倾角为80度，近地点角距为265 度，升焦点赤经100度。已知:r 0 二 6579.89967 km; (0) = -10 ；

7、 (0) =79.8489182889v(0) =7589 .30433867 m / s; (00.54681217 ; A(0) = 99 .955734求轨迹基本参数 r(t), (t) (t), v(t), (t), A(t)。二、用龙格库塔法编写程序求解此题在此函数中要求考虑非球形假设以及国际标准大气。 所以在编写此函数之前应先编写标准大气函数、阻力函数以及重力加速度函数。大气函数function Y=atmosdenty(h,vel,CL)%本函数用来求 den,Kn,ma%输入 h,vel,Cl分别表示高度，速度，特征长度R=287;%空气气体常数g0=9.806;%海平面重力加

8、速度T0=288.15;%海平面温度re=6378140;%地球半径gama=1.405;%海平面比热比b=2/re;p0=1.01325e5;%海平面大气压layers=21;%大气层数Mo=28.964;Na=6.0220978e23; sigma=3.65e-10;M=Mo 28.964 28.964 28.964 28.964 28.964 28.962 28.962 28.880 28.56 28.07 26.92 26.66 .26.5 25.85 24.69 22.66 19.94 17.94 16.84 16.17 16.17; %高度分层所得向量z=0 11.0191 20.

9、0631 32.1619 47.3501 51.4125 71.8020 86 100 110 120 150 160 170 190 .230 300 400 500 600 700 2000*1e3;%温度层向量T=T0216.65216.65228.65270.65270.65214.65186.946 210.65260.65360.65 960.65 1110.61210.651350.651550.651830.652160.652420.652590.6527002700;%热量递减率lr=-6.5e-30 1e-32.8e-3 0-2.8e-3-2e-31.693e-35e-3

10、 1e-22e-21.5e-21e-2 7e-3 5e-3 4e-33.3e-3 2.6e-3 1.7e-3 1.1e-3 0;rou0=p0/(R*T0);% p(1)=p0;T(1)=T0;rou=zeros(22);rou(1)=rou0;%首先计算出每一层端点处的密度值for i=1:layersif (lr(i)=0)rou(i+1)=rou(i)*(T(i+1)/T(i)F(-(1+(1+b* (T(i)/lr(i)-z(i)*gO/(R*(lr( i).*exp(b*g0/(R*lr(i)*(z(i+1)-z(i);elserou(i+1)=rou(i)*exp(-g0*(z(i

11、+1)-z(i)*(1-b*(z(i+1)+z(i)/2)/(R*T(i)endend% 在同温层和变温层采用不同公式计算for i=1:layersif hz(i+1)if (lr(i)=0) TM=T(i)+lr(i)*(h-z(i);den=rou(i)*(TM/T(i)*(1+(1+b* (T(i)/lr(i)-z(i)*gO/(R*lr(i) *exp(b*g0/(R*lr(i)*(h-z(i);else TM=T(i);den=rou(i)*exp(-g0*(h-z(i)*(1-b*(h+z(i)/2)/(R*(T(i);endc=sqrt(gama*R*TM);ma=vel/c;

12、MOL=M(i)+(M(i+1)-M(i)*(h-z(i)/(z(i+1)-z(i);TM=MOL*TM/Mo;Vm=sqrt(8*R*TM/pi); m=MOL*1e-3/Na; n=den/m;F=sqrt(2)*pi* n*sigmaA2*Vm;L=Vm/F;Kn=L/CL;Y=den Kn ma;return ;endend重力加速度函数function g=gravity(r,wd)%gc,gdelt 分别是径向和切向重力加速度(当地水平坐标系下的结果)； %r ， wd 分别表示距离和纬度；ywd=pi/2-wd; mu=3.9860004e14;re=6387.135e3;J2=

13、1.0826e-3;J3=2.532153e-7;J4=1.6109876e-7;gc=-(mu*(1-1.5*J2*(re/rF2*(3*cos(ywd)A2-1)-2*J3*cos(ywd)*(5*cos(y wdF2-3).*(re/r)A3-(5/8)*J4*(35*cos(ywd)A4-30*cos(ywd)A2+3)*(re/r)A4)/rA2);gdelt=-(-3*mu*si n(ywd)*cos(ywd)*(re/r)A2*(J2+0.5*J3*(5*cos(ywd)A2-1)*(re/r)/cos(ywd).+5/6*J4*(7*cos(ywd)A2-1)*(re/r)A2

14、)/rA2);g=gc,gdelt;大气阻力function cd=cdxs(mach,kn)%本函数用来计算阻力；%输入mach表示ma数，kn表示克努増数3 3.8mar=0 0.20.30.40.50.60.80.90.951.051.11.21.62.02.55 10 99;cdr=0.475475 0.475475 0.47576 0.48336 0.488965 0.508345 0.565630.618165 0.668135.1.031795 1.01707 0.990565 0.815955 0.69236 0.60971 0.54606 0.5130.494 0.48317

15、 0.48317; cdc=interp1(mar,cdr,mach);if kn14.5 cdfm=1.75+sqrt(pi)/(2*mach*sqrt(1.41/2); cd=cdfm;elsecdfm=1.75+sqrt(pi)/(2*mach*sqrt(1.41/2); cd=cdc+(cdfm-cdc)*(1/3*log10(kn/(1/2)+0.5113);end轨迹运动函数function gud=guidao(ts,tf,N)%本函数用来求解弹道式运动轨迹% N 迭代次数% ts 迭代初始时间% tf 迭代终止时间format long r=zeros(1,N+1); la=z

16、eros(1,N+1); del=zeros(1,N+1);v=zeros(1,N+1);fai=zeros(1,N+1);A=zeros(1,N+1);x1=6579899.67;x2=-10*pi/180;x3=-79.8489182889*pi/180;x4=7589.30433867;x5=0.54681217*pi/180;x6=99.955734*pi/180;r(1)=6579899.67;la(1)=-10*pi/180;del(1)=-79.8489182889*pi/180;v(1)=7589.30433867;fai(1)=0.54681217*pi/180;A(1)=9

17、9.955734*pi/180;x=zeros(1,N+1);y=zeros(1,N+1);z=zeros(1,N+1);x(1)=r(1)*cos(del(1)*cos(la(1); y(1)=r(1)*cos(del(1)*sin(la(1); z(1)=r(1)*sin(del(1);w=7.292116e-5;m=350;s=4;h=(tf-ts)/N;t=ts:h:tf;for i=1:Nre=6387135;hi=x1-re;vel=x4;CL=0.5;Y=atmosdenty(hi,vel,CL);denty=Y(1);ma=Y(3);kn=Y(2);cd=cdxs(ma,kn)

18、;wd=x3;g=gravity(x1,wd);gc=g(1);gdelt=g(2);%ki4(i=1,2,.6)k11=x4*sin(x5); k21=x4*cos(x5)*sin(x6)/(x1*cos(x3); k31=x4*cos(x5)*cos(x6)/x1; TOC o 1-5 h z k41=(-0.5*de nty*x4A2*s*cd-m*gc*s in (x5)+m*gdelt*cos(x5)*cos(x6)-m*wA2* x1*cos(x3)*(cos(x5).*cos(x6)*sin(x3)-sin(x5)*cos(x3)/m;k5 仁 1/(m*x4)*(m*x4A2/

19、x1*cos(x5)-m*gc*cos(x5)-m*gdelt*s in (x5)*cos(x6)+ m*wA2*x1*cos(x3).*(sin(x5)*cos(x6)*sin(x3)+cos(x5)*sin(x3)+2*m*w*x4*sin(x6)*cos(x3); k61=1/(m*x4*cos(x5)*(m*x4A2/x1*cos(x5)A2*sin(x6)*tan(x3)-m*gdelt*sin(x6)+m*wA2*x1.*sin(x6)*sin(x3)*cos(x3)-2*m*w*x4*(sin(x5)*cos(x6)*cos(x3)-cos(x5)*si n(x3);%ki2(i

20、=1,2,.6)re=6387135;hi=x1+h/2*k11-re;vel=x4+h/2*k41;CL=0.5;Y=atmosdenty(hi,vel,CL);denty=Y(1);ma=Y(3);kn=Y(2);cd=cdxs(ma,kn);wd=x3+h/2*k31;g=gravity(x1+h/2*k11,wd);gc=g(1);gdelt=g(2);k12=(x4+h/2*k41)*sin(x5+h/2*k51);k22=(x4+h/2*k41)*cos(x5+h/2*k51)*sin(x6+h/2*k61)/(x1+h/2*k11)*cos(x3 +h/2*k31);k32=(x

21、4+h/2*k41)*cos(x5+h/2*k51)*cos(x6+h/2*k61)/(x1+h/2*k11);k42=(-0.5sin(x6+h/2*k61)*sin(x3+h/2*k31)*cos(x3+h/2*k31)-2*m*w*(x4+h/2*k41)*( sin(x5+h/2*k51).*cos(x6+h/2*k61)*cos(x3+h/2*k31)-cos(x5+h/2*k51)*sin(x3+h/2*k31);%ki3(i=1,2,.6)re=6387135;hi=x1+h/2*k12-re;vel=x4+h/2*k42;CL=0.5;Y=atmosdenty(hi,vel,C

22、L);denty=Y(1);ma=Y(3);kn=Y(2);cd=cdxs(ma,kn);wd=x3+h/2*k32;g=gravity(x1+h/2*k12,wd);gc=g(1);de nty*(x4+h/2*k41F2*s*cd-m*gc*si n(x5+h/2*k51)+m*gdelt*cos( x5+h/2*k51).*cos(x6+h/2*k61)-m*wA2*(x1+h/2*k11)*cos(x3+h/2*k31)*(cos(x5+h/2*k51)*cos(x6+h/2*k61)*sin(x3+h/2*k31)-sin(x5+h/2*k51)*cos(x3+h/2*k31)/m;

23、k52=1/(m*(x4+h/2*k41)*(m*(x4+h/2*k41)A2/(x1+h/2*k11)*cos(x5+h/2*k51) -m*gc* . TOC o 1-5 h z cos(x5+h/2*k51)-m*gdelt*sin(x5+h/2*k51)*cos(x6+h/2*k61)+m*wA2*(x1+h/2 *k11).*cos(x3+h/2*k31).*(sin(x5+h/2*k51)*cos(x6+h/2*k61)*sin(x3+h/2*k31)+cos(x5+h/2*k51)*sin (x3+h/2*k31).+2*m*w*(x4+h/2*k41)*sin(x6+h/2*k

24、61)*cos(x3+h/2*k31);k62=1/(m*(x4+h/2*k41)*cos(x5+h/2*k51)*(m*(x4+h/2*k41)A2/(x1+h/2*k11) *cos(x5+h/2*k51)A2.*sin(x6+h/2*k61)*tan(x3+h/2*k31)-m*gdelt*sin(x6+h/2*k61)+m*wA2*(x1+h/ 2*k11).gdelt=g(2); k13=(x4+h/2sin(x6+h/2*k62)*sin(x3+h/2*k32)*cos(x3+h/2*k32)-2*m*w*(x4+h/2*k42)*( sin(x5+h/2*k52).*cos(x6

25、+h/2*k62)*cos(x3+h/2*k32)-cos(x5+h/2*k52)*sin(x3+h/2*k32);%ki4(i=1,2,.6)re=6387135;hi=x1+h*k13-re;vel=x4+h*k43;CL=0.5;Y=atmosdenty(hi,vel,CL);denty=Y(1);k42)*sin(x5+h/2*k52);k23=(x4+h/2*k42)*cos(x5+h/2*k52)*sin(x6+h/2*k62)/(x1+h/2*k12)*cos(x3 +h/2*k32);k33=(x4+h/2*k42)*cos(x5+h/2*k52)*cos(x6+h/2*k62

26、)/(x1+h/2*k12);k43=(-0.5*de nty*(x4+h/2*k42F2*s*cd-m*gc*si n(x5+h/2*k52)+m*gdelt*cos( x5+h/2*k52).*cos(x6+h/2*k62)-m*wA2*(x1+h/2*k12)*cos(x3+h/2*k32)*(cos(x5+h/2*k52)*cos(x6+h/2*k62)*sin(x3+h/2*k32)-sin(x5+h/2*k52)*cos(x3+h/2*k32)/m;k53=1/(m*(x4+h/2*k42)*(m*(x4+h/2*k42)A2/(x1+h/2*k12)*cos(x5+h/2*k52

27、) -m*gc* .cos(x5+h/2*k52)-m*gdelt*sin(x5+h/2*k52)*cos(x6+h/2*k62)+m*wA2*(x1+h/2 TOC o 1-5 h z *k12).*cos(x3+h/2*k32).*(sin(x5+h/2*k52)*cos(x6+h/2*k62)*sin(x3+h/2*k32)+cos(x5+h/2*k52)*sin(x3+h/2*k32).+2*m*w*(x4+h/2*k42)*sin(x6+h/2*k62)*cos(x3+h/2*k32);k63=1/(m*(x4+h/2*k42)*cos(x5+h/2*k52)*(m*(x4+h/2*

28、k42)A2/(x1+h/2*k12)*cos(x5+h/2*k52)A2.*sin(x6+h/2*k62)*tan(x3+h/2*k32)-m*gdelt*sin(x6+h/2*k62)+m*wA2*(x1+h/2*k12).ma=Y(3);kn=Y(2);cd=cdxs(ma,kn);wd=x3+h*k33;g=gravity(x1+h*k13,wd);gc=g(1);gdelt=g(2);k14=(x4+h*k43)*sin(x5+h*k53);k24=(x4+h*k43)*cos(x5+h*k53)*sin(x6+h*k63)/(x1+h*k13)*cos(x3+h*k33)k34=(

29、x4+h*k43)*cos(x5+h*k53)*cos(x6+h*k63)/(x1+h*k13);k44=(-0.5*de nty*(x4+h*k43f2*s*cd-m*gc*si n(x5+h*k53)+m*gdelt*cos(x5+h TOC o 1-5 h z *k53)*cos(x6+h*k63).-m*wA2*(x1+h*k13)*cos(x3+h*k33)*(cos(x5+h*k53)*cos(x6+h*k63)*sin(x3+h*k33)-sin(x5+h*k53)*cos(x3+h*k33)/m;k54=1/(m*(x4+h*k43)*(m*(x4+h*k43)A2/(x1+h

30、*k13)*cos(x5+h*k53)-m*gc*cos(x5+h*k53)-m*gdelt.*sin(x5+h*k53)*cos(x6+h*k63)+m*wA2*(x1+h*k13)*cos(x3+h*k33)*(sin(x5+h*k53)*cos(x6+h*k63)*sin(x3+h*k33)+cos(x5+h*k53)*sin(x3+h*k33)+2*m*w*(x4+h.*k43)*sin(x6+h*k63)*cos(x3+h*k33);k64=1/(m*(x4+h*k43)*cos(x5+h*k53)*(m*(x4+h*k43)A2/(x1+h*k13)*cos(x5+h*k53)A2

31、*sin(x6+h*k63).*tan(x3+h*k33)-m*gdelt*sin(x6+h*k63)+m*wA2*(x1+h*k13) *sin(x6+h*k63)*sin(x3+h*k33)*cos(x3+h*k33)-2*m*w*(x4+h*k43)*(sin(x5+h*k53)*cos(x6+h*k63).*cos(x3+h*k33)-cos(x5+h*k53)*sin(x3+h*k33);%计算下一步的值x1=x1+h*(k11+2*k12+2*k13+k14)/6;x2=x2+h*(k21+2*k22+2*k23+k24)/6;x3=x3+h*(k31+2*k32+2*k33+k34)/6;x4=x4+h*(k41+2*k42+2*k43+k44)/6;x5=x5+h*(k51+2*k52+2*k53+k54)/6;x6=x6+h*(k61+2*k62+2*k63+k64)/6;%形成向量r(i+1)=x1;la(i+1)=x2;del(i+1