3. 单元刚度矩阵的物理意义与推导

在上一节中，我们已经把一维杆的连续体离散成若干节点与单元，并明确了节点位移向量与节点力向量的组织方式。本节要回答一个更本质的问题：单元刚度矩阵 k^e 中的每一个元素究竟代表什么？它为什么能描述杆件抵抗变形的能力？我们将从能量原理（最小势能原理、虚功原理）和直接法两条路径推导一维两节点杆单元的刚度矩阵，并逐项解读矩阵元素的刚度贡献。

一、先约定符号与坐标系

取一根长度为 L、截面积为 A、弹性模量为 E 的等截面直杆作为一个单元，两端节点编号为 1、2。规定 x 轴沿杆轴由节点 1 指向节点 2，节点位移 u₁、u₂ 与节点力 F₁、F₂ 均以沿 +x 方向为正。全文与图示均采用这一约定，后续不再重复说明。单元内任一点的轴向位移用线性插值表示：

u(x) = N1(x)·u1 + N2(x)·u2
N1(x) = 1 − x/L ,  N2(x) = x/L

形函数 N₁、N₂ 在节点处分别取 1 和 0，保证插值在节点上与真实位移一致。对 u(x) 求导即得轴向应变：

ε = du/dx = (u2 − u1)/L
B = [ −1/L , 1/L ]   （应变-位移矩阵）

这一步非常关键：一维杆的应变在单元内是常数，只取决于两端位移之差。这决定了后面刚度矩阵的简洁形式。

二、基于最小势能原理的推导

单元的总势能由应变能与外力势能两部分组成。应变能密度为 ½·σ·ε = ½·E·ε²，沿体积积分得应变能：

U = ½ ∫V E·ε² dV = ½ E·A·L·[(u2 − u1)/L]²
  = (EA)/(2L)·(u2 − u1)²

外力势能为 −(F₁u₁ + F₂u₂)，于是总势能 Π = U − W。最小势能原理要求 Π 对每个节点位移取驻值，即 ∂Π/∂u₁ = 0、∂Π/∂u₂ = 0。分别求偏导：

∂Π/∂u1 = (EA/L)·(u1 − u2) − F1 = 0
∂Π/∂u2 = (EA/L)·(u2 − u1) − F2 = 0

整理成矩阵形式即得到一维杆单元的刚度方程 k^eu^e = F^e：

k^e = (EA/L) · [  1  −1 ]
                [ −1   1 ]

可以看到，刚度矩阵完全由材料与几何参数 EA/L 决定，结构上高度对称。这正是能量法推导的优雅之处：不需要分析内力分布，只需写出能量表达式并对位移求驻值。

三、基于虚功原理的推导

虚功原理给出另一种等价路径。设单元产生虚位移 δu，对应虚应变 δε = B·δu^e。内力虚功等于应力在虚应变上做的功，外力虚功等于节点力在虚位移上做的功，二者相等：

∫V δε·σ dV = δu^eT · F^e
σ = E·ε = E·B·u^e
⇒ (∫V B^T·E·B dV) · u^e = F^e

括号内的积分就是单元刚度矩阵的一般表达式 k^e = ∫_V B^TEB dV。对一维杆，B 为常数矩阵，积分退化为乘以体积 AL：

k^e = AL · [−1/L, 1/L]^T · E · [−1/L, 1/L]
    = (EA/L) · [ 1  −1 ; −1  1 ]

结果与能量法完全一致。虚功原理的价值在于它给出了可推广到二维、三维单元的统一公式，一维杆只是其中最简单的特例。

四、直接法：从弹簧类比理解刚度贡献

把杆单元想象成一根轴向刚度为 k = EA/L 的弹簧。对弹簧施加单位位移，可以逐列读出刚度矩阵元素。设 u₁ = 1、u₂ = 0，则杆被压缩（或拉伸，取决于方向），需要节点 1 施加 +k 的力、节点 2 施加 −k 的力来维持该位移状态，这正好对应矩阵第一列 [k, −k]^T。同理，令 u₂ = 1、u₁ = 0，得到第二列 [−k, k]^T。

由此可以总结出矩阵元素的物理含义：

• 对角元素 k₁₁ = k₂₂ = EA/L：节点自身单位位移所引起的该节点反力，恒为正，代表杆件抵抗轴向变形的能力。
• 非对角元素 k₁₂ = k₂₁ = −EA/L：一个节点单位位移在另一节点引起的反力，恒为负，体现两端位移的耦合与相互牵制。
• 行和为零：k₁₁ + k₁₂ = 0，说明单元作刚体平移（u₁ = u₂）时不产生任何内力，这是刚度矩阵奇异性（无约束时不可逆）的根源。

五、算例：一个具体单元的数值验证

设某钢杆单元 E = 200 GPa、A = 500 mm²、L = 400 mm。先统一单位到 N 与 mm：E = 2×10⁵ MPa，则 EA/L = 2×10⁵ × 500 / 400 = 2.5×10⁵ N/mm。单元刚度矩阵为：

k^e = 2.5e5 · [  1  −1 ]  N/mm
              [ −1   1 ]

若给定 u₁ = 0、u₂ = 0.8 mm，则节点力 F₁ = 2.5×10⁵×(0 − 0.8) = −2.0×10⁵ N，F₂ = +2.0×10⁵ N。负号表示节点 1 处需要向 −x 方向的约束反力，与杆被拉伸的物理直觉一致。再用材料力学校核：σ = E·ε = 2×10⁵ × (0.8/400) = 400 MPa，轴力 N = σA = 400 × 500 = 2.0×10⁵ N，与矩阵计算结果吻合。这说明刚度矩阵的推导与经典杆件理论是自洽的。

六、刚度矩阵性质与组装前的检查

在进入组装之前，应确认单元刚度矩阵满足以下性质，它们既是理论结论，也是实用的排错依据：

若在编程实现中发现矩阵不对称、行列式非零或行和不为零，通常意味着形函数写错、B 矩阵符号弄反，或 EA/L 的单位未统一。建议在组装前对每个单元做一次上述四项检查，这比在整体方程求解失败后再回头排查要高效得多。