趙 煒，鄒曉陽，閔占奎，王文歡，王顥鈞

（1.國網(wǎng)甘肅省電力公司電力科學研究院，甘肅 蘭州 730070；2.上海電力大學能源與機械工程學院，上海 201306；3.國網(wǎng)甘肅綜合能源服務有限公司，甘肅 蘭州 730070）

0 引言

水電是可再生清潔能源，生產(chǎn)成本大幅低于煤電，因此發(fā)展水電不僅有助于環(huán)保、減少碳排放，而且能夠減少對煤炭等不可再生資源的消耗。我國水力資源豐富，主要位于西部地區(qū)，如黃河、長江、金沙江等，水電開發(fā)潛力巨大。然而部分河流泥沙含量大，容易導致推力軸瓦燒損和事故停機，造成很大經(jīng)濟損失。2018年7月，劉家峽水電廠新投運的7、8號機組——兩臺位于“穿黃排沙”洞出口處的水輪發(fā)電機組帶額定負荷運行，短時間內(nèi)推力軸承溫度迅速上升超過報警和停機溫度，造成推力軸承燒損事故[1]。推力軸承燒損的主要原因之一是軸向水推力大，導致彈性塑料推力瓦面油膜減少或消失，推力瓦面發(fā)生干摩擦，溫度急劇升高，超出推力瓦溫度上限[1]。此外，白巖水電站、龍羊峽水電站和尼泊爾崔樹里水電站的推力軸承燒瓦事故，也均由軸向水推力過大引起[2]。

近幾年，王偉 等通過模型試驗和現(xiàn)場測試的方法研究了水推力，分析了水推力過大的原因[3]。CFD仿真技術具有工作量小、時間短、成本低、參數(shù)化分析方便等優(yōu)點，已經(jīng)廣泛應用于流體流動模擬與分析。唐聰 等基于ANSYS CFX仿真分析了古里水電站混流式水輪機真機和模型上冠流域在不同泄水減壓結構下的流動狀況、相對泄漏量和軸向水推力[4]。李浩亮 等利用CFD計算了混流式蓄能機組的軸向水推力，并與現(xiàn)場測試的結果進行比較，得到了軸向水推力隨出力和水頭的變化規(guī)律[5]。JI 等人提出水輪機主流道CFD仿真結合上冠和下環(huán)理論計算的方法來計算水輪機軸向水推力，分析了軸向水推力隨水頭和流量的變化，與實驗測試結果相符[6]?；贑FD多相流對混流式水輪機含沙水流動進行仿真的研究工作也開展了不少。黃劍峰 等、廖姣 等和孫毅 等利用ANSYS FLUENT模擬了泥沙密度、泥沙直徑、泥沙濃度3個參數(shù)中的1個或幾個變化時含沙水在水輪機主流道內(nèi)的穩(wěn)態(tài)流動狀況，分析了泥沙對水輪機的磨損情況[7-9]。然而，水輪機全流道軸向水推力隨工況參數(shù)和含沙水參數(shù)的變化還有待進一步研究。

文中針對劉家峽水電廠的新投機組[1]，分別建立水輪機主流道、上冠流道和下環(huán)流道的CFD仿真模型，采用MIXTURE多相流模型和RNG k-ε兩方程模型模擬含沙水固液兩相流，基于FLUENT求取各流道的穩(wěn)態(tài)流場，通過流道合成計算獲得水輪機全流道的軸向水推力，分析其隨工況參數(shù)和含沙水參數(shù)的變化。

1 混流式水輪機含沙水流動仿真模型

1.1 混流式水輪機全流道模型及邊界條件

根據(jù)混流式水輪機軸向水推力的研究文獻[2],軸向水推力F包含3個主要分力，分別為流體作用在轉(zhuǎn)輪內(nèi)表面的軸向分力F1、作用在上冠外表面的軸向分力F2和作用在下環(huán)外表面的軸向分力F3，即：

其中，轉(zhuǎn)輪內(nèi)表面包括葉片外表面、上冠內(nèi)表面和下環(huán)內(nèi)表面。通過建立混流式水輪機的主流道、上冠流道和下環(huán)流道的模型進行仿真計算，合成為水輪機全流道的軸向水推力。

參考劉家峽水電廠的新投機組[1]，建立水輪機主流道、上冠流道和下環(huán)流道的三維幾何模型，如圖1所示。主流道比較復雜，沿流動方向依次為蝸殼、固定導葉、活動導葉、轉(zhuǎn)輪和尾水管，輪廓尺寸約為長33 m、寬17 m、高17 m，轉(zhuǎn)輪最大外徑約為5 m，流道入口直徑約為6 m，出口為尾水管出口。上冠流道的入口為寬度很窄的圓環(huán)面，無減壓泄水口，不考慮微小泄漏，上冠流道無出口。下環(huán)流道的入口為寬度很窄的圓環(huán)面，出口為高度很小的圓柱面。

利用ICEM CFD對水輪機流道進行網(wǎng)格劃分，流道內(nèi)部劃分為六面體網(wǎng)格，流道表面劃分為四邊形網(wǎng)格，對近壁面網(wǎng)格進行加密。通過interface將蝸殼、固定導葉、活動導葉、轉(zhuǎn)輪、尾水管的網(wǎng)格組裝起來，獲得主流道網(wǎng)格，如圖2所示。主流道網(wǎng)格總數(shù)約為710.3萬，其中轉(zhuǎn)輪網(wǎng)格數(shù)約為204萬。轉(zhuǎn)輪的葉片數(shù)為15，根據(jù)轉(zhuǎn)輪的旋轉(zhuǎn)周期結構，上冠流道和下環(huán)流道均在周向上取整體的1/15進行網(wǎng)格劃分，以減小網(wǎng)格模型。上冠流道網(wǎng)格數(shù)約為17.4萬，下環(huán)流道網(wǎng)格數(shù)約為31.5萬。

圖1 水輪機全流道三維模型

圖2 主流道網(wǎng)格

主流道入口采用速度入口邊界條件，在湍流流動下，速度分布采用純經(jīng)驗冪律公式描述，入口壓力將水頭和平均流速代入伯努利方程計算得出。主流道出口為壓力出口，出口壓力設為0。采用凍結轉(zhuǎn)子方法來模擬轉(zhuǎn)輪流域與其前后的靜止流域的界面作用，轉(zhuǎn)輪壁面跟隨轉(zhuǎn)輪流域一起轉(zhuǎn)動，其余壁面靜止，壁面與流域之間無滑移。上冠流道和下環(huán)流道入口均為壓力入口，入口壓力設置為主流道穩(wěn)態(tài)流動時相應位置的平均壓力。下環(huán)出口為壓力出口，出口壓力設為0。

1.2 含沙水流動控制方程和求解

含沙水屬于固液兩相流，主相為水，次相為泥沙。文中主要關注含沙水對水輪機軸向水推力的影響，基于ANSYS FLUENT 2019平臺對含沙水流動進行仿真時，選擇MIXTURE多相流模型。MIXTURE多相流模型屬于歐拉-歐拉方法，將固體顆粒視為與連續(xù)相相互滲透的擬流體進行處理，計算量適中，計算穩(wěn)定性好，工程應用廣泛。

MIXTURE模型的連續(xù)性方程、動量方程、次相的體積分數(shù)方程分別為

其中t為時間，ρm是含沙水的平均密度，vm是質(zhì)量平均的速度矢量，p為壓強，g為重力加速度向量，F(xiàn)為體積力向量，μm為體積平均的動力粘度，n為多相流的相數(shù)，αs為第s相的體積分數(shù)，vdr,s為次相s的漂移速度向量，次相體積分數(shù)方程等號右邊第二項求和項為相間質(zhì)量轉(zhuǎn)化速率[10]，▽為向量微分算子。MIXTURE模型利用代數(shù)滑移方程計算次相相對于主相的滑移速度，進而由滑移速度計算次相的漂移速度[10]。在ANSYS FLUENT平臺上仿真時泥沙的動力粘度設置為微小量，相間的拖曳力采用Schiller-Naumann模型，相間滑移速度采用Manninen-et-al模型。

利用雷諾時均方法（Reynolds-averaged Navier-Stokes, RANS）求解含沙水固液兩相流的控制方程。雷諾時均方法將瞬時速度和瞬時壓力等效為短時間內(nèi)的平均值與脈動值之和，代入到湍流流動控制方程中，因湍流脈動出現(xiàn)的新應力項，稱為雷諾應力。目前工程上應用最廣泛的湍流模型為兩方程模型，即將雷諾應力表示為湍動粘度的函數(shù)，湍動粘度通過湍動能k和湍動耗散率ε的方程求解。標準k-ε模型適用于計算充分發(fā)展的湍流，計算穩(wěn)定性好，對計算資源要求不高，適合于較大雷諾數(shù)、低旋和弱浮力的湍流模擬。RNG k-ε模型與標準k-ε模型相似，對標準k-ε模型進行了修正，除了適合于較大雷諾數(shù)、低旋和弱浮力的湍流模擬，也能夠模擬高應變率、大曲率和低雷諾數(shù)湍流[10]。在FLUENT 2019平臺上采用RNG k-ε兩方程模型和標準壁面函數(shù)對含沙水流動進行模擬，模型參數(shù)值采用FLUENT中的默認值。

求解含沙水兩相流的控制方程時，采用壓力-速度耦合算法。壓力離散采用體積力加權法，動量離散采用二階迎風格式，體積分數(shù)、湍動能和耗散率離散均采用QUICK法。松弛因子采用比FLUENT的默認值小的值，以保證數(shù)值求解時具有良好的收斂性。

2 軸向水推力的仿真結果

2.1 額定工況下轉(zhuǎn)輪內(nèi)外流道的壁面壓力

對額定工況下水輪機全流道含沙水流動進行穩(wěn)態(tài)仿真。額定工況參數(shù)為：水頭93.6 m，流量175.4 m3/s，轉(zhuǎn)速136.4 r/min。含沙水參數(shù)為：水密度1 000 kg/m3，泥沙密度2 710 kg/m3，泥沙顆粒Φ 0.1 mm[11，12]。

圖3和圖4分別為額定工況下清水（泥沙含量0）和泥沙含量30 kg/m3的含沙水流動時轉(zhuǎn)輪、上冠和下環(huán)流道壁面壓力?？梢钥吹剑D(zhuǎn)輪壁面壓力呈圓周對稱分布，高壓區(qū)位于葉片的入流區(qū)域，低壓區(qū)位于葉片的出流區(qū)域及靠近出流邊的背壓面。葉片的工作面和背壓面有一定的壓力差，葉片表面的靜壓大體上順著流線方向遞減。上冠流域的壓力較大，壓力隨著半徑增大而增大，壓力梯度也增大。下環(huán)流域最大壓力位于入口處，且沿著間隙流道，壓力呈均勻下降趨勢。進入下環(huán)空腔，流通截面大，壓力由正變負。下環(huán)底部外側直角處出現(xiàn)了漩渦，流動較為復雜。泥沙含量30 kg/m3的含沙水流動時的壓力分布狀況與清水流動時的基本相同，轉(zhuǎn)輪壁面和上冠壁面的壓力略微增加，而下環(huán)壁面壓力略微下降。額定工況下，清水穩(wěn)態(tài)流動時的水輪機軸向水推力為9 425.8 kN，方向向下。而30 kg/m3的含沙水流動時的軸向水推力清水時的軸向水推力之比約為1.006。

2.2 軸向水推力隨工況參數(shù)的變化

以額定工況下清水穩(wěn)態(tài)流動時的水輪機軸向水推力為比較基準，定義軸向力比，分析軸向水推力隨轉(zhuǎn)速、流量、水頭等工況參數(shù)的變化。仿真時含沙水參數(shù)保持不變。

圖3 額定工況清水流動水輪機流道壁面壓力

圖4 額定工況泥沙含量30 kg/m3的含沙水流動水輪機流道壁面壓力

圖5 為水頭為93.6 m時不同轉(zhuǎn)速下水輪機軸向水推力比隨流量的變化曲線，圖6為流量175.4 m3/s時不同轉(zhuǎn)速下水輪機軸向水推力比隨水頭的變化曲線?？梢钥吹?，含沙水參數(shù)不變時，在0.85~1.15倍的額定轉(zhuǎn)速范圍內(nèi)，轉(zhuǎn)速對水輪機軸向水推力的影響不大，除了流量95.4 m3/s處軸向力比值的最大差值約為10%之外，其他流量下軸向力比值的最大差值在5%以內(nèi)，圖6很清楚地表明了這一點。水頭對軸向水推力的影響更小，圖6所示流量175.4 m3/s時相同轉(zhuǎn)速不同水頭下軸向力比值之差在1%以內(nèi)。流量對軸向水推力的影響最大，如圖5所示，流量從95.4 m3/s增加到175.4 m3/s時，軸向力比值從約0.3近似線性增加到1.0左右。

圖5 不同轉(zhuǎn)速下水輪機軸向水推力比隨流量的變化

圖6 不同轉(zhuǎn)速下水輪機軸向水推力比隨水頭的變化

2.3 軸向水推力隨含沙水參數(shù)的變化

保持工況參數(shù)與額定工況參數(shù)相同，仿真分析軸向水推力隨泥沙密度、泥沙直徑、泥沙含量等含沙水參數(shù)的變化。

圖7為泥沙含量30 kg/m3時不同泥沙密度下水輪機軸向水推力比隨泥沙顆粒直徑的變化曲線，圖8為泥沙顆粒直徑0.1 mm時不同泥沙密度下水輪機軸向水推力比隨泥沙含量的變化曲線?？梢钥吹剑嗌澈亢苄。ㄈ?0 kg/m3）時，泥沙密度在1 710~2 710 kg/m3的范圍內(nèi)變化時，對水輪機軸向水推力的影響較小。泥沙顆粒直徑對軸向水推力的影響也很小，圖7所示泥沙含量30 kg/m3時相同泥沙密度不同顆粒直徑的軸向力比之差在1%以內(nèi)。其他參數(shù)保持不變時，軸向水推力增大隨泥沙密度的增大而增大。泥沙含量越大，不同泥沙密度下的水推力比值相差越大，如圖8中泥沙含量為210 kg/m3時，2 710 kg/m3密度下的軸向水推力比與1 710 kg/m3密度下的增大5%。泥沙含量對軸向水推力的影響很大。泥沙含量從30 kg/m3增加到210 kg/m3時，軸向力比值近似線性增加，且密度越大，增加的幅度越大，密度為1 710 kg/m3時增加約7%，密度為2 710 kg/m3時增加約12%。

圖7 不同泥沙密度下水輪機軸向水推力比隨泥沙顆粒直徑的變化

圖8 不同泥沙密度下水輪機軸向水推力比隨泥沙含量的變化

根據(jù)上述的仿真結果和分析，可以看到含沙水在水輪機中流動時，對水輪機的軸向水推力會產(chǎn)生較為復雜的影響。與清水相比，含沙水增大了水輪機軸向水推力，尤其在泥沙含量很大時，將會大幅度地增加軸向水推力，這是造成過流泥沙含量大的水輪發(fā)電機組推力軸瓦燒損的原因。在工況參數(shù)中，流量對軸向水推力的影響最為顯著，減小流量可以明顯地減小水推力。文中的研究對于減小混流式水輪機軸向水推力以避免推力軸瓦燒損具有重要的指導意義。含沙水在水輪機中的流動及其引起推力軸瓦燒損是非常復雜的，今后將考慮CFD瞬態(tài)仿真分析和試驗測試等深入研究水輪機的軸向水推力。

3 結論

文中建立水輪機主流道、上冠流道和下環(huán)流道的CFD仿真模型，在ANSYS FLUNET中采用MIXTURE多相流模型和RNG k-ε兩方程模型對水輪機含沙水的湍流流動進行穩(wěn)態(tài)計算，獲得了水輪機全流道的穩(wěn)態(tài)流場，分析了水輪機全流道水推力隨工況參數(shù)和含沙水參數(shù)的變化，獲得的主要結果如下：

（1）額定工況下，清水流動和小泥沙含量（如30 kg/m3）的含沙水流動下的水輪機轉(zhuǎn)輪內(nèi)外流道壁面的壓力大小和分布接近，含沙水造成水輪機軸向水推力增大。

（2）含沙水參數(shù)不變時，轉(zhuǎn)速在額定轉(zhuǎn)速附近變化或水頭在上下限范圍內(nèi)變化時，水輪機軸向水推力的變化很小。流量從95.4~175.4 m3/s時，軸向水推力比值從約0.3近似線性增加到1.0左右。

（3）工況參數(shù)不變時，泥沙顆粒直徑在0.01~0.2 mm范圍內(nèi)變化時，軸向水推力變化很小。泥沙密度增大時，軸向水推力增大，且泥沙含量越大，增大幅度越大。泥沙含量從30~210 kg/m3時，軸向水推力近似線性增加，泥沙密度為1 710 kg/m3時增加約7%，泥沙密度為2 710 kg/m3時增加約12%。