[FEM-3]有限元分析的基本流程

---

1 问题描述

下图是一个一维三连杆结构，相关的材料参量和尺寸为

$$\begin{array}{c}\begin{array}{r}{E}_1=E_2=E_3=2\times 10^5\mathrm{Pa}3{\mathrm{A}}_1=2{\mathrm{A}}_2={\mathrm{A}}_3=0.06{\mathrm{m}}^2{\mathrm{l}}_1={\mathrm{l}}_2={\mathrm{l}}_3=0.1\mathrm{m}\end{array}\end{array}$$

图1 1D三连杆结构的受力图示

1.1 节点编号及单元划分

在杆件连接处划分节点，对应的节点和单元编号(numbering of element)如图1所示。将每一个单元分离出来，并标出每一个节点的位移和外力，如图2所示。此处位移和力的方向均以X轴正方向来标注。

图2 各个单元的节点位移和外力

1.2 计算各个单元的单元刚度方程

（1）单元①的刚度方程为

$$\begin{array}{c}\left[\begin{array}{cc}\frac{E_1{A}_1}{l_1}& -\frac{E_1{A}_1}{l_1}\\ -\frac{E_1{A}_1}{l_1}& \frac{E_1{A}_1}{l_1}\end{array}\right]\left[\begin{array}{c}{u}_1\\ u_2\end{array}\right]=\left[\begin{array}{c}{P}_1^{(1)}\\ P_2^{(1)}\end{array}\right]\end{array}$$

（2）单元②的刚度方程为

$$\begin{array}{c}\left[\begin{array}{cc}\frac{E_2{A}_2}{l_2}& -\frac{E_2{A}_2}{l_2}\\ -\frac{E_2{A}_2}{l_2}& \frac{E_2{A}_2}{l_2}\end{array}\right]\left[\begin{array}{c}{u}_2\\ u_3\end{array}\right]=\left[\begin{array}{c}{P}_2^{(2)}\\ P_3^{(2)}\end{array}\right]\end{array}$$

（3）单元③的刚度方程为

$$\begin{array}{c}\left[\begin{array}{cc}\frac{E_3{A}_3}{l_3}& -\frac{E_3{A}_3}{l_3}\\ -\frac{E_3{A}_3}{l_3}& \frac{E_3{A}_3}{l_3}\end{array}\right]\left[\begin{array}{c}{u}_3\\ u_4\end{array}\right]=\left[\begin{array}{c}{P}_3^{(3)}\\ P_4^{(3)}\end{array}\right]\end{array}$$

1.3 组装各单元刚度方程

由于整体结构是由各个单元按一定连接关系组合而成的，因此，需要按照节点的对应位置将式(2)、式(3)、式(4)进行组装，以形成整体刚度方程

$$\begin{array}{c}\left[\begin{array}{cccc}\frac{E_1{A}_1}{l_1}& -\frac{E_1{A}_1}{l_1}& 0& 0\\ -\frac{E_1{A}_1}{l_1}& \frac{E_1{A}_1}{l_1}+\frac{E_2{A}_2}{l_2}& -\frac{E_2{A}_2}{l_2}& 0\\ 0& -\frac{E_2{A}_2}{l_2}& \frac{E_2{A}_2}{l_2}+\frac{E_3{A}_3}{l_3}& -\frac{E_3{A}_3}{l_3}\\ 0& 0& -\frac{E_3{A}_3}{l_3}& \frac{E_3{A}_3}{l_3}\end{array}\right]\left[\begin{array}{c}{u}_1\\ u_2\\ u_3\\ u_4\end{array}\right]=\left[\begin{array}{c}{P}_1^{(1)}\\ P_2^{(1)}+P_2^{(2)}\\ P_3^{(2)}+P_3^{(3)}\\ P_4^{(3)}\end{array}\right]\end{array}$$

式(5)中的$P_1^{(1)}$、$P_2^{(1)}+P_2^{(2)}$、$P_3^{(2)}+P_3^{(3)}$、$P_4^{(3)}$就是节点1、2、3、4上的合成节点力，即

$$\begin{array}{c}{P}_1=P_1^{(1)}{P}_2=P_2^{(1)}+P_2^{(2)}{P}_3=P_3^{(2)}+P_3^{(3)}{P}_4=P_4^{(3)}\end{array}$$

由图1可知，$P_1=-100\mathrm{N}$、$P_2=0\mathrm{N}$、$P_3=50\mathrm{N}$，$P_4$为支座的支反力。将该结构的材料参数和几何尺寸带入式(5)，可得

$$\begin{array}{c}\left[\begin{array}{cccc}4\times 10^4& -4\times 10^4& 0& 0\\ -4\times 10^4& 1\times 10^5& -6\times 10^4& 0\\ 0& -6\times 10^4& 1.8\times 10^5& -1.2\times 10^5\\ 0& 0& -1.2\times 10^5& 1.2\times 10^5\end{array}\right]\left[\begin{array}{c}{u}_1\\ u_2\\ u_3\\ u_4\end{array}\right]=\left[\begin{array}{c}-100\\ 0\\ 50\\ P_4\end{array}\right]\end{array}$$

1.4 处理边界条件并求解

图1所示结构的位移边界条件为：$u_4=0$，代入上式可得

$$\begin{array}{c}\left[\begin{array}{cccc}4\times 10^4& -4\times 10^4& 0& 0\\ -4\times 10^4& 1\times 10^5& -6\times 10^4& 0\\ 0& -6\times 10^4& 1.8\times 10^5& -1.2\times 10^5\\ 0& 0& -1.2\times 10^5& 1.2\times 10^5\end{array}\right]\left[\begin{array}{c}{u}_1\\ u_2\\ u_3\\ 0\end{array}\right]=\left[\begin{array}{c}-100\\ 0\\ 50\\ P_4\end{array}\right]\end{array}$$

划去上述刚度矩阵的第4例和第4行，则有

$$\begin{array}{c}\left[\begin{array}{ccc}4\times 10^4& -4\times 10^4& 0\\ -4\times 10^4& 1\times 10^5& -6\times 10^4\\ 0& -6\times 10^4& 1.8\times 10^5\end{array}\right]\left[\begin{array}{c}{u}_1\\ u_2\\ u_3\end{array}\right]=\left[\begin{array}{c}-100\\ 0\\ 50\end{array}\right]\end{array}$$

求解得到

$$\begin{array}{c}\begin{array}{r}\begin{array}{cc}{u}_1& =-4.58333\times 10^{-3}\mathrm{m}\\ u_2& =-2.08333\times 10^{-3}\mathrm{m}\\ u_3& =-4.16667\times 10^{-4}\mathrm{m}\end{array}\}\end{array}\end{array}$$

1.5 求支反力

求解得到支反力$P_4$为

$$\begin{array}{c}{P}_4=-1.2\times 10^5\times u_3=50N\end{array}$$

1.6 求各个单元的其它力学量（应变、应力）

各个单元的应变为

$$\begin{array}{c}\begin{array}{rl}{\epsilon }^{(1)}& =\frac{u_2-u_1}{l_1}=\frac{(-2.08333+4.58333)\times 10^{-3}}{0.1}=2.49997\times 10^{-2}\mathrm{m}\\ {\epsilon }^{(2)}& =\frac{u_3-u_2}{l_2}=\frac{(-0.416667+2.08333)\times 10^{-3}}{0.1}=1.6667\times 10^{-2}\mathrm{m}\\ {\epsilon }^{(3)}& =\frac{u_4-u_3}{l_3}=\frac{(0+4.16667)\times 10^{-4}}{0.1}=4.16667\times 10^{-3}\mathrm{m}\end{array}\end{array}$$

各个单元的应力为

$$\begin{array}{c}\begin{array}{rl}{\sigma }^{(1)}& =E_1{\epsilon }^{(1)}=4.999\times 10^3\mathrm{Pa}\\ {\sigma }^{(2)}& =E_2{\epsilon }^{(2)}=3.3333\times 10^3\mathrm{Pa}\\ {\sigma }^{(3)}& =E_3{\epsilon }^{(3)}=8.3333\times 10^2\mathrm{Pa}\end{array}\end{array}$$

这样可以得到一种直观的有限元分析思路：

• 复杂的几何和受力对象划分为一个一个形状比较简单的标准“构件”，称为单元
• 给出单元节点的位移和受力描述，构件单元刚度矩阵
• 通过单元与单元之间的节点连接关系进行单元的组装，得到结构的整体刚度矩阵
• 根据位移约束和受力状态，处理边界条件，并进行求解

图3 有限元分析的基本流程图示

参考

1. 曾攀. 有限元分析基础教程[M]. 北京: 清华大学出版社, 2008.