付録C 安定微係数推算プログラム



In[1]:= (* 飛行速度[m/s] *)

In[2]:= u=7.5;

In[3]:= (* 飛行高度[m] *)

In[4]:= h0=4.5;

In[5]:= (* 空気密度[kg/m^3] *)

In[6]:= rho=0.119*9.8;

In[7]:= (* 全機有害抵抗係数 *)

In[8]:= cd0=0.02134;

In[9]:= (* 飛行機効率 *)

In[10]:= ef=0.8;

In[11]:= (* 機体質量[kg] *)

In[12]:= m=100;

In[13]:= (* 2次元翼の揚力傾斜[rad^-1] *)

In[14]:= a0=6.088; (* DAE11 *)

In[15]:= (* 重力加速度[m/s^2] *)

In[16]:= g=9.8;

In[17]:= (* 重心位置(mac%) *)

In[18]:= h=0.33;

In[19]:= (* 空力中心位置(mac%) *)

In[20]:= hnw=0.25;

In[21]:= (* 慣性モーメント[kg m^2] *)

In[22]:= ixx=1729;

In[23]:= iyy=117;

In[24]:= izz=1762;

In[25]:= ixz=15;

In[26]:= (********************* 主翼 *********************)

In[27]:= (* スパン [m] *)

In[28]:= bw=29;

In[29]:= (* 平面形 *)

In[30]:=cw[y_]:=Module[{t},t=Abs[y];If[t<5.7,1.102,If[t<9.9,1.102+(t-5.7)*(0.888-1.102)/4.2, 0.888+(t-9.9)*(0.651-0.888)/4.6]]]

In[31]:= (* 面積[m^2] *)

In[32]:= sw=NIntegrate[cw[y],{y,-bw/2,bw/2}]

Out[32]= 28.0001

In[33]:= (* 平均空力翼弦[m] *)

In[34]:= macw=NIntegrate[cw[y]^2,{y,0,bw/2}]*2/sw

Out[34]= 0.988678

In[35]:= (* 平均空力翼弦のスパン方向の位置[m] *)

In[36]:=macwy=Module[{t},For[t=0,t<=bw/2,t+=0.001,If[cw[t]-macw<0,Break[]]];t]

Out[36]= 7.925

In[37]:= (* アスベクト比 *)

In[38]:= arw=bw^2/sw

Out[38]= 30.0356

In[39]:= (* 地面効果を考慮した主翼アスベクト比 *)

In[40]:= areh=(1+33*(h0/bw)^1.5)/(33*(h0/bw)^1.5)*arw

Out[40]= 44.9257

In[41]:= (* 主翼揚力傾斜[rad^-1] *)

In[42]:= aw=a0/(1+a0/(Pi*areh))

Out[42]= 5.83625

In[43]:= (* 上半角[rad] *)

In[44]:= dihedral[y_]:=If[y<7.8,0.01923,0.33734];

In[45]:= (* 抵抗係数 *)

In[46]:= cdw=0.0095;

In[47]:= (******************* 水平尾翼 *******************)

In[48]:= (* スパン[m] *)

In[49]:= bh=2.9;

In[50]:= (* 平面形 *)

In[51]:= ch[y_]:=0.744+Abs[y]*(0.558-0.744)/1.45

In[52]:= (* 面積[m^2] *)

In[53]:= sh=NIntegrate[ch[y],{y,-bh/2,bh/2}]

Out[53]= 1.8879

In[54]:= (* 重心から測った位置[m] *)

In[55]:= lh=4.4;

In[56]:= (* 水平尾翼容積 *)

In[57]:= vhr=lh*sh/(macw*sw)

Out[57]= 0.300066

In[58]:= (* 平均空力翼弦[m] *)

In[59]:= mach=NIntegrate[ch[y]^2,{y,0,bh/2}]*2/sh

Out[59]= 0.655429

In[60]:= (* アスベクト比 *)

In[61]:= arh=bh^2/sh

Out[61]= 4.45469

In[62]:= (* 水平尾翼揚力傾斜[rad^-1] *)

In[63]:= ah=a0/(1+a0/(Pi*arh))

Out[63]= 4.24245

In[64]:= (******************* 垂直尾翼 *******************)

In[65]:= (* スパン[m] *)

In[66]:= bv=2.3;

In[67]:= (* 平面形 *)

In[68]:=cv[y_]:=If[y<0,0.8+y*(0.8-0.65)/0.65,0.8+y*(0.49-0.8)/1.65]

In[69]:= (* 面積[m^2] *)

In[70]:= sv=NIntegrate[cv[y],{y,-0.65,1.65}]

Out[70]= 1.5355

In[71]:= (* 重心から測った位置[m] *)

In[72]:= lv=5.3;

In[73]:= (* 垂直尾翼容積 *)

In[74]:= vvr=lv*sv/(bw*sw)

Out[74]= 0.0100223

In[75]:= (* 平均空力翼弦[m] *)

In[76]:= macv=NIntegrate[cv[y]^2,{y,-0.65,1.65}]/sv

Out[76]= 0.678951

In[77]:= (* アスベクト比 *)

In[78]:= arv=bv^2/sv

Out[78]= 3.44513

In[79]:= (* 揚力傾斜[rad^-1] *)

In[80]:= av=a0/(1+a0/(Pi*arv))

Out[80]= 3.89633

In[81]:= (* 力の作用点の位置[m] *)

In[82]:= zvr=0.5;

In[83]:= zvd=0.5;

In[84]:= (**************************************************)

In[85]:= (* 全機揚力係数 *)

In[86]:= cl=2*m*g/(rho*u^2*sw)

Out[86]= 1.06709

In[87]:= (* 全機抵抗係数 *)

In[88]:= cd=cd0+cl^2/(Pi*ef*areh)

Out[88]= 0.0314248

In[89]:= (* 必要推力[N] *)

In[90]:= nt=rho*u^2*sw*cd/2

Out[90]= 28.8601

In[91]:= (* ダウンウォッシュ特性 *)

In[92]:= deda=2*aw/(Pi*areh)

Out[92]= 0.0827026

In[93]:= (* CLα *)

In[94]:= cla=aw*(1+ah/aw*sh/sw*(1-deda))

Out[94]= 6.09864

In[95]:= (* CDα *)

In[96]:= cda=2*cl*cla/(Pi*ef*areh)

Out[96]= 0.115273

In[97]:= (* 胴体容積[m^3] *)

In[98]:= vf=1;

In[99]:= (* 胴体容積比 *)

In[100]:= vfr=vf/(macw*sw)

Out[100]= 0.0361231

In[101]:= (* サイドウォッシュ特性 *)

In[102]:= dsdb=0;

In[103]:= dsdr=0;

In[104]:= (******************* 縦の安定微係数 *****************)

In[105]:= (* 無次元 *)

In[106]:= cxu=-2*nt/(rho*sw*u^2)-2*cd

Out[106]= -0.0942743

In[107]:= czu=0;

In[108]:= cmu=0;

In[109]:= cxa=cl*(1-2*cla/(Pi*ef*areh))

Out[109]= 0.951815

In[110]:= cza=-cla

Out[110]= -6.09864

In[111]:= cma=aw*((h-hnw)-vhr*ah/aw*(1-deda)+vfr*2/aw)

Out[111]= -0.628587

In[112]:= cma'=-2*vhr*lh/macw*ah*deda

Out[112]= -0.937088

In[113]:= cxq=0;

In[114]:= czq=-2*vhr*ah

Out[114]= -2.54603

In[115]:= cmq=-2*vhr*lh/macw*ah

Out[115]= -11.3308

In[116]:= czde=-sh/sw*ah*1

Out[116]= -0.286046

In[117]:= cmde=-vhr*ah*1

Out[117]= -1.27302

In[118]:= (* 有次元 *)

In[119]:= xu=rho*u*sw/(2*m)*cxu

Out[119]= -0.11544

In[120]:= zu=rho*u*sw/(2*m)*(czu-2*cl)

Out[120]= -2.61333

In[121]:= mu=0;

In[122]:= xa=rho*u^2*sw/(2*m)*cxa

Out[122]= 8.74135

In[123]:= za=rho*u^2*sw/(2*m)*cza

Out[123]= -56.0091

In[124]:= ma=rho*u^2*sw*macw/(2*iyy)*cma

Out[124]= -4.87821

In[125]:= ma'=rho*u*sw*macw^2/(4*iyy)*cma'

Out[125]= -0.479335

In[126]:= xq=0;

In[127]:= zq=rho*u*sw*macw/(4*m)*czq

Out[127]= -1.54118

In[128]:= mq=rho*u*sw*macw^2/(4*iyy)*cmq

Out[128]= -5.79589

In[129]:= zde=rho*u^2*sw/(2*m)*czde

Out[129]= -2.62701

In[130]:= mde=rho*u^2*sw*macw/(2*iyy)*cmde

Out[130]= -9.87935

In[131]:= (**************** 横の安定微係数 ********************)

In[132]:= (* 無次元 *)

In[133]:= cyb=-sv/sw*av*(1-dsdb)

Out[133]= -0.213671

In[134]:=clb=-2*a0/(sw*bw)*NIntegrate[dihedral[y]*cw[y]*y,{y,0,bw/2}]

Out[134]= -0.314195

In[135]:= cnb=vvr*av*(1-dsdb)-2*vfr*macw/bw

Out[135]= 0.0365872

In[136]:= cyp=0;

In[137]:= clp=-4*a0/(sw*bw^2)*NIntegrate[cw[y]*y^2,{y,0,bw/2}]

Out[137]= -0.875301

In[138]:=cnp=-4/(sw*bw^2)*(cl-cda)*NIntegrate[cw[y]*y^2,{y,0,bw/2}]

Out[138]= -0.136847

In[139]:= cyr=sv/sw*av*(2*lv/bw+dsdr)

Out[139]= 0.0781004

In[140]:=clr=8/(sw*bw^2)*cl*NIntegrate[cw[y]*y^2,{y,0,bw/2}]+zvr*sv/(bw*sw)*av*(2*lv/bw+dsdr)

Out[140]= 0.308187

In[141]:=cnr=-8/(sw*bw^2)*cdw*NIntegrate[cw[y]*y^2,{y,0,bw/2}]-vvr*av*(2*lv/bw+dsdr)

Out[141]= -0.0170052

In[142]:= clda=0;

In[143]:= cnda=0;

In[144]:= cydr=sv/sw*av

Out[144]= 0.213671

In[145]:= cyda=0;

In[146]:= cldr=zvd*sv/(bw*sw)*av

Out[146]= 0.00368398

In[147]:= cndr=-vvr*av

Out[147]= -0.0390502

In[148]:= (* 有次元 *)

In[149]:= yb=rho*u^2*sw/(2*m)*cyb

Out[149]= -1.96233

In[150]:= lb=rho*u^2*sw*bw/(2*ixx)*clb

Out[150]= -4.83981

In[151]:= nb=rho*u^2*sw*bw/(2*izz)*cnb

Out[151]= 0.553027

In[152]:= yp=0;

In[153]:= lp=rho*u*sw*bw^2/(4*ixx)*clp

Out[153]= -26.0671

In[154]:= np=rho*u*sw*bw^2/(4*izz)*cnp

Out[154]= -3.99908

In[155]:= yr=rho*u*sw*bw/(4*m)*cyr

Out[155]= 1.38671

In[156]:= lr=rho*u*sw*bw^2/(4*ixx)*clr

Out[156]= 9.17805

In[157]:= nr=rho*u*sw*bw^2/(4*izz)*cnr

Out[157]= -0.496944

In[158]:= xdt=0;

In[159]:= yda=0;

In[160]:= lda=0;

In[161]:= nda=0;

In[162]:= ydr=rho*u^2*sw/(2*m)*cydr

Out[162]= 1.96233

In[163]:= ldr=rho*u^2*sw*bw/(2*ixx)*cldr

Out[163]= 0.0567474

In[164]:= ndr=rho*u^2*sw*bw/(2*izz)*cndr

Out[164]= -0.590257

In[165]:= (* primed derivative *)

In[166]:= lbp=(lb+(ixz/ixx)*nb)/(1-ixz^2/(ixx*izz))

Out[166]= -4.83537

In[167]:= lpp=(lp+(ixz/ixx)*np)/(1-ixz^2/(ixx*izz))

Out[167]= -26.1037

In[168]:= lrp=(lr+(ixz/ixx)*nr)/(1-ixz^2/(ixx*izz))

Out[168]= 9.17441

In[169]:= ldap=(lda+(ixz/ixx)*nda)/(1-ixz^2/(ixx*izz))

Out[169]= 0

In[170]:= ldrp=(ldr+(ixz/ixx)*ndr)/(1-ixz^2/(ixx*izz))

Out[170]= 0.0516305

In[171]:= nbp=(nb+(ixz/izz)*lb)/(1-ixz^2/(ixx*izz))

Out[171]= 0.511864

In[172]:= npp=(np+(ixz/izz)*lp)/(1-ixz^2/(ixx*izz))

Out[172]= -4.2213

In[173]:= nrp=(nr+(ixz/izz)*lr)/(1-ixz^2/(ixx*izz))

Out[173]= -0.418842

In[174]:= ndap=(nda+(ixz/izz)*lda)/(1-ixz^2/(ixx*izz))

Out[174]= 0

In[175]:= ndrp=(ndr+(ixz/izz)*ldr)/(1-ixz^2/(ixx*izz))

Out[175]= -0.589818

前の章 目次に戻る 次の章