Simulation method of gas and heat subnetworks of integrated energy system based on variable step-size characteristic line method

Abstract

A kind of comprehensive energy system gas, heat subnetwork simulation method based on variable step characteristic line method, including based on the network structure of gas network, heat network etc. in comprehensive energy system subnet, construct the network node file and pipe branch file of each subnet;Determine minimum time step as basic step;According to network node file and pipe branch file, construct network correlation matrix, describe network topology relationship;Using characteristic line method realizes the preprocessing to the branch of time step greater than basic step;Under each basic step, in combination with network correlation matrix and characteristic line method, realize the state quantity solution to each node;When simulation time interval reaches the simulation time step of a branch, then the intermediate node state quantity of the branch is updated using the characteristic line method, otherwise the branch is not calculated.The present application considers the length difference of pipe in gas network and heat network, greatly reduces the calculation amount in network simulation process under the premise of ensuring the simulation accuracy of each branch, and realizes high-performance simulation.

Description

Technical Field

[0001] This invention pertains to the simulation of integrated energy systems, and in particular to a simulation method for the gas and heat subnets of integrated energy systems based on the variable step size characteristic line method. Background Technology

[0002] To improve energy efficiency and reduce carbon emissions, integrated energy systems combining electricity, heat, and gas have been promoted and developed globally in recent years. The close coupling of multiple energy forms not only brings good energy and environmental benefits but also introduces complex operational characteristics and risks. Digital time-domain simulation has become an essential tool for the planning, design, and operation control of integrated energy systems. The transmission processes of each energy subgrid profoundly affect the dynamic characteristics of the system; therefore, achieving efficient simulation of the transmission of different energy flows (electricity, heat, and gas) is of great significance for studying the operational characteristics of integrated energy systems. Existing power network simulation technology is relatively mature, and related research mainly focuses on the simulation of gas and heat networks.

[0003] In response, efficient simulation methods for subnet power flow transmission dynamics have been developed both domestically and internationally, which can be broadly categorized into three main types:

[0004] 1) Finite Difference Method. This type of research utilizes numerical difference schemes to discretize the partial differential equations describing energy flow characteristics. This method possesses high numerical stability and allows for larger simulation step sizes. However, in gas network simulations, it exhibits severe numerical dispersion.

[0005] 2) Time-Frequency Transformation Method. This type of research references the methodology of deducing power networks from "fields" to "paths," guiding the analysis of gas and heat networks and constructing energy flow transfer functions for these networks. This method possesses clear physical meaning, but it requires frequency domain analysis of time-domain state variables, presenting a trade-off between computational accuracy and efficiency.

[0006] 3) Method of Characteristics. This method transforms the partial differential equations into a set of ordinary differential equations by appropriately selecting the simulation time step and spatial step. In gas network simulations, there is no numerical dispersion. However, to ensure high simulation accuracy, the simulation step size is usually small.

[0007] The above methods all start by solving the energy flow transport equation, without considering the actual characteristic of the difference in pipeline length between gas networks and heat networks.

[0008] Based on the characteristics of gas and heat networks, this invention proposes a high-performance variable step size characteristic line method, which is suitable for the simulation of gas and heat network subnetworks in integrated energy systems. Summary of the Invention

[0009] This invention provides a simulation method for gas and heat subnets of integrated energy systems based on the variable step size characteristic line method. Starting from the large differences in pipe length in actual gas and heat networks, this invention adopts the variable step size concept, which can significantly reduce the number of calculations for longer pipes in the network and improve simulation efficiency while ensuring the simulation accuracy of each branch.

[0010] The technical solution of the present invention is as follows:

[0011] A simulation method for gas and heat subnets of a comprehensive energy system based on the variable step-size characteristic line method is characterized by the following steps:

[0012] Step 1: Construct the network topology files for the gas network and heating network:

[0013] The gas network topology file consists of a gas network node file and a gas network pipeline branch file. The gas network node file is constructed using "node number, node type, and injection flow rate". The node type is divided using "0 / 1 / 2", where "0" indicates that the node is a source node, "1" indicates that the node is an intermediate node, and "2" indicates that the node is a load node. The injection flow rate represents the mass flow rate of gas flowing into the node. The gas network pipeline branch file is constructed using "branch number, starting node, ending node, pipeline length, pipeline diameter, and pipeline friction coefficient".

[0014] The heating network topology file consists of heating network node files and heating network pipeline branch files. The heating network node files are constructed using "node number, node type, and injected flow rate," and the node type is divided using "0 / 1 / 2." The heating network pipeline branch files are constructed using "branch number, starting node, ending node, pipeline length, pipeline diameter, pipeline thermal resistance, and mass flow rate," where the mass flow rate is the mass flow rate of the thermal working fluid flowing through the pipeline, and is a constant value.

[0015] Step 2: Based on the pipeline branch file, use the simulation time step of the shortest pipeline branch as the base step size Δt. b Select the simulation time step Δt and simulation space step Δx for the pipeline branch;

[0016] Step 2.1 Based on the aforementioned pipe branch file, select the shortest pipe branch in the network and determine the simulation space step size Δx for the shortest pipe branch. l The formula is as follows:

[0017]

[0018] In the formula, L l M is the length of the shortest pipe branch. l The required number of space segments;

[0019] Step 2.2 Calculate the ratio of the length of the remaining pipe branch to the length of the shortest pipe branch. The formula is as follows:

[0020]

[0021] In the formula, i is the pipe branch index, L i Let be the length of the i-th pipe branch;

[0022] Step 2.3 Set the time step of the shortest pipe branch as the base step and select the simulation space step for each pipe branch. The formula is as follows:

[0023]

[0024] Step 2.4 Set the time step of the shortest pipe branch as the base time step, and select the simulation time step for each pipe branch. The formula is as follows:

[0025]

[0026] In the formula, k 11 ~k 23 It is a constant;

[0027] Step 3: Construct the network association matrix based on the network node file and the pipeline branch file:

[0028] Step 3.1 Construct a node-inflow branch correlation matrix A using network nodes and inflow branches. in :a in,ij =1 indicates that "material flow g" flows from pipe j into node i;

[0029] Step 3.2 Construct a node-outflow branch correlation matrix A using network nodes and outflow branches. out :a out,ij =1 indicates that "material flow g" flows from node i to pipe j;

[0030] Step 3.3 Construct the branch head-node association matrix A using the branch head and network nodes. pn1 :a pn1,ij =1 indicates that the potential f at the beginning of branch i is equal to the f of node i;

[0031] Step 3.4 Construct the branch end-node association matrix A using branch ends and network nodes. pn2 :a pn2,ij =1 indicates that the potential f at the end of branch i is equal to the f at node i;

[0032] Step 4: Preprocess the pipe branches with a simulation time step larger than the base time step using the method of characteristics. This involves updating the intermediate nodes of the branches with a simulation time step larger than the base time step and calculating the coefficients γ0 and γ of the branch. M The formula is as follows:

[0033]

[0034]

[0035]

[0036] In the formula, f and g are the state variables of the network nodes, where f represents "potential" and g represents "material flow"; the subscript r is the spatial point index, and 0 and M represent the beginning and end of the branch, respectively; the superscript n indicates the nth time step; α 1~4 β 1~4 c1 and c2 are constants, determined by the branch parameters; c1 and c2 are constant coefficient matrices, determined by the branch parameters.

[0037] Step 5: Set the time variable t = 0, the branch index i = 0, and the time interval variable Δt. c = 0;

[0038] Step 6: Time variable t = t + Δt b ;

[0039] Step 7: Solve for the state variables of network nodes by combining the network incidence matrix and the characteristic line method:

[0040] In the gas network, the pressure P at the source node is known. sr and load node flow q ld And the injected traffic q of the intermediate node int If the value is 0, the unknown state variables of each node are represented as follows:

[0041]

[0042]

[0043] In the formula, P ld P represents the pressure at the load node. int The pressure at the intermediate node; q sr For source node traffic; Y 11~33 B is a constant matrix, determined by the gas network structure; int B ld It is a constant matrix, determined by the distribution of the network's state variables in the previous time step;

[0044] In a heating network, given the source node temperature, then the intermediate node temperature T is... int and load node temperature T ldRepresented as:

[0045]

[0046] In the formula, T sr Y represents the source node temperature. 21~32 B is a constant matrix, determined by the structure of the heating network; int B ld It is a constant matrix, determined by the distribution of the network's state variables in the previous time step;

[0047] Step 8: Branch index i = i + 1, the time interval variable Δt for branch i. c,i = Δt c,i + Δt b ;

[0048] Step 9: Determine Δt c,i Is it greater than or equal to the simulation time step Δt of this branch? i If Δt c,i ≤ Δt i If Δt is true, then branch i is not processed; otherwise, Δt is true. c,i = 0, branch i completes one branch state variable calculation, the calculation process is the same as step 4;

[0049] Step 10: Determine if the branch index i is greater than the total number of branches: if yes, set i = 0 and proceed to step 11; otherwise, return to step 8.

[0050] Step 11: Determine if time t is greater than the set total simulation time T. max If the result is greater than 0, the simulation ends; otherwise, return to step 6.

[0051] Compared with existing technologies, the beneficial effects of this invention are that it can improve the dynamic simulation efficiency of integrated energy systems, help to understand the dynamic operating characteristics of gas networks and heat networks, and facilitate the design, analysis and operation control of integrated energy systems. Attached Figure Description

[0052] Figure 1 This is a flowchart of the simulation method for the gas and heat subnet of a comprehensive energy system based on the variable step size characteristic line method of the present invention.

[0053] Figure 2 This is a schematic diagram of the numerical format of the method of characteristics.

[0054] Figure 3 This is a diagram describing network nodes and pipeline branches.

[0055] Figure 4 This is a schematic diagram of a 34-node natural gas network.

[0056] Figure 5This is a schematic diagram of a 51-node thermal network.

[0057] Figure 6 This is a schematic diagram of the simulation results of the gas network node flow.

[0058] Figure 7 This is a schematic diagram of the simulation results of the pressure at the nodes of the gas network.

[0059] Figure 8 This is a schematic diagram of the simulation results of the temperature at the nodes of the heating network. Detailed Implementation

[0060] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0061] Please see Figure 1 , Figure 1 The flowchart of the integrated energy system gas and heat subnet simulation method based on the variable step size characteristic line method of this invention is shown in the figure. It includes 11 steps. Step 1: Based on the network structure of the gas network, heat network and other subnets in the integrated energy system, construct the network node file and pipeline branch file of each subnet. Step 2: According to the pipeline branch file, select the simulation time Δt and spatial step size Δx of each branch, and determine the minimum time step size as the base step size Δt. b Step 3: Construct a network association matrix based on the network node file and pipeline branch file to describe the network topology. Step 4: Use the characteristic line method to preprocess branches with a time step larger than the base step. Step 5: Set the time variable t = 0, the branch index i = 0, and the time interval variable Δt. c = 0; Step 6, Time variable t = t + Δt b Step 7: In the simulation iteration, the basic step size Δt b Next, combining the network correlation matrix and the characteristic line method, the state variables of each node are solved; Step 8, branch index i = i + 1, the time interval variable Δt of branch i. c,i = Δt c,i + Δt b Step 9: Determine the time interval variable Δt for branch i. c,i If the simulation time step is greater than or equal to that of the branch, then use the method of characteristics and the obtained nodal state variables to calculate the state variables of branch i and set Δt. c,iSet to 0, otherwise do not calculate; Step 10: If the branch index i is greater than the total number of branches, set i to 0, otherwise return to step 8; Step 11: If the time variable t is greater than the set total simulation time, end the simulation, otherwise return to step 6. Among them, steps 1 to 5 constitute the initialization phase of the simulation method, which only needs to be executed once at the beginning of the program, and steps 6 to 11 constitute the main loop of the simulation method.

[0062] First, the principles of the gas network and heat network transmission equations and the method of characteristics used in this invention will be explained.

[0063] In the gas network, gas pressure P and gas flow rate q are selected as state variables, and in the heat network, temperature T and heat power Q are selected as state variables. Then, the transmission characteristics of gas and heat energy flow can be expressed as Equation (1) and Equation (2), respectively:

[0064] (1)

[0065] Where t and x are time and space variables, respectively, S is the cross-sectional area of ​​the pipe, D is the pipe diameter, and v a For the speed of sound in gas, v b Let λ be the average flow velocity in the pipe. g The coefficient of friction of the pipeline.

[0066] (2)

[0067] Where m is the mass flow rate, v is the flow velocity, and λ is the mass flow rate. h For the thermal resistance of the pipe, c w This is the specific heat capacity of the working fluid.

[0068] According to equations (1) and (2), the characteristics of gas and heat energy flow can be expressed in a unified form:

[0069] (3)

[0070] Where f and g are state variables, f represents "potential", g represents "matter flow", and k 11 ~k 23 It is a constant.

[0071] Using the method of characteristics, by appropriately selecting the time step and spatial step, the partial differential equation (3) can be transformed into a system of ordinary differential equations (4), such as... Figure 2 As shown:

[0072] (4)

[0073] Therefore, in the gas network and the heat network, the time step Δt and the spatial step Δx respectively satisfy:

[0074] (5)

[0075] In this context, the subscripts N and H represent the gas network and the heating network, respectively.

[0076] The steps of this invention will be described in further detail below.

[0077] Step 1: Forming subnet network node files and pipeline branch files

[0078] like Figure 3 As shown, network node files and pipeline branch files are used to describe the entire network topology. Gas networks and heating networks have different energy flow transmission characteristics, therefore the parameters of the constructed files differ.

[0079] In the gas network, network nodes are described using "node number, node type, and injected flow rate." "Node type" is categorized using "0 / 1 / 2." A node type of "0" indicates a source node; "1" indicates an intermediate node; and "2" indicates a load node. "Injected flow rate" represents the mass flow rate of gas flowing into that node. Pipeline branch files are described using "branch number, starting node, ending node, pipe length, pipe diameter, and pipe friction coefficient."

[0080] In a heating network, the description of network nodes is constructed using "node number, node type, and injected flow rate." Similar to gas networks, the "node type" is categorized using "0 / 1 / 2." This invention considers the heating network to operate under mass regulation, i.e., constant mass flow rate and variable temperature. Therefore, the "injected flow rate" represents the mass flow rate of the thermal medium flowing into that node, which is a constant value. The description of pipeline branch files is constructed using "branch number, starting node, ending node, pipeline length, pipeline diameter, pipeline thermal resistance, and mass flow rate." Here, "mass flow rate" refers to the mass flow rate of the thermal medium flowing through the pipeline, which is also a constant value.

[0081] Step 2: Selecting the Simulation Step Size for Branches

[0082] First, based on the pipeline branch file, select the shortest pipeline in the network and determine the simulation space step size Δx for the shortest pipeline. l :

[0083] (6)

[0084] Among them, L l M is the length of the shortest pipe. l The required number of space segments.

[0085] Then, the length ratio of the remaining branch to the shortest branch is calculated, and the spatiotemporal step size L of each branch is selected. i And the shortest branch L l The time step size is set as the base step size Δt b :

[0086] (7)

[0087] Where k is the length ratio and i is the branch index.

[0088] Step 3: Generation of the association matrix

[0089] Based on the start and end point information of each branch in the pipeline branch file, an association matrix A can be generated. in A out A pn1 A pn2 :

[0090] 1) Node-inflow branch correlation matrix A in :a in,ij =1 indicates that "material flow g" flows from pipe j into node i;

[0091] 2) Node-Outflow Branch Association Matrix A out :a out,ij =1 indicates that "material flow g" flows from node i to pipe j;

[0092] 3) Branch head-node association matrix A pn1 :a pn1,ij =1 indicates that the potential f at the beginning of branch i is equal to the f of node i;

[0093] 4) Branch end-node correlation matrix A pn2 :a pn2,ij =1 indicates that the potential f at the end of branch i is equal to the f at node i.

[0094] Step 4: Branch Preprocessing

[0095] Depend on Figure 2 It can be seen that each intermediate node of the pipeline has two characteristic lines C+ and C- passing through it. The state variables of the intermediate nodes can be obtained using equation (4). Equation (4) can be discretized using the trapezoidal rule. If the state variables of the intermediate nodes are known... and Then the state quantity at the next moment It can be represented as:

[0096] (8)

[0097] Where c1 and c2 are constant coefficient matrices.

[0098] like Figure 2 As shown, the xt plane can be divided into layers E and O based on the presence or absence of pipe boundary points. The intermediate nodes between layers E and O and back to layer E (i.e.,...) Figure 2The state variables of the black nodes in the equation are calculated and recorded as an intermediate node update (which can be calculated using equation (8)). Since only one feature line passes through the boundary point, the state variables of the boundary point cannot be directly solved. Therefore, the state variable g of the boundary point is first expressed as a function of f, and the corresponding coefficients γ0 and γ are calculated. M :

[0099] (9)

[0100] (10)

[0101] The preprocessing of branches involves updating intermediate nodes for branches whose simulation time step is greater than the base time step, and calculating the coefficients γ0 and γ of that branch. M .

[0102] Step 5: Variable Setting

[0103] Set the time variable t = 0, the branch index i = 0, and the time interval variable Δt. c = 0.

[0104] Step 6: Update the time variable

[0105] The time variable t = t + Δt b .

[0106] Step 7: Solving for network node state variables

[0107] Based on the correlation matrix, the injected traffic g of each node n It can be represented as:

[0108] (11)

[0109] Among them, g0, g M The set of flow rates at the beginning and end points of the pipeline, with a size of N. b ×1, N b This represents the number of branch roads.

[0110] The set of potentials of each node is represented as f. n Then the sets of potential values ​​at the beginning and end points of each branch, f0 and f1, are... M It can be represented as:

[0111] (12)

[0112] Combining equations (9) to (12), g n and f n satisfy:

[0113] (13)

[0114] Where Y is a constant coefficient matrix, which is only related to the network structure and branch parameters, and B is determined by the state variables of each branch at the previous time step.

[0115] Let the subscripts sr / int / ld represent the source node, intermediate node, and load node, respectively. Then, equation (13) can be rewritten as:

[0116] (14)

[0117] In a gas network, the pressure P at the source node is usually known. sr and load node flow q ld And the injected traffic q of the intermediate node int The value is 0. Therefore, the unknown state variables of each node can be expressed as:

[0118] (15)

[0119] (16)

[0120] In a quality-regulated heating network, the heat power Q and temperature T have a linear relationship, so equation (13) can be rewritten as:

[0121] (17)

[0122] Where, d sr and d ld These are the injected traffic from the source node and the load node, respectively.

[0123] Typically, if the temperature of the source node is known, the temperatures of the intermediate nodes and the load nodes can be expressed as:

[0124] (18)

[0125] In the simulation iteration, at each basic step size, the state variables of each node can be calculated using equations (15) and (18).

[0126] Step 8: Update branch index and branch time interval variables

[0127] Branch index i = i + 1, and the time interval variable Δt for branch i. c,i = Δt c,i + Δt b .

[0128] Step 9: Branch update judgment

[0129] Determine the time interval variable Δt of the branch c,i Is it greater than or equal to the simulation time step Δt of this branch? i If Δt c,i < Δt iIf Δt is true, then branch i is not processed; otherwise, Δt is true. c,i = 0, branch i completes one branch state variable calculation, the calculation process is the same as step 4.

[0130] Step 10: Determine whether all branches have been traversed.

[0131] Determine if the branch index i is greater than the total number of branches: if yes, set i = 0 and proceed to step 11; otherwise, return to step 8.

[0132] Step 11: Simulation End Judgment

[0133] If the time variable t is greater than the set total simulation time, the simulation ends; otherwise, return to step 6.

[0134] Application Examples

[0135] The following will use a 34-node natural gas network (such as...) Figure 4 (as shown) and 51-node thermal networks (such as...) Figure 5 Taking the example shown, the implementation effects of the present invention will be further illustrated, but this should not be used to limit the scope of protection of the present invention. In the results of the embodiments, M1 represents the present invention, and M2 represents the traditional feature line method.

[0136] The simulation results of the air network are as follows Figure 6 , Figure 7 As shown. Within 1-2 hours, load node N n34 The flow rate at the source node increases from 0.1 kg / s to 0.5 kg / s, while the flow rates at other load nodes remain unchanged; within 2.5 h to 3.5 h, the flow rate at source node N... n1 The air pressure rose from 200 kPa to 204.5 kPa. This shows that when the load node N... n34 When the flow rate changes, this node and the load node N n18 The pressure will change. And when the source node N... n1 When the pressure at node N changes, the pressure at the load node also changes accordingly, with a certain time delay. n18 and node N n34 With source node N n1 The time delay varies depending on the distance. The greater the distance, the longer the time delay. The simulation results for the heat network are shown in Figure 8. Due to the slower flow rate of the working fluid, the time delay in the heat network is more significant compared to that in the gas network. Furthermore, the further away from the source node, the longer the time delay.

[0137] Table 1

[0138]

[0139] Table 1 compares the implementation effects of the present invention with those of the traditional characteristic line method (including mean relative error MRE and simulation running time t). rAs can be seen, this invention has high simulation accuracy, with average relative errors all less than 0.5%. Furthermore, simulation time is significantly reduced, with gas network simulation time reduced by approximately 38.54% and heating network simulation time reduced by approximately 51.09%.

[0140] Finally, it should be noted that the described embodiments are only some, not all, of the embodiments in this application. All other embodiments obtained by those skilled in the art based on the embodiments in this application without inventive effort are within the scope of protection of this application.

Claims

1. A simulation method for gas and heat subnets of a comprehensive energy system based on the variable step-size characteristic line method, characterized in that, Includes the following steps: Step 1: Construct the network topology files for the gas network and heating network: The gas network topology file consists of a gas network node file and a gas network pipeline branch file. The gas network node file is constructed using "node number, node type, and injection flow rate". The node type is divided using "0 / 1 / 2", where "0" indicates that the node is a source node, "1" indicates that the node is an intermediate node, and "2" indicates that the node is a load node. The injection flow rate represents the mass flow rate of gas flowing into the node. The gas network pipeline branch file is constructed using "branch number, starting node, ending node, pipeline length, pipeline diameter, and pipeline friction coefficient". The heating network topology file consists of heating network node files and heating network pipeline branch files. The heating network node files are constructed using "node number, node type, and injected flow rate," and the node type is divided using "0 / 1 / 2." The heating network pipeline branch files are constructed using "branch number, starting node, ending node, pipeline length, pipeline diameter, pipeline thermal resistance, and mass flow rate," where the mass flow rate is the mass flow rate of the thermal working fluid flowing through the pipeline, and is a constant value. Step 2: Based on the pipeline branch file, use the simulation time step of the shortest pipeline branch as the base step size Δt. b Select the simulation time step Δt and simulation space step Δx for the pipeline branch; Step 2.1 Based on the aforementioned pipe branch file, select the shortest pipe branch in the network and determine the simulation space step size Δx for the shortest pipe branch. l The formula is as follows: In the formula, L l M is the length of the shortest pipe branch. l The required number of space segments; Step 2.2 Calculate the ratio of the length of the remaining pipe branch to the length of the shortest pipe branch. The formula is as follows: In the formula, i is the pipe branch index, L i Let be the length of the i-th pipe branch; Step 2.3 Set the time step of the shortest pipe branch as the base step and select the simulation space step for each pipe branch. The formula is as follows: Step 2.4 Set the time step of the shortest pipe branch as the base time step, and select the simulation time step for each pipe branch. The formula is as follows: In the formula, k 11 ~k 21 It is a constant; Step 3: Construct the network association matrix based on the network node file and the pipeline branch file: Step 3.1 Construct a node-inflow branch correlation matrix A using network nodes and inflow branches. in :a in,ij =1 indicates that "material flow g" flows from pipe branch i into node j; Step 3.2 Construct a node-outflow branch correlation matrix A using network nodes and outflow branches. out :a out,ij =1 indicates that "material flow g" flows from node j to pipe branch i; Step 3.3 Construct the branch head-node association matrix A using the branch head and network nodes. pn1 :a pn1,ij =1 indicates that the potential f at the beginning of branch i is equal to the f at node j; Step 3.4 Construct the branch end-node association matrix A using branch ends and network nodes. pn2 :a pn2,ij =1 indicates that the potential f at the end of branch i is equal to the f at node j; Step 4: Preprocess the pipe branches with a simulation time step larger than the base time step using the method of characteristics. This involves updating the intermediate nodes of the branches with a simulation time step larger than the base time step and calculating the coefficients γ0 and γ of the branch. M The formula is as follows: In the formula, f and g are the state variables of the network nodes, where f represents "potential" and g represents "material flow"; the subscript r is the spatial point index, and 0 and M represent the beginning and end of the branch, respectively; the superscript n indicates the nth time step; α 1~4 β 1~4 c1 and c2 are constants, determined by the branch parameters; c1 and c2 are constant coefficient matrices, determined by the branch parameters. Step 5: Set the time variable t = 0, the branch index i = 0, and the time interval variable Δt. c = 0; Step 6: Time variable t = t + Δt b ; Step 7: Solve for the state variables of network nodes by combining the network incidence matrix and the characteristic line method: In the gas network, the pressure P at the source node is known. sr and load node flow q ld And the injected traffic q of the intermediate node int If the value is 0, the unknown state variables of each node are represented as follows: In the formula, P ld P represents the pressure at the load node. int The pressure at the intermediate node; q sr For source node traffic; Y 11~33 B is a constant matrix, determined by the gas network structure; int B ld It is a constant matrix, determined by the distribution of the network's state variables in the previous time step; In a heating network, given the source node temperature, then the intermediate node temperature T is... int and load node temperature T ld Represented as: In the formula, T sr Y represents the source node temperature. 21~32 B is a constant matrix, determined by the structure of the heating network; int B ld It is a constant matrix, determined by the distribution of the network's state variables in the previous time step; Step 8: Branch index i = i + 1, the time interval variable Δt for branch i. c,i = Δt c,i + Δt b ; Step 9: Determine Δt c,i Is it greater than or equal to the simulation time step Δt of this branch? i If Δt c,i ≤ Δt i If Δt is true, then branch i is not processed; otherwise, Δt is true. c,i = 0, branch i completes one branch state variable calculation, the calculation process is the same as step 4; Step 10: Determine if the branch index i is greater than the total number of branches: if yes, set i = 0 and proceed to step 11; otherwise, return to step 8. Step 11: Determine if time t is greater than the set total simulation time T. max If the result is greater than 0, the simulation ends; otherwise, return to step 6.

Images