9.2 有限元方法入门

9.2 有限元方法入门

本节摘要：有限元不是黑箱，是第五章能量原理的离散化：单元刚度矩阵来自单元势能对节点位移求导，组装后解线性方程组就是找总势能最低的位移场。本节亲手装配杆与梁单元，用一致载荷向量与集中载荷向量的对账，讲清"误差从哪来、精度到哪去"。

从瑞利-里茨到单元装配

5.2 节的瑞利-里茨法有个先天短板：全局试探函数在复杂几何面前写不出来。有限元的对策朴素而致命——把结构切成小片，每片里用最简单的多项式当试探函数（杆单元线性、梁单元三次），片与片只在节点上握手。全局函数写不出，分段多项式永远写得出。总势能从"一个泛函"变成"所有单元势能之和"，对节点位移求导置零，得到 K 乘 u 等于 F 的线性方程组：K 是刚度矩阵，u 是节点位移，F 是节点力。第五章说"泛函极值缩成线性代数"，这里就是那间 algebra 车间的装配线。

单元刚度矩阵不需要死背：两节点杆单元 k 等于 EA 除以 L，梁单元的四乘四矩阵是三次埃尔米特形函数对势能求导的产物。真正要弄清的是载荷怎么进方程——同一个分布载荷，搬法不同，误差天差地别，这是本节两份对账记录的主角。

图：连续梁离散为单元与两种载荷搬运方式

对账一：受轴向均布载荷的杆

试件：长两米的杆，EA 等于二乘十的八次方牛，沿轴向均布剪来每米一万牛的载荷，左端固定。精确解是二次抛物线位移，自由端伸长 qL 平方除以二 EA，等于零点一毫米。用一个杆单元算三遍：精确载荷搬运（一致向量，按形函数把载荷分到两端各 qL/2）、粗糙搬运（全部堆到自由端）、以及解析真值。

# 杆单元：一致向量 vs 集中向量
EA, L, q = 2e8, 2.0, 1e4

k = EA / L                            # 单元刚度 N/m
u_exact = q*L**2/(2*EA)               # 解析自由端位移
u_consist = (q*L/2) / k               # 一致向量：两端各 qL/2
u_lumped  = (q*L)   / k               # 集中向量：总量堆到自由端

print(f"解析真值     = {u_exact*1e3:6.4f} mm")
print(f"一致载荷向量 = {u_consist*1e3:6.4f} mm（节点上精确）")
print(f"集中载荷向量 = {u_lumped*1e3:6.4f} mm（误差 {(u_lumped/u_exact-1)*100:.0f} %）")
​

解析真值     = 0.1000 mm
一致载荷向量 = 0.1000 mm（节点上精确）
集中载荷向量 = 0.2000 mm（误差 100 %）
​

判读：粗看像巧合，细看是定理——杆的精确位移场是二次的，而两节点杆单元内部位移是线性的，"内部"必然有误差；但节点这个位置上，线性插值恰好穿过真解，位移在节点上精确。一致载荷向量由变分原理（形函数加权）导出，把这份节点精确性完整送达；集中向量破坏了加权规则，误差翻倍。判读规则由此而来：分布载荷必须按一致向量搬运，自己手算节点载荷时别偷懒平分。

对账二：悬臂梁的两种搬运与一次精确命中

试件升级为梁：悬臂，长两米，EI 等于四点五乘十的六次方（5.2 节用过的同款数据），均布每米五千牛，外加一组对照——自由端集中十千牛。梁单元刚度矩阵四乘四，夹持端去掉两个自由度后只剩二乘二，用克莱姆法则手解。分布载荷两种搬法各算一遍，端部集中载荷也来一发——理论预言：端部载荷下梁单元位移场恰好是三次真解，一个单元就精确命中。

# 梁单元：夹持-自由，2x2 手解
import math

EI, L, q, P = 4.5e6, 2.0, 5e3, 1e4

# 自由度 (w_tip, th_tip) 的缩减刚度：EI/L3 乘 [[12, -6L], [-6L, 4L2]]
c  = EI/L**3
K11, K12, K22 = 12*c, -6*L*c, 4*L**2*c
det = K11*K22 - K12*K12

def solve(f1, f2):
    w = ( K22*f1 - K12*f2)/det
    th = (-K12*f1 + K11*f2)/det
    return w, th

w_c, _  = solve(q*L/2, -q*L**2/12)   # 一致向量（带节点力矩）
w_l, _  = solve(q*L/2, 0.0)          # 集中向量（丢力矩）
w_p, _  = solve(P, 0.0)              # 端部集中载荷

w_exact_udl = q*L**4/(8*EI)
w_exact_pt  = P*L**3/(3*EI)
print(f"均布 真值 = {w_exact_udl*1e3:6.4f} mm | 一致 {w_c*1e3:6.4f} | "
      f"集中 {w_l*1e3:6.4f}（高 {(w_l/w_exact_udl-1)*100:.1f} %）")
print(f"端载 真值 = {w_exact_pt*1e3:6.4f} mm | 单元解 {w_p*1e3:6.4f}（精确命中）")
​

均布 真值 = 2.2222 mm | 一致 2.2222 | 集中 2.9630（高 33.3 %）
端载 真值 = 5.9259 mm | 单元解 5.9259（精确命中）
​

判读这份梁记录：一致向量把均布载荷问题的端部挠度直接送到真值——位移是四次曲线，单元内部三次插值处处有误差，唯独端点被变分原理精确锚定，这叫节点的超收敛。集中向量丢了节点力矩信息，梁凭空"软"了三成三，是典型的"载荷搬错、刚度白算"。端部集中载荷下单元解一次命中：三次埃尔米特插值恰好就是真解的形状，网格再粗也无妨。这组记录解释了有限元的两个日常现象：为什么加密网格总能改善分布载荷问题（单元变小，内部误差被压掉），以及为什么节点位移总比单元内部应力先收敛——检查有限元结果，先看节点，再看导出量。

从手解到商用软件

本节的二乘二手解放大成十万自由度的方程组，就是商用有限元的内核。工程上真正要守的只剩三条判读纪律：输入侧查一遍载荷搬运（一致还是集中）、求解侧看收敛（网格加密一档，关键量变化小于百分之几）、输出侧验量纲与平衡（支反力之和等于外载，一票否决）。三条纪律之外，软件是可靠的——它只把 5.2 节的泛函极值重复执行了亿万次。

本节要点回顾

• 有限元是瑞利-里茨的分片版：单元刚度矩阵来自单元势能求导，组装即势能求和。
• 分布载荷按形函数加权搬成一致载荷向量，平分堆放会让梁凭空软三成。
• 杆与梁的节点位移有超收敛性：内部有误差，节点上精确。
• 端部集中载荷下梁单元一次命中真解：三次插值恰是真解形状。
• 有限元判读三纪律：查载荷搬运、看网格收敛、验支反力平衡。
• 应力是位移的导出量：先收敛的是节点位移，后收敛的是应力，验收标准要分档写。

数值台收工，最后回到测量本身——应变花与 DIC 的实验台。