風電機組葉片覆冰狀態與動力學特性的定量關系研究*

文 | 陳志剛，李錄平，楊波，劉勝先，封江，龔妙

我國的疆域遼闊，氣象條件復雜，在北部地區的陸地和海上，南部地區的高海拔山區等，冬季的風雪霜凍給風電機組的運行帶來了極大危害。運行在冬季濕冷地區的風電機組，在冬季經常會出現嚴重覆冰現象。工作中的風電機組葉片覆冰后，會破壞風電機組葉片的動力特性與氣動特性，冰的質量、剛度等性質會影響葉片的材料特性，造成葉片過載及載荷不均勻，使其特性與設計初衷不符，造成風電機組不能正常發電，甚至導致嚴重的安全事故。及時、準確地檢測出葉片的覆冰狀態，便于對風電機組采取及時有效的運行和維護措施，對提高風電場運行的安全可靠性和經濟性，具有重要意義。

自20世紀90年代開始，對于構件覆冰檢測的研究已經引起人們廣泛關注，研究對象主要集中在航空領域（如飛機），主要的測量方式包括：機械法、電學法、熱學法、光學法、波導法等。但是，目前直接用于風電機組葉片覆冰檢測的研究仍然比較少，能夠成熟用于工程實際的風電機組葉片覆冰檢測方法與技術仍少有資料報道。

本文探索一種利用葉片模態參數變化特征來檢測葉片覆冰狀態的方法：以翼展1120mm（不包含葉片根部）長的小風電機組葉片為研究對象，在Solidworks平臺上建立起三維實體模型。建立在振動力學、動力學分析以及柔性多體動力學理論基礎上，在ANSYS結構分析軟件中對該葉片在覆冰前后狀態下模態進行分析，得到前6階固有頻率和振型，比較覆冰前后風電機組葉片動力特性的變化情況，提取葉片覆冰后的振型曲率變化特征，獲得診斷葉片覆冰程度、覆冰位置的基本方法，為工程實際中診斷葉片覆冰故障提供理論依據。

研究對象簡介

本文以一小型風電機組葉片為研究對象，在三維制圖工具軟件Solidworks中建立風電機組葉片三維實體模型（見圖1）。如圖1所示，葉片為弦長140mm的等截面直葉片，內部為實體結構，長度為1120mm，由一種翼型界面拉伸而成。

圖1 模擬風電機組葉片三維實體模型

圖2 葉片結冰沿翼展方向分布示意

葉片覆冰狀態描述模型

風電機組葉片覆冰的空間分布描述采用如圖2所示的模型，假設葉片在某一區段附著一定厚度的冰，其覆冰狀態（Blade Icing Condition，簡稱BIC）信息用下式描述：

式中，x是葉片覆冰段沿翼展方向的起始坐標（相對坐標值）；Δx是葉片覆冰段沿翼展方向的長度（相對值）；λ是葉片覆冰段的覆冰厚度（相對坐標值）。

一、覆冰區域起始點坐標參數

定義葉片覆冰的起點坐標參數為x（見圖2），其表示的是沿翼展方向覆冰區域的起始點相對位置，參數為一個無量綱量。以葉片與輪轂的結合處為翼展方向坐標原點，定義描述覆冰起始位置的坐標參數計算公式為：

式中，L是所測風電機組葉片的翼展長度，m；X是覆冰段起始點距離葉根的坐標，m；x的取值范圍為：0～1。

二、覆冰區域長度參數

定義覆冰長度參數為Δx，其表示在翼展方向覆冰區域沿翼展方向的相對長度。設覆冰段起點端到另一端的長度為 ΔX（見圖2），定義描述覆冰區域相對長度參數的計算公式為：

式中，ΔX是覆冰段沿翼展方向的長度，m。Δx的取值范圍為：0～1。

三、覆冰厚度參數

根據風電場的風能資源條件，風電機組葉片的結構大小多種多樣，因此用冰層厚度的絕對值不能準確反映風電機組葉片覆冰厚度。因此，定義覆冰厚度相對值λ，用其描述覆冰區域平均覆冰厚度：

式中，d是所測風電機組葉片覆冰區域平均覆冰厚度，m。

覆冰葉片動力學特性計算數學模型

一、覆冰葉片等效材料特性參數計算模型

采用有限元方法計算覆冰葉片的動力學特性前，需要對比覆冰前與覆冰后葉片材料性能參數的變化情況。為此，本文在計算葉片中取中間一段180mm長的葉片段，建立葉片片段覆冰模型，覆冰相對厚度λ取值范圍為0～8.035714（1～9mm厚度）。根據所設置的葉片材料密度和冰層密度，用Solidworks軟件計算出葉片結構的質量變化情況。圖3、圖4所示為葉片段相對覆冰厚度λ分別取值為0.0、2.0時的覆冰模型。

圖3 葉片片段λ為0的模型

圖4 葉片片段λ為2的覆冰模型

圖5 葉片片段結構網格化圖

圖6 葉片片段λ為2的結構網格化圖

將覆冰前后的葉片段結構導入ANSYS中，分別進行力學分析，在覆冰葉片段組合體中定義葉片材料的彈性模量E＝4.26e10，泊松比為0.22，密度為1950kg/m3；冰的彈性模量E＝2e9，泊松比為0.1，密度為900kg/m3。由于計算對象為小葉片，所以定義葉片為Solid45實體模型，網格化以后得到圖 5－圖6。

定義葉片段的一端為全約束，在葉片的另一端Y軸反方向施加一個50N力載荷，載荷施加圖和軟件計算后的結構變形圖如圖7、圖8所示。

在葉片段的一端設置全約束，另一端施加一個力偶載荷產生5N?m力矩，通過ANSYS軟件計算得到覆冰葉片結構的扭轉變形圖，如圖9所示。通過施加力和扭矩（力偶）分別得到葉片的彎曲變形和扭轉變形，在施加外界載荷中力和扭矩不變的條件下，分別對不同覆冰相對厚度的葉片進行結構受力變形計算。

本文采用的抗彎剛度實驗測量公式為：

采用的扭轉剛度實驗測量公式為：

式中，F為結構所受力載荷，N；M為結構所受力偶載荷， N?m；y為結構在力載荷下的位移，m；θ為結構在力偶載荷下的扭轉角度，deg。

首先根據式（5）、（6）計算得到覆冰葉片片段結構的抗彎剛度和扭轉剛度，結合質量變化、抗彎剛度及扭轉剛度計算得到三項參數變化率隨覆冰厚度增加的變化關系圖，如圖10所示。

從圖10中可以看出，葉片的相對覆冰厚度λ從0到8之間變化時，葉片材料等效性能參數呈現如下變化規律：

（1）葉片覆冰后，材料等效質量、等效彎曲剛度、等效扭轉剛度均有不同程度增加；

（2）葉片的等效質量隨覆冰的相對厚度幾乎呈線性增加，當相對覆冰厚度λ達到8時，葉片的等效質量增加近35%左右；

（3）當覆冰相對厚度λ在 4.5以內時，葉片的等效彎曲剛度隨覆冰厚度的增加呈現快速增加趨勢，之后，隨著覆冰相對厚度的繼續增加，等效彎曲剛度增加的趨勢減緩，當覆冰厚度達到8時，葉片的等效彎曲剛度增加近35%左右；

（4）當覆冰相對厚度λ在4.5以內時，葉片的等效扭轉剛度隨覆冰厚度的增加呈現緩慢增加趨勢，之后，隨著覆冰相對厚度的繼續增加，等效扭轉剛度增加的趨勢加快，當覆冰厚度λ達到8.0時，葉片的等效扭轉剛度增加近12%左右。

圖7 葉片片段λ為0的彎曲變形圖

圖8 葉片片段λ為2的彎曲變形圖

圖9 葉片片段λ為2的扭轉變形圖

圖10 葉片片段在不同覆冰相對厚度下的材料性能變化圖

二、覆冰葉片模態參數變化量計算模型

（一）葉片模態頻率變化量計算模型

將葉片在揮舞方向的振動問題簡化成在根部固接的懸臂彎曲振動問題。懸臂梁的模態頻率計算問題可演化成由下式描述的特征值問題：

式中，K和M分別是葉片整體剛度矩陣和質量矩陣；φ是正則化振型；ω是模態頻率。假定風電機組葉片覆冰使葉片結構剛度矩陣和質量矩陣都產生一個細微攝動量，則對φ和ω也會產生一個細微改變量，改變后的振動方程為：

式中，ΔM，ΔK和Δφ分別為葉片結構的質量矩陣、剛度矩陣和振型的改變量。

對風電機組葉片這樣的變截面懸臂梁結構，在不同的覆冰厚度情況下，覆冰對質量矩陣的影響程度和剛度矩陣影響程度是不同的。對式（8）適當變形，得到：

在某一階振型下，有：

式（9）、（10）中，φT表示正則化振型矩陣φ的轉置矩陣。從上兩式可以分析得出，若剛度增加則葉片結構自振頻率增加，剛度降低則結構頻率會降低；而質量增加，結構頻率就會降低，反之，結構頻率會增加。

（二）葉片模態振型曲率計算模型

由于風電機組葉片為懸臂梁結構，因為揮舞方向的抗彎能力最差，所以變形主要以揮舞方向彎曲變形為主，因此采用梁結構單元進行推導，其振動微分方程為：

式中，x表示梁的軸向方向的坐標；I（x）表示梁的抗彎截面模量；E、ρ分別表示材料的楊氏模量和密度；A（x）表示梁的截面積；C（x）表示材料的阻尼；y（x,t）表示梁的變形（位移）；f（x,t）表示作用在梁上的外力。

根據模態分析理論，式（11）的解可表示為各階模態振型的疊加形式：

式中，φi（x）為位移模態振型；qi（t）為模態坐標；i為振型階數。

根據梁結構單元振動微分方程式（11）和振動位移公

式（12），將式（11）化簡得到：

式中，ms、cs、ks分別表示梁的第s階模態質量、阻尼和剛度；φs表示梁的第s階模態振型；qs表示梁的第s階模態坐標。

根據材料力學給出的梁結構靜力彎曲關系：

式中：M是葉片彎矩，ρ是曲率半徑。

設f（x,t）＝F（x）ejωt，由上兩式得到 ：

曲率是位移的二階導數，將上式離散化后得到：

矩陣形式為：

其中：

y′′的微分增量為 ：

覆冰葉片模態特性變化規律計算與分析

一、覆冰葉片模態頻率變化量計算與分析

在Solidworks中建立葉片覆冰模型。圖11－圖13分別表示在葉片不同覆冰起點位置上覆上一定長度區域的冰，并對不同覆冰狀態下的葉片進行模態分析。

采用ANSYS軟件，對不同覆冰情況下的葉片進行模態分析，得到頻率以及振型數據。葉片與冰的材料參數定義同上，將所有頻率數據整理出來，計算出相對于未覆冰狀態下葉片結構的頻率變化率（圖14）。圖14中曲線說明有三個參數，第一個參數表示覆冰起點x，第二個參數表示覆冰長度Δx，第三個參數表示相對覆冰厚度λ。

圖14中包含所測葉片結構前六階模態頻率。其中，從第一階頻率至第六階頻率分別為揮舞一階振型頻率、擺振一階振型頻率、揮舞二階振型頻率、揮舞三階振型頻率、扭轉一階振型頻率和揮舞四階振型頻率。

圖11 覆冰起點x＝0.25，覆冰長度Δx＝0.5，覆冰厚度λ＝4時葉片三維圖

圖12 覆冰起點x＝0.25，覆冰長度Δx＝0.25，覆冰厚度λ＝4時葉片三維圖

圖13 覆冰起點x＝0.75，覆冰長度Δx＝0.25，覆冰厚度λ＝2時葉片三維圖

從圖14可以發現：

（1）葉片表面覆冰后，會導致模態頻率的變化；覆冰越嚴重，模態頻率的相對變化值越大。

（2）不同覆冰狀態下各階模態頻率的變化率也是不同的。整體而言，葉片覆冰對一階模態頻率影響最大。

（3）當葉片所有表面覆冰時，葉片模態頻率變化最大且與覆冰厚度成正比關系。

（4）葉片根部覆冰后導致葉片結構剛度增加明顯，故模態頻率升高；葉片葉尖區域覆冰后導致葉片結構質量增加明顯，引起模態頻率下降。

圖14 不同覆冰狀態下葉片結構前六階模態頻率變化

圖15 葉片各點在不同覆冰條件下一階曲率變化率

二、葉片覆冰狀態與葉片模態振型的關系

由于葉片覆冰后，結構剛度顯著增大，導致的曲率變化［Δy′′］，主要由［Δφ′′r］［ΔYr］［Δφr］三者共同作用。而［Δy′′］ 與［Δφ′′r］的變化在葉片位置坐標上存在統一對應性的變化。

用ANSYS軟件仿真計算出的葉片結構固有特性參數包括模態頻率和模態振型，所獲得的模態振型主要以位移模態數據為主。根據差分法計算出葉片結構的各階曲率模態參數。這里主要分析葉片結構的揮舞方向的模態曲率，根據下式得出葉片結構的曲率變化率曲線圖，如圖15所示。

式中：ξi為i點處覆冰后的曲率變化率；φwi為i點處覆冰后曲率值；φfi為i點處無冰狀態下曲率值。

從圖15中可以看到：

（1）在葉片全覆冰情況下，覆冰2mm時曲率變化較小，曲率變化率曲線無規則；當覆冰達到4mm時曲率變化明顯增大，同時曲率變化率曲線呈現明顯的拋物線形狀。由此可見，在覆冰區域一定時，覆冰厚度越大，葉片模態振型曲率的變化率越大。

（2）覆冰起點為0.5，覆冰長度分別為0.25和0.50時，葉片對應位置的曲率變化率曲線也有凸起形狀，同時可以看出覆冰長度Δx＝0.50的曲率變化率大于覆冰長度Δx＝0.25的曲率變化率。由此可見，在覆冰起點位置相同時，覆冰區域越大，葉片模態振型曲率的變化率越大。

從上述分析結論可以看出，通過對覆冰前后葉片結構的曲率變化進行分析研究，能夠進一步對覆冰狀態進行判斷。

攝影：李博

結論

本文提出了風電機組葉片覆冰狀態的空間描述模型，給出了葉片覆冰起點參數、覆冰長度參數和覆冰厚度參數的定義；給出了覆冰葉片動力學參數變化量計算模型。以一小型風電機組葉片為例，在Solidworks三維建模軟件中建立三維實體模型，應用ANSYS結構分析軟件對風電機組葉片進行有限元模型分析，通過有限元計算，獲得葉片在不同覆冰狀態下的模態頻率變化率、模態振型曲率變化率，通過分析發現：

（1）葉片覆冰區域靠近葉尖部分時，葉片模態頻率呈現下降趨勢；葉片所有表面覆冰時，模態頻率變化率與覆冰厚度成近似正比關系。葉片根部區域覆冰后，葉片結構以剛度增加為主，故模態頻率升高。葉片葉尖區域覆冰后，葉片結構以質量分布變化為主，導致模態頻率下降。

（2）葉片模態振型曲率變化率值與葉片覆冰位置之間關系密切，根據葉片模態振型曲率變化率的不同，能夠在一定程度上對不同覆冰情況下葉片覆冰位置和覆冰厚度進行判斷。

（3）葉片模態頻率變化率、模態振型曲率變化率與葉片覆冰狀態具有明顯的定量關系。工程實際中，通過這種定量關系，可以診斷葉片的覆冰故障。