仇永興，康 順

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

基于自由渦方法的風輪偏航氣動特性研究

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

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

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

風輪；風切變；偏航風；非定常模擬；自由渦

0 引 言

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

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

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

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

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

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

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

1.1 非線性升力線方法

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

圖1 葉片微段的局部坐標系及渦系構成示意圖Fig.1 Schematic of blade element coordinate and vortex system composition

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

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

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

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

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

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

（ii）計算附著渦段中點的合速度，解出葉片攻角徑向分布；

（iii）通過葉片當?shù)匾硇偷臍鈩犹匦郧€，得到升力系數(shù)；

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

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

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

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

1.2 自由渦計算

1.2.1 渦尾跡模型

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

為了有效模擬脫體渦的運動和擾動作用，并降低尾跡計算量，建立了一種由風輪附近的自由運動的渦面（包含尾渦線和脫體渦線）和較遠下游的葉尖渦線構成的渦尾跡模型。模型假設：

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

圖2 風輪單根葉片的渦尾跡Fig.2 Vortex wake of one blade

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

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

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

為了模擬尾跡運動中的粘性效應，同時避免尾跡速度場的奇異性，在渦尾跡模型中引入Lamb-Oseen粘性渦核模型［20］：

其中：

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

1.2.2 渦尾跡控制方程及求解

渦尾跡由自由渦面和葉尖渦線構成。而構成自由渦面和葉尖渦線的是離散的渦線段。渦線段的運動控制方程為：

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

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

2 算例模型與方案介紹

為了驗證計算方法的有效性，并且研究風切變條件下風輪偏航氣動特性，選取以下兩種風力機作為算例模型：

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

研究方案具體如下：

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

式中，uhub和Hhub分別為風輪輪轂位置的風速和輪轂高度。

3 結果與分析

3.1 方法有效性驗證

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

圖3給出了均勻風條件下，NREL Phase VI風力機葉片槳矩角3°、無偏航時功率和推力計算值與實驗值的比較。圖4給出了均勻風條件下，Tjreborg風力機，葉片槳矩角0.4°，無偏航時功率計算值與實驗值［14］的比較。

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

實驗值的波動可能與實際運行中風入流存在風切變有關，下面將詳細研究。

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

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

以上計算表明，計算方法模擬風輪氣動特性是準確有效的。

3.2 風切變條件下風輪氣動特性分析

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

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

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

圖8中，葉片I運動到方位角約90°時，其轉矩為最大值；運動到約270°時為最小值。這是因為風切變條件下，風輪上半周風速較高，下半周風速較低。當風輪無偏航時，旋轉平面最高點和最低點附近分別存在葉片氣動載荷的極大、極小值。

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

圖8 風切變時風輪及葉片轉矩隨方位角的變化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 偏航條件下風輪氣動特性分析

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

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

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

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

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

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

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

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

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

圖12給出了切變風條件下，偏航30°時計算得到的風輪尾跡幾何。由頂視圖可見，由于偏航造成尾跡在來流方向上發(fā)生偏斜運動；由側視圖可見，由于切變風造成尾跡發(fā)生上揚運動。由于偏航和切變風的共同作用，形成了結構比較復雜的風輪尾跡。

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

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

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

4 結 論

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

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

［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）周文平，唐勝利.水平軸風力機穩(wěn)定偏航氣動性能計算［J］.太陽能學報，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）仇永興，康順.基于非線性升力面方法的風力機尾跡數(shù)值模擬［J］.工程熱物理學報，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）黃水林，李春華，徐國華.基于自由尾跡和升力面方法的雙旋翼懸停氣動干擾計算［J］.空氣動力學學報，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）許波峰，王同光，張震宇，等.基于自由渦尾跡和遺傳算法的葉尖小翼氣動優(yōu)化設計［J］.空氣動力學學報，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

國家自然科學基金（51176046）

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

仇永興，康順.基于自由渦方法的風輪偏航氣動特性研究［J］.空氣動力學學報，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.