直驅風力發電機主軸承載荷分布數值分析

李海江，王騰飛，郭四洲，馬帥軍，閆柯

(1.中車株洲電機有限公司，湖南 株洲 412000；2.西安交通大學 現代設計及轉子軸承系統教育部重點實驗室,西安 710049)

1 概述

風能是重要的清潔能源，大規模風能開發與利用是我國重要的能源戰略。傳統非直驅風力發電機傳動鏈由葉片、輪轂、主軸、齒輪箱和發電機組成，齒輪箱故障會導致風電裝備發生安全事故，增加維護成本[1]。近年來，大功率(目前已高達12 MW)直驅風力發電機組已成為我國風電新增容量的主力機型之一[2]。一種典型的直驅風力發電機組傳動鏈示意圖如圖1所示，發電機組由風輪(葉片+輪轂)直接驅動發電機，發電機轉軸(主軸)一般采用空心軸，發電機軸承(直驅風力發電機組主軸承)為2套單列圓錐滾子軸承。直驅風力發電機組傳動鏈中沒有齒輪箱，從根本上解決了齒輪箱的安全、維護等重大問題，但強陣風直接沖擊在主軸承上，對主軸承性能提出了極高要求。主軸承一般在低速重載工況下工作，準確分析其載荷分布，合理設計和配置主軸承，對于直驅風力發電機組安全運行至關重要。

圖1 直驅風力發電機組傳動鏈示意圖Fig.1 Transmission chain diagram of direct-driven wind turbine

直驅風力發電機主軸承(以下簡稱主軸承)受力分析是一個復雜的非線性問題，其難點在于既要考慮軸承各個方向剛度之間的耦合性與非線性，又要與主軸的受力相結合進行軸系分析，目前一般采用有限元法[3-4]。采用有限元法分析復雜軸系問題時，軸承主要有4種處理方法：

1)忽略軸承剛度將其簡化成支點，或忽略軸承剛度的耦合性和非線性將其簡化成線性彈簧，然后進行軸系分析[5]。該方法過于簡化，難以反映軸承的真實受力狀態。

2)采用全實體有限元模型，通過有限元非線性接觸算法求解軸承的受力與變形[6-7]。該方法可以考慮軸承內部結構的影響，計算結果接近真實值，但全實體有限元模型網格數量大，對計算硬件資源要求極高，耗時耗力，且不易收斂。

3)將滾子與滾道的接觸簡化為非線性彈簧，通過有限元非線性接觸算法求解軸承的受力與變形[8]。該方法減少了計算量，但非線性彈簧建模繁瑣。

4)在有限元法基礎上，考慮軸承剛度耦合性和非線性，建立軸系剛度非線性有限元模型并進行數值求解[8-10]。該方法理論基本完備，模型簡化合理，在處理復雜軸系問題時具有更高的計算效率和計算精度。

現有文獻鮮有報道考慮軸系剛度非線性的直驅風力發電機主軸承載荷分布有限元分析方法。針對直驅風力發電機主軸承的結構特點，在第4種方法的基礎上，提出一種考慮軸承剛度耦合性和非線性的有限元法，推導軸承單元剛度矩陣的解析表達式，采用考慮剪切變形的歐拉-伯努利梁單元建立空心主軸模型，構建直驅風力發電機軸系整體計算模型，采用牛頓-拉夫遜法求解極限風載作用下主軸承的載荷分布,運用自主開發的程序對某非標2TS軸承計算驗證。載荷分布計算結果與成熟商業軟件結果進行對比，并將所提計算模型的摩擦力矩計算結果與試驗結果進行對比。

2 軸承單元與空心軸單元模型

當軸系承受載荷不同時，軸承中受載滾動體數量和所受載荷大小不同，軸承變形與載荷之間存在非線性關系，剛度表現為非線性。軸承各方向剛度之間相互影響，存在耦合關系。以直驅風力發電機中常用的2TS軸承為例推導考慮剛度非線性和耦合性的剛度矩陣表達式。

對于鋼制軸承，滾子與套圈滾道的接觸載荷為[11]

Q=kδ1.11，

(1)

式中：k為滾子與套圈滾道的接觸剛度系數；δ為滾子法向接觸變形。

剛度系數的修正公式為[9]

(2)

式中：Lwe為滾子有效長度；αe，αi，αf分別為圓錐滾子與外圈、內圈及擋邊的接觸角。

在徑向載荷Fr、軸向載荷Fa及力矩M聯合作用下，圓錐滾子軸承內圈相對于外圈將產生徑向位移δr、軸向位移δa和轉角θ，如圖2所示。

圖2 圓錐滾子軸承受力示意圖Fig.2 Stress diagram of tapered roller bearing

在位置角ψj處的滾子位移為

δrj=δrcosψj，

(3)

δaj=δa，

(4)

式中：δrj，δaj分別為第j個滾子的徑向位移和軸向位移。

轉角θ會引起位置角ψj處滾子的軸向位移和徑向位移變化，但不會改變總的法向接觸變形。由轉角θ引起位置角ψj處滾子的軸向位移為

(5)

式中：Dpw為滾子組節圓直徑。

滾子總的軸向位移為

(6)

轉角θ引起的徑向位移較小，可忽略?？紤]滾子偏轉對接觸載荷及變形的影響，將滾子沿長度方向切成m片(圖2)，可得滾子與滾道接觸法線法向變形為

δnj=δajcosα+δbajsinα=

(7)

第j個滾子第i個切片所受的徑向力、軸向力分別為

(8)

(9)

第j個滾子第i個切片在軸承中心形成的力矩為

(10)

由于滾子發生偏轉，切片產生的附加力矩為

(11)

則第j個滾子第i個切片形成的總力矩為

Mji=Meji+Mθji=

(12)

式中：xi為第i個切片中心的坐標。

根據載荷及力矩平衡關系可得

(13)

(14)

(15)

由(12)—(14)式可得Fr，Fa，M與δr，δa，θ的關系，建立軸承單元剛度方程

F=Kδ，

(16)

其中，位移矢量δ=(δr,δa,θ)T，載荷矢量F=(Fr,Fa,M)T。計算得到軸承剛度矩陣為

(17)

直驅風力發電機主軸一般為變截面空心軸，異形輪轂與主軸相連，為方便編程將其等效為主軸的延伸段，軸單元采用歐拉-伯努利梁單元建模。編制程序時，在主軸上載荷(風載、重力、單邊磁拉力)作用位置及軸承位置處設置節點，將軸分為5個軸段，每一段都是一個6×6的梁單元，其單元剛度矩陣表達式見文獻[12]。主軸總體剛度矩陣由所有梁單元剛度矩陣組合得到，其單元剛度Ks組合形式為

Ks=

(18)

3 軸承載荷分布數值求解

為計算軸承載荷分布，需要將軸和軸承的剛度矩陣進行組合形成總體剛度矩陣。最終軸承節點的位移和相應軸節點處的位移必須滿足變形協調關系，軸承節點產生的力和相應軸節點處產生的力要滿足平衡條件。假定安裝軸承的2個節點分別為n1和n2，軸承節點剛度矩陣為Kb，軸剛度矩陣為Ks，只需把共節點處的軸剛度矩陣和軸承剛度矩陣疊加即生成系統總體剛度矩陣K，可表示為

K=

(19)

以此建立系統力學平衡方程

P=Kδ,

(20)

式中：P為載荷矢量，即施加于軸系的實際載荷；δ為軸變形量。

可以采用牛頓-拉夫遜法求解該非線性方程組，收斂判據為相鄰兩次迭代的節點位移或節點力小于給定小量。計算得到的節點位移矢量即為軸上各節點的位移，將軸承所在節點的位移代入軸承載荷計算公式，即可得到各軸承的載荷。將上述過程編制成VB程序，可快速求解軸承載荷分布。

4 模型驗證

以某型號直驅風力發電機2TS主軸承為例，軸系模型如圖3所示。軸系載荷中心位于圖1中的輪轂中心處。為便于計算，將輪轂簡化為空心主軸的延長，載荷中心位于圖3中空心主軸最左端軸心O處，軸承外圈固定。軸承相關參數見表1。等效載荷Fx=100 kN，Fy=960 kN，Fz=110 kN，My=-10 000 kN·m,Mz=-900 kN·m。內外圈及滾子材料為滲碳鋼，其彈性模量為206 GPa，泊松比為0.3。

圖3 軸系模型示意圖Fig.3 Diagram of shafting model

表1 圓錐滾子軸承基本參數Tab.1 Basic parameters of tapered roller bearing

4.1 載荷分布驗證

為驗證計算模型的正確性，采用上述方法計算軸承載荷分布，并與成熟商業軟件計算結果進行對比。

理論計算及商業軟件基于各自的算法與程序，得到的載荷分布如圖4所示(滾子0°位置角位于軸承正下方)，軸承承載區基本一致，前軸承最大載荷偏差為5.6%，后軸承最大載荷偏差為1%，誤差在允許范圍內，說明了計算模型的正確性。將最大載荷代入自主開發的圓錐滾子接觸分析數值計算程序，即可求得滾子接觸應力分布。綜合滾子載荷分布與接觸應力分布，對照風電行業相關設計標準，可判斷主軸承是否能安全運行。

圖4 軸承載荷分布Fig.4 Load distribution of bearings

4.2 摩擦力矩驗證

軸承摩擦力矩測量示意圖如圖5所示，發電機豎直安置于平臺上，尼龍吊帶一端固定在發電機定軸上，另一端固定在叉車上，叉車拖動定軸旋轉，測量拉力值，換算得到主軸承摩擦力矩。裝配過程對2TS軸承端蓋預緊，預緊力取650 kN。測量裝配完成后的6臺風力發電機軸承摩擦力矩(每臺測量10次，取平均值)，結果見表2。

圖5 主軸承摩擦力矩測量示意圖Fig.5 Measuring diagram of friction torque of main bearing

表2 摩擦力矩測量值Tab.2 Measured values of friction torque

由表2可知：摩擦力矩測量值在4 176～4 680 N·m范圍內，最大值與最小值之間的相對誤差為10.8%。摩擦力矩測量值波動是軸承潤滑狀態、密封狀態、旋轉精度等因素共同影響的結果，此外，使用叉車測量大尺寸電機軸承的摩擦力矩，在實際操作上也會帶來一定誤差。

采用經驗公式計算2TS主軸承的理論摩擦力矩[13]

(21)

式中：μ為摩擦因數；d為摩擦力作用內徑；P為軸承當量動載荷。

用提出的模型計算得到試驗條件下主軸承的等效靜載荷見表3，參考ISO 281：2007 Rolling bearings—dynamic load ratings and rating life，可由等效靜載荷計算得到當量動載荷P，由(21)式計算出M。

表3 軸承載荷計算結果Tab.3 Load calculation results of bearings

(21)式中摩擦因數μ與軸承運行、承載、潤滑、密封狀態等因素有關，目前難以精確測定，文獻[14]給出了摩擦因數取值范圍為0.001 8～0.002 8，將其代入(21)式可得前后軸承摩擦力，進而得到總摩擦力矩M總=M前+M后=3 817～5 938 N·m。摩擦力矩實測值在4 176～4 680 N·m范圍內波動，包含在計算值3 817～5 938 N·m范圍內，說明了計算模型的正確性。

由于(21)式中摩擦因數μ取值波動較大，最大值與最小值之間的相對誤差達27.4%，也造成摩擦力矩計算值相對誤差較大，而在表2中，摩擦力矩測量值相對誤差僅10.8%。今后需進一步研究摩擦因數μ的取值范圍，從而更好地驗證計算模型的準確性。

5 結束語

針對直驅風力發電機主軸承建立了一種考慮軸承剛度非線性和耦合性的軸承單元，推導了軸承剛度矩陣的解析表達式。針對直驅風力發電機變截面空心主軸建立了梁單元主軸模型，將軸承與主軸模型組合得到直驅風力發電機軸系有限元計算模型，并采用牛頓-拉夫遜法進行求解。通過與商用軟件計算結果及現場測試結果對比，證實了提出的模型能夠有效計算給定外載荷下直驅風力發電機主軸承的載荷分布，已成功應用于直驅風力發電機的開發設計。