有限元方法中全局刚度矩阵和负载向量的python代码

在有限元方法中，全局刚度矩阵（K矩阵）和负载向量（F向量）是求解结构力学问题的核心部分，通常用于线性方程组的建立。以下是使用Python（numpy库）创建简单二维梁元素的全局刚度矩阵和负载向量的一个基本示例：

import numpy as np

# 梁的节点数量和元素数量
n_nodes = 4  # 假设有四个节点
n_elements = 2  # 假设有两个元素

# 一维元素的刚度系数（这里假设都是单位长度，需要根据实际材料和边界条件调整）
element_stiffness = 1.0

# 全局节点编号和位置（这里仅作示例，实际应用中会根据几何布局计算）
nodes_positions = [(0, 0), (1, 0), (1, 1), (0, 1)]
node_numbers = [i for i in range(n_nodes)]

# 构建局部和全局索引
local_element_indices = [(0, 1), (1, 2)]  # 每个元素的两端节点
global_element_indices = []
for element_local in local_element_indices:
    global_element_indices.extend([node_numbers[i] for i in element_local])

# 索引转换函数
def get_global_index(local_index):
    return node_numbers.index(local_index[0]) * n_nodes + node_numbers.index(local_index[1])

# 创建局部刚度矩阵
local_dofs = [0, 1]
local_stiffness_matrix = np.zeros((len(local_dofs), len(local_dofs)))
for e in range(n_elements):
    for i in range(len(local_dofs)):
        for j in range(len(local_dofs)):
            local_stiffness_matrix[i, j] += element_stiffness * (
                1 if local_element_indices[e][i] < local_element_indices[e][j] else -1
            )

# 转换为全局刚度矩阵
K = np.zeros((n_nodes * n_dofs, n_nodes * n_dofs))
for global_index, local_index in enumerate(global_element_indices):
    K[get_global_index(local_dofs[0]), get_global_index(local_dofs)] += local_stiffness_matrix[
        :, :
    ]
    K[get_global_index(local_dofs[1]), get_global_index(local_dofs)] -= local_stiffness_matrix[
        :, :
    ]

# 创建负载向量，例如均匀分布的集中力
load_vector = np.zeros(n_nodes * n_dofs)
central_load = 10.0
load_position = node_numbers.index(0) * n_nodes + node_numbers.index(1)  # 第二个节点的第二个DOF
load_vector[load_position] = central_load

# 结果展示
print("全局刚度矩阵(K):")
print(K)
print("\n负载向量(F):")
print(load_vector)