水泵水輪機壓水充氣過程的非穩態數值模擬湍流模型比較

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

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

0 引 言

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

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

1 湍流模型

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

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

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

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

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

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

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

2 模擬方案

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

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

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

3 方案與網格無關性驗證

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

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

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

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

4 結果分析

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

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

式中：k是絕熱指數；R是氣體常數。

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

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

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

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

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

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

5 結 論

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