5 MW風力機塔影效應雙向流固耦合分析

李 麗， 潘 揚， 胡丹梅

(1.國網湖南省電力有限公司防災減災中心(電網輸變電設備防災減災國家重點實驗室)，長沙 410129；2.上海電力學院 能源與機械工程學院，上海 200090)

近年來，由于風力機單機容量不斷增大，風力機葉輪的掃風面積和工作高度也隨之提高，導致作為支撐結構的塔筒規格不斷增大，塔筒引起的塔影效應也更明顯。塔影效應是指塔架對流場產生干擾，而流場將此干擾傳遞到風力機葉片上，導致風力機的氣動特性發生變化。由于風力機塔架對空氣有阻擋和排擠作用，使得其附近的流場發生變化，空氣作用于葉片上的氣動載荷發生劇烈波動，造成葉片的形變與其應力也會發生相應變化。

對于目前的兆瓦級大型風力機而言，塔影效應已經成為一個不可忽視的影響因素。研究表明，塔影效應使單個葉片載荷波動20%～40%,風輪載荷波動6%～12%[1]，說明塔架會對風力機組的運行產生影響，因此塔影效應具有一定的研究價值。

1968年，研究人員就提出了風力機塔影效應的概念[2]。進入21世紀后，隨著風電產業的興起，風電從業者開始認識到塔架對風力機的影響不可忽略。Das等[3]建立了關于塔影效應和分切效應的時域模型，用于計算塔架造成的風力機周期波動。Sescu等[4]利用商業CFD軟件對塔影效應進行了仿真計算和分析。李少林等[5]在設計動態風力機模擬器時對塔影效應進行了模擬。

流固耦合是指流體與固體之間發生相互干涉，如果這種干涉既包括流體對固體的干涉，又包括固體對流體的干涉，則稱為雙向流固耦合。國內外學者針對風力機雙向流固耦合情況下的工作狀態進行了研究。Hsu等[6]和Bazilevs等[7-8]利用2種方法對NREL 5 MW海上風力機[9]進行了單根葉片和整機的雙向流固耦合數值計算，研究內容包括風力機固體結構和流場。潘旭[10]利用Ansys軟件對兆瓦級風力機葉片進行了雙向流固耦合數值模擬，分析了蒙皮材料與葉厚等因素對葉片工作表現的影響。何福添[11]對3 MW風力機葉片進行了雙向流固耦合數值模擬，并根據模擬結果進行了葉片設計優化。李媛[12]對二維翼型和三維葉片進行了流固耦合，并研究了不同風速條件下風力機葉片的流固耦合特性。目前，研究人員在研究風力機雙向流固耦合時往往僅關注單根葉片或單獨葉輪，對風力機整機塔影效應的研究相對較少。

筆者以NREL 5 MW風力機[9]為模型，通過數據交互模塊System Coupling將流體模擬軟件Fluent與固體分析軟件Transient Structural進行聯動，實現風力機固體結構與其周圍流場的雙向耦合，展開了關于塔影效應的流固耦合研究。

1 計算模型與網格劃分

1.1 物理模型與網格

風力機模型葉輪直徑為126 m，額定風速為11.4 m/s，額定轉速為12.1 r/min，將轉速取整為12 r/min，葉輪旋轉周期為5 s。塔筒為圓臺狀，上底面直徑為3.87 m，下底面直徑為6 m，高為87.6 m。首先利用3D建模軟件ProE建立風力機整機模型和立方體流場域模型，再根據布爾運算得到風力機的計算流場。流場分為旋轉域和靜止域，旋轉域是半徑為70 m、高為8 m的圓柱體，葉輪位于旋轉域中央，被其完全包裹；靜止域為長方體，截面高為400 m，寬為600 m，入口距葉輪旋轉平面為125 m，出口距葉輪旋轉平面為500 m，機艙與塔筒位于靜止域內，流場域整體模型如圖1所示。通過二次開發端口將模型直接導入Ansys Workbench，對風力機葉輪進行厚度賦值，使其成為厚度為8 cm的殼結構，將機艙與塔筒視為剛體。

1.1.2 劃分計算網格

利用流固耦合的Fluent流體模塊(FFF)自帶的Mesh功能對流場進行網格劃分，為能啟用Fluent中動網格的彈簧光順功能和網格重構功能，采用四面體網格，并對固體表面進行加密，網格單元個數為259萬，最大網格偏斜小于0.85，滿足計算要求。利用Workbench中的Mechanical Model對風力機葉輪進行殼單元網格劃分，生成1.6萬個六面體網格，如圖2所示。葉片采用45°鋪層各向異性玻璃纖維環氧樹脂復合材料，其材料力學性能如表1所示，其中E1、E2分別為垂直纖維方向與沿纖維方向的彈性模量，G12為剪切模量，ρ為密度，σ12為泊松比[13]。

圖2 葉輪殼網格

E1/GPaE2/GPaG12/GPaσ12ρ/(g·cm-3)398.63.80.282.1

1.1.3 可信度驗證

為驗證計算的準確性，分別利用627萬、190萬結構化網格和415萬、259萬非結構化網格在非流固耦合情況下進行計算，所得功率偏差分別為0.28%、3.6%、3.2%和5.3%。由于流固耦合計算非常占用計算機資源，為節省計算時間和提高計算機工作穩定性，最終采用259萬非結構化四面體網格，計算誤差約為5%，滿足計算要求。

1.2 計算方法與耦合方程

由于NREL 5 MW風力機的幾何尺寸較大，為了使網格數量保持在計算機可計算范圍內，需對網格的最小尺寸進行控制，這導致固體表面的網格尺寸難以模擬出邊界層。因此，湍流模型選用帶有湍流漩渦修正和低雷諾數修正的RNGk-ε兩方程模型，以達到利用標準壁面函數[14]處理邊界層的目的。研究表明，k-ε模型對于流動的預測具有一定準確性[15-16]。采用滑移網格法進行流體計算，設定旋轉域轉速為12 r/min，時間步長為0.1 s，流場入口風速為11.4 m/s，出口為自由出口，離散格式為二階迎風，利用Simple算法進行求解。對于葉輪固體部分，通過添加約束限制輪轂自由度，只保留旋轉方向自由度，并在葉輪上加載角速度與標準重力，將整個葉輪外表面設為流固耦合面。

Ansys-Fluent的流固耦合方程如下[11]:

(1)

式中：A為系數矩陣，其中下標F代表流體域，S代表固體域，I代表流固交界面；U為流場中的速度；p為流場中的壓力；δ為葉片位移；R為殘差。

利用SPSS軟件進行描述性統計分析；利用Aquachem水化學軟件繪制Piper三線圖進行水化學類型分析；選取超標較嚴重的硫酸鹽(SO2-4)、亞硝酸鹽(NO-2)、硝酸鹽(NO-3)、總硬度和溶解性總固體(TDS)共5項評價指標，利用熵權密切值法對地下水水質進行評價。評價的標準采用《地下水質量標準》(GB/T14848—93)將地下水質量劃分為5類(Ⅰ—Ⅴ)，見表1。

如圖3所示，利用Workbench中的System Coupling連接Fluent模塊和Transient Structural，并在System Coupling中設置Fluent的計算順序為1，Transient Structural的計算順序為2，使流體先作用于固體，固體再反作用于流體。將固體耦合面和流體耦合面配對，設置時間步長為0.1 s，進行耦合計算。

圖3 Workbench模塊結構圖

2 計算結果分析

為研究塔影效應對葉輪的影響，筆者在相同的工況條件下對有塔筒和無塔筒這2種情況進行計算，并進行對比分析。

2.1 風力機葉輪氣動載荷

圖4和圖5分別為風力機葉輪在25～30 s內轉矩和所受軸向推力的變化曲線。在25～30 s內葉輪旋轉了1周，轉矩和軸向推力均出現3次峰值和3次谷值，驗證了風力機的3P閃變[17]。由圖4和圖5可知，在有塔筒的情況下2種曲線的波動程度均比無塔筒時更劇烈。在風力機工況基本穩定的時間段內(25～30 s)，在有塔筒的情況下轉矩峰值約為4 030 kN·m，谷值約為3 800 kN·m，峰谷值之差為230 kN·m；無塔筒時轉矩輸出的峰值與谷值分別為4 000 kN·m和3 930 kN·m，峰谷值之差為70 kN·m。對于葉輪所受軸向推力，在有塔筒的情況下軸向推力的峰值和谷值分別約為787 kN和743 kN，最大與最小值之差為44 kN；在無塔筒的情況下軸向推力的峰值和谷值分別約為779 kN和759 kN，最大與最小值之差為20 kN。綜上，塔筒對風力機葉輪所受氣動載荷有明顯影響。以本計算為例，在塔影效應影響下葉輪所受轉矩和軸向推力的波動幅度是未受塔影效應影響時的3倍以上，說明對于兆瓦級大型風力機，塔影效應的影響不可忽視。

圖4 葉輪轉矩曲線

圖5 葉輪軸向推力曲線

2.2 葉輪形變與應力

圖6給出了模擬時間t為10 s、20 s和30 s時風力機葉輪的形變放大圖，放大倍數約為19倍。由于葉輪旋轉周期T為5 s，在時間為5 s的整數倍時，z方向的葉片位于塔筒正前方，如圖6(b)中葉片1所示。由圖6可以看出，在10 s時風力機葉輪工況尚未穩定，葉片的形變較大，位于塔筒前方葉片的軸向揮舞形變和周向擺振形變均小于其他2根葉片；30 s時葉輪工況基本穩定，3根葉片的擺振形變程度較為接近，但揮舞形變仍有一定差異，此現象與文獻[18]相符。

(a) 10 s時葉片軸向揮舞形變圖

(b) 10 s時葉片周向擺振形變圖

(c) 20 s時葉片軸向揮舞形變圖

(d) 20 s時葉片周向擺振形變圖

(e) 30 s時葉片軸向揮舞形變圖

(f) 30 s時葉片周向擺振形變圖

為直觀反映塔影效應對葉片形變的影響，表2給出了模擬時間t為10 s、20 s和30 s時在有塔筒和無塔筒情況下3根葉片尖部的位移。無塔筒時3根葉片的葉尖位移基本保持一致，偏差不超過5%，在有塔筒干涉時塔筒正前方葉片葉尖的位移明顯小于其余2根葉片，且另外2根葉片間也存在較大差異。

表2 葉片尖部位移

圖7和圖8分別是風力機葉輪的最大形變和最大等效應力曲線，在有塔筒時葉輪的形變和應力的波動程度和數值均高于無塔筒時。與無塔筒時相比，有塔筒時葉輪處于較不穩定的狀態。當經過塔筒時，葉片都會受到一次塔筒干涉，這相當于在風力機運行過程中對葉輪施加了一個周期性的氣動干擾，使葉輪運行的受力情況更為惡劣，這說明了塔影效應會對風力機組運行的穩定性產生負面影響。

圖7 葉輪最大形變曲線

圖8 葉輪最大等效應力曲線

2.3 流固耦合與非流固耦合結果對比

圖9為流固耦合和非流固耦合下葉高為0.6R處葉片附近的流線圖。由圖9可以看出，流固耦合下流線的總體趨勢與非流固耦合情況相似，但在葉片表面附近的流線變化較大。無塔筒時，在流固耦合下葉片前緣和后緣處的流線均出現明顯彎折和回流現象；有塔筒時，葉片壓力面甚至出現了渦旋，說明流固耦合下的塔影效應比非耦合情況更為明顯。

(a) 非耦合，無塔筒

(b) 耦合，無塔筒

(c) 非耦合，有塔筒

(d) 耦合，有塔筒

無論有無塔筒，在流固耦合下葉片對流場的擾動均大于非耦合情況，造成此現象的原因可能是葉片的形變與位移使得葉片表面壓力發生變化。為研究葉片表面的壓力變化，圖10給出了在耦合和非耦合情況下葉片表面的壓力分布。由圖10可知，無論是壓力面還是吸力面，在耦合情況下葉片表面壓力均略低于非耦合情況下，且在耦合情況下壓力面的壓力分布情況比非耦合情況下更復雜，這說明形變對葉片的氣動受力具有明顯影響。

3 結 論

(1)在靠近塔筒的過程中葉片的氣動載荷會發生較大突變，以本計算為例，葉輪受到的周向轉矩與軸向推力的波動幅度是無塔筒時的3倍以上。

(2)有塔筒時葉片的形變程度大于無塔筒時，且有塔筒時3根葉片的葉尖位移差異較大，在塔筒正前方葉片的葉尖位移明顯小于其他2根葉片。

(3)在風力機運行過程中塔筒會不斷對葉輪產生氣動干擾，導致整個葉片的形變和應力要劣于無塔筒時，塔影效應對風力機運行的穩定性有不利影響。

(4)考慮流固耦合的計算更符合風力機的實際工作條件。與非流固耦合相比，在流固耦合下葉片對流場的擾動較大，葉片附近流場流線的渦旋現象更明顯，且考慮流固耦合時葉片表面的壓力分布更復雜。

(5)對于減小塔影效應的對策，現在較為主流的解決方法是獨立變槳距控制，但學者們提出的控制策略大多為反饋控制，其變槳具有一定滯后性。為了探尋減弱塔影效應對風力機影響的方法，筆者曾嘗試將塔筒替換為桁架并進行了分析，結果表明桁架式風力機的氣動載荷波動程度小于塔筒式風力機。采用桁架可在一定程度上提高機組運行的穩定性與安全性，但由于桁架結構復雜，安裝和維護成本較高，其實用性仍需要進一步探討。

(a) 非耦合，壓力面

(b) 耦合，壓力面

(c) 非耦合，吸力面

(d) 耦合，吸力面

：

[1] 范忠瑤, 康順, 趙萍. 上風向風力機塔影效應的數值模擬研究[J].工程熱物理學報, 2012, 33(10): 1707-1710.

FAN Zhongyao, KANG Shun, ZHAO Ping. Numerical simulations of upwind wind turbine tower shadow effects[J].JournalofEngineeringThermophysics, 2012, 33(10): 1707-1710.

[2] CERMAK J E, HORN J D. Tower shadow effect[J].JournalofGeophysicalResearch, 1968, 73(6): 1869-1876.

[3] DAS S, KARNIK N, SANTOSO S. Time-domain modeling of tower shadow and wind shear in wind turbines[J].ISRNRenewableEnergy, 2011, 2011: 1-11.

[4] SESCU A, ANDERSEN B, AFJEH A A, et al. Computational investigation of tower shadow effects on wind turbines[C]//ASME2011InternationalMechanicalEngineeringCongressandExposition. Denver, Colorado, USA:American Society of Mechanical Engineers, 2011.

[5] 李少林, 張興, 楊淑英, 等. 風力機動態模擬器的仿真研究[J].太陽能學報, 2010, 31(10): 1366-1372.

LI Shaolin, ZHANG Xing, YANG Shuying, et al. Simulation research of dynamic wind turbine simulator[J].ActaEnergiaeSolarisSinica, 2010, 31(10): 1366-1372.

[6] HSU M C, BAZILEVS Y. Fluid-structure interaction modeling of wind turbines: simulating the fullmachine[J].ComputationalMechanics,2012, 50(6): 821-833.

[7] BAZILEVS Y, HSU M C, TAKIZAWA K, et al. ALE-VMS and ST-VMS methods for computer modeling of wind-turbine rotor aerodynamics and fluid-structure interaction[J].MathematicalModelsandMethodsinAppliedSciences, 2012, 22(S2): 1-61.

[8] BAZILEVS Y, TAKIZAWA K, TEZDUYAR T E, et al. Aerodynamic and FSI analysis of wind turbines with the ALE-VMS and ST-VMS methods[J].ArchivesofComputationalMethodsinEngineering, 2014, 21(4): 359-398.

[9] JONKMAN J, BUTTERFIELD S, MUSIAL W, et al. Definition of a 5-MW reference wind turbine for offshore system development[R]. Colorado, USA: National Renewable Energy Laboratory, 2009.

[10] 潘旭. MW級風力發電機風輪葉片流固耦合場強度分析[D]. 鄭州: 鄭州大學, 2011.

[11] 何福添. 水平軸風力機葉片的設計及流固耦合計算[D]. 武漢: 武漢理工大學, 2013.

[12] 李媛. 風力機葉片流固耦合數值模擬[D]. 北京: 華北電力大學, 2013.

[13] 胡丹梅, 張志超, 孫凱, 等. 風力機葉片流固耦合計算分析[J].中國電機工程學報, 2013, 33(17): 98-104.

HU Danmei, ZHANG Zhichao, SUN Kai, et al. Computational analysis of wind turbine blades based on fluid-structure interaction[J].ProceedingsoftheCSEE, 2013, 33(17): 98-104.

[14] 方平治, 顧明, 談建國, 等. 數值模擬大氣邊界層中解決壁面函數問題方法研究[J].振動與沖擊, 2015, 34(2): 85-90.

FANG Pingzhi, GU Ming, TAN Jianguo, et al. Method to solve the wall function problem in simulation of atmospheric boundary layer[J].JournalofVibrationandShock, 2015, 34(2): 85-90.

[15] 金新陽, 楊偉, 金海, 等. 數值風工程應用中湍流模型的比較研究[J].建筑科學, 2006, 22(5): 1-5, 28.

JIN Xinyang, YANG Wei, JIN Hai, et al. Comparative study of turbulence model on prediction of wind pressure and wind velocity around rectangular buildings[J].BuildingScience, 2006, 22(5): 1-5, 28.

[16] 張群峰, 何鴻濤. 不同湍流模型數值模擬三維軸對稱凸體分離流動的比較[J].科學技術與工程, 2009, 9(13): 3693-3697, 3703.

ZHANG Qunfeng, HE Hongtao. Numerical study on three-dimensional separated flow on an axisymmetric bump by different turbulent models[J].ScienceTechnologyandEngineering, 2009, 9(13): 3693-3697, 3703.

[17] 胡煜, 伍青安, 袁越, 等. 風電引起3p閃變的仿真分析[J].電力自動化設備, 2013, 33(3): 108-111.

HU Yu, WU Qing'an, YUAN Yue, et al. Simulative analysis of 3p voltage flicker caused by wind farm integration[J].ElectricPowerAutomationEquipment, 2013, 33(3): 108-111.

[18] 胡斌, 劉雙. 風電場風速對風機葉片性能影響的模擬分析[J].動力工程學報, 2016, 36(12): 993-999.

HU Bin, LIU Shuang. Influence of wind speed on the performance of wind turbine blades[J].JournalofChineseSocietyofPowerEngineering, 2016, 36(12): 993-999.