Previous | ToC | Next |
赤木: というわけで、プログラムとかの解説ね。この章で最後ね。
学生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 |