參數耦合下風力發電機葉片機械顫振檢測研究

張 純，孫 山，韓煜航

（1.西安熱工研究院有限公司，西安 710054；2.西安交通大學 能源與動力工程學院，西安 710049）

風力發電機是轉化氣流動能為機械能，最終轉化為電能的主要動力機械設備[1]，目前已成為不可或缺的可再生能源發電設備，在全球風電產業不斷發展下[2]，風力發電機的應用場景越發廣泛，但同時暴露出許多技術問題。葉片機械顫振不僅會影響風力發電機組的運行效率和可靠性，還可能導致嚴重的結構損壞和安全事故。因此，對于風力發電機葉片機械顫振的檢測和控制顯得尤為重要。

文獻[3]采用非線性諧波法對風力發電機葉片實行三維非定常流動分析，得到葉面脈動壓力，通過流固轉換程序獲取葉片所受氣動激振力。采用模態疊加法實現機械顫振檢測。文獻[4]通過構建風力發電機葉片三維殼體模型，采用有限元理論建立風力發電機葉片運動學方程，實現機械顫振檢測。以上方法均存在一定的局限性，導致氣動阻尼系數計算結果不理想、檢測正確率低、計算時間長的問題。

為了解決上述方法中存在的問題，提出參數耦合下風力發電機葉片機械顫振檢測方法。

1 風力發電機葉片結構動力學方程推導

基于懸臂旋轉歐拉－伯努利梁運動微分方程和哈密頓原理建立風力發電機旋轉葉片結構模型[5-6]，如下所示：

式中：U 表示應變能；T 表示動能；W 表示虛功；δ 表示變分。

考慮復合材料正交各向異性，改進風力發電機葉片揮舞彎曲和扭轉控制微分方程為如下形式：

式中：x 表示自變量；R 表示定義域；L 表示氣動載荷下風力發電機葉片氣動力[7]；m 表示單位長度葉片質量；e 表示質心偏移；w 表示揮舞位移；GJ 表示扭轉剛度；Rm表示極回轉半徑；ξ 表示扭轉位移；ρ 表示葉片密度；A 表示葉片橫截面積；ψ 表示風輪轉速；Rm1和Rm2表示繞主中心軸和垂直于彈性軸的軸質量回轉半徑；EI 表示揮舞彎曲剛度；P 表示氣動載荷下葉片變槳距[8]；′、·和″、··表示對應參數一階和二階導。

基于西奧道森理論，可得到非定態氣動載荷下葉片氣動力：

式中：VL表示升力系數；ρ∞表示空氣密度；l 表示半弦長；d 表示彈性軸與弦中點距離；v 表示來流速度；C（k）表示西奧道森函數。

離散化葉片結構模型，并采用非保守系統Lagrange方程獲取風力發電機葉片結構動力學方程[9]，如下所示：

式中：q 表示廣義坐標矢量；M、C 和K 表示結構質量、阻尼矩陣和剛度矩陣；QF表示結構外激勵。

通過振型疊加法求解結構動力學方程[10]，基于特征方程Kφi=Mλiφi求解特征值λ1，λ2，…，λM和特征向量φ1，φ2，…，φM，忽略系統高階振型，采用保留的m 個低階振型構建振型矩陣Ω＝［o1，o2，…，om］，由此將式（4）風力發電機葉片結構動力學方程近似為下式形式：

2 風力發電機葉片機械顫振檢測方法

2.1 CFD/CSD 參數耦合算法

采用非保守系統Lagrange 方程，推導出葉片的結構動力學方程，能更準確地描述風力發電機葉片的動力學行為。將非穩態氣動載荷與葉片的結構參數進行耦合求解，得到葉片的結構響應和振動情況。

在式（4）風力發電機葉片結構動力學方程中，右端項QF即為非穩態氣動載荷。將歐拉方程表示為積分方式[11]，如下所示：

為了推進非穩態時間時，引入偽時間步τ，并利用雙時間步長方法求解歐拉方程，獲得葉片的結構響應，具體表達式為

在CFD 模擬中，使用CFD 技術對風流場數值模擬，計算出時間變化的氣流場分布，并獲得每個時間點上的氣動載荷分布[12]。在CSD 模擬中，根據當前時間的氣動載荷和先前時間步的結構響應計算葉片的位移、速度和加速度等參數。利用交錯時間推進方式，將氣動載荷傳遞給結構動力學模型，反饋得到葉片的結構響應。具體流程如圖1 所示。

圖1 CFD/CSD 參數耦合算法流程Fig.1 CFD/CSD parameter coupling algorithm flow chart

為了減小標準CFD/CSD 參數耦合算法誤差，并提高計算精度，需要對該算法加以優化：在風力發電機葉片氣動彈性計算中，非穩態氣動力和結構位移通常表現出連續變化的趨勢。由此對上圖步驟（1）加以改進，基于tm-2、tm-1和tm3 個時刻氣動力實行三點拋物插值，生成tm～tm+1內氣動力關于時間的二次函數，并將其代入式（4）中，可實現對tm+1時刻的結構參數求解。通過這種優化方法，能夠有效提升CFD/CSD 參數耦合算法的計算精度。這樣的優化可以使非穩態氣動力和結構位移在時間上更加連續，從而更準確地反映風力發電機葉片在風流中的響應特性。

2.2 葉片機械顫振邊界計算

為了進一步提升CFD/CSD 參數耦合算法計算效率，提出基于沃爾泰拉級數的非穩態氣動降階模型[13]。沃爾泰拉級數為無窮級數，在系統受到小擾動時，可將非穩態氣動力轉換為二階沃爾泰拉級數形式，其在時域內離散形式如下所示：

式中：u（·）表示輸入值；n 表示離散時間變量；k、k1和k2表示時移；y（n）表示響應值；hi（·）表示系統的i階沃爾泰拉核；h0表示穩態項。

基于階躍響應對非穩態氣動力加以辨識，構建CFD/CSD 參數耦合系統階躍響應如下所示：

式中：δ（n）表示離散域內的單位脈沖信號；ζ0表示單位階躍激勵。

辨識后可得部分非對角分量和全部二階核對角分量的近似一階核（n），如下所示：

經由以上方法可生成非穩態氣動降階模型，將該模型作為廣義位移輸入，輸出廣義力，進而計算QF，并利用模態疊加法計算葉片結構動力學方程，通過調整葉片振型縮放系數，使每個振型的響應達到最大值，最終獲取到葉片機械顫振邊界。

2.3 葉片機械顫振檢測

采用能量法在能量交換的角度分析風力發電機葉片是否存在機械顫振[14-15]，由于風力發電機葉片機械結構氣動阻尼遠遠高于機械阻尼，因此，氣動力做功可作為葉片振動系統自外界獲取到的主要能量，通過檢測非穩態氣動力在一個振動周期內對風力發電機葉片的做功情況可確定是否存在葉片機械顫振。用Tb表示振動周期，（x，t）表示風力發電機葉片表面任意點，P（x，t）表示靜壓，n（x，t）表示單位外法向量，V（x，t）表示速度矢量，則周期Tb內單位面積上非穩態氣動力做功W 如下所示：

阻尼的主要作用是消耗能量，為便于計算，簡化振動系統中非線性阻尼為線形黏性阻尼。在流場繞流情況較好且系統中不存在較強非線性時，可采用黏性阻尼ζ 等效近似替代氣動阻尼，如下所示：

式中：∑W 表示非穩態氣動力做功之和；qcfd表示正則坐標下模態振幅。

基于沃爾泰拉級數的非穩態氣動降階模型和CFD/CSD 參數耦合算法實現非穩態計算后，進一步計算一個振動周期內非穩態氣動力做功和等效模態氣動阻尼比，若非穩態氣動力做功和等效模態氣動阻尼比為正，則風力發電機葉片存在機械顫振情況，需要對其進一步關注分析，以免對風力發電機安全運行產生威脅。

3 實驗與結果

為了驗證參數耦合下風力發電機葉片機械顫振檢測研究整體有效性，需要測試參數耦合下風力發電機葉片機械顫振檢測研究。

為模擬風力發電機葉片機械顫振效果，實驗所用設備及實驗流程如下所示：

（1）實驗設備：鼓風機、調速器、葉片、壓電薄 膜、LabVIEW 數據采集器、NI PXI-1042Q、上位機等。

（2）實驗流程：利用鼓風機生成氣流并通過風洞加速，采用調速器調節風速使其符合實驗要求，在風洞出口0.1 m 處固定葉片，使用壓電薄膜采集葉片信號，經由LabVIEW 數據采集器上傳至上位機，計算葉片機械顫振情況作為實驗實際值。

通過不斷調節各項條件，采集不同實驗環境下的實驗樣本，利用所提方法驗證不同實驗中葉片機械顫振情況，并與實際值加以對比，以評價方法在風力發電機葉片機械顫振檢測中的性能。

通過所提方法獲取不同葉間相角差下氣動阻尼系數，并與實際值加以比較，結果如圖2 所示。

圖2 氣動阻尼系數檢測結果Fig.2 Aerodynamic damping coefficient test results

由圖2 可以看出，所提方法得到的氣動阻尼系數整體更接近于實際值。氣動阻尼系數與風力發電機葉片機械顫振檢測息息相關，準確的氣動阻尼系數是有效識別機械顫振的基礎。

判定不同條件下機械顫振情況是風力發電機葉片機械顫振檢測的主要目的，采用檢測正確率作為指標評估所提方法在顫振孕育狀態和顫振爆發狀態中的檢測表現，結果如表1 所示。

表1 機械顫振檢測結果Tab.1 Mechanical flutter test results

由表1 可以看出，所提方法檢測正確率高達90%以上，適合用于實際檢測之中。

驗證所提方法在風力發電機葉片機械顫振檢測中的響應效率，檢測結果如表2 所示。

表2 計算時間檢測結果Tab.2 Calculate time detection results

由表2 可以看出，所提方法在風力發電機葉片機械顫振檢測中的響應效率明顯較高。所提方法基于沃爾泰拉級數構建非穩態氣動降階模型，進而有效提升了CFD/CSD 參數耦合算法計算效率。

4 結語

風力發電機旋轉葉片結構動力學方程，利用優化后CFD/CSD 參數耦合算法求解流體和結構動力學方程，建立非穩態氣動降階模型求解顫振邊界，計算風力發電機旋轉葉片動態載荷，利用能量法實現風力發電機機械顫振檢測。該方法能夠獲取到更為準確的氣動阻尼系數、檢測正確率更高、計算時間更短，為風力發電機的進一步發展和應用奠定基礎。