基于PIMPLE算法的低負(fù)荷水泵水輪機(jī)內(nèi)部流態(tài)研究

王 浩 博,周 大 慶,郭 俊 勛

(1.河海大學(xué) 水利水電學(xué)院,江蘇 南京 210098; 2.河海大學(xué) 能源與電氣學(xué)院,江蘇 南京 211100)

0 引 言

近年來,抽水蓄能機(jī)組在低負(fù)荷工況運(yùn)行時(shí)出現(xiàn)的轉(zhuǎn)輪疲勞斷裂[1-3]及運(yùn)行失穩(wěn)[4]等問題日益引起重視。一方面由于活動(dòng)導(dǎo)葉開度有限,導(dǎo)葉出口旋渦結(jié)構(gòu)[5]與無葉區(qū)動(dòng)靜干涉影響相互作用,從而產(chǎn)生異常低頻壓力脈動(dòng)現(xiàn)象[6],并伴隨著壓力脈動(dòng)頻率與葉片固有頻率發(fā)生共振的危險(xiǎn)[7-8];另一方面由于低負(fù)荷工況下,轉(zhuǎn)輪出水邊水體攜帶正向圓周速度分量,從而在尾水管中形成偏心螺旋狀渦帶,并且渦帶在旋轉(zhuǎn)過程中撞擊壁面,引起機(jī)組產(chǎn)生不規(guī)則的偏心振動(dòng),嚴(yán)重威脅著水電站安全運(yùn)行。因此,分析水泵水輪機(jī)在低負(fù)荷工況下內(nèi)部流場變化規(guī)律[9-11],歸納機(jī)組產(chǎn)生異常波動(dòng)的影響因素,對水泵水輪機(jī)穩(wěn)定運(yùn)行具有重要意義。

目前國內(nèi)外學(xué)者一般采用數(shù)值模擬和試驗(yàn)研究等方法對偏離額定工況運(yùn)行的機(jī)組穩(wěn)定性進(jìn)行探究。肖瓊[12]等通過對模型水泵水輪機(jī)泵工況進(jìn)行非定常模擬,總結(jié)出低負(fù)荷工況下無葉區(qū)內(nèi)部的渦團(tuán)分裂聚合對壓力脈動(dòng)變化具有較大影響;趙亞萍[13]等對多能互補(bǔ)系統(tǒng)中水輪機(jī)的轉(zhuǎn)輪進(jìn)行優(yōu)化設(shè)計(jì),得出葉片幾何參數(shù)優(yōu)化可以減弱低負(fù)荷運(yùn)行時(shí)尾水管中產(chǎn)生的壓力脈動(dòng)幅值;Skripkin[14]等基于模型試驗(yàn)結(jié)果,揭示機(jī)組在不同負(fù)荷區(qū)間運(yùn)行時(shí)尾水管內(nèi)部渦帶振動(dòng)分布;林亞濤[15]等對混流式水輪機(jī)進(jìn)行數(shù)值模擬,發(fā)現(xiàn)了在低負(fù)荷工況運(yùn)行時(shí),轉(zhuǎn)輪附近會(huì)出現(xiàn)更為明顯的壓力幅值波動(dòng);Li[16]等通過對原型可逆式機(jī)組進(jìn)行試驗(yàn)研究,劃分出不同負(fù)荷工況下壓力脈動(dòng)幅值分布區(qū)域,并針對其傳播機(jī)制進(jìn)行定量分析。盡管已有上述研究,但目前對多種低負(fù)荷工況水泵水輪機(jī)整體的內(nèi)部流態(tài)分析仍相對較少。

本文采用開源軟件OpenFOAM中的PIMPLE算法對不同低負(fù)荷工況下水泵水輪機(jī)運(yùn)行過程進(jìn)行數(shù)值模擬,揭示了機(jī)組在低負(fù)荷工況下內(nèi)流場特性和無葉區(qū)各頻率壓力脈動(dòng)的變化規(guī)律,探究了尾水管內(nèi)部的渦帶運(yùn)動(dòng)對進(jìn)口流態(tài)分布的影響程度。本文相關(guān)研究結(jié)果可為解決水泵水輪機(jī)運(yùn)行穩(wěn)定性問題提供參考。

1 計(jì)算域分布和網(wǎng)格劃分

本文以白山抽水蓄能電站原型機(jī)組為研究對象,通過UG12.0軟件建立水泵水輪機(jī)全流道模型,相關(guān)計(jì)算域如圖1所示。

圖1 水泵水輪機(jī)計(jì)算域Fig.1 Calculation domain of pump turbine

機(jī)組額定出力為139 MW,額定水頭為105.8 m,相關(guān)機(jī)組尺寸數(shù)據(jù)和運(yùn)行參數(shù)如表1所列。

表1 水泵水輪機(jī)參數(shù)Tab.1 Pump turbine parameters

研究對象采用ICEM軟件對轉(zhuǎn)輪、活動(dòng)導(dǎo)葉及進(jìn)出口管道延伸段進(jìn)行六面體結(jié)構(gòu)化網(wǎng)格劃分,蝸殼和固定導(dǎo)葉采用四面體網(wǎng)格進(jìn)行劃分,并基于OpenFOAM中內(nèi)置前處理工具snappyHexMesh對尾水管進(jìn)行網(wǎng)格劃分和局部加密。采用6種不同網(wǎng)格數(shù)量分配方案對計(jì)算域進(jìn)行網(wǎng)格無關(guān)性驗(yàn)證,網(wǎng)格總數(shù)量與機(jī)組效率變化趨勢如圖2所示。

圖2 網(wǎng)格無關(guān)性驗(yàn)證Fig.2 Mesh independence verification

當(dāng)網(wǎng)格數(shù)量達(dá)到600萬后,效率變化趨近于平緩,考慮到計(jì)算時(shí)間消耗和湍流場近壁面網(wǎng)格精度,最終網(wǎng)格總量確定為682.3萬。圖3為各計(jì)算區(qū)域網(wǎng)格數(shù)量分配結(jié)果:進(jìn)出口管道延伸段為25.8萬、蝸殼區(qū)域?yàn)?16.8萬、導(dǎo)葉區(qū)域?yàn)?50.9萬、轉(zhuǎn)輪區(qū)域?yàn)?82.6萬、尾水管區(qū)域?yàn)?06.2萬。其中在轉(zhuǎn)輪及活動(dòng)導(dǎo)葉等高雷諾數(shù)區(qū)域控制網(wǎng)格節(jié)點(diǎn)y+值不高于15,符合工程中湍流模型應(yīng)用要求。

2 數(shù)值求解方法

2.1 控制方程

本文采用OpenFOAM軟件對水泵水輪機(jī)多種低負(fù)荷工況進(jìn)行數(shù)值模擬,在不考慮流體壓縮性及能量方程求解情況下,選擇單相流不可壓縮求解器。因此Navier-Stokes控制方程主要由動(dòng)量守恒方程及連續(xù)性方程所組成,形式如下:

(1)

?·u=0

(2)

式中:u為流速,m/s;p為壓強(qiáng),Pa;ρ為流體密度,kg/m3;t為時(shí)間,s;?為哈密頓算子;?2為拉普拉斯算符;ν為運(yùn)動(dòng)黏度系數(shù),m2/s。

2.2 初始條件和湍流模型

本文中水泵水輪機(jī)計(jì)算域進(jìn)出口均采用壓力邊界條件,維持機(jī)組凈水頭處于105.8 m,旋轉(zhuǎn)部分設(shè)置為moving Wall Velocity,其余壁面設(shè)置為無滑移邊界條件,通過AMI(Arbitrary Mesh Interface)任意網(wǎng)格交界面插值方式實(shí)現(xiàn)接觸面之間的數(shù)據(jù)傳遞。為保證數(shù)值求解精度,采用二階迎風(fēng)格式離散對流項(xiàng)相關(guān)變量,并基于中心線性差分格式對梯度項(xiàng)和擴(kuò)散項(xiàng)進(jìn)行求解。

穩(wěn)態(tài)計(jì)算中選擇SST k-ω[17]湍流模型,瞬態(tài)計(jì)算中選擇SAS(scale-adaptive simulation)-SST[18]湍流模型。其中尺度自適應(yīng)模擬SAS方法可根據(jù)流動(dòng)不同尺度區(qū)域進(jìn)行自適應(yīng)選擇,在處理近壁面和小范圍分離流動(dòng)上采用RANS方法,而在流動(dòng)處于大尺度分離狀態(tài)時(shí),會(huì)切換為大渦模擬LES(Large eddy simulation)方法,并在一定空間尺度上進(jìn)行濾波處理,因此在處理較大范圍的分離流動(dòng)方面具有一定的優(yōu)勢。

2.3 求解算法

水泵水輪機(jī)穩(wěn)態(tài)運(yùn)行時(shí),采用OpenFOAM中的simpleFoam求解器,結(jié)合MRF(多重參考系法)實(shí)現(xiàn)穩(wěn)態(tài)計(jì)算中轉(zhuǎn)輪相對轉(zhuǎn)動(dòng),在瞬態(tài)計(jì)算中采用pimpleFoam求解器。其求解器原理主要基于PIMPLE算法,將SIMPLE[19](Semi-Implicit Method for Pressure-Linked Equations)半隱式求解方法及PISO[20](Pressure Implicit with Splitting of Operator)算法相結(jié)合,由于繼承半隱式計(jì)算中的松弛因子修正,并外置多層動(dòng)量方程循環(huán)求解,使得PIMPLE算法在采用大時(shí)間步長處理流態(tài)較復(fù)雜的內(nèi)部場時(shí),仍擁有較好的穩(wěn)定性與計(jì)算精度,相關(guān)算法設(shè)置過程如圖4所示。

圖4 PIMPLE算法計(jì)算流程Fig.4 Flow chart of PMPLE algorithm calculation

為加快計(jì)算收斂速度,將穩(wěn)態(tài)計(jì)算結(jié)果作為瞬態(tài)計(jì)算初始場,相關(guān)物理量最終收斂允許值設(shè)置為10-5,內(nèi)層循環(huán)次數(shù)和網(wǎng)格正交性修正次數(shù)均為3,外層修正次數(shù)設(shè)置為40次,時(shí)間步長設(shè)置為0.001 67 s,在機(jī)組外特性監(jiān)測參數(shù)穩(wěn)定后,額外計(jì)算至少6個(gè)轉(zhuǎn)輪周期。

3 模擬結(jié)果分析

3.1 計(jì)算精度驗(yàn)證

通過進(jìn)行多種工況點(diǎn)比較對模擬精度進(jìn)行驗(yàn)證。如圖5所示,總體來看PIMPLE算法和FLUENT軟件中的SIMPLE算法所得外特性參數(shù)結(jié)果均與試驗(yàn)數(shù)據(jù)相接近,但低負(fù)荷工況下PIMPLE算法計(jì)算結(jié)果更貼近于試驗(yàn)數(shù)據(jù),因此,可以考慮將該算法用于探究水泵水輪機(jī)低負(fù)荷工況流場計(jì)算中。

圖5 模擬精度驗(yàn)證Fig.5 Simulation accuracy verification

3.2 轉(zhuǎn)輪內(nèi)部流態(tài)分析

轉(zhuǎn)輪附近壓力脈動(dòng)監(jiān)測點(diǎn)布置情況如圖6所示,沿?zé)o葉區(qū)周向設(shè)置5個(gè)GV監(jiān)測點(diǎn),并在固定導(dǎo)葉出口處設(shè)置5個(gè)SV監(jiān)測點(diǎn),觀察無葉區(qū)壓力脈動(dòng)傳播時(shí)對上游流態(tài)的干擾程度。

圖6 壓力脈動(dòng)監(jiān)測點(diǎn)布置Fig.6 Layout of pressure pulsation monitoring points

以z=0作為基準(zhǔn)面,在無葉區(qū)靠近上冠和下環(huán)位置分別布置壓力脈動(dòng)監(jiān)測點(diǎn)GU、GD,探究不同高程無葉區(qū)截面處壓力脈動(dòng)幅值的波動(dòng)情況,并進(jìn)行監(jiān)測點(diǎn)處壓力無量綱系數(shù)及轉(zhuǎn)頻分析,其中壓力無量綱系數(shù)定義為

(3)

如圖7所示為監(jiān)測點(diǎn)在不同工況下的頻域圖變化情況。由于無葉區(qū)壓力脈動(dòng)幅值及頻率主要受到轉(zhuǎn)輪與活動(dòng)導(dǎo)葉之間的動(dòng)靜干涉作用,可觀察到GV、SV等監(jiān)測點(diǎn)均存在整數(shù)倍葉頻fn的高幅值現(xiàn)象,而在向上游傳播時(shí),固定導(dǎo)葉區(qū)域高頻幅值及低頻脈動(dòng)幅值大量減弱。并且GU、GD兩個(gè)監(jiān)測點(diǎn)壓力脈動(dòng)波動(dòng)情況與GV變化趨勢一致,在低流量工況0.43QN、0.54QN、0.66QN時(shí)均存在著以0.05fN為主的低頻脈動(dòng)現(xiàn)象,而隨著工況點(diǎn)逐漸接近于額定流量,壓力脈動(dòng)在整數(shù)倍葉頻處幅值開始減弱,無葉區(qū)內(nèi)部流態(tài)得到改善。

由于低負(fù)荷工況下轉(zhuǎn)輪進(jìn)口處易形成高速水環(huán)[21],如圖8中紅色區(qū)域所示。機(jī)組流量為0.43QN時(shí),旋渦結(jié)構(gòu)從葉片吸力面生成,并貼附于流道內(nèi),而此時(shí)內(nèi)部葉道渦結(jié)構(gòu)從進(jìn)口處發(fā)展,延伸至葉片根部堵塞進(jìn)口主流的運(yùn)動(dòng)。因此在轉(zhuǎn)輪進(jìn)口邊及出口邊均可觀察到分布較為分散的回流渦結(jié)構(gòu),并在壁面負(fù)壓區(qū)的共同作用下,轉(zhuǎn)輪進(jìn)口水體形成明顯的脫流現(xiàn)象,使得進(jìn)口主流速度方向偏離最優(yōu)進(jìn)口角度,造成更大的能量損失,以上是造成無葉區(qū)低頻脈動(dòng)現(xiàn)象的主要原因。而隨著機(jī)組過流量增大后,葉道內(nèi)流線相對平滑光順未生成旋渦結(jié)構(gòu)。

圖8 導(dǎo)葉與轉(zhuǎn)輪內(nèi)部速度流線Fig.8 Velocity streamlines diagram inside vanes and runner

3.3 尾水管內(nèi)部流態(tài)分析

如圖9所示,在尾水管進(jìn)口處及距離轉(zhuǎn)輪出口距離為0.5D1和D1的位置,分別設(shè)置監(jiān)測截面s0、s1、s2,并在尾水管進(jìn)口s0截面上設(shè)置Line1監(jiān)測線,用于分析尾水管進(jìn)口處壁面徑向距離速度分布。

圖9 尾水管監(jiān)測截面布置Fig.9 Monitoring section layout in draft tube

機(jī)組在0.43QN工況時(shí),尾水管渦帶形狀為貼近尾水管壁面的帶狀結(jié)構(gòu),如圖10所示。而隨著活動(dòng)導(dǎo)葉開度增大,尾水管內(nèi)部偏心螺旋狀渦帶[22]逐漸轉(zhuǎn)化為同心圓柱形渦帶,而在靠近額定流量時(shí)未形成明顯渦帶形狀。由于渦帶中心區(qū)壓強(qiáng)較低,且渦帶運(yùn)動(dòng)頻率與轉(zhuǎn)輪轉(zhuǎn)頻相關(guān),可發(fā)現(xiàn)隨著流量升高,尾水管不同高程截面上渦帶中心低壓區(qū)[23]分布呈現(xiàn)向中心偏移的趨勢,并且低壓區(qū)覆蓋面積也隨之減小,尾水管內(nèi)部壓強(qiáng)逐漸分布均勻。

圖10 尾水管進(jìn)口處Line1速度分量變化Fig.10 Velocity component in the Line1 of the draft tube inlet

機(jī)組處于0.43QN工況時(shí),尾水管進(jìn)口處軸向速度Ua及周向速度Ut幅值波動(dòng)較大,而徑向速度主要與機(jī)組轉(zhuǎn)速相關(guān),因此變化相對不明顯;截面中心處Ua軸向速度及Ut周向速度分量由于尾水旋渦的存在遞減至負(fù)值,且在低流量工況下各速度分量均存在著明顯的振蕩現(xiàn)象。而在0.66QN工況時(shí),轉(zhuǎn)輪出口處周向速度分量起主導(dǎo)作用,并伴隨著軸向速度分量幅值波動(dòng),進(jìn)而在該工況下可觀察到體積較大的偏心螺旋細(xì)條狀渦帶,而渦帶在運(yùn)動(dòng)過程中會(huì)不斷撞擊管壁,使得進(jìn)口處流態(tài)相對紊亂,容易產(chǎn)生流動(dòng)分離及旋渦現(xiàn)象。當(dāng)流量增加至0.8QN時(shí),尾水管內(nèi)部可觀察到圓柱形同心狀渦帶,并且此時(shí)尾水管進(jìn)口處各速度分量處于向平緩發(fā)展的趨勢,未出現(xiàn)回流結(jié)構(gòu),截面處低壓區(qū)分布也明顯改善,隨著流量增大達(dá)到0.92QN工況時(shí),尾水管內(nèi)部未發(fā)現(xiàn)渦帶存在,徑向處各速度分量幅值波動(dòng)相對平穩(wěn),尾水管出口流態(tài)受干擾較小,因此,機(jī)組在該工況下運(yùn)行時(shí)具有較高的水力效率。

4 結(jié) 論

(1) 本文采用開源軟件OpenFOAM中PIMPLE算法,對水泵水輪機(jī)在多種低負(fù)荷工況的運(yùn)行過程進(jìn)行數(shù)值模擬,得出無葉區(qū)內(nèi)部低頻脈動(dòng)主要受到轉(zhuǎn)輪進(jìn)口處高速水環(huán)及活動(dòng)導(dǎo)葉尾部流態(tài)結(jié)構(gòu)的影響,且低頻脈動(dòng)在向上游傳播后幅值大量削減,對上游固定導(dǎo)葉及蝸殼內(nèi)部流場干擾較小。

(2) 尾水管渦帶的存在形式影響著內(nèi)部流場分布。機(jī)組在低負(fù)荷工況運(yùn)行時(shí),轉(zhuǎn)輪出口水體以圓周速度分量為主,并伴隨著較為劇烈的軸向速度幅值波動(dòng),此時(shí)尾水管直錐段處易形成偏心螺旋條狀渦帶。而隨著流量逐漸增加,尾水管進(jìn)口處各速度分量相對穩(wěn)定,未形成明顯的渦帶現(xiàn)象。

(3) 由于SAS-SST湍流模型在一定范圍內(nèi)仍屬于RANS求解方法,考慮到RANS方法由于在分離流動(dòng)較大的區(qū)域內(nèi)渦識(shí)別的局限性,其計(jì)算誤差較大,后期將結(jié)合OpenFOAM軟件中不同湍流模型的應(yīng)用來獲取更精準(zhǔn)的機(jī)組運(yùn)行狀態(tài)。