导弹总体设计《防空导弹主级燃料相对质量因数计算》大作业

导弹总体设计《防空导弹主级燃料相对质量因数计算》大作业班号：_________学号：_________姓名：_________一、要求：1、根据已知条件，采用数值积分法求解相对量运动微分方程组，计算燃料相对质量因数。2、综合运用积分、插值计算等计算方法，采用C、C++，或者Matlab等语言的一种，编制计算程序。二、已知条件：1、分离条件：速度V0=500m/s，时间t0=3s，x方向位置x0=674m，y方向位置y0=329m，初始攻角01.5，初始弹道倾角026。2、气动参数阻力因数(,)xCMa插值表：（）Ma2468101.50.04300.05110.06510.08470.11202.10.03600.04360.05580.07360.09732.70.03080.03720.04810.06410.08493.30.02650.03230.04190.05600.07464.00.02220.02720.03560.04780.0644升力因数斜率(,)yCMa插值表：（）Ma12468101.50.03020.03040.03060.03090.03110.03132.00.02790.02800.02840.02860.02880.02902.50.02610.02640.02670.02690.02720.02743.00.02470.02480.02510.02540.02570.02593.50.02260.02270.02310.02330.02360.02384.00.02090.02100.02130.02160.02190.0221注：以上攻角单位均为度。3、发动机参数：比冲Is=2156/Nskg；重力加速度g=9.801；推重比2.2P；翼载205880/pNm。4、导引规律：三点法。目标匀速直线等高迎头飞行，420/TVms，15000Tym，（Ty=14200m），034200TDm（0TD=32000m）（弹目斜距）5、大气模型：参见教科书。三、计算思路①插值表的处理1.对于大气参数密度ρ、声速a，都是高度h的单值函数，本次计算采用最方便的线性插值，对于指定高度h，找出上下最接近的两点h1,h2进行计算。2.对于Cx，Cyα它由两个参数决定，采用双线性插值。具体算法如下：将Cx（α，Ma)看作是攻角和马赫数的二元点函数Z（x，y)，则对于任意点p（a，b），找出距离a最近的点A（x1，y1）、B（x2，y2）、C（x3，y3）、D（x4，y4），如上图示。先用P点的横坐标在x轴方向（A,B方向）上进行第一次线性插值，计算出E点的值，同样在C,D方向进行第二次插值求出F点的值，最后在y轴方向（E，F方向）用纵坐标进行线性插值求出P点的函数值。②解方程组1.四阶龙格库塔公式：y(i+1)=y(i)+h*(K1+2*K2+2*K3+K4)/6K1=f(x(i),y(i))K2=f(x(i)+h/2,y(i)+hK1/2)K3=f(x(i)+h/2,y(i)+hK2/2)K4=f(x(i)+h,y(i)+hK3)2.三点法导引的相对量运动微分方程组为：2020200ssin12(1)1sin1cotcotsincoscos2(1)1cotcotsxsssTTsTTsssyTRssTTIVCIIdVdPpPIVdydVVyyPgVddddIVqPgyIdyVdPgIdxVdPgIdVdPVCIIPpPgyqqIV0cotcotTTVMaCytqqV此方程组变量较多，但是只有方程1、方程2需要进行复杂的求解运算，其余的方程均可代入上一步的计算结果求得。以μ为自变量，初始μ=0，选定一个步长h，利用已知的初值条件，用龙格库塔法解方程1、方程2，方程3~8可由1、2的结果代入计算求得，经过多次运算求解。具体步骤如下：已知初始参数g=9.801m/s^2，μ=0，Is=2156N*s/kg，P=2.2，p0=5880N/m^2，α0=1.5°，θ=26°，vt=vt0=420m/s，yt0=15000m，dt0=34200m，易得xt0=31569.39345，cotq0=xt0/yt0=2.10462623，h=0.001。1.利用大气参数表计算出当前高度y=y（h)（初值y0=329m）下的密度ρ,声速a；2.利用当前速度（初值v0=500m/s）和声速解方程7计算出当前马赫数Ma；3.解方程3，4，得dx/du，dy/du；4.用龙格库塔法解方程1，得dv/du；5.用龙格库塔法解方程2，得dθ/du；6.将上述方程的结果代入方程5、6、8，解得α、q，t；7.用新值代替旧值重复1~6步骤，判断y=y0时，计算结束。四、计算结果击中目标时，μ=0.346，t=37.49648，x=16993.632，y=15058.591，v=849.348，θ=60.776°，α=5.439°，q=41.289°，xt=17080.872，Ma=2.873.五、C语言计算程序代码#includestdio.h#includestring.h#includemath.h#definePI3.1415926535//圆周率π#defineh0.001//步长#definev0500.0//分离时导弹初速度v0单位：s#definet03.0//分离时刻t0，单位：s//#definex0674.0//#definey0329.0#definealpha01.5*PI/180//攻角α，单位：rad#definetheta026.0*PI/180//俯仰角θ，单位：rad#defineis2156.0//总冲Is,单位：N*s/kg#defineg9.801//重力加速度g,单位：m/s^2#defineP2.2//推重比P#definep05880.0//翼载p0，单位：N/m^2#definevt420.0//飞机速度vt，单位：m/s#defineyt015000.0//飞机初始纵坐标yt，单位：m#definedt034200.0//弹目斜距Dt0，单位：m#definext031569.39345//飞机初始橫坐标xt0，单位：mxt0=sqrt(dt0^2-(yt0-y0)^2)+x0#definecotq02.10462623//初始目标对制导站的航向角的余切，cotq0=xt0/yt0;doubleinter_linear(doublex0,doublex1,doubley0,doubley1,doublex){//线性插值函数doublea0,a1,y;a0=(x-x1)/(x0-x1);a1=(x-x0)/(x1-x0);y=a0*y0+a1*y1;returny;}doubleinterp1(doublex[],doubley[],intn,doublea){//一维线性插值函数doublez;inti=0;for(i=0;i(n-1);i++){if((a=x[i])&&(a=x[i+1]))break;//确定a的位置}z=inter_linear(x[i],x[i+1],y[i],y[i+1],a);returnz;}doubleinterp2(doublex[],doubley[],doublez[6][6],intm,intn,doublea,doubleb){//双线性插值函数inti,j;doublew,w1,w2;for(i=0;i(m-1);i++){//确定a在x轴的位置if((a=x[i])&&(a=x[i+1]))break;}for(j=0;j(n-1);j++){//确定b在y轴的位置if((b=y[j])&&(b=y[j+1]))break;}if(ax[0])i=0;//若给定的点超出范围时，采用线性外推进行插值elseif(ax[m-1])i=m-1;if(by[0])j=0;elseif(by[n-1])j=n-1;w1=inter_linear(x[i],x[i+1],z[i][j],z[i+1][j],a);//在x轴方向进行第一次插值w2=inter_linear(x[i],x[i+1],z[i][j+1],z[i+1][j+1],a);//在x轴方向进行第二次插值w=inter_linear(y[j],y[j+1],w1,w2,b);//在y轴方向进行插值returnw;}intmain(){//*****************************计算用到的参数**************************************doublecx,cya;//阻力系数Cx，升力系数斜率Cyαdoublealpha1[5]={2,4,6,8,10},ma1[5]={1.5,2.1,2.7,3.3,4.0};//Cx插值表中的α向量,Ma向量doublealpha2[6]={1,2,4,6,8,10},ma2[6]={1.5,2.0,2.5,3.0,3.5,4.0};//Cyα插值表中的α向量，Ma向量doublecx0[6][6]={{0.0430,0.0511,0.0651,0.0847,0.1120},\{0.0360,0.0436,0.0558,0.0736,0.0973},\{0.0308,0.0372,0.0481,0.0641,0.0849},\{0.0265,0.0323,0.0419,0.0560,0.0746},\{0.0222,0.0272,0.0356,0.0478,0.0644}};//Cx插值表doublecy0[6][6]={{0.0302,0.0304,0.0306,0.0309,0.0311,0.0313},\{0.0279,0.0280,0.0284,0.0286,0.0288,0.0290},\{0.0261,0.0264,0.0267,0.0269,0.0272,0.0274},\{0.0247,0.0248,0.0251,0.0254,0.0257,0.0259},\{0.0226,0.0227,0.0231,0.0233,0.0236,0.0238},\{0.0209,0.0210,0.0213,0.0216,0.0219,0.0221}};//Cyα插值表doubleyh[23]={0,1000,2000,3000,4000,5000,6000,7000,8000,\9000,10000,11000,12000,13000,14000,15000,\16000,17000,18000,19000,20000,21000,22000};//大气高度y(h)doubleph[23]={1.22505,1.11168,1.00646,0.90913,0.81913,0.73612,0.65969,\0.58950,0.52517,0.46635,0.41270,0.36391,0.31083,0.26549,\0.22675,0.19367,0.16542,0.14128,0.12068,0.10307,0.08803,\0.07487,0.06373};//大气密度表ρ(h)doubleah[23]={340.29,336.43,332.53,328.58,324.58,320.53,316.43,312.27,\308.06,303.79,299.46,295.07,295.07,295.07,295.07,295.07,\295.07,295.07,295.07,295.07,295.07,295.75,296.43};//大气声速表a(h)inti;doublep,a,q,cotq;doublex=6