朱新宇，郭志平，楊燕昭，李慶安

(1.內(nèi)蒙古工業(yè)大學機械工程學院，內(nèi)蒙古 呼和浩特 010051；2.哈爾濱工業(yè)大學能源科學與工程學院，黑龍江 哈爾濱 150001；3.中國科學院工程熱物理研究所，北京 100190)

1 引言

隨著日益嚴峻的能源危機，新能源的開發(fā)與利用得到了足夠的重視，水平軸風力發(fā)電機也得到了快速發(fā)展。然而，水平軸風力發(fā)電機受湍流強度的影響較大，大多安裝在距離城市中心和用電工業(yè)區(qū)較遠的高山、海洋和平原等地區(qū)。所以造成輸電距離遠、建造成本高和輸電損耗大等缺點。因此，受湍流強度影響較小的垂直軸風力發(fā)電機得到科研人員的廣泛關(guān)注。在風力發(fā)電機的大家族中，垂直軸風力發(fā)電機是被最早發(fā)明的，但發(fā)展程度卻低于水平軸風力發(fā)電機，這主要由垂直軸風力發(fā)電機的相對風速和攻角隨回轉(zhuǎn)角的變化而變化，導致其流場復(fù)雜多變，其計算和研究都十分困難。在計算機性能較低的年代，僅靠人力計算其一次流場，將花費一年以上時間。因此，至今垂直軸風力發(fā)電機都沒有一套統(tǒng)一的標準。

盡管對于垂直軸風力發(fā)電機的設(shè)計還不成熟，但近些年來，在科研人員的不懈努力下也取得了一些成果。例如，Li et al.[1]在風洞試驗中，利用激光多普勒測速儀(Laser Doppler Velocimetry，LDV)通過局部風速在三個不同的葉尖速比下，測量并計算了垂直軸風力發(fā)電機的相對風速和攻角。Bhargav[2]等人應(yīng)用CFD仿真技術(shù)研究脈動風速對于垂直軸風力機的影響。楊燕昭等人[3]在CFD仿真中，采用不同湍流模型分析了垂直軸風力發(fā)電機的力矩系數(shù)和葉片壓力系數(shù)，得出k-ω(SST)湍流模型模擬效果較好。Yang Y[4]等人通過實驗與CFD仿真技術(shù)研究垂直軸風力機葉片最佳安裝角。Yang Y[5]等人基于Q準則，并使用CFD仿真技術(shù)對垂直軸風力機產(chǎn)生的渦與渦對尾流場的影響做了詳細分析。Haitian Zhu[6]等人采用CFD仿真技術(shù)研究實度對垂直軸風力發(fā)電機性能的影響。由于空氣流動對風力發(fā)電機的性能影響較大，研究人員大多關(guān)注葉片周邊的流場，為了進一步提升整體垂直軸風力發(fā)電機的氣動性能，并為建立垂直軸風力發(fā)電廠建立基礎(chǔ)，理清垂直軸軸風力發(fā)電機整體流場特性成為了研究風力發(fā)電機不可缺少的一環(huán)。

2 模型

2.1 物理模型

圖1為垂直軸風力發(fā)電機風輪的簡化模型的軸測圖，此模型包括兩枚直線翼葉片、四根葉片支架和一根回轉(zhuǎn)軸。風輪直徑為D=2.0m，翼展高度為1.2m，葉片翼型采用標準對稱翼型NACA0021，其弦長為0.265m。圖中笛卡爾坐標系為風輪的絕對坐標系，x軸為自由風速方向，y軸為垂直于來流方向，z軸為葉片翼展方向。圖2為垂直軸風力發(fā)電機風輪的俯視圖，其回轉(zhuǎn)方向為順時針，當葉片處于風輪最下端時，葉片的方位角為θ=0。自由來流風速U從左到右吹過風輪。為了后續(xù)更好的分析風輪的流場特性，將其流場分為上流域(x<0)和下流域(x>0)。

圖1 垂直軸風力發(fā)電機風輪的簡化模型

圖2 垂直軸風力發(fā)電機風輪的俯視圖

2.2 湍流模型

建立的風力發(fā)電機模型按照(1)式來計算

(1)

式中，Re為雷諾數(shù)，U0為來流風速(m/s)，c為翼弦長(m)，v為運動粘性系數(shù)(m2/s)。

可以得到：求出的雷諾數(shù)均大于臨界雷諾數(shù)。由此可知風力發(fā)電機周圍流場的流動狀態(tài)以湍流為主[7]。

因為k-ω SST湍流模型具有良好的穩(wěn)定性，收斂性和對自由來流的湍流度不敏感性，并且能夠?qū)毫μ荻攘鲃拥膶?shù)層做出精確預(yù)報等特性。所以本次采用的是雷諾時均湍流模型中的k-ω(SST)模型。

該模型假定湍流粘性為μt，則有

(2)

式中，ρ為空氣密度，k為湍流動能，ω為湍流頻率。

湍流動能k方程

(3)

式中，μ為時均速度，μ1為湍動粘度，δ是“Kronecherdelta”符號。

湍流頻率ω方程

(4)

因為Menter[8]認為原始k-ω模型沒有考慮湍流剪應(yīng)力的輸運，這會導致對于渦粘性的過分估計，因此提出對公式提出了對渦粘性進行限制的改進

(5)

2.2 數(shù)值模型

圖3(a)(b)和(c)展示了CFD仿真的數(shù)值模型，圖(a)為CFD仿真的整體模型，其模型沿x軸方向的長度為20D，沿y軸方向的寬度為10D，沿z軸方向的高度為1D。圖(b)展示了數(shù)值域風輪附近的網(wǎng)格，由于風輪附近的風速和壓力變化梯度較大，對風輪附近的網(wǎng)格進行了加密處理。此數(shù)值模型采用已被廣泛應(yīng)用在風機CFD仿真領(lǐng)域的滑移網(wǎng)格技術(shù)。因此，此數(shù)值模型分為動網(wǎng)格(move mesh)和靜網(wǎng)格(static mesh)。圖(c)展示了葉片附近的被進一步加密網(wǎng)格和邊界層網(wǎng)格，其中雷諾數(shù)Re=2.89×105，根據(jù)文獻[9]可知邊界層的y+<1，邊界層總厚度為2×10-5m，增長因子為1.25。

圖3 CFD仿真的數(shù)值模

入口采用速度入口邊界條件，入口風速U0=8m/s，湍流強度5%；出口采用壓力出口邊界條件；回轉(zhuǎn)區(qū)域交界面采用interface邊界條件，回轉(zhuǎn)區(qū)域速度滿足葉尖速比λ=(ωR)/U0=2.19；回轉(zhuǎn)軸、葉片以及葉片支架采用無滑移壁面條件。

3 結(jié)果分析

3.1 單葉片功率系數(shù)曲線

為了驗證CFD仿真的有效性，將CFD仿真所得的功率系數(shù)與來自參考文獻[10]的實驗數(shù)據(jù)進行對比。參考文獻[10]中的實驗數(shù)據(jù)是利用多點壓力測量儀在風洞實驗所得。其對比曲線如圖4所示。圖中橫坐標與縱坐標分別代表方位角和單個葉片的功率系數(shù)，其曲線為葉片回轉(zhuǎn)一周所得。在方位角為50°<θ<160°的上流域區(qū)間，由CFD仿真所得數(shù)據(jù)與風洞實驗所得數(shù)據(jù)能夠較好的擬合。然而，在其它方位角處CFD仿真所得曲線略小于由風洞實驗所得數(shù)據(jù)。其原因是此時葉片主要處于風輪的下流域，由于上流域紊亂氣流的影響使得下流域的氣流非常復(fù)雜，以目前的數(shù)值仿真模型很難對其進行精準計算。因此，當葉片轉(zhuǎn)到下流域時，CFD仿真對其計算會出現(xiàn)一定偏差。由圖4還可得，CFD仿真所得曲線峰值對應(yīng)的方位角比由風洞實驗所得大，說明CFD仿真具有一定的延時性。由上可知：CFD仿真與實際實驗是有一定偏差的，但該偏差對于本次利用CFD仿真研究垂直軸風力發(fā)電機的流場特性影響較小。

圖4 單葉片功率系數(shù)曲線

3.2 垂直于z軸的風速云圖

由于風輪的能量轉(zhuǎn)換和風輪對氣流的阻擋效應(yīng)，使得垂直軸風力發(fā)電機的下流域出現(xiàn)一個寬大的低風速區(qū)，這個低風速區(qū)的長短直接影響處于其下流域風機的工作效率，這對于集群風力發(fā)電機的密度具有重要影響。圖5描述了垂直軸風力發(fā)電機下流域風速云圖。此云圖是風輪的俯視圖，其截面位于z=0處，圖5的(a)、(b)、(c)和(d)分別代表θ=0°、θ=45°、θ=90°和θ=135°處的風速云圖。

圖5 垂直于z軸風速云圖

由圖5可知：風輪的低風速區(qū)域不是關(guān)于x軸對稱分布，其低風速區(qū)域偏向于方位角270°<θ<90°側(cè)(即圖2中y的負半軸)，出現(xiàn)這個現(xiàn)象的原因主要是風輪的回轉(zhuǎn)方向所致。沿著x軸方向，低風速區(qū)域面積逐漸增大，但是風速值大小在逐漸恢復(fù)。處于上流域的葉片對下流域的氣流產(chǎn)生很大影響，使得處于下流域的葉片周圍的風速比處于上流域葉片周圍的風速低很多，這對處于下流域葉片的氣動性能降低。在圖5中，當方位角θ=0°時風輪的尾流低風速區(qū)域最長；當方位角θ=45°時，風輪的尾流低風速區(qū)最短。說明當θ=45°時，下流域的低風速區(qū)域恢復(fù)的更快。從圖5中還可看出葉片在135°<θ<270°時葉片周圍跟隨一個低風速區(qū)域，這是由于葉片旋轉(zhuǎn)方向與風速方向相互作用，葉片尾流會跟隨葉片，這是氣動性能下降的原因之一。

3.3 垂直于x軸的風速云圖

圖6描述了垂直于x軸的風速分布云圖。其取自沿x軸方向不同位置處的風速云圖，圖中(a)、(b)、(c)、(d)分別代表了方位角θ=0°、45°、90°、135°時垂直于x軸截面的風速云圖。

圖6 垂直于x軸風速分布云圖

如圖6所示，中心區(qū)域風速值最低，隨著x/R的增大，低風速區(qū)的面積逐漸擴散，但是由于低風速區(qū)與自由流之間能量交換和相互作用，使風速值在逐漸增加。當θ=0°和θ=135°時，在x/R=1.0處出現(xiàn)風速值大于8m/s的區(qū)域，這是由于此處低風速區(qū)域還沒有得到充分擴散，同時葉片對外部氣流產(chǎn)生擠壓而引起的。由圖6還可得，所有云圖關(guān)于絕對坐標系z=0的平面對稱。

3.4 風速值曲線

為了定量的分析風數(shù)值的規(guī)律，所以建立風速值曲線。圖7中的(a)、(b)、(c)和(d)分別代表了方位角θ=0°、θ=45°、θ=90°和θ=135°時沿x軸方向風速大小的分布。x軸表示局部風速值，y軸表示沿絕對坐標系y軸的橫向方向。其帶正方形特征的曲線、帶圓形特征的曲線、帶上三角形特征的曲線、帶下三角形特征的曲線和帶菱形特征的曲線分別代表x/R=-2、0.5、1、4和8處風速大小的分布。在y軸方向上，取-2

圖7 風速值曲線圖

3.5 平均風速值曲線

為了展示了風輪下流域各方位角的平均風速值曲線，在風輪回轉(zhuǎn)過程中，每轉(zhuǎn)5°取一次風速值，再求轉(zhuǎn)過一周風速的平均值。

如圖8所示，在x/R=-2時，此處位于風輪的上流域，可以看到此時的風速略小于來流風速，產(chǎn)生這個現(xiàn)象的原因是風機回轉(zhuǎn)形成的氣流所影響的；當x/R=0.5時，在y=0處，出現(xiàn)一個風速突變值。其原因是由于風輪回轉(zhuǎn)軸對流場產(chǎn)生影響而造成的；當x/R=4.0時，風速值最??；當x/R=8時，其風速值恢復(fù)到來流風速的四分之三。

圖8 平均風速值曲線圖

4 結(jié)論

通過對垂直軸風力發(fā)電機下流域的流場特性進行分析，可得出以下結(jié)論：

1)CFD仿真所得數(shù)據(jù)與風洞試驗數(shù)據(jù)有一定的偏差，但對于本次用CFD仿真研究垂直軸風力發(fā)電機流場特性影響較小。

2)低風速區(qū)域不關(guān)于x軸對稱，其偏向于270°<θ<90°，但是關(guān)于z軸對稱，所以在考慮風力發(fā)電機分布式應(yīng)交錯排布，預(yù)計交錯式排布可以增加風場的風能利用。

3)隨著x/R的增大，低風速區(qū)的面積逐漸擴散，其值先減小后增大。在x/R=4處，風速值達到最小值，在x/R=8處，風速值恢復(fù)到初始風速的75%。所以在風場中上一臺風力發(fā)電機距下游的風力發(fā)電機間隔應(yīng)大于4R且間隔最好接近8R。