如何手动进行 FEA – 第 2 部分：2D单元！

您将学习如何手动对 2D 单元进行 FEA。但正如我总是指出的那样——要深思熟虑为什么要学习这个。

        在之前的文章中，描述了在简单一维杆单元的手工计算中使用 FEM 的方法。这些例子很好地展示了有限元法的基本概念。这就是为什么它们通常构成有关 FEA 的书籍的第一章。我们从那里学习了如何计算单元刚度矩阵，将它们组装成全局刚度矩阵，形成力和位移向量，然后求解线性代数方程以获得位移，从中我们可以获得应力和应变等量。

        然而，当您学习 FEA 的理论方面并解决此类示例时，在某些时候您会开始想知道：

微分方程在哪里？

        毕竟，有限元法是解决这些方程的一种方法！那么，我们怎样才能在不使用偏微分方程（PDE）的情况下解决任何问题呢？也许更重要的是，FEA 程序如何做到这一点？

        人们普遍误解这些程序直接使用偏微分方程求解。实际上，他们通常使用所谓的弱形式（太复杂和无聊，无法在这里详细讨论），这是从控制选定物理问题的偏微分方程获得的。更重要的是，他们也不直接使用这种弱形式。为了解释如何做到这一点，让我们回顾一下上一篇文章，其中我们解决了杆单元示例。

        您可能还记得，我们从简单的力学关系中获得了刚度矩阵特性。您可以对简单的一维元素使用这种直接方法。但最常见的 2D 和 3D 单元又如何呢？

它们的刚度矩阵无法直接获得，因此可以通过另一种方式计算。这就是需要偏微分方程的地方。当我们将偏微分方程、本构方程和形函数结合起来时，我们可以得到弱形式（如果我们使用变分方法，则称为泛函形式）。

        然后，我们提取刚度矩阵的公式以及特定单元类型的力矢量。其余的内容与我们在上一篇文章中看到的类似：

• 首先，使用公式计算所有单元的刚度矩阵

•  然后你把它们组装起来

• 接下来，计算力矢量

• 添加边界条件

• 最后，你解出方程。

这就是所有 FEA 程序在实践中的运作方式。

手头的任务…

让我们继续看示例，看看如何手动完成。我们将求解一个固定在一侧并在另一侧的单个节点处受到集中力的简单 平板（平面受力构件）：

我们将使用“结构钢”，因此杨氏模量为 210 000 MPa，泊松比为 0.3。

我们需要做的第一件事是离散化板。这意味着我们将把它分成有限单元。在这种情况下，两个三角形单元就足够了。由于矩阵的大小，如果没有计算机，更多的单元问题就很难解决。我们还对单元和节点进行编号：

数学慢慢开始

这种平面应力单元的刚度矩阵如下所示：

[k] = t A [B]^T [D] [B]

 在这里： t – 单元厚度；A——单元的面积；[B] – 应变-位移矩阵；[D] – 材料属性矩阵；上标中的 T 表示矩阵转置。

FEA 程序对单元的体积进行积分，但幸运的是，在这种具有常量参数的简单情况下，无需这样做。您可以使用节点坐标计算单元面积。幸运的是，在这个例子中，我们可以使用简单的几何关系。我们将使用如下所示的符号：

知道上述编号后，我们可以计算线性三角形单元的应变-位移矩阵：

其中 β 和 γ 是形函数系数。您可以使用节点坐标来计算它们：

这里的关键是节点总是逆时针编号。

使用相反的节点编号也可以解决该问题。但请注意，我们在这里使用的所有公式都适用于逆时针方法。FEA 软件中也进行了同样的操作。

你可能想知道，如果有多个节点编号不同于 1,2 和 3 的单元该怎么办。我们是否按照编号的顺序（例如 16,17,19）取出它们的节点？

不，我们实际上不知道。相反，我们选择任何节点作为“第一个”，然后逆时针方向将随后的两个节点分别视为上述公式中的第二个和第三个节点。

这称为全局和本地节点编号。

        全局节点数是考虑整个计算域时给出的节点数。它们对于每个节点都是唯一的（例如 16,17 和 19）。

然后，当我们计算每个单元的矩阵时，我们使用局部编号。在这种情况下，我们将数字 1、2 和 3 指定给节点，保持逆时针编号规则。

最后，我们还需要材料属性矩阵的公式。对于平面应力，由下式给出：

 在哪里：E——杨氏模量v – 是泊松比

现在我们已经有了所有必需的公式，我们可以继续解决我们的问题了。让我们从第一个单元开始，使用本地节点编号计算其矩阵：

由于上述一切，我们终于可以计算该单元的刚度矩阵：

局部刚度矩阵中的行和列对应于全局坐标中的节点自由度。

对于第一个元素的全局节点编号，按照局部编号给出的顺序，为 1、3、2，因此该刚度矩阵中包含以下自由度 (DOF)：

1（代表节点 1的x 自由度）、2（代表节点 1的y自由度）、5（ 代表节点 3的x 自由度）、6（代表节点 3的y自由度）、3（代表节点 2的x自由度）、4（ 代表节点 2的y自由度）。

好吧……我花了一些时间来真正消化这个。

如果你明白了，请跳过下一章……下面我将进一步解释刚刚发生的事情！

添加了一个小解释！

好吧，让我们一点一点地消化这些编号规则。我将简单地陈述我们已经知道的事实，并附上一些精心制作的图像，希望这能让事情更容易看到。

事实 1：我们有一个模型，它具有全局节点编号。您已经在开始时看到了这些，这里有一个简短的提醒（全局节点编号以绿色圆圈标记）：

事实 2：每个全局节点都有 2 个自由度（因为它是一个平面应力示例）。有可能会造成混淆……所以我将在全局坐标中添加一个“G”下标。整套全球坐标如下：

事实 3：在解决全局问题之前，我们将“单独”计算每个单元。为了分别求解每个单元，我们需要为每个单元分配本地节点编号。它们始终是“1、2、3”，并且按逆时针方向分配。我将使用上面的单元 1，并将本地节点编号标记为红色：

在我们继续之前……请注意一件事！当仅考虑“本地”单元 1 时，全局节点 3（左侧绿色）具有“本地”编号 2（右侧红色）。

此类编号问题将始终是一个问题。尽管只有从“1”开始对全局节点进行编号才比较棘手……因为局部节点（红色节点！）将始终按逆时针方向编号为 1、2 和 3。

事实 4：我相信您知道这是怎么回事……本地节点有其本地坐标。我们之前已经对它们进行了编号，但所有内容都集中在一处：

当然，上面的本地坐标有适合本地节点编号的下标。正如全局坐标具有适合全局节点编号的下标一样。

事实 5：现在我们准备再次看看对单元 1 的局部刚度矩阵所做的事情。当然，如果您想再次回顾一下，数学就在那里。

对我来说，棘手的部分是“行和列编号”。当然，单元刚度矩阵是在该单元的局部系统中组装的，这意味着我们最初可以这样对列进行编号：

        但这种编号方法并不是非常有用（但有序且有意义！）。它没有用，因为每个单元都有本地节点，按逆时针方向编号为 1 到 3。因此，对于每个单元，该局部刚度矩阵具有相同的行数和列数！这就是我们“返回”全局坐标的地方。我们知道全局节点号是什么，因此我们可以用“适合”全局系统的行和列“映射”上述矩阵。

        正如您在上面所看到的，在我们的示例中，节点 1 在本地坐标和全局坐标中具有相同的编号（顺便说一句，完全不必如此！）。本地节点 2（红色）是全局节点 3（绿色），本地节点 3（红色）是全局节点 2（绿色）。

这意味着，根据上面的内容，我们可以添加：

一直使用下标不太舒服。这就是为什么我们将使用“正常数字”来对全局自由度进行编号。简单的说：x ₁ 将是“1”；y ₁ 将是“2”；x ₂ 将是“3”；y ₂ 将是“4”；x ₃ 将是“5”；y ₃ 将是“6”；x ₄ 将是“7”；y ₄ 将是“8”。

以上适用于全局自由度，因此适用于上面以绿色标记的自由度。如果我们简单地用数字替换符号，我们将得到所示的刚度矩阵：

好吧……我希望现在编号已经清楚了:)

现在让我们再次计算第二个单元分配本地编号的矩阵。材料属性矩阵将与第一个元素的材料属性矩阵相同，因为它们由相同的材料制成：

再次感谢上面的内容，我们可以构建第二个单元的刚度矩阵，如下所示：

对于第二个元素，以下全局编号的节点属于该元素：1、4、3。所以这里的 DOF 是：1（代表节点 1的x 自由度）、2（ 代表节点 1的y 自由度）、7（ 代表节点 4的x 自由度）、8（ 代表节点 4的y 自由度）、5（ 代表节点 3的x 自由度）、6（代表节点节点 3的y 自由度）。

希望此时绿色数字并不令人意外！

请注意，单元 2 在 GLOBAL 节点 1 中的 LOCAL 节点编号为 1。这意味着，本地节点 2（记住，始终逆时针）将位于全局节点 4 中，本地节点 3 将位于全局节点 3 中。

其余的工作就像我之前描述的那样。

下一阶段——全局刚度矩阵

由于我们拥有两个单元的刚度矩阵，因此我们可以继续下一步 - 将这些矩阵组合成全局刚度矩阵。

我们只需将适当的项相互添加并形成一个新的 8×8 矩阵即可做到这一点。例如，从第一个局部刚度矩阵中获取 25 ^th 元素（第 2 列第 5 行 - 使用基于全局自由度的编号），并将其添加到第二个局部刚度矩阵的相应元素中以获得 25全局刚度矩阵的 ^th 元素。

这是此操作的结果：

上面的这些听起来像是咒语吗？

别担心——你已经知道需要做的一切了。但为了确定起见，我们在这里也添加一个简单的解释。

到目前为止，我们有两个局部刚度矩阵（每个元素一个）。我们已经设法将这两个矩阵中的行和列“分配”到全局坐标（绿色数字）。两者看起来都是这样的：

现在，我们需要做的就是制作一个 8×8 的全局刚度矩阵。我们的模型总共有 8 个全局自由度。这就是为什么我们的全局矩阵是 8×8。这些将是列和行中从 1 到 8 的“绿色”自由度。像这样：

很酷的是，我们有一个全局刚度矩阵。糟糕的是，它是空的！幸运的是，Jakub 已经解释了如何填充它。

您所需要做的就是从同一行和同一列中获取数字（全球坐标，所以绿色数字！）并将它们加在一起。然后将结果放入全局刚度矩阵的同一行和同一列中。

如果你制作一个图形，其实很简单：

完成此操作后，我们应该确保矩阵关于其对角线对称。然后，我们可以通过划掉与这些 DOF 相关的列和行来消除受约束的 DOF (1,2,3,4)。

我们来计算一下位移吧！

现在需要求解的方程是：

[K] {u} = {f}

 在这里：{u} – 位移向量     {f} – 外力矢量

为了求解这个方程，我们求刚度矩阵的逆矩阵：

{u}=[K]^-1 {f}

上面的结果是：

– 节点 3，x 位移：u ₅ = 0.03 mm

– 节点 3，y 位移：u ₆ = -0.154 mm

– 节点 4，x 位移：u ₇ = -0.03 mm

– 节点 4，y 位移：u ₈ = -0.153 mm

最后，应力和应变！

当然，我们可以使用这些节点位移来计算应变和应力：

{ε} = [B] {u}

{σ} = [D] {ε} = [D] [B] {u}

 应力是：

•  对于第一个单元：

•  对于第二个单元：

通过将结果与使用开源 FEA 求解器 CalculiX 获得的结果进行比较，证实了结果。他们的结果完全一致。显示位移的截图如下所示：

正如您所看到的，即使在如此简单的情况下，手动计算也确实很费力。即使有两个三角形元素，矩阵也会变得相当大。却让我们可以了解 FEA 软件在线性静态平面应力分析中的实际工作原理。

如果这是一个“设计案例”，那么这些问题迫切需要答案：

        当然，我非常清楚这一点！我完全理解这篇文章的目的不是给你这些答案......而是向你展示求解器如何进行数学计算。

        但在一个很多人都觉得“ FEA 中全是数学”的世界里，这可能是我唯一的机会向您展示，设计不仅仅是 数学问题！这就是我在这里提到它的原因！