📚 课程目录第1章杆单元轴向拉伸有限元建模与位移求解

1. 杆单元轴向拉伸有限元建模与位移求解

取一根等截面直杆，左端固定、右端受轴向拉力。设杆长 L = 500 mm，截面积 A = 200 mm²，弹性模量 E = 210 GPa，右端集中力 F = 42 kN。杆件只承受轴向载荷，可退化为一维问题。

由材料力学，轴力沿杆长为常数 N = F = 42 kN，轴向应力 σ = N/A = 210 MPa，轴向应变 ε = σ/E = 1.0×10⁻³，右端位移 u = εL = 0.5 mm。这组结果将作为有限元解的对照基准。

一维杆单元有两个节点，每个节点一个轴向自由度。单元刚度矩阵由 k = EA/L 决定：

k_e = (EA/L) * [ 1  -1
                -1   1 ]

将杆划分为 2 个等长单元，共 3 个节点，节点坐标 x = [0, 250, 500] mm。每个单元长度 L_e = 250 mm，单元刚度系数：

EA/L_e = (210000 MPa × 200 mm²) / 250 mm
       = 168000 N/mm = 1.68×10⁵ N/mm

单元①连接节点1-2，单元②连接节点2-3，两者刚度矩阵形式相同，只是装配位置不同。

整体自由度向量为 U = [u₁, u₂, u₃]ᵀ。将两个单元矩阵按节点编号叠加，得到整体刚度矩阵：

K = (EA/L_e) * [ 1  -1   0
                -1   2  -1
                 0  -1   1 ]

代入 EA/L_e = 1.68×10⁵ N/mm，整体方程 K·U = P，其中载荷向量 P = [R₁, 0, 42000]ᵀ N，R₁ 为左端支反力。

左端固定：u₁ = 0。划去第1行第1列，剩余 2×2 方程组：

1.68×10⁵ * [ 2  -1 ] [u₂]   [    0 ]
            [-1   1 ] [u₃] = [42000 ]

由第二式：1.68×10⁵(−u₂ + u₃) = 42000，得 u₃ − u₂ = 0.25 mm。
由第一式：1.68×10⁵(2u₂ − u₃) = 0，得 u₃ = 2u₂。
联立解得 u₂ = 0.25 mm，u₃ = 0.50 mm。

与解析解 u(250)=0.25 mm、u(500)=0.50 mm 完全一致，说明两单元线性杆模型对常轴力问题是精确的。

单元轴力由 N = (EA/L_e)(u_j − u_i) 计算：

单元	节点	Δu (mm)	轴力 N (kN)	应力 σ (MPa)
①	1→2	0.25	42.0	210
②	2→3	0.25	42.0	210

左端支反力 R₁ = −42 kN，与外力平衡。轴力、应力、位移三项均与解析解一致，验证通过。

杆单元对常轴力问题是精确的：只要载荷只作用在节点上、截面与材料沿杆不变，任意单元划分都会得到相同的节点位移和单元轴力。若存在分布轴向载荷或变截面，则需加密网格或改用带一致载荷向量的单元。