Fluid-solid coupling numerical simulation method, device and medium for shale gas reservoir considering proppant distribution

Abstract

The invention belongs to the technical field of reservoir numerical simulation, and discloses a shale gas reservoir fluid-solid coupling numerical simulation method, device and medium considering proppant distribution. The method comprises: combining a discrete unit model, a pore network model, a lattice Boltzmann method and an artificial neural network model to establish a proppant fracture conductivity calculation method, and obtain the law of proppant fracture conductivity changing with the increase of closure stress; using a computational fluid dynamics coupled discrete unit method to realize the two-phase bidirectional coupling of proppant flow particles and fluid in a hydraulic fracture, forming a proppant migration numerical simulation method in a hydraulic fracture, and obtaining the actual distribution characteristics of the proppant in the hydraulic fracture; combining the above-mentioned proppant migration model and the proppant fracture conductivity calculation model to establish a shale gas reservoir fluid-solid coupling integrated numerical simulation method considering the actual distribution characteristics of the proppant in the fracture and the stress sensitivity of the proppant fracture area conductivity, so as to improve the accuracy of shale gas reservoir capacity analysis.

Description

Technical Field

[0001] The present invention relates to the technical field of oil reservoir numerical simulation, and more particularly to a shale gas reservoir fluid-solid coupling numerical simulation method, device and medium taking proppant distribution into consideration. Background Art

[0002] Shale gas is a crucial component of unconventional oil and gas reserves and production increases. Due to its extremely low reservoir permeability, hydraulic fracturing is key to the effective development of shale gas reservoirs. During hydraulic fracturing, large amounts of fracturing fluid are pumped into the formation through perforations, inducing the formation and expansion of hydraulic fractures. After pumping is stopped, the fluid pressure within the wellbore and fracture network drops below that within the shale gas reservoir, causing the hydraulic fractures to close. After hydraulic fracture closure, the uneven distribution caused by proppant migration and deposition divides the fracture region into supported and unsupported areas, each with different conductivity characteristics. The supported area can maintain high conductivity for production flow. Therefore, the actual distribution of proppant within the fracture, the stress sensitivity of the conductivity of the supported fracture region, and the overall fluid-solid coupling effect all have a significant impact on the accuracy of hydraulic fracture closure and shale gas production predictions. However, existing numerical simulation methods for shale gas reservoirs have not yet comprehensively considered the above-mentioned influencing factors. It is urgent to establish an integrated numerical simulation method and technology that can comprehensively consider the actual distribution characteristics of proppants, the stress sensitivity of the proppant fracture area conductivity, and the overall fluid-solid coupling effect of shale gas reservoirs, so as to improve the accuracy of shale gas reservoir production capacity prediction. Summary of the Invention

[0003] The purpose of the present invention is to provide a shale gas reservoir fluid-solid coupling numerical simulation method, device and medium that considers proppant distribution. By combining a proppant migration model and a proppant fracture conductivity calculation model, an integrated fluid-solid coupling numerical simulation method for shale gas reservoirs is established that considers the actual distribution characteristics of proppant within the fracture and the stress sensitivity of the proppant fracture regional conductivity. This improves the accuracy of shale gas reservoir productivity analysis and can provide technical support for the integrated analysis of shale gas reservoir fracturing design and productivity evaluation.

[0004] In order to achieve the above object, the present invention adopts the following technical solutions:

[0005] According to a first aspect of the present invention, a method for fluid-solid coupling numerical simulation of shale gas reservoirs taking into account proppant distribution is provided, the method comprising:

[0006] The discrete element method was used to simulate the compaction process of the proppant layer under different closure stresses. The pore network model, lattice Boltzmann method, and artificial neural network model were combined to calculate the permeability and conductivity of the propped fractures. The relationship between the aperture and conductivity of the propped section and the closure stress was obtained.

[0007] The computational fluid dynamics coupled discrete element method is used to achieve bidirectional coupling between proppant flow particles and fluid in hydraulic fractures. Computational fluid dynamics is used to simulate the conservation of momentum and mass of the sand-carrying fluid. A discrete element model is used to simulate the flow and contact of proppant particles. Fluid forces are applied to the particle phase based on drag models to achieve numerical simulation of proppant migration in hydraulic fractures and obtain the actual distribution characteristics of the proppant in the hydraulic fractures.

[0008] For shale gas reservoirs with fractured areas formed after hydraulic fracturing, a matrix-microfracture-hydraulic fracture coupling model is used for production simulation. A single-pore medium model is used outside the fractured area, a dual-medium model is used for the bedrock and microfractures within the fractured area, and an embedded discrete fracture model is used for the hydraulic fractures. The actual distribution characteristics of the proppant and the relationship between the opening and conductivity of the propping section and the closure stress are considered within the hydraulic fracture.

[0009] According to a second aspect of the present invention, there is provided a fluid-solid coupling numerical simulation device for shale gas reservoirs taking into account proppant distribution, the device comprising:

[0010] The computational model building module is configured to use the discrete element method to simulate the compaction process of the proppant layer under different closure stresses. It combines the pore network model, the lattice Boltzmann method, and the artificial neural network model to calculate the permeability and conductivity of the propped fractures and obtain the relationship between the aperture and conductivity of the propped section and the closure stress.

[0011] The proppant distribution calculation module is configured to use computational fluid dynamics coupled discrete element method to achieve bidirectional coupling between proppant flow particles and fluid in hydraulic fractures. It uses computational fluid dynamics to simulate the conservation of momentum and mass of the sand-carrying fluid, and uses a discrete element model to simulate the flow and contact of proppant particles. It applies fluid forces to the particle phase based on drag models, etc., to achieve numerical simulation of proppant migration in hydraulic fractures and obtain the actual distribution characteristics of proppant in hydraulic fractures.

[0012] The numerical simulation module is configured to conduct production simulation for shale gas reservoirs in the fractured zone formed after hydraulic fracturing using a matrix-microfracture-hydraulic fracture coupling model. A single-pore medium model is used outside the fractured zone, a dual-medium model is used for the bedrock and microfractures within the fractured zone, and an embedded discrete fracture model is used for the hydraulic fractures. The actual distribution characteristics of the proppant and the relationship between the opening and conductivity of the propping section and the closure stress are considered within the hydraulic fracture.

[0013] According to a third aspect of the present invention, a readable storage medium is provided, wherein the readable storage medium stores one or more programs, and the one or more programs can be executed by one or more processors to implement the method described above.

[0014] The present invention has at least the following beneficial effects:

[0015] Based on the Euler particle flow model to determine the distribution information of proppant in the fracture, the present invention combines the stress-strain of the proppant layer and the stress sensitivity law of proppant fracture conductivity to successfully establish a set of fluid-solid coupling numerical simulation models suitable for shale gas reservoirs. First, a hybrid model is used to simulate the shale gas reservoir containing matrix, microcracks and hydraulic fractures, in which a single pore model is used to simulate the non-reservoir reformed area, and a multi-porosity model is used to simulate the area containing reservoir reformed. The hydraulic fractures are explicitly characterized by an embedded discrete fracture model, and the influence of proppant distribution is considered in the hydraulic fracture closure simulation. Specifically, the proppant distribution information in the fracture is calculated by the Euler particle model, and the fluid-solid coupling parameters are provided for the hydraulic fracture closure process based on the stress-sensitive proppant fracture conductivity change and the stress-strain relationship of the proppant layer. Based on the simulation method of medium matrix-microcracks-hydraulic fracture coupling, the basic characteristics of shale gas flow are taken into account to analyze the factors affecting the cumulative gas production of shale gas reservoirs. This modeling approach combines proppant distribution simulation with stress-sensitive conductivity variations to construct a fluid-solid coupling model for shale gas reservoirs that takes proppant distribution into account. This improves the accuracy of cumulative gas production simulations for shale gas reservoirs and is a key technology for developing integrated numerical simulation methods. This model establishes an integrated numerical simulation method that organically integrates proppant flow analysis, proppant fracture conductivity calculation, and shale gas reservoir productivity analysis, providing a reference simulation method and theoretical foundation for the efficient development of shale gas reservoirs and the optimization of reservoir fracturing and stimulation technologies. BRIEF DESCRIPTION OF THE DRAWINGS

[0016] Figure 1 A flow chart of a fluid-solid coupling numerical simulation method for shale gas reservoirs considering proppant distribution according to an embodiment of the present invention is shown;

[0017] Figure 2 A schematic diagram of fracture closure and proppant compaction simulation according to an embodiment of the present invention is shown;

[0018] Figure 3 A schematic diagram of the process of extracting the pore structure of the proppant layer after compaction according to an embodiment of the present invention is shown;

[0019] Figure 4 A schematic diagram showing a comparison of pore network models according to an embodiment of the present invention is shown;

[0020] Figure 5 A schematic diagram of sub-throat division according to an embodiment of the present invention is shown;

[0021] Figure 6 A schematic diagram of the distribution of conductivity within a single fracture according to an embodiment of the present invention is shown;

[0022] Figure 7A schematic diagram of fracture network geometry with different numbers of fractures according to an embodiment of the present invention is shown;

[0023] Figure 8 A schematic diagram showing the change of cumulative gas production with proppant pumping volume according to an embodiment of the present invention is shown;

[0024] Figure 9 A schematic diagram of a fracture network model according to an embodiment of the present invention is shown;

[0025] Figure 10 A schematic diagram showing the effect of branch seam positions on the cumulative gas production difference according to an embodiment of the present invention is shown;

[0026] Figure 11 A schematic diagram showing the influence of branch seam positions on stress according to an embodiment of the present invention is shown;

[0027] Figure 12 A structural diagram of a shale gas reservoir fluid-solid coupling numerical simulation device considering proppant distribution according to an embodiment of the present invention is shown. DETAILED DESCRIPTION

[0028] In order to enable those skilled in the art to better understand the technical solution of the present invention, the present invention is described in detail below with reference to the accompanying drawings and specific embodiments. The embodiments of the present invention are further described in detail below with reference to the accompanying drawings and specific embodiments, but are not intended to limit the present invention. For the various steps described herein, if there is no necessity for a contextual relationship between each other, the order in which they are described as examples herein should not be regarded as limiting, and those skilled in the art should know that they can be adjusted in order as long as the logic between them is not destroyed, resulting in the inability to implement the entire process.

[0029] The embodiment of the present invention provides a fluid-solid coupling numerical simulation method for shale gas reservoirs taking into account proppant distribution. Figure 1 The flowchart of the method is shown in FIG. The method comprises the following steps:

[0030] S1. The discrete element method is used to simulate the compaction process of the proppant layer under different closure stresses. Combined with the pore network model, the lattice Boltzmann method, and the artificial neural network model, a calculation method for the permeability and conductivity of the propped fracture is established that takes into account both computational efficiency and pore throat characteristics. The relationship between the aperture and conductivity of the propped fracture and the closure stress is obtained.

[0031] In this embodiment, when establishing a calculation model for propped fracture conductivity, the calculation efficiency of the porous media pore network is improved by improving the conductivity calculation and artificial neural network (ANN) model in the pore network model. Specifically, a discrete element method is used to simulate the compaction process of the proppant layer under different closure stresses. The pore network model and the lattice Boltzmann method are combined to improve the calculation method for propped fracture permeability and conductivity. The relationship between the opening and conductivity of the support section and the closure stress is obtained, and the stress-strain and conductivity of the propped fracture are determined.

[0032] In some embodiments, step S1 comprises the following steps:

[0033] Step S11: establishing an improved support joint conductivity calculation model.

[0034] The propped fracture area is the main gas flow channel. Accurately calculating the changes in the proppant layer's conductivity under closure can characterize the fracture material properties required for fluid-solid coupling simulation, which helps improve the accuracy of shale gas reservoir productivity analysis.

[0035] The discrete element method (DEM) was used to simulate the compaction process of the proppant layer under closure stress and calculate the change in proppant fracture aperture. Fracture closure was achieved by placing a controllable movable plate above the proppant layer and applying a constant pressure. The fracture aperture was defined as the distance between the lower fixed plate and the upper movable plate. Once the movable plate reached a stable state, the pore structure and porosity within the proppant accumulation were determined by combining discrete element calculation data, a self-developed Python image processing program, and digital core technology. The pore network was then extracted using ImageJ software and OpenPNM, and the permeability and fracture conductivity of the proppant layer were calculated based on an improved pore network model.

[0036] In this embodiment, the permeability calculation of the proppant layer is an important parameter for evaluating the conductivity and gas production performance of the hydraulic fracture. For the permeability calculation of the proppant layer, the LBM combined with the PNM method is adopted.

[0037] Step S12: Improve conductivity calculation in the pore network model.

[0038] The conventional pore network model (CPNM) retains the position and connection relationship between pores and throats, and simplifies the complex pore and throat structures in real porous media into connections of regular geometric bodies, which improves the calculation efficiency but reduces the accuracy of flow calculation. To further improve the accuracy of the calculation, the regular geometric cross-sectional shape is replaced by the real cross-sectional shape of the throat to improve the accuracy of the conductivity calculation. However, for complex real cross-sectional shapes, their conductivity cannot be obtained through empirical or theoretical formulas. Therefore, it is necessary to use direct simulation methods such as LBM to simulate the single-phase flow in the throat to determine the conductivity in the improved pore network model ( g ij ) numerical values. Another improved approach is to preserve throat structural information as much as possible, performing direct numerical flow simulation based on the actual throat structure to calculate its conductivity. However, this method requires direct numerical simulation of each throat, which is computationally intensive. Therefore, developing an accurate and efficient method for calculating throat conductivity is of great significance.

[0039] The throat of a porous medium is generally smaller in cross-section than the pores, and the flow resistance primarily originates in the throat. The pore network model employed in this paper only considers the flow resistance in the throat, while the flow resistance in the pores can be neglected. Therefore, the conductivity of the throat is a key parameter in determining flow calculations in the pore network model.

[0040] The difference between the improved pore network model and the traditional pore network model is that it uses the real cross-sectional shape of the throat to replace the regular geometric cross-sectional shape to improve the accuracy of the conductivity calculation. However, for complex real cross-sectional shapes, their conductivity cannot be obtained through experience or theoretical formulas. Therefore, the present invention uses the lattice Boltzmann method to simulate the single-phase flow in the throat to determine the conductivity value in the improved pore network model. On the other hand, in order to consider the impact of the continuous change of the throat cross-sectional shape along the pore centroid line on the conductivity calculation under actual conditions, the present invention further improves the pore network model by connecting a series of simplified throats in series. The pore network model adopted by the present invention only considers the flow resistance of the throat, and the flow resistance in the pores can be ignored.

[0041] In this embodiment, the conductivity calculation formula in the improved pore network model is:

[0042] (1)

[0043] (2)

[0044] Step S13: using an artificial neural network model (ANN) to improve the calculation efficiency of the actual model permeability and fracture conductivity.

[0045] LBM flow simulations within a single throat are highly efficient for determining its conductivity characteristics. However, due to the complex pore network structure and the large number of throats in porous media, LBM flow simulations for all throats still require significant computational resources. Therefore, ANNs can be used to summarize the relationship between throat flow characteristics and cross-sectional shape, thereby avoiding traversal calculations and improving computational efficiency.

[0046] There are three types of layers in an artificial neural network model: input layer, hidden layer, and output layer. The performance of the neural network model during training is evaluated by the mean square error (MSE), mean relative error (MARE), and coefficient of determination (DE). The final output can be calculated using formula (3):

[0047] (3)

[0048] (4)

[0049] (5)

[0050] (6)

[0051] Where MSE, MARE and DE represent mean square error, mean relative error and coefficient of determination, respectively; is an input parameter; is the output value calculated by the model; is the hyperbolic tangent function; is a linear transfer function; is the total number of evaluation parameters; i The index of the parameter to be evaluated for the input layer; j is the node index of the hidden layer; k is the node index of the output layer; and Represent the target value of the parameter and the value obtained by calculation and prediction respectively; represents the average of the target values; m is the number of nodes in the hidden layer; n is the number of nodes in the input layer; represents the bias term; is the bias term of the hidden layer; is the bias term of the output layer; is the weight coefficient; is the weight between the hidden layer and the output layer, indicating the j hidden layer nodes and the k The connection weights between the output nodes, the second weight matrix; is the weight between the input layer and the hidden layer, indicating the i The input node and jThe connection weights between hidden layer nodes, the first weight matrix.

[0052] S2. Computational fluid dynamics coupled discrete element method (CFD-DEM) is used to achieve bidirectional coupling between proppant flow particles and fluid phases within hydraulic fractures. Computational fluid dynamics is used to simulate the momentum and mass conservation of the sand-carrying fluid. A discrete element model is used to simulate the flow and contact of proppant particles. Fluid forces are applied to the particle phase based on drag models, forming a numerical simulation method for proppant migration within hydraulic fractures and obtaining the actual distribution characteristics of proppant in hydraulic fractures.

[0053] In this example, Eulerian particle simulation was used to determine proppant distribution within fractures. The fluid phase was considered the primary phase, while the proppant was treated as the secondary phase. Both phases shared a unified pressure field during the numerical calculations. By analogizing the random particle collision motion with the thermal motion of gas molecules, the proppant phase stress was derived. The concept of particle phase temperature was introduced to establish a proppant migration model, which was embedded in the proppant distribution calculations. This model determined information such as the particle-phase stress and ultimately revealed the proppant distribution in the fractures.

[0054] In some embodiments, step S2 specifically includes the following steps:

[0055] Step S21: achieving bidirectional coupling between proppant flow particles and fluid in the hydraulic fracture.

[0056] The commonly used method for achieving bidirectional coupling between proppant flow particles and the fluid phase within hydraulic fractures is the computational fluid dynamics coupled discrete element method (CFD-DEM). This requires consideration of the transport of the fluid phase relative to the proppant phase (primarily due to drag forces). Furthermore, the influence of accumulated proppant on fluid flow within narrow fractures cannot be ignored. The overall coupled computational process involves a computational fluid dynamics solver, flow field information update, particle-fluid interaction calculation, particle volume fraction calculation on the fluid grid, particle information update, and a discrete element method solver. Fluid flow information, such as velocity and pressure distribution, is obtained by solving the Navier-Stokes (NS) equations using the finite volume method, and the calculation is performed using the ANSYS Fluent solver. Contact and motion within the proppant particle population are simulated using the discrete element method, implemented using the EDEM solver. Furthermore, the forces acting on the proppant particles and the influence of accumulated proppant on fluid flow must be considered to achieve bidirectional coupled calculations. Information transfer between the computational fluid dynamics solver and the discrete element solver is a key component of the phase coupling and is implemented using an ANSYS Fluent user-defined interface program.

[0057] Specifically, at each fluid computation time step, the solid phase volume fraction is calculated within each fluid grid. This describes the effect of the cumulative change in proppant particles within the fracture region on the fluid mass conservation equation. A source term is introduced into the momentum conservation equation to reflect the proppant's reaction force on the sand-carrying fluid. The Navier-Stokes equations are then solved to obtain converged results for the fluid velocity, volume fraction, and pressure distribution within the study area. This fluid phase information is then shared with the discrete element solver, initiating proppant phase computation. The fluid forces acting on the proppant particles are primarily drag forces, which are calculated using empirical formulas, assuming that all particles within the same fluid grid experience the same fluid force. In addition to the fluid forces, particle-particle and particle-wall contact forces are calculated using the discrete element method. Based on the net force, the velocity and position of each proppant particle are updated according to Newton's second law. Next, the proppant particle velocity, diameter, and position information are passed to the fluid phase solver to calculate the particle phase volume fraction within the fluid grid and the source term in the fluid momentum conservation equation, thus accounting for the proppant's influence on the fluid phase. The process of bidirectional coupling of proppant and fluid is completed by cycling in sequence.

[0058] This coupling method requires that the fluid computational grid resolution be larger than the particle size to ensure that multiple particles are contained within each fluid grid. Furthermore, the discrete element solver time step is typically much smaller than the fluid phase computation time step. Within each fluid computation time step, the discrete element is computed an integer multiple of the times to ensure that the computational times of the two solvers match.

[0059] Step S22: using the Euler particle model to determine the proppant distribution information within the fracture.

[0060] The Euler particle model describes the flow characteristics of a proppant-fluid mixture based on a solid-phase pressure model and kinetic theory. The fluid phase is considered the primary phase, while the proppant is considered the secondary phase. Both phases share a unified pressure field in numerical calculations. The subscripts p and q refer to different phases, respectively. The mass and momentum conservation equations for the q phase are:

[0061] (7)

[0062] (8)

[0063] Where, represent q Phase volume fraction; represent q Phase velocity vector; represent p Phase velocity vector; represent q density of the phase; t Represents time; represents the gradient operator; 、 Represent phase q Relative to the phase p The relative speed, phase p Relative to the phase q relative speed; n represents the number of phases; represents the mass transfer from the p phase to the q phase; represents the mass transfer from the q phase to the p phase; is the acceleration due to gravity; K pq for p and q interphase momentum exchange coefficient; represent q External body forces of the phase; represent q Phase lift; represent q The wall lubrication force of the phase; express q Virtual mass force of phase; express q Phase drag.

[0064] is the q-corresponding force-strain tensor, which can be calculated using formula (9):

[0065] (9)

[0066] Where, represents the volume fraction; and Respectively represent q Shear viscosity and bulk viscosity of the phase; represents the unit tensor; represents the transpose of the velocity matrix;

[0067] Represents the phase interface velocity, in m·s -1 , which is determined according to the direction of mass transfer between phases:

[0068] (10)

[0069] Step S23: Determine the specific information of the particle phase.

[0070] In the Euler particle model, the stress of the proppant phase is derived by analogizing the random particle collision motion to the thermal motion of gas molecules. Introducing the concept of particle phase temperature, the particle motion fluctuation is linearly related to the temperature value, and the momentum conservation equation of the particle phase s can be written as formula (11)

[0071] (11)

[0072] Where, represents the pressure of the particle phase; t Indicates time; represents the momentum exchange coefficient between the fluid phase and the particle phase; N Indicates the total number of phases contained in the system; represents the volume fraction of the particle phase; represents the particle phase pressure; represents the particle phase velocity vector; represents the gradient operator; n represents the number of phases; express p Phase velocity vector; Represents the fluid phase l Relative to the particle phase s relative speed; Represents the particle phase s Relative to the fluid phase l relative speed; Represents the mass transfer between the fluid phase and the particle phase; Represents the mass transfer between the particle phase and the fluid phase; represents the external volume force of the particle phase; represents the lift force of the particle phase; represents the wall lubrication force of the particle phase; represents the virtual mass force of the particle phase; represents the drag force of the particle phase; is the acceleration due to gravity, represents the stress-strain tensor of the particle phase;

[0073] When the volume fraction is less than the critical value of accumulation, the particle phase pressure can be calculated by the Syamlal-O'Brien formula:

[0074] (12)

[0075] Where, represents the particle collision restitution coefficient; represents the radial distribution function; represents the particle temperature, K.

[0076] Stress-strain tensor of the particle phase It can be written as:

[0077] (13)

[0078] (14)

[0079] Where, and represent the shear viscosity and bulk viscosity of the particle phase, respectively; represents the unit tensor, represents the particle diameter; e s represents the elastic recovery coefficient of the particle phase;

[0080] When the volume fraction of the particle phase reaches the critical value of stacking, the shear viscosity (including collision part, dynamic part and friction part) can be expressed as:

[0081] (15)

[0082] Where, represents the particle diameter; represents the internal friction angle; parameter n and m is an empirical constant; Indicates the minimum value of particle phase concentration; Indicates the maximum value of the particle phase concentration.

[0083] S3. For shale gas reservoirs with fractured zones formed after hydraulic fracturing, a matrix-microfracture-hydraulic fracture coupling model is used for production simulation. A single-pore medium model is used outside the fractured zone, a dual-medium model is used for the bedrock and microfractures within the fractured zone, and an embedded discrete fracture model is used for the hydraulic fractures. Furthermore, the actual proppant distribution characteristics calculated in step 2 and the relationship between proppant section aperture and conductivity as a function of closure stress calculated in step 1 are considered within the hydraulic fractures. Based on this, an integrated fluid-solid coupling numerical simulation method for shale gas reservoirs is established, which takes into account the actual proppant distribution characteristics within the fractures and the stress sensitivity of the conductivity in the proppant fracture region.

[0084] In this embodiment, hydraulic fractures create a stimulated fracture region (SRV) around the horizontal well in a hydraulically fractured shale gas reservoir. A single-pore medium model is used outside the stimulated fracture region, while a dual-medium model is used for the bedrock and microfractures within the stimulated fracture region. An embedded discrete fracture model is used for the hydraulic fractures. This bedrock-microfracture-hydraulic fracture coupling model is established, resulting in an improved fluid-solid coupling model for shale gas reservoirs that takes into account proppant distribution. Numerical simulations are then performed based on this improved fluid-solid coupling model. For example, numerical simulations can be used to analyze closed fracture conductivity, cumulative gas production, reservoir pressure, and stress distribution patterns under different scenarios. This approach is applicable to simulations involving both single-group fractures and multi-group fractures via staged fracturing.

[0085] In some embodiments, step S3 specifically includes the following steps:

[0086] Step S31: Establishing a shale gas seepage control equation.

[0087] Based on the characteristics of shale gas reservoirs, the following assumptions are made for the gas seepage model: ① Not only free gas but also adsorbed gas exists in the shale matrix; ② adsorbed gas is not considered in microfracture and hydraulic fracture systems; ③ free gas is real gas, and shale reservoirs are slightly compressible; ④ gas adsorption obeys the Langmuir isothermal adsorption model; ⑤ gas migration in the shale matrix considers the Kinkenberg effect, and gas migration in microfracture and hydraulic fracture systems obeys Darcy's law; ⑥ shale gas adsorption, desorption, and flow are all isothermal processes.

[0088] Gas flow satisfies the law of conservation of mass, and its mathematical expression is as follows:

[0089] (16)

[0090] Where, Representation Boundary The unit normal vector of ; 、 and They are regions The mass accumulation term, mass flow term and source-sink term in . Among them, the mass accumulation term is expressed as:

[0091] (17)

[0092] Where, is the Lagrangian porosity considering the fluid-solid coupling effect; is the gas density in kg·m -3 ; represents the amount of gas adsorbed, which is a function of gas pressure and can be determined by Langmuir adsorption isotherm:

[0093] (18)

[0094] (19)

[0095] Where, is the gas pressure, in Pa; is the molar mass of the gas, in g·mol -1 ; is the gas compressibility coefficient; is the universal gas constant, which is 8.314 kJ·(kmol·K) -1 ; is the reservoir temperature in K; is the rock bulk density in kg·m -3 ; is the gas density under standard conditions, in kg·m -3 ; Langmui pressure, which represents the gas pressure when the maximum adsorption capacity is half, and the unit is Pa; is the Langmuir volume, which indicates the maximum adsorption capacity of the shale matrix, and the unit is m 3 kg -1 .

[0096] Furthermore, the mass flow term in formula (16) It can be written as:

[0097] (20)

[0098] (twenty one)

[0099] (twenty two)

[0100] (twenty three)

[0101] Where, is the absolute permeability, in m 2 ; is the gas viscosity, in Pa·s; is the acceleration due to gravity, in m·s -2 ; is the reservoir depth in m; is the Klinkenberg coefficient, which introduces the gas slip effect into the gas flow model by expressing the apparent permeability as a function of gas pressure; is the compressibility coefficient of the gas, in Pa -1 ; is the real gas Knudsen diffusion coefficient.

[0102] Step S32: Establishing the shale reservoir deformation control equation.

[0103] During the mining process, the reservoir undergoes a small static linear elastic deformation. The tensile stress is positive and the compressive stress is negative. The reservoir stress balance equation is:

[0104] (twenty four)

[0105] Where, is the total stress tensor; is the rock density, in kg·m -3 The effective stress equations of the unfractured area and the fractured area are:

[0106] (25)

[0107] Where, represents the shale matrix elastic tensor, Pa; is the Biot coefficient; is the unit tensor; and represent the equivalent elastic tensor and bulk modulus of the fracture unit in the fracturing stimulation area; b l is the coupling coefficient; l represents the first l subgrid blocks; Indicates the l The total pressure in each subunit; SRV represents the fracture stimulation area;

[0108] According to the small deformation assumption, the strain tensor ε With displacement vector u The following relations are satisfied:

[0109] (26)

[0110] Where, represents the gradient operator, the superscript T Represents matrix transpose.

[0111] Step S33: Establishing constitutive relations.

[0112] The geometric shape of the hydraulic fracture consists of two parallel plates. The degree of closure of the proppant section of the hydraulic fracture mainly depends on the deformation degree of the proppant layer. Therefore, its constitutive model can be written as formulas (27) to (29):

[0113] (27)

[0114] (28)

[0115] (29)

[0116] Where, The nonlinear stress-strain relationship that represents the deformation of the proppant layer can usually be obtained from the proppant compaction mechanical experiment. is the proppant normal strain; It represents the relative displacement difference between the two fracture surfaces of the hydraulic fracture, in m; is the initial crack opening, in m; represents the crack normal vector; represents the displacement above the crack surface; represents the displacement below the crack surface.

[0117] The closed constitutive model of the unsupported section of the hydraulic fracture satisfies the normal contact constraint of the free surface. The specific expression is as follows:

[0118] (30)

[0119] (31)

[0120] Where, is the residual hydraulic fracture aperture, in m; is the normal contact penalty parameter, with a value of 1×10 20 , this constitutive model allows a small degree of interpenetration between contacting crack surfaces.

[0121] Based on the Kozeny-Carman model, the matrix dynamic permeability is:

[0122] (32)

[0123] The dynamic permeability of natural fractures is:

[0124] (33)

[0125] Where, Represents the matrix permeability, unit is m 2 ; is the porosity; subscript 0 represents the initial state parameter; is the bulk modulus of the natural fracture, in Pa; is the volumetric strain of the mesh block.

[0126] According to the parallel plate theory and the cubic law, the dynamic permeability of an unsupported fracture can be written as:

[0127] (34)

[0128] A large number of experimental results show that the permeability of the support joint is related to the effective stress, and its generalized expression is as follows:

[0129] (35)

[0130] Where, This equation represents the relationship between fracture permeability and effective stress. It was obtained through numerical experiments on stress sensitivity analysis of fracture conductivity.

[0131] Step S34: solving the fluid-solid coupling model of the shale gas reservoir.

[0132] In this embodiment, the shale gas reservoir fluid-solid coupling model includes the governing equations and models in steps S31 to S33. The specific steps of solving the model using a hybrid numerical discretization method combining finite volume and discontinuous finite element methods include: ① First, the reservoir and hydraulic fractures are each meshed separately. Based on the coordinate information of the reservoir and fracture grids, the fracture grid system is embedded into the reservoir grid system; ② The integral form of the seepage field equation is established, numerically discretized using the finite volume method, the time term is discretized using implicit difference, and the residual form of the seepage field equation is implicitly solved using the Newton iteration method; ③ The incremental form of the stress field equation is established, numerically discretized using the discontinuous finite element method, and implicitly solved using the Newton iteration method; ④ The seepage field equation and the stress field equation are sequentially solved using a fixed stress splitting algorithm, and the overall fluid-solid coupling model is iterated repeatedly until convergence. During the iteration process, physical parameters such as reservoir porosity, permeability, fracture aperture, and fluid properties are updated in real time.

[0133] Thus, the embodiments of the present invention have established a bedrock-microfracture-hydraulic fracture coupling model for shale gas reservoirs, forming an improved fluid-solid coupling model for shale gas reservoirs that can take proppant distribution into consideration.

[0134] In order to demonstrate the feasibility and superiority of the present invention, the following examples are given. Example 1:

[0135] Set the numerical simulation parameters for shale gas reservoir production according to the following table:

[0136] Table 1. Parameter settings for numerical simulation of shale gas reservoir production

[0137] ,

[0138] First, to reduce the amount of calculation, only a portion of the horizontal well staged fracturing area was intercepted for simulation. The calculated seepage area size was 250m×120m×30m, which included a horizontal well and a hydraulic fracture, and the perforation position was located at the center of the fracture. x Direction and y The horizontal stresses in the direction are 40MPa and 35MPa respectively. z The overburden stress was 30 MPa, and the flow boundary was a closed boundary. The detailed settings for the macroscopic numerical simulation parameters are summarized in Table 1. To fully account for the impact of geostress on gas flow in the reservoir, the calculation area was expanded during the stress field simulation. The expanded stress field calculation area was 1250 m × 1120 m × 230 m.

[0139] Then, using the method described in the present invention, during the simulation process, 42 MPa high-pressure fluid was first injected into the hydraulic fracture to simulate the hydraulic fracturing process to obtain the initial fracture aperture distribution. Then, the proppant filling situation in the hydraulic fracture was initialized based on the proppant distribution data obtained from the proppant flow simulation.

[0140] Table 2 is the parameter setting table for proppant flow numerical simulation

[0141] ,

[0142] The distribution of proppant within a single fracture can be summarized into three representative profiles: cumulative distribution near the perforation location (Type 1), sand bank distribution in the middle of the fracture (Type 2), and a narrow, long, and thin sand-spreading pattern (Type 3). Using the increase in the viscosity of the sand-carrying fluid as an example, proppant profiles were simulated using the parameters in Table 2.

[0143] Figure 6 The simulation results of the conductivity after fracture closure show that the propped fracture region has significantly higher conductivity than the unsupported region, and a highly conductive channel exists at the interface between the two regions. This is because the two regions have different closing behaviors under stress. The propped region can maintain a relatively large opening, while the crack opening at the interface changes from large to small, forming a "bridging" phenomenon, resulting in higher conductivity. This is something that some traditional models cannot simulate, demonstrating the model's ability to accurately simulate fluid-structure interactions such as fracture closure and proppant support. Example 2:

[0144] 1m 3 min -1 , 3m 3 min -1 , 6 m 3 min -1 , 9 m 3 min -1 and 12 m 3 min -1 The pumping rate is 0.1, and a certain amount (15 m 3 ) Proppant enters the staged area. Figure 7 Showing 3 m 3 min -1 and 12 m 3 min -1 The fracture conductivity distribution corresponding to the rate and the proppant distribution ratio of the three groups of fractures. Among them, the fracture conductivity data is adopted The simulation in this example focuses on the influence of limited pumping volume, emphasizing the importance of reasonable distribution of proppant among multiple groups of fractures and filling of the near-wellbore area. Example 3:

[0145] Figure 8 The maximum cumulative gas production in the middle segment area increases with the number of fracture groups. However, it can be seen that when the number of groups is greater than 3 and the pumping volume is greater than 25m 3 After that, increasing the fracture density and pumping volume has little effect on improving the cumulative gas production. Therefore, based on the numerical simulation method of the present invention, reference data that is helpful for optimizing the fracturing design can be obtained. Example 4:

[0146] Figure 9 This figure shows a fracture model with three sets of hydraulic fractures and branch fractures. Cumulative gas production is compared when the branch fracture distance (X) is 25 m, 50 m, and 70 m. Branch fractures divert proppant migration within the hydraulic fractures, reducing the length or height of the proppant bank within the main fracture. Figure 10 When the middle branch fracture is located far away from the branch fracture and the horizontal well (X = 75 m), its diversion effect is smaller, and it is also beneficial to increase the conductivity of the deep area, and the cumulative gas production gradually increases. Figure 11 Shows the x-direction principal stress perpendicular to the hydraulic fracture wall under different conditions The smaller the principal stress, the smaller the normal stiffness of the fracture region, corresponding to the unsupported area of ​​the hydraulic fracture. It can be seen that the unsupported area of ​​the hydraulic fracture decreases with the increase of the distance between the branch fractures.

[0147] The embodiment of the present invention also provides a shale gas reservoir fluid-solid coupling numerical simulation device taking into account the distribution of proppant, such as Figure 12 As shown, the device includes:

[0148] The calculation model building module 121 is configured to establish a calculation model for the conductivity of the support joint and determine the stress-strain and conductivity of the support joint;

[0149] The proppant distribution calculation module 122 is configured to establish a proppant migration model to perform embedded proppant distribution calculation, determine particle response force information, and obtain fracture proppant distribution;

[0150] The numerical simulation module 123 is configured to perform numerical simulation using the improved fluid-solid coupling model of the shale gas reservoir; wherein, a single-pore medium model is used outside the fracturing zone, a dual-medium model is used for the bedrock and micro-fractures within the fracturing zone, and an embedded discrete fracture model is used for hydraulic fractures.

[0151] In some embodiments, the computational model building module is further configured to:

[0152] Establishing an improved proppant fracture conductivity calculation model, wherein the improved proppant fracture conductivity calculation model uses a discrete element method to simulate the compaction process of the proppant layer under the action of closure stress and calculates the change in proppant fracture aperture; wherein fracture closure is achieved by placing a movable plate above the proppant layer and applying a constant pressure, and the proppant fracture aperture is the distance between the lower fixed plate and the upper movable plate. When the movable plate reaches a stable state, the pore structure and porosity within the proppant accumulation are determined, the pore network is extracted, and the permeability and fracture conductivity of the proppant layer are calculated based on the improved pore network model;

[0153] The conductivity calculation formula in the improved pore network model is:

[0154] (1)

[0155] (2)

[0156] Where, is the conductivity; k is the marker symbol of the sub-throat, indicating different segments in the pore network; For the throat k Conductivity of the segment; is the length of the subthroat; is the fluid dynamic viscosity; is the displacement acceleration; is the volume flow rate; is the fluid viscosity in the LBM method; is the fluid density;

[0157] The artificial neural network model is used to improve the computational efficiency of the improved pore network model.

[0158] It should be noted that the various device structures described in this embodiment belong to the same technical concept as the method described previously, and achieve the same technical effects through the same principles, which will not be repeated here.

[0159] An embodiment of the present invention further provides a readable storage medium, which stores one or more programs. The one or more programs can be executed by one or more processors to implement the method described in any of the above embodiments.

[0160] Furthermore, although exemplary embodiments have been described herein, the scope includes any and all embodiments based on the present invention with equivalent elements, modifications, omissions, combinations (e.g., cross-cutting schemes of various embodiments), adaptations, or changes. The elements of the claims are to be interpreted broadly based on the language employed in the claims and are not limited to the examples described in this specification or during the prosecution of this application, which examples are to be construed as non-exclusive. Therefore, it is intended that the specification and examples be considered as examples only, with the true scope and spirit being indicated by the following claims and their full scope of equivalents.

[0161] The above description is intended to be illustrative and not restrictive. For example, the above examples (or one or more versions thereof) may be used in combination with each other. For example, a person of ordinary skill in the art may use other embodiments when reading the above description. In addition, in the above detailed description, various features may be grouped together to simplify the invention. This should not be interpreted as an intention that a feature of an unclaimed invention is essential to any claim. On the contrary, the subject matter of the present invention may have less than all the features of an embodiment of a particular invention. Therefore, the following claims are incorporated into the detailed description as examples or embodiments, with each claim independently serving as a separate embodiment, and it is contemplated that these embodiments may be combined with each other in various combinations or permutations. The scope of the invention should be determined with reference to the appended claims and the full scope of equivalents to which such claims are entitled.

Claims

1. A fluid-solid coupling numerical simulation method for shale gas reservoirs considering proppant distribution, characterized in that: The method comprises: The compaction process of the proppant layer under different closure stresses was simulated by using the discrete element method. The permeability and conductivity of the propping seams were calculated by combining the pore network model, the lattice Boltzmann method and the artificial neural network model. The relationship between the opening and conductivity of the propping section and the closure stress was obtained. The computational fluid dynamics coupled discrete unit method is used to achieve the two-way coupling of proppant flow particles and fluid in hydraulic fractures. The computational fluid dynamics is used to simulate the momentum and mass conservation of the sand-carrying fluid. The discrete unit model is used to simulate the flow and contact of proppant particles. Fluid forces are applied to the particle phase according to the drag model to achieve numerical simulation of proppant migration in hydraulic fractures and obtain the actual distribution characteristics of proppant in hydraulic fractures. For shale gas reservoirs with fractured areas after hydraulic fracturing, a matrix-microfracture-hydraulic fracture coupling model is used for production simulation, in which a single-pore medium model is used outside the fractured area, a dual-medium model is used for bedrock and microfractures in the fractured area, and an embedded discrete fracture model is used for hydraulic fractures; the actual distribution characteristics of proppants and the relationship between the opening and conductivity of the propping section and the closing stress are considered in the hydraulic fractures; The discrete element method is used to simulate the compaction process of the proppant layer under different closure stresses. The pore network model, lattice Boltzmann method and artificial neural network model are combined to calculate the permeability and conductivity of the propping seam, and the relationship between the opening and conductivity of the propping section and the closure stress is obtained, including: The discrete element method is used to simulate the compaction process of the proppant layer under different closing stresses and calculate the change in the proppant crack opening. The crack closure is achieved by setting a controllable movable plate above the proppant layer and applying a constant pressure. When the movable plate reaches a stable state, the pore structure and porosity inside the proppant accumulation are obtained by combining discrete element calculation data, image processing and digital core technology, and the pore network is extracted. The permeability and fracture conductivity of the proppant layer are calculated based on the improved pore network model; the conductivity calculation formula in the improved pore network model is: In the formula, g μ is the conductivity; k is the marker symbol of the subthroat; it represents different segments in the pore network; g k,μ is the conductivity of the subthroat k segment; L k is the length of the throat; μ is the fluid dynamic viscosity; a x is the displacement acceleration; q k is the volume flow rate; μ LBM is the fluid viscosity in the LBM method; ρ is the fluid density; The artificial neural network model is used to improve the computational efficiency of the improved pore network model.

2. The method according to claim 1, characterized in that The computational fluid dynamics coupled discrete element method is used to achieve the two-way coupling of proppant flow particles and fluid in hydraulic fractures. The computational fluid dynamics is used to simulate the momentum and mass conservation of the sand-carrying fluid. The discrete element model is used to simulate the flow and contact of proppant particles. The fluid force is applied to the particle phase according to the drag model to achieve the numerical simulation of proppant migration in hydraulic fractures and obtain the actual distribution characteristics of proppant in hydraulic fractures, including: In each fluid calculation time step, the solid phase volume fraction in each fluid grid is calculated to describe the influence of the cumulative change of proppant particles in the fracture area on the fluid phase mass conservation equation, and the force source term is introduced into the momentum conservation equation to obtain the NS equation. The NS equation is solved to obtain the converged fluid velocity, volume fraction and pressure distribution in the region, and the fluid phase information is shared with the discrete unit solver. The calculation of the proppant phase is started. The fluid action on the proppant particles is calculated by an empirical formula, and all particles in the same fluid grid are subjected to the same fluid force. In addition to the fluid force, the particle-particle and particle-wall contact forces are calculated based on the discrete unit method, and the velocity and position information of each proppant particle is updated according to the particle-particle and particle-wall contact forces. The velocity, diameter and position information of the proppant particles are transmitted to the fluid phase solver to calculate the particle phase volume fraction in the fluid grid and the source term of the force in the fluid momentum conservation equation, determine the influence of the proppant on the fluid phase, and repeat in sequence to complete the proppant-fluid bidirectional coupling simulation process; The Euler particle model is used to determine the proppant distribution information within the fracture; The particle stress information is determined based on the Euler particle model.

3. The method according to claim 1, characterized in that The artificial neural network model is used to improve the computational efficiency of the improved pore network model, including: The artificial neural network model was evaluated using mean square error, mean relative error and coefficient of determination; The final output of the artificial neural network model is calculated, and the mean square error, mean relative error and determination coefficient are used to evaluate the performance of the artificial neural network model during the training process: wherein the calculation formula of the final output of the artificial neural network model and the calculation formula of the mean square error, mean relative error and determination coefficient are as follows: Where MSE, MARE and DE represent mean square error, mean relative error and determination coefficient respectively; x i is the input parameter; y is the output value calculated by the model; F is the hyperbolic tangent function; l is the linear transfer function; N is the total number of evaluation parameters; i is the serial number of the input layer evaluation parameter; j is the node index of the hidden layer; k is the node index of the output layer; y i and Represent the target value of the parameter and the value obtained by calculation and prediction respectively; represents the average value of the target value; m is the number of nodes in the hidden layer; n is the number of nodes in the input layer; b represents the bias term; is the bias term of the hidden layer; is the bias term of the output layer; w is the weight coefficient; is the weight between the hidden layer and the output layer, representing the connection weight between the jth hidden layer node and the kth output node, the second weight matrix; Is the weight between the input layer and the hidden layer, representing the connection weight between the i-th input node and the j-th hidden layer node, the first weight matrix.

4. The method according to claim 1, characterized in that The Euler particle model is used to determine the proppant distribution information within the fracture, including: The Euler particle model describes the flow characteristics of the proppant and fluid mixture based on the solid phase pressure model and kinetic theory; the fluid phase is the main phase and the proppant is the secondary phase, and the main phase and the secondary phase share a unified pressure field data in the numerical calculation; the subscripts p and q are used to refer to different phases, respectively, and the mass conservation and momentum conservation equations of the q phase are expressed as: In the formula, α q represents the volume fraction of q phase; v q represents the velocity vector of the q phase; v p represents the velocity vector of the p phase; ρ q represents the density of phase q; t represents time; represents the gradient operator; v qp 、v pq They represent the relative speed of phase q relative to phase p and the relative speed of phase p relative to phase q respectively; n represents the number of phases; represents the mass transfer from the p phase to the q phase; represents the mass transfer from phase q to phase p; g is the acceleration due to gravity; K pq is the momentum exchange coefficient between p and q phases; F q represents the external volume force of phase q; F lift,q represents the lift force of phase q; F wl,q represents the wall lubrication force of phase q; F vm,q represents the virtual mass force of the q phase; F td,q represents the drag force of the q phase; is the q-corresponding force-strain tensor, calculated by formula (9): In the formula, α q represents volume fraction; μ q and λ q represent the shear viscosity and bulk viscosity of the q phase, respectively; represents the unit tensor; represents the transpose of the velocity matrix; Determine v according to the direction of mass transfer between phases pq :

5. The method according to claim 1, characterized in that Determine the particle stress information based on the Euler particle model, including: In the Euler particle model, the stress of the proppant phase is derived by analogizing the random particle collision motion to the thermal motion of gas molecules. The particle phase temperature is introduced, and the particle motion fluctuation and temperature value are linearly related. The momentum conservation equation of the particle phase s is expressed as: In the formula, p s represents the pressure of the particle phase; t represents time; K ls represents the momentum exchange coefficient between the fluid phase and the particle phase; N represents the total number of phases contained in the system; α s represents the volume fraction of particle phase; ρ s represents the particle phase pressure; v s represents the particle phase velocity vector; represents the gradient operator; n represents the number of phases; v p represents the velocity vector of the p phase; v ls represents the relative velocity of the fluid phase l relative to the particle phase s; v sl It represents the relative velocity of particle phase s relative to fluid phase l; Represents the mass transfer from the fluid phase to the particle phase; represents the mass transfer from the particle phase to the fluid phase; F s represents the external volume force of the particle phase; F lift,s represents the lift force of the particle phase; F wl,s represents the wall lubrication force of the particle phase; F vm,s represents the virtual mass force of the particle phase; F td,s represents the drag force of the particle phase; g is the acceleration due to gravity, represents the stress-strain tensor of the particle phase; When the volume fraction is less than the critical value of accumulation, the particle phase pressure is calculated by the following formula: In the formula, e ss represents the particle collision restitution coefficient; g 0,ss represents the radial distribution function; Θ s represents the particle temperature; Stress-strain tensor of the particle phase It is expressed as: In the formula, μ s and λ s represent the shear viscosity and bulk viscosity of the particle phase, respectively; represents the unit tensor, d s represents the particle diameter; e s represents the elastic recovery coefficient of the particle phase; When the volume fraction of the particle phase reaches the critical value of stacking, the shear viscosity μ s It is expressed as: In the formula, represents the internal friction angle; parameters n and m are empirical constants; α s,min Represents the minimum value of particle phase concentration; α s,max Indicates the maximum value of the particle phase concentration.

6. The method according to claim 1, characterized in that For shale gas reservoirs with fractured areas after hydraulic fracturing, a matrix-microfracture-hydraulic fracture coupling model is used for production simulation. A single-pore medium model is used outside the fractured area, a dual-medium model is used for bedrock and microfractures in the fractured area, and an embedded discrete fracture model is used for hydraulic fractures, including: A shale gas seepage control equation is established, and the shale gas seepage control equation satisfies the following conditions: There are free gas and adsorbed gas in the shale matrix; adsorbed gas is not considered in microfracture and hydraulic fracture systems; free gas is real gas, and shale reservoirs can be compressed; gas adsorption obeys the isothermal adsorption model; the migration of gas in the shale matrix considers the Kinkenberg effect, and the migration of gas in microfracture and hydraulic fracture systems obeys Darcy's law; the adsorption, desorption and flow of shale gas are all isothermal processes; Gas flow satisfies the law of conservation of mass, and the mathematical expression is as follows: Where t represents time, Γ and Ω represent the boundary and region respectively; n represents the unit normal vector of the boundary Γ; A, F and q are the mass accumulation term, mass flow term and source and sink term in the region Ω respectively; among which, the mass accumulation term is expressed as: A=fr g +(1-φ)m (17) Where φ is the Lagrangian porosity considering the fluid-solid coupling effect; ρ g is the gas density; m represents the gas adsorption amount, which is a function of the gas pressure. The gas adsorption amount and gas density are determined by the following formula: Where p is the gas pressure; M is the gas molar mass; Z is the gas compressibility factor; R is the universal gas constant; T is the reservoir temperature; ρ r is the rock bulk density; ρ gstd is the gas density under standard conditions; p L is the Langmui pressure, which represents the gas pressure when the maximum adsorption capacity is half; V L is the Langmuir volume, which indicates the maximum adsorption capacity of the shale matrix; The mass flow term F in formula (16) is expressed as: Where k is the absolute permeability; μ is the gas viscosity; g is the gravitational acceleration; D is the reservoir depth, m; b is the Klinkenberg coefficient; c is the gas compressibility coefficient; D k is the real gas Knudsen diffusion coefficient; Establish the shale reservoir deformation control equation: During the mining process, the reservoir undergoes a small static linear elastic deformation. Under the sign convention that tensile stress is positive and compressive stress is negative, the reservoir stress balance equation is expressed as: Where σ is the total stress tensor; ρ b is the rock density; is the gradient operator; The effective stress equations of the unfractured area and the fractured area are: Where C represents the shale matrix elastic tensor, Pa; ε is the strain tensor; α is the Biot coefficient; I is the unit tensor; C up and K dr are the equivalent elastic tensor and bulk modulus of the fracture unit in the fracturing stimulation area; b l is the coupling coefficient; l represents the lth subgrid block in the matrix nested grid block within each fracture unit; p l represents the total pressure in the lth subunit; SRV represents the fracturing transformation area; The strain tensor ε and the displacement vector u satisfy the following relationship: In the formula, represents the gradient operator, and the superscript T represents the matrix transpose; Establish constitutive model; When the geometry of the hydraulic fracture consists of two parallel plates and the degree of closure of the propping section of the hydraulic fracture depends on the deformation of the proppant layer, the constitutive model is expressed as: in=in + -in - (29) In the formula, f s The nonlinear stress-strain relationship representing the deformation of the proppant layer; ε s is the proppant normal strain; u represents the relative displacement difference between the two fracture surfaces of the hydraulic fracture; d HF0 is the initial crack opening; n HF represents the crack normal vector; u + represents the displacement above the crack surface; u - represents the displacement below the crack surface; The closed constitutive model of the unsupported section of the hydraulic fracture satisfies the normal contact constraint of the free surface. The specific expression is as follows: d HF =d HF0 +u·n HF (31) Where, d HF is the remaining hydraulic fracture opening; E n is the normal contact penalty parameter; The matrix dynamic permeability is determined by the following formula: In the formula, k m represents the matrix permeability, k m0 represents the matrix permeability in the initial state, φ m represents the porosity, φ m0 represents the porosity in the initial state; The dynamic permeability of natural fractures is determined by the following formula: In the formula, k f is the bulk modulus of the natural fracture, k f0 is the bulk modulus of the natural fracture in the initial state; ε V is the volume strain of the grid block; K dr is the bulk modulus in the fracture stimulation area; According to the parallel plate theory and the cubic law, the dynamic permeability k of an unsupported fracture is HF It is expressed as: Based on the correlation between the permeability of the support joint and the effective stress, the generalized expression for determining the permeability of the support joint is as follows: k HF,s =f k (p s ) (35) In the formula, f k This is the relationship between fracture permeability and effective stress.

7. A shale gas reservoir fluid-solid coupling numerical simulation device considering proppant distribution, characterized in that: The device comprises: The computational model building module is configured to use a discrete element method to simulate the compaction process of the proppant layer under different closure stresses, and to calculate the permeability and conductivity of the propping seam by combining the pore network model, the lattice Boltzmann method and the artificial neural network model, and to obtain the relationship between the opening and conductivity of the propping section and the closure stress; The support distribution calculation module is configured to use computational fluid dynamics coupled discrete element method to achieve two-way coupling of proppant flow particles and fluid in hydraulic fractures, use computational fluid dynamics to simulate the momentum and mass conservation of sand-carrying fluid, use discrete element model to simulate the flow and contact of proppant particles, apply fluid force to the particle phase according to the drag model, realize the numerical simulation of proppant migration in hydraulic fractures, and obtain the actual distribution characteristics of proppant in hydraulic fractures; The numerical simulation module is configured to use a matrix-microfracture-hydraulic fracture coupling model to simulate production in shale gas reservoirs in the fractured zone after hydraulic fracturing, wherein a single-pore medium model is used outside the fractured zone, a dual medium model is used for bedrock and microfractures in the fractured zone, and an embedded discrete fracture model is used for hydraulic fractures; wherein the actual distribution characteristics of proppants and the relationship between the opening and conductivity of the propping section and the closure stress are considered in the hydraulic fractures; The computational model building module is further configured to: The discrete element method is used to simulate the compaction process of the proppant layer under different closing stresses and calculate the change in the proppant crack opening. The crack closure is achieved by setting a controllable movable plate above the proppant layer and applying a constant pressure. When the movable plate reaches a stable state, the pore structure and porosity inside the proppant accumulation are obtained by combining discrete element calculation data, image processing and digital core technology, and the pore network is extracted. The permeability and fracture conductivity of the proppant layer are calculated based on the improved pore network model; the conductivity calculation formula in the improved pore network model is: In the formula, g μ is the conductivity; k is the marker symbol of the subthroat; it represents different segments in the pore network; g k,μ is the conductivity of the subthroat k segment; L k is the length of the throat; μ is the fluid dynamic viscosity; a x is the displacement acceleration; q k is the volume flow rate; μ LBM is the fluid viscosity in the LBM method; ρ is the fluid density; The artificial neural network model is used to improve the computational efficiency of the improved pore network model. 8 . A non-transitory computer-readable storage medium storing instructions, which, when executed by a processor, executes the method according to claim 1 .

Images