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

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

0 引言

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

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

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

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

1.1 推力軸承基本理論

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

1.2 推力瓦三維模型

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

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

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

2 流固耦合分析

2.1 流固耦合基本理論

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

以及時(shí)均動(dòng)量方程[11]：

對(duì)應(yīng)的湍動(dòng)能κ以及湍流頻率ω方程分為式（3）和式（4）：

式中：下標(biāo)f為流體，s為固體；τ為流體與固體應(yīng)力；u為位移；q為熱流量；T為溫度。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3.1.1 不同轉(zhuǎn)速下推力瓦應(yīng)變對(duì)比

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

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

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

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

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

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

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

3.2.1 不同油膜厚度下推力瓦應(yīng)變對(duì)比

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

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

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

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

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

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

4 結(jié)論

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

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

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

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