基于Rayleigh法的風機塔筒預應力基本頻率解析計算方法

雷振博， 劉 綱,2， 楊 微， 李 楊

(1.重慶大學 土木工程學院，重慶 400045; 2. 重慶大學 山地城鎮建設與新技術教育部重點實驗室，重慶 400045;3. 中國船舶重工集團海裝風電股份有限公司，重慶 401122; 4. 重慶大學 機械工程學院，重慶 400044)

風電作為應對生態環境污染和氣候變化的主力軍，逐漸得到世界各國的普遍共識和大力發展。統計表明，截至2019年我國風電裝機容量已達2.21億kW，雄踞世界第一[1]。到2050年，我國風電裝機容量更有望達到30億kW，屆時將滿足國內30%以上的電力需求[2]。風力發電技術將成為我國經濟向綠色發展轉型、實現碳中和的重要手段。

設計過程中為避免出現共振且同時考慮經濟安全等因素，塔筒第一階頻率(基頻)往往會遠離風輪旋轉頻率(1f)、過槳頻率(3f)以及主要的環境荷載頻率[3]，故對風力機自振頻率的計算精度要求較高。而主流大型風力機塔筒通常采用變截面錐形筒體形式[4]，截面剛度隨高度呈非線性變化，這給塔筒基頻的準確計算帶來較大困難。

目前風力機塔筒自振頻率計算主要有解析和數值計算兩類方法。在解析計算方面，Byrne[5]將塔筒結構等效為懸臂梁，將機艙和風輪等效為頂部集中質量，基于結構動力學原理，獲得塔筒等效單自由度體系基頻計算公式。牛文杰[6]考慮樁-彈性地基作用影響，采用柔度法對風機結構進行多自由度動力分析，通過求解微分方程確定風力機塔筒基頻。在數值計算方面，楊春寶等[7]考慮風力機塔筒變截面特性，建立了近海風力機塔筒自振頻率求解數值計算方法，并利用實際工程進行有效性驗證。曾夢偉等[8]將風輪與機艙等效為質量點，建立機艙-塔筒-法蘭剛柔多體耦合精細化有限元模型，采用數值模擬方法計算塔筒結構自振頻率。田英鵬等[9]將風力機法蘭和螺栓連接視為剛性連接，將風輪、機艙及法蘭等效為質量點，采用殼單元構建塔筒結構有限元模型進行基頻計算，并通過試驗驗證了計算精度。

總體而言，采用精細化有限元模型能夠較為準確獲得風力機塔筒的基頻，但建模工作量極大，并不適合于塔筒前期設計的簡化、快速計算。而目前解析方法較為簡略，未考慮塔筒截面變剛度特性、機艙和風輪載荷引起的塔筒預變形，計算精度不高，難以滿足工程實際需求。針對以上問題，根據風機塔筒結構構造形式及力學特征，計入機艙及風輪質量對塔筒的預應力作用，提出基于Rayleigh法的風機塔筒預應力基頻解析計算方法。同時，考慮塔筒變形形狀、頂部集中質量偏心等多種情況，構建塔筒基頻工程實用解析計算方法，為風力機塔筒設計及共振控制提供方法支撐。

1 能量法基頻計算原理

1.1 風力機塔筒模型簡化

錐臺型塔筒由數段鋼制錐筒通過法蘭盤連接而成，每段錐筒的直徑、壁厚均不一致，由底向上逐漸減小,如圖1(a)所示。大量研究、試驗和現場實測表明，塔筒主要沿順風向和橫風向振動，且往往以順風向為主，故可將空間塔筒簡化為二維模型，取塔筒底面中心點為坐標原點，沿豎直方向為z軸，沿順風向為x軸，如圖1(b)所示。

圖1 風機等效單自由度結構Fig.1 Equivalent single degree of freedom flexible structure

由圖1可知，塔筒由n段鋼制塔段及n-1個法蘭盤組成，塔段從下至上的高度分別為l1，l2，…，ln，塔筒總高為L。考慮塔筒沿高度方向具有變化的彎曲剛度EI(z)及質量密度m(z)。法蘭盤質量依下至上為m1，m2，…，mn-1，塔筒頂部機組及風輪的質量記為mn，其對塔筒的預應力為集中質量自質量產生的軸向壓力N=mng。

1.2 能量法自振頻率計算原理

為方便采用能量法建立風力機塔筒運動方程，可將塔筒近似等效為單自由度體系[10]。根據結構動力學原理，假設塔筒某階振動的形狀函數為ψ(z)，則高度z處塔筒在t時刻的位移v(z,t)為

v(z,t)=ψ(z)Z(t)

(1)

式中：Z(t)為廣義坐標，表示t時刻塔頂處x向的位移大小。

根據虛功原理可得塔筒廣義單自由度體系的運動方程為

(2)

其中，

(3)

式中：m*為廣義質量；c*為廣義阻尼；k*為廣義彎曲剛度。

結合Rayleigh能量法和式(3)，可推導得到塔筒廣義單自由度體系的基頻為[11]

(4)

2 考慮預應力的基頻解析方法

2.1 塔筒預變形

大量文獻表明，塔筒在頂部機艙、風輪載荷作用下會產生壓應力和預變形，塔筒中的壓應力會導致塔筒剛度變化，進而影響塔筒頻率[12]。因此，機艙、風輪載荷引起壓預應力效應是影響風機塔筒基頻的重要因素之一。

假定塔筒在機艙、風輪載荷下產生的預變形為e(t)(見圖1(b)),根據幾何關系有

(5)

故考慮預應力的塔筒基頻為

(6)

塔筒的廣義單自由度體系的運動方程為

(7)

(8)

式中，δe(t)為對預變形為e(t)的變分。

對比式(4)和式(6)可知，塔筒機艙、風輪載荷作用下產生的壓應力將降低塔筒基頻。

2.2 平均等效當量慣性矩及壁厚

塔筒沿高度方向具有變化的壁厚及截面直徑，故簡化基頻計算時，采用等效壁厚t*來表示塔筒壁厚的影響[13]，其表達式為

(9)

式中：li為第i段塔筒高度；ti為第i段塔筒壁厚。

針對截面變直徑引起塔筒沿高度剛度變化問題，基于位移等效原理，采用平均等效當量慣性矩I*近似計算[14]

(10)

式中：i1為系數，當考慮塔筒頂部橫向集中荷載作用時i1=3，當考慮塔筒橫向均布荷載作用時i1=4；Ii為通過積分獲得的第i段塔筒的慣性矩，其表達式為

(11)

式中：Di為第i段塔筒的底部中直徑；di為第i段塔筒的頂部中直徑。

3 多參數影響下基本頻率分析

3.1 考慮塔筒變形形狀影響

由Rayleigh法基本原理可知，基頻計算精度依賴于所取形狀函數與真實振型的相似度，對典型懸臂結構，通常以集中荷載或均布荷載變形下的撓曲線作為形狀函數。風力機塔筒結構的力學模型，如圖2所示。

圖2 風機變形影響因素Fig.2 Influence factors of fan deformation

3.1.1 塔頂集中荷載

如圖2(a)所示。根據材料力學可知，懸臂柱自由端作用水平集中荷載產生的撓曲線為[15]

(12)

將式(12)代入式(1)可得塔頂集中荷載下塔筒變形的形狀函數為

(13)

將式(13)代入式(6)可得集中荷載作用下，考慮預壓應力影響的塔筒基頻表達式為

(14)

式中：ρ為材料密度；D為整個風力機塔筒底部中直徑；d為頂部中直徑。

3.1.2 均布荷載

由材料力學可知，懸臂柱橫向作用水平均布荷載產生的撓曲線為

(15)

將式(15)代入式(1)可得均布荷載下塔筒變形的形狀函數為

(16)

將式(16)代入式(6)可得均布荷載作用下，考慮預壓應力影響的塔筒基頻表達式為

(17)

3.2 考慮頂部質量偏心影響

由于塔筒頂部的機艙、葉輪在x向有偏心，即圖2(b)中集中力N與塔筒的中心軸有偏心距Δ。塔頂質量偏心將引起附加彎矩，則塔筒廣義幾何剛度應修正為

(18)

將式(13)、式(18)代入式(6)，式(16)、式(18)代入式(6)可分別得到集中荷載作用下考慮偏心的塔筒基頻f3、均布荷載作用下考慮偏心的塔筒基頻f4為

(19)

4 數值算例

4.1 3.4 MW-100塔型

以兩個風機塔筒為例，通過有限元數值計算結果驗證所提塔筒基頻計算公式的適用性。塔筒1為某3.4 MW風機100 m高塔筒，該型風機的機艙質量為126.75 t，風輪質量為96.5 t；風機塔筒整體呈圓臺狀，塔筒計算總高度為96.7 m；塔筒分節及截面詳細尺寸，如表1所示。塔筒材料采用Q345鋼，彈性模量取為2.06×105MPa，密度取為7 850 kg/m3。塔筒頂部風輪沿x方向偏心-4.475 m，機艙沿x方向偏心0.953 m，沿y向偏心-0.038 m。

表1 塔筒幾何參數Tab.1 Tower geometry parameters

通過ANSYS軟件建立該塔筒的有限元殼模型，并計入塔頂機艙、風輪質量偏心的影響。塔筒選用SHELL63單元模擬，并考慮實際塔筒的漸變分段效果。頂部機艙和風輪以及各法蘭盤等效為點質量，采用MASS21單元模擬。塔筒底部與基礎固定連接，建立的有限元模型如圖3所示。

圖3 風機塔筒有限元模型示意圖Fig.3 Finite element model of wind turbine tower

預應力模態分析通過有預應力的模態分析方法來實現，有預應力結構的固有頻率和模態計算，需要先進行結構的靜力分析獲得應力剛度矩陣，將其代入頻率分析方程來得到考慮機艙和風輪質量預應力的模態分析[16]。通過模態分析，獲取塔筒x向和y向的基頻如圖4所示。

從圖4可知，塔筒在x向和y向第一階頻率分別為0.222 Hz和0.225 Hz。根據式(14)、式(18)和式(19)，計算得到不同條件下的風機塔筒基頻，如表2所示。

圖4 100 m風機塔筒預應力模態示意圖Fig.4 Pre-stress mode diagram of 100 m wind turbine tower

表2 100 m塔筒基頻計算結果對比Tab.2 Frequency calculation of 100 m wind turbine

4.2 3.4 MW-90塔型

塔筒2為某3.4 MW風機90 m高塔筒。該型風機的機艙質量為126.75 t，風輪質量為101.72 t；風機塔筒整體呈圓臺狀，塔筒計算總高度為86.7 m；塔筒頂部風輪沿x方向偏心-4.475 m，機艙沿x方向偏心0.953 m,沿y向偏心-0.038 m。同4.1節相同建模方法，采用ANSYS有限元軟件計算得到塔筒在x向和y向基頻分別為0.262和0.265。根據式(14)、式(18)和式(19)，計算得到不同條件下的風機塔筒基頻如表3所示。

表3 90 m塔筒頻率計算結果對比Tab.3 Frequency calculation of 90 m wind turbine

從表2、表3計算結果可知，無論x向還是y向，采用Rayleigh能量法所得基頻均高于精細化有限元模型計算值，這主要因為能量法計算所假設的變形公式式(13)、式(16)均無法與真實變形曲線完全吻合所致；相比而言，f2，f4較f1，f3更加接近于有限元數值解，因此，宜采用均布荷載作用下的形狀函數進行解析法計算。同時，因x向偏心較大，考慮頂部質量偏心對x向基頻計算精度有顯著影響，在集中荷載變形下，塔筒1和塔筒2的精度分別可提高2.8%和1.9%，均布荷載變形下分別可提高3.4%和2.6%。故當頂部質量偏心較大時，應考慮偏心對解析法求解的影響。

5 試驗驗證

5.1 模型設計

以塔筒1為原型在實驗室搭建試驗塔筒，模型幾何相似比Sl為1/20。由于原型筒壁最薄處僅14 mm，無法完全遵守幾何相似關系進行縮尺，故主要保證模型與原型動力特性相似，即保證構件截面彎曲剛度相似與質量分布相似[17-18]。通過材性試驗測得其彈性模量為1.525×105MPa，為保證加工精度可控，模型塔筒每節截面和壁厚不變，按原型塔筒相應節底截面進行相似設計，試驗塔筒模型幾何參數如表4所示。

表4 模型幾何參數Tab.4 Model geometric parameters

與原型塔筒對應，模型塔筒總共分5段，段與段之間由法蘭連接，法蘭盤上設置配重籃，塔筒頂部設置配重箱考慮機艙和風輪質量。模型由Q345熱軋鋼板焊接而成組成，底部法蘭連接600 mm×600 mm×20 mm的矩形連接板，并通過18顆螺栓和鋼板連接。試驗模型頂部配重箱沿x方向偏心6 cm。筒體與法蘭之間采用對接焊縫連接，各法蘭間采用螺栓連接，如圖5所示。

5.2 測點及加載方案設計

在各配重籃和頂部配重箱的x和y向均布設P15H-2壓電加速度傳感器(見圖5)。

圖5 試驗模型Fig.5 Experiment model

5.3 頻率測試

在塔頂施加錘擊激勵使其產生自由振動，測得塔頂加速度時程如圖6所示。通過傅里葉變換，得到塔筒兩個方向的功率譜如圖7所示。

圖6 塔頂加速度響應時程曲線Fig.6 Acceleration response at the top of tower

圖7 加速度功率譜Fig.7 Acceleration power spectrum

從圖7可看出，模型對應x向和y向基頻分別為1.18 Hz和1.22 Hz。同時，試驗中僅采用頂部激勵時，圖6的加速度響應出現典型的拍頻振動特征及正交耦合現象，表明試驗模型的振動符合高柔結構的響應特征[19-21]。加速度功率譜僅有一個峰值，表明塔筒以第一階振動為主，這與大多數文獻及試驗結果相吻合[22-24]。

5.4 有限元數值模型

與4.1節相同建模方法相同，彈性模量采用材性試驗測得的1.525×105MPa。結合表4模型幾何參數，采用ANSYS有限元軟件計算試驗縮尺模型的x向和y向基頻分別為1.203 Hz和1.212 Hz。x向和y向的基本頻率與試驗測試的基頻誤差均在3%以內，與試驗測試值吻合較好。

5.5 解析法對比驗證

結合傳統風力機塔筒頻率解析方法及式(18)，計算可得試驗模型解析法的理論值與動力測試基頻，如表5所示。

表5 理論值與實測值對比Tab.5 Comparison between calculation and experiment

從表5可知，試驗塔筒在x向實測基頻低于y向基頻，這主要是因為該塔筒在x向存在偏心，將減少塔筒在該向剛度，故實測頻率值偏小，故在計算塔筒基頻時有必要計入偏心影響。就傳統解析法和本文所提方法而言，后者較前者更加精確，即x向計算精度提高4.52%，y向計算精度提高4.50%。這主要是因為本文方法綜合考慮了機艙、風輪質量和偏心，以及塔筒不同變形形狀等情況的影響。

6 結 論

針對風力機塔筒基頻解析計算精度難以滿足工程要求的問題，提出一種基于Rayleigh法的風力機塔筒基頻快速估算方法。考慮機艙、風輪質量對塔筒壓力作用，以及塔筒不同變形形狀函數、頂部集中質量偏心等系列實際情況，基于理論推導給出了解析計算公式，并通過數值算例與縮尺模型試驗驗證了所提方法的適用性及有效性。所得結論如下：

(1) 有限元模型與試驗結果均表明采用均布荷載作用下的變形形狀曲線將得到更高的解析計算精度，兩個數值算例中，f2和f4與有限元模型計算值的最大相對誤差僅為3.47%。因此，宜采用均布荷載作用下的撓曲線作為Rayleigh法的形狀函數。

(2) 頂部質量偏心對塔筒基頻計算精度有顯著影響，特別是在偏心距較大時。采用集中荷載作用下的變形曲線時，考慮偏心情況下塔筒1和塔筒2的精度分別可提高2.8%和1.9%；采用均布荷載作用下變形曲線時，考慮偏心情況下塔筒1和塔筒2的計算精度分別提高3.4%和2.6%。因此，在計算風機塔筒基頻時，應考慮機艙、風輪質量偏心的影響。

(3) 縮尺模型試驗表明，本文方法考慮了機艙、風輪質量及偏心，以及塔筒不同變形形狀等多種情況的綜合影響，x向基頻計算精度可提高4.52%，y向基頻計算精度可提高4.50%，較傳統解析方法更加精確。