方成跃 赵观辉

中国舰船研究设计中心，湖北武汉 430064

0 引 言

压力管路中，因流速剧烈变化引起动量转换从而在管路中产生一系列急骤的压力交替变化的水力撞击现象，称为水锤现象，简称水锤。就本质而言，水锤就是管道瞬变流动中的一种压力波，其产生是由于管道中某一截面内的流速发生了改变，从而使该处的压力产生一个突然的跃升或下跌，这个压力的瞬变波就称为水锤［1］。

在管内流体不发生液柱分离且各性质均匀时，水锤波速一般取常数。但在实际工程中，很多因素都会影响到水锤波速，如管内出现气泡使流体性质发生变化，或管路中设置的波纹管和软管等特殊管路元件造成的影响等。同时，管路的支撑方式也会对波速造成一定的影响。而水锤波速的大小与压力波的幅值和周期均有直接关系，因此，管内水锤波速变化的影响分析对于水锤现象的研究具有重要意义。文献［2］针对不同关阀情况下整管波速变化的影响进行了分析，但对关阀时间与管路内水锤周期的关系没有进行系统说明。文献［3］对波纹管引起的局部波速变化情况采用取平均波速的方法进行了分析，但忽略了波的散射现象。本文将借助流体软件FLOWMASTER对以上两种典型情况进行仿真分析，在改善上述问题的前提下研究不同情况下波速变化对管路内水锤现象的影响。

1 水锤波速

1.1 波速计算方法

水锤波速理论发展至今，应用比较广泛且也比较符合实际工程计算的是弹性水锤波速理论公式［4］，具体计算公式如下：

式中，a为水锤波速；K为流体的体积弹性模量；ρ为流体密度；D为管径；E为管壁材料的弹性模量；e为管壁厚度；c1为管道支撑相关系数。

上式适用于性质均匀的单相水，若考虑到水中含气的情况，流体的体积弹性模量K值将减小［5］。当微小气泡均匀地分布于薄壁水管中时，a值可按下式计算：

式中，ρ为水的密度；ρ′为空气的密度；W 为混合液的体积；W′为混合液中气泡的体积；为空气的容积掺入率；K为水的体积弹性模量；K′为空气的体积弹性模量。

由式（1）和式（2）可知，管内水锤波速不仅与管道自身性质，如几何尺寸、材料弹性模量等有关，还与管道的支撑方式有关，同时，管内流体的属性，如含气率等也对波速具有重要影响。

1.2 波速对压力波传播特性的影响

水锤波速作为管内压力波的重要参数，对压力波的影响主要体现在压力波幅值与传播周期两个方面。根据Joukowsky的经典水锤计算公式，可得：

式中，ΔP为管内压力变化；ΔV为管内流体速度变化。

管内水锤压力波动幅值与波速成正比，因此，当波速降低时可有效缓解水锤现象，改善管路发生水锤时的管壁受力情况。

根据经典水锤理论，管内压力波动为周期变化，若管内波速恒定，其周期为［6］：

式中，L为管路长度。

2 水锤波速变化的仿真分析

FLOWMASTER是一款专业的一维工程流体管路系统分析软件，擅长于对流体管路系统进行整体分析，且带有专门的水锤分析功能。该软件具有丰富的元件库，可使仿真模型的搭建方便快捷；并且该软件能对用户关注的流量和压力等参数进行精确预测，以进一步分析管路系统的特性［7］。本文选取FLOWMASTER对不同波速变化情况进行了初步的仿真分析，并对其影响特性与规律进行了相关归纳总结。

2.1 整管水锤波速变化时的水锤现象

首先，选取最典型的情况进行分析：管内波速恒定，管路上游为恒压源，下游为阀门，通过速关阀门产生水锤现象，并比较分析不同波速下水锤压力波的传播特性。在FLOW MASTER软件中搭建如图1所示的简单管路系统，各元件具体参数如表1所示。

通过对弹性圆管中的波速进行设置，分别对不同波速下的水锤现象进行仿真。由于水锤为阀门快速关闭引起，因而选取阀前节点2为分析对象，研究其压力波动特性。不同波速下节点2的压力波动如图2所示。

图1 简单管路系统图Fig.1 Simple pipe system diagram

表1 管路系统元件参数与仿真环境变量Tab.1 Components parameters and simulation environment variables of the pipe system

从图中可以直观地看出，波速的变化对管路内的水锤压力波影响较大，各波速下压力波的具体特性如表2所示。

表2 不同水锤波速下压力波特性Tab.2 Pressure wave characteristics at different

图2 节点2的压力波动图Fig.2 Pressure wave of node 2

由表中数据可得，波速与压力波周期成反比，与式（4）吻合。但水锤波速从1000 m/s变化到500 m/s时，压力波的最大幅值几乎相等，这是因为该波速对应的压力波反射周期TR大于关阀时间 Tc=0.5 s，即该水锤现象属于“缓慢事件”［8］，故压力波幅值应使用下式求解：

根据上式，当水锤波速由1000 m/s减小到500 m/s时，对应的发射周期从0.2 s增加到了0.4 s，其他参数未变，因而其最大幅值的理论计算值不变，这与仿真结果基本相符。

当水锤波速降至100 m/s（通常发生在管内出现大量气泡的极端情况下），压力波反射周期TR小于关阀时间Tc（该水锤现象属于“快速事件”）时，此时，可采用式（4）计算对应幅值。由式（4）和式（5）计算可得，幅值应为原波速下的一半。但仿真结果却与此有较大差值，幅值与理论计算结果相比偏小，这是因为其峰值出现时间较晚，压力波能量在管内由于摩擦等因素损耗增大。

综合本小节分析结果可知，水锤波速下降时，对应的压力波周期会延长，且严格成反比关系。但是，若压力波的反射周期大于关阀时间，则波速的变化对压力波幅值基本不造成影响，当反射周期小于关阀时间时，压力幅值才会出现明显的下降趋势。

2.2 管路不同位置局部波速变化时的水锤现象

在实际管路系统中，基于某些具体设计要求，通常需要在管路中添加一些特殊的管路元件，如波纹管、软管等，同时，管路中某些位置可能会出现汽化或空气渗入，这些都会造成管路局部的波速出现较大变化。

针对波纹管引起的局部波速变化情况，文献［4］通过分别计算钢管与波纹管内的波速，然后以求平均的方法来计算整管波速，进而再使用平均波速对水锤现象进行分析。计算公式如下：

该种方法未考虑波纹管的具体位置，其将局部波速的变化转化为了全管平均波速的变化。本节将采用FLOWMASTER对波速变化的具体位置所带来的影响进行分析。

搭建类似于图1的管路系统，将单管换成由3个弹性管道元件组成的串联管路，在保持总长度不变的情况下，调整其中单个管道的波速，从而实现控制管路局部波速的目的。系统如图3所示，管路的参数大致与图1所示的系统相同，新增管路元件与变量更改信息如表3所示。

2.2.1 短管置于管路末端

首先，将短管置于阀门前，即管路末端，分别对短管内波速为1200，200 m/s的两种情况进行仿真分析，得到阀门前节点5的压力波变化情况如图4所示。

图3 局部波速变化试验管路Fig.3 Pipe system with wave velocity changes in local position

图4 短管位于管路末端时节点4的压力波动图Fig.4 Pressure wave of node 4 when short tube located at the end of the pipe

由图可见，短管波速降低后，节点5处的压力波幅值出现了一定程度的下降，而且压力波周期明显延长。若采用传统的平均波速计算方法，由式（4）和式（6）可得波速变化前后周期分别为T=0.2 s，Tave=0.24 s。与图中结果相比，发现Tave明显比仿真结果偏小，这是因为不同波速的管道接口处会发生压力波的散射，即压力波在不同波速交界处会出现反射与穿透现象。

根据文献［8］和文献［9］中的理论，不同波速管道内的流体相当于两种不同的介质，介质的特性与流体的密度以及波速有关，将其用Z表示，计算公式如下：

在如图3所示的管路系统中，流体密度保持不变，波速在短管处发生变化，因而会产生波的散射现象，其原理如图5所示。

图5 压力波散射原理图Fig.5 Scattering schematic of pressure wave

当压力波如图5所示从介质一向介质二传播时，在界面处发生散射，其穿透系数Tr与反射系数R的计算公式为：

由式（8）、式（9）可得，波的穿透系数总为正数，即当波穿越不同介质界面时，相位保持不变；而波的反射系数则可正可负，当Z2＞Z1时，R为正数，反射波相位保持不变；当Z2＜Z1时，R为负数，反射波相位出现180°偏移，即相位反转。

2.2.2 短管置于管路中段

将短管位置调节至管路中段，参考上节步骤重复仿真过程，得到新的节点压力波动，如图6所示。

图6 短管位于管路中段时节点4的压力波动图Fig.6 Pressure wave of node 4 when short tube located in the middle of the pipe

对图6进行分析可发现，将短管位置由管路末端移至中段后，压力波的传输特性出现了较大变化，波形较为复杂，且幅值变大，幅值变化速率增加，管路动态特性出现明显的恶化情况。这充分反应了波速变化位置的重要性，以及采用平均波速来分析水锤现象的不合理性。

压力波的波形出现巨大变化的主要原因仍然是由于不同界面处的散射现象，从而导致界面处产生了多个新的压力波，使波形严重复杂化，具体原理如图7所示。

由图7可知，对于位于中段的短管，由于散射会产生大量新的压力波，会对管路内的压力波产生较大干扰，进而导致压力波幅值增大，动态特性变差等不良后果。

2.2.3 短管置于管路前端

将短管位置移至管路前端，参考上节步骤重复仿真过程，得到新的节点压力波动图如图8所示。

图7 短管位于管路中段时的压力波散射原理图Fig.7 Scattering schematic of pressure wave when short tube located in the middle of the pipe

图8 短管位于管路前端时节点4的压力波动图Fig.8 Pressure wave of node 4 when short tube located in the front of the pipe

由图8可知，当短管位于管路前端时，只有一个散射截面，散射情况较简单，波速降低后波形仍然较规则。由于短管的距离较短，波在管内传播的时间也相对较短，因而压力波周期基本无变化。但是，发射波的叠加效应仍然使压力波的幅值出现了增大的情况，因而当管路前端出现波速局部下降时，不会改善水锤现象。

综上所述，由于管内不同波速界面处的散射现象，使得局部波速变化对水锤现象的影响较为复杂，需要考虑波速变化的具体位置。在关阀引起的水锤现象中，阀前位置的局部波速下降能明显改善管路水锤现象，优化管路动态特性；而在管路中段与前端出现波速下降时，改善效果不明显，甚至还会出现恶化现象，尤其是当波速下降区域位于管路中段时，效果最为恶劣。

3 结 论

水锤波速作为水锤现象中的关键参数之一，受管道材料、尺寸、支撑方式以及管内流体性质等多种因素的影响，因而分析波速变化对水锤现象带来的影响具有很好的实际工程意义。本文通过流体仿真软件FLOWMASTER，分别对典型管路系统水锤现象中整管波速变化以及不同位置的局部波速变化的影响进行了分析，得出如下结论：

1）整段管路水锤波速发生变化时，水锤压力波周期与波速大小成反比，即降低水锤波速可有效改善管路的动态特性。

2）当管路水锤为“缓慢事件”时，水锤波速的改变对压力波幅值影响较小；当管路水锤为“快速事件”时，水锤波速与压力波幅值成正比，此时，降低波速能有效缓解水锤现象。

3）管路中出现局部波速变化时，在波速变化界面处会出现压力波的散射现象，使管内压力波发生明显变化。

4）局部波速的变化对水锤现象的影响与其位置密切相关，在关阀引起的水锤现象中，阀前波速下降能显著改善管内水锤现象，而在管路前端与中段则效果有限，甚至会使管内水锤现象恶化。

［1］王树人.水击理论与水击计算［M］.北京：清华大学出版社，1981.

［2］YANG Y S，MA X Y，LI X Y.Analysis the impact of wave velocity on the calculation of pressure boost in water hammer of long-distance water pipelines［C］//Water Resource and Environmental Protection（ISWREP）.Xi’an，2011.

［3］王英浩.火箭推进剂输送系统中波纹管对水击压力的吸收分析［D］.长沙：国防科学技术大学，2004.WANG Y H.Analysis of bellows absorb water hammer in rocket propellant conveyer pipeline［D］.Changsha：National University of Defense Technology，2004.

［4］怀利E B，斯特里特V L.瞬变流［M］.北京：水利电力出版社，1983.

［5］夏小娟.有压瞬变流模型试验中水锤波速控制方法的研究［D］.南京：河海大学，2005.XIA X J.Study on methods of controlling water-hammer wave velocity in model tests for fluid transients［D］.Nanjing：Hohai University，2005.

［6］GHIDAOUI M S，ZHAO M，MCINNIS D A，et al.A review of water hammer theory and practice［J］.Applied Mechanics Reviews，2005，58（1）：49-76.

［7］王勇.供水系统水锤数值计算及动态模拟［D］.合肥：合肥工业大学，2009.WANG Y.Numerical calculation and dynamic simulation of water hammer in water supply system［D］.Hefei：Hefei University of Technology，2009

［8］VILBERG I K.Airbag for piping systems［D］.Trondheim：NTNU，2010.

［9］LI Y，ALESSANDRA B M，KIM H，et al.Variation of wave speed determined by the PU-loop with proximity to a reflection site［C］//33rd Annual International Conference of the IEEE EMBS.Boston，2011.