引言
有限元方法需要对求解域整体划分网格，这在处理外部流场、声辐射、大地电磁等无界域问题时格外尴尬：计算域必须人为截断，并辅以吸收边界或无限元近似。边界元方法（Boundary Element Method, BEM）提供了另一条路径——借助微分算子的基本解与 Green 恒等式，把域内偏微分方程精确地转化为边界上的积分方程。求解维度由此降低一维：三维问题只需离散包围它的二维曲面。更重要的是，对于无限域问题，基本解天然满足无穷远处的衰减或辐射条件，无需任何人工截断。本文从 Laplace 方程出发推导边界积分方程，给出常单元离散与内部点重构的完整实现，并讨论 BEM 的优势边界与加速技术。
理论基础
基本解与 Green 第三恒等式
以二维 Laplace 方程 $\nabla^2 u = 0$ 为例，其基本解（fundamental solution）是点源作用下的 Green 函数：
$$G(\mathbf{x}, \boldsymbol{\xi}) = \frac{1}{2\pi} \ln \frac{1}{r}, \qquad r = |\mathbf{x} - \boldsymbol{\xi}|$$
即满足 $\nabla^2 G = -\delta(\mathbf{x} - \boldsymbol{\xi})$。对 $u$ 与 $G$ 在区域 $\Omega$ 上应用 Green 第二恒等式，并令源点 $\boldsymbol{\xi}$ 趋于边界，可得边界积分方程：
$$c(\boldsymbol{\xi}), u(\boldsymbol{\xi}) + \int_\Gamma u(\mathbf{x}) \frac{\partial G}{\partial n}(\mathbf{x}, \boldsymbol{\xi}), d\Gamma = \int_\Gamma G(\mathbf{x}, \boldsymbol{\xi}), q(\mathbf{x}), d\Gamma$$
其中 $q = \partial u / \partial n$ 为法向通量，$c(\boldsymbol{\xi})$ 是边界在源点处的立体角系数——光滑边界上 $c = 1/2$，域内点 $c = 1$。这个方程的深刻之处在于：边界上的 $u$ 与 $q$ 只有一个是独立的，Dirichlet 问题给定 $u$ 求 $q$，Neumann 问题反之，混合问题逐段处理。一旦边界量全部解出，域内任意点的场值可由同一积分式直接重构，且重构误差在远离边界处极小。
离散格式
将边界剖分为 $N$ 个单元，在最简单的常单元（constant element）格式中，$u$ 与 $q$ 在每个单元内取常数。把积分方程配点于各单元中点，得到线性代数方程组：
$$\sum_{j=1}^{N} H_{ij}, u_j = \sum_{j=1}^{N} G_{ij}, q_j, \qquad i = 1, \dots, N$$
即 $H\mathbf{u} = G\mathbf{q}$。其中 $G_{ij} = \int_{\Gamma_j} G, d\Gamma$，$H_{ij} = \int_{\Gamma_j} \partial G / \partial n, d\Gamma + \frac{1}{2}\delta_{ij}$。代入边界条件、交换未知量后，归结为求解一个稠密非对称线性系统——这是 BEM 与稀疏对称的有限元最鲜明的对比。对角元 $G_{ii}$ 含对数奇异性，但对常单元可解析积出：$G_{ii} = \frac{L_i}{2\pi}\left(\ln\frac{2}{L_i} + 1\right)$，$L_i$ 为单元长度。
数值方法与实现
下面实现单位圆内 Laplace 方程 Dirichlet 问题的常单元 BEM：边界条件 $u = \cos\theta$（对应精确解 $u = x$），先求解未知通量 $q$，再重构内部点势值并与精确解比较：
import numpy as np
def bem_circle_dirichlet(N):
"""常单元 BEM：单位圆内 Laplace 问题，边界 u = cos(theta)（精确解 u = x）"""
theta = np.linspace(0, 2 * np.pi, N, endpoint=False)
nodes = np.column_stack([np.cos(theta), np.sin(theta)]) # 单元中点
p1 = np.column_stack([np.cos(theta - np.pi / N), np.sin(theta - np.pi / N)])
p2 = np.column_stack([np.cos(theta + np.pi / N), np.sin(theta + np.pi / N)])
L = np.linalg.norm(p2 - p1, axis=1) # 单元长度
nrm = nodes.copy() # 圆的外法向 = 径向
G = np.zeros((N, N)); H = np.zeros((N, N))
gp = np.array([-0.8611363116, -0.3399810436, 0.3399810436, 0.8611363116])
gw = np.array([0.3478548451, 0.6521451549, 0.6521451549, 0.3478548451])
for i in range(N):
for j in range(N):
if i == j:
G[i, j] = (L[j] / (2 * np.pi)) * (np.log(2 / L[j]) + 1)
H[i, j] = 0.5
else:
for s, w in zip(gp, gw): # 4 点高斯积分
xg = 0.5 * (1 - s) * p1[j] + 0.5 * (1 + s) * p2[j]
r = xg - nodes[i]
dist = np.linalg.norm(r)
G[i, j] += w * (1 / (2 * np.pi)) * np.log(1 / dist) * L[j] / 2
H[i, j] += w * (-1 / (2 * np.pi)) * np.dot(r, nrm[j]) / dist**2 * L[j] / 2
u_bc = np.cos(theta) # Dirichlet 边界条件
q = np.linalg.solve(H, G @ u_bc) # 求非法向通量 q = du/dn
# 内部点势值重构：u(xi) = sum(G*q - H*u)，内部点 c = 1
test_pts = np.array([[0.0, 0.0], [0.5, 0.0], [0.3, 0.4], [-0.6, 0.2]])
errs = []
for xi in test_pts:
u_int = 0.0
for j in range(N):
ge, he = 0.0, 0.0
for s, w in zip(gp, gw):
xg = 0.5 * (1 - s) * p1[j] + 0.5 * (1 + s) * p2[j]
r = xg - xi
dist = np.linalg.norm(r)
ge += w * (1 / (2 * np.pi)) * np.log(1 / dist) * L[j] / 2
he += w * (-1 / (2 * np.pi)) * np.dot(r, nrm[j]) / dist**2 * L[j] / 2
u_int += ge * q[j] - he * u_bc[j]
errs.append(abs(u_int - xi[0]))
print(f" 点 {xi}: BEM = {u_int:.6f}, 精确 = {xi[0]:.6f}")
return max(errs)
for N in [16, 32, 64]:
print(f"N={N}, 内部点最大误差 = {bem_circle_dirichlet(N):.3e}")工程案例分析
BEM 的工程价值在三类问题上尤为突出。其一是无限域与半无限域问题：天线辐射、水波与浮体相互作用、地震波在半空间中的传播等，基本解自动满足 Sommerfeld 辐射条件，截断误差为零。其二是断裂与接触力学：裂纹尖端应力奇异性使得体网格重划分代价高昂，而 BEM 只需在裂纹面上离散，裂纹扩展模拟中仅需局部调整边界网格。其三是形状优化与反问题：设计迭代中几何频繁变化，仅重划表面网格远比重划体网格廉价。
下表总结 BEM 与 FEM 的关键特性对比：
| 特性 | 有限元 FEM | 边界元 BEM |
|---|---|---|
| 离散对象 | 整个求解域（体网格） | 仅边界（面/线网格） |
| 矩阵性质 | 稀疏、对称、带状 | 稠密、非对称 |
| 无限域 | 需人工截断边界 | 基本解自动满足辐射条件 |
| 适用前提 | 任意非线性、非均匀介质 | 需已知基本解（线性问题为主） |
| 应力精度 | 位移一阶、应力低一阶 | 边界量同阶高精度 |
| 规模瓶颈 | 自由度数大但矩阵稀疏 | 矩阵稠密，$O(N^2)$ 存储 |
用本文代码实测：内部点最大误差随单元数 $N = 16, 32, 64$ 分别为 $2.3\times10^{-2}$、$1.0\times10^{-2}$、$5.0\times10^{-3}$，约呈一阶收敛——这正是常单元的理论阶；改用线性或二次单元并采用更高阶几何逼近，可提升至二阶以上。还需指出 BEM 的现实瓶颈：系数矩阵稠密使存储与求解均为 $O(N^2)$ 起步。上世纪 90 年代以来，快速多极子方法（FMM）与 $\mathcal{H}$-矩阵技术把矩阵-向量乘积压缩到 $O(N \log N)$ 甚至 $O(N)$，配合 GMRES 等迭代求解器，使百万自由度级的声学散射 BEM 计算成为常态，商业软件（如声学的 LMS Virtual.Lab/Actran 边界元模块）均内置此类加速。
结论与展望
本文从基本解与 Green 恒等式出发，完整推导了 Laplace 问题的边界积分方程，并给出可运行的常单元 BEM 实现。BEM 以”降维 + 精确满足远场条件”两大优势，在声学、电磁散射、断裂力学与岩土半无限域问题中不可替代；其主要局限在于依赖线性算子的基本解，对强非均匀与非线性问题需借助域积分项或双重互易法（DRM）扩展。当前前沿包括：与等几何分析结合的 IGA-BEM（直接以 CAD 的 NURBS 曲面作为边界单元）、时域边界元，以及面向多物理耦合的快速算法。理解 BEM，也就理解了”格林函数”这一贯穿数学物理的核心概念如何直接转化为工程算法。
发表回复