兆瓦級風(fēng)電制動器熱-結(jié)構(gòu)耦合仿真研究

沙智華，劉楠，劉宇，馬付建，尹劍，張生芳

(大連交通大學(xué) 機(jī)械工程學(xué)院，遼寧 大連 116028) *

?

兆瓦級風(fēng)電制動器熱-結(jié)構(gòu)耦合仿真研究

沙智華，劉楠，劉宇，馬付建，尹劍，張生芳

(大連交通大學(xué) 機(jī)械工程學(xué)院，遼寧 大連 116028)*

針對大兆瓦風(fēng)電制動器制動過程，考慮制動摩擦副作用區(qū)域?qū)挾燃捌溆绊懴碌木€速度徑向差異，提出速度梯度循環(huán)法，對制動過程摩擦副進(jìn)行熱-結(jié)構(gòu)耦合分析.基于ANSYS軟件建立制動器摩擦副三維瞬態(tài)傳熱有限元模型，運(yùn)用速度梯度循環(huán)法推導(dǎo)出熱流密度的加載式，用以計(jì)算制動過程中產(chǎn)生的摩擦熱流，對制動區(qū)域溫度場進(jìn)行數(shù)值模擬.從分析結(jié)果表明制動閘片摩擦區(qū)域溫度分布在制動盤徑向呈現(xiàn)弧度明顯的等溫分布，溫度梯度隨半徑增大而增大.以速度梯度循環(huán)法將熱分析結(jié)果代入結(jié)構(gòu)場對閘片摩擦區(qū)域受力及變形進(jìn)行耦合分析并預(yù)估其磨損狀況.通過與傳統(tǒng)均勻加載方法對比發(fā)現(xiàn)使用速度梯度循環(huán)法的分析結(jié)果與實(shí)際更為接近.所提出的分析方法為模擬大尺寸盤式制動器的摩擦制動過程提供參考.

風(fēng)電制動器; 熱-結(jié)構(gòu)耦合; 速度梯度循環(huán)法 ;有限元仿真

0 引言

風(fēng)力發(fā)電機(jī)組在超風(fēng)速、需要故障排除及日常維護(hù)等情況時需要停機(jī)制動，為保證在最短時間內(nèi)停機(jī)和在風(fēng)機(jī)維護(hù)時的操作安全，需要風(fēng)電制動器提供機(jī)械制動，因此風(fēng)電制動器是風(fēng)力發(fā)電必不可少的配套裝置.在制動器工作時，利用非旋轉(zhuǎn)元件與旋轉(zhuǎn)元件之間的相互摩擦來阻止轉(zhuǎn)動或轉(zhuǎn)動的趨勢，最終實(shí)現(xiàn)風(fēng)機(jī)系統(tǒng)的停機(jī).在這個過程中制動摩擦?xí)r，由于摩擦副接觸不均勻，實(shí)際接觸面積遠(yuǎn)小于名義接觸面積且接觸不連續(xù)造成接觸表面溫度分布不均勻，導(dǎo)致盤片發(fā)生熱變形，直接影響接觸狀態(tài)和接觸壓力，從而影響摩擦熱流輸入強(qiáng)度，這種耦合行為導(dǎo)致摩擦材料熱彈性不穩(wěn)定(Thermoelastic Instability, TEI)[1]，使得風(fēng)電制動器產(chǎn)生振動與噪聲[2].同時，熱應(yīng)力的作用易使摩擦材料熱衰退，產(chǎn)生初始裂紋并加劇摩擦副表面刮削現(xiàn)象[3].由此可見，在制動過程中，溫度場、應(yīng)力場對風(fēng)電制動器的制動效能有很大影響，二者的耦合分析是制動器設(shè)計(jì)的重要內(nèi)容和摩擦副材料選擇的重要理論依據(jù).

制動摩擦表面溫度模型早期由Bolk[4]提出，后來Jaeger對該理論進(jìn)行了發(fā)展，論證并提出了在半無限體表面矩形移動熱源的數(shù)學(xué)模型[5]，此后的研究工作大多集中在多熱源以及熱源之間的相互作用對摩擦表面溫度的影響.華林[6]等人基于盤式制動器制動過程中能量耗散的研究，建立了緊急制動過程中制動盤與摩擦片瞬態(tài)溫度場分析的有限元模型.王新波[7]等人，運(yùn)用循環(huán)迭代的計(jì)算方法加載熱流密度.該方法可以較好的模擬制動盤周向的應(yīng)力應(yīng)變分布，但其并未考慮制動盤各節(jié)點(diǎn)徑向旋轉(zhuǎn)線速度差異對實(shí)際制動過程的影響.

本文根據(jù)某型號大兆瓦風(fēng)電制動器的實(shí)際制動工況，并考慮了溫度場、結(jié)構(gòu)場在制動閘片上的耦合問題.采用ANSYS軟件建立制動過程的三維瞬態(tài)熱-結(jié)構(gòu)耦合模型.重點(diǎn)研究制動區(qū)域摩擦熱流的加載方式，提出速度梯度循環(huán)加載的方法，并將在笛卡爾坐標(biāo)系下的傳統(tǒng)耦合模型用柱坐標(biāo)表達(dá)，簡化了熱載荷的加載.解釋了TEI現(xiàn)象對制動閘片應(yīng)力及應(yīng)變分布的影響并預(yù)估其磨損狀況.最后通過與傳統(tǒng)均勻加載方法對比論證速度梯度循環(huán)法分析的可靠性.

1 制動器熱-結(jié)構(gòu)耦合模型的建立

1.1 熱傳導(dǎo)模型的建立

盤式制動器在進(jìn)行制動工作時，將生成大量摩擦熱，并使制動盤、閘片的溫度呈現(xiàn)梯度分布.盤式制動器的主要傳熱方式為熱傳導(dǎo)和熱對流，由于制動過程中溫度較高，也會產(chǎn)生少量熱輻射.因此在建立熱傳導(dǎo)模型時，本文做以下假設(shè)：

(1) 制動的過程中摩擦系數(shù)μ不變；

(2) 制動時，制動盤、閘片的接觸面上各個點(diǎn)的瞬時溫度對應(yīng)相等，熱流密度根據(jù)盤與閘片的材料屬性自然分配；

(3) 制動盤、閘片的接觸界面為理想平面；

(4) 制動過程中的摩擦功全部轉(zhuǎn)化為摩擦熱.因此在計(jì)算時，可將制動盤、制動閘片的熱流輸入作為邊界熱流輸入，則摩擦表面輸入熱流密度滿足.

(1)

為簡化后續(xù)熱流密度加載算法，選擇柱坐標(biāo)系，可將式(1)改寫為

(2)

其中，q(r,θ,t)為t時刻制動盤柱坐標(biāo)下(r,θ)處吸收的熱流密度，ζ為摩擦熱分配比例，μ為摩擦系數(shù)，p(r,θ,t)為摩擦副的面壓力，ω(r,θ,t)為制動盤角速度.

(5)制動盤與制動閘片的材料為各向同性材料，在制動過程中材料的熱物性參數(shù)不隨溫度變化而改變.

在制動過程中，制動盤與制動閘片之間存在自然熱傳導(dǎo)和摩擦熱流密度輸入[8].根據(jù)假設(shè)(3)制動閘片的溫度與制動盤的溫度相等，即自然熱傳導(dǎo)和摩擦熱流密度輸入關(guān)系為：

(3)

其中，λp為制動閘片熱導(dǎo)率，Tp為制動閘片的溫度.

1.2 熱應(yīng)力的計(jì)算

假設(shè)材料在制動過程中只發(fā)生線彈性變形.物體由于熱膨脹只產(chǎn)生線應(yīng)變，剪切應(yīng)變?yōu)榱?

這種由于材料熱變形產(chǎn)生的應(yīng)變可以看作是物體的初應(yīng)變.在計(jì)算應(yīng)力時包含初應(yīng)變項(xiàng)，代入熱流密度可以得出熱應(yīng)力的關(guān)系式為

(4)

其中，σ為材料應(yīng)力矩陣，D為材料彈性矩陣，ε為材料應(yīng)變矩陣，Δε為溫度變化引起的溫度應(yīng)變.

對于三維模型，綜合制動閘片溫度表達(dá)式[9]得出：

(5)

其中，α為材料應(yīng)力矩陣，Tp為制動閘片的溫度，T0為初始溫度.

1.3 邊界條件

根據(jù)制動器實(shí)際工況，可以認(rèn)為制動閘片不動，制動盤做逆時針轉(zhuǎn)動.由于制動壓力作用在制動閘片摩擦面背面，假設(shè)作用在此作用面上的壓力是均布的，但考慮制動盤、制動閘片的彈性變形.制動閘片是固接在制動連桿上并將閘片兩側(cè)凹槽安裝在導(dǎo)軌上可沿導(dǎo)軌上下滑動，故對制動閘片施加X，Z軸兩個方向的固定約束，如圖1所示.

圖1 制動摩擦副相對位置關(guān)系及位移約束示意圖

2 熱流密度速度梯度循環(huán)加載方法

熱流密度的加載是溫度場模擬的關(guān)鍵，溫度場的準(zhǔn)確模擬也是熱應(yīng)力的計(jì)算以及后續(xù)得出有效熱結(jié)構(gòu)耦合分析結(jié)果的基礎(chǔ).

為了給制動區(qū)域的每一摩擦作用單元準(zhǔn)確分配實(shí)際工況下的熱流密度，在加載式的推導(dǎo)過程中采用歸一化思想用以確定比例系數(shù).

2.1 熱流密度加載式的理論推導(dǎo)

由式(2)可知，熱流密度是關(guān)于摩擦熱分配比例、摩擦系數(shù)、摩擦副的面壓力、角速度以及徑向半徑的函數(shù)，由假設(shè)(1)、(3)可知：摩擦系數(shù)μ不變，制動閘片受均布壓力，故摩擦副的面壓力p不變，因此熱流密度與制動盤的轉(zhuǎn)速和半徑r有關(guān).且熱流密度與柱坐標(biāo)半徑(制動盤旋轉(zhuǎn)軸距制動區(qū)域的每一摩擦作用單元的距離)成正比，即

(6)

設(shè)比例系數(shù)為k1，因此公式可改寫為

(7)

在制動過程中整個制動區(qū)域每一時刻的總熱流密度

(8)

其中，n為此有限元模型中制動區(qū)域在該時刻的有效摩擦面積上的單元數(shù)；qi為制動區(qū)域任意摩擦作用單元的熱流密度.

將式(8)適當(dāng)變形從而求得比例系數(shù)

(9)

將式(9)代入式(8)最終得到每一摩擦作用單元的熱流密度

(10)

2.2 摩擦功加載式理論推導(dǎo)

基于假設(shè)(4)摩擦熱全部由摩擦力做功產(chǎn)生，所以本小節(jié)首先對摩擦力的加載進(jìn)行理論推導(dǎo).

制動區(qū)域摩擦作用點(diǎn)的摩擦功W與摩擦線速度v成正比.因?yàn)?/p>

v=ωr

則有：

W∝r

設(shè)比例系數(shù)為k2，因此公式可改寫為

W=k2r

(11)

所以只需要計(jì)算出k2值就能較為準(zhǔn)確的求解制動區(qū)域每一摩擦作用點(diǎn)的摩擦功Wi.

由于制動過程中整個制動區(qū)域每一時刻的總摩擦功

=k2r1+k2r2+…+k2ri+

(12)

其中，n為此有限元模型中制動區(qū)域在該時刻的有效摩擦面積上的節(jié)點(diǎn)數(shù).

將式(12)適當(dāng)變形從而求得比例系數(shù)

(13)

將式(13)代入(11)最終得到每一摩擦作用點(diǎn)的摩擦力

(14)

3 基于速度梯度循環(huán)法的熱-結(jié)構(gòu)耦合的算例分析

以某型號兆瓦級風(fēng)電制動器的結(jié)構(gòu)為例，制動閘片的材料為銅基粉末冶金復(fù)合材料.在有限元計(jì)算分析中，使用ANSYS中APDL模塊編寫相應(yīng)命令流，實(shí)現(xiàn)材料屬性邊界條件的施加，其中制動盤制動過程的初始轉(zhuǎn)速為2 000 r/min，施加制動壓力載荷為17 000 N.基于速度梯度循環(huán)法編寫熱流密度加載函數(shù)，在制動閘片頂面施加制動載荷，然后通過函數(shù)加載的方式施加制動閘片底面摩擦區(qū)域的摩擦熱流，設(shè)置載荷步進(jìn)行熱分析，之后將前面生成的溫度場中的熱單元轉(zhuǎn)化成結(jié)構(gòu)單元，加載摩擦載荷，最后將熱分析生成的*.rth文件作為體載荷加載到制動閘片結(jié)構(gòu)分析文件上，時間步長及其加載時間通過全局命令流控制.分析計(jì)算所需材料參數(shù)如表1所示.

圖2為制動閘片的約束條件及載荷施加示意圖，根據(jù)實(shí)際工況要求假設(shè)制動盤逆時針轉(zhuǎn)動，以仿真坐標(biāo)系作為參考坐標(biāo)系，由上一部分對邊界條件的分析可知制動閘片左側(cè)凹槽表面受X，Z方向的約束、底而受Y方向約束，并在如圖2所示表面施加壓力及模擬施加熱流密度及摩擦力.圖3表示了熱-結(jié)構(gòu)耦合理論流程.

圖2 制動閘片的約束條件及載荷示意圖

圖3 熱-結(jié)構(gòu)耦合過程示意圖

根據(jù)以上理論、邊界條件、有限元模型、計(jì)算公式，采用ANSYS求解器對制動器制動過程溫度場進(jìn)行計(jì)算，得到制動閘片的溫度場、結(jié)構(gòu)場分布.

3.1 制動器閘片的溫度場分布特性

圖4為采用速度梯度循環(huán)法加載的制動閘片在制動過程結(jié)束即19.16 s時的瞬態(tài)溫度場分布云圖，從圖4(a)可以看出，制動過程中制動閘片溫度分布不均勻，以圖中右側(cè)為Z軸正向則閘片等溫線向右側(cè)傾斜即呈“左高右低”分布且閘片溫度值也服從“左高右低”的分布規(guī)律，圖4(b)閘片底面高溫區(qū)集中在沿制動盤半徑方向的外側(cè)(圖中為Z軸負(fù)向)，溫度最大值分布在在閘片底面遠(yuǎn)離制動盤旋轉(zhuǎn)中心的端點(diǎn)處.

(a) 閘片左側(cè)面

(b) 閘片下底面

制動過程產(chǎn)生的摩擦熱流作用于閘片底面并向閘片厚度方向傳熱，但摩擦熱流在底面分布是不均勻的，由于閘片沿制動盤半徑方向的外側(cè)線速度大于內(nèi)側(cè)，而線速度大往往導(dǎo)致摩擦熱作用明顯并且產(chǎn)生熱流大，所以外側(cè)溫度高于內(nèi)側(cè)溫度.

3.2 制動器閘片的應(yīng)變分布特性對比分析

由于在建模過程中對閘片的位移約束施加在如圖2所示方向，所以圖5所示制動閘片下底面合位移分布云圖中閘片沿X軸方向的變形左側(cè)大于右側(cè)，圖5(a)為采用平均熱流密度分布法加載的制動閘片在19.16 s時的瞬態(tài)溫度場分布及閘片在制動器位移分布云圖，與圖5(b)采用速度梯度循環(huán)法的云圖對比，可以發(fā)現(xiàn)兩圖制動閘片應(yīng)變均沿制動盤相對運(yùn)動方向增大并沿制動盤徑向遞增，最大變形區(qū)域集中在沿X軸正方向的右側(cè)，但在沿制動盤半徑方向圖5(b)比圖5(a)變化趨勢明顯，并且圖5(b)應(yīng)變變化區(qū)域呈由X軸偏向Z軸方向的波浪線分布.

(a) 閘片左側(cè)面

(b) 閘片下底面

閘片在X軸方向受到位移約束，其摩擦應(yīng)變由遠(yuǎn)離約束端到靠近約束端遞減；又因?yàn)殚l片沿制動盤徑向的摩擦線速度遞增使單位面積產(chǎn)生的熱流密度遞增，閘片底面高溫區(qū)集中在沿制動盤半徑方向的外側(cè)，由于耦合作用的存在，使閘片摩擦應(yīng)變沿X軸正向偏向Z軸負(fù)向遞增.這也驗(yàn)證了使用速度梯度循環(huán)法與傳統(tǒng)均勻加載方法相比的準(zhǔn)確性.

3.3 制動器閘片的應(yīng)力分布特性對比分析

圖6為制動閘片在19.16 s時的等效應(yīng)力分布云圖，其分別采用均勻加載法如圖6(a)和速度梯度循環(huán)法如圖6(b)加載.由于閘片主要受作用于閘片頂面的壓力載荷和底面與制動盤摩擦產(chǎn)生的摩擦力，將速度梯度循環(huán)法與均勻加載方法相比較，可以看出采用速度梯度循環(huán)法加載處理后，閘片底面的等效應(yīng)力分布沿制動盤方向由內(nèi)向外遞增.由于制動器實(shí)際制動過程中制動盤逆時針旋轉(zhuǎn)而閘片左側(cè)面受到導(dǎo)軌的約束，制動盤、閘片摩擦作用產(chǎn)生的應(yīng)力集中在閘片的左側(cè)，又因?yàn)檠刂苿颖P半徑方向速度梯度的存在，從而導(dǎo)致外側(cè)摩擦作用更為明顯，而速度梯度循環(huán)法可有效模擬該工況，使仿真分析結(jié)果與實(shí)際相一致.

(a) 采用均勻加載方法

(b) 采用速度梯度循環(huán)法

結(jié)合圖4，沿制動盤半徑方向制動閘片摩擦面外側(cè)材料累積的熱量高于材料內(nèi)側(cè)，因此閘片外側(cè)溫度明顯高于內(nèi)側(cè)溫度，從而在制動盤徑向產(chǎn)生溫度變化梯度，這一溫度梯度使得制動閘片摩擦區(qū)域變形要向徑向擴(kuò)展，但受到非溫升的約束而不能自由進(jìn)行，導(dǎo)致閘片在摩擦高溫區(qū)域沿制動盤徑向產(chǎn)生了很大的壓應(yīng)力.

由于閘片沿制動盤半徑方向(圖中為Z軸負(fù)向)等效應(yīng)力分布靠向外側(cè)，又有外側(cè)線速度大于內(nèi)側(cè)，而底面高溫區(qū)又分布在外側(cè)，由摩擦學(xué)相關(guān)理論[11]可推測出制動閘片底面材料的磨損在外側(cè)要比內(nèi)側(cè)嚴(yán)重.

4 結(jié)論

本文基于大兆瓦風(fēng)電制動器制動摩擦副的結(jié)構(gòu)特性，考慮作用于制動閘片的相對線速度徑向差異及其影響下的摩擦熱流與應(yīng)力的耦合過程，提出速度梯度循環(huán)法，并推導(dǎo)出熱流密度的加載式，建立了大兆瓦風(fēng)電制動器三維熱-結(jié)構(gòu)耦合瞬態(tài)有限元模型對制動閘片熱-結(jié)構(gòu)耦合過程進(jìn)行分析.得出以下結(jié)論：

(1)在風(fēng)電制動器制動過程中，制動閘片的溫度場與應(yīng)力場之間存在耦合關(guān)系，使用速度梯度循環(huán)法可有效模擬實(shí)際工況.并得出制動閘片摩擦區(qū)域沿制動盤半徑方向呈現(xiàn)出弧度明顯的等溫帶且溫度梯度沿半徑尺寸增大方向遞增，較準(zhǔn)確的反應(yīng)制動過程中閘片溫度場與應(yīng)力應(yīng)變分布特性；

(2)制動閘片沿制動盤半徑尺寸方向外側(cè)溫度明顯高于內(nèi)側(cè)溫度，從而閘片在制動盤半徑方向產(chǎn)生溫度梯度，導(dǎo)致閘片在摩擦高溫區(qū)域沿制動盤徑向產(chǎn)生了很大的壓應(yīng)力.且通過以上分析得出制動閘片摩擦表面的磨損沿制動盤半徑增大方向外側(cè)比內(nèi)側(cè)嚴(yán)重.

[1]黃健萌, 高誠輝. 盤式制動器熱-結(jié)構(gòu)耦合的數(shù)值建模與分析[J].機(jī)械工程學(xué)報(bào), 2008(2):145-151.

[2]KWANGJIN L, BARBER JR. Frictionally excited thermoelastic instability in automotive disk brakes[J]. ASME Journal of Tribology, 1993, 115: 607-614.

[3]ZAGRODZKI P, LAM KB, BAHKALI EA, et al. Nonlinear transient behavior of a sliding system with frictionally excited thermoelastic instability[J]. ASME Journal of Tribology, 2001, 123:699-708.

[4]BOLK H. The flash temperature concept[J]. Wear, 1963, 6:483-493.

[5]EVTUSHENKO OO, IVANYK EH, HORBACHOVA NV. Analytic methods for thermal calculation of brakes(review)[J]. Materials Science, 2000, 36 (6):857-862.

[6]華林, 向上升.汽車氣壓盤式制動器瞬時溫度場研究[J].潤滑與密封, 2007, 32(5): 8-11.

[7]王新波.風(fēng)電制動器的熱-結(jié)構(gòu)耦合分析[D]大連:大連理工大學(xué), 2010.

[8]馬保吉, 朱均.摩擦制動系統(tǒng)摩擦表面的摩擦學(xué)研究[J].機(jī)械科學(xué)與技術(shù), 1999(1):14-16.

[9]胡育勇.風(fēng)機(jī)主軸制動器摩擦副熱-力耦合有限元分析[D].南昌:南昌大學(xué), 2014.

[10]劉建秀.銅基粉末冶金摩擦材料研制及其高溫疲勞磨損和沖擊性能研究[D].河南:中國工程物理研究院, 2003.

[11]溫詩鑄.摩擦學(xué)原理[M].4版，北京：清華大學(xué)出版社, 2012.

Study on Coupling Analysis of Thermal-Structure for Megawatt Wind Turbine Brakes

SHA Zhihua, LIU Nan, LIU Yu, MA Fu Jian, YIN Jian, ZHANG Shengfang

(School of Mechanical Engineering, Dalian Jiaotong University, Dalian 116028)

In view of the braking process for megawatt wind turbine brakes and considering the effects that the width of brake friction area and linear velocity radial difference, the gradient of velocities circulation method is put forward, and the thermal-structure for the brake in braking process is analyzed. First, the model of three dimensional transient heat transfers for brake friction pair is build based on ANSYS. Then heat flowing density is loaded by the gradient of velocities circulation method. At last, geometrical parameters are determined, and the loading function of heat flowing density is plugged. The result of temperature field analysis in the disc brake friction slice of temperature radius direction presents a radian obvious temperature zone. Direction of increased temperature gradient along the radius size is reduced. Using gradient of velocities circulation method integrated the results of thermal analysis generation into the structure under coupling analysis of force and deformation analysis of brake friction slice, wear condition can be forecasted. Compared with traditional uniform loading method using the analysis of the velocity gradient method is closer to the actual results. The analysis method provides the reference for simulation of the large size disc brake friction braking condition process.

wind turbine brakes; thermal-structural coupling; gradient of velocities circulation method; finite element simulation

1673- 9590(2016)06- 0066- 06

2016-06-03

國家自然科學(xué)基金資助項(xiàng)目 (51475066); 遼寧省自然科學(xué)基金資助項(xiàng)目(2015020114); 遼寧省教育廳優(yōu)秀人才計(jì)劃資助項(xiàng)目(LR2015012)

沙智華(1973-)，女，教授，博士，主要從事CAD/CAM/CAE，機(jī)械結(jié)構(gòu)優(yōu)化設(shè)計(jì)方面

A

E- mail:zhsha@djtu.edu.cn.