李 麗， 潘 揚(yáng)， 胡丹梅

(1.國網(wǎng)湖南省電力有限公司防災(zāi)減災(zāi)中心(電網(wǎng)輸變電設(shè)備防災(zāi)減災(zāi)國家重點(diǎn)實(shí)驗(yàn)室)，長沙 410129；2.上海電力學(xué)院 能源與機(jī)械工程學(xué)院，上海 200090)

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

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

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

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

筆者以NREL 5 MW風(fēng)力機(jī)[9]為模型，通過數(shù)據(jù)交互模塊System Coupling將流體模擬軟件Fluent與固體分析軟件Transient Structural進(jìn)行聯(lián)動，實(shí)現(xiàn)風(fēng)力機(jī)固體結(jié)構(gòu)與其周圍流場的雙向耦合，展開了關(guān)于塔影效應(yīng)的流固耦合研究。

1 計(jì)算模型與網(wǎng)格劃分

1.1 物理模型與網(wǎng)格

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

1.1.2 劃分計(jì)算網(wǎng)格

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

圖2 葉輪殼網(wǎng)格

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

1.1.3 可信度驗(yàn)證

為驗(yàn)證計(jì)算的準(zhǔn)確性，分別利用627萬、190萬結(jié)構(gòu)化網(wǎng)格和415萬、259萬非結(jié)構(gòu)化網(wǎng)格在非流固耦合情況下進(jìn)行計(jì)算，所得功率偏差分別為0.28%、3.6%、3.2%和5.3%。由于流固耦合計(jì)算非常占用計(jì)算機(jī)資源，為節(jié)省計(jì)算時(shí)間和提高計(jì)算機(jī)工作穩(wěn)定性，最終采用259萬非結(jié)構(gòu)化四面體網(wǎng)格，計(jì)算誤差約為5%，滿足計(jì)算要求。

1.2 計(jì)算方法與耦合方程

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

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

(1)

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

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

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

圖3 Workbench模塊結(jié)構(gòu)圖

2 計(jì)算結(jié)果分析

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

2.1 風(fēng)力機(jī)葉輪氣動載荷

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

圖4 葉輪轉(zhuǎn)矩曲線

圖5 葉輪軸向推力曲線

2.2 葉輪形變與應(yīng)力

圖6給出了模擬時(shí)間t為10 s、20 s和30 s時(shí)風(fēng)力機(jī)葉輪的形變放大圖，放大倍數(shù)約為19倍。由于葉輪旋轉(zhuǎn)周期T為5 s，在時(shí)間為5 s的整數(shù)倍時(shí)，z方向的葉片位于塔筒正前方，如圖6(b)中葉片1所示。由圖6可以看出，在10 s時(shí)風(fēng)力機(jī)葉輪工況尚未穩(wěn)定，葉片的形變較大，位于塔筒前方葉片的軸向揮舞形變和周向擺振形變均小于其他2根葉片；30 s時(shí)葉輪工況基本穩(wěn)定，3根葉片的擺振形變程度較為接近，但揮舞形變?nèi)杂幸欢ú町?，此現(xiàn)象與文獻(xiàn)[18]相符。

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

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

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

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

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

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

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

表2 葉片尖部位移

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

圖7 葉輪最大形變曲線

圖8 葉輪最大等效應(yīng)力曲線

2.3 流固耦合與非流固耦合結(jié)果對比

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

(a) 非耦合，無塔筒

(b) 耦合，無塔筒

(c) 非耦合，有塔筒

(d) 耦合，有塔筒

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

3 結(jié) 論

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

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

(3)在風(fēng)力機(jī)運(yùn)行過程中塔筒會不斷對葉輪產(chǎn)生氣動干擾，導(dǎo)致整個葉片的形變和應(yīng)力要劣于無塔筒時(shí)，塔影效應(yīng)對風(fēng)力機(jī)運(yùn)行的穩(wěn)定性有不利影響。

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

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

(a) 非耦合，壓力面

(b) 耦合，壓力面

(c) 非耦合，吸力面

(d) 耦合，吸力面

：

[1] 范忠瑤, 康順, 趙萍. 上風(fēng)向風(fēng)力機(jī)塔影效應(yīng)的數(shù)值模擬研究[J].工程熱物理學(xué)報(bào), 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] 李少林, 張興, 楊淑英, 等. 風(fēng)力機(jī)動態(tài)模擬器的仿真研究[J].太陽能學(xué)報(bào), 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級風(fēng)力發(fā)電機(jī)風(fēng)輪葉片流固耦合場強(qiáng)度分析[D]. 鄭州: 鄭州大學(xué), 2011.

[11] 何福添. 水平軸風(fēng)力機(jī)葉片的設(shè)計(jì)及流固耦合計(jì)算[D]. 武漢: 武漢理工大學(xué), 2013.

[12] 李媛. 風(fēng)力機(jī)葉片流固耦合數(shù)值模擬[D]. 北京: 華北電力大學(xué), 2013.

[13] 胡丹梅, 張志超, 孫凱, 等. 風(fēng)力機(jī)葉片流固耦合計(jì)算分析[J].中國電機(jī)工程學(xué)報(bào), 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] 方平治, 顧明, 談建國, 等. 數(shù)值模擬大氣邊界層中解決壁面函數(shù)問題方法研究[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] 金新陽, 楊偉, 金海, 等. 數(shù)值風(fēng)工程應(yīng)用中湍流模型的比較研究[J].建筑科學(xué), 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] 張群峰, 何鴻濤. 不同湍流模型數(shù)值模擬三維軸對稱凸體分離流動的比較[J].科學(xué)技術(shù)與工程, 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] 胡煜, 伍青安, 袁越, 等. 風(fēng)電引起3p閃變的仿真分析[J].電力自動化設(shè)備, 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] 胡斌, 劉雙. 風(fēng)電場風(fēng)速對風(fēng)機(jī)葉片性能影響的模擬分析[J].動力工程學(xué)報(bào), 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.