范志鋒，林巧鋒，黃葉雯，蘆 月，陳齊燈，鄭 源

（1．福建華電電力工程有限公司，福建省福州市 350001；2．中國華電集團有限公司福建分公司，福建省福州市 350001；3．河海大學，能源與電氣學院，江蘇省南京市 211100；4．河海大學，水利水電學院，江蘇省南京市 210098；5．福建華電金湖電力有限公司，福建省將樂縣 353399）

0 引言

貫流式機組在低水頭、大流量及潮汐能開發(fā)利用等方面具有其他機組無法比擬的優(yōu)勢[1-3]，其中，燈泡貫流式機組因具有結構緊密、水力損失小及建設周期短的特點，近年來發(fā)展迅猛、應用廣泛[4-6]，對機組性能和穩(wěn)定性的研究得到越來越多的重視[7-10]。目前研究貫流式機組水力性能的方法主要是試驗和數(shù)值模擬[11]，但由于貫流式機組流道較短、水體慣性占比較大，模型試驗測得的運行參數(shù)與真機誤差較大，而現(xiàn)場試驗又費時費力，無法反映機組內(nèi)部流態(tài)，因此，數(shù)值模擬是研究貫流式機組水力性能的主要手段。

此外，由于貫流式機組運行水頭低，縱橫尺度相當，水頭與轉輪直徑之比較小，所以上下游水庫中的自由液面以及水體重力因素對機組內(nèi)部流動的影響較其他類型水輪機更明顯，研究時要考慮自由液面對其水力性能帶來的影響[12-13]。在自由液面方面，國內(nèi)外學者開展了相關研究，周大慶[14]等采用VOF模型對抽水蓄能電站水輪機工況進行數(shù)值模擬，結果表明：基于VOF模型的數(shù)值計算具有較好的計算精度，能更真實地反映機組的實際運行條件；Benchikh等[15]基于VOF方法對Darrieus水平軸水輪機進行了多相數(shù)值計算，結果表明：多相模型模擬結果與試驗相比誤差更小，且能更真實地預測機組功率和轉矩系數(shù)；趙亞萍[16]等對三葉片臥式貫流式水輪機進行數(shù)值模擬，在模擬時通過求解體積函數(shù)的輸移擴散方程來確定自由水面的位置，得出需要考慮自由液面才能滿足幾何相似及流動相似的特性；邵方凱[17]對漂浮式垂直軸水輪機建立了在自由液面下周期性旋轉模型，結果表明：自由液面的存在會在低速比下增加機組功率，高速比下則會降低?？梢姡杂梢好鎸λC械數(shù)值模擬的結果存在不同程度的影響，而此類研究中對燈泡貫流式機組的探討較少，且大多沒有考慮水體重力和上下游水池帶來的影響，少數(shù)考慮到此因素的研究也只是把水池表面設置為無滑移壁面[18-21]。

基于上述研究現(xiàn)狀，本文以福建某燈泡貫流式水輪機為研究對象，構建包含上下游水池在內(nèi)的全過流系統(tǒng)三維模型，在考慮水體重力的基礎上，針對水池表面與大氣相通的特點，結合VOF模型對不考慮自由液面和考慮自由液面兩種情況分別進行數(shù)值模擬，以更符合機組實際運行狀態(tài)，并從能量特性、不同測點壓力脈動及內(nèi)部流場特性等方面對比單相工況和兩相工況的差異，揭示上下游水池自由液面對燈泡貫流式機組水力性能的影響。

1 計算模型與計算方法

1．1 機組參數(shù)與模型試驗

本文以福建某燈泡貫流式機組模型裝置為研究對象，該機組水輪機型號為GZTF07B-WP-560，額定功率為21MW，額定轉速為93.75r/min，額定水頭為9.3m，最大水頭11.9m，最小水頭4.7m，額定出力為21.875MW，額定流量為255.86m3/s，葉片數(shù)為4，活動導葉數(shù)為16。該機組的模型水輪機試驗是在哈爾濱大電機研究所的高水頭試驗I臺上進行的。試驗臺模型效率綜合測試誤差小于±0.25%，模型效率重復測試誤差小于±0.10%。模型的效率、空化、壓力脈動等試驗能在同一試驗臺同一模型上進行，試驗臺主要測量儀器設備見表1所示。

表1 模型試驗臺主要測量儀器設備Table 1 Main measuring instruments and equipment of model test bench

模型試驗中的能量試驗是在水頭為常數(shù)、無空化條件下進行的，每個工況點采集時間為60s，試驗包括整個水輪機運行范圍：模型試驗水頭為4～2.5m，模型活動導葉開度為24°～90°，槳葉開度為5°～45°，根據(jù)模型與原型換算公式，可得到原型水輪機能量試驗結果，其中，在額定工況下，即水頭9.3m、導葉開度為68.5°、槳葉開度為36.1°時，穩(wěn)定運行時對應的機組流量為260.19m3/s，轉輪力矩為2086.2 kN·m，整機效率為86.2%。

1．2 計算域及網(wǎng)格劃分

對該機組進行額定工況下三維數(shù)值模擬計算，利用UG軟件對水輪機全流道進行同等尺寸建模，三維模型包括上游進水池、進水流道、導葉、轉輪、出水流道、下游出水池等六部分構成，上下游水池高度差為額定水頭，進出水池均為三面過水，圖1為燈泡貫流式水輪機全流道結構示意圖。在考慮自由液面數(shù)值模擬時，僅在上下游水池的頂端加入一定高度的空氣域，空氣層高度與液體高度接近1：1。

圖1 燈泡貫流式水輪機全流道結構示意圖Figure 1 Structural diagram of the whole passage of bulb tubular turbine

本文利用ICEM軟件對模型進行網(wǎng)格劃分，在不考慮自由液面模擬時，采用非結構化四面體網(wǎng)格對進出水池、進出水流道進行劃分，對流態(tài)復雜的導葉和轉輪區(qū)域進行網(wǎng)格加密。在考慮自由液面模擬時，由于氣體和水體在密度上有較大的差異，為了能較為清晰地捕捉兩相介質在自由液面處的相界面，需要對上下游水池自由液面交界面進行網(wǎng)格加密，因此在上下游水池部分采取結構化網(wǎng)格進行劃分，圖2為考慮自由液面水輪機各部分網(wǎng)格劃分示意圖。

圖2 考慮自由液面水輪機各部分網(wǎng)格劃分示意圖Figure 2 Schematic diagram of grid division of each part of hydraulic turbine considering free surface

在不考慮自由液面的額定工況模擬中進行網(wǎng)格無關性驗證，所以現(xiàn)采取5種不同數(shù)量的網(wǎng)格劃分情況，總數(shù)依次為550萬、625萬、790萬、845萬以及910萬，對機組進行額定工況下的穩(wěn)態(tài)數(shù)值模擬計算，具體計算網(wǎng)格數(shù)及其效率如表2所示，從表中可以發(fā)現(xiàn)網(wǎng)格數(shù)超過845萬之后機組的效率變化已經(jīng)很小，出于對時間成本和計算精度的考慮，最終選取方案4的網(wǎng)格數(shù)進行后續(xù)計算。

表2 不同網(wǎng)格劃分方案及機組額定工況效率模擬結果Table 2 Simulation results of different grid division schemes and unit efficiency under rated working conditions

1．3 考慮自由液面VOF模型

本文基于Fluent19.2軟件平臺VOF多相流模型進行三維湍流數(shù)值模擬計算，VOF模型是建立在固定的歐拉網(wǎng)格下的表面跟蹤法，其特點在于通過建立體積函數(shù)來表征流體所占單元的體積分數(shù)，以此創(chuàng)建和追蹤自由表面。在單元中，若q相流體體積分數(shù)為αq，那么可能存在三種情況，一是αq為0即單元里不存在第q相流體，二是αq為1即單元里充滿了第q相流體，三是0＜αq＜1即單元里包含了第q相流體和其他單相或者多相流體的界面。在本文中，分別定義αa和αw為貫流式機組上游水池和下游水池內(nèi)氣體和液體所占的體積分數(shù)，則氣體的體積分數(shù)可用液體的體積分數(shù)表示，即：

其中，液體的體積率函數(shù)αw的控制方程表達式為：

引入VOF方法的k-ε湍流模型與單相流的k-ε湍流模型形式完全相同，只有密度ρ和動力黏度μ的具體表達式不同，在此系統(tǒng)中的物性參數(shù)均由體積率來控制，其表達式如下：

式中：ρw——液體相的密度；

ρa——空氣相的密度；

μw——液體相的分子黏性；

它以特殊鋁材為材料，可應用于幾乎所有行業(yè)的輸送機結構中，運輸絕大多數(shù)種類的貨物。它被組裝在OTU換向臺上，并配有電源、B L D C電機和控制器，可根據(jù)特定的應用程序進行設置。

μa——空氣相的分子黏性。

通過對液體的體積分數(shù)αw的求解，密度ρ和動力黏度μ即可通過式（3）和式（4）得出。

1．4 計算工況及邊界條件

不考慮自由液面：選取RNGk-ε湍流模型對控制方程進行封閉，在計算中方程組壓力項為二階中心差分格式，速度項、湍動能和湍動能耗散率均為二階迎風格式，求解方法選定為SIMPLEC算法，轉輪區(qū)域設置為滑移網(wǎng)格，固定旋轉速度為額定轉速93.75r/min，固體壁面設置為無滑移邊界，臨近固體壁面處采取標準壁面函數(shù)，不考慮壁面粗糙度對流場的影響。最大迭代步數(shù)為20步，收斂殘差為10-5，設置葉片旋轉每三度所經(jīng)過的時間作為步長，約為0.005s。上游水池進口設置壓力進口，下游水池出口設置壓力出口，考慮到重力場對水頭較低的貫流式機組內(nèi)部流態(tài)的影響，通過Fluent中的用戶自定義函數(shù)設置上下游水池進出口壓力沿水池高度變化。

考慮自由液面：打開VOF模型，添加明渠，上下游水池進出口均給定液面高度，水池頂部空氣域設為開放面，空氣層區(qū)域水相分數(shù)為0，即充滿空氣；其余各部分水的體積分數(shù)設置為1，其他條件與不考慮自由液面時相同，考慮自由液面后的流體初始狀態(tài)如圖3所示。

圖3 考慮自由液面水輪機初始流場Figure 3 Initial flow field of hydraulic turbine considering free surface

2 計算結果對比與分析

2．1 能量特性

對于不考慮自由液面和考慮自由液面的工況分別進行三維數(shù)值模擬，從計算開始至16s后計算收斂，各監(jiān)測數(shù)據(jù)逐漸趨于穩(wěn)定值，選取穩(wěn)定后相同計算時間內(nèi)的結果進行對比。表3所示為不考慮和考慮自由液面水輪機的能量特性對比，從表中可以發(fā)現(xiàn)，在機組水頭和轉速相同的條件下，不考慮自由液面時，機組的流量較小、轉輪力矩較大，效率達到91.1%；考慮自由液面時，機組的流量升高了1.79%，轉輪力矩減小了6.65%，效率降低了7.5%，且與模型試驗的結果更接近，誤差更小。這說明貫流式機組在實際運行時，受到重力場和自由液面的影響，機組過流量有少許增加，而轉輪出力有大幅下降，機組效率有明顯降低。這種差異會影響貫流式機組的初期設計，導致機組出力不足與設計值誤差較大的情況發(fā)生。

表3 不考慮和考慮自由液面水輪機的能量特性對比Table 3 Comparison of energy characteristics of hydraulic turbines without and without free surface

2．2 壓力脈動規(guī)律

為探究貫流式機組流道內(nèi)各位置壓力變化情況，分別在導葉前、轉輪前及尾水管進口、出口四個截面上布置若干壓力監(jiān)測點，考慮重力場影響，在每個截面外緣位置以及截面上同一縱向位置分別布置三個監(jiān)測點，監(jiān)測點布置圖如圖4所示。為了使計算結果更準確，采樣時間為12個周期，選取最后兩個周期進行壓力脈動分析。壓力脈動的時域值是通過數(shù)值計算結果直接獲得的，而頻域值則需要進行傅里葉變換才能獲取。本文取無量綱的壓力脈動系數(shù)Cp來反映各過流部位的壓力脈動特性，具體公式為：

圖4 水輪機壓力脈動監(jiān)測點布置圖Figure 4 Layout of pressure fluctuation monitoring points of hydraulic turbine

式中：pi——監(jiān)測點在某一時刻的壓力值；

通過對考慮自由液面模擬下各監(jiān)測點壓力脈動的初步分析，發(fā)現(xiàn)各監(jiān)測截面上，外緣位置監(jiān)測點的波動情況無明顯差異，因此，這里只對比截面上同一縱向位置上的三個點，圖5為考慮自由液面的兩相工況各截面縱向分布的三個點的壓力脈動頻域圖，從頻域特性上看，導葉進口處的壓力主頻均出現(xiàn)在四倍轉頻fn處（4fn），且壓力脈動系數(shù)Cp值都很小，縱向三個點中位于中間位置的A6’點Cp值最高為0.0013，其次頻Cp值為主頻幅值的38%，除此以外出現(xiàn)了很多幅值頻率不定的干擾信號，說明導葉進口在一定程度上出現(xiàn)了低頻壓力脈動，這可能是由于導葉進口處靠近進水池，水位的波動帶來的影響，但整體數(shù)值較小，流態(tài)相對平穩(wěn)；轉輪進口處的壓力主頻也在4fn處，在8fn、12fn、16fn處出現(xiàn)了不同程度的諧波，其壓力波動情況在四個截面中最強烈，三個點中外緣位置的B4′點Cp值最高為0.258，這是因為此位置距離轉輪最近，受到導葉和轉輪的動靜干擾；尾水管進口處的壓力主頻也在4fn處，外緣位置的C4′點Cp值最高為0.023，同時其高頻區(qū)域起伏也最大，倍頻信號幅值呈減小趨勢，說明水流在進入尾水管時對壁面造成了一定的沖擊，三個點中越靠近中心位置，諧波也逐漸消失，表明在尾水管進口處軸心區(qū)域流態(tài)較穩(wěn)定；尾水管出口處的壓力頻率主要在fn處、0.78fn、0.44fn處，這是由于此位置距離轉輪區(qū)域較遠，受到轉輪旋轉的影響較小，此外，該截面低頻壓力脈動占主導地位，其中D4′點主頻Cp值最高為0.056。

圖5 兩相工況各截面縱向點壓力脈動頻域圖Figure 5 Frequency domain diagram of pressure fluctuation at longitudinal points of each section under two-phase working condition

將不考慮自由液面和考慮自由液面的壓力脈動頻域進行比對，選取兩種工況下各截面上波動幅值最高的對應點進行對比，如圖6所示為不考慮和考慮自由液面各截面壓力脈動頻域對比圖，從圖中可以發(fā)現(xiàn)四個截面處單相和多相下的壓力脈動頻域變化規(guī)律基本相同，且主頻均在同一位置，不同的是，氣液兩相工況整體在頻域圖中出現(xiàn)了很多幅值頻率不等的干擾信號，雜頻較多，特別是靠近進水池的導葉進口處和靠近出水池的尾水管出口處最為明顯，這是因為考慮自由液面后機組運行受到空氣的干擾，液面的變動需要經(jīng)歷一段時間后才能逐漸趨于穩(wěn)定。此外，考慮自由液面后四個截面壓力脈動的主頻幅值均增大：導葉進口處的兩相工況主頻Cp值為單相時的1.34倍；轉輪進口處的兩相工況主頻Cp值為單相時的1.27倍，諧波幅值無明顯差異；尾水管進口處的兩相工況主頻Cp值為單相時的1.23倍，且兩相工況在4fn之前區(qū)域出現(xiàn)了部分雜頻；尾水管出口處的兩相兩相工況主頻Cp值為單相時的1.43倍，且均為低頻壓力脈動。綜合對比圖來說，考慮自由液面時，機組上下游水池液面受到空氣層的干擾，水相和氣相交互引起水池內(nèi)流態(tài)波動，這種波動隨后進入流道內(nèi)影響轉輪旋轉，繼而影響整個機組的壓力脈動狀況。

圖6 不考慮和考慮自由液面各截面壓力脈動頻域對比圖Figure 6 Frequency domain comparison of pressure pulsation of each section without and without free surface

2．3 機組內(nèi)部流態(tài)

為了進一步分析加入自由液面后內(nèi)部流場水流分布情況，對進出水池自由液面的水相進行了對比分析。圖7為兩相工況進出水池自由液面處的水相分布。從圖中可以看出，進出水池自由液面水相分布均為對稱分布，不同的是進口處水相體積范圍明顯比出口多，最高值分布在邊緣處，一定程度上呈現(xiàn)梯度變化，而出水池則沒有了規(guī)律，最高值分布在中心區(qū)域。

圖7 進出水池自由液面處水相分布圖Figure 7 Distribution of water phase at the free surface of inlet and outlet pool

如圖8所示為不考慮自由液面和考慮自由液面的機組縱向剖面壓力分布云圖?？梢园l(fā)現(xiàn)，在重力作用的影響下，進出水池的水體部分壓力均呈現(xiàn)沿水深增加的梯度變化規(guī)律，且轉輪前區(qū)域的壓力整體都比轉輪后區(qū)域的壓力數(shù)值大，在轉輪區(qū)域兩個工況壓力沿軸向都有一定的變化；考慮自由液面后，進出水池的底部高壓區(qū)域變大，且與空氣層交界處位置的壓力與不考慮自由液面時相比增大2.4%，說明液相明顯受到了氣相的影響，此外，空氣層部分壓力不沿梯度變化，與實際情況相符。這是因為不考慮自由液面時進出水池被假設為密閉容器，水池頂部設置為固體壁面，而考慮自由液面后進出水池直接與大氣相通，水池頂部設置為開放面，自由液面的存在增加了水池頂部的初始壓力，因而造成水池內(nèi)壓力隨水深的增加而增大，整體大于不考慮自由液面時的水池壓力。

圖8 縱向剖面壓力分布對比云圖Figure 8 Comparison cloud chart of pressure distribution in longitudinal section

如圖9所示為不考慮自由液面和考慮自由液面的導葉進口湍動能分布對比圖。從圖中可以發(fā)現(xiàn)兩種工況下的湍動能分布規(guī)律基本一致，均為對稱分布，且數(shù)值相差很小，兩相工況略高于單相工況；高湍動能區(qū)域都在導葉區(qū)域上部靠近燈泡體壁面的位置，低值區(qū)域均在流道的中間部位；整個截面上湍動能值從外壁面到燈泡體壁面呈現(xiàn)先減小后增大的變化規(guī)律，不考慮自由液面的湍動能等值線較為流暢、過渡均勻，而考慮自由液面后的湍動能等值線有突變情況，且在靠近外壁面處突變更明顯，這可能與考慮自由液面后，氣相和水相發(fā)生碰撞，有能量交換的情況相關。

圖9 導葉進口湍動能分布對比圖Figure 9 Comparison of turbulent kinetic energy distribution at guide vane inlet

如圖10所示分別為考慮自由液面和不考慮自由液面的轉輪葉片吸力面和壓力面的渦核區(qū)域速度分布對比圖。從圖中可以看出單相工況和兩相工況轉輪渦核區(qū)速度分布規(guī)律基本一致，壓力面幾乎沒有渦核分布，吸力面進水側分布較密集，從輪轂處向輪緣處逐漸增大，在葉片外緣側尖端速度達到最大值；考慮自由液面后，兩相工況葉片外緣側的速度值高于單相工況，且高速區(qū)域范圍明顯大于單相工況，說明加入氣相后使轉輪吸力面的渦核速度增加，這是因為自由液面部分氣相和液相的相互作用引起了液面的波動，這種波動會隨液相的流動作用在轉輪部分，而波動產(chǎn)生的能量增加了葉片區(qū)域渦核的動量，從而干擾了轉輪的旋轉，使轉輪部分的水力損失增加，這也是造成轉輪出力下降的原因。

圖10 轉輪葉片進出水面渦核區(qū)域速度對比圖Figure 10 Velocity comparison diagram of runner blade in and out of water surface vortex core area

如圖11所示為不考慮自由液面和考慮自由液面的尾水管內(nèi)速度分布對比云圖。從速度云圖中可以看出兩種工況速度都是呈現(xiàn)對稱分布的，在尾水管入口處速度最高，且最高流速相差不大，同時在尾水管中間部分形成低速區(qū)域，說明此處有同程度的旋渦產(chǎn)生，不考慮自由液面時低速區(qū)域主要集中在中間部分，而加入氣體相后低速區(qū)域范圍擴大且主要分布在尾水管的中部及靠近出口部位，藍色區(qū)域分布更為廣泛，這是由于尾水管靠近下游出水池，受到出水池空氣層的干擾，兩相工況尾水管內(nèi)流態(tài)更為混亂，容易形成擾流和回流，在實際運行時，尾水管部分區(qū)域還可能會出現(xiàn)真空情況。

圖11 尾水管內(nèi)速度分布對比云圖Figure 11 Comparison cloud chart of velocity distribution in draft tube

3 結語

本文以福建某燈泡貫流式水輪機為研究對象，在考慮水體重力的基礎上，采用VOF模型對不考慮自由液面和考慮自由液面兩種情況分別對機組設計工況下的穩(wěn)態(tài)運行過程進行數(shù)值模擬，探究加入氣相后對機組能量特性、壓力脈動及內(nèi)部流場的影響，對比結果表明，考慮自由液面后：

（1）能量特性方面，水輪機的流量比單相工況時提高了1.79%，轉輪力矩減小了6.65%，效率降低了7.5%，與模型試驗的結果更接近，這種差異會影響貫流式機組的初期設計，導致機組出力不足與設計值誤差較大的情況發(fā)生。

（2）壓力脈動方面，各截面監(jiān)測點主頻沒有發(fā)生變化，但壓力脈動幅值均大于單相工況，靠近水池的導葉進口面和尾水管出口面幅值變化較其他面更大，且出現(xiàn)了很多幅值頻率不等的干擾信號。

（3）內(nèi)部流態(tài)方面，進出水池的底部高壓區(qū)域變大，與空氣層交界處位置的壓力比單相時增大2.4%，導葉進口處湍動能等值線有明顯突變，轉輪吸力面渦核高速區(qū)域范圍明顯增大，尾水管內(nèi)低速區(qū)域范圍擴大且向靠近出口部位分布，內(nèi)部流態(tài)更混亂。

綜上所述，當貫流式機組數(shù)值模擬考慮重力項和自由液面后，更符合機組實際運行狀況，受到氣相和水相的交互影響，機組性能也會產(chǎn)生一定的差異，在設計機組和對機組運行狀態(tài)預測時需要考慮自由液面的影響。