CN105468844A

CN105468844A - Analogy method of water-gas coupling transient flow in pipeline

Info

Publication number: CN105468844A
Application number: CN201510819625.5A
Authority: CN
Inventors: 周领; 王欢; 马佳杰; 潘天文; 刘德有; 王沛; 夏林
Original assignee: Hohai University HHU
Current assignee: Hohai University HHU
Priority date: 2015-11-23
Filing date: 2015-11-23
Publication date: 2016-04-06
Anticipated expiration: 2035-11-23
Also published as: CN105468844B

Abstract

The invention discloses a simulating method of water-air coupling transient flow in a pipeline based on a finite volume method, and considers the influence of the convection item impacting the water body and the depth of tailwater during the water-air coupling process. Firstly, the calculation area is divided into three parts: impact water body, water-air interface and trapped air mass. For the water body part, firstly, the calculation area of the water body is discretized into calculation grid units, and the discrete format of the basic differential equation of transient flow is established; then, the numerical flux of the control volume of the water body unit is solved using the Godunov scheme. Runge-Kutta method is used to solve the numerical integral term of the water body discrete equation. The dynamic tracking of the water-air interface is realized by using the Godunov scheme and the ideal gas state equation simultaneously. Simulations with second-order accuracy are achieved by introducing slope limiters to suppress spurious numerical oscillations. Accurate and rigorous numerical simulation and transient analysis can be carried out on the hydraulic phenomenon of water-air coupling effect of water flow impacting trapped air mass in the water pipeline system.

Description

Simulation Method of Water-Air Coupled Transient Flow in Pipeline

技术领域technical field

本发明涉及一种基于有限体积法的管道内水-气耦合瞬变流的模拟方法，属于水电站(泵站)水力学数值计算技术领域。The invention relates to a simulating method of water-gas coupling transient flow in a pipeline based on a finite volume method, and belongs to the technical field of hydraulic numerical calculation of a hydropower station (pumping station).

背景技术Background technique

在实际水电站、泵站等输水管道系统中，其沿线的凸起段、封闭端及放空检修段等部位常可能滞留气团，在系统启停及工况转换过程中，将可能会发生水-气耦合作用的复杂水力瞬变现象，其对应的异常水锤(或称：气-水锤)可能大到足以导致管道爆破的程度。实际输水系统所出现的破坏事故，很多与这种含滞留气团瞬变流有关。然而，至今为止，实际管道系统的设计标准只考虑管道载满水的情况，并不考虑滞留气团的存在及其危害；对于含滞留气团的情况，管道设计中尚无相应的计算标准，且已有成果尚不完善。因此，为保证输水管道系统的安全运行，针对滞留气团所引起的瞬变流现象进行深入细致的研究，具有重大的现实意义。In actual water pipeline systems such as hydropower stations and pumping stations, air masses may often remain in the raised sections, closed ends, and vented maintenance sections along the lines, and water- The complex hydraulic transient phenomenon of gas coupling, and its corresponding abnormal water hammer (or: gas-water hammer) may be large enough to cause pipeline explosion. Most of the destruction accidents in the actual water delivery system are related to the transient flow with trapped air mass. However, so far, the design standards of the actual pipeline system only consider the situation that the pipeline is full of water, and do not consider the existence and harm of trapped air mass; for the situation with trapped air mass, there is no corresponding calculation standard in the pipeline design, and it has been The results are not perfect yet. Therefore, in order to ensure the safe operation of the water pipeline system, it is of great practical significance to conduct in-depth and detailed research on the transient flow phenomenon caused by the trapped air mass.

很对学者对上述复杂瞬变流进行了研究，并提出了相应的模拟方法。现有的数学模型基本上为一维模型，主要包括基于特征性法的弹性水体模型和基于刚性理论的刚性水体模型。与弹性水体模型相比，刚性水体模型具有简单、快捷等优点，但是，其仅适用于较大体积的滞留气团情况的压力预测。现有的弹性水体模型主要是通过特性线法进行求得，能够有效的预测多数简单管路等瞬变过程的压力变化。当涉及管道波速变化的管道系统时，特性线法处理起来非常麻烦。此外，对于水流冲击滞留气团的水-气耦合作用瞬变流的模拟，已有的模型均忽略了动量方程和连续方程中对流项的影响，这是因为：(1)对于流动较为缓慢的水力瞬变问题，对流项影响很小；(2)考虑对流项将大大增加了特征线法的求解难度。水流快速冲击滞留气团的复杂水-气耦合作用瞬变流涉及较高的马赫数甚至剧烈的流速变化，其对流项势必对流动有一定的影响。Many scholars have studied the complex transient flow mentioned above, and put forward corresponding simulation methods. The existing mathematical models are basically one-dimensional models, mainly including the elastic water body model based on the characteristic method and the rigid water body model based on the rigid theory. Compared with the elastic water body model, the rigid water body model has the advantages of simplicity and quickness, but it is only suitable for the pressure prediction of the larger volume of trapped air mass. The existing elastic water body model is mainly obtained by the characteristic line method, which can effectively predict the pressure change of most simple pipelines and other transient processes. The characteristic line method is cumbersome when it comes to piping systems with varying pipe wave velocities. In addition, for the simulation of the transient flow of the water-air coupling effect of water impacting the trapped air mass, the existing models have ignored the influence of the convective term in the momentum equation and the continuity equation, because: (1) For the relatively slow-flowing hydraulic For transient problems, the convection term has little influence; (2) Considering the convection term will greatly increase the difficulty of solving the characteristic line method. The complex water-air coupling transient flow in which the water flow rapidly impacts the trapped air mass involves a high Mach number or even a severe flow velocity change, and its convective term is bound to have a certain impact on the flow.

有限体积法已广泛用于双曲线系统的求解，如：气体动力学方面。该方法能很好遵守动量和质量守恒，并且对于不连续问题提供有效的模拟。已有学者尝试将有限体积法用于简单管道内水锤问题的研究。然而，对于水流快速冲击滞留气团的复杂水-气耦合作用瞬变流，以及考虑尾水情况，尚无关于有限体积法应用的报道。The finite volume method has been widely used in the solution of hyperbolic systems, such as gas dynamics. The method obeys the conservation of momentum and mass well, and provides an efficient simulation for discontinuity problems. Scholars have tried to apply the finite volume method to the study of water hammer in simple pipelines. However, there is no report on the application of the finite volume method to the complex water-air coupling transient flow where the water flow rapidly impacts the trapped air mass, and considering the tail water.

因此，如果能够充分考虑水-气耦合作用过程中对流项、尾水等因素的影响，，引入新的方法和技巧消除缺陷，就能有助于水流冲击滞留气团过程中水-气耦合作用瞬变压力的数值模拟方法的不断改进。Therefore, if the influence of factors such as convection and tail water during the water-air coupling process can be fully considered, and new methods and techniques are introduced to eliminate defects, it will be helpful for the instantaneous water-air coupling action during the process of water impacting the trapped air mass. The continuous improvement of the numerical simulation method of variable pressure.

发明内容Contents of the invention

目的：为了弥补现有弹性水体模型的插值复杂等的不足，本发明提供一种管道内水-气耦合瞬变流的模拟方法，基于有限体积法，算法简单，易于实现，可以充分考虑对流项、尾水等因素，且避免繁琐插值等问题。Purpose: In order to make up for the shortcomings of the existing elastic water body model, such as complex interpolation, the present invention provides a simulation method for water-air coupling transient flow in pipelines, based on the finite volume method, the algorithm is simple, easy to implement, and the convection item can be fully considered , tail water and other factors, and avoid problems such as cumbersome interpolation.

技术方案：为解决上述技术问题，本发明采用的技术方案为：Technical solution: In order to solve the above-mentioned technical problems, the technical solution adopted in the present invention is:

一种基于有限体积法的管道内水-气耦合瞬变流的模拟方法，采用有限体积法和Godunov格式来模拟有压管道系统中含滞留气团的水-气耦合作用水力瞬变现象，具体步骤如下：A method for simulating water-gas coupled transient flow in pipelines based on finite volume method, using finite volume method and Godunov scheme to simulate the hydraulic transient phenomenon of water-gas coupling in pressurized pipeline system with trapped air mass, specific steps as follows:

步骤1：将管道系统内瞬变流进行划分为水体、水-气交界面、滞留气团三大部分，并建立相应的控制方程，根据工程实例确定初始条件以及边界条件；Step 1: Divide the transient flow in the pipeline system into three parts: water body, water-air interface, and trapped air mass, and establish corresponding control equations, and determine the initial conditions and boundary conditions according to engineering examples;

步骤2：根据有限体积法划分计算网格，并建立离散方程；Step 2: Divide the calculation grid according to the finite volume method, and establish the discrete equation;

步骤3：采用Godunov格式求解水体的离散方程的数值通量，并取得二阶精度；Step 3: Use the Godunov scheme to solve the numerical flux of the discrete equation of the water body, and obtain the second-order accuracy;

步骤4：通过基于Runge-Kutta方法对水体的离散方程的数值积分项进行求解，从而得到二阶显式有限体积法的Godunov格式；Step 4: Solve the numerical integral term of the discrete equation of the water body based on the Runge-Kutta method to obtain the Godunov scheme of the second-order explicit finite volume method;

步骤5：给出二阶显式有限体积法的Godunov格式所满足的CFL条件(Courant-Friedrichs-Lewycriterion)；Step 5: Give the CFL condition (Courant-Friedrichs-Lewycriterion) satisfied by the Godunov scheme of the second-order explicit finite volume method;

步骤6：采用Godunov格式和理想气体状态方程联立实现水-气交界面的动态追踪。Step 6: Use the Godunov scheme and the ideal gas state equation to realize the dynamic tracking of the water-gas interface.

进一步地，步骤2中，对第i个控制单元，定义其上、下游界面编号分别为i-1/2、i+1/2。Further, in step 2, for the i-th control unit, define its upstream and downstream interface numbers as i-1/2 and i+1/2, respectively.

进一步地，步骤2中，对控制单元i，建立的流动变量u的积分方程为：Further, in step 2, for the control unit i, the integral equation of the flow variable u established is:

${U u}_{i i}^{n no + + 11} = = {U u}_{i i}^{n no} + + \frac{Δ Δ t t}{Δ Δ x x} (({f f}_{i i + + 11 / / 22} - - {f f}_{i i - - 11 / / 22})) + + \frac{Δ Δ t t}{Δ Δ x x} {&Integral; &Integral;}_{i i - - 11 / / 22}^{i i + + 11 / / 22} s the s d d x x - - - - - - ((11))$

其中，上标n和n+1分别代表t和t+Δt时步；为u在整个控制体的平均值； $u = (\begin{matrix} H \\ V \end{matrix}),$ H是测压管水头，V是平均截面速率；f为单元界面处的通量；为源项；f为达西-威斯巴哈摩阻系数；D为管径。Among them, the superscripts n and n+1 represent t and t+Δt time steps, respectively; is the average value of u in the whole control body; $u = (\begin{matrix} h \\ V \end{matrix}),$ H is the piezometric head, V is the average section velocity; f is the flux at the cell interface; is the source item; f is the Darcy-Weissbach friction coefficient; D is the pipe diameter.

进一步地，步骤3包含以下子步骤：Further, step 3 includes the following sub-steps:

步骤3.1：求解水体内部控制单元界面处通量。Step 3.1: Solve the flux at the interface of the control unit inside the water body.

首先，基于黎曼问题，根据Godunov格式，对任一内部控制单元i(1<i<N)，界面i+1/2处的通量为：First, based on the Riemann problem, according to the Godunov scheme, for any internal control unit i (1<i<N), the flux at interface i+1/2 is:

${f f}_{i i + + 11 / / 22} = = {\overset{&OverBar; &OverBar;}{A A}}_{i i + + 11 / / 22} {u u}_{i i + + 11 / / 22} ((t t)) = = \frac{11}{22} {\overset{&OverBar; &OverBar;}{A A}}_{i i + + 11 / / 22} {{(\begin{matrix} 11 & a a / / g g \\ g g / / a a & 11 \end{matrix}) {U u}_{L L}^{n no} - - (\begin{matrix} - - 11 & a a / / g g \\ g g / / a a & - - 11 \end{matrix}) {U u}_{R R}^{n no}}} - - - - - - ((22))$

其中， $\begin{matrix} u = (\begin{matrix} H \\ V \end{matrix}); & \overset{&OverBar;}{A} = (\begin{matrix} \overset{&OverBar;}{V} & a^{2} / g \\ g & \overset{&OverBar;}{V} \end{matrix}); \end{matrix}$ H是测压管水头；V是平均截面速率；是V的平均值，为一常数；为在n时步时，u分别到界面i+1/2左、右侧两侧的平均值。in, $\begin{matrix} u = (\begin{matrix} h \\ V \end{matrix}); & \overset{&OverBar;}{A} = (\begin{matrix} \overset{&OverBar;}{V} & a^{2} / g \\ g & \overset{&OverBar;}{V} \end{matrix}); \end{matrix}$ H is the piezometric head; V is the average section velocity; is the average value of V and is a constant; It is the average value of u to the left and right sides of the interface i+1/2 at n time steps.

接着，通过引入MUSCL-Hancock格式计算内部单元通量f_i+1/2，从而取得二阶精度。Then, the internal unit flux f _i+1/2 is calculated by introducing the MUSCL-Hancock scheme, so as to obtain second-order precision.

更进一步地，步骤3.1中，需选择斜率限制器，以保证解中不出现虚假振荡。Furthermore, in step 3.1, a slope limiter needs to be selected to ensure that no spurious oscillations appear in the solution.

步骤3.2：构建虚拟控制单元以求解冲击水体上下游边界控制单元界面处通量。为在边界面处也取得二阶精度，分别在起始控制单元I₁上游侧、终点控制单元N下游侧构建两个虚拟控制单元I_-1、I₀，以及I_N+1、I_N+2，并假定在虚拟单元处的流动信息与边界处是一致的，从而可求解边界黎曼问题，且相应的Godunov通量f_1/2。和f_N+1/2也可像内部单元那样进行计算。Step 3.2: Construct a virtual control unit to solve the flux at the interface of the control unit at the upstream and downstream boundaries of the impacting water body. In order to obtain second-order precision at the boundary surface, two virtual control units I _-1 , I ₀ , and I _N+1 _, I _{N+ 2} , and assuming that the flow information at the virtual unit is consistent with that at the boundary, the boundary Riemann problem can be solved, and the corresponding Godunov flux f _1/2 . and _fN+1/2 can also be calculated like internal cells.

步骤3.3：将负特征线与黎曼向量相结合以求得管道进口边界控制单元界面处通量。Step 3.3: Combine the negative characteristic line with the Riemann vector to obtain the flux at the interface of the control unit at the pipe inlet boundary.

更进一步地，步骤4中，引入源项后，基于Runge-Kutta方法对水体的离散方程的数值积分项进行求解，瞬变流基本微分方程最终解的二阶有限体积法Godunov格式为：Furthermore, in step 4, after introducing the source term, the numerical integral term of the discrete equation of the water body is solved based on the Runge-Kutta method, and the second-order finite volume method Godunov format of the final solution of the transient flow basic differential equation is:

${U u}_{i i}^{n no + + 11} = = {\overset{&OverBar; &OverBar;}{U u}}_{i i}^{n no + + 11} + + Δ Δ t t s the s (({\overset{&OverBar; &OverBar;}{\overset{&OverBar; &OverBar;}{U u}}}_{i i}^{n no + + 11})) - - - - - - ((33))$

其中，为n+1时步，控制单元i在纯对流时，流动变量u的通量；为采用时间分裂法第一次更新后的通量。in, is the n+1 time step, when the control unit i is in pure convection, the flux of the flow variable u; is the flux after the first update using the time-splitting method.

更进一步地，步骤5中，对流项满足CFL条件(Courant-Friedrichs-Lewycriterion)，进一步可推得CFL条件下的最大时间步长Δt_max,CFL：Furthermore, in step 5, the convection item satisfies the CFL condition (Courant-Friedrichs-Lewycriterion), and the maximum time step Δt _max,CFL under the CFL condition can be deduced further:

${C C}_{r r} = = \frac{m m a a x x ((\overset{&OverBar; &OverBar;}{λ λ})) \cdot &Center Dot; Δ Δ t t}{Δ Δ x x} \leq \leq 11 &DoubleRightArrow; &DoubleRightArrow; {Δt Δt}_{max max,, C C F f L L} = = \frac{Δ Δ x x}{M m a a x x ((\overset{&OverBar; &OverBar;}{λ λ}))} - - - - - - ((44))$

其中，C_r为柯朗数，为矩阵的特征值。Among them, C _r is the Courant number, for the matrix eigenvalues of .

更进一步地，步骤5中，源项满足以下稳定性约束，并可推得适用于源项的最大时间步长Δt_max,s：Furthermore, in step 5, the source term satisfies the following stability constraints, and the maximum time step Δt _max,s applicable to the source term can be deduced:

$| | \frac{{\overset{&OverBar; &OverBar;}{\overset{&OverBar; &OverBar;}{U u}}}_{i i}^{n no + + 11}}{{\overset{&OverBar; &OverBar;}{U u}}_{i i}^{n no + + 11}} | | \leq \leq 11 &DoubleRightArrow; &DoubleRightArrow; {Δt Δt}_{max max,, s the s} = = min min ((- - 44 \frac{{\overset{&OverBar; &OverBar;}{U u}}_{i i}^{n no + + 11}}{s the s (({U u}_{i i}^{n no + + 11}))},, - - 22 \frac{{U u}_{i i}^{n no + + 11}}{s the s (({\overset{&OverBar; &OverBar;}{\overset{&OverBar; &OverBar;}{U u}}}_{i i}^{n no + + 11}))})) - - - - - - ((55))$

更进一步地，包含对流项和源项的最大允许时间步长为：Furthermore, the maximum allowable time step including convective and source terms is:

Δt_max＝min(Δt_max,CFL,Δt_max,s)(6)Δt _max =min(Δt _max,CFL ,Δt _max,s )(6)

进一步地，步骤6中，采用Godunov格式和理想气体状态方程联立实现水-气交界面的动态追踪：Further, in step 6, the dynamic tracking of the water-air interface is realized by using the Godunov scheme and the ideal gas state equation simultaneously:

首先，对尾水部分进行以下假定：(1)在整个瞬变过程中，尾水深度保持不变；(2)尾水长度随上游冲击水体的运动而变化；(3)尾水部分在瞬变过程中保持静止；First, the following assumptions are made on the tailwater: (1) the depth of the tailwater remains constant during the entire transient process; (2) the length of the tailwater changes with the movement of the upstream impacting water body; (3) the tailwater part remain stationary during the change;

其次，假定滞留气团的压缩、膨胀遵守理想气体状态方程变化规律；Secondly, it is assumed that the compression and expansion of the trapped air mass follow the law of the ideal gas state equation;

最后，基于上述给定条件，将正特征线和黎曼向量相结合，再联立水-气交界面动量和连续性方程、气团的控制方程，以求得冲击水体末端边界控制单元界面处通量；从而实现水-气交界面的动态追踪。Finally, based on the above given conditions, the positive characteristic line and the Riemann vector are combined, and then the momentum and continuity equation of the water-air interface and the control equation of the air mass are combined to obtain the flow rate at the interface of the boundary control unit at the end of the impacting water body. amount; thus realizing the dynamic tracking of the water-air interface.

有益效果：本发明提供的管道内水-气耦合瞬变流的模拟方法，成功地克服了有限体积法在追踪水-气交界面的动态追踪的难题，并且该方法简单且易于实现；(2)非线性对流项很容易加入到解中，而这些项在现有的特征线法求解中通常是忽略的，因此它可将该模型的应用范围扩展到包含较大马赫数的问题；(3)该计算方法通过时间算子分裂方法，为模拟多维水锤瞬变提供了框架，而特征线法较难扩展到包含多维瞬变的流体问题。Beneficial effects: the simulation method of the water-air coupling transient flow in the pipeline provided by the present invention successfully overcomes the problem of dynamic tracking of the water-air interface by the finite volume method, and the method is simple and easy to implement; (2 ) nonlinear convection terms are easily added to the solution, and these terms are usually ignored in the existing solution of the characteristic line method, so it can extend the application range of the model to problems involving larger Mach numbers; (3 ) This calculation method provides a framework for simulating multidimensional water hammer transients through the time operator splitting method, while the characteristic line method is difficult to extend to fluid problems involving multidimensional transients.

附图说明Description of drawings

图1为本发明的基本流程图；Fig. 1 is the basic flowchart of the present invention;

图2为实施例的输水管道系统示意图；Fig. 2 is the schematic diagram of the water pipeline system of the embodiment;

图3为实施例的网格划分示意图；Fig. 3 is the grid division schematic diagram of embodiment;

图4为实施例下，水-气耦合作用瞬变过程气团的压力变化；Fig. 4 is under the embodiment, the pressure change of the air mass in the transient process of water-air coupling;

图中：1-上游水库；2-输水管；3-阀门；4-滞留气团；5-冲击水体；6-尾水。In the figure: 1-upstream reservoir; 2-water delivery pipe; 3-valve; 4-retained air mass; 5-impact water body; 6-tail water.

具体实施方式detailed description

下面结合具体实施例对本发明作更进一步的说明。The present invention will be further described below in conjunction with specific examples.

如图1所示，一种基于有限体积法的管道内水-气耦合瞬变流的模拟方法，按以下步骤进行：将管道系统内瞬变流进行划分为水体、水-气交界面、滞留气团三大部分，并建立相应的控制方程，根据工程实例确定初始条件以及边界条件；根据有限体积法划分计算网格，并建立离散方程；采用Godunov格式求解水体的离散方程的数值通量，并取得二阶精度；通过基于Runge-Kutta方法对水体的离散方程的数值积分项进行求解，从而得到二阶显式有限体积法的Godunov格式；给出二阶显式有限体积法的Godunov格式所满足的CFL条件(Courant-Friedrichs-Lewycriterion)；采用Godunov格式和理想气体状态方程联立实现水-气交界面的动态追踪。As shown in Figure 1, a simulation method of water-air coupled transient flow in pipelines based on the finite volume method is carried out as follows: divide the transient flow in the pipeline system into water body, water-gas interface, stagnant The three major parts of the air mass, and establish the corresponding control equations, determine the initial conditions and boundary conditions according to the engineering examples; divide the calculation grid according to the finite volume method, and establish the discrete equation; use the Godunov scheme to solve the numerical flux of the discrete equation of the water body, and Obtain the second-order accuracy; solve the numerical integral term of the discrete equation of the water body based on the Runge-Kutta method, so as to obtain the Godunov scheme of the second-order explicit finite volume method; give the Godunov scheme of the second-order explicit finite volume method to satisfy The CFL condition (Courant-Friedrichs-Lewycriterion) is used; the dynamic tracking of the water-air interface is realized by using the Godunov scheme and the ideal gas state equation simultaneously.

下面将结合附图和实施例对本发明技术方案做进一步的详细描述。The technical solutions of the present invention will be further described in detail below in conjunction with the accompanying drawings and embodiments.

实施例:如图2所示的实验系统来研究水流快速冲击滞留气团的瞬变流现象。整个系统由上游压力水箱，有机透明管，阀门组成。管道总长8.83m，内径40mm，管道水平，阀门距离管道末端3.25m，上游入口压力为0.16MPa，管道末端封闭。初始时，阀门全关，阀后尾水水深为20mm。实验测得水锤波速为400m/s，管道平均摩阻0.075～0.095。水-气耦合作用过程由突然开启下游阀门引起，高速高清摄像机记录其从全关至全开时间为0.07s～0.09s。Embodiment: the experimental system as shown in Figure 2 is used to study the transient flow phenomenon that the water flow rapidly impacts the stranded air mass. The whole system is composed of upstream pressure water tank, organic transparent tube and valve. The total length of the pipeline is 8.83m, the inner diameter is 40mm, the pipeline is horizontal, the valve is 3.25m away from the end of the pipeline, the upstream inlet pressure is 0.16MPa, and the end of the pipeline is closed. Initially, the valve is fully closed, and the tail water depth behind the valve is 20mm. According to the experiment, the wave velocity of water hammer is 400m/s, and the average friction resistance of the pipeline is 0.075-0.095. The water-air coupling process is caused by the sudden opening of the downstream valve, and the high-speed high-definition camera records that the time from fully closed to fully opened is 0.07s to 0.09s.

本发明的具体步骤为：Concrete steps of the present invention are:

步骤1：将管道系统内瞬变流进行划分为水体、水-气交界面、滞留气团三大部分，并建立相应的控制方程，根据工程实例确定初始条件以及边界条件。Step 1: Divide the transient flow in the pipeline system into three parts: water body, water-air interface, and trapped air mass, and establish corresponding control equations, and determine the initial conditions and boundary conditions according to engineering examples.

水体的基本微分方程为：The basic differential equation of water body is:

$\frac{\partial \partial H h}{\partial \partial t t} + + V V \frac{\partial \partial H h}{\partial \partial x x} + + \frac{{a a}^{22}}{g g} \frac{\partial \partial V V}{\partial \partial x x} = = 00 - - - - - - ((11))$

$\frac{\partial \partial V V}{\partial \partial t t} + + V V \frac{\partial \partial V V}{\partial \partial x x} + + g g \frac{\partial \partial H h}{\partial \partial x x} + + \frac{f f | | V V | | V V}{22 D D.} = = 00 - - - - - - ((22))$

其中，沿管线距离x与时间t是自变量；H(x,t)是测压管水头；V(x,t)是平均截面速率；g是重力加速度；a是波速；f为达西-威斯巴哈摩阻系数；D为管径。Among them, the distance along the pipeline x and time t are independent variables; H(x,t) is the piezometric head; V(x,t) is the average section velocity; g is the acceleration of gravity; a is the wave velocity; f is Darcy- Weissbach Mo resistance coefficient; D is the pipe diameter.

滞留气团的控制方程为：The governing equation of the trapped air mass is:

${H h}_{a a} {V V}_{a a}^{m m} = = {H h}_{a a 00} {V V}_{a a 00}^{m m} {orH or H}_{a a} {L L}_{a a}^{m m} = = {H h}_{a a 00} {L L}_{a a 00}^{m m} - - - - - - ((33))$

其中，H_a、V_a、L_a分别为气团瞬变压力、体积和长度；其对应的初始值分别为H_a0、V_a0、L_a0；m为理想气体状态方程的多变指数，由于实验中瞬变过程极为短暂，仅为数秒钟，气团的压缩膨胀可视为绝热过程，m为1.4。Among them, H _a , V _a , and L _a are the transient pressure, volume, and length of the air mass, respectively; their corresponding initial values are H _a0 , V _a0 , and L _a0 ; m is the variable exponent of the ideal gas state equation. The medium transient process is very short, only a few seconds, and the compression and expansion of the air mass can be regarded as an adiabatic process, and m is 1.4.

水-气交界面的控制方程为：The governing equation for the water-air interface is:

$\frac{{dL L}_{f f}}{d d t t} = = {V V}_{w w} - - - - - - ((44))$

H_w＝H_a(5)H _w = H _a (5)

其中，L_f为冲击水体的长度；V_w、H_w分别水-气交界面处流速和压力。Among them, L _f is the length of the impacting water body; V _w and H _w are the flow velocity and pressure at the water-air interface, respectively.

本实施例下，计算域为从上游水库出水口至阀门之间的管道；初始条件为阀门全关时，初始流速为0，阀门上游水体初始压力等于上游水库静压水头，阀门下游滞留气团初始压力为大气压力，尾水压力考虑重力作用；边界条件为：管道入口处，水库提供恒压边界，且等于上游水库静压水头；下游球阀快速开启建立水力瞬变，管道末端封闭。In this embodiment, the calculation domain is the pipeline from the water outlet of the upstream reservoir to the valve; the initial condition is that when the valve is fully closed, the initial flow rate is 0, the initial pressure of the water body upstream of the valve is equal to the static pressure head of the upstream reservoir, and the initial condition of the trapped air mass downstream of the valve is The pressure is atmospheric pressure, and the tail water pressure considers the effect of gravity; the boundary conditions are: at the inlet of the pipeline, the reservoir provides a constant pressure boundary, which is equal to the static pressure head of the upstream reservoir; the downstream ball valve opens quickly to establish a hydraulic transient, and the end of the pipeline is closed.

步骤2：根据有限体积法划分计算网格，并建立离散方程。Step 2: Divide the calculation grid according to the finite volume method, and establish discrete equations.

图3为本实施例下水锤区网格离散系统。对第i个控制体，定义其上、下游界面编号分别为i-1/2、i+1/2。Fig. 3 is the grid discrete system in the water hammer area of this embodiment. For the i-th control body, define its upstream and downstream interface numbers as i-1/2 and i+1/2, respectively.

建立离散方程的具体步骤为：The specific steps to establish a discrete equation are:

(a)微分方程(1)(2)的黎曼问题可近似为含常系数的线性双曲系统的黎曼问题：(a) The Riemann problem of differential equations (1)(2) can be approximated as the Riemann problem of a linear hyperbolic system with constant coefficients:

$\frac{\partial \partial u u}{\partial \partial t t} + + \frac{\partial \partial f f ((u u))}{\partial \partial x x} = = s the s ((u u)),, f f ((u u)) = = \overset{&OverBar; &OverBar;}{A A} u u - - - - - - ((66))$

其中 $\begin{matrix} u = (\begin{matrix} H \\ V \end{matrix}); & s = (\begin{matrix} 0 \\ - \frac{f | V | V}{2 D} \end{matrix}); & \overset{&OverBar;}{A} = (\begin{matrix} \overset{&OverBar;}{V} & a^{2} / g \\ g & \overset{&OverBar;}{V} \end{matrix}); \end{matrix}$ 是V的平均值，为一常数。in $\begin{matrix} u = (\begin{matrix} h \\ V \end{matrix}); & the s = (\begin{matrix} 0 \\ - \frac{f | V | V}{2 D.} \end{matrix}); & \overset{&OverBar;}{A} = (\begin{matrix} \overset{&OverBar;}{V} & a^{2} / g \\ g & \overset{&OverBar;}{V} \end{matrix}); \end{matrix}$ is the average value of V and is a constant.

(b)在控制单元i(从界面i-1/2到i+1/2)及时间段Δt(从t到t+Δt)内对方程(6)积分，可以得到流动变量u的离散方程：(b) Integrating equation (6) in control unit i (from interface i-1/2 to i+1/2) and time period Δt (from t to t+Δt), the discrete equation of flow variable u can be obtained :

${U u}_{i i}^{n no + + 11} = = {U u}_{i i}^{n no} + + \frac{Δ Δ t t}{Δ Δ x x} (({f f}_{i i + + 11 / / 22} - - {f f}_{i i - - 11 / / 22})) + + \frac{Δ Δ t t}{Δ Δ x x} {&Integral; &Integral;}_{i i - - 11 / / 22}^{i i + + 11 / / 22} s the s d d x x - - - - - - ((77))$

其中，上标n和n+1分别代表t和t+Δt时步。为u在整个控制体的平均值；f为单元界面处的通量；为源项；where the superscripts n and n+1 represent t and t+Δt time steps, respectively. is the average value of u in the whole control body; f is the flux at the unit interface; is the source item;

步骤3：采用Godunov格式求解水体的离散方程的数值通量，并取得二阶精度。Step 3: Use the Godunov scheme to solve the numerical flux of the discrete equation of the water body, and obtain second-order accuracy.

步骤3.1：求解水体内部控制单元界面处通量Step 3.1: Solve the flux at the interface of the control unit inside the water body

首先，根据Godunov方法，其黎曼问题为以下初值问题：First, according to the Godunov method, its Riemann problem is the following initial value problem:

$\frac{\partial \partial u u}{\partial \partial t t} + + \frac{\partial \partial f f ((u u))}{\partial \partial x x} = = 00 - - - - - - ((88))$

${u u}^{n no} ((x x)) = = \{\begin{matrix} {U u}_{L L}^{n no},, & x x < < {x x}_{i i + + 11 / / 22} \\ {U u}_{R R}^{n no},, & x x > > {x x}_{i i + + 11 / / 22} \end{matrix} - - - - - - ((99))$

其中，为在n时步时，u分别到界面i+1/2左、右侧两侧的平均值。in, It is the average value of u to the left and right sides of the interface i+1/2 at n time steps.

对矩阵进行简化计算，可以得到其特征值及特征向量分别为：pair matrix Carrying out simplified calculations, its eigenvalues and eigenvectors can be obtained as follows:

$\begin{matrix} \overset{&OverBar; &OverBar;}{{λ λ}_{11}} = = \overset{&OverBar; &OverBar;}{V V} - - a a,, & \overset{&OverBar; &OverBar;}{{λ λ}_{22}} = = \overset{&OverBar; &OverBar;}{V V} + + a a \end{matrix} - - - - - - ((1010))$

${K K}^{((11))} = = [\begin{matrix} 11 \\ - - g g / / a a \end{matrix}] a a n no d d {K K}^{((22))} = = [\begin{matrix} 11 \\ g g / / a a \end{matrix}] - - - - - - ((1111))$

因特征向量线性无关，进一步可得到：Since the eigenvectors are linearly independent, it can be further obtained:

${U u}_{L L} = = {Σ Σ}_{i i = = 11}^{22} {α α}_{i i} {K K}^{((i i))} a a n no d d {U u}_{R R} = = {Σ Σ}_{i i = = 11}^{22} {β β}_{i i} {K K}^{((i i))} - - - - - - ((1212))$

求解四个未知系数α₁，α₂，β₁，β₂，可得：Solving the four unknown coefficients α ₁ , α ₂ , β ₁ , β ₂ , we can get:

$\begin{matrix} {α α}_{11} = = \frac{11}{22} (({H h}_{L L}^{n no} - - \frac{a a}{g g} {V V}_{L L}^{n no})),, & {α α}_{22} = = \frac{11}{22} (({H h}_{L L}^{n no} + + \frac{a a}{g g} {V V}_{L L}^{n no})) \end{matrix} - - - - - - ((1313))$

$\begin{matrix} {β β}_{11} = = \frac{11}{22} (({H h}_{R R}^{n no} - - \frac{a a}{99} {V V}_{R R}^{n no})),, & {β β}_{22} = = \frac{11}{22} (({H h}_{R R}^{n no} + + \frac{a a}{99} {V V}_{R R}^{n no})) \end{matrix} - - - - - - ((1414))$

黎曼问题(方程(8)和(9))一般解的原始变量形式为：The original variable form of the general solution to the Riemann problem (equations (8) and (9)) is:

u(x,t)＝β₁K⁽¹⁾+a₂K⁽²⁾(15)u(x,t)=β ₁ K ⁽¹⁾ +a ₂ K ⁽²⁾ (15)

现在，利用方程(15)可得到变量在界面i+1/2处的精确解：Now, the exact solution of the variable at interface i+1/2 can be obtained by using equation (15):

${u u}_{i i + + 11 / / 22} ((t t)) = = \frac{11}{22} {{(\begin{matrix} 11 & a a / / g g \\ g g / / a a & 11 \end{matrix}) {U u}_{L L}^{n no} - - (\begin{matrix} - - 11 & a a / / g g \\ g g / / a a & - - 11 \end{matrix}) {U u}_{R R}^{n no}}} - - - - - - ((1616))$

从而，对t∈[tⁿ,tⁿ⁺¹]，任一内部控制单元i(1<i<N)，在界面i+1/2处的通量为：Thus, for t∈[t ⁿ ,t ⁿ⁺¹ ], for any internal control unit i(1<i<N), the flux at interface i+1/2 is:

${f f}_{i i + + 11 / / 22} = = {\overset{&OverBar; &OverBar;}{A A}}_{i i + + 11 / / 22} {u u}_{i i + + 11 / / 22} ((t t)) = = \frac{11}{22} {\overset{&OverBar; &OverBar;}{A A}}_{i i + + 11 / / 22} {{(\begin{matrix} 11 & a a / / g g \\ g g / / a a & 11 \end{matrix}) {U u}_{L L}^{n no} - - (\begin{matrix} - - 11 & a a / / g g \\ g g / / a a & - - 11 \end{matrix}) {U u}_{R R}^{n no}}} - - - - - - ((1717))$

接着，引入MUSCL-Hancock格式来计算内部单元通量f_i+1/2，取得二阶精度。具体步骤为：Then, the MUSCL-Hancock scheme is introduced to calculate the internal unit flux f _i+1/2 to obtain second-order precision. The specific steps are:

(a)数据重构：采用每个控制单元[x_i-1/2,x_i+1/2]内的分段线性函数代替数据单元平均值则在极值点处的值为：(a) Data reconstruction: use the piecewise linear function in each control unit [ _xi-1/2 , _xi+1/2 ] instead of the average value of the data unit then at the extremum point The value is:

$\begin{matrix} {U u}_{i i}^{L L} = = {U u}_{i i}^{n no} - - \frac{Δ Δ x x}{22} {Δ Δ}_{i i},, & {U u}_{i i}^{R R} = = {U u}_{i i}^{n no} + + \frac{Δ Δ x x}{22} {Δ Δ}_{i i} \end{matrix} - - - - - - ((1818))$

其中，Δ_i为选定的适度斜坡向量，用以增加计算格式的精度，并保证解中不出现虚假振荡。本发明中，Δ_i选择MINMOD限制器：where _Δi is a moderately sloped vector selected to increase the accuracy of the calculation scheme and to ensure that no spurious oscillations appear in the solution. In the present invention, _Δi selects the MINMOD limiter:

${Δ Δ}_{i i} = = M m I I N N M m O o D D. (({σ σ}_{i i}^{n no},, {σ σ}_{i i - - 11}^{n no})) = = \{\begin{matrix} {σ σ}_{i i}^{n no},, & i i f f,, | | {σ σ}_{i i}^{n no} | | < < | | {σ σ}_{i i - - 11}^{n no} | |,, a a n no d d,, {σ σ}_{i i}^{n no} {σ σ}_{i i - - 11}^{n no} > > 00 \\ {σ σ}_{i i - - 11}^{n no},, & i i f f,, | | {σ σ}_{i i}^{n no} | | > > | | {σ σ}_{i i - - 11}^{n no} | |,, a a n no d d,, {σ σ}_{i i}^{n no} {σ σ}_{i i - - 11}^{n no} > > 00 \\ 00,, & i i f f,, {σ σ}_{i i}^{n no} {σ σ}_{i i - - 11}^{n no} < < 00 \end{matrix} - - - - - - ((1919))$

其中， $\begin{matrix} σ_{i}^{n} = (U_{i + 1}^{n} - U_{i}^{n}) / Δ x, & σ_{i - 1}^{n} = (U_{i}^{n} - U_{i - 1}^{n}) / Δ x . \end{matrix}$ in, $\begin{matrix} σ_{i}^{no} = (u_{i + 1}^{no} - u_{i}^{no}) / Δ x, & σ_{i - 1}^{no} = (u_{i}^{no} - u_{i - 1}^{no}) / Δ x . \end{matrix}$

(b)推算：对每一个单元[x_i-1/2,x_i+1/2]，方程(18)中的边界外推值根据下式乘以0.5Δt推算：(b) Extrapolation: For each unit [x _i-1/2 ,x _i+1/2 ], the boundary extrapolation value in equation (18) Calculate according to the following formula multiplied by 0.5Δt:

$\begin{matrix} {\overset{&OverBar; &OverBar;}{U u}}_{i i}^{L L} = = {U u}_{i i}^{L L} + + \frac{11}{22} \frac{Δ Δ t t}{Δ Δ x x} [[f f (({U u}_{i i}^{L L})) - - f f (({U u}_{i i}^{R R}))]],, \\ {\overset{&OverBar; &OverBar;}{U u}}_{i i}^{R R} = = {U u}_{i i}^{R R} + + \frac{11}{22} \frac{Δ Δ t t}{Δ Δ x x} [[f f (({U u}_{i i}^{L L})) - - f f (({U u}_{i i}^{R R}))]] \end{matrix} - - - - - - ((2020))$

(c)黎曼问题：为计算内部界面通量f_i+1/2，需结合下面数据求解传统的黎曼问题：(c) Riemann problem: In order to calculate the internal interface flux f _i+1/2 , it is necessary to solve the traditional Riemann problem with the following data:

${U u}_{L L}^{n no} &equiv; &equiv; {\overset{&OverBar; &OverBar;}{U u}}_{i i}^{R R},, {U u}_{R R}^{n no} &equiv; &equiv; {\overset{&OverBar; &OverBar;}{U u}}_{i i}^{L L} - - - - - - ((21 twenty one))$

将方程(21)代入方程(17)，即可得到管中充满水时，通量在时间t＝[tⁿ,tⁿ⁺¹]内，所有内部单元界面i+1/2处的二阶格式。 ^Substituting Equation (21) into Equation ( ¹⁷ ), we can get the second-order Format.

步骤3.2：构建虚拟控制单元以求解冲击水体上下游边界控制单元界面处通量。为在边界面处也取得二阶精度，分别在起始控制单元I₁上游侧、终点控制单元N下游侧构建两个虚拟控制单元I_-1、I₀，以及I_N+1、I_N+2，并假定在虚拟单元处的流动信息与边界处是一致的。从而可求解边界黎曼问题，且相应的Godunov通量f_1/2。和f_N+1/2也可像内部单元那样进行计算。Step 3.2: Construct a virtual control unit to solve the flux at the interface of the control unit at the upstream and downstream boundaries of the impacting water body. In order to obtain second-order precision at the boundary surface, two virtual control units I _-1 , I ₀ , and I _N+1 _, I _{N+ 2} , and assume that the flow information at the virtual unit is consistent with that at the boundary. Thus the boundary Riemann problem can be solved, and the corresponding Godunov flux f _1/2 . and _fN+1/2 can also be calculated like internal cells.

本实施例中，上游水库恒水位边界：In this embodiment, the constant water level boundary of the upstream reservoir:

在上游边界处，与负特征线相关的黎曼不变量为：At the upstream boundary, the Riemann invariant associated with the negative characteristic line is:

${H h}_{11 / / 22} - - \frac{a a}{g g} {V V}_{11 / / 22} = = {H h}_{11}^{n no} - - \frac{a a}{g g} {V V}_{11}^{n no} - - - - - - ((22 twenty two))$

其中，H_1/2＝H_res，H_res为上游水库静压水头。Wherein, H _1/2 =H _res , and H _res is the hydrostatic head of the upstream reservoir.

从而可推得变量u_1/2(t)＝(H_1/2,V_1/2)。根据虚拟单元处的假定，临近管道进口的虚拟单元I_-1和I₀的相应值为：so that it can be deduced Variable u _1/2 (t)=(H _1/2 , V _1/2 ). According to the assumption at the virtual unit, the corresponding values of the virtual units I _-1 and I ₀ near the pipeline inlet are:

${U u}_{- - 11}^{n no + + 11} = = {U u}_{00}^{n no + + 11} = = {U u}_{11 / / 22} - - - - - - ((23 twenty three))$

步骤4：通过基于Runge-Kutta方法对水体的离散方程的数值积分项进行求解，从而得到二阶显式有限体积法的Godunov格式。Step 4: Solve the numerical integral term of the discrete equation of the water body based on the Runge-Kutta method to obtain the Godunov scheme of the second-order explicit finite volume method.

具体实现过程为：The specific implementation process is:

(a)纯对流时：(a) For pure convection:

${\overset{&OverBar; &OverBar;}{U u}}_{i i}^{n no + + 11} = = {U u}_{i i}^{n no} - - \frac{Δ Δ t t}{Δ Δ x x} (({f f}_{i i + + 11 / / 22}^{n no} - - {f f}_{i i - - 11 / / 22}^{n no})) - - - - - - ((24 twenty four))$

(b)利用源项乘以Δt/2进行更新：(b) Update by multiplying the source term by Δt/2:

${\overset{&OverBar; &OverBar;}{\overset{&OverBar; &OverBar;}{U u}}}_{i i}^{n no + + 11} = = {\overset{&OverBar; &OverBar;}{U u}}_{i i}^{n no + + 11} + + \frac{Δ Δ t t}{22} s the s (({\overset{&OverBar; &OverBar;}{U u}}_{i i}^{n no + + 11})) - - - - - - ((2525))$

(c)利用源项乘以Δt再次更新：(c) Update again by multiplying the source term by Δt:

${U u}_{i i}^{n no + + 11} = = {\overset{&OverBar; &OverBar;}{U u}}_{i i}^{n no - - 11} + + Δ Δ t t s the s (({\overset{&OverBar; &OverBar;}{\overset{&OverBar; &OverBar;}{U u}}}_{i i}^{n no + + 11})) - - - - - - ((2626))$

步骤5：给出二阶显式有限体积法的Godunov格式所满足的CFL条件(Courant-Friedrichs-Lewycriterion)。Step 5: Give the CFL condition (Courant-Friedrichs-Lewycriterion) satisfied by the Godunov scheme of the second-order explicit finite volume method.

CFL条件下的最大时间步长△t_max,CFL：The maximum time step △t _max,CFL under CFL conditions:

${C C}_{r r} = = \frac{M m a a x x ((\overset{&OverBar; &OverBar;}{λ λ})) \cdot \cdot Δ Δ t t}{Δ Δ x x} \leq \leq 11 &DoubleRightArrow; &DoubleRightArrow; {Δt Δt}_{max max,, C C F f L L} = = \frac{Δ Δ x x}{M m a a x x ((\overset{&OverBar; &OverBar;}{λ λ}))} - - - - - - ((2727))$

其中，C_r为柯朗数。Among them, C _r is the Courant number.

此外，引入的源项满足以下稳定性约束：Furthermore, the introduced source term satisfies the following stability constraints:

$| | \frac{{\overset{&OverBar; &OverBar;}{\overset{&OverBar; &OverBar;}{U u}}}_{i i}^{n no + + 11}}{{\overset{&OverBar; &OverBar;}{U u}}_{i i}^{n no + + 11}} | | \leq \leq 11 - - - - - - ((2828))$

适用于源项的最大时间步长△t_max,s：The maximum time step size △t _max,s applicable to the source term:

${Δt Δt}_{max max,, s the s} = = m m i i n no ((- - 44 \frac{{\overset{&OverBar; &OverBar;}{U u}}_{i i}^{n no + + 11}}{s the s (({\overset{&OverBar; &OverBar;}{U u}}_{i i}^{n no + + 11}))},, - - 22 \frac{{U u}_{i i}^{n no + + 11}}{s the s (({\overset{&OverBar; &OverBar;}{\overset{&OverBar; &OverBar;}{U u}}}_{i i}^{n no + + 11}))})) - - - - - - ((2929))$

可推得包含对流项和源项的最大允许时间步长为：The maximum allowable time step that can be derived to include convective and source terms is:

△t_max＝min(△t_max,cFL，△t_max,s)(30)△t _max ＝min(△t _max,cFL ，△t _max,s )(30)

首先，对尾水部分进行以下假定：(1)在整个瞬变过程中，尾水深度保持不变；(2)尾水长度随上游冲击水体的运动而变化；(3)尾水部分在瞬变过程中保持静止。其次，假定滞留气团的压缩、膨胀遵守理想气体状态方程变化规律。最后，基于上述给定条件，将正特征线和黎曼向量相结合，再联立水-气交界面动量和连续性方程、气团的控制方程，以求得冲击水体末端边界控制单元界面处通量。从而实现水-气交界面的动态追踪。First, the following assumptions are made for the tail water: (1) the depth of the tail water remains constant during the entire transient process; (2) the length of the tail water changes with the movement of the upstream impacting water body; (3) the tail water part is in the transient Stay still during the change. Secondly, it is assumed that the compression and expansion of the trapped air mass follow the law of the ideal gas state equation. Finally, based on the above given conditions, the positive characteristic line and the Riemann vector are combined, and then the momentum and continuity equation of the water-air interface and the control equation of the air mass are combined to obtain the flow rate at the interface of the boundary control unit at the end of the impacting water body. quantity. In this way, the dynamic tracking of the water-air interface can be realized.

具体步骤如下：Specific steps are as follows:

在水-气交界面处，与负特征线相关的黎曼不变量为At the water-air interface, the Riemann invariant associated with the negative characteristic line is

${H h}_{A A - - W W} + + \frac{a a}{g g} {V V}_{A A - - W W} = = {H h}_{N N f f + + 11}^{n no} + + \frac{a a}{g g} {V V}_{N N f f + + 11}^{n no} - - - - - - ((3131))$

其中，为冲击水体终点控制单元N_f中压力和流速；H_A-W、V_A-W为冲击水体末端处压力和流速in, is the pressure and flow velocity in the control unit N _f at the end of the impact water body; H _AW and V _AW are the pressure and flow velocity at the end of the impact water body

水流冲击过程中，水-气交界面处连续性方程和动量方程为：During the impact of water flow, the continuity equation and momentum equation at the water-air interface are:

A·V_A-W＝(A-A_C)·V_w(32)A·V _AW ＝(AA _C )·V _w (32)

g·(A_C·H_C-A·H_A-W)+A·V_A-W·(V_w-V_C)＝0(33)g·(A _C ·H _C -A·H _AW )+A·V _AW ·(V _w -V _C )＝0(33)

其中，A为管道截面积；A_C为尾水湿周截面面积Among them, A is the cross-sectional area of the pipeline; A _C is the cross-sectional area of the wetted tail water

水体长度变化为：The length of the water body changes as:

${L L}_{f f} = = {L L}_{f f}^{n no} + + \frac{{V V}_{w w} + + {V V}_{w w}^{n no}}{22} Δ Δ t t - - - - - - ((3434))$

联立气团控制方程，可得：Simultaneous air mass governing equations can be obtained:

${H h}_{a a} = = {H h}_{a a 00} \cdot &Center Dot; {((\frac{{L L}_{a a 00}}{L L - - {L L}_{f f}^{n no} - - \frac{{V V}_{w w} + + {V V}_{w w}^{n no}}{22} Δ Δ t t}))}^{m m} - - - - - - ((3535))$

将上述5个方程联立求解，即可得到H_A-W,V_A-W,H_A,H_C和V_w。Solve the above five equations simultaneously to get H _AW , V _AW , H _A , H _C and V _w .

根据虚拟单元处的假定，临近I_N的虚拟单元I_N+1和I_N+2处的相应值为：According to the assumption at the virtual unit, the corresponding values at the virtual units I _N+1 and I _N+2 adjacent to I _N are:

${U u}_{N N f f + + 11}^{n no + + 11} = = {U u}_{N N f f + + 22}^{n no + + 11} = = (\begin{matrix} {H h}_{A A - - W W} \\ {V V}_{A A - - W W} \end{matrix}) - - - - - - ((3636))$

在进行下一时步计算前，需对冲击水体部分的计算网格进行校核和网格重构，使得末端微段水体满足0<ΔL_f≤Δx，其中ΔL_f＝L_f-N_fΔx。当滞留气团压缩或膨胀时，可能会出现Δx<ΔL_f或ΔL_f≤0情况。为了满足0<ΔL_f≤Δx要求，需要增加或减少计算网格。Before the calculation of the next time step, it is necessary to check and reconstruct the calculation grid of the impacted water body, so that the terminal micro-segment water body satisfies 0<ΔL _f ≤ Δx, where ΔL _f =L _f -N _f Δx. When the trapped air mass is compressed or expanded, Δx < ΔL _f or ΔL _f ≤ 0 may occur. In order to meet the requirements of 0<ΔL _f ≤Δx, it is necessary to increase or decrease the calculation grid.

通过以上方法编程计算后，将计算结果均与实验数据做对比。图4给出了本实施例下水流冲击滞留气团过程中气团压力曲线图，同时也给出了忽略尾水影响的计算结果(将滞留气团等体积折算为满管滞留气团的情况，忽略尾水与冲击水体相互作用)。可以看出，本发明方法计算的压力的幅值和时间响应均与实验吻合良好。同时，可以看出忽略尾水影响会引起压力幅值和周期均有较大的计算误差。After programming and calculating by the above method, the calculation results are compared with the experimental data. Fig. 4 has provided the air mass pressure curve diagram in the process of the water flow impacting the trapped air mass under the present embodiment, and also provided the calculation result of ignoring the influence of the tail water (converting the equal volume of the trapped air mass to the situation of the full pipe trapped air mass, ignoring the tail water interact with the impacting water body). It can be seen that the amplitude and time response of the pressure calculated by the method of the present invention are in good agreement with the experiment. At the same time, it can be seen that ignoring the tail water effect will cause large calculation errors in both the pressure amplitude and period.

以上所述仅是本发明的优选实施方式，应当指出：对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。The above is only a preferred embodiment of the present invention, it should be pointed out that for those of ordinary skill in the art, without departing from the principle of the present invention, some improvements and modifications can also be made. It should be regarded as the protection scope of the present invention.

Claims

1. A simulation method of water-gas coupling transient flow in a pipeline is characterized by comprising the following steps: a finite volume method and a Godunov format are adopted to simulate the hydraulic transient phenomenon of water-gas coupling effect containing stagnant air mass in a pressure pipeline system, and the method comprises the following specific steps:

step 1: dividing transient flow in a pipeline system into three major parts, namely a water body, a water-air interface and a stagnant air mass, establishing a corresponding control equation, and determining initial conditions and boundary conditions according to an engineering example;

step 2: dividing a computational grid according to a finite volume method, and establishing a discrete equation;

and step 3: solving the numerical flux of a discrete equation of the water body by adopting a Godunov format, and obtaining second-order precision;

and 4, step 4: solving a numerical integral term of a discrete equation of the water body by a Runge-Kutta method, thereby obtaining a Godunov format of a second-order explicit finite volume method;

and 5: giving a CFL condition met by a Godunov format of a second-order explicit finite volume method;

step 6: and the Godunov format and an ideal gas state equation are adopted to realize the dynamic tracking of the water-gas interface simultaneously.

2. The method for simulating a water-gas coupled transient flow in a pipeline according to claim 1, wherein: in step 2, the integral equation of the flow variable u is established for the control unit i (from the interface i-1/2 to i +1/2) and the time period Δ t (from t to t + Δ t) as follows:

U_{i}^{n + 1} = U_{i}^{n} + \frac{Δ t}{Δ x} (f_{i + 1 / 2} - f_{i - 1 / 2}) + \frac{Δ t}{Δ x} {&Integral;}_{i - 1 / 2}^{i + 1 / 2} s d x - - - (1)

wherein, superscripts n and n +1 represent t and t + Δ t time steps, respectively;is the average value of u over the entire control volume;

u = (\begin{matrix} H \\ V \end{matrix}),

h is the piezometer tube head, V is the average section velocity; f is the flux at the cell interface;

s = (\begin{matrix} 0 \\ - \frac{f | V | V}{2 D} \end{matrix}),

is a source item; f is the Darcy-Weisbaha friction coefficient; d is the pipe diameter.

3. The method for simulating a water-gas coupled transient flow in a pipeline according to claim 1, wherein: the step 3 is specifically as follows:

step 3.1: solving flux at the interface of the control unit in the water body:

first, based on the riemann problem, according to the Godunov format, for any internal control unit i (1< i < N), the flux at interface i +1/2 is:

f_{i + 1 / 2} = {\overset{&OverBar;}{A}}_{i + 1 / 2} u_{i + 1 / 2} (t) = \frac{1}{2} {\overset{&OverBar;}{A}}_{i + 1 / 2} {(\begin{matrix} 1 & a / g \\ g / a & 1 \end{matrix}) U_{L}^{n} - (\begin{matrix} - 1 & a / g \\ g / a & - 1 \end{matrix}) U_{R}^{n}} - - - (2)

wherein,

u = (\begin{matrix} H \\ V \end{matrix}); \overset{&OverBar;}{A} = (\begin{matrix} \overset{&OverBar;}{V} & a^{2} / g \\ g & \overset{&OverBar;}{V} \end{matrix});

h is a piezometer tube water head; v is the average cross-sectional velocity;is the average value of V, which is a constant;the average value of u to the left side and the right side of the interface i +1/2 respectively when n time steps occur;

next, the internal element flux f is calculated by introducing MUSCL-Hancock format_i+1/2Thereby obtaining second-order precision;

step 3.2: constructing a virtual control unit to solve the flux at the interface of the upstream and downstream boundary control units of the impact water body: to achieve a second-order accuracy also at the boundary surfaces, respectively at the start control unit I₁Two virtual control units I are constructed on the upstream side and the downstream side of the end point control unit N_-1、I₀And I_N+1、I_N+2And assuming that the flow information at the virtual cell is consistent with that at the boundary, the boundary Riemann problem can be solved and the corresponding Godunov flux f_1/2And f_N+1/2The calculation is also performed like the internal unit;

step 3.3: and combining the negative characteristic line with the Riemann vector to obtain the flux at the interface of the pipeline inlet boundary control unit.

4. A method of simulating a water-gas coupled transient flow in a pipeline as claimed in claim 3, wherein: internal element flux f was calculated at the introduction of MUSCL-Hancock format_i+1/2In the process, the slope limiter is selected to ensure that spurious oscillations do not occur in the solution.

5. A method of simulating a water-gas coupled transient flow in a pipeline as claimed in claim 3, wherein: in the step 4, after the source term is introduced, the numerical integral term of the discrete equation of the water body is solved based on a Runge-Kutta method, and the Godunov format of the final solution of the transient flow basic differential equation is as follows:

U_{i}^{n + 1} = {\overset{&OverBar;}{U}}_{i}^{n + 1} + Δ t s ({\overset{&OverBar;}{\overset{&OverBar;}{U}}}_{i}^{n + 1}) - - - (3)

wherein,in the time step of n +1, the control unit i flows the flux of the variable u in the pure convection;the flux after the first update using time-splitting.

6. The method for simulating a water-gas coupled transient flow in a pipeline according to claim 5, wherein: in the step 5, the Godunov format of the second-order explicit finite volume method satisfies the CFL condition, and the maximum time step Δ t under the CFL condition is obtained_max,CFL：

C_{r} = \frac{Max (\overset{&OverBar;}{λ}) \cdot Δt}{Δx} \leq 1 &DoubleRightArrow; {Δt}_{\max, CFL} = \frac{Δx}{Max (\overset{&OverBar;}{λ})} - - - (4) .

7. The method for simulating a water-gas coupled transient flow in a pipeline according to claim 6, wherein: the source term satisfies the following stability constraint and can be derived as the maximum time step Δ t applicable to the source term_max,s：

| \frac{{\overset{&OverBar;}{\overset{&OverBar;}{U}}}_{i}^{n + 1}}{{\overset{&OverBar;}{U}}_{i}^{n + 1}} | \leq 1 &DoubleRightArrow; {Δt}_{\max, s} = \min (- 4 \frac{{\overset{&OverBar;}{U}}_{i}^{n + 1}}{s ({\overset{&OverBar;}{U}}_{i}^{n + 1})}, - 2 \frac{U_{i}^{n + 1}}{s ({\overset{&OverBar;}{\overset{&OverBar;}{U}}}_{i}^{n + 1})}) - - - (5) .

8. The method for simulating a water-gas coupled transient flow in a pipeline according to claim 1, wherein: the maximum allowed time step containing the convection item and the source item is:

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

9. the method for simulating a water-gas coupled transient flow in a pipeline according to claim 1, wherein: the step 6 specifically includes:

first, the following assumptions were made for the tail water fraction: (1) in the whole transient process, the tail water depth is kept unchanged; (2) the length of tail water changes along with the movement of upstream impact water body; (3) the tailwater portion remains stationary during transients;

secondly, the compression and expansion of the retained gas mass are assumed to follow the change rule of the ideal gas state equation;

finally, combining the positive characteristic line with the Riemann vector based on given conditions, and then simultaneously establishing a water-gas interface momentum and continuity equation and a control equation of an air mass to obtain the flux at the interface of the control unit of the boundary at the tail end of the impact water body; thereby realizing the dynamic tracking of the water-gas interface.