山地典型地形下的2 MW風力機仿真研究

■ 龐博鄧勝祥（1.中南大學能源科學與工程學院；2.長沙理工大學可再生能源電力技術(shù)湖南省重點實驗室）

山地典型地形下的2 MW風力機仿真研究

■ 龐博1*鄧勝祥1，2
（1.中南大學能源科學與工程學院；2.長沙理工大學可再生能源電力技術(shù)湖南省重點實驗室）

根據(jù)我國中部某實測地的山地狹管及風速數(shù)據(jù)，使用Solidworks建立2 MW風力機模型，并分別建立了風力機在平地與狹管中的計算域模型，然后通過ICEM劃分網(wǎng)格，用Fluent軟件對模型進行仿真。仿真結(jié)果表明，額定風速下，平地風力機輸出功率比理論值低，與平地相比，狹管對流體有明顯的加速作用，可使狹管風力機比同入流風速下的平地風力機功率高約5.0％ ～43.5％。研究結(jié)果為在山地狹管中布置風力機提供了參考。

復雜山地；CFD仿真；輸出功率；狹管效應(yīng)；水平軸風力機；流場分布；Wilson理論

0　引言

隨著我國風電產(chǎn)業(yè)的高速發(fā)展，風電場的選址地點逐漸從沿海與西部地區(qū)向內(nèi)陸發(fā)展，地形逐漸從平地向山地、丘陵等復雜地形轉(zhuǎn)變，而提高風電場發(fā)電效率的關(guān)鍵是在電場宏觀選址區(qū)域進行科學的微觀選址。山地狹管是山地與丘陵地帶常見的一種地形，當氣流流經(jīng)山地狹管時，由于空氣質(zhì)量不能大量堆積，其將加速流過狹管，可以顯著提高此段風力機的輸出功率[1,2]。

目前國內(nèi)外關(guān)于地形對風力機運行過程的研究主要集中在數(shù)值仿真與實驗?zāi)M兩方面。1)實驗方面：通常集中在單一風力機或純地形研究領(lǐng)域，鮮見有人將這兩方面結(jié)合研究。哈爾濱工業(yè)大學的郭文星[3]對二維山坡、二維山脊、三維軸對稱山丘的典型模型風電場進行了實驗研究，并與CFD仿真結(jié)果進行了對比，肯定了數(shù)值模擬的準確性。2)數(shù)值模擬方面：荷蘭海事研究所的Michel Make等[4]應(yīng)用試驗結(jié)合數(shù)值模擬，對比研究了兩種型號風力機表面的擾流情況，結(jié)果表明數(shù)值模擬能準確計算出風力機運行的情況。

現(xiàn)今，國內(nèi)研究者主要使用兩類國外商業(yè)軟件進行山地風力機研究與選址，一類以丹麥KAMMWAsP、美國Meso Map、加拿大WEST為代表，采用中尺度氣象模式與小尺度模式結(jié)合分析，不適合分析局部細微流場；另一類以丹麥WAsP、英國Windfarm和法國Meteodyn WT為代表，采用標準線性流動模型，但在陡峭地形往往會過高估計此地形的風速[5]。

本文利用Wilson理論建立了山地風力機模型，并利用CFD 方法對平地風力機與狹管風力機進行數(shù)值仿真，研究了兩種地形下風力機在多風速下的功率-速度曲線、額定輸出功率條件下的流場分布情況。結(jié)果表明，CFD模擬出的風輪轉(zhuǎn)矩值與風力機尾流發(fā)展結(jié)果，考慮了山體對風力機的影響、壓力與湍耗散的存在，比理論計算更精確地展示風力機運行的真實情況，而合理的功率損失也驗證了數(shù)值計算的正確性。

1　建立數(shù)學模型

1.1數(shù)學模型

仿真基于穩(wěn)態(tài)不可壓縮三維定常雷諾時均N-S方程(RANS)進行數(shù)值模擬，采用剪切壓力輸運SST k-ω湍流模型使方程組封閉；求解器采用Segregated隱式三維穩(wěn)態(tài)算法、 SIMPLE壓力-速度耦合算法；通用控制方程的離散采用有限容積法；對流項差分采用二階迎風格式[6-8]。過程如下：

不可壓縮流體連續(xù)性方程：

雷諾時均N-S方程(RANS)：

比耗散率ω方程：

湍動能k方程：

式中，ui為速度；τij為粘性切應(yīng)力；F1、F2為混合函數(shù)，F(xiàn)2用來修正F1在自由剪切流中的誤差；y 為網(wǎng)格點到最近壁面的距離；v為分子運動粘度；uj為速度分量；t為時間；P為壓力；σω、 σω2、β*為3個封閉常數(shù)；為分子動力粘度；為湍流動力粘度；Fi為質(zhì)量力；上標“ ′ ”為脈動值。

1.2Wilson葉片建模理論

Wilson理論較于Schnlitz理論，考慮了渦流損失因素；相對于Glauert理論，考慮了葉尖損失與升阻比對氣動性能的影響；且理論本身考慮了軸向與周向誘導因子，是當前最合理的葉片設(shè)計理論。其核心思想就是要使各葉素的風能利用系數(shù)Cpr達到最大[5,9]，如下：

半徑r處風能利用系數(shù)：

半徑r處尖速比：

式中，ar為半徑r處軸向誘導因子，0＆lt;ar＆lt;1；br為半徑r處周向誘導因子，0＆lt;br＆lt;1；Ftr為半徑r處葉尖損失系數(shù)，0＆lt;Ftr＆lt;1； λr為半徑r處的尖速比；Cpr為半徑r處風能利用系數(shù)；βr為半徑r處入流角；ω為葉輪轉(zhuǎn)速；B為葉片數(shù)；R為葉輪半徑；Vc為來流風速。

2　建模與網(wǎng)格劃分

2.1風力機與狹管建模

本文根據(jù)湖南某地一年風塔實測風資源數(shù)據(jù)確定風力機參數(shù)，設(shè)計風速11.5 m/s，風輪額定轉(zhuǎn)速17 r/min，風力機機電效率81％，風能利用系數(shù)Cp=0.4，尖速比為7，輪轂半徑1.5 m。使用Wilson法，通過 Matlab計算出如圖1所示的2 MW風力機沿葉片展向各葉素的原始安裝角與弦長；然后采用6階擬合，得到圖2所示的安裝角與弦長擬合曲線，擬合曲線相關(guān)系數(shù)的平方(R2)分別為0.9988與0.9971。Solidworks風力機葉片建模結(jié)果如圖3所示。

圖1　沿葉片各葉素的原始弦長與安裝角

本文使用Global Mapper提取已公開的測風塔實測地STRM高程數(shù)據(jù)，并劃分實測地狹管區(qū)域的等高線，將其導入Solidworks中，適當優(yōu)化后建立狹管模型。選取狹管漸縮段之后布置風力機，狹管與風力機模型如圖4所示，主要尺寸為入口谷寬408.9 m、谷底寬118.1 m、風力機距入口159.8 m。

圖2　沿葉片各葉素的擬合弦長與安裝角

圖3　葉片模型

圖4　狹管、風力機模型及尺寸

2.2計算域建模與網(wǎng)格劃分

本文分別布置了平地風力機與狹管風力機地形，并根據(jù)最小經(jīng)驗尺寸進行計算域劃分，仿真后確定沒有邊界對模型仿真造成影響[10]。其中，平地計算域如圖4所示，長28D、寬8D、總高5.35D、風力機距計算域前端3D。由于狹管地形尺寸較平地大，為了準確模擬此地形下游遠場尾流的發(fā)展情況，計算域擴大為長60D、寬26.08D、高10D、風力機距計算域前端9D，如圖5所示。

圖4　平地計算域模型及尺寸

圖5　狹管計算域模型及尺寸

為提高計算的準確性與效率，全部模型均采用ICEM劃分網(wǎng)格，綜合考慮分塊合理性與計算效率，除如圖6a中葉輪較小的“Y”型域采用非結(jié)構(gòu)網(wǎng)格，其他計算域均采用結(jié)構(gòu)網(wǎng)格劃分；為提高狹管地形仿真的可靠性，在如圖6b所示的山體間區(qū)域?qū)ｉT劃分了一個域提高網(wǎng)格的密度和質(zhì)量；兩個計算域均采取了合理的漸變網(wǎng)格策略以提高計算的效率[11,12]。

圖6　葉輪與狹管區(qū)域網(wǎng)格劃分

對于平地風力機計算域，本文從約120萬網(wǎng)格開始，對壓力梯度較大部分每次增加10萬～15萬網(wǎng)格，并進行額定風速11.50 m/s的計算，直至最終風輪轉(zhuǎn)矩值無明顯變化，最終選擇1726452個網(wǎng)格的方案。對于類似狹管風力機計算域的，從約200萬網(wǎng)格開始計算比較，風速設(shè)定為11.50 m/s，最終選擇2916723個網(wǎng)格的方案[13]。

3　結(jié)果與討論

3.1速度分析

本文所有云圖均以葉輪輪轂中心為(0,0,0)原點。由于來流風速的擾流、葉輪圓周運動、塔架的干擾等因素，空氣在流經(jīng)風力機時將產(chǎn)生較強的湍流，進而會在尾流中產(chǎn)生擾動，增大這些地方的速度梯度[14]。

圖7從原點后30～2230 m，每隔440 m做一個切片，共6個。切片1可發(fā)現(xiàn)由于塔架的存在，葉輪下部有一塊速度漸變的方形區(qū)域；因為風驅(qū)動風輪轉(zhuǎn)動產(chǎn)生力矩，而槳葉對氣流產(chǎn)生反轉(zhuǎn)矩作用，兩葉片間的扇形區(qū)域速度從前葉片至后葉片逐漸增大，符合風力機旋轉(zhuǎn)尾跡理論；受葉片阻擋的影響，葉片后部的半徑略大于葉輪的“Y型”區(qū)域內(nèi)速度梯度變化較大，區(qū)域直徑約為1.5D～2D；從切片2～6可發(fā)現(xiàn)，隨著距離的增加，低速尾流區(qū)域發(fā)展成了端部為圓的“Y型區(qū)域”，影響范圍擴張速度適中，至風力機下游2230 m的切面6處，速度變化微小且均恢復到10.75 m/s以上，氣流已趨于穩(wěn)定。

圖8為平地計算域風速-Y軸曲線，Y軸穿過輪轂與機艙的中軸線，可看出受風輪影響，風速在風輪前端快速下降，至輪轂表面處降至0 m/s，而在機艙后部風速迅速升至約9 m/s，接著平穩(wěn)上升，至風力機后方1700 m處已恢復至10.50 m/s。

圖8 平地風力機風速-Y軸曲線

圖9為原點上方30 m處Y切面風速云圖?？煽闯?，受風力機與山體間的擾流影響，分布明顯比平地時復雜得多，兩座山體外側(cè)形成了高速尾流帶，兩座山體與風力機后部明顯出現(xiàn)了3條低速尾流帶，同時風力機與兩側(cè)山體間則形成了2條高速尾流帶，類似于射流卷吸效應(yīng)，使風力機后的低速尾流帶在約-220 m處迅速分叉，于約-1000 m處消失。而兩座山體由于外形差異導致尾流差別很大，左側(cè)山體內(nèi)部向狹管凸出、外部有一斷崖，使得流體較難繞到山體后側(cè)，增強了其低速尾流帶，可影響到-800 m后；右側(cè)山體內(nèi)部凹陷、整體呈圓錐臺狀，流體容易繞流，使其低速尾流帶較弱，于約-400 m處消失，但中低速區(qū)域都可影響到-4500 m處。

圖9 狹管計算域Y=30 m切面風速

圖10為狹管計算域風速-Y軸曲線，Y軸穿過輪轂與機艙的中軸線?？煽闯觯L速受狹管影響加速，至風力機前方達到11.80 m/s；接著在風輪前端快速下降；受狹管內(nèi)高風速的影響，在機艙后部風速迅速升至約11.65 m/s；至狹管尾部漸擴段且受尾流影響，風速降至約6.5 m/s；接著由于受兩邊高速尾流帶影響，尾流分叉，風速快速升到約10 m/s直至末端。

圖10　沿Y軸狹管風力機風速

3.2全風速段功率分析

本文對平地風力機與狹管風力機從風力機切入速度3.00 m/s至額定功率速度11.50 m/s分別做了仿真，并與各風速的理論輸出功率值進行了對比，如圖11所示。由于考慮到實際情況并未將葉片設(shè)成光滑表面，且Wilson方法未考慮輪轂損失、尾流損失因子、葉輪受塔架回流等因素的影響，因此平地風力機輸出功率較理論值有一定減小，但變化幅度并不大。而因為狹管對來流明顯的加速能力，使狹管風力機在各個風速條件下的輸出功率都明顯高過平地與理論風力機的功率值，且狹管風力機可在風速較低時就達到風力機額定輸出功率2 MW，從而使風力機在更廣闊的風速段獲得額定發(fā)電功率。

圖11 全風速段3種風力機輸出功率

圖12展示的是理論、狹管及平地風力機三者間的差值?？蓮拇藞D直觀看出3條差值曲線的斜率呈不同程度的增大，且有關(guān)狹管風力機的2條曲線的斜率明顯更大。平地-理論差值升高較慢，最大值位于11.60 m/s附近，在140 kW以下。而有關(guān)狹管風力機的2條差值曲線均在9.60 m/ s附近達到最大值，此處平地-理論差值約為80 kW，而狹管-理論差值＆gt;750 kW，占2 MW理論功率的37.5％，是同風速下平地-理論差值的9.38倍，是平地-理論最大值的5.36倍；另外，此風速狹管-平地差值約為870 kW，為2 MW理論功率的43.5％，是同風速下平地-理論差值的10.88倍，是平地-理論最大值的6.21倍。

圖12　各風速下3種風力機功率差

4　結(jié)論

1)狹管區(qū)域?qū)α鹘?jīng)其內(nèi)的流體有明顯的加速效果，從而可明顯提高風力機的輸出功率。對于本文建立的2 MW風力機，在入口風速約9.6 m/s時，狹管對風力機功率的提升達到最大，比同風速條件下的平地風力機高870 kW，占2 MW理論功率的43.5％，是同風速下平地-理論功率差值的10.88倍，功率提升效果十分明顯。

2)從切入風速3 m/s至額定風速11.5 m/s，狹管對功率提升的效果逐漸增加，使狹管風力機在約9.60 m/s的風速下即可達到設(shè)計11.50 m/s才能達到的2 MW額定功率。

3)Wilson葉片建模法總體來說考慮比較全面，模型輸出功率符合預(yù)期，但由于考慮到實際情況仿真時并未將葉片設(shè)成光滑表面，且Wilson方法未考慮輪轂損失、尾流損失因子、葉輪受塔架回流等因素的影響，使同風速下平地風力機功率比理論值稍低。

4)狹管區(qū)域收縮段對流體加速效果明顯，但風力機后部的整體流場明顯比單機時復雜得多，且影響距離與兩座山體的外形密切相關(guān)，由于擾流等因素將產(chǎn)生多條高速與低速尾流帶，彼此相互影響。在選擇狹管區(qū)域與進行風力機排布時，需充分考慮尾流發(fā)展的情況。

[1] 沈晶, 賴旭. 峽谷地形條件下風電場風況數(shù)值模擬研究[J].水電能源科學, 2011, 29(8): 167-171.

[2] 龐加斌. 沿海和山區(qū)強風特性的觀測分析與風洞模擬研究[D]. 上海: 同濟大學, 2006.

[3] 郭文星. 復雜山地地形風場 CFD 多尺度數(shù)值模擬[D]. 哈爾濱: 哈爾濱工業(yè)大學, 2010.

[4] Michel Make, Guilherme Vaz. Analyzing scaling effects on offsho-re wind turbines using CFD[J]. Renewable Energy, 2015, (83): 1326-1340.

[5] 陳潔鷹. 2MW風力機葉片優(yōu)化設(shè)計及其關(guān)鍵性能仿真研究[D]. 沈陽: 東北大學, 2012.

[6] Carrión M, Steijl R, Woodgate M, et al. Aeroelastic analysis of wind turbines using a tightly coupled CFD–CSD method[J]. Journal of Fluids and Structures, 2014, (50): 392-415.

[7] Alexandros Makridisn, John Chick. Validation of a CFD model of wind turbine wakes with terrain effects[J]. Journal of Wind Engineering and Industrial Aerodynamics, 2013, (123): 12-29.

[8] 李少華, 王東華, 岳巍澎, 等. 雙風力機風向變化時尾流及陣列數(shù)值研究[J]. 動力工程學報, 2011, 31(10): 768-778.

[9] 汪泉. 風力機葉片氣動外形與結(jié)構(gòu)的參數(shù)化耦合設(shè)計理論研究[D]. 重慶: 重慶大學, 2013.

[10] Francesco Balduzzi, Alessandro Bianchini, Riccardo Maleci, et al. Critical issues in the CFD simulation of Darrieus wind turbines[J]. Renewable Energy, 2016, (85): 419-435.

[11] Young-Tae Lee, Hee-Chang Lim. Numerical study of the aerody-namic performance of a 500W Darrieus-type vertical-axis wind turbine[J]. Renewable Energy, 2015, (83): 407-415.

[12] Abdullah Mobin Chowdhury, Hiromichi Akimoto, Yutaka Hara-M. Comparative CFD analysis of vertical axis wind turbine in upright and tilted configuration[J]. Renewable Energy, 2016, (85): 327-337.

[13] Juliana B R Loureiro, Alexandre T P Alho, Atila P Silva Freire. The numerical computation of near wall turbulent flow over a steep hill[J]. Journal of Wind Engineering and Industrial Aerodynamics, 2008, 96(5): 540-561.

[14] Li Y, Castro A M, Sinokrot T, et al. Coupled multi-body dynamics and CFD for wind turbine simulation including explicit wind turbulence[J]. Renewable Energy, 2015, (76): 338-361.

2016-03-11

可再生能源電力技術(shù)湖南省重點實驗室(長沙理工大學)開放基金資助項目(2012ZNDL008）

龐博( 1991—)，男，碩士研究生，主要從事風能資源評估，風力機發(fā)電、計算機仿真與優(yōu)化，熱工設(shè)備檢測與熱工計算方面的研究。cspangbo@163.com