基于量子粒子群算法的風力機葉片損傷識別

孫芳錦， 盧 琛

(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檢測算法未來在實際應用工程風機損傷工況檢測領域具有廣闊的技術應用性和前景。