4.3 有限元法

4.3 有限元法

差分方法把导数换成差商，前提取决于规整网格；真实世界的机翼、发动机缸体、生物器官从不规整。有限元换一条路线：不做差商近似，而是承认解算不出来，转而在一个由分片多项式拼成的有限维空间里，寻找弱形式（2.4 节的遗产）意义下的最佳逼近。几何自由度由此彻底打开——这是它统治结构工程半个多世纪的理由。本节把一维泊松问题从头手算一遍，再看程序实现与收敛验证。

从弱形式到分片拼图

以 −u'' = f（零边值）为例。两边乘检验函数 v（边界取零）并分部积分：

∫ u'·v' dx = ∫ f·v dx。

这就是弱形式：一阶导数替代二阶导数，对 u 的光滑要求降了一级。伽辽金的思路直截了当：把 u 与 v 都限制在同一个有限维空间 V_h（下标 h 记单元尺寸），u = Σ U_j·φ_j，代入弱形式，得到线性方程组

K·U = F，其中 K_ij = ∫ φ_i'·φ_j' dx，F_i = ∫ f·φ_i dx。

K 叫刚度矩阵（结构工程的传统叫法），F 叫载荷向量。理论保证（迦辽金正交性）：这个 U 是真解在 V_h 中的能量范数最佳逼近——有限元的每一次离散都是最优投影，这是它数学上的贵族血统。

基函数的选择是艺术所在。帽形函数是最经典的一维选择：每个节点 i 上放一个"帐篷"，在本单元内线性上升、相邻单元内线性下降、其余处为零、节点处取一。任意连续分片线性函数都是它们的线性组合——解被想象成由帐篷搭出的折线，节点值 U_j 就是帐篷的高度。

图解：帽形基函数与单元剖分

帽形基函数与刚度矩阵组装示意

手算一遍：两个单元上的泊松问题

取区间 [0,1]、f = 1、零边值，剖成两个单元、三个节点（x = 0、0.5、1）。单元长度 h = 0.5。帽形基只有中间节点的 φ₂ 是"完整帐篷"，端点被边界条件清零。逐项积分（φ' 在每个单元内是常数 ±1/h）：

K = (1/h)·[2, −1; −1, 2]（缩到内部自由度后是 1x1：K = 2/h），F = f·h = 0.5。

解得 U = f·h²/2·2/2…——直接算：内部方程 (2/h)·U₂ = f·h，U₂ = f·h²/2 = 0.125。真解 u = x(1−x)/2 在 x = 0.5 处恰是 0.125——两个单元就精确命中节点值。这不是巧合：一维线性问题里，分片线性空间对二次真解的节点插值是精确的（弹性力学里的"超收敛"现象）。工程现场用几十万个单元逼近的真实结构，遵循同一个组装逻辑：单元矩阵是局部小块，共享节点处叠加——网格生成的复杂性被消化在这一步之前，组装本身是纯机械。

数值实验：程序组装与收敛验证

把两个单元的幸运推广成系统的收敛表：N 个单元线性有限元解泊松问题，与真解 u = x(1−x)/2 对表，预期误差 O(h²)（最大范数下线性元对二阶真解约为 h²/8 量级）。

import numpy as np

def fem_poisson(n):
    h = 1.0/(n+1)
    K = (2*np.eye(n) - np.eye(n, k=1) - np.eye(n, k=-1)) / h
    F = np.full(n, h)
    U = np.linalg.solve(K, F)
    x = np.linspace(0, 1, n+2)[1:-1]
    return x, U, h

for n in [7, 15, 31, 63]:
    x, U, h = fem_poisson(n)
    exact = x*(1-x)/2
    err = np.abs(U - exact).max()
    print(f"n={n:3d}  h={h:.4f}  最大误差 {err:.3e}")
​

输出里误差随 n 加倍下降约四倍——二阶收敛兑现。对照 4.2 的差分表会发现两者数字惊人接近：一维均匀网格上，线性有限元的刚度方程与三点中心差分代数上完全等价（线性元在均匀网格上的离散方程就是二阶中心差分）。区别在身份不在数值：有限元的最优性、边界条件的自然处理（诺伊曼条件自动蕴含在弱形式里，无需特殊格式）、以及向任意单元形状的推广，都是差分路线不具备的。

有限元的看家本领与代价

三条本领撑起它的行业地位。几何任意性：三角形、四面体可以贴合任何边界，曲边界用等参元映射；机翼网格生成虽贵，但可行。边界条件的自然性：狄利克雷条件强加于节点，诺伊曼与罗宾条件自动进入载荷向量——热流、面力这类"力的边界"在变分框架里是原生公民。误差可控与自适应：事后误差估计子指出误差集中区域，自动加密局部网格，把算力花在刀刃上——现代结构仿真的自适应循环（解、估、加、再解）是第七章工作流程的核心环节。代价也明确：组装与求解大型稀疏系统的成本、网格生成的人工负担、以及单元形状过畸（质量差的三角形）对精度的伤害。这些工程细节在第七章软件生态一节再收拢。

💡 关键直觉：有限元不是"把差分换成三角形"，而是换了一种提问方式——不问"导数近似得准不准"，而问"在允许的折线函数里，谁的加权平均残差最小"。提问方式一换，几何、边界、自适应三扇门同时打开。

从一维到二维：单元、质量矩阵与约束

三维世界的三件常备兵器值得预先打点。单元形状：二维三角形与四边形、三维四面体与六面体，基函数从帽形长成"金字塔"或"小山"；曲边界用等参映射处理——单元顶点贴真实边界，映射把直边单元弯成曲边，几何误差从 O(h) 压到 O(h²) 以上。单元质量有一个朴素的判据：内角不要过尖、长短比不要过大，畸形单元的误差常数会恶化到吃掉加密收益。

质量矩阵：瞬态问题（u_t 项进入弱形式）会多出一个 M_ij = ∫φ_i·φ_j dx，叫质量矩阵。它有两种取法：一致质量矩阵（如实积分，非对角元非零）与集中质量矩阵（把每行对角化，物理直觉是"把单元质量分摊到节点"）。显式步进时集中质量免去解方程，是显式动力学（碰撞、冲击）的标配；一致质量对波的传播频散更小。两者的取舍是瞬态有限元的第一道工艺选择。

约束处理：狄利克雷边界在程序里有三种写法——直接划行置单位（改矩阵结构，简单粗暴）、罚函数（把约束乘大系数加进刚度矩阵，不改结构但引入病态风险）、拉格朗日乘子（严格但多出未知量）。商业软件默认多为罚函数或划行法，理解三种写法对读懂求解日志至关重要——约束没施加上的症状（刚体漂移、解整体平移）在结构计算里是常见报障原因。

本节要点

• 弱形式路线：分部积分降光滑要求，伽辽金投影给出能量范数最优逼近；
• 帽形基与组装：分片线性搭出解空间，单元矩阵局部组装，均匀一维网格退化为三点差分；
• 手算验证：两单元泊松问题节点值精确命中，超收敛现象可亲手复算；
• 看家本领：任意几何、自然边界条件、误差估计驱动的自适应加密；
• 与差分的关系：数值相同、身份不同——最优性与几何普适性是真正的分水岭。

有限元的弱形式对平衡问题如鱼得水，但流动与激波问题最在乎另一件事——守恒。下一节的有限体积法把守恒律变成离散格式的原生属性，让每一步迭代都分毫不差地守恒。