简介《空气动力学基础北航精品课程》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 雷诺数、马赫数和克努森数的实际判断准则无量纲数的意义在于只要两个流动的无量纲控制参数一致即使尺度不同流场也是相似的。空气动力学里必须心里有数的三个数是参数表达式工程判断雷诺数 (Re)( \rho V L / \mu )判断层流/湍流影响阻力构成马赫数 (Ma)( V / a )判断压缩性是否显著克努森数 (Kn)( \lambda / L )Kn 0.01 时连续介质假设失效以标准海平面条件为例空气密度约为 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(fRe_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_panels60, chord1.0): beta np.linspace(0, 2*np.pi, n_panels*2, endpointFalse) 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_deg5.0, n_panels60): 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((n1, n1)) b np.zeros(n1) 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}{12/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)} ]其中 (\gamma1.4) 是空气的比热比。当 (Ma \to 0)这个公式做泰勒展开后正好回到不可压缩伯努利方程加上动压 (0.5\rho V^2)。你可以在下一节的脚本里对比不可压缩动压与可压缩动压的偏差趋势。5.2 面积-速度关系的三种工况判别一维定常等熵管流中截面积变化与速度变化的关系为[ \frac{dA}{A} (Ma^2 - 1)\frac{dV}{V} ]这个式子给出了三个完全不同的物理区域亚声速Ma 1时面积减小速度增大超声速Ma 1时面积增大速度继续增大声速处Ma1面积取极值。超声速风洞的拉瓦尔喷管就是利用这个原理先收缩到喉道达到声速再扩张加速到超声速。如果你的仿真模型中出现“喉道处马赫数不是 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{(\gamma1)Ma_1^2}{2(\gamma-1)Ma_1^2}\right]^{\gamma/(\gamma-1)} \left[\frac{\gamma1}{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/(gamma1) * (Ma1**2 - 1) # 总压比公式 p02_p01 ((gamma1)*Ma1**2 / (2 (gamma-1)*Ma1**2))**(gamma/(gamma-1)) * \ ((gamma1) / (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(fMa1{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, labelfMa{Ma}) ax.set_xlabel(V/V∞); ax.set_ylabel(Cp) ax.legend(); ax.grid(True) plt.show()运行后可以看到Ma0.3 时两条曲线几乎重合Ma0.6 时可压缩修正不可忽略Ma0.8 时吸力峰附近的 Cp 明显偏离不可压缩值。这提示一条工程经验当翼型表面局部速度比达到 1.3 以上时哪怕来流马赫数只有 0.5局部也可能已经接近音速必须用可压缩修正。6.2 验证热力学一致性的自查方法写这类脚本时最容易出错的点是量纲和不一致的单位制。自查顺序按三层做第一确认所有长度用米、速度用米每秒、压力用帕斯卡第二检查极限行为——比如马赫数趋近于零时可压缩 Cp 公式应该退化为不可压缩公式第三把已知实验值或 CFD 参考值拿来对比比如 NACA 0012 在 Ma0.3、攻角 0 度时的零升力条件对称翼型无论如何压缩修正Cp 分布都应该上下对称。如果不对称说明翼型几何生成或边界条件方向有误而不是公式问题。6.3 把这个流程固化成自己的“空气动力学计算工具包”建议把你写过的这些函数整理成独立的 Python 模块接口统一为“输入马赫数、雷诺数、翼型几何、攻角输出升阻力系数和边界层厚度”。这样后续遇到新的翼型或飞行工况可以直接调参数跑结果而不是每次重新翻 PDF 推导公式。对做仿真的人来说这套脚本并不是替代 CFD而是作为快速估算和 CFD 结果合理性检验的“第一道筛子”——先算一个粗值再和求解器输出对比能帮你提前发现边界条件设置错误、网格问题或湍流模型选择不当。这也正是刘沛清这门课反复强调的素养先有物理直觉再谈数值精度。本文还有配套的精品资源点击获取