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

拟合窗一开过 10/L 和 0.1/η,本来该是 −5/3 的斜率会趴到 −1.4。

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

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

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

把 time speed 滑块推高,破碎的间隔会缩短,从大涡到小涡的顺序不变。

Each eddy lives a turnover time τ ∝ ℓ^(2/3), then splits into smaller ones. The smallest eddies (gray) dissipate as heat.

time speed: 1.0×

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

三把尺子塞进同一道流场

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

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

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

量纲分析架起桥梁

Kolmogorov的假设只有两行。

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

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

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

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

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

E(k) = C\, \varepsilon^{2/3}\, k^{-5/3}

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

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

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

\frac{\eta}{L} = \mathrm{Re}^{-3/4}, \quad \frac{\lambda}{L} \sim \mathrm{Re}^{-1/2}, \quad \frac{\tau_\eta}{T} \sim \mathrm{Re}^{-1/2}

其中 $\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}$ 单元 — 当今超算亦无法企及。

拟合窗的十进制宽度由边界直接定下。 $k$ 从 $10/L$ 到 $0.1/\eta$ ,宽度是 $\log_{10}(L/\eta)-2$ 个 decade。 $\mathrm{Re}=10^{4}$ 时 $L/\eta\approx 1{,}000$ ,窗宽是 1 个 decade。Reynolds 每增大十倍, $\eta/L$ 乘上 $0.18$ 。这个收缩用两次,到 $\mathrm{Re}=10^{6}$ , $\eta/L$ 大约变成 $0.032$ 倍, $L/\eta$ 拉开到约 $3\times 10^{4}$ ,窗宽约 2.5 个 decade。滑块上 $\mathrm{Re}=10^{6}$ 的直线看起来接近 4 个 decade,是因为窗外面的峰和粘性尾巴也进了刻度。

代码里的常数取 $u'=1$ 、 $L=1$ ,则 $\varepsilon=u'^{3}/L=1$ 。乘上 $f_L$ 与 $f_\eta$ 之前,惯性段是 $E(k)=1.5\,k^{-5/3}$ 。注释里的拟合斜率 −1.661 和期望值 −1.667 相差 0.006。噪声加到 10%,只要守住窗口,三位还在。

这段间隔在圆管里第一次打开的过程,接在水流散开、边界层分成多层的临界雷诺数后面。名义 Reynolds 越大, $L$ 与 $\eta$ 之间要放的网格按 $\mathrm{Re}^{9/4}$ 增加, $10^{4}$ 的三维已经是 $10^{9}$ 个单元。

能量流动的道路

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

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

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

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

把 Reynolds 数滑块往上推,-5/3 直线段会变长。

Model energy spectrum E(k) on log-log axes (Pope 2000 form).

Reynolds number: Re = 10^5.0 ≈ 1.0e+5 show -5/3 reference line

L/η = 5.62e+3 · η = 1.78e-4 · λ (Taylor) = 1.22e-2

$\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全部解出。在哪里下刀,就是建模的判断。

把拟合窗开到两端时斜率会趴到 −1.4#

$\eta/L = \mathrm{Re}^{-3/4}$ 决定窗有多宽。掩码不放在 $10/L < k < 0.1/\eta$ 里面, $f_L$ 和 $f_\eta$ 就把两端扳弯,本该只跟通量 $\varepsilon$ 走的斜率会落到 −1.4。代码注释里的 −1.661,是守住那扇窗时的数字。

即便 $\mathrm{Re} = 10^4$ ,三维网格已经很沉, $\mathrm{Re} = 10^7$ 时格子走到 $10^{15}$ 。在那之下切开、交给亚格子时,系数漏成负数的地方写在Germano 恒等式里动态系数变成负数的位置。相信级联通量,和把切完的系数原样拿去用,是两件不同的决定。