水泵水輪機(jī)壓水充氣過(guò)程的非穩(wěn)態(tài)數(shù)值模擬湍流模型比較

張 新，趙俊龍，舒崚峰，陳順義，方 杰

（1.中國(guó)電建集團(tuán)華東勘測(cè)設(shè)計(jì)研究院有限公司，杭州 311122；2.中山大學(xué)中法核工程與技術(shù)學(xué)院，廣東珠海 519082）

0 引 言

抽水蓄能是目前電力系統(tǒng)中應(yīng)用最為廣泛、壽命周期最長(zhǎng)、容量最大、技術(shù)最成熟的一種儲(chǔ)能技術(shù)［1，2］，隨著風(fēng)電和太陽(yáng)能等隨機(jī)性間歇可再生能源裝機(jī)的快速增長(zhǎng)、核電開(kāi)發(fā)加快、超高壓遠(yuǎn)距離輸電、柔性直流電網(wǎng)等發(fā)展，抽水蓄能電站將對(duì)電網(wǎng)的儲(chǔ)能調(diào)節(jié)作用和安全運(yùn)行保障有更重要的作用［3，4］。水泵水輪機(jī)是抽水蓄能電站的核心裝置，要兼顧水泵抽水和水輪機(jī)發(fā)電兩種運(yùn)行工況。壓水氣系統(tǒng)主要用于電站的發(fā)電調(diào)相以及抽水調(diào)相工況，用于調(diào)控機(jī)組在抽水泵工況和發(fā)電水輪機(jī)工況之間的平穩(wěn)轉(zhuǎn)換［5］。如圖1所示，壓水氣系統(tǒng)的工作原理，是通過(guò)向水輪機(jī)轉(zhuǎn)輪室內(nèi)通入壓縮空氣的方法，將液面壓至轉(zhuǎn)輪以下，使轉(zhuǎn)輪在可壓縮空氣中旋轉(zhuǎn)以減小旋轉(zhuǎn)阻力和機(jī)械振動(dòng)，是保障抽水蓄能電站正常工況轉(zhuǎn)換的重要輔助系統(tǒng)［6］。

隨著水泵水輪機(jī)單機(jī)容量和額定水頭的提高，為了兼顧水泵工況和水輪機(jī)工況的空蝕要求，水泵水輪機(jī)裝機(jī)的淹沒(méi)深度越來(lái)越大，對(duì)壓水氣系統(tǒng)的供氣壓力提出了更高的要求［7-8］。高壓高速的非穩(wěn)態(tài)膨脹氣動(dòng)過(guò)程，常伴隨著劇烈的溫度變化、壓力沖擊、氣動(dòng)噪聲等現(xiàn)象，對(duì)壓水氣系統(tǒng)管路結(jié)構(gòu)的機(jī)械強(qiáng)度和熱疲勞強(qiáng)度也形成嚴(yán)峻考驗(yàn)。因此，研究壓水氣系統(tǒng)工作時(shí)的氣動(dòng)特性，對(duì)保障抽水蓄能電站的正常運(yùn)行、優(yōu)化水泵水輪機(jī)調(diào)控系統(tǒng)的設(shè)計(jì)方法，均有重要意義［9-12］。本文分別采用Spalart-Allmaras（SA）模型、Shear-Stress Transport κ-ω（SST）模型和Scale-Adaptive Simulation（SAS）模型，對(duì)壓水氣系統(tǒng)的壓水充氣工況進(jìn)行非穩(wěn)態(tài)模擬研究，對(duì)比不同湍流模型在壓水氣系統(tǒng)模擬問(wèn)題中的適用性，同時(shí)初步分析中壓充氣階段系統(tǒng)關(guān)鍵環(huán)節(jié)的流動(dòng)特征和熱力學(xué)變化。

1 湍流模型

N-S 方程是數(shù)值模擬求解流體問(wèn)題的基礎(chǔ)，但由于本構(gòu)方程中湍流應(yīng)力的求解難度極大，在實(shí)際工程應(yīng)用中通常采用時(shí)均化的方式進(jìn)行簡(jiǎn)化求解［13］。時(shí)均化是將瞬時(shí)湍流應(yīng)力分解為時(shí)均量和脈動(dòng)量、并對(duì)脈動(dòng)量進(jìn)行建模求解的方法，稱(chēng)為雷諾平均方法（Reynolds Averaging Navier-Stokes，RANS）。方程（1）是RANS方程的基本形式，其中針對(duì)項(xiàng)的各種時(shí)均化近似求解建模，稱(chēng)為RANS方法的湍流模型。

式中：t為時(shí)間；ρ為密度；u為速度；p為壓力；δ為Kronecker 符號(hào)；下標(biāo)i，j和k代表各分量，帶撇的上標(biāo)代表脈動(dòng)量。

針對(duì)可壓縮氣體充放氣問(wèn)題的模擬，SA 模型和SST κ-ω 模型是目前較為常用的湍流模型。方程（2）是SA 模型的控制方程［14］，以多組經(jīng)驗(yàn)化的代數(shù)參數(shù)計(jì)算，將應(yīng)力求解簡(jiǎn)化為關(guān)于中間變量ν~ 輸運(yùn)的單方程問(wèn)題［15］（單方程模型），降低了數(shù)值模擬的計(jì)算量。該方法的模型參數(shù)專(zhuān)門(mén)針對(duì)氣動(dòng)問(wèn)題設(shè)計(jì)，適合求解可壓縮氣體的高速流動(dòng)問(wèn)題，因此被廣泛用于燃?xì)夤苈沸孤兜雀邏悍艢膺^(guò)程的數(shù)值模擬。

式中：代表湍流運(yùn)動(dòng)黏性系數(shù)；Gν是黏性生成項(xiàng)；Yν是黏性耗散項(xiàng)。

方程（3）和方程（4）是SST 模型的控制方程，是典型的兩方程模型，通過(guò)設(shè)計(jì)混合函數(shù)F1，實(shí)現(xiàn)適用于壁面低速區(qū)的κ-ω模型與高速區(qū)的κ-ε模型的統(tǒng)一求解［16］。該模型對(duì)流速大范圍變化的問(wèn)題具有良好的適應(yīng)能力，適合氣罐排氣類(lèi)問(wèn)題的模擬研究。

式中：κ是湍動(dòng)能；ω是湍流頻率；Pκ是湍流生成速率；σκ3、σω3、α3、β3均為常數(shù)。

除了以上兩種模型，本文還考察了SAS 模型在高壓氣罐排氣問(wèn)題模擬中的表現(xiàn)。該模型是一種新近發(fā)展出的非穩(wěn)態(tài)RANS 方法，借鑒了大渦模擬（Large Eddy Simulation，LES）方法使用Kolmogorov 尺度進(jìn)行空間濾波的思想，引入馮卡門(mén)尺度識(shí)別流動(dòng)的穩(wěn)定和非穩(wěn)定區(qū)域進(jìn)行自適應(yīng)的精細(xì)化模擬。該方法除了對(duì)流速大范圍變化的問(wèn)題具有良好適應(yīng)性，還有較高的流動(dòng)結(jié)構(gòu)識(shí)別能力和噪聲分析能力，是進(jìn)行高壓排氣問(wèn)題精細(xì)化研究的潛在選擇，其具體控制方程，見(jiàn)參考文獻(xiàn)［17，18］.

2 模擬方案

壓水充氣工況中，壓水氣系統(tǒng)通過(guò)釋放儲(chǔ)氣罐中壓縮空氣的方式，自主地向水輪機(jī)轉(zhuǎn)輪室內(nèi)充氣并下壓水位，因此圖1所示充氣主管與補(bǔ)氣管的工作狀態(tài)相互獨(dú)立。本文對(duì)東南某抽水蓄能電站的壓水氣系統(tǒng)的壓水充氣工況進(jìn)行仿真，模擬所用流域建模如圖2所示。模擬采用全通道模型，保留了儲(chǔ)氣罐至轉(zhuǎn)輪室的全部管路結(jié)構(gòu)，并使用結(jié)構(gòu)化網(wǎng)格進(jìn)行空間離散，如圖3所示。

模擬過(guò)程與邊界條件設(shè)置與實(shí)際情況保持一致：壓水充氣工況中，補(bǔ)氣管路球閥關(guān)閉、充氣主管路球閥全開(kāi)；液壓球閥打開(kāi)后，初始?jí)毫? MPa的壓縮空氣，由儲(chǔ)氣罐沿管路排至轉(zhuǎn)輪室和尾水管中；忽略轉(zhuǎn)輪室和尾水管水位的下降過(guò)程，保留由尾水水位產(chǎn)生的0.9 MPa 恒定背壓條件；壓水氣罐與上部管路，裸露于室溫25 ℃的地下廠房?jī)?nèi)，應(yīng)考慮外壁面的空氣自然對(duì)流，傳熱系數(shù)約為50 W/m2·K［19］；下部管路、轉(zhuǎn)輪室和尾水管填埋于混凝土結(jié)構(gòu)內(nèi)，近似為絕熱條件。表1列出了模擬所用關(guān)鍵參數(shù)及設(shè)置。

表1 模擬參數(shù)設(shè)置列表Tab.1 The list of simulation parameters

3 方案與網(wǎng)格無(wú)關(guān)性驗(yàn)證

由于水電站混凝土澆筑建設(shè)的特點(diǎn)和安全運(yùn)行的要求，難以在建成和在建電站壓水氣系統(tǒng)內(nèi)增布測(cè)量設(shè)備，故缺乏相關(guān)電站壓水氣系統(tǒng)的運(yùn)行數(shù)據(jù)。為保證數(shù)值模擬結(jié)果的可信度，本文使用與壓水氣系統(tǒng)模擬相同的模擬方法，對(duì)文獻(xiàn)［20］中描述的高壓氣罐排氣研究（圖4）進(jìn)行數(shù)值模擬，并將模擬結(jié)果與文獻(xiàn)［20］的實(shí)驗(yàn)數(shù)據(jù)做對(duì)比。模擬使用模型及網(wǎng)格如圖5所示，模型總體積13.07 L，采用結(jié)構(gòu)化網(wǎng)格進(jìn)行空間離散，第一層網(wǎng)格高度0.01 mm，指數(shù)增長(zhǎng)率1.3，網(wǎng)格質(zhì)量0.7，計(jì)算時(shí)使用可伸縮壁面函數(shù)。邊界條件與文獻(xiàn)［20］一致，氣罐初始溫度25 ℃，初始?jí)毫?00 kPa，管路出口壓力1 atm，節(jié)流閥門(mén)通徑4 mm，氣罐與管路壁面為絕熱條件。

由于中高壓氣罐排氣過(guò)程中，氣流會(huì)在管路達(dá)到聲速并形成塞流，所以文獻(xiàn)中類(lèi)似的高速氣流問(wèn)題模擬，通常使用針對(duì)氣動(dòng)問(wèn)題優(yōu)化的SA 湍流模型［21］。本文亦使用SA 模型求解圖4模型，并分別使用了具有相同網(wǎng)格拓?fù)浜瓦吔鐚訁?shù)的多套網(wǎng)格進(jìn)行計(jì)算，網(wǎng)格單元數(shù)量分別為5×105、1×106、2×106、4×106。模擬結(jié)果與文獻(xiàn)［20］數(shù)據(jù)對(duì)比，見(jiàn)圖6。

圖6為放氣過(guò)程中罐內(nèi)壓力的變化曲線(xiàn)。曲線(xiàn)顯示，本文SA模型的模擬結(jié)果，與文獻(xiàn)［20］中絕熱條件的一維理論模擬結(jié)果和實(shí)驗(yàn)結(jié)果，在中高壓放氣階段吻合較好，但在低壓放氣階段與實(shí)驗(yàn)數(shù)據(jù)偏差較大。其原因正如文獻(xiàn)［20］所述：這是由于實(shí)驗(yàn)裝置非絕熱，在中高壓放氣階段氣體高速膨脹近似絕熱過(guò)程，因此與絕熱模擬結(jié)果一致；而低壓放氣時(shí)流速大幅下降，氣體膨脹轉(zhuǎn)為多變過(guò)程，而與絕熱模型模擬結(jié)果出現(xiàn)較大偏差。中高壓模擬結(jié)果與文獻(xiàn)數(shù)據(jù)對(duì)比的一致性表明，將本文設(shè)計(jì)的模擬方案用于模擬壓水氣系統(tǒng)的中高壓排氣工況，具有一定的合理性。

同時(shí)，圖6中不同網(wǎng)格數(shù)量的模擬結(jié)果高度重合，說(shuō)明針對(duì)13.07 L 高壓容器的放氣過(guò)程模擬，單元數(shù)量5×105的結(jié)構(gòu)化網(wǎng)格即可滿(mǎn)足網(wǎng)格無(wú)關(guān)性要求。因此，根據(jù)模擬問(wèn)題類(lèi)型的一致性，認(rèn)為壓水氣系統(tǒng)模擬在達(dá)到同等網(wǎng)格質(zhì)量、邊界層質(zhì)量和關(guān)鍵區(qū)域分辨率的情況下，可滿(mǎn)足網(wǎng)格無(wú)關(guān)性要求。最終壓水氣系統(tǒng)網(wǎng)格參數(shù)確定為：第一層網(wǎng)格高度0.01 mm，指數(shù)增長(zhǎng)率1.3，網(wǎng)格質(zhì)量0.5，高速氣流區(qū)域平均分辨率4 mm、低速氣流區(qū)域平均分辨率20 mm，網(wǎng)格單元總計(jì)6.05×106。

4 結(jié)果分析

本文對(duì)圖2所示壓水氣系統(tǒng)初始?jí)毫? MPa 的壓水充氣工況進(jìn)行數(shù)值模擬，模擬總時(shí)長(zhǎng)20 s。圖7為壓水充氣過(guò)程中，氣罐內(nèi)監(jiān)測(cè)點(diǎn)記錄的壓力、溫度和密度數(shù)據(jù)變化曲線(xiàn)。模擬結(jié)果顯示，初始?jí)毫? MPa的壓水氣系統(tǒng)，運(yùn)行20 s時(shí)氣罐壓力下降31.5%、溫度下降30 ℃，故該過(guò)程中儲(chǔ)氣罐運(yùn)行于中高壓排氣階段、罐體溫度下降顯著。在此過(guò)程中，不同位置測(cè)點(diǎn)的壓力和溫度完全相同，說(shuō)明氣罐內(nèi)部熱力學(xué)性質(zhì)均一；不同湍流模型模擬數(shù)據(jù)完全重合，說(shuō)明壓水氣系統(tǒng)壓水充氣工況的氣罐部分模擬結(jié)果，對(duì)湍流模型較不敏感。

將模擬記錄的壓力和溫度數(shù)據(jù)，按照剛性容器絕熱放氣過(guò)程理論方程（5）［22］，計(jì)算氣罐放氣過(guò)程中的密度變化率，并與模擬中氣罐中心測(cè)點(diǎn)的監(jiān)測(cè)值對(duì)比。如圖7所示，理論值與模擬值完全吻合，說(shuō)明在壓水氣系統(tǒng)壓水充氣工況中，儲(chǔ)氣罐前半段的排氣過(guò)程為絕熱過(guò)程。考慮到模擬中實(shí)際設(shè)置了氣罐壁面?zhèn)鳠嵯禂?shù)為50 W/m2·K的非絕熱條件，該模擬結(jié)果表明：壓水氣系統(tǒng)壓水充氣工況中，可以忽略氣罐外表面因自然對(duì)流形成的傳熱過(guò)程。

式中：k是絕熱指數(shù)；R是氣體常數(shù)。

圖8為壓水充氣管路入口（氣罐出口）徑向測(cè)點(diǎn)記錄的速度、壓力、溫度模擬結(jié)果。數(shù)據(jù)顯示，高壓氣體由氣罐進(jìn)入管路后迅速膨脹，空氣內(nèi)能轉(zhuǎn)化為動(dòng)能，壓力和溫度降低而速度大大提高。根據(jù)測(cè)點(diǎn)的位置關(guān)系，可知該階段壓水氣系統(tǒng)管路入口同一橫截面上壓力始終相等，速度和溫度均為對(duì)稱(chēng)分布；管路中心速度最高而溫度最低，且相對(duì)壁面附近氣體的速度差和溫度差近乎不隨時(shí)間發(fā)生變化，說(shuō)明此處膨脹為絕熱過(guò)程，動(dòng)能升高主要消耗系統(tǒng)內(nèi)能。

表2列出了本文模擬數(shù)據(jù)與寶泉電站壓水管路實(shí)測(cè)數(shù)據(jù)［23，24］的對(duì)比，在相近的氣罐壓降條件下兩者管路溫度基本一致，說(shuō)明本文各湍流模型的模擬結(jié)果具有合理性。結(jié)合圖7，可知SA、SST 和SAS 不同湍流模型計(jì)算得到的壓力和溫度標(biāo)量數(shù)據(jù)沒(méi)有區(qū)別；速度矢量數(shù)據(jù)在管路主流區(qū)內(nèi)結(jié)果相同，在壁面附近SA 模型與SAS 模型結(jié)果一致，相對(duì)SST 模型結(jié)果偏低約15%。考慮到該偏差值較小，且多數(shù)數(shù)據(jù)基本一致，可綜合認(rèn)為湍流模型在壓水氣系統(tǒng)管路入口的差異不足以干擾分析結(jié)果。

圖9為壓水氣系統(tǒng)尾水管空氣射流線(xiàn)上，測(cè)點(diǎn)記錄的速度和溫度模擬結(jié)果。數(shù)據(jù)顯示，壓縮空氣經(jīng)管路注入尾水管和轉(zhuǎn)輪室后，體積迅速膨脹，導(dǎo)致管路出口附近劇烈的速度和溫度波動(dòng)。根據(jù)測(cè)點(diǎn)的位置關(guān)系，可知隨著射流在尾水中的發(fā)展，射流氣體逐漸降速升溫，動(dòng)能重新轉(zhuǎn)換為內(nèi)能，符合等壓膨脹的一般特征。其中，高壓空氣到達(dá)管路出口時(shí)，溫度已降至0 ℃，并在20 s 內(nèi)降至-20 ℃，平均降溫速度高達(dá)-1 ℃/s，可能在管路與尾水管連接處產(chǎn)生較大的不均勻熱應(yīng)力，或產(chǎn)生局部結(jié)冰，威脅設(shè)備運(yùn)行安全。

尾水管內(nèi)的模擬結(jié)果顯示出，SA、SST和SAS不同湍流模型計(jì)算管路出口處射流時(shí)，所得溫度標(biāo)量和速度矢量沒(méi)有區(qū)別；但在射流深入尾水管而降速后，不同湍流模型對(duì)速度數(shù)據(jù)的計(jì)算，則出現(xiàn)較為明顯的差別。其中，SA 模型計(jì)算得到的射流降速幅度最大，SST 模型在射流中部計(jì)算結(jié)果與SAS 模型基本一致，在射流尾部計(jì)算結(jié)果較SAS模型略有偏低。考慮到SA模型是針對(duì)有壁面的氣動(dòng)問(wèn)題進(jìn)行設(shè)計(jì)，對(duì)無(wú)壁面的射流問(wèn)題模擬可能存在固有偏差，因此在開(kāi)展基于速度數(shù)據(jù)的壓水氣系統(tǒng)尾水部分的研究時(shí)，推薦采信SST模型和SAS模型的模擬數(shù)據(jù)；而在分析溫度等標(biāo)量數(shù)據(jù)時(shí)，SA、SST 和SAS 模型間的差異較小，不足以干擾分析結(jié)果。

通過(guò)可視化的流場(chǎng)數(shù)據(jù)分析，本文發(fā)現(xiàn)壓水氣系統(tǒng)中的節(jié)流孔板和管路出口，是壓水氣系統(tǒng)壓水充氣過(guò)程中，流動(dòng)變化最劇烈的兩處關(guān)鍵環(huán)節(jié)。圖10 顯示了節(jié)流閥在管路內(nèi)部制造的局部射流結(jié)構(gòu)。壓力云圖顯示出，該局部射流具有典型的膨脹波-壓縮波序列結(jié)構(gòu)，說(shuō)明射流為超聲速射流。溫度云圖則顯示出，射流激烈的膨脹加速過(guò)程大幅消耗氣體內(nèi)能，經(jīng)過(guò)節(jié)流閥的氣體溫度大幅下降，局部氣溫小于-50 ℃。超低溫氣流持續(xù)沖擊下流管路結(jié)構(gòu)（如止回閥），可能造成機(jī)械強(qiáng)度下降和結(jié)構(gòu)凍結(jié)等安全問(wèn)題。圖11則顯示了高壓氣體注入轉(zhuǎn)輪室-尾水管后形成的亞音速射流結(jié)構(gòu)，由于流速較慢、空間較大，該射流有充足的條件降速升溫，對(duì)尾水管壁面的沖擊較弱。

圖10 和圖11 亦對(duì)比了不同的湍流模型模擬壓水氣系統(tǒng)壓水充氣工況時(shí)，對(duì)相同流動(dòng)結(jié)構(gòu)的捕捉效果的差異。數(shù)據(jù)顯示，SST 模型和SAS 模型獲得的射流結(jié)構(gòu)幾乎相同，而SA 模型的射流則存在長(zhǎng)度偏短、膨脹-壓縮波結(jié)構(gòu)衰減偏快、亞聲速射流邊界形態(tài)模糊等問(wèn)題。圖12將不同時(shí)刻、不同湍流模型計(jì)算所得尾水管射流的等速度線(xiàn)圖像進(jìn)行對(duì)比，可知以上差異在整個(gè)非穩(wěn)態(tài)模擬過(guò)程中始終存在。因此相較于SA 模型，SST 模型和SAS 模型更適合壓水氣系統(tǒng)中高壓壓水充氣工況流動(dòng)發(fā)展過(guò)程的機(jī)理分析和研究。

5 結(jié) 論

本文使用SA 模型、SST 模型和SAS 模型，對(duì)壓水氣系統(tǒng)的壓水充氣工況進(jìn)行了全通道非穩(wěn)態(tài)數(shù)值模擬。結(jié)果表明，3 種模型均適用于壓水氣系統(tǒng)的宏觀現(xiàn)象分析和研究。其中SST和SAS模型對(duì)流動(dòng)結(jié)構(gòu)的識(shí)別能力優(yōu)于SA模型，更適合用于分析流動(dòng)機(jī)理。若考慮湍流模型計(jì)算效率的差異，SST 模型則可以更好的平衡模擬精度與模擬速度。根據(jù)模擬結(jié)果，本文發(fā)現(xiàn)壓水充氣工況的前半段整體符合絕熱過(guò)程的發(fā)展規(guī)律，但在節(jié)流閥和管路出口等局部位置，存在劇烈變化的射流結(jié)構(gòu)，可能導(dǎo)致壓水氣系統(tǒng)發(fā)生故障，有必要開(kāi)展深入研究。□