賴芬， 王鳳鳴， 朱相源， 常沛然， 李國君

(西安交通大學 熱流科學與工程教育部重點實驗室，陜西 西安 710049)

泵是把機械能轉化為流體壓力能的通用機械，其中離心泵使用最廣泛，特別是在水利水電工程中。但環(huán)境惡化引起的泥沙流失導致河流中含有泥沙，例如，中國的黃河泥沙含量為37.8 kg/m3[1]，美國的科羅拉多河泥沙含量為10 kg/m3[2]。水利水電工程中的離心泵一般按清水工況設計，但其內(nèi)部的實際流動為液固兩相流動。固體顆粒的存在不僅會降低離心泵的效率，還會對離心泵壁面造成侵蝕磨損。磨損將導致頻繁地維修和更換過流部件，造成輸運系統(tǒng)停機，輸運費用大幅度增加。Wilson等[3]指出大型礦場每小時的停機成本大約為10萬美元。因此，研究黃河和長江流域使用的含沙水離心泵內(nèi)的液固兩相流動及其磨損特性具有重要意義。

但固體顆粒引起的壁面磨損是個復雜過程，最終的磨損形態(tài)受很多因素的影響，包括顆粒特性、流場特性、壁面特性等。國內(nèi)外學者對壁面磨損的影響因素進行了大量的試驗研究。例如，Tian等[4]通過科里奧利磨損試驗臺測試了不同顆粒粒徑下壁面的磨損系數(shù)，結果表明顆粒形狀、粒徑分布等顆粒特性會對磨損系數(shù)產(chǎn)生重要的影響。陳紅生等[5]通過離心泵液固兩相流水力試驗，發(fā)現(xiàn)造成局部磨損的重要原因是葉輪出口附近的射流-尾流結構。Wiedenroth[6]采用超聲設備測試了4種葉輪葉片的磨損量，發(fā)現(xiàn)當顆粒撞擊角較大時，硬度高的壁面受到的磨損較嚴重。

隨著計算機技術的迅猛發(fā)展，數(shù)值模擬方法已成為研究各種物理現(xiàn)象的重要手段。目前關于離心泵液固兩相流的數(shù)值模擬及磨損預測，磨損模型主要選用Finnie[7]和Tabakoff[8]模型，研究對象主要針對離心泵的某一過流部件[1,9-11]。例如，Noon等[9]應用Finnie磨損模型預測了石灰漿輸送泵蝸殼的磨損形態(tài)，結果表明蝸殼的磨損率隨撞擊速度、質(zhì)量濃度、顆粒粒徑的增大而增大。劉娟等[10]應用Finnie磨損模型對低體積分數(shù)的離散相顆粒在離心泵中的運動規(guī)律及葉輪壁面的磨損特性進行了數(shù)值研究，發(fā)現(xiàn)液固相密度差距越大，固體顆粒的運動跟隨性越差，固體顆粒與過流表面發(fā)生碰撞的幾率增大，葉輪壁面的磨損強度增加。黃先北等[1]采用Tabakoff 磨損模型對不同泥沙及不同入口工況下離心泵葉輪壁面的磨損規(guī)律進行了探索，發(fā)現(xiàn)顆粒粒徑會顯著影響葉輪壁面的磨損形態(tài)和位置，顆粒在離心泵入口分布越均勻，壁面磨損越分散，磨損位置軸對稱性越明顯。何創(chuàng)新[11]采用Tabakoff 磨損模型對單級雙吸中開式離心泵葉輪壁面的磨損特性進行了分析，結果表明葉輪入口前的密封體導葉能有效地抑制葉片的磨損，合理的葉片型線設計可以顯著降低葉片壓力側的磨損。

然而，大部分磨損模型的建立和推導都是基于氣固兩相流試驗的。磨損模型能否準確預測液固兩相流下的壁面磨損特性需進一步驗證。Zhang等[12]研究表明侵蝕磨損研究中心(E/CRC)提出的磨損模型不僅適用于氣固兩相流動，而且在液固兩相流動中也能獲得準確的結果。Peng等[13]指出E/CRC磨損模型與其他磨損模型相比，在CFD軟件中更易執(zhí)行，模型中考慮了顆粒硬度和形狀因素，得到的數(shù)值結果與試驗結果更接近。但現(xiàn)有研究中未見采用E/CRC磨損模型預測離心泵壁面磨損特性的研究。因此，本文針對黃河和長江流域使用的含沙水離心泵展開了研究，建立了基于E/CRC磨損模型的離心泵壁面磨損特性預測數(shù)值方法，對各個區(qū)域的磨損形態(tài)進行了預測并對比了不同區(qū)域的最大和平均磨損率，分析了磨損率變化規(guī)律并預測了最大磨損率發(fā)生位置，另外還討論了顆粒粒徑及濃度對離心泵葉輪磨損特性的影響。

1 液固兩相流及磨損預測數(shù)學模型

本文采用的流動介質(zhì)為黃河、長江流域的含沙水，根據(jù)長江夏季泥沙濃度[14]及黃河多年平均泥沙濃度[15]，泥沙濃度大約為32 kg/m3，體積分數(shù)小于3%，泵內(nèi)流動為低濃度液固流。文獻[16]指出歐拉-拉格朗日方法的使用條件是離散相體積分數(shù)小于10%～12%。因此，本文符合歐拉-拉格朗日方法的使用條件，為了提高計算精確性，本文選用雙向耦合的歐拉-拉格朗日方法求解。計算時忽略顆粒與顆粒間的相互作用力，并假定液固相之間沒有質(zhì)量和能量交換。計算時將液體視為連續(xù)相，液體流場通過在歐拉坐標系中求解雷諾時均方程獲得；將固體顆粒視為離散相，固體顆粒運動通過在拉格朗日坐標系中求解顆粒軌跡方程獲得。液體湍流脈動引起的顆粒擴散采用隨機游走模型。顆粒撞擊壁面后發(fā)生的動量變化由Grant 和Tabakoff(G&T)[8]提出的顆粒碰撞反彈模型模擬，磨損率由E/CRC 提出的磨損模型進行計算，詳細的數(shù)學模型如下。

1.1 連續(xù)相控制方程

1) 質(zhì)量守恒方程：

(1)

2) 動量守恒方程：

(2)

(3)

式中：μt為湍流粘度；μt由標準k-ε模型確定[17-18];k為湍動能；δij為“Kronecker delta”符號。

Fi為動量交換源項分量，F(xiàn)i從以下方程獲得:

(4)

1.2 離散相控制方程

文獻[16]指出顆粒相動量源項主要來自拖曳力、重力、虛擬質(zhì)量力、壓力梯度力、熱泳力、布朗力、Saffman 升力和Magnus升力產(chǎn)生的動量。但熱泳力主要是由溫度梯度引起的，離心泵中溫度變化很小，因此，本文忽略了熱泳力；布朗力和Saffman 升力只有處理亞微觀粒子時才考慮，Magnus升力只有處理旋轉粒子時才考慮，本文所考慮的粒子不是亞微觀粒子，也不旋轉，因此，本文忽略布朗力、Saffman 升力和Magnus升力，只考慮拖曳力、重力、虛擬質(zhì)量力和壓力梯度力產(chǎn)生的動量源項。

(5)

(6)

式中：dp為顆粒直徑；μ為液體的分子粘度；Cd為拖曳系數(shù)，定義為:

(7)

式中：a1、a2、a3為由Morsi等[19]給定的常數(shù)，Re為相對雷諾數(shù)，定義為：

(8)

Fv為用于加速顆粒周邊流體的虛擬質(zhì)量力矢量，定義為：

(9)

Fp為液體中的壓力梯度力矢量，當固體顆粒通過高壓區(qū)時，它對顆粒軌跡有著重要的影響，定義為：

(10)

1.3 隨機游走模型

顆粒軌道方程中的流體速度采用瞬時速度來考慮顆粒的湍流擴散，并計算足夠多的代表性顆粒的軌跡來考慮湍流對顆粒的隨機性影響。計算時考慮顆粒與流體的離散渦之間的相互作用，流體的脈動速度假定為時間的分段常量函數(shù)，在渦團的特征生存時間內(nèi)脈動速度保持為常量，滿足高斯概率密度分布函數(shù)的隨機脈動速度u′、v′、w′為：

(11)

(12)

(13)

(14)

式中ζ為服從正態(tài)分布的隨機數(shù)。

渦團的特征生存時間定義為常量：

(15)

式中：TL為流體的拉格朗日積分時間尺度；k為湍動能；ε為湍動能耗散率。

1.4 G&T顆粒碰撞反彈模型

Peng等[13]對比了2種顆粒碰撞反彈模型和5種磨損模型預測的結果，結果表明G&T顆粒碰撞反彈模型結合E/CRC磨損模型預測的結果與試驗結果最接近，因此，本文選用G&T顆粒碰撞反彈模型。模型中，顆粒撞擊壁面后，其垂直于壁面切線方向的動量變化率為法向恢復系數(shù)Vn2/Vn1，其平行于壁面切線方向的動量變化率為切向恢復系數(shù)Vt2/Vt1，分別定義為：

Vn2/Vn1=0.993-1.76β+1.56β2-0.49β3

(16)

Vt2/Vt1=0.988-1.66β+2.11β2-0.67β3

(17)

式中：Vn、Vt分別表示顆粒撞擊速度的垂直分量和切線分量， m/s；下標1、2分別表示撞擊前和撞擊后；β為顆粒撞擊角， rad。

1.5 E/CRC磨損模型

顆粒撞擊壁面后，對壁面造成的磨損與壁面材料、顆粒特性、撞擊角等因素相關。與其他磨損模型相比，E/CRC磨損模型考慮了顆粒硬度和形狀因素，得到的數(shù)值結果與試驗結果更接近。E/CRC磨損模型中磨損率計算方程為[12]：

(18)

(19)

式中：ER為磨損率(每單位面積的質(zhì)量損失)，kg/m2；C和n分別取值2.17×10-7，2.41；BH為布氏硬度;FS為顆粒形狀系數(shù)，對于球形顆粒取值為0.2；Vp為顆粒速度，m/s;β為顆粒撞擊角， rad；Ai值見表1。

表1 Ai值Table 1 Values of Ai

1.6 數(shù)學模型驗證

由于文獻中關于離心泵壁面磨損率的試驗數(shù)據(jù)較少，因此，本文采用直徑為50 mm、曲率半徑為76.9 mm的 90°彎管試驗對數(shù)學模型進行驗證。雖然離心泵涉及旋轉部件，整體結構與彎管不一致，運動條件也有差別，但由磨損模型可知，磨損率主要與顆粒形狀、硬度、速度、撞擊角和局部壁面特性有關，與是否涉及旋轉部件和設備整體結構無關。運動條件差別將導致磨損模型中的撞擊角和撞擊速度差別，也就是說磨損模型中考慮了運動條件的差別。因此，彎管試驗是可以驗證數(shù)值模擬方法的。

彎管壁面磨損率的數(shù)值結果與Zeng等[20]采用電化學方法測量獲得的試驗結果對比如圖1所示。由圖1可知，沿著彎管曲率角逆時針方向，磨損率逐漸增加。彎管磨損率數(shù)值結果與試驗結果吻合良好，最大磨損率與最大磨損位置預測準確。因此，該數(shù)學模型可以準確地預測液固兩相流中壁面的磨損率與最大磨損位置。

圖1 彎管磨損率數(shù)值結果與試驗結果對比Fig.1 Comparison of elbow erosion rate between numerical results and experimental data

2 模型泵及數(shù)值設定

2.1 模型泵及數(shù)值計算域

本文研究對象為單級單吸離心泵，主要幾何參數(shù)如表2所示，其設計參數(shù)為：流量25 m3/h，揚程15 m，轉速2 500 r/min，比轉速135。數(shù)值計算域包括進口延伸段、葉輪、無葉擴壓器、蝸殼及出口延伸段，其中進口延伸段長度為葉輪進口直徑的3倍，出口延伸段長度為蝸殼出口直徑的3倍。

表2 模型泵主要幾何參數(shù)Table 2 Main geometric parameters of model pump

2.2 網(wǎng)格劃分

為了生成高質(zhì)量的網(wǎng)格，采用ICEM-CFD里的結構化六面體網(wǎng)格對數(shù)值計算域進行劃分。為了準確地捕捉近壁面湍流，在近壁面采用加密處理，并調(diào)整第1層網(wǎng)格高度，將壁面y+值控制在30附近，最終的網(wǎng)格劃分結果如圖2所示。

圖2 計算域網(wǎng)格劃分Fig.2 Grid generation of computational domain

2.3 數(shù)值設定

數(shù)值計算時選取25°清水為連續(xù)相，球形沙顆粒(SiO2)為離散相。采用滑移網(wǎng)格技術模擬轉子與定子間的相對運動，設置進口延伸段和擴散段與葉輪的交界面為滑移界面，葉輪計算域設在旋轉坐標系，其余計算域設在靜止坐標系。進口邊界條件設為流體速度，出口邊界條件設為自由出流。球形沙顆粒從進口處垂直于進口面恒定釋放，速度與流體速度一致。非穩(wěn)態(tài)計算時間步長設為0.000 1 s，即葉輪每旋轉1.5°求解一次，每個時間步長迭代次數(shù)設為20次。

2.4 網(wǎng)格無關性驗證

為了驗證網(wǎng)格無關性，選取7組結構化六面體網(wǎng)格對數(shù)值計算域進行劃分，并計算模型泵在設計工況下壁面的平均磨損率。平均磨損率隨網(wǎng)格數(shù)變化趨勢如圖3所示。由圖3可知，隨著網(wǎng)格數(shù)增加，平均磨損率呈減小趨勢，但當網(wǎng)格數(shù)增加至286萬后，平均磨損率的變化很小，在5%范圍內(nèi)，因此，本文選取模型泵的計算網(wǎng)格數(shù)為286萬。

圖3 網(wǎng)格無關性研究設計工況Fig.3 Grid independence study

3 計算結果與分析

為了探索離心泵磨損初期各個區(qū)域的磨損形態(tài)和磨損率以及其受顆粒粒徑及濃度的影響，本文對設計工況、不同顆粒粒徑工況、不同顆粒濃度工況下的模型泵進行了研究。由于磨損初期壁面磨損較輕，壁面變化引起的通流型線變化較小，對磨損發(fā)生的情況和磨損條件的影響較小，因此，本文忽略通流型線的變化對磨損的影響。

文獻[2]指出黃河砂中值粒徑為30 μm，多砂河常見粒徑為90 μm。雖然實際泵流動的兩相流顆粒的直徑不會是某一均勻直徑，但采用平均粒徑計算的結果與采用一定顆粒直徑范圍計算的結果差別很小。為了分析粒徑對磨損形態(tài)和磨損率的影響，本文設計工況下選取平均粒徑60 μm，選取顆粒濃度32 kg/m3；不同顆粒粒徑工況時保持其他參數(shù)不變，顆粒直徑設為30、90、120 μm；不同顆粒濃度工況時保持其他參數(shù)不變，顆粒濃度設為19.5、44.5、57 kg/m3；具體計算工況如表3所示。

表3 計算工況Table 3 Calculation conditions

3.1 設計工況下離心泵壁面的磨損特性

為了探索離心泵壁面的磨損特性，對設計工況下模型泵各個區(qū)域的磨損形態(tài)進行了預測并對比了不同區(qū)域的最大和平均磨損率，分析了葉片與后蓋板交界處磨損率隨葉片曲率角的變化規(guī)律并預測了最大磨損率發(fā)生位置。

圖4為模型泵各個區(qū)域磨損形態(tài)的時間演化。由圖4可知，與其他過流部件相比，蝸殼和出口延伸段的磨損較輕，葉輪的磨損較嚴重，尤其是葉片前緣及葉片壓力側尾緣附近區(qū)域。與葉輪前蓋板相比，葉輪后蓋板磨損區(qū)域面積更大。當磨損時間由0.5 s增至2.0 s，蝸殼及出口延伸段的磨損率變化較小，葉輪及擴壓段的磨損率及磨損區(qū)域面積增加明顯，尤其是葉片壓力側尾緣附近區(qū)域。葉片前緣的磨損區(qū)域面積增加不明顯。當磨損時間為2 s時，前蓋板靠近葉片壓力側尾緣區(qū)域、后蓋板中部區(qū)域、后蓋板外圍區(qū)域均遭受了嚴重的磨損。

圖5為不同區(qū)域的最大和平均磨損率對比圖。由圖5可知，最大磨損率最大值發(fā)生在葉片和后蓋板區(qū)域，為7.9×10-4kg/m2；平均磨損率最大值發(fā)生在后蓋板區(qū)域，為1.4×10-6kg/m2；最大和平均磨損率最小值均發(fā)生在進口延伸段，分別為1.0×10-8和 2.6×10-9kg/m2；最大和平均磨損率分布規(guī)律不一致，例如，前蓋板的最大磨損率低于葉片，但平均磨損率高于葉片；同一區(qū)域的最大和平均磨損率差距巨大，例如，葉片的最大磨損率是平均磨損率的2 079倍。由以上分析可以推測，葉輪后蓋板是磨損最嚴重的區(qū)域，最大磨損率發(fā)生在葉片與后蓋板交界處。

圖6為葉片與后蓋板交界處磨損率隨葉片曲率角的變化規(guī)律。由圖6可知，葉片壓力側與吸力側的磨損率隨葉片曲率角的變化規(guī)律不同，葉片吸力側的磨損率曲線較平緩，葉片壓力側的磨損率曲線波動較大，葉片前緣吸力側的磨損率大于壓力側，葉片尾緣壓力側的磨損率大于吸力側；但無論是葉片壓力側還是吸力側，磨損率均是在葉片前緣取得最大值，最大磨損率發(fā)生在葉片曲率角59.8°附近，為7.9×10-4kg/m2。

3.2 顆粒粒徑對葉輪磨損特性的影響

為了進一步探索離心泵壁面磨損特性的影響因素，對不同顆粒粒徑工況下(工況1、2、3)離心泵的磨損率進行了計算。由上述分析可知，葉輪是離心泵內(nèi)磨損最嚴重的過流部件，最大磨損率發(fā)生在葉片與后蓋板交界處。因此，這里分析顆粒粒徑對葉輪磨損形態(tài)及磨損率變化規(guī)律的影響，結果如圖7、圖8所示。圖7為不同顆粒粒徑工況葉輪磨損云圖變化，其中左邊為前蓋板，中間為葉片，右邊為后蓋板。由圖7可知，葉輪磨損形態(tài)受顆粒粒徑影響顯著。隨著顆粒粒徑的增大，磨損區(qū)域面積逐漸減小。這主要是由于在保持顆粒濃度不變的情況下，隨著顆粒粒徑的增大，顆粒數(shù)減小，葉輪壁面受到的撞擊次數(shù)減小。當顆粒粒徑為30 μm 時，葉輪前后蓋板中部及外圍均受到了不同程度的磨損，葉片前緣及葉片壓力側尾緣附近區(qū)域磨損嚴重；當顆粒粒徑增大至90 μm 時，前蓋板及葉片的磨損區(qū)域面積顯著減小，而后蓋板的磨損區(qū)域面積幾乎不變，但同一位置的磨損率顯著減?。划旑w粒粒徑增大至120 μm 時，前蓋板的磨損區(qū)域只有葉片壓力側尾緣附近區(qū)域，后蓋板的磨損區(qū)域仍然是中部及外圍。

圖8為不同顆粒粒徑工況葉片與后蓋板交界處磨損率隨葉片曲率角的變化規(guī)律。其中，圖8(a)為葉片壓力側，葉片曲率角在-22.5°～59.8 °變化；圖8(b)為葉片吸力側，葉片曲率角在-32.5°～59.8°變化。由圖8(a)可知，在各個粒徑工況下，葉片壓力側磨損率基本是隨葉片曲率角先減小后增大，最大峰值在葉片曲率角59.8°附近，第2峰值在葉片曲率角-22.5°附近；由圖8(b)可知，在各個粒徑工況下，葉片吸力側磨損率基本是隨葉片曲率角先小幅度變化到葉片曲率角55°附近迅速上升，在葉片曲率角59.8°附近達到最大值；隨著顆粒粒徑增大，葉片壓力側及吸力側磨損率總體上呈減小趨勢。

圖4 磨損形態(tài)的時間演化Fig.4 Time evolution of erosion pattern

圖5 設計工況t=10 s時不同區(qū)域的最大和平均磨損率對比Fig.5 Comparisons of the maximum and average erosion rates for different regions under design condition at t=10 s

3.3 顆粒濃度對葉輪磨損特性的影響

為了進一步探索離心泵壁面磨損特性的影響因素，對不同顆粒濃度工況下(工況4、5、6)離心泵的磨損率進行了計算并分析顆粒濃度對葉輪磨損形態(tài)及磨損率變化規(guī)律的影響，結果如圖9、圖10所示。圖9為不同顆粒濃度工況葉輪磨損云圖變化，其中左邊為前蓋板，中間為葉片，右邊為后蓋板。由圖9可知，葉輪磨損形態(tài)受顆粒濃度影響顯著。隨著顆粒濃度的增大，同一位置的磨損率顯著增大。這主要是由于在保持顆粒粒徑不變的情況下，隨著顆粒濃度的增大，顆粒數(shù)增大，葉輪壁面受到的撞擊次數(shù)增多。但葉輪各個部分的磨損區(qū)域面積幾乎不變，這說明改變顆粒濃度，不會改變顆粒撞擊角。當顆粒濃度為19.5 kg/m3時，嚴重磨損主要集中在葉片前緣附近區(qū)域；當顆粒濃度增大至44.5 kg/m3時，同一位置的磨損率增大，葉片壓力側尾緣附近變成嚴重磨損區(qū)域；當顆粒濃度增大至57 kg/m3時，同一位置的磨損率進一步增大，葉片壓力側尾緣附近嚴重磨損區(qū)域面積增加明顯，葉片前緣附近嚴重磨損區(qū)域面積變化較小。

圖6 磨損率隨葉片曲率角的變化規(guī)律(dp=60 μm,Cm=32 kg/m3)Fig.6 Variation of the erosion rate along the blade curvature angle(dp=60 μm,Cm=32 kg/m3)

圖7 不同顆粒粒徑葉輪磨損云圖Fig.7 Impeller erosion contours for different particle diameters

圖10為不同顆粒濃度工況葉片與后蓋板交界處磨損率隨葉片曲率角的變化規(guī)律。其中，圖10(a)為葉片壓力側，葉片曲率角在-22.5°～59.8°變化；圖10(b)為葉片吸力側，葉片曲率角在-32.5°～59.8°變化。由圖10(a)可知，在各個濃度工況下，葉片壓力側磨損率的分布基本上是葉片前緣尾緣大，中間區(qū)域小，在葉片曲率角59.8°附近達到最大值，在葉片曲率角-22.5°附近達到第2峰值；由圖10(b)可知，在各個濃度工況下，葉片吸力側磨損率的分布基本上是葉片前緣大其余區(qū)域小，也是在葉片曲率角59.8°附近達到最大值；隨著顆粒濃度增大，葉片壓力側及吸力側磨損率隨葉片曲率角的變化規(guī)律較一致，同一位置上的磨損率基本上呈增大趨勢。

圖8 不同顆粒粒徑下磨損率隨葉片曲率角的變化Fig.8 Variations of the erosion rates along the blade curvature angle for different particle diameters

圖9 不同顆粒濃度葉輪磨損云圖Fig.9 Impeller erosion contours for different particle concentrations

圖10 不同顆粒濃度下磨損率隨葉片曲率角的變化Fig.10 Variations of the erosion rates along the blade curvature angle for different particle concentrations

4 結論

1) 葉輪是離心泵內(nèi)磨損最嚴重的過流部件，最嚴重的磨損發(fā)生在葉片前緣及葉片壓力側尾緣附近區(qū)域。

2) 進口延伸段的最大和平均磨損率小于其它區(qū)域，葉片和后蓋板的最大磨損率大于其它區(qū)域，葉片與后蓋板交界處葉片曲率角59.8°附近的最大磨損率達到最大值。

3) 前蓋板及葉片的磨損區(qū)域面積隨粒徑增大顯著減小，而后蓋板的磨損區(qū)域面積隨粒徑增大幾乎不變，同一位置的磨損率隨粒徑增大呈減小趨勢。

4) 葉輪各個部分的磨損區(qū)域面積隨濃度增大幾乎不變，同一位置的磨損率隨濃度增大呈增大趨勢，葉片壓力側及吸力側磨損率隨濃度增大變化規(guī)律較一致。