基于流固耦合的抽水蓄能機組推力瓦變形以及安全性能分析

洪云來，秦 程，楊 雄，章志平，張送校，吳中華，王雪梅，張 智，鄭 源

（1.江西洪屏抽水蓄能有限公司，江西省宜春市 336000；2.河海大學能源與電氣學院，江蘇省南京市 211100；3.河海大學水利水電學院，江蘇省南京市 210098）

0 引言

近些年來，我國抽水蓄能電站建設取得了長足的進展。抽水蓄能機組啟動迅速、調(diào)節(jié)靈活、負荷變化范圍廣，在削峰填谷、調(diào)頻調(diào)壓方面有著重要作用[1]。由于抽水蓄能機組具有轉(zhuǎn)速高、工況復雜、工況切換頻繁等工作特點，軸承振動、軸承磨損等一系列關乎機組安全運行的問題變得尤為嚴峻[2，3]。因此，對推力軸承的結(jié)構、強度進行計算和分析十分必要。

目前，國內(nèi)主要采用流固耦合方法計算推力軸承強度[4]，屈波等[3]運用雙向流固耦合技術對抽水蓄能機組推力軸承進行數(shù)值模擬，研究了機組穩(wěn)態(tài)運行時推力軸承油膜的運行機制，分析總結(jié)了潤滑參數(shù)與彈性模量對軸瓦壓力分布的影響。于向軍等[5]在基于間接耦合的方法基礎上，通過流固耦合數(shù)值計算了大中型空軸式靜壓軸承，根據(jù)數(shù)據(jù)總結(jié)出軸承發(fā)生形變的規(guī)律以及影響形變的因素。孟凡明[6]通過運用非線性優(yōu)化法與影響系數(shù)法，編制了軸承流固耦合有限元程序，使用三維熱彈流分析計算了某水潤滑軸承，比較實驗結(jié)果，驗證了所編寫程序的有效與可靠性。

本文以江西洪屏抽水蓄能機組為對象，建立推力軸承流固耦合計算模型，先對油膜進行數(shù)值模擬，得到油膜的壓力分布，然后將其加載到結(jié)構分析的推力瓦模型表面，計算推力瓦的等效應力和變形，最后進行推力瓦的剛度分析，研究了機組運行整個過程中油膜運行機制對推力軸承產(chǎn)生的影響。

1 數(shù)學模型的建立

1.1 推力軸承基本理論

推力軸承在水力機械中起到了至關重要的作用，是專門承受軸向力的機械部件[7]。整個軸的軸向力通過油膜傳遞到推力軸承的推力瓦上。在機組啟動前，需通過高壓油頂起裝置讓壓力油到達推力瓦與鏡板之間，并在其之間形成靜壓油膜，隨著機組轉(zhuǎn)速不斷升高到達額定轉(zhuǎn)速時，高壓油泵停止工作，此時油膜依靠軸瓦與鏡板之間的楔形間隙形成，稱為動壓油膜[8]。江西洪屏抽水蓄能機組每個推力瓦下方采用彈簧簇支撐，采用可撕調(diào)整墊保證各推力瓦在同一水平面上，每塊推力瓦上均布置一個高壓油進油口，機組開啟時啟動高壓油頂起油泵[9]。

1.2 推力瓦三維模型

根據(jù)圖1所示江西洪屏抽水蓄能機組的瓦面尺寸圖，通過UG建立推力瓦的實體模型，其中，扇形瓦面兩邊各有一對稱的斜坡，扇瓦內(nèi)徑890mm，扇瓦外徑2010mm，額定轉(zhuǎn)速500r/min，瓦塊數(shù)為12[10]。圖2為推力瓦的實體三維模型。

圖1 瓦面結(jié)構與尺寸Figure 1 Tile surface structure and dimensions

圖2 推力瓦三維模型Figure 2 Thrust tile 3D model

2 流固耦合分析

2.1 流固耦合基本理論

流體計算時，假設所計算流體為連續(xù)介質(zhì)，均質(zhì)各向同性并且為不可壓縮流體，流體中的每個質(zhì)點均處于熱平衡狀態(tài)。考慮到推力軸承中的流體同時存在層流流動和湍流流動，因此，采用κ-ω SST湍流模型。計算時流體滿足時均連續(xù)性方程：

以及時均動量方程[11]：

對應的湍動能κ以及湍流頻率ω方程分為式（3）和式（4）：

式中：下標f為流體，s為固體；τ為流體與固體應力；u為位移；q為熱流量；T為溫度。

2.2 流固耦合計算流程

由于油膜較薄，穩(wěn)態(tài)運行時瓦面各點所受壓力會引起相應各點產(chǎn)生一定的位移量，而各點位移又會反過來作用于瓦面的壓力分布?？梢越柚鶤NSYS19.2中的Static Structural進行相應的單向流固耦合計算。

流體區(qū)域利用FLUENT對控制方程進行離散，采用SIMPLEC方法隱式求解速度和壓力項，壓力項采用二階中心差分格式，對流項、湍動能及耗散率均采用二階迎風格式[3]。進口位于兩塊軸瓦之間，因軸承流量一定，設置為速度進口條件，進油溫度為35℃，進口速度為0.5m/s；出口位于模型外側(cè)，設置為自由出流條件。鏡板面設置為wall，轉(zhuǎn)速為500r/min，設置油溫為70℃。模型的兩個側(cè)面設置為周期性邊界條件：φ1=φ2（其中：φ=u，v，w，k，∈），壁面為無滑移邊界條件。固體區(qū)域利用穩(wěn)態(tài)結(jié)構分析模塊Static Structural，巴氏合金層下表面設置固定面約束[12]。

流體域采用動網(wǎng)格技術，并將油膜層下表面與巴氏合金層上表面設置為流固耦合面。計算時先求解潤滑油流場，得到耦合面的流場壓力，通過System Coupling將壓力傳遞到結(jié)構分析模塊，完成耦合面流場壓力到固體變形的轉(zhuǎn)化，反之亦然，從而保證耦合面數(shù)據(jù)實時傳遞。

流固耦合的數(shù)值模擬流程圖如圖3所示，首先在A中進行流體區(qū)域的數(shù)值計算，計算完畢后將流體耦合面處各點的壓力傳遞到B中，然后在B中進行穩(wěn)態(tài)結(jié)構計算，經(jīng)過計算得出瓦面各點的位移變化。然后根據(jù)耦合面位移變化進行流體區(qū)域相應面網(wǎng)格的變形及體網(wǎng)格重構，接下來再進行下一輪的流體計算，之后重復前面的過程，經(jīng)過若干輪計算與傳遞后停止計算，此時瓦面各處變形值及壓力值滿足設定的收斂條件。

圖3 ANSYS流固耦合流程圖Figure 3 Flow - solid coupling flow chart of ANSYS

2.3 結(jié)構場設置

計算區(qū)域分為流體區(qū)域與固體區(qū)域。由于流體計算區(qū)域比較復雜，故采用分塊網(wǎng)格技術，各塊之間采用交界面連接，其中，油膜部分均采用結(jié)構化網(wǎng)格。經(jīng)網(wǎng)格數(shù)量無關性檢驗后確定油膜厚度方向分為6層，油膜網(wǎng)格數(shù)目為20萬；非油膜區(qū)域采用非結(jié)構化網(wǎng)格，網(wǎng)格數(shù)目為25萬；固體計算區(qū)域取巴氏合金層，采用類似油膜網(wǎng)格的劃分方法，網(wǎng)格數(shù)為7萬，見圖4。采用非結(jié)構網(wǎng)格劃分推力瓦結(jié)構。

圖4 固體計算網(wǎng)格圖Figure 4 Solid compute grid diagram

對結(jié)構體進行部分計算時，需要對推力瓦的各個方向添加約束條件。約束如圖5所示：A為對推力瓦頂部的固定約束（fixed support），固定住推力瓦瓦面避免產(chǎn)生位移。B處為重力作用約束（standard earth gravity），推力瓦在工作時受到地心引力的作用，設置重力加速度為9.8066m/s2。

圖5 推力瓦結(jié)構圖Figure 5 Thrust tile structure drawing

在數(shù)值模擬計算時對推力瓦的工作進行了簡化處理，認為每塊推力瓦的瓦面油膜分布情況與運動規(guī)律相同，且油流過軸瓦與鏡板之間對推力瓦的作用具有周期性，不同方位上推力瓦上的油膜分布情況與運動規(guī)律基本相同。本文選取一塊推力瓦作為研究對象，探究不同轉(zhuǎn)速工況以及不同油膜厚度情況下油膜對推力瓦產(chǎn)生的影響。

2.4 安全系數(shù)理論的定義

在機組運行時，推力軸承承受交變應力，因此，采用名義應力法進行分析。疲勞安全系數(shù)采用式（10）：

式中：nσ為疲勞極限；KσD為零件疲勞降低系數(shù)；σa為應力幅；[n]為許用安全系數(shù)。

3 計算結(jié)果及分析

3.1 不同轉(zhuǎn)速下的計算結(jié)果與分析

選取4種不同的機組轉(zhuǎn)速值，分別為100、350、500、725r/min（飛逸轉(zhuǎn)速），此時取得油膜厚度為0.066mm。分別對這幾個轉(zhuǎn)速工況進行計算，結(jié)果如下。

3.1.1 不同轉(zhuǎn)速下推力瓦應變對比

通過流固耦合計算可以得到推力瓦在不同轉(zhuǎn)速的總變形分布情況。圖6為機組在4個不同轉(zhuǎn)速下推力瓦固體部分總變形圖，可以看出當轉(zhuǎn)速為100r/min時［見圖6（a）］，推力瓦最大形變?yōu)?.25×10-3mm，隨著轉(zhuǎn)速不斷增加［見圖6（b）～圖6（d）］，推力瓦的形變不斷增加至2.22×10-2mm［見圖6（d）］，總體呈現(xiàn)形變增大的趨勢。對于推力瓦表面形變分布，對比4種工況最大形變位置幾乎沒有變化并且位于推力瓦中心處，從瓦面中心沿周向逐漸減小。可見，推力瓦的形變隨轉(zhuǎn)速升高而增加，同時最大形變位于推力瓦中心處。

圖6 瓦面應變分布（單位：m）Figure 6 Thrust tile surface strain distribution（unit：m）

3.1.2 不同轉(zhuǎn)速下推力瓦安全系數(shù)對比

圖7為不同轉(zhuǎn)速下推力瓦的安全系數(shù)，對比這4種轉(zhuǎn)速工況下，轉(zhuǎn)速為100r/min時［見圖7（a）］，此時的推力瓦安全系數(shù)較高，但是再隨著轉(zhuǎn)速升高，推力瓦瓦面開始出現(xiàn)安全系數(shù)下降的情況［見圖7（b）～圖7（d）］。隨著轉(zhuǎn)速升高，瓦面低安全系數(shù)區(qū)域開始從瓦面中心向周圍不斷擴大，并且逐漸向推力瓦上側(cè)靠近［見圖7（d）］。隨著轉(zhuǎn)速升高，瓦面中心安全系數(shù)降低。結(jié)合圖6可以看出，安全系數(shù)和瓦面應變分布相對應。由于轉(zhuǎn)速升高推力瓦產(chǎn)生的形變逐漸變大，因此，瓦面中心的安全系數(shù)不斷降低。為提高推力瓦的安全系數(shù)，應考慮選擇合適的金屬材料，減小推力瓦的形變，同時也應避免機組出現(xiàn)飛逸工況，以免出現(xiàn)推力瓦發(fā)生形變較大，安全系數(shù)降低的情況。

圖7 瓦面安全系數(shù)分布Figure 7 Thrust tile surface safety coefficient distribution

3.2 不同油膜厚度下的計算結(jié)果與分析

在轉(zhuǎn)速為500 r/min的情況下取4個工況點，即分別取油膜厚度為0.066、0.056、0.076、0.086mm，并對推力瓦在每個運行工況點下進行有限元求解，得到各個工況下推力瓦的應變及安全系數(shù)分布云圖。

3.2.1 不同油膜厚度下推力瓦應變對比

從圖8中可看出，當油膜厚度為0.056mm和0.066mm時，可以看出此時推力瓦受力不均勻［見圖8（a）、圖8（b）］；在這兩種工況下，形變主要集中在瓦面中心處。當厚度為0.076mm時［見圖8（c）］，發(fā)生最大形變的區(qū)域有所增大且向瓦面上側(cè)靠近，其余區(qū)域沿最大形變中心向兩邊逐漸減小。當油膜厚度增加到0.086mm時［見圖8（d）］，推力瓦應變分布與圖6（a）、圖6（b）相似，瓦面中心發(fā)生最大形變，由中心向周圍發(fā)生的應變呈梯度減小。從總體趨勢上看，隨著油膜厚度的增加，最大形變的值呈現(xiàn)先增加后急劇減少的趨勢。為減少推力瓦的形變對機組運行的影響，則需要維持合適的油膜厚度，盡量避免推力瓦出現(xiàn)較大形變。

圖8 瓦面應變分布（單位：m）（一）Figure 8 Strain distribution on tile surface（unit：m）（No.1）

圖8 瓦面應變分布（單位：m）（二）Figure 8 Strain distribution on tile surface（unit：m）（No.2）

3.2.2 不同油膜厚度下推力瓦安全系數(shù)對比

圖9為不同油膜厚度下推力瓦安全系數(shù)，可以看出4個不同工況下，推力瓦的最小安全系數(shù)出現(xiàn)在瓦面中心處。從圖中可看出，隨著厚度增加到0.076mm，瓦面中心的安全系數(shù)不斷減?。垡妶D9（a）～圖9（c）］，安全系數(shù)低的區(qū)域不斷變大。當油膜厚度為0.086 mm時［見圖9（d）］，瓦面中心的安全系數(shù)有所增加?？傮w上看，瓦面中心的安全系數(shù)隨著厚度增加先減小后又有明顯的增加。對比圖8和圖9，在油膜厚度為0.076mm時，發(fā)生最大形變的區(qū)域最大且安全系數(shù)較低的區(qū)域也最大，在機組運行時，應避免出現(xiàn)油膜厚度為0.076mm的情況或者避免長時間在該工況下運行。

圖9 瓦面安全系數(shù)分布Figure 9 Thrust tile surface safety coefficient distribution

4 結(jié)論

本文基于流固耦合方法，建立了推力瓦模型，計算并分析了推力瓦在不同轉(zhuǎn)速以及不同油膜厚度下推力瓦應力及安全系數(shù)，得出以下結(jié)論：

（1）同一油膜厚度下，推力瓦的應變隨著機組轉(zhuǎn)速升高而增大；安全系數(shù)隨著轉(zhuǎn)速升高而降低，均在飛逸時出現(xiàn)最差工況。

（2）同一轉(zhuǎn)速下，推力瓦的應變隨著油膜厚度增加，最大形變的值呈現(xiàn)先增加后急劇減少的趨勢；安全系數(shù)隨著油膜厚度增加，先減小后又有明顯的增加，均在油膜厚度為0.076mm處出現(xiàn)了最大形變區(qū)域最大和安全系數(shù)較低區(qū)域最大的情況，機組應避免在油膜厚度為0.076mm處長時間運行。

（3）在不同工況下，油膜對推力瓦均有不同程度的影響，在機組運行過程中保持穩(wěn)定的轉(zhuǎn)速以及維持合適的油膜厚度，這將有利于提高推力瓦安全系數(shù)。