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

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

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

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

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

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

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

1 基底彎矩均值解析表達(dá)式

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

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

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

圖2 來流風(fēng)速分解示意圖Fig.2 Wind speed decomposition

圖3 平均風(fēng)速分解示意圖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ù)和升力系數(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)

每個(gè)葉片阻力和升力合力的均值分別為

(10)

(11)

式中，i=1,2,3。

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

(12)

(13)

(14)

(15)

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

(16)

(17)

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

2 案例分析

以美國國家可再生能源實(shí)驗(yàn)室(NREL)5 MW陸上風(fēng)機(jī)為例，將均值解析表達(dá)式的計(jì)算結(jié)果與OpenFAST的計(jì)算結(jié)果進(jìn)行比較，分析均值解析表達(dá)式產(chǎn)生誤差的原因。NREL 5 MW陸上風(fēng)機(jī)主要參數(shù)，如表1所示。

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

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

表2 5 MW風(fēng)機(jī)葉片截面和翼型信息表Tab.2 Blade sections and airfoil information for 5 MW wind turbine

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

2.1 風(fēng) 場

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

(18)

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

(19)

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

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

(20)

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

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

(21)

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

(22)

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

a=12,

Lc=5.67×min(60 m,H)

(23)

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

(24)

2.2 攻角檢驗(yàn)

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

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

2.3 計(jì)算結(jié)果及誤差分析

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

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

表3 風(fēng)機(jī)在不同剛?cè)嵝栽O(shè)置基底彎矩誤差結(jié)果Tab.3 Error when different stiffness and flexibility

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

表4 3個(gè)葉片上風(fēng)荷載合力均值誤差表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中的計(jì)算結(jié)果發(fā)現(xiàn)，現(xiàn)有的解析表達(dá)式計(jì)算得到的基底彎矩均值均小于OpenFAST計(jì)算得到的均值，所以按照現(xiàn)有規(guī)范設(shè)計(jì)偏危險(xiǎn)，需要進(jìn)行基底彎矩的修正。

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

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

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

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

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

(25)

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

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

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

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

(26)

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

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

(27)

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

修正后半經(jīng)驗(yàn)方法計(jì)算的基底彎矩可以表示為

(28)

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

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

表6 1.5 MW和15 MW風(fēng)機(jī)半解析方法計(jì)算結(jié)果Tab.6 Calculation results of 1.5 MW and 15 MW wind turbine by semi analytical method

4 結(jié) 論

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

(1)導(dǎo)致解析表達(dá)式與OpenFAST計(jì)算結(jié)果誤差的主要原因有兩方面，其一，結(jié)構(gòu)柔性導(dǎo)致了偏心距的改變，從而產(chǎn)生了新的偏心彎矩；其二，塔筒柔性導(dǎo)致了力臂變化，而解析表達(dá)式中無法精確計(jì)算結(jié)構(gòu)變形后各點(diǎn)風(fēng)力的實(shí)際彎矩。并且風(fēng)機(jī)結(jié)構(gòu)整體越柔，結(jié)構(gòu)變形越大，力臂的改變也越大，兩種方法計(jì)算結(jié)果的誤差就越大。

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

致謝

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