物理定律、偏微分方程和数值建模   有限元法 

有限元法简介

空间和时间相关问题的物理定律通常用偏微分方程（PDE）来描述。对于绝大多数的几何结构和所面对的问题来说，可能无法求出这些偏微分方程的解析解。不过，在通常的情况下，可以根据不同的离散化 类型来构造出近似的方程，得出与这些偏微分方程近似的数值模型方程，并可以用数值方法求解。如此，这些数值模型方程的解就是相应的偏微分方程真实解的近似解。有限元法（FEM）就是用来计算出这些近似解的。

举例来说，某函数 u 可能是一个偏微分方程中的因变量（即温度、电势、压力等）。可以根据下列表达式，通过基函数的线性组合将函数 u 近似为新的函数 u_h：

⁽¹⁾

以及

⁽²⁾

在此，ψ_i 代表这些基函数，而 u_i 则代表用来对 u 进行近似的 u_h 函数中的系数。下图用一个一维问题阐明这一原理。例如，u 可以表示某一均匀加热的杆在特定长度（x）处的温度。此图中的线性基函数的值，在各自的节点处为 1，在其他节点处为 0。在这个例子中，函数 u 的定义域所在的 x-轴部分（即这根杆的长度）共有七个单元。

使用有限元法的好处之一就是该方法在离散度的选择方面提供了极大的自由（同时包括用于离散空间和离散基函数的单元的离散度选择）。比如，在上图中，这些单元均匀地分布在 x-轴上（虽然并不总会是这种情况）。在函数 u 的一个梯度较大的区域中，也可以使用较小的单元，如下所示。

这两幅图都表明，选定的线性基函数在 x-轴方向上获得的支持（仅有一个较为狭窄的非零区间）和重叠非常有限。根据手头的问题，可以选择其他类型的函数而不是线性函数。

有限元法的另一个优点是该理论已经发展得较为成熟了，原因就在于偏微分方程问题的数值表述式和弱表达式之间的密切关系（见下面的部分）。例如，当数值模型方程在计算机上求解时，该理论在误差估计或误差边界 估计方面是较为有效的。

回顾有限元方法的历史，可知该方法是在 20 世纪 40 年代初被德裔美国数学家 Richard Courant 首次提出的。虽然 Courant 报道了有限元方法在诸多问题上的应用，但几十年之后该方法才在结构力学之外的领域获得了普遍的应用，成就了现在的地位。

代数方程、常微分方程、偏微分方程和物理定律

物理定律通常使用数学语言来表达。例如，各类守恒定律（如能量守恒定律、质量守恒定律和动量守恒定律等）都可以用偏微分方程（PDE）来表达。这些定律也可以用相关变量（包括温度、密度、速度、电势以及其他因变量）的本构关系来表达。

微分方程包含有相应的表达式，可以在自变量（x, y, z, t）发生变化时确定因变量的小幅变化。这一小幅变化也被称为因变量对应于自变量的导数。

假设一个固体具有时变温度，但在空间上的变化忽略不计。在这种情况下，通过内能（热）守恒方程，就可以导出在热源 g 的作用下，随着时间的小幅变化而发生的温度变化的方程式：

⁽³⁾

在此， 表示密度，而 C_p 则代表热容量。温度 T 是因变量，时间 t 是自变量。函数 可以描述随温度和时间而变化的一个热源。方程 (3) 表明，如果温度在随着时间而变化，则它必然会由热源 所平衡（或所引起）。此方程是用一个自变量（t）的导数所表示的一个微分方程。这种微分方程被称为常微分方程（ODE）。

在某些情况下，当某一时间的温度 t₀ 为已知时（称为初始条件），即可得到方程 (3) 的一个解析解，表达式如下：

⁽⁴⁾

如此，该固体中的温度通过一个代数方程（4）来表示，其中的某个时间值 t₁ 就会有一个对应时间的温度值 T₁。

物理属性常常会随着时间和空间发生变化。例如，该固体中靠近热源处的温度可能比其他位置略高。这种温度差异进而引起该固体内部不同部分之间产生热通量。在这种情况下，根据能量守恒定律就可以导出一个传热方程，该方程同时具有时间变量和空间变量（x），如：

⁽⁵⁾

同之前一样，T 是因变量，而 x（x = (x, y, z)）和 t 则是自变量。该固体中的热通量矢量由 q = (q_x, _y, q_z) 表示，而 q 的发散 则描述了热通量沿着空间坐标的变化。在笛卡尔坐标系中，q 的发散被定义为：

⁽⁶⁾

因此，方程（5）表明，在所有方向上都有了改变时，如果净通量发生了变化，以至于 q 的发散（变化的总和）不为零，则必须有一个热源以及/或者随时间变化的温度变化来进行平衡（或引发）。

可以通过传导热通量的本构关系来描述一个固体中的热通量，这也称为傅里叶定律：

⁽⁷⁾

在上述方程式中，k 表示导热系数。方程（7）表明，在导热系数为比例常数的情况下，热通量与温度梯度成正比。方程（7）（(5) 中）给出了以下的微分方程：

⁽⁸⁾

在此，导数是以 t、x、y 和 z 表示的。在某个微分方程是用一个以上的自变量的导数来表示的情况下， 该微分方程就被称为偏微分方程（PDE），这是因为每个导数都可能代表（几个可能方向中的）某个方向上的变化。还需注意的是，常微分方程中的导数是用 d 来表示的，而偏微分方程中的导数则是用更卷曲的 ∂ 来表示的。

除了方程（8），还可以知道的就是某个时间 t₀ 上的温度或者某个位置 x₀ 上的热通量。此类知识可应用于方程（8）的初始条件和边界条件。在许多情况下，偏微分方程都无法通过解析方法来求解（即得出不同时间和位置下的因变量的值）。有时，要得到一个如下的解析表达式，可能非常困难，甚至几乎是不可能的，例如方程（8）中的：

⁽⁹⁾

在不用解析法求解偏微分方程的前提下，另一种方案就是通过寻找近似的数值解 来求解数值模型方程。有限元法正是这种类型的方法——一种求解偏微分方程的数值方法。

类似于上面提到的热能守恒方程，可以推导出动量守恒与质量守恒的方程（这两个方程构成了流体动力学的基础）。此外，亦可以推导出空变与时变问题中的电磁场和通量方程，从而得到偏微分方程组。

继续这一讨论，让我们看看如何从偏微分方程中推导出所谓的弱形式公式。

源自于弱公式的有限元法：基函数和试函数

假定正在研究的一个散热器中的温度分布由方程（8）给出，但现正处于稳定状态，这就意味着方程（8）中的温度场的时间导数为零。模型域 Ω 的域方程如下：

⁽¹⁰⁾

此外，假定沿边界（∂Ω₁）的温度已知，同时垂直于其他一些边界（∂Ω₂）的热通量的表达式也已知。在其余的边界上，热通量在向外的方向（∂Ω₃）上为零。这些边界上的边界条件就成为：

⁽¹¹⁾

⁽¹²⁾

⁽¹³⁾

其中，h 表示传热系数，T_amb 表示环境温度。边界表面上向外的单位法向矢量由 n 表示。 方程（10）至（13）描述了这一散热器的数学模型，如下图所示。

下一步是将方程（10）的两边都乘以一个试函数 φ，并在域 Ω 上积分：

⁽¹⁴⁾

试函数 φ 与方程的解 T 被假定属于希尔伯特空间（Hilbert space）。希尔伯特空间是一个具有无限维度的函数空间，并带有具备特定属性的函数。它可以被看作是具有一定属性的函数的集合；这样一来，这些函数可以同向量空间中的普通向量一样被方便地操作。例如，可以在该集合中生成函数的线性组合（这些函数有明确的长度，称为模 ），并且可以像欧几里德矢量一样测量函数之间的角度。

实际上，可以通过有限元方法简单地将这些函数转换为普通的矢量。有限元法是一种系统性的方法，将无限维函数空间中的函数转换为有限维函数空间中的一类函数，最后再转换为可以用数值方法处理的普通矢量（在某一矢量空间中）。

如果要求（14）对试函数空间中的所有试函数都成立，而不是方程（10）对 Ω 中的所有点都成立，则可以得到弱形式公式。因此，基于方程（10）的问题公式有时也称为逐点公式。在我们所说的伽辽金法 中，假设解 T 同测试函数属于相同的希尔伯特空间。这通常写为 φ ϵ h 和 T ϵ h，其中 H 表示希尔伯特空间。 使用格林第一恒等式（实质上是进行分部积分）, 就可以推出以下方程（14）：

⁽¹⁵⁾

通过要求此等式对希尔伯特空间中的所有 试函数都成立，可以实现方程（10）的弱形式公式或称为变分公式。之所以说是“弱”，是因为其放宽了（10）的要求，也就是偏微分方程的各项在每一个点上都必须被明确定义的要求。相反的是，只有在积分时才要求（14）和（15）是相等的关系。例如，弱公式化完全允许解的一阶导数不连续，因为这种情况并不妨碍积分。但是，它为二阶导数引入的分布 则并不是普通意义上的函数。因此，在不连续点上要求（10）成立是没有意义的。

有时可以对某个分布进行积分，以使（14）被明确定义。可以证明的是，弱公式化以及通过（13）得到的边界条件（11）都是与通过逐点公式化求出的解直接相关的。此外，对于解可微分 的情况（即二阶导数明确定义），这些解是相同的。这些公式化是等效的，因为从（10）推导（15）的过程依赖于格林第一恒等式，而其只有在 T 有连续的二阶导数的情况下才成立。

这是有限元公式化的第一步。利用弱公式化，就有可能对数学模型方程进行离散化，从而得到数值模型方程。可以利用伽辽金法——许多可能的有限元法公式化中的一种——来进行离散化。

首先，要实现离散化，就意味着要在希尔伯特空间 H 的有限维子空间中寻找方程（15）的近似解；如此，T ≈ T_h。这就是说，近似解被表示为一组属于子空间的基函数 ψ_i 的线性组合：

⁽¹⁶⁾

由此，对每个试函数 ψ_j 而言，方程（15）的离散化形式即变为：

⁽¹⁷⁾

这里的未知数，就是函数 T(x) 的近似解中的系数 T_i。随后，方程（17）就变成了一个方程组，该方程组与有限维函数空间拥有相同的维度。如果使用了试函数 ψ_j 中的数字 n，使 j 从 1 一直变到 n，那么就可以根据（17）得到一个方程数量为 n 的方程组。方程（16）中也有 n 个未知的系数（T_i）。

一旦体系被离散化并被施加了边界条件后， 根据以下表达式就可以得到一个方程组：

⁽¹⁸⁾

其中，T 是未知矢量，且 T _h = {T₁, .., T_i, …, T_n}；A 则是一个 nxn 的矩阵，其元素 A_ji 中的每个方程 j 都含有系数 T_i。右边是维度从 1 到 n 的矢量。A 是系统矩阵，通常称为（消除）的刚度矩阵 ——这是有限元方法的首次应用，也是其在结构力学中的用途。

如果源函数在温度方面是非线性的，或者传热系数取决于温度，那么该方程组也是非线性的，矢量 b 就成为了未知系数 T_i 的一个非线性函数。

有限元方法的优点之一是它能够选择试函数和基函数。在非常小的几何区域的支集之上，是有可能选择试函数和基函数的。这意味着，方程（17）在任意一处都为零——除非是在函数 ψ_j 和 ψ_i 重叠的非常有限的区域上，因为上面所有的积分都包括了函数 i 和 j（或它们的梯度）的乘积。很难用三维空间来描述试函数和基函数的支集，但其二维的类比却是能够被可视化的。

假设有一个二维的几何域，并且选用了 x 和 y 的线性函数，每个函数在点 i 上的值为 1，但在其他点 k 上的值为零。下一步是使用三角形对这一二维域进行离散化，并为某一三角形网格中的两个相邻节点 i 和 j 给出基函数（试函数或形函数）。

两个相邻的基函数共享两个三角形的单元。因此，两个基函数之间有一些重叠，如上所示。此外，请注意，如果 i = j，则函数之间会完全重叠。这些贡献形成了未知矢量 T 的系数，这一未知矢量与系统矩阵的对角线分量 A_jj 相对应。

比如说，假设现在这两个基函数更进一步地分开了。这两个函数不共享单元，但它们有一个共同的单元顶点。如下图所示，它们不重叠。

当这两个基函数重叠时，方程（17）具有非零值，且对系统矩阵的贡献也是非零的。当没有重叠时，积分为零，因此对系统矩阵的贡献也为零。

这意味着，在从 1 到 n 的节点上，对（17）的方程组中的每个方程来说，它们都只能从共享同一个单元的相邻节点中得到若干个非零的项。方程（18）中的系统矩阵 A 变得稀疏，而对应于重叠 ij:s 的矩阵分量才有非零项。这一代数方程组的解可以作为该偏微分方程的近似解。网格越稠密，近似解就越接近真实解。

瞬态问题（时变问题）

可以在瞬态（时变）的情况下进一步定义该散热器中的热能平衡。根据伽辽金方法，每个试函数 ψ_j 的离散弱公式化可以写作：

⁽¹⁹⁾

在此，系数 T_i 是时变函数，而基函数和试函数则仅依赖于空间坐标。再者，在时间域上的时间导数不是离散的。

一种方法是对时间域也使用有限元法，但这种做法可能会耗费大量的计算资源。经常采取的另一种方案则是通过直线法来对时间域进行独立的离散化。比如可以使用有限差分法。其最简单的形式可以用下面的差分近似法来表示：

⁽²⁰⁾

给出的是方程（19）中的两个可能有限差分逼近。当未知的系数 T_i,t 以 t + Δt 的形式表示时，就可以得到第一个式子：

⁽²¹⁾

在面对线性问题时，在每一个时间步长上都需要求解一个线性方程组。如果是非线性的问题，则必须在每个时间步长内求解相应的非线性方程组。由于在 t + Δt 处的解是被方程（21）隐含地给出的，所以这种时间推进方案被称为隐式法。

第二个公式则基于 t 处的解：

⁽²²⁾

该式表明，一旦在某一给定时间上的解（T_i,t）已知，那么方程（22）就能显式地给出在 t + Δt (T_{i, t+Δt}) 处的解。换言之，对于一个显式的时间推进方案，不需要在每个时间步长上都求解一个方程组。显式时间推进方案的缺点是它们有一个稳定性方面的时间步进限制。对于热问题来说（如此处所强调的情况），显式方法需要非常短的时间步长。隐式方案允许更大的时间步长，减少了如（22）这样的方程所需的计算资源（在每一个时间步长上都要对这些方程进行求解）。

在实践中，现代化的时间步进算法会根据具体问题自动在显式和隐式步进法之间切换。此外，方程（20）中的差分方程被替换为一个多项式，其阶次和步长可以发生变化，具体取决于所要解决的问题和求解所需的时间。现代化的时间推进方案会根据数值解的时间演化来自动地控制多项式的阶次以及步长。

下面有几个例子，对最常用的几种方法加以说明：

• 向后微分公式（BDF）法
• 广义 α 法
• 不同的 Runge-Kutta 法

不同的单元

如上所述，伽辽金法采用了与基函数和试函数相同的函数集。然而，即使是这种方法，也可以通过很多种方式（理论上是无穷多的）来定义基函数（即伽辽金有限元公式中的单元）。让我们来回顾一下最常用的几种单元。

对于二维和三维的线性函数，最常见的元素如下图所示。此图和上图 给出的是线性基函数（被定义在三角形网格中，形成了三角形的线性单元）。基函数被表示为节点位置（二维时：x 和 y；三维时：x、y 和 z）的函数。

在二维面上，矩形单元常常被用于结构力学分析。它们还可用于计算流体动力学（CFD）和传热建模中的边界层网格剖分。它们的三维类比就是所谓的六面体单元，后者也常被应用于结构力学和边界层网格剖分。在从六面体边界层单元到四面体单元的过渡中，锥体单元通常被放置在边界层单元的顶端。

下图显示的是相应的二阶单元（二次单元）。在此，面对一个域边界的边和面通常是弯曲的，而面对该域内部的边和面则是直线或平面。但是请注意，也可以将所有的边和曲都定义为是弯曲的。拉格朗日单元和巧凑边点元是二维和三维建模中最常用的单元类型。拉格朗日单元使用下面所有的节点（黑色、白色和灰色），而巧凑边点元则不使用灰色的节点。

博客“在多物理场模型中追踪单元阶次”中给出了二阶（二次）拉格朗日元的二维图形，非常漂亮。在上述单元的内部，很难用三维的形式描述这些二次基函数的基，但是可以用色块来表示单元表面的函数数值。

在讨论有限元法时，需要考虑的一个重要因素就是误差估计。原因在于，当达到估计出的误差宽容度时，就会发生收敛。请注意，这里的讨论具有更一般的性质，而不是局限于特定的有限元方法。

有限元法给出的是数学模型方程的一个近似解。数值方程的解与数学模型方程的精确解之间的差值就是误差：e = u - u_h。

在许多情况下，可以在得出数值方程的解之前就估计出误差的大小（即先验 误差估计）。先验 估计通常仅用于预测所用有限元方法的收敛阶数。例如，如果某个问题是适定的，并且相应的数值方法可以收敛，那么根据 O(h^α)（其中 α 表示收敛阶数），随着通常的单元尺寸 h 的减小，误差的模也会减小。由此可见，随着网格密度的增加，误差的模也会快速地降低。

不过，只有简单的问题才能进行先验 估计。此外，估计出的结果往往会包含不同的未知常数，从而不可能给出定量的预测。后验 估计使用的则是近似解，并结合了相关问题的其他近似，以估计出误差的模。

构造解方法

一个非常简单但却通用的误差估计方法（用于数值方法和偏微分问题），就是对问题进行略微改动——如这一篇博客文章 所述—— 使预定义的解析表达式成为改动后的问题的真实解。这种方法的优点是未对数值方法或其背后的数学问题进行过假设。此外，由于解是已知的，所以可以很容易地计算出误差的大小。通过谨慎地选择分析表达式，就可以对方法和问题的不同方面进行研究。

让我们来看一个例子，对这一点进行说明。假设有一种数值方法可以对一个单位正方形（Ω）上的泊松方程进行求解，且该正方形具有齐次边界条件

⁽²³⁾

⁽²⁴⁾

此方法可用于对改动后的问题进行求解

⁽²⁵⁾

⁽²⁶⁾

其中，

⁽²⁷⁾

这里，

是可以被自由选择的一个解析表达式。另外，如果

⁽²⁸⁾

则 是改动后的问题的精确解，此时可以直接计算出误差大小：

⁽²⁹⁾

如此，就可以为不同选择的离散化程度和 计算出误差及其模。如果改动后问题的解与未改动问题的解具有相同的特性，那么改动后问题的误差就可以用作未改动问题的近似误差。在实践中，可能很难知道是否是这种情况——这是此方法的缺点。这种方法的优点在于它的简单性和普遍性：既可以用于非线性问题和时变问题（瞬态问题），也可以用于任何的数值方法。

目标定向的误差估计

如果可以从近似解中选择出一个函数（或泛函数），并将其作为一个重点物理量来进行误差估计，那么就可以通过解析方法精准地估算出此物理量的计算误差（或界限）。此类估计依赖于对偏微分方程残差的后验 计算，以及对所谓的对偶问题 进行的近似求解。对偶问题与所选择的函数是直接相关的（并由其定义）。

这种方法的缺点在于其依赖于“对偶问题”的精确计算，而且只给出了所选函数的误差估计（而没有涉及其他物理量）。这种方法的优势在于其较高的普适性和较合理的资源消耗（用于误差计算）。

网格收敛

网格收敛是一种简单的方法，该方法比较了不同的网格剖分方案所得到的近似解。在理想情况下，一套非常精细的网格剖分方案所得出的近似解就可以作为真实解的近似了。较粗化的网格剖分方案所得出的近似解的误差，可以由下式直接估算出：

⁽³⁰⁾

在实践中，要对非常精细的网格剖分方案（比实际所需精细得多的方案）进行近似求解，其实是较困难的。因此，在习惯上会使用最精细网格的近似来达到此目的。对每个网格细化来说，也可以从所得解的变化情况中估算出收敛性。如果近似解位于一个收敛的区域之中，那么这个解的变化幅度会随着网格的细化而收窄；如此，所得的近似解也就会越来越接近于真实解。

下图显示了一个椭圆形膜的结构力学基准模型；得益于对称性，只需要对该膜的四分之一部分进行计算就可以了。载荷作用在该几何的外边缘上，而沿 x 和 y 轴的边界被认为是对称的。

对不同网格类型和单元尺寸的数值模型方程进行求解。例如，下图描述了用于二次基函数的矩形拉格朗日单元，这些基函数是根据上图得出的。

根据更早给出的这幅图，对该点上的应力和应变进行了计算。下面的图表显示的是此点上的 σ_x 所得的相对值。此值应为零，因此与零值的任何差异都是一种误差。为了得到一个相对误差，将计算出的 σ_x 除以计算出的 σ_y，以便为相对误差的估计给出正确的数量级。

上图表明，随着各单元的单元尺寸（h）的减小，相对误差也相应减小。在这种情况下，随着基函数阶次（单元阶次）的升高，收敛曲线也变得更为陡峭。不过，需要注意的是，在单元尺寸一定的情况下，随着阶次的升高，数值模型中的未知项的数量也会增加。这就意味着，当我们增加单元的阶次时，我们也要为更高的精确度而付出代价，这种代价就是计算耗时的增加。如果不使用更高阶的单元，还可以采取的另一种方法就是为较低阶次的单元选用更细化的网格。

网格自适应

在计算出了这些数值方程的解 u_h 之后，就可以用后验 局部误差估计值来创建一个密度更大的网格，该网格具有较大的误差。然后可以使用细化的网格来计算出第二个近似解。

下图描述了一个被加热的圆柱体在受到稳态流体流动作用下的温度场。对这一稳态问题进行了两次求解：一次是用基础网格，另一次是用一个细化网格（被基本网格计算出的误差估计所控制）。该细化网格在温度和热通量方面的计算精度更高，而这一点可能正是该实例所需要的。

对时变（瞬态）的对流问题来说，也可以通过前序时步的解来实现对流网格的细化。在下图给出的例子中，相场被用来计算喷墨打印机中的墨水液滴与空气之间的界面。该界面是由相场函数的等值面所给出的，其值等于 0.5。在这个界面上，相场函数的值迅速地从 1 变为 0。在此相场函数的这些陡峭梯度的周围，我们可以使用误差估计来自动完成网格细化的工作，而流场则可以用来对流网格细化，以便仅在相场等值面的面前才使用更细的网格。

其他有限元公式

在上述例子中，我们为基函数和试函数使用了相同的函数集来实现模型方程的离散化。如果一个有限元公式可以使试函数不同于基函数，则该公式称为 Petrov-Galerkin 法。这是一种常用的方法；例如，在解决对流-扩散问题的过程中，只会对流线方向进行稳定化处理。其也被称为流线迎风 /Petrov-Galerkin（SUPG）法。

在耦合方程组的求解过程中，不同的因变量可能会用到不同的基函数。一个典型的例子是纳维-斯托克斯方程的求解，其中的压力往往比速度更平滑、更易进行近似。在某类方法中，如果一个耦合方程组中不同的因变量的基函数（以及试函数）属于不同的函数空间，那么这类方法便称为混合有限元法。

发布日期：2016 年 3 月 15 日
上次修改日期：2017 年 2 月 21 日