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

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

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

0 引 言

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

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

1 湍流模型

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

式中：t為時間；ρ為密度；u為速度；p為壓力；δ為Kronecker 符號；下標i，j和k代表各分量，帶撇的上標代表脈動量。

針對可壓縮氣體充放氣問題的模擬，SA 模型和SST κ-ω 模型是目前較為常用的湍流模型。方程（2）是SA 模型的控制方程［14］，以多組經(jīng)驗化的代數(shù)參數(shù)計算，將應力求解簡化為關(guān)于中間變量ν~ 輸運的單方程問題［15］（單方程模型），降低了數(shù)值模擬的計算量。該方法的模型參數(shù)專門針對氣動問題設計，適合求解可壓縮氣體的高速流動問題，因此被廣泛用于燃氣管路泄露等高壓放氣過程的數(shù)值模擬。

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

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

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

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

2 模擬方案

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

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

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

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

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

由于中高壓氣罐排氣過程中，氣流會在管路達到聲速并形成塞流，所以文獻中類似的高速氣流問題模擬，通常使用針對氣動問題優(yōu)化的SA 湍流模型［21］。本文亦使用SA 模型求解圖4模型，并分別使用了具有相同網(wǎng)格拓撲和邊界層參數(shù)的多套網(wǎng)格進行計算，網(wǎng)格單元數(shù)量分別為5×105、1×106、2×106、4×106。模擬結(jié)果與文獻［20］數(shù)據(jù)對比，見圖6。

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

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

4 結(jié)果分析

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

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

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

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

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

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

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

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

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

5 結(jié) 論

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