风力机叶片受到阵风或随机风作用时，其会产生动载荷且振动衰减缓慢，导致叶片产生空间扭曲等变形。风力机叶片在载荷作用下的振动破坏是风力机的常见故障。振动信号采集主要采用压电传感器^[1]、PVDF压电纤维贴到风力机叶片表面检测振动^[2]，内嵌TiN涂层的光纤光栅传感器进行复合材料板的振动检测^[3]以及嵌入光纤声发射传感器(fiber optic acoustic emission sensors，FOAES)^[4-5]等方式。文献[1]采用粗纤维压电复合材料(macro fiber composite, MFC)驱动器等进行振动控制。文献[6-7]分析了转速、桨距角等对叶片固有频率等动态特性的影响，通过改变叶轮的转动方式实现叶片挥舞方向上振动的主动控制。文献[8]通过在机舱下面安装被动式隔振系统来减少风力机的动态振动响应。文献[9-11]采用弯-扭耦合(differential stiffness bend-twist coupling, DSBT)控制的方法，通过对叶片变形、攻角等的控制，实现风力机叶片在极端载荷下的自适应控制。文献[12-13]在直流高压驱动下，采用等离子体激励器(plasma actuator)改变圆柱形物体的升力和阻力系数，抑制叶片失速。文献[14-15]在叶片上安装褶翼(microtabs)进行主动力控制，并采用不同的控制算法进行对比分析，以降低叶根与叶尖的气动力集中载荷。以上研究均不同程度地抑制了叶片振动，但仍存在因需外加装置而造成结构复杂的不足。

本文针对叶片振动保护的问题，提出通过铺层设计将压电材料嵌入叶片复合材料内部以构成智能风力机叶片，该设计在一定程度上可克服现有叶片安装和制作工艺复杂的缺点。将叶片看成是复合材料层合板壳结构，建立了压电板壳式复合材料叶片结构的动力学方程，通过有限元分析研究了电压对嵌有压电材料叶片的振动抑制作用，结果表明：在风载作用下，控制电压对叶片振动有抑制作用。

1 压电板壳结构叶片设计与建模

如图 1所示，采用层合结构，将压电纤维嵌入、铺设在叶片内部，实现在叶片制作时就将压电智能材料嵌入复合材料叶片中，工艺简单，可不改变叶片的固有外部形状，构成智能叶片。

根据IEEE标准^[16]建立线性压电材料的压电本构方程：

$ \left\{ \begin{array}{l} {\varepsilon _{ij}} = S_{ijkl}^E{\sigma _{kl}} + {d_{ijk}}{E_k}\\ {D_i} = {d_{ijk}}{\sigma _{jk}} + \varepsilon _{ij}^o{E_j} \end{array} \right. $

(1)

$ \left\{ \begin{array}{l} {\sigma _{ij}} = c_{ijkl}^\mathit{{E}}{\varepsilon _{kl}} - {e_{ijk}}{\mathit{{E}}_k}\\ {D_i} = {e_{ijk}}{\varepsilon _{jk}} + \varepsilon _{ij}^t{\mathit{{E}}_j}\\ {\sigma _{ij}} = c_{ijkl}^D{\varepsilon _{kl}} - {h_{ijk}}{D_k}\\ {E_i} = - {h_{ijk}}{\varepsilon _{kl}} + \beta _{ik}^t{D_k}\\ {\varepsilon _{ij}} = s_{ijkl}^D{\sigma _{kl}} + {g_{ijk}}{D_k}\\ {E_i} = - {g_{ijk}}{\sigma _{kl}} + \beta _{ik}^\sigma {D_k} \end{array} \right. $

(2)

$ \left\{ \begin{array}{l} {d_{ijk}} = {e_{ilm}}s_{lmjk}^\mathit{{ E}} = {g_{ljk}}\varepsilon _{il}^\sigma \\ {g_{ijk}} = {d_{ljk}}\beta _{il}^\sigma = {h_{ilm}}s_{lmjk}^D\\ {e_{ijk}} = {d_{ljk}}c_{lmjk}^\mathit{{ E}} = {h_{ilk}}\varepsilon _{il}^E\\ {h_{ijk}} = {e_{ljk}}\beta _{il}^\varepsilon = {g_{ilm}}c_{lmjk}^D \end{array} \right. $

(3)

式中：ε_ij为电应变；S_ijkl^E，S_klij^D分别为材料的短路和开路弹性柔顺常数，m²/N，基于连续介质力学，材料的弹性柔顺常数与弹性刚度系数互为倒数，并与对称性有关；σ_ij为应力，N/m²；d_ijk为压电应变常数，m/V；E_i为电场强度；D_i为电位移；ε_ij^σ为应力为零时的介电常数，F/m；c_ijkl^E，c_ijkl^D分别为材料的短路和开路弹性刚度系数，N/m²；e_ijk为压电应力常数，N/Vm；g_ijkl是压电电压常数，Vm/N；h_ijk是压电劲度常数，V/m；β_ik为介质隔离率。

1.1 合成力及力矩

根据弹性扁壳内力/力矩分析以及压电本构方程，将机械能场与电场共同作用产生的压电复合材料扁壳单元中曲面单位长度上的合成力及力矩定义如下：

$ \left\{ \begin{array}{l} {N_{11}} = N_{11}^{\rm{m}} + N_{11}^{\rm{e}} = \int_\xi {}{\sigma _{11}}{\rm{d}}\xi + \int_\xi {}{d_{31}}{E_3}{\rm{d}}\xi \\ {N_{22}} = N_{22}^{\rm{m}} + N_{22}^{\rm{e}} = \int_\xi {}{\sigma _{22}}{\rm{d}}\xi + \int_\xi {}{d_{32}}{E_3}{\rm{d}}\xi \\ {N_{12}} = N_{12}^{\rm{m}} + N_{12}^{\rm{e}} = \int_\xi {}{\sigma _{12}}{\rm{d}}\xi \end{array} \right. $

(4)

$ \left\{ \begin{array}{l} {M_{11}} = M_{11}^{\rm{m}} + M_{11}^{\rm{e}} = \int_\xi {}{\sigma _{11}}{\rm{d}}\xi + \int_\xi {}{d_{31}}{E_3}{\rm{d}}\xi \\ {M_{22}} = M_{22}^{\rm{m}} + M_{22}^{\rm{e}} = \int_\xi {}{\sigma _{22}}{\rm{d}}\xi + \int_\xi {}{d_{32}}{E_3}{\rm{d}}\xi \\ {M_{12}} = M_{12}^{\rm{m}} + M_{12}^{\rm{e}} = \int_\xi {}{\sigma _{12}}{\rm{d}}\xi \end{array} \right. $

(5)

式中：N_ij^m，M_ij^m表示弹性壳体的机械力和力矩；N_ij^e，M_ij^e表示逆压电作用产生的力和力矩；上标e表示电场作用部分；上角标m表示压电壳的机械能场作用部分。

基于Hamilton原理、线性压电理论和扁壳理论基本假设推导压电板壳智能结构的系统动力学方程^[17-18]：

$ \delta \int_{{\mathit{t}_{\rm{0}}}}^{{\mathit{t}_{\rm{1}}}} {\left( {{K_{\rm{E}}}{\rm{ - }}{U_{\rm{E}}} + W} \right){\rm{d}}t = 0} $

(6)

式中：K_E为系统的动能；U_E为系统的变形能，即弹性应变能和压电应变能；W为系统的功，即边界力功W_B、外加面载荷功W_L和压电电势功W_Q；δ为变分运算算子。

1) 动能变分在时间区间[t₀, t₁]内的表达式如下：

$ \begin{array}{l} \int_{{\mathit{t}_{\rm{0}}}}^{{\mathit{t}_{\rm{1}}}} {\delta {K_{\rm{E}}}{\rm{d}}t} = {\rm{ - }}\rho h\int_{{\mathit{t}_{\rm{0}}}}^{{\mathit{t}_{\rm{1}}}} {\int_{{\mathit{a}_{\rm{1}}}}^{} {\int_{{\mathit{a}_2}}^{} {} } } ({{\mathit{\ddot u}}_1}\delta {u_1} + \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{{\mathit{\ddot u}}_2}\delta {u_2} + {{\mathit{\ddot u}}_3}\delta {u_3}){\rm{d}}{a_1}{\rm{d}}{a_2}{\rm{d}}t \end{array} $

(7)

2) 应变能变分表达式：

$ \int_{{\mathit{t}_{\rm{0}}}}^{{\mathit{t}_{\rm{1}}}} {\delta {U_{\rm{E}}}\rm{d}t} = \int_{{\mathit{t}_{\rm{0}}}}^{{\mathit{t}_{\rm{1}}}} {\int\limits_\mathit{V} {(\delta {P_m} + \delta {P_0}){\rm{d}}V{\rm{d}}t} } $

(8)

3) 系统功变分表达式：

$ \begin{array}{l} \int_{{\mathit{t}_{\rm{0}}}}^{{\mathit{t}_{\rm{1}}}} {\delta {W_{\rm{Q}}}{\rm{d}}} t = \int_{{\mathit{t}_{\rm{0}}}}^{{\mathit{t}_{\rm{1}}}} {} [\int_{{\mathit{a}_{\rm{1}}}} {\int_{{\mathit{a}_{\rm{2}}}} {({d_{31}}{s_{11}} + {\rm{ }}{d_{32}}{s_{22}} + } } \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{d_{33}}{s_{33}} + {\varepsilon _{33}}{E_3} + {Q_3})\delta \varphi {\rm{d}}{a_1}{\rm{d}}{a_2}- \\ \int_{{\mathit{a}_{\rm{1}}}} {\int_{{\mathit{a}_{\rm{2}}}} {\int_{{\mathit{a}_3}} {\frac{{\partial \left[{({d_{31}}{s_{11}} + {d_{32}}{s_{22}} + {d_{33}}{s_{33}} + {\varepsilon _{33}}{E_3}){A_1}{A_2}} \right]}}{{\partial {\mathit{a}_3}}}} } } \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left. {\delta \varphi {\rm{d}}{a_1}{\rm{d}}{a_2}{\rm{d}}{a_3}} \right]{\rm{d}}t \end{array} $

(9)

1.2 系统动力学方程

根据最终得到的压电板壳结构Hamilton方程(6) 推出压电板壳智能结构的力-电耦合系统动力学方程：

$ \begin{array}{l} \frac{{\partial \left( {N_{11}^{\rm{m}} - N_{11}^{\rm{e}}} \right)}}{{\partial {\alpha _1}}} + \frac{{\partial N_{21}^{\rm{m}}}}{{\partial {\alpha _2}}} + N_{21}^{\rm{m}}{\rm{ - }}\\ \left( {N_{22}^{\rm{m}} - N_{22}^{\rm{e}}} \right) + Q_{13}^{\rm{m}}\frac{1}{{{R_1}}} + {q_{12}} = \rho {{\mathit{\ddot u}}_1} \end{array} $

(10)

2 风载作用下振动有限元分析

根据建立的动力学方程，得到智能壳体的应力应变关系。基于机械几何方程与电学几何方程，结合本构方程及能量泛函，根据方程(10)，利用变分之任意性，得压电单元的动态有限元方程为：

$ \begin{array}{l} \left( {\mathit{\boldsymbol{M}}_m^\varepsilon + \mathit{\boldsymbol{M}}_\mathit{p}^\varepsilon } \right){{\ddot \delta }^\varepsilon } + \left( {\mathit{\boldsymbol{K}}_m^\varepsilon + \mathit{\boldsymbol{K}}_\mathit{p}^\varepsilon } \right){\delta ^\varepsilon } - \\ \mathit{\boldsymbol{K}}_{m\mathit{a}}^\varepsilon \varphi _\mathit{a}^\varepsilon - \mathit{\boldsymbol{K}}_{m\mathit{z}}^\varepsilon \varphi _\mathit{z}^\varepsilon = \mathit{\boldsymbol{F}}_1^\varepsilon + \mathit{\boldsymbol{F}}_2^\varepsilon \end{array} $

(11)

式中：

$ \begin{array}{l} \mathit{\boldsymbol{M}}_m^\varepsilon = \smallint \mathit{\boldsymbol{N}}_\delta ^{\rm{T}}{\rho _m}{\mathit{\boldsymbol{N}}_\delta }d{{\rm{V}}_{\rm{p}}}\\ \mathit{\boldsymbol{M}}_\mathit{p}^\varepsilon = \smallint \mathit{\boldsymbol{N}}_\delta ^{\rm{T}}{\rho _p}{\mathit{\boldsymbol{N}}_\delta }{\rm{d}}V\\ \mathit{\boldsymbol{K}}_m^\varepsilon = \int\limits_{{V_m}} {} \mathit{\boldsymbol{B}}_\delta ^{\rm{T}}{\mathit{\boldsymbol{C}}_m}{\mathit{\boldsymbol{B}}_\delta }{\rm{d}}{V_p}\\ \mathit{\boldsymbol{K}}_\mathit{p}^\varepsilon = \int\limits_{{V_m}} {} \mathit{\boldsymbol{B}}_\delta ^{\rm{T}}{\mathit{\boldsymbol{C}}_p}{\mathit{\boldsymbol{B}}_\delta }{\rm{d}}V \end{array} $

M_m(M_p)和K_m(K_p)分别为压电层与基体层的单元耦合质量和刚度矩阵，其中，ε为压电层的介电常数。

$ \left\{ \begin{array}{l} \mathit{\boldsymbol{F}}_{\rm{1}}^\varepsilon = \int\limits_{{V_m}} {} {\mathit{\boldsymbol{N}}_\delta }{f^a}{\rm{d}}V + \int\limits_{{S_a}} {} \mathit{\boldsymbol{N}}_\delta ^{\rm{T}}{\mathit{\boldsymbol{T}}^\varepsilon }{\rm{d}}s\\ \mathit{\boldsymbol{F}}_2^\varepsilon = \sum \int\limits_{{S_q}} {} \mathit{\boldsymbol{N}}_\delta ^{\rm{T}}Q{\rm{d}}s \end{array} \right. $

(12)

式中：F₁^ε为面和体机械力向量，F₂^ε表示电场力向量。在动力学问题中阻尼起着重要的作用，引入瑞利阻尼阵，表征压电层合壳的固有阻尼，由式(11) 得：

$ \mathit{\boldsymbol{M\ddot X}} + {\mathit{\boldsymbol{C}}_R}\mathit{\boldsymbol{\dot X}} + \mathit{\boldsymbol{KX}} = {\mathit{\boldsymbol{F}}_1} + {\mathit{\boldsymbol{F}}_a} $

(13)

$ {\mathit{\boldsymbol{C}}_R} = A\mathit{\boldsymbol{M}} + B\mathit{\boldsymbol{K}} $

(14)

式中A，B为瑞利阻尼系数。

低阶振型在工程结构的动力反应中起主导作用^[19]，通常取低阶振型来确定瑞利阻尼系数A和B，考虑叶片前2阶振型的频率，通过模态分析计算得叶片前2阶固有频率：ω₁=1.325 Hz，ω₂=2.753 Hz，联立方程组：

$ \left\{ \begin{array}{l} A/\left( {2{\omega _1}} \right) + B \times {\omega _1}/2 = \xi \\ A/\left( {2{\omega _2}} \right) + B \times {\omega _2}/2 = \xi \end{array} \right. $

将ξ=0.04代入方程组中，求得A=0.072，B=0.02。

取压电材料为PZT-4^[20]，其性能参数如下：

$ \begin{array}{l} {C_{1111}} = {C_{2222}} = 139{\rm{\; GPa}}\\ {C_{3333}} = 115{\rm{\; GPa}}, \\ {C_{1122}} = 77.8{\rm{\; GPa}}\\ {C_{1133}} = {C_{2233}} = 74.3{\rm{\; GPa}}\\ {C_{2323}} = {C_{1313}} = 25.6{\rm{\; GPa}}\\ {C_{1212}} = 30.6{\rm{\; GPa}}\\ {e_{311}} = {e_{322}} = - 5.2{\rm{ \;C/}}{{\rm{m}}^{\rm{2}}}\\ {e_{113}} = {e_{223}} = 12.7{\rm{\; C/}}{{\rm{m}}^{\rm{2}}}\\ {e_{333}} = 15.1{\rm{\; C/}}{{\rm{m}}^{\rm{2}}}\\ {\varepsilon _{11}} = {\varepsilon _{22}} = 13.06{\rm{\; nF/m}}\\ {\varepsilon _{33}} = 11.51{\rm{\; nF/m}} \end{array} $

表 1为玻璃钢材料参数，该材料具有沿纤维方向阻尼小，剪切方向阻尼大的特点。

根据表 1和压电材料的参数，所设计的风力机功率为11 kW，额定风速为8 m/s；设计的叶片长度为5.5 m。翼型选用NACA4412，其最大厚长比为12%，气动中心在距前缘30%处，叶片的几何形状如图 2所示，实际安装角θ_i如表 2所示。

叶片载荷工况取为极限工况，叶片迎面受到的均布载荷^[21-22]为：p=0.5 ρv²，其中ρ为空气密度, v为风速，以极限工况计算取v=45 m/s，ρ=1.225 kg/m³，则叶片迎风面受到的均布载荷为1.24 kPa。结合结构几何参数和荷载参数，采用MATLAB语言完成计算程序的编写和调试。给定初始条件，通过多步迭代得到最终解。

图 3、图 4和图 5分别是叶根一阶y向位移分量、叶尖二阶y向位移分量和叶尖部节点绕x轴的转角。

由于风力机振动能量主要集中在低阶，对前2阶模态实施控制，并观察控制效果，比较控制电压施加前后振动位移的变化，从图中可知位移的变化较为明显。在图 3中，控制电压施加前y向位移分量，最大振动幅值为0.4 mm，控制电压施加后，叶片振动快速衰减，相比施加电压前，8 s后已衰减80%。

仿真结果表明：施加控制电压，能够在一定程度上改善控制效果，振动抑制效果明显。

3 结论

本文将风力机叶片等效为2个矩形扁壳的结构形式，建立了压电板壳式复合材料叶片结构的机-电耦合动力学模型，利用有限元法并结合MATLAB软件模拟分析了风载作用下压电材料在智能结构中的振动抑制性能。结果表明：在风载作用下，压电材料的应用能有效地抑制了叶根和叶尖的位移，并可有效降低作用前期的叶尖部转角，控制电压对叶片振动有明显的抑制作用，显著提高了风力机叶片的气动弹性稳定性，对提高叶片使用寿命有重要意义。