word格式完美整理/數值實驗作業1實驗背景實驗名稱：鉛垂平面飛行彈道仿真及分析實驗內容與要求：根據描述飛行器在鉛垂平面內運動的數學模型,編制某導彈的鉛垂平面無控飛行彈道仿真程序,利用計算機解算初始段無控飛行彈道,對初始段彈道參數的變化規律進行分析。2建立數學模型：3計算方法研究確定數值積分方法和積分步長使用算法：四階龍格—庫塔法。使用步長：h=0.0054原始數據：1初值x=0<m>y=20.0<m>=18=18 v=20<m/s>z=0<rad/s> m=52.38<kg>2攻角與馬赫數范圍<僅用于插值計算>攻角=0~10馬赫數=0.1~0.93阻力系數表馬赫數攻角<>02468100.1.4177.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 4升力系數表馬赫數攻角<>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.39503.21624.03230.9.0000.69331.59352.47063.32734.17905推力數據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第一級工作結束時間：2.1126s,第二級工作結束時間：44.0832s6發動機質量秒流量t<s>0.2.12.10544.144.105100秒流量<kg/s>2.3622.3620.210590.210590.0.7轉動慣量t<s>.02.02.46.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.698導彈重心〔起自頭部t<s>.02.02.410.018.026.032.038.042.044.0XG<m>.9381.9095.9091.9026.8969.8928.8907.8896.8895.88969靜穩定力矩系數馬赫數攻角<>02468100.10.0000-0.0104-0.0341-0.0564-0.0771-0.09850.20.0000-0.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-0.0252-0.0425-0.0578-0.07390.90.0000-0.0053-0.0229-0.0391-0.0538-0.0693當導彈重心變化時的修正公式：10阻尼力矩導數當Xg=.9381時馬赫數攻角<>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.5090-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當Xg=.8896時馬赫數攻角<>02468100.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.76360.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其它參數特征面積S<m2>特征長度L<m>毛翼展<m>音速SONIC<m/s>大氣密度<kg/m3>0.02271.80.5343.131.225使用的插值算法：氣動數據插值——等距雙變元拋物線插值；推力、重心、轉動慣量等——不等距一元線性插值。5空氣動力和空氣動力矩表達式6編制計算程序計算機算法采用了matlab實現,源程序見附件.該程序有八個函數組成,各函數之間的調用關系如下圖所示.子函數子函數initl主函數main子函數rk_4子函數result子函數savedata子函數drawing子函數dery子函數interp子函數interp33子函數interp31子函數interp11子函數initl的功能是輸入求解導彈運動方程組所需的原始數據.子函數rk_4是四階龍格-庫塔法積分算法子函數,其中調用了子函數dery.子函數dery功能是計算微分方程組的右端函數,其中調用了子函數interp.插值子函數interp功能是計算所有需要插值的參數,其中調用了子函數interp11和interp33.子函數interp11是不等距單變元線性插值函數,主要用于轉動慣量Jz,質心位置xg,推力p等單變量的插值.子函數interp33是等距雙變元拋物線-線性插值函數,主要用于氣動力系數cx,cy,氣動力矩系數,,等雙變元參數的插值,該函數調用了子函數interp31.子函數interp31是等距單變元拋物線插值函數,被interp33調用,完成子函數interp33的等距雙變元拋物線-線性插值功能.關于程序中出現的數組和變量名,作如下說明:acx:阻力系數數組;acy:升力系數數組;amzaf:某一質心位置下的靜穩定性導數數組;amzwz:阻尼力矩系數導數數組;axg:質心位置Xg隨時間的變化規律;ajz:轉動慣量Jz隨時間的變化規律;ap:起飛,續航發動機的推力值;amc:起飛,續航發動機的燃料質量的秒流量值;agc:質心位置變化的始末值;andm:氣動數據插值所需的Ma數的最小,最大值;andaf:氣動數據插值所需的攻角的最小,最大值;y:存放積分結果的數組,該數組在程序開始時存放積分初值;dy:存放右端函數數值的數組;b:存放三個時間值的數組,其中b<2>存放起飛發動機工作結束時間;b<3>存放續航發動機工作結束時間;程序中一些主要的變量名有:L—特征長度;S—特征面積;SONIC—聲速C;RHO—大氣密度;h—積分步長;程序中其他變量都是存放中間結果的變量.7計算結果運行程序得到彈道曲線,速度曲線,攻角曲線如下:通過運用龍格庫塔法和一些插值方法求解導彈運動方程組,獲得了導彈各運動參數的變化規律.通過該算法獲得的變化規律比較接近真實情況.附件:源程序functionmain<>clearallclearglobalglobaly;globalii;ii=0;h=0.005;initl<>;whiley<7>>=0ii=ii+1;result<ii>;rk_4<8,h>;endsavedata<ii>;drawing<>;%原始數據初始化functioninitl<>globalacxacyajzamzafamzwzaxgapamcagcandmandafbLSSONICRHO;globaly;y=[020.1801802052.38];%馬赫數ma%攻角alpha%三個時間b<1>為導彈離軌時間b<2>為起飛發動機工作結束時間b<3>為續航發動機工作結束時間b=[02.112644.0832];%系數表維數n1=9;n2=6;%andm最小,最大值andm<1>=0.1;andm<2>=0.9;%andaf最小,最大值andaf<1>=0;andaf<2>=10;%阻力系數acx=[.4177.4404.5219.6603.85341.1023;.3858.4086.4903.6290.82261.0723;.3779.4007.4827.6218.81601.0666;.3785.4015.4838.6234.81841.07;.3787.4018.4846.6249.82091.0738;.3829.4062.4897.6310.82481.0835;.3855.4091.4934.6363.83581.0938;.4082.4321.5175.6621.86411.1254;.4947.5192.6073.7571.96721.2392];%升力系數acy=[.0000.64301.47582.28703.07133.8463;.0000.64541.48072.29423.09153.8731;.0000.64801.48582.30143.09153.8731;.0000.65121.49232.31073.10393.8891;.0000.65541.50072.32273.11973.9092;.0000.66171.51342.34093.14263.9401;.0000.66981.53042.36613.17753.9835;.0000.67921.55012.39503.21624.0323;.0000.69331.59352.47063.32734.1790];%推力ap=[.000.15.492.112.273.538.7825.4542.8043.6844.08;331.2614.3505.4607.848.6543.9742.0141.0040.8040.792.22];%發動機質量秒流量amc=[0.2.12.10544.144.105100;2.3622.3620.210590.210590.00000.0000];%轉動慣量ajz=[0.2.02.46.410.414.418.422.426.430.434.038.442.444.0;8.357.887.867.817.787.757.737.717.707.707.697.697.697.69];%導彈重心axg=[.02.02.410.018.026.032.038.042.044.0;.9381.9095.9091.9026.8969.892808907.8896.8895.8896];agc=[.9381.8896];%靜穩定力矩系數amzaf=[0.0000-0.0104-0.0341-0.0564-0.0771-0.0985;0.0000-0.0104-0.0341-0.0564-0.0770-0.09830.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-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];%阻尼力矩導數%當xg=.9381時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;-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];%%當xg=.8896時%其他參數S=0.0227;L=1.8;SONIC=343.13;RHO=1.225;%四階龍格-庫塔法子函數functionrk_4<n,h>globalydy;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>;fori=1:nold_y<i>=y<i>;endforj=1:3fori=1:ny1<i>=old_y<i>+a<j>*dy<i>;y<i>=y<i>+a<j+1>*dy<i>/3;enddery<y1>;endfori=1:ny<i>=y<i>+a<1>*dy<i>/3;end%右端子函數functiondery<y>globaldy;globalLSSONICRHO;globalmaabs_alphacxcymzafmzwzjzalphapmc;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>;aa<4>=cos<alpha*pi/180>;ifalpha<0cy=-cy;endxf=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;%轉換成度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;%插值子函數functioninterp<>globalacxacyajzamzafamzwzaxgapamcagcandmandafbL;globalmaabs_alphacxcymzafmzwzxgjzpmc;globaly;ify<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>>;elsexg=axg<2,10>;jz=ajz<2,14>;p=ap<2,11>;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>;ifabs_alpha~=0.mzaf=mzaf0+cy*<xg-agc<1>>/<abs_alpha*L>;%去掉了*RADelsemzaf=mzaf0;endmzwz=interp33<ma,abs_alpha,9,6,andm,andaf,amzwz>;%不等距單變元線性插值子函數functionres=interp11<yy,n,x>forj=1:<n-1>ifx<=yy<1,j+1>i=j;break;elsei=n-1;endendres=yy<2,i>+<yy<2,i+1>-yy<2,i>>*<x-yy<1,i>>/<yy<1,i+1>-yy<1,i>>;%等距雙變元拋物線線性插值子函數functionres=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><0i=1;elseif<i-n2>>=0i=n2-1;endyy1=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;%等距單變元拋物線插值子函數functionres=interp31<x,n1,i,a,yy>h=<a<2>-a<1>>/<n1-1>;i1=fix<<x-a