断面性能の計算式(7)-正多角形-

座標法と数式処理ソフト 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章 断面の主軸

 


円に内接する正多角形 です。

Mery で以下のスクリプトを ファイル名 "100-7.bat" で保存して実行してください。


/* polygon
@echo off & cls
path C:\maxima-5.47.0\bin;%path%
  type %0 | maxima --very-quiet
  pause
  goto:eof
*/

n : 5$
a : 2 * %pi / n$
/*
Pt : matrix(
  [r * (1 + cos(a*0)), r * (1 - sin(a*0))],
  [r * (1 + cos(a*1)), r * (1 - sin(a*1))],
  [r * (1 + cos(a*2)), r * (1 - sin(a*2))],
  [r * (1 + cos(a*3)), r * (1 - sin(a*3))],
  [r * (1 + cos(a*4)), r * (1 - sin(a*4))],
  [r * (1 + cos(a*0)), r * (1 - sin(a*0))]
)$
*/
listp(b: map(
  lambda([i], r * [1 + cos(a*i), 1 - sin(a*i)]), 
  endcons(0, makelist(i, i, 0, n-1)) /* [0,1,2,3,4,0] */
))$
matrixp(Pt: apply('matrix, b))$

x : col(Pt, 1)$
y : col(Pt, 2)$

Ax : sum( (s[i] : x[i+1] * y[i] - x[i] * y[i+1]) / 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]$

print("area")$
print("  Ax =", factor(trigsimp(trigreduce(Ax))))$
print("  Ax =", float(ev(Ax)))$
print("elastic center")$
print("  gx =", factor(trigsimp(trigreduce(gx))))$ print("  gy =", factor(trigsimp(trigreduce(gy))))$
print("  gx =", float(ev(gx)))$ print("  gy =", float(ev(gy)))$
print("moment of inertia of area")$
print("  Ix =", factor(trigsimp(trigreduce(Ix))))$ print("  Iy =", factor(trigsimp(trigreduce(Iy))))$
print("  Ix =", float(ev(Ix)))$ print("  Iy =", float(ev(Iy)))$
quit()$

 


計算式は

となります。

 


座標を三角関数で与えたので複雑な計算式が返ってきました。当方の環境で n=23 は解けましたが n=24 から Ax 以外は発散しました。

 


プログラムは読みやすく書くことに気をつけています。自分でもわからなくなるからです。

 

 

次回は ホームページ「Maxima/面積計算 座標法 倍横距法」 を取り上げます。このブログを書くきっかけになったホームページです。