適用于獨立變槳控制方式的風機氣動模型和風模型的研究

張 迪，王維慶，王海云，李 強

（1.新疆大學電氣工程學院，新疆 烏魯木齊 830047；2.北京金風科創風電設備有限公司，北京 100176）

0 引言

風機變槳控制分為集體變槳和獨立變槳，隨著風機的容量，葉輪和塔架的尺寸的增加，如何降低風機載荷的問題日益凸顯出來，獨立變槳同集體變槳相比的優勢在于降低風機載荷同時保證風力發電機的輸出額定功率。獨立變槳是變槳控制發展的趨勢。少數國外先進的風機制造商已經將獨立變槳應用于兆瓦級風機上，國內先進的風機制造商也加大了對獨立變槳控制技術的研究。

隨著風機變槳控制方式從集體變槳過渡到獨立變槳，由于在傳統的應用于3葉片集體變槳控制策略的風機氣動模型中，3只槳葉的槳距角對風機氣動轉矩的影響是耦合在一起的， 無法利用模型來分析單只葉片的槳距角對單只葉片的拍擊力和升力的影響，所以需要回到葉素理論中來研究適用于獨立變槳風機氣動模型的建立方法，因此需要一種新的建模方式。

1 風機氣動模型的建立

1.1 應用于集體變槳的風機氣動模型

[1]中的風機風能利用系數Cp的經驗計算公式為：

式中，λ為葉尖速比；λi為中間變量；β為槳距角；C1=0.5176， C2=116， C3=0.4， C4=5， C5=21， C6=0.006 8。

風輪吸收的功率表達式為：

式中，ρ為空氣密度；R為葉片半徑；V為風速。風輪產生的轉矩為：

式中，ω為風輪轉速。

葉尖速比λ的計算公式：

將式（4）和（5）代入式（3）中得出風機氣動轉矩 T 的表達式為：

基于式 （1）、（2）、（6） 在 simulink 中搭建模型結構，仿真條件為槳距角β取0°，空氣密度為1.29 kg/m3，葉片半徑40 m；風輪轉速ω從0 r/min變化到35 r/min， 風速從3 m/s到14 m/s。

仿真結果：①風機的輸出轉矩T變化情況如圖1中曲線所示；②根據文獻[1]，Cp最大值對應的減速比λopt為 8.1，帶入式（5）可以計算出最優轉速ωopt和與其對應的最優轉矩Topt。用光滑的曲線連接形成圖1中虛線。虛線中的低風速下段可以作為低風速下最大功率跟蹤的T-ω參考曲線。

圖1 槳距角為0°，不同風速下葉輪轉矩隨轉速變化曲線

當風速高于額定風速時，通過變槳電機驅動槳葉增大槳距角，降低了Cp，從而降低風機氣動轉矩，維持風輪轉速額定。圖2繪制了風速為12 m/s、槳距角為0°，風速為13 m/s、槳距角為1°，風速14 m/s、槳距角5°轉矩同轉速變化的對比曲線，可以看出增大槳距角對降低氣動轉矩的效果是非常明顯的。

圖2 變槳對氣動轉矩的影響

這種建模的方式可以用于對集體變槳策略的仿真驗證，但是對于獨立變槳，由于模型沒有獨立的槳距角輸入，而只能輸入3個槳葉集體的槳距角，所以無法應用于獨立變槳策略的驗證。

1.2 應用于獨立變槳的風機氣動模型

根據文獻[2]的葉素理論如圖3所示。若以葉素為參考系，風將會疊加一個同葉素轉動方向相反，大小相等的線速度-u，風速V和-u矢量合為葉素的相對風速w。w為葉素真正感受到的風速。由于相對風速w的作用，在風機葉素產生了與w同一方向的阻力Fdrag與垂直于阻力的升力Flift。

圖3 葉素受力分析

式中，Cd為阻力系數，Cl為升力系數，同Cp一樣沒有量綱，二者的取值變化決定于攻角i，影響Cd、Cl的因素有彎曲度，葉片厚度，表面粗糙度和雷諾數，Cd、Cl數據參考翼型為NACA0015；L為弦長；dr可以理解為一段很短的葉展。攻角變化范圍為0-2π。

為了計算出作用在風機主軸上的氣動轉矩Tair和作用在葉片上的拍擊力Fflap，首先需要將葉素上的阻力和推力經行軸向和切向分解與合成 （如圖4）。

每一小段葉素上產生的升力dFu和拍擊力dFa由公式（10）和公式（11）計算出來

圖4 Fdrag和Flift的分解與合成

式中，I為入流角。

為了實現仿真，將微分放大為變化量。將每只槳葉被分為均勻的10段，每段的弦長為此段兩端弦長的平均值，分別建立每段葉片的拍擊力△Fa和升力△Fu的模型。

將一只葉片上每段的拍擊力△Fai與對輪轂處的轉矩△Fuiri求和，最終計算出每只葉片上的拍擊力F（flap）j與 3只葉片總體在主軸上的氣動轉矩 Tair。

式中，i為每段葉片的標號；ri為第i段葉片中心距離輪轂的距離；j為葉片數。

基于式（12）～（17）在 simulink 中建立仿真模型，仿真條件為：空氣密度1.25 kg/m3， 葉片長度40 m，風速從3 m/s至13 m/s，風輪轉速從0到35 r/min，槳距角保持-1°。

仿真結果如圖5，主軸氣動轉矩隨著風速，風輪轉速變化曲線。由于槳距角為-1°，所以整體的主軸氣動轉矩大于圖1中的主軸氣動轉矩。

圖5中橫坐標為主軸轉速，縱坐標為氣動轉矩，曲線下方的面積為葉輪吸收的風能，根據最大面積的方法同樣可以得到最佳轉速值以及對應的最佳轉矩值作為低風速下的最大功率跟蹤T-ω參考曲線。

1.3 兩種建模方式的比較

（1）第一種風機氣動模型是否可以反映風機的氣動轉矩關鍵在于對Cp（λ，β）曲線的擬合準確度，對于不同類型的葉片需要大量的數據進行擬合，而第二種氣動模型的建立只需要知道葉片的升力系數Cl和阻力系數Cd，這兩個參數可以在對應的葉片參數中找到，能夠直接放入模型中，而且由于單只葉片模型由一段段葉片拼接而成，由此就能實現對單獨一只葉片的槳距角和風速進行獨立修改，可應用于獨立變槳中。

（2）第一種風機氣動模型無法計算單只葉片的拍擊力而第二種模型可以計算。

（3）第二種模型的分段數可以增加，將會更加準確模擬葉片，具有發展潛力。

圖5 新模型不同風速下葉輪轉矩隨轉速變化曲線

2 風剪塔影對風機氣動模型輸出的影響

2.1 風模型同風機氣動模型的影響與配合

風剪塔影最終對于風機氣動模型模型的影響包括：①在單只葉片上會產生1P的拍擊力波動，波動頻率為葉輪旋轉周期的倒數。②在風機主軸轉矩中產生3P的波動，波動頻率為1P波動的頻率的三倍。塔影引起的波動要比風剪更大。

由于風機的每只葉片被分為10段，風速模型輸出的風速也與這10段每一段中心的位置對應。即風速最終都應化為葉片位置角度θ與葉片位置點到輪轂中心距離ri的函數從而風速模型輸出30個風速變量與3只葉片中的30個葉片段位置對應。

2.2 風剪切以及與模型配合

風剪切是指當隨著高度的不同，風速隨之變化從而造成轉矩和拍擊力的周期性波動，參考文獻[3]中的風剪切公式：

式中，V（h）為某一高度處風速；V（h0）為輪轂高度處風速；α為風剪切經驗指數，取1/7。

將式（18）轉化為每段葉片位置角度的函數如式（19）：

式中，ri為第i段葉片距離輪轂的半徑長度；ω為風輪轉速；根據式（20）和式（21），作用在單只葉片某段上的風速，會疊加一個頻率為 （ω/60）Hz的分量。仿真設定條件為：葉片半徑40 m，輪轂高度70 m。轉速為19 r/min，槳距角為4°，輪轂處的風速為13 m/s，觀察與一段葉片對應的風剪模型的輸出風速如圖6。

圖6 風剪切效應下的風速波動

對風速模型輸出的數據經行快速傅里葉變化得到頻譜圖如圖7。

圖7 風速頻譜圖

由圖7可以看出，風速中除了直流分量外，還有明顯疊加了0.2 Hz諧波含量。

風剪切造成葉片拍擊力的1 P波動如圖8，頻率為0.2 Hz，即等于葉輪旋轉周期5 s的倒數，而仿真設定轉速為12 r/min，也就是0.2 r/s， 即5 s/r，即葉輪旋轉周期周期為5 s仿真結果驗證了風剪對單葉片拍擊力和氣動轉矩影響的結論。

圖8 風剪切造成葉片拍擊力的1P波動

風剪切造成主軸轉矩3P波動，頻率為0.6Hz如圖9。

2.3 塔影效應及與模型配合

塔影效應指塔架正面的風速由于塔架的作用而減小，所以當葉片經過塔架的豎直位置時，葉片上的拍擊力會減小，單只葉片旋轉一周后由風速引起的拍擊力會產生1P波動。而對于3葉片風機，葉輪旋轉一周后，每只葉片都經過一次塔架位置，所以在主軸轉矩中會產生3P波動。結合文獻[3]的塔影效應計算公式得出每段葉片對應位置的風速表達式：

圖9 風剪切造成主軸轉矩3P波動

式中，a為塔架的半徑 ；R為葉片半徑；H為塔架高度；x為葉片與塔架的距離；塔影效應只在位置角度θ為π到π內起作用。仿真的其他條件不變考慮塔影效應的情況下觀察到單只葉片的拍擊中產生了1P波動，主軸轉矩中產生3P波動，塔影效應造成葉片拍擊力的1P波動，頻率為0.2 Hz，如圖10。

圖10 塔影效應造成葉片拍擊力的1P波動

塔影效應造成主軸轉矩3P波動，頻率為0.6 Hz如圖11。

圖11 塔影效應造成主軸轉矩3P波動

2.4 風剪塔影對風機的影響

風剪塔影都會造成單只葉片的拍擊力中的1P波動，主軸轉矩中的3P波動，但是由于塔影效應對風速的影響比風剪切更加嚴重。所以塔影效應產生的波動更加嚴重。

3 獨立變槳對葉片拍擊力的影響

為了減小風對葉片拍擊力的影響，需要對每只葉片的槳距角進行獨立的控制，風速為13 m/s，轉速為12 r/min的條件下，3只槳葉的槳距角分別為2°，4°，6°，從仿真結果可以看出槳距角的增大對葉片的拍擊力的影響如圖12。

圖12 獨立變槳對葉片拍擊力的影響

4 結論

（1）綜合比較了基于Cp經驗公式與基于葉素理論的風機氣動模型，通過分析兩種模型構架思想和仿真結果，說明了基于葉素理論的風機氣動模型應用于獨立變槳的正確性。

（2）在基于葉素理論的風機氣動模型與基于風剪切塔影效應理論的風模型基礎上，通過分析二者聯合仿真的結果，驗證了風剪塔影導致單只葉片拍擊力的1P波動和3只葉片產生主軸轉矩中的3P波動。

參考文獻：

［1］ HeierS.GridIntegrationofWindEnergyConversionSystems［M］.1998.

［2］ 葉杭冶.風力發電機組的控制技術［M］.機械工業出版社2002:53-33.

［3］ 李少林，張興，謝震，等.雙饋風力發電系統3次功率脈動的研究［J］.電網技術 2010， 34（4）:38-39.

［4］ Song S H,Jeong B C,Oh J H,et al.Simulation of output characteristics of rotor blades using a hardware-in-loop wind turbine simulator［J］.IEEE Trans on Energy Conversion,2005,3（6）:1791-1796.

［5］ Anderson C G,Richon J-B,Campbell T J.An Aerodynamic Moment-Controlled Surface for Gust Load Alleviation on Wind Turbine Rotors［J］.IEEE TRANSACTIONS ON CONTROL SYSTEMS TECHNOLOGY,1998,6（5）:578-580.

［6］ Dolan D S L,Lehn P W.Simulation model of wind turbine 3p torque oscillations due to wind shear and tower shadow［J］.IEEE Trans on Energy Conversion， 2006， 21（3）:717-724.

［7］ Geyler M,Caselitz P.Individual Blade Pitch Control Design for Load Reduction on Large Wind Turbines［C］//Proceedings of European Wind Energy Conference,2007:1-4.

［8］ 董萍，吳捷，楊金明，等.風力發電機組建模研究現狀［J］.太陽能學報 2004， 25（5）:613-615.

［9］ Van Baars G E,Bongers P M M.Flexible wind turbine model validation［J］.Wind Engineering,1992,16（4）:247-256.

［10］ Muljadi E,Butterfield C P.Pitch-Controlled Variable-Speed Wind Turbine Generation［J］.IEEE Trans on Industry Applications,2001,37（1）:240-246.