Dynamic simulation method for annular pipe network and mesh pipe network based on DES (Data Encryption Standard)

Abstract

The invention discloses a dynamic simulation method for an annular pipe network and a net-shaped pipe network based on DES, and belongs to the technical field of regional heating system DH online planning and optimization, complex heating network model development and simulation and the like. An accurate and efficient variable time step model suitable for annular and netlike pipe networks is developed by using DES (Data Encryption Standard). A hydraulic model which converts a complex net-shaped structure into a tree-shaped structure easy to process and can simulate water flow in a complex annular and net-shaped heat supply network and a dynamic thermodynamic model which is expanded from one-way flow to support two-way flow are established, and the hydraulic model and the dynamic thermodynamic model are integrated in a DES frame. The model can dynamically define time and space discretization, effectively improves the accuracy and efficiency of simulation, is suitable for many application scenes such as the regional heat supply industry and energy communities, has certain potential to be integrated into an online operation optimization model of a regional heat supply system, and therefore more intelligent and low-carbon heat supply is achieved. And complex iterative calculation required by network operation is effectively managed.

Description

Technical Field

[0001] The present invention belongs to the technical fields of online planning and optimization of district heating systems (DH), development and simulation of complex heating network models, and particularly relates to the development of a multifunctional dynamic hydraulic and thermal model using a discrete event simulation (DES) method, which can simulate a district heating network with a complex mesh topology and diversified operation strategies, and realize accurate and efficient variable-step spatiotemporal simulation. Background Art

[0002] Modern district heating (DH) systems present new challenges for network modeling due to their complex topology and operation strategies. In tree and ring topology networks, the centrally concentrated heat production makes the flow direction in each pipe unidirectional. However, in mesh networks with loops and networks with multiple production sites, the flow direction in the pipes may change depending on the hydraulic conditions. Models for such networks need to support not only bidirectional flows and decentralized production, but also scalability in terms of computational speed and memory consumption. Dénarié et al. developed a thermohydraulic model in MATLAB with a fixed time step. The Lagrangian model was applied to a tree network with 485 edges and completed a 24-hour simulation in about 4 seconds, which is more than 100 times faster than the corresponding finite volume method (FVM) model. However, this Lagrangian method has not yet been extended to mesh networks, which limits its application in more complex network structures. Boghetti et al. developed an open source mesh network model and used Python for more efficient hydraulic simulations. However, the authors pointed out that computational speed is still a limiting factor and computational efficiency remains a bottleneck. In addition, choosing the right time resolution to balance computational overhead and accuracy is a key issue in energy system simulation. In a district heating system, different components, such as circulation pumps and boilers, require different time resolutions to accurately capture temperature fluctuations. Therefore, it is particularly important to use variable-step simulation to simulate the network connecting different components to ensure that the system simulation is both accurate and efficient.

[0003] In view of the above problems, a dynamic simulation method for loop pipe network and mesh pipe network based on DES is proposed. This method is suitable for grid heating network with multiple heat sources and cycles. It can efficiently simulate the operation of heating system under complex control strategies, covering the dynamic changes of water supply temperature and flow from heating plant to thermal power station, and improving simulation accuracy and efficiency. Summary of the invention

[0004] The purpose of the present invention is to provide a dynamic simulation method for a ring pipe network and a mesh pipe network based on DES.

[0005] The technical solution of the present invention:

[0006] A DES-based ring pipe network and mesh pipe network dynamic simulation method comprises the following steps:

[0007] Step S1, developing a hydraulic model capable of simulating water flow in complex annular and meshed heating networks;

[0008] S1.1 Mathematical description. Expanding from a tree network to a mesh network with one or more loops, in which hydraulic simulation and thermal simulation are decoupled, this decoupling allows the mass flow of each part of the network to be determined independently before thermal analysis;

[0009] The relationship between the number of cycles M, the number of vertices V, and the number of edges E in a mesh network is:

[0010] E=M+V-1 (1) Solving the mass flow of a mesh network requires at least E independent equations.

[0011] For a branched pipe network without loops, M = 0, and at least V-1 equations are required to solve the mass flow of the branched pipe network. Assuming no water leakage, each node maintains mass balance, expressed as:

[0012]

[0013] Among them, D j It is the water demand or supply of user and production site nodes, and the internal nodes are 0. j is a set of pipes connected to node j. i is the mass flow in pipe i, which is positive when flowing in the normal direction and negative when flowing in the reverse direction. The sign is positive for normal inflow and negative for normal outflow. Since the total mass input in the network must match the total demand ∑D j =0, thus forming V-1 linearly independent node balance equations, thus forming a definite equation system for calculating the mass flow rate of each pipeline in the tree network.

[0014] Meshes containing loops require additional equations to solve for mass flow, since each loop introduces an additional edge (variable). To solve this problem, M additional independent equations are required, which are beyond the scope of the node mass balance equation. Based on Kirchhoff's second law, which states that in any closed loop, the sum of the pressure differences along the loop must be zero, the pressure equation for each loop c is derived:

[0015]

[0016] Among them, △P i (m i ) indicates the flow rate m iThe pressure in the lower pipeline i decreases. It is positive when flowing in the normal direction and negative when flowing in the reverse direction. The sign of each pressure term depends on whether the normal direction of the pipeline is consistent with the loop direction. The pressure loss in the pipeline is due to frictional loss and local pressure loss. Since no booster pump is required in the network, the frictional loss is quantified by the Darcy - Weisbach equation, where the variables include the pipeline length L, pipeline diameter D, water density ρ, flow velocity v, and friction coefficient f:

[0017]

[0018] where m represents the flow rate and A represents the cross - sectional area of the pipeline;

[0019] The friction coefficient f depends on the flow conditions:

[0020]

[0021] ε is the absolute roughness of the pipe wall:

[0022] Re = |v|Dρ / μ (6) where μ is the dynamic viscosity of water. To simplify the calculation, a validated explicit approximation formula is chosen to replace the implicit Colebrook - White equation for the friction factor in turbulent flow (Re≥4000). In the range of transitional flow (2000 < Re < 4000), linear interpolation is used.

[0023] S1.2 Solve the mass flow rate in the tree - shaped network. For a tree - shaped network with only one production node, assume that the normal direction of the pipeline is consistent with the only flow direction, that is, in the water supply network, it flows from the production node to the user node, and in the return water network, it flows from the user node back to the production node. The mass flow rate of all pipelines is calculated in a bottom - up manner, starting from the user nodes, and aggregating the flow rates at the confluence points until reaching the production node. Topological sorting can be used to pre - determine the order of calculating the pipeline flow rates. Since the tree - shaped topology remains static throughout the simulation, this order will remain unchanged after initialization.

[0024] S1.3 Solve the mass flow rate in the mesh network. After establishing the initial water distribution, the mass flow rate of the boundary nodes is determined based on the flow rates of the connected non - loop edges, where a positive value indicates flow into the area and a negative value indicates flow out of the area. Specifically, the mass flow rate of the boundary nodes is obtained by calculating the flow rates of the non - loop edges connected to these nodes. For each boundary node, if the flow rate of the connected non - loop edge is positive, it means water flows into the block where the node is located; if the flow rate is negative, it means water flows out of the node. In this way, the flow rates of the boundary nodes can accurately reflect the input and output of the system, ensuring mass conservation and balance in the entire network.

[0025] To handle the nonlinear subproblems in each loop block, the Newton-Raphson algorithm is employed with the following iterative steps:

[0026] (1) Calculate the pressure drop △P of each pipeline i (m i ) and its derivative d△P i (m i ) / dm i .

[0027] (2) Following the above method, the pressure drop in each cycle is accumulated, and its derivative is accumulated in the same way.

[0028] (3) Evaluate the pressure difference imbalance in each cycle. If the pressure difference error P of all cycles in the block j If it is zero within the tolerance range, the system is considered to have converged. If not, the updated mass flow step length △m is calculated using the following method: (k) .

[0029] Δm (k) =m (k+1) -m (k) =-J -1 P (7)

[0030] Among them, m (k) represents the mass flow step size of the kth iteration, J represents the Jacobian matrix, and P represents the pressure error;

[0031] (4) Update the mass flow of all edges in the block based on the updated flow of the truncated edges and the initial flow of the boundary nodes.

[0032] △m (k) represents the update step size in the kth iteration. A positive step size indicates an increase in mass flow along the default circulation direction. The Jacobian matrix (J) is:

[0033]

[0034] Among them, P n Indicates pressure error; m n represents mass flow rate;

[0035] Step S2, extending the dynamic thermal model from unidirectional flow to support bidirectional flow;

[0036] S2.1 Calculate the bidirectional flow in the pipe. The temperature of water particles in the pipe decreases exponentially with the travel time τ due to heat conduction to the surroundings:

[0037]

[0038] Among them, T inis the temperature of the water when it enters the pipe, T out is the temperature of the water as it exits the pipe. k1 and k2 are coefficients for single, parallel, and double pipes. The Lagrangian DES method involves tracking the motion of the water front, which is an infinitely thin segment of water moving inside the pipe. The water front acts as a dynamic sampling point, recording various inlet parameters such as creation time, temperature, velocity, and distance to adjacent fronts. The water front is defined at the breakpoints in the temperature profile. Numerical analysis shows that by identifying these breakpoints and assuming a linear transition of temperature between adjacent points, the water temperature can be predicted with high accuracy. The bidirectional flow model follows the same principle, but the flow direction is taken into account when calculating the travel time.

[0039] The travel time of each particle is determined by the arrival time t out and entry time t in =Calculated by the difference between |v1| and |v2|. Assume that: at t1, the temperature of the current inlet (left end) changes, the first and only water front F1 is generated, and its arrival event is scheduled according to the flow rate |v1|; at t2, when the flow stops, water fronts (F2 and F3) are generated at the inlet and outlet, and the arrival event is cancelled; at t3, when the flow direction is reversed, water front F2 reaches the current target node (left end), and the arrival event of F1 is scheduled according to the flow rate |v3|. At the same time, water front F4 is generated at the current inlet (right end). Fronts F3 and F4 are extremely close, indicating a step change in temperature; at t4, front F1 reaches the current outlet, and F3 becomes the next front to arrive, and its arrival time is scheduled according to |v3|. The travel distance of particle F1 in the time period from t1 to t2 is equal to the travel distance in the time period from t3 (reverse flow—flow rate changes from 0 to v3) to t4 (particle F1 reaches the end of the pipe). The entry time and travel time of F1 can be derived from this principle and further described as follows:

[0040]

[0041] in, represents the entry time of particle F2; represents the entry time of particle F1; represents the arrival time of particle F1; represents the arrival time of particle F2; Δt represents the time interval; v1 represents the flow velocity at the time of entry; v3 represents the tassel at the time of arrival; represents the travel time of particle F1; represents the travel time of particle F2;

[0042] Likewise, the travel time and entry time of any arbitrary water front can be calculated based on the previous arrival time, the corresponding entry time, the current velocity, and the corresponding entry velocity.

[0043] Taking into account the sign of the velocity, the equations for calculating the travel time and entry time of a particle in a bidirectional flow agree with the equations for a unidirectional flow and are given by:

[0044]

[0045]

[0046] in, represents the entry time after the kth iteration; represents the arrival time after the kth iteration; v (k) represents the flow rate after the kth iteration; represents the flow rate at the time of entry after the kth iteration; τ (k) represents the travel time after the kth iteration; t (k) represents the time after the kth iteration;

[0047] There are four states of flow: positive, +0, -0, and negative. Positive and negative flow refer to the flow in the normal and reverse directions of the pipe, respectively. +0 and -0 refer to the states where the flow stops after positive or negative flow. The travel time function of a water particle through a pipe is piecewise linear and resets in two cases: when the flow direction is reversed or when the flow is temporarily stopped. This includes the following cases: When the flow direction changes, the travel time is reset to zero. During the period when the flow stops, the travel time is undefined. When the flow resumes, the reset travel time value will depend on the relative direction of the flow before and after the stop. If the flow direction is reversed, the travel time is set to the duration of the stop. If the flow resumes in the original direction, the travel time is extended by the duration of the stopped flow.

[0048] S2.2 Analyze the effect of flow direction changes on adjacent pipes. Changes in flow direction within a pipe can significantly affect the inlet temperature of adjacent pipes, so it is necessary to introduce additional water fronts to monitor these discontinuities. The mixed temperature at node j is is the outlet temperature of each pipe It is calculated by taking the weighted average of:

[0049]

[0050] in Represents the set of pipes that flow into node j.

[0051] Step S3, integrating hydraulic and dynamic thermal calculations in a discrete event simulation DES framework;

[0052] The DES model represents a parallel process as a series of events in time order that represent changes in the state of the system. During the simulation process, events can be scheduled, rescheduled, or canceled. These events are managed uniformly in an event queue, which is organized according to the scheduled activation time of the events. The discrete event simulation process involves iteratively removing the oldest event from the queue and executing the event routines associated with it until the event queue is empty. The dynamic scheduling capability of discrete event simulation makes it particularly suitable for simulations that require variable time steps.

[0053] The specific process is as follows:

[0054] Initialize the event queue with the mass flow change event and the inlet change event at the start time. If the event queue is empty, the process ends. If it is not empty, remove the first event from the event queue and process it.

[0055] If the first event is a mass flow change event, the mass flow is solved through the hydraulic model, and the node temperature is updated based on the previous mass flow. Each pipe is traversed to update the flow direction, the frontier double-ended queue is maintained, a new water flow front is inserted, and the arrival time of the event is updated. The node temperature is updated according to the updated mass flow. If the node temperature changes, an additional water flow front is inserted. If there is a next change, the mass flow change event is rescheduled. The event queue is maintained by re-minimum heap, and finally the event queue is determined again to see if it is empty, and the cycle repeats.

[0056] If the first event is a water front arrival event, the outlet temperature is updated and the arriving water front is removed from the double-ended queue. If the double-ended queue is empty, the inlet temperature change event of all downstream objects is triggered. If the double-ended queue is not empty, the arrival time needs to be rescheduled before triggering the inlet temperature change event of all downstream objects. Finally, it is determined again whether the event queue is empty, and the cycle repeats.

[0057] If the first event is an inlet temperature change event, a new water front is created in the pipeline. When the double-ended queue is not empty, the event queue is directly checked to see if it is empty. If the double-ended queue is empty, an arrival event needs to be dispatched, and then the event queue is checked to see if it is empty, and the cycle repeats.

[0058] Beneficial effects of the present invention: A DES-based dynamic simulation method for annular pipe networks and mesh pipe networks, which converts complex mesh structures into easy-to-handle tree structures, allows nonlinear sub-problems to be solved independently, and effectively solves the distribution problem of the hydraulic system. The dynamic thermal model is extended from unidirectional flow to support bidirectional flow, and a double-ended queue is used to allow elements to be inserted and deleted from both ends, which effectively adapts to the dynamic characteristics of water front management in bidirectional flow simulation. This accurate and efficient variable time step model for mesh networks uses the DES method, which can dynamically define time and space discretization, which is crucial to improving the accuracy and efficiency of temperature simulation. BRIEF DESCRIPTION OF THE DRAWINGS

[0059] Figure 1 It is a mesh network decomposed into cyclic blocks and tree structures.

[0060] Figure 2 It is the flow chart of the meshed DH network hydraulic model.

[0061] Figure 3 It is a schematic diagram of bidirectional flow.

[0062] Figure 4 It is the water front produced in adjacent pipes when the direction of water flow changes.

[0063] Figure 5 It is a flow chart of mesh DH network modeling based on DES.

[0064] Figure 6 It is a mesh DH network layout.

[0065] Figure 7 It is the comparison between the measured and simulated values ​​of the node temperature.

[0066] Figure 8 is the relationship between simulation time and water front number.

[0067] The purpose, features and advantages of the present invention will be further described with reference to the accompanying drawings and examples. DETAILED DESCRIPTION

[0068] The following will be combined with the drawings in the examples of the present invention to clearly and completely describe the technical solutions in the examples of the present invention. Obviously, the examples described are only part of the examples of the present invention, not all of the examples. Based on the examples in the present invention, all other examples obtained by ordinary technicians in this field without making creative work are within the scope of protection of the present invention. It should be noted that the examples of the present invention use simulations to test the gain effect of the present invention, and experiments can also be used to test the gain effect of the present invention.

[0069] The present invention proposes a DES-based dynamic simulation method for annular pipe networks and mesh pipe networks, establishes a hydraulic model that can simulate water flow in complex annular and mesh heating networks and a dynamic thermal model that is expanded from unidirectional flow to support bidirectional flow, and integrates them in the DES framework, effectively improving the accuracy and efficiency of the simulation.

[0070] The present invention provides a DES-based dynamic simulation method for annular pipe networks and mesh pipe networks, which expands the hydraulic-thermal network model from a tree network to a mesh network with one or more loops, and the water distribution in the tree component is not affected by the flow in the loop, such as Figure 1 As shown. Three types of events are defined: mass flow change events, water flow front arrival events, and inlet temperature change events. Mass flow change events are critical for water distribution management, and factors that trigger these events include changes in water supply from multiple heat source plants, changes in consumer water demand, and adjustments to pressure and flow control devices (such as valves or pumps). Each pipe has a water flow front arrival event and an inlet temperature change event. The water flow front arrival events of each pipe are managed by its own double-ended queue, ensuring local processing. Inlet temperature change events are triggered by arrival events from upstream pipes, or by direct temperature adjustments at upstream nodes, such as heat source plants (in supply networks) or users (in return networks). Figure 2 The hydraulic simulation process of a mesh network is outlined. The mesh network is effectively simplified to a tree structure, the spanning tree, by selecting a cutoff edge for each loop. After initializing the mass flow rates on these cutoff edges, the remaining distribution is resolved using the established procedure for tree networks.

[0071] Figure 3 A bidirectional flow model is shown. The key events of the extended DES method applied to the pipeline include: at t1, the temperature of the current inlet (left end) changes, generating the first and only water front F1, and scheduling its arrival event according to the flow rate |v1|; at t2, when the flow stops, water fronts (F2 and F3) are generated at the inlet and outlet, and the arrival event is cancelled; at t3, when the flow direction is reversed, water front F2 reaches the current target node (left end), and the arrival event of F1 is scheduled according to the flow rate |v3|. At the same time, water front F4 is generated at the current inlet (right end). Fronts F3 and F4 are extremely close, indicating a step change in temperature; at t4, front F1 reaches the current outlet, and F3 becomes the next arriving front, and its arrival time is scheduled according to |v3|.

[0072] Figure 4Demonstrates the creation of additional water fronts in adjacent pipes when the flow direction in a pipe changes. The flow direction in pipe 1 is reversed. This reversal requires the creation of additional water fronts not only in the pipe with the reverse flow, but also in the previous downstream pipe (Pipe 3) and the new downstream pipe (Pipe 0) due to the change in mixing temperature at the node.

[0073] Figure 5 The DES flow chart for simulating a meshed heating network is shown. Unlike the simulation of a tree network, the calculation method for mass flow change events has been redesigned. In the revised model, the water temperature of the entire network is initially calculated based on the previous hydraulic conditions. Subsequent updates to individual pipes require recalculation of the water temperature to reflect the updated hydraulic conditions. This adjustment introduces a new heat front, resulting from the change in flow direction. The main advantage of this approach is its flexibility: heat updates can be performed in any order, avoiding the constraints imposed by a fixed heat propagation order when the flow direction changes. Due to the large number of arrival events that need to be rescheduled, an event queue is maintained at the end of each mass flow change event to improve simulation efficiency.

[0074] Reference Figure 6 In one embodiment of the present invention, the dynamic simulation method of the ring pipe network and the mesh pipe network based on DES simulates the supply side of the main district heating network in City A for 85 days, from December 1, 2021 to February 24, 2022. The pipe diameters of the heating network range from DN 80 to DN 900, most of which are prefabricated insulated stainless steel pipes, wrapped with a high-density polyethylene outer layer and a rigid polyurethane foam layer inside. The supply and return pipes are installed in parallel underground, and their depth and spacing meet the construction standard of "CJJ / T81-2013".

[0075] The energy supply of the network in this embodiment comes from two heat sources: a main heat source with a capacity of 200MW, which operates continuously throughout the heating season; the other is a peak boiler with a maximum output of 350MW, which is enabled during peak load periods. The main heat source is located at node 0 on the east side of the network, while the peak boiler is located at node 184 on the west side. The outdoor temperature varies between -12.5℃ and 12.3℃. The peak boiler operated for a total of 487 hours, from 8 pm on December 10 to 2 am on December 31, during the so-called "peak period". The network includes 86 heat exchange stations (nodes 1 to 86), of which 85 heat exchange stations (except node 63) operated normally during this period.

[0076] In this embodiment, after merging pipes of the same diameter, the network consists of 186 pipes and 185 nodes, forming a single circulation block containing two independent circulations. There are 47 pipes arranged along these two circulations in the network. The water pump is installed at the heat source or heat exchange station, and no pump or check valve is set in the circulation block. The throttle valve in the circulation block remains fully open throughout the heating season.

[0077] This embodiment uses a long short-term memory (LSTM) model trained on measurement data from the 2023-2024 heating season for missing measurement data from two operating heat exchange stations, node 48 and node 36. The model captures the impact of outdoor air temperature, weekly patterns, and daily rhythms on mass flow. The mass flow of the abandoned heat exchange station (node ​​63) is zero. For the remaining heat exchange stations, any missing data for short time periods are filled by linear interpolation.

[0078] In this embodiment, the total mass flow rate of the heat exchange station is used as a reference for the water supply mass flow rate, and the correction coefficients a and b determined by the least squares method are used to adjust the mass flow rate of the heat source:

[0079] 2. main +bm peak =m sub (16)

[0080] Among them, m main and m peak represent the measured mass flow rates of the primary heat source and the peak boiler, respectively.

[0081] This embodiment uses the hourly water supply temperatures of the two heat sources and the mass flow rates of all heat exchange station nodes as input data. Figure 6 The cycle basis selected for this block is shown, where the two thinnest pipes are selected as cutoff pipes, indicated by dashed lines.

[0082] In this embodiment, the physical properties of water, including density, specific heat capacity and dynamic viscosity, are set according to 80°C conditions. Based on the average outdoor air temperature of -1.2°C, the ground temperature is set to 0°C and the thermal conductivity of the soil is 2.5W / (m·K). It is assumed that all pipes are parallel single pipes, and the heat loss on the water supply side is affected by the return water temperature. The distance between the water supply pipe and the return water pipe varies from 0.43 meters to 1.47 meters depending on the pipe diameter. When calculating the heat loss of the water supply pipe, the return water temperature is assumed to be 35°C, which is a uniform value based on the measured average return water temperature. The thermal conductivity of the pipe insulation layer ranges from 0.025 to 0.05W / (m·K), depending on the construction time. The burial depth from the ground surface to the top of the outer layer of the pipe varies from 0.7 meters to 1.1 meters, depending on the road type and pipe diameter. It is assumed that the absolute roughness of all pipes is 0.045 mm.

[0083] The temperature deviation in this embodiment is defined as the difference between the time-weighted average temperature of the simulated result and the measured value at each substation. A positive deviation means that the simulated temperature is higher than the measured value, and vice versa. These deviations can arise from a number of factors, including incorrect calibration of the heat meter, inaccurate pipeline information, or operational problems such as incorrect connection of the supply and return pipes or water leaks. The calibration method can correct the simulated temperature by directly subtracting a constant offset, effectively offsetting the simulated temperature vertically.

[0084] This embodiment selects three representative nodes: two are located near the heat source (node ​​25 and node 66), and one is located at the farthest end of the long branch (node ​​84). Their hourly temperatures are compared as follows: Figure 7 These results show that the DES model achieves very high accuracy in predicting temperature at different locations.

[0085] This example performs a computational speed test over a period of 85 days using five different models. The baseline model is the original DES model without any further optimization. The speed-up effect of eliminating redundant water fronts through tolerance threshold technology is compared. If the effect of the water front on the temperature of the simulated target node is less than the specified threshold (unit: °C), the water front will be eliminated. The temperature thresholds for the four enhanced models are set between 0.00001 and 0.01 °C. When two water fronts arrive at the end of the pipe at the same time, it is counted as only one arrival event because the event is only scheduled once. The reported computation time is the average of 10 test runs.

[0086] In this example, 99.9% of the events are water front arrival events. The most intensive water front arrival events occur in the area where substations 51-56 are located, which is the deepest branch during off-peak hours. By adopting a tolerance threshold method of 0.01°C, the total number of arrival events can be reduced by 83%. This method is particularly effective in deep branches of the network, and the arrival events in a single pipe can be reduced by up to 93%. Therefore, the average number of arrival events per pipeline per day is reduced to 85. This significant reduction is due to the elimination of redundant water fronts, which reduces the number of water fronts throughout the simulation. A threshold of 0.00001°C can halve the number of water fronts. Since eliminating water fronts also suppresses the generation of new water fronts in downstream pipes, the total number of water fronts created and eliminated will be lower after applying the tolerance threshold.

[0087] In this embodiment, Figure 8The relationship between the simulation time and the number of water fronts is shown. For the benchmark calculation, the calculation time is about 1 second. The simulation time decreases linearly with the number of water fronts. Setting the threshold to 0.01°C reduces the calculation time by 72% compared to the benchmark model, and the average error difference between the substations is ±0.005°C. In addition, there is no significant difference in the temperature pattern. Therefore, it is recommended to use the threshold model. The regression line in the figure does not pass through the origin because there is an additional overhead in the hydraulic calculation itself, which is independent of the number of water fronts. The hydraulic simulation time remains stable in the five models, with an average of about 0.07 seconds. The hydraulic calculation requires four Newton-Raphson iterations in the initial step to converge, while in the subsequent mass flow change event, only 1.8 iterations are required on average because the iteration starts from the solution of the previous step. The calculation time of the hydraulic simulation in the 0.01°C threshold model accounts for 25.7%. It can be seen that the dynamic simulation method of annular pipe network and mesh pipe network based on DES proposed in the present invention can effectively improve the simulation accuracy and efficiency.

[0088] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. It should be pointed out that a person skilled in the art can make several improvements and modifications without departing from the technical principles of the present invention, and these improvements and modifications should also be regarded as within the scope of protection of the present invention.

Claims

1. A dynamic simulation method for a ring pipe network and a mesh pipe network based on DES, characterized in that: The steps include: Step S1, developing a hydraulic model capable of simulating water flow in complex annular and meshed heating networks; Step S2, extending the dynamic thermal model from unidirectional flow to support bidirectional flow; Step S3: Integrate hydraulic and dynamic thermal calculations in the discrete event simulation (DES) framework.

2. A DES-based dynamic simulation method for a ring pipe network and a mesh pipe network according to claim 1, characterized in that: The step S1 specifically includes: S1.1 Mathematical description; Expanding from a tree network to a mesh network with one or more loops, in which hydraulic simulation and thermal simulation are decoupled. This decoupling allows the mass flow of each part of the network to be determined independently before thermal analysis; The relationship between the number of cycles M, the number of vertices V, and the number of edges E in a mesh network is: E=M+V-1(1) Solving the mass flow rate of a mesh network requires at least E independent equations; For a branched pipe network without loops, M = 0. At least V-1 equations are needed to solve the mass flow of the branched pipe network. Assuming no water leakage, each node maintains mass balance, which can be expressed as: Among them, D j is the water demand or supply of the user and production site nodes, and the internal nodes are 0; E j is a set of pipes connected to node j; m i is the mass flow in pipe i, which is positive when flowing in the normal direction and negative when flowing in the reverse direction; the sign is positive for normal inflow and negative for normal outflow; since the total mass input in the network must match the total demand ∑D j =0, thus forming V-1 linearly independent node balance equations, thus forming a definite equation system for calculating the mass flow rate of each pipeline in the tree network; Meshes containing loops require additional equations to solve for mass flow, since each loop introduces an additional edge. To solve this problem, M additional independent equations are required, which are beyond the scope of the node mass balance equation. Based on Kirchhoff's second law, which states that in any closed loop, the sum of the pressure differences along the loop must be zero, the pressure equation for each loop c is derived: Among them, △P i (m i ) indicates the flow rate m i The pressure drop in the lower pipe i is positive when flowing in the normal direction and negative when flowing in the reverse direction. The sign of each pressure term depends on whether the normal direction of the pipe is consistent with the direction of the loop. The pressure loss in the pipe is due to friction loss and local pressure loss. Since booster pumps are not required in the network, friction loss is quantified by the Darcy-Weisbach equation, where the variables include the pipe length L, pipe diameter D, water density ρ, flow velocity v and friction coefficient f: Among them, m represents the flow rate, and A represents the cross-sectional area of ​​the pipe; The friction coefficient f depends on the flow conditions: ε is the absolute roughness of the tube wall: Re=|v|Dρ / μ(6) Where μ is the dynamic viscosity of water. To simplify the calculation, the proven explicit approximate formula is selected to replace the implicit Colebrook-White equation for the friction factor of turbulent flow. In the transitional flow range, the linear interpolation method is used for processing. S1.2 Solve for mass flow in a tree network; for a tree network with only one production node, assume that the normal direction of the pipe is consistent with the only direction of flow, that is, in the supply network, the flow is from the production node to the user node, and in the return network, the flow is from the user node back to the production node; the mass flow of all pipes is calculated in a bottom-up manner, starting from the user node, summarizing the flow at the confluence point until it reaches the production node; topological sorting can be used to predetermine the order in which the pipe flows are calculated; since the tree topology remains static throughout the simulation, this order will always remain unchanged after initialization; S1.3 solves the mass flow in the mesh network; after the initial water distribution is established, the mass flow of the boundary nodes is determined based on the flow of the connected non-loop edges, where positive values ​​indicate flow into the area and negative values ​​indicate flow out of the area; specifically, the mass flow of the boundary nodes is obtained by calculating the flow of the non-loop edges connected to these nodes; for each boundary node, if the flow of the connected non-loop edge is positive, it means that water flows into the block where the node is located; and if the flow is negative, it means that water flows out of the node; in this way, the flow of the boundary nodes can accurately reflect the input and output of the system, ensuring the conservation and balance of mass of the entire network; To handle the nonlinear subproblems in each loop block, the Newton-Raphson algorithm is employed with the following iterative steps: (1) Calculate the pressure drop △P of each pipeline i (m i ) and its derivative d△P i (m i ) / dm i ; (2) accumulating the pressure drop in each cycle according to the above method, and accumulating its derivative in the same way; (3) Evaluate the pressure difference imbalance in each cycle; if the pressure difference error P of all cycles in the block j If it is zero within the tolerance range, the system considers that it has converged; if not, the updated mass flow step length △m is calculated using the following method (k) ; Δm (k) =m (k+1) -m (k) =-J -1 P(7) Among them, m (k) represents the mass flow step size of the kth iteration, J represents the Jacobian matrix, and P represents the pressure error; (4) Update the mass flow of all edges in the block according to the updated flow of the truncated edges and the flow initially set at the boundary nodes; △m (k) represents the update step size in the kth iteration; a positive step size indicates an increase in mass flow along the default circulation direction; the Jacobian matrix J is: Among them, P n Indicates pressure error; m n Indicates mass flow rate.

3. A DES-based ring pipe network and mesh pipe network dynamic simulation method according to claim 1, characterized in that: The step S2 specifically includes: S2.1 calculates the bidirectional flow in the pipe; the temperature of water particles in the pipe decreases exponentially with the travel time τ due to heat conduction to the surroundings: Among them, T in is the temperature of the water when it enters the pipe, T out is the temperature of the water as it exits the pipe; k1 and k2 are coefficients for single, parallel, and double pipes; the Lagrangian DES method involves tracking the motion of the water front, which is an infinitely thin section of water moving inside the pipe; the water front acts as a dynamic sampling point that records various inlet parameters; the water front is defined at the breakpoints in the temperature profile; numerical analysis shows that by identifying these breakpoints and assuming a linear transition of temperature between adjacent points, the water temperature can be predicted with high accuracy; the bidirectional flow model follows the same principle, but the flow direction is taken into account when calculating the travel time; The travel time of each particle is determined by the arrival time t out and entry time t in =Calculate the difference between |v1| and |v2|; Assume that: at t1, the current inlet temperature changes, generating the first and only water front F1, and scheduling its arrival event according to the flow rate |v1|; at t2, when the flow stops, water fronts F2 and F3 are generated at the inlet and outlet, and the arrival event is cancelled; at t3, when the flow direction reverses, water front F2 reaches the left end of the current target node, and schedules the arrival event of F1 according to the flow rate |v3|; at the same time, water front F4 is generated at the right end of the current inlet; fronts F3 and F4 are extremely close, indicating a step change in temperature; at t4, front F1 reaches the current outlet, and F3 becomes the next arriving front, and its arrival time is scheduled according to |v3|; the travel distance of particle F1 in the time period from t1 to t2 is equal to the travel distance in the time period from t3 to t4; the entry time and travel time of F1 can be derived from this principle and further described: in, represents the entry time of particle F2; represents the entry time of particle F1; represents the arrival time of particle F1; represents the arrival time of particle F2; Δt represents the time interval; v1 represents the flow velocity at the time of entry; v3 represents the tassel at the time of arrival; represents the travel time of particle F1; represents the travel time of particle F2; Likewise, the travel time and entry time of any arbitrary water front can be calculated from the previous arrival time, the corresponding entry time, the current velocity, and the corresponding entry velocity; Taking into account the sign of the velocity, the equations for calculating the travel time and entry time of a particle in a bidirectional flow agree with the equations for a unidirectional flow and are given by: in, represents the entry time after the kth iteration; represents the arrival time after the kth iteration; v (k) represents the flow rate after the kth iteration; represents the flow rate at the time of entry after the kth iteration; τ (k) represents the travel time after the kth iteration; t (k) represents the time after the kth iteration; There are four states of water flow: positive flow, +0, -0 and negative flow; positive flow and negative flow refer to the flow of water in the normal direction and opposite direction along the pipe respectively; +0 and -0 refer to the state where the water stops flowing after positive flow or negative flow; the travel time function of water particles through the pipe is piecewise linear and resets in two cases: when the flow direction is reversed or the flow is temporarily stopped; specifically, the following situations are included: when the flow direction changes, the travel time is reset to zero; during the flow stop, the travel time is undefined; when the flow resumes, the reset travel time value will depend on the relative direction of the flow before and after the stop; if the flow direction is reversed, the travel time is set to the duration of the stop; if the flow returns to its original direction, the travel time will be extended by the duration of the stopped flow; S2.2 Analyze the effect of flow direction changes on adjacent pipes; changes in flow direction within a pipe can significantly affect the inlet temperature of adjacent pipes, so additional water fronts need to be introduced to monitor these discontinuities; the mixing temperature at node j is the outlet temperature of each pipe It is calculated by taking the weighted average of: in Represents the set of pipes that flow into node j.

4. A DES-based dynamic simulation method for a ring pipe network and a mesh pipe network according to claim 1, characterized in that: The step S3 specifically includes: The DES model represents parallel processes as a series of events arranged in time order, which represent changes in the system state; during the simulation process, events can be scheduled, rescheduled, or canceled; these events are uniformly managed in an event queue, which is organized according to the scheduled activation time of the events; the discrete event simulation process involves iteratively removing the earliest event from the queue and executing the event routines associated with it until the event queue is empty; the dynamic scheduling capability of discrete event simulation makes it particularly suitable for simulations that require variable time steps; The specific process is as follows: Initialize the event queue with mass flow change event and inlet change event at the start time; if the event queue is empty, the process ends; if it is not empty, remove the first event from the event queue and process it; If the first event is a mass flow change event, the mass flow is solved through the hydraulic model, and the node temperature is updated based on the previous mass flow. Each pipe is traversed to update the flow direction, the frontier double-ended queue is maintained, a new water flow front is inserted, and the arrival time of the event is updated; the node temperature is updated according to the updated mass flow, and if the node temperature changes, an additional water flow front is inserted; if there is a next change, the mass flow change event is rescheduled; the event queue is maintained by re-minimum heap, and finally it is determined again whether the event queue is empty, and the cycle repeats; If the first event is a water front arrival event, the outlet temperature is updated and the arriving water front is removed from the double-ended queue; if the double-ended queue is empty, the inlet temperature change event of all downstream objects is triggered; if the double-ended queue is not empty, the arrival time needs to be rescheduled before triggering the inlet temperature change event of all downstream objects; finally, it is determined again whether the event queue is empty, and the cycle repeats; If the first event is an inlet temperature change event, a new water front is created in the pipeline. When the double-ended queue is not empty, a check is made directly to see if the event queue is empty. If the double-ended queue is empty, an arrival event is scheduled, and then a check is made to see if the event queue is empty, and the cycle repeats.

Images