MMS：不可压缩纳维-斯托克斯方程的验证
理论与物理
MMS是什么
老师，我听说MMS（制造解法）是“制造假解”，为什么要特意制造假解呢？这有意义吗？
问得好。MMS（Method of Manufactured Solutions）是代码验证的最强工具。关键点不是“假解”，而是“自己选择的解”。
像纳维-斯托克斯方程这样复杂的方程，通常不存在解析解。但使用MMS，通过将任意光滑函数“指定为解”，可以严格检查求解器的实现是否在数学上正确。
啊，“指定为解”是什么意思？可以随便说一个不满足方程的函数是解吗？
这就是MMS的窍门。例如，将纳维-斯托克斯方程的微分算子记为 $\mathcal{L}$，原本是寻找满足
$$ \mathcal{L}(u) = 0 $$的 $u$。但在MMS中，我们选择自己喜欢的 $u_\text{mfg}$，并将其代入方程。显然 $\mathcal{L}(u_\text{mfg}) \neq 0$，所以将其残差作为源项 $f$ 加到右边：
$$ \mathcal{L}(u_\text{mfg}) = f $$也就是说，人工构造一个“在存在 $f$ 这个源项的世界里，$u_\text{mfg}$ 是精确解”的问题。让求解器解这个问题，并确认其解是否与 $u_\text{mfg}$ 一致。
原来如此！是通过反算源项来使逻辑自洽啊。但是，这样真的能找到bug吗？
能找到。因为要确认收敛阶数（order of accuracy）。如果用二阶精度格式进行离散化，那么网格宽度 $h$ 减半时，误差应该变为 $1/4$。如果这个规律被破坏，就说明代码某处有bug。
例如，有汽车制造商的外饰CFD代码中，对流项的离散化存在符号错误。在通常的基准测试（如Lid-driven cavity等）中，结果看起来“大致正确”，但用MMS发现收敛阶数降到了一阶，从而找到了bug。
不可压缩纳维-斯托克斯方程
那么，请告诉我应用MMS的目标方程。不可压缩纳维-斯托克斯方程是什么形式？
不可压缩NS方程由动量守恒方程和连续性方程（质量守恒）两者组成。二维形式如下：
动量守恒方程（$x$ 方向）:
动量守恒方程（$y$ 方向）:
连续性方程:
各项的物理意义
- 非定常项 $\partial u / \partial t$：流速的时间变化率。定常问题中为零。
- 对流项 $(u \cdot \nabla)u$：流体自身输运动量的非线性效应。高雷诺数时占主导。例如河流汇合点处流速加快的现象。
- 压力梯度项 $-\nabla p / \rho$：压力差推动流体的力。是泵的驱动力或翼面升力的来源。
- 粘性项 $\nu \nabla^2 u$：分子粘性引起的动量扩散。像蜂蜜那样缓慢扩散的流动。
- 源项 $f$：重力或外力。在MMS中，这里会加入“使逻辑自洽”的项。
假设条件与适用范围
- 不可压缩（马赫数 $\ll$ 0.3）：密度 $\rho$ 恒定
- 牛顿流体：粘性应力与应变速率线性正比
- 等温：忽略温度变化引起的物性变化
- 连续介质假设：Knudsen数 $\ll$ 1（分子平均自由程远小于特征尺寸）
制造解的构建
将MMS应用于NS方程时，制造解怎么选？随便什么都可以吗？
理论上什么都可以，但实用上有必须遵守的规则：
- 满足连续性方程：必须满足 $\nabla \cdot \mathbf{u}_\text{mfg} = 0$，否则会给不可压缩求解器提供矛盾的输入
- 足够光滑：如果不可微，则无法验证高阶精度格式
- “激活”所有项：例如，如果制造解是常数，那么对流项和粘性项都为零，测试就失去了意义
具有代表性的是Taylor-Green涡。这实际上也是NS方程的精确解（无需源项即可成立），但非常适合用来解释MMS。
Taylor-Green涡型制造解（2D）:
如何确认这个制造解满足连续性方程呢？
实际计算偏导数即可：
$$ \frac{\partial u_\text{mfg}}{\partial x} = \pi\sin(\pi x)\sin(\pi y)\,e^{-2\pi^2\nu t} $$ $$ \frac{\partial v_\text{mfg}}{\partial y} = -\pi\sin(\pi x)\sin(\pi y)\,e^{-2\pi^2\nu t} $$相加为零。完美满足连续性方程。这个制造解中 $u$ 和 $v$ 的结构是“符号反转+sin/cos互换”的设计，所以总能满足不可压缩条件。
源项的推导
那么，终于要计算源项了。是把制造解代入NS方程对吧？
没错。依次计算各项。以 $x$ 方向的动量方程为例说明。
首先是非定常项：
$$ \frac{\partial u_\text{mfg}}{\partial t} = 2\pi^2\nu\cos(\pi x)\sin(\pi y)\,e^{-2\pi^2\nu t} $$对流项（$u \cdot \partial u/\partial x + v \cdot \partial u/\partial y$）：
$$ u\frac{\partial u}{\partial x} = -\pi\cos(\pi x)\sin(\pi y)\cdot\pi\sin(\pi x)\sin(\pi y)\cdot e^{-4\pi^2\nu t} $$ $$ v\frac{\partial u}{\partial y} = \pi\sin(\pi x)\cos(\pi y)\cdot(-\pi)\cos(\pi x)\cos(\pi y)\cdot e^{-4\pi^2\nu t} $$粘性项：
$$ \nu\nabla^2 u = \nu\cdot 2\pi^2\cos(\pi x)\sin(\pi y)\,e^{-2\pi^2\nu t} $$非定常项与粘性项相互抵消，对于Taylor-Green涡，源项 $f_x = 0$。也就是说这个制造解同时也是NS方程的精确解。
诶？源项为零？那还能作为MMS使用吗？
可以使用。即使 $f=0$，“已知精确解”这一点没有改变，所以收敛阶数的验证完全可以进行。不过在实际工作中，更常使用更一般的制造解。例如：
$$ u_\text{mfg} = \sin(ax)\sin(by) + c $$ $$ v_\text{mfg} = \cos(ax)\cos(by) + d $$选择这种不满足连续性方程的形式，通过调整系数 $a, b$ 也是一种方法。这种情况下源项非零，因此可以同时测试源项的实现。
这里出于教学目的使用Taylor-Green涡，但在实际工作中，使用更复杂的制造解（多项式+三角函数的组合）来激活所有离散化项是铁则。
MMS步骤总结：
- 选择制造解 $u_\text{mfg}$（设计成满足连续性方程）
- 代入NS方程的微分算子 $\mathcal{L}$ 计算源项 $f = \mathcal{L}(u_\text{mfg})$
- 将源项 $f$ 添加到求解器右边，并将 $u_\text{mfg}$ 设为边界条件
- 在多种网格分辨率下求解，计算与 $u_\text{mfg}$ 的误差
- 确认收敛阶数 $p$ 是否与离散化格式的理论值一致
收敛阶数的评估
收敛阶数具体是怎么得出数字的？
系统地细化网格，计算每个网格下的误差 $e$。两个网格之间的观测收敛阶数为：
其中 $e_1, e_2$ 分别是网格宽度 $h_1, h_2$ 下的误差范数。
能给我看个具体例子吗？比如二阶精度格式会怎么样？
来看实际的数值例子。将网格细化为 $h$, $h/2$, $h/4$, $h/8$ 的情况：
| 网格 | 单元数 | $h$ | $L_2$ 误差范数 | 观测阶数 $p$ | 判定 |
|---|---|---|---|---|---|
| Level 1 | $16 \times 16$ | $1/16$ | $3.87 \times 10^{-3}$ | -- | -- |
| Level 2 | $32 \times 32$ | $1/32$ | $9.72 \times 10^{-4}$ | 1.99 | PASS |
| Level 3 | $64 \times 64$ | $1/64$ | $2.43 \times 10^{-4}$ | 2.00 | PASS |
| Level 4 | $128 \times 128$ | $1/128$ | $6.08 \times 10^{-5}$ | 2.00 | PASS |
判定标准: 观测阶数在理论值（2.0）的 $\pm 0.1$ 以内则为 PASS。若系统性偏低，则怀疑有bug。
也写下 $L_2$ 范数的计算公式：
离散形式为：
其中 $V_i$ 是单元 $i$ 的体积（面积）。
数值解法与实现
有限体积法离散化
CFD的话，主流是有限体积法（FVM）吧。MMS的源项在FVM中怎么处理？
FVM处理的是单元平均值，所以源项也要在单元体积上积分。单元 $i$ 的离散化动量方程为：
$$ \sum_f \phi_f u_f - \sum_f \nu (\nabla u)_f \cdot \mathbf{S}_f = -\frac{1}{\rho}\sum_f p_f \mathbf{S}_f + \int_{V_i} f \, dV $$其中 $\phi_f = \mathbf{u}_f \cdot \mathbf{S}_f$ 是面通量，$\mathbf{S}_f$ 是面法向量。源项的体积积分 $\int_{V_i} f \, dV$ 近似为单元中心值 $f_i \cdot V_i$（单元中心单点近似）。
单点近似精度够吗？
问得好。单元中心单点近似是二阶精度，所以对于二阶精度格式的验证没有问题。但是，如果要验证三阶及以上精度，则需要用高斯求积对源项进行高精度积分。这里是容易忽略的点，请注意。
源项的实现
源项怎么添加到求解器里？
主要有两种方法：
- 方法1: 显式添加源项 -- 在动量方程右边添加用户定义的源项函数。OpenFOAM用
fvOptions，Fluent用UDF的DEFINE_SOURCE宏。 - 方法2: 作为体积力添加 -- 以重力项相同的形式给出源项。实现简单，但有时需要另外处理压力项的源项。
实际工作中方法1更普遍。将源项函数硬编码，或实现为以坐标和时间作为参数的函数。
误差范数的计算
最好计算 $L_1$, $L_2$, $L_\infty$ 三种误差范数：
计算三种有什么意义？只用 $L_2$ 不行吗？
不能说不行，但存在 $L_\infty$（最大误差）范数崩溃的情况。例如，仅在边界附近精度下降的bug，在 $L_2$ 范数中容易被平均化而难以发现，但用 $L_\infty$ 范数就能立刻发现。如果三种范数都得出正确的收敛阶数，那么代码的可靠性就相当高了。
使用SymPy自动推导源项
源项的手工计算，感觉在对流项的微分上很容易出错，有点害怕...
实际工作中，用SymPy（Python的符号计算库）自动推导是铁则。手工计算仅用于验算。特别是非线性对流项 $(u \cdot \nabla)u$ 的展开中符号错误频发，强烈建议用符号计算交叉检查。参考代码：
import sympy as sp
x, y, t, nu = sp.symbols('x y t nu')
pi = sp.pi
# 制造解的定义
u = -sp.cos(pi*x) * sp.sin(pi*y) * sp.exp(-2*pi**2*nu*t)
v = sp.sin(pi*x) * sp.cos(pi*y) * sp.exp(-2*pi**2*nu*t)
p = -sp.Rational(1,4)*(sp.cos(2*pi*x)+sp.cos(2*pi*y))*sp.exp(-4*pi**2*nu*t)
# 连续性方程的确认
div = sp.diff(u, x) + sp.diff(v, y)
print("div(u) =", sp.simplify(div)) # =>なった
詳しく
報告