基于Godunov 格式的输水管道阻塞数值模拟与识别

胡朝仲，杨旭，李赟杰，白光国，周领，4

（1.云南水利水电职业学院水利工程学院，云南昆明 650499；2.中水北方勘测设计研究有限责任公司，天津 300222；3.河海大学水利水电学院，江苏南京 210098；4.河海大学长江保护与绿色发展研究院，江苏南京 210098）

0 引言

引调水工程和城市供水管网建设的发展，大大满足了生活和生产用水，同时也暴露出了一些问题，如管道泄漏、阻塞等。在实际输水管道中，生物膜堆积、腐蚀和沉积过程会形成堵塞，且堵塞的大小和严重程度通常随着时间的推移而增加［1］。管道阻塞会导致额外的能量损失和水质问题，因此管道阻塞的准确模拟和精准检测十分重要。

以往的管道瞬变流数值模拟大都基于特征线法求解，该方法简单、高效，常用于各类瞬变流现象模拟［2，3］。然而，特征线法在求解复杂管道模型时需要插值近似或调整波速，容易产生数值耗散或色散。基于此，Zhao等［4］提出一种有限体积法Godunov 求解格式，该求解格式受库朗特数的影响较小。Zhou等［5，6］基于该求解格式，纳入了不同类型的动态摩阻，结果表明Brunone 摩阻模型结合Godunov 求解格式在低库朗特数工况条件下数值耗散小，能更好地计算复杂管道模型。管道阻塞分为局部阻塞和管段阻塞，局部阻塞即尺度小于最小波长，可以视为一个不连续点，如一个未完全关闭的阀门；管段阻塞的尺度大于最小波长，可以视为一段直径较小的管段串联在管道中［7］。对于管道阻塞，特别是管段阻塞，在数值模拟时可能涉及到库朗特数小于1的情况，此时特征线法求解格式会有较大误差。

本文提出了一种基于Godunov 求解格式的管道阻塞模型。首先，给出含Brunone 动态摩阻的控制方程和边界条件，以此建立相应的管道阻塞模型。然后，分别对局部阻塞和管段阻塞进行模型验证，并根据波传播理论分析不同阻塞工况压力波动曲线突变原因，并给出阻塞参数识别方法。通过参数敏感性分析研究了局部阻塞和管段阻塞的影响因素和结果。

1 数学模型及求解

1.1 控制方程

管道瞬变流的基本控制方程为［2，8，9］：

式中：H为测压管水头；V为管道断面的平均流速；a为波速；g为重力加速度；x为沿管道轴向长度；t为时间；D为管道内径；ρ为水体密度；f为达西-维斯巴赫系数；SIGN（V）为符号函数；k3为Brunone摩阻系数［10］。k3通常采用Vardy等［11］所给的系数公式：

1.2 边界条件

对于边界条件，通过在边界两端各增加两个虚拟控制体，可以得到对边界控制体和内部控制体的统一计算模拟，并达到二阶精度要求［12］。结合黎曼问题在边界处的解析解、上下游边界处的动量和质量方程可得到边界处的未知量。

对于上游恒定水库水位HR，流速可以表示为：

对于下游阀门-水库边界条件，有：

对于管道中阀门或串联节点边界条件，主要考虑连接边界处的水头损失，即：

1.3 Godunov求解格式

对于有限体积法Godunov 求解格式，首先将瞬变流基本方程以黎曼问题的求解格式表示为矩阵形式［4-6］：

式中：u为求解向量；f（u）为控制体通量矢量；Aˉ为矢量f中各变量对u中各求解变量的偏导数矩阵；s（u）为源项，代表了管道的摩阻项。

图1 沿管道轴向的计算网格

式中：fi-1/2表示i-1/2界面的数值通量；fi+1/2表示i+1/2界面的数值通量。因此，在计算下一时步的向量U时，需涉及当前时步的控制体两边界数值通量和源项。

针对Godunov 格式，控制体边界数值通量可由局部黎曼问题得到，可描述为初值问题：

对于和，采用分段多项式重构近似，其精度取决于多项式的形式。对于二阶Godunov 近似，通过引入MUSCLHancock 格式进行线性插值重构。因为二阶格式会在高梯度附近产生假振荡，采用MINMOD 斜率限制器削减这种数值震荡［6，13］。因此，可以得到控制体积各界面的通量向量。源项采用显式的二阶龙格-库塔离散化方法，考虑到Courant-Friedrichs-Lewy（CFL）定律和龙格库塔离散稳定性条件，可得到源项的允许时间步长［6］。

2 模型验证

针对局部阻塞、管段阻塞两种堵塞类型对本文所建模型进行验证、讨论和分析。

2.1 局部阻塞模型验证

将模拟结果与Lee等［14］2015 年中的论文数据作对比，该系统上游为水库，下游阀门，管道中间某处有一未完全关闭阀门。具体参数为：管道长度L=2 000 m，上游恒定水位HR=50 m，中间阀门位于距上游水库0.3L处，稳态损失系数K=1 890，波速a=1 000 m/s，管道直径D=0.5 m，初始流量Q0=0.02 m3/s，稳态摩阻系数f=0.02。阀门瞬间关闭，模型模拟结果对比如图2 所示。图2中，对纵坐标压力进行标准化处理，即H*=（H-H0）/（aV0/g），其中H*为无量纲水头，H0为稳态压力，V0为稳态流速。

由图2 可知，本文提出的基于Godunov 求解格式的管道局部阻塞模型与Lee 等2015 年模型模拟结果基本重合，验证了局部阻塞模型的正确性。此外，根据模拟结果可知，对于管道局部阻塞，其压力波动曲线在第一峰值处会突然增加，这是由于压力波从末端阀门传播到阻塞处时会发生反射，反射的正波与末端阀门突然关闭产生的主波叠加，导致反射波传播到阀门末端时压力骤增，由此可以根据压力骤增的时间判断局部阻塞发生位置。在本算例中，第一峰值压力骤增的时刻为t=2.8 s，根据波速a=1 000 m/s，可以得到阻塞点距下由阀门位置xl=at/2=1 400 m，与实际阻塞点位置一致。

图2 局部阻塞模型验证对比图

2.2 管段阻塞模型验证

将模拟结果与Duan等［15］实验数据作对比，该系统上游为水库，下游阀门，管道分为三段，中间段连接直径较小的管道。具体参数为：管道总长度L=41.5 m，其中L1=8.8 m，L2=12.1 m，L3=20.6 m；上游恒定水位HR=38.2 m，波速a=1 180 m/s，管道分别为直径D1=0.072 4 m，D2=0.022 25 m，D3=0.072 4 m，初始流量Q0=9.7×10-5m3/s，稳态摩阻系数f=0.045 9，阀门关闭时间为0.01 s。模型与实验对比结果如图3 所示，图3（a）和（b）分别为稳态摩阻与动态摩阻模型模拟结果。

由图3 可知，提出的管段阻塞模型能较好地复现压力波动曲线周期和衰减。图3（a）所示，仅考虑稳态摩阻时，在第一峰值处与实验结果基本重合，但随着时间变化，最大与最小压力值和实验结果偏差越大。图3（b）所示，考虑动态摩阻模型，随着时间变化仍能较好的模拟实验最大最小压力值，仅会发生微小的相位偏移。与局部阻塞一致，末端阀门产生的压力波在管段阻塞处也会发生反射。当管径突然变小时，反射的压力波为正波，管径突然变大时，反射的压力波为负波。当正波或负波传播到末端阀门时，会产生相应的压力增加或减小，因此在本结果前半个周期中会出现压力波先增加后减小的现象。同理，据此可以计算管段阻塞的位置和长度。以本工况为例，前半个周期内压力波骤增和骤减的时间分别为t1=0.035 s，t2=0.055 s，根据波速a=1 180 m/s，可以得到阻塞管段末端距下游阀门20.65 m，阻塞长度为11.8 m，这与实际阻塞位置和长度基本一致。

图3 管段阻塞模型验证对比图

3 参数敏感性分析

3.1 局部阻塞

采用2.1 节局部阻塞工况分别对管道初始流量Q0、局部阻塞损失系数K和局部阻塞位置xl进行参数敏感性分析，结果如图4所示。

对管道初始流量Q0参数敏感性分析，分别采用0.02、0.04和0.08 m3/s，结果如图4（a）所示。图4（a）表明，阻塞大小不变的情况下，随着初始流量的增加，在第一峰值由阻塞引起的反射压力幅值随之增大，即反射系数随初始流量增大而增大。Yan［16］理论推导了管道局部阻塞引起的反射系数解析解，反射系数r可以表示为：

式中：A为管道横截面积。

由解析解可知：其他条件不变的情况下，随着管道初始流量的增大，反射系数也会增大，且增大的幅度与流量呈正比。这也证明了模拟的反射系数与解析解结果一致。此外，随着流量的增加，压力曲线的峰谷压力也会增大，衰减速度更快。

对局部阻塞损失系数K进行参数敏感性分析，K值分别取1 890、3 000 和5 000，结果如图4（b）所示。所得规律与管道初始流量类似，具体表现为：在其他条件不变的情况下，随着K值的增大（即局部阻塞程度增加），峰值反射增大，且反射系数与K呈正比（与解析解一致）；峰谷压力增大，衰减速度加快。

对局部阻塞位置xl进行参数敏感性分析，局部阻塞位置xl分别距上游水库0.3L、0.5L和0.7L，结果如图4（c）所示。由图4（c）可知，不同局部阻塞位置在前半个周期内表现为压力上升时间不同，即发生反射的时间不同。因此，可根据压力波反射的时间估计局部阻塞位置，具体见3.1节。

图4 局部阻塞参数敏感性分析

3.2 管段阻塞

对管段阻塞进行参数敏感性分析，系统上游为水库，下游阀门，管道分为三段，中间段连接直径较小的管道。具体参数为：管道总长度L=300 m，其中L1、L2、L3均为100 m；上游恒定水位HR=50 m，波速a=1 000 m/s，管道直径分别为D1=0.2 m，D2=0.14 m，D3=0.2 m，初始流量Q0=0.01 m3/s，稳态摩阻系数f=0.02，阀门瞬间关闭。分别对管段阻塞程度和阻塞长度进行参数敏感性分析，结果如图5所示。

对管段阻塞程度进行参数敏感性分析，保持其他参数不变，中间管道直径D2分别选取0.14、0.16、0.18和0.2 m 进行计算分析，结果如图5（a）所示。图5（a）表明，随着阻塞管段直径的减小，阻塞反射幅值增加，表现为第一压力峰值增大，且系统最大峰值压力不再为第一峰值压力。对于随时间变化的压力波动曲线，阻塞管段直径减小使得系统压力波动周期增大，这是由于布拉格共振效应引起的［17，18］。

对管段阻塞长度进行参数敏感性分析，保持其他参数不变，第三段管道长度恒定100 m，第一段管道长度和中间管道长度之和保持200 m 不变，中间管道长度L2分别选取20、60 和100 m 进行计算分析，结果如图5（b）所示。图5（b）表明，随着管段阻塞长度增加，第一峰值压力下降时间延长，这是因为瞬态波穿过阻塞管段时间变短。波动周期随管段阻塞长度增加而逐渐增大，与阻塞管段直径变化引起的相移一致，这也是由布拉格共振效应引起的。

图5 管段阻塞参数敏感性分析

4 结论

本文首先提出了基于Godunov 求解格式的管道阻塞模型，并针对局部阻塞和管段阻塞分别进行模型验证和结果分析，主要结论有：①本文提出的基于Godunov 求解格式的管道阻塞模型精确度高，能较好的模拟管道阻塞压力轨迹；②根据压力波传播理论，可根据压力的骤增或骤减计算管道中局部阻塞的位置以及管段阻塞的位置和长度；③局部阻塞压力峰值大小受管道初始流量、局部阻塞损失系数的影响，压力峰值出现的时间受局部阻塞位置的影响；④管段阻塞程度影响压力峰值和波动周期大小，管段阻塞长度影响压力峰值出现的时间和波动周期。

同时应当注意到，在计算阻塞位置或长度时，压力变化点的时间不精确会导致误差，利用一些信号处理技术如小波变换等可能会提高计算精度。