上学期为了赶课题进度我开始啃《计算流体力学大串讲》这本教材。说实话第一次翻开的时候我是有点懵的前一页还在推导连续方程后面突然跳到各种差分格式再过几页又是一堆稳定性曲线。我当时的第一反应是把它当字典用哪个公式用到就去翻哪段结果越翻越乱做题照样卡壳。后来换了个思路把前3章当成一条完整的逻辑链来读——方程、离散、稳定性这三件事其实是同一个故事的不同侧面。这篇笔记就是按这个思路整理出来的记录了我读1-3章的推导过程、阅读路线和踩过的坑。它不是标准答案只是一个过来人把知识“串”起来的方式希望对正在入门计算流体力学的朋友有点参考价值。1. 为什么我决定先用“串讲”的方式读前3章1.1 大部分CFD教材把知识切成了孤岛市面上主流的计算流体力学教材几乎都有一个共同问题章节之间天然是断开的。讲方程的章节拼命堆数学讲离散的章节又假设你已经吃透了偏微分方程等到讲稳定性的时候前面的格式又已经忘得差不多了。再加上有限差分、有限体积、有限元这几大流派并存初学者经常陷入一种“单章能看懂合上书全忘”的状态。《计算流体力学大串讲》这个名字吸引我的地方就是它试图把这些章节之间的“接口”给补上。前3章的内容在我的理解里基本是围绕“连续方程怎么转变到代数方程代数方程又要满足什么条件才能用”这条主线展开的。如果只按目录顺序读跟在知乎上收藏一堆零零散散的公式没有区别。1.2 我的阅读策略一章三遍主线优先我给自己定了一个比较笨但很有效的策略每一章读三遍。第一遍只看标题、公式、图表目的是搞清楚这一章在回答什么问题第二遍拿着笔从第一个公式推导到最后一个公式推不动就停下来查资料第三遍合上书凭记忆把这一章的核心逻辑写在一张A4纸上。配合这个策略我给自己拉出了三条贯穿前3章的主线方程线N-S方程怎么被简化成模型方程为什么简化了还能代表真实流动。离散线偏微分方程怎么变成差分方程截断误差从哪里冒出来格式精确性由什么决定。稳定性线为什么格式看似合理跑出来的结果却是乱的CFL条件到底在管什么。这三个词贯穿了我整个阅读过程。接下来就按这条路线把前3章里我认为最核心的内容拆开讲。2. 第一根主线从N-S方程到模型方程是怎么一步一步“降维”的2.1 N-S方程到底重在哪里计算流体力学的物理起点不用多说就是纳维-斯托克斯方程。不可压缩流动的基本形式可以写成$$ \frac{\partial u_i}{\partial x_i} 0 $$$$ \rho\left(\frac{\partial u_i}{\partial t} u_j\frac{\partial u_i}{\partial x_j}\right) -\frac{\partial p}{\partial x_i} \mu \frac{\partial^2 u_i}{\partial x_j \partial x_j} f_i $$这个方程看着不复杂实际上藏着三个大麻烦。第一对流项 $u_j \partial u_i / \partial x_j$ 是非线性的速度的乘积让方程在数学上极难处理第二压力 $p$ 和速度场是耦合在一起的没有关于压力的独立方程需要靠连续性方程去约束第三真实流动里不同的尺度跨度极大从机翼绕流到管道湍流涡的尺度可以差好几个数量级。这就是为什么教材的前几章一定会花大量篇幅讲方程的性质而不是直接给格式。因为你在设计数值方法之前必须知道待求解的方程是双曲型、抛物型还是椭圆型的不同数学性质的方程对离散格式的要求完全不同。这一点我在第一次读的时候完全没意识到后来做稳定性分析才回过头来补课。2.2 模型方程不是偷懒是“小白鼠”直接拿N-S方程去验证一个新格式成本太高而且因为方程太复杂出了问题你很难判断是格式本身的问题还是方程耦合带来的问题。所以教材在第2章左右几乎都会引入几个经典模型方程模型方程数学形式对应的物理特征典型应用场景线性对流方程$\dfrac{\partial u}{\partial t} a\dfrac{\partial u}{\partial x} 0$波形以速度 $a$ 平移不衰减激波捕捉、边界条件研究扩散方程$\dfrac{\partial u}{\partial t} \alpha \dfrac{\partial^2 u}{\partial x^2}$梯度抹平信息向全空间传播热传导、粘性扩散线性对流扩散方程$\dfrac{\partial u}{\partial t} a\dfrac{\partial u}{\partial x} \alpha \dfrac{\partial^2 u}{\partial x^2}$对流与扩散同时存在粘性流体的边界层Burgers方程$\dfrac{\partial u}{\partial t} u\dfrac{\partial u}{\partial x} 0$保留非线性的最小模型可能产生间断激波、交通流这四个方程每一个都像一只“小白鼠”。线性对流方程保留了对流特性但扔掉了黏性和非线性扩散方程保留了抛物型特性但扔掉了波传播方向。你在这几个方程上验证过的格式很多性质可以直接迁移到N-S方程的求解中。比如迎风格式的耗散特性、中心差分格式的振荡问题在模型方程上表现得特别清晰放在N-S方程里反而容易被各种耦合现象掩盖。我强烈建议如果你在读前3章不要跳过错模型方程。我见过太多人直接跳过这方面内容结果后面写代码算激波管问题时连数值振荡都不知道是格式本身的属性还以为是边界条件没设对。2.3 无量纲化不是形式主义是帮你“换一副眼镜看流动”前3章还有一个看似枯燥、其实非常重要的环节无量纲化。把变量换成特征尺度以后N-S方程会变成类似这样的形式$$ \frac{\partial u_i^}{\partial x_i^}0 $$$$ \frac{\partial u_i^}{\partial t^} u_j^\frac{\partial u_i^}{\partial x_j^} -\frac{\partial p^}{\partial x_i^} \frac{1}{Re}\frac{\partial^2 u_i^}{\partial x_j^* \partial x_j^*} $$无量纲化之后出现了一个关键参数雷诺数 $Re$。它实质上是惯性力与粘性力的比。$Re$ 很小扩散项占主导方程行为趋于抛物型$Re$ 很大对流项占主导方程行为趋于双曲型。佩克莱特数 $Pe Re \cdot Pr$ 在对流扩散方程里也是同样的道理表征对流输运和扩散输运的相对强度。看懂了这一点再往后读稳定性分析就会顺很多。因为很多格式的适用条件是跟这些无量纲参数绑定的比如后面对流项的CFL数、扩散项的网格傅里叶数本质上都在描述“离散世界里信息能不能在物理时间尺度内正确传播”。无量纲化让我意识到计算流体力学不只是数值技巧它更像是在离散世界里重建一个缩小版的物理系统。这也是前3章给我的最大观念转变。3. 第二根主线有限差分是怎么把“微积分”翻译成“加减乘除”的3.1 泰勒展开是一切差分格式的起点计算机能处理的只有加减乘除没法直接求导数。有限差分的核心思想就是把导数近似成相邻网格点上函数值的线性组合。这个近似过程的理论依据就是泰勒展开。以一阶导为例。把 $u(x\Delta x)$ 在 $x$ 处泰勒展开$$ u_{i1} u_i \Delta x \left.\frac{\partial u}{\partial x}\right|_i \frac{\Delta x^2}{2}\left.\frac{\partial^2 u}{\partial x^2}\right|_i \frac{\Delta x^3}{6}\left.\frac{\partial^3 u}{\partial x^3}\right|_i \cdots $$稍微整理一下就能得到前向差分的表达式$$ \left.\frac{\partial u}{\partial x}\right|i \frac{u{i1} - u_i}{\Delta x} - \frac{\Delta x}{2}\left.\frac{\partial^2 u}{\partial x^2}\right|_i O(\Delta x^2) $$也就是说用 $\frac{u_{i1}-u_i}{\Delta x}$ 近似导数时主要误差来自 $\frac{\Delta x}{2}u$ 这一项步长越大误差线性增长所以这是一阶精度格式。用同样的方法处理 $u_{i-1}$就得到后向差分也是一阶精度。如果我们把 $u_{i1}$ 的展开式减去 $u_{i-1}$ 的展开式偶数阶项全部消掉得到的就是二阶中心差分$$ \left.\frac{\partial u}{\partial x}\right|i \frac{u{i1} - u_{i-1}}{2\Delta x} - \frac{\Delta x^2}{6}\left.\frac{\partial^3 u}{\partial x^3}\right|_i O(\Delta x^4) $$首项误差从 $O(\Delta x)$ 变成了 $O(\Delta x^2)$精度上升了一个台阶。二阶导数也类似通过对 $u_{i1}$ 和 $u_{i-1}$ 的展开式相加得到$$ \left.\frac{\partial^2 u}{\partial x^2}\right|i \frac{u{i1} - 2u_i u_{i-1}}{\Delta x^2} O(\Delta x^2) $$我在初学的时候有一个误区以为“格式精度越高结果就一定越好”。这个想法被后面的稳定性分析彻底纠正了高阶精度的格式可能带有更严重的数值振荡低阶格式反而常具备单调性不易产生非物理的波动。所以第一步是把“精度”和“稳定性”这两个概念分开来看。3.2 把离散格式写进方程一维线性对流方程实例前3章最好玩的部分就是把前面那些差分公式组合起来真正去离散一个偏微分方程。拿最简单的一维线性对流方程$$ \frac{\partial u}{\partial t} a\frac{\partial u}{\partial x} 0, \quad a 0 $$假设空间用一阶后向差分也就是迎风方向格式时间用一阶前向差分得到$$ \frac{u_j^{n1} - u_j^n}{\Delta t} a\frac{u_j^n - u_{j-1}^n}{\Delta x} 0 $$定义网格傅里叶数或CFL数 $\nu a\Delta t / \Delta x$方程可以整理成$$ u_j^{n1} (1-\nu)u_j^n \nu u_{j-1}^n $$这个形式特别直观$u_j^{n1}$ 本质上是当前步 $u_j^n$ 和上游 $u_{j-1}^n$ 的线性加权平均。系数 $1-\nu$ 和 $\nu$ 都是非负的这意味着新值不会跳出旧值范围格式能做到“有界”。物理意义也很清楚当 $a0$ 时信息是从左往右传的所以计算某个点的新值时应该看它的左边邻居这就是迎风的思想。如果我们反过来用空间中心差分配合时间前向$$ \frac{u_j^{n1} - u_j^n}{\Delta t} a\frac{u_{j1}^n - u_{j-1}^n}{2\Delta x} 0 $$这个格式的截断误差二阶比一阶迎风“看起来更高端”但它天然带锯齿型振荡的隐患。原因在第4章稳定性分析里会揭露这个格式对任何 $\Delta t$ 都是不稳定的换句话说根本不能用。只看截断误差会误导你必须配合稳定性分析才能判断格式好坏这是我读前3章最深的体会之一。3.3 修正方程截断误差背后的物理面孔第3章后半部分有一个非常精彩的概念修正方程。差分格式的截断误差并不仅仅是一个“数学误差”把截断误差项凑回原方程你会得到另一个偏微分方程它才是差分格式真正在求解的方程。以一阶迎风格式为例。它的修正方程长这样$$ \frac{\partial u}{\partial t} a\frac{\partial u}{\partial x} -\frac{a\Delta x}{2}\left(1-\nu\right)\frac{\partial^2 u}{\partial x^2} \cdots $$右边多了一项二阶导数项数值上相当于给原方程额外添加了一个“人工粘性”。这个人工粘性会让波形变钝、峰值降低这就是人们常说的数值耗散。如果你用的是中心差分格式修正方程里的主误差项就会变成三阶导数项对应一种色散效应波形不会变钝但会产生尾随的振荡这叫数值频散。这个认识对我后来看各种高级格式非常有帮助。比如TVD格式、WENO格式本质就是通过限制器控制这种隐含的人工粘性和频散让格式在间断附近不振荡、在光滑区域保持高精度。前3章能把修正方程理解透后面学这些高级格式就轻松很多。4. 第三根主线稳定性分析才是数值格式的“安全线”4.1 为什么要做稳定性分析真的只是怕发散吗很多人刚接触稳定性分析时以为不稳定的格式运行之后会直接计算溢出变成NaN只要避开就好。实际上即使是一个“线性稳定”的格式在真实计算时仍然可能因为边界处理、非线性效应等因素出现局部振荡只是振荡不会无限增长而已。稳定性分析的真正价值是给我们一个步长选择的量化依据告诉我们 $\Delta t$ 和 $\Delta x$ 之间必须满足什么约束格式才能保证误差不增长。数值计算里有两种误差离散误差和舍入误差。舍入误差总是存在如果格式在离散意义上是稳定的舍入误差的增长就会被控制在一个有限倍数内如果格式不稳定舍入误差会在时间推进中逐渐放大最终把真实解淹没。所以稳定性分析本质上回答的是舍入误差会不会在迭代过程中被“放大器”一样放大。4.2 von Neumann稳定性分析手把手推导前3章里的稳定性分析大多采用von Neumann方法这个方法推导起来并不复杂但每一步的物理含义都值得停下来想一想。把差分解写成单个傅里叶模的叠加$$ u_j^n \lambda^n e^{ikj\Delta x} $$这里 $k$ 是波数$\lambda$ 是增长因子。把它代入迎风格式 $u_j^{n1} (1-\nu)u_j^n \nu u_{j-1}^n$两边同时除以 $\lambda^n e^{ikj\Delta x}$得到$$ \lambda 1 - \nu(1 - e^{-ik\Delta x}) $$利用欧拉公式 $e^{-i\theta} \cos\theta - i\sin\theta$取增长因子的模方$$ |\lambda|^2 (1-\nu\nu\cos\theta)^2 (\nu\sin\theta)^2 $$化简之后是一个极其经典的结果$$ |\lambda|^2 1 - 2\nu(1-\nu)(1-\cos\theta) $$因为 $1-\cos\theta \geq 0$要让 $|\lambda|^2 \leq 1$ 对任意波数成立必须满足$$ 0 \le \nu \le 1 $$这就是一维线性对流方程迎风格式的稳定性条件通常写作$$ \nu \frac{a\Delta t}{\Delta x} \le 1 $$这一串推导下来CFL条件不再是“教材上背出来的数字”而是一个可以自己用手推出来的结论。我在推导之后自己写了个几行Python的小算例去验证 $\nu0.8$ 和 $\nu1.2$ 两种情况的差别。import numpy as np import matplotlib.pyplot as plt def advection_upwind(N100, cfl0.8, steps200): x np.linspace(0, 1, N) dx x[1] - x[0] u np.exp(-((x - 0.25) / 0.05) ** 2) for _ in range(steps): u_new u.copy() u_new[1:] u[1:] - cfl * (u[1:] - u[:-1]) u_new[0] 0.0 u u_new return x, u x, u advection_upwind(cfl0.8) print(u.max())实测下来$\nu0.8$ 时波形虽然会有轻微抹平但整体形状在$\nu1.2$ 时很快就出现高频锯齿振荡然后增长到失控。亲手跑一遍比看十遍公式都更有说服力。4.3 CFL条件的物理本质信息传播速度的“限速牌”稳定性分析算完以后教材通常会单独拿出CFL条件来解释一下物理含义。我这里用自己的话说一遍。假设某一层网格上的物理信息以速度 $a$ 传播在一个时间步 $\Delta t$ 内它能传播的距离是 $a\Delta t$。一个合理的数值格式必须保证信息在离散网格上的传播距离不能超过一个网格间距 $\Delta x$。如果 $\Delta t$ 太大$a\Delta t \Delta x$真实物理上信息已经跨过了好几层网格但格式只往相邻网格传了一个格子数值解的依赖域覆盖不了物理解的依赖域公式上就表现为不稳定。所以CFL条件本质上是一个“限速牌”。它把时间步长和空间步长绑定在一起$$ \Delta t \le \frac{\Delta x}{a} $$实际工程中我们通常不取在边界上。我自己做问题时一般把CFL数控制在0.5左右既保证了稳定性又不至于太保守导致计算量浪费。这里还要补充一点显式扩散格式也有类似的限制二阶导数项带来的限制是 $\alpha\Delta t/\Delta x^2 \le 0.5$因为扩散信息在网格间以二阶方式传播对时间步长的限制更苛刻。这也是为什么很多商业CFD软件对流项和扩散项会采用不同的离散策略扩散项经常用隐式格式就是为了绕开这个跟 $\Delta x^2$ 成反比的致命限制。5. 读前3章时踩过的坑以及我调整后的学法5.1 坑之一抄公式一时爽合上书本全忘光我第一次读前3章的时候最大的错误就是“抄公式”。看到推导过程觉得自己明白了拿起笔原样抄一遍抄完还有种莫名的成就感。结果到了要写程序的时候脑子里只剩下“CFL小于1”这一个结论其他细节全都模糊了。后来我换成了最小算例法。每次读完一个格式就手动取5个网格点把离散方程一列一列写出来。就拿迎风格式来说我会从 $j1$ 开始把 $u_1^{n1}$、$u_2^{n1}$ 直到 $u_5^{n1}$ 全部展开写成向量形式$$ \mathbf{u}^{n1} A\mathbf{u}^n $$然后观察矩阵 $A$ 的非零元素分布。这一下就看明白了两件事第一时间推进本质上是一个矩阵迭代第二迎风格式的信息传递方向在矩阵里一目了然。这个方法比单纯看公式有效得多可能多花15分钟却能让你真正从“看过”变成“会用”。5.2 坑之二迷信高阶格式忽略了格式的有界性读前3章时我有一段时间刻意追求高阶格式。原因很简单截断误差更低听起来更厉害。但后来被一组对比实验泼了冷水。中心差分在光滑区域确实精度高但一旦遇到间断或者大梯度区域就会在间断附近的上下游产生非物理振荡。更严重的是线性高阶格式几乎都会出现这个问题——Godunov定理说得非常清楚一个线性单调格式最高只能有一阶精度。之后我就把观念调整过来了实际开发中很少只用一个格式走天下。常用的做法是光滑区域用高阶格式间断附近自动切换成一阶迎风或带限制器的高阶格式这正是TVD/WENO类方法的核心思路。前3章教会我讨论格式精度不能脱离“解在哪个区域”这个前提。5.3 坑之三把“CFL小于1”当作万能口诀之前很长一段时间我遇到计算发散第一反应就是“把CFL调小”。这个做法有时候有效却治标不治本。前3章读完之后我才意识到一个发散的计算背后可能有完全不同的原因如果是双曲型问题先检查CFL数是否大于1如果是扩散主导问题先检查 $\alpha\Delta t/\Delta x^2$ 是否大于0.5如果都不是再检查边界条件有没有给对迎风方向是否跟波传播方向一致。我甚至整理了一张简易排查表每次程序发散了就按顺序核对。这个过程看起来繁琐实际上比盲目调小步长更省时间因为很多发散问题的根源根本不在时间步长上而是在离散方向或者边界处理上。前3章如果学透了你会有能力自己判断问题出在哪一环节而不是靠猜。5.4 接下来第4章以后的阅读计划读完前3章之后我的整体框架已经比较清晰了方程是描述物理的“神”差分是逼近物理的“形”稳定性是约束数值方法的“规矩”。第4章开始大概率会进入更具体的离散方法对比比如有限体积与有限差分的关系或者压力耦合算法。这一部分我打算采用同样的思路先理清这一章要解决什么问题再动手推公式然后写最小算例验证最后回到主线框架里给新的知识点找到位置。如果你也正在读《计算流体力学大串讲》或者类似的CFD教材我的建议是不要被前几章的数学量吓退。第一次读只需要抓住“方程线、离散线、稳定性线”这三条线其他细节都可以暂时跳过。读到后面突然卡住的时候再回过头补前3章的细节效率会高很多。最后说一个我自己的习惯每读完一章我会在笔记本上画一张一页纸的思维导图把这一章跟前面章节的联系用箭头标出来。这张图不追求完整只记录“这个知识点是干嘛用的”“它卡在哪个环节”。前3章读完后我的三张图从N-S方程一路指到了CFL条件形成了一条闭链路。现在做项目遇到数值不稳定问题我脑子里会直接浮现出这条链路对应的位置而不是像以前一样乱试参数。这就是我写这篇笔记最大的收获。