CN112257355A - Modeling method for medium-low pressure gas distribution pipe network of hydrogen-doped natural gas - Google Patents

Modeling method for medium-low pressure gas distribution pipe network of hydrogen-doped natural gas Download PDF

Info

Publication number
CN112257355A
CN112257355A CN202011191416.8A CN202011191416A CN112257355A CN 112257355 A CN112257355 A CN 112257355A CN 202011191416 A CN202011191416 A CN 202011191416A CN 112257355 A CN112257355 A CN 112257355A
Authority
CN
China
Prior art keywords
hydrogen
pipeline
natural gas
node
formula
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.)
Pending
Application number
CN202011191416.8A
Other languages
Chinese (zh)
Inventor
周苏洋
孙锴宇
顾伟
陈晓刚
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Southeast University
State Grid Zhejiang Electric Power Co Ltd
Original Assignee
Southeast University
State Grid Zhejiang Electric Power Co Ltd
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Southeast University, State Grid Zhejiang Electric Power Co Ltd filed Critical Southeast University
Priority to CN202011191416.8A priority Critical patent/CN112257355A/en
Publication of CN112257355A publication Critical patent/CN112257355A/en
Pending legal-status Critical Current

Links

Images

Classifications

    • 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
    • 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
    • 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/08—Thermal analysis or thermal optimisation
    • 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)
  • Geometry (AREA)
  • Computer Hardware Design (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Evolutionary Computation (AREA)
  • General Engineering & Computer Science (AREA)
  • Computational Mathematics (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • Algebra (AREA)
  • Computing Systems (AREA)
  • Fluid Mechanics (AREA)
  • Mathematical Physics (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

The invention provides a method for building a medium-low pressure gas distribution pipe network of hydrogen-doped natural gas, which comprises the following steps: s1, establishing a steady-state isothermal model of the hydrogen-doped natural gas network pipeline, wherein the steady-state isothermal model comprises parameter settings of density, Reynolds number, dynamic viscosity and friction coefficient of the hydrogen-doped natural gas and mathematical model descriptions of nodes, pipelines and loops of the hydrogen-doped natural gas network; s2, based on the mathematical model description in the step S1, establishing a solving method of the medium-low pressure gas distribution pipe network of the hydrogen-doped natural gas by adopting an iteration method; s3, according to the solving method in the step S2, mathematical solving models of node pressure and pipeline flow of the network are respectively established by adopting a node method, a ring network method and a ring energy method, internal relations among the mathematical solving models obtained by comparing the three methods are analyzed, characteristics and application occasions of the mathematical solving models are compared and analyzed, and finally the modeling method of the medium-low distribution gas pipeline network of the hydrogen-doped natural gas is realized.

Description

Modeling method for medium-low pressure gas distribution pipe network of hydrogen-doped natural gas
Technical Field
The invention belongs to the technical field of mathematical modeling of energy systems, and particularly relates to a modeling method for a medium-low pressure gas distribution pipe network of hydrogen-doped natural gas.
Background
With the continuous development of society, the standards of end users for the quality of energy are continuously improved, and further requirements for the high efficiency and cleanness of energy are provided. In recent years, with the vigorous popularization of electric energy and new energy technologies in China, especially the proposal of policies of ' two substitution ' -electric energy substitution ' and ' clean substitution ', energy transformation is accelerated step by step, so that how to improve the consumption level of renewable energy and how to improve the energy use structure are important research directions of future energy networks.
On one hand, with the increasing attention of people to the environmental pollution problem and the practical problems of energy storage and the like, in the traditional fossil energy, the consumption proportion of coal and petroleum is continuously reduced, and the consumption proportion of natural gas as clean fossil energy is gradually increased; on the other hand, with the development of clean energy utilization technology, hydrogen as another clean energy is also considered to be an important component in the future clean energy consumption structure. However, hydrogen gas, as an ideal clean energy source, has the advantages of high energy density, wide flammability range, low minimum ignition energy, high combustion rate, no pollution emission generated by combustion and the like under the same weight, but also has the problems of influencing unit efficiency, being difficult to store and transport and the like. The hydrogen can be transported only after being compressed, the transportation difficulty is high, and the cost is high. And hydrogen is easy to gasify or fire, the technical development of the pure hydrogen combustion engine is immature, and the development of hydrogen fuel is greatly limited due to the lack of corresponding fuel infrastructure and fuel supply equipment. The hydrogen is doped into the natural gas to prepare the hydrogen-doped natural gas, so that the hydrogen-doped natural gas is a feasible mode for transition to pure hydrogen energy, the hydrogen-doped natural gas has higher combustion rate and less pollutant emission compared with the natural gas, and compared with the hydrogen, the hydrogen-doped natural gas can be transported by utilizing the existing natural gas pipe network, the use limit is smaller, the safety coefficient is higher, and the hydrogen-doped natural gas is a novel mixed clean energy which is more in line with the current actual requirement.
Meanwhile, the hydrogen and the natural gas have different physical and chemical properties, the transportation conditions of the hydrogen and the natural gas in the existing gas distribution pipe network are different, the physical model for the transportation of the hydrogen-doped natural gas pipe section cannot be accurately and effectively reflected by the original natural gas pipe network model, and the research on the transportation characteristics of the hydrogen-doped natural gas is very important work. Although the existing literature demonstrates the feasibility of transporting the hydrogen-doped natural gas by using the existing natural gas pipeline, the mathematical modeling of the medium-and-low-pressure gas distribution pipe network for transporting the hydrogen-doped natural gas is still in an exploration stage, so that the key of the problem lies in providing an effective and accurate mathematical modeling method for the medium-and-low-pressure gas distribution pipe network of the hydrogen-doped natural gas.
Disclosure of Invention
In order to effectively depict the pipe section transportation characteristics of the hydrogen-doped natural gas, the invention provides a modeling method of a medium-low pressure gas distribution pipe network of the hydrogen-doped natural gas, and provides the following technical scheme:
a method for building a medium-low distribution pipe network of hydrogen-doped natural gas comprises the following steps:
s1, establishing a steady-state isothermal model of the hydrogen-doped natural gas network pipeline, wherein the steady-state isothermal model comprises parameter settings of density, Reynolds number, dynamic viscosity and friction coefficient of the hydrogen-doped natural gas and mathematical model descriptions of nodes, pipelines and loops of the hydrogen-doped natural gas network;
s2, based on the mathematical model description in the step S1, establishing a solving method of the medium-low pressure gas distribution pipe network of the hydrogen-doped natural gas by adopting an iteration method;
s3, according to the solving method in the step S2, mathematical solving models of node pressure and pipeline flow of the network are respectively established by adopting a node method, a ring network method and a ring energy method, internal relations among the mathematical solving models obtained by comparing the three methods are analyzed, characteristics and application occasions of the mathematical solving models are compared and analyzed, and finally the modeling method of the medium-low distribution gas pipeline network of the hydrogen-doped natural gas is realized.
Further, the specific process of setting the parameters of the density, Reynolds number, kinematic viscosity and molar group coefficient of the hydrogen-doped natural gas is as follows:
the density rho of the hydrogen-doped natural gas is relative density delta and air density rhoAir (a)The specific expression is as follows:
Figure BDA0002752856330000021
in formula (1): rho0The dry air density in the standard state (273.15K, 101.3kPa) was 1.293kg/m3(ii) a p represents the pressure intensity of the hydrogen-doped natural gas and has the unit of MPa; t represents the temperature of the hydrogen-doped natural gas and has the unit of K;
the density ρ of the hydrogen-doped natural gas can be obtained by the formula (2), that is:
Figure BDA0002752856330000022
in the formula (2), the reaction mixture is,
Figure BDA0002752856330000023
represents CH4The density of (a) of (b),
Figure BDA0002752856330000024
represents CH4The volume of (a) to (b),
Figure BDA0002752856330000025
represents H2The density of (a) of (b),
Figure BDA0002752856330000026
represents H2The volume of (a);
the hydrogen volume ratio σ% was set to 10%, that is:
Figure BDA0002752856330000027
the characterization formula of the Reynolds number Re of the hydrogen-doped natural gas is as follows:
Figure BDA0002752856330000028
in formula (4): gamma denotes the kinematic viscosity of the gas in m2S; mu is dynamic viscosity in the unit of N.s/m2(ii) a v represents the gas flow rate in m/s; v represents the gas volume flow rate in m3S; m represents the gas mass flow rate, and the unit is kg/s; d represents the inner diameter of the pipeline and has the unit of m;
the calculation expression for the kinematic viscosity μ is as follows:
Figure BDA0002752856330000031
in formula (5): Δ represents the relative density of the hydrogen-doped natural gas under the standard state, namely the ratio of the air density under the standard state; rho is the gas density in kg/m3(ii) a T represents the gas temperature in K;
the friction coefficient lambda of the hydrogen-doped natural gas is calculated by adopting a Colebrook equation, namely:
Figure BDA0002752856330000032
in formula (6): k represents the absolute roughness of the inner wall of the tube in mm; d represents the inner diameter of the pipe, mm; re represents the Reynolds number of the hydrogen-loaded natural gas.
Further, in step S1, the mathematical model description process of the nodes, pipes, and loops of the hydrogen-loaded natural gas network includes:
two basic equations are constructed to describe the network state: a node-pipeline correlation matrix A and a loop-pipeline correlation matrix BETA;
element a of the node-pipe incidence matrix AijDeducing according to a kirchhoff first law, wherein the construction condition is that the pressure of a reference node is known;
Figure BDA0002752856330000033
element b of said loop-pipeline association matrix BETAnjDeducing according to a kirchhoff second law, and assuming that the reference flow direction of the loop is clockwise;
Figure BDA0002752856330000041
simulating a node voltage solution and a loop current solution in the power system, wherein the node-pipeline correlation matrix A and the loop-pipeline correlation matrix BETA are both matrices describing the state of the pipe network; according to the graph theory, there is an orthogonal relationship between node-pipeline correlation matrix a and loop-pipeline correlation matrix beta, that is, the inner product of the row vector of node-pipeline correlation matrix a and the row vector of loop-pipeline correlation matrix beta is always equal to zero, and the specific expression is as follows:
A×BT＝B×AT＝0 (9)
for the hydrogen-doped natural gas network, assuming that the number of nodes is i (including reference nodes), the number of pipelines is j, and the number of loops is n, the corresponding relationship among the number of nodes i, the number of pipelines j, and the number of loops n can be obtained as follows:
j＝n+i-1 (10)
aiming at a hydrogen-doped natural gas network, the following equation is used for description;
obtaining a node flow equation according to a kirchhoff first law, wherein the flow of the inflow node is equal to the flow of the outflow node, and the pressure of a reference node is known:
A(i-1)×jVj×1＝q(i-1)×1 (11)
in the formula (11), q(i-1)×1Representing the traffic of the egress node, A(i-1)×jRepresenting a loop-pipe correlation matrix, V, with reference nodes removedj×1Representing the flow vector of each pipe section;
the loop pressure drop equation is derived from kirchhoff's second law, with a pressure drop along any one closed loop of 0:
Bn×jΔpj×1＝0n×1 (12)
in the formula (12), Bn×jRepresenting a loop-pipe correlation matrix with n loops and j pipes, Δ pj×1Representing the pressure difference between the upper end and the lower end of each pipe section;
according to the pressure drop equation of the pipeline, the relation between the pressure drop of the pipeline and the pressure difference of the upper end and the lower end of the pipeline is as follows:
Figure BDA0002752856330000042
in formula (13): p is a radical of(i-1)×1Expressed as the pressure of each node relative to a reference node;
due to A x BT＝B×AT0, so
Figure BDA0002752856330000051
Therefore, the formula can be expressed by one of the two formulas;
when the pipe diameter is known, each pipeline has two unknowns of pressure drop delta p and flow, namely 2j unknowns, and the listed equation number is j + (i-1) + n ═ 2j for solving;
secondly, according to the corresponding quantitative relation between the pipeline pressure drop delta p and the pipeline flow, the equation of the pipeline flow pressure drop is obtained:
Δp＝G|V|·V (14)
in formula (14): g represents a diagonal element matrix of the flow resistance coefficient of the pipeline, namely the flow resistance coefficient of each pipeline section is on the corresponding diagonal line, and the rest elements of the matrix are 0; v represents the calculated flow vector of the pipeline in m3/h；
For medium-low pressure pipe network, the pipeline pressure drop delta p between the connection node a and the connection node babAnd the pipe flowThe calculation formula of the dynamic resistance coefficient G is as follows:
Figure BDA0002752856330000052
in formula (15): Δ pabRepresents the pressure drop of the pipe connecting nodes a and b in kPa2；paAnd pbRepresenting the pressure at nodes a and b in kPa; lambda represents the hydraulic friction coefficient of the pipeline; znRepresenting the compression factor of the gas under the standard condition, and taking 1; d represents the inner diameter of the pipeline and has the unit of mm; rho represents the density of the natural gas doped with hydrogen and has the unit of kg/m3；TnThe standard state absolute temperature is represented, namely 273.15K; t represents the temperature of the pipeline and has the unit of K; l represents the length of the pipeline, and the unit is km; z represents the compression factor of the hydrogen-doped natural gas, and the expression of Z is as follows:
Figure BDA0002752856330000053
in formula (16): t iscRepresents the critical temperature of the gas, in K; p is a radical ofcRepresents the critical pressure of the gas in Pa; p is a radical ofavDenotes the mean pressure of the pipe mn in Pa, pavThe expression is as follows:
Figure BDA0002752856330000054
in the formula (17), pmRepresenting the nodal pressure, p, of the beginning m of the pipe sectionnRepresenting the node pressure at the end n of the pipe section;
let element S in diagonal matrix S of j x jjjSatisfies the following conditions:
Figure BDA0002752856330000061
then, according to the equations (11) to (14) and (18), the following relationships can be obtained:
Figure BDA0002752856330000062
ASATp＝q (20)
B(G|V|V)＝0 (21)
further, in step S2, the solving method for establishing the medium-low pressure distribution pipe network in the hydrogen-doped natural gas by using the iterative method specifically includes:
suppose there are N variables x1,x2,L,xNThe equation is known:
Fi(x1,x2,L,xN)＝0,i＝1,2,L,N (22)
let variable matrix X ═ X1 x2 … xi]', unbalance matrix δ X ═ δ X1 δx2 … δxi]', function matrix F ═ F1 F2 … Fi]', the square matrix of unbalance amount deltaX2＝[(δx1)2 (δx2)2 … (δxi)2]', then each FiThe Taylor series expansion is adopted to obtain:
Figure BDA0002752856330000063
let the partial derivatives form a Jacobian matrix J, where element J is number (i, J)ijIs composed of
Figure BDA0002752856330000064
Then:
F(X+δX)＝F(X)+J×δX+O(δX2) (25)
neglecting O (δ X)2) And let F (X + δ X) ═ 0, a linear equation system for the unbalance amount δ X can be obtained:
J×δX＝-F(X) (26)
in formulae (25) and (26): deltax lets each function approach zero at the same time,
the amount of unbalance is added to the solution vector to continuously correct the initial value,
Xk+1＝Xk+δXk (27)
in formula (27): xk+1Represents the result of the iteration of the (k + 1) th solution vector, XkRepresents the result of the k-th iteration of the solution vector, deltaXkCorrection values representing the solution vector for the k-th iteration.
Further, in step S3, the specific process of establishing the mathematical solution model of the node pressure and the pipeline flow of the network by using the node method is as follows:
transform equation (20) into:
ASATp-q＝0 (28)
equation (28) is written as:
F(p)＝ASATp-q＝0 (29)
firstly, giving an initial estimation value of each node pressure, and secondly, obtaining the unbalance of the node pressure as an error correction quantity by continuously correcting the estimation value:
Jp×δpk＝-F(pk) (30)
in formula (30): j. the design is a squarepRepresenting a Jacobian matrix formed by partial derivation of the node pressure p; p is a radical ofkRepresenting the node pressure value of the kth iterative solution;
Figure BDA0002752856330000071
the iterative process for correcting the node pressure estimate is:
pk+1＝pk+δpk (32)
in the formula (32), pk+1Representing the node pressure value of the (k + 1) th iteration solution; p is a radical ofkRepresenting the node pressure value of the kth iterative solution; δ pkRepresenting the node pressure correction value of the kth iterative solution;
the iteration termination condition is as follows:
|δpk|≤ε (33)
in formula (33): ε represents the calculation accuracy.
In each iteration, the value of S varies according to equations (13), (14) and (18), and the sign of V is continuously corrected according to the sign of Δ p.
Further, in step S3, the specific process of establishing the mathematical solution model of the node pressure and the pipeline flow of the network by using the ring network method is as follows:
according to equation (21), the equation can be written as:
F(V)＝BG|V|V＝0 (34)
equation (34) can be converted to:
BG|V+BTe|(V+BTe)＝0 (35)
in the equation (35), e represents a correction amount added to the branch flow rate to correct the branch flow rate estimate;
then:
F(e)＝BG|V+BTe|(V+BTe)＝0 (36)
firstly, giving an initial value of loop flow, and then continuously iterating until a real solution is reached; in the iterative solution process of the ring network method, f (e) is not equal to zero, that is, a loop flow error is introduced into each loop, and accordingly, the unbalance amount of the loop flow can be obtained, namely, the unbalance amount can be regarded as an error correction amount:
Je×δek＝-F(ek) (37)
in formula (37): j. the design is a squareeRepresenting a Jacobian matrix formed by performing partial derivation on the flow correction value e of the branch circuit; e.g. of the typekRepresenting a branch flow correction value of the kth iterative solution; delta ekRepresenting the corrected value of the branch flow correction value of the kth iterative solution;
Figure BDA0002752856330000081
the iterative process for correcting the loop flow estimate is:
ek+1＝ek+δek (39)
the iteration termination condition is as follows:
|δek|≤ε (40)
in formula (40): epsilon is the calculation precision;
in the process JeThe binding e is changed continuously.
Further, the specific process of establishing the mathematical solution model of the node pressure and the pipeline flow of the network by using the ring energy method is as follows:
according to equations (11) and (21), the matrix form can be written:
Figure BDA0002752856330000091
the iteration method is 'calculating twice', and the solution is realized through iterative approximation between a true value and an estimated value:
Figure BDA0002752856330000092
in formula (42): v'kRepresenting the estimated value of the pipe section flow obtained by the k-th calculation; vkThe true value of the pipe section flow of the kth calculation is represented; vk+1Represents the true value of the pipe section flow calculated at the (k + 1) th time; omega is a weight coefficient;
the iteration termination condition is as follows:
|Vk+1-Vk|≤ε (43)
in formula (43): ε represents the calculation accuracy.
Advantageous effects
The medium-low pressure gas distribution pipe network modeling method of the hydrogen-doped natural gas fully considers the physical characteristic difference of the hydrogen-doped natural gas and the pure natural gas, establishes the medium-low pressure gas distribution pipe network modeling method of the hydrogen-doped natural gas, utilizes an iterative solution principle through three different methods of a node method, a ring network method and a ring energy method, contrasts and analyzes respective characteristics and operation results, fully considers the difference of physicochemical properties of the hydrogen-doped natural gas and the requirements of actual engineering, enables pipe network modeling to better meet the requirements of the actual engineering under the condition of ensuring simulation accuracy, and achieves the medium-low pressure gas distribution pipe network modeling of the hydrogen-doped natural gas which better meets the requirements of the actual engineering under the condition of ensuring the simulation accuracy.
Drawings
FIG. 1 is a flow chart of a method implemented in an embodiment of the invention.
Fig. 2 is a diagram of an 11-node hydrogen-loaded natural gas pipeline network in the embodiment of the invention.
Fig. 3 is a graph of a simulation result of a steady-state node method of an 11-node hydrogen-doped natural gas pipe network in the embodiment of the present invention.
Fig. 4 is a comparison graph of node pressure simulation results of three calculation methods in the embodiment of the present invention.
Fig. 5 is a comparison graph of pipe segment flow simulation results of three calculation methods in the embodiment of the present invention.
Detailed Description
The accompanying drawings, which are incorporated in and constitute a part of this specification, illustrate an embodiment of the invention and, together with the description, serve to explain the invention and not to limit the invention.
The simulated hydrogen-doped natural gas pipe network system in the embodiment of the invention is shown in figure 2, the system is a medium-pressure hydrogen-doped natural gas pipe network in a certain area, the system comprises 11 nodes, 12 sections of pipe networks, the gas source node is number 11 node, the reference pressure is 1.8MPa, the length, the inner diameter, the rated temperature and the relative density of fluid in the pipe of each pipe network are shown in table 1, and the load flow of each node is shown in table 2.
Table 1. data of hydrogen-blended natural gas pipe network in the calculation example
Serial number Head and tail ends Length/km Inner diameter/mm temperature/K Relative density
1 11-1 12.3 495 289.15 0.494
2 1-2 7.6 495 289.00 0.494
3 2-3 8.1 495 288.85 0.494
4 1-4 11.1 495 288.65 0.494
5 4-5 8.9 495 288.75 0.494
6 5-6 14.4 495 288.85 0.494
7 7-6 10.2 495 288.65 0.494
8 11-7 7.8 495 288.95 0.494
9 7-8 13.7 495 288.35 0.494
10 8-9 11.5 495 288.55 0.494
11 10-9 13.8 495 288.25 0.494
12 10-5 12.6 495 288.15 0.494
TABLE 2 node load parameters in the example
Serial number 1 2 3 4 5 6 7 8 9 10
Load flow/(m)3/h) 0 4500 4200 3630 4700 5500 2300 1150 3930 1600
The method for building the medium-low pressure gas distribution pipe network of the hydrogen-loaded natural gas, which is provided by the embodiment of the invention, is shown by referring to fig. 1, and comprises the following steps:
s1, establishing a steady-state isothermal model of the hydrogen-doped natural gas network pipeline, wherein the model mainly comprises 2 aspects: firstly, setting parameters of gas density, Reynolds number, dynamic viscosity and friction coefficient of the hydrogen-doped natural gas, and secondly, describing mathematical models of nodes, pipelines and loops of a hydrogen-doped natural gas network;
s2, establishing a solving method of the medium-low pressure gas distribution pipe network in the hydrogen-doped natural gas by adopting an iteration method based on the mathematical model description of the steady-state isothermal model of the network pipeline;
s3, according to the solving method, a node method, a ring network method and a ring energy method are adopted to establish a mathematical solving model of node pressure and pipeline flow of the network, the internal relation between the mathematical solving models obtained by comparing the three methods is analyzed, the characteristics and the application occasions of the mathematical solving models are compared and analyzed, and finally the modeling method of the medium-low distribution pipe network of the hydrogen-doped natural gas is realized.
In step S1, the specific process of setting the parameters of the density, reynolds number, kinematic viscosity, and molar group coefficient of the hydrogen-doped natural gas is as follows:
1) density of natural gas mixed with hydrogen
The density rho of the hydrogen-doped natural gas is relative density delta and air density rhoAir (a)The specific expression is as follows:
Figure BDA0002752856330000111
in formula (1): rho0The dry air density in the standard state (273.15K, 101.3kPa) was 1.293kg/m3(ii) a p represents the pressure intensity of the hydrogen-doped natural gas and has the unit of MPa; t represents the temperature of the hydrogen-doped natural gas and has the unit of K.
The natural gas density of the hydrogen-doped gas is
Figure BDA0002752856330000112
In the formula (2), the reaction mixture is,
Figure BDA0002752856330000113
represents CH4The density of (a) of (b),
Figure BDA0002752856330000114
represents CH4The volume of (a) to (b),
Figure BDA0002752856330000115
represents H2The density of (a) of (b),
Figure BDA0002752856330000116
represents H2The volume of (a).
Considering the transportation safety of pipeline facilities and related transformation engineering of end equipment, the hydrogen loading ratio of the hydrogen-loaded natural gas is often limited, and is usually limited to the range of 5-20%, so the hydrogen volume ratio sigma% is uniformly set to 10% by the method, namely:
Figure BDA0002752856330000117
the density of the natural gas is the product of the relative density delta and the air density, and then the relative density of the natural gas is as follows:
Figure BDA0002752856330000118
2) reynolds number Re
The reynolds number is a similar criterion number for characterizing viscous influence in fluid mechanics, and a smaller number means that the viscous force influence is more significant, and a larger number means that the inertia of the fluid has a larger influence on the fluid flow in the circular tube, and the reynolds number is characterized by the following formula:
Figure BDA0002752856330000121
in formula (4): gamma denotes the kinematic viscosity of the gas, m2S; μ denotes dynamic viscosity, N.s/m2(ii) a v represents the gas flow rate, m/s; v represents the gas volume flow, m3S; m represents gas mass flow, kg/s; d represents the inner diameter of the pipe, m.
3) Dynamic viscosity mu
The dynamic viscosity represents the internal friction force generated by the interaction of the fluids, can be obtained by the ratio of stress to strain rate, and can be obtained by a calculation expression of the dynamic viscosity according to the pressure, the temperature and the relative density of the hydrogen-doped natural gas under a standard state:
Figure BDA0002752856330000122
in formula (5): Δ represents the relative density of the hydrogen-doped natural gas under the standard state, namely the ratio of the air density under the standard state; rho is gas density, kg/m3(ii) a T represents the gas temperature, K.
4) Coefficient of friction resistance lambda
The friction coefficient of the hydrogen-doped natural gas pipeline is generally calculated by adopting a Colebrook equation, namely:
Figure BDA0002752856330000123
in formula (6): k represents the absolute roughness of the inner wall of the pipe, and is mm, the thickness of the steel pipe is generally 0.1-0.2 mm, and the thickness of the polyethylene pipe is generally 0.01 mm; d represents the inner diameter of the pipe, mm; re represents the Reynolds number.
The parameters of the simulated hydrogen-doped natural gas pipe network system in the example of the present invention were calculated according to the equations (1) to (6) as shown in table 3.
Table 3 data of the hydrogen-blended natural gas pipeline network in the calculation example
Figure BDA0002752856330000124
Figure BDA0002752856330000131
In step S1, the mathematical model description process of the nodes, pipes, and loops of the hydrogen-doped natural gas network includes:
for a hydrogen loaded natural gas network, two basic equations can first be constructed to describe the network state: a node-pipeline correlation matrix A and a loop-pipeline correlation matrix BETA.
Element a of the node-pipe incidence matrix AijThe method can be deduced according to a kirchhoff first law, and the construction condition is that the pressure of the reference node is known.
Figure BDA0002752856330000132
Figure BDA0002752856330000133
Element b of loop-pipeline association matrix BETAnjIt can be deduced from kirchhoff's second law, assuming that the reference flow direction of the loop is clockwise.
Figure BDA0002752856330000141
Figure BDA0002752856330000142
The two matrixes logically describe the relationship among the nodes, the pipelines and the loops, not only show the relationship among the pipelines and the loops, but also realize the mutual conversion among the node pressure, the pipeline pressure drop and the loop pressure drop, and by analogy with a node voltage solution and a loop current solution in a power system, the node-pipeline correlation matrix A and the loop-pipeline correlation matrix BETA can be considered to be matrixes describing the state of a pipe network. According to the graph theory, there is an orthogonal relationship between node-pipeline correlation matrix a and loop-pipeline correlation matrix beta, i.e. the inner product of the row vectors of node-pipeline correlation matrix a and the row vectors of loop-pipeline correlation matrix beta is always equal to zero:
A×BT＝B×AT＝0 (9)
for the hydrogen-doped natural gas network, assuming that the number of nodes is i (including reference nodes), the number of pipelines is j, and the number of loops is n, the relationship between the three can be obtained as follows:
j＝n+i-1 (10)
for a loaded natural gas network, the following equations can be used for description.
Obtaining a node flow equation according to a kirchhoff first law, wherein the flow of the inflow node is equal to the flow of the outflow node (the pressure of a reference node is known):
A(i-1)×jVj×1＝q(i-1)×1 (11)
in the formula (11), q(i-1)×1Representing the traffic of the egress node, A(i-1)×jRepresenting a loop-pipe correlation matrix, V, with reference nodes removedj×1Representing the flow vector of each pipe section.
The loop pressure drop equation is derived from kirchhoff's second law, with a pressure drop along any one closed loop of 0:
Bn×jΔpj×1＝0n×1 (12)
in the formula (12), Bn×jRepresenting a loop-pipe correlation matrix with n loops and j pipes, Δ pj×1Representing the pressure difference between the upper end and the lower end of each pipe section.
According to the pressure drop equation of the pipeline, the pressure drop of the pipeline is related to the pressure difference of the upper end and the lower end of the pipeline:
Figure BDA0002752856330000151
in formula (13): p is a radical of(i-1)×1Expressed as the pressure of each node relative to a reference node.
Due to A x BT＝B×AT0, so
Figure BDA0002752856330000154
Therefore, the formula can be expressed by one of the above formulae.
When the pipe diameter is known, each pipeline has two unknowns of pressure drop and flow, namely 2j unknowns in total, and the listed equation number is j + (i-1) + n-2 j, so that the solution can be realized.
Secondly, the corresponding quantity relation exists between the pipeline pressure drop and the pipeline flow, namely the pipeline flow pressure drop equation:
Δp＝G|V|V (14)
in formula (14): g represents a diagonal element matrix of the flow resistance coefficient of the pipeline, namely the flow resistance coefficient of each pipeline section is on the corresponding diagonal line, and the rest elements of the matrix are 0; v represents the calculated flow vector of the pipeline in m3/h。
For medium and low pressure pipe networks (low pressure: 0-0.01 MPa; medium pressure pipeline: 0.01-0.4 MPa), the following are adopted in a unified way:
Figure BDA0002752856330000152
in formula (15): Δ pabRepresenting the pressure drop, kPa, of the pipe connecting nodes a and b2；paAnd pbIndicating the pressure at nodes a and b, kPa; lambda represents the hydraulic friction coefficient of the pipeline and is related to the Reynolds number Re, the relative density delta of the hydrogen-doped natural gas and the like; znRepresenting the compression factor of the gas under the standard condition, and taking 1; d represents the inner diameter of the pipeline, mm; rho represents the density of the natural gas doped with hydrogen, kg/m3；TnRepresents the standard state absolute temperature, 273.15K; t represents the pipe temperature, K; l represents the length of the pipeline, km;
table 4. coefficient of flow resistance of natural gas pipeline in the calculation
Sequence number of pipe sections 1 2 3 4 5 6 7 8 9 10 11 12
Coefficient of flow resistance 0.4051 1.1007 2.1351 1.1090 3.5608 1.7439 0.8519 0.5530 0.8568 1.2508 4.6096 11.8671
Z represents the compression factor of the hydrogen-doped natural gas, and the formula is as follows:
Figure BDA0002752856330000153
in formula (16): t iscRepresents the critical temperature of the gas, K; p is a radical ofcRepresents the critical pressure of the gas, Pa;
TABLE 5 compression factor of the pipe network of the exemplary embodiment of the Natural gas Loading
Sequence number of pipe sections 1 2 3 4 5 6 7 8 9 10 11 12
Compression factor 0.9722 0.9722 0.9721 0.9720 0.9721 0.9721 0.9720 0.9721 0.9719 0.9720 0.9718 0.9718
pavRepresents the average pressure of the pipe mn, Pa, expressed as:
Figure BDA0002752856330000161
in the formula (17), pmRepresenting the nodal pressure, p, of the beginning m of the pipe sectionnRepresenting the node pressure at the end n of the pipe section.
Let element S in diagonal matrix S of j x jjjSatisfies the following conditions:
Figure BDA0002752856330000162
then, according to the equations (11) to (14) and (18), the following relationships can be obtained:
Figure BDA0002752856330000163
ASATp＝q (20)
B(G|V|V)＝0 (21)
therefore, the node method, the ring network method and the ring energy method can be respectively applied to solving.
Further, in step S2, the process of establishing the solution method for the medium-low pressure distribution pipe network of the hydrogen-doped natural gas is as follows:
the common solution method is to introduce the amount of unbalance in iteration, perform taylor series expansion, and limit the amount of unbalance below a specified value, which can be regarded as an accurate solution.
Suppose there are N variables x1,x2,L,xNThe equation is known:
Fi(x1,x2,L,xN)＝0,i＝1,2,L,N (22)
let variable matrix X ═ X1 x2 … xi]', unbalance matrix δ X ═ δ X1 δx2 … δxi]', function matrix F ═ F1 F2 … Fi]', square of unbalance amountMatrix delta X2＝[(δx1)2 (δx2)2 … (δxi)2]', then each FiThe Taylor series expansion is adopted to obtain:
Figure BDA0002752856330000171
let the partial derivatives form a Jacobian matrix J, where element J is number (i, J)ijComprises the following steps:
Figure BDA0002752856330000172
then:
F(X+δX)＝F(X)+J×δX+O(δX2) (25)
neglecting O (δ X)2) And let F (X + δ X) ═ 0, a linear equation system for the unbalance amount δ X can be obtained:
J×δX＝-F(X) (26)
in formulae (25) and (26): δ X may let each function approach zero at the same time.
Finally, the unbalance amount is added into the solution vector to continuously correct the initial value.
Xk+1＝Xk+δXk (27)
In formula (27): xk+1Represents the result of the iteration of the (k + 1) th solution vector, XkRepresents the result of the k-th iteration of the solution vector, deltaXkCorrection values representing the solution vector for the k-th iteration.
The key to the solution is to construct an invertible Jacobian square matrix.
Further, in step S3, the specific process of establishing the mathematical solution model of the node pressure and the pipeline flow of the network by using the node method is as follows:
according to equation (20), it can be transformed into:
ASATp-q＝0 (28)
equation (28) is written as:
F(p)＝ASATp-q＝0 (29)
and (3) solving by using a node method, firstly providing an initial estimation value of the pressure of each node, and secondly continuously correcting the estimation value until a final result is obtained. Because the node pressure value is only an estimated value, and F (p) is actually not equal to zero in the iteration process, the unbalance amount of the node pressure can be obtained, and the node pressure value can be regarded as an error correction amount:
Jp×δpk＝-F(pk) (30)
in formula (30): j. the design is a squarepRepresenting a Jacobian matrix formed by partial derivation of the node pressure p; p is a radical ofkAnd representing the node pressure value solved by the kth iteration.
Figure BDA0002752856330000181
The iterative process for correcting the node pressure estimate is:
pk+1＝pk+δpk (32)
in the formula (32), pk+1Representing the node pressure value of the (k + 1) th iteration solution; p is a radical ofkRepresenting the node pressure value of the kth iterative solution; δ pkAnd representing the node pressure correction value solved by the kth iteration.
The iteration termination condition is as follows:
|δpk|≤ε (33)
in formula (33): ε is the calculation accuracy and is small enough.
In each iteration, the value of S varies according to equations (13), (14) and (18), and the sign of V is continuously corrected according to the sign of Δ p.
Further, in step S3, the specific process of establishing the mathematical solution model of the node pressure and the pipeline flow of the network by using the ring network method is as follows:
according to equation (21), the equation can be written as:
F(V)＝BG|V|V＝0 (34)
the ring network method needs to define a group of closed loops in a pipe network, give an initial estimation value of branch flow of each pipeline and ensure the balance of the flow at each node. Given each pipe branch flow, a loop flow also needs to be introduced. The loop flow is a correction e added to the branch flow to correct the branch flow estimate, in the same direction as the loop reference direction.
A particular natural gas pipeline may be subject to more than one loop flow rate, and in general, the pipeline flow rate is a function of the initial estimate and all loop flow rates, i.e., equation (34) can be converted to:
BG|V+BTe|(V+BTe)＝0 (35)
in the equation (35), e represents a correction amount added to the branch flow rate to correct the branch flow rate estimate;
then:
F(e)＝BG|V+BTe|(V+BTe)＝0 (36)
first an initial value of the loop flow is given and then iterations continue until a true solution is reached. In the iterative solution process of the ring network method, F (e) is not equal to zero, namely, a loop flow error is introduced into each loop. The unbalance of the loop flow can be obtained according to the above, and can be regarded as an error correction amount:
Je×δek＝-F(ek) (37)
in formula (37): j. the design is a squareeRepresenting a Jacobian matrix formed by performing partial derivation on the flow correction value e of the branch circuit; e.g. of the typekRepresenting a branch flow correction value of the kth iterative solution; delta ekAnd representing the corrected value of the branch flow correction value solved by the kth iteration.
Figure BDA0002752856330000191
The iterative process for correcting the loop flow estimate is:
ek+1＝ek+δek (39)
the iteration termination condition is as follows:
|δek|≤ε (40)
in formula (40): ε is the calculation accuracy and is small enough.
In the process JeThe binding e is changed continuously.
The specific process of establishing the mathematical solving model of the node pressure and the pipeline flow of the network by adopting the ring energy method is as follows:
according to equations (11) and (21), the matrix form can be written:
Figure BDA0002752856330000201
the iteration method is 'calculating twice', and the solution is realized through iterative approximation between a true value and an estimated value:
Figure BDA0002752856330000202
in formula (42): v'kRepresenting the estimated value of the pipe section flow obtained by the k-th calculation; vkThe true value of the pipe section flow of the kth calculation is represented; vk+1Represents the true value of the pipe section flow calculated at the (k + 1) th time; ω is a weight coefficient, and ω is generally 0.5.
The iteration termination condition is as follows:
|Vk+1-Vk|≤ε (43)
in the formula: ε is the calculation accuracy and is small enough.
As shown in fig. 3, a graph of the simulation result of the steady-state node method of the 11-node hydrogen-doped natural gas pipe network is shown. 4-5, which are graphs comparing the results of the node pressure and pipe section flow simulation for the three methods.
The initial value of the node method is the node pressure, the requirement on the initial value is not high, the equation number is 10 nodes with unknown pressure, the iteration times are the most, the iteration times are 49 times, and the convergence speed is the slowest. Because the pipeline flow is decided by the pressure difference, when the pipe diameter of the pipeline is large, the pressure difference is small and the pipeline has different sensitivity to the influence of the node pressure, the precision requirement of the pipeline flow can not be ensured, thereby influencing the convergence speed. According to the simulation result, the calculation precision and complexity of the node method are between those of the ring network method and the ring energy method. The initial value of the ring network method is the pipeline flow, and the requirement is met by a node flow balance equation. The number of equations is minimum, equal to 2 network loops, the iteration times are minimum, 4 times, and the solving speed is fastest. However, the algorithm has extremely high requirements on the selection of the initial pipeline flow, the selection of the initial value directly influences the solving result, and the error fluctuation of the result is large. The initial value of the ring energy method is the pipeline flow, the requirement on the initial value is not high, the equation number is the sum of the node number with unknown network pressure and the loop number, namely the number of the pipelines is 12, and the iteration number is 45. The loop energy method needs to consider the node flow and the loop pressure drop equation at the same time, the obtained most accurate result is obtained, and the result of the algorithm is generally taken as a standard value in literature research. However, the number of the equations of the ring energy method is increased by 2 compared with the node method and 10 compared with the ring network method, and if large-scale urban complex natural gas distribution network load flow calculation is involved, the calculation complexity is greatly improved.
As shown in tables 6 and 7, the simulation results of the node pressure and the pipe flow of the three algorithms are compared with the algorithm complexity and the error.
TABLE 6 comparison of simulation results of three algorithms
Figure BDA0002752856330000211
TABLE 7 comparison of complexity and error for three algorithms
Figure BDA0002752856330000212
Figure BDA0002752856330000221
As can be seen from the data in table 7, the number of equations of the ring network method is the minimum, the number of iterations is the minimum, but the simulation result is the worst, and there is a larger relative error compared with the ring network method and the node method. The ring energy method and the ring network method directly carry out iterative analysis on the flow, and can directly identify the correctness of the initial pipeline flow direction. Although the node method cannot directly solve the pipeline flow, the correctness of the initial pipeline flow direction can be well identified through the sign of delta p. On the calculation precision, the results of the node method and the ring energy method have no obvious error, and higher precision is obtained. Meanwhile, for the operation control platform of the gas distribution network, node pressure monitoring points are generally dense and accurate, the number of pipeline flow monitoring points is small, and the accuracy is low, so that the method for carrying out load flow calculation by adopting a node method is more advantageous. From the analysis of application occasions, the node pressure intensity is generally monitored by a gas distribution pipe network control platform, so that the node method has more practical application advantages.
It should be noted that the various features of the embodiments described in the above embodiments, as well as other nodes added to the system, may be combined in any suitable manner without contradiction. In order to avoid unnecessary repetition, other possible combinations of the present invention will not be described.

Claims (7)

1. A medium and low distribution pipe network modeling method of hydrogen-doped natural gas is characterized by comprising the following steps:
s1, establishing a steady-state isothermal model of the hydrogen-doped natural gas network pipeline, wherein the steady-state isothermal model comprises parameter settings of density, Reynolds number, dynamic viscosity and friction coefficient of the hydrogen-doped natural gas and mathematical model descriptions of nodes, pipelines and loops of the hydrogen-doped natural gas network;
s2, based on the mathematical model description in the step S1, establishing a solving method of the medium-low pressure gas distribution pipe network of the hydrogen-doped natural gas by adopting an iteration method;
s3, according to the solving method in the step S2, mathematical solving models of node pressure and pipeline flow of the network are respectively established by adopting a node method, a ring network method and a ring energy method, internal relations among the mathematical solving models obtained by comparing the three methods are analyzed, characteristics and application occasions of the mathematical solving models are compared and analyzed, and finally the modeling method of the medium-low distribution gas pipeline network of the hydrogen-doped natural gas is realized.
2. The modeling method for medium-low gas distribution pipe networks of hydrogen-loaded natural gas according to claim 1, wherein in step S1, the specific process for setting the parameters of density, reynolds number, kinematic viscosity and molar group coefficient of hydrogen-loaded natural gas is as follows:
the density rho of the hydrogen-doped natural gas is relative density delta and air density rhoAir (a)The specific expression is as follows:
Figure FDA0002752856320000011
in formula (1): rho0The dry air density in the standard state (273.15K, 101.3kPa) was 1.293kg/m3(ii) a p represents the pressure intensity of the hydrogen-doped natural gas and has the unit of MPa; t represents the temperature of the hydrogen-doped natural gas and has the unit of K;
the density ρ of the hydrogen-doped natural gas can be obtained by the formula (2), that is:
Figure FDA0002752856320000012
in the formula (2), the reaction mixture is,
Figure FDA0002752856320000013
represents CH4The density of (a) of (b),
Figure FDA0002752856320000014
represents CH4The volume of (a) to (b),
Figure FDA0002752856320000015
represents H2The density of (a) of (b),
Figure FDA0002752856320000016
represents H2The volume of (a);
the hydrogen volume ratio σ% was set to 10%, that is:
Figure FDA0002752856320000017
the characterization formula of the Reynolds number Re of the hydrogen-doped natural gas is as follows:
Figure FDA0002752856320000018
in formula (4): gamma denotes the kinematic viscosity of the gas in m2S; mu is dynamic viscosity in the unit of N.s/m2(ii) a v represents the gas flow rate in m/s; v represents the gas volume flow rate in m3S; m represents the gas mass flow rate, and the unit is kg/s; d represents the inner diameter of the pipeline and has the unit of m;
the calculation expression for the kinematic viscosity μ is as follows:
Figure FDA0002752856320000021
in formula (5): Δ represents the relative density of the hydrogen-doped natural gas under the standard state, namely the ratio of the air density under the standard state; rho is the gas density in kg/m3(ii) a T represents the gas temperature in K;
the friction coefficient lambda of the hydrogen-doped natural gas is calculated by adopting a Colebrook equation, namely:
Figure FDA0002752856320000022
in formula (6): k represents the absolute roughness of the inner wall of the tube in mm; d represents the inner diameter of the pipe, mm; re represents the Reynolds number of the hydrogen-loaded natural gas.
3. The modeling method for medium-low distribution networks of hydrogen-loaded natural gas according to claim 2, wherein in step S1, the mathematical model description of the nodes, pipes and loops of the hydrogen-loaded natural gas network comprises:
two basic equations are constructed to describe the network state: a node-pipeline correlation matrix A and a loop-pipeline correlation matrix BETA;
element a of the node-pipe incidence matrix AijDeducing according to a kirchhoff first law, wherein the construction condition is that the pressure of a reference node is known;
Figure FDA0002752856320000023
element b of said loop-pipeline association matrix BETAnjDeducing according to a kirchhoff second law, and assuming that the reference flow direction of the loop is clockwise;
Figure FDA0002752856320000031
simulating a node voltage solution and a loop current solution in the power system, wherein the node-pipeline correlation matrix A and the loop-pipeline correlation matrix BETA are both matrices describing the state of the pipe network; according to the graph theory, there is an orthogonal relationship between node-pipeline correlation matrix a and loop-pipeline correlation matrix beta, that is, the inner product of the row vector of node-pipeline correlation matrix a and the row vector of loop-pipeline correlation matrix beta is always equal to zero, and the specific expression is as follows:
A×BT＝B×AT＝0 (9)
for the hydrogen-doped natural gas network, assuming that the number of nodes is i (including reference nodes), the number of pipelines is j, and the number of loops is n, the corresponding relationship among the number of nodes i, the number of pipelines j, and the number of loops n can be obtained as follows:
j＝n+i-1 (10)
aiming at a hydrogen-doped natural gas network, the following equation is used for description;
obtaining a node flow equation according to a kirchhoff first law, wherein the flow of the inflow node is equal to the flow of the outflow node, and the pressure of a reference node is known:
A(i-1)×jVj×1＝q(i-1)×1 (11)
in the formula (11), q(i-1)×1Representing the amount of traffic that flows out of the node,A(i-1)×jrepresenting a loop-pipe correlation matrix, V, with reference nodes removedj×1Representing the flow vector of each pipe section;
the loop pressure drop equation is derived from kirchhoff's second law, with a pressure drop along any one closed loop of 0:
Bn×jΔpj×1＝0n×1 (12)
in the formula (12), Bn×jRepresenting a loop-pipe correlation matrix with n loops and j pipes, Δ pj×1Representing the pressure difference between the upper end and the lower end of each pipe section;
according to the pressure drop equation of the pipeline, the relation between the pressure drop of the pipeline and the pressure difference of the upper end and the lower end of the pipeline is as follows:
Figure FDA0002752856320000032
in formula (13): p is a radical of(i-1)×1Expressed as the pressure of each node relative to a reference node;
due to A x BT＝B×AT0, so
Figure FDA0002752856320000041
Therefore, the formula can be expressed by one of the two formulas;
when the pipe diameter is known, each pipeline has two unknowns of pressure drop delta p and flow, namely 2j unknowns, and the listed equation number is j + (i-1) + n ═ 2j for solving;
secondly, according to the corresponding quantitative relation between the pipeline pressure drop delta p and the pipeline flow, the equation of the pipeline flow pressure drop is obtained:
Δp＝G|V|·V (14)
in formula (14): g represents a diagonal element matrix of the flow resistance coefficient of the pipeline, namely the flow resistance coefficient of each pipeline section is on the corresponding diagonal line, and the rest elements of the matrix are 0; v represents the calculated flow vector of the pipeline in m3/h；
For medium-low pressure pipe network, the pipeline pressure drop delta p between the connection node a and the connection node babMeter (2)The calculation formula and the calculation formula of the pipeline flow resistance coefficient G are as follows:
Figure FDA0002752856320000042
in formula (15): Δ pabRepresents the pressure drop of the pipe connecting nodes a and b in kPa2；paAnd pbRepresenting the pressure at nodes a and b in kPa; lambda represents the hydraulic friction coefficient of the pipeline; znRepresenting the compression factor of the gas under the standard condition, and taking 1; d represents the inner diameter of the pipeline and has the unit of mm; rho represents the density of the natural gas doped with hydrogen and has the unit of kg/m3；TnThe standard state absolute temperature is represented, namely 273.15K; t represents the temperature of the pipeline and has the unit of K; l represents the length of the pipeline, and the unit is km; z represents the compression factor of the hydrogen-doped natural gas, and the expression of Z is as follows:
Figure FDA0002752856320000043
in formula (16): t iscRepresents the critical temperature of the gas, in K; p is a radical ofcRepresents the critical pressure of the gas in Pa; p is a radical ofavDenotes the mean pressure of the pipe mn in Pa, pavThe expression is as follows:
Figure FDA0002752856320000044
in the formula (17), pmRepresenting the nodal pressure, p, of the beginning m of the pipe sectionnRepresenting the node pressure at the end n of the pipe section;
let element S in diagonal matrix S of j x jjjSatisfies the following conditions:
Figure FDA0002752856320000051
from the equation (11) -and, the following relationship is obtained:
Figure FDA0002752856320000052
ASATp＝q (20)
B(G|V|V)＝0 (21)
4. the modeling method for the medium-low pressure distribution pipe network of the hydrogen-loaded natural gas according to claim 3, wherein in the step S2, the solving method for establishing the medium-low pressure distribution pipe network of the hydrogen-loaded natural gas by using the iteration method specifically comprises the following steps:
suppose there are N variables x1,x2,L,xNThe equation is known:
Fi(x1,x2,L,xN)＝0,i＝1,2,L,N (22)
let variable matrix X ═ X1 x2…xi]', unbalance matrix δ X ═ δ X1 δx2…δxi]', function matrix F ═ F1F2…Fi]', the square matrix of unbalance amount deltaX2＝[(δx1)2 (δx2)2…(δxi)2]', then each FiThe Taylor series expansion is adopted to obtain:
Figure FDA0002752856320000053
let the partial derivatives form a Jacobian matrix J, where element J is number (i, J)ijIs composed of
Figure FDA0002752856320000054
Then:
F(X+δX)＝F(X)+J×δX+O(δX2) (25)
neglecting O (δ X)2) And make F (X +δ X) is 0, a linear equation set for the unbalance amount δ X can be obtained:
J×δX＝-F(X) (26)
in formulae (25) and (26): deltax lets each function approach zero at the same time,
the amount of unbalance is added to the solution vector to continuously correct the initial value,
Xk+1＝Xk+δXk (27)
in formula (27): xk+1Represents the result of the iteration of the (k + 1) th solution vector, XkRepresents the result of the k-th iteration of the solution vector, deltaXkCorrection values representing the solution vector for the k-th iteration.
5. The modeling method for the medium-low gas distribution pipe network of the hydrogen-loaded natural gas as claimed in claim 3, wherein in the step S3, the concrete process of building the mathematical solution model of the node pressure and the pipeline flow of the network by using the node method is as follows:
transform equation (20) into:
ASATp-q＝0 (28)
the equation is written as:
F(p)＝ASATp-q＝0 (29)
firstly, giving an initial estimation value of each node pressure, and secondly, obtaining the unbalance of the node pressure as an error correction quantity by continuously correcting the estimation value:
Jp×δpk＝-F(pk) (30)
in formula (30): j. the design is a squarepRepresenting a Jacobian matrix formed by partial derivation of the node pressure p; p is a radical ofkRepresenting the node pressure value of the kth iterative solution;
Figure FDA0002752856320000061
the iterative process for correcting the node pressure estimate is:
pk+1＝pk+δpk (32)
in the formula (32), pk+1Representing the node pressure value of the (k + 1) th iteration solution; p is a radical ofkRepresenting the node pressure value of the kth iterative solution; δ pkRepresenting the node pressure correction value of the kth iterative solution;
the iteration termination condition is as follows:
|δpk|≤ε (33)
in formula (33): ε represents the calculation accuracy.
In each iteration, the value of S varies according to equations (13), (14) and (V), and the sign of the desired V needs to be continuously corrected according to the sign of Δ p.
6. The method for solving the medium-low pressure distribution pipe network in the hydrogen-loaded natural gas as claimed in claim 3, wherein in the step S3, the specific process for establishing the mathematical solution model of the node pressure and the pipeline flow of the network by using the ring network method is as follows:
according to the formula, the formula can be written as:
F(V)＝BG|V|V＝0 (34)
the formula can be converted into:
BG|V+BTe|(V+BTe)＝0 (35)
in the equation (35), e represents a correction amount added to the branch flow rate to correct the branch flow rate estimate;
then:
F(e)＝BG|V+BTe|(V+BTe)＝0 (36)
firstly, giving an initial value of loop flow, and then continuously iterating until a real solution is reached; in the iterative solution process of the ring network method, f (e) is not equal to zero, that is, a loop flow error is introduced into each loop, and accordingly, the unbalance amount of the loop flow can be obtained, namely, the unbalance amount can be regarded as an error correction amount:
Je×δek＝-F(ek) (37)
in formula (37): j. the design is a squareeIndicating the flow of the opposite branchCalculating a Jacobian matrix formed by partial derivatives according to the correction value e; e.g. of the typekRepresenting a branch flow correction value of the kth iterative solution; delta ekRepresenting the corrected value of the branch flow correction value of the kth iterative solution;
Figure FDA0002752856320000081
the iterative process for correcting the loop flow estimate is:
ek+1＝ek+δek (39)
the iteration termination condition is as follows:
|δek|≤ε (40)
in formula (40): epsilon is the calculation precision;
in the process JeThe binding e is changed continuously.
7. The method for solving the medium-low pressure distribution pipe network in the hydrogen-loaded natural gas according to claim 3, wherein the specific process for establishing the mathematical solution model of the node pressure and the pipeline flow of the network by adopting the ring energy method is as follows:
according to the formula and the formula, the matrix form can be written:
Figure FDA0002752856320000082
the iteration method is 'calculating twice', and the solution is realized through iterative approximation between a true value and an estimated value:
Figure FDA0002752856320000083
in formula (42): v'kRepresenting the estimated value of the pipe section flow obtained by the k-th calculation; vkThe true value of the pipe section flow of the kth calculation is represented; vk+1True representing the (k + 1) th calculated spool piece flowA value; omega is a weight coefficient;
the iteration termination condition is as follows:
|Vk+1-Vk|≤ε (43)
in formula (43): ε represents the calculation accuracy.
CN202011191416.8A 2020-10-30 2020-10-30 Modeling method for medium-low pressure gas distribution pipe network of hydrogen-doped natural gas Pending CN112257355A (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202011191416.8A CN112257355A (en) 2020-10-30 2020-10-30 Modeling method for medium-low pressure gas distribution pipe network of hydrogen-doped natural gas

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202011191416.8A CN112257355A (en) 2020-10-30 2020-10-30 Modeling method for medium-low pressure gas distribution pipe network of hydrogen-doped natural gas

Publications (1)

Publication Number Publication Date
CN112257355A true CN112257355A (en) 2021-01-22

Family

ID=74268286

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202011191416.8A Pending CN112257355A (en) 2020-10-30 2020-10-30 Modeling method for medium-low pressure gas distribution pipe network of hydrogen-doped natural gas

Country Status (1)

Country Link
CN (1) CN112257355A (en)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113987943A (en) * 2021-10-29 2022-01-28 国家石油天然气管网集团有限公司 Natural gas pipe network gas calorific capacity prediction method suitable for mixed gas source
CN114386206A (en) * 2022-01-14 2022-04-22 国网江苏省电力有限公司经济技术研究院 Modeling method of hydrogen-mixed natural gas network optimization scheduling model

Citations (13)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101333961A (en) * 2008-08-07 2008-12-31 清华大学 Hydrogen gas natural gas mixed fuel engine optimizing method
US20150261893A1 (en) * 2014-04-22 2015-09-17 China University Of Petroleum - Beijing Method and apparatus for determining pipeline flow status parameter of natural gas pipeline network
CN106777708A (en) * 2016-12-21 2017-05-31 天津大学 Steady state analysis method of electric power-natural gas regional comprehensive energy system
CN107291990A (en) * 2017-05-24 2017-10-24 河海大学 Energy stream emulation mode based on electrical interconnection integrated energy system transient Model
US9897260B1 (en) * 2017-04-18 2018-02-20 Air Products And Chemicals, Inc. Control system in an industrial gas pipeline network to satisfy energy consumption constraints at production plants
US9915399B1 (en) * 2017-04-18 2018-03-13 Air Products And Chemicals, Inc. Control system in a gas pipeline network to satisfy demand constraints
CN109697308A (en) * 2018-11-30 2019-04-30 天津大学 A kind of natural gas transmission systematic steady state modeling method considering pipe network time-delay characteristics
CN110717643A (en) * 2019-08-07 2020-01-21 重庆大学 Natural gas network gas storage configuration method based on global sensitivity analysis
CN110765622A (en) * 2019-10-28 2020-02-07 南方电网科学研究院有限责任公司 Energy flow obtaining system, equipment and medium of natural gas pipeline model
CN110796295A (en) * 2019-10-15 2020-02-14 西安交通大学 Energy internet air network transmission optimization method
CN111046594A (en) * 2020-01-09 2020-04-21 东南大学 Hot water heating network dynamic simulation method based on cross iteration principle
CN111259547A (en) * 2020-01-16 2020-06-09 清华大学 Natural gas path modeling method for operation control of comprehensive energy system
CN111563315A (en) * 2020-04-08 2020-08-21 重庆大学 Topological analysis-based steady-state energy flow calculation method for electricity-gas comprehensive energy system

Patent Citations (13)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101333961A (en) * 2008-08-07 2008-12-31 清华大学 Hydrogen gas natural gas mixed fuel engine optimizing method
US20150261893A1 (en) * 2014-04-22 2015-09-17 China University Of Petroleum - Beijing Method and apparatus for determining pipeline flow status parameter of natural gas pipeline network
CN106777708A (en) * 2016-12-21 2017-05-31 天津大学 Steady state analysis method of electric power-natural gas regional comprehensive energy system
US9897260B1 (en) * 2017-04-18 2018-02-20 Air Products And Chemicals, Inc. Control system in an industrial gas pipeline network to satisfy energy consumption constraints at production plants
US9915399B1 (en) * 2017-04-18 2018-03-13 Air Products And Chemicals, Inc. Control system in a gas pipeline network to satisfy demand constraints
CN107291990A (en) * 2017-05-24 2017-10-24 河海大学 Energy stream emulation mode based on electrical interconnection integrated energy system transient Model
CN109697308A (en) * 2018-11-30 2019-04-30 天津大学 A kind of natural gas transmission systematic steady state modeling method considering pipe network time-delay characteristics
CN110717643A (en) * 2019-08-07 2020-01-21 重庆大学 Natural gas network gas storage configuration method based on global sensitivity analysis
CN110796295A (en) * 2019-10-15 2020-02-14 西安交通大学 Energy internet air network transmission optimization method
CN110765622A (en) * 2019-10-28 2020-02-07 南方电网科学研究院有限责任公司 Energy flow obtaining system, equipment and medium of natural gas pipeline model
CN111046594A (en) * 2020-01-09 2020-04-21 东南大学 Hot water heating network dynamic simulation method based on cross iteration principle
CN111259547A (en) * 2020-01-16 2020-06-09 清华大学 Natural gas path modeling method for operation control of comprehensive energy system
CN111563315A (en) * 2020-04-08 2020-08-21 重庆大学 Topological analysis-based steady-state energy flow calculation method for electricity-gas comprehensive energy system

Non-Patent Citations (6)

* Cited by examiner, † Cited by third party
Title
GIULIO GUANDALINI，等: "Dynamic modeling of natural gas quality within transport pipelines in presence of hydrogen injections", APPLIED ENERGY, no. 185, 31 December 2017 (2017-12-31), pages 1712 - 1723 *
SUYANG ZHOU: "Design and Evaluation of Operational Scheduling Approaches for HCNG Penetrated Integrated Energy System", VOLUME, 26 June 2019 (2019-06-26), pages 87792 - 87801 *
吕勃蓬，等: "天然气管网稳态分析方法研究现状与展望", 油气田地面工程, no. 11, 1 November 2012 (2012-11-01), pages 20 - 22 *
吴嫦: "天然气掺混氢气使用的可行性研究", 中国硕士论文全文数据库, no. 4, 15 April 2019 (2019-04-15), pages 038 - 1454 *
田贯三，等: "城镇天然气管网水力分析数学模型与计算方法", 天然气工业, no. 03, 28 May 2002 (2002-05-28), pages 96 - 98 *
黄明，等: "利用天然气管道掺混输送氢气的可行性分析", 煤气与热力, no. 04, 15 April 2013 (2013-04-15), pages 39 - 42 *

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113987943A (en) * 2021-10-29 2022-01-28 国家石油天然气管网集团有限公司 Natural gas pipe network gas calorific capacity prediction method suitable for mixed gas source
CN114386206A (en) * 2022-01-14 2022-04-22 国网江苏省电力有限公司经济技术研究院 Modeling method of hydrogen-mixed natural gas network optimization scheduling model
CN114386206B (en) * 2022-01-14 2023-11-21 国网江苏省电力有限公司经济技术研究院 Modeling method of hydrogen-mixed natural gas network optimization scheduling model

Similar Documents

Publication Publication Date Title
CN109344436B (en) Online simulation method for large complex natural gas pipe network system
CN112257355A (en) Modeling method for medium-low pressure gas distribution pipe network of hydrogen-doped natural gas
CN111563315B (en) Topology analysis-based steady-state energy flow calculation method for electric-gas comprehensive energy system
CN109800968B (en) Power-gas interconnection system probability energy flow analysis method considering natural gas system thermodynamic process
CN109726483A (en) A kind of radial heat supply network model of electric heating interconnection integrated energy system and its system
CN110866213A (en) Multi-network steady-state energy flow analysis method and device for electricity-gas integrated energy system
CN112069692A (en) Optimization solving method for natural gas pipe network transmission difference calculation
CN110765622A (en) Energy flow obtaining system, equipment and medium of natural gas pipeline model
Sheng et al. Two-stage state estimation approach for combined heat and electric networks considering the dynamic property of pipelines
CN110728032A (en) Quick power flow calculation method for electricity-heat interconnection comprehensive energy system considering ring network
CN108389136B (en) Natural gas probability-fuzzy energy flow analysis method considering multiple uncertainties
CN110532642A (en) A kind of calculation method that integrated energy system probability can flow
CN111414675A (en) Double-layer robust state estimation method and system for electric heating comprehensive energy system
CN115688617A (en) Method, system and equipment for calculating energy flow of hydrogen-doped natural gas pipeline and application of method, system and equipment
Ma et al. A novel analytical unified energy flow calculation method for integrated energy systems based on holomorphic embedding
CN110796295B (en) Energy Internet air network transmission optimization method
CN111783309A (en) Dynamic simulation method of steam heating network based on internal conservation
CN108964143B (en) Natural gas network static equivalent model of electricity-gas comprehensive energy system
CN110717643A (en) Natural gas network gas storage configuration method based on global sensitivity analysis
Xianxi et al. Modeling and simulation of steam pipeline network with multiple supply sources in iron& steel plants
CN114117819B (en) Steady-state simulation method for hot steam network
CN114880817A (en) Gas transmission pipeline dynamic hydraulic prediction method and device, electronic equipment and storage medium
CN109002419B (en) Dynamic analysis method and device for natural gas pipe network
CN105160062A (en) Hydraulic check method for reverse return pipe network
Song et al. Dynamic analysis and modeling of the natural gas pipeline using the electrical analogy

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