弹道计算大作业Word文件下载.docx
《弹道计算大作业Word文件下载.docx》由会员分享,可在线阅读,更多相关《弹道计算大作业Word文件下载.docx(17页珍藏版)》请在冰豆网上搜索。
通过模型计算可得
m/s2
零升阻力系数
0.02
/rad
密度
kg/m3
1.2仿真要求
请使用Simulink或Buildfly完成以下仿真任务:
(1)请完成该导弹的无控飞行弹道仿真;
(2)请完成该导弹的平衡滑翔方案飞行弹道仿真;
(3)请完成该导弹的最大升阻比滑翔飞行弹道仿真;
二、模型的建立
2.1升力和阻力模型
已知展弦比为
的飞行器的升力线斜率为:
根据飞行力学相关知识,飞行器的升力系数和阻力系数为:
其中,升力线斜率由
(1)式可得;
为效率系数:
。
由升力系数和阻力系数,得到导弹的升力和阻力为:
2.2大气和重力加速度模型
在计算过程中,大气密度采用如下模型:
其中,
为海平面的大气密度;
重力加速度采用如下模型:
,
为地球半径;
为飞行器距离地面的高度。
2.3无控飞行
假设导弹的运动始终在铅垂平面,根据飞行力学知识,得到导弹无控飞行时的运动学和动力学方程为:
在上述模型中,假设俯仰角
为0。
2.4平衡滑翔
所谓的“平衡”可以理解为垂直于速度方向受力平衡,即
因此得到平衡滑翔时的导弹运动学和动力学方程:
由于弹道倾角的变化率为常数,方程组中的第二个方程等于0。
这个方程可以用来求攻角
2.5最大升阻比滑翔飞行弹道
联立
(1)式、
(2)式可得升阻比的表达式为:
从上式可以看出,由于展弦比
、零升阻力系数
为常数,因此升阻比只和攻角有关,是关于攻角的函数。
因此要使升阻比达到最大,须使
得到
因此,以最大升阻比滑翔时导弹运动学和动力学方程为:
三、仿真结果
3.1无控飞行弹道仿真
根据无控弹道模型,写出s函数,搭建的仿真模块如下图所示:
图1无控飞行仿真模块
由于初始条件给定,因此模块没有输入;
输出有六个,分别为导弹的射程变化、高度变化、速度变化、弹道倾角变化、攻角变化以及密度变化。
模块的仿真时间由高度变化决定,当高度降为0(导弹落到地面上)时仿真结束。
导出数据后画图如下:
图2无控飞行时各参数变化
3.2平衡滑翔弹道仿真
平衡滑翔弹道仿真模块如下图所示:
图3平衡滑翔模块
取仿真时间为150s,无输入,输出分别为:
导弹的射程变化、高度变化、速度变化、弹道倾角变化、攻角变化以及密度变化。
得到各参量时间变化图如下:
图4平衡滑翔飞行时各参数变化
3.3最大升阻比滑翔弹道仿真
按最大升阻比飞行时弹道仿真模块如下图所示:
图5最大升阻比飞行模块
取仿真时间为180s,无输入,输出分别为:
图4最大升阻比飞行时各参数变化
附录
附表1无控弹道飞行时完整的s函数
无控弹道
function[sys,x0,str,ts,simStateCompliance]=trace2(t,x,u,flag)
switchflag,
case0,
[sys,x0,str,ts,simStateCompliance]=mdlInitializeSizes;
case1,
sys=mdlDerivatives(t,x,u);
case2,
sys=mdlUpdate(t,x,u);
case3,
sys=mdlOutputs(t,x,u);
case4,
sys=mdlGetTimeOfNextVarHit(t,x,u);
case9,
sys=mdlTerminate(t,x,u);
otherwise
DAStudio.error('
Simulink:
blocks:
unhandledFlag'
num2str(flag));
end
function[sys,x0,str,ts,simStateCompliance]=mdlInitializeSizes
sizes=simsizes;
sizes.NumContStates=4;
sizes.NumDiscStates=0;
sizes.NumOutputs=5;
sizes.NumInputs=0;
sizes.DirFeedthrough=0;
sizes.NumSampleTimes=1;
sys=simsizes(sizes);
x0=[0;
2000;
100;
-5/180*pi];
str=[];
ts=[00];
simStateCompliance='
UnknownSimState'
;
functionsys=mdlDerivatives(t,x,u)
S=1.7;
%参考面积,m^2
AR=0.86;
%展弦比
e=0.9;
%效率因子;
m=115;
%质量,kg
g0=9.8;
%海平面重力加速度,m/s^2
Rd=6371000;
%地球半径
r=Rd/(Rd+x
(2));
g=g0*r^2;
%飞行器所在高度的重力加速度
rho0=1.225;
%海平面大气密度,kg/m^3
T0=288.15;
rho=rho0*(1-0.0065*x
(2)/T0)^4.2288;
%飞行器所在高度的大气密度
alpha=-x(4);
%无控飞行时
CLa=3.141592*AR/(1+sqrt(1+(AR/2)^2));
%升力线斜率,/rad
CDo=0.02;
%零升阻力系数
epsilon=1/(pi*e*AR);
%诱导阻力因子
CL=CLa*alpha;
%升力系数
CD=CDo+epsilon*CL^2;
%阻力系数
X=CD*1/2*rho*x(3)^2*S;
Y=CL*1/2*rho*x(3)^2*S;
%以下为飞行器在铅垂平面的运动方程
dx=x(3)*cos(x(4));
dy=x(3)*sin(x(4));
dv=-X/m-g*sin(x(4));
dtheta=Y/(m*x(3))-g*cos(x(4))/x(3);
sys=[dx;
dy;
dv;
dtheta];
functionsys=mdlUpdate(t,x,u)
sys=[];
functionsys=mdlOutputs(t,x,u)
y1=x
(1);
y2=x
(2);
y3=x(3);
sys=[x
(1)x
(2)x(3)x(4)rho];
functionsys=mdlGetTimeOfNextVarHit(t,x,u)
sampleTime=1;
%Example,setthenexthittobeonesecondlater.
sys=t+sampleTime;
functionsys=mdlTerminate(t,x,u)
附表2平衡滑翔飞行部分代码
平衡滑翔飞行
alpha=2*m*g*cos(x(4))/(rho*x(3)^2*S*CLa);
X=CD*1/2*rho*sqrt(x(3)^2)*S;
Y=CL*1/2*rho*sqrt(x(3)^2)*S;
dtheta=0;
dy;
dv;
rho0=1.225;
T0=288.15;
rho=rho0*(1-0.0065*x
(2)/T0)^4.2288;
S=1.7;
AR=0.86;
m=115;
g0=9.8;
Rd=6371000;
r=Rd/(Rd+x
(2));
g=g0*r^2;
%飞行器所在高度的重力加速度
CLa=pi*AR/(1+sqrt(1+(AR/2)^2));
alpha=2*m*g*cos(x(4))/(rho*x(3)^2*S*CLa);
y
(1)=x
(1);
y
(2)=x
(2);
y(3)=x(3);
y(4)=x(4);
y(5)=alpha;
y(6)=rho;
sys=[y];
附表3最大升阻比飞行部分代码
最大升阻比飞行
alpha=sqrt(CDo*pi*AR*e)/CLa;
%零升阻力系数
sys=[x
(1)x
(2)x(3)x(4)alpharho];