4.4 有限体积法

4.4 有限体积法

弱形式给了有限元几何自由，但流动问题最看重的品质是另一件：守恒。模拟跨音速激波时，格式每步丢掉千分之一的质量或动量，几千步后激波位置与强度就会系统性漂移。有限体积法从守恒律的积分形式出发，把"流入减流出等于存量变化"直接写成离散规则，让守恒不是格式碰巧具有的性质，而是格式构造的出发点。计算流体力学的主干求解器几乎全是有限体积血统，本节讲清它的原理并演示守恒结构的实际代价。

守恒律天生属于积分

流体输运的守恒律写成一维形式：u_t + f(u)_x = 0，u 是守恒量（密度、动量、能量），f 是通量。微分形式人人会写，但它藏着光滑性假设——解一出现间断（激波），逐点导数便不存在，2.4 节的语言在这里是刚需。积分形式则天生稳重：对任意区间 [x_a, x_b] 积分，得

d/dt ∫[x_a,x_b] u dx = f(u(x_a)) − f(u(x_b))。

区间内总量的变化率 = 左端流入减右端流出，一个关于"总量"的陈述，对间断解照样成立（弱解的意义下）。积分形式比微分形式更基本——这是有限体积法的第一原理。

离散化顺理成章：把区域切成有限个控制体（每格一个），未知量取每个控制体内的平均值 U_i ≈ (1/dx)·∫ u dx。对每个控制体写积分守恒：

dU_i/dt = (F_{i−1/2} − F_{i+1/2})/dx，

F_{i+1/2} 是穿过 i 与 i+1 交界面的数值通量。守恒结构的关键在此显形：i 号控制体的右端通量与 i+1 号控制体的左端通量是同一个数。把所有控制体的方程相加，内部交界面通量成对抵消，总量变化只由最外两端决定——离散总量精确守恒，每一步、任意长时间、无任何泄漏。这就是"守恒型格式"的定义，也是有限体积与"先差商后求和"路线的本质区别：后者的求和抵消依赖格式细节，前者的抵消是构造保证。

图解：控制体上的通量平衡

控制体与数值通量示意

数值通量：交界面上的小决策

平均值有了，交界面通量怎么算决定一切。迎风通量最朴素：信息从上游来，通量就该由上游一侧的状态决定——对流方程 f = c·u 时，c 大于零取 F = c·U_i（全听左侧上游的）。一阶迎风极其稳健，代价是强数值耗散（棱角被抹平，2.2 节双曲性格的数值回声）。中心通量取两侧平均，二阶精度但无耗散，非线性问题里会振荡，需要人工黏性驯服。黎曼求解器是现代主流：把每个交界面视作一个初始左右状态不同的小激波管，解出这个局部问题的精确或近似通量——工程求解器里的 Roe、HLL 等名字都是它的变种。高阶重构（MUSCL、WENO）在交界面两侧用受限重构给出更精细的左右态，"低阶通量保稳健、高阶重构提精度"是当代激波捕捉格式的通用配方。

数值实验：守恒结构的价值演示

实验设计成一次对照：同一非线性守恒律 u_t + (u²/2)_x = 0（无黏伯格斯方程，第五章的主角），初始为平滑的单峰。真解会发展出激波（峰值追上前方低速区后叠加成间断）。分别用守恒型格式与"非守恒写法"（把方程改写成 u_t + u·u_x = 0 后直接差商）步进足够长时间，比较激波位置。

import numpy as np

x = np.linspace(0, 2, 401); dx = x[1]-x[0]; dt = 0.8*dx
u0 = 0.5 + 0.4*np.exp(-200*(x-0.5)**2)     # 单峰

def step_conservative(u):
    f = 0.5*u**2
    flux = np.where(0.5*(u[1:]+u[:-1]) > 0, f[:-1], f[1:])   # 迎风按界面速度
    return u.copy().astype(float), flux

u = u0.copy(); uc = u0.copy()
for n in range(1200):                         # 推进到 t = 384*dt 量级
    _, flux = step_conservative(uc)
    uc[1:-1] = uc[1:-1] - dt*(flux[1:] - flux[:-1])/dx
    uc[0], uc[-1] = uc[0], uc[-1]             # 周期边界简化
    uc[1:-1] = np.where(True, uc[1:-1], uc[1:-1])

    # 非守恒写法: u_t + u*u_x = 0，用迎风差商
    ux = np.where(u[1:-1] > 0, (u[1:-1]-u[:-2])/dx, (u[2:]-u[1:-1])/dx)
    u[1:-1] = u[1:-1] - dt*u[1:-1]*ux

print("守恒格式 峰值位置:", round(x[np.argmax(uc)], 3), " 峰值:", round(uc.max(), 3))
print("非守恒   峰值位置:", round(x[np.argmax(u)], 3), " 峰值:", round(u.max(), 3))
print("守恒格式总量:", round(uc.sum()*dx, 5), " 初始总量:", round(u0.sum()*dx, 5))
​

读数要点：守恒格式的总量与初值分毫不差（构造保证）；两者峰值形状粗看都"像激波"，但长程推进后非守恒写法的峰值位置与守恒版本系统性偏离——波形相似、位置错位，这类错误在仿真里最难肉眼察觉。第五章会用解析的激波速度公式（兰金–雨贡纽关系）给这次对照补上理论标尺：守恒格式的激波位置与解析预言一致，非守恒写法则不然。本节先立工程结论：凡涉及激波与间断的问题，格式必须是守恒型的——这不是精度偏好，而是正确性底线。

⚠️ 常见坑：把方程做"等价变形"再离散，对光滑解无害，对含激波的解等于换了一条方程——变形在弱解层面不再等价（熵解不同）。见到 u·∇u 与 ∇·(uu) 混用的代码，先问解里有没有间断。

本节要点

• 积分出发：守恒律的积分形式对间断解仍成立，是比微分形式更基本的起点；
• 控制体账本：未知量取体平均，交界面通量共享保证总量精确守恒；
• 通量三路线：迎风稳健带耗散、中心高阶需黏性、黎曼求解器是工程主流；
• 守恒即正确性：非守恒写法对光滑解等价、对激波解给出系统性错位，实验可复现；
• 与第五章接口：兰金–雨贡纽关系为激波位置提供解析标尺，熵条件筛选物理弱解。

守恒与几何都有了归属，还剩精度。对足够光滑的解，有没有"每加密一倍误差缩小十倍以上"的方法？下一节的谱方法把这个问题推到指数收敛的极致，也把第一章傅里叶的遗产变成计算引擎。