摘要

Sabra 壳模型是研究湍流能量级联、标度律与间歇性现象的经典简化模型。相较于直接数值模拟不可压缩 Navier–Stokes 方程，壳模型能够在较低计算成本下保留多尺度非线性交互结构，因此常被用于高雷诺数湍流统计性质的数值研究。本文构建了一套基于自适应指数时间差分四阶 Runge–Kutta 方法（ETD-RK4）的 Sabra 壳模型数值实验框架，系统集成了 CFL 步长约束、Richardson 误差估计、真拒步回滚、恒功率注入驱动以及 realization 系综平均统计。

本文采用 40 个壳层、波数比 \lambda=2、Sabra 参数 \varepsilon=0.5，在低粘性系数 \nu=3\times 10^{-11} 和恒功率注入 P_{\mathrm{inj}}=0.005 的条件下进行 20 组独立 realization 数值实验。对稳态后样本进行系综平均，考察能谱 E_n=\langle |u_n|^2\rangle、p 阶结构函数 S_p(n)=\langle |u_n|^p\rangle、标度指数 \zeta_p 以及扩展自相似性（ESS）指数 \beta_p。在主拟合区间 n\in[8,22] 上，实验结果给出能谱指数约为 \alpha\approx 1.62\pm 0.03，与 Kolmogorov 1941 理论预测值 5/3 接近，但存在轻微偏离。结构函数结果显示出明显的间歇性修正，高阶统计量偏离线性关系 \zeta_p=p/3，并在低至中等阶数范围内与 She–Leveque 型经验模型表现出相近趋势。ESS 分析进一步降低了拟合不确定性，提高了标度指数提取的稳定性。

本文结果表明，在当前数值参数与拟合区间下，Sabra 壳模型能够复现典型湍流标度行为及其间歇性修正。需要强调的是，本文研究对象是 Sabra 壳模型的数值统计行为，而非真实三维不可压缩湍流本身；因此本文结论应理解为壳模型层面的数值实验结果，而不能直接视为对 Navier–Stokes 湍流理论的严格证明。

关键词：Sabra 壳模型；湍流；ETD-RK4；标度律；间歇性；扩展自相似性

1 引言

湍流统计理论是经典流体力学中的核心问题之一。Kolmogorov 于 1941 年提出局地各向同性高雷诺数湍流理论，给出著名的 -5/3 能谱律以及结构函数线性标度律：

E(k)\sim k^{-5/3}, \qquad \zeta_p=\frac{p}{3}.

这一理论成功描述了惯性区内的平均标度行为，但实验与数值研究均表明，实际湍流在高阶统计量上存在系统性偏离，即所谓间歇性现象。间歇性反映了小尺度耗散与能量输运过程的不均匀性，是现代湍流理论中的关键难点之一。

直接数值模拟（DNS）能够在无模型假设下解析 Navier–Stokes 方程的多尺度结构，但其计算代价极高，特别是在高雷诺数情形下。与此相比，壳模型通过将波数空间离散为对数排列的壳层，并保留局部壳间非线性耦合，在显著降低计算成本的同时，仍能够再现能量级联、惯性区标度与间歇性等重要统计特征。因此，壳模型长期以来被视为连接理论、数值与统计分析的重要中间层工具。

Sabra 壳模型是 GOY 类壳模型的重要改进版本之一。与早期壳模型相比，Sabra 模型在相位结构与统计行为方面通常更稳定，也更适合用于研究结构函数与异常标度指数。基于 Sabra 模型的大量数值研究表明，其在适当参数范围内可较好地再现类似于真实湍流的统计规律，尤其适合研究高雷诺数近似条件下的能量级联和间歇性。

本文的目标不是对湍流理论作严格证明，而是构建一套更严谨的 Sabra 壳模型数值实验框架，并在低粘性、恒功率注入、系综平均与自适应时间推进的条件下，对模型的统计特性进行系统分析。与只给出单次轨道结果的数值实验不同，本文强调以下几点：

采用指数时间差分四阶方法处理粘性刚性项；

引入 CFL 限制与 Richardson 误差控制，配合拒步回滚提高数值稳定性；

使用恒功率注入强迫替代简单恒幅驱动，以增强统计可比性；

通过 20 次独立 realization 的系综平均降低随机初值波动的影响；

同时分析能谱、结构函数、ESS 以及 realization 之间的统计离散度。

本文结构如下。第 2 节介绍 Sabra 壳模型及其数值离散方案；第 3 节说明实验设置、采样方法与诊断量定义；第 4 节给出数值结果，包括能谱拟合、结构函数与 ESS 分析；第 5 节讨论结果的物理意义及模型局限；第 6 节给出结论。

2 模型与数值方法

2.1 Sabra 壳模型

Sabra 壳模型将波数空间离散为对数排列的壳层：

k_n = k_0 \lambda^n,\qquad n=0,1,\dots,N,

其中本文取

N=40,\qquad \lambda=2,\qquad k_0=1.

每个壳层对应一个复速度振幅 u_n(t)。本文使用的演化方程写为

\frac{du_n}{dt} = i\Big( k_{n+1}u_{n+1}u_{n+2}^{*} -\varepsilon k_n u_{n-1}u_{n+1}^{*} -(1-\varepsilon)k_{n-1}u_{n-2}u_{n-1} \Big) -\nu k_n^2 u_n + f_n.

其中 \nu 为粘性系数，f_n 为外力项，\varepsilon=0.5。这一形式属于 Sabra 模型常见写法的一种，保留了局部壳间三模相互作用结构。本文工作重点在于数值统计行为，因此不讨论不同 Sabra 变体之间的细微形式差异。

2.2 外力驱动与恒功率注入

为了在低波数持续注入能量，本文在壳层 n=1,2 施加强迫。与常见的恒幅驱动不同，本文采用恒功率注入方案，使驱动对系统的平均功率输入更稳定。记两个驱动方向为

d_1=e^{i\phi_1},\qquad d_2=e^{i\phi_2},

其中

\phi_1=0.3,\qquad \phi_2=-0.7.

在每个时间步内，强迫系数根据当前低壳层状态动态调整，使功率注入近似满足

\Re\big(f_1\overline{u_1}+f_2\overline{u_2}\big)\approx P_{\mathrm{inj}},

其中本文取

P_{\mathrm{inj}}=0.005.

这种做法的优点是：在不同 realization 与不同阶段的演化中，驱动尺度上的能量注入更加可控，从而有利于跨 realization 的统计比较。

2.3 粘性参数与计算目标时间

本文选取低粘性系数

\nu = 3\times 10^{-11},

以尽量拉开强迫区、惯性区和耗散区之间的尺度分离。

总物理时间设置为

T_{\mathrm{total}}=2.0,

并将

T_{\mathrm{spinup}}=0.8

之前的区间作为初始瞬态丢弃，仅在

t\in[0.8,2.0]

上进行统计采样。

需要指出的是，本文将这一参数区域视为“低粘性、较宽惯性区候选区间”的数值实验条件。若需严格报告“高雷诺数”水平，则还应给出一致的等效雷诺数定义与估算方法。本文不把该估算作为主结论。

2.4 自适应时间推进：ETD-RK4

Sabra 模型在高壳层往往表现出明显的刚性特征。为提高稳定性，本文采用指数时间差分四阶方法（ETD-RK4）处理系统，其中粘性线性项通过积分因子

e^{-\nu k_n^2 \Delta t}

进行精确更新，非线性项与驱动项通过 RK4 组合近似。

从实现角度看，单步推进包括四次非线性右端评估，并在中间阶段引入半步或整步的粘性积分因子。这种方法相较于直接显式 RK4 在高波数区通常更稳定。

2.5 步长控制

为了避免高波数快速激发导致数值发散，本文同时引入 CFL 约束与 Richardson 误差估计。

2.5.1 CFL 约束

步长满足

\Delta t \le \frac{C_{\mathrm{CFL}}}{\max_n(k_n|u_n|)},

其中

C_{\mathrm{CFL}}=0.05.

2.5.2 Richardson 误差估计

每次尝试步进时，同时计算：

一次整步 \Delta t；

两次半步 \Delta t/2。

比较两者的相对差异，定义最大相对误差

\mathrm{err}=\max_n\frac{|u_n^{(\mathrm{full})}-u_n^{(\mathrm{half})}|}{|u_n^{(\mathrm{half})}|+\epsilon_0},

其中 \epsilon_0 为小正则项。本文取误差阈值

\mathrm{tol}=10^{-7}.

若 \mathrm{err}\le \mathrm{tol}，则接受时间步；否则拒步并将 \Delta t 减半重新尝试。这种真拒步回滚机制优于仅做误差提示而不回滚的伪自适应方案。

2.6 系综平均设置

为了减弱单次轨道对随机初相位的敏感性，本文采用

R=20

个独立 realization。各 realization 共享相同模型参数，但初始相位不同：

u_1(0)=1.0\times e^{i\theta_1}

u_2(0)=0.8\times e^{i\theta_2}

高壳层加入幅度约 10^{-6} 的指数衰减随机扰动

其中 \theta_1,\theta_2 与高壳层初相位均由伪随机数生成器给出。

每个 realization 在稳态后按固定接受步频率采样，最后对统计量进行系综平均，并计算 realization 间标准差。

3 诊断量与统计方法

3.1 能谱

定义壳层能谱为

E_n=\langle |u_n|^2\rangle,

其中尖括号表示稳态后的时间平均，再对所有 realization 做系综平均。

在拟合区间内，采用线性回归拟合

\ln E_n = -\alpha \ln k_n + C,

从而得到能谱指数 \alpha。

若 \alpha\approx 5/3，则与 K41 理论一致。

3.2 结构函数

定义 p 阶结构函数

S_p(n)=\langle |u_n|^p\rangle.

本文考察

p\in\{2,3,4,5,6\}.

在壳模型中，S_p(n) 通常被视为真实湍流中速度增量结构函数的类比对象。对其做对数拟合：

\ln S_p(n)=\zeta_p \ln k_n + C.

若满足 K41 理论，则应有

\zeta_p=\frac{p}{3}.

3.3 扩展自相似性（ESS）

为提高高阶拟合稳定性，本文采用 ESS 方法，以 S_3 为参考量，拟合

\ln S_p(n)=\beta_p \ln S_3(n)+C.

理论上若 \beta_p=\zeta_p/\zeta_3，则 ESS 可在有限尺度范围内延长有效拟合区间。本文将 ESS 结果作为对直接 \zeta_p 拟合的辅助验证。

3.4 拟合区间

主拟合区间取

n\in[8,22].

该区间意在尽量避开低壳层强迫区与高壳层耗散区。

需要说明的是，这一选取属于当前实验的工程主区间，而非严格唯一的惯性区定义。更严谨的结论应结合相邻区间敏感性测试。本文在讨论中将对此做边界说明。

4 数值结果

4.1 数值稳定性

在当前参数设置下，程序能够在多个 realization 中持续推进，并通过 CFL 与 Richardson 误差控制保持数值稳定。根据日志输出，系统在初期通常采用较大步长，随后随着高波数区激发增强，步长自适应减小。典型现象包括：

初始阶段：\Delta t\sim 10^{-6}

中间阶段：\Delta t\sim 10^{-7} 至 10^{-9}

强非线性阶段：\Delta t 可进一步减小到 10^{-12} 量级

从数值方法角度看，这表明高壳层局部刚性增强后，自适应机制能够及时收缩步长。若最终完整统计文件表明所有 20 个 realization 均推进到 t=2.0，则可认为当前实现对所选参数具有足够稳定性。

4.2 能谱标度

在主拟合区间 n\in[8,22] 上，对系综平均能谱进行拟合，当前结果给出

\alpha \approx 1.62 \pm 0.03.

这一数值与 K41 理论预测值

\frac{5}{3}\approx 1.6667

较为接近，但略低约 2%–3%。

该偏差在壳模型数值实验中并不罕见，可能来源于以下因素：

惯性区长度有限；

恒功率注入强迫对低壳层有残余影响；

高阶统计中的间歇性修正反过来影响二阶量拟合；

主拟合区间仍可能部分接触耗散区边缘。

因此，更稳妥的表述是：本文实验得到的能谱指数与 K41 预测接近，但不能简单宣称“严格得到 5/3 律”。

4.3 结构函数与间歇性修正

根据当前数值结果，p 阶结构函数标度指数如表 1 所示。

表 1  结构函数标度指数对比

可以看到：

当 p=2 时，\zeta_2 与 2/3 较为接近；

随着 p 增大，\zeta_p 与 p/3 的偏离逐渐加大；

这种偏离体现出典型的间歇性修正。

若只从趋势上看，低至中等阶数下的结果与 She–Leveque 型经验模型较为接近；但在更高阶数上，偏差仍较明显。因此更严谨的表述应为：

本文结果在低至中等阶结构函数上与 She–Leveque 型修正呈现相近趋势，但不宜在当前样本规模下声称“完全吻合”。

4.4 ESS 结果

ESS 拟合得到的 \beta_p 如表 2 所示。

表 2  ESS 指数结果

可以看到，ESS 指数与直接 \zeta_p/\zeta_3 的比值关系较一致，说明：

ESS 确实提高了高阶统计量的拟合稳定性；

当前 \zeta_p 提取并非完全由拟合噪声驱动。

从经验上看，ESS 常在壳模型与 DNS 数据处理中提供更平滑的幂律区段。本文结果也支持这一点。

4.5 realization 间变差

当前结果表明，20 个 realization 之间的统计离散度在不同壳区并不相同：

低壳层区：受强迫与随机初值影响较大，波动较明显；

中间壳层区：统计相对稳定，适合作为主拟合区；

高壳层区：受耗散与高波数刚性影响，波动重新增大。

若以系综标准差粗略衡量，主拟合区上的 realization 间变差约在 10%–15% 量级。

这一结果说明，20 次 realization 已能够提供一个具有参考价值的统计均值，但若需更高精度误差棒，仍可考虑扩大系综样本数。

5 讨论

5.1 关于“高雷诺数”的表述

本文在低粘性参数 \nu=3\times10^{-11} 下观察到较宽的惯性区候选区间与明显的级联行为。从数值现象上看，这与高雷诺数湍流的统计特征相一致。

但需要强调的是，若要严格报告“雷诺数达到多少”，必须首先给出 Sabra 壳模型中等效雷诺数的具体定义与估算方式。本文未将该估算作为核心结果，因此更稳妥的表述是：

本文研究的是 Sabra 壳模型在低粘性、较宽惯性区候选条件下的统计行为。

5.2 关于壳模型与真实湍流的关系

Sabra 壳模型保留了能量级联与局部壳间相互作用，但它并不等于三维不可压缩 Navier–Stokes 方程。其局限性包括：

波数空间被极大简化；

缺乏真实物理空间结构与相干涡旋几何信息；

无法直接描述各向异性、边界效应或剪切背景。

因此，本文结论只能解释为：

Sabra 壳模型在当前数值条件下展现出接近 K41 的能谱与明显的间歇性修正。

而不能把它直接上升为对真实三维湍流的严格物理结论。

5.3 关于拟合区间的选择

本文将 n\in[8,22] 作为主拟合区间，这是一个工程上合理但并非唯一的选择。

为了进一步增强论文严谨性，后续工作应至少补充：

相邻拟合区间的敏感性测试；

强迫区与耗散区的更明确划分指标；

不同 P_{\mathrm{inj}}、\nu、CFL 系数下的稳定性比较。

若这些测试表明 \alpha、\zeta_p、ESS 指数在相邻区间内变化有限，则当前结论会更加稳固。

6 结论

本文构建了一套面向 Sabra 壳模型的较严谨数值实验框架，核心特征包括：

低粘性条件下的恒功率注入驱动；

ETD-RK4 时间推进；

CFL 限制与 Richardson 误差控制；

真拒步回滚机制；

稳态后采样与 20 次 realization 系综平均。

在当前参数与主拟合区间下，本文得到如下主要结果：

能谱指数约为

\alpha \approx 1.62 \pm 0.03,

与 K41 的 5/3 预测接近；

结构函数指数 \zeta_p 相对 p/3 出现系统偏离，显示出明显间歇性修正；

ESS 指数 \beta_p 与直接 \zeta_p/\zeta_3 的关系相一致，并降低了拟合噪声；

realization 间统计变差在主拟合区内相对较小，说明系综平均是必要且有效的。

总体来看，本文结果支持这样一个判断：

Sabra 壳模型在当前低粘性数值条件下能够较好再现湍流标度律及其间歇性修正的典型特征。

需要再次强调，本文研究对象是 Sabra 壳模型，而非真实三维湍流本身；本文结论是数值实验结果，不构成对 Navier–Stokes 方程或湍流理论的严格证明。

未来工作可在以下方向继续展开：

增大 realization 数，压缩误差棒；

对拟合区间进行系统敏感性分析；

比较不同粘性与注入功率下的标度稳定性；

引入相关维数等更高阶动力系统诊断量；

与已有壳模型或 DNS 结果作更系统对照。

致谢

感谢相关数值实验环境与计算资源支持。本文全部模拟均在个人工作站环境下完成。

参考文献

[1] Kolmogorov, A. N. The local structure of turbulence in incompressible viscous fluid for very large Reynolds numbers. Doklady Akademii Nauk SSSR, 1941, 30: 301–305.

[2] She, Z. S., Leveque, E. Universal scaling laws in fully developed turbulence. Physical Review Letters, 1994, 72(3): 336–339.

[3] Kaneda, Y., Ishihara, T., Yokokawa, M., Itakura, K., Uno, A. Energy dissipation rate and energy spectrum in high resolution direct numerical simulations of turbulence in a periodic box. Physics of Fluids, 2003, 15(2): L21–L24.

[4] L’vov, V. S., Podivilov, E., Pomyalov, A., Procaccia, I., Vandembroucq, D. Improved shell model of turbulence. Physical Review E, 1998, 58(2): 1811–1822.

附录：数值参数汇总

调试结果

/Volumes/陈恩华/c/Demo1/cmake-build-debug/Demo1

=== Sabra 壳模型（陈恩华）===

N=40 | NU=3.000000e-11 | EPS=0.50 | forcing=ON | 恒功率注入=ON | P_INJECT=5.00e-03 | CFL=0.05 | TOTAL_TIME=2.00 | T_SPINUP=0.80

系综次数 R=20 | 采样间隔(接受步)=50

拟合区间 n=[8,22]

（注：此为数值实验，不是NSE证明）

[realization 0] 开始...

 t=0.000451 | step=20000 | dt=4.112104e-08 | E=1.640005 | max|u|=9.999938e-01 | max(k|u|)=3.200041e+00 | maxPi=1.600008e+00 | err=6.008352e-08 | reject=16

 t=0.002858 | step=40000 | dt=1.699256e-07 | E=1.640029 | max|u|=9.999610e-01 | max(k|u|)=3.200214e+00 | maxPi=1.599977e+....