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

Abstract

The invention discloses an analogy method of water-gas coupling transient flow in a pipeline on the basis of a finite volume method, and considers the influence of the convective term and the tailwater depth of an impact water body in a water-gas coupling function process. Firstly, a computational domain is divided into three parts including the impact water body, a water-gas interface and a retention air mass. Firstly, the water body computational domain is firstly dispersed into a computation grid unit by aiming at a water body part, a discrete format of a transient flow basic differential equation is established, and then, a Godunov format is adopted to solve the numerical flux of a water body unit control volume. A Runge-Kutta method is adopted to solve the numerical value integral term of a water body discrete equation. A Godunov format and an ideal gas state equation are simultaneous to realize the dynamic tracking of a water-gas interface. A slope restrictor is introduced to restrict false numerical value oscillation so as to realize two-order precision simulation. The water-gas coupling function hydraulic phenomenon of a stream impact trapped air mass in a water pipeline system can be subjected to strict and accurate numerical value analogy and transient analysis.

Description

Simulation method of water-gas coupling transient flow in pipeline

Technical Field

The invention relates to a finite volume method-based simulation method of water-gas coupling transient flow in a pipeline, belonging to the technical field of hydraulic numerical calculation of hydropower stations (pump stations).

Background

In water pipe systems such as actual hydropower stations, pump stations and the like, air masses may be retained at positions such as a convex section, a closed end, an emptying maintenance section and the like along the water pipe systems, a complex hydraulic transient phenomenon of water-air coupling may occur in the process of starting and stopping the systems and converting working conditions, and corresponding abnormal water hammers (or called as gas-water hammers) may be large enough to cause pipe blasting. Many of the destructive incidents of actual water delivery systems are related to such transient flows containing stagnant air pockets. However, to date, the design criteria for practical piping systems have only considered the case of a full pipe, and not the presence of stagnant air pockets and their hazards; for the condition of containing the retention air mass, no corresponding calculation standard exists in the pipeline design, and the existing results are not complete. Therefore, in order to ensure the safe operation of the water pipeline system, the intensive research is carried out on the transient flow phenomenon caused by the retained air mass, and the method has great practical significance.

The scholars have studied the complex transient flow and put forward a corresponding simulation method. The existing mathematical model is basically a one-dimensional model and mainly comprises an elastic water body model based on a characteristic method and a rigid water body model based on a rigidity theory. Compared with an elastic water body model, a rigid water body model has the advantages of simplicity, quickness and the like, but is only suitable for pressure prediction of a large-volume trapped air mass situation. The existing elastic water model is mainly obtained by a characteristic line method, and the pressure change of transient processes such as most simple pipelines can be effectively predicted. When a pipeline system in which the wave velocity of the pipeline varies is involved, the characteristic line method is very troublesome to handle. In addition, for the simulation of the water-gas coupling transient flow of the water flow impact retention air mass, the existing models neglect the influence of the momentum equation and the continuous equation on the convection term, because: (1) for the problem of hydraulic transient with slower flow, the influence on the flow term is small; (2) the solution difficulty of the characteristic line method is greatly increased by considering the convection term. The complex water-gas coupling transient flow of the water flow fast impacting the stagnant air mass relates to high Mach number and even violent flow velocity change, and the convection term of the transient flow has certain influence on the flow.

Finite volume methods have been widely used for solving hyperbolic systems, such as: aerodynamic aspects. The method is well able to respect momentum and mass conservation and provides an effective simulation of discontinuity problems. There have been attempts by scholars to use the limited volume method for the study of the problem of water hammer in simple pipes. However, there is no report on the application of the limited volume method for complex transient flows of water-air coupling where water flow rapidly impacts stagnant air mass, and considering the tail water situation.

Therefore, if the influence of convection terms, tail water and other factors in the water-gas coupling process can be fully considered, and new methods and skills are introduced to eliminate defects, the continuous improvement of the numerical simulation method of the transient pressure of the water-gas coupling in the process of water flow impacting the trapped air mass can be facilitated.

Disclosure of Invention

The purpose is as follows: in order to make up for the defects of complex interpolation and the like of the existing elastic water body model, the invention provides a simulation method of water-gas coupling transient flow in a pipeline, which is based on a finite volume method, has simple algorithm and easy realization, can fully consider the factors such as convection terms, tail water and the like, and avoids the problems of fussy interpolation and the like.

The technical scheme is as follows: in order to solve the technical problems, the technical scheme adopted by the invention is as follows:

a simulation method of water-gas coupling transient flow in a pipeline based on a finite volume method adopts the finite volume method and Godunov format to simulate the hydraulic transient phenomenon of water-gas coupling effect containing stagnant air mass in a pressure pipeline system, and 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 the CFL condition (Courant-Friedrichs-Lewycrion) satisfied by the Godunov format of the 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.

Furthermore, in step 2, the numbers of the upstream and downstream interfaces of the ith control unit are defined as i-1/2 and i +1/2 respectively.

Further, in step 2, for the control unit i, the integral equation of the flow variable u is established as follows:

$U_i^{n+1}=U_i^n+\frac{\mathrm{\Δ}t}{\mathrm{\Δ}x}(f_{i+1/2}-f_{i-1/2})+\frac{\mathrm{\Δ}t}{\mathrm{\Δ}x}{\∫}_{i-1/2}^{i+1/2}sdx---\left(1\right)$

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=\left(\begin{array}{c}H\\ V\end{array}\right),$ h is the piezometer tube head, V is the average section velocity; f is the flux at the cell interface;is a source item; f is the Darcy-Weisbaha friction coefficient; d is the pipe diameter.

Further, step 3 comprises the following substeps:

step 3.1: and solving the 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}={\stackrel{\&OverBar;}{A}}_{i+1/2}{u}_{i+1/2}\left(t\right)=\frac{1}{2}{\stackrel{\&OverBar;}{A}}_{i+1/2}\{\left(\begin{array}{cc}1& a/g\\ g/a& 1\end{array}\right)U_L^n-\left(\begin{array}{cc}-1& a/g\\ g/a& -1\end{array}\right)U_R^n\}---\left(2\right)$

wherein, $\begin{array}{cc}u=\left(\begin{array}{c}H\\ V\end{array}\right);& \stackrel{\&OverBar;}{A}=\left(\begin{array}{cc}\stackrel{\&OverBar;}{V}& a^2/g\\ g& \stackrel{\&OverBar;}{V}\end{array}\right);\end{array}$ h is a piezometer tube water head; v is the average cross-sectional velocity;is the average value of V, which is a constant;is the average value of u to the left and right sides of the interface i +1/2 at n time steps.

Next, the internal element flux f is calculated by introducing MUSCL-Hancock format_i+1/2Thereby achieving second order accuracy.

Further, in step 3.1, a slope limiter is selected to ensure that no spurious oscillations occur in the solution.

Step 3.2: and constructing a virtual control unit to solve the flux at the interface of the control unit of the upstream and downstream boundaries 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/2. And f_N+1/2The calculation can also be done like an 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.

Furthermore, in step 4, after the source term is introduced, the numerical integral term of the discrete equation of the water body is solved based on the 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}={\stackrel{\&OverBar;}{U}}_i^{n+1}+\mathrm{\&Delta;}ts\left({\stackrel{\&OverBar;}{\stackrel{\&OverBar;}{U}}}_i^{n+1}\right)---\left(3\right)$

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.

Further, in step 5, the convection term satisfies the CFL condition (Courant-Friedrichs-Lewycriterion), and the maximum time step Δ t under the CFL condition can be further deduced_max,CFL：

$C_r=\frac{max\left(\stackrel{\&OverBar;}{\mathrm{\&lambda;}}\right)\&CenterDot;\mathrm{\&Delta;}t}{\mathrm{\&Delta;}x}\&le;1\&DoubleRightArrow;{\mathrm{\&Delta;t}}_{\max ,CFL}=\frac{\mathrm{\&Delta;}x}{Max\left(\stackrel{\&OverBar;}{\mathrm{\&lambda;}}\right)}---\left(4\right)$

Wherein, C_rIs a function of the number of the korangian,is a matrixThe characteristic value of (2).

Further, in step 5, the source term satisfies the following stability constraint, and the maximum time step Δ t applicable to the source term can be derived_max,s：

$\left|\frac{{\stackrel{\&OverBar;}{\stackrel{\&OverBar;}{U}}}_i^{n+1}}{{\stackrel{\&OverBar;}{U}}_i^{n+1}}\right|\&le;1\&DoubleRightArrow;{\mathrm{\&Delta;t}}_{\max ,s}=\min \left(-4\frac{{\stackrel{\&OverBar;}{U}}_i^{n+1}}{s\left(U_i^{n+1}\right)},-2\frac{U_i^{n+1}}{s\left({\stackrel{\&OverBar;}{\stackrel{\&OverBar;}{U}}}_i^{n+1}\right)}\right)---\left(5\right)$

Further, the maximum allowed time step including the convection item and the source item is:

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

further, in step 6, the dynamic tracking of the water-gas interface is realized by adopting Godunov format and ideal gas equation of state in a simultaneous manner:

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, based on the given conditions, combining the positive characteristic line with the Riemann vector, and then combining the momentum and continuity equation of the water-gas interface and the control equation of the 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.

Has the advantages that: the simulation method of the water-gas coupling transient flow in the pipeline provided by the invention successfully overcomes the difficulty of dynamic tracking of a finite volume method on a tracking water-gas interface, and is simple and easy to realize; (2) nonlinear convection terms are easily added into the solution, and the terms are usually ignored in the existing solution of the characteristic line method, so that the application range of the model can be expanded to the problem containing large Mach number; (3) the calculation method provides a frame for simulating the multidimensional water hammer transient through a time operator splitting method, and the characteristic line method is difficult to expand to the problem of the fluid containing the multidimensional transient.

Drawings

FIG. 1 is a basic flow diagram of the present invention;

FIG. 2 is a schematic view of an embodiment of a water conduit system;

FIG. 3 is a schematic diagram of an embodiment of meshing;

FIG. 4 is a graph showing the pressure change of an air mass during a water-air coupling transient process, according to an embodiment;

in the figure: 1-an upstream reservoir; 2-water conveying pipe; 3-a valve; 4-retention of gas mass; 5-impacting the water body; 6-tail water.

Detailed Description

The present invention will be further described with reference to the following examples.

As shown in fig. 1, a finite volume method-based simulation method of water-gas coupling transient flow in a pipeline is performed according to the following steps: 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; dividing a computational grid according to a finite volume method, and establishing a discrete equation; solving the numerical flux of a discrete equation of the water body by adopting a Godunov format, and obtaining second-order precision; 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; giving the CFL condition (Courant-Friedrichs-Lewycrion) satisfied by the Godunov format of the second order explicit finite volume method; and the Godunov format and an ideal gas state equation are adopted to realize the dynamic tracking of the water-gas interface simultaneously.

The technical solution of the present invention will be further described in detail with reference to the accompanying drawings and examples.

Example an experimental system as shown in figure 2 was used to study transient flow phenomena where water flow rapidly impacts a stagnant air mass. The whole system consists of an upstream pressure water tank, an organic transparent pipe and a valve. The total length of the pipeline is 8.83m, the inner diameter is 40mm, the pipeline is horizontal, the distance between the valve and the tail end of the pipeline is 3.25m, the upstream inlet pressure is 0.16MPa, and the tail end of the pipeline is closed. At the beginning, the valve is closed completely, and the water depth of tail water after the valve is 20 mm. The water hammer wave velocity is 400m/s, and the average friction resistance of the pipeline is 0.075-0.095. The water-gas coupling process is caused by suddenly opening a downstream valve, and the high-speed high-definition camera records that the time from full closing to full opening is 0.07 s-0.09 s.

The method comprises the following specific steps:

step 1: dividing transient flow in a pipeline system into three 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.

The basic differential equation of a water body is:

$\frac{\&part;H}{\&part;t}+V\frac{\&part;H}{\&part;x}+\frac{a^2}{g}\frac{\&part;V}{\&part;x}=0---\left(1\right)$

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

wherein, along the pipeline, the distance x and the time t are independent variables; h (x, t) is the piezometer tube head; v (x, t) is the average cross-sectional velocity; g is the acceleration of gravity; a is the wave velocity; f is the Darcy-Weisbaha friction coefficient; d is the pipe diameter.

The control equation for trapped air mass is:

$H_a{V}_a^m=H_{a0}{V}_{a0}^m{\mathrm{orH}}_a{L}_a^m=H_{a0}{L}_{a0}^m---\left(3\right)$

wherein H_a、V_a、L_aAir mass transient pressure, volume and length, respectively; their corresponding initial values are respectively H_a0、V_a0、L_a0(ii) a m is the polytropic index of the ideal gas equation of state, and since the transient process in the experiment is very short and is only a few seconds, 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-gas interface is:

$\frac{{\mathrm{dL}}_f}{dt}=V_w---\left(4\right)$

H_w＝H_a(5)

wherein L is_fIs the length of the impinging water body; v_w、H_wFlow rate and pressure at the water-gas interface, respectively.

In this embodiment, the calculation domain is a pipeline from the water outlet of the upstream reservoir to the valve; when the initial condition is that 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 water head of the upstream reservoir, the initial pressure of the retained air mass downstream of the valve is atmospheric pressure, and the gravity action is considered in the tail water pressure; the boundary conditions are as follows: at the inlet of the pipeline, the reservoir provides a constant pressure boundary which is equal to the static pressure water head of the upstream reservoir; the downstream ball valve is rapidly opened to establish hydraulic transient and the end of the pipeline is closed.

Step 2: and dividing the computational grid according to a finite volume method, and establishing a discrete equation.

Fig. 3 shows the discrete grid system of the launching hammer zone in the present embodiment. For the ith control body, the serial numbers of the upstream and downstream interfaces are defined as i-1/2 and i +1/2 respectively.

The specific steps for establishing the discrete equation are as follows:

(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{\&part;u}{\&part;t}+\frac{\&part;f\left(u\right)}{\&part;x}=s\left(u\right),f\left(u\right)=\stackrel{\&OverBar;}{A}u---\left(6\right)$

wherein $\begin{array}{ccc}u=\left(\begin{array}{c}H\\ V\end{array}\right);& s=\left(\begin{array}{c}0\\ -\frac{f\left|V\right|V}{2D}\end{array}\right);& \stackrel{\&OverBar;}{A}=\left(\begin{array}{cc}\stackrel{\&OverBar;}{V}& a^2/g\\ g& \stackrel{\&OverBar;}{V}\end{array}\right);\end{array}$ Is the average value of V, which is a constant.

(b) Integrating equation (6) over control unit i (from interface i-1/2 to i +1/2) and over time period Δ t (from t to t + Δ t) yields a discrete equation for flow variable u:

$U_i^{n+1}=U_i^n+\frac{\mathrm{\&Delta;}t}{\mathrm{\&Delta;}x}(f_{i+1/2}-f_{i-1/2})+\frac{\mathrm{\&Delta;}t}{\mathrm{\&Delta;}x}{\&Integral;}_{i-1/2}^{i+1/2}sdx---\left(7\right)$

wherein the superscripts n and n +1 represent t and t + Δ t time steps, respectively.Is the average value of u over the entire control volume; f is the flux at the cell interface;is a source item;

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

Further, step 3 comprises the following substeps:

step 3.1: solving flux at control unit interface inside water body

First, according to the Godunov method, its riemann problem is the following initial value problem:

$\frac{\&part;u}{\&part;t}+\frac{\&part;f\left(u\right)}{\&part;x}=0---\left(8\right)$

$u^n\left(x\right)=\left\{\begin{array}{cc}{U}_L^n,& x<x_{i+1/2}\\ U_R^n,& x>x_{i+1/2}\end{array}\right.---\left(9\right)$

wherein,is the average value of u to the left and right sides of the interface i +1/2 at n time steps.

For matrixAnd simplifying calculation to obtain characteristic values and characteristic vectors which are respectively as follows:

$\begin{array}{cc}\stackrel{\&OverBar;}{{\mathrm{\&lambda;}}_1}=\stackrel{\&OverBar;}{V}-a,& \stackrel{\&OverBar;}{{\mathrm{\&lambda;}}_2}=\stackrel{\&OverBar;}{V}+a\end{array}---\left(10\right)$

$K^{\left(1\right)}=\left[\begin{array}{c}1\\ -g/a\end{array}\right]and{K}^{\left(2\right)}=\left[\begin{array}{c}1\\ g/a\end{array}\right]---\left(11\right)$

because the characteristic vectors are linearly independent, further the following can be obtained:

$U_L=\underset{i=1}{\overset{2}{\&Sigma;}}{\mathrm{\&alpha;}}_i{K}^{\left(i\right)}and{U}_R=\underset{i=1}{\overset{2}{\&Sigma;}}{\mathrm{\&beta;}}_i{K}^{\left(i\right)}---\left(12\right)$

solving for four unknown coefficients α₁，α₂，β₁，β₂The following can be obtained:

$\begin{array}{cc}{\mathrm{\&alpha;}}_1=\frac{1}{2}(H_L^n-\frac{a}{g}{V}_L^n),& {\mathrm{\&alpha;}}_2=\frac{1}{2}(H_L^n+\frac{a}{g}{V}_L^n)\end{array}---\left(13\right)$

$\begin{array}{cc}{\mathrm{\&beta;}}_1=\frac{1}{2}(H_R^n-\frac{a}{9}{V}_R^n),& {\mathrm{\&beta;}}_2=\frac{1}{2}(H_R^n+\frac{a}{9}{V}_R^n)\end{array}---\left(14\right)$

the original variable form of the general solution to the riemann problem (equations (8) and (9)) is:

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

now, using equation (15), the exact solution of the variable at interface i +1/2 can be obtained:

$u_{i+1/2}\left(t\right)=\frac{1}{2}\{\left(\begin{array}{cc}1& a/g\\ g/a& 1\end{array}\right)U_L^n-\left(\begin{array}{cc}-1& a/g\\ g/a& -1\end{array}\right)U_R^n\}---\left(16\right)$

thus, for t ∈ [ tⁿ,tⁿ⁺¹]Any internal control unit i (1)<i<N), the flux at interface i +1/2 is:

$f_{i+1/2}={\stackrel{\&OverBar;}{A}}_{i+1/2}{u}_{i+1/2}\left(t\right)=\frac{1}{2}{\stackrel{\&OverBar;}{A}}_{i+1/2}\{\left(\begin{array}{cc}1& a/g\\ g/a& 1\end{array}\right)U_L^n-\left(\begin{array}{cc}-1& a/g\\ g/a& -1\end{array}\right)U_R^n\}---\left(17\right)$

next, MUSCL-Hancock format was introduced to calculate the internal element flux f_i+1/2And second-order precision is obtained. The method comprises the following specific steps:

(a) and (3) data reconstruction: using each control unit [ x ]_i-1/2,x_i+1/2]Inner piecewise linear function instead of data cell meanThen at the extreme pointThe values of (A) are:

$\begin{array}{cc}{U}_i^L=U_i^n-\frac{\mathrm{\&Delta;}x}{2}{\mathrm{\&Delta;}}_i,& U_i^R=U_i^n+\frac{\mathrm{\&Delta;}x}{2}{\mathrm{\&Delta;}}_i\end{array}---\left(18\right)$

wherein, Delta_iIs a moderate slope vector selected to increase the accuracy of the calculation format and to ensure that no spurious oscillations occur in the solution. In the present invention,. DELTA._iSelecting a MINMOD limiter:

${\mathrm{\&Delta;}}_i=MINMOD\left({\mathrm{\&sigma;}}_i^n,{\mathrm{\&sigma;}}_{i-1}^n\right)=\left\{\begin{array}{cc}{\mathrm{\&sigma;}}_i^n,& if,\left|{\mathrm{\&sigma;}}_i^n\right|<\left|{\mathrm{\&sigma;}}_{i-1}^n\right|,and,{\mathrm{\&sigma;}}_i^n{\mathrm{\&sigma;}}_{i-1}^n>0\\ {\mathrm{\&sigma;}}_{i-1}^n,& if,\left|{\mathrm{\&sigma;}}_i^n\right|>\left|{\mathrm{\&sigma;}}_{i-1}^n\right|,and,{\mathrm{\&sigma;}}_i^n{\mathrm{\&sigma;}}_{i-1}^n>0\\ 0,& if,{\mathrm{\&sigma;}}_i^n{\mathrm{\&sigma;}}_{i-1}^n<0\end{array}\right.---\left(19\right)$

wherein, $\begin{array}{cc}{\mathrm{\&sigma;}}_i^n=(U_{i+1}^n-U_i^n)/\mathrm{\&Delta;}x,& {\mathrm{\&sigma;}}_{i-1}^n=(U_i^n-U_{i-1}^n)/\mathrm{\&Delta;}x.\end{array}$

(b) and (4) calculating: for each unit [ x_i-1/2,x_i+1/2]Boundary extrapolation value in equation (18)Calculated by multiplying 0.5 Δ t by:

$\begin{array}{c}{\stackrel{\&OverBar;}{U}}_i^L=U_i^L+\frac{1}{2}\frac{\mathrm{\&Delta;}t}{\mathrm{\&Delta;}x}\&lsqb;f\left(U_i^L\right)-f\left(U_i^R\right)\&rsqb;,\\ {\stackrel{\&OverBar;}{U}}_i^R=U_i^R+\frac{1}{2}\frac{\mathrm{\&Delta;}t}{\mathrm{\&Delta;}x}\&lsqb;f\left(U_i^L\right)-f\left(U_i^R\right)\&rsqb;\end{array}---\left(20\right)$

(c) riemann problem: for calculating internal interface flux f_i+1/2The following data are combined to solve the conventional riemann problem:

$U_L^n\&equiv;{\stackrel{\&OverBar;}{U}}_i^R,U_R^n\&equiv;{\stackrel{\&OverBar;}{U}}_i^L---\left(21\right)$

by substituting equation (21) into equation (17), the flux at time t ═ t when the tube is filled with water can be obtainedⁿ,tⁿ⁺¹]Second order format at all internal cell interfaces i + 1/2.

Step 3.2: and constructing a virtual control unit to solve the flux at the interface of the control unit of the upstream and downstream boundaries 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 assume that the flow information at the virtual cell is consistent with that at the boundary. The boundary Riemann problem can thus be solved, and the corresponding Godunov flux f_1/2. And f_N+1/2The calculation can also be done like an 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.

In this embodiment, the constant water level boundary of the upstream reservoir:

at the upstream boundary, the riemann invariants associated with the negative feature line are:

$H_{1/2}-\frac{a}{g}{V}_{1/2}=H_1^n-\frac{a}{g}{V}_1^n---\left(22\right)$

wherein H_1/2＝H_res，H_resIs the upstream reservoir static pressure water head.

Thereby can be pushedVariable u_1/2(t)＝(H_1/2,V_1/2). Virtual unit I adjacent to the inlet of the duct according to the assumptions at the virtual unit_-1And I₀The corresponding values of (a) are:

$U_{-1}^{n+1}=U_0^{n+1}=U_{1/2}---\left(23\right)$

and 4, step 4: the Godunov format of the second-order explicit finite volume method is obtained by solving a numerical integral term of a discrete equation of the water body based on the Runge-Kutta method.

The specific implementation process is as follows:

(a) pure convection:

${\stackrel{\&OverBar;}{U}}_i^{n+1}=U_i^n-\frac{\mathrm{\&Delta;}t}{\mathrm{\&Delta;}x}(f_{i+1/2}^n-f_{i-1/2}^n)---\left(24\right)$

(b) update with the source term multiplied by Δ t/2:

${\stackrel{\&OverBar;}{\stackrel{\&OverBar;}{U}}}_i^{n+1}={\stackrel{\&OverBar;}{U}}_i^{n+1}+\frac{\mathrm{\&Delta;}t}{2}s\left({\stackrel{\&OverBar;}{U}}_i^{n+1}\right)---\left(25\right)$

(c) update again with the source term multiplied by Δ t:

$U_i^{n+1}={\stackrel{\&OverBar;}{U}}_i^{n-1}+\mathrm{\&Delta;}ts\left({\stackrel{\&OverBar;}{\stackrel{\&OverBar;}{U}}}_i^{n+1}\right)---\left(26\right)$

and 5: the CFL condition (Courant-Friedrichs-Lewycriton) satisfied by the Godunov format of the second order explicit finite volume method is given.

Maximum time step △ under CFL conditionst_max,CFL：

$C_r=\frac{Max\left(\stackrel{\&OverBar;}{\mathrm{\&lambda;}}\right)\&CenterDot;\mathrm{\&Delta;}t}{\mathrm{\&Delta;}x}\&le;1\&DoubleRightArrow;{\mathrm{\&Delta;t}}_{\max ,CFL}=\frac{\mathrm{\&Delta;}x}{Max\left(\stackrel{\&OverBar;}{\mathrm{\&lambda;}}\right)}---\left(27\right)$

Wherein, C_rIs the Koran number.

Furthermore, the introduced source terms satisfy the following stability constraints:

$\left|\frac{{\stackrel{\&OverBar;}{\stackrel{\&OverBar;}{U}}}_i^{n+1}}{{\stackrel{\&OverBar;}{U}}_i^{n+1}}\right|\&le;1---\left(28\right)$

maximum time step △ t for source item_max,s：

${\mathrm{\&Delta;t}}_{\max ,s}=min(-4\frac{{\stackrel{\&OverBar;}{U}}_i^{n+1}}{s\left({\stackrel{\&OverBar;}{U}}_i^{n+1}\right)},-2\frac{U_i^{n+1}}{s\left({\stackrel{\&OverBar;}{\stackrel{\&OverBar;}{U}}}_i^{n+1}\right)})---\left(29\right)$

The maximum allowable time step containing the convection item and the source item can be derived as:

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

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.

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 stagnant air mass are assumed to follow the change rule of the ideal gas state equation. And finally, combining the positive characteristic line with the Riemann vector based on the given conditions, and then combining the momentum and continuity equation of the water-gas interface and the control equation of the 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.

The method comprises the following specific steps:

at the water-gas interface, the Riemann invariant associated with the negative characteristic line is

$H_{A-W}+\frac{a}{g}{V}_{A-W}=H_{Nf+1}^n+\frac{a}{g}{V}_{Nf+1}^n---\left(31\right)$

Wherein,for impacting water body terminal control unit N_fMedium pressure and flow rate; h_A-W、V_A-WFor impacting pressure and flow velocity at the end of the body of water

In the water flow impact process, a continuity equation and a momentum equation at a water-gas interface are as follows:

A·V_A-W＝(A-A_C)·V_w(32)

g·(A_C·H_C-A·H_A-W)+A·V_A-W·(V_w-V_C)＝0(33)

wherein A is the sectional area of the pipeline; a. the_CIs the cross-sectional area of the wet periphery of tail water

The water body length change is as follows:

$L_f=L_f^n+\frac{V_w+V_w^n}{2}\mathrm{\&Delta;}t---\left(34\right)$

simultaneous air mass control equation, we can get:

$H_a=H_{a0}\&CenterDot;{\left(\frac{L_{a0}}{L-L_f^n-\frac{V_w+V_w^n}{2}\mathrm{\&Delta;}t}\right)}^m---\left(35\right)$

the 5 equations are solved simultaneously to obtain H_A-W,V_A-W,H_A,H_CAnd V_w。

Proximity I according to assumptions at the virtual cell_NVirtual unit I of_N+1And I_N+2The corresponding values of (a) are:

$U_{Nf+1}^{n+1}=U_{Nf+2}^{n+1}=\left(\begin{array}{c}{H}_{A-W}\\ V_{A-W}\end{array}\right)---\left(36\right)$

before the next time step of calculation, checking and reconstructing the calculation grids of the water body impacting part to ensure that the tail end micro-section water body meets 0<ΔL_fΔ x, wherein Δ L_f＝L_f-N_fΔ x. When the trapped air mass is compressed or expanded, Δ x may occur<ΔL_fOr Δ L_fUnder the condition of less than or equal to 0. To satisfy 0<ΔL_fAnd the requirement of less than or equal to deltax requires increasing or decreasing the calculation grids.

After the calculation is programmed by the method, the calculation results are compared with the experimental data. Fig. 4 shows a pressure curve of the air mass in the process of impacting the stagnant air mass by the water flow in the embodiment, and also shows a calculation result of neglecting the influence of the tail water (converting the equal volume of the stagnant air mass into the case of full-tube stagnant air mass, and neglecting the interaction between the tail water and the impacting water body). It can be seen that both the amplitude and the time response of the pressure calculated by the method of the present invention are in good agreement with the experiment. Meanwhile, it can be seen that neglecting the influence of the tail water causes large calculation errors in both the pressure amplitude and the period.

The above description is only of the preferred embodiments of the present invention, and it should be noted that: it will be apparent to those skilled in the art that various modifications and adaptations can be made without departing from the principles of the invention and these are intended to be within the scope of the invention.

Claims (9)

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{\mathrm{\&Delta;}t}{\mathrm{\&Delta;}x}(f_{i+1/2}-f_{i-1/2})+\frac{\mathrm{\&Delta;}t}{\mathrm{\&Delta;}x}{\&Integral;}_{i-1/2}^{i+1/2}sdx---\left(1\right)$ 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=\left(\begin{array}{c}H\\ V\end{array}\right),$ h is the piezometer tube head, V is the average section velocity; f is the flux at the cell interface; $s=\left(\begin{array}{c}0\\ -\frac{f\left|V\right|V}{2D}\end{array}\right),$ 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}={\stackrel{\&OverBar;}{A}}_{i+1/2}{u}_{i+1/2}\left(t\right)=\frac{1}{2}{\stackrel{\&OverBar;}{A}}_{i+1/2}\{\left(\begin{array}{cc}1& a/g\\ g/a& 1\end{array}\right)U_L^n-\left(\begin{array}{cc}-1& a/g\\ g/a& -1\end{array}\right)U_R^n\}---\left(2\right)$

wherein, $u=\left(\begin{array}{c}H\\ V\end{array}\right);\stackrel{\&OverBar;}{A}=\left(\begin{array}{cc}\stackrel{\&OverBar;}{V}& a^2/g\\ g& \stackrel{\&OverBar;}{V}\end{array}\right);$ 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}={\stackrel{\&OverBar;}{U}}_i^{n+1}+\mathrm{\&Delta;}ts\left({\stackrel{\&OverBar;}{\stackrel{\&OverBar;}{U}}}_i^{n+1}\right)---\left(3\right)$

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{\mathrm{Max}\left(\stackrel{\&OverBar;}{\mathrm{\&lambda;}}\right)\&CenterDot;\mathrm{\&Delta;t}}{\mathrm{\&Delta;x}}\&le;1\&DoubleRightArrow;{\mathrm{\&Delta;t}}_{\max ,\mathrm{CFL}}=\frac{\mathrm{\&Delta;x}}{\mathrm{Max}\left(\stackrel{\&OverBar;}{\mathrm{\&lambda;}}\right)}---\left(4\right).$

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：

$\left|\frac{{\stackrel{\&OverBar;}{\stackrel{\&OverBar;}{U}}}_i^{n+1}}{{\stackrel{\&OverBar;}{U}}_i^{n+1}}\right|\&le;1\&DoubleRightArrow;{\mathrm{\&Delta;t}}_{\max ,s}=\min (-4\frac{{\stackrel{\&OverBar;}{U}}_i^{n+1}}{s\left({\stackrel{\&OverBar;}{U}}_i^{n+1}\right)},-2\frac{U_i^{n+1}}{s\left({\stackrel{\&OverBar;}{\stackrel{\&OverBar;}{U}}}_i^{n+1}\right)})---\left(5\right).$

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.