座標法と数式処理ソフト maxima で断面性能の計算式を求めるやり方を紹介します。
座標法でも 形状寸法を変数とすれば 計算式を 求めることができます。多角形に限られますが 積分をする必要はありません。
数式処理ソフトを利用すれば 多角形の頂点の座標を定義するだけで 断面積から主軸の角度まで算定することができます。いまのところ、円弧や放物線は 折れ線に置き換える必要があるのですが、今後の研究によって 合理的な手法が提案されてくると思われます。
多角形の断面性能の計算式は
となります。座標に数値を代入すれば 断面性能が算定できます。
ただし
座標 = [x1, y1], [x2, y2], [x3, y3], …, [xn, yn]
Ax = S, Ix = Jx, Iy = Jy, Ixy = Jxy
図心回りの2次モーメント
Ixc = Ix - Ax * gy ** 2
Iyc = Iy - Ax * gx ** 2
図心回りの相乗モーメント
Ixyc = Ixy - Ax * gx * gy
Xc軸からの主軸の角度 α
α = atan2(2 * Ixyc, Iyc - Ixc) / 2
Σ は i = 1, n-1
座標は閉鎖形の頂点でそれを結ぶ線が交差しないようにする
※参考文献
○材料力学本論 / 前澤・吉峯訳 / コロナ社 付録2 平面図形の諸元
○建造力学 / 谷口著 / 裳華房 第20章 断面の主軸
さて
[ Maxima/面積計算 座標法 倍横距法 / http://ja.wikibooks.org/wiki/Maxima/%E9%9D%A2%E7%A9%8D%E8%A8%88%E7%AE%97_%E5%BA%A7%E6%A8%99%E6%B3%95_%E5%80%8D%E6%A8%AA%E8%B7%9D%E6%B3%95 ]
上記 ホームページに
多角形の面積計算
座標(1,1),(3,1),(3,3),(1,3)の多角形の面積を求める。maxima のプログラムがあります。
(オリジナル)
/*size:行列の寸法関数 area:面積計算関数*/
size(Z):=block(fn(f):=1,C:matrixmap(fn,col(Z,1)),R:matrixmap(fn,row(Z,1)),
return(matrix([transpose(C).C,R.transpose(R)])))$
area(B):=block(B:addcol(B,col(B,1)),B:addcol(B,col(B,2)),
return(sum(B[1,i]*(B[2,i+1]-B[2,i-1]),i,(1+1),(size(B)[1,2]-2+1))/2))$
AT:matrix([1,1],
[3,1],
[3,3],
[1,3]
);
print("AT=",AT,"S=",area(transpose(AT)))$
参考にさせていただいたプログラムですが 当方には解読不能でした。それなりにアレンジしたものを下記に示します。
(オリジナルを改変したもの)
/*area:面積計算関数*/
area(B):=block(x:col(B,1),y:col(B,2),n:length(B)-1,
return(sum(x[i]*y[i+1]-x[i+1]*y[i],i,1,n)/2)[1])$
AT:matrix([1,1],
[3,1],
[3,3],
[1,3],[1,1]
)$
print("AT=",AT,"S=",area(AT))$
本稿のプログラムは これをベースにしています。
本稿のやり方で処理します。
Mery で以下のスクリプトを ファイル名 "100-8.bat" で保存して実行してください。
/* regular square
@echo off & cls
path C:\maxima-5.47.0\bin;%path%
type %0 | maxima --very-quiet
pause
goto:eof
*/
Pt : matrix([1, 1], [3, 1], [3, 3], [1, 3], [1, 1])$
n : length(Pt)-1$
x : col(Pt, 1)$
y : col(Pt, 2)$
Ax : sum( (s[i] : x[i] * y[i+1] - x[i+1] * y[i]) / 2, i, 1, n)[1]$
gx : sum(s[i] * (x[i] + x[i+1]) / (6 * Ax), i, 1, n)[1]$
gy : sum(s[i] * (y[i] + y[i+1]) / (6 * Ax), i, 1, n)[1]$
Ix : sum(s[i] * ( (y[i] + y[i+1]) ^ 2 + y[i] ^ 2 + y[i+1] ^ 2) / 24, i, 1, n)[1]$
Iy : sum(s[i] * ( (x[i] + x[i+1]) ^ 2 + x[i] ^ 2 + x[i+1] ^ 2) / 24, i, 1, n)[1]$
Ixy : sum(s[i] * ( (x[i] + x[i+1]) * (y[i] + y[i+1]) + x[i] * y[i] + x[i+1] * y[i+1]) / 24,\
i, 1, n)[1]$
Ixc : Ix - Ax * gy ^ 2$
Iyc : Iy - Ax * gx ^ 2$
Ixyc : Ixy - Ax * gx * gy$
alpha : if Ixyc = 0 then 0 else atan2(2 * Ixyc, Iyc - Ixc) / 2$
print("area")$
print(" Ax =", float(ev(Ax)))$
print("elastic center")$
print(" gx =", float(ev(gx)))$ print(" gy =", float(ev(gy)))$
print("moment of inertia of area")$
print(" Ix =", float(ev(Ix)))$ print(" Iy =", float(ev(Iy)))$
print("centrifugal moment of area")$
print(" Ixy =", float(ev(Ixy)))$
print("moment of inertia of area by elastic center")$
print(" Ixc =", float(ev(Ixc)))$ print(" Iyc =", float(ev(Iyc)))$
print("centrifugal moment of area by elastic center")$
print(" Ixyc =", float(ev(Ixyc)))$
print("principal axis of area")$
print(" alpha =", float(ev(alpha)))$
quit()$
計算結果は
となります。
あらかじめ関数を作っておけば ruby なら
pt = [[1, 1], [3, 1], [3, 3], [1, 3], [1, 1]]
puts "area"
puts " Ax = %s" % a = pt.ax
puts "elastic center"
puts " gx = %s" % gx = pt.gx; puts " gy = %s" % gy = pt.gy
ix = pt.ix + a * gy ** 2
iy = pt.iy + a * gx ** 2
ixy = pt.jxy + a * gx * gy
puts "moment of inertia of area"
puts " Ix = %s" % ix ; puts " Iy = %s" % iy
puts "centrifugal moment of area"
puts " Ixy = %s" % ixy
puts "moment of inertia of area by elastic center"
puts " Ixc = %s" % pt.ix; puts " Iyc = %s" % pt.iy
puts "centrifugal moment of area by elastic center"
puts " Ixyc = %s" % pt.jxy
puts "principal axis of area"
puts " alpha = %s" % pt.alpha
__END__
計算結果は(外部変形用ライブラリ jw.rb / 未公開版 を利用)
となります。
数式処理ソフト maxima で座標を変数とすれば 多角形の断面性能の計算式を求めることがわかりました。
次回は ruby/座標法/断面性能 を取り上げます。木材や鋼材の呼称を座標扱いしてみます。