俞蕓蕓,周大慶,于安,劉佳佳

(1. 南京航空航天大學(xué)金城學(xué)院，江蘇 南京 211100；2. 河海大學(xué)能源與電氣學(xué)院，江蘇 南京 211100; 3. 南水北調(diào)東線江蘇水源有限責(zé)任公司揚州分公司, 江蘇 揚州 225000)

軸流泵是一種低揚程、大流量且應(yīng)用廣泛的泵型，近年來其受到越來越多的關(guān)注.空化是評價泵性能的三大指標(biāo)之一，與效率、穩(wěn)定性的研究相比，空化的研究還不是十分的充分，空化問題制約著現(xiàn)代水泵行業(yè)的發(fā)展，也是國內(nèi)外學(xué)者研究的熱點之一[1].軸流泵空化特性的研究一般采用數(shù)值模擬和試驗驗證等方法.隨著計算機科學(xué)的快速發(fā)展，CFD技術(shù)被廣泛應(yīng)用于數(shù)值模擬中，國內(nèi)外學(xué)者對算法的改進也使得數(shù)值模擬的精確度越來越高，數(shù)值模擬已經(jīng)成為研究空化現(xiàn)象的重要手段之一[2].

近年來，用于空化計算的湍流模型得到了進一步修正，其中一種由標(biāo)準(zhǔn)k-ε模型改進而來的局部時均化模型(partially-averaged Navier-Stokes model, PANS模型)被越來越多地應(yīng)用于空化流動計算中[3-4].GIRIMAJI等[5]認為PANS模型是雷諾平均RANS和直接數(shù)值模擬DNS之間的橋梁模型，并提出其控制參數(shù)是未分解湍動能比值fk以及未分解耗散率比值fε，并對該模型進行了初步驗證.LAKSHMIPATHY等[6]利用PANS模型對圓柱繞流進行了模擬，并采用不同的參數(shù)fk和fε值研究了PANS模型的可實施性和準(zhǔn)確性.羅大海等[7]采用PANS湍流模型對超聲速斜面空腔流動進行了數(shù)值模擬，并取定值fε=1，而fk值隨流場中的時間和空間變化，結(jié)果表明，相比于取定值，采用動態(tài)的fk可以更加準(zhǔn)確地模擬出三維流動的結(jié)構(gòu).黃彪等[8]采用基于k-ε模型改進的PANS湍流模型進行水翼的空化流動計算，證明了fk值對空化計算的精度有較大影響.施衛(wèi)東等[9]采用PANS模型對繞二維Clark-Y型水翼進行了非定常計算，通過對fk取不同的值研究了其對非定?？栈挠绊?石磊等[10]采用PANS模型對軸流泵葉頂間隙空化進行了數(shù)值模擬和試驗驗證，準(zhǔn)確地模擬出了葉頂區(qū)的空化特性.HUANG等[11]采用修正的PANS模型對繞臺階流模型進行數(shù)值計算，證明了其在流場壓力、速度分布和旋渦形態(tài)上具有更準(zhǔn)確的預(yù)測能力.

基于前人研究，文中利用計算軟件CFX二次開發(fā)技術(shù)對標(biāo)準(zhǔn)PANS模型中的fk值(未分解湍動能占總湍動能的比值)進行動態(tài)定義，形成修正的PANS模型.利用該湍流模型對軸流泵全流道進行空化數(shù)值模擬和試驗測試，驗證修正PANS模型在軸流泵空化特性模擬上的可靠性.

1 數(shù)值計算理論及方法

1.1 基本控制方程

CFD計算方法中涉及的流體動力學(xué)控制方程是基于質(zhì)量守恒定律、動量守恒定律、能量守恒定律形成的3個控制方程，根據(jù)常規(guī)水力機械運行條件，流動介質(zhì)為水，且常溫常壓，不考慮熱量交換，一般認為其內(nèi)部流場為等溫不可壓縮的三維流動，即不涉及能量交換，只包含質(zhì)量和動量守恒定律[12]，因此，基本方程包括連續(xù)性方程和Navier-Stockes方程.

連續(xù)性方程[13]

(1)

Navier-Stockes方程[14]

(2)

式中：ρ為密度；p為流體壓力；ui，uj為不同方向流速的瞬時值；μ為動力黏性系數(shù)；μt為湍流黏性系數(shù).

1.2 湍流模型修正

文中采用的湍流模型是修正的PANS模型，即基于標(biāo)準(zhǔn)PANS模型，利用CFX的二次開發(fā)技術(shù)編寫CCL語言，對PANS模型中的fk值進行動態(tài)定義.

由于軸流泵內(nèi)湍流流動的高度復(fù)雜性，在數(shù)值模擬計算中需要對湍流模擬方法進行研究和分析，根據(jù)不同的情況選擇不同的湍流模型進行數(shù)值模擬，以取得與實際流態(tài)更為吻合的結(jié)果.數(shù)值模擬方法分為直接模擬(DNS)和非直接模擬，非直接模擬應(yīng)用最為廣泛的是雷諾時均法(RANS)，雷諾時均法的核心是如何準(zhǔn)確合理地求解時均化的雷諾方程，從而取代直接求解瞬時Navier-Stokes方程的方法[15].常見的雷諾時均法包括零方程模型、一方程模型、k-ε模型、k-ω模型等.文中采用的修正PANS模型是基于k-ε模型改變而來，因此對該模型進行詳細介紹.

k-ε模型在計算時不直接計算雷諾應(yīng)力，而是通過湍動黏度(渦黏系數(shù))對湍流應(yīng)力進行參數(shù)表達，即根據(jù)Boussinesq渦黏假定得出[16]

(3)

式中：δij的取值方式中，當(dāng)i=j時，δij=1；當(dāng)i≠j時，δij=0；k為湍動能，即

(4)

由式(3)可知，k-ε模型的關(guān)鍵在于求解μt.在式(4)關(guān)于湍動能k的方程基礎(chǔ)上，引入湍動能耗散率ε，形成k-ε模型，ε計算式為

(5)

可得

(6)

式中：Cμ為經(jīng)驗系數(shù)，根據(jù)LAUNDER的研究，Cμ取值為0.09[17].

k-ε模型比零方程和一方程模型在計算精度和計算成本上有了很大改進，且具有廣泛的適用性，但是k-ε模型在推導(dǎo)過程中仍進行了一定的假設(shè)和近似，比如它是針對湍流發(fā)展充分的高湍流區(qū)域建立的，在近壁區(qū)等湍流發(fā)展不充分的地方可能存在失真，因此需要采用壁面函數(shù)法等措施進行修正；此外，k-ε模型在強旋流、壁面彎曲流動時計算存在一定誤差.

文中基于k-ε模型，利用局部時均法引入fk和fε對運輸方程進行修改，形成PANS模型.fk，fε定義如下

(7)

式中：k，ε分別為總湍動能和總耗散率；下標(biāo)u代表PANS模型中未分解的物理量.

PANS模型將標(biāo)準(zhǔn)k-ε模型的運輸方程改為

(8)

(9)

并針對式(9)中的以下系數(shù)進行了修正

(10)

式中：Cε1=1.44；Cε2=1.92；σk=1.0；σε=1.3.

PANS模型是可以從直接模擬模型(DNS)過渡到雷諾時均化模型(RANS)的計算模型，這個過程通過調(diào)整模型中的fk和fε進行控制，fk和fε的值取[0,1].根據(jù)大量文獻，fε在高雷諾數(shù)時取值為1更為合適[18].而fk的取值對計算的結(jié)果影響較大，當(dāng)fk=1時，PANS控制方程將復(fù)原到RANS模型；當(dāng)fk=0時，PANS控制方程直接對N-S方程求解，即DNS模型，在fk值由小變大的過程中，PANS模型從DNS平滑過渡到RANS模型.當(dāng)fk取值較小時，雖然可以使得求解結(jié)果更為精確，但對應(yīng)的網(wǎng)格數(shù)量增多，且對計算機能力要求也大大增加，因此在PANS模型計算中fk的取值尤為關(guān)鍵.

目前應(yīng)用較為廣泛的取值方法是將fk值取為[0,1]的一個數(shù)或幾個數(shù)，然后進行計算,從而選出較為適合的值，此方法是在整個計算域取統(tǒng)一的值，需要較大的計算量進行對比且無法得到不同區(qū)域內(nèi)最合適的fk值.

文中利用CFX軟件二次開發(fā)技術(shù)中的CCL語言，對fk值進行了動態(tài)定義，即空間上瞬時地根據(jù)當(dāng)?shù)鼐W(wǎng)格條件和湍流長度尺度定義fk的值，形成修正的PANS(modified partially-averaged Navier-Stokes， MPANS)模型，fk公式為[19]

fk=min[1,3(Δ/l)2/3],

(11)

式中：l為當(dāng)?shù)赝牧鏖L度尺度，l=k1.5/ε；Δ為當(dāng)?shù)鼐W(wǎng)格尺度，Δ=(ΔxΔyΔz)1/3.

由式(11)可知，fk的取值可以根據(jù)當(dāng)?shù)鼐W(wǎng)格和湍流尺度進行動態(tài)修改，當(dāng)湍流尺度較小時，fk=1，此時計算方法改為RANS法；而當(dāng)湍流尺度較大時，fk<1，此時計算方法趨向于直接計算湍流方程[20]，可更為準(zhǔn)確地模擬出該區(qū)域的流動特性.

修正的PANS模型解決了標(biāo)準(zhǔn)PANS模型無法精確取到合適的fk值，以及無法在不同區(qū)域取不同fk值的問題；與普通雷諾時均化模型相比，修正PANS模型在高湍流流動區(qū)域內(nèi)數(shù)值模擬的精度更高；與直接模擬相比，在湍流尺度小的區(qū)域進行了時均化處理，減少了計算成本，提高了網(wǎng)格尺寸的適應(yīng)性.因此，在湍流計算中選取修正的PANS模型更為合適.

2 數(shù)值模擬

2.1 計算域及網(wǎng)格劃分

文中研究對象為立式軸流泵裝置，根據(jù)工程要求對全流道進行三維建模，模型由進水流道、葉輪、后置導(dǎo)葉體、出水流道構(gòu)成，模型額定轉(zhuǎn)速為1 367 r/min，設(shè)計流量為0.375 m3/s，葉輪直徑為300 mm，輪轂直徑為120 mm，葉輪葉片數(shù)為4，葉頂間隙為0.3 mm.利用ICEM進行六面體結(jié)構(gòu)化網(wǎng)格劃分，為了準(zhǔn)確反映邊界層內(nèi)的流動情況，精確模擬靠近葉片表面的空化形態(tài)，對網(wǎng)格的Y+值進行了控制，通過調(diào)整固體壁面第1層網(wǎng)格的高度來改變Y+值[21]，將其控制在30以內(nèi)，流體域整體網(wǎng)格數(shù)量約為300萬.

2.2 求解方法設(shè)置

利用CFX軟件進行軸流泵全流道空化定常計算，將修正的PANS模型導(dǎo)入湍流控制方程中，計算離散方法采用有限體積法[22]，湍流模型對流項采用高階精度、中心差分格式，收斂控制設(shè)置為最少迭代步數(shù)不少于5 000步，殘差收斂精度設(shè)置為10-5.空化模型采用Zwar-Gerber-Belamri模型[23]，氣液兩相之間的傳輸采用基于Rayleigh-Plesset方程[24-25]的計算方法.進口邊界條件采用總壓進口，以控制在不同汽蝕余量下的空化計算，進口處氣相體積分數(shù)為0，液相體積分數(shù)為1，出口設(shè)置為質(zhì)量流量出口.流體域介質(zhì)設(shè)定為“Water at 25 ℃”，“Water Vapour at 25 ℃”，水在25 ℃時對應(yīng)的飽和蒸汽壓力為3 175 Pa.

2.3 計算結(jié)果分析

根據(jù)基于修正PANS模型的定常計算結(jié)果，對泵的空化性能進行綜合分析.為了準(zhǔn)確對比不同工況下葉輪內(nèi)同一位置處的性能，在葉輪內(nèi)進行分層設(shè)置，在葉輪內(nèi)定義某個徑向截面，引入量綱為一的變量Span,

Span=r/R，

(12)

式中：r為該截面距離輪轂的半徑長；R為葉輪室距離輪轂的半徑長.Span=0時，該截面為輪轂面；Span=1時，該截面為葉輪室壁面.

2.3.1 不同汽蝕余量下的空化性能

選取不同汽蝕余量即NPSH為6.81，5.18和4.52 m時的計算結(jié)果，進行空泡體積分數(shù)δ、渦量Ω分布、壓力p分布的對比，在葉輪內(nèi)取截面Span=0.95，Span=0.50并進行平鋪展開，對比結(jié)果見圖1，2.

圖1 NPSH=6.81，5.18，4.52 m時Span=0.95處空化特性對比

圖2 NPSH=6.81, 5.18, 4.52 m時Span=0.50處空化特性對比

根據(jù)圖1，2的對比可知，從葉輪輪緣至葉輪中部的徑向分布上不同汽蝕余量下表現(xiàn)的規(guī)律一致，隨著汽蝕余量的減小，各截面上空泡體積分數(shù)變大，速度渦量變大，壓力逐漸減小，而這些變化在葉片的吸力面表現(xiàn)得更為明顯，這與空泡的主要發(fā)生位置相對應(yīng).對比截面Span=0.95和Span=0.50處的結(jié)果可知，葉片上的空泡在靠近輪緣處相較于葉片中部更多，葉片輪緣處空化相對較嚴重，且輪緣處的空泡主要分布在葉片的前端，而葉輪中部的空泡則在葉片的尾部較多；葉輪中部速度渦量值相比于輪緣處較小，且整體分布規(guī)則，整體大渦量區(qū)域減小，葉片表面壓力由中部向輪緣處逐漸增大，葉片中部由于空化較弱，壓力相差值減小，壓力分布較輪緣處均勻.通過以上的分析可知，葉片表面壓力和速度的變化與空化發(fā)生的位置和區(qū)域大小息息相關(guān)，不同汽蝕余量下各變量也表現(xiàn)出不同的規(guī)律.

2.3.2 不同流量下的空化性能

在不同流量工況下進行了數(shù)值計算，分別選取流量為0.87Q，0.93Q，1.00Q，1.07Q，1.13Q，計算得出對應(yīng)的臨界汽蝕余量分別為4.97，5.16，5.37，5.72 和6.40 m.隨著流量的增大，揚程降低，臨界汽蝕余量增大，說明泵是否發(fā)生空化不僅與進口壓力有關(guān)，也與流量大小有關(guān)，在某一固定的汽蝕余量下，如5.37 m時，泵在設(shè)計流量下運行時處于臨界空化狀態(tài)，而在大流量工況下已發(fā)生嚴重空化，小流量工況下則尚未發(fā)生空化.

分析流量對空泡的影響，選取相同汽蝕余量，對比分析流量對軸流泵空化性能的影響.當(dāng)汽蝕余量NPSH=4.87 m時，對不同流量工況下泵葉輪內(nèi)的空化狀態(tài)進行分析，由圖3各流量下的空泡分布可知，空泡體積分數(shù)為0.15.

圖3 不同流量下的空泡分布圖

由圖3可知，在相同汽蝕余量下，不同流量工況表現(xiàn)出不同的空泡分布特點.隨著流量的增大，空泡分布的整體量在變大，空化逐漸嚴重，葉片吸力面的空泡分布由進水邊逐漸向出水邊轉(zhuǎn)移，而葉片壓力面進水邊區(qū)域的空泡出現(xiàn)了從無到有的變化，當(dāng)流量達到1.13Q時，葉片吸力面空泡與壓力面空泡的產(chǎn)生和潰滅連為一體，使空泡不只是附著于葉片，而是充斥在葉間流道內(nèi)，此時空化十分嚴重.

分析相同汽蝕余量(4.87 m)下不同流量工況時的空化對泵內(nèi)流態(tài)的影響，取葉輪和導(dǎo)葉中間截面(Span=0.50)并對其展開，各流量對應(yīng)的流線分布如圖4所示.

圖4 NPSH=4.87 m時不同流量下的流線分布

由圖4可知，在不同流量工況下，從流線的整體分布來看，葉輪內(nèi)部的流線變化不大，導(dǎo)葉內(nèi)流態(tài)變化較大.在小流量工況0.87Q下，導(dǎo)葉吸力面出現(xiàn)較大區(qū)域的旋渦，隨著流量的增加，旋渦逐漸減小直至消失；在大流量工況1.13Q下，導(dǎo)葉吸力面流線較為平穩(wěn)，但是在導(dǎo)葉的進口處出現(xiàn)扭曲.隨著空化逐漸嚴重，導(dǎo)葉體內(nèi)的流態(tài)逐漸平穩(wěn)，說明空化未對導(dǎo)葉體內(nèi)流態(tài)產(chǎn)生大的影響，流量的增加使得導(dǎo)葉流道內(nèi)的水流具有一定的夾持作用，能夠減小水流的旋渦.

為了分析不同流量下空化對葉輪內(nèi)流態(tài)穩(wěn)定性的影響，得到圖5所示葉輪內(nèi)不同截面上平均湍動能值曲線.由圖可知，在相同汽蝕余量(4.87 m)下，不同流量工況時葉輪內(nèi)部的湍動能分布不同，隨著流量的增大，湍動能值變大，湍流耗散變嚴重，這與空化逐漸嚴重有關(guān)，空化的發(fā)生導(dǎo)致葉輪內(nèi)部流動的不穩(wěn)定性加劇.從湍動能在葉輪的徑向分布特點來看，各流量工況下湍動能值具有較為明顯的特點，其值整體上隨著流量的增加而增大.對于同一流量工況下，湍動能值在Span=0.05即靠近輪轂處最大，沿徑向先減小至最小值，而后在Span=0.95即靠近輪緣處增大，說明在輪轂和輪緣處湍流耗散較嚴重，水流穩(wěn)定性降低，與空泡的邊緣處重合.大流量下空泡體積分數(shù)較大，其潰滅對水流湍動能的影響更大，這導(dǎo)致了流量1.13Q下葉輪出口和葉輪內(nèi)部的湍動能值均比其他工況下大.綜上可知，軸流泵在不同流量下發(fā)生空化所表現(xiàn)的性能有所不同，空泡的分布、流態(tài)的變化以及湍動能值的大小均隨著流量的改變呈現(xiàn)出規(guī)律性變化.

圖5 不同流量下葉輪內(nèi)各截面平均湍動能值曲線

3 試驗驗證

在對軸流泵進行空化定常計算并得出軸流泵空化特性以后，對軸流泵模型進行了試驗，并對比數(shù)值計算和試驗結(jié)果中的空化特性曲線和空泡分布，驗證文中所采取的計算方法的可靠性.

本次空化試驗在河海大學(xué)水力機械多功能試驗臺上進行，試驗中為了驗證數(shù)值計算的準(zhǔn)確性，對軸流泵的外特性和空化特性均進行了測試，主要涉及的參數(shù)有以下幾個：揚程、流量、轉(zhuǎn)矩、真空度、壓力等.空化試驗采用能量法，保持設(shè)計流量為0.375 m3/s，通過系統(tǒng)回路內(nèi)抽真空減小試驗系統(tǒng)的空化余量值，從而減小有效汽蝕余量.在抽取不同真空度值時測試軸流泵對應(yīng)的揚程，并繪制揚程和汽蝕余量的曲線，根據(jù)揚程下降2%確定臨界汽蝕余量.繪制試驗汽蝕余量-揚程曲線，并與數(shù)值模擬計算結(jié)果進行對比，如圖6所示.

圖6 試驗與數(shù)模空化曲線對比

由圖6可知，隨著汽蝕余量的減小，泵的揚程逐漸減小，空化越來越嚴重，試驗與數(shù)值模擬結(jié)果曲線趨勢相同，但試驗測得的揚程較數(shù)值計算結(jié)果略大，兩者誤差小于2.9%，這是由試驗條件及試驗系統(tǒng)引起的誤差，且在合理范圍內(nèi).圖中點a為試驗所測得的臨界汽蝕余量，其值為5.68 m；點A為數(shù)值模擬計算所得臨界汽蝕余量，其值為5.37 m，試驗測得的臨界汽蝕余量較數(shù)值計算結(jié)果偏大，空化性能較數(shù)值結(jié)果略差.這是由于實際液體中包含某些不可溶氣體，形成了空化核，空化核的存在將會降低液體的抗拉強度，試驗時若對水體雜質(zhì)的排除不夠徹底，將會增加水體中空化核的含量，而數(shù)值模擬時默認的空化核數(shù)量較少，這將導(dǎo)致試驗空化性能比數(shù)值計算結(jié)果差.另外，試驗時葉輪體各壁面具有一定的粗糙度，當(dāng)泵內(nèi)水流流經(jīng)某些粗糙度較大的地方時，將會引起速度和壓力的突變，增加附近區(qū)域的湍流強度，導(dǎo)致在粗糙區(qū)域內(nèi)產(chǎn)生局部空化，也在一定程度上降低了泵的空化性能.

試驗中除對效率、揚程等數(shù)據(jù)進行采集外，還在葉輪室壁面上安裝了透視窗口，利用閃頻儀對葉輪內(nèi)的空泡進行拍攝，拍攝結(jié)果如圖7所示.

選取試驗和數(shù)值計算中空化開始階段和空化較為嚴重的2個工況，即圖6中的點B，b和點C，c，分別進行空泡分布對比，數(shù)值計算結(jié)果按照空泡體積分數(shù)0.15以上顯示，如圖7所示.

圖7 空泡分布對比圖

由圖7可看出，數(shù)值模擬下的空泡與試驗所拍攝的結(jié)果基本一致，隨著汽蝕余量的減小，試驗和數(shù)值模擬空化逐漸嚴重，空泡逐漸變多，說明數(shù)值計算能夠較為準(zhǔn)確地模擬泵內(nèi)的空化性能以及空泡的分布情況，修正PANS模型在軸流泵的空化模擬中可信度較高.

4 結(jié) 論

1) 基于標(biāo)準(zhǔn)PANS模型，利用CFX軟件二次開發(fā)技術(shù)中的CCL語言，對fk值進行了動態(tài)定義，使其可以瞬時地根據(jù)當(dāng)?shù)鼐W(wǎng)格條件和湍流長度尺度修改其值，形成了修正PANS模型，解決了標(biāo)準(zhǔn)PANS模型在計算中無法精確取到合適的fk值以及在不同區(qū)域取不同fk值的問題.

2) 建立軸流泵全流道模型進行數(shù)值計算，引入量綱為一的變量Span對葉輪內(nèi)部進行分層，對比不同工況下葉輪內(nèi)部的空化特性.分析得出隨著汽蝕余量的減小，葉輪內(nèi)部的空泡體積分數(shù)和速度渦量逐漸變大，葉片上壓力和速度的變化與空化發(fā)生的位置相對應(yīng).

3) 對不同流量工況進行數(shù)值計算，通過監(jiān)測揚程繪制NPSH-H曲線，對比分析了在不同流量下泵發(fā)生空化的臨界點，以及流量對泵葉輪內(nèi)部的空泡分布、流線分布和湍動能分布的影響.

4) 經(jīng)過試驗驗證可知，隨著汽蝕余量的減小，揚程逐漸減小，空化逐漸嚴重，試驗與數(shù)值計算結(jié)果趨勢相同.利用閃頻儀對葉輪內(nèi)的空泡進行拍攝，所得空泡照片與數(shù)值計算結(jié)果基本一致，說明修正PANS模型能夠較為準(zhǔn)確地模擬泵內(nèi)的空化性能以及空泡的分布情況，在軸流泵的空化模擬中可信度較高.