有限元法 (FEM)

有限元法 (Finite Element Method)

有限元法（Finite Element Method, FEM）是一种求解偏微分方程（Partial Differential Equations, PDEs）边值问题和初边值问题的数值方法。其核心思想是将连续求解域离散为有限个互不重叠的简单子区域（称为"单元"），在每个单元上构造近似函数，再通过变分原理或加权残量法将原PDE转化为代数方程组进行求解。有限元法因其对复杂几何形状和非齐次材料的强适应性，成为工程分析和科学计算中应用最广泛的数值方法之一，覆盖固体力学、流体力学、热传导、电磁场和声学等领域。

有限元法的发展可追溯至20世纪40年代。1943年，Richard Courant在研究扭转问题的变分方法时首次提出了在三角形子区域上使用分片连续函数近似解的思想，被视为有限元法的雏形。1956年，M. J. Turner和R. W. Clough等在飞机结构分析中首次系统性地使用了三角形和矩形单元，Clough于1960年正式提出了"有限单元法"这一术语。20世纪60至70年代，O. C. Zienkiewicz和J. T. Oden等为有限元法建立了严格的数学基础，并将其从固体力学推广到其他物理领域。此后，随着计算机技术的迅猛发展，FEM成为现代工程设计不可或缺的工具。

基本原理

有限元法的数学基础建立在变分形式（也称弱形式）之上。考虑一个典型的椭圆型边值问题：

其中$\Omega\subset\mathbb{R}^d$为有界区域，$a(\mathbf{x})$和$b(\mathbf{x})$为已知系数，$f(\mathbf{x})$为源项。对方程两边乘以检验函数$v$并在$\Omega$上积分，利用格林公式降低导数阶数，得到弱形式：

其中$H_0^1(\Omega)$为满足边界条件的Sobolev空间。这一变分表述将原始方程的要求从"每一点满足"弱化为"在加权积分意义下满足"，为使用分片低阶多项式逼近奠定了基础。

有限元法的核心步骤分为三步。

第一，离散化：将区域$\Omega$划分为有限个简单单元。一维问题使用线段单元，二维问题使用三角形或四边形单元，三维问题使用四面体或六面体单元。相邻单元在共有的节点处连接，所有单元的集合构成求解域的一个网格。单元尺寸越小，数值解越精确，但计算代价也越大。

第二，单元分析：在每个单元上构造形函数（shape functions）。设单元$e$上有$n_e$个节点，节点值为$u_i$，则在单元内构造插值近似：

其中$N_i^{(e)}(\mathbf{x})$为形函数，具有"紧支性"——仅在当前单元内非零，并在节点$i$处取值为1、在其他节点处取值为0。形函数的构造方式决定了单元的精度和收敛性质。最常用的是Lagrange型形函数，对应线性或高次多项式插值；此外还有Hermite型形函数，在节点处同时匹配函数值和导数值。

将单元近似代入弱形式，得到单元级的线性代数方程（称为单元刚度方程）：

其中$\mathbf{K}^{(e)}$为单元刚度矩阵，其元素为$K_{ij}^{(e)}=\int_{e}(a\nabla N_i\cdot\nabla N_j+bN_iN_j)d\mathbf{x}$；$\mathbf{f}^{(e)}$为单元载荷向量。

第三，总体组装与求解：将各单元的刚度矩阵按节点的全局编号"组装"成总体刚度矩阵$\mathbf{K}$——这是一个大规模的稀疏矩阵，因为每个单元只与其相邻的少数单元通过共享节点产生耦合。组装的本质是将每个单元的局部自由度贡献叠加到全局自由度上，对应数学上的"直接求和"操作。引入边界条件后，得到总体线性系统$\mathbf{K}\mathbf{u}=\mathbf{f}$。由于$\mathbf{K}$的对称正定性（对自伴椭圆问题），通常采用共轭梯度法、Cholesky分解或多重网格法高效求解。

收敛性与误差估计

有限元解的误差来自两个方面：一是离散误差（因有限维逼近代替连续问题所致），二是数值积分误差（因使用数值积分近似单元矩阵积分所致）。对线性单元，误差的先验估计为：

其中$h$为单元尺寸（最大单元直径），$C$为与网格形状正则性有关的常数。该估计表明：当网格加密（$h\to0$）时，解在$H^1$范数下线性收敛；若使用$p$次多项式单元，收敛阶为$O(h^p)$。在$L^2$范数下，通过Aubin-Nitsche对偶论证可得到更高一阶的收敛速率$O(h^{p+1})$。

后验误差估计通过计算当前解与某种增强逼近之间的局部差异来估计误差，是实现自适应网格加密（adaptive mesh refinement）的理论基础。常用的策略是：在误差较大的区域局部加密网格，而在解足够光滑的区域保持粗网格，从而在给定自由度总数下获得最优的计算精度。

重要扩展

混合有限元法：对涉及多个物理场的耦合问题（如达西渗流中的流速和压力、斯托克斯流中的速度和压力），将不同的物理量用不同的形函数空间逼近，可避免应力/通量的后处理误差并满足inf-sup条件（Ladyzhenskaya-Babuška-Brezzi条件）。

非协调有限元：放松单元间解的连续性要求，允许单元交界面上的位移具有微小跳跃。Crouzeix-Raviart单元是典型的非协调线性三角形单元，在处理不可压缩约束时特别有效。

等参单元：通过相同的形函数同时描述几何坐标和物理场值，使单元边界的曲率可以精确匹配曲线/曲面边界。等参映射将物理空间中的一般形状单元映射到参数空间中的标准单元（如单位正方形或标准三角形），大大提高了处理复杂几何的能力。

hp-有限元法：同时加密网格（$h$-加密）和提升单元多项式阶次（$p$-加密）。在解具有局部奇异性的区域使用$h$-加密，在光滑区域使用高次单元，可实现指数收敛，是求解高精度问题（如断裂力学中的应力强度因子计算）的有力工具。

扩展有限元法（XFEM）：通过在形函数中嵌入富集函数（如阶跃函数和裂尖渐近函数），允许裂纹等强不连续面独立于网格存在，避免了传统FEM中需沿裂纹面重新划分离散网格的繁琐操作。

时域有限元：对时间域的离散，可将有限元法扩展至抛物型和双曲型问题。典型策略包括空间使用有限元离散、时间使用有限差分（如Newmark-β法）或Rothe方法（先时后空），以及时空全有限元法。

计算实现与软件

有限元法的实际计算步骤通常为：前处理（几何建模、网格划分、材料属性定义）→单元刚度计算与组装→施加边界条件和载荷→求解稀疏线性系统→后处理（可视化、误差分析、导出应力/应变/通量等导数量）。主流商业软件包括ABAQUS、ANSYS、COMSOL Multiphysics和NASTRAN；开源替代方案有deal.II、FEniCS、LibMesh和FreeFEM。这些软件通常提供丰富的单元库、材料模型和多物理场耦合接口。

在大规模并行计算中，FEM面临的关键挑战包括：非结构化网格的区域分解与负载均衡（使用ParMETIS或Zoltan进行图划分）、高效通信下的共轭梯度法并行化，以及可扩展的预条件子（如域分解预条件和代数多重网格）。当代FEM计算已可处理包含数十亿自由度的大规模工程问题。

经济学与社会科学中的应用

虽然FEM起源于工程领域，但其方法论已被引入经济学研究。在空间经济学中，FEM被用于求解连续空间上的人口迁移模型和区域经济增长模型，其中的扩散-聚集方程本质上是反应-扩散方程（属于PDE）。在宏观经济学中，异质性代理人模型有时需求解关于财富或收入分布的函数方程，FEM提供的分片多项式逼近是替代有限差分法和投影法的有效方案。在环境经济学中，FEM被用于求解地下水污染扩散方程和大气污染物传输的数值模拟。此外，最优控制问题中的Hamilton-Jacobi-Bellman方程也可通过FEM进行数值逼近，从而求解连续时间的动态优化问题。

局限性与前沿方向

传统FEM的主要限制包括：对于对流占优问题（如高Péclet数流动）的数值不稳定，需引入迎风格式或SUPG（流线迎风Petrov-Galerkin）稳定化；处理高速冲击和接触碰撞等强非线性问题时，网格可能严重扭曲导致求解失败，此时无网格法（如SPH）或物质点法（MPM）更具优势；对于高频波动问题（如声学散射），每个波长内需要足够多的单元，计算代价随频率三次方增长，广义有限元法（GFEM）和超网格（megagrid）方法正在改善这一瓶颈。

当前FEM的前沿方向包括：等几何分析（使用CAD中的NURBS基函数直接进行有限元分析，消除网格生成瓶颈）、增减量有限元（基于形函数导数计算极值量）、多尺度有限元法（MsFEM，在粗网格上构造捕捉细观尺度信息的基函数）、不确定性量化（将随机场参数化输入FEM并计算输出的概率分布）以及基于深度学习的FEM加速（用神经网络作为PDE求解器或替代传统预条件子）。