流体力学基础

外流空气动力学：原理与诊断验证

外流建模的核心是判断哪些物理可以省略：压缩性、三维诱导阻力、边界层转捩与远场截断各有量化门槛。本文给出马赫数、普朗特-格劳厄特修正、展弦比与升力线斜率的核算方法及越界信号。

林工程师更新于 2026/9/22
CFD结果诊断与可信度验证物理建模与适用边界流体力学基础外流空气动力学阻塞修正升阻力系数压缩性修正诱导阻力进阶

外流空气动力学：原理与诊断验证

外流建模的难点不在求解器，而在决定哪些物理可以省略：来流是否可压缩、机翼是否必须按三维处理、边界层在何处转捩、远场截断放在多远。每一项都有可计算的判据，本文把它们换算成具体数字，并说明越界后会出现什么可观测的偏差，便于在建模阶段就确定模型层级。外流算例出现 20% 量级的气动力偏差时，湍流模型通常不是首因。计算域截断、风洞阻塞、首层网格落在缓冲层、以及把非定常尾迹当成定常解，这四类问题足以解释大部分与实验的差距。本文把每条判据换算成具体数字，并给出可以逐步执行的排查顺序。

亚声速压缩性对压力的放大

在临界马赫数以下，可用普朗特-格劳厄特修正把不可压压力系数换算到可压：

Cp=Cp,01−M∞2C_p = \frac{C_{p,0}}{\sqrt{1-M_\infty^{2}}}

修正因子随马赫数迅速增长：M=0.3M = 0.3 时 1/0.91=1.0481/\sqrt{0.91} = 1.048（+4.8%），M=0.5M = 0.5 时 1/0.75=1.1551/\sqrt{0.75} = 1.155（+15.5%），M=0.7M = 0.7 时 1/0.51=1.4001/\sqrt{0.51} = 1.400（+40%）。这说明即使整体不可压，只要局部加速到 M>0.3M > 0.3，翼型上表面的吸力峰就会被低估十几个百分点。对 NACA 0012 这类对称翼型，阻力发散马赫数在 $C_L \approx 0$ 附近约为 0.74，越过它激波阻力急剧上升，任何基于不可压假设的结论都失效。

三维效应与诱导阻力

有限展弦比机翼的升力线斜率低于二维值，并产生诱导阻力：

a=a01+a0π e AR,CD,i=CL2π e ARa = \frac{a_0}{1+\dfrac{a_0}{\pi\,e\,AR}}, \qquad C_{D,i} = \frac{C_L^{2}}{\pi\,e\,AR}

二维升力线斜率取 a0=2π=6.28 rad−1a_0 = 2\pi = 6.28\,\mathrm{rad^{-1}}，Oswald 效率 e=0.85e = 0.85：

  • AR=8AR = 8：a=6.28/(1+6.28/21.36)=4.85 rad−1=0.085 deg−1a = 6.28/(1+6.28/21.36) = 4.85\,\mathrm{rad^{-1}} = 0.085\,\mathrm{deg^{-1}}
  • AR=4AR = 4：a=6.28/(1+6.28/10.68)=3.95 rad−1=0.069 deg−1a = 6.28/(1+6.28/10.68) = 3.95\,\mathrm{rad^{-1}} = 0.069\,\mathrm{deg^{-1}}

同样的 $C_L = 0.5，诱导阻力分别是，诱导阻力分别是C_{D,i} = 0.25/21.36 = 0.0117（（AR = 8）与）与0.0234（（AR = 4）。若总阻力系数为）。若总阻力系数为C_D = 0.032，把，把AR = 4的机翼按二维算，会把诱导阻力当成型阻，的机翼按二维算，会把诱导阻力当成型阻，C_D低估0.023即70低估 0.023 即 70% 以上。判断方法：只要AR < 10且且C_L > 0.3$，就必须按三维建模或至少施加诱导阻力修正。

边界层转捩与湍流模型的匹配

低雷诺数翼型（RecRe_c 在 105∼10610^{5}\sim10^{6} 量级）上，层流分离泡与转捩位置直接决定阻力和失速特性。转捩对应的动量厚度雷诺数 Reθ,t=ρUθ/μRe_{\theta,t} = \rho U\theta/\mu 随来流湍流度 TuTu 显著下降，工程上常用的对应关系为：

风洞湍流度从 1% 升到 3%，转捩所需的 ReθRe_\theta 下降约 60%，转捩点明显前移，同一翼型在两个风洞里的阻力可以差 20%。平板算例中 U∞=30 m/sU_\infty = 30\,\mathrm{m/s}、x=0.25 mx = 0.25\,\mathrm{m} 给出 Rex=5×105Re_x = 5\times10^{5}，对应 Reθ,t≈510Re_{\theta,t} \approx 510，与 Tu≈1%Tu \approx 1\% 一行吻合，可作为入口湍流度取值的自检。

模型选择的判据：Rec>5×106Re_c > 5\times10^{6} 且来流湍流度正常时，全湍流 kk-ω\omega SST 足够；RecRe_c 在 105∼10610^{5}\sim10^{6} 且需要准确阻力时，应使用 γ\gamma-ReθRe_\theta 转捩模型或直接做 LES；把全湍流模型套在低雷诺数翼型上，计算阻力通常比实验高 30% 以上。

TuTu Reθ,tRe_{\theta,t}
0.1% 约 1070
1% 约 530
3% 约 210
6% 约 165

远场与近壁的分辨率要求

外流对首层网格的要求与近壁处理绑定。以 Rec=3×106Re_c = 3\times10^{6} 为例，估算壁面量并反算首层高度：

import math
rho, nu, U, c = 1.225, 1.5e-5, 50.0, 1.0
Re = U * c / nu
Cf = 0.0592 / Re**0.2          # 湍流平板
tau_w = 0.5 * Cf * rho * U**2  # Pa
u_tau = math.sqrt(tau_w / rho) # m/s
for yplus in (1, 30, 50):
    print(f"y+={yplus:2d}  y1={yplus*nu/u_tau*1e6:6.1f} um")
print(f"Re={Re:.2e}  Cf={Cf:.4f}  tau_w={tau_w:.2f} Pa  u_tau={u_tau:.2f} m/s")

输出给出 y+=1y^{+} = 1 对应 y1≈7.9 μmy_1 \approx 7.9\,\mu\mathrm{m}，y+=50y^{+} = 50 对应 395 μm395\,\mu\mathrm{m}。两种策略之间没有中间地带：用壁函数就必须保证首层落在 30∼30030\sim300，否则改用全解析近壁。远场方面，压力远场或自由来流边界应放在 20c20c 以上，入口 5c5c、出口 20c20c，否则升力斜率与阻力都会被域截断系统性地拉偏。

建模决策顺序

  1. 由 M∞M_\infty 与局部最大马赫数判断是否需要可压缩求解器。
  2. 由 ARAR 与 $C_L$ 判断是否需要三维建模或诱导阻力修正。
  3. 由 RecRe_c 与来流湍流度判断转捩是否重要，据此选湍流模型。
  4. 由 RecRe_c 反算 y+y^{+} 目标与首层高度，二选一确定近壁策略。
  5. 按 5c5c / 10c10c / 20c20c 布置入口、侧向与出口边界，并记录域尺寸。
  6. 把上述五项写入算例说明，后续任何一项改变都重做对应的量级核算。

什么时候不能再用不可压假设

可压缩性的强弱由马赫数决定，密度变化可由等熵关系估算：

M=U∞a∞,ρ0ρ=(1+0.2M2)2.5M = \frac{U_\infty}{a_\infty}, \qquad \frac{\rho_0}{\rho} = \left(1+0.2M^{2}\right)^{2.5}

海平面声速 a∞=340.3 m/sa_\infty = 340.3\,\mathrm{m/s}。核算三档：

  • M=0.3M = 0.3（U∞=102 m/sU_\infty = 102\,\mathrm{m/s}）：ρ0/ρ=(1.018)2.5=1.046\rho_0/\rho = (1.018)^{2.5} = 1.046，密度变化 4.6%，不可压求解器仍可用。
  • M=0.5M = 0.5（U∞=170 m/sU_\infty = 170\,\mathrm{m/s}）：密度变化 14%，已超出 5% 的经验门槛。
  • M=0.7M = 0.7（U∞=238 m/sU_\infty = 238\,\mathrm{m/s}）：密度变化 29%，必须用可压缩求解器并考虑激波。

首层网格：壁函数与缓冲层的错配

对 Rec=3.3×106Re_c = 3.3\times10^{6} 的翼型，用湍流平板估算壁面摩擦：

Cf≈0.0592Rex0.2,τw=12CfρU∞2,uτ=τw/ρ,y1=y+νuτC_f \approx \frac{0.0592}{Re_x^{0.2}}, \qquad \tau_w = \frac{1}{2}C_f\rho U_\infty^{2}, \qquad u_\tau = \sqrt{\tau_w/\rho}, \qquad y_1 = \frac{y^{+}\nu}{u_\tau}

代入 Rex=3.3×106Re_x = 3.3\times10^{6}，Rex0.2=20.2Re_x^{0.2} = 20.2，Cf=0.0029C_f = 0.0029，τw=0.5×0.0029×1.225×2500=4.4 Pa\tau_w = 0.5\times0.0029\times1.225\times2500 = 4.4\,\mathrm{Pa}，uτ=4.4/1.225=1.90 m/su_\tau = \sqrt{4.4/1.225} = 1.90\,\mathrm{m/s}。若目标 y+=1y^{+} = 1（低雷诺/全解析近壁），首层中心高度 y1=1×1.5×10−5/1.90=7.9 μmy_1 = 1\times1.5\times10^{-5}/1.90 = 7.9\,\mu\mathrm{m}；若目标 y+=50y^{+} = 50（标准壁函数），y1=395 μmy_1 = 395\,\mu\mathrm{m}。

真正的坑是"两头不靠"：用标准壁函数却把首层放在 y+=5∼10y^{+} = 5\sim10，壁面剪切被高估 20%～40%，阻力随之虚高。诊断办法是计算后按面积统计壁面 y+y^{+}，确认最小值、最大值与分布都落在所选处理方式的区间内。

故障模式与判定试验

现象 根因 判定试验
吸力峰 CpC_p 比实验低 15% 局部 M>0.3M>0.3 却用不可压求解 计算表面最大马赫数，超 0.3 换可压求解器
阻力在 M=0.75M=0.75 附近突增而计算平滑 求解器未捕捉激波，或网格过粗 加密激波附近网格并检查是否出现 MM 数跃变
AR=4AR=4 机翼的 $C_D$ 明显偏低 按二维处理，漏掉诱导阻力 用 CD,i=CL2/(πeAR)C_{D,i}=C_L^2/(\pi e AR) 估算漏掉多少
低雷诺数翼型阻力高 30% 以上 全湍流模型，无转捩 换 γ\gamma-ReθRe_\theta 模型，比较转捩位置
升力斜率对域尺寸敏感 远场太近或用了速度入口 域扩到 20c20c 并改压力远场复算
首层 y+y^{+} 落在 5～10 壁函数与网格不匹配 统计 y+y^{+} 分布，二选一重做近壁策略
$C_D$ 比实验高 25% 且随域缩小升高 计算域上下边界压缩流管，等效阻塞 把域扩到 10c10c 重算，若 $C_D$ 回落超过 2% 则确认
升力曲线斜率比薄翼理论低 15% 域太近、或远场边界用速度入口而非压力远场 换用压力远场/无反射边界复算，比较斜率
壁面 CfC_f 偏高但压差阻力正常 首层 y+y^{+} 落在 5～10 的缓冲层 统计 y+y^{+} 面积分布，改成全解析或放回对数区
$C_L$ 随时间周期性振荡 ±0.05 涡脱落未做时均统计 记录 $C_L$ 时间序列，取 10 个周期以上的均值
残差降到 10−510^{-5} 但气动力仍漂移 迭代/统计未收敛，或欠松弛过大 监控 $C_L、、C_D$ 的历史曲线是否出现平台
尾迹总压亏损与阻力积分不一致 后处理平面位置在分离区或回流区内部 把积分平面移到尾迹充分发展处，重算动量亏损

验收时要能回答的问题

  1. $C_L、、C_D$ 的参考面积与动压是否与实验口径一致？
  2. 域尺寸是否满足入口 5c5c、上下 10c10c、出口 20c20c 的最低要求？
  3. 壁面 y+y^{+} 分布是否与所选近壁处理匹配，且给出最大最小值？
  4. 非定常算例的统计窗口是否覆盖 10 个以上的脱落周期？
  5. 表面压力积分与尾迹动量积分给出的阻力是否在 5% 以内一致？
  6. 网格加密后气动力变化是否小于 2%，观测阶是否合理？

计算域尺寸与阻塞效应

风洞实验中模型占截面越大，来流被"堵"得越厉害，实测阻力偏高。二维模型在闭式风洞中的固体阻塞修正为：

ϵs=π248(ch)2,Uc=U∞(1+ϵs),CDcorr=CDmeas(1+ϵs)−2\epsilon_s = \frac{\pi^{2}}{48}\left(\frac{c}{h}\right)^{2}, \qquad U_c = U_\infty\left(1+\epsilon_s\right), \qquad C_D^{corr} = C_D^{meas}\left(1+\epsilon_s\right)^{-2}

hh 为风洞高度或 CFD 计算域高度。核算两档弦高比：

  • c/h=0.2c/h = 0.2：ϵs=0.2056×0.04=0.82%\epsilon_s = 0.2056\times0.04 = 0.82\%，阻力修正约 1.6%，可忽略。
  • c/h=0.5c/h = 0.5：ϵs=0.2056×0.25=5.1%\epsilon_s = 0.2056\times0.25 = 5.1\%，阻力修正约 10%，不修正就会把阻塞误读成"模型偏保守"。

CFD 中并不存在真实的壁面阻塞，但计算域上下边界给成滑移/对称面时同样会压缩流管，效果与阻塞等价。稳妥做法是把上下边界放到 10c10c 以外，入口 5c5c、出口 20c20c 以外，并在验收表中记录域尺寸；若把域缩小一半后 $C_D$ 上升超过 2%，说明域尺寸已经在污染结果。

非定常尾迹的统计窗口

钝体或大攻角翼型会周期性脱落涡，瞬时升力永远在振荡。以圆柱绕流为例：

St=fDU∞St = \frac{fD}{U_\infty}

亚临界圆柱 St≈0.2St \approx 0.2。取 D=0.1 mD = 0.1\,\mathrm{m}、U∞=50 m/sU_\infty = 50\,\mathrm{m/s}，脱落频率 f=0.2×50/0.1=100 Hzf = 0.2\times50/0.1 = 100\,\mathrm{Hz}，周期 10 ms。要得到稳定时均量，至少需要 10 个周期即 0.1 s 的物理时间，时间步应小于 1/(20f)=0.5 ms1/(20f) = 0.5\,\mathrm{ms}。若报告里给的是某一时刻的瞬时 $C_L，或采样窗口不足5个周期，，或采样窗口不足 5 个周期，C_L$ 的标准差会达到均值的 10% 以上。

用尾迹动量收支做独立核对

远场尾迹中的动量亏损应等于模型受到的阻力，这是与表面压力积分完全独立的第二条证据：

D=ρ∞U∞∫wake(U∞−u) dyD = \rho_\infty U_\infty \int_{wake} \left(U_\infty - u\right)\,\mathrm{d}y

在翼型后方 1c1c 处取积分平面，若表面压力积分给出 $C_D = 0.0060，而尾迹积分给出，而尾迹积分给出0.0052$（低 13%），常见原因是积分平面落在分离泡内部或网格在尾迹方向过粗。两个数字一致到 5% 以内，才说明阻力不是后处理口径的产物。

OpenFOAM 中可用的力系数函数对象：

functions
{
    forceCoeffs
    {
        type            forceCoeffs;
        libs            (forces);
        patches         (airfoil);
        rho             rhoInf;
        rhoInf          1.225;
        liftDir         (0 1 0);
        dragDir         (1 0 0);
        CofR            (0.25 0 0);   // 四分之一弦
        pitchAxis       (0 0 1);
        magUInf         50;
        lRef            1.0;          // 弦长
        Aref            1.0;          // 二维: c * 1
        writeControl    timeStep;
        timeInterval    1;
    }
}

参考资料

升阻力系数的定义本身就是一个误差来源： $$ C_L = \frac{L}{\frac{1}{2}\rho_\infty U_\infty^{2} A}, \qquad C_D = \frac{D}{\frac{1}{2}\rho_\infty U_\infty^{2} A} 二维翼型取 A=c×1A = c\times1（弦长乘单位展长），三维机翼取参考面积（通常为机翼平面投影面积），二者混用会让 $C_D差出一个展弦比。取NACA0012、差出一个展弦比。取 NACA 0012、c = 1,\mathrm{m}、、U_\infty = 50,\mathrm{m/s}、海平面空气、海平面空气\rho = 1.225,\mathrm{kg/m^3}、、\nu = 1.5\times10^{-5},\mathrm{m^2/s}$：

  • Rec=50×1/1.5×10−5=3.3×106Re_c = 50\times1/1.5\times10^{-5} = 3.3\times10^{6}
  • 动压 q∞=12×1.225×502=1531 Paq_\infty = \frac{1}{2}\times1.225\times50^{2} = 1531\,\mathrm{Pa} 若测得的阻力为 D=9.2 ND = 9.2\,\mathrm{N}（单位展长），则 $C_D = 9.2/1531 = 0.0060。对比。对比Re = 3.3\times10^{6}下NACA0012的公开阻力极曲线（下 NACA 0012 的公开阻力极曲线（C_D \approx 0.0065，含尾缘修正），偏差约8，含尾缘修正），偏差约 8%，属于可接受范围；若算出C_D = 0.020$，必须先怀疑参考面积与动压口径，而不是模型。
  1. Glauert H., The Elements of Aerofoil and Airscrew Theory, Cambridge University Press, 1926.
  2. Abbott I.H., von Doenhoff A.E., Theory of Wing Sections, Dover, 1959.
  3. Michel R., “Etude de la transition sur les profils d’aile,” ONERA Report 1/1578A, 1951.
  4. Abu-Ghannam B.J., Shaw R., “Natural Transition of Boundary Layers — The Effects of Turbulence, Pressure Gradient, and Flow History,” Journal of Mechanical Engineering Science, 1980.
  5. Maskell E.C., “A Theory of the Blockage Effects on Bluff Bodies and Stalled Wings in a Closed Wind Tunnel,” Aeronautical Research Council R&M 3400, 1963.
  6. Barlow J.B., Rae W.H., Pope A., Low-Speed Wind Tunnel Testing, 3rd ed., Wiley, 1999.
  7. Roshko A., “On the Development of Turbulent Wakes from Vortex Streets,” NACA Report 1191, 1954.
  8. Menter F.R., “Two-Equation Eddy-Viscosity Turbulence Models for Engineering Applications,” AIAA Journal, 1994.