Multi-scale modeling and fracture simulation method for compressed air energy storage caprock integrity assessment

Abstract

The application discloses a kind of multi-scale modeling and fracture simulation method for compressed air energy storage cover layer integrity evaluation, steps include: constructing multi-scale heterogeneous geomechanics model;Establish the fatigue damage constitutive relationship under rock loading and unloading condition;Establish fracture simulation control equation;Multi-scale calculation strategy and local adaptive grid are constructed;Cyclic loading condition simulation and integrity quantitative evaluation are carried out.The application is designed for the core damage mode of CAES gas storage low pressure fatigue, compared with general structure analysis software, the pertinence and accuracy of cover layer integrity evaluation are higher, the heterogeneous geological characteristics of cover layer are truly reflected, the process of crack initiation and propagation in cover layer under cyclic injection and production conditions is accurately simulated, and the advanced prediction and quantitative evaluation of cover layer integrity failure risk are realized.

Description

Technical Field

[0001] This invention relates to a multi-scale modeling and fracture simulation method for assessing the integrity of compressed air energy storage caprock, belonging to the interdisciplinary technical field of deep earth engineering and computational mechanics. Background Technology

[0002] Compressed air energy storage is a key supporting technology for building new power systems and realizing the large-scale, high-efficiency consumption of renewable energy. Its core infrastructure consists of gas storage facilities constructed using underground salt caverns, abandoned mines, depleted oil and gas reservoirs, or aquifers. During operation, the gas storage facility undergoes periodic pressurization and depressurization, resulting in significant pressure fluctuations within the storage tank. The caprock covering the gas storage facility acts as a natural barrier against high-pressure gas leakage. The integrity of the caprock directly determines the safety, sealing, and service life of the gas storage facility. Under long-term cyclic loading, existing micro-cracks within the caprock may expand, leading to internal lining failure and the formation of gas leakage channels. Therefore, accurate assessment of the caprock's integrity is crucial.

[0003] Currently, the assessment of caprock integrity in engineering mainly relies on empirical formulas and analytical solutions, as well as traditional numerical simulation methods. These two methods each have the following limitations:

[0004] Empirical formulas and analytical solutions: Based on elasticity mechanics or porosity elasticity theory, the caprock is simplified into a homogeneous, isotropic flat plate or thick-walled cylindrical model, and its safety factor is calculated using analytical formulas. These methods ignore the actual heterogeneity of the caprock, the distribution of natural fractures, and the influence of geological structures. Their assessment results are overly idealistic, often too conservative, or unable to identify local weak points, failing to meet the needs of refined design and risk assessment.

[0005] Traditional numerical simulation methods: use commercial finite element software to establish relatively simplified two-dimensional or three-dimensional models for stability analysis. Although it is an improvement over analytical methods, the following problems still exist: (1) Coarse modeling: geological models are mostly based on the assumption of layered homogeneity, and fail to make full use of data from boreholes, logging, geophysical exploration, etc. to construct a real three-dimensional heterogeneous field that can reflect the spatial variability of rock mechanical parameters; (2) Simple constitutive relations: linear elastic or elastoplastic constitutive relations are usually used, which cannot describe the progressive damage and failure process of rocks under cyclic loading, especially the complete mechanism of microcrack initiation, stable propagation and even unstable connection; (3) Contradiction between computational efficiency and accuracy: if crack propagation is to be simulated in detail, an extremely dense grid needs to be arranged in the potential failure area, resulting in a huge amount of computation for the whole model, which is difficult to apply to multiple cyclic simulations on a real engineering scale. Summary of the Invention

[0006] This invention provides a multi-scale modeling and fracture simulation method for the integrity assessment of compressed air energy storage caprock. This method can realistically reflect the heterogeneous geological characteristics of the caprock, accurately simulate the crack initiation and propagation process in the caprock under cyclic injection and production conditions, and realize the advanced prediction and quantitative assessment of the caprock integrity failure risk.

[0007] To achieve the above objectives, this invention provides a multi-scale modeling and fracture simulation method for integrity assessment of compressed air energy storage caprock, comprising the following steps:

[0008] S1. Construct a multi-scale heterogeneous geomechanical model;

[0009] S2. Establish the constitutive relationship of fatigue damage under rock loading and unloading conditions;

[0010] S3. Establish the control equations for fracture simulation;

[0011] S4. Construct a multi-scale computation strategy and a local adaptive mesh;

[0012] S5. Perform cyclic load condition simulation and quantitative integrity assessment.

[0013] Furthermore, the specific process of S1 is as follows:

[0014] S1.1 Generation of Random Parameter Field: Using the uniaxial compressive strength data of rock obtained from boreholes as the basic data, the expected strength values ​​of unsampled points in the entire caprock area are estimated using the Kriging interpolation method, forming a trend field of strength. To characterize the spatial variability of the parameters, a random perturbation field is superimposed on the trend field. The perturbation field is generated using a sequential Gaussian simulation method, the core of which is to assume that the parameters follow a Gaussian distribution and, based on the spatial correlation structure obtained from the borehole data variability function analysis using a spherical model, to sequentially generate random values ​​that conform to this distribution and correlation at the grid nodes. Finally, the material parameters of each grid cell are expressed as follows:

[0015] ;

[0016] in, This represents material parameters, such as uniaxial compressive strength, tensile strength, and elastic parameters; Trend is the trend value estimated by Kriging, which can be directly used through the predict method of the sklearn.gaussian_process.GaussianProcessRegressor instance in Python; Residual is the random residual value generated by the sequential Gaussian simulation, which is obtained through the instantiation of the gstools.SRF class in Python. Before the call, the variogram model is preset to a spherical model. Through this process, a three-dimensional random field reflecting the non-homogeneous distribution of the parameter space is obtained.

[0017] S1.2 Deterministic Macroscopic Structural Modeling: For macroscopic geological structures clearly identified through geological surveys, geophysical exploration (such as seismic and well logging), or remote sensing interpretation, including but not limited to faults, fold axial planes, large joints, lithological interfaces, and unconformities, a deterministic geometric modeling method is used for precise characterization. The process is as follows: First, based on the interpretation results, discrete control points or control lines of the structure in space are obtained. These data define the basic morphology and spatial distribution of the structure. Then, according to the geological genesis and geometric characteristics of the structure (such as the cross-sectional morphology of faults and the geometric pattern of folds), an appropriate surface construction algorithm is selected. For bedding structures, a continuous surface is generated using discrete point-based triangulation and surface interpolation algorithms. For linear structures, a spline curve fitting method is used to generate their centerlines. Finally, the generated geometric entities, i.e., surfaces or curves, are used as deterministic boundary conditions and embedded into the three-dimensional geological model. This allows for strict control of the geometric morphology and spatial location of macroscopic structures in the numerical model.

[0018] S1.3, Fusion Modeling: The geometric information of the deterministic structure is used as a cutting condition and imported into the 3D volume mesh model with generated random parameter fields. The processing method is as follows: First, all mesh elements are traversed through the spatial geometric intersection algorithm to accurately determine the spatial relationship between each element and the deterministic structural surface (such as the fault surface). For elements that are crossed by the surface, the program will perform mesh cutting and reclassification operations to generate new nodes and elements at the interface between the surface and the element, ensuring that the structural surface becomes the internal geometric boundary of the model. Subsequently, based on the signed distance function of the element centroid relative to the structural surface, the system automatically divides all elements, including newly generated elements, into different material regions or contact pairs: elements located on the positive side of the surface normal are marked as region A, elements on the negative side are marked as region B, and elements within a certain threshold, i.e., the influence bandwidth, are separately marked as transition region C, i.e., the fault influence zone.

[0019] After completing the region division, the parameter field is remapped and reassigned. For each grid cell, the final material parameters are no longer directly interpolated from the original random field at the cell's center point. Instead, a preset attribute correction function is applied based on the region identifier to which the cell belongs. This correction function is expressed as follows: the original random parameter field values ​​are... Multiply by a region-dependent correction factor ,Right now:

[0020] ;

[0021] in, The center point of the unit The region identification function, Representing intact rock mass regions A and B, Represents the region within the fault zone. Represents region C, which is influenced by the fault. For the corresponding material properties, a reduction or correction factor can be set (e.g., a setting can be made). , , This rule systematically characterizes the weakening or altering effect of deterministic structures on the material properties of their neighboring regions;

[0022] Ultimately, the resulting integrated three-dimensional geomechanical fundamental model is mathematically expressed as a piecewise function field controlled by region identifiers:

[0023] ;

[0024] in, The model is at a point in space. The set of mechanical properties, It is encapsulated The property correction function of the rule; this model is the input for subsequent simulation of fracture and seepage processes; This refers to the original random parameter field generated in step S1.1; I(x,y,z) is an indicator function used to determine the position of a point (x,y,z) relative to the deterministic structure. It is directly called through the query method of the Python library scipy.spatial.cKDTree to calculate the shortest distance from the point to the triangular mesh, i.e., the deterministic structure; f is an attribute correction function that adjusts P based on the value of I. random After correction, the function for f can be written in Python as follows:

[0025] def f(I_value, P_random_value):

[0026] correction_factors = { 0: 1.0, # Intact rock mass

[0027] 1: 0.3, # Within the fault zone

[0028] 2: 0.7 # Within the fault influence zone};

[0029] factor = correction_factors.get(I_value, 1.0);

[0030] return P_random_value * factor.

[0031] Furthermore, the specific process of S2 is as follows:

[0032] S2.1. An internal state variable, the cumulative fatigue damage factor Df, is introduced. This factor is related not only to the stress amplitude Δσ of the current load cycle, but also to the unloading path of the current cycle and the existing damage state d of the material. Its evolution equation is defined as:

[0033] ;

[0034] Among them, D f The current cumulative fatigue damage factor is given by N, where N is the number of cycles, C and m are material fatigue parameters, and σ0 is the reference stress. The unloading rate is represented by α and β, which are coupling coefficients. The core innovation of this model lies in the fact that the unloading term (1 + α *|dσ_unload / dt|) accelerates the fatigue damage accumulation rate during rapid unloading due to inertial friction and local stress redistribution at the microcrack surface; the damage coupling term (1 + β * d) reflects the weakening of the material's fatigue resistance by existing damage, i.e., the more severe the damage, the easier it is for fatigue cracks to propagate in subsequent cycles. This constitutive relation will be embedded as a source term in the subsequent phase-field evolution equation.

[0035] Furthermore, the specific process of S3 involves establishing governing equations applicable to porous media capping layers based on continuum mechanics:

[0036] Rock phase-field damage evolution equation: The phase-field fracture method is used to describe the damage and fracture process of the caprock material. A continuous scalar field variable d is introduced, and its evolution is driven by the principle of energy minimization. The rock fatigue damage effect established in step S2 is dynamically embedded into the phase-field evolution equation. The standard phase-field evolution equation is:

[0037] ;

[0038] Where H is the historical strain energy density, The critical energy release rate. The feature length is 3 to 5 times the minimum grid size. The Laplace operator is used, where d represents the current damage state of the material; the driving term H is expanded to include fatigue-related historical variables Hf. fat :

[0039] ;

[0040] Where, ψ elastic Where γ is the current elastic strain energy, D is the fatigue damage driving coefficient, and D is the current elastic strain energy. f As the current cumulative fatigue damage factor, It is the first This calculation step ensures that the damage driving force is not only related to the current and historical stress states, but also directly related to the fatigue damage D quantified from the cyclic load history.f Related. Therefore, even during the unloading phase of a single cycle or at lower stress levels, as long as fatigue damage D... f Through continuous accumulation, d may still evolve, thus achieving a precise description of the physical mechanism of stable microcrack propagation under cyclic loading, especially during the unloading stage.

[0041] Furthermore, S4, which balances the computational accuracy and efficiency of engineering-scale simulation, employs a global-local multi-scale computational framework and adaptive crack front densification. The specific process is as follows:

[0042] S4.1 Global Scale Analysis: An initial computational mesh is established on the entire 3D model of the cap layer, with a global mesh size L. global Determined based on the characteristic geometric dimensions of the model, it is defined as follows:

[0043] ;

[0044] Where H is the average thickness of the cap layer, and L c Let be the diameter of the gas storage cavity; on this grid, solve the coupled control equations established in step S3, that is, simultaneously solve the mechanical equilibrium equations and the modified phase field evolution equations:

[0045] ;

[0046] ;

[0047] in, Represents physical strength;

[0048] By discretizing and iteratively solving using the existing standard finite element method, the displacement field u at each global grid node is finally obtained. global Stress field σ global ε strain field global and damage field d global ; where d global It provides an intuitive representation of the damage state of materials, and its distribution cloud map is used to identify potential damage areas;

[0049] S4.2 Local Sub-model Intervention Criteria and Establishment: Setting the Damage Threshold d crit =0.25; when the global calculation results display the d of a node in a certain region. global Value reaches or exceeds d crit When an area is identified as a potentially critical damage zone, the system automatically creates a geometric boundary centered on that area, extending at least 3*L beyond the damaged area. globalThe local sub-model; the geometry of the local sub-model is inherited from the geometry of the global model in that region, and its material property field is obtained by mapping the heterogeneous field of the global model in that region to a finer mesh through bilinear interpolation; the mesh size within the local sub-model is L. global / 100;

[0050] S4.3 Coupling between Local Sub-model and Global Model: The local sub-model and global model are strongly coupled through displacement boundary conditions. The specific process is as follows: From the completed global model results, extract the displacement values ​​u of all nodes on the sub-model boundary at the current load step. global The displacement value is used as a mandatory displacement boundary condition and directly applied to the corresponding boundary nodes of the local sub-model. This method ensures that the deformation at the boundary of the sub-model is completely consistent with the overall deformation response of the global model, realizing the information transfer from the global to the local.

[0051] S4.4 Adaptive Mesh Refinement at the Crack Lead-Outside: In the high damage gradient region of the local sub-model or global model, adaptive mesh refinement based on the phase field gradient mode is adopted. The specific process is as follows: After each calculation step, the phase field gradient vector on all elements is calculated. Thus, its modulus field is obtained. ,Right now:

[0052] ;

[0053] Module length field It is obtained by solving the shape function derivatives of the phase field variable d within the element; the gradient threshold g is set. crit =0.99; During the calculation process, the value of each unit is checked in real time. Value, if If the element is located at the leading edge of a sharp crack change, the mesh of that element and its adjacent layer is immediately subdivided (e.g., the hexahedral element is divided into eight parts); when the crack propagates through, the area... The value fell back to g crit In this case, the mesh can be automatically merged and coarsened. This strategy ensures that computational resources are always focused on the crack tip regions that require the finest resolution.

[0054] Furthermore, the specific process of S5 is as follows:

[0055] S5.1 Boundary and Initial Condition Settings:

[0056] S5.1-1 Model Boundary Setting: The geometric boundary of the model is determined based on the actual engineering geological conditions. The bottom boundary is set in a geologically stable rock layer at twice the height of the gas storage cavity, which constrains the displacement in all directions, i.e., the fixed boundary.

[0057] S5.1-2 Lateral boundaries are set according to the regional geostress field characteristics: stress boundary conditions with linear gradient distribution are applied in the direction of maximum horizontal principal stress, and displacement constraints, i.e., normal fixed, are applied in the direction of minimum horizontal principal stress to simulate far-field stress constraints; the top boundary is the ground surface, set as a free boundary or an equivalent overlying stratum pressure is applied.

[0058] S5.1-3, Initial Geostress Field Generation: The initial geostress field is obtained through self-weight stress equilibrium calculation. Specifically, after applying the above boundary conditions, only gravity load is applied, and a static mechanical analysis is run once until the model reaches mechanical equilibrium. The stress field generated in the model at this time is the initial geostress field σ formed by the rock mass self-weight and boundary constraints. initial This stress field will serve as the initial state for subsequent cyclic load analysis;

[0059] S5.2 Load Application: A periodic pressure load P(t) conforming to the design operating law is applied to the inner wall surface of the gas storage cavity to simulate long-term injection-production cycles. P(t) is a value at the minimum operating pressure P min and maximum operating pressure P max Functions that change cyclically between:

[0060] ;

[0061] in, Average pressure; This refers to the pressure amplitude. The initial phase angle;

[0062] S5.3 Coupled Solution and Result Output: In the time domain, the coupled mechanical equilibrium equations and the modified phase field evolution equations are iteratively solved to simulate the response of the cap layer under multiple cyclic loads. Based on the simulation, the following key evaluation results are output:

[0063] S5.3-1, Damage / Crack Evolution Dynamic Diagram: Visually displays the entire process of fatigue damage accumulation and macroscopic crack initiation and propagation;

[0064] S5.3-2, Integrity Failure Risk Map: The system automatically identifies and marks whether the damage zone (d > 0.99) forms a continuous through path extending upward from the cavity boundary. If it forms, it is judged as a potential failure of sealing integrity.

[0065] S5.3-3 Fatigue Life Prediction Curve: Based on the simulation results, plot the curve of damage state d at key locations (such as the weakest point) of the caprock as a function of the number of cycles N. Define a critical damage value (such as d). failure = 0.8), to predict the fatigue life of the caprock under given operating conditions, i.e., the number of failure cycles;

[0066] S5.3-4, Recommendation for safe operating pressure window: Different pressure amplitudes, i.e., P, should be studied through parametric analysis. max With P min The combination of factors affects fatigue life. By plotting the life-pressure amplitude relationship curve, the maximum allowable pressure fluctuation range that ensures the integrity of the caprock within the target life is obtained.

[0067] This invention constructs a multi-scale heterogeneous geomechanical model, establishes the constitutive relationship of fatigue damage under rock loading and unloading conditions, establishes the control equation for fracture simulation, constructs a multi-scale calculation strategy and a local adaptive mesh, and performs cyclic loading simulation and quantitative integrity assessment. It combines an unloading path-dependent fatigue damage model with the phase-field fracture method, overcoming the limitation of traditional models that only consider damage in the loading segment. This more realistically reveals the progressive failure process of caprock under cyclic loading from a physical mechanism perspective. Simultaneously, the three-dimensional heterogeneous parametric field model constructed based on geostatistics transforms the identification of weak zones from empirical judgment to data-driven scientific prediction, significantly improving the model's realism. Furthermore, the combination of a global-local multi-scale framework and adaptive crack front densification technology achieves a holistic understanding and local focus, enabling large-scale, multi-cycle engineering problem simulation while ensuring high accuracy in crack propagation simulation. Finally, it outputs quantitative indicators directly applicable to engineering design, such as fatigue life prediction curves and safety pressure windows, achieving a leap from post-judgment to pre-judgment and scheme optimization. This invention is designed for the core failure mode of low-pressure fatigue in CAES gas storage facilities. Compared with general structural analysis software, it has higher pertinence and accuracy in assessing caprock integrity, truly reflects the heterogeneous geological characteristics of the caprock, and accurately simulates the process of crack initiation and propagation in the caprock under cyclic injection and production conditions, thus realizing advanced prediction and quantitative assessment of caprock integrity failure risk. Attached Figure Description

[0068] Figure 1 This is a flowchart of the process of the present invention. Detailed Implementation

[0069] The invention will now be further described with reference to the accompanying drawings.

[0070] like Figure 1 As shown, a multi-scale modeling and fracture simulation method for integrity assessment of compressed air energy storage caprock includes the following steps:

[0071] S1. Construct a multi-scale heterogeneous geomechanical model;

[0072] S2. Establish the constitutive relationship of fatigue damage under rock loading and unloading conditions;

[0073] S3. Establish the control equations for fracture simulation;

[0074] S4. Construct a multi-scale computation strategy and a local adaptive mesh;

[0075] S5. Perform cyclic load condition simulation and quantitative integrity assessment.

[0076] Example: Taking the Yunlong Lake Laboratory Compressed Air Energy Storage Project as an example, the process is as follows:

[0077] (1) Integrate the logging data of 3 wells, use sequential Gaussian simulation to generate a three-dimensional heterogeneous field of caprock elastic modulus and tensile strength, and deterministically model the top morphology of the gas storage tank;

[0078] (2) The fatigue parameters C, m, α, β and phase field parameters Gc, l0 of the rock were calibrated by indoor cyclic loading and unloading tests;

[0079] (3) Establish a global model and set a local sub-model with a mesh size of 0.5 meters to be started in the damaged area when d is greater than 0.25; apply a load with 500 cycles per year and pressure varying between 5 and 18 MPa to simulate 10 loading and unloading cycles; the simulation results show that damage concentration occurs in the local area after 15 years of operation, but no through cracks are formed until the 10th loading and unloading cycle.

[0080] (4) Output the dN curve of the region. It is predicted that under the current working conditions, it will take more than 10,000 loading and unloading cycles to reach the critical damage. The integrity of the cap layer is evaluated to ensure that the design requirements are met. Through simulation by adjusting the pressure amplitude, it is found that reducing the upper limit pressure to 10 MPa can extend the predicted life to more than 10,000,000 loading and unloading cycles, which can be regarded as an infinite life. This provides a quantitative basis for optimizing the operation plan.

Claims

1. A multi-scale modeling and fracture simulation method for assessing the integrity of compressed air energy storage caprock, characterized in that, Includes the following steps: S1. Construct a multi-scale heterogeneous geomechanical model; S2. Establish the constitutive relationship of fatigue damage under rock loading and unloading conditions; S3. Establish the control equations for fracture simulation; S4. Construct a multi-scale computation strategy and a local adaptive mesh; S5. Perform cyclic load condition simulation and quantitative integrity assessment.

2. The multi-scale modeling and fracture simulation method for integrity assessment of compressed air energy storage caprock as described in claim 1, characterized in that, The specific process of S1 is as follows: S1.1 Generation of Random Parameter Field: Using the uniaxial compressive strength data of rock obtained from boreholes as the basic data, the expected strength values ​​of unsampled points in the entire caprock area are estimated using the Kriging interpolation method, forming a trend field of strength. To characterize the spatial variability of the parameters, a random perturbation field is superimposed on the trend field. The perturbation field is generated using a sequential Gaussian simulation method, the core of which is to assume that the parameters follow a Gaussian distribution and, based on the spatial correlation structure obtained from the borehole data variability function analysis using a spherical model, to sequentially generate random values ​​that conform to this distribution and correlation at the grid nodes. Finally, the material parameters of each grid cell are expressed as follows: ; in, The material parameters are represented by Trend, which is the trend value estimated by Kriging, and Residual is the random residual value generated by the sequential Gaussian simulation. The variogram model is preset to a spherical model before calling. Through this process, a three-dimensional random field reflecting the non-homogeneous distribution of the parameter space is obtained. S1.2 Deterministic Macroscopic Structural Modeling: For macroscopic geological structures clearly identified through geological surveys, geophysical exploration, or remote sensing interpretation, including faults, fold axial planes, large joints, lithological interfaces, and unconformities, a deterministic geometric modeling method is used for precise characterization. The process is as follows: First, based on the interpretation results, discrete control points or control lines of the structure in space are obtained. These data define the basic morphology and spatial distribution of the structure. Then, according to the geological genesis and geometric characteristics of the structure, an appropriate surface construction algorithm is selected. For bedding structures, a continuous surface is generated using discrete point-based triangulation and surface interpolation algorithms. For linear structures, a spline curve fitting method is used to generate their centerlines. Finally, the generated geometric entities, i.e., surfaces or curves, are used as deterministic boundary conditions and embedded into the three-dimensional geological body model. S1.3, Fusion Modeling: The geometric information of the deterministic structure is used as a cutting condition and imported into the 3D volumetric mesh model with a generated random parameter field. The processing method is as follows: First, all mesh elements are traversed through the spatial geometric intersection algorithm to accurately determine the spatial relationship between each element and the deterministic structure surface. For elements that are crossed by the surface, the program performs mesh cutting and reclassification operations to generate new nodes and elements at the interface between the surface and the element, ensuring that the structure surface becomes the internal geometric boundary of the model. Subsequently, based on the signed distance function of the element centroid relative to the structure surface, the system automatically divides all elements, including newly generated elements, into different material regions or contact pairs: elements located on the positive side of the surface normal are marked as region A, elements on the negative side are marked as region B, and elements within a certain threshold, i.e., the influence bandwidth, are marked as transition region C, i.e., the fault influence zone. After completing the region division, the parameter field is remapped and reassigned. For each grid cell, the final material parameters are no longer directly interpolated from the original random field at the cell's center point. Instead, a preset attribute correction function is applied based on the region identifier to which the cell belongs. This correction function is expressed as follows: the original random parameter field values ​​are... Multiply by a region-dependent correction factor ,Right now: ； in, The center point of the unit The region identification function, Representing intact rock mass regions A and B, Represents the region within the fault zone. Represents region C, which is influenced by the fault. This refers to the corresponding material property reduction or correction factor; Ultimately, the resulting integrated three-dimensional geomechanical fundamental model is mathematically expressed as a piecewise function field controlled by region identifiers: ； in, The model is at a point in space. The set of mechanical properties, It is encapsulated The property correction function of the rule; this model is the input for subsequent simulation of fracture and seepage processes; This refers to the original random parameter field generated in step S1.1; I(x,y,z) is an indicator function used to determine the position of a point (x,y,z) relative to the deterministic structure. It is directly called through the query method of the Python library scipy.spatial.cKDTree to calculate the shortest distance from the point to the triangular mesh, i.e., the deterministic structure; f is an attribute correction function that adjusts P based on the value of I. random After correction, the function for f can be written in Python as follows: def f(I_value, P_random_value): correction_factors = { 0: 1.0, # Intact rock mass 1: 0.3, # Within the fault zone 2: 0.7 # Within the fault influence zone}; factor = correction_factors.get(I_value, 1.0); return P_random_value * factor.

3. The multi-scale modeling and fracture simulation method for integrity assessment of compressed air energy storage caprock as described in claim 1, characterized in that, The specific process of S2 is as follows: S2.

1. An internal state variable, the cumulative fatigue damage factor Df, is introduced. This factor is related not only to the stress amplitude Δσ of the current load cycle, but also to the unloading path of the current cycle and the existing damage state d of the material. Its evolution equation is defined as: ； Among them, D f The current cumulative fatigue damage factor is given by N, where N is the number of cycles, C and m are material fatigue parameters, and σ0 is the reference stress. The unloading rate is represented by α and β, which are coupling coefficients.

4. The multi-scale modeling and fracture simulation method for integrity assessment of compressed air energy storage caprock according to claim 3, characterized in that, The specific process of S3 involves establishing governing equations applicable to porous media capping layers based on continuum mechanics: Rock phase-field damage evolution equation: The phase-field fracture method is used to describe the damage and fracture process of the caprock material. The rock fatigue damage effect established in step S2 is dynamically embedded into the phase-field evolution equation. The standard phase-field evolution equation is as follows: ； Where H is the historical strain energy density, The critical energy release rate. The feature length is 3 to 5 times the minimum grid size. The Laplace operator is used, where d represents the current damage state of the material; the driving term H is expanded to include fatigue-related historical variables Hf. fat : ； Where, ψ elastic Where γ is the current elastic strain energy, D is the fatigue damage driving coefficient, and D is the current elastic strain energy. f As the current cumulative fatigue damage factor, It is the first One calculation step.

5. The multi-scale modeling and fracture simulation method for integrity assessment of compressed air energy storage caprock according to claim 4, characterized in that, S4, which balances the computational accuracy and efficiency of engineering-scale simulation, employs a global-local multi-scale computational framework and adaptive crack front densification. The specific process is as follows: S4.1 Global Scale Analysis: An initial computational mesh is established on the entire 3D model of the cap layer, with a global mesh size L. global Determined based on the characteristic geometric dimensions of the model, it is defined as follows: ； Where H is the average thickness of the cap layer, and L c Let be the diameter of the gas storage cavity; on this grid, solve the coupled control equations established in step S3, that is, simultaneously solve the mechanical equilibrium equations and the modified phase field evolution equations: ； ； in, Represents physical strength; By discretizing and iteratively solving using the existing standard finite element method, the displacement field u at each global grid node is finally obtained. global Stress field σ global ε strain field global and damage field d global ; where d global It provides an intuitive representation of the damage state of materials, and its distribution cloud map is used to identify potential damage areas; S4.2 Local Sub-model Intervention Criteria and Establishment: Setting the Damage Threshold d crit When the global calculation results display the d of a node in a certain region... global Value reaches or exceeds d crit When an area is identified as a potentially critical damage zone, the system automatically creates a geometric boundary centered on that area, extending at least 3*L beyond the damaged area. global The local sub-model; the geometry of the local sub-model is inherited from the geometry of the global model in that region, and its material property field is obtained by mapping the heterogeneous field of the global model in that region to a finer mesh through bilinear interpolation; the mesh size within the local sub-model is L. global / 100; S4.3 Coupling between Local Sub-model and Global Model: The local sub-model and global model are strongly coupled through displacement boundary conditions. The specific process is as follows: From the completed global model results, extract the displacement values ​​u of all nodes on the sub-model boundary at the current load step. global The displacement value is then applied directly to the corresponding boundary nodes of the local sub-model as a mandatory displacement boundary condition. S4.4 Adaptive Mesh Refinement at the Crack Lead-Outside: In the high damage gradient region of the local sub-model or global model, adaptive mesh refinement based on the phase field gradient mode is adopted. The specific process is as follows: After each calculation step, the phase field gradient vector on all elements is calculated. Thus, its modulus field is obtained. ,Right now: ； Module length field It is obtained by solving the shape function derivatives of the phase field variable d within the element; the gradient threshold g is set. crit During the calculation process, the performance of each unit is checked in real time. value, If the element is located at the forefront of a sharp crack change, the mesh of that element and its adjacent layer is immediately subdivided; when the crack propagates through, the area... The value fell back to g crit In the following cases, the mesh is automatically merged and coarsened.

6. The multi-scale modeling and fracture simulation method for integrity assessment of compressed air energy storage caprock as described in claim 5, characterized in that, The specific process of S5 is as follows: S5.1 Boundary and Initial Condition Settings: S5.1-1 Model Boundary Setting: The geometric boundary of the model is determined based on the actual engineering geological conditions. The bottom boundary is set in a geologically stable rock layer at twice the height of the gas storage cavity, which constrains the displacement in all directions, i.e., the fixed boundary. S5.1-2 Lateral boundaries are set according to the regional geostress field characteristics: stress boundary conditions with linear gradient distribution are applied in the direction of maximum horizontal principal stress, and displacement constraints, i.e., normal fixed, are applied in the direction of minimum horizontal principal stress to simulate far-field stress constraints; the top boundary is the ground surface, set as a free boundary or an equivalent overlying stratum pressure is applied. S5.1-3, Initial Geostress Field Generation: The initial geostress field is obtained through self-weight stress equilibrium calculation. Specifically, after applying the above boundary conditions, only gravity load is applied, and a static mechanical analysis is run once until the model reaches mechanical equilibrium. The stress field generated in the model at this time is the initial geostress field σ formed by the rock mass self-weight and boundary constraints. initial This stress field will serve as the initial state for subsequent cyclic load analysis; S5.2 Load Application: A periodic pressure load P(t) conforming to the design operating law is applied to the inner wall surface of the gas storage cavity to simulate long-term injection-production cycles. P(t) is a value at the minimum operating pressure P min and maximum operating pressure P max Functions that change cyclically between: ； in, Average pressure; This refers to the pressure amplitude. The initial phase angle; S5.3 Coupled Solution and Result Output: In the time domain, the coupled mechanical equilibrium equations and the modified phase field evolution equations are iteratively solved to simulate the response of the cap layer under multiple cyclic loads. Based on the simulation, the following key evaluation results are output: S5.3-1, Damage / Crack Evolution Dynamic Diagram: Visually displays the entire process of fatigue damage accumulation and macroscopic crack initiation and propagation; S5.3-2, Integrity Failure Risk Map: The system automatically identifies and marks whether the damage zone forms a continuous through path extending upward from the cavity boundary. If it forms, it is determined to be a potential failure of the sealing integrity. S5.3-3 Fatigue life prediction curve: Based on the simulation results, plot the curve of the damage state d at the key location of the caprock as a function of the number of cycles N. By defining the critical damage value, predict the fatigue life of the caprock under a given operating condition, i.e. the number of failure cycles. S5.3-4, Recommendation for safe operating pressure window: Different pressure amplitudes, i.e., P, should be studied through parametric analysis. max With P min The combination of factors affects fatigue life. By plotting the life-pressure amplitude relationship curve, the maximum allowable pressure fluctuation range that ensures the integrity of the caprock within the target life is obtained.

Images