引言

涡激振动是一种典型的流致振动现象，流体绕结构产生周期性的旋涡脱落，流体的涡脱频率接近结构的自振频率时，质量轻柔、阻尼比低的结构极易发生的一种限幅振动.如大跨度桥梁、高层建筑以及海洋立管等结构就经常在极低的风速下发生涡振，影响结构的疲劳寿命^[1].

矩形断面作为桥梁结构以及高层建筑结构的典型代表性断面，其流场复杂性和气动力的不稳定性，引起了广泛的关注.与经典的圆形或方形圆柱体直接由上下剪切层之间的相互作用触发，之后形成交替脱落的旋涡不同，矩形断面是由于旋涡从前缘分离的剪切层向上卷起，并周期性地重新附着在圆柱体侧面上并滑移，之后与尾流干涉失稳，最终在剪切层尾端形成交替脱落的Karman涡，即所谓的“撞击剪切层不稳定性”^[2-4].根据Matsumot等^[5]的研究，竖向涡振的运动诱导涡旋激励是剪切层不稳定的结果，它在横截面前缘产生涡旋，涡旋在体表面向后缘移动，速度约等于迎面流速度的60%.在一个振荡周期中，前缘涡到达后缘并与该处产生的二次涡结合.此外，运动诱导涡似乎触发扭转颤振，在Karman涡街存在时，扭转颤振的临界风速会提高.在2008 年的第六届国际钝体空气动力学和应用学术讨论会（BBAA6）上^[6]，学者们正式定下了5∶1 矩形断面的空气动力学基准性研究（BARC），之后有大量风洞试验和数值模拟对其开展研究相关研究.Zhang等^[7]通过3D LES研究了展现长度和网格分辨率对矩形柱流场的影响，Nguyen等^[8]通过风洞试验结合CFD计算分析了5∶1矩形圆柱周围的表面压力分布和相关性，揭示了涡激振动是由运动诱导的前缘涡引发的，接近后缘的强相关流动特征是结构响应增加的原因.但是现有研究大都直接对结构某一阶物理量进行分析，如速度、压力、涡量等，而流场的细节往往非常复杂，包含了大量的细节信息.如果直接对这些高维度的流场数据进行处理，会使得流场的分析和控制变得非常困难.

由于在非常复杂的流动中也存在特定模式，为了进一步了解矩形复杂的流动结构以及其流态演变机理，我们通过流场降阶模型来从高维的非线性流场数据中提取低维流动特征.近些年来，随着高精度数值模拟技术的迅猛发展，低维降阶技术在复杂流动分析中的应用日益广泛^[9，10]，对于流场降阶最常用的为本征正交分解（POD）与动力学模态分解（DMD）两种方法，叶坤等^[11]发现相比于POD，DMD更加适用于圆柱绕流的稳定性分析.寇家庆等^[12，13]利用DMD对跨声速抖振进行了模态分析，并对其在流体力学中的应用进行了综述.

已有研究表明，二维RANS模型可能会高估分离泡的大小和侧向面上的压力波动^[14，15]，但从Bruno等^[16]通过POD对 LES模拟得到的压力场分解结果来看，虽然流场的三维特征不可忽视，但气动力主要还是由二维模态决定的，其频率基本不会有大的改变.

因此，考虑到实际的计算资源有限，本研究采用二维雷诺平均（RANS）SST k-ω湍流模型对5:1矩形断面的涡激振动进行了二维数值模拟计算.为进一步了解矩形断面流固耦合过程中的旋涡脱落机理，低维模态空间的角度探究了静/涡振态矩形断面的绕流模态的特征，并从流场模态演化的角度探究涡激振动内在的演变机制，通过动力学模态分解（dynamic mode decomposition，DMD）方法对其静态绕流场和涡振态流场进行了模态分解，基于能量占比的大小，分析涡振最大无量纲振幅处各阶模态表征的物理意义，并基于DMD的主导模态建立降阶模型，重构了矩形断面的涡量场.

1 数值方法

1.1 基本控制方程

为了获得结构周围流体的流动状态，需要对二维矩形断面周围流体进行数值模拟.流体状态可视为二维不可压缩黏性流体，忽略能量方程，雷诺平均后，得到连续性方程与动量方程可以表达为:

其中，t为时间，p为平均压力，ρ和v分别为空气的密度和运动黏性系数，下标（i，j=1-2）为坐标轴的两个方向，${\stackrel{-}{u}}_i$和${\stackrel{-}{u}}_{}^{\mathrm{\text{'}}}$表示沿着i轴的平均速度和脉动速度，$\partial {\stackrel{-}{u}}_i{\stackrel{-}{u}}_j$为雷诺应力项.本研究控制方程的离散和求解基于商业软件Fluent实现，其中湍流模型采用对分离流动适用性良好的SST k-ω模型，选择有限体积法求解控制方程，空间离散采用二阶迎风格式，时间离散采用二阶隐式时间积分，速度-压力耦合方程则采用SIMPLEC算法来求解，计算收敛残差精度设置为1E-6.

1.2 计算模型

基于文献^[17]的风洞试验参数，对宽高比为5:1的矩形断面涡激振动进行二维数值模拟，其中截面宽度B=0.3m，高度D=0.06m，采用自编译UDF程序嵌入Fluent软件，实现矩形竖向单自由度的流固耦合计算.矩形断面的具体结构参数如表1所示.

矩形断面可简化为竖向自由度弹簧-质量-阻尼系统，矩形断面振动模型如图1所示.采用Newmark-β法求解结构的动力学方程:

式中，y（t）、$\dot{y}（t）$、$\ddot{y}（t）$分别为结构瞬时的竖向位移、速度、加速度; m、c和k分别为结构单位长度的质量、阻尼及刚度; F_L（t）为结构在竖向方向的瞬时升力，可以通过Fluent内置的UDF宏命令提取得到.

1.3 网格及边界条件设置

模型的计算域及边界条件设置如图2所示，上、下侧边界距离模型形心均为10B，上游入口边界及下游出口边界距离模型形心分别为10B、20B，整体阻塞率为1%.来流风攻角为0°攻角，上游边界采用速度入口（Velocity-inlet），下游边界采用压力出口（Pressure-outlet），上下侧边界均采用对称边界（Symmetry）.经过网格无关性验证及时间步长无关性验证后，采用网格数量为60400的二维数值网格和时间步长t=0.001s进行矩形断面涡激振动数值模拟.

整体网格示意图如图3所示，网格计算域划分为刚性域+动网格域，刚性域中结构附近边界层网格为20层，第一层网格高度为1E-4 m，动网格域采用四边形网格，整体网格增长率不超过1.05，在最高计算风速（U₀=4.43 m/s）下的网格整体YPlus均值小于1，符合湍流模型的基本要求^[17].动网格技术（Dynamic mesh）采用弹簧光顺（Smoothing-spring），为保证动网格更新迭代趋于收敛，迭代步数设置为30步.

1.4 DMD方法

动力学模态分解是对原始数值模拟的结果进行采样，得到N个时刻的流场快照信息，即$X_1^N=\left\{x_1，x_2，x_3，\cdots ，x_{N-1}\right\}$，记录以时间为间隔的流场数据，并建立以下两个快照矩阵.

假设上个时刻流场x_i和下个时刻x_i+1的变化为线性关系，即

根据上式（4）的假设可得X₂=AX₁，其中为A系数矩阵.只要确定了A矩阵或者获得了它的特征值和特征向量，就可获悉该系统的主要特征和演变规律，而DMD的核心就是对上述快照进行处理，提取出系统主导的特征值和主要模态.

对矩形X₁进行奇异值分解（SVD）:

其中矩阵S为由包含r个奇异值组成的对角矩阵，U，V为酉矩阵，满足${\mathit{U}}^{*}\mathit{U}=\mathit{I}，{\mathit{V}}^{*}\mathit{V}=\mathit{I}，{\mathit{U}}^{*}，{\mathit{V}}^{*}$分别为U，V的复共轭转置矩阵.

将式（7）带入X₂=AX₁可得:

由于此时获得的A矩阵维度很高，其求解计算量很大.因此，可以提供一个矩阵$\tilde{\mathit{A}}$来代替高维矩阵A，且这两个矩阵相似，将A投影至X₁的正交基U上进行降阶，得到A的相似矩阵:

用的$\tilde{\mathit{A}}$特征值来近似估计主要的A的特征值，并计算出相似矩阵$\tilde{\mathit{A}}$的特征值和特征向量.记A的第j个特征值为μ_j，特征向量为A_j.第j阶DMD模态为:

各阶DMD模态的流动增长率（衰减率）g_i和振荡频率f_j为:

任意时刻流场可表示为:

其中α_j表示第j阶模态

2 模型验证

矩形断面风洞试验与二维数值模拟的涡激振动结果对比如图5所示，矩形断面的涡振折减风速（V=U/Df_h）区间为V=6.2~11.5，折算风速V=11时达到了最大无量纲位移，h/D=0.064.对比风洞试验和数值模拟结果，其涡振发生区间基本吻合，其涡振位移响应也在可接受的范围内.

图6展现了各工况下升力主频随折减风速变化的曲线，其涡振折减风速区间为V=8~11，流体涡脱频率发生明显的锁定现象，此时升力主频与结构固有频率之比接近1，当在锁定区外，涡振时的升力主频与静止状态下的涡脱频率相同，验证了此数值模拟方法的准确性.

图7给出了无量纲振幅最大（V=11）时的位移时程图，把这个阶段涡振响应分为初始段、发展段和稳定段三个部分，把无量纲时间（0~1500）划分为初始段，此时处于涡振最开始，流固耦合作用剧烈，流场极其不稳定; 把（1500~3000）划分为发展段，位移响应发展迅速; 3000之后划分为稳定段，此时涡振状态呈极限环振动.并后续对这三个部分分别进行DMD分析.

此时，涡振态矩形稳定振动时的频率为6.13Hz，对应的无量纲频率St=0.091，而该风速下静态矩形的涡脱频率为7.86Hz，其对应的无量纲频率St=0.116.

3 计算结果及分析

3.1 静态矩形的DMD分解

为了探究矩形断面涡振发生时的旋涡脱落机理，本文采用DMD算法对V=11（此时无量纲振幅最大）处以及对应风速下的静态矩形涡量场进行分析.首先提取静态矩形流场稳定后的10个周期的非定常涡量场数据，每个周期采样40次，定义模态的二范数为该阶模态的能量^[8，19]，之后根据模态能量的大小进行排序.

对提取的流场快照进行DMD分解，图8为静态矩形DMD模态Ritz特征值的分布情况，其中颜色的深浅对应各阶模态能量的大小，颜色越深表示该阶模态所对应的能量越大，模态的特征值基本分布在单位圆上，说明此时对应的静态矩形绕流场的频率趋于稳定.同时，表2给出了各阶模态对应的增长率与折减频率的关系.流动状态稳定后各阶模态的增长率都很小，说明此时前5阶模态均属于稳定模态.

DMD模态能量与频率的关系如图9所示，其中能量占比最大的为第1阶模态，其频率为0，增长率也接近为0，表征着流场的平均模态.其余流场的振荡模态均成对出现，对应特征值的共轭复数，因此只分析其中的一阶模态.

除开平均模态，流场的能量主要分布在2~4阶模态上，第2阶模态的频率与静态矩形涡脱主频相同，是稳定流动阶段流场的主导模态，模态3为静态矩形旋涡脱落的二倍频，是尾缘涡交替脱落的模态，模态4为静态矩形旋涡脱落的三倍频，是流场的高阶谐波模态，是对流场的补充.

流场振荡的频率除开涡脱主频及其倍频成分，其余能量占比都很小.选取其前4阶DMD模态云图进行分析，如图10所示，第1阶模态为平均模态，通过数值计算可知其与时均化流场基本相同，第2阶模态可以看出矩形整体的旋涡脱落结构为正对称的，其3阶模态的整体旋涡脱落结构为反对称的，4阶模态整体旋涡脱落结构为正对称的.对于静态矩形而言，前缘涡分离后，受到撞击剪切层不稳定性的作用不断发展，但发展到下游后，断面周围的旋涡尺度偏小，在与尾涡的合并中占弱势地位，其尾缘脱落的Karman涡强劲，主要的旋涡脱落频率由尾涡来控制，从整体上来看，随着模态阶数的增加，断面周围涡的尺度都在相应地降低，各阶模态形成的尾流都有很强的相关性.

3.2 涡振态矩形的DMD分解

涡振态矩形最大无量纲振幅为V=11处，把该风速下整个涡振态分为初始段、发展段以及稳定段三个阶段，每个阶段都提取12个周期的非定常涡量场数据，每个周期采样40次，对提取的流场快照进行DMD分解，以此来分析矩形涡振在振幅发生、发展到稳定的过程中，流场模态的变化趋势.

图11为涡振态矩形DMD模态Rize特征值的分布，初始阶段有模态分布在单位圆外，表明该模态不稳定，其余模态分布在单位圆上或单位圆内，表明该阶模态为周期性模态或稳定模态.对于涡振各个阶段，第一阶对应的平均流场模态都为周期性的稳定模态.

在涡振发展的不同阶段，主导流场的模态频率会有所不同，表3给出了各阶段前5阶DMD模态及其对应的增长率与频率，对于初始段，此时的流场还没有完全稳定，主导此时流场振动的模态为2阶模态，其对应St频率为0.117，与静态流场的涡脱频率接近，其增长率为-0.028，即该模态为不断衰减的稳定模态，与之相对应的2倍频模态3同样也为不断衰减的稳定模态.其3阶模态对应St频率为0.091，和结构固有频率接近，并且增长率为0.111，为不稳定模态，表明该流体模态正处于和结构模态耦合的过程中，还在不断发展.4阶模态频率相对较小，这可能是由剪切层的不稳定性造成的自持性振荡^[18]，增长率为-0.048说明该模态为强烈衰减的稳定模态.到了发展阶段，模态2~4的频率均表现为结构主频及其倍频，对应的模态也均为稳定模态，其中3阶模态与4阶模态频率都与结构3倍频接近，但3阶模态衰减率很大.对于稳定阶段，主要的流场模态频率已被结构频率及其倍频锁定，各阶模态的增长率都接近于0，对应的模态均为周期性的稳定模态.

图12给出了涡振态矩形DMD模态能量与频率的关系，从能量的角度来看，初始段能量分布于各阶模态，初始能量都很低，能量的分布较为平均，此时的各阶模态还属于竞争状态，到了发展阶段，流场模态的能量逐渐被结构主频模态及其倍频模态捕获，此时高能量的模态都集中于结构主频模态及其倍频模态周围，但此时并不是单一模态，而是每一阶结构频率附近同时存在几阶流场模态互相竞争，对于稳定阶段，流场模态竞争已经完成，流场模态已完全被结构频率模态捕获，此时流场达到了新的周期性平衡.

图13为涡振稳定段前4阶DMD模态云图，此时的流场涡脱频率被结构频率锁定，相较于静态流场模态，涡振态流场主要旋涡脱落模式与静态流场相似，但是对于涡振态，由于由剪切层不稳定引起的冲击前缘涡合与运动诱导涡合并向下运输，其能量对于尾缘涡来说更加占据主要地位，运动诱导涡与结构振动频率同步，从而形成了尾涡的锁定，而由于截面的振荡，其各阶模态的尾流相关性较静态流场差，尾涡分布也更为分散.

3.3 非定常流场分析

图14为静态矩形2阶模态一个周期内流场旋涡演化，对于静态矩形而言，其涡脱主要由上游造成，形成了交替脱落的前缘涡，由于剪切层的K-H不稳定性，旋涡在剪切层周围不断发展，形成大尺度旋涡，之后再附都矩形为的尾部，与尾缘涡形成干涉，从而从尾缘产生交替脱落的Karman涡.

图15为涡振态矩形2阶模态一个周期内流场旋涡演化，相较于静态矩形，涡振态矩形前缘处脱落的旋涡明显更加强烈，其旋涡的尺度更大，这是由于除冲击前缘形成的涡脱外，涡振态矩形还会产生运动诱导的旋涡^[5，20]，由于本文数值模拟采用二维的雷诺平均SST k-ω湍流模型，其小尺度的旋涡会直接数值平均为大尺度旋涡，从而未能直接从旋涡的形态上直接捕捉到运动诱导涡与前缘分离涡的区别，但是准确捕捉到了其旋涡强度的变化.相比与静态，其涡振态旋涡再附点也会提前，这进一步缩小了前缘涡与尾缘涡的相位差，从而导致涡振态矩形的尾缘涡脱频率会低于静态矩形，形成了独特的涡振的锁频现象.

4 结论

本文基于数值模拟方法对长宽比为5∶1的矩形断面在竖向自由度的涡激振动进行了计算，并通过DMD方法，从低维模态空间的角度探究了静/涡振态矩形断面的绕流模态的特征，最后通过稳定态下主要流体振荡模态的非定常流场，对矩形的旋涡脱落模式进行了分析，结论如下:

1）基于DMD方法准确地捕捉到了流体模态的特征，各阶模态流场存在显著差异.流体状态稳定后，流场主要由频率为0的模态，与涡脱主频相等模态以及其对应的高阶谐波模态组成.

2）静/涡振态矩形的流场模态特征存在显著区别.相较于静态流场模态，涡振态流场主要旋涡脱落模式与静态流场相似，但其模态频率更低，且结构剪切层周围的旋涡强度更高，与尾涡而合并后在尾缘脱落的Karman涡相关性也越差.

3）从流场模态演化的角度阐明了涡振从产生、发展到稳定极限环运动过程中内在的演变机制.在涡振发展的不同阶段，主导流场的模态频率会有所不同.初始段，主要模态的能量集中于静态流场的涡脱频率上; 发展段，流体频率会向结构的固有频率及其倍频靠近，此时一个结构频率周围可能存在多个与之接近的流体模态.稳定段，流体模态频率完全被结构所捕获，流场趋于新的平衡.

4）从低维流场模态角度解释了涡振的锁频现象的原因.在振动状态下，未能直接捕捉到运动诱导涡，由其周围旋涡的强度推测，其与剪切层不稳定引起的冲击前缘涡合并后，在结构尾缘时又与尾涡合并，导致前缘涡与尾缘涡的相位差减小，从而使得涡振态矩形的尾缘涡脱频率会低于静态矩形，形成了独特的涡振的锁频现象.

参考文献