一种基于气液交界面耦合的调压室通气洞风速模拟方法与流程

本发明属于水利水电工程技术领域，具体的说涉及一种基于气液交界面耦合的调压室通气洞风速模拟方法。

背景技术：

随着我国水电建设事业的深入，西南山区一批大型、特大型水电站正在修建。西南地区水力资源具有高落差，河道陡峭、山谷狭窄等普遍特点，基于此，在综合考虑地形、地质、降低工程投资、安全性等多方面因素后，水电站多采用地下式厂房，这些电站普遍的工程特点是机组装机容量大且采用地下式厂房，为了减小管道内的水击压力，通常要修建调压室。对于这类地下式水电站来说，设置的调压室通常是修建在山体中的开敞式调压室，为保证调压室水面为大气压，该类调压室必须通大气，实现通风、减压等功能。地下式水电站调压室的通风是通过修建通气洞来实现的，如图1所示。在过渡过程中，不同工况转换时，调压室上方气体压力脉动是剧烈的，低压可低于大气压力，当负荷变化很大时，甚至可以减少到汽化压力，在这种正负压力剧烈交变的过程中，会对调压室内壁混凝土表面进行剥蚀，造成钢筋外露，强度降低，对水电站及建筑物构成威胁，通气洞的设置可有效的消除这些问题。

在水电站过渡过程中，机组甩负荷、增负荷都将引起调压室水位上下波动，导致调压室内气体体积、压强、流速等物理量发生改变(即处于非恒定状态)，进而调压室内非恒定状态的气体作为扰动源引起通气洞内气体的非恒定流动，可能产生较大的风速。风速的大小关系到通气洞的安全使用和工程投资，因此需对风速进行精确模拟，据此认识清楚通气洞在过渡过程中的风速特征与变化机理，为通气洞的设计、优化等提供依据。

对于调压室通气洞风速的模拟研究，前人尝试过试验和CFD模拟，但这两种方法都极为耗时耗力，应用起来很不方便与经济。一维的数值模拟方法可以快速高效的模拟气体的瞬变流过程，但是对于通气洞来说，其一端与调压室相连，存在气(通气洞气体)、液(调压室水面)的交界面，如何对该交界面的气液作用进行模拟是风速模拟的关键环节。另外，通气洞接大气端边界条件的模拟、气体瞬变流方程的建立、模拟程序的构建等等问题，都是整个风速模拟的重要环节。

技术实现要素：

本发明所要解决的技术问题是为克服现有技术的不足，提供一种基于气液交界面耦合的调压室通气洞风速模拟方法，实现调压室通气洞在过渡过程中的风速分布与变化的精确模拟。该方法建立了气液耦合模型，基于该气液耦合模型，进而建立了包含引水隧洞、调压室和通气洞在内的液体-气体运动整体模型，可以同时计算得到调压室水位波动变化过程和通气洞内各断面风速变化过程，反映了通气洞与调压室间的相互作用机理及通气洞内风速的波动本质。

本发明解决其技术问题所采用的技术方案是：

一种基于气液交界面耦合的调压室通气洞风速模拟方法，包括以下步骤：

步骤1，建立通气洞内的气体运动基本方程，建立引水隧洞内的液体运动基本方程，建立调压室涌浪基本方程，推求气体密度与压力及压缩性系数的关系，建立上游水库端、压力管道机组侧、通气洞接大气端及串联管连接处的边界条件，建立引水隧洞、调压室及通气洞内的流量和压力的初始条件；

步骤2，基于调压室涌浪基本方程和通气洞内的气体运动基本方程，推导气液交界面耦合方程，建立调压室水位液面与通气洞气体分界面的气液耦合模型；利用所述气液耦合模型将通气洞内的气体运动基本方程、引水隧洞内的液体运动基本方程、调压室涌浪基本方程、气体密度与压力及压缩性系数的关系耦合在一起，建立包含引水隧洞、调压室和通气洞在内的液体-气体运动整体模型；

步骤3，依据建立的包含引水隧洞、调压室和通气洞在内的液体-气体运动整体模型及特征线解法，结合针对设有通气洞的水电站引水隧洞与调压室建立的边界条件和初始条件，并对边界条件中的压力管道机组侧流量边界取不同流量变化过程，并采用计算机编程建立模拟程序，该模拟程序可以同时计算调压室水位波动和通气洞过渡过程中的风速。

而且，对于所述气液耦合模型：在推导气液交界面耦合方程时，要用到调压室涌浪基本方程和通气洞内的气体运动基本方程，两组方程根据一定的假设建立联系，可以得到气液交界面耦合方程，然后与气体、液体压力管道的特征线解法联立，可以实现对气液耦合模型各计算参数的模拟计算。

而且，所述调压室内涌浪水位基本方程为：

1)动量方程

和

2)连续性方程

Q_P1-Q_P2＝Q_S，

式中：A为引水管道断面积，m²；L为引水管道长度,m；Z为调压室内的水位,m，以向上为正；k为阻力系数；Q_P1为1号管段末端面流量，m³/s；Q_P2为2号断面首断面流量，m³/s；Q_S为进出调压室的流量，m³/s，以进入调压室内为正。

而且，所述气体管道流量边界方程为：

1)连续方程：

和

2)动量方程：

相关参数的定义如下：

(a)M为质量流量，p为气体绝对压力，θ为管道断面形心连线与水平面的夹角，α为惯性因子，B为波速，g为重力加速度，f为Darcy-Weisbach摩阻系数，M＝ρvA；

(b)波速可由状态方程计算得到：其中：ρ为气体密度、z为压缩性系数、R为气体常数、T为绝对温度；

(c)当通气洞断面中心连线沿+x方向升高时，θ取正值。

而且，所述假设具体为：

假定在每个计算时长Δt内，调压室内水体体积的变化量等于G断面气体体积的变化量。

而且，所述气液交界面耦合方程为：

和

式中：

R_S、K为中间变量，水头损失系数K＝Δt/(2F)，L为引水管道长度，F为调压室断面面积，R为管道水力半径，C为摩擦因子，ζ为局部回头损失系数，g为重力加速度；ρ为液体密度，ρ_气为气体密度；C_P，B₁，C_M，B₂为特征线解法中间变量；C_PP为气体管道G断面的中间变量；C_BB为参数；P为t时刻交界断面处的压强；P_-Δt为t-Δt时刻交界断面处的压强；Z_S0为t-Δt时刻调压室水位；Q_S0为t-Δt时刻调压室内的流量；H_P为t时刻P点的压力水头；Q_S为t时刻进出调压室内的流量。

而且，所述模拟程序的计算过程如下：

第1步：将引水隧洞和通气洞均划分为若干个计算子管段，同时形成计算断面，并选取计算时间步长；

第2步：计算引水隧洞内液体流量和压力水头、通气洞管道内部气体流量和压力；

第3步：计算上游边界结点流量和压力；

第4步：计算下游边界结点流量和压力；

第5步：在某一时刻每一网格结点上的变量值以及上、下游边界变量值计算已知后，再进行下一时刻变量的计算，并以此类推，计算至所需时长为止。

本发明的有益效果是：此种风速模拟方法理论依据充分，实现方式简单。通过气液交界面的耦合，实现引水隧洞液体与通气洞气体的同时、联合模拟，能够真实反映管道液体及调压室水位波动对通气洞风速的影响，得到的风速模拟结果真实可信、精度高，并且计算稳定、速度快，可以方便用于不同体型的通气洞在不同运行工况下的风速的模拟，模拟结果可直接用于指导调压室与通气洞的设计，为通气洞的体型设计与优化提供了可靠的保障。

附图说明

图1设置通气洞的带调压室水电站引水发电系统示意图；

图2气液耦合模型；

图3特征网格与特征线；

图4阻抗式调压室示意图；

图5气液耦合模型各管道断面示意图；

图6串联管边界条件；

图7通气洞典型断面布置图；

图8调压室涌浪对比图；

图9 A断面风速对比图；

图10 B断面风速对比图；

图11 C断面风速对比图；

图12 D断面风速对比图；

图13本发明方法的流程框图。

具体实施方式

下面对本发明内容作进一步详细描述。

气液耦合模型如图2所示。本发明建立通气洞气体和引水隧洞液体的瞬变流基本方程与特征线解法，然后再建立气液耦合模型与其他边界条件，实现通气洞气体和引水隧洞液体过渡过程中的非恒定流运动的联合模拟。

1通气洞内的气体运动基本方程及特征线解法

(1)基本方程

对于地下式水电站调压室通气洞中空气流动速度的模拟属于气体管道瞬变流的范畴，所采用的计算方法与液体管道瞬变流类似，同样需要首先建立瞬变流基本方程(连续方程、动量方程)。在方程建立之前，针对气体管道瞬变流的特殊性做如下假设：

(a)流动是等温的；

(b)管壁的弹性可以忽略；

(c)流动是一维的；

(d)摩擦系数是壁面粗糙度及雷诺数的函数，在瞬变流计算中采用定常运动的摩擦系数；

(e)沿着管道的动能的变化可以忽略。

1)连续方程：

2)动量方程：

相关参数的定义如下：

(a)A为管道断面面积，x为距离，t为时间，D为管道断面直径，M为质量流量，p为气体绝对压力，θ为管道断面形心连线与水平面的夹角，α为惯性因子，B为波速，g为重力加速度，f为Darcy-Weisbach摩阻系数，M＝ρvA，ρ为气体密度，v为流速；

(b)波速可由状态方程计算得到：其中：z为压缩性系数、R为气体常数、T为绝对温度；

(c)当交通洞断面中心连线沿+x方向升高时，θ取正值。

(2)特征线解法

连续性方程可写成：

动量方程可写成：

以上I₁、I₂方程可以合并成：

式中λ是一个未知参数，λ任意取两个不同的实数会得到两个不同的方程，以上将I₁、I₂方程可以合并得到I方程的目的是为了寻求是否存在两个特殊的值，使得上述方程可以变换得到常微分方程。

由微分学知识可以知道，若x是t的某一未知函数，则M(x,t)可表示如下式：

若上式中则式(5)的第一个括号里的项可简化为

同理对于来说，若则式(5)的第二个括号里的项可简化为在同一个方程中的值应该是一样的，则应满足：

即：

此时，方程I变形如下式：

方程(9)成立的条件是：

方程(10)中并不是质点速度，而是为了使方程I,成立，而形成的x、t的一种特定关系，将这种关系定义为特征线解法，便得到特征线方程组：

将C⁺方程(11)中各项均乘以dx＝Bdt/α，采用定时间间隔法，并沿着特征线A→P进行积分(如图3所示其中A、P、B为网格点)，积分得到有限差分的形式如下：

方程(15)中，M_P、M_A分别为P、A点的质量流量，P_P、P_A分别为P、A点的压力。方程(15)的前两项为完全差分形式。方程(15)由积分给定的后两项，必须沿着特征线A→P进行积分，使得方程(15)符合定常状态的方程

的解，因为定常流动是不定常流动状态的特殊情况，其中Δx为空间步长，s＝(2gΔx sinθ)/B²，P₁、P₂为Δx两侧的压力。于是得到C⁺方程：

同理，C^-方程(13)中各项均乘以dx＝Bdt/α，采用定时间间隔法，沿着特征线B→P进行积分，可以得到C^-方程：

方程(16)(17)中，利用了摩阻的二阶估算，把P点的流动与A点的流动或与B点的流动加以平均。其中绝对值符号的引用，是为了使方程对于相反的流动也是同样成立的。

C⁺、C^-方程中s＝(2gΔx sinθ)/B²，对于水平管线管s＝0，(e^s-1)/s＝1，利用沿特征线积分(C⁺：A→P，C^-：B→P)的方法可得C⁺方程、C^-方程：

C⁺：p_P＝C_P-C_BM_P

C^-：p_P＝C_M+C_BM_P

C_P、C_B、C_M为中间变量。

进而可得P点的压力和质量流量：

式中：

2引水隧洞内的液体运动基本方程及特征线解法

(1)基本方程

引水隧洞在电站中起着输送水流的作用，从输水系统的水流来看，流速从v₀变成v₀+Δv，引用流量从Q₀变成Q₀+ΔQ。由于水流惯性的作用，引水隧洞中各点的压力随着水流惯性的变化而变化，有着明显的升高和降低，以弹性波的形式沿管线传播，也即大家所熟悉的水击现象；在无压隧洞或明渠中各点的水深随着水流的惯性变化而变化，以重力波的形式沿渠道传播；调压室水位也发生上下波动，这种往返于引水隧洞和调压室内的水流运动是一种质量波。所以输水系统中的水流从初始的恒定状态过渡到另一种恒定状态的过程中呈非恒定流的特征。此过程即为引水隧洞的瞬变过程。

引水隧洞液体瞬变流基本方程即为水锤基本方程(动量方程、连续性方程)，其方程的建立基于以下假设：

(a)管道内的液体为一维流动，且流速沿横截面呈均匀分布；

(b)管壁和流体是线弹性的，应力与应变成比例；

(c)管道中恒定流的阻力计算公式对于瞬变流同样适用。

1)动量方程：

2)连续性方程：

式中：V为管道中流速，以向下游为正；H为压力水头；x为距离，t为时间，g为重力加速度，以管道进口为原点，向下游为正；c为水锤波速；d为管道直径；为纵坡。

(2)特征线解法

水击基本方程可写成如下形式：

用未知因子λ对Y₁、Y₂线性组合得：

若H＝H(x,t),V＝V(x,t)是上式方程的解，根据微分法则，全导数可写成：

很明显若令

Y式可写为：

由公式(26)解得的两个特征值

将带入方程(27)，得到特征线方程组：

将上述方程沿特征线C⁺A→P和C^-B→P积分(如图3所示)，其中摩阻损失项采用二阶精度数值积分，并用流量代替断面流速，整理可得：

C⁺:H_Pi＝C_A-B_AQ_Pi

C^-:H_Pi＝C_B-B_BQ_Pi

式中：C_A＝H_Ai-1+B_AQ_Ai-1-RQ_Ai-1|Q_Ai-1|，C_B＝H_Bi-1-B_BQ_Bi-1+RQ_Bi-1|Q_Bi-1|，

A_A、A_B为特征线上A、B两点所在断面的断面积；

H_Ai-1、Q_Ai-1为A节点(即i-1节点)在t时刻的压力水头和流量；

H_Bi-1、Q_Bi-1为B节点(即i+1节点)在t时刻的压力水头和流量。

由以上公式p_P＝C_P-C_BM_P和p_P＝C_M+C_BM_P可以解得在t+Δt时刻i节点的压力水头H_pi、流量Q_pi：

为了满足库朗条件及保证计算结果收敛，计算时采用的时间步长Δt和管道划分的空间步长Δx须满足条件

3调压室涌浪基本方程及涌浪求解

(1)调压室涌浪基本方程

调压室水位波动变化过程的模拟要基于以下假设条件：

(a)液体属于不可压缩性流体，引水隧洞和压力管道管壁是刚性的。即液体在引水隧洞和压力管道中是刚性运动的；

(b)调压室内的水体惯性比引水隧洞和压力管道内的水体惯性要小，可以忽略不计；

(c)沿程水力损失忽略不计；

(d)计算水头损失时，稳态条件下的相应流速公式同样适用于过渡过程中。

基于以上假设，可列出如图4所示的调压室底部动量方程和连续性方程如式(35)、(36)：

1)动量方程

2)连续性方程

Q_P1-Q_P2＝Q_S (36)

式中：Q为流量，t为时间，A为引水管道断面积，m²；L为引水管道长度,m；Z为调压室水位,m，以向上为正；k为阻力系数；Q_P1为1号管段末端面流量，m³/s；Q_P2为2号断面首断面流量，m³/s；Q_S为进出调压室的流量，m³/s，以进入调压室内为正。

(2)调压室水位与流量公式求解

如图4所示的阻抗式调压室，忽略了调压室内水体的惯性以及底部进口处的局部水头损失，1号管道末断面，2号管道首断面以及调压室底部断面，三个断面处的水头相等；即满足：

H₁＝H₂＝H_P (37)

由液体压力管道的特征线解法，可以得到1号管道末断面C⁺方程和2号管段首断面C^-方程，即：

H_P＝C_P-B₁Q_P1 (38)

H_P＝C_M+B₂Q_P2 (39)

方程(38)、(39)中，C_P，B₁，C_M，B₂为特征线解法中间变量。

连续性方程满足：

Q_P1-Q_P2＝Q_S (40)

Z与P点水压力满足：

H_P＝Z+R_SQ_S|Q_S| (41)

调压室水位与流量的关系为：

方程(42)中，F为调压室断面面积。

写成差分形式方程：

方程(43)中，Δt为时间步长，Z和Z_S0分别为t时刻和t-Δt时刻的调压室水位。

综合以上方程可得：

式中：水头损失系数K＝Δt/(2F)，F、L意义同前，R为管道水力半径，C为摩擦因子，R_S、K为中间变量，ζ为局部回头损失系数，g为重力加速度。

Q_P1为t-Δt时刻1号管段末断面的流量，m³/s；

Q_P2为t-Δt时刻2号管段首断面的流量，m³/s；

Z_S0为t-Δt时刻调压室水位，m；

Q_S0为t-Δt时刻调压室内的流量，m³/s；

H_P为t时刻P点的压力水头，m；

Q_S为t时刻进出调压室内的流量，m³/s。

(3)经典四阶龙格-库塔法

调压室涌浪基本方程的连续性方程和动量方程是两个一阶常微分方程，可以使用经典的四阶龙格-库塔法求解。将连续性方程和动量方程变换形式得到：

其中：t为时间，V为引水隧洞水流速度，Z为调压室水位，f为引水隧洞断面面积，g为重力加速度，L为引水隧洞长度，F为调压室断面面积，Q_t为压力管道流量，h_w为引水隧洞水头损失，h_r为调压室底部水头损失，f₁、f₂为函数符号。

根据经典的四阶龙格-库塔公式得到：

K₁＝Δtf₁(t,V,Z)

K₂＝Δtf₁(t+Δt/2,V+E₁/2,Z+K₁/2)

K₃＝Δtf₁(t+Δt/2,V+E₂/2,Z+K₂/2)

K₄＝Δtf₁(t+Δt,V+E₃,Z+K₃)

E₁＝Δtf₂(t,V,Z)

E₂＝Δtf₂(t+Δt/2,V+E₁/2,Z+K₁/2)

E₃＝Δtf₂(t+Δt/2,V+E₂/2,Z+K₂/2)

E₄＝Δtf₂(t+Δt,V+E₃,Z+K₃)

其中：t为时间，V为引水隧洞水流速度，Z为调压室水位，f₁、f₂为函数符号，K₁，K₂，K₃，K₄，E₁，E₂，E₃，E₄为中间变量，Z_t、V_t为t时刻的调压室水位与引水隧洞水流速度，Z_t+Δt、V_t+Δt为t+Δt时刻的调压室水位与引水隧洞水流速度，Δt为时间步长，根据组合工况的初始工况提供的调压室水位Z₀和引水隧洞水流速度V₀初始条件，计算初始工况下调压室水位。随后我们将初始工况计算所得调压室水位Z和引水隧洞水流速度V作为组合工况下调压室涌浪计算的初始条件。甩负荷工况下，将引水隧洞中的引用流量为零作为组合工况下调压室水位的极值发生的条件；增负荷工况下将引水隧洞中的引用流量等于压力管道中引用流量，作为组合工况下调压室水位的极值发生的条件。

4气体密度与压力、压缩性系数的关系推求

针对模拟的调压室交通洞的瞬变流过程，根据瞬变过程快速的特点，需考虑气体可压缩性的影响。故进行调压室交通洞过渡过程中的气体密度与压力、压缩性系数的关系推求。

气体状态方程：

p＝zρRT (46)

式(46)中，R为气体常数，T为绝对温度。虑气体可压缩性，将式(46)左边对压强p、右边对气体压缩系数z和密度ρ求全微分，得：

dp＝ρRTdz+zRTdρ (47)

结合气体压缩性系数的表达式，其中V为气体比容，将式(47)进行离散化处理，可得任一时刻t的压力、密度(分别用p、ρ表示)与其前一计算时刻t-Δt(分别用p_-Δt、ρ_-Δt表示)间的递推关系：

p-p_-Δt＝2ρpV-ρp_-ΔtV_-Δt-ρ_-ΔtpV (48)

整理式(48)可得：

求解式(49)，可得ρ的计算式：

其中：当断面压力大于前一时刻压力时，采用密度递推公式(50)计算；当断面压力小于前一时刻压力时，采用密度递推公式(51)计算。

5气液交界面耦合方程的推导

在推导气液交界面耦合方程时，要用到调压室水位基本方程和气体管道流量边界方程，两个方程根据一定的假设建立联系，可以得到气液交界面耦合方程，然后与气体、液体压力管道特征线解法联立，可以实现对气液耦合模型各计算参数的模拟计算。

调压室内涌浪水位基本方程具体为：方程(35)和(36)；

气体管道流量边界方程具体为：方程(1)和(2)；

一定的假设具体为：假定在每个计算时长Δt内，调压室内水体体积的变化量等于G断面气体体积的变化量；

气液交界面耦合方程具体为：方程(56)和(57)；

气体、液体压力管道特征线解法分别为通气洞内的气体运动基本方程特征线解法和引水隧洞内的液体运动基本方程特征线解法。

气液交界面耦合方程：如图5，假定在每个计算时长Δt内，调压室内水体体积的变化量等于G断面气体体积的变化量，即调压室内气体的变化量等于通气洞内的气体变化量，则调压室与通气洞连接断面G的气体流量M_J满足下式：

M_J＝-ρ_气Q_S (52)

由前文中边界条件可知通气洞与调压室连接处G断面的气体压强P满足：

P＝C_PP-C_BBM_J (53)

假定此时调压室水体上方气体压强与G断面气体压强相等，则液体压力管道R点的压力水头H_P与调压室水位Z以及水体上方的气体相对压强P-P_-Δt满足以下关系：

P＝C_PP+C_BBρ_气Q_S (55)

以上式中：Q_S为进出调压室的流量，m³/s；ρ为水体密度，kg/m³；C_PP为气体管道G断面的中间变量；C_BB为参数；P为t时刻交界断面处的压强，单位Pa；P_-Δt为t-Δt时刻交界断面处的压强，单位Pa。

将以上公式与四阶龙格-库塔法中公式联立，可以得到R点的压力水头H_P及进出调压室的流量Q_S公式：

公式(56)与公式(57)联立可以计算得到每一个Δt时刻的调压室水位Z以及G断面处的气体压强P，然后与气体管道特征线方程联立，进而求得每一时刻各断面的风速及压力变化。

6边界条件

(1)上游水库端

在过渡过程中，上游水库的水位可以认为保持不变，若忽略管道进口处的水头损失，则：

H_P＝H_RES (58)

上游水库端满足C^-特征线方程：

C^-:H_P＝C_M+B_PQ_P (59)

以上公式中：H_RES为上游水库的水位，以下游水面高程为基准，m；H_P为水库与管道连接点处的压力水头，m；Q_P为边界处管道进口流量，m³/s；C_M为中间变量。

(2)压力管道机组侧流量边界

本分发明气液耦合模型省略了水轮机组方程、下游压力管道方程以及下游水库方程，通过给定压力管道与机组进口端处某一线性流量变化过程Q_t作为边界条件，然后建立起C⁺边界方程：

C⁺：H_P＝C_P-B_PQ_P (60)

式中：Q_P＝Q_t；C_P为中间变量；B_P＝a/(gA)，a为波速，A为管道断面积。

(3)交通洞接大气端

直接与大气相接的边界，其压力恒为大气压力，即p_p＝p₀。应用C^-方程，可得此边界结点的未知量M_p：

(4)串联管

串联管连接处应满足质量流量连续条件，同时忽略连接处的局部损失，则还应满足压力相等条件。运用C⁺方程和C^-方程，可得到串联管边界条件：

7初始条件

对于定常等温流动，M为常数，此时将动量方程从x＝0、p＝p₁到x＝Δx、p＝p₂积分，整理可得：

式(64)中：s＝(2gΔx sinθ)/B²。此式给出定常状态的抛物线压强梯度，管道初始恒定状态各点的压强值都要满足此式，称为定常状态方程。对于水平管线，s＝0、e^s＝1、(e^s-1)/s＝1，此时式(64)可以变为：

对于初始状态管道内的气体处于静止的情况，各断面的流量均为零，管道内各点的压强处处等于外界大气的压强值。

本发明的工作原理与过程如下：

本发明通过建立调压室水位液面与通气洞气体分界面的气液耦合模型，将通气洞内的气体与引水隧洞内的液体耦合在一起，再通过特征线解法，实现两者在过渡过程中非恒定流动的联合模拟，精确地得到通气洞内风速的变化过程。具体工作过程如下：

依据以上建立的包含引水隧洞、调压室和通气洞在内的液体-气体运动整体模型及特征线求解方法，结合针对本发明的研究对象而提出的边界条件与初始条件，对气液耦合模型的压力管道流量边界取不同流量变化过程，并采用C++语言进行计算机编程，建立可以同时计算调压室水位波动和通气洞过渡过程中的风速模拟程序。模拟方法与模拟程序的计算过程如下：

(1)将引水隧洞和通气洞均划分为若干个计算子管段，子管段划分方法：各管段长度ΔL＝B*Δt，Δt为计算时间步长，长度为L的通气洞可以划分为n＝L/ΔL段，同时形成n+1个断面；

(2)引水隧洞流量和压力水头的计算：对管道内部任一计算断面，根据t＝0时刻的边界条件与初始条件，计算中间变量C_P、C_M、C_B的值，由式(33)、(34)计算t+Δt时刻的流量Q_P和压力水头H_P；通气洞管道内部气体流量和压力的计算：对管道内部任一计算断面，根据t＝0时刻的边界条件与初始条件，计算中间变量C_P1、C_M1、C_B1的值，由式(18)、(19)计算t+Δt时刻的质量流量M_P和压力P_p；

(3)上游边界结点流量和压力的计算：对于气液耦合模型，上游边界(上游水库端)，外无其它结点，边界断面处压力水头近乎恒定值，满足C^-方程；

(4)下游边界结点流量和压力的计算：对于液体压力管道，流量边界为给定流量变化过程，满足C⁺方程；同样对于气体管道，上游边界结点(通大气端)外无结点，因此只能采用C⁺方程，大气压力作为上游已知边界条件；

(5)在t+Δt时刻每一网格结点上的变量值以及上、下游边界变量值计算已知后，再进行t+2Δt时刻变量的计算，并以此类推，计算至所需时长为止。

下面结合实施例对本发明作进一步详细描述，但是本发明不限于给出的例子。

某如图2所示的设置通气洞的带调压室水电站，参数如下：L₁＝1000m，A₁＝93.986m²；L₂＝500m，A₂＝93.986m²；调压室断面积F＝500m²；气体压力管道LL₁＝680m，AA₁＝20m²；LL₁₀＝85m，AA₁₀＝20m²。其中，L₁、A₁分别为引水隧洞的长度、断面积；L₂、A₂分别为压力管道的长度、断面积，LL₁、AA₁分别为通气洞水平段的长度、面积，LL₁₀、AA₁₀分别为通气洞竖直段的长度、断面积。

以下内容是根据前面的方程及边界条件，通过模拟方法编制相应的计算程序，然后对调压室与通气洞连接断面处在耦合与不耦合两种情况下进行了对比分析(包括涌浪波动变化过程、典型断面风速变化过程)。

注：耦合情况是指利用特征线解法及耦合断面方程，可以同时计算得到每一时刻的调压室水位以及各断面风速变化过程；不耦合情况是指先利用液体特征线方程计算得到调压室水位波动变化过程，然后建立调压室水位变化速率与连接断面处的气体流量之间的联系，得到气体管道流量边界方程，然后再利用气体管道特征线方程对管道内各断面风速进行数值模拟。

利用上面的参数及相应程序，计算得到了耦合与不耦合两种情况下的调压室涌浪波动过程及各典型断面风速波动过程，见图8～图12：

由图8～图12可知：

(1)在体型参数及边界条件不变的情况下，耦合与不耦合情况计算得到的调压室涌浪值，其涌浪变化规律基本相同，涌浪差值(耦合条件下涌浪值减去不耦合情况下涌浪值)相差很小，差值的最大值为-0.082m，此时不耦合计算得到的涌浪值为8.39m,差值对涌浪的影响比例为0.977％，差值产生的原因是由于耦合情况下气液交界断面处的压强是不断变化的，其相对压强不等于零，而不耦合情况其相对压强恒等于零，由公式(4-3)可知，该断面处的气体压强会对调压室水位产生一定的影响；

(2)耦合与不耦合情况下，各断面风速值变化规律一致；调压室水位变化速度较大时，风速值由明显的弹性波和质量波叠加而成，而调压室水位变化速度较小时，载波即弹性波的影响很小，可以忽略，风速波动近乎为光滑的曲线，由此可以得知，调压室水位变化速率也会对载波产生影响，原因在于载波与气体惯性成正比，惯性越大，载波影响越显著，而惯性与质量成正比，当调压室水位波动变化减弱时，每个Δt时段内通气洞内运动的气体流量减小，惯性作用减弱，因而载波的影响也会减弱；

(3)从A、B、C、D各断面风速变化规律可以看出，从A-D断面(即通气洞接大气端至通气洞与调压室连接端)，载波对风速变化曲线影响越来越小，至D断面风速为光滑的曲线，前面已经介绍了风速波动曲线由弹性波和质量波叠加而成，质量波由调压室水位变化速率决定，而弹性波的影响是气体惯性决定的，运动的气体遇到出口处的大气被反射，在靠近大气端的断面处，气体惯性最大，因而载波对该断面处的风速影响也最大，靠近调压室端则与之相反；

(4)耦合与不耦合情况下，两者计算得到的风速值存在一定的偏差，最大差值为-3.488m/s，此时不耦合计算得到的风速值为-25.554m/s，差值对风速的影响比例为13.650％。