Python 有限元从零开始

简介

学习有限元的最佳路径是自己动手实现一遍。本教程不依赖任何现成 FEM 框架，从三角单元出发，用纯 NumPy/SciPy 构建 2D 线弹性有限元求解器，覆盖从网格生成到后处理的完整流程。适合想深入理解 FEM 底层原理的工程师与学生。

学习路径

章节 1：三角单元与形状函数（2–3 天）

目标：理解 FEM 的核心——从弱形式到单元矩阵

实现：

import numpy as np

# 三角单元形状函数 N1, N2, N3
# 面积坐标 (L1, L2, L3)
def shape_functions(xi, eta):
    return [1 - xi - eta, xi, eta]

# 形函数导数（Jacob 变换）
def dN_dXY(xi, eta, coords):
    # coords: 3x2 节点坐标
    J = np.array([[coords[1,0]-coords[0,0], coords[2,0]-coords[0,0]],
                  [coords[1,1]-coords[0,1], coords[2,1]-coords[0,1]]])
    detJ = np.linalg.det(J)
    # ... 组装 B 矩阵

章节 2：刚度矩阵组装（2–3 天）

目标：掌握全局刚度矩阵的组装逻辑

关键：

• 单元刚度矩阵 K_e = ∫ B^T D B dΩ
• 高斯积分（2D 3点/7点积分）
• 全局矩阵组装（稀疏矩阵 COO 格式）

章节 3：边界条件与求解（1–2 天）

目标：掌握 Dirichlet/Neumann 边界条件的处理

关键：

• 硬边界（直接置零）/ 惩罚法
• SciPy 稀疏线性求解器
• 结果回代与应力恢复

章节 4：网格生成与后处理（1–2 天）

目标：自己写 Delaunay 三角剖分 + 结果可视化

工具：

• scipy.spatial.Delaunay 做三角剖分
• matplotlib 等值线 / 伪彩色图

适合人群

• 机械/土木/航空航天专业学生，理解 FEM 公式推导但想动手验证
• 有编程基础，想从第一性原理理解 FEM 而非仅会操作软件
• 为学习 FEniCS / Abaqus 打基础（先自己实现再用框架会理解更深）