CN113656926B - 基于Schohl卷积近似的管道瞬变流模拟方法 - Google Patents

Abstract

本发明公开了一种基于Schohl卷积近似的管道瞬变流模拟方法，包括：构建含Zielke非恒定摩阻项和粘弹性项的管道***控制方程；将控制方程以黎曼问题的求解格式表示为矩阵形式，将非恒定摩阻项和粘弹性项放入源项；建立有限体积法求解体系下的计算网格，通过二阶Godunov格式计算网格边界通量；对源项中的Zielke非恒定摩阻项和粘弹性项进行Schohl卷积近似，将源项合并到矩阵方程中，并进行稳定性约束；结合实际工程给出计算初始条件，得出计算结果。本发明结合有限体积法Godunov格式，对控制方程中的非恒定摩阻项和粘弹性项进行Schohl卷积近似，可解决现有一维瞬变流模拟方法依赖库朗特数，涉及到非恒定摩阻和粘弹性项计算时间长、求解复杂等问题。

Description

基于Schohl卷积近似的管道瞬变流模拟方法

技术领域

本发明属于计算水力学领域，具体涉及一种基于Schohl卷积近似的管道瞬变流非恒定摩阻和粘弹性模拟方法。

背景技术

输水管网***由于阀门突然地启闭，会产生瞬时高压流动，称为瞬变流。如果压力过高则会造成管道***的泄漏或损坏，甚至引起爆管，严重威胁着饮用水安全供给。因此，探明管道***液体流动状态，对保障输水管网安全稳定运行至关重要。

管道***的一维瞬变流模拟是探明管道***流动状态的有效手段。以往的一维瞬变流模型中多使用Zielke原始加权函数非恒定摩阻模型，结合广义K-V粘弹性模型，利用特征线法求解瞬变管道压力和流量变化。该方法存在以下缺陷：1)Zielke原始加权函数模型需要存储历史时刻流速变化，导致计算时间较长；2)原始K-V粘弹性模型求解方法复杂，且使用一阶近似；3)特征线法对于单管较为方便，但在复杂管道模拟时需要插值近似，容易产生数值耗散。

因此，亟需一种精确、高效、简便的输水管网瞬变流模拟方法，对于输水管道***安全运行、管道状态实时监测具有重要的理论意义和实际应用价值。

发明内容

发明目的：为了克服现有技术的一维瞬变流模型中存在的计算时间长、求解复杂、通用性差以及难以实时监测等问题，提供一种基于Schohl卷积近似的管道瞬变流非恒定摩阻和粘弹性非恒定摩阻和粘弹性模拟方法。

技术方案：为实现上述目的，本发明提供一种基于Schohl卷积近似的管道瞬变流非恒定摩阻和粘弹性模拟方法，包括如下步骤：

S1：构建含Zielke非恒定摩阻项和粘弹性项的管道***控制方程；

S2：将控制方程以黎曼问题的求解格式表示为矩阵形式，将非恒定摩阻项和粘弹性项放入源项；

S3：建立有限体积法求解体系下的计算网格，通过二阶Godunov格式计算网格边界通量；

S4：对源项中的Zielke非恒定摩阻项和粘弹性项进行Schohl卷积近似，将源项合并到矩阵方程中，并进行稳定性约束；

S5：结合实际工程给出计算初始条件，利用Fortran编程得出计算结果。

进一步地，所述步骤S1中管道***控制方程的方程式为：

式中，H为管道某截面处压力；V为管道某截面的平均流速；ε_r为管道迟滞应变；a为有压管道水中波速；g为重力加速度；x为沿管道轴线的距离；t为时间；D为管道内径；ρ为流体密度；f为达西-维斯巴赫系数；v为水的动力粘度；u为微分变量；W为历史流速变化的加权函数。

进一步地，所述步骤S2具体为：

其中，

式中，u为求解向量，f(u)为控制体通量矢量，为矢量f中各变量对u中各求解变量的偏导数矩阵，s(u)为源项，分别代表了管道的粘弹性项和摩阻项(显式)，若s(u)＝0，即不考虑管道粘弹性和摩阻项，则该***为常系数线性齐次***。

进一步地，所述步骤S3中计算网格的建立方法为：

在有限体积法求解体系下，将管段分为N个控制体，建立计算网格，每个控制体长Δx，计算时间步长为Δt，对于第i(1≤i≤N)个控制体，其左右边界分别记为i-1/2和i+1/2；为求解该黎曼问题，将矩阵方程沿x方向在第i个控制体两边界界面之间积分得到以下公式：

式中，f_i-1/2表示i-1/2界面的数值通量；f_i+1/2表示i+1/2界面的数值通量；向量u可由控制体i内的平均值表示，故上式可整理为：

式中，上标n表示t时刻，即当前计算时步；上标n+1表示t+Δt时刻，即下一计算时步。因此，在计算下一时步的向量U时，需涉及当前时步的控制体两边界数值通量和源项。

进一步地，所述步骤S3中网格边界通量的计算方法为：

针对有限体积法Godunov求解格式，控制体边界的数值通量可由界面处的局部黎曼问题得到：对于一个常系数线性齐次双曲***，黎曼问题可描述为下列初值问题：

式中，第n个计算时步中向量u在i+1/2界面左侧的平均值，/>和/>是向量/>的向量元素；/>第n个计算时步中向量u在i+1/2界面右侧的平均值，/>和是向量/>的向量元素；x_i+1/2为第i个控制体右侧边界沿管道轴线的距离；

通过求解特征方程能够得到两个相互独立的特征值：

式中，I为单位矩阵，为要求解的特征值，j＝1～2；根据Rankine-Hugoniot条件 在t∈[t，t+Δt]的计算时间内，内部边界i+1/2处的通量值/>可表示为：

其中，

对于和/>采用分段多项式重构近似，其精度取决于多项式的形式。

进一步地，所述步骤S3中对于和/>采用二阶Godunov近似，通过引入MUSCL-Hancock格式进行线性插值重构，具体的步骤如下：

A1：数据重构：

通过构造分段线性函数U_i(x)来获取第i个控制体左右边界极值和/>

其中，Ω_i为在第i个计算单元的斜率向量，适当选取Ω_i的值可以在确保结果没有虚假振荡的同时提高方案的精度。本文选择MINMOD(流速限制函数)斜率限制器：

式中，

A2：时间进化计算：

对于每个控制体[x_i-1/2，x_i+1/2]，左右边界极值按0.5Δt时间进化得到：

式中，分别为时间推进后的左右边界极值，/>和/>分别为左右边界极值/>处的通量，

A3：黎曼问题的近似求解：

对于边界条件，通过在两端边界各增加两个虚拟控制体-1、0和N+1、N+2，可以得到对边界控制体和内部控制体的统一计算模拟，并达到二阶精度要求；

根据黎曼问题在边界处的解析解：

可推得与正负特征线相关的黎曼不变量方程：

上游水库边界，有

且代入上式得

其虚拟单元满足：

下游边界条件，有

式中，C_V＝(Q₀τ)²/2ΔH₀，Q₀和ΔH₀分别为初始时刻阀门处的流量和水头损失，τ为阀门开度；H_RD为下游水库水头；当C_P-H_RD＞0时，S＝+1，当C_P-H_RD＜0时，S＝-1；

对于下游阀门+水库边界

根据末端阀门孔口方程，可知

带入阀门正特征线黎曼不变量方程，得

其虚拟单元满足：

进一步地，所述步骤S4中对源项中的Zielke非恒定摩阻项和粘弹性项进行Schohl卷积近似的具体方法为：

对于Zielke非恒定摩阻项，在公式(2)中，令

式中，

其中，W_app(τ)是W(τ)的近似函数，由W_i(τ)的和函数组成，i＝1～5；m_i和n_i分别是加权函数的系数，对于Schohl卷积近似摩阻，m_i＝(1.051；2.358；9.021；29.47；79.55)；n_i＝(26.65；100；669.6；6497；57990)。令

则有，

故

对于粘弹性项，有

式中，ε_r(t)为ε_r的时间函数；α为泊松比；e为管道壁厚；J(t*)为蠕变柔量，t*为微分变量；对于广义K-V模型，蠕变柔量函数可以表示为

式中，J_k和τ_k均为粘弹性力学模型参数，可根据蠕变实验数据校核得到；M为校核得到的J_k的个数；对J(t*)求偏导数可得

根据卷积的微分性质，有

令

即

根据卷积的交换性可得

令

则

根据Schohl卷积近似，对于t+Δt时刻，有

故

式中，

CO_A3＝CO_A1·CO_A2 (45)

CO_A4＝CO_A1(1-CO_A2) (46)。

进一步地，所述步骤S4中将源项合并到矩阵方程中，并进行稳定性约束的具体方法为：

在步骤S3中未考虑源项的影响，可采用时间***法将源项引入非线性***求解，求解精度取决于对源项的积分方式；对于二阶精度，采用显式二阶Runge-Kutta离散化方法求解源项s(u)及其空间积分，通过两次更新源项，得到：

式中，

稳定性约束主要包括对流项的CFL(Courant-Friedrich-Lewy)条件以及源项两次更新时的约束；对于CFL条件，有

式中，C_r为库朗数，Δt_max，CFL为CFL条件下的最大时间步长；

利用二阶Runge-Kutta格式引入源项，约束条件为

因此可以得到

因此，计算所允许的最大时间步长需满足

Δt≤min(Δt_max，CFL，Δt_max，s) (54)

有益效果：本发明与现有技术相比，具备如下优点：

1、对Zielke原始加权函数模型进行Schohl卷积近似，降低了计算时长。

2、本发明提供的模型提出了新的基于Schohl卷积近似的广义K-V模型求解方法，简化了求解难度，提升了计算精度。

3、使用二阶Godunov格式求解控制方程，降低了复杂管道模拟时的数值耗散程度。

附图说明

图1是本发明的实现流程图；

图2是沿管道轴向的计算网格示意图；

图3是本发明反法模拟结果与实验对比图。

具体实施方式

下面结合附图和具体实施例，进一步阐明本发明，应理解这些实施例仅用于说明本发明而不用于限制本发明的范围，在阅读了本发明之后，本领域技术人员对本发明的各种等价形式的修改均落于本申请所附权利要求所限定的范围。

如图1所示，本发明提供一种基于Schohl卷积近似的管道瞬变流非恒定摩阻和粘弹性模拟方法，包括如下步骤：

S1：构建含Zielke非恒定摩阻项和粘弹性项的管道***控制方程：

式中，H为管道某截面处压力；V为管道某截面的平均流速；ε_r为管道迟滞应变；a为有压管道水中波速；g为重力加速度；x为沿管道轴线的距离；t为时间；D为管道内径；ρ为流体密度；f为达西-维斯巴赫系数；ν为水的动力粘度；u为微分变量；W为历史流速变化的加权函数。

S2：将控制方程以黎曼问题的求解格式表示为矩阵形式，将非恒定摩阻项和粘弹性项放入源项：

其中，

式中，u为求解向量，f(u)为控制体通量矢量，为矢量f中各变量对u中各求解变量的偏导数矩阵，s(u)为源项，分别代表了管道的粘弹性项和摩阻项(显式)，若s(u)＝0，即不考虑管道粘弹性和摩阻项，则该***为常系数线性齐次***。

S3：一、建立有限体积法求解体系下的计算网格：

如图2所示，在有限体积法求解体系下，将管段分为N个控制体，建立计算网格，每个控制体长Δx，计算时间步长为Δt，对于第i(1≤i≤N)个控制体，其左右边界分别记为i-1/2和i+1/2；为求解该黎曼问题，将矩阵方程沿x方向在第i个控制体两边界界面之间积分得到以下公式：

式中，f_i-1/2表示i-1/2界面的数值通量；f_i+l/2表示i+1/2界面的数值通量；向量u可由控制体i内的平均值表示，故上式可整理为：

式中，上标n表示t时刻，即当前计算时步；上标n+1表示t+Δt时刻，即下一计算时步。因此，在计算下一时步的向量U时，需涉及当前时步的控制体两边界数值通量和源项。

二、通过二阶Godunov格式计算网格边界通量：

针对有限体积法Godunov求解格式，控制体边界的数值通量可由界面处的局部黎曼问题得到：对于一个常系数线性齐次双曲***，黎曼问题可描述为下列初值问题：

式中，第n个计算时步中向量u在i+1/2界面左侧的平均值，/>和/>是向量/>的向量元素；/>第n个计算时步中向量u在i+1/2界面右侧的平均值，/>和/>是向量/>的向量元素；x_i+1/2为第i个控制体右侧边界沿管道轴线的距离；

通过求解特征方程能够得到两个相互独立的特征值：

式中，I为单位矩阵，为要求解的特征值，j＝1～2；根据Rankine-Hugoniot条件 在t∈[t，t+Δt]的计算时间内，内部边界i+1/2处的通量值/>可表示为：

其中，

对于和/>采用分段多项式重构近似，其精度取决于多项式的形式。

对于和/>采用二阶Godunov近似，通过引入MUSCL-Hancock格式进行线性插值重构，具体的步骤如下：

A1：数据重构：

通过构造分段线性函数u_i(x)来获取第i个控制体左右边界极值和/>

其中，Ω_i为在第_i个计算单元的斜率向量，适当选取Ω_i的值可以在确保结果没有虚假振荡的同时提高方案的精度。本文选择MINMOD(流速限制函数)斜率限制器：

式中，

A2：时间进化计算：

对于每个控制体[x_i-1/2，x_i+1/2]，左右边界极值按0.5Δt时间进化得到：

式中，分别为时间推进后的左右边界极值，/>和/>分别为左右边界极值/>处的通量，

A3：黎曼问题的近似求解：

对于边界条件，通过在两端边界各增加两个虚拟控制体-1、0和N+1、N+2，可以得到对边界控制体和内部控制体的统一计算模拟，并达到二阶精度要求；

根据黎曼问题在边界处的解析解：

可推得与正负特征线相关的黎曼不变量方程：

上游水库边界，有

且代入上式得

其虚拟单元满足：

下游边界条件，有

式中，C_V＝(Q₀τ)²/2ΔH₀，Q₀和ΔH₀分别为初始时刻阀门处的流量和水头损失，τ为阀门开度；H_RD为下游水库水头；当C_P-H_RD>0时，S＝+1，当C_P-H_RD<0时，S＝-1；

对于下游阀门+水库边界

根据末端阀门孔口方程，可知

带入阀门正特征线黎曼不变量方程，得

其虚拟单元满足：

S4：一、对源项中的Zielke非恒定摩阻项和粘弹性项进行Schohl卷积近似：

对于Zielke非恒定摩阻项，在公式(2)中，令

式中，

其中，W_app(τ)是W(τ)的近似函数，由W_i(τ)的和函数组成，i＝1～5；m_i和n_i分别是加权函数的系数，对于Schohl卷积近似摩阻，m_i＝(1.051；2.358；9.021；29.47；79.55)；n_i＝(26.65；100；669.6；6497；57990)。令

则有，

故

对于粘弹性项，有

式中，ε_r(t)为ε_r的时间函数；α为泊松比；e为管道壁厚；J(t*)为蠕变柔量，t*为微分变量；对于广义K-V模型，蠕变柔量函数可以表示为

式中，J_k和τ_k均为粘弹性力学模型参数，可根据蠕变实验数据校核得到；M为校核得到的J_k的个数；对J(t*)求偏导数可得

根据卷积的微分性质，有

令

即

根据卷积的交换性可得

令

则

根据Schohl卷积近似，对于t+Δt时刻，有

故

式中，

CO_A3＝CO_A1·CO_A2 (45)

CO_A4＝CO_A1(1-CO_A2) (46)

二、将源项合并到矩阵方程中，并进行稳定性约束：

在步骤S3中未考虑源项的影响，可采用时间***法将源项引入非线性***求解，求解精度取决于对源项的积分方式；对于二阶精度，采用显式二阶Runge-Kutta离散化方法求解源项s(u)及其空间积分，通过两次更新源项，得到：

式中，

/>

稳定性约束主要包括对流项的CFL(Courant-Friedrich-Lewy)条件以及源项两次更新时的约束；对于CFL条件，有

式中，C_r为库朗数，Δt_max，CFL为CFL条件下的最大时间步长；

利用二阶Runge-Kutta格式引入源项，约束条件为

因此可以得到

因此，计算所允许的最大时间步长需满足

Δt≤min(Δt_max，CFL，Δt_max，s) (54)

S5：结合实际工程给出计算初始条件，利用Fortran编程得出计算结果：

带入实际工况中的初始条件，包括上游水库或水箱水头Hr，管道直径D，根据管道长度所划分的控制体的个数N，以及每段控制体的长度Δx，管道波速a，雷诺数Re等，可求得管道中各个控制体的测压管水头H和流速V，从而实现了管道水力瞬变的模拟。

基于上述方案，本实施例中为了验证本发明方法的效果，选取Covas于2003年设计搭建的粘弹性管道水力瞬变实验装置***用于验证本发明方法的有效性，实验和模拟对比结果如图3所示，结果表明本发明所提新的方法能够很好的模拟管道水力瞬变的压力峰值、衰减和相位。

Claims (6)

1.一种基于Schohl卷积近似的管道瞬变流模拟方法，其特征在于，包括如下步骤：

S1：构建含Zielke非恒定摩阻项和粘弹性项的管道***控制方程；

S2：将控制方程以黎曼问题的求解格式表示为矩阵形式，将非恒定摩阻项和粘弹性项放入源项；

S3：建立有限体积法求解体系下的计算网格，通过二阶Godunov格式计算网格边界通量；

S4：对源项中的Zielke非恒定摩阻项和粘弹性项进行Schohl卷积近似，将源项合并到矩阵方程中，并进行稳定性约束；

S5：结合实际工程给出计算初始条件，得出计算结果；

所述步骤S4中对源项中的Zielke非恒定摩阻项和粘弹性项进行Schohl卷积近似的具体方法为：

对于Zielke非恒定摩阻项，在公式(2)中，令

式中，

其中，W_app(τ)是W(τ)的近似函数，由W_i(τ)的和函数组成；m_i和n_i分别是加权函数的系数，对于Schohl卷积近似摩阻，令

则有，

故

对于粘弹性项，有

式中，ε_r(t)为ε_r的时间函数；α为泊松比；e为管道壁厚；J(t*)为蠕变柔量，t*为微分变量；对于广义K-V模型，蠕变柔量函数可以表示为

式中，J_k和τ_k均为粘弹性力学模型参数，可根据蠕变实验数据校核得到；M为校核得到的J_k的个数；对J(t*)求偏导数可得

根据卷积的微分性质，有

令

即

根据卷积的交换性可得

令

则

根据Schohl卷积近似，对于t+Δt时刻，有

故

式中，

CO_A3＝CO_A1·CO_A2 (45)

CO_A4＝CO_A1(1-CO_A2) (46)

所述步骤S4中将源项合并到矩阵方程中，并进行稳定性约束的具体方法为：

采用时间***法将源项引入非线性***求解，求解精度取决于对源项的积分方式；对于二阶精度，采用显式二阶Runge-Kutta离散化方法求解源项s(u)及其空间积分，通过两次更新源项，得到：

式中，

稳定性约束包括对流项的CFL条件以及源项两次更新时的约束；对于CFL条件，有

式中，C_r为库朗数，Δt_max,CFL为CFL条件下的最大时间步长；

利用二阶Runge-Kutta格式引入源项，约束条件为

因此可以得到

因此，计算所允许的最大时间步长需满足

Δt≤min(Δt_max,CFL,Δt_max,s) (54)。

2.根据权利要求1所述的一种基于Schohl卷积近似的管道瞬变流模拟方法，其特征在于，所述步骤S1中管道***控制方程的方程式为：

式中，H为管道某截面处压力；V为管道某截面的平均流速；ε_r为管道迟滞应变；a为有压管道水中波速；g为重力加速度；x为沿管道轴线的距离；t为时间；D为管道内径；ρ为流体密度；f为达西-维斯巴赫系数；ν为水的动力粘度；u为微分变量；W为历史流速变化的加权函数。

3.根据权利要求1所述的一种基于Schohl卷积近似的管道瞬变流模拟方法，其特征在于，所述步骤S2具体为：

其中，

式中，u为求解向量，f(u)为控制体通量矢量，为矢量f中各变量对u中各求解变量的偏导数矩阵，s(u)为源项，分别代表了管道的粘弹性项和摩阻项，若s(u)＝0，即不考虑管道粘弹性和摩阻项，则该***为常系数线性齐次***。

4.根据权利要求1所述的一种基于Schohl卷积近似的管道瞬变流模拟方法，其特征在于，所述步骤S3中计算网格的建立方法为：

在有限体积法求解体系下，将管段分为N个控制体，建立计算网格，每个控制体长Δx，计算时间步长为Δt，对于第i，1≤i≤N，个控制体，其左右边界分别记为i-1/2和i+1/2；为求解该黎曼问题，将矩阵方程沿x方向在第i个控制体两边界界面之间积分得到以下公式：

式中，f_i-1/2表示i-1/2界面的数值通量；f_i+1/2表示i+1/2界面的数值通量；向量u可由控制体i内的平均值表示，故上式可整理为：

式中，上标n表示t时刻，即当前计算时步；上标n+1表示t+Δt时刻，即下一计算时步。

5.根据权利要求4所述的一种基于Schohl卷积近似的管道瞬变流模拟方法，其特征在于，所述步骤S3中网格边界通量的计算方法为：

针对有限体积法Godunov求解格式，控制体边界的数值通量可由界面处的局部黎曼问题得到：对于一个常系数线性齐次双曲***，黎曼问题可描述为下列初值问题：

式中， 和/>是向量/>的向量元素；/> 和/>是向量/>的向量元素；x_i+1/2为第i个控制体右侧边界沿管道轴线的距离；

通过求解特征方程能够得到两个相互独立的特征值：

式中，I为单位矩阵，为要求解的特征值；根据Rankine-Hugoniot条件/> 在t∈[t，t+Δt]的计算时间内，内部边界i+1/2处的通量值/>可表示为：

其中，

对于和/>采用分段多项式重构近似。

6.根据权利要求5所述的一种基于Schohl卷积近似的管道瞬变流模拟方法，其特征在于，所述步骤S3中对于和/>采用二阶Godunov近似，通过引入MUSCL-Hancock格式进行线性插值重构，具体的步骤如下：

A1：数据重构：

通过构造分段线性函数U_i(x)来获取第i个控制体左右边界极值和/>

其中，Ω_i为在第i个计算单元的斜率向量，选择流速限制函数MINMOD斜率限制器：

式中，

A2：时间进化计算：

对于每个控制体[x_i-1/2，x_i+1/2]，左右边界极值按0.5Δt时间进化得到：

式中，分别为时间推进后的左右边界极值，/>和/>分别为左右边界极值处的通量，

A3：黎曼问题的近似求解：

对于边界条件，通过在两端边界各增加两个虚拟控制体-1、0和N+1、N+2，可以得到对边界控制体和内部控制体的统一计算模拟，并达到二阶精度要求；

根据黎曼问题在边界处的解析解：

可推得与正负特征线相关的黎曼不变量方程：

上游水库边界，有

且代入上式得

其虚拟单元满足：

下游边界条件，有

式中，C_V＝(Q₀τ)²/2ΔH₀，Q₀和ΔH₀分别为初始时刻阀门处的流量和水头损失，τ为阀门开度；H_RD为下游水库水头；当C_P-H_RD>0时，S＝+1，当C_P-H_RD<0时，S＝-1；

对于下游阀门+水库边界

根据末端阀门孔口方程，可知

带入阀门正特征线黎曼不变量方程，得

其虚拟单元满足：