復雜大氣入流下海上風機力學特性研究1)

徐順 趙偉文 萬德成

(上海交通大學船舶海洋與建筑工程學院,船海計算水動力學研究中心(CMHL),上海 200240)

引言

近年來,傳統化石能源的不可再生性以及對環境造成的污染,使其難以滿足當今社會對能源的巨大需求.風能由于其無污染、可再生及儲量大等特點逐漸成為最具發展潛力的新能源之一[1].2021 年全球風能報告指出[2],2020 年新增裝機容量高達93 GW,同比增長53%.風能技術的不斷進步,促使風機尺寸逐漸向大型化發展,并導致大氣邊界層(atmospheric boundary layer,ABL)復雜來流對風機的運行性能、疲勞載荷等方面產生惡劣影響.因此,研究大氣邊界層內風機的運行性能和力學特性具有重要意義.

目前,生成滿足大氣邊界層復雜來流統計特性的方法主要有兩種[3]:人工合成和域前模擬.人工合成的基本思想是在入口平面處基于嚴格的數學推導給出入流風速的計算公式.寧旭[4]采用諧波合成法生成了符合高雷諾數湍流能譜分布的隨機風速序列,并將其應用于風機數值模擬.Huang等[5]提出了“DSRFG”方法,該方法生成的脈動風速場嚴格滿足連續性條件且能夠并行計算生成不同空間位置的風速時程.Castro 和 Paz[6]在“DSRFG”方法的基礎上,通過引入時間尺度參數考慮了脈動風速場的時間相關性.但是,人工合成法給定的風速不能嚴格滿足NS 方程,需沿下游發展一段距離才能獲得符合統計特性的大氣湍流風場.Xie等[7]對傳統2D 渦方法進行了參數化研究,提出了改進的2D 渦方法,將符合湍流統計特性的下游自適應距離由原來的5h減小到3h(h為計算域的高度).

域前模擬的基本做法是預先生成一個符合大氣湍流統計特性的風場,并將其作為計算區域的入口邊界條件.在域前模擬中,生成大氣湍流風場的方法主要有兩種:利用湍流激發裝置以及依靠表面粗糙度.在實驗中經常利用湍流激發裝置(如粗糙元、網格板等)來促使大氣湍流風場的生成[8-10],后來有部分學者將其應用于數值模擬中.Phuc等[11]利用尖劈和粗糙元構建了數值風洞,其數值模擬出的湍流風場統計特性與實驗值符合良好.周桐等[12]利用尖劈和粗糙元等湍流激發裝置數值模擬了大氣邊界層入口湍流,與“CDRFG”方法對比,該方法生成的流場結構更加合理.但是,該方法模擬出的湍流風場特性高度依賴于湍流激發裝置附近的網格分辨率,且激發裝置與流場相互作用機理機制尚不明確.

另一種域前模擬方法是直接依靠計算域底部的粗糙度,無需進行復雜的網格劃分,直接長時間大范圍數值模擬出大氣邊界層湍流風場,如美國國家可再生能源實驗室研發的風機氣動性能數值計算軟件“SOWFA”便基于此方法生成復雜大氣入流[13].Lee等[14]研究了不同表面粗糙度和穩定性的大氣入流及風機尾流對下游風機疲勞載荷的影響,結果表明表面粗糙度和上游風機尾流對風機疲勞載荷影響顯著.Ning 和Wan[15]基于LES 結合致動線模型,研究了不同大氣穩定性下尾流彎曲效應及其對下游風機氣動性能的影響.白鶴鳴等[16]基于“SOWFA”軟件,對對流大氣邊界層下錯列排布三風機進行數值模擬,發現垂向錯列可以有效提升后排風機的功率,且不會惡化其氣動性能.李德順等[17]利用LES 結合致動線模型的方法對一臺處于中性大氣入流下的外場試驗機組進行了葉根載荷分析,發現葉根揮舞載荷對大氣中的湍流結構響應明顯.Johlas等[18-19]通過結合“SOWFA”和“OpenFAST”[20],對大氣復雜入流下浮式風機的尾流場進行了分析,并對比了不同浮式支撐平臺對風機尾流場和平臺六自由度運動的影響.

以上研究中,大都采用域前模擬中的依靠地表粗糙度的方式來生成風機的大氣復雜入流,但其中卻少有文獻重點關注風機的力學特性響應.因此,本文采用域前模擬的方法,通過定義典型海面粗糙度,直接長時間大范圍數值模擬真實環境下中性大氣邊界層復雜入流.并結合LES 和致動線模型,對中性大氣復雜入流下的海上固定式風機進行數值模擬,并以均勻入流計算結果作為參考基準,定量分析大氣復雜入流對風機氣動性能和葉根載荷的影響.

1 數值方法

1.1 致動線模型

風機葉片利用致動線模型進行模擬,該方法無需求解葉片表面邊界層,能大幅降低計算成本并獲得較為精確的結果.其最早由S?rensen 和Shen[21]提出,基本思想是利用致動線上的致動點代表沿著徑向離散化的風機葉片,然后基于葉素理論計算致動點的體積力,并將其反作用于流場以體現風機葉片對流場的作用.致動點上的體積力可采用下式計算

式中,L和D分別代表葉片半徑r處葉素的升力和阻力,ρ 為空氣密度,Urel為相對入流風速,c為葉片半徑r處二維翼型的弦長,dr為葉素的寬度,CL和CD分別為升力系數和阻力系數,可根據攻角確定,eL和eD為升力和阻力方向的單位矢量.

圖1 顯示了葉片半徑處二維翼型的速度矢量,其中相對入流風速可由下式確定

圖1 二維翼型的速度矢量Fig.1 Velocity vectors of two-dimensional airfoil

式中,Uz和Uθ分別為入流風速的軸向速度和切向速度,Ω 為轉子角速度.

另外,若直接將致動點上的體積力作用于流場,會造成數值奇異性,在此采用高斯核函數進行體積力光順投影,光順后的體積力表達式如下

式中,N為葉片上致動點的數目,(xi,yi,zi) 表示第i個致動點的位置,di為致動點距離投影點的距離,ε 代表投影的寬度,本文取 ε ≈2Δx[22],Δx為葉片附近的網格尺寸.

1.2 控制方程

本文采用大渦模擬方法,對中性復雜大氣入流下的海上固定式風機進行數值模擬,空間過濾后的流場的控制方程如下

式中,動量方程右端第一項為修正壓力梯度項,第二項為背景壓力梯度項,用以克服計算域底部表面粗糙度并驅動大氣湍流的生成,第三項為考慮地球自轉的科氏力項,第四項代表溫度引起的浮力項,第五項為流體應力張量項,采用“Smagorinsky”亞格子模型封閉,第六項為風機葉片體積力源項,由致動線模型計算得出,在域前模擬階段忽略此項.另外,還需要求解位溫輸運方程以獲得位溫場

式中,熱通量qj由下式計算

式中,Prt為普朗特數,在中性和對流大氣邊界層中取 1/3,關于控制方程更詳細的說明可查閱參考文獻[23].對于均勻入流下的風機,其流場的控制方程未考慮科氏力項和溫度項的影響,也不需要背景壓力梯度驅動流體流動,在此不再重復贅述.

1.3 局部坐標系

在本文的結果分析中,需要對風機轉子氣動力和力矩以及葉片根部的力和力矩進行結果分析,兩者所采用的坐標系不同,在此進行簡要說明.

風機轉子氣動力和力矩采用的是輪轂局部坐標系,如圖2(a)所示,其坐標原點為轉子軸與旋轉平面的交點,x軸為沿著輪轂中心線指向下游方向,z軸垂直于輪轂中心線并與葉片具有相同的方位角,y軸根據右手螺旋定則確定.

圖2 局部坐標系示意圖Fig.2 Local coordinate systems

葉片根部的力和力矩基于葉片局部坐標系,如圖2(b)所示,其中i代表葉片編號.y軸平行于弦線并指向葉片后緣,z軸為沿著槳距軸指向葉尖,x軸根據右手螺旋定則確定.另外,兩個局部坐標系均會隨著葉片旋轉.

2 算例設置

2.1 風力機參數

本文所采用的研究對象為美國國家可再生能源實驗室(National Renewable Energy Laboratory,NREL)開發的5 MW 風機,其是一種傳統迎風向可變槳距和轉速的3 葉片風力機.該風機的額定風速為 1 1.4 m/s,對應的額定功率和額定轉速分別為5 MW 和12.1 r/min,該風機的主要參數如下表1 所示,關于該風機更詳細的數據可參考文獻[24].

表1 NREL 5 MW 風機主要參數Table 1 Main parameters of NREL 5 MW wind turbine

2.2 大氣入流算例設置

考慮真實復雜大氣入流的海上固定式風機數值模擬需要分兩階段進行:域前模擬和主模擬.在域前模擬階段,需要長時間大范圍的依靠底面粗糙度來生成真實的大氣湍流風場.如圖3 所示,計算域的四周側壁均采用周期性邊界條件,以達到通過減小計算域尺寸降低計算量的目的,計算域頂部設置為滑移邊界條件,底面采用“Moeng”壁面模型[25],表面粗糙度設置為0.001,對應典型低粗糙度的海面情況.計算域的長度、寬度和高度分別設置為3000 m×3000 m×1000 m,三個方向的網格分辨率為10 m×10 m×10 m,對應的網格數量為900 萬.輪轂高度處的風速設置為風機的額定風速 1 1.4 m/s,入流風向設置為240°(西南方向入流)以避免在周期性邊界條件中同一高度處的渦結構幾乎同時到達下游邊界時又同時傳入上游邊界時“卡死”現象的出現.

圖3 域前模擬計算域及邊界條件(單位:km)Fig.3 Calculation region and boundary conditions of precursor simulation (unit:km)

對于初始位溫的設置,從底部至700 m 高度設置300 K 均勻分布,700 m 至800 m 高度設置逆溫反轉層,其中的溫度線性增長到308 K,800 m 至1000 m 高度的溫度以 0 .003 K/m 的速率線性增長.為了使生成的大氣湍流達到準平衡態,域前模擬階段的計算時間為19 000 s,計算步長為0.4 s,保留最后1000 s 的計算數據作為主模擬階段的輸入.

主模擬階段的計算域和背景網格的分辨率與域前模擬階段相同,但北面和東面側壁邊界修改為零梯度邊界條件以防止下游風機尾流循環進入上游邊界,如圖4 所示.另外,為了更加精確的捕捉風機尾流區域的流場細節,對風機附近的區域進行了二級網格加密,加密后風機附近的網格尺寸為2.5 m×2.5 m×2.5 m,總網格數量約為1400 萬.時間步長的設置需要滿足CFL 收斂性條件

圖4 主模擬計算域布置 (單位:m)Fig.4 Successor simulation calculation region (unit:m)

式中,U=11.4 m/s,為輪轂高度入流風速;R=63 m,為轉子半徑;Ω=12.1 r/min,為入流風速下的轉子轉速;Δx=2.5 m,為風機葉片附近的網格分辨率.代入上式,可得 Δt＜0.031 s,因此,本文的時間步長設置為0.02 s.計算時長為1000 s,對應于域前模擬階段的流場數據保存時間,并利用最后350 s 的計算數據進行結果分析.本課題組在風機數值模擬方向具有扎實的研究基礎,關于風機附近網格劃分策略的詳細細節可參考文獻[26-28].

2.3 均勻入流算例設置

如圖5 所示,均勻入流算例計算域的長度、寬度和高度分別設置為 1 0D×3D×3D,D為轉子直徑126 m.風機距離上游邊界為 2D,同樣,與主模擬做法相一致,背景網格的尺寸為 1 0 m×10 m×10 m,并對風機葉片附近區域采取二級網格加密策略,加密后的網格尺寸為 2 .5 m×2.5 m×2.5 m,網格數量約為320 萬.另外,上游邊界為速度入口邊界條件,對應于11.4 m/s 的均勻入流,入流風向為270°(西側方向入流),下游邊界為自由流出邊界條件,四周壁面為周期性邊界條件.時間步長為0.02 s,計算時間為400 s,并利用最后350 s 的計算數據進行結果分析.

圖5 均勻入流計算域布置Fig.5 Calculation region of uniform inflow

3 結果分析

3.1 氣動功率

圖6 顯示了計算結果最后350 s 中性大氣入流和均勻入流下海上固定式風機的功率時歷曲線.可以看出,均勻入流下的轉子氣動功率輸出較為穩定,但大氣復雜入流下轉子功率輸出變化十分復雜,轉子功率在t=25 s 附近開始存在一個較低值,且持續時間超過50 s,這是由大氣復雜入流中的大尺度低速氣流導致.另外,隨著時間的增加,氣動功率值迅速增加至6 MW 附近,在t=90 s 附近又存在較大的下降梯度,隨后又逐漸增加,并在均勻入流風機的氣動功率值附近振蕩.

圖6 功率對比Fig.6 Comparison of rotor power

在數據統計特性方面,分別統計了氣動功率的最大值、最小值、平均值、均方根和標準差,如表2所示.可以看出,兩種工況的氣動功率均方根值差異不大,分別為5.22 MW 和5.30 MW,因為氣動功率的均方根值主要由平均風速決定,兩種計算工況的平均風速均在11.4 m/s 附近.但是,大氣入流中的復雜湍流會引起氣動功率振蕩幅值的極大增加,進而引起標準差的增加,其值約為均勻入流氣動功率標準差的4 倍.

表2 氣動功率統計值Table 2 Aerodynamic power statistics

3.2 軸向推力和偏航力矩

軸向推力和偏航力矩屬于風機轉子的氣動力和力矩,其局部坐標系為輪轂坐標系,其中軸向推力沿輪轂坐標系x軸,偏航力矩繞z軸.圖7(a) 和圖7(b)分別顯示了大氣入流和均勻入流兩種工況下風機轉子的軸向推力和偏航力矩對比曲線.與氣動功率輸出類似,均勻入流下風機的軸向推力穩定在600 kN 附近,而復雜大氣入流下的軸向推力輸出無論是變化幅度還是變化梯度均較大,這同樣歸因于復雜大氣入流下的大尺度湍流來流.

圖7 軸向推力和偏航力矩對比Fig.7 Comparisons of rotor thrust and yaw moment

當風機轉子的其中兩個葉片位于塔架同一側時會導致盤面受力的橫向不對稱性,并由此產生偏航力矩.因此,即便在均勻入流下,風機偏航力矩也會隨葉片旋轉作周期性變化,變化周期對應于葉片旋轉周期.然而,復雜大氣入流下的流場擾動會使得風機的偏航力矩隨時間呈現明顯區別于均勻入流工況的復雜變化,明顯加劇偏航力矩的周期性變化幅度,并顯著增加其均方根值,如圖7(b)所示.因此,在風機偏航優化控制中需重點關注復雜大氣入流帶來的不利影響[29].

表3 顯示了兩種入流工況下軸向推力和偏航力矩的統計值.雖然軸向推力的均方根值差異不大,但復雜大氣入流下的軸向推力標準差與均勻入流相比增大了約53 倍.另外,針對偏航力矩,復雜大氣入流的流場擾動極大的增加了其幅值,偏航力矩的最大值和最小值分別增加到均勻入流的10 倍和8 倍,并導致均方根和標準差的比值為均勻入流的4.3 倍,相關文獻[4]也表明復雜入流下各類結構載荷的標準差值較均勻入流時均有大幅增加.

表3 軸向推力和偏航力矩統計值Table 3 Statistics of rotor thrust and yaw moment

3.3 擺振剪力和彎矩

擺振剪力和彎矩的參考坐標系為風機葉片局部坐標系,其中擺振剪力沿x軸,擺振彎矩繞y軸.圖8分別顯示了擺振剪力和彎矩的對比曲線由圖8(a)可以看出,均勻入流下葉根部位的擺振剪力隨葉片旋轉作周期性變化,變化周期對應于葉片旋轉周期.另外,大氣入流下的擺振剪力也近似作周期性變化,但是,由于復雜大氣入流下的高強度湍流結構及速度沿垂向分布的不均勻性,導致擺振剪力的變化幅度較為劇烈且隨機.擺振彎矩主要由擺振剪力引起,兩者的變化趨勢基本一致,在此不再重復贅述.

圖8 擺振剪力和彎矩對比Fig.8 Comparisons of flapwise shear force and bending moment

表4 給出了擺振剪力和彎矩的統計值,可以看出,復雜大氣入流和均勻入流工況下擺振剪力和彎矩的均方根值差異不大.但是,由于復雜大氣入流的流向風速在垂向上存在速度梯度,當風機葉片旋轉至輪轂平面上方時,擺振剪力和彎矩值較均勻入流增大;反之,較均勻入流減小,這導致了擺振剪力和彎矩標準差值的顯著增加.

表4 擺振剪力和彎矩統計值Table 4 Statistics of flapwise shear force and bending moment

3.4 速度云圖

圖9 顯示了復雜大氣入流下某時刻的流場速度云圖.由圖9(a)可以看出,輪轂高度處的來流速度同時存在大尺度的高速和低速的氣流團結構,并且在風輪旋轉平面處存在明顯的不對稱性,使得偏航力矩振蕩幅值急劇增加.由于風機吸收了入流風的能量,導致風機盤面后方的尾流區域存在明顯的速度損失.本文并未對輪轂進行建模,因此在輪轂中心處存在一條狹長的高速氣流,該氣流隨著向下游傳播迅速和低速風機尾流相互摻混.另外,尾流膨脹和蜿蜒效應也清晰可見.

圖9 瞬時速度尾流場Fig.9 Instantaneous velocity wake field

圖9(b)展示的是240°入流風垂直面的瞬時尾流速度云圖,可以看出,由于海面的阻滯作用,在靠近計算域底部高度處存在著大尺度低速氣流結構.另外,在輪轂高度附近存在大尺度低速羽流結構,其與計算域底部的大尺度低速氣流團相互作用會顯著增加葉片擺振剪力和彎矩的振蕩幅值.在風輪旋轉平面下游,可以明顯看出風機尾流與外部大氣之間的摻混作用.但由于海面低粗糙度的作用,使得中性復雜大氣入流的環境湍流強度較高地表粗糙度時低,導致風機尾流速度隨下游流向距離恢復較慢.

3.5 尾渦結構

第三代渦識別方法能夠合理回答渦的六大要素問題,并定量化表示渦的特性.因此,本文采用第三代“Liutex”渦識別方法[30],對中性復雜大氣入流下海上固定式風機的尾渦結構進行可視化識別.根據圖10 可以看出,由于中性大氣入流的復雜湍流場,導致在風機上游就已形成多尺度湍流渦結構.風機葉片的旋轉會誘導出清晰可見的葉尖渦,該葉尖渦與復雜大氣入流渦結構的相互作用引起葉尖渦的迅速破碎,并與小尺度環境渦結構相互摻混向下游傳播.另外,風機尾渦的膨脹效應也清晰可見.

圖10 尾渦結構(| Liutex|=0.18)Fig.10 Vortex structures (| Liutex|=0.18)

4 結論

本文采用基于大渦模擬的域前模擬方法,生成了真實環境下中性大氣邊界層復雜入流,并利用致動線模型對風機葉片進行數值模擬,定量研究了中性復雜大氣入流下海上固定式風機的力學特性響應,重點分析了風機的氣動性能以及轉子和葉片根部的力和力矩,通過與均勻入流工況進行對比,得出了以下結論.

(1)氣動功率的均方根值由平均入流風速決定,但是由于中性大氣來流的復雜性,導致氣動功率輸出振蕩幅值加劇.另外,大氣入流中的大尺度低速氣流團使得氣動功率輸出值在較長一段時間內過低.

(2)中性大氣復雜入流使得風機軸向推力的標準差相較與均勻入流急劇增加,風機盤面前的復雜來流擾動引起偏航力矩最值、均方根值以及標準差值的較大增加.

(3) 中性復雜大氣入流流向速度沿垂向的梯度、流場的復雜性以及輪轂高度處大尺度低速羽流結構的共同作用,加劇了葉片根部擺振剪力和彎矩的響應幅值和標準差值.

隨著風機逐漸向大型化發展,復雜大氣入流會進一步惡化風機葉片的入流條件,增加風機葉片的力學響應幅值,并引起標準差值的極大增加.因此,復雜大氣入流下風機葉片的力學響應特性需得到更加廣泛關注.