弹道仿真.docx

弹道仿真.docx

《弹道仿真.docx》由会员分享，可在线阅读，更多相关《弹道仿真.docx（13页珍藏版）》请在冰豆网上搜索。

弹道仿真

导弹无控弹道仿真

蒋洋

（北京理工大学宇航学院）

摘要：

在进行导弹总体设计时，会进行无控和有控弹道仿真来验证导弹的设计是否合理，利用matlab对导弹运动方程组进行求解可以得到导弹的无控弹道曲线等。

关键词：

无控弹道；matlab；导弹运动方程组

Abstract:

Inthedesignofthemissile,wewilldouncontrolledandcontrolledballistictrajectorysimulationtoverifywhetherthedesignisreasonableornot.Withthehelpofmatlab,wecaneasilysolvetheequationsofmotionofamissiletogetuncontrolledtrajectorycurve.

Keywords:

Uncontrolledtrajectory;matlab;theequationsofmotionofamissile

1、数学模型及公式

由于四阶龙格-库塔法精度高且易于编写程序，所以在本次仿真实验中求解导弹运动方程组将使用龙格-库塔法，龙格-库塔法的运算公式如下：

在实验中使用的数学模型如下方程所述

以上数学模型满足以下假设：

1）侧向运动参数

及舵偏角

都比较小。

这样可以令：

且略去小量的乘积

……以及参数

对阻力X的影响。

2）导弹基本上在铅垂面内飞行，即其弹道与铅垂面弹道差别不大。

3）俯仰操纵机构的偏转仅取决于纵向运动参数；而偏航、倾斜操纵机构的偏转仅取决于侧向运动参数。

其他相关公式：

空气动力阻力

升力

空气动力矩

其中，

2、原始数据

1）初值

X=0（m）y=20（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

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

Jy（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）阻尼力矩导数

当

=0.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

当

=0.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）其他参数

特征面积（m2）

特征长度（m）

毛翼展（m）

音速（m/s）

大气密度（kg/m3）

0.0227

1.8

0.5

343.13

1.225

3、Matlab程序编写

1）主程序main

functionmain

h=input（'输入h='）;

t=input（'输入t='）;%输入步长h和时间t

x=[20;0.314;0;0.314;0;20;52.38];

%按照导弹运动方程顺序x=[v,,z,,x,y,m]

n=fix（t/h）;

fori=1:

n

RES（:

i）=x;%将数据存入RES

x=rk4（h*i,x,h）;

ifRES（6,end）<0

n=i;

break

end

end

time=num2str（i*h）;

disp（['飞行时间=',time,'s']）;

leng=num2str（RES（5,end））;

disp（['飞行距离=',leng,'m']）;%输出运算结果

picture（RES,h,n）%画图

2）画图程序picture

functionpicture（RES,h,n）

figure

（1）

plot（RES（5,:

）,RES（6,:

）,'Linewidth',1）

title（'弹道曲线','FontSize',15）;ylim（[0160]）;

xlabel（'x/m','FontSize',15）,ylabel（'y/m','FontSize',15）

grid%利用RES中的x，y的数据画出弹道曲线

figure

（2）

plot（0:

h:

（n-1）*h,（RES（4,:

）-RES（2,:

））*57.3,'Linewidth',1）;

title（'攻角时间曲线','FontSize',15）;

xlabel（'t/s','FontSize',15）,ylabel（'{\alpha}/（\circ）','FontSize',15）;

grid%利用RES中的alpha数据画出以时间t为横轴的曲线

figure（3）

plot（0:

h:

（n-1）*h,RES（1,:

）,'Linewidth',1）

title（'速度时间曲线','FontSize',15）;

xlabel（'t/s','FontSize',15）,ylabel（'V/（m*s^{-1}）','FontSize',15）

grid%利用RES中速度V的数据画出速度时间曲线

figure（4）

plot（0:

h:

（n-1）*h,RES（2,:

）,'Linewidth',1）

title（'弹道倾角曲线','FontSize',15）;

xlabel（'t/s','FontSize',15）,ylabel（'{\theta}/（rad）','FontSize',15）

grid%利用RES中的

的数据画出弹道倾角曲线

figure（5）

plot（0:

h:

（n-1）*h,RES（4,:

）,'Linewidth',1）

title（'俯仰角曲线','FontSize',15）;

xlabel（'t/s','FontSize',15）,ylabel（'{\vartheta}/（rad）','FontSize',15）

grid%利用RES中的

数据画出俯仰角曲线

figure（6）

plot（0:

h:

（n-1）*h,RES（3,:

）,'Linewidth',1）

title（'俯仰角速度曲线','FontSize',15）;

xlabel（'t/s','FontSize',15）,ylabel（'{\omega}_{z}/（rad*s^{-1}）','FontSize',15）

grid%利用RES中的

数据画出俯仰角速度曲线

3）四阶龙格-库塔法程序rk4

functionxadd=rk4（t,x,h）

K1=f（t,x）;

K2=f（t+h/2,x+h/2*K1）;

K3=f（t+h/2,x+h/2*K2）;

K4=f（t+h,x+h*K3）;

xadd=x+h/6*（K1+2*K2+2*K3+K4）;

4）插值以及运算程序f

functiondy=f（t,x）

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

（以矩阵方式输入原始数据，由于数据太多，此处省略）

P1=interp1（tP（1,:

）,tP（2,:

）,t,'linear'）;

tmc=interp1（tMc（1,:

）,tMc（2,:

）,t,'linear'）;

jy=interp1（tJy（1,:

）,tJy（2,:

）,t,'linear'）;

XG1=interp1（tXG（1,:

）,tXG（2,:

）,t,'linear'）;

%使用interp1函数对数据进行线性插值

alpha=abs（180*（x（4）-x

（2））/pi）;Ma=x

（1）/SONIC;

Cx=interp2（0:

2:

10,0:

0.1:

0.9,aCx,alpha,Ma,'cubic'）;

Cy=interp2（0:

2:

10,0:

0.1:

0.9,aCy,alpha,Ma,'cubic'）;

%使用interp2函数对升力和阻力系数矩阵插值

ifx（4）-x

（2）<0%判断攻角是否为负来对升力系数和静稳定力矩进行变换正负处理

Cy=-Cy;

mza=-mza;

end

T=1/2*RHO*x

（1）^2*S;

X=Cx*T;

Y=Cy*T;%计算空气动力

P=P1;

mc=tmc;

Jz=jy;

mza0=interp2（0:

2:

10,0:

0.1:

0.9,mza,alpha,Ma,'cubic'）;

ifalpha~=0

mza=mza0+Cy*（XG1-0.9381）/（L）;

else

mza=mza0;

end

mzw1=interp2（0:

2:

10,0:

0.1:

0.9,mzw9381,alpha,Ma,'cubic'）;

mzw2=interp2（0:

2:

10,0:

0.1:

0.9,mzw8896,alpha,Ma,'cubic'）;

%使用interp2函数对阻尼力矩导数矩阵插值

mzw=interp1（[0.9381,0.8896],[mzw1,mzw2],XG1,'linear'）;

dy=zeros（7,1）;

dy

（1）=（P*g*cos（x（4）-x

（2））-X-x（7）*g*sin（x

（2）））/x（7）;

dy

（2）=（P*g*sin（x（4）-x

（2））+Y-x（7）*g*cos（x

（2）））/（x（7）*x

（1））;

dy（3）=（mza/Jz+mzw/Jz*x（3）*L/x

（1））*T*L;

dy（4）=x（3）;

dy（5）=x

（1）*cos（x

（2））;

dy（6）=x

（1）*sin（x

（2））;

dy（7）=-mc;%导弹运动方程组

4、运行结果

1）数值结果

输入h=0.005

输入t=12

飞行时间t=11.535s

飞行距离L=2658.2788m

2）图象