有限元程式碼

redufa發表於2024-10-13

import numpy as np

class FiniteElementAnalysis:
    def __init__(self, nodes, elements, material_properties):
        self.nodes = nodes
        self.elements = elements
        self.material_properties = material_properties
        self.global_stiffness_matrix = None

    def shape_functions(self, xi, eta, zeta):
        # 二次等參單元的形函式
        N = [
            (1 - xi) * (1 - eta) * (1 - zeta) / 8,
            (1 + xi) * (1 - eta) * (1 - zeta) / 8,
            (1 + xi) * (1 + eta) * (1 - zeta) / 8,
            (1 - xi) * (1 + eta) * (1 - zeta) / 8,
            (1 - xi) * (1 - eta) * (1 + zeta) / 8,
            (1 + xi) * (1 - eta) * (1 + zeta) / 8,
            (1 + xi) * (1 + eta) * (1 + zeta) / 8,
            (1 - xi) * (1 + eta) * (1 + zeta) / 8
        ]
        return N

    def compute_element_stiffness_matrix(self, element):
        # 計算單元剛度矩陣
        node_ids = element['nodes']
        E = self.material_properties['E']
        nu = self.material_properties['nu']
        t = self.material_properties['t']
        Ke = np.zeros((24, 24))

        # 高斯積分點
        gauss_points = [(0.57735, 0.57735, 0.57735), (-0.57735, 0.57735, 0.57735), ...]  # 更多積分點

        for xi, eta, zeta in gauss_points:
            N = self.shape_functions(xi, eta, zeta)
            # 計算雅可比矩陣和其逆
            J = np.array([[...], [...], [...]])  # 雅可比矩陣
            det_J = np.linalg.det(J)
            inv_J = np.linalg.inv(J)

            # 計算B矩陣
            B = np.array([[...], [...], [...]])  # B矩陣

            Ke += B.T @ np.array([[E / (1 - nu**2) * np.array([[1, nu, 0], [nu, 1, 0], [0, 0, (1 - nu) / 2]])]]) @ B * det_J * t

        return Ke

    def assemble_global_stiffness_matrix(self):
        num_nodes = len(self.nodes)
        self.global_stiffness_matrix = np.zeros((3 * num_nodes, 3 * num_nodes))

        for element in self.elements:
            Ke = self.compute_element_stiffness_matrix(element)
            for i, node_i in enumerate(element['nodes']):
                for j, node_j in enumerate(element['nodes']):
                    self.global_stiffness_matrix[3*node_i:3*node_i+3, 3*node_j:3*node_j+3] += Ke[3*i:3*i+3, 3*j:3*j+3]

    def apply_boundary_conditions(self, bc):
        for node_id, dof in bc.items():
            for d in dof:
                self.global_stiffness_matrix[3*node_id+d, :] = 0
                self.global_stiffness_matrix[:, 3*node_id+d] = 0
                self.global_stiffness_matrix[3*node_id+d, 3*node_id+d] = 1

    def solve(self, forces):
        displacements = np.linalg.solve(self.global_stiffness_matrix, forces)
        return displacements

# 示例資料
nodes = np.array([[0, 0, 0], [1, 0, 0], [1, 1, 0], [0, 1, 0], [0, 0, 1], [1, 0, 1], [1, 1, 1], [0, 1, 1]])
elements = [{'nodes': [0, 1, 2, 3, 4, 5, 6, 7]}]
material_properties = {'E': 210e9, 'nu': 0.3, 't': 1}

# 建立有限元分析物件
fea = FiniteElementAnalysis(nodes, elements, material_properties)
fea.assemble_global_stiffness_matrix()

# 應用邊界條件
bc = {0: [0, 1, 2], 1: [0, 1, 2]}  # 例如，節點0和1的x, y, z方向固定
fea.apply_boundary_conditions(bc)

# 施加外力
forces = np.zeros(24)
forces[3*2+2] = -1000  # 在節點2的z方向施加1000N的力

# 求解
displacements = fea.solve(forces)
print("Displacements:", displacements)

相關文章