CN105512363B - 基于Godunov格式的有压管道中水柱分离的模拟方法 - Google Patents

Abstract

本发明公开了一种基于Godunov格式的有压管道中水柱分离的模拟方法，该方法首先将计算域离散成计算网格单元，建立瞬变流基本微分方程的离散格式；然后采用时空均为二阶精度的Godunov格式求解纯对流时，单元控制体的数值通量；接着，基于二阶Runge‑Kutta离散格式的时间算子分裂法，在纯对流控制方程的解中引入源项，从而得到最终解的二阶显式有限体积法(FVM,Finite Volume Method)Godunov格式；最后，根据算得的压力判定是否发生水柱分离：若已形成，则进行水柱分离模块计算；反之，进入下一时步计算。此外，本发明通过引入斜率限制器来抑制虚假的数值振荡。本发明可对输水管道系统中瞬变蒸汽空化现象进行准确、严格的数值模拟和瞬变分析，从而达到规避工程风险，保护水利系统的安全的目的。

Description

基于Godunov格式的有压管道中水柱分离的模拟方法

技术领域

本发明涉及一种基于Godunov格式的有压管道中水柱分离的模拟方法，属于水电站(泵站)水力学数值模拟计算技术领域。

背景技术

当输水管道系统中液体压力降低到汽化压力时，可能发生瞬变蒸汽空化现象。这些空化泡可能发育长大到占据整个管道截面的程度，破坏液柱的连续性。Bergant和Simpson等学者在文献中指出分离的液柱再弥合可能产生突升的高压脉冲，导致水力系统破坏。因此，为避免水柱分离现象，需对该现象进行严格的瞬变分析。

几十年来，很多学者采用各种方法来预测水柱分离过程中的瞬变压力，应用最广的是由Streeter和Wylie基于特征线法(MOC，Method of Characteristics)提出的离散蒸汽—空穴模型(DVCM,discrete vapor-cavity model)，该模型允许汽穴在压力达到或低于汽化压力的任何计算节点形成，并且作为一内部边界占据整个管道截面，因该模型易于嵌入到标准的水锤商业包中，且再现了水柱分离过程中很多的物理特性，成为目前模拟水柱分离现象最常用的模型。

经典MOC-DVCM的主要缺点是当采用相对细密的网格时，计算结果会出现不真实的压力脉冲。Chaiko在FVM方法基础上采用均质连续模型来描述水柱分离过程中的汽穴过程，然而，其算法相比MOC而言过于复杂，实现困难。因此，如果能够合理保留各种方法的优点，引入新的方法和技巧消除缺陷，就能有助于水柱分离过程中瞬变压力的数值模拟方法的不断改进。

发明内容

发明目的：为弥补现有MOC-DVCM方法存在的虚假震荡等不足，本发明基于FVM，提供了一种算法简单，易于实现的模拟方法，可以满足计算网格细密时不会出现人为的脉冲峰值。

技术方案：一种基于Godunov格式的有压管道中水柱分离的模拟方法，具体步骤如下：

步骤1：建立瞬变流基本微分方程，根据工程实例确定计算域、初始条件以及边界条件；

步骤2：根据FVM划分计算网格，并建立离散方程；

步骤3：采用Godunov方法求解纯对流时的离散方程，得到纯对流时控制单元界面处的数值通量，并取得二阶精度；

步骤4：通过基于二阶Runge-Kutta离散格式的时间算子分裂法，在纯对流控制方程的解中引入源项，从而得到最终解的二阶显式FVM-Godunov格式；

步骤5：根据计算得到的压力判定是否形成水柱分离：若已形成，则进行水柱分离模块计算；反之，进入下一时步计算。

进一步地，步骤2中，将空间域x离散为长度为Δx，时间域t离散为间隔为Δt的控制单元，数量为N。对第i个控制单元，定义其上、下游界面编号分别为i-1/2、i+1/2。

进一步地，步骤2中，对控制单元i，建立的流动变量u的积分方程为：

其中，上标n和n+1分别代表t和t+Δt时步；为u在整个控制体的平均值；H是测压管水头，V是平均截面速率；f为单元界面处的通量；为源项；f为达西-威斯巴哈摩阻系数；D为管径。

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

步骤3.1：求解内部控制单元界面处通量。

首先，基于黎曼问题，根据Godunov格式，对任一内部控制单元i(1<i<N)，界面i+1/2处的通量为：

其中，H是测压管水头；V是平均截面速率；是V的平均值，为一常数；为在n时步时，u分别到界面i+1/2左、右侧两侧的平均值。

接着，通过引入MUSCL-Hancock格式计算内部单元通量f_i+1/2，从而取得空间和时间上的二阶精度。

更进一步地，步骤3.1中，在引入MUSCL-Hancock格式计算内部单元通量f_i+1/2过程中，需选择斜率限制器，以保证解中不出现虚假振荡。

步骤3.2：求解边界控制单元界面处通量。为在边界面处也取得二阶精度，分别在起始控制单元1上游侧、终点控制单元N下游侧构建两个虚拟控制单元I_-1、I₀，以及I_N+1、I_N+2，并假定在虚拟单元处的流动信息与边界处是一致的。从而可求解边界黎曼问题，且相应的Godunov通量f_1/2和f_N+1/2也可像内部单元那样进行计算。

更进一步地，步骤3中，纯对流时离散方程中的对流项满足CFL条件(Courant-Friedrichs-Lewy criterion)，进一步可推得CFL条件下的最大时间步长Δt_max,CFL：

其中，C_r为柯朗数，为矩阵的特征值。

进一步地，步骤4中，引入源项后，瞬变流基本微分方程解的二阶FVM-Godunov格式为：

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

更进一步地，步骤4中，源项满足以下稳定性约束，并可推得适用于源项的最大时间步长Δt_max,s：

更进一步地，包含对流项和源项的最大允许时间步长为：

Δt_max＝min(Δt_max,CFL,Δt_max,s) (6)

进一步地，步骤5中，水柱分离计算模块中，假定汽穴仅在控制体的中部形成，为计算汽穴体积，将每个控制体等分为两个相同的控制体。当为纯单相流动时，第i个控制体内两部分的压力和流量均由FVM计算。一旦任一小控制体压力达到或低于液体汽化压力，则两部分压力设为汽化压力水头，计算体内流量则依然采用Godunov格式求得。汽穴体积由汽穴上下游侧的流量差值计算：

其中，分别为n+1时步时，第i个控制体内，上、下游小控制体内的流量。

更进一步地，步骤5中，一旦汽穴体积小于或等于0，则汽穴溃灭，流动恢复为纯单相流动,压力和流量再次均由二阶FVM-Godunov格式计算。

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

(1)本发明提供的FVM-DVCM方法成功地克服了FVM方法在追踪汽穴并预测其发育、溃灭方面的难题，并且同MOC-DVCM方法一样简单且易于实现；(2)当采用相对细密的网格时，结果中不会出现人为的脉冲峰值；(3)非线性对流项很容易加入到解中，而这些项在MOC-DVCM中通常是忽略的，因此它可将该模型的应用范围扩展到包含马赫数的问题；(4)该计算方法通过时间算子分裂方法，为模拟多维水锤瞬变提供了框架，而MOC解较难扩展到包含多维瞬变的流体问题。

附图说明

图1为本发明的基本流程图；

图2为实施例的水库-管道-阀门-水库系统示意图；

图3为实施例的水锤区的网格离散系统图；

图4为实施例的水柱分离区的网格离散系统图；

图5为实施例下，网格数为32时，阀门末端的压力曲线图；

图6为实施例下，网格数为256时，阀门末端的压力曲线图；

图2中：1-上游水库；2-输水管；3-球阀；4-下游水库。

具体实施方式

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

基于Godunov格式的有压管道中水柱分离的模拟方法，按以下步骤进行：建立瞬变流基本微分方程，根据工程实例确定计算域、初始条件以及边界条件；根据FVM划分计算网格，并建立离散方程；采用时空均为二阶精度的Godunov格式求解纯对流时，单元控制体的数值通量；接着，基于二阶Runge-Kutta离散格式的时间算子分裂法，在纯对流控制方程的解中引入源项，从而得到最终解的二阶显式FVM-Godunov格式；最后，根据算得的压力判定是否发生水柱分离：若已形成，则进行水柱分离模块计算；反之，进入下一时步计算。图1给出了本发明的基本流程图。

下面将结合附图和实施例对本发明技术方案做进一步的详细描述。

实施例:Bergant和Simpson于1995年设计搭建了如图2所示的实验系统来研究水柱分离-弥合水锤现象。整个系统由上游压力水箱，铜管，下游球阀，下游压力水箱组成。铜管总长37.2m，内径22.1mm，管道倾角3.2°。水锤波速为1319m/s。伴随水柱分离的水锤事件由突然关闭下游球阀引起。实验工况参数见表1。

表1 Bergant和Simpson实验工况参数

本发明的具体步骤为：

步骤1：建立瞬变流基本微分方程，根据工程实例确定计算域、初始条件以及边界条件。

瞬变流基本微分方程为：

其中，沿管线距离x与时间t是自变量；H(x,t)是测压管水头；V(x,t)是平均截面速率；g是重力加速度；a是波速；f为达西-威斯巴哈摩阻系数；D为管径。

本实施例下，计算域为从上游水库出水口至阀门之间的管道；初始条件为球阀全开时，初始流速为V₀，初始压力根据上游水库静压水头减去相应的稳态摩阻损失得到；边界条件为：管道入口处，水库提供恒压边界，且等于上游水库静压水头；下游球阀快速关闭建立水锤瞬变。

步骤2：根据FVM划分计算网格，并建立离散方程。

图3为本实施例下水锤区网格离散系统。对第i个控制体，定义其上、下游界面编号分别为i-1/2、i+1/2。

建立离散方程的具体步骤为：

(a)微分方程(1)(2)的黎曼问题可近似为含常系数的线性双曲系统的黎曼问题：

其中是V的平均值，为一常数。

(b)在控制单元i(从界面i-1/2到i+1/2)及时间段Δt(从t到t+Δt)内对方程(3)积分，可以得到流动变量u的离散方程：

其中，上标n和n+1分别代表t和t+Δt时步。为u在整个控制体的平均值；f为单元界面处的通量；为源项；

步骤3：采用Godunov方法求解纯对流时的离散方程，得到纯对流时控制单元界面处的数值通量，并取得二阶精度。

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

步骤3.1：计算纯对流时，内部控制单元界面处通量

首先，根据Godunov方法，纯对流时，其黎曼问题为以下初值问题：

其中，为在n时步时，u分别到界面i+1/2左、右侧两侧的平均值。

对矩阵进行简化计算，可以得到其特征值及特征向量分别为：

因特征向量线性无关，进一步可得到：

求解四个未知系数α₁，α₂，β₁，β₂，可得：

黎曼问题(方程(5)和(6))一般解的原始变量形式为：

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

现在，利用方程(12)可得到变量在界面i+1/2处的精确解：

从而，对t∈[tⁿ,tⁿ⁺¹]，任一内部控制单元i(1<i<N)，在界面i+1/2处的通量为：

接着，引入MUSCL-Hancock格式来计算内部单元通量f_i+1/2，取得空间和时间上的二阶精度。具体步骤为：

(a)数据重构：采用每个控制单元[x_i-1/2,x_i+1/2]内的分段线性函数代替数据单元平均值则在极值点处的值为：

其中，Δ_i为选定的适度斜坡向量，用以增加计算格式的精度，并保证解中不出现虚假振荡。本发明中，Δ_i选择MINMOD限制器：

其中，

(b)推算：对每一个单元[x_i-1/2,x_i+1/2]，方程(15)中的边界外推值根据下式乘以0.5Δt推算：

(c)黎曼问题：为计算内部界面通量f_i+1/2，需结合下面数据求解传统的黎曼问题：

将方程(18)代入方程(14)，即可得到管中充满水时，通量在时间t＝[tⁿ,tⁿ⁺¹]内，所有内部单元界面i+1/2处的二阶格式。

步骤3.2：边界条件的实现

为取得二阶精度，分别在起始控制单元I₁上游侧、终点控制单元I_N下游侧构建两个虚拟控制单元I_-1、I₀，以及I_N+1、I_N+2，并假定在虚拟单元处的流动信息与边界处是一致的。从而可求解边界黎曼问题，且相应的Godunov通量f_1/2。和f_N+1/2也可像内部单元那样进行计算。

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

(a)上游水库恒水位边界：

在上游边界处，与负特征线相关的黎曼不变量为：

其中，H_1/2＝H_res，H_res为上游水库静压水头。

从而可推得变量u_1/2(t)＝(H_1/2,V_1/2)。根据虚拟单元处的假定，临近管道进口的虚拟单元I_-1和I₀的相应值为：

(b)下游阀门边界：

在下游边界处，与正特征线相关的黎曼不变量为：

下游边界条件为：阀门在阀门关闭时间T_c内关闭。水头、流速之间的关系为：

V_N+1/2＝0,T_c<t (23)

其中，H_D为出口压力水头，ζ_v为阀门处损失系数，变量u_N+1/2(t)＝(H_N+1/2,V_N+1/2)。根据虚拟单元处的假定，临近I_N的虚拟单元I_N+1和I_N+2处的相应值为：

此外，步骤3中，对流项满足CFL条件(Courant-Friedrichs-Lewy criterion)，进一步可推得CFL条件下的最大时间步长Δt_max,CFL：

其中，C_r为柯朗数。

步骤4：通过基于二阶Runge-Kutta离散格式的时间算子分裂法，在纯对流控制方程的解中引入源项，从而得到最终解的二阶显式FVM-Godunov格式。

具体实现过程为：

(a)纯对流时：

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

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

此外，引入的源项满足以下稳定性约束：

将方程(27)、(28)分别代入方程(29)，可以得到适用于源项的最大时间步长Δt_max,s：

结合方程(25)和(30)，可推得包含对流项和源项的最大允许时间步长为：

Δt_max＝min(Δt_max,CFL,Δt_max,s) (31)

步骤5：根据计算得到的压力判定是否形成水柱分离：若已形成，则进行水柱分离模块计算；反之，进入下一时步计算。

图4为本实施例下水柱分离区的网格离散系统图。假定汽穴仅在控制体的中部形成，为计算汽穴体积，将每个控制体等分为两个相同的控制体。当为纯单相流动时，第i个控制体内两部分的压力和流量均由FVM计算。一旦任一小控制体压力达到或低于液体汽化压力，则两部分压力设为汽化压力水头，计算体内流量则依然采用Godunov格式求得。汽穴体积由汽穴上下游侧的流量差值计算：

其中，Q_ui、Q_i分别为第i个控制体内，上、下游小控制体内的流量。

一旦汽穴体积小于或等于0，则汽穴溃灭，流动恢复为纯单相流动,压力和流量再次均由二阶FVM-Godunov格式计算。

通过以上方法编程计算后，将FVM-DVCM和MOC-DVCM计算结果均与实验数据做对比。计算中采用的网格数为32、256。图5-6分别给出了本实施例下阀门处的压力曲线图。可以看出，FVM-DVCM计算的压力的幅值和时间响应均与实验吻合良好，即使在计算单元数量少至32的情形下。当网格数量增加时，计算的压力波动强度会有一定程度的增大，但其响应规律保持不变。而在经典MOC-DVCM计算结果中，当网格数为32时，压力波动曲线与实验结果吻合良好。一旦网格数增加，就会出现不真实的压力脉冲和高频振荡。针对水柱分离现象，相比经典的MOC-DVCM，本发明提出的FVM-DVCM在细密网格下，可以获得更准确、稳定的结果。

Claims (8)

1.一种基于Godunov格式的有压管道中水柱分离的模拟方法，其特征在于，采用Godunov格式来模拟输水管道系统中水柱分离-弥合水锤现象，具体步骤如下：

步骤1：建立瞬变流基本微分方程，根据工程实例确定计算域、初始条件以及边界条件；

步骤2：根据FVM划分计算网格，并建立离散方程；

步骤3：采用Godunov方法求解纯对流时的离散方程，得到纯对流时控制单元界面处的数值通量，并取得二阶精度；

步骤4：通过基于二阶Runge-Kutta离散格式的时间算子分裂法，在纯对流控制方程的解中引入源项，从而得到最终解的二阶显式FVM-Godunov格式；

步骤5：根据计算得到的压力判定是否形成水柱分离：若已形成，则进行水柱分离模块计算；反之，进入下一时步计算；

步骤2中，对控制单元i，建立的流动变量u的积分方程为：

其中，上标n和n+1分别代表t和t+Δt时步；为u在整个控制体的平均值；H是测压管水头，V是平均截面速率；Δx为控制单元尺寸；dx为积分变量；f为单元界面处的通量；为源项；f为达西-威斯巴哈摩阻系数；D为管径；

步骤3包含以下子步骤：

步骤3.1：求解内部控制单元界面处通量；

首先，基于黎曼问题，根据Godunov格式，对任一内部控制单元i，1<i<N，界面i+1/2处的通量为：

其中，H是测压管水头；V是平均截面速率；是V的平均值，为一常数；a为波速；g是重力加速度；为在n时步时，u分别到界面i+1/2左、右侧两侧的平均值；

接着，通过引入MUSCL-Hancock格式计算内部单元通量f_i+1/2，从而取得空间和时间上的二阶精度；

步骤3.2：求解边界控制单元界面处通量；为在边界面处也取得二阶精度，分别在起始控制单元1上游侧、终点控制单元N下游侧构建两个虚拟控制单元I_-1、I₀，以及I_N+1、I_N+2，并假定在虚拟单元处的流动信息与边界处是一致的；从而可求解边界黎曼问题，且相应的Godunov通量f_1/2和f_N+1/2也可像内部单元那样进行计算。

2.如权利要求1所述的基于Godunov格式的有压管道中水柱分离的模拟方法，其特征在于，步骤3.1中，在引入MUSCL-Hancock格式计算内部单元通量f_i+1/2过程中，需选择斜率限制器，以保证解中不出现虚假振荡。

3.如权利要求1所述的基于Godunov格式的有压管道中水柱分离的模拟方法，其特征在于，对流项满足CFL条件，进一步可推得CFL条件下的最大时间步长Δ_tmax,CFL：

其中，C_r为柯朗数，为矩阵的特征值。

4.如权利要求1所述的基于Godunov格式的有压管道中水柱分离的模拟方法，其特征在于，步骤4中，引入源项后，瞬变流基本微分方程解的二阶FVM-Godunov格式为：

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

5.如权利要求1所述的基于Godunov格式的有压管道中水柱分离的模拟方法，其特征在于，源项满足以下稳定性约束，并可推得适用于源项的最大时间步长Δ_tmax,s：

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

6.如权利要求1所述的基于Godunov格式的有压管道中水柱分离的模拟方法，其特征在于，包含对流项和源项的最大允许时间步长为：

Δt_max＝min(Δt_max,CFL,Δt_max,s) (6)

其中，Δt_max,CFL为CFL条件下所允许的最大时间步长；Δt_max,s为采用Runge-Kutta格式引入源项时，源项满足约束性条件时所允许的最大时间步长。

7.如权利要求1所述的基于Godunov格式的有压管道中水柱分离的模拟方法，其特征在于，步骤5中，水柱分离计算模块中，假定汽穴仅在控制体的中部形成，为计算汽穴体积，将每个控制体等分为两个相同的控制体；当为纯单相流动时，第i个控制体内两部分的压力和流量均由FVM计算；一旦任一小控制体压力达到或低于液体汽化压力，则两部分压力设为汽化压力水头，计算体内流量则依然采用Godunov格式求得；汽穴体积由汽穴上下游侧的流量差值计算：

其中，分别为第i个控制体在n、n+1时步的汽穴体积；Q_ui、Q_i分别为第i个控制体内，上、下游小控制体内的流量。

8.如权利要求7所述的基于Godunov格式的有压管道中水柱分离的模拟方法，其特征在于，步骤5中，一旦汽穴体积小于或等于0，则汽穴溃灭，流动恢复为纯单相流动,压力和流量再次均由二阶FVM-Godunov格式计算。