飞行器建模与仿真.docx
- 文档编号:6603905
- 上传时间:2023-01-08
- 格式:DOCX
- 页数:20
- 大小:283.40KB
飞行器建模与仿真.docx
《飞行器建模与仿真.docx》由会员分享,可在线阅读,更多相关《飞行器建模与仿真.docx(20页珍藏版)》请在冰豆网上搜索。
飞行器建模与仿真
飞行器建模与仿真作业
班级:
飞设12
姓名:
潘周周
学号:
2110702015
作业一
一、题目
已知微分方程组
用四阶龙哥库塔法求解该方程组。
二、程序
1.编写rk4f.m函数
functionr=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=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;
fori=1:
1:
N;
%ki1(i=1,2,...6)
k11=x2;
k21=x3;
k31=x4;
k41=x5;
k51=x6;
k61=w0^6/a6*(xs-x1-a1*x2/w0-a2*x3/w0^2-a3*x4/w0^3-a4*x5/w0^4-a5*x6/w0^5);
%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=w0^6/a6*(xs-(x1+h/2*k11)-a1*k12/w0-a2*k22/w0^2-a3*k32/w0^3-a4*k42/w0^4-a5*k52/w0^5);
%%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=w0^6/a6*(xs-(x1+h/2*k12)-a1*k13/w0-a2*k23/w0^2-a3*k33/w0^3-a4*k43/w0^4-a5*k53/w0^5);
%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=w0^6/a6*(xs-(x1+h*k13)-a1*k14/w0-a2*k24/w0^2-a3*k34/w0^3-a4*k44/w0^4-a5*k54/w0^5);
%计算下一步的值
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;
%形成向量
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)
plot(t,aa,'-')
subplot(2,3,2)
plot(t,ab,'-')
subplot(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)
三、运行结果
作业2
一、题目
弹道式再入轨迹仿真。
已知太空舱质量为350Kg,从近地点高度为200Km的轨道返回。
偏心率为0.2,轨道倾角为80度,近地点角距为265度,升焦点赤经100度。
已知:
求轨迹基本参数
。
2、用龙格库塔法编写程序求解此题
在此函数中要求考虑非球形假设以及国际标准大气。
所以在编写此函数之前应先编写标准大气函数、阻力函数以及重力加速度函数。
(1)大气函数
functionY=atmosdenty(h,vel,CL)
%本函数用来求den,Kn,ma
%输入h,vel,Cl分别表示高度,速度,特征长度
R=287;%空气气体常数
g0=9.806;%海平面重力加速度
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=[Mo28.96428.96428.96428.96428.96428.96228.96228.88028.5628.0726.9226.66...
26.525.8524.6922.6619.9417.9416.8416.1716.17];
%高度分层所得向量
z=[011.019120.063132.161947.350151.412571.802086100110120150160170190...
2303004005006007002000]*1e3;
%温度层向量
T=[T0216.65216.65228.65270.65270.65214.65186.946210.65260.65360.65960.651110.6...
1210.651350.651550.651830.652160.652420.652590.6527002700];
%热量递减率
lr=[-6.5e-301e-32.8e-30-2.8e-3-2e-31.693e-35e-31e-22e-21.5e-21e-27e-35e-34e-3...
3.3e-32.6e-31.7e-31.1e-30];
rou0=p0/(R*T0);
%p
(1)=p0;
T
(1)=T0;
rou=zeros(22);
rou
(1)=rou0;
%首先计算出每一层端点处的密度值
fori=1:
layers
if~(lr(i)==0)
rou(i+1)=rou(i)*(T(i+1)/T(i))^(-(1+(1+b*(T(i)/lr(i)-z(i)))*g0/(R*(lr(i)))))......
*exp(b*g0/(R*lr(i))*(z(i+1)-z(i)));
else
rou(i+1)=rou(i)*exp(-g0*(z(i+1)-z(i))*(1-b*(z(i+1)+z(i))/2)/(R*T(i)));
end
end
%在同温层和变温层采用不同公式计算
fori=1:
layers
ifh 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)))*g0/(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)))); end c=sqrt(gama*R*TM); ma=vel/c; 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*sigma^2*Vm; L=Vm/F; Kn=L/CL; Y=[denKnma]; return; end end (2)重力加速度函数 functiong=gravity(r,wd) %gc,gdelt分别是径向和切向重力加速度(当地水平坐标系下的结果); %r,wd分别表示距离和纬度; ywd=pi/2-wd; mu=3.9860004e14; re=6387.135e3; J2=1.0826e-3; J3=2.532153e-7; J4=1.6109876e-7; gc=-(mu*(1-1.5*J2*(re/r)^2*(3*cos(ywd)^2-1)-2*J3*cos(ywd)*(5*cos(ywd)^2-3)... *(re/r)^3-(5/8)*J4*(35*cos(ywd)^4-30*cos(ywd)^2+3)*(re/r)^4)/r^2); gdelt=-(-3*mu*sin(ywd)*cos(ywd)*(re/r)^2*(J2+0.5*J3*(5*cos(ywd)^2-1)*(re/r)/cos(ywd)... +5/6*J4*(7*cos(ywd)^2-1)*(re/r)^2)/r^2); g=[gc,gdelt]; (3)大气阻力 functioncd=cdxs(mach,kn) %本函数用来计算阻力; %输入mach表示ma数,kn表示克努増数 mar=[00.20.30.40.50.60.80.90.951.051.11.21.62.02.533.851099]; cdr=[0.4754750.4754750.475760.483360.4889650.5083450.565630.6181650.668135... 1.0317951.017070.9905650.8159550.692360.609710.546060.5130.4940.483170.48317]; cdc=interp1(mar,cdr,mach); ifkn<0.0146 cd=cdc; elseifkn>14.5 cdfm=1.75+sqrt(pi)/(2*mach*sqrt(1.41/2)); cd=cdfm; else cdfm=1.75+sqrt(pi)/(2*mach*sqrt(1.41/2)); cd=cdc+(cdfm-cdc)*(1/3*log10(kn/(1/2))+0.5113); end (4)轨迹运动函数 functiongud=guidao(ts,tf,N) %本函数用来求解弹道式运动轨迹 %N迭代次数 %ts迭代初始时间 %tf迭代终止时间 formatlong r=zeros(1,N+1); la=zeros(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)=99.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; fori=1: N re=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); 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; k41=(-0.5*denty*x4^2*s*cd-m*gc*sin(x5)+m*gdelt*cos(x5)*cos(x6)-m*w^2*x1*cos(x3)*(cos(x5)... *cos(x6)*sin(x3)-sin(x5)*cos(x3)))/m; k51=1/(m*x4)*(m*x4^2/x1*cos(x5)-m*gc*cos(x5)-m*gdelt*sin(x5)*cos(x6)+m*w^2*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*x4^2/x1*cos(x5)^2*sin(x6)*tan(x3)-m*gdelt*sin(x6)+m*w^2*x1... *sin(x6)*sin(x3)*cos(x3)-2*m*w*x4*(sin(x5)*cos(x6)*cos(x3)-cos(x5)*sin(x3))); %ki2(i=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=(x4+h/2*k41)*cos(x5+h/2*k51)*cos(x6+h/2*k61)/(x1+h/2*k11); k42=(-0.5*denty*(x4+h/2*k41)^2*s*cd-m*gc*sin(x5+h/2*k51)+m*gdelt*cos(x5+h/2*k51)... *cos(x6+h/2*k61)-m*w^2*(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; k52=1/(m*(x4+h/2*k41))*(m*(x4+h/2*k41)^2/(x1+h/2*k11)*cos(x5+h/2*k51)-m*gc*... cos(x5+h/2*k51)-m*gdelt*sin(x5+h/2*k51)*cos(x6+h/2*k61)+m*w^2*(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*k61)*cos(x3+h/2*k31)); k62=1/(m*(x4+h/2*k41)*cos(x5+h/2*k51))*(m*(x4+h/2*k41)^2/(x1+h/2*k11)*cos(x5+h/2*k51)^2... *sin(x6+h/2*k61)*tan(x3+h/2*k31)-m*gdelt*sin(x6+h/2*k61)+m*w^2*(x1+h/2*k11)... *sin(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,CL); 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); gdelt=g (2); k13=(x4+h/2*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)/(x1+h/2*k12); k43=(-0.5*denty*(x4+h/2*k42)^2*s*cd-m*gc*sin(x5+h/2*k52)+m*gdelt*cos(x5+h/2*k52)... *cos(x6+h/2*k62)-m*w^2*(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)^2/(x1+h/2*k12)*cos(x5+h/2*k52)-m*gc*... cos(x5+h/2*k52)-m*gdelt*sin(x5+h/2*k52)*cos(x6+h/2*k62)+m*w^2*(x1+h/2*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*k42)^2/(x1+h/2*k12)*cos(x5+h/2*k52)^2... *sin(x6+h/2*k62)*tan(x3+h/2*k32)-m*gdelt*sin(x6+h/2*k62)+m*w^2*(x1+h/2*k12)... *sin(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+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); 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=(x4+h*k43)*cos(x5+h*k53)*cos(x6+h*k63)/(x1+h*k13); k44=(-0.5*denty*(x4+h*k43)^2*s*cd-m*gc*sin(x5+h*k53)+m*gdelt*cos(x5+h*k53)*cos(x6+h*k63)... -m*w^2*(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)^2/(x1+h*k13)*cos(x5+h*k53)-m*gc*cos(x5+h*k53)-m*gdelt... *sin(x5+h*k53)*cos(x6+h*k63)+m*w^2*(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)^2/(x1+h*k13)*cos(x5+h*k53)^2*sin(x6+h*k63)... *tan(x3+h*k33)-m*gdelt*sin(x6+h*k63)+m*w^2*(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)... *co
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 飞行器 建模 仿真