鉛垂平面飛行彈道仿真及分析

1、數(shù)值實(shí)驗(yàn)作業(yè)1 實(shí)驗(yàn)背景實(shí)驗(yàn)名稱(chēng)：鉛垂平面飛行彈道仿真及分析實(shí)驗(yàn)內(nèi)容與要求：根據(jù)描述飛行器在鉛垂平面內(nèi)運(yùn)動(dòng)的數(shù)學(xué)模型，編制某導(dǎo)彈的鉛垂平面無(wú)控飛行彈道仿真程序，利用計(jì)算機(jī)解算初始段無(wú)控飛行彈道，對(duì)初始段彈道參數(shù)的變化規(guī)律進(jìn)行分析。2 建立數(shù)學(xué)模型：3計(jì)算方法研究確定數(shù)值積分方法和積分步長(zhǎng)使用算法：四階龍格庫(kù)塔法。使用步長(zhǎng)：h=0.0054 原始數(shù)據(jù)：1）初值x=0(m)y=20.0(m) J=18 q=18v=20(m/s)wz=0(rad/s)m=52.38(kg)2）攻角與馬赫數(shù)范圍(僅用于插值計(jì)算)攻角=010 馬赫數(shù)=0.1 0.93）阻力系數(shù)表馬赫數(shù)攻角()02468100.1.41

2、77.4404.5219.6603.85341.10230.2.3858.4086.4903.6290.82261.07230.3.3779.4007.4827.6218.81601.06660.4.3785.4015.4838.6234.81841.07000.5.3787.4018.4846.6249.82091.07380.6.3829.4062.4897.6310.82841.08350.7.3855.4091.4934.6363.83581.09380.8.4082.4321.5175.6621.86411.12540.9.4947.5192.6073.7571.96721.2392

3、4）升力系數(shù)表馬赫數(shù)攻角()02468100.1.0000.64301.47582.28703.07133.84630.2.0000.64541.48072.29423.08143.85980.3.0000.64801.48582.30143.09153.87310.4.0000.65121.49232.31073.10393.88910.5.0000.65541.50072.32273.11973.90920.6.0000.66171.51342.34093.14363.94010.7.0000.66981.53042.36613.17753.98350.8.0000.67921.55012

4、.39503.21624.03230.9.0000.69331.59352.47063.32734.17905）推力數(shù)據(jù)t(s).000.15.492.112.273.538.7825.4542.8043.6844.08P(kgf)331.2614.3505.4607.848.6543.9742.0141.0040.8040.792.22第一級(jí)工作結(jié)束時(shí)間：2.1126s，第二級(jí)工作結(jié)束時(shí)間：44.0832s 6）發(fā)動(dòng)機(jī)質(zhì)量秒流量t(s)0.2.12.10544.144.105100秒流量(kg/s)2.3622.3620.210590.210590.0. 7）轉(zhuǎn)動(dòng)慣量t(s).02.02.4

5、6.410.414.418.422.426.430.434.038.442.444.0Jz(kgms)8.357.887.867.817.787.757.737.717.707.707.697.697.697.69 8）導(dǎo)彈重心（起自頭部）t(s).02.02.410.018.026.032.038.042.044.0XG(m).9381.9095.9091.9026.8969.8928.8907.8896.8895.88969）靜穩(wěn)定力矩系數(shù)馬赫數(shù)攻角()02468100.10.0000 -0.0104 -0.0341 -0.0564 -0.0771 -0.09850.20.0000 -0.

6、0104 -0.0341 -0.0564 -0.0770 -0.09830.30.0000 -0.0104 -0.0341 -0.0564 -0.0769 -0.09820.40.0000 -0.0105 -0.0342 -0.0564 -0.0768 -0.09790.50.0000 -0.0104 -0.0339 -0.0560 -0.0761 -0.09690.60.0000 -0.0093 -0.0314 -0.0521 -0.0708 -0.09030.70.0000 -0.0080 -0.0286 -0.0477 -0.0650 -0.08290.80.0000 -0.0065 -

7、0.0252 -0.0425 -0.0578 -0.07390.90.0000 -0.0053 -0.0229 -0.0391 -0.0538 -0.0693當(dāng)導(dǎo)彈重心變化時(shí)的修正公式：10）阻尼力矩導(dǎo)數(shù)當(dāng)Xg=.9381時(shí)馬赫數(shù)攻角()02468100.1-0.4686 -0.4829 -0.4982 -0.5130 -0.5272 -0.54090.2-0.4707 -0.4850 -0.5003 -0.5150 -0.5292 -0.54290.3-0.4744 -0.4886 -0.5039 -0.5186 -0.5327 -0.54640.4-0.4797 -0.4939 -0.50

8、90 -0.5237 -0.5378 -0.55140.5-0.4882 -0.5022 -0.5173 -0.5318 -0.5458 -0.55930.6-0.5089 -0.5227 -0.5376 -0.5520 -0.5658 -0.57910.7-0.5366 -0.5502 -0.5649 -0.5790 -0.5927 -0.60580.8-0.5738 -0.5871 -0.6014 -0.6153 -0.6287 -0.64150.9-0.6272 -0.6407 -0.6553 -0.6694 -0.6830 -0.6960當(dāng)Xg=.8896時(shí)馬赫數(shù)攻角()0246810

9、0.1-0.6179 -0.6384 -0.6600 -0.6805 -0.6999 -0.71820.2-0.6207 -0.6410 -0.6626 -0.6830 -0.7024 -0.72070.3-0.6253 -0.6455 -0.6670 -0.6874 -0.7067 -0.72490.4-0.6319 -0.6521 -0.6734 -0.6937 -0.7129 -0.73100.5-0.6424 -0.6624 -0.6835 -0.7036 -0.7226 -0.74060.6-0.6669 -0.6866 -0.7074 -0.7272 -0.7459 -0.7636

10、0.7-0.6997 -0.7190 -0.7395 -0.7589 -0.7774 -0.79480.8-0.7435 -0.7624 -0.7824 -0.8014 -0.8194 -0.83650.9-0.8069 -0.8266 -0.8474 -0.8672 -0.8859 -0.903511）其它參數(shù)特征面積S(m2)特征長(zhǎng)度L(m)毛翼展(m)音速SONIC(m/s)大氣密度(kg/m3)0.02271.80.5343.131.225使用的插值算法：氣動(dòng)數(shù)據(jù)插值等距雙變?cè)獟佄锞€(xiàn)插值；推力、重心、轉(zhuǎn)動(dòng)慣量等不等距一元線(xiàn)性插值。5 空氣動(dòng)力和空氣動(dòng)力矩表達(dá)式6 編制計(jì)算程序計(jì)算機(jī)算法

11、采用了matlab實(shí)現(xiàn),源程序見(jiàn)附件.該程序有八個(gè)函數(shù)組成,各函數(shù)之間的調(diào)用關(guān)系如下圖所示.子函數(shù)initl主函數(shù)main子函數(shù)rk_4子函數(shù)result子函數(shù)savedata子函數(shù)drawing子函數(shù)dery子函數(shù)interp子函數(shù)interp33子函數(shù)interp31子函數(shù)interp111) 子函數(shù)initl的功能是輸入求解導(dǎo)彈運(yùn)動(dòng)方程組所需的原始數(shù)據(jù).2) 子函數(shù)rk_4是四階龍格-庫(kù)塔法積分算法子函數(shù),其中調(diào)用了子函數(shù)dery.3) 子函數(shù)dery功能是計(jì)算微分方程組的右端函數(shù),其中調(diào)用了子函數(shù)interp.4) 插值子函數(shù)interp功能是計(jì)算所有需要插值的參數(shù),其中調(diào)用了子函數(shù)i

12、nterp11和interp33.5) 子函數(shù)interp11是不等距單變?cè)€(xiàn)性插值函數(shù),主要用于轉(zhuǎn)動(dòng)慣量Jz, 質(zhì)心位置xg, 推力p等單變量的插值.6) 子函數(shù)interp33是等距雙變?cè)獟佄锞€(xiàn)-線(xiàn)性插值函數(shù),主要用于氣動(dòng)力系數(shù)cx, cy, 氣動(dòng)力矩系數(shù),等雙變?cè)獏?shù)的插值,該函數(shù)調(diào)用了子函數(shù)interp31.7) 子函數(shù)interp31是等距單變?cè)獟佄锞€(xiàn)插值函數(shù),被interp33調(diào)用,完成子函數(shù)interp33的等距雙變?cè)獟佄锞€(xiàn)-線(xiàn)性插值功能.關(guān)于程序中出現(xiàn)的數(shù)組和變量名,作如下說(shuō)明:acx:阻力系數(shù)數(shù)組;acy: 升力系數(shù)數(shù)組;amzaf:某一質(zhì)心位置下的靜穩(wěn)定性導(dǎo)數(shù)數(shù)組;amzw

13、z:阻尼力矩系數(shù)導(dǎo)數(shù)數(shù)組;axg:質(zhì)心位置Xg隨時(shí)間的變化規(guī)律;ajz:轉(zhuǎn)動(dòng)慣量Jz隨時(shí)間的變化規(guī)律;ap:起飛,續(xù)航發(fā)動(dòng)機(jī)的推力值;amc: 起飛,續(xù)航發(fā)動(dòng)機(jī)的燃料質(zhì)量的秒流量值;agc:質(zhì)心位置變化的始末值;andm:氣動(dòng)數(shù)據(jù)插值所需的Ma數(shù)的最小,最大值;andaf: 氣動(dòng)數(shù)據(jù)插值所需的攻角的最小,最大值;y:存放積分結(jié)果的數(shù)組,該數(shù)組在程序開(kāi)始時(shí)存放積分初值;dy:存放右端函數(shù)數(shù)值的數(shù)組;b:存放三個(gè)時(shí)間值的數(shù)組,其中b(2)存放起飛發(fā)動(dòng)機(jī)工作結(jié)束時(shí)間;b(3)存放續(xù)航發(fā)動(dòng)機(jī)工作結(jié)束時(shí)間;程序中一些主要的變量名有:L特征長(zhǎng)度;S特征面積;SONIC聲速C;RHO大氣密度;h積分步長(zhǎng);程

14、序中其他變量都是存放中間結(jié)果的變量.7 計(jì)算結(jié)果運(yùn)行程序得到彈道曲線(xiàn), 速度曲線(xiàn), 攻角曲線(xiàn)如下:通過(guò)運(yùn)用龍格庫(kù)塔法和一些插值方法求解導(dǎo)彈運(yùn)動(dòng)方程組,獲得了導(dǎo)彈各運(yùn)動(dòng)參數(shù)的變化規(guī)律.通過(guò)該算法獲得的變化規(guī)律比較接近真實(shí)情況.附件:源程序function main()clear allclear globalglobal y;global ii;ii=0;h=0.005;initl();while y(7)=0 ii=ii+1; result(ii); rk_4(8,h);endsavedata(ii);drawing();%原始數(shù)據(jù)初始化function initl()global acx a

15、cy ajz amzaf amzwz axg ap amc agc andm andaf b L S SONIC RHO;global y;y=0 20. 18 0 18 0 20 52.38;%馬赫數(shù) ma%攻角 alpha%三個(gè)時(shí)間 b(1)為導(dǎo)彈離軌時(shí)間 b(2)為起飛發(fā)動(dòng)機(jī)工作結(jié)束時(shí)間 b(3)為續(xù)航發(fā)動(dòng)機(jī)工作結(jié)束時(shí)間b=0 2.1126 44.0832;%系數(shù)表維數(shù)n1=9;n2=6;%andm最小,最大值andm(1)=0.1;andm(2)=0.9;%andaf最小,最大值andaf(1)=0;andaf(2)=10;%阻力系數(shù)acx=.4177 .4404 .5219 .660

16、3 .8534 1.1023; .3858 .4086 .4903 .6290 .8226 1.0723; .3779 .4007 .4827 .6218 .8160 1.0666; .3785 .4015 .4838 .6234 .8184 1.07; .3787 .4018 .4846 .6249 .8209 1.0738; .3829 .4062 .4897 .6310 .8248 1.0835; .3855 .4091 .4934 .6363 .8358 1.0938; .4082 .4321 .5175 .6621 .8641 1.1254; .4947 .5192 .6073 .7

17、571 .9672 1.2392;%升力系數(shù)acy=.0000 .6430 1.4758 2.2870 3.0713 3.8463; .0000 .6454 1.4807 2.2942 3.0915 3.8731; .0000 .6480 1.4858 2.3014 3.0915 3.8731; .0000 .6512 1.4923 2.3107 3.1039 3.8891; .0000 .6554 1.5007 2.3227 3.1197 3.9092; .0000 .6617 1.5134 2.3409 3.1426 3.9401; .0000 .6698 1.5304 2.3661 3.

18、1775 3.9835; .0000 .6792 1.5501 2.3950 3.2162 4.0323; .0000 .6933 1.5935 2.4706 3.3273 4.1790;%推力ap=.000 .15 .49 2.11 2.27 3.53 8.78 25.45 42.80 43.68 44.08; 331.2 614.3 505.4 607.8 48.65 43.97 42.01 41.00 40.80 40.79 2.22;%發(fā)動(dòng)機(jī)質(zhì)量秒流量amc=0. 2.1 2.105 44.1 44.105 100; 2.362 2.362 0.21059 0.21059 0.0000

19、 0.0000;%轉(zhuǎn)動(dòng)慣量ajz=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; 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;%導(dǎo)彈重心axg=.0 2.0 2.4 10.0 18.0 26.0 32.0 38.0 42.0 44.0; .9381 .9095 .9091 .9026 .8969 .8928 08907 .8896 .8895 .8896;agc=.9381 .8896;%靜穩(wěn)定力矩系數(shù)amzaf=0

20、.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

21、 -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;%阻尼力矩導(dǎo)數(shù)%當(dāng)xg=.9381時(shí)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;

22、-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;%當(dāng)xg=.8

23、896時(shí)%其他參數(shù)S=0.0227; L=1.8; SONIC=343.13; RHO=1.225;%四階龍格-庫(kù)塔法子函數(shù)function rk_4(n,h)global y dy;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);for i=1:n old_y(i)=y(i);end for j=1:3 for i=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 for

24、i=1:n y(i)=y(i)+a(1)*dy(i)/3; end%右端子函數(shù)function dery(y)global dy ;global L S SONIC RHO;global ma abs_alpha cx cy mzaf mzwz jz alpha p mc;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);a

25、a(4)=cos(alpha*pi/180);if alpha0 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;%轉(zhuǎn)換成度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;%插值子函數(shù)function i

26、nterp()global acx acy ajz amzaf amzwz axg ap amc agc andm andaf b L ;global ma abs_alpha cx cy mzaf mzwz xg jz p mc;global y;if y(1)b(3) %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);else xg=axg(2,10); jz=ajz(2,14); p=ap(2,11);

27、 mc=amc(2,6);endcx=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);if abs_alpha=0. mzaf=mzaf0+cy*(xg-agc(1)/(abs_alpha*L);%去掉了*RADelse mzaf=mzaf0;endmzwz=interp33(ma,abs_alpha,9,6,andm,andaf,amzwz);%不等距單變?cè)€(xiàn)性插值子函

28、數(shù)function res=interp11(yy,n,x)for j=1:(n-1) if x=yy(1,j+1) i=j;break; else i=n-1; endend res=yy(2,i)+(yy(2,i+1)-yy(2,i)*(x-yy(1,i)/(yy(1,i+1)-yy(1,i);%等距雙變?cè)獟佄锞€(xiàn)線(xiàn)性插值子函數(shù)function res=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=n2-1;end yy1=interp31(x,n1,i,a,yy);yy2=interp31(x,n1,i+1,a,yy);res=(qq-(i-1)*h-bb(1)*yy2-(qq-i*h-bb(1)*yy1)/h;%等距單變?cè)獟佄锞€(xiàn)插值子函數(shù)function res=interp31(x,n1,i,