为 EPANET 增加代数水锤水力求解器以模拟配水系统中的非稳态流

《Journal of Water Process Engineering》：Empowering EPANET with an algebraic water hammer hydraulic solver to model unsteady flow in water distribution systems

　　

•代数水锤（AWH）已集成到开源的EPANET 3中，可在同一次模拟中完成稳态与非稳态流分析。•准自适应求解器可检测阀门和泵事件，并从G-GGA格式切换至AWH格式。•每根管道为一个计算单元，因此求解器不需要对管道进行离散化。•两个网络上的四种瞬态工况与Bentley HAMMER匹配，NSE最高达0.90。•在包含1489根管道的网络上，运行时间比Bentley HAMMER短26%。引言给水分配系统中非稳态流条件产生的极端压力对基础设施完整性和公共健康构成了严重威胁[1], [2], [3]。当发生阀门突然关闭、泵故障或需求快速变化时，由此产生的剧烈压力瞬变可能对管道造成破坏性影响[4], [5], [6], [7]。虽然不如损坏和无法使用的管道那样明显，但水锤事件可能导致污染物通过受损的管壁和管件渗入系统[8], [9]。特别是当压力达到负值时，受污染的土壤水和病原体可能通过微观裂缝侵入给水分配系统[10]。在这种情况下，公共健康可能受到严重威胁[11]。大量研究表明，非稳态流引起的极端压力事件会导致水质下降，并引发影响数千名消费者的疾病[12], [13], [14]。因此，对给水分配系统的非稳态流进行水力模拟对于提高现有基础设施的可持续性和保护公共健康至关重要[15]。除了直接后果之外，水锤事件的破坏性影响还会影响供水企业在管理和规划饮用水分配网络时面临的许多决策。在基础设施方面，反复的涌浪载荷会加速老旧主干管的疲劳，扩大已有的泄漏，松动接头，并缩短管道、阀门和管件的有效使用寿命[6], [7]。在运营方面，泵停机或快速阀门关闭后的负压阶段会在地下缺陷处产生侵入通道，这与有记录的水源性疾病暴发存在关联[10], [11], [12]。这些风险直接转化为实际的设计与规划问题，例如水锤罐、空气室或止回阀应安装在何处；如何安排泵的加速和减速序列；在操作规程中应采用何种阀门关闭规律；压力释放装置应如何确定规格；以及在老旧区域中应优先在哪里进行修复。要在事故发生之前（而非之后）回答这些问题，需要使用能够再现实际运行条件下瞬态压力幅值、时序和空间分布的水力模型。因此，非稳态流水力求解器的实用价值不仅在于再现单个涌浪事件，还在于使供水企业能够比较不同的保护策略、筛查脆弱位置，并将瞬态安全裕度纳入日常运营和资本规划决策中。如果流动特性随时间变化，则该流动被归类为非稳态流。管道中的非稳态流通常使用两种建模方法之一来处理，具体取决于是否考虑流体的可压缩性：刚体水柱法（RWC）和弹性水柱法（EWC）。在RWC方法中，假设流体不可压缩且管壁刚体。流动和扬程由保留惯性效应但忽略波浪传播的普通微分方程控制[16], [17]。该假设仅在流动特性逐渐变化时适用，例如缓慢阀门操作或日需量变化期间。当发生水锤事件时，不可压缩假设不再成立，压力波在流域中传播；因此，瞬态事件只能使用EWC方法一致地进行建模。在EWC方法中，流体被视为弹性且因此可压缩，管道流动的特征参数（如流量和压力）通过使用一维双曲型偏微分方程来确定[18], [19], [20]。因此，EWC方法是模拟给水分配系统中突发性水力事件唯一物理一致的选项[21]，本研究采用了该方法。在过去的三十年中，EPANET已成为最广泛用于给水分配系统水力和水质分析的开源软件[22], [23]。该程序中的稳态和扩展期模拟使用全局梯度算法（GGA）进行，其中Newton-Raphson方法应用于网络的流量守恒和能量守恒方程[24], [25], [26]。由于其卓越的收敛能力和计算效率，GGA是一种稳健的方法，构成了大多数商业水力建模软件的核心计算引擎[27], [28], [29]。EPANET 3.0版本由独立的科研人员社区和程序员开发，使用面向对象的编程语言C++生成（https://github.com/OpenWaterAnalytics/epanet-dev）。该版本包含了压力依赖需求建模[30]、动态水箱液位计算[31]和基于扬程的收敛准则[32]。其模块化设计使其成为与水力建模相关扩展的适宜平台[16], [33], [34], [35], [36]。自GGA最初引入以来，已对其进行了许多改进。广义GGA（G-GGA）方法可作为一个例子。在该方法中，水箱节点的流量计算得到改进，并消除了经典GGA在水箱靠近时表现出的不稳定性[37]。随后，通过保持GGA的收敛特性和矩阵结构，稳态流水力求解器被赋予了解不可压缩非稳态流水力问题的能力。用于不可压缩非稳态流条件的水力求解器称为RWC-GGA[16], [17], [38]。最近，特征线法（MOC）的相容方程通过代数水锤（AWH）方法被转化为GGA风格的代数方程系统。因此，通过利用GGA方法的高收敛性能和效率，可以在给水分配系统中进行可压缩瞬态分析[39]。AWH之前的所有阶段——GGA、G-GGA和RWC-GGA——提供了一种将属于GGA算法族的AWH方法用于水锤分析并集成到EPANET中的自然方式。尽管存在这种算法渊源，但至今尚未有开源的瞬态求解器被直接集成到EPANET的计算核心中。现有的基于Python的工具（如TSNet[40]及其并行化后继PTSNet[41]）作为外部MOC引擎运行，而非EPANET原生求解器的扩展。此外，这些软件包存在拓扑约束，例如无法将阀门和泵连接到水库，而在EPANET中不存在此类约束。需要强调的是，将瞬态求解器集成到EPANET 3中并非易事。EPANET的稳态水力求解器（G-GGA）假设每根管道有一个流量值，而在AWH方法中，必须为每根管道的上游端和下游端分别定义流量值。由于必须求解每根管道的特征方程，必须利用这些管道中的回溯流量和扬程。这只能通过在每个时间步持续记录节点扬程和管道流量并在需要时检索记录数据来实现。本质上，通过为每根管道定义两个流量值并持续记录和利用过去的扬程和流量数据，AWH方法可以在不对管道进行分段的情况下对给水分配系统的非稳态流进行水力分析[39]。基于Nault等人[39]开发的AWH方法，EPANET-AWH是由将AWH方法直接集成到面向对象的EPANET 3中而开发的开源准自适应扩展。在正常运行条件下，求解器以稳态G-GGA模式运行。当检测到突发性水力事件时，会自动切换至AWH求解器，剩余模拟以可压缩非稳态模式进行。此外，与MOC将管道分为多个分段不同，将每根管道作为单个计算单元处理，显著减少了计算量。与Nault等人（2016年）不同，当前实现还包含了瞬时加速度基（IAB）非稳态摩擦项。由于该项，压力波得以如预期般耗散。得益于EPANET-AWH，当前EPANET中用于稳态建模的算法框架已扩展至包含可压缩瞬态水力分析。因此，面向工程师和研究人员的开源非稳态水力建模软件得以向公众提供，无需使用商业水力建模程序。章节片段代数水锤在维非稳态可压缩流动中，质量和动量必然守恒。质量守恒由连续性方程给出，如公式（1）所示：c2/(gA)·?Q/?x + ?H/?t = 0动量守恒由公式（2）表达：1/(gA)·?Q/?t + ?H/?x + J_Sf + J_Uf = 0J_Sf = fQ|Q|/(2gDA2)J_Uf = (k_UF/gA)·(?Q/?t + c*·sign(Q)·?Q/?x)其中Q为流量（m3/s），H为扬程（m），c为波速（m/s），g为重力加速度（m/s2），A为管道截面积（m2），D为管道直径（m），J_Sf为稳态摩擦单位EPANET-AWH的结构EPANET 3当前水力求解器可以通过采用递归和准稳态的G-GGA方法来分析给水分配系统的水力特性。Ko?ucu等人[16]将用于给水分配系统缓变非稳态流水力分析的RWC-GGA方法应用到了EPANET 3中。虽然RWC-GGA方法能够模拟缓慢变化（如渐进阀门关闭）的水力效应，但无法模拟突发性事件（如突然的工况案例本研究认为，使用两个案例分析该软件的输出以验证EPANET-AWH功能是否正确是必要的。在此范围内，使用EPANET-AWH和Bentley HAMMER对两个不同的给水分配网络进行了非稳态流水力分析。第一个网络受到突然阀门关闭的影响，第二个网络受到泵停机、泵启动和管段处突然抽取50 L/s水量的影响。参考节点扬程的一致性讨论将代数水锤方法集成到EPANET 3中，使得在同一模拟中对给水分配系统的稳态和非稳态流水力进行准自适应分析成为可能。EPANET 3的新水力求解器基于GGA具有稳健的收敛能力，为进行可压缩非稳态流分析提供了有效且有用的框架。对工程师而言，直接结果是涌浪评估不再需要导出到外部工具。结论本研究提出了EPANET-AWH，这是一个将代数水锤方法集成到EPANET 3中的开源扩展。该求解器在两个网络和覆盖阀门关闭、泵停机、泵启动和突然抽取的四种瞬态工况上进行了测试。在事件窗口内，与Bentley HAMMER的吻合度达到NSE 0.90，在较大网络上控制涌浪保护设计的峰值扬程差异小于2 m。在Nault等人（2018年）的1489根管道网络上，符号说明符号描述单位A管道截面积m2A??结合管道特征和非管道水头损失项的对角算子矩阵–A??, A??连接未知扬程节点的管道和非管道连接的关联矩阵–A??A, A??B, A??C管道起始节点、管道末端节点和非管道连接的A??分区–A??连接已知扬程（水库）节点的管道和非管道连接的关联矩阵–A??A, A??B, A??C起始、末端和非管道连接的A??分区–A??对角CRediT作者贡献声明Mehmet Melih Ko?ucu：撰写——审阅与编辑，撰写——原稿，可视化，验证，软件，资源，方法论，调查，形式分析，数据管理，概念化。利益冲突声明作者声明不存在已知的财务利益或个人关系可能对本文工作产生影响的冲突。致谢本研究由伊斯坦布尔理工大学科学研究项目部门在项目编号MAB-2024-45834下提供资助。Mehmet Melih Ko?ucu