Branching pipe network thermal dynamic modeling method based on transfer function

Abstract

The application provides a kind of branch pipe network thermal dynamic modeling method based on transfer function, belongs to the technical field of electric digital data processing.The application first establishes the advection diffusion-heat loss control equation of single pipeline, and deduces the single pipe frequency domain transfer function through Laplace transform;Then, combined with theoretical analysis and high-precision numerical simulation, the turbulent temperature diffusion coefficient α and temperature loss coefficient γ are accurately calibrated and corrected, which significantly improves the physical consistency of the model;For the topology structure of branch pipe network, it is decomposed into multiple independent branches according to the heat user load node, and the frequency domain expression of the overall pipe network return water temperature is constructed based on the flow weighting principle;Finally, the time domain response is obtained by inverse Laplace transform.The application effectively solves the contradiction between the insufficient accuracy of traditional node method and the complexity of high-precision numerical method, balances the modeling accuracy and calculation efficiency, and is suitable for real-time optimization scheduling of heating pipe network in integrated energy system.

Description

Technical Field

[0001] This invention belongs to the field of electrical digital data processing technology, specifically relating to a method for thermal dynamic modeling of branched pipe networks based on transfer functions. Background Technology

[0002] Against the backdrop of the low-carbon transformation of energy systems, district heating networks, as the core carrier of integrated energy systems (IES), achieve cascaded energy utilization through the coupling of multiple energy flows including electricity, heat, and gas, becoming a key pathway to improve system flexibility and reduce carbon emissions. [1] .

[0003] Thermal transients refer to the dynamic process of temperature propagation in a heating network caused by heat source regulation, load changes, or start-up and shutdown operations, with a time scale often ranging from tens of minutes to several hours. [5][6] In the dynamic adjustment of heating pipe networks, thermal transient phenomena are prevalent. These phenomena are essentially the spatiotemporal evolution of temperature waves propagating along the pipes and have multiple impacts on energy dispatch. Adjusting the operating state of the pipe network will trigger two key phenomena. [7][8] (1) Temperature front diffusion: When temperature fluctuations caused by heat source regulation or load changes propagate in the pipeline, due to fluid turbulent mixing (causing the step temperature signal to diffuse during transmission) and radial heat transfer through the pipe wall, the temperature distribution exhibits a gradient smoothing characteristic, rather than an ideal step change (the user-side temperature response deviates from the ideal step); (2) Dynamic heat loss: Pipeline radial heat exchange (continuous heat exchange between the pipeline and the environment) is coupled with fluid axial transport, forming unsteady heat loss, resulting in irreversible energy attenuation during transmission. In the dynamic regulation process of the heating network, this thermal transient phenomenon is common. The superposition of the two key phenomena mentioned above causes the user-side temperature response to lag significantly behind the heat source regulation, and the response curve deviates from the ideal model. If the thermal transient effect is ignored, it may lead to unmet heat user demand, equipment overload, or system safety risks.

[0004] Branched pipe networks are widely used as the mainstream topology for urban district heating and industrial steam transmission due to their advantages of short total pipeline length and low construction cost. However, the flow distribution coupling effect caused by their tree-like topology and the thermal inertia of long pipelines pose unique challenges to thermal dynamic modeling: the heat flow propagation path is nonlinearly distributed, the flow distribution of the main pipe and branch pipes affects each other, and heat source regulation or load changes need to be transmitted through multiple branches, resulting in the thermal response process having both axial delay and radial diffusion characteristics. Although thermal transients are a common characteristic of heating pipe networks, the topological characteristics of branched pipe networks make their thermal dynamic behavior more complex, requiring higher modeling accuracy and efficiency.

[0005] Existing methods for modeling and solving heating network systems mainly include the following categories:

[0006] (1) Steady-state thermal model and lumped model

[0007] Steady-state thermal models for district heating (DH) networks have been widely used in the simulation and analysis of DH networks. [9]

[10] and operation optimization

[11]

[12] However, steady-state models cannot capture the thermal transients of DH networks.

[0008] Kunz's early research on thermal transient models of DH pipe networks focused on lumped modeling of DH pipes.

[13] In this model, each pipe is set as a lumped heat capacity, and the water temperature is set as the average of the inlet and outlet temperatures. The drawback of this method is that a lumped model with only one heat capacity cannot capture the distributed parameter characteristics of a DH network.

[0009] Chinese patent CN111310343B provides a method for modeling the thermal circuit of a heating network for integrated energy system scheduling. Based on the energy conservation equation of the thermal pipes in the heating network, a lumped parameter model of the thermal pipes in the frequency domain is established. Based on this model, the thermal branch characteristics of the heating network are established. Thermal topological constraints of the heating network are established for the confluence and divergence processes at each node of the heating network. A thermal circuit model of the heating network is established by integrating the thermal branch characteristics and thermal topological constraints. The shortcomings of this method are that it does not correct key parameters (turbulent temperature diffusion coefficient α, heat loss coefficient γ) and does not consider temperature front diffusion and dynamic heat loss during thermal transient processes.

[0010] (2) Node method

[0011] A further development of techniques for simulating thermal transients in DH networks is the nodal method.

[14] Because the velocity in the DH pipeline is continuously changing, the time delay of the temperature wave between the two nodes is also constantly changing. The principle of the nodal method is to track the time delay by calculating the time required for water mass to move from one node to another. Based on the temperature history time series of the pipeline inlet temperature, the outlet temperature can be calculated according to the changes in heat loss and heat capacity of the pipeline. In this way, the new temperature of all nodes in the network at each time step can be calculated.

[0012] gabrielaitienkv et al. applied the nodal method to the DH system, focusing on the dynamic performance of temperature distribution distortion.

[15]

[16] The key advantage of the nodal method is its low time complexity, making it highly suitable for online optimization. Expressing nodal temperatures through linear equations simplifies the description of the dynamic characteristics of the supply and return heating network, leading to its widespread use in integrated energy system optimization modeling. [2][4]

[17] While the nodal method provides an effective approach for thermal transient modeling, it still suffers from two main drawbacks: (a) the nodal method does not address the propagation of temperature waves along the pipe, as the diffusion and smearing of temperature waves are predictable. (b) the nodal method only calculates the nodal temperatures, not the temperature distribution along the pipe.

[18] The accuracy is limited when the node method is applied to mesh DH networks or DH networks with multiple heat sources.

[0013] Chinese patent CN110135013B provides an integrated method for constructing a dynamic thermal model of a heating network. This method utilizes a node-based approach to establish a dynamic thermal model of pipes and nodes. It integrates all nodes and pipes to build the dynamic thermal model of the heating network by constructing an inflow pipe matrix, an inflow pipe starting node matrix, an inflow pipe quantity matrix, and a pipe starting node matrix. The drawback of this method is that for large-scale urban heating networks (containing tens of thousands of nodes and pipes), the matrices generated by this method become extremely large and sparse. Solving such a high-dimensional system of differential-algebraic equations requires enormous computational resources.

[0014] (3) Numerical methods

[0015] Numerical methods can provide accurate predictions and sufficient information on thermal transients. However, a drawback of numerical methods is their time consumption, especially when applied to the online operational optimization of large-scale DH networks.

[0016] The element method based on finite element analysis discretizes space and time, decomposing partial differential equations into several algebraic problems that can be solved using numerical schemes.

[19]

[20]

[21]

[22]

[23] This method typically offers high accuracy, but careful parameter tuning is required to avoid numerical diffusion and excessive computation time.

[24]

[25] .

[0017] The unique aspect of the Lagrange method lies in the observer's perspective; the observer moves with the flow rather than remaining in a fixed position.

[26]

[27] It can overcome the limitations of Euler-based discretization.

[26] Numerical diffusion with inherent correlation. To improve the accuracy and speed of the solution, scholars have conducted extensive research on solution methods. Stevanovic et al. proposed a highly efficient numerical method for predicting thermal transients in DH networks. [5]

[28] Employing efficient numerical methods

[28] The flow rate of the DH pipeline was calculated using a third-order Lagrange interpolation polynomial, and the energy equation was solved using the method of characteristics. The predicted temperature front propagation showed good agreement with measured data. Wang et al. then proposed a new third-order numerical solution to solve the hydraulic-thermal coupling model.

[22] This method not only preserves the sharp temperature fronts during heat propagation but also achieves good computational accuracy even with coarse meshes. However, higher-order schemes typically require more stringent stability conditions, which directly means that smaller time steps must be used to maintain computational stability. In dynamic simulations, this can lead to a sharp increase in the total number of computation steps, significantly extending the simulation time and making it difficult to meet the demands of real-time applications.

[0018] Chinese patent CN112100778B provides a scheduling method and system for an integrated energy system. This patent uses differential processing to obtain a heat network pipeline transmission model from the one-dimensional thermal energy flow conservation equation. Furthermore, by transforming the difficult-to-solve one-dimensional thermal energy flow conservation partial differential equation into an ordinary differential equation through differential algorithms, it effectively reduces the computational load and simplifies the calculation of the dynamic outlet water temperature of the heat network pipeline. However, its complex differential processing leads to an extremely high computational load.

[0019] (4) Other modeling methods

[0020] Besides nodal methods and finite difference numerical methods, scholars have also tried to use other methods to model and solve heating systems. For example, Wang used an autoregressive moving average (ARMA) model to describe the temperature dynamics of heating systems.

[29] Xu et al. used a quasi-dynamic pipe model to describe the transport delay effect, that is, they used the inlet temperature and heat loss model to estimate the outlet temperature. [3] Analytical solutions require that the relevant equations be explicitly solvable so that the system of equations and its solutions can be represented in matrix form. This method is significantly more computationally efficient than numerical methods, and analytical solutions hold great potential in the simulation of integrated energy systems (IES) because they facilitate the integration of district heating networks with other networks within the system.

[30] However, these models can only handle quasi-dynamic or steady-state simulations. Gou used the RC model to establish a heat flow model for the DHS, modeling the transmission delay and heat loss of the pipe using both heat potential and heat flow sources.

[31] The heat flow model accurately describes the heat migration and transfer process, and its analogy to electrical conductivity facilitates the application of mature electrical simulation software for DHS simulation. Denarie proposed a novel numerical method to simulate heat transfer in long pipes.

[32] This model employs a characterization method to partition the heat capacity of water between the turbulent core and the boundary layer. Compared to the finite volume method and the nodal method, this model has lower computational cost, provides more accurate results, and does not require artificial smoothing of temperature waves. K. Sartor demonstrated the limitations of the one-dimensional finite volume method in pipe heat transfer simulation and proposed an alternative method that considers heat loss and pipe inertia.

[33] The results show that the thermal inertia of the pipeline has a significant impact on the temperature distribution at the outlet.

[0021] As mentioned above, it can be seen that among the current modeling methods coupled with integrated energy systems, there are almost no pipeline modeling methods that can take into account the thermal transient process of the pipeline network. For pipeline systems with long transmission times, time delay can play an important role in flexible control. However, at this time, radial heat transfer occurs drastically, and thermal transient changes become non-negligible. Therefore, the coupling of IES optimization considering pipeline characteristics faces two huge gaps: (1) The non-negligibility of thermal transient phenomena: Pipelines can provide great energy flexibility in energy system optimization, but frequent heat source actions make thermal transient phenomena non-negligible. Ignoring thermal transient phenomena will lead to unmet heat user needs; (2) The complexity of solving thermal transient processes: Large pipeline networks contain hundreds of pipes. High-precision solution methods require several hours to solve a single transient, which cannot meet the needs of online optimization.

[0022] Summarizing existing technologies, we can identify the core contradiction in current pipeline network modeling in energy systems: increasing the range of pipeline network movements to enhance flexibility exacerbates thermal transient effects, while existing models cannot balance accuracy and efficiency.

[0023] In summary, existing heating network modeling methods have significant limitations in dealing with thermal transient analysis: (1) Traditional node methods calculate time delays through mass tracking, but can only characterize the macroscopic delay of axial heat flow transmission and cannot analyze the microscopic mechanism of temperature front diffusion, resulting in a decrease in the prediction accuracy of the model in complex topologies; (2) High-precision numerical methods (such as the third-order finite difference method) can analyze the spatiotemporal evolution of thermal transients, but require discretization of the network for solution, and the computational complexity increases exponentially with the network size, making it difficult to meet the real-time requirements of online optimization; (3) Existing technologies mostly focus on simple networks or steady-state models, and do not adequately support modular decomposition and parallel computing of branched network topologies, resulting in limited model scalability and inability to be directly applied to practical scenarios; (4) Existing technologies do not fully consider the impact of temperature front diffusion and dynamic heat loss on the thermal dynamic response of the network, resulting in reduced accuracy; high-precision models have high computational complexity and are difficult to meet the needs of online optimization; low-order models are simple to calculate, but lack accuracy.

[0024] In response, the transfer function method is widely used in the field of automatic control to quickly solve system responses, providing a new approach for modeling virtual energy storage in pipeline networks and helping to improve the response speed and accuracy of scheduling models: a frequency domain transfer function method is proposed.

[34] This method is used for thermal characterization and energy flexibility design in areas with circulating radiative heating / cooling systems and significant thermal mass. Simultaneously, the method of establishing system transfer functions has also been applied in the field of heat transfer modeling with good results, providing a new approach to solving the real-time state of pipe networks under heat source operation.

[35] . Summary of the Invention

[0025] This invention is based on the above background technology and aims to provide a method for thermal dynamic modeling of branched pipe networks based on transfer functions.

[0026] This invention provides a method for thermal dynamic modeling of branched pipe networks based on transfer functions. The branched pipe network includes multiple independent branches and has the following characteristics: S10, establishing the advection diffusion-heat loss control equation for a single pipe. Later on time t Perform a Laplace transform and derive the single-tube frequency domain transfer function: ,in, Indicates the axial position of the pipe x Location, Time t The difference between the fluid temperature and the ambient temperature at that time. V Indicates the average flow velocity of the fluid. Represents the convection term. This represents the heat loss term. c Indicates the temperature loss coefficient. Indicates the diffusion term. α Represents the turbulent temperature diffusion coefficient, customizable. , Represents the Laplace transform. s S20, considering the thermal inertia of the pipe wall and the inner diameter, the temperature loss coefficient is determined through a finite difference numerical simulation experiment. c Perform precise calibration and correction; S30, apply Reynolds averaging to the advection diffusion-heat loss control equation, and determine the turbulent temperature diffusion coefficient by combining the gradient diffusion assumption and the Prandtl number of turbulence. α The associated physical quantities were then used to simulate the pipeline using the SST turbulent kinetic energy-specific dissipation rate turbulence model. The internal flow field information was compressed into one-dimensional information using a velocity-weighted averaging method. Finally, the accurate turbulent temperature diffusion coefficient was obtained by fitting the one-dimensional information with the associated physical quantities. α S40, the branched pipe network is decomposed into multiple independent branches according to the heat user load nodes, and the independent branches are decomposed into multiple pipes, based on the calibrated and corrected temperature loss coefficient. c and the turbulent temperature diffusivity α S50, and construct the frequency domain expression of the return water temperature of the branched pipe network based on the flow weighting principle; S50, perform inverse Laplace transform on the frequency domain expression of the return water temperature of the branched pipe network to obtain the return water temperature response in the time domain so as to realize the thermal dynamic modeling of the branched pipe network.

[0027] The heat dynamics modeling method for branched pipe networks based on transfer functions provided by this invention may also have the following feature: step S10 includes the following sub-steps: S11, establishing the advection diffusion-heat loss control equation for a single pipe: S12, Custom Then, the Laplace transform of the advection diffusion-heat loss control equation is used to obtain the ordinary differential equation in the frequency domain: S13, Solve the ordinary differential equation to obtain the single-tube frequency domain transfer function: .

[0028] The method for thermal dynamic modeling of branched pipe networks based on transfer functions provided by this invention may also have the following feature: wherein, in step S12, it is assumed that the solution form of the ordinary differential equation is: , l Let be a constant, and the characteristic equation be: The final root is: .

[0029] The heat dynamics modeling method for dendritic networks based on transfer functions provided by this invention may also have the following feature: step S20 includes the following sub-step: S21, associating the heat loss term based on heat transfer theory. Theoretical value and introduce the inner diameter of the pipe. D , and thus c Theoretical value: In the above formula, r Indicates fluid density, c p Indicates the specific heat capacity of a fluid. R This represents the thermal resistance per unit length of the pipe. and D and R Related; S22, define a custom energy retention rate: S23. After establishing a finite difference model of pipe heat loss considering the thermal inertia of the pipe wall as a reference model, a temperature pulse is applied as input at the pipe inlet of the reference model. The temperature response at the pipe outlet is calculated by running the reference model. The energy retention rate of the reference model is obtained by spatial integration of the temperature response based on the user-defined energy retention rate, and denoted as the numerical simulation solution. S24, using the analytical solution of the advection diffusion-heat loss control equation. Solving for the custom energy retention rate yields the energy retention rate of the theoretical model, denoted as . ,in, Indicates the pulse half-width. This indicates the initial height of the temperature pulse. S25 represents the error function; for several groups D and R The combination of will have different D The combinations are divided into different classes, and the fitting method is used to find... The value makes and The best match is used to determine the correction factor. τ is the coefficient obtained from the fitting. k D To take into account the inner diameter of the pipe D The correction factor is used to finally obtain: .

[0030] The heat dynamics modeling method for branched pipe networks based on transfer functions provided by this invention may also have the following feature: step S21 includes the following sub-step: S211, the heat loss per unit length of the pipe is given by Newton's law of cooling: S212, based on the heat capacity per unit length of the fluid. C = pAc p The relationship between the rate of decrease in fluid temperature due to heat loss per unit length is determined as follows: ,in, A This represents the cross-sectional area of ​​the fluid inside the pipe; S213, derived from heat transfer theory: S214, the heat loss term This represents the contribution of temperature decay, therefore its theoretical value is expressed as: S215, after sorting, we get: S216, Introducing the inner diameter of the pipe D Substitution A From: .

[0031] The thermal dynamic modeling method for branched pipe networks based on transfer functions provided by this invention may also have the following feature: wherein, in step S23, in the reference model... x Apply a persistent [condition] at =0. Q The temperature in seconds is T 输入 The rectangular temperature pulse is taken as input, assuming ,but: Then, the boundary value problem in the above equation is changed to an initial value problem, that is, the action of the heat source in time can be regarded as the generation of a rectangular temperature pulse in space: In step S24, the convection-diffusion-attenuation equation within the pipe is given using the Green's function method. Green's function in the form of point source solutions: The solution to the rectangular temperature pulse is the convolution of the Green's function: ,in, Let the integral variable be represented. Substitute the Green's function into the convolution formula, and obtain the analytical solution of the advection diffusion-heat loss control equation through variable substitution and integration: .

[0032] The heat dynamic modeling method for branched pipe networks based on transfer functions provided by this invention may also have the following feature: wherein step S30 includes the following sub-steps: S31, ... V Decomposed into ,Will Decomposed into ,in, and It is the pulsation quantity in turbulent flow. and It is the Reynolds average quantity in turbulent flow. , S32, Applying Reynolds average to the advection diffusion-heat loss governing equation: S33, Turbulent heat flux based on the gradient diffusion assumption Modeling Substituting this into the advection diffusion-heat loss control equation after applying Reynolds average, we obtain: ,in, express α The theoretical value, α 水 The thermal conductivity of water is represented by... α 水 The order of magnitude is extremely small, so use express α Theoretical value; S34, using the Prandtl number of turbulence. Will With turbulent viscosity Related: Where β represents turbulent kinetic energy, oh This represents the specific dissipation rate, since and ,so , D S35 represents the inner diameter of the pipe; the SST turbulent kinetic energy-specific dissipation rate turbulence model is used for... D and V After performing computational fluid dynamics simulation on the combined pipes to calculate its internal complete flow field, the internal complete flow field is simplified in one dimension using the velocity-weighted averaging method to obtain: ,in, A This represents the cross-sectional area of ​​the fluid inside the pipe. The cross-sectional average velocity, r Represents radial coordinates, Represents velocity distribution. Represents radial coordinates r Turbulent temperature diffusivity at point S36; select several different D and V By combining the elements and repeating step S35, several different combinations are obtained. Data points, and their relationship with D × V Perform linear fitting to determine the correction factor k α ,but: .

[0033] The heat dynamic modeling method for branched pipe networks based on transfer functions provided by this invention may also have the following feature: step S40 includes the following sub-steps: S41, decomposing the branched pipe network according to heat user load nodes. n Each independent branch contains a complete flow path from the heat source to the corresponding load node and back to the heat source; S42, the components of the independent branch are divided into upstream segments. Load nodes and downlink segments symmetrical to the uplink segments. , among which, in the first k Temperature drop at the load node of each of the independent branches , k =1,2,3,……, n Each of the aforementioned independent branches is 2 m It consists of several pipes, including the upstream section. m The pipeline and the downstream section mentioned above m The frequency domain transfer function of each of the aforementioned pipes is as follows: , , ,in, i The number of the pipe representing a single, independent branch. i =1,2,3,……, m , L i , G i , α i , c i 、V i These are, respectively, the pipe length corresponding to the pipeline of the independent branch, the single-pipe frequency domain transfer function, and the corrected... αCorrected c as well as V S43, the total return water temperature is recorded as Then it flows into the first k The water supply from the aforementioned independent branch flows through the upstream section. The load node and the downlink segment The next k The return water temperature of each of the aforementioned independent branches is: ,in, S44: Calculate the flow weight of each independent branch according to the flow rate ratio, representing the change in heat source temperature; w k Later Result obtained by flow weighting: S45, Decoupling and rearranging the formula in step S44 yields the frequency domain expression for the return water temperature of the branched pipe network: Among them, the heat source transfer function , No. k Load transfer function of each load node .

[0034] The heat dynamics modeling method for branched pipe networks based on transfer functions provided by this invention may also have the following feature: wherein, in step S44, the flow weight of each of the independent branches... , Q k Indicates the first k The flow rate of an independent branch, Q i Indicates the first i The flow rate of the pipe.

[0035] The method for thermal dynamic modeling of branched pipe networks based on transfer function provided by the present invention may also have the following feature: wherein, in step S50, the inverse fast Fourier transform algorithm is used to perform the inverse Laplace transform.

[0036] The beneficial effects of this invention are:

[0037] (1) Transfer function modeling method based on parameter correction: Based on the one-dimensional advection diffusion-heat loss control equation, the single-tube frequency domain transfer function model is derived through Laplace transform. Combined with theoretical analysis and high-precision numerical simulation, the corrected turbulent temperature diffusion coefficient is calibrated. α ) and heat loss coefficient ( c The physical correlation expression of the pipeline material properties, fluid velocity and ambient temperature is used for the first time. This expression significantly improves the physical consistency of thermal transient modeling of long pipelines and reduces the temperature distribution prediction error by more than 50% compared with the traditional nodal method.

[0038] (2) Modular decomposition and parallel computing strategy: For branched pipe network topology, a modular decomposition strategy is proposed to decouple the complex pipe network into independent branch combinations and construct a network model using flow weights and branch transfer functions. This strategy supports parallel computing, avoids solving large simultaneous equations, and reduces the single transient simulation time by 80% compared with the third-order finite difference method, while maintaining millisecond-level response accuracy. Attached Figure Description

[0039] Figure 1 This is a flowchart of a method for thermal dynamic modeling of branched pipe networks based on transfer functions in this embodiment.

[0040] Figure 2 This is the relationship between the temperature loss coefficient and the thermal resistance per unit length for different pipe diameters in the test examples of this invention.

[0041] Figure 3 This is the dimensionless flow velocity distribution curve as a function of dimensionless radius in the CFD simulation results of the test example of this invention.

[0042] Figure 4 This is the distribution curve of thermal diffusivity as a function of dimensionless radius in the CFD simulation results of the test example of this invention.

[0043] Figure 5 This is the fitting result of the turbulent temperature diffusion coefficient in the test example of the present invention.

[0044] Figure 6 The test examples of this invention use experimental test data and literature. [6] The single-pipe frequency domain transfer function of the embodiment is verified by testing and simulating the temperature changes at the pipe inlet and outlet.

[0045] Figure 7 This is a comparison of the pipeline network modeling method in the test examples of this invention. Detailed Implementation

[0046] To make the technical means, creative features, objectives and effects of this invention easier to understand, the following embodiments, in conjunction with the accompanying drawings, specifically illustrate a method for thermal dynamic modeling of branched pipe networks based on transfer functions.

[0047] Example

[0048] Before elaborating on the heat dynamics modeling method for dendritic networks based on transfer functions in this embodiment, the following basic assumptions are established for heat transfer dendritic networks:

[0049] (1) The fluid flow is an incompressible, uniform and fully developed steady-state flow.

[0050] (2) The temperature and velocity distribution inside the pipe is axially symmetric.

[0051] (3) The thermal properties of the pipe material are constant and independent of time.

[0052] (4) The hydraulic conditions of the heat transfer branch network are stable, and the heat carrier is a single-phase fluid (water) filling the pipes. Extreme cases such as water evaporation or two-phase flow are ignored.

[0053] (5) The heat source regulation mode is quality regulation, that is, adjusting the outlet water temperature of the heat source under stable hydraulic conditions.

[0054] This embodiment provides a method for thermal dynamic modeling of branched pipe networks based on transfer functions, wherein the branched pipe network includes multiple independent branches, and each independent branch includes multiple pipes.

[0055] Figure 1 This is a flowchart of a method for thermal dynamic modeling of branched pipe networks based on transfer functions, as described in this embodiment. Figure 1 As shown, the thermal dynamic modeling method for dendritic networks based on transfer functions in this embodiment includes the following steps S10~S50:

[0056] S10, Derive the single-transistor frequency domain transfer function, including the following sub-steps S11~S13:

[0057] S11, Ideally, the temperature pulse undergoes diffusion and deformation during propagation, similar to the transmission and diffusion of a pulse signal. To explore its dynamic laws, the pipe is considered as a one-dimensional space, and a temperature pulse thermal convection-diffusion equation considering heat loss is established.

[0058] In steady turbulent flow, the dynamics of the thermal field are governed by the following advection diffusion-heat loss control equation:

[0059] .

[0060] in, Indicates the axial position of the pipe x Location, Time t The difference between the fluid temperature and the ambient temperature at that time. V Indicates the average flow velocity of the fluid. Represents the convection term. This represents the heat loss term. c Indicates the temperature loss coefficient. Indicates the diffusion term. α This represents the turbulent temperature diffusion coefficient.

[0061] S12, Custom Then, by performing a Laplace transform on the advection diffusion-heat loss governing equation, we obtain the ordinary differential equation in the frequency domain:

[0062] .

[0063] in, Represents the Laplace transform. s Represents complex frequency variable

[0064] S13, Assume the solution to the ordinary differential equation in step S12 is in the form of... The characteristic equation is Then its final root is ,in, l It is a constant.

[0065] By selecting a physically appropriate attenuation root, the single-tube frequency domain transfer function is obtained:

[0066] .

[0067] From the above equation, we can see that the core parameters of the single-tube frequency domain transfer function are: α and c .literature [9] Point out α Mainly related to the average fluid velocity V and the inner diameter of the pipe D Related; c Except for the average fluid velocity V and the inner diameter of the pipe D In addition to being related to the working fluid properties (heat transfer coefficient) k Kinematic viscosity of working fluid n (Related to) Literature [6] Based on the phenomenology of turbulence, the following equation is derived: , .

[0068] However, despite the literature [6] right α and c Rigorous derivation was performed, but its estimation based on turbulence phenomenology lacked sufficient experimental verification. Although the derived relationships can reflect the influence trend of the parameters, the calculated results do not fit the experimental data well.

[0069] Therefore, it is necessary to provide this α and c A more accurate calibration method is sought for the two parameters, as shown in steps S20 and S30 below.

[0070] S20, for temperature loss coefficient c Perform precise calibration and correction (this step assumes the ambient temperature is constant and uniform), including the following sub-steps S21~S25:

[0071] S21, based on heat transfer theory, related heat loss terms Theoretical value And introduce the inner diameter of the pipe.D , and thus c Theoretical value Specifically, it includes the following sub-steps S211~S216:

[0072] S211, the heat loss per unit length of the pipe is given by Newton's law of cooling:

[0073] .

[0074] in, R The thermal resistance per unit length of a pipe depends on the pipe geometry (such as pipe diameter). (and insulating materials.) Heat loss causes the fluid temperature to drop, and its energy balance is based on the heat capacity per unit length of fluid.

[0075] S212, based on the heat capacity per unit length of the fluid C = pAc p Determine heat loss per unit length Q loss The relationship that leads to the rate of decrease in fluid temperature is as follows:

[0076] .

[0077] in, A This represents the cross-sectional area of ​​the fluid inside the pipe. r Indicates fluid density, c p This indicates the specific heat capacity of the fluid.

[0078] S213, derived from heat transfer theory, is as follows:

[0079] .

[0080] S214, Now, it is necessary to correlate the heat loss term of the advection diffusion-heat loss control equation in step S11. In the advection diffusion-heat loss governing equation, the heat loss term is: This represents the contribution of temperature decay, which can be expressed as:

[0081] .

[0082] S215, because heat transfer in the branched pipe network is always stable, ,so:

[0083] .

[0084] S216, in order to be applicable to actual engineering simulations, introduces the inner diameter of the pipe. D Substitute into step S215 A From:

[0085] .

[0086] The above step S216 gives c Theoretical value The expression for heat loss and its key influencing factors are discussed. However, in actual branched pipe networks, factors such as the thermal inertia of the pipe wall cause the theoretical expression to fail to accurately reflect actual heat loss.

[33] For verification and correction c An accurate baseline model for heat loss needs to be established. Subsequent iterations will utilize the literature...

[33] The proposed pipeline heat loss model (considering the influence of pipe wall thermal inertia, with good agreement between simulation and experimental results) serves as a correction benchmark, as detailed below:

[0087] S22, define a custom energy retention rate:

[0088] .

[0089] S23, establish a finite difference model of pipe heat loss considering the thermal inertia of the pipe wall as the benchmark model.

[36] Then, a temperature pulse is applied as input at the pipe inlet of the benchmark model, and the temperature response at the pipe outlet is calculated by running the benchmark model, based on a custom energy retention rate. The energy retention rate of the baseline model is obtained by spatial integration of the temperature response, denoted as the numerical simulation solution. .

[0090] Specifically, in this embodiment, in the baseline model x Apply a persistent [condition] at =0. Q The temperature in seconds is T 输入 A rectangular temperature pulse is used as input, assuming an initial temperature ,but:

[0091] .

[0092] However, since the analytical solution for the time boundary value is too complex, the boundary value problem in the above equation is changed to an initial value problem, that is, the action of the heat source in time can be regarded as the generation of a rectangular temperature pulse in space:

[0093] ,in, This represents the initial height of the temperature pulse. It is the pulse half-width.

[0094] S24, Calculate the energy retention rate of the theoretical model. This includes the following sub-steps S241~S245:

[0095] S241, using Green's function method, derives the convection-diffusion-attenuation equations within the pipe. Green's function in the form of point source solutions: .

[0096] S242, the solution for the rectangular temperature pulse is the convolution of the aforementioned Green's function:

[0097] .

[0098] in, This represents the integral variable.

[0099] S243, Substitute the Green's function from step S241 into the convolution formula of step S242, and obtain the analytical solution of the advection diffusion-heat loss control equation from step S11 through variable substitution and integration:

[0100] .

[0101] S244, Solve for the custom energy retention rate using the analytical solution of the advection diffusion-heat loss control equation from step S11. The energy retention rate of the theoretical model is obtained, denoted as . .

[0102] in, This represents the error function.

[0103] S25, for several groups D and R The combination of will have different D The combinations are divided into different classes, and the fitting method is used to find... The value makes and The best match is used to determine the correction factor. τ is the coefficient obtained from the fitting. k D To take into account the inner diameter of the pipe D The correction factor is used to finally obtain:

[0104] .

[0105] S30, temperature diffusivity of turbulent flow α Perform precise calibration / correction (this step assumes the ambient temperature is constant and uniform), including the following sub-steps S21~S25:

[0106] S31, In turbulent flow, physical quantities can be decomposed into average and fluctuating components, thus... V Decomposed into ,Will Decomposed into .

[0107] in, and It is the pulsation quantity in turbulent flow. and It is the Reynolds average quantity in turbulent flow. , .

[0108] S32, Apply Reynolds average to the advection diffusion-heat loss control equation in step S11:

[0109] .

[0110] S33, using the gradient diffusion assumption to assess turbulent heat flux Modeling Substituting this into the advection diffusion-heat loss control equation after applying Reynolds average, we obtain:

[0111] .

[0112] in, express α The theoretical value, α 水 The thermal conductivity of water is represented by... α 水 The order of magnitude is extremely small (typically on the order of 10). -6 Therefore, use express α The theoretical value.

[0113] S34, due to It is not a constant; it depends on the turbulent characteristics of the flow (such as turbulence intensity and scale), so the Prandtl number of turbulence is used. Will With turbulent viscosity Related:

[0114] .

[0115] Where β represents turbulent kinetic energy, oh This represents the specific dissipation rate, since and ,so (The vortex frequency is inversely proportional to the characteristic time scale, which is consistent with the literature) [6] (Consistent with the conclusions in the text) D Indicates the inner diameter of the pipe. It is a dimensionless parameter (typically ranging from 0.7 to 2.0, depending on the fluid and flow conditions).

[0116] S35, the flow inside the pipe is turbulent, and cavitation effects are not considered. SST turbulent kinetic energy-specific dissipation rate turbulence model.

[37] This model exhibits high computational accuracy for attached boundary layer turbulence and moderately separated flows, making it suitable for fluid calculations in this embodiment. It combines the high accuracy of the standard turbulent kinetic energy-specific dissipation rate turbulence model in the near-wall region with its robustness in the mainstream region.

[38] It is particularly suitable for scenarios with adverse pressure gradients and moderately separated flows (such as the boundary layer on the pipe wall), and can automatically switch between near-wall and far-field turbulence calculations through a mixing function.

[0117] CFD simulation can obtain the complete flow field, but the resulting turbulent temperature diffusion coefficient is a radial distribution function. The governing equations need to describe the one-dimensional axial diffusion effective coefficient. Given that axial temperature diffusion in the pipe is mainly dominated by the high-velocity mainstream region, this step proposes a velocity-weighted averaging method (convection-dominated averaging) for one-dimensional simplification, specifically in this embodiment:

[0118] The SST turbulent kinetic energy-specific dissipation rate turbulence model was adopted.

[37] right D and V After performing computational fluid dynamics simulation on the combined pipe to calculate its internal complete flow field, the internal complete flow field is simplified to one dimension using the velocity-weighted averaging method to obtain:

[0119] .

[0120] in, The cross-sectional average velocity, r Represents radial coordinates, This indicates the velocity distribution (in fully developed pipe turbulence, the high-speed region, i.e. the pipe center, contributes more to heat transport and is given a higher weight). Represents radial coordinates r The turbulent temperature diffusion coefficient at that location.

[0121] S36, Select several different groups D and V By combining the elements and repeating step S35, several different combinations are obtained. Data points, and their relationship with D × V Perform linear fitting to determine the correction factor k α ,but:

[0122] .

[0123] S40, Construct the frequency domain expression for the return water temperature of the branched pipe network. In this embodiment, the branched pipe network adopts a quality regulation mode (constant flow rate, regulating supply water temperature). Under this mode, it is assumed that the hydraulic conditions are stable. This step specifically includes the following sub-steps S41~S45:

[0124] S41, decompose the branched pipe network according to the heat user load nodes. n Each independent branch contains a complete flow path from the heat source to the corresponding load node and back to the heat source.

[0125] S42 divides the composition of independent branches into upstream segments. Load nodes and the downlink segment symmetrical to the uplink segment. , among which, in the first k Temperature drop at the load node of each independent branch , k =1,2,3,……, n .

[0126] Each independent branch has 2 m It consists of several pipes, including the upstream section. m The pipeline and the downstream section m A pipeline.

[0127] The single-pipe frequency domain transfer function of each pipe is:

[0128] .

[0129] Transfer function for the uplink segment:

[0130] .

[0131] Transfer function for the downlink segment:

[0132] .

[0133] in, i This indicates the pipe number of a separate, independent branch. i =1,2,3,……, m , L i , G i , α i , c i 、V i These are the pipe length of the independent branch, the single-pipe frequency domain transfer function, and the corrected... α Corrected c as well as V .

[0134] S43, the total return water temperature is recorded as Then it flows into the first k The water supply from each independent branch line flows through the upstream section. Load nodes and downlink The next k The return water temperature of each independent branch is:

[0135] .

[0136] in, This indicates the change in the temperature of the heat source.

[0137] S44, calculate the flow weight of each independent branch according to the flow ratio. Later Result obtained by flow weighting:

[0138] .

[0139] in, Q k Indicates the first k The flow rate of an independent branch, Q i Indicates the first i The flow rate of the pipe.

[0140] S45, Decoupling and rearranging the formula in step S44 yields the frequency domain expression for the return water temperature of the branched pipe network:

[0141] .

[0142] Among them, the heat source transfer function , No. k Load transfer function of each load node .

[0143] S50 uses the inverse fast Fourier transform algorithm to perform an inverse Laplace transform on the frequency domain expression of the return water temperature of the branched pipe network, thereby obtaining the return water temperature response in the time domain to realize the thermal dynamic modeling of the branched pipe network.

[0144] Test case

[0145] This test case uses the transfer function-based thermal dynamic modeling method for branched pipe networks from the embodiments for actual testing.

[0146] (1) Correction of temperature loss coefficient

[0147] In step S20, for c The experimental parameters were adjusted and set as shown in Table 1 below: The insulation material was polyurethane foam, with different pipe diameters... DCorresponding to different insulation layer thicknesses and thermal resistance per unit length R Typical thermal resistance of DN900 pipe R The range is 0.4 m / W to 1.0 K·m / W, 0.5 K·m / W to 1.2 K·m / W for DN500, and 1.2 K·m / W to 2.3 K·m / W for DN250.

[0148] Table 1 (Experimental Parameter Settings)

[0149]

[0150] Verification of theoretical values ​​in step S216 To ensure universality, multiple sets of simulation experiments with different thermal resistances were conducted on the same pipe diameter, and the baseline model was solved using FDM.

[36] The inlet boundary condition is a rectangular pulse ( T 输入 =10℃, Q =3600s). The change of pipe outlet temperature over time was obtained through numerical calculation. Numerical dissipation does not affect the total energy, so the energy retention rate can be calculated by the ratio of the temperature-time integral area to the initial pulse area. , and then establish D , R and c The relationship. Based on energy retention rate. To establish a bridge, a fitting regression method is used. c and D , R The relationship. The fitting results are as follows: Figure 2 As shown.

[0151] Figure 2 This is the relationship between the temperature loss coefficient and the thermal resistance per unit length for different pipe diameters in the test examples of this invention.

[0152] like Figure 2 As shown, under commonly used insulation materials, the heat loss coefficient of the advection diffusion-heat loss control equation for a branched pipe network is... c With pipe diameter D and thermal resistance R Significantly correlated. Within the typical thermal resistance range, larger diameter pipes exhibit lower [correlation]. c This does not mean that larger diameter pipes have less heat loss; on the contrary, larger diameter pipes usually have a larger total heat loss, but their corresponding mass flow rate is larger, so more heat loss is required per unit temperature drop, resulting in a relatively lower temperature loss coefficient.

[0153] During the fitting process, it was found that classifying different pipe diameters into different groups for fitting yielded better results. The fitting correction results indicate that, considering the thermal inertia of the pipe wall, the actual heat loss is higher than the theoretical calculation (correction factor). kγ >1), and literature

[22] The conclusions are consistent: during heating, the thermal inertia of the pipe wall reduces the rate of temperature rise, and during cooling, it slows down the rate of temperature drop. Optimal correction factors for different pipe diameters. k γ Differences exist: the larger the pipe diameter, the larger the correction factor. k γ The smaller the value, the weaker the influence of pipe wall thermal inertia on water temperature propagation, which is consistent with the literature.

[33]

[34] The results obtained are consistent.

[0154] Finally, the advection diffusion-heat loss control equation derived from this test example, considering the pipe's thermal inertia and diameter variation, t =1.673, k D =-0.6475, the expression for the temperature loss coefficient is: .

[0155] (2) Correction of turbulent temperature diffusivity

[0156] In step S30, in the fully developed turbulent flow of the circular pipe, the fully developed section is extracted ( x / D =80) speed With turbulent diffusion coefficient Radial distribution. Specific CFD simulation references.

[37] The final CFD simulation results are as follows: Figure 3 and Figure 4 As shown.

[0157] Figure 3 This is the dimensionless flow velocity distribution curve as a function of dimensionless radius in the CFD simulation results of the test example of this invention.

[0158] like Figure 3 As shown, the velocity distribution calculated by CFD is as expected, exhibiting a parabolic shape, consistent with the literature.

[39] The results were consistent. As the pipe diameter increased, the difference in turbulent velocities increased, with the flow velocity at the center of the DN900 pipe approaching 1.2 times the average flow velocity.

[0159] Figure 4 This is the distribution curve of thermal diffusivity as a function of dimensionless radius in the CFD simulation results of the test example of this invention.

[0160] like Figure 4 As shown, The distribution is more uneven, spanning approximately 10 kilometers from the pipe wall to the center. 3 The magnitude is due to the weak turbulent diffusion caused by the low flow velocity in the near-wall region.

[0161] After processing the obtained distribution data using the velocity averaging method, the relationship between the temperature diffusion coefficient and the pipe diameter and velocity (which corresponds to the cross-sectional average velocity in CFD simulation) can be obtained.

[0162] Figure 5 This is the fitting result of the turbulent temperature diffusion coefficient in the test example of the present invention.

[0163] like Figure 5 As shown, the overall one-dimensional temperature diffusivity of the pipeline increases with both velocity and pipe diameter, which is basically consistent with the theoretically derived expression and the literature. [6] The conclusions drawn are consistent with this trend. The fitting results are shown in Table 2:

[0164] Table 2 (Fitting results for thermal diffusivity)

[0165]

[0166] final k α =0.854, the fitting expression for pipe diameter and velocity is: .

[0167] (3) Verification of single-tube frequency domain transfer function

[0168] Select total length L =1 km, inner diameter D The test was conducted using a straight pipe with a diameter of 500 mm. The water flow velocity was assumed to be... V =1 m / s, initial steady-state inlet temperature 70℃. At the 10th minute, the inlet temperature increases by 2℃. The experimentally measured pipe outlet temperature is as follows: Figure 6 Midpoint shown (1 minute interval). Calculate based on pipe diameter and flow velocity. c =1.62×10 -6 (s) -1 ), α =0.427 m 2 / s. Analytical solution of the transfer function (time-domain response) as follows: Figure 6 As shown by the solid red line in the middle.

[0169] The results show that the transfer function solution accurately predicts the temperature transfer delay time. An offset occurs in the initial stage of temperature rise due to experimental measurement fluctuations, but the offset is relatively small. After the temperature stabilizes, the temperature reduction (heat loss effect) predicted by the transfer function solution is basically consistent with the experimental results, verifying the accuracy of the prediction. c The accuracy of the correction is confirmed. The temperature diffusion slope (front smoothness) presented by the transfer function solution is consistent with the experimental results, verifying the accuracy of the correction. α The accuracy of the correction.

[0170] Figure 6The solid blue line represents data based on literature. [6] The analytical solution calculated shows that the predicted temperature loss is significantly overestimated, resulting in a temperature drop of nearly 0.5℃ during the approximately 20-minute transfer process. This is because the literature... [6] The RKOK turbulence model used is applicable to the high Reynolds number region, where convective heat transfer is enhanced, leading to an overestimation of the γ value. However, the flow velocity in typical heating pipe networks rarely reaches the high Reynolds number region; therefore, the literature... [6] The recommended expression given is not accurate in describing this experiment. By comparing with experimental and literature results, the accuracy of the corrected parameter expression in step S42 of this embodiment was verified.

[0171] Thus far, a heat transfer function model for a single heating pipe, combining physical accuracy and computational efficiency, has been successfully established and validated. However, actual district heating systems are networks composed of numerous pipes connected by complex topologies (branched, looped). The dynamic behavior of a single pipe in the network is coupled and superimposed due to flow branching, merging, and flow distribution. Therefore, extending the validated single-pipe transfer function model to the entire heating network and constructing a unified networked dynamic model is a crucial next step in modeling heating networks to support the collaborative optimization of integrated energy systems.

[0172] (4) Comparison of pipeline network modeling methods

[0173] Two methods were used to construct transfer functions for the example network, and both methods employed IFFT to transform the frequency domain to the time domain. The final results are as follows: Figure 7 As shown, the transfer function obtained by the global matrix method for pipe networks may contain high-order poles due to the presence of high-order rational fractions in the denominator, leading to pole sensitivity. Using the IFFT method may result in numerical overflow. To ensure stability, only the Stehfest method can be used for time-frequency conversion of the results. In contrast, the transfer function obtained by the branch-superposition method for pipe networks is in pure exponential form and can therefore be considered unconditionally stable, allowing for time-frequency conversion using the IFFT method. Ultimately, it was found that the transfer function obtained by the branch method accurately reflects the different delay times and different contributions to return water temperature corresponding to the two branches, while the result of the global matrix transfer function of the pipe network analyzed by the Stehfest method only roughly reflects the overall return water temperature of the pipe network, lacking accuracy. However, this does not mean that the solution obtained by the global matrix transfer function method for pipe networks is inaccurate; it simply means that during the time-frequency conversion process, the Stehfest method leads to variations in the result when dealing with non-smooth, rapidly changing functions with steep gradients. (See the literature...)

[40] The details are explained in the text. The branching method decomposes the complex transfer function into an IFFT-friendly exponential form through topological decoupling, avoiding the numerical ill-conditioning problem of the global matrix transfer function method for pipeline networks. However, its cost is that it ignores the thermal coupling between branches, and may fail in scenarios with sudden changes in traffic or strong interactions.

[0174] [1]QIN X, SUN H, SHEN X, et al. A generalized quasi-dynamic model forelectric-heat coupling integrated energy system with distributed energyresources[J]. Applied Energy, 2019, 251: 113270.DOI: 10.1016 / j.apenergy.2019.05.073。

[0175] [2]WANG Y, DONG H, MA K, et al. Multi frequency stabilityoptimization of integrated energy systems considering virtual energy storagecharacteristics ofheating networks[J]. Applied Thermal Engineering, 2024,257: 124254.DOI:10.1016 / j.applthermaleng.2024.124254。

[0176] [3]XU X, LYU Q, QADRAN M, et al. Quantification of Flexibility of aDistrict Heating System for the Power Grid[J]. IEEE Transactions onSustainable Energy, 2020, 11(4): 2617-2630.DOI: 10.1109 / TSTE.2020.2968507。

[0177] [4]CHEN F, DENG H, CHEN Y, et al. Distributed robust cooperativescheduling of multi-region integrated energy system considering dynamiccharacteristics ofnetworks[J]. International Journal of Electrical Power&Energy Systems, 2023, 145: 108605.DOI: 10.1016 / j.ijepes.2022.108605。

[0178] [5]ZIVKOVIC B, PRICA S, MASLOVARIC B, et al. Prediction of thermaltransients in district heating systems[J]. Energy Conversion and Management,2009, 50(9): 2167-2173.DOI: 10.1016 / j.enconman.2009.04.034。

[0179] [6]CHERTKOV M, NOVITSKY N. Thermal Transients in District HeatingSystems[J].Energy, 2019, 184: 22-33.DOI: 10.1016 / j.energy.2018.01.049。

[0180] [7]VAN DER HEIJDE B, FUCHS M, RIBAS C, et al. Dynamic equation-basedthermo-hydraulic pipe model for district heating and cooling systems[J].Energy Conversion and Management, 2017, 151: 158-169.DOI:10.1016 / j.enconman.2017.08.072。

[0181] [8]SALAVATI S, REES S J. Dynamic thermal response modelling ofturbulent fluid flow through pipelines with heat losses[J]. InternationalJournal of Heat and Mass Transfer, 2020, 151: 119440.DOI: 10.1016 / j.ijheatmasstransfer.2020.119440。

[0182] [9]LI X Z, WU J Z, JENKINS N, et al. Combined analysis of electricityand heat networks[J]. Applied Energy, 2016, 162: 1238-1250.DOI: 10.1016 / j.apenergy.2015.01.102。

[0183]

[10] PENG Z G, GUO Q L, SUN H B. Interactions of district electricityand heating systems considering time-scale characteristics based on quasi-steadymulti-energy flow[J]. Applied Energy, 2016, 167: 230-243.DOI:10.1016 / j.apenergy.2015.10.099。

[0184]

[11] SARTOR K, QUOILIN S, DEWALLEF P. Simulation and optimization of aCHP biomass plant and district heating network[J]. Applied Energy, 2014, 130:474-483.DOI: 10.1016 / j.apenergy.2014.02.039。

[0185]

[12] FANG T T, LAHDELMA R. Genetic optimization of multi-plant heatproduction in district heating networks[J]. Applied Energy, 2015, 159: 610-619.DOI: 10.1016 / j.apenergy.2015.08.125。

[0186]

[13] KUNZ V J, HALDI P A, SARLOS G. Dynamic behaviour of districtheatingnetworks[J]. Fernwaerme International, 1991, 20: 104-119。

[0187]

[14] PARK Y S, KIM W T, KIM B K, et al. Simple models for operationaloptimization, chapter 2: state of the art[M] / / BØHM B. IEA District Heatingand Cooling, Annex VI, Report 2002: S1. Kgs. Lyngby: Technical University ofDenmark, 2002:7-10。

[0188]

[15] DE-CARVALHO R, VALENTE R A F, ANDRADE-CAMPOS A. Dynamictemperaturesimulation in district heating systems in Denmark regardingpronounced transient behaviour[J]. Journal of Civil Engineering andManagement, 2011, 17(1): 79-88。

[0189]

[16] GABRIELAITIENE I. Numerical simulation of a district heatingsystem with emphases on transient temperature behaviour[D]. Vilnius: VilniusGediminas Technical University, 2011。

[0190]

[17] LI Z, WU W, SHAHIDEHPOUR M, et al. Combined Heat and PowerDispatch Considering Pipeline Energy Storage of District Heating Network[J].IEEE Transactions on Sustainable Energy, 2016, 7(1): 12-22.DOI: 10.1109 / TSTE.2015.2467383。

[0191]

[18] STEVANOVIC V D, ZIVKOVIC B, PRICA S, et al. Prediction of thermaltransients in district heating systems[J]. Energy Conversion and Management,2009, 50(9): 2167-2176.DOI: 10.1016 / j.enconman.2009.04.034。

[0192]

[19] GABRIELAITIENE I, KACIANAUSKAS R, SUNDEN B. Thermo-hydraulicfinite element modelling of district heating network by the uncoupledapproach[J]. Journal of Civil Engineering and Management, 2003, 9(2): 153-162.DOI:10.1080 / 13923730.2003.10531321。

[0193]

[20] WANG H, MENG H. Improved thermal transient modeling with new 3-order numerical solution for a district heating network with consideration ofthe pipe wallsthermal inertia[J]. Energy, 2018, 160: 171-183.DOI: 10.1016 / j.energy.2018.06.214。

[0194]

[21] BETANCOURT SCHWARZ M, MABROUK M T, SANTO SILVA C, et al. Modifiedfinite volumes method for the simulation of dynamic district heating networks[J]. Energy, 2019, 182: 954-964.DOI: 10.1016 / j.energy.2019.06.038。

[0195]

[22] ARABKOOHSAR A, KHOSRAVI M, ALSAGRI A S. CFD analysis of triple-pipesfor a district heating system with two simultaneous supply temperatures[J].International Journal of Heat and Mass Transfer, 2019, 141: 432-443.DOI:10.1016 / j.ijheatmasstransfer.2019.06.101。

[0196]

[23] CAPONE M, GUELPA E, VERDA V. Numerical model for the analysis ofthermal transients in district heating networks[C] / / E3S Web of Conferences.2020, 197: 01004.DOI: 10.1051 / e3sconf / 202019701004。

[0197]

[24] VIVIAN J, DE URIBARRI P M Á, EICKER U, et al. The effect ofdiscretization on the accuracy of two district heating network models basedon finite-differencemethods[J]. Energy Procedia, 2018, 149: 625-634.DOI:10.1016 / j.egypro.2018.08.227。

[0198]

[25] WANG Y, YOU S, ZHANG H, et al. Thermal transient prediction ofdistrict heating pipeline: optimal selection of the time and spatial stepsfor fast and accuratecalculation[J]. Applied Energy, 2017, 206: 900-910.DOI:10.1016 / j.apenergy.2017.08.061。

[0199]

[26] OPPELT T, URBANECK T, GROSS U, et al. Dynamic thermo-hydraulicmodel of district cooling networks[J]. Applied Thermal Engineering, 2016,102: 336-345.DOI: 10.1016 / j.applthermaleng.2016.03.168。

[0200]

[27] SOLAR ENERGY LABORATORY. TRNSYS Manual 17 - Volume 4 -Mathematical Reference[R]. Madison: University of Wisconsin-Madison, 2009。

[0201]

[28] STEVANOVIC V, PRICA S, MASLOVARIC B, et al. Efficient numericalmethod for district heating system hydraulics[J]. Energy Conversion andManagement, 2007, 48(5): 1536-1543.DOI: 10.1016 / j.enconman.2006.11.020。

[0202]

[29] WANG Y, LI Y, ZHANG Y, et al. Optimized operation of integratedenergy systems accounting for synergistic electricity and heat demandresponse under heat loadflexibility[J]. Applied Thermal Engineering, 2024,243: 122640.DOI: 10.1016 / j.applthermaleng.2024.122640。

[0203]

[30] LIU X, MANCARELLA P, WU J. A tool for integrated analysis ofmulti-vector district energy networks[C] / / 2015 IEEE Eindhoven PowerTech.Eindhoven: IEEE, 2015: 1-6.DOI: 10.1109 / PTC.2015.7232647。

[0204]

[31] GOU X, CHEN Q, HE K. Real-time quantification for dynamic heatstorage characteristic of district heating system and its application indispatch ofintegrated energy system[J]. Energy, 2022, 259: 124960.DOI:10.1016 / j.energy.2022.124960。

[0205]

[32] DÉNARIÉ A, APRILE M, MOTTA M. Heat transmission over long pipes:New model for fast and accurate district heating simulations[J]. Energy,2019, 166: 267-276.DOI: 10.1016 / j.energy.2018.09.186。

[0206]

[33] SARTOR K, DEWALLEF P. Experimental validation of heat transportmodelling in district heating networks[J]. Energy, 2017, 137: 961-968.DOI:10.1016 / j.energy.2017.02.161。

[0207]

[34] ZHAO T, LI H, LI X, et al. A frequency domain dynamic simulationmethod for heat exchangers and thermal systems[J]. Energy, 2024, 286:129552.DOI: 10.1016 / j.energy.2023.129552。

[0208]

[35] LIN X, ZHANG N, LUO Z, et al. Balanced operation strategies ofdistrict heating systems based on dynamic hydraulic-thermal modeling[J].Energy and Built Environment, 2025, 6(3): 466-483.DOI: 10.1016 / j.enbenv.2024.01.005。

[0209]

[36] WANG H, MENG H, ZHU T. New model for onsite heat loss stateestimation of general district heating network with hourly measurements[J].Energy Conversion and Management, 2018, 157: 71-85.DOI: 10.1016 / j.enconman.2017.11.062。

[0210]

[37] MENTER F R. Two-equation eddy-viscosity turbulence models forengineering applications[J]. AIAA Journal, 1994, 32(8): 1598-1605.DOI:10.2514 / 3.12149。

[0211]

[38] BELHOCINE A, IBRAHEEM O F. Numerical simulation of thermallydeveloping turbulent flow through a cylindrical tube[J]. The InternationalJournal of Advanced Manufacturing Technology, 2019, 102(5-8): 2001-2012.DOI:10.1007 / s00170-019-03315-y。

[0212]

[39] SARMA P K, CHADA K, DHARMA RAO V, et al. Evaluation of momentumand thermal eddy diffusivities for turbulent flow in tubes[J]. InternationalJournal of Heat and Mass Transfer, 2010, 53(5-6): 1237-1242.DOI:10.1016 / j.ijheatmasstransfer.2009.11.023。

[0213]

[40] KUO L. Review of inverse Laplace transform algorithms for Laplace-space numerical approaches[J]. Numerical Algorithms, 2012, 63(2):339-355.DOI: 10.1007 / s11075-012-9625-3.

[0214] Those skilled in the art should understand that this invention is not limited to the above embodiments. The embodiments and descriptions in the specification are merely illustrative of the principles of the invention. Various changes and modifications can be made to this invention without departing from its spirit and scope, and all such changes and modifications fall within the scope of the invention as claimed. The scope of protection of this invention is defined by the appended claims and their equivalents.

Claims

1. A method for thermal dynamic modeling of a branched pipe network based on transfer function, wherein the branched pipe network includes multiple independent branches, and each independent branch includes multiple pipes. Its features are, Includes the following steps: S10, Establish the advection diffusion-heat loss control equation for a single pipe. Later on time t Perform a Laplace transform and derive the single-tube frequency domain transfer function: ， in, Indicates the axial position of the pipe x Location, Time t The difference between the fluid temperature and the ambient temperature at that time. V Indicates the average flow velocity of the fluid. Represents the convection term. This represents the heat loss term. γ Indicates the temperature loss coefficient. Represents the diffusion term. α Indicates the turbulent temperature diffusivity. Custom , Represents the Laplace transform. s Represents a complex frequency variable; S20, considering the thermal inertia of the pipe wall and the inner diameter, the temperature loss coefficient is determined through a finite difference numerical simulation experiment. γ Perform precise calibration and correction; S30, after applying Reynolds averaging to the advection diffusion-heat loss control equation, the turbulent temperature diffusion coefficient is determined by combining the gradient diffusion assumption and the Prandtl number of turbulence. α The associated physical quantities Subsequently, computational fluid dynamics simulations of the pipeline were performed using the SST turbulent kinetic energy-specific dissipation rate turbulence model, and the internal flow field information related to the simulation was compressed into one-dimensional information using the velocity-weighted averaging method. Finally, by combining the one-dimensional information with the associated physical quantity, the accurate turbulent temperature diffusion coefficient is obtained. α ; S40, the branched pipe network is decomposed into multiple independent branches according to the heat user load nodes, and the independent branches are decomposed into multiple pipes, based on the calibrated and corrected temperature loss coefficient. γ and the turbulent temperature diffusivity α And based on the flow weighting principle, a frequency domain expression for the return water temperature of the branched pipe network is constructed; S50, perform an inverse Laplace transform on the frequency domain expression of the return water temperature of the branched pipe network to obtain the return water temperature response in the time domain, thereby realizing the thermal dynamic modeling of the branched pipe network.

2. The method for thermal dynamic modeling of branched pipe networks based on transfer function according to claim 1, characterized in that: in, Step S10 includes the following sub-steps: S11, Establish the advection diffusion-heat loss control equation for a single pipe: ； S12, Custom Then, the Laplace transform of the advection diffusion-heat loss control equation is used to obtain the ordinary differential equation in the frequency domain: ； S13, Solve the ordinary differential equation to obtain the single-tube frequency domain transfer function: 。 3. The method for thermal dynamic modeling of branched pipe networks based on transfer function according to claim 2, characterized in that: in, In step S12, it is assumed that the solution to the ordinary differential equation is in the form of: , λ It is a constant. The characteristic equation is: , The final root is: .

4. The method for thermal dynamic modeling of branched pipe networks based on transfer function according to claim 1, characterized in that: in, Step S20 includes the following sub-steps: S21, based on heat transfer theory, the heat loss term is associated with... Theoretical value and introduce the inner diameter of the pipe. D , and thus γ Theoretical value: ， In the above formula, ρ Indicates fluid density, c p Indicates the specific heat capacity of a fluid. R This represents the thermal resistance per unit length of the pipe. and D and R Related; S22, define a custom energy retention rate: ； S23. After establishing a finite difference model of pipe heat loss considering the pipe wall thermal inertia as a baseline model, a temperature pulse is applied as input at the pipe inlet of the baseline model. The temperature response at the pipe outlet is calculated by running the baseline model. The energy retention rate of the baseline model is obtained by spatial integration of the temperature response based on the user-defined energy retention rate, denoted as the numerical simulation solution. ; S24, using the analytical solution of the aforementioned advection diffusion-heat loss control equation. Solving for the custom energy retention rate yields the energy retention rate of the theoretical model, denoted as . , in, Indicates the pulse half-width. This indicates the initial height of the temperature pulse. Represents the error function; S25, for several groups D and R The combination of will have different D The combinations are divided into different classes, and the fitting method is used to find... The value makes and The best match is used to determine the correction factor. τ is the coefficient obtained from the fitting. k D To take into account the inner diameter of the pipe D The correction factor The final result is: .

5. The method for thermal dynamic modeling of branched pipe networks based on transfer function according to claim 4, characterized in that: in, Step S21 includes the following sub-steps: S211, The heat loss per unit length of the pipe is given by Newton's law of cooling: ； S212, based on the heat capacity per unit length of the fluid C = ρAc p The relationship between the rate of decrease in fluid temperature due to heat loss per unit length is determined as follows: ， in, A This represents the cross-sectional area of ​​the fluid inside the pipe; S213, derived from heat transfer theory, is as follows: ； S214, the heat loss term This represents the contribution of temperature decay, therefore its theoretical value is expressed as: ； S215, after sorting: ； S216, Introducing the inner diameter of the pipe D Substitution A From: 。 6. The method for thermal dynamic modeling of branched pipe networks based on transfer function according to claim 4, characterized in that: in, In step S23, in the reference model x Apply a persistent [condition] at =0. Q The temperature in seconds is T 输入 The rectangular temperature pulse is used as input. Assumption ,but: , The boundary value problem in the above equation is then transformed into an initial value problem, that is, the action of the heat source in time can be regarded as the generation of a rectangular temperature pulse in space: ； In step S24, the convection-diffusion-attenuation equation within the pipe is given using the Green's function method. Green's function in the form of point source solutions: ， The solution to the rectangular temperature pulse is the convolution of the Green's function: ， in, Represents the integral variable. Substituting the Green's function into the convolution formula, and obtaining the analytical solution of the advection diffusion-heat loss control equation through variable substitution and integration: 。 7. The method for thermal dynamic modeling of branched pipe networks based on transfer function according to claim 1, characterized in that: in, Step S30 includes the following sub-steps: S31, will V Decomposed into ,Will Decomposed into , in, and It is the pulsation quantity in turbulent flow. and It is the Reynolds average quantity in turbulent flow. , ; S32, Applying Reynolds average to the advection diffusion-heat loss control equation: ； S33, using the gradient diffusion assumption to assess turbulent heat flux Modeling Substituting this into the advection diffusion-heat loss control equation after applying Reynolds average, we obtain: ， in, express α The theoretical value, α 水 The thermal conductivity of water is represented by... α 水 The order of magnitude is extremely small, so use express α Theoretical value; S34, using the Prandtl number for turbulence. Will With turbulent viscosity Related: ， Where β represents turbulent kinetic energy, ω This represents the specific dissipation rate, since and ,so , D This indicates the inner diameter of the pipe; S35, using the SST turbulent kinetic energy-specific dissipation rate turbulence model for... D and V After performing computational fluid dynamics simulation on the combined pipes to calculate its internal complete flow field, the internal complete flow field is simplified in one dimension using the velocity-weighted averaging method to obtain: ， in, A This represents the cross-sectional area of ​​the fluid inside the pipe. The cross-sectional average velocity, r Represents radial coordinates, Represents velocity distribution. Represents radial coordinates r Turbulent temperature diffusivity at that location; S36, Select several different groups D and V By combining the elements and repeating step S35, several different combinations are obtained. Data points, and their relationship with D × V Perform linear fitting to determine the correction factor k α ,but: 。 8. The method for thermal dynamic modeling of branched pipe networks based on transfer function as described in claim 1, Its features are: in, Step S40 includes the following sub-steps: S41, decompose the branched pipe network according to the heat user load nodes. n Each independent branch contains a complete flow path from the heat source to the corresponding load node and back to the heat source; S42, the composition of the independent branch is divided into uplink segments. Load nodes and downlink segments symmetrical to the uplink segments. , among which, in the first k Temperature drop at the load node of each of the independent branches , k =1,2,3,……, n ， Each of the aforementioned independent branches is 2 m It consists of several pipes, including the upstream section. m The pipeline and the downstream section mentioned above m The pipeline mentioned above, The single-tube frequency domain transfer function for each of the aforementioned pipes is: ， ， ， in, i The number of the pipe representing a single, independent branch. i =1,2,3,……, m ， L i , G i , α i , γ i 、V i These are, respectively, the pipe length corresponding to the pipeline of the independent branch, the single-pipe frequency domain transfer function, and the corrected... α Corrected γ as well as V ; S43, the total return water temperature is recorded as Then it flows into the first k The water supply from the aforementioned independent branch flows through the upstream section. The load node and the downlink segment The next k The return water temperature of each of the aforementioned independent branches is: ， in, Indicates the change in the temperature of the heat source; S44, Calculate the flow weight of each independent branch according to the flow ratio. w k Later Result obtained by flow weighting: ； S45, Decoupling and rearranging the formula in step S44 yields the frequency domain expression for the return water temperature of the branched pipe network: ， Among them, the heat source transfer function , No. k Load transfer function of each load node .

9. The method for thermal dynamic modeling of branched pipe networks based on transfer function according to claim 8, characterized in that: in, In step S44, the flow weight of each of the independent branches , Q k Indicates the first k The flow rate of an independent branch, Q i Indicates the first i The flow rate of the pipe.

10. The method for thermal dynamic modeling of branched pipe networks based on transfer function according to claim 1, characterized in that: in, In step S50, the inverse fast Fourier transform algorithm is used to perform the inverse Laplace transform.

Images