基于非定常面元法的海上浮式風機氣動性能研究

劉利琴，陳迪郁，沈文君，李 昊，鄧萬如

(1. 天津大學 水利工程仿真與安全國家重點實驗室，天津 300072； 2. 交通運輸部天津水運工程科學研究所 工程泥沙交通行業重點實驗室，天津 300456)

葉片是風力機的關鍵部件，其表面形狀對整個風機的空氣動力學性能和風力機的風能利用率都會產生重要的影響，因此葉片的氣動性能一直都是風力機的研究重點。為評估風力發電機效率和結構安全，需要正確預估不同來流條件下葉片上的風載荷。特別是在研究海上浮式風機時，自身的繞流情況會變得更為復雜。為了提高風能的利用效率，必須要準確快速地解決風機面臨的空氣動力學問題。

目前計算風機氣動載荷的常用方法主要有葉素動量理論(BEM)、三維勢流理論、計算流體力學(CFD)方法。其中，葉素動量理論基于工程經驗提出，計算速度快，但不能很好地考慮風力機流場尾流特性；計算流體力學(CFD)方法可以準確地考慮流場黏性及尾流，但計算速度非常慢，難以用于大規模工程計算；基于三維勢流理論的面元法計算風力機氣動載荷逐漸受到人們的重視。該方法最早應用于模擬飛機機翼的空氣動力學，Hess和Smith[1]首次采用面元法來計算無升力體的流體動力，為面元法理論奠定了基礎。與葉素動量法和計算流體力學方法相比，面元法能在較少的時間內獲得較高的計算精度[2]。作為一種數值求解方法，面元法可以真實地模擬葉片的幾何形狀，并在表面上滿足邊界條件，將問題簡化為計算物體表面奇點單元的強度，不需要對整個流場進行求解，從計算的角度來看具有更高的經濟性和實用性[3]。近年來有不少學者將面元法拓展到風力機葉片的氣動性能計算上，并且成功得到了應用。Bermúdez等[4]開發了基于三維低階非定常速度勢的面元法求解器，對NREL水平軸風力機進行了數值模擬，和試驗測試的結果對比表明該求解器可以很好地模擬非定常效應并作為工程設計軟件使用。Tescione等[5]基于自由尾渦/面元法模擬了垂直軸風力機的三維氣動性能，通過和垂直軸風機試驗數據對比，驗證了自由渦尾流模型的正確性，并且該數值模型在已經觀察到的中心收縮和側向膨脹的情況下，可以很好地捕獲水平面中渦旋結構的動力學行為以及尖端渦旋的垂直運動。劉洋等[6]基于三維速度勢面元法來估算風機的三維非定常氣動力，計算值和試驗數據對比良好，證明該方法在貼體流動和小分離工況下具有較好的適用性。周捍瓏等[7]以面元法理論為基礎，MEXICO風輪模型和試驗數據為支撐，對風力機的偏航性能進行了預測。近年來，海上浮式風力機發展迅速，并逐步商業化。由于海洋環境復雜，且浮式基礎運動將影響風力機的氣動載荷，因而海上浮式風力機氣動載荷研究更加重要。Netzband等[8]將面元法應用到海上浮式風力發電機的氣動力和水動力計算，針對浮式基礎的運動展開研究，風機的定常以及非定常模擬結果和3種不同的雷諾平均N-S方程(RANS)計算結果具有一致性。

基于非定常面元法，重點研究海上浮式風力機氣動載荷特性，開發了水平軸風機定常及非定常氣動載荷的計算程序，考慮剪切風、塔影效應以及浮式基礎運動，進一步研究了浮式基礎運動對風力機氣動載荷的影響。

1 面元法基本理論

1.1 基本方程

風力機葉片的表面邊界SB是已知的，如圖1所示，當風力機處于流場中時，總的速度勢符合拉普拉斯方程，用不可壓縮無旋連續性方程來表示：

圖1 計算空間的流場示意Fig. 1 Potential flow of computational domain

(1)

根據格林公式，在物面SB上分布一系列的源σ和偶極子μ來構造求解方程。源用于模擬風力機葉片厚度帶來的影響，而反對稱項例如偶極子(或渦)則用于解決升力問題。對于一般的三維勢流情況，還需要考慮尾流的建模、尾流脫落初始方向和幾何形狀。因此在尾流上也布置偶極子，這樣可以將方程轉化成積分的形式[9]：

(2)

式中:n為指向物體內部的法向量；r為點(x,y,z)到奇點面元的距離；Φ∞為無窮來流的速度勢。

內外速度勢之間的差異定義為偶極子的強度：

(3)

內外速度勢在法向上求導的差值定義為源的強度：

(4)

1.2 邊界條件

直接在表面邊界上指定法向速度為0的方程稱為Neumann邊界條件，而在邊界上給定速度勢的大小間接使得法向速度分量等于0，這樣的邊界條件稱為Dirichlet邊界條件?；谒俣葎莸拿嬖ㄟM行分析，故采用Dirichlet邊界條件，物面SB內部的速度勢是常量：

(5)

物體內部存在的速度勢為：

(6)

上述方程是間接滿足邊界條件的基礎(其中SW表示尾流區域)，但是該常數取不同的值也會產生不同的解法，例如假設該常數等于0，速度勢函數將會包括Φ∞這一項，奇點的強度就會特別大。這里不妨假設該常數等于Φ∞，這樣結合方程(5)和方程(6)，問題簡化為擾動產生的速度勢為0。

(7)

σ=n·V∞

(8)

為了進一步求解SB上未知的偶極子強度和尾流上的偶極子分布，將物體表面SB離散成若干個單元，在離散時，采用四邊形面元，每個面元上均勻分布常值面源和面偶極子，取面元的中心為控制點(見圖2)，物面上的面元數量設為N，尾流的面元數量設為NW，在每一時刻方程(7)離散為一組關于μ和μl的代數方程：

圖2 葉片表面的面元分布Fig. 2 Panel distribution on blade surface

(9)

其中，Ck、Cl和Bk稱為影響系數，采用面元局部坐標系(圖3)進行計算，使用它們與強度的乘積來表示奇點產生的速度勢：

圖3 面元坐標系示意Fig. 3 Schematic diagram of panel coordinate system

(10)

(11)

2 數值離散模型

2.1 葉片及尾流模型

選取NREL 5 MW風機[10]作為計算對象，該風機是三葉片水平軸風力機，葉型采用DU和NACA系列翼型。根據葉片的幾何參數進行建模并劃分網格，為了提高計算的準確性，在弦長方向上劃分30個網格，展長方向網格數為50，最后一共得到1 500個面元，如圖4所示。

圖4 葉片網格劃分示意Fig. 4 Schematic diagram of blade meshing

假設風速為V0，風機的旋轉速度為Ω，r表示計算點到風機旋轉中心的距離，則葉片上的來流速度為[11]：

V∞=V0-Ω×r

(12)

由于影響系數Ck和Bk只和葉片的表面形狀有關，因此不需要重復計算，但是在每個時間步都要重新計算來流速度V∞，根據式(8)更新表面的源強度大小。

針對非定常面元法，采用自由尾流模型，假定尾流是自由運動的，偶極子單元按照當地的流動速度進行移動，尾流上的偶極子分布[12]如圖5所示。因此需要在每個時間步計算尾流處的誘導速度，從而確定尾流的幾何形狀，單元的位移增量為：

圖5 尾流上的偶極子分布[12]Fig. 5 Doublet distribution on wake

(Δx,Δy,Δz)=(u,v,w)·Δt

(13)

同時在每個時刻脫落偶極子的強度不同，每經過一個時間步長，某一尾流面元上的偶極子向后一個面元處轉移，每條尾渦帶上第一個面元處的偶極子強度又等于葉片的定常尾流速度勢。根據該理論可以獲得葉片當前時刻步和下個時刻步的物理量變化關系[13]。圖6給出了采用自由尾流模型，當浮式基礎做縱搖運動時生成的尾流形狀。

圖6 自由尾流模型Fig. 6 Free wake model

每一時間步自動生成尾流面網格形式，確定了尾流面的形狀后尾流影響系數也隨之確定，利用Morino庫塔條件來得到尾流面上的偶極子強度的初值，每條尾渦帶上第一個偶極子單元的大小可以用未知量的形式來表示：

ΓW(j)=μI B,j-μ1,j

(14)

式中:j為不同的尾渦條帶序號，IB為弦向網格數量(面元控制點編號從葉片下翼面后緣變到前緣再到上翼面后緣逐漸增加),μI B,j和μ1,j分別表示葉片后緣上、下表面偶極子強度。

2.2 等壓庫塔條件及黏性修正

采用Morino庫塔條件只是得到計算的初值，但是對風機葉片這種三維翼型來說會產生較大的誤差，還需要逐步迭代滿足等壓庫塔條件，保證尾緣上下表面的壓力相等。

Δp=pupper-plower=0

(15)

在每個面元上的壓力系數為：

(16)

式中:v和p分別表示當地的速度和壓力；vref和pref分別表示來流速度和壓力；μ為葉片表面的速度勢。

進一步對壓力積分可以得到作用在風機上的推力和轉矩，面積為ΔS的單個面元上沿nk方向上受到的力為：

(17)

由于面元法研究的是理想流體的勢流問題，無法考慮流體的黏性產生的阻力,通過Prandtl-Schlichting表面摩擦阻力公式[14]對推力和轉矩進行修正，該經驗公式可以計入形狀的影響，提高計算的準確度。葉片表面的摩擦阻力系數為：

(18)

式中:t表示葉片截面處的最大厚度，C表示葉片截面處翼型的弦長，Re表示雷諾數。推力和轉矩黏性修正項的表達式分別為：

(19)

(20)

式中:Vm表示面元m上的表面切向速度，在葉片坐標系下的3個分量分別為(Vm,i,Vm,j,Vm,k)；Cf,m表示面元m上的摩擦阻力系數；NC和NR分別表示弦長方向和展長方向的網格劃分數量；Z表示葉片數量；(x,y,z)代表面元m控制點的坐標；ΔSm表示面元m的面積。

3 浮式基礎運動影響

當浮式基礎產生6自由度運動時，不僅會改變作用在葉片上的風速方向，平臺自身的運動速度也會對來流速度做出貢獻，因此需要定義坐標轉換關系來描述平臺運動帶來的影響。在描述過程中一共有3個坐標系，分別是：大地坐標系∑o1x1y1z1(坐標系1)，固定于平臺并隨平臺轉動的參考坐標系∑o2x2y2z2(坐標系2)，隨風機葉片轉動的葉片局部坐標系∑o3x3y3z3(坐標系3)。坐標系示意如圖7所示。

圖7 坐標系示意Fig. 7 Coordinate system diagram

根據坐標的轉換關系，如浮式基礎繞重心旋轉(φ,θ,ψ)，坐標系1到坐標系2的轉換矩陣形式則為：

(21)

坐標系1到坐標系2的轉動角度即為橫搖角θroll、縱搖角θpitch和艏搖角θyaw。將風速和平臺的運動速度疊加后得到任一時刻葉片在坐標系2下的來流速度：

Vrel2(t)=f(θroll,θpitch,θyaw)vwind(t)+vP(t)+ωP×rp

(22)

其中，vwind為風速，vP和ωP分別為平臺的平動和轉動速度，rp為葉片面元到平臺重心的矢量。

由于葉片在旋轉過程中方位角的變化以及在安裝時也會導致坐標系2和坐標系3之間存在角度的偏轉，設2、3坐標系之間旋轉的角度為(p,q,r)，經過坐標變換后，來流速度變為：

Vrel3(t)=f(p,q,r)Vrel2(t)+ωo×ro

(23)

其中，ωo表示葉片的旋轉速度，ro表示葉片面元到旋轉中心的矢量，將坐標系3下得到的速度代入式(8)中即可得到每一時刻步面元上的源強度大小。通過該過程將平臺的運動和風速共同處理為非定常來流進行分析，建立了海上浮式風機的非定常面元法計算模型，為下一步的非定常氣動性能分析奠定了基礎。

4 剪切風和塔影效應

采用指數率來描述風速隨高度的變化規律，考慮剪切風的影響時，風速的表達式為[15]：

(24)

其中，v(z)表示高度z位置的風速，vh表示輪轂處的風速，h表示輪轂高度，α表示風剪切指數，取值范圍0.10～0.25之間。在每一時刻計算葉片旋轉到不同位置處每個面元的高度z，然后將vshear賦予式(22)中的vwind，采用非定常面元法也能適用于求解復雜來流的情況。

由于塔架的存在會導致塔架后面的風速降低，采用潛流模型[16]來描述塔影效應，方位角為120°～240°之間的區域影響因子為：

(25)

其中，D表示面元所在高度的塔柱直徑，x和y分別表示面元到塔柱中心的坐標分量?？紤]塔影效應的情況下，將風速乘以影響因子賦予式(22)中的vwind。

浮式基礎不運動時，葉片上風速大小和方向的變化只和方位角有關，不同葉片的功率之間只存在相位的差距[17]，因此多葉片風機不需要做重復的計算，對于三葉片風機，總的功率為；

P=P1+P2+P3=P1(θ)+P1(θ+120°)+P1(θ+240°)

(26)

5 計算程序及驗證

根據非定常面元法的基本理論和計算模型，編寫了海上浮式風機的氣動載荷計算程序，計算流程圖如圖8所示。

圖8 計算流程圖Fig. 8 Flow diagram of calculation

該程序同樣可以用來處理非定常流。在有限元軟件中對葉片進行網格劃分，根據上一時刻計算的誘導速度獲取尾流面的位置信息，然后采用Newton-Rapshon方法進行迭代直到壓力符合等壓庫塔條件，隨后輸出氣動載荷，然后進入下一時間步的計算。程序對葉片的形狀沒有要求，不依賴已有翼型的氣動參數，在后處理中可以顯示整個葉片表面的壓力場分布，并且迭代次數少，采用的算例在每個時間步迭代若干次后即可收斂，時間成本低。

為了驗證面元法計算程序的合理性，給定定常風速，考慮葉片在安裝時的錐角以及傾角，高風速下引進變槳策略，計入槳距角的影響。模擬了不同風速和轉速下，風機的功率和推力，并和文獻[10]中給出的NREL的設計值、文獻[18]中的CFD計算結果以及文獻[19]中的BEM/CFD結果進行對比，如圖9所示。在引入黏性修正后，在額定工況下的輸出功率為5.30 MW，和NREL的設計值(5.29 MW)十分接近。

圖9表明，面元法計算趨勢是正確的，功率計算值偏大推力計算值偏小，推力值和NREL設計值之間存在一些差別，但是和BEM理論的計算值較為接近，主要的原因是文中和文獻[19]都還未考慮塔架的存在。

圖9 功率和推力計算對比Fig. 9 Power and thrust calculation comparison diagram

圖10給出了風速為9 m/s時，在0.45R、0.60R、0.80R、0.90R這幾個截面處的壓力系數分布。從圖10中可以看出，在同一截面上，來流速度以一定的角度流過葉片，然后在前緣點發生分離，一部分沿下表面流動，另一部分繞過前緣點后再沿上表面流動，最后在尾流點匯集，在后緣點上下表面的壓力系數值相等滿足等壓庫塔條件。葉片壓力系數較大的區域集中在葉片0.3R長度附近，并且在吸力面前緣會出現較大的峰值，因為吸力面前緣是葉片表面形狀變化最劇烈的區域，面元法的速度是通過速度勢的局部求導獲得的，曲率過大會導致計算得到的速度也偏大，而在真實的流場中，速度過大時會在黏性作用下被很快的耗散掉[20]，這也是計算功率偏大的主要原因。

圖10 不同截面處的壓力系數和整個葉片上的壓力系數分布云圖Fig. 10 Pressure coefficient at different sections and cloud diagram of pressure coefficient distribution on the entire blade

6 風力機氣動性能分析

6.1 剪切風和塔影效應對風力機功率的影響

對浮式風力機來說，隨著方位角的變化，由于剪切風和塔影效應的影響，導致風速的大小和方向隨時間產生周期性變化，使風機的輸出功率產生波動。

圖11給出了在額定風速11.4 m/s下單個葉片的輸出功率，從圖中可以看出剪切風會導致功率出現簡諧波動，波動的范圍在1.42～2.01 MW之間，而塔影效應會導致在方位角達到180°時降到最低值1.55 MW，兩者聯合作用時，下降的幅度更大，最小的功率值為1.26 MW。相比剪切風，塔影效應對風機輸出功率的穩定性影響更大，當風機轉到塔架位置時會導致功率突然降低隨后又突然上升，呈V字形變化，而總的輸出頻率是風機旋轉頻率的三倍(見圖11)。圖12總結了幾種不同模型下總功率的變化趨勢。以上分析表明：剪切風的影響機理是使功率產生簡諧波動并且降低輸出功率，塔影效應產生的主要影響是當葉片旋轉到塔柱位置時產生功率的驟降，兩者共同作用下，波動幅值比塔影效應單獨作用時要小，剪切風能適當降低塔影效應帶來的不穩定波動，氣動載荷出現風輪旋轉的3P頻率成分。從結果來看非定常面元法能很好地模擬風機的非定常效應，反映入流風速變化對氣動力帶來的影響。

圖11 剪切風和塔影效應對功率的影響Fig. 11 Effect of shear wind and tower shadow on power

圖12 總功率隨方位角的變化Fig.12 Variation of total power with azimuth angle

6.2 浮式基礎運動對氣動載荷的影響

為了研究浮式基礎運動對風機氣動載荷的影響，且為了避免變槳工況，設定風機在風速9 m/s下運行，轉速為10.33 r/min，計算時間步長為0.1 s。浮式基礎縱蕩速度和風速在一個方向上，直接影響來流風速的大小，而縱搖和艏搖會改變風速的方向，使風機處于偏航狀態，因此選取了對風機氣動性能影響顯著的3個自由度的運動：縱蕩、縱搖和艏搖，將浮式基礎這3個自由度的運動簡化為正弦運動，運動周期和幅值見表1，數值的選取參考了真實海況下浮式基礎的運動[21-23]。當葉片旋轉到不同的角度時，由于浮式基礎的轉動角度和提供的附加速度大小都不相同，不能簡單的轉換相位進行疊加，采用非定常面元法的計算思路是在計算不同的葉片時改變初始方位角，具體體現在式(23)中的角度p上，從而導致轉換矩陣也各不相同，計算得到不同葉片上的氣動載荷隨后進行疊加，最后得到總的功率和推力隨時間變化的結果，如圖13～17所示。

表1 平臺運動工況

圖13 浮式基礎縱蕩對總功率和推力的影響Fig. 13 The influence of surge on floating foundation on total power and thrust

以上計算表明：浮式基礎的縱蕩和縱搖對風力機氣動載荷的影響較大，當浮式基礎在這兩個自由度做正弦運動時，風力機總功率和推力的變化周期和浮式基礎的運動周期一致。浮式基礎發生艏搖運動時，風力機的總功率變化周期和葉片的旋轉周期一致，變化較為平穩，但是單個葉片上氣動載荷的波動較為復雜，每個葉片上氣動載荷的變化周期等于艏搖運動的周期，同時每個葉片之間不止存在相位角的差距，氣動載荷的幅值和變化趨勢都不相同。

此外，計算結果還表明，浮式基礎縱搖和縱蕩對葉片氣動載荷的波動影響最大，不考慮平臺運動時算得的功率值為2.78 MW，而縱搖會使得風機功率的范圍從1.15 MW變化到4.84 MW。因此，浮式基礎縱搖顯著影響風機輸出功率的穩定性，實際中應對其加以必要的控制。在浮式基礎縱蕩的影響下，文中算例中，使得風機功率的變化范圍從2.04 MW變化到3.43 MW，雖然縱蕩運動不會影響來流風速的方向，但是會提供額外的附加速度，當平臺縱蕩運動速度過大時會顯著改變風輪輸出功率的波動范圍；浮式基礎艏搖對風機總功率波動范圍的影響非常小，風機功率的波動范圍在0.05 MW以內。

圖14 浮式基礎縱搖時的總功率和推力Fig. 14 Total power and thrust of pitch

圖15 浮式基礎艏搖時總功率和推力Fig.15 Total power and thrust of yaw

圖16 浮式基礎縱搖時各個葉片上的轉矩Fig. 16 Moment on each blade of pitch

圖17 浮式基礎艏搖時各個葉片上的轉矩Fig. 17 Moment on each blade of yaw

根據以上分析可知，適當的浮式基礎運動會提高風機的輸出功率，但是會導致風機輸出功率的不穩定，而葉片之間相位角的差距能適當彌補這種波動帶來的不穩定性。

圖18進一步給出了計算穩定后，不同時刻風機葉片表面的壓力場分布云圖?？梢钥闯?，在0.5R～0.9R長度范圍內，風機葉片吸力面前緣部分壓力最大，壓力的最大絕對值先變大再減??；壓力面上的壓力隨時間變化不明顯，同時經過非線性迭代后，在后緣上下表面的壓力大小相同，并且在該區域范圍內壓力過渡平滑，說明滿足非定常等壓庫塔條件。

7 結 語

基于三維速度勢的非定常面元法，開發了浮式水平軸風力機的氣動載荷計算程序，以NREL 5 MW風機為例，分析了剪切風、塔影效應以及浮式基礎運動對風機氣動載荷的影響。研究結論如下：

1) 研究了剪切風和塔影效應對風機氣動性能的影響。結果表明，考慮剪切風和塔影效應時，單個葉片上氣動載荷的變化頻率等于風機的轉動頻率，塔影效應對風機功率的穩定性影響顯著，風機整體在方位角位于60°、180°和300°時會發生輸出功率的突然下降。

2) 分析了浮式基礎縱蕩、縱搖和艏搖運動對風力機氣動性能的影響。結果表明，總的輸出功率和推力的變化周期以浮式基礎運動的周期為主。浮式基礎縱蕩和縱搖對風機整體功率波動幅值的影響很大，不利于海上風機的穩定運行；而對單個葉片而言，其氣動載荷對浮式基礎的艏搖會更加敏感。

3) 分析了風機葉片壓力分布特征。結果表明，葉片壓力系數較大的區域集中在葉片0.3R長度附近，壓力較大且變化比較劇烈的區域在0.5R～0.9R展長處的葉片前緣附近，在風機葉片結構設計時應對該區域的強度給與額外的關注。