一管雙機式混流式水輪機壓力脈動預測與振動分析

楊夢嬌，陳　端，張昌兵，李金遙

(四川大學水利水電學院， 成都 610065)

0　引　言

水力發電機組安全、穩定運行一直是廣大專家和學者重點關注的問題，尤其是機組運行過程中的水力振動問題[1-3]。引起水力振動的原因是多方面的，包括葉輪機械固有的轉動部件與靜止部件間的動靜干擾，葉片繞流產生的卡門渦列，偏離最優工況下尾水管中發展出的旋轉渦帶[4]等。由于水輪機流道的復雜性，采用實驗方法測量其內部流場比較困難，且成本高、周期長。隨著計算機和計算流體力學(CFD)的不斷發展與完善，應用數值模擬方法計算水輪機內部的復雜流場，預測壓力脈動成為一種重要工具，并取得了一系成果[5-8]。楊建明在文獻[9]中指出尾水管渦帶是造成水力機械在部分負荷時引起機組振動和出力擺動的主要根源。文獻[10]中說明水電站過渡過程也可以引起振動。基于振動帶來的危害，蘇華山[11]等人研究了振動條件下液力耦合器葉輪內部兩相流動特性，探索了振動對液力耦合器內外特性影響。曾立飛[12]等人針對調節閥振動問題，考慮閥桿系統振動對閥內流場的影響，進一步揭示由流體誘發振動的原因。龐立軍[13]等人針對三峽右岸部分機組振動、產生異常噪音的問題，模擬了固定導葉開展出水邊處的渦街振動，確定了產生異常噪聲的激振源，并提出有效優化方案。

對于引水管路較長的水電站，為了節約成本，往往采用一根引水管帶兩臺、甚至多臺機組的方式。但機組之間不可避免地會出現水力干擾的問題，嚴重時可能影響機組的穩定運行。文獻[14]報道了某電站一根總管帶四臺水輪機組，機組的振動問題十分突出，大修周期約為2～3年，平均每兩年就要更換一個轉輪，這嚴重影響了電站的正常運行及其綜合效率。文獻[15]對其進行了數值分析，得出了振動主要是由于轉輪引起的結論。文獻[16]對水輪機進行了尾水管補氣數值模擬，有效的補氣方式能有效減輕水輪機的水力振動。

由于一根總管帶多臺水輪發電機組更容易造成水力振動。但到目前為止，水輪機的數值模擬主要是針對單臺機組，沒有考慮機組之間的水力干擾導致的水力因素，因此，有必要對其水力振動原因作進一步分析。研究表明[17]，標準k-ε模型用于應變率特別大時可能導致負的正應力，為了模擬的流動更加符合湍流的物理定律，Realizablek-ε模型將原來標準k-ε模型中不變的黏度系數Cμ與應變率聯系起來，引入與旋轉和曲率有關的內容，因而Realizablek-ε模型能更好地模擬射流撞擊、分離流、二次流、旋流和圓口射流等復雜的流動問題，且計算結果較準確，精確較高[18]。因此，本文以一管兩機混流式水輪機為例，采用時均N-S控制方程組結合Realizablek-ε模型對其水力振動進行數值分析，為機組的穩定運行提供科學依據。

1　數學模型

水輪機的水流按不可壓縮流處理。其中岔管、蝸殼、尾水管的水流采用絕對坐標系下控制方程；轉輪內的水流采用的旋轉坐標系下的控制方程。

1.1　在絕對坐標系統下的控制方程

水輪機的水流視為不可壓縮流動，通過岔管、蝸殼、尾水管的水的連續性方程和動量方程如下：

▽V=0

(1)

(2)

式中：V是水流的絕對速度，m/s；t是時間，s；ρ0是水的密度，kg/m3；p是壓力Pa；ν和λ分別是第一和第二流體黏度系數。

1.2　在旋轉坐標系統下的控制方程

在旋轉坐標系下，通過轉輪的水的連續性方程和動量方程如下：

▽W=0

(3)

(4)

式中：W是水流的相對速度，m/s；ω是旋轉角速度，rad/s；r是半徑，m；p是壓力Pa。

1.3　Realizable k-ε模型

Realizablek-ε模型的方程如下[19]：

湍動能k輸運方程：

湍流耗散率ε輸運方程：

(6)

式中：ρ為流體密度；t為時間；k為湍動強度；ε為耗散率；μ為動力黏度系數；v為運動黏度系數；μt為湍流黏性系數；ui、uj分別為流體在xi、xj方向的瞬時流速；Gk為平均速度梯度引起的湍動能k的產生項；S為平均應變率張量系數；G2為常數；σk、σε分別為k、ε的普朗特數。

湍流黏性系數計算：

(7)

上述公式中，取C2=1.9，σk=1.0，σε=1.2。

2　計算方法

2.1　計算域與網格

計算域包括一根岔管和兩臺相同型號的混流式水輪機。水輪機包括蝸殼、轉輪和尾水管，如圖1所示。水輪機型號為HLD126-LJ145，固定導葉數8個，活動導葉數Z1=16個，轉輪葉片數Z2=14，尾水管形狀為彎肘型，其參數如表1所示。

表1　水輪機參數Tab.1 Turbine parameters

本文對水輪機從岔管入口至尾水管出口的所有部件進行全流道三維計算，運用基于壓力法的Pressure-Based求解器，動量方程中壓力與速度的耦合采用廣泛應用的SIMPLE[19]算法。采用非結構化網格，在空間離散上，動量方程采用二階迎風分格式，湍流量輸運方程采用一階迎風分格式， 壓力項離散具有二階精度。在時間上，采用一階隱式格式進行離散。

由于結構比較復雜，采用適用性強的非結構化網格進行劃分，如圖1。轉動部件和非轉動部件之間采用的交界面(interface)模式為：動靜部件耦合，采用瞬態計算的方法。

圖1　計算區域及網格Fig.1 Computational domain and mesh

為了保證計算模型能真實模擬水流內部的流動過程，且滿足精度要求，在比較復雜的部件如轉輪處采用細網格，在相鄰部件連接處進行了均勻細化處理，使得網格變化均勻過渡。網格劃分總數為250萬。

2.2　邊界條件

流體狀態：流體選用液態水，狀態為湍流，采用不可壓縮黏性模型，密度為998.2 kg/m3，黏度為0.001 kg/(m·s)。

入口條件：入口采用速度進口邊界條件，采用平均流速。

出口條件：出口采用壓力出口條件。根據水輪機的吸出高度Hs=-8.0 m，設定尾水管出口壓力為8.0 m。

壁面條件：固壁面采用無滑移邊界條件(wall)，近壁區采用標準壁面函數來確定其附近區域的流動。

“壁面函數”適用于實體邊界：

(8)

絕對和旋轉坐標系之間的界面條件：

V=W+ωr

(9)

式中：ωr是圓周速度, m/s.

2.3　計算工況

為了研究不同開度和水頭下機組振動的情況，對表2所示工況進行計算。

表2　計算工況Tab.2 Computational operation conditions

為了更好地研究內部水流的變化情況，通過在各過流部件內設置監測點，計算其在不同工況下的壓力脈動情況。監測點的位置如圖2所示。圖2(a)中P1～P5為Z=0平面上的點；圖2(b)中P6～P10為1號機組上Z=-0.6 m面上分布的點，其中P6為該截面上的中點，P7～P10為繞P6周向分布、距離中心0.6 m的點；P16為Z=-1.2 m平面上的中點。(P11～P15、P17為2號機組上與1號機組對應位置上的監測點)

圖2　監測點位置Fig.2 The location of monitoring points

3　結果與分析

先進行定常計算，當計算收斂時，將定常計算的結果作為非定常計算的初始條件。非定常計算時，依據收斂條件，時間步長應滿足Courant-Friedrichs-Lewy(CFL)條件：CCFL=uΔt/Δx, 對于水輪機模擬計算，取CCFL≤1，并由最小網格尺寸和最大流速，時間步長應小于5×10-4s，計算時取1×10-4s。計算過程中，對選定測點位置的壓力脈動監測，并進行頻譜分析，分析水力振動的原因。

3.1　轉輪葉片背面壓力分析

以1號轉輪為例，其葉片背面壓力分布如圖3所示。

由圖3可知，從葉片進水邊至出水邊，壓力呈減小趨勢，轉輪葉片靠近上冠處的壓力分布比靠近下環位置的壓力分布更均勻。葉片背面進水邊靠近下環位置處壓力梯度較大，由圖3(f)可見，該處產生了負壓，表明在大開度、高水頭工況下，該處水流流速變化大，導致較大的壓力梯度，對葉片造成沖擊與振動，引起疲勞破壞。下環靠近出水邊存在大范圍的負壓區，該區域容易產生氣蝕破壞，影響機組的穩定運行。計算結果表明，在同一水頭下，開度越大，下環出水邊負壓區的范圍越大，其表面最大壓力值越大。在同一開度下，水頭越高，表面的最大負壓值越大。2號轉輪的葉片壓力分布規律與1號轉輪葉片壓力分布規律大致相同，但2號機組負壓區的范圍更大，表明兩臺水輪機的工作狀態存在差異。

圖3　轉輪葉片背面壓力分布圖Fig.3 Pressure contour on blade back surface of runner

3.2　壓力脈動

為了進一步分析壓力脈動的變化情況，對各個監測點壓力隨時間的變化進行快速傅里葉變換(FFT)。通過比較，選擇蝸殼進口處的P4、P5，截面Z=-0.6 m 的中心P6、P11和在該截面上離中心0.6 m的點P7、P12幾個典型監測點進行分析。圖4為工況1下的壓力脈動計算結果和對應的頻譜分布情況，脈動壓力幅值及對應的功率譜見表3。

由圖4中壓力脈動變化過程和表3中壓力脈動幅值可見。蝸殼中的和尾水管中心線的壓力脈動周期較短，而尾水管渦帶掠過的區域，壓力脈動周期相對較長。由表3可見，尾水管中心線處的壓力脈動幅值較小，而渦帶掠過的區域壓力脈動幅值很大，其值是蝸殼里壓力脈動幅值的5倍左右，是尾水管中心線處的20多倍。表中監測點4、6、7(屬于1號機組)的幅值和功率密度值分別大于對應監測點5、11、12(屬于2號機組)。說明1號機組的振動強于2號機組，且振動產生的能量(功率密度的大小代表能量高低)也高于2號機組，這與實際相符合。并且1號機組P6的功率密度的值是對應2號機組P11的值的數倍，在轉輪出口處兩臺機組產生了較大差距。這可能是由于在轉輪內，兩臺機組的水流的流態不同，產生的壓力脈動不同，導致了進入尾水管的水流流態不同，也就造成了二者在尾水管能量耗散的差異性。壓力脈動的不穩定，引起水流的沖擊與振動，故而導致廠房出現振動，產生異常噪聲。這正是引起兩臺機組振動不同的原因。

圖4　壓力脈動時域圖和功率密度頻域圖Fig.4 Time process of the pressure pulsation and power spectrum density (PSD)

表3　壓力脈動的幅值及最大功率譜密度Tab.3 Amplitude of pressure pulsation and maximum power spectral density

注：壓力脈動幅值表示峰-峰值。

由頻譜分析結果可知，發現兩臺機組對應監測點的主頻基本一致，說明相同位置受影響的主要因素相同。蝸殼壓力脈動的主頻為140.0 Hz，這與水輪機轉輪葉片的固有頻率f2=nZ2/60=140 Hz相同，表明蝸殼里壓力脈動主要受水輪機轉輪葉片動靜干涉的影響，另外也包含一些低頻成分；尾水管中心線處壓力脈動頻率主要為10.0 Hz和140.0 Hz，與水輪機固有頻率f1=n/60=10 Hz和轉輪葉片固有頻率相同， 主要受到轉輪的影響。在尾水管渦帶掠過的區域，壓力脈動的主頻主要為10.0 Hz及其倍頻。

由表3中的最大功率譜密度值來看，尾水管中心線位置的功率譜密度值很小，表明該處的壓力脈動所含能量較小。蝸殼里的最大功率譜密度值在10.0 MPa2·s左右，而尾水管渦帶區域的最大功率譜密度值很大，達到180.0 MPa2·s，是蝸殼里最大功率譜密度值的18倍，表明尾水管渦帶單位時間所含振動能量較大，這是造成機組振動的主要因素。

3.3　尾水管渦帶

幾種典型工況下的尾水管渦帶如圖5所示，圖5中左邊為1號機組，右邊為2號機組。

圖5　穩定時刻尾水管內三維渦帶圖Fig.5 Three-dimensional Vorex rope in the draft tube at a certain stabilization time

當機組偏離最優工況區運行時，由于轉輪葉片出口產生圓周方向的速度分量，會引起尾水管內水流的旋轉，產生渦帶，進而產生能量損失，影響機組的穩定工作[20]。由圖5可知，工況不同，尾水管渦帶的形狀有差別。同一臺機組，在同一開度下，隨著水頭的增加，流量增大，渦帶逐漸減小，渦帶強度減弱，能量耗散減小，振動減小。在設計工況(即工況7)下無漩渦。比較兩臺水輪機，在相同水頭、相同開度下，尾水管進口里的渦帶是相似的，對應的功率譜密度值也基本相等。

因此可以發現，兩臺不同的機組，在相同水頭、相同開度下，兩臺聯合運行機組的尾水管進口處的渦帶差別并不明顯，可以推測尾水管處渦帶產生的振動并不是兩臺機組產生差異的癥結所在。

4　結論與展望

本文通過對一根引水管帶兩臺水輪機的非恒定流場進行了三維數值預測，可以得到如下結論：

(1)采用不可壓縮流N-S方程結合real Realizablek-ε模型較好地模擬了一管帶雙水輪機的流場，得到了水輪機的壓力分布和渦帶分布。

(2)轉輪葉片背面靠近下環出水邊出現負壓區，同一水頭下，開度越大，最大負壓值越大；同一開度下，水頭越高，最大負壓值越大，越容易產生空蝕破壞，也是引起機組產生振動的重要原因之一。計算結果表明，兩臺機組轉輪葉片負壓區的面積和最大負壓值有差異，導致兩臺機的空蝕破壞和振動程度不相同。

(3)由壓力脈動的頻譜分析可知，蝸殼里壓力脈動主頻為140.0 Hz，與水輪機和轉輪葉片固有頻率相同，主要受轉輪葉片的動靜干擾的影響；尾水管里壓力脈動的影響因素較多，除了和轉輪葉片頻率影響外，還包括水輪機固有頻率10.0 Hz及其倍頻、水輪機導葉頻率160.0 Hz等其他因素的影響。影響最大的還是尾水管渦帶掠過的區域，其最大功率譜密度值是蝸殼里的18倍，表明尾水管蝸帶單位時間所含振動能量較大，這是造成機組振動的主要因素。

(4)在設計開度下，沒有尾水渦帶出現。隨著活動導葉開度減小, 尾水管渦帶由轉輪泄水錐底部產生，并逐漸發展，在25°開度下發展為相對于尾水管中心線對稱的三根渦帶，并且渦帶強度也逐漸增大。在活動導葉開度從25°減小到15°的過程中，尾水渦帶又逐漸過渡到一根渦帶，形成偏心渦帶，渦帶強度進一步增大。模擬結果與實際情況吻合。

(5)由于兩臺機組之間尾水管進口振動因子存在差異化，說明兩臺機組之間振動原因的不同，振動結果不同，導致1號、2號機組的振動情況不同；同一臺機組，在相同開度下，隨著水頭增加，尾水管處渦帶減小，振動減小，但在相同工況下，兩臺聯合運行機組的尾水管進口處的渦帶差別卻不大。說明引起兩臺機組振動的差異的位置主要發生在轉輪處而非尾水管的渦帶，但振動的主要因素卻是由尾水管渦帶引發的。

(6)盡管采用Real Realizablek-ε模型比較成功地模擬了一管多機式混流式水輪機的渦帶和壓力脈動特性，但由于該模型本身的時均性，兩機振動的差異并不明顯。在本文的基礎上，有必要進一步對模型進行研究和對比，以及開展流-固耦合的多物理場的研究，以便能更好地模擬出機組之間振動的差異和規律，為機組穩定運行和維護提供科學依據。

□

參考文獻：

[1]陳婧,馬震岳,劉志明,等.水輪機壓力脈動誘發廠房振動分析[J].水力發電,2004,30(5):24-27.

[2]肖若富,韋彩新,韓風琴,等.混流式水輪機轉輪的動力學研究[J].大電機技術,2001:41-43.

[3]王正偉,周凌九,黃源芳.尾水管渦帶引起的不穩定流動計算與分析[J].清華大學學報 (自然科學版), 2002,42(12):1 647-1 650.

[4]劉宇,楊建明,戴江,等.混流式水輪機三維非定常湍流計算[J].水力發電學報,2004,23(4):102-105.

[5]ZUO Zhi-gang, LIU Shu-hong, LIU De-min et al.Numerical analyses of pressure fluctuations induced by interblade vortices in a model Francis turbine[J].Journal of Hydrodynamics, 2015,27(4):513-521.

[6]XIAO Ye-xiang, WANG Zheng-wei, ZHANG Jin, et al.Numerical predictions of pressure pulses in a Francis pump turbine with misaligned guide vanes[J].Journal of Hydrodynamics, 2014,26(2):250-256.

[7]QIAN Zhong-dong, YANG Jian-dong, HUAI Wen-xin.Numerical simulation and analysis of pressure pulsation in Francis Hydraulic Turbine with air admission[J].Journal of Hydrodynamics, 2007,19(4):467-472.

[8]桂中,唐澍,潘羅平.混流式水輪機尾水管非定常流動模擬及不規則壓力脈動預測[J].中國水利水電科學研究院學報, 2006, 4(1):68-73.

[9]楊建明.水輪機尾水管和轉輪中湍流計算研究[D].北京：清華大學,1999.

[10]李文鋒,馮建軍,羅興锜,等.基于動網格技術的混流式水輪機轉輪內部瞬態流動數值模擬[J].水力發電學報,2015,34(7):64-73.

[11]蘇華山,陳從平,趙美云,等.泵輪軸向振動條件下高速液力耦合器特性[J].2017,33(7):51-57.

[12]曾立飛,劉觀偉,毛靖儒,等.調節閥振動對閥內流場影響的數值模擬[J].中國電機工程學報, 2015,35(8):1 977-1 982.

[13]龐立軍,呂桂萍,鐘蘇,等.水輪機固定導葉的渦街模擬與振動分析[J].機械工程學報,2011,4(22):159-166.

[14]陳世云.大橋電廠水輪機過流部件汽蝕分析及處理[J].四川水力發電,2006,25(3):88-90.

[15]唐巍,張昌兵,李雷.中低比轉速混流式水輪機空化與水力振動分析[J].水電能源科學,2014,32(8):145-147.

[16]李雷,張昌兵,唐巍.混流式水輪機水力振動及補氣減振研究[J].中國農村水利水電,2015,(2):132-135.

[17]Shih T H, Lion W W, Shabbir A, et al.A newk-εeddy viscosity model for high Reynolds number turbulent flow[J].Computer Fluids, 1995,24(3):227-238.

[18]張師帥.CFD技術原理與應用[M].武漢:華中科技大學出版社,2016:84-85.

[19]王福軍.計算流體力學分析——CFD軟件原理與應用[M].北京:清華大學出版社,2004:114-141.

[20]符杰,曾永忠,宋文武,等.貫流混流式水輪機尾水管壓力脈動分析[J].設計與研究,2005,6(32):5-6,24.