0 引言

高雷诺壁湍流的研究近年来成为湍流研究的热点。在高雷诺数下，流场内流动结构的尺度层级更加丰富，不同尺度间相互作用的形式更加复杂，并且出现了不同于低雷诺数边界层的统计特性以及物理结构，比如大/超大尺度结构。这种外区的大/超大尺度结构会影响湍流边界层的湍动能的生成和分布，甚至会对近壁区的小尺度结构所对应的雷诺应力和湍动能产生影响。为了量化湍流边界层多尺度结构的运动学和动力学特性以及它们之间的相互作用，需要将大小尺度结构进行合理的分离^[1]。在高雷诺数的情况下，耗散尺度和含能尺度充分分离，因此分离尺度相对比较明显。常见的方法是通过设定某个截断波长，对谱空间进行截断滤波。例如，Mathis等^[2]推荐边界层厚度δ作为流向截断尺度，即λ_x/δ=1，用于研究大尺度结构对小尺度脉动的振幅调制。这一标准在湍流边界层的研究中被广泛采用^[3-5]。此外，Zhang等^[6]提出了分别在流向、展向和时间方向上的内尺度截断尺度λ_x⁺=1000π, λ_z⁺=125π和λ_t⁺=260，这在某种意义上是“最优”的截断尺度，因为所分离出来的小尺度湍流脉动接近于准定常准均匀假设(Qausi-steady qausi-homogeneous, QSQH)，即小尺度湍流脉动与大尺度脉动产生的壁面剪切之间相互独立。

尽管谱分析在尺度分离中有广泛的应用，但它也有一定的局限性。首先，傅立叶变换隐含地假设被分析信号达到统计稳态。对于通过单点测量(例如热线风速仪)获得的速度脉动的时间序列，可以通过增加采样的持续时间来满足这一条件。然而，对于像粒子图像测速(PIV)这样的场测量，视野范围(FOV)通常受到测试仪器的限制不可能太大，因此瞬时脉动速度在空间上的分布，如流向速度u(x)，往往不能满足统计稳态条件。这使得针对短信号的傅里叶分析无法准确解析大尺度(低波数)处的能谱。其次，在中低Re的情况下，尺度分离并不明显^{[4, 7]}，这使得选择合适的分离尺度相当困难。另外，人工给定一个大、小尺度之间的固定的截断阈值，在物理上与湍流结构的时变特性相矛盾。

经验模态分解(Empirical Mode Decomposition, EMD)^[8]及其变体是另外一种尺度分离选择，其本身具有处理非稳态信号的能力。EMD是通过一系列的筛选(sifting)过程将时间序列信号分解为有限数量的本征模态函数(Intrinsic Mode Functions, IMFs，在下文中也简称为模态)，而不依赖于预先确定的模态形式, 即可以自适应地捕捉到信号的局部尺度。在湍流尺度分析中得到了较为广泛的应用^[9-10]，Cheng^[11]等最近对通过DNS模拟的槽道湍流速度场进行了二维EMD (B-EMD)分析, 得到了流场中的条带结构和流向涡结构，并证明这些结构的尺度沿法向高度自相似增长的特性，这与传统的附着涡假设所得到的结论是一致的。

由Dragomiretskiy和Zosso^[12]提出的变分模态分解(VMD)是传统EMD的最新升级。VMD方法的基本思想是使所分解出的模态应该具有最紧致的带宽。在此基础上，VMD用严格的数学定义明确地构造了一个凸优化问题。由于最小带宽规则的限制，不同VMD模态在尺度空间中相互重叠最小，这避免了传统EMD的模态掺混问题。这种优势使得VMD在信号处理方面具有巨大的潜力。在此基础上，Wang^[13-14]的工作对VMD方法在多维多分量信号分解上的应用进行了探索，提出了多分量准二维VMD(MC-QBVMD)方法并用于湍流边界层拟序结构的分析，验证了这种方法在湍流多尺度分析中的有效性。然而这种自适应尺度分解方法得到的模态所对应的物理结构的本质仍不清楚。

本文将利用MC-QBVMD方法，对直接数值模拟(DNS)得到的湍流边界层进行分解。重点对分解得到的各阶模态中的拟序结构的尺度和几何特性进行分析，揭示其物理本质。本文分为以下几个部分，第1节将介绍变分模态分解方法及其衍生的MC-QBVMD方法；第2节将介绍所使用的DNS边界层数据；第3节将对分解的到的结果进行分析讨论；第4节给出结论。

1 变分模态分解方法 1.1 一维变分模态分解

在这里我们简要的介绍一维VMD算法的核心思想和实现方法，以方便后面的讨论。关于算法本身的更多细节可以参考文献[12]。

VMD算法的核心是将一列时序信号f(t)分解成一系列本征模态函数(Intrinsic mode function, IMF)的线性叠加，即：

$ f(t) = \sum\limits_{k = 1}^K {{u_k}} $

(1)

其中u_k为IMF，其特征是极值点的数量与过零点的数量相等或仅相差一个。模态数量K需要提前指定。VMD算法的核心目的是要使得各阶模态在尺度上充分分离，这也就引出了VMD算法的核心优化函数，即：单个模态的带宽要最窄。这一最窄带宽准则可以表述成以下的优化函数，即：

$ \begin{array}{l} \mathop {{\mathop{\rm argmin}\nolimits} }\limits_{{u_k}, {\omega _k}} \left\{ {\sum\limits_k {\left\| {{\partial _t}\left[ {\left( {\delta (t) + \frac{{\rm{i}}}{{\pi t}}} \right){u_k}(t)} \right]{{\rm{e}}^{ - {\rm{i}}{\omega _k}t}}} \right\|_2^2} } \right\}{\rm{ }}\\ {\rm{s}}{\rm{.t}}{\rm{. }}\quad f = \sum\limits_k {{u_k}} \end{array} $

(2)

上式目标函数中的L2范数‖ ·‖²是各阶模态的带宽的估计值。其中$ \left(\delta(t)+\frac{\mathrm{i}}{\pi t}\right) u_{k}(t)$实际上是模态u_k的分析函数a_k(t)：

$ a_{k}(t)=u_{k}(t)+\mathrm{i} H u_{k}(t) $

(3)

公式(3)中，H为希尔伯特变换(Hilbert Transform)。分析函数a_k(t)的傅里叶变换是一个只在正频率有幅值的函数。${{\rm{e}}^{ - {\rm{i}}\omega kt}} $的作用是将整个信号的频域上移动ω_k，其中ω_k可以认为是模态u_k的中心频率，或者载波频率。将经过平移的分析信号关于t求偏导数，再求L2范数，相当于求a_k(t)的频谱关于ω_k的转动惯量，因此可以反应u_k的带宽。由此可见，公式(2)所表达的优化函数的含义是：求一个原始信号f的模态的集合{u_k}和中心频率ω_k的组合，使得所有模态带宽相加之和最小。对该目标函数的约束条件是，所得全部模态得线性叠加能够完全还原原始信号。

求解公式(2)中得优化问题即可得到原始信号f的模态，下文称为VMD模态。可以使用变向乘子法(Alternate direction method of multipliers, ADMM)来求解这一优化问题。为了简明起见不再给出推导求解的过程，仅仅给出Dragomiretskiy等^[12]给出的迭代求解公式，如下所示：

$ \hat u_k^{n + 1} = \left( {{{\hat f}_m} - \sum\limits_{j \ne k} {{{\hat u}_j}} + \frac{{\hat \lambda }}{2}} \right)\frac{1}{{1 + 2\alpha {{\left( {\omega - {\omega _k}} \right)}^2}}} $

(4)

$ \omega _k^{n + 1} = \frac{{\int_0^\infty \omega {{\left| {{{\hat u}_k}} \right|}^2}{\rm{d}}\omega }}{{\int_0^\infty {{{\left| {{{\hat u}_k}} \right|}^2}} {\rm{d}}\omega }} $

(5)

其中尖帽符号$ {\hat \cdot }$代表傅里叶变换，λ为拉格朗日乘子，α是用以加速迭代收敛的权重系数。λ和α需要提前设定。设定分解模态的数量和残差的阈值之后即可不断的迭代更新公式(4)和(5)，直至结果收敛。从公式(4)和(5)中可以看到，模态u_k在每一步都迭代更新成其他模态残差的Wiener滤波，而中心频率ω_k则更新成$ \hat{u}_{k}$的重心。

1.2 多分量准二维变分模态分解

瞬时湍流速度场是三维三分量的，但是湍流边界层的流场本身的特殊性，比如脉动速度沿法向的分布具有不均匀性、流场中的拟序结构在不同方向上的尺度差别巨大等等导致同时在两个以上的维度上进行分解非常困难。比如二维FFT无法对x-y平面的速度场进行分解，因为法向不满足周期性条件；二维EMD分解x-z平面流场时往往无法有效的提取其中的条带状结构等等^[14]。

针对这个问题, Wang等^[13]提出了多分量准二维变分模态分解(MC-QBVMD)算法来分解二维多分量数据。这个方法的思路是，对于一个二维多分量的数据信号$ f(x, y)=\left\{f_{m}(x, y), m=1: M\right\}$，假设x方向为分解方向，y方向为辅助方向，M为分量个数，则f(x, y)可以看成一系列x方向上的空间信号在y方向上的组合。假设由于拟序结构的存在，同一阶模态中主方向上的尺度沿辅助方向没有变化。因此在其迭代求解的过程中，可以先在主方向上进行尺度分解，然后在辅助方向上对不同位置处的中心尺度进行重新平衡，使得每阶模态的主方向上的尺度在辅助方向上保持一致。对于不同的速度分量，MC-QBVMD要求不同分量之间的尺度保持一致，这是因为对于模态u_k中某个拟序结构引起的脉动，其在不同分量上的投影的尺度应该是一致的。

MC-QBVMD的目标函数可以表示为：

$ \begin{array}{l} \mathop {\min }\limits_{{u_{m, k}}, {\omega _k}} \left\{ {\sum\limits_{m, k} {\left\| {{\partial _x}\left[ {\left( {\delta (x) + \frac{{\rm{i}}}{{\pi x}}} \right){u_{m, k}}(x, y)} \right]{{\rm{e}}^{ - {\rm{i}}{\omega _k}x}}} \right\|_F^2} } \right\}\\ {\rm{ s}}{\rm{.t}}{\rm{. }}\quad {f_m} = \sum\limits_k {{u_{m, k}}} , \quad \forall m \in [1, M] \end{array} $

(6)

上式假设x为主分解方向。由于f为二维模态，因此目标函数中的L2范数被替换成Frobenius范数‖ ·‖_F来估计二维模态的带宽。对于一个二维矩阵u_ij，对比公式(2)与MC-QBVMD的目标函数可以发现，MC-QBVMD的目标函数相当于不同分量和不同y处的一维空间信号的带宽之和。

由于MC-QBVMD的目标函数相当于不同空间位置处的多分量一维信号的带宽的线性叠加，因此其求解的步骤与VMD非常相似。使用ADMM求解式(6)，可以得到迭代解为：

$ \begin{array}{*{20}{c}} {\hat u_{m \cdot k}^{n + 1}(\omega , y) = }\\ {\frac{{\left( {{{\hat f}_m}(\omega , y) - \sum\limits_{j \ne k} {{{\hat u}_{m \cdot j}}} (\omega , y) + \hat \lambda (\omega , y)/2} \right)}}{{1 + 2\alpha {{\left( {\omega - {\omega _k}} \right)}^2}}}} \end{array} $

(7)

$ \omega _k^{n + 1} = \frac{{\int\limits_{D(y)} {} \int_0^\infty \omega \sum\limits_m {\left( {{{\left| {{{\hat u}_{m, k}}(\omega , y)} \right|}^2}} \right)} {\rm{d}}\omega {\rm{d}}y}}{{\int\limits_{D(y)} {} \int_0^\infty {\sum\limits_m {\left( {{{\left| {{{\hat u}_{m, k}}(\omega , y)} \right|}^2}} \right)} } {\rm{d}}\omega {\rm{d}}y}} $

(8)

公式(7)说明在MC-QBVMD中，u_k的更新步相当于对f(x, y)在每个y位置和每个分量进行一维的VMD分解，而ω_k的更新步不但在各个分量之间进行平均，同时在不同的y位置进行平衡，使得每一阶一维模态u_{m, k}在整个y域上具有相同的中心频率。关于MC-QBVMD算法的更多细节可参照Wang等^[13]的工作。

2 DNS湍流边界层数据

我们将在后文中分析马德里理工大学(Universidad Politécnica de Madrid) Javier Jimenez教授课题组公开的DNS湍流边界层数据(https://torroja.dmt.upm.es/turbdata/blayers)。该数据为体发展边界层流动，雷诺数跨度Re_τ=400~2000，模拟了光滑平板之上的一个六面体内的速度场，展向采用的是周期性边界条件，流向和法向为非周期性边界条件。在流向和法向采用交错网格三点紧致差分方法。展向上使用傅里叶谱进行计算，并使用2/3律来消除混淆误差。求解压力使用的Poisson方程中采用的是二阶有限差分格式。时间推进使用的是半隐式三阶Runge-Kutta方法。关于这两套DNS数据的更多细节可以参考原始文献[15]。

在当前工作中，我们选取6个特征雷诺数下的流向-法向平面作为研究对象，流向上选取各个雷诺数对应的位置上下游各6δ(δ为特征雷诺数当地的边界层厚度)的范围，法向上选取y≤δ的区域，基本涵盖了边界层中所有流向尺度的结构。在当前选择的流向范围内，雷诺数的变化在±3%以内，因此后面将采用流向范围中心位置处的边界层参数作为全局参数进行无量纲。各个工况的流向平均速度型U和流向脉动速度型σ_u如图 1所示。不同工况下的边界层参数如表 1所示。

3 变分模态分解结果 3.1 瞬时流场结构分解

使用MC-QBVMD对湍流边界层中流向-法向(x-y)平面内的多分量瞬时速度场进行尺度分离。

需要说明的是，湍流边界层中脉动速度的统计量在平行壁面平面内是准均匀的，在垂直壁面方向是非均匀的。这种情况适合应用准二维分解方法进行分解，分解的主方向是流向，辅助方向是法向。而传统的二维分解方法，例如二维FFT、或者二维EMD，由于在两个方向上均进行分解，所以并不适用于流向-法向平面内的二维数据集。参考Wang的工作^[13]，模态的阶数定为K=2, 罚系数定为α=100，λ=0。

图 2展示了湍流边界层(Re_τ=1750)某单帧速度场以及MC-VMD分解得到的两阶模态。这里我们仅仅展示原视野范围的一小部分，以便更加清晰地显示模态中的流动结构。第一阶模态(图 2(b))很好的捕捉到了原始流场中的大尺度结构，即沿流向排布的交替排列的高低速流体，其流向尺度达到了边界层厚度的量级。这种结构与前人观察到的大尺度结构在x-y平面的是一致的^[16-17]，其特征都是向流向倾斜，一端与壁面接触，另一端向对数区延伸。由于u、v两个分量是同时分解的，因此可以看到这些附着于壁面的结构与大尺度的扫掠喷射事件密切相关：即高速区对应大尺度的Q4事件，低速区对应大尺度的Q2事件。

第二阶模态(图 2(c))显示了瞬时的小尺度模态，其特征尺度比第一阶模态小一个量级，可用内尺度进行标度(λ_x⁺~300)。小尺度模态中可以观察到近壁面(y⁺ < 50)频繁发生的小尺度喷射和扫掠事件。这些事件在流向交替排列，并且与大尺度结构一样与流向方向存在倾角。这些小尺度结构可能与湍流边界层内区的自循环机制有关^[18]。

3.2 大小尺度模态结构特征

在观察MC-QBVMD模态的瞬态结构(如图 2)时，我们发现两阶模态均表征出相类似的结构，即：从壁面倾斜延伸向外的、沿流向交替排列的高低速条状结构。为了进一步揭示这种结构的形态学特征，我们统计了大/小尺度模态中不同速度分量的倾斜结构的结构倾角随高度的变化。计算结构倾角的方法如下，首先计算各个速度分量在不同法向高度处的两点相关，即：

$ \begin{array}{c} R_{u u}\left(r_{x}, r_{y} ; y\right)= \\ \frac{\left\langle u\left(x+r_{x}, y+r_{y}, t\right) u(x, y, t)\right\rangle_{x, t}}{\sigma_{u}^{2}(y)} \end{array} $

(9)

其中:u代表VMD瞬时模态的三个脉动速度分量(u, v, w)中的一个，〈·〉代表在x、t两个维度上做系综平均。然后求R_uu峰值在不同r_y处的最大值对应的流向偏移量r_x^max。通过对y高度附近的(r_x^max, r_y)进行线性拟合即可得到当地结构的倾斜角θ，表示倾斜结构与壁面的夹角。在当前的计算中，我们选取位置y及其附近的21点作为拟合所用的点集。一个典型高度处(y⁺=30)的大/小尺度MC-QBVMD模态的两点相关云图R_uu、R_vv和R_ww如图 3所示，图 3中也显示了检测高度附近的(r_x^max, r_y)线性拟合得到的斜率。图 4给出大/小尺度VMD模态的不同速度分量的结构倾角随高度的变化，为便于比较，原始流场的结构倾角随高度的变化也一并给出。首先，不同雷诺数下原始流场以及MC-QBVMD模态的结构倾角的变化趋势基本一致。但是大小尺度模态的几何倾角随高度的变化有明显的区别，不同分量的结构倾角也各不相同。对于u分量，原始流场的结构倾角随高度增加逐渐增加，在y⁺≈200之后稳定在60°左右。大尺度u分量的结构倾角与原始流场基本相同，区别是在y⁺≈200之后，其倾角继续缓慢增长。小尺度u分量的倾角则始终高于大尺度模态的u, 最高可以到达80°(y⁺≈200)，之后缓慢减小。v分量的倾角远远高于u分量，原始的v分量在y⁺≈60处就达到了90°倾角，即其结构与壁面垂直。小尺度v分量的结构倾角的变化趋势与原始流场的v基本一致。大尺度模态v分量倾角的增长趋势较慢，在y⁺≈200处达到90°。w分量的倾角相比其他两个分量的倾角要小，说明其结构更加向壁面倾斜。原始流场中的w分量的倾角在y⁺≈200稳定于50°左右。而大尺度w分量和小尺度w分量的结构倾角均随法向高度的增加而持续增加，但前者结构倾角的数值始终小于后者。

对图 4显示出的结果可以做以下三点更加深入的讨论。第一，仔细分析结构倾角的雷诺数影响可以发现，原始速度场的结构倾角呈现出较为明显的雷诺数依赖性，但在将原始速度场分解后，大/小尺度模态的结构倾角的雷诺数依赖性都有不同程度的减弱，尤其是在大尺度u和w分量(图 4(a，c))上。这似乎暗示，对于平行壁面的速度分量，其倾斜的大尺度结构的几何特征在所研究的雷诺数范围内具有“普适性”，该推测需要在更高的雷诺数下去做进一步的验证。第二，从象限分析的角度出发，v分量的结构倾角达到90°这一现象，可以理解为Q2/Q4事件对v分量湍流脉动的贡献达到均衡。图 4(b)显示，相比于大尺度而言，小尺度v分量的结构倾角都在一个更低的流层上收敛至90°，且不同雷诺数之间的差异很小，这与小尺度的内区自维持循环只限于近壁流区、且具有雷诺数无关性的传统观点相吻合。另一方面，大尺度v分量的结构倾角收敛至90°的法向高度，似乎随着雷诺数的增大而增大。这也许可以说明，以Q2喷射事件为特征的大尺度低速结构是更为常见的近壁大尺度结构；或者可以说明，壁面条件对大尺度v分量的影响延伸至很高的流区，且与雷诺数相关。第三，附着涡模型(attached eddy model)假设附着涡在平行于壁面的速度分量(u和w)上的流向尺度正比于附着涡所处的法向高度^[16]，即：l_x~y，这一假设暗示了附着涡在不同的流层上具有相同的结构倾角。图 4(a, c)显示，原始流场u、w分量的结构倾角确实在一定的y⁺以上不再改变，但大/小两种尺度的u、w分量的结构倾角并未呈现出类似的收敛趋势。因此，不能简单认为在大尺度模态中发现的从壁面延伸向外区的条状结构就是附着涡的物理表征。

4 结论

本工作应用MC-QBVMD对湍流边界层x-y平面脉动速度场进行了分解。MC-QBVMD的大尺度模态对应的是与LSM相似的结构。小尺度模态则比较复杂，在近壁区倾斜的Q2和Q4事件比较相似，而在对数区则近似于垂直的涡结构。MC-QBVMD小尺度模态的在流向尺度和形态特征上与流向涡的尺度非常接近，且两者在不同雷诺数的工况下均有较好的对应关系，表明x-y平面的MC-QBVMD小尺度模态应该对应于原始流场中的小尺度流向涡结构在x-y平面的投影。

在当前工作中我们只将原始流场分解为大/小两种尺度的模态，这样做虽然便于表征和分析，但可能并不能够最好的反映湍流的多尺度特性。更多关于模态数量以及更高维度的分解的工作需要在后续的工作中开展。