有限单元方法（Finite Element Method, FEM）把复杂计算区域划分为有限个互不重叠、通过节点连接的单元，在每个单元内用低阶基函数近似未知场，再把所有单元方程组装成总体代数方程。与规则网格上的有限差分方法相比，有限元对复杂几何、非均匀网格和自然边界条件具有更强的适应性。

有限单元法的基本思想

有限元计算通常包括以下步骤：

1. 把偏微分方程写成适合离散的弱形式；
2. 将计算区域剖分为一维线段、二维三角形或四边形等单元；
3. 在每个单元中选取节点和形函数；
4. 用节点未知量的线性组合近似单元内部的场；
5. 计算单元矩阵和单元载荷向量；
6. 按局部节点与总体节点的对应关系进行组装；
7. 施加边界条件并求解总体线性方程组；
8. 根据节点解计算单元内部的场、梯度、应力等派生量。

有限元近似并不要求控制方程在每个点严格成立，而是要求残差对选定的权函数在积分意义下为零。

加权残差法

考虑边值问题

$$Lu=f\quad\text{in }\Omega,$$

以及边界条件。用有限项展开

$$u_h(\mathbf x)=\sum_{j=1}^{N}U_jN_j(\mathbf x)$$

近似精确解，其中 $N_j$ 是基函数或形函数，$U_j$ 是待求节点值。代入原方程产生残差

$$R=Lu_h-f.$$

加权残差法要求

$$\int_\Omega w_iR\,d\Omega=0, \qquad i=1,2,\ldots,N.$$

若取 $w_i=N_i$，便得到 Galerkin 有限元方法。

强形式与弱形式

一维 Poisson 问题

以一维方程为例：

$$-\frac{d}{dx}\left(k\frac{du}{dx}\right)=f, \qquad x\in(0,L).$$

在边界 $\Gamma_D$ 上给定 $u=\bar u$，在边界 $\Gamma_N$ 上给定通量

$$k\frac{du}{dx}n=\bar q.$$

强形式要求 $u$ 具有足够高的光滑性，使二阶导数在区域内有定义。将方程乘以检验函数 $w$ 并积分：

$$-\int_0^Lw\frac{d}{dx} \left(k\frac{du}{dx}\right)dx =\int_0^Lwf\,dx.$$

分部积分得

$$\int_0^Lk\frac{dw}{dx}\frac{du}{dx}\,dx -\left[wk\frac{du}{dx}\right]_0^L =\int_0^Lwf\,dx.$$

结合 Neumann 边界条件，可以写成弱形式：寻找满足 Dirichlet 边界条件的 $u$，使对所有在 $\Gamma_D$ 上为零的检验函数 $w$，都有

$$\int_\Omega k\nabla w\cdot\nabla u\,d\Omega =\int_\Omega wf\,d\Omega +\int_{\Gamma_N}w\bar q\,d\Gamma.$$

分部积分把 $u$ 的最高导数由二阶降为一阶，因此弱形式对近似函数光滑性的要求更低，也使分片多项式基函数成为可能。

本质边界与自然边界

• Dirichlet 边界条件直接规定未知量 $u$，称为本质边界条件，需要在试探空间或总体方程中显式施加；
• Neumann 边界条件规定通量或牵引，经过分部积分自然出现在边界积分中，称为自然边界条件。

这一区别是有限元边界处理的核心。

单元、节点与形函数

将区域划分为 $N_e$ 个单元：

$$\Omega=\bigcup_{e=1}^{N_e}\Omega^e, \qquad \Omega^e\cap\Omega^r=\varnothing\quad(e\ne r)$$

其中相邻单元只在公共边界或节点处相交。

在单元 $e$ 内，近似解写为

$$u_h^e(\mathbf x) =\sum_{a=1}^{n_e}N_a^e(\mathbf x)U_a^e =\mathbf N^e\mathbf U^e.$$

形函数应满足节点插值性质

$$N_a^e(\mathbf x_b)=\delta_{ab},$$

以及单位分解

$$\sum_{a=1}^{n_e}N_a^e(\mathbf x)=1.$$

单位分解保证常数场能够被精确表示。相邻单元的公共节点共享同一个总体未知量，从而保证 $u_h$ 在单元之间连续。

一维线性单元

物理坐标中的形函数

设单元 $e$ 的两个节点为 $x_1^e,x_2^e$，长度

$$h_e=x_2^e-x_1^e.$$

一维线性形函数为

$$N_1(x)=\frac{x_2^e-x}{h_e}, \qquad N_2(x)=\frac{x-x_1^e}{h_e}.$$

单元近似解为

$$u_h^e(x)=N_1(x)U_1^e+N_2(x)U_2^e.$$

其导数为常数：

$$\frac{du_h^e}{dx} =\begin{bmatrix}-1/h_e&1/h_e\end{bmatrix} \begin{bmatrix}U_1^e\\U_2^e\end{bmatrix}.$$

单元刚度矩阵

对于 Poisson 型方程，单元刚度矩阵为

$$K_{ab}^e =\int_{x_1^e}^{x_2^e} k\frac{dN_a}{dx}\frac{dN_b}{dx}\,dx.$$

若单元内 $k$ 为常数，则

$$\mathbf K^e =\frac{k}{h_e} \begin{bmatrix} 1&-1\\ -1&1 \end{bmatrix}.$$

单元载荷向量为

$$F_a^e=\int_{x_1^e}^{x_2^e}N_a f\,dx.$$

若 $f$ 在单元内为常数，则

$$\mathbf F^e =\frac{fh_e}{2} \begin{bmatrix}1\\1\end{bmatrix}.$$

一维二次单元

三节点二次单元在自然坐标 $\xi\in[-1,1]$ 上取节点 $-1,0,1$，形函数为

$$N_1(\xi)=\frac12\xi(\xi-1),$$

$$N_2(\xi)=1-\xi^2,$$

$$N_3(\xi)=\frac12\xi(\xi+1).$$

二次单元可以精确表示二次变化的场，通常比相同单元数的线性单元精度更高，但每个单元自由度更多，数值积分和组装也更复杂。

自然坐标与等参变换

在标准单元中定义形函数，再把标准单元映射到物理单元：

$$x(\xi)=\sum_{a=1}^{n_e}N_a(\xi)x_a^e.$$

若几何映射和场变量插值使用同一组形函数，称为等参单元。Jacobian 为

$$J=\frac{dx}{d\xi}.$$

积分和导数按

$$dx=J\,d\xi,$$

$$\frac{dN_a}{dx}=\frac{dN_a}{d\xi}\frac{d\xi}{dx} =J^{-1}\frac{dN_a}{d\xi}$$

进行转换。二维问题中 Jacobian 是矩阵，其行列式控制面积变换；若 $\det\mathbf J\leq0$，说明单元节点顺序错误或单元发生翻转。

单元组装

设单元 $e$ 的局部节点编号对应总体节点

$$\mathcal I^e=(I_1^e,I_2^e,\ldots,I_{n_e}^e).$$

单元矩阵元素按映射累加到总体矩阵：

$$K_{I_a^eI_b^e}\mathrel{+}=K_{ab}^e,$$

$$F_{I_a^e}\mathrel{+}=F_a^e.$$

组装后得到

$$\mathbf K\mathbf U=\mathbf F.$$

每个单元只连接少数相邻节点，因此总体矩阵通常是大型稀疏矩阵。有限元程序不应按一般稠密矩阵存储和求解，而应使用稀疏矩阵格式和相应线性求解器。

边界条件的施加

Dirichlet 条件

若总体自由度 $r$ 满足 $U_r=\bar U_r$，常见处理方式包括：

1. 消元：把已知值对应的项移到右端，并删除相应行列；
2. 行列修正：把第 $r$ 行和列清零、对角元置 1，右端置为 $\bar U_r$，同时修正其他方程右端；
3. 罚函数法：在对应对角元加入很大的罚参数。

直接消元或一致的行列修正能够保留矩阵的对称性；只修改一行而不处理相应列可能破坏对称结构。

Neumann 条件

通量或牵引通过边界积分进入载荷向量：

$$F_a^{\Gamma} =\int_{\Gamma_N}N_a\bar q\,d\Gamma.$$

若给定零通量，则该边界积分为零，不需要额外修改总体矩阵。

数值积分

单元矩阵通常需要在标准单元上积分。Gauss 求积写为

$$\int_{-1}^{1}g(\xi)\,d\xi \approx\sum_{q=1}^{n_q}w_qg(\xi_q).$$

$n_q$ 点 Gauss-Legendre 求积可以精确积分不超过 $2n_q-1$ 次的多项式。二维四边形单元常使用一维 Gauss 点的张量积；三角形单元使用相应的三角形积分公式。

积分点过少会产生欠积分，可能引起零能模式；积分点过多虽然准确，但增加计算成本。应根据形函数阶次、材料参数变化和积分项次数选择求积阶数。

二维三角形线性单元

三节点三角形单元内采用线性插值：

$$u_h(x,y)=N_1U_1+N_2U_2+N_3U_3.$$

面积坐标形函数可以写为

$$N_i(x,y)=\frac{a_i+b_ix+c_iy}{2A},$$

其中 $A$ 是三角形面积，系数由三个节点坐标确定。形函数梯度在单元内为常数：

$$\nabla N_i=\frac{1}{2A} \begin{bmatrix}b_i\\c_i\end{bmatrix}.$$

对标量 Poisson 方程，若 $k$ 在单元内为常数，则

$$K_{ij}^e =\int_{\Omega^e}k\nabla N_i\cdot\nabla N_j\,d\Omega =\frac{k}{4A}(b_ib_j+c_ic_j).$$

线性三角形能够灵活逼近复杂边界，但单元内梯度为常数。要提高精度，可以加密网格或使用高阶三角形单元。

二维四边形单元

四节点双线性单元在自然坐标 $(\xi,\eta)\in[-1,1]^2$ 上的形函数为

$$N_1=\frac14(1-\xi)(1-\eta),$$

$$N_2=\frac14(1+\xi)(1-\eta),$$

$$N_3=\frac14(1+\xi)(1+\eta),$$

$$N_4=\frac14(1-\xi)(1+\eta).$$

几何映射为

$$x(\xi,\eta)=\sum_{a=1}^{4}N_ax_a, \qquad y(\xi,\eta)=\sum_{a=1}^{4}N_ay_a.$$

Jacobian 矩阵为

$$\mathbf J= \begin{bmatrix} \partial x/\partial\xi&\partial y/\partial\xi\\ \partial x/\partial\eta&\partial y/\partial\eta \end{bmatrix}.$$

形函数对物理坐标的梯度由

$$\begin{bmatrix} \partial N_a/\partial x\\ \partial N_a/\partial y \end{bmatrix} =\mathbf J^{-1} \begin{bmatrix} \partial N_a/\partial\xi\\ \partial N_a/\partial\eta \end{bmatrix}$$

求得，面积元为

$$d\Omega=|\det\mathbf J|\,d\xi d\eta.$$

四边形单元适合结构化或半结构化网格，但严重扭曲会降低插值精度和矩阵条件。

含时问题的有限元离散

以波动方程

$$\rho\frac{\partial^2u}{\partial t^2} -\nabla\cdot(c\nabla u)=s$$

为例，对空间作 Galerkin 离散得到

$$\mathbf M\ddot{\mathbf U} +\mathbf K\mathbf U =\mathbf F(t),$$

其中质量矩阵和刚度矩阵分别为

$$M_{ij}=\int_\Omega\rho N_iN_j\,d\Omega,$$

$$K_{ij}=\int_\Omega c\nabla N_i\cdot\nabla N_j\,d\Omega.$$

若考虑阻尼，则得到

$$\mathbf M\ddot{\mathbf U} +\mathbf C\dot{\mathbf U} +\mathbf K\mathbf U =\mathbf F(t).$$

一致质量矩阵与集中质量矩阵

按形函数积分得到的 $\mathbf M$ 称为一致质量矩阵，一般不是对角矩阵。将每行质量适当集中到对角线上，可得到集中质量矩阵（lumped mass matrix）。对角质量矩阵使显式时间推进无需每步求解线性方程组，因此在大规模波动模拟中非常重要。

采用中心差分时，

$$\ddot{\mathbf U}^n \approx\frac{\mathbf U^{n+1}-2\mathbf U^n+\mathbf U^{n-1}}{\Delta t^2},$$

从而

$$\mathbf U^{n+1} =2\mathbf U^n-\mathbf U^{n-1} +\Delta t^2\mathbf M^{-1} \left(\mathbf F^n-\mathbf K\mathbf U^n\right).$$

若 $\mathbf M$ 已集中为对角矩阵，$\mathbf M^{-1}$ 只对应逐节点除法。

稳定时间步长受离散系统最高固有频率限制：

$$\Delta t\leq\frac{2}{\omega_{\max}},$$

其中 $\omega_{\max}^2$ 是广义特征值问题

$$\mathbf K\boldsymbol\phi =\omega^2\mathbf M\boldsymbol\phi$$

的最大特征值。

一维算例：Poisson 方程

考虑

$$-u''(x)=f(x),\qquad x\in(0,1),$$

$$u(0)=u(1)=0.$$

使用均匀线性单元，单元刚度矩阵为

$$\mathbf K^e=\frac1h \begin{bmatrix}1&-1\\-1&1\end{bmatrix}.$$

组装后内部节点方程为

$$\frac1h \left(-U_{i-1}+2U_i-U_{i+1}\right)=F_i.$$

它在形式上与有限差分离散非常相似，但右端 $F_i$ 是源函数与形函数的积分，而不是简单的点值。两种方法在均匀一维线性问题中可能得到相同或相近的代数结构，但它们的出发点不同：有限差分直接近似微分算子，有限元则离散弱形式。

有限元与有限差分的比较

特征	有限差分法	有限单元法
基本出发点	用差商近似导数	对弱形式作分片基函数近似
网格	规则网格最方便	非结构网格和复杂几何更灵活
边界	复杂边界处理较困难	自然适应曲折边界和 Neumann 条件
系数矩阵	规则、模板固定	稀疏但需单元积分和组装
局部加密	实现相对困难	可通过单元尺寸和阶次灵活控制
编程工作	简单问题较直接	网格、映射、积分和组装更复杂
波动模拟	显式格式高效	质量集中后也可高效显式推进

有限元并非在所有情况下都优于有限差分。规则矩形区域中的大规模波传播问题，有限差分往往更简单高效；复杂地形、曲线界面或局部网格加密问题，有限元通常更有优势。

版权所有

版权归属：Staaaaaaaaar

许可证：署名-非商业性-相同方式共享 4.0 国际 (CC-BY-NC-SA-4.0)

编辑此页

贡献者: Staaaaaaaaar

上一页二维地球介质中地震波传播计算