胡丹梅， 潘 揚， 張建平

(上海電力學院 能源與機械工程學院， 上海 200090)

塔影效應是指由于塔架對流場的干涉，導致塔架附近葉片的工作特性發(fā)生變化。當水平軸風力機的葉片經過塔架附近時，空氣作用于葉片上的氣動載荷會隨風速變化而變化，造成單個葉片的受力狀況與其他葉片之間存在較大差異，這種差異對于風力機的正常運行是較為不利的。

塔影效應的研究起源于國外，其概念早在1968年便被提出，但當時的風力機技術處于起步階段，作為非主要影響因素的塔影效應并未受到重視。進入21世紀后，隨著風電產業(yè)的興起，風電從業(yè)者開始認識到塔架對風力機的影響不可忽略，因此對塔影效應的研究逐漸增多。2001年，Wang等[1]建立了基于規(guī)定尾流渦模型和近尾流動力學模型的高分辨塔影效應模型，其計算結果在一定程度上與實際測量結果相符。Das等[2]為計算塔架造成的風力機周期波動，建立了有關塔影效應和風剪切效應的時域模型。對于3葉片風力機來說，一定時間內產生的振蕩次數是旋轉頻率P的3倍，這種周期性波動稱為3P效應。我國早期的研究著重于塔影效應對電力系統(tǒng)方面的影響。張玉良等[3]研究了塔影效應對風功率的影響，導出了下風向風力機的風功率計算式。李少林等[4]在設計動態(tài)風力機模擬器時涉及了對塔影效應的模擬，并計算出受塔影效應影響的等效風速，用以替代空間平均風速。由于葉片經過塔架時等效風速會發(fā)生脈動，可用等效風速計算出塔影效應引起的轉矩脈動。張雪等[5]和董升等[6]均進行了基于等效風速法的風力機塔影效應和風剪切效應仿真研究，并得出了相似的結果。范忠瑤等[7]以東方汽輪機有限公司的2.5 MW上風向型DF90為模型，進行了塔影效應的數值模擬，研究表明塔影效應使單個葉片載荷波動20%～40%，風輪載荷波動6%～12%，且風輪的旋轉作用導致塔前來流始終發(fā)生偏轉，引起塔架兩側受力不均。

目前，國內外對于塔影效應的研究尚未成熟，特別是對兆瓦級大型風力機的數值模擬相對較少。筆者采用Fluent軟件對美國可再生能源實驗室(NREL)的5 MW風力機進行模擬，分析在不同來流風速下塔影效應的表現，為后續(xù)減弱塔影效應的研究提供參考。

1 風力機3D建模及網格劃分

1.1 建立模型

選取NREL 5 MW風力機作為研究對象，風力機為3葉片上風向，葉輪直徑為126 m，葉片長度為61.5 m[8]。將文獻[8]中葉片17個截面處的弦長和翼型數據導入ProE軟件中，通過連接翼型生成葉片模型。將輪轂簡化為球體，機艙簡化為長方體，將塔筒設為上底面直徑為3.87 m、下底面直徑為6 m的圓臺，將包括葉片在內的各組件進行組裝，得到整機模型。為了模擬出風力機運行時葉輪轉動的工作狀態(tài)，需將風力機的外部流場(計算域)分為2個部分，即將葉輪完全包圍并隨之旋轉的圓盤狀旋轉域和旋轉域外圍不隨葉輪旋轉的靜止域，如圖1所示。旋轉域半徑為70 m，深度為8 m，關于旋轉平面對稱，旋轉中心距地面高度為90 m；靜止域入口距葉輪旋轉平面為250 m，出口距葉輪旋轉平面距離為1 000 m，寬為600 m，高為400 m，靜止域為長方體。

圖1 計算域模型

1.2 網格劃分

采用ICEM-CFD網格劃分軟件對旋轉域和靜止域分別進行結構化網格劃分。劃分旋轉域網格時，先劃分出三分之一旋轉域網格，再沿旋轉軸進行周期性旋轉，得到整個旋轉域的計算網格，如圖2所示。

圖2 旋轉域劃分網格

對葉片外圍流場進行O剖分，并單獨劃分流動邊界層，如圖3所示。由于風力機模型的葉輪直徑達126 m，幾何尺寸很大，為保證網格數量在可計算范圍之內，沿葉片展向的網格節(jié)點數較少，導致邊界層網格的長和寬較大。為防止網格長寬比過大，影響計算精度[9]，采用壁面函數法處理邊界層；為滿足Y+要求[10]，第1層網格高度為1.5 mm，網格增長比為1.1，共14層。

圖3 不同葉高處的網格分布

Fig.3 Grid distribution at different blade heights

表1 不同網格數的計算結果比較

選取額定工況進行網格無關性驗證，結果如表1所示?？紤]到計算的準確性，選擇網格數為626.75萬，其中旋轉域的網格數為346.37萬，靜止域的網格數為280.38萬。

2 計算方法

2.1 控制方程

基于N-S方程進行數值模擬，考慮到湍流模型需與壁面函數法兼容，以及葉輪的旋轉效應，選用帶有湍流漩渦修正的RNGk-ε模型，離散格式為2階迎風，利用Simple算法進行求解[11]。計算采用滑移網格法，旋轉域設為以葉輪旋轉軸為軸進行旋轉，風力機部件視為剛性，不考慮流固耦合。

RNGk-ε模型的控制方程為[12]：

(1)

(2)

式中：ρ為空氣密度；k為湍流脈動動能；ε為脈動耗散率；ρε為耗散項的密度；t為時間；xi為空間某一方向(i、j=1，2，3)；ui為速度在某一方向上的分量(i、j=1，2，3)；μeff為有效黏性系數；Gk為產生的湍流動能；Gb為浮力產生的湍流動能；YM為可壓縮湍流中過渡擴散產生的波動；αk為湍流動能普朗特數的倒數；αε為湍流耗散率普朗特數的倒數；C1ε、C2ε和C3ε均為模型默認常數；Rε、Sk和Sε均為用戶自定義項。

2.2 塔影效應

風力機葉輪的轉矩為[13]：

(3)

(4)

式中：λ為葉尖速比；R為葉片旋轉半徑，m；ω為葉片旋轉角速度，rad/s；Cp為功率系數；β為槳距角，(°)；Tm為轉矩，N·m；v為運動黏度，m2/s。

由式(3)可知，葉片槳距角一定時，葉輪轉矩與葉輪上游處的風速有關。

3 模擬結果與分析

3.1 計算可信度驗證

3.1.1 功率曲線驗證

為驗證數值模擬的準確性，選取25 m/s、20 m/s、15 m/s、11.4 m/s(額定風速)、9 m/s、7 m/s和5 m/s作為來流風速進行模擬。表2為各工況下的葉片槳距角[8]。如圖4所示，模擬值和設計值的功率曲線在額定工況下符合情況較好，最大相對誤差為4.6%，可認為模型和計算方法具有可信度。

表2 不同工況下的葉片槳距角

圖4 風力機功率曲線

3.1.2 3P效應

設置11.4 m/s的均勻來流，旋轉域的旋轉角速度為12.1 rad/min，采用滑移網格法，設時間步長為0.027 548 2 s(對應旋轉域旋轉角度為2°)，進行非定常計算。計算過程中對整個風力機葉輪和3個葉片分別進行轉矩監(jiān)測，待計算時間足夠長，轉矩趨于穩(wěn)定后，記錄葉輪旋轉1周過程中的轉矩，如圖5所示。由圖5可以看出，在葉輪旋轉1周的過程中，相比于無塔筒，有塔筒時轉矩出現3次大幅下降，且方位角相互間隔120°，符合3P效應。

圖5 1個旋轉周期內整個葉輪的轉矩

3.2 模擬結果分析

3.2.1 塔影效應對葉片工作特性的影響

圖6和圖7分別為在11.4 m/s、9 m/s、7 m/s和5 m/s的來流風速下1根葉片在旋轉1周過程中轉矩和軸向推力的變化曲線，可較直觀地反映塔筒對葉片工作特性的影響。葉片尖部垂直向上時方位角為0°，葉片在塔筒正前方時方位角為180°。

由圖6(a)可知，在葉片的1個旋轉周期內，轉矩出現了1次谷值和2次峰值，谷值出現在方位角180°處(C點)，2個峰值分別出現在方位角140°(B點)和210°處(A點)，且210°處轉矩大于140°處轉矩。由圖6可知，轉矩的峰值和谷值出現的方位角不隨來流風速的變化而改變。

(a) 來流風速為11.4 m/s,葉輪轉速為12.1 rad/min

(b) 來流風速為9 m/s,葉輪轉速為10.2 rad/min

(c) 來流風速為7 m/s,葉輪轉速為8.3 rad/min

(d) 來流風速為5 m/s,葉輪轉速為7.3 rad/min

將最大轉矩和最小轉矩的差值與葉片旋轉1周的平均轉矩的比值定義為轉矩相對峰谷差，可表示為：

(5)

式中：TA為較大轉矩峰值；TC為最小轉矩；Tave為葉片旋轉1周的平均轉矩。

轉矩相對峰谷差的大小代表輸出轉矩的相對波動幅度，其值越大，轉矩的相對波動幅度越大。在4個來流風速下轉矩相對峰谷差分別為12.85%、12.88%、13.25%和14.93%。

將轉矩2個峰值之間的差值與葉片旋轉1周的平均轉矩的比值定義為轉矩相對峰峰差，可表示為：

(6)

式中：TB為較小轉矩峰值。

轉矩相對峰峰差反映了葉片轉矩波動規(guī)律的不對稱性，其值越大，說明葉片在塔筒兩側時工況的差異越大。在4個來流風速下轉矩相對峰峰差分別為1.36%、1.38%、1.83%和3.31%。綜上，來流風速越小，塔影效應引起的轉矩相對峰谷差越大，轉矩相對峰峰差越大。

對比圖6和圖7可知，葉輪所受軸向推力的變化同樣呈現峰-谷-峰規(guī)律。由圖7(a)可知，峰值分別出現在方位角140°(B點)和210°處(A點)，谷值出現在方位角180°處(C點)。如圖7所示，峰谷值出現的方位角與圖6中峰谷值的方位角一致，但140°處的峰值卻高于210°處的峰值。與轉矩相似，葉片所受軸向推力的峰值和谷值的方位角不隨來流風速的變化而變化。因此可認為塔影效應對葉片的氣動影響主要集中在葉片方位角140°～210°的時域內。

將最大軸向推力與最小軸向推力的差值與葉片旋轉1周的平均軸向推力的比值定義為推力相對峰谷差，其表達式為：

(7)

式中：NB為較大軸向推力峰值；NC為最小軸向推力；Nave為葉片旋轉1周的平均軸向推力。

推力相對峰谷差越大，說明葉片經過塔筒時受到的軸向推力變化越劇烈。在4個來流風速下軸向推力相對峰谷差分別為6.65%、6.16%、6.33%和5.68%。

將軸向推力2個峰值間的差值與葉片旋轉1周的平均軸向推力的比值定義為推力相對峰峰差，其表達式為：

(8)

式中：NA為較大軸向推力峰值。

推力相對峰峰差越大，說明葉片在塔筒兩側時所受推力情況差異越大。在4個來流風速下推力相對峰峰差分別為1.22%、1.02%、0.62%和0.20%。

與轉矩的變化特性相反，總體上葉片所受軸向推力的相對波動幅度隨來流風速的減小而減小。

(a) 來流風速為11.4 m/s

(b) 來流風速為9 m/s

(c) 來流風速為7 m/s

(d) 來流風速為5 m/s

葉片轉矩和軸向推力在方位角180°附近發(fā)生劇烈變化，說明葉片在塔筒附近時，氣動特性受塔筒影響會發(fā)生較大改變。圖8給出了塔影效應方程坐標，葉片的氣動特性在塔筒兩側表現出一定程度的不對稱性，這種不對稱性在根據等效風速法計算出的結果中并未體現[4-6]，而在同樣使用數值模擬方法的結果中卻有所體現[7,14]。

等效風速法[15]為：

v=vH+vtower

(9)

(10)

式中：v為等效風速；vH為輪轂高度處的風速；vtower為塔影效應下的風速；v0為空間平均風速；a為塔筒半徑；y為葉片微元到塔筒中線的水平距離；x為葉輪旋轉平面到塔筒中線的水平距離。

由式(10)可知，在等效風速法中，vtower在y方向上是關于塔筒對稱的函數，且未考慮葉輪旋轉造成的旋轉效應，因此推測這種不對稱性來源于旋轉效應對風速的影響。為了進一步分析這種不對稱性，圖9給出了額定風速下葉片在方位角140°、150°、180°和210°處、水平高度為77.4 m平面的風速云圖。

由圖9可知，與普通的圓柱繞流相比，塔筒附近的速度分布沿順時針方向有微小偏斜，偏斜方向與葉輪旋轉方向一致，說明旋轉效應造成了塔筒附近流場的不對稱性；葉片尾部有一個低速區(qū)，反映了風經過葉片后動能的損失。按風速損失程度由大到小排序分別為圖9(d)、圖9(a)、圖9(b)和圖9(c)，與轉矩出現的峰-谷-峰變化規(guī)律相符，也解釋了2個峰值在數值上的不對稱性。

(a) 方位角140°處

(b) 方位角150°處

(c) 方位角180°處

(d) 方位角210°處

為尋找能減弱塔影效應的有效對策，筆者在保持其他計算條件不變的情況下，將塔筒替換為桁架，并進行了模擬計算與分析。限于文章篇幅，僅列出分析結果，即與塔筒式風力機相比，桁架式風力機的葉輪和塔架所受氣動載荷的波動程度均明顯減小。

3.2.2 塔影效應下風力機流場及尾跡分析

計算出額定工況下風力機整機的周圍流場，并與無塔筒時的流場進行對比。圖10給出了葉片尖端指向地面時(對應葉片方位角180°)有、無塔筒情況下在葉片高度20%、40%、60%和80%處的流線及風速云圖。由圖10可知，與無塔筒相比，有塔筒時葉片周圍近距離處流線的變化不大，塔筒對葉片的影響主要體現在速度場方面。在有塔筒的情況下，葉片周圍的速度場受到明顯“擠壓”，葉片前緣處的高速區(qū)和尾緣處的低速區(qū)范圍均小于無塔筒的情況。塔筒下游出現了不同程度的渦流，渦流大小隨葉高的增大而減小，在20%和80%葉高處出現了渦流分離現象。對于風力機單機，塔影效應影響的關鍵區(qū)域在靠近葉根、弦長較大處，此處葉片近流場變化最為明顯。

塔筒除了會干涉葉片周圍的流場，也會對風力機下游的尾跡流場產生一定的影響。圖11給出了在額定風速下20%、40%、60%和80%葉高處風力機尾跡流場的風速云圖，發(fā)現塔筒對尾跡流場有明顯干擾。無塔筒時，風力機尾流速度場基本關于葉輪中軸線呈對稱分布。有塔筒時，除20%葉高處之外，其他3個截面處塔筒下游均出現了一個長度約為250 m(2D)的尾跡，且該尾跡在葉高較大處更為清晰。

(a) 無塔筒、20%葉高處

(b) 有塔筒、20%葉高處

(c) 無塔筒、40%葉高處

(d) 有塔筒、40%葉高處

(e) 無塔筒、60%葉高處

(f) 有塔筒、60%葉高處

(g) 無塔筒、80%葉高處

(h) 有塔筒、80%葉高處

4 結 論

(1) 當來流風為均勻來流且風向一定時，不同來流風速下葉片在1個旋轉周期內轉矩和軸向推力的峰值、谷值出現在相同的葉片方位角， 140°～210°方位角是塔影效應的主要影響時域。當來流風速減小時，在葉片掃掠過塔筒的過程中，其轉矩的相對波動幅度會減小，葉片所受軸向推力的相對波動幅度會增大。

(2) 在大型風力機的運行過程中，葉輪的旋轉效應對塔筒引起的風力機工況變化特性有較明顯的影響。其主要表現為葉片在靠近、遠離塔筒的過程中，葉片轉矩和軸向推力的變化規(guī)律具有不對稱性，因此在大型風力機塔影效應的研究中有必要考慮旋轉效應。

(a) 無塔筒、20%葉高處

(b) 有塔筒、20%葉高處

(c) 無塔筒、40%葉高處

(d) 有塔筒、40%葉高處

(e) 無塔筒、60%葉高處

(f) 有塔筒、60%葉高處

(g) 無塔筒、80%葉高處

(h) 有塔筒、80%葉高處

(3) 塔筒對風力機流場的速度分布有顯著影響。對于葉片附近流場，葉高越小處塔筒造成的影響越大，因此靠近葉根的部分是塔影效應影響風力機單機的關鍵區(qū)域；對于遠尾跡流場，葉高越大處塔筒造成的尾跡越明顯，塔筒尾跡長度約為葉輪直徑的2倍。

[1] WANG Tongguang, COTON F N. A high resolution tower shadow model for downwind wind turbines[J].JournalofWindEngineeringandIndustrialAerodynamics, 2001, 89(1): 873-892.

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

[3] 張玉良, 李仁年, 巫發(fā)明. 下風向風力機的塔影效應研究[J].山東建筑大學學報, 2008, 23(3): 243-246.

ZHANG Yuliang, LI Rennian, WU Faming. The research of tower shade effect of down wind turbine[J].JournalofShandongJianzhuUniversity, 2008, 23(3): 243-246.

[4] 李少林, 張興, 楊淑英, 等. 風力機動態(tài)模擬器的仿真研究[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.

[5] 張雪, 王維慶, 王海云, 等. 考慮風切和塔影效應的風力機風速模型[J].電測與儀表, 2015, 52(8): 56-60.

ZHANG Xue, WANG Weiqing, WANG Haiyun, et al. Model of wind speed considered wind shear and tower shadow effect on wind turbine[J].ElectricalMeasurement&Instrumentation, 2015, 52(8): 56-60.

[6] 董升, 周彪, 韓肖清. 基于剪切和塔影效應的風力機轉矩脈動模型仿真研究[J].山西電力, 2016(1): 1-5.

DONG Sheng, ZHOU Biao, HAN Xiaoqing. Simulation of wind turbine torque pulsation model based on shear and tower shadow effect[J].ShanxiElectricPower, 2016(1): 1-5.

[7] 范忠瑤, 康順, 趙萍. 上風向風力機塔影效應的數值模擬研究[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.

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

[9] 漆文邦, 鄭超. 基于單元不同長寬比網格劃分的有限元誤差分析研究[J].中國農村水利水電, 2010(2): 108-110, 113.

QI Wenbang, ZHENG Chao. A study of sensitivity to mesh distortion of isoparameter elements[J].ChinaRuralWaterandHydropower, 2010(2): 108-110, 113.

[10] 汪志成, 紀李文, 周書民, 等. 壁面參數對垂直軸風力發(fā)電機流體仿真的效果分析[J].現代制造工程, 2015(12): 43-46.

WANG Zhicheng, JI Liwen, ZHOU Shumin, et al. The analysis of numerical simulation performance for the vertical axis wind turbine on wall-parameter[J].ModernManufacturingEngineering, 2015(12): 43-46.

[11] 李少華, 匡青峰, 吳殿文, 等. 1.2 MW風力機整機流場的數值模擬[J].動力工程學報, 2011, 31(7): 551-556.

LI Shaohua, KUANG Qingfeng, WU Dianwen, et al. Numerical simulation on flow field of a 1.2 MW wind turbine[J].JournalofChineseSocietyofPowerEngineering, 2011, 31(7): 551-556.

[12] 任年鑫, 歐進萍. 大型海上風力機尾跡區(qū)域風場分析[J].計算力學學報, 2012, 29(3): 327-332.

REN Nianxin, OU Jinping. Numerical analysis for the wake zone of large offshore wind turbine[J].ChineseJournalofComputationalMechanics, 2012, 29(3): 327-332.

[13] 張弘鯤, 孟祥星. 塔影效應引起的風電機組輸出功率波動問題[J].東北電力技術, 2011, 32(4): 33-35.

ZHANG Hongkun, MENG Xiangxing. On output power fluctuation of wind generators caused by tower shadow effect[J].NortheastElectricPowerTechnology, 2011, 32(4): 33-35.

[14] 田仁斌, 李學敏, 何成明, 等. 上風型水平軸風力機塔架干涉的數值模擬[J].工程熱物理學報, 2014, 35(12): 2417-2420.

TIAN Renbin, LI Xuemin, HE Chengming, et al. Numerical simulation of upwind-type horizontal axis wind turbine blade-tower interaction[J].JournalofEngineeringThermophysics, 2014, 35(12): 2417-2420.

[15] DOLAN D S L, LEHN P W. Simulation model of wind turbine 3p torque oscillations due to wind shear and tower shadow[J].IEEETransactionsonEnergyConversion, 2006, 21(3): 717-724.