梁单元

梁单元是结构分析中的基本组成部分，旨在模拟主要通过弯曲承受荷载的构件。与仅处理轴向力的桁架单元不同，梁单元可以支持横向荷载，从而产生内部剪力和弯矩。

1. 理论基础：欧拉-伯努利梁理论

标准“经典”梁单元的公式基于欧拉-伯努利梁理论。该理论提供了将单元的变形与其内部应变联系起来的基本运动学关系。该理论建立在几个关键假设之上：

平截面保持平面且垂直：这是最关键的假设。它指出，最初平坦且垂直于梁中心轴的横截面，在梁变形后仍将保持平坦并垂直于该轴。这一假设的一个重要结果是认为剪切变形可以忽略不计。
小位移和小转角：假设梁的挠度和转角很小。这允许进行线性近似，从而极大地简化了数学计算。
线弹性材料：假设梁的材料是各向同性的，并遵循胡克定律，即应力与应变成正比。

基于这些假设，我们推导出了公式推导中最重要的方程。这个方程将梁的物理弯曲与其材料纤维所经历的内部应变联系起来。

为了理解这一点，想象梁弯曲成一个光滑的弧线。沿其长度方向上的任意一点，我们可以用“最佳拟合”圆来近似这段弧线。这个圆的半径称为曲率半径，用 ρ (rho) 表示。

一个非常平缓、轻微的弯曲对应的是一个半径非常大的圆，因此 ρ 值很大。
一个急剧、紧密的弯曲对应的是一个半径很小的圆，因此 ρ 值很小。

现在，考虑距离梁中性轴（位于该圆的中心）距离为 y 的一根纤维上的应变。根据弧的几何形状，轴向应变 (ε_xx) 由下式给出： $ϵ_{x x} = - \frac{y}{ρ}$

这个简单而强大的关系表明，应变在中性轴处 (y=0) 为零，并且随着远离中性轴而线性增加。

为了结构力学中的数学便利性，通常使用一个称为曲率的量，用 κ (kappa) 表示（参见此处）。曲率就是曲率半径的倒数： $κ = \frac{1}{ρ}$ 这意味着平缓的弯曲（大 ρ）具有小的曲率，而急剧的弯曲（小 ρ）具有大的曲率。利用这一定义，我们的应变方程变为： $ϵ_{x x} = - y κ$

最后关键的一步是将这个物理曲率与梁的横向位移函数 v(x) 联系起来。对于欧拉-伯努利理论中假定的小挠度，存在一个直接的数学近似： $κ = \frac{\frac{d^{2} v}{d x^{2}}}{{[1 + {(\frac{d v}{d x})}^{2}]}^{\frac{3}{2}}} \approx \frac{d^{2} v}{d x^{2}}$

将这些关系结合起来，我们得到了最终的、将内部应变与位移场的二阶导数联系起来的本质方程。这个方程构成了我们有限元推导的基础：

$ϵ_{x x} (x, y) = - y \frac{d^{2} v}{d x^{2}}$

2. 单元定义与自由度 (DOFs)

为了捕捉弯曲行为，我们定义一个两节点梁单元。每个节点需要两个自由度 (DOFs) 来描述其状态：

一个横向位移，v。
一个转角，θ。

对于小角度，转角等于位移曲线的斜率，因此 θ = dv/dx。

这意味着一个梁单元总共有四个自由度。这些自由度被收集在单元的广义节点位移向量 q 中：

3. 位移场与埃尔米特形函数的必要性

有限元方法的核心是使用有限数量的节点参数来近似连续的位移场 v(x)。简单的线性函数是不够的，因为它的二阶导数始终为零，这意味着曲率为零，因此弯曲应变为零。

能表示非零、变化曲率的最简单多项式是三次多项式，它有四个未知系数 (a₀ 到 a₃)。这恰好与我们四个节点自由度相匹配。 $v (x) = a_{0} + a_{1} x + a_{2} x^{2} + a_{3} x^{3}$

虽然我们可以使用 a_i 系数，但用形函数来表述问题要有效得多，形函数将位移 v(x) 直接与向量 q 中的节点自由度联系起来。为此使用的特定形函数称为埃尔米特多项式。

与只插值数值的更简单的拉格朗日形函数不同，埃尔米特形函数同时在节点处插值数值及其一阶导数。这正好是我们处理位移 (v) 和转角 (θ = dv/dx) 所需要的。最重要的结果是，埃尔米特多项式保证了相邻单元之间的C¹ 连续性，确保跨节点的斜率连续，变形后的形状是光滑的。

4. 埃尔米特形函数的逐步推导

我们将推导四个形函数，N₁(x) 到 N₄(x)，使得：

每个

N_{i} (x)

都是通过假定

v (x) = a_{0} + a_{1} x + a_{2} x^{2} + a_{3} x^{3}

并求解四个系数得到的，使得

q_{i} = 1

且所有其他自由度 (

q_{j} = 0, j \neq i

) 为零。

我们来推导 N₁(x)，找到一个唯一的、在节点1处对应单位位移 (v₁=1) 而所有其他节点自由度为零的三次多项式。

步骤1：定义 N₁(x) 的节点条件 该函数必须满足： 1. N₁(0) = 1 （节点1处单位位移） 2. N₁'(0) = 0 （节点1处零斜率） 3. N₁(L) = 0 （节点2处零位移） 4. N₁'(L) = 0 （节点2处零斜率）

步骤2：将条件应用于一般三次多项式 令 N₁(x) = a x³ + b x² + c x + d。其导数为 N₁'(x) = 3a x² + 2b x + c。

应用这四个条件： 1. N₁(0) = 1 => a(0) + b(0) + c(0) + d = 1 => d = 1 2. N₁'(0) = 0 => 3a(0) + 2b(0) + c = 0 => c = 0 3. N₁(L) = 0 => aL³ + bL² + 0*L + 1 = 0 => aL³ + bL² = -1 4. N₁'(L) = 0 => 3aL² + 2bL + 0 = 0 => 3aL² + 2bL = 0

步骤3：求解系数 a 和 b 由条件 (4)，我们得到 b = - (3/2)aL。将此代入条件 (3)： aL³ + (-3/2 * aL)L² = -1 aL³ - 3/2 * aL³ = -1 -1/2 * aL³ = -1 => a = 2/L³

现在求 b： b = -3/2 * (2/L³) * L => b = -3/L²

步骤4：组装形函数 N₁(x) 将系数 a, b, c, d 代回多项式形式：

$N_{1} (x) = (\frac{2}{L^{3}}) x^{3} - (\frac{3}{L^{2}}) x^{2} + 1$

对其他三种单位情况遵循相同的过程，我们推导出所有四个埃尔米特形函数。

5. 推导单元刚度矩阵

最后一步是利用这些形函数，从基本积分推导单元刚度矩阵： $𝐊 = \int_{Ω} 𝐁^{𝖳} 𝐄 𝐁 d V$

1. 定义 B 矩阵： 应变为 $ϵ_{x x} = - y \frac{d^{2} v}{d x^{2}}$ 且 $v = ⟨ N_{1} N_{2} N_{3} N_{4} ⟩ 𝐪$ 。曲率 d²v/dx² 为：因此，应变-位移矩阵 B 为：

2. 建立刚度积分： 体积微分为 dV = dA dx。对于一维应力状态，弹性矩阵 E 仅为杨氏模量 E。 $𝐊 = \int_{0}^{L} \int_{A} (𝐁)_{4 \times 1}^{𝖳} E (𝐁)_{1 \times 4} d A d x$ 代入 B 的表达式：不依赖于面积 (A) 的项可以移到内积分之外：

3. 最终结果： 我们认识到 $\int_{A} y^{2} d A$ 是截面二次矩 I 的定义。这极大地简化了积分： $𝐊 = \int_{0}^{L} {[\frac{d^{2} 𝐍}{d x^{2}}]}^{𝖳} E I [\frac{d^{2} 𝐍}{d x^{2}}] d x$ 计算四个埃尔米特形函数的二阶导数，形成矩阵乘积，并对 x 从 0 到 L 进行积分，最终得到经典的4x4 欧拉-伯努利梁刚度矩阵：

* 有时，我们会加一个下标 e（如 $𝐊_{e}$ ）来表示上述矩阵是单元的刚度矩阵。

任意方向梁与变换矩阵

如果我们考虑6个自由度（考虑梁柱单元），梁也可以被拉伸或压缩，那么单元的刚度矩阵为 $𝐊_{e}^{local} = [\begin{matrix} \frac{E A}{L} & 0 & 0 & - \frac{E A}{L} & 0 & 0 \\ 0 & \frac{12 E I}{L^{3}} & \frac{6 E I}{L^{2}} & 0 & - \frac{12 E I}{L^{3}} & \frac{6 E I}{L^{2}} \\ 0 & \frac{6 E I}{L^{2}} & \frac{4 E I}{L} & 0 & - \frac{6 E I}{L^{2}} & \frac{2 E I}{L} \\ - \frac{E A}{L} & 0 & 0 & \frac{E A}{L} & 0 & 0 \\ 0 & - \frac{12 E I}{L^{3}} & - \frac{6 E I}{L^{2}} & 0 & \frac{12 E I}{L^{3}} & - \frac{6 E I}{L^{2}} \\ 0 & \frac{6 E I}{L^{2}} & \frac{2 E I}{L} & 0 & - \frac{6 E I}{L^{2}} & \frac{4 E I}{L} \end{matrix}]$

现在考虑一根与全局 $x$ 轴成 $α$ 角的梁。在这种情况下，我们必须将全局坐标系 (x, y) 旋转角度 $α$ 以到达局部坐标系 (x', y')。

自由度 $u$ 和 $v$ 的变换与桁架单元中 u 和 v 的变换相同。由于自由度 $θ$ 在两个坐标系中相同，变换矩阵变为

$𝐓 = [\begin{matrix} c & s & 0 & 0 & 0 & 0 \\ - s & c & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & c & s & 0 \\ 0 & 0 & 0 & - s & c & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 \end{matrix}]$ 其中 $c = \cos α, s = \sin α .$ 全局坐标系下的刚度矩阵变为 $𝐊_{e}^{global} = 𝐓^{𝖳} 𝐊_{e}^{local} 𝐓$