A large-diameter artificial compressed air energy storage chamber optimization design method and system

Abstract

The application relates to the technical field of underground geotechnical engineering, and discloses a large-diameter artificial compressed air energy storage chamber optimization design method and system. The method comprises the following steps: constructing a three-dimensional geometric model comprising surrounding rock, a concrete lining and a steel sealing layer, and dividing time sequence grid layers according to layered pouring height; sequentially activating the time sequence grid layers, calculating thermal shrinkage displacement, and performing normal accumulation on a contact interface to generate an initial emptying defect distribution; applying an alternating load to the contact interface to calculate fatigue damage, and updating the normal contact stiffness distribution; taking a required compensation stiffness distribution as a design variable, taking the stress variance of the steel sealing layer as a target to construct an optimization function, and solving the required compensation stiffness distribution; and discretizing the required compensation stiffness distribution to output grouting control parameters and mesh screen layout area parameters. The application can dynamically quantify interface emptying and stiffness degradation, and convert theoretical stiffness requirements into on-site construction guidance indicators, thereby guaranteeing the long-term operation safety of the chamber structure.

Description

Technical Field

[0001] This invention relates to the field of underground geotechnical engineering technology, specifically to an optimized design method and system for large-diameter artificial compressed air energy storage chambers. Background Technology

[0002] Large-diameter artificial compressed air energy storage chambers, as a type of efficient large-scale energy storage infrastructure, typically employ a composite load-bearing and sealing structure consisting of surrounding rock, concrete lining, and an internal steel sealing layer to cope with frequent filling and releasing of high-pressure gases. However, in actual engineering construction and long-term service, the interfacial stress state of this composite structure is extremely complex, and existing engineering design and evaluation methods still have significant limitations.

[0003] During the construction phase of the tunnel, the layered pouring of large-volume concrete generates a heat of hydration effect. Existing conventional numerical simulation and design methods typically employ static boundary condition settings, making it difficult to dynamically adjust the thermodynamic parameters of the upper concrete layer by incorporating the actual construction sequence and the temperature evolution of the already poured lower layer. This static and singular calculation model cannot accurately reproduce the true thermal shrinkage behavior of concrete, resulting in an inability to accurately predict the initial void defects caused by thermal shrinkage at the interface between the surrounding rock and the concrete lining. Furthermore, the geometric gaps at this interface are currently difficult to quantify effectively.

[0004] Once the energy storage chamber is operational, frequent inflation and deflation operations will subject the chamber structure to severe alternating loads. Initial interface voids left over from the initial construction will, under the continuous action of alternating loads, induce fatigue damage and propagation at the defect edges, resulting in severe degradation of local support stiffness. Current structural assessment systems generally lack the ability to quantitatively assess this dynamic degradation process of interface contact stiffness. Due to uneven distribution of support stiffness, the internal steel sealing layer, directly bearing the gas pressure, will experience localized stress concentration, increasing the risk of fatigue yielding and tearing of the sealing layer.

[0005] To eliminate the hidden dangers caused by local voids, grouting reinforcement is usually required in engineering projects. However, existing anti-void support parameter designs rely heavily on macroscopic engineering experience and lack a reverse optimization mechanism oriented towards structural stress safety. Existing technologies cannot derive the required interface compensation stiffness from the perspective of minimizing the stress variance of the steel sealing layer, nor can they scientifically map the theoretical compensation stiffness distribution into specific physical construction parameters such as discrete grouting hole locations, foundation grouting volume, target grouting pressure, and dense mesh layout area. This severe disconnect between theoretical mechanical analysis and precise on-site support operations leads to the blindness of grouting reinforcement, failing to form an engineering closed loop from defect prediction to anti-void support, and ultimately making it difficult to guarantee the safe and stable operation of large-diameter compressed air energy storage chambers over a long lifespan. Summary of the Invention

[0006] To address the shortcomings of existing technologies, this invention provides an optimized design method and system for large-diameter artificial compressed air energy storage chambers. This method solves the problems of existing design methods being unable to dynamically quantify the voids and stiffness degradation at the rock-lining interface caused by thermal shrinkage during construction and alternating loads during operation in large-diameter compressed air energy storage chambers. Furthermore, it lacks a mechanism for physical stiffness compensation and reverse optimization of grouting reinforcement parameters to address void defects, which leads to uneven stress on the steel sealing layer and difficulty in ensuring operational safety.

[0007] To achieve the above objectives, the present invention provides the following technical solution: The first aspect of this invention provides an optimized design method for a large-diameter artificial compressed air energy storage chamber, comprising the following steps: A three-dimensional geometric model of the chamber, including the surrounding rock, the concrete lining area and the steel sealing layer, is constructed and divided into finite element meshes. The concrete lining area is divided into multiple temporal mesh layers according to the layered pouring height. The time-series grid layers are activated sequentially according to the construction timeline. The thermal shrinkage displacement of each time-series grid layer is calculated and accumulated normally at the contact interface between the surrounding rock and the concrete lining area to generate an initial void defect distribution. An alternating load is applied to the contact interface, fatigue damage at the edge of the initial void defect distribution is calculated, and the normal contact stiffness distribution of the contact interface is updated based on the fatigue damage. Extract the stress field of the steel sealing layer under the alternating load, use the required stiffness distribution of the contact interface as the design variable, construct an optimization function with the goal of minimizing the stress variance of the steel sealing layer, and solve for the required stiffness distribution. The stiffness distribution requiring compensation is discretized, and the grouting control parameters and the parameters of the dense mesh layout area are output.

[0008] The specific innovative principles and implementation details of the above-mentioned optimized design method for large-diameter artificial compressed air energy storage chambers are explained below: To address the computational resource consumption problem caused by the large scale of large-diameter chamber models, this invention establishes a global and local collaborative grid processing mechanism.

[0009] By constructing a global three-dimensional geometric model containing the surrounding rock, the concrete lining area and the steel sealing layer and dividing it into finite element meshes, the stress gradient distribution under static analysis is extracted. The region where the stress gradient distribution exceeds a preset threshold is taken as a local sub-model, and the displacement field at the intersection node with the local sub-model is mapped to the intercept boundary of the local sub-model as the displacement boundary condition. Subsequently, based on the layered pouring height, the concrete lining area in the local sub-model is divided into multiple temporal mesh layers. This mechanism confines high-precision damage and thermal shrinkage calculations to high stress gradient regions, reducing the overall computational dimensionality while ensuring the accuracy of mechanical response calculations.

[0010] During the defect generation stage, this invention reconstructs the thermodynamic boundary conditions by combining the actual construction timeline. Specifically, it activates the current time-series grid layer to be calculated and obtains the temperature rise time series data of the activated time-series grid layers below it. Based on the temperature rise time series data, the specific heat capacity and heat generation rate coefficient of the current time series grid layer to be calculated are updated using a Bayesian inference algorithm. This step utilizes the measured or a priori temperature evolution law of the bottom layer of poured concrete to dynamically correct the constitutive thermodynamic parameters of the upper structure, and then calculates the temperature field and thermal shrinkage strain to generate the thermal shrinkage displacement.

[0011] The thermal shrinkage displacement is then accumulated normally on the contact interface. When the absolute value of the accumulation is greater than the void tolerance threshold, a geometric gap is generated. The geometric gaps across the entire domain constitute the initial void defect distribution.

[0012] This process achieves a quantitative mapping from the evolution of hydration thermophysics during construction to the geometric gaps in the structure.

[0013] During the operational damage assessment phase, this invention transforms macroscopic structural fatigue into the degradation of contact surface properties. By establishing a cohesive model at the contact interface and applying the alternating load, the cumulative fatigue damage variable of the interface after a preset number of cycles is extracted as fatigue damage. By utilizing the undamaged state proportion coefficient, the initial normal contact stiffness of the contact interface is transformed into a degraded normal contact stiffness matrix, which serves as the updated normal contact stiffness distribution. This mechanism accurately reflects the deteriorating effect of the gas chamber inflation / deflation cyclic load on the expansion of the vacancy edge.

[0014] To eliminate stress concentration in the steel sealing layer caused by local voids and stiffness degradation, this invention constructs a reverse optimization system guided by stiffness compensation. The adjoint variable method is used to derive the partial derivatives of the equivalent stress of the steel sealing layer with respect to the required stiffness compensation at each node of the contact interface, and a sensitivity Jacobian matrix is ​​constructed to obtain the nodal stress sensitivity.

[0015] To avoid singularities in the numerical solution process and accelerate convergence, the stress sensitivity of the nodes is sorted in descending order and the percentage is accumulated. The nodes in the latter column are filtered out according to the cutoff rate threshold, and the stiffness of the nodes in the latter column that need to be compensated is constrained to zero in subsequent iterations.

[0016] Subsequently, within the reduced-dimensional sensitive subspace, a sequential quadratic programming algorithm was used to solve the optimization function, yielding the required stiffness distribution with the goal of minimizing the stress variance of the steel sealing layer. Thus, the optimal support stiffness required to maintain uniform stress on the steel sealing layer was determined through theoretical numerical analysis.

[0017] Finally, this invention establishes a physical conversion model from theoretical stiffness to actual engineering construction parameters. By employing a spatial clustering algorithm, the non-zero nodes in the stiffness distribution requiring compensation are divided into multiple discrete grouting target areas, and the centroid coordinates of each discrete grouting target area are extracted as the target borehole coordinates. Further extract the volume of the initial void defect at the corresponding position of each discrete grouting target area as the basic grouting volume, establish a mapping model between the target grouting pressure and stiffness increment to solve the target grouting pressure, and generate grouting control parameters to guide on-site construction by combining the target hole position coordinates and the basic grouting volume.

[0018] In addition, the continuous void area of ​​each discrete grouting target area is calculated. If the continuous void area is greater than the critical safety area threshold, the coordinates of the corresponding boundary nodes are extracted to generate the parameters of the dense mesh layout area, thus realizing the engineering logic closed loop from defect prediction, stiffness solution to anti-void support parameter output.

[0019] A second aspect of the present invention provides an optimized design system for a large-diameter artificial compressed air energy storage chamber, which executes the optimized design method for a large-diameter artificial compressed air energy storage chamber provided in the first aspect of the present invention, including: The modeling and discretization module is used to construct a three-dimensional geometric model of the chamber, which includes the surrounding rock, the concrete lining area and the steel sealing layer, and to divide it into finite element meshes. The concrete lining area is divided into multiple temporal mesh layers according to the layered pouring height. The defect generation module is used to activate the time-series grid layers sequentially according to the construction time axis, calculate the thermal shrinkage displacement of each time-series grid layer, and perform normal accumulation on the contact interface between the surrounding rock and the concrete lining area to generate an initial void defect distribution. The stiffness update module is used to apply alternating loads to the contact interface, calculate the fatigue damage at the edge of the initial void defect distribution, and update the normal contact stiffness distribution of the contact interface based on the fatigue damage. The optimization solution module is used to extract the stress field of the steel sealing layer under the alternating load, construct an optimization function with the required stiffness distribution of the contact interface as the design variable and the goal of minimizing the stress variance of the steel sealing layer, and solve for the required stiffness distribution. The parameter output module is used to discretize the stiffness distribution that needs to be compensated and output grouting control parameters and dense mesh layout area parameters.

[0020] In one specific implementation, the defect generation module is specifically used to: activate the current time-series grid layer to be calculated, and obtain the temperature rise time series data of the activated time-series grid layer below the time-series grid layer; Based on the temperature rise time series data, the specific heat capacity and heat generation rate coefficient of the current time series grid layer to be calculated are updated using a Bayesian inference algorithm, and the temperature field and thermal contraction strain are calculated accordingly to generate the thermal contraction displacement. The thermal contraction displacement is accumulated normally on the contact interface, and a geometric gap is generated when the absolute value of the accumulation is greater than the void tolerance threshold. The initial void defect distribution is formed by the geometric gaps in the entire domain.

[0021] In one specific implementation, the optimization solution module is specifically used to: derive the partial derivative of the equivalent stress of the steel sealing layer with respect to the stiffness to be compensated at each node of the contact interface using the adjoint variable method, and construct a sensitivity Jacobian matrix to obtain the nodal stress sensitivity. The stress sensitivity of the nodes is sorted in descending order and the percentage is accumulated. The nodes in the latter column are filtered out according to the cutoff rate threshold, and the stiffness of the nodes in the latter column that needs to be compensated is constrained to zero in subsequent iterations. Within the dimensionality-reduced sensitive subspace, the optimization function is solved using a sequential quadratic programming algorithm to obtain the stiffness distribution that needs compensation.

[0022] In one specific implementation, the parameter output module is specifically used to: use a spatial clustering algorithm to divide the non-zero nodes in the stiffness distribution to be compensated into multiple discrete grouting target areas, and extract the centroid coordinates of each discrete grouting target area as the target hole position coordinates; The volume of the initial void defect at the corresponding position of each discrete grouting target area is extracted as the basic grouting volume. A mapping model between the target grouting pressure and the stiffness increment is established to solve for the target grouting pressure. The grouting control parameters are generated by combining the target hole position coordinates and the basic grouting volume. Calculate the continuous void area of ​​each discrete grouting target area. If the continuous void area is greater than the critical safety area threshold, extract the coordinates of the corresponding boundary nodes to generate the parameters of the dense mesh layout area.

[0023] This invention provides an optimized design method and system for a large-diameter artificial compressed air energy storage chamber. It offers the following advantages: 1. This invention constructs a three-dimensional geometric model of the chamber that coordinates global and local factors, and uses a Bayesian inference algorithm to dynamically correct thermodynamic parameters by combining the bottom temperature rise time series data. Then, the thermal contraction displacement is accumulated normally on the contact interface to generate the initial void defect distribution. This allows for accurate quantification of the mapping relationship from the evolution of hydration thermophysics during construction to void defects at the interface between the surrounding rock and the lining, while reducing the dimensionality of finite element calculation.

[0024] 2. This invention establishes a cohesive model at the contact interface to extract fatigue damage under alternating loads to update the normal contact stiffness distribution. By constructing a sensitivity Jacobian matrix, it inversely solves the stiffness distribution that minimizes the stress variance of the steel sealing layer in the dimension-reduced sensitive subspace. Therefore, it can dynamically evaluate the contact surface support degradation effect caused by the gas filling and defilling cycle, thereby eliminating the hidden danger of stress concentration in the steel sealing layer caused by local voids.

[0025] 3. This invention uses a spatial clustering algorithm to transform the stiffness distribution to be compensated into the target hole coordinates of the discrete grouting target area, and combines the initial void defect volume and stiffness increment mapping model to solve the grouting control parameters including the target grouting pressure and the dense mesh layout area. This allows the theoretical virtual stiffness requirements to be directly transformed into physical operation indicators that guide the refined construction on site, thereby establishing an engineering closed-loop mechanism from defect prediction to void prevention support. Attached Figure Description

[0026] Figure 1 Flowchart of the optimization design method provided in the embodiments of the present invention Figure 2 Flowchart of multi-scale modeling and temporal mesh reconstruction provided for embodiments of the present invention Figure 3 A flowchart of the defect dynamic inference driven by measured data provided in the embodiments of the present invention. Figure 4 A flowchart of cyclic service degradation evolution based on a cohesion model provided in this embodiment of the invention. Figure 5 A flowchart of the dimensionality reduction inverse optimization solution based on stress sensitivity provided in this embodiment of the invention. Figure 6 Flowchart for Discretization Mapping and Engineering Process Parameter Generation Provided in Embodiments of the Invention Figure 7 This is a diagram of the system architecture for optimizing and compensating for service degradation of the sealing structure of an artificial chamber, provided in an embodiment of the present invention. Detailed Implementation

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

[0028] See attached document Figure 1The optimized design method for large-diameter artificial compressed air energy storage chambers provided in this embodiment relies on a computer-aided engineering simulation platform and runs on a computer system equipped with a multi-core processor, memory, and data communication bus. To ensure the consistency of multi-source data, the computer system performs condition alignment processing based on a unified timestamp when receiving geological survey parameters, structural dimensions, and measured sensor data from the engineering site through the data interface. The processor executes the program code stored in the memory to achieve numerical simulation and inverse parameter optimization of the large-span, multi-layered coupled structure.

[0029] For physical scenarios involving large spans, deep burial, and complex multi-layered media coupling, conventional single-layer analytical solutions are insufficient to reflect the true nonlinear contact state. Based on the general finite element simulation principle, this embodiment constructs a three-dimensional geometric model of the chamber, including the surrounding rock, concrete lining area, and steel sealing layer, and divides it into finite element meshes. The concrete lining area in the three-dimensional geometric model is further divided into multiple temporal mesh layers according to the layered pouring height. This spatial modeling and temporal slicing discretization provides the geometric and mesh foundation for subsequently tracking the evolution of the physical field during complex construction processes.

[0030] In actual construction, the tunnel is usually constructed using a bottom-up, layered pouring process. The heat loss from the hydration of large volumes of concrete can cause volume shrinkage, resulting in non-uniform gaps at the interfaces between different media. To quantify this initial physical damage, the temporal grid layers are activated sequentially according to the construction timeline. The thermal shrinkage displacement of each temporal grid layer is calculated, and the thermal shrinkage displacement is accumulated normally at the contact interface between the surrounding rock and the concrete lining area to generate the initial void defect distribution.

[0031] At this stage, thermal shrinkage displacement was chosen as the primary control parameter because the temperature gradient is the direct physical driver of early interface delamination. By tracking the temperature field changes and shrinkage effects over time, the non-ideal contact boundary state of the chamber at the end of the construction period was established.

[0032] During its service life, the chamber undergoes compressed air filling and releasing operations, and the alternating temperature and pressure environment promotes the further evolution of early defects. Based on this physical process, alternating loads are applied to the contact interface to calculate the fatigue damage at the edge of the initial void defect distribution, and the normal contact stiffness distribution of the contact interface is updated based on the fatigue damage.

[0033] This process transforms the propagation and evolution of interfacial cracks at the microscopic level into a stiffness degradation index in macroscopic continuity mechanics, thus restoring the true mechanical support environment of underground chambers after long-term service.

[0034] Local loss of interface support stiffness can easily lead to stress imbalance in the internal steel sealing layer that provides a seal. Traditional single-point extremum determination is prone to getting stuck in local optima and cannot take into account the overall deformation coordination of thin-walled containers.

[0035] Based on this, the stress field of the steel sealing layer under alternating load is extracted. The required stiffness distribution at the contact interface is used as the design variable, and the optimization function is constructed with the goal of minimizing the stress variance of the steel sealing layer. The required stiffness distribution is then solved.

[0036] By establishing a reverse optimization mechanism guided by stress homogenization, the physical stiffness required to be added at the soil-rock interface to maintain the overall mechanical equilibrium of the system was quantified.

[0037] The continuous stiffness distribution field obtained from theoretical calculations is an abstract mechanical quantity that cannot be directly read by the mechanical equipment on the construction site. To achieve the connection between theoretical calculations and actual engineering, the stiffness distribution to be compensated is discretized, and grouting control parameters and dense mesh layout area parameters are output.

[0038] At the level of output result judgment, the system combines the upper limit of the consolidation physics of the grouting material for comprehensive judgment, avoiding process failure caused by relying solely on a single stiffness gap value. At the same time, it introduces the judgment logic of dense mesh auxiliary reinforcement for extremely weakened areas, thus completing the data closed loop from theoretical design to construction intervention.

[0039] See attached document Figure 2 The tunnel engineering encompasses both macroscopic rock mass distribution and local nonlinear contact interfaces. To resolve the contradiction between the redundancy of computational degrees of freedom brought about by full-scale spatial modeling and the accuracy requirements of local damage evolution, this scheme, based on finite element sub-model technology, performs bottom-level algorithm-level topological reconstruction and spatiotemporal decoupling on the technical features of constructing a three-dimensional geometric model of the tunnel and dividing it into finite element meshes and multiple temporal mesh layers as described in the claims. The specific implementation process includes steps S201 to S204.

[0040] Step S201: Construct a global three-dimensional geometric model including the surrounding rock, concrete lining area and steel sealing layer, and divide the global finite element mesh.

[0041] The computer system imports the geological outline of the chamber from on-site surveys using a computer-aided design interface, and establishes spatial geometric topological relationships by combining this with internal structural design drawings. During the mesh discretization stage, the system calls a mesh generator to discretize the solid region into three-dimensional continuous elements.

[0042] Since deep underground engineering is controlled by a complex geostress field, the system applies a three-dimensional initial geostress field and gravity load to the far-field outer boundary of the global finite element mesh, and runs a static analysis algorithm to solve the initial stress tensor in the entire computational domain. Then, by differentiating the spatial coordinates, the stress gradient distribution of the global three-dimensional geometric model under static analysis is extracted.

[0043] For the generation rules of the basic finite element mesh and the matrix assembly and solution of the linear elastic statics equations, those skilled in the art can use the existing finite element theory system.

[0044] In step S202, the stress distribution in the surrounding rock is non-uniform due to the influence of complex geological structures. If the computational resources are spread evenly across the entire area, it will not only be inefficient, but also easily mask the damage evolution characteristics of local high gradient areas.

[0045] To identify damage-sensitive areas, the system performs spatial filtering based on the calculated stress gradient distribution, extracting regions where the stress gradient distribution exceeds a preset threshold as local sub-models. In this process, a single strength extremum determination ignores the risk of complex failure under complex three-dimensional stress states. Therefore, the preset threshold is selected not only based on the standard value of the tensile strength of the surrounding rock material but also by incorporating a stress tolerance coefficient based on the Von Mises equivalent yield criterion.

[0046] Specifically, the preset threshold is set as the mathematical product of the standard value of tensile strength and the stress tolerance coefficient; among which, in order to balance the sensitivity of damage evolution capture and the computational efficiency of algorithm dimensionality reduction, the dimensionless value range of the stress tolerance coefficient is set to 0.65 to 0.85.

[0047] Specifically, the system uses the rock mass integrity index in the engineering survey report as a decay multiplier, multiplies it by the stress tolerance coefficient, and then applies it to the tensile strength standard value to obtain the final stress gradient determination threshold.

[0048] The system iterates through all cell integration points in the global 3D geometric model, and topologically merges the sets of cells whose stress gradient scalars exceed a preset threshold to generate independent local sub-model geometric boundaries. This truncation operation strips away the background rock mass, reducing the degree of freedom scale of subsequent multiphysics coupling iterations.

[0049] In step S203, after the local sub-model is truncated, its original mechanical equilibrium state is physically interrupted. In order to realistically reconstruct the physical constraints and load transfer effects of the external unmodeled rock mass on this local area, the system extracts the displacement field on the set of nodes at the intersection with the local sub-model in the global three-dimensional geometric model, and maps it as displacement boundary conditions to the truncated boundary of the local sub-model.

[0050] The system identifies all slave nodes on the outer boundary of the local sub-model and searches for the master elements that enclose these slave nodes in the global 3D geometric model, then uses finite element shape functions for spatial displacement interpolation. Considering that local mesh distortion may cause singularities in the Jacobian matrix of the shape function, leading to interpolation calculation failure, the system adds a tolerance check before interpolation.

[0051] When the determinant of the Jacobian matrix of the principal element approaches zero, the system skips shape function calculations and uses a weighted average method of nearest nodes within the local search radius to perform a fallback substitution for boundary displacements. The mathematical calculation logic for interpolation mapping when no singularities occur is as follows: ; in, The spatial coordinates on the boundary of the local sub-model are represented as The interpolated displacement vector obtained at the node; The representation encloses the coordinate point in the global 3D geometric model. The total number of nodes in the corresponding main unit is limited by the topology type of the main unit. In this embodiment, when a hexahedral unit is used, the value is 8. Characterizing the first principal unit The shape function of each node in coordinates The value at that location; Characterizing the principal element extracted in the global statics solution The global displacement vector of each node.

[0052] Through this shape function space mapping algorithm, the system transforms macroscopic deformation into the first type of forced constraint condition of local boundary, ensuring the mechanical equivalence of dimensionality reduction calculation.

[0053] Step S204: After establishing the high-precision spatial physical boundary, a time dimension needs to be further introduced to match the actual project progress.

[0054] Due to the constraints of large-volume concrete temperature control and formwork trolley technology on the engineering site, the concrete lining of the tunnel exhibits a bottom-up, progressive pouring characteristic. To accurately reproduce the time span effect of this process in numerical space, the system divides the concrete lining area in the local sub-model into multiple temporal grid layers based on the construction layer pouring height in the engineering design.

[0055] The system establishes elevation coordinate axes in the opposite direction of gravitational acceleration in the local sub-model, and uses the set layered pouring height values ​​as step sizes to perform Boolean slicing on the 3D mesh of the concrete lining area. Each slice generated by the segmentation corresponds to a specific pouring batch on site and is defined as an independent temporal mesh layer in terms of physical properties.

[0056] Each temporal grid layer is assigned a specific activation timestamp, thereby binding the static spatial geometric grid with the dynamic construction timeline, providing underlying grid support for subsequent tracking of the spatiotemporal evolution of the hydration heat gradient.

[0057] See attached document Figure 3In this embodiment, the heat of hydration release and subsequent volume shrinkage during the solidification and hardening process of large-volume concrete are the direct physical causes of the initial void defect between the surrounding rock and the concrete lining. Due to fluctuations in ambient temperature, humidity, and material composition at the construction site, using fixed theoretical thermodynamic parameters will lead to deviations between the simulation results and the actual deformation state.

[0058] To quantify this initial physical damage and improve the boundary realism of numerical simulation, the system activates the temporal grid layer sequentially according to the construction time axis and performs dynamic parameter inversion in combination with measured data. The specific implementation process includes steps S301 to S303.

[0059] Step S301: In order to synchronously reproduce the layered pouring process on site in the numerical model, the system uses the element birth and death technique to sequentially activate the current time-series grid layer to be calculated along the preset construction time axis.

[0060] When the current temporal grid layer to be calculated is activated and given an initial unhardened physical state, the temporal grid layer below it has already been cast and is in the hydration and heat release stage in actual engineering. The heat from the bottom layer will directly change the temperature evolution boundary of the upper material through interlayer conduction.

[0061] Based on this causal relationship, the computer system obtains the temperature rise time series data from the pre-embedded fiber optic temperature sensor within the activated time series grid layer located below the currently computed time series grid layer through a data interface. However, due to the inconsistency between the fixed sampling frequency of the field sensors and the adaptive calculation step size of the simulation program, direct data mapping can lead to abrupt changes in boundary conditions.

[0062] The system employs a cubic spline interpolation algorithm to resample and align the temperature rise time series data with the current calculation step size along the time axis. After resampling, a moving average filtering algorithm is used to filter out abnormal discrete points caused by high-frequency electrical noise from the sensors, thereby forming a continuous and smooth observation data vector in the time domain.

[0063] In step S302, after obtaining accurate observation data vectors, the system uses a Bayesian inference algorithm to update the specific heat capacity and heat generation rate coefficient of the current time-series grid layer to be calculated.

[0064] Given the highly nonlinear coupling characteristics of heat transfer and heat generation behavior in concrete, conventional gradient optimization algorithms can get stuck in local extrema. Bayesian inference can combine prior engineering knowledge to find the parameter combination that best conforms to the laws of physical evolution in the global probability space.

[0065] The specific heat capacity and heat generation rate coefficient are used as a joint parameter vector to be estimated, and the actual temperature rise response that has occurred in the lower layer is used to constrain the thermodynamic parameters of the upper region to be calculated.

[0066] The system constructs a likelihood function based on the Gaussian error distribution and, combined with the prior distribution of parameters, establishes a posterior probability distribution model. The mathematical calculation logic of this posterior probability distribution model is as follows: ; in, A joint parameter vector representing the specific heat capacity and the heat generation rate coefficient; The characterization includes Temperature rise time series data observation vector at each time sampling point; The posterior probability density function characterizing the joint parameter vector given the observation vector; The prior probability density function characterizing the joint parameter vector is specifically expressed as a multidimensional uniform distribution determined according to the concrete mix design specification. Characterization in the The actual observed temperature value fed back by the sensor at each sampling time; Representation in the current parameter vector The output of the lower finite element solver is the first The theoretically calculated temperature value corresponding to each time point; The statistical variance characterizing the observed data; The sign indicates proportionality.

[0067] To ensure the mathematical completeness of the algorithm when the observed data has minimal fluctuations or variances approaching zero, a non-zero, minimal positive real number is introduced into the denominator. As a numerical stability penalty term, its value range is set to 10. -6 ~10 -4 This is used to constrain the lower bound of numerical calculations and avoid division by zero anomalies.

[0068] The system employs a Markov chain Monte Carlo sampling algorithm to sample the aforementioned posterior probability distribution in a high-dimensional space. To ensure the global reliability of the sampling results in a statistical sense, the system computes multiple Markov chains in parallel and uses the Gelman-Rubin statistic to determine convergence.

[0069] When the Gelman-Rubin statistic is less than the set threshold of 1.2, it is determined that the sampling process has reached a stationary distribution. Subsequently, the maximum a posteriori probability estimate after convergence is extracted and used to replace the initial nominal specific heat capacity and heat generation rate coefficient of the current time-series grid layer to be calculated.

[0070] For the specific state transition criteria of the Markov chain Monte Carlo sampling algorithm, those skilled in the art can implement it using standard statistical theory.

[0071] In step S303, relying on the updated high-precision specific heat capacity and heat generation rate coefficient, the system inputs these parameters into the transient thermo-mechanical coupling solver to calculate the temperature field evolution of the current time-series grid layer to be calculated within the set hardening and curing cycle.

[0072] The change in temperature gradient causes the internal volume expansion and contraction of the structure. The system combines the linear expansion coefficient of the material with the creep constitutive model to transform the transient temperature field into thermal contraction strain at each grid node.

[0073] Constrained by external rock mass and interlayer bonding, thermal shrinkage strain ultimately manifests as non-uniform displacement at the contact interface. The system extracts the spatial node displacement vectors of the time-series grid layers at the end of the corresponding calculation step and calculates their thermal shrinkage displacement.

[0074] To eliminate slip interference parallel to the interface and accurately locate the void state, the system projects the thermal shrinkage displacement vector onto the normal direction of each local node at the contact interface between the surrounding rock and the concrete lining. Since the normal directions of the three-dimensional discrete mesh are discontinuous, the system uses the weighted average of the normal vectors of adjacent element surfaces sharing the same node as the transient reference normal vector for that local contact point, and performs a dot product projection operation.

[0075] The normal displacement projection values ​​generated at the contact interface by all sequentially activated and calculated temporal mesh layers are algebraically summed. To avoid false void detection caused by computer floating-point truncation errors, the system introduces an interface void tolerance threshold, which is based on the minimum feature size of the finite element mesh multiplied by 10. -3 Sure.

[0076] When the accumulated normal displacement scalar shows a separation trend away from the surrounding rock surface, and its absolute value is greater than the interface void tolerance threshold, the scalar value is recorded and a corresponding geometric gap is generated on the interface unit. The set of these gaps in the entire interface space constitutes the initial void defect distribution.

[0077] This distribution, as an initial boundary condition with physical memory, is directly transferred to the cyclic load calculation module during the subsequent service life.

[0078] See attached document Figure 4 In this embodiment, the artificial chamber undergoes a continuous cycle of high-temperature, high-pressure gas filling and low-temperature, low-pressure gas release during its decades-long service life. This alternating temperature and pressure load causes the interface material to be continuously subjected to alternating stress impacts, leading to the initiation and propagation of microcracks at the edges of the initial void defect distribution. Uniaxial static analysis cannot capture this type of damage effect that accumulates over time.

[0079] Based on the principle of continuous damage mechanics, the system introduces a cohesive force model to quantify the fatigue cumulative damage at the interface, transforming the physical degradation process at the micro level into a macroscopically calculable stiffness reduction index. The specific implementation process includes steps S401 to S403.

[0080] In step S401, the system embeds zero-thickness cohesive elements in batches along a continuous spatial grid at the interface between the surrounding rock and the concrete lining area to establish an interface cohesive model.

[0081] The cohesive model adopts the bilinear traction-separation rule as the constitutive control logic. Based on geological survey data and concrete material mechanical test reports, the system assigns initial normal contact stiffness and tangential contact stiffness to the cohesive unit.

[0082] After establishing the basic constitutive relations, the system applies a pressure and temperature alternating cyclic load spectrum for the compressed air energy storage equipment during its operation to the local sub-model. Due to the cyclic action of the alternating load and the geometrical abrupt changes caused by the initial voiding defect distribution, the cohesive elements at the defect geometric edges continuously bear stress singularity impacts, becoming regions with highly concentrated strain energy and sensitive areas for fatigue crack propagation.

[0083] For the topology generation algorithm of zero-thickness cohesive elements and the stiffness matrix assembly of the bilinear traction-separation law, those skilled in the art can use standard nonlinear finite element contact mechanics methods.

[0084] Step S402: Based on the acquired local high-frequency alternating stress state, the system calculates the cumulative fatigue damage variable of each cohesive integration point after a preset number of cycles.

[0085] Since the use of cycle-by-cycle direct time integration in long-cycle fatigue calculations leads to an exponential increase in computational cost, the system employs an envelope extrapolation algorithm based on cyclic cycles and an energy dissipation criterion for damage evolution calculation.

[0086] The evolution rate of the cumulative fatigue damage variable is driven by the variation in the strain energy release rate at the interface within each load cycle, and its mathematical calculation logic is as follows: ; in, Characterizing cumulative fatigue damage variables With load cycle number The rate of evolution of partial derivatives; and The fatigue evolution characteristic constants of the interface material are obtained by fitting the values ​​through standard interface fatigue shear and pull-out physical tests. In this embodiment... The value range is set to 10. -5 ~10 -3 , The value range is set to 1.5 to 3.0; The amplitude scalar of strain energy release rate, which is actually calculated from the work done by external force at the cohesive integral point under a single stress cycle spectrum, represents the strain energy release rate converted from the work done by external force at the corresponding cohesive integral point. Characterizes the critical fracture strain energy release rate required for the contact interface to reach a complete fracture state. To avoid variable amplitude scalars under extreme loads... Approaching the critical value indefinitely When the denominator approaches zero, it causes the numerical integral to diverge. The system introduces a minimal positive real number into the denominator. As the energy cutoff tolerance, its value is set to 10 based on the solver's floating-point precision. -6 .

[0087] The cumulative fatigue damage variable is obtained by integral calculation using the above formula. It is a dimensionless scalar with a value between 0 and 1. When the value is zero, it indicates that there is no physical damage to the contact interface. When the value reaches one, it indicates that the interface has completely debonded and broken.

[0088] Step S403: After obtaining the cumulative fatigue damage variable distribution of the entire interface space, the system performs physical reduction on the support stiffness matrix of the contact interface based on the variable.

[0089] Contact stiffness directly determines the interface's ability to transmit support reactions to the internal structure. The system utilizes the undamaged state proportion coefficient to transform the initial normal contact stiffness into a degraded, continuously distributed matrix of normal contact stiffness in space. The mathematical mapping logic for its physical reduction is as follows: ; in, Characterizing the contact interface in spatial coordinates The normal contact stiffness matrix at the node after cyclic service degradation; This represents the cumulative fatigue damage variable obtained by calculating the integral at this coordinate node through step S402; The initial normal contact stiffness matrix is ​​defined at this coordinate node; This represents the residual stiffness tolerance coefficient introduced to prevent numerical singularities in the overall system stiffness matrix. When the accumulated fatigue damage variables approach one, i.e., when physical cracking occurs at the interface, direct multiplication will cause the local stiffness matrix elements to completely degenerate to zero, thus leading to the failure of inverting the global stiffness matrix in the finite element method.

[0090] To maintain the numerical stability of the algorithm, the system introduces a residual stiffness tolerance coefficient for fallback limiting, with a fixed value of 10. -5The updated stiffness distribution matrix accurately reflects the physical state of local softening and support loss at the chamber interface during the later stages of service, providing an accurate mechanical boundary for subsequent assessment of the stress evolution of the steel sealing layer.

[0091] See attached document Figure 5 In this embodiment, the degradation of local interface stiffness directly alters the original load transfer path of the system, leading to uneven stress distribution and localized high stress concentration in the steel sealing layer, which bears the main sealing and seepage prevention function. To maintain the overall structural safety and stability over a long lifespan, the system establishes a reverse optimization mechanism guided by stress homogenization. By deriving the stress sensitivity matrix to interface stiffness, the system reduces the dimensionality of high-dimensional design variables and solves for the distribution of stiffness requiring compensation. The specific implementation process includes steps S501 to S504.

[0092] In step S501, the system uses the degraded spatial continuous distribution matrix of normal contact stiffness as the basic support boundary, runs the finite element solver, extracts the three-dimensional stress tensor at each mesh integration point of the steel sealing layer under alternating load, and transforms it into an equivalent stress field according to the VonMises yield criterion.

[0093] Compared to the one-sided approach of relying on the maximum stress at a single point in conventional design, controlling only local extreme values ​​in the steel sealing layer of a thin-walled large cylindrical structure will destroy the overall coordinated deformation capability of the structure.

[0094] To avoid local optima and achieve uniform stress distribution across the entire system, the system constructs an optimization objective function with the core objective of minimizing the spatial variance of the equivalent stress across the entire steel sealing layer. This objective drives the stress distribution on the pressure-bearing thin-shell material to tend towards uniformity.

[0095] Meanwhile, the system sets the stiffness distribution vector of each node at the interface between the surrounding rock and the concrete lining as a design variable.

[0096] In step S502, since the number of contact interface nodes in the large artificial chamber model is extremely large, if the stiffness of all nodes to be compensated is directly used as an independent variable in the optimization solution, it will cause the curse of dimensionality and lead to non-convergence of the calculation.

[0097] To reduce the optimization scale, the system introduces the adjoint variable method for stress sensitivity analysis. The core principle of the adjoint variable method is to construct virtual adjoint loads to transform the implicit derivative of the objective function with respect to design variables into an explicit algebraic equation for solution, thereby avoiding the large amount of repetitive finite element calculations generated by the finite difference method.

[0098] The system derives the partial derivatives of the equivalent stress of the steel sealing layer with respect to the compensable stiffness of local nodes at the contact interface, and constructs the sensitivity Jacobian matrix. The mathematical calculation logic is as follows: ; in, The sensitivity Jacobian matrix for the first The stress sensitivity column vector of each interface node; The equivalent stress vector characterizing the steel sealing layer; Characterizing the first The stiffness scalar that needs to be compensated at each interface node; The characterization is a transpose operator matrix that maps the differentiation of three-dimensional stress components to the equivalent stress gradient; Characterizing the displacement-strain geometry matrix in finite element analysis; Characterizes the global displacement vector of the system under external loads. This is because partial derivatives with respect to the displacement vector are required. The underlying algorithm must call the inverse of the global stiffness matrix. To prevent the global stiffness matrix from becoming ill-conditioned or even singular due to mesh distortion or stiffness degradation, the system injects a small Tikhonov regularization parameter (on the order of 10) into the diagonal elements of the global stiffness matrix before the inversion operation. -8 ).

[0099] By introducing adjoint vectors to solve the pseudo-loads of the global equilibrium equations for the design variables, the system obtains the sensitivity scalar values ​​for each node of the interface. For the pseudo-load assembly and partial derivative chain derivation of the finite element adjoint variable method, those skilled in the art can utilize existing structural sensitivity analysis theories.

[0100] Step S503: After obtaining the stress sensitivity scalar values ​​of all nodes, the system performs absolute value taking and sorts the set of values ​​in descending order.

[0101] If a small boundary stiffness perturbation cannot cause a change in the stress of the internal steel sealing layer, it indicates that the corresponding node contributes very little to the global objective optimization.

[0102] To avoid engineering applicability bias caused by relying on a single fixed percentage truncation, the system adopts a weighted judgment logic based on cumulative contribution. The system successively accumulates the sensitivity scalars after descending order and calculates its percentage of the total sensitivity of the entire domain.

[0103] The system filters out nodes whose stress sensitivity scalars are located in the lower columns of the sorted list based on a set cutoff rate threshold, and forces the corresponding stiffness scalars requiring compensation to be zero in subsequent iterations. This cutoff rate threshold is determined based on a balance between controlling resource consumption and retaining critical optimization paths. In this embodiment, the cutoff rate threshold is dynamically set to the boundary point where the cumulative sensitivity reaches 85%–90%. This cutoff operation maps the high-dimensional physical optimization space to a low-dimensional sensitivity subspace, reducing the iterative degrees of freedom in the optimization solution.

[0104] In step S504, within the reduced-dimensional sensitive subspace, the system combines the stress variance objective function of the steel sealing layer with the constraints of the overall mechanical equilibrium equation of the system, and uses a sequential quadratic programming algorithm to perform inverse optimization iteration.

[0105] In each optimization iteration, the system updates the global stiffness matrix of the finite element based on the current stiffness scalar that needs to be compensated, recalculates the stress field of the steel sealing layer, and extracts the current stress variance.

[0106] When the relative rate of change of stress variance output in two consecutive iterations is less than the set convergence tolerance threshold (which is 10 in this embodiment), -4 When the number of iterations reaches the preset maximum number of iterations (set to 200 in this embodiment), the system determines that the optimization process has reached steady-state convergence.

[0107] The system stops iterating and outputs the optimal compensation stiffness value for the high-sensitivity nodes in the current state. Simultaneously, it combines the zero values ​​of the filtered nodes to reconstruct and generate a compensation stiffness distribution covering the entire contact interface. This distribution precisely quantifies the amount of physical support compensation required at specific coordinate locations in space to offset service degradation effects.

[0108] See attached document Figure 6 In this embodiment, the stiffness distribution to be compensated, obtained through inverse optimization, is presented as a continuous set of node vectors in the numerical space. Due to limitations in the operational accuracy of on-site construction equipment and engineering costs, this continuous ideal stiffness compensation cannot be accurately reproduced point by point in the actual physical space.

[0109] In order to seamlessly integrate the theoretical calculation results into the on-site construction process, the system establishes a discretized mapping model based on spatial clustering and porous media mechanics theory, which transforms the continuous stiffness compensation requirements into specific grouting process parameters. The specific implementation process includes steps S601 to S603.

[0110] Step S601: The system extracts the three-dimensional spatial coordinates and corresponding stiffness scalar values ​​of the non-zero nodes in the stiffness distribution to be compensated, and constructs a feature dataset.

[0111] Grouting reinforcement operations at the engineering site spread outwards in a spherical or ellipsoidal shape from the borehole. Based on this, the system introduces a spatial clustering algorithm to divide the feature dataset into regions. Since the distribution of physical defects at the underground chamber interface exhibits a highly irregular spatial morphology, the system employs a density-based spatial clustering algorithm to effectively identify high-density stiffness compensation zones of arbitrary shapes. To ensure the stability and engineering feasibility of the clustering algorithm in complex 3D meshes, the system uses the spatial neighborhood search radius and the minimum number of core nodes as key constraints.

[0112] Specifically, the search radius of the spatial neighborhood is determined based on the average characteristic edge length of the finite element mesh of the contact surface, and is taken as 1.5 to 2.0 times the average characteristic edge length; the minimum core node number threshold is set to 3 to 5 based on the typical node density within the borehole influence range.

[0113] Because of the difference in physical dimensions between spatial geometric coordinates and stiffness gradient change rate, the system performs range standardization on the feature dataset before clustering, scaling all features to a dimensionless range of zero to one.

[0114] By calculating the Euclidean distance and stiffness gradient change rate between nodes in the normalized feature dataset, the system divides the global stiffness distribution requiring compensation into multiple independent discrete grouting reinforcement target areas. After division, the system extracts the geometric centroid coordinates of each discrete grouting reinforcement target area and directly maps them to the target borehole coordinates of the grouting boreholes in the engineering site.

[0115] In step S602, after determining the target borehole coordinates, the system needs to further quantify the grouting process parameters for each target area. The improvement in interfacial contact stiffness essentially depends on the physical filling of initial voids by the grout and its penetrating and compacting effect on the relaxed rock mass.

[0116] Before delving into specific mathematical derivations, the system constructs an analytical framework based on the elastic expansion mechanism in pore medium mechanics. Since the injected grout needs to overcome the initial geostress and pore water pressure of the rock mass, and achieves mechanical compression of the surrounding medium through volume expansion, this physical filling of the microscopic pore dimension is equivalent to an increase in the stiffness of the contact interface on a macroscopic level.

[0117] The system extracts the geometric gap volume of each discrete grouting reinforcement target area at its corresponding spatial location and uses it as the basic grouting volume of the corresponding borehole. In order to achieve the total required stiffness compensation output by reverse optimization, the injected grout needs to be split or penetrate into the microfractures of the rock mass under specific pressure.

[0118] The system, based on the theory of porous media mechanics, establishes a mathematical mapping model between grouting pressure and stiffness increment. The mathematical calculation logic of this mapping model is as follows: ; in, Characterizes the target grouting pressure for the corresponding discrete grouting reinforcement target area; Characterizes the initial pore water pressure of the rock mass at this spatial depth location; The volume compressibility coefficient characterizes the grouting medium, and its value is determined based on the specific grout material composition properties used in the project. In this embodiment, the value range is set to 10. -4 ~10 -3 MPa -1 ; Characterizes the set of interface nodes contained in the current discrete grouting reinforcement target area; Characterized by the algebraic sum of the stiffness scalars that need to be compensated for at all nodes within the target area; Characterizes the initial reference contact stiffness of the local interface.

[0119] Considering the possibility of extreme physical states where local regions may completely detach, causing the initial stiffness to approach zero, the system introduces a minimal positive real number into the denominator. As a numerical stability tolerance, its fixed value is 10. -6 This is to avoid program calculation crashes caused by division by zero exceptions.

[0120] This nonlinear logarithmic mapping relationship accurately reveals the physical attenuation characteristics between grouting pressure and stiffness increase. That is, as the interface stiffness increases, the medium pores gradually become saturated, and the additional grouting pressure required to achieve the same stiffness increment increases exponentially.

[0121] By solving the above model, the system transforms the abstract finite element stiffness compensation into a fluid dynamics pressure control index that can be directly read by on-site construction equipment.

[0122] In step S603, based on the above discretization mapping and mechanical equivalent transformation, the system obtains the grouting hole coordinates, foundation grouting volume, and target grouting pressure covering the entire chamber.

[0123] To address extreme conditions where local stiffness gaps are too large and grouting alone cannot effectively restore the continuity of the overall structure, the system further calculates the continuous void area of ​​each discrete grouting reinforcement target area on its two-dimensional projection plane. When this continuous void area exceeds a set critical safety area threshold, the system extracts the coordinates of the outermost boundary node of the target area and generates them as parameters for the dense mesh layout area.

[0124] The system rearranges these discrete process parameters (including grouting indicators and the range of dense mesh netting) according to the bottom-up spatial elevation sequence and procedural logic of on-site construction, generating a structured engineering process parameter matrix. This matrix contains a unique number for each borehole and its corresponding construction indicator parameters, and is imported via a data bus into the on-site 3D visualization construction management platform or the control terminal of the automated grouting trolley to guide precise targeted reinforcement operations on-site.

[0125] This data-driven closed-loop feedback mechanism ensures the stiffness coordination and stress safety of the artificial grouting chamber sealing structure throughout its long service life. For the data communication protocol of the automated grouting equipment and the front-end rendering mechanism of the visualization management platform, those skilled in the art can utilize existing industrial IoT communication standards and computer graphics interfaces.

[0126] See attached document Figure 7 In this embodiment, to ensure that the aforementioned degradation optimization compensation method based on measured data and cohesion model can be reliably and efficiently executed in the industrial field, the system has carried out targeted modular physical deployment on the underlying architecture of the computing device.

[0127] The various logical functional units, supported by the underlying hardware, interact with each other via a high-speed internal bus. This hardware-level data flow mechanism is designed to meet the stringent computational requirements of large-scale three-dimensional finite element nonlinear solutions. The specific construction and operation mechanism of this system includes steps S701 to S704.

[0128] In step S701, the system allocates an independent address mapping area within the server's motherboard storage space to instantiate the data acquisition and 3D modeling module. This module establishes a hard connection with the field sensor network via a standard industrial Ethernet interface to capture temperature and displacement time-series vectors in real time during the construction process.

[0129] To address the data discrepancy issue caused by inconsistent sampling frequencies from multiple sensors in the field, this module introduces a high-precision clock synchronization mechanism at the hardware level. Specifically, the system uses the lowest frequency sensor acquisition cycle as a benchmark and performs downsampling and timestamp alignment operations on the high-frequency time-series vectors to ensure strict matching of multiphysics boundary conditions in the time domain. After completing external data preprocessing, the integrated mesh generation engine within the module generates a continuous medium mesh and interface cohesive elements containing temporal activation information based on preset geometric feature parameters. The generated global mesh topology matrix and node coordinate matrix are serialized and written to the system's shared memory pool for subsequent computation modules to access without delay.

[0130] In step S702, to address the massive matrix operations involved in the defect dynamic simulation and service degradation evolution process, the system deploys a core mechanics solution engine within a heterogeneous computing architecture. This solution engine is logically divided into a defect dynamic simulation module and a service degradation calculation module. The defect dynamic simulation module is allocated to a central processing unit with multi-threading capabilities, and performs parallel sampling of thermodynamic parameters in the posterior probability space based on a Bayesian inference instruction set, thereby calculating the initial delamination geometric gap.

[0131] Based on the microscopic defect data obtained from the above deduction, the service degradation calculation module calls the floating-point acceleration core of the graphics processor and performs integral updates of the state variables of a large number of cohesive elements according to the bilinear traction-separation law and fatigue energy dissipation criterion. Since the calculation results of the two modules are physically causally related, the system adopts a cross-device memory access protocol to ensure the synchronization and data consistency of the stiffness matrix degradation reduction operation.

[0132] In step S703, after obtaining the degraded interface stiffness boundary, the system starts the inverse optimization and discrete mapping module residing in the cache. To overcome the computational bottleneck of the high-dimensional optimization space, the inverse optimization module incorporates a sensitivity analysis solver based on the adjoint variable method. This solver constructs and solves a large sparse Jacobian matrix, and performs physical memory-level forced zero-setting truncation on stiffness compensation nodes with extremely low sensitivity. Since the massive number of mesh nodes causes the Jacobian matrix to exhibit extremely sparse or even locally singular characteristics, the underlying computation unit uses a compressed row storage format to allocate memory. Meanwhile, the system detects that the determinant of the stiffness matrix is ​​close to zero (in this embodiment, the judgment tolerance is set to be less than 10). −12 When the problem occurs, a hardware fault-tolerant branch based on Singular Value Decomposition (SVD) to find the pseudo-inverse is automatically triggered to prevent the solution thread from crashing due to a serious exception.

[0133] Subsequently, the invoked sequential quadratic programming algorithm outputs a continuous vector of ideal stiffness distribution within this reduced-dimensional subspace. The discrete mapping module, through its internally embedded density space clustering logic algorithm unit, transforms the aforementioned continuous vector into a discrete engineering process parameter matrix composed of grouting hole coordinates, grouting volume, and target grouting pressure, and sends it to the on-site construction terminal via digital input / output ports.

[0134] Step S704: Based on the above-described modular deployment architecture, the present invention also provides a computer system and a matching computer-readable storage medium.

[0135] The computer system includes at least one processor and a memory communicatively connected to the processor. The memory, as a non-volatile storage medium, contains computer-executable instructions. When the processor reads and executes these instructions, the on / off states of the transistors in the hardware circuitry change periodically, thereby physically driving the underlying computing unit to precisely execute the entire process of mathematical operations and control operations covered in the aforementioned embodiments, including model building, parameter inversion, fatigue damage evolution, stiffness inverse optimization, and discrete mapping of grouting parameters.

[0136] For the selection of electronic components and circuit integration of the underlying hardware such as the central processing unit, volatile random access memory, and serial peripheral interface bus involved in the computer system, those skilled in the art can use existing general standards for computer architecture.

Claims

1. A method for optimizing the design of a large-diameter artificial compressed air energy storage chamber, characterized in that, Includes the following steps: A three-dimensional geometric model of the chamber, including the surrounding rock, the concrete lining area and the steel sealing layer, is constructed and divided into finite element meshes. The concrete lining area is divided into multiple temporal mesh layers according to the layered pouring height. The time-series grid layers are activated sequentially according to the construction timeline. The thermal shrinkage displacement of each time-series grid layer is calculated and accumulated normally at the contact interface between the surrounding rock and the concrete lining area to generate an initial void defect distribution. An alternating load is applied to the contact interface, fatigue damage at the edge of the initial void defect distribution is calculated, and the normal contact stiffness distribution of the contact interface is updated based on the fatigue damage. Extract the stress field of the steel sealing layer under the alternating load, use the required stiffness distribution of the contact interface as the design variable, construct an optimization function with the goal of minimizing the stress variance of the steel sealing layer, and solve for the required stiffness distribution. The stiffness distribution requiring compensation is discretized, and the grouting control parameters and the parameters of the dense mesh layout area are output.

2. The optimized design method for a large-diameter artificial compressed air energy storage chamber according to claim 1, characterized in that, The steps of constructing a three-dimensional geometric model of the chamber, including the surrounding rock, the concrete lining area, and the steel sealing layer, and dividing it into finite element meshes, and dividing the concrete lining area into multiple temporal mesh layers according to the layered pouring height, include: A global three-dimensional geometric model containing the surrounding rock, the concrete lining area and the steel sealing layer is constructed and divided into finite element meshes. The stress gradient distribution under static analysis is extracted. The region where the stress gradient distribution exceeds a preset threshold is taken as a local sub-model, and the displacement field at the intersection node with the local sub-model is mapped to the intercepted boundary of the local sub-model as the displacement boundary condition. Based on the layered pouring height, the concrete lining area in the local sub-model is divided into multiple temporal grid layers.

3. The optimized design method for a large-diameter artificial compressed air energy storage chamber according to claim 1, characterized in that, The steps of activating the temporal grid layers sequentially according to the construction timeline, calculating the thermal shrinkage displacement of each temporal grid layer, and performing normal accumulation at the contact interface between the surrounding rock and the concrete lining area to generate the initial void defect distribution include: Activate the current time-series grid layer to be calculated, and obtain the temperature rise time series data of the activated time-series grid layers below the current time-series grid layer; Based on the temperature rise time series data, the specific heat capacity and heat generation rate coefficient of the current time series grid layer to be calculated are updated using the Bayesian inference algorithm, and the temperature field and thermal contraction strain are calculated accordingly to generate the thermal contraction displacement. The thermal shrinkage displacement is accumulated normally on the contact interface. When the absolute value of the accumulation is greater than the void tolerance threshold, a geometric gap is generated. The geometric gaps over the entire domain constitute the initial void defect distribution.

4. The optimized design method for a large-diameter artificial compressed air energy storage chamber according to claim 1, characterized in that, The steps of applying alternating loads to the contact interface, calculating fatigue damage at the edge of the initial void defect distribution, and updating the normal contact stiffness distribution of the contact interface based on the fatigue damage include: A cohesive model is established at the contact interface and the alternating load is applied. The cumulative fatigue damage variable of the interface after a preset number of cycles is extracted as the fatigue damage. Using the undamaged state proportion coefficient, the initial normal contact stiffness of the contact interface is transformed into a degraded normal contact stiffness matrix, which serves as the updated normal contact stiffness distribution.

5. The optimized design method for a large-diameter artificial compressed air energy storage chamber according to claim 1, characterized in that, The steps of extracting the stress field of the steel sealing layer under alternating load, constructing an optimization function with the required stiffness distribution at the contact interface as the design variable and minimizing the stress variance of the steel sealing layer as the objective, and solving for the required stiffness distribution include: The partial derivatives of the equivalent stress of the steel sealing layer with respect to the required compensation stiffness of each node at the contact interface are derived using the adjoint variable method, and a sensitivity Jacobian matrix is ​​constructed to obtain the nodal stress sensitivity. The stress sensitivity of the nodes is sorted in descending order and the percentage is accumulated. The nodes in the latter column are filtered out according to the cutoff rate threshold. The stiffness of the nodes in the latter column that needs to be compensated is constrained to zero in subsequent iterations. Within the dimensionality-reduced sensitive subspace, the optimization function is solved using a sequential quadratic programming algorithm to obtain the stiffness distribution that needs compensation.

6. The optimized design method for a large-diameter artificial compressed air energy storage chamber according to claim 1, characterized in that, The step of discretizing the stiffness distribution requiring compensation and outputting grouting control parameters and dense mesh layout area parameters includes: A spatial clustering algorithm is used to divide the non-zero nodes in the stiffness distribution to be compensated into multiple discrete grouting target areas, and the centroid coordinates of each discrete grouting target area are extracted as the target hole position coordinates. The volume of the initial void defect at the corresponding position of each discrete grouting target area is extracted as the basic grouting volume. A mapping model between the target grouting pressure and the stiffness increment is established to solve for the target grouting pressure. The grouting control parameters are generated by combining the target hole position coordinates and the basic grouting volume. Calculate the continuous void area of ​​each discrete grouting target area. If the continuous void area is greater than the critical safety area threshold, extract the coordinates of the corresponding boundary nodes to generate the parameters of the dense mesh layout area.

7. An optimized design system for a large-diameter artificial compressed air energy storage chamber, characterized in that, The method for optimizing the design of a large-diameter artificial compressed air energy storage chamber as described in any one of claims 1-6 includes: The modeling and discretization module is used to construct a three-dimensional geometric model of the chamber, which includes the surrounding rock, the concrete lining area and the steel sealing layer, and to divide it into finite element meshes. The concrete lining area is divided into multiple temporal mesh layers according to the layered pouring height. The defect generation module is used to activate the time-series grid layers sequentially according to the construction time axis, calculate the thermal shrinkage displacement of each time-series grid layer, and perform normal accumulation on the contact interface between the surrounding rock and the concrete lining area to generate an initial void defect distribution. The stiffness update module is used to apply alternating loads to the contact interface, calculate the fatigue damage at the edge of the initial void defect distribution, and update the normal contact stiffness distribution of the contact interface based on the fatigue damage. The optimization solution module is used to extract the stress field of the steel sealing layer under the alternating load, construct an optimization function with the required stiffness distribution of the contact interface as the design variable and the goal of minimizing the stress variance of the steel sealing layer, and solve for the required stiffness distribution. The parameter output module is used to discretize the stiffness distribution that needs to be compensated and output grouting control parameters and dense mesh layout area parameters.

8. The optimized design system for a large-diameter artificial compressed air energy storage chamber according to claim 7, characterized in that, The defect generation module is specifically used for: Activate the current time-series grid layer to be calculated, and obtain the temperature rise time series data of the activated time-series grid layers below the current time-series grid layer; Based on the temperature rise time series data, the specific heat capacity and heat generation rate coefficient of the current time series grid layer to be calculated are updated using the Bayesian inference algorithm, and the temperature field and thermal contraction strain are calculated accordingly to generate the thermal contraction displacement. The thermal shrinkage displacement is accumulated normally on the contact interface. When the absolute value of the accumulation is greater than the void tolerance threshold, a geometric gap is generated. The geometric gaps over the entire domain constitute the initial void defect distribution.

9. The optimized design system for a large-diameter artificial compressed air energy storage chamber according to claim 7, characterized in that, The optimization solution module is specifically used for: The partial derivatives of the equivalent stress of the steel sealing layer with respect to the required compensation stiffness of each node at the contact interface are derived using the adjoint variable method, and a sensitivity Jacobian matrix is ​​constructed to obtain the nodal stress sensitivity. The stress sensitivity of the nodes is sorted in descending order and the percentage is accumulated. The nodes in the latter column are filtered out according to the cutoff rate threshold. The stiffness of the nodes in the latter column that needs to be compensated is constrained to zero in subsequent iterations. Within the dimensionality-reduced sensitive subspace, the optimization function is solved using a sequential quadratic programming algorithm to obtain the stiffness distribution that needs compensation.

10. The optimized design system for a large-diameter artificial compressed air energy storage chamber according to claim 7, characterized in that, The parameter output module is specifically used for: A spatial clustering algorithm is used to divide the non-zero nodes in the stiffness distribution to be compensated into multiple discrete grouting target areas, and the centroid coordinates of each discrete grouting target area are extracted as the target hole position coordinates. The volume of the initial void defect at the corresponding position of each discrete grouting target area is extracted as the basic grouting volume. A mapping model between the target grouting pressure and the stiffness increment is established to solve for the target grouting pressure. The grouting control parameters are generated by combining the target hole position coordinates and the basic grouting volume. Calculate the continuous void area of ​​each discrete grouting target area. If the continuous void area is greater than the critical safety area threshold, extract the coordinates of the corresponding boundary nodes to generate the parameters of the dense mesh layout area.

Images