直驅(qū)風(fēng)力發(fā)電機(jī)主軸承載荷分布數(shù)值分析

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

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

1 概述

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

圖1 直驅(qū)風(fēng)力發(fā)電機(jī)組傳動(dòng)鏈?zhǔn)疽鈭DFig.1 Transmission chain diagram of direct-driven wind turbine

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

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

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

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

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

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

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

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

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

Q=kδ1.11，

(1)

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

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

(2)

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

在徑向載荷Fr、軸向載荷Fa及力矩M聯(lián)合作用下，圓錐滾子軸承內(nèi)圈相對(duì)于外圈將產(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個(gè)滾子的徑向位移和軸向位移。

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

(5)

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

滾子總的軸向位移為

(6)

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

δnj=δajcosα+δbajsinα=

(7)

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

(8)

(9)

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

(10)

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

(11)

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

Mji=Meji+Mθji=

(12)

式中：xi為第i個(gè)切片中心的坐標(biāo)。

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

(13)

(14)

(15)

由(12)—(14)式可得Fr，F(xiàn)a，M與δr，δa，θ的關(guān)系，建立軸承單元?jiǎng)偠确匠?/p>

F=Kδ，

(16)

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

(17)

直驅(qū)風(fēng)力發(fā)電機(jī)主軸一般為變截面空心軸，異形輪轂與主軸相連，為方便編程將其等效為主軸的延伸段，軸單元采用歐拉-伯努利梁?jiǎn)卧?。編制程序時(shí)，在主軸上載荷(風(fēng)載、重力、單邊磁拉力)作用位置及軸承位置處設(shè)置節(jié)點(diǎn)，將軸分為5個(gè)軸段，每一段都是一個(gè)6×6的梁?jiǎn)卧?，其單元?jiǎng)偠染仃嚤磉_(dá)式見(jiàn)文獻(xiàn)[12]。主軸總體剛度矩陣由所有梁?jiǎn)卧獎(jiǎng)偠染仃嚱M合得到，其單元?jiǎng)偠菿s組合形式為

Ks=

(18)

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

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

K=

(19)

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

P=Kδ,

(20)

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

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

4 模型驗(yàn)證

以某型號(hào)直驅(qū)風(fēng)力發(fā)電機(jī)2TS主軸承為例，軸系模型如圖3所示。軸系載荷中心位于圖1中的輪轂中心處。為便于計(jì)算，將輪轂簡(jiǎn)化為空心主軸的延長(zhǎng)，載荷中心位于圖3中空心主軸最左端軸心O處，軸承外圈固定。軸承相關(guān)參數(shù)見(jiàn)表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àn)證

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

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

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

4.2 摩擦力矩驗(yàn)證

軸承摩擦力矩測(cè)量示意圖如圖5所示，發(fā)電機(jī)豎直安置于平臺(tái)上，尼龍吊帶一端固定在發(fā)電機(jī)定軸上，另一端固定在叉車(chē)上，叉車(chē)拖動(dòng)定軸旋轉(zhuǎn)，測(cè)量拉力值，換算得到主軸承摩擦力矩。裝配過(guò)程對(duì)2TS軸承端蓋預(yù)緊，預(yù)緊力取650 kN。測(cè)量裝配完成后的6臺(tái)風(fēng)力發(fā)電機(jī)軸承摩擦力矩(每臺(tái)測(cè)量10次，取平均值)，結(jié)果見(jiàn)表2。

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

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

由表2可知：摩擦力矩測(cè)量值在4 176～4 680 N·m范圍內(nèi)，最大值與最小值之間的相對(duì)誤差為10.8%。摩擦力矩測(cè)量值波動(dòng)是軸承潤(rùn)滑狀態(tài)、密封狀態(tài)、旋轉(zhuǎn)精度等因素共同影響的結(jié)果，此外，使用叉車(chē)測(cè)量大尺寸電機(jī)軸承的摩擦力矩，在實(shí)際操作上也會(huì)帶來(lái)一定誤差。

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

(21)

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

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

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

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

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

5 結(jié)束語(yǔ)

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