具有粘弹性的三稳态能量采集系统的动力学分析

Kinetics Analysis of a Tri-Stable Energy Harvesting System with Viscoelastic Properties

摘要: 本文提出了一种具有粘弹性项的三稳态能量采集系统，并深入研究了该系统的主共振响应特性。首先，利用多尺度法得到了主共振附近稳态振幅的解析解，并通过矩方法得到了系统的稳态响应矩。其次，通过蒙特卡罗仿真验证了所用方法的正确性。进一步通过数值分析研究随机跳变、随机分岔和有效带宽扩展等现象，最后分析了系统参数对振幅的二阶矩和平均输出功率的影响，发现对于每个时间延迟都存在一个最佳的反馈增益值使能量采集系统的平均输出功率最大，能量采集效果最佳。这些研究结果为优化三稳态能量采集系统的性能提供了理论支持。

Abstract: In this paper, a tri-stable energy harvesting system with viscoelastic term is proposed, and the main resonance response characteristics of the system are studied in depth. Firstly, the analytical solution of the steady-state amplitude near the principal resonance is obtained by using the multi-scale method, and the steady-state response moment of the system is obtained by the moment method. Secondly, the correctness of the proposed method is verified by Monte Carlo simulation. Further, the phenomena of random jumps, random bifurcations and effective bandwidth expansion are investigated by numerical analysis, and finally, the influence of system parameters on the second-order moment of the amplitude and the average output power is analyzed, and it is found that there is an optimal feedback gain value for each time delay, so that the average output power of the energy harvesting system is the largest and the energy harvesting effect is the best. These results provide theoretical support for optimizing the performance of tri-stable energy harvesting systems.

1. 引言

随着全球经济的快速发展和人口的持续增长，能源需求不断攀升，传统能源如煤炭、石油和天然气的过度开采与利用导致了严重的环境污染和生态破坏，加剧了全球气候变化和资源枯竭问题[1]。在这一背景下，寻找和开发清洁、可再生的新能源成为全球关注的焦点。能量采集技术作为一种新兴的能源解决方案，旨在从环境中采集分散且未被充分利用的能源(如机械振动能、热能、太阳能和风能等)，并将其转化为电能，为低功耗电子设备提供可持续的能源支持[2]。在众多能量采集技术中，基于机械振动能的采集系统因其能量来源广泛[3]、环境适应性强而备受关注。例如，机械设备运行中的振动[4]、人类行走或车辆行驶引发的环境震动[5]，甚至自然界的风能和波浪能，都可以通过能量采集技术转化为电能。

能量采集技术主要可分为电容式、压电式和电磁式三种类型[6]。其中，压电能量采集技术通过压电片的正压电效应将环境中的机械振动能转化为电能，为微电子设备提供电力支持。该技术能够为电池充电或直接为设备供电，从而有效延长无线传感器网络和物联网等设备的工作时间，并显著提升其运行可靠性。同时，压电能量采集系统因其结构简单、能量转换效率高、响应速度快等优势[7]，成为研究的热点之一。

压电能量采集技术的发展可追溯至1880年，法国著名物理学家Jacques Curie和Pierre Curie发现了压电效应[8]，为利用压电式转换机制实现振动能量采集奠定了理论基础。Sodano等[9]提出了一种线性压电悬臂梁结构并推导了机电耦合方程，通过分析与验证，讨论了系统重要参数与输出电压的影响。线性压电能量采集系统在特定情况下性能良好，但有效工作频率范围相对较窄，严重影响了其在实际应用中的普及性和实用性。为了克服这一局限，许多学者开始研究在能量采集器中加入非线性因素[10]-[13]。Gafforelli等[14]研究了单稳态能量采集系统的特性，通过实验发现非线性使系统的工作频带变宽。Erturk等[15]设计了一种双稳态能量采集器，并通过一系列实验研究了采集器在受到外部激励时的响应特性，研究结果表明，该模型在较宽的频率范围内能够有效地提取并转化大量能量，从而实现宽带能量的高效收集。然而，由于双稳态能量采集器在低频激励下产生电能的效率较低，国内外学者逐渐将研究焦点转向三稳态能量采集系统。其中，Zhou等[16]-[18]设计了一种三稳态压电能量采集器，其结构包含悬臂梁，悬臂梁的底端固定一尖端磁体，这一尖端磁体在两个固定磁体之间振荡，非常类似于双稳态采集器的早期设计。Panyam等[19]揭示了三稳态能量采集系统的非线性响应机制，结果表明系统具有多解特性，且能够实现高能阱间振荡。

在这些能量采集器的基础上，学者们开始研究各种因素对其的影响。Yang等[20]研究了一种在窄带随机激励下具有时间延迟的新型粘弹性能量收集器，并通过理论分析得到系统的稳态响应。与粘性物体或弹性物体不同，粘弹性介质表现出对应变速率敏感的粘性行为，又表现出类似弹簧的弹性性质。这种“中间”特征似乎比粘性和弹性这两种理想状态更适合于模拟实际问题。

在旋转机械领域，涡轮叶片受周期性气流扰动产生的气动激励，其主导频率成分集中在转子通过频率及其谐波附近；土木工程中斜拉桥索结构的风致振动，当风速处于锁定区间时，其涡脱频率将稳定在结构固有频率邻域；甚至生物医学工程领域，人工心脏瓣膜启闭过程承受的血流脉动压力，其频谱能量在心动周期特征频率处呈现显著聚集。这些实例印证了：当系统响应受特定共振机制主导时，外界激励的有效成分往往自动窄带化。本文采用的窄带激励假设，正是捕捉了此类物理本质特征，为揭示系统在特征频段内的非线性动力学行为提供了有效的分析框架。

本文的其余部分组织如下：第2节是三稳态能量采集系统的模型建立以及理论分析；第3节通过数值模拟分析了能量采集系统各参数对其采集性能的影响。

2. 三稳态能量采集系统的模型建立以及理论分析

2.1. 用多尺度方法求解系统的稳态响应

Figure 1. Model of a tri-stable energy harvesting system: (a) Equivalent schematic diagram of the collector; (b) Piezoelectric mechanism with time delay

图1. 三稳态能量采集系统的模型：(a) 采集器的等效原理图；(b) 带有时间延迟的压电机制

考虑一个粘弹性的三稳态能量采集系统如图1所示，该模型由一个与电路耦合的压电器件及机械振荡器组成。通过牛顿第二定律与基尔霍夫律得到系统的运动方程为：

$\{\begin{array}{l}M\ddot{\overline{X}}+\overline{\rho }Z\left(\delta ,t,\overline{X}\right)+\overline{C}\dot{\overline{X}}\left(t\right)+\text{d}\overline{U}\left(\overline{X}\right)/\text{d}\overline{X}-\chi \overline{V}\left(t\right)=-M{\ddot{\overline{X}}}_b\\ C_p\dot{\overline{V}}\left(t\right)+V\left(t\right)/R+\overline{k}\dot{\overline{X}}\left(t\right)=\overline{\beta }\overline{V}\left(t-\tau \right)\end{array}$ (1)

其中，$X$ 表示磁铁$M$ 的位移，$\dot{\overline{X}}$ 为磁铁$M$ 的速度，$\overline{\rho }$ 表示粘弹性项，$X_b$ 表示基本加速度，$\overline{C}$ 表示线性粘性阻尼系数，势函数为$U\left(\overline{X}\right)=k_1{\overline{X}}^2/2+k_2{\overline{X}}^4/4+k_3{\overline{X}}^6/6$ ，$\chi$ 表示线性机电耦合系数，$C_p$ 表示压电元件的有效电容，$\overline{V}$ 表示感应电压，$R$ 表示等效电阻，即$R=R_l{R}_p\left(R_l+R_p\right)$ ，其中$R_l$ 为负载电阻，$R_p$ 为压电电阻。$\overline{k}$ 为电路中的压电耦合项，$\beta$ 和$\tau$ 分别为电路中的反馈增益和时间延迟。

为了使式(1)无量纲化，引入以下参数变换，

$\begin{array}{l}x=\frac{X}{l},x_b=\frac{{\overline{X}}_b}{l},V=\frac{C_p}{\overline{\chi }l},{\omega }_n=\sqrt{\frac{k_1}{M}},C=\frac{\overline{C}}{M},{\widehat{\alpha }}_1=\frac{k_2{l}^2}{M},\\ {\widehat{\alpha }}_2=\frac{k_3{l}^2}{M},\stackrel{⌢}{\rho }=\frac{\overline{\rho }}{M},\lambda =\frac{1}{R{C}_p},\widehat{\chi }=\frac{{\chi }^{2}l}{M{C}_p},k=\frac{k}{\chi },\beta =\frac{\overline{\beta }}{C_p}\end{array}$

由此，具有粘弹性的无量纲方程可以表示为：

$\{\begin{array}{l}\ddot{x}\left(t\right)+\stackrel{⌢}{\rho }Z\left(\beta ,t,x\right)+Cx\left(t\right)+\text{d}U\left(x\right)/\text{d}x-\widehat{\chi }v\left(t\right)=-{\ddot{x}}_b,\\ \dot{v}\left(t\right)+\lambda v\left(t\right)+k\dot{x}\left(t\right)=\beta v\left(t-\tau \right)\end{array}$ (2)

式中，$x\left(t\right)$ 为质量$M$ 的无量纲位移，$C=2\mu {\omega }_n$ 为无量纲线性粘性阻尼系数，$v\left(t\right)$ 为通过负载电阻的无量纲电压，$U\left(x\right)={\omega }_n^2{x}^2/2+{\alpha }_1^4/4+{\alpha }_2^6/6$ 为无量纲势函数，${\omega }_n$ 为系统固有频率，${\alpha }_1,{\alpha }_2$ 为无量纲非线性刚度系数，$\stackrel{⌢}{\chi }$ 为无量纲机电耦合系数，$\lambda$ 为无量纲时间常数比，k为电路中的无量纲压电耦合项，$\beta$ 和$\tau$ 分别为电路中的无量纲反馈增益和时间延迟。

为了应用多尺度方法，本节引入小参数$\epsilon$ 来得到无量纲参数：$\stackrel{⌢}{\rho }=\epsilon \rho ,\widehat{\mu }=\epsilon \mu ,{\widehat{\alpha }}_1=\epsilon {\alpha }_1,{\widehat{\alpha }}_2=\epsilon {\alpha }_2,\widehat{\chi }=\epsilon \chi ,$ $\widehat{F}=\epsilon F$ 。式(2)可改写为：

$\{\begin{array}{l}\ddot{x}\left(t\right)+\epsilon \rho Z\left(\beta ,x,t\right)+2\mu {\omega }_n\dot{x}\left(t\right)+{\omega }_n^{2}x\left(t\right)+{\alpha }_1{x}^3\left(t\right)+{\alpha }_2{x}^3\left(t\right)-\chi v\left(t\right)=F\cos \left(\Omega t+\gamma W\left(t\right)\right),\\ \dot{v}\left(t\right)+\lambda v\left(t\right)+\kappa x\left(t\right)=\beta v\left(t-\tau \right).\end{array}$ (3)

根据文献[21] [22]，粘弹性项可以表示为：

$Z\left(\delta ,x,t\right)={\displaystyle {\int }_0^t\frac{1}{\delta }}\exp \left(-\frac{t-s}{\beta }\right)x\left(s\right)\text{d}s$ (4)

根据多尺度法，$x\left(t\right)$ 和$v\left(t\right)$ 的近似解可设为：

$\{\begin{array}{l}x\left(t\right)=x_0\left(T_0,T_1\right)+\epsilon x_1\left(T_0,T_1\right)+O\left({\epsilon }^2\right),\\ v\left(t\right)=v_0\left(T_0,T_1\right)+\epsilon v_1\left(T_0,T_1\right)+O\left({\epsilon }^2\right).\end{array}$ (5)

引入微分算子

$\{\begin{array}{l}\frac{\text{d}}{\text{d}t}=D_0+\epsilon D_1+O\left({\epsilon }^2\right),\\ \frac{{\text{d}}^2}{\text{d}{t}^2}=D_0^2+2\epsilon D_0{D}_1+O\left({\epsilon }^2\right),\end{array}$ (6)

其中，$D_0=\partial /\partial T_0,D_1=\partial /\partial T_1$ 。

把式(5)、(6)代入式(3)中，取$\epsilon$ 的相似幂次项，并分别令它们的系数为零得到如下的方程：

$\begin{array}{l}O\left({\epsilon }^0\right):\\ \{\begin{array}{l}{D}_0^2{x}_0\left(T_0,T_1\right)+{\omega }_n^2{x}_0\left(T_0,T_1\right)=0,\\ D_0{V}_0\left(T_0,T_1\right)+\lambda V_0\left(T_0,T_1\right)+k{D}_0{x}_0\left(T_0,T_1\right)=\beta v_0\left(T_0-\tau ,T_1\right).\end{array}\end{array}$ (7)

$\begin{array}{l}O\left({\epsilon }^1\right):\\ \{\begin{array}{l}{D}_0^2{x}_1\left(T_0,T_1\right)+{\omega }_n^2{x}_1\left(T_0,T_1\right)=-2D_0{D}_1{x}_0\left(T_0,T_1\right)-2\mu {\omega }_n{D}_0{x}_0\left(T_0,T_1\right)-{\alpha }_2{x}_0^5\left(T_0,T_1\right)-{\alpha }_1{x}_0^3\left(T_0,T_1\right)\\ +\chi v_0\left(T_0,T_1\right)+F\cos \left(\Omega T_0+\gamma W\left(T_1\right)\right)-\rho Z\left(\beta ,x,t\right)\\ D_0{v}_1\left(T_0,T_1\right)+\lambda v_1\left(T_0,T_1\right)=-D_1{v}_0\left(T_0,T_1\right)-k{D}_0{x}_0\left(T_0,T_1\right)+\beta v_1\left(T_0-\tau ,T_1\right).\end{array}\end{array}$ (8)

方程(3)的通解可表示为：

$\{\begin{array}{l}{x}_0\left(T_0,T_1\right)=A_1\left(T_1\right){\text{e}}^{i\omega T_0}+cc,\\ v_0\left(T_0,T_1\right)=B_1\left(T_1\right){\text{e}}^{i\omega T_0}+cc.\end{array}$ (9)

式中，$A_1\left(T_1\right),B\left(T_1\right)$ 为复函数，$cc$ 表示前面各项的共轭，将式(9)代入式(7)可得：

$i{\omega }_n{B}_1+\lambda B_1+i{\omega }_{n}k{A}_1=\beta B_1{\text{e}}^{i{\omega }_n\tau }.$ (10)

由式(10)得到$B_1\left(T_1\right)$ 为：

$B_1\left(T_1\right)=\frac{-i{\omega }_{n}k}{i{\omega }_n+\lambda -\beta {\text{e}}^{-i{\omega }_n\tau }}A\left(T_1\right).$ (11)

将式(9)代入式(8)，得到：

$\begin{array}{c}{D}_0^2{x}_1+{\omega }_n^2{x}_1=\left(-2i{\omega }_n{A^{\prime }}_1-2i{\omega }_n^2\mu A_1-3{\alpha }_1\overline{A_1}{A}_1^2-10{\alpha }_2{\overline{A_1}}^2{A}_1^3\right){\text{e}}^{i{\omega }_n{T}_0}-\frac{i\omega k\chi A_1}{i\omega +\lambda -\beta {\text{e}}^{i{\omega }_n\tau }}{\text{e}}^{i{\omega }_n{T}_0}\\ -\left({\alpha }_1{A}_1^3+5a_2{A}_1^4\overline{A_1}\right){\text{e}}^{i{\omega }_n{T}_0}-\alpha A_1^5{\text{e}}^{5i{\omega }_n{T}_0}+\frac{F}{2}{\text{e}}^{i\left(\Omega T_0+\gamma W\left(T_1\right)\right)}-\rho Z\left(\delta ,x_0,T_0\right)+cc.\end{array}$ (12)

根据式(4)和式(9)得到粘弹性项可以表示成如下形式：

$\begin{array}{l}Z\left(\delta ,x_0,T_0\right)=\frac{1}{\delta }{\displaystyle {\int }_0^{T_0}\exp \left(-\frac{T_0-s}{\delta }\right)x_0\left(s\right)\text{d}s}\\ =\frac{A_1}{\delta }\exp \left(-\frac{T_0}{\delta }\right){\displaystyle {\int }_0^{T_0}\exp \left(\left(\frac{1}{\delta }+i{\omega }_n\right)s\right)\text{d}s}+cc\\ =\frac{A_1}{1+i{\omega }_n}\exp \left(-\frac{T_0}{\delta }\right)\left(\exp \left(\left(\frac{1}{\delta }+i{\omega }_n\right)T_0\right)-1\right)+cc\\ =\frac{1-{\omega }_n\delta i}{1+{\omega }_n{\delta }^2}{A}_1\exp \left(-\frac{T_0}{\delta }\right)+cc\end{array}$ (13)

消除久期项得到：

$-2i{\omega }_n{A^{\prime }}_1-\rho \frac{1-\omega \delta i}{1+{\omega }^2{\delta }^2}{A}_1-2i{\omega }_n^2\mu A_1-3{\alpha }_1\overline{A_1}{A}_1^2-10{\alpha }_2{\overline{A_1}}^2{A}_1^3-\frac{i{\omega }_{n}k\chi A_1}{i{\omega }_n+\lambda -\beta {\text{e}}^{-i{\omega }_n\tau }}+\frac{F}{2}{\text{e}}^{i\left(\Omega T_0-{\omega }_n{T}_0+\gamma W\left(T_1\right)\right)}=0.$ (14)

为了研究式(3)的主共振，可以引入失谐频率$\sigma$ 来表示主共振为：

$\Omega ={\omega }_n+\epsilon \sigma .$ (15)

将式(15)代入(14)中可得：

$-2i{\omega }_n{A^{\prime }}_1-A_1\rho \frac{1-\omega \delta i}{1+{\omega }^2{\delta }^{{}^2}}-2i{\omega }_n^2\mu A_1-3{\alpha }_1\overline{A_1}{A}_1^2-10{\alpha }_2{\overline{A_1}}^2{A}_1^3-\frac{i{\omega }_{n}k\chi A_1}{i{\omega }_n+\lambda -\beta {\text{e}}^{-i{\omega }_n\tau }}+\frac{F}{2}{\text{e}}^{i\left(\sigma T_1+\gamma W\left(T_1\right)\right)}=0$ (16)

由于$A_1$ 是复数，可设$A_1\left(T_1\right)=\left(1/2\right)a\left(T_1\right){\text{e}}^{i\phi \left(T_1\right)}$ ，同时引入$A_1\left(T_1\right)=\left(1/2\right)a\left(T_1\right){\text{e}}^{i\phi \left(T_1\right)}$ 。$a$ 代表系统振动的振幅。将二者都代入式(16)，再将实部与虚部分离得到：

$\{\begin{array}{l}{a}^{\prime }=-{\mu }_{1}a+\frac{\rho \delta }{2\left(1+{\omega }_n^2\delta \right)}a+\frac{F}{2{\omega }_n}\sin \eta ,\\ a{\eta }^{\prime }={\sigma }_{1}a-\frac{\rho }{2{\omega }_n\left(1+{\omega }_n^2\delta \right)}a-\frac{3{\alpha }_1{a}^3}{8{\omega }_n}-\frac{5{\alpha }_2{a}^5}{16{\omega }_n}+\frac{F}{2{\omega }_n}\cos \eta +\gamma a{W}^{\prime }\left(T_1\right).\end{array}$ $$ (17)

其中，$\eta \left(T_1\right)=\sigma T_1+\gamma W\left(T_1\right)-\phi \left(T_1\right),{\mu }_1=\mu {\omega }_n+\chi k\left(\lambda -\beta \cos \left({\omega }_n\tau \right)\right)/H,{\sigma }_1=\sigma -\chi k\left({\omega }_n+\beta \sin \left(\omega \tau \right)\right)/H,$

$H=2\left[{\left({\omega }_n+\beta \sin \left({\omega }_n\tau \right)\right)}^2+{\left(\lambda -\beta \cos \left({\omega }_n\tau \right)\right)}^2\right].$

$$ 由式(17)可以得到$a$ 和$\eta$ ，同时推导出位移和电压的一阶近似解如下：

$x\left(t\right)=a\left(\epsilon t\right)\cos \left(\Omega t-\eta \left(\epsilon t\right)\right)+{\rm O}\left(\epsilon \right),$ (18)

$v\left(t\right)=\frac{-i{\omega }_{n}k}{i{\omega }_n+\lambda -\beta {\text{e}}^{i{\omega }_n\tau }}a\left(\epsilon t\right)\cos \left(\Omega t-\eta \left(\epsilon t\right)+\theta \right)+{\rm O}\left(\epsilon \right),$ (19)

其中，$\theta =\arctan \left[\left(\lambda -\beta \cos \left({\omega }_n\tau \right)\right)/\left({\omega }_n+\beta \sin \left({\omega }_n\tau \right)\right)\right]$ 。

2.2. 用矩方法求解系统的稳态响应矩

由式(16)可知，稳态解的随机微分方程为：

$\{\begin{array}{l}\text{d}a=\left[-{\mu }_{1}a+\frac{\rho \delta }{2\left(1+{\omega }_n^2{\delta }^2\right)}+\frac{F}{2{\omega }_n}\sin \eta \right]\text{d}{T}_1,\\ \text{d}\eta =\left[{\sigma }_1-\frac{\rho \delta }{2{\omega }_n\left(1+{\omega }_n^2{\delta }^2\right)}-\frac{3{\alpha }_1{a}^2}{8{\omega }_n}-\frac{5{\alpha }_2{a}^4}{16{\omega }_n}+\frac{F}{2{\omega }_{n}a}\cos \eta \right]\text{d}{T}_1+\gamma \text{d}W\left(T_1\right)\text{.}\end{array}$ (20)

当$\gamma =0$ 时，式(20)不存在随机激励，则设$a^{\prime }=0,{\eta }^{\prime }=0$ 可得到稳态解$a=a_0,\eta ={\eta }_0$ 满足如下方程：

$\{\begin{array}{l}{\mu }_1{a}_0-\frac{\rho \delta a_0}{2\left(1+{\omega }_n^2{\delta }^2\right)}=\frac{F}{2{\omega }_n}\sin {\eta }_0,\\ {\sigma }_1{a}_0-\frac{\rho \delta a_0}{2{\omega }_n\left(1+{\omega }_n^2{\delta }^2\right)}-\frac{3a_1{a}^3}{8{\omega }_n}-\frac{5a_2{a}^5}{16{\omega }_n}=-\frac{F}{2{\omega }_n}\cos {\eta }_0.\end{array}$ (21)

简化式(21)，可以得到：

${\left({\sigma }_1{a}_0-\frac{\rho \delta a_0}{2{\omega }_n\left(1+{\omega }_n^2{\delta }^2\right)}-\frac{3{\alpha }_1{a}_0^3}{8{\omega }_n}-\frac{5{\alpha }_2{a}_0^4}{16{\omega }_n}\right)}^2+{\left({\mu }_1{a}_0-\frac{\rho \delta a_0}{2\left(1+{\omega }^2{\delta }^2\right)}\right)}^2=\frac{F^2}{4{\omega }_n^2}.$ (22)

当$\gamma \ne 0$ 时，式(20)的非平凡解设为：$a=a_0+a_1,\eta ={\eta }_0+{\eta }_1$ ，其中$a_0,{\eta }_0$ 满足(20)，$a_1,{\eta }_1$ 是微扰项，将式(22)代入式(20)，得到式(20)在$\left(a_0,{\eta }_0\right)$ 处的线性化方程：

$\{\begin{array}{l}{a^{\prime }}_1=\left(-{\mu }_1+\frac{\rho \delta }{2\left(1+{\omega }_n^2{\beta }^2\right)}\right)a_1+\frac{F}{2{\omega }_n}\cos {\eta }_0\cdot {\eta }_1,\\ {{\eta }^{\prime }}_1=\left(\frac{{\sigma }_1}{a_0}-\frac{\rho \delta }{2\left(1+{\omega }_n^2{\beta }^2\right)}-\frac{9\alpha a_0}{8{\omega }_n}-\frac{25{\alpha }_1{a}_0^3}{16{\omega }_n}\right)a_1-\frac{F}{2a_0{\omega }_n}\sin {\eta }_0\cdot {\eta }_1+\gamma W^{\prime }\left(T_1\right).\end{array}$ (23)

根据式(21)，方程(23)可改写为：

$\{\begin{array}{l}{a^{\prime }}_1=-H_1{a}_0-G_1{\eta }_1,\\ {{\eta }^{\prime }}_1=G_2{a}_1-H_1{\eta }_1+\gamma W^{\prime }\left(T_1\right).\end{array}$ (24)

其中，$H_1={\mu }_1-\frac{\rho \delta }{2\left(1+{\omega }_n^2{\delta }^2\right)}$ ，$G_1={\sigma }_1{a}_0-\frac{\rho \delta }{2\left(1+{\omega }_n^w{\delta }^2\right)}-3{\alpha }_1{a}_0^3/8{\omega }_n-5{\alpha }_2{a}_0^5/16{\omega }_n$ ，$G_2={\sigma }_1/a_0-\frac{\rho \delta }{2{\omega }_n{a}_0\left(1+{\omega }_n^2{\delta }^2\right)}-9{\alpha }_1{a}_0/8{\omega }_n-{\alpha }_2{a}_0^{3}75/16{\omega }_n$ 。

式(23)的二阶矩满足下列等式：

$\{\begin{array}{l}\frac{\text{d}E{a}^2}{\text{d}{T}_1}=-2H_{1}E{a}_1^2-2G_{1}E{a}_1{\eta }_1,\\ \frac{\text{d}E{a}_1{\eta }_1}{\text{d}{T}_1}=G_{2}E{a}_1^2-2H_{1}E{a}_1{\eta }_1-2\mu E{\eta }_1^2,\\ \frac{\text{d}E{\eta }_1^2}{\text{d}{T}_1}=2G_{2}E{a}_1{\eta }_1-2H_{1}E{\eta }_1^2+{\gamma }^2.\end{array}$ (25)

令$\frac{\text{d}E{a}_1^2}{\text{d}{T}_1}=\frac{\text{d}E{a}_1{\eta }_1}{\text{d}{T}_1}=\frac{\text{d}E{\eta }_1^2}{\text{d}{T}_1}=0$ ，可以得到二阶稳态矩：

$\{\begin{array}{l}E{a}_1^2=\frac{G_1^2{\gamma }^2}{4H_1\left(H_1^2+G_1{G}_2\right)},\\ E{a}_1{\eta }_1=\frac{-G_1^2{\gamma }^2}{4H_1\left(H_1^2+G_1{G}_2\right)},\\ E{\eta }_1^2=\frac{\left(2H_1^2+G_1{G}_2\right){\gamma }^2}{4H_1\left(H_1^2+G_1{G}_2\right)}.\end{array}$ (26)

取方程(22)两边的期望，可得到解的一阶和二阶稳态矩为：

$\{\begin{array}{l}Ea=a_0,E{a}^2=a_0+E{a}_1^2,\\ E\eta ={\eta }_0,E{\eta }^2={\eta }_0^2+E{\eta }_1^2.\end{array}$ (27)

同理，输出电压的一阶和二阶稳态矩为：

$Ev=\frac{{\omega }_{n}k}{\sqrt{H/2}}{a}_0,$ (28)

$E{v}^2=\frac{2{\omega }_n^2{k}^2}{H}\left(a_0^2+E{a}_1^2\right),$ (29)

平均输出功率为：

$EP=\lambda \chi E{v}^2.$ (30)

3. 数值分析

3.1. 随机跳变现象

为了确保理论解与数值解之间具有良好的一致性，本节通过蒙特卡罗数值模拟直观地展示了系统的响应特性。图2展示了系统的数值模拟结果。其中，通过多尺度方法得到的振幅一阶和二阶稳态矩的理论结果(式(27)和(28))以实线形式呈现。同时，采用蒙特卡罗方法对原始系统(4)进行数值求解，其结果以圆圈标记。从图中可以看出，理论结果与数值模拟结果高度一致，表明理论分析的准确性。

图3展示了振幅的一阶矩和输出电压的一阶矩随激励频率变化的关系曲线。从图中可以观察到，这两条曲线均呈现出向右弯曲的共振特性。随着系统非线性刚度系数的增加，曲线向右弯曲的趋势变得更加显著。这一现象表明，通过引入非线性刚度，可以有效扩展三稳态能量采集器的工作带宽，从而使其在更宽的频率范围内实现高效的能量转换。因此，非线性刚度的合理设计对于提升三稳态能量采集器的性能具有重要意义。

Figure 2. Diagram of the first-order steady-state moment of amplitude and voltage as a function of detuning parameter σ: (a) First-order steady-state moment of amplitude; (b) First-order steady-state moment of the voltage. Solid line: Theoretical solution; Circle: Numerical solution

图2. 振幅与电压的一阶稳态矩随失谐参数σ变化的函数图：(a) 振幅的一阶稳态矩；(b) 电压的一阶稳态矩。实线：理论解；圆圈：数值解

Figure 3. Diagram of the first-order steady-state moment of amplitude and voltage as a function of detuning parameters: (a) First-order steady-state moment of amplitude; (b) First-order steady-state moment of the voltage

图3. 振幅与电压的一阶稳态矩随失谐参数变化的函数图：(a) 振幅的一阶稳态矩；(b) 电压的一阶稳态矩

3.2. 系统参数对能量采集系统的影响

在接下来的分析中，依次引入了噪声强度、反馈增益、时滞、粘弹性项、压电耦合项和时间常数比，研究了这些因素对振幅的二阶稳态矩$E{a}^2$ 和平均输出功$EP$ 率的影响。

图4描述了振幅的二阶稳态矩$E{a}^2$ 和平均功率$EP$ 与反馈增益$\beta$ 和时间延迟$\tau$ 之间的关系。结果显示，随着反馈增益$\beta$ 的不断增加，振幅的二阶稳态矩$E{a}^2$ 首先下降，达到一个最小值后整体呈上升趋势。相比之下，平均功率$EP$ 则表现出不同的变化趋势：随着反馈增益$\beta$ 的增加，平均功率$EP$ 首先上升，达到一个最大值后整体呈下降趋势。因此，可以推测对于特定的反馈增益$\beta$ ，存在最优的时间延迟$\tau$ ，使得能量采集器的性能达到最佳。因此，能够推测出，针对特定的反馈增益$\beta$ ，存在一个最优时间延迟$\tau$ ，可让能量采集器性能达到最佳状态。在实际系统场景中，我们能够通过探索，找到一组最优的反馈增益与时间延迟组合，借此让能量采集系统的工作效能，在所处环境与条件限制下，达到近乎理想的水平，实现能量采集最大化，或者使采集效率和稳定性等关键性能指标得到显著提升。

Figure 4. Three-dimensional plot of the second-order steady-state moment of the amplitude as a function of feedback gain and time delay

图4. 振幅的二阶稳态矩与随反馈增益和时间延迟变化的三维图

图5展示了噪声强度对振幅的二阶稳态矩和平均功率的影响。分析结果表明，随着噪声强度$\gamma$ 的增加，振幅的二阶稳态矩$E{a}^2$ 的最小值和平均功率$EP$ 的最大值均呈现增大的趋势。这表明，图中出现的极限环现象是由噪声强度的增大所诱导的。此外，当反馈增益$\beta <-0.6$ 时，平均功率先持续增大，随后迅速下降。这一现象揭示了反馈增益对能量采集器性能的显著影响。

Figure 5. Diagram of the second-order steady-state moment of the amplitude and the average output power as a function of feedback gain at different noise intensities

图5. 不同噪声强度下振幅的二阶稳态矩与平均输出功率随反馈增益的变化图

综合分析结果表明，当反馈增益约为0.6时，能量采集器的性能达到最佳状态。此外，通过对比不同噪声强度下的平均功率，我们发现较大的随机噪声水平可以在一定程度上提高能量采集器的性能。这可能是因为适当的噪声扰动能够激发系统的非线性特性，从而增强其能量转换效率。因此，在实际应用中，合理控制噪声强度对于优化能量采集器的性能至关重要。

4. 结束语

本文研究了窄带随机激励下具有粘弹性项的三稳态能量采集系统的响应特性。通过多尺度方法，提出了一种求解系统解析解的有效方法。首先，利用该方法推导了系统一阶矩和二阶矩的解析表达式，为系统的统计特性分析提供了理论基础。其次，对所建立的模型进行验证，结果表明理论解析解与数值模拟结果高度一致，验证了模型的准确性和可靠性。进一步分析发现，在参数空间中，系统存在随机跳跃现象，这些现象揭示了系统在特定参数范围内的复杂非线性动力学行为。最后，讨论了系统参数对均方电压和平均输出功率的影响。结果表明，随着压电耦合系数的增大，平均输出功率不断增大。对于每个时间延迟都存在一个最佳的反馈增益值使能量采集系统的平均输出功率最大，能量采集效果最佳。同时，对于较大的随机噪声水平可以在一定程度上提高能量采集器的性能。