在工程湍流模型中, 湍流积分尺度作为表征湍流大涡运动的特征长度^[1], 是构建湍流模型和优化工业设备的关键参数. 例如积分尺度决定了温度脉动的初始能量分布尺度, 当通过改变加热方式调整初始温度脉动的积分尺度时, 其温度脉动衰减率会随之变化, 显著影响燃烧稳定性^[2]; 在风力机尾流演化模拟中, 积分尺度的低估会错误预测湍流动能的横向扩散速率, 进而影响风电场布局优化^[3]. 此外, 积分尺度也决定着风荷载的脉动特性谱, 是工程抗风设计的关键依据^[4–6].

湍流积分尺度的演化规律不仅深刻反映了能量从大尺度向小尺度的传递过程, 更是构建湍流封闭模型的基础. 尽管其工程应用的重要性已得到广泛认可, 但积分尺度的演化规律仍存在分歧. 自Taylor^[7]首次提出积分尺度的统计定义以来, 其时间依赖性一直是湍流研究中的争议焦点. 早期理论研究中, Kolmogorov和Lumley^[8]基于能量级串的局部均衡假设, 推导出积分尺度随时间的幂律增长关系$ L\propto {t}^{2/7} $, 而Dryden^[9]通过相似性分析提出$ L\propto {t}^{1/2} $, Saffmann^[10]又通过引入Birkhoff不变量, 在低波数湍动能谱的幂律指数$ \sigma =2 $的条件下推导出$ L\propto {t}^{2/5} $, 这几种理论预测的显著分歧, 本质上源于湍动能谱幂律指数的具体值存在争议, 而该指数与空间衰减或能谱的低波数依赖性密切相关^[11]. 实验研究虽试图验证这些理论, 但受限于风洞尺寸与测量技术, 难以完全捕捉能谱低波数区的数据, 从而导致积分尺度测量产生偏差^[12,13], 这一矛盾在DNS直接数值模拟兴起后更为凸显. 因此, 尽管DNS能够精确解析高波数区的耗散过程^[14], 但周期性边界条件限制了湍动能谱中的最小可解析波数$ {k}_{L} $^[15,16], 使得大尺度涡结构的演化无法完整捕捉.

在湍流衰减动力学研究中, 耗散系数$ {C}_{\varepsilon } $作为表征能量耗散与湍动能输运均衡状态的关键参数, 其非定常特性直接影响着湍流封闭模型的构建. 而从耗散系数$ {C}_{\varepsilon } $定义式$ {C}_{\varepsilon }=\varepsilon L/{u}^{3} $($ \varepsilon $为湍动能耗散率, $ u $为均方根速度)可以看出, 积分尺度的准确测量会直接影响耗散标度. 经典Kolmogorov理论基于局部均衡假设推导出$ {C}_{\varepsilon } $的常数特性, 但该结论在有限雷诺数低波数衰减湍流中面临挑战^[18]. 实验与DNS研究表明^[19], 初始条件下大尺度运动的非均衡演化将导致$ {C}_{\varepsilon }\sim{Re}_{\lambda }^{-1} $的耗散标度律, 这区别于经典的Kolmogorov耗散规律. 然而, 由于能谱低波数缺失的影响, 数值模拟中$ {C}_{\varepsilon } $往往无法真实反映湍流系统在尺度相互作用中的非均衡特性, 这种能谱低波数的缺失掩盖了大尺度脉动对能量级串的调节作用.

不同于Kolmogorov理论的局部相似性假设^[8], George^[20]通过全局自相似性理论框架重新审视了衰减湍流的全局动力学, 认为湍动能谱$ E(k, t) $和湍流传输谱$ T(k, t) $均可以泰勒尺度$ \lambda $为唯一长度尺度进行标度化, 指出积分尺度$ L $与泰勒尺度$ \lambda $的比值$ L/\lambda $应保持恒定, 这一结论似乎与经典的Kolmogorov理论$ L/\lambda \propto {Re}_{\lambda }^{1} $^[8,21]的形成鲜明对比. DNS实测数据^[22]揭示, 在衰减湍流演化过程中, 湍动能谱峰值对应的波数$ {k}_{p} $会随时间向低波数迁移, 为了更多地捕捉到低波数信息, de Bruyn Kops和Riley^[23]提出$ {k}_{L}/{k}_{p} $应小于0.3, 但现有修正模型多局限于均衡湍流假设, 难以刻画湍动能谱低波数区的大尺度运动.

针对上述问题, 本研究基于广义von Kármán谱模型^[24]的修正框架, 将湍动能谱的大尺度结构考虑在内, 实现对积分尺度计算误差的量化评估. 通过对各项均匀同性湍流的DNS数据验证, 证实修正后的积分尺度演化遵循Saffmann理论预测的时间幂律关系$ L\propto {t}^{2/5} $, 同时揭示了湍流系统在引入大尺度结构后的非均衡特征.

本文所研究的对象是自由衰减的各向均匀同性湍流, 采用DNS直接数值模拟的方法生成数据, 对不可压缩纳维-斯托克斯(Navier-Stokes)方程在全时空尺度上进行数值求解.

其中$ \boldsymbol{u} $是无散度速度场($ \nabla \boldsymbol{u}=0 $); $ \boldsymbol{P} $是压力; $ \rho $是流体密度; $ \nu $是运动黏度. DNS算例采用Rogallo^[25]提出的伪谱法求解, 并用半隐式法处理黏性项, 时间方向上采用四阶显式龙格-库塔(Runge-Kutta)方法进行积分. 本研究计算区域为一个边长为$ 4{\pi } $的周期性立方体, 计算网格数为$ {384}^{3} $, 空间分辨率$ {k}_{\mathrm{m}\mathrm{a}\mathrm{x}}\eta =1.65 $($ {k}_{\mathrm{m}\mathrm{a}\mathrm{x}} $为最大波数, $ \eta $为耗散尺度). 该算例从一个完全发展的自由衰减湍流场开始计算, 且该自由衰减流场是由随机初始场计算得到的^[26]. 随机场的能谱分布与Comte-Bellot和Corrsin^[12]的实验能谱一致. 模拟时间t的区间为[0.24350, 0.30325], 该区间内的湍流已经经历了长时间的衰减, 符合各向同性湍流的均衡特征.

在均匀各向同性湍流研究中, 积分尺度与能量的准确计算依赖于能谱在全部波数范围内的完整描述. 然而, 受限于实验测量或数值模拟的物理条件(如计算域尺寸或网格分辨率), 实际可获取的能谱往往缺失部分低波数信息. 这一问题对积分尺度的测量影响尤为显著.

湍动能$ E $与能谱$ E\left(k\right) $之间的关系为

其中$ u $为均方根脉动速度; 尖括号表示求系综平均; $ k $为波数. 谱空间中湍流积分尺度求解如下:

但在实际DNS模拟实验中, 我们无法得到波数低于$ {k}_{L} $的波数值, 因此实际值为

其中$ {u}_{\mathrm{m}} $, $ {L}_{\mathrm{m}} $为未修正的均方根脉动速度和积分尺度(公式中下标$ \mathrm{m} $代表未修正值, 未标注下标的$ u $和$ L $为真实值), 本研究中DNS可取得的最小波数值$ {k}_{L}=0.52 $. 由于低波数段部分缺失, 基于波数范围$ [{k}_{L}, +\infty ) $的积分计算将导致积分尺度实际值$ {L}_{\mathrm{m}} $显著低于真实值$ L $.

图1展示了不同时间下能谱$ E\left(k\right) $随波数的变化, 从图中可以看出本研究算例中最小可解析波数与能谱峰值波数之比$ {k}_{L}/{k}_{p}=1 $, 该值大于de Bruyn Kops和Riley^[23]提出的$ {k}_{L}/{k}_{p}=0.3 $, 这意味着数值模拟未能完全解析湍流大尺度结构的动态演化, 导致大涡运动信息部分丢失. 事实上, 随着能谱峰值波数$ {k}_{p} $向低波数方向迁移^[22], $ {k}_{L}/{k}_{p} $比值呈现单调递增趋势, 使得低波数缺失现象愈发严重, 此时若仍采用缺失后的能谱计算积分尺度, 误差会随时间的推移被放大.

本研究采用von Kármán谱模型^[24]的广义形式作为理论框架, 其函数形式可表示为

其中$ B $为调整幅值的无量纲归一化系数; $ x=k/{k}_{\mathrm{e}} $, $ {k}_{\mathrm{e}} $为特征波数, 用于划分谱模型中低波数到高波数的过渡. 需要指出的是, 该谱模型假设的对象是低波数区, 能谱按$ {k}^{4} $增长, 而惯性子区按$ {k}^{-5/3} $衰减. 在有限雷诺数下能谱惯性子区的行为会轻微偏离经典$ {k}^{-5/3} $律^[17], 因此Gamard和George^[27]在保留(6)式框架的基础上, 将有限雷诺数和低波数缺失的影响考虑在内, 通过引入湍动能谱幂律指数$ p $, 创建了一个涵盖所有波数范围的广义模型:

其中系数$ {C}_{p} $和特征波数$ {k}_{e} $通过能量和积分尺度的积分约束确定^[28,29], 参数$ p $根据George^[20]的相似关系确定. 基于此, 我们通过修正的能谱模型来补偿数值模拟中的低波数缺失效应. 具体实施时, 将积分下限为$ {k}_{L} $的湍动能积分值与积分下限为0的理论积分值作比, 再进行简单的变量替换并转换为不完全贝塔函数:

其中变量$ z $由以下公式定义:

$ {I}_{z} $为不完全贝塔函数:

其中贝塔函数:

同理我们将(7)式谱模型进行积分尺度积分作比, 可以得到(12)式:

要对$ L $进行迭代求解, 需要将变量$ z $用准确的积分尺度$ L $表示, 在Wang和George^[28]的附录A中:

为了将$ L $与$ z $建立联系, 我们进行$ \dfrac{{k}_{L}}{{k}_{p}}=\dfrac{{k}_{L}L}{{k}_{p}L} $的参数转换. (13)式可以看出, 当参数$ p $确定时, $ {k}_{p}L $为常数. (12)式左侧的$ {L}_{{\mathrm{m}}} $为已知的实际值, 右侧的变量$ z $依赖于$ L $, 因此方程可以表示为$ {L}_{\mathrm{m}}/L=f\left(L\right) $, 要想得到修正后的$ L $需要提供$ {u}_{\mathrm{m}}^{2} $, $ {L}_{\mathrm{m}} $和最小波数$ {k}_{L} $的值, 然后进行隐式迭代得到真实值.

需特别指出, 以上公式均含有的参数是指数$ p $, 因此正确选择$ p $对修正的准确性至关重要. 图2(a)和图2(b)分别展示了参数$ p $和低波数缺失程度$ {k}_{L}/{k}_{p} $对修正后能量比$ {u}_{\mathrm{m}}^{2}/{u}^{2} $及积分尺度比$ {L}_{\mathrm{m}}/L $的影响, 横坐标$ {k}_{L}/{k}_{p} $作为独立变量进行赋值.

从图2(a)可见, $ {u}_{\mathrm{m}}^{2}/{u}^{2} $比对$ p $值的敏感性较弱. 例如当$ {k}_{L}/{k}_{p}=0.5 $时, $ p=1 $和$ p=2 $对应的能量损失差异不足5%. 而图2(b)中积分尺度比$ {L}_{\mathrm{m}}/L $对$ p $和$ {k}_{L}/{k}_{p} $的综合影响显著敏感. 当$ p=2 $且$ {k}_{L}/ {k}_{p}=0.5 $时, $ {L}_{\mathrm{m}}/L $的偏差为16.91%, 这是由于积分尺度计算中低波数区域$ E\left(k\right)/k $的权重更高((3)式), 而低波数缺失直接削弱了其对积分尺度的贡献. 正如第3节提到的, 随着$ {k}_{L}/{k}_{p} $增大, 能谱峰值会向更低波数迁移, 不同$ p $值对应的$ {L}_{\mathrm{m}}/L $曲线呈现显著发散. 例如, 当$ {k}_{L}/{k}_{p}=1.0 $且$ p=2 $时, 积分尺度偏差差异可达39.7%, 这一显著差异揭示参数$ p $的选择需与动力学约束严格匹配. George^[20]首次证明, 当所研究数据遵循$ {u}^{2}\propto {t}^{n} $的幂律衰减特性时, 存在$ p=-2 n-1 $的关系, 其中$ n $为湍动能衰减指数, 该关联确保了修正过程中指数$ n $与积分尺度演化的自洽性.

在George^[20]的自相似性假设下, 泰勒尺度$ \lambda $为单一长度标度并在全波数范围内有效, 而泰勒尺度的平方在该假设下将呈线性增长^[10,20,30]:

其中$ {t}_{0} $为虚拟原点, 选择$ {t}_{0} $的作用在于调整时间轴的起点, 在模拟中, 实际的时间起点可能与理论模型存在偏差, $ {t}_{0} $可以通过平移时间轴来消除这种偏差. 本研究在确定(14)式中$ n $的值后得到$ {t}_{0}= +0.00387 $, $ n $值详见第5节. 泰勒尺度$ \lambda $的定义:

同理(8)式的推导思路, 我们得到耗散率$ \varepsilon $的修正公式:

根据(16)式要求, $ p $必须满足$ p < 5/3 $, 在修正过程中$ p $接近$ 5/3 $时(16)式右侧趋于1, 即耗散率受低波数的影响可以忽略, 我们的DNS计算结果也表明, 利用谱方法$ {\varepsilon }_{\mathrm{m}}=2\nu \displaystyle\int_{{k}_{L}}^{\infty}{k}^{2}E\left(k\right)\mathrm{d}k $与梯度法$ {\varepsilon }_{\mathrm{m}}=\nu \left\langle{\dfrac{\partial {u}_{i}}{\partial {x}_{k}}\dfrac{\partial {u}_{i}}{\partial {x}_{k}}}\right\rangle $得到的计算误差小于$ {10}^{-10} $, 所以当$ p\geqslant 5/3 $时取$ {\varepsilon }_{\mathrm{m}}=\varepsilon $. 根据(15)式和(16)式可以得到$ \lambda $的修正值, 通过对(14)式求导, 可以更方便地得到指数$ n $:

因此, 我们利用关系式$ p=-2 n-1 $可以得到修正后的$ p $值. 当$ p $的初始假设值与其修正后的值相差小于$ 1{0}^{-3} $时, 迭代过程截止.

(14)式中$ n $值的确定依赖于泰勒尺度的动态演化特征. 图3(a)对比了修正前后泰勒尺度平方的时间导数$ \mathrm{d}{\lambda }^{2}/\nu \mathrm{d}t $的演化趋势, 在$ 0.29725 <t < 0.30275 $的区域内, 导数曲线趋于水平, 对应指数$ n=-1.40 $. 此时修正前后导数值的拟合误差小于5%, 有效验证了(17)式的理论预测. 泰勒尺度$ \lambda $的修正依赖于均方根脉动速度$ u $与耗散率$ \varepsilon $的共同作用, 图2(a)中可以看出, 修正前后$ {u}^{2} $的误差小于5%, 而耗散率$ \varepsilon $的修正贡献在参数$ p\geqslant 5/3 $时可以忽略, 这表明$ \lambda $的修正主要源于$ u $的微小调整, 这在图3(b)的$ {\lambda }^{2} $演化曲线中也得到体现. 值得注意的是, 当时间$ t > 0.27850 $时, 修正后的$ {\lambda }^{2} $曲线与(14)式理论预测的拟合误差小于0.4%, 存在$ {\lambda }^{2}\sim{t}^{1} $, 且实际的$ {\lambda }_{\mathrm{m}}^{2} $曲线也接近该线性关系. 湍动能随时间的演化如图4所示, 可以看出修正前后的数据均遵循$ {u}^{2}\propto {t}^{-1.40} $的标度关系, 这与Saffmann和Oberlack^[11,12]在指数$ \sigma =2 $下得到的$ \overline{u}' \propto {t}^{-6/5} $框架吻合, 其中$ \overline{u}' $为湍动能的统计平均值. 图5为修正前后积分尺度$ L $随时间的演化, 通过参数$ p $的迭代优化, 模型恢复了低波数区$ E\left(k\right)/k $对积分尺度的贡献, 使得图5中$ L $呈现出Saffmann和Oberlack^[11,12]的理论演化特性$ L\propto {t}^{2/5} $. 需要指出的是, 本研究受限于低波数的缺失($ {k}_{L}/{k}_{p}=1 $), 导致了对积分尺度的低估, 这本质上反映了低波数信息缺失对能量级串的影响.

此外, 本文还将谱空间求得的积分尺度与物 理空间利用速度相关函数计算得到的积分尺度$ \varLambda $进行了对比. 在物理空间, 积分尺度的定义为$ \varLambda =\displaystyle\int_{0}^{{r}_{0}}{\rho }_{uu}\left(r\right)\mathrm{d}r $, 其中, $ {\rho }_{uu}(r) = \dfrac{\langle{u (x)u(x+r)}\rangle}{\langle{{u}^{2}}\rangle}$是纵向脉动速度相关, 积分上限$ {r}_{0} $取$ {\rho }_{uu}\left(r\right) $第一次为零值的距离. 图6展示了两种计算方法(物理空间与谱空间)计算得到的积分尺度比值随时间的演化. 可以看到, 物理空间积分尺度与修正后谱空间积分尺度的比值($ \varLambda /L $)小于其与未修正谱空间积分尺度的比值($ \varLambda /{L}_{\mathrm{m}} $), 且两比值都大于1. 这说明在考虑了低波数的能谱信息之后, 谱空间积分尺度明显变大, 但其数值仍小于利用速度相关函数计算得到的积分尺度. 造成这一现象的主要原因有以下两点: 1)本文的计算域不够大, 有研究表明计算域的长度至少在8倍积分尺度$ \varLambda $的情况下^[1,17], 利用速度相关函数计算得到的积分尺度才能有较好的精度, 在本文中计算域的长度约为积分尺度$ \varLambda $的6倍; 2)本文虽然对谱空间进行了修正, 但本文算例中网格数较少, 因此谱空间中的波数也较少, 这会导致谱空间下积分尺度的数值仍然偏小(在与本文算例计算域长度相同, 但网格数为$ {640}^{3} $的另一个算例中, $ \varLambda /L $比本文更小, 约为1.54). 理论上, 在网格数足够的情况下, 物理空间积分尺度与修正后谱空间积分尺度的比值($ \varLambda /L $)接近于1.

与格栅湍流中非均衡性的研究类似, 我们通过考察耗散系数$ {C}_{\varepsilon } $随时间的演化来判断整体能量传输与耗散之间是否均衡. 自由衰减湍流中耗散系数$ {C}_{\varepsilon } $随时间的演化特性如图7所示, 修正前的$ {C}_{\varepsilon } $稳定, 符合统计均衡态特征, 这一结果在高雷诺数下的自由衰减湍流中被普遍接受. 值得注意的是, 修正前的$ {C}_{\varepsilon } $维持小于1的恒定值, 这归因于能谱低波数缺失导致的积分尺度低估. 而修正后的结果呈现$ {C}_{\varepsilon }\sim{Re}_{\lambda }^{-1} $的标度关系, 表现出能量过度耗散特征, 这与Goto和Vassilicos^[19]在湍流衰减初期阶段(未达到统计平衡)得到的高雷诺数实验结果一致. Bos和Rubinstein^[31]的理论推导进一步支持了这种非均衡耗散标度的合理性. 耗散系数出现差异的原因是: 尽管能谱$ E(k, t) $在中间波数区呈现$ {k}^{-5/3} $局部平稳特征, 但系统整体未达到平稳, 谱模型(7)式纳入大尺度动力学后, 这会导致能量输运与耗散的动态失衡, 进而得到非均衡统计特性. Steiros^[32]关于各项同性衰减湍流的研究也表明, 有限雷诺数下大尺度区的非均衡特性不能保证Kolmogorov耗散标度$ {C}_{\varepsilon }=1 $. DNS和EDQNM模拟^[33,34]证实: 当演化超过大尺度结构主导的能量集中区域时(表现为$ kL\approx 2-10 $区间的谱峰值), 黏性耗散逐渐消退, 系统便从新耗散尺度$ {C}_{\varepsilon }\sim {Re}_{\lambda }^{-1} $向经典耗散尺度$ {C}_{\varepsilon }=1 $转变. 因此随湍流的演化, 低波数缺失会进一步加剧, 使得Kolmogorov尺度占据主导地位, 湍流系统趋于均衡.

两个特征尺度的比值$ L/\lambda $与$ {Re}_{\lambda } $的关系也反映了这种非均衡特性, 如图8所示. 修正前的特征尺度比值存在$ L/\lambda \sim{Re}^{1}_{\lambda} $关系, 这符合传统的Kolmogorov局部均衡理论预测^[8]. 修正后的$ L/\lambda $呈现趋稳现象, 这源于在初始大尺度的扰动下, 耗散过程持续遵循新的耗散规律^[35]$ \varepsilon \propto {u}^{2}/{L}^{2} $, 耗散系数呈现文献中报道的非均衡标度$ {{C}_{\varepsilon }\sim Re}_{\lambda }^{-1} $, 这与格栅湍流的研究结果一致^[36–39].

本研究基于广义von Kármán谱模型^[24], 系统地分析了低波数缺失对各向均匀同性湍流统计参数的影响. 受周期性边界条件和计算域尺寸限制, 实际DNS直接数值模拟方法因未能完全解析湍动能谱低波数区的大尺度运动, 导致积分尺度实际值$ {L}_{\mathrm{m}} $较理论真实值$ L $显著偏低, 且误差随能谱峰值向低波数迁移而显著放大. 我们通过引入湍动能谱修正模型, 恢复了低波数区域$ E\left(k\right)/k $的贡献. 修正后的积分尺度$ L $随时间的变化规律呈现$ L\propto {t}^{2/5} $, 这符合Saffmann理论预测的幂律关系; 修正后耗散系数$ {C}_{\varepsilon } $的演化满足湍流非均衡耗散规律$ {C}_{\varepsilon }\sim{Re}_{\lambda }^{-1} $, 这与Goto和Vassilicos^[19]在湍流衰减初期阶段得到的高雷诺数实验结果一致. 研究通过重构数值模拟中缺失的能谱大尺度信息, 揭示了湍流系统的非平衡态特性, 证实了大尺度涡结构对湍流能量耗散过程具有显著调控作用. 进一步分析表明, 在有限雷诺数或者受初始条件影响较大的湍流流动中, 大尺度结构对流动的影响显著, 湍流无法在全尺度实现均衡. 这一发现突破了Kolmogorov的局部平衡假设, 将为湍流多尺度相互作用中大尺度结构的统计建模提供新范式.