Rock crack evolution numerical simulation method based on data analysis

Abstract

The invention discloses a rock crack evolution numerical simulation method based on data analysis, and relates to the technical field of rock mechanics engineering numerical simulation, and the method comprises the following steps: S1, obtaining microstructure data of a target rock, including a CT scanning image, an SEM image and mineral energy spectrum data; and S2, constructing a spatial distribution dynamic mechanical parameter field based on the microstructure data. According to the rock crack evolution numerical simulation method based on data analysis, a dynamic mechanical parameter field driven by a microstructure is constructed, so that the inherent defect that rock is simplified into a homogeneous material in a traditional method is overcome; microscopic characteristics such as mineral hardness gradient, grain boundary strength distribution and pore morphology are effectively quantified into an elastic modulus field, a tensile strength field and a fracture energy field which are spatially distributed, so that the numerical prediction of crack initiation and expansion truly reflects the inherent heterogeneity of the rock for the first time.

Description

Technical Field

[0001] The present invention relates to the technical field of numerical simulation of rock mechanics engineering, and in particular to a numerical simulation method of rock crack evolution based on data analysis. Background Art

[0002] In the field of rock mechanics and engineering, accurately predicting the initiation, propagation and evolution of rock cracks is crucial for assessing rock stability. Traditional numerical simulation methods, such as the finite element method, discrete element method or extended finite element method, face significant limitations in simulating the evolution of rock cracks. These methods usually assume that rock is a continuous homogeneous material or rely on simplified macroscopic mechanical parameters, failing to fully consider the complex microscopic heterogeneity inherent in rock. Rocks are essentially heterogeneous bodies composed of mineral particles, cements, natural pores and microcracks. The spatial distribution of their mineral components, differences in grain boundary strength and the randomness of micro-defects fundamentally dominate the initial position, propagation path and speed of cracks.

[0003] However, existing mainstream technologies have difficulty effectively converting real microstructural features, such as mineral phase distribution and pore geometry, into key input parameters to drive the simulation process, and are even more unable to quantify the feedback effects of these heterogeneous factors in real time during the dynamic evolution of cracks. This leads to systematic deviations between simulation results and laboratory observations or engineering practice when predicting crack network morphology, fracture energy consumption, and final failure mode, especially under complex multi-field coupling conditions. Therefore, there is an urgent need for a numerical simulation method that can deeply integrate the real microstructural information of rocks and dynamically reflect its influence on crack evolution. Summary of the Invention

[0004] To achieve the above objectives, the present invention is implemented through the following technical solutions: a numerical simulation method of rock crack evolution based on data analysis, comprising the following steps: S1: Acquire the microstructural data of the target rock, including CT scan images, SEM images, and mineral energy spectrum data; further, first obtain the microstructural multimodal data of the target rock sample, including three-dimensional tomographic images acquired by medical CT equipment, multi-field high-resolution microscopic images acquired by scanning electron microscopes, and mineral composition distribution data generated by electron probe microanalyzers. For special rock samples with low permeability or high water content, pre-processing in a vacuum environment or freezing conditions is required to maintain the original structural morphology. After all data are unified to the same coordinate system through a spatial registration algorithm, they are input into a deep convolutional neural network to segment the mineral phase boundaries and pore areas. When encountering grain adhesion or sub-resolution microcracks, the edge enhancement algorithm and three-dimensional continuity verification are activated.

[0005] S2: Based on the aforementioned mesostructural data, a spatially distributed dynamic mechanical parameter field is constructed. Furthermore, based on the segmentation results and mineral energy spectrum data, a heterogeneous parameter transfer function system is designed. This includes directly invoking standard parameters from the micromechanics experimental database for single-phase mineral regions, calculating equivalent properties weighted by elemental proportions in multiphase mineral intergrowth regions, and introducing stress concentration correction factors in porous regions based on 3D morphology analysis. The transfer process strictly adheres to mesostructural constraints, implementing an exponential smoothing algorithm at mineral phase transition boundaries. The final output is a spatially discretized dynamic mechanical parameter field, including the elastic modulus field, the tensile strength field, and the fracture energy field.

[0006] S3: The extended finite element method (EFM) was used to drive the numerical simulation of crack evolution, integrating microfracture acoustic emission monitoring data in real time during the simulation. Furthermore, the EFM was used to drive the crack evolution simulation while simultaneously incorporating microfracture signals collected by an acoustic emission sensor array. The raw waveforms were bandpass filtered and subjected to time-frequency analysis to extract the energy release rate and source coordinates, which were then synchronously correlated with the numerical simulation time step.

[0007] S4: When the deviation between the monitored data and the simulated path exceeds a set threshold, the dynamic mechanical parameter field is dynamically corrected. Furthermore, when the monitored acoustic emission energy trend continuously deviates from the simulated crack tip mechanical response, a dynamic ellipsoidal influence domain is defined centered on the crack front, and a particle filter inversion mechanism is initiated. This involves randomly generating a swarm of particles representing the fracture energy field distribution, calculating the theoretical response of each particle under an extended finite element method, updating particle weights and resampling based on Bayes' theorem, and outputting a corrected fracture energy field distribution. The inversion process is constrained by mineral phase boundaries, limiting the correction range in hard mineral regions and allowing greater freedom along the connectivity direction in porous regions.

[0008] For high-stress hard rock environments, an event-driven emergency response model is employed. This includes: acoustic emission mainshock events trigger updates to the full particle set, with the radius of the affected region expanding exponentially with the crack acceleration. For soft rock shear conditions, a progressive optimization approach is employed: an energy accumulation window is set, and a small step inversion is performed each time a threshold energy unit is reached within the window. The affected region is stretched along the potential slip surface. Fluid parameter constraints are introduced for hydraulic fracturing conditions: the fracturing fluid viscosity limits the particle search space, and the proppant coordinate set serves as the basis for topological node generation.

[0009] S5: The validity of the simulation results is determined through a dual verification mechanism of physical morphology and mathematical topology. Furthermore, a dual verification is performed after the simulation is completed. This includes: During the physical morphology verification phase, the simulated crack network is spatially aligned with the actual fracture point cloud obtained by 3D laser scanning. Normal sampling points are set along the main path to calculate the maximum distance extreme value, which must be less than the product of the characteristic mineral size and a preset coefficient. During the mathematical topology verification phase, the topological invariants of the simulated and measured crack networks are extracted. A life cycle barcode is generated using a persistent homology algorithm. Topological similarity is determined when the overlap length of the stable intervals of the main ring structure and the deviation in the number of branches meet the set criteria.

[0010] If any verification fails, the coordinates of the failure region are located, and the parameter field optimization cycle is restarted back to the corresponding simulation time step. At this time, the adaptation strategy is based on the failure type, including: geometric path deviation focuses on correcting the fracture energy field gradient, and topological structure mismatch adjusts the intensity field anisotropy. Among them, three core mechanisms work together throughout the entire process: parameter field construction pre-embeds physical boundaries, dynamic correction of response acoustic emission event characteristics, and verification results drive targeted optimization cycles to form a closed-loop control system.

[0011] Preferably, in step S2, a deep learning model is used to segment the mineral phases, pores, and microcracks in the CT scan image and SEM image to generate a microstructure distribution map. Furthermore, when there is a systematic deviation between the mineral phase segmentation results and the energy spectrum data, a local re-acquisition process is initiated: the field of view density of the scanning electron microscope is increased in the deviation area, and the electron probe point scanning frequency is simultaneously increased. If grain boundary identification is affected by insufficient contrast in the microscopic image, a chemical dye is injected to enhance mineral phase differentiation. The dye should be selected to avoid chemical reactions with rock components.

[0012] Preferably, in step S2, the construction of the dynamic mechanical parameter field includes: designing a heterogeneous parameter transfer function based on the mineral energy spectrum data and the micromechanical experimental calibration results; and mapping the mineral hardness and grain boundary strength into a spatially distributed elastic modulus field, tensile strength field and fracture energy field through the transfer function.

[0013] Preferably, the dynamic correction of step S4 includes: coupling the particle filter algorithm with the extended finite element method to establish a dynamic association between the energy release rate of the acoustic emission event and the simulated crack tip stress intensity factor; using the acoustic emission energy release rate as the observation variable and the crack tip stress intensity factor as the state variable, and inverting the local fracture energy field through particle filter weight update.

[0014] Preferably, the inversion is performed only on a local region of the crack growth front, and the update frequency is synchronized with the acoustic emission event rate.

[0015] Preferably, the physical morphology verification of step S5 includes: aligning the simulated crack network topology with the 3D laser scanning point cloud of the actual damage section, calculating the Hausdorff distance between the two, and when the Hausdorff distance is less than a preset ratio of the characteristic length, determining that the physical verification is passed.

[0016] Preferably, the mathematical topology verification includes: extracting the Betti numbers and Euler characteristic numbers of the simulated and experimental crack networks, calculating the topological similarity through persistent homology analysis, and determining that the mathematical verification is passed when the similarity meets a predetermined topological matching standard.

[0017] Preferably, the microfracture acoustic emission monitoring data is derived from acoustic emission signals collected in a laboratory or on-site, and the energy release rate and event location coordinates are extracted after time-frequency analysis.

[0018] Preferably, if either the physical morphology verification or the mathematical topology verification fails in step S5, the process automatically returns to step S4 to restart the parameter field optimization cycle.

[0019] Preferably, the dynamic mechanical parameter field construction, dynamic correction and dual verification mechanism work together in the entire crack evolution process to quantify the dynamic impact of microscopic heterogeneity. Furthermore, the initialization phase inherits the topological relationship of the mineral phase to define the correction boundary; the extended calculation period adopts event-triggered resource allocation, including: real-time correction is performed in the area associated with the acoustic emission event, and the parameter field is frozen in the non-associated area; a reverse diagnostic path is established in the verification backtracking period: when the topological invariant anomaly points to a specific mineral phase, the parameter transfer function is automatically triggered to be recalibrated. Under extreme working conditions, the degradation mode is activated, including: the fine topological analysis of deep rock bursts is closed and switched to a rapid verification based on the total energy of acoustic emissions; when continuous optimization fails, the current optimal solution is output and the risk area is marked.

[0020] The present invention provides a numerical simulation method for rock crack evolution based on data analysis. It has the following beneficial effects: This data-analysis-based numerical simulation method for rock crack evolution overcomes the inherent defect of traditional methods that simplify rocks into homogeneous materials by constructing a dynamic mechanical parameter field driven by microstructures. It effectively quantifies microscopic characteristics such as mineral hardness gradients, grain boundary strength distributions, and pore morphology into spatially distributed elastic modulus fields, tensile strength fields, and fracture energy fields, enabling the numerical prediction of crack initiation and propagation to truly reflect the inherent heterogeneity of rocks for the first time. Combining acoustic emission monitoring data with a particle filter inversion algorithm, real-time calibration of parameter fields during crack evolution is achieved, improving the accuracy of crack network morphology predictions under complex stress paths. Compared with traditional numerical methods, the crack path prediction error has been reduced, and it remains stable and reliable under extreme working conditions such as high-in-situ stress rockbursts and shear slip of layered rock masses.

[0021] This data-analysis-based numerical simulation method for rock crack evolution establishes a closed-loop quality control system through a dual verification mechanism. Physical morphology verification ensures geometric consistency between simulated cracks and real fracture surfaces, while mathematical topology verification quantifies network connectivity characteristics at the structural level. This dual-track screening minimizes false positives. This creates a full-cycle data closed loop, from microscopic characterization to macroscopic prediction. Parameter field construction embeds physical boundary constraints, corrects responses to real-time monitoring signals, and verification results drive a targeted optimization cycle, reducing the cost of manual trial and error. BRIEF DESCRIPTION OF THE DRAWINGS

[0022] Figure 1 Schematic diagram of a framework of a numerical simulation method for rock crack evolution based on data analysis according to the present invention; Figure 2 This is a control logic timing diagram of a numerical simulation method for rock crack evolution based on data analysis of the present invention. DETAILED DESCRIPTION

[0023] The following will clearly and completely describe the technical solutions in the embodiments of the present invention in conjunction with the accompanying drawings. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative efforts are within the scope of protection of the present invention.

[0024] See also Figure 1 and Figure 2 The present invention provides a technical solution: a numerical simulation method for rock crack evolution based on data analysis, comprising the following steps: S1: Acquire microstructural data of the target rock, including CT scans, SEM images, and mineral energy spectrum data. It should be noted that, during the specific implementation process, a medical CT device is first used to obtain three-dimensional tomographic images at submicron resolution for the target rock sample, ensuring that the scan layer thickness is no greater than one-third of the rock's smallest mineral grain size. For micron-scale pore structure, a scanning electron microscope (SEM) is used in backscattering mode to acquire high-resolution SEM images of multiple fields of view, with each field covering at least 10 representative mineral grains. A full-scale microstructural map is constructed using image stitching techniques. Simultaneously, an electron probe microanalyzer is used to perform mineral energy spectrum scanning, generating a mineral phase distribution map at corresponding locations in the SEM image. If carbonate-soluble minerals or clay minerals are detected, additional energy spectrum point scanning is performed to quantify the cement composition. For low-permeability shales or tight sandstones, the samples are pretreated in a vacuum environment to prevent microcrack closure. For soft rocks with a water content exceeding 5%, cryo-fixation is used to preserve the original pore morphology. After all data acquisition is completed, the CT volume data, SEM two-dimensional image sequence and energy spectrum data are unified into the same coordinate system through a feature point-based spatial registration algorithm. If the registration error exceeds 20% of the pixel size, the local area data is re-collected until the accuracy requirements are met. S2: Constructing a spatially distributed dynamic mechanical parameter field based on mesostructural data. It should be further explained that, in the specific implementation process, based on the registered mesostructural data, a deep convolutional neural network is first used to segment mineral phase boundaries and pore regions. For rock samples with significant grain size differences, a multi-scale feature fusion module is activated to ensure the continuous identification of microcracks. Based on the compositional distribution calibrated by the mineral energy spectrum data, a mapping library of mineral type and microhardness is established. If amorphous cement is detected, its interface strength properties are supplemented by nanoindentation experimental data. Subsequently, a heterogeneous parameter transfer function is designed. This includes: for hard mineral regions such as quartz and feldspar, the local elastic modulus gradient distribution is calculated based on the lattice orientation; for weak interlayers such as mica and clay, anisotropic strength attenuation characteristics are assigned based on the lamellar structural orientation in the SEM image; and for pore regions, the fracture energy parameter is dynamically reduced based on the results of three-dimensional morphological analysis. When the pore is flat and the long axis is parallel to the principal stress direction, an additional stress concentration correction factor is activated. All parameter fields are spatially discretized using voxels as units. A transition smoothing algorithm is used at the boundaries of mineral phase transitions to avoid mechanical mutations. If the difference in elastic modulus between adjacent voxels exceeds a set threshold, subvoxel interpolation and reconstruction are initiated. The final output is a dynamic mechanical parameter field that strictly corresponds to the microstructure, including but not limited to the elastic modulus field, tensile strength field, and fracture energy field. S3: An extended finite element method (EFM) is used to drive numerical simulations of crack evolution, integrating microfracture acoustic emission monitoring data in real time during the simulation process. It should be noted that, in the specific implementation, an extended finite element model is established based on a dynamic mechanical parameter field, driving crack evolution calculations under preset stress boundary conditions. The raw waveform signals collected by the AE sensor array are subjected to wavelet transform to extract characteristic frequency components, identify effective microfracture events, and locate the source coordinates. These signals are then synchronously mapped to corresponding spatial locations in the simulation domain, establishing a synchronous association between the AE event sequence and the numerical simulation time step. During crack propagation, the AE energy release rate and the simulated crack tip stress intensity factor are compared in real time. If the two deviate in the opposite direction within three consecutive event cycles, a significant deviation is identified and the degree of deviation is recorded. When the accumulated deviation exceeds a critical threshold, a dynamic correction command is triggered. For deep hard rock high-in-situ stress conditions, event clustering analysis is used to select the AE signals corresponding to the dominant crack cluster and filter out background rock noise. In the case of anisotropic cracking in layered rock, the monitoring data weight coefficient is adjusted based on the bedding plane azimuth to ensure physical consistency of the data fusion. All fusion operations are completed in independent data assimilation modules to avoid interfering with the parallel computing architecture of the core solver; S4: When the deviation between the monitored data and the simulated path exceeds a set threshold, the dynamic mechanical parameter field is dynamically corrected. It should be noted that, in specific implementation, when the cumulative deviation between the monitored acoustic emission energy release rate and the simulated crack tip stress intensity factor exceeds a critical threshold, a particle filter-driven local parameter field correction process is initiated. An ellipsoidal influence domain is defined, centered at the crack front of the previous simulation step. The current fracture energy field data for all voxels within the domain are extracted as the initial state of the particle swarm. Using the acoustic emission event energy as the observation value and the crack tip stress intensity factor as the state variable, the fracture energy field correction is inverted using a particle weight update mechanism. The inversion process strictly adheres to mineral phase boundary constraints, including limiting the fracture energy field to a narrow range in hard mineral regions and allowing greater correction freedom in weak interlayer regions. Connectivity detection is activated when encountering densely porous zones to prevent the correction from leading to unphysical mechanical isolation zones. Once the correction is complete, the fracture energy parameters of the voxels within the influence domain are immediately updated, and the associated coupling terms of the elastic modulus and tensile strength fields are simultaneously optimized, while retaining the original parameters of regions unaffected by crack propagation. For sudden splitting of hard rock under high ground stress, a post-event retrospective correction mode is used to complete parameter field iteration within three calculation steps after the acoustic emission main shock. For progressive shear failure of layered rock masses, a continuous small-step incremental update strategy is activated to ensure real-time correction. S5: The validity of the simulation results is determined through a dual verification mechanism of physical morphology and mathematical topology. It should be further explained that, in the specific implementation process, after completing the crack evolution simulation, the damaged rock sample is first subjected to 3D laser scanning to obtain real fracture point cloud data. The simulated crack network topology is then spatially aligned with the measured point cloud using a feature point matching algorithm. The physical morphology verification phase calculates the Hausdorff distance between the two. This involves setting equally spaced sampling points along the main crack path and calculating the minimum distance extreme value in the normal direction of each point. Verification is considered successful if all sampling point extreme values ​​are less than the product of the characteristic length and a preset scaling factor. The mathematical topology verification phase extracts the topological invariants of the simulated and measured crack networks. This includes generating a skeletonized structure diagram based on the crack branch points and endpoints, calculating the Betti number to describe the number of independent loop holes, and quantifying connectivity differences using Euler characteristic numbers. A lifecycle barcode of the topological features is constructed using a persistent homology algorithm. Topological similarity is determined when the overlap of the barcode's stable intervals exceeds a set standard. In the case of complex networks with multiple cracks, the system prioritizes segmentation into sub-crack systems for block verification. For tree-like bifurcated cracks generated by hydraulic fracturing, topological similarity is calculated based on the consistency of bifurcation angles. If any verification fails, the coordinates of the failed region are automatically marked, triggering a parameter field backtracking optimization instruction.

[0025] In step S2, a deep learning model is used to segment the mineral phases, pores, and microcracks in the CT scan and SEM images to generate a microstructural distribution map. It should be noted that, in the specific implementation, the spatially registered CT volume data and SEM image sequences are fed into a pre-trained deep convolutional neural network, which employs an encoder-decoder architecture and embeds a multi-scale feature fusion module. The CT data are first sliced ​​along three orthogonal planes to identify the three-dimensional topological continuity of mineral phase boundaries and pore space. Microcrack morphology and grain distribution characteristics in the SEM images are simultaneously analyzed, and a sub-pixel edge enhancement algorithm is activated when grain sizes below the scanning resolution are encountered. For hard rocks such as granite, the weight of quartz feldspar grain boundary segmentation is prioritized. For porous rocks such as sandstone, the pore connectivity detection branch is activated to distinguish isolated pores from seepage channels. If mineral phase adhesion is observed in the output segmentation map, physical constraints are applied based on the mineral composition distribution of the energy spectrum data. These constraints include forced segmentation in regions with significant energy spectrum differences between adjacent minerals, while retaining morphological smoothing in gradient transition zones. The final generated microstructure distribution map needs to pass the grain boundary closure check, and the shortest distance interpolation is used to repair the crack paths with breakpoints to ensure the topological integrity of the subsequent parameter field construction.

[0026] In step S2, the construction of the dynamic mechanical parameter field includes: designing a heterogeneous parameter transfer function based on the mineral energy spectrum data and the micromechanical experimental calibration results; mapping the mineral hardness and grain boundary strength into a spatially distributed elastic modulus field, tensile strength field and fracture energy field through the transfer function. It should be further explained that in the specific implementation process, a transfer function system from mineral type to mechanical parameters is constructed based on the spatial distribution of components calibrated by the mineral energy spectrum data and the micromechanical experimental database. For the single-phase mineral-dominated area, the stress-strain curve characteristic values ​​of the calibration experiment are directly called to generate basic parameters; for the multi-phase mineral coexistence area, the equivalent properties are calculated by weighting according to the proportion of energy spectrum elements. If a calcite-dolomite gradient transition zone is detected, the lattice mismatch correction factor is activated. The transfer function design adheres to microstructural constraints. These include assigning a positively correlated elastic modulus distribution to hard mineral regions such as quartz based on CT grayscale gradients, and assigning anisotropic tensile strength attenuation to layered mica schist structures based on SEM-identified schistosity dips. In pore regions, based on 3D topography analysis, a stress concentration factor compensation mechanism is activated when the pore aspect ratio exceeds a critical value. Physical consistency checks are performed on all parameter fields after generation. This includes reconstructing transition zones using an exponential smoothing algorithm based on grain boundary distances if the elastic modulus jump between adjacent voxels exceeds the mineral phase transition threshold. Furthermore, the fracture energy field is constructed by correlating pore connectivity paths to prevent the formation of non-physical crack barriers in isolated high-energy regions.

[0027] The dynamic correction in step S4 involves coupling a particle filter algorithm with the extended finite element method to establish a dynamic correlation between the energy release rate of acoustic emission events and the simulated crack tip stress intensity factor. Using the acoustic emission energy release rate as the observation variable and the crack tip stress intensity factor as the state variable, the local fracture energy field is inverted through particle filter weight updating. It should be further explained that in the specific implementation, during the crack propagation simulation, the energy release rate of microfracture events captured by the acoustic emission sensor is used as the observation sequence, and the crack tip stress intensity factor is used as the state variable to construct a dynamic correlation model for the particle filter system. An ellipsoidal influence domain is established with the current crack front as the origin, and a randomly generated particle swarm characterizes the distribution of fracture energy field parameters within the domain. Each particle represents a set of fracture energy field correction hypotheses, and the theoretical response of the crack tip stress intensity factor under these hypotheses is calculated using the extended finite element method. When the cumulative residual between the acoustic emission energy observation and the particle prediction exceeds an adaptive threshold, a Bayesian inversion mechanism is triggered, which includes updating particle weights based on the observation likelihood and resampling high-weight particles to generate a corrected fracture energy field distribution. The inversion process imposes mineral phase constraints, including: fracture energy corrections in hard mineral regions must not exceed the fluctuation range of micromechanical experimental calibration values, and in the porous region, greater degrees of freedom are allowed for corrections in the direction of the connected pore network. If the particle weight entropy fails to converge after three consecutive inversions, it is determined to be a sudden change in rock mass heterogeneity, and the influence domain expansion process is initiated.

[0028] The inversion is performed only on a localized region at the crack propagation front, and the update frequency is synchronized with the acoustic emission event rate. It should be further explained that, in the specific implementation, the stress concentration area at the crack tip is defined as a dynamic influence domain of an ellipsoid, with its principal axis oriented along the current crack propagation vector and its semi-axis length adaptively adjusted by the spatial distribution density of the acoustic emission events. When the acoustic emission sensor array captures a valid microfracture event, the extended finite element calculation process is immediately frozen and a local inversion procedure is initiated, including: the influence domain range is dynamically expanded and contracted with the crack velocity, the number of voxels in the domain is compressed to the typical mineral grain coverage for brittle cracking in hard rock, and creep failure in soft rock is extended to the potential shear band width.

[0029] The update frequency is strictly tied to the characteristics of the acoustic emission event chain. For mainshock-aftershock rupture sequences, inversion is completed within three simulation time steps after the mainshock event. Successive microseismic swarms are triggered in batches when the integrated energy of the event reaches a single correction threshold. In the case of rockbursts with multiple crack branches under high in-situ stress, an independent influence domain and particle swarm are assigned to each active crack tip, and parallel inversion is achieved using the event source location coordinates. After inversion, only the fracture energy parameters of voxels within the influence domain are updated, and the correction timestamps are recorded for backtracking topology verification.

[0030] The physical morphology verification of step S5 includes: aligning the simulated crack network topology with the 3D laser scanning point cloud of the actual damage section, calculating the Hausdorff distance between the two, and judging that the physical verification is passed when the Hausdorff distance is less than 5% of the preset ratio of the characteristic length. It should be further explained that in the specific implementation process, after completing the crack evolution simulation, a blue light 3D scanner is used to obtain the surface point cloud data of the damaged sample with a resolution not less than the average grain size of the rock, and the simulated crack topology is mapped to the measured coordinate system through the characteristic mineral corner matching algorithm. Equally spaced verification points are set along the main path of the simulated crack, and the distance to the nearest measured point is searched at each point along the normal direction of the section, and the maximum distance extreme value of all verification points is recorded as the Hausdorff distance reference value. When the rock sample contains multiple secondary cracks, the main crack path is aligned first and then the branch network is processed; if a melt fracture is generated in a high temperature and high pressure environment, the anti-thermal deformation calibration program is activated to correct the point cloud deformation error. Distance determination uses a relative standard, including using the rock's characteristic mineral dimensions as a reference length. Physical morphology verification is considered passed when the maximum distance is less than the product of this reference length and a preset coefficient. Brittle fractures in hard rock are also verified for curvature consistency along the crystalline fracture path, while plastic fractures in soft rock focus on shear band width matching.

[0031] High-temperature melt fractures require thermal deformation compensation, including correcting point cloud coordinates based on the rock sample's thermal expansion coefficient and filtering out transient topological features below the glass transition temperature. Anisotropic rock mass verification utilizes directional weighting, including increasing the distance tolerance for sampling points in the bedding direction to twice that of the vertical direction and compressing scale coordinates along the bedding plane using topological barcodes. In-situ engineering fractures utilize multi-source data fusion, including LiDAR point cloud compensation for obscured areas and core sampling data to assist in verifying hidden crack branching.

[0032] Mathematical topological verification involves extracting the Betti numbers and Euler characteristic numbers of the simulated and experimental crack networks, calculating topological similarity through persistent homology analysis, and passing the mathematical verification when the similarity meets the predetermined topological matching criteria, i.e., when the similarity is greater than 90%. It should be further explained that, in the specific implementation process, the simulated and measured crack networks, after physical registration, are first skeletonized to generate a simplified topological structure diagram. Key topological invariants are extracted, including calculating the zero-dimensional Betti number to characterize the number of independent crack branches and the one-dimensional Betti number to describe the total amount of closed loop holes; the Euler characteristic number is also simultaneously solved to quantify the global connectivity of the network. A lifecycle barcode of the topological features is constructed through a persistent homology algorithm, including progressively filtering the crack network with a distance function, recording the scale range of each ring structure from formation to disappearance, and generating a multidimensional barcode dataset. The similarity determination phase compares the distribution of stable intervals between the simulated and measured barcodes. Topological similarity is determined when the overlap length of the birth-death scale intervals of the main ring holes exceeds a critical ratio and the deviation in the number of independent branches is less than the allowable error. For tree-like bifurcations formed by hydraulic fracturing, a histogram of the bifurcation order distribution is additionally calculated. For cracks with a cyclic loading history, the topological evolution trajectory is verified by separating subnetworks at each loading stage.

[0033] Microfracture acoustic emission monitoring data is derived from acoustic emission signals collected in the laboratory or in the field. After time-frequency analysis, the energy release rate and event location coordinates are extracted. It should be noted that in the laboratory, a vacuum-coupled piezoelectric sensor array is attached to the rock sample surface, while in the field, a three-component acoustic emission probe set is drilled through a borehole. The raw waveform signal is first bandpass filtered to eliminate interference from the equipment's resonant frequency. A high-frequency, narrow-window capture mode is used for brittle fracture events in hard rock, while a low-frequency, wide-window capture mode is used for shear slip events in soft rock. The time-frequency analysis phase combines short-time Fourier transform and wavelet packet decomposition techniques. These techniques include extracting the energy integral of the signal's main frequency band as a baseline for the energy release rate and inverting the source coordinates using the time of arrival of the sensor array. When encountering in-situ tectonic noise, a rock mass characteristic frequency matching filter algorithm is used to retain only event signals that match the theoretical spectral characteristics of rock fracture. If there are differences in wave velocity between multiple layers of media within the monitoring area, ray tracing is used to correct for location offset errors. The final output is standardized event sequence data, including event timestamps, energy release rate scalar values, and three-dimensional spatial coordinates, which are then used by the dynamic correction module.

[0034] If either the physical morphology verification or the mathematical topology verification fails in step S5, the system automatically returns to step S4 to restart the parameter field optimization loop. It should be further explained that, in the specific implementation process, when the Hausdorff distance of the physical morphology verification exceeds the standard or the mathematical topology similarity is insufficient, the system automatically locates the spatial coordinates of the failure region and backtracks to the corresponding simulation time step. The characteristics of the acoustic emission event sequence at that moment are extracted, including: if there is a continuous deviation trend between the event energy release rate and the stress intensity factor, the size of the local inversion particle swarm is increased; if the event location shows multiple crack interference, the influence domain is expanded to cover the interference area.

[0035] When restarting the optimization cycle, a failure type adaptation strategy is applied. This includes prioritizing correction of the fracture energy field gradient distribution for geometric path deviations (physical verification failures) and focusing on adjusting the intensity field anisotropy at the crack bifurcation point for topological structure mismatches (mathematical verification failures). The previous parameter field version is retained as a baseline after each restart. When the number of consecutive restarts reaches the upper limit corresponding to the rock mass heterogeneity complexity level, a mesoscopic data re-acquisition process is forced. For verification failures in high-temperature and high-pressure environments, the thermomechanical coupling parameter field compensation module is additionally activated before re-entering the optimization cycle.

[0036] Dynamic mechanical parameter field construction, dynamic correction, and dual verification mechanisms synergize throughout the crack evolution process to quantify the dynamic impact of microscopic heterogeneity. It should be further explained that, during the specific implementation process, three types of mechanisms form a closed-loop control system throughout the crack evolution cycle: during the parameter field construction phase, mineral phase constraint boundaries are embedded to set the physically feasible domain for dynamic correction; during inversion optimization, the particle swarm sampling strategy is adaptively adjusted based on the characteristics of real-time acoustic emission events, and parameter field version tags are recorded for verification and backtracking; the dual verification results not only determine the final validity but also drive the restart of the targeted optimization cycle by locating the failure area.

[0037] For the brittle cracking conditions of hard rock, the collaborative system focuses on the response calibration of the stress intensity factor at the crack tip, and the parameter field update frequency is strictly synchronized with the acoustic emission main shock event; for the shear failure of soft rock, it switches to a long-term progressive optimization mode, and adds the consistency check of the shear band inclination in the topology verification stage.

[0038] When faced with complex hydraulic fracturing conditions, the collaborative system automatically activates the fluid-solid coupling pathway, including: constraining the particle swarm search space based on fracturing fluid viscosity parameters, and topologically verifying pore processing rules based on proppant distribution data. Throughout the process, all technical modules share a unified spatiotemporal coordinate system, ensuring that mineral phase boundaries, crack path coordinates, and acoustic emission source locations operate under the same datum.

[0039] By constructing a dynamic mechanical parameter field driven by microstructures, the inherent defect of traditional methods that simplify rocks into homogeneous materials is overcome, and microscopic characteristics such as mineral hardness gradients, grain boundary strength distributions, and pore morphology are effectively quantified into spatially distributed elastic modulus fields, tensile strength fields, and fracture energy fields, so that the numerical prediction of crack initiation and propagation can truly reflect the inherent heterogeneity of rocks for the first time. Combining acoustic emission monitoring data with the particle filter inversion algorithm, real-time calibration of parameter fields during crack evolution is achieved, improving the accuracy of crack network morphology prediction under complex stress paths. Compared with traditional numerical methods, the crack path prediction error has been reduced, and it still maintains stability and reliability in extreme working conditions such as high-in-situ stress rock bursts and shear slip of layered rock masses. A closed-loop quality control system is established through a dual verification mechanism. Physical morphology verification ensures geometric consistency between simulated cracks and real fracture surfaces, while mathematical topology verification quantifies network connectivity characteristics at the structural level. This dual-track screening minimizes false positives. This creates a full-cycle data closed loop, from microscopic characterization to macroscopic prediction. Parameter field construction embeds physical boundary constraints, corrects responses to real-time monitoring signals, and verification results drive targeted optimization cycles, reducing the cost of manual trial and error.

[0040] It should be noted that, in this document, relational terms such as first and second, etc., are used only to distinguish one entity or operation from another entity or operation, and do not necessarily require or imply any actual relationship or order between these entities or operations. Moreover, the terms "comprises," "comprising," or any other variations thereof are intended to cover non-exclusive inclusion, such that a process, method, article, or device comprising a series of elements includes not only those elements, but also other elements not explicitly listed, or elements inherent to such process, method, article, or device. In the absence of further limitations, an element defined by the phrase "comprising a..." does not exclude the presence of additional identical elements in the process, method, article, or device comprising the element.

[0041] While embodiments of the present invention have been shown and described, it will be appreciated by those skilled in the art that various changes, modifications, substitutions, and variations may be made to these embodiments without departing from the principles and spirit of the invention, and that the scope of the invention is defined by the appended claims and their equivalents.

Claims

1. A numerical simulation method for rock crack evolution based on data analysis, characterized in that: The steps include: S1: Obtain microstructural data of the target rock, including CT scan images, SEM images, and mineral energy spectrum data; S2: constructing a spatially distributed dynamic mechanical parameter field based on the mesoscopic structural data; S3: Use the extended finite element method to drive the numerical simulation of crack evolution, and integrate the microfracture acoustic emission monitoring data in real time during the simulation process; S4: When the deviation between the monitoring data and the simulation path exceeds a set threshold, dynamically correcting the dynamic mechanical parameter field; S5: The validity of the simulation results is determined through a dual verification mechanism of physical morphology and mathematical topology.

2. The method for numerical simulation of rock crack evolution based on data analysis according to claim 1, characterized in that: In step S2, a deep learning model is used to segment the mineral phases, pores and microcracks in the CT scan image and the SEM image to generate a mesostructure distribution map.

3. The method for numerical simulation of rock crack evolution based on data analysis according to claim 2, characterized in that: In step S2, the construction of the dynamic mechanical parameter field includes: designing a heterogeneous parameter transfer function based on the mineral energy spectrum data and the micromechanical experimental calibration results; mapping the mineral hardness and grain boundary strength into a spatially distributed elastic modulus field, tensile strength field and fracture energy field through the transfer function.

4. The method for numerical simulation of rock crack evolution based on data analysis according to claim 1, characterized in that: The dynamic correction of step S4 includes: coupling the particle filter algorithm with the extended finite element method to establish a dynamic correlation between the energy release rate of the acoustic emission event and the simulated crack tip stress intensity factor; using the acoustic emission energy release rate as the observation variable and the crack tip stress intensity factor as the state variable, and inverting the local fracture energy field through particle filter weight update.

5. The method for numerical simulation of rock crack evolution based on data analysis according to claim 4, characterized in that: The inversion is performed only on a local region of the crack growth front, and the update frequency is synchronized with the acoustic emission event rate.

6. The method for numerical simulation of rock crack evolution based on data analysis according to claim 1, characterized in that: The physical morphology verification of step S5 includes: aligning the simulated crack network topology with the 3D laser scanning point cloud of the actual failure section, calculating the Hausdorff distance between the two, and determining that the physical verification has passed when the Hausdorff distance is less than a preset ratio of the characteristic length.

7. The method for numerical simulation of rock crack evolution based on data analysis according to claim 6, characterized in that: The mathematical topology verification includes: extracting Betti numbers and Euler characteristic numbers of the simulated and experimental crack networks, calculating topological similarity through persistent homology analysis, and determining that the mathematical verification is passed when the similarity meets a predetermined topological matching standard.

8. The method for numerical simulation of rock crack evolution based on data analysis according to claim 1, characterized in that: The micro-fracture acoustic emission monitoring data is derived from acoustic emission signals collected in a laboratory or on-site, and the energy release rate and event location coordinates are extracted after time-frequency analysis.

9. The method for numerical simulation of rock crack evolution based on data analysis according to claim 1, characterized in that: If either the physical form verification or the mathematical topology verification fails in step S5, the process automatically returns to step S4 to restart the parameter field optimization cycle.

10. A method for numerical simulation of rock crack evolution based on data analysis according to any one of claims 3 to 7, characterized in that: The dynamic mechanical parameter field construction, dynamic correction and dual verification mechanism work synergistically throughout the crack evolution process to quantify the dynamic impact of mesoscopic heterogeneity.

Images