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

(1.中車株洲電機有限公司，湖南 株洲 412000；2.西安交通大學(xué) 現(xiàn)代設(shè)計及轉(zhuǎn)子軸承系統(tǒng)教育部重點實驗室,西安 710049)

1 概述

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

圖1 直驅(qū)風(fēng)力發(fā)電機組傳動鏈示意圖Fig.1 Transmission chain diagram of direct-driven wind turbine

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

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

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

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

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

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

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

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

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

Q=kδ1.11，

(1)

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

剛度系數(shù)的修正公式為[9]

(2)

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

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

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

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

δrj=δrcosψj，

(3)

δaj=δa，

(4)

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

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

(5)

式中：Dpw為滾子組節(jié)圓直徑。

滾子總的軸向位移為

(6)

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

δnj=δajcosα+δbajsinα=

(7)

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

(8)

(9)

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

(10)

由于滾子發(fā)生偏轉(zhuǎn)，切片產(chǎn)生的附加力矩為

(11)

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

Mji=Meji+Mθji=

(12)

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

根據(jù)載荷及力矩平衡關(guān)系可得

(13)

(14)

(15)

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

F=Kδ，

(16)

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

(17)

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

Ks=

(18)

3 軸承載荷分布數(shù)值求解

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

K=

(19)

以此建立系統(tǒng)力學(xué)平衡方程

P=Kδ,

(20)

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

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

4 模型驗證

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

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

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

4.1 載荷分布驗證

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

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

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

4.2 摩擦力矩驗證

軸承摩擦力矩測量示意圖如圖5所示，發(fā)電機豎直安置于平臺上，尼龍吊帶一端固定在發(fā)電機定軸上，另一端固定在叉車上，叉車拖動定軸旋轉(zhuǎn)，測量拉力值，換算得到主軸承摩擦力矩。裝配過程對2TS軸承端蓋預(yù)緊，預(yù)緊力取650 kN。測量裝配完成后的6臺風(fēng)力發(fā)電機軸承摩擦力矩(每臺測量10次，取平均值)，結(jié)果見表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范圍內(nèi)，最大值與最小值之間的相對誤差為10.8%。摩擦力矩測量值波動是軸承潤滑狀態(tài)、密封狀態(tài)、旋轉(zhuǎn)精度等因素共同影響的結(jié)果，此外，使用叉車測量大尺寸電機軸承的摩擦力矩，在實際操作上也會帶來一定誤差。

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

(21)

式中：μ為摩擦因數(shù)；d為摩擦力作用內(nèi)徑；P為軸承當(dāng)量動載荷。

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

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

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

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

5 結(jié)束語

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