多工況下兆瓦級海上風電齒輪箱均載性能優(yōu)化設(shè)計

陳巖松，朱才朝，譚建軍，宋朝省，宋海藍

(重慶大學 機械傳動國家重點實驗室，重慶 400044)

為了應(yīng)對日益嚴重的世界能源危機，風能作為一種最具發(fā)展?jié)摿Φ木G色可再生能源，引起了世界各國廣泛關(guān)注。風電齒輪箱是支撐風電機組高效開發(fā)風能資源的關(guān)鍵傳動部件，但在時變風速作用下，百米級大尺寸、大柔性葉片與塔筒的耦合變形會使風電齒輪箱的輸入載荷產(chǎn)生明顯的隨機特性，容易惡化齒輪箱行星級均載性能[1-2]，增大故障失效率[3]，因此開展多工況下風電齒輪箱行星級均載性能優(yōu)化研究對提高風電機組運行安全具有重要的意義。

1 風電齒輪箱系統(tǒng)動力學建模

1.1 風電齒輪箱結(jié)構(gòu)及傳動原理

圖1 風電機組結(jié)構(gòu)圖Fig. 1 Structure of wind turbine

表1 風電齒輪箱傳動結(jié)構(gòu)參數(shù)

1.2 風電齒輪箱系統(tǒng)動力學模型

1.2.1 外部激勵

風速頻率表示在風電機組20年運行期間內(nèi)，相同平均風速總時長與測量時長的百分比值，是反映風場平均風速統(tǒng)計特性的重要參數(shù)。根據(jù)某風場風速統(tǒng)計數(shù)據(jù)[15]，采用雙參數(shù)威布爾分布對風速頻率進行擬合，其概率密度函數(shù)為

(1)

式中，k=8.426，c=1.708。

對應(yīng)的分布函數(shù)為

(2)

式中：u為風速；k為形狀參數(shù)；c為尺度參數(shù)。

考慮風切變效應(yīng)，輪轂高度處的風速

(3)

式中：zhub為輪轂處高度；uhub為輪轂處風速；z為欲求的風速離地面高度。

采用Kaimal模型計算湍流風，湍流風在空間3個分量的譜為

(4)

(5)

式中，Λu為湍流尺度系數(shù)。

空間中k,h2點處風速相關(guān)模型表示為

(6)

式中：PCoh(f)為空間相干度；Sk,k(f)和Sh,h(f)分別為k和h處的功率譜；Sk,h(f)為空間k和h處的互功率譜；f為頻率。

1.2.2 內(nèi)部激勵

齒輪嚙合是齒輪箱內(nèi)部激勵的主要來源之一。如圖2所示，筆者采用齒輪分片法，將齒輪沿齒寬方向劃分為若干薄片，用彈簧單元表示一對薄片齒的嚙合，其嚙合剛度通過ISO 6336標準[17]計算。齒輪副法向嚙合剛度函數(shù)為拋物線模型，最大值在節(jié)點處，計算公式為

圖2 輪齒分片示意圖Fig. 2 Schematic diagram of gear slicing

cmax=c′CR，

(7)

式中：c′為單齒嚙合剛度；CR為輪齒結(jié)構(gòu)系數(shù)。

單齒嚙合剛度的計算公式為

c′=cthCMCBcosβ，

(8)

式中：cth為理論單齒嚙合剛度；CM為理論修正系數(shù)；CB為基本齒條系數(shù)；β為螺旋角。

理論單齒嚙合剛度cth定義為

(9)

(10)

式中：zn1和zn2分別為主動輪和從動輪的當量齒數(shù)；x1和x2分別為主動輪和從動輪的變位系數(shù)。

齒輪時變嚙合剛度函數(shù)可以表示為

(11)

傳動誤差計算公式為

(12)

式中：φ10和φ20分別為主動輪和從動輪進入嚙合時刻的初始轉(zhuǎn)角；φ1和φ2分別為主動輪和從動輪轉(zhuǎn)動的位移角；N1和N2分別為主動輪和從動輪的齒數(shù)。

1.2.3 齒輪箱系統(tǒng)動力學模型

采用SIMPACK建立風電齒輪箱系統(tǒng)動力學模型如圖3所示。采用模態(tài)縮減法計算部件柔性變形[18-19]，其中，柔性體節(jié)點平移速度為

圖3 SIMPACK齒輪箱系統(tǒng)動力學模型Fig. 3 Gearbox dynamic model of SIMPACK

(13)

式中：r為柔性體節(jié)點位移；GAB為局部參考到慣性參考的轉(zhuǎn)換矩陣；s和u分別為柔性體變形前和變形后的節(jié)點位移；I為節(jié)點慣性張量；Φ為節(jié)點模態(tài)矩陣；ξ為節(jié)點廣義坐標矩陣。

柔性體節(jié)點角速度為

GωJ=GωB+GωP，

(14)

式中：GωB為物體的角速度；GωP為由于物體變形引起的角速度。

基于公式(13)和式(14)，可得柔性體動能表達式為

(15)

式中，ρ,V,m和m分別為部件的密度、體積、節(jié)點質(zhì)量和質(zhì)量矩陣。

柔性體勢能包括重力和彈性勢能，其表達式為

(16)

式中：x為部件參考系相對于慣性參考系的位移；A為轉(zhuǎn)換矩陣；g為重力加速度；q為模態(tài)坐標系；K為廣義剛度矩陣。

基于式(15)和式(16)，可得風電機組齒輪箱系統(tǒng)動力學方程為

(17)

式中：D為模態(tài)阻尼矩陣；fg為廣義重力；ψ為代數(shù)約束方程；λ為拉格朗日算子；Q為廣義力。

2 風電齒輪箱均載優(yōu)化

2.1 設(shè)計變量

齒輪修形一般包括齒廓修形和齒向修形。如圖4(a)所示，針對齒廓修形，采用齒頂修緣，包括最大修形量e、修形長度λ和修形曲線3個要素。其中，修形曲線采用拋物線。如圖4(b)所示，齒向修形則采用修鼓形方式。由于風電齒輪箱第一級行星級容易受到非扭載荷的作用，造成其動態(tài)性能較差[20]，因此重點優(yōu)化第一級行星級齒輪修形參數(shù)，共計9個優(yōu)化變量，如表2所示。

圖4 齒輪修形示意圖Fig. 4 Schematic diagram of gear modification

表2 優(yōu)化變量參數(shù)表

2.2 約束條件

基于式(18)～式(19)計算齒廓修形參數(shù)[21]

(18)

式中：最大修形量ek=fKT+fm；fKT為彈性變形，加工誤差fm=fpb+1/3ff，fpb為基節(jié)誤差，ff為齒形誤差。l為沿單對齒嚙合的上界點至嚙合始點長度。x為嚙合位置坐標(原點在界點處)。

λ=Pb(εα-1)，

(19)

式中：λ為修形長度；Pb為齒輪基節(jié)；εα為齒輪端面重合度。

基于式(20)計算齒向修形參數(shù)[21]

(20)

式中：δ為鼓形量；Ft為嚙合圓周力；b為齒寬。

分別選取式(18)～式(20)計算值的70%和130%作為表2中修形參數(shù)的最小值和最大值。

2.3 優(yōu)化目標

如式(21)所示，對n種風速工況下行星級均載系數(shù)求和，并將其最小值作為優(yōu)化目標。

minf∑=k1f1+k2f2+…+knfn，

(21)

式中，ki(i∈[1,n])為第i種工況對應(yīng)的權(quán)重值，基于式(22)計算得到。fi(i∈[1,n])為第i種工況下行星級均載系數(shù)。

(22)

式中，f(u)為風速概率密度函數(shù)，可由公式(1)計算得到。ai和bi分別為該風速工況下的風速區(qū)間的下界和上界。

采用拉丁超立方抽樣方法，對表2中9個優(yōu)化變量進行抽樣組合，然后將每一種抽樣組合作為齒輪箱系統(tǒng)動力學模型(式(17))中第一級行星級齒輪修形參數(shù)的輸入，并計算行星級均載系數(shù)[22-23]。最后，采用支持向量機回歸(SVR)[24]重構(gòu)第i種風速工況下齒輪修形參數(shù)與均載系數(shù)之間的映射函數(shù)，即可得式(21)中fi，其中fi(i∈[1,n])。

假設(shè)x是第i種風速工況下的9個齒輪修形參數(shù)，f(x)則是對應(yīng)的均載系數(shù)，利用SVR可以構(gòu)建x和f(x)之間的非線性關(guān)系為

f(x)=wTφ(x)+b，

(23)

式中：f(x)表示預(yù)測值；φ(x)表示非線性映射函數(shù)；w和b是系數(shù)[25]。

圖5所示為在多風速工況下風電齒輪箱行星級均載特性優(yōu)化流程。

圖5 風電機組多工況優(yōu)化流程Fig. 5 Flow chart of multi-operating condition optimization of wind turbine

首先，根據(jù)風速頻率曲線(式(1))生成具有n種風速工況，并作為風電機組整機模型的外部激勵，進而計算主軸與輪轂連接處六自由度載荷(Fx，F(xiàn)y，F(xiàn)z，Mx，My，Mz)；然后，將每一種風速工況下計算得到的主軸與輪轂連接處六自由度載荷作為齒輪箱系統(tǒng)動力學模型的載荷輸入，計算不同齒輪修形參數(shù)組合(表2)下的行星級均載系數(shù)，并利用SVR構(gòu)建修形參數(shù)與均載系數(shù)之間的映射關(guān)系；最后，將n種風速工況下行星級均載系數(shù)進行加權(quán)求和，以其最小值作為優(yōu)化目標，采用遺傳算法求解最優(yōu)的齒輪修形參數(shù)。

在一臺CPU主頻為2.9G Hz的計算機上使用SVR進行齒輪修形參數(shù)尋優(yōu)所需要的計算時間約10 min。而采用SIMPACK動力學模型進行一次計算需耗時20 min。由此可見，文中提出的優(yōu)化方法可以大幅提高優(yōu)化效率。

3 結(jié)果與討論

3.1 外部激勵分析

如圖6(a)所示，從風電機組切入至切出風速區(qū)間5～23 m/s，間隔3 m/s，共計取n=7組風速工況，并利用式(22)計算對應(yīng)的風速頻率，作為式(21)中工況的權(quán)重值。參照IEC 61400-1標準，湍流強度設(shè)置為0.14[16]，可得額定風速(11 m/s)下u、v和w3個方向的湍流風速，如圖6(b)所示。

圖6 風速頻率及時域風速Fig. 6 Wind speed frequency and time domain wind speed

圖7所示為7種風速工況下主軸與輪轂連接處的六自由度載荷。從圖中可以看出，在額定風速以下，隨著風速的增加，氣動轉(zhuǎn)矩也逐漸增大，但非扭載荷與氣動轉(zhuǎn)矩的比值顯著下降；在額定風速以上時，由于變槳系統(tǒng)開始運行，此時主軸處氣動轉(zhuǎn)矩近似恒定，而非扭載荷與氣動轉(zhuǎn)矩的比值逐漸增大。

圖7 不同風速下主軸載荷Fig. 7 Mainshaft load under different wind speeds

3.2 內(nèi)部激勵分析

筆者僅針對第一級行星級齒輪開展修形優(yōu)化設(shè)計，圖8所示為非扭載荷對齒輪箱行星級齒輪副傳動誤差的影響。從圖中可以看出，隨著風速的增加，傳動誤差幅值也逐漸增大；當考慮非扭載荷時，傳動誤差幅值略有增大。

圖8 非扭載荷對傳動誤差影響Fig. 8 Influence of non-torque load on transmission error

圖9 不同風速下優(yōu)化前后行星級傳動誤差Fig. 9 Planetary gear transmission error before and after optimization under different wind speeds

圖10所示為額定風速工況下行星級齒輪副時變嚙合剛度優(yōu)化前后對比。從圖中可以看出，優(yōu)化后的齒輪副嚙合剛度激勵(波動)顯著減小。圖11所示為不同風速工況下行星級齒輪副嚙合剛度激勵優(yōu)化前后對比。從圖中可以看出，齒輪副嚙合剛度激勵受風速工況的影響不大，整體變化較為平緩。但優(yōu)化后的行星級嚙合剛度激勵顯著降低。

圖10 行星級優(yōu)化前后齒輪副時變嚙合剛度對比Fig. 10 Comparison of time-varying meshing stiffness of planetary gears before and after optimization

圖11 不同風速下行星級齒輪副嚙合剛度激勵Fig. 11 Meshing stiffness of planetary gears under different wind speeds

3.3 齒面載荷分析

圖12 額定風速下優(yōu)化前后齒面載荷分布Fig. 12 Load distribution of tooth surface before and after optimization under rated wind speed

3.4 均載性能分析

通過行星輪和內(nèi)齒圈動態(tài)嚙合力計算均載系數(shù)k1，同時利用行星輪軸承力計算均載系數(shù)k2[24-25]，分別如式(24)和式(25)所示。

(24)

(25)

圖13 額定風速下行星輪-內(nèi)齒圈動態(tài)嚙合力及均載系數(shù)k1Fig. 13 Dynamic meshing force and load sharing coefficient of planet-ring gear under rated wind speed

圖14 額定風速下行星輪軸承力及均載系數(shù)k2Fig. 14 Dynamic bearing force and load sharing coefficient of planet gear under rated wind speed

利用式(26)計算行星輪上、下風向軸承力差值與軸承力額定值的百分比p。

(26)

式中：Fupwind和Fdownwind分別為行星輪上、下風向軸承力；Frated為額定的行星輪軸承力。

圖15所示為第1個行星輪上、下風向軸承力差值。從圖中可以看出，當不考慮非扭載荷時，行星輪1上、下風向軸承力差值較小，而考慮非扭載荷后，行星輪1上、下風向軸承力差值顯著增大，但隨著風速的增加，其差值又逐漸減小。其主要原因是在低風速工況時，氣動轉(zhuǎn)矩較小，行星輪軸承力也較小，此時行星輪軸承力易受到非扭載荷的影響[2,19]，而在額定風速及以上時，氣動轉(zhuǎn)矩為額定值，此時行星輪軸承以承擔扭矩載荷為主。優(yōu)化后的行星輪上、下風向軸承力差值明顯減小，偏載現(xiàn)象得到明顯改善，與圖12結(jié)果相符。

圖15 優(yōu)化前后行星輪1上下風向軸承力Fig. 15 Upwind and downwind bearing force of planet gear 1 before and after optimization

圖16所示為優(yōu)化前后的行星級均載系數(shù)k1和k2對比。從圖中可以看出，不考慮非扭載荷時，行星級均載系數(shù)隨著風速的增加而增大，在達到額定風速后保持不變，但考慮非扭載荷后，隨著風速的增大，行星級均載系數(shù)先減小再增大。結(jié)合圖7可知，在額定風速工況以下時，隨著風速的增加，氣動轉(zhuǎn)矩也逐漸增大，但非扭載荷與氣動轉(zhuǎn)矩的比值顯著下降，因此行星級均載系數(shù)逐漸減小，均載性能提高。當達到額定風速后，主軸處氣動轉(zhuǎn)矩近似恒定，但非扭載荷與氣動轉(zhuǎn)矩的比值逐漸增大，導(dǎo)致行星級均載系數(shù)增大，均載性能降低。優(yōu)化后的行星級均載性能在不同風速工況下均有所提升，在額定風速及以上時均載優(yōu)化效果較好。

圖16 優(yōu)化前后行星級均載系數(shù)對比Fig. 16 Comparison of load sharing coefficient of planetary stage before and after optimization

4 結(jié) 論

該文考慮不同風速工況下風電齒輪箱時變輸入載荷的影響，結(jié)合風電齒輪箱系統(tǒng)動力學模型，構(gòu)建了多風速工況下風電齒輪箱行星級均載性能優(yōu)化模型，以齒輪箱行星級均載系數(shù)綜合最小為優(yōu)化目標對行星級齒輪副進行修形優(yōu)化設(shè)計，對比了不同風速工況下優(yōu)化前后風電齒輪箱運行性能，得出以下結(jié)論：

1) 在額定風速以下時，隨著風速的增加，非扭載荷與氣動轉(zhuǎn)矩比值顯著減小，而在額定風速及以上時，在變槳控制作用下其比值逐漸增大。較大的非扭載荷占比會造成行星級齒輪和軸承產(chǎn)生明顯的偏載。

2) 不同風速工況下，優(yōu)化后的行星級齒輪副傳動誤差和嚙合剛度波動幅值均顯著降低，并且齒面載荷分布均勻，行星輪上、下風向軸承的偏載現(xiàn)象得到有效改善。

3) 隨著風速的增加，非扭載荷占比變化會使行星級均載系數(shù)先減小再增大，并且優(yōu)化后的行星級均載性能在不同風速工況下均有所提升，在額定風速及以上時均載優(yōu)化效果較好。