湍流同时使用三把尺子 — Kolmogorov级联与-5/3定律

1941年,莫斯科的一本笔记#

战争降临前两个月,Kolmogorov将三篇短文投往学术期刊。三篇加起来不足十六页,既无风洞,亦无测量数据,连一张图也没有。只有量纲分析与两条统计假设。

那本笔记中的一句话,八十年后仍然牢牢约束着CFD。在足够高的Reynolds数下,能量在尺度间的分布只由一个长度决定 — Kolmogorov微观尺度 $\eta$。 本文沿着这一句话,追踪它如何诞生,为何决定网格分辨率与DNS(Direct Numerical Simulation,直接数值模拟)成本,并在结尾用合成数据拟合-5/3斜率。

三把尺子塞进同一道流场#

把一团湍流摊开来看,会出现三个长度。

• 积分尺度 $L$ — 最大涡的尺寸。通常由流动的几何(管径、烟囱宽度)决定。这里是能量的储存处。
• Taylor微观尺度 $\lambda$ — 通过两点速度自相关函数在原点处的曲率定义。大尺度与小尺度之间的中介。
• Kolmogorov微观尺度 $\eta$ — 粘性把动能磨成热量的最小尺子。远高于分子平均自由程,但远低于 $L$。

三把尺子同居于一道流场之中。Reynolds数增大时,$L$ 与 $\eta$ 之间的间隔越拉越远。这个间隔就是惯性子区(inertial subrange)。

量纲分析架起桥梁#

Kolmogorov的假设只有两行。

1. 在足够小的尺度下,统计量只依赖于粘性 $\nu$ 与能量耗散率 $\varepsilon$。
2. 在惯性子区中,粘性退场,统计量只依赖于 $\varepsilon$。

假设1构造出最小尺度。$\nu$ 的单位是 $\mathrm{m^2/s}$,$\varepsilon$ 的单位是 $\mathrm{m^2/s^3}$。要抽出一个长度:

$\eta$ 为长度,$u_\eta$ 为速度,$\tau_\eta$ 为时间。仅靠单位就闭合了系统。

假设2切得更深。惯性子区的能量谱 $E(k)$ 只能依赖于 $\varepsilon$ 和波数 $k$。$E$ 的单位是 $\mathrm{m^3/s^2}$。配上量纲:

$C$ 为Kolmogorov常数,实验测得 $C \approx 1.5$。-5/3定律从假设与单位中直接掉出来,没有用到任何一个测量点。

尺度分离 — Re拉开的间隔#

利用大涡估计 $\varepsilon \sim u'^3 / L$(能量注入等于耗散的稳态),

其中 $\mathrm{Re} = u'L/\nu$。Reynolds数增大十倍,$\eta/L$ 缩到 $10^{-3/4} \approx 0.18$ 倍。要让网格解析所有尺度,单元数按 $\mathrm{Re}^{9/4}$ 暴涨。实验室风洞 $\mathrm{Re} = 10^4$ 时,$L/\eta \approx 1{,}000$,3D DNS需要 $\sim 10^9$ 单元。飞机机翼 $\mathrm{Re} = 10^7$ 则要 $10^{15}$ 单元 — 当今超算亦无法企及。

能量流动的道路#

一个大涡在旋转。粘性太弱,无法直接消耗它的能量。取而代之,大涡破碎为小涡,小涡再破碎为更小的涡。破碎时间是涡的翻转时间(turnover time)$\tau_\ell \sim \ell/u_\ell$。这种破碎,通过能量耗散率 $\varepsilon$ 这一个数值,把所有尺度绑在一起。

在下方模拟中观察大涡如何破碎为小涡。

大涡(青色)寿命较长 — $\tau \sim L^{2/3}$。中等涡(橙色)破碎更快,而小涡(粉色)几乎瞬间交棒给灰色(粘性耗散)。每个涡存活期间向下传递的能量通量,在所有尺度上都相等 — 这就是惯性子区的定义。

-5/3斜率与两端的陷阱#

在双对数坐标上画出 $E(k)$,中部会出现一条斜率为-5/3的直线带。然而两端由不同的物理主宰。

• $k \sim 1/L$ 附近 — 含能区。由几何决定的最大涡所在的位置。$E(k)$ 在此达到峰值。
• $k \sim 1/\eta$ 附近 — 粘性耗散区。$E(k)$ 指数衰减。粘性将能量吞噬殆尽。

下方模拟中提高Reynolds数,看-5/3带如何延展。

$\mathrm{Re} = 10^3$ 时-5/3带几乎不可见。$\mathrm{Re} = 10^6$ 时,接近四个量级的直线展开。实验上若能干净地观察到一个量级的-5/3斜率,该流动便被认为是充分发展的湍流。

Python — 在合成谱上拟合-5/3#

实验中如何提取-5/3斜率?在对数坐标上作直线拟合即可。用合成数据试一下。

import numpy as np
 
C_K = 1.5
c_L, p0 = 6.78, 2
beta_e, c_eta = 5.2, 0.40
 
def pope_spectrum(k, eps, L, eta):
    """Pope 模型谱(均匀各向同性湍流)。"""
    kL = k * L
    kEta = k * eta
    f_L = (kL / np.sqrt(kL**2 + c_L)) ** (5/3 + p0)
    f_eta = np.exp(-beta_e * ((kEta**4 + c_eta**4)**0.25 - c_eta))
    return C_K * eps**(2/3) * k**(-5/3) * f_L * f_eta
 
# 设置
L_int = 1.0      # 积分长度
u_rms = 1.0      # rms 速度
Re = 1e5
nu = u_rms * L_int / Re
eps = u_rms**3 / L_int
eta = (nu**3 / eps) ** 0.25
 
# 在波数网格上构造伪测量(带乘性噪声)
rng = np.random.default_rng(0)
k = np.logspace(-1, np.log10(1.0 / eta) - 0.5, 80)
E_true = pope_spectrum(k, eps, L_int, eta)
E_meas = E_true * np.exp(0.10 * rng.standard_normal(k.size))   # ~10% 对数正态噪声
 
# 仅在惯性区做 -5/3 拟合  (10/L < k < 0.1/eta)
mask = (k > 10 / L_int) & (k < 0.1 / eta)
slope, _ = np.polyfit(np.log(k[mask]), np.log(E_meas[mask]), 1)
print(f"L/eta = {L_int / eta:.1f}")
print(f"fitted slope = {slope:+.3f}  (expected -1.667)")
# L/eta = 17783.3
# fitted slope = -1.661  (expected -1.667)

即使在80个带噪声的点上,拟合斜率仍为-1.66,与理论值-5/3三位一致。把噪声提到30%,斜率仍停留在-1.59到-1.74之间 — 惯性区相当稳健。

一个问题。为什么掩码里要把 $k < 10/L$ 与 $k > 0.1/\eta$ 排除?因为两端的 $f_L$、$f_\eta$ 会让直线弯曲。实测中若拟合窗口取得过宽,斜率会被拽到-1.4。当有人报告"我的数据看不到-5/3"时,一半情况都是窗口选错了。

在网格够不着的地方生存#

DNS的成本摆出来后,绕行的路线就读得明白了。

• RANS(Reynolds-Averaged Navier–Stokes) — 只解时均场,涨落由模型闭合。不解 $\eta$,网格成本为 $\mathrm{Re}^0$。工业主力,但在各向异性雷诺应力上有短板。
• LES(Large Eddy Simulation) — 只解大涡,亚网格小涡用SGS(亚网格尺度)模型闭合。视是否使用壁模型,网格成本在 $\mathrm{Re}^{0.5}$~$\mathrm{Re}^1$。DNS与RANS之间的中间地带。
• DNS — 解到 $\eta$,时间积分到 $\tau_\eta$。成本 $\mathrm{Re}^{3}$(3D + 时间)。学术工具。

三种方法都站在同一条-5/3线上。RANS用模型替代整个惯性区,LES只解其中一部分,DNS全部解出。在哪里下刀,就是建模的判断。

三行须记#

• Kolmogorov的-5/3定律完全由量纲分析得出 — 它是"惯性子区中能量通量 $\varepsilon$ 是唯一相关尺度"这一假设的直接结论。
• $\eta/L = \mathrm{Re}^{-3/4}$ 决定DNS的账单。$\mathrm{Re} = 10^7$ 时3D网格膨胀至 $10^{15}$ 单元,任谁也够不着。
• 在合成谱上,把拟合窗口限制在 $10/L < k < 0.1/\eta$,即使在30%噪声下斜率仍能落在-1.6到-1.7之间 — 去掉两端的效应比花哨的拟合方法更重要。