Gas storage volume design method for large-scale compressed air energy storage power station

Abstract

This invention discloses a method for designing the volume of a large-scale compressed air energy storage power station gas storage reservoir. This method constructs a memory-corrected state equation that reflects the impact of historical operations, establishes a dynamic heat transfer analysis that simulates energy transmission along discrete paths, and uses multivariable coupled iterative solution with the ultimate optimization goal of minimizing the system's total entropy production rate to determine the optimal gas storage reservoir volume. This method synergistically couples a fluid memory effect model, a path topology dynamic heat transfer model, and an entropy production rate optimization model to achieve global, dynamic optimization of the gas storage reservoir volume, significantly improving the accuracy, safety, and cost-effectiveness of large-scale gas storage reservoir designs.

Description

Technical Field

[0001] The present invention relates to a compressed air energy storage technology, in particular to a method for designing the volume of a gas storage reservoir in a large-scale compressed air energy storage power station. Background Art

[0002] Million-kilowatt-class large-scale compressed air energy storage (CAES) power plants are key technologies for supporting high-proportion renewable energy integration and ensuring the safe and stable operation of power systems. Large underground gas storage facilities, as core equipment, are crucial for long-term, safe, and efficient operation. During gas storage operation, compressed air inevitably contains water vapor. Improper humidity control can cause serious problems on multiple levels. High-pressure, wet air not only reacts with surrounding rocks like rock salt and sandstone, impacting the long-term geological stability of the gas storage, but also forms electrolyte solutions on the surfaces of metal components. Combining with carbon dioxide and sulfides in the air, they create an acidic environment, causing severe electrochemical corrosion and seriously threatening the service life and safety of valves, pipelines, and critical pressure-bearing equipment. Furthermore, during rapid energy release, the gas undergoes adiabatic expansion, causing a sudden drop in temperature. Supersaturated water vapor or condensed water droplets can form hydrates or ice on throttle valves or turbine blades, causing system blockages and even equipment damage.

[0003] Currently, engineering practices primarily focus on controlling compressed air humidity by combining source dehydration with process monitoring. Technically, source dehydration generally employs established industrial drying technologies. For example, cooling dehydration utilizes a freeze dryer to forcibly cool compressed air to below a preset pressure dew point, condensing most of the water vapor into liquid water that is removed by a separator. For scenarios requiring extremely low humidity, adsorption drying is further employed. This utilizes the physical adsorption properties of desiccants such as silica gel, activated alumina, or molecular sieves to deeply capture water molecules in the air, reducing the pressure dew point to -60°C or even lower. For process monitoring, industrial-grade online humidity monitors or dew point meters are typically installed at key points within the gas storage facility, such as injection ports, production outlets, and specific locations within the facility, to provide real-time, fixed-point, and periodic monitoring of gas humidity. Operationally, front-end drying equipment is often manually adjusted based on seasonal changes in ambient humidity, combined with monitoring data. Alternatively, a fixed, conservative drying target is maintained year-round to ensure condensation-free operation under all operating conditions.

[0004] However, while existing technologies can provide basic humidity assurance, they are difficult to apply to large-scale, multi-million-kilowatt gas storage facilities with complex and dynamic operating conditions. For example, existing methods treat the gas storage facility as an ideal container with uniform humidity (i.e., a lumped parameter model). In reality, the internal temperature and flow field are extremely non-uniform, easily leading to the formation of localized high-humidity zones or condensation-dominated paths. Due to the lack of accurate prediction capabilities for these localized condensation risk points, existing technologies can only adopt a global over-drying strategy, deeply dehydrating all gases, which results in energy waste. Summary of the Invention

[0005] The purpose of the invention is to provide a method for designing the volume of a gas storage reservoir for a large-scale compressed air energy storage power station, in order to solve the technical problems existing in the prior art.

[0006] Technical solution: A method for designing the volume of a large-scale compressed air energy storage power station gas storage reservoir, including:

[0007] Receive design requirement data and determine initial design parameters accordingly;

[0008] Based on historical operating data and initial design parameters, a fluid memory effect model is constructed to generate a memory-corrected state equation for subsequent entropy production rate calculation;

[0009] Based on the initial design parameters, a path topology dynamic heat transfer model is established to obtain the time-varying equivalent heat transfer coefficient for subsequent entropy production rate calculation;

[0010] The total entropy production rate of the system, which represents the thermodynamic irreversibility of the system, is calculated by coupling the memory-corrected state equation with the time-varying equivalent heat transfer coefficient.

[0011] Iterative optimization is performed to search for the gas storage volume that minimizes the total entropy production rate of the system, and this volume is determined as the final design volume.

[0012] Beneficial effect: The present invention synergistically couples the fluid memory effect model, the path topology dynamic heat transfer model and the entropy production rate optimization model to achieve global and dynamic optimization of the gas storage volume, significantly improving the accuracy, safety and economy of the design of large-scale gas storage. BRIEF DESCRIPTION OF THE DRAWINGS

[0013] Figure 1 It is a flow chart of the present invention.

[0014] Figure 2 It is a flow chart of establishing a path topology dynamic heat transfer model of the present invention.

[0015] Figure 3 It is a flow chart of the present invention for generating a heat transfer path network.

[0016] Figure 4 It is a flow chart of forming a dynamic path weight vector according to the present invention.

[0017] Figure 5 It is a flow chart of the present invention for obtaining a set of path-specific heat transfer characteristic data. DETAILED DESCRIPTION

[0018] In order to make the purpose, technical solutions and advantages of the present invention clearer, the following Figures 1 to 5 , and specific embodiments are provided to further describe the present invention in detail. It should be noted that the specific embodiments herein are only used to explain the present invention and are not intended to limit the scope of protection of the present invention. Any modifications, equivalent substitutions, improvements, etc. made within the spirit and principles of the present invention shall be included in the scope of protection of the present invention.

[0019] During the research process, the applicant found that the main problems are as follows:

[0020] Existing humidity control methods suffer from two major issues: First, they cannot accurately predict the distribution of local condensation, forcing them to resort to a global over-dehydration strategy, which results in excessive energy consumption. Specifically, existing methods treat the gas storage as a uniform container and employ a global average humidity control strategy. However, in reality, the temperature and flow fields within the gas storage are extremely non-uniform, making localized areas prone to high humidity and even condensation. Existing methods cannot accurately predict the spatial and temporal distribution of these localized condensation points. To mitigate this risk, they resort to over-dehydration of the entire gas, resulting in unnecessary high energy consumption.

[0021] Second, the humidity memory effect is ignored, making it impossible to adjust strategies based on historical conditions, resulting in insufficient control precision. This makes it difficult to strike a balance between safety and economy. Specifically, existing control strategies are static responses that ignore the humidity memory effect of the gas storage reservoir. Residual moisture in the system will re-enter the gas phase in subsequent cycles, changing the initial humidity conditions. However, existing methods cannot proactively adjust drying strategies based on historical moisture conditions, resulting in insufficient control precision and potential condensation risks.

[0022] Example 1: A method for designing the volume of a gas storage reservoir applicable to a million-kilowatt compressed air energy storage power station is provided. When executed on a computing device (such as a server or workstation), the method includes the following steps:

[0023] Step S101: receiving design requirement data and determining initial design parameters accordingly.

[0024] In this embodiment, the design requirement data is the input basis for starting the entire design process, which is mainly derived from the overall design objectives of the power station. Specifically, it at least includes the design power generation power P of the power station. design (e.g. 1000MW), rated energy storage duration t storage (e.g. 8 hours), and geological exploration data of the gas storage location data(For example, the target formation is a deep salt cavern.) Based on these inputs, combined with engineering experience and relevant design specifications, the type of gas storage can be preliminarily determined (for example, an underground gas storage using an abandoned salt mine), its basic shape can be selected as a nearly vertical cylinder, and a series of initial design parameters can be set to form an initial set, such as: the initial temperature T0 of the gas storage (for example, 300K), the minimum operating pressure p min (such as 5MPa), maximum operating pressure p max (such as 10MPa), temperature range T min =a℃,T max = b ° C. The shape factor k, etc., characterizes the geometric characteristics of the gas storage reservoir.

[0025] Step S102: constructing a fluid memory effect model based on historical operating data and initial design parameters to generate a memory-corrected state equation for subsequent entropy production rate calculation.

[0026] In this embodiment, this step is primarily used to modify the traditional gas equation of state to account for the historical dependence of compressed air during the cycle, known as the memory effect. The basic principle is that after a gas undergoes multiple non-ideal compressions and expansions, its intermolecular forces and energy distribution leave traces, causing its behavior at the same PVT point to differ from that of gas that has never undergone a cycle, affecting its current thermodynamic state. The model analyzes historical operating data (time series of pressure and temperature) to quantify this historical influence and express it as a correction term, ultimately generating a more accurate memory-corrected equation of state that incorporates historical information. This equation is used to calculate the pressure drop entropy production rate because it more accurately predicts the actual pressure loss caused by throttling and friction.

[0027] Historical operating data is used to capture deviations between real-world gas behavior and the ideal model and serves as the source of "memory." This is achieved by collecting long-term series of pressure (p_raw(t)), temperature (T_raw(t)), and other data from gas storage facilities with similar operating conditions or high-fidelity CFD simulation databases.

[0028] Step S103: Based on the initial design parameters, a path topology dynamic heat transfer model is established to obtain a time-varying equivalent heat transfer coefficient for subsequent entropy production rate calculation.

[0029] In this embodiment, in order to solve the problem of overly rough heat transfer coefficient values ​​in traditional designs, the complex heat transfer process inside the gas storage reservoir is deconstructed into a network consisting of a large number of discrete, parallel heat transfer paths. Each path has independent, calculable heat transfer properties. More importantly, the importance (i.e., weight) of the path will be dynamically adjusted as the flow field (such as pressure and velocity distribution) inside the gas storage reservoir changes in real time. By dynamically weighting and summing the heat transfer contributions of all paths, a time-varying, high-precision equivalent heat transfer coefficient h is finally obtained.eff (t) for subsequent calculation of heat transfer entropy production rate. These data can accurately reflect the real, unsteady heat transfer process between the gas storage wall and the gas.

[0030] Different from the traditional lumped parameter method or full-scale CFD, the novel heat transfer analysis model, the path topology dynamic heat transfer model can capture the main contradictions and time-varying characteristics of the heat transfer process at a lower computational cost.

[0031] The time-varying equivalent heat transfer coefficient h output by the model eff (t), which replaces the fixed empirical value in traditional design and will serve as the key basis for subsequent entropy production rate calculations (especially heat transfer entropy production calculations).

[0032] Step S104 , coupling the memory-corrected state equation with the time-varying equivalent heat transfer coefficient to calculate the total entropy production rate of the system that characterizes the thermodynamic irreversibility of the system.

[0033] In this embodiment, this step is the core of thermodynamic optimization. Entropy production is a measure of the energy dissipation and irreversibility of the system. The total entropy production rate S total It consists of three parts: flow entropy production rate S flow (caused by fluid viscous dissipation), heat transfer entropy production rate S heat (caused by finite temperature difference heat transfer) and pressure drop entropy production rate S pressure (Caused by pressure drop due to throttling, friction, etc.) Calculate S pressure When calculating S heat When , a time-varying equivalent heat transfer coefficient is needed to obtain a more accurate heat flux. In this way, the two core models are tightly coupled to the calculation of entropy production rate.

[0034] Step S105 , searching for the gas storage volume that minimizes the total entropy production rate of the system through an iterative optimization process, and determining this volume as the final design volume.

[0035] In this embodiment, in order to find the optimal gas storage volume V final The principle is that the volume of the gas storage V is proportional to the total entropy production rate S of the system. total There is a nonlinear constraint relationship between them: if the volume is too small, the flow rate is fast, and the flow and pressure drop entropy generation dominate; if the volume is too large, the wall heat exchange area is large, the heat transfer entropy generation may increase, and the construction cost rises sharply. Therefore, there must be an optimal volume that minimizes the total entropy generation rate. By establishing S total The functional relationship between V and the volume is obtained, and a numerical optimization algorithm (such as gradient descent method or genetic algorithm) is used to iteratively solve the problem while satisfying all engineering constraints (such as power, safety, and cost) to find the volume value corresponding to the minimum entropy production rate.

[0036] Through the above steps, this method organically combines the fluid memory effect, dynamic heat transfer and the second law of thermodynamics, breaking through the limitations of over-simplified models and empirical dependence of parameters in traditional design methods, thereby being able to design a million-kilowatt gas storage volume that achieves the best balance between safety, economy and efficiency.

[0037] In other words, steps S104 and S105 together form a closed optimization loop based on the second law of thermodynamics. The underlying principle is that any real thermodynamic process is accompanied by an increase in entropy (i.e., entropy production), and the magnitude of entropy production is a measure of the system's energy dissipation and irreversible losses. The optimal system design should minimize the total entropy production rate while still achieving its intended function.

[0038] The specific calculation process can be as follows: starting with a candidate volume V, using the outputs of steps S102 and S103, calculate the total entropy generation rate of the system within a complete operating cycle at that volume. An optimization algorithm (such as gradient descent or genetic algorithm) is used to adjust the volume V, and the entropy generation rate is repeatedly calculated until the volume value that minimizes the entropy generation rate is found. This is the final design volume sought by the present invention. Because it is obtained by comprehensively considering the fluid memory effect, dynamic heat transfer effects, and the system's global thermodynamic optimality, it is more accurate and efficient.

[0039] In summary, this embodiment can significantly improve the accuracy and reliability of the volume design of large-scale gas storage facilities, avoiding the huge design margins and potential economic and safety issues caused by the traditional method's reliance on static empirical values.

[0040] Example 2: As a refinement and optimization of step S102 in Example 1, this example describes in detail the specific process of constructing a fluid memory effect model and generating a memory-corrected state equation.

[0041] Step S201 : Processing historical operation data and constructing a standardized state trajectory database containing multiple complete inflation and deflation cycles.

[0042] In this embodiment, the historical operation data can be obtained from the actual operation records of similar gas storage facilities, or generated through high-fidelity simulation. For example, the pressure p is collected for one month with a sampling frequency of 1 Hz. raw (t), temperature T raw (t) and density ρ raw (t) Data.

[0043] The standardized state trajectory refers to the curve formed in the state space (such as the PT coordinate system) after normalizing the pressure, temperature and other data of the complete inflation and deflation cycle (from the start of inflation to the end of deflation).

[0044] First, a data cleaning algorithm (such as the 3σ criterion or isolation forest) is applied to remove outliers and noise points from the raw data. Then, by analyzing the sign of the first-order derivative of pressure p(t) (dp / dt), the start and end times of the inflation phase (dp / dt>ε), storage phase (|dp / dt|≤ε), and deflation phase (dp / dt<-ε) are automatically identified. Finally, each identified cycle data is normalized (for example, mapped to the [0, 1] interval) to form a standardized state trajectory, Trajectory_i, which is stored in the database.

[0045] Intermediate results: The output of this step is a structured historical trajectory feature library, which contains standardized trajectory data of N complete cycles, providing a basis for subsequent memory feature quantification.

[0046] Step S202 calculates the cumulative impact of historical cycles on the current gas state based on the standardized state trajectory database and constructs this into a total memory function. Prior to this, a step of assigning weights to each historical cycle is preferably included. In this embodiment, this step is the core of quantitative memory.

[0047] Step S202a, generating a comprehensive memory weight vector that assigns weights to each historical cycle.

[0048] Specifically, it can be as follows: based on the time interval between each historical cycle and the current moment, a time weight that decreases as the time distance increases is calculated for each historical cycle through a time decay function;

[0049] Calculate the geometric similarity between the state trajectory of each historical cycle and the current state trajectory, and assign trajectory similarity weights to each historical cycle based on the geometric similarity;

[0050] The time weight and trajectory similarity weight are weighted and fused to generate a comprehensive weight for each historical cycle. The comprehensive weights of all historical cycles together constitute a comprehensive memory weight vector.

[0051] Since not all historical cycles are equally important, intuitively, cycles that are more recent in time and have more similar operating conditions should have a greater impact. Each element in the generated comprehensive memory weight vector represents the influence of a historical cycle on the current state.

[0052] Based on the time interval between each historical cycle and the current moment, the time decay function w time (i)=exp(-λ*(t current -t cycle_i )) / τ decay ) to calculate the time weight. Among them, λ is the attenuation coefficient (for example, 0.5), τ decayis the characteristic decay time (e.g. 7 days). This weight decreases exponentially as the time distance increases. In some embodiments, the decay coefficient λ can be set to 0.1, and the characteristic decay time τ decay Can be set to 1 day (86400 seconds), which means that loops from a day ago will have significantly less influence.

[0053] Calculate the state trajectory of each history cycle i Trajectory with current (or target) state current Specifically, the Euclidean distance or dynamic time warping (DTW) distance between the two curves can be used and substituted into the Gaussian kernel function w similarity (i)=exp(-||T rajectorycurrent -T rajectoryi ∣∣2 2 / σ 2 similarity ) to calculate the trajectory similarity weight. The higher the similarity, the smaller the distance, the greater the weight.

[0054] Perform weighted fusion on the above two weights to obtain the comprehensive weight w i =α×w time (i)+β×w similarity (i), where α and β are fusion coefficients (e.g., α = 0.4, β = 0.6, indicating that the influence of time recency is slightly greater than that of morphological similarity). The combined weights of all historical cycles together form the comprehensive memory weight vector W memory .

[0055] Step S202b: Constructing a total memory function. Specifically, based on the standardized state trajectory database, constructing a molecular interaction memory kernel function for characterizing intermolecular interaction memory, a compression history memory function for characterizing compression history memory, and a temperature history memory function for characterizing temperature history memory;

[0056] The total memory function is formed by nonlinearly weighted coupling of the molecular interaction memory kernel function, the compression history memory function and the temperature history memory function.

[0057] In this embodiment, a dual-factor weighting approach is used to identify historical cycles that have the greatest impact on the current state, typically those that are close (time weight) and similar (similarity weight). This intelligently filters the most valuable information from massive amounts of historical data, improving the relevance and accuracy of the memory model.

[0058] Under some working conditions, the specific calculation process is as follows:

[0059] Constructing molecular interaction memory kernel function K molecular: This function characterizes the interaction memory of gas molecules under non-ideal conditions. It is constructed by weighted summing the product of the deviations between the actual state and the ideal state over the history cycle. It is used to characterize the molecular force memory caused by the change in molecular spacing due to history compression.

[0060] Constructing the compressed history memory function M compression This function characterizes the impact of a history of severe pressure fluctuations on the gas state. It is constructed by weighted summing characteristic parameters such as the maximum pressure ratio and pressure change rate during the historical cycles. It is used to characterize the effect of cyclic pressure intensity on gas mechanical fatigue.

[0061] Construct temperature history memory function M thermal This function characterizes the impact of the temperature field change history. It is constructed by weighted summing characteristic parameters such as the maximum temperature difference and temperature change rate in the historical cycle. It is used to characterize the impact of cyclic temperature shocks on the thermal history of the gas.

[0062] Nonlinear weighted coupling: The three sub-functions mentioned above are combined into a total memory function through nonlinear coupling: M total (p history , T history )=γ1×K molecular +γ2×M compression +γ3×M thermal The coupling coefficients γ1, γ2, and γ3 can be obtained by fitting the experimental data.

[0063] Step S203 : The total memory function is used as a correction factor and introduced into the reference gas state equation to form a memory correction state equation.

[0064] The traditional ideal gas state equation is p×V=n×R×T. The total memory function M obtained in step S202b is total Introduced as a correction term, the memory-corrected state equation is obtained: p×V=n×R×T×[1+α memory ×M total (p history , T history )] Among them, α memory It is the overall strength coefficient of the memory effect calibrated through experiments.

[0065] This embodiment more accurately describes the thermodynamic behavior of gases in large-capacity gas storage facilities under complex operating conditions than traditional equations. For example, after multiple consecutive high-intensity inflations, even if the current pressure and temperature are the same as after a lower-intensity inflation, the equation can predict different gas density and enthalpy values ​​due to the memory function. This lays the foundation for accurate energy calculations and entropy generation analysis.

[0066] Example 3: As a refinement and optimization of step S103 in Example 1, this example describes in detail the specific process of establishing a path topology dynamic heat transfer model and obtaining a time-varying equivalent heat transfer coefficient.

[0067] Step S301 : Based on the geometry of the gas storage defined by the initial design parameters, the physical gradient field of the wall is analyzed, and a heat transfer path network consisting of multiple discrete paths is generated accordingly.

[0068] Specifically, this can be achieved by: in the physical gradient field, selecting grid points where both the temperature gradient and the height gradient exceed respective preset thresholds, and aggregating them into a set of path starting points;

[0069] Taking each point in the path starting point set as the source point, an iterative tracking algorithm is executed along the composite direction of the temperature gradient and the height gradient to generate a coordinate sequence of the path point by point until the preset termination condition is met; the coordinate sequences of multiple paths together constitute a heat transfer path network.

[0070] The heat transfer path network is a virtual topological structure used to describe how energy is transferred within a gas storage reservoir. It consists of paths representing energy transfer channels and nodes at their intersections. Specifically, it is represented as a discretized data structure that describes the primary pathways for energy transfer from high-temperature to low-temperature zones within the gas storage reservoir. It primarily consists of a sequence of node coordinates representing the paths. In this embodiment, the physical gradient field refers to the temperature gradient gradT and the surface geometric height gradient gradh at each point on the gas storage reservoir wall grid. This is described in detail below.

[0071] In one working condition, the following can be done: Based on the volume V and shape coefficient k in the initial design parameters, a 3D geometric model of the gas storage reservoir is established. Then, the inner wall surface of the gas storage reservoir is divided into an M×N grid.

[0072] Calculate the physical gradients at each grid point, mainly the temperature gradient gradT (obtained from CFD simulation or sensor data) and the geometric height gradient gradh (representing the direction of gravity). Filter out those points whose gradient modulus is greater than a preset threshold (e.g. |gradT|>T threshold And |gradh|>h threshold ), as the seed point where energy transfer is most likely to occur, forming a set of path starting points.

[0073] Starting from each seed point, along the composite direction of temperature gradient and height gradient Direction=wT×gradT+w h×gradh performs path tracing. The algorithm iteratively calculates the coordinates of the next point with a set step size δs until the path reaches the gas storage boundary, the gradient disappears, or the preset maximum length is reached, thereby generating a coordinate sequence for each path. The coordinate sequences of all paths together constitute the heat transfer path network.

[0074] In a certain scenario, the above process is as follows:

[0075] First, based on the geometry of the gas storage reservoir (for example, a cylindrical salt cavern with a diameter of 20 meters and a height of 100 meters), the finite element method is used to divide its three-dimensional wall into a uniform grid of 100x200. Then, the physical gradient field is calculated at each grid point.

[0076] Among all grid points, those whose temperature gradient |gradT| exceeds the preset threshold of 0.5K / m and whose height gradient |∇h| exceeds the preset threshold of 0.01 are selected, and these points are aggregated into a set of path starting points.

[0077] Taking each point in the set as the source point, the composite direction of the temperature gradient and the height gradient (for example, Direction=0.7×gradT norm +0.3×gradh norm ) performs an iterative tracking algorithm. This algorithm generates a coordinate sequence for each path, point by point, with a step size of 0.1 meters, until the path reaches the gas storage boundary or the gradient falls below a termination threshold. The coordinate sequences of all paths together constitute the heat transfer path network.

[0078] Through a dual-gradient field-based path generation method, advantageous heat transfer pathways typically emerge in areas of greatest temperature difference (dominated by gradT), while also being guided and constrained by wall geometry (e.g., ridges and depressions, dominated by gradh). Intelligently identifying these energy highways, driven by both thermodynamics and geometric constraints, significantly reduces computational complexity compared to global CFD calculations while preserving critical physical information.

[0079] In step S302 , the geometric and thermodynamic properties of each path in the heat transfer path network are independently calculated to obtain a set of path-specific heat transfer characteristic data.

[0080] Specifically, it can be: analyzing the geometric characteristics of each path in the heat transfer path network, and classifying each path into one of the predefined path types according to a preset geometric classification criterion;

[0081] For each path, the heat transfer correlation formula that matches the predefined path type of the path is analyzed and adopted from the heat transfer correlation formula, and combined with the internal flow field state, the heat transfer properties of the path are calculated;

[0082] The heat transfer properties of all paths are combined into path-specific heat transfer characteristic data.

[0083] Predefined path types are categorized based on their geometric characteristics (e.g., curvature), such as direct, reflective, and bypass paths. Heat transfer correlations are empirical or semi-empirical formulas that describe the relationship between the Nusselt number (Nu) and dimensionless parameters such as the Reynolds number (Re) and the Prandtl number (Pr).

[0084] In some embodiments, the process is as follows: for each generated path Path i , calculate its geometric parameters, such as the actual length L actual_i , the straight-line distance L from the starting point to the end point straight_i , mean curvature k avg_i Based on these geometric features, paths can be automatically classified into different predefined types, such as: Direct path: small curvature, close to a straight line. Reflective path: one or more sudden changes in direction occur at the wall. Circumferential path: large curvature, winding path.

[0085] For example, calculate the mean curvature k of the path avg , if k avg <0.05m -1 , it is classified as a “direct path”; if k avg >0.2m -1 , it is classified as a bypass path; if it is between the two and there is a direction mutation point, it is classified as a reflection path.

[0086] For each path, according to its type, the best matching formula is selected from the heat transfer correlation library to calculate its heat transfer properties. For example, the local Reynolds number Re of the path is calculated. path_i ,Then:

[0087] For the direct path, the Nusselt number Nu can be calculated using the classic Dittus-Boelter formula: direct .

[0088] For reflection and bypass paths, Nu direct Multiply the correction factor that takes into account the reflection enhancement or turbulence enhancement on the basis of Nu reflection or Nu tortuosity .

[0089] Finally, the heat transfer coefficient h of each path is calculated path_i and thermal resistance R thermal_i .

[0090] For example, the calculation process can be:

[0091] Direct Path Nu direct =0.023×Re 0.8 ×Pr0.4 ×(1+path correction factor), where the path correction factor can be related to the wall roughness.

[0092] Reflection Path Nu reflection =Nu direct ×(1+reflection enhancement factor×number of reflections), where the reflection enhancement factor can be set to 0.1.

[0093] Flow path Nu tortuosity =Nu direct ×(L actual / L straight ) 0.6 , where L is the path length. Through the above calculations, a set of path-specific heat transfer characteristic data is finally obtained.

[0094] The output path feature database contains the type, geometric parameters and thermodynamic properties of each path.

[0095] This differentiated classification-before-calculation approach takes advantage of the significant differences in heat transfer mechanisms (such as laminar flow, turbulent flow, and boundary layer separation and reattachment) between different path types. By matching specific correlations to each path type, the approach improves the physical realism and accuracy of heat transfer calculations, representing a significant technological advancement compared to the one-size-fits-all lumped parameter approach.

[0096] Step S303 , in response to the real-time changes in the flow field state inside the gas storage, a dynamically changing weight is assigned to each path to form a dynamic path weight vector.

[0097] Specifically, the process is as follows: according to the internal flow field state, the local Reynolds number and the pressure difference at both ends of each path are evaluated, and based on this, the path activation probability is calculated to indicate whether the path is activated;

[0098] Extract the effective heat transfer cross-sectional area and thermal resistance of each path from the path-specific heat transfer characteristic data, and calculate its area weight and thermal resistance weight accordingly.

[0099] The path activation probability, area weight and thermal resistance weight are fused to generate a comprehensive path weight and normalize it to form a dynamic path weight vector.

[0100] During gas storage operation, gas flow is not uniform; some areas experience more active flow and heat transfer. The weighting mechanism is designed to quantify this real-time variation in activity.

[0101] In this embodiment, the path activation probability is used to describe the probability or flag of whether a certain path becomes an effective heat transfer channel under the current flow field conditions.

[0102] In some working conditions, the calculation process is as follows:

[0103] Calculate the path activation probability w activation_i :Calculate the local Reynolds number Re of each path based on real-time flow field data path_i and the pressure difference Δp i When the Reynolds number exceeds the critical value (indicating that the flow is turbulent) and the pressure difference is large enough (indicating that there is a driving force), the path is considered to be activated and its activation probability P active_i Set to 1, otherwise 0.

[0104] For example, if Re path If the Reynolds number is greater than the critical one (e.g. 2300), the path activation probability of the path is 1, indicating that the channel is in a turbulent state and heat transfer is active; otherwise, it is 0.1 (a weak heat conduction effect is retained).

[0105] Calculate the area weight w area_i : Effective heat transfer cross-sectional area A based on the path eff_i Calculation, the larger the area, the higher the weight. area_i =A eff_i / ∑[A eff_j ].

[0106] Calculate thermal resistance weight w resistance_i :Based on path thermal resistance R thermal_i Calculation, the smaller the thermal resistance, the easier the heat transfer, the higher the weight. resistance_i =(1 / R thermal_i ) / ∑[1 / R thermal_j ].

[0107] It should be noted that from the path-specific heat transfer characteristic data, the area weight of each path is proportional to the effective cross-sectional area, and the thermal resistance weight is proportional to the inverse of the thermal resistance.

[0108] Fusion generates comprehensive path weight: Fusion of the above three to obtain comprehensive weight W path_i =(w α area_i *w β resistance_i *w γ activation_i ) / normalization factor. The weights of all paths constitute the dynamic path weight vector. The weight exponent can be calibrated based on experimental data, for example, α=0.5, β=1.0, and γ=1.0.

[0109] In some embodiments, the path activation probability can be a continuous function between 0 and 1, rather than just a binary judgment. For example, a Sigmoid function P can be used. active =1 / (1+exp(-(Re-Re critical ) / scale)) to smooth the transition.

[0110] By weighted summing the heat transfer coefficients of all paths and their corresponding dynamic weights, the time-varying equivalent heat transfer coefficient h is finally obtained. eff (t).

[0111] In step S304, the path-specific heat transfer characteristic data is weighted and integrated in combination with the dynamic path weight vector to form a time-varying equivalent heat transfer coefficient.

[0112] Specific implementation method: The dynamic weight vector W obtained in step S303 is path_i (t) and the heat transfer coefficient h of each path obtained in step S302 path_i (t) Perform weighted summation to obtain the instantaneous equivalent heat transfer coefficient: h eff_instant (t)=∑[W path_i (t)*h path_i (t)].

[0113] In order to avoid computational noise, the instantaneous value can be smoothed by time to obtain a more stable h eff_filtered (t).

[0114] Through the above steps, h is based on the physical model and changes in real time with the working conditions. eff (t) replaces the fixed empirical value in traditional design (such as 10W / (m 2 •K), underground 5W / (m 2 •K)). This improves the accuracy of predictions of the thermodynamic behavior of gas storage, especially during the dynamic process of rapid conversion between filling and degassing.

[0115] Example 4: As a refinement of steps S104 and S105 in Example 1, the specific process of calculating the total entropy production rate of the system and performing iterative optimization to determine the final design volume is described in detail.

[0116] In step S401, based on the memory-corrected state equation and the time-varying equivalent heat transfer coefficient, the flow entropy generation rate, heat transfer entropy generation rate, and pressure drop entropy generation rate caused by fluid viscous dissipation, wall heat transfer, and inlet and outlet pressure drops are calculated respectively, and the three are summed to obtain the total entropy generation rate of the system.

[0117] In this embodiment, this step decomposes the total irreversibility of the system into three main sources and performs accurate calculations using the refined model obtained in the previous embodiment.

[0118] Specifically, the process is as follows: when calculating the flow entropy production rate and the heat transfer entropy production rate, the topological structure of the heat transfer path network is used as the basis, and the irreversibility of the fluid dissipation and heat transfer processes along the discrete paths of the network are quantified respectively;

[0119] When calculating the pressure drop entropy generation rate, the memory-corrected state equation is used to calculate the throttling and friction losses, and the memory effect of historical operation is reflected in the entropy generation analysis.

[0120] Both the flow entropy generation rate and the heat transfer entropy generation rate are based on the topological structure of the heat transfer path network generated in Example 2. Specifically, the flow entropy generation rate is obtained by calculating the viscous dissipation caused by the velocity gradient along each discrete path; the heat transfer entropy generation rate is obtained by calculating the irreversible heat transfer loss caused by the temperature difference along each path. This path-based calculation method can more accurately locate and quantify the main sources of entropy generation. Specifically as follows:

[0121] In other words, based on the topological structure of the heat transfer path network, the flow entropy generation rate S is calculated. flow For each heat transfer path, the viscous dissipation rate is calculated based on the velocity gradient of the flow field around it, and then the flow entropy generation rate S of the path is obtained. flow_i At the same time, the velocity field interference between adjacent paths is considered and the coupling entropy production term is calculated. Finally, the path weight W is used path_i The entropy production of each path and the coupled entropy production are weightedly integrated to obtain the total flow entropy production rate S flow_total .

[0122] Similarly, based on the heat transfer path network, the heat transfer entropy production rate S is calculated heat For each path, according to its heat transfer q heat_i and the temperature T at both ends of the path hot_i , T cold_i , calculate the heat transfer entropy production S heat_i =q heat_i *(1 / T cold_i -1 / T hot_i . Heat transfer q heat_i The calculation of depends on the path heat transfer characteristics and time-varying equivalent heat transfer coefficient obtained in step S302. Finally, the dynamic weight W path_i (t) Weighted summation to obtain the total heat transfer entropy production rate S heat_total (t).

[0123] The pressure drop entropy production rate is calculated using the memory-corrected state equation generated in Example 3. Specifically, when calculating the entropy production due to throttling and friction loss, the fluid's density ρ, enthalpy change Δh and other thermodynamic parameters are all solved by the corrected equation that includes the memory effect. This allows the entropy production analysis to take into account the impact of historical operation, which is unprecedented. The three entropy production rate components mentioned above are summed to obtain the total entropy production rate of the system S system_total .

[0124] Whether calculating the entropy generation caused by valve throttling or pipe friction, accurate pressure drop values ​​and gas density are required. By using the memory-corrected state equation to calculate these values, the memory effect of historical operation can be reflected in the entropy generation analysis, resulting in a more accurate pressure drop entropy generation rate S than traditional methods. pressure_total .

[0125] Sum S total (t)=S flow_total +S heat_total (t)+S pressure_total By integrating the time over a complete cycle of inflation and deflation, the total entropy production of a single cycle can be obtained.

[0126] Step S402: Establish a functional relationship between the total entropy production rate of the system and the volume of the gas storage reservoir, and solve the volume that minimizes the functional relationship while satisfying preset engineering constraints, and use it as the final design volume.

[0127] In one embodiment, the implementation process is as follows:

[0128] Establish an entropy production rate-volume relationship: Select a series of volume values ​​within a range around the current volume V (e.g., V ± 20%). Repeat step S401 for each volume value to obtain the corresponding total entropy production rate of the system. Using these (V, S_total) data points, we can fit the functional relationship Stotal = f(V).

[0129] Constructing a multi-objective optimization function: Actual engineering design should not only consider thermodynamic efficiency, but also economic and performance. Therefore, construct a multi-objective optimization function: F(V)=w1*(S total (V) / S ref )+w2*(Cost(V) / Cost ref )+w3*(Performance_penalty(V)); where Cost(V) is the construction cost function related to the volume, and Performance_penalty(V) is the penalty function when certain performance indicators are violated.

[0130] Set constraints: including power constraints (P out ≥P design ), safety pressure constraint (p min ≤p≤p max ) and investment cost constraint (Cost(V)≤Budget).

[0131] Step S403: executing an iterative optimization process, using a coupled iterative solution algorithm with multi-variable synchronous updating.

[0132] Multivariable synchronous updating refers to an iterative strategy for solving strongly coupled problems. In one iteration, all interdependent variables are updated sequentially instead of optimizing a single variable independently.

[0133] For example, within an iterative loop, the four core coupling variables, namely the gas storage volume V, memory function M, heat transfer coefficient h, and entropy production rate S, are updated in a specific logical order rather than independently and in parallel.

[0134] Specifically, in a single iteration, the memory-corrected state equation, the time-varying equivalent heat transfer coefficient, and the total entropy production rate of the system are updated in sequence according to the current iteration's gas storage volume. Based on the updated total entropy production rate, the new gas storage volume for the next iteration is solved.

[0135] The single iteration cycle is repeated until the changes in the gas storage volume, the memory-corrected state equation, the time-varying equivalent transfer coefficient, and the total entropy production rate of the system all converge to their respective preset convergence criteria.

[0136] In one embodiment, the process is as follows:

[0137] Start with an initial volume V(0) (eg, a volume calculated by conventional methods).

[0138] A single iteration loop (from nth to n+1th) proceeds as follows:

[0139] Based on the current volume V(n) and historical trajectory, the memory function is recalculated to obtain the updated memory-corrected state equation (i.e., the updated M memory (n+1)), that is, the operation of updating the memory effect is completed.

[0140] Based on the current volume V(n) (affecting geometry) and the updated memory effect (affecting flow field), the heat transfer path network is reconstructed and the updated time-varying equivalent heat transfer coefficient h is calculated. eff (n+1).

[0141] Based on the updated M memory (n+1) and h eff (n+1), recalculate the total entropy production rate S of the system total (n+1).

[0142] Solve the optimization problem minF(V), where the entropy production term uses the S just calculated total (n+1), and obtain a new gas storage volume V(n+1) for the next iteration.

[0143] Repeat the single iteration loop until the changes of all core variables are less than their respective preset convergence criteria, such as satisfying |V(n+1)-V(n)|<εV, |Mmemory (n+1)-M memory (n)∣<εM,∣h eff (n+1)-h eff (n)∣<εh and∣S total (n+1)-S total (n)∣<εS.

[0144] The final output volume V(n+1) after iterative convergence is the final design volume V that minimizes the total entropy production rate of the system. final .

[0145] In other words, in a loop n, the algorithm updates other variables in turn according to the current volume V_n:

[0146] M n+1 =f(Vn, historical data);h n+1 =f(Vn,M n+1 );S n+1 =f(Vn,M n+1 , h n+1 ); Based on the updated S n+1 , using quasi-Newton methods such as BFGS, a new gas storage volume V_{n+1} is solved.

[0147] Repeat the above cycle until the changes of all core variables converge to their respective preset acceptance criteria. For example: relative change in volume |V n+1 -V n | / V n <0.1%; relative change in entropy production rate |S n+1 -S n | / S n <0.1%; when all criteria are met at the same time, the iteration terminates, and the volume V_{n+1} at this time is the final design volume.

[0148] This embodiment enables the four mutually coupled aspects of volume, gas state, heat transfer characteristics and system efficiency to reach the optimal state in a coordinated and synchronous manner, avoiding the problems of local optimality and global suboptimality caused by ignoring the coupling effect in traditional step-by-step design, and realizing true system-level global optimization.

[0149] Example 5: The process of constructing the fluid memory effect model involved in Example 2 is described in detail, especially the determination method and physical connotation of key parameters.

[0150] Step S501 : calibrate and interpret key parameters of the memory effect model.

[0151] The memory decay coefficient λ is a parameter used to characterize the rate at which historical cycle information decays over time. Its physical meaning is related to the thermal inertia of the gas storage surrounding rock, that is, the time scale required for the surrounding rock to forget a thermal shock.

[0152] The coupling weight coefficients γ1, γ2, and γ3 respectively represent the contribution weights of the three effects of molecular action memory, compression history memory, and temperature history memory to the total influence of the current gas state.

[0153] The value of the memory decay coefficient λ can be determined by fitting the historical temperature response data of the gas storage. Specifically, after a significant injection / production event, the temperature recovery curve at a certain point on the gas storage wall is monitored over time. This curve can generally be approximated as an exponential decay. By fitting this curve to the decay function e-λt, the value of λ can be determined. For example, for large salt rock gas storage, which has a long thermal relaxation time, the value of λ is typically between 0.1 and 0.3.

[0154] The coupling weight coefficients γ1, γ2, and γ3 can be calibrated using supervised learning methods. First, a set of input-output experimental data sets containing various working conditions is prepared, where the input is the historical state trajectory and the output is the measured gas storage pressure. Then, a loss function is defined, such as the mean square error between the predicted pressure and the measured pressure: L(γ1, γ2, γ3) = ∑ j=1 M (p predicted，j(γ1，γ2，γ3) -p actual，j ) 2 Finally, a nonlinear optimization algorithm (such as the Levenberg-Marquardt algorithm) is used to find the value of (γ1, γ2, γ3) that minimizes the loss function L. This process ensures that the value of the weight coefficient is data-driven rather than subjective.

[0155] The necessity of triple coupling is explained as follows: three memory functions are needed because they correspond to different physical mechanisms. molecular ) corresponds to the historical dependence of the non-ideal nature of gas under high pressure; the compressed historical memory (M compression ) corresponds to the mechanical behavior memory of the surrounding rock under cyclic load; temperature history memory (M thermal ) corresponds to the heat storage and release effects of the surrounding rock. The lack of any one of them will lead to deviations in the model's predictions for specific working conditions.

[0156] Step S502: Refine the calculation method of trajectory similarity.

[0157] In the second embodiment, the core of the dynamic time warping (DTW) distance is to calculate the cumulative distance matrix D. For two time series trajectories Ta (length m) and Tb (length n), the element D(i, j) of the matrix D is calculated using the following recursive formula: D(i, j) = d(Ta(i), Tb(j)) + min{D(i-1, j), D(i, j-1), D(i-1, j-1)}, where d(Ta(i), Tb(j)) is the Euclidean distance between points Ta(i) and Tb(j). The recursive calculation starts from D(1, 1) and ends at D(m, n), and the final value is the DTW distance between the two trajectories. It can effectively measure the similarity between trajectories that are similar in shape but may be scaled or translated along the time axis.

[0158] Step S503: verifying the long-term stability of the memory function model.

[0159] To ensure that the memory function does not diverge over long, multi-loop iterations, numerical methods can be used to assess its stability. For example, the recursive process of the memory function can be treated as a discrete dynamic system. By calculating its state evolution after multiple cycles, the maximum Lyapunov exponent of the system can be evaluated. A negative value of this exponent indicates that the system is stable and that historical perturbations converge over time. A positive value indicates that chaotic or divergent behavior may exist, and the model parameters or structure need to be re-examined.

[0160] Example 6: The process of establishing the path topology dynamic heat transfer model involved in Example 3, especially its judgment criteria and parameter selection basis, is explained.

[0161] Step S601: Select a gradient threshold and a path termination condition.

[0162] Gradient threshold T threshold , h threshold It is not based on empirical settings, but is determined through sensitivity analysis. Specifically, the final calculated full-cycle average equivalent heat transfer coefficient h -eff As the goal, analyze its effect on T threshold and h threshold Plot h -eff The optimal threshold value is the inflection point where the slope of the curve decreases from large to small and then flattens out. This point means that the contribution of gradients below this threshold to the overall heat transfer is negligible, while maintaining computational efficiency.

[0163] The path termination condition is further as follows: the gradient modulus | | Direction | | is less than a minimum value εg (such as machine precision). The coordinates of the path tracking point enter the predefined exit or wall absorption area. The actual path length L actual Exceed the characteristic size of the gas reservoir by a preset multiple (e.g. 3) to prevent infinite loops.

[0164] Step S602: Constructing a more refined path classification criterion, including:

[0165] Introducing path tortuosity τ path : Defined as the actual length of the path L actual The straight-line distance L between its two end points straight The ratio of τ path =L actual / L straight .

[0166] Introducing local curvature k local The curvature of each point on the path. New classification criteria, including:

[0167] Direct path: τ path <1.1. This type of path has direct energy transfer and low loss.

[0168] Reflection path: There is at least one point on the path with local curvature k local Exceeding a large reflection threshold k reflect , and the path as a whole τ path Smaller.

[0169] Flow path: τ path >1.5. This type of path usually occurs near obstacles or geometric dead corners, and the heat transfer efficiency is low.

[0170] Hybrid paths: Other paths that do not fall into the above categories. This classification method is based on quantitative indicators, eliminates ambiguity, and provides a basis for subsequent matching of different heat transfer models (such as the reflection path model that considers local enhanced heat transfer at the inflection point).

[0171] Step S603: Optimize the public weight index.

[0172] The exponents α, β, and γ in the path weight formula are similarly obtained using a parameter calibration method similar to that described in Example 5. Prepare a set of experimental or high-precision simulation results containing precisely measured local heat flux density data. With the goal of minimizing the error between the predicted and measured heat fluxes, construct a loss function and use a global optimization algorithm (such as particle swarm optimization or genetic algorithm) to search for the optimal exponent combination (α, β, γ).

[0173] Example 7: The calculation method of the total entropy production rate of the system involved in Example 4 is explained, especially how to perform the calculation based on a discrete path network and how to determine the multi-objective optimization weights.

[0174] Step S701 : providing a complete calculation method for flow and heat transfer entropy production based on a path network.

[0175] The entropy production rate is essentially the spatial integral of a field quantity. To combine it with a discrete path network model, it is necessary to transform the continuous integral operation into a summation operation over discrete paths. Specifically, the following steps are involved:

[0176] Calculate the flow entropy production rate: For each path P in the heat transfer path network athi , which consists of a series of discrete points {P1, P2, ..., P k}. At each point P j The velocity vector vj can be obtained from the CFD flow field data.

[0177] The components of the velocity gradient tensor, such as dydvx, can be approximated by the central difference scheme: dydvx≈2Δyvx(y+Δy)-vx(y-Δy), where the required velocity value is obtained by interpolating the flow field data around the path point Pj.

[0178] Substituting all the calculated velocity gradient components into the formula of the viscous dissipation function Φv, the viscous dissipation rate at point Pj is obtained.

[0179] Path P athi The total flow entropy generation rate S flow_i By numerically integrating the dissipation rate of all points on the path along the path length (e.g. trapezoidal rule), we can obtain: S flow_i =∫ Path_i (μ / T local )Φ v ds≈∑ j=1 k-1 (μ / T j )}Φ v d v,j δ j .

[0180] Finally, the total flow entropy production rate is obtained by weighted summing the contributions of all paths: S flow_total =∑ i W path_ i S flow_i .

[0181] Similarly, the total heat transfer entropy production rate can also be calculated by taking the weighted sum of the heat transfer entropy production rates of each path. athi Heat transfer entropy generation S heat_i It can be obtained by integrating the local heat transfer entropy production at each point on the path. The local heat transfer entropy production density is k(gradT) 2 / T 2 .

[0182] Step S702: Determine the weight of the multi-objective optimization function, specifically:

[0183] To avoid subjectivity of weights w1, w2, and w3, this embodiment uses the analytic hierarchy process (AHP) to determine them.

[0184] Invite multiple domain experts to conduct pairwise comparisons on the three objectives of minimizing entropy production (efficiency), minimizing costs (economy), and minimizing performance penalties (safety and reliability) to construct a judgment matrix. For example, if an expert believes that efficiency is slightly more important than economy, then fill the corresponding position with 3; otherwise, fill the corresponding position with 1 / 3.

[0185] Calculate the maximum eigenvalue λ of the judgment matrix max and its corresponding normalized eigenvector, which is the required weight vector (w1, w2, w3) T .

[0186] Calculate the consistency ratio CR=CI / RI. Where CI=(λ max -N) / (N-1), where N is the target number (here, 3) and RI is the preset average random consistency index. CR is usually required to be less than 0.1. If it is not met, the experts need to readjust the judgment matrix to ensure the logical consistency of the decision.

[0187] Step S703: designing a penalty function to handle engineering constraints.

[0188] For the constraints in the optimization problem, such as the safety pressure constraint p≤p max , an external penalty function can be designed. When p>p max When the penalty term P(V)=R×(pp max ) 2 , where R is a large positive number (penalty factor). This causes the objective function value corresponding to any solution that violates the constraint to increase dramatically, and thus be naturally excluded during the optimization process.

[0189] Example 8: The numerical implementation, convergence guarantee and robustness measures of the coupled iterative optimization process described in Example 4 are described.

[0190] Step S801: providing an initial value selection and divergence processing mechanism for the iterative algorithm.

[0191] The initial volume V(0) for the coupled iterations should not be chosen randomly. A preferred strategy is to first perform a quick calculation using a simplified conventional model (e.g., an ideal gas model that ignores memory effects and heat transfer evolution) to obtain a preliminary volume estimate. This value is then used as the initial value V(0) for the complex coupled iterations. This ensures that the starting point of the iteration is within the attracting domain of the optimal solution, accelerating convergence.

[0192] During the iteration process, a monitor is set to track the objective function value F(Vk) in real time. If F(Vk) appears for more than a preset number of times (for example, 3 times) k )>F(V k-1 ), then the iteration is judged to be divergent. At this time, the divergence processing mechanism is triggered: the current step size η k Forced halving, i.e. η k ←η k / 2.

[0193] Abandon the currently calculated divergent state X(k) and roll back the system state to the previous stable point X(k-1).

[0194] Restart iteration from X(k-1) using the halved step size. If rollback is triggered multiple times in a row, the user will be notified of convergence failure and advised to check the model parameters.

[0195] In step S802, a fixed step size may cause oscillation near the optimal point or slow convergence. An adaptive step size can dynamically adjust the step size according to the progress of the iteration.

[0196] In the optimization subproblem of solving the new volume, a line search method with the Armijo criterion is used to determine the step size η k At the kth iteration, starting from an initial step size (e.g. η0 = 1), we multiply it by a reduction factor (e.g. 0.5) until we find the first step size η that satisfies the following inequality: k :

[0197] F(V k -η k gradF(V k ))≤F(V k )-c·η k ∣∣gradF(V k )∣∣ 2 ; where c is usually 10 -4 This criterion ensures that the objective function can be sufficiently reduced at each iteration, thus theoretically ensuring the convergence of the algorithm to the local optimal solution.

[0198] Step S803: sensitivity analysis of the final volume.

[0199] After obtaining the final optimized volume V final Finally, a sensitivity analysis is performed to assess its robustness to the uncertainty of key input parameters. Key input parameters are identified, such as the thermal conductivity k of the rock mass in the geological data. rock , the power generation power P in the initial design requirements designFor each key parameter, perturb it within its possible value range (for example, ±10%). After each perturbation, rerun the entire optimization process to obtain a new optimal volume.

[0200] The final volume V is calculated by final The partial derivative (or coefficient of variation) with respect to the change of each parameter quantifies its sensitivity. For example, dV final / dk rock The larger the absolute value of , the higher the accuracy requirement for the rock mass thermal conductivity in the final design. This analysis result can guide the aspects that require key surveys and precise designs in engineering practice.

[0201] In some embodiments of the present application, the physical mechanism of the memory effect is further elaborated. Molecular interaction memory can be traced back to the Virial expansion of the real gas state equation in thermodynamic statistics. The drastic changes in historical pressure and temperature change the average distance and energy distribution between gas molecules, so that the Virial coefficients (such as the second Virial coefficient B(T)) that characterize the intermolecular attraction and repulsion show dependence on the historical path, rather than just a function of the current temperature. The compression and temperature history memory is related to the thermal-mechanical coupling behavior of the surrounding rock of the gas storage reservoir. As a huge heat capacity and elastic body, the rock mass has a hysteresis and relaxation effect in its response to the gas temperature and pressure. This macroscopic hysteresis effect is fed back to the gas state equation through the memory function.

[0202] It should be noted that the memory effect model described in this paper is primarily applicable to non-ideal gas environments with high pressures (e.g., greater than 5 MPa) and wide temperature ranges (e.g., cyclic temperature differences greater than 50 K). Under these conditions, intermolecular interactions and heat exchange effects with the surrounding rock are significant. Under low-pressure, constant-temperature conditions close to those of an ideal gas, the value of the memory function approaches zero, and the model automatically degenerates into a traditional equation of state.

[0203] Specifically, the calculation process of the Lyapunov exponent is as follows: First, the iterative process of the memory function is represented as a discrete mapping of the state vector Xk+1=F(Xk). Then, a small perturbation δX0 is imposed on the reference trajectory. The evolution of the perturbation is calculated through multiple iterations: δX k =JF(X k-1 )·JF(X k-2 )…JF(X0)·δX0, where JF is the Jacobian matrix of the mapping F. The maximum Lyapunov exponent L is defined as: L=limk→∞k1ln|||δX0||||||δX k In numerical calculations, this exponent can be stably solved by performing QR decomposition on the product of Jacobian matrices.

[0204] Furthermore, the reflection threshold k for path classification reflectThe characteristic hydraulic diameter D of the gas storage reservoir h When the local curvature radius of the path is less than D h When the flow state changes significantly like hitting the wall, it can be considered that reflection occurs. reflect Can be set to 2 / D h .

[0205] It should be noted that the maximum path length is set to three times the characteristic size. The theoretical basis for this is that, for diffusion or convection processes within a closed or semi-enclosed space, if a particle's trajectory is too long to reach the boundary, it usually means that it has been trapped in an eddy or local circulation zone. In heat transfer analysis, such overly circuitous paths contribute minimally to overall heat transport over long distances. Setting this upper limit is a computational trade-off, intended to prioritize the capture of the main paths that are decisive for global heat transfer within limited computing resources, avoiding endless and inefficient tracking in local dead zones.

[0206] For paths classified as mixed, segmentation or weighted averaging can be employed. For example, a mixed path can be decomposed into several subpaths with distinct characteristics (e.g., approximately straight segments and curved segments). The corresponding heat transfer correlation is then calculated for each segment and integrated along the path. Alternatively, the equivalent heat transfer coefficient can be estimated by taking a weighted average of the heat transfer correlations for the direct and bypass paths based on the overall path geometric parameters (e.g., mean curvature and tortuosity).

[0207] Optionally, the coupling entropy production between paths can be calculated by introducing an influence function. athj The mutual interference intensity of the velocity field or temperature field can be modeled as a function of the path spacing d ij The attenuation function f(dij). The coupling entropy production term can be approximately expressed as Scoupling ij =C·f(d ij )·∣∣gradvi-gradv j ∣∣ 2 , where C is the coupling coefficient. The total coupling entropy production is the sum of this term for all adjacent path pairs.

[0208] During the optimization process, an adaptive strategy can be used to select the penalty factor R. In the early stages of an iteration, R can be kept small to allow the optimization algorithm to explore a wider range of space. As the iterations progress, if the solution continues to wander outside the constraint bounds, the value of R is gradually increased (for example, by multiplying it by a factor greater than 1, such as 1.2, with each iteration) to impose a stronger penalty and force the solution to converge to the feasible region.

[0209] After the judgment matrices of multiple experts are integrated for group decision-making, if the overall consistency does not meet the requirements (CR>0.1), the following method can be used to handle it: first, identify the judgment element (or elements) that cause the greatest inconsistency; then, feedback the relevant information of the judgment element to each expert, and organize a round of back-to-back re-scoring or group discussion until the revised judgment matrix meets the consistency requirements.

[0210] Setting the divergence criterion to three consecutive times is an engineering heuristic that strikes a balance between sensitivity and stability. Setting it to just one time is overly sensitive and may misinterpret normal oscillations during the iteration process as divergence. Setting it too many times (e.g., 10 times) may cause the algorithm to stray too far from the optimal solution, making it difficult to recover. Three times is generally considered a reasonable observation window, sufficient to distinguish a persistent deteriorating trend from temporary random fluctuations.

[0211] If the coupled iterations ultimately fail to converge, the system should provide a fail-safe design solution. This solution is based on the most conservative physical model calculations. For example, it completely ignores memory effects (using the ideal gas equation of state) and uses the most conservative heat transfer coefficient value (i.e., the most intense heat transfer) in existing design specifications. Although this solution may be less economical, it ensures absolute design safety and can be submitted as a final alternative.

[0212] Example 9 provides a comprehensive design case from beginning to end that integrates all of the above technologies (including the content in the supplementary explanation) to fully demonstrate the complete application process, theoretical support and technical advantages of the present invention in engineering practice. It should be noted that the content that has been described in detail in the above examples is omitted here.

[0213] Design of the core salt rock gas storage for a planned 1000MW / 8000MWh compressed air energy storage power station in western China.

[0214] Step S901: Initial parameter determination and data preprocessing

[0215] Receiving design requirements: Rated power generation P design =1000MW, energy storage duration t storage =8h, cycle efficiency target>60%.

[0216] Geological exploration data revealed a roof depth of 1,000 meters and a thickness of 150 meters within the target salt layer. Kriging interpolation was performed on the borehole sampling data to generate a three-dimensional distribution field of rock mass mechanical parameters (elastic modulus, Poisson's ratio) and thermophysical parameters (thermal conductivity, specific heat capacity) within the gas storage area. Data integrity testing confirmed that the sampling density in key areas met design requirements.

[0217] Combined with geological conditions and safety regulations, the pressure range is determined to be p min =6MPa, p max =12MPa; considering the initial ground temperature of the surrounding rock, the initial temperature T0 is set to 320K. The gas storage reservoir is initially designed as a vertical cylinder with a height-to-diameter ratio of 3.

[0218] Step S902: Preliminary calculation of the volume of the fused memory effect

[0219] The high-pressure air and surrounding salt rock are considered as a tightly coupled non-ideal system. The van der Waals effect of the gas and the thermal hysteresis of the surrounding rock together constitute the memory of the system.

[0220] Using historical data of similar salt rock gas storage, the memory function coupling weights are calibrated to γ1=0.25, γ2=0.4, and γ3=0.35 using the least squares method described in Example 5.

[0221] Collect historical data and build a standardized status trajectory database.

[0222] According to the method of embodiment 5, the comprehensive memory weight vector and the total memory function M are calculated. total .

[0223] Establish memory correction state equation: p×V=nZ eff RT, where the equivalent compression factor Z eff =1+α memory M total .

[0224] By solving the revised energy balance equation, the initial volume V1=2.8×105m3 considering the memory effect is obtained.

[0225] Step S903: Analyze the dynamic heat transfer process based on the path topology

[0226] This step transforms traditional lumped parameter-based heat transfer analysis into a structure- and flow-based, path-topology analysis with a clear, physical image. Its theoretical basis for heat transfer is that, in complex internal convection, the total heat transfer can be viewed as the superposition of energy transport paths driven by a series of coherent structures of varying scales and shapes.

[0227] The gas storage reservoir with a volume of V1 is subjected to non-uniform meshing, and the mesh quality (such as aspect ratio and orthogonality) is checked using the method described in Example 6 to ensure high-quality mesh.

[0228] Based on the temperature field and velocity field obtained by the preliminary CFD simulation, the heat transfer path network is generated and classified using the quantitative criteria in Example 6 (such as tortuosity and local curvature).

[0229] For each path, the heat transfer characteristics are calculated by selecting the best matching formula from a preset heat transfer correlation library containing heat transfer formulas under different Reynolds numbers and geometric conditions.

[0230] The dynamic path weight is calculated by the method of Example 3 and finally integrated into the time-varying equivalent heat transfer coefficient h eff (t). Analysis shows that in the initial stage of gas injection, h eff (t) can be as high as 20W / (m 2 ·K), and at the end of gas storage it drops to 3W / (m 2 ·K).

[0231] h eff (t) Substitute into the heat transfer evolution differential equation and calculate the heat transfer corrected volume V2 = 3.1×10 5 m 3 .

[0232] Step S904: System optimization based on minimization of entropy production

[0233] For an optimal engineering system, the energy loss (i.e., total entropy generation) caused by irreversibility should be minimized. To this end, the calculation process is as follows:

[0234] A volume parameter sweep is performed with V2 as the center. For each volume, the memory effect model and the path heat transfer model are called simultaneously.

[0235] The path network-based numerical integration method described in Example 7 is used to calculate the flow entropy generation, heat transfer entropy generation, and pressure drop entropy generation respectively.

[0236] The AHP method was used to determine the weights of each optimization objective. After scoring and consistency testing by the expert group, the weight vectors were obtained as follows: efficiency (0.55), economy (0.30), and safety (0.15).

[0237] A multi-objective optimization function including a penalty function was constructed to find the minimum entropy production rate. The Pareto frontier solution set obtained during the optimization process showed that the system's overall performance reached its optimal value when the volume was V3 = 3.0×105m3.

[0238] Step S905: Final volume determination and verification

[0239] Physical processes such as gas non-ideality, wall heat transfer, and flow dissipation are coupled and influence each other. Through a synchronous iterative framework, a self-consistent solution for the entire coupled system is sought, and its mathematical convergence can be proven using fixed point theories such as the compression image principle. The details are as follows:

[0240] By conducting numerical perturbation experiments near the optimization point, the normalized sensitivity coefficients between the core variables are calculated and used as coupling strength parameters.

[0241] With V3 as the initial value, the multivariable synchronous update algorithm with adaptive step size and divergence processing mechanism described in the eighth embodiment is started.

[0242] After 15 iterations, the relative changes in volume, memory function, heat transfer coefficient and entropy production rate are all less than 10-4, meeting the convergence criterion. The final design volume V is determined. final =3.02×10 5 m 3 .

[0243] V final The data was input into independent third-party geomechanics and thermodynamics simulation software. The results showed that during the simulated 10-year operation period, the gas storage's circulating pressure and temperature remained within the design envelope, and the maximum deformation of the surrounding rock was less than the safety threshold, verifying the effectiveness of the design.

[0244] Compared with the design scheme using the industry's common simplified method (ideal gas + fixed heat transfer coefficient), the volume obtained by the present invention is reduced by about 13%, directly saving tens of millions of yuan in cavity manufacturing costs, while the predicted value of the energy storage cycle efficiency is increased by 2.5 percentage points.

[0245] By constructing and applying a fluid memory effect model, the accuracy of the prediction of the true thermodynamic state of gas in the gas storage can be significantly improved, especially in complex operating scenarios with long cycles and multiple operating conditions. Specifically, by constructing a standardized state trajectory database containing multiple complete charging and discharging cycles, a rich historical information foundation is provided for the model; then, by constructing a triple memory function containing molecular interactions, compression history, and temperature history, and combining it with dynamically assigned comprehensive memory weights, the model can quantify and distinguish the cumulative impact of different historical events on the current state; finally, this total memory function is introduced as a correction factor into the benchmark gas state equation, so that the state equation is no longer an isolated instantaneous relationship, but contains historical dependence in the time dimension. For the specific scenario of large underground gas storage, the huge heat capacity characteristics of the surrounding rock make the thermal memory effect particularly prominent. The traditional method ignores this effect and the energy estimation error can reach 5%-10%. By precisely accounting for this memory effect, this method can more accurately predict key parameters such as gas density and enthalpy, thereby providing a solid foundation for preventing unexpected phase changes (such as ice blockage) during rapid energy release, accurately calculating energy storage efficiency, and subsequently accurately calculating the pressure drop entropy production rate. This solves the problems of insufficient control accuracy and potential safety risks caused by the memoryless assumption in existing technologies.

[0246] By establishing a path topology dynamic heat transfer model, the existing problem of crude heat transfer coefficient values ​​that fail to reflect the spatial heterogeneity within the gas storage reservoir can be resolved, enabling accurate prediction of localized overheating or undercooling areas. The approach primarily deconstructs the complex overall heat transfer problem into multiple parallel, physically distinct subprocesses. A heat transfer path network consisting of multiple discrete paths is generated based on the reservoir geometry and physical gradient field, providing a structural framework for energy transfer. The geometric and thermodynamic properties of each path in the network are independently calculated, forming a set of path-specific heat transfer characteristic data. In response to real-time changes in the internal flow field state, a dynamically changing weight vector is assigned to each path, and ultimately, through weighted integration, the time-varying equivalent heat transfer coefficient is fused. In large-scale, megawatt-class gas storage reservoirs, due to their enormous spatial scale and complex geometry, dead zones and convection-dominated channels inevitably exist. These are condensation-dominated paths that cannot be captured by traditional lumped parameter models. By dynamically changing the path weights, the heat transfer intensity in these critical areas can be quantitatively identified in real time, providing unprecedentedly precise guidance for avoiding catastrophic local condensation risks. This allows designers to move away from the "global over-drying strategy" of over-drying 100% of the gas to address the 1% risk zone. While ensuring safety, it significantly reduces unnecessary energy consumption, offering significant economic and engineering value.

[0247] By deeply coupling entropy generation calculations with two physical models—fluid memory effects and path topology heat transfer—this method establishes a more accurate and reliable system for thermodynamic irreversibility assessment based on real physical processes rather than simplified assumptions. Specifically, the calculations of flow and heat transfer entropy generation rates are based on the topological structure of the heat transfer network, quantifying the macroscopic entropy generation down to each discrete energy transfer path, thereby identifying the primary energy dissipation regions and causes. When calculating pressure drop entropy generation, a memory-corrected equation of state is used to account for throttling and friction losses. This means that the impact of historical memory effects on current energy dissipation is accurately accounted for for the first time. For compressed air energy storage systems, the efficiency bottleneck lies precisely in the irreversible losses during compression, heat transfer, and throttling. By more accurately quantifying these losses, subsequent optimization (i.e., finding the volume that minimizes the total entropy generation rate) is based on a more realistic foundation. It is no longer a rough, experience-based optimization, but a rigorous optimization process guided by the second law of thermodynamics and supported by precise physical models. The design volume ultimately found can fundamentally reduce the ineffective dissipation of energy throughout the life cycle of the gas storage facility and improve the overall operating efficiency of the energy storage system.

[0248] By employing a coupled iterative solution algorithm with simultaneous updates of multiple variables, this method can obtain a globally optimal and physically self-consistent gas storage design volume, effectively avoiding the inherent local optima and global suboptimal drawbacks of traditional sequential design methods. Specifically, within a single iteration, rather than optimizing each variable independently, the associated memory-corrected equation of state, time-varying equivalent heat transfer coefficient, and system total entropy generation rate are sequentially updated based on the current iteration's gas storage volume. Based on the updated system total entropy generation rate, the new gas storage volume for the next iteration is then solved. This forms a closed-loop feedback and correction mechanism, acknowledging and addressing the strong coupling between various physical fields (flow, temperature, and pressure) and geometric parameters (volume) in gas storage design. This coupling effect is particularly pronounced in the design of large-scale energy storage power plants: small changes in volume can lead to significant changes in flow velocity, which in turn affects heat transfer and flow dissipation; and gas history, in turn, affects its density, altering the pressure distribution and flow pattern. The coupled iterative algorithm ensures that the final output design volume is one of the few, or even the only, solutions that can simultaneously satisfy the three interrelated and complex constraints of fluid memory effect, dynamic heat transfer characteristics, and minimization of entropy production rate, thus ensuring the overall coordination, robustness and ultimate optimal performance of the design scheme.

Claims

1. A method for designing the volume of a large-scale compressed air energy storage power station gas storage reservoir, characterized in that: include: Receive design requirement data and determine initial design parameters accordingly; Based on historical operating data and initial design parameters, a fluid memory effect model is constructed to generate a memory-corrected state equation for subsequent entropy production rate calculation; Based on the initial design parameters, a path topology dynamic heat transfer model is established to obtain the time-varying equivalent heat transfer coefficient for subsequent entropy production rate calculation; The total entropy production rate of the system, which represents the thermodynamic irreversibility of the system, is calculated by coupling the memory-corrected state equation with the time-varying equivalent heat transfer coefficient. Iterative optimization is performed to search for the gas storage volume that minimizes the total entropy production rate of the system, and this volume is determined as the final design volume; Among them, the path topology dynamic heat transfer model is established, including: based on the geometric shape of the gas storage defined by the initial design parameters, the physical gradient field of the wall is analyzed, and a heat transfer path network composed of multiple discrete paths is generated accordingly; for each path in the heat transfer path network, its geometric and thermodynamic properties are independently calculated to obtain a set of path-specific heat transfer characteristic data; in response to the real-time changes in the flow field state inside the gas storage, a dynamically changing weight is assigned to each path to form a dynamic path weight vector; combined with the dynamic path weight vector, the path-specific heat transfer characteristic data is weighted and integrated to form a time-varying equivalent heat transfer coefficient; Generating a heat transfer path network includes: selecting grid points in a physical gradient field where both temperature gradients and height gradients exceed respective preset thresholds, and aggregating these points into a set of path starting points; using each point in the set of path starting points as a source point, executing an iterative tracking algorithm along the composite direction of the temperature gradient and height gradient to generate a coordinate sequence of the path point by point until a preset termination condition is met; the coordinate sequences of the multiple paths together constitute the heat transfer path network; Each path is assigned a dynamically changing weight to form a dynamic path weight vector. This includes: evaluating the local Reynolds number and pressure difference at both ends of each path based on the internal flow field state, and calculating the path activation probability based on this to indicate whether the path is activated; extracting the effective heat transfer cross-sectional area and path thermal resistance of each path from the path-specific heat transfer characteristic data, and calculating its area weight and thermal resistance weight accordingly; fusing the path activation probability, area weight, and thermal resistance weight to generate a comprehensive path weight and normalizing it to form a dynamic path weight vector; Among them, obtaining a set of path-specific heat transfer characteristic data includes: analyzing the geometric characteristics of each path in the heat transfer path network, and classifying each path into one of the predefined path types based on preset geometric classification criteria; for each path, analyzing and adopting the heat transfer correlation formula that matches the predefined path type of the path from the heat transfer correlation formula, and calculating its heat transfer properties for the path in combination with the internal flow field state; and combining the heat transfer properties of all paths into path-specific heat transfer characteristic data.

2. The method according to claim 1, characterized in that Construct a fluid memory effect model and generate a memory-corrected equation of state, including: Process historical operating data and construct a standardized state trajectory database containing multiple complete charging and discharging cycles; Based on the standardized state trajectory database, the cumulative impact of historical cycles on the current gas state is calculated and constructed into a total memory function; The total memory function is used as a correction factor and introduced into the reference gas state equation to form the memory-corrected state equation.

3. The method according to claim 2, characterized in that Construct a total memory function, including: Based on the standardized state trajectory database, a molecular interaction memory kernel function is constructed to represent the intermolecular interaction memory, a compressed history memory function is constructed to represent the compressed history memory, and a temperature history memory function is constructed to represent the temperature history memory. The total memory function is formed by nonlinearly weighted coupling of the molecular interaction memory kernel function, the compression history memory function and the temperature history memory function.

4. The method according to claim 3, characterized in that Before constructing the total memory function, it also includes generating a comprehensive memory weight vector that assigns weights to each historical cycle. Specifically: Based on the time interval between each historical cycle and the current moment, a time weight that decreases as the time distance increases is calculated for each historical cycle through a time decay function; Calculate the geometric similarity between the state trajectory of each historical cycle and the current state trajectory, and assign trajectory similarity weights to each historical cycle based on the geometric similarity; The time weight and trajectory similarity weight are weighted and fused to generate a comprehensive weight for each historical cycle. The comprehensive weights of all historical cycles together constitute a comprehensive memory weight vector.

5. The method according to claim 1, wherein Calculate the total entropy production rate of the system and search for the final design volume, including: Based on the memory-corrected state equation and the time-varying equivalent heat transfer coefficient, the flow entropy generation rate, heat transfer entropy generation rate, and pressure drop entropy generation rate caused by fluid viscous dissipation, wall heat transfer, and inlet and outlet pressure drop are calculated respectively, and the total entropy generation rate of the system is obtained by summing the three. A functional relationship between the total entropy production rate of the system and the volume of the gas storage reservoir is established, and the volume that minimizes the functional relationship while satisfying the preset engineering constraints is solved as the final design volume.

6. The method according to claim 5, characterized in that The calculation of the flow entropy production rate and the heat transfer entropy production rate are based on the topological structure of the heat transfer path network, and the irreversibility of the fluid dissipation and heat transfer along the discrete paths of the network is quantified respectively. When calculating the pressure drop entropy generation rate, the memory-corrected state equation is used to calculate the throttling and friction losses, and the memory effect of historical operation is reflected in the entropy generation analysis.

7. The method according to claim 5, characterized in that The iterative optimization process uses a coupled iterative solution algorithm with synchronous updates of multiple variables, specifically: In a single iteration, the memory-corrected state equation, the time-varying equivalent heat transfer coefficient, and the total entropy production rate of the system are updated in sequence according to the gas storage volume of the current iteration. Based on the updated total entropy production rate of the system, the new gas storage volume for the next iteration is solved. The single iteration cycle is repeated until the changes in the gas storage volume, the memory-corrected state equation, the time-varying equivalent transfer coefficient, and the total entropy production rate of the system all converge to their respective preset convergence criteria.

Images