弹道计算大作业之欧阳道创编.docx
《弹道计算大作业之欧阳道创编.docx》由会员分享，可在线阅读，更多相关《弹道计算大作业之欧阳道创编.docx（12页珍藏版）》请在冰点文库上搜索。
弹道计算大作业之欧阳道创编
弹道计算大作业
时间：
2021.03.06
创作：
欧阳道
一、初始条件和要求
1.1初始条件
已知给定的初始条件如下：
表1初始条件
名称
符号
给定值
单位
参考面积
1.7
m2
展弦比
0.86
/
效率因子
0.9
/
质量
115
kg
重力加速度
通过模型计算可得
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）;
rho0=1.225;
T0=288.15;
rho=rho0*（1-0.0065*x
（2）/T0）^4.2288;
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）
sys=[];
附表2平衡滑翔飞行部分代码
平衡滑翔飞行
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;%飞行器所在高度的大气密度
CLa=3.141592*AR/（1+sqrt（1+（AR/2）^2））;%升力线斜率,/rad
CDo=0.02;%零升阻力系数
epsilon=1/（pi*e*AR）;%诱导阻力因子
alpha=2*m*g*cos（x（4））/（rho*x（3）^2*S*CLa）;
CL=CLa*alpha;%升力系数
CD=CDo+epsilon*CL^2;%阻力系数
X=CD*1/2*rho*sqrt（x（3）^2）*S;
Y=CL*1/2*rho*sqrt（x（3）^2）*S;
dx=x（3）*cos（x（4））;
dy=x（3）*sin（x（4））;
dv=-X/m-g*sin（x（4））;
dtheta=0;
sys=[dx;dy;dv;dtheta];
functionsys=mdlOutputs（t,x,u）
rho0=1.225;
T0=288.15;
rho=rho0*（1-0.0065*x
（2）/T0）^4.2288;
S=1.7;%参考面积,m^2
AR=0.86;%展弦比
m=115;%质量,kg
g0=9.8;%海平面重力加速度,m/s^2
Rd=6371000;%地球半径
r=Rd/（Rd+x
（2））;
g=g0*r^2;%飞行器所在高度的重力加速度
CLa=pi*AR/（1+sqrt（1+（AR/2）^2））;%升力线斜率,/rad
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最大升阻比飞行部分代码
最大升阻比飞行
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;%飞行器所在高度的大气密度
CLa=3.141592*AR/（1+sqrt（1+（AR/2）^2））;%升力线斜率,/rad
CDo=0.02;%零升阻力系数
epsilon=1/（pi*e*AR）;%诱导阻力因子
alpha=sqrt（CDo*pi*AR*e）/CLa;
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=mdlOutputs（t,x,u）
rho0=1.225;
T0=288.15;
rho=rho0*（1-0.0065*x
（2）/T0）^4.2288;
AR=0.86;%展弦比
e=0.9;%效率因子;
CLa=3.141592*AR/（1+sqrt（1+（AR/2）^2））;%升力线斜率,/rad
CDo=0.02;%零升阻力系数
alpha=sqrt（CDo*pi*AR*e）/CLa;
sys=[x
（1）x
（2）x（3）x（4）alpharho];
时间：
2021.03.06
创作：
欧阳道
升级会员 