劉宇寧，王秀勇，劉志遠，白小榜

（1.蘭州理工大學 能源與動力工程學院，蘭州 730050；2.重慶水泵廠有限責任公司國家級企業(yè)技術(shù)中心，重慶 400030）

0 引言

多級離心泵在工農(nóng)業(yè)生產(chǎn)中的應用極為廣泛，由于多級泵水力模型的開發(fā)制造成本較高，因而在設計過程中通常借助于CFD技術(shù)對其水力性能進行預測，并對泵內(nèi)的流場結(jié)構(gòu)進行分析研究，為水力性能的進一步優(yōu)化提供有效的參考信息［1］。雖然CFD技術(shù)在流體機械行業(yè)中已被廣泛應用，但數(shù)值計算精度問題依然是當前研究的重點內(nèi)容。由于目前關于多級離心泵的數(shù)值模擬研究多集中于泵內(nèi)流動特征的分析，而關于提高水力性能預測精度的研究還很少，因而在既有理論基礎上，有必要對影響多級離心泵數(shù)值計算精度的因素進行研究，從而為提高水力性能的預測精度提供可靠而又可行的數(shù)值模擬方法。

研究結(jié)果表明，影響數(shù)值計算精度的主要因素有流動狀態(tài)是否定常、湍流模型的選取以及近壁面網(wǎng)格尺度等［2-5］。在工程應用中，關于流動控制方程的求解通常采用雷諾時均化（RANS）方法，該方法具有較強的計算經(jīng)濟性，對計算機資源要求不高，在離心泵數(shù)值模擬中的應用最為廣泛［6-7］，但由于時均化過程會導致部分湍流脈動信息丟失，流場保真性較差；大渦模擬（LES）方法在反映湍流脈動特征方面具有RANS無可比擬的優(yōu)勢，常用于捕捉泵內(nèi)部的渦核分布［8］，計算精度較高，但對計算資源的要求也很高；近年來，結(jié)合了LES與RANS各自優(yōu)點的分離渦模擬（DES）方法在流體機械的數(shù)值模擬研究中取得了良好的效果［9-10］，但計算過程中可能出現(xiàn)模化應力不足問題，而延遲分離渦模擬（DDES）方法可以對該問題進行修正，并且在離心泵數(shù)值模擬中的應用效果比較明顯［11］，有望成為今后離心泵數(shù)值模擬所選用湍流模型的主流。

由于目前探討離心泵數(shù)值計算精度影響因素的研究對象主要集中于單級單吸結(jié)構(gòu)，為了系統(tǒng)地研究各因素對多級離心泵數(shù)值計算精度的影響，本文以某型號的多級離心泵為研究對象，采用ANSYS FLUENT計算軟件，首先對基于RNGk-ε湍流模型的定常和非定常計算的結(jié)果進行對比；然后比較分析RNGk-ε和DDES兩種湍流模型在多級離心泵數(shù)值模擬中的應用效果；最后探討葉輪近壁面網(wǎng)格尺度的不同對數(shù)值計算結(jié)果的影響。通過上述研究，為多級離心泵的數(shù)值模擬方法提供有價值的參考。

1 數(shù)值計算方法

1.1 計算模型

本文以并聯(lián)結(jié)構(gòu)的節(jié)段式多級泵為研究對象，結(jié)合試驗泵的實際情況，級數(shù)取為兩級，其設計參數(shù)為：流量Q=2 400 m3/h，揚程H=190 m，轉(zhuǎn)速n=1 480r/min，首級和次級葉輪的葉片數(shù)均為6枚，導葉的葉片數(shù)為8枚。由于并聯(lián)式結(jié)構(gòu)關于末級壓出室中心面兩側(cè)對稱，為減少數(shù)值計算的工作量，計算域在此取其中一側(cè)，包括吸入室、首級葉輪、次級葉輪、導葉、壓出室、葉輪前后腔以及口環(huán)間隙，整體幾何模型如圖1所示。

應用ICEM CFD進行網(wǎng)格劃分。由于泵體結(jié)構(gòu)復雜，若進行六面體結(jié)構(gòu)化網(wǎng)格劃分，則工作量巨大，耗時很長，技術(shù)要求極高，難以達到在一般工程應用中進行普遍推廣的目的，因而對吸入室、導葉和壓出室采用容易實施的四面體非結(jié)構(gòu)網(wǎng)格劃分方案。葉輪是影響泵性能好壞的核心部件，葉輪內(nèi)流場結(jié)構(gòu)的模擬精度是決定泵整體性能計算精度的關鍵，為了研究葉輪近壁面處的網(wǎng)格尺度對數(shù)值計算精度的影響，對葉輪采用兩種網(wǎng)格劃分方案，一種是進行四面體非結(jié)構(gòu)化網(wǎng)格劃分，另一種是進行六面體結(jié)構(gòu)化網(wǎng)格劃分并加邊界層。葉輪的網(wǎng)格劃分情況如圖2所示。

圖1 幾何模型

圖2 葉輪網(wǎng)格劃分

根據(jù)網(wǎng)格無關性驗證結(jié)果，當網(wǎng)格總數(shù)增加到1 100×104后，揚程和效率的計算值趨于穩(wěn)定，此時吸入室、導葉和壓出室的網(wǎng)格尺度設定為12 mm；葉輪采用四面體網(wǎng)格劃分時設定葉片表面網(wǎng)格尺度為6 mm，全局網(wǎng)格尺度為8 mm；葉輪采用六面體網(wǎng)格劃分時，設定近壁面第一層網(wǎng)格的高度為0.05 mm，網(wǎng)格尺度增長比率為1.2。

1.2 湍流模型

RNGk-ε湍流模型［12］采用重正化群（RNG）的統(tǒng)計技術(shù)，它會根據(jù)不同強度的旋流來適當?shù)馗淖兺牧黟ざ?，進而對旋流流動進行修正，同時它也兼顧了平均流動中的旋轉(zhuǎn)及旋流情況。

DDES湍流模型［13］是在DES的基礎上，將延遲函數(shù)引入到方程當中，混合長度尺度被重新構(gòu)建，從而解決了DES模型中由網(wǎng)格誘導產(chǎn)生的分離流動（Grid Induced Separation，GIS）問題，進而減弱了DES對網(wǎng)格的依賴性。SST-DDES對長度尺度的修正公式為：

式中l(wèi)DDES——DDES混合長度尺度，m；

lRANS——RANS混合長度尺度，m；

fd——DDES模型的延遲函數(shù)；

CDES——DES函數(shù)校準常數(shù)，CDES=0.61；

Δmax——本地網(wǎng)格單元最大尺度，m；

Cd1，Cd2——經(jīng)驗常數(shù)，Cd1=20，Cd2=3；

rd——當?shù)赝牧鞒叨扰c到壁面距離的比值；

v——運動黏度，m2/s；

vt——旋渦運動黏度，m2/s；

k——卡門通用常數(shù)，k=0.412；

d——網(wǎng)格單元到壁面距離，m；

S——應變率張量，s-1；

Ω——渦量張量，s-1。

1.3 計算方法

采用有限體積法對控制方程進行離散。壓力和速度的耦合采用SIMPLEC算法。定常計算時，采用動參考系模型；非定常計算時，采用滑移網(wǎng)格模型，并使用前面的定常計算結(jié)果作為初始流場。壁面采用無滑移邊界條件；入口給定速度邊界條件；出口給定自由出流邊界條件；靜止域與旋轉(zhuǎn)域之間采用Interface邊界條件進行數(shù)據(jù)傳遞。將葉輪每轉(zhuǎn)過3°圓周角所需要的時間定義為時間步長，其值為3.378 38×10-4s。定常計算時，當泵出口總壓趨于穩(wěn)定時，可視為計算結(jié)果收斂；非定常計算時，當泵出口總壓呈現(xiàn)周期性波動時，取物理量在最后一個旋轉(zhuǎn)周期內(nèi)的時均值進行性能計算。

2 計算結(jié)果分析

2.1 定常與非定常計算結(jié)果對比

圖3示出全計算域采用四面體網(wǎng)格時，在全流量工況范圍內(nèi)，RNGk-ε湍流模型分別采用定常與非定常2種計算方法預測得到的揚程、軸功率和效率的計算結(jié)果。

圖3 RNG k-ε湍流模型定常與非定常計算結(jié)果比較

由圖3（a）所示的揚程計算結(jié)果來看，定常與非定常計算在全流量工況范圍內(nèi)的結(jié)果有明顯的不同，其中非定常計算的揚程在各個工況點均高于試驗值，但揚程計算曲線隨流量增加時的變化趨勢與試驗曲線吻合良好；而定常計算的揚程在（0～0.3）Qd以及（1.1～1.2）Qd工況范圍內(nèi)高于試驗值，在（0.3～1.1）Qd工況范圍內(nèi)低于試驗值，揚程的計算曲線與試驗曲線存在兩個交點，導致?lián)P程計算曲線隨流量增加時的變化趨勢與試驗曲線不完全一致，在（0～0.8）Qd工況范圍內(nèi)下降速度較快，在（0.8～1.2）Qd工況范圍內(nèi)下降速度又相對緩慢。從揚程相對計算誤差大小的角度來看，在（0.2～0.5）Qd以及（0.9～1.2）Qd工況范圍內(nèi)，定常計算的揚程計算誤差小于非定常計算，原因是這兩個工況范圍處于定常計算曲線與試驗曲線的交點位置附近，當偏離這兩個工況范圍時，定常計算的相對計算誤差大于非定常計算；定常計算的揚程相對計算誤差最大值在關死工況點，為12.23%，其次是0.8Qd工況點，為9.10%，其他工況點的計算誤差位于2.54%～6.61%之間；非定常計算的揚程相對計算誤差最大值在1.2Qd工況點，為10.47%，其他工況點的計算誤差相差不大，均位于4.13%～7.32%之間。

由圖3（b）所示的軸功率計算結(jié)果來看，定常與非定常計算的結(jié)果差異明顯，軸功率的非定常計算值在全流量工況點都高于試驗值，但其計算曲線隨流量增加時的變化趨勢與試驗曲線吻合良好，能夠準確反映出軸功率與流量之間的特性關系；而軸功率的定常計算值在0.6Qd工況點之后小于試驗值，之前大于試驗值，并且在關死工況點附近遠大于試驗值，導致定常計算曲線的變化趨勢僅在（0.6～1.2）Qd工況范圍內(nèi)與試驗曲線基本一致，但在（0～0.6）Qd工況范圍內(nèi)隨流量的增加呈現(xiàn)下降的趨勢，與試驗曲線的變化趨勢相反，不能正確反映軸功率在小流量工況區(qū)間內(nèi)的變化特征。從軸功率相對計算誤差大小的角度來看，定常與非定常計算均在關死工況點存在最大計算誤差，其中定常計算的相對計算誤差為53.36%，非定常計算為22.80%；定常計算的最小誤差在0.6Qd工況點，僅有0.09%，該工況點剛好是軸功率的定常計算曲線與試驗曲線的交點；非定常計算的誤差隨著流量的增加逐漸減小，在1.2Qd工況點達到最小值，為5.06%。

由圖3（c）所示的效率計算結(jié)果來看，定常與非定常計算的效率計算值均在設計工況點之前小于試驗值，之后大于試驗值；在（0～0.8）Qd工況范圍內(nèi)，效率的非定常計算曲線與試驗曲線的吻合度明顯高于定常計算；在（0.8～1.2）Qd工況范圍內(nèi)，效率的定常與非定常計算結(jié)果比較接近，隨著流量的增加效率的計算值也在增大，而試驗值則變化平緩，效率的計算曲線與試驗曲線在設計工況附近有明顯的區(qū)別。從效率相對計算誤差大小的角度來看，定常與非定常計算的最大誤差均在0.2Qd工況點，其中定常計算的相對誤差大小為23.72%，非定常計算為7.36%；定常計算誤差的次高點在0.3Qd工況點，為14.89%，其他工況點相對計算誤差的大小位于1.95%～7.25%之間；非定常計算誤差的次高點在1.2Qd工況點，為6.22%，其他工況點相對計算誤差的大小位于1.98%～3.97%之間。

總體來看，在全流量工況范圍內(nèi)，多級泵各性能參數(shù)的非定常計算曲線與試驗曲線的吻合度都遠高于定常計算，尤其是效率的非定常計算結(jié)果，無論是在計算精度方面還是在計算曲線隨流量的變化趨勢方面，都具有定常計算不可比擬的明顯優(yōu)勢。因此，對多級泵進行性能預測時，采用非定常計算的方法可以獲得更精確的計算結(jié)果。

2.2 RNG k-ε與DDES湍流模型計算結(jié)果對比

圖4示出全計算域采用四面體網(wǎng)格時，在全流量工況范圍內(nèi)進行非定常計算的前提下，RNGk-ε與DDES 2種湍流模型關于揚程、軸功率和效率的計算結(jié)果。

圖4 RNG k-ε與DDES湍流模型非定常計算結(jié)果比較

由圖4（a）所示的揚程計算結(jié)果來看，在全流量工況范圍內(nèi)，DDES湍流模型的揚程計算值明顯低于RNGk-ε湍流模型，其計算曲線與試驗曲線之間具有極高的吻合度，計算精度遠高于RNGk-ε湍流模型。DDES的相對計算誤差最大值在1.2Qd工況點，為2.48%，其次是0.4Qd工況點，為2.09%，在其他工況點的相對計算誤差值均小于1.30%。

由圖4（b）所示的軸功率計算結(jié)果來看，在全流量工況范圍內(nèi)，DDES湍流模型的軸功率計算值明顯低于RNGk-ε湍流模型，并且也都低于試驗值，但二者的軸功率計算曲線隨流量增加時的變化趨勢仍然保持較高的相似度。從軸功率相對計算誤差大小的角度來看，DDES湍流模型在全流量工況范圍內(nèi)的計算誤差都明顯低于RNGk-ε湍流模型，其相對計算誤差的最小值在關死工況點，為0.05%，其次是1.2Qd工況點，為1.95%，在（0.2～1.0）Qd工況范圍內(nèi)的相對計算誤差值比較穩(wěn)定，基本維持在4.22%左右。

由圖4（c）所示的效率計算結(jié)果來看，在全流量工況范圍內(nèi)，DDES湍流模型的效率計算值明顯高于RNGk-ε湍流模型，并且也都高于試驗值，但DDES湍流模型的效率計算曲線在設計工況點附近更為平坦一些，可以更準確地反映效率隨流量增加時的變化趨勢。從效率相對計算誤差大小的角度來看，在全流量工況范圍內(nèi)，DDES湍流模型的計算誤差值均高于RNGk-ε湍流模型，其計算誤差的平均值比RNGk-ε高出1.7%左右。

總體來看，當計算域全局采用四面體網(wǎng)格時，DDES湍流模型在多級泵的揚程和軸功率的計算精度上比RNGk-ε湍流模型具有明顯的優(yōu)勢，雖然在預測效率計算曲線隨流量增加時的變化趨勢上也具有一定的優(yōu)勢，但在計算精度上與RNGk-ε湍流模型相比卻處于明顯的劣勢狀態(tài)。實際上，DDES湍流模型對近壁面的網(wǎng)格尺度具有較高的要求，而四面體網(wǎng)格劃分的方法遠達不到其要求，因而下面將進一步探討近壁面網(wǎng)格尺度的大小對DDES湍流模型計算結(jié)果的影響。

2.3 葉輪近壁面網(wǎng)格尺度對計算結(jié)果的影響

吸入室、導葉和壓出室均采用四面體網(wǎng)格，葉輪分別采用四面體（無邊界層）和六面體（有邊界層）網(wǎng)格進行數(shù)值計算。根據(jù)數(shù)值計算的統(tǒng)計結(jié)果，葉輪近壁面網(wǎng)格的y+值分布情況如圖5所示，從圖可以看出，網(wǎng)格無邊界層時的y+值集中分布在250～1 000的區(qū)間內(nèi)，有邊界層時的y+值集中分布在10～30的區(qū)間內(nèi)。雖然葉輪的兩種網(wǎng)格劃分方案的y+值均不滿足DDES湍流模型y+＜5的要求，但仍可以用來定性討論分析近壁面網(wǎng)格尺度的變化對計算結(jié)果的影響，更細致的研究將在后期展開。

圖5 不同網(wǎng)格劃分方案近壁面y+值分布情況

圖6示出采用非定常計算時，DDES湍流模型在葉輪有無邊界層網(wǎng)格的情況下關于揚程、軸功率和效率的計算結(jié)果。

圖6 DDES湍流模型有無邊界層計算結(jié)果比較

由圖6（a）所示的揚程計算結(jié)果來看，當在葉輪壁面附近加入邊界層使y+值大幅減小時，揚程的計算曲線與試驗曲線幾乎完全重合，與無邊界層的計算結(jié)果相比，加入邊界層后揚程的計算精度有了明顯提高。葉輪網(wǎng)格加入邊界層后，揚程的最大相對計算誤差在1.2Qd工況點，僅有1.43%，其次是設計工況點，為1.06%，其他工況點的相對計算誤差值均在0.50%左右。

由圖6（b）所示的軸功率計算結(jié)果來看，在葉輪壁面附近加入邊界層時，軸功率的計算值比無邊界層時有所提高，但在0.2Qd工況點之后仍小于試驗值。在考慮了邊界層對計算結(jié)果的影響后，除了在關死點軸功率的相對計算誤差比無邊界層時增加了1.17%外，在其他工況點的相對計算誤差均明顯減小，最低可以減小0.55%，最高減小了2.87%，平均能夠減小1.20%。

由圖6（c）所示的效率計算結(jié)果來看，有邊界層時的效率計算曲線與試驗曲線的吻合度進一步提高，在全流量工況范圍內(nèi)比無邊界層時的計算結(jié)果更精確，除了關死工況點的效率在理論上不存在計算誤差外，其他工況點有邊界層時的相對計算誤差均比無邊界層時減小了2.50%左右。

總體來看，當減小葉輪近壁面處的網(wǎng)格尺度時，DDES湍流模型關于多級泵各性能參數(shù)的計算精度均有明顯提高，并且可以推測，當葉輪尤其是所有過流部件近壁面網(wǎng)格尺度的y+值完全滿足DDES湍流模型的要求時，其計算精度會有進一步的提高。

2.4 各數(shù)值計算方法計算誤差綜合分析

為了綜合比較上述各數(shù)值計算方法的優(yōu)劣，在此對各性能參數(shù)在全流量工況范圍內(nèi)相對計算誤差大小的平均值及相對計算誤差的均方差進行統(tǒng)計，如表1所示。

表1 全流量工況點計算誤差分析 （%）

相對計算誤差大小的平均值可以反映計算值與試驗值之間的偏離程度，而均方差則可以反映計算曲線與試驗曲線在曲線形狀上的相似度；平均值越小，則計算值越接近試驗值；均方差越小，則計算曲線隨流量增加時的變化趨勢越接近試驗曲線。

由表1所示各計算方法相對計算誤差的對比來看，在全流量工況點，定常計算方法的計算誤差最大，計算曲線與試驗曲線在曲線形狀的相似度上也最低，而非定常計算方法在上述兩個問題上都有明顯改善。計算域全局采用四面體網(wǎng)格并且進行非定常計算時，雖然DDES湍流模型關于效率的計算誤差比RNGk-ε湍流模型高出了1.72%，但關于揚程、軸功率的計算誤差卻都比RNGk-ε湍流模型低了4%，并且各性能參數(shù)相對計算誤差的均方差也都比RNGk-ε湍流模型低很多，也就是DDES湍流模型在預測性能曲線隨流量增加時的變化趨勢上比RNGk-ε湍流模型要精確很多，所以綜合來看，在相同的計算方法下，DDES湍流模型比RNGk-ε湍流模型更適合于預測多級泵的性能參數(shù)。當采用DDES湍流模型進行非定常計算時，在葉輪近壁面區(qū)域加入邊界層網(wǎng)格后，各性能參數(shù)的計算誤差進一步減小，計算曲線與試驗曲線的吻合度進一步提高。

綜上所述，對多級泵的性能參數(shù)進行預測時，采用DDES湍流模型和非定常計算方法，能夠獲得較高的計算精度，尤其是對葉輪的近壁面網(wǎng)格進行加密處理時，可以使計算精度進一步提高。在該數(shù)值計算方法中，網(wǎng)格劃分容易實現(xiàn)，網(wǎng)格數(shù)量容易控制，非定常計算所需要的時間也僅為定常計算的4～5倍，但計算精度卻比常規(guī)的采用RNGk-ε湍流模型和定常計算的方法高出很多，適合于一般工程應用。

3 結(jié)論

（1）在全流量工況點，非定常計算方法不僅計算誤差比定常計算方法小，并且性能計算曲線的形狀也更接近于試驗曲線；定常計算時在個別工況點的計算精度極高具有一定的偶然性，與該工況點距離計算曲線和試驗曲線交點位置的遠近有關。

（2）全計算域采用四面體網(wǎng)格和非定常計算的條件下，DDES湍流模型關于效率的計算精度稍低于RNGk-ε湍流模型，但關于揚程和軸功率的計算精度卻遠高于RNGk-ε湍流模型，并且各性能參數(shù)計算曲線的形狀更接近于試驗曲線。

（3）當對近壁面網(wǎng)格進行加密處理時，DDES湍流模型的計算精度可以進一步提高。