导弹总体设计《防空导弹主级燃料相对质量因数计算》
一、要求：
1、 根据已知条件，采用数值积分法求解相对量运动微分方程组，计算燃料相对质量因数。 2、 综合运用积分、插值计算等计算方法，采用C、C++，或者Matlab等语言的一种，编制
计算程序。 二、已知条件：
1、分离条件：速度V0=500m/s，时间t0=3s，x方向位置x0=674m，y方向位置y0=329m，初始攻角?0?1.5，初始弹道倾角?0?26。 阻力因数Cx(Ma,?)插值表：
?（） Ma 1.4 1.5 2.1 2.7 3.3 4.0 1 0.0401 0.0389 0.0322 0.0276 0.0236 0.0197 2 0.0442 0.0430 0.0360 0.0308 0.0265 0.0222 4 0.0524 0.0511 0.0436 0.0372 0.0323 0.0272 6 0.0667 0.0651 0.0558 0.0481 0.0419 0.0356 8 0.0866 0.0847 0.0736 0.0641 0.0560 0.0478 10 0.1145 0.1120 0.0973 0.0849 0.0746 0.0644 升力因数斜率C?y(Ma,?)插值表：
?（） Ma 1.4 1.5 2.0 2.5 3.0 3.5 4.0
注：以上攻角单位均为度。
1 0.0307 0.0302 0.0279 0.0261 0.0247 0.0226 0.0209 2 0.0309 0.0304 0.0280 0.0264 0.0248 0.0227 0.0210 4 0.0310 0.0306 0.0284 0.0267 0.0251 0.0231 0.0213 6 0.0314 0.0309 0.0286 0.0269 0.0254 0.0233 0.0216 8 0.0316 0.0311 0.0288 0.0272 0.0257 0.0236 0.0219 10 0.0318 0.0313 0.0290 0.0274 0.0259 0.0238 0.0221 3、 发动机参数：比冲Is=2156N?s/kg；重力加速度g=9.801；推重比P?2.2；翼载。
p0?5880N/m2
4、 导引规律：三点法。目标匀速直线等高迎头飞行，VT?420m/s，yT?15000m，
DT0?34200m（弹目斜距）。
注：因初始马赫数为1.4706，初始攻角为0.0262，为便于程序编写，故先采用双线性插值法对原题目阻力因数、升力因数斜率数据扩充。
基础一：所采用插值方法是双线性差值，对其原理简单介绍如下：
又称为双线性内插，是有两个变量的插值函数的线性插值扩展，其核心思想是在两个方向分别进行一次线性插值。首先在x轴方向上，对R1和R2两个点进行插值，然后根据R1和R2对P点进行插值，公式如下：
f(x,y)??f(Q11)f(Q21)(x2?x)(y2?y)?(x?x1)(y2?y)(x2?x1)(y2?y1)(x2?x1)(y2?y1)f(Q12)f(Q11)(x2?x)(y?y1)?(x?x1)(y?y1)(x2?x1)(y2?y1)(x2?x1)(y2?y1)
基础二：考虑如下一阶微分方程初值问题
?dy??f(x,y),a?x?b?dx ??y(a)?y0由显式单步法的四级四阶经典Runge?Kutta公式，可迭代求解上述微分方程，
h??yn?1?yn?6(K1?2K2?2K3?K4)??K1?f(xn,yn)?hh?K?f(x?,y?K1)?2nn22?hh?K?f(x?,y?K2)nn?322???K4?f(xn?h,yn?hK3)具体求解：以下说明对本题中相对量运动微分方程组的求解：（使用C++编写程序，具体
程序见.cbp文件）
?dvIs?v2CxIsIs???sin??d?1??2Pp(1??)P0??vTIs1dydv[?(v?y)]?d?yPgvsin?d?d???T?IsvT?d?1?(cotq?)cot?0?PgyT?Is?dy?vsin???d?Pg??dx?Isvcos??d?Pg?d?Iscos??v??d?P????v2C?IyTRIs??s?2Pp0(1??)1????Ma?v?c?dvd?dxdyf?,g?,b?,m?Runge?Kutta 模块1：函数实现，其中d?d?d?d?void RK4(double(*f)(double u,double v,double theta,double y),double(*g)(double u,double v,double theta,double y),double(*b)(double v,double theta),double(*m)(double v,double theta),double initial[5],double resu[5],double h) {
double f1,f2,f3,f4,g1,g2,g3,g4,b1,b2,b3,b4,m1,m2,m3,m4; double u0,v0,theta0,x0,y0,v1,theta1,x1,y1;
u0=initial[0];v0=initial[1];theta0=initial[2];x0=initial[3];y0=initial[4]; f1=f(u0,v0,theta0,y0); g1=g(u0,v0,theta0,y0); b1=b(v0,theta0); m1=m(v0,theta0);
f2=f(u0+h/2,v0+h*f1/2,theta0+h*g1/2,y0+h*m1/2); g2=g(u0+h/2,v0+h*f1/2,theta0+h*g1/2,y0+h*m1/2); b2=b(v0+h*f1/2,theta0+h*g1/2); m2=m(v0+h*f1/2,theta0+h*g1/2);
f3=f(u0+h/2,v0+h*f2/2,theta0+h*g2/2,y0+h*m2/2); g3=g(u0+h/2,v0+h*f2/2,theta0+h*g2/2,y0+h*m2/2); b3=b(v0+h*f2/2,theta0+h*g2/2); m3=m(v0+h*f2/2,theta0+h*g2/2);
f4=f(u0+h,v0+h*f3,theta0+h*g3,y0+h*m3); g4=g(u0+h,v0+h*f3,theta0+h*g3,y0+h*m3); b4=b(v0+h*f3,theta0+h*g3); m4=m(v0+h*f3,theta0+h*g3);
v1=v0+h*(f1+2*f2+2*f3+f4)/6;
theta1=theta0+h*(g1+2*g2+2*g3+g4)/6; x1=x0+h*(b1+2*b2+2*b3+b4)/6;
y1=y0+h*(m1+2*m2+2*m3+m4)/6; }
模块2：对相对量运动微分方程组代值及函数实现：
double f(double u,double v,double theta,double y){ double dv;
dv=2156/(1-u)-(2156/25872)*Cx*1.225*pow(1-y/44330.8,4.2559)*v*v/(1-u)-980*sin(theta);//大气密度=1.225*pow(1-y/4.3308,4.2559) return(dv); }
double g(double u,double v,double theta,double y){ double num,den;
den=1+(2.1057455-2.7997*u)/tan(theta);//第二方程的分母项
num=2*420*2156/(15000*2.2*9.801)-(420/15000)*y*(1/v*v)*(1/sin(theta))*(2156/(1-u)-(2156/25872)*Cx*1.225*pow(1-y/44330.8,4.2559)*v*v/(1-u)-980*sin(theta));//第二方程的分子项 dth=num/den; return(dth); }
double b(double v,double theta){ double dx;
dx=v*cos(theta)*(2156/(2.2*9.801));//取重力加速度g=9.801 return(dx); }
double m(double v,double theta){ double dy;
dy=v*sin(theta)*(2156/(2.2*9.801)); return(dy); }
模块3：包含对阻力因数、升力因数斜率的插值求解；sa[7][6]存放Cx;sx[6][6存放cya；
for(j=0;j<7;j++){
if(sa_[j][0]>Ma) break;} sa_[0][0]=1/57.3; for(k=0;k<6;k++){
if(sa_[0][k]>arf) break;}
cya=sa[j-1][k-1]*(sa_[j][0]-Ma)*(sa_[0][k]-arf)+sa[j][k-1]*(Ma-sa_[j-1][0])*(sa_[0][k]-arf)+sa[j-1][k]*(sa_[j][0]-Ma)*(arf-sa_[0][k-1])+sa[j][k]*(Ma-sa_[j-1][0])*(arf-sa_[0][k-1]); //双线性插值求升力系数斜率cya arf1=resu[1]*dth+2156*cos(resu[2])/2.2;
arf2=(1.225*pow(1-resu[4]/44330.8,4.2559))*resu[1]*resu[1]*(2156/25872)*cya/(1-resu[0])+2156/(1-resu[0]);
arf=arf1/arf2;
//求攻角arf,arf=arf1/arf2; for(j=0;j<6;j++){
if(sx_[j][0]>Ma) break;} sx_[0][0]=1/57.3; for(k=0;k<6;k++){
if(sx_[0][k]>arf) break;}
Cx=sx[j-1][k-1]*(sx_[j][0]-Ma)*(sx_[0][k]-arf)+sx[j][k-1]*(Ma-sx_[j-1][0])*(sx_[0][k]-arf)+sx[j-1][k]*(sx_[j][0]-Ma)*(arf-sx_[0][k-1])+sx[j][k]*(Ma-sx_[j-1][0])*(arf-sx_[0][k-1]); } //双线性插值求升力系数斜率Cx
模块4：函数主体，对Runge?Kutta和相对量微分方程的反复调用和迭代：
运行程序后按照示例输入数据 会得到如下结果：
导弹满足与目标的的遭遇条件时，即导弹高度达到15000米时，燃料相对质量因数约为29.4%，即为所求，且此时v?1072.5m/s,??57.13o。
西工大导弹总体设计主级燃料相对质量因数计算大作业
