CN113656926B - Pipeline transient flow simulation method based on Schohl convolution approximation - Google Patents
Pipeline transient flow simulation method based on Schohl convolution approximation Download PDFInfo
- Publication number
- CN113656926B CN113656926B CN202110986943.6A CN202110986943A CN113656926B CN 113656926 B CN113656926 B CN 113656926B CN 202110986943 A CN202110986943 A CN 202110986943A CN 113656926 B CN113656926 B CN 113656926B
- Authority
- CN
- China
- Prior art keywords
- term
- pipeline
- schohl
- boundary
- calculation
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 84
- 230000001052 transient effect Effects 0.000 title claims abstract description 22
- 238000004088 simulation Methods 0.000 title claims abstract description 20
- 238000004364 calculation method Methods 0.000 claims abstract description 44
- 230000004907 flux Effects 0.000 claims abstract description 26
- 239000011159 matrix material Substances 0.000 claims abstract description 19
- 239000013598 vector Substances 0.000 claims description 27
- 230000008569 process Effects 0.000 claims description 18
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 claims description 16
- 238000002474 experimental method Methods 0.000 claims description 4
- 238000011144 upstream manufacturing Methods 0.000 claims description 4
- 230000001133 acceleration Effects 0.000 claims description 3
- 239000012530 fluid Substances 0.000 claims description 3
- 230000005484 gravity Effects 0.000 claims description 3
- 238000012886 linear function Methods 0.000 claims description 3
- 230000010354 integration Effects 0.000 claims 1
- 230000008859 change Effects 0.000 description 2
- 238000012544 monitoring process Methods 0.000 description 2
- 230000010355 oscillation Effects 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000001514 detection method Methods 0.000 description 1
- 239000003651 drinking water Substances 0.000 description 1
- 235000020188 drinking water Nutrition 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 239000007788 liquid Substances 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/10—Geometric CAD
- G06F30/18—Network design, e.g. design based on topological or interconnect aspects of utility systems, piping, heating ventilation air conditioning [HVAC] or cabling
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
- G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
- G06F30/28—Design optimisation, verification or simulation using fluid dynamics, e.g. using Navier-Stokes equations or computational fluid dynamics [CFD]
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2111/00—Details relating to CAD techniques
- G06F2111/04—Constraint-based CAD
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2111/00—Details relating to CAD techniques
- G06F2111/10—Numerical modelling
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2113/00—Details relating to the application field
- G06F2113/08—Fluids
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2113/00—Details relating to the application field
- G06F2113/14—Pipes
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
- G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Mathematical Physics (AREA)
- Geometry (AREA)
- Mathematical Analysis (AREA)
- Mathematical Optimization (AREA)
- Pure & Applied Mathematics (AREA)
- General Engineering & Computer Science (AREA)
- Evolutionary Computation (AREA)
- Computer Hardware Design (AREA)
- Computational Mathematics (AREA)
- Algebra (AREA)
- Data Mining & Analysis (AREA)
- Computing Systems (AREA)
- Fluid Mechanics (AREA)
- Computer Networks & Wireless Communication (AREA)
- Operations Research (AREA)
- Databases & Information Systems (AREA)
- Software Systems (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
The invention discloses a pipeline transient flow simulation method based on Schohl convolution approximation, which comprises the following steps: constructing a pipeline system control equation containing Zielke non-constant friction terms and viscoelasticity terms; expressing a control equation in a matrix form in a solution format of a Riemann problem, and putting a non-constant friction term and a viscoelastic term into a source term; establishing a calculation grid under a finite volume method solving system, and calculating grid boundary flux through a second-order Godunov format; performing Schohl convolution approximation on a Zielke non-constant friction term and a viscoelasticity term in a source term, merging the source term into a matrix equation, and performing stability constraint; and (5) giving out calculation initial conditions by combining with actual engineering to obtain a calculation result. The invention combines the Godunov format of the finite volume method to carry out Schohl convolution approximation on the non-constant friction term and the viscoelasticity term in the control equation, and can solve the problems that the existing one-dimensional transient flow simulation method depends on the Brownet number, and involves long calculation time, complex solution and the like of the non-constant friction term and the viscoelasticity term.
Description
Technical Field
The invention belongs to the field of computational hydraulics, and particularly relates to a transient flow non-constant friction resistance and viscoelasticity simulation method for a pipeline based on Schohl convolution approximation.
Background
The water pipe network system can generate instant high-pressure flow, called instant flow, due to the sudden opening and closing of the valve. If the pressure is too high, the pipeline system can be leaked or damaged, even a pipe burst is caused, and the safe supply of the drinking water is seriously threatened. Therefore, the detection of the liquid flow state of the pipeline system is important to ensuring the safe and stable operation of the water delivery pipeline network.
One-dimensional transient flow simulation of a piping system is an effective means of ascertaining the flow state of the piping system. In the past one-dimensional transient flow model, a Zielke original weighting function non-constant friction model is mostly used, and a generalized K-V viscoelastic model is combined, so that the pressure and flow change of the transient pipeline are solved by a characteristic line method. The method has the following defects: 1) The Zielke original weighting function model needs to store the flow velocity change at the historical moment, so that the calculation time is longer; 2) The original K-V viscoelastic model solving method is complex and uses first order approximation; 3) The characteristic line method is convenient for a single tube, but interpolation approximation is needed in complex pipeline simulation, and numerical dissipation is easy to generate.
Therefore, an accurate, efficient and simple transient flow simulation method for the water delivery pipe network is needed, and the method has important theoretical significance and practical application value for safe operation of the water delivery pipe system and real-time monitoring of the pipe state.
Disclosure of Invention
The invention aims to: in order to solve the problems of long calculation time, complex solution, poor universality, difficulty in real-time monitoring and the like in a one-dimensional transient flow model in the prior art, the method for simulating the transient flow non-constant friction resistance and the viscoelasticity of the pipeline based on Schohl convolution approximation is provided.
The technical scheme is as follows: to achieve the above object, the present invention provides a method for simulating transient flow non-constant friction and viscoelasticity of a pipeline based on Schohl convolution approximation, comprising the steps of:
s1: constructing a pipeline system control equation containing Zielke non-constant friction terms and viscoelasticity terms;
s2: expressing a control equation in a matrix form in a solution format of a Riemann problem, and putting a non-constant friction term and a viscoelastic term into a source term;
s3: establishing a calculation grid under a finite volume method solving system, and calculating grid boundary flux through a second-order Godunov format;
s4: performing Schohl convolution approximation on a Zielke non-constant friction term and a viscoelasticity term in a source term, merging the source term into a matrix equation, and performing stability constraint;
s5: and (3) giving calculation initial conditions by combining with actual engineering, and obtaining a calculation result by utilizing Fortran programming.
Further, the equation of the control equation of the pipeline system in the step S1 is:
wherein H is the pressure at a certain section of the pipeline; v is the average flow velocity of a certain section of the pipeline; epsilon r Is the hysteresis strain of the pipeline; a is the wave velocity in the water of the pressurized pipeline; g is gravity acceleration; x is the distance along the axis of the pipe; t is time; d is the inner diameter of the pipeline; ρ is the fluid density; f is the darcy-visbach coefficient; v is the dynamic viscosity of water; u is a differential variable; w is a weighted function of historical flow rate variation.
Further, the step S2 specifically includes:
wherein,
where u is the solution vector, f (u) is the control volume flux vector,for the partial derivative matrix of each solving variable in each variable pair u in the vector f, s (u) is a source term and represents the viscoelasticity term and the friction term (explicit) of the pipeline respectively, and if s (u) =0, i.e. the viscoelasticity and the friction term of the pipeline are not considered, the system is a constant coefficient linear homogeneous system.
Further, the method for establishing the computing grid in the step S3 is as follows:
under a finite volume method solving system, dividing a pipe section into N control bodies, establishing a calculation grid, wherein the length delta x of each control body is calculated, the calculation time step length is delta t, and for the i (i is more than or equal to 1 and less than or equal to N) control bodies, the left and right boundaries of the i control bodies are respectively marked as i-1/2 and i+1/2; to solve the Riemann problem, the matrix equation is integrated along the x-direction between the two boundary interfaces of the ith control body to obtain the following formula:
wherein f i-1/2 A numerical flux representing the i-1/2 interface; f (f) i+1/2 A numerical flux representing the i+1/2 interface; vector u may be derived from the average value within control volume iThe above formula can be organized as:
in the formula, the superscript n represents the time t, namely the current calculation time step; the superscript n+1 indicates the time t+Δt, i.e., the next calculation time step. Thus, when calculating the vector U of the next time step, the two boundary value fluxes and the source term of the control body of the current time step are needed to be related.
Further, the calculation method of the grid boundary flux in the step S3 is as follows:
for the finite volume method Godunov solution format, the numerical flux of the control volume boundary can be derived from the local Riemann problem at the interface: for a constant coefficient linear homogeneous hyperbolic system, the Riemann problem can be described as the following initial problem:
in the method, in the process of the invention,average value of vector u on the left side of i+1/2 interface in the nth calculation step,/>And->Is vector->Vector elements of (a); />Average value of vector u on the right side of i+1/2 interface in the nth calculation step,/>Andis vector->Vector elements of (a); x is x i+1/2 The distance along the axis of the pipeline is the right side edge of the ith control body;
by solving characteristic equationsTwo mutually independent characteristic values can be obtained:
wherein I is an identity matrix,j=1 to 2 for the eigenvalue to be solved; according to the Rankine-Hugonitot condition At t E [ t, t+Deltat]Flux values at the inner boundary i+1/2 +.>Can be expressed as:
wherein,
for the followingAnd->The approximation is reconstructed using a piecewise polynomial, the accuracy of which depends on the form of the polynomial.
Further, in the step S3, forAnd->Linear interpolation by introducing MUSCL-Hancock format using second order Godunov approximationThe reconstruction comprises the following specific steps:
a1: and (3) data reconstruction:
by constructing piecewise linear function U i (x) To obtain the boundary extremum of the left and right of the ith control bodyAnd->
Wherein Ω i Is thatIn the slope vector of the ith calculation unit, properly selecting Ω i The values of (2) can improve the accuracy of the scheme while ensuring that the result is free of spurious oscillations. Here, a MINMOD (flow rate limiting function) slope limiter is selected:
in the method, in the process of the invention,
a2: time evolution calculation:
for each control volume [ x ] i-1/2 ，x i+1/2 ]Left and right boundary extremumEvolution at 0.5 Δt time gave:
in the method, in the process of the invention,left and right boundary extremum after time advance, < ->And->Left and right boundary extremum->The flux at the point(s) is (are),
a3: approximation solution of Riemann problem:
for boundary conditions, two virtual control bodies-1, 0 and N+1, N+2 are respectively added at the boundaries of two ends, so that unified calculation simulation of the boundary control body and the internal control body can be obtained, and the second-order precision requirement is met;
from the analytical solution of the Riemann problem at the boundary:
the Riemann invariant equation associated with positive and negative feature lines can be deduced:
an upstream water reservoir boundary is provided with
And is also provided withSubstituting into the above to obtain
The virtual units thereof satisfy:
downstream boundary conditions, with
Wherein C is V ＝(Q 0 τ) 2 /2ΔH 0 ，Q 0 And DeltaH 0 Respectively the flow and head loss at the valve at the initial moment, and tau is the valve opening; h RD Is a downstream reservoir head; when C P -H RD At > 0, s= +1, when C P -H RD When < 0, s= -1;
for downstream valve + reservoir boundaries
From the end valve orifice equation, it can be seen that
Is brought into a valve positive characteristic line Riemann invariant equation to obtain
The virtual units thereof satisfy:
further, the specific method for performing Schohl convolution approximation on the Zielke non-constant friction term and the viscoelastic term in the source term in the step S4 is as follows:
for Zielke non-constant friction term, in equation (2), let
In the method, in the process of the invention,
wherein W is app (τ) is an approximation of W (τ), by W i A sum function composition of (τ), i=1 to 5; m is m i And n i The coefficients of the weighting functions, respectively, approximate friction, m for Schohl convolution i ＝(1.051；2.358；9.021；29.47；79.55)；n i = (26.65; 100;669.6;6497; 57990). Order the
There is a case where the number of the group,
therefore, it is
For the viscoelastic term, there are
Wherein ε r (t) is ε r Is a function of time of (2); alpha is poisson's ratio; e is the wall thickness of the pipeline; j (t) is creep compliance, t is differential variable; for the generalized K-V model, the creep compliance function can be expressed as
Wherein J is k And τ k Are all viscoelastic mechanical model parameters and can be obtained by checking according to creep experiment data; m is J obtained by checking k Is the number of (3); obtaining partial derivative of J (t)
According to the differential nature of convolution, there are
Order the
I.e.
Is available from the interchangeability of convolution
Order the
Then
According to the Schohl convolution approximation, for the t+Δt time instant, there is
Therefore, it is
In the method, in the process of the invention,
CO A3 ＝CO A1 ·CO A2 (45)
CO A4 ＝CO A1 (1-CO A2 ) (46)。
further, in the step S4, the specific method for merging the source term into the matrix equation and performing stability constraint is as follows:
in the step S3, the influence of the source item is not considered, the source item can be introduced into a nonlinear system to be solved by adopting a time division method, and the solving precision depends on an integral mode of the source item; for second-order precision, solving a source term s (u) and spatial integral thereof by adopting an explicit second-order Runge-Kutta discretization method, and obtaining the source term by updating the source term twice:
in the method, in the process of the invention,
stability constraints mainly include the CFL (Courant-Friedrich-Lewy) condition of the stream item and the constraint when the source item is updated twice; for CFL conditions, there are
Wherein C is r For the Brownian number, Δt max，CFL Maximum time step under CFL conditions;
the source term is introduced by using the second-order Runge-Kutta format, and the constraint condition is that
Thus can obtain
Therefore, the maximum allowable time step is calculated to be satisfied
Δt≤min(Δt max，CFL ，Δt max，s ) (54)
The beneficial effects are that: compared with the prior art, the invention has the following advantages:
1. schohl convolution approximation is carried out on the Zielke original weighting function model, so that the calculation time is reduced.
2. The model provided by the invention provides a novel generalized K-V model solving method based on Schohl convolution approximation, so that the solving difficulty is simplified, and the calculation accuracy is improved.
3. And a second-order Godunov format is used for solving a control equation, so that the numerical dissipation degree in complex pipeline simulation is reduced.
Drawings
FIG. 1 is a flow chart of an implementation of the present invention;
FIG. 2 is a schematic illustration of a computational grid along the axial direction of a pipeline;
FIG. 3 is a graph showing the simulation results of the inverse method of the present invention compared with the experiment.
Detailed Description
The present invention is further illustrated in the accompanying drawings and detailed description which are to be understood as being merely illustrative of the invention and not limiting of its scope, and various modifications of the invention, which are equivalent to those skilled in the art upon reading the invention, will fall within the scope of the invention as defined in the appended claims.
As shown in FIG. 1, the invention provides a method for simulating transient flow non-constant friction and viscoelasticity of a pipeline based on Schohl convolution approximation, which comprises the following steps:
s1: constructing a pipeline system control equation containing Zielke non-constant friction terms and viscoelasticity terms:
wherein H is the pressure at a certain section of the pipeline; v is the average flow velocity of a certain section of the pipeline; epsilon r Is the hysteresis strain of the pipeline; a is the wave velocity in the water of the pressurized pipeline; g is gravity acceleration; x is the distance along the axis of the pipe; t is time; d is the inner diameter of the pipeline; ρ is the fluid density; f is the darcy-visbach coefficient; v is the dynamic viscosity of water; u is a differential variable; w is a weighted function of historical flow rate variation.
S2: the control equation is expressed in a matrix form in a solution format of the Riemann problem, and a non-constant friction term and a viscoelastic term are put into a source term:
wherein,
where u is the solution vector, f (u) is the control volume flux vector,for the partial derivative matrix of each solving variable in each variable pair u in the vector f, s (u) is a source term and represents the viscoelasticity term and the friction term (explicit) of the pipeline respectively, and if s (u) =0, i.e. the viscoelasticity and the friction term of the pipeline are not considered, the system is a constant coefficient linear homogeneous system.
S3: 1. establishing a calculation grid under a finite volume method solving system:
as shown in FIG. 2, under a finite volume method solving system, a pipe section is divided into N control bodies, a calculation grid is established, the length Deltax of each control body is calculated, the calculation time step length is Deltat, and for the ith (i is more than or equal to 1 and less than or equal to N) control bodies, the left and right boundaries of the ith control body are respectively marked as i-1/2 and i+1/2; to solve the Riemann problem, the matrix equation is integrated along the x-direction between the two boundary interfaces of the ith control body to obtain the following formula:
wherein f i-1/2 A numerical flux representing the i-1/2 interface; f (f) i+l/2 A numerical flux representing the i+1/2 interface; vector u may be derived from the average value within control volume iThe above formula can be organized as:
in the formula, the superscript n represents the time t, namely the current calculation time step; the superscript n+1 indicates the time t+Δt, i.e., the next calculation time step. Thus, when calculating the vector U of the next time step, the two boundary value fluxes and the source term of the control body of the current time step are needed to be related.
2. Grid boundary flux is calculated by a second order Godunov format:
for the finite volume method Godunov solution format, the numerical flux of the control volume boundary can be derived from the local Riemann problem at the interface: for a constant coefficient linear homogeneous hyperbolic system, the Riemann problem can be described as the following initial problem:
in the method, in the process of the invention,average value of vector u on the left side of i+1/2 interface in the nth calculation step,/>And->Is vector->Vector elements of (a); />Average value of vector u on the right side of i+1/2 interface in the nth calculation step,/>And->Is vector->Vector elements of (a); x is x i+1/2 The distance along the axis of the pipeline is the right side edge of the ith control body;
by solving characteristic equationsTwo mutually independent characteristic values can be obtained:
wherein I is an identity matrix,j=1 to 2 for the eigenvalue to be solved; according to the Rankine-Hugonitot condition At t E [ t ]，t+Δt]Flux values at the inner boundary i+1/2 +.>Can be expressed as:
wherein,
for the followingAnd->The approximation is reconstructed using a piecewise polynomial, the accuracy of which depends on the form of the polynomial.
For the followingAnd->The linear interpolation reconstruction is carried out by adopting second-order Godunov approximation and introducing MUSCL-Hancock format, and the specific steps are as follows:
a1: and (3) data reconstruction:
by constructing piecewise linear function u i (x) To obtain the boundary extremum of the left and right of the ith control bodyAnd->
Wherein Ω i Is thatIn the first place i Slope vectors of the calculation units are properly selected to be omega i The values of (2) can improve the accuracy of the scheme while ensuring that the result is free of spurious oscillations. Here, a MINMOD (flow rate limiting function) slope limiter is selected:
in the method, in the process of the invention,
a2: time evolution calculation:
for each control volume [ x ] i-1/2 ，x i+1/2 ]Left and right boundary extremumEvolution at 0.5 Δt time gave:
in the method, in the process of the invention,left and right boundary extremum after time advance, < ->And->Left and right boundary extremum->The flux at the point(s) is (are),
a3: approximation solution of Riemann problem:
for boundary conditions, two virtual control bodies-1, 0 and N+1, N+2 are respectively added at the boundaries of two ends, so that unified calculation simulation of the boundary control body and the internal control body can be obtained, and the second-order precision requirement is met;
from the analytical solution of the Riemann problem at the boundary:
the Riemann invariant equation associated with positive and negative feature lines can be deduced:
an upstream water reservoir boundary is provided with
And is also provided withIs substituted intoObtaining
The virtual units thereof satisfy:
downstream boundary conditions, with
Wherein C is V ＝(Q 0 τ) 2 /2ΔH 0 ，Q 0 And DeltaH 0 Respectively the flow and head loss at the valve at the initial moment, and tau is the valve opening; h RD Is a downstream reservoir head; when C P -H RD >When 0, S= +1, when C P -H RD <At 0, s= -1;
for downstream valve + reservoir boundaries
From the end valve orifice equation, it can be seen that
Is brought into a valve positive characteristic line Riemann invariant equation to obtain
The virtual units thereof satisfy:
s4: 1. schohl convolution approximation was performed on the Zielke non-constant friction term and the viscoelastic term in the source term:
for Zielke non-constant friction term, in equation (2), let
In the method, in the process of the invention,
wherein W is app (τ) is an approximation of W (τ), by W i A sum function composition of (τ), i=1 to 5; m is m i And n i The coefficients of the weighting functions, respectively, approximate friction, m for Schohl convolution i ＝(1.051；2.358；9.021；29.47；79.55)；n i = (26.65; 100;669.6;6497; 57990). Order the
There is a case where the number of the group,
therefore, it is
For the viscoelastic term, there are
Wherein ε r (t) is ε r Is a function of time of (2); alpha is poisson's ratio; e is the wall thickness of the pipeline; j (t) is creep compliance, t is differential variable; for generalized K-V model, vermicularThe variable compliance function may be expressed as
Wherein J is k And τ k Are all viscoelastic mechanical model parameters and can be obtained by checking according to creep experiment data; m is J obtained by checking k Is the number of (3); obtaining partial derivative of J (t)
According to the differential nature of convolution, there are
Order the
I.e.
Is available from the interchangeability of convolution
Order the
Then
According to the Schohl convolution approximation, for the t+Δt time instant, there is
Therefore, it is
In the method, in the process of the invention,
CO A3 ＝CO A1 ·CO A2 (45)
CO A4 ＝CO A1 (1-CO A2 ) (46)
2. incorporating the source term into the matrix equation and performing stability constraints:
in the step S3, the influence of the source item is not considered, the source item can be introduced into a nonlinear system to be solved by adopting a time division method, and the solving precision depends on an integral mode of the source item; for second-order precision, solving a source term s (u) and spatial integral thereof by adopting an explicit second-order Runge-Kutta discretization method, and obtaining the source term by updating the source term twice:
in the method, in the process of the invention,
/>
stability constraints mainly include the CFL (Courant-Friedrich-Lewy) condition of the stream item and the constraint when the source item is updated twice; for CFL conditions, there are
Wherein C is r For the Brownian number, Δt max，CFL Maximum time step under CFL conditions;
the source term is introduced by using the second-order Runge-Kutta format, and the constraint condition is that
Thus can obtain
Therefore, the maximum allowable time step is calculated to be satisfied
Δt≤min(Δt max，CFL ，Δt max，s ) (54)
S5: giving calculation initial conditions by combining actual engineering, and obtaining calculation results by utilizing Fortran programming:
the initial conditions brought into the actual working condition comprise an upstream reservoir or water tank water head Hr, a pipeline diameter D, the number N of control bodies divided according to the pipeline length, the length Deltax of each section of control body, a pipeline wave speed a, a Reynolds number Re and the like, and the pressure measuring pipe water head H and the flow velocity V of each control body in the pipeline can be obtained, so that the simulation of the pipeline hydraulic transient is realized.
Based on the above scheme, in this embodiment, in order to verify the effect of the method of the present invention, a viscoelastic pipeline hydraulic transient experimental apparatus system designed and built by Covas in 2003 is selected to verify the effectiveness of the method of the present invention, and experimental and simulation comparison results are shown in fig. 3, which shows that the new method of the present invention can well simulate the pressure peak value, attenuation and phase of the pipeline hydraulic transient.
Claims (6)
1. A method for transient flow simulation of a pipeline based on Schohl convolution approximation, comprising the steps of:
s1: constructing a pipeline system control equation containing Zielke non-constant friction terms and viscoelasticity terms;
s2: expressing a control equation in a matrix form in a solution format of a Riemann problem, and putting a non-constant friction term and a viscoelastic term into a source term;
s3: establishing a calculation grid under a finite volume method solving system, and calculating grid boundary flux through a second-order Godunov format;
s4: performing Schohl convolution approximation on a Zielke non-constant friction term and a viscoelasticity term in a source term, merging the source term into a matrix equation, and performing stability constraint;
s5: giving out calculation initial conditions by combining with actual engineering to obtain a calculation result;
the specific method for performing Schohl convolution approximation on the Zielke non-constant friction term and the viscoelasticity term in the source term in the step S4 is as follows:
for Zielke non-constant friction term, in equation (2), let
In the method, in the process of the invention,
wherein W is app (τ) is an approximation of W (τ), by W i A sum function of (τ); m is m i And n i The coefficients of the weighting functions, respectively, approximate friction for Schohl convolution, let
There is a case where the number of the group,
therefore, it is
For the viscoelastic term, there are
Wherein ε r (t) is ε r Is a function of time of (2); alpha is poisson's ratio; e is the wall thickness of the pipeline; j (t) is creep compliance, t is differential variable; for the generalized K-V model, the creep compliance function can be expressed as
Wherein J is k And τ k Are all viscoelastic mechanical model parameters and can be obtained by checking according to creep experiment data; m is J obtained by checking k Is the number of (3); obtaining partial derivative of J (t)
According to the differential nature of convolution, there are
Order the
I.e.
Is available from the interchangeability of convolution
Order the
Then
According to the Schohl convolution approximation, for the t+Δt time instant, there is
Therefore, it is
In the method, in the process of the invention,
CO A3 ＝CO A1 ·CO A2 (45)
CO A4 ＝CO A1 (1-CO A2 ) (46)
in the step S4, the specific method for merging the source term into the matrix equation and performing stability constraint is as follows:
introducing a source term into a nonlinear system to solve by adopting a time division method, wherein the solving precision depends on an integration mode of the source term; for second-order precision, solving a source term s (u) and spatial integral thereof by adopting an explicit second-order Runge-Kutta discretization method, and obtaining the source term by updating the source term twice:
in the method, in the process of the invention,
stability constraints include constraints on CFL conditions for a flow item and on two updates to a source item; for CFL conditions, there are
Wherein C is r For the Brownian number, Δt max,CFL Maximum time step under CFL conditions;
the source term is introduced by using the second-order Runge-Kutta format, and the constraint condition is that
Thus can obtain
Therefore, the maximum allowable time step is calculated to be satisfied
Δt≤min(Δt max,CFL ,Δt max,s ) (54)。
2. The method for transient flow simulation of a pipeline based on Schohl convolution approximation according to claim 1, wherein the equation of the control equation of the pipeline system in step S1 is:
wherein H is the pressure at a certain section of the pipeline; v is the average flow velocity of a certain section of the pipeline; epsilon r Is the hysteresis strain of the pipeline; a is the wave velocity in the water of the pressurized pipeline; g is gravity acceleration; x is the distance along the axis of the pipe; t isTime; d is the inner diameter of the pipeline; ρ is the fluid density; f is the darcy-visbach coefficient; v is the dynamic viscosity of water; u is a differential variable; w is a weighted function of historical flow rate variation.
3. The method for transient flow simulation of a pipeline based on Schohl convolution approximation according to claim 1, wherein the step S2 is specifically:
wherein,
where u is the solution vector, f (u) is the control volume flux vector,the partial derivative matrix of each solving variable in each variable pair u in the vector f is obtained, s (u) is a source term and represents the viscoelasticity term and the friction resistance term of the pipeline respectively, and if s (u) =0, i.e. the viscoelasticity and the friction resistance term of the pipeline are not considered, the system is a constant coefficient linear homogeneous system.
4. The method for simulating pipeline transient flow based on Schohl convolution approximation according to claim 1, wherein the method for establishing the computing grid in step S3 is as follows:
under a finite volume method solving system, dividing a pipe section into N control bodies, establishing a calculation grid, wherein the length delta x of each control body is calculated, the calculation time step length is delta t, and for the ith control body, i is more than or equal to 1 and less than or equal to N, the left and right boundaries of the control bodies are respectively marked as i-1/2 and i+1/2; to solve the Riemann problem, the matrix equation is integrated along the x-direction between the two boundary interfaces of the ith control body to obtain the following formula:
wherein f i-1/2 A numerical flux representing the i-1/2 interface; f (f) i+1/2 A numerical flux representing the i+1/2 interface; vector u may be derived from the average value within control volume iThe above formula can be organized as:
in the formula, the superscript n represents the time t, namely the current calculation time step; the superscript n+1 indicates the time t+Δt, i.e., the next calculation time step.
5. The method for transient flow simulation of a pipeline based on Schohl convolution approximation according to claim 4, wherein the calculation method of the grid boundary flux in step S3 is as follows:
for the finite volume method Godunov solution format, the numerical flux of the control volume boundary can be derived from the local Riemann problem at the interface: for a constant coefficient linear homogeneous hyperbolic system, the Riemann problem can be described as the following initial problem:
in the method, in the process of the invention, and->Is vector->Vector elements of (a); /> And->Is vector->Vector elements of (a); x is x i+1/2 The distance along the axis of the pipeline is the right side edge of the ith control body;
by solving characteristic equationsTwo mutually independent characteristic values can be obtained:
wherein I is an identity matrix,is a feature value of the solution; according to the Rankine-Hugonitot condition +.> At t E [ t, t+Deltat]Flux values at the inner boundary i+1/2 +.>Can be expressed as:
wherein,
for the followingAnd->The approximation is reconstructed using a piecewise polynomial.
6. The method for transient flow simulation of a pipeline based on Schohl convolution approximation according to claim 5, wherein in step S3, forAnd->The linear interpolation reconstruction is carried out by adopting second-order Godunov approximation and introducing MUSCL-Hancock format, and the specific steps are as follows:
a1: and (3) data reconstruction:
by constructing piecewise linear function U i (x) To obtain the boundary extremum of the left and right of the ith control bodyAnd->
Wherein Ω i Is thatAt the slope vector of the ith calculation unit, a flow rate limiting function MINMOD slope limiter is selected:
in the method, in the process of the invention,
a2: time evolution calculation:
for each control volume [ x ] i-1/2 ，x i+1/2 ]Left and right boundary extremumEvolution at 0.5 Δt time gave:
in the method, in the process of the invention,left and right boundary extremum after time advance, < ->And->Respectively, left and right boundary extreme valuesThe flux at the point(s) is (are),
a3: approximation solution of Riemann problem:
for boundary conditions, two virtual control bodies-1, 0 and N+1, N+2 are respectively added at the boundaries of two ends, so that unified calculation simulation of the boundary control body and the internal control body can be obtained, and the second-order precision requirement is met;
from the analytical solution of the Riemann problem at the boundary:
the Riemann invariant equation associated with positive and negative feature lines can be deduced:
an upstream water reservoir boundary is provided with
And is also provided withSubstituting into the above to obtain
The virtual units thereof satisfy:
downstream boundary conditions, with
Wherein C is V ＝(Q 0 τ) 2 /2ΔH 0 ，Q 0 And DeltaH 0 Respectively the flow and head loss at the valve at the initial moment, and tau is the valve opening; h RD Is a downstream reservoir head; when C P -H RD >When 0, S= +1, when C P -H RD <At 0, s= -1;
for downstream valve + reservoir boundaries
From the end valve orifice equation, it can be seen that
Is brought into a valve positive characteristic line Riemann invariant equation to obtain
The virtual units thereof satisfy:
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110986943.6A CN113656926B (en) | 2021-08-26 | 2021-08-26 | Pipeline transient flow simulation method based on Schohl convolution approximation |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110986943.6A CN113656926B (en) | 2021-08-26 | 2021-08-26 | Pipeline transient flow simulation method based on Schohl convolution approximation |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN113656926A CN113656926A (en) | 2021-11-16 |
| CN113656926B true CN113656926B (en) | 2024-03-26 |
Family
ID=78492897
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202110986943.6A Active CN113656926B (en) | 2021-08-26 | 2021-08-26 | Pipeline transient flow simulation method based on Schohl convolution approximation |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN113656926B (en) |
Families Citing this family (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111767684B (en) * | 2020-06-30 | 2024-04-19 | 西安理工大学 | Modeling method of optimized friction source item implicit format two-dimensional shallow water equation |
| CN114896908B (en) * | 2022-05-19 | 2023-03-28 | 四川大学 | Water drop flow field based on flux limiter and water drop collection rate calculation method |
| CN115391069B (en) * | 2022-10-27 | 2023-02-03 | 山东省计算中心（国家超级计算济南中心） | Parallel communication method and system based on ocean mode ROMS |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109918787A (en) * | 2019-03-08 | 2019-06-21 | 河海大学 | The analogy method of aqueous vapor two-phase homogeneous flow in aqueduct based on finite volume method |
| CN111339701A (en) * | 2020-02-26 | 2020-06-26 | 河海大学 | Pipeline leakage characteristic Godunov simulation method based on Brunone dynamic friction resistance |
| CN111695307A (en) * | 2020-05-21 | 2020-09-22 | 河海大学 | Water hammer finite volume simulation method considering dynamic friction resistance explicitly |
| CN111985166A (en) * | 2020-07-31 | 2020-11-24 | 河海大学 | Pipeline hydraulic transient simulation method and storage medium with implicit consideration of dynamic friction resistance |
-
2021
- 2021-08-26 CN CN202110986943.6A patent/CN113656926B/en active Active
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109918787A (en) * | 2019-03-08 | 2019-06-21 | 河海大学 | The analogy method of aqueous vapor two-phase homogeneous flow in aqueduct based on finite volume method |
| CN111339701A (en) * | 2020-02-26 | 2020-06-26 | 河海大学 | Pipeline leakage characteristic Godunov simulation method based on Brunone dynamic friction resistance |
| CN111695307A (en) * | 2020-05-21 | 2020-09-22 | 河海大学 | Water hammer finite volume simulation method considering dynamic friction resistance explicitly |
| CN111985166A (en) * | 2020-07-31 | 2020-11-24 | 河海大学 | Pipeline hydraulic transient simulation method and storage medium with implicit consideration of dynamic friction resistance |
Also Published As
| Publication number | Publication date |
|---|---|
| CN113656926A (en) | 2021-11-16 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN113656926B (en) | Pipeline transient flow simulation method based on Schohl convolution approximation | |
| Ramos et al. | Surge damping analysis in pipe systems: Modelling and experiments | |
| CN109918787B (en) | Finite volume method based simulation method for gas-liquid two-phase homogeneous mass flow in water delivery pipeline | |
| CN103353908B (en) | A kind of pipe resistance coefficient Method for Accurate Calculation based on numerical computations | |
| CN111339701B (en) | Pipeline leakage characteristic Godunov simulation method based on Brunone dynamic friction resistance | |
| Weinerowska-Bords | Alternative approach to convolution term of viscoelasticity in equations of unsteady pipe flow | |
| CN111985166A (en) | Pipeline hydraulic transient simulation method and storage medium with implicit consideration of dynamic friction resistance | |
| Zhou et al. | An accurate and efficient scheme involving unsteady friction for transient pipe flow | |
| CN113553737A (en) | Valve flow prediction method based on valve pressure difference | |
| CN110174237A (en) | The experiment porch of fluid state in a kind of measurement oil pipe | |
| CN111695307B (en) | Water hammer finite volume simulation method considering dynamic friction resistance explicitly | |
| Hullender | Alternative approach for modeling transients in smooth pipe with low turbulent flow | |
| CN106777770B (en) | The analogy method of hole stream in aqueduct based on finite volume method | |
| M GAD et al. | Impact of pipes networks simplification on water hammer phenomenon | |
| Pal et al. | Efficient approach toward the application of the Godunov method to hydraulic transients | |
| CN114811448B (en) | Method for pipeline leakage detection, leakage flow velocity estimation and leakage positioning under flowing condition | |
| Franchini et al. | Use of Torricelli's equation for describing leakages in pipes of different elastic materials, diameters and orifice shape and dimensions | |
| Li et al. | 1D-3D coupling investigation of hydraulic transient for power-supply failure of centrifugal pump-pipe system | |
| CN113468832B (en) | Two-dimensional simulation method for shock tube diaphragm rupture process based on overlapped moving grids | |
| Bazhlekov et al. | Numerical investigation of the dynamic influence of the contact line region on the macroscopic meniscus shape | |
| CN110608891B (en) | Cold flow test system of liquid rocket engine and parallel storage tank propellant conveying balance performance test method | |
| Koren | Euler flow solutions for transonic shock wave–boundary layer interaction | |
| CN112163722A (en) | Method and device for predicting gas supply state of natural gas pipe network | |
| Ozawa et al. | An indirect measurement method of transient pressure and flow rate in a pipe using steady state kalman filter | |
| CN110895202B (en) | Test device for researching influence of environmental factors on flow measurement |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |