飞行方案设计大作业1.docx
《飞行方案设计大作业1.docx》由会员分享，可在线阅读，更多相关《飞行方案设计大作业1.docx（15页珍藏版）》请在冰点文库上搜索。
飞行方案设计大作业1
航天飞行动力学大作业
韩谨阳
2015300464
1、方案飞行2、弹道设计3、卫星摄动与机动
第一部分
飞行方案
卫星的摄动与机动
第三部分
弹道设计
第二部分
飞行方案大作业
一、问题描述
在已知导弹质量、转动惯量、发动机推力等参数的情况下，导弹分为三个飞行方案，即三个阶段飞行。
阶段一：
飞行距离在
，采用追踪法，其中方案高度与距离的关系、方案弹道倾角与高度的关系如下：
（1）
阶段二：
飞行距离在
，采用追踪法，其中方案高度与距离的关系、方案弹道倾角与高度的关系、导弹因燃料消耗而质量改变参数如下：
（2）
（3）
阶段三：
飞行方案
，而最终目标位置为
采用比例导引法
（4）
要求：
1）计算纵向理想弹道，给出采用瞬时平衡假设
时所有纵向参数随时间的变化曲线。
2）不考虑气动力下洗影响，计算飞行器沿理想弹道飞行时，你认为可以作为特性点的5个以上点处的纵向短周期扰动运动的动力系数，并分析其在特性点处的自由扰动的稳定性，以及计算在各个特性点处弹体传递函数
。
二、建立模型
基于“瞬时平衡”假设，导弹在铅垂平面内运动的质心运动方程组为：
（5）
因为阶段一不考虑导弹质量随时间的变化，因此阶段一的模型需要联立公式
（1）、公式（5）；
其中攻角
可根据瞬时平衡假设
从而可得到导弹攻角与弹道倾角之间的关系
（6）
其中
（7）
其中假设公式
（1）的
中的
又因为阶段二需要考虑导弹质量随时间的变化，因此阶段二的模型需要联立公式
（2）公式（5）、公式（6）、公式（7）
最后一阶段，因为利用了比例导引法
公式（4）的k=2，可得导弹到达目标的相对微分方程为
而导引率
、其中k=2；
因为第三阶段的初始参数及终点坐标均为直角坐标系，由下图可知将
代入到公式（4），得到直角坐标系下的微分方程组
另外补充方程法向平衡方程：
三、算法实现
编程使用MATLAB软件，并运用欧拉方程解微分方程，即ode45函数；
四、程序源代码
*************************阶段一******************************
functiondy=jieduan1（t,y）
dy=zeros（4,1）;
m=320;
g=9.8;
P=2000;
q=0.5*1.2495*（（288.15-0.0065*y（4））/288.15）.^4.2558*y
（1）.^2;
k=-9;
dk=-0.5;
Hi=2000*cos（0.000314*1.1*y（3））+5000;
dHi=-2000*0.000314*1.1*sin（y（3））;
delta=k*（y（4）-Hi）+dk*（dy（3）-dHi）;
alpha=0.34*delta;
Xb=（0.2+0.005*alpha^2）*q*0.45;
Yb=（0.25*alpha+0.05*delta）*q*0.45;
dy=zeros（4,1）;
dy
（1）=P*cos（alpha）/m-Xb/m-g*sin（y
（2））;
dy
（2）=P*sin（alpha）/m/y
（1）+Yb/m/y
（1）-g*cos（y
（2））/y
（1）;
dy（3）=y
（1）*cos（y
（2））;
dy（4）=y
（1）*sin（y
（2））;
end
******************************阶段二******************************
functiondy=jieduan2（t,y）
dy=zeros（4,1）;
m=320-0.46*t;
g=9.8;
P=2000;
q=0.5*1.2495*（（288.15-0.0065*y（4））/288.15）.^4.2558*y
（1）.^2;
k=-0.25;
Hi=3050;
delta=k*（y（4）-Hi）;
alpha=0.34*delta;
Xb=（0.2+0.005*alpha^2）*q*0.45;
Yb=（0.25*alpha+0.05*delta）*q*0.45;
dy
（1）=P*cos（alpha/180*pi）/m-Xb/m-g*sin（y
（2）/180*pi）;
dy
（2）=P*sin（alpha/180*pi）/m/y
（1）+Yb/m/y
（1）-g*cos（y
（2）/180*pi）/y
（1）;
dy（3）=y
（1）*cos（y
（2）/180*pi）;
dy（4）=y
（1）*sin（y
（2）/180*pi）;
end
*******************************阶段三********************************
functiondy=jieduan3（t,y）
v=y（4）;
k=10;
m=285.04-0.46*t;
q0=-atan（3050/6000）;
g=9.8;
q1=0.5*1.2495*（（288.15-0.0065*y
（2））/288.15）.^4.2558*y（4）.^2;
k1=10;
dk1=0.05;
dy=zeros（4,1）;
r=sqrt（y
（1）^2+y
（2）^2）;
q=atan（y
（2）/（y
（1）-30000））;
elta=q-y（3）;
dr=-v*cos（elta）;
tht=q0+k*（q-q0）;
dq=v/r*sin（elta）;
dtht=k*dq;
delta=k1*（y（3）-tht）+dk1*（dy（3）-dtht）;
alpha=0.34*delta;
dy
（1）=-dr*cos（q）+r*sin（q）*dq;
dy
（2）=-dr*sin（q）-r*cos（q）*dq;
Yb=（0.25*alpha+0.05*delta）*q1*0.45;
dy（3）=（2000*sin（alpha）/m+Yb/m-g*cos（y（3）））/v;
y（4）=v;
end
***********************************main函数************************************
m
（1）=287.2204;%导弹质量
P=2000;%发动机推力
g=9.8;
k=5;
det
（1）=0.045;
a
（1）=0.6186;
sit
（1）=-0.000002024;
V
（1）=217.2867;%初始速度
x
（1）=24000;%初始位置
H
（1）=3071;%初始高度
H1
（1）=3050;
S=0.45;%参考面积
L=2.5;%参考长度
k1=-0.14;
k2=-0.06;
sit1
（1）=sit
（1）;
p0=1.2495;
T0=288.15;
T
（1）=T0-0.0065*H
（1）;
p
（1）=p0*（T
（1）/T0）^4.25588;
q
（1）=1/2*p
（1）*V
（1）^2;%大气密度计算公式
Cx
（1）=0.2+0.005*a
（1）^2;
Cy
（1）=0.25*a
（1）+0.05*det
（1）*180/pi;%升力系数
Y
（1）=Cy
（1）*q
（1）*S;
X
（1）=Cx
（1）*q
（1）*S;
SIT
（1）=（P*sind（a
（1））+（Y
（1）-m
（1）*g*cos（sit
（1））））/m
（1）/V
（1）;
Q
（1）=atan（-H
（1）/（30000-x
（1）））+pi;
r
（1）=6708.2039;
R
（1）=-V
（1）*cos（Q
（1））;
n
（1）=Q
（1）+pi;
SIT1
（1）=k/r
（1）*（V
（1）*sin（n
（1）））;
mza=-0.1;%俯仰力矩系数对攻角的偏导数
mzdet=0.024;%俯仰力矩系数对舵偏角的偏导数
t=0;
i=0;
dt=0.01;
ms=0.46;%质量秒消耗量
whileH>0&H1>0%运用迭代法求解
i=i+1;
t=t+dt;
det（i+1）=k1*（sit（i）-sit1（i））+k2*（SIT（i）-SIT1（i））;
a（i+1）=-mzdet/mza*det（i）*180/pi;
Cy（i+1）=0.25*a（i）+0.05*det（i）*180/pi;
Cx（i+1）=0.2+0.005*a（i）^2;
Y（i+1）=Cy（i）*q（i）*S;
X（i+1）=Cx（i）*q（i）*S;
m（i+1）=m（i）-ms*dt;
sit（i+1）=sit（i）+（P*sind（a（i））+（Y（i）-m（i）*g*cos（sit（i））））/m（i）/V（i）*dt;
V（i+1）=V（i）+（P*cosd（a（i））-（X（i）+m（i）*g*sin（sit（i））））/m（i）*dt;
x（i+1）=x（i）+V（i）*cos（sit（i））*dt;
H（i+1）=H（i）+V（i）*sin（sit（i））*dt;
Q（i+1）=atan（-H（i）/（30000-x（i）））+pi;
sit1（i+1）=k*（Q（i）-Q
（1））;
H1（i+1）=H（i）+V（i）*sin（sit1（i））;
SIT（i+1）=（sit（i+1）-sit（i））/dt;
r（i+1）=（H（i）^2+（30000-x（i））^2）^（1/2）;
R（i+1）=（r（i+1）-r（i））/dt;
n（i+1）=acos（-R（i）/V（i））+pi;
SIT1（i+1）=k/r（i）*（V（i）*sin（n（i）））;
T（i+1）=T0-0.0065*H（i+1）;
p（i+1）=p0*（T（i+1）/T0）^4.25588;
q（i+1）=1/2*p（i+1）*V（i+1）^2;
end
plot（x,H）;
holdon
[t,y]=ode45（'jieduan1',[039.0564],[250007000]）;
plot（y（:
3）,y（:
4））;
holdon
[t,y]=ode45（'jieduan2',[39.0564115],[192.768-0.00991002998.71]）;
plot（y（:
3）,y（:
4））;
其中每一段的初始值，均为上阶段的结束值
所以每一阶段计算结束后，需要再给出所有数据的结果，找到每一段距离相对应的数据，即为初始值。
五、结果分析
制出导弹三个阶段的飞行轨迹如图
（1）
图
（1）
图
（2）是第一阶段纵向参数随时间的变化曲线；
图
（2）
图（3）时第二阶段纵向飞行参数随时间的变化曲线
由图
（1）导弹在第一阶段，从初始高度7000m,开始下降飞行，在距离9100m时，开始变为登高飞行，距离达到24000m至目标30000m这一阶段为导弹的下降寻找目标阶段；
由图
（2）得，第二阶段的飞行速度先增加后减小，在第一阶段末尾阶段速度减小至192.768m/s；
弹道倾角先减小后增加，海拔高度随时间的增加而减小；
由图（3）得，第三阶段为登高飞行，所以弹道倾角和海拔高度分别在0度和3050m之间振荡，而速度也基本在140m/s至150m/s之间徘徊；
六、特性点的动力系数、传函
分别取特性点1：
x=0时；
特性点2：
x=9100时；
特性点3：
x=24000时；
特性点4：
x=30000时
由纵向自由扰动的稳定性条件
即纵向自由扰动运动稳定。
根据以下公式：
得到以下值：
特性点1的传递函数：
特征点2的传递函数：
特性点3的传递函数：
特征点4的传递函数：
升级会员 