張自超 李延頻 陳德新

(1.華北水利水電大學電力學院，鄭州 450045；2.華北水利水電大學河南省流體機械工程技術研究中心，鄭州 450045；3.華北水利水電大學烏拉爾學院，鄭州 450045)

0 引言

廣泛應用于黃河沿岸灌溉泵站中的雙吸離心泵，各部件存在明顯的磨損問題。葉輪作為高速旋轉的部件，其磨損問題更為嚴重，對性能影響更顯著。因此，研究雙吸離心泵葉輪的磨損特性具有重要的意義。

目前，數(shù)值模擬是研究水力機械泥沙磨損的主要手段。較多學者采用定常計算研究水力機械的磨損特性[1-8]。然而，泵內存在明顯的動靜干涉作用，必然導致固液兩相流的內流場具有非定常特性。如壓力脈動就是動靜干涉引起的，這說明動靜干涉對流場的影響不可忽視。

越來越多的學者關注了泵內固液兩相流的非定常特性[9-13]。然而，水泵的磨損與固相速度、固相濃度和固相與壁面的碰撞角等流場參數(shù)有關，因此，水泵的磨損也必然具有非定常特性。

葉輪作為泵內受動靜干涉作用影響最大的部件，其磨損的非定常特性必然顯著。而現(xiàn)有研究關注較多的是揚程、功率、壓力脈動等的非定常特性，并未對磨損的非定常特性進行過多研究，尤其是對雙吸離心泵葉輪內的泥沙磨損非定常特性的相關研究較少。本文以雙吸離心泵為研究對象，采用歐拉-歐拉方法對其內部固液兩相流進行非定常計算，研究不同工況下，葉輪壁面上多個監(jiān)測點的固相速度、固相濃度和磨損率的非定常特性，并與定常計算結果進行對比。

1 歐拉-歐拉方法

1.1 控制方程

歐拉-歐拉方法的控制方程為

(1)

(2)

其中

αl+αs=1

式中k——相類別(l為液相，s為固相)

α——體積分數(shù)μt——湍動粘度

v——時均速度p——壓強

Fi——相間作用力gi——體積力

κ——泥沙擴散系數(shù)ρ——密度

t——時間μk——動力粘度

λ——應力xi、xj——坐標分量

vki、vkj——速度梯度

其中，相間阻力的計算體現(xiàn)在式(2)的相間作用力項Fi中，而含泥沙擴散系數(shù)項是根據(jù)梯度輸移假定，通過泥沙擴散系數(shù)建立各相濃度、速度脈動關聯(lián)項與濃度梯度的關系得到的，具體可參考文獻[14]。

1.2 磨損模型

圖爾薩大學的磨損模型是目前泥沙兩相流磨損計算的經(jīng)典模型[15]，考慮了顆粒沖擊速度、沖擊角和壁面材料的硬度對磨損的影響，即

(3)

其中

F(θ)=5.4θ-10.11θ2+10.93θ3-6.33θ4+1.42θ5

(4)

ER=WRmsand/Acell

(5)

式中WR——磨損失重率，壁面磨損失去的質量與顆粒質量的比值，kg/kg

BH——壁面材料的布氏硬度

Fs——顆粒形狀系數(shù)

vs——顆粒入射速度，m/s

θ——顆粒沖擊角，rad

a、C——經(jīng)驗常數(shù)，取2.41、2.17×10-7

ER——磨損率，kg/(m2·s)

msand——固相質量流量，kg/s

Acell——計算單元壁面網(wǎng)格面積

棱角顆粒Fs=1，半球形顆粒Fs=0.53，球形顆粒Fs=0.2，本文取Fs=0.2。

2 計算對象和設置

2.1 計算對象

本文的研究對象為山西尊村引黃泵站的雙吸離心泵，主要幾何參數(shù)：額定流量Qn為3 m3/s，額定轉速n為490 r/min，設計揚程H為32 m，比轉數(shù)ns為162，葉輪直徑D為1.12 m，葉片數(shù)z為6，泵進口直徑D1為1.1 m，泵出口直徑D2為0.9 m。

計算含沙條件為：含沙質量濃度5 kg/m3、體積分數(shù)0.188%時，顆粒密度2 650 kg/m3；沙粒粒徑ds為25 μm的懸移質典型粒徑。

2.2 網(wǎng)格及邊界條件

雙吸離心泵的計算網(wǎng)格由Gambit生成，經(jīng)過網(wǎng)格無關性檢查，最終確定各部分網(wǎng)格單元數(shù)分別為：葉輪603 620；吸水室564 882；壓水室400 447；葉輪前后腔水體64 288，口環(huán)間隙較小，計算復雜，而且對葉片磨損影響不大，因此本文不考慮口環(huán)間隙；進口延長段44 100；出口延長段56 202；總計網(wǎng)格數(shù)為1 733 539。網(wǎng)格無關性檢查如圖1所示，雙吸離心泵計算域及網(wǎng)格劃分情況如圖2所示。

圖1 網(wǎng)格無關性檢查Fig.1 Grid-independent verification

圖2 計算域及網(wǎng)格劃分Fig.2 Computational domain and grid

應用Fluent進行數(shù)值模擬，運用歐拉-歐拉方法對雙吸泵進行固液兩相流數(shù)值計算。采用Phase Coupled SIMPLE算法求解二階迎風格式的離散方程。湍流模型采用RNGk-ε模型，因為該模型可以更好地處理高應變率及流線彎曲程度較大的流動[16]；考慮到相間阻力和滑移速度的影響，湍流多相流模型采用Dispersed turbulence模型，該模型是目前應用較多的一種多相流湍流模型[17-18]；泥沙擴散系數(shù)模型采用Diffusion-in-VOF模型進行計算。壁面磨損計算采用圖爾薩大學的磨損模型[15]。

計算域進口采用速度進口，分別給定固、液兩相的速度，并給定固相體積分數(shù)；出口采用自由出流；過流部件內壁面，不考慮粗糙度，設為光滑壁面，對液相采用無滑移壁面邊界條件，對固相采用自由滑移壁面邊界條件，近壁區(qū)采用標準壁面函數(shù)。

以定常計算的結果為初始條件進行非定常計算，共計算10個葉輪旋轉周期，每個旋轉周期計算120個時間步，即每旋轉3°為一個時間步。

2.3 阻力模型選擇

相間阻力模型是固液兩相流計算中重要的相間作用因素，合理地選擇阻力模型可以提高固液兩相流的計算精度。

常用的固液相間阻力模型有Syamlal-O’Brien模型[19]、Wen-Yu模型[20]、Gidaspow模型[21]和Huilin-Gidaspow模型[22]等。

Syamlal-O’Brien模型是基于流化床中顆粒的沉降得到的，因為模型中考慮了固體剪切粘度參數(shù)，隨著固相濃度趨于零，固相剪切粘度也趨于零[23],并不適用于稀疏固液兩相流的計算[24]。Wen-Yu模型被推薦用于固相體積分數(shù)小于20%的情況[17,25]，而Gidaspow模型和Huilin-Gidaspow模型是Wen-Yu模型和Ergun模型的組合，這兩種模型可以用于更高固相濃度的兩相流計算。但是，本文含沙體積分數(shù)遠小于20%，此時，Wen-Yu模型、Gidaspow模型和Huilin-Gidaspow模型的本質均為Wen-Yu模型，尤其是在固相體積分數(shù)小于1%的情況下。因此，本文選用Wen-Yu模型作為相間阻力模型進行計算。

2.4 葉片監(jiān)測點布置

為了監(jiān)測葉輪內部固相參數(shù)的脈動情況，需在葉片壁面上設置監(jiān)測點，各監(jiān)測點位置如圖3所示。

圖3 葉片上監(jiān)測點設置示意圖Fig.3 Monitoring locations on blade

分別在葉片吸力面和壓力面上各設置3個監(jiān)測點，沿葉片吸力面中間流線從進口到出口分別設置3個監(jiān)測點，分別為：S1位于葉片吸力面進口處，S2位于葉片吸力面中間位置，S3位于葉片吸力面出口處，同樣對應壓力面設置3個監(jiān)測點P1、P2、P3。6個監(jiān)測點位于葉片上并隨葉片一起旋轉。

各監(jiān)測點監(jiān)測的固相參數(shù)為固相速度、固相濃度和監(jiān)測點處的壁面磨損率和沖擊角，各監(jiān)測值只記錄參數(shù)的大小，不監(jiān)測方向。

3 數(shù)值計算驗證

3.1 外特性驗證

為了驗證數(shù)值計算的合理性，將數(shù)值計算得到的雙吸離心泵揚程和效率與試驗值作對比，并分析誤差，如圖4所示。計算工況分別為0.4Qn、0.6Qn、0.8Qn、Qn、1.2Qn。

圖4 雙吸離心泵外特性計算值和試驗值對比Fig.4 Comparison of calculated value and experimental value for external characteristics of double suction centrifugal pump

由圖4可知，計算得到的效率和揚程與試驗值均較為接近，但有一定的誤差。揚程的計算值高于試驗值，最大誤差為0.8Qn處的3.6%；而效率計算值較接近試驗值，最大誤差為Qn處的2.1%。因此，外特性計算值與試驗值較吻合，計算誤差在可接受的范圍內，計算結果合理、可靠，可用于進一步的分析和研究。

3.2 固液兩相流固相體積分數(shù)分布驗證

為了進一步對固液兩相流流場進行驗證，選取文獻[26]中的一個圓管為計算對象，如圖5所示。采用與上述雙吸離心泵完全相同的計算模型進行固液兩相流計算，將計算得到的固相體積分數(shù)分布與試驗值進行對比驗證。

圖5 計算圓管示意圖Fig.5 Sketch of calculated circular pipe

圓管直徑Dp為0.15 m，長度L為3 m。進口流速vl為3 m/s，進口平均含沙體積分數(shù)αvf為15%。顆粒粒徑ds為90 μm。

圓管水平布置，坐標原點位于圓管底部，y軸豎直向上，線1為驗證固相體積分數(shù)分布的位置，距離進口3L/4且豎直向上。

計算結果如圖6所示。圖中y*=y/Dp，αs為圓管中線1上各點的固相體積分數(shù)，αvf為圓管進口平均固相體積分數(shù)。

圖6 固相體積分數(shù)計算值和試驗值對比Fig.6 Calculated solid concentration distribution compared with experimental data

由圖6可知，本文采用的兩相流計算模型得到的固相體積分數(shù)的計算值與試驗值分布趨勢一致，結果較為吻合，相差不大，最大誤差為9.5%，計算誤差在可接受的范圍內。由此證明，本文采用的固液兩相流計算模型合理、可靠，可用于流場的進一步計算和分析。

4 結果分析

4.1 葉輪磨損特性

4.1.1葉片上磨損率脈動特性分析

圖7給出了各監(jiān)測點不同工況時磨損率非定常和定常計算結果對比，圖中，T為葉輪旋轉周期。

由圖7可知，葉輪表面各監(jiān)測點的磨損率具有周期性，其波動周期等于葉輪旋轉周期，這與動靜干涉引起的葉輪內固液兩相流流場具有脈動特性有關，這說明葉輪壁面的磨損具有脈動特性。同時，最大磨損率與最小磨損率之間相差較大，脈動特性明顯。隨著流量的增大，磨損率逐漸增大，脈動更加劇烈。在葉片吸力面，葉片進口處的磨損率較小，出口處的磨損率較大；而在葉片壓力面，葉片中部的磨損率較小，最大磨損率在葉片出口處。葉片壓力面出口處為葉輪磨損最大位置，而壓力面中部的磨損率最小，這可能與該處流態(tài)變化有關。

圖7 各監(jiān)測點不同工況時磨損率非定常和定常結果對比Fig.7 Comparisons for unsteady and steady results of erosion rate under different operation conditions at each monitoring point

由圖7可以看出，定常計算得到的磨損率遠小于非定常的結果，為了更清楚地對比二者的計算結果，表1給出了各監(jiān)測點不同工況時定常計算得到磨損率和非定常計算得到的最大磨損率的對比。由表1可以看出，定常計算得到的磨損率遠小于非定常的結果，這說明采用定常方法計算雙吸離心泵的壁面磨損率時，誤差較大。而隨著流量的增大，定常和非定常計算得到磨損率的誤差增大。這說明定常計算并不能準確反映葉輪壁面上磨損的脈動特性。

表1 定常和非定常計算磨損率對比Tab.1 Comparison of unsteady and steady results of erosion rate kg/(m2·s)

4.1.2葉片上磨損分布特性分析

為了分析葉片表面上的磨損分布情況，圖8給出了不同工況下，葉片吸力面和壓力面的磨損云圖。

為了對計算得到的磨損分布進行驗證，圖9給出了離心泵葉片表面的實際磨損情況，其中圖9a為文獻[27-28]得到的離心泵葉片表面磨損形貌試驗結果示意圖，圖9b為雙吸離心泵葉片現(xiàn)場磨損圖。由于本文研究的雙吸離心泵缺少磨損試驗，圖9b為同類型其它的雙吸離心泵的實際磨損情況，根據(jù)磨損規(guī)律，其磨損情況可以定性地驗證本文的數(shù)值計算得到的磨損分布情況。

由圖8可知，葉片表面的磨損率隨著流量的增大而增大，磨損面積也逐漸增大。各工況下，磨損最大的部位為頭部和尾部，這與圖7中非定常計算結果一致。吸力面中部的磨損面積大于壓力面，而在頭部和尾部，壓力面的磨損大于吸力面，這與圖9a的試驗結果一致。同時根據(jù)圖7得到的結論和圖8的分布情況，可以看出，壓力面尾部磨損最為嚴重，這與圖9b中的雙吸離心泵葉片現(xiàn)場磨損情況一致。同時也說明本文的計算結果可靠。

圖8 不同工況下葉片表面的磨損分布云圖Fig.8 Calculated erosion distributions of blade surface under different operation conditions

圖9 離心泵葉片實際磨損情況Fig.9 Actual erosion conditions of centrifugal pump blade

4.2 固液兩相流場特性分析

由式(3)～(5)可知，對于確定的計算對象和條件，壁面材料硬度、顆粒形狀系數(shù)等均為定值，影響磨損率的主要流場參數(shù)為沖擊角、固相體積分數(shù)和固相速度。為了進一步分析葉輪壁面磨損特性及兩相流場特性，對各監(jiān)測點在不同工況時的沖擊角、固相體積分數(shù)和固相速度的脈動特性進行分析。

4.2.1沖擊角及沖擊角函數(shù)特性

沖擊角對磨損率的影響，通過沖擊角函數(shù)F(θ)來反映，沖擊角函數(shù)F(θ)的具體形式見式(4)。圖10為壁面沖擊角的示意圖，沖擊角為固相速度與壁面切線的夾角，其大小為0°～90°，由圖可知，固相顆粒在撞擊壁面之后反射射出。

圖10 壁面沖擊角示意圖Fig.10 Diagram of wall impact angle

圖11為沖擊角函數(shù)隨沖擊角變化規(guī)律，由圖可知，隨著沖擊角的增大，沖擊角函數(shù)先增大后減小，其最大值出現(xiàn)在50°附近，這說明存在一個使磨損最大的沖擊角，此時的切削作用最顯著。

圖11 沖擊角函數(shù)變化曲線Fig.11 Law of impact angle function

當沖擊角為20°～80°時，沖擊角函數(shù)大于1，說明沖擊角增強了磨損損失；當沖擊角小于20°或大于80°時，沖擊角函數(shù)小于1，說明沖擊角減弱了磨損損失；當沖擊角達到50°附近時，沖擊角函數(shù)最大，沖擊角對磨損損失的影響最大。

圖12給出了各監(jiān)測點不同工況時，沖擊角及其函數(shù)的波動規(guī)律。由圖12可知，沖擊角和沖擊角函數(shù)隨葉輪旋轉呈周期性變化，其脈動曲線與圖7中磨損率的脈動曲線較相似，這說明沖擊角對磨損的影響至關重要，其決定了磨損率的脈動特性，以及最大磨損率和最小磨損率。沖擊角具有增強或減弱磨損的作用。同一監(jiān)測點處，不同工況下，沖擊角大小有所不同，這與工況改變后，流場的流態(tài)有所改變，導致沖擊角發(fā)生變化有關。

圖12 各監(jiān)測點不同工況時沖擊角及其函數(shù)波動規(guī)律對比Fig.12 Comparisons for fluctuation of impact angle and its function under different operation conditions at each monitoring point

由此可知，沖擊角對磨損的影響為增強或減弱磨損效果的作用，其決定了磨損率的脈動特性及最大磨損率、最小磨損率。

4.2.2固相體積分數(shù)特性

圖13給出了各監(jiān)測點不同工況時，固相體積分數(shù)的變化規(guī)律。由圖13可知，固相體積分數(shù)隨著葉輪旋轉也具有周期性。同一監(jiān)測點處，隨著流量的增大，固相體積分數(shù)有所增加；而隨著流量減小，固相體積分數(shù)隨時間的波動更為顯著，其原因可能是小流量時，偏離最優(yōu)工況，流態(tài)變差，二次流和漩渦增多，流場波動明顯。相同工況下，吸力面和壓力面對應監(jiān)測點的固相體積分數(shù)分布相差不大。吸力面和壓力面從進口到出口，固相體積分數(shù)分布波動逐漸減小，這是因為葉片進口處的流態(tài)較為復雜，容易出現(xiàn)二次流和漩渦，在葉片出口處，流態(tài)基本較為穩(wěn)定。

圖13 各監(jiān)測點不同工況時固相體積分數(shù)波動規(guī)律對比Fig.13 Comparison for solid volume fraction under different operation conditions at each monitoring point

通過與圖7對比，發(fā)現(xiàn)固相體積分數(shù)的脈動曲線與磨損率脈動曲線差別較大，這說明固相體積分數(shù)并不能決定磨損率的脈動特性，只是對其大小有所影響。

4.2.3固相速度特性

圖14給出了各監(jiān)測點不同工況時，固相速度的變化規(guī)律。由圖可知，隨著葉輪旋轉，固相速度的變化具有周期性。同一監(jiān)測點處，隨著流量增大，固相速度逐漸增大。通過對比圖7可以發(fā)現(xiàn)，最大固相速度出現(xiàn)在葉片壓力面出口處(P3)，最小速度位于壓力面中間(P2)，這與圖7中磨損率最大和最小的位置完全相同，這說明圖7中壓力面中間(P2)處出現(xiàn)磨損率最小主要是由該處的固相速度最小引起的，而葉片壓力面出口(P3)處磨損率最大主要是由該處固相速度最大引起的。

圖14 各監(jiān)測點不同工況時固相速度波動規(guī)律對比Fig.14 Comparisons for solid velocity under different operation conditions at each monitoring point

對比圖14和圖7可以看出，磨損率與固相速度的脈動規(guī)律非常吻合，這說明固相速度對磨損大小的影響更為顯著。

5 結論

(1)葉輪表面磨損率分布具有周期性，其脈動周期為葉輪旋轉周期，其脈動特性明顯。隨著流量增大，磨損率增大，脈動特性更明顯。葉片壓力面出口處的磨損率最大。

(2)定常計算的磨損率遠小于非定常結果，隨著流量增大，二者誤差增大。因此，定常計算結果不能準確預測葉輪壁面上的磨損損失。

(3)磨損最大的部位為頭部和尾部。吸力面中部的磨損面積大于壓力面，而在頭部和尾部，壓力面的磨損大于吸力面。

(4)葉輪壁面沖擊角、沖擊角函數(shù)、固相體積分數(shù)和固相速度分布均具有周期性。沖擊角函數(shù)的脈動曲線與磨損率脈動曲線較相似，沖擊角對磨損損失具有增強或減弱的作用，存在一個使磨損率最大的沖擊角，其決定了磨損率的脈動特性。

(5)在小流量時，固相體積分數(shù)脈動明顯。固相體積分數(shù)的脈動曲線與磨損率脈動曲線差別較大，其對磨損率脈動特性影響較小，對磨損大小有所影響。固相速度脈動規(guī)律及固相速度最大、最小值所處位置均與磨損率脈動特性相似，固相速度對磨損率影響顯著。