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

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

0 引 言

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

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

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

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

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

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

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

1 計(jì)算模型與方法

1.1 計(jì)算模型

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

1.2 網(wǎng)格劃分

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

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

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

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

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

1.3 控制方程

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

(1)

(2)

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

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

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

2.1 風(fēng)機(jī)氣動(dòng)力特性

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

風(fēng)機(jī)葉片的葉尖線速度與對(duì)應(yīng)風(fēng)速之比稱為翼尖速比，所研究風(fēng)機(jī)模型的翼尖速比為

(3)

式中：ω為風(fēng)機(jī)轉(zhuǎn)速；R為風(fēng)機(jī)半徑；U為來(lái)流風(fēng)速。

推力是指風(fēng)沿來(lái)流方向?qū)︼L(fēng)機(jī)表面的作用力，推力系數(shù)CT為

(4)

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

風(fēng)能利用系數(shù)Cp用于表征風(fēng)輪對(duì)風(fēng)能的捕獲能力，其表達(dá)式為

(5)

式中：M為風(fēng)力扭矩。

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

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

表1 風(fēng)機(jī)翼尖速比為8時(shí)不同計(jì)算方法推力系數(shù)和風(fēng)能利用系數(shù)對(duì)比

2.2 尾渦分布

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

(6)

式中：Ω為旋渦強(qiáng)度；S為剪切應(yīng)變率。圖4中Q= 0.01，并用風(fēng)速方向速度分量值染色。

圖4 風(fēng)機(jī)尾渦結(jié)構(gòu)示例

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

圖5 不同截面流場(chǎng)速度軸向分布圖

2.3 瞬時(shí)流線圖

2.3.1 葉片表面正負(fù)壓分界線

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

圖6 葉片表面正負(fù)壓分界線示例

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

圖7 葉片表面流線圖

2.3.2 來(lái)流方向流線

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

圖8 來(lái)流方向上不同截面的瞬時(shí)流線

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

2.3.3 翼型截面流線

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

圖9 5 MW風(fēng)機(jī)葉片5個(gè)典型截面

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

圖10 不同翼型截面的瞬時(shí)流線

3 結(jié) 論

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

數(shù)值計(jì)算結(jié)果的對(duì)比說(shuō)明，該數(shù)值計(jì)算方法對(duì)于計(jì)算風(fēng)機(jī)氣動(dòng)力特性和流場(chǎng)真實(shí)情形的模擬較為準(zhǔn)確，計(jì)算結(jié)果吻合較好，可作為今后風(fēng)機(jī)研究的方法和手段。