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

Abstract

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

Description

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

技术领域

本发明涉及一种基于有限体积法的输水管道内水气两相均质流的模拟方法，属于水电站(泵站)水力学数值模拟计算技术领域。

背景技术

在输水管道系统中，阀门或机组的突然启闭可能导致管道内的压力突变，当管道内的压力变化超过管壁材料的承受能力时，可能会造成严重的管道系统破坏，甚至会威胁人身安全。近些年，管道内的两相流问题一直是国内外研究的热门课题，不同于一般的单相流问题，两相流需要考虑不同相间的耦合作用。在两相流问题中，均质流模型被认为是应用最广的模型。在输水管道系统中，若有少量的自由气体均匀的分布在水体中，且假设水气两相不会随着流动发生相互移动，这样的模型被成为两相均质流模型。

目前，针对输水管道内的气-液两相均质流问题，其模拟的方法主要为特征线法(MOC,Method of Characteristics)以及有限差分法(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)水锤的基本方程：

其中，沿管线距离x与时间t是自变量；ρ(x,t)是平均截面流体密度；v(x,t)是平均截面速率；g是重力加速度；J为单位管长的摩擦力；ρ_l为水体的密度；

(2)气-液两相流水锤波速与压力的关系：

其中，c_m为混合流体密度；C为纯液体下水锤波速；P_ref为参考压力；ρ_fref为参考压力下两相流的密度；ψ_ref为参考压力下混合流体内气体的初始体积分数；

(3)气-液两相流压力与密度的关系：

其中，ρ为混合流体的密度，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)转化为拟线性的非守恒向量格式：

其中 D为管道内径，f为达西韦斯巴赫系数； P为管道内平均截面的压力；F(W)为非守恒系统向量形式的通量；

(b)忽略初始常量，对方程(5)中的通量进行积分：

(c)在从界面i-1/2到i+1/2的单元i及从t到t+Δt的时间段Δt内对方程(5)积分，得到流动变量W的离散方程：

其中，上标n，n+1/2和n+1分别代表t、t+1/2Δt和t+Δt时步；为W在整个单元处流体参数的平均值；F为界面处的通量。

进一步地，步骤3中，通过使用TVD版本的MUSCL-Hancock格式，利用斜率限制器得到二阶精度的格式的方法为：

(a)引入斜率限制器

其中，为斜率限制器参数；

(b)非守恒变量的重建

(c)进一步重建

(d)黎曼问题

进一步地，步骤4中，求解黎曼问题，并通过黎曼解得到模型内部界面通量的方法为：

(a)一般的双曲系统的黎曼问题也就是初始值问题：

W_t+F_x＝0 (12)

这里和均为常数值，F为界面处通量；对方程(6)在单元处进行积分，得到：

(b)在有限体积法中，穿过断点处，利用Rankine-Hugoniot条件ΔF＝λΔW，得到：

(c)利用压力与波速的一般方程推导出：

(d)通过一阶的线性近似：

其中，并且当网格足够精细时，线性近似足够保证精度；

(e)对微分方程(17)在单元处进行积分：

(f)求解黎曼问题中中间处参数值：

(g)求解单元界面通量值：

其中，根据求出的变量ρ_*的值，通过方程(4)，利用牛顿迭代法可求解得到方程(20)中的压力P。

进一步地，步骤5中，模型中二阶边界的求解方法为：

(a)以左侧为例，当已知边界处的压力值时：

在边界处，沿着特征线dx/dt＝v-c_m，满足为了得到二阶精度的边界解，边界处引入一个虚拟网格，在左侧单元内从(1/2+,tⁿ)到(0,t^n+1/2)对微分关系式进行积分，得到：

其中，和为已知量，为已知的边界条件，未知量可以通过以下方程求解得出：

结合方程(21)，待求值可以被推导为：

(b)以模型左侧为例，当已知边界处的流量值时：

模型左侧边界处的流速可由得到，利用方程(3)，(4)通过牛顿迭代法进一步求得未知量首先假设变量 然后将该值带入方程(4)计算得到将计算得到的带入方程(3)计算得到将这些参数带入方程(25)进一步得到利用新的到的参数得到重复以上步骤直到计算收敛；一旦得到和的值,就可以利用公式计算得出：

进一步地，步骤6中，引入源项通过二阶龙格库塔法求解离散方程：

其中，为n+1时步，控制单元i在纯对流时，流动变量W的通量；为采用时间分裂法第一次更新后的通量。

进一步地，步骤7中，建立有效的稳定性约束条件，更新初始变量并开展下一时步计算的方法为：

(a)由于采用了显式二阶龙格-库塔离散法将S引入求解，稳定性约束不仅要包括对流部分的CFL准则，而且还要包括源项的约束；由CFL得到:

(b)显式二阶龙格-库塔离散化约束：

(c)由于对流项和源项采用相同的时间步长Δt,因而采用而不是最后给出包括对流部分和源项的最大允许时间步长:

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

(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)建立基本控制方程：

瞬变流基本微分方程：

其中，沿管线距离x与时间t是自变量；ρ(x,t)是平均截面流体密度；v(x,t)是平均截面速率；g是重力加速度；J为单位管长的摩擦力；ρ_l为水体的密度。

气-液两相流水锤波速与压力的关系方程：

其中，c_m为混合流体密度；C为纯液体下水锤波速；P_ref为参考压力；ρ_fref为参考压力下两相流的密度；ψ_ref为参考压力下混合流体内气体的初始体积分数；

气-液两相流压力与密度的关系方程：

其中，ρ为混合流体的密度，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)转化为拟线性的非守恒向量格式：

其中 (D为管道内径，f为达西韦斯巴赫系数)； (P为管道内平均截面的压力)；F(W)为非守恒系统向量形式的通量。

(b)忽略初始常量，对方程(5)中的通量进行积分：

(c)在单元i(从界面i-1/2到i+1/2)及时间段Δt(从t到t+Δt)内对方程(5)积分，得到流动变量W的离散方程：

其中，上标n，n+1/2和n+1分别代表t、t+1/2Δt和t+Δt时步；为W在整个单元处流体参数的平均值；F为界面处的通量。

步骤3：通过使用Total Variation Diminishing(TVD)版本的MUSCL-Hancock格式，利用斜率限制器得到二阶精度同时避免了虚假震荡。

进一步的，步骤3包含以下子步骤：

步骤3.1：引入斜率限制器

其中，为斜率限制器参数。

步骤3.2：非守恒变量的重建

步骤3.3：进一步重建

步骤3.4：黎曼问题

步骤4：求解黎曼问题，并通过黎曼解得到模型单元内部界面通量。

根据，对对任一控制体i(1≤i<N)，界面i+1/2处黎曼解的流动变量值为：

模型单元界面通量值为：

其中，根据求出的变量ρ_*的值，通过方程(4)，利用牛顿迭代法可求解得到方程(13)中的压力P。

步骤5：模型中二阶边界的求解。

本实施例中，包含两种边界条件：

(a)已知上游水库边界处的压力值

在边界处，沿着特征线dx/dt＝v-c_m，满足为了得到二阶精度的边界解，边界处引入一个虚拟网格，如图2所示，1为自由气体，在左侧单元内从(1/2⁺,tⁿ)到(0,t^n+1/2)对微分关系式进行积分，得到：

其中，和为已知量，为已知的边界条件，未知量可以通过以下方程求解得出：

结合方程(14)，待求值可以被推导为：

(b)已知下游水库边界处的流量值

模型右侧边界处的流速可由得到，利用方程(3)，(4)通过牛顿迭代法进一步求得未知量首先假设变量(例如)然后将该值带入方程(4)计算得到将计算得到的带入方程(3)计算得到将这些参数带入方程(18)进一步得到利用新的到的参数得到重复以上步骤直到计算收敛。一旦得到和的值,就可以利用公式计算得出。

步骤6:通过基于二阶Runge-Kutta离散格式的时间算子分裂法，引入源项，求得离散方程解的二阶显式戈杜诺夫格式。

具体实现过程为：

(a)纯对流时：

(b)利用源项乘以Δt/2进行更新：

(c)利用源项乘以Δt再次更新：

步骤7:建立有效的稳定性约束条件，更新初始变量并开展下一时步计算。

(a)由于采用了显式二阶龙格-库塔离散法将S引入求解，稳定性约束不仅要包括对流部分的Courant-Friedrichs-Lewy(CFL)准则，而且还要包括源项的约束。由CFL得到:

(b)显式二阶龙格-库塔离散化约束：

(c)由于对流项和源项采用相同的时间步长Δt,因而采用而不是最后给出包括对流部分和源项的最大允许时间步长:

通过以上方法编程计算后，将非守恒形式的模型计算结果与实验数据作对比。图4-6分别给出了本实施例下网格数N＝200，库朗特数Cr_max＝0.95时，压力传感器1、2处的压力曲线图。可以看出，基于本模型下的预测值能够很好的拟合实验检测值，本发明的有效性以及准确性得到了很好的验证。同时，针对输水管道中气-液两相均质流，相比MOC类方法以及现有文献中的守恒形式的方法，本发明提出的非守恒有限体积戈杜诺夫格式，在模型预测方面更加的准确且更加的有效率。

Claims (10)

1.一种基于有限体积法的输水管道内水气两相均质流的模拟方法，其特征在于，采用了二阶非守恒变量的有限体积法戈杜诺夫格式来模拟管道系统内均质两相流瞬变过程，具体步骤如下：

步骤1：在非守恒形式的欧拉体系下，构建含自由气体的两相均质流基本控制方程，根据模拟工况确定计算域、初始条件以及边界条件；

步骤2：通过有限体积法戈杜诺夫格式划分计算网格，并建立离散方程；

步骤3：通过含有TVD形式的MUSCL-Hancoke格式，得到重建的流动变量，使其具有二阶精度；

步骤4：求解模型内部界面黎曼解，得到控制单元界面通量；

步骤5：根据已有的边界条件，通过牛顿迭代法利用黎曼不变量得到二阶的边界值；

步骤6：引入源项，通过基于二阶龙格库塔法离散格式的时间分裂法，求得时间梯度上二阶FVM-Godunov离散方程；

步骤7：建立稳定性约束条件，更新初始变量进行下一时步计算。

2.根据权利要求1所述的基于有限体积法戈杜诺夫格式的输水管道中水气两相均质流的模拟方法，其特征在于，步骤1中，在非守恒形式的欧拉体系下，构建含自由气体的两相均质流基本控制方程，需在水锤问题的基础上假定：(a)管道内自由气体体积含量很低(<≈1％)且气体在水体中均匀分散，气-水两相流以等效单相流体处理；(b)瞬变过程时间尺度很小，忽略整个过程中气体的吸收与释放；(c)管道内流体为无粘性流动且整个过程发生在等温条件下不考虑热传递；(d)管道和水体均为刚性，管道截面积在整个瞬变流过程中不发生变化。

3.根据权利要求1所述的基于有限体积法的输水管道中气-液两相均质流的模拟方法，其特征在于，步骤1中，构建的含自由气体均质两相流基本控制方程包含：

(1)水锤的基本方程：

其中，沿管线距离x与时间t是自变量；ρ(x,t)是平均截面流体密度；v(x,t)是平均截面速率；g是重力加速度；J为单位管长的摩擦力；ρ_l为水体的密度；

(2)气-液两相流水锤波速与压力的关系：

其中，c_m为混合流体密度；C为纯液体下水锤波速；P_ref为参考压力；ρ_fref为参考压力下两相流的密度；ψ_ref为参考压力下混合流体内气体的初始体积分数；

(3)气-液两相流压力与密度的关系：

其中，ρ为混合流体的密度，P为流体的压力。

4.根据权利要求1所述的基于有限体积法的输水管道内水气两相均质流的模拟方法，其特征在于，步骤2中，在Godunov格式下建立均质两相流的计算网格的方法为：

(a)建立初始网格：

将空间域X离散为N个长度为Δx、时间域T上离散为间隔为Δt的控制单元，在空间处，第i个网格以i-1/2、i+1/2为边界，编号为i，用其代表该单元处流体参数的平均值；

(b)建立虚拟网格：

为了方便模型计算，在模型两边边界处各建立两个虚拟单元，分别编号为-1、0、N+1、N+2并假设虚拟单元处的流动参数与边界处的参数相同。

5.根据权利要求3所述的基于有限体积法的输水管道内水气两相均质流的模拟方法，其特征在于，步骤2中，在管道内两相流离散模型系统基础上，建立离散方程的步骤为：

(a)将微分方程(1)(2)转化为拟线性的非守恒向量格式：

其中D为管道内径，f为达西韦斯巴赫系数； P为管道内平均截面的压力；F(W)为非守恒系统向量形式的通量；

(b)忽略初始常量，对方程(5)中的通量进行积分：

(c)在从界面i-1/2到i+1/2的单元i及从t到t+Δt的时间段Δt内对方程(5)积分，得到流动变量W的离散方程：

其中，上标n，n+1/2和n+1分别代表t、t+1/2Δt和t+Δt时步；为W在整个单元处流体参数的平均值；F为界面处的通量。

6.根据权利要求1所述的基基于有限体积法的输水管道内水气两相均质流的模拟方法，其特征在于，步骤3中，通过使用TVD版本的MUSCL-Hancock格式，利用斜率限制器得到二阶精度的格式的方法为：

(a)引入斜率限制器

其中，为斜率限制器参数；

(b)非守恒变量的重建

(c)进一步重建

(d)黎曼问题

7.根据权利要求1所述的基于基于有限体积法的输水管道内水气两相均质流的模拟方法，其特征在于，步骤4中，求解黎曼问题，并通过黎曼解得到模型内部界面通量的方法为：

(a)一般的双曲系统的黎曼问题也就是初始值问题：

W_t+F_x＝0 (12)

这里和均为常数值，F为界面处通量；对方程(6)在单元处进行积分，得到：

(b)在有限体积法中，穿过断点处，利用Rankine-Hugoniot条件ΔF＝λΔW，得到：

(c)利用压力与波速的一般方程推导出：

(d)通过一阶的线性近似：

其中，并且当网格足够精细时，线性近似足够保证精度；

(e)对微分方程(17)在单元处进行积分：

(f)求解黎曼问题中中间处参数值：

(g)求解单元界面通量值：

其中，根据求出的变量ρ_*的值，通过方程(4)，利用牛顿迭代法可求解得到方程(20)中的压力P。

8.根据权利要求1所述的基于有限体积法的输水管道中气-液两相均质流的模拟方法，其特征在于，步骤5中，模型中二阶边界的求解方法为：

(a)以左侧为例，当已知边界处的压力值时：

在边界处，沿着特征线dx/dt＝v-c_m，满足为了得到二阶精度的边界解，边界处引入一个虚拟网格，在左侧单元内从(1/2⁺,tⁿ)到(0,t^n+1/2)对微分关系式进行积分，得到：

其中，和为已知量，为已知的边界条件，未知量可以通过以下方程求解得出：

结合方程(21)，待求值可以被推导为：

(b)以模型左侧为例，当已知边界处的流量值时：

模型左侧边界处的流速可由得到，利用方程(3)，(4)通过牛顿迭代法进一步求得未知量首先假设变量然后将该值带入方程(4)计算得到将计算得到的带入方程(3)计算得到将这些参数带入方程(25)进一步得到利用新的到的参数得到重复以上步骤直到计算收敛；一旦得到和的值,就可以利用公式计算得出：

9.根据权利要求1所述的基于有限体积法的输水管道中气-液两相均质流的模拟方法，其特征在于，步骤6中，引入源项通过二阶龙格库塔法求解离散方程：

其中，为n+1时步，控制单元i在纯对流时，流动变量W的通量；为采用时间分裂法第一次更新后的通量。

10.根据权利要求1所述的基于有限体积法的输水管道中气-液两相均质流的模拟方法，其特征在于，步骤7中，建立有效的稳定性约束条件，更新初始变量并开展下一时步计算的方法为：

(a)由于采用了显式二阶龙格-库塔离散法将S引入求解，稳定性约束不仅要包括对流部分的CFL准则，而且还要包括源项的约束；由CFL得到:

(b)显式二阶龙格-库塔离散化约束：

(c)由于对流项和源项采用相同的时间步长Δt,因而采用而不是最后给出包括对流部分和源项的最大允许时间步长: