A simulation method for the deformation and fragmentation of quasi-brittle materials under the action of supercritical CO2

Abstract

The present invention discloses a simulation method for the deformation - fragmentation of quasi - brittle materials under the action of supercritical CO2, which includes the following steps: statistically analyzing the distribution of natural fractures in coal - rock masses, determining geometric parameters, and writing a fracture generation program; establishing a numerical calculation model, dividing entity units, embedding fracture units at the junctions of natural fractures and entity units and automatically generating a set; successively constructing mechanical constitutive relations reflecting the plastic deformation of coal - rock matrix, the mixed - mode fracture of non - penetrating fractures, and the separation / extrusion / compression - shear friction between blocks after the complete fracture of coal - rock; determining the evolution law of coal - rock mechanical parameters with the action time of supercritical CO2, inputting material parameters and setting initial conditions; setting the increment step size, output parameters, and state variables at the moment when fracture units fail and contact pairs are activated. The present invention provides an effective simulation tool for the evolution of the sealing performance of caprocks under the corrosion of CO2 and the action of complex stress coupling under the conditions of CO2 geological storage.

Description

Technical Field

[0001] The present invention relates to the technical fields of CO2 geological storage and carbon emission reduction, and particularly relates to a simulation method for the deformation and fragmentation of quasi-brittle materials under the action of supercritical CO2. Background Art

[0002] CCUS (Carbon Capture, Utilization and Storage) technology is a key technology for reducing CO2 emissions during fossil energy power generation and industrial production processes, and is also a backup technology for China to achieve carbon neutrality. CO2 geological storage is the core component of CCUS technology, which determines the development potential and direction of CCUS technology. The key to the long-term and safe storage of CO2 is no leakage in the limited geological reservoir environment, and the core element is the change characteristics of the deformation and fragmentation mechanical properties of the caprock under the long-term action of CO2 and the coupling action of deep stress. In deep coal seams, CO2 easily reaches the supercritical state conditions (7.38 MPa, 31.4 °C), and this corrosive fluid will have a significant impact on the properties and mechanical properties of the caprock. Using a reasonable numerical simulation method to predict the development height of caprock deformation and fragmentation under the coupling action of supercritical CO2 and stress is of great significance for quantitatively analyzing the CO2 leakage law, optimizing CO2 injection parameters, storage location, etc.

[0003] In recent years, domestic and foreign scholars have mainly used finite element and discrete element methods to carry out simulation studies on the deformation and fracture of quasi-brittle media such as CO2 geological storage and coal and rock. Based on the finite element method, using elastic-plastic mechanics theory and damage mechanics theory, study the influence range of stress, displacement, plastic zone, damage zone, etc. of the caprock caused by CO2 injection into the formation. Based on the discrete element method, coal and rock are simplified into an aggregate composed of basic components such as masonry and spheres. The intact coal and rock are rigid bodies or elastic bodies, and the Coulomb slip model is mostly used for the fracture part, and then the fragmentation process of coal and rock under complex stress states is numerically calculated. However, the fractures in the caprock are randomly distributed. Under the combined influence of in-situ stress and CO2 injection pressure, the above-mentioned intact coal rock mass and fractures are often subjected to complex stress, resulting in complex response processes such as plastic deformation, mixed fracture mode, separation / extrusion / shear friction, etc.; more importantly, supercritical CO2 will corrode components such as quartz, calcite, and clay minerals in the rock, resulting in a decrease in its strength, fracture energy, and shear mechanical properties. Ignoring these two factors will seriously affect the accuracy of the numerical simulation results of CO2 geological storage. Summary of the Invention

[0004] To solve the above problems, the present invention provides a simulation method for the deformation - fragmentation of quasi - brittle materials under the action of supercritical CO2, which can analyze the deformation - fragmentation mechanism of the caprock under the long - term coupling action of deep stress and supercritical CO2 corrosion, and the leakage law of CO2 along the caprock and its internal fractures, providing an effective simulation means for optimizing the parameters of CO2 geological storage.

[0005] To achieve the above object, the technical solution adopted by the present invention is as follows:

[0006] A simulation method for the deformation - fragmentation of quasi - brittle materials under the action of supercritical CO2, comprising the following steps:

[0007] S1. Statistically analyze the distribution of natural fractures in coal - rock mass, determine geometric parameters, and write a fracture generation program;

[0008] S2. Establish a numerical calculation model, divide entity units, embed fracture units at the junctions of natural fractures and entity units, and automatically generate a set;

[0009] S3. Sequentially construct mechanical constitutive relations reflecting the plastic deformation of coal - rock matrix, the mixed - mode fracture of non - penetrating fractures, and the separation / extrusion / compression - shear friction between blocks after the complete fracture of coal - rock; compile a finite - element - discrete - element FDEM calculation program to realize the whole process of plastic deformation - mixed - mode fracture - fragmentation friction of coal - rock;

[0010] S4. Experimentally determine the evolution law of coal - rock mechanical parameters with the action time of supercritical CO2, input material parameters, and set initial conditions;

[0011] S5. Set the increment step size, output parameters, and state variables at the moment of fracture unit failure and contact pair activation, and calculate the model responses before and after unit failure through finite - element and discrete - element calculation units respectively;

[0012] S6. Analyze the deformation - fragmentation results of the loaded coal - rock specimens under the action of supercritical CO2 for different times.

[0013] Furthermore, in step S1, the CT scanning method is used to determine the distribution of main fractures in coal - rock, statistically analyze the frequency of fractures appearing at different azimuth angles, and statistically analyze the average spacing and its variance of adjacent fractures; on this basis, a main - fracture automatic generation program is compiled. Specifically: according to the average spacing and its variance of the main fractures determined by scanning, randomly generate discrete points in MATLAB, determine the coordinates of the discrete points, connect the scattered points into Delaunay triangles, enclose Voronoi polyhedra by the perpendicular bisectors of the adjacent sides of the triangles, generate multiple groups of Voronoi polyhedra, and then select a set of Voronoi polyhedra whose azimuth angle results are close to those obtained by CT scanning; simulate the spatial distribution of main fractures in coal - rock by each face of the above - mentioned Voronoi polyhedra.

[0014] Further, the step S2 includes the following steps:

[0015] S2.1. Establish a numerical calculation model and divide the entity unit mesh;

[0016] Based on step S1, further determine the shape and size of the research object according to the purpose of numerical simulation, and use tools such as "cutting" and "rounding" to cut the Voronoi polyhedron model into the required numerical model;

[0017] Before mesh division, determine the maximum element size L according to the following formula:

[0018] L = 9πEF e / [32(1 - μ 2 )σ p

[0019] In the formula, E is the elastic modulus; μ is the Poisson's ratio; F e is the fracture energy; σ p is the tensile strength;

[0020] Seed the numerical model through the "Seed edge" of the Mesh module, with the seed spacing being L, and then divide the entity unit mesh of the numerical model through the "Seed part" command;

[0021] S2.2. Embed crack elements at the junction of natural fractures and entity units and automatically divide the set.

[0022] (1) Update the node numbers of the entity units;

[0023] After mesh division, generate an ".inp file" through CAE, copy all the numbers below the keywords "*nset" and "*elset" and above "*asssembly" in the file to Excel, and find the maximum node number n max and the maximum element number e max ;

[0024] Keep all node coordinates unchanged, copy the nth i (0 < i ≤ n max ) node a times, and at the same time keep the node number of the element with the smallest element number as n i,min unchanged, and increase the node numbers of the remaining elements to [(a - 1)×10 + n i,min , where a is the number of elements shared by a certain node;

[0025] (2) Create the node numbers and element numbers of the crack elements;

[0026] Taking the quadrilateral entity unit as an example, denote the node numbers after the above re - numbering as N 01 ,N​02 , N 03 , N 04 ; The adjacent elements sharing the n 0i -th node are numbered as e 01 , e 02 , e 03 , e 04 . Combine the new node numbers and the corresponding elements into an ordered array:

[0027] Array = (e 01 , e 02 , e 03 , e 04 , N 01 , N 02 , N 03 , N 04 )

[0028] Reverse the order of N 01 to N 04 , and the node numbers of the crack elements between the quadrilateral solid elements e 01 , e 02 , e 03 , e 04 can be obtained;

[0029] The smallest number of the crack element is determined by the integer power of 10 obtained by adding 1 to the digit where the largest number of the solid element is located; for example, if the largest node number of the solid element is 756, the smallest number of the crack element is 1000.

[0030] (3) Embed crack elements at the main fissures;

[0031] Add the keyword "*Elset = cohelemAll" after the keyword "*Part" and before the keyword "*Assembly" in the ".inp file"; copy the corresponding e i and N i numbers in the order that the first column is the crack element number and the 2nd to 5th columns are the node numbers, so as to embed zero-thickness crack elements at all main fissures and all boundaries of the solid elements;

[0032] (4) Automatically divide the crack element set;

[0033] Copy the crack element numbers e i and their node coordinates N i (x i , y i , z i ) to Excel, and use the inner product method to solve the direction vector d coh = (d1, d2, d3) of all crack elements, and the direction vectors d of each face of the Voronoi polyhedronip , if the two are collinear, then calculate the distance from any point on the crack element to the Voronoi surface by using the spatial vector method. If it is 0, it is determined that the crack element is located on the surface of the Voronoi polyhedron; copy the numbers of all crack elements and node numbers that meet the above three conditions into the.inp file, and add the keyword "*Elset=cohElem01" after the keyword "*Part" and before the keyword "*Assembly", so as to form a set named "cohelem01" for all crack elements on the main fissures;

[0034] On this basis, obtain the set of crack elements "cohelem02" at the boundaries of all solid elements through Boolean operations.

[0035] Further, the step S3 includes the following steps:

[0036] S3.1. Elastic-plastic deformation of intact coal and rock mass;

[0037] In the elastic-plastic deformation stage, each element obeys the generalized Hooke's law:

[0038] σ = E0:(ε - ε p )

[0039] where σ is the stress, E0 is the elastic stiffness matrix, and ε is the total strain. The plastic strain is determined by the yield function and the plastic potential function. The expression of the yield function F is:

[0040]

[0041] where F is the yield function, <g>is the Macaulay symbol, p = trace(σ) / 3, γ = 3(1 - K c ) / (2K c - 1), q = [2(S:S) / 3] 1 / 2 , S = σ + pI where σ i,max (i = 1, 2, 3) are the maximum effective principal stresses, is the ratio of biaxial to uniaxial compressive yield stress, σ t and σ c (ε p ) are the tensile and compressive stresses respectively; K c is the parameter controlling the shape of the yield surface in the deviatoric plane, and the value for rock-like materials is 0.667;

[0042] G = [(δσ t tanψ) 2 + q 2 - ptanψ

[0043] where δ is the cusp curvature of the meridian of the plastic potential function in the tensile section in the q - p plane, generally taking a value of 0.1; ψ is the dilation angle;

[0044] S3.2. Mixed-mode fracture of non-penetrating cracks;

[0045] The rock fracture process is divided into two stages: elastic deformation and mixed-mode fracture;

[0046] The constitutive relationship in the elastic stage is:

[0047] σ c = D 0,c ε c

[0048] Once the following conditions are reached, the ductile fracture process begins:

[0049]

[0050] where D 0,c is the elastic stiffness matrix, σ c,n , σ c,s , σ c,t are the normal and two tangential stresses (peak stresses), and ε c is the strain;

[0051] To derive the constitutive equation for mixed-mode ductile fracture, first establish the fracture energies G c,n , G c,S (G c,S = G c,s + G c,t , G c,s and G c,t are the fracture energies in two tangential directions), G c,m The expression is:

[0052]

[0053] where j represents the n, S, m (i.e., tensile, shear, mixed-mode) fracture modes;

[0054] The expression for the fracture energy mixing ratio is:

[0055] ξ = G c,S / (G c,n + G c,S )

[0056] The criterion for complete fracture of coal and rock is:

[0057]

[0058] where χ are the tensile and shear fracture energies and material parameters at complete fracture respectively;

[0059] Under tensile, shear, and mixed-mode fracture modes, the constitutive equations of each element are:

[0060]

[0061] where d c,j , D c,j , σ c,j , S c,j are the damage variable, elastic modulus, stress, plastic displacement, and total displacement respectively, and there are:

[0062]

[0063] where S c,n and S c,S are the displacements under pure tensile and shear conditions respectively;

[0064] Assume that the plastic displacement under the mixed-mode fracture mode is equal to the shear plastic displacement, i.e., there is

[0065]

[0066] where the unknown parameter d c,j is determined through the following steps: ① Use three-point bending tests and loading-unloading direct shear tests to obtain the σ c,m -S c,m curve, σ c,n -S c,n curve, σ c,S -S c,S Curve, calculate the corresponding damage evolution curve and fracture energy ② The "coal-rock complete fracture criterion expression" was used to analyze the multiple groups of experimental The data are fitted to obtain the fracture energy mixing ratio ξ under the same material parameter χ; ③ On this basis, the "peak stress equation" is combined to obtain the peak load under a certain ξ The corresponding σ c,n With σ c,S The peak displacement component S is calculated according to the elastic expression c,n ,S c,S ; ④ Combine the "fracture energy expression" and "constitutive equation" to obtain any σ c,m The corresponding σ c,n With σ c,S Component and the corresponding displacement component S c,n ,S c,S ⑤Replace the stress and displacement in the mixed mode with the mechanical parameters of pure tension and pure shear, and calculate σ by the "constitutive equation" c,n , σ c,S and S c,n ,S c,S Plastic strain and tensile / shear damage variables under the conditions; ⑦ Calculated by the "constitutive equation" to obtain d under any mixing ratio c,m -S c,m Relationship; thereby determining the mixed ductile fracture constitutive equation;

[0067] S3.3, separation / compression / shear friction of through-fractures or blocks;

[0068] Once the fracture energy satisfies the fracture energy expression, the rock is completely fractured, and the zero-thickness crack unit on the fracture surface is deleted in the numerical calculation to eliminate its mechanical influence;

[0069] On this basis, the criterion for the relationship between separation, compression and shear friction between blocks is first established:

[0070] ① Separation: When the distance between any two nodes of adjacent entity elements is greater than 0, separation occurs;

[0071] ② Extrusion: l = 0 and there is compressive stress between the two;

[0072] ③ Shear friction: if l = 0 and there is shear stress along the structural surface between the two;

[0073] If the adjacent blocks are separated (l>0), there is no interaction force between them, and the motion of the blocks obeys Newton's second law;

[0074] If extrusion occurs, the mechanical constitutive relationship in the normal direction of the structural surface is:

[0075] σ n = D n N max n / (N max - n)

[0076] where σ n is the compressive stress; N max is the maximum closure of the structural plane, determined by three-dimensional topography scanning; n is the closure of the structural plane, a variable; D n is the normal modulus of the structural plane, and its value is the parameter of the adjacent rock block;

[0077] If shear friction occurs, the constitutive equation of compression-shear friction for the rough structural plane:

[0078]

[0079] where σ S (σ S,p ) is the (peak) shear stress; D S is the shear stiffness; S(s p ) is the (peak) shear displacement; the parameter p is determined by σ S,p , S p , the residual shear stress σ S,r and the residual shear displacement S r ;

[0080] In step S3.4, the finite element method (FEM) and the discrete element method (DEM) are coupled for solution;

[0081] Taking the deleted zero-thickness fracture element as the boundary, regarding its internal part as a whole, according to the method in step S3.3, the nodal forces of all discrete blocks are solved by the DEM method; taking these nodal forces as the boundary conditions, according to the method in step S3.1, the nodal forces of the internal intact blocks are solved by the FEM method;

[0082] The force field data is transmitted through the shared nodes of the solid element and the fracture element. On this basis, the plastic strain and stress σ response of the solid element are calculated according to step S3.1; the damage, displacement, and stress σ c response of the fracture element are calculated according to step S3.2; if there is fracture energy and at the same time |σ - σ c | ≤ 10 -5 , then the stress, displacement, etc. responses of the solid element and the fracture element are directly transmitted to the next calculation step;

[0083] If the fracture energy indicates that the fracture element is completely fractured, at this time, the zero-thickness fracture element is deleted, and at the same time, the separation / extrusion / shear friction between the coal and rock blocks is calculated according to step S3.3;

[0084] If at the same time |σ - σ c | > 10 -5 Then, the additional nodal force is taken as |σ - σ c |, and the nodal forces of the solid elements and crack elements are recalculated according to steps S3.1 and S3.2 until |σ - σ c | ≤ 10 -5 is satisfied; then, the responses such as the stresses and displacements of the solid elements and crack elements are transmitted to the next calculation step, and the above calculation process is repeated.

[0085] Furthermore, in step S4, the coal-rock specimens are immersed in a supercritical carbon dioxide (ScCO2) preparation device for 0 - 60 days, and then for the intact coal-rock specimens, cylindrical specimens with V-shaped notches, and specimens with rough through-crack, the "triaxial pressure testing machine", "electronic universal testing machine", and "multi-functional rock triaxial dynamic shear testing machine" are respectively used to obtain the mechanical parameters of each specimen under different ScCO2 immersion times;

[0086] Based on the above experimental results, 1 field variable is set under the "Property" module, and the field variable value variables = the number of days the specimen is immersed in ScCO2. At the same time, the corresponding parameters such as "elastic modulus", "Poisson's ratio", "tensile / shear peak stress", and "frictional displacement" are set; according to the experimental conditions or engineering conditions, the corresponding stress boundaries, displacement boundaries, initial stress fields, etc. are set under the "Load" module.

[0087] Furthermore, in step S5, in the "verify mesh" of the "Mesh" module, the minimum stable increment step of the solid elements is checked to determine the minimum increment time step of the discrete element method (DEM);

[0088] The variables at the failure of the crack elements and the activation of the contact pairs are set under the Element type in the "Mesh" module;

[0089] The variables to be output are set in the create field output under the "Step" module, including stress, strain, plastic strain, damage variable, fracture energy, frictional stress, frictional displacement, etc.;

[0090] On this basis, the numerical calculation is submitted.

[0091] The present invention has the following beneficial effects:

[0092] 1) The present invention provides an effective simulation tool for the sealing performance evolution of caprocks under the corrosion of CO2 and the coupling action of complex stresses in the condition of CO2 geological storage. By using the mechanical constitutive relationship provided by the present invention, as well as the method for obtaining the mechanical parameters of coal and rock after being soaked in supercritical CO2 (ScCO2) and the parameter input method, it is possible to achieve the numerical simulation of the whole process of elastic-plastic deformation → damage → fracture → separation / extrusion / shear friction between blocks of quasi-brittle materials such as coal and rock under load conditions at any ScCO2 corrosion time.

[0093] 2) This method deletes crack elements according to the fracture energy. On this basis, the spatial regions for discrete element calculation with explicit time integration and finite element calculation with implicit time integration are demarcated, which can greatly improve the convergence of numerical solutions and accelerate the calculation efficiency at the same time.

[0094] 3) The global embedded crack element method can realize the automatic grouping of natural large-scale cracks and micro-cracks inside coal and rock blocks. Combining with the finite discrete element FDEM method, it can ultimately realize the arbitrary expansion and random fragmentation of cracks in coal and rock under the corrosion of CO2 and the coupling action of complex stresses.

[0095] 4) The present invention can be used to analyze the laws of crack initiation, growth, and propagation evolution of CO2 reservoirs and caprocks under the influence of different CO2 injection pressures and different geological parameters, so as to optimize the CO2 geological storage location and injection pressure and achieve the long-term safe storage of CO2. BRIEF DESCRIPTION OF THE DRAWINGS

[0096] Other features, purposes, and advantages of the present invention will become more obvious by reading the detailed description of non-limiting embodiments with reference to the following drawings:

[0097] Figure 1 It is a schematic diagram of the numerical simulation method flow.

[0098] Figure 2 It is the experimental result of CT scanning of coal and rock.

[0099] Figure 3 It is the numerical calculation model and the distribution of through and non-through cracks inside it;

[0100] In the figure: 1 - crack element at the main crack; 2 - crack element between solid elements.

[0101] Figure 4 It is a schematic diagram of the mechanical response of the mixed fracture mode.

[0102] Figure 5 It is the mechanical experiment diagram of three-point bending and through shear fracture.

[0103] In the figure: (a) Type I fracture specimen and its dimensions; (b) Type II fracture specimen and its dimensions; (c) Tensile damage evolution curve; (d) Shear damage evolution curve.

[0104] Figure 6 It is a three-dimensional morphology of a rough section and its direct shear test diagram.

[0105] Figure 7 It is the numerical calculation result of supercritical CO2 immersion;

[0106] In the figure: (a) is the numerical calculation result of supercritical CO2 immersion for 0 d; (b) is the numerical calculation result of supercritical CO2 immersion for 15 d; (c) is the numerical calculation result of supercritical CO2 immersion for 60 d. Specific implementation manners

[0107] The present invention will be described in detail below in conjunction with specific embodiments. The following embodiments will help those skilled in the art to further understand the present invention, but do not limit the present invention in any form. It should be noted that for those of ordinary skill in the art, without departing from the concept of the present invention, several modifications and improvements can still be made. These all belong to the protection scope of the present invention.

[0108] As Figure 1 shown, a simulation method for deformation and fragmentation of quasi-brittle materials under the action of supercritical CO2 according to an embodiment of the present invention includes the following steps:

[0109] S1. Statistically analyze the distribution of natural fractures in coal and rock masses, determine geometric parameters, and write a fracture generation program; specifically, use the CT scanning method to determine the distribution of main fractures in coal and rock, count the frequency of fractures appearing at different azimuth angles, count the average spacing and its variance of adjacent fractures; on this basis, compile an automatic main fracture generation program, specifically: according to the average spacing and its variance of the main fractures determined by scanning, randomly generate discrete points in MATLAB, determine the coordinates of the discrete points, connect the scattered points into Delaunay triangles, enclose the Voronoi polyhedron by the perpendicular bisectors of the adjacent sides of the triangles, generate multiple groups of Voronoi polyhedrons, and then select a set of Voronoi polyhedrons that are close to the azimuth angle results of the main fractures obtained by CT scanning; simulate the spatial distribution of the main fractures in coal and rock by each face of the above Voronoi polyhedron;

[0110] S2. Establish a numerical calculation model, divide the entity units, and embed fracture units at the junction of natural fractures and entity units and automatically generate a set;

[0111] S2.1. Establish a numerical calculation model and divide the entity unit mesh;

[0112] On the basis of step S1, further determine the shape and size of the research object according to the purpose of numerical simulation, and use tools such as "cutting" and "rounding" to cut the Voronoi polyhedron model into the required numerical model;

[0113] Before mesh generation, the maximum element size L is determined according to the following formula:

[0114] L = 9πEF e / [32(1 - μ 2 )σ p

[0115] where E is the elastic modulus; μ is the Poisson's ratio; F e is the fracture energy; σ p is the tensile strength;

[0116] Use the "Seed edge" in the Mesh module to seed the numerical model with a seed spacing of L, and then use the "Seed part" command to divide the solid element mesh of the numerical model;

[0117] S2.2. Embed crack elements at the intersections of natural cracks and solid elements and automatically divide the sets.

[0118] (1) Update the node numbers of the solid elements;

[0119] After mesh generation, generate an ".inp file" through CAE. Copy all the numbers below the keywords "*nset" and "*elset" and above "*assembly" in the file into Excel, and find the maximum node number n max , the maximum element number e max ;

[0120] Keep all node coordinates unchanged. Duplicate the n i (0 < i ≤ n max )th node a times, and at the same time keep the node number n i,min of the element with the smallest element number unchanged. Increase the node numbers of the remaining elements to [(a - 1) × 10 + n i,min , where a is the number of elements sharing a certain node;

[0121] (2) Create the node numbers and element numbers of the crack elements;

[0122] Taking the quadrilateral solid element as an example, denote the node numbers after the above re - numbering as N 01 , N 02 , N 03 , N 04 ; The numbers of the adjacent elements sharing the n 0i th node are e 01 , e 02 , e 03 , e 04 . Combine the new node numbers and the corresponding elements into an ordered array:

[0123] ​Array = (e 01 , e 02 , e 03 , e 04 , N 01 , N 02 , N 03 , N 04 )

[0124] Reverse the order of N 01 ~N 04 , and the node numbers of the crack elements between the quadrilateral solid elements e 01 , e 02 , e 03 , e 04 can be obtained;

[0125] The smallest number of the crack element is determined by the integer power of 10 obtained by adding 1 to the digit where the largest number of the solid element is located; for example, if the largest node number of the solid element is 756, the smallest number of the crack element is 1000.

[0126] (3) Embed crack elements at the main fissures;

[0127] Add the keyword "*Elset = cohelemAll" after the keyword "*Part" and before the keyword "*Assembly" in the ".inp file"; copy the corresponding e i and N i numbers in the order that the first column is the crack element number and the 2nd - 5th columns are the node numbers, so as to embed zero - thickness crack elements at all main fissures and all boundaries of the solid elements;

[0128] (4) Automatically divide the crack element set;

[0129] Copy the crack element number e i and its node coordinates N i (x i , y i , z i ) to Excel, and use the inner - product method to solve the direction vector d coh =(d1, d2, d3) of all crack elements, and the direction vectors d ip , if the two are collinear, then calculate the distance from any point on the fracture element to the Voronoi surface by the method of spatial vectors. If it is 0, it is determined that the fracture element is located on the surface of the Voronoi polyhedron; copy the numbers of all fracture elements and node numbers that meet the above three conditions into the.inp file, and add the keyword "*Elset=cohElem01" after the keyword "*Part" and before the keyword "*Assembly", so as to form a set named "cohelem01" of fracture elements on all main fissures;

[0130] On this basis, obtain the set "cohelem02" of fracture elements at the boundaries of all solid elements through Boolean operations.

[0131] S3. Sequentially construct the mechanical constitutive relations reflecting the plastic deformation of coal-rock matrix, the mixed-mode fracture of non-penetrating fractures, and the separation / extrusion / compressive-shear friction between blocks after the complete fracture of coal-rock; compile the finite element-discrete element FDEM calculation program to realize the whole process of coal-rock plastic deformation-mixed-mode fracture-crushing friction;

[0132] S3.1. Elastic-plastic deformation of intact coal-rock mass;

[0133] In the elastic-plastic deformation stage, each element obeys the generalized Hooke's law:

[0134] σ = E0:(ε - ε p )

[0135] where σ is the stress, E0 is the elastic stiffness matrix, and ε is the total strain. The plastic strain is determined by the yield function and the plastic potential function. The expression of the yield function F is:

[0136]

[0137] where F is the yield function, <g>is the Macaulay symbol, p = trace(σ) / 3, γ = 3(1 - K c ) / (2K c - 1), q = [2(S:S) / 3] 1 / 2 , S = σ + pI where σ i,max (i = 1, 2, 3) are the maximum effective principal stresses, is the ratio of the biaxial to uniaxial compressive yield stresses, σ t and σ c (ε p ) are the tensile and compressive stresses respectively; K c is the parameter controlling the shape of the yield surface in the deviatoric plane, and the value for rock-like materials is 0.667;

[0138] G = [(δσ t tanψ) 2 + q 2 - ptanψ

[0139] where δ is the cusp curvature of the meridian of the plastic potential function in the tensile section in the q - p plane, generally taking a value of 0.1; ψ is the dilation angle;

[0140] S3.2, Mixed - mode fracture of non - penetrating cracks;

[0141] The rock fracture process is divided into two stages: elastic deformation and mixed - mode fracture;

[0142] The constitutive relationship in the elastic stage is:

[0143] σ c = D 0,c ε c

[0144] Once the following conditions are reached, the ductile fracture process begins:

[0145]

[0146] where D 0,c is the elastic stiffness matrix, σ c,n , σ c,s , σ c,t are the normal and two tangential stresses (peak stresses), ε c is the strain;

[0147] To derive the constitutive equation for mixed - mode ductile fracture, first establish the fracture energies G c,n , G c,S (G c,S = G c,s + G c,t , G c,s and G c,t are the fracture energies in two tangential directions), G c,m The expression is:

[0148]

[0149] where j represents the n, S, m (i.e., tensile, shear, mixed-mode) fracture modes;

[0150] The expression for the fracture energy mixing ratio is:

[0151] ξ = G c,S / (G c,n + G c,S )

[0152] The criterion for complete fracture of coal and rock is:

[0153]

[0154] where χ are the tensile and shear fracture energies and material parameters at complete fracture, respectively;

[0155] Under the tensile, shear, and mixed fracture modes, the constitutive equations of each element are:

[0156]

[0157] where d c,j , D c,j , σ c,j , S c,j are the damage variable, elastic modulus, stress, plastic displacement, and total displacement, respectively, and there are:

[0158]

[0159] where S c,n and S c,S are the displacements under pure tensile and shear conditions, respectively;

[0160] Assume that the plastic displacement under the mixed fracture mode is equal to the shear plastic displacement, i.e., there is

[0161]

[0162] where the unknown parameter d c,j is determined through the following steps: ① Conduct three-point bending tests and loading-unloading direct shear tests to obtain the σ c,m -S c,m curve, σ c,n -S c,n curve, σ c,S -S c,S Curve, calculate the corresponding damage evolution curve and fracture energy ② The "coal-rock complete fracture criterion expression" was used to analyze the multiple groups of experimental The data are fitted to obtain the fracture energy mixing ratio ξ under the same material parameter χ; ③ On this basis, the "peak stress equation" is combined to obtain the peak load under a certain ξ The corresponding σ c,n With σ c,S The peak displacement component S is calculated according to the elastic expression c,n ,S c,S ; ④ Combine the "fracture energy expression" and "constitutive equation" to obtain any σ c,m The corresponding σ c,n With σ c,S Component and the corresponding displacement component S c,n ,S c,S ⑤Replace the stress and displacement in the mixed mode with the mechanical parameters of pure tension and pure shear, and calculate σ by the "constitutive equation" c,n , σ c,S and S c,n ,S c,S Plastic strain and tensile / shear damage variables under the conditions; ⑦ Calculated by the "constitutive equation" to obtain d under any mixing ratio c,m -S c,m Relationship; thereby determining the mixed ductile fracture constitutive equation;

[0163] S3.3, separation / compression / shear friction of through-fractures or blocks;

[0164] Once the fracture energy satisfies the fracture energy expression, the rock is completely fractured, and the zero-thickness crack unit on the fracture surface is deleted in the numerical calculation to eliminate its mechanical influence;

[0165] On this basis, the criterion for the relationship between separation, compression and shear friction between blocks is first established:

[0166] ① Separation: When the distance between any two nodes of adjacent entity elements is greater than 0, separation occurs;

[0167] ② Extrusion: l = 0 and there is compressive stress between the two;

[0168] ③ Shear friction: if l = 0 and there is shear stress along the structural surface between the two;

[0169] If the adjacent blocks are separated (l>0), there is no interaction force between them, and the motion of the blocks obeys Newton's second law;

[0170] If extrusion occurs, the mechanical constitutive relationship in the normal direction of the structural surface is:

[0171] σ n = D n N max n / (N max - n)

[0172] where σ n is the compressive stress; N max is the maximum closure of the structural plane, determined by three-dimensional topography scanning; n is the closure of the structural plane, a variable; D n is the normal modulus of the structural plane, and its value is the parameter of the adjacent rock block;

[0173] If shear friction occurs, the constitutive equation of compression-shear friction for rough structural planes:

[0174]

[0175] where σ S (σ S,p ) is the (peak) shear stress; D S is the shear stiffness; S(s p ) is the (peak) shear displacement; the parameter p is determined by σ S,p , S p , the residual shear stress σ S,r and the residual shear displacement S r ;

[0176] In step S3.4, the finite element FEM and the discrete element DEM are coupled for solution;

[0177] Taking the deleted zero-thickness crack element as the boundary, regarding its internal part as a whole, according to the method in step S3.3, the nodal forces of all discrete blocks are solved by the DEM method; taking these nodal forces as the boundary conditions, according to the method in step S3.1, the nodal forces of the internal intact blocks are solved by the FEM method;

[0178] The force field data is transmitted through the shared nodes of the solid element and the crack element. On this basis, according to step S3.1, the plastic strain and stress σ response of the solid element are calculated; according to step S3.2, the damage, displacement, and stress σ c response of the crack element are calculated; if there is fracture energy and at the same time |σ - σ c | ≤ 10 -5 , then directly transmit the responses such as the stress and displacement of the solid element and the crack element to the next calculation step;

[0179] If the fracture energy indicates that the crack element is completely fractured, at this time, the zero-thickness crack element is deleted, and at the same time, the separation / extrusion / shear friction between the coal and rock blocks is calculated according to step S3.3;

[0180] If at the same time |σ - σ c | > 10 -5 , the additional nodal force is taken as |σ - σ c |, and the nodal forces of the solid elements and crack elements are recalculated according to steps S3.1 and S3.2 until |σ - σ c | ≤ 10 -5 holds; then the responses such as the stresses and displacements of the solid elements and crack elements are transferred to the next calculation step, and the above calculation process is repeated

[0181] S4. Experimentally determine the evolution law of coal and rock mechanical parameters with the action time of supercritical CO2, input material parameters and set initial conditions; specifically, soak the coal and rock specimens in a supercritical carbon dioxide (ScCO2) preparation device for 0 - 60 days, and then use a "triaxial pressure testing machine", an "electronic universal testing machine", and a "multi-functional rock triaxial dynamic shear testing machine" for intact coal and rock specimens, cylindrical specimens with V-shaped notches, and specimens with rough through cracks respectively to obtain the mechanical parameters of each specimen under different ScCO2 soaking times; based on the above experimental results, set 1 field variable under the "Property" module, and the field variable value variables = the soaking days of the specimen in ScCO2. At the same time, set corresponding parameters such as "elastic modulus", "Poisson's ratio", "tensile / shear peak stress", and "frictional displacement"; according to the experimental conditions or engineering conditions, set the corresponding stress boundary, displacement boundary, initial stress field, etc. under the "Load" module;

[0182] S5. Set the increment step size, output parameters, and state variables at the moment of crack element failure and contact pair activation, and calculate the model responses before and after element failure through finite element and discrete element methods respectively; in the "verify mesh" of the "Mesh" module, check the minimum stable increment step of the solid elements, and thus determine the minimum increment time step of the discrete element DEM;

[0183] Set the variables at the moment of crack element failure and contact pair activation under the Element type in the "Mesh" module;

[0184] Set the variables to be output in the create field output under the "Step" module, including stress, strain, plastic strain, damage variable, fracture energy, frictional stress, frictional displacement, etc.;

[0185] On this basis, submit the numerical calculation.

[0186] S6. Analyze the deformation - fragmentation results of the loaded coal and rock specimens under different action times of supercritical CO2.

[0187] Example 1

[0188] S1. Process the coal and rock specimens from the Jincheng mining area into specimens with dimensions of 100 mm × 100 mm × 50 mm. Obtain the distribution of natural fractures inside them through CT scanning, as Figure 2 shown. Statistically analyze geometric parameters such as the spacing of natural fractures and the angle with the short side direction, and compile a program for automatically generating natural fractures.

[0189] S2. In the "Part module", generate multiple Voronoi polyhedra according to the results of Step S1. Based on the CT scanning results, select the closest Voronoi polyhedron, and establish a numerical calculation model of a cylinder with a diameter of 50 mm and a height of 100 mm through a cutting method. The maximum size of the fracture element is calculated to be 6 mm. Under the "Mesh module", globally divide the numerical model into tetrahedral solid element meshes. Then, through the method of globally embedding fracture elements, embed three-dimensional 6-node fracture elements at the natural fractures and the boundaries of solid elements in the numerical model respectively, and form fracture element sets of "cohelem01" and "cohelem02" respectively, as Figure 3 shown.

[0190] S3. According to the mechanical constitutive relations of plastic deformation of coal and rock matrix, mixed-mode fracture of non-penetrating fractures, and separation / extrusion / compressive-shear friction between blocks after complete fracture of coal and rock, compile the corresponding Fortran calculation program for the finite element-discrete element FDEM to realize the whole process of plastic deformation - mixed-mode fracture - fragmentation friction of coal and rock.

[0191] S4. Place cylindrical specimens, cylindrical specimens with V-notch, cylindrical specimens with circular-ring notch, and cubic specimens with rough fracture surfaces in a high-pressure reactor containing supercritical CO2 and let them stand for 0 d to 60 d. After taking them out, conduct uniaxial / triaxial compression experiments, three-point bending experiments, through shear experiments, shear friction tests, and three-dimensional topography scans of rough fracture surfaces respectively to obtain the elastic modulus, Poisson's ratio, total strain, plastic strain, ratio of triaxial peak stress to uniaxial peak stress, and dilatancy angle of coal and rock blocks, and assign the above parameters to solid elements; the tensile / shear / mixed-mode moduli, tensile / shear peak stresses, and material parameters of non-penetrating structural surfaces of coal and rock, and assign the above parameters to fracture elements; the natural spacing of rough fracture surfaces; parameters such as shear modulus, peak load, peak displacement, and residual displacement during compressive-shear friction, and assign the above parameters to the contact pairs activated after the failure of fracture elements. The parameters corresponding to different supercritical CO2 soaking times are shown in Table 1.

[0192] Table 1 Mechanical parameters

[0193]

[0194] S5. Generate a disc-shaped discrete rigid body component with a diameter of 60 mm under the Part module to simulate the indenter load during the experimental process. Then, copy 2 discrete rigid body components under the Assembly module, and set the center of the component to coincide with the axis of the cylindrical coal-rock specimen and be in contact with the upper and lower bottom surfaces of the cylinder respectively to simulate the "indenter - specimen - base" system. Set the loading speed of the upper rigid body to 0.002 mm / min in the Load module, and constrain all displacements of the lower rigid body.

[0195] Determine that the minimum stable increment step of the tetrahedral solid element is 10 according to verify mesh under the Mesh module -8 , and set an analysis step in the Step module with a time increment step of 10 -8 . Set the output variables as stress, strain, plastic strain, fracture energy, frictional displacement, frictional stress, state variable STATUS, etc.

[0196] S6. Create a new task under the Job module, and input the prepared FDEM subroutine in the User subroutine file tab; set the number of parallel computing cores to 90 cores in the Edit Job tab and submit the calculation.

[0197] Thus, the numerical calculation results of the action of supercritical CO2 for 0 - 60 d are obtained, as Figure 7 shown. The results show that the longer the action time of supercritical CO2, the higher the degree of coal-rock fragmentation under load and the lower the long-term sealing performance of the caprock.

[0198] Compared with the existing coal-rock failure models under load, on the one hand, the present invention takes into account the corrosion and strength reduction effects of supercritical CO2 on coal-rock; on the other hand, it considers the elastic-plastic deformation of intact coal-rock mass, the mixed-mode fracture of non-penetrating cracks, and the separation / extrusion / shear friction effects of penetrating cracks. On this basis, the whole process of coal-rock mass deformation - fragmentation is simulated by the FDEM numerical method, providing a powerful method for the long-term safe and stable geological storage of supercritical CO2.

[0199] The specific embodiments of the present invention have been described above. It should be understood that the present invention is not limited to the above specific embodiments, and those skilled in the art can make various deformations or modifications within the scope of the claims, which do not affect the essence of the present invention.< / g> < / g>

Claims

1. A simulation method for the deformation and fragmentation of quasi-brittle materials under the action of supercritical CO2, characterized in that: It includes the following steps: S1. Statistically analyze the distribution of natural fractures in coal and rock mass, determine geometric parameters, and write a fracture generation program; S2. Establish a numerical calculation model, divide solid elements, embed fracture elements at the junction of natural fractures and solid elements, and automatically generate a set; S3. Sequentially construct mechanical constitutive relations reflecting the plastic deformation of coal and rock matrix, the mixed-mode fracture of non-penetrating fractures, and the separation / extrusion / compressive-shear friction between blocks after the complete fracture of coal and rock; Compile a finite element-discrete element FDEM calculation program to realize the whole process of plastic deformation-mixed-mode fracture-crushing friction of coal and rock; S4. Experimentally determine the evolution law of mechanical parameters of coal and rock with the action time of supercritical CO2, input material parameters, and set initial conditions; S5. Set the increment step size, output parameters, and state variables at the moment of fracture element failure and contact pair activation, and calculate the model response before and after element failure through finite element and discrete element calculation units respectively; S6. Analyze the deformation and fragmentation results of loaded coal and rock specimens under the action of supercritical CO2 for different times.

2. The simulation method for the deformation and fragmentation of quasi-brittle materials under the action of supercritical CO2 according to claim 1, characterized in that: In step S1, the CT scanning method is used to determine the distribution of main fractures in coal and rock, count the frequency of fractures appearing at different azimuth angles, and count the average spacing and its variance of adjacent fractures; On this basis, a main fracture automatic generation program is compiled. Specifically: According to the average spacing and its variance of the main fractures determined by scanning, randomly generate discrete points in MATLAB, determine the coordinates of the discrete points, connect the scattered points into Delaunay triangles, and enclose Voronoi polyhedra by the perpendicular bisectors of the adjacent sides of the triangles. Generate multiple groups of Voronoi polyhedra, and then select a set of Voronoi polyhedra close to the azimuth angle results of the main fractures obtained by CT scanning; The spatial distribution of the main fractures in coal and rock is simulated by each face of the above Voronoi polyhedra.

3. The simulation method for the deformation and fragmentation of quasi-brittle materials under the action of supercritical CO2 according to claim 1, characterized in that: Step S2 includes the following steps: S2.

1. Establish a numerical calculation model and divide the solid element mesh; On the basis of step S1, further determine the shape and size of the research object according to the purpose of numerical simulation, and use "cutting" and "rounding" tools to cut the Voronoi polyhedron model into the required numerical model; Before mesh division, determine the maximum element size L according to the following formula: L = 9πEF e / [32(1 - μ 2 )σ p ​ where E is the elastic modulus; μ is the Poisson's ratio; F e is the fracture energy; σ p is the tensile strength; Use the "Seed edge" in the Mesh module to seed the numerical model, with the seed spacing being L, and then use the "Seed part" command to divide the solid element mesh of the numerical model; S2.

2. Embed fracture elements at the junction of natural fractures and solid elements, and automatically divide the set; (1) Update the node numbers of solid elements; After mesh generation, generate a ".inp file" through CAE. Copy all the numbers below the keywords "*nset" and "*elset" and above "*assembly" in the file into Excel, and find the maximum node number n max , the maximum element number e max ; Keep all node coordinates unchanged, and copy the n i (0 < i ≤ n max )-th node a times. At the same time, keep the node number n with the smallest element number unchanged, and increase the node numbers of the remaining elements to [(a - 1) × 10 + n i,min , where a is the number of elements shared by a certain node; i,min ​ (2) Create the node numbers and element numbers of fracture elements; In the quadrilateral solid element, denote the node number after updating the node numbers of the solid element as N 01 , N 02 , N 03 , N 04 ; The numbers of the adjacent elements sharing the n 0i -th node are e 01 , e 02 , e 03 , e 04 , and combine the new node numbers and the corresponding elements into an ordered array: Array=(e 01 , e 02 , e 03 , e 04 , N 01 , N 02 , N 03 , N 04 ) Reverse the order of N 01 to N 04 to obtain the node numbers of the crack elements between the quadrilateral solid elements e 01 , e 02 , e 03 , e 04 ; The smallest number of fracture elements is determined by the integer power of 10 obtained by adding 1 to the digit where the largest number of solid elements is located; (3) Embed fracture elements at the main fractures; Add the keyword "*Elset=cohelemAll" after the keyword "*Part" and before the keyword "*Assembly" in the ".inp file"; copy it into the corresponding e in the order that the first column is the crack element number and the 2nd to 5th columns are the node numbers. i and N i numbers, so as to embed zero-thickness crack elements at all major fissures and all solid element boundaries. (4) Automatically divide the fracture element set; Copy the crack element number e i and its node coordinates N i (x i , y i , z i ) into Excel, and use the inner product method to solve the direction vector d of all crack elements coh =(d1, d2, d3), and the direction vector d of each face of the Voronoi polyhedron ip . If the two are collinear, then calculate the distance from any point on the crack element to the Voronoi face through the spatial vector method. If it is 0, it is determined that the crack element is located on the face of the Voronoi polyhedron; copy the numbers of all crack elements and node numbers that meet the above three conditions into the.inp file, and add keywords after the keyword "*Part" and before "*Assembly" "*Elset=cohElem01”, so as to form a set named "cohelem01" by all the fracture elements on the main fractures; On this basis, the crack element set "cohelem02" at the boundaries of all solid elements is obtained through Boolean operations.

4. The simulation method for the deformation and fragmentation of quasi-brittle materials under the action of supercritical CO2 according to claim 1, characterized in that: The step S3 includes the following steps: S3.1, elastoplastic deformation of intact coal and rock mass; In the elastoplastic deformation stage, each element obeys the generalized Hooke's law: σ = E0:(ε - ε p ) Among them, σ is the stress, E0 is the elastic stiffness matrix, ε is the total strain, and the plastic strain is determined by the yield function and the plastic potential function. The expression of the yield function F is as follows: where F is the yield function and <g> is the Macaulay symbol, p = trace(σ) / 3, γ = 3(1 - K c ) / (2K c - 1), q = [2(S:S) / 3] 1 / 2 , S = σ + pI where σ i,max (i = 1, 2, 3) is the maximum effective principal stress, is the ratio of the biaxial to uniaxial compressive yield stress, σ t and σ c (ε p ) are the tensile and compressive stresses respectively; K c is the parameter controlling the shape of the yield surface on the deviatoric plane, and the value for rock-like materials is 0.667; G = [(δσ t tan ψ) 2 + q 2 - p tan ψ Among them, δ is the cusp curvature of the meridian of the plastic potential function in the q-p plane in the tension section, generally taking a value of 0.1; ψ is the dilation angle; S3.2, non-penetrating crack mixed-mode fracture; The rock fracture process is divided into two stages: elastic deformation and mixed-mode fracture; The constitutive relationship in the elastic stage is: σ c = D 0,c ε c Once the following conditions are met, the ductile fracture process is entered: where D 0,c is the elastic stiffness matrix, and σ c,n , σ c,s , σ c,t are the normal and two tangential stresses (peak stresses), and ε c is the strain; To derive the constitutive equation for mixed-mode ductile fracture, the fracture energy G under tensile, shear, and mixed modes is first established c,n , G c,S (G c,S = G c,s + G c,t , G c,s and G c,t are the fracture energies in two tangential directions), G c,m The expression is: Among them, j represents the n, S, m (i.e., tension, shear, mixed) fracture modes; The expression of the fracture energy mixing ratio is: ξ = G c,S / (G c,n + G c,S ) The criterion for complete fracture of coal and rock is: Among them, χ is the tensile, shear fracture energy and material parameters at complete fracture, respectively; In the tension, shear, and mixed fracture modes, the constitutive equations of each element are: where d c,j , D c,j , σ c,j , S c,j are the damage variable, elastic modulus, stress, plastic displacement and total displacement respectively, and there are: Among them, S c,n and S c,S are displacements under pure tension and shear conditions, respectively. Assume that the plastic displacement in the mixed fracture mode is equal to the shear plastic displacement, that is, Among them, the unknown parameters d c,j It is determined by the following steps: ① Using three-point bending test and loading and unloading direct shear test to obtain σ under tension, shear and mixed modes c,m -S c,m Curve, σ c,n -S c,n Curve, σ c,S -S c,S Curve, calculate the corresponding damage evolution curve and fracture energy ② The "coal-rock complete fracture criterion expression" was used to analyze the multiple groups of experimental The data are fitted to obtain the fracture energy mixing ratio ξ under the same material parameter χ; ③ On this basis, the "peak stress equation" is combined to obtain the peak load under a certain ξ The corresponding σ c,n With σ c,S The peak displacement component S is calculated according to the elastic expression c,n ,S c,S ; ④ Combine the "fracture energy expression" and "constitutive equation" to obtain any σ c,m The corresponding σ c,n With σ c,S Component and the corresponding displacement component S c,n ,S c,S ⑤Replace the stress and displacement in the mixed mode with the mechanical parameters of pure tension and pure shear, and calculate σ by the "constitutive equation" c,n , σ c,S and S c,n ,S c,S Plastic strain and tensile / shear damage variables under the conditions; ⑦ Calculated by the "constitutive equation" to obtain d under any mixing ratio c,m -S c,m Relationship; thereby determining the mixed ductile fracture constitutive equation; S3.3, separation / extrusion / shear friction of penetrating cracks or blocks; Once the fracture energy satisfies the fracture energy expression, the rock is completely fractured, and the zero-thickness crack elements on the fracture surface generated by the complete fracture of the rock are deleted in the numerical calculation to eliminate their mechanical effects; On this basis, the criteria for the separation, compression, and shear friction relationships between blocks are first established: ① Separation: When the distance l between any two nodes of adjacent solid elements is > 0, separation occurs; ② Extrusion: l = 0 and there is compressive stress between them; ③ Shear friction: If l = 0 and there is shear stress along the structural plane between them; If adjacent blocks are separated (l > 0), there is no mutual force between them, and the motion of the blocks obeys Newton's second law; If extrusion occurs, the constitutive relationship in the normal direction of the structural plane is: σ n = D n N max n / (N max - n) Among them, σ n is the compressive stress; N max is the maximum closure of the structural plane, determined by three-dimensional topography scanning; n is the closure of the structural plane, a variable; D n is the normal modulus of the structural plane, and its value is the parameter of the adjacent rock block; If shear friction occurs, the constitutive equation of the compression-shear friction of the rough structural plane: where, σ S (σ S,p ) is the (peak) shear stress; D S is the shear stiffness; S(s p ) is the (peak) shear displacement; the parameter p is determined by σ S,p , S p , the residual shear stress σ S,r and the residual shear displacement S r ; In step S3.4, the finite element FEM and the discrete element DEM are coupled for solution; Taking the deleted zero-thickness crack element as the boundary, regarding its interior as a whole, according to the method of step S3.3, the nodal forces of all discrete blocks are solved by the DEM method; taking this nodal force as the boundary condition, according to the method of step S3.1, the nodal forces of the intact blocks inside are solved by the FEM method; Transfer the force field data through the shared nodes of the solid element and the crack element. On this basis, calculate the plastic strain and stress σ response of the solid element according to step S3.1; calculate the damage, displacement, and stress σ of the crack element according to step S3.2 c response; if there is fracture energy G c ≤G c F , and at the same time |σ - σ c | ≤ 10 -5 , then directly transfer the stress and displacement responses of the solid element and the crack element to the next calculation step; If the fracture energy indicates that the crack element is completely fractured, the zero-thickness crack element is deleted at this time, and the separation / extrusion / shear friction between coal and rock blocks is calculated according to step S3.3; If at the same time |σ - σ c | > 10 -5 , then the additional nodal force is taken as |σ - σ c |, and the nodal forces of the solid elements and crack elements are recalculated according to steps S3.1 and S3.2 until |σ - σ c | ≤ 10 -5 holds; then the stresses and displacement responses of the solid elements and crack elements are transferred to the next calculation step, and the above calculation process is repeated.

5. A simulation method for the deformation and fragmentation of quasi-brittle materials under the action of supercritical CO2 as claimed in claim 1, characterized in that: In the step S4, the coal and rock specimens are immersed in a supercritical carbon dioxide (ScCO2) preparation device for 0 - 60 d, and then for the intact coal and rock specimens, cylindrical specimens with V-notch, and specimens with rough penetrating cracks, the "triaxial pressure testing machine", "electronic universal testing machine", and "multi-functional rock triaxial dynamic shear testing machine" are respectively used to obtain the mechanical parameters of each specimen under different ScCO2 immersion times; Set 1 field variable in the "Property" module, and the field variable value variables = the number of days the specimen is immersed in ScCO2, and at the same time set the corresponding "elastic modulus", "Poisson's ratio", "tensile / shear peak stress", and "friction displacement"; according to the experimental conditions or engineering conditions, set the corresponding stress boundary, displacement boundary, and initial stress field in the "Load" module.

6. A simulation method for the deformation and fragmentation of quasi-brittle materials under the action of supercritical CO2 as claimed in claim 1, characterized in that: In step S5, in "verify mesh" of the "Mesh" module, check the minimum stable increment step of the solid elements to determine the minimum increment time step of the discrete element DEM. Under the "Element type" of the "Mesh" module, set the variables when the crack element fails and the contact pair is activated. Under the "create field output" of the "Step" module, set the variables to be output, including stress, strain, plastic strain, damage variable, fracture energy, frictional stress, and frictional displacement. On this basis, submit the numerical calculation.

Images