| 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 |