使用有限单元法对简单结构进行分析的方法
最近，我的参赛作品涉及对简单结构进行应力分析。虽然有限元分析是一种有效的应力分析工具，但市面上大多数有限元分析软件都是闭源的，且价格不菲。如果自己动手实现，就需要一定的工程力学知识。因此，我尝试用计算机专业同学易于理解的方式，介绍一个有限元分析的数学模型。本人是一名非土木等工科专业的学生，如果有错误，还请专业人士不吝指正。
一.平面弹性力学中的传统应力分析方法
该部分主要介绍应力分析的基础知识，如果理解有难度的话可跳过。在弹性力学中，主要有三大基本方程，分别是平衡微分方程、几何方程和物理方程。这三大方程构成了平面弹性力学的基础，通常被称为“弹性力学的基本方程”，它们共同描述了材料在受力情况下的行为。
1.常见的物理量
要理解三大基本方程，需要先认识几个物理量
(1) 体力
- 定义：体力是作用于单位体积上的力，通常指重力、惯性力等作用于整个体积的外力。
- 符号：常用 f 表示，单位为 N/m^2
(2) 面力
- 定义：面力是作用于单位面积上的力，通常指外部施加在物体表面上的力，如压力、摩擦力等。
- 符号：常用 \bar{f} 表示，单位为 $N/m^2$（即帕斯卡）。
(3) 内力
- 定义：内力是物体内部各部分之间的相互作用力，通常是由于材料的应力和应变引起的。内力维持物体的结构稳定。
- 符号：常用 N 表示，具体形式可以用应力张量 σ 或者简单的 F 来表示，单位为 N。
(4) 应力
- 定义：应力是单位面积上承受的内力，描述材料内部的力分布。它反映了材料在外力作用下的反应。
- 符号：常用 \sigma , \tau 表示，分别作为正应力和切应力 ,单位为 $N/m^2$（即帕斯卡）。
(5) 应变
-
定义：应变是材料在外力作用下的变形程度，表示原始长度的相对变化。它是无量纲的，通常用百分比或小数表示。
-
符号：常用 \varepsilon , \gamma 表示，分别代表正应变和切应变 ,定义为：
\varepsilon =\frac{\varDelta L}{L_0}
其中，ΔL 是长度的变化量，$L_0$是原始长度。
(6) 位移
-
定义：位移是物体从原始位置到新位置的变化量，描述了物体在外力作用下的移动情况。
-
符号：常用 u , v , w 表示，分别作为不同坐标下的位移 ,单位为米（m）。
2.平面问题三大方程
了解这些物理量后，可以给出平面弹性力学的三个方程
(1)平面问题的平衡微分方程：
\left\{ \begin{array}{l}
\frac{\partial \sigma _x}{\partial x}+\frac{\partial \tau _{yx}}{\partial y}+f_x=0\\
\frac{\partial \sigma _y}{\partial y}+\frac{\partial \tau _{xy}}{\partial x}+f_y=0\\
\end{array} \right.
(2)平面问题的几何方程：
\left\{ \begin{array}{l}
\varepsilon _x=\frac{\partial u}{\partial x}\\
\varepsilon _y=\frac{\partial v}{\partial y}\\
\gamma _{xy}=\frac{\partial v}{\partial x}+\frac{\partial u}{\partial y}\\
\end{array} \right.
(3)物理方程：
\left\{ \begin{array}{l}
\varepsilon _x=\frac{1}{E}\left( \sigma _x-\mu \sigma _y-\mu \sigma \right) ,\ \gamma _{yz}=\frac{1}{G}\tau _{yz}\\
\varepsilon _y=\frac{1}{E}\left( \sigma _y-\mu \sigma _z-\mu \sigma _x \right) ,\ \gamma _{yz}=\frac{1}{G}\tau _{zx}\\
\varepsilon _z=\frac{1}{E}\left( \sigma _z-\mu \sigma _x-\mu \sigma _y \right) ,\ \gamma _{xy}=\frac{1}{G}\tau _{xy}\\
\end{array} \right.
其中E为弹性模量，$\mu$为泊松比
3.位移和应力边界条件
我们在学习微积分时接触过微分方程，可以将这里的边界条件理解为用于消除解微分方程时不定积分所带来的常数项的影响。通过带入边界条件，可以使方程具有唯一解。从物理意义上理解，边界条件表示在边界上的位移与约束，或应力与面力之间的关系。
位移边界条件：
设在$S_u$部分边界上给定位移分量$\bar{u}(s)和\bar{v}(s)$，则有
\left( u_s \right) =\bar{u}\left( s \right) ,\left( v_s \right) =\bar{v}\left( s \right)
应力边界条件：
\left\{ \begin{array}{l}
\bar{f}_x\left( s \right) =\left( l\sigma _x+m\tau _{yx} \right) _s\\
\bar{f}_y\left( s \right) =\left( m\sigma _y+l\tau _{xy} \right) _s\\
\end{array} \right.
这里的 l , m 为方向余弦
通过建立平面弹性力学的三大方程，理论上只要给出物体结构的应力、位移约束条件和体力分布，就可以计算出任意一点的应力、应变和位移。然而，由于解微分方程的困难，实际上我们能通过这种方式解决的问题非常有限，而且难以使用计算机程序进行应力分析。而有限元法有效地解决了这一问题，使得复杂结构的分析变得更加可行。
二. 有限单元法
1.有限单元法概述
有限单元法（Finite Element Method, FEM）是一种数值技术，用于求解复杂工程和物理问题中的偏微分方程。它通过将连续的物理区域离散化为有限数量的小部分（即“单元”），并在每个单元内建立简单的数学模型，从而近似整个区域的行为。
有限单元法相比于传统方法具有多个优势。首先，它能够处理复杂几何形状和边界条件，适用于非线性和大变形问题，而传统方法通常只适用于简单结构。其次，有限单元法将整个结构离散化为小单元，使得问题的求解过程可以通过数值方法进行，这极大地简化了求解复杂微分方程的过程。此外，有限单元法能够提供局部应力和应变的详细分布信息，而传统方法往往只能给出全局平均值或近似解。最后，随着计算机技术的发展，有限单元法的计算效率大幅提高，使得在工程实践中广泛应用成为可能。这些优势使得有限单元法成为现代工程分析中不可或缺的工具。
2.有限元分析的一般过程
有限元分析所涉及的内容十分复杂，在这里我们仅分析如图所示由大小和形状完全相同的正方体组成的结构。
(1) 结构离散化
从 z=0 层开始，将结构按坐标分割成正六面体，给每个正六面体及顶点编号，顶点编号顺序如图。定点编号包括局部节点编号及整体节点编号。
(2) 数据准备
输入弹性模量E，泊松比 \mu，材料密度 \rho，位移边界条件（z=0处位移为0），力边界条件（物体底部有约束力）
边界条件在后续计算中会体现
(3) 计算单元刚度矩阵
中学物理告诉我们，弹簧的长度变化量与受力成正比。在有限单元法中对每个单元进行分析的时候，我们可以将单元类比弹簧，得到公式
K^ea=P^e
这里的a为单元八个结点的位移向量组成的列阵，其形状为24*1，P为结点载荷，可以理解为八个结点在x，y，z方向上的受力组成的列阵，其形状为24*1。K为单位刚度矩阵，与弹簧的弹性系数类似，其形状为24*24
a. Gauss积分法
要计算K和P，需要引入Gauss积分法
目标：设定要计算的的积分形式为
I=\int_a^b{f\left( x \right) dx}
变换区间：为了方便起见，通常将积分区间从 [a,b] 转换为标准区间 [−1,1]。这一变换可以通过线性变换实现：
x=\frac{b-a}{2}\xi +\frac{a+b}{2}
其中 \xi \in \left[ -1,1 \right]
变量代换：相应的，dx的变化量为
dx=\frac{b-a}{2}d\xi
因此，积分变为：
I=\int_{-1}^1{f\left( \frac{b-a}{2}\xi +\frac{a+b}{2} \right) \frac{b-a}{2}d\xi}=\int_{-1}^1{g\left( \xi \right) d\xi}=\sum_{i=1}^n{H_ig\left( \xi _i \right)}
此处的 H_i 为积分权系数，\xi_i 为积分点，对于正六面体等参单元，取 \xi_1=0.557350269189626,\xi_2=-0.557350269189626 两个积分点，H_1=H_2=1
通过引入Gauss积分法，我们得以将积分运算转换为矩阵运算
b. 三维单元插值函数矩阵
如图所示的六面体8结点一次单元的插值函数为
N_i=\frac{1}{8}\left( 1+\xi _i\xi \right) \left( 1+\eta _i\eta \right) \left( 1+\zeta _i\zeta \right) \left( i=1,2,\cdots ,8 \right)
式中 \left( \xi _i,\eta _i,\zeta _i \right) 为结点i的坐标：
1（-1，-1，-1）；2（+1，-1，-1）；3（+1，+1，-1）；4（-1，+1，-1）
5（-1，+1，-1）；6（+1，-1，+1）；7（+1，+1，+1）；8（-1，+1，+1）
设
I=\left[ \begin{matrix}
1& 0& 0\\
0& 1& 0\\
0& 0& 1\\
\end{matrix} \right]
那么插值函数矩阵
N=\left[ \begin{matrix}
IN_1& IN_2& IN_3\\
\end{matrix}\ \cdots \ IN_8 \right]
c. 应变矩阵
三维等参单元的应变子矩阵为
B_i=\left[ \begin{matrix}{}
\frac{\partial N_i}{\partial x}& 0& 0\\
0& \frac{\partial N_i}{\partial y}& 0\\
0& 0& \frac{\partial N_i}{\partial z}\\
\frac{\partial N_i}{\partial y}& \frac{\partial N_i}{\partial x}& 0\\
0& \frac{\partial N_i}{\partial z}& \frac{\partial N_i}{\partial y}\\
\frac{\partial N_i}{\partial z}& 0& \frac{\partial N_i}{\partial x}\\
\end{matrix} \right]
令 Jacobi 矩阵
J=\left[ \begin{matrix}{}
\frac{\partial x}{\partial \xi}& \frac{\partial y}{\partial \xi}& \frac{\partial z}{\partial \xi}\\
\frac{\partial x}{\partial \eta}& \frac{\partial y}{\partial \eta}& \frac{\partial z}{\partial \eta}\\
\frac{\partial x}{\partial \zeta}& \frac{\partial y}{\partial \zeta}& \frac{\partial z}{\partial \zeta}\\
\end{matrix} \right] =\left[ \begin{matrix}{}
\sum_{i=1}^n{\frac{\partial N_i}{\partial \xi}x_i}& \sum_{i=1}^n{\frac{\partial N_i}{\partial \xi}y_i}& \sum_{i=1}^n{\frac{\partial N_i}{\partial \xi}z_i}\\
\sum_{i=1}^n{\frac{\partial N_i}{\partial \eta}x_i}& \sum_{i=1}^n{\frac{\partial N_i}{\partial \eta}y_i}& \sum_{i=1}^n{\frac{\partial N_i}{\partial \eta}z_i}\\
\sum_{i=1}^n{\frac{\partial N_i}{\partial \zeta}x_i}& \sum_{i=1}^n{\frac{\partial N_i}{\partial \zeta}y_i}& \sum_{i=1}^n{\frac{\partial N_i}{\partial \zeta}z_i}\\
\end{matrix} \right]
则
\left[ \begin{array}{c}
\frac{\partial N_i}{\partial x}\\
\frac{\partial N_i}{\partial y}\\
\frac{\partial N_i}{\partial z}\\
\end{array} \right] =J^{-1}\left[ \begin{array}{c}
\frac{\partial N_i}{\partial \xi}\\
\frac{\partial N_i}{\partial \eta}\\
\frac{\partial N_i}{\partial \zeta}\\
\end{array} \right]
将八个结点的应变子矩阵横向拼接，得到应变矩阵
B=\left[ \begin{matrix}{}
B_1& B_2& B_3& B_4& B_5& B_6& B_7& B_8\\
\end{matrix} \right]
d. 材料刚度矩阵
D=\frac{E}{\left( 1+\nu \right) \left( 1-2\nu \right)}\left[ \begin{matrix}{l}
1-\nu& \nu& \nu& 0& 0& 0\\
\nu& 1-\nu& \nu& 0& 0& 0\\
\nu& \nu& 1-\nu& 0& 0& 0\\
0& 0& 0& \frac{1-2\nu}{2}& 0& 0\\
0& 0& 0& 0& \frac{1-2\nu}{2}& 0\\
0& 0& 0& 0& 0& \frac{1-2\nu}{2}\\
\end{matrix} \right]
这里的E为弹性模量， \nu 为泊松比
e. 单元刚度矩阵
K^e=\int_{-1}^1{\int_{-1}^1{\int_{-1}^1{B^TDB\left| J \right|d\xi d\eta d\zeta}}}
=\sum_{p=1}^{n_p}{\sum_{q=1}^{n_q}{\sum_{r=1}^{n_r}{H_pH_qH_rB^T\left( \xi _r,\eta _q,\zeta _p \right) DB\left( \xi _r,\eta _q,\zeta _p \right) \left| J\left( \xi _r,\eta _q,\zeta _p \right) \right|}}}
此处的积分点取 \xi_1=0.557350269189626,\xi_2=-0.557350269189626 ，H取1
（4）计算结构等效结点载荷列向量
a.体力产生的等效结点载荷的Gauss积分
P_{b}^{e}=\int_{-1}^1{\int_{-1}^1{\int_{-1}^1{N^Tf\left| J \right|}d\xi d\eta d\zeta}}
=\sum_{p=1}^{n_p}{\sum_{q=1}^{n_q}{\sum_{r=1}^{n_r}{H_pH_qH_rN^T\left( \xi _r,\eta _q,\zeta _p \right) f\left| J\left( \xi _r,\eta _q,\zeta _p \right) \right|}}}
其中 f 为体力向量，这里将其视为重力向量，形状为3*1
b.作用于底面上的面力产生的等效结点载荷的积分
P_{b}^{e}=\int_{-1}^1{\int_{-1}^1{N^T\bar{t}Ld\eta d\zeta}}
=\sum_{p=1}^{n_p}{\sum_{q=1}^{n_q}{H_pH_qN^T\left( \xi _p,\eta _q,C \right) \bar{t}L\left( \xi _p,\eta _q,C \right)}\ \left( \zeta =C\ ,C=1,-1 \right)}
其中 \bar{t} 为面力向量，这里指底部的约束力，形状为3*1。L为面积分变换系数，是一个标量
L=\left[ \left( \frac{\partial y}{\partial \eta}\frac{\partial z}{\partial \zeta}-\frac{\partial y}{\partial \zeta}\frac{\partial z}{\partial \eta} \right) ^2+\left( \frac{\partial z}{\partial \eta}\frac{\partial x}{\partial \zeta}-\frac{\partial z}{\partial \zeta}\frac{\partial x}{\partial \eta} \right) ^2+\left( \frac{\partial x}{\partial \eta}\frac{\partial y}{\partial \zeta}-\frac{\partial x}{\partial \zeta}\frac{\partial y}{\partial \eta} \right) ^2 \right] ^{1/2}\ \left( \xi =\pm 1 \right)
代入 Jacobi 矩阵即可
（5）组集单位刚度矩阵和结构等效结点载荷列向量
计算出各个单元的刚度矩阵和等效节点载荷后，我们需要进行组集，这个过程可以理解为将所有单元相互联系起来，使其相互影响。
a.平面问题的组集
为了描述方便，这里仅对平面三角形单元的组集方式进行说明
对于平面问题，n个结点总共有2n个自由度，离散结构所有结点的位移列向量为
a=\left[ u_1\ v_1\ u_2\ v_2\ \cdots \ u_n\ v_n \right] ^T
=\left[ \boldsymbol{a}_1\ \boldsymbol{a}_2\ \cdots \ \boldsymbol{a}_{\boldsymbol{n}} \right] ^T
若某单元3个结点1，2，3对应的总体结点编号为i，m，j，且总体节点编号i<m<j，则该单元的结点位移可以表达为
a^e=Ga
其中
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \begin{matrix}{}
1& \cdots& i& \cdots& j& \cdots& m& \cdots& n\\
\end{matrix}
G_{6*2n}=\left( \begin{matrix}{}
0& \cdots& I& \cdots& 0& \cdots& 0& \cdots& 0\\
0& \cdots& 0& \cdots& 0& \cdots& I& \cdots& 0\\
0& \cdots& 0& \cdots& I& \cdots& 0& \cdots& 0\\
\end{matrix} \right)
I=\left[ \begin{matrix}
1& 0\\
0& 1\\
\end{matrix} \right]
令
K=\sum_{e=1}^{M_e}{\left( G^TK^eG \right)}
P=\sum_{e=1}^{M_e}{\left[ G^T\left( P_{b}^{e}+P_{s}^{e} \right) \right]}
这里的K和P分别称为结构总体刚度矩阵和结构总体结点载荷列向量
可得整体结构的有限元平衡方程
Ka=P
b.立体问题的组集
对于立体问题，n个结点总共有3n个自由度，离散结构所有结点的位移列向量为
a=\left[ u_1\,\,v_1\,\,w_1\ u_2\,\,v_2\ w_2\,\,\cdots \,\,u_n\,\,v_n\ w_n \right] ^T
=\left[ \boldsymbol{a}_1\ \boldsymbol{a}_2\ \cdots \ \boldsymbol{a}_{\boldsymbol{n}} \right] ^T
每个单元有八个结点，按平面问题的组集方式操作即可
(6) 引入位移边界条件
对于上文提到的结构，只需要固定底部，即z=0处的位移为0
如果给定 j 结点自由度的位移为零，即位移边界条件 a_j=0 ，则可以采用对角元素改1的方法。具体方法是将结构有限元平衡方程改写成如下形式
\begin{matrix}{}
1& \ \ \ 2& \ \cdots& j& \cdots& \ \ n\\
\end{matrix}
\begin{array}{l}
1\\
2\\
\vdots\\
j\\
\vdots\\
n\\
\end{array}\left[ \begin{matrix}{}
K_{11}& K_{12}& \cdots& 0& \cdots& K_{1n}\\
K_{21}& K_{22}& \cdots& 0& \cdots& K_{2n}\\
\vdots& \vdots& & \vdots& & \vdots\\
0& 0& \cdots& 1& \cdots& 0\\
\vdots& \vdots& & \vdots& & \vdots\\
K_{n1}& K_{n2}& \cdots& 0& \cdots& K_{nn}\\
\end{matrix} \right] \left[ \begin{array}{l}
a_1\\
a_2\\
\vdots\\
a_j\\
\vdots\\
a_n\\
\end{array} \right] =\left[ \begin{array}{l}
P_1\\
P_2\\
\vdots\\
0\\
\vdots\\
P_n\\
\end{array} \right]
式中 n 为离散结构的结点自由度总数。将刚度矩阵的对角元素 K_{jj} 改为1,而将第 j 行和第 j 列的其他所有元素改为零，同时将载荷列向量中的 P_j 改为零即可
(7) 求解线性方程组
求解方程组
Ka=P
求解线性方程组的方法有很多，如Gauss消元法，迭代法等，这里不再赘述
(8) 计算单元应变和应力
计算出位移列阵 a 后，自然可以得知各单元结点位移 a^e
a. 计算单元中任意一点的应变
\left[ \begin{array}{l}
\varepsilon _x\left( x_q,y_q,z_q \right)\\
\varepsilon _y\left( x_q,y_q,z_q \right)\\
\varepsilon _z\left( x_q,y_q,z_q \right)\\
\gamma _{xy}\left( x_q,y_q,z_q \right)\\
\gamma _{yz}\left( x_q,y_q,z_q \right)\\
\gamma _{zx}\left( x_q,y_q,z_q \right)\\
\end{array} \right] =\left[ \begin{matrix}{}
B_1\left( \xi _q,\eta _q,\zeta _q \right)& B_2\left( \xi _q,\eta _q,\zeta _q \right)& \cdots& B_8\left( \xi _q,\eta _q,\zeta _q \right)\\
\end{matrix} \right] \left[ \begin{array}{c}
a_{1}^{e}\\
a_{2}^{e}\\
a_{2}^{e}\\
a_{2}^{e}\\
\end{array} \right]
其中
x_q=\sum_{i=1}^8{N_i\left( \xi _q,\eta _q,\zeta _q \right) x_i},y_q=\sum_{i=1}^8{N_i\left( \xi _q,\eta _q,\zeta _q \right) y_i},z_q=\sum_{i=1}^8{N_i\left( \xi _q,\eta _q,\zeta _q \right) z_i}
进一步按下式可计算单元中任意一点的应力：
\left[ \begin{array}{l}
\sigma _x\left( x_q,y_q,z_q \right)\\
\sigma _y\left( x_q,y_q,z_q \right)\\
\sigma _z\left( x_q,y_q,z_q \right)\\
\tau _{xy}\left( x_q,y_q,z_q \right)\\
\tau _{yz}\left( x_q,y_q,z_q \right)\\
\tau _{zx}\left( x_q,y_q,z_q \right)\\
\end{array} \right] =\left[ \begin{array}{l}
\varepsilon _x\left( x_q,y_q,z_q \right)\\
\varepsilon _y\left( x_q,y_q,z_q \right)\\
\varepsilon _z\left( x_q,y_q,z_q \right)\\
\gamma _{xy}\left( x_q,y_q,z_q \right)\\
\gamma _{yz}\left( x_q,y_q,z_q \right)\\
\gamma _{zx}\left( x_q,y_q,z_q \right)\\
\end{array} \right]
一般只需计算Gauss积分点和结点的应变和应力
三. 总结
这篇文章介绍了一种使用有限元法对物体结构进行分析的方法，其中的数学原理和物理定律不再详细说明。如果您对这些内容感兴趣或想处理更复杂的应力分析问题，可以参考严波老师的《有限元法基础》。未来如果有机会，我会补充有限单元法的Python代码。作为一名非土木等工科专业的计算机专业学生，我对一些专业知识的理解还不够深入，如果文章中有错误，欢迎专业人士指正。
参考文献：
严波 有限单元法基础 北京 高等教育出版社, 2022. Print.
王彪 工程力学 合肥 中国科学技术大学出版社, 2008. Print.
徐芝纶 弹性力学简明教程 #31532;5 & #29256; 北京 高等教育出版社, 2018. Print.
