彭緒意 徐利君 胥千鑫 章志平 劉 澤 王雪梅 張 智 張玉全 鄭 源

(1.江西洪屏抽水蓄能有限公司 江西宜春 330603；2.河海大學(xué)能源與電氣學(xué)院 江蘇南京 211100；3.河海大學(xué)水利水電學(xué)院 江蘇南京 210098)

抽水蓄能電站在電力系統(tǒng)中具有調(diào)峰填谷、調(diào)頻、調(diào)相、緊急事故備用和黑啟動(dòng)等多種功能， 成為現(xiàn)代電力系統(tǒng)中有效的、不可或缺的調(diào)節(jié)工具[1-2]。抽水蓄能機(jī)組在工作時(shí)啟停頻繁、工作轉(zhuǎn)速高，在不同工況時(shí)需要調(diào)整旋轉(zhuǎn)運(yùn)行方向，同時(shí)由于其結(jié)構(gòu)鋼強(qiáng)度低易受水力、機(jī)械、電磁等復(fù)雜因素的影響，機(jī)組振動(dòng)、軸承磨損問題突出[3]。因此研究抽水蓄能機(jī)組運(yùn)行過程中推力軸承的應(yīng)變和疲勞壽命，以提高機(jī)組運(yùn)行穩(wěn)定性顯得十分必要。

近些年，CFD數(shù)值模擬技術(shù)得到快速發(fā)展，已廣泛應(yīng)用于流體計(jì)算的各個(gè)領(lǐng)域，基于流固耦合理論研究推力軸承也受到普遍關(guān)注[4]。翟黎明[3]運(yùn)用流、固、熱耦合的方法研究了推力軸承熱彈流潤滑特性，分析了軸承的靜壓潤滑特性。屈波等人[4]運(yùn)用雙向流固耦合，通過CFD技術(shù)對(duì)抽水蓄能機(jī)組推力軸承進(jìn)行數(shù)值模擬，研究了瓦面傾角、油膜厚度以及瓦面彈性模量對(duì)其潤滑性能的影響。李永海等[5]依據(jù)傳熱學(xué)、摩擦學(xué)等相關(guān)理論，研究了推力瓦表面溫度和對(duì)流換熱系數(shù)的施加方法以及推力瓦的熱變形分布。王青華等[6]利用計(jì)算流體力學(xué)方法，研究了油膜厚度和瓦面傾角對(duì)油膜承載力、等效應(yīng)力的影響。梁興鑫等[7]設(shè)計(jì)階梯橡膠支撐的水潤滑可傾瓦推力軸承，并基于流固耦合理論計(jì)算了橡膠墊不同尺寸和瓦面材料對(duì)推力盤的影響。唐乾皓等[8]基于數(shù)值模擬，通過優(yōu)化推力瓦幾何參數(shù)解決了推力瓦溫度過高的問題，同時(shí)也證明了推力瓦幾何參數(shù)是影響推力軸承油系統(tǒng)的重要參數(shù)。以上研究均是對(duì)推力軸承中推力瓦的分析，對(duì)于推力軸承中推力頭的研究較少。然而推力頭在機(jī)組啟停頻繁的情況下，承受交變應(yīng)力，容易產(chǎn)生疲勞磨損，這也成為機(jī)組安全運(yùn)行的一大隱患。

本文作者在分析江西洪屏抽水蓄能電站推力軸承結(jié)構(gòu)及運(yùn)行數(shù)據(jù)的基礎(chǔ)上[9]，建立了推力頭流固耦合模型，運(yùn)用CFD技術(shù)，將油膜的壓力分布加載到鏡板面上，分析了推力頭在不同運(yùn)行參數(shù)下應(yīng)變以及壽命情況，為工程應(yīng)用提供指導(dǎo)意義。

1 模型的建立

1.1 推力軸承結(jié)構(gòu)

推力軸承作為水力機(jī)械的重要部件之一，承擔(dān)機(jī)組的軸向載荷，保證機(jī)組的安全穩(wěn)定運(yùn)轉(zhuǎn)。一般情況下，推力軸承由推力頭、鏡板、推力瓦、油箱、支撐系統(tǒng)等組成，如圖1所示為推力軸承結(jié)構(gòu)示意圖。推力軸承的推力頭下方和鏡板連接，隨機(jī)組一同旋轉(zhuǎn)，推力瓦則為靜止部件。文中模擬的抽水蓄能電站的推力軸承共有12塊扇形瓦塊，繞主軸均勻分布。瓦塊被油槽內(nèi)的潤滑油完全浸沒，可以保證瓦面和鏡板面間的潤滑，避免干摩擦。當(dāng)機(jī)組啟動(dòng)前，壓力油在軸瓦表面形成靜壓油膜。在機(jī)組運(yùn)行過程中，經(jīng)冷卻后的潤滑油流入瓦間油槽，通過推力頭將冷卻的潤滑油帶入推力瓦與鏡板之間，形成動(dòng)壓油膜。

圖1 推力軸承結(jié)構(gòu)示意

1.2 三維模型的建立

根據(jù)某抽水蓄能電站推力頭平面尺寸圖，通過UG建立推力頭的三維實(shí)體模型。如圖2所示為推力頭的1/12實(shí)體模型及尺寸。

圖2 推力頭三維模型及尺寸(mm)

為便于分析，將計(jì)算部分分為2部分：流體區(qū)域和固體區(qū)域[4]。流體部分中的油膜部分采用結(jié)構(gòu)化網(wǎng)格，在保證計(jì)算精確性且不浪費(fèi)資源的情況下，經(jīng)網(wǎng)格數(shù)量無關(guān)性檢驗(yàn)后確定油膜厚度方向分為6層，油膜網(wǎng)格數(shù)目保持在62萬；除流體區(qū)域的油膜部分，其余流體部分均采用非結(jié)構(gòu)網(wǎng)格，網(wǎng)格數(shù)目為93萬；固體計(jì)算區(qū)域?yàn)橥屏︻^，可通過ANSYS19.2中的Static Structural中的mesh確定網(wǎng)格單元尺寸為0.01 m，生成非結(jié)構(gòu)網(wǎng)格,見圖3。

圖3 推力頭網(wǎng)格劃分

1.3 數(shù)學(xué)模型的建立

流體流動(dòng)必須遵循基本的物理守恒定律，包括質(zhì)量守恒定律、動(dòng)量守恒定律和能量守恒定律[10]。基于N-S方程通過式(1)—(3)控制方程描述流場(chǎng)。

質(zhì)量守恒方程：

(1)

動(dòng)量守恒方程：

(2)

式中：ft為體積力矢量；tf為剪切力張量。

總能量方程：

vρff+SE

(3)

式中：htot為總晗；λ為導(dǎo)熱系數(shù)；SE為能量源項(xiàng)；ΔT為溫度變化量。

在進(jìn)行流體計(jì)算時(shí)考慮到，推力軸承中的流體既有層流流動(dòng)又有湍流流動(dòng)，采用κ-ωSST模型可以很好地處理低雷諾數(shù)到高雷諾數(shù)的變化，對(duì)應(yīng)的湍動(dòng)能κ以及湍流頻率ω方程如公式(4)和(5)所示。

(4)

β″ω2+Dω

(5)

在流固耦合計(jì)算中，需要通過交接面進(jìn)行流體和固體的數(shù)據(jù)傳導(dǎo)，數(shù)據(jù)傳導(dǎo)需要滿足壓力、位移、熱流以及溫度守恒，如式(6)[10]所示。

(6)

式中：下標(biāo)f為流體；s為固體。

2 計(jì)算過程及理論

2.1 流體部分的計(jì)算方法

機(jī)組運(yùn)行時(shí)，推力頭主要會(huì)受到離心力和壓力，在對(duì)推力頭進(jìn)行靜力學(xué)分析時(shí)，需要獲得油膜的壓力分布。由于k-omega SST對(duì)邊界有較好的處理方式，故采用計(jì)算結(jié)果較好的k-omega SST湍流模型[11]。

模型邊界條件設(shè)置為：進(jìn)口為推力瓦進(jìn)油側(cè)所對(duì)應(yīng)的模型位置，設(shè)為速度進(jìn)口并且設(shè)置進(jìn)油溫度為35 ℃，進(jìn)口速度為0.5 m/s；出口為模型的另一側(cè)，設(shè)置為自由出口；鏡板面設(shè)置為旋轉(zhuǎn)壁面，旋轉(zhuǎn)速度為500 r/min并設(shè)置溫度為70 ℃。壁面采用無滑移邊界條件[12]。采用SIMPLEC對(duì)速度壓力進(jìn)行耦合求解。其中壓力項(xiàng)采用二階迎風(fēng)格式，對(duì)流項(xiàng)、湍動(dòng)能及耗散率均采用二階迎風(fēng)格式。

2.2 流固耦合計(jì)算流程

在穩(wěn)態(tài)運(yùn)行時(shí)，推力頭與油膜接觸面各點(diǎn)所受壓力會(huì)引起相應(yīng)各點(diǎn)產(chǎn)生一定的位移量，該問題屬于流固耦合問題，可通過ANSYS19.2中的Static Structural進(jìn)行相應(yīng)的計(jì)算[13]。

圖4所示為基于流固耦合理論的計(jì)算流程。首先進(jìn)行流體區(qū)域的數(shù)值計(jì)算，然后將計(jì)算得到的流固耦合面處各點(diǎn)的壓力傳遞到固體區(qū)域進(jìn)行穩(wěn)態(tài)結(jié)構(gòu)計(jì)算，經(jīng)過計(jì)算得出各網(wǎng)格結(jié)點(diǎn)的位移變化，最后根據(jù)耦合面位移變化進(jìn)行流體區(qū)域相應(yīng)面網(wǎng)格的變形及體網(wǎng)格重構(gòu)，再進(jìn)行下一輪的流體計(jì)算。重復(fù)前面的計(jì)算步驟，經(jīng)過若干輪計(jì)算與傳遞后停止計(jì)算，此時(shí)各處變形值及壓力值滿足設(shè)定的收斂條件。

圖4 流固耦合計(jì)算流程

2.3 結(jié)構(gòu)場(chǎng)約束設(shè)置

在對(duì)結(jié)構(gòu)體進(jìn)行分析計(jì)算時(shí)，對(duì)結(jié)構(gòu)體添加了如圖5所示的約束條件。其中：A為角加速度(Rotational velocity)，推力頭在旋轉(zhuǎn)過程中受到離心力作用，設(shè)定旋轉(zhuǎn)速度為500 r/min。B為施加集中應(yīng)力(Force)，此處設(shè)定為推力頭的一個(gè)面上。C和D均為圓柱約束(Cylindrical support)，設(shè)定推力頭的切向轉(zhuǎn)動(dòng)。

2.4 材料的疲勞計(jì)算

根據(jù)相關(guān)手冊(cè)可以查到轉(zhuǎn)輪材料的疲勞極限值，對(duì)于簡化的S-N曲線可以采用冪函數(shù)繪制，如公式(7)所示[14]。

SαN=C

(7)

式中：α、C為材料常數(shù)；S為應(yīng)力；N為應(yīng)力循環(huán)次數(shù)。

式(7)兩邊取對(duì)數(shù)得：

lgN=lgC-αlgS

(8)

令a=lgC，b=-α，得lgN=a+blgS。

在Workbench中的Engineering Data中，獲得材料的S-N曲線。如圖6所示。

圖6 金屬的S-N曲線

在Workbench中采用Fatigue tools對(duì)推力頭進(jìn)行壽命計(jì)算，載荷類型設(shè)置為脈動(dòng)循環(huán)載荷，載荷比例由0到1按正弦比例隨時(shí)間變化如圖7所示。在計(jì)算疲勞壽命時(shí)用Goodman修正公式更為保守，在工程應(yīng)用中一般采用保守的計(jì)算方式[12]。因此選擇Goodman為推力頭疲勞壽命估算模型[13]。

圖7 交變循環(huán)應(yīng)力

3 計(jì)算結(jié)果與分析

3.1 油膜壓力分析

由于油膜極薄的特性，通過物理實(shí)驗(yàn)采集油膜參數(shù)是極其困難的，多數(shù)通過數(shù)值模擬進(jìn)行研究。為對(duì)推力頭進(jìn)行靜力學(xué)分析，計(jì)算了油膜在不同運(yùn)行參數(shù)下的壓力分布情況，以及油膜最大壓力和平均壓力隨轉(zhuǎn)速和油膜厚度的變化情況，分別如圖8、9所示。圖8(a)所示為機(jī)組在不同轉(zhuǎn)速下，沿旋轉(zhuǎn)半徑方向且過某一鏡板面中心的油膜壓力變化情況?？梢姴煌D(zhuǎn)速下的油膜壓力變化趨勢(shì)比較相近，都出現(xiàn)了沿半徑方向先增大后減小的趨勢(shì)，且均在半徑0.75 m左右出現(xiàn)壓力最高點(diǎn)。結(jié)合圖9(a)可知，隨著轉(zhuǎn)速的不斷升高，油膜表面的最大壓力和平均壓力均不斷增加。由于軸向推力與軸瓦和鏡板之間的承載動(dòng)壓、相對(duì)運(yùn)動(dòng)以及潤滑油的黏滯作用這些因素有關(guān)，隨著轉(zhuǎn)速的增加即相對(duì)運(yùn)動(dòng)速度增加得比較迅速，因此軸向推力不斷增加，相對(duì)應(yīng)的壓力也隨轉(zhuǎn)速不斷增加。

圖8 油膜壓力隨轉(zhuǎn)速和油膜厚度變化

圖8(b)所示為機(jī)組處于額定轉(zhuǎn)速500 r/min時(shí)，不同的油膜厚度下沿旋轉(zhuǎn)半徑方向且過某一鏡板面中心的油膜壓力情況。不同油膜厚度下沿半徑方向壓力分布呈現(xiàn)出和不同轉(zhuǎn)速下相同的變化趨勢(shì)，且在半徑0.75 m左右處出現(xiàn)壓力最高點(diǎn)。結(jié)合圖9(b)可知，當(dāng)油膜厚度從55 μm增加至76 μm時(shí)，油膜的壓力隨其厚度的增加而增加，而當(dāng)油膜厚度繼續(xù)增加時(shí)，油膜的壓力急劇下降。當(dāng)油膜厚度增加到76 μm時(shí)，最大壓力和平均壓力之間的差值最大，說明在機(jī)組運(yùn)行過程中，當(dāng)油膜在達(dá)到76 μm附近時(shí)，油膜表面將會(huì)產(chǎn)生劇烈的壓力梯度變化。這一方面說明由于變工況情況導(dǎo)致的油膜厚度變化，會(huì)嚴(yán)重影響油膜表面參數(shù)變化；另一方面說明，油膜厚度不會(huì)長久保持在76 μm，運(yùn)行期間油膜厚度達(dá)到76 μm附近時(shí)會(huì)短暫影響這一時(shí)期的機(jī)組安全穩(wěn)定運(yùn)行，后果嚴(yán)重的將導(dǎo)致燒瓦現(xiàn)象。油膜壓力變化趨勢(shì)與章志平等[14]的數(shù)值模擬結(jié)果相符。

圖9 油膜最大壓力和平均壓力隨轉(zhuǎn)速和油膜厚度變化

3.2 推力頭應(yīng)變分析

在機(jī)組啟動(dòng)過程中，選取轉(zhuǎn)速分別為100、250、350、500 r/min，分析了推力頭下的總變形分布情況，結(jié)果如圖10所示?？梢姴煌D(zhuǎn)速下對(duì)應(yīng)的推力頭應(yīng)變分布情況基本一致。推力頭表面的最大形變分布主要集中在其外緣，從推力頭外緣向內(nèi)側(cè)形變逐漸減小，可能是由于外緣遠(yuǎn)離主軸，外緣受到的力矩大于靠近主軸處所受到的力矩，因此推力頭的外緣形變比較大。圖11示出了不同轉(zhuǎn)速下推力頭最大、最小和平均形變。整體上隨著機(jī)組轉(zhuǎn)速的升高推力頭的最大形變值有所增大，且在轉(zhuǎn)速100～250 r/min之間推力頭的應(yīng)變?cè)黾虞^快，在轉(zhuǎn)速超過250 r/min后推力頭的形變?cè)黾铀俣戎饾u放緩。轉(zhuǎn)速100 r/min下的最大形變與轉(zhuǎn)速500 r/min下的最大形變相差近65%。雖然不同轉(zhuǎn)速下的最大形變相差較大，但是最大形變相較于推力頭尺寸都屬于小范圍變形[4]。如推力頭外緣的半徑為2 125 mm，而500 r/min所對(duì)應(yīng)的最大變形量僅為1.266 8 mm。因此，不同轉(zhuǎn)速下引起的推力頭形變對(duì)機(jī)組安全穩(wěn)定運(yùn)行影響較小。

圖10 不同轉(zhuǎn)速下推力頭應(yīng)變分布

圖11 不同轉(zhuǎn)速下推力頭最大、最小和平均形變分布

在轉(zhuǎn)速為500 r/min時(shí)取不同油膜厚度，通過有限元分析得到不同油膜厚度下推力頭應(yīng)變的分布云圖，如圖12所示。4種油膜厚度下的推力頭應(yīng)變均沿半徑方向由內(nèi)向外逐漸增大，在遠(yuǎn)離主軸位置處產(chǎn)生最大應(yīng)變，即推力頭表面應(yīng)變從靠近主軸處到遠(yuǎn)離主軸處逐漸變大。結(jié)合鏡板面應(yīng)變分析(見圖13)，推力頭呈現(xiàn)內(nèi)徑下沉外徑上翹的形狀。由于研究的4種不同油膜厚度下推力頭的最大形變量均大于油膜的厚度，因此推力頭的形變?cè)谝欢ǔ潭壬细淖兞擞湍さ暮穸确植?。油膜厚度?6～76 μm時(shí)，油膜的厚度分布受推力頭形變影響而產(chǎn)生的不均勻情況不斷加劇，當(dāng)厚度增加到86 μm時(shí)，油膜厚度分布的不均勻程度又有所改善。圖14示出了最大、最小和平均形變與油膜厚度的關(guān)系。當(dāng)油膜厚度為56～76 μm時(shí)，推力頭表面的最大形變值不斷增加，當(dāng)油膜厚度繼續(xù)增大時(shí)，推力頭表面的形變值急劇減小。圖14中形變隨油膜厚度的變化趨勢(shì)和圖9(b)中油膜表面壓力隨油膜厚度變化的趨勢(shì)相同，說明油膜表面壓力對(duì)推力頭的形變有一定程度的影響?？傮w而言，油膜厚度達(dá)到一定值時(shí)會(huì)對(duì)推力頭造成較大的形變，從而影響油膜的分布，這對(duì)推力軸承的潤滑和散熱都是不利的，在機(jī)組運(yùn)行時(shí)應(yīng)盡可能減小油膜厚度和推力頭最大形變之差，保證油膜厚度分布盡可能均勻。

圖12 不同油膜厚度下推力頭應(yīng)變分布

圖13 油膜厚度66 μm下鏡板面應(yīng)變

圖14 不同油膜厚度下推力頭最大、最小和平均形變分布

3.3 推力頭壽命分析

根據(jù)S-N曲線對(duì)不同運(yùn)行參數(shù)推力頭進(jìn)行疲勞壽命計(jì)算，得到的循環(huán)次數(shù)的最大和最小值如圖15所示。

圖15 不同轉(zhuǎn)速和油膜厚度下推力頭應(yīng)力循環(huán)次數(shù)

圖15(a)所示為轉(zhuǎn)速為500 r/min，不同油膜厚度下的推力頭應(yīng)力循環(huán)次數(shù)，隨著油膜厚度的增加，在油膜厚度為76 μm處出現(xiàn)了谷值，可能是因?yàn)樵谟湍ず穸葹?6 μm時(shí)所對(duì)應(yīng)的壓力出現(xiàn)峰值,加劇了應(yīng)力集中的程度?？傮w而言，機(jī)組在運(yùn)行過程中產(chǎn)生了作用于推力頭上的交變應(yīng)力，并出現(xiàn)不同程度的應(yīng)力集中，應(yīng)考慮對(duì)推力頭應(yīng)力集中處進(jìn)行表面強(qiáng)化或者改善其結(jié)構(gòu)。

圖15(b)所示為油膜厚度為65 μm時(shí)推力頭在不同轉(zhuǎn)速下的應(yīng)力循環(huán)次數(shù)。不同轉(zhuǎn)速下推力頭的最大應(yīng)力循環(huán)次數(shù)為1×106；在轉(zhuǎn)速從100 r/min增加到250 r/min的過程中出現(xiàn)了推力頭最小應(yīng)力循環(huán)次數(shù)急劇降低的情況，在轉(zhuǎn)速從250 r/min增加到500 r/min時(shí)最小應(yīng)力循環(huán)次數(shù)繼續(xù)小幅下降并趨于穩(wěn)定。

從圖16所示的額定轉(zhuǎn)速推力頭應(yīng)力循環(huán)次數(shù)云圖可看出，推力頭出現(xiàn)了局部壽命較低的情況(圖中圓圈部分)，所以從100 r/min開始出現(xiàn)最小應(yīng)力循環(huán)次數(shù)降低的情況。從圖17所示推力頭應(yīng)力分布云圖，可發(fā)現(xiàn)推力頭出現(xiàn)壽命較低的位置均出現(xiàn)了局部應(yīng)力集中，這可能是因?yàn)樵撎幫屏︻^外形尺寸發(fā)生突變從而引起局部范圍的應(yīng)力集中，使該處壽命較低。

圖16 額定轉(zhuǎn)速推力頭應(yīng)力循環(huán)次數(shù)云圖

圖17 額定轉(zhuǎn)速推力頭應(yīng)力分布云圖

4 結(jié)論

運(yùn)用流固耦合理論，通過Workbench對(duì)推力頭進(jìn)行了應(yīng)變和壽命分析，并對(duì)比了不同運(yùn)行參數(shù)下的油膜壓力分布，得到如下結(jié)論：

(1)隨著機(jī)組轉(zhuǎn)速的升高，油膜壓力也隨之升高；隨著油膜厚度的增加，油膜壓力先增加后減小。

(2)推力頭的最大形變主要集中在其外緣，隨著轉(zhuǎn)速的升高最大形變隨之增加?？傮w上形變量和推力頭的尺寸相比屬于小變形范圍，對(duì)機(jī)組的安全穩(wěn)定運(yùn)行影響較??；隨著油膜厚度的增加，不同油膜厚度下發(fā)生的最大形變均大于油膜厚度，在一定程度上，推力頭的形變影響了油膜的均勻分布，在實(shí)際運(yùn)行中，油膜的厚度和推力頭最大形變應(yīng)盡量接近，以保證機(jī)組安全穩(wěn)定的運(yùn)行。

(3)基于S-N曲線，分別模擬了不同轉(zhuǎn)速和不同油膜厚度下的推力頭疲勞壽命。隨著轉(zhuǎn)速升高，應(yīng)力循環(huán)次數(shù)在應(yīng)力集中處不斷降低；隨著油膜厚度增加，在油膜厚度為76 μm處出現(xiàn)了谷值。因此應(yīng)考慮通過對(duì)推力頭應(yīng)力集中處進(jìn)行表面強(qiáng)化或者改善其結(jié)構(gòu)，延長推力頭壽命。