数值抛物型方程的有限体积法
好的，我们开始学习“数值抛物型方程的有限体积法”。这个方法结合了有限差分法的直观和有限元法处理复杂几何区域的优势，是计算流体力学和传热学等领域中求解抛物型方程（如热传导方程、对流-扩散方程）的核心技术。
第一步：理解物理背景与控制方程
我们从一个典型的抛物型方程——瞬态热传导方程开始：
\[\frac{\partial u}{\partial t} = \alpha \nabla^2 u + S \]
这里：
- \(u\) 是待求解的物理量（如温度）。
- \(t\) 是时间。
- \(\alpha\) 是扩散系数（大于0）。
- \(\nabla^2\) 是拉普拉斯算子，代表扩散项。
- \(S\) 是源项（如内部热源）。
这个方程描述了物理量（如热量）随时间扩散的过程。有限体积法的核心思想不是直接近似这个微分方程，而是忠实于其背后的物理守恒定律。
第二步：有限体积法的核心思想——积分守恒
有限体积法的基本出发点是将计算区域划分为许多互不重叠的控制体积（或称为网格单元）。对于每个控制体积，我们要求所关注的物理量（如能量）满足守恒定律。
- 空间离散： 将求解域 \(\Omega\) 划分为 \(N\) 个控制体积 \(\Omega_i\) (\(i = 1, 2, ..., N\))。
- 积分方程： 对上述热传导方程在任意一个控制体积 \(\Omega_i\) 上积分，并应用高斯散度定理：
\[ \int_{\Omega_i} \frac{\partial u}{\partial t} dV = \alpha \int_{\Omega_i} \nabla^2 u dV + \int_{\Omega_i} S dV \]
\[ \int_{\Omega_i} \frac{\partial u}{\partial t} dV = \alpha \oint_{\partial \Omega_i} (\nabla u) \cdot \mathbf{n} dA + \int_{\Omega_i} S dV \]
这里：
-
\(\partial \Omega_i\) 是控制体积 \(\Omega_i\) 的边界。
-
\(\mathbf{n}\) 是边界上指向外侧的单位法向量。
-
\(dA\) 是面积微元。
这个积分形式的物理意义非常明确：
-
左边：控制体内物理量 \(u\) 随时间的变化率。
-
右边第一项：通过控制体边界 \(\partial \Omega_i\) 的净通量（例如，净热流）。\(\nabla u \cdot \mathbf{n}\) 是法向梯度，与流量相关。
- 右边第二项：控制体内的源项总和。
因此，方程表达了：控制体内物理量的增加 = 流入的净通量 + 内部产生的量。这才是最根本的守恒形式。
第三步：一维情况下的详细推导
为了更具体，我们考虑一维无源项的热传导方程：
\[\frac{\partial u}{\partial t} = \alpha \frac{\partial^2 u}{\partial x^2} \]
计算域为 \([a, b]\)，将其划分为多个控制体积。对于第 \(i\) 个控制体积（一个小区间），其中心点（网格点）为 \(x_i\)，左右边界分别为 \(x_{i-1/2}\) 和 \(x_{i+1/2}\)。
- 对控制体积积分：
\[ \int_{x_{i-1/2}}^{x_{i+1/2}} \frac{\partial u}{\partial t} dx = \alpha \int_{x_{i-1/2}}^{x_{i+1/2}} \frac{\partial^2 u}{\partial x^2} dx \]
- 处理各项：
- 时间项：假设 \(u\) 在控制体内平均值 \(\bar{u}_i\) 代表该控制体的值，则左边近似为 \(\frac{d\bar{u}_i}{dt} \Delta x_i\)，其中 \(\Delta x_i = x_{i+1/2} - x_{i-1/2}\)。
- 扩散项（通量项）：利用牛顿-莱布尼茨公式，右边积分可精确求出：
\[ \alpha \int_{x_{i-1/2}}^{x_{i+1/2}} \frac{\partial^2 u}{\partial x^2} dx = \alpha \left[ \frac{\partial u}{\partial x} \bigg|_{x_{i+1/2}} - \frac{\partial u}{\partial x} \bigg|_{x_{i-1/2}} \right] \]
这正好是通过右边界的通量减去左边界的通量，即**净通量**。
- 得到半离散方程：
\[ \Delta x_i \frac{d\bar{u}_i}{dt} = \alpha \left( \frac{\partial u}{\partial x} \bigg|_{i+1/2} - \frac{\partial u}{\partial x} \bigg|_{i-1/2} \right) \]
现在，关键问题是如何计算界面 \(x_{i±1/2}\) 处的导数（即通量）\(\frac{\partial u}{\partial x}\)。
第四步：界面通量的计算——格式的核心
界面上的通量需要由相邻控制体的值来构造。最常用的方法是中心格式。
- 在界面 \(i+1/2\) 处，我们近似：
\[ \frac{\partial u}{\partial x} \bigg|_{i+1/2} \approx \frac{u_{i+1} - u_i}{x_{i+1} - x_i} = \frac{u_{i+1} - u_i}{\Delta x} \]
（这里假设均匀网格，\(\Delta x\) 为网格点间距）。
- 同理，界面 \(i-1/2\) 处：
\[ \frac{\partial u}{\partial x} \bigg|_{i-1/2} \approx \frac{u_i - u_{i-1}}{\Delta x} \]
- 得到全离散方程：
将通量表达式代入半离散方程：
\[ \Delta x \frac{d u_i}{dt} = \alpha \left( \frac{u_{i+1} - u_i}{\Delta x} - \frac{u_i - u_{i-1}}{\Delta x} \right) \]
\[ \frac{d u_i}{dt} = \alpha \frac{u_{i+1} - 2u_i + u_{i-1}}{\Delta x^2} \]
你会发现，这正好是用二阶中心差分离散扩散项 \(\frac{\partial^2 u}{\partial x^2}\) 得到的结果。这说明在一维扩散问题上，有限体积法（中心格式）和有限差分法是等价的。
第五步：时间离散
上面得到的是一个关于时间的常微分方程组 \(\frac{d\mathbf{u}}{dt} = \mathbf{F}(\mathbf{u})\)。我们需要用时间积分方法求解。
- 显式欧拉法：最简单，但有稳定性限制（CFL条件）。
\[ u_i^{n+1} = u_i^n + \frac{\alpha \Delta t}{\Delta x^2} (u_{i+1}^n - 2u_i^n + u_{i-1}^n) \]
稳定性要求 \(\frac{\alpha \Delta t}{\Delta x^2} \le \frac{1}{2}\)。
- 隐式方法（如Crank-Nicolson或全隐式）：无条件稳定，但需要求解线性方程组。例如，全隐式欧拉法：
\[ u_i^{n+1} = u_i^n + \frac{\alpha \Delta t}{\Delta x^2} (u_{i+1}^{n+1} - 2u_i^{n+1} + u_{i-1}^{n+1}) \]
这需要联立所有网格点的方程，形成一个三对角线性方程组，可以用托马斯算法高效求解。
第六步：推广到多维与复杂情况
有限体积法的优势在多维和复杂边界问题上更为突出。
- 多维： 思想完全一样。在二维或三维中，对控制体积（矩形、三角形、多边形等）积分，净通量是通量向量在所有边界面（或棱）上的积分之和。通量计算需要处理每个面的法向梯度。
- 非结构网格： 有限体积法可以很自然地应用于三角形、四面体等非结构网格，非常适合复杂几何形状。
- 非线性与源项： 对于非线性问题或复杂源项，通常采用“算子分裂”或“延迟修正”等策略，将非线性项或源项当作一个独立的子问题处理。
总结
数值抛物型方程的有限体积法，其精髓在于：
- 物理守恒优先：从积分形式的守恒定律出发，保证了数值解即使在粗网格下也满足离散意义上的守恒。
- 灵活的几何适应性：通过定义控制体积和计算界面通量，可以处理各种复杂的计算区域。
- 模块化：时间离散、空间离散（通量计算）、源项处理等相对独立，便于构建复杂问题的求解器。
这种方法因其严格的守恒性和对复杂几何的良好适应性，成为计算物理和工程领域中不可或缺的数值工具。