劉 毅,王秀勇,2*,董 峰,杜永峰

(1.蘭州理工大學 能源與動力工程學院,甘肅 蘭州 730050;2.甘肅省流體機械及系統(tǒng)重點實驗室,甘肅 蘭州 730050)

0 引 言

近年來,應用計算流體動力學(computational fluid dynamics,CFD)技術對泵進行水力性能預測及優(yōu)化,已經(jīng)成為泵類產(chǎn)品前期研發(fā)的重要方法[1,2]。CFD的計算精度一直是被關注的重點,同時也是計算流體動力學的核心問題。

在計算域模型表達、網(wǎng)格劃分方式、壁面函數(shù)應用、湍流模型選取、邊界條件設置、流動狀態(tài)是否定常等方面,國內外學者已經(jīng)對影響數(shù)值計算精度的因素進行了相關研究[3,4]。李曉俊等人[5]研究了網(wǎng)格質量對離心泵數(shù)值模擬精度的影響,發(fā)現(xiàn)湍流模型的選取與邊界層網(wǎng)格質量密切相關。劉宇寧等人[6]對多級離心泵水力性能數(shù)值模擬精度的影響因素進行了研究,發(fā)現(xiàn)非定常計算的精度高于定常計算。

在研究泵的非定常流動特性時,時間步長是影響其計算精度的一個重要因素,通常要將葉輪每旋轉1°～3°所需要的時間設置為時間步長[7-9]。朱榮生等人[10]對核主泵小流量工況下壓力脈動特性進行研究,將葉輪每旋轉2°作為時間步長進行了非定常計算。王東偉等人[11]在研究離心泵非定?？栈鲌鰰r,將葉輪每旋轉3°作為一個時間步長。但目前并沒有研究對非定常計算過程中如何確定合理時間步長這一問題進行系統(tǒng)分析。

眾所周知,減小時間步長能夠提高非定常計算的精度,而目前的CFD技術還無法完全真實還原泵內的實際流場結構,也就是說,數(shù)值計算始終會存在一定的誤差[12]。

為了研究時間步長與計算精度之間的關系,首先需要明確泵內與時間密切關聯(lián)的主要流動現(xiàn)象。由于各種頻率下的湍流脈動是影響核主泵內流場分布特征的主要因素,由此可以推斷得出,時間步長所對應的響應頻率對湍流脈動頻率譜的覆蓋率,才是影響非定常計算精度的主要因素。因此,有必要對非定常計算精度隨時間步長減小時的變化規(guī)律展開探討研究。

筆者以核主泵為研究對象,對全計算域進行六面體結構化網(wǎng)格劃分,應用FLUENT軟件,選擇基于剪切應力輸運(shear stress transport,SST)k-ω模型的延遲分離渦模擬(delayed detached-eddy simulation,DDES)進行非定常計算[13],分析時間步長變化對核主泵水力性能的計算精度、瞬時性能曲線的分布特征及旋渦場的影響,并在此基礎上,探討時間步長對核主泵內復雜流動特征的解析效應,從而為確定合理的時間步長提供理論參考。

1 計算模型及數(shù)值計算方法

1.1 計算模型

筆者以某型號核主泵原型樣機為研究對象,其設計工況點參數(shù)為:流量Qd=630 m3/h,揚程H=35 m,轉速n=990 r/min,比轉速ns=105,葉輪的葉片數(shù)為6枚,導葉的葉片數(shù)為10枚。

計算域包括進水段、葉輪、導葉、壓水室、前腔、口環(huán)間隙及后腔等過流區(qū)域。

整體幾何模型如圖1所示。

圖1 計算域模型分解圖

1.2 網(wǎng)格劃分

筆者應用ICEM軟件對計算域進行六面體結構化網(wǎng)格劃分。計算域網(wǎng)格劃分情況如圖2所示。

圖2 計算域網(wǎng)格劃分情況

此時,計算域網(wǎng)格總數(shù)量達到1.5×107,完全滿足當模型泵的網(wǎng)格數(shù)為9.42×106時,即可實現(xiàn)網(wǎng)格無關性的要求[14]。

近壁面第一層網(wǎng)格尺度用無量綱壁面距離y+來表示。

各過流部件近壁面y+值的分布情況,如圖3所示。

圖3 壁面y+值分布圖

其中,y+≤5的網(wǎng)格數(shù)量占全計算域近壁面網(wǎng)格總數(shù)的65.77%;510的網(wǎng)格數(shù)量為16.76%。

由于該泵的三維模型在結構上比較復雜,在應用ICEM劃分網(wǎng)格時,若近壁面第一層網(wǎng)格的尺度在現(xiàn)在的基礎上進一步減小,則網(wǎng)格質量變得較差,不能滿足計算對網(wǎng)格質量的要求。

受軟件網(wǎng)格劃分功能的限制,暫且還沒有更好的辦法來進一步降低y+值,但目前已經(jīng)具有的y+值分布也基本能夠滿足DDES模型的要求。

1.3 邊界條件設置

筆者采用動量方程及連續(xù)性方程對整個流場進行求解,湍流模型選擇基于SSTk-ω的DDES,進口采用速度進口,出口采用自由出流邊界條件,各流體域之間的耦合面設置為interface。

筆者采用有限體積法對控制方程進行離散,用SIMPLEC算法實現(xiàn)壓力與速度之間的耦合,壁面采用無滑移壁面邊界條件。

進行定常計算時,動靜域交界面之間采用多重參考坐標系模型;非定常計算時,以定常計算結果作為初始化流場,動靜域交界面采用滑移網(wǎng)格模型。

筆者分別將葉輪每旋轉Δα=1°、1.5°、2°和3°所需要的時間設置為時間步長,以Δα=1°為時間步長完成整個計算所需要的時間,比Δα=1.5°所需的時間增加約50%,近似為Δα=2°所需時間的2倍,Δα=3°所需時間的3倍;

每個時間步長設置迭代20次,收斂精度為10-5;在葉輪旋轉第5圈時,泵出口處總壓的時均值不再發(fā)生明顯變化,采用葉輪旋轉第5圈的數(shù)據(jù)進行后處理分析。

2 水力性能計算結果分析

2.1 計算結果與試驗結果對比

在設計工況點下,4種時間步長的水力性能計算結果與試驗結果對比情況,如表1所示。

表1 各時間步長非定常計算結果與試驗結果的對比

表1中,各性能參數(shù)的計算結果都為時均值,即葉輪旋轉一周的過程中,各性能參數(shù)在每個時間步長下監(jiān)測結果的平均值。

(1)由揚程的計算結果與試驗結果對比可知:當Δα≤1.5°時,揚程的計算值大于試驗值;Δα≥2.0°時,則揚程的計算值小于試驗值;并且,當Δα=1.5°時,揚程的計算值最大,之后隨著時間步長的增大,揚程的計算值持續(xù)減小;Δα≤2°時,揚程在各時間步長下的相對計算誤差大小差異不大,其中,Δα=1.5°的計算誤差相對較大,為0.40%;Δα=3.0°時,揚程的相對計算誤差最大,為1.14%;

(2)由軸功率的計算結果與試驗結果對比可知:Δα=1.0°時的計算值略高于試驗值,從Δα=1.5°開始計算值低于試驗值,此時的相對計算誤差最小,為0.1%;之后,隨著時間步長的增大,軸功率的計算值持續(xù)減小,而相對計算誤差則持續(xù)增大,在Δα=3°時最大,為1.02%;

(3)由效率的計算結果與試驗結果對比可知:Δα=1.5°和Δα=2.0°時的計算值高于試驗值,而Δα=1.0°和Δα=3.0°時的計算值則略低于試驗值(其中,Δα=1.5°時效率的計算值最大,之后隨著時間步長的增大,效率的計算值逐漸減小;Δα=1.5°時效率的相對計算誤差最大,為0.48%,其他3個時間步長下,效率的相對計算誤差位于0.12%～0.18%之間)。

相對計算誤差大小如圖4所示。

圖4 不同時間步長相對計算誤差

由圖4綜合對比來看:隨著時間步長的減小,各性能參數(shù)的相對計算誤差沒有明確的變化規(guī)律,也就是單個性能參數(shù)的計算精度并不一定隨著時間步長的減小而呈現(xiàn)持續(xù)提高的趨勢;但Δα≤2°時,揚程和軸功率的計算誤差明顯低于Δα=3.0°,而效率在所有時間步長下的計算誤差均較小;

因此,綜合各性能參數(shù)的計算誤差可以判斷:在對泵進行非定常計算時,時間步長不大于葉輪每旋轉2°所需要的時間,即可獲得較高的計算精度;這意味著時間步長的大小存在臨界值,當時間步長小于該臨界值時,可以在總體上獲得較高的計算精度。

揚程、軸功率和效率在各時間步長下的相對計算誤差綜合分析情況,如表2所示。

表2 各時間步長相對計算誤差分析

由表2總體來看:時間步長越小,揚程、軸功率和效率三者計算誤差的平均值越小,誤差的離散度也越小;即時間步長越小,各性能參數(shù)的綜合計算精度越高;

根據(jù)各性能參數(shù)計算誤差的平均值和均方差在不同時間步長之間的差異來看:Δα≤2°時,計算誤差的平均值及均方差都明顯小于Δα=3.0°;其中,Δα=2°時,計算誤差的平均值比Δα=3°減小0.43%,均方差減小0.21%;而計算誤差的平均值和均方差在Δα=1°、Δα=1.5°和Δα=2°之間的差異較小,平均值最多只相差0.08%,均方差最多只相差0.09%。

由此可知:當Δα≤2°時,雖然時間步長變小會在總體上提高計算的精度,但提高幅度不明顯;

上述分析結果進一步表明:Δα=2°是核主泵非定常計算精度發(fā)生質變的一個關鍵臨界值;時間步長大于葉輪每旋轉2°所需要的時間時,水力性能的計算精度會明顯降低;

由上述分析可知:當時間步長小于某個臨界值時,隨著時間步長的減小,單個性能參數(shù)的計算精度不一定會持續(xù)提高。雖然所有性能參數(shù)的綜合計算精度會有所提高,但提高幅度較小。由于完成整個非定常計算所需要的時間與時間步長成反比的關系,當時間步長減為原來的一半,完成整個計算所需要的時間比原來增加幾乎一倍。

綜合考慮計算精度和計算所消耗時間這兩個因素,筆者將葉輪每旋轉Δα=2°所需要的時間設置為非定常計算的時間步長,這樣既可以滿足一般工程應用的計算精度要求,又可以大幅節(jié)約計算所消耗的時間。

2.2 水力性能時域分析

筆者在設計工況點進行非定常計算時,采用了4種不同的時間步長,但從計算誤差的平均值和均方差在不同時間步長之間的對比分析來看,Δα=1.5°和Δα=2.0°二者之間具有很高的相似度。故在后續(xù)的研究內容中,筆者僅選取Δα=1°、Δα=2°和Δα=3°的計算結果進行分析。

不同時間步長下水力性能計算值對應的時域圖,如圖5所示。

圖5 不同時間步長下水力性能時域圖T—葉輪旋轉360°所用時間;t—葉輪旋轉60°所用時間

圖5中,從水力性能是否為定常的角度來看:

在葉輪旋轉的過程中,泵的揚程、軸功率和效率并非保持為常數(shù)狀態(tài),而是在一定的區(qū)間范圍內持續(xù)脈動(例如,當Δα=3°時,揚程的脈動幅度為時均值的±1.7%左右,軸功率為±2.8%左右,效率為±1.2%左右;當Δα=2°時,各水力性能的脈動幅度與Δα=3°相接近;而當Δα=1°時,揚程和軸功率的脈動幅度較前二者明顯偏小,但效率略微偏大,其中,揚程在時均值附近的脈動幅度為±1.1%左右,軸功率為±2.1%左右,效率為±1.4%左右);

由于水力性能并非定常狀態(tài),因而對泵的水力性能進行試驗測試時,針對一個工況點進行多次采樣時,取平均值是非常重要的。

由各性能參數(shù)在不同時間步長之間的對比來看:

在各時間步長下,揚程和軸功率的脈動周期均不明顯,而效率的脈動周期相對明顯一些,波峰和波谷分別交替出現(xiàn)30次,為葉輪和導葉的干涉周期;揚程和軸功率的脈動頻率隨時間步長的減小而明顯加快,脈動幅度始終在變,同時曲線分布重心上移,也就是時均值隨時間步長的減小而增大;效率的脈動頻率在各時間步長下總體都保持一致,曲線分布重心也基本相同;只是,當Δα=1°時,效率的瞬時分布曲線出現(xiàn)了更多的微小波動特征。

總體來看:當時間步長減小時,各性能曲線在分布特征上能表現(xiàn)出更多的、小尺度的波動特性,這也正是核主泵內流動特征復雜的細微體現(xiàn);導葉和壓水室內旋渦運動尤其明顯,非定常流動特征加強,從而對核主泵的水力性能產(chǎn)生明顯的影響,導致水力性能周期性變化的特征受到干涉,脈動周期不再明顯,該現(xiàn)象是由泵的內因引起的,并非數(shù)值計算方法的問題。

由于時間步長減小后各水力性能的瞬時曲線在分布特征上會出現(xiàn)更多的小幅度的波動特性,這意味著此時的時間步長對核主泵內更細微的非定常流動特征的分辨能力有所提升,能夠識別出更高頻率下的湍流脈動對流動的影響,因而計算精度會更高一些;這也正是Δα=2°時,各水力性能的綜合計算精度明顯高于Δα=3°的原因。

由于湍流的脈動頻率主要集中在某段區(qū)間內(例如自由射流的湍流脈動頻率譜主要分布在3 500 Hz以內,主頻率在400 Hz左右[15，16]),當時間步長的響應頻率(f=6n/Δα,n為核主泵的轉速,r/min)大于該區(qū)間的上限時,并不會使計算精度得到明顯提高;正如Δα=1°(f=5 940 Hz)時,各水力性能的綜合計算精度略高于Δα=2°(f=2 970 Hz),但提高幅度卻沒有Δα=2°相對于Δα=3°(f=1 980Hz)那么明顯;因為此時Δα=2°所對應時間步長的響應頻率正好位于湍流脈動頻率的上限附近,也就是非定常計算的時間步長存在臨界值,當時間步長小于該臨界值后,對泵內非定常流動特征的響應頻率即可滿足高精度計算的要求。但隨時間步長的繼續(xù),減小計算精度的提高幅度卻不再那么明顯。

3 內流場分析

不同時間步長下動靜葉柵內的旋渦分布狀態(tài),如圖6所示。

圖6 動靜葉柵內旋渦分布圖

不同時間步長下壓水室內的旋渦分布狀態(tài),如圖7所示。

圖7 壓水室內旋渦分布圖

此處,筆者采用Q準則對泵內的渦結構進行提取[17，18],并選取葉輪完成第5個旋轉周期末的結果,對流場進行定性分析。

由圖6可知:葉輪內的旋渦運動在不同時間步長下的分布情況非常相似,流場結構相對簡單,流動的核心區(qū)域幾乎沒有旋渦運動,只在葉片背面及葉輪外緣附近存在兩條細長的帶狀渦(其中,緊靠葉片背面的帶狀渦較窄,而貼近葉輪流道出口外緣的帶狀渦相對較寬);受導葉葉片的干涉作用,帶狀渦在導葉的葉片入口處開始破裂,并流入導葉內部;由于葉輪和導葉的動靜干涉具有徑向對稱的效果,葉輪內的旋渦分布狀態(tài)也呈現(xiàn)出良好的徑向對稱性。

導葉內的流場結構比較復雜,整個流道空間幾乎都被大量形狀各異、分布不均的小尺度旋渦占據(jù);旋渦在導葉流道入口處呈現(xiàn)為帶狀結構,沿著流動方向,帶狀渦逐漸發(fā)展,破碎后成為各種形狀的小尺度旋渦;受葉輪干涉作用時序效應的影響,導葉相鄰兩流道內的旋渦場結構并不完全相同;雖然葉輪和導葉的動靜干涉具有徑向對稱的效果,但導葉內的旋渦分布狀態(tài)并沒有呈現(xiàn)出良好的對稱性(原因在于壓水室內的流場結構是非徑向對稱的,如圖7所示,導葉每個流道的出口邊界條件均不相同,使各個導葉流道內的旋渦場結構出現(xiàn)差異)。

總體來看:導葉內的旋渦場在不同時間步長下的宏觀結構比較相似,但隨著時間步長的減小,捕捉到的渦結構尺度更小、數(shù)量更多,且分布更為離散;與Δα=3°的渦結構相比,Δα=2°的渦結構明顯與Δα=1°的相似度更高,由此也可表明Δα=2°時的時間步長具有足夠精度來解析湍流脈動對導葉內流場的影響。

由圖7壓水室內的旋渦分布可知:由于葉輪沿逆時針方向旋轉,壓水室內液體的循環(huán)方向也為逆時針,內循環(huán)的液體對由導葉左側射入壓水室的液體產(chǎn)生較強的沖擊作用,從而使壓水室左側的旋渦運動強于右側;

Δα=1°與Δα=2°時間步長下的旋渦分布十分相似,且都比Δα=3°時的旋渦數(shù)量多,強度也略大,這進一步表明,較小的時間步長能夠捕捉到更細致的旋渦運動特征,同時也表明,Δα=2°時的時間步長具有足夠的精度來解析壓水室內的湍流運動情況。

湍流脈動是由各種不同尺度和不同旋向的旋渦相互疊加共同作用的結果,大尺度旋渦主要引起低頻脈動,小尺度旋渦主要引起高頻脈動,所以時間步長能夠反映的頻率譜越寬,就能涉及更多小尺度旋渦對流動的影響,數(shù)值計算的精度也就越高;

由不同時間步長下,各過流部件內旋渦場分布情況對比可知:相較于Δα=3°,Δα=2°與Δα=1°所對應的時間步長,均能分辨出更小尺度的旋渦運動,這意味著較小的時間步長能夠解析更高頻率下的湍流脈動情況,因而其計算精度也更高;由于Δα=1°時所捕捉的旋渦場中小尺度旋渦的數(shù)量比Δα=2°略多一些,因而其計算精度也比Δα=2°略高一些;同時可以進一步明確,Δα=2°所對應的時間步長是核主泵非定常計算精度明顯提高的臨界值[19]。

綜上所述,在對核主泵展開非定常計算研究時,時間步長確實存在一個臨界值,該臨界值對應的響應頻率能夠覆蓋核主泵內湍流脈動頻率譜的絕大多數(shù)頻段,從而使計算精度明顯提高;

當時間步長小于該臨界值時,響應頻率雖然會進一步提高,但由于湍流脈動高頻段所占的比例極低,所以其計算精度提升的幅度有限。

4 結束語

筆者通過改變非定常計算過程中的時間步長,對比分析了核主泵水力性能在計算精度、水力性能的脈動特性,及旋渦場分布特征等方面的變化異同點,并探討研究了時間步長對核主泵非定常計算精度的影響機理,確定了時間步長的合理取值范圍。

研究結果如下:

(1)隨時間步長的減小,單項水力性能的計算精度不一定能提高,但各水力性能的綜合計算精度能得以提高;時間步長對應的葉輪旋轉角度Δα=2°是非定常計算精度能夠明顯提高的臨界值;當Δα≤2°時,隨時間步長的減小,非定常計算精度的提高幅度較小;

(2)在葉輪旋轉的過程中,各水力性能的時域圖呈現(xiàn)出脈動的特性;揚程和軸功率的脈動周期不規(guī)律,而效率則呈現(xiàn)出周期性脈動的特征;隨著時間步長的減小,揚程和軸功率的脈動頻率加快,效率的脈動頻率始終保持不變;當Δα≤2°時,各水力性能的瞬時曲線在幾何特征上出現(xiàn)了小尺度波動的形狀,表明此時的時間步長能夠反映流場中更細微的流動特征;

(3)時間步長減小能夠捕捉到更小尺度的旋渦運動,反映出更高頻率的湍流脈動對流動的影響;當Δα=2°時,導葉和壓水室內的渦量場結構與Δα=1°時相似,小尺度渦的數(shù)量均比Δα=3°多,并且分布狀態(tài)更加離散;

(4)時間步長的響應頻率所能覆蓋的湍流脈動頻率譜的范圍,決定了非定常計算的精度,覆蓋范圍越寬,計算精度也越高;以葉輪每旋轉Δα=2°所需要的時間作為時間步長,既可以獲得較高的非定常計算精度,又可以節(jié)約計算所消耗的時間。

在后續(xù)的工作中,筆者將對核主泵內影響湍流脈動頻率的主要因素,以及泵內湍流脈動頻率區(qū)間的確定展開深入研究,為準確定位非定常計算的精度提供理論支撐。