黃國慶，袁 酊，劉 敏，許 楠

(1.重慶大學 土木工程學院，重慶 400045；2.中國廣核新能源控股有限公司，北京 100000)

風能是一種可再生、無污染的綠色能源，儲量豐富。我國陸地風力發(fā)電儲量合計10 億kW[1]。風電行業(yè)發(fā)展迅速，截至2019年底，我國陸上風電累計裝機容量位列全球排名首位，累計陸上裝機容量占全球陸上風電裝機容量的37%。為提高風機的經濟效益，單機容量變得越來越大，目前正向5 MW以上發(fā)展[2]。

當前風機的選型主要按照國際電工委員會發(fā)布的標準IEC 61400-1[3]進行。根據(jù)風機安裝場地的風速和湍流度將場地分為Ⅰ級、Ⅱ級、Ⅲ級和S級以便選擇風機葉片型號。但這種場地分級方法只能大致將實際場地的風速和湍流參數(shù)與上述標準分類匹配，而實際風場與規(guī)范規(guī)定風場存在的差異無法考慮。因此在精細化設計，特別是塔筒設計中還需要大量分析。在設計塔筒時，可以用數(shù)值仿真軟件進行估計來獲得比較準確的響應，但是數(shù)值仿真需要進行耗時的時程分析才能確定動力響應的統(tǒng)計特性，計算效率較低，難以適應風機快速、精細選型的需求。因此急需更加高效的風機塔筒響應計算方法。

目前設計規(guī)范中常采用等效靜風荷載將復雜的動力分析問題轉換為靜力問題來評估風致結構響應(如構件的彎矩、軸力和變形等)。借鑒這一思想，Tempel[4]將風力機系統(tǒng)等效為一個單自由度體系，通過理論方法計算獲得簡化的風振動力響應。賀廣零[5]提出了兩自由度體系陣風荷載因子，其安全性較單自由度體系陣風荷載因子更高。柯世堂等[6]通過建立一體化有限元模型，揭示了風力發(fā)電塔輪系統(tǒng)風振響應背景與共振部分的分布特性和模態(tài)間耦合效應的作用機理。Xu等[7]推導了風機停機狀態(tài)下順風向和橫風向基底彎矩及其極值的解析表達式。該解析表達式由均值和脈動兩部分組成。在Xu等的研究中，基底彎矩極值的解析表達式是在準定常假設的基礎上推導的，需要將葉片上各個截面的風荷載在葉片長度方向進行積分，形式復雜，不便于直接應用于設計。為了便于工程應用，Xu等[8]對多個不同型號的風機進行數(shù)值仿真后，通過參數(shù)擬合獲得了簡化的基底彎矩極值解析表達式，并將簡化公式計算的響應與數(shù)值仿真得到的響應進行了對比[9]，結果表明在400 kW的小型風機上該簡化公式的適用性較好。該簡化公式可快速簡便地獲得響應預估結果，被日本《風力發(fā)電設備塔架結構設計指南及解說》[10]所采納。

基底彎矩的均值在基底彎矩極值中占比可達50%，所以準確地估計其均值是計算基底彎矩極值的基礎。需要指出的是在Xu等的研究中，計算基底彎矩均值時不考慮結構變形對基底彎矩的影響。該假設對于小型風機是合理的，但對于大型兆瓦級風機，由于葉片長度和塔筒高度都顯著增加，結構整體剛度明顯降低，葉片變形和塔筒變形相比小型風機明顯增加，導致更大的基底彎矩等響應。因此基于柔性塔筒和剛性葉片假設的均值解析表達式對5 MW及以上的大型風機的適用性還有待進一步驗證。

限于篇幅，本文以塔筒基底彎矩的均值為研究對象。針對風機的不利工況之一的停機狀態(tài)，本文研究了塔筒響應均值計算解析表達式。首先基于以往學者的工作重新推導了風機基底彎矩均值解析表達式，然后以5 MW的變槳距控制風機為例，將均值解析表達式的計算結果與OpenFAST軟件的計算結果進行比較，探討了誤差產生的原因。最后提出了考慮葉片和塔筒變形的基底彎矩均值半解析方法，以期減小解析公式計算誤差，為風機塔筒響應均值的快速計算提供更精確的方法。

1 基底彎矩均值解析表達式

根據(jù)日本《風力發(fā)電設備塔架結構設計指南及解說》，對于變槳距控制風力機，當停機狀態(tài)下作用于塔筒的最大風荷載發(fā)生在偏航角為90°，槳距角為90°，葉片1方位角為90°時，如圖1所示。

圖1 順槳狀態(tài)下變槳距控制風力機的姿勢Fig.1 Parking state of pitch control wind turbine

來流風速包括水平方向分量平均風速U0和脈動風速u0以及豎向、橫向的脈動風速w0和v0，如圖2所示。由于各翼型的升阻力系數(shù)是基于垂直于葉片長度方向的截面定義的，因此計算葉片上作用的風荷載時需要將來流風速分解到沿葉片長度方向和垂直葉片長度方向。橫向脈動風速v0始終垂直于葉片長度方向，因此v0不需要分解。將順風向風速U0+u0與豎向脈動風速w0的矢量組合速度分解為沿葉片長度的分量W1和垂直葉片長度方向的分量V1。在圖2中，U0+u0,w0,W1和V1處于同一平面內，其中分量W1引起沿葉片長度方向的荷載幾乎為0，V1由均值部分U和脈動部分u組成，即V1=U+u。V1的均值部分U由平均風速U0的分解得到，如圖3所示。葉片2相對于風向是傾斜的，葉片上平行于橫截面的平均風速為U=U0cosβ，其中，β為葉片的方位角。V1的脈動部分u由脈動風速u0和w0分解得到，u=u0cosβ-w0sinβ。葉片3上風速分解的方法與葉片2的分解方法類似，不予贅述。

圖2 來流風速分解示意圖Fig.2 Wind speed decomposition

圖3 平均風速分解示意圖Fig.3 Mean wind speed decomposition

圖4 葉片2截面所受升阻力示意圖Fig.4 Lift force and drag force of a section of blade 2

(1)

(2)

式中：ρ為空氣密度；CD(θ+θ′,r)和CL(θ+θ′,r)分別為在距離葉根r處的某截面在攻角為θ+θ′時的阻力系數(shù)和升力系數(shù)；c(r)為在距離葉根r處葉片的弦長。

(3)

將D和L分解到x,y方向，并在θ處泰勒展開，可得阻力fD為

fD(r)=D(θ+θ′,r)cosθ′-L(θ+θ′,r)sinθ′≈

(4)

相似的，升力fL可以表示為

(5)

(6)

(7)

(8)

(9)

每個葉片阻力和升力合力的均值分別為

(10)

(11)

式中，i=1,2,3。

計算塔筒基底彎矩時，先將分布在葉片上的風荷載簡化為作用在輪轂高度處的集中荷載與附加彎矩，輪轂處集中荷載可以表示為

(12)

(13)

(14)

(15)

在實際情況中，塔筒底部的彎矩還包括由于風機質量的不對稱分布導致的偏心彎矩，但本文中所述基底彎矩只考慮由等效到輪轂處的葉片總風荷載和塔筒上風荷載導致的彎矩。塔筒底部由順風向風荷載和橫風向風荷載產生的彎矩均值分別為

(16)

(17)

式中：H為輪轂高度；d(z)為塔筒直徑；CD(z)為塔筒截面阻力系數(shù)；CL(z)為塔筒截面升力系數(shù)；U(z)為沿高度變化的平均風速。

2 案例分析

以美國國家可再生能源實驗室(NREL)5 MW陸上風機為例，將均值解析表達式的計算結果與OpenFAST的計算結果進行比較，分析均值解析表達式產生誤差的原因。NREL 5 MW陸上風機主要參數(shù)，如表1所示。

表1 5 MW風機計算中所用主要參數(shù)Tab.1 Some parameters of 5 MW wind turbine

5 MW風機的每個葉片有19個控制截面，對應8個不同的翼型，各控制截面與翼型的對應關系，如表2所示。其中部分翼型的升力系數(shù)和阻力系數(shù)隨攻角的變化，如圖5所示。

表2 5 MW風機葉片截面和翼型信息表Tab.2 Blade sections and airfoil information for 5 MW wind turbine

圖5 風機葉片部分截面升力系數(shù)和阻力系數(shù)隨攻角變化Fig.5 Lift coefficient and drag coefficient for different airfoils

2.1 風 場

平均風速剖面采用指數(shù)率，U(z)=UH(z/H)α1，H為輪轂高度，UH為輪轂高度處平均風速，本案例中設置為30 m/s；α1為地面粗糙度指數(shù)。湍流度設置為0.18，不考慮湍流度沿高度的變化。脈動風速u0,v0,w03個方向的自功率譜均采用Kaimal譜，由式(18)給出

(18)

式中:k=u0，v0,w0；f為頻率；Lk為各脈動風速成分的積分尺度，其取值由IEC規(guī)范給出

(19)

式中，ΛU為湍流積分尺度，其定義為

ΛU=0.7·min(60 m,H)

(20)

式中，min(60 m,H)為取輪轂高度H與60 m中較小者。標準差之間的關系定義為

σv0=0.8σu0,σw0=0.5σu0

(21)

相干函數(shù)cohi,j采用風力發(fā)電機設計要求(IEC規(guī)范)中u0方向風速分量的相干函數(shù)

(22)

式中：d為兩點間的向量在與平均風向垂直的平面上的投影線的長度；a為衰減系數(shù)；Lc為相干尺度參數(shù)。在IEC規(guī)范中，參數(shù)a和Lc的取值為

a=12,

Lc=5.67×min(60 m,H)

(23)

IEC規(guī)范中沒有指定橫風向風速分量和豎向風速分量的相干函數(shù)。本文采用OpenFAST進行風場模擬時采用的橫風向風速分量和豎向風速分量的相干函數(shù)[12]

(24)

2.2 攻角檢驗

從圖5中可發(fā)現(xiàn)升力系數(shù)和阻力系數(shù)均對攻角敏感，對攻角的校核是風荷載計算的關鍵一步，故有必要對比OpenFAST和解析表達式計算的攻角均值。圖6給出了葉片2由兩種方法計算的葉片19個控制截面的平均攻角對比情況?？梢钥吹絻煞N方法得到的各截面攻角吻合良好，說明了采用解析表達式計算的攻角是正確的。

圖6 OpenFAST和解析表達式計算的葉片截面的攻角θ均值對比Fig.6 Comparison of mean value of blade attack angle calculated by OpenFAST and analytic method

2.3 計算結果及誤差分析

該停機工況在OpenFAST中運行一條600 s時程需要約10 min，為得到準確的統(tǒng)計量需要計算大量時程。使用解析表達式編程后大幅提升效率。

為了探究風機結構剛度對基底彎矩均值計算的影響，通過調整葉片和塔筒的自由度，分別設置了如表3所示的4個計算工況。在使用OpenFAST計算時，通過TurbSim模塊模擬20個風場時程，每個風速時程長度為630 s，然后計算風機響應時程。計算基底彎矩均值時，為了避免瞬態(tài)效應，前30 s未使用，將這20個時程的均值再取平均后作為該工況的響應統(tǒng)計量，計算結果如表3所示。

表3 風機在不同剛柔性設置基底彎矩誤差結果Tab.3 Error when different stiffness and flexibility

從表3中可得出以下結論。①當塔筒設置為剛性時，比較工況1和工況2，葉片的剛柔性對風荷載引起的基底彎矩的影響不大。②與工況1剛性塔筒剛性葉片相比，工況3和工況4中柔性塔筒工況下，由于塔筒頂部在風作用下產生位移，導致實際的偏心距與剛性塔筒工況中的初始偏心距略有差別，從而產生新的偏心彎矩值。而實際情況中的基底總彎矩中包含了風力彎矩、偏心彎矩和附加彎矩三部分(OpenFAST獲得的彎矩包含各種荷載作用的響應)。在從中得到風力彎矩時，由于偏心彎矩發(fā)生變化，可能導致風力彎矩計算有所偏差。③在保持葉片為剛性的前提下，塔筒的剛度變化對結果有一定的影響，見工況1和工況3；計算中發(fā)現(xiàn)在剛性塔筒和柔性塔筒工況中，葉片上升力和阻力的合力相差不大(見表4所示工況1和工況3的結果對比)，說明塔筒變形導致的力臂變化可能是導致基底彎矩差別的主要原因。④當塔筒設置為柔性時，如工況3和工況4所示，柔性葉片工況產生了比剛性葉片工況更大的彎矩，可能的原因是：風機整體結構越柔，變形越大，力臂的變化也越大(圖7顯示葉片2各截面沿U方向的位移均值)。

表4 3個葉片上風荷載合力均值誤差表Tab.4 Error of mean value of wind load on 3 blades

圖7 葉片2各截面沿U方向的位移均值Fig.7 Mean displacement of blade 2 along the U direction

3 基底彎矩均值修正系數(shù)

從表3中的計算結果發(fā)現(xiàn)，現(xiàn)有的解析表達式計算得到的基底彎矩均值均小于OpenFAST計算得到的均值，所以按照現(xiàn)有規(guī)范設計偏危險，需要進行基底彎矩的修正。

根據(jù)IEC規(guī)范，輪轂高度處湍流度的變化范圍約為0.11～0.18。柔性塔筒柔性葉片工況下，輪轂高度處平均風速UH=30 m/s，風剖面指數(shù)α1=0.2，湍流度分別為0.10和0.18時順風向和橫風向基底彎矩的均值，如表5所示。

表5 風機在不同湍流度下基底彎矩均值的相對差值Tab.5 Error of response mean under different turbulence

基底彎矩的誤差與輪轂高度處平均風速、風剖面指數(shù)和湍流度相關，但從表5中發(fā)現(xiàn)，在不同湍流度下基底彎矩的差別不大。因此，基底彎矩均值的修正系數(shù)為輪轂處平均風速UH和風剖面指數(shù)α1的函數(shù)。

IEC規(guī)范將地面粗糙度分成了4個等級，分別對應的風剖面指數(shù)為0.10，0.15，0.20，0.27。由于在平均風速大于等于25 m/s時風機進入停機狀態(tài)，同時在IEC規(guī)范中，以輪轂高度處的平均風速將風機的建設等級分為了3類，其中Ⅰ類風機的參考平均風速為50 m/s，因此將風速的計算范圍確定為25～50 m/s。本文計算了在α1分別為0.10，0.15，0.20和0.27時的不同地面粗糙度場地上，輪轂高度處湍流度為0.18，輪轂處平均風速從25 m/s變化至50 m/s時，解析表達式與OpenFAST的順風向和橫風向基底風力彎矩均值誤差。

定義順風向和橫風向基底彎矩均值修正系數(shù)分別為ΔD和ΔL，表示為

(25)

在不同的平均風速下順風向和橫風向基底彎矩修正系數(shù)，如圖8和圖9所示。

圖8 順風向風力彎矩修正系數(shù)Fig.8 Correction factor of along-wind bending moment

圖9 橫風向風力彎矩修正系數(shù)Fig.9 Correction factor of across-wind bending moment

順風向平均風力彎矩修正系數(shù)的表達式為

(26)

橫風向平均風力彎矩修正系數(shù)的表達式為

ΔL(α1,UH)=-0.002 5UH+0.03α1+1.22

(27)

式中，α1為0.10,0.15,0.20,0.27。

修正后半經驗方法計算的基底彎矩可以表示為

(28)

圖10 修正后順風向基底彎矩誤差Fig.10 Error of along-wind bending moment after correction

圖11 修正后橫風向基底彎矩誤差Fig.11 Error of across-wind bending moment after correction

表6 1.5 MW和15 MW風機半解析方法計算結果Tab.6 Calculation results of 1.5 MW and 15 MW wind turbine by semi analytical method

4 結 論

本文首先基于以往學者的工作，重新推導了停機狀態(tài)下風機塔筒底部彎矩均值的解析表達式，然后以美國國家可再生能源實驗室的5 MW風機為例，對解析表達式和OpenFAST軟件的計算結果進行了對比，檢驗了解析表達式對大型風機的適用性；最后針對大型風機基底彎矩均值的計算提出了修正方法。主要有以下結論：

(1)導致解析表達式與OpenFAST計算結果誤差的主要原因有兩方面，其一，結構柔性導致了偏心距的改變，從而產生了新的偏心彎矩；其二，塔筒柔性導致了力臂變化，而解析表達式中無法精確計算結構變形后各點風力的實際彎矩。并且風機結構整體越柔，結構變形越大，力臂的改變也越大，兩種方法計算結果的誤差就越大。

(2)與OpenFAST計算的結果相比，解析表達式計算的基底彎矩總是偏小，直接使用該解析表達式偏危險。因此本文提出了基底彎矩均值計算的半解析方法，提高了計算精度。本文中用5 MW風機的另一計算工況和1.5 MW，15 MW風機驗證了所提出修正公式的可行性。

致謝

感謝國家自然科學基金(51778546)和111引智基地項目(B18062)對本研究的支持。