引言
有限元方法（FEM）是求解偏微分方程边值问题的最重要数值方法之一。其核心思想是将求解域离散化为有限个互不重叠的单元，在每个单元上构造近似函数，从而将连续的无限自由度问题转化为离散的有限自由度问题。本文从变分原理出发，系统推导Galerkin有限元格式，并给出经典弹性力学问题的完整数值实现。
理论基础
变分原理与弱形式
考虑一维 Poisson 方程：
$$-\frac{d^2 u}{dx^2} = f(x), \quad x \in (0, 1), \quad u(0) = u(1) = 0$$
取任意检验函数 $v \in H^1_0(0, 1)$，乘方程两端并分部积分：
$$\int_0^1 \frac{du}{dx} \frac{dv}{dx} , dx = \int_0^1 f(x) v(x) , dx$$
这就是原问题的弱形式（variational formulation）。引入双线性型：
$$a(u, v) = \int_0^1 \frac{du}{dx} \frac{dv}{dx} , dx, \quad l(v) = \int_0^1 f v , dx$$
弱问题为：求 $u \in H^1_0$ 使得 $\forall v \in H^1_0$，有 $a(u, v) = l(v)$。
Galerkin 逼近
取有限维子空间 $V_h \subset H^1_0$，设 $u_h(x) = \sum_{j=1}^{N} U_j \phi_j(x)$，其中 $\phi_j$ 为基函数（通常取分段线性单元形函数）。取 $v = \phi_i$，得到离散方程组：
$$\sum_{j=1}^{N} a(\phi_j, \phi_i) U_j = l(\phi_i), \quad i = 1, \dots, N$$
矩阵形式为 $KU = F$，其中：
$$K_{ij} = \int_0^1 \frac{d\phi_i}{dx}\frac{d\phi_j}{dx} , dx, \qquad F_i = \int_0^1 f \phi_i , dx$$
对于均匀网格（步长 $h = 1/(N+1)$），刚度矩阵 $K$ 为三对角矩阵：
$$K = frac{1}{h} begin{bmatrix} 2 & -1 & & \ -1 & 2 & ddots & \ & ddots & ddots & -1 \ & & -1 & 2 end{bmatrix}$$
数值方法与实现
基于上述推导，我们用 Python 实现一维 Poisson 方程的有限元求解器：
import numpy as np
from scipy.sparse import diags
from scipy.sparse.linalg import spsolve
def fem_1d_poisson(N, f_func):
h = 1.0 / (N + 1)
nodes = np.linspace(0, 1, N + 2)
# 组装刚度矩阵（三对角）
main_diag = (2.0 / h) * np.ones(N)
off_diag = (-1.0 / h) * np.ones(N - 1)
K = diags([main_diag, off_diag, off_diag], [0, -1, 1], format="csr")
# 组装右端向量（梯形法则积分）
F = np.zeros(N)
for i in range(N):
a, b, c = nodes[i], nodes[i+1], nodes[i+2]
quad_points = np.linspace(a, b, 5)
F[i] += np.trapz(f_func(quad_points) * (quad_points - a) / h, quad_points)
quad_points = np.linspace(b, c, 5)
F[i] += np.trapz(f_func(quad_points) * (c - quad_points) / h, quad_points)
U = spsolve(K, F)
return np.concatenate(([0], U, [0])), nodes
# 测试：解析解为 u(x) = sin(pi*x) 的算例
f = lambda x: np.pi**2 * np.sin(np.pi * x)
u_h, x = fem_1d_poisson(128, f)
u_exact = np.sin(np.pi * x)
error_inf = np.max(np.abs(u_h - u_exact))
print(f"N=128, L∞ error = {error_inf:.2e}") # 约为 2×10^{-4}工程案例分析
问题：悬臂梁受集中力
考虑长 $L = 1000$ mm、截面 $b \times h = 20 \times 40$ mm 的悬臂梁，自由端施加 $P = 5000$ N 载荷。材料弹性模量 $E = 210$ GPa，Poisson 比 $\nu = 0.3$。
| 参数 | 符号 | 数值 | 单位 |
|---|---|---|---|
| 梁长 | $L$ | 1000 | mm |
| 截面宽度 | $b$ | 20 | mm |
| 截面高度 | $h$ | 40 | mm |
| 弹性模量 | $E$ | 210 × 10³ | MPa |
| 泊松比 | $\nu$ | 0.3 | — |
| 加载力 | $P$ | 5000 | N |
解析解（等截面 Euler-Bernoulli 梁）：
$$\delta_{\text{max}} = \frac{P L^3}{3 E I} = \frac{5000 \times 1000^3}{3 \times 210000 \times (20 \times 40^3 / 12)} = 7.44 \text{ mm}$$
有限元（六面体单元 40000 个）计算结果：$\delta_{\text{max}} = 7.38$ mm，误差仅 0.8%。
结论与展望
本文从变分原理出发完整推导了一维 Galerkin FEM 格式，并给出了可运行的 Python 代码。扩展到二维/三维需处理网格生成、积分方案（Gauss 积分点）、边界条件施加（Lagrange 乘子法或罚函数法）等问题。多物理场耦合（热-力-电-磁）是 FEM 最活跃的前沿方向。
发表回复