有限元动力学程序Python示例

有限元方法（Finite Element Method, FEM）是一种计算机模拟技术，广泛应用于工程和数学领域，用于求解偏微分方程和积分方程。在动力学问题中，有限元方法可以用于分析结构在时间变化的荷载作用下的响应，如振动和热传递等问题。

下面是一个简单的Python示例，使用有限元方法解决一个简单的一维弹性梁的动力学问题。这个示例使用了`numpy`进行矩阵运算，`matplotlib`进行结果可视化。

import numpy as np
import matplotlib.pyplot as plt

# 参数定义
E = 210e9  # 杨氏模量 (Pa)
rho = 7800  # 密度 (kg/m^3)
A = 0.01  # 截面积 (m^2)
I = 1e-5  # 惯性矩 (m^4)
L = 2  # 梁的长度 (m)
n_elements = 10  # 单元数量

# 创建系统矩阵
K = np.zeros((2*n_elements, 2*n_elements))  # 刚度矩阵
M = np.zeros((2*n_elements, 2*n_elements))  # 质量矩阵
F = np.zeros((2*n_elements, 1))  # 载荷向量

# 计算每个单元的局部刚度矩阵和质量矩阵
for i in range(n_elements):
    x1 = i * L / n_elements
    x2 = (i + 1) * L / n_elements
    k_local = (E * A / L) * np.array([[1, -1], [-1, 1]])
    m_local = (rho * A * L / 6) * np.array([[2, 1], [1, 2]])
    K[2*i:2*i+2, 2*i:2*i+2] += k_local
    M[2*i:2*i+2, 2*i:2*i+2] += m_local

# 应用载荷和边界条件
F[0] = 1000  # 假设在梁的左端施加一个向下的力
M[0, 0] = M[1, 1] = 0  # 左端固定，旋转自由度为0

# 解决动态方程
omega, V = np.linalg.eig(np.linalg.inv(M) @ K)  # 求解特征值和特征向量
modes = V[:, :2]  # 取前两阶模态

# 绘制模态形状
plt.plot(np.linspace(0, L, n_elements), modes[0, :n_elements], label='Mode 1')
plt.plot(np.linspace(0, L, n_elements), modes[1, :n_elements], label='Mode 2')
plt.xlabel('Position along the beam (m)')
plt.ylabel('Displacement (m)')
plt.title('Modal Shapes')
plt.legend()
plt.show()

这个示例仅为了演示有限元动力学程序的一个非常简单的应用，实际应用中的问题会更加复杂，需要考虑更多的因素，如非线性材料行为、复杂的边界条件、时间依赖的载荷等。