A discrete element method combined with finite element method artificial crack dynamic crack width calculation method

Abstract

This invention discloses a method for calculating the dynamic width of artificial fractures using a combination of discrete element method (DEM) and finite element method (FEM). The method applies the DEM to determine the constitutive characteristic parameters of proppant particles and then uses these parameters in a finite element numerical simulation of fracture width closure, achieving efficient prediction of the dynamic width of artificial fractures during mining. This invention fully considers the mechanical characteristics of proppant particle packs, using the DEM to obtain the constitutive characteristic parameters of the proppant packs. This eliminates the reliance on triaxial experiments of the particles and removes the influence of the non-homogeneity of the proppant particle size, effectively reducing costs and improving efficiency. Compared to traditional methods, this invention considers both the bulk mechanical characteristics of proppant particle packs and the solution efficiency of engineering-scale models, making it more suitable for field applications.

Description

Technical Field

[0001] This invention belongs to the field of hydraulic fracturing technology for oil and gas field development, specifically involving a method for calculating the dynamic fracture width of artificial fractures using a combination of discrete element method and finite element method. Background Technology

[0002] Hydraulic fracturing is an important method for enhancing oil and gas well production and injection, and it is now an essential technology for the efficient development of low-permeability and shale reservoirs. However, after fracturing and commissioning of oil and gas wells, the support width of artificially filled fractures will close under the effective stress caused by the fluid-structure interaction effect in the reservoir, and the degree of closure varies with the development dynamics of the reservoir. Especially for reservoirs with high closure pressure, the closure of artificial fractures caused by stress disturbance is more obvious. However, current reservoir numerical simulations still set the fracture width as a constant value and do not consider the loss of fracture conductivity caused by fracture closure. This not only affects the accuracy of production prediction, but also hinders the selection of wells and reservoirs for repeated fracturing.

[0003] Due to limitations in geometric dimensions and time scales, the dynamic fracture width closure under long-term mining conditions is difficult to obtain experimentally. Fracture width damage obtained using fracture conductivity meters primarily originates from proppant embedding and cannot account for fracture closure under the influence of reservoir fluid-structure interaction-induced stress. Numerical simulation is an effective means to address this problem, but current dynamic fracture width closure models do not consider the constitutive characteristics of propped fractures, and the selection of parameters such as Young's modulus of propped fractures in the models lacks a basis. Since the single-particle elastic parameters of commonly used proppants (such as ceramsite, coated sand, and quartz sand) are difficult to obtain, some researchers have proposed the concept of apparent Young's modulus to characterize the constitutive behavior of proppant particle packs and have developed a method to obtain the apparent Young's modulus under closure pressure loading conditions through steel plate conductivity testing. However, the constitutive characteristics of particle packs are affected by complex factors such as particle size, particle properties, particle shape, proppant concentration, and closure pressure. Experimental simulation testing is costly, time-consuming, complex in sample preparation, and subject to stringent conditions. Currently, there is no usable apparent Young's modulus prediction chart, which restricts the field application of this method.

[0004] Given the complexity of the mechanical characteristics of granular aggregates, conventional numerical simulation methods such as the finite element method, finite difference method, and boundary element method all have limitations in analyzing the deformation behavior of granular masses. However, the discrete element method (DEM) can reveal granular mechanical phenomena that conventional mechanics theories cannot explain. It is not limited by the amount of deformation and can effectively simulate the discontinuous mechanical phenomena of granular aggregate media, making it the most mature numerical simulation method for the mechanical characteristics of granular masses. Therefore, there is an urgent need to propose a method that uses the DEM to determine the constitutive characteristic parameters of proppant particles and applies it to the finite element numerical simulation of fracture width closure, in order to achieve efficient prediction of the dynamic fracture width of artificial fractures during mining. Summary of the Invention

[0005] This invention addresses the lack of existing methods for predicting the dynamic width of artificial fractures after fracturing at the engineering scale. It proposes a method for calculating the dynamic width of artificial fractures using a combination of discrete element and finite element methods. This method is applicable to calculations under engineering scale and long-cycle mining conditions, meeting the needs of on-site fracturing design and production capacity prediction.

[0006] This invention is achieved using the following technical solution: a method for calculating the dynamic crack width of artificial cracks using a combination of discrete element method and finite element method, comprising the following steps:

[0007] Step A: Obtain the particle size distribution data of the proppant used on site, perform a three-dimensional scan on a single particle, and obtain the geometry of the single particle;

[0008] Step B: Construct a single-particle discrete element model and obtain the discrete element modeling parameters for a single particle;

[0009] Step C: Combine the discrete element modeling parameters obtained in Step B to establish a discrete element model of the proppant particle pack and obtain the elastoplastic constitutive mechanical parameters of the proppant particle pack.

[0010] Step D: Apply the elastoplastic constitutive mechanical parameters of the proppant mass to the dynamic fracture width closed finite element calculation during the artificial fracture mining process to realize the calculation of the dynamic fracture width of the artificial fracture.

[0011] Furthermore, step B is implemented in the following manner:

[0012] Based on the particle size distribution data and geometry obtained in step A, a single-particle discrete element model is constructed. A single-particle uniaxial compression experiment is conducted, and the same pressure conditions as those in the single-particle uniaxial compression experiment are applied to the single-particle discrete element model. The uniaxial compression stress-strain curve of the single particle during loading is obtained. The effective modulus, stiffness ratio, and friction coefficient of the discrete element linear model are adjusted, and the stress-strain curve obtained from the simulation is fitted with the results of the single-particle uniaxial compression experiment to determine the constitutive characteristic parameters of the single particle, thereby obtaining the discrete element modeling parameters of the single particle.

[0013] Furthermore, step C specifically includes the following steps:

[0014] Step C1: Establish a discrete element model of the proppant particle pack. The particle size in this model is set according to the particle size distribution in step A. The effective modulus, stiffness ratio, and friction coefficient of the particles are the same as the fitting results in step B.

[0015] Step C2: Apply closed pressure to the discrete element model of proppant particle packing to simulate the uniaxial compression process of proppant particles under confining pressure; record the stress-strain curves during the compression process, and obtain the elastic modulus of the proppant particles in the elastic stage, the stress yield strength in the plastic stage, and the tangential modulus parameters based on these curves.

[0016] Furthermore, step D specifically includes the following steps:

[0017] Step D1: Establish a three-dimensional model of the reservoir with artificial fractures at the engineering scale. The geometric dimensions of the fractures and reservoir in the model are set according to the field fracture monitoring results, and its elastic-plastic constitutive parameters are determined by the discrete element simulation results in step C2.

[0018] Step D2: Solve the solid mass balance equation, fluid mass balance equation, Biot effective stress equation, and Kozeny-Carman permeability response equation during reservoir exploitation using the finite element method; calculate the fracture width closure of artificial fractures under fluid-solid coupling during fractured reservoir exploitation, and analyze the main factors affecting fracture width closure characteristics.

[0019] Furthermore, in step D1, the three-dimensional model of the reservoir with artificial fractures includes a simulated wellbore, simulated fractures, and reservoir matrix. Initial values ​​and boundary conditions are assigned to the three-dimensional model of the reservoir with artificial fractures, wherein the simulated fractures are set as elasto-plastic materials, and the reservoir matrix is ​​set as an ideal elastic material, the parameters of which are obtained from the physical property parameter tests of the reservoir core; the simulated fractures and the reservoir matrix are set as a deformable and permeable interface, and the reservoir matrix interface is set as a pressure interface.

[0020] Furthermore, step D2 is specifically solved based on the following principle:

[0021] After initializing the model, the seepage field during reservoir production is first solved based on the fluid mass conservation equation to obtain the reservoir pore pressure distribution law; the effective stress change induced by pore pressure change is solved based on the Biot effective stress principle to solve the continuity equation and calculate the deformation equation of the reservoir solid phase unit; the relationship between strain and porosity and permeability is established according to the Kozeny-Carmen equation, and the reservoir porosity and permeability data are updated after each time step; after checking the model convergence, the calculation of the next time step is carried out, and after the calculation set time step is reached, the closure of the support fracture width is checked.

[0022] Compared with the prior art, the advantages and positive effects of the present invention are as follows:

[0023] The proposed discrete element method (DEM) combined with finite element method (FEM) model for predicting the dynamic width of artificial cracks fully considers the mechanical characteristics of proppant particle packs. The DEM is used to obtain the constitutive characteristic parameters of the proppant packs, eliminating reliance on triaxial experiments and removing the influence of particle size heterogeneity. This effectively reduces costs and improves efficiency. Compared to existing crack width prediction techniques, the numerical simulation method combining DEM and FEM considers both the bulk mechanical characteristics of proppant particle packs and the solution efficiency of engineering-scale models, making it more suitable for practical field applications. Attached Figure Description

[0024] Figure 1 This is a schematic diagram of the proppant particle shape, scanning three-dimensional model, and single-particle discrete element model as described in an embodiment of the present invention;

[0025] Figure 2 This is a schematic diagram showing the fitting of the uniaxial compression experiment and discrete element simulation results of the proppant particles described in the embodiments of the present invention;

[0026] Figure 3 In the diagram, (a) is a schematic diagram of the discrete element model of the proppant particle pack, and (b) is a schematic diagram of servo loading.

[0027] Figure 4 This is a schematic diagram of the stress-strain curve of the proppant particle pack as described in an embodiment of the present invention;

[0028] Figure 5 This is a schematic diagram showing the changes in elastic modulus and tangential modulus of the proppant particle pack as a function of confining pressure according to an embodiment of the present invention.

[0029] Figure 6 The following is a schematic diagram of the finite element geometric model of the reservoir with artificial fractures according to an embodiment of the present invention: (a) is a schematic diagram of the fracture geometry; (b) is a schematic diagram of the three-dimensional model of the reservoir.

[0030] Figure 7 This is a schematic diagram of the solution process for the multi-field coupling model of reservoirs with artificial fractures as described in an embodiment of the present invention;

[0031] Figure 8 The curves showing the dynamic crack width closure of the artificial crack as described in the embodiments of the present invention are as follows: (a), (b), (c), and (d) correspond to the curves for 10 days, 20 days, 30 days, and 40 days, respectively. Detailed Implementation

[0032] To better understand the above-described objects, features, and advantages of the present invention, the present invention will be further described below in conjunction with the accompanying drawings and embodiments. Many specific details are set forth in the following description to provide a thorough understanding of the present invention; however, the present invention may be practiced in other ways than those described herein, and therefore, the present invention is not limited to the specific embodiments disclosed below.

[0033] This embodiment discloses a method for calculating the dynamic fracture width of artificial fractures after hydraulic fracturing based on a combination of discrete element and finite element methods. The method determines the constitutive characteristic parameters of proppant particles using the discrete element method and applies them to the finite element numerical simulation of fracture width closure. This method is applicable to artificial fractures formed after various reservoir stimulation processes and can be used in hydraulic fracturing scheme design, well and layer selection for repeated fracturing, and post-fracturing production prediction. Specifically, it includes the following steps:

[0034] Step A: Obtain the particle size distribution data of the proppant used on site, perform a three-dimensional scan on a single particle, and obtain the geometry of the single particle;

[0035] Step B: Construct a single-particle discrete element model and obtain the discrete element modeling parameters for a single particle;

[0036] Step C: Combine the discrete element modeling parameters obtained in Step B to establish a discrete element model of the proppant particle pack and obtain the elastoplastic constitutive mechanical parameters of the proppant particle pack.

[0037] Step D: Apply the constitutive characteristic parameters of the proppant deposit to the dynamic fracture width closed finite element calculation during the artificial fracture mining process to realize the calculation of the dynamic fracture width of the artificial fracture.

[0038] To better understand the present invention, the following detailed description is provided with reference to specific examples:

[0039] I. Obtaining the particle size distribution of proppant particles:

[0040] In step A, a particle size analyzer is used to obtain data such as the particle size distribution of the proppant used on site. Specifically, commonly used proppants on site, such as quartz sand, can be used. In this embodiment, the proppant particle sizes include 8 / 16 mesh (1.18~2.36mm), 20 / 40 mesh (0.425~0.85mm), and 40 / 70 mesh (0.212~0.425mm). The particle size analyzer is used to obtain the particle size distribution of proppant particles of each group.

[0041] II. Obtaining Single-Particle Discrete Element Modeling Parameters

[0042] Based on the particle size and geometry obtained from the particle size analyzer in step A, a single-particle discrete element model is constructed. A servo method is applied to the discrete element model with the same pressure conditions as in the single-particle uniaxial compression experiment. The uniaxial compression stress-strain curve of the particle during loading is obtained. The effective modulus, stiffness ratio, and friction coefficient of the discrete element linear model are adjusted. The simulated stress-strain curve is then fitted with the experimental results to determine the constitutive characteristic parameters of the single particle. The specific steps include:

[0043] (1) Perform a three-dimensional scan of the proppant particles to obtain their geometric morphology and establish a corresponding three-dimensional model file. Based on this three-dimensional model file, construct a single-particle discrete element model. The proppant particle diagram, the scanned three-dimensional model, and the single-particle discrete element model are as follows: Figure 1 As shown;

[0044] (2) Conduct uniaxial compression tests on single proppant particles, record the deformation of each particle during pressure loading, and use the following formula to convert the force F into stress σ to obtain the stress-strain curve of a single particle:

[0045]

[0046] In the formula, σ is the compressive stress on the particle, MPa; F is the loading force on the particle, N; and a is the effective contact area at the bottom of the particle, m². 2 After the experiment, the particle strain under different loading pressure conditions was recorded.

[0047] (3) The single-particle discrete element model described above was used to simulate the uniaxial failure process of a single particle. A servo method was used to apply the same uniaxial pressure as in the experiment to the model. During this process, the strain of the single-particle discrete element model under different stresses was monitored and recorded, forming the stress-strain curve of the single particle under compression. The discrete element particles in the model were selected as linear elastic materials. Parameters such as the effective modulus, stiffness ratio, and friction coefficient of the particles in the model were adjusted to fit the simulation results of the uniaxial compression experiment. Specifically, as follows... Figure 2 As shown, the discrete element modeling parameters for a single particle are obtained.

[0048] III. Obtaining the elastoplastic constitutive mechanical parameters of the proppant particle packing includes the following steps:

[0049] Step C1: Establish a discrete element model of the proppant particle pack. The particle size in the model is set according to the proppant particle size distribution obtained from particle size analysis tests. The effective modulus, stiffness ratio, and friction coefficient of the particles are the same as the fitting results in Step B. Apply confining pressure to the discrete element model of the proppant particle pack using a servo method to ensure close contact between the particles in the model. Specifically:

[0050] (1) First, a group of loose particles is generated in the computational domain according to the test results of the particle size analyzer. To restore the particle packing state under the condition of ground stress compaction, the particle packing unit is established by closing the loading plate after randomly generating particles. In this embodiment, a cubic geometric space with a side length of 3.6 cm is first established, and then rigid cluster particles for simulating proppant are generated in it. Each rigid cluster is composed of several sub-particles, and the sub-particles constituting the cluster are in rigid contact with each other. The particle size range in this embodiment is 20 / 40 mesh; the single particle parameters (effective modulus, stiffness ratio, friction coefficient) in the discrete element model of proppant particle packing are selected according to the fitting results in step B.

[0051] (2) For discrete element problems with stress boundaries, after the granular particle group is generated, a servo process is needed to adjust the contact between the particle system to ideal conditions in order to achieve a force equilibrium state of the particles. In this embodiment, as shown... Figure 3 As shown, to simulate the confining pressure conditions of a real reservoir environment, stress boundary conditions are set for the discrete element model. Constraint forces are applied to the boundary by applying a small velocity. The confining pressure is loaded through a servo process to form a proppant particle accumulation.

[0052] Step C2: Apply a closed pressure to the discrete element model of the proppant particle packing to simulate the uniaxial compression process of the proppant particles under confining pressure; record the stress-strain curves during the compression process, and obtain the elastic modulus E in the elastic stage and the stress yield strength σ in the plastic stage of the proppant particles based on these curves. y Tangent modulus E t These parameters, namely the elastoplastic constitutive mechanical parameters of granular aggregates, specifically include:

[0053] (1) After the proppant particles reach equilibrium, uniaxial compression of the sample is started to obtain the stress-strain curve during the rock sample compression process. Figure 4 The table shows the stress-strain curves of the particles under different confining pressures. It can be seen that the stress-strain curves of the particles under triaxial confining pressure conform to the constitutive characteristics of typical elastoplastic materials. When the stress applied by the loading plate is lower than the material's elastic limit stress σ... y0 At this stage, the material is in the linear elastic stage, during which the strain and the applied stress satisfy Hooke's law, and the slope of the stress-strain curve is equal to the elastic modulus of the material. When the stress exceeds the elastic limit, the proppant particles undergo uniform plastic deformation, and as the amount of plastic deformation increases, the material's resistance to plastic deformation gradually increases. This stage is the strain hardening stage, until the plastic stage yield strength σ is reached. y It then manifests as an ideal plastic state. The σ of the proppant particles... y like Figure 4 As shown in the image.

[0054] (2) The overall constitutive model of the proppant particles exhibits linearly reinforced elastoplastic characteristics. Its elastoplastic constitutive mechanical behavior can be described by two parameters: elastic modulus and tangential modulus. In this embodiment, the elastic modulus E and tangential modulus E of the proppant particles were simulated. t like Figure 5 As shown in the image.

[0055] IV. Calculating the dynamic fracture width closure during artificial fracture mining includes the following steps:

[0056] Step D1: Establish a 3D model of the reservoir at engineering scale, including artificial fractures. The geometric dimensions of the fractures and reservoir in the model are set according to the field fracture monitoring results. Specifically:

[0057] (1) Establish a three-dimensional model of the reservoir containing artificial fractures. This model includes three parts: simulated wellbore, simulated fractures, and reservoir matrix, such as... Figure 6 As shown;

[0058] (2) Initial values ​​and boundary conditions are assigned to the three-dimensional model of the reservoir containing artificial fractures. The fractures are set as elastoplastic materials, and their elastoplastic constitutive parameters are determined by discrete element simulation of the constitutive characteristic behavior of proppant accumulation in step C2. The reservoir part in the model is set as an ideal elastic material, and its parameters are obtained from the physical property parameter test of the reservoir core. The model uses the sweep method to create a structured grid, and the grid in the fracture area is locally refined to improve the calculation accuracy. The fracture and matrix are set as deformable and permeable interfaces, and the matrix interface is set as a pressure interface.

[0059] Step D2: Apply the finite element method to solve the solid mass balance equation, fluid mass balance equation, Biot effective stress equation, and Kozeny-Carman permeability response equation during reservoir exploitation; the solution process of the model is as follows: Figure 7 As shown in the image.

[0060] The basic principle is as follows: After assigning initial values ​​to the model, the seepage field during reservoir production is first solved based on the fluid mass conservation equation to obtain the reservoir pore pressure distribution law; the effective stress change induced by pore pressure change is solved based on the Biot effective stress principle to solve the solid continuity equation and calculate the deformation characteristics of the reservoir solid unit; the relationship between strain and porosity and permeability is established according to the Kozeny-Carmen equation, and the reservoir porosity and permeability data are updated after each time step; after checking the model convergence, the calculation for the next time step is carried out. During the model calculation, the deformation behavior prediction of the fracture body is based on the elastoplastic mechanical parameters of the proppant particle pack. After reaching the set calculation time step, the fracture width closure distribution characteristics of the propped fractures are obtained, specifically:

[0061] (1) The finite element method was used to solve the established fractured reservoir model. Both the reservoir matrix and the artificial fractures are porous materials, and their solid phase mass balance equation is shown in Equation 2 below:

[0062]

[0063] In the formula, ρ is the density, kg·m -3 φ is porosity, dimensionless; ν is velocity, m·s -1 t represents time, s; the superscript S indicates a solid.

[0064] Based on the vector integral equation:

[0065]

[0066] Substitute equation 3 into equation 2, ignoring (1-φ)ρ. S The gradient of the term over time is calculated and divided by ρ on both sides of the equation. S Equation 2 can be rewritten in the following form:

[0067]

[0068] Equation 5 is the mass balance equation for the fluid phase, where the superscript and subscript π represent the fluid phase:

[0069]

[0070] Superimposing Equations 5 and 4 yields the system mass balance equation:

[0071]

[0072] According to Biot's effective stress principle, the expression for the effective stress of the system is as follows:

[0073]

[0074] In the formula, σ is the stress, MPa; α is the Biot coefficient, dimensionless.

[0075] The system continuity equation is as follows:

[0076]

[0077] In the formula, L is the stress differential operator, which is dimensionless; σ is the stress on the solid element, in N, as shown in equation (6-13); ρ is the density in equation (6-7), in kg·m³. -3 .

[0078] (2) The porosity stress response relationship between the reservoir and the fracture is shown below:

[0079]

[0080] In the formula, ε b φ is the volumetric strain in the solid phase of the porous medium, which is dimensionless; φ0 is the initial porosity, which is dimensionless.

[0081] The relationship between pore deformation and permeability is established using the Kozeny-Carman formula, which is as follows:

[0082]

[0083] In the formula, k is the permeability, m 2 ε represents the porosity of the proppant particle packing, dimensionless; KC is the Kozeny-Carman constant, dimensionless; S b The specific surface area of ​​the proppant particle pack (mp ... -1 .

[0084] Step D3: Using the finite element model of the reservoir with artificial fractures established in step D1 and the model solution method proposed in step D2, calculate the fracture width closure of the artificial fractures under the fluid-solid coupling effect during the mining of the fractured reservoir, and analyze the main factors affecting the fracture width closure characteristics.

[0085] Figure 8 The figure shows the variation of fracture width closure value with production time under different initial fracture widths obtained from numerical simulation in this embodiment. The flow process of fluid in artificial fractures will significantly disturb the pore pressure, and the larger the high conductivity fracture region, the higher the induced stress value, thus forming a greater degree of fracture deformation. When the well has been producing for 10 days, due to its proximity to the bottom of the well, the fluid velocity at the fracture opening is relatively high, and the fracture surface displacement gradually decreases along the fracture propagation direction; however, as reservoir production continues and the production-induced stress increases, the supporting fractures near the wellbore gradually enter the plastic deformation stage, the rate of increase of fracture width closure value at the fracture opening gradually slows down, while the deformation in the middle and end regions of the fracture, which are in the elastic stage, continues to increase. During the production time of 20 to 40 days, the deformation along the fracture surface gradually tends to level off. Figure 8 It can also be seen that the increase in crack width effectively increases the interaction area between the fluid and the skeleton in the crack, significantly increases the effective stress generated by the fluid-structure interaction effect in the crack, and triggers a greater degree of support crack width closure.

[0086] This invention combines two numerical simulation methods, discrete element method and finite element method, to obtain the mechanical behavior of proppant particle packing and further conduct finite element numerical simulation of fracture width closure under reservoir mining conditions, thereby obtaining the fracture width damage law of the elastic deformation segment and plastic deformation segment of the proppant particle body.

[0087] The above description is merely a preferred embodiment of the present invention and is not intended to limit the present invention in any other way. Any person skilled in the art may make changes or modifications to the above-disclosed technical content to create equivalent embodiments for application in other fields. However, any simple modifications, equivalent changes, and modifications made to the above embodiments based on the technical essence of the present invention without departing from the scope of the present invention shall still fall within the protection scope of the present invention.

Claims

1. A method for calculating the dynamic crack width of artificial cracks using a combination of discrete element method and finite element method, characterized in that, Includes the following steps: Step A: Obtain the particle size distribution data of the proppant used on site, perform a three-dimensional scan on a single particle, and obtain the geometry of the single particle; Step B: Construct a single-particle discrete element model and obtain the discrete element modeling parameters for a single particle; Step C: Combine the discrete element modeling parameters obtained in Step B to establish a discrete element model of the proppant particle pack and obtain the elastoplastic constitutive mechanical parameters of the proppant particle pack. Step D: Apply the elastoplastic constitutive mechanical parameters of the proppant mass to the dynamic fracture width closed finite element calculation during the artificial fracture mining process to realize the calculation of the dynamic fracture width of the artificial fracture. The finite element calculation includes: establishing a three-dimensional model of the reservoir with artificial fractures at the engineering scale. The geometric dimensions of the fractures and reservoir in the model are set according to the on-site fracture monitoring results, and its elastoplastic constitutive parameters are determined by the discrete element simulation results in step C; and applying the finite element method to solve the solid phase mass balance equation, fluid mass balance equation, Biot effective stress equation, and Kozeny Carman permeability response equation during the reservoir exploitation process. Calculate the fracture width closure of artificial fractures under the fluid-solid coupling effect during the exploitation of fractured reservoirs, and analyze the main factors affecting the fracture width closure characteristics.

2. The method for calculating the dynamic crack width of artificial cracks using the discrete element method combined with the finite element method according to claim 1, characterized in that: Step B is implemented in the following manner: Based on the particle size distribution data and geometry obtained in step A, a single-particle discrete element model is constructed. A single-particle uniaxial compression experiment is conducted, and the same pressure conditions as those in the single-particle uniaxial compression experiment are applied to the single-particle discrete element model to obtain the uniaxial compression stress-strain curve of the single particle during the loading process. The effective modulus, stiffness ratio, and friction coefficient of the discrete element linear model are adjusted, and the stress-strain curve obtained from the simulation is fitted with the uniaxial compression experimental results of a single particle to determine the constitutive characteristic parameters of the single particle, thereby obtaining the discrete element modeling parameters of the single particle.

3. The method for calculating the dynamic crack width of artificial cracks using the discrete element method combined with the finite element method according to claim 1, characterized in that: Step C specifically includes the following steps: Step C1: Establish a discrete element model of the proppant particle pack. The particle size in this model is set according to the particle size distribution in step A. The effective modulus, stiffness ratio, and friction coefficient of the particles are the same as the fitting results in step B. Step C2: Apply closed pressure to the discrete element model of proppant particle packing to simulate the uniaxial compression process of proppant particles under confining pressure; record the stress-strain curves during the compression process, and obtain the elastic modulus of the proppant particles in the elastic stage, the stress yield strength in the plastic stage, and the tangential modulus parameters based on these curves.

4. The method for calculating the dynamic crack width of artificial cracks using the discrete element method combined with the finite element method according to claim 1, characterized in that: In step D, the three-dimensional model of the reservoir with artificial fractures includes a simulated wellbore, simulated fractures, and reservoir matrix. Initial values ​​and boundary conditions are assigned to the three-dimensional model of the reservoir with artificial fractures. The simulated fractures are set as elastoplastic materials, and the reservoir matrix is ​​set as an ideal elastic material. The parameters are obtained from the physical property parameter test of the reservoir core. The simulated fractures and the reservoir matrix are set as a deformable and permeable interface, and the reservoir matrix interface is set as a pressure interface.

5. The method for calculating the dynamic crack width of artificial cracks using the discrete element method combined with the finite element method according to claim 1, characterized in that: Step D is specifically solved based on the following principles: After initializing the model, the seepage field during reservoir production is first solved based on the fluid mass conservation equation to obtain the reservoir pore pressure distribution law; the effective stress change induced by pore pressure change is solved based on the Biot effective stress principle to solve the continuity equation and calculate the deformation equation of the reservoir solid phase unit; the relationship between strain and porosity and permeability is established according to the Kozeny-Carmen equation, and the reservoir porosity and permeability data are updated after each time step; after checking the model convergence, the calculation of the next time step is carried out, and after the calculation set time step is reached, the closure of the support fracture width is checked.

Images