【python】【材料力学】计算任意界面的质心坐标及惯性矩

背景

etabs自定义截面

公式

This code uses the following formulas:

Area of the polygon:

$$A=\frac{1}{2}\sum _{i=1}^n(x_i{y}_{i+1}-x_{i+1}{y}_i)$$

Centroid coordinates:

$$C_x=\frac{1}{6A}\sum _{i=1}^n(x_i+x_{i+1})(x_i{y}_{i+1}-x_{i+1}{y}_i)$$

$$C_y=\frac{1}{6A}\sum _{i=1}^n(y_i+y_{i+1})(x_i{y}_{i+1}-x_{i+1}{y}_i)$$

Moment of inertia:

$$Ix=\frac{1}{12}\sum {i=0}^{n-1}(y{i}^2+yiyi+1+y{i+1}^2)(xiyi+1-xi+1yi)$$
$$Iy=\frac{1}{12}\sum {i=0}^{n-1}(x{i}^2+xixi+1+x{i+1}^2)(xiyi+1-xi+1yi)$$
$$Ixy=\frac{1}{24}\sum {i=0}^{n-1}(2xiyi+xi+1yi+2xi+1yi+1+xiyi+1)(xiyi+1-xi+1yi)$$
其中，$n$ 是多边形的顶点数，$(xi,yi)$ 表示多边形中第 $i$ 个顶点的坐标，$(xi+1,yi+1)$ 表示多边形中第 $i+1$ 个顶点的坐标，$Ix$ 代表多边形相对于 $x$ 轴的惯性矩，$Iy$ 代表多边形相对于 $y$ 轴的惯性矩，$I_{xy}$ 代表多边形相对于 $x$ 和 $y$ 轴的二次惯性矩。

代码

import matplotlib.pyplot as plt

def polygon_area(vertices):
    n = len(vertices)
    area = 0
    for i in range(n):
        area += (vertices[i][0] * vertices[(i + 1) % n][1] - vertices[(i + 1) % n][0] * vertices[i][1])
    return 0.5 * abs(area)

def centroid(vertices):
    n = len(vertices)
    cx, cy = 0, 0
    area = polygon_area(vertices)
    for i in range(n):
        factor = (vertices[i][0] * vertices[(i + 1) % n][1] - vertices[(i + 1) % n][0] * vertices[i][1])
        cx += (vertices[i][0] + vertices[(i + 1) % n][0]) * factor
        cy += (vertices[i][1] + vertices[(i + 1) % n][1]) * factor
    return (1 / (6 * area)) * cx, (1 / (6 * area)) * cy

def moment_of_inertia(vertices):
    n = len(vertices)
    Ix, Iy, Ixy = 0, 0, 0
    area = polygon_area(vertices)
    for i in range(n):
        factor = (vertices[i][0] * vertices[(i + 1) % n][1] - vertices[(i + 1) % n][0] * vertices[i][1])
        Ix += (vertices[i][1]**2 + vertices[i][1]*vertices[(i + 1) % n][1] + vertices[(i + 1) % n][1]**2) * factor
        Iy += (vertices[i][0]**2 + vertices[i][0]*vertices[(i + 1) % n][0] + vertices[(i + 1) % n][0]**2) * factor
        Ixy += (vertices[i][0]*vertices[(i + 1) % n][1] + 2*vertices[i][0]*vertices[i][1] + 2*vertices[(i + 1) % n][0]*vertices[(i + 1) % n][1] + vertices[(i + 1) % n][0]*vertices[i][1]) * factor
    return (1 / 12) * abs(Ix), (1 / 12) * abs(Iy), (1 / 24) * abs(Ixy)

def plot_polygon(vertices, centroid):
    polygon = plt.Polygon(vertices, fill=None, edgecolor='r', linewidth=2)
    plt.gca().add_patch(polygon)
    plt.scatter(*centroid, color='blue', marker='o', s=100, label='Centroid')
    plt.axis('equal')
    plt.legend()
    plt.show()

# 可增加，从cad中读取
vertices = [(0, 0), (30, 0), (30, 30), (0, 30), (-20, 20)]

area = polygon_area(vertices)
centroid_coordinates = centroid(vertices)
Ix, Iy, Ixy = moment_of_inertia(vertices)

print("Area:", area)
print("Centroid:", centroid_coordinates)
print("Moment of Inertia: Ix =", Ix, "Iy =", Iy, "Ixy =", Ixy)

plot_polygon(vertices, centroid_coordinates)

总结

常用的markdown编辑器可以渲染latex格式代码
常用公式可以编程