桥为什么会那样下垂 — 桁架有限元的直接刚度法

1956年，波音的 M. J. Turner 用手算后掠翼应力，算到一半卡住了。构件有几百根，平衡方程没法一根一根解。他和同事给出的办法是直接刚度法。先做单元刚度，再在自由度重叠的位置相加，然后解一次 $Ku = f$。

先看角度。杆放平，竖直那一格就死了。把角度滑条拉到0，4×4刚度里竖直自由度的行和列会一起变成0。

看 $\theta$ 为0时哪些格子空掉。杆只沿着自己的轴硬。角度抬高以后，那些0会拆进 $x$ 和 $y$。

刚度从一根弹簧说起

杆单元是只沿轴向受力的弹簧。长度 $L$、截面积 $A$、弹性模量 $E$，轴向刚度就是 $k = EA/L$。每拉长一单位长度所需要的力。两端轴向位移 $u_1, u_2$ 对应的节点力是：

$f_1, f_2$ 是两端节点力，括号里是局部刚度矩阵 $k_e$。行和为0。两个节点一起挪的刚体运动不需要力。所以拿掉边界条件，矩阵就是奇异的。

从局部到整体 — 旋转变换

构件朝各个方向歪着。局部轴向位移得写成整体的 $x, y$，才组得起来。设构件与 $x$ 轴的夹角为 $\theta$，$c = \cos\theta$，$s = \sin\theta$。轴向位移是整体位移的投影，$u_{\text{axial}} = c\,u_x + s\,u_y$。这个变换从两边乘上去，就得到4×4的整体单元刚度。

每一项连接节点1的 $x, y$ 和节点2的 $x, y$。$\theta = 0$ 时 $s = 0$、$c = 1$，上面的矩阵乘上 $EA/L$ 之后，第一行是1, 0, −1, 0，第二行全是0。第三行是第一行变号，第四行也是0。2×2块秩为1，就是这些0造成的。坐标一歪，刚度张量跟着转。那种转动写在用度量把本构张量转过去里。

把重叠的自由度相加 — 整体组装

组装就是在重叠的自由度上把刚度加上。节点 $i$ 的自由度放在整体指标 $2i, 2i+1$ 上，先 $x$ 后 $y$。把单元的4个局部自由度送到这些指标上，再把单元刚度累进去。

$\mathbf{L}_e$ 是把单元自由度送到整体自由度的选择矩阵，$K$ 是整体刚度矩阵。代码并不去乘选择矩阵，而是按指标直接加。一个节点上接了好几根构件，贡献就堆在那个对角块里。

边界条件：抹去支座

刚组装完的 $K$ 是奇异的。结构还飘着，刚体运动没被挡住。支座上的位移必须收成0，才解得开。把自由度分成自由集合 $f$ 和固定集合 $c$。

$u_c = 0$ 时只留下上面一行，$K_{ff}\,u_f = F_f$。这个缩减系统非奇异。反力事后用 $R_c = K_{cf}\,u_f$ 收回来。$R_c$ 是支座还给结构的力。

把荷载滑条和EA一起调高，端部下垂会怎么减小，画面上直接看得到。

EA变大，同样的荷载下垂得更少。刚度矩阵就是从荷载到位移的那个映射。荷载不动，只把EA加倍，看端部位移是不是往一半走。

Python — 求解12杆桁架#

墙边固住的三跨悬臂。8个节点，12根杆。坐标、连接、荷载进去，位移和轴力出来。

import numpy as np
 
nodes = np.array([[0,0],[0,1],[1,0],[1,1],[2,0],[2,1],[3,0],[3,1]], float)
members = [(0,2),(2,4),(4,6),(1,3),(3,5),(5,7),
           (2,3),(4,5),(6,7),(1,2),(3,4),(5,6)]
EA = 8.0e6           # 轴向刚度 EA [N]
fixed = [0, 1]       # 固定在墙上的节点
 
def bar_stiffness(p1, p2, EA):
    d = p2 - p1
    L = np.hypot(*d)
    c, s = d / L
    k = EA / 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]])
    return k, L, (c, s)
 
ndof = nodes.shape[0] * 2
K = np.zeros((ndof, ndof))
geom = []
for a, b in members:                       # 组装整体刚度
    k, L, cs = bar_stiffness(nodes[a], nodes[b], EA)
    dof = [2*a, 2*a+1, 2*b, 2*b+1]
    K[np.ix_(dof, dof)] += k
    geom.append((L, cs))
 
F = np.zeros(ndof)                          # 自由端（节点6,7）施加24 kN荷载
for n in (6, 7):
    F[2*n+1] -= 12.0e3
 
fixed_dof = [d for n in fixed for d in (2*n, 2*n+1)]
free_dof = [d for d in range(ndof) if d not in fixed_dof]
 
u = np.zeros(ndof)                          # 缩减系统 K_ff u_f = F_f
u[free_dof] = np.linalg.solve(K[np.ix_(free_dof, free_dof)],
                              F[free_dof])
 
for m, (a, b) in enumerate(members):        # 杆件轴力 N = (EA/L)[-c,-s,c,s]·u
    L, (c, s) = geom[m]
    dof = [2*a, 2*a+1, 2*b, 2*b+1]
    N = EA / L * np.array([-c, -s, c, s]) @ u[dof]
    print(f"member {a}-{b}: N = {N/1e3:+7.2f} kN "
          f"({'tension' if N > 0 else 'compression'})")
 
print(f"free-end drop = {u[2*6+1]*1e3:.3f} mm")

np.ix_ 把单元刚度加进整体的格子。这个循环和商业FEA的组装是同一种结构。

节点间距是1 m。9根水平杆和竖杆的长度是 $L = 1$。代码里的EA是 $8.0\times 10^6$ N，于是 $k = EA/L = 8.0\times 10^6$ N/m。端部一个节点的竖向荷载是 $12\times 10^3$ N，所以只推这一根杆，位移是 $1.5\times 10^{-3}$ m，也就是1.5 mm。桁架端部不是这个除法。杆件把荷载分开扛，真正求解的是缩减后的 $K_{ff}$。

斜杆有三根。连接(0, 1)和(1, 0)的那根长度是 $\sqrt{2}$，刚度是水平杆的 $1/\sqrt{2}$。EA相同，斜杆更软。这根杆的角度是45°。$c^2$ 和 $s^2$ 都是 $1/2$，$cs$ 的绝对值也是 $1/2$。4×4里不为0的格子全是 $EA/L$ 的一半，长度是 $\sqrt{2}$，格子的大小就是 $EA/(2\sqrt{2})$。水平杆的竖直格子空着为0，斜杆则四个方向都按这个大小填满。同一个EA，只靠角度就分成这样。8个节点有16个自由度。墙上的节点0和1约束了其中4个，所以求解的是12×12。二维至少要约束3个自由度，这个例子约束了4个。墙上一个节点的 $y$ 若放开，个数就贴着下限。两个都放开，$K_{ff}$ 会再次奇异。

当刚度矩阵变奇异时

现场最常见的失败就是奇异矩阵。分叉有三条。

边界条件漏了，刚体运动还在，$K_{ff}$ 就是奇异的。二维至少约束3个自由度，三维至少6个。

没有分成三角形的四边形节间，杆不够，是机构。桁架得用三角形填满，不然会晃。

长度是0，或者节点叠在一起，$L = 0$ 会让除法炸掉。坐标相同的点，合并之前先滤掉。

少固定一个墙节点，组装就停在这里

水平杆的1.5 mm说的是单根杆。12杆悬臂的端部位移要解 $K_{ff}\,u_f = F_f$ 才有，墙上那4个自由度不拿掉，这个式子就不成立。支承少一个，刚体模态还在，位移就不是一个数，而是零空间。

杆换成面单元之后，同一种组装会因为别的理由变硬。薄壳上那种锁死怎么出现，写在剪切闭锁，以及MITC怎样把刚度还回来里。桁架上先核对4个支承。只有模型变薄了，再转到那一页。