基于PIMPLE算法的低負荷水泵水輪機內部流態研究

王 浩 博,周 大 慶,郭 俊 勛

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

0 引 言

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

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

本文采用開源軟件OpenFOAM中的PIMPLE算法對不同低負荷工況下水泵水輪機運行過程進行數值模擬,揭示了機組在低負荷工況下內流場特性和無葉區各頻率壓力脈動的變化規律,探究了尾水管內部的渦帶運動對進口流態分布的影響程度。本文相關研究結果可為解決水泵水輪機運行穩定性問題提供參考。

1 計算域分布和網格劃分

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

圖1 水泵水輪機計算域Fig.1 Calculation domain of pump turbine

機組額定出力為139 MW,額定水頭為105.8 m,相關機組尺寸數據和運行參數如表1所列。

表1 水泵水輪機參數Tab.1 Pump turbine parameters

研究對象采用ICEM軟件對轉輪、活動導葉及進出口管道延伸段進行六面體結構化網格劃分,蝸殼和固定導葉采用四面體網格進行劃分,并基于OpenFOAM中內置前處理工具snappyHexMesh對尾水管進行網格劃分和局部加密。采用6種不同網格數量分配方案對計算域進行網格無關性驗證,網格總數量與機組效率變化趨勢如圖2所示。

圖2 網格無關性驗證Fig.2 Mesh independence verification

當網格數量達到600萬后,效率變化趨近于平緩,考慮到計算時間消耗和湍流場近壁面網格精度,最終網格總量確定為682.3萬。圖3為各計算區域網格數量分配結果:進出口管道延伸段為25.8萬、蝸殼區域為116.8萬、導葉區域為150.9萬、轉輪區域為282.6萬、尾水管區域為106.2萬。其中在轉輪及活動導葉等高雷諾數區域控制網格節點y+值不高于15,符合工程中湍流模型應用要求。

2 數值求解方法

2.1 控制方程

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

(1)

?·u=0

(2)

式中:u為流速,m/s;p為壓強,Pa;ρ為流體密度,kg/m3;t為時間,s;?為哈密頓算子;?2為拉普拉斯算符;ν為運動黏度系數,m2/s。

2.2 初始條件和湍流模型

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

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

2.3 求解算法

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

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

為加快計算收斂速度,將穩態計算結果作為瞬態計算初始場,相關物理量最終收斂允許值設置為10-5,內層循環次數和網格正交性修正次數均為3,外層修正次數設置為40次,時間步長設置為0.001 67 s,在機組外特性監測參數穩定后,額外計算至少6個轉輪周期。

3 模擬結果分析

3.1 計算精度驗證

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

圖5 模擬精度驗證Fig.5 Simulation accuracy verification

3.2 轉輪內部流態分析

轉輪附近壓力脈動監測點布置情況如圖6所示,沿無葉區周向設置5個GV監測點,并在固定導葉出口處設置5個SV監測點,觀察無葉區壓力脈動傳播時對上游流態的干擾程度。

圖6 壓力脈動監測點布置Fig.6 Layout of pressure pulsation monitoring points

以z=0作為基準面,在無葉區靠近上冠和下環位置分別布置壓力脈動監測點GU、GD,探究不同高程無葉區截面處壓力脈動幅值的波動情況,并進行監測點處壓力無量綱系數及轉頻分析,其中壓力無量綱系數定義為

(3)

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

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

圖8 導葉與轉輪內部速度流線Fig.8 Velocity streamlines diagram inside vanes and runner

3.3 尾水管內部流態分析

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

圖9 尾水管監測截面布置Fig.9 Monitoring section layout in draft tube

機組在0.43QN工況時,尾水管渦帶形狀為貼近尾水管壁面的帶狀結構,如圖10所示。而隨著活動導葉開度增大,尾水管內部偏心螺旋狀渦帶[22]逐漸轉化為同心圓柱形渦帶,而在靠近額定流量時未形成明顯渦帶形狀。由于渦帶中心區壓強較低,且渦帶運動頻率與轉輪轉頻相關,可發現隨著流量升高,尾水管不同高程截面上渦帶中心低壓區[23]分布呈現向中心偏移的趨勢,并且低壓區覆蓋面積也隨之減小,尾水管內部壓強逐漸分布均勻。

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

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

4 結 論

(1) 本文采用開源軟件OpenFOAM中PIMPLE算法,對水泵水輪機在多種低負荷工況的運行過程進行數值模擬,得出無葉區內部低頻脈動主要受到轉輪進口處高速水環及活動導葉尾部流態結構的影響,且低頻脈動在向上游傳播后幅值大量削減,對上游固定導葉及蝸殼內部流場干擾較小。

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

(3) 由于SAS-SST湍流模型在一定范圍內仍屬于RANS求解方法,考慮到RANS方法由于在分離流動較大的區域內渦識別的局限性,其計算誤差較大,后期將結合OpenFOAM軟件中不同湍流模型的應用來獲取更精準的機組運行狀態。