曹 宇， 韓兆龍,b， 周 岱,b,c， 雷 航

(上海交通大學 a. 船舶海洋與建筑工程學院； b. 海洋工程國家重點實驗室上海； c. 高新船舶與深海開發(fā)裝備協同創(chuàng)新中心， 上海 200240)

垂直軸風力機因控制系統(tǒng)簡單、整體受力合理等優(yōu)勢在海上風電產業(yè)中應用廣泛[1].但是，復雜的海洋風環(huán)境，如風速過大、風向多變會導致垂直軸風力機的轉子產生極大的轉矩，造成底部浮式平臺失穩(wěn).對此，學者們對海上垂直軸風力機的穩(wěn)定性進行了研究.徐應瑜等[2]在單柱式平臺的底部附加3根約束錨鏈，研究了10 MW浮式風機在極端風浪條件下的運動特性.劉中柏等[3]通過控制垂蕩板的邊長和面積，研究了海上風電浮式基礎的運動特性.但是，以上方法均會增加海上風電的工程成本.同軸對轉式垂直軸風力機是在同一支撐軸上設置上、下共2個轉動方向完全相反的風力機轉子，通過結構優(yōu)化，在保證浮式平臺穩(wěn)定性的基礎上，控制工程的成本.

目前，針對同軸對轉式垂直軸風力機的研究較少，本文的部分研究借鑒Lam等[4]關于同平面對轉式垂直軸風力機的成果：該類型風力機的結構的穩(wěn)定性更高，尾端脫落渦消散更快，風能利用率更高.基于該成果，采用計算流體力學(CFD)數值模擬方法，對比在不同的葉尖風速比(TSR)下，對轉式與獨立式垂直軸風力機的空氣動力響應情況，包括轉子的轉矩和風力機的風能利用效率以及渦量場分布情況，進而研究對轉式垂直軸風力機的氣動和穩(wěn)定性能.

1 模型和數值方法

1.1 計算模型

基于計算精確度和計算成本的要求，選擇雷諾時均剪切應力傳輸RANS SSTk-ω模型作為湍流模型模擬對轉式垂直軸風力機的氣動性能.RANS SSTk-ω模型可以通過時均化處理將求解復雜的Navier-Stokes方程轉化為求解簡單的平均流方程.此外，在求解渦旋黏性時，用比耗散率ω代替耗散率ε，使得模型在低雷諾數、過渡段和高雷諾數時均有較好的結果[5]，保證了結果的精確度.RANS SSTk-ω湍流模型由湍動能k和ω兩項輸運方程組成.

(1)

(2)

(3)

(4)

(5)

(6)

式中：vx、vy、vz分別為x、y、z方向的速度分量；t為來流時間；μ為運動黏度；μT為渦旋黏度；σk為k的湍流普朗特數；Pk為k方程的分項；β*、α和β為參數項，且β*=0.09，α=1，β=1.5；S為廣義源項；σω為ω的湍流普朗特數；F1為ω方程的協調方程.方程的具體求解過程見文獻[6].

圖1 對轉式垂直軸風力機三維模型圖Fig.1 A 3D model sketch of counter-rotating vertical axis wind turbine

圖2 對轉式垂直軸風力機轉子的平面俯視圖Fig.2 Top view of rotor of counter-rotating vertical axis wind turbine

1.2 幾何模型

圖1所示為對轉式垂直軸風力機的三維幾何模型，風力機包括同軸的上、下轉子，上轉子順時針轉動，下轉子逆時針轉動，并利用中心支撐將葉片與底部浮式平臺連接.圖2為上、下轉子的平面俯視圖，葉片的安裝角(弦線與切線的夾角)β′=6°.轉子直徑D=2 m，高度為1.32 m；H型葉片翼型采用NACA0021，葉片高度H=1.2 m，弦長為0.265 m.由于連接桿和中心支撐對風力機氣動和力學性能的影響有限，所以本文模型不包括風力機連接桿和中心支撐[7].

1.3 計算域

利用CFD數值模擬方法和SSTk-ω湍流模型，研究在風場環(huán)境下的對轉式垂直軸風力機的空氣動力響應，并建立計算域，如圖3所示.圖3(a)中，計算域前端為速度進口，來流速度v∞=8 m/s；后端為壓力出口，壓力大小為一個標準大氣壓.由于滑移壁面對流場的摩擦阻力為0，流速和渦量分布固定不變，考慮數值模擬中計算域的尺寸遠小于場地試驗的流場環(huán)境，且計算域外的風速和渦量分布與壁面的情況不同，所以將計算域的4個壁面均設置為滑移壁面，可以更真實地反映實際場地試驗的流場情況，保證模擬的準確性.圖3(b)和(c)中，計算域包括旋轉域和外流場域.旋轉域的直徑為3.2 m，高度為 1.32 m，上、下轉子的間距為0.265 m，外流場域的尺寸為 31 m×12 m×10 m.計算模型中用重疊網格模擬轉子的旋轉特性，外流場域為背景區(qū)域，旋轉域為重疊區(qū)域，兩者之間的邊界為重疊界面.每1個時間步內，背景區(qū)域和重疊區(qū)域通過界面交換數據，

圖3 計算域和邊界的參數說明Fig.3 Parameters of computational domain and boundary

因此，重疊界面上2個區(qū)域的網格始終一致，確保了數據交換的準確性.旋轉中心軸在x方向距離速度進口6 m，距離壓力出口25 m，以保證模擬的準確性和尾渦擴散的完整度[8].

1.4 網格劃分

計算網格的拓撲結構如圖4所示.為了加快界面數據的交換速度，提高數值模擬的效率，模型的網格劃分采用切割體網格形式：圖4(a)中，重疊網格背景區(qū)域的外流場域基礎網格尺寸設置為 0.5 m，最小網格尺寸設置為0.25 m；重疊網格重疊區(qū)域的旋轉域基礎網格尺寸設置為0.05 m，最小網格尺寸設置為0.025 m；重疊邊界的網格尺寸設置為 0.05 m.外流場域與旋轉域的網格尺寸相差較大，為了保證2個區(qū)域網格的平穩(wěn)過渡，在外流場域和旋轉域之間外加1個網格尺寸為0.2 m的圓柱.同時，為了研究旋轉域后側遠端的尾渦發(fā)展過程，對外流場的尾渦區(qū)域進行加密處理，加密區(qū)的網格尺寸為 0.2 m.圖4(b)為葉片的網格劃分情況，旋轉域網格的基礎尺寸為0.05 m，葉片邊緣的網格尺寸為 0.01 m.旋轉域到葉片邊緣的網格呈逐漸加密的趨勢，且由于葉片前緣的曲率較大，為保證計算結果的精確性，前緣網格尺寸更小，為0.005 m.為了精確分析風力機的氣動性能和葉片上渦的情況，在模型的葉片表面附加棱柱層網格，如圖4(c)和(d)所示.棱柱層的總厚度為0.005 m，增長率為120%，層數為30.基于模型的雷諾數Re=1.293×105，模型首層棱柱層厚度為3×10-5m，首層網格與壁面間距的無量綱參數

y+=u*y/ν<1

式中：u*為近壁面摩擦速度；y為首層網格與壁面的間距；ν為運動黏度.以保證葉片邊緣附著渦的穩(wěn)定性和模擬的準確性[9].

圖4 計算網格的拓撲結構圖Fig.4 Topology structure of computational grids

1.5 求解器設置

采用隱式非定常流動方法求解連續(xù)性方程，采用二階時間差分法求解動量方程.在非定常計算中，將求解器的時間步長設置為0.002 s，保證數值模擬在轉子每轉動2° 時進行1個時間步計算；1次時間步包含10次迭代，保證模擬的湍動能殘差始終在10-5左右波動[10].整個模擬過程設置為10個周期，總時長為3.74 s，以保證模擬結果收斂完全.

2 模型驗證和網格獨立性分析

2.1 驗證方法

通過對比數值模擬與試驗得到的風能利用效率(Cp)，驗證風力機模擬的精確性.Cp是衡量風力機氣動性能的關鍵因素，計算公式為

其中：Q為風力機轉矩；Ω為轉子角速度；ρ為來流密度.

圖5為不同TSR時，CFD數值模擬與風洞試驗[11]的Cp對比曲線和相對誤差曲線.其中：

TSR=ΩD/2v∞

圖5(a)中，當TSR<1.80時，風洞試驗與數值模擬的Cp值的相似度較高；當TSR>1.80時，Cp值的相似度較低，在曲線的最高點處，數值模擬曲線偏左的誤差率達6.7%.這是由于本模型的CFD數值模擬不考慮包含在風洞試驗中的塔架、支撐等連接件，同時模擬的結果忽略了機械摩擦耗能的影響，所以數值模擬的Cp值較大.即在相同的Cp下，CFD數值模擬對應的TSR更小，對應的曲線偏左.當 TSR=2.10時，數值模擬和風洞試驗的Cp值均達到最大，約為0.18.圖5(b)中，當TSR>1.30時，Cp的相對誤差均小于12%， 模擬結果與試驗結果吻合較好；當TSR=1.05時，Cp相對誤差最大，約為28.5%.這是由于當TSR=1.05時，風力機處于動態(tài)失速狀態(tài)，葉片邊緣出現了大量脫落渦，而風力機葉片、連接桿和中心支撐等構件上的脫落渦會導致葉片升力和風力機轉矩的降低，從而引起Cp值的降低.CFD數值模擬忽略了連接件、中心支撐的摩擦耗能的影響，未考慮忽略構件上的脫落渦影響，導致Cp值偏大，精度出現較大偏差.綜上可知，當TSR=2.10時，Cp值最大，且Cp相對誤差約為0.因此，后文數值模擬的分析均基于TSR=2.10的情況.

圖5 CFD數值模擬與風洞試驗的Cp值對比Fig.5 Comparison of Cp values of CFD numerical simulation and wind tunnel test

2.2 網格獨立性分析

表1為TSR=2.10時， 基于Cp相對誤差的網格獨立性驗證.隨著網格由粗糙變精細(網格尺寸由大變小)，Cp值增大，Cp相對誤差減小，邊界層的y+值減小.在精細網格中，風洞試驗的Cp相對誤差和y+值均最小.網格劃分越細，Cp相對誤差越小，y+越小，模擬結果的收斂性和精確性越好[12].但是，隨著網格由中等變精細，總網格數增加了37%，Cp相對誤差僅降低了2%.因此，綜合考慮模型的計算成本和精確度要求，后文研究均基于中等網格的情況.

表1 網格獨立性驗證Tab.1 Grid independence verification

3 分析與討論

3.1 轉矩對比

圖6為TSR=2.10時，對轉式與獨立式垂直軸風力機的總轉矩(Qt)隨來流時間的對比關系.綠色框中為一個周期內總轉矩的變化規(guī)律[13]，當風力機的葉片位于上風區(qū)且方位角θ=60°～120° 時，葉片升力和風力機轉矩取得最大值.而本文的研究對象為兩葉片垂直軸風力機，在同一周期內，2個葉片會分別轉動到上風區(qū)θ=60°～120° 的位置.因此，1個周期內風力機轉矩會取得2次最大值，且兩值的大小相近.由圖可知，獨立式風力機轉矩的最大值為30 N·m，最小值為-10 N·m，兩值相差40 N·m.而對轉式風力機的總轉矩曲線在Qt=0附近小幅度波動，最大值與最小值相差小于3 N·m.可知，同等風場條件下，對轉式風力機傳遞到浮式平臺的的總轉矩更小，結構的穩(wěn)定性更強.這是由于對轉式風力機存在完全相同的上、下轉子，在對轉過程中，轉子傳遞到浮臺的轉矩近似為大小相等、方向相反，即轉矩幾乎抵消，所以總轉矩在0附近波動且波動幅度很小.

圖6 風力機的轉矩對比Fig.6 Comparison of torques of different wind turbines

3.2 氣動性能分析

圖7為不同TSR的對轉式與獨立式垂直軸風力機的Cp對比關系.其中，Cp1、Cp2分別為對轉式風力機上、下轉子的風能利用效率，Cps為獨立式風力機的風能利用效率.圖中，Cp1和Cp2曲線基本重合，相對誤差小于1%.這是由于對轉式風力機的轉子為反對稱結構，上、下轉子的氣動性能相似，流場分布和脫渦情況也基本一致，所以Cp值基本相同.而Cp1、Cp2與Cps的差異明顯，相對誤差大于10%.當TSR=1.30時，對轉式與獨立式風力機的Cp曲線出現1個交叉點，交叉點兩側的Cp特征完全相反.當TSR<1.30時，Cps大于Cp1和Cp2，最大差值點在TSR=1.05處，相對誤差約為55.8%.當TSR>1.30時，Cps小于Cp1和Cp2，最大差值點在TSR=2.10處，相對差值約為10.85%.

圖7 轉子的Cp值對比Fig.7 Comparison of Cp values of different rotors

圖8為不同TSR的葉片攻角α與θ的關系曲線，圖9為獨立式和對轉式垂直軸風力機的三維渦量圖.基于圖8和9，對圖7中TSR<1.30時，Cps大于Cp1和Cp2的情況進行分析(取TSR=1.05時的風況).

圖9 三維渦量圖Fig.9 3D vorticity diagrams

圖8中α與θ為正弦關系，隨著TSR的增大，α曲線的幅值減小.風力機的失速情況在α>15° 時發(fā)生[14]，圖中水平虛線為α=15° 的情況，點A、A′為TSR=1.05時水平虛線與關系曲線的交點，點B、B′為TSR=2.40時水平虛線與關系曲線的交點.可以看出，線段AA′明顯長于線段BB′，即當TSR=1.05、α>15° 時，風力機對應的方位角的范圍更大，失速時間更長.因此，TSR較小時，風力機會提前進入失速狀態(tài)，且失速時間更長，渦脫落現象更嚴重.圖9中，黑色虛線框內為風力機轉子內部和葉片邊緣的渦流情況.當TSR<1.30時，相較于獨立式風力機，對轉式風力機的葉片邊緣和轉子內部的脫渦現象更嚴重.由于葉片邊緣過多的脫渦會導致葉片升力降低、阻力上升，使得風力機的升阻比減小[15]，所以當TSR<1.30時，對轉式風力機的Cp值更低.

圖10為TSR=2.10時，對轉式與獨立式風力機的轉矩系數(CM)隨θ變化的關系.其中，

圖11為獨立式風力機順時針轉動的俯視圖.基于圖10和11，對圖7中TSR>1.30時，Cps小于Cp1、Cp2的情況進行分析(取TSR=2.10時的風況).

圖10中，在θ=0°～360° 的1個周期內，獨立式與對轉式風力機的CM隨θ的變化過程相似.CM在θ=60°，270°處存在極大值，且在θ=270° 處取得最大值.當θ=60° 時，對轉式風力機的CM大于獨立式風力機的CM，相對誤差為54%，且對轉式的CM曲線整體高于獨立式風力機的CM曲線.這是由于對轉式風力機轉子內脫落渦的擾動更強烈，所以轉矩系數更高[16].圖11中，當θ=60° 時，葉片位于下風區(qū)；當θ=270° 時，葉片位于上風區(qū).因此，當θ=270° 時，風力機的轉矩系數更大.

圖10 風力機的轉矩系數對比Fig.10 Comparison of torque coefficients of different wind turbines

圖12 轉子中平面渦量圖Fig.12 Vorticity in middle plane of rotors

圖12為TSR=2.10、θ=60° 時，對轉式風力機上轉子與獨立式風力機轉子的中平面渦量對比關系.其中，中平面指位于風力機轉子高度1/2處的平面.基于圖12，對圖10中θ=60° 時CM的情況進行分析.圖12(a)中，對轉式風力機的遠端渦流長度較短且脫渦的強度較低，由外流場域前端進入的風能更多地被轉子葉片吸收，轉子后端的損失風能較少.這是由于風力機遠端脫渦的比率和強度對轉子的轉矩系數有一定影響.風力機來流方向的遠端渦流長度較短且脫渦強度較低，說明風力機尾部渦脫落較少，能量散失較低，從而引起風力機轉矩系數增加[17].而對轉式風力機尾部脫渦較少，這可能是由于在風力機對轉過程中，上、下轉子間的脫渦交錯會產生渦流互擾，一定程度上阻礙流體通過風力機轉子向流場域尾部流動，導致更多的風能被風力機轉子吸收，使得尾部脫渦較少.此外，由于風能利用效率與轉矩系數存在正比關系，即

Cp=(DΩ/2v∞)CM

所以，當TSR>1.30時，對轉式風力機的風能利用效率更高.

4 結論

針對海上漂浮獨立式垂直軸風力機的失穩(wěn)問題，提出同軸對轉式垂直軸風力機的設計理念.基于計算流體力學理論，借助RANS SSTk-ω湍流模型對對轉式風力機進行數值模擬，并對比對轉式與獨立式垂直軸風力機的轉矩和風能利用效率.結合渦流理論，進一步分析不同TSR時，對轉式與獨立式風力機的氣動性能差異性的原因.主要結論如下：

(1) 同等外流場條件下，相較于獨立式風力機，對轉式風力機的轉子傳遞到浮式平臺的的總轉矩在0附近波動且波動幅度很小，其浮式平臺的穩(wěn)定性更強.

(2) 當TSR<1.30時，風力機葉片進入失速狀態(tài)的時間長、程度深.在長時間的失速狀態(tài)下，對轉式風力機的葉片邊緣和轉子內部脫渦更嚴重，風能利用效率更低.當TSR>1.30時，相較于獨立式風力機，對轉式風力機的遠端渦流長度較短且脫渦強度較低，風能利用效率更高.

(3) 本文提出的新型結構理念和分析研究方法可以為海上漂浮垂直軸風力機的氣動和穩(wěn)定性能的優(yōu)化提供一定的參考價值.