【全网最详细FEM介绍】有限元基础教程 | 平面高阶单元
线性等参单元在受弯过程中，常常出现“剪切自锁”现象，结果精度极度下降，除了前文提及的增加形函数二次项的方法（非协调单元）、选择性减缩积分技术、增加单元密度，采用高阶单元也是常用的方法。
本文分享的是平面高阶单元的单元刚度计算过程，并结合具体数值案例进一步分析，与前几种方法进行精度对比。由于Matlab的编程门槛较低，本文基于Matlab平台进行数值实现，读者可在此基础上使用其他编程语言数值实现。
主要内容有：
6节点二次三角形等参单元 单元介绍 有限元格式 单元刚度矩阵计算 8节点四边形等参单元 单元介绍 有限元格式 单元刚度矩阵计算 数值案例 案例一 案例二 附录 常见单元局部编码参考
6节点二次三角形单元LST
单元介绍
如下图所示的6节点二次三角形单元又称线性应变三角形单元(Linear Strain Triangle, LST)，是在常应变三角形基础上每条边的中点增加了一个节点，可以使单元边界发生弯曲，单元边界可以是直线，也可以是曲线。图中为了表达单元可以表示弯曲特征，特意画的曲线边界。
 |  |
|---|---|
| 物理坐标系下的LST单元 | 自然坐标系下的LST单元 |
LST单元共有6个节点，每个节点共有2个平动自由度，记作 ，相应的有12个节点力，记作 ：
有限元格式
接下来会带着大家再次熟悉有限元基本格式，在以后的有限元问题描述上我都会不断的串讲有限元基本格式，加深印象。因为格式真的很固定，特别适合计算机处理，当然这不代表有限元简单，只不过很多内容不需要我们去推倒，直接拿来用就行（应用型教学）。以下的理论读者不要觉得枯燥，是和编程直接相关的理论基础，本书将与编程无关的理论均剔除，在后面的程序中，以下的理论公式均有所体现。
形函数
单元节点的局部编码顺序影响着形函数的结构，以下给出两种常见的有限元节点局部编码形式，如下图所示。
 |  |
|---|---|
| LST单元局部编码顺序1 | LST单元局部编码顺序2 |
左图表示的局部编码顺序为逆时针逐个节点节点进行编码，对应的形函数为：
右图的局部编码顺序为逆时针先对角点的节点进行编码再对单元边界的中点进行编码，对应的形函数为：
Abaqus采用的局部编码顺序为右图所示，本次的程序也按照顺序2进行编排。
常见单元局部编码顺序可参考Gmsh帮助文档Section10.2 Node ording的或者Abaqus帮助文档Theory-Elements-Continuum elements，文末的附录中给出常见单元的局部编码形式。
位移场
单元内某一点的位移，可通过形函数与节点位移描述。单元的位移场可以近似表达为：
矩阵表达形式为：
简写：
形函数 也可以反映母单元和子单元的映射关系，设母单元内部一点的坐标为 ，可由形函数插值得到：
雅可比矩阵可表示为：
不要被上面的表达式所吓倒哈！其实在数值编程的过程中，可以把它换成另一种形式：
如此以来是不是简单许多，左侧项为形函数的偏导矩阵，右侧项为节点坐标。
function[JacobianMatrix,XYDerivatives] = Jacobian(nodeCoordinates,naturalDerivatives)
JacobianMatrix = nodeCoordinates'*naturalDerivatives;
XYDerivatives = naturalDerivatives/JacobianMatrix;
end
几何方程
单元内某一点的应变，可通过几何矩阵B与节点位移描述。
其中：
再次展开：
以上的公式中体现的是形函数对物理坐标轴的偏导，这时候就可以用雅可比矩阵进行偏导矩阵转换：
进一步变换：
本构方程
针对线弹性各向同性材料，平面应力状态下：
平面应变状态：将上式中 的换为 即可。
单元刚度矩阵
在对三角形单元区域内进行积分时，可采用Hammer数值积分技术，单元体积上的积分可表示为：
上式中， 为单元厚度， 为Hammer积分点的数量， 为积分点的权重系数。常用Hammer积分点位置及权重系数汇总见下表：
| 阶次 | 积分点数目 | \xi | \eta | W |
|---|---|---|---|---|
| 1 | 1 | 1/3 | 1/3 | 1/2 |
| 2 | 3 | 1/2 | 1/2 | 1/6 |
| 0 | 1/2 | 1/6 | ||
| 1/2 | 0 | 1/6 | ||
| 2 | 3 | 1/6 | 1/6 | 1/6 |
| 2/3 | 1/6 | 1/6 | ||
| 1/6 | 2/3 | 1/6 |
下图中示意了三角形单元中Hammer积分点的位置，二阶单元中有两种积分点分布形式，本文取后者的分布形式。
 |
|---|
| 三角形单元Hammer积分点的位置 |
8节点二次四边形单元Q8
单元介绍
如下图所示的8节点二次四边形单元又称Q8单元，在Abaqus的单元号中有CPS8、CPS8R、CPE8、CPE8R都表示Q8单元。和LST单元同样的道理，在Q4单元的基础上每条边的中点增加一个节点，可以使得单元边界发生弯曲，在模拟有曲线边界的结构时非常有优势。
下图给出了Q8单元分别在物理坐标系和自然坐标系下的单元形态，为了可以表示弯曲特征，特意将单元边界更加曲线化。
 |  |
|---|---|
| 物理坐标系下的Q8单元 | 自然坐标系下的Q8单元 |
Q8单元共有8个节点，每个节点共有2个平动自由度，记作 ，相应的有16个节点力，记作 ：
有限元格式
接下来会带着大家再次熟悉有限元基本格式，在以后的有限元问题描述上我都会不断的串讲有限元基本格式，加深印象。因为格式真的很固定，特别适合计算机处理，当然这不代表有限元简单，只不过很多内容不需要我们去推倒，直接拿来用就行（应用型教学）。以下的理论读者不要觉得枯燥，是和编程直接相关的理论基础，本书将与编程无关的理论均剔除，在后面的程序中，以下的理论公式均有所体现。
形函数
Q8单元的局部编码如下图所示，对应的形函数为：
 |
|---|
| Q8单元局部节点编码 |
在Q8单元的基础上还可以在单元中心额外增加一个节点，转变为Q9单元，具体的形函数可以表示为：
位移场
单元内某一点的位移，可通过形函数与节点位移描述。单元的位移场可以近似表达为：
简写：
形函数 也可以反映母单元和子单元的映射关系，设母单元内部一点的坐标为 ，可由形函数插值得到：
雅可比矩阵可表示为：
不要被上面的表达式所吓倒哈！其实在数值编程的过程中，可以把它换成另一种形式：
如此以来是不是简单许多，左侧项为形函数的偏导矩阵，右侧项为节点坐标。
function[JacobianMatrix,XYDerivatives] = Jacobian(nodeCoordinates,naturalDerivatives)
JacobianMatrix = nodeCoordinates'*naturalDerivatives;
XYDerivatives = naturalDerivatives/JacobianMatrix;
end
几何方程
单元内某一点的应变，可通过几何矩阵B与节点位移描述。
其中：
再次展开：
以上的公式中体现的是形函数对物理坐标轴的偏导，这时候就可以用雅可比矩阵进行偏导矩阵转换：
进一步变换：
本构方程
针对线弹性各向同性材料，平面应力状态下：
平面应变状态：将上式中 的换为 即可。
单元刚度矩阵
在对单元区域内进行积分时，可采用高斯——勒让德数值积分技术，至于为什么选用这个积分方法，我的猜想是与勒让德多项式积分区间在[-1,1]之内有关，个人猜想哈。
上式中， 为单元厚度， 为高斯积分点的数量， 为高斯点的权重系数。常用高斯积分点位置及权重系数汇总见下表：
| 积分点 | 积分点坐标 | 权重系数 |
|---|---|---|
| 1 | 0.00000000000000 | 2.00000000000000 |
| 2 | 0.577350269189626 | 1.00000000000000 |
| 3 | 0.774596669241483 | 0.55555555555556 |
| 0.00000000000000 | 0.88888888888889 |
代码修改
单元刚度矩阵遍历方式
在前面章节的单元刚度计算中，给出了每个积分点的具体 位置及对应的权重系数，采用一层循环即可计算出单元刚度矩阵，这样的处理方法更为直观，但是在积分点的函数中需要写入每一个积分点的位置，平面单元还哈，如果到了空间单元中，积分点数量更多，写起来也是比较繁琐的，如果可以在网上直接copy的话那更好。
对于高斯积分点循环过程中，可以采用双层循环（平面单元）或三层循环（空间单元），这样的处理方式中，积分点函数可以写的更为简洁一些，只需写入一维积分的情况即可。
functionsamp = gauss(ngp)
% 一维情况下的高斯积分点位置
samp=zeros(ngp,2);
if ngp==1
samp=[0.2];
elseif ngp==2
samp=[1./sqrt(3) 1.;
-1./sqrt(3) 1.];
elseif ngp==3
samp= [.2*sqrt(15.) 5./9;
0.8./9.;
-.2*sqrt(15.) 5./9];
elseif ngp==4
samp= [0.8611363115940530.347854845137454;
0.3399810435848560.652145154862546;
-0.3399810435848560.652145154862546;
-0.8611363115940530.347854845137454];
end
end
function[weights,locations] = gauss(option)
% 平面单元所有的积分点位置
switch option
case'third'
locations = [-0.774596669241483-0.774596669241483;
0.-0.774596669241483;
0.774596669241483-0.774596669241483;
-0.7745966692414830.;
0.0.;
0.7745966692414830.;
-0.7745966692414830.774596669241483;
0.0.774596669241483;
0.7745966692414830.774596669241483];
weights = [0.555555555555556*0.555555555555556;
0.555555555555556*0.888888888888889;
0.555555555555556*0.555555555555556;
0.888888888888889*0.555555555555556;
0.888888888888889*0.888888888888889;
0.555555555555556*0.888888888888889;
0.555555555555556*0.555555555555556;
0.555555555555556*0.888888888888889;
0.555555555555556*0.555555555555556];
case'complete'
locations = [ -0.577350269189626-0.577350269189626;
0.577350269189626-0.577350269189626;
0.5773502691896260.577350269189626;
-0.5773502691896260.577350269189626];
weights = [1;1;1;1];
case'reduced'
locations = [00];
weights = 4;
end
end
后处理代码修改
由于单元类型发生了转变，在vtk输出函数中也应做出相应改变，参考VTK语法对于二次单元的单元号定义：
% Q8单元
fprintf(vtk_file, 'CELLS %d %d\n', num_elements, num_elements * 9);
fori = 1:num_elements
fprintf(vtk_file, '8 %d %d %d %d %d %d %d %d\n', element_info(i, 1) - 1, element_info(i, 2) - 1,element_info(i, 3) - 1,element_info(i, 4) - 1, ...
element_info(i, 5) - 1, element_info(i, 6) - 1,element_info(i, 7) - 1,element_info(i, 8) - 1);
end
fprintf(vtk_file, 'CELL_TYPES %d\n', num_elements);
fori = 1:num_elements
fprintf(vtk_file, '23\n');
end
% LST单元
fprintf(vtk_file, 'CELLS %d %d\n', num_elements, num_elements * 7);
fori = 1:num_elements
fprintf(vtk_file, '6 %d %d %d %d %d %d\n', element_info(i, 1) - 1, element_info(i, 2) - 1,element_info(i, 3) - 1,element_info(i, 4) - 1, ...
element_info(i, 5) - 1, element_info(i, 6) - 1);
end
fprintf(vtk_file, 'CELL_TYPES %d\n', num_elements);
fori = 1:num_elements
fprintf(vtk_file, '22\n');
end
数值案例
案例一
考虑下图所示的悬臂梁模型，模型右端下方节点固定，上方节点仅约束 方向的位移，在模型左端下方节点施加竖直向下的集中力，弹性模量 ，泊松比 ，单元厚度为2，按照平面应力状态处理，单元细长比取10。
此时的有限元解精度很差，出现"剪切自锁"现象；使用力偶荷载构造纯弯载荷工况可以自行试验，规律是一致的，都容易出现"剪切自锁"。
材料参数与约束信息：
*Material
10000., 0.0,2.0,1
*Load
1,2,-0.1
*Constr
11,1,0.
11,2,0.
12,1,0.
33,1,0.
34,1,0.
精度对比
分别使用LST单元和Q8单元进行有限元分析，厚度方向上均使用一个单元，长度方向使用10个单元进行网格划分。
LST单元
将LST单元的位移场结果与Abaqus位移场结果进行对比，吻合度极高，验证了程序的正确性。
 |  |
|---|---|
| 自研程序-Umag | Abaqus-Umag |
 |  |
| 自研程序-U1 | Abaqus-U1 |
 |  |
| 自研程序-U2 | Abaqus-U2 |
Q8单元
将Q8单元的位移场与Abaqus的位移场进行对比，吻合度极高，验证了程序的正确性。在应力场对比时，需要将Abaqus的应力磨平选项中的AVG参数调整至100%，自研程序与Abaqus的应力场吻合度极高。
 |  |
|---|---|
| 自研程序-Umag | Abaqus-Umag |
 |  |
| 自研程序-Sxx | Abaqus-Sxx |
 |  |
| 自研程序-Syy | Abaqus-Syy |
 |  |
| 自研程序-Sxy | Abaqus-Sxy |
总结
围绕着上述有限元模型，将CPS4、CPS4I（非协调单元）、CPS4SR（选择性减缩积分）、CPS6、CPS8单元计算的结果整理于下表：
| 单元类型 | CPS4 | CPS4I | CPS4SR | CPS6 | CPS8 |
|---|---|---|---|---|---|
| 最大挠度 | 0.1130 | 0.7392 | 0.7392 | 0.7397 | 0.7411 |
从表中可以看出，在线性等参单元出现剪切自锁时，可采取措施：
在原有形函数基础上增加满足弯曲特征的二次项，采用非协调单元 单元刚度矩阵计算时采用“选择性减缩积分技术” 采用高阶单元
当然也可以适度增加线性单元的网格密度，效果不一定显著。
案例二
现考虑下图所示的4点弯曲简支深梁进行有限元分析，根据对称性，可以对模型的一半结构进行分析。网格离散形式为 ，单元类型采用Q8单元。该算例取自Amar教授的经典著作：《Introduction to finite element analysis using matlab and abaqus》
 |
|---|
| 数值模型几何参数和边界条件（长度单位：毫米） |
 |
|---|
| 模型有限元分析（长度单位：毫米） |
材料属性及边界条件信息：
*Material
弹性模量，泊松比，厚度，平面应力标志（平面应变为2）
40000., 0.17,100.0,1
*Load
节点号，自由度号，加载力
42,2,-30000
*Constr
节点号，自由度号，施加位移
2,2,0.0
9,1,0.0
18,1,0.0
27,1,0.0
36,1,0.0
45,1,0.0
69,1,0.0
86,1,0.0
103,1,0.0
120,1,0.0
精度对比
将Q8单元的位移场与Abaqus的位移场进行对比，吻合度极高，验证了程序的正确性。在应力场对比时，需要将Abaqus的应力磨平选项中的AVG参数调整至100%，自研程序与Abaqus的应力场吻合度较高。
 |  |
|---|---|
| 自研程序-Umag | Abaqus-Umag |
 |  |
| 自研程序-Sxx | Abaqus-Sxx |
 |  |
| 自研程序-Syy | Abaqus-Syy |
 |  |
| 自研程序-Sxy | Abaqus-Sxy |
对于本次案例中的深梁，按照工程中常见的细长梁理论来计算跨中挠度，跨中挠度计算公式为：
根据以上公式，计算得到的跨中挠度为0.12mm，小于有限元程序中跨中处的竖向位移0.1551mm。在工程梁理论的计算公式中，不考虑剪切应力对梁变形的影响，导致略小于程序中的计算结果。
应力云图显示梁的上部受压，下部受拉，中心各种处几乎不受力；剪应力 云图中可以看出在铰支座附近和集中荷载之间的区域内剪应力较大，其他地方的剪应力可以忽略不计，符合力学常识。
附录
杆系单元
Line: Line3: Line4:
v
^
|
|
0-----+-----1 --> u 0----2----1 0---2---3---1
平面三角单元
Triangle: Triangle6: Triangle9/10:
v
^
|
2 2 2
|`\ |`\ | \
| `\ | `\ 7 6
| `\ 5 `4 | \
| `\ | `\ 8 (9) 5
| `\ | `\ | \
0----------1--> u 0-----3----1 0---3---4---1
Triangle12/15:
v
^
|
2
| \
9 8
| \
10 (14) 7
| \
11 (12) (13) 6
| \
0---3---4---5---1--> u
平面四边形单元
Quadrangle: Quadrangle8: Quadrangle9:
v
^
|
3-----------2 3-----6-----2 3-----6-----2
| | | | | | |
| | | | | | |
| +---- | --> u 7 5 7 8 5
| | | | | |
| | | | | |
0-----------1 0-----4-----1 0-----4-----1
四面体单元
Tetrahedron: Tetrahedron10:
v
.
,/
/
2 2
,/|`\ ,/|`\
,/ | `\ ,/ | `\
,/ '. `\ ,6 '. `5
,/ | `\ ,/ 8 `\
,/ | `\ ,/ | `\
0-----------'.--------1 --> u 0--------4--'.--------1
`\. | ,/ `\. | ,/
`\. | ,/ `\. | ,9
`\. '. ,/ `7. '. ,/
`\. |/ `\. |/
`3 `3
`\.
` w
六面体单元
Hexahedron: Hexahedron20: Hexahedron27:
v
3----------2 3----13----2 3----13----2
|\ ^ |\ |\ |\ |\ |\
| \ | | \ | 15 | 14 |15 24 | 14
| \ | | \ 9 \ 11 \ 9 \ 20 11 \
| 7------+---6 | 7----19+---6 | 7----19+---6
| | +-- |-- | -> u | | | | |22 | 26 | 23|
0---+---\--1 | 0---+-8----1 | 0---+-8----1 |
\ | \ \ | \ 17 \ 18 \ 17 25 \ 18
\ | \ \ | 10 | 12| 10 | 21 12|
\| w \| \| \| \| \|
4----------5 4----16----5 4----16----5
参考文献
https://gmsh.info/doc/texinfo/gmsh.html#High_002dorder-elements https://help.3ds.com/2024/english/dssimulia_established/SIMACAETHERefMap/simathe-c-tritetwedge.htm?contextscope=all&id=3947a13adffd4720ae102aa71bcf54f0 Ferreira A J M. MATLAB codes for finite element analysis[M]. Amsterdam: Springer Netherlands, 2009. Khennane A. Introduction to finite element analysis using MATLAB® and abaqus[M]. CRC Press, 2013.
