基于增量谐波平衡法的涡激振动抑制系统的周期解及其动力学分析

Periodic Solutions and Dynamic Analysis of a Vortex⁃Induced Vibration Suppression System Based on the Incremental Harmonic Balance Method

摘要

本文提出了应用于流体环境的圆柱体结构压电能量收集吸振器， 研究在涡激振动下实现结构振动抑制和压电能量收集一体化策略。首先采用数值方法研究了减振器对涡激振动系统的被动振动抑制效果， 并采用增量谐波平衡法获得了涡激振动能量采集器的半解析周期解， 依据Floquet理论确定了周期解的稳定性。研究发现， 在涡激振动系统的锁定区内， 非线性流固耦合效应可诱导多解共存非线性动力学行为。当系统配置线性吸振装置时， 系统的多解共存现象消失， 振动主频不再固定于结构固有频率， 其值随风速呈现逐渐递增的趋势。另外， 分析了线性减振器的尖端钝体质量， 以及内梁的长度、 宽度和厚度等参数对主结构幅值变化吸振器电压输出的影响规律， 系统对比评价了相关参数对减振效果的作用， 明确了线性减振器高效设计的参数可选区间。

Abstract

This study proposed a novel piezoelectric energy harvesting dynamic vibration absorber （DVA） integrated into cylindrical structures subjected to fluid flow， and aimed to achieve simultaneous vortex-induced vibration suppression and piezoelectric energy harvesting. Firstly， numerical simulations were employed to investigate the passive vibration suppression effectiveness of the dynamic vibration absorber on the vortex-induced vibration system. Subsequently， the incremental harmonic balance （IHB） method was adopted to derive semi-analytical periodic solutions for the vibration energy harvester， with the solution stability analyzed via Floquet theory. The study reveals that within the lock-in range of the vortex-induced vibration system， nonlinear fluid-structure interaction （FSI） effects can induce coexisting multiple solutions and complex nonlinear dynamic behaviors. Notably， when the system is equipped with a linear vibration absorber， the system exhibits eliminated coexisting solutions and a dominant vibration frequency that monotonically increases with wind speed， rather than remaining locked to the structural natural frequency. Furthermore， a parametric study evaluates key design variables such as the tip mass of the bluff body and beam dimensions on vibration mitigation performance. This work quantifies both the primary structure’s amplitude attenuation and the DVA’s voltage output， and establishes optimal parameter ranges for high-efficiency linear DVA design.

Graphical abstract

关键词

涡激振动 / 能量收集 / 增量谐波平衡法 / 周期解 / 频率锁定

Key words

vortex-induced vibration / energy harvesting / incremental harmonic balance method / periodic solutions / frequency lock-in

0 引　言

涡激振动（Vortex Induced Vibration， VIV）是一种具有自限性振幅特性的自激振动现象， 广泛存在于航空航天、 海洋工程、 土木工程等领域^［1‐5］。当柔性或弹性支撑的钝体受来流作用时， 其尾部交替脱落的漩涡可诱发大幅振动响应^［6‐7］。在一定的流速范围内， 涡流脱落频率与结构固有频率同步， 结构会在锁定或同步状态下表现出类似共振的振荡运动^［8］。但是， 这种特殊的涡激振动状态容易引起噪声和疲劳裂纹， 从而影响结构安全^［9］。因此， 对于这类结构而言， 要实现高效设计， 精准预测其在频率锁定状态下的动态响应就显得尤为关键， 同时还需及时采取有效措施控制振动响应。当前对于这类振动的抑制方法可以分为主动抑制和被动抑制两种。Zhang等^［10］通过数值模拟研究了VIV的电磁控制（洛伦兹力控制）。Akhtar等^［11］在一个固定圆柱上使用了一对吸力制动器， 并观察到波动力的抑制和平均阻力系数的减小。此外， Dai等^［12］利用非线性时滞反馈来抑制弹性安装的圆柱体的振幅， 并认为它在控制空气动力系统的响应方面是有效的。Mehmood等^［13］研究了线性和非线性速度反馈控制器在降低锁定状态下弹性安装刚性圆柱体的VIV方面的有效性。这些控制方法主要是使用主动方法进行抑制， 相较于主动抑制方法， 被动式压电能量收集减振器在成本、 可靠性和适用场景方面的优势更明显。被动式压电减振器无需复杂传感器、 控制器及持续能源供应， 可以大幅降低设备购置、 安装与维护成本； 同时因其无复杂的电子元件和活动部件， 故障点减少， 受环境因素影响较小， 因而使用寿命较长； 另外， 在偏远无稳定供电地区、 强电磁干扰环境， 或对控制实时性要求不高的一般工程场景， 被动式压电能量收集减振器无需复杂控制算法和能源输入， 能够稳定发挥减振作用。

已有的被动控制方案中， 大多数振动能量被阻尼所耗散， 如果能将冗余的振动能量收集起来， 就可以在一些恶劣条件下为传感器元件持续稳定地供电。同时， 相对于自然界中其他能量， 振动能量具有较高的能量密度， 更易捕获， 具有更高的能量转换效率。因此， 研究人员开始致力于如何捕获环境中的振动能并将其转化为电能， 以实现低功耗电子设备的自给供能^［14］。压电式振动俘能器装置利用压电效应^［15］将机械能转化为电能， 其环境适应性强， 并且能与微电子技术较好地融合。然而， 目前基于涡激振动的能量收集侧重于直接收集振动能量， 要求主结构具有较大的振动幅值， 但这会危害结构的安全性。

涡激振动系统是一种高维强非线性振动系统， 在其求解过程中， 数值方法虽然简单， 精度高， 但难以揭示多解共存、 突弹跳跃等非线性现象， 需要借助解析或半解析方法^{［16‐17］}。传统的解析方法（正规摄动法、 Lindstedt-Poincaré法、 平均法、 多尺度法和KBM法等）难以求解其高阶周期解^{［18‐19］}。Zhou等^［20］使用谐波平衡分析发现， 在高激励水平下， 非线性涡激振动系统存在多解共存现象。然而， 目前通过解析或半解析方法研究涡激振动抑制及能量收集系统的周期解的工作相对较少。

增量谐波平衡法（Incremental Harmonic Balance Method， IHBM）由Lau等^［21］在1981年首次提出， 其将非线性弹性系统周期运动的幅度增量变分原理与Rayleigh-Ritz法相结合， 形成了增量谐波平衡法。随后， Lau等^［22］提出了将增量谐波平衡法和时间、 空间域的有限条法相结合的增量时间-空间有限元法， 并进一步将增量谐波平衡法与有限法相结合， 用于分析弹性薄板的非线性振动问题。增量谐波平衡方法结合了增量迭代法与谐波平衡法的特点， 其灵活的收敛精度控制能力使其成为目前求解强非线性问题应用较为广泛的方法之一。

综上所述， 目前涡激振动和能量收集的研究较多集中于单一功能， 对于同时实现振动抑制以及能量收集的研究较少。此外， 对于安装减振装置后的涡激振动抑制系统， 尚未开展解析解、 多解和跳跃等非线性动力学行为研究， 其对于完善振动抑制以及能量收集效果至关重要。针对以上不足， 本文的创新性贡献如下： 提出压电能量收集吸振器用于涡激振动系统， 力求实现被动抑制振动和收集能量的双重功能； 应用IHB方法求解了涡激振动系统的周期解， 介绍了具有线性吸振器的涡激振动系统的模型设计、 工作原理及IHB解的推导； 描述了吸振器参数对主结构减振效果的影响， 揭示了压电能量收集吸振器抑制涡激振动的内在机制。

1 模型及系统方程

本文研究的系统由安装在悬臂梁自由端的圆柱体及其内部集成的压电能量收集吸振器组成。该吸振器为等宽矩形截面的压电悬臂梁， 其自由端附加质量块， 通过正压电效应将机械能转化为电能。此设计既能耗散主结构振动能量， 又可实现能量回收。图 1 所示为压电能量收集吸振器示意图， 该结构受到均匀流速的影响。图 1（a） 中w_{\mathrm{1}}(x,t)表示主结构的绝对位移， 而w_{\mathrm{2}}(x,t)表示吸振器与主结构的相对位移。图 1（b） 所示为压电悬臂梁的截面图。L_{bn}和L_b分别表示压电层和外梁长度， b_n表示压电层的宽度， h_s和h_p分别表示梁和压电层的厚度。在低风速激励环境下， 悬臂梁结构的振动能量主要集中于基频附近， 其一阶模态因固有频率与激励频带高度匹配而占据主导地位， 而高阶模态对应的振动形态贡献度较低。尤其对于本文所研究的低风速环境下的悬臂式振动结构而言， 高阶模态引发的振动响应幅度极小， 可忽略其干扰影响。因此， 在后续关于外部支撑弹性梁与减振器中压电梁的动力学分析中， 均仅考虑一阶模态的作用， 以简化理论模型并突出主要振动特征。

早期研究表明， Van der Pol型和 Rayleigh型自激振子可以用来模拟波动尾流行为^{［23‐25］}。后续研究又对经典尾流振子模型提出了几种修正^{［26‐27］}。Facchinetti等^［28］对基于位移、 速度和加速度耦合的低阶现象模型进行了分析， 发现考虑加速度耦合的模型能够定性和定量地描述涡激振荡现象的持续锁定区域。因此， 本文将使用上述加速度耦合的尾流振荡器模型进行分析。需要指出， 这些降阶模型仅适用于低雷诺数（Re）区域（300Re1.5×10⁵）。 本文研究的涡激力应为

\begin{array}{c}{F}_{\mathrm{V}\mathrm{I}\mathrm{V}}=\frac{C_{L\mathrm{0}}{\rho }_{\mathrm{0}}D{U_{\mathrm{0}}}^{\mathrm{2}}{L}_{\mathrm{0}}}{\mathrm{4}}q-\frac{C_d{\rho }_{\mathrm{0}}D{U}_{\mathrm{0}}{L}_{\mathrm{0}}}{\mathrm{2}}\cdot \\ \left[{\varphi }_{\mathrm{1}}(L_b)+\frac{D}{\mathrm{2}}{{\varphi }_{\mathrm{1}}}^{\mathrm{\text{'}}}(L_b)\right]\dot{\eta },\end{array}（1）

式中： C_d为阻力系数， 考虑为1.2； C_{L\mathrm{0}}为名义升力系数， 确定为0.3； U_{\mathrm{0}}为气流速度； {\rho }_{\mathrm{0}}为空气的密度； L_{\mathrm{0}}为圆柱体的长度； {\varphi }_{\mathrm{1}}(x)为主结构的一阶模态振型； \eta (t)为第一个模态的模态坐标。

研究发现， 考虑加速度耦合的模型定性和定量地捕获了涡激振动现象的持续锁定区域， 本文将使用上述加速度耦合的尾流振荡器模型进行分析。尾流振子的方程形式可以写为

\begin{array}{c}\ddot{q}+\mathrm{\Lambda }{w}_s(q^{\mathrm{2}}-\mathrm{1})\dot{q}+{w_s}^{\mathrm{2}}q=\\ \frac{P}{D}\left[({\varphi }_{\mathrm{1}}(L_b)+\frac{D}{\mathrm{2}}{{\varphi }_{\mathrm{1}}}^{\mathrm{\text{'}}}(L_b)\right]\ddot{\eta },\end{array}（2）

式中： P和\mathrm{\Lambda }是实验确定的常数， 分别为12和0.24； w_s为涡流脱落频率， 定义为w_s=\frac{\mathrm{2}\mathrm{\pi }{S}_t{U}_{\mathrm{0}}}{D}。斯特劳哈尔数S_t是雷诺数和圆柱体横截面的函数， 在此分析中数值为0.2。因此， 系统的运动方程可以表示为

\begin{array}{c}\ddot{\eta }+\mathrm{2}{\zeta }_{\mathrm{11}}{w}_{\mathrm{11}}\dot{\eta }+{w_{\mathrm{11}}}^{\mathrm{2}}\eta +{\lambda }_{Lb}{\lambda }_{Lbn}\ddot{r}+m{{\lambda }_{Lb}}^{\mathrm{2}}\ddot{\eta }={\lambda }_{Lb}{F}_{viv},\\ \ddot{r}+\mathrm{2}{\zeta }_{\mathrm{21}}{w}_{\mathrm{21}}\dot{r}+{w_{\mathrm{21}}}^{\mathrm{2}}r+{\lambda }_{Lbn}{\lambda }_{Lb}\ddot{\eta }-\theta V=\mathrm{0},\\ \ddot{q}+\mathrm{\Lambda }{w}_s(q^{\mathrm{2}}-\mathrm{1})\dot{q}+{w_s}^{\mathrm{2}}q=\frac{P}{D}{\lambda }_{Lb}\ddot{\eta },\\ C_P\dot{V}+\frac{V}{R}+\theta \dot{r}=\mathrm{0},\end{array}（3）

式中： {\lambda }_{Lb}={\varphi }_{\mathrm{1}}(L_b)+\frac{D}{\mathrm{2}}{{\varphi }_{\mathrm{1}}}^{\mathrm{\text{'}}}(L_b)； {\lambda }_{Lbn}=M_t{\psi }_{\mathrm{1}}(L_{bn})+\frac{B}{\mathrm{2}}{M}_t{{\psi }_{\mathrm{1}}}^{\mathrm{\text{'}}}(L_{bn})+m_{\mathrm{2}}{\int }_{\mathrm{0}}^{L_{bn}}{\psi }_{\mathrm{1}}\mathrm{d}x； m=M_t+m_{\mathrm{2}}{L}_{bn}； C_p=\frac{\mathrm{2}{e}_{\mathrm{33}}{b}_n{L}_{bn}}{h_p}为压电能量收集吸振器的电容； \theta =E_P{d}_{\mathrm{31}}{b}_n(h_p+h_s){{\psi }_{\mathrm{1}}}^{\mathrm{\text{'}}}(L_{bn})为压电耦合项， {\psi }_{\mathrm{1}}(x)为吸振器的一阶模态振型； r(t)为第一个模态的模态坐标； {\zeta }_{\mathrm{11}}和w_{\mathrm{11}}分别为VIV系统外部支撑梁的机械阻尼比和一阶外梁固有频率； {\zeta }_{\mathrm{21}}和w_{\mathrm{21}}分别为内部压电梁的机械阻尼比和一阶内梁固有频率； 考虑外梁的第一阶模态频率w_{\mathrm{11}}=24.4441 rad/s。

2 周期解及其稳定性分析

2.1 增量谐波平衡法

为了揭示VIV系统的非线性动力学行为， 采用IHB方法推导了其高阶近似周期解。首先， 将等式（3）重新写成状态方程

\left[\begin{array}{cccc}\mathrm{1}+m{{\lambda }_{Lb}}^{\mathrm{2}}& {\lambda }_{Lb}{\lambda }_{Lbn}& \mathrm{0}& \mathrm{0}\\ {\lambda }_{Lb}{\lambda }_{Lbn}& \mathrm{1}& \mathrm{0}& \mathrm{0}\\ -{\beta }_m{\lambda }_{Lb}& \mathrm{0}& \mathrm{1}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}\right]\left[\begin{array}{c}\ddot{\eta }\\ \ddot{r}\\ \ddot{q}\\ \ddot{v}\end{array}\right]+\left[{\begin{array}{cccc}\mathrm{2}{\zeta }_{\mathrm{11}}{w}_{\mathrm{11}}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{2}{\zeta }_{\mathrm{21}}{w}_{\mathrm{21}}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \theta & \mathrm{0}& C\end{array}}_P\right]\left[\begin{array}{c}\dot{\eta }\\ \dot{r}\\ \dot{q}\\ \dot{v}\end{array}\right]+

\left[\begin{array}{cccc}{{\lambda }_{Lb}}^{\mathrm{2}}{C}_m& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& -{\beta }_m& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}\right]\left[\begin{array}{c}\dot{\eta }\\ \dot{r}\\ \dot{q}\\ \dot{v}\end{array}\right]U+\left[\begin{array}{cccc}\mathrm{0}& \mathrm{0}& -C_L{\lambda }_{Lb}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& {{\beta }_s}^{\mathrm{2}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}\right]\left[\begin{array}{c}\eta \\ r\\ q\\ v\end{array}\right]U^{\mathrm{2}}+

\left[\begin{array}{cccc}{w_{\mathrm{11}}}^{\mathrm{2}}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& {w_{\mathrm{21}}}^{\mathrm{2}}& \mathrm{0}& -\theta \\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \frac{\mathrm{1}}{R}\end{array}\right]\left[\begin{array}{c}\eta \\ r\\ q\\ v\end{array}\right]+\left[\begin{array}{c}\mathrm{0}\\ \mathrm{0}\\ {\beta }_w\\ \mathrm{0}\end{array}\right]U{q}^{\mathrm{2}}\dot{q}=\mathrm{0},（4）

式中： {\beta }_w=\mathrm{2}\mathrm{\pi }{S}_t\mathrm{\Lambda }{D}^{-\mathrm{1}}， {\beta }_s=\mathrm{2}\mathrm{\pi }{S}_t{D}^{-\mathrm{1}}， {\beta }_m=\frac{P}{D}， C_m=\frac{\mathrm{1}}{\mathrm{2}}{C}_d{\rho }_{\mathrm{0}}{L}_{\mathrm{0}}D。

引入时间尺度变化\tau =\Omega t（\Omega是响应频率）以及状态变量

\mathit{Y}={[\begin{array}{cccc}{y}_{\mathrm{1}}& y_{\mathrm{2}}& y_{\mathrm{3}}& y_{\mathrm{4}}\end{array}]}^{\mathrm{T}}={[\begin{array}{cccc}\eta & r& q& v\end{array}]}^{\mathrm{T}}，

等式（4）写为

\begin{array}{c}{\Omega }^{\mathrm{2}}\mathit{M}{\mathit{Y}}^{″}+\Omega \mathit{C}{\mathit{Y}}^{\text{'}}+\Omega {\mathit{C}}_{u}U{\mathit{Y}}^{\text{'}}+\mathit{K}\mathit{Y}+\\ {\mathit{K}}_u{U}^{\mathrm{2}}\mathit{Y}+\Omega U{\mathit{N}}_{\mathrm{1}}(\mathit{Y},{\mathit{Y}}^{\text{'}})=\mathrm{0},\end{array}（5）

式中： {()}^{\text{'}}表示\frac{\partial ()}{\partial \tau }。方程（5）中各个系数矩阵和非线性参数向量表示为

\mathit{M}=\left[\begin{array}{cccc}\mathrm{1}+m{{\lambda }_{Lb}}^{\mathrm{2}}& {\lambda }_{Lb}{\lambda }_{Lbn}& \mathrm{0}& \mathrm{0}\\ {\lambda }_{Lb}{\lambda }_{Lbn}& \mathrm{1}& \mathrm{0}& \mathrm{0}\\ -{\beta }_m{\lambda }_{Lb}& \mathrm{0}& \mathrm{1}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}\right]，

\mathit{C}=\left[{\begin{array}{cccc}\mathrm{2}{\zeta }_{\mathrm{11}}{w}_{\mathrm{11}}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{2}{\zeta }_{\mathrm{21}}{w}_{\mathrm{21}}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \theta & \mathrm{0}& C\end{array}}_P\right]，

{\mathit{C}}_u=\left[\begin{array}{cccc}{{\lambda }_{Lb}}^{\mathrm{2}}{C}_m& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& -{\beta }_m& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}\right]，

\mathit{K}=\left[\begin{array}{cccc}{w_{\mathrm{11}}}^{\mathrm{2}}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& {w_{\mathrm{21}}}^{\mathrm{2}}& \mathrm{0}& -\theta \\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \frac{\mathrm{1}}{R}\end{array}\right]，

{\mathit{K}}_u=\left[\begin{array}{cccc}\mathrm{0}& \mathrm{0}& -C_L{\lambda }_{Lb}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& {{\beta }_s}^{\mathrm{2}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}\right]，

{\mathit{N}}_{\mathrm{1}}(\mathit{Y},{\mathit{Y}}^{\text{'}})=\left[\begin{array}{c}\mathrm{0}\\ \mathrm{0}\\ {\beta }_w\\ \mathrm{0}\end{array}\right]y_{\mathrm{3}}^{\mathrm{2}}{\dot{y}}_{\mathrm{3}}\mathrm{。}

IHB方法的第一步是增量过程， 即使用Newton-Raphson 程序将微分方程线性化。假设振动状态的响应频率和风速为{\Omega }_{\mathrm{0}}和U_{\mathrm{0}}， 相邻状态可以通过添加相应的小增量定义为\Omega ={\Omega }_{\mathrm{0}}+\mathrm{\Delta }\Omega和U=U_{\mathrm{0}}+\mathrm{\Delta }U。

IHB 方法的第二步是谐波平衡过程， 即 Galerkin 过程。近似周期解Y_{\mathrm{0}}和相关的微小增量\mathrm{\Delta }\mathit{Y}用前n项傅里叶级数表示为

\begin{array}{l}{y}_{i\mathrm{0}}\approx a_{i\mathrm{0}}+{\displaystyle \sum _{j=\mathrm{1}}^n}(a_{ij}cos(j\tau )+b_{ij}sin(j\tau ))={\mathit{C}}_s{\mathit{a}}_i^{\mathrm{T}},\\ \mathrm{\Delta }{y}_i\approx \mathrm{\Delta }{a}_{i\mathrm{0}}+{\displaystyle \sum _{j=\mathrm{1}}^n}\left(\Delta \right.a_{ij}cos(j\tau )+\Delta b_{ij}sin(j\tau ))={\mathit{C}}_s\mathrm{\Delta }{\mathit{a}}_i^{\mathrm{T}},\end{array}（6）

式中： i=\mathrm{1,2},\mathrm{3,4}。j表示近似周期解的调和阶数。

{\mathit{C}}_s=[\begin{array}{cccccc}\mathrm{1}& \mathrm{c}\mathrm{o}\mathrm{s}(\tau )& \mathrm{s}\mathrm{i}\mathrm{n}(\tau )& \cdots & \mathrm{c}\mathrm{o}\mathrm{s}(n\tau )& \mathrm{s}\mathrm{i}\mathrm{n}(n\tau )\end{array}]，

{\mathit{a}}_i=[\begin{array}{cccccc}{a}_{i\mathrm{0}}& a_{i\mathrm{1}}& b_{i\mathrm{1}}& \cdots & a_{in}& b_{in}\end{array}]，

以及

\mathrm{\Delta }{\mathit{a}}_i=[\begin{array}{cccccc}\mathrm{\Delta }{a}_{i\mathrm{0}}& \mathrm{\Delta }{a}_{i\mathrm{1}}& \mathrm{\Delta }{b}_{i\mathrm{1}}& \cdots & \mathrm{\Delta }{a}_{in}& \mathrm{\Delta }{b}_{in}\end{array}]。

记 \mathit{S}=\left[\begin{array}{cccc}{\mathit{C}}_s& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& {\mathit{C}}_s& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& {\mathit{C}}_s& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& {\mathit{C}}_s\end{array}\right]，

响应的系数向量\mathit{A}={[\begin{array}{cccc}{\mathit{a}}_{\mathrm{1}}& {\mathit{a}}_{\mathrm{2}}& {\mathit{a}}_{\mathrm{3}}& {\mathit{a}}_{\mathrm{4}}\end{array}]}^{\mathrm{T}}， 小增量的系数向量\mathrm{\Delta }\mathit{A}={[\mathrm{\Delta }\begin{array}{cccc}{\mathit{a}}_{\mathrm{1}}& \mathrm{\Delta }{\mathit{a}}_{\mathrm{2}}& \mathrm{\Delta }{\mathit{a}}_{\mathrm{3}}& \mathrm{\Delta }{\mathit{a}}_{\mathrm{4}}\end{array}]}^{\mathrm{T}}。

使用Galerkin平均法平衡方程（5）中的谐波项。

\begin{array}{c}\delta \mathrm{\Delta }{\mathit{A}}^{\mathrm{T}}\left\{{\int }_{\mathrm{0}}^{\mathrm{2}\pi }{\mathit{S}}^T[{{\Omega }_{\mathrm{0}}}^{\mathrm{2}}\mathit{M}{\mathit{S}}^{″}+{\mathit{C}}_c{\mathit{S}}^{\text{'}}+{\mathit{K}}_k\mathit{S}]\mathrm{d}\tau \right\}\mathrm{\Delta }\mathit{A}=\\ -\delta \mathrm{\Delta }{\mathit{A}}^{\mathrm{T}}{\int }_{\mathrm{0}}^{\mathrm{2}\mathrm{\pi }}{\mathit{S}}^{\mathrm{T}}\mathit{R}\mathrm{d}\tau -\delta \mathrm{\Delta }{\mathit{A}}^{\mathrm{T}}{\int }_{\mathrm{0}}^{\mathrm{2}\mathrm{\pi }}{\mathit{S}}^{\mathrm{T}}{\mathit{\Omega }}_{\mathrm{\Omega }}\mathrm{\Delta }\Omega \mathrm{d}\tau -\\ \delta \mathrm{\Delta }{\mathit{A}}^{\mathrm{T}}{\int }_{\mathrm{0}}^{\mathrm{2}\mathrm{\pi }}{\mathit{S}}^{\mathrm{T}}{\mathit{U}}_U\mathrm{\Delta }U\mathrm{d}\tau ,\end{array}（7）

其中，

{\mathit{C}}_c={\Omega }_{\mathrm{0}}\mathit{C}+{\Omega }_{\mathrm{0}}{U}_{\mathrm{0}}{\mathit{C}}_u+{\Omega }_{\mathrm{0}}{U}_{\mathrm{0}}{\mathit{N}}_{\mathrm{13}},

{\mathit{K}}_k=\mathit{K}+{U_{\mathrm{0}}}^{\mathrm{2}}\mathit{K}{K}_u+{\Omega }_{\mathrm{0}}{U}_{\mathrm{0}}{\mathit{N}}_{\mathrm{12}},

\mathit{R}={\Omega }_{\mathrm{0}}^{\mathrm{2}}\mathit{M}{{\mathit{Y}}_{\mathrm{0}}}^{\mathrm{\text{'}}\mathrm{\text{'}}}+{\Omega }_{\mathrm{0}}\mathit{C}{{\mathit{Y}}_{\mathrm{0}}}^{\mathrm{\text{'}}}+{\Omega }_{\mathrm{0}}{U}_{\mathrm{0}}{\mathit{C}}_u{{\mathit{Y}}_{\mathrm{0}}}^{\mathrm{\text{'}}}+\mathit{K}{Y}_{\mathrm{0}}+{U_{\mathrm{0}}}^{\mathrm{2}}{\mathit{K}}_u{\mathit{Y}}_{\mathrm{0}}+{\Omega }_{\mathrm{0}}{U}_{\mathrm{0}}{\mathit{N}}_{\mathrm{11}},

{\mathit{U}}_U={\Omega }_{\mathrm{0}}{\mathit{C}}_u{{\mathit{Y}}_{\mathrm{0}}}^{\mathrm{\text{'}}\mathrm{\text{'}}}+\mathrm{2}{U}_{\mathrm{0}}{\mathit{K}}_u{\mathit{Y}}_{\mathrm{0}}+{\Omega }_{\mathrm{0}}{\mathit{N}}_{\mathrm{11}},

{\mathit{\Omega }}_{\mathrm{\Omega }}=\mathrm{2}{\Omega }_{\mathrm{0}}\mathit{M}{{\mathit{Y}}_{\mathrm{0}}}^{\mathrm{\text{'}}\mathrm{\text{'}}}+\mathit{C}{{\mathit{Y}}_{\mathrm{0}}}^{\mathrm{\text{'}}}+U_{\mathrm{0}}{\mathit{C}}_u{{\mathit{Y}}_{\mathrm{0}}}^{\mathrm{\text{'}}}+U_{\mathrm{0}}{\mathit{N}}_{\mathrm{11}},

{\mathit{N}}_{\mathrm{11}}=\left[\begin{array}{c}\mathrm{0}\\ \mathrm{0}\\ {\beta }_w({y_{\mathrm{3}}}^{\mathrm{2}}{y_{\mathrm{3}}}^{\text{'}}{)}_{\mathrm{0}}\\ \mathrm{0}\end{array}\right]，

{\mathit{N}}_{\mathrm{12}}=\left[\begin{array}{cccc}\mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{2}{\beta }_w(y_{\mathrm{3}}{y_{\mathrm{3}}}^{\text{'}}{)}_{\mathrm{0}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}\right]，

{\mathit{N}}_{\mathrm{13}}=\left[\begin{array}{cccc}\mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& {\beta }_w({y_{\mathrm{3}}}^{\mathrm{2}}{)}_{\mathrm{0}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}\right]\mathrm{。}

则方程（7）可以写为

{\mathit{J}}_a\mathrm{\Delta }\mathit{A}=-\overline{\mathit{R}}-{\mathit{J}}_{\mathrm{\Omega }}\mathrm{\Delta }\Omega -{\mathit{J}}_U\mathrm{\Delta }U，

{\mathit{J}}_a={\int }_{\mathrm{0}}^{\mathrm{2}\mathrm{\pi }}{\mathit{S}}^T[{\Omega }_{\mathrm{0}}^{\mathrm{2}}\mathit{M}{\mathit{S}}^{″}+{\mathit{C}}_c{\mathit{S}}^{\text{'}}+{\mathit{K}}_k\mathit{S}]\mathrm{d}\tau ,（8）

\overline{\mathit{R}}={\int }_{\mathrm{0}}^{\mathrm{2}\mathrm{\pi }}{\mathit{S}}^T\mathit{R}\mathrm{d}\tau， {\mathit{J}}_{\Omega }={\int }_{\mathrm{0}}^{\mathrm{2}\mathrm{\pi }}{\mathit{S}}^T{\mathit{\Omega }}_{\mathrm{\Omega }}\mathrm{d}\tau，

{\mathit{J}}_U={\int }_{\mathrm{0}}^{\mathrm{2}\mathrm{\pi }}{\mathit{S}}^T{\mathit{U}}_{\mathrm{U}}\mathrm{d}\tau \mathrm{。}

方程（8）中的未知数数量大于方程数量， 需要固定一个不参与迭代的参数。在本文中， 方程（6）中的b_{\mathrm{11}}=\mathrm{0}， \mathrm{\Delta }{b}_{\mathrm{11}}=\mathrm{0}， 位移的谐波相位被取作基准相位。那么， 方程（8）可以改写为

{\tilde{\mathit{J}}}_a\mathrm{\Delta }\tilde{\mathit{A}}=\overline{\mathit{R}}-{\mathit{J}}_U\mathrm{\Delta }U,（9）

其中， {\overline{\mathit{J}}}_a=[\begin{array}{cccccc}{\mathit{J}}_{a,\{\mathrm{1}\}}& {\mathit{J}}_{a,\{\mathrm{2}\}}& {\mathit{J}}_{\Omega }& {\mathit{J}}_{a,\{\mathrm{4}\}}& \cdots & {\mathit{J}}_{a,\{\mathrm{8}n+\mathrm{4}\}}\end{array}]，

\mathrm{\Delta }\tilde{A}={[\begin{array}{cccccc}\mathrm{\Delta }{\mathit{A}}_{\{\mathrm{1}\}}& \mathrm{\Delta }{\mathit{A}}_{\{\mathrm{2}\}}& \mathrm{\Delta }\Omega & \mathrm{\Delta }{\mathit{A}}_{\{\mathrm{4}\}}& \cdots & \mathrm{\Delta }{\mathit{A}}_{\{\mathrm{8}n+\mathrm{4}\}}\end{array}]}^{\mathrm{T}}\mathrm{。}

非线性系统频率响应中容易在某临界处发生突然跳跃， 形成尖峰点， 增量谐波平衡法迭代在这些点处往往难以收敛。为了提升其计算效率及其在尖锋点处的收敛性， 本文采用弧长延拓方法结合三次外插技术， 对下一迭代点的初始解进行预测。设x_{\mathrm{0}}， x_{\mathrm{1}}， x_{\mathrm{2}}和x_{\mathrm{3}}为已求得的解， 则两点之间的弧长可以表示为

s_i=\left|x_i-x_{i-\mathrm{1}}\right|=\sqrt[]{{\displaystyle \sum _{i=\mathrm{1}}^{\mathrm{6}n+\mathrm{1}}}[x_i(j)-x_{j-\mathrm{1}}{(j)]}^{\mathrm{2}}},(i=\mathrm{1,2},\mathrm{3,4},\cdots )。设下一迭代点x_{\mathrm{4}}的弧长增量为\mathrm{\Delta }s， 则下一迭代点的初始解可表示为

x_{\mathrm{4}}={{\displaystyle \sum _{i=\mathrm{0}}^{\mathrm{3}}}\left({\displaystyle \prod _{\begin{array}{l}i=\mathrm{0}\\ j\ne i\end{array}}^{\mathrm{3}}}\frac{t_{\mathrm{4}}-t_i}{t_i-t_j}\right)x}_i,

其中， t_{\mathrm{0}}=\mathrm{0}， t_{\mathrm{1}}=t_{\mathrm{0}}+s_{\mathrm{1}}， t_{\mathrm{2}}=t_{\mathrm{1}}+s_{\mathrm{2}}， t_{\mathrm{3}}=t_{\mathrm{2}}+s_{\mathrm{3}}， t_{\mathrm{4}}=t_{\mathrm{3}}+\mathrm{\Delta }s。应用弧长延拓法结合三次外插， 在每次增量谐波平衡法计算前给出周期解的预测值， 并将其作为Newton-Raphson 迭代的初始解， 可以保持并追踪算法的迭代轨迹， 同时缩短迭代的收敛时间。

2.2 周期解的稳定性

本文使用Floquet理论^［29］来确定周期解的稳定性。记\mathit{Z}=［\mathrm{\Delta }{y}_{\mathrm{1}}\mathrm{\Delta }{y}_{\mathrm{2}}\mathrm{\Delta }{y}_{\mathrm{3}}\mathrm{\Delta }{y_{\mathrm{1}}}^{\mathrm{\text{'}}}\mathrm{\Delta }{y_{\mathrm{2}}}^{\mathrm{\text{'}}}\mathrm{\Delta }{y_{\mathrm{3}}}^{\mathrm{\text{'}}}\mathrm{\Delta }{y}_{\mathrm{4}}］ 为方程（5）周期解附近的扰动向量， 残差向量\mathit{R}等于零。通过在方程周期解附近进行泰勒展开， 得到了方程的周期解， 忽略小扰动的非线性项， 可以得到关系式

{\mathit{Z}}^{\text{'}}=\mathit{Q}(\tau )\mathit{Z},（10）

其中，

\begin{array}{c}\mathit{Q}(\tau )=\\ \left[\begin{array}{ccc}{\mathrm{0}}_{\mathrm{3}\times \mathrm{3}}& {\mathit{I}}_{\mathrm{3}\times \mathrm{3}}& {\mathrm{0}}_{\mathrm{3}\times \mathrm{1}}\\ {\left(-\frac{\mathrm{1}}{{\Omega }_{\mathrm{0}}^{\mathrm{2}}}{\mathit{M}}_{\mathrm{0}}^{-\mathrm{1}}{\mathit{K}}_{\mathrm{0}}\right)}_{\mathrm{3}\times \mathrm{3}}& {\left(-\frac{\mathrm{1}}{{\Omega }_{\mathrm{0}}^{\mathrm{2}}}{\mathit{M}}_{\mathrm{0}}^{-\mathrm{1}}{\mathit{C}}_{\mathrm{0}}\right)}_{\mathrm{3}\times \mathrm{3}}& -\frac{\mathrm{1}}{{\Omega }_{\mathrm{0}}^{\mathrm{2}}}{\mathit{M}}_{\mathrm{0}}^{-\mathrm{1}}{\mathit{B}}_{\mathrm{0}}\\ {\mathrm{0}}_{\mathrm{1}\times \mathrm{3}}& \begin{array}{ccc}\mathrm{0}& -\theta /C_p& \mathrm{0}\end{array}& -\frac{\mathrm{1}}{{\Omega }_{\mathrm{0}}R{C}_P}\end{array}\right],\end{array}

{\mathit{M}}_{\mathrm{0}}=\left[\begin{array}{ccc}\mathrm{1}+m{{\lambda }_{Lb}}^{\mathrm{2}}& {\lambda }_{Lb}{\lambda }_{Lbn}& \mathrm{0}\\ {\lambda }_{Lb}{\lambda }_{Lbn}& \mathrm{1}& \mathrm{0}\\ -{\beta }_m{\lambda }_{Lb}& \mathrm{0}& \mathrm{1}\end{array}\right]， {\mathit{B}}_{\mathrm{0}}=\left[\begin{array}{r}\mathrm{0}\\ -\theta \\ \mathrm{0}\end{array}\right],

{\mathit{K}}_{\mathrm{0}}=\left[\begin{array}{ccc}{w}_{\mathrm{11}}^{\mathrm{2}}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& w_{\mathrm{21}}^{\mathrm{2}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}\right]+U_{\mathrm{0}}^{\mathrm{2}}\left[\begin{array}{ccc}\mathrm{0}& \mathrm{0}& -C_L{\lambda }_{Lb}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& {\beta }_s^{\mathrm{2}}\end{array}\right]+{\Omega }_{\mathrm{0}}{U}_{\mathrm{0}}\left[\begin{array}{ccc}\mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{2}{\beta }_w(y_{\mathrm{3}}{y_{\mathrm{3}}}^{\mathrm{\text{'}}}{)}_{\mathrm{0}}\end{array}\right],

\begin{array}{c}{\mathit{C}}_{\mathrm{0}}={\Omega }_{\mathrm{0}}\left[\begin{array}{ccc}\mathrm{2}{\zeta }_{\mathrm{11}}{w}_{\mathrm{11}}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{2}{\zeta }_{\mathrm{21}}{w}_{\mathrm{21}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}\right]+\\ {\Omega }_{\mathrm{0}}{U}_{\mathrm{0}}\left[\begin{array}{ccc}{{\lambda }_{Lb}}^{\mathrm{2}}{C}_m& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& -{\beta }_m\end{array}\right]+{\Omega }_{\mathrm{0}}{U}_{\mathrm{0}}\left[\begin{array}{ccc}\mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& {\beta }_w({y_{\mathrm{3}}}^{\mathrm{2}}{)}_{\mathrm{0}}\end{array}\right]\mathrm{。}\end{array}

通过在［0， 2π］上对方程（10）进行积分并以单位矩阵的不同列作为初始状态， 可以获得转移矩阵。转移矩阵的特征值被称为Floquet乘子。当所有的Floquet乘子都在单位圆内时， 获得稳定周期解， 反之得到不稳定解。

3 系统动力学的研究

本节讨论线性压电能量收集吸振器在控制涡激振动方面的有效性， 以及线性吸振器的参数对振动抑制和能量收集的影响。

3.1 线性吸振器控制动力系统的有效性

为了突出减振效果， 首先针对动力学控制方程（3）， 使用Runge-Kutta （ode45）数值方法分别求出有无吸振器下VIV系统的数值解， 然后根据表 1 数据进行仿真， 结果如图 2 所示。由图 2 可以看出， 无吸振器时， 结构固有频率由于只在锁频区域（4~6.2 m/s）的风速下与漩涡脱落频率同步， 因而产生了大幅振荡。为了在抑制结构大幅振荡的同时收集一些耗散的能量， 本文设计了一种压电能量收集吸振器放置于主结构的内部。由图 2结果可以看出， 添加压电能量收集吸振器后， 主结构的振幅被显著降低。

图3（a）为无减振时5.5 m/s风速下主结构的位移响应时域图， 图3（b）为有减振时5.5 m/s风速下主结构和压电能量收集吸振器的位移响应时域图， 虽然图中主结构的振幅降低， 但压电能量收集吸振器却有着较大的振幅， 同时压电能量收集吸振器也会因正压电效应而产生一定电能， 从而实现了抑制振动和收集能量的双重功能。

3.2 周期解、 多解和频率锁定现象

通过采用增量谐波平衡法（IHB）对涡激振动（VIV）系统的动力学方程进行求解， 能够有效揭示其多解共存特性。结果表明， 在锁频（Lock-in）区域中， 系统存在显著的多解共存现象， 这些解在相同参数下并存。当来流速度处于3.89~4.02 m/s和6.13~6.32 m/s两个临界区间时， 同一流速下可对应两种不同振幅的稳态解， 如图 4（a） 所示。虚线表示不稳定解， 实线表示稳定解。多解现象的出现是由系统的非线性引起的， 激励频率和振幅的变化会导致系统发生鞍结分岔， 使原本单一解分裂为多个解， 利用IHB法求解可以看出多解共存的情况。当初始点位于某吸引盆内时， 轨迹将渐近趋向对应吸引子， 且相空间中各吸引子的稳定流形形成吸引盆边界， 反映了不同初始条件（如初始位移、 速度）下系统最终收敛的稳态。

图 5 给出了速度为3.96和6.2 m/s时未装减振器情况下的吸引盆。位移和速度初始值分别在［-0.04，0.04］， ［-2，2］中选取， 电压的初值固定为0， 取最后10%时间段的响应来判断稳态输出特性。

图 5 中， 在蓝色区域内将收敛到较大振幅的周期解， 而在红色区域内将收敛到较小振幅周期解。

用增量谐波平衡法求出安装减振器后VIV系统的周期解， 结果如图 6（a） 所示， 对比数值方法和IHB方法的解可以看出， 两条曲线走势较为相似， 都呈现双峰特性， 两者最大相对误差小于10%。图 6（b） 显示了风速5.5 m/s时的位移时域曲线， 可以看出位移曲线均呈周期形式。

图 7（a） 中给出了通过数值方法和IHB法获得的频率特性， 可以看出两者吻合较好。未安装吸振器时， 涡激振动会发生锁频现象， 而安装线性吸振器后， 由于减振结构与其他结构间的相互作用， 响应频率不会被“锁定”到固定频率， 所以仍然呈现出随风速的递增而递增的变化趋势， 如图 7（b） 所示。

3.3 不同参数下的减振效果对比

本节研究线性减振器的各种参数， 即尖端钝体质量、 内部压电梁的长度、 宽度和厚度的影响规律， 主要从主结构的幅值和吸振器的电压方面来比较对减振效果的影响。已有研究表明， 吸振器末端附加质量在实现振动抑制效能中具有重要的影响， 但其减振性能并非随着质量的增加而单调提升， 只有当尖端质量与内梁的结构参数相匹配时， 方能起到良好的抑振效果。

图 8（a） 给出了压电梁长0.18 m， 宽0.02 m， 厚0.000 5 m时不同尖端质量情况下主结构的振幅。由于图中各曲线之间的差异并不十分显著， 难以仅凭直观观察明确判断出哪种情况下的减振效果更好。因此， 有必要借助具体的计算分析手段来判断， 本文引入平均位移思想对不同M_{\mathrm{t}}值下的曲线进行量化处理， 从而更为精准地评估各曲线所代表的减振效果。具体而言， 可按以下步骤运用平均位移法进行分析： 首先， 确定风速的计算区间， 图中风速范围为3.4~9.0 m/s， 将此作为统一的分析区间然后， 获取各曲线在该区间内不同风速点对应的位移数据， 并计算平均位移。设M_{\mathrm{t}}取不同值时对应的曲线位移数据分别为y_{\mathrm{1}i},y_{\mathrm{2}i},y_{\mathrm{3}i},\dots（i表示风速点序号）， 对于某一M_{\mathrm{t}}值对应的曲线， 其平均位移\overline{y}=\frac{\mathrm{1}}{n}{\displaystyle \sum _{i=\mathrm{1}}^n}{y}_i（n为该区间内风速点的数量）。平均位移值越小， 表明在该M_{\mathrm{t}}取值下， 结构在整个风速区间内的总体振动位移响应越小， 也就意味着对应的减振效果越好。

图 8（b） 给出了压电梁长0.18 m， 宽0.02 m， 厚0.000 5 m时不同尖端质量情况下吸振器输出电压的大小。采用同样的方法判断能量收集的效果， 平均电压值越大， 表明在该M_{\mathrm{t}}取值下， 结构在整个风速区间内的总体电压越大， 也就意味着对应的能量收集效果越好。本研究以涡激振动抑制为核心目标， 故优先致力于主结构振幅的最小化调控。在确保涡激振动抑制效果显著的前提下， 进一步探索能量收集潜力， 以实现压电元件输出电压的最大化。计算结果表明， M_{\mathrm{t}}=0.0388 87 kg时的减振效果较好， 同时也有一定的电压输出， 为了计算方便， 在进一步的研究中考虑M_{\mathrm{t}}=0.04 kg。

压电能量收集吸振器的梁长度、 宽度和厚度对减振和能量输出的影响结果如图 9 所示。

尖端质量M_{\mathrm{t}}为0.04 kg， 内梁宽 0.02 m， 厚 0.000 5 m时， 不同内梁长度下的位移幅值和电压输出分别如图 9（a） 和图 9（b） 所示， 可以看出， 随着梁长度的增加， 吸振器减振效果会有一定提升， 但吸振器电压输出也会降低。这是因为长度增大导致吸振器的频率趋向涡激振动装置的一阶固有频率， 在锁频范围内产生较好的振动抑制效果。但梁长度并不是越大越好， 长度较大的梁会导致吸振器的大幅相对位移， 增大了碰到主结构内壁的可能， 从而影响结构动力学响应。在进行振动抑制参数调控时， 应当合理选择吸振器的长度， 在抑制涡激振动主结构振幅的同时， 合理保留必要的振动能量用于机电转换， 形成“振动抑制优先、 能量收集协同”的双目标优化框架。图 9（c） 给出了内梁长度0.18 m， 厚度0.000 5 m， M_{\mathrm{t}}=0.04 kg时不同内梁宽度下的位移幅值， 可以看出， 在低风速时减振效果较差， 高风速时减振效果较好。这是因为脱涡频率与来流风速之间存在正相关关系， 低风速时的脱涡频率较低， 而较小的内梁宽度也导致减振器较低的固有频率， 两者之间存在频率同步的条件。相反， 较大的内梁宽度也导致减振器较高的固有频率， 只能在较大的风速下才会产生较好的能量收集效果。同样， 当减振器的厚度增大时， 也会显著增大固有频率， 使得系统在低风速时的减振效果较差， 而高风速时的减振效果较好， 如图 9（e） 所示。图 9（d） 给出了内梁长0.18 m， 厚0.0005 m， M_{\mathrm{t}}为0.04 kg时不同内梁宽度的电压幅值。可以看出， 在低风速下， 较小的宽度导致大幅的响应输出； 在高风速下， 较大的宽度导致较大的响应输出。图 9（f） 给出了不同内梁厚度的电压幅值， 相应输出随着内梁厚度的增大而增大。从振动抑制和能量收集双重目标来看， 适当增加减振压电梁的厚度， 有利于降低特定风速下的振动幅值， 同时能够提高电压输出。

4 结　论

本文提出了一种通过压电能量收集吸振器实现涡激振动抑制的被动控制方案， 并且采用IHB法求解了具有线性减振器的涡激振动的半解析周期解， 进一步绘制了响应位移幅值、 电压幅值和响应频率的曲线。研究表明： 在涡激振动的频率锁定区域中， 系统存在显著的多解共存现象； 与传统涡激振动不同， 具有线性吸振器的涡激振动系统的响应频率不会被“锁定”到固定频率， 而是会呈现为单调变化的频率曲线。通过量化分析减振器参数（尖端钝体质量、 内梁的长度、 宽度和厚度）对主结构振幅及输出电压的影响， 本文明确了减振器的最优参数区间， 为工程设计提供了理论依据。此外， 本文将增量谐波平衡法应用于涡激振动减振模型， 有助于更好地理解涡激振动的多解和频率锁定问题。后续将进一步探索非线性减振器的涡激振动抑制问题。

参考文献

[1] 马烨璇， 宋志友， 徐万海. 基于能量传递规律的海洋立管涡激振动抑制研究［J］. 力学学报， 2022， 54（4）： 901-911.

[2] MA Yexuan， SONG Zhiyou， XU Wanhai. Study on vortex-induced vibration suppression of marine riser based on energy transfer［J］.Chinese Journal of Theoretical and Applied Mechanics， 2022， 54（4）： 901-911. （in Chinese）

[3] 李佳鸿， 李正良， 王涛. 基于神经网络的输电塔钢管构件涡激振动幅值预测方法［J］. 工程力学， 2024， 41（1）： 64-75.

[4] LI Jiahong， LI Zhengliang， WANG Tao. Prediction method for vortex-induced vibration amplitude of steel tubes in transmission towers based on neural network［J］. Engineering Mechanics， 2024， 41（1）： 64-75.（in Chinese）

[5] 肖志鹏， 仇翯辰， 周磊. 复合材料支撑机翼撑杆位置与结构综合优化设计［J］. 工程力学， 2019， 36（9）： 213-220.

[6] XIAO Zhipeng， QIU Hechen， ZHOU Lei. Integrated optimization design of strut location and structure for composite strut-braced wing［J］. Engineering Mechanics， 2019， 36（9）： 213-220. （in Chinese）

[7] QIN G F， ZHANG H X， LI D. Numerical study on vortex induced vibration of hydrofoils with trailing-edge truncation［J］. Ocean Engineering， 2023， 275： 114083.

[8] ZHANG J K， HUANG B， ZHANG M J， et al. Investigation of vortex-induced vibration characteristics around the NACA0009 hydrofoil with a focus on the lock-in phenomenon［J］. Ocean Engineering， 2024， 312： 119082.

[9] BEARMAN P W. Circular cylinder wakes and vortex-induced vibrations［J］. Journal of Fluids and Structures， 2011， 27（5/6）： 648-658.

[10] 周佳豪， 马文勇， 黄伯城. 临界雷诺数区光滑圆柱振动与气动力研究［J］. 工程力学， 2019， 36（S1）： 306-310.

[11] ZHOU Jiahao， MA Wenyong， HUANG Bocheng. Study on vibration and aerodynamic forces of a smooth cylinder in the critical Reynolds number range［J］. Engineering Mechanics， 2019， 36（S1）： 306-310. （in Chinese）

[12] LAI Z H， WANG S B， ZHU L K， et al. A hybrid piezo-dielectric wind energy harvester for high-performance vortex-induced vibration energy harvesting［J］. Mechanical Systems and Signal Processing， 2021 150（1）： 107212.

[13] 胡志康， 鲁丽， 王琳， 等. 错列双管涡激振动风洞实验研究［J］. 化工机械， 2018， 45（6）： 675-679.

[14] HU Zhikang， LU Li， WANG Lin， et al. Wind tunnel experiment on vortex-induced vibrationof staggered double circular tube［J］. Chemical Machinery， 2018， 45（6）： 675-679. （in Chinese）

[15] ZHANG H， FAN B C， CHEN Z H， et al. Numerical study of the suppression mechanism of vortex-induced vibration by symmetric Lorentz forces［J］. Journal of Fluids and Structures， 2014， 48： 62-80.

[16] AKHTAR I， NAYFEH A. On controlling the bluff body wake using a reduced-order model［C］//4th Flow Control Conference， 2008： 4189.

[17] DAI H L， ABDELKEFI A， WANG L， et al. Time-delay feedback controller for amplitude reduction in vortex-induced vibrations［J］. Nonlinear Dynamics， 2015， 80（1）： 59-70.

[18] MEHMOOD A， ABDELKEFI A， AKHTAR I， et al. Linear and nonlinear active feedback controls for vortex-induced vibrations of circular cylinders［J］. Journal of Vibration and Control， 2014， 20（8）： 1137-1147.

[19] ZHOU S， CAO J， WANG W， et al. Modeling and experimental verification of doubly nonlinear magnet-coupled piezoelectric energy harvesting from ambient vibration［J］. Smart Materials and Structures， 2015， 24（5）： 055008.

[20] 袁江波， 谢涛， 单小彪， 等. 压电俘能技术研究现状综述［J］. 振动与冲击， 2009， 28（10）： 36-42.

[21] YUAN Jiangbo， XIE Tao， SHAN Xiaobiao， et al. A review of current situation for piezoelectric energy harvesting［J］. Journal of Vibration and Shock， 2009， 28（10）： 36-42. （in Chinese）

[22] LIU G， LU Z R， LI W， et al. A new semi-analytical technique for nonlinear systems based on response sensitivity analysis［J］. Nonlinear Dynamics， 2021， 103（2）： 1529-1551.

[23] LOGHMAN E， KAMALI A， BAKHTIARI-NEJAD F， et al. On the combined shooting-pseudo-arclength method for finding frequency response of nonlinear fractional-order differential equations［J］. Journal of Sound and Vibration， 2022， 516： 116521.

[24] PANYAM M， DAQAQ M F. Characterizing the effective bandwidth of tri-stable energy harvesters［J］. Journal of Sound and Vibration， 2017， 386： 336-358.

[25] NIU J C， LI X F， XING H J. Superharmonic resonance of fractional-order mathieu-duffing oscillator［J］. Journal of Computational and Nonlinear Dynamics， 2019， 14（7）： 071005.

[26] ZHOU S X， CAO J Y， LIN J. Theoretical analysis and experimental verification for improving energy harvesting performance of nonlinear monostable energy harvesters［J］. Nonlinear Dynamics， 2016， 86（3）： 1599-1611.

[27] LAU S L， CHEUNG Y K. Amplitude incremental variational principle for nonlinear vibration of elastic systems［J］. Journal of Applied Mechanics， 1981， 48（4）： 959-964.

[28] LAU S L， CHEUNG Y K， WU S Y. Incremental harmonic balance method with multiple time scales for aperiodic vibration of nonlinear systems［J］. Journal of Applied Mechanics， 1983， 50（4a）： 871-876.

[29] SKOP R A， GRIFFIN O M. A model for the vortex-excited resonant response of bluff cylinders［J］. Journal of Sound and Vibration， 1973， 27（2）： 225-233.

[30] IWAN W D， BLEVINS R D. A model for vortex-induced oscillation of structures［J］. Journal of Applied Mechanics， 1974， 41（3）： 581-586.

[31] HARTLEN R T， CURRIE I G. Lift-oscillator model of vortex-induced vibration［J］. Journal of the Engineering Mechanics Division， 1970， 96（5）： 577-591.

[32] DECUYPER J， DE TROYER T， RUNACRES M C， et al. Nonlinear state-space modelling of the kinematics of an oscillating circular cylinder in a fluid flow［J］. Mechanical Systems and Signal Processing， 2018， 98： 209-230.

[33] HUYNH B H， TJAHJOWIDODO T， ZHONG Z W， et al. Design and experiment of controlled bistable vortex induced vibration energy harvesting systems operating in chaotic regions［J］. Mechanical Systems Signal Processing， 2018， 98（1）： 1097-1115.

[34] FACCHINETTI M L， DE LANGRE E， BIOLLEY F. Coupling of structure and wake oscillators in vortex-induced vibrations［J］. Journal of Fluids and Structures， 2004， 19（2）： 123-140.

[35] LIU G， LÜ Z R， LIU J K， et al. Quasi-periodic aeroelastic response analysis of an airfoil with external store by incremental harmonic balance method［J］. International Journal of Non-Linear Mechanics， 2018， 100： 10-19.

基金资助

国家自然科学基金资助项目(12372026)

山西省回国留学人员科研资助项目(2022-150)

山西省自然科学基金资助项目(202203021211086)