王新茹， 陳 剛， 肖龍飛， 溫斌榮， 田新亮

(上海交通大學 a.海洋工程國家重點實驗室； b.高新船舶與深海開發(fā)裝備協(xié)同創(chuàng)新中心， 上海 200240)

0 引 言

近年來，隨著海上風電的快速發(fā)展，海上風能的利用已成為可再生能源發(fā)展的重要方向[1]。風機的氣動力性能直接影響其發(fā)電效率和使用壽命。因此，對風機周圍流場的精細化模擬和分析作為海上風機組設計和研究的基礎，逐漸成為人們關注的關鍵技術之一[2-3]。

隨著計算機性能的提高和數(shù)值計算方法的發(fā)展，人們對風機周圍流場的數(shù)值分析方法進行了大量的研究和探索。

INGRAM[4]在一維動量理論的基礎上，提出葉素動量理論，該理論是計算風機葉片氣動力性能的經(jīng)典方法。但是，葉素動量理論在不同雷諾數(shù)和攻角下對二維翼型升力系數(shù)和阻力系數(shù)的依賴使其難以精確模擬復雜工況下風機的三維旋轉(zhuǎn)運動。為減小誤差，人們引入葉尖損失修正模型、動態(tài)失速修正模型、輪轂損失修正模型和偏航尾跡修正模型等諸多修正方法[5]。這些方法雖基本實現(xiàn)了快速高效地計算風機載荷情況，但在風機周圍流場分布方面卻顯得無能為力。為進一步精確計算風機氣動力性能，并提供風機周圍流場信息，KEITH等[6]和LANDAHI等[7]分別提出各種氣動力學勢流模型，如面源法、升力面理論等，對風機的三維旋轉(zhuǎn)進行較為準確的近似模擬。這類勢流方法在大型風機非定常氣動分析方面取得了較好的結(jié)果[8-10]。然而，這類方法通過工程模型的方式近似模擬葉片周圍的流體黏性和流動分離等關鍵現(xiàn)象，不能提供葉片周邊精密細致的流場分布。

近幾年，CFD方法應運而生，基于雷諾平均的Navier-Stokes方程，選取適當?shù)耐牧髂Ｐ蛯︼L機葉片周圍的三維流場進行數(shù)值模擬，通過計算機求解代數(shù)方程得到結(jié)果，最終還可將計算結(jié)果通過圖形或圖像進行可視化顯示，以便更好地理解和分析問題。但是，由于風機的翼型結(jié)構及其繞流場形狀十分復雜，風機的實型建模和高質(zhì)量的網(wǎng)格劃分變得困難，因此前期的CFD計算大多采用SORENSEN等[11]提出的制動線和制動面模型，將風機的實型簡化成原本并不存在的線或圓盤并離散成多個單元，求解所得的單位葉片所受升力和阻力便通過這些虛擬載體代入葉素理論從而獲得整個風機所受的升力和阻力。這種近似模型計算方法雖在計算精度上有所改善，但忽略了風機的結(jié)構對其氣動力性能和周圍流場影響，使結(jié)果與現(xiàn)實情況相比仍存在一定的差距。

ZHOU等[12]對風機模型進行原型建模，采用滑移網(wǎng)格方法考慮風機旋轉(zhuǎn)，探究風機結(jié)構細節(jié)及葉片旋轉(zhuǎn)對其周圍流場產(chǎn)生的影響。TRAN等[13]、LIU 等[14]和CAI 等[15]對完全考慮結(jié)構特點的NREL 5 MW風機進行實尺度建模，通過重疊網(wǎng)格的方法對旋轉(zhuǎn)風機進行數(shù)值模擬，探究風機的流場特征及其在風作用下的運動響應。但是，模型所采用網(wǎng)格數(shù)量較少，對流場細節(jié)的捕捉還不夠精細，使其在揭示物理現(xiàn)象和物理機制時略粗糙。

本文在對NREL 5 MW參考模型進行實尺度精細化建模的基礎上，采用極大規(guī)模的網(wǎng)格處理模型特征，還原葉片扭轉(zhuǎn)角、翼型過渡和葉尖等部分結(jié)構細節(jié)，采用OpenFOAM計算軟件，結(jié)合自主開發(fā)的旋轉(zhuǎn)動網(wǎng)格處理方法，利用pimpleDyMFoam求解器對風機葉片的氣動力性能和周圍流場進行全方位、多角度的精細化模擬和分析。

圖1 計算模型示例(非正確比例)

1 計算模型與方法

1.1 計算模型

所采用的風機模型為NREL 5 MW參考風機[16]，這是一款上風向三葉片水平軸風機。出于計算簡化考慮，未考慮風機塔柱的影響。風機盤面直徑D=126 m，其所處圓柱型流域如圖1所示。風機旋轉(zhuǎn)平面到入口的距離約2D，距出口為10D。來流風速取10 m/s，坐標軸x的正方向與風向一致，風機轉(zhuǎn)速為順時針(+x)12 r/min。大氣參數(shù)分別取：密度ρ=1.29 kg/m3，運動學黏性系數(shù)ν=14.8×10-6m2/s。

1.2 網(wǎng)格劃分

CFD的本質(zhì)是將連續(xù)介質(zhì)進行區(qū)域化離散，通過插值運算將物理場中的連續(xù)量近似到各離散點上，再求解各離散點上的物理量從而獲得所求流場中各物理量的連續(xù)性分布。也就是說，計算是在生成的各離散網(wǎng)格點上進行的。因此，將一個給定形狀的物體在所需研究的流場內(nèi)進行合適的網(wǎng)格劃分十分重要。

本文首先通過ANSYS中的網(wǎng)格劃分軟件ICEM將圓柱型流域內(nèi)部嵌套一個O型網(wǎng)格用于模擬加密的葉片網(wǎng)格，對風機葉片尾流部分的網(wǎng)格進行適當加密，以期捕捉到更真實、精細的尾渦結(jié)構。而后，利用snappyHexMesh工具對網(wǎng)格進行捕捉、分割、移動，以及添加邊界層網(wǎng)格。由于風機為實體，因此在流場計算中需將其內(nèi)部的網(wǎng)格掏空，并增加相應的邊界層。整個計算域的網(wǎng)格劃分結(jié)果如圖2所示，風機葉片及邊界層劃分如圖3所示，網(wǎng)格總數(shù)量為1 520萬個。

圖2 全計算域網(wǎng)格劃分結(jié)果

圖3 風機葉片附近網(wǎng)格及邊界層

同時，基于OpenFOAM 1.6版本，在求解器中引入課題組研發(fā)的整體旋轉(zhuǎn)網(wǎng)格的動網(wǎng)格方法[17]，避免滑移網(wǎng)格法信息交互誤差，使整個計算網(wǎng)格隨風機葉片一起旋轉(zhuǎn)，保證長時間高精度高效率的計算。

1.3 控制方程

風機外流場的馬赫數(shù)一般小于0.3，故可將流體視為牛頓不可壓縮流體。非定常流場的控制方程(N-S方程)[18]為

(1)

(2)

式(1)和式(2)中：p為壓力；ρ為流場密度；xi、xj為笛卡爾坐標系下點的坐標(i、j=1、2、3，分別表示x、y、z方向)；ui為速度在xi方向上的分量；t為時間。

計算流體動力學開源軟件OpenFOAM采用有限體積法對N-S方程進行離散化，基于張量法和面向?qū)ο蠹夹g[19]，用于解決連續(xù)介質(zhì)中的流體力學問題。采用PIMPLE(PISO和SIMPLE混合算法)的離散格式求解速度與壓力的耦合，用于插值、梯度、拉普拉斯和散度計算的空間格式分別為線性、高斯線性、高斯線性修正和高斯線性。對流項采用二階迎風格式，時間項采用二階Crank-Nicholson格式。以0.000 7 s步長在上海交通大學超級計算中心平臺中以64核計算兩周得到計算結(jié)果。

2 計算結(jié)果分析

2.1 風機氣動力特性

風機的基本工作原理是利用空氣流經(jīng)風機葉片時葉片所產(chǎn)生的升力或阻力推動葉片轉(zhuǎn)動，通過傳動系統(tǒng)將轉(zhuǎn)矩傳到發(fā)電機，將風能轉(zhuǎn)化為電能。因此，風機作為一種流體機械存在效率問題。通常，主要通過推力系數(shù)和風能利用系數(shù)評價風機的氣動力性能[20]。

風機葉片的葉尖線速度與對應風速之比稱為翼尖速比，所研究風機模型的翼尖速比為

(3)

式中：ω為風機轉(zhuǎn)速；R為風機半徑；U為來流風速。

推力是指風沿來流方向?qū)︼L機表面的作用力，推力系數(shù)CT為

(4)

式中：T為風向方向的風載荷分量；ρ0為空氣密度。

風能利用系數(shù)Cp用于表征風輪對風能的捕獲能力，其表達式為

(5)

式中：M為風力扭矩。

通過對OpenFOAM中翼尖速比為8時計算得到的風機所受推力和風能利用系數(shù)與其他計算方法在相同翼尖速比時所得結(jié)果[8， 21-22]進行對比，分析風機氣動力性能細節(jié)以及所采用模擬方法的實用性，如表1所示。

觀察可知，風機所受推力系數(shù)CT和風能利用系數(shù)CP的數(shù)值計算結(jié)果與其他計算方法所得結(jié)果基本吻合，間接證明所采用數(shù)值模擬方法的實用性。

表1 風機翼尖速比為8時不同計算方法推力系數(shù)和風能利用系數(shù)對比

2.2 尾渦分布

當水平軸風機的葉片被限制在一定長度時，葉片旋轉(zhuǎn)會導致在風機下游部分產(chǎn)生螺旋形渦。產(chǎn)生的渦主要存在于兩個區(qū)域：風機的葉尖處，也稱葉尖渦；風機下游旋轉(zhuǎn)軸的水平方向處。根據(jù)距離風機葉片的遠近，一般將風機尾渦區(qū)劃分為近尾渦區(qū)(x/D≤ 1)和遠尾渦區(qū)(x/D>1)，近尾渦區(qū)的各氣動力學參數(shù)會對風機自身的運行成本、發(fā)電功率和使用壽命產(chǎn)生較大影響，而遠尾渦區(qū)則會影響風電場中其他風機的氣動力性能。通過渦的可視化體現(xiàn)，可清楚地了解風機的尾渦結(jié)構及其下游流場形態(tài)。采用Q準則對尾渦區(qū)的渦結(jié)構進行表示，Q的定義為

(6)

式中：Ω為旋渦強度；S為剪切應變率。圖4中Q= 0.01，并用風速方向速度分量值染色。

圖4 風機尾渦結(jié)構示例

圖5給出了流場中順風向距離風機葉片x/D=-1.0、0.5、1.0、1.5、2.0時，風速在z軸方向上的對比曲線。z/R用于表征z方向距離與風機半徑的比值，同時引入α=u/U因子表征風機葉片不同截面高度處流場的速度損失情況?？梢钥闯?，流場的風速分布關于x軸基本對稱。在遠尾渦區(qū)，流場速度受風機輪轂和葉片阻隔的影響而有所下降，且距離越近，速度下降越明顯；在近尾渦區(qū)，由于風機輪轂和葉片(尤其是葉尖處)的旋轉(zhuǎn)帶動作用，流場風速在z= 0 m和z>40 m處已基本恢復。

圖5 不同截面流場速度軸向分布圖

2.3 瞬時流線圖

2.3.1 葉片表面正負壓分界線

風流過旋轉(zhuǎn)的風機后，其周圍流場在風機的前后表面會形成正負壓區(qū)域，圖6所示的即為風機表面流場正負壓分界線情況。圖中黑色表示正壓，灰色表示負壓。由觀察可知，風機葉片的壓力面成正負壓交替分布狀，而在其吸力面，流場壓力分布以負壓為主。

圖6 葉片表面正負壓分界線示例

在葉片表面用三維流線表征流場的流動特性，清晰可見在分界線周圍流線分支、匯聚，且在風機輪轂后方形成渦，如圖7所示。

圖7 葉片表面流線圖

2.3.2 來流方向流線

為更形象細致地研究風機對于風機流場的干擾作用，沿來流方向分別取x/D= -1.0、0.5、1.0、1.5、2.0處的流場截面，背景以壓力染色，分別觀察不同截面處流線與壓力特征，圓形表征風機盤面大小，y、z軸為各截面處真實坐標值，如圖8所示。

圖8 來流方向上不同截面的瞬時流線

由觀察可得：在葉片截面處，葉片周圍流場的壓力增大，瞬時流線在葉片吸力面處形成旋渦；距離風機x/D= 0.5處的流場最為明顯的特征是截面流線在旋轉(zhuǎn)軸中心處形成旋渦，流場壓力降為負值，且越靠近旋轉(zhuǎn)軸中心，三葉片投影區(qū)域的反向壓力最大，說明該處截面受旋轉(zhuǎn)葉片的影響顯著；在x/D= 1處，流場壓力仍為負值，但流線與之前不同，呈類海螺狀，旋渦形狀與之前相比變大；近尾渦區(qū)(見圖8(d))之后，壓力負值逐漸減小，流線螺旋態(tài)減小，風機葉片對流場的影響減弱。

2.3.3 翼型截面流線

選取5個典型截面的壓強系數(shù)，截面分布(r/R=0.30、0.47、0.63、0.80、0.95)如圖9所示，分析翼型截面的流場流線分布。

圖9 5 MW風機葉片5個典型截面

如圖10所示，流場壓力在翼型處存在明顯分界，壓力面壓力增大，吸力面壓力減少，過渡區(qū)與翼型處于同一截面，壓力成階梯狀由大至小分布，流線在壓力面和引力面分別形成渦，隨著翼型截面的減小，產(chǎn)生的渦也相應減小直至消失。

圖10 不同翼型截面的瞬時流線

3 結(jié) 論

基于OpenFOAM，采用整體旋轉(zhuǎn)動網(wǎng)格處理方法對大型風機周圍三維流場進行了模擬與分析。從數(shù)值模擬的流場來看，由于翼型的存在，風機周圍流場結(jié)構復雜，葉尖渦與尾流區(qū)分界明顯，來流方向上不同截面的流場特性隨著其與風機表面距離遠近的變化而變化，渦的產(chǎn)生對風機自身部件和風電場其他風機都會產(chǎn)生一定的作用。

數(shù)值計算結(jié)果的對比說明，該數(shù)值計算方法對于計算風機氣動力特性和流場真實情形的模擬較為準確，計算結(jié)果吻合較好，可作為今后風機研究的方法和手段。