Rock nanoindentation numerical simulation method based on pfc discrete element

Abstract

The application discloses a rock nanoindentation numerical simulation method based on PFC discrete elements, and relates to unconventional natural gas reservoir fracturing technology. The method comprises the following steps: based on the XRD-SEM-EDS mineral spatial distribution and the micro-parameters of single mineral nanoindentation experimental data of a rock sample, a load-displacement curve is dynamically loaded and iteratively matched; combined with the mineral spatial distribution data, a multi-mineral heterogeneous model is constructed in stages and dynamically loaded and simulated; based on the crack number development, force chain distribution and particle displacement vector characteristics, the regulation and control mechanism of the mineral interface coupling effect and the pore stress concentration on the crack propagation path is quantitatively analyzed. The application can analyze the heterogeneous evolution mechanism of the force chain development and crack propagation in the loading, holding and unloading stages, reveal the dynamic mechanism that the mineral components and pore structure cooperatively regulate and control the crack propagation path and nonlinear load response, and provide a solution path for rock micro-damage mechanism prediction and heterogeneous material discrete element numerical simulation modeling.

Description

Technical Field

[0001] This invention relates to the field of unconventional natural gas reservoir fracturing technology, specifically to a numerical simulation method for rock nanoindentation based on PFC discrete element method. Background Technology

[0002] In the field of unconventional natural gas utilization, the strategic importance of unconventional natural gas as a clean energy source is becoming increasingly prominent. Traditional hydraulic fracturing technology relies on the construction of macroscopic fracture networks, but due to the complex mechanisms of the microscopic heterogeneity of unconventional natural gas reservoirs, existing theoretical systems are unable to effectively reveal the cross-scale correlation between mineral composition, pore structure, and crack propagation.

[0003] Current research primarily employs nanoindentation technology to characterize microscopic mechanical parameters, supplemented by continuous media simulation methods such as the finite element method (FEM). Nanoindentation uses nanoscale probes to penetrate the rock surface and inverts microscopic mechanical parameters (such as elastic modulus and hardness) through load-displacement curves. The "pop-in phenomenon" (sudden load drop) observed in experiments corresponds to fracture at mineral particle interfaces, while the "elbow feature" reflects the creep characteristics of clay minerals. Traditional methods such as the finite element method (FEM) assume a continuous medium, making it difficult to characterize contact failures between mineral particles (such as frictional sliding and cementation fracture). Therefore, they face theoretical bottlenecks in analyzing contact failures and crack dynamic evolution paths in multiphase media. Especially for discontinuous mechanical behaviors such as weak cementation effects at mineral interfaces and pore-induced stress concentration, traditional numerical methods struggle to achieve high-precision simulations of multiphase coupling effects.

[0004] In recent years, the Particle Flow Code (PFC) method has demonstrated groundbreaking advantages. Compared to continuous medium methods, PFC can accurately characterize the contact constitutive relationships between rock and mineral phases by constructing multi-scale particle systems. Its explicit dynamic algorithm iteratively calculates particle trajectories using Newton's laws, dynamically capturing the entire process of force chain reconstruction and crack network evolution during loading and unloading. Studies have confirmed that PFC simulations have significant advantages in reproducing key features of nanoindentation experiments, such as the pop-in phenomenon, elbow, and creep displacement response. However, the traditional Grain-Based Model (GBM) model discretizes rocks into aggregates of contacting mineral particles and simulates rock failure processes by defining parameters such as interparticle bonding strength and friction coefficient. This method suffers from technical bottlenecks such as mineral particle size homogenization, and past nanoindentation PFC simulations have exhibited non-physical oscillations, limiting the reliability of the simulation results.

[0005] Patent document CN116911040A discloses a method for calibrating discrete element parameters of granite based on nanoindentation testing. This document constructs a heterogeneous model using Voronoi polygons and combines nanoindentation data to calibrate mesoscopic parameters, improving the scientific accuracy of the parameter calibration. However, its use of homogenized Voronoi polygon modeling technology fails to consider the spatial differentiation characteristics of multiple mineral phases in shale, leading to deviations in mineral grain size distribution from the actual geological structure and hindering high-precision simulation of the microscopic damage mechanism of heterogeneous shale. Summary of the Invention

[0006] To overcome the limitations of existing technologies, this invention proposes a numerical simulation method for rock nanoindentation based on PFC discrete element method. A heterogeneous model system containing pore defects is constructed using a force servo control dynamic loading algorithm and multi-mineral phase separation modeling technology. By integrating XRD-SEM-EDS multimodal mineral distribution data, the spatial heterogeneity of the main mineral mechanical properties of the rock is accurately characterized, revealing the evolution law of crack network and the nonlinear load response characteristics synergistically regulated by mineral interface coupling effects and pore structure.

[0007] The technical solution adopted in this invention is as follows:

[0008] A numerical simulation method for rock nanoindentation based on PFC discrete element method includes the following steps:

[0009] Step 1: Based on the spatial distribution of minerals in the rock by XRD-SEM-EDS and the experimental results of single mineral nanoindentation, the microscopic parameters of the single mineral homogeneous model of the rock are calibrated. The dynamic loading algorithm of force servo control is used for loading. The load-displacement curve is iteratively matched by trial and error method to achieve a high degree of agreement between the simulated and experimental load-displacement curves.

[0010] Step 2: Using the calibrated microscopic parameters and XRD-SEM-EDS mineral spatial distribution data, a multi-mineral heterogeneous model is constructed using a multi-mineral phase modeling method, and a force servo control dynamic loading algorithm is used for dynamic loading simulation.

[0011] Step 3: Monitor the loading process data, and quantitatively analyze the control mechanism of mineral interface coupling effect and pore stress concentration on crack propagation path based on crack number development, force chain distribution and particle displacement vector characteristics.

[0012] Furthermore, step 1 specifically includes:

[0013] Step 1.1: Based on the joint characterization results of mineral spatial distribution by XRD-SEM-EDS, the indentation point data of large-particle minerals in the nanoindentation experiment are selected as the basis for calibration of single-mineral homogeneous model parameters. Based on the load-displacement curve, the elastic modulus and hardness of the mineral sample can be calculated.

[0014] Step 1.2: Establish a PFC 2D single-mineral homogeneous model, optimize the size and total number of particles, balance the stress through servo control and establish a parallel cementation system, including four steps: sample preparation, pre-compression, cementation, and pressure head.

[0015] Step 1.3: Loading is performed using a force servo control dynamic loading algorithm. Loading and unloading are controlled separately, and the process follows unidirectional control logic to eliminate reverse loading disturbances. The Y-axis data source is switched from real-time load to target load, and the target load is recorded synchronously to reproduce the quasi-static response of the test.

[0016] Step 1.4: Adjust the micro-parameters (including particle size, modulus, cementation strength, etc.) in reverse through trial and error iteration until the error between the simulated and experimental load-displacement curves is less than a set percentage (e.g., ≤10%), and the simulated and experimental indentation failure characteristics are similar, that is, the crack propagation mode matches the experiment, and obtain the micro-parameters of the single mineral homogeneous model.

[0017] Furthermore, the specific process of step 2 is as follows:

[0018] Step 2.1: Based on the XRD-SEM-EDS joint characterization data, use AutoCAD to perform geometric topological simplification of the spatial distribution structure of mineral components.

[0019] Step 2.2: Construct five types of feature models, including clay macroporosity, fissure defects, blocky distribution, layered distribution, and mosaic distribution models. Here, certain modeling principles must be followed, taking into account factors such as model size, porosity defects, and mineral distribution, and selecting appropriate regions as the basis for modeling.

[0020] Step 2.3: Use the multiphase separation-recombination modeling method to establish a porous multi-mineral heterogeneous discrete model.

[0021] Step 2.4: Import the single mineral micro parameters obtained from the calibration in Step 1, use a force servo control head to dynamically load the material, and obtain the simulation experiment results.

[0022] Furthermore, in step 3, during the loading process, the load-displacement curve of the model, the number, location, inclination angle, and frequency of different types of cracks, the size and evolution law of force chains, and the particle displacement vector are monitored and analyzed.

[0023] The beneficial technical effects brought about by this invention are as follows.

[0024] This invention proposes a numerical simulation method for rock nanoindentation based on PFC discrete element method, and constructs a three-dimensional quantitative analysis system for the synergistic effect of mineral spatial distribution, crack dynamic evolution and pore structure, which breaks through the limitations of traditional two-dimensional experimental observation.

[0025] This invention utilizes a force servo control dynamic loading algorithm combined with XRD-SEM-EDS multimodal mineral characterization data for multi-mineral phase modeling, constructing a discrete element model (DEM) of nanoindentation in heterogeneous rocks with porosity defects. This overcomes the limitations of the traditional GBM modeling method's assumption of mineral grain size homogeneity, accurately revealing the regulatory mechanism of mineral grain size spatial heterogeneity and pore distribution on crack initiation and propagation, achieving high-precision simulation of the microscopic damage mechanism of heterogeneous rocks. Furthermore, the force servo control dynamic loading algorithm effectively suppresses the abnormal oscillation phenomenon of the load-displacement curve in previous nanoindentation PFC discrete element numerical simulations, enabling high-precision reproduction of nanoindentation experiments.

[0026] This invention can analyze the heterogeneous evolution mechanism of force chain development and crack propagation at each stage of loading, holding and unloading, and can reveal the dynamic mechanism of the synergistic regulation of crack propagation path and nonlinear load response by mineral composition and pore structure, providing a solution for predicting rock micro-damage mechanism and modeling discrete element numerical simulation of heterogeneous materials. Attached Figure Description

[0027] Figure 1 This is a flowchart of the simulation method described in this invention;

[0028] Figure 2 This is a flowchart for constructing a homogeneous model of a single mineral;

[0029] Figure 3 This is a flowchart for establishing a PFC 2D single-mineral homogeneous model;

[0030] Figure 4 This is a flowchart for constructing a multi-mineral heterogeneous model;

[0031] Figure 5 This is a flowchart for establishing a PFC 2D multi-mineral heterogeneous model;

[0032] Figure 6 It shows the load-displacement curves and failure mode diagrams of a single-mineral homogeneous quartz model in numerical simulation.

[0033] Figure 7 This is a diagram of damage analysis in a numerical simulation of a multi-mineral heterogeneous model. Detailed Implementation

[0034] The present invention will be further described in detail below with reference to the accompanying drawings and specific embodiments:

[0035] This embodiment proposes a numerical simulation method for rock nanoindentation based on PFC discrete element method. Through cross-scale coupling analysis of nanoindentation experiments and discrete element simulation, combined with XRD-SEM-EDS multimodal mineral characterization data, a multi-mineral heterogeneous model considering pore defects is constructed. Using a force servo control dynamic loading algorithm, the method elucidates the regulation mechanism of rock micro-heterogeneity on mechanical behavior and reveals the intrinsic relationship between the mechanical properties of mineral components and the dynamic evolution of cracks.

[0036] See Figure 1 This embodiment uses a numerical simulation method for rock nanoindentation based on PFC discrete element method, taking shale material as an example, and mainly includes the following three steps:

[0037] Step 1: Based on the XRD-SEM-EDS and single-mineral nanoindentation experimental results of shale samples, the microscopic parameters of quartz, calcite, and clay were calibrated. Dynamic loading was performed using a force servo control dynamic loading algorithm. Through trial and error iterative matching, the simulated and experimental load-displacement curves were highly consistent, and a homogeneous single-mineral nanoindentation model of shale was constructed.

[0038] Step 2: Using the microscopic parameters calibrated in Step 1, combined with XRD-SEM-EDS mineral spatial distribution data, multi-mineral phase modeling is performed to construct a shale nanoindentation multi-mineral heterogeneous model, and dynamic loading simulation is conducted.

[0039] Step 3: Based on crack quantity development, force chain distribution, and particle displacement vector characteristics, quantitatively analyze the control mechanism of mineral interface coupling effect and pore stress concentration on crack propagation path. During loading, monitor and analyze data such as the model's load-displacement curve, the number, location, inclination angle, and frequency of different types of cracks, force chain size and evolution law, and particle displacement vector.

[0040] like Figure 2 As shown, the single-mineral homogeneous model construction process in step 1 specifically includes the following steps:

[0041] Step 1.1: Based on the combined XRD-SEM-EDS characterization results, the spatial distribution of three main minerals—quartz, calcite, and clay—was determined. Data from the indentation points of large-particle minerals in the nanoindentation experiment were selected as the basis for calibrating the parameters of the single-mineral homogeneous model. Based on the load-displacement curve, the elastic modulus and hardness of the mineral sample can be calculated. The specific equations are as follows:

[0042]

[0043] Where H is hardness, Pmax is maximum load, Ac is contact area, Er is reduction modulus, S = dP / dh is contact stiffness, α = 1.034 is the shape factor of the glass indenter, E and v are the elastic modulus and Poisson's ratio of the shale sample, respectively, Ei and Vi are the elastic modulus and Poisson's ratio of the indenter, respectively, and Ei = 1140 GPa and vi = 0.07 for the glass indenter.

[0044] Step 1.2: Construct a PFC 2D single-mineral homogeneous model, optimize size and total number of particles, balance stress through servo control, and establish a parallel cementation system, including sample preparation, pre-compression, cementation, indentation, and loading. The specific process is as follows: Figure 3 As shown:

[0045] A two-dimensional discrete element model of the shale sample was established, and a 10μm equally spaced indentation array was constructed to eliminate boundary effects based on the principle of experimental comparability. While ensuring computational accuracy, the total number of particles was optimized by setting the sample aspect ratio to 0.5 (10μm×5μm), which effectively improved computational efficiency.

[0046] A servo control mechanism was used to adjust the displacement of the confining pressure baffle, achieving initial stress equilibrium of the specimen based on the particle elastic modulus threshold. The equilibrium stress threshold was determined according to the particle elastic modulus. This process allowed for a small amount of elastic deformation in the specimen, which was verified not to affect subsequent experimental results.

[0047] A linear parallel bonded system was constructed, where the interparticle contact was modeled using a parallel bond, while the particle baffle contact maintained a linear relationship. The formulas for calculating the total force and total moment are as follows:

[0048]

[0049] Among them, F c For the total force, F l For linear force, F d For damping force, The adhesive force at the interface. M c For the total torque, It represents torque.

[0050] The upper baffle was removed, and the remaining three side baffles were fixed. After the sample was balanced, a two-dimensional model of the Berkovich indenter was created. Based on previous research, the three-dimensional structure of the Berkovich indenter was equivalently transformed into a two-dimensional conical indenter with a 70.3° half-apex angle. This geometric transformation maintained a constant projected area-indentation depth ratio, ensuring the equivalence of the load-displacement response.

[0051] During the loading phase, since PFC has a built-in time stepping concept, the loading speed in the simulation should be much higher than the loading speed in the actual experiment, without affecting the authenticity of the results, in order to shorten the time required for a single simulation.

[0052] Step 1.3: Dynamic loading is performed using a force servo control dynamic loading algorithm. The specific process is as follows:

[0053] Establish time-varying target load function F target1 = K × t. K is the proportionality coefficient, representing the rate of load change. When the real-time load F of the pressure head... actual <F target At that time, accelerated downlink loading is executed, with a loading rate v1 = K × (F) target -F actual When F actual ≥F target At that time, the loading rate immediately returns to zero, and the loading time parameter is set to obtain the final load F. final .

[0054] Reaching F final Then, maintain the target load F target =F final Perform closed-loop force control and set the load holding time parameter.

[0055] Constructing the time-varying unloading function F target2 =F final -K×t. When F actual >F target2 At that time, accelerated uplink unloading is executed, and the unloading rate v2 = K × (F) actual -F target2 When F actual ≤F target2 At that time, the unloading rate immediately returns to zero. The unloading time is designed to be symmetrical with the loading time, and the final load converges to zero.

[0056] Simultaneously record target load and real-time displacement data to construct load-displacement characteristic curves.

[0057] Here, a force servo control dynamic loading algorithm is adopted, which strictly follows unidirectional control logic during loading and unloading. During loading, forward loading occurs when the real-time load is lower than the target load, and the speed becomes zero when it is higher than the target load; the unloading process is the opposite. This eliminates reverse loading disturbances and accurately reproduces the mechanical path of the nanoindentation experiment. Switching the Y-axis data source from real-time load to target load effectively suppresses signal oscillations caused by the elastic collision of discrete particles, which is more in line with the data output logic of the nanoindenter. Compared with the traditional loading mode, the advantages are: the uniform speed loading mode does not conform to the force control loading logic of the nanoindentation experimenter; in the ordinary force servo control mode, when the sample particles detach from the cementation constraint, the bidirectional movement of the indenter causes repeated collisions between the detached particles and the indenter, resulting in high-frequency fluctuations in the measured load. The force servo control dynamic loading algorithm in this invention avoids reverse motion disturbances by adopting a unidirectional force servo control strategy during loading / unloading, and replaces the Y-axis data source with the target load value instead of the real-time load measurement value, avoiding interference from discrete particle collisions and successfully eliminating abnormal fluctuations. Numerical verification based on quartz samples shows that the improved loading mode successfully reproduces the quasi-static load response characteristics of the indoor nanoindentation experiment, and the selection of the loading mode has a decisive influence on the PFC simulation of nanoindentation experiment.

[0058] Step 1.4: The main microscopic parameters, such as particle size, elastic modulus, and cementation strength, are adjusted using a trial-and-error iterative method. Through parameter inverse identification, the relative error between the simulated load-displacement curve and experimental data is made ≤10%, while the crack propagation mode closely matches the physical experiment, verifying the model's effectiveness and obtaining the microscopic parameters of the single-mineral homogeneous model. Load-displacement curve comparison and indentation failure characteristics are shown below. Figure 6 As shown.

[0059] See further Figure 4 The specific process for constructing the multi-mineral heterogeneous model in step 2 is as follows:

[0060] Step 2.1: Based on the combined XRD-SEM-EDS characterization data, AutoCAD is used to simplify the two-dimensional spatial distribution structure of mineral components. XRD analysis revealed that the shale sample contained three main minerals: quartz, calcite, and clay. The distribution of Al, Si, and Ca elements in the EDS analysis results can be used to briefly delineate the distribution areas of clay, quartz, and calcite. Combined with the SEM results, a simplified mineral component diagram can be obtained using AutoCAD.

[0061] Here, the XRD-SEM-EDS joint characterization data is a mineralogy dataset obtained through the collaborative acquisition of three analytical techniques, including XRD (X-ray diffraction), SEM (scanning electron microscopy), and EDS (energy dispersive spectroscopy). It aims to achieve a comprehensive characterization of mineral composition, morphology, and elemental distribution, breaking through the limitations of single techniques and revealing mineral heterogeneity (such as the distribution patterns of clay, quartz, and pyrite in shale) through data fusion.

[0062] Step 2.2: Construct five types of feature models: clay macropores, fissure defects, blocky distribution, layered distribution, and mosaic distribution models.

[0063] To ensure that the model reflects the heterogeneous characteristics of shale, the following modeling guidelines must be followed:

[0064] First, since the longitudinal and transverse distribution characteristics of shale material are similar at the same point, XY plane cross-sectional SEM data can be used as the basis for longitudinal cross-section modeling. Second, since the spacing between indentation points needs to be greater than three times the indentation radius to eliminate mutual influence between indentation points, the distance from the indentation center to the boundary should be ≥5μm to eliminate the boundary effect. Third, the selected area should contain ≥2 mineral phases and consider the pore distribution to meet the research needs of multi-mineral heterogeneity of shale.

[0065] Based on the above modeling criteria and the simplified mineral composition diagram obtained from XRD-SEM-EDS results, five types of feature models can be selected: clay macropore, fracture defect, massive distribution, layered distribution, and mosaic distribution model.

[0066] Step 2.3: To overcome the particle size homogenization defect in traditional GBM methods for micron-scale modeling, this invention proposes a multiphase separation-recombination modeling method to construct a multi-mineral heterogeneous model containing voids. This mainly includes phase separation construction, particle embedding, and weak cementation treatment, such as... Figure 5 As shown.

[0067] The specific method is as follows:

[0068] Based on the mineral composition of the sample, single-mineral rectangular models of quartz, calcite, and clay were constructed by phase separation. The particles in the distribution areas of each single mineral were retained, and the particles in redundant areas were deleted. Multi-phase particle mosaic merging was performed to obtain a multi-mineral heterogeneous model of shale nanoindentation, which made all parameters of quartz, calcite, and clay controllable and showed multi-scale particle size distribution characteristics.

[0069] A servo control mechanism was used to adjust the displacement of the confining pressure baffle to achieve initial stress equilibrium in the specimen. The equilibrium stress threshold was determined based on the elastic modulus of the particles. This process allowed for a small amount of elastic deformation in the specimen, which was verified not to affect subsequent experimental results.

[0070] The LPBM model was used for both intragranular cementation and intergranular cementation. The intergranular cementation strength was appropriately reduced to simulate the weak cementation characteristics between different minerals. The mineral boundary contact used a linear contact model. The pore system was selectively removed by regional particles after cementation.

[0071] Remove the upper baffle, fix the remaining three side baffles, and after the sample is balanced, perform a two-dimensional model of the Berkovich indenter. The three-dimensional structure of the Berkovich indenter is equivalently transformed into a two-dimensional conical indenter with a 70.3° half-apex angle. This geometric transformation maintains a constant projected area-to-depth ratio, ensuring the equivalence of the load-displacement response.

[0072] Step 2.4: Based on the multi-mineral heterogeneous model obtained in the previous step, import the calibration micro-parameters of each single mineral obtained in Step 1. The indenter insertion position can be controlled by methods such as left and right translation or deletion of surface region particles. Force servo control mode is adopted during the loading stage. The 8mN load insertion result of the multi-mineral heterogeneous model is shown below. Figure 7 As shown, the upper figure represents the distribution of cracks and particle displacements caused by pressure on the model, while the lower figure represents the distribution of the pressure chain.

[0073] Step 3: During the loading process, data is collected through monitoring. Based on the crack number development, force chain distribution, and particle displacement vector characteristics, the regulation mechanism of mineral interface coupling effect and pore stress concentration on crack propagation path is quantitatively analyzed.

[0074] The data that needs to be monitored and obtained includes: the model's load-displacement curves, the number, location, inclination angle, and frequency of different types of cracks, the size and evolution of force chains, and particle displacement vectors. The specific process is as follows:

[0075] Real-time monitoring of target load and indenter displacement is used to construct load-displacement characteristic curves.

[0076] Based on different contact types, servo functions are used to classify cracks and monitor their respective quantities, including intragranular and intergranular cracks, tensile cracks and shear cracks.

[0077] Based on the different types of crack generation, a servo function is used to statistically monitor the number of each type of crack generated at regular intervals from the start to the end of loading, so as to obtain the generation frequency of different types of cracks.

[0078] By utilizing the built-in monitoring functions of PFC 2D, such as ball, fracture dfn name, fracture rosette, contact force, and ball displacement, the location of crack initiation, crack inclination angle, force chain size and evolution, and particle displacement vector can be monitored in real time.

[0079] As can be seen from the above embodiments, this invention constructs a multi-mineral heterogeneous model considering pore defects by combining cross-scale coupling analysis of nanoindentation experiments and discrete element simulations with XRD-SEM-EDS multimodal mineral characterization data. It utilizes a force servo control dynamic loading algorithm to reveal the regulatory mechanism of rock microscopic heterogeneity on crack evolution. This breakthrough overcomes the limitations of traditional two-dimensional experimental observations and the bottleneck of nanoindentation numerical simulation, systematically analyzing the heterogeneous evolution law of the force chain throughout the loading-holding-unloading cycle. It establishes a quantitative mapping relationship between the mechanical properties of mineral components and crack propagation paths, elucidating the synergistic damage mechanism of pore structure-induced stress concentration and mineral interface cementation strength. It achieves accurate simulation of the three-dimensional dynamic evolution of multiphase medium crack networks, revealing the energy conversion mechanism of elastic strain energy release driving the initiation of circumferential secondary cracks. It accurately reproduces the differential deformation characteristics of quartz dominated by elasticity and clay dominated by plasticity, ensuring a high degree of agreement between the simulation results and actual nanoindentation experimental conditions. This provides theoretical basis and methodological support for optimizing fracturing parameters in unconventional natural gas reservoirs and for discrete element modeling of heterogeneous materials.

Claims

1. A numerical simulation method for rock nanoindentation based on PFC discrete element method, characterized in that, include: Step 1: Based on the XRD-SEM-EDS mineral spatial distribution and single-mineral nanoindentation experimental results of the rock sample, the microscopic parameters of the single-mineral homogeneous model are calibrated. The force servo control dynamic loading algorithm is used for loading. The loading, holding, and unloading processes follow unidirectional control logic to eliminate reverse loading disturbances. The Y-axis data source is switched from real-time load to target load. The target load is recorded synchronously to reproduce the experimental quasi-static response. The load-displacement curve is matched by trial and error iteration method to achieve a high degree of agreement between the simulated and experimental load-displacement curves. Step 2: Using the calibrated microscopic parameters and XRD-SEM-EDS mineral spatial distribution data, multi-mineral phase modeling is performed to construct a multi-mineral heterogeneous model, and a force servo control dynamic loading algorithm is used for dynamic loading simulation. Specifically, this includes: Step 2.1, based on the XRD-SEM-EDS joint characterization data, performing geometric topological simplification on the spatial distribution structure of mineral components; Step 2.2: Construct five types of heterogeneous feature models, including clay macropore, fracture defects, blocky distribution, layered distribution, and mosaic distribution models; Step 2.3: Use the multiphase separation-recombination modeling method to establish a multi-mineral heterogeneous discrete model with pores: specifically, use the multiphase separation-recombination modeling method to construct single-mineral models by phase separation, and generate a multi-mineral heterogeneous pore model by particle mosaic merging to achieve multi-scale particle size distribution parameter control; based on the particle elastic modulus threshold, use a servo control mechanism to adjust the displacement of the confining pressure baffle to achieve the initial stress balance of the sample; unify the intracrystalline cementation and intercrystalline cementation LPBM models, appropriately reduce the intercrystalline cementation strength, use a linear contact model for mineral boundary contact, and generate the pore system through selective particle deletion; delete the upper baffle, fix the remaining three side baffles, and perform two-dimensional modeling of the glass indenter after the sample is balanced; Step 2.4: Import the microscopic parameters obtained from the calibration in Step 1, use force servo control to dynamically load the indenter, and obtain the simulation experimental results; Step 3: Monitor the loading process data, and quantitatively analyze the control mechanism of mineral interface coupling effect and pore stress concentration on crack propagation path based on crack number development, force chain distribution and particle displacement vector characteristics.

2. The numerical simulation method for rock nanoindentation based on PFC discrete element method according to claim 1, characterized in that: Step 1 includes: Step 1.1: Based on the joint characterization results of mineral spatial distribution by XRD-SEM-EDS, the indentation point data of large-particle minerals in the nanoindentation experiment are selected as the basis for calibrating the micro-parameters of the homogeneous single-mineral model, and the elastic modulus and hardness of the single mineral are calibrated. Step 1.2: Establish a PFC 2D single-mineral homogeneous model; Step 1.3: Apply the load using a force servo control dynamic loading algorithm; Step 1.4: Adjust the micro parameters through trial and error iteration until the error between the simulated and experimental load-displacement curves is less than the set percentage, so that the simulated and experimental indentation failure characteristics match, and obtain the micro parameters of the single mineral homogeneous model.

3. The numerical simulation method for rock nanoindentation based on PFC discrete element method according to claim 2, characterized in that: In step 1.1, the elastic modulus and hardness of the rock sample are calculated based on the load-displacement curve, and the specific equations are as follows: (1) (2) (3) Where H is hardness, Pmax is maximum load, Ac is contact area, Er is reduction modulus, S=dP / dh is contact stiffness, α is the shape factor of the glass indenter, E and v are the elastic modulus and Poisson's ratio of the rock sample, respectively, and Ei and vi are the elastic modulus and Poisson's ratio of the indenter, respectively.

4. The numerical simulation method for rock nanoindentation based on PFC discrete element method according to claim 2, characterized in that: The construction process in step 1.2 includes: A two-dimensional discrete element model of the rock sample was established, with the aspect ratio of the sample set to 10 μm × 5 μm. A servo control mechanism is used to adjust the displacement of the confining pressure baffle to achieve initial stress balance of the specimen; Constructing a linear parallel cementation system between particles; Remove the upper baffle, fix the remaining three side baffles, and after the sample is balanced, perform a two-dimensional model of the glass indenter.

5. The numerical simulation method for rock nanoindentation based on PFC discrete element method according to claim 4, characterized in that: In step 1.2, a linear parallel bonding model is used to construct a linear parallel bonding system. The formulas for calculating the total force and total moment are as follows: (4) (5) Among them, F c For the total force, F l For linear force, F d For damping force, For the bonding force on the interface, M c For the total torque, It represents torque.

6. The numerical simulation method for rock nanoindentation based on PFC discrete element method according to claim 1, characterized in that: The modeling in step 2.2 follows three criteria: based on the principle of orthogonal anisotropy of materials, XY plane cross-sectional data is selected; the center distance of the indentation from the boundary is ≥5 μm to avoid boundary effects. The selected area must contain at least two mineral phases and take into account the porosity distribution.

7. The numerical simulation method for rock nanoindentation based on PFC discrete element method according to any one of claims 1-6, characterized in that: The monitoring in step 3 includes: monitoring the load-displacement curve of the model during loading, monitoring the number, location, tilt angle, and frequency of different types of cracks, monitoring the force chain size and evolution law, and particle displacement vector.

8. The numerical simulation method for rock nanoindentation based on PFC discrete element method according to claim 7, characterized in that: During the monitoring in step 3: Real-time monitoring of target load and indenter displacement to construct load-displacement characteristic curves; Based on different contact types, servo functions are used to classify cracks and monitor their respective quantities, including intragranular and intergranular cracks, tensile cracks and shear cracks. Based on the occurrence of different types of cracks, a servo function is used to statistically monitor the number of each type of crack generated at intervals from the start to the end of loading, so as to obtain the frequency of different types of cracks. Using the monitoring function of PFC 2D, the location of crack initiation, crack inclination angle, force chain size and evolution, and particle displacement vector are monitored in real time.

Images