水泵水轮机压水充气过程的非稳态数值模拟湍流模型比较

张 新，赵俊龙，舒崚峰，陈顺义，方 杰

（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 模型则可以更好的平衡模拟精度与模拟速度。根据模拟结果，本文发现压水充气工况的前半段整体符合绝热过程的发展规律，但在节流阀和管路出口等局部位置，存在剧烈变化的射流结构，可能导致压水气系统发生故障，有必要开展深入研究。□