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

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

Info

Publication number
CN113656926A
CN113656926A CN202110986943.6A CN202110986943A CN113656926A CN 113656926 A CN113656926 A CN 113656926A CN 202110986943 A CN202110986943 A CN 202110986943A CN 113656926 A CN113656926 A CN 113656926A
Authority
CN
China
Prior art keywords
term
pipeline
schohl
formula
boundary
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202110986943.6A
Other languages
English (en)
Other versions
CN113656926B (zh
Inventor
周领
李赟杰
薛子剑
陆燕清
冯瑞琳
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Hohai University HHU
Original Assignee
Hohai University HHU
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Hohai University HHU filed Critical Hohai University HHU
Priority to CN202110986943.6A priority Critical patent/CN113656926B/zh
Publication of CN113656926A publication Critical patent/CN113656926A/zh
Application granted granted Critical
Publication of CN113656926B publication Critical patent/CN113656926B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/10—Geometric CAD
    • G06F30/18—Network design, e.g. design based on topological or interconnect aspects of utility systems, piping, heating ventilation air conditioning [HVAC] or cabling
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/28—Design optimisation, verification or simulation using fluid dynamics, e.g. using Navier-Stokes equations or computational fluid dynamics [CFD]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/04—Constraint-based CAD
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/10—Numerical modelling
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2113/00—Details relating to the application field
    • G06F2113/08—Fluids
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2113/00—Details relating to the application field
    • G06F2113/14—Pipes
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • General Engineering & Computer Science (AREA)
  • Geometry (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Physics (AREA)
  • Pure & Applied Mathematics (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Computational Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Algebra (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • Operations Research (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • Computing Systems (AREA)
  • Fluid Mechanics (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

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中管道***控制方程的方程式为：
Figure BDA0003231020760000021
Figure BDA0003231020760000022
式中，H为管道某截面处压力；V为管道某截面的平均流速；εr为管道迟滞应变；a为有压管道水中波速；g为重力加速度；x为沿管道轴线的距离；t为时间；D为管道内径；ρ为流体密度；f为达西-维斯巴赫系数；v为水的动力粘度；u为微分变量；W为历史流速变化的加权函数。
进一步地，所述步骤S2具体为：
Figure BDA0003231020760000023
其中，
Figure BDA0003231020760000024
Figure BDA0003231020760000025
式中，u为求解向量，f(u)为控制体通量矢量，
Figure BDA0003231020760000026
为矢量f中各变量对u中各求解变量的偏导数矩阵，s(u)为源项，分别代表了管道的粘弹性项和摩阻项(显式)，若s(u)＝0，即不考虑管道粘弹性和摩阻项，则该***为常系数线性齐次***。
进一步地，所述步骤S3中计算网格的建立方法为：
在有限体积法求解体系下，将管段分为N个控制体，建立计算网格，每个控制体长Δx，计算时间步长为Δt，对于第i(1≤i≤N)个控制体，其左右边界分别记为i-1/2和i+1/2；为求解该黎曼问题，将矩阵方程沿x方向在第i个控制体两边界界面之间积分得到以下公式：
Figure BDA0003231020760000027
式中，fi-1/2表示i-1/2界面的数值通量；fi+1/2表示i+1/2界面的数值通量；向量u可由控制体i内的平均值
Figure BDA0003231020760000028
表示，故上式可整理为：
Figure BDA0003231020760000029
式中，上标n表示t时刻，即当前计算时步；上标n+1表示t+Δt时刻，即下一计算时步。因此，在计算下一时步的向量U时，需涉及当前时步的控制体两边界数值通量和源项。
进一步地，所述步骤S3中网格边界通量的计算方法为：
针对有限体积法Godunov求解格式，控制体边界的数值通量可由界面处的局部黎曼问题得到：对于一个常系数线性齐次双曲***，黎曼问题可描述为下列初值问题：
Figure BDA0003231020760000031
式中，
Figure BDA0003231020760000032
第n个计算时步中向量u在i+1/2界面左侧的平均值，
Figure BDA0003231020760000033
和
Figure BDA0003231020760000034
是向量
Figure BDA0003231020760000035
的向量元素；
Figure BDA0003231020760000036
第n个计算时步中向量u在i+1/2界面右侧的平均值，
Figure BDA0003231020760000037
和
Figure BDA0003231020760000038
是向量
Figure BDA0003231020760000039
的向量元素；xi+1/2为第i个控制体右侧边界沿管道轴线的距离；
通过求解特征方程
Figure BDA00032310207600000310
能够得到两个相互独立的特征值：
Figure BDA00032310207600000311
式中，I为单位矩阵，
Figure BDA00032310207600000312
为要求解的特征值，j＝1～2；根据Rankine-Hugoniot条件
Figure BDA00032310207600000313
Figure BDA00032310207600000314
在t∈[t，t+Δt]的计算时间内，内部边界i+1/2处的通量值
Figure BDA00032310207600000315
可表示为：
Figure BDA00032310207600000316
其中，
Figure BDA00032310207600000317
对于
Figure BDA00032310207600000318
和
Figure BDA00032310207600000319
采用分段多项式重构近似，其精度取决于多项式的形式。
进一步地，所述步骤S3中对于
Figure BDA00032310207600000320
和
Figure BDA00032310207600000321
采用二阶Godunov近似，通过引入MUSCL-Hancock格式进行线性插值重构，具体的步骤如下：
A1：数据重构：
通过构造分段线性函数Ui(x)来获取第i个控制体左右边界极值
Figure BDA00032310207600000322
和
Figure BDA00032310207600000323
Figure BDA00032310207600000324
其中，Ωi为
Figure BDA00032310207600000325
在第i个计算单元的斜率向量，适当选取Ωi的值可以在确保结果没有虚假振荡的同时提高方案的精度。本文选择MINMOD(流速限制函数)斜率限制器：
Figure BDA0003231020760000041
式中，
Figure BDA0003231020760000042
A2：时间进化计算：
对于每个控制体[xi-1/2，xi+1/2]，左右边界极值
Figure BDA0003231020760000043
按0.5Δt时间进化得到：
Figure BDA0003231020760000044
Figure BDA0003231020760000045
式中，
Figure BDA0003231020760000046
分别为时间推进后的左右边界极值，
Figure BDA0003231020760000047
和
Figure BDA0003231020760000048
分别为左右边界极值
Figure BDA0003231020760000049
处的通量，
Figure BDA00032310207600000410
Figure BDA00032310207600000411
A3：黎曼问题的近似求解：
Figure BDA00032310207600000412
对于边界条件，通过在两端边界各增加两个虚拟控制体-1、0和N+1、N+2，可以得到对边界控制体和内部控制体的统一计算模拟，并达到二阶精度要求；
根据黎曼问题在边界处的解析解：
Figure BDA00032310207600000413
可推得与正负特征线相关的黎曼不变量方程：
Figure BDA00032310207600000414
上游水库边界，有
Figure BDA00032310207600000415
且
Figure BDA00032310207600000416
代入上式得
Figure BDA00032310207600000417
其虚拟单元满足：
Figure BDA00032310207600000418
下游边界条件，有
Figure BDA00032310207600000419
式中，CV＝(Q0τ)2/2ΔH0，Q0和ΔH0分别为初始时刻阀门处的流量和水头损失，τ为阀门开度；HRD为下游水库水头；当CP-HRD＞0时，S＝+1，当CP-HRD＜0时，S＝-1；
对于下游阀门+水库边界
Figure BDA0003231020760000051
根据末端阀门孔口方程，可知
Figure BDA0003231020760000052
带入阀门正特征线黎曼不变量方程，得
Figure BDA0003231020760000053
其虚拟单元满足：
Figure BDA0003231020760000054
进一步地，所述步骤S4中对源项中的Zielke非恒定摩阻项和粘弹性项进行Schohl卷积近似的具体方法为：
对于Zielke非恒定摩阻项，在公式(2)中，令
Figure BDA0003231020760000055
式中，
Figure BDA0003231020760000058
其中，Wapp(τ)是W(τ)的近似函数，由Wi(τ)的和函数组成，i＝1～5；mi和ni分别是加权函数的系数，对于Schohl卷积近似摩阻，mi＝(1.051；2.358；9.021；29.47；79.55)；ni＝(26.65；100；669.6；6497；57990)。令
Figure BDA0003231020760000056
则有，
Figure BDA0003231020760000057
故
Figure BDA0003231020760000061
对于粘弹性项，有
Figure BDA0003231020760000062
式中，εr(t)为εr的时间函数；α为泊松比；e为管道壁厚；J(t*)为蠕变柔量，t*为微分变量；对于广义K-V模型，蠕变柔量函数可以表示为
Figure BDA0003231020760000063
式中，Jk和τk均为粘弹性力学模型参数，可根据蠕变实验数据校核得到；M为校核得到的Jk的个数；对J(t*)求偏导数可得
Figure BDA0003231020760000064
根据卷积的微分性质，有
Figure BDA0003231020760000065
令
Figure BDA0003231020760000066
即
Figure BDA0003231020760000067
根据卷积的交换性可得
Figure BDA0003231020760000071
令
Figure BDA0003231020760000072
则
Figure BDA0003231020760000073
根据Schohl卷积近似，对于t+Δt时刻，有
Figure BDA0003231020760000074
故
Figure BDA0003231020760000075
式中，
Figure BDA0003231020760000076
Figure BDA0003231020760000077
COA3＝COA1·COA2 (45)
COA4＝COA1(1-COA2) (46)。
进一步地，所述步骤S4中将源项合并到矩阵方程中，并进行稳定性约束的具体方法为：
在步骤S3中未考虑源项的影响，可采用时间***法将源项引入非线性***求解，求解精度取决于对源项的积分方式；对于二阶精度，采用显式二阶Runge-Kutta离散化方法求解源项s(u)及其空间积分，通过两次更新源项，得到：
Figure BDA0003231020760000078
式中，
Figure BDA0003231020760000081
Figure BDA0003231020760000082
Figure BDA0003231020760000083
稳定性约束主要包括对流项的CFL(Courant-Friedrich-Lewy)条件以及源项两次更新时的约束；对于CFL条件，有
Figure BDA0003231020760000084
式中，Cr为库朗数，Δtmax，CFL为CFL条件下的最大时间步长；
利用二阶Runge-Kutta格式引入源项，约束条件为
Figure BDA0003231020760000085
因此可以得到
Figure BDA0003231020760000086
因此，计算所允许的最大时间步长需满足
Δt≤min(Δtmax，CFL，Δtmax，s) (54)
有益效果：本发明与现有技术相比，具备如下优点：
1、对Zielke原始加权函数模型进行Schohl卷积近似，降低了计算时长。
2、本发明提供的模型提出了新的基于Schohl卷积近似的广义K-V模型求解方法，简化了求解难度，提升了计算精度。
3、使用二阶Godunov格式求解控制方程，降低了复杂管道模拟时的数值耗散程度。
附图说明
图1是本发明的实现流程图；
图2是沿管道轴向的计算网格示意图；
图3是本发明反法模拟结果与实验对比图。
具体实施方式
下面结合附图和具体实施例，进一步阐明本发明，应理解这些实施例仅用于说明本发明而不用于限制本发明的范围，在阅读了本发明之后，本领域技术人员对本发明的各种等价形式的修改均落于本申请所附权利要求所限定的范围。
如图1所示，本发明提供一种基于Schohl卷积近似的管道瞬变流非恒定摩阻和粘弹性模拟方法，包括如下步骤：
S1：构建含Zielke非恒定摩阻项和粘弹性项的管道***控制方程：
Figure BDA0003231020760000095
Figure BDA0003231020760000096
式中，H为管道某截面处压力；V为管道某截面的平均流速；εr为管道迟滞应变；a为有压管道水中波速；g为重力加速度；x为沿管道轴线的距离；t为时间；D为管道内径；ρ为流体密度；f为达西-维斯巴赫系数；ν为水的动力粘度；u为微分变量；W为历史流速变化的加权函数。
S2：将控制方程以黎曼问题的求解格式表示为矩阵形式，将非恒定摩阻项和粘弹性项放入源项：
Figure BDA0003231020760000091
其中，
Figure BDA0003231020760000092
Figure BDA0003231020760000093
式中，u为求解向量，f(u)为控制体通量矢量，
Figure BDA0003231020760000094
为矢量f中各变量对u中各求解变量的偏导数矩阵，s(u)为源项，分别代表了管道的粘弹性项和摩阻项(显式)，若s(u)＝0，即不考虑管道粘弹性和摩阻项，则该***为常系数线性齐次***。
S3：一、建立有限体积法求解体系下的计算网格：
如图2所示，在有限体积法求解体系下，将管段分为N个控制体，建立计算网格，每个控制体长Δx，计算时间步长为Δt，对于第i(1≤i≤N)个控制体，其左右边界分别记为i-1/2和i+1/2；为求解该黎曼问题，将矩阵方程沿x方向在第i个控制体两边界界面之间积分得到以下公式：
Figure BDA0003231020760000101
式中，fi-1/2表示i-1/2界面的数值通量；fi+l/2表示i+1/2界面的数值通量；向量u可由控制体i内的平均值
Figure BDA0003231020760000102
表示，故上式可整理为：
Figure BDA0003231020760000103
式中，上标n表示t时刻，即当前计算时步；上标n+1表示t+Δt时刻，即下一计算时步。因此，在计算下一时步的向量U时，需涉及当前时步的控制体两边界数值通量和源项。
二、通过二阶Godunov格式计算网格边界通量：
针对有限体积法Godunov求解格式，控制体边界的数值通量可由界面处的局部黎曼问题得到：对于一个常系数线性齐次双曲***，黎曼问题可描述为下列初值问题：
Figure BDA0003231020760000104
式中，
Figure BDA0003231020760000105
第n个计算时步中向量u在i+1/2界面左侧的平均值，
Figure BDA0003231020760000106
和
Figure BDA0003231020760000107
是向量
Figure BDA0003231020760000108
的向量元素；
Figure BDA0003231020760000109
第n个计算时步中向量u在i+1/2界面右侧的平均值，
Figure BDA00032310207600001010
和
Figure BDA00032310207600001011
是向量
Figure BDA00032310207600001012
的向量元素；xi+1/2为第i个控制体右侧边界沿管道轴线的距离；
通过求解特征方程
Figure BDA00032310207600001013
能够得到两个相互独立的特征值：
Figure BDA00032310207600001014
式中，I为单位矩阵，
Figure BDA00032310207600001015
为要求解的特征值，j＝1～2；根据Rankine-Hugoniot条件
Figure BDA00032310207600001016
Figure BDA00032310207600001017
在t∈[t，t+Δt]的计算时间内，内部边界i+1/2处的通量值
Figure BDA00032310207600001018
可表示为：
Figure BDA0003231020760000111
其中，
Figure BDA0003231020760000112
对于
Figure BDA0003231020760000113
和
Figure BDA0003231020760000114
采用分段多项式重构近似，其精度取决于多项式的形式。
对于
Figure BDA0003231020760000115
和
Figure BDA0003231020760000116
采用二阶Godunov近似，通过引入MUSCL-Hancock格式进行线性插值重构，具体的步骤如下：
A1：数据重构：
通过构造分段线性函数ui(x)来获取第i个控制体左右边界极值
Figure BDA0003231020760000117
和
Figure BDA0003231020760000118
Figure BDA0003231020760000119
其中，Ωi为
Figure BDA00032310207600001110
在第i个计算单元的斜率向量，适当选取Ωi的值可以在确保结果没有虚假振荡的同时提高方案的精度。本文选择MINMOD(流速限制函数)斜率限制器：
Figure BDA00032310207600001111
式中，
Figure BDA00032310207600001112
A2：时间进化计算：
对于每个控制体[xi-1/2，xi+1/2]，左右边界极值
Figure BDA00032310207600001113
按0.5Δt时间进化得到：
Figure BDA00032310207600001114
Figure BDA00032310207600001115
式中，
Figure BDA00032310207600001116
分别为时间推进后的左右边界极值，
Figure BDA00032310207600001117
和
Figure BDA00032310207600001118
分别为左右边界极值
Figure BDA00032310207600001119
处的通量，
Figure BDA00032310207600001120
Figure BDA00032310207600001121
A3：黎曼问题的近似求解：
Figure BDA00032310207600001122
对于边界条件，通过在两端边界各增加两个虚拟控制体-1、0和N+1、N+2，可以得到对边界控制体和内部控制体的统一计算模拟，并达到二阶精度要求；
根据黎曼问题在边界处的解析解：
Figure BDA0003231020760000121
可推得与正负特征线相关的黎曼不变量方程：
Figure BDA0003231020760000122
上游水库边界，有
Figure BDA0003231020760000123
且
Figure BDA0003231020760000124
代入上式得
Figure BDA0003231020760000125
其虚拟单元满足：
Figure BDA0003231020760000126
下游边界条件，有
Figure BDA0003231020760000127
式中，CV＝(Q0τ)2/2ΔH0，Q0和ΔH0分别为初始时刻阀门处的流量和水头损失，τ为阀门开度；HRD为下游水库水头；当CP-HRD>0时，S＝+1，当CP-HRD<0时，S＝-1；
对于下游阀门+水库边界
Figure BDA0003231020760000128
根据末端阀门孔口方程，可知
Figure BDA0003231020760000129
带入阀门正特征线黎曼不变量方程，得
Figure BDA00032310207600001210
其虚拟单元满足：
Figure BDA00032310207600001211
S4：一、对源项中的Zielke非恒定摩阻项和粘弹性项进行Schohl卷积近似：
对于Zielke非恒定摩阻项，在公式(2)中，令
Figure BDA00032310207600001212
式中，
Figure BDA0003231020760000137
其中，Wapp(τ)是W(τ)的近似函数，由Wi(τ)的和函数组成，i＝1～5；mi和ni分别是加权函数的系数，对于Schohl卷积近似摩阻，mi＝(1.051；2.358；9.021；29.47；79.55)；ni＝(26.65；100；669.6；6497；57990)。令
Figure BDA0003231020760000131
则有，
Figure BDA0003231020760000132
故
Figure BDA0003231020760000133
对于粘弹性项，有
Figure BDA0003231020760000134
式中，εr(t)为εr的时间函数；α为泊松比；e为管道壁厚；J(t*)为蠕变柔量，t*为微分变量；对于广义K-V模型，蠕变柔量函数可以表示为
Figure BDA0003231020760000135
式中，Jk和τk均为粘弹性力学模型参数，可根据蠕变实验数据校核得到；M为校核得到的Jk的个数；对J(t*)求偏导数可得
Figure BDA0003231020760000136
根据卷积的微分性质，有
Figure BDA0003231020760000141
令
Figure BDA0003231020760000142
即
Figure BDA0003231020760000143
根据卷积的交换性可得
Figure BDA0003231020760000144
令
Figure BDA0003231020760000145
则
Figure BDA0003231020760000146
根据Schohl卷积近似，对于t+Δt时刻，有
Figure BDA0003231020760000147
故
Figure BDA0003231020760000148
式中，
Figure BDA0003231020760000151
Figure BDA0003231020760000152
COA3＝COA1·COA2 (45)
COA4＝COA1(1-COA2) (46)
二、将源项合并到矩阵方程中，并进行稳定性约束：
在步骤S3中未考虑源项的影响，可采用时间***法将源项引入非线性***求解，求解精度取决于对源项的积分方式；对于二阶精度，采用显式二阶Runge-Kutta离散化方法求解源项s(u)及其空间积分，通过两次更新源项，得到：
Figure BDA0003231020760000153
式中，
Figure BDA0003231020760000154
Figure BDA0003231020760000155
Figure BDA0003231020760000156
稳定性约束主要包括对流项的CFL(Courant-Friedrich-Lewy)条件以及源项两次更新时的约束；对于CFL条件，有
Figure BDA0003231020760000157
式中，Cr为库朗数，Δtmax，CFL为CFL条件下的最大时间步长；
利用二阶Runge-Kutta格式引入源项，约束条件为
Figure BDA0003231020760000158
因此可以得到
Figure BDA0003231020760000161
因此，计算所允许的最大时间步长需满足
Δt≤min(Δtmax，CFL，Δtmax，s) (54)
S5：结合实际工程给出计算初始条件，利用Fortran编程得出计算结果：
带入实际工况中的初始条件，包括上游水库或水箱水头Hr，管道直径D，根据管道长度所划分的控制体的个数N，以及每段控制体的长度Δx，管道波速a，雷诺数Re等，可求得管道中各个控制体的测压管水头H和流速V，从而实现了管道水力瞬变的模拟。
基于上述方案，本实施例中为了验证本发明方法的效果，选取Covas于2003年设计搭建的粘弹性管道水力瞬变实验装置***用于验证本发明方法的有效性，实验和模拟对比结果如图3所示，结果表明本发明所提新的方法能够很好的模拟管道水力瞬变的压力峰值、衰减和相位。

Claims (8)

1.一种基于Schohl卷积近似的管道瞬变流模拟方法，其特征在于，包括如下步骤：
S1：构建含Zielke非恒定摩阻项和粘弹性项的管道***控制方程；
S2：将控制方程以黎曼问题的求解格式表示为矩阵形式，将非恒定摩阻项和粘弹性项放入源项；
S3：建立有限体积法求解体系下的计算网格，通过二阶Godunov格式计算网格边界通量；
S4：对源项中的Zielke非恒定摩阻项和粘弹性项进行Schohl卷积近似，将源项合并到矩阵方程中，并进行稳定性约束；
S5：结合实际工程给出计算初始条件，得出计算结果。
2.根据权利要求1所述的一种基于Schohl卷积近似的管道瞬变流模拟方法，其特征在于，所述步骤S1中管道***控制方程的方程式为：
Figure FDA0003231020750000011
Figure FDA0003231020750000012
式中，H为管道某截面处压力；V为管道某截面的平均流速；εr为管道迟滞应变；a为有压管道水中波速；g为重力加速度；x为沿管道轴线的距离；t为时间；D为管道内径；ρ为流体密度；f为达西-维斯巴赫系数；ν为水的动力粘度；u为微分变量；W为历史流速变化的加权函数。
3.根据权利要求1所述的一种基于Schohl卷积近似的管道瞬变流模拟方法，其特征在于，所述步骤S2具体为：
Figure FDA0003231020750000013
其中，
Figure FDA0003231020750000014
Figure FDA0003231020750000015
式中，u为求解向量，f(u)为控制体通量矢量，
Figure FDA0003231020750000016
为矢量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个控制体两边界界面之间积分得到以下公式：
Figure FDA0003231020750000021
式中，fi-1/2表示i-1/2界面的数值通量；fi+1/2表示i+1/2界面的数值通量；向量u可由控制体i内的平均值
Figure FDA0003231020750000022
表示，故上式可整理为：
Figure FDA0003231020750000023
式中，上标n表示t时刻，即当前计算时步；上标n+1表示t+Δt时刻，即下一计算时步。
5.根据权利要求4所述的一种基于Schohl卷积近似的管道瞬变流模拟方法，其特征在于，所述步骤S3中网格边界通量的计算方法为：
针对有限体积法Godunov求解格式，控制体边界的数值通量可由界面处的局部黎曼问题得到：对于一个常系数线性齐次双曲***，黎曼问题可描述为下列初值问题：
Figure FDA0003231020750000024
式中，
Figure FDA0003231020750000025
Figure FDA0003231020750000026
和
Figure FDA0003231020750000027
是向量
Figure FDA0003231020750000028
的向量元素；
Figure FDA0003231020750000029
Figure FDA00032310207500000210
和
Figure FDA00032310207500000211
是向量
Figure FDA00032310207500000212
的向量元素；xi+1/2为第i个控制体右侧边界沿管道轴线的距离；
通过求解特征方程
Figure FDA00032310207500000213
能够得到两个相互独立的特征值：
Figure FDA00032310207500000214
式中，I为单位矩阵，
Figure FDA00032310207500000215
为要求解的特征值；根据Rankine-Hugoniot条件
Figure FDA00032310207500000216
Figure FDA00032310207500000217
在t∈[t，t+Δt]的计算时间内，内部边界i+1/2处的通量值
Figure FDA00032310207500000218
可表示为：
Figure FDA00032310207500000219
其中，
Figure FDA0003231020750000031
对于
Figure FDA0003231020750000032
和
Figure FDA0003231020750000033
采用分段多项式重构近似。
6.根据权利要求5所述的一种基于Schohl卷积近似的管道瞬变流模拟方法，其特征在于，所述步骤S3中对于
Figure FDA0003231020750000034
和
Figure FDA0003231020750000035
采用二阶Godunov近似，通过引入MUSCL-Hancock格式进行线性插值重构，具体的步骤如下：
A1：数据重构：
通过构造分段线性函数Ui(x)来获取第i个控制体左右边界极值
Figure FDA0003231020750000036
和
Figure FDA0003231020750000037
Figure FDA0003231020750000038
其中，Ωi为
Figure FDA0003231020750000039
在第i个计算单元的斜率向量，选择MINMOD(流速限制函数)斜率限制器：
Figure FDA00032310207500000310
式中，
Figure FDA00032310207500000311
A2：时间进化计算：
对于每个控制体[xi-1/2，xi+1/2]，左右边界极值
Figure FDA00032310207500000312
按0.5Δt时间进化得到：
Figure FDA00032310207500000313
Figure FDA00032310207500000314
式中，
Figure FDA00032310207500000315
分别为时间推进后的左右边界极值，
Figure FDA00032310207500000316
和
Figure FDA00032310207500000317
分别为左右边界极值
Figure FDA00032310207500000318
处的通量，
Figure FDA00032310207500000319
Figure FDA00032310207500000320
A3：黎曼问题的近似求解：
Figure FDA00032310207500000321
对于边界条件，通过在两端边界各增加两个虚拟控制体-1、0和N+1、N+2，可以得到对边界控制体和内部控制体的统一计算模拟，并达到二阶精度要求；
根据黎曼问题在边界处的解析解：
Figure FDA0003231020750000041
可推得与正负特征线相关的黎曼不变量方程：
C±:
Figure FDA0003231020750000042
上游水库边界，有
Figure FDA0003231020750000043
且
Figure FDA0003231020750000044
代入上式得
Figure FDA0003231020750000045
其虚拟单元满足：
Figure FDA0003231020750000046
下游边界条件，有
Figure FDA0003231020750000047
式中，CV＝(Q0τ)2/2ΔH0，Q0和ΔH0分别为初始时刻阀门处的流量和水头损失，τ为阀门开度；HRD为下游水库水头；当CP-HRD>0时，S＝+1，当CP-HRD<0时，S＝-1；
对于下游阀门+水库边界
Figure FDA0003231020750000048
根据末端阀门孔口方程，可知
Figure FDA0003231020750000049
带入阀门正特征线黎曼不变量方程，得
Figure FDA00032310207500000410
其虚拟单元满足：
Figure FDA00032310207500000411
7.根据权利要求2所述的一种基于Schohl卷积近似的管道瞬变流模拟方法，其特征在于，所述步骤S4中对源项中的Zielke非恒定摩阻项和粘弹性项进行Schohl卷积近似的具体方法为：
对于Zielke非恒定摩阻项，在公式(2)中，令
Figure FDA00032310207500000412
式中，
Figure FDA0003231020750000051
其中，Wapp(τ)是W(τ)的近似函数，由Wi(τ)的和函数组成；mi和ni分别是加权函数的系数，对于Schohl卷积近似摩阻，令
Figure FDA0003231020750000052
则有，
Figure FDA0003231020750000053
故
Figure FDA0003231020750000054
对于粘弹性项，有
Figure FDA0003231020750000055
式中，εr(t)为εr的时间函数；α为泊松比；e为管道壁厚；J(t*)为蠕变柔量，t*为微分变量；对于广义K-V模型，蠕变柔量函数可以表示为
Figure FDA0003231020750000056
式中，Jk和τk均为粘弹性力学模型参数，可根据蠕变实验数据校核得到；M为校核得到的Jk的个数；对J(t*)求偏导数可得
Figure FDA0003231020750000057
根据卷积的微分性质，有
Figure FDA0003231020750000061
令
Figure FDA0003231020750000062
即
Figure FDA0003231020750000063
根据卷积的交换性可得
Figure FDA0003231020750000064
令
Figure FDA0003231020750000065
则
Figure FDA0003231020750000066
根据Schohl卷积近似，对于t+Δt时刻，有
Figure FDA0003231020750000067
故
Figure FDA0003231020750000068
式中，
Figure FDA0003231020750000069
Figure FDA0003231020750000071
COA3＝COA1·COA2 (45)
COA4＝COA1(1-COA2) (46)。
8.根据权利要求7所述的一种基于Schohl卷积近似的管道瞬变流模拟方法，其特征在于，所述步骤S4中将源项合并到矩阵方程中，并进行稳定性约束的具体方法为：
采用时间***法将源项引入非线性***求解，求解精度取决于对源项的积分方式；对于二阶精度，采用显式二阶Runge-Kutta离散化方法求解源项s(u)及其空间积分，通过两次更新源项，得到：
Figure FDA0003231020750000072
式中，
Figure FDA0003231020750000073
Figure FDA0003231020750000074
Figure FDA0003231020750000075
稳定性约束包括对流项的CFL(Courant-Friedrich-Lewy)条件以及源项两次更新时的约束；对于CFL条件，有
Figure FDA0003231020750000076
式中，Cr为库朗数，Δtmax,CFL为CFL条件下的最大时间步长；
利用二阶Runge-Kutta格式引入源项，约束条件为
Figure FDA0003231020750000077
因此可以得到
Figure FDA0003231020750000081
因此，计算所允许的最大时间步长需满足
Δt≤min(Δtmax,CFL,Δtmax,s) (54)。
CN202110986943.6A 2021-08-26 2021-08-26 基于Schohl卷积近似的管道瞬变流模拟方法 Active CN113656926B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110986943.6A CN113656926B (zh) 2021-08-26 2021-08-26 基于Schohl卷积近似的管道瞬变流模拟方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110986943.6A CN113656926B (zh) 2021-08-26 2021-08-26 基于Schohl卷积近似的管道瞬变流模拟方法

Publications (2)

Publication Number Publication Date
CN113656926A true CN113656926A (zh) 2021-11-16
CN113656926B CN113656926B (zh) 2024-03-26

Family

ID=78492897

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110986943.6A Active CN113656926B (zh) 2021-08-26 2021-08-26 基于Schohl卷积近似的管道瞬变流模拟方法

Country Status (1)

Country Link
CN (1) CN113656926B (zh)

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111767684A (zh) * 2020-06-30 2020-10-13 西安理工大学 一种优化的摩阻源项隐式格式二维浅水方程建模方法
CN114896908A (zh) * 2022-05-19 2022-08-12 四川大学 一种基于通量限制器的水滴流场及水滴收集率计算方法
CN115391069A (zh) * 2022-10-27 2022-11-25 山东省计算中心（国家超级计算济南中心） 基于海洋模式roms的并行通讯方法及***
CN116663146B (zh) * 2023-05-30 2023-11-17 西安理工大学 非圆形管道沿程阻力的计算方法

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109918787A (zh) * 2019-03-08 2019-06-21 河海大学 基于有限体积法的输水管道内水气两相均质流的模拟方法
CN111339701A (zh) * 2020-02-26 2020-06-26 河海大学 基于Brunone动态摩阻的管道泄漏特性Godunov模拟方法
CN111695307A (zh) * 2020-05-21 2020-09-22 河海大学 显式考虑动态摩阻的水锤有限体积模拟方法
CN111985166A (zh) * 2020-07-31 2020-11-24 河海大学 隐式考虑动态摩阻的管道水力瞬变模拟方法和存储介质

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109918787A (zh) * 2019-03-08 2019-06-21 河海大学 基于有限体积法的输水管道内水气两相均质流的模拟方法
CN111339701A (zh) * 2020-02-26 2020-06-26 河海大学 基于Brunone动态摩阻的管道泄漏特性Godunov模拟方法
CN111695307A (zh) * 2020-05-21 2020-09-22 河海大学 显式考虑动态摩阻的水锤有限体积模拟方法
CN111985166A (zh) * 2020-07-31 2020-11-24 河海大学 隐式考虑动态摩阻的管道水力瞬变模拟方法和存储介质

Cited By (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111767684A (zh) * 2020-06-30 2020-10-13 西安理工大学 一种优化的摩阻源项隐式格式二维浅水方程建模方法
CN111767684B (zh) * 2020-06-30 2024-04-19 西安理工大学 一种优化的摩阻源项隐式格式二维浅水方程建模方法
CN114896908A (zh) * 2022-05-19 2022-08-12 四川大学 一种基于通量限制器的水滴流场及水滴收集率计算方法
CN114896908B (zh) * 2022-05-19 2023-03-28 四川大学 一种基于通量限制器的水滴流场及水滴收集率计算方法
CN115391069A (zh) * 2022-10-27 2022-11-25 山东省计算中心（国家超级计算济南中心） 基于海洋模式roms的并行通讯方法及***
CN116663146B (zh) * 2023-05-30 2023-11-17 西安理工大学 非圆形管道沿程阻力的计算方法

Also Published As

Publication number Publication date
CN113656926B (zh) 2024-03-26

Similar Documents

Publication Publication Date Title
CN113656926A (zh) 基于Schohl卷积近似的管道瞬变流模拟方法
Coulaud et al. Numerical modelling of nonlinear effects in laminar flow through a porous medium
Alamian et al. A state space model for transient flow simulation in natural gas pipelines
Dupuis et al. Coupling lattice Boltzmann and molecular dynamics models for dense fluids
Massari et al. Is the leak head–discharge relationship in polyethylene pipes a bijective function?
CN111339701B (zh) 基于Brunone动态摩阻的管道泄漏特性Godunov模拟方法
CN103729505B (zh) 一种基于cfd的阀门当量长度计算方法
LaViolette On the history, science, and technology included in the Moody diagram
CN111985166A (zh) 隐式考虑动态摩阻的管道水力瞬变模拟方法和存储介质
CN109855691B (zh) 一种差分式层流流量测量方法及装置
CN109117579B (zh) 一种多孔孔板流量计的设计计算方法
Zhou et al. An accurate and efficient scheme involving unsteady friction for transient pipe flow
CN113553737B (zh) 一种基于阀门压差的阀门流量预测方法
CN107014451A (zh) 基于广义回归神经网络推测超声波流量传感器系数的方法
CN110426085A (zh) 一种节流式流量测量装置流量算法
Hullender Alternative approach for modeling transients in smooth pipe with low turbulent flow
CN111695307B (zh) 显式考虑动态摩阻的水锤有限体积模拟方法
Kane Sphere drag data at supersonic speeds and low Reynolds numbers
Pal et al. Efficient approach toward the application of the Godunov method to hydraulic transients
CN116734968A (zh) 一种流量系数的校验方法、装置、设备及存储介质
Li et al. 1D-3D coupling investigation of hydraulic transient for power-supply failure of centrifugal pump-pipe system
Al-Khomairi USE OF THE STEADY-STATE ORIFICE EQUATION IN THE COMPUTATION OF TRANSIENT FLOW THROUGH PIPE LEAKS.
Shermukhamedov et al. Mathematical model of dynamic characteristics of a hydraulic drive with distributed parameters taking into account the influence of the temperature factor
CN110608891B (zh) 一种液体火箭发动机冷流试验***及并联贮箱推进剂输送均衡性能试验方法
CN210400858U (zh) 一种液体火箭发动机冷流试验***

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant