周 胡，趙文超，萬(wàn)德成

(上海交通大學(xué) 船舶海洋與建筑工程學(xué)院 海洋工程國(guó)家重點(diǎn)實(shí)驗(yàn)室 高新船舶與海洋開(kāi)發(fā)裝備協(xié)同創(chuàng)新中心，上海 200240)

在可再生能源領(lǐng)域中，風(fēng)能的開(kāi)發(fā)和利用技術(shù)是最為成熟也是最具大規(guī)模商業(yè)開(kāi)發(fā)前景的，受到了世界各國(guó)的廣泛關(guān)注。特別是最近幾年海上風(fēng)能發(fā)電的興起，再一次掀起了研究和利用風(fēng)能的熱潮。風(fēng)力機(jī)將風(fēng)能轉(zhuǎn)化為機(jī)械能，風(fēng)資源的品質(zhì)對(duì)風(fēng)力機(jī)的功率、壽命和運(yùn)行等起著重要影響。風(fēng)力機(jī)的來(lái)流風(fēng)由于受地球表面的影響，最顯著的特征是在時(shí)間或空間上的分布不均。這種非均勻的來(lái)流風(fēng)可以用風(fēng)速廓線的形式予以量化表示。由于風(fēng)速廓線引起的風(fēng)切變效應(yīng)的存在，葉片在旋轉(zhuǎn)過(guò)程中將經(jīng)歷風(fēng)速的周期性變化，導(dǎo)致風(fēng)輪受到風(fēng)力載荷的不均，這種不均會(huì)對(duì)風(fēng)力機(jī)的使用壽命和運(yùn)行安全產(chǎn)生深遠(yuǎn)影響。同時(shí)，在未來(lái)海上風(fēng)電發(fā)展過(guò)程中，為了降低成本，葉輪直徑會(huì)隨著額定功率的增加而增加，風(fēng)切變對(duì)風(fēng)力機(jī)功率的影響將更為重要，必須在設(shè)計(jì)時(shí)引起注意。

隨著計(jì)算機(jī)性能和計(jì)算方法的飛速提高和發(fā)展，使用計(jì)算流體力學(xué)理論數(shù)值模擬風(fēng)力機(jī)空氣動(dòng)力學(xué)性能的方法日益受到了重視，數(shù)值的方法不僅能夠準(zhǔn)確預(yù)測(cè)風(fēng)力機(jī)性能，還能夠清楚地觀察流場(chǎng)細(xì)節(jié)，深化對(duì)繞流場(chǎng)流動(dòng)問(wèn)題本質(zhì)的理解。使用基于雷諾時(shí)均Navier－Stokes 方程求解風(fēng)力機(jī)繞流場(chǎng)的研究已有很多［1－4］，但是使用數(shù)值方法模擬非均勻風(fēng)或剪切風(fēng)作用下風(fēng)力機(jī)空氣動(dòng)力性能的工作還不是很多［5－6］，這也正是本文的立足點(diǎn)。

Nilay Sezer－Uzol 等［5］使用基于自由渦流理論的勢(shì)流方法研究了穩(wěn)態(tài)和瞬態(tài)風(fēng)剪切對(duì)水平軸風(fēng)力機(jī)的性能特別是渦結(jié)構(gòu)的影響，文中著重研究了三種典型工況，分別是均勻風(fēng)、指數(shù)率的風(fēng)速廓線及瞬態(tài)極限剪切風(fēng)，研究表明風(fēng)速廓線的存在對(duì)風(fēng)力機(jī)的尾渦結(jié)構(gòu)和葉片表面的壓力將產(chǎn)生很大非對(duì)稱影響。廖明夫等［7］基于對(duì)數(shù)律分布的風(fēng)切變模型從理論上對(duì)非均勻風(fēng)對(duì)風(fēng)力機(jī)功率影響進(jìn)行了詳細(xì)的研究和分析，得出了非均勻風(fēng)引起的風(fēng)力機(jī)設(shè)計(jì)功率的損失與風(fēng)力機(jī)葉輪直徑的大小相關(guān)等有用結(jié)論。劉磊等［6］依據(jù)指數(shù)率的風(fēng)速分布曲線，使用了CFD 方法對(duì)非均勻風(fēng)作用下風(fēng)力機(jī)三維非定常氣動(dòng)特性進(jìn)行了研究，同時(shí)還研究了不同風(fēng)切變指數(shù)對(duì)風(fēng)力機(jī)載荷波動(dòng)振幅的影響。文中基于OpenFOAM 開(kāi)源平臺(tái)，選用任意網(wǎng)格界面元法瞬態(tài)求解器進(jìn)行非均勻風(fēng)作用下風(fēng)力機(jī)空氣動(dòng)力性能預(yù)報(bào)方法的研究。

1 OpenFOAM 簡(jiǎn)介

OpenFOAM 全稱為Open Field Operation and Manipulation，即場(chǎng)的操作和處理的開(kāi)源計(jì)算平臺(tái)。它是一個(gè)基于有限體積法，可用于對(duì)連續(xù)介質(zhì)力學(xué)問(wèn)題進(jìn)行數(shù)值計(jì)算的面向?qū)ο蟮腃+ +庫(kù)，同時(shí)還提供了許多預(yù)編譯好的求解器、輔助工具和模型庫(kù)等。由于使用了許多C+ +語(yǔ)言的高級(jí)工具，例如模板類(lèi)、操作符重載、多態(tài)等［8］，OpenFOAM 具有強(qiáng)大的可定制性和可拓展性，使用者和開(kāi)發(fā)者可以根據(jù)所求問(wèn)題的特殊性編寫(xiě)自己的求解器，重點(diǎn)關(guān)注求解的流程，而不需要關(guān)注離散和求解的最底層知識(shí)，這是許多商業(yè)軟件無(wú)法比擬的。將OpenFOAM 作為底層庫(kù)類(lèi)來(lái)構(gòu)建自己的求解器是許多研究者選擇OpenFOAM 作為求解器開(kāi)發(fā)平臺(tái)的重要原因。

上海交通大學(xué)萬(wàn)德成課題組已經(jīng)在OpenFOAM 的開(kāi)發(fā)和使用上做了大量豐富工作。查晶晶、曹洪建等［9－10］基于OpenFOAM 求解器interDyMFoam，開(kāi)發(fā)實(shí)現(xiàn)了數(shù)值粘性水池造波和阻尼消波，對(duì)圓柱波浪爬高等進(jìn)行了廣泛研究。沈志榮等［11］將六自由度運(yùn)動(dòng)模塊植入OpenFOAM，開(kāi)發(fā)了naoe－FOAM－SJTU 求解器，實(shí)現(xiàn)了對(duì)船舶在波浪上運(yùn)動(dòng)的數(shù)值模擬。王強(qiáng)、周胡等［12－14］基于MRFSimpleFoam 和pimpleDyMFoam求解器對(duì)風(fēng)力機(jī)風(fēng)輪的氣動(dòng)力性能、風(fēng)力機(jī)葉片與支撐塔架相互作用的耦合流場(chǎng)進(jìn)行了數(shù)值模擬。

2 計(jì)算模型與網(wǎng)格生成

Phase VI 實(shí)驗(yàn)風(fēng)力機(jī)兩葉片模型是美國(guó)國(guó)家新能源實(shí)驗(yàn)室于2000年在美國(guó)國(guó)家航空航天局Ames 120×80 ft 風(fēng)洞中進(jìn)行系列實(shí)驗(yàn)的模型。由于該模型具有詳細(xì)的實(shí)驗(yàn)數(shù)據(jù)［15］，所以模型常被用來(lái)進(jìn)行風(fēng)力機(jī)葉片空氣動(dòng)力特性數(shù)值模擬結(jié)果的驗(yàn)證。通過(guò)求解RANS 方程結(jié)合k－ω SST 模型對(duì)實(shí)驗(yàn)風(fēng)力機(jī)模型外流場(chǎng)非定常流動(dòng)進(jìn)行數(shù)值計(jì)算。計(jì)算模型為上風(fēng)向風(fēng)力機(jī)，從x 軸正向看，風(fēng)輪逆時(shí)針轉(zhuǎn)動(dòng)，轉(zhuǎn)速為72 轉(zhuǎn)/分鐘。風(fēng)力機(jī)直徑10.058 m，輪轂高度為12.192 m，葉片翼型為NREL S809，詳細(xì)模型數(shù)據(jù)可以參考文獻(xiàn)［16］。

網(wǎng)格的生成流程大致可以分為兩過(guò)程，首先借助ICEM CFD 網(wǎng)格軟件劃分背景網(wǎng)格，再使用OpenFOAM自帶的SnappyHexMesh 工具通過(guò)調(diào)整system 文件夾下的snappyHexMeshDict 字典參數(shù)自動(dòng)生成最終網(wǎng)格。最終劃分的網(wǎng)格情況如圖1 所示，風(fēng)力機(jī)距入口5 m，距出口20 m。在來(lái)流風(fēng)向，靠近風(fēng)力機(jī)前后0.8 m 范圍內(nèi)進(jìn)行了加密，同時(shí)為了捕捉尾渦，風(fēng)力機(jī)后方5 m 內(nèi)也進(jìn)行了加密。最終網(wǎng)格量大概為80 萬(wàn)左右?；平缑?AMI)網(wǎng)格如圖2 所示，交接面的具體生成過(guò)程可以參考下文任意網(wǎng)格界面元法的介紹。

計(jì)算工況中的空氣密度、運(yùn)動(dòng)粘度、轉(zhuǎn)速和槳矩角參考Sequence S 系列實(shí)驗(yàn)，選取輪轂處5、10、15 和25 m/s四個(gè)風(fēng)速進(jìn)行數(shù)值計(jì)算，四個(gè)風(fēng)速下的初始參數(shù)設(shè)置參如表1 所示。

圖1 整體網(wǎng)格情況Fig.1 Global mesh

圖2 AMI 交界面網(wǎng)格Fig.2 Mesh of AMI Interface

表1 實(shí)驗(yàn)各工況參數(shù)Tab.1 Experiment parameter of different computation cases

3 數(shù)值模擬

3.1 控制方程

在慣性參考系中，風(fēng)力機(jī)風(fēng)輪以設(shè)定的角速度繞固定軸轉(zhuǎn)動(dòng)，流場(chǎng)為非定常。由于風(fēng)力機(jī)外流場(chǎng)的風(fēng)速較低，繞流場(chǎng)可以看成不可壓。其時(shí)均控制方程可以表述為:

3.2 湍流模型

為了計(jì)算雷諾應(yīng)力項(xiàng)，需要引入k－ω SST 模型［17］使控制方程封閉。其中k 和ω 分別代表湍動(dòng)能和渦量脈動(dòng)強(qiáng)度。兩者的輸運(yùn)方程為:

3.3 離散方法

OpenFOAM 軟件包基于有限體積法，采用空間網(wǎng)格的形式將計(jì)算區(qū)域劃分為若干的控制體，在每個(gè)控制體上分別求解連續(xù)性方程、動(dòng)量方程和能量方程。

在控制方程各項(xiàng)離散中，時(shí)間項(xiàng)采用隱式歐拉方法，對(duì)流項(xiàng)采用一階高斯迎風(fēng)格式，擴(kuò)散項(xiàng)采用修正的高斯線性方法?？刂品匠讨械乃俣群蛪毫Φ慕怦钍褂肞IMPLE 算法(PISO(pressure－implicit split－operator)和SIMPLE(semi－implicit method for pressure－linked equations)混合算法)，PIMPLE 算法的主要結(jié)構(gòu)是從PISO 中繼承的，主要的區(qū)別表現(xiàn)在對(duì)PISO 算法的每個(gè)時(shí)間步再用SIMPLE 算法求解以得到更穩(wěn)定的解。速度與壓力解耦算法的深入解析可以參考文獻(xiàn)［18］。

3.4 邊界條件

3.4.1 入口邊界條件

風(fēng)速廓線能確定風(fēng)速沿高度的變化規(guī)律，是確定給定高度處風(fēng)力機(jī)輸出功率的前提。風(fēng)速沿高度的變化是風(fēng)的重要特性之一。風(fēng)速廓線一般有兩種描述方式:一種是按邊界層理論得到的對(duì)數(shù)風(fēng)速廓線;一種是按實(shí)測(cè)結(jié)果得到的指數(shù)風(fēng)速廓線。文中采用指數(shù)風(fēng)速廓線來(lái)描述入口的速度分布。

風(fēng)速廓線的指數(shù)分布規(guī)律可以近似表示為［19］:

圖3 輪轂處兩種風(fēng)速下的入口風(fēng)速廓線Fig.3 Wind profile of 5 m/s and 15 m/s at hub location

在OpenFOAM 中，沒(méi)有現(xiàn)成的指數(shù)分布速度邊界條件，需要自己編寫(xiě)。一般有兩種邊界條件編寫(xiě)的方法:其一是將邊界條件加到某一個(gè)特定求解器，在使用求解器的時(shí)候就可以使用編譯好的新的邊界條件;其二是將需要添加的邊界條件編譯成一個(gè)動(dòng)態(tài)鏈接庫(kù)，需要使用該邊界條件的時(shí)候，只需要在OpenFOAM的system 文件夾里的controlDict 中添加相關(guān)語(yǔ)句即可。文中將使用一個(gè)名為windprofileVelocity 類(lèi)型的邊界條件，通過(guò)編寫(xiě)動(dòng)態(tài)鏈接庫(kù)的方法添加。

在windprofileVelocity 邊界條件編譯好后，邊界條件的使用就像使用其它已經(jīng)編譯好的邊界條件一樣。0文件夾下的U 文件里新的邊界設(shè)置如下:

其中，inlet 是入口邊界的名字，n 是來(lái)流風(fēng)速的方向，meanValue 是輪轂處的風(fēng)速，value 為入口速度初始值。最終計(jì)算時(shí)入口的速度分布如圖4 所示，與預(yù)定義的速度廓線一致，說(shuō)明入口邊界條件編寫(xiě)正確。

圖4 輪轂處5 m/s 和15 m/s 速度計(jì)算模型入口速度分布Fig.4 Inlet velocity distribution at speeds of 5 m/s and 15 m/s

入口處的壓力設(shè)為零梯度，湍動(dòng)能k 和渦量脈動(dòng)強(qiáng)度ω 設(shè)為固定值，固定值的取法可根據(jù)經(jīng)驗(yàn)公式求得，渦粘度也通過(guò)計(jì)算求得。

3.4.2 出口邊界條件出口處的速度邊界條件為零梯度，壓力為固定值0，湍動(dòng)能k 和渦量脈動(dòng)強(qiáng)度ω 設(shè)為固定值，湍動(dòng)能k 和渦量脈動(dòng)強(qiáng)度ω 設(shè)為零梯度，渦粘度通過(guò)計(jì)算求得。

3.4.3 壁面條件

計(jì)算域外圍的速度設(shè)為固定值0;壓力設(shè)為零梯度;湍動(dòng)能k、渦量脈動(dòng)強(qiáng)度ω 和渦粘度采用壁面函數(shù)。

3.4.4 風(fēng)力機(jī)模型邊界條件

風(fēng)力機(jī)的塔架和不動(dòng)的輪轂的速度設(shè)為固定值0，葉片和轉(zhuǎn)動(dòng)的輪轂設(shè)為movingWallVelocity;整個(gè)風(fēng)力機(jī)的壓力邊界條件設(shè)為零梯度;湍動(dòng)能k、渦量脈動(dòng)強(qiáng)度ω 和渦粘度采用壁面函數(shù)。

3.5 任意網(wǎng)格界面元法

在OpenFOAM 新的版本中，采用任意網(wǎng)格界面元法來(lái)處理旋轉(zhuǎn)的交界面。任意網(wǎng)格界面元法的本質(zhì)是一種滑移網(wǎng)格技術(shù)。任意網(wǎng)格界面元法和之前版本中的通用網(wǎng)格界面法原理基本類(lèi)似，但是也存在一些細(xì)微區(qū)別。任意網(wǎng)格界面元法主要是通過(guò)插值實(shí)現(xiàn)動(dòng)靜區(qū)域流場(chǎng)參數(shù)和信息的交換。為描述方便，界面兩邊運(yùn)動(dòng)和靜止的面分別定義為主面(master)和從面(slave)。其求解過(guò)程主要步驟如下:

圖5 主從面的關(guān)系示意Fig.5 Sketch map of relation between rotor and stator

3.5.1 創(chuàng)建AMI 交界面

在OpenFOAM 中，創(chuàng)建AMI 交界面可以通過(guò)網(wǎng)格的拓?fù)溥\(yùn)算完成，例如使用topoSet 工具來(lái)選取相應(yīng)的面域。之后需要把該交界面復(fù)制一份，使得一個(gè)面屬于主面，一個(gè)屬于從面，如圖5 所示，圖中rotor side 代表主面，stator side 代表從面。這一步操作可以使用OpenFOAM 中的createBaffles 和splitOrMergeBaffles兩個(gè)工具完成。

3.5.2 尋找相鄰面

當(dāng)AMI 交界面生成后，就可以開(kāi)始計(jì)算了。計(jì)算的第一步便是尋找相鄰面。尋找相鄰面的過(guò)程是已知主面尋找相鄰從面的過(guò)程。

如圖6 所示，運(yùn)動(dòng)部分的主面O 與靜止部分的從面a、b、c 三個(gè)面相鄰，故可得這三個(gè)面為主面O 的相鄰面。具體搜尋相鄰面的算法在此不作介紹。

3.5.3 計(jì)算權(quán)重

權(quán)重的計(jì)算是相互的，即需要計(jì)算從主面到從面的權(quán)重，也需要計(jì)算從從面到主面的權(quán)重。權(quán)重由相重疊區(qū)域占該面的面積來(lái)表示:

主面到從面的權(quán)重:

圖6 相互滑移網(wǎng)格交界面的處理［21］Fig.6 Interaction of different sliding meshes

從面到主面的權(quán)重:

式中:i 為某一個(gè)主面的第i 個(gè)相鄰從面;j 為某一個(gè)從面的第j 個(gè)相鄰主面;為第j 個(gè)主面和從面重疊區(qū)域面積大小，第i 個(gè)從面和主面重疊面積大小;為某一從面第j 個(gè)主面的面積和某一主面第i 個(gè)從面的面積。

例如對(duì)于圖6，想要計(jì)算從面b 到主面O 的權(quán)重，那么可以通過(guò)以上公式計(jì)算得到權(quán)重。由此可以得到與主面O 相鄰的所有從面到主面的權(quán)重。

3.5.4 差值

有了權(quán)重，那么就可以通過(guò)帶權(quán)重差值的方式把速度場(chǎng)、壓力場(chǎng)等在滑移面兩側(cè)實(shí)現(xiàn)數(shù)據(jù)交互。

由主面中變量差值得出從面中變量值:

ΦS= ∑WMj_to_SΦMj

由從面中變量差值得出主面中變量值:

ΦM= ∑WSi_to_MΦSi

式中:ΦMj，ΦSi為某一從面第j 個(gè)主面中Φ 值和某一主面第i 個(gè)從面中的Φ 值。

通過(guò)這種帶權(quán)重的差值，就實(shí)現(xiàn)了網(wǎng)格相對(duì)運(yùn)動(dòng)時(shí)主面和從面之間數(shù)據(jù)的交互。

4 計(jì)算結(jié)果與分析

4.1 推力的時(shí)歷曲線

風(fēng)力機(jī)在非均勻風(fēng)的作用下，葉片在旋轉(zhuǎn)過(guò)程處于不同的位置將經(jīng)歷風(fēng)速的周期性變化，這種變化會(huì)導(dǎo)致葉片在整個(gè)的受風(fēng)面上風(fēng)載荷的不均勻，會(huì)對(duì)風(fēng)力機(jī)的使用年限和運(yùn)行安全產(chǎn)生很大影響。圖7(a)～(d)分別表示一圈內(nèi)不同風(fēng)速均勻風(fēng)和剪切風(fēng)作用下推力隨方位角的變化對(duì)比圖。為方便比較，此處的推力為單個(gè)葉片所受的推力，非均勻風(fēng)下的葉片經(jīng)歷了從高風(fēng)速區(qū)到低風(fēng)速區(qū)再回到高風(fēng)速區(qū)的過(guò)程。從圖中可以明顯觀察到推力也經(jīng)歷了一個(gè)從高到低再到高的過(guò)程，并且推力在風(fēng)速最低處(對(duì)應(yīng)方位角180°)出現(xiàn)了最低值，這與理論分析完全一致。經(jīng)過(guò)計(jì)算發(fā)現(xiàn)一圈內(nèi)的葉片所受推力的平均值變化不超過(guò)5%，但是由于非均勻風(fēng)引起的脈動(dòng)值卻差距很大。例如10 m/s 風(fēng)速下，均勻風(fēng)下推力的脈動(dòng)量占平均推力的7%左右，但是非均勻風(fēng)作用下推力的脈動(dòng)量卻占了平均推力的14.3%，幾乎是均勻風(fēng)的兩倍。對(duì)于25 m/s 風(fēng)速，均勻風(fēng)和非均勻風(fēng)對(duì)推力脈動(dòng)量的影響差異將更大，均勻風(fēng)影響下推力的脈動(dòng)量約是平均推力的4.2%，非均勻風(fēng)影響下推力的脈動(dòng)量約是平均推力的11.1%，約是均勻風(fēng)的三倍。

從上面的分析可以看出，非均勻風(fēng)影響下風(fēng)力機(jī)葉片所受推力的脈動(dòng)值會(huì)比均勻風(fēng)影響下的脈動(dòng)值大很多，這對(duì)于結(jié)構(gòu)的疲勞將產(chǎn)生非常不利的影響。所以在風(fēng)力機(jī)的疲勞設(shè)計(jì)階段，需要根據(jù)具體風(fēng)場(chǎng)的風(fēng)速廓線情況，對(duì)非均勻風(fēng)所引起的疲勞載荷予以考慮。

圖7 四種不同風(fēng)速下剪切風(fēng)與均勻風(fēng)推力時(shí)歷曲線Fig.7 Comparison of time histories of thrust for different uniform and nonuniform winds

4.2 渦的分布

尾渦區(qū)的紊流結(jié)構(gòu)以及葉尖渦的捕獲是風(fēng)力機(jī)氣動(dòng)力學(xué)研究所感興趣的，使用速度梯度張量的二階不變量Q 對(duì)尾流場(chǎng)的渦結(jié)構(gòu)進(jìn)行可視化處理，Q 的定義為:

其中，Ωij代表漩渦的強(qiáng)度，Sij表示剪切應(yīng)變率，分別是速度梯度反對(duì)稱和對(duì)稱分離，可以理解為剪切應(yīng)變率和渦量間的局部平衡量［22］。

2.6 s 輪轂處風(fēng)速為5、10、15 和25 m/s 時(shí)，非均勻風(fēng)影響下風(fēng)力機(jī)尾渦情況如圖8 所示。圖中Q 使用上面定義的速度梯度張量的二階不變量的等值線表達(dá)，并且使用速度染色。風(fēng)輪上部區(qū)域渦的顏色要比下部渦的顏色深，這與非均勻風(fēng)入口設(shè)置是一致的，也驗(yàn)證了文中計(jì)算模型的正確性。從圖中可以清晰看到葉尖渦和葉片根部過(guò)渡區(qū)的渦的脫落，同時(shí)在高風(fēng)速下還能看到葉片處渦的分離。此外，塔架和葉尖渦相互作用也可以明顯看出，而塔架下半部分沒(méi)有明顯的渦生成，可能是下半部分網(wǎng)格變稀疏的緣故造成的。

圖8 不同風(fēng)速非均勻風(fēng)影響下風(fēng)力機(jī)尾渦情況Fig.8 Vortex structure for uniform and nonuniform winds at different speeds

4.3 尾流場(chǎng)分析

尾跡區(qū)分析是風(fēng)力機(jī)空氣動(dòng)力學(xué)分析的重要部分。圖9(a)給出了輪轂處5 m/s 風(fēng)速非均勻風(fēng)作用下不同順風(fēng)向位置中剖面處的風(fēng)速曲線，分別截取了距葉片4、8、12 和16 m 處與塔架平行的直線上風(fēng)速隨高度的變化關(guān)系。從圖中可以清晰的觀察到，距離葉片越近尾流場(chǎng)越紊亂，在距離葉片4 m 縱坐標(biāo)為0 時(shí)，風(fēng)力機(jī)輪轂的地方風(fēng)速有急劇的下降，這主要是由于輪轂的遮蔽作用，同時(shí)也注意到在距離葉片16 m 處風(fēng)速基本恢復(fù)且接近于入口處的非均勻分布。圖9(b)給出了輪轂處10 m/s 風(fēng)速非均勻風(fēng)下不同順風(fēng)向位置中剖面處的風(fēng)速曲線。與圖9(a)一樣，在距離輪轂較近的尾流區(qū)風(fēng)速驟降，隨著據(jù)葉片距離的增大，風(fēng)速逐漸恢復(fù)但是10 m/s 的恢復(fù)比5 m/s 風(fēng)速下要慢很多，可以推出高風(fēng)速下尾流區(qū)的影響更長(zhǎng)，需要更長(zhǎng)的距離進(jìn)行恢復(fù)。在進(jìn)行風(fēng)場(chǎng)風(fēng)力機(jī)布置時(shí)，一定要著重考慮風(fēng)力機(jī)之間的間距以減小上一風(fēng)力機(jī)尾流場(chǎng)對(duì)下一風(fēng)力機(jī)功率的影響。

圖9 輪轂處5 m/s 和10 m/s 橫風(fēng)向尾跡區(qū)風(fēng)速曲線Fig.9 Cross－wake speed distribution at 5 m/s and 10 m/s

為了進(jìn)一步探討非均勻風(fēng)影響下塔架與葉片的相互作用，以10 m/s 風(fēng)速為例，取2.5 s 時(shí)風(fēng)力機(jī)下方葉片r/R=0.15 和r/R=0.95 位置的水平截面進(jìn)行研究，如圖10 所示。圖中背景以速度值染色，渦由Q 等值面表達(dá)并由速度染色，白色圓圈代表塔架，白色翼型代表的是葉片。首先觀察速度場(chǎng)，可以看到塔筒前后的速度值有明顯下降，而塔筒后的速度場(chǎng)有明顯的波浪形變化，這是典型的圓柱繞流的速度場(chǎng)。從圖中還可以看出葉片尖部(z= －4.78 m)的攻角要比葉片根部(z= －1.51 m)的攻角小，葉片尖部泄渦更為明顯。葉片尖部的泄渦與塔架的尾渦相互作用可以從圖10(b)中看出，而葉根處塔架與葉片尾渦相互作用的過(guò)程卻不是特別明顯。在輪轂10 m/s 非均勻風(fēng)速下葉片分離渦中的葉尖渦與塔架尾渦相互作用的過(guò)程更明顯，影響更大。

圖10 10 m/s 非均勻風(fēng)速下風(fēng)力機(jī)下葉片不同截面塔架與葉片相互作用尾渦情況Fig.10 Blade tip vortices and tower vortices interaction for uniform and nonuniform winds at two different sections(cut plane at z= －1.51 m and －4.78 m)

4.4 葉片截面的壓力系數(shù)

通過(guò)將均勻風(fēng)和非均勻風(fēng)作用下葉片各個(gè)截面壓力分布情況進(jìn)行對(duì)比可以進(jìn)一步探討非均勻風(fēng)對(duì)風(fēng)力機(jī)空氣動(dòng)力性能的影響的細(xì)節(jié)。選取了2.5 s 時(shí)(即風(fēng)力機(jī)正好旋轉(zhuǎn)3 圈時(shí))對(duì)四個(gè)不同截面(r/R =0.466，0.633，0.8，0.95)的壓力情況做對(duì)比研究。文中用到的壓力系數(shù)定義如下:

其中，P0為葉片表面附近的壓力值;P∞為無(wú)窮遠(yuǎn)處的壓力值，文中取0;U 代表風(fēng)速;ω 代表風(fēng)輪角速度;r 代表截面距轉(zhuǎn)動(dòng)中心的距離。

從圖11 中能夠清晰的看出在5 m/s 低風(fēng)速時(shí)，非均勻風(fēng)的影響很小，葉片截面的壓力系數(shù)只有細(xì)微差別，同時(shí)隨著截面從根部向葉尖處發(fā)展，這種影響會(huì)越來(lái)越小，這從10 m/s 風(fēng)速下的不同截面的壓力系數(shù)分布能更清晰的看出。

圖11 5 m/s 和10 m/s 均勻風(fēng)和非均勻風(fēng)葉片各截面壓力系數(shù)分布曲線Fig.11 Pressure coefficient distribution for uniform and nonuniform winds at four sections at 5，10 m/s

5 結(jié) 語(yǔ)

在考慮非均勻風(fēng)影響的情況下，對(duì)風(fēng)力機(jī)整體結(jié)構(gòu)的三維非定常氣動(dòng)問(wèn)題進(jìn)行了數(shù)值模擬，得出以下結(jié)論:

1)在非均勻風(fēng)影響下，風(fēng)力機(jī)所受推力的脈動(dòng)值會(huì)比均勻風(fēng)影響下的脈動(dòng)值大很多，這對(duì)結(jié)構(gòu)的疲勞將產(chǎn)生非常不利的影響。所以在風(fēng)力機(jī)的疲勞設(shè)計(jì)階段，需要根據(jù)具體風(fēng)場(chǎng)的風(fēng)速廓線情況，對(duì)非均勻風(fēng)所引起的疲勞載荷予以考慮。

2)非均勻風(fēng)的存在會(huì)引起尾渦結(jié)構(gòu)的非對(duì)稱，同時(shí)對(duì)尾流場(chǎng)的分析中可以發(fā)現(xiàn)尾流區(qū)風(fēng)速分布非常紊亂，且高風(fēng)速下風(fēng)速的恢復(fù)更慢。

3)非均勻風(fēng)影響下的葉片截面壓力系數(shù)在低風(fēng)速時(shí)與均勻風(fēng)相比變化不大，但是高風(fēng)速時(shí)有很大差距，且隨著截面從葉片根部到葉片尖部發(fā)展，這種影響會(huì)逐漸變小。

［1］S?RENSEN N N，MICHELSEN J，SCHRECK S.Navier－Stokes predictions of the NREL phase VI rotor in the NASA Ames 80 ft× 120 ft wind tunnel［J］.Wind Energy，2002，5(2－3):151－169.

［2］李宇紅，張慶麟.風(fēng)力機(jī)葉片三維流動(dòng)特性與氣動(dòng)性能的數(shù)值分析［J］.太陽(yáng)能學(xué)報(bào)，2008，29(9):1172－1176.(LI Yuhong，ZHANG Qinglin.Numerical simulation of flow field and aerodynamic performance of a wind turbine blade［J］.Acta Energiae Solaris Sinica，2008，29(9):1172－1176.(in Chinese))

［3］張果宇，蔣勁，劉長(zhǎng)陸.風(fēng)力發(fā)電機(jī)整機(jī)氣動(dòng)性能數(shù)值模擬計(jì)算與仿真研究［J］.華東電力，2009，37(3):449－452.(ZHANG Guoyu，JIANG Jin，LIU Changlu.Numerical simulation of aerodynamic performance for wind turbines［J］.East China Electric Power，2009，37(3):449－452.(in Chinese))

［4］任年鑫，歐進(jìn)萍.大型海上風(fēng)力機(jī)尾跡區(qū)域風(fēng)場(chǎng)分析［J］.計(jì)算力學(xué)學(xué)報(bào)，2012，29(3):327－331.(REN Nianxin，OU Jinping.Numerical analysis for the wake zone of large offshore wind turbine［J］.Chinese Journal of Computational Mechanics，2012，29(3):327－331.(in Chinese))

［5］SEZER－UZOL N，UZOL O.Effect of steady and transient wind shear on the wake structure and performance of a horizontal axis wind turbine rotor［J］.Wind Energy，2013，16(1):1－17.

［6］劉磊，石可重，楊科，等.風(fēng)切變對(duì)風(fēng)力機(jī)氣動(dòng)載荷的影響［J］.工程熱物理學(xué)報(bào)，2010，31(10):1667－1670.(LIU Lei，SHI Kezhong，YANG Ke，et al.Effect of wind shear on the aerodynamic load of wind turbine［J］.Journal of Engineering Thermophysics，2010，31(10):1667－1670.(in Chinese))

［7］廖明夫，徐可，吳斌，等.風(fēng)切變對(duì)風(fēng)力機(jī)功率的影響［J］.沈陽(yáng)工業(yè)大學(xué)學(xué)報(bào)，2008，30(2):163－167.(LIAO Mingfu，XU Ke，WU Bin，et al.Effect of wind shear on wind turbine power［J］.Journal of Shengyang University of Technology，2008，30(2):163－167.(in Chinese))

［8］JASAK H，JEMCOV A，TUKOVIC Z.Openfoam:A c + + library for complex physics simulations［C］//Proceedings of International Workshop on Coupled Methods in Numerical Dynamics.2007.

［9］查晶晶，萬(wàn)德成.用OpenFOAM 實(shí)現(xiàn)數(shù)值水池造波和消波［J］.海洋工程，2011，29(3):1－12.(CHA Jingjing ，WAN Decheng.Numerical wave generation and absorption based on OpenFOAM［J］.The Ocean Engineering，2011，29(3):1－12.(in Chinese))

［10］CAO H J，WAN D C.Development of multidirectional nonlinear numerical wave tank by naoe－FOAM－SJTU solver［J］.International Journal of Ocean System Engineering，2014，4(1):52－59.

［11］SHEN Z R，WAN D C.RANS computations of added resistance and motions of ship in head waves［J］.International Journal of Offshore and Polar Engineering，2013，23(4):263－271.

［12］周胡，王強(qiáng)，萬(wàn)德成.風(fēng)機(jī)葉片三維繞流場(chǎng)數(shù)值模擬［C］//第十一屆全國(guó)水動(dòng)力學(xué)學(xué)術(shù)會(huì)議暨第二十四屆全國(guó)水動(dòng)力學(xué)研討會(huì)并周培源教授誕辰110 周年紀(jì)念大會(huì).2012:627－636.(ZHOU Hu，WANG Qiang，WAN Decheng.Numerical Simulation of the 3D viscous flow field over wind turbine Blade［C］//Proceedings of 11th National Conference of Hydrodynamics.2012:627－636.(in Chinese))

［13］WANG Q，ZHOU H，WAN D C.Numerical simulation of wind turbine blade－tower interaction［J］.Journal of Marine Science and Application，2012，11(3):321－327.

［14］周胡，萬(wàn)德成.下風(fēng)向風(fēng)力機(jī)塔影效應(yīng)的非定常數(shù)值模擬［C］//第二十五屆全國(guó)水動(dòng)力學(xué)研討會(huì)暨第十二屆全國(guó)水動(dòng)力學(xué)學(xué)術(shù)會(huì)議.2013:277－283.(ZHOU Hu，WAN Decheng.Unsteady numerical simulation of tower shadow of downwind wind turbine［C］//Proceedings of 12th National Conference of Hydrodynamics.2013:277－283.(in Chinese))

［15］FINGERSH L J，SIMMS D，HAND M，et al.Wind tunnel testing of NREL’s unsteady aerodynamics experiment［J］.AIAA paper，2001:35－46.

［16］HAND M M，SIMMS D，F(xiàn)INGERSH L，et al.Unsteady aerodynamics experiment phase v:test configuration and available data campaigns［R］.National Renewable Energy Laboratory，2001.

［17］MENTER F R.Two －equation eddy －viscosity turbulence models for engineering applications［J］.AIAA Journal，2012，32(8):76－82.

［18］JASAK H.Error analysis and estimation for the finite volume method with applications to fluid flows［D］.London:University of London，1996.

［19］賀德馨.風(fēng)工程與工業(yè)空氣動(dòng)力學(xué)［M］.北京:國(guó)防工業(yè)出版社，2006.(HE Dexin.Wind engineering and industrial aerodynamics［M］.Beijing:National Defense Industry Press，2006.(in Chinese))

［20］GB.建筑結(jié)構(gòu)荷載規(guī)范［S］.2002.(GB.Load code for the design of building structures［S］.2002.(in Chinese))

［21］JASAK H，BEAUDOIN M.OpenFOAM TURBO TOOLS:from general purpose CFD to turbomachinery simulations［C］//Proceedings of ASME－JSME－KSME Joint Fluids Engineer－ing Conference (AJK2011－FED).2011.

［22］JEONG J，HUSSAIN F.On the identification of a vortex［J］.Journal of Fluid Mechanics，1995，285(69):69－94.