飞行力学弹道设计与分析.

飞行力学弹道设计与分析
一．题目重述 (1)
二．第一次作业 (2)
1.求解思路 ........................................................................................................................................... 2 2. 计算结果：（在MATLAB 环境下编程求解） (3)
三．第二次作业 (4)
1.要求 (4)
2.计算结果展示及问题分析 (5)
3.精度分析 (7)
四．第三次作业 (8)
1.作业要求 (8)
2.简单思路分析 (8)
3.计算结果及简单问题分析 (8)
4.下面我们讨论步长DT 对结果的影响 (10)
(10)
五．第四次作业 (11)
1.作业要求 (11)
2.简要分析与计算结果 (12)
一．题目重述
某一型号导弹，给定参数如下
0m =320kg; P=2000N; s m =0.46kg/s;
0v =250m/s; 0x =0m; 0H =7000m
0θ=0o ; 0ϕ=0o ; 0α=0o ;
Sref=0.452m ; Lref=2.5m;
0.250.05Cy z αδ=+; 20.20.005Cx α=+;
0.10.024z m z αδ=-+; (气动计算里的角度单位都为度)
y Y C qSref =; x X C qSref =; z Mz m qSrefLref =; 动压 20.5q v ρ=;
其中0 1.2495ρ= 0288.15T =
00.0065T T H =-; 00
=T T ρρ 4.25588（）; 现要求设计,K K 两个参数
在瞬时平衡假设下满足一下飞行方案(飞行中舵偏角015z δ<)
1. 9100x m < 0.0/s m kg s =
*2000cos(0.000314 1.1)5000H x =⨯⨯+
**()()z K H H K H H δ=-+-
2. 910024000m x m <<
*3000H m =(原题3050m 感觉不太合理，这里改为3000m)
*()z K H H KH δ=-+
3. 24000m x < ,0y >
目标位置 30000,0m m x m y m ==
**()z K K θθ
δθθθθ=-+-（） 采用比例导引法攻击目标 二．第一次作业
1.求解思路
选取一个微小的时间步长dt 从初始的飞行状态（可看作初
始点）的0v =250m/s; 0θ=0o ;可以计算出dt 时刻后下一点的x,y.从而可以计算出下一点的H 和H ，带入式*()z K H H KH δ=-+中求出z δ，又根据瞬时平衡假设求出攻角α，再计算出升力，阻力带入导弹纵向平面的运动方程利用欧拉法求解,v θ的微分方程组就可以求出第二个点的,v θ。

然后这样迭代把所有点连线得导弹飞行轨迹。

我们可以用图表示
11111,,,,n n n n n n n v x y H z Yf Xf v θδαθ+++++→→→→→→ 通过此迭代过程我们可以得最后的轨迹
2. 计算结果：（在matlab 环境下编程求解）
K =-0.08，K =-0.032；
其他参数曲线
如下
我们可以看到，计算的结果满足要求
三．第二次作业
1.要求
在前面2段方案飞行的基础上，这次我们加上了第三段的比例导引飞行（攻击目标段）
即要设计
*
b
d dq
K
dt dt
θ
=中的比例系
b
K(这里
b
K为了与前面一
二段的方案飞行的K区别)
和**()z K K θθ
δθθθθ=-+-（）中的K θ，K θ三个参数 2.计算结果展示及问题分析
首先我们选取了==-0.055 =-0.782.5 b K K K θθ，，
（K θ，K θ两个参数取变量对结果影响不大，故这里取常值） 得到如下结果
我们不难发现在舵偏角的变化曲线中在方案飞行与导引飞行衔接处（24000,3000x m y m ==）会发生4度左右的突变，这样的阶跃信号会引起调整时的振动。

是我们不希望看到的 。

而这样的突变靠改变后面的K θ，K θ两个参数（即使使它们随位移或时间变化）是很难实现的。

这里我们经过matlab 环境中反复的实验给出了一种解决方案，我认为这里的根本问题在于方案飞行和导引飞行的舵偏
角控制方程不一样所引发的突变。

所以我们给出了以下的平滑过渡方案
再第三段飞行中给方案飞行的的舵偏角控制方程和导引飞行的舵偏角乘上相应的系数12ηη和，即：
**1()()z K H H K H H δ=-+-（方案飞行的控制方程）
**2()z K K θθ
δθθθθ=-+-（） （导引飞行本来的控制方程） 修正后的导引的舵偏方程为
121+2z z z δηδηδ=⨯⨯
其中sin(2*pi/(4*6000)*X(i 1=1-))+η
sin(2*pi/(4*6000)*X 2=(i))η
12ηη和趋势如图所示
经过以上过度处理后我们可以得到如下结果
我们可以清楚地看到这次弹道倾角的变化趋于平滑。

但在方案飞行中的一段和二段之间舵偏角因发动机突然施加推力且控制方程相同还尚未找到合适的办法使之平滑过渡。

3.精度分析
对攻击目标时做精度分析
取弹道中后面的几个点分析
X=[2.9996 2.9999 3.0001 3.0003]*4
10m Y=[ 1.7486 0.6614 -0.4246 -1.5096]m
取(X(2)，Y(2)),(X(3),Y(3))两个点做线化分析可以得到距目标点（30000，0）的误差为0.12m。

四．第三次作业
1.作业要求
在原有的基础上去掉瞬时平衡假设，考虑导弹的动态特性。

进一步得到新的弹道，与原有的理想弹道做对比。

2.简单思路分析
在去掉瞬时平衡假设后，攻角的计算由原来的通过力矩平衡由舵偏角计算z δα→变为 =-αϕθ
其中θ的计算和原来的一样，而ϕ可由z ω得出
3.计算结果及简单问题分析
经过稍微的调整系数后K =-0.03，K =-0.032
dt=0.01s
可以看到考虑了动态特性后，导弹在比例导引段与原有的理想弹道相差较大而方案飞行段与理想弹道吻合度很高。

分析其他参数
我们可以看到在比例导引阶段有三次明显的震荡。

4.下面我们讨论步长dt对结果的影响
当dt>0.024s时弹道发散不能形成完整弹道图。

当dt=0.0224s
dt=0.0215s
dt=0.0212
我们可以看到当dt 取得越小时弹道拟合的越好（dt<0.024s ），这也一定程度上说明了欧拉法的选择收敛性，而非无条件收敛，所以最终我们选取dt=0.01s （在与瞬时平衡弹道对比图中dt=0.01s ）就能很好的拟合。

五．第四次作业 1.作业要求
在前面的基础上，求出5个动力系数，选取合适的特征点，并分别求出舵偏角对弹道倾角，及高度的传递函数（分别几位,H G G δθδ）。

2.简要分析与计算结果
1）特征点选取
我们共选取了4个特征点，主要的依据是画出5个动
力系数的图后，在图中选择一些拐点，驻点，最值点。

他们
所对应的为特征点。

下图为5个动力系数的图
我们选取如下特征点
特征点X a1 a2 a3 a4 a5 注
1 9101.2m 109.40 5.47 -1.31 0.064 0.0059 发动机开始工作
2 30000m 298.34 14.92 -3.58 0.095 0.013
3 攻击目标时
3 23999m 156.32 7.91 -1.88 0.072 0.0079 转为导引飞行时
4
6349.7m 150.90 7.55 -1.81 0.033 0.0067 第一段飞行拐点
（a1,a2,a3,a4,a5的单位在上图中有说明） 2）传递函数的求解
在忽略长周期干扰的影响下我们利用公式
25153425
32
123
()=-=()a s a a s a a a a s G z s s A s A s A δθθδ-+-∆∆+++ （1） 其中 114A a a =+;2214A a a a =+;30A ≈ 可以求出对θ的传函 再利用公式
00(sin )(cos )d y
V V dt
θθθ∆=∆+∆ ( 2 ) 式中0(sin )0V θ∆≈两遍同取拉式变换把（1）带入可得
25153425
32
123
cos =H a s a a s a a a a V G s s A s A s A θθ-+-+++ ( 3) 由上面的公式（1），（3）可以得到传函,H G G δθδ。

下面分别给出4个特征点的传函 特征点1
2132
0.00590.64830.1165
=109.512.47s s G s s s
δθ--++ 2132
181.7(0.00590.64830.1165)
=(109.512.47)
H s s G s s s s θ⨯--++ 特征点2
2232
0.01328 3.9620.5371
=298.443.16s s G s s s
δθ--++
2232
234.7(0.01328 3.9620.5371)
=(298.443.16)
H s s G s s s s θ⨯--++ 特征点3
2332
0.007894 1.2340.1959
=156.419.0s s G s s s
δθ--++ 2H 332
217.2(0.007894 1.2340.1959)
=(156.419.0)s s G s s s s θ⨯--++ 特征点4
2432
0.0067 1.0060.1106
=150.912.57s s G s s s
δθ--++ 2432
194.3(0.0067 1.0060.1106)
=(150.912.57)
H s s G s s s s θ⨯--++ (以上的传函与上表中的特征点一一对应)。