为什么三角形单元你的应力分析总是不准？Abaqus二阶三角形单元揭秘，精度直接翻倍！
在本文中，我们将详细介绍静态结构分析中使用的二阶三角形单元的公式。此类单元通常用于有限元法 (FEM) 中来解决工程问题。
我们将考虑以下假设：
平面应力或平面应变条件
线弹性、各向同性、均质材料。
单元仅在XY平面上加载。
小变形条件。
介绍
二阶三角单元的每个单元具有6个节点和每个节点2个自由度（DOF），总共12 doF。因此，我们已经可以确保该元素的刚度矩阵将是12×12矩阵。该矩阵也将是对称和奇异的。
用于结构分析的元素的自由度是位移。对于二阶三角元素，每个节点有两个可能的位移：一个在x方向上，一个在y方向上，如下所示：
所需的材料特性仅弹性模量（E）和泊松比（ν）。
单元几何和形状函数
对于二阶三角形单元，节点位于边的顶点和中点。xy平面上的节点坐标可以表示如下：
制定任何有限元的初始步骤之一是定义计算单元内任何位置处的位移分量的函数。最简单的函数类型是多项式。要定义多项式，必须确定其所需的系数数量。确定每个多项式所需系数数量的方法如下所示。
为了帮助构造多项式，我们将使用帕斯卡三角形。
因此，位移函数U如下：
可以使用矩阵表示法将其写为：
或者以更紧凑的方式：
并且，位移函数v如下：
同样可以写为
或者以更紧凑的方式：
请注意，以上方程包括二阶（二次）项。这就是为什么此单元称为二阶有限单元的原因。
为了确定系数 ，节点位移可以用作边界条件，从而形成12个代数方程。由于有12个系数可以计算，因此该系统将具有独特的解决方案。这些方程式使用矩阵表示法表示。
上面的矩阵方程可以更简洁地表达如下：
为了计算系数 ，应在上面的方程式中隔离列矩阵{a}。因此：
因此，位移函数可以计算如下：
可以执行上面的矩阵运算，生成以下位移函数表达式：
其中 是插值函数或形状函数。注意：
其中矩阵[v]包含变量x和y，矩阵[p]包含节点坐标，因此其逆也由节点坐标组成。因此，形状函数将由节点坐标（创建网格后建立）和变量x和y定义。
应变-位移矩阵（B 矩阵）
对于平面应力或平面应变条件，应变的三个分量可以计算如下：
现在可以使用先前开发的位移函数方程来计算这些量。请注意，由于变量x和y仅出现在形状函数 中，因此只需要区分这些函数。
上述表达式可以写成矩阵形式：
B矩阵是形状函数相对于X和Y的微分构成：
关于应变-位移矩阵（B 矩阵）的重要结论
了解应变分量如何随单元区域变化对于任何有限元公式都至关重要。
例如，在单元的整个区域中，线性三角单元（也称为恒定应变三角形（CST））具有恒定应变和应力。在压力和应变梯度高的区域中，这种特征是不可取的，因为它需要使用许多单元，从而大大增加了计算成本。在当前的公式中，我们具有二阶有限元。我们已经开发了方程来计算应变成分作为形状函数的微分。但是，形状函数并未明确得出。仅显示用于得出它的表达式。因此，无法得出结论如何变化。
因此，我们现在将以不同的方式计算应变分量，以说明它们在单元区域内的变化情况。
使用初始形式的位移函数，我们有以下表达式：
通过对这些函数进行微分，我们可以计算出应变分量，如下所示：
这些表达式的关键结论是，二阶三角单元的应变成分在x和y方向上线性变化。与恒定应变三角形相比，这是一个显着改善。这种线性变化允许在压力和应变高梯度的区域使用较少的单元，从而降低计算成本。
本构矩阵（D 矩阵）
二阶三角形单元的应力分量可以使用以下表达式计算：
对于平面应力条件，本构矩阵D定义为：
对于平面应变条件，本构矩阵D定义为：
应该注意的是，平面应力和平面应变条件之间二阶三角单元的唯一差异在于本构矩阵D的计算。
另一个需要注意的要点是，由于应力是通过将本构矩阵 D 乘以应变矢量来计算的，并且应变分量在 x 和 y 方向上线性变化，因此应力分量也将在 x 和 y 方向上线性变化。
刚度矩阵 [k]
现在我们已经定义了应变位移矩阵（B 矩阵）和本构矩阵（D 矩阵），我们可以评估二阶三角形单元的刚度矩阵。首先，我们需要计算单元的总应变能。
总应变能
计算总应变能很重要，因为它使我们能够使用以下两种方法之一来评估单元刚度矩阵：
卡斯蒂利亚诺第一定理
总势能最小原理
这两种方法都需要单元的总应变能，可以使用以下表达式计算：
总势能
最小总势能的原理将用于评估刚度矩阵K。因此，我们首先需要计算单元的总势能，这可以使用以下表达式完成：
这里，W 表示与单元内所有可能的节点力相关的功。这些力如下图所示：
因此，W可以计算为：
这里
因此，单元的总势能可以计算如下：
总势能最小原理
在稳定平衡条件下，物体的总势能最小化。为了找到最小值，我们将总势能函数求微分并使之为零。
对于二阶三角单元，我们需要根据总势能函数对节点位移求微分，并将这些微分等于零，如下所示。
这些导数将产生 12 个代数方程，可以用矩阵表示法表示如下：
其中 是二阶三角单元的刚度矩阵。
重要的是，这组12个代数方程表示x和y方向的每个节点处的平衡条件。例如，第一个方程对应于x方向上节点1的静态平衡，而第二个方程与y方向上节点1的静态平衡有关，依此类推。
可以证明，刚度矩阵是使用以下表达式计算的，该表达式是从这 12 个导数导出的。
正如预期的那样，刚度矩阵是一个 12×12 矩阵，是奇异且对称的。
为了评估这种积分，通常使用数值集成，特别是在单元区域上应用的高斯正交过程。
结论
这篇博文详细介绍了静态结构分析中使用的二阶三角形有限元的公式。有限元法 (FEM) 中使用的单元假设平面应力或应变条件、线弹性、各向同性、均质材料、xy 平面载荷和小位移。
要点
元素特征：
6 个节点，每个节点 2 个自由度 (DOF)，总共 12 个 DOF。
刚度矩阵： ，对称奇异。
位移功能：
使用多项式定义 x 和 y 形函数。
以矩阵形式表示的函数：
应变 - 置换矩阵（B-Matrix）：
从位移函数导出的应变分量：
以矩阵形式表示：
应变变化：
二阶单元提供应变的线性变化，与恒定应变三角形相比提高了计算效率。
本构矩阵（D-matrix）：
适用于平面应力和平面应变条件。
使用： 计算的应力分量
刚度矩阵（[k]）：
通过总应变能或总势能计算。
积分表达式：
使用数值积分求解，通常是高斯积分法。
等参公式
等级公式也可以应用于此元素。这种方法的优点是它使用相同的形状函数来插值元素的几何形状和位移。
由于二阶三角元素的形状函数是二次函数，因此使用等法公式可以使三角形的侧面弯曲，符合二阶多项式。
该特征对于表示域的弯曲部分（例如孔）特别有用。准确地表示几何形状，从而降低计算成本所需的元素更少。
下图说明了使用等参公式制定的元素与按照此处描述的方式制定的元素之间的差异。
