固体火箭弹道

最新推荐文章于 2024-09-12 08:51:28 发布

___xh___

最新推荐文章于 2024-09-12 08:51:28 发布

阅读量1.2k

点赞数 31

文章标签： matlab

本文链接：https://blog.csdn.net/djzdym/article/details/137874146

版权

close all;clear all;clc;clf;
global m0 dand phi1k th t1 t2 ta1 t1k ta2 t2k ta3 t3k t4h t4k;
th=200;t1=5;t2=8;ta1=86.4;t1k=ta1+5;ta2=t1k+81.7;t2k=ta2+5;ta3=t2k+70.5;t3k=ta3+5;
t4h=t3k+th;t4k=t4h+176;
options=odeset(‘MaxStep’,0.5,‘RelTol’,1e-4,‘AbsTol’,1e-6);
tspan=(0:0.001:620) ; %t区间步长
% Y= 【v，theat，x，y mass】
m0=(170+60+18+3.5)1000;%初始质量
y0=[1 pi/2 0 5 m0];%初始条件
[t,y]=ode45(@funt4,tspan,y0,options);
%%
% y2=[t,ma,betak,Theatk57.3,rk,ny,alpha];
% figure(‘numbertitle’,‘off’,‘name’,‘y–x曲线’)
% plot(y(:,3),y(:,4),‘b-’,‘lineWidth’,1.5);grid on;xlabel(‘x/m’);ylabel(‘y/m’);title(‘y–x曲线’);
figure(‘numbertitle’,‘off’,‘name’,‘y–t曲线’)
plot(t,y(:,4),‘b-’,‘lineWidth’,1.5);grid on;xlabel(‘t/s’);ylabel(‘y/m’);title(‘y–t曲线’);
figure(‘numbertitle’,‘off’,‘name’,‘v–t曲线’)
plot(t,y(:,1)/1000,‘b-’,‘lineWidth’,1.5);grid on;xlabel(‘t/s’);ylabel(‘v/km/s’);title(‘y–t曲线’);

figure(‘numbertitle’,‘off’,‘name’,‘theat–t曲线’)
plot(t,y(:,2)*57.3,‘b-’,‘lineWidth’,1.5);grid on;xlabel(‘t/s’);ylabel(‘theat/°’);
hold on;
plot(dand(:,1),dand(:,9),‘r-’,‘lineWidth’,1.5);grid on;xlabel(‘t/s’);ylabel(‘phi/°’);hold on;
plot(dand(:,1),dand(:,8),‘b+’,‘lineWidth’,1.5);grid on;xlabel(‘t’);ylabel(‘a/°’);title(‘pta–t曲线’);
legend(“theat”,“phi”,“a”)

figure(‘numbertitle’,‘off’,‘name’,‘x–y曲线’)
plot(y(:,3)/1e3,y(:,4)/1e3,‘b-’,‘lineWidth’,1.5);grid on;xlabel(‘x’);ylabel(‘y’);title(‘x–y曲线’);

figure(‘numbertitle’,‘off’,‘name’,‘masss–t曲线’)
plot(t,y(:,5),‘b-’,‘lineWidth’,1.5);grid on;xlabel(‘t/s’);ylabel(‘mass’);title(‘mass–t曲线’);
% y2=[t,y(1),rk,Theatk,ma,betak,ny,alpha57.3,phi];
figure(‘numbertitle’,‘off’,‘name’,‘ny–t曲线’)
plot(dand(:,1),dand(:,7),‘b-’,‘lineWidth’,1.5);grid on;xlabel(‘t/s’);ylabel(‘ny’);title(‘ny–t曲线’);
figure(‘numbertitle’,‘off’,‘name’,‘a–t曲线’)
plot(dand(:,1),dand(:,8),‘b-’,‘lineWidth’,1.5);grid on;xlabel(‘t’);ylabel(‘a/°’);title(‘a–t曲线’);
figure(‘numbertitle’,‘off’,‘name’,‘thdd–t曲线’)
plot(dand(:,1),dand(:,4)53.7,‘b-’,‘lineWidth’,1.5);grid on;xlabel(‘t’);ylabel(‘thdd/°’);title(‘thdd–t曲线’);
function dydt=funt4(t,y)
% 发射坐标系火箭弹道主动段计算程序
% 近似处理：地球为均质圆球，忽略地球扁率及gφ，we=0;忽略附加力；力矩瞬时平衡；忽略二阶小量；
% 状态变量5个： Y= 【v，theat，x，y mass】 dydt=Y’
dydt=zeros(5,1); %状态向量
%%
format long;
global dand rk3 vk3 tvlpk3 phi1k m0 th t1 t2 ta1 t1k ta2 t2k ta3 t3k t4h t4k phik2;
%% 上游专业输入值
% 总体专业输入火箭参数
S=3.143.35^2/4;%特征面积
R=6371110;
pid=57.295779513082323;
m0=(170+60+18+3.5)1000;%初始质量
% 控制专业输入
a0phi=1;
% 动力专业输入火箭动力数据动力专业根据级确定推力
% Pe= 45e39.8;%发动机等效推力，P-X1C
% 气动专业输入气动力源数据，Cx=f(ma)=f1|alpha=0+f2|alpha^2
% Cy_a=f(ma)
ma0=[0.6,0.8,0.9,1.0,50];
Cx10=[-0.2000,-0.2210,-0.2315,-0.2420,-0.2500];
Cx20=[-0.0005,-0.0005,-0.0005,-0.0005,-0.0005,];
Cya0=[0.2505,0.253,0.254,0.255,0.256];
Cydetaz0=[0.04,0.047,0.043,0.045,0.045];
mza0=[-0.0290,-0.0298,-0.0300,-0.0302,-0.0302];
mzdetaz0=[-0.0220,-0.0240,-0.0232,-0.0230,-0.0234];
%
h=y(4);%计算当前高度→马赫数→气动插值→气动力→方程组；
fqx=airCond(y(4));
cair=fqx(4); %计算当前大气声速
rhoair=fqx(3); %计算当前大气密度
ma=y(1)/cair; %马赫数 y1 是速度
q=0.5rhoairy(1)y(1); %动压
% 程序角优化参数 phi2k二级末俯仰角度，phi3k三级末俯仰角，th四级初始滑行时间。
% am最大攻角，a-am时间调节因子
% 飞行时序，受限于发动机
am=6.87/pid;a=0.023;phi2k=12/pid;phi3k=10/pid;
%计算程序攻角
if t<t1 % 垂直段|10
alpha=0;
X=0;
Y=0;
phi=(alpha+y(2));
dydt(5)=-170/0.920.081000/86.4;
Pe=4525e3;
elseif t<=t2 %负攻角转弯段|0.8ma
alpha=4amexp(a*(t1-t))(exp(a(t1-t))-1);
phi=(alpha+y(2));
dydt(5)=-1700.921000/86.4;
Pe=4525e3;
elseif t<=ta1 %弹道重力转弯段ta1发动机主动段关机
alpha=0;
phi=(alpha+y(2));
dydt(5)=-1700.921000/86.4;
Pe=4525e3;
elseif t<=t1k %分离段，t1k一级分离，t3为一级关机点
alpha=0;
phi=(alpha+y(2));
if abs(t-t1k)<0.5
phi1k=(alpha+y(2));
y(5)=m0-170e3;
end
dydt(5)=-1700.081000/5;
Pe=0;
elseif t<=ta2 %二级主动段
phi=phi1k+(phi2k-phi1k)(t-t1k)/(ta2-t1k);
phi=y(2);
alpha=phi-y(2);
dydt(5)=-600.92/81.71e3;
Pe=1925e3;
elseif t<=t2k %二级分离滑行段
phi=y(2);
if abs(t-t2k)<0.5
phik2=y(2);
end
dydt(5)=-600.08/51e3;
Pe=0;
alpha=phi-y(2);
elseif t<=ta3 %三级主动段
phi=y(2);
alpha=(phi-y(2));
dydt(5)=-180.9/70.51e3;
Pe=655e3;
elseif t<=t3k %三级分离
phi=y(2);
alpha=phi-y(2);
if abs(t3k-t)<0.5
rk3=R+y(4);
vk3=y(1);
tvlpk3=y(2)+atan2(y(3),R+y(4));
end
Pe=0;
dydt(5)=-180.1/5*1e3;
elseif t<=t4h %四级滑行段
% [Et4,vt4,rt4,tlvpt4]= kpl(t3k,rk3,vk3,tvlpk3,t);
Pe=0;
X=0;
Y=0;
dydt(5)=0;
phi=y(2);
alpha=phi-y(2);

elseif t<=t4k %四级主动段
Pe=125e3;
dydt(5)=-3.50.9/71.81e3;
X=0;
Y=0;
phi=y(2);
alpha=phi-y(2);
end
%% 气动计算
if y(4)<80e3

[cl,cd,cma]=ClCdCma(ma,alpha);
X=q*S*cd/5;
%升力计算
Y=q*S*cl;

else
X=0
Y=0;
end
%%
rk=sqrt((R+y(4))^2+y(3)y(3));
f=6.67e-11;%引力常数
M=5.965e24;%地球质量
g=-fM/rk/rk; %计算当前重力加速度，引力y轴假设
% 积分计算
dydt(1)=(Pe-X)/y(5)+gsin(y(2));
if 190<t && t<=250
dydt(2)=(Pealpha+Y)/y(5)/y(1)+g/y(1)cos(y(2))-0.1/57;
else
dydt(2)=(Pealpha+Y)/y(5)/y(1)+g/y(1)cos(y(2));
end
dydt(3)=y(1)cos(y(2));
dydt(4)=y(1)sin(y(2));
betak=atan2(y(3),R+y(4));
Theatk=y(2)+betak;%y(2)+atan2(y(3),R+y(4))
rk=sqrt((R+y(4))^2+y(3)y(3));
ny=(Pealpha)/y(5)/9.8;
phi=phi57.3;
y2=[t,y(1),rk,Theatk,ma,betak,ny,alpha57.3,phi];
dand=[dand;y2];
end
function [ y ] = airCond(h)
%UNTITLED5 此处显示有关此函数的摘要
% 输入 h，m；输出t p rho a
h=h/1000;
R0=6371110/1000;
rhosl=1.225;
psl=101325;
z=1/(1/h-1/R0);
if z<=11.0191
w=1-h/44.3308;
T=288.15w;
p=pslw^5.2559;
rho=rhoslw^4.2559;
elseif z<=20.0631
w=exp((14.9647-h)/6.3416);
T=216.65;
p=pslw0.11953;
rho=rhoslw0.15898;
elseif z<=32.1619
w=1+(h-24.9021)/221.552;
T=221.552w;
p=pslw^(-34.1629)0.025158;
rho=rhoslw^(-35.1629)0.032722;
elseif z<=47.3501
w=1+(h-39.7499)/89.4107;
T=250.350w;
p=pslw^(-12.2011)0.0028338;
rho=rhoslw^(-13.2011)0.0032618;
elseif z<=51.4125
w=exp((48.6252-h)/7.9223);
T=270.650;
p=pslw0.00089155;
rho=rhoslw0.0009492;
elseif z<=71.8020
w=1-(h-59.439)/88.2218;
T=247.021w;
p=pslw^12.20110.00021671;
rho=rhoslw^11.20110.00025280;
elseif z<=86
w=1-(h-78.0303)/100.2950;
T=200.59w;
p=pslw^17.08160.000012274;
rho=rhoslw^16.08160.000017632;
elseif z<=91
w=exp((87.2848-h)/5.47);
T=186.87;
p=psl*(2.273+1.042e-3h)1e-6w;
rho=rhoslw3.6411e-6;
else
T=0;
p=psl0;
rho=rhosl0;
end
a=20.0468T^0.5;
y=[T,p,rho,a];
end
function[cl,cd,cma]=ClCdCma(ma,alpha)
% qdmat=xlsread(“C:\Users\47382\Desktop\dat1.xlsx”);
alphaqd=[-10 -6 -2 2 6 10];
maqd=[0.1 1.5 3 4.5 6 7.5];
% clqd=qdmat(:,3);
% cdqd=qdmat(:,4);
% cmaqd=qdmat(:,7);
% clqd=reshape(clqd,[length(alphaqd),length(maqd)]);
% clqd=clqd’;
% cdqd=reshape(cdqd,[length(alphaqd),length(maqd)]);
% cdqd=cdqd’;
%
% cmaqd=reshape(cmaqd,[length(alphaqd),length(maqd)]);
% cmaqd=cmaqd’;
cdqd=[0.234800000000000,0.151800000000000,0.116700000000000,0.116700000000000,0.151800000000000,0.234800000000000;0.440000000000000,0.319100000000000,0.268100000000000,0.268100000000000,0.319100000000000,0.440000000000000;0.447700000000000,0.316500000000000,0.269000000000000,0.269000000000000,0.316500000000000,0.447700000000000;0.466500000000000,0.338100000000000,0.292000000000000,0.292000000000000,0.338100000000000,0.466500000000000;0.490600000000000,0.363600000000000,0.318400000000000,0.318400000000000,0.363600000000000,0.490600000000000;1.46050000000000,1.36320000000000,1.33310000000000,1.33310000000000,1.36320000000000,1.46050000000000];
clqd=[-0.793500000000000,-0.420600000000000,-0.117900000000000,0.117900000000000,0.420600000000000,0.793500000000000;-0.918800000000000,-0.445700000000000,-0.102100000000000,0.102100000000000,0.445700000000000,0.918800000000000;-0.900300000000000,-0.409200000000000,-0.0899000000000000,0.0899000000000000,0.409200000000000,0.900300000000000;-0.878000000000000,-0.395100000000000,-0.0852000000000000,0.0852000000000000,0.395100000000000,0.878000000000000;-0.865400000000000,-0.386900000000000,-0.0823000000000000,0.0823000000000000,0.386900000000000,0.865400000000000;-0.689900000000000,-0.278600000000000,-0.0457000000000000,0.0457000000000000,0.278600000000000,0.689900000000000];
cmaqd=[0.00120000000000000,0.00220000000000000,0.00360000000000000,0.00360000000000000,0.00220000000000000,0.00120000000000000;-0.00190000000000000,0.00210000000000000,0.00530000000000000,0.00530000000000000,0.00210000000000000,-0.00190000000000000;0.00230000000000000,0.00550000000000000,0.00830000000000000,0.00830000000000000,0.00550000000000000,0.00230000000000000;0.00310000000000000,0.00640000000000000,0.00920000000000000,0.00920000000000000,0.00640000000000000,0.00310000000000000;0.00340000000000000,0.00680000000000000,0.00970000000000000,0.00970000000000000,0.00680000000000000,0.00340000000000000;0.00350000000000000,0.00700000000000000,0.0100000000000000,0.0100000000000000,0.00700000000000000,0.00350000000000000];

cl=interp2(maqd,alphaqd,clqd,ma,alpha,‘spline’);

cd=interp2(maqd,alphaqd,cdqd,ma,alpha,‘spline’);

cma=interp2(maqd,alphaqd,cmaqd,ma,alpha,‘spline’);
end
%%
% vk=dand(end,2);
% rk=dand(end,3);
% Thk=dand(end,4);
% Thk=pi/4;
% f=6.67e-11;%引力常数
% M=5.965e24;%地球质量
% u=fM;%引力常数
% %g=fM/r/r;
% uk=vkvk/(u/rk);
% P=rkukcos(Thk)cos(Thk);
% e=sqrt(1+uk(uk-2)cos(Thk)cos(Thk))
% a=-urk/(rkvkvk-2*u)%长轴跟机械能有关系
% %