流体力学的基础 — 从纳维-斯托克斯方程到湍流模型的完整解释

分类：基础理论 / 流体力学 | 2026-03-25 | 阅读时间：约30分钟

流体力学是描述液体和气体流动所遵循物理规律的学科体系。在CAE中的流体分析（CFD：Computational Fluid Dynamics）广泛应用于汽车空气动力学、飞机翼型气流、热交换器设计、建筑物风荷载评估、电子设备冷却设计等现代工程的各个领域。本文从流体力学的数学基础开始，系统讲解实务中必需的湍流模型选择、CFD求解器的内部算法等内容。

CAE visualization for fluid dynamics - technical simulation diagram
流体力学

1. 流体的基本性质

连续体假设与流体质点

🧑‍🎓

老师，流体就是"水分子聚集在一起"对吧。CFD是不是在计算每一个分子？

🎓

完全同意。流体力学采用连续体假设。1 cm³空气中有约2.7×10¹⁹个分子，但从工程尺度来看，"连续分布的物质"这个近似就足够精确了。CFD网格单元即使只有1 mm³，其中也包含着数万亿分子。

不过在微流体装置（通道宽度数μm以下）或接近真空的极低压下，连续体假设会失效。此时需要通过Knudsen数 $Kn = \lambda/L$（$\lambda$：平均自由程, $L$：特征长度）检验。当 $Kn > 0.01$ 时需要特别注意。

🧑‍🎓

普通的汽车和飞机的流动可以用连续体吗？"流体质点"是什么，和分子不一样吗？

🎓

流体质点是"无穷小但包含足够多分子的微小体积元素"。比分子大，比流动结构小。在流体质点上可以定义位置、速度、压力、温度等状态量。这是连续体力学的出发点。

粘性的物理意义

粘性是"流体层间传递动量"的性质。对于牛顿流体，剪切应力 $\tau$ 与剪切速率 $du/dy$ 成正比（牛顿粘性律）：

$$\tau = \mu \frac{du}{dy}$$

其中 $\mu$ 是动力粘性系数（粘度）[Pa·s]，$u$ 是流速，$y$ 是距壁面的距离。运动粘性系数 $\nu = \mu/\rho$ [m²/s]（$\rho$：密度）在CFD中也经常用到。温度依赖性很强，液体中粘度随温度升高而降低，气体中粘度随温度升高而升高。

流体动力粘性系数 μ [Pa·s]（20℃）运动粘性系数 ν [m²/s]备注
空气1.81×10⁻⁵1.51×10⁻⁵温度变化影响大
水1.00×10⁻³1.00×10⁻⁶作为基准易于使用
轻油（20℃）3.0×10⁻³3.4×10⁻⁶温度依赖性很强
甘油（20℃）1.491.19×10⁻³粘度是水的1490倍
发动机油 SAE30（40℃）~0.1~1.1×10⁻⁴高温时快速下降
水银（20℃）1.53×10⁻³1.13×10⁻⁷密度大，所以 ν 小

可压缩性 vs 不可压缩性

🧑‍🎓

空气是可以压缩的吧。但我看到CFD解"不可压缩性"，什么时候能用？

🎓

实务中的经验是Mach数 $Ma < 0.3$。流速在音速（空气中约340 m/s）的30%以下时，密度变化不到1%，可视为不可压缩。时速100 km的车（约28 m/s，$Ma \approx 0.08$）完全没问题。反之战斗机巡航（$Ma \approx 0.85$）或超音速火箭必须考虑压缩性。

不可压缩模型的优点是能独立求解能量方程，密度是常数，计算简化很多。SIMPLE法这样的基于压力的求解器也是这样才能用的。

Mach数定义为 $Ma = V/a$（$V$：流速, $a = \sqrt{\gamma RT}$：音速）。$\gamma$ 是比热比（空气：1.4），$R$ 是气体常数（空气：287 J/(kg·K)），$T$ 是绝对温度。

Mach数范围流动名称密度变化推荐方法代表例
$Ma < 0.3$不可压缩区<2%不可压缩求解器（SIMPLE等）汽车、建筑、低速风扇
$0.3 < Ma < 0.8$亚音速可压缩2〜20%可压缩求解器或密度修正客机巡航、高速风扇
$0.8 < Ma < 1.2$跨音速很大、可能有冲击波必须用可压缩求解器翼型跨音速设计
$Ma > 1.2$超音速冲击波、膨胀波密度基求解器（Roe/AUSM）火箭、超音速飞机

牛顿流体 vs 非牛顿流体

满足 $\tau = \mu \, du/dy$ 的流体称为牛顿流体（水、空气、大多数低分子有机溶剂）。血液、聚合物溶液、润滑脂等许多实用流体表现出非牛顿性质。

模型本构式代表例CFD处理
牛顿$\tau = \mu \dot{\gamma}$水、空气、乙醇标准NS方程直接应用
幂律（Power Law）$\tau = K \dot{\gamma}^n$（$n<1$：假塑性，$n>1$：膨胀性）番茄酱、聚合物熔体代入有效粘度 $\mu = K\dot{\gamma}^{n-1}$
Bingham塑性体$\tau = \tau_0 + \mu_p \dot{\gamma}$（当 $\tau > \tau_0$ 时流动）润滑脂、混凝土、泥流考虑屈服应力 $\tau_0$ 的模型
Herschel-Bulkley$\tau = \tau_0 + K\dot{\gamma}^n$血液、发胶幂律+屈服应力的组合
Carreau-Yasuda$\mu = \mu_\infty + (\mu_0-\mu_\infty)[1+(\lambda\dot{\gamma})^a]^{(n-1)/a}$血液、牙膏Ansys Fluent等直接输入可用

2. 连续方程（质量守恒）

积分形与微分形的推导

对任意检查体积 $V$ 应用质量守恒。检查体积内质量变化率 = 流入质量通量 − 流出质量通量：

$$\frac{\partial}{\partial t}\int_V \rho \, dV + \oint_S \rho \mathbf{v} \cdot \hat{n} \, dS = 0$$

应用Gauss发散定理（$\oint_S \mathbf{F}\cdot\hat{n}\,dS = \int_V \nabla\cdot\mathbf{F}\,dV$）得到微分形式：

$$\frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{v}) = 0$$

展开后：$\frac{\partial \rho}{\partial t} + \rho\nabla\cdot\mathbf{v} + \mathbf{v}\cdot\nabla\rho = 0$，用物质导数表示为：

$$\frac{D\rho}{Dt} + \rho\nabla\cdot\mathbf{v} = 0$$

不可压缩流体（$\rho = \text{常数}$，即 $D\rho/Dt = 0$）时：

$$\nabla \cdot \mathbf{v} = 0 \quad \Leftrightarrow \quad \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} + \frac{\partial w}{\partial z} = 0$$
🧑‍🎓

$\nabla \cdot \mathbf{v} = 0$ 经常出现。具体表示什么？

🎓

粗俗说法："任何地方都不积累流体，任何地方都不冒出来"。想象一根断面积变化的管子。变细时流速加快，变粗时流速变慢。这个关系用 $A_1 v_1 = A_2 v_2$（连续方程的积分形式）表示，微分形式就是 $\nabla \cdot \mathbf{v} = 0$。

在CFD中，如果这个条件不满足，就会"质量凭空消失或凭空产生"导致出错。continuity residual（连续残差）很大时，通常是网格或边界条件设置的问题。

管道截面变化与流速的关系（实例）

直径 $D_1 = 50$ mm的管道缩小到 $D_2 = 25$ mm时，入口流速 $v_1 = 2$ m/s，出口流速为：

$$v_2 = v_1 \frac{A_1}{A_2} = v_1 \left(\frac{D_1}{D_2}\right)^2 = 2 \times \left(\frac{50}{25}\right)^2 = 8 \text{ m/s}$$

截面积减小到1/4时，流速增加4倍。结合伯努利定理（$p + \frac{1}{2}\rho v^2 + \rho g z = \text{常数}$），可求出压力下降 $\Delta p = \frac{1}{2}\rho(v_2^2 - v_1^2) = \frac{1}{2} \times 1000 \times (64-4) = 30\,000$ Pa ≈ 0.3 atm。

物质导数的概念

🧑‍🎓

NS方程里有 $D/Dt$ 这个记号。和普通偏导数 $\partial/\partial t$ 有什么区别？

🎓

这是流体力学初学者最容易卡的地方。$\partial/\partial t$ 是"看固定位置随时间的变化"（欧拉描述）。$D/Dt$ 是"坐在流体质点上随时间的变化"（拉格朗日描述）。

比如河流水温：$\partial T/\partial t$ 是"在桥上测量同一地点的水温变化"，$DT/Dt$ 是"坐橡皮艇跟着漂流感受水温变化"。两者的关系是 $\frac{D}{Dt} = \frac{\partial}{\partial t} + (\mathbf{v} \cdot \nabla)$，其中 $(\mathbf{v} \cdot \nabla)$ 项表示"对流"。

🧑‍🎓

那CFD代码里如何离散化 $D/Dt$ 呢？

🎓

在FVM中，$\partial/\partial t$ 部分用时间差分（向后差分、Crank-Nicolson等）离散化，对流项 $(\mathbf{v} \cdot \nabla)\phi$ 作为单元边界面的通量求值。如何评估这个对流通量就是"对流格式"，有一阶迎风、QUICK、高阶带限制器的格式等多种。后面会详细说明。

$$\frac{D}{Dt} = \frac{\partial}{\partial t} + u\frac{\partial}{\partial x} + v\frac{\partial}{\partial y} + w\frac{\partial}{\partial z} = \frac{\partial}{\partial t} + (\mathbf{v} \cdot \nabla)$$

NS方程的推导与物理意义

对流体质点应用牛顿第二定律（$F = ma$）。作用在流体质点上的力包括：(1) 压力梯度力、(2) 粘性力（剪切应力的发散）、(3) 体积力（重力等）。用牛顿流体的应力张量：

$$\sigma_{ij} = -p\delta_{ij} + \mu\left(\frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i}\right) - \frac{2}{3}\mu\delta_{ij}\frac{\partial u_k}{\partial x_k}$$

代入动量守恒式，假设不可压缩（$\nabla\cdot\mathbf{v}=0$）和常数粘度，得到纳维-斯托克斯方程：

$$\rho \frac{D\mathbf{v}}{Dt} = \rho\left(\frac{\partial \mathbf{v}}{\partial t} + (\mathbf{v}\cdot\nabla)\mathbf{v}\right) = -\nabla p + \mu \nabla^2 \mathbf{v} + \rho \mathbf{g}$$

展开的x分量形式：

$$\rho \left(\frac{\partial u}{\partial t} + u\frac{\partial u}{\partial x} + v\frac{\partial u}{\partial y} + w\frac{\partial u}{\partial z}\right) = -\frac{\partial p}{\partial x} + \mu\left(\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} + \frac{\partial^2 u}{\partial z^2}\right) + \rho g_x$$

整理各项的物理意义：

项公式物理意义CFD注意事项
非定常项$\rho \partial\mathbf{v}/\partial t$流体质点的动量时间变化定常分析中省略为零
惯性项（对流项）$\rho(\mathbf{v} \cdot \nabla)\mathbf{v}$流动引起的动量传输非线性！数值不稳定的主要原因
压力梯度项$-\nabla p$压力差引起的驱动力不可压缩时需单独求解压力方程
粘性项$\mu \nabla^2 \mathbf{v}$粘性引起的动量扩散（平滑化）高Re时较小，注意数值扩散
体积力$\rho \mathbf{g}$重力、浮力（热对流中重要）Boussinesq近似或完全可压缩处理
🧑‍🎓

强调惯性项是"非线性"，具体什么问题？

🎓

非线性就是有 $\mathbf{v}$ 与 $\mathbf{v}$ 的乘积。线性方程可以用叠加原理，很容易求解。但非线性就不行了。这个惯性项正是产生湍流、涡流、Karman涡街这些复杂现象的根本原因。CFD需要迭代计算（欠松弛）就是因为非线性对流项的处理。

反过来，Re数非常小（$Re \ll 1$）的Stokes流中可以忽略惯性项，变成线性方程"Stokes方程"，能得到精确解。微流路生物流体分析和滑动轴承膜等用这个近似。

无量纲化与雷诺数的推导

用代表长度 $L$、代表速度 $U$、动压 $\rho U^2$ 进行无量纲化（带星号量为无量纲）：

$$\frac{\partial \mathbf{v}^*}{\partial t^*} + (\mathbf{v}^* \cdot \nabla^*)\mathbf{v}^* = -\nabla^* p^* + \frac{1}{Re}\nabla^{*2}\mathbf{v}^*$$

这里出现了唯一的无量纲参数——雷诺数：

$$Re = \frac{\rho U L}{\mu} = \frac{UL}{\nu} = \frac{\text{惯性力}}{\text{粘性力}}$$

无量纲化的重要启示：若两个流动的 Re 数相同且几何相似，则流动模式也相似。这是风洞试验用缩小模型的理论基础。例如实际模型的1/10缩小模型用10倍流速就能得到相同的Re数，压力系数等无量纲量完全一致。

4. 能量方程

在可压缩流和流热耦合问题中，需要将能量方程与连续式、动量式联立求解。用全焓 $H = h + V^2/2$（$h$：比焓, $V = |\mathbf{v}|$）的守恒形式：

$$\frac{\partial(\rho H)}{\partial t} + \nabla \cdot (\rho H \mathbf{v}) = \frac{\partial p}{\partial t} + \nabla \cdot (k \nabla T) + \nabla \cdot (\boldsymbol{\tau} \cdot \mathbf{v}) + \dot{Q}$$

在不可压缩、常物性、忽略粘性耗散的情况下，简化为纯温度的对流扩散方程：

$$\rho c_p \left(\frac{\partial T}{\partial t} + \mathbf{v} \cdot \nabla T\right) = k \nabla^2 T + \dot{Q}$$

其中 $c_p$ 是定压比热 [J/(kg·K)]，$k$ 是热传导系数 [W/(m·K)]。此方程中温度场单向依赖于速度场（速度→温度），但在自然对流或Boussinesq近似中，密度对温度的依赖会产生双向耦合。

🧑‍🎓

能量方程是CFD必须解的吗？不解也可以吗？

🎓

实务中经常做"等温分析"把温度方程关掉。比如汽车车身空气动力学，如果只想知道风压、阻力系数，用不可压缩等温就够了。但电子设备散热（CPU排热）、发动机燃烧室、超音速流动必须求解能量方程。

要注意的是Ansys Fluent等软件中"Energy: OFF"忘记关掉的话，温度会变成乱七八糟的值但计算还能跑。温度与密度有函数关系（理想气体 $\rho = p/(RT)$）时更要当心。汽车发动机舱热流动这种"高温气体+强制对流"的问题，必须把能量方程打开。

粘性耗散项的重要性

在高速流和高粘流体中，粘性耗散（$\Phi = \boldsymbol{\tau}:\nabla\mathbf{v}$）变得重要。

$$\Phi = 2\mu\left[e_{xx}^2 + e_{yy}^2 + e_{zz}^2 + e_{xy}^2 + e_{yz}^2 + e_{zx}^2\right] - \frac{2}{3}\mu(\nabla\cdot\mathbf{v})^2$$

其中 $e_{ij}$ 是应变速率张量的分量。高粘高分子聚合物挤出成形或高Mach数超音速流中粘性耗散导致的温度上升（气动加热）不能忽视。大气圈再入航天器之所以高温，正是这个现象。

5. 雷诺数与流动尺度

层流、过渡、湍流的判定

流动状态圆管内（水力直径 $D_h$）平板上（距离 $x$）特征
层流$Re_D < 2300$$Re_x < 5 \times 10^5$流线整齐，压力损失∝流速一次方
过渡流$2300 \sim 4000$$5\times10^5 \sim 10^6$不稳定，难以分析预测
湍流$Re_D > 4000$$Re_x > 10^6$混乱涡运动，压力损失∝流速1.75〜2次方
🧑‍🎓

实际的机械设备一般是什么Re数范围呢？有点摸不着边。

🎓

举例说明：

  • 毛细血管（$D \approx 10\,\mu\text{m}$）: $Re \approx 0.001$ — 完全层流，惯性可忽略的Stokes流
  • 淋浴喷头喷孔: $Re \approx 500\sim2000$ — 层流~过渡区
  • 家庭自来水管（φ13 mm, 2 L/min 流量）: $Re \approx 3000$ — 过渡~湍流
  • 汽车车身（100 km/h, 特征长度 4 m）: $Re \approx 1.1\times10^7$ — 完全湍流
  • 大型客机翼（巡航速度, 弦长 7 m）: $Re \approx 5\times10^8$ — 极高Re湍流
  • 泵的叶轮: $Re \approx 10^6\sim10^7$ — 湍流

CFD用得最多的就是湍流区，所以湍流模型选择至关重要。

边界层发展

平板上边界层厚度 $\delta(x)$ 随流向距离 $x$ 的增长规律：

$$\delta_{\text{层流}} \approx \frac{5.0\,x}{\sqrt{Re_x}}, \qquad \delta_{\text{湍流}} \approx \frac{0.37\,x}{Re_x^{1/5}}$$

壁面摩擦系数 $C_f$（局部）：

$$C_{f,\text{层流}} = \frac{0.664}{\sqrt{Re_x}}, \qquad C_{f,\text{湍流}} \approx \frac{0.0594}{Re_x^{1/5}}$$

湍流边界层生长更快，壁面剪切应力更大，摩擦阻力更大。但湍流边界层对逆压力梯度（adverse pressure gradient）的抗分离能力强。高尔夫球表面的小凹陷就是为了主动形成湍流边界层，让分离点向后移，降低阻力。

6. 湍流基础与湍流模型

雷诺平均（RANS）与湍流脉动

湍流中速度随时间脉动。将变量分解为时间平均值与脉动分量（雷诺分解）：

$$\mathbf{v} = \overline{\mathbf{v}} + \mathbf{v}', \quad p = \bar{p} + p' \quad (\text{其中}\; \overline{\mathbf{v}'} = 0,\; \overline{p'} = 0)$$

代入NS方程并取时间平均，得到RANS方程：

$$\rho \left(\frac{\partial \overline{u}_i}{\partial t} + \overline{u}_j\frac{\partial \overline{u}_i}{\partial x_j}\right) = -\frac{\partial \bar{p}}{\partial x_i} + \frac{\partial}{\partial x_j}\left[\mu\frac{\partial \overline{u}_i}{\partial x_j} - \rho\overline{u'_i u'_j}\right]$$

$-\rho\overline{u'_i u'_j}$ 称为雷诺应力张量。对称张量，独立分量6个（三维）。这成为新的未知量，需要建立额外的"湍流模型"来闭合方程组（闭合问题）。

🧑‍🎓

"闭合问题"是什么意思？未知数多了，没办法吗？

🎓

正好。RANS方程出现了雷诺应力这个新的未知量（三维6个分量），但方程个数没增加。这叫"方程未闭合"。湍流模型的作用就是用已知的平均量近似雷诺应力。

最简单的近似是涡粘性假设（Boussinesq近似）：$-\rho\overline{u'_i u'_j} \approx \mu_t\left(\partial\overline{u}_i/\partial x_j + \partial\overline{u}_j/\partial x_i\right) - \frac{2}{3}\rho k \delta_{ij}$，加入湍流粘性系数 $\mu_t$，用额外的输运方程求解。$k = \frac{1}{2}\overline{u'_i u'_i}$ 是湍流动能。

湍流能量级联与Kolmogorov尺度

按Kolmogorov理论，湍流从大尺度（注入能量尺度 $L$）逐步向小尺度传递能量（能量级联）。

$$\eta = \left(\frac{\nu^3}{\varepsilon}\right)^{1/4}, \quad \tau_\eta = \left(\frac{\nu}{\varepsilon}\right)^{1/2}, \quad u_\eta = (\nu\varepsilon)^{1/4}$$

其中 $\varepsilon$ 是湍流能量耗散率 [m²/s³]。高Re流中尺度比 $L/\eta \sim Re^{3/4}$，完全求解所有尺度（DNS）需要计算量 $\propto Re^{9/4}$，代价巨大。

方法求解尺度计算代价适用Re数
DNS（直接数值模拟）直到Kolmogorov尺度 $\eta$$\propto Re^{9/4}$（极高）$Re < 10^4$（基础研究）
LES（大涡模拟）惯性小子区以上的尺度$\propto Re^{13/7}$（高）$Re \sim 10^5\sim10^6$（详细分析）
RANS（雷诺平均）所有尺度平均化低~中任意（实用）
DES（分离涡模拟）壁近处RANS + 自由区LES中~高$Re \sim 10^6\sim10^7$

主要湍流模型比较

模型输运方程擅长的流动不擅长的流动应用场景
Spalart-Allmaras（SA）1个（$\tilde{\nu}$）翼面附着边界层分离、逆压梯度飞机外部空气动力
标准 k-ε2个（$k$, $\varepsilon$）充分发展管内流、喷流壁近处、分离建筑通风、烟羽
RNG k-ε2个（RNG推导）有涡、曲率的流动强分离混合罐、室内通风
可实现 k-ε2个（可实现条件）喷流、剪切流、通道流强旋转涡轮机初步检查
标准 k-ω2个（$k$, $\omega$）壁近处、逆压梯度自由流依赖性（外部流问题大）边界层详细解析
k-ω SST (Menter)2个（混合函数）壁近+自由流兼顾，分离流强浮力效应（自然对流）最通用。汽车、翼型、涡轮
RSM（雷诺应力）7个（全应力分量）强各向异性、旋流收敛慢、复杂旋转燃烧器、旋风分离器
LES（Smagorinsky等）非定常+SGS模型大涡详细结构、声学壁近处高Re（成本过高）噪声、燃烧详细过程
DES/DDESSA or SST 基础大规模分离（车尾流等）收敛判断难汽车尾流、钝头体
🧑‍🎓

为什么k-ω SST被称"最通用"？k-ε不行吗？

🎓

k-ε的弱点在壁面附近。要精确模拟，需要依赖壁函数（wall function），y+得在30〜300范围才行。k-ω则擅长壁面，但自由流（远离壁面）对初始湍流条件过度敏感。

SST（Shear Stress Transport）用混合函数"壁面近k-ω，自由流ε-like"，两者长处结合，还加了湍流剪切应力过估修正。自1994年Menter提出后，航空、汽车、涡轮机就都用这个。"不知道选什么"的第一选择就是SST。

y⁺（y-plus）的含义与壁面网格设计

$y^+$ 是用"粘性子层尺度"对壁面距离进行无量纲化：

$$y^+ = \frac{y \cdot u_\tau}{\nu}, \qquad u_\tau = \sqrt{\frac{\tau_w}{\rho}} \text{（摩擦速度）}$$

壁面附近速度分布（对数律）：

$$u^+ = \frac{u}{u_\tau} = \begin{cases} y^+ & (y^+ < 5,\; \text{粘性子层}) \\ \frac{1}{\kappa}\ln(y^+) + B & (y^+ > 30,\; \text{对数律层, } \kappa \approx 0.41, B \approx 5.2) \end{cases}$$
湍流模型/处理方式推荐 y⁺壁面处理备注
k-ε（标准壁函数）30〜300壁函数补间粗网格可用。计算成本低
k-ε（增强壁处理）y⁺ ≈ 1 也可自动切换Ansys Fluent推荐设置
k-ω SST（低Re求解）y⁺ ≈ 1粘性子层求解精度优先。计算代价高
k-ω SST（壁函数）y⁺ = 30〜300壁函数大规模分析的折衷方案
LES / DESy⁺ ≈ 1〜5壁面求解网格极细，成本高

若y⁺ 在推荐范围外，壁面剪切应力、传热系数会产生大误差。分析前必须用 y⁺计算器 检查第一层网格高度。

7. 可压缩流与冲击波

等熵流与驻点条件

无摩擦、无热交换的等熵流中，完全气体的状态量仅由Mach数决定：

$$\frac{T_0}{T} = 1 + \frac{\gamma-1}{2}Ma^2$$
$$\frac{p_0}{p} = \left(1 + \frac{\gamma-1}{2}Ma^2\right)^{\gamma/(\gamma-1)}, \qquad \frac{\rho_0}{\rho} = \left(1 + \frac{\gamma-1}{2}Ma^2\right)^{1/(\gamma-1)}$$

下标0表示驻点（全量）条件，无下标表示静量。$\gamma = 1.4$（空气），$Ma = 1$ 时喉部条件：$T^*/T_0 = 2/(\gamma+1) = 0.833$，$p^*/p_0 = 0.528$。

德·拉瓦尔喷管与面积比

要在喷管内产生超音速，需要先收缩后膨胀的形状（德·拉瓦尔喷管、收缩膨胀喷管）。喉部面积 $A^*$ 与某点面积 $A$ 的比与Mach数的关系：

$$\frac{A}{A^*} = \frac{1}{Ma}\left[\frac{2}{\gamma+1}\left(1 + \frac{\gamma-1}{2}Ma^2\right)\right]^{(\gamma+1)/[2(\gamma-1)]}$$
🧑‍🎓

火箭喷管就是这个形状吧。但为什么截面扩大就超音速了？亚音速里"细就快"，怎么反了？

🎓

"音速超过后连续方程反向工作"的现象。亚音速 $d\rho/\rho \ll dA/A$，管细了流速快（体积守恒≈断面守恒）。超音速后 $d\rho/\rho$ 比 $dA/A$ 大得多——密度下降快，截面虽扩大，连续性让流速必须更快。

喉部达到 $Ma = 1$ 时出现"堵塞（choke）"：上游压力再提也通过喉部的质量流量不增。火箭发动机喇叭形喷管正是利用这原理把燃烧气加速到超音速，获得最大推力。

垂直冲击波关系式（Rankine-Hugoniot）

垂直冲击波前后（下标1：上游, 2：下游）的关系式：

$$Ma_2^2 = \frac{(\gamma-1)Ma_1^2 + 2}{2\gamma Ma_1^2 - (\gamma-1)}$$
$$\frac{p_2}{p_1} = \frac{2\gamma Ma_1^2 - (\gamma-1)}{\gamma+1}, \qquad \frac{\rho_2}{\rho_1} = \frac{(\gamma+1)Ma_1^2}{(\gamma-1)Ma_1^2+2}$$
$$\frac{T_2}{T_1} = \frac{[2\gamma Ma_1^2 - (\gamma-1)][(\gamma-1)Ma_1^2 + 2]}{(\gamma+1)^2 Ma_1^2}$$

冲击波是几个平均自由程厚度的薄不连续面，压力、密度、温度陡增，熵增加（不可逆过程）。CFD中捕捉冲击波需要控制数值扩散，密度基求解器（Roe法、AUSM法、JST法等）加TVD格式必不可少。

$Ma_1$（上游）$Ma_2$（下游）$p_2/p_1$$T_2/T_1$$\rho_2/\rho_1$
1.01.0001.0001.0001.000
1.50.7012.4581.3201.862
2.00.5774.5001.6882.667
3.00.47510.332.6793.857
5.00.41529.005.8005.000

8. 计算流体力学（CFD）概论

有限体积法（FVM）的离散化

CFD最常用的是有限体积法（FVM：Finite Volume Method）。计算域分割成单元（控制体积），对每个单元应用保存律的积分形式。在单元中心 $P$ 与相邻单元 $N$ 的边界面处评估通量。

$$\frac{\partial}{\partial t}(\rho\phi V_P) + \sum_{\text{faces}} \dot{m}_f \phi_f = \sum_{\text{faces}} \Gamma_\phi A_f \left(\frac{\partial \phi}{\partial n}\right)_f + S_\phi V_P$$

$\phi$ 是守恒变量（速度分量、温度、$k$、$\varepsilon$ 等），$\dot{m}_f$ 是边界面的质量通量，$\Gamma_\phi$ 是扩散系数，$A_f$ 是面积，$S_\phi$ 是源项。FVM的特点是严格满足每个单元的守恒律（守恒性好）。

压力速度耦合：SIMPLE法

🧑‍🎓

不可压缩的NS方程没有压力方程吧？CFD怎么求压力？

🎓

问得好！不可压缩中压力是"满足连续方程 $\nabla \cdot \mathbf{v} = 0$ 的约束"而确定。这就是压力速度耦合问题。SIMPLE（半隐式压力关联方程法）的流程：

  1. 用假定的压力场 $p^*$ 解动量方程，得到假的速度 $\mathbf{v}^*$
  2. $\mathbf{v}^*$ 不满足连续方程，所以建立压力修正量 $p'$ 的Poisson方程求解
  3. $p = p^* + \alpha_p p'$、$\mathbf{v} = \mathbf{v}^* + \mathbf{v}'$ 用欠松弛系数 $\alpha_p$ 更新
  4. 重复直到残差足够小

进阶有SIMPLEC（系数更精确）、非定常用PISO（Pressure-Implicit Splitting Operators）。OpenFOAM的非定常分析常用PIMPLE（SIMPLE+PISO混合），效率好。

对流格式的比较

格式精度阶稳定性数值扩散振荡用途
一阶迎风（UDS/UD）一阶◎（无条件稳定）大无收敛检查、初值、稳健性优先
中心差分（CDS）二阶△（高Pe数不稳定）无有LES、DNS（数值扩散最小需要）
线性迎风（LUDS）二阶○中小定常RANS标准选择
QUICK（三阶迎风加权）三阶○小略有精度优先的RANS
TVD格式（MUSCL+限制器）二~三阶○小抑制冲击波、密度不连续流
WENO（高阶权重本质无振荡）五阶以上○极小无超音速、爆炸波、精密LES

实务中常见做法是先用一阶迎风收敛，再切换为二阶精度格式提高精度的"两阶段法"。直接用高阶格式容易发散。

收敛判定与残差读法

🧑‍🎓

OpenFOAM计算时会输出各种"残差"。怎样才算"收敛"？残差要多小？

🎓

一般目标是：

  • 速度分量 ($U_x, U_y, U_z$): $10^{-4}$ 以下最好，$10^{-3}$ 时工程上有时可接受
  • 连续残差 (continuity): $10^{-4}$ 以下。大的话质量守恒破坏
  • 能量（温度）: $10^{-7}$ 以下（值小易收敛）
  • 湍流变量 ($k, \omega, \varepsilon$): $10^{-5}$ 以下

但只看残差是危险的。监控关心的物理量（阻力系数 $C_D$、压力损失 $\Delta p$、热流 $q$）是否不再变化。残差下降但 $C_D$ 还飘忽不定，说明没真正收敛。定量监控点设置必须做，这样才能判断答案稳定了。

9. CFD实务注意事项

边界条件设置错误与对策

CFD分析中最常见的失败就是边界条件误设。虽然计算能跑，但物理上无意义的边界条件会导致完全错误的结果。

边界常见错误正确设置与检查要点
入口（Inlet）湍流强度和长度尺度用默认值（TI=5%）TI = 1〜5%（内部流路）或0.1〜1%（外部流）。长度尺度约为水力直径的0.07倍
出口（Outlet）预期有反流处设压力出口，反流条件未设反流可能的地方设backflow湍流条件。计算域延伸到出口后方也有效
壁面（Wall）旋转壁的速度错误设为零（涡轮机械）使用Moving Wall / MRF（多参考框架）/ Sliding Mesh
对称面（Symmetry）对非对称流（Karman涡街）用对称条件Karman涡、非对称分离不能用对称。半模型也须用非定常
热边界外壁设绝热（adiabatic）做热分析。全温与静温混淆依据测量值或规格明确设热通量或对流传热系数+外气温
周期边界周期长度与物理周期不符压力降指定型（Fluent的Periodic BC）与流量指定型的区分

网格品质与数值扩散

🧑‍🎓

"数值扩散"经常听到，但它具体是什么问题？和物理扩散怎么区分？

🎓

数值扩散是"数值误差表现得像物理扩散"的现象。一阶迎风格式的截断误差展开，实际上加进了伪扩散系数 $\nu_{\text{num}} \approx U\Delta x/2$。这项大的话，本来尖锐的温度梯度或涡就会"模糊"。特别是网格不对齐流向时更坏。

辨别方法是加细网格看结果是否变化。粗网格和细网格差异大，说明数值扩散在主导。这叫网格无关性验证（Grid Independence Study）。实务中至少比较两个网格级别（粗、细）。

品质指标良好范围问题值影响
纵横比（Aspect Ratio）< 100（非边界层）> 1000收敛差、精度低
斜度（Skewness）< 0.5> 0.85离散误差增大、易发散
正交性（Non-orthogonality）< 40° (OpenFOAM定义)> 70°压力速度耦合误差大
体积比（Volume Change Ratio）< 10> 100体积剧变导致局部误差

物性值的数据源与注意

物性值错误是"计算能跑但答案全错"的常见原因。遵守以下原则：

🧑‍🎓

最后问一个：做CFD最容易犯的三个错误是什么？

🎓

①不验证就直接切细网格：应先用粗网格跑通，检查"物理是否合理"，再加细网格。一开始就10倍网格，答案不对时修改成本巨大。

②湍流强度用默认值（5%）不改：入口条件的湍流强度（Turbulent Intensity）对结果影响大，却经常被忽视。应依实验或规范设定。

③只看残差判断收敛：监控感兴趣的物理量（$C_D$、$\Delta p$、传热系数），不能只看残差。另外用基准（