CN111339701A - 基于Brunone动态摩阻的管道泄漏特性Godunov模拟方法 - Google Patents

Abstract

本发明公开了一种基于Brunone动态摩阻的管道泄漏特性Godunov模拟方法，首先构建了含Brunone动态摩阻模型的非恒定管流控制方程；然后根据有限体积法Godunov格式建立泄漏管道计算网格系统，并将控制方程转化为黎曼问题；接着采用二阶Godunov格式计算网格交界面通量；随后处理包含泄漏处的边界条件，计算控制体内部流体变量和控制体泄漏边界处流体变量。本发明可以解决MOC格式中非常困难和复杂的动态摩阻模型的二阶精度问题，从而更加准确、高效的模拟泄漏管道水力瞬变特性，确保管道系统的泄漏检测和安全运行。

Description

基于Brunone动态摩阻的管道泄漏特性Godunov模拟方法

技术领域

本发明涉及一种基于Brunone动态摩阻的管道泄漏特性Godunov模拟方法，属于水电站(泵站)管道检测技术领域。

背景技术

管道漏损检测是全球供水行业的热点和难点问题，管道泄漏不仅造成巨额的经济损失，同时可能造成一些次生灾害，如水质受到影响、建筑基础的损坏等，因此寻求准确有效的泄漏检测方法是一个重要研究课题。传统的稳态或准稳态摩阻近似不能在快速瞬变流动中再现实验测量的压力阻尼，因此动态摩阻模型被引入来解释压力耗散，其结果证明更加贴近实验。特征线法(MOC)凭借其简单、准确的优点经常被用来计算管道瞬变流产生的压力。由于公式的复杂性，基于卷积的动态摩阻模型通常采用显式固定网格MOC解的一阶近似，Brunone模型可以很容易地用隐式或显式MOC格式求解。

特征线法在网格固定时，可以很理想的满足库朗数条件，然而实际工程中的管道系统往往会出现不同的波速或长度，所以可能需要改变波速或者网格长度来满足库朗数条件，这就让计算变得繁琐而且可能导致计算结果产生较大的误差。对于简单水锤问题，一阶有限体积法Godunov格式(GTS,Godunov type scheme)格式和固定网格MOC格式给出了相同的结果，但当库朗数小于1时，两者都会导致很强的数值耗散。二阶Godunov格式即使在库朗数明显小于1时也具有很好的鲁棒性和稳定性。

发明内容

针对在管道泄漏水力瞬变数值模拟时，对于实际管道检测中的库朗数小于1带来的数值耗散问题、鲁棒性和稳定性问题以及动态摩阻的问题，本发明提出了一种基于Brunone动态摩阻的管道泄漏特性Godunov模拟方法，可以解决MOC格式中非常困难和复杂的动态摩阻模型的二阶精度问题，从而更加准确、高效的模拟泄漏管道水力瞬变特性，确保管道系统的泄漏检测和安全运行。

为达到上述目的，本发明采用的技术方案如下：

基于Brunone动态摩阻的管道泄漏特性Godunov模拟方法，包括：

将改进的Brunone动态摩阻模型引入非恒定管道控制方程；

在有限体积法体系下将泄漏管道分为N个控制体，将非恒定管道控制方程转化为黎曼问题；

采用二阶Godunov格式求解黎曼问题，得到控制体内部边界通量；

根据控制体内部边界通量计算下一时刻控制体内部流体变量；

根据控制体内部流体变量计算管道两端边界和控制体泄漏边界处流体变量。

进一步的，所述将改进的Brunone动态摩阻模型引入非恒定管道控制方程，包括：

考虑动态摩阻的非恒定管道控制方程为：

其中，H为测压水头，V为平均流速，a为声波在水中传播的速度，g为重力加速度，x为距离，t为时间，D为管径，f_q为准稳态摩擦系数，J_u为动态摩阻项；

动态摩阻项J_u为：

其中，k为摩阻系数。

进一步的，所述摩阻系数确定如下：

采用试验法进行经验预测；

或者根据下式确定：

其中，C^*为Vardy剪切衰减系数，

对于层流：C^*＝0.00476；

对于紊流：

Re为雷诺数。

进一步的，所述在有限体积法体系下建立泄漏管道计算网格，将非恒定管道控制方程转化为黎曼问题，包括：

假定泄漏点发生在第i个控制体的右边界；

将考虑动态摩阻的非恒定管道控制方程写成黎曼问题的形式：

其中，

对于第i个控制体，将

沿控制体左边界i-1/2和右i+1/2积分得到：

其中，Δx为控制体长度，

表示n+1时刻流体变量u在第i个控制体的平均值，

表示n时刻流体变量u在第i个控制体的左边界通量，

表示n时刻流体变量u在第i个控制体的右边界通量。

进一步的，所述采用二阶Godunov格式求解黎曼问题，得到控制体内部边界通量，包括：

在FVM-Godunov格式中，黎曼问题是以下初值问题：

其中，

为n时刻流体变量u在i+1/2界面左控制单元平均值，

为n时刻变量u在i+1/2界面右控制单元平均值，x_i+1/2表示i+1/2边界处的x值，uⁿ(x)为n时刻x位置处流体变量；

通过求解特征方程

得到两个相互独立的特征值：

分别对应以下两个特征向量：

由于特征向量K⁽¹⁾，K⁽²⁾线性独立，所以

和

表示为：

其中，

由此得出系数α₁、α₂、β₁和β₂：

则黎曼问题的通解为：

u(x,t)＝β₁K⁽¹⁾+α₂K⁽²⁾

进而得到流体变量u在i+1/2界面的解析解：

在t∈[tⁿ，tⁿ⁺¹]的计算时间内，流体变量u在第i个控制体的右边界通量表示为：

同理，得到流体变量u在第i个控制体的左边界通量

进一步的，所述根据控制体内部边界通量计算下一时刻控制体内部流体变量，包括：

通过构造分段线性函数来代替单元平均值

其中，

分别为流体变量u在第i个控制体的左、右边界值，Δx为网格单元宽度，

为n时刻流体变量u在第i个控制体内的平均值，

Ω_i是斜率矢量，上角标L表示x→x_i-1/2且x>x_i-1/2，上角标R表示x→x_i+1/2且x＜x_i+1/2；

对于每个网格[x_i-1/2，x_i+1/2]，边界值

跟据分段线性函数按0.5Δt进化得到：

最终得到Riemann问题的近似求解：

进一步的，所述斜率矢量选择MINMOD斜率限制器：

其中，

进一步的，所述根据控制体内部流体变量计算管道两端边界和控制体泄漏边界处流体变量，包括：

在Godunov格式中，存在：

其中，

表示

的特征值；

即：

其中，V_i+1/2表示i+1/2边界处流速，H_i+1/2表示i+1/2边界处测压水头；

V_R，H_R分别表示i+1/2界面右控制单元的流速和水头平均值，取值如下：

V_L，H_L分别表示i+1/2界面左控制单元的流速和水头平均值，取值如下：

管道两端边界满足：

上游恒定水位H_r，H_1/2＝H_r，下游阀门处流速为0，V_N+1/2＝0；

对于上下游边界条件，流体变量分别由上式得到；

控制体泄漏边界处存在：

其中，V_leak和H_leak分别为泄漏处流速和测压水头，C_dA_g为管道泄漏系数，A为管道横截面积；

根据连续性方程：

其中，

和

分别代表泄漏处左、右边界流速；

假定泄漏左右界面压力与泄漏处压力相等，即：

式中，

和

分别代表泄漏处左、右边界测压水头；

联立上述公式，即求得管道泄漏边界处流速和测压水头。

与现有技术相比，本发明提供的基于Brunone动态摩阻模型的管道泄漏特性Godunov模拟方法，具有如下优点：

(1)二阶Godunov格式能解决实际管道泄漏检测中库朗数小于1带来的数值耗散问题，且具有很好的鲁棒性和稳定性；

(2)在二阶Godunov格式中考虑动态摩阻的影响，与实际情况更接近；

(3)考虑Brunone动态摩阻的管道泄漏特性Godunov模拟方法可以很好解决MOC格式中非常困难和复杂的动态摩阻模型的二阶精度问题。

附图说明

图1是本发明方法流程图；

图2是本发明泄漏管道的FVM网格；

图3是本发明管道泄漏处的边界条件图示；

图4是实施例1实验值和考虑Brunone动态摩阻的二阶Godunov格式模拟的管道末端压力的比较；

图5是实施例2泄漏情况下二阶Godunov格式与Lee(2015)MOC模拟的管道末端压力标准化值的比较。

具体实施方式

下面对本发明作进一步描述。以下实施例仅用于更加清楚地说明本发明的技术方案，而不能以此来限制本发明的保护范围。

本发明提供一种基于Brunone动态摩阻的管道泄漏特性Godunov模拟方法，利用改进的Brunone动态摩阻模型的二阶Godunov格式来模拟输水管道泄漏水力瞬变问题，参见图1，具体步骤如下：

1.将改进的Brunone动态摩阻模型引入非恒定管道控制方程；具体为：

步骤1：将改进的Brunone动态摩阻模型引入非恒定管道控制方程；

考虑动态摩阻的非恒定管道控制方程为：

式中，H为测压水头；V为平均流速；a为声波在水中传播的速度；g为重力加速度；x为距离；t为时间；D为管径；f_q为准稳态摩擦系数，由Hagen-Poiseuille定律和Colebrook-White公式确定；J_u为动态摩阻项。

本发明将改进的Brunone动态摩阻模型引入非恒定管道控制方程，具体如下：

式中，k为摩阻系数，可由两种方法确定：一是用试验法进行经验预测；二是根据Vardy剪切衰减系数

其值大小依赖于水流流态和雷诺数的大小。对于层流：C^*＝0.00476；对于紊流：

Re为雷诺数。

2.在FVM(有限体积法)体系下，建立泄漏管道计算网格，将控制方程转化为黎曼问题；具体为：

参见图2，建立泄漏管道FVM体系的计算网格，将泄漏管道分为N个控制体，每个控制体长Δx，控制体中点编号为i(1≤i≤N)，第i个控制体左右边界分别为i-1/2和i+1/2，假定泄漏点发生在第i个控制体的右边界，参见图3。

将步骤1中的控制方程写成黎曼问题的形式：

式中，

流体变量H和V用它们在每个控制体积内的平均值表示，对于第i个控制体，将方程(4)沿控制体边界i-1/2和i+1/2积分得到：

其中，

表示n+1时刻流体变量u在第i个控制体的平均值，

表示n时刻流体变量u在第i个控制体的左边界通量，

表示n时刻流体变量u在第i个控制体的右边界通量。

3.采用二阶Godunov格式计算控制体边界通量；具体为，

(3a)计算控制体之间网格界面的质量和动量通量，包括：

在FVM-Godunov格式中，数值通量是通过求解每个单元界面上的局部Riemann问题来确定的，一般双曲组的Riemann问题是以下初值问题：

式中，

为n时刻流体变量u在i+1/2边界左侧的平均值，

为n时刻变量u在i+1/2边界右侧的平均值，x_i+1/2表示i+1/2边界处的x值。

通过求解特征方程

能够得到两个相互独立的特征值：

分别对应以下两个特征向量：

由于特征向量K⁽¹⁾，K⁽²⁾线性独立，所以

和

可表示为：

可以得出系数α₁、α₂、β₁和β₂：

则

分别是矢量

和

里面的压力和流速，即：

则上述黎曼问题的通解为：

u(x,t)＝β₁K⁽¹⁾+α₂K⁽²⁾ (13)

由式(13)得到控制体界面i+1/2处的解析解：

由此，在t∈[tⁿ，tⁿ⁺¹]的计算时间内，内部边界i+1/2处的通量值可表示为：

同理，求解

内部边界i+1/2指的是内部网格点的i+1/2边界；

控制体界面i+1/2既包含内部边界，也包含实际外部边界。

二阶Godunov格式会在高梯度附近产生虚假振荡，本发明构造了一种非线性方法，通过MUSCL-Hancock格式进行线性插值重构，从而取得具有二阶精度的计算格式，基于Godunov格式引入MUSCL-Hancock格式进行线性重构的具体步骤如下：

第一步：数据重构

通过构造分段线性函数来代替单元平均值

其中，

分别为流体变量u在第i个控制体的左、右边界值，Δx为网格单元宽度，

为n时刻流体变量u在第i个控制体内的平均值。

Ω_i是适当选择的斜率矢量，以在确保结果没有虚假振荡的同时提高方案的精度，为此本发明选择了MINMOD斜率限制器：

式中，上角标L表示x→x_i-1/2且x>x_i-1/2；上角标R表示x→x_i+1/2且x＜x_i+1/2；

第二步：时间进化计算

对于每个网格[x_i-1/2，x_i+1/2]，边界值

跟据公式(18)按0.5Δt进化得到：

第三步：Riemann问题的近似求解：

将公式(21)带入到公式(17)中即可得到具有二阶精度的Godunov格式每个内部计算单元界面i+1/2处的数值通量，进一步可以得出下一时刻每个控制体的流体变量H和V。

4.计算边界单元：管道两端边界及泄漏处边界；

在Godunov格式中：

其中，

表示

的特征值。

即：

其中，V_R，H_R分别表示i+1/2界面右控制单元的流速和水头平均值，取值如下：

V_L，H_L分别表示i+1/2界面左控制单元的流速和水头平均值，取值如下：

上游恒定水位H_r，H_1/2＝H_r，下游阀门处流速为0，V_N+1/2＝0。

所以，对于上下游边界条件，流体变量可分别由式(23)和式(24)得到。

对于管道泄漏处(V_i+1/2,H_i+1/2)：

式中，V_leak和H_leak分别为泄漏处流速和压力水头；C_dA_g为管道泄漏系数；A为管道横截面积。

根据连续性方程：

式中，

和

分别代表了泄漏左右界面流速。

假定泄漏左右界面压力与泄漏处压力相等，即

式中，

和

分别代表了泄漏左右界面压力水头。

联立公式(23)～(27)，可以求得泄漏边界处的流速和水头。

5.处理计算结果，与已发表的实验结果和模拟结果对比验证。

实施例

为了验证并分析本发明考虑Brunone动态摩阻模型的管道泄漏特性Godunov模拟方法的模拟效果，实施例1选取Bergant于2001年设计搭建的管道水力瞬变实验装置系统用于验证本发明Godunov-BM模拟方法的有效性，实施例2选取Lee2015年管道泄漏的模拟结果验证本发明基于Godunov格式的泄漏模型的有效性。实施例1种参数：管道长37.23m，直径0.0221m，上游水箱恒定水头32m，波速1319m/s，管道流速V₀＝0.3m/s，雷诺数Re＝5600，阀门瞬间关闭。实施例2中参数：管道长2000m，波速1000m/s，管道直径为0.5m，上游恒定水头50m，初始流量Q₀＝0.02m³/s，雷诺数为50000，恒定摩阻系数f＝0.02，不考虑动态摩阻的影响，泄漏孔口在距上游水库300m处，初始泄漏量为0.2Q₀。

采用本发明方法编程计算，将考虑Brunone动态摩阻模型的二阶Godunov格式计算结果与实施例1实验结果对比，压力曲线如图4所示；将考虑管道泄漏的二阶Godunov格式计算结果与实施例2Lee2015年模拟结果对比，压力曲线如图5所示。由图4可知，与不考虑动态摩阻的模型相比，动态摩阻模型能精准描述瞬变管道压力幅值和衰减，随着周期数增加仍能较好的吻合，而传统准恒定摩阻模型模拟的压力波在两个周期后误差逐渐增大，相位偏移加大。由图5可知，管道在发生泄漏时末端压力会出现局部的畸变点，且畸变量与时间呈正向关系，二阶Godunov格式模拟泄漏管道末端压力结果与MOC方法所得结果一致。因此，考虑Brunone动态摩阻的管道泄漏特性Godunov模拟方法可以很好解决管道泄漏模拟过程中MOC格式非常困难和复杂的动态摩阻模型的二阶精度问题，能更好地再现实验压力振荡和准确预测泄漏管道瞬变的压力阻尼。

以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明技术原理的前提下，还可以做出若干改进和变形，这些改进和变形也应视为本发明的保护范围。

Claims (8)

1.基于Brunone动态摩阻的管道泄漏特性Godunov模拟方法，其特征在于，包括：

将改进的Brunone动态摩阻模型引入非恒定管道控制方程；

在有限体积法体系下将泄漏管道分为N个控制体，将非恒定管道控制方程转化为黎曼问题；

采用二阶Godunov格式求解黎曼问题，得到控制体内部边界通量；

根据控制体内部边界通量计算下一时刻控制体内部流体变量；

根据控制体内部流体变量计算管道两端边界和控制体泄漏边界处流体变量。

2.根据权利要求1所述的基于Brunone动态摩阻的管道泄漏特性Godunov模拟方法，其特征在于，所述将改进的Brunone动态摩阻模型引入非恒定管道控制方程，包括：

考虑动态摩阻的非恒定管道控制方程为：

其中，H为测压水头，V为平均流速，a为声波在水中传播的速度，g为重力加速度，x为距离，t为时间，D为管径，f_q为准稳态摩擦系数，J_u为动态摩阻项；

动态摩阻项J_u为：

其中，k为摩阻系数。

3.根据权利要求2所述的基于Brunone动态摩阻的管道泄漏特性Godunov模拟方法，其特征在于，所述摩阻系数确定如下：

采用试验法进行经验预测；

或者根据下式确定：

其中，C^*为Vardy剪切衰减系数，

对于层流：C^*＝0.00476；

对于紊流：

Re为雷诺数。

4.根据权利要求2所述的基于Brunone动态摩阻的管道泄漏特性Godunov模拟方法，其特征在于，所述在有限体积法体系下建立泄漏管道计算网格，将非恒定管道控制方程转化为黎曼问题，包括：

假定泄漏点发生在第i个控制体的右边界；

将考虑动态摩阻的非恒定管道控制方程写成黎曼问题的形式：

其中，

对于第i个控制体，将

沿控制体左边界i-1/2和右i+1/2积分得到：

其中，Δx为控制体长度，

表示n+1时刻流体变量u在第i个控制体的平均值，

表示n时刻流体变量u在第i个控制体的左边界通量，

表示n时刻流体变量u在第i个控制体的右边界通量。

5.根据权利要求4所述的基于Brunone动态摩阻的管道泄漏特性Godunov模拟方法，其特征在于，所述采用二阶Godunov格式求解黎曼问题，得到控制体内部边界通量，包括：

在FVM-Godunov格式中，黎曼问题是以下初值问题：

其中，

为n时刻流体变量u在i+1/2界面左控制单元平均值，

为n时刻变量u在i+1/2界面右控制单元平均值，x_i+1/2表示i+1/2边界处的x值，uⁿ(x)为n时刻x位置处流体变量；

通过求解特征方程

得到两个相互独立的特征值：

分别对应以下两个特征向量：

由于特征向量K⁽¹⁾，K⁽²⁾线性独立，所以

和

表示为：

其中，

由此得出系数α₁、α₂、β₁和β₂：

则黎曼问题的通解为：

u(x,t)＝β₁K⁽¹⁾+α₂K⁽²⁾

进而得到流体变量u在i+1/2界面的解析解：

在t∈[tⁿ，tⁿ⁺¹]的计算时间内，流体变量u在第i个控制体的右边界通量表示为：

同理，得到流体变量u在第i个控制体的左边界通量

6.根据权利要求5所述的基于Brunone动态摩阻的管道泄漏特性Godunov模拟方法，其特征在于，所述根据控制体内部边界通量计算下一时刻控制体内部流体变量，包括：

通过构造分段线性函数来代替单元平均值

其中，

分别为流体变量u在第i个控制体的左、右边界值，Δx为网格单元宽度，

为n时刻流体变量u在第i个控制体内的平均值，

Ω_i是斜率矢量，上角标L表示x→x_i-1/2且x>x_i-1/2，上角标R表示x→x_i+1/2且x＜x_i+1/2；

对于每个网格[x_i-1/2，x_i+1/2]，边界值

跟据分段线性函数按0.5Δt进化得到：

最终得到Riemann问题的近似求解：

7.根据权利要求6所述的基于Brunone动态摩阻的管道泄漏特性Godunov模拟方法，其特征在于，所述斜率矢量选择MINMOD斜率限制器：

其中，

8.根据权利要求5所述的基于Brunone动态摩阻的管道泄漏特性Godunov模拟方法，其特征在于，所述根据控制体内部流体变量计算管道两端边界和控制体泄漏边界处流体变量，包括：

在Godunov格式中，存在：

其中，

表示

的特征值；

即：

其中，V_i+1/2表示i+1/2边界处流速，H_i+1/2表示i+1/2边界处测压水头；

V_R，H_R分别表示i+1/2界面右控制单元的流速和水头平均值，取值如下：

V_L，H_L分别表示i+1/2界面左控制单元的流速和水头平均值，取值如下：

管道两端边界满足：

上游恒定水位H_r，H_1/2＝H_r，下游阀门处流速为0，V_N+1/2＝0；

对于上下游边界条件，流体变量分别由上式得到；

控制体泄漏边界处存在：

其中，V_leak和H_leak分别为泄漏处流速和测压水头，C_dA_g为管道泄漏系数，A为管道横截面积；

根据连续性方程：

其中，

和

分别代表泄漏处左、右边界流速；

假定泄漏左右界面压力与泄漏处压力相等，即：

式中，

和

分别代表泄漏处左、右边界测压水头；

联立上述公式，即求得管道泄漏边界处流速和测压水头。

Images