尾緣襟翼偏轉角對不同翼型的垂直軸風力機氣動影響研究

戴孟祎， 張志豪， 涂佳黃， 韓兆龍， 周 岱， 朱宏博

(1. 上海交通大學 船舶海洋與建筑工程學院， 上海 200240；2. 湘潭大學 土木工程學院， 湖南 湘潭 411105)

近年來，風能以其清潔的特性、豐富的儲量以及日趨成熟的風力發電技術，成為未來最具前景的新能源形式之一[1].風力機作為捕獲風能的重要裝置，其構型對于提高風力機的風能轉化率、降低生產和運行維護成本具有重要意義[2].與傳統水平軸風力機相比，垂直軸風力機具有設計簡單、對風向不敏感等優點[3]，但垂直軸風力機較低的風能利用率成為了其大規模商業化應用的瓶頸.

在旋轉過程中，垂直軸風力機葉片的攻角會出現周期性變化，當葉片的攻角超過失速攻角時，葉片表面氣體就會流動分離，從而使得葉片的氣動性能急劇下降[4].因此，研究人員提出利用不同的流動控制手段改善葉片的氣動性能.流動控制手段主要可分為被動控制[5-7]和主動控制[8-10]兩大類.其中，可變尾緣襟翼技術作為一種常見的主動控制手段，可通過改變翼型改變葉片的氣動性能.

Chen等[8]通過對比固定漿距和變槳距垂直軸風力機的性能，發現葉片在旋轉過程中做適當幅度的俯仰運動，能夠有效抑制功率輸出、轉速和轉矩輸出的波動，并提高風力機的功率系數.Li等[9]為探索最優槳距控制空間解，搭建由遺傳算法和計算流體力學(CFD)仿真模塊組成的變槳距自動優化平臺，得到能在較廣的葉尖速比范圍內使風力機平均功率系數均得以提升的優化葉片漿距.向斌等[10]為減小動態失速對垂直軸風力機氣動性能的影響，提出在翼型尾緣處布置動態格尼襟翼的主動控制方法.研究結果表明，采用主動式格尼襟翼提高了垂直軸風力機的運行穩定性，并增大風力機在低風速下的啟動力矩.繆維跑等[11]對NACA0012兩段式襟翼翼型進行數值模擬，利用弦線變換得到襟翼擺角與攻角的關系，發現襟翼翼型在靜態條件下由擺角引起的攻角遷移現象.研究結果表明，隨襟翼擺角增大，有效攻角范圍減小，翼型攻角產生遷移，受力亦發生變化.祖紅亞等[12]以NACA0018為基準翼型，采用數值模擬方法研究尾緣襟翼相對長度對翼型氣動性能的影響.研究結果表明，襟翼對翼型周圍主渦發展和變化具有影響，不僅改善了翼型的失速特性，同時也提高了翼型的氣動性能，襟翼翼型的失速攻角在研究范圍內均大于基準翼型.

上述研究均是針對某一特定翼型下襟翼對葉片或風力機氣動性能的影響所展開.然而，對于不同的翼型，襟翼對葉片和風力機氣動性能的影響存在一定差異.因此，針對不同翼型下帶有尾緣襟翼的風力機進行氣動性能分析和優化設計具有重要工程意義.

本文采用計算流體動力學數值模擬方法，建立添加分離式尾緣襟翼的H型垂直軸風力機二維模型，基于NACA0018、NACA0021和NACA0024共3種基礎翼型與-16°、-8°、0°、8°、16°共5組襟翼偏轉角參數，從風能利用率、葉片氣動力系數、風力機荷載和流場分布等方面對比分析了尾緣襟翼對不同翼型的H型垂直軸風力機的氣動性能影響，為進一步研究和提高H型垂直軸風力機氣動性能提供參考.

1 計算模型和方法

1.1 幾何模型

采用文獻[13]研究的三葉片H型垂直軸風力機作為原始模型，具體模型如圖1所示.該垂直軸風力機使用NACA0021對稱翼型，弦長c=85.8 mm，風輪旋轉直徑D=1 030 mm，轉軸直徑為c/6，自由來流速度U∞=9 m/s，葉片安裝角為0°.3個葉片在同一個平面內，葉片長度l=1.456 4 m, 互為120°夾角.風力機模型的展弦比(葉片長度與弦長之比)接近17，相對較大.根據Rezaeiha等[14]的結論，該風力機葉片的三維效應對葉片中部截面影響不顯著，適用于二維CFD模擬.因此，為減小計算量，本研究使用二維CFD模擬.為提高計算效率，忽略支撐桿對風力機整體性能的影響.

圖1 原始H型垂直軸風力機3D模型Fig.1 3D model of original H-type VAWT

圖2為風力機的旋轉示意圖.其中，α為葉片相對局部速度與弦線之間的夾角，即攻角；R為風輪旋轉半徑；ω′為角速度；FT為葉片所受切向力；FN為法向力；FL為升力；FD為阻力.假定順風區的自由來流速度與逆風區的自由來流速度相同[15]，即不考慮流動誘導的影響，相對局部速度可表示為

圖2 風力機旋轉示意圖Fig.2 Schematic diagram of wind turbine rotation

(1)

式中：λ為葉尖速比；θ為葉片所在方位角.

攻角可表示為

(2)

根據葉片所在方位角θ的不同，將風輪旋轉區域劃分為逆風區(45°≤θ<135°)、背風區(135°≤θ<225°)、順風區(225°≤θ<315°)和向風區(315°≤θ<45°)4個區域.功率系數Cp可由單葉片的彎矩系數Cm決定.

Cp=CmNλ

(3)

(4)

式中：N為葉片個數；M為風力機彎矩；ρ為空氣密度；AC為參考面積.

此外，選取切向力系數CT、法向力系數CN、推力系數Cth和側向力系數Cla共4個關鍵氣動力參數，對作用在葉片上的氣動力進行分析.

(5)

(6)

(7)

(8)

式中：H為葉片沿展向長度；Fth為葉片所受推力；Fla為側向力；A為風力機的掃風面積.在二維計算中，FT、FN、Fth和Fla可根據計算軟件監測到的葉片x和y方向上的力Fx和Fy與方位角θ的關系計算得到.

在原始模型的基礎上，建立添加分離式尾緣襟翼的垂直軸風力機幾何模型，葉片結構如圖3所示.葉片弦長與原始模型的弦長相同，為c=85.8 mm.將葉片前緣至襟翼旋轉中心的距離定義為主翼長度c1，襟翼旋轉中心至葉片尾緣的距離定義為襟翼長度c2，主翼與襟翼的間隙定義為翼縫寬度d.參照祖紅亞等[16-17]的工作，設置襟翼長度為0.3倍弦長，即c2=0.3c，設置翼縫寬度d=1.7%c.為了研究不同襟翼偏轉角對垂直軸風力機氣動性能的影響，襟翼偏轉角β可繞其旋轉中心在 ±16° 范圍內偏轉，定義偏轉角為正表示襟翼向上偏轉，偏轉角為負表示襟翼向下偏轉.

圖3 分離式尾緣襟翼葉片示意圖Fig.3 Schematic diagram of blade with separated trailing edge flap

1.2 數值方法

運用以有限體積法為基礎的商業軟件Fluent對垂直軸風力機氣動性能和流場特性進行分析.考慮到風力機的性能在很大程度上取決于葉片邊界層的發展，因此對轉捩點的準確預測至關重要.采用4方程轉捩剪切應力輸運(SST)湍流模型求解納維-斯托克斯(N-S)方程.轉捩SST湍流模型除了使用k-ωSST湍流模型中采用的湍流動能和耗散率方程外，還結合另外兩個考慮間歇因子γ和當地轉捩雷諾數Reθ t的輸運方程[18-19].其優點在于，一方面避免了一般情況下對平均場進行積分的過程，計算周期短且計算要求較低；另一方面對自由來流湍流度、分離和壓力梯度等影響轉捩的因素敏感，可以更準確地捕捉層流到湍流之間的過渡[20-21].

間歇因子表示在轉捩區某一固定點流動處于有脈動和無脈動狀態的時間比例，其輸運方程為

(12)

式中：t為時間尺度；μ為動力黏度系數；Uj為y軸方向的速度分量；xj為y軸方向的位移分量；μt為轉捩常數；S為應變率；Fle為轉捩長度函數；Fon和Ftu為轉捩控制函數；Ω為渦量；ca1、ca2、ce2和σγ為轉捩常數.

(14)

式中：混合函數Fθ t為邊界層中關閉源項；cθ t和σθ t為轉捩常數.

在垂直軸風力機非定常氣動特性的計算過程中，計算時間步長設置為T/1 440，其中，T為一個旋轉周期的時長.確保數值模擬過程中風力機轉子每轉動0.25°進行一個時間步計算.單次時間步長包含20次迭代，保證模擬的湍動能殘差ξ始終控制在ξ=10-6附近波動.基于Rezaeiha等[14]的研究結果，湍動能殘差控制為ξ=10-6可以兼顧模擬的時間復雜度和結果精確度.整個模擬過程設置為20個周期，總時長Tmax=2.72 s，前10個周期的初始化及流場發展可以保證模擬結果的充分收斂，第11～20個周期的模擬用于分析風力機氣動性能的相關參數.

2 計算域及有效性驗證

2.1 計算域和網格分布

計算流體動力學分析要求高精度流場，本文建立的二維模擬區域保證流場內流體能夠充分發展.為減小進出口邊界和側邊界對計算精度的影響[22]，將風力機布置在靠近上游的位置，入口到旋轉中心的距離為10D，旋轉域直徑為1.5D，計算域沿入流方向的長度為30D，寬度為20D.模型的阻塞比(在二維計算中定義為風輪直徑與計算域寬度之比)為5%，由對稱邊界條件導致的流動加速非常微弱，對風力機計算空氣動力學性能的影響可以忽略[23].在計算中，忽略連接桿件的影響，將整個流場分成4個區域：外流域F、轉軸控制域Z1、環形旋轉域Z2和葉片控制域Z3，各區域之間設置交界面，通過滑移網格技術實現風輪旋轉，具體計算域劃分如圖4所示.考慮葉片表面的粗糙度和擾流效應，葉片布置為無滑移壁面.邊界設置為前端速度入口U∞=9 m/s.為了統一模型參數的計算與分析，湍流強度設置為3%.空氣密度ρ=1.225 kg/m3，動力黏度μ=1.789 4×10-5Pa·s，壓力出口為0 Pa，基于弦長的雷諾數Re=5.28×104.

圖4 計算流場區域劃分Fig.4 Division of computational flow field region

為保證旋轉域與外流場的良好過渡和下游流場的充分發展與穩定，對外流域的網格以旋轉軸為中心進行十字形加密處理，加密寬度為3D，促進實現流場的進一步精確模擬.分別采用二維結構網格與非結構網格進行拓撲，轉軸控制域Z1和外流場域F采用結構網格，環形旋轉域Z2和葉片控制域Z3采用非結構網格進行網格加密，貼近壁面處使用邊界層網格，首層網格高度設置為 0.01 mm，網格增長率為1.2%，使葉片壁面處y+≈1，以滿足準確刻畫邊界層流動的要求.詳細網格布置如圖5所示，單個葉片控制域網格數約為17萬，總網格數約為68萬.所有仿真計算都在配置為Intel(R) Xeon(R) CPU E5-2676 V3@2.4 Hz的小型服務器上進行，一個計算周期大約需要消耗1 h，約20 h完成一個算例.

圖5 計算域的網格拓撲Fig.5 Mesh topology of computational domain

2.2 網格獨立性測試

網格質量通常對模擬結果的準確性有顯著影響.因此，在正式計算前，進行網格獨立性測試.改變葉片的最小表面尺寸，共設置3類網格方案，分別為粗網格、中網格和細網格.對葉尖速比λ=2.64時不同網格方案計算的功率系數進行比較，具體網格設置和結果如表1所示.

表1 基于功率系數誤差的網格獨立性測試Tab.1 Mesh independence test based on errors of power coefficient

結果表明，細網格和中網格方案的功率系數差別非常小，分別為0.354和0.350，但兩者的總網格數量相差較大，約19萬網格.另一方面，粗網格方案雖然網格數更少，但計算得到的功率系數的相對誤差較大.需要注意的是，網格數量越多并不意味著結果越準確，因為還會受網格形狀和分布的影響.因此，綜合考慮計算時長與計算精度的要求，本研究基于中網格方案來執行其他仿真.

2.3 模型有效性驗證

為評估CFD計算結果的準確性，選用文獻[13]的實驗結果與當前數值模擬結果進行對比驗證.該實驗數據已被廣泛用于研究垂直軸風力機的氣動性能[24-25]，計算結果與實驗數據對比如圖6所示.

圖6 現有結果與實驗和數值數據對比Fig.6 Comparison of current results with experimental and numerical data

可知，在葉尖速比λ達到額定葉尖速比2.64前，功率系數Cp隨λ增大而增大；而λ達到2.64后，功率系數Cp隨λ增加呈減小態勢.在低葉尖速比下，數值模擬的結果較好地表現了功率系數Cp隨葉尖速比λ變化的趨勢，而在高葉尖速比時，數值模擬結果對風力機的功率系數存在一定的高估，這與文獻[25]中的數值計算結果類似.在高葉尖速比下，Cp被高估主要有兩個原因：一是在二維CFD模擬中忽略了葉片高度，無法描述葉片兩端的局部流動特征，不存在葉片尖端損失；二是由于CFD模擬過程中的幾何簡化與建模，例如葉片輻條和連接至風力機塔架的支撐桿會導致較大的阻力，并降低風洞測量時的風力機性能，而這些在本研究的數值模擬中被忽略.相關研究表明，由于支撐臂摩擦而導致的風力機性能損失可以達到20%左右，在高葉尖速比下，這一損失將更明顯[13, 15].對于低葉尖速比(存在動態失速)，由于幾何簡化，同樣存在類似高估的問題.而在中等葉尖速比下，可能由于二維轉捩SST湍流模型無法準確模擬風力機氣動效應的固有三維流動復雜性，使這種預期的高估被抵消一部分.

總之，本文所采用的數值模型能夠較好地反映垂直軸風力機在風洞中的氣動特性，準確地再現實驗結果.因此，目前的數值模型可以作為后續研究的一種可靠模擬方法.

3 結果與分析

NACA對稱翼型族具有較高的升阻比和良好的失速特性，在H型垂直軸風力機的葉片選型上得到了廣泛應用.本文選用3種不同厚度的基礎翼型NACA0018、NACA0021和NACA0024，作為不同H型垂直軸風力機的葉片主體，葉片弦長與原始風力機模型的弦長一致.為分析尾緣襟翼對不同垂直軸風力機氣動性能的影響規律，針對添加分離式尾緣襟翼的3種基礎翼型NACA0018、NACA0021和NACA0024，分別選取襟翼偏轉角為 -16°、-8°、0°、8°和16°的5種情況，取原風力機模型的額定葉尖速比2.64進行計算模擬，討論風力機葉片的氣動性能與葉片附近區域的流場特性.

3.1 氣動性能與氣動荷載

葉片彎矩系數Cm是決定風能利用系數Cp的關鍵參數.圖7給出了翼型為NACA0018、NACA0021和NACA0024共3種垂直軸風力機在不同襟翼偏轉角下單個葉片的彎矩系數隨方位角的變化.可知，在逆風區，正向襟翼偏轉角能夠有效提高葉片的彎矩系數，且β=16°時，葉片的瞬時彎矩系數達到最優，此時NACA0018風力機的最大彎矩系數由0.179提高至0.246，增長了37.4%，NACA0021和NACA0024風力機的最大彎矩系數也分別增長了25.7%和30.4%.原因是在此方位角區間內，正向偏轉的襟翼增大了葉片的有效攻角，使葉片升阻比得到提升.同樣，在順風區，負向襟翼偏轉角對葉片的瞬時彎矩系數產生有利影響，且在β=-16° 時彎矩系數達到最優.考慮到中等葉尖速比下的垂直軸風力機本身氣動性能較優，并不存在較嚴重的動態失速問題[26]，因此可知，在其他區域，襟翼偏轉角的提升效果有限，甚至低于無偏轉時的風力機模型，對葉片彎矩系數造成不利影響.

圖7 不同襟翼偏轉角下風力機的單葉片彎矩系數Fig.7 Moment coefficient of blade at different flap angles

單個葉片的切向力系數和法向力系數對比分別如圖8和圖9所示.切向力方向與葉片繞旋轉軸旋轉的切線方向相同為正；法向力垂直于葉片繞轉軸旋轉的切線，指向風輪旋轉中心為正.與單葉片彎矩系數隨方位角的變化相似，逆風區的正向襟翼偏轉角和順風區的負向偏襟翼轉角可以使葉片的切向力系數得到提高，且β=16°時，葉片切向力系數在θ=90°處達到最大值.此時NACA0018風力機的最大切向力系數由1.757提高至2.532，增長了44.1%，NACA0021和NACA0024風力機的最大切向力系數也分別增長了30.6%和37.9%.

圖8 不同襟翼偏轉角下風力機的單葉片切向力系數Fig.8 Tangential force coefficient of blade with different flap angle

圖9 不同襟翼偏轉角下風力機的單葉片法向力系數Fig.9 Normal force coefficient of blade at different flap angles

從葉片的法向力系數對比可知，在整個旋轉周期范圍內，負向襟翼偏轉角能夠有效降低法向力系數.在β=-16° 時，葉片法向力系數在θ=90°處達到最大值.此時NACA0018風力機的最大法向力系數由10.310下降到7.202，降低了30.1%，NACA0021和NACA0024風力機的最大法向力系數也分別降低了24.6%和27.2%，表示風力機在旋轉過程中支撐桿受到的荷載可以得到大幅降低.

不同垂直軸風力機在一個穩定的旋轉周期內平均功率系數和平均彎矩系數的對比如圖10(a)和圖10(b)所示.顯然，分離式尾緣襟翼的偏轉角對風力機的動力性能有重要影響.對3種垂直軸風力機而言，在某些方位角區間內尾緣襟翼的偏轉能夠使風力機葉片的瞬時彎矩系數得到提高，但風力機整機的平均功率系數和平均彎矩系數均隨襟翼偏轉角的絕對值的增大而減小，說明簡單地設置固定偏轉襟翼并不能使整機動力性能得到明顯改善，這與圖7分析結果一致.

對于大型垂直軸風機的設計與運用，葉片周期性的轉動將導致整機受到周期性的疲勞荷載，這些疲勞荷載是需要重點考量的因素.風力機的氣動載荷是衡量其氣動性能的主要參數，圖10(c)和10(d)分別為不同垂直軸風力機在一個穩定的旋轉周期內平均推力系數和平均側向力系數的對比，推力和側向力分別指風力機整體受到的沿來流方向的力和垂直于來流方向的力.可知，負向襟翼偏轉角會在一定程度上降低風力機的推力系數；而偏轉角為正時，推力系數隨偏轉角度的增大而增大.對于側向力系數，整體呈現為隨偏轉角絕對值的增大而增大的趨勢.說明單一地改變襟翼偏轉角對改善風機整體荷載的作用效果不突出，與圖8和圖9的分析結果相符.

圖10 不同垂直軸風力機的對比Fig.10 Comparison of different VAWTs

為進一步分析尾緣襟翼對不同厚度翼型的垂直軸風力機動力性能的影響，將各風力機在不同襟翼偏轉角下的功率系數與β=0° 時的功率系數作差得到ΔCp，并求功率系數的變化率，即ΔCp/Cp(β=0°)，可直觀得到不同風力機ΔCp隨尾緣襟翼的變化分布，如圖11所示.對任意風力機而言，襟翼偏轉角的絕對值越大，功率系數變化越明顯，且變化速率逐漸增大.對比3種垂直軸風力機，當襟翼偏轉角為負時，NACA0024風力機的功率系數隨偏轉角變化最明顯，結合圖10(a)可知，偏轉角從0°減小至 -16°，其功率系數由0.301降低至0.129，下降了57.1%，NACA0021和NACA0018風力機的功率系數也分別下降了50.3%和39.2%；當襟翼偏轉角為正時，NACA0018風力機的功率系數隨偏轉角變化最明顯，偏轉角從0° 增大至16°，其功率系數由0.352降低至0.054，下降了84.7%，NACA0021和NACA0018風力機的功率系數也分別下降了65.3%和55.3%.由此可見，在負向襟翼偏轉角下，風力機功率系數受偏轉角影響的程度與翼型厚度呈正相關；而在正向襟翼偏轉角下，功率系數受偏轉角影響的程度與翼型厚度呈負相關.

圖11 不同垂直軸風力機ΔCp隨尾緣襟翼變化分布Fig.11 Distribution of ΔCp with trailing edge flap deflection angle for different airfoils

3.2 流場分布

動態失速與葉片表面發生的流動分離直接相關[27].圖12給出了NACA0018和NACA0024風力機在不同襟翼偏轉角下翼型周圍的渦量分布對比，以觀察流動分離的形成和發展，并進一步解釋風力機氣動性能特性的內在機理.

圖12 不同風力機翼型周圍的渦量分布對比Fig.12 Comparison of vorticity distributions around airfoil for different VAWTs

在方位角θ=60° 時，對于β=0° 的模型，NACA0018風力機的氣流基本附著在葉片表面，未發生流動分離現象，而NACA0024風力機尾緣下游區域存在顯著的周期性旋渦結構；對于β=16° 的模型，兩種不同翼型工況下，襟翼的偏轉導致襟翼處產生明顯的流動分離，尾緣處形成周期性的卡門渦街.當葉片旋轉到θ=90° 位置時，對于β=0° 的模型，由于攻角的增大，尾緣產生輕微流動分離，形成周期性尾渦.同時由于翼縫的存在，襟翼靠近翼縫的位置出現小尺度的渦；對于β=16° 的模型，襟翼的偏轉延緩了尾緣處的流動分離，并使翼縫處的渦得到小范圍擴散，解釋了圖7所示的葉片產生更大彎矩系數的結果.當葉片旋轉到θ=120° 位置時，由于攻角增大，葉片表面的氣流高度分離.后緣渦和前緣渦的強度大大增加，表明在這一階段發生了嚴重的動態失速，對應圖7中彎矩系數的急劇下降.在θ=150° 位置，由于攻角逐漸減小，從前緣產生的氣流分離逐漸受到抑制，后緣產生的旋渦結構向下游處進一步發展.對于β=16° 的工況，葉片襟翼的偏轉加劇了后緣附近處的流動分離.當葉片旋轉到θ=180° 位置時，從葉片前緣和后緣脫落的渦逐漸消散，攻角在此位置再次歸零，氣流逐漸重新附著到葉片表面，流動分離強度減弱，因而導致圖7中彎矩系數激增.對比不同方位角下，NACA0018和NACA0024風力機有無襟翼偏轉工況下的渦量圖可知，β=16° 的偏轉襟翼對NACA0018風力機流場分布的影響比對NACA0024風力機的影響更顯著，尤其在θ=120° 和θ=150° 的位置，這與圖11得到的分析結果一致.

4 結論

針對添加尾緣襟翼的垂直軸風力機的空氣動力學問題，采用計算流體動力學轉捩SST湍流模型，開展3種不同分離式尾緣襟翼的翼型(NACA0018、NACA0021和NACA0024)葉片的H型垂直軸風力機氣動性能的數值研究.對比不同襟翼偏轉角 -16°、-8°、0°、8°、16° 下3種風力機模型的風能利用率、葉片氣動力系數、風力機荷載和流場分布等結果，主要結論如下：

(1) 襟翼偏轉角度變化能夠導致風力機氣動性能變化.在本文的研究模型中，與尾緣襟翼不發生偏轉的模型相比，風力機的平均功率系數隨襟翼偏轉角的絕對值增大而減小.

(2) 在逆風區(45°≤θ<135°)，正向襟翼偏轉角可以有效提高葉片的彎矩系數，當β=16° 時，葉片的瞬時彎矩系數達到最優，NACA0018、NACA0021和NACA0024風力機的彎矩系數分別增長了37.4%、25.7%和30.4%；在順風區(225°≤θ< 315°)，負向襟翼偏轉角對葉片的瞬時彎矩系數產生有利影響，在β=-16° 時彎矩系數達到最優.

(3) 對不同翼型的H型垂直軸風力機，當襟翼偏轉角為負時，風能利用率受偏轉角影響的程度與翼型厚度呈正相關，即NACA0024風力機受影響程度最大，在β=-16° 時風能利用率下降了57.1%；當襟翼偏轉角為正時，風能利用率受偏轉角影響的程度與翼型厚度呈負相關，即NACA0018風力機受影響程度最大，在β=16° 時風能利用率下降了84.7%.