[FEM-5]局部坐标系中的杆单元描述

---

单元的描述包括单元的几何和节点描述、位移场、应变场、应力场、势能，也就是要充分描述问题的三大类变量(位移、应变、应力)以及三大类方程(平衡方程、几何方程、物理方程)来计算单元的势能，然后通过最小势能原理或虚功原理来得到单元的方程，

实际上，单元内位移场的描述就是其试函数的选取

1 单元的几何和节点描述

图1 局部坐标系中的杆单元

上图为一个在局部坐标系(local coordinate system)中的杆单元(bar element)，由于有两个端节点(Node 1和Node 2)，则基本变量为节点位移(向量)矩阵 ${\mathbf{q}}^e$

$$\begin{array}{c}\begin{array}{r}{\mathbf{q}}^e=\left[\begin{array}{c}{u}_1\\ u_2\end{array}\right]={\left[\begin{array}{cc}{u}_1& u_2\end{array}\right]}^T\end{array}\end{array}$$

将每一个描述物体位置状态的独立变量叫做一个自由度(DOF, degree of freedom)

上述的节点位移为两个自由度。节点力(向量)列阵 ${\mathbf{P}}^e$

$$\begin{array}{c}\begin{array}{r}{\mathbf{P}}^e=\left[\begin{array}{c}{P}_1\\ P_2\end{array}\right]={\left[\begin{array}{cc}{P}_1& P_2\end{array}\right]}^T\end{array}\end{array}$$

若该单元承受有沿轴向的分布载荷，可以将其等效到节点上。

利用函数插值、几何方程、物理方程以及势能计算公式，可以将单元的所有力学参量(即场变量：$u(x)$、$\epsilon (x)$、$\sigma (x)$和${\prod }^e$)用节点位移列阵${\mathbf{q}}^e$及相关的插值函数表示。

2 单元位移场的表达

设该单元的位移场为$u(x)$，由Taylor级数，其可以表示为

$$\begin{array}{c}\begin{array}{r}u(x)=a_0+a_{1}x+a_2{x}^2+\cdots \end{array}\end{array}$$

该函数将由两个端节点的位移$u_1$和$u_2$来进行插值确定，可取式(3)的前两项来作为该单元的位移插值模式(interpolation model)

$$\begin{array}{c}\begin{array}{r}u(x)=a_0+a_{1}x\end{array}\end{array}$$

式中：$a_0$和$a_1$为待定系数(unknows)。

单元节点条件为

$$\begin{array}{c}\begin{array}{rl}u(x){|}_{x=0}& =u_1\\ u(x){|}_{x=l^e}& =u_2\end{array}\end{array}$$

代入式(4)中可求得$a_0$和$a_1$为

$$\begin{array}{c}\begin{array}{rl}{a}_0& =u_1\\ a_1& =\frac{u_2-u_1}{l^e}\end{array}\end{array}$$

即有

$$\begin{array}{c}\begin{array}{rl}u(x)& =u_1+\frac{u_2-u_1}{l^e}x\\ & =(1-\frac{x}{l^e})u_1+\frac{x}{l^e}{u}_2\\ & =\mathbf{N}(x)\cdot {\mathbf{q}}^e\end{array}\end{array}$$

式中：$\mathbf{N}(x)$为形状函数矩阵(shape function matrix)，为

$$\begin{array}{c}\begin{array}{r}\mathbf{N}(x)=\left[\begin{array}{cc}1-\frac{x}{l^e}& \frac{x}{l^e}\end{array}\right]\end{array}\end{array}$$

3 单元应变场的表达

由弹性力学的几何方程，1D问题的应变为

$$\begin{array}{c}\begin{array}{rl}\epsilon (x)& =\frac{\mathrm{d}u(x)}{\mathrm{d}x}\\ & =\left[\begin{array}{cc}-\frac{1}{l^e}& \frac{1}{l^e}\end{array}\right]\left[\begin{array}{c}{u}_1\\ u_2\end{array}\right]\\ & =\mathbf{B}(x)\cdot {\mathbf{q}}^e\end{array}\end{array}$$

其中

$$\begin{array}{c}\begin{array}{rl}\mathbf{B}(x)& =\frac{\mathrm{d}\mathbf{N}(x)}{\mathrm{d}x}\\ & =\left[\begin{array}{cc}-\frac{1}{l^e}& \frac{1}{l^e}\end{array}\right]\end{array}\end{array}$$

称为应变矩阵(strain-displacement matrix)。

4 单元应力场的表达

由弹性力学的物理方程，1D问题的应力为

$$\begin{array}{c}\begin{array}{rl}\sigma (x)& =E^e\epsilon (x)\\ & =E^e\cdot \mathbf{B}(x)\cdot {\mathbf{q}}^e\\ & =\mathbf{S}(x)\cdot {\mathbf{q}}^e\end{array}\end{array}$$

其中

$$\begin{array}{c}\begin{array}{rl}\mathbf{S}(x)& =E^e\cdot \mathbf{B}(x)\\ & =\left[\begin{array}{cc}-\frac{E^e}{l^e}& \frac{E^e}{l^e}\end{array}\right]\end{array}\end{array}$$

称为应力矩阵(stress-displacement matrix)。

5 单元势能的表达

基于式(9)和式(10)，计算得到单元势能的表达式为

$$\begin{array}{c}\begin{array}{rl}\prod ^e& =U^e-W^e\\ & =\frac{1}{2}{\int }_{{\Omega }^e}\sigma (x)\cdot \epsilon (x)\mathrm{d}\Omega -(P_1^e\cdot u_1+P_2^e\cdot u_2)\\ & =\frac{1}{2}{\int }_0^{l^e}\mathbf{S}(x)\cdot {\mathbf{q}}^e\cdot \mathbf{B}(x)\cdot {\mathbf{q}}^e\cdot A^e\mathrm{d}x-(P_1^e\cdot u_1+P_2^e\cdot u_2)\\ & =\frac{1}{2}{\int }_{{\Omega }^e}{\mathbf{q}}^{eT}\cdot {\mathbf{S}}^T(x)\cdot \mathbf{B}(x)\cdot {\mathbf{q}}^e\cdot A^e\mathrm{d}x-(P_1^e\cdot u_1+P_2^e\cdot u_2)\\ & =\frac{1}{2}\left[\begin{array}{cc}{u}_1& u_2\end{array}\right]\left[\begin{array}{c}-\frac{E^e}{l^e}\\ \frac{E^e}{l^e}\end{array}\right]\left[\begin{array}{cc}-\frac{1}{l^e}& \frac{1}{l^e}\end{array}\right]\left[\begin{array}{c}{u}_1\\ u_2\end{array}\right]\cdot A^e\cdot l^e-\left[\begin{array}{cc}{P}_1^e& P_2^e\end{array}\right]\left[\begin{array}{c}{u}_1\\ u_2\end{array}\right]\\ & =\frac{1}{2}\left[\begin{array}{cc}{u}_1& u_2\end{array}\right]\left[\begin{array}{cc}\frac{E^e{A}^e}{l^e}& -\frac{E^e{A}^e}{l^e}\\ -\frac{E^e{A}^e}{l^e}& \frac{E^e{A}^e}{l^e}\end{array}\right]\left[\begin{array}{c}{u}_1\\ u_2\end{array}\right]-\left[\begin{array}{cc}{P}_1^e& P_2^e\end{array}\right]\left[\begin{array}{c}{u}_1\\ u_2\end{array}\right]\\ & =\frac{1}{2}{\mathbf{q}}^{eT}\cdot {\mathbf{K}}^e\cdot {\mathbf{q}}^e-{\mathbf{P}}^{eT}\cdot {\mathbf{q}}^e\end{array}\end{array}$$

上式中第3行到第4行的推导并不是运用的矩阵转置公式$ (M \cdot N)^T = N^T \cdot M^T $，而是因为$ \mathbf{S}(x) $和$ \mathbf{q}(x)^e $均为单位向量，可以直接相互转置和调换顺序

式中：${\mathbf{K}}^e$为单元刚度矩阵(stiffness matrix of element)，即

$$\begin{array}{c}\begin{array}{rl}{\mathbf{K}}^e& =\left[\begin{array}{cc}\frac{E^e{A}^e}{l^e}& -\frac{E^e{A}^e}{l^e}\\ -\frac{E^e{A}^e}{l^e}& \frac{E^e{A}^e}{l^e}\end{array}\right]\\ & =\frac{E^e{A}^e}{l^e}\left[\begin{array}{cc}1& -1\\ -1& 1\end{array}\right]\end{array}\end{array}$$

${\mathbf{P}}^e$为节点力列阵(nodal force vector)，即

$$\begin{array}{c}\begin{array}{rl}{\mathbf{P}}^e& =\left[\begin{array}{c}{P}_1^e\\ P_2^e\end{array}\right]\end{array}\end{array}$$

6 单元刚度方程

对式(13)取极小值，可以得到单元刚度方程(stiffness equation of element)

$$\begin{array}{c}\begin{array}{r}\underset{(2\times 2)}{{\mathbf{K}}^e}\cdot \underset{(2\times 1)}{{\mathbf{q}}^e}=\underset{(2\times 1)}{{\mathbf{P}}^e}\end{array}\end{array}$$

7 变截面杆单元的推导

图2 变截面杆件及单元

上图为一受轴载荷的线性变截面杆件，两端的截面积为$A_1$和$A_2$，长度为$l$，材料的弹性模量为$E$

建立的变截面杆单元如图2(b)所示，其中横截面积为

$$\begin{array}{c}\begin{array}{r}A(x)=(1-\frac{x}{l})A_1+\frac{x}{l}{A}_2\end{array}\end{array}$$

其节点位移列阵为${\mathbf{q}}^e={\left[\begin{array}{cc}{u}_1& u_2\end{array}\right]}^T$，节点力列阵为${\mathbf{P}}^e={\left[\begin{array}{cc}{P}_1& P_2\end{array}\right]}^T$

由于为两个节点，故单元的位移模式取为

$$\begin{array}{c}\begin{array}{r}u(x)=a_0+a_{1}x=\mathbf{N}(x)\cdot {\mathbf{q}}^e\end{array}\end{array}$$

得到单元的形状函数矩阵为

$$\begin{array}{c}\begin{array}{r}\mathbf{N}(x)=\left[\begin{array}{cc}(1-\frac{x}{l})& \frac{x}{l}\end{array}\right]\end{array}\end{array}$$

由几何方程，得到该单元的应变矩阵为

$$\begin{array}{c}\begin{array}{r}\mathbf{B}(x)=\left[\begin{array}{cc}-\frac{1}{l}& \frac{1}{l}\end{array}\right]\end{array}\end{array}$$

该单元的应力矩阵为

$$\begin{array}{c}\begin{array}{r}\mathbf{S}(x)=E\cdot \mathbf{B}(x)=\left[\begin{array}{cc}-\frac{E}{l}& \frac{E}{l}\end{array}\right]\end{array}\end{array}$$

该单元刚度方程为

$$\begin{array}{c}\begin{array}{r}\underset{(2\times 2)}{{\mathbf{K}}^e}\cdot \underset{(2\times 1)}{{\mathbf{q}}^e}=\underset{(2\times 1)}{{\mathbf{P}}^e}\end{array}\end{array}$$

其中单元刚度矩阵为

$$\begin{array}{c}\begin{array}{rl}\underset{(2\times 2)}{{\mathbf{K}}^e}& ={\int }_{\Omega }{\mathbf{S}}^T(x)\cdot \mathbf{B}(x)\cdot \mathrm{d}\Omega \\ & ={\int }_0^l\left[\begin{array}{c}-\frac{E}{l}\\ \frac{E}{l}\end{array}\right]\left[\begin{array}{cc}-\frac{1}{l}& \frac{1}{l}\end{array}\right]\cdot A(x)\mathrm{d}x\\ & ={\int }_0^l\frac{E}{l^2}\left[\begin{array}{cc}1& -1\\ -1& 1\end{array}\right]\cdot [(1-\frac{x}{l})A_1+\frac{x}{l}{A}_2]\mathrm{d}x\\ & =\frac{E}{l^2}\left[\begin{array}{cc}1& -1\\ -1& 1\end{array}\right]{\int }_0^l[(1-\frac{x}{l})A_1+\frac{x}{l}{A}_2]\mathrm{d}x\\ & =\frac{E}{l^2}\left[\begin{array}{cc}1& -1\\ -1& 1\end{array}\right][(x-\frac{x^2}{2l})A_1+\frac{x^2}{2l}{A}_2]{|}_0^l\\ & =\frac{E}{l^2}\left[\begin{array}{cc}1& -1\\ -1& 1\end{array}\right][(l-\frac{l}{2})A_1+\frac{l}{2}{A}_2]\\ & =\frac{E(A_1+A_2)}{2l}\left[\begin{array}{cc}1& -1\\ -1& 1\end{array}\right]\end{array}\end{array}$$

可以看到，与等截面杆单元相比(式(14))，该线性变截面单元的截面取了平均面积。

8 多连杆的整体刚度矩阵

图3 局部坐标系中一维三连杆结构

上图为局部坐标系中的一维三连杆结构，故每个杆件的单元刚度方程分别为

$$\begin{array}{c}\begin{array}{rl}{\mathbf{K}}^{(1)}\cdot {\mathbf{q}}^{(1)}& ={\mathbf{P}}^{(1)}\\ \left[\begin{array}{cc}{K}_{11}^{(1)}& K_{12}^{(1)}\\ K_{21}^{(1)}& K_{22}^{(1)}\end{array}\right]\cdot \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}\end{array}$$

$$\begin{array}{c}\begin{array}{rl}{\mathbf{K}}^{(2)}\cdot {\mathbf{q}}^{(2)}& ={\mathbf{P}}^{(2)}\\ \left[\begin{array}{cc}{K}_{22}^{(2)}& K_{23}^{(2)}\\ K_{32}^{(2)}& K_{33}^{(2)}\end{array}\right]\cdot \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}\end{array}$$

$$\begin{array}{c}\begin{array}{rl}{\mathbf{K}}^{(3)}\cdot {\mathbf{q}}^{(3)}& ={\mathbf{P}}^{(3)}\\ \left[\begin{array}{cc}{K}_{33}^{(3)}& K_{34}^{(3)}\\ K_{43}^{(3)}& K_{44}^{(3)}\end{array}\right]\cdot \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}\end{array}$$

组装成整体刚度方程为

$$\begin{array}{c}\begin{array}{rl}{\mathbf{K}}^{All}\cdot {\mathbf{q}}^{All}& ={\mathbf{P}}^{All}\\ \left[\begin{array}{cccc}{K}_{11}^{(1)}& K_{12}^{(1)}& 0& 0\\ K_{21}^{(1)}& K_{22}^{(1)}+K_{22}^{(2)}& K_{23}^{(2)}& 0\\ 0& K_{32}^{(2)}& K_{33}^{(2)}+K_{33}^{(3)}& K_{34}^{(3)}\\ 0& 0& K_{43}^{(3)}& K_{44}^{(3)}\end{array}\right]\cdot \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]\\ \left[\begin{array}{cccc}{K}_{11}^{(1)}& K_{12}^{(1)}& 0& 0\\ K_{21}^{(1)}& K_{22}^{(1)}+K_{22}^{(2)}& K_{23}^{(2)}& 0\\ 0& K_{32}^{(2)}& K_{33}^{(2)}+K_{33}^{(3)}& K_{34}^{(3)}\\ 0& 0& K_{43}^{(3)}& K_{44}^{(3)}\end{array}\right]\cdot \left[\begin{array}{c}{u}_1\\ u_2\\ u_3\\ u_4\end{array}\right]& =\left[\begin{array}{c}{P}_1\\ P_2\\ P_3\\ P_4\end{array}\right]\end{array}\end{array}$$

整体刚度矩阵的组合方式如下图所示

图4 局部坐标系中一维多连杆结构整体刚度矩阵组装方式

参考

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