有压管道系统非恒定流模拟研究进展

王云莉 ，胡 宏

（1.重庆交通大学 西南水运工程科学研究所，重庆 400016；2.内河航道整治技术交通行业重点实验室，重庆 400016）

科研与管理

有压管道系统非恒定流模拟研究进展

王云莉1，2，胡 宏1，2

（1.重庆交通大学 西南水运工程科学研究所，重庆 400016；2.内河航道整治技术交通行业重点实验室，重庆 400016）

介绍了国内外有压管道系统非恒定流数值模拟研究的发展现状，现阶段有压管道系统非恒定流水力特性数值模拟中存在主要问题，提出了解决这些问题的主要方法及数值模拟的发展趋势，可为类似研究提供参考。

数值模拟；有压管道系统；非恒定流；研究进展

1 发展现状

20世纪70年代，计算机的计算速度与效率大大提高，数值模拟方法模拟水击成为了现实。数值模拟主要的优势在于对边界条件的模拟，在水击与调压室涌浪的联合计算中有明显优势，解析法和图解法慢慢被淘汰。一维数值模拟中运用最广最多的是电算法，是以特征线法为基础发展起来的。但随着科学技术的不断发展，人们越来越重视对流场内部的研究，水击的研究方向也逐渐由一维向三维模拟转变，三维数值模拟是目前水击的发展方向，三维数值模拟与成熟的特征线解法比较起来，更能准确地模拟出流场结构与阀门各种特性在关闭过程中的变化。因此，三维数值模拟也逐渐成为研究有压管道系统非恒定流水动力特性的主要模拟方法［1-3］。

2 有压管道系统的模拟

在水利枢纽等有压管道系统的设计和运用中，水击是要考虑的重要因素之一，而为了减小水击的危害，常需要在有压管道中设置调压室。当电站处于某一固定负荷时，机组的引用流量固定，调压室内的水位为某一固定水位。随着电站负荷的改变，机组引用流量也相应发生改变，此时调压室内发生质量波动，波动的稳定是调压室研究中的重要内容。若调压室水位过低，空气从调压室进入压力钢管形成含气水锤；若调压室水位过高，水流则从调压室溢出，两者都可能导致重大的安全事故，所以工程上对调压室的最低、最高涌浪水位非常重视［4-6］。

此外，在初期水电站设计中，一般是把调压室涌浪和压力管道水击分开进行计算，忽略其相互影响。伴随着计算机技术的发展与成熟，可以将水击与调压井进行联合计算。数值模拟是进行有压管道非恒定流计算的常用方法。目前，有压管道系统数值模拟的方法有很多，主要分为一维数值模拟和三维数值模拟。

2.1 一维数值模拟

2.1.1 控制方程

2.1.1.1 水击基本方程

水击基本方程主要有：运动方程、连续方程。以管道进口为原点，指向下游为正，则有：

运动方程：

运动方程是从管道水体中选取隔离体，应用牛顿第二定律建立的。

连续方程：

考虑到管道中水体的可压缩性及管壁的弹性，对非恒定流连续方程做适当变换而得到连续方程。式中 g为重力加速度（m/s2）；H为压力水头（m）；v为流速（m/s）；x为距离（m）；t为时间（s）；f为摩阻系数；D为管道内径（m）；θ为管道轴线与水平面夹角（°）；a为水击波的传播速度（m/s）。

2.1.1.2 调压井涌浪方程

运动方程：

连续方程：

式中 z为库水位与调压室水位之差；hf为引水隧洞总的水头损失；L为引水隧洞总的长度；f为引水隧洞断面面积；F为调压室断面面积；v为引水隧洞中的流速；Q为引水隧洞中的流量；Qt为流入或流出调压室的流量；Q0为一台水轮机的引用流量；m为机组的数量。

2.1.2 水击和调压室涌浪联合运行的一维数学模型

联合运行的数学模型中，通常以水击的运动方程和连续方程作为基本方程，将调压室的方程组作为边界条件引入计算。通过方程的离散、网格划分、边界条件分析、时间步长的设定等进行计算。将方程组离散、求全倒数、整理得特征线方程：

2.2 三维数值模拟

2.2.1 控制方程

2.2.1.1 连续方程

连续方程是任何数值计算必须遵守的方程，即质量守恒方程：

式中 V为控制体；A为控制面，首项为流入控制体内质量增量，第二项为流出控制体质量差。其微分表达式为：

对于不可压缩流体，ρ为常数。

2.2.1.2 N-S方程

N-S方程是大多流体流动必须遵守的方程，即为Navier-Stokes方程，也称为动量守恒方程。

式中 fx，fy，fz分别为单位流体质量力在不同方向上的分量；P为流体内应力的张量；sx，sy，sz是流体内的广义源项。对于黏性的可压缩流体，s，f为零。N-S方程则比较系统地描述了流体的实际流动形态，所有的流动问题的数值解法都是围绕该方程组进行求解的。

2.2.2 方程求解

随着人们对自然界中紊流运动从整体不断向微观发展，数值分析在求解明渠水流流动问题的精细化水流模拟要求也越来越高。随着研究者不断深入的探讨研究与发展，有限差分法、有限单元法、有限体积法、边界单元法等越来越成熟，在实际过程中发挥了重要作用。

2.2.2.1 有限差分法

从实际问题引出的微分方程出发，将连续的定解区域离散为合理分布的网格，把微分方程及定解条件中的微商用相应的差商代替，对得到的差分代数方程组进行求解，国外学者已从数学上严格证明了差分方程收敛性相容性、稳定性及其定解条件，其结果就是原方程在定解区域上的近似解。有限差分法有较成熟的理论基础，具体算法简单易行。

2.2.2.2 有限单元法

将流动问题计算域剖分成各个小区域，对这些小区域运用不同的插值函数选择区域基函数，将区域结果整合至整个计算域。运用有限差分和变分原理，对不同形状的单元选择逼近的插值函数进行积分处理，适用于边界复杂的问题。

2.2.2.3 有限体积法

使单元中每个节点周围都包含一个控制体，假设网格节点上的因变量在网格节点之间的分布规律，将控制方程对每个控制体体积分，组成一系列离散方程，再结合边界条件进行求解。牛顿守恒定律自然表现形式，适用于任何类型网格单元。有限体积法可以认为是在有限元法基础上改进的有限差分法，结合了两者思想，使计算形式更加统一。目前在很多领域，该方法都得到了广泛应用。

2.2.2.4 边界元法

它不需要对整个问题域进行划分，只需要将问题域的边界划分为一系列的单元，边界微分方程划分为积分方程，再在边界上转化为代数方程而进行求解。该方法可以降低求解问题的维数，对于明渠流自由表面的处理较简单，且精度较高。

2.2.3 紊流数值模拟理论

2.2.3.1 紊流数值模拟

对于雷诺数高于2300的流动称为紊流，自然界中常见的水流问题多数为紊流，一般的流体力学往往不能有效地解决紊流问题，而常常结合一定的紊流模型才能达到求解问题的精度。

2.2.3.2 紊流数值模拟方法

紊流运动是非常复杂的非线性流动，数值模拟常需要很高的精度才能达到流体运动的实际状态，经过近些年研究学者的努力，对紊流模型的研究取得了非常显著成就。研究方法分类如图1。

图1 紊流模型分类形式

由于本研究主要应用Reynolds平均法（RANS），现对其进行详细介绍。

（1）Reynolds平均法（RANS）。紊流运动是非常复杂的非线性流动，对于直接求解N-S方程很难达到其精确的解析解，研究者往往研究紊流的瞬时性，所以常把奈流的瞬时值看作是时间平均值和脉动值的叠加。Reynolds平均法在研究三维紊流问题时具有很高的精度，并得到了广泛应用。在Reynolds平均法中任一瞬时变量的时间平均值定义为：

在Reynolds平均法中用时均值和脉动值代替瞬时值后，可得到Reynolds时均的紊流运动控制方程。

连续方程：

N-S方程：

由于时均方程中出现了雷诺应力项uiuj导致方程组不封闭，故引入涡黏性假设。

式中 μt为涡黏性系数；ui为i时刻的时均速度；k为紊动动能。

根据对Reynolds应力作出的假设或处理方式不同，目前常用的模型有两大类：Reynolds应力模型和涡粘模型。Reynolds应力模型方法中，直接构建表示Reynolds应力的方程，然后联立求解及建立新的Reynolds应力方程。通常情况下Reynolds应力方程是微分形式，称为Reynolds应力方程模型。若将Reynolds应力方程的微分形式简化为代数应力方程，则称为代数应力方程模型［12］。

（2）涡粘模型。主要分为：零方程模型、一方程模型和二方程模型。其中：

零方程模型：最出名的是混掺长度模型推导出μt，混合长度lm和时均速度u关系式，则：

一方程模型：为了弥补混合长度的局限性，建议在湍流时均方程和Reynolds方程的基础上，建立μt，l和k的关系式，则：

二方程模型：随着研究者对于求解精度的需要，转而寻求通用性更高、方法更为精细的紊流模型，更有效地确定特征长度分布。

2.2.4 自由水面的模拟

在水利水电工程、船舶、海洋工程中，水流的自由表面流的模拟具有着重大的现实应用价值。目前，常采用欧拉法和拉格朗日法处理自由液面流动问题。欧拉模型主要体现为在固定的坐标系下或坐标系以一定的方式移动来适应自由表面的移动，而拉格朗日模型则是坐标会随着流体域的变化而变化，该方法能精确地定位自由面的形状变化，但是算法极不稳定。

随着自由表面模拟的不断深入发展，越来越多的技术被应用到实际工程中来，目前应用最多的有对称性边界条件、标高函数法、MAC法、体积分辨率法等。

2.2.4.1 对称性边界

对称性边界条件俗称“刚盖假定”，是最早应用于自由水面模拟的技术，对称性边界条件即假设通过此边界的所有物理梯度量为零，假设自由液面边界不会随时间变动而变动，该方法通常应用于自由表面比较规则且变化幅度不是太大的流动问题。当处理泄洪、水跃、溃坝水流、波浪等水流问题时，显然该法的精度不能满足问题需要。

2.2.4.2 标高函数法

假设水深是关于x，y，z的单值函数，以此建立自由表面的运动方程：

方程组达到封闭，即可求解表面任一处的水深h。但该方法也有一定局限性，对于复杂的自由表面流问题，即该法很难有效解决水深A不是关于x，y，t的单值函数的水流问题。

2.2.4.3 MAC法

MAC法的显著特点就是能有效解决水深的多值坐标函数问题。主要表现为在计算域空间的网格内设置如同质点的“标记点”，这些点只是用来记录追踪流体的运动轨迹，不参与计算，根据流场的计算，可以得到每个标记点的速度，以此确定被标记的网格单元的状态，从而确定出下一时刻自由面的位置。

2.2.4.4 体积分率法（VOF）

体积分率法又称为体积追踪法，是将流体在一个网格单元中所占的体积比率当作指标，可以求解空气与水两相介质的交互界面，计算体积比率的改变来追踪分离界面的运动，利用体积分率来标记某计算网格为水、空气、自由液面，再利用几何观念求解此流体体积方程式。釆用线性方法，根据网格内水与空气的比例来代表自由液面的斜率，表现出网格内流体间的界面，进而决定自由液面的形状。如何决定控制体积内空气与水物理性质所占的比例，则需通过体积分率来决定。表达形式为：

式中 δΩcell为控制体积的大小；δΩwater为控制体内水的体积。当F=0时，表示控制体内充满空气；当F=1时，表示控制体内充满水；当0＜F＜1时，表示控制体内同时存在水与空气，也就是自由液面的位置。

2.2.5 三维数学模型

水击和调压室涌浪联合运行三维数学模型中，运用FLUENT软件进行计算。边界条件主要包括：固壁边界条件、进出口边界条件。

对于本模型而言，给出的具体边界条件如图2。

图2 有压管道系统三维数学模型示意图

进口边界类型：进口1为速度进口；出口边界类型：出口1为压力出口，采用VOF技术处理水/空气自由面波动问题；出口2为OUTFLOW。当出口2有闸门时，要模拟闸门的动态启闭和有压管道系统的模拟时，就要引入UDF函数和动网格技术。对各种闸门的启闭过程进行动态仿真时，动网格用于调整和更新闸门运动后的网格，以确保计算网格质量；UDF用于定义阀门的开关规律，使阀门实现连续的启闭过程。为了便于计算，本研究中闸门是线性启闭。

3 结语

（1）有压管道系统非恒定流一维数值模拟已相对成熟，其具有代表性的特征线法有计算速度快，效率高，对于复杂性管道的适应性很强的优点，流行了一段时间。但是伴随特征线法的广泛使用，其多维模拟能力较差的弊端也越来越明显，无法得到有压管道系统流场分布和压力分布，而三维数值模拟刚好可以满足这一点，有压管道系统非恒定流三维数值模拟将是今后的发展方向。

（2）有压管道系统非恒定流三维数值模拟中涉及到空气和水流两种介质不断交换问题，可采用成熟的VOF模型进行处理。VOF模型是为解决多相流问题而发展起来的流体体积模型。其原理是根据网格单元内流体体积的变化来追踪自由面的变化。近年来，VOF模型在计算两相流中得到应用。技术已趋于成熟。喻宇［14］运用VOF理论成果模拟出有压管道水流冲击气团的局部三维流场。

（3）国内外的研究现状来看，已有调压室与水击联合进行的一维数值模拟和有压管道的局部三维数值模拟，但尚无一套完整理论进行调压室与水击联合进行的三维数值模拟，这是一个巨大的创新点，建议开展此研究，为类似工程提供理论参考，为有压管道系统数值模拟提供研究方向。

［1］段昌国.水电站引水管道水击的解析与计算［J］.科学通报，1980（3）：134-137.

［2］石宝珍，王学芳，叶宏开.长输管线水击计算方法的研究［J］.清华大学学报（自然科学版），1986（1）.

［3］张晓宏，金钟元.碧口电站压力引水系统非恒定流的分析与计算［J］.陕西水力发电，1992（4）.

［4］张师华.水击压力和调压室涌浪联合计算的程序编制和分析［J］.武汉水利电力学院学报，1985（1）：23-26.

［5］张健，索丽生.阻抗式调压室甩负荷涌浪计算显式公式［J］.水力发电学报，1999（2）：70-74.

［6］彭天玫，杨效亮.阻抗式调压井穿井压力计算及最优阻抗孔口的研究［J］.水力发电学报，1986（3）：61-77.

［7］是勋刚.湍流直接数值模拟的进展与前景［J］.水动力学研究与进展（A辑），2001（4）：21-27.

［8］罗坤，刘小云，樊建人，等.三维混合层的直接数值模拟和实验验证［J］.中国电机工程学报，2005（6）：13-18.

［9］吴腾虎.管道内颗粒悬浮流中大颗粒和湍流相互作用的直接数值模拟研究［D］.杭州：浙江大学，2011.

［10］李斌，颂平.基于湍流尺度的混合RANS/LES模型［J］.北京航空航天大学学报，2008（1）：31-35.

［11］刘小兵，曾永忠.混流式水轮机中旋涡流的LES法预测［J］.大电机技术，2009（2）：22-27.

［12］马雷.边界层理论在低比转速离心泵叶片设计中的应用［J］.兰州理工大学学报，2009（2）：41-45.

［13］严继松，廖国玲.有压管道充水过程水力特性三维数值模拟［J］.水利水电技术，2015， 11（8）： 219-224.

［14］喻宇.有压管道水流冲击气团的局部三维流场数值模拟研究［D］.哈尔滨：哈尔滨工业大学，2011.

Research progress in numerical simulation method of unsteady flow in pressure pipeline system

WANG Yun-li1，2，HU-Hong1，2
（1.Key Laboratory of Inland Waterway Regulation Engineering， Ministry of Transport， PRC，Chongqing 400016，China；2.Southwest Hydraulic Institute for Waterways，Chongqing Jiaotong University，Chongqing 400016，China）

This paper introduces the development of numerical simulation method of unsteady flow in pressure pipeline system at home and abroad，puts forward the existing problems in the simulation of pressure pipeline system at the present stage， and analyzes the possible methods to solve these problems， which can provide a reference for similar research.

numerical simulation；pressure pipe system；unsteady flow；research progress

TV131.4 文献标识码：B 文章编号：1672-9900（2017）06-0003-05

2017-09-15

内河航道整治技术交通行业重点实验室开放基金项目（NHHD-201512）

王云莉（1974-），女（汉族），云南永善人，副研究员，主要从事港口及航道专业研究，（Tel）13271829317。

王艳肖）