抽水蓄能機組推力頭的強度與壽命分析*

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

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

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

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

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

1 模型的建立

1.1 推力軸承結構

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

圖1 推力軸承結構示意

1.2 三維模型的建立

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

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

為便于分析，將計算部分分為2部分：流體區域和固體區域[4]。流體部分中的油膜部分采用結構化網格，在保證計算精確性且不浪費資源的情況下，經網格數量無關性檢驗后確定油膜厚度方向分為6層，油膜網格數目保持在62萬；除流體區域的油膜部分，其余流體部分均采用非結構網格，網格數目為93萬；固體計算區域為推力頭，可通過ANSYS19.2中的Static Structural中的mesh確定網格單元尺寸為0.01 m，生成非結構網格,見圖3。

圖3 推力頭網格劃分

1.3 數學模型的建立

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

質量守恒方程：

(1)

動量守恒方程：

(2)

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

總能量方程：

vρff+SE

(3)

式中：htot為總晗；λ為導熱系數；SE為能量源項；ΔT為溫度變化量。

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

(4)

β″ω2+Dω

(5)

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

(6)

式中：下標f為流體；s為固體。

2 計算過程及理論

2.1 流體部分的計算方法

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

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

2.2 流固耦合計算流程

在穩態運行時，推力頭與油膜接觸面各點所受壓力會引起相應各點產生一定的位移量，該問題屬于流固耦合問題，可通過ANSYS19.2中的Static Structural進行相應的計算[13]。

圖4所示為基于流固耦合理論的計算流程。首先進行流體區域的數值計算，然后將計算得到的流固耦合面處各點的壓力傳遞到固體區域進行穩態結構計算，經過計算得出各網格結點的位移變化，最后根據耦合面位移變化進行流體區域相應面網格的變形及體網格重構，再進行下一輪的流體計算。重復前面的計算步驟，經過若干輪計算與傳遞后停止計算，此時各處變形值及壓力值滿足設定的收斂條件。

圖4 流固耦合計算流程

2.3 結構場約束設置

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

2.4 材料的疲勞計算

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

SαN=C

(7)

式中：α、C為材料常數；S為應力；N為應力循環次數。

式(7)兩邊取對數得：

lgN=lgC-αlgS

(8)

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

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

圖6 金屬的S-N曲線

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

圖7 交變循環應力

3 計算結果與分析

3.1 油膜壓力分析

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

圖8 油膜壓力隨轉速和油膜厚度變化

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

圖9 油膜最大壓力和平均壓力隨轉速和油膜厚度變化

3.2 推力頭應變分析

在機組啟動過程中，選取轉速分別為100、250、350、500 r/min，分析了推力頭下的總變形分布情況，結果如圖10所示。可見不同轉速下對應的推力頭應變分布情況基本一致。推力頭表面的最大形變分布主要集中在其外緣，從推力頭外緣向內側形變逐漸減小，可能是由于外緣遠離主軸，外緣受到的力矩大于靠近主軸處所受到的力矩，因此推力頭的外緣形變比較大。圖11示出了不同轉速下推力頭最大、最小和平均形變。整體上隨著機組轉速的升高推力頭的最大形變值有所增大，且在轉速100～250 r/min之間推力頭的應變增加較快，在轉速超過250 r/min后推力頭的形變增加速度逐漸放緩。轉速100 r/min下的最大形變與轉速500 r/min下的最大形變相差近65%。雖然不同轉速下的最大形變相差較大，但是最大形變相較于推力頭尺寸都屬于小范圍變形[4]。如推力頭外緣的半徑為2 125 mm，而500 r/min所對應的最大變形量僅為1.266 8 mm。因此，不同轉速下引起的推力頭形變對機組安全穩定運行影響較小。

圖10 不同轉速下推力頭應變分布

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

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

圖12 不同油膜厚度下推力頭應變分布

圖13 油膜厚度66 μm下鏡板面應變

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

3.3 推力頭壽命分析

根據S-N曲線對不同運行參數推力頭進行疲勞壽命計算，得到的循環次數的最大和最小值如圖15所示。

圖15 不同轉速和油膜厚度下推力頭應力循環次數

圖15(a)所示為轉速為500 r/min，不同油膜厚度下的推力頭應力循環次數，隨著油膜厚度的增加，在油膜厚度為76 μm處出現了谷值，可能是因為在油膜厚度為76 μm時所對應的壓力出現峰值,加劇了應力集中的程度?？傮w而言，機組在運行過程中產生了作用于推力頭上的交變應力，并出現不同程度的應力集中，應考慮對推力頭應力集中處進行表面強化或者改善其結構。

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

從圖16所示的額定轉速推力頭應力循環次數云圖可看出，推力頭出現了局部壽命較低的情況(圖中圓圈部分)，所以從100 r/min開始出現最小應力循環次數降低的情況。從圖17所示推力頭應力分布云圖，可發現推力頭出現壽命較低的位置均出現了局部應力集中，這可能是因為該處推力頭外形尺寸發生突變從而引起局部范圍的應力集中，使該處壽命較低。

圖16 額定轉速推力頭應力循環次數云圖

圖17 額定轉速推力頭應力分布云圖

4 結論

運用流固耦合理論，通過Workbench對推力頭進行了應變和壽命分析，并對比了不同運行參數下的油膜壓力分布，得到如下結論：

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

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

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