MMS: 非圧縮性方程式的验证

分类: V&V / 制造解法 | 统一版 2026-04-12

MMS verification of incompressible Navier-Stokes equations showing Taylor-Green vortex velocity field and convergence rate plot

制造解法(MMS)对Navier-Stokes方程求解器的验证 -- Taylor-Green涡的制造解和收敛阶评估

MMS的理论基础

MMS是什么

🧑🎓

老师，我听说MMS（制造解法）是"制造虚假的解"，但我不太理解为什么要制造虚假解呢？这样做有意义吗？

🎓

很好的问题。MMS（Method of Manufactured Solutions）是代码验证的最强工具。关键是理解"虚假解"这个概念——它实际上是"自己选择的解"。

对于Navier-Stokes方程这样的复杂方程，通常没有分析解。但MMS允许我们选择任意光滑函数作为"解"，然后反推源项，以严格数学方法验证求解器的实现是否正确。

🧑🎓

等等，"将函数视为解"是什么意思？不满足方程的函数也能称为解吗？

🎓

这就是MMS的妙处。设微分算子为 $\mathcal{L}$，通常我们要找满足以下条件的 $u$ ：

$$ \mathcal{L}(u) = 0 $$

但在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。如果这个关系破坏了，就说明代码中存在错误。

我记得一个例子：一个汽车制造商的外观CFD代码中对流项离散化有符号错误。在标准基准测试（如Lid-driven cavity）中，结果看起来"还不错"，但MMS检验发现收敛阶下降到一阶，因此成功发现了这个bug。

非压缩性NS方程

🧑🎓

那么，MMS要应用的目标方程是什么？非压缩性Navier-Stokes方程的形式如何？

🎓

非压缩性NS方程由动量守恒方程和连续性方程（质量守恒）组成。在二维中写成如下形式：

动量守恒方程（$x$ 方向）：

$$ \frac{\partial u}{\partial t} + u\frac{\partial u}{\partial x} + v\frac{\partial u}{\partial y} = -\frac{1}{\rho}\frac{\partial p}{\partial x} + \nu\left(\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2}\right) + f_x $$

动量守恒方程（$y$ 方向）：

$$ \frac{\partial v}{\partial t} + u\frac{\partial v}{\partial x} + v\frac{\partial v}{\partial y} = -\frac{1}{\rho}\frac{\partial p}{\partial y} + \nu\left(\frac{\partial^2 v}{\partial x^2} + \frac{\partial^2 v}{\partial y^2}\right) + f_y $$

连续性方程：

$$ \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} = 0 $$

制造解的构造

🧑🎓

在NS方程中应用MMS时，制造解该如何选择？可以随意选择吗？

🎓

理论上可以任意选择，但实务中有需要遵守的规则：

满足连续性方程：必须 $\nabla \cdot \mathbf{u}_\text{mfg} = 0$，否则会对非压缩求解器提供自相矛盾的输入
充分光滑：函数必须可微，否则无法验证高阶精度方案
激活所有项：例如，常数制造解会使对流项和粘性项都为零，失去测试意义

典型的例子是Taylor-Green涡。它实际上是NS方程的严格解（无源项成立），是MMS演示的理想选择。

Taylor-Green涡型制造解（2D）：

$$ u_\text{mfg}(x,y,t) = -\cos(\pi x)\sin(\pi y)\,e^{-2\pi^2 \nu t} $$

$$ v_\text{mfg}(x,y,t) = \sin(\pi x)\cos(\pi y)\,e^{-2\pi^2 \nu t} $$

$$ p_\text{mfg}(x,y,t) = -\frac{1}{4}\bigl(\cos(2\pi x) + \cos(2\pi y)\bigr)\,e^{-4\pi^2 \nu t} $$

🧑🎓

这个制造解如何验证满足连续性方程？

🎓

直接计算偏导数即可：

$$ \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 $$

这些不满足连续性方程，导致非零源项，同时也测试了源项实现。

教学目的上我们用Taylor-Green涡，但实务中应该使用更复杂的制造解（多项式+三角函数组合），激活所有离散化项，这是铁则。

收敛阶的评估

🧑🎓

收敛阶具体如何计算出来？

🎓

系统地细化网格，计算各网格的误差 $e$。两个网格间的观测收敛阶为：

$$ p = \frac{\log(e_1 / e_2)}{\log(h_1 / h_2)} $$

其中 $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	通过
Level 3	$64 \times 64$	$1/64$	$2.43 \times 10^{-4}$	2.00	通过
Level 4	$128 \times 128$	$1/128$	$6.08 \times 10^{-5}$	2.00	通过

判定标准：观测阶在理论值（2.0）的 $\pm 0.1$ 范围内为通过。系统性低于理论值时怀疑存在bug。

🎓

$L_2$ 范数的计算公式也写下来：

$$ \|e\|_{L_2} = \sqrt{\frac{1}{\Omega}\int_\Omega \bigl(u_\text{numerical} - u_\text{mfg}\bigr)^2 \, d\Omega} $$

离散形式：

$$ \|e\|_{L_2} \approx \sqrt{\frac{\sum_{i=1}^{N} (u_i - u_{\text{mfg},i})^2 \cdot V_i}{\sum_{i=1}^{N} V_i}} $$

其中 $V_i$ 是单元 $i$ 的体积（或面积）。

MMS的数值计算方法

有限体积法的离散化

🧑🎓

在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$ ：

$$ \|e\|_{L_1} = \frac{1}{N}\sum_{i=1}^{N} |u_i - u_{\text{mfg},i}| $$

$$ \|e\|_{L_\infty} = \max_{i} |u_i - u_{\text{mfg},i}| $$

🧑🎓

为什么要计算三种？只用 $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))  # => 0

# x方向源项: f_x = du/dt + u*du/dx + v*du/dy + dp/dx - nu*(d2u/dx2 + d2u/dy2)
fx = (sp.diff(u,t) + u*sp.diff(u,x) + v*sp.diff(u,y)
      + sp.diff(p,x) - nu*(sp.diff(u,x,2) + sp.diff(u,y,2)))
print("f_x =", sp.simplify(fx))  # => 0 (Taylor-Green涡情况)

MMS的实务应用

分步骤操作指南

🧑🎓

实际执行MMS时，从头到尾需要做什么？

🎓

按以下7个步骤：

选定制造解：选择满足连续性方程的光滑函数。确认所有方程项都被激活
推导源项：用SymPy等符号计算。与手工计算交叉检查
设置计算区域和边界条件：边界处使用制造解值作为Dirichlet条件
系统地生成网格：至少4层网格（如16x16, 32x32, 64x64, 128x128）。网格比最好是2倍
各网格上运行计算：相同求解器设置。残差充分收敛
计算误差范数：$L_1$, $L_2$, $L_\infty$ 三种。分别对速度 $u$, $v$ 和压力 $p$ 计算
判断收敛阶：$p = \log(e_1/e_2) / \log(h_1/h_2)$ 是否与理论值一致

🧑🎓

边界条件全部用Dirichlet吗？压力的边界条件怎么处理？

🎓

速度用Dirichlet（固定为制造解值）。压力需要小心：非压缩求解器通常在一点固定压力（参考压力），或全边界上使用Neumann条件（从制造解计算压力梯度）。全边界Dirichlet会导致over-constrain，要查看求解器文档。

OpenFOAM实现例

🧑🎓

OpenFOAM中具体怎么做？

🎓

OpenFOAM中最简单的方法是用 icoFoam（非压缩层流求解器）加 fvOptions 添加源项。步骤：

blockMeshDict 生成 $[0,1] \times [0,1]$ 正方形区域
0/U 用 codedFixedValue 设置制造解初始条件和边界条件
constant/fvOptions 用 codedSource 定义源项
4层网格运行，postProcess -func 'error' 收集误差

Taylor-Green涡时源项为零，fvOptions 可不用。但通用制造解时源项实现是必需的。

Ansys Fluent实现

🎓

Fluent中使用UDF（用户定义函数）：

DEFINE_SOURCE(x_mom_source, ...)：$x$ 动量源项
DEFINE_PROFILE(u_bc, ...)：边界条件的制造解速度值
DEFINE_INIT(mms_init, ...)：初始条件设为制造解

Fluent用C语言写UDF，源项数式可用SymPy的 sp.ccode() 函数直接转换为C代码。

品质保证检查表

🧑🎓

执行MMS时应检查哪些要点？

🎓

以下项全部通过就可判定MMS检验合格：

#	检查项目	合格标准
1	制造解满足连续性方程	$\nabla \cdot \mathbf{u}_\text{mfg} = 0$（符号验证）
2	源项推导正确	SymPy和手工计算一致
3	边界条件与制造解一致	边界上 $u_\text{num} = u_\text{mfg}$
4	各网格残差充分收敛	残差 $< 10^{-10}$（离散化误差远大于舍入误差）
5	$L_2$ 范数收敛阶与理论值一致	$\|p_\text{obs} - p_\text{theory}\| < 0.1$
6	$L_\infty$ 范数收敛阶也检验	同上
7	速度 $u$, $v$ 和压力 $p$ 全部检验	所有变量均合格

MMS的软件比较

CFD求解器MMS支持情况

🧑🎓

不同求解器之间MMS的易用性有差异吗？

🎓

差异很大。关键在于能否灵活添加源项。

求解器	源项添加方法	MMS友好特性	难度
OpenFOAM	`fvOptions` / `codedSource`	完整源代码访问，自定义求解器	低
Ansys Fluent	UDF (`DEFINE_SOURCE`)	C语言UDF，Scheme脚本	中
STAR-CCM+	Field Function / User Code	Java API，表格输入	中
COMSOL	方程形式建模	任意PDE右侧直接输入源项	低
FEniCS/Firedrake	弱形式直接定义	Python符号定义，UFL记法	低
SU2	源代码修改	内置MMS检验框架	中

🧑🎓

OpenFOAM和COMSOL最简单啊。商用求解器UDF开发环境成了障碍...

🎓

完全同意。不过Fluent 2024R2以后开始支持Python UDF，门槛降低了不少。SU2是开源的，MMS检验甚至被纳入官方CI/CD回归测试。重视代码验证文化的求解器，MMS支持往往更完善。

MMS的前沿研究

三维MMS的扩展

🧑🎓

到目前为止都是二维。实务的CFD是三维呀。三维扩展难吗？

🎓

原理相同，但连续性方程 $\nabla \cdot \mathbf{u} = 0$ 要在三个分量上满足。三维Taylor-Green涡为：

$$ u = A\cos(ax)\sin(by)\sin(cz) $$ $$ v = B\sin(ax)\cos(by)\sin(cz) $$ $$ w = C\sin(ax)\sin(by)\cos(cz) $$

其中系数满足 $Aa + Bb + Cc = 0$。三维时源项推导的项数爆炸增长，SymPy几乎成为必须工具。

非定常问题的时间精度验证

🎓

时间积分方案精度也能同步验证。步骤：

空间网格固定在足够细的等级（空间误差不主导）
系统变化时间步长 $\Delta t$ 为 $\Delta t$, $\Delta t/2$, $\Delta t/4$...
在最终时刻计算误差，检验时间方向收敛阶

Euler隐格式期望一阶，Crank-Nicolson期望二阶，BDF2期望二阶。空间和时间分开验证很关键，同时细化会混淆是哪项精度降低。

湍流模型的应用

🧑🎓

RANS湍流模型（比如Spalart-Allmaras）也能用MMS？

🎓

完全可以。事实上湍流模型正是MMS大显身手的地方。因为湍流模型的实现中源项（生成项、消散项）众多，容易藏bug。

Spalart-Allmaras模型情形下，对湍流粘性 $\tilde{\nu}$ 的输运方程也设置制造解，导出结合NS方程的源项。Rumsey（NASA LaRC）的公开验证数据是实际标准。

MMS的故障排除

收敛阶低于理论值

🧑🎓

老师，我的二阶方案收敛阶只有1.5，这是bug吗？

🎓

可能是bug，但先检查以下几点：

网格是否足够细？ 过粗的网格不会进入asymptotic range，高阶截断误差项主导，收敛阶不符理论。提高分辨率重试
迭代残差收敛充分吗？ 残差 $10^{-6}$ 时，细网格上离散误差和迭代误差同级，破坏收敛阶。降至 $10^{-10}$ 以下
源项实现有舍入误差？ 用 double 而非 float
边界处理精度如何？ 内部二阶、边界一阶会降低总体精度

🧑🎓

啊，我残差设的 $10^{-6}$ 啦...那可能就是这个原因了。

🎓

在MMS中"离散化误差 >> 迭代误差"这条件是关键。通常CFD工程计算 $10^{-6}$ 足够，但MMS验证必须极限收敛。初学者这点最容易忘。

源项的不匹配

🧑🎓

源项实现后收敛阶变成零（网格细化完全没用），完全坏了...

🎓

零阶收敛="网格细化误差不减"，几乎肯定是源项计算错误。常见原因：

符号错误：对流项 $(u \cdot \nabla)u$ 展开中 $+/-$ 反号
密度处理错误：假设 $\rho = 1$，但源项里乘了 $\rho$（或反过来）
坐标系不一致：源项的 $(x,y)$ 与求解器内部坐标系错位
时刻评估错位：非定常问题中源项的时刻是时间步起点还是终点

排查方法：先用Taylor-Green涡（零源项、已知严格解）验证收敛阶。如果正确，再加源项。这样就能定位源项实现的问题。

边界条件的陷阱

🧑🎓

边界条件还有其他需要注意的？

🎓

有几个典型陷阱：

周期边界条件：Taylor-Green涡与周期BC契合，但制造解必须真的周期。否则边界不连续，收敛阶崩溃
压力参考点：非压缩求解器中压力有常数不确定性。MMS中计算值和制造解要先减去各自平均值再比较
壁函数干扰：湍流计算若用壁函数，则壁附近模型化与MMS严格解矛盾。应用"分辨到壁"（resolve to wall）验证

🧑🎓

MMS真是深啊...但只要这样扎实做，CFD代码就能获得很高的信度了！

🎓

就是这样。MMS是"代码正确的数学证明"最接近的方法。ASME V&V 20和AIAA V&V指南都把MMS列为推荐的代码验证方法。新开发求解器，或给现有求解器加新功能，都必须用MMS验证。这习惯一旦养成，就能显著提升数值模拟的信度。