CN109918787B - 基于有限体积法的输水管道内气液两相均质流的模拟方法 - Google Patents

基于有限体积法的输水管道内气液两相均质流的模拟方法 Download PDF

Info

Publication number
CN109918787B
CN109918787B CN201910173947.5A CN201910173947A CN109918787B CN 109918787 B CN109918787 B CN 109918787B CN 201910173947 A CN201910173947 A CN 201910173947A CN 109918787 B CN109918787 B CN 109918787B
Authority
CN
China
Prior art keywords
equation
flow
pipeline
gas
model
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.)
Active
Application number
CN201910173947.5A
Other languages
English (en)
Other versions
CN109918787A (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 CN201910173947.5A priority Critical patent/CN109918787B/zh
Publication of CN109918787A publication Critical patent/CN109918787A/zh
Application granted granted Critical
Publication of CN109918787B publication Critical patent/CN109918787B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Landscapes

  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明公开了一种基于有限体积法的输水管道内水气两相均质流的模拟方法，这种方法基于非守恒的形式；通过有限体积法对模型管道在时间和空间上进行离散；接着，通过含有TVD的MUSCL‑Hancoke格式，得到重建的流动变量，使其具有二阶精度；接着，求解模型内部界面处黎曼解进而得到计算界面通量值；接着，针对模型的边界条件，使用牛顿迭代法利用黎曼不变量等到二阶的边界值；再接着，在模型中进入了源项，通过数值方法进行离散，使其具有二阶精度；最后，基于初始给定的参数(自由气体含量、常压下纯水水锤波速等)，模拟出任意时刻管道内的压力波动。结果表明，相较于传统的守恒形式均质流模型，该模型的计算精度和计算效率更高。

Description

基于有限体积法的输水管道内气液两相均质流的模拟方法
技术领域
本发明涉及一种基于有限体积法的输水管道内气-液两相均质流的模拟方法，属于水电站(泵站)水力学数值模拟计算技术领域。
背景技术
在输水管道系统中，阀门或机组的突然启闭可能导致管道内的压力突变，当管道内的压力变化超过管壁材料的承受能力时，可能会造成严重的管道系统破坏，甚至会威胁人身安全。近些年，管道内的两相流问题一直是国内外研究的热门课题，不同于一般的单相流问题，两相流需要考虑不同相间的耦合作用。在两相流问题中，均质流模型被认为是应用最广的模型。在输水管道系统中，若有少量的自由气体均匀的分布在水体中，且假设气-液两相不会随着流动发生相互移动，这样的模型被成为两相均质流模型。
目前，针对输水管道内的气-液两相均质流问题，其模拟的方法主要为特征线法(MOC,Method ofCharacteristics)以及有限差分法(FDM,Finite Difference Method)。MOC由于其计算简单，可以较好的模拟管道内压力波动而被广泛应用。然而，由于均质流模型变波速、变密度的问题，MOC需要做插值运算而造成模型预测与实际相差较大。FDM则避免了线性插值问题，在计算精度上相较于MOC提升了不少。随着均质流模型的愈发成熟，寻求更高精度的模型成为了学术界和工业界的目标。最近，León等人提出了基于有限体积法戈杜诺夫格式的二阶高分辨率均质流模型，该模型采用了TVD(Total VariationDiminishing总变差减小)格式，同时在对流部分采用CFL(Courant-Friedrichs-Lewy)约束，在他们的模型中，黎曼问题的求解只需简单的代数过程不需要通过迭代，这大大的提高了模型的计算效率。尽管如此，León等人提出的模型相较于MOC等模型，仍然存在计算效率低下的问题，因此，高效率高精度模型的研究依然是国内外一个重要的课题。
发明内容
发明目的：为弥补现有技术在模拟输水管道中气-液两相均质流时存在计算精度低等不足，本发明基于有限体积法戈杜诺夫格式，通过采用非守恒形式的控制方程，提供了一种算法简单，易于实现的模拟方法，以在高计算效率的前提下得到高精度解。
技术方案：为解决上述问题，本发明采用以下技术手段：
一种基于有限体积法的输水管道内气-液两相均质流的模拟方法，其特征在于，采用了二阶非守恒变量的有限体积法戈杜诺夫格式来模拟管道系统内均质两相流瞬变过程，具体步骤如下：
步骤1：在非守恒形式的欧拉体系下，构建含自由气体的两相均质流基本控制方程，根据模拟工况确定计算域、初始条件以及边界条件；
步骤2：通过有限体积法戈杜诺夫格式划分计算网格，并建立离散方程；
步骤3：通过含有TVD形式的MUSCL-Hancoke格式，得到重建的流动变量，使其具有二阶精度；
步骤4：求解模型内部界面黎曼解，得到控制单元界面通量；
步骤5：根据已有的边界条件，通过牛顿迭代法利用黎曼不变量得到二阶的边界值；
步骤6：引入源项，通过基于二阶龙格库塔法离散格式的时间分裂法，求得时间梯度上二阶FVM-Godunov离散方程；
步骤7：建立稳定性约束条件，更新初始变量进行下一时步计算。
进一步地，步骤1中，在非守恒形式的欧拉体系下，构建含自由气体的两相均质流基本控制方程，需在水锤问题的基础上假定：(a)管道内自由气体体积含量很低(<≈1％)且气体在水体中均匀分散，气-水两相流以等效单相流体处理； (b)瞬变过程时间尺度很小，忽略整个过程中气体的吸收与释放；(c)管道内流体为无粘性流动且整个过程发生在等温条件下不考虑热传递；(d)管道和水体均为刚性，管道截面积在整个瞬变流过程中不发生变化。
进一步地，步骤1中，构建的含自由气体均质两相流基本控制方程包含：
(1)水锤的基本方程：
Figure GDA0002982439880000021
Figure GDA0002982439880000022
其中，沿管线距离x与时间t是自变量；ρ(x,t)是平均截面流体密度；v(x,t)是平均截面速率；g是重力加速度；J为单位管长的摩擦力；ρl为水体的密度；
(2)气-液两相流水锤波速与压力的关系：
Figure GDA0002982439880000031
其中，cm为混合流体密度；C为纯液体下水锤波速；Pref为参考压力；ρfref为参考压力下两相流的密度；ψref为参考压力下混合流体内气体的初始体积分数；
(3)气-液两相流压力与密度的关系：
Figure GDA0002982439880000032
其中，ρ为混合流体的密度，P为流体的压力。
进一步地，步骤2中，在Godunov格式下建立均质两相流的计算网格的方法为：
(a)建立初始网格：
将空间域X离散为N个长度为Δx、时间域T上离散为间隔为Δt的控制单元，在空间处，第i个网格以i-1/2、i+1/2为边界，编号为i，用其代表该单元处流体参数的平均值；
(b)建立虚拟网格：
为了方便模型计算，在模型两边边界处各建立两个虚拟单元，分别编号为-1、 0、N+1、N+2并假设虚拟单元处的流动参数与边界处的参数相同。
进一步地，步骤2中，在管道内两相流离散模型系统基础上，建立离散方程的步骤为：
(a)将微分方程(1)(2)转化为拟线性的非守恒向量格式：
Figure GDA0002982439880000033
其中
Figure GDA0002982439880000034
D为管道内径，f为达西韦斯巴赫系数；
Figure GDA0002982439880000035
Figure GDA0002982439880000036
P为管道内平均截面的压力；F(W)为非守恒系统向量形式的通量；
(b)忽略初始常量，对方程(5)中的通量进行积分：
Figure GDA0002982439880000041
(c)在从界面i-1/2到i+1/2的单元i及从t到t+Δt的时间段Δt内对方程(5) 积分，得到流动变量W的离散方程：
Figure GDA0002982439880000042
其中，上标n，n+1/2和n+1分别代表t、t+1/2Δt和t+Δt时步；
Figure GDA0002982439880000043
为W在整个单元处流体参数的平均值；F为界面处的通量。
进一步地，步骤3中，通过使用TVD版本的MUSCL-Hancock格式，利用斜率限制器得到二阶精度的格式的方法为：
(a)引入斜率限制器
Figure GDA0002982439880000044
其中，
Figure GDA0002982439880000045
为斜率限制器参数；
(b)非守恒变量的重建
Figure GDA0002982439880000046
(c)进一步重建
Figure GDA0002982439880000047
(d)黎曼问题
Figure GDA0002982439880000048
进一步地，步骤4中，求解黎曼问题，并通过黎曼解得到模型内部界面通量的方法为：
(a)一般的双曲系统的黎曼问题也就是初始值问题：
Wt+Fx＝0 (12)
Figure GDA0002982439880000051
这里
Figure GDA0002982439880000052
和
Figure GDA0002982439880000053
均为常数值，F为界面处通量；对方程(6)在单元处进行积分，得到：
Figure GDA0002982439880000054
(b)在有限体积法中，穿过断点处，利用Rankine-Hugoniot条件ΔF＝λΔW，得到：
Figure GDA0002982439880000055
(c)利用压力与波速的一般方程
Figure GDA0002982439880000056
推导出：
Figure GDA0002982439880000057
(d)通过一阶的线性近似：
Figure GDA0002982439880000058
其中，
Figure GDA0002982439880000059
并且当网格足够精细时，线性近似足够保证精度；
(e)对微分方程(17)在单元处进行积分：
Figure GDA00029824398800000510
(f)求解黎曼问题中中间处参数值：
Figure GDA00029824398800000511
(g)求解单元界面通量值：
Figure GDA00029824398800000512
其中，
Figure GDA00029824398800000513
根据求出的变量ρ*的值，通过方程(4)，利用牛顿迭代法可求解得到方程(20)中的压力P。
进一步地，步骤5中，模型中二阶边界的求解方法为：
(a)以左侧为例，当已知边界处的压力值
Figure GDA0002982439880000061
时：
在边界处，沿着特征线dx/dt＝v-cm，满足
Figure GDA0002982439880000062
为了得到二阶精度的边界解，边界处引入一个虚拟网格，在左侧单元内从(1/2+,tn)到(0, tn+1/2)对微分关系式
Figure GDA0002982439880000063
进行积分，得到：
Figure GDA0002982439880000064
其中，
Figure GDA0002982439880000065
和
Figure GDA0002982439880000066
为已知量，
Figure GDA0002982439880000067
为已知的边界条件，未知量
Figure GDA0002982439880000068
可以通过以下方程求解得出：
Figure GDA0002982439880000069
Figure GDA00029824398800000610
结合方程(21)，待求值
Figure GDA00029824398800000611
可以被推导为：
Figure GDA00029824398800000612
(b)以模型左侧为例，当已知边界处的流量值
Figure GDA00029824398800000613
时：
模型左侧边界处的流速可由
Figure GDA00029824398800000614
得到，利用方程(3)，(4)通过牛顿迭代法进一步求得未知量
Figure GDA00029824398800000615
首先假设变量
Figure GDA00029824398800000616
Figure GDA00029824398800000617
然后将该值带入方程(4)计算得到
Figure GDA00029824398800000618
将计算得到的
Figure GDA00029824398800000619
带入方程(3)计算得到
Figure GDA00029824398800000620
将这些参数带入方程(25)进一步得到
Figure GDA00029824398800000621
利用新的到的参数得到
Figure GDA00029824398800000622
重复以上步骤直到计算收敛；一旦得到
Figure GDA00029824398800000623
和
Figure GDA00029824398800000624
的值,
Figure GDA00029824398800000625
就可以利用公式计算得出：
Figure GDA00029824398800000626
进一步地，步骤6中，引入源项通过二阶龙格库塔法求解离散方程：
Figure GDA0002982439880000071
其中，
Figure GDA0002982439880000072
为n+1时步，控制单元i在纯对流时，流动变量W的通量；
Figure GDA0002982439880000073
为采用时间分裂法第一次更新后的通量。
进一步地，步骤7中，建立有效的稳定性约束条件，更新初始变量并开展下一时步计算的方法为：
(a)由于采用了显式二阶龙格-库塔离散法将S引入求解，稳定性约束不仅要包括对流部分的CFL准则，而且还要包括源项的约束；由CFL得到:
Figure GDA0002982439880000074
(b)显式二阶龙格-库塔离散化约束：
Figure GDA0002982439880000075
(c)由于对流项和源项采用相同的时间步长Δt,因而采用
Figure GDA0002982439880000076
而不是
Figure GDA0002982439880000077
最后给出包括对流部分和源项的最大允许时间步长:
Figure GDA0002982439880000078
有益效果：与现有技术相比，本发明具有如下优点：
(1)本发明提供的基于非守恒形式有限体积法气-液两相均质流模型收敛到了正确结果并取得二阶精度，同时如MOC类方法一样简单且易于实现；(2)相较于León等人提出的二阶精度模型，在达到同一准确度的条件下，本发明的计算时长更短，计算效率更高；(3)相较于León等人提出的二阶精度模型，在同一计算时长的条件下，本发明的计算精确度更高。
附图说明
图1为本发明的基本流程图；
图2为实施例下模型的初始网格左侧边界离散系统图；
图3为实施例对比实验装置示意图；
图4为实施例下，压力修正系数C-ap＝1.0时，阀门末端的压力曲线图；
图5为实施例下，压力修正系数C-ap＝0.9时，阀门末端的压力曲线图；
图6为实施例下，压力修正系数C-ap＝0.5时，阀门末端的压力曲线图。
具体实施方式
下面结合具体实施例，进一步阐明本发明，应理解这些实施例仅用于说明本发明而不用于限制本发明的范围，在阅读了本发明之后，本领域技术人员对本发明的各种等价形式的修改均落于本申请所附权利要求所限定的范围。
基于有限体积法的输水管道中气-液两相均质流的模拟方法，按以下步骤进行：1.在非守恒形式的欧拉体系下，构建含自由气体的两相均质流基本控制方程，根据模拟工况确定计算域、初始条件以及边界条件；2.通过戈杜诺夫格式划分计算网格，并建立离散方程；3.通过含有Total Variation Diminishing(TVD) 形式的MUSCL-Hancoke格式，得到重建的流动变量，使其具有二阶精度；4.求解模型内部界面黎曼解，得到控制单元界面通量；5.根据已有的边界条件，通过牛顿迭代法利用黎曼不变量得到二阶的边界值；6.引入源项，通过基于二阶龙格库塔法离散格式的时间分裂法，求得时间梯度上二阶FVM-Godunov离散方程；7.建立稳定性约束条件，更新初始变量进行下一时步计算。图1给出了本发明的基本流程图。
下面将结合附图和实施例对本发明技术方案做进一步的详细描述。
实施例:
为了验证并分析本发明提供的非守恒格式有限体积法均质两相流模型的模拟效果，选取Chaudry于1990年设计搭建的均质流实验装置系统用于验证本发明方法的有效性。整个系统由上游水箱，管道，下游球阀，下游水箱组成。管道总长30.6m，内径0.026mm，管道水平。达西韦斯巴赫系数为0.0195，纯水水锤波速为715m/s。瞬变过程由突然关闭下游球阀引起。实验工况参数为：初始流速2.940m/s，上游水库静压水头21.700m，稳定时气体质量流动率1.15× 10-5kg/s，下游水体中气体体积含量0.0053。
本发明的具体步骤为：
步骤1：在非守恒形式的有限体积法体系下，构建含自由气体的瞬变流基本控制方程，根据工程实际确定计算域、初始条件以及边界条件。
(a)建立基本控制方程：
瞬变流基本微分方程：
Figure GDA0002982439880000091
Figure GDA0002982439880000092
其中，沿管线距离x与时间t是自变量；ρ(x,t)是平均截面流体密度；v(x,t)是平均截面速率；g是重力加速度；J为单位管长的摩擦力；ρl为水体的密度。
气-液两相流水锤波速与压力的关系方程：
Figure GDA0002982439880000093
其中，cm为混合流体密度；C为纯液体下水锤波速；Pref为参考压力；ρfref为参考压力下两相流的密度；ψref为参考压力下混合流体内气体的初始体积分数；
气-液两相流压力与密度的关系方程：
Figure GDA0002982439880000094
其中，ρ为混合流体的密度，P为流体的压力。
(b)确定计算域、初始条件及边界条件
计算域：从上游水库出水口至阀门之间的管道；
初始条件：球阀全开时，初始流速为2.940m/s，各节点初始压力根据上游水库静压水头减去相应的稳态摩阻损失得到；
边界条件为：管道入口处，水库提供恒压边界，且等于上游水库静压水头；下游球阀快速关闭引发空穴流瞬变，此处压力变化(图示4)由压力传感器3收集得到。
步骤2：在下建立均质两相流的计算网格并建立离散模型：
将空间域X离散为N个长度为Δx、时间域T上离散为间隔为Δt的控制单元，在空间处，第i个网格以i-1/2、i+1/2为边界，编号为i，用其代表该单元处流体参数的平均值。
为了方便模型计算，在模型两边边界处各建立两个虚拟单元，分别编号为-1、 0、N+1、N+2并假设虚拟单元处的流动参数与边界处的参数相同。
在管道内两相流离散模型系统基础上，建立离散方程的步骤为：
(a)将微分方程(1)(2)转化为拟线性的非守恒向量格式：
Figure GDA0002982439880000101
其中
Figure GDA00029824398800001011
(D为管道内径，f为达西韦斯巴赫系数)；
Figure GDA0002982439880000102
Figure GDA0002982439880000103
(P为管道内平均截面的压力)；F(W)为非守恒系统向量形式的通量。
(b)忽略初始常量，对方程(5)中的通量进行积分：
Figure GDA0002982439880000104
(c)在单元i(从界面i-1/2到i+1/2)及时间段Δt(从t到t+Δt)内对方程 (5)积分，得到流动变量W的离散方程：
Figure GDA0002982439880000105
其中，上标n，n+1/2和n+1分别代表t、t+1/2Δt和t+Δt时步；
Figure GDA0002982439880000106
为W在整个单元处流体参数的平均值；F为界面处的通量。
步骤3：通过使用TotalVariationDiminishing(TVD)版本的MUSCL-Hancock 格式，利用斜率限制器得到二阶精度同时避免了虚假震荡。
进一步的，步骤3包含以下子步骤：
步骤3.1：引入斜率限制器
Figure GDA0002982439880000107
其中，
Figure GDA0002982439880000108
为斜率限制器参数。
步骤3.2：非守恒变量的重建
Figure GDA0002982439880000109
步骤3.3：进一步重建
Figure GDA00029824398800001010
步骤3.4：黎曼问题
Figure GDA0002982439880000111
步骤4：求解黎曼问题，并通过黎曼解得到模型单元内部界面通量。
根据，对对任一控制体i(1≤i<N)，界面i+1/2处黎曼解的流动变量值为：
Figure GDA0002982439880000112
模型单元界面通量值为：
Figure GDA0002982439880000113
其中，
Figure GDA0002982439880000114
根据求出的变量ρ*的值，通过方程(4)，利用牛顿迭代法可求解得到方程(13)中的压力P。
步骤5：模型中二阶边界的求解。
本实施例中，包含两种边界条件：
(a)已知上游水库边界处的压力值
Figure GDA0002982439880000115
在边界处，沿着特征线dx/dt＝v-cm，满足
Figure GDA0002982439880000116
为了得到二阶精度的边界解，边界处引入一个虚拟网格，如图2所示，1为自由气体，在左侧单元内从(1/2+,tn)到(0,tn+1/2)对微分关系式
Figure GDA0002982439880000117
进行积分，得到：
Figure GDA0002982439880000118
其中，
Figure GDA0002982439880000119
和
Figure GDA00029824398800001110
为已知量，
Figure GDA00029824398800001111
为已知的边界条件，未知量
Figure GDA00029824398800001112
可以通过以下方程求解得出：
Figure GDA00029824398800001113
Figure GDA00029824398800001114
结合方程(14)，待求值
Figure GDA00029824398800001115
可以被推导为：
Figure GDA0002982439880000121
(b)已知下游水库边界处的流量值
Figure GDA0002982439880000122
模型右侧边界处的流速可由
Figure GDA0002982439880000123
得到，利用方程(3)，(4)通过牛顿迭代法进一步求得未知量
Figure GDA0002982439880000124
首先假设变量
Figure GDA0002982439880000125
(例如
Figure GDA0002982439880000126
)然后将该值带入方程(4)计算得到
Figure GDA0002982439880000127
将计算得到的
Figure GDA0002982439880000128
带入方程(3)计算得到
Figure GDA0002982439880000129
将这些参数带入方程(18)进一步得到
Figure GDA00029824398800001210
利用新的到的参数得到
Figure GDA00029824398800001211
重复以上步骤直到计算收敛。一旦得到
Figure GDA00029824398800001212
和
Figure GDA00029824398800001213
的值,
Figure GDA00029824398800001214
就可以利用公式计算得出。
Figure GDA00029824398800001215
步骤6:通过基于二阶Runge-Kutta离散格式的时间算子分裂法，引入源项，求得离散方程解的二阶显式戈杜诺夫格式。
具体实现过程为：
(a)纯对流时：
Figure GDA00029824398800001216
(b)利用源项乘以Δt/2进行更新：
Figure GDA00029824398800001217
(c)利用源项乘以Δt再次更新：
Figure GDA00029824398800001218
步骤7:建立有效的稳定性约束条件，更新初始变量并开展下一时步计算。
(a)由于采用了显式二阶龙格-库塔离散法将S引入求解，稳定性约束不仅要包括对流部分的Courant-Friedrichs-Lewy(CFL)准则，而且还要包括源项的约束。由CFL得到:
Figure GDA00029824398800001219
(b)显式二阶龙格-库塔离散化约束：
Figure GDA0002982439880000131
(c)由于对流项和源项采用相同的时间步长Δt,因而采用
Figure GDA0002982439880000132
而不是
Figure GDA0002982439880000133
最后给出包括对流部分和源项的最大允许时间步长:
Figure GDA0002982439880000134
通过以上方法编程计算后，将非守恒形式的模型计算结果与实验数据作对比。图4-6分别给出了本实施例下网格数N＝200，库朗特数Crmax＝0.95时，压力传感器1、2处的压力曲线图。可以看出，基于本模型下的预测值能够很好的拟合实验检测值，本发明的有效性以及准确性得到了很好的验证。同时，针对输水管道中气-液两相均质流，相比MOC类方法以及现有文献中的守恒形式的方法，本发明提出的非守恒有限体积戈杜诺夫格式，在模型预测方面更加的准确且更加的有效率。

Claims (6)

1.一种基于有限体积法的输水管道内气-液两相均质流的模拟方法，其特征在于，采用了二阶非守恒变量的有限体积法戈杜诺夫格式来模拟管道系统内均质两相流瞬变过程，具体步骤如下：
步骤1：在非守恒形式的欧拉体系下，构建含自由气体的两相均质流基本控制方程，根据模拟工况确定计算域、初始条件以及边界条件；
步骤1中，在非守恒形式的欧拉体系下，构建含自由气体的两相均质流基本控制方程，需在水锤问题的基础上假定：(a)管道内自由气体体积含量很低且气体在水体中均匀分散，气-液两相流以等效单相流体处理；(b)瞬变过程时间尺度很小，忽略整个过程中气体的吸收与释放；(c)管道内流体为无粘性流动且整个过程发生在等温条件下不考虑热传递；(d)管道和水体均为刚性，管道截面积在整个瞬变流过程中不发生变化；
步骤1中，构建的含自由气体均质两相流基本控制方程包含：
(1)水锤的基本方程：
Figure FDA0002982439870000011
Figure FDA0002982439870000012
其中，沿管线距离x与时间t是自变量；ρ(x,t)是平均截面流体密度；v(x,t)是平均截面速率；g是重力加速度；J为单位管长的摩擦力；ρl为水体的密度；
(2)气-液两相流水锤波速与压力的关系：
Figure FDA0002982439870000013
其中，cm为混合流体密度；C为纯液体下水锤波速；Pref为参考压力；ρfref为参考压力下两相流的密度；ψref为参考压力下混合流体内气体的初始体积分数；
(3)气-液两相流压力与密度的关系：
Figure FDA0002982439870000014
其中，ρ为混合流体的密度，P为管道内平均截面的压力；
步骤2：通过有限体积法戈杜诺夫格式划分计算网格，并建立离散方程；
步骤2中，在Godunov格式下建立均质两相流的计算网格的方法为：
(a)建立初始网格：
将空间域X离散为N个长度为Δx、时间域T上离散为间隔为Δt的控制单元，在空间处，第i个网格以i-1/2、i+1/2为边界，编号为i，用其代表该单元处流体参数的平均值；
(b)建立虚拟网格：
为了方便模型计算，在模型两边边界处各建立两个虚拟单元，分别编号为-1、0、N+1、N+2并假设虚拟单元处的流动参数与边界处的参数相同；
步骤2中，在管道内两相流离散模型系统基础上，建立离散方程的步骤为：
(a)将微分方程(1)(2)转化为拟线性的非守恒向量格式：
Figure FDA0002982439870000021
其中
Figure FDA0002982439870000022
D为管道内径，f为达西韦斯巴赫系数；
Figure FDA0002982439870000023
Figure FDA0002982439870000024
P为管道内平均截面的压力；F(W)为非守恒系统向量形式的通量；
(b)忽略初始常量，对方程(5)中的通量进行积分：
Figure FDA0002982439870000025
(c)在从界面i-1/2到i+1/2的单元i及从t到t+Δt的时间段Δt内对方程(5)积分，得到流动变量W的离散方程：
Figure FDA0002982439870000026
其中，上标n，n+1/2和n+1分别代表t、t+1/2Δt和t+Δt时步；
Figure FDA0002982439870000027
Figure FDA0002982439870000028
为W在整个单元处流体参数的平均值；F为界面处的通量；
步骤3：通过含有TVD形式的MUSCL-Hancoke格式，得到重建的流动变量，使其具有二阶精度；
步骤4：求解模型内部界面黎曼解，得到控制单元界面通量；
步骤5：根据已有的边界条件，通过牛顿迭代法利用黎曼不变量得到二阶的边界值；
步骤6：引入源项，通过基于二阶龙格库塔法离散格式的时间分裂法，求得时间梯度上二阶FVM-Godunov离散方程；
步骤7：建立稳定性约束条件，更新初始变量进行下一时步计算。
2.根据权利要求1所述的基于有限体积法的输水管道内气-液两相均质流的模拟方法，其特征在于，步骤3中，通过使用TVD版本的MUSCL-Hancock格式，利用斜率限制器得到二阶精度的格式的方法为：
(a)引入斜率限制器
Figure FDA0002982439870000031
其中，
Figure FDA0002982439870000032
为斜率限制器参数；
(b)非守恒变量的重建
Figure FDA0002982439870000033
(c)进一步重建
Figure FDA0002982439870000034
(d)黎曼问题
Figure FDA0002982439870000035
3.根据权利要求1所述的基于有限体积法的输水管道内气-液两相均质流的模拟方法，其特征在于，步骤4中，求解黎曼问题，并通过黎曼解得到模型内部界面通量的方法为：
(a)一般的双曲系统的黎曼问题也就是初始值问题：
Wt+Fx＝0 (12)
Figure FDA0002982439870000036
这里
Figure FDA0002982439870000041
和
Figure FDA0002982439870000042
均为常数值，F为界面处通量；对方程(6)在单元处进行积分，得到：
Figure FDA0002982439870000043
(b)在有限体积法中，穿过断点处，利用Rankine-Hugoniot条件ΔF＝λΔW，得到：
Figure FDA0002982439870000044
(c)利用压力与波速的一般方程
Figure FDA0002982439870000045
推导出：
Figure FDA0002982439870000046
(d)通过一阶的线性近似：
Figure FDA0002982439870000047
其中，
Figure FDA0002982439870000048
并且当网格足够精细时，线性近似足够保证精度；
(e)对微分方程(17)在单元处进行积分：
Figure FDA0002982439870000049
(f)求解黎曼问题中中间处参数值：
Figure FDA00029824398700000410
(g)求解单元界面通量值：
Figure FDA00029824398700000411
其中，
Figure FDA00029824398700000412
根据求出的变量ρ*的值，通过方程(4)，利用牛顿迭代法可求解得到方程(20)中的压力P。
4.根据权利要求1所述的基于有限体积法的输水管道内气-液两相均质流的模拟方法，其特征在于，步骤5中，模型中二阶边界的求解方法为：
(a)以左侧为例，当已知边界处的压力值
Figure FDA0002982439870000051
时：
在边界处，沿着特征线dx/dt＝v-cm，满足
Figure FDA0002982439870000052
为了得到二阶精度的边界解，边界处引入一个虚拟网格，在左侧单元内从(1/2+,tn)到(0,tn+1/2)对微分关系式
Figure FDA0002982439870000053
进行积分，得到：
Figure FDA0002982439870000054
其中，
Figure FDA0002982439870000055
和
Figure FDA0002982439870000056
为已知量，
Figure FDA0002982439870000057
为已知的边界条件，未知量
Figure FDA0002982439870000058
可以通过以下方程求解得出：
Figure FDA0002982439870000059
Figure FDA00029824398700000510
结合方程(21)，待求值
Figure FDA00029824398700000511
可以被推导为：
Figure FDA00029824398700000512
(b)以模型左侧为例，当已知边界处的流量值
Figure FDA00029824398700000513
时：
模型左侧边界处的流速可由
Figure FDA00029824398700000514
得到，利用方程(3)，(4)通过牛顿迭代法进一步求得未知量
Figure FDA00029824398700000515
首先假设变量
Figure FDA00029824398700000516
然后将该值带入方程(4)计算得到
Figure FDA00029824398700000517
将计算得到的
Figure FDA00029824398700000518
带入方程(3)计算得到
Figure FDA00029824398700000519
将这些参数带入方程(25)进一步得到
Figure FDA00029824398700000520
利用新的到的参数得到
Figure FDA00029824398700000521
重复以上步骤直到计算收敛；一旦得到
Figure FDA00029824398700000522
和
Figure FDA00029824398700000523
的值,
Figure FDA00029824398700000524
就可以利用公式计算得出：
Figure FDA00029824398700000525
5.根据权利要求1所述的基于有限体积法的输水管道内气-液两相均质流的模拟方法，其特征在于，步骤6中，引入源项通过二阶龙格库塔法求解离散方程：
Figure FDA0002982439870000061
其中，
Figure FDA0002982439870000062
为n+1时步，控制单元i在纯对流时，流动变量W的通量；
Figure FDA0002982439870000063
为采用时间分裂法第一次更新后的通量。
6.根据权利要求1所述的基于有限体积法的输水管道内气-液两相均质流的模拟方法，其特征在于，步骤7中，建立有效的稳定性约束条件，更新初始变量并开展下一时步计算的方法为：
(a)由于采用了显式二阶龙格-库塔离散法将S引入求解，稳定性约束不仅要包括对流部分的CFL准则，而且还要包括源项的约束；由CFL得到:
Figure FDA0002982439870000064
(b)显式二阶龙格-库塔离散化约束：
Figure FDA0002982439870000065
(c)由于对流项和源项采用相同的时间步长Δt,因而采用
Figure FDA0002982439870000066
而不是
Figure FDA0002982439870000067
最后给出包括对流部分和源项的最大允许时间步长:
Figure FDA0002982439870000068
CN201910173947.5A 2019-03-08 2019-03-08 基于有限体积法的输水管道内气液两相均质流的模拟方法 Active CN109918787B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910173947.5A CN109918787B (zh) 2019-03-08 2019-03-08 基于有限体积法的输水管道内气液两相均质流的模拟方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910173947.5A CN109918787B (zh) 2019-03-08 2019-03-08 基于有限体积法的输水管道内气液两相均质流的模拟方法

Publications (2)

Publication Number Publication Date
CN109918787A CN109918787A (zh) 2019-06-21
CN109918787B true CN109918787B (zh) 2021-05-11

Family

ID=66963810

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910173947.5A Active CN109918787B (zh) 2019-03-08 2019-03-08 基于有限体积法的输水管道内气液两相均质流的模拟方法

Country Status (1)

Country Link
CN (1) CN109918787B (zh)

Families Citing this family (12)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110569541A (zh) * 2019-08-01 2019-12-13 天津大学 管道水锤分析方法
CN110633502B (zh) * 2019-08-20 2021-02-02 中南大学 一种考虑织物透气性的超声速降落伞数值模拟方法
CN111046567B (zh) * 2019-12-18 2020-07-31 中国水利水电科学研究院 一种基于Godunov格式的城市排水管网水流数值模拟方法
CN111414679B (zh) * 2020-03-12 2022-06-03 河海大学 一种岩塞爆破水气过渡过程水力特性的计算方法
CN111414683B (zh) * 2020-03-16 2021-10-29 河海大学 一种考虑动态摩阻的水气耦合瞬变流的模拟方法
CN112861263B (zh) * 2021-02-22 2024-02-13 西北工业大学 一种适用于可压缩两相流的计算模拟方法
CN113094917A (zh) * 2021-04-21 2021-07-09 电子科技大学成都学院 一种高压油管的单向阀开启计算方法
CN113435136B (zh) * 2021-06-24 2023-06-20 河海大学 耦合能量方程的输水管道气-液两相均质流的模拟方法
CN113361217B (zh) * 2021-07-07 2022-10-11 中国海洋大学 一种高效的两相流无网格数值模型实施方法、装置
CN113656926B (zh) * 2021-08-26 2024-03-26 河海大学 基于Schohl卷积近似的管道瞬变流模拟方法
CN114254572B (zh) * 2021-12-16 2024-01-02 西北工业大学太仓长三角研究院 考虑污染物沉积的航发压气机流场性能预测方法及系统
CN114896908B (zh) * 2022-05-19 2023-03-28 四川大学 一种基于通量限制器的水滴流场及水滴收集率计算方法

Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5880378A (en) * 1996-08-19 1999-03-09 Southwest Research Institute Critical flow venturi with variable and continuous range
US20080103742A1 (en) * 2005-08-17 2008-05-01 Jiun-Der Yu Coupled Algorithms on Quadrilateral Grids for Generalized Axi-Symmetric Viscoelastic Fluid Flows
CN103226641A (zh) * 2013-05-10 2013-07-31 中国石油大学(华东) 深水气液两相流循环温度压力耦合计算方法
CN105512363A (zh) * 2015-11-25 2016-04-20 河海大学 基于Godunov格式的有压管道中水柱分离的模拟方法
CN106503396A (zh) * 2016-11-14 2017-03-15 中国电建集团昆明勘测设计研究院有限公司 基于有限差分法与有限体积法耦合的多维水力系统瞬变模拟方法
CN106777770A (zh) * 2017-01-09 2017-05-31 河海大学 基于有限体积法的输水管道中空穴流的模拟方法
CN107038295A (zh) * 2017-04-06 2017-08-11 中国水利水电科学研究院 一种水锤泵内部流道评价及优化方法

Patent Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5880378A (en) * 1996-08-19 1999-03-09 Southwest Research Institute Critical flow venturi with variable and continuous range
US20080103742A1 (en) * 2005-08-17 2008-05-01 Jiun-Der Yu Coupled Algorithms on Quadrilateral Grids for Generalized Axi-Symmetric Viscoelastic Fluid Flows
CN103226641A (zh) * 2013-05-10 2013-07-31 中国石油大学(华东) 深水气液两相流循环温度压力耦合计算方法
CN105512363A (zh) * 2015-11-25 2016-04-20 河海大学 基于Godunov格式的有压管道中水柱分离的模拟方法
CN106503396A (zh) * 2016-11-14 2017-03-15 中国电建集团昆明勘测设计研究院有限公司 基于有限差分法与有限体积法耦合的多维水力系统瞬变模拟方法
CN106777770A (zh) * 2017-01-09 2017-05-31 河海大学 基于有限体积法的输水管道中空穴流的模拟方法
CN107038295A (zh) * 2017-04-06 2017-08-11 中国水利水电科学研究院 一种水锤泵内部流道评价及优化方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
Godunov’s Method for Simulatinons of Fluid-Structure Interaction in Piping Systems;Janez Gale;《Journal of Pressure Vessel Technology》;20080831;全文 *
基于有限体积法 Godunov 格式的水锤计算模型;赵越;《水利水电科技进展》;20190131;全文 *

Also Published As

Publication number Publication date
CN109918787A (zh) 2019-06-21

Similar Documents

Publication Publication Date Title
CN109918787B (zh) 基于有限体积法的输水管道内气液两相均质流的模拟方法
Zhou et al. Investigation of hydraulic transients of two entrapped air pockets in a water pipeline
Fuertes-Miquel et al. Transient phenomena during the emptying process of a single pipe with water–air interaction
CN106844913B (zh) 一种基于三维cfd的滞留气团热力学特性模拟方法
Zhou et al. Phenomenon of white mist in pipelines rapidly filling with water with entrapped air pockets
Wang et al. Water hammer simulation using explicit–implicit coupling methods
Leon Improved modeling of unsteady free surface, pressurized and mixed flows in storm-sewer systems
CN105468844A (zh) 管道内水-气耦合瞬变流的模拟方法
Yang et al. Dynamic analysis of the pump system based on MOC–CFD coupled method
Zhou et al. An accurate and efficient scheme involving unsteady friction for transient pipe flow
Vasconcelos et al. Innovative simulation of unsteady low-pressure flows in water mains
Ferrari et al. A 5‐equation, transient, hyperbolic, 1‐dimensional model for slug capturing in pipes
Tamhankar et al. Experimental and CFD analysis of flow through venturimeter to determine the coefficient of discharge
Zhou et al. Energy dissipation in a rapid filling vertical pipe with trapped air
Modisette Pipeline thermal models
CN106777770B (zh) 基于有限体积法的输水管道中空穴流的模拟方法
Chao et al. Study on the algorithm for solving two-fluid seven-equation two-pressure model
CN111695307B (zh) 显式考虑动态摩阻的水锤有限体积模拟方法
CN113435136B (zh) 耦合能量方程的输水管道气-液两相均质流的模拟方法
Li et al. 1D-3D coupling investigation of hydraulic transient for power-supply failure of centrifugal pump-pipe system
Cheng et al. Simulation of complex seepage field of a gravity dam foundation using a CFD-based approach
Liu et al. A discrete unified gas kinetic scheme for simulating transient hydrodynamics in porous media with fractures
CN109670265B (zh) 高压加热器建模方法、装置和计算机设备
CN105973319A (zh) 一种控制棒驱动机构排污系统水力特性计算方法
Hou et al. Discussion of “rigid water column model for simulating the emptying process in a pipeline using pressurized air” by Oscar E. Coronado-Hernández, Vicente S. Fuertes-Miquel, Pedro L. Iglesias-Rey, and Francisco J. Martínez-Solano

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