2009302026王勋

航空外弹道学课程设计姓名：王勋班级： 08030902班学号： 2009302026时间： 2012.7.12一、已知条件及题目要求1.查表可知，标准下落时间()21.12S s s2．气象：760on h mmHg =288.4on K τ=29.27R =温度梯度35.86210G -=⨯度/米8.4on mmHg =e3. 弹丸参数：216.5q kg =0.299d m =弹长 2.11L m =9.806g =4. 空气阻力系数：00.160x C =5. 初始条件：400/u m s =0w =0p =2000m H =要求：列出弹道参数，并画出炸弹弹道曲线。

二、题目分析对于所给题目，取直角坐标系Oxy ，坐标原点取在投弹高度为H 的投弹点O 上，x 轴取在飞机投弹瞬间速度1v →的铅垂面（投弹面）内的水平方向，y 轴铅直向下，如图1所示。

转角方向规定顺时针为正。

设弹道上任一点M 速度向量v →在x 轴与y 轴上的分量分别为u 与w 。

图1 水平轰炸图由题目给出的条件，P=0，W=0可以知道，所要设计的题目类型为无外力的水平轰炸，由于飞机速度V 在X 轴方向上，因而初始条件为：当0t =时，u v,w 0,x=0,y 0,0.θ====建立方程组如下：()G(v )1=()G(v )+cos equations(1)sin du CH H y dt d CH H y g dt dx u v dt dy v dt v τττττωωθωθ⎧=--⎪⎪⎪--⎪⎪⎪==⎨⎪⎪==⎪⎪⎪⎪⎩（）（2）（3）（4）（5） 三、相关参数计算(1)弹道系数C 的计算 0.2460.650.16C C i x0x0===2233id 0.2460.299C 10100.10158216.5q ⨯=⨯=⨯= 其中d 为弹直径，q 为炸弹的重量。

(2)288.4G*(2000)y τ=+-，其中310*862.5G -=(3)查表得)Ma (C x0，其中：v Ma a v *(288.40.005862*2000-0.005862y)^0.520.074v *(300.124-0.005862y)^0.520.074==+= Ma 可作如下近似)Ma (C x0=（4）(v )G τ的计算onH τγγ 其中0()x v C a由（3）可求得，τ由（2）计算可求得。

（5）(y)H τ的计算onH τγγ (6)气压的计算：1/RG ON ON Gh h (1*(H-y))τ=- 其中-3G 5.86210=⨯(7)音速的计算：a =(8)密度的计算：ON 1.225γ=τγR h 6.13=ON )(H γγ=y其中γ为空气比重，ON γ为地面空气比重标准值。

四、方程变形由如下方程：()G(v )1=()G(v )+cos equations(2)sin du CH H y dt d CH H y g dt dx u v dt dy v dt v τττττωωθωθ⎧=--⎪⎪⎪--⎪⎪⎪==⎨⎪⎪==⎪⎪⎪⎪⎩（）（2）（3）（4）（5） 令w P u= 把上述方程组改变为：()() equations(3)du CH H y G v udt dP gdt u dx udt dy uPdt v ττ⎧=--⎪⎪⎪=⎪⎪⎪=⎨⎪⎪=⎪⎪⎪=⎪⎩积分的初始条件：当0t =时，1,0,0,0u v P x y ====带入上述转换化简得：1404.73710*[1()]*()* equations(4)ON RG x ON du G C H y C Ma uv dtdP g dt u dy uP dt dx u dtτττ-⎧=-⨯--⎪⎪⎪=⎪⎪⎨⎪=⎪⎪⎪=⎪⎩ （1）（2）（3）（4）五、运用RungeKutta 法求解RungeKutta 法的主要优点是计算精度较高，能满足通常的计算要求；每次计算，只用到前一步的计算结果；计算过程中改变步长不受限制。

求解C 代码如下：#include <math.h>#include <stdio.h>#include<stdlib.h>#define H 2000#define G 5.682e-3#define R 29.27#define tao0 288.4#define C 0.10165#define U 400#define g 9.806#define dt 0.01//求相关参数；double tao(double y2) //求τ{ double tao_;tao_=tao0-G*(H-y2);return(tao_);}double Hg(double y2) //比重函数Hg{ double Hg_;Hg_=pow((1-G/tao0*(H-y2)),1/(R*G))*sqrt(tao0/tao(y2));return(Hg_);}double V(double y0,double y1) //速度{ double V_;V_=sqrt(y0*y0+y1*y1);return(V_);}double Cx0(double y0,double y1,double y2) //弹丸阻力系数Cx0{ double a;double Cx0_;a=20.074*sqrt(tao(y2));//音速为aif(V(y0,y1)/a<0.8)Cx0_=0.6;else Cx0_=V(y0,y1)/a-0.1;return(Cx0_);}double Gt(double y0,double y1,double y2) //空气阻力函数{ double Gt_;Gt_=4.737e-4*V(y0,y1)*Cx0(y0,y1,y2)*sqrt(tao0/tao(y2));return(Gt_);}//x速度u，y为P，z为高度double f1(double x, double y, double z) //方程组（1）式{return (-1)*C*Hg(z)*Gt(x,x*y,z)*x;}double f2( double x) //方程组（2）式{return g/x;}double f3( double x, double y) //方程组（3）式{return x*y;}double f4(double x) //方程组（4）式{return 1.0*x;}double RungeKutta(double x0, double y0, double z0,double m0, double a, double b, int n,double x[], double y[], double z[],double m[], double(*f1)(double, double, double),double (*f2)(double), double (*f3)(double, double), double (*f4)(double )) //RungeKutta法{double k1, k2, k3, k4 ,m1, m2, m3, m4, n1, n2, n3, n4, c1, c2, c3, c4 ,h,t; int count;count =0;h =(b-a)/n;x[0] = x0; //初值y[0] = y0;z[0] = z0;m[0] = m0;FILE *file= fopen(“result1.txt”,”w”);for(int i = 1;i<=n;i++){t = a+(i-1)*h;k1 = f1(x[i-1], y[i-1], z[i-1]);m1 = f2(x[i-1]);n1 = f3(x[i-1], y[i-1]);c1 = f4(x[i-1]);k2 = f1(x[i-1]+k1*h/2.0, y[i-1]+m1*h/2.0, z[i-1]+n1*h/2.0);m2 = f2(x[i-1]+k1*h/2.0);n2 = f3(x[i-1]+k1*h/2.0, y[i-1]+m1*h/2.0);c2 = f4(x[i-1]+k1*h/2.0);k3 = f1(x[i-1]+k2*h/2.0, y[i-1]+m2*h/2.0, z[i-1]+n2*h/2.0);m3 = f2(x[i-1]+k2*h/2.0);n3 = f3(x[i-1]+k2*h/2.0, y[i-1]+m2*h/2.0);c3 = f4(x[i-1]+k2*h/2.0);k4 = f1(x[i-1]+k3*h/2.0, y[i-1]+m3*h/2.0, z[i-1]+n3*h/2.0);m4 = f2(x[i-1]+k3*h/2.0);n4 = f3(x[i-1]+k3*h/2.0, y[i-1]+m3*h/2.0);c4 = f4(x[i-1]+k3*h/2.0);x[i] = x[i-1]+h*(k1+2*k2+2*k3+k4)/6.0;y[i] = y[i-1]+h*(m1+2*m2+2*m3+m4)/6.0;z[i] = z[i-1]+h*(n1+2*n2+2*n3+n4)/6.0;m[i] = m[i-1]+h*(c1+2*c2+2*c3+c4)/6.0;++count;fprintf(file,"当t=%4.2f时，速度x= %f m/s p=w/u= %f 高度h= %f 水平距离m= %f\n",t+dt,x[i],y[i],z[i],m[i]);}fprintf(file,"计算次数为 count=%d\n", count);fclose(file);return 0;}void main(){double x[210];double y[210];double z[210];double m[210];RungeKutta(U,0.0,0.0,0.0,0.0,22.12,200,x,y,z,m, f1, f2, f3, f4); }运行结果截图：图2 运行结果参数六、运用MATLAB绘制弹道图利用离散点绘制出弹道曲线后，考虑使用插值法求得落点坐标，并在弹道曲线上标示如图。

Matlab代码如下：[x y] =textread('D:\result.txt','%f%f');% Set up fittype and options.ft = 'linearinterp';opts = fitoptions( ft );opts.Normalize = 'on';% Fit model to data.fitresult = fit( y, x, ft, opts );yy1 = 2000; % y = 2000，要求的落点xx1 = fitresult( yy1 ); % 与y对应的x值% %%% 画点标注plot(x,y,xx1,yy1,'rs','MarkerEdgeColor','k','MarkerFaceColor','g','Ma rkerSize',6)text(xx1,yy1,[num2str(xx1,'落点坐标\n(%5.2f'),',',num2str(yy1,'%5.2f'),')'],'EdgeColor','red','Backgroun dColor',[.7 .9 .7],'VerticalAlignment','bottom');title('水平轰炸弹道曲线');set(gca,'ydir','reverse','xaxislocation','top');运行结果如下图：图3 水平轰炸弹道曲线七、结论通过编程计算出弹道后，运用matlab 绘图得到的弹道曲线，使用龙格库塔法计算弹道误差较小，得到曲线很好。