張 純，孫 山，韓煜航

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

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

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

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

1 風(fēng)力發(fā)電機葉片結(jié)構(gòu)動力學(xué)方程推導(dǎo)

基于懸臂旋轉(zhuǎn)歐拉－伯努利梁運動微分方程和哈密頓原理建立風(fēng)力發(fā)電機旋轉(zhuǎn)葉片結(jié)構(gòu)模型[5-6]，如下所示：

式中：U 表示應(yīng)變能；T 表示動能；W 表示虛功；δ 表示變分。

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

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

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

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

離散化葉片結(jié)構(gòu)模型，并采用非保守系統(tǒng)Lagrange方程獲取風(fēng)力發(fā)電機葉片結(jié)構(gòu)動力學(xué)方程[9]，如下所示：

式中：q 表示廣義坐標(biāo)矢量；M、C 和K 表示結(jié)構(gòu)質(zhì)量、阻尼矩陣和剛度矩陣；QF表示結(jié)構(gòu)外激勵。

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

2 風(fēng)力發(fā)電機葉片機械顫振檢測方法

2.1 CFD/CSD 參數(shù)耦合算法

采用非保守系統(tǒng)Lagrange 方程，推導(dǎo)出葉片的結(jié)構(gòu)動力學(xué)方程，能更準(zhǔn)確地描述風(fēng)力發(fā)電機葉片的動力學(xué)行為。將非穩(wěn)態(tài)氣動載荷與葉片的結(jié)構(gòu)參數(shù)進行耦合求解，得到葉片的結(jié)構(gòu)響應(yīng)和振動情況。

在式（4）風(fēng)力發(fā)電機葉片結(jié)構(gòu)動力學(xué)方程中，右端項QF即為非穩(wěn)態(tài)氣動載荷。將歐拉方程表示為積分方式[11]，如下所示：

為了推進非穩(wěn)態(tài)時間時，引入偽時間步τ，并利用雙時間步長方法求解歐拉方程，獲得葉片的結(jié)構(gòu)響應(yīng)，具體表達式為

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

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

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

2.2 葉片機械顫振邊界計算

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

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

基于階躍響應(yīng)對非穩(wěn)態(tài)氣動力加以辨識，構(gòu)建CFD/CSD 參數(shù)耦合系統(tǒng)階躍響應(yīng)如下所示：

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

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

經(jīng)由以上方法可生成非穩(wěn)態(tài)氣動降階模型，將該模型作為廣義位移輸入，輸出廣義力，進而計算QF，并利用模態(tài)疊加法計算葉片結(jié)構(gòu)動力學(xué)方程，通過調(diào)整葉片振型縮放系數(shù)，使每個振型的響應(yīng)達到最大值，最終獲取到葉片機械顫振邊界。

2.3 葉片機械顫振檢測

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

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

式中：∑W 表示非穩(wěn)態(tài)氣動力做功之和；qcfd表示正則坐標(biāo)下模態(tài)振幅。

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

3 實驗與結(jié)果

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

為模擬風(fēng)力發(fā)電機葉片機械顫振效果，實驗所用設(shè)備及實驗流程如下所示：

（1）實驗設(shè)備：鼓風(fēng)機、調(diào)速器、葉片、壓電薄 膜、LabVIEW 數(shù)據(jù)采集器、NI PXI-1042Q、上位機等。

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

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

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

圖2 氣動阻尼系數(shù)檢測結(jié)果Fig.2 Aerodynamic damping coefficient test results

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

判定不同條件下機械顫振情況是風(fēng)力發(fā)電機葉片機械顫振檢測的主要目的，采用檢測正確率作為指標(biāo)評估所提方法在顫振孕育狀態(tài)和顫振爆發(fā)狀態(tài)中的檢測表現(xiàn)，結(jié)果如表1 所示。

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

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

驗證所提方法在風(fēng)力發(fā)電機葉片機械顫振檢測中的響應(yīng)效率，檢測結(jié)果如表2 所示。

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

由表2 可以看出，所提方法在風(fēng)力發(fā)電機葉片機械顫振檢測中的響應(yīng)效率明顯較高。所提方法基于沃爾泰拉級數(shù)構(gòu)建非穩(wěn)態(tài)氣動降階模型，進而有效提升了CFD/CSD 參數(shù)耦合算法計算效率。

4 結(jié)語

風(fēng)力發(fā)電機旋轉(zhuǎn)葉片結(jié)構(gòu)動力學(xué)方程，利用優(yōu)化后CFD/CSD 參數(shù)耦合算法求解流體和結(jié)構(gòu)動力學(xué)方程，建立非穩(wěn)態(tài)氣動降階模型求解顫振邊界，計算風(fēng)力發(fā)電機旋轉(zhuǎn)葉片動態(tài)載荷，利用能量法實現(xiàn)風(fēng)力發(fā)電機機械顫振檢測。該方法能夠獲取到更為準(zhǔn)確的氣動阻尼系數(shù)、檢測正確率更高、計算時間更短，為風(fēng)力發(fā)電機的進一步發(fā)展和應(yīng)用奠定基礎(chǔ)。