航天飛行動力學大作業(yè)

韓謹陽

2015300464

1、方案飛行2、彈道設(shè)計3、衛(wèi)星攝動與機動

第一部分

飛行方案

衛(wèi)星的攝動與機動

彈道設(shè)計

飛行方案大作業(yè)

一、問題描述

在已知導彈質(zhì)量、轉(zhuǎn)動慣量、發(fā)動機推力等參數(shù)的狀況下，導彈分為三個飛行方案，即三個

階段飛行。

階段一：

飛行距離在xv9IOO，〃，采用追蹤法，其中方案高度與距離的關(guān)系、方案彈道傾角與高度

的關(guān)系如下：

H*=2000xcos(0.000314x1.1xx)+5000

階段二：

飛行距圖在24000/〃>x>9100，〃，采用途標法，其中方案舟度與距崗的關(guān)系、方案異道

傾角與高度的關(guān)系、導彈因燃料消托而質(zhì)量變化參數(shù)如下：

H'=3050m

...(2)

&)+勺〃

mt=0.46^/s(3)

階段三：

飛行方案x>24000/7：&&y>0,而豉終目的位JE為為”=3()()()()〃？

采用比例導引法

r—=Vxsin7-V7sinrj

（ItMr

MF

dO,,dq

-----=k-L（4）

dtdt

。'-On=k（q-q。）

S.=%（。-〃）+（（6-?！?/p>

規(guī)定：

1）計算縱向理想彈道，給出采用瞬時平衡假設(shè)=0時所有縱6參數(shù)隨時間的

變化曲線。

2）不考慮氣動力下沈影響，計算飛行器沿理想彈道飛行時，你認為可以作為特性點的5個

以上點處的縱向短周期擾動運動的動力系數(shù)，并分析其在特性點處的自由擾動的穩(wěn)定

性，以及計算在各個特性點處彈體傳遞函數(shù)I曝（S），W,6（S）,1囁（S）。

二、建立模型

基于“瞬時平衡”假設(shè)，導彈在船垂平面內(nèi)運動的質(zhì)心運動方程組為：

dV.

m—=Pcosa,-Xb-mgsin0

mV-=Ps\x\a+Y-mgcos0

戊hh⑸

—=Vcos0

—=rsin0

由于階段一不考慮導彈質(zhì)量應時間的變化，因此階段一的模型需要聯(lián)立公式（1）、公式（5）：

其中攻角a可根據(jù)瞬時平衡假設(shè)

nfa+=0

從而可得到導彈攻角與彈道傾角之間的關(guān)系

a-(6)

2c#/

其中

丫"Cg

其中假設(shè)公式(1)的用=A〃(夕夕)此(夕一。)中的勺=9勺=05：

又由于階段二需要考慮導彈質(zhì)量防葉間的變化，因此階段二的模型需安聯(lián)立公式(2)公式

(5)、公式(6)、公式(7)

最終一階段，由于運用了比例導引法

公式(4)的k=2,可得導彈抵達目的的相對微分方程為

dr

—=-Kcosrj

r—=Vsinrj

dt

而導引率

(10,dq

——=k-L、其中k=2;

dtdt

由于第三階段的初始參數(shù)及終點坐標均為直角坐標系，由下圖可知將

x=30000-rcosq

y=rsinq

代入到公式(4),得到直角坐標系下的微分方榴組

內(nèi)

打

心?

-s4s1n。

五

五

瓦co

附

力

力

--

加

五

-加s1ngs1ng

此外撲充方程法向平衡方程:

mV—=PsinaY-mgcos0

dt

三、算法實現(xiàn)

編杼使用MATLAB軟件，并運用歐拉方代解微分方侵，印。de45畫數(shù)：

四、程序源代碼

*************************階段-***********************

fundiondy=jieduanl(t,y)

dy=zeros(4,1)；

m=320：

g=9.8；

P=;

q=0.5*1.2495*((288.15-0.0065*y⑷)/288,15)<4.2558*y(1)<2;

k=-9；

dk=-0.5；

Hi=*cos(0.000314*1.1*y(3))+5000:

dHi=-*0.000314*1.1*sin(y⑶)；

delta=k*(y(4)-Hi)+dk*(dy(3)-dHi);

alpha=0.34*deIta;

Xb=(0.2+0.005*alpha"2)*q*0.45；

Yb=(0.25*alpha+0.05*deIta)*q*0.45；

dy=zeros(4.1);

dy(1)=P*cos(aIpha)/m-Xb/m-g*sin(y(2))；

dy(2)=P*sin(aIpha)/m/y(1)+Yb/1n/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=；

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*deIta;

Xb=(0.2+0.005*aIpha'2)*q*0.45：

Yb=(0.25*alpha+0.05*deIta)*q*0.45；

dy(1)=P*cos(alpha/180*pi)/nr-Xb/m-?*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

*******************************汾段三*at：******************************

functiondy=jieduan3(t,y)

v=y(4)；

k=10;

m=285.04-0.46*t；

qO=-atan(3050/6000)；

g=9.8;

q1=0.5*1.2495*((288.15-0.0065*y(2))/288.15)44.2558*y(4)/2;

k1=10;

dk1=0.05;

dy=zeros(4,1)；

r=sqrt(y⑴-2+y⑵-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；

deIta=k1*(y(3)-tht)+dk1*(dy(3)-dtnt);

alpha=0.34*deIta;

dy(1)=-dr*cos(q)+r*sin(q)*dq；

dy(2)=-dr*sin(q)-r*cos(q)*dq;

Yb=(0.25*alpha+0.05*deIta)*q1*0.45;

dy(3)=(*sin(aIpha)/m+Yb/m-g+cos(y(3)))/v；

v(4)=v;

end

***********************************main函數(shù)

************************************

m(1)=287.2204：%導彈質(zhì)量

P=;$發(fā)動機拉力

g=9.8；

k=5;

det(1)=0.045;

a(1)=0.6186；

sit(1)=-0.;

V(1)=217.2867；%初始速度

x(1)=24000；先初始性五

H(1)=3071;%初始高度

H1(1)=3050;

S=0.45；%參照面積

L=2.5：氣賽照長度

k1=-0.14；

k2=-O.06；

sit1(1)=sit(1);

p0=1.2495;

T0=288.15;

T(1)=T0-0.0065*H(1);

p(1)=pO*(T(1)/TO)-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;%升力系數(shù)

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⑴))+pi;

r(1)=6708.2039;

R(1)=-V(1)*cos(0(1));

n(1)=Q(1)+pi;

SIT1(1)=k/r(1)*(V(1)*sin(n(1)))；

mza=-0.1;%俯仰力矩系數(shù)對攻角的偏導數(shù)

mzdet=0.024；%俯仰力矩系數(shù)對舵偏角的偏導數(shù)

t=0;

i=0;

dt=0.01;

ms=0.46；$質(zhì)量秒消耗量

whileH>0&H1>0%運用迭代法求解

i=i+1;

t=t+dt;

det(i+1)=k1*(sit(i)-sit1(i))+<2*(SIT(i)-SIT1(i));

a(i+1)=-mzdet/mza*det(i)*180/oi;

Cy(i+1)=0.25*a(i)+0.05*det⑴*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))+1(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*(0(i)-Q(1));

HI(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)A(1/2)；

R(i+1)=(r(i+1)-r(i))/dt;

n(i+1)=acos(-R⑴/V⑴)+pi;

SIT1(i+1)=k/r(i)*(V(i)*sin(n(i)))；

T(i+1)=T0-0.0065*H(i+1)；

p(i+1)=pO*(T(i+1)/TO)N.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))；

其中每一段的初始值.均為上階段的結(jié)束值

因此每一階段計算結(jié)束后,需要再給出所有效據(jù)的成果，找到每一段距離相對應的效據(jù)，即

為初始值。

五、成果分析

制出導彈三個階段的飛行軌跡如圖(1)

圖（1）

圖（2）是第一階段縱向參數(shù)隨時間的變化曲線;

圖（3）時第二階段縱向飛行參數(shù)隨時間的變化曲線

工Hv銳使化曲線；Jut,T”曲變ii1曲線□

由圖(1)導彈在第一階段，從初始高度7000m,開始下降飛行，在距離9100m時，

開始變?yōu)榈菎u飛行，龍離到達24000m至目的30000m這一階段為導彈的下降尋找目的階

段；

由圖(2)得.第二階段的飛行速度先增長后減小.在第一階段末尾階段速度減小至

192.768m/s：

彈道傾角先減小后增長，海拔高度版時間的增長而減?。?/p>

由圖(3)得，第三階段為登高飛行，因此彈道傾角和海拔高度分別在0度和3050m之

間振蕩，而速度也基本在140m/s至150m/s之間徘徊；

六、特性點的動力系數(shù)、傳函

分別取特性點1:x=0時；

特性點2：x=9100時：

特性點3:x=24000時:

特性點4:x=30000時

由縱向自由擾動的穩(wěn)定性條件％I+>0即縱向自由擾動運動椅定。

根據(jù)如下公式:

K=45%—%%

%+%%

61

Ml+丹2%

%+%+%.

+%%

(嵋)。=(m^5￡)0

=（m%S&

（M泣

%

（%）。

J,。

（V）。

%

J,。

（0+尸’）

（叫

（3

（叫

得到如下值:

M：。M6

?nPnVaNO

構(gòu)性點170000.60122502113.912507.33

特性點291000.9273192-3676.815882.4357

將件點3210000.92311206-1971.115?173.86694

特性點4300001.2195195-2672.563641.4152

r

。22?24?25-4?35Kalfl

特性點10.006349216.7108-1.51010.051420.0052818-0.01760.3860.0115

特性點20.00631911.672-2.30140.070960.009191-0.026220.29270.01131

特性點30.0063496.26808-1.534340.049680.004936-0.016680.399410.011189

特性點40.006349218.18431-20360.0735810.0084-0.026080.31330.01372

特性點1的傳遞函數(shù):

v(s)=-0.01