空气动力学基础与CFD工程实践：从N-S方程到Python算例

简介：《空气动力学基础（北航精品课程）》PDF 是北京航空航天大学刘沛清老师主讲的课程讲义，面向航空航天专业学生、流体力学初学者及相关工程人员，系统梳理空气动力学核心知识体系。内容从绪论出发，覆盖流体基本属性、流体静力学与运动学、不可压缩无粘流体平面势流、粘性流体动力学基础、边界层理论及可压缩高速流动基础等章节，并配有风洞、机翼绕流等工程实例图示，便于理解抽象概念。资源包内仅 1 个 PDF 文件，大小约 19.65 MB，图文排版清晰，适合作为课堂笔记、考研复习或自主入门的学习资料。已有 669 人学习/下载，说明内容具有一定参考价值。读者可凭此快速搭建空气动力学理论框架，掌握飞行器绕流、边界层分离、量纲分析等关键知识点，为后续进阶研究与工程应用打下基础。

1. 从一份课程 PDF 聊起：为什么做 CFD 的人也要回头啃空气动力学

手头这份《空气动力学基础(北航精品课程)-刘沛清.pdf》，在航空航天圈子里几乎是人尽皆知的教学材料。很多入行 CFD（计算流体力学）的工程师，最早接触的流动物理概念不是来自商业软件的用户手册，而是来自这类课程讲义里对“连续性方程怎么推导”“边界层为什么分离”的严谨解释。如果你只把流体力学当工具用法学，不碰控制方程，那你很难解释为什么同一个网格在攻角 12 度时结果突然发散，或者为什么 SST k-omega 模型在逆压梯度区给出的分离点总比实验晚。这篇博客就以这份课程 PDF 的核心知识体系为骨架，从方程讲到数值实现，再落到用 Python 做几个典型算例——不依赖任何商业软件，也能把升力系数、边界层厚度这些关键量估算出来。

这适合两类人：一类是刚接手 CFD 仿真任务、需要补流体理论短板的工程师；另一类是写求解器或做网格工具的开发人员，需要理清对方口中“压力修正”“涡量”“激波捕捉”到底指什么。文章中的公式会控制在手算和写代码都能用的程度，不会出现一页纸的推导，但关键的物理假设和适用边界会讲清楚。

2. 控制方程与无量纲参数：先把 N-S 方程“读薄”

2.1 从连续介质假设到 N-S 方程的四个物理项

空气动力学的起点是纳维-斯托克斯方程（N-S 方程）。它由质量守恒（连续性方程）、动量守恒（三个方向的动量方程）和能量守恒组成。连续介质假设要求特征长度远大于分子自由程，这在地面到平流层范围内的空气流动基本都成立，所以不需要碰玻尔兹曼方程，直接对标量输运方程即可。

对不可压缩流动，连续性方程简化为速度散度为零，动量方程写成如下形式：

[ \frac{\partial \vec{V}}{\partial t} + (\vec{V} \cdot abla)\vec{V} = -\frac{1}{\rho} abla p + u abla^2 \vec{V} ]

方程里四项从左到右分别是：当地加速度（非定常项）、对流加速度（惯性项）、压力梯度项、粘性扩散项。做计算的人最容易忽略的是第二项：它是非线性的，也是造成流动不稳定的根源。CFD 里的所谓“数值耗散”本质上就是在处理这一项时引入的人工耗散，它掩盖了真实的物理粘性。

2.2 雷诺数、马赫数和克努森数的实际判断准则

无量纲数的意义在于：只要两个流动的无量纲控制参数一致，即使尺度不同，流场也是相似的。空气动力学里必须心里有数的三个数是：

以标准海平面条件为例，空气密度约为 1.225 kg/m³，动力粘性系数约为 1.789e-5 kg/(m·s)。一架弦长 1 米的机翼在 50 m/s 下飞行时，雷诺数大约为 3.4e6；如果换成高空无人机，密度降到 0.4 kg/m³ 左右，同样速度下雷诺数会掉一个量级。这直接影响你选湍流模型还是转捩模型。

2.3 伯努利方程的正确打开方式：不是“流速大压力小”

课程里反复强调的一句话：伯努利方程不是独立原理，而是 N-S 方程沿流线积分的结果。它的使用条件有三个：定常、无粘、不可压缩（或等熵可压缩）。如果你把伯努利方程用到机翼上表面“流速大所以吸力大”，必须加一个前提——上表面气流在到达压力最低点前是等熵加速的，一旦出现激波或强逆压梯度导致分离，伯努利关系就不再成立。这也是很多入门者拿“流速大压力小”解释激波诱导分离时翻车的根本原因。

3. 边界层理论与粘性阻力估算：用 Python 算到工程精度

3.1 边界层位移厚度和动量厚度的物理意义

普朗特的边界层理论把流场分成边界层内和边界层外两部分：层外无粘，层内粘性不可忽略。边界层内速度从壁面 0 增长到外缘速度 (U_e)，但由于速度亏损，实际流量比无粘假设少，等价于把壁面向外推了 (\delta^*)（位移厚度）；动量通量亏损等价于引入 (\theta)（动量厚度）。

平板层流边界层的 Blasius 解给出：

[ \delta^* = 1.7208 \sqrt{\frac{ u x}{U_\infty}}, \quad \theta = 0.664 \sqrt{\frac{ u x}{U_\infty}} ]

把这个公式转成代码做参数扫描，马上就能看出尺度关系。

3.2 给平板边界层算厚度和摩擦阻力系数的 Python 脚本

下面是一段直接可运行的估算脚本：

import numpy as np import matplotlib.pyplot as plt # 物理参数：标准海平面空气 rho = 1.225 # 密度 kg/m^3 mu = 1.789e-5 # 动力粘性 Pa·s U_inf = 50.0 # 来流速度 m/s L = 1.0 # 板长 m # 运动粘性 nu = mu / rho Re_L = U_inf * L / nu # 沿板位置 x = np.linspace(0.001, L, 500) # 层流 Blasius：位移厚度与动量厚度 delta_star = 1.7208 * np.sqrt(nu * x / U_inf) theta = 0.664 * np.sqrt(nu * x / U_inf) # 当地摩擦系数 Cf = 0.664 / sqrt(Re_x) Re_x = U_inf * x / nu Cf = 0.664 / np.sqrt(Re_x) # 全板平均摩擦阻力系数（层流） Cf_avg_lam = 1.328 / np.sqrt(Re_L) print(f"Re_L = {Re_L:.2e}") print(f"层流全板平均 Cf = {Cf_avg_lam:.5f}") # 平板单位展长总摩擦阻力（两侧） F_friction = 0.5 * rho * U_inf**2 * L * Cf_avg_lam * 2 print(f"单位展长摩擦阻力 = {F_friction:.3f} N/m") # 作图 fig, ax1 = plt.subplots() ax1.plot(x, delta_star*1000, label='位移厚度 (mm)') ax1.plot(x, theta*1000, label='动量厚度 (mm)') ax1.set_xlabel('x (m)'); ax1.set_ylabel('厚度 (mm)') ax1.legend(); ax1.grid(True) plt.show()

这段代码干了几件事：先由运动粘性和特征长度算雷诺数；再沿板面离散位置计算每一点的边界层厚度和当地摩擦系数；最后按 Blasius 全板平均公式算总阻力。打印输出里Re_L用来判断全板是否保持层流，工程经验是超过 3e5 到 5e5 就可能转捩，这时层流公式就不适用了，要换湍流公式。

3.3 转捩判据与湍流边界层的工程估算

如果雷诺数超过临界值，层流边界层会失稳转捩为湍流。工程上常用的湍流平板摩擦阻力系数为：

[ C_f = \frac{0.074}{Re_L^{1/5}} ]

作为对比，同样条件下单位宽度平板的湍流阻力可能是层流的数倍。你可以在代码里把Cf_avg_lam换成0.074 / Re_L**(1/5)再算一遍，观察阻力增量。很多机翼设计都尽量保持层流区域，原因就在这里：层流区的摩擦阻力能比湍流区低一个量级。

关于转捩位置，工程上常用 (Re_{x,tr} \approx 5 \times 10^5)（自然转捩）做粗略估计，更多时候结合粗糙度修正。表面污染物、铆钉凸起都会提前触发转捩，这一点在飞行器表面涂装设计里也极为关键。

4. 翼型升力产生的物理本质与库塔条件：从环量到升力线

4.1 库塔-茹科夫斯基定理与环量

绕翼型流动与绕圆柱流动的关键区别在于：翼型尖后缘通过库塔条件确定环量大小。库塔条件表述为：对于给定攻角的翼型，流动会在后缘光滑脱体，上下表面流速在下游趋于一致。这一条看似经验规则，实际上是由粘性边界层在后缘的“截止”作用决定的——无粘解有无数个，粘性解只有一个。

环量 (\Gamma) 与升力的关系由库塔-茹科夫斯基定理给出：

[ L' = \rho V_\infty \Gamma ]

其中 (L') 是单位展长升力。实际工程中更常用升力系数 (C_l)，它与环量的关系为：

[ C_l = \frac{2\Gamma}{V_\infty c} ]

其中 (c) 是翼型弦长。薄翼理论给出理想攻角范围内的线性关系：(C_l = 2\pi\alpha)（(\alpha) 用弧度表示）。真实粘性流动中，升力线斜率略低于此值，并且在大攻角时因流动分离而下降。

4.2 用涡面法算 NACA 0012 的升力系数

下面用一个简化但物理概念完整的涡面法（Vortex Panel Method）代码，对对称翼型 NACA 0012 做升力估算。这个方法不需要求解 N-S 方程，而是将翼型表面离散成若干涡面片，通过满足物面不可穿透条件和库塔条件求解涡强分布：

import numpy as np import math # NACA 0012 翼型坐标生成（对称翼型） def naca0012_coords(n_panels=60, chord=1.0): beta = np.linspace(0, 2*np.pi, n_panels*2, endpoint=False) x = 0.5 * chord * (1 - np.cos(beta)) # 余弦加密，前缘后缘点距小 y = 0.6 * chord * (0.2969*np.sqrt(x/chord) - 0.1260*(x/chord) - 0.3516*(x/chord)**2 + 0.2843*(x/chord)**3 - 0.1036*(x/chord)**4) # 只取上表面（从后缘到前缘）和下表面（从前缘到后缘） x_upper = x[:n_panels]; y_upper = y[:n_panels] x_lower = x[n_panels:][::-1]; y_lower = y[n_panels:][::-1] # 逆时针排列：下表面后缘->前缘，再上表面前缘->后缘 xs = np.concatenate([x_lower[::-1], x_upper[::-1]]) ys = np.concatenate([y_lower[::-1], y_upper[::-1]]) return xs, ys # 涡面法求解升力 def vortex_panel_lift(alpha_deg=5.0, n_panels=60): xs, ys = naca0012_coords(n_panels) n = len(xs) - 1 # 面板数 alpha = math.radians(alpha_deg) V_inf = 1.0 # 面板几何 Xc = (xs[:-1] + xs[1:]) / 2.0 Yc = (ys[:-1] + ys[1:]) / 2.0 dx = xs[1:] - xs[:-1] dy = ys[1:] - ys[:-1] S = np.hypot(dx, dy) tx = dx / S # 切向单位向量 ty = dy / S nx = ty # 法向（逆时针旋转90度） ny = -tx # 构建影响系数矩阵 A = np.zeros((n+1, n+1)) b = np.zeros(n+1) for i in range(n): for j in range(n): rx = Xc[i] - Xc[j] ry = Yc[i] - Yc[j] # 第j个涡面的诱导速度在i点处（解析解） # 由涡片诱导速度公式推导 # 这里用简化版：远场近似 + 奇异性处理 r = math.hypot(rx, ry) + 1e-8 # 涡片段引起的法向速度系数 A[i, j] = ry / (2*math.pi*r**2) * (-S[j]) # 简化示意 A[i, n] = 1.0 # 常数项对应来流法向分量 b[i] = V_inf * (math.sin(alpha)*nx[i] - math.cos(alpha)*ny[i]) # 库塔条件：后缘上下表面压力相等 # 简化为涡强之和为零 A[n, :n] = 1.0 A[n, n] = 0.0 b[n] = 0.0 # 解线性方程组 gamma = np.linalg.solve(A, b) # 总环量 Gamma = np.sum(gamma[:-1] * S) # 弦长 c = max(xs) - min(xs) Cl = 2 * Gamma / (c * V_inf) return Cl for a in [0, 2, 4, 6, 8, 10]: cl = vortex_panel_lift(a) print(f"攻角 {a:2d}°: C_l = {cl:.4f}")

说明一下：这是一个教学级简化实现，实际上的涡面法对影响系数的推导比这完整得多，包含面板自身的诱导速度解析式。但它展示了一条完整的路径：几何建模（余弦加密离散翼型）→ 奇点法布涡面 → 物面边界条件（法向速度为 0）→ 库塔条件封闭方程 → 解线性系统求环量 → 用库塔-茹科夫斯基定理算升力。

如果你把这段代码的A[i,j]部分换成完整影响系数公式，得到的曲线将非常接近薄翼理论值 (2\pi\alpha)。这正对应课程里强调的核心逻辑：升力不是“伯努利效应”的简单产物，而是粘性导致的后缘库塔条件决定环量，环量再决定升力。

4.3 有限翼展与升力线理论的失速边界

无限翼展的二维翼型没有翼尖涡，升力线斜率是每弧度 (2\pi)。真实机翼有展弦比限制，翼尖涡会诱导下洗速度，等效于减小了有效攻角。升力线理论给出：

[ C_L = \frac{2\pi\alpha}{1+2/AR} ]

其中 AR 是展弦比（展长的平方除以机翼面积）。AR 越小，升力线斜率越低。这对无人机设计很重要：多旋翼的旋翼叶片 AR 只有 3-5，升力效率远低于 AR 为 9-12 的固定翼。值得注意的是，升力线理论假设环量沿展向为椭圆分布，在失速前有效；一旦翼根或翼尖先失速，展向环量不再是椭圆分布，升力系数曲线就开始弯曲，最终到最大升力系数 (C_{L,max}) 后急剧下降。

5. 可压缩流动：激波、膨胀波与面积-速度关系

5.1 从不可压缩到可压缩：密度不再是常数

当马赫数超过 0.3 时，密度变化对流动的影响开始超过 5%，必须考虑压缩性。这一章在课程中占据重要位置，因为跨声速飞行器的设计瓶颈几乎都出在“局部马赫数超过 1”后出现的激波和激波-边界层干扰上。

等熵可压缩流的压力-速度关系由欧拉方程积分得到：

[ \frac{p_0}{p} = \left(1 + \frac{\gamma-1}{2}Ma^2\right)^{\gamma/(\gamma-1)} ]

其中 (\gamma=1.4) 是空气的比热比。当 (Ma \to 0)，这个公式做泰勒展开后正好回到不可压缩伯努利方程加上动压 (0.5\rho V^2)。你可以在下一节的脚本里对比不可压缩动压与可压缩动压的偏差趋势。

5.2 面积-速度关系的三种工况判别

一维定常等熵管流中，截面积变化与速度变化的关系为：

[ \frac{dA}{A} = (Ma^2 - 1)\frac{dV}{V} ]

这个式子给出了三个完全不同的物理区域：亚声速（Ma < 1）时面积减小速度增大；超声速（Ma > 1）时面积增大速度继续增大；声速处（Ma=1）面积取极值。超声速风洞的拉瓦尔喷管就是利用这个原理，先收缩到喉道达到声速，再扩张加速到超声速。如果你的仿真模型中出现“喉道处马赫数不是 1”，说明边界条件或网格质量有问题，排查方向通常在进出口压力的设定上。

5.3 正激波关系与总压损失的快速计算脚本

正激波是超声速来流减速为亚声速的最简单模型，正激波前后的马赫数关系为：

[ Ma_2^2 = \frac{Ma_1^2 + 2/(\gamma-1)}{2\gamma Ma_1^2/(\gamma-1) - 1} ]

总压比（衡量激波不可逆损失的重要指标）为：

[ \frac{p_{02}}{p_{01}} = \left[\frac{(\gamma+1)Ma_1^2}{2+(\gamma-1)Ma_1^2}\right]^{\gamma/(\gamma-1)} \left[\frac{\gamma+1}{2\gamma Ma_1^2 - (\gamma-1)}\right]^{1/(\gamma-1)} ]

写成 Python 来计算不同马赫数下的激波损失：

import numpy as np gamma = 1.4 def normal_shock(Ma1): """输入激波前马赫数，返回波后马赫数、静压比和总压比""" Ma2 = np.sqrt((Ma1**2 + 2/(gamma-1)) / (2*gamma*Ma1**2/(gamma-1) - 1)) p2_p1 = 1 + 2*gamma/(gamma+1) * (Ma1**2 - 1) # 总压比公式 p02_p01 = ((gamma+1)*Ma1**2 / (2 + (gamma-1)*Ma1**2))**(gamma/(gamma-1)) * \ ((gamma+1) / (2*gamma*Ma1**2 - (gamma-1)))**(1/(gamma-1)) return Ma2, p2_p1, p02_p01 for Ma in [1.2, 1.5, 2.0, 2.5, 3.0]: Ma2, p2p1, p02p01 = normal_shock(Ma) print(f"Ma1={Ma:.1f}: Ma2={Ma2:.3f}, p2/p1={p2p1:.2f}, 总压比={p02p01:.4f}")

从输出可以看到：来流马赫数从 1.2 升到 3.0，波后总压比从约 0.99 掉到约 0.33。这就是为什么超声速飞行器设计要尽量避免强激波——每道激波都在消耗总压，总压损失最终体现为阻力增加。超声速进气道的多级斜激波设计，就是把一道强正激波拆成多道弱斜激波加一道更弱的正激波，总压恢复率可以显著提升。

6. 用 Python 做课程公式的参数扫描：把 PDF 里的关系变成能用的工具

6.1 一口气对比理想气体、不可压缩与可压缩的压力系数

这一节给一个综合算例：以 NACA 0012 翼型在来流马赫数 0.3/0.6/0.8 三种条件下的表面压力分布估算为场景，对比不可压缩伯努利和可压缩等熵关系给出的压力系数差异。压力的无因次系数定义为：

[ C_p = \frac{p - p_\infty}{\frac{1}{2}\rho_\infty V_\infty^2} ]

不可压缩条件下，给定速度比 (V/V_\infty) 后，(C_p = 1 - (V/V_\infty)^2)。可压缩等熵条件下推导可得：

[ C_p = \frac{2}{\gamma Ma_\infty^2}\left{\left[1+\frac{\gamma-1}{2}Ma_\infty^2\left(1-\left(\frac{V}{V_\infty}\right)^2\right)\right]^{\gamma/(\gamma-1)} - 1\right} ]

import numpy as np import matplotlib.pyplot as plt gamma = 1.4 V_ratio = np.linspace(0.5, 1.5, 100) # 速度比 V/V_inf def cp_incompressible(Vr): return 1.0 - Vr**2 def cp_compressible(Vr, Ma): term = 1 + (gamma-1)/2 * Ma**2 * (1 - Vr**2) return 2/(gamma*Ma**2) * (term**(gamma/(gamma-1)) - 1) fig, ax = plt.subplots(figsize=(8,5)) ax.plot(V_ratio, cp_incompressible(V_ratio), 'k--', label='不可压缩') for Ma in [0.3, 0.6, 0.8]: cp = [cp_compressible(vr, Ma) for vr in V_ratio] ax.plot(V_ratio, cp, label=f'Ma={Ma}') ax.set_xlabel('V/V∞'); ax.set_ylabel('Cp') ax.legend(); ax.grid(True) plt.show()

运行后可以看到：Ma=0.3 时两条曲线几乎重合，Ma=0.6 时可压缩修正不可忽略，Ma=0.8 时吸力峰附近的 Cp 明显偏离不可压缩值。这提示一条工程经验：当翼型表面局部速度比达到 1.3 以上时，哪怕来流马赫数只有 0.5，局部也可能已经接近音速，必须用可压缩修正。

6.2 验证热力学一致性的自查方法

写这类脚本时，最容易出错的点是量纲和不一致的单位制。自查顺序按三层做：第一，确认所有长度用米、速度用米每秒、压力用帕斯卡；第二，检查极限行为——比如马赫数趋近于零时，可压缩 Cp 公式应该退化为不可压缩公式；第三，把已知实验值或 CFD 参考值拿来对比，比如 NACA 0012 在 Ma=0.3、攻角 0 度时的零升力条件，对称翼型无论如何压缩修正，Cp 分布都应该上下对称。如果不对称，说明翼型几何生成或边界条件方向有误，而不是公式问题。

6.3 把这个流程固化成自己的“空气动力学计算工具包”

建议把你写过的这些函数整理成独立的 Python 模块，接口统一为“输入马赫数、雷诺数、翼型几何、攻角，输出升阻力系数和边界层厚度”。这样后续遇到新的翼型或飞行工况，可以直接调参数跑结果，而不是每次重新翻 PDF 推导公式。对做仿真的人来说，这套脚本并不是替代 CFD，而是作为快速估算和 CFD 结果合理性检验的“第一道筛子”——先算一个粗值，再和求解器输出对比，能帮你提前发现边界条件设置错误、网格问题或湍流模型选择不当。这也正是刘沛清这门课反复强调的素养：先有物理直觉，再谈数值精度。

本文还有配套的精品资源，点击获取