用Python手搓一个有限元分析器：从5节点三角形单元到云图可视化（附完整代码）
用Python手搓一个有限元分析器：从5节点三角形单元到云图可视化（附完整代码）
有限元分析（FEA）作为工程仿真领域的核心技术，广泛应用于结构力学、热传导、流体动力学等场景。对于工程师和科研人员而言，理解有限元背后的数学原理固然重要，但亲手实现一个迷你有限元程序更能加深对算法细节的掌握。本文将带你用Python从零构建一个支持5节点三角形单元的平面应力分析器，并实现位移/应力云图的可视化。
1. 环境准备与基础理论
1.1 工具链选择
我们选择Python生态中的两个核心库：
- NumPy：处理矩阵运算和线性代数计算
- Matplotlib：实现计算结果的可视化
安装依赖只需一行命令：
pip install numpy matplotlib1.2 三角形单元理论基础
三节点三角形单元（常应变单元）是最基础的有限元单元类型，其刚度矩阵推导基于以下关键公式：
应变-位移关系： $$ \mathbf{\varepsilon} = \mathbf{Bu} $$
应力-应变关系（平面应力）： $$ \mathbf{\sigma} = \mathbf{D\varepsilon} $$
单元刚度矩阵： $$ \mathbf{K^e} = tA\mathbf{B^TDB} $$
其中：
- $t$为厚度
- $A$为单元面积
- $\mathbf{B}$为应变-位移矩阵
- $\mathbf{D}$为弹性矩阵
2. 代码实现步骤
2.1 网格定义与预处理
首先定义节点坐标和单元连接关系：
import numpy as np # 材料参数 E = 210e9 # 弹性模量(Pa) nu = 0.2 # 泊松比 t = 0.01 # 厚度(m) # 5节点示例网格 nodes = np.array([ [0, 1], # 节点0 [1, 1], # 节点1 [0, 0], # 节点2 [1, 0], # 节点3 [2, 0] # 节点4 ]) # 3个三角形单元 elements = np.array([ [0, 2, 3], # 单元0 [1, 3, 4], # 单元1 [3, 1, 0] # 单元2 ])2.2 单元刚度矩阵计算
实现三角形单元刚度矩阵计算函数：
def compute_element_stiffness(E, nu, t, node_coords): """计算三节点三角形单元刚度矩阵""" # 计算单元面积 mat = np.array([ [1, node_coords[0,0], node_coords[0,1]], [1, node_coords[1,0], node_coords[1,1]], [1, node_coords[2,0], node_coords[2,1]] ]) A = 0.5 * abs(np.linalg.det(mat)) # 弹性矩阵(平面应力) D = (E/(1-nu**2)) * np.array([ [1, nu, 0], [nu, 1, 0], [0, 0, (1-nu)/2] ]) # 应变-位移矩阵B x, y = node_coords[:,0], node_coords[:,1] B = np.zeros((3,6)) for i in range(3): j, m = (i+1)%3, (i+2)%3 B[0, 2*i] = y[j]-y[m] B[1, 2*i+1] = x[m]-x[j] B[2, 2*i] = x[m]-x[j] B[2, 2*i+1] = y[j]-y[m] B = B / (2*A) return t * A * (B.T @ D @ B)2.3 全局刚度矩阵组装
将单元刚度矩阵组装到全局矩阵：
n_nodes = len(nodes) K_global = np.zeros((2*n_nodes, 2*n_nodes)) for elem in elements: K_elem = compute_element_stiffness(E, nu, t, nodes[elem]) # 将单元刚度矩阵添加到全局矩阵 for i in range(3): for j in range(3): K_global[2*elem[i]:2*elem[i]+2, 2*elem[j]:2*elem[j]+2] += K_elem[2*i:2*i+2, 2*j:2*j+2]3. 边界条件与求解
3.1 施加边界条件
设置固定边界和载荷：
# 固定节点3,4,5的x,y位移 (对应索引2,3,4) fixed_dofs = [4,5,6,7,8,9] # 节点0施加x方向载荷 force = np.zeros(2*n_nodes) force[0] = 1e6 # 1MN的x方向力 # 处理边界条件 free_dofs = [i for i in range(2*n_nodes) if i not in fixed_dofs] K_free = K_global[np.ix_(free_dofs, free_dofs)] force_free = force[free_dofs]3.2 求解位移
求解线性方程组得到自由度的位移：
displacements = np.zeros(2*n_nodes) displacements[free_dofs] = np.linalg.solve(K_free, force_free) print("节点位移：") for i in range(n_nodes): print(f"节点{i}: ux={displacements[2*i]:.2e}, uy={displacements[2*i+1]:.2e}")4. 结果可视化
4.1 位移云图绘制
使用Matplotlib绘制位移场：
import matplotlib.pyplot as plt import matplotlib.tri as tri def plot_displacement(nodes, elements, displacements, component='x'): """绘制位移云图""" disp = displacements[::2] if component == 'x' else displacements[1::2] fig, ax = plt.subplots(figsize=(8,6)) triang = tri.Triangulation(nodes[:,0], nodes[:,1], elements) # 绘制填充云图 tcf = ax.tricontourf(triang, disp, levels=20, cmap='jet') fig.colorbar(tcf, label=f'Displacement ({component})') # 叠加网格线 ax.triplot(triang, 'k-', lw=0.5, alpha=0.5) ax.set_title(f'{component.upper()}方向位移云图') ax.set_aspect('equal') plt.show() plot_displacement(nodes, elements, displacements, 'x') plot_displacement(nodes, elements, displacements, 'y')4.2 应力计算与可视化
计算单元应力并绘制云图：
def compute_stress(E, nu, node_coords, displacements): """计算单元应力""" # 与刚度矩阵相同的B矩阵计算 x, y = node_coords[:,0], node_coords[:,1] A = 0.5 * abs(x[0]*(y[1]-y[2]) + x[1]*(y[2]-y[0]) + x[2]*(y[0]-y[1])) B = np.zeros((3,6)) for i in range(3): j, m = (i+1)%3, (i+2)%3 B[0, 2*i] = y[j]-y[m] B[1, 2*i+1] = x[m]-x[j] B[2, 2*i] = x[m]-x[j] B[2, 2*i+1] = y[j]-y[m] B = B / (2*A) D = (E/(1-nu**2)) * np.array([ [1, nu, 0], [nu, 1, 0], [0, 0, (1-nu)/2] ]) return D @ B @ displacements # 计算所有单元应力 stresses = [] for elem in elements: elem_disp = np.concatenate([displacements[2*elem[i]:2*elem[i]+2] for i in range(3)]) stresses.append(compute_stress(E, nu, nodes[elem], elem_disp)) stresses = np.array(stresses) # 绘制应力云图 def plot_stress(nodes, elements, stresses, stress_type='x'): """绘制应力云图""" if stress_type == 'x': data = stresses[:,0] # σxx elif stress_type == 'y': data = stresses[:,1] # σyy else: data = stresses[:,2] # τxy fig, ax = plt.subplots(figsize=(8,6)) triang = tri.Triangulation(nodes[:,0], nodes[:,1], elements) # 三角形常应力绘图 tpc = ax.tripcolor(triang, facecolors=data, edgecolors='k', cmap='jet') fig.colorbar(tpc, label=f'Stress ({stress_type})') ax.set_title(f'应力云图 (σ{stress_type})') ax.set_aspect('equal') plt.show() plot_stress(nodes, elements, stresses, 'x') plot_stress(nodes, elements, stresses, 'y') plot_stress(nodes, elements, stresses, 'xy')5. 扩展与优化建议
5.1 性能优化技巧
对于大规模问题，稀疏矩阵存储可以显著减少内存占用：
from scipy.sparse import lil_matrix # 使用稀疏矩阵存储 K_global = lil_matrix((2*n_nodes, 2*n_nodes)) for elem in elements: K_elem = compute_element_stiffness(E, nu, t, nodes[elem]) for i in range(3): for j in range(3): K_global[2*elem[i]:2*elem[i]+2, 2*elem[j]:2*elem[j]+2] += K_elem[2*i:2*i+2, 2*j:2*j+2]5.2 高阶单元实现
将三节点单元扩展为六节点二次单元，只需修改B矩阵的计算方式：
def compute_quadratic_B_matrix(node_coords): """六节点三角形单元的B矩阵计算""" # 包含中点节点的形函数导数计算 # 具体实现取决于采用的形函数形式 pass5.3 结果验证方法
通过与解析解或商业软件对比验证结果可靠性：
# 示例：与理论解对比 theory_solution = {...} # 填入理论解 error = np.linalg.norm(displacements - theory_solution) print(f"相对误差: {error:.2e}")在实际项目中，建议先用简单模型验证代码正确性，再逐步扩展到复杂问题。对于非线性或动力学问题，可以考虑在现有框架上增加时间积分和非线性求解器。