基于非定常面元法的海上浮式風(fēng)機(jī)氣動(dòng)性能研究

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

(1. 天津大學(xué) 水利工程仿真與安全國家重點(diǎn)實(shí)驗(yàn)室，天津 300072； 2. 交通運(yùn)輸部天津水運(yùn)工程科學(xué)研究所 工程泥沙交通行業(yè)重點(diǎn)實(shí)驗(yàn)室，天津 300456)

葉片是風(fēng)力機(jī)的關(guān)鍵部件，其表面形狀對整個(gè)風(fēng)機(jī)的空氣動(dòng)力學(xué)性能和風(fēng)力機(jī)的風(fēng)能利用率都會(huì)產(chǎn)生重要的影響，因此葉片的氣動(dòng)性能一直都是風(fēng)力機(jī)的研究重點(diǎn)。為評估風(fēng)力發(fā)電機(jī)效率和結(jié)構(gòu)安全，需要正確預(yù)估不同來流條件下葉片上的風(fēng)載荷。特別是在研究海上浮式風(fēng)機(jī)時(shí)，自身的繞流情況會(huì)變得更為復(fù)雜。為了提高風(fēng)能的利用效率，必須要準(zhǔn)確快速地解決風(fēng)機(jī)面臨的空氣動(dòng)力學(xué)問題。

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

基于非定常面元法，重點(diǎn)研究海上浮式風(fēng)力機(jī)氣動(dòng)載荷特性，開發(fā)了水平軸風(fēng)機(jī)定常及非定常氣動(dòng)載荷的計(jì)算程序，考慮剪切風(fēng)、塔影效應(yīng)以及浮式基礎(chǔ)運(yùn)動(dòng)，進(jìn)一步研究了浮式基礎(chǔ)運(yùn)動(dòng)對風(fēng)力機(jī)氣動(dòng)載荷的影響。

1 面元法基本理論

1.1 基本方程

風(fēng)力機(jī)葉片的表面邊界SB是已知的，如圖1所示，當(dāng)風(fēng)力機(jī)處于流場中時(shí)，總的速度勢符合拉普拉斯方程，用不可壓縮無旋連續(xù)性方程來表示：

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

(1)

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

(2)

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

內(nèi)外速度勢之間的差異定義為偶極子的強(qiáng)度：

(3)

內(nèi)外速度勢在法向上求導(dǎo)的差值定義為源的強(qiáng)度：

(4)

1.2 邊界條件

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

(5)

物體內(nèi)部存在的速度勢為：

(6)

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

(7)

σ=n·V∞

(8)

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

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

(9)

其中，Ck、Cl和Bk稱為影響系數(shù)，采用面元局部坐標(biāo)系(圖3)進(jìn)行計(jì)算，使用它們與強(qiáng)度的乘積來表示奇點(diǎn)產(chǎn)生的速度勢：

圖3 面元坐標(biāo)系示意Fig. 3 Schematic diagram of panel coordinate system

(10)

(11)

2 數(shù)值離散模型

2.1 葉片及尾流模型

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

圖4 葉片網(wǎng)格劃分示意Fig. 4 Schematic diagram of blade meshing

假設(shè)風(fēng)速為V0，風(fēng)機(jī)的旋轉(zhuǎn)速度為Ω，r表示計(jì)算點(diǎn)到風(fēng)機(jī)旋轉(zhuǎn)中心的距離，則葉片上的來流速度為[11]：

V∞=V0-Ω×r

(12)

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

針對非定常面元法，采用自由尾流模型，假定尾流是自由運(yùn)動(dòng)的，偶極子單元按照當(dāng)?shù)氐牧鲃?dòng)速度進(jìn)行移動(dòng)，尾流上的偶極子分布[12]如圖5所示。因此需要在每個(gè)時(shí)間步計(jì)算尾流處的誘導(dǎo)速度，從而確定尾流的幾何形狀，單元的位移增量為：

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

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

(13)

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

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

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

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

(14)

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

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

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

Δp=pupper-plower=0

(15)

在每個(gè)面元上的壓力系數(shù)為：

(16)

式中:v和p分別表示當(dāng)?shù)氐乃俣群蛪毫Γ籿ref和pref分別表示來流速度和壓力；μ為葉片表面的速度勢。

進(jìn)一步對壓力積分可以得到作用在風(fēng)機(jī)上的推力和轉(zhuǎn)矩，面積為ΔS的單個(gè)面元上沿nk方向上受到的力為：

(17)

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

(18)

式中:t表示葉片截面處的最大厚度，C表示葉片截面處翼型的弦長，Re表示雷諾數(shù)。推力和轉(zhuǎn)矩黏性修正項(xiàng)的表達(dá)式分別為：

(19)

(20)

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

3 浮式基礎(chǔ)運(yùn)動(dòng)影響

當(dāng)浮式基礎(chǔ)產(chǎn)生6自由度運(yùn)動(dòng)時(shí)，不僅會(huì)改變作用在葉片上的風(fēng)速方向，平臺(tái)自身的運(yùn)動(dòng)速度也會(huì)對來流速度做出貢獻(xiàn)，因此需要定義坐標(biāo)轉(zhuǎn)換關(guān)系來描述平臺(tái)運(yùn)動(dòng)帶來的影響。在描述過程中一共有3個(gè)坐標(biāo)系，分別是：大地坐標(biāo)系∑o1x1y1z1(坐標(biāo)系1)，固定于平臺(tái)并隨平臺(tái)轉(zhuǎn)動(dòng)的參考坐標(biāo)系∑o2x2y2z2(坐標(biāo)系2)，隨風(fēng)機(jī)葉片轉(zhuǎn)動(dòng)的葉片局部坐標(biāo)系∑o3x3y3z3(坐標(biāo)系3)。坐標(biāo)系示意如圖7所示。

圖7 坐標(biāo)系示意Fig. 7 Coordinate system diagram

根據(jù)坐標(biāo)的轉(zhuǎn)換關(guān)系，如浮式基礎(chǔ)繞重心旋轉(zhuǎn)(φ,θ,ψ)，坐標(biāo)系1到坐標(biāo)系2的轉(zhuǎn)換矩陣形式則為：

(21)

坐標(biāo)系1到坐標(biāo)系2的轉(zhuǎn)動(dòng)角度即為橫搖角θroll、縱搖角θpitch和艏搖角θyaw。將風(fēng)速和平臺(tái)的運(yùn)動(dòng)速度疊加后得到任一時(shí)刻葉片在坐標(biāo)系2下的來流速度：

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

(22)

其中，vwind為風(fēng)速，vP和ωP分別為平臺(tái)的平動(dòng)和轉(zhuǎn)動(dòng)速度，rp為葉片面元到平臺(tái)重心的矢量。

由于葉片在旋轉(zhuǎn)過程中方位角的變化以及在安裝時(shí)也會(huì)導(dǎo)致坐標(biāo)系2和坐標(biāo)系3之間存在角度的偏轉(zhuǎn)，設(shè)2、3坐標(biāo)系之間旋轉(zhuǎn)的角度為(p,q,r)，經(jīng)過坐標(biāo)變換后，來流速度變?yōu)椋?/p>

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

(23)

其中，ωo表示葉片的旋轉(zhuǎn)速度，ro表示葉片面元到旋轉(zhuǎn)中心的矢量，將坐標(biāo)系3下得到的速度代入式(8)中即可得到每一時(shí)刻步面元上的源強(qiáng)度大小。通過該過程將平臺(tái)的運(yùn)動(dòng)和風(fēng)速共同處理為非定常來流進(jìn)行分析，建立了海上浮式風(fēng)機(jī)的非定常面元法計(jì)算模型，為下一步的非定常氣動(dòng)性能分析奠定了基礎(chǔ)。

4 剪切風(fēng)和塔影效應(yīng)

采用指數(shù)率來描述風(fēng)速隨高度的變化規(guī)律，考慮剪切風(fēng)的影響時(shí)，風(fēng)速的表達(dá)式為[15]：

(24)

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

由于塔架的存在會(huì)導(dǎo)致塔架后面的風(fēng)速降低，采用潛流模型[16]來描述塔影效應(yīng)，方位角為120°～240°之間的區(qū)域影響因子為：

(25)

其中，D表示面元所在高度的塔柱直徑，x和y分別表示面元到塔柱中心的坐標(biāo)分量。考慮塔影效應(yīng)的情況下，將風(fēng)速乘以影響因子賦予式(22)中的vwind。

浮式基礎(chǔ)不運(yùn)動(dòng)時(shí)，葉片上風(fēng)速大小和方向的變化只和方位角有關(guān)，不同葉片的功率之間只存在相位的差距[17]，因此多葉片風(fēng)機(jī)不需要做重復(fù)的計(jì)算，對于三葉片風(fēng)機(jī)，總的功率為；

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

(26)

5 計(jì)算程序及驗(yàn)證

根據(jù)非定常面元法的基本理論和計(jì)算模型，編寫了海上浮式風(fēng)機(jī)的氣動(dòng)載荷計(jì)算程序，計(jì)算流程圖如圖8所示。

圖8 計(jì)算流程圖Fig. 8 Flow diagram of calculation

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

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

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

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

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

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

6 風(fēng)力機(jī)氣動(dòng)性能分析

6.1 剪切風(fēng)和塔影效應(yīng)對風(fēng)力機(jī)功率的影響

對浮式風(fēng)力機(jī)來說，隨著方位角的變化，由于剪切風(fēng)和塔影效應(yīng)的影響，導(dǎo)致風(fēng)速的大小和方向隨時(shí)間產(chǎn)生周期性變化，使風(fēng)機(jī)的輸出功率產(chǎn)生波動(dòng)。

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

圖11 剪切風(fēng)和塔影效應(yīng)對功率的影響Fig. 11 Effect of shear wind and tower shadow on power

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

6.2 浮式基礎(chǔ)運(yùn)動(dòng)對氣動(dòng)載荷的影響

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

表1 平臺(tái)運(yùn)動(dòng)工況

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

以上計(jì)算表明：浮式基礎(chǔ)的縱蕩和縱搖對風(fēng)力機(jī)氣動(dòng)載荷的影響較大，當(dāng)浮式基礎(chǔ)在這兩個(gè)自由度做正弦運(yùn)動(dòng)時(shí)，風(fēng)力機(jī)總功率和推力的變化周期和浮式基礎(chǔ)的運(yùn)動(dòng)周期一致。浮式基礎(chǔ)發(fā)生艏搖運(yùn)動(dòng)時(shí)，風(fēng)力機(jī)的總功率變化周期和葉片的旋轉(zhuǎn)周期一致，變化較為平穩(wěn)，但是單個(gè)葉片上氣動(dòng)載荷的波動(dòng)較為復(fù)雜，每個(gè)葉片上氣動(dòng)載荷的變化周期等于艏搖運(yùn)動(dòng)的周期，同時(shí)每個(gè)葉片之間不止存在相位角的差距，氣動(dòng)載荷的幅值和變化趨勢都不相同。

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

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

圖15 浮式基礎(chǔ)艏搖時(shí)總功率和推力Fig.15 Total power and thrust of yaw

圖16 浮式基礎(chǔ)縱搖時(shí)各個(gè)葉片上的轉(zhuǎn)矩Fig. 16 Moment on each blade of pitch

圖17 浮式基礎(chǔ)艏搖時(shí)各個(gè)葉片上的轉(zhuǎn)矩Fig. 17 Moment on each blade of yaw

根據(jù)以上分析可知，適當(dāng)?shù)母∈交A(chǔ)運(yùn)動(dòng)會(huì)提高風(fēng)機(jī)的輸出功率，但是會(huì)導(dǎo)致風(fēng)機(jī)輸出功率的不穩(wěn)定，而葉片之間相位角的差距能適當(dāng)彌補(bǔ)這種波動(dòng)帶來的不穩(wěn)定性。

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

7 結(jié) 語

基于三維速度勢的非定常面元法，開發(fā)了浮式水平軸風(fēng)力機(jī)的氣動(dòng)載荷計(jì)算程序，以NREL 5 MW風(fēng)機(jī)為例，分析了剪切風(fēng)、塔影效應(yīng)以及浮式基礎(chǔ)運(yùn)動(dòng)對風(fēng)機(jī)氣動(dòng)載荷的影響。研究結(jié)論如下：

1) 研究了剪切風(fēng)和塔影效應(yīng)對風(fēng)機(jī)氣動(dòng)性能的影響。結(jié)果表明，考慮剪切風(fēng)和塔影效應(yīng)時(shí)，單個(gè)葉片上氣動(dòng)載荷的變化頻率等于風(fēng)機(jī)的轉(zhuǎn)動(dòng)頻率，塔影效應(yīng)對風(fēng)機(jī)功率的穩(wěn)定性影響顯著，風(fēng)機(jī)整體在方位角位于60°、180°和300°時(shí)會(huì)發(fā)生輸出功率的突然下降。

2) 分析了浮式基礎(chǔ)縱蕩、縱搖和艏搖運(yùn)動(dòng)對風(fēng)力機(jī)氣動(dòng)性能的影響。結(jié)果表明，總的輸出功率和推力的變化周期以浮式基礎(chǔ)運(yùn)動(dòng)的周期為主。浮式基礎(chǔ)縱蕩和縱搖對風(fēng)機(jī)整體功率波動(dòng)幅值的影響很大，不利于海上風(fēng)機(jī)的穩(wěn)定運(yùn)行；而對單個(gè)葉片而言，其氣動(dòng)載荷對浮式基礎(chǔ)的艏搖會(huì)更加敏感。

3) 分析了風(fēng)機(jī)葉片壓力分布特征。結(jié)果表明，葉片壓力系數(shù)較大的區(qū)域集中在葉片0.3R長度附近，壓力較大且變化比較劇烈的區(qū)域在0.5R～0.9R展長處的葉片前緣附近，在風(fēng)機(jī)葉片結(jié)構(gòu)設(shè)計(jì)時(shí)應(yīng)對該區(qū)域的強(qiáng)度給與額外的關(guān)注。