基于致動(dòng)線方法的風(fēng)力機(jī)氣動(dòng)數(shù)值模擬

朱 翀，王同光，鐘 偉

（南京航空航天大學(xué) 江蘇省風(fēng)力機(jī)設(shè)計(jì)高技術(shù)研究重點(diǎn)實(shí)驗(yàn)室，江蘇 南京 210016）

0 引 言

為了獲得風(fēng)力機(jī)葉片的氣動(dòng)性能，目前使用的最主要的工具是結(jié)合二維翼型性能數(shù)據(jù)的葉素動(dòng)量（BEM）理論。由于其簡(jiǎn)便性，工業(yè)應(yīng)用中使用BEM理論設(shè)計(jì)新的風(fēng)力機(jī)葉片。然而，為了獲得風(fēng)力機(jī)周圍流場(chǎng)的精確信息，為了獲得風(fēng)力機(jī)或風(fēng)力機(jī)組更精確的載荷信息和發(fā)電功率等，需要更精致的技術(shù)，比如計(jì)算流體力學(xué)（CFD）方法［1－3］。目前 CFD 方法已經(jīng)成為研究風(fēng)力機(jī)空氣動(dòng)力學(xué)的常用工具。然而在使用CFD方法對(duì)風(fēng)力機(jī)尾流流場(chǎng)進(jìn)行研究時(shí)，為了同時(shí)滿足風(fēng)力機(jī)葉片附近區(qū)域和尾流區(qū)域的數(shù)值精度需要大量網(wǎng)格和很高的計(jì)算成本，這限制了CFD方法在風(fēng)力機(jī)尾流流場(chǎng)研究中的應(yīng)用。結(jié)合BEM理論的簡(jiǎn)便性與CFD方法對(duì)流場(chǎng)的精確模擬性，Shen等學(xué)者提出了致動(dòng)線（AL）方法［4］。Troldborg在其博士論文中使用致動(dòng)線方法對(duì)Tj?reborg葉片與NM80葉片進(jìn)行了數(shù)值模擬［5］；Shen使用致動(dòng)線方法對(duì)MEXICO實(shí)驗(yàn)葉片進(jìn)行了數(shù)值模擬［6］，均取得了較理想的結(jié)果。而在國(guó)內(nèi)，致動(dòng)線方法尚未有較系統(tǒng)的研究。

致動(dòng)線方法是融合了引入沿葉片展向分布的體積力代表風(fēng)力機(jī)葉片的致動(dòng)線技術(shù)與三維Navier－Stokes方程的一種方法。風(fēng)力機(jī)葉片繞軸旋轉(zhuǎn)，產(chǎn)生氣動(dòng)力，從流體微團(tuán)的角度來(lái)看，流體微團(tuán)感受到的是動(dòng)量的變化，故通過(guò)在Navier－Stokes方程中引入動(dòng)量源項(xiàng)（即體積力）來(lái)表示葉片的這種作用。因此致動(dòng)線方法的優(yōu)勢(shì)在于使用體積力來(lái)代替風(fēng)力機(jī)葉片，其所需要的網(wǎng)格單元數(shù)遠(yuǎn)遠(yuǎn)小于為了捕捉葉片幾何細(xì)節(jié)所耗費(fèi)的網(wǎng)格單元數(shù)，而且在計(jì)算過(guò)程中也避免了花費(fèi)大量資源去計(jì)算葉片附面層，可以把更多的網(wǎng)格與計(jì)算資源投入到風(fēng)力機(jī)的尾流區(qū)域，最終的計(jì)算結(jié)果中包含的流場(chǎng)信息（比如速度、壓力、渦量等各種參數(shù)）也不比CFD方法獲得的少，所以致動(dòng)線方法是研究風(fēng)力機(jī)尾流流場(chǎng)的很好的工具，國(guó)外的學(xué)者在使用致動(dòng)線方法研究風(fēng)力機(jī)尾流流場(chǎng)方面已經(jīng)取得極大的進(jìn)展［4－6］。致動(dòng)線方法與CFD方法都是通過(guò)Navier－Stokes方程組來(lái)求解流場(chǎng)信息的，所以在數(shù)值精度方面應(yīng)該是完全一樣的。在計(jì)算結(jié)果的準(zhǔn)確性方面，雖然致動(dòng)線方法使用體積力代替風(fēng)力機(jī)葉片，無(wú)法獲得葉片附面層流場(chǎng)的信息，但其使用了風(fēng)洞實(shí)驗(yàn)測(cè)得的翼型氣動(dòng)力數(shù)據(jù)，一定程度上彌補(bǔ)了以上簡(jiǎn)化造成的計(jì)算結(jié)果準(zhǔn)確性的降低。而B(niǎo)EM理論雖然能夠快速得到葉片的載荷、功率等信息，但是完全無(wú)法獲得相關(guān)的流場(chǎng)信息。因此致動(dòng)線方法在風(fēng)力機(jī)研究方面有其獨(dú)特的優(yōu)勢(shì)。但是由于代替葉片的體積力是由二維翼型性能數(shù)據(jù)計(jì)算而得，故致動(dòng)線方法非常依賴于所提供的二維翼型性能數(shù)據(jù)的準(zhǔn)確性。

本文使用致動(dòng)線方法對(duì)南京航空航天大學(xué)自主研發(fā)的NH1500葉片進(jìn)行數(shù)值模擬，并將致動(dòng)線方法的數(shù)值模擬結(jié)果與BEM理論、CFD方法的數(shù)值模擬結(jié)果以及風(fēng)洞實(shí)驗(yàn)［7］結(jié)果進(jìn)行了系統(tǒng)的比較，來(lái)突顯致動(dòng)線方法區(qū)別于BEM理論與CFD方法獨(dú)特的優(yōu)勢(shì)。

1 致動(dòng)線方法

1.1 控制方程

本文數(shù)值模擬使用基于有限體積法的解算器。將風(fēng)力機(jī)周圍流場(chǎng)看作不可壓有粘流動(dòng)，故求解的控制方程組為旋轉(zhuǎn)參考系下的不可壓雷諾平均Navier－Stokes方程組，其表達(dá)式為：

式中：Ω為以S為邊界的控制體；W 表示守恒變量；Fc表示對(duì)流通量；Fv表示粘性通量；Q表示體積力源項(xiàng)。W，F(xiàn)c，F(xiàn)v，Q 的表達(dá)式如下：

式中：n為單位法向量；ρ為空氣密度；u，v和w 為速度的三個(gè)分量；p為壓強(qiáng)為剪切應(yīng)力；fe為參考系旋轉(zhuǎn)導(dǎo)致的體積力加速度；fε為代表葉片作用的體積力，將在后面進(jìn)行詳細(xì)的推導(dǎo)；

逆變速度V的表達(dá)式為：

方程組中各個(gè)變量都是定義在旋轉(zhuǎn)參考系下的，該參考系的轉(zhuǎn)軸與風(fēng)力機(jī)風(fēng)輪轉(zhuǎn)軸重合，轉(zhuǎn)速與風(fēng)輪轉(zhuǎn)速相同，是一個(gè)非慣性參考系。為了修正因參考系旋轉(zhuǎn)而附加給氣流的離心力和科式力，在方程組中定義體積力加速度為：

式中：ω為參考系的旋轉(zhuǎn)角速度；r為位置矢量。

本文采用對(duì)邊界層湍流和自由剪切湍流都有較好模擬效果的SST k－ω模型［8］來(lái)計(jì)算湍流粘性。

1.2 葉片氣動(dòng)模型

代表風(fēng)力機(jī)葉片的體積力是基于葉素假設(shè)，根據(jù)當(dāng)?shù)氐牧鲌?chǎng)信息與已知的二維翼型性能計(jì)算所得。圖1為葉片的一個(gè)翼型截面圖，其中x方向?yàn)閬?lái)流方向。

圖1 葉片的翼型截面Fig.1 Cross－sectional airfoil element of blade

葉片翼型截面的當(dāng)?shù)厮俣葹椋?/p>

其中，Ω為角速度，r為葉片翼型截面徑向位置，Vx和Vθ分別為葉片翼型截面的軸向速度和切向速度。

軸向速度

切向速度

本文計(jì)算使用坐標(biāo)系為笛卡爾坐標(biāo)系，葉片位于平面x＝0內(nèi)。u、v、w分別為流場(chǎng)中葉片翼型截面x、y、z方向的速度，Δy和Δz為葉片翼型截面與葉根在y和z方向上的距離。

葉片翼型截面的當(dāng)?shù)厮俣扰c風(fēng)力機(jī)葉片平面之間的夾角為：

葉片翼型截面的當(dāng)?shù)赜菫椋?/p>

其中γ為當(dāng)?shù)貥嘟恰?/p>

在確定了葉片翼型截面當(dāng)?shù)厮俣扰c當(dāng)?shù)赜侵?，即可以得到單位展長(zhǎng)的升力與阻力：

其中CL＝CL（α，Re）和CD＝CD（α，Re）分別為升力系數(shù)和阻力系數(shù)，Re為基于當(dāng)?shù)叵议L(zhǎng)c的雷諾數(shù)，eL和eD為升力和阻力方向的單位向量。

1.3 高斯分布

為了使代表葉片的體積力在流場(chǎng)中光順的分布，將體積力在流場(chǎng)中以致動(dòng)線為中心線呈三維高斯分布［4］。

經(jīng)過(guò)高斯分布的單位體積力為：

分布因子：

其中d＝｜l－sei｜為網(wǎng)格單元中心點(diǎn)與致動(dòng)線上第i個(gè)點(diǎn)之間的距離，ε為控制氣動(dòng)力分布的分布因子。

最后，將分布后的單位體積力fε（l）作為體積力源項(xiàng)加載到CFD計(jì)算的動(dòng)量方程中，即為方程（1）中的體積力源項(xiàng)Q中的fε。而其他的物理量（如速度、壓力等）則通過(guò)求解方程（1）而得。

2 計(jì)算模型與計(jì)算條件

計(jì)算模型選用南京航空航天大學(xué)自主研發(fā)的NH1500風(fēng)力機(jī)，該風(fēng)力機(jī)為三葉片式，葉片長(zhǎng)度為40.5m。選擇4個(gè)具有代表性的運(yùn)行狀態(tài)進(jìn)行計(jì)算，具體計(jì)算狀態(tài)見(jiàn)表1，其中Vin為入流風(fēng)速，Ω為風(fēng)輪轉(zhuǎn)速，λ為葉尖速比。

表1 計(jì)算狀態(tài)Table 1 Computational states

致動(dòng)線方法的計(jì)算域結(jié)構(gòu)如圖2所示，速度入口距離致動(dòng)線為5倍風(fēng)輪半徑，壓力出口距離致動(dòng)線為20倍風(fēng)輪半徑，遠(yuǎn)場(chǎng)邊界距離中心線為9倍風(fēng)輪半徑。針對(duì)NH1500風(fēng)力機(jī)三葉片的軸對(duì)稱特點(diǎn)，應(yīng)用了120°旋轉(zhuǎn)周期邊界條件，使得計(jì)算域只有實(shí)際區(qū)域的三分之一，從而使網(wǎng)格單元數(shù)和計(jì)算量減少。周期邊界上一側(cè)的網(wǎng)格節(jié)點(diǎn)和另一側(cè)的網(wǎng)格節(jié)點(diǎn)一一對(duì)應(yīng)，緊鄰周期邊界一側(cè)的計(jì)算域外的“鏡像單元”信息由緊鄰另一側(cè)的計(jì)算域內(nèi)的單元提供。為了充分利用網(wǎng)格，將計(jì)算域分為內(nèi)場(chǎng)和外場(chǎng)，內(nèi)場(chǎng)為致動(dòng)線附近的區(qū)域，對(duì)該區(qū)域的網(wǎng)格進(jìn)行加密處理，而外場(chǎng)區(qū)域內(nèi)的網(wǎng)格則較為稀疏。實(shí)際計(jì)算域網(wǎng)格拓?fù)浣Y(jié)構(gòu)如圖3所示。整個(gè)計(jì)算域采用六面體網(wǎng)格，網(wǎng)格單元數(shù)為1.5×106。

圖2 計(jì)算域結(jié)構(gòu)圖Fig.2 Structure diagram of computational domain

圖3 計(jì)算域網(wǎng)格拓?fù)浣Y(jié)構(gòu)Fig.3 Topology of computational domain

為了驗(yàn)證致動(dòng)線方法的可行性，除了使用致動(dòng)線方法，還使用了BEM理論與CFD方法對(duì)NH1500葉片進(jìn)行了數(shù)值模擬作為參照。其中BEM理論［9－10］主要是以成熟商業(yè)軟件BLADED作為參照而編寫(xiě)的，考慮了經(jīng)典的普朗特葉尖損失修正［11］，計(jì)算結(jié)果與BLADED軟件計(jì)算結(jié)果相一致。CFD方法則是對(duì)NH1500葉片的附面層進(jìn)行了精確的數(shù)值模擬，計(jì)算區(qū)域及網(wǎng)格拓?fù)浣Y(jié)構(gòu)與致動(dòng)線方法網(wǎng)格基本一致，網(wǎng)格單元數(shù)為4×106。

用于研究NH1500葉片氣動(dòng)性能的風(fēng)洞實(shí)驗(yàn)是在中國(guó)氣動(dòng)研究與發(fā)展中心12m×16m的低速風(fēng)洞中完成的。

3 分布因子ε的選取

分布因子ε是致動(dòng)線方法中非常重要的一個(gè)參數(shù)，它關(guān)系著體積力的分布范圍，是能否合理表現(xiàn)出葉片在流場(chǎng)中作用的關(guān)鍵。在此對(duì)不同分布因子ε的數(shù)值模擬結(jié)果進(jìn)行比較分析。一般情況下，分布因子ε的取值都是網(wǎng)格中致動(dòng)線部分的網(wǎng)格單元尺寸的倍數(shù)，即ε＝n×Δs，其中n為倍數(shù)，Δs為致動(dòng)線網(wǎng)格單元的尺寸。如果分布因子ε選取過(guò)小，則體積力在流場(chǎng)中作用范圍過(guò)于集中，與實(shí)際葉片對(duì)流場(chǎng)的影響不相符，導(dǎo)致計(jì)算結(jié)果與真實(shí)情況相差較大，同時(shí)還可能會(huì)影響計(jì)算的收斂性。如果分布因子ε選取過(guò)大，則會(huì)使體積力在流場(chǎng)中作用范圍過(guò)大，影響葉尖流場(chǎng)的分布，這當(dāng)然也是與實(shí)際葉片對(duì)流場(chǎng)的影響不相符的。

以入流風(fēng)速V0＝9m／s的情況為例進(jìn)行分析，以n作為變量，分別取1，1.5，2，2.5，3，3.5，4。圖4為不同分布因子ε下沿葉片展向分布的切向力因子Ct，圖中還增加了CFD方法數(shù)值模擬的結(jié)果作為參照。由圖4可以看出，當(dāng)n取1的時(shí)候，計(jì)算所得的切向力因子Ct明顯過(guò)小。這說(shuō)明n取1時(shí)，因體積力的分布范圍過(guò)于集中而不能真實(shí)合理的反映葉片在流場(chǎng)中的作用，從而該計(jì)算結(jié)果錯(cuò)誤。而當(dāng)n取1～4時(shí)，計(jì)算所得的切向力因子Ct在葉片大部分區(qū)域基本相同，只在葉尖處有區(qū)別；n取4時(shí)，葉尖處的切向力因子Ct有明顯的突變。這說(shuō)明n取4時(shí)，體積力的分布范圍過(guò)大，影響了葉尖的流場(chǎng)分布。綜合比較得出當(dāng)n取2時(shí)，體積力的分布較為合理，能夠真實(shí)反應(yīng)葉片對(duì)周圍流場(chǎng)區(qū)域的影響。

圖4 沿葉片展向分布的切向力因子CtFig.4 Tangential force coefficient Ct distributed along the blade

因此分布因子的選擇將會(huì)影響葉片氣動(dòng)力的計(jì)算以及代表葉片作用的體積力在流場(chǎng)中的作用范圍，選擇合適的分布因子才能夠獲得較準(zhǔn)確的葉片流場(chǎng)。

4 數(shù)值模擬結(jié)果與分析

4.1 載荷分布

圖5為入流風(fēng)速8m／s下致動(dòng)線方法、BEM理論、CFD方法計(jì)算出的沿葉片展向的法向力系數(shù)Cn和切向力系數(shù)Ct分布，其中橫坐標(biāo)為葉片當(dāng)?shù)卣瓜蛭恢门c葉片長(zhǎng)度的比值，縱坐標(biāo)分別為法向力系數(shù)Cn和切向力系數(shù)Ct。由于三種方法中對(duì)葉根處的處理方法并不相同，導(dǎo)致葉根處結(jié)果相差較大，不具有參考價(jià)值；且葉根處的發(fā)電量在整個(gè)葉片發(fā)電量中所占比例極小，故在此不考慮葉根處的情況。

圖5 沿葉片展向分布的法向力因子Cn與切向力因子CtFig.5 Normal force coefficient Cn and tangential force coefficient Ct distributed along the blade

總的來(lái)說(shuō)，圖5中三種方法計(jì)算所得的沿葉片展向的法向力系數(shù)Cn和切向力系數(shù)Ct分布趨勢(shì)大致相同。但是致動(dòng)線方法與BEM理論的計(jì)算結(jié)果更為接近，這是因?yàn)橹聞?dòng)線方法與BEM理論都是通過(guò)二維翼型性能插值計(jì)算得到各個(gè)截面的氣動(dòng)力，進(jìn)而計(jì)算出其它相關(guān)信息；而CFD方法則是完全通過(guò)對(duì)流場(chǎng)的數(shù)值模擬得到所有的信息。致動(dòng)線方法把得到的各個(gè)截面的氣動(dòng)力通過(guò)體積力的方式作用到流場(chǎng)中，再通過(guò)求解Navier－Stokes方程得到葉片性能及周圍流場(chǎng)的信息；BEM理論則是將葉片離散成一個(gè)個(gè)截面，通過(guò)一維動(dòng)量守恒定律進(jìn)行分析，相鄰截面互不干擾。致動(dòng)線方法中葉片不同截面之間是有相互影響的，而B(niǎo)EM理論中不同截面之間是互不干擾的，這也是致動(dòng)線方法與BEM理論的一個(gè)區(qū)別。

除了所使用的計(jì)算方法對(duì)結(jié)果造成的影響，致動(dòng)線方法與BEM理論計(jì)算所需的二維翼型性能數(shù)據(jù)對(duì)計(jì)算結(jié)果有更大的影響。入流風(fēng)速8m／s下，致動(dòng)線方法計(jì)算所得的沿葉片展向的法向力系數(shù)Cn和切向力系數(shù)Ct要大于BEM理論的，而B(niǎo)EM理論的大于CFD方法的。這是因?yàn)橹聞?dòng)線方法與BEM理論計(jì)算所用的二維翼型風(fēng)洞實(shí)驗(yàn)性能數(shù)據(jù)是雷諾數(shù)3×106下的，遠(yuǎn)遠(yuǎn)大于入流風(fēng)速8m／s時(shí)葉片各截面真實(shí)的雷諾數(shù)，過(guò)大的雷諾數(shù)導(dǎo)致計(jì)算出的氣動(dòng)力偏高，從而使致動(dòng)線方法與BEM理論的計(jì)算結(jié)果較CFD方法要偏大。而且葉片展向各個(gè)翼剖面的雷諾數(shù)也并不相同，使用同一雷諾數(shù)的翼型氣動(dòng)力數(shù)據(jù)來(lái)計(jì)算，也會(huì)使結(jié)果產(chǎn)生一定的誤差。

而在葉尖處CFD方法計(jì)算所得的法向力系數(shù)Cn和切向力系數(shù)Ct卻大于致動(dòng)線方法與BEM理論計(jì)算所得的，這是因?yàn)槿~片葉尖由于葉尖繞流效應(yīng)的存在而使周圍流場(chǎng)三維效應(yīng)顯著。雖然BEM理論已經(jīng)考慮了經(jīng)典的普朗特葉尖損失修正，但是畢竟普朗特葉尖損失修正是經(jīng)驗(yàn)公式，不可能完全體現(xiàn)出葉尖效應(yīng)；而CFD方法是對(duì)葉尖直接進(jìn)行模擬，故應(yīng)該可以較好的體現(xiàn)出葉尖效應(yīng)，但是由于葉尖處并不是常規(guī)的翼型形狀，在葉尖采集處理所得的法向力系數(shù)Cn和切向力系數(shù)Ct會(huì)有一些誤差，故CFD方法中的Cn、Ct只采集到葉片的展向94%位置處。而致動(dòng)線方法在沒(méi)有使用任何經(jīng)驗(yàn)公式修正的情況下，計(jì)算所得葉尖處的Cn、Ct與BEM理論使用了普朗特葉尖損失修正的計(jì)算所得的Cn、Ct已經(jīng)比較接近了，與CFD方法計(jì)算所得有一些小差距。這是因?yàn)殡m然說(shuō)致動(dòng)線方法是通過(guò)直接求解Navier－Stokes方程來(lái)獲得葉尖附近的流場(chǎng)信息，能夠通過(guò)代表葉片作用的體積力來(lái)模擬葉尖繞流效應(yīng)，不同于無(wú)粘假設(shè)的BEM理論，不需要使用傳統(tǒng)的普朗特葉尖損失修正公式來(lái)對(duì)葉尖進(jìn)行修正。但是由于葉尖繞流造成的三維效應(yīng)，靠近葉尖處的翼型截面性能不再是二維的，而致動(dòng)線方法所使用的翼型性能數(shù)據(jù)卻是二維的。

4.2 功率系數(shù)

風(fēng)力機(jī)功率系數(shù)Cp是風(fēng)力機(jī)葉片的主要?dú)鈩?dòng)性能參數(shù)，體現(xiàn)了風(fēng)力機(jī)葉片的發(fā)電效率，是衡量風(fēng)力機(jī)葉片氣動(dòng)性能優(yōu)劣的重要指標(biāo)。圖6為NH1500葉片在不同葉尖速比下風(fēng)洞實(shí)驗(yàn)、CFD方法、致動(dòng)線方法所得的功率系數(shù)。從圖6中可以看出無(wú)論是風(fēng)洞實(shí)驗(yàn)還是CFD方法或致動(dòng)線方法，都是在葉尖速比為9.5左右取得最大功率系數(shù)Cp。風(fēng)洞實(shí)驗(yàn)的最大功率系數(shù)Cp為0.492，CFD方法的為0.505，致動(dòng)線方法的為0.521。不僅僅是風(fēng)洞實(shí)驗(yàn)獲得的最大功率系數(shù)Cp要略小于CFD方法與致動(dòng)線方法計(jì)算所得的最大功率系數(shù)Cp，而且所有實(shí)驗(yàn)獲得的功率系數(shù)Cp都略小于CFD方法與致動(dòng)線方法相對(duì)應(yīng)葉尖速比下的功率系數(shù)Cp。由于風(fēng)洞實(shí)驗(yàn)中所使用的為1／16的模型，其雷諾數(shù)必然是小于CFD方法與致動(dòng)線方法中的雷諾數(shù)，造成了風(fēng)洞實(shí)驗(yàn)中葉片的升力偏低而阻力偏高，這也就解釋了CFD方法與致動(dòng)線方法所得的功率系數(shù)Cp要略大于風(fēng)洞實(shí)驗(yàn)的功率系數(shù)Cp。而致動(dòng)線方法所得的功率系數(shù)Cp則要比CFD方法的略高一些，是因?yàn)橹聞?dòng)線方法中所使用的二維翼型性能數(shù)據(jù)為雷諾數(shù)3×106下的，與CFD方法中的雷諾數(shù)相比偏高；且在高葉尖速比，也即低風(fēng)速下，兩種方法中雷諾數(shù)的差距更大，所得的功率系數(shù)Cp相對(duì)而言相差也更大一些。

圖6 功率系數(shù)CpFig.6 The power coefficient Cp

5 尾流區(qū)域的簡(jiǎn)要分析

葉片壓力面的壓力高，吸力面的壓力低，壓力差推動(dòng)氣流繞過(guò)葉尖，形成葉尖渦。風(fēng)力機(jī)的葉尖渦在來(lái)流風(fēng)速的疊加下，不容易從跡線或速度矢量上直接顯示出漩渦的位置，因而引入渦量來(lái)顯示渦。圖7是使用致動(dòng)線方法計(jì)算的入流風(fēng)速9m／s下拖出的葉尖渦的渦量等值面，左右兩個(gè)平面上顯示的是軸向速度等值線。從圖7中可以清晰的看出從葉尖拖出的葉尖渦在尾流區(qū)域呈螺旋結(jié)構(gòu)向下游發(fā)展，其誘導(dǎo)速度以一定的周期性規(guī)律作用于尾流流場(chǎng)。

圖7 尾流的渦量等值面和軸向速度等值線圖Fig.7 Isosurface of the vorticity and contours of the axial velocity in the wind turbine wake flow filed

為了更好的觀測(cè)尾流區(qū)域流場(chǎng)的變化，取葉片所在的軸向平面的軸向速度等值線圖，如圖8所示。圖8中，通過(guò)葉尖渦對(duì)尾流流場(chǎng)的誘導(dǎo)作用可以清楚的看出葉尖渦在尾流區(qū)域是呈膨脹趨勢(shì)向下游發(fā)展。這一系列的結(jié)論都是與CFD方法對(duì)尾流區(qū)域葉尖渦發(fā)展趨勢(shì)的研究結(jié)論［12－15］相一致的。這也驗(yàn)證了使用致動(dòng)線方法研究尾流流場(chǎng)的可行性。

圖8 軸向速度等值線圖Fig.8 Contours of the axial velocity

為了準(zhǔn)確的分析尾流區(qū)域流場(chǎng)的計(jì)算精度，選取軸向誘導(dǎo)因子這一典型流場(chǎng)參數(shù)進(jìn)行研究。圖9為葉片所在的軸向平面的軸向誘導(dǎo)因子等值線圖，橫軸為風(fēng)力機(jī)風(fēng)輪下游位置x與風(fēng)輪半徑R的比值，縱軸為風(fēng)輪徑向位置y與風(fēng)輪半徑R的比值。圖中橫坐標(biāo)為0處黑線代表葉片的位置。從圖9中可以看出葉片和葉片周圍的軸向誘導(dǎo)因子接近1／3，這與經(jīng)典動(dòng)量理論中氣動(dòng)性能最優(yōu)的風(fēng)力機(jī)是一致的。尾流區(qū)域的外部和中間的誘導(dǎo)因子都要小于尾流區(qū)域的主體部分。很明顯的是，軸向誘導(dǎo)因子從葉片位置向下游方向不斷增加，即是軸向速度不斷減小。在葉片下游遠(yuǎn)方，軸向誘導(dǎo)因子最終達(dá)到了0.5，是葉片和葉片周圍誘導(dǎo)因子的兩倍，這也驗(yàn)證了經(jīng)典動(dòng)量葉素理論的假設(shè)。

圖9 軸向誘導(dǎo)因子等值線圖Fig.9 Contours of the axial induced factor

6 結(jié) 論

本文以南京航空航天大學(xué)自主研發(fā)的NH1500葉片為計(jì)算模型，使用致動(dòng)線方法、BEM理論與CFD方法分別進(jìn)行了數(shù)值模擬。并且與風(fēng)洞實(shí)驗(yàn)進(jìn)行了比較。通過(guò)系統(tǒng)的比較，在對(duì)流場(chǎng)細(xì)節(jié)的捕捉方面，致動(dòng)線方法的計(jì)算精度已經(jīng)達(dá)到了CFD方法的計(jì)算精度，與風(fēng)洞實(shí)驗(yàn)吻合。為了捕捉風(fēng)力機(jī)葉片的附面層，CFD方法耗費(fèi)了大量網(wǎng)格與計(jì)算資源，而致動(dòng)線方法則不需要去求解風(fēng)力機(jī)葉片附面層，可以把更多的網(wǎng)格與計(jì)算資源投入到風(fēng)力機(jī)尾流流場(chǎng)。因此致動(dòng)線方法相比于CFD方法，節(jié)約了大量的網(wǎng)格與計(jì)算資源，卻能夠達(dá)到CFD方法對(duì)風(fēng)力機(jī)尾流流場(chǎng)模擬同樣的精度，最終的計(jì)算結(jié)果中包含的流場(chǎng)信息也不比CFD方法獲得的少。而B(niǎo)EM理論則完全不能夠獲得風(fēng)力機(jī)葉片周圍的流場(chǎng)信息。因此致動(dòng)線方法是研究風(fēng)力機(jī)尾流流場(chǎng)的有力工具，特別是對(duì)于風(fēng)場(chǎng)中多臺(tái)風(fēng)力機(jī)之間的相互干擾的研究，致動(dòng)線方法的優(yōu)勢(shì)相比于CFD方法更加明顯。但是致動(dòng)線方法中使用的二維翼型性能數(shù)據(jù)的準(zhǔn)確性在一定程度上也制約著致動(dòng)線方法的發(fā)展。研究三維效應(yīng)對(duì)二維翼型性能的影響，使翼型性能數(shù)據(jù)滿足致動(dòng)線方法計(jì)算的需求，將有利于提高致動(dòng)線方法的準(zhǔn)確性。

［1］SANDERSE B，Van Der PIJL S P，KOREN B.Review of computational fluid dynamics for wind turbine wake aerodynamics［EB／OL］. ［2011－02－04］.http：／／onlinelibrary.wiley.com／doi／10.1002／We.458／abstract

［2］VERMEER L J，S?RENSEN J N，CRESPO A.Wind turbine wake aerodynamics［J］.Progress in Aerospace Sciences，2003，39（6／7）：467－510.

［3］王琦峰，陳伯君，安亦然.大型風(fēng)力機(jī)三維空氣動(dòng)力學(xué)數(shù)值模擬［J］.空氣動(dòng)力學(xué)學(xué)報(bào)，2011，29（6）：810－814.

［4］S?RENSEN J N，SHEN W Z.Numerical modeling of wind turbine wakes［J］.Journal of Fluids Engineering，2002，124：393－399.

［5］TROLDBORG N.Actuator line modeling of wind turbine wakes［D］.Denmark：Technical University of Denmark，2008.

［6］SHEN W Z，ZHU W J，S?RENSEN J N.Actuator line／Navier－Stokes computations for the MEXICO rotor：comparison with detailed measurements［EB／OL］.［2011－10－03］.http：／／onlinelibrary.wiley.com／doi／10.1002／we.510／abstract

［7］肖京平，陳立，許波峰，武杰.1.5MW風(fēng)力機(jī)氣動(dòng)性能研究［J.空氣動(dòng)力學(xué)學(xué)報(bào)，2011294529－534.

［8］MENTER F R.Two－equation eddy－viscosity turbulence models for engineering applications［J］.AIAA Journal，1994，32（8）：1598－1605.

［9］GLAUERT H.Aerodynamic theory，vol.4［M］.Berlin，Germany：Julius Springer：1935：169－360.

［10］MARTIN O L，Hansen.Aerodynamics of wind turbines（2nd Edition）［M］.UK：Earthscan，2008，45－62.

［11］PRANDTL L，BETZ A.Vier Abhandlungen zur hydrodynamik und aerodynamik［C］.G?ttingen，city，1927：88－92.

［12］鐘偉，王同光.風(fēng)力機(jī)葉尖渦的數(shù)值模擬［J］.南京航空航天大學(xué)學(xué)報(bào)，2011435640－644.

［13］PORTE－AGEL F，WU Y T，LU H，et al.Large－eddy simulation of atmospheric boundary layer flow through wind turbines and wind farms［J］.Journal of Wind Engineering and Industrial Aerodynamics，2011，99（4）：154－168.

［14］JOHANSEN J，S?RENSEN N N，MICHELSEN J A，et al.Detached－eddy simulation of flow around the NREL phase VI blade［J］.Wind Energy，2002，5：185－197.

［15］楊瑞.水平軸風(fēng)力機(jī)尾流結(jié)構(gòu)及空氣動(dòng)力學(xué)性能的研究［D］.中國(guó)：蘭州理工大學(xué)，2009.