回答这个问题的最大问题是，在你看到方程完全发展之前，它并不像你想象的那么简单，直到你看到它们完全发展，即使那样你也在尽你最大的努力来分配一个含义。作为一名机械工程专业的研究生，这对我来说很难，因为我习惯于用我的物理直觉来引导我穿过数学森林。

对于应用于模态分析的 FEM 问题，您基本上从一个强形式方程开始。例如，弹性杆的自由振动由下式给出： $$ \lambda \rho u + (Eu,_x),_x = 0,\ x \in (0, L) $$ 这里，$\rho$ 是密度，$E$ 是弹性模量，$u$ 是位移。让我们应用一个固定边界条件（$u(L) = 0$）和一个自由边界条件（例如 $-Eu,_x(0)=0$）。您可以使用加权函数 $w$ 将其转换为等效的弱形式： $$ \int_0^Lw,_xEu,_xdx\ -\lambda\int_0^Lw\rho u\, dx = 0 $$数学最终失去了它的物理解释。从这里开始（并且在这个答案中为了简洁而跳过了很多严谨性），我们基本上会将 $w$ 限制在 Sobolev 空间的一个子集中，通过离散化创建有限维权重函数 $w^h$ 和解 $u^h$ 的域来应用 Galerkin 近似，使 $u^h$ 成为 Sobolev 函数 $v^h$ 和给定函数 $ 的总和g^h$，然后将函数 $w^h$ 和 $v^h$ 定义为线性独立基函数的有限和（即 $v^h = \sum_{A=1}^N c_A N_A$ 和$w^h = \sum_{B=1}^N d_A N_A$)。通过应用一些线性代数定理并大量重新排列，最终得到如下等式： $$ (K-\lambda^h_kM)\psi_k = 0 $$ 这里，$\psi_k$ 是第 k 个特征向量，代表位移。很难对我掩盖的数学表示赞赏。这本书 然后将函数 $w^h$ 和 $v^h$ 定义为线性独立基函数的有限和（即 $v^h = \sum_{A=1}^N c_A N_A$ 和 $w^h = \sum_{B=1}^N d_A N_A$)。通过应用一些线性代数定理并大量重新排列，最终得到如下等式： $$ (K-\lambda^h_kM)\psi_k = 0 $$ 这里，$\psi_k$ 是第 k 个特征向量，代表位移。很难对我掩盖的数学表示赞赏。这本书 然后将函数 $w^h$ 和 $v^h$ 定义为线性独立基函数的有限和（即 $v^h = \sum_{A=1}^N c_A N_A$ 和 $w^h = \sum_{B=1}^N d_A N_A$)。通过应用一些线性代数定理并大量重新排列，最终得到如下等式： $$ (K-\lambda^h_kM)\psi_k = 0 $$ 这里，$\psi_k$ 是第 k 个特征向量，代表位移。很难对我掩盖的数学表示赞赏。这本书 很难对我掩盖的数学表示赞赏。这本书 很难对我掩盖的数学表示赞赏。这本书有限元方法对它是如何开发的有完整的解释（第 1 章介绍了该方法的基础知识，第 7.3 节开发了频率分析部分）。

为了澄清起见，我们简单地将这些矩阵称为 K 和 M，方程的形式类似于梁的简单模态分析，因此我们简单地将其指定为物理解释。实际上，K 和 M 矩阵高度依赖于网格以及元素的连接方式。这些大矩阵由网格中的每个元素构成。

离散化真的很重要。从物理上讲，对于这个问题，节点代表我们正在计算位移的点，并且连接到该节点的元素在施加的力下限制其运动。你也可以把这个问题的网格想象成有$n$个相互连接的杆。这样，每个杆件只会真正影响它所连接的那些。每个元素由一个小矩阵表示。元素矩阵的计算如下： $$ m_{ab}^e = \int_{\Omega^e} N_a \rho N_b\,d\Omega \\ k_{ab}^e = \int_{\Omega^e } N_{a,x} E N_{b,x}\,d\Omega $$ 这里，$\Omega^e$ 是元素本身的域。当元素矩阵被组装到全局矩阵中时，大矩阵中的一些元素将从几个较小的矩阵中添加项。

$N_a$ 和 $N_b$ 符号是跨越我们选择的 Sobolev 空间子集的基函数。这些非常重要，因为不同的基函数会改变解。我们可以选择使用线性或抛物线函数，甚至是不同类型的函数，如拉格朗日甚至 B 样条基函数。

同样，我不知道 OP 有多少 FEM 理论，而且很难总结这一切是如何工作的。这可能还不够，但希望它有助于突出 FE 方法背后的一些理论。

有限元分析是一种在工程师中得到广泛应用的数学工具。然而，在对有限元分析在其中发挥如此重要作用的计算机模拟主题进行了一年多的研究之后，我还没有找到一个令人满意的解释来说明它们是如何真正起作用的……

FEM的主要背景是60年代的结构工程。工程师是非常实际的人，因此他们最初设计了一个系统，允许他们设置代数方程，其中可以设置其结构的不同点或“节点”之间的关系。这些代数关系被进一步证明是许多不同类型的，不仅是结构关系，还可以制定热关系和许多其他关系：需要分析的只是空间的适当离散化，以网格的形式与节点相关对彼此。

对于我们特定的结构工程案例，该过程可以概括为三个主要步骤，因此主要思想以图形的形式出现：

第 1 步：离散化 在这个基本但经典的示例中，我选择将我的样本结构仅划分为 5 个元素，其中节点在梁之间的交汇处清晰可辨。有四个节点：

（自由程度）

每个节点将有 6DOF（自由度）：三个用于每个轴上的线性位移 (X,Y,Z)，三个用于围绕每个轴旋转 (X,Y,Z)，因为我们正在处理 3D。许多可用的示例提供了更“简单”的 2D 情况，但在我看来，这只会使事情变得更加复杂。

一旦找到了节点并且它们之间存在相互关联的网络，我们就可以认为我们有一个网格。在我们的示例中，表中描述了相关性： 元素用于将节点“链接”到彼此。该表为节点建立了一个“拓扑”。

第 2 步：元素表征和形状函数 在这一步中，FEM 公式和文献变得非常尴尬和讨厌。事实上，这是一切的核心，也是 FEM 与其他求解 PDE 方法的不同之处。

实际上，尽管数学复杂，但我们正在寻找一种表征元素材料特性的方法。为此，该方法要求节点之间的那些链接的行为服从某个公式。这个公式是实际的形状函数。事实上，形状函数可以是任何数学公式，它可以帮助我们对没有定义网格的点进行插值。这个出现在节点之间的“幽灵”实体实际上就是有限元。实际上，作为工程师，我们对 FEM 的实现更感兴趣，而不是它的公式，所以重要的是要了解，对于不同的形状函数，我们会获得不同的元素矩阵。

2DOF 梁单元矩阵 3DOF 梁单元矩阵 3DOF Timoshenko 梁单元矩阵 6DOF Timoshenko 梁单元矩阵

根据选择的公式，我们有不同程度的插值，因此可能更高或更低的精度。此外，根据选择的公式，我们可能有不同的方式来定位和关联我们的节点。

对于我们的示例，我们可以选择文献中提供的任何公式（以上是结构工程中最常用的公式）。重要的是要注意这些元素矩阵的内部结构，它们是对称的，并且清楚地划分为部分，每个部分对应于位于元素边界上的节点（梁的情况下为 2 个节点 - 矩阵中的 4 个象限）。

第 3 步：矩阵组装和解决方案

因为节点之间的关系需要同时完成，所以我们必须将所有方程设置为构成一个代数方程组。我们要求解的矩阵方程（至少在静力学中）如下： \begin{equation} [F]= [Kg]·[u], \end{equation}

其中[F]是施加外力的矢量，[Kg]是系统的整体刚度矩阵，[u]是施加外力导致的位移矢量。[Kg]、[F]和[u]的大小是自由度数乘以节点数的大小，即[Kg]的平方和[F]和[u]一维的。为了生成力向量，我们需要做的就是收集施加的力（线性力和力矩），并根据它们所施加的节点的索引对它们进行排序。对于全局刚度矩阵，需要一个更费力的过程，通过它我们遍历每个元素的特定刚度矩阵。在每一个中，我们只得到与我们存储在矩阵中的节点位置相对应的部分，并将其添加到可能来自其他元素的并发数据中。警告：在输入全局刚度矩阵之前，我们必须将局部坐标转换为全局坐标！

为了得到位移矢量，首先需要强制约束，然后求解得到的代数方程组。要做到这一点，有两种经典方法：惩罚法和拉格朗日乘数法，但我们不会在这里讨论...

全局矩阵组装

之后，我们需要做的就是使用任何可用的代数方程求解器（LU分解是最扩展的一种），并获得系统的解。

第 4 步：计算矩阵的特征值

一旦您组装了一个矩阵并表示您的物理模型，您就可以继续对其进行其他数学运算，例如计算其特征值。特征值具有“神奇”特性，即刚度矩阵是位移矢量的一种变换形式，表示位移最大的形式。这些位移与结构系统的物理特性及其自由振动周期直接相关。

如果您熟悉标准 FEM 分析工作，那么模态分析的想法很简单。在标准 FEM 分析中，您传递瞬态弹性波方程（暂时忽略阻尼）$$ \ddot{r}(t) + A r(t) = f(t) \tag{1}\label{1 },$$ 是使用微分算子 $A$ 描述位移行为 $r(t)$ 的数学模型（其中包括材料定律；在线性弹性的情况下，$Ar = \mu \Delta r + (\mu+\lambda)\nabla(\nabla\cdot r)$，其中 $\lambda$ 和 $\mu$ 是 Lamé 参数) 到线性方程组 $$ M \ddot{\vec{r} }(t) + K \vec r(t) = M\vec f(t),$$ 其中 $K$ 是刚度矩阵，$M$ 是质量矩阵，正如您在问题中所写。Galerkin 方法），在有限元方法中，它被选为网格每个单元格上的分段多项式的跨度。在这个网站上已经有很多关于它是如何工作的问题的答案（例如，这个），所以我不打算在这里详细介绍；我要指出的是，在分段线性函数的特定情况下，\eqref{1} 解的有限维逼近是由其在网格节点处的值唯一确定的，这正是 $\ vec r$ 表示。

同样，特征模态——它们是静止状态！-- 满足共振方程$$ Ar = \lambda r,$$ 其中$\lambda$ 是特征值。对于时间相关的情况 \eqref{1}，你会得到有限维（广义）特征值问题 $$ K\vec r = \lambda M\vec r.$$ 这现在是一个标准的数值线性代数，存在几种方法（应该指出，这些方法主要通过计算特征向量来工作，从中可以很容易地获得特征值）。特别是，由于 $K$ 和 $M$ 又大又稀疏，所以像 Arnoldi 这样的 Krylov 方法将是最有效的。这正是 COMSOL 所做的（使用实现这些方法的 FORTRAN 库 ARPACK）。

（你说得对，特征向量不是唯一的；通常——尤其是在 Krylov 方法中——你希望将它们归一化为 $r^TMr=1$，正如 Bill Greene 所说，请参阅COMSOL 手册的第 55 页；这只留下符号的选择，它是任意的，通常甚至是随机的（由于方法的初始化）。）