有限单元法

有限单元法（英语：Finite element method），即使用有限单元分析物理现象，是一种用于求解微分方程组或积分方程组数值解的数值方法。

在解偏微分方程的过程中，主要难点是如何构造一个方程来逼近原本研究的方程，并且该过程还需要保持数值稳定性。目前有许多处理的方法，他们各有利弊。当区域改变时（就像一个边界可变的固体），当需要的精确度在整个区域上变化，或者当解缺少光滑性时，有限元方法是在复杂区域（像汽车、船体结构、输油管道）上解偏微分方程的一个很好的选择。

为了解决问题，有限单元法将大型物理系统细分为更小、更简单的部分，称为有限元（英文：finite element）。这是通过在空间维度上进行特定的空间离散化来实现的，该离散化是通过构建对象的网格实现的：解决方案的数值域具有有限数量的点。边值问题的有限单元法公式化最终形成了一个代数方程组。该方法在域上近似未知函数^[1]。然后，将对这些有限元建模的简单方程式组合成一个对整个问题进行建模的较大方程式系统。然后，有限单元法通过最小化关联的误差函数，使用来自变异演算的变异方法来近似求解。

由研究人员建立的有限单元网格，然后使用套装软件找到解决磁性问题的方法。颜色表示分析人员已为每个区域设置了材料属性，在此图中，浅蓝色的铁磁成分（可能是铁），浅紫色的是空气。

有限单元法解决左侧问题的方法，涉及电磁屏蔽。磁性圆筒形遮蔽罩保护内部区域不受外部磁场影响。如插入图例中的比例尺所示，颜色表示磁场的强度，红色为高强度磁场。圆柱体内部的区域是低强度的（深蓝色，具有宽间隔的磁通量线），这表明遮蔽罩的性能达到了设计目标。

将整个物理系统细分为更简单的部分具有以下优点^[2]：

精确表示复杂的几何形状。
可以描述多样的材料特性。
轻松表示整体解决方案。
精确描述局部现象。

该方法的工作流程包括

（1）将问题的域划分为子域的集合，每个子域由一组单元方程表示为原始问题，然后（2）系统地将所有单元方程组重组为用于最终计算的全域方程组。

在上面的第一步中，单元方程是简化过的方程，可以局部地近似要研究的原始复杂方程组，其中原始方程通常是偏微分方程。为了求此方程式的近似解，通常将有限单元法作为伽辽金法的特例来处理。用数学语言来说，该过程是将残差和加权函数取内积，并将该积分设为零。简而言之，它是通过将试验函数拟合到偏微分方程中来最小化近似误差的过程。残差是由试验函数引起的误差，权重函数是投影残差的多项式逼近函数。该过程消除了偏微分方程中的所有空间导数，从而使偏微分方程局部近似为一组稳态问题的代数方程，或是一组用于瞬态问题的常微分方程。如果基础偏微分方程是线性的，则单元方程也是线性的，反之亦然。稳态问题中出现的代数方程组，便利用数值线性代数方法求解，而瞬态问题中出现的常微分方程组则使用其他数值方法（例如欧拉方法或Runge-Kutta法）通过数值积分来求解。

历史

有限元法最初起源于土木工程和航空工程中的弹性和结构分析问题的研究。它的发展可以追溯到Alexander Hrennikoff（1941）和Richard Courant (1942)的工作。这些先驱者使用的方法具有很大的差异，但是他们具有共同的本质特征：利用网格离散化将一个连续区域转化为一族离散的子区域，通常叫做元.Hrennikoff的工作离散用类似于格子的网格离散区域; Courant的方法将区域分解为有限个三角形的子区域，用于求解来源于圆柱体转矩问题的二阶椭圆偏微分方程. Courant的贡献推动了有限元的发展，绘制了早期偏微分方程的研究结果。

有限元方法的发展开始于五十年代中后期使用在机身框架和结构分析上，并于六十年代通过斯图加特大学的John Argyris（英语：John Argyris）和柏克莱加州大学的Ray W. Clough（英语：Ray W. Clough）在土木工程中的应用工作中积累经验。

基于五十年代至六十年代大型水坝计算研究的实践经验，1965年，中国计算数学专家冯康发表了《基于变分原理的差分格式》一文，奠定了有限元计算方法的严格数学理论，为后世有限元计算方法的实际应用提供了理论保证。且冯康教授的“有限元法”严密理论体系是先于西方的，是国际公认的当代计算数学的一项重大成就，不同的是冯康教授只是从数学方面提出有限元法的。^[3]^[4]

有限元概念

单元

单元（Element）是由节点组成的几何体，如三角形单元，四面体单元等。

节点

节点（Node）是单元几何体的端点、顶点或特定点，单元的各物理量变化均体现在节点上，例如在弹性力学问题中，一个有两个节点的线单元的质量集中在两个节点上，受力也只能作用在节点上，变形也用节点的位移表示。

自由度

节点自由度（Degree of Freedom，简写 DoF），是节点上变量的个数，例如用位移法解结构问题时节点自由度为3，表示单个节点上三个坐标方向上的位移，又例如热分析时节点自由度为1，表示某个节点处的温度值。

网格

网格（Mesh）是由多个单元通过共用节点组成的单元网络，用以表示待解问题域。

分析方法

以下用有限元分析解决两个简单问题，更一般的问题可以类似的推导出来。

P1是一个较简单的一维问题

{\mbox{ P1 }}:{\begin{cases}u''(x)=f(x){\mbox{ in }}(0,1),\\u(0)=u(1)=0,\end{cases}}

其中 $f$ 是已知函数, $u$ 是关于 $x$ 的未知函数, $u''$ 是 $u$ 对 $x$ 的二阶导数。

二维比较简单的问题是狄利克雷问题

{\mbox{P2 }}:{\begin{cases}u_{xx}(x,y)+u_{yy}(x,y)=f(x,y)&{\mbox{ in }}\Omega ,\\u=0&{\mbox{ on }}\partial \Omega ,\end{cases}}

其中 $\Omega$ 是 $(x,y)$ 平面上的连通开区域，它的边界 $\partial \Omega$ 是良好的（例如，光滑流形或多边形）, $u_{xx}$ 和 $u_{yy}$ 分别表示 $x$ 和 $y$ 的二阶导数。问题P1能够通过计算不定积分而直接解决。然而，解决边值问题的这一方法只有在空间维数为1时才可用，并且不能推广到高维问题以及形如 $u+u''=f$ 的问题。出于这种考虑，我们将用有限元方法解决P1并将其推广至问题P2.

我们的描述分为两步，每步都反映了用有限元解决边值问题的本质。

将原问题描述为它的弱形式，或变分形式。这一步很少或不需要计算。
离散化，将弱形式在有限维空间离散化。

这两步之后，我们可以构造一个大型有限维线性方程，线性方程的解就是原边值问题的逼近解。然后，这一有限维问题由计算机求解。

弱解形式

P1的弱解形式

第一步是将问题P1和P2转化为他的等价变分形式，或弱解形式。

如果 $u$ 是问题P1的解，那么对任何满足边界条件的光滑函数 $v$ ,有

（1） $\int _{0}^{1}f(x)v(x)\,\mathrm {d} x=\int _{0}^{1}u''(x)v(x)\,\mathrm {d} x$

相反如果 $u$ 对任何光滑函数 $v(x)$ 满足 $u(0)=u(1)=0$ 和（1），

那么 $u$ 是P1的解。对于二次可导函数 $u$ 证明这一点是非常容易的（利用中值定理）。

通过对（1）的右侧使用分部积分，可以得到

（2） ${\begin{aligned}\int _{0}^{1}f(x)v(x)\,\mathrm {d} x&=\int _{0}^{1}u''(x)v(x)\,\mathrm {d} x\\&=u'(x)v(x)|_{0}^{1}-\int _{0}^{1}u'(x)v'(x)\,\mathrm {d} x\\&=-\int _{0}^{1}u'(x)v'(x)\,\mathrm {d} x=-\phi (u,v)\end{aligned}}$

其中假设 $v(0)=v(1)=0$ 。

P2的弱解形式

f当我们使用格林恒等式来表示式（2）， P2可以 $u$ 的积分型式表示，在此定义 $\phi (u,v)$

\int _{\Omega }fv\,ds=-\int _{\Omega }\nabla u\cdot \nabla v\,ds\equiv -\phi (u,v),

此处 $\nabla$ 代表梯度，即为二维平面上的内积。另外 $\,\!\phi$ 可以转为内积空间 $H_{0}^{1}(\Omega )$ ，且 $\Omega$ 的一次微分函数 $\partial \Omega$ 为零。我们也可以假设 $v\in H_{0}^{1}(\Omega )$ （详见索伯列夫空间）也可以显示解的存在性和唯一性。

证明解的存在性和唯一性

离散化

A function in

H_{0}^{1},

with zero values at the endpoints (blue), and a piecewise linear approximation (red)

P1 和 P2 通过上述过程被离散化，并简化为子问题 (3)。基本思路是将无限维线性问题替换掉：

找到

u\in H_{0}^{1}

使

\forall v\in H_{0}^{1},\;-\phi (u,v)=\int fv

表示唯有限维度的形式：

子问题(3) Find

u\in V

such that

\forall v\in V,\;-\phi (u,v)=\int fv

此处 $V$ 是 $H_{0}^{1}$ 的一个有限维度线性子空间。

$V$ 有许多可能的形式，但对于有限单元而言，在此将假定 $V$ 存在于分段多项式函数的空间中。

对于 P1

在此在区间 $(0,1)$ 之中选择 $n$ 个 $x$ 的可能值 $0=x_{0}<x_{1}<\cdots <x_{n}<x_{n+1}=1$ 接着定义 $V$ 为：

V=\{v:[0,1]\rightarrow \mathbb {R} \;:v{\mbox{ is continuous, }}v|_{[x_{k},x_{k+1}]}{\mbox{ is linear for }}k=0,\dots ,n{\mbox{, and }}v(0)=v(1)=0\}

令e $x_{0}=0$ 和 $x_{n+1}=1$ 。观察到在 $V$ 之中的函数根据微积分的基本定义是不可微分的。当然，当 $v\in V$ ，则通常不定义 $x=x_{k}$ ， $k=1,\ldots ,n$ 的导数，但导数事实上存在于每一个 $x$ 的位置，并可以利用这些导数来进行部分积分运算。

两种不同维度（二维及三维空间）的分段线性函数

对于 P2

$V$ 是属于 $\Omega$ 的一系列函数。在右图中，图片下半部是一个15边形的平面 $\Omega$ 的三角分割（英语：Polygon triangulation），以及该多边形的分段线性函数（图片上半部彩色部分），即在三角分割（英语：Polygon triangulation）所形成的每个三角形上呈线性；空间 $V$ 则由在所在的三角分割（英语：Polygon triangulation）的每个三角形上的函数线性组合而成。

我们希望当下面的三角形网格变得越来越精细，离散子问题（3）的解在某种意义上将收敛到原始边界值问题P2的解。为了测量此网格的细度，三角分割（英语：Polygon triangulation）由一很小的实数 $h>0$ 所表示。此参数将与三角分割（英语：Polygon triangulation）中最大或平均三角形的大小有关。当我们提高三角分割（英语：Polygon triangulation）的精度时（分割出更多三角形），分段线性函数的空间 $V$ 应会随 $h$ 变动。因此，在某些文献中会以 $V_{h}$ 来代表。

参考文献

^ Logan, Daryl L. A first course in the finite element method 5th ed. Stamford, CT: Cengage Learning. 2012. ISBN 978-0-495-66825-1. OCLC 664675951.
^ Reddy, J. N., An Introduction to the Finite Element Method (Third ed.). An Introduction to the Finite Element Method (Third ed.). McGraw-Hill. 2006. ISBN 9780071267618. 请检查|isbn=值 (帮助).
^ 冯, 康. 基于变分原理的差分格式. 应用数学与计算数学. 1965.
^ 实践出题、直觉判断、求异思维——冯康的创新要诀. 中国科学院数学与系统科学研究院. [2010-07-26]. （原始内容存档于2013-07-10）.

外部链接

MIT Video Lecture on the Finite Element Method
Multiphysics Glossary （页面存档备份，存于互联网档案馆）（Glossary of Multiphysics and Finite Element Modeling terms by COMSOL）
NAFEMS（页面存档备份，存于互联网档案馆） -- The International Association for the Engineering Analysis Community
IFER （页面存档备份，存于互联网档案馆） -- Internet Finite Element Resources - an annotated list of FEA links and programs
Workshop "The Finite Element Method in Biomedical Engineering, Biomechanics and Related Fields" （页面存档备份，存于互联网档案馆）
Finite Element Analysis Resources- Finite Element news, articles and tips
COMSOL Multiphysics Finite Element Analysis Software （页面存档备份，存于互联网档案馆） - Official site
CAD, Finite Element Analysis（Abaqus,Ansys）, CAE, Programming （页面存档备份，存于互联网档案馆）- FEM, CAD, Programming, discussion forums
Finite Element Books- books bibliography
Mathematics of the Finite Element Method （页面存档备份，存于互联网档案馆）
Finite Element Methods for Partial Differential Equations （页面存档备份，存于互联网档案馆）
FEM AVI-gallery at CompMechLab site, St.Petersburg State Polytechnical University, Russia
Intro to FEA
Introduction to FEA for EM modeling（includes list of currently available software）
Finite Element modeling of light feapower^{[永久失效链接]}
propagation （页面存档备份，存于互联网档案馆）
World Association of Fatigue, Durability and Fracture Mechanics - Fatigue for Finite Element Models （页面存档备份，存于互联网档案馆）

[1] Logan, Daryl L. A first course in the finite element method 5th ed. Stamford, CT: Cengage Learning. 2012. ISBN 978-0-495-66825-1. OCLC 664675951.

[2] Reddy, J. N., An Introduction to the Finite Element Method (Third ed.). An Introduction to the Finite Element Method (Third ed.). McGraw-Hill. 2006. ISBN 9780071267618. 请检查|isbn=值 (帮助).

[FKOrig-3] 冯, 康. 基于变分原理的差分格式. 应用数学与计算数学. 1965.

[4] 实践出题、直觉判断、求异思维——冯康的创新要诀. 中国科学院数学与系统科学研究院. [2010-07-26]. （原始内容存档于2013-07-10）.

[1]

[2]

[3]

[4]