有限元分析及运用课程笔记第二章：基于直接刚度法的杆系有限元方法

本笔记主要来源为清华大学公开课《有限元分析及应用》（曾攀、雷丽萍)。非常感谢两位老师深入浅出、循序渐进的讲解，部分图片、公式也出自该公开课。

我的专业是自动化控制，在空余时间学习这门课程并记录课程笔记只是出于兴趣。
我的评注内容都放在这种的引用块内，方便浏览时理解，如果存在错误欢迎与各位探讨。
笔记只供快速浏览或复习，实际学习建议观看原视频或参考相关资料，机械工程、力学有些问题问我也不会。

2.1 弹簧的力学分析原理

弹簧力和位移的关系:F=kδF=k\deltaF=kδ
其中，kkk为弹簧的刚度系数,FFF为作用在端点的力，δ\deltaδ为端点位移量

更一般地，对于弹簧整体受力分析有：

{k(u2−u1)=F2F1+F2=0\left\{ \begin{aligned} k(u_2-u_1) = F_2 \\ F_1+F_2 = 0 \end{aligned} \right. {k(u2−u1)=F2F1+F2=0
变换得到：
{k(u1−u2)=F1k(u2−u1)=F2(1)\left\{ \begin{aligned} k(u_1-u_2) = F_1 \\ k(u_2-u_1) = F_2 \end{aligned} \right. \tag{1} {k(u1−u2)=F1k(u2−u1)=F2(1)

我们可以将公式(1)改写为矩阵形式（弹簧的平衡方程）:

[k−k−kk][u1u2]=[F1F2](2)\begin{bmatrix} k & -k \\ -k & k \end{bmatrix} \begin{bmatrix} u_1 \\ u_2 \end{bmatrix} = \begin{bmatrix} F_1 \\ F_2 \end{bmatrix} \tag{2} [k−k−kk][u1u2]=[F1F2](2)

同理，对于组装的弹簧在B点施加力F2F_2F2，我们有:

[k1−k1−k1k1][u1u2]=[F1F2−F12](3)\begin{bmatrix} k_1 & -k_1 \\ -k_1 & k_1 \end{bmatrix} \begin{bmatrix} u_1 \\ u_2 \end{bmatrix} = \begin{bmatrix} F_1 \\ F_2 - F_{12} \end{bmatrix} \tag{3} [k1−k1−k1k1][u1u2]=[F1F2−F12](3)

[k2−k2−k2k2][u2u3]=[F2+F12F3](4)\begin{bmatrix} k_2 & -k_2 \\ -k_2 & k_2 \end{bmatrix} \begin{bmatrix} u_2 \\ u_3 \end{bmatrix} = \begin{bmatrix} F_2 + F_{12} \\ F_3 \end{bmatrix} \tag{4} [k2−k2−k2k2][u2u3]=[F2+F12F3](4)

这里F12F_{12}F12指的是弹簧之间的相互作用力

将公式（3）与公式（4）相加可以得到:

[k1−k10−k1k1+k2−k20−k2k2][u1u2u3]=[F1F2F3](5)\begin{bmatrix} k_1 & - k_1 & 0 \\ - k_1 & k_1 + k_2 & -k_2 \\ 0 & -k_2 & k_2 \end{bmatrix} \begin{bmatrix} u_1 \\ u_2 \\ u_3 \end{bmatrix}= \begin{bmatrix} F_1 \\ F_2 \\ F_3 \end{bmatrix} \tag{5} ⎣⎡k1−k10−k1k1+k2−k20−k2k2⎦⎤⎣⎡u1u2u3⎦⎤=⎣⎡F1F2F3⎦⎤(5)
由于端点固定，u1=0,u3=0u_1 = 0,u_3 = 0u1=0,u3=0,由此可以得到：
(k1+k2)u2=F2(k_1+k_2)u_2=F_2 (k1+k2)u2=F2
我们就可以得到u2u_2u2了，并由此可以将F1、F2F_1、F_2F1、F2求出。

2.2 弹簧单元与杆单元的比较

我们知道杠单元的平衡关系为：
平衡方程：σ(x)=F/A平衡方程：\sigma(x) = F/A 平衡方程：σ(x)=F/A
物理方程：σ=Eε物理方程：\sigma = E \varepsilon 物理方程：σ=Eε
几何方程：ε(x)=δ/l几何方程：\varepsilon(x) = \delta/l 几何方程：ε(x)=δ/l
可以的得到:
F=AEl⋅δ(6)F = \frac{AE}{l} \cdot \delta \tag{6} F=lAE⋅δ(6)

令 k=AElk=\frac{AE}{l}k=lAE 就与弹簧关系一致！

2.3 杠单元的坐标变换

上面已经得到了局部坐标系下的杆单元描述
KeUe=Fe(7)K^eU^e=F^e\tag{7}KeUe=Fe(7)
其中
Ke=AEl[1−1−11]K^e=\frac{AE}{l}\begin{bmatrix}1&-1\\-1&1\end{bmatrix}Ke=lAE[1−1−11]
Ue=[u1u2]U^e=\begin{bmatrix}u_1\\u_2\end{bmatrix}Ue=[u1u2]
Fe=[F1F2]F^e=\begin{bmatrix}F_1\\F_2\end{bmatrix}Fe=[F1F2]

这里的 eee 是 element 的缩写，只是标记这是在局部坐标系，不是 eee 次方的意思

如果要在世界坐标系下处理，只要变换到世界坐标系下即可！

假定从世界坐标变换为局部坐标的坐标变换矩阵为 Tw→eT^{w \rightarrow e}Tw→e，那么有：
Ue=Tw→eUw(8)U^e = T^{w \rightarrow e} U^ w \tag{8} Ue=Tw→eUw(8)

我这里使用 www 标记为世界坐标，避免和 eee 混淆所以与课程标注不一致，但是推导思路和逻辑会完全保持一致

将 (8) 式代入 (7) 式，我们可以得到:

KeTw→eUw=FeK^e T^{w \rightarrow e} U^w = F^e KeTw→eUw=Fe

为了得到世界坐标下的力，我们只要左右同时乘以 Tw→e⊤{T^{w \rightarrow e}}^\topTw→e⊤（其实也就是Te→wT^{e \rightarrow w}Te→w）,有
Tw→e⊤KeTw→e‾Uw=Tw→e⊤Fe‾KwUw=Fw(9)\underline{{T^{w \rightarrow e}}^\top K^e T^{w \rightarrow e}} U^w = \underline{{T^{w \rightarrow e}}^\top F^e} \\ K^w U^w = F ^ w \tag{9} Tw→e⊤KeTw→eUw=Tw→e⊤FeKwUw=Fw(9)

在世界坐标系下单元刚度矩阵为：
Kw=Tw→e⊤KeTw→eK^w = {T^{w \rightarrow e}}^\top K^e T^{w \rightarrow e} Kw=Tw→e⊤KeTw→e

世界坐标系下的节点力矩阵为：
Fw=Tw→e⊤FeF^w = {T^{w \rightarrow e}}^\top F^e Fw=Tw→e⊤Fe

我认为这里的变换其实只要得到 Tw→eT^{w \rightarrow e}Tw→e 就可以求到所有矩阵的实际值了。由于局部坐标系中杠的自由度为2(u1、u2u_1、u_2u1、u2)，三维中杠的自由度为6（两个端点，每个端点三个坐标值），所以 Tw→eT^{w \rightarrow e}Tw→e 其实是一个 2×62 \times 62×6 的矩阵。

2.4 一个四杆结构的实例分析

题目如下：

解题思路为：

结构离散化与标号

一共有四个单元，四个节点，确定坐标轴，确定节点的世界坐标，得到每个单元的世界位置矩阵 UwU^wUw
各个单元矩阵描述

根据 Kw=Tw→e⊤KeTw→eK^w = {T^{w \rightarrow e}}^\top K^e T^{w \rightarrow e}Kw=Tw→e⊤KeTw→e ，得到每个单元的刚度矩阵,需要注意刚度矩阵中节点对应关系（即方程（9）KwUw=FwK^w U^w = F ^ wKwUw=Fw 的每项的关系）
建立世界坐标系刚度方程

类似方程（5）的获取方式按照对应节点编号进行组装，根据输入和约束可以确定节点力矩阵 FwF^wFw
求解刚度方程

我们已经得到了刚度方程的所有项，求解矩阵中的未知数。

2.5 软件模拟求解

略过