基于自由渦方法的風(fēng)輪偏航氣動(dòng)特性研究

仇永興，康 順

（1.華北電力大學(xué)電站設(shè)備狀態(tài)檢測與控制教育部重點(diǎn)實(shí)驗(yàn)室，北京 102206；2.西安現(xiàn)代控制技術(shù)研究所，陜西西安 710065）

基于自由渦方法的風(fēng)輪偏航氣動(dòng)特性研究

仇永興1，＊，康 順1，2

（1.華北電力大學(xué)電站設(shè)備狀態(tài)檢測與控制教育部重點(diǎn)實(shí)驗(yàn)室，北京 102206；2.西安現(xiàn)代控制技術(shù)研究所，陜西西安 710065）

采用非線性升力線方法和時(shí)間精確自由渦方法，對風(fēng)切變條件下風(fēng)輪的偏航氣動(dòng)特性進(jìn)行了非定常數(shù)值模擬。為了改善非線性升力線方法的計(jì)算收斂性，提出了基于Newton-Raphson方法的環(huán)量方程迭代算法。建立了一種包含粘性渦核并由自由渦面和葉尖渦線構(gòu)成的改進(jìn)尾跡模型，可以提高尾跡的計(jì)算精度和效率。以NREL Phase VI風(fēng)力機(jī)和Tjreborg風(fēng)力機(jī)為算例，分別對均勻風(fēng)和切變風(fēng)條件下、風(fēng)輪處于不同偏航角時(shí)的氣動(dòng)力進(jìn)行計(jì)算，并與測量值進(jìn)行了比較。分析了風(fēng)切變和偏航風(fēng)對風(fēng)輪氣動(dòng)特性產(chǎn)生影響的主要因素及風(fēng)輪氣動(dòng)特性的變化規(guī)律。結(jié)果表明：計(jì)算方法可以有效模擬風(fēng)切變條件下風(fēng)輪偏航氣動(dòng)特性；在風(fēng)切變條件下，偏航會增加葉片扭矩的波動(dòng)，當(dāng)風(fēng)輪處于正偏航角時(shí)，葉片揮舞力矩的波動(dòng)減弱，處于負(fù)偏航角時(shí)，揮舞力矩波動(dòng)增強(qiáng)。

風(fēng)輪；風(fēng)切變；偏航風(fēng)；非定常模擬；自由渦

0 引 言

風(fēng)力機(jī)運(yùn)行在風(fēng)速、風(fēng)向不斷變化的風(fēng)場中。由于大氣邊界層的作用，來流風(fēng)存在垂直風(fēng)切變。在處于偏航和風(fēng)切變條件下，風(fēng)輪葉片會受到較復(fù)雜的非定常氣動(dòng)載荷，影響葉片疲勞壽命與安全。對風(fēng)輪偏航氣動(dòng)特性的研究，?；诰鶆蝻L(fēng)假設(shè)［1-3］。因此，研究風(fēng)切變條件下風(fēng)輪偏航氣動(dòng)特性具有十分重要且實(shí)際的工程意義。

基于URANS方程的CFD（Computational Fluid Dynamics）方法用于研究復(fù)雜來流條件下風(fēng)輪的氣動(dòng)載荷模擬，仍然存在計(jì)算成本過高的問題。當(dāng)前工業(yè)領(lǐng)域應(yīng)用最為廣泛的方法仍然是基于葉素動(dòng)量理論的相關(guān)方法。這些方法的主要缺點(diǎn)是假設(shè)條件多，大量依靠經(jīng)驗(yàn)?zāi)Ｐ?。特別是使用動(dòng)量尾跡模型過度簡化了風(fēng)輪尾跡運(yùn)動(dòng)，不能準(zhǔn)確模擬復(fù)雜來流條件下尾跡的運(yùn)動(dòng)，從而造成誘導(dǎo)速度計(jì)算不準(zhǔn)確，嚴(yán)重影響計(jì)算精度［4］。

基于勢流假設(shè)的自由渦方法，使用環(huán)量沿徑向連續(xù)變化的附著渦線模擬葉片的氣動(dòng)效應(yīng)；用沿著葉片徑向分布的尾渦線和脫體渦線模擬葉片的尾跡.隨著風(fēng)輪的轉(zhuǎn)動(dòng)，形成螺旋形的葉片尾跡。假設(shè)尾跡渦線沿著流線方向自由運(yùn)動(dòng)，其空間位置是流場解的一部分，由來流條件和渦線的相互擾動(dòng)作用共同決定。自由渦方法已經(jīng)在直升機(jī)旋翼數(shù)值模擬領(lǐng)域取得豐富成果［5-6］，并在風(fēng)力機(jī)數(shù)值模擬領(lǐng)域快速發(fā)展［7-10］。其優(yōu)勢在于假設(shè)條件少，能夠模擬葉片和尾跡的運(yùn)動(dòng)過程［11］，并解出流場和葉片的非定常氣動(dòng)載荷。相較于葉素動(dòng)量理論，自由渦方法更接近物理實(shí)際，同時(shí)在計(jì)算量和計(jì)算時(shí)間上又顯著少于CFD方法。

自由渦方法包括葉片氣動(dòng)計(jì)算和渦尾跡計(jì)算兩部分。葉片氣動(dòng)計(jì)算常采用非線性升力線方法。由于使用翼型氣動(dòng)實(shí)驗(yàn)數(shù)據(jù)計(jì)算附著渦環(huán)量，計(jì)入了葉片邊界層及分離流動(dòng)的粘性效應(yīng)，因此計(jì)算更加準(zhǔn)確。由于該方法常使用簡單迭代算法求解環(huán)量的非線性方程［7］，導(dǎo)致環(huán)量迭代不易收斂。

渦尾跡計(jì)算是通過求解渦線運(yùn)動(dòng)的控制方程，得到渦尾跡的瞬時(shí)空間位置。常用的方法有二階龍格庫塔方法［7，10］（2nd-order Runge-Kutta scheme，RK2）和預(yù)測-校正二階向后差分（Predictor-Corrector with 2nd-ordor Backward difference scheme，PC2B）的時(shí)間精確解法。Leishman等的研究表明，PC2B方法在計(jì)算精度和數(shù)值穩(wěn)定性上優(yōu)于其他方法［12］。但是該方法的計(jì)算量大、效率低。因此，需要簡化渦尾跡模型，降低計(jì)算量。常將風(fēng)力機(jī)尾跡簡化為由葉尖渦線構(gòu)成的尾跡［13］，使尾跡模型中不包含脫體渦，降低了尾跡的計(jì)算精度。

為了解決上述問題，改進(jìn)了自由渦的求解方法。為了改善采用非線性升力線方法中環(huán)量方程迭代收斂性，提出了一種基于Newton-Raphson方法的非線性環(huán)量方程迭代算法。為了提高時(shí)間精確自由渦方法的計(jì)算精度和效率，提出了一種由自由渦面（包含尾渦線和脫體渦線）和葉尖渦線構(gòu)成的簡化渦尾跡模型。以NREL Phase VI兩葉片失速型風(fēng)力機(jī)［14-15］和Tjreborg三葉片兆瓦級風(fēng)力機(jī)［16］為例，分別對均勻來流和風(fēng)切變條件下，風(fēng)輪處于不同風(fēng)速和偏航角時(shí)的氣動(dòng)力進(jìn)行了計(jì)算，并與測量值進(jìn)行比較。分析了風(fēng)切變和偏航對風(fēng)輪氣動(dòng)特性產(chǎn)生影響的主要因素。研究了風(fēng)切變條件下，風(fēng)輪不同偏航角時(shí)，氣動(dòng)特性的變化規(guī)律。結(jié)果表明，所建立的改進(jìn)自由渦方法具有較高的計(jì)算精度和計(jì)算效率，能夠有效模擬風(fēng)切變條件下風(fēng)輪偏航氣動(dòng)特性。該方法對于推動(dòng)風(fēng)力機(jī)的工程計(jì)算及優(yōu)化設(shè)計(jì)都具有重要而且實(shí)際的意義。

1 研究方法與數(shù)值模型

1.1 非線性升力線方法

葉片的氣動(dòng)效應(yīng)使用非線性升力線方法模擬。將葉片簡化為一根環(huán)量階躍變化的附著渦線（升力線）位于葉片翼型的1/4弦長處。由于氣動(dòng)力沿葉片徑向的變化，環(huán)量也相應(yīng)地變化。依據(jù)Helmholtz定律和Kelvin定律，升力線的環(huán)量沿葉片徑向變化會產(chǎn)生尾渦線，隨時(shí)間發(fā)生變化會產(chǎn)生脫體渦線。尾渦和脫體渦從葉片尾緣流向下游。因此，在葉片上附著渦線、尾渦線和脫體渦線構(gòu)成渦環(huán)，如圖1所示。

圖1 葉片微段的局部坐標(biāo)系及渦系構(gòu)成示意圖Fig.1 Schematic of blade element coordinate and vortex system composition

依據(jù)Joukowski升力定理，葉片微段上附著渦線段產(chǎn)生的升力為：

也可由翼型升力系數(shù)計(jì)算葉片微元段的升力，即：

其中，ГB是葉片微段的附著渦環(huán)量；u是葉片微段的合速度，包括來流風(fēng)速、葉片旋轉(zhuǎn)線速度、葉片附著渦段中點(diǎn)的擾動(dòng)速度。擾動(dòng)速度通過Biot-Savart定律解出，方法參見文獻(xiàn)［10］；c是葉片微段的當(dāng)?shù)叵议L；CL是翼型升力系數(shù)，可由實(shí)驗(yàn)或CFD數(shù)據(jù)獲得，隱含了邊界層及邊界層分離的影響。聯(lián)立式（1）和式（2）可構(gòu)建求解附著渦環(huán)量ГB的方程，常采用線性迭代法求解［7］。

由于環(huán)量方程是非線性方程，使用線性迭代法往往造成解的震蕩甚至不收斂。為了提高解的精度和數(shù)值穩(wěn)定性，構(gòu)建了基于Newton-Raphson方法的非線性方程迭代算法。將環(huán)量方程寫為：

其中，f是一個(gè)趨近于0的小量；n表示附著渦微段序號。環(huán)量計(jì)算過程為：

（i）給定葉片和尾跡渦的初始環(huán)量和初始尾跡；

（ii）計(jì)算附著渦段中點(diǎn)的合速度，解出葉片攻角徑向分布；

（iii）通過葉片當(dāng)?shù)匾硇偷臍鈩?dòng)特性曲線，得到升力系數(shù)；

（iv）解方程組（5），得到葉片附著渦環(huán)量的增量；

判斷新舊環(huán)量值的最大相對誤差小于0.001為收斂，若未收斂則用下式更新環(huán)量值，并返回步（ii）繼續(xù)計(jì)算：

上式中，松弛因子ω推薦取0.5。方程（5）的Jacobi矩陣中的偏微分運(yùn)算，使用數(shù)值微分求解。

解出附著渦環(huán)量就得到了整個(gè)渦系的環(huán)量值，同時(shí)獲得葉片合速度、攻角以及氣動(dòng)載荷等沿徑向的分布。此外，以上計(jì)算均假設(shè)葉片各微元段的繞流是二維流動(dòng)，為了計(jì)入三維效應(yīng)引入了Du-Selig三維失速延遲模型，見文獻(xiàn)［17］。

1.2 自由渦計(jì)算

1.2.1 渦尾跡模型

使用時(shí)間精確算法求解尾跡渦線運(yùn)動(dòng)，計(jì)算精度高但是效率低。若完全滿足Helmholtz定律和Kelvin定律，則需要計(jì)算尾跡中所有尾渦線和脫體渦線的運(yùn)動(dòng)，計(jì)算量過大，顯著降低了自由渦方法的效率優(yōu)勢。因此，有必要簡化尾跡模型。一般的方法是使用葉尖尾渦線代替整個(gè)自由渦，忽略了脫體渦的運(yùn)動(dòng)，降低了尾跡模型的精度。

為了有效模擬脫體渦的運(yùn)動(dòng)和擾動(dòng)作用，并降低尾跡計(jì)算量，建立了一種由風(fēng)輪附近的自由運(yùn)動(dòng)的渦面（包含尾渦線和脫體渦線）和較遠(yuǎn)下游的葉尖渦線構(gòu)成的渦尾跡模型。模型假設(shè)：

（i）葉尖渦是由尾渦構(gòu)成的集中渦。渦面運(yùn)動(dòng)到風(fēng)輪下游一定距離后，全部卷入葉尖渦。葉尖渦渦線從渦面位于葉尖區(qū)域的卷曲中心引出，如圖2所示；

圖2 風(fēng)輪單根葉片的渦尾跡Fig.2 Vortex wake of one blade

（ii）尾跡中的脫體渦運(yùn)動(dòng)到下游一定距離后，對葉片流場的擾動(dòng)作用忽略不計(jì)。所以，脫體渦只包含在自由渦面模型中；

（iii）每根渦線在運(yùn)動(dòng)過程中，環(huán)量值保持不變，但是渦核半徑會隨著運(yùn)動(dòng)時(shí)間和渦線拉伸率發(fā)生變化。

前兩條假設(shè)，是基于流動(dòng)顯示實(shí)驗(yàn)的觀察結(jié)果［18-19］，最后一條假設(shè)是基于Helmholtz渦守恒定律。為了保證流場中尾渦環(huán)量的連續(xù)性，渦面的尾渦環(huán)量需要傳遞給葉尖渦線。具體方法是，在兩者的連接邊界處，將渦面中與葉尖渦同向的尾渦環(huán)量沿著尾跡徑向積分，該積分值為葉尖渦線的環(huán)量值。

為了模擬尾跡運(yùn)動(dòng)中的粘性效應(yīng)，同時(shí)避免尾跡速度場的奇異性，在渦尾跡模型中引入Lamb-Oseen粘性渦核模型［20］：

其中：

式中，a＝1.25643；渦線拉伸率ε是相鄰兩個(gè)時(shí)刻渦線長度的相對變化率；δV為湍流粘性系數(shù)，文獻(xiàn)［10］指出對于小型風(fēng)力機(jī)該參數(shù)一般取100，大型風(fēng)力機(jī)取1000較好；Sc為渦線運(yùn)動(dòng)的初始時(shí)間，一般取0.5～1保證初始渦核半徑rc約等于葉片平均厚度的一半。引入渦核模型后的擾動(dòng)速度計(jì)算參見文獻(xiàn)［10］。

1.2.2 渦尾跡控制方程及求解

渦尾跡由自由渦面和葉尖渦線構(gòu)成。而構(gòu)成自由渦面和葉尖渦線的是離散的渦線段。渦線段的運(yùn)動(dòng)控制方程為：

采用預(yù)測-校正二階向后差分法（PC2B）求解控制方程（7）。式中Ψ為葉片旋轉(zhuǎn)角度；ζ為風(fēng)輪相對坐標(biāo)系下尾跡節(jié)點(diǎn)的旋轉(zhuǎn)幅角；Ω為風(fēng)輪旋轉(zhuǎn)角速度；uw為尾跡節(jié)點(diǎn)處流速，是風(fēng)速和渦擾動(dòng)速度的合速度。方程（7）離散后，得到方程（8）有：

式中，j為Ψ增加方向序號；k為ζ增加方向序號。方程（8）在求解過程采用了偽隱式解法，具體解法參見文獻(xiàn)［21］。

2 算例模型與方案介紹

為了驗(yàn)證計(jì)算方法的有效性，并且研究風(fēng)切變條件下風(fēng)輪偏航氣動(dòng)特性，選取以下兩種風(fēng)力機(jī)作為算例模型：

NREL Phase VI風(fēng)力機(jī)：這是一個(gè)失速型兩葉片風(fēng)力機(jī)，風(fēng)輪直徑10.046m，葉片截面采用S809翼型。該風(fēng)力機(jī)由美國可再生能源實(shí)驗(yàn)室設(shè)計(jì)，在NASA AMES研究中心的大型風(fēng)洞中進(jìn)行實(shí)驗(yàn)，由于實(shí)驗(yàn)全面且數(shù)據(jù)精確，一直作為驗(yàn)證風(fēng)力機(jī)氣動(dòng)計(jì)算方法的重要參考［15］。

研究方案具體如下：

（1）方法有效性驗(yàn)證。NREL Phase VI風(fēng)力機(jī)，風(fēng)速5～12m/s，葉片槳矩角3°，無偏航，計(jì)算風(fēng)力機(jī)功率和推力；風(fēng)速7m/s，葉片槳矩角3°，給定風(fēng)輪偏航角為10°、30°、45°和60°，計(jì)算葉片旋轉(zhuǎn)一周的平均氣動(dòng)力。Tjreborg風(fēng)力機(jī)，風(fēng)速5～13m/s，葉片槳矩角0.4°，無偏航，計(jì)算風(fēng)力機(jī)功率。

式中，uhub和Hhub分別為風(fēng)輪輪轂位置的風(fēng)速和輪轂高度。

3 結(jié)果與分析

3.1 方法有效性驗(yàn)證

圖3 功率和推力隨風(fēng)速的變化曲線對比，NREL Phase VIFig.3 Comparison of the computed power and the thrust with respect to the wind speed，NREL Phase VI

圖3給出了均勻風(fēng)條件下，NREL Phase VI風(fēng)力機(jī)葉片槳矩角3°、無偏航時(shí)功率和推力計(jì)算值與實(shí)驗(yàn)值的比較。圖4給出了均勻風(fēng)條件下，Tjreborg風(fēng)力機(jī)，葉片槳矩角0.4°，無偏航時(shí)功率計(jì)算值與實(shí)驗(yàn)值［14］的比較。

圖3中可見，風(fēng)速小于10m/s時(shí)，功率和推力的計(jì)算值與實(shí)驗(yàn)值較好吻合。當(dāng)風(fēng)速大于10m/s后，功率計(jì)算值與實(shí)驗(yàn)值誤差變大，但推力仍較好吻合。主要原因是葉片已經(jīng)處于失速狀態(tài)，吸力面存在較強(qiáng)的三維分離流動(dòng)［17］。計(jì)算方法受限于二維流動(dòng)假設(shè)和三維失速延遲模型（Du-Selig三維失速延遲模型［17］）的計(jì)算能力，模擬葉片大分離工況存在較大誤差。圖4中可見，計(jì)算值與實(shí)驗(yàn)值很好吻合。

實(shí)驗(yàn)值的波動(dòng)可能與實(shí)際運(yùn)行中風(fēng)入流存在風(fēng)切變有關(guān)，下面將詳細(xì)研究。

圖4 功率隨風(fēng)速的變化曲線對比，TjreborgFig.4 Variation of the computed power with respect to the wind speed，Tjreborg

圖5給出了均勻風(fēng)條件下，NREL Phase VI風(fēng)力機(jī)，葉片槳矩角3°，風(fēng)速7m/s時(shí)，風(fēng)輪處于不同偏航角時(shí)葉片法向和切向載荷系數(shù)沿徑向的分布。法向和切向載荷系數(shù)的定義見文獻(xiàn)［14］?？梢?，當(dāng)偏航角為10°、30°和45°時(shí)，計(jì)算結(jié)果與實(shí)驗(yàn)數(shù)據(jù)很好吻合。當(dāng)風(fēng)輪處于60°偏航角時(shí)，葉片存在較大分離流動(dòng)和徑向二次流動(dòng)，影響了計(jì)算的準(zhǔn)確性，但是計(jì)算值在趨勢上仍然與實(shí)驗(yàn)值保持一致，且切向力系數(shù)值與實(shí)驗(yàn)值吻合良好。

以上計(jì)算表明，計(jì)算方法模擬風(fēng)輪氣動(dòng)特性是準(zhǔn)確有效的。

3.2 風(fēng)切變條件下風(fēng)輪氣動(dòng)特性分析

為了便于分析，圖7中定義了葉片方位角和偏航角。

圖5 給定偏航角的葉片氣動(dòng)力系數(shù)沿徑向分布對比，NREL phase VIFig.5 Comparison of the radial distributions of the computed aerodynamic force coefficients under the yaw conditions，NREL phase VI

圖6 風(fēng)速沿高度方向的分布對比Fig.6 Comparison of the wind speed distributions in the altitudinal direction

圖8中，葉片I運(yùn)動(dòng)到方位角約90°時(shí)，其轉(zhuǎn)矩為最大值；運(yùn)動(dòng)到約270°時(shí)為最小值。這是因?yàn)轱L(fēng)切變條件下，風(fēng)輪上半周風(fēng)速較高，下半周風(fēng)速較低。當(dāng)風(fēng)輪無偏航時(shí)，旋轉(zhuǎn)平面最高點(diǎn)和最低點(diǎn)附近分別存在葉片氣動(dòng)載荷的極大、極小值。

圖7 葉片方位角偏航角定義示意圖Fig.7 Schematic of definitions of the blade azimuth angle and yaw angle

圖8 風(fēng)切變時(shí)風(fēng)輪及葉片轉(zhuǎn)矩隨方位角的變化Fig.8 Variation of the computed shaft torque of the rotor and blades with respect to the azimuth angle under shear wind condition

3.3 偏航條件下風(fēng)輪氣動(dòng)特性分析

定義葉片氣動(dòng)力的兩個(gè)分量作為分析對象。沿著葉片徑向分布的氣動(dòng)力，在風(fēng)輪旋轉(zhuǎn)方向的分量是切向力，表示為Ftangent；在風(fēng)輪軸向的分量是推力，表示為Taxial。沿著葉片徑向，切向力和推力分別向風(fēng)輪旋轉(zhuǎn)中心取矩，其合力矩分別是葉片的扭矩和揮舞力矩。定義葉片切向力系數(shù)為：葉片推力系數(shù)為：

式中，ρ為空氣密度；uhub為風(fēng)輪輪轂位置風(fēng)速，R為風(fēng)輪掃風(fēng)面積。由式（9）、式（10）可見，葉片上切向力系數(shù)和推力系數(shù)的分布表示了葉片氣動(dòng)力在旋轉(zhuǎn)方向和推力方向的分布。

圖9 不同偏航角時(shí)葉片切向力系數(shù)在旋轉(zhuǎn)平面內(nèi)的分布云圖Fig.9 Distributions of the computed tangent force coefficient of the blade on the rotation plane under yaw conditions

圖10 不同偏航角時(shí)葉片推力系數(shù)在旋轉(zhuǎn)平面內(nèi)的分布云圖Fig.10 Distributions of the computed axial thrust coefficient of the blade on the rotation plane under yaw conditions

圖9和圖10分別給出了均勻風(fēng)和切變風(fēng)條件下，輪轂風(fēng)速10.9m/s、Tjreborg風(fēng)力機(jī)處于不同偏航角時(shí)，單個(gè)葉片的切向力系數(shù)和推力系數(shù)在旋轉(zhuǎn)平面內(nèi)的分布云圖。圖中可見，偏航角相反時(shí)，均勻風(fēng)條件下，氣動(dòng)力系數(shù)的分布呈180°中心對稱；切變風(fēng)條件下分布不具有對稱性。隨著葉片旋轉(zhuǎn)，切向力系數(shù)的波動(dòng)在大部分葉片半徑范圍內(nèi)都存在，而推力系數(shù)的波動(dòng)主要集中在葉片60%半徑之外。

圖9（b）可見，切變風(fēng)條件下，正偏航時(shí)，葉片切向力系數(shù)的分布與圖9（a）相比更加不均勻。其中，風(fēng)輪上半周切向力系數(shù)極大值和極值區(qū)域都有所增大，風(fēng)輪下半周切向力系數(shù)極小值更小，極小值區(qū)域更大；負(fù)偏航時(shí)，旋轉(zhuǎn)平面內(nèi)切向力系數(shù)分布相較均勻風(fēng)時(shí)，極大值區(qū)域減小。當(dāng)葉片運(yùn)行在方位角約-120°～90°區(qū)間內(nèi)，切向力系數(shù)隨相位角的變化比較平緩，沿徑向的分布比較均勻。

圖10（b）可見，切變風(fēng)條件下，正偏航時(shí)，旋轉(zhuǎn)平面內(nèi)推力系數(shù)的分布與圖10（a）相比波動(dòng)減小，推力系數(shù)的極值區(qū)域和變化幅度都有所減?。回?fù)偏航時(shí)，推力系數(shù)的波動(dòng)增大。當(dāng)葉片運(yùn)行在方位角約45°～225°區(qū)間時(shí)，葉片半徑60%以外的區(qū)域，推力系數(shù)保持在極大值。

圖11給出了均勻風(fēng)和切變風(fēng)條件下，不同偏航角時(shí)，單個(gè)葉片的扭矩和揮舞力矩隨方位角的變化對比。圖中可見，與均勻風(fēng)條件相比，切變風(fēng)條件下，風(fēng)輪正偏航時(shí)，葉片扭矩隨方位角的波動(dòng)明顯增大。其中，扭矩極大值由方位角0°附近移動(dòng)至60°附近，極小值由方位角200°附近移動(dòng)至240°附近。揮舞力矩的波動(dòng)明顯減??；負(fù)偏航時(shí)，扭矩的波動(dòng)幅值沒有明顯變化，但是極大值所在的方位角大約從210°變化至150°附近。揮舞力矩的波幅明顯增大，極值所在方位角沒有明顯變化。

對圖11中給出的葉片扭矩和揮舞力矩的波動(dòng)規(guī)律進(jìn)行定量比較：相對于均勻風(fēng)條件下葉片扭矩和揮舞力矩的波幅，在切變風(fēng)條件下，正偏航時(shí)，葉片扭矩的波幅增加了約57.7%（26.7kN·m），揮舞力矩波幅減小了37.8%（77.1kN·m）；負(fù)偏航時(shí)，扭矩波幅增加了10.8%（5.0kN·m），揮舞力矩波幅增加了56.6%（114kN·m）?？梢姡陲L(fēng)切變條件下，偏航風(fēng)增強(qiáng)了葉片扭矩的波動(dòng)；當(dāng)風(fēng)輪正偏航時(shí)，葉片揮舞力矩波動(dòng)減弱，負(fù)偏航時(shí)波動(dòng)增強(qiáng)。

圖12給出了切變風(fēng)條件下，偏航30°時(shí)計(jì)算得到的風(fēng)輪尾跡幾何。由頂視圖可見，由于偏航造成尾跡在來流方向上發(fā)生偏斜運(yùn)動(dòng)；由側(cè)視圖可見，由于切變風(fēng)造成尾跡發(fā)生上揚(yáng)運(yùn)動(dòng)。由于偏航和切變風(fēng)的共同作用，形成了結(jié)構(gòu)比較復(fù)雜的風(fēng)輪尾跡。

圖11 不同偏航角時(shí)葉片轉(zhuǎn)矩和揮舞力矩隨方位角的變化對比Fig.11 Variation of the computed shaft torque and flap moment of the blade with respect to the azimuth angle under yaw conditions

圖12 風(fēng)輪尾跡形狀，偏航30°，切變風(fēng)Fig.12 Shape of the rotor wake under the condition of 30°of yaw angle with shear wind

以上研究均假設(shè)風(fēng)速是定常的，風(fēng)力機(jī)實(shí)際運(yùn)行中，風(fēng)速是非定常的。由于所研究的模擬方法屬于時(shí)間步進(jìn)法，因此可以將非定常風(fēng)速作為數(shù)值模擬的邊界條件進(jìn)行研究。相關(guān)研究將會在今后的研究工作中開展。

4 結(jié) 論

基于自由渦方法，建立了一個(gè)高效模擬水平軸風(fēng)力機(jī)處于風(fēng)切變條件下，偏航氣動(dòng)特性分析的方法，并對兩種不同尺度的風(fēng)力機(jī)和不同運(yùn)行工況進(jìn)行了計(jì)算，在與實(shí)驗(yàn)測試數(shù)據(jù)比較的基礎(chǔ)上，分析了風(fēng)輪處于風(fēng)切變和不同偏航時(shí)，葉片氣動(dòng)力和扭矩的變化規(guī)律。

結(jié)果表明，所建立的計(jì)算方法可以較好的模擬風(fēng)切變條件下，風(fēng)輪偏航時(shí)葉片氣動(dòng)載荷的變化。分析發(fā)現(xiàn)，在風(fēng)切變條件下，偏航會增加葉片扭矩的波動(dòng)，當(dāng)風(fēng)輪處于正偏航角時(shí)，葉片揮舞力矩的波動(dòng)減弱，處于負(fù)偏航角時(shí)，揮舞力矩波動(dòng)增強(qiáng)。

［1］ Tongchitpakdee C，Benjanirat S，Sankar L N.Numerical simulation of the aerodynamics of horizontal axis wind turbines under yawed flow conditions［J］.Journal of Solar Energy Engineering，2005，127（4）：464-475.

［2］ Zhou Wenping，Tang Shengli.The aerodynamic performance computation of HAWT in steady yaw［J］.Acta Energiae Solaris Sinica，2011，32（9）：1315-1320.（in Chinese）周文平，唐勝利.水平軸風(fēng)力機(jī)穩(wěn)定偏航氣動(dòng)性能計(jì)算［J］.太陽能學(xué)報(bào)，2011，32（9）：1315-1320.

［3］ Qiu Yongxing，Kang Shun.Numerical simulation of wind turbine wake based on nonlinear lifting surface method［J］.Journal of Engineering Thermophysics，2012，33（12）：2072-2075.（in Chinese）仇永興，康順.基于非線性升力面方法的風(fēng)力機(jī)尾跡數(shù)值模擬［J］.工程熱物理學(xué)報(bào)，2012，33（12）：2072-2075.

［4］ Gupta S，Leishman J G.Comparison of momentum and vortex methods for the aerodynamic analysis of wind turbines［C］.43rd AIAA Aerospace Sciences Meeting and Exhibit，Reno：2005.

［5］ Huang Shuilin，Li Chunhua，Xu Guohua.An analytical method for aerodynamic interaction of twin rotors based on free-vortex and lifting-surface models［J］.Acta Aerodynamica Sinica，2007，25（3）：390-395.（in Chinese）黃水林，李春華，徐國華.基于自由尾跡和升力面方法的雙旋翼懸停氣動(dòng)干擾計(jì)算［J］.空氣動(dòng)力學(xué)學(xué)報(bào)，2007，25（3）：390-395.

［6］ Xu Bofeng，Wang Tongguang，Zhang Zhenyu，et al.Aerodynamic optimal design of winglet based on free vortex wake method and genetic algorithm［J］.Acta Aerodynamica Sinica，2013，31（1）：132-136.（in Chinese）許波峰，王同光，張震宇，等.基于自由渦尾跡和遺傳算法的葉尖小翼氣動(dòng)優(yōu)化設(shè)計(jì)［J］.空氣動(dòng)力學(xué)學(xué)報(bào)，2013，31（1）：132-136.

［7］ Sebastian T，Lackner M A.Development of a free vortex wake method code for offshore floating wind turbines［J］.Renewable Energy，2012，46（2012）：269-275.

［8］ Sebastian T.Understanding the unsteady aerodynamics and near wake of an offshore floating horizontal axis wind turbine［D］.［Ph D Thesis］.University of Massachusetts Amherst，2011.

［9］ Shen Xin，Zhu Xiaocheng，Du Zhaohui.Wind turbine aerodynamics and loads control in wind shear flow［J］.Energy，2011，36（2011）：1424-1434.

［10］Sant T.Improving BEM-based aerodynamic models in wind turbine design codes［D］.［Ph D Thesis］.Delft University of Technology，2007.

［11］Nilay Sezer Uzol，Oguz Uzol.Effect of steady and transient wind shear on the wake structure and performance of a horizontal axis wind turbine rotor［J］.Wind Energy，2011，16（1）：1-17.

［12］Leishman J G，Bagai A，Bhagwat M J.Free-vortex filament methods for the analysis of helicopter rotor wakes［J］.Journal of Aircraft，2002，39（5）：759-775.

［13］Leishman J G.Challenges in modeling the unsteady aerodynamics of wind turbines［J］.Wind Energy，2002，5（2-3）：85-132.

［14］Hand M M，Simms D A，F(xiàn)ingersh L J，et al.Unsteady aerodynamics experiment phase vi：wind tunnel test configurations and available data campaign［R］.Technical report NREL/TP-500-29955，NREL，2001.

［15］Simms D，Schreck S，Hand M，et al.NREL unsteady aerodynamics experiment in the NASA-Ames wind tunnel：a comparison of predictions to measurements［R］.Technical report NREL/TP-500-29494，NREL，2001.

［16］PF，BEP，KSH.The Tjreborg wind turbine loads during normal operation mode for CEC［R］.DG XII.DK.ELSAMPROJEKT A/S，EP94/456.1994.

［17］Yu Guohua，Shen Xin，Zhu Xiaocheng，et al.An insight into the separate flow and stall delay for HAWT［J］.Renewable Energy，2011，36（1）：69-76.

［18］Landgrebe A J.The wake geometry of a hovering rotor and its influence on rotor performance［J］.Journal of the American Helicopter Society，1972，17（4）：2-15.

［19］Massouh F，Dobrev I.Exploration of the vortex wake behind of wind turbine rotor［J］.Journal of Physics：Conference Series，2007，75：012036.

［20］Vatistas G H，Kozel V.A simpler model for concentrated vortices［J］.Experiments in Fluids，1991，11（1）：73-76.

［21］Bagai A，Leishman J G.Rotor free-wake modeling using apseudoimplicit relaxation algorithm［J］.Journal of Aircraft，1995，32（6）：1276-1285.

Aerodynamic characteristics of HAWT rotor in crosswind based on free vortex method

Qiu Yongxing1，Kang Shun1，2
（1.North China Electric Power University，Key Laboratory of CMCPPE Ministry of Education，Beijing 102206，China；2.Xi′an Institute of Modern Control Technology，Xi′an 710065，China）

In order to study the aerodynamic characteristics of HAWT rotor under yawed conditions with shear wind，an unsteady numerical simulation method is developed.It consists of a nonlinear lifting line method and a time-accurate free vortex method.To improve the convergence of the nonlinear lifting line method，an iterative algorithm，based on Newton-Raphson method，is used to solve the circulation equations.In order to improve the accuracy and the efficiency of the algorithm for rotor wake simulation，an improved wake vortex model is developed，which includes a viscous vortex core model，a vortex sheet model and a tip vortex filaments model.The NREL Phase VI and Tjreborg wind turbines are taken as the testing examples.The aerodynamic loads of both the rotor and the blades are computed under different yawed conditions with uniform wind or the shear wind.The results are compared to the experimental data.The key factors and the rules of the aerodynamic characteristics of the rotor and the blades are analyzed.The results show that the method developed in this paper can simulate the aerodynamic characteristics of HAWT effectively under the yawed conditions with shear wind.The fluctuations of shaft torque of blades are increased；the fluctuations of flap moment of blades are decreased at the plus yaw angles and increased at the minus yaw angles.

HAWT rotor；shear wind；yaw wind；unsteady numerical simulation；free vortex method

TK83

：Adoi：10.7638/kqdlxxb-2013.0025

0258-1825（2015）02-0246-08

2013-03-05；

：2013-08-07

國家自然科學(xué)基金（51176046）

仇永興＊（1985-），男，甘肅蘭州人，博士研究生，研究方向：風(fēng)力機(jī)空氣動(dòng)力學(xué)數(shù)值模擬研究.E-mail：qiu＿yx＠163.com

仇永興，康順.基于自由渦方法的風(fēng)輪偏航氣動(dòng)特性研究［J］.空氣動(dòng)力學(xué)學(xué)報(bào)，2015，33（2）：246-253.

10.7638/kqdlxxb-2013.0025 Qiu Y X，Kang S.Aerodynamic characteristics of HAWT rotor in crosswind based on free vortex method［J］.Acta Aerodynamica Sinica，2015，33（2）：246-253.