铅垂平面飞行弹道仿真及分析.docx

铅垂平面飞行弹道仿真及分析.docx

《铅垂平面飞行弹道仿真及分析.docx》由会员分享，可在线阅读，更多相关《铅垂平面飞行弹道仿真及分析.docx（18页珍藏版）》请在冰豆网上搜索。

铅垂平面飞行弹道仿真及分析

数值实验作业

1实验背景

实验名称：

铅垂平面飞行弹道仿真及分析

实验内容与要求：

根据描述飞行器在铅垂平面内运动的数学模型，编制某导弹的铅垂平面无控飞行弹道仿真程序，利用计算机解算初始段无控飞行弹道，对初始段弹道参数的变化规律进行分析。

2建立数学模型：

3计算方法研究

确定数值积分方法和积分步长

使用算法：

四阶龙格—库塔法。

使用步长：

h=0.005

4原始数据：

1）初值

x=0（m）y=20.0（m）=18=18v=20（m/s）z=0（rad/s）m=52.38（kg）

2）攻角与马赫数范围（仅用于插值计算）

攻角=0~10马赫数=0.1~0.9

3）阻力系数表

马赫数

攻角（）

0

2

4

6

8

10

0.1

.4177

.4404

.5219

.6603

.8534

1.1023

0.2

.3858

.4086

.4903

.6290

.8226

1.0723

0.3

.3779

.4007

.4827

.6218

.8160

1.0666

0.4

.3785

.4015

.4838

.6234

.8184

1.0700

0.5

.3787

.4018

.4846

.6249

.8209

1.0738

0.6

.3829

.4062

.4897

.6310

.8284

1.0835

0.7

.3855

.4091

.4934

.6363

.8358

1.0938

0.8

.4082

.4321

.5175

.6621

.8641

1.1254

0.9

.4947

.5192

.6073

.7571

.9672

1.2392

4）升力系数表

马赫数

攻角（）

0

2

4

6

8

10

0.1

.0000

.6430

1.4758

2.2870

3.0713

3.8463

0.2

.0000

.6454

1.4807

2.2942

3.0814

3.8598

0.3

.0000

.6480

1.4858

2.3014

3.0915

3.8731

0.4

.0000

.6512

1.4923

2.3107

3.1039

3.8891

0.5

.0000

.6554

1.5007

2.3227

3.1197

3.9092

0.6

.0000

.6617

1.5134

2.3409

3.1436

3.9401

0.7

.0000

.6698

1.5304

2.3661

3.1775

3.9835

0.8

.0000

.6792

1.5501

2.3950

3.2162

4.0323

0.9

.0000

.6933

1.5935

2.4706

3.3273

4.1790

5）推力数据

t（s）

.000

.15

.49

2.11

2.27

3.53

8.78

25.45

42.80

43.68

44.08

P（kgf）

331.2

614.3

505.4

607.8

48.65

43.97

42.01

41.00

40.80

40.79

2.22

第一级工作结束时间：

2.1126s，第二级工作结束时间：

44.0832s

6）发动机质量秒流量

t（s）

0.

2.1

2.105

44.1

44.105

100

秒流量（kg/s）

2.362

2.362

0.21059

0.21059

0.

0.

7）转动惯量

t（s）

.0

2.0

2.4

6.4

10.4

14.4

18.4

22.4

26.4

30.4

34.0

38.4

42.4

44.0

Jz（kgms）

8.35

7.88

7.86

7.81

7.78

7.75

7.73

7.71

7.70

7.70

7.69

7.69

7.69

7.69

8）导弹重心（起自头部）

t（s）

.0

2.0

2.4

10.0

18.0

26.0

32.0

38.0

42.0

44.0

XG（m）

.9381

.9095

.9091

.9026

.8969

.8928

.8907

.8896

.8895

.8896

9）静稳定力矩系数

马赫数

攻角（）

0

2

4

6

8

10

0.1

0.0000

-0.0104

-0.0341

-0.0564

-0.0771

-0.0985

0.2

0.0000

-0.0104

-0.0341

-0.0564

-0.0770

-0.0983

0.3

0.0000

-0.0104

-0.0341

-0.0564

-0.0769

-0.0982

0.4

0.0000

-0.0105

-0.0342

-0.0564

-0.0768

-0.0979

0.5

0.0000

-0.0104

-0.0339

-0.0560

-0.0761

-0.0969

0.6

0.0000

-0.0093

-0.0314

-0.0521

-0.0708

-0.0903

0.7

0.0000

-0.0080

-0.0286

-0.0477

-0.0650

-0.0829

0.8

0.0000

-0.0065

-0.0252

-0.0425

-0.0578

-0.0739

0.9

0.0000

-0.0053

-0.0229

-0.0391

-0.0538

-0.0693

当导弹重心变化时的修正公式：

10）阻尼力矩导数

当Xg=.9381时

马赫数

攻角（）

0

2

4

6

8

10

0.1

-0.4686

-0.4829

-0.4982

-0.5130

-0.5272

-0.5409

0.2

-0.4707

-0.4850

-0.5003

-0.5150

-0.5292

-0.5429

0.3

-0.4744

-0.4886

-0.5039

-0.5186

-0.5327

-0.5464

0.4

-0.4797

-0.4939

-0.5090

-0.5237

-0.5378

-0.5514

0.5

-0.4882

-0.5022

-0.5173

-0.5318

-0.5458

-0.5593

0.6

-0.5089

-0.5227

-0.5376

-0.5520

-0.5658

-0.5791

0.7

-0.5366

-0.5502

-0.5649

-0.5790

-0.5927

-0.6058

0.8

-0.5738

-0.5871

-0.6014

-0.6153

-0.6287

-0.6415

0.9

-0.6272

-0.6407

-0.6553

-0.6694

-0.6830

-0.6960

当Xg=.8896时

马赫数

攻角（）

0

2

4

6

8

10

0.1

-0.6179

-0.6384

-0.6600

-0.6805

-0.6999

-0.7182

0.2

-0.6207

-0.6410

-0.6626

-0.6830

-0.7024

-0.7207

0.3

-0.6253

-0.6455

-0.6670

-0.6874

-0.7067

-0.7249

0.4

-0.6319

-0.6521

-0.6734

-0.6937

-0.7129

-0.7310

0.5

-0.6424

-0.6624

-0.6835

-0.7036

-0.7226

-0.7406

0.6

-0.6669

-0.6866

-0.7074

-0.7272

-0.7459

-0.7636

0.7

-0.6997

-0.7190

-0.7395

-0.7589

-0.7774

-0.7948

0.8

-0.7435

-0.7624

-0.7824

-0.8014

-0.8194

-0.8365

0.9

-0.8069

-0.8266

-0.8474

-0.8672

-0.8859

-0.9035

11）其它参数

特征面积S（m2）

特征长度L（m）

毛翼展（m）

音速SONIC（m/s）

大气密度

（kg/m3）

0.0227

1.8

0.5

343.13

1.225

使用的插值算法：

气动数据插值——等距双变元抛物线插值；

推力、重心、转动惯量等——不等距一元线性插值。

5空气动力和空气动力矩表达式

6编制计算程序

计算机算法采用了matlab实现,源程序见附件.

该程序有八个函数组成,各函数之间的调用关系如下图所示.

1）子函数initl的功能是输入求解导弹运动方程组所需的原始数据.

2）子函数rk_4是四阶龙格-库塔法积分算法子函数,其中调用了子函数dery.

3）子函数dery功能是计算微分方程组的右端函数,其中调用了子函数interp.

4）插值子函数interp功能是计算所有需要插值的参数,其中调用了子函数interp11和interp33.

5）子函数interp11是不等距单变元线性插值函数,主要用于转动惯量Jz,质心位置xg,推力p等单变量的插值.

6）子函数interp33是等距双变元抛物线-线性插值函数,主要用于气动力系数cx,cy,气动力矩系数

等双变元参数的插值,该函数调用了子函数interp31.

7）子函数interp31是等距单变元抛物线插值函数,被interp33调用,完成子函数interp33的等距双变元抛物线-线性插值功能.

关于程序中出现的数组和变量名,作如下说明:

acx:

阻力系数数组;

acy:

升力系数数组;

amzaf:

某一质心位置下的静稳定性导数数组;

amzwz:

阻尼力矩系数导数数组;

axg:

质心位置Xg随时间的变化规律;

ajz:

转动惯量Jz随时间的变化规律;

ap:

起飞,续航发动机的推力值;

amc:

起飞,续航发动机的燃料质量的秒流量值;

agc:

质心位置变化的始末值;

andm:

气动数据插值所需的Ma数的最小,最大值;

andaf:

气动数据插值所需的攻角

的最小,最大值;

y:

存放积分结果的数组,该数组在程序开始时存放积分初值;

dy:

存放右端函数数值的数组;

b:

存放三个时间值的数组,其中b

（2）存放起飞发动机工作结束时间;b（3）存放续航发动机工作结束时间;

程序中一些主要的变量名有:

L—特征长度;

S—特征面积;

SONIC—声速C;

RHO—大气密度

;

h—积分步长;

程序中其他变量都是存放中间结果的变量.

7计算结果

运行程序得到弹道曲线,速度曲线,攻角曲线如下:

通过运用龙格库塔法和一些插值方法求解导弹运动方程组,获得了导弹各运动参数的变化规律.通过该算法获得的变化规律比较接近真实情况.

附件:

源程序

functionmain（）

clearall

clearglobal

globaly;

globalii;

ii=0;

h=0.005;

initl（）;

whiley（7）>=0

ii=ii+1;

result（ii）;

rk_4（8,h）;

end

savedata（ii）;

drawing（）;

%原始数据初始化

functioninitl（）

globalacxacyajzamzafamzwzaxgapamcagcandmandafbLSSONICRHO;

globaly;

y=[020.1801802052.38];

%马赫数ma

%攻角alpha

%三个时间b

（1）为导弹离轨时间b

（2）为起飞发动机工作结束时间b（3）为续航发动机工作结束时间

b=[02.112644.0832];

%系数表维数

n1=9;

n2=6;

%andm最小,最大值

andm

（1）=0.1;

andm

（2）=0.9;

%andaf最小,最大值

andaf

（1）=0;

andaf

（2）=10;

%阻力系数

acx=[.4177.4404.5219.6603.85341.1023;

.3858.4086.4903.6290.82261.0723;

.3779.4007.4827.6218.81601.0666;

.3785.4015.4838.6234.81841.07;

.3787.4018.4846.6249.82091.0738;

.3829.4062.4897.6310.82481.0835;

.3855.4091.4934.6363.83581.0938;

.4082.4321.5175.6621.86411.1254;

.4947.5192.6073.7571.96721.2392];

%升力系数

acy=[.0000.64301.47582.28703.07133.8463;

.0000.64541.48072.29423.09153.8731;

.0000.64801.48582.30143.09153.8731;

.0000.65121.49232.31073.10393.8891;

.0000.65541.50072.32273.11973.9092;

.0000.66171.51342.34093.14263.9401;

.0000.66981.53042.36613.17753.9835;

.0000.67921.55012.39503.21624.0323;

.0000.69331.59352.47063.32734.1790];

%推力

ap=[.000.15.492.112.273.538.7825.4542.8043.6844.08;

331.2614.3505.4607.848.6543.9742.0141.0040.8040.792.22];

%发动机质量秒流量

amc=[0.2.12.10544.144.105100;

2.3622.3620.210590.210590.00000.0000];

%转动惯量

ajz=[0.2.02.46.410.414.418.422.426.430.434.038.442.444.0;

8.357.887.867.817.787.757.737.717.707.707.697.697.697.69];

%导弹重心

axg=[.02.02.410.018.026.032.038.042.044.0;

.9381.9095.9091.9026.8969.892808907.8896.8895.8896];

agc=[.9381.8896];

%静稳定力矩系数

amzaf=[0.0000-0.0104-0.0341-0.0564-0.0771-0.0985;

0.0000-0.0104-0.0341-0.0564-0.0770-0.0983

0.0000-0.0104-0.0341-0.0564-0.0769-0.0982;

0.0000-0.0105-0.0342-0.0564-0.0768-0.0979;

0.0000-0.0104-0.0339-0.0560-0.0761-0.0969;

0.0000-0.0093-0.0314-0.0521-0.0708-0.0903;

0.0000-0.0080-0.0286-0.0477-0.0650-0.0829;

0.0000-0.0065-0.0252-0.0425-0.0578-0.0739;

0.0000-0.0053-0.0229-0.0391-0.0538-0.0693];

%阻尼力矩导数

%当xg=.9381时

amzwz=[-0.4686-0.4829-0.4928-0.5130-0.5272-0.5409;

-0.4707-0.4850-0.5003-0.5150-0.5292-0.5429;

-0.4744-0.4886-0.5039-0.5186-0.5327-0.5464;

-0.4797-0.4939-0.5090-0.5237-0.5378-0.5514;

-0.4882-0.5022-0.5173-0.5318-0.5458-0.5593;

-0.5089-0.5227-0.5376-0.5520-0.5658-0.5791;

-0.5366-0.5502-0.5649-0.5790-0.5927-0.6058;

-0.5738-0.5871-0.6014-0.6153-0.6287-0.6415;

-0.6272-0.6407-0.6553-0.6694-0.6830-0.6960];

%%当xg=.8896时

%其他参数

S=0.0227;L=1.8;SONIC=343.13;RHO=1.225;

%四阶龙格-库塔法子函数

functionrk_4（n,h）

globalydy;

dy=zeros（n,1）;

old_y=zeros（1,n）;

y1=zeros（1,n）;

a

（1）=h/2;

a

（2）=h/2;

a（3）=h;

a（4）=h;

dery（y）;

fori=1:

n

old_y（i）=y（i）;

end

forj=1:

3

fori=1:

n

y1（i）=old_y（i）+a（j）*dy（i）;

y（i）=y（i）+a（j+1）*dy（i）/3;

end

dery（y1）;

end

fori=1:

n

y（i）=y（i）+a

（1）*dy（i）/3;

end

%右端子函数

functiondery（y）

globaldy;

globalLSSONICRHO;

globalmaabs_alphacxcymzafmzwzjzalphapmc;

aa=zeros（1,4）;

q=RHO*y

（2）*y

（2）/2;

ma=y

（2）/SONIC;

alpha=y（5）-y（3）;

abs_alpha=abs（alpha）;

interp（）;

aa

（1）=sin（y（3）*pi/180）;aa

（2）=cos（y（3）*pi/180）;

aa（3）=sin（alpha*pi/180）;aa（4）=cos（alpha*pi/180）;

ifalpha<0

cy=-cy;end

xf=cx*q*S;yf=cy*q*S;

wzt=y（4）*L/y

（2）;

dy

（1）=1;

dy

（2）=（p*aa（4）-xf-9.81*y（8）*aa

（1））/y（8）;

dy（3）=（p*aa（3）+yf-9.81*y（8）*aa

（2））/（y

（2）*y（8））/pi*180;%转换成度

dy（4）=（mzaf*alpha+mzwz*wzt）*q*S*L/jz;

dy（5）=y（4）/pi*180;

dy（6）=y

（2）*aa

（2）;

dy（7）=y

（2）*aa

（1）;

dy（8）=-mc;

%插值子函数

functioninterp（）

globalacxacyajzamzafamzwzaxgapamcagcandmandafbL;

globalmaabs_alphacxcymzafmzwzxgjzpmc;

globaly;

ify

（1）

%a=interp11（tt,3,y

（1））;

xg=interp11（axg,10,y

（1））;

jz=interp11（ajz,14,y

（1））;

p=interp11（ap,11,y

（1））;

mc=interp11（amc,6,y

（1））;

elsexg=axg（2,10）;

jz=ajz（2,14）;

p=ap（2,11）;

mc=amc（2,6）;

end

cx=interp33（ma,abs_alpha,9,6,andm,andaf,acx）;

cy=interp33（ma,abs_alpha,9,6,andm,andaf,acy）;

mzaf0=interp33（ma,abs_alpha,9,6,andm,andaf,amzaf）;

ifabs_alpha~=0.

mzaf=mzaf0+cy*（xg-agc

（1））/（abs_alpha*L）;%去掉了*RAD

elsemzaf=mzaf0;

end

mzwz=interp33（ma,abs_alpha,9,6,andm,andaf,amzwz）;

%不等距单变元线性插值子函数

functionres=interp11（yy,n,x）

forj=1:

（n-1）

ifx<=yy（1,j+1）

i=j;break;

elsei=n-1;

end

end

res=yy（2,i）+（yy（2,i+1）-yy（2,i））*（x-yy（1,i））/（yy（1,i+1）-yy（1,i））;

%等距双变元抛物线线性插值子函数

functionres=interp33（x,qq,n1,n2,a,bb,yy）

h=（bb

（2）-bb

（1））/（n2-1）;

i=fix（（qq-bb

（1））/h+1）;

if（i-1）<0

i=1;

elseif（i-n2）>=0

i=n2-1;

end

yy1=interp31（x,n1,