CN104050341B - 基于液相与气相耦合的调压室水位波动模拟方法及系统 - Google Patents

Abstract

本发明公开了一种基于液相与气相耦合的调压室水位波动模拟方法及系统，包括步骤：步骤1，采用Pressimann空间四点隐格式分别离散液体管道非恒定流控制方程和气体管道非恒定流控制方程；步骤2，将调压室和通气廊道管道化，构建垂直管道液气耦合模型；步骤3，采用追赶法的前扫描过程求解垂直管道液气耦合模型，获取水气分界面高程；步骤4，根据获得的水气分界面高程调整气体管道网格数和液体管道网格数，并保证水气分界面位于气体管道末网格线和液体管道末网格线间。本发明方法计算简便，模拟准确度可达到三维数值模拟法准确度，可广泛运用于水电站引水发电系统的工程设计，特别是带通气廊道的调压室的工程设计。

Description

基于液相与气相耦合的调压室水位波动模拟方法及系统

技术领域

本发明属于水电站数值模拟技术领域，特别涉及一种基于液相与气相耦合的调压室水位波动模拟方法及系统。

背景技术

随着水电建设发展，调压室应用越来越普遍，规模越来越庞大。一方面，水电站流道的水力设计、结构设计等与调压室位置、形式和尺寸密切相关；另一方面，一些关键的调节保证参数，比如蜗壳压力、尾水管压力、调压室最高涌浪水位、调压室最低涌浪水位等均受调压室参数影响，调压室对调节系统稳定性和调节品质等都具有非常大甚至决定性的影响。

工程设计时，一方面要使设计的调压室的最小面积满足压力、涌浪极值以及系统稳定的要求，并留有一定的安全系数；另一方面设计的调压室的最大面积受到地形地质、围岩稳定等诸多因素限制，同时，调压室投资占流道系统比例相当大，应此在满足调节保证和机组运行稳定的前提下，应尽量减少工程量，满足经济合理原则。

调压室的工作原理是利用扩大的断面和自由反射水击波，减小水击压力并改善机组的运行条件。设计阶段，需要采用数值模拟技术对整个引水发电系统的动态特性及瞬变过程进行仿真模拟，以确保各部分设计均满足要求。调压室作为引水发电系统的重要组成部分，其水位波动模拟的准确性直接影响引水发电系统仿真模拟的可靠性，因此数值模拟中调压室边界条件的选择极其重要。

调压室的型式多种多样，比如开敞式调压室、气垫式调压室等。在地下式水电站中，为减少调压室的高度和工程量，采用通气廊道将调压室与外界大气相连，因此调压室水位波动模拟同时涉及液相与气相。另外有一种能减小地下式水电站调压室工程量的方法，采用一面积较小的竖直连接管连接压力管道和调压室底部，由于竖直连接管面积小流速大，在数值模拟过程中则不能忽略竖直连接管内的水流惯性。

已有的调压室水位波动模拟方法主要有解析解法和数值解法。解析解法由于仅考虑简单调压室系统的水位波动，只能在诸多简化条件下单纯的计算水位波动，而不能与引水发电系统中机组等边界条件联合求解，在调压室设计中已经很少被采用。

调压室水位波动的数值解法有一维数值模拟法和三维数值模拟法。现有的一维数值模拟法中，把调压室当作蓄水池考虑，调压室内液体满足静压分布，不考虑水在调压室内的流动，忽略调压室边壁与水流间的摩擦阻力；然后，结合管道非恒定流基本方程和引水发电系统中边界条件模拟调压室内水位波动，一维数值模拟法的模拟过程中常常忽略气体管道对水位波动的影响。

三维数值模拟法更准确，采用三维数值模拟法不仅能获得一些宏的物理参数，比如调压室内水位波动随时间的变化规律、阻抗式调压室阻抗系数等，同时还能模拟诸如流场、漩涡等一维数值模拟法无法获得的微观现象。三维计数值模拟法能考虑调压室体型等复杂流道边界条件，并采用考虑了液相与气相间耦合作用的两相流模型。

现有的调压室水位波动的数值解法中，一维数值模拟法是假设调压室内水流满足静压分布，因此未考虑调压室内水流惯性，对于带有长连接管的调压室，连接管内水流惯性对水位波动和其他调保参数的影响非常明显，如不考虑调压室内水流惯性则会带来较大的计算误差，影响模拟准确性。另一方面，对于与通气廊道相连的调压室，一维数值模拟法忽略了通气廊道内气体流动对水位波动产生的影响，同时，已有的一维数值模拟法也不能准确模拟通气廊道内风速等重要参数。

三维数值模拟法虽然准确，但其计算一个工况所需时间很长，而所要计算的工况常常很多，因此在工程设计阶段完全采用三维数值模拟法并不现实；再者，由于水电站引水发电系统过渡过程的模拟包含水力-机械-电气三方面耦合，三维数值模拟法的优势在于流体方面，在机械控制、电气调节等方面还没有很好的模拟方法。因此对整个引水发电系统模拟时，三维数值模拟法存在一定局限性。

目前国内外对液相与气相耦合作用下的调压室水位波动的研究不多，传统模拟方法不能满足带通气廊道调压室的水位波动模拟需要。

发明内容

针对现有技术在模拟带通气廊道的调压室内水位波动时存在的不足，本发明利用液体管道非恒定流控制方程和气体管道非恒定流控制方程，提出了一种基于液相与气相耦合的调压室水位波动模拟方法及系统。

为解决上述技术问题，本发明采用如下的技术方案：

基于液相与气相耦合的调压室水位波动模拟方法，包括步骤：

步骤1，采用Pressimann空间四点隐格式分别离散液体管道非恒定流控制方程和气体管道非恒定流控制方程，获得对应的离散方程；

步骤2，将调压室和通气廊道管道化，以调压室内的水气交界面为界划分液体管道和气体管道，并以水气交界面为气体管道末断面和液体管道末断面，耦合液体管道与气体管道非恒定流，构建气体管道末断面质量流量和液体管道末断面体积流量的关系、气体管道末断面压力和液体管道末断面水头的关系以及液体管道末断面体积流量和水气分界面高程的关系，即获得了垂直管道液气耦合模型；

步骤3，基于步骤1获得的离散方程，采用追赶法前扫描过程获得液体管道末断面体积流量和水头的关系以及气体管道末断面质量流量和压力的关系，所述的液体管道的前扫描方向为调压室底部向水气交界面方向，所述的气体管道的前扫描方向为通气廊道出口向水气交界面方向；

步骤4，联合垂直管道液气耦合模型、液体管道末断面体积流量和水头的关系以及气体管道末断面质量流量和压力的关系，获得液体管道末断面的体积流量和水头、气体管道末断面的质量流量和压力、以及水气分界面高程；

步骤5，根据液体管道末断面的体积流量和水头以及气体管道末断面的质量流量和压力，采用追赶法后扫描过程获得液体管道各断面的体积流量和水头以及气体管道各断面的质量流量和压力；所述的液体管道的后扫描方向为水气交界面向调压室底部方向，所述的气体管道的后扫描方向为水气交界面向通气廊道出口方向；

步骤6，根据获得的水气分界面高程调整气体管道网格数和液体管道网格数，并保证水气分界面位于气体管道末网格线和液体管道末网格线间。

步骤1进一步包括子步骤：

1.1采用Pressimann空间四点隐格式差分格式离散液体管道非恒定流控制方程，获得液体管道相邻断面的水头增量和体积流量增量间的关系，即液体管道非恒定流控制方程的离散方程；所述的水头增量指计算时刻液体管道断面的水头较上一时刻水头的增量；所述的体积流量增量指计算时刻液体管道断面的体积流量较上一时刻体积流量的增量；

1.2采用Pressimann空间四点隐格式差分格式离散气体管道非恒定流控制方程，获得气体管道相邻断面的压力增量和质量流量增量间的关系，即气体管道非恒定流控制方程的离散方程；所述的压力增量指计算时刻气体管道断面的压力较上一时刻压力的增量；所述的质量流量增量指计算时刻气体管道断面的质量流量较上一时刻质量流量的增量。

步骤2中所述的气体管道末断面质量流量和液体管道末断面体积流量的关系为：其中，M_U为气体管道末断面的质量流量，Q_N为液体管道末断面的体积流量，ρ_气为空气密度。

步骤2中所述的气体管道末断面压力和液体管道末断面水头的关系为：

其中，H_N为液体管道末断面的水头，Z为调压室水气分界面高程，p_U为气体管道末面的压力，ρ_气为空气密度，g为重力加速度。

步骤2中所述的液体管道末断面体积流量和水气分界面高程的关系为：

其中，Z为调压室水气分界面高程，Q_N为液体管道末断面的体积流量；Z⁰为前一时刻调压室水气分界面高程，为前一时刻液体管道末断面的体积流量；Δt为时间步长，F为调压室横截面积。

步骤3中所述的液体管道末断面体积流量和水头的关系的获得进一步包括子步骤：

3a-1基于调压室与上游管道交界处特征线法的C⁺方程、调压室与下游管道交界处特征线法的C^-方程、以及调压室与上游管道、下游管道交界处的体积流量连续性方程和能量守恒方程，获得液体管道首断面的体积流量增量和水头增量的关系；所述的液体管道首断面指液体管道与上游管道、下游管道相连的断面；所述的体积流量增量为液体管道首断面的体积流量较上一时刻体积流量的增量，所述的水头增量为液体管道首断面的水头较上一时刻水头的增量；

3a-2根据液体管道首断面的体积流量增量和水头增量的关系以及液体管道非恒定流控制方程的离散方程，获得液体管道各断面的体积流量增量和水头增量的递推关系；

3a-3根据液体管道各断面的体积流量增量和水头增量的递推关系，采用追赶法的前扫描过程递推出液体管道末断面的体积流量增量和水头增量的关系；

且，步骤3中所述的气体管道末断面质量流量和压力的关系的获得具体为：

3b-1根据气体管道非恒定流控制方程的离散方程以及气体管道首断面的压力增量，获得气体管道各断面的质量流量增量和压力增量的关系；所述的气体管道首断面即气体管道与大气相连的断面；所述的压力增量为气体管道断面的压力较上一时刻压力的增量；所述的质量流量增量为气体管道断面的质量流量较上一时刻质量流量的增量；

3b-2根据气体管道各断面的质量流量增量和压力增量的关系，采用追赶法的前扫描过程递推出气体管道末断面的质量流量增量和压力增量的关系。

步骤5进一步包括子步骤：

5.1根据液体管道末断面的体积流量获得液体管道末断面的体积流量增量，根据追赶法前扫描过程中获得的液体管道各断面体积流量增量和水头增量的关系式，采用追赶法后扫描过程获得液体管道各断面的体积流量和水头；

5.2根据气体管道末断面的质量流量获得气体管道末断面的质量流量增量，根据追赶法前扫描过程中获得的气体管道各断面质量流量增量和压力增量的关系式，采用追赶法后扫描过程获得气体管道各断面的质量流量和压力。

步骤6进一步包括子步骤：

6.1若当前水气分界面位于气体管道末网格线和液体管道末网格线间，则保持当前气体管道网格数和当前液体管道网格数不变；

6.2若当前水气分界面升高且超过了气体管道末网格线，则减少气体管道网格数，增加液体管道网格数，使得当前水气分界面位于气体管道末网格线和液体管道末网格线间；

6.3若当前水气分界面降低且低于液体管道末网络线，则减少液体管道网格数，增加气体管道网格数，使得当前水气分界面位于气体管道末网格线和液体管道末网格线间。

本发明还提供了一种基于液相与气相耦合的调压室水位波动模拟系统，包括：

离散模块，用来采用Pressimann空间四点隐格式分别离散液体管道非恒定流控制方程和气体管道非恒定流控制方程，获得对应的离散方程；

垂直管道液气耦合模型构建模块，用来将调压室和通气廊道管道化，以调压室内的水气交界面为界划分液体管道和气体管道，并以水气交界面为气体管道末断面和液体管道末断面，耦合液体管道与气体管道非恒定流，构建气体管道末断面质量流量和液体管道末断面体积流量的关系、气体管道末断面压力和液体管道末断面水头的关系以及液体管道末断面体积流量和水气分界面高程的关系，即获得了垂直管道液气耦合模型；

前扫描模块，用来基于步骤1获得的离散方程，采用追赶法的前扫描过程获得液体管道末断面体积流量和水头的关系以及气体管道末断面质量流量和压力的关系，所述的液体管道的前扫描方向为调压室底部向水气交界面方向，所述的气体管道的前扫描方向为通气廊道出口向水气交界面方向；

垂直管道液气耦合模型求解模块，用来联合垂直管道液气耦合模型、液体管道末断面体积流量和水头的关系以及气体管道末断面质量流量和压力的关系，获得液体管道末断面的体积流量和水头、气体管道末断面的质量流量和压力、以及水气分界面高程；

后扫描模块，根据液体管道末断面的体积流量和水头以及气体管道末断面的质量流量和压力，采用追赶法的后扫描过程获得液体管道各断面的体积流量和水头以及气体管道各断面的质量流量和压力；所述的液体管道的后扫描方向为水气交界面向调压室底部方向，所述的气体管道的后扫描方向为水气交界面向通气廊道出口方向；

动网格调整模块，用来根据获得的水气分界面高程调整气体管道网格数和液体管道网格数，并保证水气分界面位于气体管道末网格线和液体管道末网格线间。

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

1、考虑了调压室内的水流惯性和边壁摩擦阻力影响，考虑了气相与液相的耦合作用对调压室内水位波动的影响，在获得调压室内水位波动变化过程的同时，还可以获得通气廊道内风速、压力等参数的变化过程。

2、本发明为一维模拟方法，将复杂矩阵方程求解简化为代数方程的求解，编程简单，高效省时，模拟准确度可达到三维数值模拟法的准确度，满足工程设计需要。

3、适用范围广，可模拟各种类型调压室水位波动。

4、液体管道和气体管道的空间步长一致，但时空间步长与时间步长相互独立，在空间步长取较大情况下，仍可以取较小空间步长，有利于采用动网格追踪水气分界面。

附图说明

图1为Preissmann空间四点隐格式差分格式示意图；

图2为调压室与通气廊道相连示意图；

图3为水气交界面及垂直管道液气耦合模型的示意图；

图4为本发明方法流程示意图。

图中，1-上游管道；2-调压室下游管道；3-水流方向；4-调压室；5-液体管道网格线；6-水气交界面；7-气体管道网格线；8-气流方向；9-通气廊道；10-通气廊道出口。

具体实施方式

本发明基于液体管道与气体管道的非恒定流耦合思想，将调压室当成管道，采用管道非恒定流控制方程求解调压室内各网格的流量和水头，采用相同的差分格式求解气体管道瞬变过程，并将液体与气体的控制方程耦合在同一模型中，该模型不仅考虑了调压室内水流惯性影响，而且考虑了液体与气体间的相互作用；同时，提出了动网格技术追踪水位波动，实现液体网格与气体网格的自动增减。

下面将结合附图和具体实施方式对本发明技术方案做进一步说明。

本发明方法的具体步骤如下：

步骤1，采用Pressimann四点空间隐格式分别离散液体管道非恒定流控制方程和气体管道非恒定流控制方程。

(1)液体管道非恒定流控制方程和气体管道非恒定流控制方程

描述一维非棱柱体的有压非恒定流控制方程包括连续性方程(1)和动量方程(2)：

$V\frac{\∂H}{\∂x}+\frac{\∂H}{\∂t}+\frac{a^2}{g}\frac{\∂V}{\∂x}+\frac{a^{2}V}{\mathrm{gA}}\frac{\∂A}{\∂x}-\sin \mathrm{\β}\·V=0---\left(1\right)$

$g\frac{\∂H}{\∂x}+V\frac{\∂V}{\∂x}+\frac{\∂V}{\∂t}+f\frac{V\left|V\right|}{2D}=0---\left(2\right)$

式(1)～(2)中：

x表示距离，单位：m；

t表示时间，单位：s；

V为管道断面平均流速，单位：m/s；

H为测压管水头，单位：m；

a为水击波波速，单位：m/s；

A为管道过流断面面积，单位：m²；

β为管轴线与水平面的夹角，单位：rad；

g为重力加速度，单位：m/s²；

f为Darcy-Weisbach摩阻系数，无单位；

D为管道内直径，单位：m。

气体是流体的一种，气体管道的非恒流控制方程与液体管道类似。一维气体管道非恒定流控制方程仍包括连续性方程(3)和动量方程(4)：

$\frac{B^2}{A}\frac{\∂M}{\∂x}+\frac{\∂p}{\∂t}=0---\left(3\right)$

$\frac{\∂p}{\∂x}+\frac{1}{A}\frac{\∂M}{\∂t}+\frac{\mathrm{pg}}{B^2}\sin \mathrm{\β}+\frac{f{B}^2{M}^2}{2D{A}^{2}p}=0---\left(4\right)$

式(3)～(4)中：

x表示距离，单位：m；

t表示时间，单位：s；

M为管道断面质量流量，单位：kg/s；

p为管道断面压力，单位：Pa；

B为声波波速，单位：m/s；

A为管道过流断面面积，单位：m²；

β为管轴线与水平面的夹角，单位：rad；

g为重力加速度，单位：m/s²；

f为Darcy-Weisbach摩阻系数，无单位；

D为管道内直径，单位：m。

(2)液体管道非恒定流控制方程和气体管道非恒定流控制方程的离散

图1和式(5)为考虑权重系数θ的Preissmann空间四点隐格式差分格式，图1中，i+1、i表示相邻网格的编号，n+1、n表示相邻时刻的编号，Δx为空间步长，Δt为时间步长。采用该差分格式离散公式(1)～(2)所示的液体管道非恒定流控制方程，并以增量方式表示水头和体积流量，获得式(6)和式(7)所示的离散方程：

$\begin{array}{c}f(x,t)=\frac{\mathrm{\θ}}{2}(f_{i+1}^{n+1}+f_i^{n+1})+\frac{1-\mathrm{\θ}}{2}(f_{i+1}^n+f_i^n)\\ \frac{\∂f}{\∂x}=\mathrm{\θ}\frac{f_{i+1}^{n+1}-f_i^{n+1}}{\mathrm{\Δx}}+(1-\mathrm{\θ})\frac{f_{i+1}^n-f_i^n}{\mathrm{\Δx}}\\ \frac{\∂f}{\∂t}=\frac{f_{i+1}^{n+1}-f_{i+1}^n+f_i^{n+1}-f_i^n}{2\mathrm{\Δt}}\end{array}---\left(5\right)$

式(5)中，x表示距离，t表示时间；Δx为空间步长，Δt为时间步长；i+1、i表示相邻网格的编号，即相邻断面编号；n+1、n表示相邻时刻的编号；权重系数θ取值范围为0.5～1.0，根据经验进行取值，权重系数θ取值越大，Preissmann空间四点隐格式越稳定，但准确性降低。

对于液体管道，用Q/A来表示管道断面平均流速V，Q表示管道断面体积流量，以增量方式表示水头和体积流量，获得式(6)和式(7)的离散方程：

A1·ΔH_i+1+B1·ΔQ_i+1＝C1·ΔH_i+D1·ΔQ_i+F1         (6)

A2·ΔH_i+1+B2·ΔQ_i+1＝C2·ΔH_i+D2·ΔQ_i+F2         (7)

式(6)～(7)中，系数A1、A2、B1、B2、C1、C2、D1、D2、F1、F2均为已知值，ΔH表示计算时刻管道断面水头较上一时刻的增量，ΔQ表示计算时刻管道断面体积流量较上一时刻的增量，即，ΔH_i+1和ΔH_i分别表示液体管道断面i+1、i水头较上一时刻的增量，ΔQ_i+1表示液体管道断面i+1、i体积流量较上一时刻的增量。

对于气体管道非恒定流控制方程，以增量方式表示压力和质量流量，可得到式(8)和式(9)所示的离散方程：

A1'·Δp_i+1+B1'·ΔM_i+1＝C1'·Δp_i+D1'·ΔM_i+F1'         (8)

A2'·Δp_i+1+B2'·ΔM_i+1＝C2'·Δp_i+D2'·ΔM_i+F2'         (9)

式(8)～(9)中，系数A1'、A2'、B1'、B2'·、C1'、C2'、D1'、D2'、F1'、F2'均为已知值，ΔM表示计算时刻气体管道断面质量流量较上一时刻的增量，Δp表示计算时刻气体管道断面压力较上一时刻的增量，即，ΔM_i+1表示气体管道断面i+1的质量流量较上一时刻的增量，ΔM_i表示气体管道断面i的质量流量较上一时刻的增量，Δp_i+1表示气体管道断面i+1的压力较上一时刻的增量，Δp_i表示气体管道断面i的压力较上一时刻的增量。

公式(6)～(7)中各系数值如下：

$A1=1+\mathrm{\θ}\frac{\mathrm{\Δt}}{\mathrm{\Δx}}(\frac{Q_{i+1}}{A_{i+1}}+\frac{Q_i}{A_i});$

$B1=\frac{a^2}{g{A}_{i+1}}\frac{2\mathrm{\θ\Δt}}{\mathrm{\Δx}};$

$C1=-1+\mathrm{\θ}\frac{\mathrm{\Δt}}{\mathrm{\Δx}}(\frac{Q_{i+1}}{A_{i+1}}+\frac{Q_i}{A_i});$

$D1=\frac{a^2}{g{A}_i}\frac{2\mathrm{\θ\Δt}}{\mathrm{\Δx}};$

$A2=g\frac{\mathrm{\θ\Δt}}{\mathrm{\Δx}};$

$B2=\frac{1}{2A_{i+1}}[1+\frac{\mathrm{\θ\Δt}}{\mathrm{\Δx}}(\frac{Q_{i+1}}{A_{i+1}}+\frac{Q_i}{A_i})];$

$C2=g\frac{\mathrm{\θ\Δt}}{\mathrm{\Δx}};$

$D2=\frac{1}{2A_i}[1-\frac{\mathrm{\θ\Δt}}{\mathrm{\Δx}}(\frac{Q_{i+1}}{A_{i+1}}+\frac{Q_i}{A_i})];$

$\begin{array}{c}F1=-\frac{\mathrm{\Δt}}{\mathrm{\Δx}}(\frac{Q_{i+1}}{A_{i+1}}+\frac{Q_i}{A_i})(H_{i+1}-H_i)-\frac{a^2}{g}\frac{2\mathrm{\Δt}}{\mathrm{\Δx}}(\frac{Q_{i+1}}{A_{i+1}}-\frac{Q_i}{A_i})-\frac{a^2}{g}\frac{\mathrm{\Δt}}{\mathrm{\Δx}}(\frac{Q_{i+1}}{{A_{i+1}}^2}+\frac{Q_i}{{A_i}^2})(A_{i+1}-A_i)\\ +\mathrm{\Δ}t\sin \mathrm{\β}(\frac{Q_{i+1}}{A_{i+1}}+\frac{Q_i}{A_i})\end{array}$

$F2=-\frac{\mathrm{\Δt}}{2\mathrm{\Δt}}(\frac{{Q_{i+1}}^2}{{A_{i+1}}^2}-\frac{{Q_i}^2}{{A_i}^2})-g\frac{\mathrm{\Δt}}{\mathrm{\Δx}}(H_{i+1}-H_i)-\frac{\mathrm{\Δt}}{4}(\frac{f{Q}_i\left|Q_i\right|}{A_i^2{D}_i}+\frac{f{Q}_{i+1}\left|Q_{i+1}\right|}{A_{i+1}^2{D}_{i+1}})$

其中，Q_i+1和Q_i分别表示计算时刻液体管道断面i+1、i的体积流量，A_i+1和A_i分别表示计算时刻液体管道断面i+1、i的过流断面面积。

公式(8)～(9)中各系数值如下：

${A1}^{\′}=\frac{1}{2\mathrm{\Δt}};$

${B1}^{\′}=\frac{B^2\mathrm{\θ}}{\mathrm{A\Δx}};$

${C1}^{\′}=-\frac{1}{2\mathrm{\Δt}};$

${D1}^{\′}=\frac{B^2\mathrm{\θ}}{\mathrm{A\Δx}};$

${F1}^{\′}=-\frac{B^2}{\mathrm{A\Δx}}(M_{i+1}-M_i);$

${A2}^{\′}=\frac{\mathrm{\θ}}{\mathrm{\Δx}}+\frac{\mathrm{g\θ}\sin \mathrm{\β}}{2B^2}-\frac{f{B}^2(M_i-M_{i+1})|M_i+M_{i+1}|}{2D{A}^2{(p_i+p_{i+1})}^2};$

${B2}^{\′}=\frac{1}{2\mathrm{A\Δt}}+\frac{f{B}^2\left(\right|M_i+M_{i+1}\left|\right)}{2{\mathrm{DA}}^2(p_i+p_{i+1})};$

${C2}^{\′}=\frac{\mathrm{\θ}}{\mathrm{\Δx}}-\frac{\mathrm{g\θ}\sin \mathrm{\β}}{2B^2}+\frac{f{B}^2(M_i+M_{i+1})|M_i+M_{i+1}|}{2D{A}^2{(p_i+p_{i+1})}^2};$

${D2}^{\′}=-\frac{1}{2\mathrm{A\Δt}}-\frac{f{B}^2\left(\right|M_i+M_{i+1}\left|\right)}{2{\mathrm{DA}}^2(p_{i+1}+p_i)};$

${F2}^{\′}=-\frac{p_{i+1}-p_i}{\mathrm{\Δx}}-\frac{g\sin \mathrm{\β}}{2B^2}(p_{i+1}+p_i)-\frac{{\mathrm{fB}}^2(M_i+M_{i+1})|M_i+M_{i+1}|}{4D{A}^2(p_i+p_{i+1})};$

其中，M_i+1和M_i分别表示计算时刻气体管道断面i+1、i的质量流量，p_i+1和p_i分别表示计算时刻气体管道断面i+1、i的压力。

步骤2，将调压室和通气廊道管道化，耦合液体管道非恒定流与气体管道非恒定流，获得垂直管道液气耦合模型，见公式(10)～(12)。

图2所示为某地下电站调压室(4)，其通过通气廊道(9)与大气相连，由于调压室面积较大，通气廊道面积相对较小，调压室内水位波动将会引起通气廊道内产生很大风速，耦合液体管道非恒定流与气体管道非恒定流，获得模拟调压室在气体作用下的水位波动。

将调压室和通气廊道管道化，以调压室内水气交界面为界划分液体管道和气体管道，并按预设空间步长分别对液体管道和气体管道进行网格划分。液体管道的首网格线对应液体管道与管道(1,2)相连的断面，液体管道的末网格线对应调压室的水气交界面；气体管道的首网格线对应气体管道与大气相连的断面，气体管道的末网格线同样对应调压室的水气交界面。见图2，网格线1～N表示液体管道各断面，网格线1为液体管道的首网格线，网格线1～U表示气体管道各断面，调压室水位即水气交界面(6)高程。水气交界面及本发明提出的垂直管道液气耦合模型见图3，气体管道末断面采用气体管道断面U表示，液体管道末断面采用液体管道断面N表示，Z为水气分界面高程。

假设计算时间步长Δt内，流入或流出液体管道末断面的液体体积等于流出或流入气体管道末断面的气体体积，即存在式(10)：

液体管道中水流流向是从调压室底部到气液分界面，气体管道中气体流向是从通气廊道出口到气液分界面，由于水流流向和气体流向相反，所以式(10)中存在负号。

同时，在水气交界面处还满足压力方程式(11)和水位流量关系方程式(12)：

$Z=Z^0+\frac{(Q_N^0+Q_N)\mathrm{\Δt}}{2F}---\left(12\right)$

压力方程式(11)反映了液体管道末断面水头和气体管道末断面压力间的关系，水位流量关系方程式(12)反映了调压室内液体管道末断面体积流量和水气分界面高程的关系。

式(10)～(12)中：M_U为计算时刻气体管道末断面的质量流量，Q_N为计算时刻液体管道末断面的体积流量，p_U为计算时刻气体管道末断面的压力，H_N为计算时刻液体管道末断面的水头，Z为计算时刻调压室水气分界面高程，Z⁰为前一时刻的调压室水气分界面高程，为前一时刻液体管道末断面上的体积流量，F为调压室横截面积，ρ_气为空气密度，Δt为时间步长，g为重力加速度。

在气体管道末断面及液体管道末断面间区域，一共含有M_U、p_U、Q_N，H_N和Z五个未知数，上述垂直管道液气耦合模型仅包括三个方程，因此本发明增加液体管道末断面体积流量和水头的关系以及气体管道末断面质量流量和压力的关系使调压室内垂直管道液气耦合模型封闭，液体管道末断面体积流量和水头的关系以及气体管道末断面质量流量和压力的关系通过步骤3中的前扫描过程获得。

步骤3，采用追赶法法求解垂直管道液气耦合模型。

(1)采用追赶法的前扫描过程获得液体管道末断面体积流量和水头的关系

在调压室与管道(1，2)交界处，存在上游管道末断面特征线法的C⁺方程(13)、下游管道首断面的C^-方程(14)、以及体积流量连续性方程(15)和能量守恒方程(16)：

Q_sn＝Q_CP-C_QP·H_sn         (13)

Q_x1＝Q_CM+C_QM·H_x1         (14)

Q_sn＝Q₁+Q_x1(15)

$H_{\mathrm{sn}}=H_{x1}=H_1+\mathrm{\α}{Q}_1^2---\left(16\right)$

式(13)～(16)中，Q和H分别为管道断面的体积流量和水头，下标表示管道断面编号；Q_CP和C_QM为特征线法系数，α为调压室阻抗系数。

联合方程(13)～(16)求解液体管道首断面体积流量和水头的关系，并以增量形式表示为：

ΔQ₁＝EE₁·ΔH₁+FF₁         (17)

式(17)中，ΔQ₁为液体管道首断面体积流量较上一时刻的增量，ΔH₁为液体管道首断面水头较上一时刻的增量，系数EE₁和FF₁均为已知值。

式(17)即为由显格式特征线方程、连续性方程和能量守恒方程获得的隐格式管道首断面的边界条件，将其带入式(6)和(7)，获得隐格式液体管道内部断面的递推关系：

ΔQ_i＝EE_i·ΔH_i+FF_i         (18)

ΔH_i＝L_i·ΔH_i+1+M_i·ΔQ_i+1+N_i         (19)

式(18)～(19)中各系数只与前时刻各断面函数值H_i、Q_i有关，为已知。断面i上的EE_i、FF_i仅依赖于前一断面i-1上的EE_i-1、FF_i-1，i＝1断面上的EE₁、FF₁已知，通过联合求解公式(13)～(16)获得。已知EE₁、FF₁依次可推算所有的EE_i、FF_i。另外，L_i、M_i、N_i取决于EE_i和FF_i。在液体管道的前扫描过程中，首先，确定EE₁、FF₁；然后，依次推算液体管道各断面的EE_i、FF_i、L_i、M_i、N_i；最后，获得液体管道末断面的EE_N和FF_N。

以从调压室底部到水气交界面方向为前扫描方向，采用追赶法的前扫描过程得到隐格式液体管道末断面水头和体积流量的关系：

ΔQ_N＝EE_N·ΔH_N+FF_N         (20)

式(20)即为求解垂直管道液气耦合模型所需的液体管道末断面体积流量和水头的关系式，其中，ΔQ_N为液体管道末断面体积流量较上一时刻的增量，ΔH_N为液体管道末断面水头较上一时刻的增量。

(2)采用追赶法的前扫描过程获得气体管道末断面流量和压力的关系

气体管道的前扫描方向为从通气廊道出口向水气交界面方向。

由于气体管道首断面直接与大气相连，其出口压力保持不变，因此气体管道首断面压力增量始终为0，即：

Δp₁＝0         (21)

在气体管道断面1和断面2间，存在离散化后的连续性方程和动量方程的增量表达式：

A1·Δp₁+B1·ΔM₁＝C1·Δp₂+D1·ΔM₂+F1         (22)

A2·Δp₁+B2·ΔM₁＝C2·Δp₂+D2·ΔM₂+F2         (23)

公式(22)～(23)中，Δp₁为气体管道断面1的压力较上一时刻压力的增量，Δp₂为气体管道断面2的压力较上一时刻的增量；ΔM₁为气体管道断面1的质量流量较上一时刻质量流量的增量，ΔM₂为气体管道断面2的质量流量较上一时刻质量流量的增量。

联合方程(21)～(23)得到气体管道断面2的质量流量和压力的关系：

ΔM₂＝GG₂·Δp₂+KK₂         (24)

式(24)中，系数GG₂和KK₂均为已知值。

将式(24)代入气体管道非恒定流控制方程的离散方程(8)～(9)，得到气体管道各断面质量流量和压力的关系：

ΔM_i＝GG_i·Δp_i+KK_i         (25)

Δp_i＝U_i·Δp_i+1+V_i·ΔM_i+1+W_i         (26)

式(25)～(26)中，系数GG_i、KK_i、U_i、V_i、W_i均为已知值，为气体管道的前扫描过程中可以确定的中间变量。

气体管道末断面的质量流量和压力的关系式如下：

ΔM_U＝GG_U·Δp_U+KK_U         (27)

式(27)即为求解垂直管道液气耦合模型所需的气体管道末断面质量流量和压力的关系式。

(3)求解调压室内水位波动模型

联合垂直管道液气耦合模型(10)～(12)、液体管道末断面体积流量和水头的关系式(20)以及气体管道末断面质量流量和压力的关系式(27)，求解液体管道末断面的体积流量Q_N和水头H_N、气体管道末断面的质量流量M_U和压力p_U、以及水气分界面高程Z。

采用追赶法的后扫描过程分别获得液体管道和气体管道各断面的参数：

根据液体管道末断面体积流量Q_N获得对应的体积流量增量ΔQ_N，将体积流量增量ΔQ_N代入式(19)求得液体管道断面N-1的水头增量ΔH_N-1，将水头增量ΔH_N-1代入公式(18)求得液体管道断面N-1的体积流量增量ΔQ_N-1，以此类推，利用式(18)和(19)可以求解液体管道各断面上的体积流量增量和水头增量，从而获得液体管道各断面的体积流量和水头。

同理，对于气体管道，根据气体管道末断面的质量流量M_U获得对应的质量流量增量ΔM_U，将质量流量增量ΔM_U代入式(26)，求得气体管道断面(U-1)的压力增量Δp_U-1，将压力增量Δp_U-1代入公式(25)求解气体管道断面(U-1)的质量流量增量ΔM_U-1。利用式(25)和式(26)间的递推关系求解气体管道断面2～U上的质量流量和压力，再利用公式(21)～(22)求解气体管道首断面的质量流量和压力。

步骤4，采用动网格技术跟踪调压室液面波动。

由于调压室内水气交界面高程Z随时间变化，因此需要随时调整液体管道网格数和气体管道网格数，保持空间步长Δx不变。

网格数的调整分以下三种情况进行：

情况1：水气交界面高程Z在液体管道末网格线N和气体管道末网格线U间波动，但未越过网格线N和U。

这种情况下，液体管道网格数和气体管道网格数与上一时刻相同，无需要增减网格即直接进入下一时刻的计算。

当水气交界面高程Z满足不等式(28)时，说明水气交界面高程Z未越过网格线N和U：

$\left\{\begin{array}{c}(N-2)\&CenterDot;\mathrm{\&Delta;x}\&le;Z^0<(N-1)\&CenterDot;\mathrm{\&Delta;x}\\ (N-2)\&CenterDot;\mathrm{\&Delta;x}\&le;Z(N-1)\&CenterDot;\mathrm{\&Delta;x}\end{array}\right.---\left(28\right)$

式(28)中，Z⁰为上一时刻水气交界面高程，Z为当前水气交界面高程；Δx为空间步长，即相邻网格线间距离。

情况2：水气交界面高程Z升高，并越过气体管道末网格线U。

这种情况下，减少气体管道网格数，增加液体管道网格数。本具体实施中，调压室液体管道网格和气体管道网格的空间步长均为0.5m，时间步长为0.01s，因此在一个时间步长内，水气交界面高程Z的变化最多跨越一条网格线。

当水气交界面高程Z满足不等式(29)时，说明水气交界面高程升高，且越过气体管道网格线U：

$\left\{\begin{array}{c}(N-2)\&CenterDot;\mathrm{\&Delta;x}\&le;Z^0<(N-1)\&CenterDot;\mathrm{\&Delta;x}\\ (N-1)\&CenterDot;\mathrm{\&Delta;x}\&le;Z<N\&CenterDot;\mathrm{\&Delta;x}\end{array}\right.---\left(29\right)$

式(29)中，Z⁰为上一时刻水气交界面高程，Z为当前水气交界面高程；Δx为空间步长，即相邻网格线间距离。

此时，增加液体管道网格数，减少气体管道的网格数减小：

N＝N+1         (30)

U＝U-1         (31)

液体管道新增网络线断面上流量和压力水头初始值为：

$Q_{N+1}^0=Q_N---\left(32\right)$

$H_{N+1}^0=Z---\left(33\right)$

情况3：水气交界面高程Z降低，并越过液体管道末网格线N。

这种情况下，增加气体管道网格数，减少液体管道网格数。本具体实施中，调压室液体管道网格和气体管道网格的空间步长均为0.5m，时间步长为0.01s，因此在一个时间步长内，水气交界面高程Z的变化最多跨越一条网格线。

当水气交界面高程Z满足不等式(34)时，说明水气交界面高程Z降低，且越过液体管道网格线N：

$\left\{\begin{array}{c}(N-2)\&CenterDot;\mathrm{\&Delta;x}\&le;Z^0<(N-1)\&CenterDot;\mathrm{\&Delta;x}\\ (N-3)\&CenterDot;\mathrm{\&Delta;x}\&le;Z(N-2)\&CenterDot;\mathrm{\&Delta;x}\end{array}\right.---\left(34\right)$

式(34)中，Z⁰为上一时刻水气交界面高程，Z为当前水气交界面高程；Δx为空间步长，即相邻网格线间距离。

此时，减少液体管道网格数，增加气体管道网格数：

N＝N-1         (35)

U＝U+1         (36)

气体管道新增网络线断面上流量和压力初始值为：

$p_{U+1}^0=p_U---\left(38\right)$

其中，p_U为气体管道末断面压力。

通过以上方法调整相邻时刻间网格数后，再采用步骤3继续跟踪下一时刻调压室水气交界面高程(即调压室水位)。

Claims (9)

1.基于液相与气相耦合的调压室水位波动模拟方法，其特征在于，包括步骤：

步骤1，采用Pressimann空间四点隐格式分别离散液体管道非恒定流控制方程和气体管道非恒定流控制方程，获得对应的离散方程；

步骤2，将调压室和通气廊道管道化，以调压室内的水气交界面为界划分液体管道和气体管道，并以水气交界面为气体管道末断面和液体管道末断面，耦合液体管道与气体管道非恒定流，构建气体管道末断面质量流量和液体管道末断面体积流量的关系、气体管道末断面压力和液体管道末断面水头的关系以及液体管道末断面体积流量和水气分界面高程的关系，即获得了垂直管道液气耦合模型；

步骤3，基于步骤1获得的离散方程，采用追赶法前扫描过程获得液体管道末断面体积流量和水头的关系以及气体管道末断面质量流量和压力的关系，所述的液体管道的前扫描方向为调压室底部向水气交界面方向，所述的气体管道的前扫描方向为通气廊道出口向水气交界面方向；

步骤4，联合垂直管道液气耦合模型、液体管道末断面体积流量和水头的关系以及气体管道末断面质量流量和压力的关系，获得液体管道末断面的体积流量和水头、气体管道末断面的质量流量和压力、以及水气分界面高程；

步骤5，根据液体管道末断面的体积流量和水头以及气体管道末断面的质量流量和压力，采用追赶法后扫描过程获得液体管道各断面的体积流量和水头以及气体管道各断面的质量流量和压力；所述的液体管道的后扫描方向为水气交界面向调压室底部方向，所述的气体管道的后扫描方向为水气交界面向通气廊道出口方向；

步骤6，根据获得的水气分界面高程调整气体管道网格数和液体管道网格数，并保证水气分界面位于气体管道末网格线和液体管道末网格线间。

2.如权利要求1所述的基于液相与气相耦合的调压室水位波动模拟方法，其特征在于：

步骤1进一步包括子步骤：

1.1采用Pressimann空间四点隐格式差分格式离散液体管道非恒定流控制方程，获得液体管道相邻断面的水头增量和体积流量增量间的关系，即液体管道非恒定流控制方程的离散方程；所述的水头增量指计算时刻液体管道断面的水头较上一时刻水头的增量；所述的体积流量增量指计算时刻液体管道断面的体积流量较上一时刻体积流量的增量；

1.2采用Pressimann空间四点隐格式差分格式离散气体管道非恒定流控制方程，获得气体管道相邻断面的压力增量和质量流量增量间的关系，即气体管道非恒定流控制方程的离散方程；所述的压力增量指计算时刻气体管道断面的压力较上一时刻压力的增量；所述的质量流量增量指计算时刻气体管道断面的质量流量较上一时刻质量流量的增量。

3.如权利要求1所述的基于液相与气相耦合的调压室水位波动模拟方法，其特征在于：

所述的气体管道末断面质量流量和液体管道末断面体积流量的关系为：其中，M_U为气体管道末断面的质量流量，Q_N为液体管道末断面的体积流量，ρ_气为空气密度。

4.如权利要求1所述的基于液相与气相耦合的调压室水位波动模拟方法，其特征在于：

所述的气体管道末断面压力和液体管道末断面水头的关系为：其中，H_N为液体管道末断面的水头，Z为调压室水气分界面高程，p_U为气体管道末面的压力，ρ_气为空气密度，g为重力加速度。

5.如权利要求1所述的基于液相与气相耦合的调压室水位波动模拟方法，其特征在于：

所述的液体管道末断面体积流量和水气分界面高程的关系为：

其中，Z为调压室水气分界面高程，Q_N为液体管道末断面的体积流量；Z⁰为前一时刻调压室水气分界面高程，为前一时刻液体管道末断面的体积流量；Δt为时间步长，F为调压室横截面积。

6.如权利要求1所述的基于液相与气相耦合的调压室水位波动模拟方法，其特征在于：

步骤3中所述的液体管道末断面体积流量和水头的关系的获得进一步包括子步骤：

3a-1基于调压室与上游管道交界处特征线法的C⁺方程、调压室与下游管道交界处特征线法的C-方程、以及调压室与上游管道、下游管道交界处的体积流量连续性方程和能量守恒方程，获得液体管道首断面的体积流量增量和水头增量的关系；所述的液体管道首断面指液体管道与上游管道、下游管道相连的断面；所述的体积流量增量为液体管道首断面的体积流量较上一时刻体积流量的增量，所述的水头增量为液体管道首断面的水头较上一时刻水头的增量；

3a-2根据液体管道首断面的体积流量增量和水头增量的关系以及液体管道非恒定流控制方程的离散方程，获得液体管道各断面的体积流量增量和水头增量的递推关系；

3a-3根据液体管道各断面的体积流量增量和水头增量的递推关系，采用追赶法的前扫描过程递推出液体管道末断面的体积流量增量和水头增量的关系；

且，步骤3中所述的气体管道末断面质量流量和压力的关系的获得具体为：

3b-1根据气体管道非恒定流控制方程的离散方程以及气体管道首断面的压力增量，获得气体管道各断面的质量流量增量和压力增量的关系；所述的气体管道首断面即气体管道与大气相连的断面；所述的压力增量为气体管道断面的压力较上一时刻压力的增量；所述的质量流量增量为气体管道断面的质量流量较上一时刻质量流量的增量；

3b-2根据气体管道各断面的质量流量增量和压力增量的关系，采用追赶法的前扫描过程递推出气体管道末断面的质量流量增量和压力增量的关系。

7.如权利要求1所述的基于液相与气相耦合的调压室水位波动模拟方法，其特征在于：

步骤5进一步包括子步骤：

5.1根据液体管道末断面的体积流量获得液体管道末断面的体积流量增量，根据追赶法前扫描过程中获得的液体管道各断面体积流量增量和水头增量的关系式，采用追赶法后扫描过程获得液体管道各断面的体积流量和水头；

5.2根据气体管道末断面的质量流量获得气体管道末断面的质量流量增量，根据追赶法前扫描过程中获得的气体管道各断面质量流量增量和压力增量的关系式，采用追赶法后扫描过程获得气体管道各断面的质量流量和压力。

8.如权利要求1所述的基于液相与气相耦合的调压室水位波动模拟方法，其特征在于：

步骤6进一步包括子步骤：

6.1若当前水气分界面位于气体管道末网格线和液体管道末网格线间，则保持当前气体管道网格数和当前液体管道网格数不变；

6.2若当前水气分界面升高且超过了气体管道末网格线，则减少气体管道网格数，增加液体管道网格数，使得当前水气分界面位于气体管道末网格线和液体管道末网格线间；

6.3若当前水气分界面降低且低于液体管道末网络线，则减少液体管道网格数，增加气体管道网格数，使得当前水气分界面位于气体管道末网格线和液体管道末网格线间。

9.基于液相与气相耦合的调压室水位波动模拟系统，其特征在于，包括：

离散模块，用来采用Pressimann空间四点隐格式分别离散液体管道非恒定流控制方程和气体管道非恒定流控制方程，获得对应的离散方程；

垂直管道液气耦合模型构建模块，用来将调压室和通气廊道管道化，以调压室内的水气交界面为界划分液体管道和气体管道，并以水气交界面为气体管道末断面和液体管道末断面，耦合液体管道与气体管道非恒定流，构建气体管道末断面质量流量和液体管道末断面体积流量的关系、气体管道末断面压力和液体管道末断面水头的关系以及液体管道末断面体积流量和水气分界面高程的关系，即获得了垂直管道液气耦合模型；

前扫描模块，用来基于离散模块获得的离散方程，采用追赶法的前扫描过程获得液体管道末断面体积流量和水头的关系以及气体管道末断面质量流量和压力的关系，所述的液体管道的前扫描方向为调压室底部向水气交界面方向，所述的气体管道的前扫描方向为通气廊道出口向水气交界面方向；

垂直管道液气耦合模型求解模块，用来联合垂直管道液气耦合模型、液体管道末断面体积流量和水头的关系以及气体管道末断面质量流量和压力的关系，获得液体管道末断面的体积流量和水头、气体管道末断面的质量流量和压力、以及水气分界面高程；

后扫描模块，根据液体管道末断面的体积流量和水头以及气体管道末断面的质量流量和压力，采用追赶法的后扫描过程获得液体管道各断面的体积流量和水头以及气体管道各断面的质量流量和压力；所述的液体管道的后扫描方向为水气交界面向调压室底部方向，所述的气体管道的后扫描方向为水气交界面向通气廊道出口方向；

动网格调整模块，用来根据获得的水气分界面高程调整气体管道网格数和液体管道网格数，并保证水气分界面位于气体管道末网格线和液体管道末网格线间。