非均勻風影響下風力機三維氣動粘性流場的數(shù)值模擬

周 胡，趙文超，萬德成

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

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

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

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

1 OpenFOAM 簡介

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

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

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

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

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

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

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

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

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

3 數(shù)值模擬

3.1 控制方程

在慣性參考系中，風力機風輪以設定的角速度繞固定軸轉動，流場為非定常。由于風力機外流場的風速較低，繞流場可以看成不可壓。其時均控制方程可以表述為:

3.2 湍流模型

為了計算雷諾應力項，需要引入k－ω SST 模型［17］使控制方程封閉。其中k 和ω 分別代表湍動能和渦量脈動強度。兩者的輸運方程為:

3.3 離散方法

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

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

3.4 邊界條件

3.4.1 入口邊界條件

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

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

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

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

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

其中，inlet 是入口邊界的名字，n 是來流風速的方向，meanValue 是輪轂處的風速，value 為入口速度初始值。最終計算時入口的速度分布如圖4 所示，與預定義的速度廓線一致，說明入口邊界條件編寫正確。

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

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

3.4.2 出口邊界條件出口處的速度邊界條件為零梯度，壓力為固定值0，湍動能k 和渦量脈動強度ω 設為固定值，湍動能k 和渦量脈動強度ω 設為零梯度，渦粘度通過計算求得。

3.4.3 壁面條件

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

3.4.4 風力機模型邊界條件

風力機的塔架和不動的輪轂的速度設為固定值0，葉片和轉動的輪轂設為movingWallVelocity;整個風力機的壓力邊界條件設為零梯度;湍動能k、渦量脈動強度ω 和渦粘度采用壁面函數(shù)。

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

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

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

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

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

3.5.2 尋找相鄰面

當AMI 交界面生成后，就可以開始計算了。計算的第一步便是尋找相鄰面。尋找相鄰面的過程是已知主面尋找相鄰從面的過程。

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

3.5.3 計算權重

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

主面到從面的權重:

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

從面到主面的權重:

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

例如對于圖6，想要計算從面b 到主面O 的權重，那么可以通過以上公式計算得到權重。由此可以得到與主面O 相鄰的所有從面到主面的權重。

3.5.4 差值

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

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

ΦS= ∑WMj_to_SΦMj

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

ΦM= ∑WSi_to_MΦSi

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

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

4 計算結果與分析

4.1 推力的時歷曲線

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

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

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

4.2 渦的分布

尾渦區(qū)的紊流結構以及葉尖渦的捕獲是風力機氣動力學研究所感興趣的，使用速度梯度張量的二階不變量Q 對尾流場的渦結構進行可視化處理，Q 的定義為:

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

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

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

4.3 尾流場分析

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

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

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

圖10 10 m/s 非均勻風速下風力機下葉片不同截面塔架與葉片相互作用尾渦情況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ù)

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

其中，P0為葉片表面附近的壓力值;P∞為無窮遠處的壓力值，文中取0;U 代表風速;ω 代表風輪角速度;r 代表截面距轉動中心的距離。

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

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

5 結 語

在考慮非均勻風影響的情況下，對風力機整體結構的三維非定常氣動問題進行了數(shù)值模擬，得出以下結論:

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

2)非均勻風的存在會引起尾渦結構的非對稱，同時對尾流場的分析中可以發(fā)現(xiàn)尾流區(qū)風速分布非常紊亂，且高風速下風速的恢復更慢。

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

［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］李宇紅，張慶麟.風力機葉片三維流動特性與氣動性能的數(shù)值分析［J］.太陽能學報，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］張果宇，蔣勁，劉長陸.風力發(fā)電機整機氣動性能數(shù)值模擬計算與仿真研究［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］任年鑫，歐進萍.大型海上風力機尾跡區(qū)域風場分析［J］.計算力學學報，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］劉磊，石可重，楊科，等.風切變對風力機氣動載荷的影響［J］.工程熱物理學報，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］廖明夫，徐可，吳斌，等.風切變對風力機功率的影響［J］.沈陽工業(yè)大學學報，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］查晶晶，萬德成.用OpenFOAM 實現(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］周胡，王強，萬德成.風機葉片三維繞流場數(shù)值模擬［C］//第十一屆全國水動力學學術會議暨第二十四屆全國水動力學研討會并周培源教授誕辰110 周年紀念大會.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］周胡，萬德成.下風向風力機塔影效應的非定常數(shù)值模擬［C］//第二十五屆全國水動力學研討會暨第十二屆全國水動力學學術會議.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］賀德馨.風工程與工業(yè)空氣動力學［M］.北京:國防工業(yè)出版社，2006.(HE Dexin.Wind engineering and industrial aerodynamics［M］.Beijing:National Defense Industry Press，2006.(in Chinese))

［20］GB.建筑結構荷載規(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.