航天飞行动力学远程火箭弹道设计

航天飞行动力学远程火箭弹道设计

已知火箭纵向运动方程式如公式错误!未找到引用源。所示。

()0sin 1

cos cos sin e e pr P v g m

v P g m x v y v m m m t

A ϕθθαθθ

θαϕθ⎫=+⋅⎪⎪

⎪

=⋅⋅+⋅⎪⎪

⎬=⋅⎪

=⋅⎪

⎪=-⋅⎪

⎪=⋅-⎭

其中，0,,m ,,,,e v P x y θα分别为火箭飞行速度、发动机推力、火箭初始质量、弹道倾角、攻角、水平位移和飞行高度；A ϕ为角度增益系数，t 为火箭飞行时间，m 为火箭质量。仿真初始条件如表1和表2所示。

表1初始状态 序号 变量名 变量值 物理意义及单位 0 t

火箭飞行时间，s 1 θ /2π 初始弹道倾角，弧度 2 v 0 火箭初始速度， /m s

4 x 0 火箭在地面发射坐标系下的初始水平位置,m

5 y

火箭在地面发射坐标系下的初始高度,m 表2 有关参数

序号 变量名 变量值 物理意义及单位 0 0m

8000 起飞质量 kg

1 m 28.57 单位时间燃料质量消耗， /kg s

2 g

9.8 重力加速度常数， 2/N s 3

A ϕ

35 角度增益系数 4 e P

200 发动机推力，KN 5

w

7000

发动机排气速度，/m s

飞行程序角pr ϕ随火箭飞行时间的关系如公式错误!未找到引用源。：

1

21112212123

12302222= =10s, 130, 150 s

60

pr

t t t t t t fig t t t t t t t fig t t t fig t t s t πππϕπ

⎧

≤<⎪⎪⎡⎤⎪⎛⎫

--⎪⎛⎫⎢⎥=+-⋅-⋅

≤<⎨ ⎪ ⎪--⎝⎭⎢⎥⎝⎭⎪⎣⎦⎪

≤≤⎪⎪⎩

==，

问题：（1）请根据如上已知条件，完成火箭纵向运动仿真；

（2）验证齐奥尔科夫斯基公式。

解：

（1）将主动段分为三个阶段：垂直段、转弯段、瞄准段完成火箭纵向运动仿真。 垂直段：10t t <≤，

2

π

ϕ=

pr ，火箭纵向运动方程可进一步化简为：

⎪

⎪⎪⎩⎪⎪

⎪⎨

⎧==⋅+-⋅-⨯⨯=⋅+⋅⋅=+-⨯=⋅+=θθθθπθαθθθsin cos cos 8.9)2()75.288000(3510200cos sin 8.975.28800010200sin 33v y

v x v v t v g v m p t g m p v

e e 初始条件已由题干给出。

function chuizhiduan t0=0;tf=10;

[t,y]=ode45(@diyiduan,[t0 tf],[0 pi/2 0 0]); plot(y(:,3),y(:,4));

%y(1)、y(2)、y(3)、y(4)分别代表速度v,弹道倾角θ，坐标x 及y ，下同。 function dy=diyiduan(t,y) dy=zeros(4,1);

dy(1)=200000/(8000-28.57*t)-9.8*sin(y(2));

dy(2)=7000000*(pi/2-y(2))/(8000-28.57*t)+9.8*cos(y(2)); dy(3)=y(1)*cos(y(2)); dy(4)=y(1)*sin(y(2));

转弯段：21t t t <≤，⎥⎥⎦⎤

⎢⎢⎣

⎡--⎪⎭⎫ ⎝⎛-+=⎥⎥⎦

⎤⎢⎢⎣⎡--⋅-⎪⎪⎭⎫ ⎝⎛--⋅-+=601012010602922)2(221212

121

t t t t t t t t t t fig pr ππππϕ

火箭纵向运动方程可进一步化简为：

⎪

⎪⎪⎩⎪⎪

⎪⎨

⎧==⋅+-⋅-⨯⨯=⋅+⋅⋅=+-⨯=⋅+=θθθθϕθαθθθsin cos cos 8.9)()75.288000(3510200cos sin 8.975.28800010200sin 33

v y

v x v v t v g

v m p t g m p v pr e e 初始条件为垂直段的末端参数：m y m x s m v 2.775,0,2

,/6.156===

=π

θ

function zhuanwanduan t0=10;tf=130;

[t,y]=ode45(@dierduan,[t0 tf],[156.6 pi/2 0 775.2]) plot(y(:,3),y(:,4));

function dy=dierduan(t,y) dy=zeros(4,1);

dy(1)=200000/(8000-28.57*t)-9.8*sin(y(2));

pr=pi/2+(29*pi/60)*((t-10)^2/14400-(t-10)/60);

dy(2)=(7000000/(8000-28.57*t)*(pr-y(2))+9.8*cos(y(2)))/y(1); dy(3)=y(1)*cos(y(2)); dy(4)=y(1)*sin(y(2));

瞄准段：32t t t ≤≤，60

π

ϕ=

=fig pr ，火箭纵向运动方程可进一步化简为：

⎪

⎪⎪⎩⎪⎪

⎪⎨

⎧==⋅+-⋅-⨯⨯=⋅+⋅⋅=+-⨯=⋅+=θθθθπθαθθθsin cos cos 8.9)60()75.288000(3510200cos sin 8.975.28800010200sin 33v y

v x v v t v g

v m p t g m p v e e 初始条件为转弯段的末端参数：km y km x s m v 06.56,177,rad 05949.0,/3702====θ

function miaozhunduan t0=130;tf=150;

[t,y]=ode45(@dierduan,[t0 tf],[3702 0.05949 177000 56060]); plot(y(:,3),y(:,4));

function dy=dierduan(t,y) dy=zeros(4,1);

dy(1)=200000/(8000-28.57*t)-9.8*sin(y(2)); pr=pi/60;

dy(2)=(7000000/(8000-28.57*t)*(pr-y(2))+9.8*cos(y(2)))/y(1); dy(3)=y(1)*cos(y(2)); dy(4)=y(1)*sin(y(2));

在]s 150,0[s t ∈，火箭速度、弹道倾角、在地面发射坐标系下的横纵坐标变化如下：