??+????
????????+?????????????? ?????=???????????????+??????
????
??????=??????+????????+???????????????????????????
????
???=?????????? ???=??????????
??=??(?????)
??????
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.0000 % 用四阶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)
远程火箭飞行轨迹设计
![](/skin/haowen/images/icon_star.png)
![](/skin/haowen/images/icon_star.png)
![](/skin/haowen/images/icon_star.png)
![](/skin/haowen/images/icon_star.png)