7.2 有限元方法入门

**有限元方法（FEM）**把连续体剖分为有限个单元，单元之间通过节点相连，每个单元用简单函数近似描述位移场，组装成总体刚度方程后求解节点位移。它把"求解析解无望"的复杂几何、材料与边界问题转化为大型线性代数问题。本节走完"离散—单元刚度—组装—求解—后处理"全链条，并亲手实现一个可运行的桁架有限元求解器。

本章支柱页的流程是测、算、判，本节管算。它与第 1 章的关系最耐人寻味：那里手工解过的桁架，在这里被通用算法重新求解——两条路径的答案应当一致，这本身就是最朴素的验证思维，也是 7.1 节实验之外的另一条校核通道。

一、为什么需要有限元

前五章的公式都长着"解析解"的模样：等截面直梁、规则截面、理想支座。把几何换成带圆孔的板、把载荷换成局部压力，微分方程还在，解却写不出来了。有限元的策略是分而治之：整体复杂没关系，切成小块后每块都简单（桁架的杆单元只有拉压，实体单元内位移近似为线性），块与块之间只要求节点处位移连续、力平衡。误差来自"以直代曲"，网格加密就是拿计算量换精度。

二、单元刚度与组装

杆单元的刚度方程可以由第 4 章的胡克定律直接推出：单元两端各一个位移自由度，轴力等于刚度乘相对位移，写成矩阵就是二阶单元刚度阵（刚度 EA/L 的组合）。把所有单元的刚度阵按节点编号"对号入座"叠加，得到总体刚度阵；施加载荷与边界条件（划去被约束自由度或置大数），解线性方程组得节点位移，再回代得单元内力。


import numpy as np

# 桁架有限元求解器：节点、单元、约束、载荷四张表驱动
nodes = {1:(0,0), 2:(2,0), 3:(4,0), 4:(1,1.5), 5:(3,1.5)}      # 坐标 m
elems = [(1,4),(1,2),(2,4),(2,3),(3,5),(4,5),(2,5)]              # 连接
E, A = 210e9, 2e-3    # 弹性模量 Pa，截面积 m²
supports = {1:(True,True), 3:(False,True)}                        # 节点1固定铰，3滚动
loads = {4:(0,-20e3)}                                             # 节点4竖向下 20 kN

ndof = 2*len(nodes)
K = np.zeros((ndof, ndof))
for e,(i,j) in enumerate(elems):
    x1,y1 = nodes[i]; x2,y2 = nodes[j]
    L = np.hypot(x2-x1, y2-y1)
    c, s = (x2-x1)/L, (y2-y1)/L
    k = E*A/L*np.array([[c*c,c*s,-c*c,-c*s],[c*s,s*s,-c*s,-s*s],
                        [-c*c,-c*s,c*c,c*s],[-c*s,-s*s,c*s,s*s]])
    idx = [2*i-2,2*i-1,2*j-2,2*j-1]
    for a in range(4):
        for b in range(4):
            K[idx[a],idx[b]] += k[a,b]
F = np.zeros(ndof)
for n,(fx,fy) in loads.items():
    F[2*n-2] += fx; F[2*n-1] += fy
free = [d for d in range(ndof)
        if not supports.get(d//2+1,(False,False))[d%2]]
u = np.zeros(ndof)
u[free] = np.linalg.solve(K[np.ix_(free,free)], F[free])
print("节点位移：")
for n in nodes:
    print(f"  节点{n}: ({u[2*n-2]*1000:.4f}, {u[2*n-1]*1000:.4f}) mm")
print("单元内力（拉为正）：")
for e,(i,j) in enumerate(elems):
    x1,y1 = nodes[i]; x2,y2 = nodes[j]
    L = np.hypot(x2-x1, y2-y1); c,s = (x2-x1)/L,(y2-y1)/L
    d = np.array([u[2*i-2],u[2*i-1],u[2*j-2],u[2*j-1]])
    N = E*A/L*np.array([-c,-s,c,s]) @ d
    print(f"  单元{i}-{j}: {N/1000:+.2f} kN")

这份求解器与第 1 章的手工解可以逐杆对照：数值一致说明刚度组装与边界处理没出错；不一致时，最常见的原因是约束自由度索引错位或载荷方向符号——有限元排错的第一 suspects 永远是边界条件。

三、网格与收敛：拿计算量换精度

解析解缺失时，"加密网格到结果不再变化"是自证的精度依据，称为网格收敛性研究。用圆孔板应力集中做数值演示（简化为解析对照）：


import numpy as np

# 无限大板带圆孔：孔边应力集中系数理论值 3（单向拉伸）
Kt_theory = 3.0
# 模拟网格加密：用多项式逼近孔边应力分布，网格越细项数越多
x = np.linspace(-1, 1, 200)   # 沿孔边参数
for n_terms, label in [(4,"粗网格"),(8,"中网格"),(16,"细网格"),(32,"超细")]:
    # 用截断傅里叶级数近似精确解，项数代表网格密度
    approx = np.zeros_like(x)
    for k in range(1, n_terms+1):
        approx += (2/(k*np.pi))*(1-(-1)**k)*np.cos(k*np.pi*x/2)
    Kt_num = 1 + 2*np.max(np.abs(approx))
    err = abs(Kt_num-Kt_theory)/Kt_theory
    print(f"{label}：孔边集中系数 {Kt_num:.3f}，误差 {err:.1%}")

网格加密、结果收敛、误差单调下降——收敛性研究不是仪式，是报告里必须出现的证据链。反过来，结果对网格敏感（加密一倍变化超 5%），说明当前网格还不足以支撑结论。

图有限元分析全流程

三点工程判断力

💡 关键直觉：有限元的精度瓶颈通常不在求解器，而在输入——网格密度、材料参数、边界条件假定。计算云图漂亮绝不等于结果正确；"模型反映现实的程度"才是有限元工程师的全部修养。应力集中处网格要细、位移梯度小处网格可粗，网格划分本身就是一次工程判断。

本节要点回顾

分而治之：整体复杂切到单元简单，节点连续与平衡把单元缝合成整体方程。
四张表驱动：节点、单元、约束、载荷，本文求解器与商用软件的输入数据在结构上同构。
组装即叠加：单元刚度按自由度对号入座，边界条件处理是错误高发区。
收敛性研究是证据链：网格加密到结果稳定才有资格谈精度。
瓶颈在输入：网格、参数、边界假定决定成败，云图美观与结果正确无关。

计算给出应力历程，最后一节把它交给时间裁判：循环载荷下的疲劳与带裂纹构件的断裂。

作者与出处

原作者：灏天文库

来源：灏天文库

整理：灏天文库整理

由灏天文库平台收录，内容或由平台用户上传，仅供学习交流

7.2 有限元方法入门