Rock dynamic fracture type-i stress intensity factor calculation method, system and device

Abstract

The application discloses a rock dynamic fracture I-type stress intensity factor calculation method, system and device, relates to the rock dynamics technical field, and the method comprises the following steps: collecting an image sequence of a granite test piece with a straight cut groove in a Brazilian disc through a dynamic impact test on the granite test piece; performing semantic segmentation on the image by using a Mask2former model, identifying a crack region, and automatically positioning a crack tip position at each moment based on a segmentation result; processing the image sequence based on a DIC method to obtain displacement field data of a test piece surface; determining a displacement field sub-region participating in calculation in a Williams equation based on the crack tip position and the displacement field data; and iteratively solving the I-type dynamic stress intensity factor based on a least square method by adjusting the shape and size of the displacement field sub-region, so that the measurement accuracy and efficiency of the dynamic stress intensity factor are significantly improved.

Description

Technical Field

[0001] This invention relates to the field of rock dynamics technology, specifically to a method, system, and apparatus for calculating the type I stress intensity factor of dynamic rock fracture. Background Technology

[0002] Rocks are natural objects formed from minerals and rock fragments under complex geological conditions, giving them characteristics such as heterogeneity, discontinuity, and anisotropy. Impact loads, as a significant factor in rock mass failure, pose a safety hazard in many underground and mining engineering projects. The impact failure of rock is often accompanied by the dynamic propagation of cracks. Fracture mechanics is a crucial theoretical foundation for studying material fracture. Some scholars, combining a series of solutions to the stress and displacement fields of cracks at the edge of infinite plates, have proposed the concept of stress intensity factors (SIFs). SIFs reflect the strength of the elastic stress field near the crack tip, thereby assessing crack initiation and propagation, and ultimately ensuring the safety of engineering projects.

[0003] Currently, dynamic solid-state fractures (SIFs) can be calculated during rock crack propagation by combining the DIC method with the Williams equation to determine the displacement field at the crack tip. Digital image correlation (DIC) calculates the displacement field on the object surface based on the changes in speckle characteristics of randomly distributed grids before and after deformation using a correlation criterion function. It offers advantages such as non-contact measurement and full-field strain measurement, and is widely used in deformation detection of various materials. In capturing the dynamic fracture of SHPB impact brittle rock using a high-speed camera, the location and shape of the rock cracks continuously change. Therefore, after obtaining the displacement field using the DIC method, finding the location of the crack tip and selecting the sub-regions involved in the equation calculation become key steps in calculating the SIF. Regarding crack tip identification, existing research has attempted to use deep learning models based on displacement gradient changes and U-Net for crack detection and localization. For example, the U-Net model is used for automatic identification of sandstone crack tips. Besides the U-Net model, various network structures have emerged in the field of deep learning segmentation, including SegNet, Mask R-CNN, DeeplabV3+, and Mask2former, all of which have achieved remarkable results on public datasets. On the other hand, regarding the selection of the calculation region for the Williams equation, a nonlinear fracture process zone (FPZ) is generated due to the singularity of the stress field around the crack tip in brittle materials. The FPZ cannot be analyzed using linear elastic fracture mechanics (LEFM). Most studies only determine the displacement point regions involved in the calculation through experience or limited-range testing, lacking systematic quantitative analysis.

[0004] In summary, although existing methods have achieved the measurement of dynamic stress intensity factor to some extent, due to the variable crack morphology and complex image features during the dynamic fracture process of rocks, traditional methods such as tip identification based on displacement gradient changes are easily interfered with in dynamic multi-crack propagation scenarios, resulting in large crack tip positioning errors and difficulty in determining the calculation area of ​​the Williams equation, thereby reducing the calculation efficiency and accuracy of stress intensity factor. Summary of the Invention

[0005] To address the shortcomings of existing technologies, such as large crack tip positioning errors and difficulty in determining the calculation region of the Williams equation, which reduces the efficiency and accuracy of stress intensity factor calculation, this invention proposes a method, system, and device for calculating Type I stress intensity factor in dynamic rock fracture, thereby solving the problems existing in the prior art.

[0006] A method for calculating the Type I stress intensity factor of dynamic rock fracture includes the following steps: Collect image sequences of rock specimens during the dynamic fracture process; The image sequence is input into a pre-trained crack semantic segmentation model based on the Mask2former architecture to perform semantic segmentation on the image sequence to identify the crack region; and the crack tip position in each frame image is located based on the crack region. Based on the digital image correlation (DIC) method, the subpixel displacement of the speckle pattern on the surface of the rock specimen in the image sequence before and after deformation is calculated to obtain displacement field data covering the surface of the rock specimen; and the nonlinear fracture process zone in front of the crack tip is identified based on the region with increased displacement gradient in the displacement field data. In the displacement field data, the nonlinear fracture process region and the crack region in front of the crack tip are screened out with the crack tip position as the center to determine the displacement field sub-regions that participate in the calculation of the Williams displacement field equation. The displacement field data of the displacement field sub-region and the polar coordinates determined by the crack tip position are input into the Williams displacement field equation for fitting and solving, and the type I stress intensity factor of rock dynamic fracture is calculated.

[0007] Furthermore, the step of inputting the image sequence into a pre-trained crack semantic segmentation model based on the Mask2former architecture to perform semantic segmentation on the image sequence to identify crack regions, and automatically locating the crack tip position in each frame of the image based on the crack regions, specifically includes the following steps: A semantic segmentation model is constructed by adding feature layers to the Pixel Decoder and Transformer Decoder of the Mask2former model; The image sequence of the rock specimen during the dynamic fracture process is input into the semantic segmentation model, and the segmentation result of each frame is output; the segmentation result includes a binary mask corresponding to the background, rock and crack respectively; Extract the binary mask representing the crack from the segmentation result, and find the coordinate extreme points along the crack propagation direction within the crack mask region on each side of the pre-existing crack. The extreme coordinate point found is determined as the position of the crack tip at that moment.

[0008] Furthermore, the digital image correlation (DIC) method calculates the subpixel-level displacement of the speckle pattern on the surface of the rock specimen in the image sequence before and after deformation to obtain displacement field data covering the surface of the rock specimen. This specifically includes the following steps: The image before the impact load was applied and before the specimen was deformed was selected from the image sequence as the reference image; On the reference image, the entire surface of the rock specimen or the region of interest containing the crack propagation path is used as the computational mesh; where each node of the mesh represents a displacement point to be calculated. For each current image in the image sequence, iterative calculations are performed to track changes in displacement points. Specifically, this includes: for each calculation node in the current image, a subset of pixels is selected centered on its corresponding position in the reference image; an iterative optimization algorithm is used to find the candidate region in the current image that is most similar to the pixel subset; with the candidate region as the initial value, an iterative algorithm is used to perform sub-pixel level optimization to obtain the displacement vector of each node. By calculating the displacement vector of each node on the computational grid, the displacement field data of the rock specimen surface at the corresponding time is obtained.

[0009] Furthermore, the identification of the nonlinear fracture process zone in front of the crack tip based on the region of increased displacement gradient in the displacement field data specifically includes the following steps: Select a region of interest (ROI) within the rock specimen sample area, and arrange several vertical displacement measurement lines on both sides of the pre-fabricated crack in the ROI. The vertical displacements on each displacement measurement line are extracted from the displacement contour lines to obtain the vertical displacement-vertical distance displacement curve on the displacement measurement line. Identify the boundaries of regions where displacement gradients abruptly occur based on the displacement curve of vertical displacement versus vertical distance; The boundaries of the abrupt change regions identified by each displacement measurement line are connected to determine the spatial range of the nonlinear fracture process zone in front of the crack tip, thereby obtaining the nonlinear fracture process zone in front of the crack tip.

[0010] Furthermore, the shape of the displacement field sub-region participating in the calculation of the Williams displacement field equation includes a circle centered at the crack tip with a radius of... R ( n A circular region, with the crack tip as the centroid of a square and a side length of ). D ( n A square region or a triangle with the crack tip as its centroid, with a side length of . T ( n Furthermore, the line connecting the apex of the equilateral triangle near the edge of the rock is parallel to the horizontal direction of the image.

[0011] Furthermore, it also includes optimizing the Williams displacement field equations by adjusting the shape and size of the calculated displacement field sub-regions, and iteratively solving for the Type I dynamic stress intensity factor based on the least squares method, specifically including the following steps: Multiple candidate computational sub-regions of different sizes are defined with the determined crack tip location as the center; For each displacement data point within the candidate calculation sub-region, its global coordinates and crack tip coordinates are converted to polar coordinates; and the displacement component of each displacement data point in the direction perpendicular to the crack surface is extracted to generate displacement observation values. Substituting the polar coordinates and displacement observations of each data point into the Williams displacement expansion based on linear elastic fracture mechanics, for the data point containing... M A computational sub-region containing data points is constructed. M One equation, regarding the unknown parameter [ , , The overdetermined system of equations; For each candidate computational subregion, perform the following operation: Set the number of truncation terms for the Williams series. N The initial value; the overdetermined system of equations is solved iteratively using the least squares method to obtain the initial value. N Unknown parameter series under value Gradually increase the number of items N ,when When the value of approaches a stable value as N increases, the Williams equation is considered to have converged. Among all candidate computational subregions, the region shape and size that enables the Williams equation to converge earliest are selected as the optimal computational region. The optimal computational domain is obtained under convergent conditions. Value, substitute into the formula The type I stress intensity factor of dynamic fracture in rock was calculated. .

[0012] Furthermore, dynamic impact tests were conducted on rock specimens with straight grooves using a split Hopkinson bar system, and an ultra-high-speed camera was used to acquire image sequences of the rock specimens during the dynamic fracture process.

[0013] This invention also includes a system for calculating the type I stress intensity factor of dynamic rock fracture, comprising: The acquisition module is used to acquire image sequences of rock specimens during the dynamic fracture process; The segmentation module is used to input image sequences into a pre-trained crack semantic segmentation model based on the Mask2former architecture, perform semantic segmentation on the image sequences to identify crack regions, and locate the crack tip position in each frame of the image based on the crack regions. The analysis module is used to calculate the subpixel displacement of the speckle pattern on the surface of the rock specimen in the image sequence before and after deformation based on the digital image correlation (DIC) method, so as to obtain displacement field data covering the surface of the rock specimen; and to identify the nonlinear fracture process zone in front of the crack tip based on the region of increased displacement gradient in the displacement field data. The computational region determination module is used to filter out the nonlinear fracture process region and the crack region in front of the crack tip in the displacement field data, with the crack tip location as the center, so as to determine the displacement field sub-region that participates in the calculation of the Williams displacement field equation. The solver module is used to input the displacement field data of the displacement field sub-region and the polar coordinates determined by the crack tip position into the Williams displacement field equation for fitting and solving, and to calculate the rock dynamic fracture type I stress intensity factor.

[0014] The present invention also includes a computer device for calculating the type I stress intensity factor of dynamic rock fracture, comprising: a memory, a processor, and a computer program stored in the memory, wherein the processor executes the computer program to implement the steps of the method for calculating the type I stress intensity factor of dynamic rock fracture.

[0015] The present invention also includes a readable storage medium storing a computer program, the computer program including program instructions, which, when executed by a processor, are used to perform the steps of the method for calculating the type I stress intensity factor of dynamic rock fracture.

[0016] This invention provides a method for calculating the Type I stress intensity factor of dynamic fracture in rocks, which has the following beneficial effects: This invention combines the advanced deep learning instance segmentation model Mask2former with digital image correlation (DIC) full-field measurement technology for dynamic rock fracture analysis. The Mask2former model solves the problem of automatic and accurate identification of dynamic crack tip positions, overcoming the shortcomings of traditional manual interpretation, which is subjective and inefficient. The DIC technology provides full-field, continuous displacement data, laying the data foundation for accurate calculations based on Williams equations. The two work together to achieve full-process automation and intelligence from image recognition to parameter calculation, significantly improving the measurement accuracy and efficiency of dynamic stress intensity factor. By fusing the DIC displacement field and Mask2former segmentation results, the interference of the fracture process zone FPZ and crack body on linear elastic theory is quantified and eliminated, purifying the calculation data source. This allows for the identification of the displacement field sub-regions participating in the calculation of the Williams displacement field equation. Using the convergence of the Williams equation series solution as an objective criterion, the optimal calculation sub-region is automatically determined, ensuring accurate convergence of the stress intensity factor while effectively reducing the computational load and avoiding the problem of solution divergence or excessive error caused by improper region selection. Attached Figure Description

[0017] Figure 1 This is a schematic diagram of the SHPB impact test system in an embodiment of the present invention; Figure 2 In the two-dimensional problem of the embodiments of the present invention q A pyramid-shaped diagram of a function; Figure 3 This is a coordinate diagram of the node plane isoparametric unit in an embodiment of the present invention; Figure 4 This is a diagram of the Mask2former network architecture with added feature layers in an embodiment of the present invention; Figure 5 This is a schematic diagram of the segmentation result between the original experimental image and Mask2former in Experiment BX4 of this embodiment of the invention; Figure 6 This is a schematic diagram of the Mask2former segmentation results, training labels, and strain field at the crack tip in an embodiment of the present invention. Figure 7 Figure 1 illustrates the distance between MCTP and LCTP and the two tips on both sides of the pre-cracks in experiments BX1 and BX2 during the crack initiation to arrest process in an embodiment of the present invention. Figure 8 Figure 2 illustrates the distance between MCTP and LCTP and the two tips on both sides of the pre-cracks in experiments BX1 and BX2 during the crack initiation to arrest process in this embodiment of the invention. Figure 9 This is a schematic diagram of the region of interest (ROI) and displacement measurement line in an embodiment of the present invention; Figure 10 This is a schematic diagram of the vertical displacement (V) - vertical distance (Y) curve on the displacement measurement line in an embodiment of the present invention; Figure 11 This is a schematic diagram illustrating the selection of non-computational regions in an embodiment of the present invention; Figure 12 This is a schematic diagram of the calculation region selection in an embodiment of the present invention; Figure 13 This is a schematic diagram illustrating the variation of the stress intensity factor in different range regions as a function of the higher-order term N (complete circular region) in an embodiment of the present invention. Figure 14 Figure 1 illustrates the test results of the variation of stress intensity factor in different regions and the variation of higher-order term N (circular region) in an embodiment of the present invention. Figure 15 Figure 2 illustrates the test results of the variation of stress intensity factor in different regions and the variation of higher-order term N (circular region) in an embodiment of the present invention. Figure 16 Figure 1 illustrates the test results of the variation of stress intensity factor in different regions and the variation of higher-order term N (square region) in an embodiment of the present invention. Figure 17 Figure 2 illustrates the test results of the variation of stress intensity factor in different regions and the variation of higher-order term N (square region) in an embodiment of the present invention. Figure 18 Figure 1 illustrates the test results of the variation of stress intensity factor in different regions and the variation of higher-order term N (equilateral triangle region) in an embodiment of the present invention. Figure 19 Figure 2 illustrates the test results of the stress intensity factor variation and the variation of the higher-order term N (equilateral triangle region) in different ranges in an embodiment of the present invention. Figure 20 This is a schematic diagram of the geometric model and mesh generation in an embodiment of the present invention; Figure 21 This is a schematic diagram of strain verification in numerical simulation in an embodiment of the present invention; Figure 22 This is a schematic diagram of the segmentation result of Mask2former and the damage process of the numerical simulation unit in an embodiment of the present invention; Figure 23 These are schematic diagrams of the Y-direction displacement based on DIC technology and the Y-direction displacement cloud based on numerical simulation, as shown in this embodiment of the invention. Figure 24 This refers to the calculation based on DIC-Mask2former in the SHPB experiment of this invention. A schematic diagram of the numerical results obtained by the J-integral method in numerical simulation; Figure 25 This is a schematic diagram of a method for calculating the type I stress intensity factor of dynamic rock fracture in an embodiment of the present invention. Detailed Implementation

[0018] The technical solutions of the present invention will be clearly and completely described below with reference to the accompanying drawings of the embodiments of the present invention. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments.

[0019] In previous studies, various methods have been used to find the crack tips of various materials. For example, some scholars have detected crack tips in metallic materials based on changes in displacement gradients; some scholars have automatically identified crack tips in sandstone in triaxial experiments based on the U-net model; and some scholars have used Williams equations in finite element software to find crack tips in a manner similar to Pac-Man. With the development of deep learning, methods for detecting and segmenting crack images using deep learning have been gradually applied to various engineering scenarios, such as cracks in bridges, tunnels, and road surfaces.

[0020] However, the deep learning models in these studies have certain limitations in segmenting and detecting rock fracture behavior, and their applicability is restricted by factors such as engineering background, crack texture, and temporal domain selection. Currently, there are few studies on real-time detection of dynamic cracks in rocks based on various experimental rock images as datasets. Some scholars have proposed a deep convolutional neural network model to extract microcracks using CT image slices of coal and rock as datasets. Other scholars have introduced the CrackU-Net model based on U-Net and proposed a novel spatiotemporal Bayesian inference (STBI) method for segmenting and quantifying cracks generated during the dynamic fracture process of coal and rock, achieving good results. In addition to the U-Net model, various network structures have emerged in the field of deep learning segmentation, including SegNet, Mask R-CNN, DeeplabV3+, and Mask2former. These models have all achieved remarkable results on public datasets. On the other hand, regarding the selection of the computational region for the Williams equation, a nonlinear fracture process zone (FPZ) is generated due to the singularity of the stress field around the crack tip in brittle materials. The FPZ cannot be analyzed using linear elastic fracture mechanics (LEFM). This region has been analyzed and discussed in many studies. Even if the influence of FPZ on the calculation of SIFs values ​​is excluded, the size and shape of the selected calculation region still deserve further discussion.

[0021] Based on this, this invention proposes a method for calculating the Type I dynamic stress intensity factor (SIF) of rock dynamic fracture. The aim is to quantitatively study the influence of SHPB impact loading on the Type I dynamic SIF at the crack tip in granite CSTBD specimens during dynamic fracture failure, using the DIC-Mask2former optimization model as a foundation. Factors affecting computational accuracy: This invention utilizes the SHPB experimental system and ultra-high-speed photography technology to conduct a Type I dynamic fracture test on granite CSTBD specimens, capturing images of the dynamic fracture process in real time and annotating the captured digital images. Based on the Mask2former model segmentation results, the crack tip location is identified. On this basis, the displacement field is obtained using DIC technology and discussed in two categories according to linear elastic fracture theory: non-computational region and computational region. The focus in the non-computational region is on quantizing FPZ to improve computational accuracy, while the computational region discusses the impact of region shape and size on computational accuracy. The influence of the value was investigated, and the above process was finally simulated and verified based on ANSYS / LS-DYNA.

[0022] like Figure 25 As shown, the method specifically includes the following steps: SHPB impact test of S1 and CSTBD granites: This invention uses granite from a certain region as the test material. To ensure the repeatability of the experiment, the samples were divided into 8 groups, denoted as Experiments BX(1-8). The basic mechanical parameters of the rock samples are shown in Table 1. Due to the extremely high hardness and brittleness of granite, its fracture characteristics are obvious, which is helpful for studying the rock fracture process and mechanism.

[0023] Table 1. Basic mechanical parameters of the rocks used in the experiment This invention utilizes an experimental setup with a Hopkinson bar (SHPB) as follows: Figure 1 As shown in (a), the system includes components such as a nitrogen cylinder, a high-pressure gas gun, a punch, a pulse shaper, a PMMA plate, an incident rod, a transmission rod, strain gauges, an oscilloscope, a damper, a high-speed camera system, a flash, and a computer. The video acquisition system uses a Kirana ultra-high-speed camera from Specialized Imaging, ultimately obtaining 8 sets of experimental images, with 180 images in each set, totaling 1440 images.

[0024] The use of straight-grooved Brazilian disc (CSTBD) specimens to study the fracture process of rocks under high strain loading is now recognized by the International Society for Rock Mechanics (ISRM). Figure 1 (b) is a model diagram showing the detailed dimensions of the rock sample. The diameter D of the disc is 50 mm, and the pre-fabricated notch length I, width w, and depth B are 10 mm, 1 mm, and 25 mm, respectively. The notch tip is polished with a diamond wire saw to ensure the generation of a typical Type I crack path under impact loading. Finally, a series of random speckles are generated using spray painting. Figure 1 (c)

[0025] S2, DIC-based Calculation: The digital image correlation method divides the image into several equally sized grids, and calculates the displacement field of the object surface by means of a correlation criterion function based on the change of speckle characteristics of the randomly distributed grids before and after deformation; the experiment of this invention uses VIC-2D software to measure the experimental images.

[0026] In fracture mechanics, the relationship between the displacement field at the crack tip and the stress intensity factor can be expressed by the Williams equation, which is solved using information about the tip displacement field. In the experiment, the rectangular coordinates of the crack and the surrounding area in the rock fracture image need to be converted into polar coordinates using formula (3):

[0027] In formula (3) and These are the coordinates of the displacement point at the crack tip; (Note: superscript) w (This is a annotation of the Williams equations for solving the SIF series formulas) and This represents the coordinates of the crack tip. r and are the polar coordinate components of the crack tip region, respectively.

[0028] The test of this invention is a type I fracture of a rock specimen. The displacement of the crack tip region is perpendicular to the crack surface, that is, the displacement field equations can be effectively solved by the displacement component in the Y direction; the simplified formula (4) is obtained, and the rigid body displacement parameter of the specimen is introduced. Data points at the tip of the crack k , ( k =1, 2, ... As shown in formula (5): In formula (4) , These represent the displacement components in the X and Y directions, respectively; G is the shear modulus of the material, which is used in plane stress problems. =(3- ) / (1+ In plane strain problems =3-4 ; Poisson's ratio; and The coefficients of the series terms, subscript It represents a series sequence.

[0029] Formula (5) - A set of 2 Group about unknown parameters , , , and Nonlinear equations ( It is the number of equations in the DIC software. T y The error has been calculated; iterative solutions were obtained using the least squares method. and ; and Same as Type I and Type II stress intensity factors and The relationship between the stress intensity factor and the stress intensity factor can be expressed as formula (6): S3, based on numerical simulation calculate: To further verify the method of the present invention ,based on J The principle of integral numerical calculation method is solved in finite element method. And verify the experimental results. Some scholars have proposed the equivalent region integration method to calculate... J Integrals, above J The integral can be expressed as an equivalent surface integral over a finite region near the crack tip, as shown in formula (7):

[0030] (7) For the plane crack problem, the above equation can be expanded into formula (8): (8) For a two-dimensional problem, the integration region A needs to be determined first. The inner boundary Γ0 of the integration region A is usually chosen as the crack tip. In this case, the integration region A is the region enclosed by the outer boundary Γ1 of the integration region. The outer boundary Γ1 of the integration region usually coincides with the element edge.

[0031] function q ( x , y This is only to facilitate numerical computation of the integral expression. However, at all element nodes in the integration region, the function... q ( x , y All of these must have definite values. It has already been proven that... J The calculated value of the integral is relative to the assumption. q The form of the function is not sensitive. q The function can have any form, but it must be within the boundary of the integration region. q The correct values ​​must be taken. For example, for plane stress or plane strain problems, there is Γ0. q=1, which typically corresponds to the crack tip, while on the outer boundary Γ1 of the integration region... q =0.

[0032] Figure 2 (a) in the figure gives the common application in two-dimensional problems q The pyramid form of the function Figure 2 (b) in the middle is Figure 2 (a) is projected onto the x,y plane where the function is... q ( x , y The value is 1 at the crack tip, then linearly varies to the outer boundary, where it becomes 0. The value at any location within a given element... q The value of the function can be obtained by interpolating the unitary function to obtain formula (9):

[0033] (9) in n This refers to the number of nodes in the unit. q I For the node q The value that the function can take. N I It is a unit shape function. Therefore q The partial derivatives of the function can be expressed as formula (10):

[0034] (10) in The coordinates are the isoparametric coordinates of the element.

[0035] In seeking q After taking the partial derivative of the function, only the displacement vector needs to be determined. u and v coordinates x The partial derivatives, as well as the stress components and strain energy density. w The value of can be solved using the Gaussian integral method. J Integration. The values ​​of these variables at Gaussian points can also be obtained through interpolation using element-shape functions.

[0036] The four-node isoparametric element of a two-dimensional plane problem is as follows: Figure 3 As shown. Figure 3 In (a), s and t are the isoparametric coordinates of the element. Figure 3 (b) x and y The coordinates of the element are global coordinates. Variables within the element can be obtained by interpolating the values ​​of variables at the element nodes using shape functions. Considering any point within the element, its coordinates can be expressed using nodal coordinates and shape functions, as shown in formulas (10) and (11):

[0037] (10) (11) The displacement component at this point is also approximated by nodal displacement and shape function in formula (12): (12) The point corresponding to q The value of the function can also be expressed as formula (13): (13) In the above formulas, N For shape function vectors, x and y The coordinate vector at the element node. u and v This represents the displacement vector at the element node. q At the unit node q The vectors formed by the values ​​can be expressed as formula (14): (14) The variables at the element nodes in the above formula can all be directly calculated using finite element software. The key step below is to determine the shape function based on the element type. Figure 3 As can be seen from the definition of the 4-node planar isoparametric element in (a), the isoparametric coordinate system st is located at the exact center of the element, and the isoparametric coordinate system is a natural coordinate system. The coordinate axes s and t do not need to be orthogonal, nor do they need to be parallel to the global coordinate axes x and y. By defining the isoparametric coordinates as +1 or -1, the edges and corners of the isoparametric element can be completely restricted. Assume... Figure 3 The malformed unit in (b) became as follows Figure 3 The regular cell shown in (a) allows the st coordinate system to be associated with the global coordinates x and y as shown in formulas (15) and (16):

[0038] (15) (16) in x e and y e Let x and y be the global coordinates at the center of the unit. We assume that the global coordinates x and y can be expressed in terms of isoparametric coordinates s and t as shown in formulas (17) and (18):

[0039] (17) (18) Substitute node coordinates x 1. x 2. x 3. x 4. y 1. y 2. y 3. y 4 and the corresponding values ​​of s and t can be used to... a i use s and t To express this, we finally obtain formula (19): (19) The shape functions in formula (19) are represented by formula (20): (20) The derivative of the shape function with respect to the isoparametric coordinates can be obtained as formula (21): (twenty one) Therefore, we have formula (22): (twenty two) The Jacobian matrix is ​​given by formula (23): (twenty three) The derivative of the shape function with respect to the global coordinates is given by formula (24): (twenty four) Therefore, the strain can be expressed as formula (25)-formula (27): (25) (26) (27) Similarly, it can be obtained q The differential formula for the function (28): (28) After obtaining the stress, the stress can be obtained through the linear elastic constitutive relation, which is expressed as formula (29) for plane stress problems: (29) After obtaining the stress and strain, the strain energy density w This can be expressed as formula (30): (30) Finally, the J integral within a certain unit can be expressed as formula (31): (31) In formula (31), the function I can be expressed as in formula (32): (32) This surface integral can be solved using the Gaussian integral method as shown in formula (33): (33) The coordinates of the four integration points in the formula are selected as shown in formula (34): , (34) Repeat the above process to calculate the J integral over all elements within the region, and sum the results to obtain the J integral value at the crack tip. From this, the formula for the Type I stress intensity factor (35) is obtained:

[0040] (35)

[0041] S4. Using Mask2former as the framework for crack semantic segmentation, traditional neural network models typically include an input layer, an output layer, and multiple hidden layers. Mask2former effectively combines feature pyramids and Transformers on top of traditional neural networks; high-resolution features optimize the segmentation effect for small targets; the order of self-attention and cross-attention mechanisms in the TransformerDecoder is changed, accelerating model convergence and improving performance; Mask loss only calculates sampled random points instead of the entire image, saving three times the training memory; based on the Mask cross-attention mechanism, the model only focuses on attention between foreground elements, without requiring background involvement, significantly shortening training time. Figure 4The image illustrates a Transformer-based mask classification and segmentation framework, similar in structure to Mask2former. Its main process is as follows: The backbone network extracts multi-scale features (e.g., ResNet outputs feature maps at different resolutions: 8×, 16×, 32×, 64×, 128×, 256×, 512× corresponding to the number of channels). The pixel decoder fuses and upsamples the features before feeding them into the Transformer decoder. The Transformer decoder includes: a self-attention module (modeling the relationships between queries); a masked attention module (combining object queries with image features, guiding attention through the predicted mask); each submodule is followed by "residual connections + layer normalization" (Add & Norm); and an FFN further processes the features. Finally, it outputs a set of prediction results consisting of "class labels + masks," achieving panoramic, instance, or semantic segmentation. The backbone is a pixel decoder that extracts low-resolution features from an image. The Pixel Decoder progressively upsamples the low-resolution features from the backbone output to generate high-resolution per-pixel embeddings. Finally, the Transformer Decoder is a transforming decoder that manipulates image features to process object queries. The final binary annotation prediction is decoded using the per-pixel embeddings of the object queries.

[0042] This invention adds new feature layers to the Pixel Decoder and Transformer Decoder based on the original Mask2former, in order to extract features at lower resolutions and make the model more sensitive to small targets. Figure 4 ResNet_d64 was used as the backbone layer. The ResNet_d64 layer consists of a 2D convolutional layer with 3 input channels and 64 output channels, a 7x7 kernel size, a 2x2 stride, and 3x3 padding to extract feature map information. After processing through five ResNet_d64 residual blocks, five feature levels were extracted. Four lower-resolution features were then output to a Transformer Decoder for decoding. The experimental images were labeled and trained using the network model to obtain the segmentation results.

[0043] The segmentation results from Mask2former contain three categories: background, rocks, and cracks. Figure 5Since the propagation direction of a Type I crack can be approximated as a ray, and the crack is considered as a single main crack propagating to both sides in this experiment, the crack tip (CTP) is represented in the image as the pixels with non-zero values ​​in the blue channel at the top and bottom.

[0044] Figure 6 This image shows the crack tip effect identified by Mask2former from the crack initiation to arrest process on the left side of the pre-existing crack in experiment BX1, along with a cloud map of the high-strain zone near the crack tip and a comparison of the reference strain map with the original experimental image. Figure 6 It can be intuitively seen that in experiment BX1, with the concentration and expansion of the strain field, the crack tip position (MCTP) identified by Mask2former and the marked crack tip position (LCTP) almost simultaneously start to move from the pre-existing crack in the negative X-axis direction and stop moving at T=18µs. This shows that the MCTP and the actual crack tip position in the experiment have a good match in both time and space. Figure 7 , Figure 8 The distance between MCTP and LCTP and the two tips is measured on both sides of the pre-cracked sides of experimental BX1 and BX2 from crack initiation to crack arrest.

[0045] Figure 7 , Figure 8 During the crack initiation stage (T=0-5µs), the displacements of MCTP and LCTP are small, and the difference in distance between them is less than 2mm. During the crack instability propagation stage (T=5-18µs), the displacements of MCTP and LCTP in the X direction increase sharply, and the difference in distance between their tips also increases. During the crack arrest stage (T=18-21µs), the displacements of MCTP and LCTP decrease, and the difference in distance between them is less than 3.5mm. This indicates that Mask2former is more effective at identifying crack tips during the crack initiation and arrest stages than during the crack instability propagation stage. Furthermore, the main direction affecting the distance between MCTP and LCTP is also the X direction. In other words, the X-axis coordinates are not accurate during the crack instability propagation stage of MCTP. One reason may be that the initiation of multiple cracks makes it difficult for the model to capture the features of the image. Another reason may be that crack generation is not a single two-dimensional problem. Fracture also occurs in the direction perpendicular to the circular surface of the specimen, causing the propagation process from the pre-existing crack to both positive and negative X-axis directions to become discontinuous in the image, which in turn leads to inaccurate identification of crack tips.

[0046] The error caused during the unstable crack propagation stage will not affect subsequent solutions. The measurement process, due to various influencing factors, only allows for the measurement of the crack initiation stage. value. Figure 7In (a), the three-dimensional plot shows the spatial coordinate changes of the labeled crack tip (LCTP) and the model-identified crack tip (MCTP) at the left end of the pre-existing crack during rock fracture in experiment BZ1, while the two-dimensional plot shows the distance between the two tips at the same time. Figure 7 (b) Figure 8 In the diagram, (c) and (d) represent the spatial information of the right end of the pre-crack in BZ1, the left end of the pre-crack in BZ2, and the LCTP and MCTP of the right end of the pre-crack in BZ2, respectively.

[0047] S5. Selection of the computational region The impact.

[0048] S5.1 Quantization and Pairing of Non-Computational Regions (FPZ and Cracks) The impact.

[0049] The dynamic propagation of cracks and the formation of the FPZ region in front of the crack tip provide a basis for accurate measurement. This presents a significant challenge, as numerous research reports indicate that the fracture toughness measured using SCB specimens is far lower than that measured using other ISRM recommended methods, attributing this large difference to the presence of the FPZ region near the crack tip. Therefore, FPZ region quantification directly affects the solution. The accuracy of the value.

[0050] Based on DIC technology, there are two commonly used methods to identify and characterize FPZ in brittle materials. One method is to directly identify FPZ based on the high-strain region at the crack tip, and the other method is based on the evolution of the displacement gradient region. In this invention, the displacement gradient method is used, and the crack initiation displacement contour map in BX2 is used as a calculation sample to illustrate the determination of FPZ in CSTBD granite samples. Calculation steps:

[0051] Step 1: Select a Region of Interest (ROI) within the sample area, and then arrange several vertical displacement measurement lines on both sides of the pre-existing crack in the ROI, such as... Figure 9 (x=14 to x=24.5 and x=38.5 to x=49), where the Y-axis represents the direction of the column vectors in the image matrix from top to bottom, and the X-axis represents the direction of the row vectors in the image matrix from left to right.

[0052] Step 2: Extract the vertical displacement (V) from the displacement contour lines on the displacement measurement line.

[0053] Step 3: Obtain as follows Figure 10 (a) The displacement curve of vertical displacement (V) - vertical distance (Y) on the displacement measurement line, by capturing the high displacement gradient area in the DIC displacement contour map, such as... Figure 10 (c) Obtain the FPZ region. Figure 10 (b) is the curve that does not produce a sudden change in displacement gradient.

[0054] Numerous studies have shown that the ends of the FPZ region are always after the crack tip, and the maximum width of the FPZ at different times approaches a constant greater than the crack width. Based on this theory, the non-computational region is selected as the area of ​​the rectangle enclosed by parallel lines at the x-coordinates of the FPZ ends and parallel lines at the y-coordinate of the maximum FPZ width on the experimental image of this invention. To ensure the rigor of the non-computational region, the area of ​​the crack identified by the Mask2former model is added as the final non-computational region. Figure 11 . Figure 12 The non-calculation region in the model completely includes FPZ and cracks, which can also avoid edge point calculation failure caused by DIC decorrelation to a certain extent.

[0055] S5.2, The shape and size of the displacement field calculation region affect The impact.

[0056] The choice of N terms in the Williams sequence in formula (5) will affect The calculation results. It varies with the number of higher-order terms N. Existing studies have shown that the DIC method can be used to test SCB specimens in triaxial experiments. The calculation area is filled when calculating the value (the calculation point covers the entire specimen). It converges and achieves good results. For example... Figure 13 As shown in the figure, it can be seen that during the full-area displacement field measurement, the higher-order term N increases. The value also fluctuates accordingly. When N > 35, N increases. It tends to a stable value. Figure 13 The left and right ends of the middle Both converged to 0.249 and 0.258 (MPa·m1 / 2), respectively. This invention uses this as a criterion to determine the correctness of the selected displacement field for measurement. Different shapes and sizes of displacement field regions involved in the calculation were tested to find the optimal solution that satisfies the convergence condition of the Williams equation while simultaneously reducing the computational domain. Figure 13 The ordinate, Stress intensity factor, represents the stress intensity factor.

[0057] This invention first determines the non-computational region, and then... Figure 12Three shapes and sizes of the computational region were tested, including (1) a circle (a circular region with radius R(n) centered at the crack tip); (2) a square (a square region with side length D(n) centered at the crack tip); and (3) an equilateral triangle (an equilateral triangle with side length T(n) centered at the crack tip, and the line connecting the tip to the vertex of the equilateral triangle near the rock edge is parallel to the Y-axis). The values ​​of R(n), D(n), and T(n) and the corresponding areas of the computational regions are shown in Table 2. Figure 14 , Figure 15 Figures 1 and 2 illustrate the test results of the variation of stress intensity factor in different regions and the variation of the higher-order term N (circular region), as shown in Figures 1 and 2. Figure 16 , Figure 17 Figures 1 and 2 illustrate the test results of the variation of stress intensity factor in different regions and the variation of the higher-order term N (circular region), as shown in Figures 1 and 2. Figure 18 , Figure 19 Figures 1 and 2 illustrate the test results of the variation of stress intensity factor in different regions and the variation of the higher-order term N (circular region). Figure 14 (a)-(d) Figure 15 (e)-(h) Figure 16 (a)-(d) Figure 17 (e)-(h) Figure 18 (a)-(d) Figure 19 In the equation, (e)-(g) represent circular, rectangular, and equilateral triangular regions of different areas, respectively, and are the dynamic stress intensity factors in the Williams equation. As N iterates, test results show that all three shapes can eventually satisfy the convergence condition by expanding the displacement field calculation region, and the test results also show... All converge at N≥30, and at N=30, the largest circular region is n=8 (left end: 0.228MPa·m1 / 2 / right end: 0.292MPa·m). 1 / 2 The largest square region n=8 (left end: 0.241MPa·m) 1 / 2 / Right end: 0.296MPa·m 1 / 2 The largest equilateral triangle region n=7 (left end: 0.188MPa·m) 1 / 2 / Right end: 0.190 MPa·m 1 / 2 ).

[0058] A circular region of 444.9 mm² is required for convergence. The value is close to the result obtained from the full-area displacement field calculation. The solution for the square region is correct. The required area is slightly larger than a circle, while the area obtained is an equilateral triangle. The result is smaller than that obtained from the displacement field calculation of the entire region, and the selected region needs to be large enough for convergence.

[0059] Table 2. Calculation Region Shape and Size Parameters Numerical simulation: This invention uses ANSYS / LS-DYNA software to simulate and analyze the above-mentioned SHPB impact experiment.

[0060] Model Establishment: To ensure calculation accuracy, Hopkinson's rods and CSTBD rock samples were created in the software using stretching operations in the X and Z directions, referencing the shape and size of actual experimental equipment and specimens. Based on this, a hexahedral mesh with 1,663,990 elements was created. Figure 20 The model consists of three parts: an incident rod, a rock specimen, and a transmission rod, all coaxially aligned. After establishing the geometric model, the material parameters and constitutive model for numerical simulation need to be determined. The Hopkinson rod material model is defined as an elastic model by adding the keyword *MAT_ELASTIC(001). In this invention, the Riedel-Hiermaier-Thoma (RHT) model in ANSYS / LS-DYNA and the ADD_EROSION keyword are used together to describe the failure behavior of the rock specimen under impact load. According to previous literature, this constitutive model has been successfully applied to study the failure and fracturing behavior of rock materials under dynamic load, and its effectiveness has been well verified. The main parameters of the rock specimen and the Hopkinson rod model are shown in Table 3.

[0061] To simulate the loading conditions of the specimen in the impact test, boundary conditions were set in LS-PrePost. First, the keyword *BOUNDARY_SPC_SET was added to restrict the rotation of the front end of the incident rod and the end of the transmission rod. Then, *CONTACT_ERODING_SURFACE_TO_SURFACE was set to define the contact surfaces between the incident rod and the rock and between the rock and the transmission rod. Finally, the incident wave was loaded using the keyword Load.

[0062] Table 3 Parameter Table Numerical simulation verification: Taking experiment BX4 as an example, this invention verifies the results by comparing experimental and simulation results from strain waveforms. Figure 21 Damage mode Figure 22 Y-direction displacement Figure 23 Comparison, Figure 21 This represents the changes in the values ​​of incident wave, transmitted wave, and reflected wave during rock failure, and verifies the accuracy of the numerical model by measuring the equilibrium state of the three waves. Figure 22In the image, (a) represents the result of the Mask2former model segmenting the original experimental image during the rock failure process. Figure 22 (b) in the figure represents the rock failure element (crack) propagation process in the numerical model; Figure 23 In the diagram, (a) represents the vertical displacement (V) contour plot obtained by DIC after processing the original experimental image during the rock failure process. Figure 23 (b) in the figure represents the vertical displacement (V) variation of the rock element in the numerical model. Verification of the numerical simulation: Under the same time axis, the strain curves measured by selecting one element from the incident rod and the transmission rod in the numerical simulation are compared with the strain curves measured by strain gauges attached to the incident rod and the transmission rod in the experiment. It can be concluded that the strain peak values ​​and waveform trends of the two are roughly the same.

[0063] Using the pre-existing crack initiation time (T=0µs) as a baseline, the results of segmenting experimental images based on the Mask2former model were compared with the failure elements in the numerical simulation. In the Mask2former model segmentation results, the crack initiates from both ends of the pre-existing crack, then the cracks at the left and right ends proliferate unstably along the negative and positive X-axis directions respectively, and finally the crack extends to the rock edge to arrest, accompanied by the initiation and coalescence of multiple cracks. Similarly, in the numerical simulation, the crack also initiates near the pre-existing crack, accompanied by the propagation of multiple cracks. The results show that the failure models of the numerical simulation and the impact experiment are almost identical, and both exhibit the characteristics of Type I fracture in the crack initiation stage.

[0064] Using the pre-crack initiation time as a baseline (T=0µs), and selecting the same displacement unit (min=-0.02mm, max=0.02mm), the Y-direction displacement contour maps obtained based on the DIC method and experimental images were compared with those from numerical simulation. The results show that before the pre-crack initiation, the Y-direction displacement is concentrated on the upper and lower sides of the pre-crack, with opposite values, and is axially symmetrical about the X-axis of the specimen center. During the crack propagation stage (T=0-21µs), the Y-direction displacement concentration in both numerical simulation and impact experiment diffuses from the upper and lower sides of the pre-crack center towards the edges of the specimen. Finally, the Y-direction displacement values ​​on the upper and lower sides of the pre-crack are opposite and overflow the upper and lower halves of the specimen, forming a semi-circular shape. This sufficiently demonstrates the consistency between the changes in Y-axis displacement in numerical simulation and SHPB impact experiment.

[0065] In the above process, the crack initiation stage can be regarded as Type I fracture. Figure 24 (a) shows the impact experiments of BX1-BX4 based on the Williams equation and the DIC-Mask2former optimization method. The values ​​are obtained by solving the displacement and coordinate data of elements and nodes exported in the finite element software LS-prepost based on the equivalent region integration method. Values. Most of the dynamics of experiment BX(1,3,4) The value is between 0.25 and 1 (MPa·m 1 / 2 The experimental values ​​for BX2 showed a slight increasing trend between 3 and 6 µs. The downward trend in the value may be caused by the simultaneous initiation of multiple cracks, such as... Figure 24 As shown in the screenshot (a), the propagation of a set of parallel microcracks eventually leads to mutual antagonism at the crack tips, thus affecting the output of the DIC displacement field and consequently the... The values ​​are unstable. During the same period, the numerical simulation... The value showed an increasing trend from 0.72 to 2.73 (MPa·m). 1 / 2 This may be because the material in the numerical simulation is more homogeneous than real rock, causing the rock failure and crack trajectory to more closely resemble the ideal state. Figure 24 (b) shows each stage of the four experimental groups. The mean of the values ​​was compared with those in numerical simulation. The difference was analyzed using a histogram, showing a value of 0.16 at T=0µs. The difference gradually increased with time, reaching 1.00, 1.33, and 1.48 at T=3, 6, and 9µs, respectively. In summary, the numerical simulation... The value is slightly larger than that in the impact test. The values, and this numerical difference is likely caused by the damage and non-homogeneity within the rock material, in the simulation. Values ​​obtained from experimental measurements The magnitude of the value difference is within a certain range and the overall trend is roughly the same.

[0066] This invention optimizes the Williams equation solution for dynamic fracture of granite in SHPB impact tests based on deep learning and DIC technology. The process was studied, and the experimental results were verified and analyzed based on ANSYS / LS-DYNA finite element software. The specific conclusions are as follows: (1) The present invention uses the Mask2former model with added feature layers as the object for training and segmentation of experimental images and captures the real-time crack tip position in the results. (2) The study found that the size and shape of the computational region affect the convergence of the equation. The minimum computational area required for the convergence of the circular region equation is 444.9 mm. 2 FPZ was quantized based on the displacement method, and FPZ enhancements were removed from the computational domain. The solution accuracy. (3) The above experiment was numerically simulated in the finite element software ANSYS / LS-DYNA, and the corresponding solution was obtained at the crack tip according to the equivalent region integration method. The results show that the numerical simulation and experimental results are consistent. The differences are within a certain range and the overall trend is roughly the same.

[0067] Based on the same inventive concept, this invention also proposes a system for calculating the Type I stress intensity factor of dynamic rock fracture, comprising: The acquisition module is used to acquire image sequences of rock specimens during the dynamic fracture process.

[0068] The segmentation module is used to input image sequences into a pre-trained crack semantic segmentation model based on the Mask2former architecture, perform semantic segmentation on the image sequences to identify crack regions, and automatically locate the crack tip position in each frame of the image based on the crack regions.

[0069] The analysis module is used to analyze image sequences based on the Digital Image Correlation (DIC) method to obtain displacement field data on the surface of rock specimens; and to identify the nonlinear fracture process zone in front of the crack tip based on the region of increased displacement gradient in the displacement field data.

[0070] The computational region determination module is used to filter out the nonlinear fracture process region and the crack region in front of the crack tip in the displacement field data, with the crack tip location as the center, so as to determine the displacement field sub-regions involved in the calculation of the Williams displacement field equation.

[0071] The solver module is used to input the displacement field data of the displacement field sub-region and the polar coordinates determined by the crack tip position into the Williams displacement field equation for fitting and solving, and to calculate the rock dynamic fracture type I stress intensity factor.

[0072] The present invention also proposes a computer device for calculating the type I stress intensity factor of dynamic rock fracture, comprising: a memory, a processor, and a computer program stored in the memory, wherein the processor executes the computer program to implement the steps of the method for calculating the type I stress intensity factor of dynamic rock fracture.

[0073] The present invention also proposes a readable storage medium storing a computer program, the computer program including program instructions, which, when executed by a processor, are used to perform the steps of a method for calculating the type I stress intensity factor of dynamic rock fracture.

[0074] The above description is only a preferred embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any equivalent substitutions or modifications made by those skilled in the art within the scope of the technology disclosed in the present invention, based on the technical solution and inventive concept of the present invention, should be covered within the scope of protection of the present invention.

Claims

1. A method for calculating the Type I stress intensity factor of dynamic rock fracture, characterized in that, Includes the following steps: Collect image sequences of rock specimens during the dynamic fracture process; The image sequence is input into a pre-trained crack semantic segmentation model based on the Mask2former architecture to perform semantic segmentation on the image sequence in order to identify the crack region. And based on the crack region, the location of the crack tip in each frame of the image is determined; Based on the digital image correlation (DIC) method, the subpixel-level displacement of the speckle pattern on the surface of the rock specimen in the image sequence before and after deformation is calculated to obtain displacement field data covering the surface of the rock specimen. Furthermore, the nonlinear fracture process zone in front of the crack tip was identified based on the region of increased displacement gradient in the displacement field data; In the displacement field data, the nonlinear fracture process region and the crack region in front of the crack tip are screened out with the crack tip position as the center to determine the displacement field sub-regions that participate in the calculation of the Williams displacement field equation. The displacement field data of the displacement field sub-region and the polar coordinates determined by the crack tip position are input into the Williams displacement field equation for fitting and solving, and the type I stress intensity factor of rock dynamic fracture is calculated.

2. The method for calculating the type I stress intensity factor of dynamic rock fracture according to claim 1, characterized in that, The process of inputting the image sequence into a pre-trained crack semantic segmentation model based on the Mask2former architecture to perform semantic segmentation on the image sequence to identify crack regions and automatically locate the crack tip position in each frame of the image based on the crack regions specifically includes the following steps: A semantic segmentation model is constructed by adding feature layers to the Pixel Decoder and Transformer Decoder of the Mask2former model; The image sequence of the rock specimen during the dynamic fracture process is input into the semantic segmentation model, and the segmentation result of each frame is output; the segmentation result includes a binary mask corresponding to the background, rock and crack respectively; Extract the binary mask representing the crack from the segmentation result, and find the coordinate extreme points along the crack propagation direction within the crack mask region on each side of the pre-existing crack. The extreme coordinate point found is determined as the position of the crack tip at that moment.

3. The method for calculating the type I stress intensity factor of dynamic rock fracture according to claim 1, characterized in that, The digital image correlation (DIC) method calculates the subpixel-level displacement of the speckle pattern on the surface of the rock specimen in the image sequence before and after deformation to obtain displacement field data covering the surface of the rock specimen. Specifically, it includes the following steps: The image before the impact load was applied and before the specimen was deformed was selected from the image sequence as the reference image; On the reference image, the entire surface of the rock specimen or the region of interest containing the crack propagation path is used as the computational mesh; where each node of the mesh represents a displacement point to be calculated. For each current image in the image sequence, iterative calculations are performed to track changes in displacement points. Specifically, this includes: for each calculation node in the current image, a subset of pixels is selected centered on its corresponding position in the reference image; an iterative optimization algorithm is used to find the candidate region in the current image that is most similar to the pixel subset; with the candidate region as the initial value, an iterative algorithm is used to perform sub-pixel level optimization to obtain the displacement vector of each node. By calculating the displacement vector of each node on the computational grid, the displacement field data of the rock specimen surface at the corresponding time is obtained.

4. The method for calculating the type I stress intensity factor of dynamic rock fracture according to claim 1, characterized in that, The identification of the nonlinear fracture process zone in front of the crack tip based on the region of increased displacement gradient in the displacement field data specifically includes the following steps: Select a region of interest (ROI) within the rock specimen sample area, and arrange several vertical displacement measurement lines on both sides of the pre-fabricated crack in the ROI. The vertical displacements on each displacement measurement line are extracted from the displacement contour lines to obtain the vertical displacement-vertical distance displacement curve on the displacement measurement line. Identify the boundaries of regions where displacement gradients abruptly occur based on the displacement curve of vertical displacement versus vertical distance; The boundaries of the abrupt change regions identified by each displacement measurement line are connected to determine the spatial range of the nonlinear fracture process zone in front of the crack tip, thereby obtaining the nonlinear fracture process zone in front of the crack tip.

5. The method for calculating the type I stress intensity factor of dynamic rock fracture according to claim 1, characterized in that, The shape of the displacement field sub-region involved in the calculation of the Williams displacement field equation includes a circle centered at the crack tip with a radius of [missing information]. R ( n A circular region, with the crack tip as the centroid of a square and a side length of ). D ( n A square region or a triangle with the crack tip as its centroid, with a side length of . T ( n Furthermore, the line connecting the apex of the equilateral triangle near the edge of the rock is parallel to the horizontal direction of the image.

6. The method for calculating the type I stress intensity factor of dynamic rock fracture according to claim 1, characterized in that, It also includes optimizing the Williams displacement field equations by adjusting the shape and size of the calculated displacement field sub-regions, and iteratively solving for the Type I dynamic stress intensity factor based on the least squares method, specifically including the following steps: Multiple candidate computational sub-regions of different sizes are defined with the determined crack tip location as the center; For each displacement data point within the candidate calculation sub-region, its global coordinates and crack tip coordinates are converted to polar coordinates; and the displacement component of each displacement data point in the direction perpendicular to the crack surface is extracted to generate displacement observation values. Substituting the polar coordinates and displacement observations of each data point into the Williams displacement expansion based on linear elastic fracture mechanics, for the data point containing... M A computational sub-region containing data points is constructed. M One equation, regarding the unknown parameter [ , , The overdetermined system of equations; For each candidate computational subregion, perform the following operation: Set the number of truncation terms for the Williams series. N The initial value; the overdetermined system of equations is solved iteratively using the least squares method to obtain the initial value. N Unknown parameter series under value Gradually increase the number of items N ,when When the value of approaches a stable value as N increases, the Williams equation is considered to have converged. Among all candidate computational subregions, the region shape and size that enables the Williams equation to converge earliest are selected as the optimal computational region. The optimal computational domain is obtained under convergent conditions. Value, substitute into the formula The type I stress intensity factor of dynamic fracture in rock was calculated. .

7. The method for calculating the type I stress intensity factor of dynamic rock fracture according to claim 1, characterized in that, Dynamic impact tests were conducted on rock specimens with straight grooves using a split Hopkinson bar system, and an ultra-high-speed camera was used to acquire image sequences of the rock specimens during the dynamic fracture process.

8. A system for calculating the type I stress intensity factor of dynamic rock fracture, characterized in that, include: The acquisition module is used to acquire image sequences of rock specimens during the dynamic fracture process; The segmentation module is used to input image sequences into a pre-trained crack semantic segmentation model based on the Mask2former architecture, and to perform semantic segmentation on the image sequences to identify crack regions. And based on the crack region, the location of the crack tip in each frame of the image is determined; The analysis module is used to calculate the subpixel-level displacement of the speckle pattern on the surface of the rock specimen in the image sequence before and after deformation based on the digital image correlation (DIC) method, so as to obtain displacement field data covering the surface of the rock specimen. Furthermore, the nonlinear fracture process zone in front of the crack tip was identified based on the region of increased displacement gradient in the displacement field data; The computational region determination module is used to filter out the nonlinear fracture process region and the crack region in front of the crack tip in the displacement field data, with the crack tip location as the center, so as to determine the displacement field sub-region that participates in the calculation of the Williams displacement field equation. The solver module is used to input the displacement field data of the displacement field sub-region and the polar coordinates determined by the crack tip position into the Williams displacement field equation for fitting and solving, and to calculate the rock dynamic fracture type I stress intensity factor.

9. A computer device for calculating the Type I stress intensity factor of dynamic rock fracture, characterized in that, include: The memory, the processor, and the computer program stored in the memory, wherein the processor executes the computer program to implement the steps of the method for calculating the type I stress intensity factor of dynamic rock fracture as described in any one of claims 1-7.

10. A readable storage medium, characterized in that, The readable storage medium stores a computer program, which includes program instructions that, when executed by a processor, perform the steps of the method for calculating the type I stress intensity factor of dynamic rock fracture as described in any one of claims 1-7.

Images