8. 「とんでる力学」2005年 2 月分

8.1. プログラムと解説

赤木: というわけで、プログラムとかの解説ね。この章が終わればあ と一つだからさくさくいきましょ。

学生S: はあ、、、(プログラム書いたのも中身の説明するのも僕なの に何がさくさくだよこのやろう、とか思ってるかどうかは内緒)

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

学生S: これですね:

 # february.rb

 require 'vector.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
 def cds(alpha)
   alpha = fabs(alpha) % PI
   alpha = PI-alpha if alpha > PI/2
   x = alpha*180/PI
   0.01+x*x*0.0005
 end

 def alphaexact(beta)
   alpha=-beta
   100.times{
     x = alpha+beta
     direction=atan2(cl(x),cd(x))
     alpha -= (direction-alpha)*0.003
     print( alpha, direction, "\n")

   }
   alpha
 end

 class Airplane

   attr_accessor :m, :sw, :ss, :lw, :ls, :moi,
     :x, :v, :omega, :psi, :beta, :gamma

   def theta
     atan2(@v[1],@v[0])
   end

   def alphaw
     @psi - theta
   end

   def alphas
     @beta + alphaw - @gamma*cl(alphaw)+@omega*@ls/@v.abs
   end

   def fwl
     0.5* $rho*v.sqr*@sw*cl(alphaw)
   end
   def fwd
     0.5* $rho*v.sqr*@sw*cd(alphaw)
   end
   def fsl
     0.5* $rho*v.sqr*@ss*cl(alphas)
   end
   def fsd
     0.5* $rho*v.sqr*@ss*cds(alphas)
   end

   def ev
     @v*(1.0/@v.abs)
   end

   def nv
     e=ev
     [-ev[1],ev[0]].to_v
   end

   def fw
     nv*fwl-ev*fwd
   end

   def fs
     nv*fsl-ev*fsd
   end

   def fg
     [0,-1].to_v* $g*@m
   end
   def vdot
     (fg+fw+fs)*(1.0/@m)
   end

   def lwvector
     [cos(@psi),sin(@psi)].to_v*@lw
   end
   def lsvector
     [-cos(@psi),-sin(@psi)].to_v*@ls
   end
   def omegadot
     ((lwvector^fw) + (lsvector^fs))*(1.0/@moi)
   end
 end

 $rho=1
 $g=10

 def solve(dt,tend,glider)
   t = 0

   ta=[]
   xa=[]
   ya=[]
   vxa=[]
   vya=[]
   psia=[]
   alphaa=[]
   omegaa=[]
   while t < tend
     gtmp = glider.clone
     h = dt*0.5
     gtmp.x     += glider.v*h
     gtmp.v     += glider.vdot*h
     gtmp.psi   += glider.omega*h
     gtmp.omega += glider.omegadot*h

     glider.x     += gtmp.v*dt
     glider.v     += gtmp.vdot*dt
     glider.psi   += gtmp.omega*dt
     glider.omega += gtmp.omegadot*dt

     t += dt

     #  print t, " ", glider.x.join(" ")," ",  glider.v.join(" ")," ",
     #    glider.psi," ",  glider.omega,"\n"
     ta.push(t)
     xa.push(glider.x[0])
     ya.push(glider.x[1])
     vxa.push(glider.v[0])
     vya.push(glider.v[1])
     psia.push(glider.psi)
     omegaa.push(glider.omega)
     alphaa.push(glider.alphaw)
   end
  $result=[ta,xa,ya,vxa,vya,psia,omegaa,alphaa]

 end

 glider = Airplane.new
 glider.m=0.02
 glider.sw = 0.04
 glider.ss = glider.sw*0.1
 glider.lw = 0.0
 glider.ls = 0.3
 glider.moi = 2e-4
 glider.beta= -0.05
 glider.gamma= 0.1

 glider.x = [0,0].to_v
 glider.v = [7,0].to_v
 glider.psi = 0
 glider.omega=0

 default_glider= glider.clone

 unless File.exist?("february-fig2.ps")
   glider=default_glider.clone
   printer february-fig2.ps/vcps
   square
   viewport 0.1 1.1 0.8 1
   limit 0 70 -6 0
   expand 1.5
   box
   relocate 35 -8.5
   label X
   relocate -9 -3
   label Y
   solve(1.0/256, 10,glider)
   xdata $result[1]
   ydata $result[2]
    conn

   viewport 0 1 -1.5 -0.5
   limit 0 10 -0.1 0.1
   box
   relocate 5 -0.18
   label T
   relocate -1.8 0
   label \\gq,\\ga
   xdata $result[0]
   ydata $result[7]
   conn
   ydata $result[5]
   conn

   pgend
   psfix
 end

 unless File.exist?("february-fig2b.ps")
   glider=default_glider.clone
   printer february-fig2b.ps/vcps

   def cd(alpha)
     0
   end
   def cds(alpha)
     0
   end

   square
   viewport 0.1 1.1 0.8 1
   limit 0 70 -2 2
  expand 1.5
   box
   relocate 35 -4
   label X
   relocate -9 1
   label Y
   solve(1.0/256, 10,glider)
   xdata $result[1]
   ydata $result[2]
    conn

   viewport 0 1 -1.5 -0.5
   limit 0 10 -0.05 0.15
   box
   relocate 5 -0.13
   label T
   relocate -1.9 0.1
   label \\gq,\\ga
   xdata $result[0]
   ydata $result[7]
   conn
   ydata $result[5]
   conn

   pgend
   psfix
 end

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

学生S: クラスの説明とかしてもあれですよね？とりあえずグラフ書く ところからにします。まず

 glider = Airplane.new
 glider.m=0.02
 glider.sw = 0.04
 glider.ss = glider.sw*0.1
 glider.lw = 0.0
 glider.ls = 0.3
 glider.moi = 2e-4
 glider.beta= -0.05
 glider.gamma= 0.1

 glider.x = [0,0].to_v
 glider.v = [7,0].to_v
 glider.psi = 0
 glider.omega=0

 default_glider= glider.clone

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

 unless File.exist?("february-fig2.ps")
   glider=default_glider.clone
   printer february-fig2.ps/vcps
   square
   viewport 0.1 1.1 0.8 1
   limit 0 70 -6 0
   expand 1.5
   box
   relocate 35 -8.5
   label X
   relocate -9 -3
   label Y
   solve(1.0/256, 10,glider)
   xdata $result[1]
   ydata $result[2]
    conn

   viewport 0 1 -1.5 -0.5
   limit 0 10 -0.1 0.1
   box
   relocate 5 -0.18
   label T
   relocate -1.8 0
   label \\gq,\\ga
   xdata $result[0]
   ydata $result[7]
   conn
   ydata $result[5]
   conn

   pgend
   psfix
 end

solve 以外は全部グラフの枠とかプロットとかしてるだけなので、 solve の ほうをみてみます。

 def solve(dt,tend,glider)
   t = 0
   ta=[]
   xa=[]
   ya=[]
   vxa=[]
   vya=[]
   psia=[]
   alphaa=[]
   omegaa=[]
   while t < tend
     gtmp = glider.clone
     h = dt*0.5
     gtmp.x     += glider.v*h
     gtmp.v     += glider.vdot*h
     gtmp.psi   += glider.omega*h
     gtmp.omega += glider.omegadot*h

     glider.x     += gtmp.v*dt
     glider.v     += gtmp.vdot*dt
     glider.psi   += gtmp.omega*dt
     glider.omega += gtmp.omegadot*dt
     t += dt
     ta.push(t)
     xa.push(glider.x[0])
     ya.push(glider.x[1])
     vxa.push(glider.v[0])
     vya.push(glider.v[1])
     psia.push(glider.psi)
     omegaa.push(glider.omega)
     alphaa.push(glider.alphaw)
   end
  $result=[ta,xa,ya,vxa,vya,psia,omegaa,alphaa]
 end

このプログラムでは、微分方程式を解きながら線を引くのではなくて、解いた 答を配列に全部しまってそれを呼び出し側に返す、というふうにしてみました。 何故そうしたかは良く憶えてないですが、多分、線で引こうとするとプログラ ムがちょっと面倒になると思ったからです。線を引くには前の値をとっておい て、 2 点を結ぶ必要があるのですが、パネル何枚も使うとかするとその辺が ややこしいですから。

データファイルに落としてもいいのですが、そんなに計算時間がかかるわけで もないし、データの大きさもしれてるので配列にそのままいれました。配列に しまってるのは

     xa.push(glider.x[0])

とかのところです。数値積分は2次のルンゲクッタで

の形の公式を使いました。 ではなくて の値を gtmp という Airplane クラスのインスタンスにしまうようにしてます。

赤木: なるほど。 2 次にしたのはなんか理由があるの？

学生S: 1次だとさすがに精度がアレだし、4次とかは書くのがちょっと面倒になるので。

赤木: まあ、そんなに長い積分するわけでもないから、これくらいで 十分よね。微分方程式の右辺の計算は？

学生S: これはまあ式の通りなので、説明することもないです。

赤木: んなこといっても読者にはわかんないわよ。 vdot は？

学生S: えーと

  def vdot
    (fg+fw+fs)*(1.0/@m)
  end

ですね。 fg は重力で

  def fg
    [0,-1].to_v* $g*@m
  end

だから、下向きの というだけです。 fw は

  def fw
    nv*fwl-ev*fwd
  end

で、 nv, fwl は

  def nv
    e=ev
    [-ev[1],ev[0]].to_v
  end

あれ、 e って使ってないですね。 ev は速度方向の単位ベクトルで

  def ev
    @v*(1.0/@v.abs)
  end

fwl に戻ると

  def fwl
    0.5* $rho*v.sqr*@sw*cl(alphaw)
  end

で、 で、揚力です。抗力のほうも同じようなものです。

赤木: 関数が大体式1つに対応してるのね。まあ、わかりやすいかも。 でも、なんか無駄が多くない？ cl とか cd とか何回も何回も計算されるし。

学生S: 無駄は多いと思います。とりあえず正しく動くプログラムを早 く作るほうを優先にして速度はまあ、、、というつもりで作ったので。色々す れば 10 倍くらいは速くできるかもしれませんが、まあ、見てるうちに計算終 わるので。