有限元分析及应用学习笔记
到底什么是有限元，有限元到底有什么用
有限元分析讲的其实上不全是有限元分析，前面用了很大的篇幅讲述了材料力学的内容，而有限元只是处理相关问题的一个方法，即化大为小，将复杂的东西进行标准的、参数化的几何划分，来进行单元研究，最后进行集合组装得到我们整体物体的性质分析。
有限元分析通过针对复杂几何域进行分片，然后对于对象的材料，我们一般简化材料的行为，对于变形体材料有以下基本假设：
- 连续性：无空隙，可以采用连续函数描述
- 均匀性：物质的特征不随位置变化
- 各向同性：物质特征与方向性无关
- 线性弹性：变形与外力作用的关系是线性的，使得描述材料性质的方程是线性的
- 小变形：变形远小于物体的几何尺寸，可以忽略二阶及以上的高阶小量
基本名词释义
常见研究单元
- 杆：杆只能承受轴向的拉压力，只有一个方向的受力以及变形；
- 平面纯弯梁：只有弯曲形变,不受轴向伸缩形变的梁，一般来说会考虑弯曲的转角以及挠度；
- 一般梁：包含杆单元和平面纯弯梁的特性，会考虑转角、挠度和轴向变化
- 平面三角形、矩形等平面单元及四面体、正六面体单元等空间单元
常见物理量及符号
- 位移（挠度） $u$
- 应力：$\sigma_{xx}=\frac{F}{A}$ 表示 x 方向的正应力（$A$表示横截面积）；$\tau_{xy}$表示力的作用面的法向量方向是 y 坐标轴方向，力的方向为 x 坐标轴方向的剪应力（剪力）
- 应变：正应变$\epsilon_{xx}=\frac{u}{L}$，剪应变$\gamma_{xy}=\alpha+\beta=\frac{\partial{u}}{\partial{y}}+\frac{\partial{v}}{\partial{x}}$
- 弹性模量：$E=\frac{\sigma}{\epsilon}$
- 载荷：节点承载的负荷（力）
- 支反力：在被固定的点，用于维持点被固定住所施加的外力
- 刚度$K$，满足$Kq=F$（$q$为每个节点位移的列向量，$F$为每个节点力的列向量）
- 泊松比$\mu$：材料在单向受拉或受压时，横向正应变$\epsilon_{y}$与轴向正应变$\epsilon_x$的绝对值的比值的相反数$\mu=-\frac{\epsilon_y}{\epsilon_x}$
- 均布力$\overline{P}$：在外表面上均匀分布的外力
- 算子$L$：拉普拉斯算子，对位移求二阶导
- 残差函数$\Re$、形状函数矩阵$N$、几何矩阵$B$、应力矩阵$S$
三大类变量三大类方程
- 三大类变量：位移、应力、应变。其中应力与应变可以由位移表示。
- 三大类方程：平衡方程、几何方程、物理方程
- 其实我们在用的时候，也会用到两类边界条件：位移边界、力边界
处理数据的方法
微分方程求解的方法
- 解析法：对于简单的微分方程常用，但是对于比较复杂的方程难以适用。
- 差分法：对整体进行分段，用差商来代替微分，对每段列差分方程，得到线性方程组，就可以得到解在每段的近似值
- 试函数方法：设置一个带着参数，且满足边界条件的试函数（$R(x)$），通过让残差函数在试函数的限定区域内的累计最小（$\int_{area} \phi(x)R(x)dx=0$），解得试函数的参数。
- 试函数的精度取决于试函数设置的形式，如果这个形式包括解析解，则可以求出解析解；但是如果不包括，只能在试函数设置形式这个范围内，获得一个最优解。
- 试函数不一定非要是在全域上讨论，可以对全域进行分段，每一段用不同的试函数。但是这样的坏处在于这样的一系列试函数在高阶导数上不连续（在后面我们要至少保证二阶导数连续）
坐标变换的方法
$T^e$是坐标变换矩阵，以平面的杆单元（坐标位移$\begin{bmatrix} \overline{u}_1 & \overline{v}_1 & \overline{u}_2 & \overline{v}_2\end{bmatrix}$）为例，$T^e=$$$\begin{bmatrix} \cos{\alpha} & \sin{\alpha} & 0 & 0 \\ 0 & 0 & \cos{\alpha} & \sin{\alpha}\end{bmatrix}$$，我们对位移进行坐标变换$q^e=T^e\cdot q^e$，并带入原刚度方程$K^eq^e=F^e$，并在等式两边同时左乘$T^{Et}$,得到方程$T^{eT}K^e(T^eq^e)=T^{eT}F^e$，并简记为$\overline{K}^e\cdot \overline{q}^e = \overline{F}^e$。此时的新刚度矩阵的大小对应着新坐标系位移的个数，即为$4\times4$。
指标记法
- 自由指标：在每一项（加减号中间的那部分）中只出现一次的下标叫做自由指标，可以自由变化。
- 哑指标：意味着求和（爱因斯坦求和约定）
$$a_{ij}x_j=b_i \Leftrightarrow \sum_{j=1}^3 a_{ij}x_j=b_i \Leftrightarrow \begin{cases}{}
a_{11}x_1 + a_{12}x_2 + a_{13}x_3 = b_1 \\
a_{21}x_1 + a_{22}x_2 + a_{23}x_3 = b_2 \\
a_{31}x_1 + a_{32}x_2 + a_{33}x_3 = b_3
\end{cases} $$
Voigt 标记
将高阶自有指标的张量写成低阶张量形式的过程叫做 Voigt 标记，规则叫做 Voigt 规则（这个规则好神奇百度谷歌都查不到）。
用途：优化存储，减少对称矩阵的对称部分，降低矩阵大小
方法简述：从$a_{11}$开始，先沿着矩阵对角线向右下方移动，走到$a_{nn}$后，按逆时针顺序向上写（再向上，走到头再向左，再右下循环）。
方法举例：$\sigma_{ij}=\begin{bmatrix} \sigma_{11} & \sigma_{12} & \sigma_{13} \\ \sigma_{21} & \sigma_{22} & \sigma_{23} \\ \sigma_{31} & \sigma_{32} & \sigma_{33} \end{bmatrix}$，可以得到一维 Voigt 矩阵$\sigma_p=\begin{bmatrix}\sigma_{11} \\ \sigma_{22} \\ \sigma_{33} \\ \sigma_{23} \\ \sigma_{13} \\ \sigma_{12} \end{bmatrix}$
试函数法求解变形体力学方程
讲到的试函数法有加权残值法、虚功原理法、最小势能原理法、变分方法。
我们要设置满足边界条件的解，带入原控制方程，处理误差求得待定系数。
试函数法的要点
将一组满足所有边界条件试函数（基底函数）$\phi_i(x)$的线性组合$$\hat{v}(x)=c_1\phi_1(x)+c_2\phi_2(x)+\cdots+c_n\phi_n(x)$$（$c_i$为待定系数）代回边界条件$\begin{cases}{}BC(u):g_u(\hat{v}(x)=0\quad on\ S_u \\ BC(p):g_p(\hat{v}(x))=0\quad on\ S_p\end{cases}$，并通过指标处理残值函数$\Re(x)=L(\hat{v}(x))+\overline{b}\neq0$。
Galerkin 加权残值法
Galerkin 加权残值法将试函数取的基底函数作为权重，对残值函数$\Re(x,y,z)$进行加权（几个参数用几个基底函数，每次加权用一个，然后把 n 个积分=0 联立起来），在整个域中积分。
$$\int_\Omega w_{t_i}\Re(\mathrm{x,y,z})\mathrm{d}\Omega=0,\quad w_{t_i}=\phi_i$$
残值最小二乘法
残值最小二乘法调整待定系数使得残值平方和积分最小（一般此时积分许多参数都与积分微元无关，所以类似于幂函数的积分）。我们求得$E_{rr}$后，分别对 n 个参数求偏导$\frac{\partial{E_{rr}}}{\partial{c_i}}=0,\quad i=1,2\cdots N$来求得参数。
$$\underset{c_1,c_2,\cdots,c_n}{min}\left[E_n=\int_\Omega w_t\Re^2(\mathrm{x,y,z})\mathrm{d}\Omega\right], w_t=1$$
配点法
相当于简单强迫余量在域内的 N 个点上为 0（一般 n 个参数取 n 个点），准不准的就看造化了。
$$\int\Re(x_i,y_i,z_i)\mathrm{d}\Omega=0\quad i=1,2\cdots N$$
降低试函数高阶导数及边界条件要求
将应力应变均用位移表达后，可以得到基于位移的控制方程（平衡方程）
$$\begin{cases}{}
\frac{E}{1-\mu^2}\left(\frac{\partial^2{u}}{\partial{x^2}}+\frac{1-\mu}{2}\frac{\partial^2{u}}{\partial{y^2}}+\frac{1+\mu}{2}\frac{\partial^2{v}}{\partial{x}\partial{v}}\right)+\overline{b}_x=0\\
\frac{E}{1-\mu^2}\left(\frac{\partial^2{v}}{\partial{y^2}}+\frac{1-\mu}{2}\frac{\partial^2{v}}{\partial{x^2}}+\frac{1+\mu}{2}\frac{\partial^2{u}}{\partial{x}\partial{v}}\right)+\overline{b}_y=0
\end{cases}$$
以及基于位移的力边界
$$\begin{cases}{}
\frac{E}{1-\mu^2}\left[n_x(\frac{\partial{u}}{\partial{x}}+\mu\frac{\partial{v}}{\partial{y}}+n_y(\frac{1-\mu}{2})(\frac{\partial{u}}{\partial{y}}+\frac{\partial{v}}{\partial{x}})\right]=\overline{p}_x\\
\frac{E}{1-\mu^2}\left[n_y(\frac{\partial{v}}{\partial{y}}+\mu\frac{\partial{u}}{\partial{x}}+n_x(\frac{1-\mu}{2})(\frac{\partial{v}}{\partial{x}}+\frac{\partial{vu}}{\partial{y}})\right]=\overline{p}_y
\end{cases}$$
知对导数要求可以降低到一半。
虚功原理与最小势能原理
虚功原理：$\delta U=\delta W$，即$\int_\Omega \sigma_{ij}\cdot\delta\epsilon_{ij}\cdot \mathrm{d}\Omega=\int_\Omega\overline{b}_i\delta u_i\mathrm{d}\Omega+\int_{S_p}\overline{p}_i\delta u_i\mathrm{d}A$，其中虚位移$\delta c_i$具有任意性，不恒为零，与其相乘为 0 则另一个恒为零。
最小势能原理：总势能$\prod(\hat{v}(x))=U-W=\frac{1}{2}\int_{\Omega}\sigma_x\epsilon_x\mathrm{d}\Omega-\left(\int_\Omega\overline{b}_i\delta u_i\mathrm{d}\Omega+\int_{S_p}\overline{p}_i\delta u_i\mathrm{d}A\right)$，然后分别对参数求导=0，可以解得参数$c_i$。
虚功原理与最小势能原理等价，都需要全域的试函数$BC(u)$和低阶导数要求。最小势能原理与加权残值法在经过分部积分处理后是等价的。
有限元方法
- 几何域的离散化：将复杂的几何域切分为参数化、标准化、规范化的单元，方便对这些单元一一研究。
- 单元研究：通过每个单元的能量原理、虚功原理、最小势能原理等研究单元，得到$K^eq^e=P^e$
- 单元的集成装配：部分坐标系与整体坐标系的转换，尤其是旋转的变换。
- 处理位移边界条件：讲力$P$拆为已知的外力$\overline{P}_F$和未知的支反力$P_R$，将位移拆为未知的节点位移$q_u$和已知的节点位移$\overline{q}_k$，可以得到矩阵等式$\begin{bmatrix}K_1 & K_2 \\K_3 & K_4\end{bmatrix}\begin{bmatrix}q_u\\ \overline{q}_k \end{bmatrix}=\begin{bmatrix}\overline{P}_F \\ P_R\end{bmatrix}$
- 计算支反力等物理量：通过上一步，的方程$\begin{cases}K_1q_u+K_2\overline{q}_k=\overline{P}_F\\ K_3q_u+K_4\overline{q}_k=P_R\end{cases}$，求未知节点位移$q_u=K_1^{-1}(\overline{P}_F-K_2\overline{q}_k)$，与支反力$P_R=K_3q_u+K_4\overline{q}_k=K_3K_1^{-1}(\overline{P}_F-K_2\overline{q}_k)+K_4\overline{q}_k$，通过几何方程物理方程可以求得其他物理量。
重要的定理公式
- 剪应力互等定理$\tau_{xy}=\tau_{yx}$，由力矩平衡导出
- 平面问题平衡方程 $\sigma_{ij,j}+\overline{b}_i=0 \Leftrightarrow \begin{cases}{} \frac{\partial{\sigma_{xx}}}{\partial{x}}+ \frac{\partial{\sigma_{xy}}}{\partial{y}}+\overline{b}_x=0 \\ \frac{\partial{\sigma_{yy}}}{\partial{y}}+ \frac{\partial{\sigma_{xy}}}{\partial{x}}+\overline{b}_y=0 \end{cases}$，其中指标形式的逗号表示对逗号后面的指标值求偏导
- 平面问题几何方程：相对伸长量（应变）$\epsilon_{xx}=\frac{\partial{u}}{\partial{x}}$、$\epsilon_{yy}=\frac{\partial{v}}{\partial{y}}$，剪应变$\gamma_{xy}=\frac{\partial{u}}{\partial{y}}+\frac{\partial{v}}{\partial{x}}$。
- 使用中常记为$\epsilon_{ij}=\frac{1}{2}\left(u_{i,j}+u_{j,i}\right)$，逗号表示同上，而且在$i\neq j$时，$\epsilon_{ij} = \frac{1}{2}\gamma_{ij}$（据弹幕说这是爱因斯坦的记法，这个奇怪的记法就这么流传下来了）
- 位移分量$u,v$可以唯一确定三个应变分量$\epsilon_{xx},\epsilon_{yy},\gamma_{xy}$，但是应变不能唯一确定位移分量，需要满足变形协调条件$\frac{\partial^2{\epsilon_{xx}}}{\partial{y^2}}+\frac{\partial^2{\epsilon_{yy}}}{\partial{x^2}}=\frac{\partial^2{\gamma_{xy}}}{\partial{x}\partial{y}}$。其代表着材料变形的时候不撕裂、不重叠。
- 平面问题物理方程：根据泊松效应，一个方向的正应变不仅与该方向的应力有关，还受到其他方向应力的影响。
- 方程表述$\begin{cases}{}
\epsilon_{xx}=\frac{1}{E}\left(\sigma_{xx}-\mu\sigma_{yy}\right) \\
\epsilon_{yy}=\frac{1}{E}\left(\sigma_{yy}-\mu\sigma_{xx}\right) \\
\gamma_{xy}=\frac{1}{G}\tau_{xy}
\end{cases}$，其中$G = \frac{E}{2\left(1+\mu\right)}$。 - 也常用其逆形式$\begin{cases}{}
\sigma_{xx}=\frac{E}{1-\mu^2}\left(\epsilon_{xx}+\mu\epsilon_{yy}
\right) \\
\sigma_{yy}=\frac{E}{1-\mu^2}\left(\epsilon_{yy}+\mu\epsilon_{xx}
\right) \\ \tau_{xy}=G\gamma_{xy}
\end{cases}$。 - 按照指标规则和 Voigt 规则，可以将指标降阶，得到物理方程的矩阵形式$\begin{bmatrix}\sigma_{xx} \\ \sigma_{yy} \\ \tau_{xy}\end{bmatrix} = \begin{bmatrix}
\frac{E}{1-\mu^2} & \frac{E\mu}{1-\mu^2} & 0 \\
\frac{E\mu}{1-\mu^2} & \frac{E}{1-\mu^2} & 0 \\
0 & 0 & G \
\end{bmatrix} \cdot \begin{bmatrix} \epsilon_{xx} \\ \epsilon_{yy} \\ \gamma_{xy}\end{bmatrix}$
- 方程表述$\begin{cases}{}
- 位移边界$S_u$：$u_i = \overline{u}_i \quad on\ S_u$（$\overline{u}_i$表示给定的位移约束）
- 力边界$S_p$：$\sigma_{ij}n_{j}=\overline{p}_i \quad on\ S_p$（$n_j$表示边界一点上外发现的方向余弦，$\overline{p}_i$表示给定的分布力）
- 刚度方程：$K^eq^e=F^e$
- 分布力的处理：
- 长度$L$梁单元
- 中点受到力$P$：$R_A=R_B=-\frac{P}{2}$，$M_A=-\frac{PL}{8}$；
- 在距左长度$a$处（距右长度$b=L-a$处）受到力$P$：$R_A=-\frac{Pb^2}{L}\cdot (3a+b)$，$R_B=-\frac{Pa^2}{L}\cdot (a+3b)$，$M_A=-\frac{Pab^2}{L^2}$，$M_B=-\frac{Pa^2b}{L^2}$；
- 均布载荷$P_0$：$R_A=R_B=-\frac{p_0L}{2}$，$M_A=-\frac{p_0L^2}{12}$，$M_B=\frac{p_0L^2}{12}$；
- 受到直角三角形（$P_A$最小，$P_B$最大）载荷$P_0$：$R_A=-\frac{3p_0L}{20}$，$R_B=-\frac{7p_0L}{20}$，$M_A=-\frac{p_0L^2}{30}$，$M_B=\frac{p_0L^2}{20}$；
- 在距左长度为 a 的部分受到均布载荷$P_0$：$R_A=-\frac{p_0a}{2L^3}\cdot(a^3-2a^2L+2L^3)$，$R_B=-\frac{p_0a^3}{2L^3}\cdot (2L-a)$，$M_A=-\frac{p_0a^2}{12L^2}\cdot (3a^2-8aL+6L^2)$，$M_B=\frac{p_0a^3}{12L^2}\cdot (4L-3a)$；
- 受到等腰三角形均布载荷$P_0$：$R_A=R_B=-\frac{p_0L}{4}$，$M_A=-\frac{5p_0L^2}{96}$，$M_B=\frac{5p_0L^2}{96}$；
- 在距左长度$a$处（距右长度$b=L-a$处）受到弯矩$M_0$：$R_A=-\frac{6M_0ab}{L^3}$，$R_B=\frac{6M_0ab}{L^3}$，$M_A=-\frac{M_0b}{L^2}(3a-L)$，$M_B=-\frac{M_0a}{L^2}(3b-L)$。
- 平面 3 节点三角形单元
- 单元自重：三个节点竖直方向平摊，水平方向无荷载；
- 三角形一条边均布侧压（垂直于边）：在 12 节点垂直收到均布载荷$p_0$，1、2 号节点受水平方向载荷$\frac{1}{2}p_0t(y_1-y_2)$，竖直方向载荷$\frac{1}{2}p_0t(x_2-x_1)$；
- 三角形一条边 x 方向均布侧压：边的两端点 x 方向均摊；
- 三角形一条边 x 方向三角形侧压：受力大的那边占$\frac{2}{3}$，受力小的那边$\frac{1}{4}$。
- 长度$L$梁单元
一些经常讨论的情况
平面应力问题
条件
一块很薄的等厚度板（厚度为 z 轴方向）
结果
$\sigma_{zz}=0, \tau_{xz}=0, \tau_{yz}=0, \gamma_{xz}=0, \gamma_{yz}=0$
由物理方程知$\epsilon_{zz}=-\frac{\mu}{E}\left(\sigma_{xx}+\sigma_{yy}\right)\neq0$，但$\sigma_{zz}=0$
平面应变问题
条件
一个无限长的等截面柱体（无限长的是 z 轴方向）
结果
- z 轴方向的位移$w(x,y,z)$应变$\epsilon_{zz}$均为 0，与 z 相关的应变$\gamma_{xz},\gamma_{yz}$也为 0
- 同样 z 方向的应力$\sigma_{zz}$，与 z 相关的应力$\tau_{yz}, \tau_{xz}$也为 0
- 有物理方程知$\sigma_{zz}=\mu\left(\sigma_{xx}+\sigma_{yy}\right) \neq 0$，但$\epsilon_{zz}=0$（与平面应力相反）
比较
| 种类 | 平面应力问题 | 平面应变问题 | |
|---|---|---|---|
| 应力 | $\sigma_{xx}(x,y,z)$ | ||
| $\sigma_{yy}(x,y,z)$ | |||
| $\sigma_{zz}(x,y,z)$ | $0$ | $\mu(\sigma_{xx}+\sigma_{yy})$ | |
| $\tau_{xy}(x,y,z)$ | |||
| $\tau_{xz}(x,y,z)$ | $0$ | ||
| $\tau_{yz}(x,y,z)$ | $0$ | ||
| 应变 | $\epsilon_{xx}(x,y,z)$ | $\frac{1}{E}\left(\sigma_{xx}-\mu\sigma_{yy}\right)$ | $\frac{1-\mu^2}{E}\left[\sigma_{xx}-\left(\frac{\mu}{1-\mu}\right)\sigma_{yy}\right]$ |
| $\epsilon_{yy}(x,y,z)$ | $\frac{1}{E}\left(\sigma_{yy}-\mu\sigma_{xx}\right)$ | $\frac{1-\mu^2}{E}\left[\sigma_{yy}-\left(\frac{\mu}{1-\mu}\right)\sigma_{xx}\right]$ | |
| $\epsilon_{zz}(x,y,z)$ | $-\frac{\mu}{E}\left(\sigma_{xx}+\sigma_{yy}\right)$ | $0$ | |
| $\gamma_{xy}(x,y,z)$ | $\frac{1}{G}\tau_{xy}$ | ||
| $\gamma_{xz}(x,y,z)$ | $0$ | ||
| $\gamma_{yz}(x,y,z)$ | $0$ | ||
\begin{table}[] |
平面刚体位移
条件
物体内不会产生任何应变，即$\epsilon_{xx}=\frac{\partial{u}}{\partial{x}}=0$，$\epsilon_{yy}=\frac{\partial{v}}{\partial{y}}=0$，$\gamma_{xy}=\frac{\partial{v}}{\partial{x}}+\frac{\partial{u}}{\partial{y}}=0$。
结果
$\begin{cases}{}u(x,y)=-\omega_0y+u_0\\v(x,y)=\omega_0x+v_0 \end{cases}$，其中$u_0,v_0,\omega_0$表征平面刚体位移常数，$u_0,v_0$为整个物体在 x,y 方向的刚体平移量，$\omega_0$为整个物体缸体转动的角度。
简单拉杆问题
变量与方程
- 变量：位移$u(x)$，应力$\sigma_x(x)$，应变$\epsilon_x(x)$；
- 方程：平衡方程$\frac{\mathrm{d}\sigma_x}{\mathrm{d}x}=0$；几何方程$\epsilon_x=\frac{\mathrm{d}u}{\mathrm{d}x}$；物理方程$\epsilon_x=\frac{\sigma_x}{E}$;
- 边界方程：位移方程$BC(u)$：$u(x)|_{x=0}=0$；力边界$BC(p)$：$\sigma_x(x)|_{x=l}=\frac{P}{A}=\overline{p}_x$。
求解结果
应力$\sigma_x(x)=\frac{P}{A}$；应变$\epsilon_x(x)=\frac{P}{EA}$；位移$u(x)=\frac{P}{EA}x$
常见单元的自由度
- 杆单元：$u_1,u_2$两个
- 平面纯弯梁：$v_1,\theta_1,v_2,\theta_2$四个
- 一般平面梁（杆+平面纯弯梁）：$u_1,v_1,\theta_1,u_2,v_2,\theta_2$六个
- 空间梁单元（三个方向弯曲）：梁的扭转$\theta_{x1},\theta_{x2}$，沿轴向拉伸$u_1,u_2$，xOy 平面弯曲$v_1,\theta_{z1},v_2,\theta_{z2}$，xOz 平面弯曲$w_1,\theta_{y1},w_2,\theta_{y2}$共 12 个
- 平面 3 节点三角形：三个节点每个方向分别的位移$u_1,v_1,u_2,v_2,u_3,v_3$6 个
- 平面 4 节点矩形单元：三个节点每个方向分别的位移$u_1,v_1,u_2,v_2,u_3,v_3,u_4,v_4$8 个
- 轴对称体元：$u_r,w$两个
- 3 节点三角形环形单元：$u_{r1},w_1,u_{r2},w_2,u_{r3},w_3$6 个
- 4 节点矩形环形单元：$u_{r1},w_1,u_{r2},w_2,u_{r3},w_3,u_{r4},w_4$8 个
放在最后
感谢彭昊学长对本人学习有限元分析的指导与帮助，希望大家也可以更高效的学会有限元分析。
另推荐大家看几篇总结：
