基于量子粒子群算法的风力机叶片损伤识别

孙芳锦， 卢 琛

(1.桂林理工大学 土木与建筑工程学院 广西 桂林 541004； 2.广西嵌入式技术与智能系统重点实验室， 广西 桂林 541006)

最近几十年来,受自然灾害和养护不善、爆炸、冲击等人为因素的影响,建筑结构发生不同程度破坏和一定程度功能丧失,结构的健康监测显得格外关键。用于结构损伤鉴定的重要技术指标有:固有频率、模态结构振型、柔度矩阵的变化、频响函数、神经网络等。大部分结构破坏的本质是结构刚度的降低，柔度和刚度成倒数关系，可以用作结构鉴定损伤的指标。

近年来，柔度矩阵与智能算法相结合的结构损伤识别技术成为国内外科学研究的热点。许多学者使用柔度矩阵差、柔度差变化率和柔度曲率差来定位结构的损伤[1-2]，但利用柔度灵敏度[3]和广义柔度灵敏度[4]对结构损伤进行定量识别时，存在稳定性差、误差大的缺点。Zhao等[5]研究频率、振型和柔度的灵敏度,证明了柔度矩阵比频率和振型更为灵敏。

本文结合量子粒子群优化算法和广义柔度矩阵提出了一种新的损伤识别方法。

1 量子粒子群算法

Eberhart与Kennedy在1995年首先提出了粒子群算法(PSO)[6]，粒子群算法的数学模型是假设N维解空间结构中存在M个微粒。第m个粒子在t时间上位移和移动速度分别为

Xm(t)=(Xm,1(t),...Xm,n(t),...Xm,N(t)),Vm(t)=(Vm,1(t),...Vm,n(t),...Vm,N(t)),m=1,2,...M,n=1,2,...N。

当前时间粒子群中的局部区域最优预测位置和全域最优预测位置分别为

Pm(t)=(Pm,1(t),...Pm,n(t),...Pm,N(t)),G(t)=(G1(t),...Gn(t),...GN(t))。

其中

(1)

G(t)=Pg(t)=(Pg,1(t),...Pg,n(t),...Pg,N(t))

(2)

式中：f[·]表示构建的目标函数适应值g=argmin{f[Pm(t)]},1≤m≤M；g为一个位于全局最优位置粒子对应的下标。

随着时间t的增长，在粒子群算法中M个微粒的位置变化和速率更新，分别为：

Vm,n(t+1)=Vm,n(t)+c1u1(t)[Pm,n(t)-Xm,n(t)]+c2u2(t)[Gn(t)-Xm,n(t)]

(3)

Xm,n(t+1)=Xm,n(t)+Vm,n(t+1)

(4)

式中：c1和c2分别表示为局部最优和全局最优加速因子,通常取(0,2)中的一个值；u1和u2是均匀分布在区间(0,1)中的随机数。为了防止粒子在迭代过程中飞出搜索空间，必须依次确定Vm,n(t)、Xm,n(t)搜索的上、下界。

粒子群优化算法收敛速度慢，搜索模式单一，空间受限，易局限于局部优化。为了解决的上述缺点，科学家们不断探索。Sun在量子力学的理论角度上给出了量子粒子群算法(QPSO)[7]，融入了波函数与δ势阱的有机组合理念等，确保了搜索空间可以充分涵盖整个区域。QPSO搜索算法把粒子在空间出现的机率以波函数模的平方表示，从而保证以全机率寻找到全局的最优解，通过Monte Carlo方法得

(5)

式中：p是势阱吸引子；L是δ势阱的特征长度；u是均匀分布在区间(0,1)中的随机数，即u∈U(0,1)。

随着时间t的增加，QPSO算法中M个粒子的位置和速度更新如下：

Xm,n(t+1)=pm,n(t)±α*|Cn(t)-Xm,n(t)|*ln(1/um,n(t))

(6)

pm,n(t)=βn(t)*Pm,n(t)+(1-βn(t))*Gn(t)

(7)

mbest=(C1(t),C2(t),...CN(t))

(8)

QPSO算法的步骤：

步骤1：初始化粒子个体位置，并设置初始粒子个体的最优值等；

步骤2：计算各粒子的适应度值；

步骤3：比较各粒子的适应度值大小，然后更新个体最优与全局最优；

步骤4：运用公式(5)～公式(7)更新粒子的位置和速度，得到最新的位置和速度；

步骤5：重复步骤2～步骤4，直至达到停止的条件。

2 广义柔度矩阵的基本原理

无阻尼振动条件下，n个自由度结构在自由振动时的方程为

[M]x″+[k]x=0

(9)

式中：[K]、[M]分别表示结构的刚度矩阵和质量矩阵；x″、x分别代表结构的加速度响应和位移响应。

由(9)式的运动方程，得到特征方程为：

([K]-λ[M])φ=0

(10)

式中：λ、φ分别表示方程特征值、特征向量；λ表示结构各阶固有频率的平方。

将式(10)以矩阵的形式表示为：

[K][Ψ]=[M][Ψ][Λ]

(11)

式中：[Ψ]=[φ1,φ2,...φn]为n阶模态振型矩阵；[Λ]=daig(λ1,λ2,...λn)为n阶对角矩阵。

由振型的质量归一化特性得

[Ψ]T[M][Ψ]=[I]

(12)

式中：[I]为单位矩阵。

将式(11)两边同时乘[Ψ]T得

[Ψ]T[K][Ψ]=[Ψ]T[M][Ψ][Λ]

(13)

式(12)代入式(11)的右侧得

[Ψ]T[K][Ψ]=[Λ]

(14)

式(14)两侧同时求逆，整理得

[K]-1=[Ψ][Λ]-1[Ψ]T

(15)

其中柔度矩阵是刚度矩阵的逆矩阵，即：

(16)

从方程(16)可以看出，柔度矩阵是关于结构振动频谱和强度振型的向量函数，与频率倒数的平方成正比。随着模态阶数的增加，模态参数对柔度矩阵的贡献越来越小。模态柔度仅通过前几阶模态就能更好地估计结构柔度矩阵[2]，有效地克服了工程设计中仅检测前几阶模态和高阶模态不准确的问题。

为了引入广义柔度矩阵的概念[4]，在柔度矩阵基础上进行改进，具体函数为：

(17)

式中：Q=0,1,2,...。当Q=0时，[G]=[F]，广义柔度矩阵与一般柔度矩阵相等；当Q=1时，[G]=[Ψ][Λ]-2[Ψ]T。

3 目标函数的建立

根据式(17)建立结构损伤前后广义柔度矩阵差为

(18)

式中：[Gu]、[Gd]分别代表结构损伤前后广义柔度矩阵；φru、φrd分别代表结构破坏损伤前后模态振型；ωru、ωrd分别代表结构破坏损伤前后模态固有频率；m为小于n的正整数。

(19)

式中：θi为第i个单元损伤参数因子。

令

(20)

(21)

构造目标函数

Z(θ)=A-B

(22)

式中：θ=(θ1,θ2,...θn)，θ的取值范围为0≤θ≤1。考虑到QPSO算法粒子搜索的最佳取值范围和结构损伤大于90%时将丧失使用功能，将θ的取值范围调整为0.1≤θ≤0.9。

变刚度法是将损伤部位的各单元的刚度减小，再与其他常规单元共同进行计算。根据圣维南定律，这种损伤只能影响到周围的区域，而不会影响到更远的地方。如果单元选择得足够大，则损伤部位的单元刚度仅会影响到损伤部位的单元刚度，其他单元的刚度不会受到影响。用变刚度法模拟了叶片的损伤，在损伤单元的纵向、横向和剪切模量均减少的情况下，利用各部位的刚度变化来模拟损伤部位的损伤，以同一部位的刚度变化来反映损伤的严重性。假定E和E*分别为叶片完好和损伤两种状态下的弹性模量，则反映损伤程度的单元刚度折减系数x可表示为：

(23)

4 基于QPSO算法对风力机叶片的损伤识别

叶片是风力机的主要部件之一，很容易恶化和磨损，影响使用寿命。进行定期检查，在早期阶段发现故障，能减轻这些损失。这里采用一种以损伤前后结构的广义柔度矩阵差构造目标函数来识别结构损伤的方法。

选取20 kW风力机单独的一个5.532 m叶片如图1所示,将该叶片划分为92个单元如图2所示，叶片材料属性与参数如表1所示，实验数据[8]如表2所示。

图1 20 kW风力机叶片模型

图2 20 kW风力机叶片网格划分示意图

表1 20 kW风力机叶片材料属性与参数

表2 20 kW风力机叶片实验数据

迭代次数=评价次数/种群数量[9-11]。在相同的种群数量下，评价次数与迭代次数成正比，达到相同的适应度值，看所需的评价次数来评估算法的计算效率。即评价次数越少，迭代次数越少，耗时越少，计算速度越快。为了证明研究的实效性以及计算的优越性，分别按变刚度法折减刚度为30%、50%、70%设置损伤工况，利用ABAQUS软件提取单元刚度矩阵、整体刚度矩阵、质量矩阵代入构造的目标函数计算得到相应损伤工况的目标函数，再用MATLAB得到每个损伤工况的单元损伤参数、评价次数与适应度值的折线图。分别用PSO、QPSO算法来进行3种工况的损伤识别，如图3～图5所示。

(a)PSO算法 (b)QPSO算法图3 工况1下20 kW风力机叶片在不同算法下的收敛折线图

(a)PSO算法 (b)QPSO算法图4 工况2下20 kW风力机叶片在不同算法下的收敛折线图

(a)PSO算法 (b)QPSO算法图5 工况3下20 kW风力机叶片在不同算法下的收敛折线图

选取5 MW风力机单独的一个61.5 m叶片如图6所示，将叶片划分为55 366个单元如图7所示，结构物理参数如表3所示。

图6 5 MW风力机叶片模型

图7 5 MW风力机叶片网格划分示意图

表3 5 MW风力机叶片材料属性与参数

同20 kW风力机叶片的计算方法，得到识别结果如图8～图10所示。

(a)PSO算法 (b)QPSO算法图8 工况1下5 MW风力机叶片在不同算法下的收敛折线图

(a)PSO算法 (b)QPSO算法图9 工况2下5 MW风力机叶片在不同算法下的收敛折线图

(a)PSO算法 (b)QPSO算法图10 工况3下5 MW风力机叶片在不同算法下的收敛折线图

5 结果分析

20 kW和5 MW风力机叶片在不同工况下两种算法的收敛坐标如表4和表5所示。

表4 20 kW风力机叶片不同算法在不同工况下的收敛坐标

表5 5 MW风力机叶片不同算法在不同工况下的收敛坐标

在MATLAB中查看结果，每个工况每个单元损伤情况一致，只需设置一个单元损伤参数，对20 kW和5 MW风力机叶片由两种算法在不同工况下单元损伤参数如表6和表7所示。

表6 20 kW风力机叶片不同算法在不同工况下的单元损伤参数

表7 5 MW风力机叶片不同算法在不同工况下的单元损伤参数

由图3～图5、图8～图10和表4～表7分析得：

(1)QPSO算法比PSO算法识别损伤程度更精确，精度提高接近2%。

(2)由图3可知刚开始适应度值达到0.01时，2种算法的评价次数相差不大，随着适应度值达到0.001时，PSO算法评价次数达到QPSO算法的7倍，即迭代次数达到QPSO算法的7倍。收敛时的适应度值大概为0.000 1左右，PSO算法评价次数达到QPSO算法的6倍。在图8中适应度值达到0.001时，两种算法的评价次数相差不大，适应度值达到0.000 1时，PSO算法评价次数达到QPSO算法的10倍。说明当达到相同的适应度值时，QPSO算法比PSO算法收敛更快，计算效率更快。

(3)由图4可知随着风力机叶片损伤程度加大后，当适应度值为0.001左右时，PSO算法评价次数达到QPSO算法的3倍，即迭代次数达到QPSO算法的3倍。在图9中当适应度值为0.001左右时，PSO算法评价次数达到QPSO算法的2倍，即迭代次数达到QPSO算法的2倍，可以看出QPSO算法计算效率仍然大于PSO算法。

(4)随着叶片损伤程度达到70%，当适应度值为0.01左右时，图5中PSO算法评价次数达到QPSO算法的1倍左右，即迭代次数达到QPSO算法的1倍左右。图10中当适应度值为0.01左右时，PSO算法评价次数达到QPSO算法的11倍左右，即迭代次数达到QPSO算法的11倍左右。说明随着叶片损伤程度增加，QPSO算法对于大型风力机叶片的计算效率远大于PSO算法。

6 结 语

通过算法检测不同风机损伤程度工况的同型风力机不同叶片,且虽然风力机不同叶片分别划分的检测单元多与杂,具有一定的不准确性,但采用QPSO检测算法与PSO检测算法都仍然能得到较好的准确识别,表明采用QPSO检测算法目前具有良好的鲁棒性,QPSO检测算法未来在实际应用工程风机损伤工况检测领域具有广阔的技术应用性和前景。