海上浮式風機氣動性能數值模擬

吳 俊，丁金鴻，何炎平，趙永生，劉亞東

(上海交通大學 船舶海洋與建筑工程學院 海洋工程國家重點實驗室，上海 200240)

海上浮式風機氣動性能數值模擬

吳 俊，丁金鴻，何炎平，趙永生，劉亞東

(上海交通大學 船舶海洋與建筑工程學院 海洋工程國家重點實驗室，上海 200240)

采用計算流體動力學(CFD)方法，基于RANS方程和SST k-ω湍流模型，對OC3-Hywind Spar浮式平臺支撐的NREL 5 MW風機進行氣動性能模擬。對固定式風機的數值模型做網格無關性驗證，同時考慮垂直風切變的影響，并將數值結果與NREL設計數據進行對比以驗證模型的有效性。在FLUENT軟件中,設定嵌入式滑移網格和用戶定義程序(UDF)來模擬風機葉輪隨平臺的周期運動，分別研究浮式平臺的縱蕩、縱搖和首搖運動對風機氣動荷載的影響。數值結果表明平臺的縱蕩和縱搖運動對輸出功率影響較大，且平臺運動幅值越大周期越低，其氣動荷載變化越劇烈。合理控制平臺的運動幅值對提高浮式風機的發電性能和疲勞強度有很大作用。

海上浮式風機；計算流體動力學；氣動性能；六自由度運動

Abstract:Computational fluid dynamic (CFD) method is used to simulate the aerodynamic performance of floating offshore wind turbine (FOWT) experiencing platform motions.Based on the RANS equations and the k-ω shear-stress transport (SST) turbulence model,aerodynamic simulations for NREL 5MW turbine combined with OC3-Hywind spar buoy are studied.Grid independence studies are performed in fixed type offshore wind turbine,considering the effect of vertical wind shear,and the simulation results are compared with NREL published data to verify the numerical model.Furthermore,embedded moving mesh technique and user-defined functions (UDF) are adopted in FLUENT software to model the periodical motions of the rotor,and the effects of typical surge,pitch and yaw motions of the platform on turbine aerodynamic loads are investigated.The simulation results reveal that the surge and pitch motions have a greater impact on the fluctuation of the output power,and more significant variation of the aerodynamic loads would be introduced by larger amplitude and shorter period of the motion.The aerodynamic performance of the turbine and the fatigue strength of the structures would be significantly improved by reasonable controlling of the platform motions.

Keywords:floating offshore wind turbine (FOWT); computational fluid dynamics; aerodynamic performance; six-degrees-of-freedom motions

海上浮式風機在風浪流的聯合作用下產生六自由度運動，使風機周圍的流場發生劇烈變化，其復雜程度遠遠超出了陸上和離岸固定式風機的流場[1]。其中平臺的縱搖和首搖運動將導致葉輪附近出現非軸向或扭曲的流場，且風機葉片也會與尾流相互作用，使葉片的氣動荷載計算變得更加復雜。同時平臺的縱搖運動使葉輪產生有效風剪切，故葉輪表面均勻風速的假設不再適用[2]。縱蕩和縱搖運動引起的附加風速與入流風速疊加將改變葉片的相對入流速度以及葉輪-尾流之間的耦合作用。

目前大型風機氣動性能計算方法主要有：葉素動量理論；渦格法(VLM)和CFD方法。葉素動量理論將葉素切片理論與簡單動量方法相結合，并引入一些修正因子，包括動態失速效應，Glauert修正系數以及Prandtl葉尖損失因子等，該方法廣泛地用到葉片設計和氣動性能計算中。但對于浮式風機，由于高聳于平臺上部的風機隨平臺產生劇烈運動，使風機處于風車和螺旋槳的瞬態變化狀態中[3]，這種狀態無法用傳統的BEM理論來模擬。特別當風機往下風向縱搖時，葉輪與尾部湍流相互作用將產生垂直于葉輪的螺旋形回流，Sebastian[3]將這種瞬態氣動效應定義為渦環狀態 (VRS)。當風機處于渦環狀態時，BEM理論中的動量方程將不再守恒[2]。還有一種方法稱為渦格法，其是基于勢流理論提出的，并在直升機旋翼和螺旋槳分析中得到廣泛應用，Jeon等[4]曾使用該方法研究浮式風機在平臺縱搖運動下的氣動性能響應，但目前該方法在風機氣動性能模擬方面應用較少。而CFD作為風機氣動性能分析的新一代工具，能精確模擬浮式風機由于平臺運動引起的復雜流場。目前，如Chow[5],Bazilevs[6]和劉強等[7]成功運用CFD方法研究了NREL 5 MW大型風機的氣動性能和尾流特性。Tran等[8]用UBEM、CFD和FAST三種方法模擬了浮式平臺的縱搖運動對風機氣動荷載的影響并對三種方法進行對比。同時，斯圖加特大學的研究團隊利用FLOWer軟件和重疊網格技術研究NREL 5 MW浮式風機在指定的縱搖運動下的非穩態氣動響應[9]。任年鑫等[10]基于Fluent軟件中的滑移網格技術，實現了浮式風機葉片在平臺縱蕩下與周圍流場的耦合模擬；同時利用VOF方法和滑移網格，Ren等[11]成功建立了浮式風機與風-浪耦合流場之間的流固耦合模型。

這里利用嵌入式滑移網格和UDF模擬NREL 5 MW風機的旋轉及隨浮式平臺的典型周期性運動。首先對模型做網格無關性驗證，考慮垂直風切變的影響，并將無平臺運動的數值結果與NREL數據對比以驗證數值模型。分別模擬浮式風機在平臺縱蕩、縱搖和首搖運動下的氣動荷載響應，并探討不同運動幅值和周期對風機發電量和推力的影響，得出了平臺運動對浮式風機氣動性能的影響規律。

1 數值計算方法

1.1控制方程

基于雷諾時均的三維黏性不可壓Navier-Stokes方程的連續性方程和動量方程分別為：

1.2湍流模型

采用SST k-ω湍流模型，該模型結合了k-ω模型在近壁區域計算的魯棒性和精確性等優點以及k-ε模型在遠場計算的優點，其輸送方程可寫為：

2 計算模型與網格

2.1浮式風機模型

研究對象選用美國可再生能源實驗室(NREL)所開發的5 MW離岸風機[12]，具體幾何參數如表1所示。沿葉片展向，截面的翼型不同葉厚的DU翼型以及NACA-64618翼型。下部支撐平臺為Offshore Code Comparison Collaboration (OC3)研究計劃所選用的OC3-Hywind Spar式浮式平臺。

表1 NREL 5 MW海上風機總體參數Tab.1 The overall parameters of the NREL 5 MW offshore wind turbine

采用實尺度模型計算風機的氣動性能，為提高計算效率，模型不包括塔架和浮式平臺。整個計算域劃分為外流場域、球形旋轉域和葉片旋轉域，如圖1所示。其中球形旋轉域用于模擬在平臺轉動下整個葉片旋轉域的轉動，而葉片旋轉域用來定義葉輪相對球形旋轉域的自轉。每一時間步迭代前，都需用UDF定義旋轉域的運動速度，由于葉片旋轉域存在嵌套轉動，故需要使用嵌入式滑移網格以定義其運動。兩個區域彼此相互接觸的面都設為交界面，交界面之間存在相對滑移并傳遞流場信息。葉片表面設為無滑移壁面，外部流場分別設置為速度入口和壓力出口，其他四周均為對稱壁面。

圖1 計算模型示意Fig.1 Sketch of the whole computational domain

2.2網格劃分

采用結構化網格，加密旋轉域和葉片附近的網格，且對葉片賦予一定厚度的邊界層。圖2所示為葉片表面與周圍流場網格分布以及葉片展向2/3R處的網格截面圖。

圖2 網格劃分示意Fig.2 Illustration of the mesh

3 模擬結果與分析

采用FLUENT商業軟件，選取SST k-ω湍流模型，采用壓力基求解器，壓力速度的耦合選取PISO算法，壓力項采用PRESTO！格式進行離散。擴散項為中心差分格式，對流項、湍流動能方程以及湍流耗散率均采用二階迎風格式。

3.1網格無關性驗證

為確保數值模擬的精確性，需進行網格無關性驗證。其中計算域寬度與高度均取10倍葉輪半徑，上風向長度為6倍半徑，為滿足壓力出口條件，下風向取20倍半徑。對三種不同網格方案模型，在額定風速下進行氣動荷載模擬，并將結果與設計值對比，如表2所示。由表可知方案1的網格數量最少，但模擬結果偏差較大；而方案2的網格數量較適中，且功率和推力的偏差與方案3較接近，為了提高計算效率，取方案2，最小單元尺寸為2 cm，網格總量為832萬。

表2 不同網格精度的對比Tab.2 The comparison of different grid sizes

3.2考慮垂直風切變的影響

5 MW風機尺度巨大，輪轂高度距海平面90 m，沿高度方向的風速存在較大梯度，故需考慮垂直風切變的影響。根據API-RP-2A-WSD規范，采用對數風剖面模擬海上風速沿高度的變化：

其中，U0為海平面10 m高度處的1小時平均風速(m/s)。

將風速梯度編制成UDF，初始化流場速度入口，以模擬在額定風速下風機氣動荷載的變化。由于越往高處，風速越大，轉至高處的葉片受到的推力與扭矩比處于低處的葉片要大，而我們所關心的風機總體氣動荷載是對三枚葉片求和得到的。如圖3所示，雖然葉輪輸出功率和推力存在一定波動，以葉輪轉1/3圈為周期類似簡諧曲線變化，但總體波動很小，在0.2%～0.4%之間，其中平均功率相對均勻風的工況降低了1.2%。總體來說，垂直風切變對風機總體氣動性能影響較小，為了簡化，以下工況均取速度入口為均勻風。

圖3 垂直風切變對風機總體氣動性能的影響Fig.3 The effect of vertical wind shear on overall aerodynamic performance of wind turbine

圖4 風機后方不同位置截面的流場云圖Fig.4 Contour of velocity from different sections behind the rotor

3.3數值模型驗證

根據NREL 5 MW風機轉速隨來流風速的變化情況，分別模擬三個低風速工況、一個額定風速以及兩個高風速工況。首先采用旋轉參考坐標系法(MRF)進行穩態模擬，當收斂精度達到10-4的且氣動荷載變化較小后，即認為流場穩定。而瞬態計算是在穩態計算的基礎上采用滑移網格法進行模擬，當葉輪旋轉4圈后，推力和扭矩基本穩定。圖4顯示了風機后方穩定流場截面的速度云圖，當來流剛經過風機后，存在明顯的渦流，隨著流動距離不斷加大，漩渦逐漸消失。將以上的數值結果與NREL的設計數據進行對比，由圖5可知，CFD的模擬結果略小于NREL設計數據。當風速較小時，CFD結果與設計值相差很小，當風速高于額定風速時，相差越大，其中當風速達到15 m/s時，水平推力相差9.4%，而輸出功率相差9.1%。這主要是由于當風速較大時，葉片中部翼型截面處的流場逐漸進入動態失速狀態，故在有限的網格數量和所選湍流模型下很難精確模擬出該復雜流場[10]。鑒于模擬結果的誤差在可接受范圍內，從一定程度上也驗證了該數值模型模擬海上風機氣動荷載的有效性。

圖5 風機推力和輸出功率的對比Fig.5 Comparisons of fixed wind turbine thrust and rotor power

3.4平臺縱蕩的影響

當浮式平臺發生縱蕩時，風機與平臺以相同的速度運動。為了簡單起見，這里將縱蕩運動簡化為一定周期和幅值的簡諧運動，如式(7)所示。對上式求導即為縱蕩速度，當運動開始時，該速度從零開始增加，故流場不會因為平臺的運動而發生突變。在開始加載縱蕩運動前，先計算600步左右的穩態模擬和4個葉片旋轉周期的瞬態模擬，以確保流場已充分發展并消除瞬態波動。模擬工況選取額定風速為11.4 m/s，葉片旋轉速度為12.0 rpm，浮式平臺運動周期分別為12 s和20 s，運動幅值為2 m和4 m，葉輪推力及輸出功率的變化如圖6和圖7所示。

圖6 不同運動幅值下浮式風機氣動性能比較Fig.6 Comparison of blade aerodynamic performance under different surge amplitudes

圖7 不同運動周期下浮式風機氣動性能比較Fig.7 Comparison of blade aerodynamic performance under different surge periods

由圖6和圖7可知，風機氣動荷載的響應曲線也近似為簡諧曲線，與平臺的運動規律基本一致。由于平臺的縱蕩運動產生的附加速度，使風機與入流風之間的相對速度改變，從而導致葉片各截面翼型的攻角、升力系數和阻力系數發生周期性變化，最終反映到風機推力和功率的周期性變化。在縱蕩運動下，推力和輸出功率的平均值與固定式風機相比降低1%左右，主要原因是氣動荷載的簡諧曲線關于中心軸不是對稱的，即風機往下風向運動所減小的氣動荷載值要略高于往上風向運動所增加的值。對于周期相同而幅值為2 m和4 m的兩種縱蕩工況，風機推力變化幅度分別為10%和18%，而功率變化為22%和68%。該結果表明在縱蕩運動下，平臺運動幅值越大，周期越小，氣動荷載響應的變化越劇烈，且輸出功率變化較為敏感，故合理控制平臺運動幅值對保證風力發電的穩定以及風機結構包括塔架、葉片和減速齒輪的疲勞壽命都有重要意義。

下一步分析縱蕩運動對每個葉片的影響。考慮到葉輪的對稱性，只取一個葉片作為研究對象，平臺運動選取周期為12 s而幅值為4 m的工況。圖8給出了從縱蕩運動開始到氣動荷載達到最大值期間，推力和轉矩沿葉片展向的分布情況。當t=20 s時，此時縱蕩速度為零，故氣動性能的分布曲線處于中間位置；而當t=23 s時，即風機的下風向速度達到最大，此時轉矩和推力的分布都處于最低位置；當t=26 s時，風機的縱蕩速度又回到零，故與t=20 s的曲線幾乎重合；而當t=29 s時，風機的上風向速度達到最大，故葉片氣動性能分布曲線也處于最高位置。由圖8可知，葉尖和靠近輪轂部分的氣動荷載變化較小，而變化最劇烈的范圍集中在0.5R～0.9R之間，故此處需著重考慮葉片的結構和疲勞強度。

圖8 在縱蕩運動下氣動荷載沿葉片展向的分布Fig.8 The spanwise distributions of aerodynamic loads along the blade experiencing surge motion

3.5平臺縱搖的影響

將平臺的縱搖角位移簡化為簡諧運動，根據OC3-Hywind spar結構，取輪轂以下90 m作為縱搖中心。對于錨泊式浮式風機，其典型縱搖運動的頻率接近入射波浪的頻率，故這里選取周期為10 s且角位移幅值分別為2°和4°的兩種縱搖工況進行分析。圖9顯示了浮式風機縱搖過程的尾流變化情況，當風機縱搖處于平衡位置時，葉片與后方流場相互干擾較小，故尾流相對穩定。而當風機縱搖至上風向和下風向時，由于葉片的一部分受到尾流的影響，同時由于葉片的存在對流場產生一定擾動，故葉輪與尾流之間的耦合作用使尾流變化更劇烈。圖10給出了在兩種工況下風機氣動荷載的變化。由圖10可知，縱搖角度越大所引起氣動荷載的波動也越劇烈。與無平臺運動的情況相比，在幅值為2°和4°的縱搖運動下推力的平均值分別降低2.7%和5.2%，而兩者的功率平均值基本一致，均降低4.2%。其原因可能是縱搖運動使風機在入流方向的投影面減小，縱搖幅值越大投影面越小，從而受到的水平推力也就越低。對于縱搖角度為2°和4°的運動，水平推力的變化幅度分別為17.5%和35.3%，而輸出功率變化幅度分別為23.8%和65.9%。與縱蕩運動類似，縱搖運動產生的氣動荷載的劇烈波動也會對浮式風機的機艙控制系統和結構疲勞造成一定影響。

當風機發生縱搖時，風輪盤面存在的附加速度差會引起有效風剪切[3]，如圖11所示。為了研究這種效應對風機氣動荷載的影響，分別計算每個葉片上的荷載響應。圖12給出了三個葉片在縱搖運動下的水平推力的響應曲線，其中葉片#1的初始位置在Z軸正方向，葉片#2和#3與Z軸分別成±120°。由于葉輪的旋轉作用，三個葉片的氣動荷載曲線存在相位差，其中葉片#1比#3超前2/3π相位，比#2落后2/3π相位。由于葉片#1的截面到縱搖中心的垂直距離最長，相應其附加速度的變化也是最大的，故如圖12所示，作用在葉片#1上的推力的變化幅值比其他葉片要大。然而這種差別低于5%，故由于縱搖引起的風剪切效應對風機葉片氣動荷載的影響并不大。

圖9 浮式風機縱搖至不同位置的尾流速度圖Fig.9 Velocity contours of wake flow during different pitch positions of floating wind turbine

圖10 不同縱搖幅值下風機的氣動性能對比Fig.10 Comparison of blade aerodynamic performance under different pitch amplitudes

圖11 縱搖運動所引起的風輪盤面的有效風剪切Fig.11 Effective wind shear across the rotor due to pitching motion

圖12 作用在每個葉片上的水平推力Fig.12 The aerodynamic thrust force for each blade

3.6平臺首搖的影響

與縱搖類似，平臺首搖運動也會在葉輪盤面處產生風剪切，并改變相對入流速度，同時也造成葉片槳距角的改變，即葉片不同截面翼型處的攻角也相應發生變化。在附加速度的影響以及葉片槳距角的改變的共同作用，分析風機在首搖運動下的氣動荷載是非常復雜的。文中研究了風機在周期為10 s、角幅值為2°的首搖運動下的氣動性能，圖13給出了風機在該運動下的總體氣動荷載變化以及葉片#1的推力和扭矩的變化。由圖可知，除了初始階段的數值波動外，風機的總體性能相對較穩定，而單獨每個葉片的氣動荷載存在較大的波動。為了研究葉片翼型截面在首搖作用下的氣動性能，借助BEM理論進行定性分析，如圖14所示。以葉片#1為例，當首搖運動開始時(t=10 s)，葉片繞Z軸的旋轉力臂為零，即槳距角和附加速度都從零開始增加，故此時葉片的每個截面的攻角都沒有改變。顯然，附加速度與首搖角速度和轉動力臂成正比，附加速度增加會導致翼型攻角和氣動荷載的增大，而增大逆時針方向的槳距角會降低翼型攻角和氣動荷載。當t處于10～11 s之間時，雖然由于槳距角的增大會減小攻角，但附加速度的增加占主導地位，故此階段推力和轉矩都增大。過了11 s后，由于首搖角速度減小使附加速度的增幅放慢，而此時槳距角的增大占主導地位，故葉片氣動荷載下降。到了1/4T(t=12.5 s)時，首搖角速度減為零，即附加速度為零，同時逆時針方向的槳距角也達到了最大值，故此時葉片扭矩降到最低值，而由于升力在推力方向上的分量增加所帶來的補償，故水平推力并沒有降到最低值而是略低于正常工況。

圖13 在首搖運動下風機的氣動性能Fig.13 The aerodynamic performance of the rotor experiencing yaw motion

圖14 利用BEM理論分析平臺首搖的影響Fig.14 Analysis of the effect of yaw motion using BEM theory

4 結 語

1)在模型的網格無關性驗證基礎上，模擬垂直風切變的工況，結果表明風機總體氣動性能變化很小可以忽略。同時不同風速工況的模擬結果表明CFD的模擬結果略小于NREL設計數據。在低風速下，兩者相差很小，而在高風速下，相差略大，然而誤差在可接受范圍內，一定程度上驗證了該數值模型的有效性。

2)在平臺縱蕩運動下，浮式風機的推力和輸出功率的平均值相比固定式減少了1%。平臺運動幅值越大，周期越小，氣動荷載變化越劇烈，同時功率的變化較推力值更為敏感。當平臺縱蕩時，葉片氣動荷載變化最劇烈的部位集中在0.5R～0.9R之間。

3)平臺縱搖角度越大，其氣動荷載波動也越劇烈，且輸出功率的變化比推力更敏感，同時由于縱搖引起的風剪切效應可忽略不計。在縱搖運動下，浮式風機的推力和輸出功率的均值較固定式風機要低，且運動幅值越大推力的均值越低。

4)在平臺首搖運動下，風機的總體性能相對穩定，而每個葉片的氣動荷載波動較大。由于首搖角速度產生的附加速度以及葉片槳距角的改變，兩者的共同作用使風機的氣動性能預報變得更加復雜。

[1] SEBASTIAN T,LACKNER M.Analysis of the induction and wake evolution of an offshore floating wind turbine[J].Energies,2012,5(12):968-1000.

[2] SEBASTIAN T,LACKNER M A.Characterization of the unsteady aerodynamics of offshore floating wind turbines[J].Wind Energy,2013,16(3):339-352.

[3] SEBASTIAN T,LACKNER M.A comparison of First-order aerodynamic analysis methods for floating wind turbines[C]//Proceedings of the 48th AIAA Aerospace Sciences Meeting.2010.

[4] JEON M,LEE S,LEE S.Unsteady aerodynamics of offshore floating wind turbines in platform pitching motion using vortex lattice method[J].Renewable Energy,2014,65:207-212.

[5] CHOW R,VAN DAM C P.Verification of computational simulations of the NREL 5 MW rotor with a focus on inboard flow separation[J].Wind Energy,2012,15(8):967-981.

[6] BAZILEVS Y,HSU M C,AKKERMAN I,et al.3D simulation of wind turbine rotors at full scale.Part I:Geometry modeling and aerodynamics[J].International Journal for Numerical Methods in Fluids,2011,65(1-3):207-235.

[7] 劉強,楊科,黃宸武,等.5 MW 大型風力機氣動特性計算及分析[J].工程熱物理學報,2012,33(7):1155-1159.(LIU Qiang,YANG Ke,HUANG Chengwu,et al.Simulation and analysis of the aerodynamic characteristics of a 5 MW wind turbine[J].Journal of Engineering Thermophysics,2012,33(7):1155-1159.(in Chinese) )

[8] TRAN T,KIM D,SONG J.Computational fluid dynamic analysis of a floating offshore wind turbine experiencing platform pitching motion[J].Energies,2014,7(8):5011-5026.

[9] MATHA D,SCHLIPF M,CORDLE A,et al.Challenges in simulation of aerodynamics,hydrodynamics,and mooring-line dynamics of floating offshore wind turbines[R].National Renewable Energy Laboratory,US Department of Energy,Office of Energy Efficiency and Renewable Energy,2011.

[10] 任年鑫,李玉剛,歐進萍.浮式海上風力機葉片氣動性能的流固耦合分析[J].計算力學學報,2014,31(1):91-95.(REN Nianxin,LI Yugang,OU Jinping.The fluid-structure interaction analysis of aerodynamic performance of floating offshore wind turbine blade[J].Chinese Journal of Computational Mechanics,2014,31(1):91-95.(in Chinese) )

[11] REN N X,LI Y G,OU J P.Coupled wind-wave time domain analysis of floating offshore wind turbine based on computational fluid dynamics method[J].Journal of Renewable and Sustainable Energy,2014,6(2):23106.

[12] JONKMAN J,BUTTERFIELD S,MUSIAL W,et al.Definition of a 5 MW reference wind turbine for offshore system development[R].National Renewable Energy Laboratory Golden,CO,2009.

Numerical analysis of aerodynamic performance of floating offshore wind turbine

WU Jun,DING Jinhong,HE Yanping,ZHAO Yongsheng,LIU Yadong

(State Key Laboratory of Ocean Engineering,School of Naval Architecture,Ocean and Civil Engineering,Shanghai Jiao Tong University,Shanghai 200240,China)

P753; TM614

A

10.16483/j.issn.1005-9865.2016.03.005

1005-9865(2016)03-0038-09

2015-04-30

國家重點基礎研究發展計劃(973 計劃)項目(2014CB046200)；高等學校博士學科點專項科研基金新教師類資助課題(20120073120014)

吳 俊(1990-)，男，浙江湖州人，碩士，主要從事海上風機的氣動荷載模擬。E-mail:wujun2009110@sjtu.edu.cn

丁金鴻。E-mail：ahaha@sjtu.edu.cn