2.2　有压管道输水过程模拟

2.2.1　复杂输水工程的常态有压管道输水过程

常态有压管道输水为恒定输水工况，水流恒定，因此工况可以直接采用初始条件的方法进行仿真计算。首先由上游端水位开始往下游端迭代至末端面，采用达西-威斯巴赫公式计算沿程水头损失，而局部水头损失计算则按控制建筑物特性分别计算。正常输水时的稳定输水情况，此时上、下游端水位都近似保持恒定，其输水流量除去分水外沿线保持恒定，而水头则应为阶梯状直线。输水管线恒定流下沿程水头线及流量如图2.2-1和图2.2-2所示。

对于任意段管道微段，断面i与i+1之间由恒定总流方程有

其中

式中:H为测压管水头。

暂不考虑局部水头损失，管道断面管径不变时(变管径时分段计算)，按均匀流达西-威斯巴赫公式计算沿程水头损失有

对于泵站系统，恒定流主要分析供水泵站在水泵的不同组合和进、出水池水位不断变化的条件下水泵的工作点，以及相对应的泵站系统的各种水力要素，以确定水泵的稳定工作范围，泵站的流量、扬程和效率等，尽量使水泵保持在高效区运行，以实现节能降耗。

2.2.1.1　泵基本特性曲线

在一定的转速下，根据试验测算出水泵多个流量Q对应的扬程H、效率η、功率N和汽蚀余量NPSH值，绘出水泵的基本特性曲线即Q-H曲线、Q-η曲线、Q-N曲线和Q-NPSH曲线。由Q-H曲线可以看出水泵的压头在较大流量范围内随流量增大而减小。Q-η曲线表示水泵的效率，在该曲线上有一最高点，该最高点被称为设计点，无论流量过小或过大，效率都会降低。一般高效区为不低于最高效率的92%的范围。由Q-N曲线可知泵的功率一般随流量增大而增大。

2.2.1.2　管路特性曲线

泵站需要扬程为

式中:H净为净扬程，即上、下水面间的垂直高度，m;H损为通过管路液体单位液重的水力损失水头，m;hf为管路沿程损失水头，m;hj为局部损失水头，m;A为管路截面面积，m2;f为管路摩阻系数;∑ξ为各局部阻力系数之和;K为管路特性系数，对已知管路系统为一常数。

因为一般H净=const，所以曲线Q-H需是在纵坐标为H净的一条水平线叠加上对应流量Q时管路损失水头H损，所得曲线为管路特性曲线。

2.2.1.3　单泵单管泵运行工作点的确定

泵站的需要扬程曲线Q-H需与水泵的流量扬程曲线Q-H的交点就是水泵运行的工作点。

泵站的需要扬程等于水泵净扬程与管道阻力损失之和，若管道的阻力参数为S，则泵站的需要扬程为

水泵的特性曲线可以根据生产厂家提供的实测数据通过最小二乘法拟合而得到，即

能量平衡方程为

由式(2.2-7)可求出流量，将求出的流量代入式(2.2-6)可求得扬程。同理可拟合出效率曲线，求出水泵的工作效率。

2.2.1.4　复杂泵系统运行工作点的确定

在实际供水工程中，往往需要将多台同型号水泵或者不同型号水泵并联运行。

2.2.1.4.1　同型号水泵并联运行时水泵工作点的确定

当n台同型号水泵并联运行时，相当于水泵扬程H不变，流量扩大n倍，计算原理同单泵单管。此时泵站的需要扬程为

结合水泵扬程式(2.2-6)，可得能量平衡方程为

将求出的流量代入式(2.2-9)可求得扬程。

泵站流量:

泵站扬程等于水泵扬程。泵站效率为

式中:η泵为水泵的效率;η传为传动装置的效率;η电为电动机的效率;η管为管道的效率;η池为进、出水池的效率。

2.2.1.4.2　不同型号水泵并联运行时水泵工作点的确定

当不同型号水泵并联运行时，分别对不同型号水泵建立能量方程，其方法和单泵单管相似，得到的能量平衡方程为

式中:H1为大泵扬程;H2为小泵扬程;Q1为大泵流量;Q2为小泵流量;n1为大泵台数;n2为小泵台数;S1为大泵进出水管损失系数;S2为小泵进出水管损失系数;S为并联管路损失系数。

联立式(2.2-12)和式(2.2-13)求解可得Q1和Q2。

然后利用拟合得到的Q1-H1和Q2-H2曲线，分别求得大小泵的扬程H1和H2。同理可得到大小泵的效率η1和η2。

泵站流量:

泵站扬程:

泵站效率:

式中:m为并联运行水泵的台数;i为各台水泵的序号(i=1，2，3，…，m)。

2.2.2　有压管道输水水力过渡过程

解决水资源短缺问题最直接、最有效的方式就是修建调水工程，把水从水资源丰富的地区输送到水资源短缺的地区。这是解决水资源地区分布不均与供需矛盾的一项有效措施，也是促进缺水地区经济发展与水资源综合开发利用的重要途径之一。

因为工程地质条件和地形条件的不同，调水方式和结构也呈现多样化。既有渠道、管道、隧洞、天然河道及其相互连接的多种复杂结构，又有自流、加压提水以及自流与加压提水相结合的调水方式，在加压提水的调水系统中，还有大流量、低扬程泵站或小流量、高扬程泵站等形式。由此在调水工程中衍生出了许多的技术难题，管道输水系统在运行调节过程中发生的水力过渡过程就是其中之一。所谓的管道输水系统的水力过渡过程就是指由于流量的调节或事故的发生等原因，系统的边界条件发生变化，原有的恒定流条件遭到破坏，从而在系统内引发非恒定流的现象。此时系统沿线的压力、水位、流量、流速等水力参数随时间发生急剧的变化，往往会超出恒定流的设计范围，对输水过程造成不利影响，如有压管道的水击、管道负压、水体脱空、奎水溢流、水流振荡等现象。

因此，在工程的设计和优化运行当中必须对管道输水系统的水力过渡过程进行计算，确定在各种情况下引发水力过渡过程时沿线的水力要素的变化情况，从而为工程的设计和管理运行提供依据和建议。随着大量输水工程的兴建，人们对输水系统发生的过渡过程越来越重视，对输水系统水力过渡过程控制的研究已成为极具实际意义的课题。

2.2.2.1　瞬变流基本理论

瞬变流是指压力管道中的流体因某种原因而产生流速的骤然变化时，在流体的惯性作用下而导致管内流体的压力急剧变化的现象。在水力学中，瞬变流也称为水击或水锤。水击可分为刚性水击和弹性水击。在刚性水击中，水体和压力管道的管壁被作为刚体考虑。在弹性水击中，主要考虑水体的惯性和可压缩性，压力管道的管壁被看作是弹性的。通常由于弹性水击更能反映水击的实际情况，因此相关研究大多针对弹性水击进行阐述。

瞬变流的基本微分方程主要有运动方程和连续方程。其理论基础是水流运动的力学规律和连续性原理。它是全面反映有压管流中非恒定流流动规律的数学表达式，包含水流流速和水头的变化规律。瞬变流基本微分方程是研究水力过渡过程的基础。

2.2.2.1.1　瞬变流运动方程

流体的运动方程也称为动量方程，它可以用动量守恒的基本原理进行推导。为了推导瞬变流的运动方程，在发生状态改变的水流中取如图2.2-3所示的微小控制体。控制体的长度为d x，其断面面积为A，管道的倾角为α。假设x坐标的方向与原来的水流方向一致，上、下游断面的水压力为P1、P2，控制体周边界面上的阻力为P3，侧水压力为P4，重力为mg。若上游m-m断面的密度为ρ，过水断面面积为A，湿周为X，压强为P，则下游n-n断面相应各量分别为图2.2-3中θ为控制体侧壁与管线夹角，一般取cosθ=1。假设τ为控制体周边平均阻力，沿x轴方向的流速为v，则有压管道水体的加速度

根据动量定理，作用在控制体上的力的总和应当等于水体质量乘以加速度，即

因为v是时间t和坐标x的函数，所以可得

将式(2.2-18)代入式(2.2-17)，并取整理可得

因为测压管水头线f为恒定流时的阻尼系数，所以式(2.2-19)可简化为

式(2.2-20)即为有压管道瞬变流的运动偏微分方程。

2.2.2.1.2　瞬变流连续方程

对于图2.2-4所示的连续方程控制体，设上游m-m断面的流体密度为ρ，面积为A，沿x轴方向的流速为v，根据质量守恒定律可以叙述为:流入、流出控制体的水体质量之差等于控制体水体质量随时间的变化，即

整理并简化可得

将式(2.2-22)展开得

因为所以式(2.2-23)可简化为

此方程由质量守恒定律得来，适用于任何形式的流体，适用于有压流也适用于明渠水流。

管道中的水击波速公式为

式中:K为水的体积弹性模量;E为管壁材料的纵向弹性模量;g为重力加速度;D为管道内径;δ为管壁厚度。

如果忽略了管道轴向变形的影响，那么

式中:εT为管道环向应变;σT为应力。

又所以式(2.2-25)可以简化为

根据流体的体积弹性模量的定义有下式成立:

又由压力与测压管的关系P=(H-Z)γ可得

联立式(2.2-27)～式(2.2-29)可得

整理后可得

式(2.2-31)与式(2.2-24)联立并整理可得

式(2.2-32)就是分析有压管道水力过渡过程的连续方程，式中的第二、四项与第一、三项相比较小，应用时往往忽略。

2.2.2.2　瞬变流计算数学模型

管道输水工程过渡过程计算涉及的数学模型有管道瞬变流计算模型、前池边界计算模型、阀门边界计算模型、泵边界计算模型和空气阀边界计算模型。

2.2.2.2.1　管道瞬变流计算方程

对于管道中的瞬变流(非定常流动)，其计算微分方程组如下。

运动方程:

连续方程:

式中:H为沿程水头;v为平均流速;g为重力加速度;f为达西-威斯巴赫摩擦系数;α为管道中心线与水平线的夹角;D为管道直径;c为波速。

2.2.2.2.2　方程求解

水锤问题计算方法通常有解析法、图解法、特征线方法和隐格式差分方法。在求解运用实际中，特征线方法理论严谨，原理简单，求解水锤问题精度高。因此，模型选用特征线方法进行求解水锤微分方程组。

根据特征线理论，式(2.2-33)和式(2.2-34)是一对准线性双曲线型偏微分方程，存在两条特征线，可以采用特征线方法进行求解。

现对这两个偏微分方程分别用L1和L2标示:

令某实数λ乘以L1再与L2相加，其和仍应等于零:

若采用两个不同的λ值代入式(2.2-37)，则可得到包含有未知数H和v的两个新方程式，这两个新方程式是由式(2.2-35)和式(2.2-36)组成的，因此可以等价替换原有的两个方程式，选择两个特殊的λ值，目的是使新的方程组能改造为常微分方程的形式。

根据复合函数的微分运算规则:

分别令

当整个式(2.2-37)为全微分方程时，将代入式(2.2-39)和式(2.2-40)可得

将得到的两个值先后代入式(2.2-37)，可得到与式(2.2-35)和式(2.2-36)等价的两个常微分方程组为

以上就是管内流态暂态的特征线方程，特征线求解示意如图2.2-5所示。图中A点和B点代表地点x和时间t已给定的两个点，它们的H和v是已知的，曲线C+和C-称为特征线，式(2.2-42)和式(2.2-43)称为相容性方程，由于原始的基本微分方程式并没有加进任何运算上的近似，故相容性方程的解，也就是原始基本微分方程式的解。从特征线可以看出，从初始边或边界边出发，沿C+和C-两个方向可以唯一确定管道任意位置x时刻t的未知量和值，如图2.2-5所示。

大多数工程管道为刚性壁管材料制成，波速c比v大得多，式(2.2-42)和式(2.2-43)中的流速项可以忽略不计。这时的特征线C+和C-都将成直线形状，在x-t坐标图上，这两条直线的斜率分别为和式(2.2-43)简化为

由于仿真模拟往往只需要确定特定位置、特定时间或者说是要确定离散位置、离散时间的输水管线压力水头与流量大小，因此可以直接沿特征线方向离散微分方程进行差分求解。如图2.2-6所示，将管路划分为多个Δx的步段，断面排列号用i表示，计算时段P点就可以通过离散点A与B进行求解。式(2.2-44)和式(2.2-45)中流速v用流量Q代替。特征线差分计算网格如图2.2-6所示。

其中

式中:B、R分别为计算常数。

若将各参数和计算时段初始已知的参数组合在一起，则式(2.2-46)和式(2.2-47)可紧缩改写为

C+:C-:

其中

综合参数CP和CM，可在各计算时段开始时先算出。式(2.2-48)和式(2.2-49)化简为

式(2.2-52)和式(2.2-53)便是适宜于编入计算机程序的相容性方程，应当指出的是:特征线的交点，只包括图2.2-6中的矩形网格内结点，两端断面上的参数，则还必须通过各个瞬时的边界控制条件才能确定。各边界控制条件包括阀门边界计算模型、泵边界计算模型、空气阀边界计算模型、前池边界计算模型、管道连接点边界方程和初值计算。

2.2.2.3　边界条件的计算方程

单一管路的两个端部结点都只有一个相容性方程可以利用。对于上游端如图2.2-7(a)所示，只存在特征线C-，可利用的相容性方程为式(2.2-52);对于下游端如图2.2-7(b)所示，只存在特征线C+，可利用的相容性方程为式(2.2-53)。

为了要确定边界上的两个控制参数HP和QP，必须再补充一个边界条件方程。常见的边界条件方程有以下几种条件。

(1)边界上的HP或QP是独立于管道控制系统的控制参数，当管路上、下游为水位恒定的水池时，边界结点的HP是固定常数，相容性方程可用来求解QP。

(2)边界上的HP和QP之间存在一定的函数关系，当边界上有正常运转中的水泵时，泵的性能曲线规定了HP=f(QP)，它与相容性方程联立可解出边界上的HP和QP。

(3)边界上的HP和QP还与其他的边界参数有关，如发生事故停电时泵的性能曲线与其瞬时转速有关系。由于增加了瞬时泵转速这个新变量，就需要再多增加一个边界条件方程。

2.2.2.3.1　前池边界方程

对于大面积的泵站前池边界，在短促的暂态期间，水位可认为是恒定的。设前池的水位为HR，则管路上游结点1的水头HP1=HR。

当水池的面积F较小时，水池水位在暂态过程中也会发生变化，其结点水头为

式中:HR为计算时段开始时的水池水位;Q1为计算时段开始时的流量值。

HP1确定后，可采用相容性方程式(2.2-55)计算结点流量:

2.2.2.3.2　阀门边界方程

一个阀门装在已知管线内或两根不同的管线之间，阀门的孔口方程要和管段的端部条件联合处理，而且应当考虑流动反向的可能性。忽略通过阀门的流体加速或减速的影响，即假定存在于阀门中的流体体积不发生变化，则瞬变过程中仍然可以使用定常态的阀门孔口方程:

或

式中:ζ为阀门的阻力系数;C为阀门开启面积乘以流量系数。

ζ和C两者关系为

式(2.2-56)或式(2.2-57)和阀门上、下游侧的相容性方程式(2.2-52)和式(2.2-53)联立，可解得3个未知数:流量QP和阀门上、下游两侧的两个HP值。

式中:HP1、QP1分别为阀门前截面下一时刻的压头和流量;HP2、QP2分别为阀门后截面下一时刻的压头和流量。

对于在管线中的阀门，阀门前方管道中的参数用下标1表示，阀门后方管道中的参数用下标2表示。

2.2.2.3.3　泵边界方程

对泵进行流动暂态分析，首先要对叶片泵全面性能曲线进行改造，以便计算机在计算过程中取值。在这个问题解决后，计算机才能对泵处的两个边界条件方程——水头平衡方程和转速变化方程进行运算。

(1)叶片泵全面性能曲线的改造。离心泵有4个特性参数:扬程H、流量Q、轴转矩M和转速N。这4个量中只有两个是独立的。若以Q和N作为纵横两坐标参数，则全面性能曲线包含一组等H值曲线和另一组等M值曲线。

利用全面性能曲线进行水锤计算时，采用了以下两个简化的假定作为前提。

1)假定在恒定(稳态)流动条件下实测所得的曲线，可以反映暂态过程中各参数之间的关系。

2)相似关系成立。

现用一组无因次参数表示工作参数与额定参数的相对比值:

式(2.2-60)中带角标r的值表示额定值，将同一台泵在不同转速条件下的性能曲线相互叠合在一起，转化为无因次参数新坐标的性能曲线:

采用式(2.2-61)中的新坐标，可将原有全面性能曲线所包含的两组等值曲线(等H值和等M值曲线组)改绘为两条无因次曲线，从而大大便利于给参数赋值供计算调用。

但是由于都是可正可负，当然也可能出现为零值的情况，故各坐标参数的变化幅度均为-∞～+∞，以致实际上不可能绘出整个曲线。为此，麦切尔(Marchal)等提出了巧妙的方法，将式(2.2-60)中的无因次坐标参数再改造为以下的无因次新坐标参数:

式(2.2-62)中横坐标x以弧度计，用下式表示:

式(2.2-62)中纵坐标用下式表示:

横坐标的变化幅度为0～2π，纵坐标的变化幅度也有限，这就可在有限的坐标系统中绘出整个曲线。

WH(x)和WM(x)的曲线形状比较复杂，很难用简单的数学公式对它有比较满意的近似描述。WH(x)和WM(x)用系列离散数据表示，从x=0至x=2π，以等分间距取下89个WH(x)和89个WM(x)离散值，按x的次序排列在表中。

计算过程中出现的值不一定正好等于89个x的离散值之一(即x不一定正好等于Δx的整倍数)。因此，还需要通过线性内插，才能确定WH(x)和WM(x)值。内插方法如下。

令I表示的整数部分，即

实际x值位于(I-1)Δx和IΔx之间，从离散数据表中可以查到x=(I-1)Δx时所对应的WH(I)和WM(I)值，x=IΔx所对应的WH(I+1)和WM(I+1)值。用通过结点参数值的直线近似地代表两结点间的微段曲线，则该微段曲线可以近似用下列直线公式替代:

其中

对于每个计算时段，泵处的边界条件方程，都要先依据v和β值算出x值，并通过上述方法内插，计算相应的WH(x)和WM(x)值。

(2)水头平衡方程:

式中:HS为泵在吸水管一侧的测管水头，m;HP为泵在出水管一侧的测管水头，m;ΔHpump为泵的扬程，m;ΔHF为阀门引起的水头损失，m。

根据在叶片泵全面性能曲线中的讨论，可得到泵扬程为

由式(2.2-56)，阀门水头损失可写为

其中

式中:ΔH0为阀门全开时的水头损失，即设计工况下的水头损失。

在阀门水头损失式(2.2-72)中，将v2写成v|v|是为了反映水头损失的正负性质:当正向流动时，v>0，ΔHF>0;当反向流动时，v<0，ΔHF<0。

吸水管末段和出水管首段的特征线方程可按照式(2.2-48)和式(2.2-49)改写为

式中:CP、BS为吸水管中的已知参数;CM、BP为出水管中的已知参数，分别按式(2.2-50)和式(2.2-51)计算。

将式(2.2-71)～式(2.2-74)一齐代入水头平衡关系式(2.2-70)，可得

式(2.2-75)中流量为(单泵单管系统)

为了简化起见，将已知参数项归并为

则形式较为紧凑的水头平衡方程可改写为

由于式(2.2-78)中Hr、A0、A1、Δh0都为常数，亦为计算时段开始瞬间的已知参数，在密云水库调蓄工程中，为两泵并联工作向一条出水总管供水，则式(2.2-75)中的QP=2v Qr，又因为吸水管很短，故可忽略其中的水锤效应，BS=0，此时。

为解得v和β这两个未知数，除用水头平衡方程式(2.2-70)外，还必须建立泵的另一个边界条件方程——转速变化方程式。

(3)转速变化方程。泵的转速变化是由所加的不平衡转矩引起的。当启停泵发生时，由轴输入的转矩为零，作用在叶轮上的转矩及其他摩擦力矩将使转速发生变化，故有转矩方程:

其中

式中:GD2为水泵机组回转部分的飞轮力矩;M为反力矩，可取每一计算时段初始瞬刻的力矩M0(已知)和结束瞬刻的力矩M(待求值)的算术平均值;ω为泵轮的旋转角速度。

改写式(2.2-70)，得

根据全面性能曲线的结论可得

故转速变化方程的最后形式可表述为

式(2.2-84)中B0、B1、C0、m0、β0都是常数或为已知参数，只有β、v是独立的未知数。将式(2.2-79)和式(2.2-80)联立，可解得计算时段结束瞬刻的泵处状态参数β和v。

整理泵处的两个边界条件方程式(2.2-78)和式(2.2-84)，即

通过式(2.2-85)可解出β和v。但由于式(2.2-85)为非线性的超越方程，只能采用迭代方法求解其近似值。牛顿-莱福生(Newton-Raphson)迭代公式形式为

2.2.2.3.4　空气阀边界方程

空气阀按其低于大气压系统进气，高于大气压系统排气的物理特点建立模型，用差分法求解。

对空气阀作如下假定。

(1)空气等熵地流入流出阀门。

(2)管内的空气温度始终不变。

(3)流进管内的空气仍留在空气进排气阀附近。

由以上假定，要求空气阀安装在管线顶点并被选为计算截面，如图2.2-8所示。

由图2.2-8，对于进入(或排出)管道的空气部分，根据气体状态方程有

式中:P为断面处的绝对压强;V为空穴部分体积;M为空气质量;M0为空气摩尔质量，标准状况下为29g/mol;R为摩尔气体常数，一般为8.314J/(mol·K);T为热力学温度。

流过阀门的空气质量主要影响因素有管外大气的绝对温度T0、绝对压力P0和管内的绝对压力P、绝对温度T。分以下4种情况进行考虑。

空气以临界流速流进:

式中:M为空气质量;C1为进气时阀的流量系数;A1为进气时阀的流通面积;R为气体常数;P0为绝对压力;T0为管外大气的绝对温度。

空气以亚音速流进:

式中:ρ0为大气密度;P为管内的绝对压力;其他符号意义同前。

空气以临界流速流出:

式中:A2、C2分别为排气时阀的流通面积和流量系数;其他符号意义同前。

空气以亚音速流出:

式中:符号意义同前。

根据式(2.2-91)可知需要补充关于空穴体积V与质量M的表达式，根据流体连续性，对于微段，其空穴体积的增加等于流体体积的减少，即

式中:V、V0分别为当前时刻和前一时刻的空穴体积;Δt为仿真时段步长。

另外有

式中:ρAir为空气密度，标准状况下约为1.29kg/m3，一般情况下可以采用此值近似估算，在气穴时需要依据压强P与体积V来确定。

同时根据测压管水头定义，断面出水头HP满足

式中:Pa为标准大气压;z为管道底部高程;D为管道直径。

根据实际进排气情况，将式(2.2-88)～式(2.2-93)代入式(2.2-94)得

另外有

联立式(2.2-95)和式(2.2-96)以及上、下游连接管道特征线方程可以实现空气阀的仿真计算。需要说明的是ρAir和P与V有关，可以采用前一时段的压强P0和气穴体积V0推求并进行迭代计算。

2.2.2.3.5　管道连接点边界方程

管道连接(包括串联、并联、分叉和汇合)的边界条件应满足连续方程和能量方程。以如图2.2-9所示的情况为例，管道1、2串联。

由流量连续有

根据能量方程有

式(2.2-98)左边第一、二项分别表示相应断面的测压管水头和速度水头，第三项表示水头损失，K为局部损失系数，绝对值保证了流动反向时也成立;右边第一、二项分别表示相应断面的测压管水头和速度水头。对多管的分叉和汇合所形成的边界条件，同样可以按连续方程和能量方程来给出。

2.2.2.3.6　初值计算

虚拟阻抗流体网络计算法计算简单，计算时间短，可以计算复杂网络结构系统的初值，本书中对于输水工程的过渡过程计算，其初值用虚拟阻抗流体网络计算法，即采用流体网络计算方法对系统进行初值计算，将管道、水泵看作阻抗。

流体网络可以看作是由一系列的单元组成，许多单元之间以一定数量的结点相连。对如图2.2-10所示管元，规定水从k流向j时管道流量为正。设分别为单元i连接于结点k、j的结点水头，分别为单元i连接于结点k、j的结点流量，规定从结点流出流量为正，由能量守恒定律可得

式(2.2-99)中的Ki通过损失系数求得，令管道损失系数为Si，则有

式(2.2-99)用矩阵形式表示为

简写为

单元结点流量矢量为

单元特征矩阵为

单元结点水头矢量为

在管网中任一结点上需满足流量连续方程:连接此结点的所有单元从结点流出流量之和等于输入该结点的流量(若是消耗则为负值)，对任一结点有

式(2.2-103)中的“∑”为对结点m有贡献的所有N单元求和。由各结点的流量平衡方程可将各管元的能量平衡方程合并成总体方程组，简写为

式中:为管网水头矢量，由管网各结点的水头组成;为管网的输入矢量，由管网各结点处的输入流量(若为损耗则为负)组成;为管网特性矩阵，由管元特性矩阵扩充而成。

每个管元有两个连接点，由管元特性方程可知，在管网特性矩阵中每个管元有4个基值，对如图2.2-10所示的管元，管元i只对j、k两点有贡献，对其他结点的贡献为零，将管元的系数Ki叠加到管网特性矩阵的(k，k)、(j，j)位置上，将系数-Ki叠加到(k，j)、(j，k)位置上，考虑所有的管元即形成管网特性矩阵。由矩阵形成过程可知，矩阵为带状对称稀疏矩阵，解方程组时可利用该特性来节省计算机内存和计算时间。

对于管网总体方程组，必须引进适当的边界条件才能求解，对任一结点m有两种可能的边界条件，或者规定结点输入(或消耗)Cm，或者规定水头值Hm，在管网计算时，为了求解方程组，结点边界条件必须至少规定有一个水头已知。由于管网特性方程组中的系数Ki与流量有关，计算中先给定管道的流量，通过迭代计算，使得前后两次的流量计算值相差在给定误差范围内，来求得非线性方程组的解。

实际的泵站输水系统中，除了有纯管道外，还可能会有泵、阀门、调压井、空气阀、空气罐等非管道元件，需将其看作特殊管网来计算。调压井、空气阀、空气罐这些元件在系统稳定状态下不起作用，在系统初值计算时可以先不考虑其作用，系统各管道的初值计算出后，用其所在管道结点的初值计算出这些元件的相应初值，以此作为系统动态仿真计算的初值。对于水泵、阀门，考虑到在实际管网计算中，每个管道的作用相当于阻抗，所以水泵、阀门的作用也可虚拟成阻抗，只不过阀门的作用相当于正阻抗，而水泵的作用相当于负阻抗。

采用该方法，在不改变程序的条件下，对于泵站水系统的不同布置方案均可以灵活计算其初值，为泵站的过渡过程计算实现灵活组态打下了坚实的基础。

2.2.3　管道充、排水过渡过程

有压输水方式输水响应快，是现代输水工程中最为常用的输水方式。但是由于有压管道水击波速大，加之输水系统的输水线路长，有压输水系统在过渡过程中引发的水力波动也较为剧烈，管理运行的控制难度较大。特别是有压管道的充、排水操作工况，其管理控制更为复杂。首先输水管线长，管道的容水体积巨大，充、排水时间较长;再者，充、排水引发的过渡过程涉及的水流形式多，流态复杂。这些因素都影响着输水系统在充、排水操作时的快速、安全、稳定的充、排水目标的实现，因此对输水系统充、排水工况的研究具有十分重要的意义。本章对管道充、排水过程所涉及的水流特征进行了总结介绍，对输水管道一维非恒定流过程的基本控制方程及其求解方法进行了阐述，同时对模拟气液二相瞬变流最为常用的两种数学模型——自由气体-液体离散模型和变波速模型进行了介绍，从而为管道输水系统充、排水水力过渡问题模拟研究提供一定的理论依据和借鉴。

2.2.3.1　输水管道充、排水过程的水流基本特征

就输水管道中任意一个单独的管道而言，管道中的水流并不是恒定不变，而是一直随时间不规律地变化的，因此管道中的水流状态为非恒定非均匀流，它主要有以下两个特点。

(1)不同时刻流动状态的复杂性。因为受到管道进口、出口变化的入流和出流的影响，管道中的流动也会呈现出不同的状态，速度上可能是急流或者缓流，压力上可能是明流、满流或者明满流，变化频率上可能是急变流或者渐变流，相性上可能是单相流，也可能因为掺入了空气而形成水气两相流，在管道的水力条件剧烈变化时甚至有可能受到上、下游边界条件的影响而形成反向流。

(2)不同流态之间转换的复杂性。因为管道中水流状态的不固定性，管道中可能存在着不同状态的水流之间的转换，如缓流和急流的转换、明流和满流的转换等。

缓流和急流之间的转换主要表现为缓流到急流之间的水跌或是急流到缓流之间的水跃。在某一特定的管道中，受到上、下游边界条件的影响，有可能上游是急流、下游是缓流，也可能正好相反，并且上、下游条件是不断变化的，因此急流和缓流之间的过渡面也是不断变化的，这也给计算增加了难度。

明流和满流之间的转换，分界面的形式主要包括正界面和负界面两种。在不同的限制条件下，有可能是明流过渡为满流，也有可能是满流过渡为明流。另外，明流和满流之间的过渡有可能在管道的入口或者出口产生，也有可能在管道中间的某处产生，产生的具体位置由管道中的水气条件决定。

2.2.3.2　气液两相流模拟技术

有压管路水击波的传播速度是分析管道水击问题的一个十分重要和敏感的参数，它对水击压力计算值的影响很大，水击波速的大小以及变化过程将直接影响有压过渡流的压力波动过程的幅值和周期。压力管道的水击波速计算公式为

式中:K为水的体积弹性模量;ρ为水的密度;D为管道直径;E为管材的弹性模量;e为管壁厚度。

可见，水击波速的大小涉及流体和管道的许多影响因素，其值部分取决于水的体积弹性模量，部分取决于管道的膨胀能力。如果输水管道中不含有自由空气或是蒸汽，水击波速就保持常数。但是如果水中存在有被截留的气体团，或是夹杂有自由气泡，以及由于压力降低而形成的空气或是蒸汽泡，那么水的有效体积弹性模量就会发生很大变化，水击波速将变成压力和时间的函数，加之水中含有的气泡或气体本身所具有的不确定性，所以二相流的瞬变过程是十分复杂的。

自由气体在液体中的含量通常用空穴比表示，其中Vg为自由气体的体积，V为含有气泡的液体的体积。

长期以来，对于气液两相瞬变流的模拟主要采用两种数学模型:一种是自由气体-液体离散模型，该方法将压力管道中的气泡或气团限制在固定的截面上，而认为各气泡截面之间不含有气体，仍然采用固定的水击波速对瞬变流进行计算，通过这种假设，来等效地模拟水击波速变化的情况;另一种是变波速模型，该方法假定液体中的气泡均匀地分布在整个管段，在水力过渡流的计算过程中需要根据水击压力以及空穴率的变化来计算输水系统不同截面的水击波速变化，并采用变化的水击波速来进行水力过渡过程的计算。这两种方法虽然假设条件和处理方式各不相同，但是都在一定的条件下反映出了水中含有的气体与水力过渡过程之间的相互联系和相互影响，气泡的含量和体积影响着水力过渡过程的变化，反过来水力要素的变化又引起了水中气泡含量和体积发生变化，因此这两种方法在处理二相流问题上得到了广泛的应用。本章通过这两种数学模型对气液二相流的水力过渡过程进行了模拟，对气液二相瞬变流的水力特性进行了研究，并对这两种模型进行了比较，且对模型的适用性进行了探讨。

2.2.3.2.1　自由气体-液体离散模型

自由气体-液体离散模型不把注意力集中在两相流的气泡动力学上，而是从宏观着眼来处理含有气团或自由气泡的水力过渡过程。该模型把气团或气泡限制在相应的固定计算断面上，气团或气泡的体积随水击压力的变化规律符合等温的完善气体状态方程，而认为压力管道的断面与断面之间不存在气体，仍然采用固定水击波速进行计算。通过该假设，可以等效地模拟气泡均匀分布的气液二相流的水力瞬变，当然也适用于管道中的高位部分有气团聚结或是只有管系的一部分存在气团或自由气泡的情况。

根据液体中存在气体的百分比，将计算断面j和j+1之间的管涵中的自由气体都集中在断面j+1上。根据等温的完善气体状态方程，有以下关系:

式中:V0，j+1为固定在j+1断面的气团在大气压下的体积;Pa为大气压力;Pg为箱涵中的气团压力;Vj+1为压力Pg下的气团体积。

气团压力与测压管水头的关系为

式中:γa为气体容重;HP为测压管水头;Z为该截面的管顶高程;Ha为气压计压头。

采用如图2.2-11所示的矩形网格，其特征线的相容性方程为

C+:

C-:

含有气团的断面处的连续方程可以表示成如下关系式:

式中:QP、QPU分别为t+d t时刻，即待求时刻计算结点的来流和出流流量;Q、QU分别为t时刻，即已知时刻计算结点的来流和出流流量;VgP、Vg分别为t+d t和t时刻计算结点处气团的体积;φ为加权因子，将其引入是为了控制数值求解过程中的数值振荡，0.5≤φ≤1。

将式(2.2-106)～式(2.2-110)联立即可求得含有气泡的有压管涵的水力过渡过程。

自由气体-液体离散模型还可以处理水力过渡过程中发生的液柱分离及再弥合现象。在水力过渡过程中，输水管道中的某点的压力降低至水的汽化压力，当这种压力降低持续相当长的时间时，管道的水流将因为水的汽化而形成气泡使得水柱分离，而在压力升高时，被分离的水柱将再度弥合互相撞击，形成压力急剧上升，这时将对输水管道产生很大的破坏作用。在对液柱分离及再弥合现象进行处理时，令Hwi为输水管道i断面的绝对压力，Hv为水的汽化压力，如果Hwi≤Hv，则判断产生液柱分离。此时该界面就被当作特殊的内部边界条件来处理，其压力保持为水的汽化压力，则有

式中:HP为测压管水头;Zi为该截面顶点到基准面的高度;Ha为大气压水柱高度。

将式(2.2-108)～式(2.2-111)联立即可求得液柱分离过程中的水柱分离的体积、水力要素等的变化过程。

若计算过程中液柱分离的体积等于(或小于)零，则认为水柱在该瞬时已经发生弥合相互撞击，此时可以认为液柱分离的体积为零，且QP=QPU，则由式(2.2-108)和式(2.2-109)即可求得液柱分离及其再弥合的水力参数变化过程。

2.2.3.2.2　变波速模型

变波速模型假定液体中的气泡均匀地分布在整个管段，或者把输水管道分成若干管段，每个管段拥有着各自的气体含量，在任何Δx小段上以自由气泡形式存在的气体的体积取决于该小段上的绝对压力，所以每小段上液体的有效体积弹性模量会随压力而变化，而管中的波速又随着液体的体积弹性模量而变化，从而每小段的波速就会随该段的压力而变化，反过来水击波速的变化又会影响水力过渡过程中的输水沿线压力的变化。在水力过渡流的计算过程中需要根据水击压力以及空穴率的变化来计算输水系统不同截面的水击波速变化，并采用变化的水击波速来进行水力过渡过程的计算。

假设水中一小部分以自由气泡形式存在的气体体积百分比为ε(%)，气液的总体积为V，水的体积弹性模量为K，那么使水的压力增加ΔP时，水的体积变为V1:

假设把气体的体积分配到小气泡中去，并假设气体体积的变化是等温的，那么在增加压力ΔP时，气体的体积Vg为

所以气液混合物的体积VZ为

略去很小的项，所以式(2.2-114)可以简化为

体积应变可以简化为

从而可以得到水的体积弹性模量为

所以把式(2.2-117)代入式(2.2-118)可得

由于考虑管道弹性的水的有效体积弹性模量K1存在以下关系:

把式(2.2-119)代入式(2.2-117)，并忽略掉自由气体的质量，则得到考虑气泡影响的可变波速公式为

式中:P为水的绝对压力;γw为水的重度。

式(2.2-120)适用于气泡在水中只占很小的一部分，且均匀地分散在水中的气液两相流。水中掺入不同的空气含量，处于不同的管道压力下就会得到不同的水击波速。在管道压力一定的情况下，根据水中不同的气体含量(包括溶解气和自由气泡)就可以计算出相应的水击波速值，从而就可以绘制出在压力一定的情况下水击波速与水中含气量之间的关系曲线。卡普勒斯、小堀等对此进行了研究，发现了一个十分有意义的现象:管路内的掺气量超过30%时，水击波速一般就降到约100m/s，且波速曲线具有一个宽而平的最低值区;在气体含量占30%～70%之间时，在一定管道压力下，水击波速基本为一定值，且这一最低值区随着管路压力的降低而降低;但气体含量为0%和100%的两个极限值与压力无关。水击波速与气体/液体体积比之间的关系如图2.2-12所示，从图上可以看到曲线与管道的压力有关，不同的压力有不同的曲线。据水中含气量与水击波速的关系，在一定压力下随着含气量的增加波速曲线降低很快，最终进入一个宽而平的最低值区，据可变波速公式计算出的波速值一旦落到这一区域就可采用这一最低波速值。

对于变波速模型可以用本书第2章介绍的插值特征线法进行求解，在差分格网中，时间步长Δt的选取应满足Courant-Lewy稳定性准则:

这样可以避免由于流速v和波速a的变化引起的数值的不稳定，式(2.2-121)中的波速a可以取管路中的任何一处的最大波速，根据经验，Δt可以取为0.95×乘以系数0.95是为了保证波速变化时数值计算的稳定性。

2.2.3.2.3　模型适用条件

变波速模型和自由气体-液体离散模型的理论基础仍然是前文所介绍的有压过渡流的控制方程，因此在采用这两个模型对气泡均匀分布的气液二相瞬变流进行模拟时，必须要满足自由气泡在水中只占很小的一部分的假设。因为如果水中的气体百分比很大，那么有压非恒定流的控制方程以及变波速方程中的某些假设就会为气液二相瞬变流的模拟带来较大的误差。

因此，变波速模型和自由气体-液体离散模型对于气液二相瞬变流的模拟，只适合于气体含量较小的情况，在气体含量较小时，两者计算的结果符合得较好，采用这两种方法进行气液二相瞬变流的模拟均可得到比较满意的结果，对于气体含量较大，且均匀分布于水中的情况来说，使用这两种模型均会带来较大的误差。

另外，如果气体在水中并不是均匀分布，而是有一些大的气泡处于管路的某些位置上。这种情况下，变波速模型就不再适用，因为变波速模型是基于气泡在水中均匀分布的假设，而推导出的变波速公式。而自由气体-液体离散模型正是基于气泡位于断面上的假设，所以对于这种情况仍然适用。因此，对于管路中的高位部分聚结有蒸汽穴、液柱分离及再弥合问题等，自由气体-液体离散模型的假设是合理的，因此可以采用该模型进行水力过渡过程的求解。

这两种模型均是从宏观的观点来研究问题，而不是把注意力集中在气液二相流的气泡动力学上，因此这两种模型的处理均是基于一定假设的。至于气液二相瞬变问题求解的精确性尚有以下问题需要确定。

(1)管路中的自由气泡含量的确定。

(2)作为压力和时间函数的气泡从水中释放和再吸收的规律性。

(3)管路中的气泡的运动规律和分布。

这些问题均涉及气泡动力学的理论知识，而且有些问题尚无确切的定论，还有待于进一步的研究。