Non-isothermal steady-state hydraulic calculation method, equipment, storage medium and program product for urban gas pipeline network based on Broyden iteration method
By combining the Broyden iteration method with the Jacobi matrix approximation matrix, the problem of temperature and flow direction changes not being considered in steady-state hydraulic calculations is solved, and high-precision and efficient calculation of gas pipeline network simulation is achieved.
Patent Information
- Application Number
- CN202310288866.6
- Authority / Receiving Office
- CN · China
- Patent Type
- Patents(China)
- Current Assignee / Owner
- Filing Date
- 2023-03-22
- Publication Date
- 2025-09-16
- Estimated Expiration
- 2043-03-22
AI Technical Summary
The existing steady-state hydraulic calculation method fails to fully consider the changes in temperature and flow direction, resulting in a deviation between the single-pipe hydraulic calculation model and the actual flow phenomenon, affecting the accuracy of gas pipeline network simulation.
The Broyden iteration method combined with the Jacobi matrix approximation matrix is used to perform hydraulic calculations on a single pipe, gradually updating the gas composition and temperature. A non-isothermal steady-state hydraulic calculation method is constructed, and the fourth-order Runge-Kutta method is used for single-pipe hydraulic calculations, and the relevant parameters are updated during the iterative process.
The calculation accuracy and efficiency of gas pipeline network simulation are improved, the accuracy and stability of the calculation results are ensured, and it is suitable for non-isothermal steady-state hydraulic calculations of urban gas pipeline networks.
Smart Images

Figure CN116305696B_ABST
Abstract
Description
Technical Field
[0001] The present application relates to the field of hydraulic calculation technology, and in particular to non-isothermal steady-state hydraulic calculation of urban gas pipeline networks. Background Art
[0002] Urban gas pipeline network simulation uses theoretical modeling and the high computing power of computers to replicate network operating conditions and guide network management. Gas pipeline network simulation primarily includes modules such as hydraulic calculation, resistance identification, leak and blockage diagnosis, and optimized scheduling. Hydraulic calculation is the foundation for all other modules, and its accuracy directly impacts the effectiveness of pipeline network simulation.
[0003] Hydraulic calculation of urban gas pipeline networks is an important part of gas pipeline network design and simulation. It has two important functions: (1) During the design process, hydraulic calculation of gas pipeline networks can determine the pipe diameter in combination with economic friction resistance, and further determine the preliminary budget for gas pipeline network construction; (2) During operation, hydraulic calculation of urban gas pipeline networks can determine the operating conditions of the entire gas pipeline network based on boundary conditions, providing support for leakage diagnosis, intelligent regulation, etc. of the pipeline network.
[0004] Gas pipeline network hydraulic calculations can be divided into steady-state and transient hydraulic calculations, depending on whether or not the impact of load instability is considered. Transient hydraulic calculations account for the effects of load fluctuations, introducing a time term to reflect the impact of load fluctuations during the modeling process. In contrast, steady-state hydraulic calculations ignore the impact of load fluctuations, resulting in a simpler hydraulic model and significantly improving simulation efficiency.
[0005] The steps of pipe network hydraulic calculation can be roughly divided into: determining the single-pipe hydraulic calculation model; constructing the pipe network equation group; and solving the equation group using numerical solutions.
[0006] At present, steady-state hydraulic calculation is still an important means of hydraulic calculation for urban gas pipeline networks. At present, the single-pipe hydraulic calculation model for steady-state hydraulic calculation mostly adopts analytical solutions. It is combined with other pipe network equations to construct a pipe network equation group, and then the pipe network equation group is solved by numerical solutions such as Newton iteration method to realize hydraulic calculation.
[0007] However, the disadvantages of the prior art are:
[0008] Single-tube analytical solution The momentum equation It is simplified and does not consider the impact of changes in temperature, elevation, etc. along the flow direction on gas parameters. Therefore, there will be a large deviation between the single-pipe hydraulic calculation model and the actual flow phenomenon. When combined with other equations, this deviation will also be brought in, ultimately resulting in insufficient accuracy of the hydraulic calculation. Summary of the Invention
[0009] The purpose of the present invention is to solve the problem of low accuracy of existing hydraulic calculations and to provide a non-isothermal steady-state hydraulic calculation method, equipment, storage medium and program product for urban gas pipeline networks based on the Broyden iteration method.
[0010] The present invention is achieved through the following technical solutions. In one aspect, the present invention provides a non-isothermal steady-state hydraulic calculation method for a town gas pipeline network based on the Broyden iteration method, the method comprising:
[0011] Step 1: Read the topology and working condition information of the pipe network;
[0012] Step 2: Set the initial values of the gas composition and node temperature of the pipe network;
[0013] Step 3: numbering the unknown quantities;
[0014] Step 4: Set the initial value of the unknown quantity iteration of the pipe network equation group;
[0015] Step 5: Perform hydraulic calculations on a single pipe;
[0016] Step 6: Based on the hydraulic calculation results of a single pipe obtained in step 5, list the pressure drop equation of the pipe section and combine it with the node flow balance equation to construct the equation group of the pipe network;
[0017] Step 7: Determine an approximate matrix of the Jacobi matrix of the pipe network equations according to the numbers of the unknown quantities;
[0018] Step 8: Based on the gas composition and temperature settings, the Broyden iteration method is used to solve the equations of the pipe network. During the iteration process, step 5 needs to be repeated to obtain the pressure of each node and the flow of each pipe section.
[0019] Step 9: Update the gas composition and temperature of the pipe network;
[0020] Step 10: Determine whether the errors after the two gas quality updates and temperature updates are within the set allowable values. If so, it is considered that the hydraulic calculation has reached the required accuracy and the calculation is stopped; otherwise, return to step 8 and re-calculate the hydraulic calculation.
[0021] Furthermore, step 2 includes:
[0022] When the hydraulic leveling of the pipe network is first performed, the gas quality information of each node is considered to be consistent, and it is preliminarily assumed that the gas composition of each node is either methane or the average of the gas quality components of each gas source;
[0023] The initial temperature of each node is set to the ambient temperature.
[0024] Furthermore, step 4 specifically includes:
[0025] Set the initial pressure value of the node and the initial flow value of the pipe section;
[0026] By default, the initial pressure value of all nodes is equal to the maximum pressure value of all nodes with known pressure;
[0027] The default initial flow rate is 0.01 times the minimum known flow rate of the pipe section. Assume that the minimum known flow rate of the pipe section is 0.1Nm 3 / s;
[0028] It is considered that the flow direction is from the point with smaller node number to the point with larger node number;
[0029] The sign of the flow at each node is defined as: negative for inflow to the node and positive for outflow from the node.
[0030] Furthermore, step 5 specifically includes:
[0031] Step 5.1: Determine the hydraulic calculation step length for a single pipe;
[0032] Step 5.2: Divide the pipe section into several micro-element segments according to the hydraulic calculation step length of the single pipe;
[0033] Step 5.3: Perform hydraulic calculations on each micro-element segment, using the starting point parameters of the pipe segment as the starting point parameters of the first micro-element segment, and update the gas parameters, including pressure, flow rate, and temperature, as conditions for the next calculation;
[0034] Step 5.4: Use the end parameter of the previous micro-element segment as the starting parameter of the next micro-element segment, and use the same method as step 5.3 to calculate the end parameter of each micro-element segment in turn;
[0035] Step 5.5: Use the end parameters of the last infinitesimal segment as the end parameters of the pipe segment, and the single-pipe hydraulic calculation is completed.
[0036] Furthermore, step 5.3 further includes:
[0037] Step 5.3.1. For the medium-pressure network, estimate the time s for the gas to return to ambient temperature after flowing out of the source;
[0038] Step 5.3.2: For pipe sections with a flow distance from the starting point of the pipe section to the gas source greater than s, the temperature gradient is no longer calculated during hydraulic calculations. Instead, the temperature is set to the ambient temperature.
[0039] Furthermore, step 7 specifically includes:
[0040] According to the number of the unknown quantity, use the formula:
[0041]
[0042] Where A0 is a high-order matrix, which is usually solved by solving the matrix equation A0X=E. The solved X is the numerical solution of B0. The initial matrix B0 is the inverse matrix of the approximate matrix A0 of the Jacobi matrix of the equation system; f i is the residual of the i-th equation in the system of equations for the pipe network;
[0043] Determine the approximate matrix of the Jacobi matrix of the pipe network equations.
[0044] Furthermore, step 9 specifically includes:
[0045] Starting from the gas source point, along the direction of fluid flow in the pipe, the gas composition and gas temperature of each node are gradually updated.
[0046] In a second aspect, the present invention provides a computer device comprising a memory and a processor, wherein the memory stores a computer program, and when the processor runs the computer program stored in the memory, the steps of a non-isothermal steady-state hydraulic calculation method for a town gas pipeline network based on the Broyden iteration method as described above are executed.
[0047] In a third aspect, the present invention provides a computer-readable storage medium, in which a plurality of computer instructions are stored, and the plurality of computer instructions are used to enable a computer to execute a non-isothermal steady-state hydraulic calculation method for a town gas pipeline network based on the Broyden iteration method as described above.
[0048] In a fourth aspect, the present invention provides a computer program product, which, when executed by a processor, implements the non-isothermal steady-state hydraulic calculation method for urban gas pipeline networks based on the Broyden iteration method as described above.
[0049] Beneficial effects of the present invention:
[0050] The present invention proposes to apply the single-pipe numerical solution method to the pressure drop equation, and adopts the Broyden iteration method with superlinear convergence to solve the pipe network equations. The pipe network equations and unknown variables are numbered, and the approximate matrix of the Jacobi matrix is determined. While ensuring the calculation efficiency, the calculation accuracy is taken into account at the same time, further improving the application value of steady-state hydraulics in gas pipe network simulation.
[0051] The present invention is applicable to non-isothermal steady-state hydraulic calculation of urban gas pipeline networks. BRIEF DESCRIPTION OF THE DRAWINGS
[0052] In order to more clearly illustrate the technical solution of the present application, the following is a brief introduction to the drawings required for use in the embodiments. Obviously, for ordinary technicians in this field, other drawings can be obtained based on these drawings without any creative work.
[0053] Figure 1 This is the flow chart of non-isothermal steady-state hydraulic calculation;
[0054] Figure 2 is the relative pressure drop error of the high-pressure pipe network (pipe diameter 0.7m, starting pressure 3.56MPa);
[0055] Figure 3 is the relative pressure drop error of the medium-high pressure pipe network (pipe diameter 0.7m, pressure 1.356MPa);
[0056] Figure 4 is the relative pressure drop error of the medium-pressure pipe network (pipe diameter 0.3m, pressure 0.356MPa);
[0057] Figure 5 Schematic diagram of fluid temperature change along the flow direction;
[0058] Figure 6 This is an example of a medium-pressure pipe network topology diagram (circles represent pipe segment numbers, and uncircled numbers represent node numbers);
[0059] Figure 7 Determine the process for single-pipe hydraulic calculation step length;
[0060] Figure 8 Update the topology graph for the parameters. DETAILED DESCRIPTION
[0061] The embodiments of the present invention are described in detail below, examples of which are shown in the accompanying drawings, wherein the same or similar reference numerals throughout represent the same or similar elements or elements having the same or similar functions. The embodiments described below with reference to the accompanying drawings are exemplary and are intended to be used to explain the present invention, but are not to be construed as limiting the present invention.
[0062] Embodiment 1: A non-isothermal steady-state hydraulic calculation method for a town gas pipeline network based on a Broyden iteration method, the method comprising:
[0063] Step 1: Read the topology and working condition information of the pipe network;
[0064] Step 2: Set the gas composition of the pipe network and the temperature of the nodes;
[0065] Step 3: Number the unknown quantities;
[0066] Step 4: Set the initial value of the unknown quantity iteration of the pipe network equation group;
[0067] Step 5: Use the fourth-order Runge-Kutta method to perform hydraulic calculations on a single pipe;
[0068] Step 6: Based on the hydraulic calculation results of a single pipe obtained in step 5, list the pressure drop equation of the pipe section and combine it with the node flow balance equation to construct the equation group of the pipe network;
[0069] Step 7: Determine an approximate matrix of the Jacobi matrix according to the number of the unknown quantity;
[0070] Step 8: Based on the gas quality parameters and temperature settings, the Broyden iteration method is used to solve the equations of the pipe network. During the iteration process, step 5 needs to be repeated to obtain the pressure of each node and the flow rate of each pipe section.
[0071] Step 9: Use the BFS breadth-first search algorithm to update the gas quality and temperature of the pipe network. The BFS breadth-first search algorithm is used in the hydraulic calculation of the gas pipe network and can be used to update the gas quality and temperature information of the pipe network to achieve higher hydraulic calculation accuracy.
[0072] Step 10: Determine whether the errors after the two gas quality updates and temperature updates are within the set allowable values. If so, it is considered that the hydraulic calculation has reached the required accuracy and the calculation is stopped; otherwise, return to step 8 and re-calculate the hydraulic calculation.
[0073] In this implementation, it is proposed to apply the single-pipe numerical solution method to the pressure drop equation, and use the Broyden iteration method with superlinear convergence to solve the pipe network equations. This method not only ensures computational efficiency but also takes into account computational accuracy, further improving the application value of steady-state hydraulics in gas pipe network simulation.
[0074] Implementation method 2: This implementation method further limits the non-isothermal steady-state hydraulic calculation method for a town gas pipeline network based on the Broyden iteration method described in implementation method 1. In this implementation method, step 2 is further limited, specifically including:
[0075] When the hydraulic leveling of the pipe network is first performed, the gas quality information of each node is considered to be consistent, and it is preliminarily assumed that the gas composition of each node is methane or the average value of the gas quality components of each gas source;
[0076] The initial temperature of each node is set to the ambient temperature, and the temperature of each node is assumed to be consistent with the ambient temperature to quickly complete the first hydraulic leveling of the pipe network.
[0077] In this embodiment, considering that the pipeline network may be a multi-gas source pipeline network, the gas quality information of each gas source is different. However, when the pipeline network is hydraulically leveled for the first time, the gas quality information of each pipe section is considered to be consistent, and the setting is made based on the main component of the gas being methane (or the average value of the gas quality components of each gas source); based on the characteristic of the pipeline network with a small flow rate, the temperature is set, and the initial setting of the temperature of each node is the ambient temperature.
[0078] Implementation method 3: This implementation method further limits the non-isothermal steady-state hydraulic calculation method for a town gas pipeline network based on the Broyden iteration method described in implementation method 1. In this implementation method, step 4 is further limited, specifically including:
[0079] Set the initial pressure value of the node and the initial flow value of the pipe section;
[0080] Used as the initial value for solving the pipe network equations using the Broyden iteration method.
[0081] By default, the initial pressure value of all nodes is equal to the maximum pressure value of all nodes with known pressure; (at least one node pressure parameter is known in the working condition)
[0082] The pipe networks in the same area have the same pressure level system and the pressure values of each node are very close. Using the known maximum pressure value as the initial value of the unknown pressure can make the initial pressure value close to the true value, thereby improving the convergence and convergence rate of the Broyden iteration method.
[0083] By default, all flow rates are 0.01 times the known minimum pipe flow rate. Since there are usually no flow measurement points in the pipe network (except for node flow measurement points), it is assumed that the known minimum pipe flow rate is 0.1Nm 3 / s;
[0084] This is mainly because in the pipeline network, there are always many users who do not use gas at the same time, making the flow of many pipe sections 0, so a smaller initial flow value is set.
[0085] It is considered that the flow direction is from the point with smaller node number to the point with larger node number;
[0086] This identification method is convenient for calculating uniformity (it can also be changed to flow from the starting point to the end point of the pipe segment)
[0087] The sign of the flow at each node is defined as: negative for inflow to the node and positive for outflow from the node.
[0088] It is mainly used to facilitate the column node flow balance equation.
[0089] Implementation method 4: This implementation method further limits the non-isothermal steady-state hydraulic calculation method for a town gas pipeline network based on the Broyden iteration method described in implementation method 1. In this implementation method, step 5 is further limited, specifically including:
[0090] Step 5.1: Determine the hydraulic calculation step length for a single pipe;
[0091] Step 5.2: Divide the pipe section into several micro-element segments according to the hydraulic calculation step length of the single pipe;
[0092] Step 5.3: Use the fourth-order Runge-Kutta method to perform hydraulic calculations on each micro-element segment, use the relevant parameters of the pipe segment starting point as the starting point parameters of the first micro-element segment, and update the relevant gas parameters (such as pressure, flow rate, temperature, etc.) as the conditions for the next calculation;
[0093] Step 5.4: Use the end parameter of the previous micro-element segment as the starting parameter of the next micro-element segment, and use the same method as step 5.3 to calculate the end parameter of each micro-element segment in turn;
[0094] Step 5.5: Use the end parameters of the last infinitesimal segment as the end parameters of the pipe segment, and the single-pipe hydraulic calculation is completed.
[0095] In this embodiment, the fourth-order Runge-Kutta method is used. This method continuously updates relevant parameters such as gas pressure, temperature, density, flow rate, viscosity, etc. during each iterative solution process to ensure the accuracy of single-pipe hydraulic calculations.
[0096] The step size in this embodiment is set to a value based on testing, and can be optimally set to 250 meters. With this step size, the microelement segmentation can effectively utilize the fourth-order Runge-Kutta method to improve the accuracy of non-isothermal steady-state hydraulic calculations for urban gas pipeline networks while also balancing computational efficiency.
[0097] Implementation method 5: This implementation method further limits the non-isothermal steady-state hydraulic calculation method for a town gas pipeline network based on the Broyden iteration method described in implementation method 4. In this implementation method, step 5.3 is further limited, specifically including:
[0098] Step 5.3 also includes:
[0099] Step 5.3.1. For the medium-pressure network, use the fourth-order Runge-Kutta method or the Sukhov formula to estimate the process s for the gas to return to the ambient temperature after flowing out from the source;
[0100] Step 5.3.2: For pipe sections with a flow rate greater than s from the starting point to the gas source, the temperature gradient is no longer calculated when the fourth-order Runge-Kutta method is used for hydraulic calculation. Instead, set the temperature to ambient temperature.
[0101] In most cases, the temperature of the medium-pressure pipeline gas source is consistent with the ambient temperature. However, there are also cases where LNG vaporization is used to achieve pipeline peak regulation, which will cause the gas source supply temperature to deviate from the ambient temperature. The medium-pressure pipeline flow is small. When there is a heat transfer temperature difference, the use of the fourth-order Runge-Kutta iteration process for a single-tube hydraulic meter is prone to the temperature dropping below absolute zero, which will make it impossible to update the parameters. Therefore, in order to avoid this situation, for pipe sections where there is a heat transfer temperature difference between the fluid and the environment, the calculation step size needs to be reduced; for pipe sections where the fluid temperature is consistent with the ambient temperature, the temperature gradient is no longer solved. The default calculation step size can be used. The step size for single-pipe hydraulic calculations is determined by assuming that the flow from the start point of the pipe segment to the gas source is greater than s. This ensures both feasibility and efficiency of the calculation.
[0102] Implementation method 6: This implementation method further limits the non-isothermal steady-state hydraulic calculation method of a town gas pipeline network based on the Broyden iteration method described in any one of implementation methods 1 to 5. In this implementation method, step 7 is further limited, specifically including:
[0103] Step 7 specifically includes:
[0104] According to the number of the unknown quantity, use the formula:
[0105]
[0106] Among them, A0 is a high-order matrix, which is usually solved by solving the matrix equation A0X=E. The solved X is the numerical solution of B0. The initial matrix B0 is the inverse matrix of the approximate matrix A0 of the Jacobi matrix of the equation system; f i is the residual of the i-th equation in the system of equations for the pipe network;
[0107] Determine the approximate matrix of the Jacobi matrix of the pipe network equations.
[0108] The determination method of this embodiment determines the elements of the matrix by this method. Under normal circumstances, it can ensure that the Jacobi approximate matrix A0 and the exact Jacobi matrix have more than 50% of the same non-zero elements, and the positive and negative signs of the remaining non-zero elements are also consistent with the corresponding elements of the Jacobi matrix, ensuring the convergence and convergence efficiency of the Broyden iteration method.
[0109] Implementation method 7: This implementation method further limits the non-isothermal steady-state hydraulic calculation method for a town gas pipeline network based on the Broyden iteration method described in implementation method 1. In this implementation method, step 9 is further limited, specifically including:
[0110] Step 9 specifically includes:
[0111] Starting from the gas source point, along the direction of fluid flow in the pipe, the gas composition and gas temperature of each node are gradually updated.
[0112] The updating method of this embodiment can update the gas composition and gas temperature at each point in the pipe network when calculating the non-isothermal steady-state hydraulics of the urban gas pipe network, so that it continuously approaches the actual gas composition and node temperature, further improving the accuracy of the hydraulic calculation.
[0113] The following is the working principle of the non-isothermal steady-state hydraulic calculation method for urban gas pipeline networks based on the Broyden iteration method as described above:
[0114] 2.2.1 Introduction to Single-Pipe Hydraulic Calculation Model
[0115] The single-pipe hydraulic calculation adopts the currently more mature numerical solution method - the fourth-order Runge-Kutta method. This method continuously updates the relevant parameters such as gas pressure, temperature, density, flow rate, viscosity, etc. during each iterative solution process to ensure the accuracy of the single-pipe hydraulic calculation. The specific steps are as follows:
[0116] By combining the continuity equation, momentum equation, and energy equation, we can obtain the following ordinary differential equation:
[0117]
[0118] in:
[0119]
[0120] Where: m is the mass flow rate of gas in the pipeline, kg / s; A is the cross-sectional area of the pipeline, m2; T e is the temperature of the environment where the pipeline is located, K, x is the flow distance, m; K is the heat transfer coefficient between the fluid and the environment, w / (m 2 ·K); p is pressure, pa; ρ is density, m3 / kg; T is fluid temperature, K; h is enthalpy, j / kg; θ is the angle between the pipe section and the horizontal plane, rad; λ is the friction coefficient, dimensionless; d is the pipe diameter, m; g is the gravity constant, m / s 2 ;π is the circumference ratio, dimensionless.
[0121] If the starting pressure p0, temperature T0, and mass flow rate m0 of a pipe section are known, the pipe section can be divided into N micro-units. Here, we take the i-th micro-unit pipe section as an example.
[0122] (1) Using the initial pressure p of the i-th pipe section i , temperature T i and flow m i, we can find the coefficients of the equation group (1.1), and then substitute them into the equation group (1.1) to solve and Then solve according to the following formula and (where Δx is the length of the infinitesimal element).
[0123]
[0124] (2) The newly obtained and Substitute into equation (1.1) and solve Use the following formula to calculate and
[0125]
[0126] (3) Using the newly acquired and Solve Solve using the following formula and
[0127]
[0128] (4) Using the newly obtained and Solve Then calculate the end pressure p of the infinitesimal segment i+1 , terminal temperature T i+1 and terminal flow m i+1 The calculation formula of Runge-Kutta method is as follows:
[0129]
[0130] The starting parameters of the first micro-element segment are the starting parameters of the pipe segment. The starting parameters of each subsequent micro-element segment are the end parameters of the previous micro-element segment. Starting from the first micro-element segment of the pipeline, each micro-element segment uses the process of steps (1) to (4) to calculate the gradients of pressure, temperature, and flow within each step. The calculated gradients are then substituted into formula (1.5) to complete the non-isothermal steady-state hydraulic calculation of a single step. The end parameters of each step are used as the starting parameters of the next step. The hydraulic calculation of a single pipe can be completed using the same steps as above.
[0131] The fourth-order Runge-Kutta method is used for calculation only when the gas flows from the starting point to the end point of the pipeline. In a ring network, due to the influence of working conditions, the gas flow direction of the same pipeline under different working conditions may be different, that is, the gas may flow from the end of the pipeline to the starting point.
[0132] When the flow flows from the starting point to the end point, use formulas (1.2) to (1.5) for calculation. When the flow flows from the end point (or end) to the starting point, use formulas (1.6) to (1.10) for calculation.
[0133] Therefore, a new parameter symbol is introduced to characterize the flow direction of the gas, which is defined as follows:
[0134]
[0135] If M represents the conditions at the beginning of a micro-element segment and N represents the conditions at the end of the micro-element segment, then the hydraulic calculation method for the micro-element segment in the ring network is as follows:
[0136] (1) Using the initial pressure p of the i-th pipe section i , temperature T i and flow m i , we can find the coefficients of the equation group (1.1), and then substitute them into the equation group (1.1) to solve and Then solve according to the following formula and (where Δx is the length of the infinitesimal element).
[0137]
[0138] (2) The newly obtained and Substitute into equation (1.1) and solve and Use the following formula to calculate and
[0139]
[0140] (3) Using the newly acquired and Solve and Solve using the following formula and
[0141]
[0142] (4) Using the newly obtained and Solve and Then calculate the end pressure p of the infinitesimal segment i+1 , terminal temperature T i+1 and terminal flow m i+1 The calculation formula of Runge-Kutta method is as follows:
[0143]
[0144] Considering that the gas pipeline network satisfies the continuity equation, that is, under specific working conditions, the mass flow rate of the pipe section does not change along the way, therefore, for a pipe section with a starting point 1 and an end point 2, if the pressure, temperature and flow rate at the starting point 1 of the pipe section are known to be P1, T1 and m1 respectively, the pressure at the end point 2 of the pipe section can be solved by formulas (1.1), (1.6) to (1-10) to obtain p2 = f p (p1, T1, m1), temperature T2 = f T (p1, T1, m1), p2 = f p (p1, T1, m1) and T2 = f T (p1, T1, m1) represents the function mapping relationship.
[0145] 2.2.2 Pipeline network equations
[0146] Taking into account the diversity of topology (branched, ring, single gas source and multiple gas sources, etc.), the pipe network equations are constructed by combining the pipe section pressure drop equation and the node flow balance equation.
[0147] If the starting point of the pipe section k is known to be i and the end point is j, the fourth-order Runge-Kutta method can be used to obtain the starting pressure p i , temperature T i and flow m i Calculate the end pressure f pk (p i ,T i ,m i ), then the pressure drop equation of the pipe section is:
[0148] (p i -p j )-(p i -f pk (p i ,T i ,m i ))=0 (1.11)
[0149] If the flow rates of the pipe sections connected to node i are known to be G1, G2, ...G n (Outflow is positive, inflow is negative), the node flow is Q i Then the node flow balance equation of this node is
[0150]
[0151] For a pipe network topology with N+1 nodes and B pipe sections, the parameters to be determined in the network equations are the nodal pressures of the N nodes (for a gas network, there must be one constant-pressure node, so the maximum number of unknown pressures is N) and the flow rates of the B pipe sections, resulting in a total of N+B unknowns. For this network topology, we can list B linearly independent pressure drop equations for the pipe sections, as shown in Equation (1.11), and N linearly independent nodal flow balance equations, as shown in Equation (1.12), for a total of N+B equations. In summary, the number of parameters to be determined in the network equations is equal to the number of equations, and the equations are solvable.
[0152] 2.2.3 Solving the Equations
[0153] For the pipe network equations listed in Section 2.2.2, considering that the hydraulic calculation of a single pipe in the pipe pressure drop equations is solved numerically, the Jacobi matrix of the pipe network equations cannot be solved. Therefore, the traditional Newton method is no longer applicable. The Broyden iteration method does not require an accurate Jacobi matrix and effectively solves the difficulties in numerical solutions in the equations.
[0154] The Broyden iteration method also has the characteristics of superlinear convergence and high computational efficiency. Its iteration format is:
[0155] x i+1 =x i -B i F(x i )
[0156]
[0157] where δ i+1 =x i -x i-1 ,Δ i =F(x i )-F(x i-1 ), x0=initial vector, B0=initial matrix. i represents the result of the i-th iteration, F(x i ) represents the variable x i The residuals of the equations are all purely mathematical definitions here and have no actual physical meaning.
[0158] Use equation (1.13) to solve the pipe network equations listed in Section 2.2.2.
[0159] The initial matrix B0 is the inverse matrix of the approximate matrix A0 of the Jacobi matrix of the equation system. The confirmation method of the inverse matrix elements adopts the qualitative and quantitative combination method to determine the determination method as follows:
[0160]
[0161] By determining the elements of the matrix in this way, it can usually be guaranteed that the Jacobi approximate matrix A0 has more than 50% of the same non-zero elements as the exact Jacobi matrix, ensuring the convergence and convergence efficiency of the Broyden iterative method.
[0162] In large-scale pipeline networks, A0 is a high-order matrix that cannot be directly inverted. It is usually solved by solving the matrix equation A0X=E (E is the unit matrix). The solved X is the numerical solution of B0.
[0163] 2.2.4 Non-isothermal steady-state hydraulic calculation process
[0164] After each hydraulic adjustment of the pipe network, the non-isothermal steady-state hydraulic calculation uses the BFS breadth-first search algorithm, starting from the gas source, to traverse each node in turn, updating the gas composition and temperature of the node as the parameters required for the next hydraulic adjustment of the pipe network. The gas composition and temperature errors of the two times are controlled within the allowable range. The pipe network parameters (node pressure, temperature, gas composition, etc.; flow rate of each pipe section) at this time can be considered to be true values. The non-isothermal steady-state hydraulic calculation flow chart is as follows: Figure 1 .
[0165] 2.2.5 Determination of spatial calculation step size
[0166] The fourth-order Runge-Kutta method is used for single-pipe hydraulic calculations, and the selection of the step size should take into account both calculation accuracy and calculation efficiency.
[0167] In the high-pressure pipe network, under normal circumstances, the flow rate of the high-pressure pipe network does not exceed 2×10 5 Nm 3 / h, that is 55.55Nm 3 / s; The relative error of the calculated results under different step lengths is compared with the pressure drop error when the step length is 0.1m. Figure 2 When the step length is controlled within 250m, the relative error of the pressure drop is small, so the calculation step length can be 250m.
[0168] In the secondary high pressure pipe network, the flow range is selected to be 1~10Nm 3 / s was tested, and the calculation results were compared with the pressure drop error when the step length was 0.1m. The relative error is as follows: Figure 3 As shown. Except the flow rate is 2Nm 3 / s, the error fluctuates greatly. When the step length is controlled within 250m, the relative error of the pressure drop is small. Therefore, the calculation step length can be 250m.
[0169] In the medium pressure pipe network, the flow range is selected to be 0.1~1Nm 3 / s was tested, and the calculation results were compared with the pressure drop error when the step length was 0.1m. The relative error is as follows: Figure 4 As shown. It is not difficult to find that the relative error is controlled within 10 -9 % and no obvious trend of change. If hydraulic calculations are performed separately for a medium-pressure network, the calculation step size can be determined based on the accuracy requirements. Considering that different levels of pipe networks may be calculated together during hydraulic calculations, the calculation step size can also be set to 250m for medium-pressure networks.
[0170] It should be noted that in medium-pressure pipe networks, when using the fourth-order Runge-Kutta method for single-pipe hydraulic calculations, when the flow rate is small, the temperature gradient during the calculation process is large, causing the calculation to fluctuate violently. It is possible that the temperature drops below zero K during a certain calculation. There are currently two ways to deal with this situation:
[0171] (1) Reduce the calculation step size. This method can effectively reduce the amplitude of the oscillation, but due to the large range of flow changes in the pipe section, when the flow rate changes from 10 -1 Nm 3 / s reduced to 10 -5 Nm 3 / s, the calculation step size needs to be reduced by 0.1m from the default calculation step size of 250m, which will inevitably lead to a significant increase in calculation time and affect the applicability of the method.
[0172] (2) Set the adaptive step size according to the topological position of the pipe section. In theory, the closer the pipe is to the gas source, the greater the temperature difference between the fluid in the pipe and the environment. When the flow rate is small, after flowing out of the gas source, the fluid temperature will be consistent with the ambient temperature after a small process. For a pipe with a diameter of 0.45m, a length of 150m, a starting temperature of 303.15K, an ambient temperature of 283.15K, and a starting pressure of 0.356MPa, its temperature along the flow direction is as follows: Figure 5
[0173] As shown, it is not difficult to find that for a flow rate of 0.01Nm 3 / s、0.001Nm 3 / s and 0.0001Nm 3 / s, its temperature will return to the ambient temperature within 40m of the process, and the flow rate is 0.1Nm 3 / s, the heat exchange temperature difference between its temperature and the ambient temperature is only about 3K after 150m of the process.
[0174] Therefore, when the gas source temperature is higher than the ambient temperature, the fourth-order Runge-Kutta method or the Sukhov formula can be used to first estimate the flow rate s for the gas to return to the ambient temperature after flowing out from the source (that is, the flow distance when the temperature drops to the ambient temperature after flowing out of the gas source). The hydraulic calculation step length of the pipeline network is determined as follows: ① For pipe sections where the flow rate from the starting point of the pipe section to the gas source is less than s, the step length adopts an adaptive step length (that is, the step length is determined based on temperature, flow rate, and pipe section parameters). ② For pipe sections where the flow rate from the starting point of the pipe section to the gas source is greater than s, the step length adopts the default calculation step length of 250m for the medium-pressure pipeline network, and when the fourth-order Runge-Kutta method is used for hydraulic calculation, the temperature gradient is no longer solved. Instead, the temperature is considered to be constant (i.e., the ambient temperature).
[0175] Using the above method (2) can not only minimize the impact of small flow rate and large heat transfer temperature difference on the calculation efficiency, but also reduce the negative impact of temperature oscillation on pressure solution to a certain extent, thereby improving the accuracy of the calculation.
[0176] Implementation 8: This implementation is an example of a non-isothermal steady-state hydraulic calculation method for a town gas network based on the Broyden iteration method as described above, specifically including:
[0177] for Figure 6 In the pipe network topology diagram, there are 35 pipe sections and 28 nodes. The gas source nodes are 0, 18, and 26. The gas supply temperature is 288.15K. The gas composition is node 0 (95% CH4, 5% C2H6), node 18 (93% CH4, 6% C2H6, C3H8), and node 26 (95% CH4, 5% C2H6).
[0178] If the pressure of gas source 26 is known, the flow rates of gas source 0 and gas source 18 are known, the node flow rates of all user nodes (8, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 27) are known, and the node flow rates of the middle nodes (nodes not directly connected to users and gas sources, nodes 1, 2, 3, 4, 5, 6, 7, 9, 10, 11, 12, 13, 14, 15) are all 0, then in this pipe network, there are 27 nodes (except gas source node 26) with unknown pressures and 35 pipes with unknown pressures. The flow rate of each segment is unknown (no flow measurement points exist on any segment), resulting in a total of 62 unknown quantities. By listing the pressure drop equation for each segment, we can list 35 independent equations. By listing the node flow balance equations for all nodes except node 26, we can list 27 independent equations, resulting in a total of 62 network equations. It is not difficult to find that the number of unknown quantities is equal to the number of equations, and the equations are independent. Therefore, the system of equations can be solved, and non-isothermal steady-state hydraulic calculation can be performed. The specific steps of the ring network hydraulic calculation process are as follows:
[0179] (1) Read the topological structure and operating condition information of the pipeline network.
[0180] (2) Assume the gas composition of each node (or each pipe section) and the temperature of the node. Considering that the pipeline network is a multi-gas source network, the gas composition of each gas source is different. However, when the pipeline network is hydraulically leveled for the first time, the gas composition of each pipe section is considered to be consistent. Since the main component of the gas is methane, it can be preliminarily assumed that the gas composition of each node and each pipe section is methane. Due to the small flow rate of the pipeline network, it is initially assumed that the temperature of each node is the ambient temperature of 283.15K, that is, T0 = T 18 =T 26 =T1=T2=……T 17 =T 19 =……T 25 =T 27 =283.15K.
[0181] (3) Number the unknown quantities. Figure 6 In [x0,x1,x2,…,x 61 ]=[P0, P1……P 25 、P 27 ,G0=G1=G2=……=G 34 ]. That is, assigning the pipe network variables to the variables of the equation group.
[0182] (4) Determine the initial value of the unknown quantity iteration of the equation group. Assume that the pressure of the node P0, P1...P 25 、P 27 , the flow rate of the pipe section G1, G2, G3...G 34 , by default, the pressure of all nodes is equal to the maximum value of all known nodes, that is, P0=P1=……P 25 =P 27 =P 26 By default, all flow rates are 0.01 times the minimum known flow rate of the pipe section. Since there is no flow measurement point in the pipe network, the minimum flow rate is assumed to be 0.1Nm 3 / s, that is, G0=G1=G2=……=G 34 =0.001Nm 3 / s, and the flow direction is assumed to be from the node with a smaller node number to the node with a larger node number. For example, pipe segment 15 flows from node 12 to node 13. In addition, the sign of the flow rate at each node is defined as negative for inflow into the node and positive for outflow from the node.
[0183] (5) Single-pipe hydraulic calculation. Single-pipe hydraulic calculation is performed using the single-pipe hydraulic calculation steps. According to the pressure of the nodes P0, P1...P 25 、P 27 And the flow rate of the pipe section G0, G1, G2...G 34Perform hydraulic calculation on a single pipe and solve the end pressure of each pipe section as follows:
[0184] The single-pipe hydraulic calculation step specifically includes:
[0185] This example uses Figure 6 The pipe network topology shown in the figure explains the steps of non-isothermal steady-state hydraulic calculation. Non-isothermal steady-state hydraulic calculation includes two parts: single-pipe hydraulic calculation and ring network hydraulic calculation.
[0186] Single pipe hydraulic calculation is to use the fourth-order Runge-Kutta method to perform hydraulic calculation on a single pipe. The flow chart is as follows Figure 7 The detailed steps are as follows:
[0187] (5.1) Determine the hydraulic calculation step size step for a single pipe. Considering that this example only performs hydraulic calculations on a medium-pressure pipe network, the default maximum calculation step size can be selected based on the actual accuracy requirements. Let's set the default maximum calculation step size for a single pipe as max_step. First, use the fourth-order Runge-Kutta method or Sukhov formula to estimate the flow rate s of the gas after it flows out from the source to the ambient temperature (that is, the flow distance when the temperature drops to the ambient temperature after flowing out of the gas source). Then, based on Figure 7 The flowchart shown determines the calculation step size of the pipe segment.
[0188] (5.2) Divide the pipe segment into N micro-segments based on the calculation step size of the single pipe. Considering that in most cases the pipe length Length cannot be divided evenly by the step size step, it is necessary to divide the pipe segment into a series of micro-segments based on the actual situation.
[0189] If the pipe segment is an integer multiple of the step length, the number of micro-element segments N = Length / step can be directly calculated; if the pipe length is not an integer multiple of the step length, the pipe segment is divided into N segments, where the length of the first N-1 segments is step, and the length of the Nth segment is Length-step*(N-1).
[0190] (5.3) Use formulas (1.1) to (1.10) to perform hydraulic calculations on the first microelement of the pipeline. Use the starting parameters of the pipe segment as the starting parameters of the first microelement, and calculate P each time. i j 、T i j 、m i j (j=1,2,3) After that, the gas density ρ, enthalpy h, friction λ and other physical parameters are updated as the conditions for the next calculation.
[0191] (5.4) Take the end parameter of the previous micro-element segment as the starting parameter of the next micro-element segment, and use the same method as (3) to calculate the end parameter of each micro-element segment in turn.
[0192] (5.5) The end parameters of the last infinitesimal segment are used as the end parameters of the pipe segment, and the single-pipe hydraulic calculation is completed.
[0193] (6) List the equations for the pipe network. Based on the results of the single-pipe hydraulic calculation in step (5), list the pipe section pressure drop equations. Combined with the node flow balance equations, the equations for the pipe network are constructed:
[0194]
[0195] (7) Determine the Jacobi matrix approximation matrix A0. According to the numbering of the unknowns in step (4), use formula (1.14) to determine each element in A0, for example, [A0] 00 =1, [A0] 01 =-1, [A0] 27 =-1.
[0196] (8) Hydraulic balance of the pipe network. Based on the assumptions of gas quality parameters and temperature in step (2), Broyden iteration method is used to solve equation (3.1) to solve the pressure P0, P1...P of each node. 25 、P 27 And the flow rate of each pipe section G0, G1, G2...G 34 .
[0197] (9) Use the BFS breadth-first search algorithm to update the gas quality and temperature of the pipeline network. Starting from the gas source point, along the direction of fluid flow in the pipe, gradually update the gas quality parameters and gas temperature of each node.
[0198] The gas quality update method is: start from the gas source point, the gas quality component y of each node is the weighted average value of the gas quality parameters of the pipe section flowing into the node with respect to the pipe section flow rate. Figure 8 In the example, if the upstream pipe sections of node 1 are known to be pipe section 1 and pipe section 2, and the flow rates of pipe section 1 and pipe section 2 are G Pipe1 and G pipe2 Then the temperament component of node 1 is The gas in each pipe section has the same gas composition as that of its upstream node. Figure 8 middle, y Pipe3 =y Pipe4 =y Pipe5 =y Node1 .
[0199] The temperature update method is: the temperature of each node is the weighted average of the end temperature of all pipes flowing into the node with respect to the flow rate. Figure 8In the example, if the end temperatures of pipe section 1 and pipe section 2 are T Pipe1 、T Pipe2 , the temperature of node 1 is The terminal temperatures of pipe sections 3, 4, and 5 need to be updated using a fourth-order Runge-Kutta method.
[0200] (10) Determine whether the error after two gas updates (or checks) and temperature updates (or checks) is within the set allowable value. The flow chart is as follows: Figure 1 . First, determine the number of times step (9) is performed. If the number of times step (9) is performed is 1, considering that the gas quality parameters and temperature are updated for the first time, and the gas quality and temperature initially assumed are not reference-based, then return to step (5) and re-calculate the hydraulics. If the number of times step (9) is performed is greater than 1, it is necessary to compare the gas quality error and temperature error before and after to determine whether they are within the set allowable error; if the temperature error of a node in the two calculations before and after exceeds the set allowable temperature error r T Or the temperament error exceeds the allowable temperament error r y , then return to step (5) and recalculate the hydraulics; if the temperature error of all nodes exceeds the set allowable temperature error r after two calculations, T Or the temperament error exceeds the allowable temperament error r y , it is considered that the hydraulic calculation has reached the required accuracy and the calculation is stopped.
Claims
1. A non-isothermal steady-state hydraulic calculation method for urban gas pipeline network based on Broyden iteration method, characterized by: The method comprises: Step 1: Read the topology and working condition information of the pipe network; Step 2: Set the initial values of the gas composition and node temperature of the pipe network; Step 3: Number the unknown quantities; Step 4: Set the initial value of the unknown quantity iteration of the pipe network equation group; Step 5: Perform hydraulic calculations on a single pipe; Step 6: Based on the hydraulic calculation results of a single pipe obtained in step 5, list the pressure drop equation of the pipe section and combine it with the node flow balance equation to construct the equation group of the pipe network; Step 7: Determine an approximate matrix of the Jacobi matrix of the pipe network equations according to the numbers of the unknown quantities; Step 8: Based on the gas composition and temperature settings, the Broyden iteration method is used to solve the equations of the pipe network. During the iteration process, step 5 needs to be repeated to obtain the pressure of each node and the flow of each pipe section. Step 9: Update the gas composition and temperature of the pipe network; Step 10: Determine whether the errors between the two gas quality updates and the temperature update are within the set allowable value. If so, it is considered that the hydraulic calculation has reached the required accuracy and the calculation is stopped; otherwise, return to step 8 and re-calculate the hydraulic calculation; The approximate matrix method for determining the Jacobi matrix of the pipe network equations in step 7 is: According to the number of the unknown quantity, use the formula: Where, Is a high-order matrix, usually solved by the matrix equation The solution is that is The numerical solution of the initial matrix is the approximate matrix of the Jacobi matrix of the equation system The inverse matrix of is the residual of the i-th equation in the system of equations for the pipe network; Determine the approximate matrix of the Jacobi matrix of the pipe network equations.
2. The non-isothermal steady-state hydraulic calculation method for a town gas network based on the Broyden iteration method according to claim 1 is characterized in that: Step 2 includes: When the hydraulic leveling of the pipe network is first performed, the gas composition of each node is considered to be consistent. It is preliminarily assumed that the gas composition of each node is either methane or the mean of the gas composition of each gas source. The initial temperature of each node is set to the ambient temperature.
3. The non-isothermal steady-state hydraulic calculation method for a town gas network based on the Broyden iteration method according to claim 1 is characterized in that: Step 4 specifically includes: Set the initial pressure value of the node and the initial flow value of the pipe section; By default, the initial pressure value of all nodes is equal to the maximum pressure value of all nodes with known pressure; The default initial flow rate is 0.01 times the minimum known flow rate of the pipe section. Assume that the minimum known flow rate of the pipe section is 0.1Nm 3 / s; It is considered that the flow direction is from the point with smaller node number to the point with larger node number; The sign of the flow at each node is defined as: negative for inflow to the node and positive for outflow from the node.
4. The non-isothermal steady-state hydraulic calculation method for a town gas network based on the Broyden iteration method according to claim 1 is characterized in that: Step 5 specifically includes: Step 5.1: Determine the hydraulic calculation step length for a single pipe; Step 5.2: Divide the pipe section into several micro-element segments according to the hydraulic calculation step length of the single pipe; Step 5.3: Perform hydraulic calculations on each micro-element segment, using the starting point parameters of the pipe segment as the starting point parameters of the first micro-element segment, and update the gas parameters, including pressure, flow rate, and temperature, as conditions for the next calculation; Step 5.4: Use the end parameter of the previous micro-element segment as the starting parameter of the next micro-element segment, and use the same method as step 5.3 to calculate the end parameter of each micro-element segment in turn; Step 5.5: Use the end parameters of the last infinitesimal segment as the end parameters of the pipe segment, and the single-pipe hydraulic calculation is completed.
5. The non-isothermal steady-state hydraulic calculation method for a town gas network based on the Broyden iteration method according to claim 4 is characterized in that: Step 5.3 also includes: Step 5.3.1: For the medium pressure network, estimate the process of gas returning to ambient temperature after flowing out from the source ; Step 5.3.2: For the flow from the start point of the pipe section to the gas source, For pipe sections with , when performing hydraulic calculations, the temperature gradient is no longer solved, but the temperature is set to the ambient temperature.
6. The non-isothermal steady-state hydraulic calculation method for a town gas network based on Broyden iteration method according to claim 1, characterized in that: Step 9 specifically includes: Starting from the gas source point, along the direction of fluid flow in the pipe, the gas composition and gas temperature of each node are gradually updated.
7. A computer device comprising a memory and a processor, wherein a computer program is stored in the memory, wherein: When the processor runs the computer program stored in the memory, the steps of the method according to any one of claims 1 to 6 are performed.
8. A computer-readable storage medium, characterized in that The computer-readable storage medium stores a plurality of computer instructions, and the plurality of computer instructions are used to enable a computer to execute the method according to any one of claims 1 to 6.
9. A computer program product, characterized in that When the computer program is executed by a processor, the method according to any one of claims 1 to 6 is implemented.
Citation Information
Patent Citations
Implementation method of dynamic simulation process of natural gas pipeline network conveying system
CN102779199A
Flow simulation and transient well analysis method based on generalized pipe flow seepage coupling
CN113826099A