閉じた節点情報から断面二次モーメントを出力するプログラムを考えます。
これをもとに断面一次モーメント、断面積、断面係数を算出することも容易です。
定義式
断面二次モーメントは、
$$ I_z = \int_A y^2 dA \tag{1} $$
で、与えられます。断面内のあらゆる点について $y^2$を計算し、それの総和、ということですね。
コードを書く
今回引数は頂点情報の配列で渡すことを前提とします。例えば下記。
多角形なので、N番目の頂点と、N+1の頂点とそれらの垂線で構成される台形のメッシュで分割することができます。
上の図の場合は、青エリアは正数、赤エリアでは負数で返してやればOK。
台形の場合の式(1)は愚直に解けて、斜辺の一次関数を$y = ax + b$とすれば、
I_z = \int_{x_i}^{x_{i+1}} \int_{y_i}^{y_{i+1}} y^2 dy dx
= \frac{1}{12a} \left\{ ( ax_{i+1} + b )^4 - ( ax_i + b)^4 \right\} \tag{2a}
ただし、$a=0$ （つまり長方形）のときは、
I_z = \frac{1}{3} b^3 \left(x_{i+1} - x_i \right) \tag{2b}
となります。
コード
(2a)と(2b)をそのまま式に起こしてあげれば良い。引数のnodeは、
node = [[x1, y1], [x2, y2], ...]
の形式を期待。でも、あまりうまくいっていないのでnumpy的な方法にした方がいいと思います。
def moment_of_inertia_of_area(nodes):
nodes = move_x_to_positive(nodes) # 後述
func = reformat_nodes_to_func(nodes) # 後述
ans = 0.0
for f in func:
x1, x2 = f[0], f[1]
a, b = f[2][0], f[2][1]
if a == 0:
inte = (b ** 3) * (x2 - x1) / 3
else:
inte = ((a * x2 + b) ** 4 - (a * x1 + b) ** 4) / (12 * a)
ans += inte
return ans
途中で出てくるmove_x_to_positive(node)は、xを0以上に移動するための関数。
reformat_nodes_to_funcは、[x1, x2, a, b]の配列にして返す関数。この辺は引数の取り方にもよるので、適宜調整必要かと思います。
def move_x_to_positive(nodes):
min_x = min(list(map(lambda node: node[0], nodes)))
nodes = list(map(lambda node: [node[0] - min_x, node[1]], nodes))
return nodes
def reformat_nodes_to_func(nodes):
func_list = [
[nodes[i][0], nodes[i + 1][0],
calc_slope_and_intercept(
nodes[i][0], nodes[i][1], nodes[i + 1][0], nodes[i + 1][1])]
for (i, x) in enumerate(nodes[:-1])
]
return func_list
def calc_slope_and_intercept(xi, yi, xj, yj):
try:
slope = (yi - yj) / (xi - xj)
y_intercept = yi - slope * xi
except ZeroDivisionError:
slope, y_intercept = 0, 0
return slope, y_intercept
この辺汚いですな。2年くらい前のコードなので....
まあ、動けばOK！！！！
以上です。ライブラリなしで計算できるのはメリット。