Coupling calculation method of heat transfer in chamber and heat transfer in chamber wall based on one-dimensional heat conduction

Abstract

The scheme designs a kind of coupling calculation method of cavern thermodynamics and cavern wall heat conduction based on one-dimensional heat conduction, according to the equality of the convective heat transfer power of cavern air and cavern wall and the surface heat conduction power of cavern wall obtained based on temperature distribution, the sealing layer temperature T m is obtained m , and the sealing layer temperature T m is taken as the boundary condition, the cavern air thermodynamics calculation and the one-dimensional heat conduction calculation of cavern wall are carried out respectively, the cavern air temperature and pressure are obtained through the cavern air thermodynamics calculation, the temperature distribution of cavern wall along the radial direction is obtained through the one-dimensional heat conduction calculation of cavern wall until the end of the set calculation time, the change of surrounding rock temperature in the whole process is obtained, so as to realize the coupling of the thermodynamic process of air in cavern and the heat transfer process of surrounding rock of cavern wall.

Description

Technical Field

[0001] This invention patent belongs to the field of energy storage and relates to the calculation of thermodynamic processes during the filling and releasing of underground artificial gas storage chambers in compressed air energy storage power stations, as well as the heat conduction process of the chamber walls caused by changes in air temperature inside the chamber. Specifically, it is a coupled calculation method of chamber thermodynamics and chamber wall heat conduction based on one-dimensional heat conduction. Background Technology

[0002] Compressed air energy storage technology is a large-scale, long-term energy storage and power generation technology. Its technical principle is as follows: Figure 3 As shown. During periods of wind and solar power curtailment and off-peak electricity consumption, compressed air energy storage power stations use electricity to drive compressors to compress air, which is then sent to underground gas storage facilities for storage. When peak electricity demand arrives, the high-pressure air in the storage facilities is heated by a heat exchanger or combustion chamber and then sent to an expander to expand and perform work, driving a generator to generate electricity, thereby achieving peak shaving and valley filling functions for electrical energy.

[0003] Commercial compressed air energy storage power plants currently primarily utilize salt cavern gas storage. These are formed by injecting fresh water into thick underground salt layers or salt domes using a water-soluble extraction method, dissolving the salt rock layer, and then draining the saturated or near-saturated brine. Salt cavern gas storage offers advantages such as good sealing and high stability. However, salt cavern gas storage is significantly limited by geographical conditions; areas where compressed air energy storage power plants need to be built may lack salt rock, making the construction of salt cavern gas storage impossible.

[0004] Besides salt cavern gas storage facilities, underground artificial caverns have also become an important option for compressed air energy storage power plants in recent years. These are gas storage structures with a certain volume, artificially excavated from hard underground rock. Based on the shape of the gas storage chamber, they can be divided into two types: large tank type and tunnel type (see...). Figures 4-5 ).

[0005] Regardless of the type of underground artificial reservoir, its wall structure, from the inside out, mainly consists of three parts: a sealing layer, a concrete lining, and surrounding rock. The airtightness of this type of gas storage is provided by sealing materials such as steel, while the stability and deformation control of the reservoir are provided by the concrete lining.

[0006] The advantages of artificial storage chambers mainly lie in their strong engineering controllability and flexible construction according to needs. Underground artificial storage chambers are currently a key research and demonstration area for compressed air energy storage projects. Therefore, in terms of storage methods, future compressed air energy storage power plants will primarily use artificial storage chambers, supplemented by other storage methods.

[0007] Compressed air energy storage power stations have been validated through the construction of experimental power stations ranging from 5kW to 1.5MW. Currently, most operational projects have reached the 60MW to 100MW level, and 300MW-level projects are in the planning and design phase. The gas storage capacity of a 300MW-level power station is typically 100,000 m³. 3 Even hundreds of thousands of m 3 These engineering gas storage facilities typically operate under standard cycle conditions with pressure variations ranging from several megapascals, and their basic pressure is generally several megapascals, or even as high as ten megapascals or more.

[0008] During such a high-pressure differential variable-amplitude inflation and deflation process, the temperature of the air inside the chamber changes significantly with the inflation process, rising during inflation and decreasing during deflation. The temperature change inside the chamber is controlled not only by the thermodynamic properties of the air but also by the coupled control of convective heat transfer between the air and the chamber walls, as well as the heat conduction processes of the various media layers within the chamber walls.

[0009] During the engineering design phase, understanding the changes in air pressure and temperature within the tunnel is fundamental to the generator set design. On the other hand, as the heat conduction process occurs within the tunnel walls, changes in the temperature of the sealing layer, concrete lining, and surrounding rock significantly impact their stress, strain, and stability. Therefore, performing coupled calculations of tunnel thermodynamics, tunnel wall convective heat transfer, and tunnel wall heat conduction is an important task.

[0010] Establishing energy conservation equations and mass conservation equations based on the conservation of gas energy and mass within the chamber, and then establishing heat conduction equations for the multi-layered media in the chamber wall based on heat conduction, along with their respective boundary conditions, would be extremely complex.

[0011] Existing calculation methods can only calculate the charging and pumping process under standard cycle conditions, and cannot calculate the charging and pumping process starting from atmospheric pressure. This process plays a controlling role in the design, so it is necessary to find a more applicable solution method. Summary of the Invention

[0012] This patent proposes a coupling method based on one-dimensional heat conduction, which can conveniently couple the thermodynamics of the cavern with the heat conduction process of the cavern wall.

[0013] The technical means employed in this invention are as follows:

[0014] The coupled calculation method of tunnel thermodynamics and tunnel wall heat conduction based on one-dimensional heat conduction includes the following steps:

[0015] S1: Determine the geometric dimensions and operating parameters of the chamber, provide the basic thermodynamic and mechanical calculation parameters for the sealing layer, lining, and surrounding rock, and determine the total calculation time t. n The time domain to be calculated is discretized, the time step length Δt is divided, and initial temperature values ​​are assigned to the air and surrounding rock in the cave.

[0016] S2: Calculate the air inflation and deflation rates, inflation temperature, and air density in the chamber at each time step;

[0017] S3: Based on the convective heat transfer between the tunnel wall and the air and the conservation of heat conduction energy in the tunnel wall, and based on the air temperature and tunnel wall temperature distribution at time step i, calculate the temperature T of the sealing layer at time step i+1. m ;

[0018] S4: Assuming the temperature of the sealing layer remains constant in the (i+1)th time step, calculate the temperature of the air inside the chamber in the (i+1)th time step, and use a one-dimensional heat conduction model to calculate the heat conduction of the chamber wall to obtain the temperature distribution of the chamber wall.

[0019] S5: Repeat steps S3-S4 until the total computation time is greater than the total computation time.

[0020] Preferably, in each time step of the cyclic calculation, the sealing layer temperature T is calculated based on the fact that the convective heat transfer power between the air and the wall of the chamber is equal to the heat conduction power on the surface of the wall obtained based on the temperature distribution. m and with the sealing layer temperature T m As boundary conditions, thermodynamic calculations of the air in the cave and one-dimensional heat conduction calculations of the cave wall are performed separately. The temperature and pressure of the air in the cave are obtained through thermodynamic calculations of the air in the cave, and the radial temperature distribution of the cave wall is obtained through one-dimensional heat conduction calculations of the cave wall.

[0021] As a preferred option

[0022] Discretize the time domain to be calculated and select a time step Δt;

[0023] Let coordinate axis r extend outward from the center of the cave, and let the discrete interval of coordinate axis r be Δr;

[0024] Assuming the temperature T of the sealing layer in the time-step chamber is... m Constant and does not change with the thickness or location of the sealing layer;

[0025] This is used to calculate the gas temperature T inside the chamber and the temperature gradient of the surrounding rock at the end of the first time step. Right now Substituting it into equation (1) yields equation (2):

[0026]

[0027] In the formula:

[0028] k is the thermal conductivity of the concrete lining, W·(m·K). -1 ;

[0029] h c The convective heat transfer coefficient between the air inside the chamber and the chamber wall is expressed in W / (m²). 2·K);

[0030] T u The temperature of the surrounding rock;

[0031]

[0032] In the formula:

[0033] h c (TT m The value represents the convective heat transfer power of the surrounding rock at the end of the previous time interval, in W.

[0034] The thermal conductivity of the surrounding rock is expressed in W.

[0035] The required sealing layer temperature T for the next time interval can be further calculated using equation (2). m .

[0036] As a preferred approach, when calculating the fundamental solution for the gas temperature T inside the chamber, it is assumed that the chamber wall temperature is constant, the air density inside the chamber is constant, the filling or pumping rate is constant, and the air is an ideal gas. Furthermore, the energy conservation equation in equation (3) is transformed into an ordinary differential equation.

[0037]

[0038] In the formula:

[0039] V is the volume of the chamber, in meters. 3 ;

[0040] ρ is the density of the air inside the chamber, in kg / m³. 3 ;

[0041] c v The isochoric specific heat of air is 717 J / (kg·K);

[0042] T is the temperature of the air inside the chamber, in K;

[0043] t represents time, in seconds;

[0044] m in (t) is the inflation rate function;

[0045] m ex (t) is the pumping rate function;

[0046] h i To accommodate the enthalpy of air, J;

[0047] h is the enthalpy of the air inside the chamber, in J;

[0048] Z is the air compression coefficient;

[0049] R is the air constant, 286.7 J / (kg·K);

[0050] u is the internal energy of the air inside the cavern, J;

[0051] The heat exchange rate between the air inside the chamber and the chamber wall is expressed in J / s.

[0052] Combining the mass conservation equation in equation (4) and the convective heat transfer equation in equation (5), we obtain the analytical solution (6) for the gas temperature inside the chamber.

[0053]

[0054]

[0055] In the formula:

[0056] h c The convective heat transfer coefficient between the air inside the chamber and the chamber wall is expressed in W / (m²). 2 ·K);

[0057] A c The surface area of ​​the cavern is m. 2 ;

[0058] r0 is the radius of the cavern, in meters;

[0059]

[0060] In formula (6):

[0061]

[0062]

[0063]

[0064] in:

[0065] ρ av To calculate the average air density (kg / m³) within a time step Δt. 3 ;

[0066] ρ0 is the average air density inside the tunnel during the previous time period, in kg / m³. 3 ;

[0067] T0 is the air temperature inside the chamber at the previous time, in K;

[0068] m i The inflation rate for this stage is kg / s;

[0069] m e The pumping rate for this stage is kg / s;

[0070] V is the volume of the chamber, in meters. 3 ;

[0071] T m The temperature of the sealing layer is K.

[0072] As a preferred approach, when calculating the fundamental solution for one-dimensional heat conduction in the surrounding rock, the radius of the cave is set to r1, the radius of the surrounding rock within the temperature influence range is set to r3, the temperature of the outer boundary of the one-dimensional heat conduction is constant, and an x-axis is defined extending outward from a point on the cave wall. The boundary value problem under the one-dimensional heat conduction condition is as follows:

[0073]

[0074] T = T m x = 0, t > 0 (11)

[0075] T = T R0 x = L, t > 0 (12)

[0076] T = f(x) t = 0 (13)

[0077] In the formula, T R0 R is the initial temperature of the surrounding rock; L is the distance of one-dimensional heat conduction, L = r3 - r1.

[0078] As a preferred approach, the above boundary value problem is decomposed into the following two problems:

[0079] First, there is the problem of steady-state nonhomogeneity:

[0080]

[0081] T w =T m x = 0 (15)

[0082] T w =T R0 x=L (16)

[0083] The solution to this definite problem is:

[0084]

[0085] Second, unsteady homogeneous problems:

[0086]

[0087] T f =0 x=0, t>0 (19)

[0088] T f =0 x=L,t>0 (20)

[0089] T f =f(x)-T w (x) t=0 (21)

[0090] The solution to this definite problem is:

[0091]

[0092] The solutions to the original problem equations (10) to (13) are:

[0093] T(x,t)=T w (x)+T f (x, t) (23)

[0094] Substituting equations (17) and (22) into equation (23) and rearranging, we get:

[0095]

[0096] In the formula,

[0097] As a preferred option, the maximum value of m is 20.

[0098] Compared with the prior art, the beneficial technical effects of the present invention are as follows: the coupling method performs the thermodynamic calculation of the air inside the tunnel and the thermal conduction of the surrounding rock separately based on the temperature of the sealing layer as the boundary condition. Compared with the previous method of establishing the thermal conduction equation of the multi-layer medium of the tunnel wall and the respective boundary conditions for coupling solution based on thermal conduction, this method is simpler and more efficient. Moreover, this scheme abandons the more difficult axisymmetric thermal conduction calculation method and adopts a one-dimensional thermal conduction model, which can ensure the calculation accuracy while greatly simplifying the calculation difficulty. Attached Figure Description

[0099] Figure 1 The calculation flowchart for the gas storage tank, which takes into account the thermodynamic processes inside the tunnel and the heat conduction through the tunnel walls, is shown in this scheme.

[0100] Figure 2 for Figure 1 The calculation flowchart for each loop step.

[0101] Figure 3 This is a schematic diagram of a compressed air energy storage power station.

[0102] Figure 4 The figures show a longitudinal section (a) and a cross section (b) of a tunnel-type underground artificial cavern.

[0103] Figure 5 This is a schematic diagram of a large-tank compressed air energy storage sealed chamber.

[0104] Figure 6 A simplified model diagram of heat transfer in the sealing layer, lining, and surrounding rock.

[0105] Figure 7 This is a schematic diagram of the inflation and deflation rate functions.

[0106] Figure 8 This is a schematic diagram of the time discretization method.

[0107] Figure 9 This is a schematic diagram of a one-dimensional heat conduction model.

[0108] Figure 10 This is a curve showing the temperature change of the chamber over time in one embodiment.

[0109] Figure 11 for Figure 10 The curve showing the change of air pressure inside the chamber over time in the embodiment. Detailed Implementation

[0110] like Figure 1 As shown, a coupled calculation method for the thermodynamics of a tunnel and the thermal conduction of the tunnel wall based on one-dimensional heat conduction is proposed, including the following steps:

[0111] S1: Determine the geometric dimensions and operating parameters of the chamber, provide the basic thermodynamic and mechanical calculation parameters for the sealing layer, lining, and surrounding rock, and determine the total calculation time t. n The time domain to be calculated is discretized, the time step length Δt is divided, and initial temperature values ​​are assigned to the air and surrounding rock in the cave.

[0112] S2: Calculate the air inflation and deflation rates, inflation temperature, and air density in the chamber at each time step;

[0113] S3: Based on the convective heat transfer between the tunnel wall and the air and the conservation of heat conduction energy in the tunnel wall, and based on the air temperature and tunnel wall temperature distribution at time step i, calculate the temperature T of the sealing layer at time step i+1. m ;

[0114] S4: Assuming the temperature of the sealing layer remains constant in the (i+1)th time step, calculate the temperature of the air inside the chamber in the (i+1)th time step, and use a one-dimensional heat conduction model to calculate the heat conduction of the chamber wall to obtain the temperature distribution of the chamber wall.

[0115] S5: Repeat steps S3-S4 until the total computation time is greater than the total computation time.

[0116] Preferably, in each time step of the cyclic calculation, the sealing layer temperature T is calculated based on the fact that the convective heat transfer power between the air and the wall of the chamber is equal to the heat conduction power on the surface of the wall obtained based on the temperature distribution. m and with the sealing layer temperature T mAs boundary conditions, thermodynamic calculations of the air in the cave and one-dimensional heat conduction calculations of the cave wall are performed separately. The temperature and pressure of the air in the cave are obtained through thermodynamic calculations of the air in the cave, and the radial temperature distribution of the cave wall is obtained through one-dimensional heat conduction calculations of the cave wall.

[0117] Specifically, regarding the physical model used in the calculation:

[0118] Since the length of a chamber is generally much greater than its transverse diameter, the heat conduction process of the chamber can be considered as an axisymmetric problem.

[0119] For chambers using steel lining as the sealing layer, the thickness of the steel lining is generally 10mm to 20mm, which is relatively thin. Considering that its thermal conductivity is generally more than 10 times that of the lining and surrounding rock, the heat resistance of the steel lining can be ignored, and it can be assumed that the temperature at any point inside the steel lining is the same.

[0120] Considering that the thermal conductivity and isobaric specific heat of the concrete lining [2.94 W / (m·K), 0.96 kJ / (kg·K)] are close to those of the surrounding rock [2–5 W / (m·K), 0.75–1.14 kJ / (kg·K)], the concrete lining and the surrounding rock are considered as the same heat transfer medium with the same thermophysical parameters. In this scheme, both are collectively referred to as the surrounding rock. Furthermore, considering that the conduction depth of air temperature changes within the tunnel is limited in both the lining and the surrounding rock, the parameters of the lining material have a greater impact on the temperature distribution. Therefore, the parameters of both the concrete lining and the surrounding rock are taken from those of the concrete lining.

[0121] like Figure 6 As shown, the outer radius of the steel lining is denoted as r1, and the radius of the temperature or pressure influence range is denoted as r3. The temperature of the surrounding rock remains constant beyond radius r3. The thermophysical properties of the thick-walled cylinder between r1 and r3 are taken from those of the concrete lining. For the sake of simplicity, Figure 6 The steel lining is not marked in the model.

[0122] And the mathematical model for the calculation:

[0123] The existing methods all employ simultaneous solutions to solve the following energy conservation equations, mass conservation equations, convective heat transfer equations, boundary conditions for heat conduction between the air inside the tunnel and the tunnel wall, heat conduction equations of the surrounding rock, and generalized gas state equations.

[0124] The energy conservation equation is as follows:

[0125]

[0126] In the formula:

[0127] V is the volume of the chamber, in meters. 3 ;

[0128] ρ is the density of the air inside the chamber, in kg / m³. 3 ;

[0129] c v The isochoric specific heat of air is 717 J / (kg·K);

[0130] T is the temperature of the air inside the chamber, in K;

[0131] t represents time, in seconds;

[0132] m in (t) is the inflation rate function, see Figure 7 As shown;

[0133] m ex (t) is the pumping rate function, see Figure 7 As shown;

[0134] h i To accommodate the enthalpy of air, J;

[0135] h is the enthalpy of the air inside the chamber, in J;

[0136] Z is the air compression coefficient;

[0137] R is the air constant, 286.7 J / (kg·K);

[0138] u is the internal energy of the air inside the cavern, J;

[0139] The heat exchange rate between the air inside the chamber and the chamber wall is expressed in J / s.

[0140] In equation (3), the inflation rate is positive and the deflation rate is negative. For specific definitions, see [link to equation (3)]. Figure 7 . Figure 7 It also shows that the rate is constant during both the inflation and deflation phases.

[0141] There is also the mass conservation equation:

[0142]

[0143] Convection heat transfer equation:

[0144]

[0145] In the formula:

[0146] h c The convective heat transfer coefficient between the air inside the chamber and the chamber wall is expressed in W / (m²). 2 ·K);

[0147] A c The surface area of ​​the cavern is m. 2 ;

[0148] r0 is the radius of the cavern, in meters.

[0149] Boundary conditions for heat conduction between the air inside the tunnel and the tunnel wall:

[0150]

[0151] T u =T Rw r→∞ (25)

[0152] In the formula:

[0153] k is the thermal conductivity of the concrete lining, W·(m·K). -1 .

[0154] The heat conduction equation of the surrounding rock:

[0155]

[0156] In the formula:

[0157] k u is the thermal conductivity of the medium in the tunnel wall, W / (m·K);

[0158] ρ u The density of the medium in the tunnel wall is kg / m³. 3 ;

[0159] c pu The isobaric specific heat of the medium in the tunnel wall is expressed in J / (kg·K).

[0160] Generalized gas law:

[0161] p=ZρRT (27)

[0162] In the formula:

[0163] p is the gas pressure inside the chamber, in kPa;

[0164] Z is the air compressibility coefficient, which is dimensionless;

[0165] R is the air constant, 286.7 J / (kg·K).

[0166] In the previous solutions, due to the complexity of the energy conservation equation, the solution process generally assumed that the air density remained constant over a period of time. When the working pressure of the chamber varies greatly, this will lead to significant calculation errors. Existing solutions cannot calculate the changes in air within the chamber during the initial inflation process, and the temperature variation during the initial inflation is large due to the large changes in air pressure and density, which is a controlling factor in engineering design.

[0167] Therefore, to address the shortcomings of the existing calculation methods, this solution adopts the following approach:

[0168] Discretize the time domain to be calculated, and select a time step Δt, such as... Figure 8As shown, within each time period, it is assumed that the air density inside the chamber remains constant, and the inner surface temperature of the sealing layer remains constant; it is assumed that the temperature T of the chamber sealing layer within each time step is constant. m Constant and does not change with the thickness or location of the sealing layer;

[0169] In actual calculations, coordinate axes r are set outward from the center of the cave, and the coordinate axes r are discretized, with the discretization interval of the coordinate axes r set as Δr.

[0170] Using these as boundary conditions, calculate the gas temperature T inside the chamber and the surrounding rock temperature gradient at the end of a time step. Right now Substituting this into the boundary condition equation for heat conduction between the air inside the tunnel and the tunnel wall, i.e., equation (1), we obtain the following equation (2):

[0171]

[0172] In the formula:

[0173] k is the thermal conductivity of the concrete lining, W·(m·K). -1 ;

[0174] h c The convective heat transfer coefficient between the air inside the chamber and the chamber wall is expressed in W / (m²). 2 ·K);

[0175] T u The temperature of the surrounding rock;

[0176]

[0177] In the formula:

[0178] h c (TT m The value represents the convective heat transfer power of the surrounding rock at the end of the previous time interval, in W.

[0179] The thermal conductivity of the surrounding rock is expressed in W.

[0180] The required sealing layer temperature T for the next time interval can be further calculated using equation (2). m Next, the temperature T of the sealing layer will be used as the basis for further analysis. m Perform thermodynamic calculations of the air in the cave and thermal conduction calculations of the cave wall for the boundary conditions.

[0181] In calculating the fundamental solution of the gas temperature T inside the chamber, it is assumed that the temperature of the chamber wall is constant, the air density inside the chamber is constant, the filling rate or pumping rate is constant, and the air is an ideal gas (Z=1). The energy conservation equation in equation (3) is further transformed into an ordinary differential equation.

[0182]

[0183] In the formula:

[0184] V is the volume of the chamber, in meters. 3 ;

[0185] ρ is the density of the air inside the chamber, in kg / m³. 3 ;

[0186] c v The isochoric specific heat of air is 717 J / (kg·K);

[0187] T is the temperature of the air inside the chamber, in K;

[0188] t represents time, in seconds;

[0189] m in (t) is the inflation rate function;

[0190] m ex (t) is the pumping rate function;

[0191] h i To accommodate the enthalpy of air, J;

[0192] h is the enthalpy of the air inside the chamber, in J;

[0193] Z is the air compression coefficient;

[0194] R is the air constant, 286.7 J / (kg·K);

[0195] u is the internal energy of the air inside the cavern, J;

[0196] The heat exchange rate between the air inside the chamber and the chamber wall is expressed in J / s.

[0197] Combining the mass conservation equation in equation (4) and the convective heat transfer equation in equation (5), we obtain the analytical solution (6) for the gas temperature inside the chamber.

[0198]

[0199]

[0200] In the formula:

[0201] h c The convective heat transfer coefficient between the air inside the chamber and the chamber wall is expressed in W / (m²). 2 ·K);

[0202] A c The surface area of ​​the cavern is m. 2 ;

[0203] r0 is the radius of the cavern, in meters;

[0204]

[0205] In formula (6):

[0206]

[0207]

[0208]

[0209] in:

[0210] ρ av To calculate the average air density (kg / m³) within a time step Δt. 3 ;

[0211] ρ0 is the average air density inside the tunnel during the previous time period, in kg / m³. 3 ;

[0212] T0 is the air temperature inside the chamber at the previous time, in K;

[0213] m i The inflation rate for this stage is kg / s;

[0214] m e The pumping rate for this stage is kg / s;

[0215] V is the volume of the chamber, in meters. 3 ;

[0216] T m The temperature of the sealing layer is K.

[0217] Furthermore, when calculating the fundamental solution for heat conduction in surrounding rock, using axisymmetric heat conduction will significantly increase the computational difficulty. The following is an example of an axisymmetric heat conduction calculation method:

[0218] When the temperature of the sealing layer is not equal to the temperature of the air inside the chamber, convective heat transfer occurs, and the heat conduction equation of the surrounding rock is as follows (26):

[0219]

[0220] In the formula:

[0221] k u is the thermal conductivity of the medium in the tunnel wall, W / (m·K);

[0222] ρ u The density of the medium in the tunnel wall is kg / m³. 3 ;

[0223] c pu The isobaric specific heat of the medium in the tunnel wall, J / (kg·K);

[0224] The initial and boundary conditions for solving equation (10) are as follows:

[0225] T u (r, t) = f(r), t = 0 (27)

[0226] T u (r, t) = T R0 , r = r3 (28)

[0227] T u (r, t) = T m , r = r1 (29)

[0228] In the formula:

[0229] r1 is the outer radius of the sealing layer; r3 is the radius of the temperature or pressure influence range; the thermophysical parameters of the thick-walled cylinder between r1 and r3 are taken from the parameters of the concrete lining;

[0230] f(r) is a function of the surrounding rock temperature at a certain calculation start time point, which is determined by the calculation results of the previous time period.

[0231] T R 0 represents the temperature of the rock outside the radius r3, which is a constant in K.

[0232] T m The temperature of the sealing layer of the chamber is assumed to be constant within a short calculation step Δt, in K.

[0233] Applying the Laplace transform to both sides of equation (26), we get:

[0234]

[0235] In the formula:

[0236] s is the Laplace operator;

[0237] a u It is a constant.

[0238] T u (r, 0) represents the initial temperature distribution of the surrounding rock. From equation (27), we know that T u (r, 0) = f(r);

[0239] For T u The Laplace transform of is a function of r and s. For T u The relationship is shown in the following equation (32):

[0240]

[0241] Rearranging equation (30) yields:

[0242]

[0243] Equation (33) is a zero-order non-homogeneous Bessel equation with imaginary principal variables, and the general solution of the corresponding homogeneous equation is:

[0244]

[0245] In the formula:

[0246] C1 and C2 are integration constants;

[0247] I0 and K0 are Bessel functions of the first and second kind of zero-order imaginary principal variables, respectively;

[0248] Using the method of variation of constants, a particular solution of equation (33) can be obtained as follows:

[0249]

[0250] Therefore, the general solution of equation (33) is:

[0251]

[0252] After performing a Laplace transform on the boundary conditions (28) and (29), substituting them into equation (36) and rearranging, we get:

[0253]

[0254] In the formula:

[0255]

[0256]

[0257]

[0258]

[0259]

[0260]

[0261] Thus, the Plas transform of the temperature distribution function in the surrounding rock under the convective heat transfer mode is obtained.

[0262] Plas transform of temperature distribution function in surrounding rock The original function T can be obtained by further performing the inverse Plas transform of equation (38). u (r,t):

[0263]

[0264] Considering the relatively short time step in the temperature function calculation, equation (38) is solved numerically using the following equation (39):

[0265]

[0266] In the formula,

[0267] In equation (39), in principle, the larger the value of N, the more accurate the calculation. However, due to the influence of rounding errors, N is generally taken as an even number between 8 and 20 when solving equation (39). When N is 18, the accuracy of the calculation result reaches 10. -5 ~10 -7 It is orders of magnitude faster and has a relatively fast computation speed.

[0268] However, in reality, compressed air energy storage chambers are generally quite large in diameter. For tunnel-type chambers, the diameter is typically 10–14 meters, while for large tank-type chambers, the diameter can reach 40 meters. Because the thermal conductivity of the surrounding rock is relatively low, the temperature change inside the chamber has a limited impact on the surrounding rock temperature. Therefore, for such large-diameter chambers, a one-dimensional heat conduction model can be used for heat conduction in the chamber walls. This not only greatly simplifies the calculations and ensures the stability of the numerical calculations but also guarantees sufficient accuracy. The calculation model is as follows: Figure 9 As shown in the figure. In this model, the temperature of the one-dimensional heat conduction outer boundary is isothermal, that is, the initial temperature T of the surrounding rock. R0 .

[0269] In calculating the fundamental solution of one-dimensional heat conduction in the surrounding rock, let the radius of the cave be r1, the radius of the surrounding rock within the temperature influence range be r3, the temperature of the outer boundary of the one-dimensional heat conduction be constant, and let the x-axis be extended outward from a point on the cave wall as the starting point; the boundary value problem under the one-dimensional heat conduction condition is as follows:

[0270]

[0271] T = T m x = 0, t > 0 (11)

[0272] T = T R0 x = L, t > 0 (12)

[0273] T = f(x) t = 0 (13)

[0274] In the formula, T R0 R is the initial temperature of the surrounding rock; L is the distance of one-dimensional heat conduction, L = r3 - r1.

[0275] Specifically, the above boundary value problem can be broken down into the following two problems:

[0276] First, there is the problem of steady-state nonhomogeneity:

[0277]

[0278] T w =T m x = 0 (15)

[0279] T w =T R0 x=L (16)

[0280] The solution to this definite problem is:

[0281]

[0282] Second, unsteady homogeneous problems:

[0283]

[0284] T f =0 x=0, t>0 (19)

[0285] T f =0 x=L,t>0 (20)

[0286] T f =f(x)-T w (x) t=0 (21)

[0287] The solution to this definite problem is:

[0288]

[0289] The solutions to the original problem equations (10) to (13) are:

[0290] T(x,t)=T w (x)+T f (x, t) (23)

[0291] Substituting equations (17) and (22) into equation (23) and rearranging, we get:

[0292]

[0293] In the formula,

[0294] Similar to the one-dimensional finite-thickness outer boundary adiabatic model, the solution to the boundary value problem contains an infinite series summation operation. In actual calculations, the cutoff value can be determined through trial and error. In this scheme, the maximum value of m is taken as 20.

[0295] Furthermore, in a preferred embodiment, a shorter time step results in a longer total computation time and is more prone to numerical computation difficulties; conversely, a longer time step leads to a larger computation error but a shorter computation time. Actual calculations show that when the time step length Δt is between 10s and 600s, the accuracy of the calculation results meets the requirements.

[0296] Here is a specific example:

[0297] In this example, the initial inflation is performed at a predetermined inflation rate, from one atmosphere to the set minimum pressure of 8.5 MPa, followed by a 3-day rest period before entering the standard cycle condition. Through trial calculations, the initial inflation time under this mode is found to be 18.32 hours. Other calculation parameters are shown in Table 1.

[0298]

[0299] Table 1

[0300] Based on the above parameters, the temperature variations over time for the air inside the chamber, the sealing layer, and two points on the chamber wall are calculated. (See attached curves.) Figure 10 The results showed that the temperature rose rapidly during the initial inflation phase, especially at the beginning. During the lull phase after the initial inflation, the temperature gradually decreased over time. After entering the standard cycle, the air temperature varied with the inflation-high-pressure storage-evacuation-low-pressure storage cycle. The temperature rose sharply during the inflation phase, decreased slowly during the high-pressure storage phase, decreased sharply during the evacuation phase, and rose slowly during the low-pressure storage phase. The temperature trend after entering the standard cycle showed a general downward trend with increasing cycle count, initially rapid but slowing down towards a relatively stable state. The further away from the chamber wall, the more pronounced the temperature lag effect.

[0301] The curve showing the change in air pressure inside the chamber over time is shown below. Figure 11 The results showed that during the inflation phase, the pressure inside the chamber increased linearly over time. During the 3-day rest period, the pressure decreased slightly. After entering the standard cycle, the air pressure varied with the cycles of inflation, high-pressure storage, extraction, and low-pressure storage. The pressure increased sharply during the inflation phase, decreased slowly during the high-pressure storage phase, decreased sharply during the extraction phase, and increased slowly during the low-pressure storage phase. The temperature trend after entering the standard cycle showed that the overall pressure decreased with increasing cycle count, initially rapidly, then slowing down until reaching a relatively stable stage.

Claims

1. A coupled calculation method for the thermodynamics of a cavern and the thermal conduction of its walls based on one-dimensional heat conduction, characterized in that... Includes the following steps: S1: Determine the geometric dimensions and operating parameters of the chamber, provide the basic thermodynamic and mechanical calculation parameters for the sealing layer, lining, and surrounding rock, and determine the total calculation time t. n The time domain to be calculated is discretized, the time step length Δt is divided, and initial temperature values ​​are assigned to the air and surrounding rock in the cave. S2: Calculate the air inflation and deflation rates, inflation temperature, and air density in the chamber at each time step; S3: Based on the convective heat transfer between the tunnel wall and the air and the conservation of heat conduction energy in the tunnel wall, and based on the air temperature and tunnel wall temperature distribution at time step i, calculate the temperature T of the sealing layer at time step i+1. m ; S4: Assuming the temperature of the sealing layer remains constant in the (i+1)th time step, calculate the temperature of the air inside the chamber in the (i+1)th time step, and use a one-dimensional heat conduction model to calculate the heat conduction of the chamber wall to obtain the temperature distribution of the chamber wall. S5: Repeat steps S3-S4 until the total computation time is greater than the total calculation time; In each time step of the calculation, the sealing layer temperature T is calculated based on the fact that the convective heat transfer power between the air and the wall of the chamber is equal to the heat conduction power on the surface of the wall obtained based on the temperature distribution. m and with the sealing layer temperature T m As boundary conditions, thermodynamic calculations of the air in the chamber and one-dimensional heat conduction calculations of the chamber wall are performed separately. The temperature and pressure of the air in the chamber are obtained through thermodynamic calculations of the air in the chamber, and the radial temperature distribution of the chamber wall is obtained through one-dimensional heat conduction calculations of the chamber wall. In calculating the fundamental solution of one-dimensional heat conduction in the surrounding rock, let the radius of the cave be r1, the radius of the surrounding rock within the temperature influence range be r3, the temperature of the outer boundary of the one-dimensional heat conduction be constant, and let the x-axis be extended outward from a point on the cave wall as the starting point; the boundary value problem under the one-dimensional heat conduction condition is as follows: In the formula, T R0 R is the initial temperature of the surrounding rock; L is the distance of one-dimensional heat conduction, L = r3 - r1.

2. The method for coupled calculation of tunnel thermodynamics and tunnel wall heat conduction based on one-dimensional heat conduction according to claim 1, characterized in that: Discretize the time domain to be calculated and select a time step Δt; Let coordinate axis r extend outward from the center of the cave, and let the discrete interval of coordinate axis r be Δr; Assuming the temperature T of the sealing layer in the time-step chamber is... m Constant and does not change with the thickness or location of the sealing layer; This is used to calculate the gas temperature T inside the chamber and the temperature gradient of the surrounding rock at the end of the first time step. ,Right now Substituting it into equation (1) yields equation (2): In the formula: k is the thermal conductivity of the concrete lining, W·(m·K). -1 ; h c The convective heat transfer coefficient between the air inside the chamber and the chamber wall is expressed in W / (m²). 2 ·K); T u The temperature of the surrounding rock; In the formula: The convective heat transfer power of the surrounding rock at the end of the previous time interval, in W; The thermal conductivity of the surrounding rock is expressed in W. The required sealing layer temperature T for the next time interval is further calculated using equation (2). m .

3. The coupled calculation method for tunnel thermodynamics and tunnel wall heat conduction based on one-dimensional heat conduction according to claim 2, characterized in that, When calculating the fundamental solution for the gas temperature T inside the chamber, it is assumed that the chamber wall temperature is constant, the air density inside the chamber is constant, the filling rate or pumping rate is constant, and the air is an ideal gas. The energy conservation equation in equation (3) is further transformed into an ordinary differential equation. In the formula: V is the volume of the chamber, in meters. 3 ; ρ is the density of the air inside the chamber, in kg / m³. 3 ; c v The isochoric specific heat of air is 717 J / (kg·K); T is the temperature of the air inside the chamber, in K; t represents time, in seconds; m in (t) is the inflation rate function; m ex (t) is the pumping rate function; h i To accommodate the enthalpy of air, J; h is the enthalpy of the air inside the chamber, in J; Z is the air compression coefficient; R is the air constant, 286.7 J / (kg·K); u is the internal energy of the air inside the cavern, J; The heat exchange rate between the air inside the chamber and the chamber wall is expressed in J / s. Combined with the mass conservation equation in equation (4) and the convective heat transfer equation in equation (5), we obtain the analytical solution equation (6) for the gas temperature inside the chamber. In the formula: h c The convective heat transfer coefficient between the air inside the chamber and the chamber wall is expressed in W / (m²). 2 ·K); A c The surface area of ​​the cavern is m. 2 ; r0 is the radius of the cavern, in meters; In formula (6): in: ρ av To calculate the average air density (kg / m³) within a time step Δt. 3 ; ρ0 is the average air density inside the tunnel during the previous time period, in kg / m³. 3 ; T0 is the air temperature inside the chamber at the previous time, in K; m i The inflation rate for this stage is expressed in kg / s. m e The pumping rate for this stage is kg / s; V is the volume of the chamber, in meters. 3 ; T m The temperature of the sealing layer is K.

4. The method for coupled calculation of tunnel thermodynamics and tunnel wall heat conduction based on one-dimensional heat conduction according to claim 1, characterized in that, The above boundary value problem can be broken down into the following two problems: First, there is the problem of steady-state nonhomogeneity: The solution to this definite problem is: Second, unsteady homogeneous problems: The solution to this definite problem is: The solutions to the original problem equations (10) to (13) are: Substituting equations (17) and (22) into equation (23) and rearranging, we get: In the formula, .

5. The coupled calculation method for tunnel thermodynamics and tunnel wall heat conduction based on one-dimensional heat conduction according to claim 4, characterized in that, The maximum value of m is 20.

Images