3. 有限体积离散的基本框架
有限体积法(Finite Volume Method, FVM)的核心思想可以用一句话概括:把守恒律在每一个控制体积上做积分,再用界面通量把相邻单元耦合起来。与有限差分法直接在节点上近似导数、有限元法在单元上做加权残量不同,有限体积法从积分形式的守恒方程出发,天然保证离散后的通量在内部界面上成对抵消,从而在离散层面严格守恒。本节围绕这一思想,把"从连续方程到半离散常微分方程组"的构造过程拆解为四个关键环节:控制体积划分、网格与单元、界面通量、半离散格式组装。
一、从积分守恒律出发
以一维标量守恒律为例,连续形式为
∂u/∂t + ∂f(u)/∂x = 0, f(u) 为通量函数
在任意控制体积 Vi = [xi−1/2, xi+1/2] 上对空间积分,并利用高斯散度定理把通量的散度项化为边界积分,得到
d/dt ∫_{V_i} u dx + f(u_{i+1/2}) − f(u_{i−1/2}) = 0
这一步是有限体积法的"根":方程左边第一项是控制体积内守恒量的变化率,第二、三项是穿过左右界面的净通量。注意此时还没有做任何数值近似——它是精确的积分关系。接下来的全部工作,就是把积分量、界面通量分别用离散未知量表示出来。
二、控制体积划分与网格、单元
控制体积是守恒律被施加的"账本单元"。一维情形下,最常用的是单元中心型(cell-centered)划分:把计算域剖分为 N 个互不重叠的区间,每个区间既是一个网格单元(cell),也是一个控制体积,未知量存放在单元中心 xi。与之相对的是节点型(vertex-centered)划分,未知量在网格节点上,控制体积围绕节点构造,界面落在相邻节点中点。二维、三维情形同理,控制体积可以是矩形、多边形或任意多面体。
设单元 i 的长度为 Δxi,定义单元平均值
ū_i(t) = (1/Δx_i) ∫_{V_i} u(x,t) dx
把积分守恒式两边除以 Δxi,得到
dū_i/dt = −(1/Δx_i) [ f(u_{i+1/2}) − f(u_{i−1/2}) ]
这就是有限体积法的半离散形式:时间导数保留连续,空间方向已经离散为一个由界面通量构成的右端项。所谓"半离散",正是指只离散了空间、时间仍待积分。这里要特别注意符号约定:本文统一取 x 轴向右为正,通量 f 的正方向为 +x,因此单元右界面的流出通量取正号、左界面的流入通量取负号,二者之差构成净流出。若采用相反约定,全篇公式与图示必须同步翻转,不能混用。
三、界面通量的数值近似
半离散式中唯一尚未确定的是界面通量 f(ui+1/2)。界面两侧各有一个单元,物理上通量应由两侧状态共同决定,记作数值通量
F_{i+1/2} = F̂(ū_i, ū_{i+1})
数值通量必须满足两个基本性质:相容性 F̂(a, a) = f(a),即两侧状态相同时退化为物理通量;守恒性,即同一界面在左右两个单元的方程中取同一个 F̂ 值,符号相反,从而在求和时内部界面通量自动抵消。正是守恒性保证了有限体积格式的全局守恒。
常用的数值通量包括:中心平均 F̂ = (fi + fi+1)/2,简单但无耗散、易产生非物理振荡;迎风通量,按特征方向取上游值,稳定但一阶精度耗散较大;以及 Godunov 型、Roe 型、HLL 型等近似 Riemann 解通量,在双曲守恒律问题中兼顾稳定与精度。选择哪种通量,直接决定格式的稳定性、精度和适用场景,是有限体积法中最需要工程判断的一环。
四、半离散格式的构造步骤
把上述环节串起来,一个完整的有限体积半离散格式按以下顺序构造:
- 确定守恒方程与通量函数:明确守恒量 u、通量 f(u) 及其特征速度,判断方程类型(双曲、抛物或混合)。
- 划分网格与控制体积:生成单元、节点、界面拓扑,记录每个单元的体积(一维为 Δxi)、每个界面的左右邻居编号。对非结构网格,还需存储界面法向与面积。
- 在控制体积上积分:把守恒律写成积分形式,用散度定理化为界面通量之和,得到每个单元的积分平衡关系。
- 定义单元平均值并引入数值通量:用 ūi 表示积分量,用 F̂ 表示界面通量,得到 dūi/dt 的表达式。
- 选择数值通量格式:根据方程类型与精度要求选取迎风、中心、Roe、HLL 等格式,并检查相容性与守恒性。
- 组装右端项:对每个单元累加所有界面的通量贡献,形成全局半离散系统 dū/dt = R(ū),其中 R 为残差算子。
- 施加边界条件:在边界界面上用虚拟单元或指定通量封闭方程,注意边界通量的符号与内部界面一致。
以一个具体的小例子说明组装过程。设一维域 [0, 1] 均匀剖分为 4 个单元,Δx = 0.25,初始 ū = [1.0, 0.8, 0.5, 0.2],通量 f(u) = u²/2,采用一阶迎风通量(特征速度 u > 0,取左侧值)。对单元 2(ū2 = 0.8),左界面通量 F̂3/2 = f(ū1) = 0.8²/2 = 0.32,右界面通量 F̂5/2 = f(ū2) = 0.8²/2 = 0.32,于是 dū2/dt = −(1/0.25)(0.32 − 0.32) = 0。这说明在该初始分布下单元 2 暂时处于局部平衡;而单元 3(ū3 = 0.5)的右端项为 −(1/0.25)(f(0.5) − f(0.8)) = −4×(0.125 − 0.32) = 0.78 > 0,守恒量将增加。这一结果符合物理直觉:高值区向低值区输运。需要强调,以上数值仅为教学示例,用于演示组装与符号判断,并非实际计算结果。
五、常见易错点
初学者在构造半离散格式时最容易在三个地方出错。其一是符号约定不统一:右界面流出取正、左界面流入取负,一旦某个界面写反,守恒性立刻被破坏。其二是界面通量重复计算:同一界面在左右两个单元中各出现一次,必须保证取同一个 F̂ 值,不能各自用不同的近似。其三是单元体积与界面面积混淆:在非均匀或曲线网格上,除以的是单元体积而非界面面积,二者量纲不同,混用会导致格式量纲不一致。把这三条作为组装后的自检清单,可以避免绝大多数低级错误。
至此,连续守恒律已经被转化为一个关于单元平均值的半离散常微分方程组。空间离散的全部工作已经完成,剩下的问题是如何对这个方程组做时间积分——这正是后续时间推进格式要解决的问题。