引言
结构失效并不总是因为应力超过强度。细长受压构件在平均应力远低于屈服极限时，会突然发生侧向弯曲并丧失承载能力，这就是屈曲（失稳）。1907 年加拿大魁北克大桥在悬臂拼装阶段因下弦压杆失稳而整体坍塌，是工程史上由稳定性问题导致灾难的著名案例。与强度问题不同，稳定性研究的是平衡构形本身能否维持：当载荷达到某一临界值时，原有的直线平衡路径发生分岔，结构转向新的弯曲平衡位形，或在缺陷诱导下经极值点突然丧失刚度。本文从 Euler 理想压杆出发，梳理弹性屈曲的临界载荷理论、柔度判据与缺陷修正，然后给出基于有限差分离散的特征值数值解法，并结合工程算例讨论其实用边界。
理论基础
Euler 理想压杆
理想压杆假定：杆件完全平直、载荷严格通过截面形心、材料线弹性、挠度远小于杆长。取两端铰支压杆，挠曲后任意截面的弯矩为 $M = -Pv$（$v$ 为侧向挠度），代入梁的挠曲线近似微分方程 $EI,v” = -M$，得
$$EI \frac{d^2 v}{dx^2} + P v = 0$$
令 $k^2 = P/(EI)$，通解为 $v = C_1 \sin kx + C_2 \cos kx$。代入边界条件 $v(0) = v(L) = 0$，非零解要求 $\sin kL = 0$，即 $kL = n\pi$，于是临界载荷
$$P_{cr} = \frac{n^2 \pi^2 EI}{L^2}, \quad n = 1, 2, 3, \dots$$
$n = 1$ 对应最低阶屈曲模态（半波正弦），即工程上关心的 Euler 临界载荷。对其他支承形式，可引入有效长度系数 $\mu$ 统一表达：
$$P_{cr} = \frac{\pi^2 EI}{(\mu L)^2}$$
$\mu$ 的物理意义是把各种端部约束折算成等效铰支杆长度：悬臂杆 $\mu = 2$，两端固定 $\mu = 0.5$。除以截面面积 $A$ 得临界应力，并定义柔度（长细比）$\lambda = \mu L / i$（$i = \sqrt{I/A}$ 为回转半径）：
$$\sigma_{cr} = \frac{P_{cr}}{A} = \frac{\pi^2 E}{\lambda^2}$$
Euler 公式成立的前提是 $\sigma_{cr}$ 不超过材料比例极限 $\sigma_p$，即 $\lambda \ge \lambda_p = \pi\sqrt{E/\sigma_p}$。以 Q235 钢为例 $\lambda_p \approx 100$，只有大柔度杆才发生弹性屈曲。
弹塑性屈曲、缺陷与稳定判据
中柔度杆进入弹塑性范围，Engesser 提出以切线模量 $E_t = d\sigma/d\varepsilon$ 代替 $E$ 的切线模量理论；Shanley 指出加载与卸载并存时应采用折减模量。真实杆件存在初弯曲、初偏心、残余应力等缺陷，分岔点退化为极值点失稳，临界载荷明显低于 Euler 值，因此各国规范（如 GB 50017 的 a/b/c/d 四类柱子曲线）均以缺陷敏感的经验曲线替代理想解。壳体对缺陷更敏感，Koiter 的初始后屈曲理论给出了缺陷幅值与临界载荷折减的渐近关系。从能量观点看，稳定性由总势能 $\Pi$ 的二阶变分判定：$\delta^2\Pi > 0$ 稳定，$\delta^2\Pi = 0$ 即临界状态，这也是 Rayleigh–Ritz 近似解法的基础。
数值方法与实现
离散化之后，线性屈曲问题归结为广义特征值问题
$$(\mathbf{K} + \lambda \mathbf{K}_g),\boldsymbol{\phi} = \mathbf{0}$$
其中 $\mathbf{K}$ 为弹性刚度矩阵，$\mathbf{K}_g$ 为几何（初应力）刚度矩阵，$\lambda$ 为屈曲载荷因子，$\boldsymbol{\phi}$ 为屈曲模态。对铰支压杆的简化模型 $-EI,v” = Pv$，用二阶中心差分离散后可化为对称矩阵的标准特征值问题，最小特征值即临界载荷。下面给出一个可直接运行的完整算例。
工程案例分析
某平台钢立柱两端铰支，承受轴心压力，截面与材料参数如下表。要求计算其弹性临界载荷并校核是否处于弹性屈曲范围。
| 参数 | 符号 | 数值 | 单位 |
|---|---|---|---|
| 杆长 | $L$ | 3.0 | m |
| 弹性模量 | $E$ | 210 | GPa |
| 截面宽×高 | $b \times h$ | 50×100 | mm |
| 弱轴回转半径 | $i$ | 14.4 | mm |
| 材料屈服强度 | $f_y$ | 235 | MPa |
import numpy as np
from scipy.linalg import eigh
# 铰支-铰支压杆: EI v'' + P v = 0, v(0)=v(L)=0 -> 离散为特征值问题
L = 3.0 # 杆长 m
E = 210e9 # 弹性模量 Pa
b, h = 0.05, 0.10 # 矩形截面宽、高 m
I = b * h**3 / 12 # 弱轴截面二次矩 m^4
N = 200 # 内部节点数
dx = L / (N + 1)
main = 2.0 * np.ones(N); off = -1.0 * np.ones(N - 1)
D2 = (np.diag(main) + np.diag(off, 1) + np.diag(off, -1)) / dx**2
A = E * I * D2 # -EI v'' = P v -> A v = P v
eigvals, eigvecs = eigh(A) # 对称矩阵标准特征值问题
P_num = eigvals[0]
P_euler = np.pi**2 * E * I / L**2
print(f"数值解 P_cr = {P_num/1e3:.1f} kN, Euler 解 = {P_euler/1e3:.1f} kN, "
f"误差 {abs(P_num-P_euler)/P_euler*100:.3f}%")
A_sec = b * h
sigma_cr = P_euler / A_sec
i = np.sqrt(I / A_sec); lam = L / i
print(f"临界应力 {sigma_cr/1e6:.1f} MPa, 柔度 lambda = {lam:.1f}")运行结果：数值解与 Euler 解析解同为 959.5 kN，相对误差仅 0.002%，验证了差分离散的正确性。临界应力 $\sigma_{cr} = 191.9$ MPa，柔度 $\lambda = 103.9$，略大于 Q235 的界限柔度 $\lambda_p \approx 100$，处于弹性屈曲的边缘。值得强调的是，实际设计不能直接采用该值：按 GB 50017 中 b 类截面柱子曲线，$\lambda \approx 104$ 对应的稳定系数仅约 0.5，即考虑初弯曲与残余应力后许用稳定应力约为屈服强度的一半，比 Euler 理想值低得多。对于框架柱，还应注意计算长度系数需按节点约束刚度查取，而非简单套用理论 $\mu$ 值；对于薄壁开口截面，局部屈曲与畸变屈曲往往先于整体屈曲发生，Euler 公式不再适用。
结论与展望
Euler 临界载荷给出了弹性压杆稳定承载力的上限，特征值屈曲分析则是现代有限元软件评估结构稳定性的第一步。但工程结构的实际失稳几乎总是缺陷与材料非线性共同驱动的极值点问题，必须借助几何非线性全过程分析（如 Riks 弧长法）跟踪后屈曲路径，并以规范柱子曲线或直接强度法（DSM）进行设计折减。在数值层面，本文的差分特征值格式可以自然推广到变截面、弹性支承与 Winkler 地基上的压杆，此时解析解不再可得，特征值方法的优势会更加明显。
发表回复