近十年来，利用风致振动的能量收集技术，为微型无人飞行器或低功耗无线传感器网络供电，受到了国内外学者的关注^[1-6]，具有巨大的工程应用前景。Erturk等利用颤振机理设计的能量收集器^[7]，在风速为9.3 m/s时发生极限环振动(Limit Cycle Oscillations, LCOs)，向负载电路输出电功率达到10.7 mW。Abdelkefi和Ghommem对扑翼飞行器的能量收集进行了数值研究，输出电功率超过50 mW，能够满足机载摄像镜头的供电需求^[8]。Bryant和Garcia将二元翼段铰接在压电梁上，设计了微型低风速颤振能量收集器^[9]，颤振始发风速仅为1.9 m/s，在2~9 m/s的风速区间收集到的电功率达到0.85~2.2 mW。Xiang等研究了压电机翼在突风载荷作用下的能量收集性能，数值分析了压电片的位置、厚度、机翼刚性轴位置等参数对输出电功率的影响^[10]。在已有的这些研究中，虽然不同形式的气动弹性能量收集器被提出，但普遍基于线性气动弹性理论来设计，这使得根据特定工作风速设计的能量收集器，当环境风速发生改变或结构存在不确定性时，系统难以持续工作在最佳的能量收集状态。

为解决上述问题，近年来，逐渐有学者借鉴振动能量收集中采用的非线性技术^[11-13]，来提升气动弹性能量收集系统的工作性能。Li等设计了一种磁力作用下的单稳态颤振能量收集器^[14]，能够降低系统的颤振始发风速，在研究的风速范围内，能量收集性能提升了一倍。Alhadidi等^[15]利用磁力产生双稳态，进行了尾涡驰振的双稳态能量收集研究，拓宽了结构对尾涡频率的响应带宽，从而降低能量收集器对风速变化的敏感度。Zhang等实验验证了磁力双稳态对涡激振动能量收集器输出电功率的增强作用^[16]，然而双稳态同时引起了系统工作风速的提高，作者对此没有进行深入的研究。Zhou等研究了带有末端分叉的压电梁在气流扰动下的双稳态能量收集^[17]，发现在不同的风速下，磁力双稳态可能激发能量收集器在势阱间(inter-well)跃迁，也可能将系统束缚在其中一个势阱内(intra-well)，从而对能量收集性能造成截然相反的影响。总的来看，与振动能量收集相比，非线性气动弹性能量收集的研究仍处于起步阶段，存在着一定的研究空白，非线性多稳态的影响机制缺乏深入的探索，并需要对多稳态带来的弊端提出改进措施。

本文设计了一类新型的磁耦合双稳态颤振能量收集器(Bi-stable Flutter Energy Harvester, BFEH)，对双稳态的影响机理以及存在的问题进行了探讨，从而提出变势能阱双稳态的改进策略(Variable-potential-well Bi-stable Flutter Energy Harvester, VBFEH)。建立了磁耦合气动弹性能量收集系统的动力学控制方程，利用非线性磁偶极模型和平衡点稳定性理论，对出现双稳态构型的磁力参数范围进行了求解，研究了双稳态构型对系统的颤振特性以及输出电功率的影响规律。本文提出的VBFEH通过弹性支撑下外部磁铁的运动来自适应调节系统的势能阱，激发双稳态势阱间运动，达到了同时拓宽能量收集器工作风速范围和提升输出电功率的目的，并解释了变势能阱技术对颤振能量收集性能的增强机理，进一步提高了系统对风场环境的适应性。

1 VBFEH的组成及工作原理

图 1为本文提出的变势能阱双稳态颤振能量收集器的结构示意图，其主要由2根在根部粘贴有压电片的悬臂梁和在梁梢部通过刚性轴铰接的二元翼段构成。图中, U_∞为来流风速；u(t)为外部磁铁的位移；α(t)为二元翼段的俯仰位移；R_L为负载电阻; h(t)为沉浮模态位移；V⁺和V^-分别为电压片上、下表面电极上的电势。与传统的单稳态颤振能量收集器^{[7, 9]}不同，本文提出的VBFEH在悬臂梁的梢部带有永磁铁，外部磁极方向相反的永磁铁安装在支撑弹簧上，通过内外部磁铁的磁场作用，对颤振能量收集器施加非线性磁力。根据文献[18]可知，当内外部磁铁间距D小于临界值时，排斥性磁力将抵抗悬臂梁的弹性恢复力，使能量收集器发生静态失稳，出现双稳态构型，如图 2所示，w(x, t)为梁的横向振动位移；k_e为外部磁铁支撑弹簧刚度。

在双稳态构型下，原始的零平衡点变得不稳定，而在两侧对称位置出现新的稳定平衡点，发生颤振后，能量收集器将围绕2个新的平衡位置，发生势能阱间的大幅值跃迁运动。外部磁铁在弹性支撑下的运动将实时改变与内部磁铁之间的相对位置，起到调节双稳态系统势能阱深度的作用^[19]，从而实现可变势能阱的双稳态颤振能量收集。显然，磁耦合双稳态、变势能阱双稳态的颤振能量收集特性，将不同于传统的颤振能量收集系统。本文重点研究了磁铁的间距、磁偶极矩以及外部磁铁的弹性运动对颤振能量收集系统工作风速范围、输出电功率的影响规律，旨在增强颤振能量收集性能、提高系统对风场环境的适应性。

2 磁力-压电-气动弹性耦合系统建模

本文提出的VBFEH的力学分析模型可等效为图 2所示，压电换能电路与气动弹性系统的沉浮自由度发生耦合。为了研究颤振能量收集器的动态特性，本节将根据Euler-Bernoulli梁的弯曲大变形理论，通过Hamilton变分原理、压电本构方程以及Kirchhoff电学基本定律，建立双稳态颤振压电能量收集系统的磁力-压电-气动弹性耦合的动力学控制方程。

2.1 理论模型

本文对带动二元翼段沉浮运动的2根压电梁采用Euler-Bernoulli梁理论进行建模，记压电梁的横向振动位移为w(x, t)，考虑梁在颤振运动中弯曲大变形引起的非线性刚度，根据几何非线性关系^[14]可得梁横截面上任意一点的应变为

$ {S_1} = - z\frac{{{\partial ^2}w}}{{\partial {x^2}}}\left[ {1 + \frac{1}{2}{{\left( {\frac{{\partial w}}{{\partial x}}} \right)}^2}} \right] $	（1）

式中：S₁为梁的弯曲应变；z为梁截面上的点到中性轴的距离。梁根部粘贴的压电片伴随梁发生弯曲变形，其弯曲本构方程为

$ \left\{ {\begin{array}{*{20}{l}} {{T_1} = c_{11}^E{S_1} - {e_{31}}{E_3}}\\ {{D_3} = {e_{31}}{S_1} + \varepsilon _{33}^S{E_3}} \end{array}} \right. $	（2）

式中：T₁为弯曲轴向应力；c₁₁^E为压电材料弹性模量；e₃₁为压电应力常数；E₃为电场强度，对于压电片结构, $ {E_3} = \frac{{V(t)}}{{{h_{\text{p}}}}}$，V(t)为压电片两端的电压；h_p为压电片厚度；D₃为电位移；ε₃₃^S为介电常数。利用压电本构方程式(2)，在2根压电梁上积分，颤振能量收集系统的势能可以写作：

$ \begin{array}{*{20}{l}} {{E_{\rm{p}}} = \frac{1}{2}\int\limits_{{V_{\rm{p}}}} {c_{11}^E} S_1^2{\rm{d}}{V_{\rm{p}}} - \frac{1}{2}\int\limits_{{V_{\rm{p}}}} {{S_1}} {e_{31}}{E_3}{\rm{d}}{V_{\rm{p}}} + }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \frac{1}{2}\int\limits_{{V_{\rm{b}}}} {{E_{\rm{b}}}} S_1^2{\rm{d}}{V_{\rm{b}}}} \end{array} $

（3）

式中：V_p、V_b代表压电片和悬臂梁的积分体积; E_b为梁的弹性模量。第1项为压电片的弹性应变势能，第2项代表压电片机电耦合势能，第3项为悬臂梁的弹性应变势能。压电片储存的内部电能^[3]为

$ {E_{{\rm{ie}}}} = \frac{1}{2}\int_{{V_{\rm{p}}}} {{D_3}} {E_3}{\rm{d}}{V_{\rm{p}}} $	（4）

考虑二元翼段的沉浮运动、绕刚性轴的俯仰运动，以及2根压电梁的振动，系统的总动能表达式为

$ \begin{array}{*{20}{l}} {{E_{\rm{k}}} = \frac{1}{2}\int_{{V_{\rm{b}}}} {{\rho _{\rm{b}}}} {{\dot w}^2}{\rm{d}}{V_{\rm{b}}} + \frac{1}{2}\int_{{V_{\rm{p}}}} {{\rho _{\rm{p}}}} {{\dot w}^2}{\rm{d}}{V_{\rm{p}}} + \frac{1}{2}{I_\alpha }{{\dot \alpha }^2} + }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {m_{\rm{F}}}{x_\alpha }b\dot \omega \delta (x - {L_{\rm{b}}})\dot \alpha + }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \frac{1}{2}({m_{\rm{F}}} + {m_{\rm{M}}}){{\dot w}^2}\delta (x - {L_{\rm{b}}})} \end{array} $

（5）

式中：ρ_p和ρ_b分布为压电片和梁的密度；m_F为二元翼段的质量；b为翼段的半弦长；x_α为重心在刚性轴后的无量纲距离；I_α为绕刚性轴的质量惯性矩；m_M为梁梢部磁铁的质量；δ(x)为狄拉克函数。

颤振能量收集器通过压电效应向负载电路输送电能，负载电路的电阻记为R_L。每根悬臂梁两侧粘贴的压电片首先经过串联连接，然后将2根压电梁的电极并联，将输出电压施加在负载2端，如图 1所示。系统受到的外界非保守力做功包括气动载荷对二元翼段做功W_a、外部磁场对梁梢部磁铁做功W_m以及压电片向负载电路输出的电能W_e，分别写作：

$ {{W_{\rm{a}}} = - {Q_{\rm{L}}}w\delta (x - {L_{\rm{b}}}) + {Q_{\rm{M}}}\alpha } $	（6）

$ {{W_{\rm{m}}} = {F_{\rm{m}}}w\delta (x - {L_{\rm{b}}})} $	（7）

$ {{W_{\rm{e}}} = QV} $	（8）

式中：气动载荷的作用点在刚性轴位置，Q_L和Q_M表示气动升力和气动力矩，气动升力的方向向上为正，气动力矩使翼段抬头为正方向；F_m为梁梢部磁铁受到的磁力；Q为通过负载电路的总电荷。根据广义Hamilton变分原理^[9]：

$ \int_{{t_1}}^{{t_2}} \delta ({E_{\rm{k}}} - {E_{\rm{p}}} + {E_{{\rm{ie}}}} + {W_{\rm{a}}} + {W_{\rm{m}}} + {W_{\rm{e}}}){\rm{d}}t = 0 $	（9）

将式(3)~式(8)代入式(9)，首先经过体积积分，再分别对广义自由度w、α和V进行变分运算，可以得到系统的动力学控制方程：

$ \begin{array}{l} \rho A\ddot w + ({m_{\rm{F}}} + {m_{\rm{M}}})\ddot w\delta (x - {L_{\rm{b}}}) + {m_{\rm{F}}}{x_\alpha }b\ddot \alpha + \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \begin{array}{*{20}{l}} {YI({w^{\prime \prime \prime \prime }} + {w^{\prime \prime }}^3 + 4{w^\prime }{w^{\prime \prime }}{w^{\prime \prime \prime }} + {w^{\prime 2}}{w^{\prime \prime \prime \prime }}) + }\\ {{{({z_{\rm{m}}}{e_{31}}{b_{\rm{p}}}\delta (x - {L_{\rm{p}}}))}^\prime }V = ( - {Q_{\rm{L}}} + {F_{\rm{m}}}) \cdot }\\ {\delta (x - {L_{\rm{b}}})} \end{array} \end{array} $

（10）

$ {{m_{\rm{F}}}{x_\alpha }b\ddot w\delta (x - {L_{\rm{b}}}) + {I_\alpha }\ddot \alpha = {Q_{\rm{M}}}} $	（11）

$ {\frac{{\varepsilon _{33}^S}}{{{h_{\rm{p}}}}}{b_{\rm{p}}}{L_{\rm{p}}}\dot V + \frac{V}{{{R_{\rm{L}}}}} - 2\int_0^{{L_{\rm{p}}}} {\frac{{{z_{\rm{m}}}{b_{\rm{p}}}{e_{31}}}}{{{h_{\rm{p}}}}}} {{\dot \omega }^{\prime \prime }}{\rm{d}}x = 0} $	（12）

式中：ρA表示梁单位长度的质量；YI为梁的弯曲截面惯性矩；ẇ、w′分别表示梁的位移对时间t和空间x求偏导；b_p、L_p为压电片的宽度和长度；z_m为压电片厚度中心到悬臂梁中性轴的距离。

能量收集系统的等效接口电路^[20]可以表示为图 3所示的电路图，记每个压电片输出的总电流为i_p(t)，压电片等效电容为C_p，等效电容两端通过的电流为i_c(t)。负载电阻的电压、电流分别为V(t)、i_L(t)。根据电路Kirchhoff定律可得：

$ {2({i_{\rm{p}}}(t) - {i_{\rm{c}}}(t)) = \frac{{V(t)}}{{{R_{\rm{L}}}}}} $	（13）

$ {\frac{{2\int {{i_{\rm{c}}}} (t){\rm{d}}t}}{{{C_{\rm{p}}}}} = V(t)} $	（14）

联立式(13)和式(14)可得系统的等效电路方程：

$ {C_{\rm{p}}}\dot V(t) + \frac{{V(t)}}{{{R_{\rm{L}}}}} - 2{i_{\rm{p}}}(t) = 0 $	（15）

利用Hamilton原理和Kirchhoff定律得到的电学控制方程等价，对比式(12)和式(15)可知：

$ {{C_{\rm{p}}} = \frac{{\varepsilon _{33}^S}}{{{h_{\rm{p}}}}}{b_{\rm{p}}}{L_{\rm{p}}}} $	（16）

$ {{i_{\rm{p}}}(t) = \int_0^{{L_{\rm{p}}}} {\frac{{{z_{\rm{m}}}{b_{\rm{p}}}{e_{31}}}}{{{h_{\rm{p}}}}}} {{\dot w}^{\prime \prime }}{\rm{d}}x} $	（17）

本文采用模态截断法^[9]对颤振能量收集系统进行模型降阶，考虑悬臂梁的一阶弯曲模态，设梁的横向振动位移为$w(x, t) \doteq {\phi _1}\left( x \right)h(t), {\phi _1}\left( x \right) $为梁的一阶模态函数，考虑结构中存在的阻尼，将梁的位移解以及式(16)和式(17)代入式(10)~式(12)，可得：

$ \begin{array}{l} M\ddot h + {m_{\rm{F}}}{x_\alpha }b\ddot \alpha + {c_h}\dot h + {K_{h1}}h + {K_{h3}}{h^3} + \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \theta V = - {Q_{\rm{L}}} + {F_{\rm{m}}} \end{array} $

（18）

$ {{m_{\rm{F}}}{x_\alpha }b\ddot h + {I_\alpha }\ddot \alpha + {c_\alpha }\dot \alpha = {Q_{\rm{M}}}} $	（19）

$ {{C_{\rm{p}}}\dot V + \frac{V}{{{R_{\rm{L}}}}} - 2\theta \dot h = 0} $	（20）

式中：梁的模态函数${\phi _1}\left( x \right) $以梁的梢部模态位移进行归一化，即$ {\phi _1}\left( {{L_{\text{b}}}} \right) = 1$；M为系统的沉浮质量；K_h1、K_h3表示沉浮自由度的线性和非线性刚度；c_h、c_α为沉浮和俯仰阻尼系数；θ=z_me₃₁b_p· ϕ′₁x(L_p)为系统的机电耦合系数^[7]。

2.2 气动力模型

在颤振能量收集器的动力学控制方程中，气动力模型采用Theodorsen非定常气动力的Wagner函数近似表达式，用来描述颤振能量收集器所受到的时域气动载荷。当风速为U时，二元翼段在任意运动形式下所受的气动升力、气动力矩可以表示为

$ \begin{array}{*{20}{c}} {{Q_{\rm{L}}}(t) = \pi \rho {b^2}s[\ddot h(t) + U\dot \alpha (t) - ba\ddot \alpha (t)] + }\\ {2\pi \rho Ubs {Q_{3/4}}(t)(1 - {A_1} - {A_2}) + }\\ {2\pi \rho {U^2}bs[{A_1}{b_1}{w_1}(t) + {A_2}{b_2}{w_2}(t)]} \end{array} $	（21）

$ \begin{array}{l} \begin{array}{*{20}{c}} {{Q_{\rm{M}}}(t) = \pi {\rho ^2}s\left[ {ba\ddot h(t) - Ub\left( {\frac{1}{2} - a} \right)\dot \alpha (t) - } \right.}\\ {\left. {{b^2}\left( {\frac{1}{8} + {a^2}} \right)\ddot \alpha (t)} \right] + 2\pi \rho Ubs\left( {\frac{1}{2} + a} \right) \cdot } \end{array}\\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \begin{array}{*{20}{l}} {{Q_{\frac{3}{4}}}(t)(1 - {A_1} - {A_2}) + 2\pi \rho {U^2}bs\left( {\frac{1}{2} + a} \right) \cdot }\\ {[{A_1}{b_1}{w_1}(t) + {A_2}{b_2}{w_2}(t)]} \end{array} \end{array} $

（22）

式中：s为二元翼段的展长；${Q_{\frac{3}{4}}}(t) = U\alpha (t) + \dot h(t) + b\left( {\frac{1}{2} - a} \right)\dot a(t) $为翼面上3/4弦长处的下洗速度；a为刚心位于翼段弦长中点后的无量纲距离，有量纲距离为ab; A₁和A₂为Wagner中的常数，有A₁=0.165，A₂=0.335；w₁(t)和w₂(t)为Wagner函数引入的空气动力状态变量，它们满足如下微分方程：

$ \begin{array}{l} \left[ {\begin{array}{*{20}{c}} {{{\dot w}_1}}\\ {{{\dot w}_2}} \end{array}} \right] = \frac{U}{b}\left[ {\begin{array}{*{20}{c}} { - {b_1}}&0\\ 0&{ - {b_2}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{w_1}}\\ {{w_2}} \end{array}} \right] + \frac{U}{b}\left[ {\begin{array}{*{20}{c}} 0&{1/b}\\ 0&1 \end{array}} \right]\left\{ {\begin{array}{*{20}{l}} h\\ \alpha \end{array}} \right\} + \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \left[ {\begin{array}{*{20}{c}} {\frac{1}{b}}&{\frac{{1/2 - a}}{b}}\\ 1&{\frac{1}{2} - a} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {\dot h}\\ {\dot \alpha } \end{array}} \right] \end{array} $

（23）

其中:b₁和b₂为Wagner函数中的常数，有b₁=0.045 5, b₂=0.30。

2.3 磁力建模与系统的双稳态构型分析

本文采用磁偶极理论对非线性磁力进行建模，来计算压电梁梢部磁铁以及外部磁铁受到的磁力，Dipole-dipole磁偶极模型为^[21]：

$ \begin{array}{*{20}{l}} {{\mathit{\boldsymbol{F}}_{{\rm{AB}}}} = - \frac{{{\mu _0}}}{{4\pi }}\nabla ((\nabla \frac{{{\mathit{\boldsymbol{m}}_{\rm{A}}} \cdot \mathit{\boldsymbol{r}}}}{{|\mathit{\boldsymbol{r}}{|^3}}}) \cdot {\mathit{\boldsymbol{m}}_{\rm{B}}}) = \frac{{3{\mu _0}}}{{4\pi |\mathit{\boldsymbol{r}}{|^5}}} \cdot }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} [(\mathit{\boldsymbol{r}} \times {\mathit{\boldsymbol{m}}_{\rm{A}}}) \times {\mathit{\boldsymbol{m}}_{\rm{B}}} + (\mathit{\boldsymbol{r}} \times {\mathit{\boldsymbol{m}}_{\rm{B}}}) \times {\mathit{\boldsymbol{m}}_{\rm{A}}} - }\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} 2\mathit{\boldsymbol{r}}({\mathit{\boldsymbol{m}}_{\rm{A}}} \cdot {\mathit{\boldsymbol{m}}_{\rm{B}}}) + 5\mathit{\boldsymbol{r}}(\mathit{\boldsymbol{r}} \times {\mathit{\boldsymbol{m}}_{\rm{A}}}) \cdot (\mathit{\boldsymbol{r}} \times {\mathit{\boldsymbol{m}}_{\rm{B}}})]} \end{array} $

（24）

式中：m_A、m_B为压电梁梢部和外部2个磁铁的磁偶极矩矢量，本文中内外部磁铁的磁偶极矩矢量方向相反，即m_A/m_A=-m_B/m_B，它们的大小(m_A和m_B)与磁铁的磁极强度、体积成正比；μ₀为真空中的磁场常数；▽(·)为方向梯度算子。为简化磁力模型，不考虑梁弯曲变形引起的梢部磁铁的转动，在颤振能量收集器的运动过程中，2个磁铁的几何位置构型如图 2所示。记弹簧支撑的外部磁铁的水平位移为u(t); 2个磁铁间的水平间距为D; 磁铁空间距离矢量为r，它的模可表示为$ \left| \mathit{\boldsymbol{r}} \right| = \sqrt {h{{(t)}^2} + {{\left( {D + u(t)} \right)}^2}} $，r与水平方向的夹角为$ \varphi = {\text{arctan}}\left( {\frac{{{h}(t)}}{{D + u(t)}}} \right)$。将2个磁铁的几何位置关系代入磁力模型，可得到颤振能量收集器在沉浮自由度上的磁力F_m以及在外部磁铁受到的磁力F_me，分别为

$ {F_{\rm{m}}}(h,u,D) = \frac{{3{\mu _0}{m_{\rm{A}}}{m_{\rm{B}}}(4{\rm{sin}}{\kern 1pt} {\kern 1pt} \varphi - 5{{({\rm{sin}}{\kern 1pt} {\kern 1pt} \varphi )}^3})}}{{4\pi {{(h{{(t)}^2} + {{(D + u(t))}^2})}^2}}} $	（25）

$ \begin{array}{*{20}{l}} {{F_{{\rm{me}}}}(h,u,D) = }\\ {\quad \frac{{3{\mu _0}{m_{\rm{A}}}{m_{\rm{B}}}(2{\rm{cos}}{\kern 1pt} {\kern 1pt} \varphi - 5{{({\rm{sin}}{\kern 1pt} {\kern 1pt} \varphi )}^2}{\rm{cos}}{\kern 1pt} {\kern 1pt} \varphi )}}{{4\pi {{(h{{(t)}^2} + {{(D + u(t))}^2})}^2}}}} \end{array} $	（26）

弹性支撑的外部磁铁可以简化为一个单自由度振子，记质量、阻尼、刚度分别为m_e、c_e、k_e，运动方程如下：

$ {m_{\rm{e}}}\ddot u + {c_{\rm{e}}}\dot u + {k_{\rm{e}}}u = {F_{{\rm{me}}}}(h,u,D) $	（27）

从式(25)和式(26)可以看出，颤振能量收集器和外部磁铁之间通过磁力发生相互作用，它们受到的磁力均与二者的相对运动有关。

当外部磁铁固定时(u=0)，在系统的沉浮运动过程中，磁力的变化曲线如图 4所示。从图中可以看出，能量收集器受到的磁力在零点处为零，在两侧达到最大值后随距离增大而逐渐减小；外部磁铁的磁力在零点处最大，两侧逐渐减小。随着沉浮位移继续增大，内外磁铁部彼此远离，二者所受的磁力最终衰减至零，磁力与沉浮位移之间呈现强非线性关系。如图中所示，当2个磁铁的间距D减小时，磁力逐渐增强，磁力曲线与能量收集器的弹性恢复力曲线可以出现3个交点，且零点两侧的交点逐渐偏离零点。

由动力学控制方程式(18)可知，能量收集器的双稳态构型体现在沉浮自由度上。令能量收集器的速度和加速度分别为零，可以得到系统的平衡点应满足方程：

$ {K_{h1}}h + {K_{h3}}{h^3} = {F_{\rm{m}}}(h,D) $	（28）

联立求解代数超越方程式(25)和式(28)，可得颤振能量收集器的平衡点，对应图 4(a)中磁力曲线与弹性恢复力曲线交点的横坐标值。当磁铁间距较大时，磁力相对微弱，式(28)只有零解，能量收集器保持单稳态构型。

当间距小于临界值(D < D_cr)时，式(28)存在3个解{0，±h_s}，根据平衡点的稳定性理论可知，两侧的非零平衡点是稳定的，而原来无磁力系统的零平衡点由于磁力等效负刚度的存在变得不稳定，颤振能量收集器将出现双稳态构型。系统的总势能可表示为弹性势能和磁场势能的叠加，双稳态构型的势能曲线在零点两侧出现对称的势能阱，如图 5所示，在零风速状态下，颤振能量收集系统将稳定在其中一个势能阱内。

2.4 状态空间方程

为了方便进行系统的稳定性分析以及动响应求解，记[q₁  q₂  q₃  q₄  q₅]^T=[h  α  ḣ  ȧ  V]^T为系统的状态变量，将式(21)和式(22)、式(25)和式(26)代入系统的控制方程式(18)~式(20)，可化为状态向量的形式：

$ \begin{array}{l} \left[ {\begin{array}{*{20}{l}} {{q_1}}\\ {{{\dot q}_2}}\\ {{{\dot q}_3}}\\ {{{\dot q}_4}}\\ {{{\dot q}_5}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{q_3}}\\ {{q_4}}\\ { - {\mathit{\boldsymbol{M}}^{ - 1}}\mathit{\boldsymbol{K}}\left[ {\begin{array}{*{20}{l}} {{q_1}}\\ {{q_2}} \end{array}} \right] - {\mathit{\boldsymbol{M}}^{ - 1}}\mathit{\boldsymbol{C}}\left[ {\begin{array}{*{20}{l}} {{q_3}}\\ {{q_4}} \end{array}} \right] - {\mathit{\boldsymbol{M}}^{ - 1}}\left[ {\begin{array}{*{20}{l}} \theta \\ 0 \end{array}} \right]{q_5}}\\ {\frac{\theta }{{{C_{\rm{p}}}}}{q_3} - \frac{{2\theta }}{{{C_{\rm{p}}}{R_{\rm{L}}}}}{q_5}} \end{array}} \right] - \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\mathit{\boldsymbol{M}}^{ - 1}}\left[ {\begin{array}{*{20}{c}} {{K_{h3}}}\\ 0 \end{array}} \right]q_1^3 + {\mathit{\boldsymbol{M}}^{ - 1}}\left( {\left[ {\begin{array}{*{20}{c}} { - {Q_{\rm{L}}}}\\ {{Q_{\rm{M}}}} \end{array}} \right] + \left[ {\begin{array}{*{20}{c}} {{F_{\rm{m}}}}\\ 0 \end{array}} \right]} \right) \end{array} $

（29）

式中：M、K、C分别为沉浮和俯仰自由度的质量、刚度、阻尼矩阵。

定义颤振能量收集系统向负载电路输出的平均电功率为

$ P = \frac{1}{T}\int_0^T {\frac{{{V^2}}}{{{R_{\rm{L}}}}}} {\rm{d}}t $	（30）

式中:T为稳态极限环振荡的周期。

利用颤振能量收集系统的理论模型式(29)，采用四阶Runge-Kutta法对状态空间方程进行求解，计算能量收集器的结构响应以及输出电压、电功率。下文将对传统的无磁耦合颤振能量收集器(Flutter Energy Harvester, FEH)、磁耦合双稳态颤振能量收集器(BFEH)以及变势能阱双稳态颤振能量收集器(VBFEH)进行仿真对比研究，分析双稳态构型以及变势能阱双稳态技术对颤振能量收集性能的影响。

3 数值仿真验证 3.1 模型参数

仿真模型中的几何和物理参数如表 1所示，参考了本文作者在文献[14]中的模型参数，其中，悬臂梁材料为钢，压电片材料为PZT-5A。

3.2 无磁耦合颤振能量收集器(FEH)

当没有外部磁铁时，退化为传统的颤振能量收集器，将其作为能量收集性能的参照基准。根据压电-气动弹性系统的稳定性理论^{[7, 22]}，随着风速增大，当无磁耦合系统式(29)的特征根出现正实部时，能量收集器将发生动态失稳，此时的风速U就是能量收集系统工作的启动风速。

图 6所示为系统的特征根轨迹曲线，从图中可以得到无磁耦合系统的颤振风速为7.6 m/s，能量收集系统只有在高于该风速时才能开始工作，阻碍了其在低风速环境中的应用。

图 7所示为负载电阻对颤振能量收集器(FEH)的影响。其中，图 7(a)展示了负载电阻与颤振始发风速的关系，随着负载电阻的增大，能量收集器的颤振风速提高，达到最大值后逐渐减低。利用压电梁理论^[7]计算得到的本文能量收集器的机电耦合系数为θ=0.05 mN/V，是弱机电耦合系统，因此，负载电阻对颤振稳定性的改变并不明显。但从图 7(b)中可以看出，负载电阻对输出电功率的性能却有着显著的影响。系统的输出电压随负载增大而逐渐增大后趋于稳定值，而输出电功率的变化规律则与颤振风速相似，出现了先增大后降低的趋势，存在最优的负载电阻使得颤振风速和输出电功率几乎同时达到最大值，仿真分析得到的最优负载电阻约为R_{L, opt}=2 MΩ，下文将针对该最优负载进行双稳态系统能量收集性能的理论研究。

3.3 双稳态颤振能量收集性能(BFEH)

首先，对外部磁铁固定安装(u=0)的双稳态颤振能量收集器的性能进行分析。由式(28)可得，出现双稳态构型的磁铁临界间距为D_cr=24.6 mm，小于临界值时系统出现双稳态。

图 8所示为磁铁间距D=20 mm时，颤振能量收集器进入工作前后的动态响应时域对比。风速为U=4 m/s时，颤振能量收集器受扰动后的振动逐渐衰减，最终被束缚在其中一个势能阱内(h_s=-11.4 mm)。当风速为U=5.5 m/s时，系统在气流的激励下能够发生极限环振动，围绕2个平衡点在势能阱间往复跃迁。显然，双稳态颤振能量收集器在低于无磁耦合系统的线性颤振风速下，便能够有效地工作，用于收集电能。在U=5.5 m/s的风速下，双稳态颤振能量收集器的稳态输出电压的幅值为72 V，平均输出电功率达到1.3 mW。

3.3.1 内外部磁铁间距的影响

图 9为选取不同磁铁间距D时，双稳态颤振能量收集系统(BFEH)的输出电功率随风速变化的分叉图。从结果中可以看出，FEH达到线性颤振风速后才能发生极限环振动，输出电功率的幅值随着风速的增加而连续增长，这是气动弹性系统的超临界颤振特性。而BFEH在线性颤振风速之前，就可以利用双稳态势能阱间的跳跃行为(snap-through)，使系统提前进入极限环振动状态来收集电能。输出电功率的幅值在颤振发生前后存在着突变，气动弹性系统的动力学行为由超临界颤振转变为亚临界颤振。在文献[23]中，作者对磁耦合双稳态系统的亚临界颤振特性进行了风洞试验研究，亚临界颤振特性的拓扑改变与本节的仿真结果一致，体现了本文建模与仿真分析的正确性。

借助双稳态构型，颤振能量收集器能够在较低的风速环境下开始工作，意味着能量收集器的有效工作风速范围得到了拓宽。以磁铁间距D=23 mm为例，BFEH的颤振始发风速仅为3.2 m/s，比FEH降低了4.4 m/s，工作风速向低速区域延拓了58%。即使在FEH的线性颤振风速范围内(U >7.6 m/s)，BFEH的输出电功率也仍然相对较高，BFEH相比FEH在整个风速区间内都展现出优越的能量收集性能。

随着磁铁间距减小，双稳态能量收集器的输出电功率幅值有所提升，但同时也导致极限环振动的始发风速随之提高，系统的工作风速范围逐渐变窄。因此，拓宽能量收集系统工作的风速范围与增大输出电功率，这两方面的性能需要折中考虑，如图 9所示，无法通过调整磁铁间距实现性能的同步增强。

3.3.2 磁偶极矩的影响

从式(25)和式(26)可知，改变磁偶极矩矢量将直接影响非线性的磁力的大小，影响系统的动力学特性。以D=20 mm时的双稳态构型为例，分别取外部磁铁1倍、1.5倍、2倍的磁偶极矩进行仿真，结果如图 10所示，提高磁铁的磁偶极矩，能量收集系统受到的磁力增大，输出电能功率得到了有效地提升，但与减小磁铁间距的结果类似，增大磁偶极矩同样提高了颤振始发风速，缩小了能量收集系统工作的风速范围。

3.4 变势能阱技术的性能增强机理(VBFEH)

为了进一步增强颤振能量收集系统的性能，解决在双稳态颤振能量收集系统中，工作风速范围和输出电功率无法同步提升的局限，下文继续讨论弹簧支撑外部磁铁的运动对双稳态颤振能量收集性能的影响。

图 11所示为外部磁铁固定时(D=20 mm)的双稳态颤振能量收集(BFEH)，与弹簧支撑情况下(D=18.5 mm)的变势能阱双稳态颤振能量收集器(VBFEH)的性能对比。变势能阱双稳态能量收集器的参数见表 1。从结果中可以发现VBFEH的极限环始发风速比BFEH降低了1.6 m/s，能量收集系统的工作风速相对拓宽1/3，并且输出电功率没有像BFEH那样出现下降的趋势，反而在整个风速区间内得到了提升。以风速U=6 m/s时为例，BFEH的输出电功率为1.65 mW，而VBFEH的输出电功率达到2.49 mW，同比增加了约51%。

双稳态势能阱的深度可定义为不稳定平衡点和稳定平衡点之间的势能差值，它反映了系统初始动能为零时，从稳定点的势能阱底部上升至不稳定点处势能阱峰值处，继而完成势能阱间跃迁所需要的最少的非保守力做功。增强磁力可以增大平衡点间的距离，扩大双稳态的运动幅值，但也同时加深了势能阱的深度，束缚气动弹性系统在势阱间的跃迁，导致颤振风速升高。

图 12为VBFEH与BFEH的势能阱曲线对比，当外部磁铁固定时，系统的势阱深度为19.36 mJ；而对于外部磁铁有弹簧支撑的情况，当系统朝着沉浮位移零点运动时，外部磁铁受到排斥性磁力后逐渐远离内部磁铁，系统的势能曲线随着沉浮运动而实时改变，此时VBFEH的势能阱深度仅为2.45 mJ，能量收集器更容易在低风速下完成阱间跃迁运动。同时，磁场引起的磁势能在不稳定平衡点附近被释放，转换为颤振能量收集系统的动能，使系统发生绕2个平衡点的更大幅值的极限环振动。在相同风速下(以U=6 m/s为例)，VBFEH和BFEH的沉浮运动相图如图 13所示，VBFEH的相图更加平滑，在不稳定点处没有出现明显的凹陷，当系统穿越不稳定点时，VBFEH的沉浮运动速度较快，最终的极限环幅值也更大，从而提升了输出电压和电功率。

在随机激励的双稳态振动能量收集中，Lan和Qin^[19]通过附加吸引性磁铁来降低不稳定点处的势能，从而将势能阱的深度变浅，辅助双稳态系统在较低的振动水平下完成势阱间的跃迁。本文利用外部磁铁的弹性运动，以机械的方式实现了类似的势能阱调控，降低颤振能量收集器对环境风速的敏感度，研究结果证明，变势能阱技术能够在降低颤振始发风速的前提下，进一步提高能量收集的性能，为多稳态气动弹性能量收集系统的设计提供了参考。

4 结论

1) 引入非线性磁力和变势能阱技术，设计了一类新型变势能阱双稳态颤振能量收集器，建立了磁力-压电-气动弹性耦合的颤振能量收集系统理论分析模型。利用磁偶极模型和平衡点的稳定性理论，确定了颤振能量收集系统出现双稳态构型的内外部磁铁间距的参数范围。

2) 数值研究了双稳态构型的磁铁间距、磁偶极矩对能量收集性能的影响规律。在双稳态构型下，原始无磁力系统的超临界颤振变为亚临界颤振，低于线性颤振风速时系统便能够发生幅值稳定的极限环振动，在低风速下的能量收集性能显著提高。

3) 利用弹性支撑的外部磁铁的运动实现了变势能阱技术，解决了双稳态能量收集系统的工作风速范围与输出电功率无法同步提升的局限，从而进一步增强了双稳态颤振能量收集系统的性能。利用势能阱理论解释了变势能阱技术对颤振能量收集性能的增强机理，提高了系统对风场环境的适应性。