基于等效线性化方法的非线性振动能量采集器功率分析

李佳诚， 王志霞, 2， 王 炜, 2， 王 辰

(1. 天津大学 机械工程学院，天津 300350； 2. 天津市非线性动力学与控制重点实验室， 天津 300350； 3. 香港理工大学 土木与环境工程系， 香港 999077)

近年来，以便携式通信系统、无线传感器网络为代表的低功耗电子技术取得了长足的发展，在人们的日常生产活动中发挥着日益重要的作用。与此同时，考虑到很多无线传感器都工作在传统电力网络不能直接提供电能的区域，而采用电池供电又无法实现“永久化”的长效运行模式，进而成为制约这些设备进一步广泛应用的瓶颈。因此，如何有效地从环境中获取能量，摆脱传统供给方式的束缚，为众多无线传感器和低功耗电子设备提供可靠、安全且免维护的电力来源，就成为科研人员关注的焦点问题。

作为环境能量的主要形式之一，振动能不受温度、尺度等条件的制约，具备成为可靠能量来源的基本要素。目前的振动能量采集器按照工作原理不同，可分为电容式、压电式和电磁式三种；其中电磁式振动能量采集器(electromagnetic vibration energy harvester, EMH)具有低频性能好、发电量大、无需驱动电源等特点，特别适用于低频环境下(小于150 Hz)能量采集。

在工作原理方面，EMH主要是利用法拉第电磁感应定律，结构中包括永磁体和感应线圈绕组。在外界激励作用下，永磁体与线圈之间产生相对运动，导致线圈中的磁通量发生变化，从而产生感生电动势。然而，在实际应用过程中，研究人员却发现，为确保EMH实现最佳的能量转化效果，一方面需要使其尽可能地工作在振子的谐振频率附近，另一方面还需要对结构参数进行优化，只有使振动、感生、储能各环节相互协调，才能达到最佳的能量采集效果。因此，如何寻求有效手段来提高采集器的频带宽度，同时实施全面的参数优化工作，成为了EMH设计、应用过程中亟待解决的核心课题。

具体而言，传统的EMH以线性振子为主，如：Williams等[1]就利用薄膜振动的原理，设计了最初的EMH模型，还基于线性振子的思想结合有限元方法对结构进行了优化；王佩红[2]应用Maxwell较早开展了EMH的设计优化工作，研究了永磁体及各种结构参数与最大输出电压和输出功率之间的关系，从而为实际器件的加工制作提供了可靠的参考依据；李志宏等[3]设计了一种复合式能量采集模型，利用传递函数推导出了压电式和电磁式复合采集器的功率表达式，并研究了负载参数变化对于功率的影响。

虽然线性振子比较易于实现系统的功率分析、确定最优参数的取值范围，但是由于其工作带宽有限，当外界激励远离谐振频率时，振幅显著减弱，并伴随着输出电压陡然下降，难以起到有效的能量采集效果。因此，需要采用非线性手段达到增加带宽、提高振幅、改善输出效果的目的。在非线性能量采集器方面，武丽森等[4]提出了一种新型的非线性压电-电磁复合式结构，其基本原理在振动磁铁的上下各固定一个磁铁，利用磁铁相互作用的非线性力来改善俘能效果，试验结果显示当在磁铁间距为2.5 mm时，电磁单元最佳负载为18 Ω，3 dB带宽为15πrad/s；代显智等[5]利用等效磁荷理论，计算了结构振动时受到的磁力，并采用数值手段对其进行参数拟合，由此研究了采集器的非线性振动特性，并且讨论了定量分析的误差产生原因。

虽然在采集器中引入非线性因素可以达到增加带宽、改善输出效果的目的[6-16]，但随着非线性手段的引入，系统的控制方程也变得更为复杂，这无疑增大了功率分析、参数优化的难度。另一方面，学者们提出了获取最大电流或电压两种不统一的功率优化方案。如：Yan等[17]建立了双稳态多永磁电磁能量采集器理论模型，并用谐波平衡法求出解析解，再通过数值模拟和试验研究了激励频率、激励幅度、势阱电阻和势阱形状对最大电压的影响，但所采用的谐波平衡法的精度不够高，造成理论分析结果和数值模拟结果误差较大。Kecik[18]研究了非线性伪磁悬浮能量采集系统，利用Matlab自带ODE程序求解非线性方程并详细介绍了伪磁悬浮物理参数(振幅激励、负载电阻和耦合系数)对位移响应和最大电流的影响，但由于非线性因素的存在，Kecik只选取几个特定的参数值进行研究。

为深入研究上述复杂非线性因素对采集器振动特性的影响，进而形成更为便捷，精度更高的功率优化方案，本文提出了一种基于等效线性化思想[19]的非线性EMH稳态响应计算、采集器功率分析的方法。首先，将动态频率法[20]拓展至1.5自由度强非线性复杂振动系统，考虑复杂非线性因素、周期激励作用下，计算出系统的动态频率与稳态响应。文献[20]表明，较之传统的谐波平衡法，动态频率方法可以更为高效、精确地处理强非线性复杂振动问题。其次，在系统的能量方程中，利用动态频率进行Fourier级数展开，将第一阶谐波成分之外的部分视为与其等效线性化系统相对应的高频外激励，借助谐波平衡的观点，确定相关谐波项的系数；该等效线性化方程的解则对应于原有非线性系统的解，与基于弱非线性系统平均法获得的线性化系统相比，其求解精度更高。最后，利用传递函数开展针对EMH等效线性化系统的输出功率分析，讨论系统负载等关键参数变化对于最大功率的影响，从根本上克服了非线性因素存在对于优化策略的负面影响。实现了引入非线性成分、有效扩大带宽，与降低参数优化复杂性、充分发挥结构设计优势之间的平衡。

1 EMH的数学模型

如图1所示，采用悬臂梁结构EMH，整体上由振动环节和电路环节两个部分组成[21-23]，其中：振动环节采集环境中的机械能，电路环节将机械能转化为电能并加以储存，以达到向外界输出电能的功能。以下分别从上述两方面建立EMH的通用化数学模型。

1.1 振动环节的数学模型

如图2所示，EMH的拾振部分可以简化为一个机械弹簧-质量块-阻尼器组合而成的物理模型，此处的非线性因素包括磁力非线性和阻尼非线性[24]，采用Ansoft Maxwell进行磁力仿真，如图3和图4所示，七次多项式就可以很好地拟合仿真数据

(1)

式中：u为质量块位移；bn为磁力系数。

(a) 三维物理模型示意图

(b) 俯视图图1 EMHFig.1 EMH

图2 EMH拾振环节力学简化模型Fig.2 Simplified mechanical model of the EMH

图3 非线性磁场力拟合曲线Fig.3 Nonlinear magnetic force fitting curve

图4 拟合余差图Fig.4 Fitting residual diagram

(2)

考虑式(1)可得系统通用化的非线性振动方程

(3)

式中：my″(t)=AΩ2cosΩt，Fe=βI，分别为外界激振力与电磁线圈产生的恢复力；除此之外，m为模型质量，c1为阻尼系数，c2为非线性阻尼系数，k为刚度系数，β为机电耦合系数，I为流过负载的电流。

对式(3)进行化简，得到如下的简化方程

(4)

1.2 电路环节的数学模型

图5 EMH等效电路模型Fig.5 Equivalent circuit model of the EMH

根据基尔霍夫电流定律建立上述等效电路的0.5自由度方程

(5)

并进行简化，可得：

(6)

综合方程式(4)和(6)，可得EMH机电耦合系统1.5自由度非线性形式控制方程

(7)

其中：较为丰富的非线性形式，为拓展后续定量分析提供了通用化的研究基础。

2 基于动态频率方法的等效线性化分析

为准确把握非线性因素对振动特性的影响，本文采用动态频率方法计算其稳态响应。该方法最初用于单自由度强非线性系统的稳态响应计算与全局动力学分析，其总体求解思路是：在系统的机械能表达式之中引入一个待定的动态频率项，考虑系统中非线性因素对周期响应频率的影响，而后利用能量平衡的方法得到一系列包含未知量的代数方程组。与常规的非线性振动问题分析方法不同，此处得到的动态频率为一个关于时间变量t的函数，既有效地体现了系统中强扰动量的影响，又简化了复杂系统渐近解的求解过程，且极易实现程序化。在此基础上开展的等效线性化分析，将成为本文非线性EMH功率优化、参数分析的主要途径。

2.1 动态频率方法求解系统的稳态渐近解(Ω≈ω1，0)

为获取方程(4)的稳态响应，本文首先将文献[20]中提出的动态频率方法拓展至与式(7)类似的1.5自由度系统。以系统在平衡位置的周期运动为研究对象，给出方程(4)解的表达式

(8)

谐波分析过程中基频项在解的傅里叶函数展开式中占主要部分，考虑到高阶谐波项前的系数相比于基频项为小量，且将非线性因素对基频的影响以动态频率ω(t)的形式给出，因此式(8)可以表示为

(9)

将式(9)代入到式(6)求解出I的表达式

(10)

将式(4)转换为能量方程

(11)

式中，E0表示一个周期内的平均能量。

将式(9)和式(10)代入式(11)，积分并且整合方程中低于p阶的δ幂级数，可以得到相应的第p阶能量方程

O(δk+1)

(12)

其中：对于一阶近似，可以令式(12)中的p=1。

将方程左右两侧展开成为关于三角函数(sinT, cosT)的多项式,并将cosT的高阶项通过三角变化为sinT的高阶项，只保留含有cosT的一阶项，然后平衡方程两侧的三角函数同类项，便可得到一个关于未知变量的代数方程组，具体的平衡过程包括如下六个步骤：

(13)

联立上述代数方程组可以确定式(9)中的未知量。同时，上述平衡过程，并不会因为系统中非线性项的差异而发生变化，因而相较于谐波平衡法或者其它定量分析方法具有比较明显的应用优势。

2.2 等效线性化系统

对非线性系统式(7)进行等效线性化处理，采用与文献[19]类似的方法，引入高阶谐波组合项替代方程中的复杂非线性项形式，同时保留以待定固有频率ω1,0为基础的系统线性项部分，可得：

(14)

式中：γ1是等效阻尼系数；γ2是等效机电耦合系数；γ3是等效激励幅值；γ4是平衡常数；w是等效电流；Γi,0和Γ0,i为谐波项系数。显然，式(14)中的线性项系数较之原有非线性系统具有更为直观的物理意义，也便于经由传递函数开展后续功率分析。

以能量方程为基础给出系统式(14)的能量方程

(15)

3 EMH输出功率分析

等效线性化处理之后，EMH机电耦合系统可视为如下1.5自由度线性黏性阻尼系统

(16)

其中：激励成分f″(t)源自于等效后的外界激励与谐波项

(17)

(18)

对式(18)中的电流方程作拉氏变换

sW(s)+ReW(s)+kesX(s)=0

(19)

可解电流W(s)的表达式

(20)

对式(18)中的振动方程做拉氏变换并代入式(20)，可得：

(21)

简化式(21)，可得：

X(s)=H(s)F″(s)

(22)

式中，传递函数H(s)为

(23)

阻抗Z(s)为

(24)

由于形式上系统式(18)是线性方程，满足线性叠加原理，因此将输入f″(t)分解为各阶谐波激励，再利用传递函数获取位移响应X(t)

(25)

其中： 系数γ1,γ2,γ3,γ4的解析表达式比较复杂，此处并未列出。

为最大化系统的输出功率，首先从式(25)的线性系统位移响应出发计算系统的最大位移，而后代入式(20)确定最大电流值Wm

(26)

并由此获取系统的最大平均输出功率Pam表达式

(27)

其中：式(27)忽略了高阶谐波成分的影响。结果与文献[28]得到的功率表达式类似。

4 数值仿真分析

为验证本文方法的有效性，此处采用数值方法对原有非线性系统式(7)与等效线性化系统式(16)进行对比，其中数值模拟的参数大小如表1所示。结果如图6，7所示：等效线性化方程的数值解与原方程数值解吻合得较好，且一阶近似的等效线性化方程与原方程的余差已经很小；同时，二阶近似结果可以显著提高等效线性化系统的拟合精度。

表1 等效线性化方法的模拟参数

图6 系统(7)相图对比Fig.6 Comparison of system (7) phase diagram

此时传递函数H(s)的波特图与瞬态响应曲线如图8和图9所示，该系统是稳定的正阻尼系统，且其无阻尼固有圆频率数值为1.56，共振圆频率数值是1.47。

利用式(27)画出最大平均功率Pam的幅频曲线如图11所示，从图中可以看出功率最大点的固有圆频率数值大约为1.38与图8中得到的系统共振圆频率数值1.47相差0.09。故可将共振点圆频率ω1,0作为功率最大点便于后续计算，式(27)可化为

(a) u

图7 等效线性化方程解与原方程解的余差曲线

图8 传递函数H(s)波特图Fig.8 Transfer function H(s) Bode plot

图9 系统瞬态运动曲线Fig.9 System transient motion curve

图10 系统(16)相图对比Fig.10 Phase diagram comparison of system (16)

(28)

5 试验参数理论分析

5.1 试验装置

试验流程如图12所示。信号发生器产生的余弦信号经过功率放大器及SPEKTRA激振器作用于样机，迫使悬臂梁带动磁铁振动，从而磁场能够切割线圈，产生电流。

图13为EMH样机和试验装置示意图。样机主要包括NdFeB35磁铁，线圈由高导电的漆包线组成(N=480圈)，悬臂梁由铍青铜材料组成。试验装置主要由信号发生器，信号功率放大器，激振器，位移传感器，加速度传感器，数据处理仪等组成。

图11 固有圆频率Ω与最大平均功率Pam的关系曲线Fig.11 Relation curve between natural circular frequency Ω and maximum average power Pam

图12 试验流程Fig.12 The test procedure

1. 能量采集器； 2. 激振器； 3. 激光位移传感器； 4. 交流电阻箱； 5. 计算机； 6. 信号分析仪； 7. 电荷放大器； 8. 信号发生器；9. 功率放大器图13 试验装置和EMH样机图Fig.13 Test setup and prototype of the EMH

5.2 试验结果

由功率表达式(28)可以计算电磁单元的最大平均输出功率Pam。这里采用的非线性能量采集器的结构参数如表2所示。此处主要讨论负载电阻,机电耦合系数与系统最大平均功率Pam之间的关系。

如图14所示，结构动态响应的非线性和激励强度密切相关。激励加速度增大的过程中，结构的非线性现象越明显。从试验频响曲线可知，在0.2g激励条件下，结构响应接近线性结构响应；0.5g激励条件下，产生明显的软特性响应，即幅频曲线向低频偏移，同时响应幅值也会增加；1g激励条件下，幅频曲线进一步向低频偏移，幅值进一步增加，非线性现象更明显。当外界激励为1g时，此时的谐振频率大约为11 Hz,基于式(29)

ω=2πf

(29)

得到系统的谐振圆频率ω≈69.12 rad/s，再根据式(13)算出的ω1，0为71.48 rad/s，与实验测试相差2.36 rad/s。

根据式(26)可以求出谐振时EMH收集到最大电流为2.03 mA，而试验测得最大电压为31.04 mV，如图15所示，转换成最大电流为1.94 mA。由于振动时磁铁运动路径中磁场分布不均匀，所以输出电压呈现非线性的特征。

表2 非线性能量采集器试验参数

图14 不同加速度下样机的幅频曲线Fig.14 Amplitude frequency curves of the prototype under different accelerations

图15 谐振时电压输出波形Fig.15 Resonance voltage output waveform

当负载电阻趋近线圈电阻时，功率越接近最大值。设置外界激励圆频率为69.12 rad/s，线圈电阻为16 Ω，通过调节电阻箱改变负载电阻，此时负载R与最大平均功率Pam的关系如图16所示，负载的最大功率点大约16 Ω左右，此时最大平均功率为31.01 μW。

图16 负载R与最大平均功率Pam的关系曲线Fig.16 Relation curve between load R and maximum average power Pam

机电耦合系数β取决于线圈结构和磁通量密度[29]，如式(30)所示

β=NBL

(30)

式中：N为线圈匝数；B为磁通量密度；L为线圈长度。试验通过增加线圈匝数来增大机电耦合系数，由图17可知，随着机电耦合系数β的增大，功率也随之增大，且对功率的影响程度很大。

图17 机电耦合系数β与最大平均功率Pam的关系曲线

为了验证样机的输出特性，将其与近几年报道的电磁式振动能量采集器相对比。如表3所示，相较于其它的采集器，本文搭建的电磁式能量采集器能够在低频环境下，产生较满意的能量输出。

表3 当前文献报道的电磁式能量采集器性能比较

6 结 论

本文提出了一种基于等效线性化方法的非线性EMH功率优化分析策略。该方法有效地克服非线性因素对于EMH输出功率分析的负面影响，从而有效地提高了分析问题的效率，并且采用数值模拟以及试验手段验证了相关分析结果的有效性。本文的主要工作和结论如下：

(1) 以动态频率方法为基础可以开展通用化复杂非线性问题的研究，且此方法具有较高的求解精度；与此同时其分析结果又可以成为后续定量研究的理论基础，因而具有广泛的适用性。

(2) 利用等效线性化手段，可以弱化系统中的非线性因素，在近线性结构中更为直观地讨论关键参数变化对于系统振动特性的影响，将复杂的非线性分析简单化。本文据此开展了非线性EMH的功率分析工作，结合传递函数，有助于形成统一、便捷的非线性振动能量采集器功率研究途径。

(3) 等效线性化处理可以提高复杂非线性问题的分析效率，所得谐波项系数，在低阶近似时可以获得明确的解析表达式，也便于讨论其与原有系统参数之间的关系；而取高阶近似时，该系数就会变得比较复杂，因此如何在提高精度的同时降低复杂性仍是后续值得关注的科学问题。另一方面本文仅以特定的采集器试验模型为载体进行验证，并未采用更为小型化和高功率的采集器结构，后续需借助等效线性化理论制定优化策略来提高能量采集器效能。