桁架单元

我们来推导简单两节点桁架（杆）单元的著名刚度矩阵。即，考虑两端两个自由度 $q_{1}$ 和 $q_{2}$ 及其相应的力 $P_{1}$ 和 $P_{2}$ 。

源自标准矩阵分析：

力的平衡要求 $P_{1} = - P_{2}$ 。由材料力学可知 $\frac{P_{1}}{A} = E \frac{q_{1} - q_{2}}{L}, \frac{P_{2}}{A} = E \frac{q_{2} - q_{1}}{L},$ 以上方程可写为 $P_{1} = \frac{E A}{L} (q_{1} - q_{2}), P_{2} = \frac{E A}{L} (q_{2} - q_{1})$

根据结构分析，对于等截面面积 A、杨氏模量 E 和长度 L 的杆，其刚度矩阵为： $𝐊 = \frac{E A}{L} [\begin{matrix} 1 & - 1 \\ - 1 & 1 \end{matrix}]$

源自 FEM 第一性原理： 现在，利用 FEM 积分 $𝐊 = \int 𝐁^{𝖳} 𝐄 𝐁 d V$ 来推导。

位移场： 沿杆任意点的轴向位移 u(x) 可通过线性形函数从节点位移 q₁ 和 q₂ 插值得到。若 $q_{1} = 1$ 且 $q_{2} = 0$ ，则 $u (x) = 1 - \frac{x}{L}$ 若 $q_{1} = 0$ 且 $q_{2} = 1$ ，则 $u (x) = \frac{x}{L}$

利用叠加原理，得到 $u (x) = (1 - \frac{x}{L}) q_{1} + (\frac{x}{L}) q_{2} = 𝐍 (x) {\begin{matrix} q_{1} \\ q_{2} \end{matrix}}$ 其中 $𝐍 = ⟨ 1 - \frac{x}{L} \frac{x}{L} ⟩$ 为形函数矩阵。

应变场： 轴向应变 ε 是位移的导数。 $ϵ = \frac{d u}{d x} = \frac{d ٔ 𝐍}{d x} 𝐪 = ⟨ - \frac{1}{L} \frac{1}{L} ⟩ {\begin{matrix} q_{1} \\ q_{2} \end{matrix}}$
应变-位移矩阵 (B)： 由以上可知，对于此简单单元，B 矩阵为常数： $𝐁 = [- \frac{1}{L} \frac{1}{L}]$
材料矩阵 (E)： 对于一维轴向应力，材料矩阵 E 就是标量杨氏模量 *E*。
积分： 现在计算单元体积上的刚度矩阵积分（ $d V = A d x$ ）。 $𝐊 = \int_{0}^{L} 𝐁^{T} E 𝐁 A d x$
由于积分内各项相对于 x 为常数： $𝐊 = \frac{E A}{L^{2}} [\begin{matrix} 1 & - 1 \\ - 1 & 1 \end{matrix}] [x]_{0}^{L} = \frac{E A}{L} [\begin{matrix} 1 & - 1 \\ - 1 & 1 \end{matrix}]$
这成功利用 FEM 基本原理复现了已知结果。

变截面杆分析

考虑一根长度为 L 的杆，一端固定，另一端受集中荷载 P₁。其横截面面积线性变化： $A (x) = A_{0} (1 - \frac{x}{2 L}) .$

1. 解析解

通过沿长度对应变进行积分，可求得“精确”位移。

由于没有分布力，各点的内力 $A (x) E ϵ (x) = A (x) E \frac{d u}{d x}$ 必须为常数，等于轴向力 $P$ 。因此， $A (x) E \frac{d u}{d x} = P$ 从而 $\int_{0}^{L} \frac{d u}{d x} d x = \int_{0}^{L} \frac{P}{E A (x)} d x$

但 $\int_{0}^{L} \underset{d u}{\underset{⏟}{\frac{d u}{d x} d x}} = u (L) - u (0) = q_{2} - q_{1}$ 且 $\int_{0}^{L} \frac{P}{E A (x)} d x = \frac{P}{E} \int_{0}^{L} \frac{1}{1 - \frac{x}{2 L}} d x = \frac{P}{E A_{0}} L \ln 4$ 由于 $P = P_{2} = - P_{1}$ ， ${\begin{matrix} P_{1} \\ P_{2} \end{matrix}} = \frac{E A_{0}}{L \ln 4} [\begin{matrix} 1 & - 1 \\ - 1 & 1 \end{matrix}] {\begin{matrix} q_{1} \\ q_{2} \end{matrix}}$ 因此精确刚度矩阵为 $𝐊 = 0.7213 \frac{E A_{0}}{L} [\begin{matrix} 1 & - 1 \\ - 1 & 1 \end{matrix}]$

2. FEM 解（单个线性单元）

现用单个两节点有限单元对同一根杆建模。采用与桁架示例相同的线性形函数，这意味着 B 矩阵仍为 $𝐁 = \frac{1}{L} [\begin{matrix} 1 & - 1 \end{matrix}]$

关键区别在于面积 A(x) 现在位于刚度积分内：

$𝐊 = \int_{0}^{L} 𝐁^{T} E 𝐁 A (x) d x = 𝐁^{T} E 𝐁 \int_{0}^{L} A (x) d x$

$\int_{0}^{L} A (x) d x = \int_{0}^{L} A_{0} (1 - \frac{x}{2 L}) d x = A_{0} {[x - \frac{x^{2}}{4 L}]}_{0}^{L} = 0.75 A_{0} L$

代回 K 的表达式：

$𝐊 = (\frac{E}{L^{2}} [\begin{matrix} 1 & - 1 \\ - 1 & 1 \end{matrix}]) (0.75 A_{0} L) = \frac{0.75 E A_{0}}{L} [\begin{matrix} 1 & - 1 \\ - 1 & 1 \end{matrix}]$ 此结果为近似值。线性位移场假设（u(x) = Nq）导致恒定应变场（ε = Bq），不能表示变截面杆中真实的变化应变。这种差异导致误差（此处单元刚度过大）。

在此方法中，似乎已将杆替换为具有等于平均横截面面积 $\frac{3}{4} A$ 的等截面杆。误差仅约 4%。

3. 提高 FEM 精度：变截面杆的两单元模型

用两个长度为 L/2 的线性单元对变截面杆建模。可用各单元中点处的值将面积近似为常数。

单元 1（x = 0 到 L/2）： 中点在 x=L/4。 A₁ = A₀(1 - (L/4)/2L) = (7/8)A₀。
单元 2（x = L/2 到 L）： 中点在 x=3L/4。 A₂ = A₀(1 - (3L/4)/2L) = (5/8)A₀。

刚度矩阵为： $𝐊^{(1)} = \frac{E A_{1}}{L / 2} {[\begin{matrix} 1 & - 1 \\ - 1 & 1 \end{matrix}]}_{1, 2} 𝐊^{(2)} = \frac{E A_{2}}{L / 2} {[\begin{matrix} 1 & - 1 \\ - 1 & 1 \end{matrix}]}_{2, 3}$

组装： 通过叠加每个自由度（节点）的贡献，将其组合成 3×3 整体刚度矩阵 K_global。 $𝐊_{g l o b a l} = \frac{2 E}{L} [\begin{matrix} A_{1} & - A_{1} & 0 \\ - A_{1} & A_{1} + A_{2} & - A_{2} \\ 0 & - A_{2} & A_{2} \end{matrix}]$ $𝐊_{g l o b a l} = \frac{E A_{0}}{4} [\begin{matrix} 7 & - 7 & 0 \\ - 7 & 7 + 5 & - 5 \\ 0 & - 5 & 5 \end{matrix}]$

静态凝聚：

通常，我们只关心外部自由度（节点 1 和 3）间的关系，而非内部节点（2）。静态凝聚是一种用于消除内部自由度的矩阵缩减技术。分块整体系统为： $[\begin{matrix} 𝐏_{e} \\ 𝐏_{i} \end{matrix}] = [\begin{matrix} 𝐊_{e e} & 𝐊_{e i} \\ 𝐊_{i e} & 𝐊_{i i} \end{matrix}] [\begin{matrix} 𝐪_{e} \\ 𝐪_{i} \end{matrix}]$ 可写为 $𝐏_{e} = 𝐊_{e e} 𝐪_{e} + 𝐊_{e i} 𝐪_{i}$ $𝐏_{i} = 𝐊_{i e} 𝐪_{e} + 𝐊_{i i} 𝐪_{i}$ 若内部节点无外力（P_i = 0），可求解 q_i 并代回，得到仅关联外部自由度的凝聚刚度矩阵 K_condensed。 $0 = 𝐊_{i e} 𝐪_{e} + 𝐊_{i i} 𝐪_{i}$ $𝐪_{i} = - 𝐊_{i i}^{- 1} 𝐊_{i e} 𝐪_{e}$

$𝐊_{c o n d e n s e d} = 𝐊_{e e} - 𝐊_{e i} 𝐊_{i i}^{- 1} 𝐊_{i e}$

将其应用于两单元模型 $𝐊_{e e} = \frac{E A_{0}}{4} [\begin{matrix} 7 & 0 \\ 0 & 5 \end{matrix}], 𝐊_{e i} = \frac{E A_{0}}{4} [\begin{matrix} - 7 \\ - 5 \end{matrix}]$ $𝐊_{e i} = \frac{E A_{0}}{4} [\begin{matrix} - 7 & - 5 \end{matrix}], 𝐊_{i i} = \frac{12 E A_{0}}{4}$ 得到 2×2 矩阵 $𝐊_{c o n d e n s e d} = 0.7229 \frac{E A_{0}}{L} [\begin{matrix} 1 & - 1 \\ - 1 & 1 \end{matrix}]$

它提供了杆刚度的更精确结果，显著减小了误差（约 0.2%）。

4.2. 高阶单元

我们可使用单个更复杂的单元，而非更多简单单元。例如，二次单元在其中点有第三个节点，并采用二次多项式形函数。

对于 3 节点杆单元（节点位于 x=0, L, L/2），位移场为： $u (x) = ⟨ \begin{matrix} N_{1} (x) & N_{2} (x) & N_{3} (x) \end{matrix} ⟩ {\begin{matrix} q_{1} \\ q_{2} \\ q_{3} \end{matrix}}$ 其中 N₁, N₂, N₃ 为二次函数： $N_{1}^{(2)} (x) = \frac{(x - L / 2) (x - L)}{L^{2} / 2},$ $N_{2}^{(2)} (x) = - \frac{x (x - L)}{L^{2} / 4},$ $N_{3}^{(2)} (x) = \frac{x (x - L / 2)}{L^{2} / 2} .$

这导致应变场 ε(x) 线性变化，对变截面杆是更好的近似。 $ϵ = \frac{d u}{d x} = \underset{𝐁}{\underset{⏟}{⟨ \begin{matrix} \frac{d N_{1}}{d x} & \frac{d N_{2}}{d x} & \frac{d N_{3}}{d x} \end{matrix} ⟩}} {\begin{matrix} q_{1} \\ q_{2} \\ q_{3} \end{matrix}}$ 计算 3×3 刚度矩阵

然后利用静态凝聚得到 2×2 外部刚度矩阵，结果精度很高（示例误差约 0.12%）： $𝐊_{2 \times 2}^{cond} = \frac{13}{18} \frac{E A_{0}}{L} [\begin{matrix} 1 & - 1 \\ - 1 & 1 \end{matrix}] \approx 0.7222 \frac{E A_{0}}{L} [\begin{matrix} 1 & - 1 \\ - 1 & 1 \end{matrix}] .$

若采用线性形函数

则结果与选择两个单元的情形相同。