Previous ToC Next

9. 「とんでる力学」2005年 3 月分

9.1. プログラムと解説

赤木: というわけで、プログラムとかの解説ね。この章で最後ね。

学生S: ですね。

赤木: とりあえずリスト出して。

学生S: これですね:


 # march.rb

 require 'vector3.rb'
 include Math
 def fabs(x)
  x>0 ? x : -x
 end

 def cl(alpha)
   stallangle=15
   maxangle=45
   clcoef = 0.05
   alpha = atan2(sin(alpha),cos(alpha))
   x = alpha*180/PI
   xabs = fabs(x)
   if xabs <= stallangle
     y=xabs*clcoef
   elsif xabs < maxangle
     y=stallangle*clcoef*0.4
   else
     y=stallangle*clcoef*0.4*(90-x)/(90-maxangle)
   end
   x> 0 ? y : -y
 end
 def cd(alpha)
   alpha = fabs(alpha) % PI
   alpha = PI-alpha if alpha > PI/2
   x = alpha*180/PI
   0.01+x*x*0.0005
 end

 class Airplane
  x, v, u is always on ground coordinate

  all else are in plane coordinate
   attr_accessor :m, :sw, :ss, :sv, :lw, :ls, :beta,
     :x, :v, :omega, :u, :span, :i0, :dihedral

   def to_plane(x,u)
     y=[]
     x.each_index{|k| y[k]=x*u[k]}
     y.to_v
   end
   def to_ground(x,u)
     y=[]
     ut=transpose(u)
     x.each_index{|k| y[k]=x*ut[k]}
     y.to_v
   end

   def normal_vector(side)
     if side== "left"
       nl=[0,-sin(@dihedral),cos(@dihedral)]
     elsif  side== "right"
       nl=[0, sin(@dihedral),cos(@dihedral)]
     elsif  side== "horizontal"
       nl=[sin(@beta),0,cos(@beta)]
     elsif  side== "vertical"
       nl=[0,1,0]
     else
       raise "normal_vector: unknown parameter #{side}"
     end
     nl.to_v
   end
   def moment_arm(side)
     if side== "left"
       ls=[0,@span*0.25,0]
     elsif  side== "right"
       ls=[0,-@span*0.25,0]
     elsif  side== "horizontal"
       ls=[-@ls,0,0]
     elsif  side== "vertical"
       ls=[-@ls,0,0]
     else
       raise "moment_arm: unknown parameter #{side}"
     end
     ls.to_v
   end

   def area(side)
     if (side== "left") or (side== "right")
       x=@sw/2
     elsif  side== "horizontal"
       x=@ss
     elsif  side== "vertical"
       x=@sv
     else
       raise "area: unknown parameter #{side}"
     end
     x
   end
   def alphas(side,v,u,omega)
     nl=normal_vector(side)
     ls=moment_arm(side)
     vrot = omega^ls
     vrot *= 1.5 if ((side == "right") or (side == "left") )
     vreal=to_plane(v,u)+vrot
     -asin(vreal.normalize*nl)
   end

  force on a wing in plane coordinate
   def fs(side,v,u,omega)
     alpha=alphas(side,v,u,omega)
     nl=normal_vector(side)
     ls=moment_arm(side)
     vrot = omega^ls
     vrot *= 1.5 if ((side == "right") or (side == "left") )
     vreal=to_plane(v,u)+vrot
     fl= 0.5* $rho*vreal.sqr*area(side)*cl(alpha)
     fd= 0.5* $rho*vreal.sqr*area(side)*cd(alpha)
     evp=vreal.normalize
     el=normal_vector(side)-evp*(normal_vector(side)*evp)
     el*fl-evp*fd
   end

   def fg
     [0,0,-1].to_v* $g*@m
   end
   def vdot
     fair = fs("left",@v,@u,@omega)+
       fs("right",@v,@u,@omega)+
       fs("horizontal",@v,@u,@omega)+
       fs("vertical",@v,@u,@omega)
     (fg+to_ground(fair,@u))*(1.0/@m)
   end
   def xdot
     @v
   end

   def omegadot
     torque = @omega^([@i0[0]*@omega,@i0[1]*@omega,@i0[2]*@omega].to_v)
     ["left","right","horizontal","vertical"].each{|x|
       torque += fs(x,@v,@u,@omega)^moment_arm(x)
     }
     torque.each_index{|i| torque[i]= -torque[i]/@i0[i][i]}
     torque.to_v
   end
   def udot
     omegag = to_ground(omega,@u)
     [omegag^u[0],omegag^u[1],omegag^u[2]]
   end

   def pp
     print("x =",@x.join(" "),"\n")
     print("v =",@v.join(" "),"\n")
     print("om=",@omega.join(" "),"\n")
     print("do=",omegadot.join(" "),"\n")
     3.times{|k|
       print("u[#{k}]=",@u[k].join(" "),"\n")
     }
     3.times{|k|
       print("udot[#{k}]=",udot[k].join(" "),"\n")
     }
     print("vinp=",to_plane((@v.normalize),@u).join(" "),"\n")
     print("uy*v=",@u[1]*(@v.normalize),"\n")
     print("alpawl=", alphas("left",@v,@u,@omega),
         " alpawr=", alphas("right",@v,@u,@omega),"\n")
     print("alpas=", alphas("horizontal",@v,@u,@omega),"\n")
     print("alpav=", alphas("vertical",@v,@u,@omega),"\n")
     print("Fwl=", fs("left",@v,@u,@omega).join(" "),"\n")
     print("Fwr=", fs("right",@v,@u,@omega).join(" "),"\n")
   end

   def pp2
     print("x =",@x.join(" "),"\n")
     print("v =",@v.join(" "),"\n")
     print("u[1]=",@u[1].join(" "),"\n")
   end
 end
 $rho=1
 $g=10

 def solve(dt,tend,glider)
   t = 0
   $ta=[]
   $xa=[]
   $va=[]
   $ua=[]
   $omegaa=[]
   while t < tend
     gtmp = glider.clone
     gtmp.u = to_rot_mat(gtmp.u)
     h = dt*0.5
     gtmp.x     += glider.v*h
     gtmp.v     += glider.vdot*h
     ud = glider.udot
     gtmp.u[0] += ud[0]*h
     gtmp.u[1] += ud[1]*h
     gtmp.u = to_rot_mat(gtmp.u)
     gtmp.omega += glider.omegadot*h
     h=dt

     gt2=gtmp
     gtmp=glider
     glider=gt2
     gtmp.x     += glider.v*h
     gtmp.v     += glider.vdot*h
     ud = glider.udot
     gtmp.u[0] += ud[0]*h
     gtmp.u[1] += ud[1]*h
     gtmp.u = to_rot_mat(gtmp.u)
     gtmp.omega += glider.omegadot*h

     glider=gtmp
     t += dt

     print("t=",t,"\n")
     glider.pp2
     $ta.push(t)
     $xa.push(glider.x)
     $va.push(glider.v)
     $ua.push(glider.u)
     $omegaa.push(glider.omega)
   end
 end

 glider = Airplane.new
 glider.m=0.02
 glider.span = 0.5
 glider.sw = 0.04
 glider.ss = glider.sw*0.1
 glider.sv = glider.ss*0.2
 glider.lw = 0.0
 glider.ls = 0.3
 glider.dihedral = 0.15
 glider.i0=[[2e-4,0,0].to_v,[0,2e-4,0].to_v,[0,0,2e-4].to_v]
 glider.beta= 0.08
 glider.x = [0,0,0].to_v
 glider.v = [6.583,0.0,-0.621].to_v
 glider.omega=[0,0,0].to_v
 theta=0.08
 glider.u=to_rot_mat([[1,0,-0.0143].to_v,[0,1,0.05].to_v,[0,0,0].to_v])

 default_glider= glider.clone
 unless File.exist?("march-fig3.ps")
   glider=default_glider.clone
   printer march-fig3.ps/vcps
   square
   viewport 0.1 1.1 0.7 1
   limit 0 100 -29 1
   expand 1.5
   lw 2
   box BCTSV BCNTSV
   relocate -15 -3
   label Y
   glider=default_glider.clone
   solve(1.0/256, 15,glider)
   xdata $xa.collect{|x|x[0]}
   ydata $xa.collect{|x|x[1]}
   conn
   viewport 0 1 -0.6666666666666666 0
   limit 0 100 -20 0
   expand 1.5
   box BCNTSV  BCNTSV
   relocate 50 -26
   label X
   relocate -15 -3
   label Z
   xdata $xa.collect{|x|x[0]}
   ydata $xa.collect{|x|x[2]}
   conn

   viewport 0 1 -1.5 -0.5
   limit 0 15 -0.1 0.1
   expand 1.5
   box BCNTSV  BCNTSV
   relocate 7.5 -0.2
   label T
   relocate -2.5 0
   label U\\d23
   xdata $ta
   ydata $ua.collect{|x|x[1][2]}
   conn
 end
 unless File.exist?("march-fig4.ps")
   glider=default_glider.clone
   glider.dihedral=0.02
   printer march-fig4.ps/vcps
   square
   viewport 0.1 1.1 0.7 1
   limit 0 100 -29 1
   expand 1.5
   lw 2
   box BCTSV BCNTSV
   relocate -15 -3
   label Y
   solve(1.0/256, 15,glider)
   xdata $xa.collect{|x|x[0]}
   ydata $xa.collect{|x|x[1]}
   conn
   viewport 0 1 -0.6666666666666666 0
   limit 0 100 -20 0
   expand 1.5
   box BCNTSV  BCNTSV
   relocate 50 -26
   label X
   relocate -15 -3
   label Z
   xdata $xa.collect{|x|x[0]}
   ydata $xa.collect{|x|x[2]}
   conn

   viewport 0 1 -1.5 -0.5
   limit 0 15 -0.05 0.15
   expand 1.5
   box BCNTSV  BCNTSV
   relocate 7.5 -0.15
   label T
   relocate -2.5 0
   label U\\d23
   xdata $ta
   ydata $ua.collect{|x|x[1][2]}
   conn
 end
 unless File.exist?("march-fig5.ps")
   glider=default_glider.clone
   glider.sv = glider.ss*0.01
   printer march-fig5.ps/vcps
   square
   viewport 0.1 1.1 0.7 1
   limit 0 100 -20 10
   expand 1.5
   lw 2
   box BCTSV BCNTSV
   relocate -15 -3
   label Y
   solve(1.0/256, 15,glider)
   xdata $xa.collect{|x|x[0]}
   ydata $xa.collect{|x|x[1]}
   conn
   viewport 0 1 -0.6666666666666666 0
   limit 0 100 -20 0
   expand 1.5
   box BCNTSV  BCNTSV
   relocate 50 -26
   label X
   relocate -15 -3
   label Z
   xdata $xa.collect{|x|x[0]}
   ydata $xa.collect{|x|x[2]}
   conn

   viewport 0 1 -1.5 -0.5
   limit 0 15 -0.1 0.1
   expand 1.5
   box BCTSV  BCNTSV
   relocate -2.5 0
   label U\\d23
   xdata $ta
   ydata $ua.collect{|x|x[1][2]}
   conn
   viewport 0 1 -1 0
   limit 0 15 -0.1 0.1
   expand 1.5
   box BCNTSV  BCNTSV
   relocate 7.5 -0.2
   label T
   relocate -2.5 0
   label U\\d12
   xdata $ta
   ydata $ua.collect{|x|x[0][1]}
   conn
 end

赤木: で、どこから話を始める?

学生S: 前回と同じく、グラフ書くところからにします。まず

 glider = Airplane.new
 glider.m=0.02
 glider.span = 0.5
 glider.sw = 0.04
 glider.ss = glider.sw*0.1
 glider.sv = glider.ss*0.2
 glider.lw = 0.0
 glider.ls = 0.3
 glider.dihedral = 0.15
 glider.i0=[[2e-4,0,0].to_v,[0,2e-4,0].to_v,[0,0,2e-4].to_v]
 glider.beta= 0.08

 glider.x = [0,0,0].to_v
 glider.v = [6.583,0.0,-0.621].to_v
 glider.omega=[0,0,0].to_v
 theta=0.08
 glider.u=to_rot_mat([[1,0,-0.0143].to_v,[0,1,0.05].to_v,[0,0,0].to_v])
 default_glider= glider.clone

ここまでは前回と同じく初期設定です。変数名は基本的には本文のと同じで、 それぞれ Airplane クラスの インスタンス変数で、外からアクセスできるよ うにしたものです。で、これを使い回すのでコピーを default_glider にとっ ておきます。以下が本体です。

 unless File.exist?("march-fig3.ps")
   glider=default_glider.clone
   printer march-fig3.ps/vcps
   square
   viewport 0.1 1.1 0.7 1
   limit 0 100 -29 1
   expand 1.5
   lw 2
   box BCTSV BCNTSV
   relocate -15 -3
   label Y
   glider=default_glider.clone
   solve(1.0/256, 15,glider)
   xdata $xa.collect{|x|x[0]}
   ydata $xa.collect{|x|x[1]}
   conn
   viewport 0 1 -0.6666666666666666 0
   limit 0 100 -20 0
   expand 1.5
   box BCNTSV  BCNTSV
   relocate 50 -26
   label X
   relocate -15 -3
   label Z
   xdata $xa.collect{|x|x[0]}
   ydata $xa.collect{|x|x[2]}
   conn

   viewport 0 1 -1.5 -0.5
   limit 0 15 -0.1 0.1
   expand 1.5
   box BCNTSV  BCNTSV
   relocate 7.5 -0.2
   label T
   relocate -2.5 0
   label U\\d23
   xdata $ta
   ydata $ua.collect{|x|x[1][2]}
   conn
 end

solve 以外は全部グラフの枠とかプロットとかしてるだけなのも先月と同じで すが、ちょっと違うのは空間が 3 次元なので結果の位置とかも 3 次元ベクト ルの配列になってることです。 x, y 要素をプロットするのに

   xdata $xa.collect{|x|x[0]}
   ydata $xa.collect{|x|x[2]}

と collect を使って x 座標、 y 座標を切りだしてます。solve のほうにい きます。

 def solve(dt,tend,glider)
   t = 0
   $ta=[]
   $xa=[]
   $va=[]
   $ua=[]
   $omegaa=[]
   while t < tend
     gtmp = glider.clone
     gtmp.u = to_rot_mat(gtmp.u)
     h = dt*0.5
     gtmp.x     += glider.v*h
     gtmp.v     += glider.vdot*h
     ud = glider.udot
     gtmp.u[0] += ud[0]*h
     gtmp.u[1] += ud[1]*h
     gtmp.u = to_rot_mat(gtmp.u)
     gtmp.omega += glider.omegadot*h

     h=dt
     gt2=gtmp
     gtmp=glider
     glider=gt2

     gtmp.x     += glider.v*h
     gtmp.v     += glider.vdot*h
     ud = glider.udot
     gtmp.u[0] += ud[0]*h
     gtmp.u[1] += ud[1]*h
     gtmp.u = to_rot_mat(gtmp.u)
     gtmp.omega += glider.omegadot*h
     glider=gtmp

     t += dt
     print("t=",t,"\n")
     glider.pp2

     $ta.push(t)
     $xa.push(glider.x)
     $va.push(glider.v)
     $ua.push(glider.u)
     $omegaa.push(glider.omega)
   end
 end

2 次のルンゲクッタというのも先月と同じです。 x, v u は速度、位置、それ から方向を示す回転行列で、それらを直接積分するだけですね。omega は角速度ベクトルで、問題は vdot や udot です。vdot は

  def vdot
    fair = fs("left",@v,@u,@omega)+
      fs("right",@v,@u,@omega)+
      fs("horizontal",@v,@u,@omega)+
      fs("vertical",@v,@u,@omega)
    (fg+to_ground(fair,@u))*(1.0/@m)
  end

えーと、なんでしょこれ?

赤木: って、君が書いたわけだし。

学生S: ですね、、、 fs を見ると

  def fs(side,v,u,omega)
    alpha=alphas(side,v,u,omega)
    nl=normal_vector(side)
    ls=moment_arm(side)
    vrot = omega^ls
    vrot *= 1.5 if ((side == "right") or (side == "left") )
    vreal=to_plane(v,u)+vrot
    fl= 0.5* $rho*vreal.sqr*area(side)*cl(alpha)
    fd= 0.5* $rho*vreal.sqr*area(side)*cd(alpha)
    evp=vreal.normalize
    el=normal_vector(side)-evp*(normal_vector(side)*evp)
    el*fl-evp*fd
  end

えーと、まあ、その、これ結局本文の定義通りで、、、、

赤木: omegadot は?

学生S: これもなんか式の通りですね。

  def omegadot
    torque = @omega^([@i0[0]*@omega,@i0[1]*@omega,@i0[2]*@omega].to_v)
    ["left","right","horizontal","vertical"].each{|x|
      torque += fs(x,@v,@u,@omega)^moment_arm(x)
    }
    torque.each_index{|i| torque[i]= -torque[i]/@i0[i][i]}
    torque.to_v
  end

udot も、機体座標系から omega を地面系に戻して外積計算するだけで

  def udot
    omegag = to_ground(omega,@u)
    [omegag^u[0],omegag^u[1],omegag^u[2]]
  end

です。

赤木: これではいくらなんでも説明が簡単すぎよね。もうちょっとな んとかして。

学生S: 今日はちょっと、、、また後でにしませんか?

Previous ToC Next