An underground powerhouse cavern group multi-scale coupling adaptive numerical simulation method

Abstract

The application discloses a kind of underground powerhouse cavern group multiscale coupling adaptive numerical simulation methods, first, establish macro continuum model, calculate stress field, displacement field and plastic damage zone under excavation unloading;Then build local refinement trigger index including displacement anomaly, stress concentration, plastic damage and structure surface sensitivity, identify trigger area and extract local submodel;Subsequently, using non-matching grid boundary mapping method, boundary displacement, boundary stress and principal stress direction are transferred, and discrete fracture network model is constructed to represent fault, joint and fracture distribution;Key cracks are screened through crack propagation potential evaluation function, and crack initiation, propagation and penetration are calculated by embedding expansion finite element model;Finally, the equivalent stiffness degradation parameter is calculated, the deformation modulus of the macro model is updated by feedback, and the stability is recalculated and the convergence is judged;The application can realize the two-way coupling of macro response and local crack damage evolution, thereby improving the accuracy, efficiency and reliability of stability analysis.

Description

Technical Field

[0001] This invention relates to the field of numerical calculation and stability evaluation of rock mechanics in underground engineering, and in particular to a multi-scale coupled adaptive numerical simulation method for underground powerhouse cavern groups. Background Technology

[0002] As underground hydropower projects extend to deeper levels, the scale and complexity of underground powerhouse cavern complexes are continuously increasing. Influenced by high ground stress, high osmotic pressure, and complex geological structures, after excavation and unloading, structural surfaces such as faults, joints, and fissures in the surrounding rock are prone to stress concentration and damage evolution. This leads to the initiation, propagation, and eventual connection of cracks, weakening the overall stability of the local surrounding rock and the cavern complex. Therefore, conducting highly reliable simulations of the mechanical response of complex fractured rock masses is of great significance for the design, construction control, and stability evaluation of underground powerhouse cavern complexes.

[0003] Currently, numerical simulation methods for the stability of underground cavern groups mainly include continuum numerical methods, discontinuum numerical methods, and macro-local refinement coupling methods. However, when facing ultra-large-scale, complex fractured rock mass underground engineering projects, existing methods still have the following shortcomings:

[0004] (1) Traditional finite element method (FEM) or finite difference method (FDM) is usually based on the assumption of continuous medium and is suitable for the analysis of the overall stress field, displacement field and plastic damage zone of underground cavern groups. However, real rock masses have obvious heterogeneity, anisotropy and discontinuity. Traditional continuum models are difficult to explicitly characterize the crack initiation, propagation, penetration and local block instability process under the control of faults, joints and fissures, resulting in insufficient ability to identify the risks of local block shedding, fissure penetration and brittle failure.

[0005] (2) Discontinuous methods such as the Discrete Element Method (DEM), Discrete Fracture Network (DFN) model, and Extended Finite Element Method (XFEM) can describe the discontinuous deformation and crack propagation process of fractured rock masses. However, for ultra-large-scale underground powerhouse cavern groups, if a large number of fractures are explicitly constructed and crack propagation calculations are carried out in the entire domain, the computational scale and solution cost will be significantly increased, making it difficult to meet the needs of multi-condition and multi-stage analysis of large-scale projects. In addition, existing rock mass defect modeling methods mostly focus on the geometric expression of fractures, joints, or pores, and have not yet effectively solved the problems of adaptive identification of local dangerous areas, macro-local boundary condition transfer, and back-feedback of local crack propagation results for overall stability evaluation.

[0006] (3) Although existing macroscopic calculation and local detailed analysis methods attempt to balance efficiency and accuracy, the scope of local detailed analysis often relies on engineering experience or is determined solely based on single indicators such as displacement and stress. For underground powerhouse cavern groups affected by high ground stress, excavation unloading, plastic damage, and structural surfaces, a single criterion is insufficient to accurately identify hidden damage zones, structural surface control zones, and potential fracture penetration points. This can easily lead to deviations in the delineation of local detailed modeling areas, thereby affecting the accuracy of local rock damage and overall stability evaluation.

[0007] (4) Although existing multi-scale or continuous-discontinuous coupled models can simulate the fracturing process of the surrounding rock in the cavern to a certain extent, they still often face problems such as different mesh scales, inconsistent node positions, and non-coincident boundary nodes when transferring data between the macroscopic continuum model and the local refined model. If the boundary displacement, boundary stress, or principal stress direction is not accurately transferred, it can easily lead to distortion of the local model's stress environment and affect the reliability of the crack propagation calculation results. At the same time, in the process of connecting the discrete fracture network model and the extended finite element model, if all fractures are directly input into the extended finite element model, the calculation scale will be too large; if the initial crack is selected based solely on manual experience or fracture length, factors such as the distance to the cavern, stress driving, fracture connectivity, and macroscopic damage background may be ignored, resulting in insufficient identification of key damaging fractures.

[0008] Therefore, it is necessary to develop a multi-scale coupled adaptive numerical simulation method for underground powerhouse cavern groups to solve the above problems. Summary of the Invention

[0009] The purpose of this invention is to design a multi-scale coupled adaptive numerical simulation method for underground powerhouse cavern groups in order to solve the above-mentioned problems.

[0010] The present invention achieves the above objectives through the following technical solutions:

[0011] A multi-scale coupled adaptive numerical simulation method for underground powerhouse cavern groups includes the following steps:

[0012] S1. Establish a macroscopic continuum model of the underground powerhouse cavern group and calculate the overall stress field, displacement field and plastic damage zone under the excavation unloading condition;

[0013] S2. Based on the output data of step S1, construct a local refined triggering index that includes displacement anomaly, stress concentration, plastic damage and structural surface sensitivity, identify the triggering area, and extract local sub-models;

[0014] S3. Based on the local sub-model extracted in step S2, the non-matching mesh boundary mapping method is used to map the boundary displacement, boundary stress and principal stress direction spatial interpolation obtained from the macro model to the boundary of the local refined sub-model.

[0015] S4. Construct a DFN model within the local sub-model extracted in step S2 to explicitly characterize the spatial distribution of faults, joints, and fractures;

[0016] S5. Construct a crack propagation potential evaluation function, select key cracks from the DFN model generated in step S4, and identify cracks with crack propagation potential index higher than a preset threshold as key cracks, and embed them into the XFEM model as initial cracks.

[0017] S6. Based on the stress boundary conditions of the local sub-model in step S3, perform numerical calculations of the XFEM crack initiation, propagation, and penetration evolution process, and output the local crack propagation characteristics.

[0018] S7. Extract the local crack propagation features from step S6, calculate the equivalent stiffness degradation parameters, and feed back the equivalent stiffness degradation parameters to update the deformation modulus of the corresponding region of the macroscopic continuum model in step S1.

[0019] S8. Based on the updated macroscopic continuum model from step S7, recalculate the overall stress field, displacement field, and plastic damage zone of the underground powerhouse cavern group, and determine whether the convergence condition is met. If not, return to step S2 to re-identify and iteratively calculate. If yes, output the final stability evaluation and failure prediction results.

[0020] Specifically, in step S1, the macroscopic continuum model of the underground powerhouse cavern group is established based on the engineering layout of the underground powerhouse cavern group, the cavern excavation outline, the surrounding rock classification results, the geological structure distribution, the initial geostress field, the model boundary conditions, and the excavation unloading conditions. The macroscopic continuum model divides the surrounding rock of the underground powerhouse cavern group into continuous medium calculation units according to the surrounding rock classification results. Each calculation unit is assigned a partition value according to the surrounding rock category corresponding to different surrounding rock classifications. Each type of surrounding rock is assigned corresponding rock mechanics parameters, including deformation modulus, Poisson's ratio, cohesion, internal friction angle, tensile strength, compressive strength, and natural density.

[0021] Specifically, step S2 includes:

[0022] The macroscopic continuum model of the underground powerhouse cavern group is divided into several computational regions, and the local refined triggering index is calculated for each computational region.

[0023] The local refinement triggering index is obtained by coupling the macroscopic mechanical response index and the structural surface sensitivity index, and its expression is:

[0024] ;

[0025] in, Let be the local refinement trigger index for the i-th computational region. Let be the macroscopic mechanical response index of the i-th computational region. Let i be the structural surface sensitivity index for the i-th computational region. This is the structural surface sensitivity amplification factor;

[0026] The macroscopic mechanical response index is used to characterize the deformation anomalies, stress concentrations, and plastic damage levels in the computational domain under excavation unloading. Its expression is:

[0027] ;

[0028] in, Let be the displacement anomaly factor for the i-th computational region. Let be the stress concentration factor for the i-th computational region. Let be the plastic damage factor for the i-th computational region. , , The weighting coefficients for the displacement anomaly factor, stress concentration factor, and plastic damage factor are, in order, and satisfy the following:

[0029] ;

[0030] The displacement anomaly factor is determined by both the displacement amplitude and the displacement gradient, and its expression is as follows:

[0031] ;

[0032] in, Let i be the displacement value of the i-th calculation region. This represents the maximum displacement value across all computational domains in the macroscopic continuum model. Let be the displacement gradient of the i-th computational region. This represents the maximum displacement gradient across all computational domains in the macroscopic continuum model. and These are the weighting coefficients for the displacement magnitude and the displacement gradient, respectively, and they satisfy:

[0033] ;

[0034] The stress concentration factor is determined by both the maximum principal stress and the deviatoric stress, and its expression is:

[0035] ;

[0036] in, The maximum principal stress in the i-th computational region is... This represents the maximum principal stress across all computational domains in the macroscopic continuum model. For the deviatoric stress in the i-th calculation region, This represents the maximum deviatoric stress across all computational domains in the macroscopic continuum model. and These are the weighting coefficients for the maximum principal stress and the deviatoric stress, respectively, and they satisfy the following:

[0037] ;

[0038] The plastic damage factor is determined by the damage variable output from the macroscopic continuum model, and its expression is:

[0039] ;

[0040] in, Let i be the damage variable for the i-th computational region. The maximum damage variable across all computational domains in the macroscopic continuum model;

[0041] When the macroscopic continuum model does not output continuous damage variables, the plastic damage factor is determined by the plastic state of the computational domain, and its expression is:

[0042] ;

[0043] in, Let be the set of computational regions determined to be in a plastic state in the macroscopic continuum model; if the i-th computational region belongs to the plastic region, then =1; if the i-th computational region does not belong to the plastic region, then =0;

[0044] The structural surface sensitivity index is used to characterize the control effect of faults, joints, and fractures on local damage. Its expression is:

[0045] ;

[0046] in, Let be the structural surface density factor of the i-th computational region. Let be the near-hole distance factor for the i-th computational region. Let be the intersection factor of the structural surfaces in the i-th computational region. Let be the stress direction unfavorable factor for the i-th computational region. , , , These are the weighting coefficients for the structural surface density factor, near-hole distance factor, structural surface intersection factor, and stress direction unfavorable factor, respectively, and they satisfy the following:

[0047] ;

[0048] The structural surface density factor is used to characterize the development of faults, joints, and fractures within the computational region, and its expression is:

[0049] ;

[0050] in, Let be the surface density of the structure within the i-th computational region. The maximum structural surface density across all computational regions;

[0051] The proximity factor is used to characterize the spatial proximity between the structural surface and the excavation boundary of the underground powerhouse cavern. Its expression is as follows:

[0052] ;

[0053] in, Let be the minimum distance from the main structural surface to the excavation boundary of the cavern within the i-th calculation region. The maximum control distance at which the structural face affects the stability of the surrounding rock of the cavern; when Greater than season =0;

[0054] The structural plane intersection factor is used to characterize the degree to which multiple sets of structural planes intersect within the computational region to form potential block boundaries or through-failure pathways. Its expression is:

[0055] ;

[0056] in, Let i be the number of structural surfaces intersecting within the i-th computational region. The maximum number of structural surface intersections across all computational regions;

[0057] The stress direction adverse factor is used to characterize the probability of shear slip or tensile failure occurring on a structural surface after being subjected to the maximum principal stress. Its expression is:

[0058] ;

[0059] in, is the angle between the normal of the main structural surface and the direction of the maximum principal stress in the i-th calculation region;

[0060] After obtaining the local refinement trigger index for each computational region, the trigger threshold is further determined. This threshold is determined by the statistical characteristics of the local refinement trigger index of all computational regions in the current excavation stage, and its expression is as follows:

[0061] ;

[0062] in, This is the local fine-tuning trigger threshold for the current excavation stage. This represents the average value of the local refinement trigger index for the entire calculation area during the current excavation phase. The standard deviation of the local refinement trigger index for the entire calculation area during the current excavation phase. This is the threshold control coefficient;

[0063] The local sub-model is determined based on the comparison between the local refinement triggering index and the triggering threshold: If If , then the i-th computational region is identified as the local refinement trigger region; if If so, the region will remain in a macroscopic continuum calculation state.

[0064] Specifically, step S3 includes:

[0065] First, determine the spatial extent and set of boundary nodes of the local sub-model extracted in step S2; for any boundary node on the boundary of the refined local sub-model, denoted as... In the macroscopic continuum model, a macroscopic computational unit containing the spatial location of the boundary node is searched, denoted as... Where j is the boundary node number of the local refined sub-model, and k is the computational unit number in the macroscopic continuum model;

[0066] When boundary nodes Located in the macroscopic computing unit When internal, extract the node coordinates, node displacements, node stresses, and principal stress directions of the macroscopic calculation unit; when boundary nodes... When located at the common boundary of multiple macroscopic computing units, the macroscopic computing unit that is closest to the boundary node and whose inclusion or projection relationship meets the preset conditions is selected as the mapping unit;

[0067] After determining the mapping unit, the shape function interpolation method is used to establish the mapping relationship between the macroscopic computing unit nodes and the boundary nodes of the local refined sub-model; let the macroscopic computing unit be... With n nodes, the displacement vector of the r-th node is The stress tensor is The direction vector of the maximum principal stress is The shape function corresponding to this node is Then the boundary nodes of the locally refined sub-model The mapped displacement at is:

[0068] ;

[0069] in, For the boundary nodes of the local refined sub-model Mapped displacement at that point For the r-th node of the macroscopic computing unit at the boundary node The shape function value at that location, Let be the displacement vector of the r-th node in the macroscopic computing unit, and n be the number of nodes in the macroscopic computing unit.

[0070] Local refinement sub-model boundary nodes The mapped stress at that point is:

[0071] ;

[0072] in, For the boundary nodes of the local refined sub-model The mapped stress tensor at that point, Let be the stress tensor of the r-th node in the macroscopic calculation unit;

[0073] Local refinement sub-model boundary nodes The direction of the maximum principal stress at a given location is determined by the following formula:

[0074] ;

[0075] in, For the boundary nodes of the local refined sub-model The unit vector of the direction of the maximum principal stress at that location. The maximum principal stress direction vector of the r-th node in the macroscopic calculation unit is given by the denominator, which is used to normalize the direction vector obtained by interpolation.

[0076] After completing the above mapping, , and They serve as boundary nodes of the local refined sub-model. The displacement boundary conditions, stress boundary conditions, and principal stress direction parameters are defined; among them, the mapped displacement is used to constrain the boundary deformation of the local refined sub-model; the mapped stress is used to characterize the external load state of the local sub-model boundary; and the direction of the maximum principal stress is used for the evaluation of the key crack propagation potential in step S5 and the determination of the XFEM crack propagation direction in step S6.

[0077] Specifically, the DFN model construction region in step S4 is consistent with the spatial range of the local sub-model extracted in step S2, and the local coordinate system, cavern excavation boundary and principal stress direction obtained in step S3 are used as modeling references; thus, the crack geometry information in the DFN model is consistent with the stress boundary conditions, cavern spatial location and subsequent XFEM crack propagation analysis of the local sub-model.

[0078] The DFN model is established based on geological structural data and structural surface statistical parameters within the local sub-model area. The geological structural data includes fault location, fault attitude, joint set distribution, fracture development area, and spatial relationship between structural surfaces and cavern excavation boundaries. The structural surface statistical parameters include structural surface dip, dip angle, length, spacing, density, aperture, and connectivity.

[0079] When constructing the DFN model, the structural surfaces within the local sub-model are divided into deterministic structural surfaces and stochastic structural surfaces;

[0080] Deterministic structural surfaces are faults, large joints, and through fractures that have been clearly identified in field geological surveys, borehole exposure, tunnel face sketches, or 3D scanning; stochastic structural surfaces are joints and fractures generated within the local sub-model based on the statistical parameters of the structural surface group.

[0081] For deterministic structural surfaces, the corresponding structural surface geometry is directly established in the local sub-model based on its spatial location, dip, dip angle, length, and extension range; for stochastic structural surfaces, the fracture geometry is generated in the local sub-model based on the dip distribution, dip angle distribution, length distribution, spacing distribution, and density parameters of the structural surface group, and each fracture is assigned corresponding fracture parameters.

[0082] The m-th fracture in the DFN model is represented as:

[0083] ;

[0084] in, For the m-th crack, The coordinates of the center point of the crack are: The crack length is... It is prone to fractures. The angle of the fracture dip. For crack aperture, Number the structural surface group to which the crack belongs;

[0085] After the DFN model was constructed, the cracks in the DFN model were structurally identified to obtain the spatial relationship between each crack and the excavation boundary of the cavern, adjacent cracks, and the direction of local principal stress.

[0086] Calculate the minimum distance from each crack to the excavation boundary of the tunnel to determine whether the crack is located within the influence range of the tunnel boundary; identify the intersection and proximity relationships between cracks and the excavation boundary of the tunnel to determine whether cracks may participate in local tunnel wall failure; identify the intersection and potential connectivity relationships between cracks to determine whether multiple cracks may form a through-failure channel; calculate the angle between the crack normal and the direction of the maximum principal stress obtained from step S3 to determine the degree of stress disadvantage of the crack under the current excavation unloading stress environment.

[0087] For any two cracks and Calculate the minimum distance between the two. ;when Less than or equal to the preset connection distance At that time, determine the crack and cracks It has a potential connectivity relationship, and its determination formula is:

[0088] ;

[0089] in, This represents the minimum distance between the m-th fracture and the n-th fracture. This is the threshold for the distance between the fractures.

[0090] Specifically, step S5 includes:

[0091] For the m-th crack generated in step S4 Constructing a fracture propagation potential index Its expression is:

[0092] ;

[0093] in, Let m be the fracture propagation potential index of the m-th fracture. The fracture scale factor. As a factor affecting the near-cavity, As a stress driving factor, The fracture connectivity factor. As a background factor for macroscopic damage, , , , , These are the weight coefficients of the corresponding factors, and they satisfy:

[0094] ;

[0095] The crack scale factor is used to characterize the influence of the crack's own geometric scale on crack propagation, and its expression is as follows:

[0096] ;

[0097] in, Let m be the length of the m-th crack. This represents the maximum crack length among all cracks within the current local sub-model. The larger the crack length, the higher the likelihood of stress concentration at the crack tip and through-crack failure, and the larger the corresponding crack scale factor.

[0098] The proximity factor is used to characterize the spatial proximity between the fracture and the excavation boundary of the underground powerhouse cavern, and its expression is:

[0099] ;

[0100] in, Let m be the minimum distance from the m-th fissure to the excavation boundary of the cavern. The maximum control distance at which the fissure affects the local damage to the surrounding rock of the cavern; when > season =0; The closer the fissure is to the excavation boundary of the cavern, the higher the possibility that it will participate in cavern wall collapse, sidewall cracking or local through-damage.

[0101] The stress driving factor is used to characterize the driving force for shear slip or tensile propagation of the m-th crack under the local stress state obtained by mapping in step S3. Its expression is:

[0102] ;

[0103] in, Let be the tangential stress on the surface of the m-th crack. This represents the maximum tangential stress on all crack surfaces within the current local sub-model. Let be the tension normal stress on the m-th crack surface. This represents the maximum tensile normal stress on all crack surfaces within the current local sub-model. To prevent the stability coefficient from being zero, the same preset small positive number is used in the calculation of each normalized or relative change. and These are the weighting coefficients for shear-driven and tension-driven operations, respectively, and they satisfy the following:

[0104] ;

[0105] The tangential stress and tensile normal stress on the fracture surface are jointly determined by the local stress tensor obtained from step S3 and the fracture surface normal; let the unit normal vector of the m-th fracture be... The local stress tensor at the location of the crack is Then the normal stress on the crack surface is expressed as:

[0106] ;

[0107] The tangential stress on the crack surface is expressed as:

[0108] ;

[0109] in, This is the transpose of the unit normal vector of the m-th crack;

[0110] Tensioning normal stress The normal stress component that promotes crack opening and propagation; when the normal stress exhibits tensile force, take To correspond to the tensile normal stress value; when the normal stress exhibits compressive action, let =0;

[0111] The fracture connectivity factor is used to characterize the probability that the m-th fracture will form a potential through path with its adjacent fractures. Its expression is as follows:

[0112] ;

[0113] in, This represents the number of potential connected fractures within a preset connected distance range for the m-th fracture. This represents the maximum number of potentially connected fractures among all fractures in the current local sub-model; if the minimum distance between the m-th fracture and its adjacent fractures satisfies the potential connectivity criterion in step S4, then the adjacent fractures are included. ;

[0114] The macroscopic damage background factor is used to characterize the degree of damage development in the region containing the m-th crack within the macroscopic continuum model, and its expression is:

[0115] ;

[0116] in, Let m be the damage variable in the macroscopic computational region where the m-th crack is located. The maximum damage variable in the macroscopic region corresponding to the current local sub-model; when the macroscopic continuum model does not output continuous damage variables, the macroscopic damage background factor is determined according to the plastic damage factor in step S2;

[0117] After calculating the fracture propagation potential index of each fracture, a key fracture screening threshold is determined. The key fracture screening threshold is determined based on the statistical characteristics of the propagation potential indices of all fractures within the current local sub-model, and its expression is as follows:

[0118] ;

[0119] in, The threshold for screening key fractures. This represents the average value of all fracture propagation potential indices within the current local sub-model. This is the control coefficient for the crack screening threshold. This represents the standard deviation of the total crack propagation potential index within the current local sub-model.

[0120] Critical fractures are determined based on the comparison between the fracture propagation potential index and the critical fracture screening threshold; when the following conditions are met: When the m-th fracture is identified as a critical fracture, and the following condition is met: When the m-th crack is not included in the XFEM crack propagation calculation;

[0121] For the key fractures selected, the coordinates of the fracture center point, fracture length, fracture dip, fracture inclination angle, fracture aperture, fracture endpoint coordinates, and fracture surface normal vector are extracted, and the above fracture geometric parameters are mapped to the local XFEM model;

[0122] The spatial geometric position of the critical crack in the DFN model is converted into the initial crack geometry in the XFEM model coordinate system. The crack surface of the critical crack is taken as the initial crack surface of XFEM, and the end of the critical crack is taken as the crack tip position of XFEM. The local boundary displacement, boundary stress and principal stress direction obtained in step S3 are used as the boundary and loading conditions for the subsequent propagation analysis of the initial crack.

[0123] When multiple key cracks satisfy the potential connectivity relationship in step S4, the multiple key cracks are input into the XFEM model as a crack cluster, and the spatial relative positional relationship between each key crack is retained in the XFEM model for subsequent determination of whether a through-type failure channel is formed after crack propagation.

[0124] Specifically, step S6 includes:

[0125] The key cracks selected in step S5 are used as the initial cracks in the XFEM model; the spatial location, crack length, crack dip, crack dip angle, crack endpoint coordinates, and crack surface normal vector of the initial crack are provided by the DFN model constructed in step S4; the boundary displacement, boundary stress, and principal stress direction of the local XFEM model are provided by the non-matching mesh boundary mapping method in step S3.

[0126] In the crack propagation assessment process, the equivalent crack propagation driving force at the tip of the m-th critical crack is calculated. Its expression is:

[0127] ;

[0128] in, The equivalent crack propagation driving force at the tip of the m-th critical crack is... This is the Type I stress intensity factor at the tip of the m-th critical fracture, used to describe the degree to which the fracture tip continues to crack when the rock masses on both sides of the fracture open up under tensile stress. This is the Type II stress intensity factor at the tip of the m-th critical fracture, used to describe the degree to which the fracture tip continues to propagate when the rock masses on both sides of the fracture slide relative to each other along the fracture surface under shear stress. This represents the local rock mass deformation modulus. The Poisson's ratio of the local rock mass;

[0129] Equivalent crack propagation driving force Critical fracture energy of rock mass materials Compare, when the following conditions are met: When, determine that the m-th critical crack has expanded in the current calculation step; when When the m-th critical crack does not propagate in the current calculation step, it is determined that the m-th critical crack will not propagate; where, The critical fracture energy of rock mass materials;

[0130] When the critical crack meets the propagation condition, the crack propagation direction is determined based on the stress state at the crack tip; the crack propagation angle is determined using the maximum circumferential stress criterion, and the crack propagation angle satisfies:

[0131] ;

[0132] in, Let be the propagation angle of the m-th critical fracture relative to the original fracture direction;

[0133] The position of the crack tip and the crack propagation path are updated based on the crack propagation angle. If the propagating crack intersects with adjacent cracks, cavern excavation boundaries or other propagating cracks, or if the minimum distance between them is less than the preset penetration distance, it is determined that a crack penetration or potential penetration failure channel has been formed.

[0134] Specifically, step S7 includes:

[0135] The local crack propagation results output in step S6 include the cumulative crack propagation length, crack continuity, and local damage area range. For the i-th macroscopic calculation region, the propagating cracks calculated by XFEM within this region are extracted, and the number of propagating cracks located within this region is recorded as follows: ;

[0136] Based on the cumulative crack propagation length obtained in step S6, calculate the local crack propagation damage variable for the i-th macroscopic calculation region:

[0137] ;

[0138] in, Let be the local crack propagation damage variable for the i-th macroscopic computational region. Let i be the number of expanding cracks in the i-th macroscopic calculation region. Let m be the cumulative length of the m-th crack. The characteristic crack length threshold for the i-th macroscopic computational region;

[0139] When a crack penetration or potential penetration failure path exists in the i-th macroscopic calculation region as determined in step S6, the local crack propagation damage variable is corrected for penetration:

[0140] ;

[0141] in, This is the crack penetration correction factor, used to characterize the amplification effect of crack penetration on the degree of local damage; when there is no crack penetration or potential penetration failure path, =1; when there is a through crack or a potential through-path for failure. , This is the upper limit of the preset penetration correction coefficient;

[0142] Based on the local crack propagation damage variables, calculate the equivalent stiffness degradation parameters for the i-th macroscopic computational region:

[0143] ;

[0144] in, Let be the equivalent stiffness degradation parameter for the i-th macroscopic computational region. This is the lower limit of the equivalent stiffness degradation parameter, used to prevent excessive stiffness degradation in the corresponding region of the macroscopic continuum model, which would affect computational stability.

[0145] Based on the equivalent stiffness degradation parameters, the deformation modulus of the corresponding region in the macroscopic continuum model is updated using feedback:

[0146] ;

[0147] in, To provide feedback on the updated deformation modulus, To provide feedback on the deformation modulus before the update;

[0148] For macroscopic calculation regions that are not identified as local refinement trigger regions in step S2, their deformation modulus remains unchanged; for macroscopic calculation regions that have been identified as local refinement trigger regions but have not experienced crack propagation in step S6, their local crack propagation damage variable is zero, and their equivalent stiffness degradation parameter is 1. That is, the deformation modulus after feedback update is the same as the deformation modulus before update, and no stiffness degradation processing is performed.

[0149] Specifically, step S8 includes:

[0150] The results of the r-th iteration are compared with those of the (r-1)-th iteration, and the changes in displacement field, stress field, plastic damage zone, and deformation modulus are calculated respectively.

[0151] The change in displacement field is expressed as:

[0152] ;

[0153] in, Let be the change in displacement field in the r-th iteration relative to the (r-1)-th iteration. The macroscopic displacement field is obtained from the r-th iteration. The macroscopic displacement field is obtained from the (r-1)th iteration. To prevent the stability coefficient from being zero, the same preset small positive number is used in the calculation of each normalized or relative change.

[0154] The change in stress field is expressed as:

[0155] ;

[0156] in, Let be the change in stress field in the r-th iteration relative to the (r-1)-th iteration. The macroscopic stress field is obtained from the r-th iteration calculation. This represents the macroscopic stress field obtained from the (r-1)th iteration calculation;

[0157] The change in the plastic damage zone is expressed as:

[0158] ;

[0159] in, This represents the change in the plastic damage region in the r-th iteration relative to the (r-1)-th iteration. The volume of the plastic damage zone is calculated in the r-th iteration. The volume of the plastic damage zone is calculated in the (r-1)th iteration.

[0160] Introducing the deformation modulus update as a convergence criterion:

[0161] The change in deformation modulus is expressed as:

[0162] ;

[0163] in, Let be the deformation modulus update amount in the r-th iteration relative to the (r-1)-th iteration. The deformation modulus field of the macroscopic model after the r-th iteration update. The deformation modulus field of the macroscopic model after the (r-1)th iteration update;

[0164] When both conditions are met: , , , When the multi-scale coupled computation satisfies the convergence condition, it is determined that the computation is successful.

[0165] in, This is the convergence threshold of the displacement field. The stress field convergence threshold, This represents the convergence threshold of the plastic damage zone. Update the convergence threshold for deformation modulus;

[0166] If the convergence condition is not met, return to step S2, recalculate the local fine-tuning trigger index based on the updated macroscopic continuum model, re-identify the local fine-tuning trigger region, and continue executing steps S3 to S8 until the convergence condition is met. If the convergence condition is met, output the final stability evaluation result and the local crack failure prediction result. The final stability evaluation result includes the overall stress field, overall displacement field, distribution of plastic damage zone, deformation modulus degradation region, and overall stability state. The local crack failure prediction result includes the key crack propagation path, crack penetration state, potential failure channel location, and distribution of local dangerous areas.

[0167] Preferably, the macroscopic continuum model uses FLAC3D to calculate the overall stress field, displacement field, and plastic damage zone; the DFN model uses 3DEC to generate the crack spatial distribution and identify the crack spatial relationships; the XFEM model uses ABAQUS to calculate the crack evolution process; the macroscopic continuum model, DFN model, and XFEM model are coupled in stages through a unified coordinate system, boundary condition mapping, crack geometric parameter transfer, and equivalent stiffness degradation parameter feedback.

[0168] The beneficial effects of this invention are as follows:

[0169] (1) This invention constructs a local refined triggering index that includes displacement anomaly, stress concentration, plastic damage and structural surface sensitivity, and quantitatively evaluates the computational region in the macroscopic continuum model. It can comprehensively consider the macroscopic mechanical response caused by excavation unloading and the control effect of structural surfaces such as faults, joints and cracks, and avoid the identification bias caused by relying solely on human experience, a single displacement threshold or a single stress threshold to determine the local analysis region.

[0170] (2) This invention first obtains the overall stress field, displacement field and plastic damage zone of the underground powerhouse cavern group through a macroscopic continuum model, and then only carries out DFN model construction and XFEM crack propagation analysis for local areas that meet the triggering conditions. This avoids the problem of excessive calculation scale caused by fine modeling across the entire domain, and improves the accuracy of local crack failure analysis while ensuring overall calculation efficiency.

[0171] (3) This invention constructs a DFN model within a locally refined triggering region, transforming faults, joints, and fractures, which are equivalently represented in the macroscopic continuum model, into an explicit structural surface network. Based on fracture scale, proximity to the cavity, stress, connectivity characteristics, and macroscopic damage, the risk of fracture propagation and penetration is evaluated, and key fractures are selected as XFEM analysis objects. This reduces the amount of detailed calculation required for low-risk fractures and improves the targeting and computational efficiency of local fracture damage analysis.

[0172] (4) The present invention adopts the non-matching mesh boundary mapping method to map the boundary displacement, boundary stress and principal stress direction spatial interpolation obtained by the macroscopic continuum model to the boundary of the local refined sub-model, which solves the problem of boundary condition transfer caused by non-coincident nodes and inconsistent scales between the macroscopic coarse-scale mesh and the local fine-scale mesh, so that the local DFN-XFEM analysis can inherit the real mechanical environment after the overall excavation and unloading.

[0173] (5) The present invention embeds the key cracks obtained by screening into the XFEM model as the initial cracks, and determines whether the cracks are propagating based on the crack propagation driving force and critical fracture energy. It can obtain the crack initiation location, propagation direction, propagation length, penetration state and local damage area range, make up for the shortcomings of the traditional macroscopic continuum model in directly describing the local crack evolution path and crack penetration process, and provide a calculation basis for subsequent equivalent stiffness degradation feedback.

[0174] (6) This invention converts the crack propagation length, crack penetration state, and local damage area obtained from XFEM calculations into equivalent stiffness degradation parameters, and feeds back to update the deformation modulus of the corresponding region in the macroscopic continuum model; then, it recalculates the overall stress field, displacement field, and plastic damage zone, and judges whether the convergence condition is met by the changes in displacement field, stress field, plastic damage zone, and deformation modulus. This allows the results of local crack propagation to participate in the overall stability evaluation of underground powerhouse cavern groups, rather than just focusing on the output of local crack paths, thus improving the accuracy, reliability, and engineering applicability of the stability evaluation. Attached Figure Description

[0175] Figure 1 This is a flowchart of the present invention;

[0176] Figure 2 This is a schematic diagram illustrating the identification of hazardous areas in the underground powerhouse cavern complex of the present invention;

[0177] Figure 3 This is a partial schematic diagram of the non-matching mesh boundary mapping between the macroscopic model and the local sub-model of the present invention; where A is the local region of the macroscopic continuum, B is the non-matching mesh boundary mapping, and C is the local refined sub-model.

[0178] Explanation of reference numerals in the attached figures:

[0179] 1-Structural surface, 2-Main plant arch seat, 3-Main plant, 4-Main plant side wall, 5-Surrounding rock, 6-Main transformer room arch shoulder, 7-Main transformer room, 8-Connecting tunnel, 9-Boundary between main transformer room and connecting tunnel, 10-Potential risk area, 11-General surrounding rock area, 12-Local refined triggering area. Detailed Implementation

[0180] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, 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, not all, of the embodiments of the present invention. The components of the embodiments of the present invention described and shown in the accompanying drawings can generally be arranged and designed in various different configurations.

[0181] Therefore, the following detailed description of the embodiments of the invention provided in the accompanying drawings is not intended to limit the scope of the claimed invention, but merely to illustrate selected embodiments of the invention. All other embodiments obtained by those skilled in the art based on the embodiments of the invention without inventive effort are within the scope of protection of the invention.

[0182] It should be noted that similar labels and letters in the following figures indicate similar items. Therefore, once an item is defined in one figure, it does not need to be further defined and explained in subsequent figures.

[0183] In the description of this invention, it should be understood that the terms "upper," "lower," "inner," "outer," "left," "right," etc., indicate the orientation or positional relationship based on the orientation or positional relationship shown in the accompanying drawings, or the orientation or positional relationship commonly used when the product of this invention is in use, or the orientation or positional relationship commonly understood by those skilled in the art. They are only used to facilitate the description of this invention and to simplify the description, and are not intended to indicate or imply that the device or element referred to must have a specific orientation, or be constructed and operated in a specific orientation. Therefore, they should not be construed as limitations on this invention.

[0184] Furthermore, the terms "first," "second," etc., are used only to distinguish descriptions and should not be interpreted as indicating or implying relative importance.

[0185] In the description of this invention, it should also be noted that, unless otherwise explicitly specified and limited, terms such as "set" and "connection" should be interpreted broadly. For example, "connection" can be a fixed connection, a detachable connection, or an integral connection; it can be a mechanical connection or an electrical connection; it can be a direct connection or an indirect connection through an intermediate medium; it can be a connection within two components. Those skilled in the art can understand the specific meaning of the above terms in this invention according to the specific circumstances.

[0186] The specific embodiments of the present invention will now be described in detail with reference to the accompanying drawings.

[0187] like Figure 1 As shown, a multi-scale coupled adaptive numerical simulation method for underground powerhouse cavern groups includes the following steps:

[0188] S1. Establish a macroscopic continuum model of the underground powerhouse cavern group and calculate the overall stress field, displacement field and plastic damage zone under the excavation unloading condition;

[0189] In the preferred embodiment, in step S1, a macroscopic continuum model of the underground powerhouse cavern group is established based on the engineering layout of the underground powerhouse cavern group, the cavern excavation outline, the surrounding rock classification results, the geological structure distribution, the initial geostress field, the model boundary conditions, and the excavation unloading conditions.

[0190] The macroscopic continuum model divides the surrounding rock of the underground powerhouse cavern complex into continuous medium calculation units, assigning values ​​to different rock types. Each type of surrounding rock is assigned corresponding rock mechanics parameters, including deformation modulus, Poisson's ratio, cohesion, internal friction angle, tensile strength, compressive strength, and natural density. Specifically, deformation modulus and Poisson's ratio describe the deformation characteristics of the surrounding rock under excavation unloading; cohesion and internal friction angle describe the shear strength characteristics; tensile strength describes the failure threshold of the surrounding rock under tensile stress; compressive strength describes the strength characteristics of the surrounding rock under compressive stress and serves as a criterion for rock strength or a basis for parameter verification; and natural density is used to calculate the self-weight stress of the surrounding rock.

[0191] The rock mass mechanical parameters are determined by field survey data, indoor rock mechanics tests, in-situ test results and engineering analogy data, and are assigned values ​​in different zones according to the surrounding rock classification results, fault influence range and different lithological distributions.

[0192] After establishing the macroscopic continuum model, a phased excavation and unloading condition was set according to the actual excavation sequence of the underground powerhouse cavern group, and the overall stress field, displacement field, and plastic damage zone were calculated under different excavation stages. Among them, the overall stress field, displacement field, and plastic damage zone were used to identify stress concentration areas, abnormal deformation areas, and yield damage areas, respectively, and together served as input data for the calculation of the local refinement trigger index in step S2.

[0193] S2. Based on the output data of step S1, construct a local refined triggering index that includes displacement anomaly, stress concentration, plastic damage and structural surface sensitivity, identify the triggering area, and extract local sub-models;

[0194] In the preferred embodiment, in step S2, a localized refined triggering index is constructed based on the overall stress field, displacement field, plastic damage zone, and spatial distribution information of the structural surface output in step S1.

[0195] like Figure 2 As shown, the underground powerhouse cavern complex includes the main powerhouse 3, the main transformer room 7, and the connecting cave 8, all of which are located in the surrounding rock 5. Structural surfaces 1 are developed within the surrounding rock 5, including faults, joints, and fissures. Figure 2 In the diagram, the gray area represents the general surrounding rock area 11, the yellow area represents the potential risk area 10, and the red area represents the local refined triggering area 12.

[0196] In step S2, based on the overall stress field, displacement field, and plastic damage zone obtained in step S1, and combined with the spatial distribution of structural surface 1, local refined triggering indices are calculated for the main plant arch shoulder 2, the main plant side wall 4, the junction of the connecting hole 8 and the main plant 3, the main transformer room arch shoulder 6, and the junction of the main transformer room and the connecting hole 9. When the local refined triggering index of the corresponding region is greater than or equal to the preset triggering threshold, it is identified as a local refined triggering region 12 and used as the object for subsequent local sub-model extraction and refined analysis.

[0197] The local refinement trigger index is used to determine whether the computational domain in the macroscopic continuum model needs to enter the local DFN-XFEM refinement analysis.

[0198] Specifically, the macroscopic continuum model of the underground powerhouse cavern complex is divided into several computational regions, and a local refinement triggering index is calculated for each computational region.

[0199] The local refinement triggering index is obtained by coupling the macroscopic mechanical response index and the structural surface sensitivity index, and its expression is:

[0200] ;

[0201] in, Let be the local refinement trigger index for the i-th computational region. Let be the macroscopic mechanical response index of the i-th computational region. Let i be the structural surface sensitivity index for the i-th computational region. This is the structural surface sensitivity amplification factor.

[0202] The macroscopic mechanical response index is used to characterize the deformation anomalies, stress concentrations, and plastic damage levels in the computational domain under excavation unloading. Its expression is:

[0203] ;

[0204] in, Let be the displacement anomaly factor for the i-th computational region. Let be the stress concentration factor for the i-th computational region. Let be the plastic damage factor for the i-th computational region. , , These are the weighting coefficients for the displacement anomaly factor, stress concentration factor, and plastic damage factor, respectively, and they satisfy the following:

[0205] ;

[0206] The displacement anomaly factor is determined by both the displacement amplitude and the displacement gradient, and its expression is as follows:

[0207] ;

[0208] in, Let i be the displacement value of the i-th calculation region. This represents the maximum displacement value across all computational domains in the macroscopic continuum model. Let be the displacement gradient of the i-th computational region. This represents the maximum displacement gradient across all computational domains in the macroscopic continuum model. and These are the weighting coefficients for the displacement magnitude and the displacement gradient, respectively, and they satisfy:

[0209] ;

[0210] The stress concentration factor is determined by both the maximum principal stress and the deviatoric stress, and its expression is:

[0211] ;

[0212] in, The maximum principal stress in the i-th computational region is... This represents the maximum principal stress across all computational domains in the macroscopic continuum model. For the deviatoric stress in the i-th calculation region, This represents the maximum deviatoric stress across all computational domains in the macroscopic continuum model. and These are the weighting coefficients for the maximum principal stress and the deviatoric stress, respectively, and they satisfy the following:

[0213] ;

[0214] The plastic damage factor is determined by the damage variable output from the macroscopic continuum model, and its expression is:

[0215] ;

[0216] in, Let i be the damage variable for the i-th computational region. This represents the maximum damage variable across all computational domains in the macroscopic continuum model.

[0217] When the macroscopic continuum model does not output continuous damage variables, the plastic damage factor is determined by the plastic state of the computational domain, and its expression is:

[0218] ;

[0219] in, This is the set of computational regions determined to be in a plastic state in the macroscopic continuum model. If the i-th computational region belongs to the plastic region, then... =1; if the i-th computational region does not belong to the plastic region, then =0.

[0220] The structural surface sensitivity index is used to characterize the control effect of faults, joints, and fractures on local damage. Its expression is:

[0221] ;

[0222] in, Let be the structural surface density factor of the i-th computational region. Let be the near-hole distance factor for the i-th computational region. Let be the intersection factor of the structural surfaces in the i-th computational region. Let be the stress direction unfavorable factor for the i-th computational region. , , , These are the weighting coefficients for the structural surface density factor, near-hole distance factor, structural surface intersection factor, and stress direction unfavorable factor, respectively, and they satisfy the following:

[0223] ;

[0224] The structural surface density factor is used to characterize the development of faults, joints, and fractures within the computational region, and its expression is:

[0225] ;

[0226] in, Let be the surface density of the structure within the i-th computational region. This represents the maximum structural surface density across all computational regions.

[0227] The proximity factor is used to characterize the spatial proximity between the structural surface and the excavation boundary of the underground powerhouse cavern. Its expression is as follows:

[0228] ;

[0229] in, Let be the minimum distance from the main structural surface to the excavation boundary of the cavern within the i-th calculation region. This represents the maximum control distance at which the structural plane affects the stability of the surrounding rock of the cavern. When Greater than season =0.

[0230] The structural plane intersection factor is used to characterize the degree to which multiple sets of structural planes intersect within the computational region to form potential block boundaries or through-failure pathways. Its expression is:

[0231] ;

[0232] in, Let i be the number of structural surfaces intersecting within the i-th computational region. This represents the maximum number of structural plane intersections across all computational regions.

[0233] The stress direction adverse factor is used to characterize the probability of shear slip or tensile failure occurring on a structural surface after being subjected to the maximum principal stress. Its expression is:

[0234] ;

[0235] in, is the angle between the normal of the main structural surface and the direction of the maximum principal stress in the i-th calculation region.

[0236] After obtaining the local refinement trigger index for each computational region, the trigger threshold is further determined. This threshold is determined by the statistical characteristics of the local refinement trigger index of all computational regions in the current excavation stage, and its expression is:

[0237] ;

[0238] in, This is the local fine-tuning trigger threshold for the current excavation stage. This represents the average value of the local refinement trigger index for the entire calculation area during the current excavation phase. The standard deviation of the local refinement trigger index for the entire calculation area during the current excavation phase. This is the threshold control coefficient.

[0239] The local sub-model is determined based on the comparison between the local refinement triggering index and the triggering threshold: If If , then the i-th computational region is identified as the local refinement trigger region; if If so, the region will remain in a macroscopic continuum calculation state.

[0240] S3. Based on the local sub-model extracted in step S2, the non-matching mesh boundary mapping method is used to map the boundary displacement, boundary stress and principal stress direction spatial interpolation obtained from the macro model to the boundary of the local refined sub-model.

[0241] like Figure 3 As shown, a coarse mesh is used for the overall mechanical response calculation of the local region of the macroscopic continuum, while a fine mesh is used for local fine analysis in the local refined sub-model. Due to the differences in mesh scale, non-coincident node positions, and inconsistent number of boundary nodes between the coarse and fine meshes, a non-matching mesh boundary mapping method is used to transfer boundary conditions. Figure 3 In the diagram, solid black nodes represent macroscopic model boundary nodes, hollow white nodes represent local sub-model boundary nodes, and dashed lines represent the shape function interpolation relationship between them. Through shape function interpolation, the boundary displacements, boundary stresses, and principal stress directions calculated from the macroscopic model are mapped to the boundaries of the refined local sub-model.

[0242] In the preferred embodiment, in step S3, the spatial extent and set of boundary nodes of the local sub-model extracted in step S2 are first determined. For any boundary node on the boundary of the refined local sub-model, denoted as... In the macroscopic continuum model, a macroscopic computational unit containing the spatial location of the boundary node is searched, denoted as... Where j is the boundary node number of the local refined sub-model, and k is the computational unit number in the macroscopic continuum model.

[0243] When boundary nodes Located in the macroscopic computing unit When internal, extract the node coordinates, node displacements, node stresses, and principal stress directions of the macroscopic calculation unit; when boundary nodes... When located at the common boundary of multiple macroscopic computing units, the macroscopic computing unit that is closest to the boundary node and whose inclusion or projection relationship meets the preset conditions is selected as the mapping unit.

[0244] After determining the mapping unit, shape function interpolation is used to establish the mapping relationship between the macroscopic computing unit nodes and the boundary nodes of the local refined sub-model. Let the macroscopic computing unit be... With n nodes, the displacement vector of the r-th node is The stress tensor is The direction vector of the maximum principal stress is The shape function corresponding to this node is Then the boundary nodes of the locally refined sub-model The mapped displacement at is:

[0245] ;

[0246] in, For the boundary nodes of the local refined sub-model Mapped displacement at that point For the r-th node of the macroscopic computing unit at the boundary node The shape function value at that location, Let be the displacement vector of the r-th node in the macroscopic computing unit, and n be the number of nodes in the macroscopic computing unit.

[0247] Local refinement sub-model boundary nodes The mapped stress at that point is:

[0248] ;

[0249] in, For the boundary nodes of the local refined sub-model The mapped stress tensor at that point, Let be the stress tensor of the r-th node in the macroscopic calculation unit.

[0250] Local refinement sub-model boundary nodes The direction of the maximum principal stress at a given location is determined by the following formula:

[0251] ;

[0252] in, For the boundary nodes of the local refined sub-model The unit vector of the direction of the maximum principal stress at that location. The maximum principal stress direction vector of the r-th node in the macroscopic calculation unit is given, and the denominator is used to normalize the direction vector obtained by interpolation.

[0253] After completing the above mapping, , and They serve as boundary nodes of the local refined sub-model. The displacement boundary conditions, stress boundary conditions, and principal stress direction parameters.

[0254] Among them, the mapped displacement is used to constrain the boundary deformation of the local refined sub-model; the mapped stress is used to characterize the external load state of the boundary of the local sub-model; and the direction of the maximum principal stress is used for the evaluation of the key crack propagation potential in subsequent step S5 and the determination of the XFEM crack propagation direction in step S6.

[0255] S4. Construct a DFN model within the local sub-model extracted in step S2 to explicitly characterize the spatial distribution of faults, joints, and fractures;

[0256] In the preferred embodiment, in step S4, a DFN model is constructed within the local sub-model extracted in step S2 to explicitly characterize the spatial distribution of faults, joints, and fractures within the local refined triggering region, and to provide fracture geometric parameters and structural surface spatial relationship parameters for the calculation of the fracture propagation potential evaluation function in step S5.

[0257] Specifically, the DFN model construction region in step S4 is consistent with the spatial range of the local sub-model extracted in step S2, and the local coordinate system, cavern excavation boundary, and principal stress direction obtained in step S3 are used as modeling references. This ensures that the crack geometry information in the DFN model is consistent with the stress boundary conditions, cavern spatial location, and subsequent XFEM crack propagation analysis of the local sub-model.

[0258] The DFN model is established based on geological structural data and structural surface statistical parameters within the local sub-model area. The geological structural data includes fault location, fault attitude, joint set distribution, fracture development area, and the spatial relationship between structural surfaces and cavern excavation boundaries; the structural surface statistical parameters include structural surface dip, dip angle, length, spacing, density, aperture, and connectivity.

[0259] When constructing the DFN model, the structural surfaces within the local sub-model are divided into deterministic structural surfaces and stochastic structural surfaces.

[0260] Deterministic structural surfaces are faults, large joints, and through fractures that have been clearly identified in field geological surveys, borehole exposures, face sketches, or 3D scans; stochastic structural surfaces are joints and fractures generated within the local sub-model based on the statistical parameters of the structural surface group.

[0261] For deterministic structural surfaces, the corresponding structural surface geometry is directly established in the local sub-model based on its spatial location, dip direction, dip angle, length, and extent. For stochastic structural surfaces, fracture geometry is generated in the local sub-model based on the dip direction distribution, dip angle distribution, length distribution, spacing distribution, and density parameters of the structural surface group, and each fracture is assigned corresponding fracture parameters.

[0262] The m-th fracture in the DFN model is represented as:

[0263] ;

[0264] in, For the m-th crack, The coordinates of the center point of the crack are: The crack length is... It is prone to fractures. The angle of the fracture dip. For crack aperture, The structural surface group to which the crack belongs is numbered.

[0265] After the DFN model was constructed, the cracks in the DFN model were structurally identified to obtain the spatial relationship between each crack and the excavation boundary of the cavern, adjacent cracks, and the direction of local principal stress.

[0266] Specifically, the minimum distance from each crack to the excavation boundary of the tunnel is calculated to determine whether the crack is located within the influence range of the tunnel boundary; the intersection and proximity relationships between cracks and the excavation boundary of the tunnel are identified to determine whether cracks may participate in local tunnel wall failure; the intersection and potential connectivity relationships between cracks are identified to determine whether multiple cracks may form a through-damage channel; and the angle between the crack normal and the direction of the maximum principal stress obtained by mapping in step S3 is calculated to determine the degree of stress disadvantage of the crack under the current excavation unloading stress environment.

[0267] For any two cracks and Calculate the minimum distance between the two. .when Less than or equal to the preset connection distance At that time, determine the crack and cracks It has a potential connectivity relationship, and its determination formula is:

[0268] ;

[0269] in, This represents the minimum distance between the m-th fracture and the n-th fracture. This is the threshold for the distance between the fractures.

[0270] S5. Construct a crack propagation potential evaluation function, select key cracks from the DFN model generated in step S4, and identify cracks with crack propagation potential index higher than a preset threshold as key cracks, and embed them into the XFEM model as initial cracks.

[0271] In the preferred embodiment, in step S5, based on the DFN model constructed in step S4 and the local stress boundary conditions obtained by mapping in step S3, a crack propagation potential evaluation function is constructed, key cracks are screened from the crack set in the DFN model, and key cracks that meet the preset screening conditions are embedded into the XFEM model as initial cracks.

[0272] The crack propagation potential evaluation function is used to characterize the probability that a single crack in the DFN model will initiate, propagate, or connect with adjacent cracks under the current excavation unloading stress environment. Compared with directly inputting all DFN cracks into the XFEM model, this step reduces the scale of local fine-grained calculations and improves the relevance and computational efficiency of XFEM crack propagation analysis by quantitatively screening key cracks.

[0273] Specifically, for the m-th crack generated in step S4 Constructing a fracture propagation potential index Its expression is:

[0274] ;

[0275] in, Let m be the fracture propagation potential index of the m-th fracture. The fracture scale factor. As a factor affecting the near-cavity, As a stress driving factor, The fracture connectivity factor. As a background factor for macroscopic damage, , , , , These are the weight coefficients of the corresponding factors, and they satisfy:

[0276] ;

[0277] The crack scale factor is used to characterize the influence of the crack's own geometric scale on crack propagation, and its expression is as follows:

[0278] ;

[0279] in, Let m be the length of the m-th crack. This represents the maximum crack length among all cracks within the current local sub-model. A larger crack length increases the likelihood of stress concentration at the crack tip and through-crack failure, and consequently, the larger the crack scale factor.

[0280] The proximity factor is used to characterize the spatial proximity between the fracture and the excavation boundary of the underground powerhouse cavern, and its expression is:

[0281] ;

[0282] in, Let m be the minimum distance from the m-th fissure to the excavation boundary of the cavern. This represents the maximum control distance at which the fissure affects the localized damage to the surrounding rock of the cavern. When > season =0. The closer the fissure is to the excavation boundary of the tunnel, the higher the likelihood that it will contribute to the collapse of the tunnel wall, cracking of the sidewall, or localized through-damage.

[0283] The stress driving factor is used to characterize the driving force for shear slip or tensile propagation of the m-th crack under the local stress state obtained by mapping in step S3. Its expression is:

[0284] ;

[0285] in, Let be the tangential stress on the surface of the m-th crack. This represents the maximum tangential stress on all crack surfaces within the current local sub-model. Let be the tension normal stress on the m-th crack surface. This represents the maximum tensile normal stress on all crack surfaces within the current local sub-model. To prevent the stability coefficient from being zero, the same preset small positive number is used in the calculation of each normalized or relative change. and These are the weighting coefficients for shear-driven and tension-driven operations, respectively, and they satisfy the following:

[0286] ;

[0287] The tangential stress and tensile normal stress on the fracture surface are jointly determined by the local stress tensor obtained from step S3 and the fracture surface normal. Let the unit normal vector of the m-th fracture be... The local stress tensor at the location of the crack is Then the normal stress on the crack surface is expressed as:

[0288] ;

[0289] The tangential stress on the crack surface is expressed as:

[0290] ;

[0291] in, It is the transpose of the unit normal vector of the m-th crack.

[0292] Tensioning normal stress This refers to the normal stress component that promotes crack opening and propagation. When the normal stress exhibits a tensile effect, it is taken as... To correspond to the tensile normal stress value; when the normal stress exhibits compressive action, let =0.

[0293] The fracture connectivity factor is used to characterize the probability that the m-th fracture will form a potential through path with its adjacent fractures. Its expression is as follows:

[0294] ;

[0295] in, This represents the number of potential connected fractures within a preset connected distance range for the m-th fracture. This represents the maximum number of potentially connected fractures among all fractures within the current local sub-model. If the minimum distance between the m-th fracture and its adjacent fractures satisfies the potential connectivity criterion in step S4, then the adjacent fractures are included. .

[0296] The macroscopic damage background factor is used to characterize the degree of damage development in the region containing the m-th crack within the macroscopic continuum model, and its expression is:

[0297] ;

[0298] in, Let m be the damage variable in the macroscopic computational region where the m-th crack is located. This represents the maximum damage variable in the macroscopic region corresponding to the current local sub-model. When the macroscopic continuum model does not output continuous damage variables, the macroscopic damage background factor is determined based on the plastic damage factor in step S2.

[0299] After calculating the fracture propagation potential index of each fracture, a key fracture screening threshold is determined. The key fracture screening threshold is determined based on the statistical characteristics of the propagation potential indices of all fractures within the current local sub-model, and its expression is:

[0300] ;

[0301] in, The threshold for screening key fractures. This represents the average value of all fracture propagation potential indices within the current local sub-model. This is the control coefficient for the crack screening threshold. This represents the standard deviation of the total crack propagation potential index within the current local sub-model.

[0302] Critical fractures are determined based on a comparison between the fracture propagation potential index and the critical fracture screening threshold. When the following conditions are met: When the m-th fracture is identified as a critical fracture, and the following condition is met: When the m-th crack is not included in the XFEM crack propagation calculation.

[0303] For the key fractures selected, the coordinates of the fracture center point, fracture length, fracture dip, fracture angle, fracture aperture, fracture endpoint coordinates, and fracture surface normal vector are extracted, and the above fracture geometric parameters are mapped into the local XFEM model.

[0304] Specifically, the spatial geometric position of the critical crack in the DFN model is converted into the initial crack geometry in the XFEM model coordinate system. The crack surface of the critical crack is taken as the initial crack surface of XFEM, the end of the critical crack is taken as the crack tip position of XFEM, and the local boundary displacement, boundary stress and principal stress direction obtained in step S3 are taken as the boundary and loading conditions for the subsequent propagation analysis of the initial crack.

[0305] When multiple key cracks satisfy the potential connectivity relationship in step S4, the multiple key cracks are input into the XFEM model as a crack cluster, and the spatial relative positional relationship between each key crack is retained in the XFEM model for subsequent determination of whether a through-type failure channel is formed after crack propagation.

[0306] S6. Based on the force boundary conditions of the local sub-model in step S3 (force boundary conditions refer to the three types of mechanical boundary information mapped to the local sub-model boundary in step S3 and used for XFEM crack calculation), perform numerical calculations of the XFEM crack initiation, propagation and penetration evolution process.

[0307] In the preferred embodiment, in step S6, a local XFEM crack propagation model is established based on the key cracks screened in step S5 and the local boundary conditions mapped in step S3, and numerical calculations are performed on the initiation, propagation and penetration evolution process of the key cracks.

[0308] Specifically, the key cracks selected in step S5 are used as the initial cracks in the XFEM model. The spatial location, crack length, crack dip, crack inclination angle, crack endpoint coordinates, and crack surface normal vector of the initial crack are provided by the DFN model constructed in step S4; the boundary displacement, boundary stress, and principal stress direction of the local XFEM model are provided by the non-matching mesh boundary mapping method in step S3.

[0309] In the crack propagation assessment process, the equivalent crack propagation driving force at the tip of the m-th critical crack is calculated. Its expression is:

[0310] ;

[0311] in, The equivalent crack propagation driving force at the tip of the m-th critical crack is... This is the Type I stress intensity factor at the tip of the m-th critical fracture, used to describe the degree to which the fracture tip continues to crack when the rock masses on both sides of the fracture open up under tensile stress. This is the Type II stress intensity factor at the tip of the m-th critical fracture, used to describe the degree to which the fracture tip continues to propagate when the rock masses on both sides of the fracture slide relative to each other along the fracture surface under shear stress. This represents the local rock mass deformation modulus. The local Poisson's ratio is given.

[0312] Equivalent crack propagation driving force Critical fracture energy of rock mass materials Compare, when the following conditions are met: When, determine that the m-th critical crack has expanded in the current calculation step; when When the m-th critical crack does not propagate in the current calculation step, it is determined that the m-th critical crack will not propagate. It is the critical fracture energy of rock mass materials.

[0313] When the critical crack meets the propagation condition, the crack propagation direction is determined based on the stress state at the crack tip. Preferably, the maximum circumferential stress criterion is used to determine the crack propagation angle, which satisfies the following:

[0314] ;

[0315] in, Let be the expansion angle of the m-th critical fracture relative to the original fracture direction.

[0316] The crack tip position and crack propagation path are updated based on the crack propagation angle. If the propagating crack intersects with adjacent cracks, cavern excavation boundaries, or other propagating cracks, or if the minimum distance between them is less than the preset penetration distance, it is determined that a crack penetration or potential penetration failure path has been formed.

[0317] S7. Extract the local crack propagation features from step S6, calculate the equivalent stiffness degradation parameters, and feed back the equivalent stiffness degradation parameters to update the deformation modulus of the corresponding region of the macroscopic model in step S1.

[0318] In the preferred embodiment, in step S7, based on the local crack propagation results obtained in step S6, the local crack propagation damage variable and the equivalent stiffness degradation parameter are calculated, and the equivalent stiffness degradation parameter is fed back to update the deformation modulus of the corresponding region of the macroscopic continuum model in step S1.

[0319] Specifically, the local crack propagation results output in step S6 include the cumulative crack propagation length, crack continuity, and the extent of the local damage area. For the i-th macroscopic calculation region, the propagating cracks calculated by XFEM within this region are extracted, and the number of propagating cracks located within this region is recorded as follows: .

[0320] Based on the cumulative crack propagation length obtained in step S6, calculate the local crack propagation damage variable for the i-th macroscopic calculation region:

[0321] ;

[0322] in, Let be the local crack propagation damage variable for the i-th macroscopic computational region. Let i be the number of expanding cracks in the i-th macroscopic calculation region. Let m be the cumulative length of the m-th crack. is the characteristic crack length threshold for the i-th macroscopic calculation region.

[0323] In a preferred embodiment, when a crack penetration or potential penetration failure path exists in the i-th macroscopic calculation region as determined in step S6, a penetration correction is applied to the local crack propagation damage variable:

[0324] ;

[0325] in, This is the crack penetration correction factor, used to characterize the amplification effect of crack penetration on the degree of local damage. When there is no crack penetration or a potential penetration failure path, =1; when there is a through crack or a potential through-path for failure. , This is the upper limit of the preset penetration correction coefficient.

[0326] Based on the local crack propagation damage variables, calculate the equivalent stiffness degradation parameters for the i-th macroscopic computational region:

[0327] ;

[0328] in, Let be the equivalent stiffness degradation parameter for the i-th macroscopic computational region. This is the lower limit of the equivalent stiffness degradation parameter, used to prevent excessive stiffness degradation in the corresponding region of the macroscopic continuum model, which would affect computational stability.

[0329] Based on the equivalent stiffness degradation parameters, the deformation modulus of the corresponding region in the macroscopic continuum model is updated using feedback:

[0330] ;

[0331] in, To provide feedback on the updated deformation modulus, This is to provide feedback on the deformation modulus before the update.

[0332] For macroscopic calculation regions not identified as local refinement trigger regions in step S2, their deformation modulus remains unchanged. For macroscopic calculation regions identified as local refinement trigger regions but where crack propagation did not occur in step S6, their local crack propagation damage variable is set to zero, and their equivalent stiffness degradation parameter is set to 1. That is, the deformation modulus after feedback update is the same as the deformation modulus before update, and no stiffness degradation processing is performed.

[0333] S8. Based on the updated macroscopic model from step S7, recalculate the overall stress field, displacement field, and plastic damage zone of the underground powerhouse cavern group, and determine whether the convergence condition is met. If not, return to step S2 to re-identify and iteratively calculate. If yes, output the final stability evaluation and failure prediction results.

[0334] In the preferred embodiment, in step S7, the local crack propagation results have been converted into equivalent stiffness degradation parameters and fed back to update the deformation modulus of the corresponding region in the macroscopic continuum model. Therefore, in step S8, the macroscopic continuum is recalculated using the updated deformation modulus to obtain the overall displacement field, overall stress field, and plastic damage zone distribution under the new iteration.

[0335] The results of the r-th iteration are compared with those of the (r-1)-th iteration, and the changes in displacement field, stress field, plastic damage zone, and deformation modulus are calculated respectively.

[0336] The change in displacement field is expressed as:

[0337] ;

[0338] in, Let be the change in displacement field in the r-th iteration relative to the (r-1)-th iteration. The macroscopic displacement field is obtained from the r-th iteration. The macroscopic displacement field is obtained from the (r-1)th iteration. To prevent the stability coefficient from being zero, the same preset small positive number is used in the calculation of each normalized or relative change.

[0339] The change in stress field is expressed as:

[0340] ;

[0341] in, Let be the change in stress field in the r-th iteration relative to the (r-1)-th iteration. The macroscopic stress field is obtained from the r-th iteration calculation. This represents the macroscopic stress field obtained from the (r-1)th iteration.

[0342] The change in the plastic damage zone is expressed as:

[0343] ;

[0344] in, This represents the change in the plastic damage region in the r-th iteration relative to the (r-1)-th iteration. The volume of the plastic damage zone is calculated in the r-th iteration. The volume of the plastic damage zone is calculated in the (r-1)th iteration.

[0345] Furthermore, since step S7 mainly updates the deformation modulus in the macroscopic model through the equivalent stiffness degradation parameter, the deformation modulus update amount is introduced as a convergence judgment index in step S8.

[0346] The change in deformation modulus is expressed as:

[0347] ;

[0348] in, Let be the deformation modulus update amount in the r-th iteration relative to the (r-1)-th iteration. The deformation modulus field of the macroscopic model after the r-th iteration update. Let be the deformation modulus field of the macroscopic model after the (r-1)th iteration update.

[0349] When both conditions are met: , , , When the multi-scale coupled computation satisfies the convergence condition, it is determined that the computation meets the convergence condition.

[0350] in, This is the convergence threshold of the displacement field. The stress field convergence threshold, This represents the convergence threshold of the plastic damage zone. Update the convergence threshold for deformation modulus.

[0351] If the convergence condition is not met, return to step S2, recalculate the local refinement triggering index based on the updated macroscopic continuum model, re-identify the local refinement triggering region, and continue executing steps S3 to S8 until the convergence condition is met. If the convergence condition is met, output the final stability evaluation result and the local crack failure prediction result. The final stability evaluation result includes the overall stress field, overall displacement field, distribution of the plastic damage zone, deformation modulus degradation region, and overall stability state; the local crack failure prediction result includes the key crack propagation path, crack penetration state, potential failure channel location, and distribution of local hazardous areas.

[0352] In one specific implementation,

[0353] The macroscopic continuum model uses the finite difference method, finite element method, or finite volume method to calculate the overall stress field, displacement field, and plastic damage zone. In practical implementation, it can be built using FLAC3D to calculate the overall stress field, displacement field, and plastic damage zone under excavation unloading conditions. The DFN model generates the spatial distribution of cracks and identifies the spatial relationships of cracks through a discrete crack network generation module. In practical implementation, it can be built using 3DEC to generate the spatial distribution of faults, joints, and cracks, and identify crack intersection relationships, potential connectivity relationships, and the spatial relationship between cracks and the cavern excavation boundary. The XFEM model uses the extended finite element method to calculate the crack evolution process. In practical implementation, it can be built using ABAQUS to use key cracks as initial cracks and calculate the crack initiation, propagation, and penetration evolution process. The macroscopic continuum model, DFN model, and XFEM model are coupled in stages through a unified coordinate system, boundary condition mapping, crack geometric parameter transfer, and equivalent stiffness degradation parameter feedback. The aforementioned software platforms are used to perform macroscopic overall response calculation, local crack network construction, and key crack propagation analysis, respectively.

[0354] The boundary displacements, boundary stresses, and principal stress directions output by the macroscopic continuum model can be spatially interpolated and then transferred to the local sub-model; the key crack geometry parameters output by the DFN model can be transferred to the XFEM model; the crack propagation length, penetration state, and local damage range calculated by the XFEM model can be converted into equivalent stiffness degradation parameters and used to update the deformation modulus of the corresponding region of the macroscopic continuum model.

[0355] The above are merely preferred embodiments of the present invention. It should be noted that those skilled in the art can make various improvements and modifications without departing from the technical principles of the present invention, and these improvements and modifications should also be considered within the scope of protection of the present invention.

Claims

1. A multi-scale coupled adaptive numerical simulation method for underground powerhouse cavern groups, characterized in that, Includes the following steps: S1. Establish a macroscopic continuum model of the underground powerhouse cavern group and calculate the overall stress field, displacement field and plastic damage zone under the excavation unloading condition; S2. Based on the output data of step S1, construct a local refined triggering index that includes displacement anomaly, stress concentration, plastic damage and structural surface sensitivity, identify the triggering area, and extract local sub-models; S3. Based on the local sub-model extracted in step S2, the non-matching mesh boundary mapping method is used to map the boundary displacement, boundary stress and principal stress direction spatial interpolation obtained from the macro model to the boundary of the local refined sub-model. S4. Construct a DFN model within the local sub-model extracted in step S2 to explicitly characterize the spatial distribution of faults, joints, and fractures; S5. Construct a crack propagation potential evaluation function, select key cracks from the DFN model generated in step S4, and identify cracks with crack propagation potential index higher than a preset threshold as key cracks, and embed them into the XFEM model as initial cracks. S6. Based on the stress boundary conditions of the local sub-model in step S3, perform numerical calculations of the XFEM crack initiation, propagation, and penetration evolution process, and output the local crack propagation characteristics. S7. Extract the local crack propagation features from step S6, calculate the equivalent stiffness degradation parameters, and feed back the equivalent stiffness degradation parameters to update the deformation modulus of the corresponding region of the macroscopic continuum model in step S1. S8. Based on the updated macroscopic continuum model from step S7, recalculate the overall stress field, displacement field, and plastic damage zone of the underground powerhouse cavern group, and determine whether the convergence condition is met. If not, return to step S2 to re-identify and iteratively calculate. If yes, output the final stability evaluation and failure prediction results.

2. The multi-scale coupled adaptive numerical simulation method for underground powerhouse cavern groups according to claim 1, characterized in that, In step S1, the macroscopic continuum model of the underground powerhouse cavern group is established based on the engineering layout of the underground powerhouse cavern group, the cavern excavation outline, the surrounding rock classification results, the geological structure distribution, the initial geostress field, the model boundary conditions, and the excavation unloading conditions. Based on the surrounding rock classification results, the macroscopic continuum model divides the surrounding rock of the underground powerhouse cavern group into continuous medium calculation units; and assigns values ​​to each calculation unit according to the surrounding rock category corresponding to different surrounding rock classifications; each type of surrounding rock is assigned corresponding rock mass mechanical parameters, including deformation modulus, Poisson's ratio, cohesion, internal friction angle, tensile strength, compressive strength and natural density.

3. The multi-scale coupled adaptive numerical simulation method for underground powerhouse cavern groups according to claim 2, characterized in that: Step S2 includes: The macroscopic continuum model of the underground powerhouse cavern group is divided into several computational regions, and the local refined triggering index is calculated for each computational region. The local refinement triggering index is obtained by coupling the macroscopic mechanical response index and the structural surface sensitivity index, and its expression is: ； in, Let be the local refinement trigger index for the i-th computational region. Let be the macroscopic mechanical response index of the i-th computational region. Let i be the structural surface sensitivity index for the i-th computational region. This is the structural surface sensitivity amplification factor; The macroscopic mechanical response index is used to characterize the deformation anomalies, stress concentrations, and plastic damage levels in the computational domain under excavation unloading. Its expression is: ； in, Let be the displacement anomaly factor for the i-th computational region. Let be the stress concentration factor for the i-th computational region. Let be the plastic damage factor for the i-th computational region. , , The weighting coefficients for the displacement anomaly factor, stress concentration factor, and plastic damage factor are, in order, and satisfy the following: ； The displacement anomaly factor is determined by both the displacement amplitude and the displacement gradient, and its expression is as follows: ； in, Let i be the displacement value of the i-th calculation region. This represents the maximum displacement value across all computational domains in the macroscopic continuum model. Let be the displacement gradient of the i-th computational region. This represents the maximum displacement gradient across all computational domains in the macroscopic continuum model. and These are the weighting coefficients for the displacement magnitude and the displacement gradient, respectively, and they satisfy: ； The stress concentration factor is determined by both the maximum principal stress and the deviatoric stress, and its expression is: ； in, The maximum principal stress in the i-th computational region is... This represents the maximum principal stress across all computational domains in the macroscopic continuum model. For the deviatoric stress in the i-th calculation region, This represents the maximum deviatoric stress across all computational domains in the macroscopic continuum model. and These are the weighting coefficients for the maximum principal stress and the deviatoric stress, respectively, and they satisfy the following: ； The plastic damage factor is determined by the damage variable output from the macroscopic continuum model, and its expression is: ； in, Let i be the damage variable for the i-th computational region. The maximum damage variable across all computational domains in the macroscopic continuum model; When the macroscopic continuum model does not output continuous damage variables, the plastic damage factor is determined by the plastic state of the computational domain, and its expression is: ； in, Let be the set of computational regions determined to be in a plastic state in the macroscopic continuum model; if the i-th computational region belongs to the plastic region, then =1; if the i-th computational region does not belong to the plastic region, then =0; The structural surface sensitivity index is used to characterize the control effect of faults, joints, and fractures on local damage. Its expression is: ； in, Let be the structural surface density factor of the i-th computational region. Let be the near-hole distance factor for the i-th computational region. Let be the intersection factor of the structural surfaces in the i-th computational region. Let be the stress direction unfavorable factor for the i-th computational region. , , , These are the weighting coefficients for the structural surface density factor, near-hole distance factor, structural surface intersection factor, and stress direction unfavorable factor, respectively, and they satisfy the following: ； The structural surface density factor is used to characterize the development of faults, joints, and fractures within the computational region, and its expression is: ； in, Let be the surface density of the structure within the i-th computational region. The maximum structural surface density across all computational regions; The proximity factor is used to characterize the spatial proximity between the structural surface and the excavation boundary of the underground powerhouse cavern. Its expression is as follows: ； in, Let be the minimum distance from the main structural surface to the excavation boundary of the cavern within the i-th calculation region. The maximum control distance at which the structural face affects the stability of the surrounding rock of the cavern; when Greater than season =0; The structural plane intersection factor is used to characterize the degree to which multiple sets of structural planes intersect within the computational region to form potential block boundaries or through-failure pathways. Its expression is: ； in, Let i be the number of structural surfaces intersecting within the i-th computational region. The maximum number of structural surface intersections across all computational regions; The stress direction adverse factor is used to characterize the probability of shear slip or tensile failure occurring on a structural surface after being subjected to the maximum principal stress. Its expression is: ； in, is the angle between the normal of the main structural surface and the direction of the maximum principal stress in the i-th calculation region; After obtaining the local refinement trigger index for each computational region, the trigger threshold is further determined. This threshold is determined by the statistical characteristics of the local refinement trigger index of all computational regions in the current excavation stage, and its expression is as follows: ； in, This is the local fine-tuning trigger threshold for the current excavation stage. This represents the average value of the local refinement trigger index for the entire calculation area during the current excavation phase. The standard deviation of the local refinement trigger index for the entire calculation area during the current excavation phase. This is the threshold control coefficient; The local sub-model is determined based on the comparison between the local refinement triggering index and the triggering threshold: If If , then the i-th computational region is identified as the local refinement trigger region; if If so, the region will remain in a macroscopic continuum calculation state.

4. The multi-scale coupled adaptive numerical simulation method for underground powerhouse cavern groups according to claim 3, characterized in that: Step S3 specifically includes: First, determine the spatial extent and set of boundary nodes of the local sub-model extracted in step S2; for any boundary node on the boundary of the refined local sub-model, denoted as... In the macroscopic continuum model, a macroscopic computational unit containing the spatial location of the boundary node is searched, denoted as... Where j is the boundary node number of the local refined sub-model, and k is the computational unit number in the macroscopic continuum model; When boundary nodes Located in the macroscopic computing unit When internal, extract the node coordinates, node displacements, node stresses, and principal stress directions of the macroscopic calculation unit; when boundary nodes... When located at the common boundary of multiple macroscopic computing units, the macroscopic computing unit that is closest to the boundary node and whose inclusion or projection relationship meets the preset conditions is selected as the mapping unit; After determining the mapping unit, the shape function interpolation method is used to establish the mapping relationship between the macroscopic computing unit nodes and the boundary nodes of the local refined sub-model; let the macroscopic computing unit be... With n nodes, the displacement vector of the r-th node is The stress tensor is The direction vector of the maximum principal stress is The shape function corresponding to this node is Then the boundary nodes of the locally refined sub-model The mapped displacement at that point is: ； in, For the boundary nodes of the local refined sub-model Mapped displacement at that point For the r-th node of the macroscopic computing unit at the boundary node The shape function value at that location, Let be the displacement vector of the r-th node in the macroscopic computing unit, and n be the number of nodes in the macroscopic computing unit. Local refinement sub-model boundary nodes The mapped stress at that point is: ； in, For the boundary nodes of the local refined sub-model The mapped stress tensor at that point, Let be the stress tensor of the r-th node in the macroscopic calculation unit; Local refinement sub-model boundary nodes The direction of the maximum principal stress at a given location is determined by the following formula: ; in, For the boundary nodes of the local refined sub-model The unit vector of the direction of the maximum principal stress at that location. The maximum principal stress direction vector of the r-th node in the macroscopic calculation unit is given by the denominator, which is used to normalize the direction vector obtained by interpolation. After completing the above mapping, , and They serve as boundary nodes of the local refined sub-models. The displacement boundary conditions, stress boundary conditions, and principal stress direction parameters are defined; among them, the mapped displacement is used to constrain the boundary deformation of the local refined sub-model; the mapped stress is used to characterize the external load state of the local sub-model boundary; and the direction of the maximum principal stress is used for the evaluation of the key crack propagation potential in step S5 and the determination of the XFEM crack propagation direction in step S6.

5. The multi-scale coupled adaptive numerical simulation method for underground powerhouse cavern groups according to claim 4, characterized in that: The DFN model construction region in step S4 is consistent with the spatial range of the local sub-model extracted in step S2, and the local coordinate system, cavern excavation boundary and principal stress direction obtained in step S3 are used as modeling references; thus, the crack geometry information in the DFN model is consistent with the stress boundary conditions, cavern spatial location and subsequent XFEM crack propagation analysis of the local sub-model. The DFN model is established based on geological structural data and structural surface statistical parameters within the local sub-model area. The geological structural data includes fault location, fault attitude, joint set distribution, fracture development area, and spatial relationship between structural surfaces and cavern excavation boundaries. The structural surface statistical parameters include structural surface dip, dip angle, length, spacing, density, aperture, and connectivity. When constructing the DFN model, the structural surfaces within the local sub-model are divided into deterministic structural surfaces and stochastic structural surfaces; Deterministic structural surfaces are faults, large joints, and through fractures that have been clearly identified in field geological surveys, borehole exposure, tunnel face sketches, or 3D scanning; stochastic structural surfaces are joints and fractures generated within the local sub-model based on the statistical parameters of the structural surface group. For deterministic structural surfaces, the corresponding structural surface geometry is directly established in the local sub-model based on its spatial location, dip, dip angle, length, and extension range; for stochastic structural surfaces, the fracture geometry is generated in the local sub-model based on the dip distribution, dip angle distribution, length distribution, spacing distribution, and density parameters of the structural surface group, and each fracture is assigned corresponding fracture parameters. The m-th fracture in the DFN model is represented as: ； in, For the m-th crack, The coordinates of the center point of the fracture are: The crack length is... It is prone to fractures. The angle of the fracture dip. For crack aperture, Number the structural surface group to which the crack belongs; After the DFN model was constructed, the cracks in the DFN model were structurally identified to obtain the spatial relationship between each crack and the excavation boundary of the cavern, adjacent cracks, and the direction of local principal stress. Calculate the minimum distance from each crack to the excavation boundary of the tunnel to determine whether the crack is located within the influence range of the tunnel boundary; identify the intersection and proximity relationships between cracks and the excavation boundary of the tunnel to determine whether cracks may participate in local tunnel wall failure; identify the intersection and potential connectivity relationships between cracks to determine whether multiple cracks may form a through-failure channel; calculate the angle between the crack normal and the direction of the maximum principal stress obtained from step S3 to determine the degree of stress disadvantage of the crack under the current excavation unloading stress environment. For any two cracks and Calculate the minimum distance between the two. ;when Less than or equal to the preset connection distance At that time, determine the crack and cracks It has a potential connectivity relationship, and its determination formula is: ； in, This represents the minimum distance between the m-th fracture and the n-th fracture. This is the threshold for the distance between the fractures.

6. The multi-scale coupled adaptive numerical simulation method for underground powerhouse cavern groups according to claim 5, characterized in that: Step S5 includes: For the m-th crack generated in step S4 Constructing a fracture propagation potential index Its expression is: ； in, Let m be the fracture propagation potential index of the m-th fracture. The fracture scale factor. As a factor affecting the near-cavity, As a stress driving factor, For the fracture connectivity factor, As a background factor for macroscopic damage, , , , , These are the weight coefficients of the corresponding factors, and they satisfy: ； The crack scale factor is used to characterize the influence of the crack's own geometric scale on crack propagation, and its expression is as follows: ； in, Let m be the length of the m-th crack. This represents the maximum crack length among all cracks within the current local sub-model. The larger the crack length, the higher the likelihood of stress concentration at the crack tip and through-crack failure, and the larger the corresponding crack scale factor. The proximity factor is used to characterize the spatial proximity between the fracture and the excavation boundary of the underground powerhouse cavern, and its expression is: ； in, Let m be the minimum distance from the m-th fissure to the excavation boundary of the cavern. This represents the maximum control distance at which the fissure affects the localized damage to the surrounding rock of the cavern; when > season =0; The closer the fissure is to the excavation boundary of the cavern, the higher the possibility that it will participate in cavern wall collapse, sidewall cracking or local through-damage. The stress driving factor is used to characterize the driving force for shear slip or tensile propagation of the m-th crack under the local stress state obtained by mapping in step S3. Its expression is: ； in, Let be the tangential stress on the surface of the m-th crack. This represents the maximum tangential stress on all crack surfaces within the current local sub-model. Let be the tension normal stress on the m-th crack surface. This represents the maximum tensile normal stress on all crack surfaces within the current local sub-model. To prevent the stability coefficient from being zero, the same preset small positive number is used in the calculation of each normalized or relative change. and These are the weighting coefficients for shear-driven and tension-driven operations, respectively, and they satisfy the following: ； The tangential stress and tensile normal stress on the fracture surface are jointly determined by the local stress tensor obtained from step S3 and the fracture surface normal; let the unit normal vector of the m-th fracture be... The local stress tensor at the location of the crack is Then the normal stress on the crack surface is expressed as: ； The tangential stress on the crack surface is expressed as: ； in, This is the transpose of the unit normal vector of the m-th crack; Tensioning normal stress The normal stress component that promotes crack opening and propagation; when the normal stress exhibits tensile force, take To correspond to the tensile normal stress value; when the normal stress exhibits compressive action, let =0; The fracture connectivity factor is used to characterize the probability that the m-th fracture will form a potential through path with its adjacent fractures. Its expression is as follows: ； in, This represents the number of potential connected fractures within a preset connected distance range for the m-th fracture. This represents the maximum number of potentially connected fractures among all fractures in the current local sub-model; if the minimum distance between the m-th fracture and its adjacent fractures satisfies the potential connectivity criterion in step S4, then the adjacent fractures are included. ; The macroscopic damage background factor is used to characterize the degree of damage development in the region containing the m-th crack within the macroscopic continuum model, and its expression is: ； in, Let m be the damage variable in the macroscopic computational region where the m-th crack is located. The maximum damage variable in the macroscopic region corresponding to the current local sub-model; when the macroscopic continuum model does not output continuous damage variables, the macroscopic damage background factor is determined according to the plastic damage factor in step S2; After calculating the fracture propagation potential index of each fracture, a key fracture screening threshold is determined. The key fracture screening threshold is determined based on the statistical characteristics of the propagation potential indices of all fractures within the current local sub-model, and its expression is as follows: ； in, The threshold for screening key fractures. This represents the average value of all fracture propagation potential indices within the current local sub-model. This is the control coefficient for the crack screening threshold. This represents the standard deviation of the total crack propagation potential index within the current local sub-model. Critical fractures are determined based on the comparison between the fracture propagation potential index and the critical fracture screening threshold; when the following conditions are met: When the m-th fracture is identified as a critical fracture, and the following condition is met: When the m-th crack is not included in the XFEM crack propagation calculation; For the key fractures selected, the coordinates of the fracture center point, fracture length, fracture dip, fracture inclination angle, fracture aperture, fracture endpoint coordinates, and fracture surface normal vector are extracted, and the above fracture geometric parameters are mapped to the local XFEM model; The spatial geometric position of the critical crack in the DFN model is converted into the initial crack geometry in the XFEM model coordinate system. The crack surface of the critical crack is taken as the initial crack surface of XFEM, and the end of the critical crack is taken as the crack tip position of XFEM. The local boundary displacement, boundary stress and principal stress direction obtained in step S3 are used as the boundary and loading conditions for the subsequent propagation analysis of the initial crack. When multiple key cracks satisfy the potential connectivity relationship in step S4, the multiple key cracks are input into the XFEM model as a crack cluster, and the spatial relative positional relationship between each key crack is retained in the XFEM model for subsequent determination of whether a through-type failure channel is formed after crack propagation.

7. The multi-scale coupled adaptive numerical simulation method for underground powerhouse cavern groups according to claim 6, characterized in that: Step S6 includes: The key cracks selected in step S5 are used as the initial cracks in the XFEM model; the spatial location, crack length, crack dip, crack dip angle, crack endpoint coordinates, and crack surface normal vector of the initial crack are provided by the DFN model constructed in step S4; the boundary displacement, boundary stress, and principal stress direction of the local XFEM model are provided by the non-matching mesh boundary mapping method in step S3. In the crack propagation assessment process, the equivalent crack propagation driving force at the tip of the m-th critical crack is calculated. Its expression is: ； in, The equivalent crack propagation driving force at the tip of the m-th critical crack is... This is the Type I stress intensity factor at the tip of the m-th critical fracture, used to describe the degree to which the fracture tip continues to crack when the rock masses on both sides of the fracture open up under tensile stress. This is the Type II stress intensity factor at the tip of the m-th critical fracture, used to describe the degree to which the fracture tip continues to propagate when the rock masses on both sides of the fracture slide relative to each other along the fracture surface under shear stress. This represents the local rock mass deformation modulus. The Poisson's ratio of the local rock mass; Equivalent crack propagation driving force Critical fracture energy of rock mass materials Compare, when the following conditions are met: When, determine that the m-th critical crack has expanded in the current calculation step; when When the m-th critical crack does not propagate in the current calculation step, it is determined that the m-th critical crack will not propagate; where, The critical fracture energy of rock mass materials; When the critical crack meets the propagation condition, the crack propagation direction is determined based on the stress state at the crack tip; the crack propagation angle is determined using the maximum circumferential stress criterion, and the crack propagation angle satisfies: ； in, Let be the propagation angle of the m-th critical fracture relative to the original fracture direction; The position of the crack tip and the crack propagation path are updated based on the crack propagation angle. If the propagating crack intersects with adjacent cracks, cavern excavation boundaries or other propagating cracks, or if the minimum distance between them is less than the preset penetration distance, it is determined that a crack penetration or potential penetration failure channel has been formed.

8. The multi-scale coupled adaptive numerical simulation method for underground powerhouse cavern groups according to claim 7, characterized in that: Step S7 includes: The local crack propagation results output in step S6 include the cumulative crack propagation length, crack continuity, and local damage area range. For the i-th macroscopic calculation region, the propagating cracks calculated by XFEM within this region are extracted, and the number of propagating cracks located within this region is recorded as follows: ; Based on the cumulative crack propagation length obtained in step S6, calculate the local crack propagation damage variable for the i-th macroscopic calculation region: ； in, Let be the local crack propagation damage variable for the i-th macroscopic computational region. Let i be the number of expanding cracks in the i-th macroscopic calculation region. Let m be the cumulative length of the crack. The characteristic crack length threshold for the i-th macroscopic computational region; When a crack penetration or potential penetration failure path exists in the i-th macroscopic calculation region as determined in step S6, the local crack propagation damage variable is corrected for penetration: ； in, This is the crack penetration correction factor, used to characterize the amplification effect of crack penetration on the degree of local damage; when there is no crack penetration or potential penetration failure path, =1; when there is a through crack or a potential through-path for failure. , This is the upper limit of the preset penetration correction coefficient; Based on the local crack propagation damage variables, calculate the equivalent stiffness degradation parameters for the i-th macroscopic computational region: ； in, Let be the equivalent stiffness degradation parameter for the i-th macroscopic computational region. This is the lower limit of the equivalent stiffness degradation parameter, used to prevent excessive stiffness degradation in the corresponding region of the macroscopic continuum model, which would affect computational stability. Based on the equivalent stiffness degradation parameters, the deformation modulus of the corresponding region in the macroscopic continuum model is updated using feedback: ； in, To provide feedback on the updated deformation modulus, To provide feedback on the deformation modulus before the update; For macroscopic calculation regions that are not identified as local refinement trigger regions in step S2, their deformation modulus remains unchanged; for macroscopic calculation regions that have been identified as local refinement trigger regions but have not experienced crack propagation in step S6, their local crack propagation damage variable is zero, and their equivalent stiffness degradation parameter is 1. That is, the deformation modulus after feedback update is the same as the deformation modulus before update, and no stiffness degradation processing is performed.

9. The multi-scale coupled adaptive numerical simulation method for underground powerhouse cavern groups according to claim 8, characterized in that: Step S8 includes: The results of the r-th iteration are compared with those of the (r-1)-th iteration, and the changes in displacement field, stress field, plastic damage zone, and deformation modulus are calculated respectively. The change in displacement field is expressed as: ； in, Let be the change in displacement field in the r-th iteration relative to the (r-1)-th iteration. The macroscopic displacement field is obtained from the r-th iteration. The macroscopic displacement field is obtained from the (r-1)th iteration. To prevent the stability coefficient from being zero, the same preset small positive number is used in the calculation of each normalized or relative change. The change in stress field is expressed as: ； in, Let be the change in stress field in the r-th iteration relative to the (r-1)-th iteration. The macroscopic stress field is obtained from the r-th iteration calculation. This represents the macroscopic stress field obtained from the (r-1)th iteration calculation; The change in the plastic damage zone is expressed as: ； in, This represents the change in the plastic damage region in the r-th iteration relative to the (r-1)-th iteration. The volume of the plastic damage zone is calculated in the r-th iteration. The volume of the plastic damage zone is calculated in the (r-1)th iteration. Introducing the deformation modulus update as a convergence criterion: The change in deformation modulus is expressed as: ； in, Let be the deformation modulus update amount in the r-th iteration relative to the (r-1)-th iteration. The deformation modulus field of the macroscopic model after the r-th iteration update. The deformation modulus field of the macroscopic model after the (r-1)th iteration update; When both conditions are met: , , , When the multi-scale coupled computation satisfies the convergence condition, it is determined that the computation is successful. in, The threshold for convergence of the displacement field. The stress field convergence threshold, This represents the convergence threshold of the plastic damage zone. Update the convergence threshold for deformation modulus; If the convergence condition is not met, return to step S2, recalculate the local fine-tuning trigger index based on the updated macroscopic continuum model, re-identify the local fine-tuning trigger region, and continue executing steps S3 to S8 until the convergence condition is met. If the convergence condition is met, output the final stability evaluation result and the local crack failure prediction result. The final stability evaluation result includes the overall stress field, overall displacement field, distribution of plastic damage zone, deformation modulus degradation region, and overall stability state. The local crack failure prediction result includes the key crack propagation path, crack penetration state, potential failure channel location, and distribution of local dangerous areas.

10. The multi-scale coupled adaptive numerical simulation method for underground powerhouse cavern groups according to claim 6, characterized in that: The macroscopic continuum model uses FLAC3D to calculate the overall stress field, displacement field, and plastic damage zone; the DFN model uses 3DEC to generate the spatial distribution of cracks and identify the spatial relationships of cracks; the XFEM model uses ABAQUS to calculate the crack evolution process; the macroscopic continuum model, DFN model, and XFEM model are coupled in stages through a unified coordinate system, boundary condition mapping, crack geometric parameter transfer, and equivalent stiffness degradation parameter feedback.

Images