張立軍 趙昕輝 馬東辰 米玉霞 王旱祥 姜 浩

中國(guó)石油大學(xué)(華東)機(jī)電工程學(xué)院，青島，266580

0 引言

垂直軸風(fēng)力機(jī)可以接收來(lái)自任何方向的風(fēng)，其增速齒輪箱和發(fā)電機(jī)可以安裝在地面，運(yùn)行維修方便，但風(fēng)力機(jī)運(yùn)行時(shí)，其葉片常發(fā)生動(dòng)態(tài)失速(低葉尖速比時(shí)尤為明顯)[1]，因此準(zhǔn)確計(jì)算動(dòng)態(tài)失速下的葉片氣動(dòng)力系數(shù)是分析與設(shè)計(jì)垂直軸風(fēng)力機(jī)的關(guān)鍵。動(dòng)態(tài)失速是指葉片攻角發(fā)生周期性或非定常變化時(shí)，翼型的失速攻角比靜態(tài)失速攻角要大得多，且翼型氣動(dòng)特性曲線(通常為法向力系數(shù)和切向力系數(shù)隨攻角的變化曲線)要明顯滯后于靜態(tài)曲線的現(xiàn)象。翼型發(fā)生動(dòng)態(tài)失速時(shí)測(cè)得的動(dòng)態(tài)氣動(dòng)力系數(shù)與翼型靜止時(shí)的靜態(tài)氣動(dòng)力系數(shù)相差較大。

目前研究翼型動(dòng)態(tài)失速的方法主要有三種：①以Navier-Stokes方程為基礎(chǔ)的CFD數(shù)值方法[2-4]；②基于面元法和邊界層理論的黏性與無(wú)黏耦合算法[5]；③基于實(shí)驗(yàn)數(shù)據(jù)建立的半經(jīng)驗(yàn)動(dòng)態(tài)失速模型方法[6-7]。其中，半經(jīng)驗(yàn)動(dòng)態(tài)失速模型計(jì)算效率高且通用性好，通過(guò)適當(dāng)?shù)男拚憧捎糜诖怪陛S風(fēng)力機(jī)翼型的動(dòng)態(tài)氣動(dòng)系數(shù)計(jì)算?，F(xiàn)有的半經(jīng)驗(yàn)動(dòng)態(tài)失速模型主要有以下幾類(lèi)[7-11]：基于對(duì)動(dòng)態(tài)失速延遲數(shù)值修正的B-V模型和MIT模型；用三階微分方程描述氣動(dòng)系數(shù)的ONERA模型；基于翼型繞流流動(dòng)特性的B-L動(dòng)態(tài)失速模型、丹麥Ris國(guó)家實(shí)驗(yàn)室建立的用于評(píng)價(jià)水平軸風(fēng)力機(jī)氣動(dòng)特性的Ris模型。國(guó)內(nèi)外一些科研單位對(duì)上述模型進(jìn)行適當(dāng)修正，并用于各自研究領(lǐng)域[12-14]。這些模型均是針對(duì)航空翼型或水平軸風(fēng)力機(jī)翼型的，在計(jì)算垂直軸風(fēng)力機(jī)葉片動(dòng)態(tài)失速氣動(dòng)性能時(shí)，需要結(jié)合其實(shí)際工作特性對(duì)這些模型進(jìn)行修正，但相關(guān)研究較少。

由于B-L模型和MIT更多地考慮了運(yùn)動(dòng)翼型的繞流物理特性，能更好地模擬翼型動(dòng)態(tài)失速特性且實(shí)用性較強(qiáng)，因此本文首先介紹了目前常用的B-L模型和MIT模型的計(jì)算方法，并基于垂直軸風(fēng)力機(jī)工作時(shí)的非定常氣動(dòng)特性對(duì)這兩種模型進(jìn)行修正。然后結(jié)合雙致動(dòng)盤(pán)多流管理論，用修正后模型分別求解Sandia實(shí)驗(yàn)室17 m垂直軸風(fēng)力機(jī)葉片的動(dòng)態(tài)切向力系數(shù)和動(dòng)態(tài)法向力系數(shù)，將計(jì)算結(jié)果與實(shí)驗(yàn)數(shù)據(jù)進(jìn)行對(duì)比，分析這兩種修正模型對(duì)不同風(fēng)區(qū)的預(yù)測(cè)精度。最終得到一種在整個(gè)風(fēng)區(qū)都能保持較高計(jì)算精度的垂直軸風(fēng)力機(jī)葉片動(dòng)態(tài)氣動(dòng)力系數(shù)計(jì)算模型。

1 B-L模型修正

B-L模型[9]屬于半經(jīng)驗(yàn)?zāi)Ｐ?，它雖依賴(lài)經(jīng)驗(yàn)常數(shù)，但能部分預(yù)測(cè)動(dòng)態(tài)失速現(xiàn)象。B-L模型將坐標(biāo)系OXY固定在翼型上，如圖1所示。翼型的升力系數(shù)CL和阻力系數(shù)CD向弦線方向(X方向)和弦線法向(Y方向)投影，可得切向力系數(shù)CT、法向力系數(shù)CN。圖1中，W是作用在翼型上的合成風(fēng)速，可通過(guò)雙致動(dòng)盤(pán)多流管理論求得；α為翼型攻角，是合成風(fēng)速方向與翼型弦線方向的夾角。

圖1 翼型的氣動(dòng)力系數(shù)Fig.1 Aerodynamic coefficients of airfoils

B-L模型主要是確定非定常附著流、分離流動(dòng)和動(dòng)態(tài)失速渦對(duì)翼型氣動(dòng)力系數(shù)的影響，在此基礎(chǔ)上得出翼型在動(dòng)態(tài)失速下的動(dòng)態(tài)氣動(dòng)力系數(shù)。

1.1 非定常附著流

(1)

αE=φC(t)Δα+αn-1

φC(t)=A0-A1exp(-b1s)-A2exp(-b2s)

Δα=αn-αn-1

式中，CNα為對(duì)應(yīng)翼型的靜態(tài)法向力系數(shù)曲線線性段的斜率；αE為B-L模型中定義的有效攻角；φC(t)為指數(shù)形式的延遲函數(shù)；t為翼型的厚度；s為量綱一時(shí)間，s=2vt/c；v為平均來(lái)流風(fēng)速；c為翼型弦長(zhǎng)；αn-1為翼型在n-1時(shí)刻的攻角；常數(shù)A0=1.0，A1=0.30，A2=0.70，b1=0.35，b2=0.68[7]。

原始的B-L模型采用階躍函數(shù)來(lái)描述翼型攻角隨時(shí)間的變化規(guī)律。對(duì)于垂直軸風(fēng)力機(jī)而言，其翼型攻角α的變化規(guī)律為類(lèi)正弦曲線：

(2)

式中，θ為翼型當(dāng)前所處的方位角；λ為風(fēng)力機(jī)局部葉尖速比，λ=V/(ωR)；V為流體誘導(dǎo)速度；ω為風(fēng)輪旋轉(zhuǎn)角速度；R為風(fēng)輪旋轉(zhuǎn)半徑。

本文使用式(2)代替原B-L模型的階躍函數(shù)，既符合垂直軸風(fēng)力機(jī)工作的實(shí)際情況，也使攻角變化的連續(xù)性更強(qiáng)。

脈沖分量反映了翼型運(yùn)動(dòng)引起的表面法向速度的氣動(dòng)響應(yīng)，可以由活塞理論進(jìn)行求解：

(3)

式中，Ma為馬赫數(shù)。

對(duì)垂直軸風(fēng)力機(jī)而言，Ma?1，式(3)的計(jì)算結(jié)果接近于0，故可以忽略附著流的脈沖效應(yīng)[12-13]，進(jìn)而可對(duì)現(xiàn)有的B-L模型進(jìn)行簡(jiǎn)化。因此非定常附著流的法向力系數(shù)為

(4)

切向力系數(shù)為[7]

(5)

式中，CLα為對(duì)應(yīng)翼型的靜態(tài)升力系數(shù)的曲線線性段的斜率。

1.2 分離流動(dòng)

對(duì)分離流動(dòng)的描述是B-L模型中至關(guān)重要的一部分，而描述分離流動(dòng)最關(guān)鍵的是確定流動(dòng)分離點(diǎn)的位置。通過(guò)Kirchhoff流動(dòng)定理可以得到靜態(tài)流動(dòng)分離點(diǎn)與翼型在定常情況下的法向力系數(shù)、切向力系數(shù)：

(6)

(7)

式中，fN、fT分別為流動(dòng)分離點(diǎn)在翼型表面法向和切向的位置，fN=x/c，fT=x/c；x為流動(dòng)分離點(diǎn)距后緣的長(zhǎng)度；α0為零升攻角。

將式(6)、式(7)重新整理可得對(duì)應(yīng)的法向力系數(shù)和切向力系數(shù)的分離點(diǎn)位置：

(8)

(9)

在小攻角變化情況下，fN和fT的值非常接近，故在原B-L模型中，令流動(dòng)分離點(diǎn)f=fN，其中，fN為式(8)的計(jì)算結(jié)果，并將f作為計(jì)算CT和CN共同的流動(dòng)分離點(diǎn)參數(shù)，再通過(guò)函數(shù)擬合方法得到分離點(diǎn)f與攻角α的連續(xù)函數(shù)關(guān)系。這樣對(duì)于任意的攻角α就都可以得到一個(gè)流動(dòng)分離點(diǎn)f，進(jìn)而通過(guò)式(6)、式(7)得到對(duì)應(yīng)攻角下的靜態(tài)法向力系數(shù)和切向力系數(shù)。對(duì)于攻角變化范圍較大的垂直軸風(fēng)力機(jī)翼型而言，通過(guò)上述方法得到的法向力系數(shù)和切向力系數(shù)與原始數(shù)據(jù)不能很好地吻合，文獻(xiàn)[12]對(duì)這一現(xiàn)象做了相應(yīng)的解釋。為此，本文在將B-L模型用于垂直軸風(fēng)力機(jī)翼型動(dòng)態(tài)失速性能研究時(shí)，獨(dú)立計(jì)算靜態(tài)法向力分離點(diǎn)fN和靜態(tài)切向力分離點(diǎn)fT，然后分別代入式(6)、式(7)計(jì)算法向力系數(shù)CN和切向力系數(shù)CT。

(10)

(11)

(12)

(13)

式中，Tf是B-L模型中定義的半經(jīng)驗(yàn)時(shí)間常數(shù)。

(14)

(15)

1.3 動(dòng)態(tài)失速渦

(16)

(17)

求解該式便可得到渦流法向力系數(shù)。

文獻(xiàn)[13]指出，翼型在大攻角情況下，渦流對(duì)切向力的影響較大，需要加入B-L模型加以考慮。渦流對(duì)切向力系數(shù)的影響可表示為

(18)

式中，τV為量綱一渦時(shí)間。

(19)

(20)

2 MIT模型修正

對(duì)于垂直軸風(fēng)力機(jī)，翼型攻角、氣動(dòng)力和風(fēng)速關(guān)系如圖2所示，V是誘導(dǎo)速度，W是誘導(dǎo)速度V和切向速度ωR的合成風(fēng)速，F(xiàn)L、FD分別為葉素受到的升力和阻力，切向力FT、法向力FN分別為FL和FD在翼型弦線方向和法向上的分量。

圖2 翼型氣動(dòng)參數(shù)分析示意圖Fig.2 Sketch map of airfoil aerodynamic parameters analysis

MIT模型將攻角α的變化過(guò)程劃分為4個(gè)階段，分別用不同的公式計(jì)算當(dāng)前攻角下翼型的動(dòng)態(tài)切向力系數(shù)和動(dòng)態(tài)法向力系數(shù)，具體計(jì)算方法如下：

(1)翼型攻角α從0開(kāi)始增大且低于靜態(tài)失速攻角αss時(shí)，動(dòng)態(tài)升力系數(shù)和阻力系數(shù)均取靜態(tài)值，這些值可通過(guò)查詢(xún)翼型升/阻力系數(shù)數(shù)據(jù)庫(kù)得到。數(shù)據(jù)庫(kù)未給出的值，可通過(guò)線性插值的方法求得。

(2)當(dāng)攻角繼續(xù)增大，超過(guò)靜態(tài)失速攻角αss但還未到達(dá)動(dòng)態(tài)失速攻角αds時(shí)，動(dòng)態(tài)阻力系數(shù)仍采用對(duì)應(yīng)攻角的靜態(tài)值，而升力系數(shù)為

(21)

式中，CLss為靜態(tài)失速攻角αss對(duì)應(yīng)的靜態(tài)升力系數(shù)。

MIT模型用誘導(dǎo)速度V作為翼型的特征風(fēng)速。對(duì)于垂直軸風(fēng)力機(jī)而言，作用在葉片上的風(fēng)速為誘導(dǎo)速度V和切向速度ωR的合成風(fēng)速W(其大小為W)。動(dòng)態(tài)失速攻角αds定義為

(22)

(23)

CDmax=CLmaxsinα

直到攻角開(kāi)始減小。這種假設(shè)極大地簡(jiǎn)化了計(jì)算，但對(duì)模型的計(jì)算精度有一些影響。

式中，CLs、CDs分別為翼型在當(dāng)前攻角下的升力系數(shù)和阻力系數(shù)的靜態(tài)值，均以指數(shù)方式衰減至其靜態(tài)值。

當(dāng)攻角再次增大時(shí)，重復(fù)上述過(guò)程，便可得到各個(gè)攻角下的動(dòng)態(tài)升阻力系數(shù)。通過(guò)下式可得到動(dòng)態(tài)切向力系數(shù)和動(dòng)態(tài)法向力系數(shù)：

(26)

3 修正模型中風(fēng)速參數(shù)計(jì)算

B-L修正模型和MIT修正模型都需要求解各個(gè)方位下作用在葉片上的合成風(fēng)速。本文中，合成風(fēng)速采用雙致動(dòng)盤(pán)多流管理論計(jì)算得出。雙致動(dòng)盤(pán)多流管模型用2個(gè)串聯(lián)的致動(dòng)盤(pán)分別表示上風(fēng)區(qū)的半個(gè)轉(zhuǎn)子掃掠面和下風(fēng)區(qū)的半個(gè)轉(zhuǎn)子掃掠面，且假設(shè)各致動(dòng)盤(pán)的誘導(dǎo)速度不同，如圖3所示。

風(fēng)速為v∞的無(wú)窮遠(yuǎn)來(lái)流經(jīng)垂直軸風(fēng)輪時(shí)會(huì)經(jīng)歷速度衰減過(guò)程，即流經(jīng)上風(fēng)區(qū)時(shí)，速度變成上風(fēng)區(qū)誘導(dǎo)速度V；流經(jīng)轉(zhuǎn)軸處時(shí)，速度進(jìn)一步降低至均衡誘導(dǎo)速度Ve；流經(jīng)下風(fēng)區(qū)時(shí)，速度為下風(fēng)區(qū)誘導(dǎo)速度V′，三者速度之間的關(guān)系為

V=uv∞

(27)

V′=u′(2u-1)v∞

(28)

式中，u、u′分別為上風(fēng)區(qū)誘導(dǎo)因子和下風(fēng)區(qū)誘導(dǎo)因子，u<1，u′<1。

對(duì)于垂直軸風(fēng)力機(jī)， 采取分風(fēng)區(qū)討論的方法

圖3 雙致動(dòng)盤(pán)多流管模型Fig.3 Model double-multiple stream-tubes

來(lái)研究葉片受力情況。如圖4所示，方位角θ位于0°～180°的區(qū)域稱(chēng)為上風(fēng)區(qū)；方位角θ位于180°～360°的區(qū)域稱(chēng)為下風(fēng)區(qū)。

圖4 葉片受力分析圖Fig.4 Analysis of blade stress

葉片所受的切向力FT和法向力FN分別為

(29)

式中，ρ為空氣密度；H為葉片高度。

當(dāng)葉片處于上風(fēng)區(qū)時(shí)，上風(fēng)函數(shù)為

(30)

式中，N為葉片數(shù)目；n為雙致動(dòng)盤(pán)多流管理論中的流管數(shù)目。

通過(guò)fupu=1-u可迭代計(jì)算出上風(fēng)區(qū)誘導(dǎo)因子u，進(jìn)而可計(jì)算出上風(fēng)區(qū)誘導(dǎo)速度、合成風(fēng)速，以及各方位下的實(shí)際攻角值，最后得到葉片在上風(fēng)區(qū)各個(gè)方位下的動(dòng)態(tài)氣動(dòng)力系數(shù)。

同理，當(dāng)葉片處于下風(fēng)區(qū)時(shí)，計(jì)算方法與上風(fēng)區(qū)類(lèi)似，下風(fēng)函數(shù)為：

(31)

將上風(fēng)區(qū)迭代得到的各個(gè)流管的誘導(dǎo)因子作為下風(fēng)區(qū)迭代公式fdwu′=1-u′的迭代初值進(jìn)行計(jì)算，便可得到下風(fēng)區(qū)各流管的誘導(dǎo)速度、合成風(fēng)速，以及各方位下的實(shí)際攻角值，最終得到葉片在下風(fēng)區(qū)各個(gè)方位下的動(dòng)態(tài)氣動(dòng)力系數(shù)。

4 算例求解與分析

以美國(guó)Sandia實(shí)驗(yàn)室17 m垂直軸風(fēng)力機(jī)[14-16]為算例進(jìn)行計(jì)算，其基本參數(shù)如表1所示。

表1 風(fēng)力機(jī)參數(shù)

將表1中的數(shù)據(jù)代入式(30)、式(31)，計(jì)算該垂直軸風(fēng)力機(jī)的葉尖速比分別為2.33和3.09時(shí)，葉片的動(dòng)態(tài)切向力系數(shù)和動(dòng)態(tài)法向力系數(shù)。

對(duì)比圖5可以發(fā)現(xiàn)，當(dāng)葉片處于上風(fēng)區(qū)且B-L修正模型和MIT修正模型在方位角0°～60°時(shí)，計(jì)算結(jié)果與實(shí)驗(yàn)數(shù)據(jù)吻合度較好；方位角60°～180°時(shí)，B-L模型的計(jì)算結(jié)果急劇減小，而MIT修正模型的計(jì)算結(jié)果基本反映實(shí)驗(yàn)數(shù)據(jù)的變化趨勢(shì)。當(dāng)葉片處于下風(fēng)區(qū)時(shí)，B-L修正模型的計(jì)算結(jié)果與實(shí)驗(yàn)數(shù)據(jù)相差不大，而MIT修正模型的計(jì)算結(jié)果除在方位角210°～240°、330°～360°能吻合外，其他范圍與實(shí)驗(yàn)結(jié)果的一致性較差，不推薦使用。

對(duì)比圖6可以發(fā)現(xiàn)，這兩種修正模型在計(jì)算動(dòng)態(tài)法向力系數(shù)時(shí)，總體上可以反映實(shí)驗(yàn)結(jié)果的變化規(guī)律，但在一些方位角下，兩種修正模型各有優(yōu)缺點(diǎn)：對(duì)于上風(fēng)區(qū)，B-L修正模型的總體趨勢(shì)預(yù)測(cè)精度要高于MIT修正模型；對(duì)于下風(fēng)區(qū)，B-L修正模型在方位角210°附近的預(yù)測(cè)值與實(shí)驗(yàn)結(jié)果相差較大，而MIT修正模型與實(shí)驗(yàn)結(jié)果的一致性較好。

(a)=2.33

(b)=3.09圖5 兩種葉尖速比的動(dòng)態(tài)切向力系數(shù)Fig.5 Dynamic tangential force coefficientsof 2 tip speed ratios

(a)=2.33

(b)=3.09圖6 兩種葉尖速比的動(dòng)態(tài)法向力系數(shù)Fig.6 Dynamic normal force coefficients of 2 tip speed ratios

基于上述對(duì)比分析，結(jié)合雙致動(dòng)盤(pán)多流管理論可以得到一種改進(jìn)的垂直軸風(fēng)力機(jī)葉片動(dòng)態(tài)失速性能計(jì)算模型。該模型計(jì)算動(dòng)態(tài)切向力系數(shù)時(shí)，上風(fēng)區(qū)使用MIT修正模型，下風(fēng)區(qū)使用B-L修正模型；計(jì)算動(dòng)態(tài)法向力系數(shù)時(shí)，上風(fēng)區(qū)使用B-L修正模型，下風(fēng)區(qū)使用MIT修正模型。修正模型的計(jì)算流程如圖7所示。

圖7 垂直軸風(fēng)力機(jī)葉片動(dòng)態(tài)失速性能計(jì)算修正模型Fig.7 Correction model for dynamic stall performance calculation of vertical axis wind turbine blades

5 結(jié)論

(1)結(jié)合垂直軸風(fēng)力機(jī)的非定常氣動(dòng)特性，對(duì)B-L和MIT動(dòng)態(tài)失速模型進(jìn)行了修正：簡(jiǎn)化了原B-L模型中的馬赫數(shù)相關(guān)項(xiàng)及翼型非繞流升力部分，改進(jìn)了流動(dòng)分離點(diǎn)的計(jì)算方法，將原模型中的攻角變化規(guī)律修正為垂直軸風(fēng)力機(jī)運(yùn)行時(shí)攻角的實(shí)際變化規(guī)律；將MIT模型中的特征風(fēng)速定義為垂直軸風(fēng)力機(jī)工作時(shí)的實(shí)際風(fēng)速，對(duì)動(dòng)態(tài)失速攻角及最大升阻力系數(shù)的計(jì)算公式進(jìn)行了修正，使該模型更加適合垂直軸風(fēng)力機(jī)葉片的實(shí)際工況。

(2)鑒于垂直軸風(fēng)輪流場(chǎng)的復(fù)雜性，采用分風(fēng)區(qū)的方法研究葉片的動(dòng)態(tài)受力。對(duì)比2種修正模型的計(jì)算結(jié)果和Sandia實(shí)驗(yàn)室17 m垂直軸風(fēng)力機(jī)實(shí)驗(yàn)數(shù)據(jù)發(fā)現(xiàn)：MIT修正模型對(duì)上風(fēng)區(qū)的切向力系數(shù)和下風(fēng)區(qū)的法向力系數(shù)的預(yù)測(cè)精度較高；B-L修正模型對(duì)上風(fēng)區(qū)法向力系數(shù)和下風(fēng)區(qū)切向力系數(shù)的計(jì)算結(jié)果與實(shí)驗(yàn)數(shù)據(jù)吻合良好。

(3)該模型在計(jì)算動(dòng)態(tài)切向力系數(shù)時(shí)，上風(fēng)區(qū)使用MIT修正模型，下風(fēng)區(qū)使用B-L修正模型；在計(jì)算動(dòng)態(tài)法向力系數(shù)時(shí)，上風(fēng)區(qū)使用B-L修正模型，下風(fēng)區(qū)使用MIT修正模型。