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

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

0 引言

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

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

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

1 混流式水輪機(jī)含沙水流動仿真模型

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

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

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

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

利用ICEM CFD對水輪機(jī)流道進(jìn)行網(wǎng)格劃分，流道內(nèi)部劃分為六面體網(wǎng)格，流道表面劃分為四邊形網(wǎng)格，對近壁面網(wǎng)格進(jìn)行加密。通過interface將蝸殼、固定導(dǎo)葉、活動導(dǎo)葉、轉(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)周期結(jié)構(gòu)，上冠流道和下環(huán)流道均在周向上取整體的1/15進(jìn)行網(wǎng)格劃分，以減小網(wǎng)格模型。上冠流道網(wǎng)格數(shù)約為17.4萬，下環(huán)流道網(wǎng)格數(shù)約為31.5萬。

圖1 水輪機(jī)全流道三維模型

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

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

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

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

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

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

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

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

2 軸向水推力的仿真結(jié)果

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3 結(jié)論

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

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

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

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