極端服役環境下的風電機組塔架結構參數優化研究

戴巨川 劉 旋 楊書儀 文澤軍 沈祥兵

1.湖南科技大學，湘潭，411201 2.湘電風能有限公司，湘潭，411000

0 引言

風能是一種可再生的清潔能源，是各國新能源開發和利用的戰略選擇，近年來對它的開發利用發展迅速。塔架作為風電機組的支撐結構，其穩定性與可靠性直接影響整個機組的運行安全，塔架成本昂貴，制造成本占整機成本的15%～20%。因此，近年來越來越多的學者對塔架設計優化進行了研究：Negm等［1］以橫截面面積、壁厚和高度為變量，最小塔架質量、最大剛度和最大自然頻率為目標進行了塔架優化設計；Horvath等［2］以塔架等應力分布為目標，采用有限元方法（FEM）和計算流體力學（CFD）方法探討了塔架高度和橫截面參數取值范圍；Uys等［3］同時考慮材料成本與制造成本，以綜合成本最小為目標，獲得了塔架壁厚、加強筋數目和直徑等的優化值；Clifton－smith等［4］運用有限元法研究了塔架的屈曲穩定性，對屈曲應力計算公式進行了修正；Hassan等［5］以塔架固有頻率最大為目標，塔架質量為約束，對塔架的壁厚變化方式進行了優化設計；Yildirim等［6］以塔架質量最小為目標，采用遺傳算法進行了優化求解；趙吉文等［7－8］以固有頻率最大為目標，采用支持向量機建立了塔架結構參數與固有頻率之間的快速計算模型，引入粒子群優化、混沌搜索等方法進行優化求解。盡管如此，針對風電機組塔架的優化設計仍有許多問題需要解決，如獲取塔架直徑、壁厚等結構參數與應力、固有頻率及壽命等的關聯特性，尋求合理的替代模型解決優化求解時需迭代計算帶來的耗時長、效率低等問題。

本文采用一種壁厚分段式線性變化塔架結構模型，借助有限元數值仿真分析結構參數與應力、固有頻率間的關聯特性，采用支持向量機構建結構參數與應力、固有頻率間的快速計算模型。

1 塔架幾何模型

大型風電機組塔架大多采用按一定錐度變化的變截面圓筒形結構。為方便運輸，通常將塔架分為幾節，安裝時用法蘭將其連接成整體。考慮到塔架應力分布與圓筒半徑、壁厚等為非線性關系，為盡可能地改善塔架應力分布，本文采用分段式線性變化塔架結構，分為三節，如圖1所示。塔架主要結構參數包括底端直徑、頂端直徑、塔架總高度、各節塔架高度及厚度。由于風電機組在設計時一般以輪轂高度處的風速為設計依據，而塔架總高度主要由風電機組容量、風場風速和所處地理位置等決定，故進行塔架設計優化時該值視為常量；同時考慮到塔架頂端與機艙配合，頂端直徑值也可視為常量。

圖1 塔架結構示意圖

2 風電機組載荷分析

對風電機組進行載荷分析，獲得準確的塔架載荷是進行塔架受力分析、優化設計的基礎。風電機組承受的載荷特別是氣動載荷非常復雜，由于來流的時變性，作用在風輪上的氣動載荷在風電機組運行過程中是時刻變化的，是影響風電機組塔架載荷變化的主要因素。為便于分析，本文采用風電機組坐標系如圖2所示，葉片上某一截面處的載荷分解為坐標系xeyeze中沿xe軸方向分力Fxe和繞xe軸的力矩Mxe（揮舞力矩），沿ye軸方向分力Fye和繞ye軸的力矩Mye（擺振力矩），沿ze軸方向分力Fze和繞ze軸的力矩Mze（變槳距力矩）。

圖2 風電機組坐標系［9］

圖2 中，γ為風輪錐角，τ為風能仰角。坐標系xgygzg為慣性坐標系，固定于塔架底部；坐標系xnynzn固定在機艙上，其yn軸在風電機組主軸軸線上；坐標系xsyszs固定在主軸上，與主軸同步旋轉，轉速為n；xpypzp為假定無風能仰角時，固定于旋轉主軸上的坐標系；坐標系xeyeze固定在旋轉的葉片1上，ze軸沿葉片變槳軸線方向，ye軸垂直于ze軸指向塔架（上風向風電機組），xe軸垂直于葉片變槳軸線和主軸軸線。

目前，進行風電機組氣動載荷分析有多種方法，如葉素—動量（BEM）理論、基于動態入流理論的GDW方法、CFD等。在實際工程中，由于BEM理論簡便可信，故是常用的方法。以下僅對基于BEM理論的氣動載荷計算方法進行簡約介紹。在BEM理論中，取半徑r處長度為dr的葉素為對象，氣流速度與空氣動力分量如圖3所示［9－10］。

圖3 葉素上氣流速度與空氣動力分量

圖3中，A為前緣點，C為壓力中心，B為后緣點。v0為軸向風速，v1為氣流相對于葉素的速度；ω為風輪角頻率；a為軸向誘導速度系數，a′為切向誘導速度系數；α為攻角，β為節距角，φ為來流角度；dR為葉素上的氣動力，dL為dR的升力分量，dD為dR的阻力分量，dQ為氣動力切向分量（與旋轉面相切），dN為氣動力軸向分量；ve－Tip和ve－Top分別是塔架振動引起的x方向和y方向速度，ve－Bip和ve－Bop分別是葉片振動引起的x方向和y方向速度。對于葉片上不同半徑r處的葉素而言，ve－Tip和ve－Top不變，而ve－Bip和ve－Bop隨當前葉素半徑r的變化而變化。

作用在元素上的升力和阻力分別為

式中，c為半徑r處的葉片弦長；Cl、Cd分別為升力系數和阻力系數；ρ為空氣密度。

dL和dD的合力dR可以分解為軸向力dN和切向力dQ：

產生的擺振方向彎矩為

產生的揮舞方向彎矩為

式中，Ct、Cn分別為切向力系數和法向力系數。

設v為風輪中心高度h0處的來流風速，如果風速分布服從赫爾曼Holman規律［11］，則有

式中，r′為風速廓線指數，與地面粗糙度有關；ψ為葉片方位角。

沿葉片方向將整個積分空間［0，R］分成若干個子空間［ri，ri＋1］，i＝1，2，…，N，葉片產生的擺振方向總彎矩可寫為

式中，ωBx、ωBy分別為x、y方向葉片的振動角速度。

同樣，葉片產生的揮舞方向總彎矩可以寫為

葉片擺振方向剪力為

葉片揮舞方向剪力為

同一空間向量，在不同坐標系中坐標值不同，通過變換矩陣可以實現相互轉換，變換矩陣的獲取可采用方向余弦或歐拉方法。這樣就可以把風輪上的氣動載荷轉化至塔架上。

除了風輪上氣動載荷以外，塔架自身上風載和風電機組運行載荷（變槳距）、慣性載荷等對塔架也有影響，本文對此未做介紹，相關理論可參考有關文獻。在具體確定設計載荷時，《風力發電機組規范》（中國船級社，2008）規定以具體的裝配、吊裝、維修、運行狀態或設計工況同外部條件的組合為依據，應考慮到以合理概率出現的相關載荷，如極端外部條件與故障工況相關，兩者應組合為一種載荷情況。對于某2MW風電機組，作用在塔架基礎上的極限載荷按表1給定（由相關廠家提供）。

表1 塔架極限載荷

3 結構參數與應力、固有頻率關聯特性分析

塔架結構參數對應力、固有頻率、屈曲穩定性、疲勞壽命和動態響應特性等都有直接影響，獲取結構參數與這些屬性的關聯特性是塔架優化設計的前提。一般來說，塔架在極端服役環境下發生破壞的可能性較大，因此本文選取塔架結構參數與應力、固有頻率間的關系進行分析，不考慮其疲勞壽命和動態響應的影響。

3.1 結構參數與應力關聯特性

在有限元分析軟件ANSYS中采用殼單元進行塔架建模，考慮到葉片、輪轂和機艙等部件質量影響塔架固有頻率，將葉片、輪轂和機艙簡化為一集中質量點，即mass21單元節點，用mpc184無質量剛性梁單元連接mass21單元與塔架頂端節點，材料屬性的彈性模量定義為205GPa，泊松比為0.3，材料密度為7.85×103kg／m3，塔架壁厚值大小通過設置殼單元實常數實現，塔架底部節點全約束作為邊界條件。塔架基礎上的載荷是塔架頂部存在的推力產生的，塔架頂部推力大小可由基礎彎矩除以塔架高度得到，在mass21節點上施加等效載荷后可求解塔架的最大應力。考慮設計的具體情況，本文選取上段壁厚δ1變化范圍為12～18mm，中段壁厚δ2變化范圍為18～26mm，下段壁厚δ3變化范圍為26～45mm，下段塔架長度H1變化范圍為15～20m，中段塔架長度H2變化范圍為20～35m，塔架底端直徑D變化范圍為3.8～4.4m，由于不同結構參數的數量級不同，為便于比較各結構參數與塔架應力的關系，對結構參數進行歸一化處理，得到如圖4所示的塔架結構參數與最大應力的關系曲線。為觀察不同塔架結構參數組合下最大應力點位置，選取4組結構參數求解其最大應力，最大應力及其節點位置如表2所示。

圖4 塔架最大應力與結構參數的變化關系

表2 最大應力值與節點位置

從圖4可知，不同結構參數對塔架的影響程度不同，在圖4a中，隨著δ1增大，最大應力值顯著減小，但是當δ1增加到一定程度以后，對應力的影響逐漸放緩；δ2和δ3的變化對最大應力值沒有顯著影響，總體趨勢比較平緩。從圖4b可以看出，隨著塔架單段長度、塔架底端直徑增大，最大應力值總體趨勢均是逐漸減小。在表2中，4組不同塔架結構參數對應的最大應力值分別為153MPa、195MPa、177MPa和209MPa，分別位于距塔架底端48m、38m、18.5m和6m處。不難看出，塔架結構參數組合不但影響塔架的最大應力，而且影響了最大應力點的位置。

3.2 結構參數與固有頻率關聯特性

風電機組在運行過程中，葉輪、塔架等相互激勵，有可能發生共振引發機組倒塌事故，因此在設計時，應充分考慮各子系統固有頻率，避免共振的發生。接近風輪旋轉頻率的塔架一階固有頻率最易導致共振，對其采用有限元模態分析法進行求解。需要指出的是，在進行模態分析時，葉片、輪轂和機艙等部件質量不能忽略，采用集中質量的方式簡化，建模時在集中質量重心處設置一個質量節點單元mass21，用MPC184無質量剛性梁單元連接質點單元與塔頂節點，以底部節點全約束為邊界條件求解塔架固有頻率。為便于比較各結構參數與塔架固有頻率的關系，對結構參數進行歸一化處理。圖5所示為結構參數與塔架一階固有頻率的關系曲線，表3所示為4組結構參數不同時塔架的固有頻率。

圖5 塔架一階固有頻率與結構參數關系

表3 塔架固有頻率

從圖5可知，隨著各結構參數的增大，固有頻率均呈增大趨勢；不同結構參數對塔架固有頻率的影響程度不同，在所分析的結構參數中上段壁厚δ1和下段直徑D的變化對塔架的固有頻率影響較大，下段壁厚δ3和中段塔架長度H2影響較小。從表3中可以看出，不同結構參數組合下，塔架固有頻率值差異比較明顯，表中4組不同結構參數下的塔架一階固有頻率變化范圍為0.372 09～0.532 96Hz，而風電機組運行時風輪最大旋轉頻率為0.375Hz左右，第3組、4組參數下塔架易發生共振。

4 結構參數與應力、固有頻率快速計算模型

塔架結構參數與應力、固有頻率之間是一種復雜非線性關系，近似的數學表達式計算誤差大，應該采用有限元計算結果。由于在優化設計時需要反復多次計算，并且優化算法存在不確定性，一次優化結果一般不能作為最終結果，需要多次優化后進行結果比較確定，直接采用有限元分析軟件計算效率較低，尋求合理的替代模型可以有效解決這一問題。在各種人工智能方法中，支持向量機在解決小樣本、非線性及高維模式識別中表現出許多特有的優勢，并能推廣到函數擬合等其他機器學習問題中。本文采用支持向量機建模方法建立塔架結構參數與應力、固有頻率間的快速計算替代模型。

4.1 樣本數據獲取

為了構建基于支持向量機的結構參數與應力、固有頻率的快速計算模型，首先需要獲得足夠的樣本數據。采用正交試驗設計與隨機組合方法相結合的混合試驗設計方法確定參數組合，通過有限元數值仿真試驗獲得樣本值。根據塔架設計中可變的結構參數確定正交試驗的6個因素，即因素A（塔架上段壁厚δ1）、因素B（塔架中段壁厚δ2）、因素C（塔架下段壁厚δ3）、因素D（塔架下段長度H1）、因素E（塔架中段長度H2）、因素F（塔架底端直徑D）；根據塔架結構參數與最大應力、固有頻率的變化關系確定每個因素的5個水平值。試驗因素水平表見表4，按照正交試驗方法需要進行L25（56）＝25組試驗。除正交試驗獲得的25個樣本數據外，再從整個排列組合156 25種試驗方案中隨機選取200組進行試驗，共獲得225個樣本數據。

表4 塔架結構參數因素水平表

4.2 快速計算模型的構建

樣本空間為 ｛（x1，z1），（x2，z2），…，（x225，z225）｝，其中，xi＝（δ1i，δ2i，δ3i，H1i，H2i，Di），為塔架結構參數；zi＝（σi，fi），為塔架應力和固有頻率；i＝1，2，…，225。支持向量機中采用的方法是將輸入向量xi通過某種非線性關系映射到高維空間，則在高維空間回歸函數為

式中，φ為一個從輸入空間到高維特征空間的映射函數；b為偏差值；w為權向量。

式（12）即為待求的塔架結構參數與應力或固有頻率替代模型。設擬合精度為ε，求解最佳回歸函數轉換為固定間隔下尋找最小的‖w‖問題。采用ε－SVR型支持向量，考慮樣本的容錯性，給擬合精度ε這一硬性值加一松弛變量ξ，則有［12］：

式中，C為懲罰因子為松弛變量。

使用Lagrange乘子法解決上述問題，得到原問題的對偶問題為

對偶問題完全根據樣本數據表達，由樣本數據求得最優的Lagrange乘子后即可得到問題的最優解。則式（12）可寫為

將通過數值試驗獲得的225組樣本數據分為兩部分，200組用于訓練，其余用于測試。支持向量機類型選用ε－SVR類型，核函數設置類型選用徑向基核函數（RBF）。在MATLAB中分別用網格搜索法、粒子群方法和遺傳算法函數求解核函數參數，選取平均平方誤差（MSE）最小、平方相關系數最接近1的核參數，最終應力模型中取核參數c＝68.3341，g＝0.083 923，固有頻率模型中取核參數c＝24.5844，g＝0.003 814 7。

基于上述方法得到的應力回歸模型由179個支持向量組成，偏差值為－0.5034，固有頻率回歸模型由153個支持向量組成，偏差值為0.0192。表5給出了5組模型輸出與有限元仿真結果，可以看出構建的模型最大應力輸出值與有限元仿真結果誤差為－0.004%～4.7%，一階固有頻率輸出值與有限元仿真結果誤差為－0.005%～0.3%。結果表明，構建的塔架結構參數與最大應力、一階固有頻率的替代模型可用于工程實際。

表5 模型輸出值與有限元仿真結果

5 塔架結構參數優化

5.1 參數優化數學模型

選取塔架的結構參數為設計變量為

以塔架總質量最小為目標，即

式中，mtower為塔架總質量，mi為每段塔架的質量（可由塔架結構參數求得）。

約束條件包括：

（1）強度約束。為保證塔架的安全，外載荷作用下的塔架最大應力σ應小于材料的屈服應力，即：

式中，nst為安全系數；σ為材料屈服應力。

（2）固有頻率約束。風電機組運行時，如果風輪旋轉頻率或葉片通過頻率接近塔架的固有頻率，則會出現共振現象，引發嚴重的安全問題。某型風電機組運行時風輪的最大轉速為22.5r／min，風輪旋轉頻率和葉片通過頻率分別為0.375Hz和1.125Hz，根據GB18451.1－2001規定，塔架工作狀態下的彎曲固有頻率f0，n和激振頻率fR、fR，m之間應有適當的間隔，按下述公式計算：

式中，fR為正常運行范圍內風輪的最大旋轉頻率；f0，1為塔架（整機狀態下）的一階固有頻率；fR，m為m個風輪葉片的通過頻率；f0，n為塔架（整機狀態下）的n階固有頻率。

由式（19）得到塔架一階固有頻率約束為

（3）邊界約束。根據工程具體情況，塔架的結構參數在一定范圍內取值，這里選取變量的范圍為

本文采用遺傳算法進行參數優化求解，遺傳算法是一種基于生物遺傳和進化機制的自適應概率優化算法，具有很好的全局尋優能力。對于一些非線性、多模型、多目標的函數優化問題用其他優化方法較難求解，而遺傳算法可以方便地得到較好的結果。種群中的個體數目設為30，采用浮點編碼方式，交叉概率為0.95，變異概率為0.08，最大迭代次數為250次。

5.2 優化結果

基于上述方法，得到塔架結構優化參數如表6所示，同時給出了優化前后塔架最大應力和固有頻率的有限元數值計算結果，優化后的塔架總質量為100.4t，相比優化前的104t，用鋼量減少3.5%（3.6t），一階塔架固有頻率從0.4788Hz變為0.4793Hz，最大應力從124MPa降為120MPa。圖6給出了優化前后塔架應力分布曲線，優化后應力變化更平緩，塔架的安全性有所提高。

表6 塔架結構參數優化結果

圖6 塔架應力分布曲線

6 結論

（1）通過試驗設計和有限元數值仿真分析了各參數變化對應力、固有頻率的影響趨勢，獲得了塔架結構參數與應力、固有頻率的關聯特性，在塔架設計時可為各結構參數的分布提供一定的參考。

（2）采用支持向量機理論分別構建了塔架結構參數與應力、塔架結構參數與固有頻率間的快速計算模型，計算結果顯示，替代模型輸出結果能與有限元仿真結果較好吻合，表明本文構建的快速模型是可行的。

（3）構建了某2MW風電機組以塔架質量最小為目標的塔架結構參數優化設計模型，采用遺傳算法進行優化求解，優化后的塔架質量減少3.6t，減少用鋼量3.5%，為風電機組塔架實際工程設計優化提供了一種參考方法。

［1］Negm H M，Maalawi K Y.Structural Design Optimization of Wind Turbine Towers［J］.Computers ＆Structures，2000，74（6）：649－666.

［2］Horvath G，Toth L.New Methods in Wind Turbine Tower Design［J］.Wind Engineering，2001，25（3）：171－178.

［3］Uys P，Farkas J，Jarmai K，et al.Optimisation of a Steel Tower for a Wind Turbine Structure［J］.Engineering Structures，2007，29（7）：1337－1342.

［4］Clifton－smith M，Wood D.Optimisation of Selfsupporting Towers for Small Wind Turbines［J］.Wind Engineering，2010，34（5）：561－578.

［5］Hassan A，Maalawi K，Negm H.Frequency Optimization of Conical Shells Under Mass Equality Constraint［J］.JASR，2006.

［6］Yildirim S，?zkol I·.Wind Turbine Tower Optimization under Various Requirements by Using Genetic Algorithm［J］.Engineering，2010，2（8）：641－647.

［7］趙吉文，謝芳，劉永斌，等.600kW風力機塔架結構參數粒子群優化設計［J］.中國機械工程，2007，18（15）：1776－1779.Zhao Jiwen，Xie Fang，Liu Yongbin，et al.Optimization Design of Tower Structure of 600kW－Wind Turbine Based on Particls Swarm［J］.China Mechanical Engineering，2007，18（15）：1776－1779.

［8］趙吉文，劉永斌，陳軍寧，等.基于混沌搜索的600kW風力機塔架結構優化設計［J］.系統仿真學報，2007，19（11）：2425－2428.Zhao Jiwen，Liu Yongbin，Chen Junning，et al.Optimization Design of Tower Structure of 600kWWind Turbine Based on Chaotic Search［J］.Journal of System Simulation，2007，19（11）：2425－2428.

［9］Dai Juchuan，Hu Yanping，Liu Deshun，et al.Aerodynamic Loads Calculation and Analysis for Large Scale Wind Turbine Based on Combining BEM Modified Theory with Dynamic Stall Model［J］.Renewable Energy，2011，36（3）：1095－1104.

［10］戴巨川，胡燕平，劉德順，等.MW 級變槳距風電機組葉片轉矩計算與特性分析［J］.太陽能學報，2010，31（8）：1030－1036.Dai Juchuan，Hu Yanping，Liu Deshun，et al.Torque Calculation and Characteristics Analysis of Blade for MW Adjustable－ptich Wind Turbine［J］.Acta Energiae Solaris Sinica，2010，31（8）：1030－1036，2（11）：821－829.

［11］Gourieres D L.Wind Power Plants：Theory and Design［M］.Oxford：Pergamon Press，1982.

［12］Vapnik V.Statistical Learning Theory［M］.New York：Wiley，1998.