飛行器建模與仿真

1、飛行器建模與仿真作業(yè)班級(jí)：姓名：飛設(shè)12潘周周學(xué)號(hào)：2110702015作業(yè)一、題目已知微分方程組廠為=X2X2=X3X3=X4X4=X5XX6X6皿a6-XXi-a X2Wo2X3_ 2Woa3X4a4X53_4WoWoa5X6_5Wo用四階龍哥庫(kù)塔法求解該方程組。程序1.編寫(xiě)rk4f.m函數(shù)fun ctio nr=rk4f(x1,x2,x3,x4,x5,x6,xs,ts,tf,N)%四階龍格庫(kù)塔法計(jì)算微分方程組程序%x1,x2,x3,x4,x5,x6分別是初值%xs是方程組中的X%ts，tf代表代表左端點(diǎn)和右端點(diǎn)%N是迭代步數(shù)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);%計(jì)算下一步的值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%畫(huà)出圖型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編寫(xiě)主函數(shù)調(diào)用該函數(shù)rk4f(0,0,0,0,0,0,1,0,0.5,100)三、運(yùn)行結(jié)果x1與t的關(guān)系0 6010.5、 題目0.50.5蒞 o0tLilttn作業(yè)2彈道式再入軌跡仿真。已知太空艙質(zhì)量為 350Kg,從近地點(diǎn)高度為200Km勺軌道返回。偏心率為0.2，軌道傾角為80度，近地點(diǎn)角距為265 度，升焦點(diǎn)赤經(jīng)100度。已知:r 0 二 6579.89967 km; (0) = -10 ；

7、 (0) =79.8489182889v(0) =7589 .30433867 m / s; (00.54681217 ; A(0) = 99 .955734求軌跡基本參數(shù) r(t), (t) (t), v(t), (t), A(t)。二、用龍格庫(kù)塔法編寫(xiě)程序求解此題在此函數(shù)中要求考慮非球形假設(shè)以及國(guó)際標(biāo)準(zhǔn)大氣。 所以在編寫(xiě)此函數(shù)之前應(yīng)先編寫(xiě)標(biāo)準(zhǔn)大氣函數(shù)、阻力函數(shù)以及重力加速度函數(shù)。大氣函數(shù)function Y=atmosdenty(h,vel,CL)%本函數(shù)用來(lái)求 den,Kn,ma%輸入 h,vel,Cl分別表示高度，速度，特征長(zhǎng)度R=287;%空氣氣體常數(shù)g0=9.806;%海平面重力加

8、速度T0=288.15;%海平面溫度re=6378140;%地球半徑gama=1.405;%海平面比熱比b=2/re;p0=1.01325e5;%海平面大氣壓layers=21;%大氣層數(shù)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;%首先計(jì)算出每一層端點(diǎn)處的密度值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% 在同溫層和變溫層采用不同公式計(jì)算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重力加速度函數(shù)function g=gravity(r,wd)%gc,gdelt 分別是徑向和切向重力加速度(當(dāng)?shù)厮阶鴺?biāo)系下的結(jié)果)； %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)%本函數(shù)用來(lái)計(jì)算阻力；%輸入mach表示ma數(shù)，kn表示克努増數(shù)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軌跡運(yùn)動(dòng)函數(shù)function gud=guidao(ts,tf,N)%本函數(shù)用來(lái)求解彈道式運(yùn)動(dòng)軌跡% N 迭代次數(shù)% ts 迭代初始時(shí)間% tf 迭代終止時(shí)間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);%計(jì)算下一步的值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