远程火箭飞行轨迹设计

由天下分享时间：2024/9/11 3:17:13 加入收藏我要投稿点赞

??+????

????????+?????????????? ?????=???????????????+??????

????

??????=??????+????????+???????????????????????????

????

???=?????????? ???=??????????

??=??(?????)

??????

222

??=√(??+??0??)+(??+??0??)+??0??

???????? ??= √????2??????2??+????2??????2??

?=????? ??=??0?????

????2

(2????????0) ??0=???????????? ????

???? ??=????????????√????2+????2{

4.3程序设计 4.3.1 符号定义

符号 ???? ???? ???? ???? ?? ?? ???? ???? ?? ?????? ?? 定义发机推力火箭参考面积阻力因数升力系数发射坐标系下横坐标发射坐标系下纵坐标地心坐标系下横坐标地心坐标系下纵坐标火箭总质量程序角符号 ?? ?? ?? ?? ??0 ?? ?? ? ???? ?? ???? 定义攻角侧滑角滚动角偏航角地心纬度轨道半径地球半径轨道高度引力加速度分量地球引力常数椭球扁率弹道倾角

9 / 18

注：在程序实现中，用alpha代替α，用theta代表θ，Phi_pr代表??????，Rho代表大气密度。 4.3.2 函数表

函数名

function phi=Phi(t) function rho=Density(h) function F1=F_v (Pe,m,q,theta,x,y,r) function F2=F_theta (Pe,alpha,theta,x,y,r,m,v) function F3=F_x (v,theta) function F4=F_y (v,theta) 功能

火箭所设定的程序角分段函数??????(??) USSA1976大气模型的分段密度函数三自由度仿真模型中 ??? 三自由度仿真模型中 ?? 三自由度仿真模型中 ??? 三自由度仿真模型中 ??? 4.3.3 程序框图开始根据实验对象及椭球模型等设定全局参数设定各变量的初值计算设定的程序角四阶龙格库塔算法迭代计算弹道参数根据工作时间计算火箭整体质量、改变燃料质量消耗速率、发动机推力根据所在高度计算大气密度更新计算地心失径、高度、攻角否发动机停机？是绘制图像，结果分析结束 4.3.4 程序代码

function rocket_1 clear all;close all;clc;

10 / 18

step=0.1; %仿真步长 time=186.4; t = 0:step:time;

totalnum =time/step+1; %总采样点数

%程序角参数

t1=10; t2=173; t3=186.4; fig=pi/11;

% 发动机参数 t_Eegine1=61.6; t_Eegine2=65.2; t_Eegine3=59.6;

m0=35.4*1000; %起飞质量 22.68吨 m_E1=22.78*1000;

m_E1_1=20.8*1000/61.6 ; % 一级发动机单位时间燃料质量消耗 m_E2=7.05*1000;

m_E2_1=6.25*1000/65.2 ; % 二级发动机单位时间燃料质量消耗 m_E3=3.65*1000;

m_E3_1=3.32*1000/59.6 ; % 二级发动机单位时间燃料质量消耗 Pe_1=912*1000; % 单位：N Pe_2=270*1000; % 单位：N Pe_3=155*1000; % 单位：N %方程参数

Cx=0.2; %阻力系数 g=-9.80665;

A_phi=25; %角度增益系数 S_M=1/4*pi*(1.67^2); R_Earth=6371*1000; %椭球模型

ae=6378245; %地球长轴 be=6356863; %地球短轴

A0=0.349; %发射方位角，弧度

B0=0.7176; %发射点地理纬度，弧度，酒泉

Phi_0=atan(be^2*tan(B0)/ae^2); %计算发射点地心纬度

R0=ae*be/sqrt(ae^2*(sin(Phi_0))^2+be^2*(cos(Phi_0))^2); %计算发射点地球半径 mu0=B0-Phi_0; %发射点地理纬度和地心纬度差 %计算发射点地球半径在方向投影 R_0x=-R0*sin(mu0)*cos(A0); R_0y=R0*cos(mu0);

R_0z=R0*sin(mu0)*sin(A0); % 初始状态 theta(1)=pi/2;

v(1)=0; x(1)=0; y(1)=0; r(1)=R_Earth; h(1)=0; m(1)=m0; Pe(1)=Pe_1; q(1)=0;

alphaa(1)=0;

Phi_pr(1)=pi/2; %程序角 v_dao(1)=0; theta_dao(1)=0; %大气模型

rho_s=1.2250; % kg/m3 海平面标准密度 rho(1)=rho_s;

% 火箭纵向运动方程

function F1=F_v (Pe,m,q,theta,x,y,r)

F1=( Pe- Cx*q*S_M+ m*g*(y+R_Earth)/r*sin(theta)+m*g*x/r*cos(theta) ) /m; end

function F2=F_theta (Pe,alpha,theta,x,y,r,m,v) if v==0

F2=0;

else F2=( Pe* alpha +m*g*(y+R_Earth)/r*cos(theta) -m*g*x/r*sin(theta) )/v/m;

11 / 18

end end

function F3=F_x (v,theta) F3=v*cos(theta); end

function F4=F_y (v,theta) F4=v*sin(theta); end

% 飞行程序角分段函数 function phi=Phi(t) if t

phi=pi/2; else if t

phi=pi/2+(pi/2-fig)*( (t-t1)^2 / (t2-t1)^2- 2*(t-t1)/(t2-t1) ); else

phi=fig; end end end

% 计算程序角

for i=1:1:(totalnum-1)

Phi_pr(i+1)=Phi(i*step); end

function rho=Density(h) %密度函数 H=h/1000; if H<=0

rho=(1-9.6005765*0.01*H)*rho_s; elseif (0

rho=rho_s*W^4.25588;

elseif (11.0191

elseif (20.0631

rho=3.2722*0.01*(W^(-35.1629))*rho_s; elseif (32.1619

rho=3.2618*0.001*(W^(-13.2011))*rho_s; elseif (47.3501

rho=2.5280*0.0001*(W^(11.2011))*rho_s; elseif (71.8020

rho=1.7632*0.00001*(W^(16.0816))*rho_s; elseif (86.0000

rho=3.6411*0.000001*W*rho_s; elseif (91.0000200 rho=0; end end

% 用四阶Runge-Kutta求解纵向运动方程 for i=1:1:(totalnum-1)

a1=F_v (Pe(i),m(i),q(i),theta(i),x(i),y(i),r(i));

b1=F_theta (Pe(i),alphaa(i),theta(i),x(i),y(i),r(i),m(i),v(i)); c1=F_x (v(i),theta(i)); d1=F_y (v(i),theta(i)); v_dao(i+1)=a1;

12 / 18

theta_dao(i+1)=b1;

a2=F_v (Pe(i),m(i),q(i),theta(i),x(i),y(i),r(i));

b2=F_theta (Pe(i),alphaa(i),theta(i)+step/2*b1,x(i),y(i),r(i),m(i),v(i)); c2=F_x (v(i)+step/2*a1,theta(i)+step/2*b1); d2=F_y(v(i)+step/2*a1,theta(i)+step/2*b1);

a3=F_v (Pe(i),m(i),q(i),theta(i),x(i),y(i),r(i));

b3=F_theta (Pe(i),alphaa(i),theta(i)+step/2*b2,x(i),y(i),r(i),m(i),v(i)); c3=F_x (v(i)+step/2*a2,theta(i)+step/2*b2); d3=F_y (v(i)+step/2*a2,theta(i)+step/2*b2);

a4=F_v (Pe(i),m(i),q(i),theta(i),x(i),y(i),r(i));

b4=F_theta (Pe(i),alphaa(i),theta(i)+step/2*b3,x(i),y(i),r(i),m(i),v(i)); c4=F_x (v(i)+step/2*a3,theta(i)+step/2*b3); d4=F_y(v(i)+step/2*a3,theta(i)+step/2*b3);

v(i+1)=v(i)+step*(a1+2*a2+2*a3+a4)/6;

theta(i+1)=theta(i)+step*(b1+2*b2+2*b3+b4)/6; x(i+1)=x(i)+step*(c1+2*c2+2*c3+c4)/6; y(i+1)=y(i)+step*(d1+2*d2+2*d3+d4)/6;

if (i*step)<=t_Eegine1

m(i+1)=m0-m_E1_1*(i*step); Pe(i+1)=Pe_1;

elseif ((i*step)>t_Eegine1) && ( (i*step)<(t_Eegine1+t_Eegine2) ) m(i+1)=m0-m_E1-m_E2_1*(i*step-t_Eegine1); Pe(i+1)=Pe_2; elseif ( (i*step)

m(i+1)=m0-m_E1-m_E2-m_E3_1*(i*step-t_Eegine1-t_Eegine2); Pe(i+1)=Pe_3; elseif i==(totalnum-1)

m(i+1)=m0-m_E1-m_E2-m_E3; Pe(i+1)=0; end

rho(i+1)=Density(h(i)); q(i+1)=1/2*rho(i)*v(i)^2;

r(i+1)=sqrt( (x(i)+R_0x)^2+(y(i)+R_0y)^2+(0+R_0z)^2); h(i+1)=r(i)-R_Earth;

alphaa(i+1)=A_phi*(Phi_pr(i+1) - theta(i+1) ); end

set(0,'defaultfigurecolor','w'); %设置图形背景为白色 figure

plot(t,v,'black','LineWidth',1.0); xlabel('t /s');