安装有空气阀的输水管路系统空管充水过程瞬态分析

王 玲，王福军，2，黄 靖，罗建群，桂新春，谢爱华

（1.中国农业大学 水利与土木工程学院，北京 100083；2.北京市供水管网系统安全与节能工程技术研究中心，北京 100083；3.株洲南方阀门股份有限公司，湖南 株洲 412007）

1 研究背景

空管充水过程对于新建或更新改造的输水管路系统在投入正常运行之前是不可或缺的环节［1-3］。曲折管线进行空管充水工作时，易在管路的局部高点处产生水柱分离现象，随之而来的水柱弥合会产生剧烈的压力波动［4］，因此，有必要在管线布置一些进排气阀实现对管路的压力防护［5-7］。为了有效预测空气阀在管路系统中的水锤防护特性，通常需要进行较为准确的数值模拟和选型计算，但现有空气阀的应用研究多集中于满管流状态［8-10］，对于空管充水阶段空气阀的水锤防护研究较少。同时，目前用于求解空气阀数学模型的常用方法主要有直接求解法和牛顿迭代法。其中，直接求解法是指Wylie和Streeter提出的方法［11］，其将以亚声速流入与流出和以临界速度流入与流出空气阀的4种情况下的空气质量流量方程离散化，并用一系列分段方程近似，结合空气阀边界的特征线方程以及空气阀内的气体状态方程，可转换成关于空气阀内气体压力的一元二次方程，然后通过判断解的存在区域并求出下一时刻的空气阀内的气体压力；牛顿迭代法［12］是指根据当前空气阀内的气体压力判断选择通过空气阀的空气质量流量方程，然后与空气阀边界的特征线方程和空气阀内的气体方程组成非线性方程组，通过牛顿迭代法去求解下一时刻的空气阀内的气体压力。相比较而言，直接求解法具有求解速度快、计算精度高等优点，但只适用于满管流状态；而牛顿迭代法计算过程简单，虽然对于满管流计算和空管充水过程计算均适用，但计算精度易受收敛条件的限制，并且可能存在计算得到的下一时刻的空气阀内的气体压力不在当前时刻判断的气体进出流压力区域的问题。

为了准确模拟空气阀在空管充水过程中的瞬态特性，本文首先对求解空气阀模型的牛顿迭代法进行改进，然后结合空气阀模型的直接求解法，确定一种基于改进牛顿迭代法和直接求解法相结合的求解空气阀模型的计算方法。在此基础上，对曲折管线的空管充水过程进行瞬态分析，并对3种类型的空气阀在空管充水过程中的水锤防护特性进行比较，最后确定空气阀的选取类型。

2 数学模型

2.1 空管充水模型为了追踪空管充水过程中水气交界面的移动位置，本文首先假设：（1）水气交界面垂直于管轴线；（2）管道出口连通大气，充水过程中，管道内气体压力假定为恒定的大气压值。其中，对于第1条假设，在管径较小、充水水头足够大或充水流速足够快的管路系统中，水气交界面的形状是比较陡的，可认为是垂直的水气交界面。另外，由于水气交界面在整个充水水柱中所占比例较小，水气交界面的形态对充水瞬变流的整体影响并不十分明显［14-15］。对于第2条假设，由于本文所研究的管道出口是敞开式的，管内气体几乎无阻碍排出，即使管内大气受到充水水柱的压缩，其压力波动也是极其微弱的。因此，待充水管道内的气体压力假设为大气压是可以接受的。基于以上假设，采用Malekpour［13］建立的界面追踪法求解空管充水瞬变流问题。界面追踪法是指在模拟空管充水过程中，时刻计算与判断水气交界面的位置，这样不同的流动区域有各自对应的控制方程。

如图1所示，对于当前时刻已充满的计算单元，建立特征线方程［16］：

式中：B为计算常数，B=a（/gA）；H为测压管水头；Q为流量；CP为正特征线方程系数；CM为负特征线方程系数；a为水锤波速；g为重力加速度；A为管道横截面面积；下标i和t+Δt分别为管道分段节点和下一时刻点。

对于正在充水的计算单元，建立水柱前缘对应的一维连续方程和动量方程［13］：

式中：S为能量坡度，ψ为权重系数；Δt为时间步长；δx为当前时间步长内水柱行进距离，Lf为当前时间充水单元内水柱长度；S0为管道坡度；Z为管道节点高程；f为管道摩擦系数；D为管道直径；下标t为当前时刻点。

对于待填充的计算单元，由于连通大气，故管内压力假定为大气压。

为了简化上游水库与控制充水流速阀门之间的关系，把两者当成一个整体作为统一的上游边界，对应的边界控制方程为［13］：

式中：KS为上游水库至阀门间水柱的水头损失系数，KS=fLS/D；KV为阀门损失系数；LS为上游水库至阀门间的水柱长度；HR为上游水库水位。

考虑到空管充水过程中可能伴随的水柱分离现象，本文采用离散气体空穴模型（DGCM）［17-18］进行水柱分离现象的预测与计算，空穴内的气体控制方程和连续方程为

式中：ρ为水的密度；HV为水的表蒸汽压；Vg为发生水柱分离时的空穴体积；P0*为初始参考压力；α0为初始体积分数；VR为管道计算单元体积；下标iu和id分别为节点i的进、出口标号。

至此，对于已充满水柱的管道内部节点，联立式（1）、式（2）、式（6）和式（7），通过转换并求解关于气穴内气体压力的一元二次方程，计算得到节点的H，然后反代求解Qiu、Qid和Vg；对于正在充水的管道节点，联立式（1）、式（3）、式（4）、式（6）和式（7），采用牛顿迭代法，计算得到节点的H、Qiu、Qid、Vg和δx。

2.2 空气阀模型及求解方法在建立空气阀模型时，通常需要考虑如下假设［19-20］：（1）气体等熵流入流出空气阀；（2）空气阀内的气体遵守等温规律；（3）进入管内的气体集聚在空气阀位置处；（4）液体表面高度基本不变，而气体的体积相比管段内液体体积很小。同时，空气阀的状态根据空管充水的进展过程，如图1所示，可以分为3种情况：（a）充水水柱未到达空气阀位置处，空气阀的作用是排掉管路内的滞留气体，考虑到管路出口连通大气，因此，可认为空气阀处的压力基本保持为大气压值；（b）充水水柱越过空气阀所在位置，并且充水水柱前缘位于靠近空气阀出口管道的第一个计算单元；（c）充水水柱早已越过空气阀所在位置，即靠近空气阀进、出口管道的计算单元皆处于满管流状态。

图1 安装有空气阀的空管充水计算网格

对于图1（b）和图1（c）两种情况，当空气阀内不存在气体且水压高于大气压时，空气阀接头处的边界条件就是一般内截面解。当水压降到大气压以下时，空气阀开启，气体流入，在气体排出之前，气体满足定温气体方程和连续方程［11］：

式中：Pair为空气阀内的气体绝对压力；Vair为空气阀内的气体体积；m为空气阀内的气体质量；ṁ为空气阀内的气体质量流量；R为气体常数；T为空气阀内的气体绝对温度；下标P0和Pn+1分别为空气阀的进、出口标号。

同时，空气阀处的气体压力与水柱的测压管水头满足如下关系

式中：γ为水的比重；Hatm为绝对大气压头。

根据气体是否为音速或亚声速流入流出空气阀，对应的空气质量流量可分为如下4种情况［16］：

（1）气体以临界速度流入空气阀，

（2）气体以亚声速流入空气阀，

（3）气体以亚声速流出空气阀，

（4）气体以临界速度流出空气阀，

式中：Cin为气体流入系数；Cout为气体流出系数；Ain为气体流入孔口面积；Aout为气体流出孔口面积；ρair为气体密度；Patm为绝对大气压。

这样，对于进、出口计算单元已充满水柱的空气阀节点，采用Wylie等提出的直接求解法求解空气阀内的气体压力［11］，即首先对气体质量方程式（12）和式（13）进行离散，形成一系列关于空气阀内气体压力的二次方程，然后将离散化的式（12）和式（13），以及式（11）和式（14）与式（1）、式（2）、式（8）、式（9）和式（10）联立，形成最终的关于空气阀内气体压力的一元二次方程，通过判断解的存在区域，从而求得空气阀内气体压力Pair，最终反代求解HP0、HPn+1、QP0、QPn+1和Vair；对于进口计算单元已充满水柱而出口计算单元正在充水的空气阀节点，联立式（1）、式（3）、式（4）、式（8）、式（9）、式（10）和式（11）到式（14）中之一，采用改进牛顿迭代法，计算得到节点的HP0、QP0、QPn+1、Pair、Vair、ṁ和δx。基于改进牛顿迭代法和直接求解法相结合的求解空气阀模型的计算流程图如图2所示。

其中，空气阀模型的原始牛顿迭代法是指根据当前时刻计算得到的气体压力，判断气体质量流量的进出状态，根据式（11）—（14）确定所在状态的气体质量流量控制方程，然后与气体方程和特征线方程等组成非线性方程组，通过迭代求解，获得下一时刻的气体压力。空气阀模型的改进牛顿迭代法则是在获取由原始牛顿迭代法计算的下一时刻的气体压力后，接着根据所得到的下一时刻的气体压力再次判断气体质量流量的进出状态。若前后两次判断的气体质量流量的进出状态相同，则所计算的下一时刻的气体压力为正解；若前后两次判断的气体质量流量的进出状态不同，则所计算的下一时刻的气体压力将作为当前时刻的气体压力重复以上步骤进行迭代求解。当完成对4种气体质量流量的进出状态判断后仍不满足要求，则认为空气阀内无气体，空气阀节点压力按有压内截面进行求解。

例如，若当前时刻计算的空气阀内的气体压力（Pair）t=0.5Patm，则由此气体压力确定的气体质量流量控制方程为式（11），与空气阀边界的其它方程组成非线性方程组，通过原始牛顿迭代法求解获得的下一时刻的空气阀内的气体压力可能为（Pair）t+∆t=0.8Patm＞0.528Patm，这就导致采用式（11）计算的（Pair）t+∆t不在式（11）应用的压力范围内。而当采用改进的牛顿迭代法后，在获得下一时刻的空气阀内的气体压力可能为（Pair）t+∆t=0.8Patm时，则令（Pair）t=0.8Patm，重新确定的气体质量流量控制方程为式（12），通过非线性方程组的再次迭代计算，求解获得的下一时刻的空气阀内的气体压力可能为0.528Patm＜（Pair）t+∆t=0.85Patm＜Patm，这说明采用式（12）计算的（Pair）t+∆t在式（12）应用的压力范围内，故（Pair）t+∆t=0.85Patm为正解。由此可见，相比于空气阀模型的原始牛顿迭代法，改进牛顿迭代法可有效避免气体压力计算前后判断的气体质量流量进出状态不一致的问题。

2.3 空气阀模型求解方法的验证为了说明本文提出的基于改进牛顿迭代法和直接求解法相结合的求解空气阀模型的计算方法的可靠性，采用Apollonio等［21］的带有空气阀的空管充水测试系统进行数值验证。该测试系统主要由供压水塔、蝶阀和管径为DN75的U-PVC管组成。其中，供压水塔提供的充水水头为16.3 m，向上倾斜的管长为5.43 m，向下倾斜的管长为5.90 m，管道坡度为30°。口径为6.4 mm的空气阀安装在管路系统的高点位置。

图2 求解空管充水过程中空气阀模型的计算流程

图3 空气阀位置的管内压力试验值与模拟值对比

图3为空气阀位置处管内压力p的试验结果与数值结果。对比发现，由试验测得的压力变化和数值模拟获得的压力变化整体趋势较为一致，特别是压力波动稳定后的压力曲线基本是重合的。但试验结果与数值结果较大的差异发生在约2.0 s时刻，具体表现为由原始牛顿迭代法获得的压力出现数值振荡，最大压力达到了1.52 bar，远大于最大试验压力值0.76 bar，最小压力达到了负压-0.31 bar，而试验压力无负压；由改进牛顿迭代法与直接求解法相结合的计算方法获得的压力波动缓慢，最大压力为0.48 bar，比试验值稍小，这可能是由于本文的空管充水模型忽略了充水过程中空管内气体的压缩性导致的，而在实际测试过程中，管内水柱因气堵作用会有压升的波动过程。因此，在忽略管道内气体压缩性的条件下，由改进牛顿迭代法与直接求解法相结合的计算方法获得的压力变化更接近于测试结果。

3 研究对象

现以一曲折管线系统为例，采用所建立的空管充水模型进行充水瞬变流的计算与分析。如图4所示，该管路系统主要包括1个定水位的上游水库，其水位HR为102.0 m；1个安装于管道进口的控制阀，其全开损失系数为0.8；3条等管径管线，其直径为DN100，摩擦系数f为0.02，长度L分别为300 m、100 m和150 m。管线出口通大气，管路的节点高程Z见图4中所示。初始时刻，控制阀处于关死状态，阀前为上游水库静水压，阀后为大气压。当控制阀快速开启后，在上游水库压差作用下，管线开始充水。空管充水瞬变流计算时，水锤波速a为1000 m/s。控制阀快速开启后，充水水柱的长度Lw和流速vw变化如图5所示，管线膝点B处发生水柱分离时的压力PB和气穴体积Vg变化如图6所示，从充水开始至整个管线充满的过程中，充水瞬变流引起的测压管水头包络线如图7所示。

图4 待充水的曲折管线示意图

图5 管线内充水水柱长度和流速变化

图6 管线膝点B处压力和气穴体积变化

图7 管线测压管水头包络线

图5至图7的计算结果表明，整个充水管线在258 s时刻完成充水工作，在控制阀开启的瞬间达到最大充水流速4.0 m/s；充水水柱越过膝点B时，会产生水柱拉断，形成空穴，空穴最大体积可达0.022 m3，随之而来的水柱弥合引发的高压（测压管水头减去高程）水锤为27.1 m，大约为管路进口淹没水深的13.6倍；从测压管水头包络线可以看出，由水柱分离再弥合引发的高压水锤可在水锤波的作用下，传至几乎整个管线，较大的压力变化易对管路产生冲击，并引发管路振动，从而威胁了管路的安全运行。

4 3种类型空气阀的空管充水瞬变流分析

为了预防或降低曲折管线中快速充水伴随的水柱分离及再弥合水锤危害，在图4所示的管线膝点B处布置了一个DN20的空气阀。根据空气阀的类型不同，设置了3种类型空气阀的布置方案，分别为进气阀、普通进排气阀和进气微排阀。其中，进气阀是指只进气不排气的空气阀；普通进排气阀是指等口径（20 mm）进排气的空气阀，进气微排阀是指大口（20 mm）进气，小口（2 mm）微量排气的空气阀。在所建立的含有空气阀模块的空管充水模型的基础上，分别采用原始牛顿迭代法和改进牛顿迭代法与直接求解法相结合的计算方法获得了3种类型空气阀的水锤防护特性的数值模拟结果。

图8 安装进气阀的膝点B处压力和气体体积变化

图9 安装进气阀的测压管水头包络线

图10 安装普通进排气阀的膝点B处压力和气体体积变化

图11 安装普通进排气阀的测压管水头包络线

安装进气阀后，从图8的节点压力和进气阀内气体体积变化可以看出，膝点B在进气阀的补气作用下，压力维持在了-2.0 m以上，并且进气后气体体积最终维持在了0.36 m3；从图9的测压管水头包络线可以发现，安装进气阀后，可有效避免管线中汽化压力的出现，系统出现的最小压力极值为-2.3 m，最大压力极值相比于不安装空气阀的管线充水过程有较大程度的缓解。同时，还可以看出采用原始牛顿迭代法和改进牛顿迭代法与直接求解法相结合的计算方法对进气阀模型求解时，其计算结果基本无差异。可见，对于空管充水伴随的水柱分离及再弥合水锤防护，只进气不排气是最好的解决方案。但这样的做法在实践中似乎并不可行，因为补进去的空气最终需要慢慢排出，否则会影响管路的正常输水能力，增加水力损失，消耗能源。

安装普通进排气阀后，从图10的膝点B处压力和普通进排气阀内气体体积变化可以看出，当节点压力下降到大气压以下时，普通进排气阀开启进行补气，当节点压力回升到大气压以上时，普通进排气阀进入排气阶段，而在快速排气结束时刻，气穴两侧的水柱撞击产生高压水锤，采用原始牛顿迭代法得到的初次弥合高压为227.8 m，采用改进牛顿迭代法与直接求解法相结合的计算方法得到的初次弥合高压为110.5 m，均远高于不安装空气阀时膝点B处水柱分离引发的弥合高压水锤；从图11的测压管水头包络线也可以看出，采用原始牛顿迭代法和改进牛顿迭代法与直接求解法相结合的计算方法预测的整个管线均出现汽化压力和弥合高压。可见，采用等口径的进排气阀，并不能有效消除负压，并且大口径排气还会导致水柱弥合高压水锤，并传播至整个管线。

另外，从图10（a）中的通过原始牛顿迭代法获得的膝点B处的压力变化可知，当空气阀内气体体积排尽为零后，仍有极端的高压值生成，这并不是真实有效的物理值，而是由原始牛顿迭代法带来的虚拟数值振荡，并在压力波的传播下，导致图11（a）中的测压管水头包络线同样存在不真实的压力极值。然而，从图10（b）中的通过改进牛顿迭代法与直接求解法相结合的计算方法获得的膝点B处的压力变化可知，空气阀排尽气体时产生弥合高压，然后在阻尼的作用下压力峰值削弱，这与实际物理过程更为相符，由此形成的图11（b）的测压管水头包络线也更为接近实际。

安装注气微排阀后，从图12的膝点B处压力和注气微排阀内气体体积变化可以看出，在微量排气阀的作用下，在膝点B处由水柱分离引发的弥合水锤高压，采用原始牛顿迭代法计算的数值为10.7 m，采用改进牛顿迭代法与直接求解法相结合的计算方法计算的数值为12.3 m，相比于普通进排气阀均大大降低；由图13的测压管水头包络线还可以看出，在进气微排阀的补气作用下，不仅弥合水锤高压在整个管线大为降低，同时负压防护有较大缓解，最小压力极值由原始牛顿迭代法计算的数值为-4.1 m，而由改进牛顿迭代法与直接求解法相结合的计算方法计算的数值为-2.6 m。这是由于选取微小排气孔的空气阀，可使排气过程更缓慢，从而在管道内产生弹性空气气囊，延长水柱弥合时间，进而减小水流速度的变化率，最终有效的降低压力波动。同时，通过数值对比可以发现，采用原始牛顿迭代法和改进牛顿迭代法与直接求解法相结合的计算方法对注气微排阀模型求解时，对于空气阀内气体排尽后的压力波动，虽然数值上存在较小差异，但均没有出现虚假数值振荡问题，可认为由两种方法预测的压力值均较为合理。

图12 安装注气微排阀的膝点B处压力和气体体积变化

图13 安装注气微排阀的测压管水头包络线

5 结论

（1）改进了求解空气阀模型的牛顿迭代法过程，结合直接求解法，确立了空管充水过程中空气阀模型的求解方法，提高了空管充水过程中空气阀模型的计算精度；（2）为预防空管充水过程中出现的水柱分离及再弥合水锤，对3种类型空气阀的水锤防护特性进行了数值研究，通过计算发现：进气阀可有效缓解空管充水过程中出现的水柱分离及再弥合高压水锤，但滞留系统的气体，易降低管路系统的输水能力；普通进排气阀虽然可以排尽管路残留气体，但对水柱分离和弥合高压的防护并不理想；进气微排阀不仅在大口进气的过程中缓解管路负压，而且在微量排气过程中抑制弥合高压水锤的形成。因此，进气微排阀对于空管充水过程中水柱分离及再弥合高压水锤的防护最为有效。

参 考 文 献：

［1］ 杨开林，郭永鑫，付辉，等.管道充水水力瞬变模型相似律［J］.水利学报，2012，43（10）：1188-1193.

［2］ 郭永鑫，吴昊潼，杨开林，等.管道充水过程气液两相流动力学特性分析［J］.南水北调与水利科技，2013，11（5）：65-69.

［3］ 曹命凯，郑源，严继松，等.长距离有压输水管道充水过程的研究［J］.长江科学院院报，2009，26（12）：54-57.

［4］ MALEKPOUR A，KARNEY B.Profile-induced column separation and rejoining during rapid pipeline fillin［J］.Journal of Hydraulic Engineering，2014，140（11）：1-12.

［5］ BALACCO G，APOLLONIO C，PICCINNI A F.Experimental analysis of air valve behavior during hydraulic tran⁃sients［J］.Journal of Applied Water Engineering&Research，2015，3（1）：3-11.

［6］ 戚兰英，石维新，张大成.南水北调中线工程PCCP及暗涵输水系统空气阀型式研究［J］.水利水电技术，2009，40（11）：40-44.

［7］ 郑源，薛超，周大庆.设有复式空气阀的管道充、放水过程［J］.排灌机械工程学报，2012，30（1）：91-96.

［8］ 刘志勇，刘梅清.空气阀水锤防护特性的主要影响参数分析及优化［J］.农业机械学报，2009，40（6）：85-89.

［9］ 张健，朱雪强，曲兴辉，等.长距离供水工程空气阀设置理论分析［J］.水利学报，2011，42（9）：1025-1033.

［10］ 李小周，朱满林，陶灿.空气阀型式对压力管道水锤防护的影响［J］.排灌机械工程学报，2015，33（7）：599-605.

［11］ 怀利E B，斯特里特V L.瞬变流［M］.清华大学流体传动与控制教研组译.北京：水利电力出版社，1983.

［12］ 杨开林.电站与泵站中的水力瞬变及调节［M］.北京：中国水利水电出版社，2000.

［13］ MALEKPOUR A.Analysis of rapid pipeline filling including column separation&entrapped air effects［D］.To⁃ronto：University of Toronto，2014.

［14］ ZHOU L，LIU D，KARNEY B.Investigation of hydraulic transients of two entrapped air pockets in a water pipe⁃line［J］.Journal of Hydraulic Engineering，2013，139（9）：949-959.

［15］ HOU Q，TIJSSELING A S，LAANEARU J，et al.Experimental investigation on rapid filling of a large-scale pipeline［J］.Journal of Hydraulic Engineering，2014，140（11）：1-14.

［16］ WYLIE E B，STREETER V L，SUO L.Fluid transients in systems［M］.Englewood Cliffs，NJ：Prentice-Hall，1993.

［17］ SIMPSON A R，BERGANT A.Numerical comparison of pipe-column-separation models［J］.Journal of Hydrau⁃lic Engineering，1994，120（3）：361-377.

［18］ BERGANT A，SIMPSON A R，Tijsseling A S.Water hammer with column separation：A historical review［J］.Journal of Fluids and Structures，2006，22（2）：135-171.

［19］ 陈乃祥.水利水电工程的水力瞬变仿真与控制［M］.北京：中国水利水电出版社，2005.

［20］ MALEKPOUR A，KARNEY B W.A quasi-two-phase flow model for calculating filling in pipelines［C］//Pro⁃ceedings of Pipelines 2014：from Underground to the Forefront of Innovation and Sustainability，Portland，Ore⁃gon，2014.

［21］ APOLLONIO C，BALACCO G，GIUGNI M，et al.Hydraulic transients caused by air expulsion during rapid fill⁃ing of undulating pipelines［J］.Water，2016，8（1）：1-12.