用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 (i1)%3, (i2)%3 B[0, 2*i] y[j]-y[m] B[1, 2*i1] x[m]-x[j] B[2, 2*i] x[m]-x[j] B[2, 2*i1] 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*i2, 2*j:2*j2]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*i1]:.2e})4. 结果可视化4.1 位移云图绘制使用Matplotlib绘制位移场import matplotlib.pyplot as plt import matplotlib.tri as tri def plot_displacement(nodes, elements, displacements, componentx): 绘制位移云图 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, levels20, cmapjet) fig.colorbar(tcf, labelfDisplacement ({component})) # 叠加网格线 ax.triplot(triang, k-, lw0.5, alpha0.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 (i1)%3, (i2)%3 B[0, 2*i] y[j]-y[m] B[1, 2*i1] x[m]-x[j] B[2, 2*i] x[m]-x[j] B[2, 2*i1] 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_typex): 绘制应力云图 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, facecolorsdata, edgecolorsk, cmapjet) fig.colorbar(tpc, labelfStress ({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*i2, 2*j:2*j2]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})在实际项目中建议先用简单模型验证代码正确性再逐步扩展到复杂问题。对于非线性或动力学问题可以考虑在现有框架上增加时间积分和非线性求解器。