有限体积法
一、直观理解:每个小格子都做一本守恒账
有限体积法(Finite Volume Method, FVM)从积分守恒律出发:把区域划成控制体(cell),对每个控制体写
也就是:体内累积 = 源 − 净流出。只要相邻控制体共用面上的通量取值一致(左细胞的出就是右细胞的入),离散后就自动保持局部守恒,从而也保持全局守恒。
这与差分“点上逼近导数”、有限元“函数空间投影”视角不同。FVM 更像会计分录:不管你怎么估计界面上的“流速×浓度”,账本在界面两侧必须对得上。
因此,有限体积在计算流体、交通流、浅水波、守恒型对流扩散中非常自然,尤其擅长处理通量与间断(激波)相关问题。
二、适用场景与优缺点
| 类型 | 说明 |
|---|---|
| 适用场景 | 守恒律形式的 PDE;对流占优、可能出现间断;需要严格保证质量/动量/能量离散守恒;非结构网格上的工程流场;交通流、污染物输运、气体动力学简化模型。 |
| 核心优点 | 局部守恒天生满足;物理含义清晰(通量);可在一般多边形/多面体网格上构造;与 Riemann 求解器等工具结合可捕获激波。 |
| 局限与缺点 | 高精度重构与限制器设计复杂;扩散项与复杂边界实现比结构差分繁琐;分析收敛阶不如规则网格差分直观;手写二维以上代码量大。 |
使用前应确认:
- 方程能否写成守恒形式 ;
- 控制体与界面几何量(体积、面积、法向)正确;
- 数值通量满足一致性与(需要时)单调/熵条件;
- 时间步满足 CFL,并检查总质量误差是否在机器精度量级。
三、核心原理:界面数值通量
1. 一维守恒律
单元 ,单元平均 。半离散:
称为数值通量。
2. 常见数值通量
- 迎风(线性对流 ):;
- Lax–Friedrichs:加数值黏性以增强稳定;
- Godunov / Riemann 解:根据左右状态解局部波结构取通量。
一致性要求:。
3. 时间离散
可用显式 Euler / RK:
并受 CFL 限制。对扩散项可隐式处理。
4. 守恒恒等式
对所有单元求和,内部通量抵消,总变化只取决于边界通量与源——这是 FVM 的结构优势。
四、关键决策
-
通量选择
平滑流可用简单迎风或中心+人工黏性;间断问题需要更可靠的 Riemann/限制器策略。 -
CFL
线性对流常取 (格式相关)。超过后轻则振荡重则爆炸。 -
高阶重构
片常数一阶耗散大;线性重构 + minmod 等限制器可提高分辨率并抑制新极值。 -
验证清单
- 常值解保持;
- 总质量与边界通量账一致;
- 光滑区网格收敛;
- 激波位置与速度是否符合 Rankine–Hugoniot。
-
与 FDM/FEM 分工
规则网格扩散问题可用 FDM;复杂几何椭圆问题 FEM;守恒律与通量问题 FVM。
五、MATLAB 示例
%% 1) 一维线性对流:u_t + c u_x = 0,周期边界,迎风有限体积
N = 100;
L = 1;
dx = L/N;
x = ((1:N) - 0.5) * dx; % 单元中心
c = 1;
CFL = 0.8;
dt = CFL * dx / abs(c);
tEnd = 0.5;
nStep = round(tEnd / dt);
dt = tEnd / nStep; % 精确落在 tEnd
CFL = c * dt / dx;
% 初值:方波
u = zeros(1, N);
u(x >= 0.2 & x <= 0.4) = 1;
for n = 1:nStep
% 周期:右界面通量用左侧状态(c>0 迎风)
uL = u;
uR = [u(2:end), u(1)]; % 仅用于中心型;迎风只需 u
flux = c * u; % F_{i+1/2} = c * U_i (c>0)
fluxLeft = [flux(end), flux(1:end-1)];
u = u - dt/dx * (flux - fluxLeft);
end
% 真解:初值平移 c*tEnd
u0 = zeros(1, N);
u0(x >= 0.2 & x <= 0.4) = 1;
% 周期平移
shift = mod(c * tEnd, L);
xBack = mod(x - shift, L);
uTrue = zeros(1, N);
uTrue(xBack >= 0.2 & xBack <= 0.4) = 1;
figure;
stairs(x, u, "LineWidth", 1.2); hold on; grid on;
stairs(x, uTrue, "--");
ylim([-0.2, 1.2]);
legend("迎风 FVM", "平移真解", "Location", "best");
title(sprintf("线性对流:CFL=%.2f,可见数值耗散使方波变钝", CFL));
xlabel("x"); ylabel("u");
fprintf("总质量: num=%.6f, true=%.6f\n", sum(u)*dx, sum(uTrue)*dx);
%% 2) 守恒检查:随机初值 + 周期边界,总质量应几乎不变
u = rand(1, N);
mass0 = sum(u) * dx;
for n = 1:200
flux = c * u;
fluxLeft = [flux(end), flux(1:end-1)];
u = u - dt/dx * (flux - fluxLeft);
end
mass1 = sum(u) * dx;
fprintf("周期对流质量: 初始=%.12e, 最终=%.12e, 差=%.3e\n", ...
mass0, mass1, mass1 - mass0);
%% 3) 带源项的泄漏模型:u_t + c u_x = -lambda u
lambda = 0.5;
u = exp(-20*(x-0.3).^2);
for n = 1:nStep
flux = c * u;
fluxLeft = [flux(end), flux(1:end-1)];
u = u - dt/dx * (flux - fluxLeft) - dt * lambda * u;
end
figure;
plot(x, u, "LineWidth", 1.2); grid on;
title("对流 + 线性衰减源项");
xlabel("x"); ylabel("u");
%% 4) 与中心通量对比(说明为何需要迎风/数值黏性)
% 中心通量在纯对流显式下往往振荡,仅作警示实验
对 Burgers 方程 可用局部 Lax–Friedrichs:
N = 200; L = 1; dx = L/N; x = ((1:N)-0.5)*dx;
u = sin(2*pi*x) + 0.5; % 初值
dt = 0.4 * dx / max(abs(u)); % 粗估
tEnd = 0.2; nStep = ceil(tEnd/dt); dt = tEnd/nStep;
F = @(w) 0.5*w.^2;
for n = 1:nStep
a = max(abs(u)); % 全局波速上界(可改为逐界面)
uR = [u(2:end), u(1)];
uL = u;
flux = 0.5*(F(uL)+F(uR)) - 0.5*a*(uR - uL); % LF
fluxLeft = [flux(end), flux(1:end-1)];
u = u - dt/dx * (flux - fluxLeft);
dt = 0.4 * dx / max(1e-8, max(abs(u))); % 简单自适应
end
figure; plot(x, u, "LineWidth", 1.2); grid on;
title("Burgers:Lax-Friedrichs 有限体积(示意)");
六、如何解释结果
有限体积报告应强调守恒与通量,而不仅是曲线形状:
- 守恒律形式与源项含义;
- 网格(单元数、)、时间步与 CFL;
- 数值通量类型与边界处理;
- 总质量/总动量随时间的变化表;
- 与精确解或参考解的误差,以及数值耗散/色散表现;
- 间断位置、高度是否物理合理;
- 若做工程结论,说明网格收敛后关键通量是否稳定。
方波被迎风“抹圆”是数值黏性,不是物理扩散,除非模型本身含扩散项——解释时必须分开。
七、常见误区
1. 界面通量左右定义不一致
同一界面两侧用了不同 ,守恒立即破坏。这是实现中最致命的 bug 之一。
2. 把数值耗散当成物理扩散去标定系数
迎风抹平宽度随网格变化;加密后变窄。用其拟合真实扩散会得到网格依赖的假参数。
3. 忽视 CFL
“空间很细”反而迫使时间步更小;只加密空间不改 会不稳定。
4. 高阶重构不加限制器
在间断附近产生过冲/新极值,浓度变负。
5. 非守恒形式离散却宣称有限体积
对非线性通量,非守恒型差分可能导致错误激波速度。
八、小结
有限体积法的核心思想可以概括为:
在控制体上积分守恒律,把问题转化为界面数值通量的计算,从而在离散层面保持局部守恒。
它是守恒律与通量型问题的自然选择。建模中使用 FVM 时,优先保证通量一致性、CFL 与全局守恒检查,再追求高阶与复杂重构;并清楚区分数值黏性与物理扩散,以免误解释机制。