Taylor–Green 涡:结果诊断与可信度验证
Taylor–Green 涡是少数在黏性条件下有严格解析解、又能同时检验对流与耗散的无网格基准。它的可信度判据不是"涡看起来圆不圆",而是三条可量化的线:动能衰减率是否等于 4νk2、速度场 L2 误差是否随分辨率按预期阶下降、以及反推出的等效数值黏性是否随时间步收敛。本文给出这三条线的公式、阈值与一次完整手算。
解析衰减律是判断耗散是否过量的标尺
在 2π×2π 的周期方腔内,二维 Taylor–Green 涡的速度场为
ux=Ue−2νk2tsin(kx)cos(ky),uy=−Ue−2νk2tcos(kx)sin(ky),
其中 U 为初始速度幅值(m/s),ν 为运动黏度(m2/s),波数取 k=1 m−1(对应域长 2π m)。因为 u∼e−2νk2t,体积平均动能按 u2 衰减:
Ek(t)=Ek(0)e−4νk2t,λexact=4νk2.
λexact 的单位是 s−1,它是判断数值耗散是否过量的唯一标尺。取 U=1 m/s、ν=0.01 m2/s、k=1 m−1,则 λexact=4×0.01×1=0.04 s−1,对应雷诺数 Re=UL/ν=1×1/0.01=100。在 t=1 s 时动能比为 e−0.04=0.9608,t=5 s 时降为 e−0.2=0.8187。若数值曲线在 t=5 s 只剩初始动能的 0.70,说明多出的耗散已经把物理衰减淹没了。
反推等效数值黏性
从数值结果拟合对数动能斜率,得到总衰减率 λnum,定义等效数值黏性为
νnum=4k2λnum−λexact.
该式把"看起来衰减快一点"翻译成有量纲的黏性增量,便于与 ν 比较。若 λnum=0.046 s−1,则 νnum=(0.046−0.04)/(4×1)=1.5×10−3 m2/s,即数值黏性占物理黏性的 15%。工程判据是 νnum/ν≤0.1;超过这一比例时,改变 ν 得到的结论已不可信。
速度场 L2 误差与观测收敛阶
在固定时刻(建议取 t=2 s,此时涡尚未衰减到数值噪声水平)对速度场做加权 L2 误差:
EL2=(∑iVi∣uiexact∣2∑iVi∣ui−uiexact∣2)1/2,
Vi 为粒子体积。三档分辨率下的 $E_1,E_2,E_3$ 给出观测收敛阶
p=ln(Δx2/Δx1)ln(E2/E1).
一次可核对的手算:域长 2π m,三档分辨率分别取 642,1282,2562 粒子,对应 Δx=0.0982,0.0491,0.0245 m。在 t=2 s 测得 $E_1=0.0042、E_2=0.0119、E_3=0.0312$。则
p12=ln2ln(0.0119/0.0042)=0.6931.041=1.50,p23=ln2ln(0.0312/0.0119)=0.6930.964=1.39.
两段阶次落在 1.4~1.5,说明二阶格式在该分辨率区间被边界与初始化的低阶误差稀释;若 p12<1,优先检查周期边界是否真正闭合,以及时间步是否小到让时间离散误差退出主导。
时间步与分辨率必须分开扫
Taylor–Green 涡的诊断最容易出错的地方,是把空间加密和时间步缩小同时做。若 Δt 由 CFL 数固定,则 Δt∝Δx,空间与时间误差一起变小,反推出的 p 会被系统性抬高。正确做法是先固定 Δx,把 Δt 从 0.01 s 逐步缩到 0.001 s,确认 EL2 变化小于 1%,再固定 Δt 扫 Δx。只有在这两步都稳定后,上式的 p 才有解释力。
诊断脚本
import numpy as np
nu, k, U, T = 0.01, 1.0, 1.0, 2.0 # m^2/s, 1/m, m/s, s
lam_exact = 4 * nu * k**2 # 0.04 1/s
def numerical_viscosity(lam_num, k=1.0, lam_exact=lam_exact):
"""由实测衰减率反推等效数值黏性"""
return (lam_num - lam_exact) / (4 * k**2)
def fit_decay(t, Ek):
"""对 ln(Ek) 做线性拟合得到斜率 -lambda"""
slope, _ = np.polyfit(t, np.log(Ek), 1)
return -slope
def order(E1, E2, dx1, dx2):
return np.log(E2 / E1) / np.log(dx2 / dx1)
print("Ek(t=1)/Ek0 =", np.exp(-lam_exact * 1.0)) # 0.9608
print("nu_num =", numerical_viscosity(0.046)) # 1.5e-3 m^2/s
print("p12 =", order(0.0042, 0.0119, 0.0982, 0.0491)) # 1.50
失败模式与判定试验
| 现象 |
根因 |
判定试验 |
| 动能衰减明显快于 e−0.04t |
人工黏性或核函数耗散过强 |
拟合 λnum 并算 νnum/ν,看是否超 0.1 |
| 涡在 t=3 s 后停止衰减 |
空间分辨率不足,涡核被抹成常数场 |
比较 1282 与 2562 的 t=5 s 动能 |
| p12 高于理论阶 |
Δt 随 Δx 同步缩小,时间误差混入 |
固定 Δx 单独缩 Δt,看 EL2 是否已平台 |
| 速度剖面出现棋盘振荡 |
粒子无序初始化未做正则化 |
用规则格点与扰动格点各跑一次比较 EL2 |
| 动能曲线有周期性抖动 |
采样频率与涡对流频率接近,混叠 |
把输出间隔从 0.02 s 加密到 0.005 s |
| 反推黏性为负 |
拟合窗口太短或包含初始化瞬态 |
把拟合起点从 t=0 后移到 t=0.5 s |
复核与参考文献
判定通过需同时满足 νnum/ν≤0.1、p∈[1.3,2.2] 且单调、固定 Δx 下 EL2 对 Δt 的平台化。参考文献:
- Taylor, G. I., Green, A. E., “Mechanism of the production of small eddies from large ones,” Proceedings of the Royal Society A, 158(895), 1937, pp. 499–521.
- Ethier, C. R., Steinman, D. A., “Exact fully 3D Navier–Stokes solutions for benchmarking,” International Journal for Numerical Methods in Fluids, 19(5), 1994, pp. 369–375.
- Chorin, A. J., “Numerical solution of the Navier–Stokes equations,” Mathematics of Computation, 22(104), 1968, pp. 745–762.
- Brachet, M. E., et al., “Small-scale structure of the Taylor–Green vortex,” Journal of Fluid Mechanics, 130, 1983, pp. 411–452.
- Monaghan, J. J., “Smoothed Particle Hydrodynamics,” Annual Review of Astronomy and Astrophysics, 30, 1992, pp. 543–574.
- Violeau, D., Fluid Mechanics and the SPH Method: Theory and Applications, Oxford University Press, 2012.
- Morris, J. P., Fox, P. J., Zhu, Y., “Modeling low Reynolds number incompressible flows using SPH,” Journal of Computational Physics, 136(1), 1997, pp. 214–226.