Calculation method for migration of proppant in fracturing crack of horizontal well and global flow conductivity of crack

Abstract

The invention discloses a method for calculating the migration of a proppant in a horizontal well fracturing crack and the global flow conductivity of the crack, and relates to the field of unconventional oil and gas reservoir yield increase renovation, and the method comprises the following steps: (1) establishing a reservoir geologic model according to on-site logging exploration data; (2) establishing a hydraulic fracturing crack propagation fluid-solid coupling model; (3) establishing a horizontal well fracturing inter-cluster flow distribution model; (4) constructing a propping agent migration liquid-solid two-phase flow model in the hydraulic fracture; (5) solving the fluid-solid coupling model and the proppant migration model through loop iteration, and predicting the concentration distribution of the proppant in the hydraulic fracture; and (6) on the basis of the crack width and the proppant distribution result, calculating the hydraulic fracturing crack global flow conductivity. According to the method, the proppant migration model is solved by adopting a high-order discrete format, dynamic coupling simulation of hydraulic fracture propagation and proppant migration under the mine field scale is realized, the global flow conductivity of the hydraulic fracture is calculated, and reference is provided for evaluation and optimization of unconventional reservoir fracturing transformation.

Description

Technical Field

[0001] The present invention relates to the field of reservoir production enhancement and transformation, and in particular to a method for calculating proppant migration in a horizontal well hydraulic fracture and the global conductivity of the fracture. Background Art

[0002] Horizontal well drilling and hydraulic fracturing are key technologies for achieving efficient development of unconventional oil and gas resources. Fracture formation and proppant migration are two important processes in horizontal well fracturing. Fracture formation is the basis for reservoir production enhancement, and proppant transportation is the key to high and stable production. By predicting the distribution of proppant in the horizontal well fracturing cracks after fracturing construction, evaluating the sand transport results and conductivity of the fracturing network, and optimizing the parameters involved in fracturing formation and sand transport, the cost reduction and efficiency improvement of unconventional oil and gas resource development can be achieved.

[0003] Horizontal well fracturing mainly includes the stages of pre-fluid fracturing, sand-carrying fluid sand transportation, and displacement fluid anti-sinking. The production dynamics of the production well after fracturing mainly depends on the migration and placement characteristics of the proppant in the fracturing fracture. The migration of proppant in the fracturing fracture is a complex liquid-solid two-phase flow problem. At the same time, there is a significant coupling effect between the expansion of the fracturing fracture and the migration of proppant. Since indoor experiments cannot describe the dynamic migration process of proppant during the dynamic extension of fractures at the mineral scale, the simulation of proppant migration in the fracturing fracture is mostly carried out through numerical simulation. At present, the numerical simulation of proppant migration in hydraulic fractures can be divided into two categories: Euler-Lagrangian method and Euler-Euler method. The Euler-Lagrangian method is suitable for the study of proppant migration mechanism at the indoor scale. It is a problem of huge computational complexity when analyzing the migration of proppant in multi-fracture expansion of horizontal well fracturing in actual engineering. The Euler-Euler method is suitable for studying the transport and placement of proppant in hydraulic fractures at an engineering scale. However, the existing models are limited to simple planar hydraulic fractures or bifurcated fractures with fixed widths. There is a lack of high-precision numerical calculation methods for proppant migration models, making it difficult to predict the proppant migration and placement characteristics in the hydraulic fracture network of horizontal wells in fractured reservoirs with developed natural fractures.

[0004] Therefore, it is urgent to establish a method for calculating the migration of proppant in horizontal well fractures and the global conductivity of fractures. By establishing and coupling the fracture proppant expansion model and the fracture proppant migration model, a high-order numerical discrete format is introduced to solve the nonlinear hyperbolic equation of proppant migration, and the dynamic simulation of horizontal well fracturing and proppant migration is realized. According to the proppant distribution results in the fracture and the calculation model of the conductivity of the sand-filled fracture, the global conductivity of the fracture is evaluated, which provides a new idea for simulating the migration of proppant in the fracture of horizontal well fractures and optimizing the design of fracturing construction parameters. Summary of the invention

[0005] In order to overcome the defects existing in the prior art, the purpose of the present invention is to provide a method for calculating the migration of proppant in horizontal well fractures and the global conductivity of fractures, so as to solve the problems of simple fracture morphology, difficult coupling of fracture extension and proppant migration, and poor accuracy of model calculation results in the prior art, thereby accurately simulating the migration results of proppant in fractures in complex reservoirs, and optimizing the fracturing construction parameters according to the global conductivity evaluation results, so as to improve the reservoir production transformation effect and realize the cost-saving and efficiency-enhancing development of unconventional reservoir oil and gas resources.

[0006] The present invention provides a method for calculating proppant migration in horizontal well hydraulic fractures and the global conductivity of fractures. The technical solution is as follows:

[0007] S01. Establish reservoir geological model based on field logging exploration data;

[0008] S02. Establish a fluid-solid coupling model for hydraulic fracturing crack expansion:

[0009] S03, establishing a flow distribution model between horizontal well fracturing clusters;

[0010] S04. Construct a liquid-solid two-phase flow model for proppant migration in hydraulic fractures;

[0011] S05. Solve the fluid-solid coupling model and proppant migration model to predict the proppant concentration distribution in the hydraulic fracture;

[0012] S06. Calculate the global conductivity of hydraulic fracturing fractures based on the fracture width and proppant distribution results.

[0013] Furthermore, the step S01 establishes a reservoir geological model based on the on-site logging exploration data, specifically including:

[0014] (1) Analyze the well logging exploration data of the target block, determine the reservoir geological characteristics in the horizontal well fracturing section, and judge whether natural fractures are developed;

[0015] (2) If natural fractures are not developed in the target block, the model is simplified to a homogeneous reservoir geological model. If natural fractures are developed in the target block, it is necessary to characterize and model the natural fractures. A natural fracture distribution model in the reservoir is established by using statistical methods and geological exploration data. Natural fractures are treated as vertical plane fractures. The Monte Carlo method is used to determine the geometric center distribution of natural fractures in the reservoir plane, and the surface density is introduced to calculate the number of natural fractures:

[0016]

[0017] The distribution of natural fracture lengths is described by the normal distribution probability density function:

[0018]

[0019] Fisher distribution is used to characterize the trend characteristics of natural fractures in reservoirs:

[0020]

[0021] In the above formula, ρ nf represents the surface density of natural fractures; N nf Indicates the number of natural fractures; A res represents the reservoir area; l nf represents the distribution of natural fracture length; θ nf represents the distribution of natural fracture strike angles; l ave represents the average length of natural fractures; σ l represents the standard deviation of natural fracture length; θ ave represents the average strike angle of natural fractures; K θ represents the Fisher distribution constant; K rand Represents a random number between 0 and 1;

[0022] (3) Obtain the geological parameters, completion parameters, construction parameters and material parameters required for fracturing construction design.

[0023] Furthermore, the step S02 establishes a fluid-solid coupling model for hydraulic fracturing crack extension, specifically including:

[0024] (1) The relationship between hydraulic fracture unit displacement and stress is established based on the three-dimensional displacement discontinuity method:

[0025]

[0026] In the above formula, G represents the shear modulus; v represents the Poisson's ratio; σ represents the stress acting on the hydraulic fracture unit; D represents the displacement of the hydraulic fracture unit; I represents the kernel function; the subscripts N, H, and L represent the normal, dip, and strike directions along the hydraulic fracture unit, respectively; represents the local coordinates of the hydraulic fracture unit;

[0027] (2) The flow of fracturing fluid in hydraulic fractures is described by Poiseuille's law. Considering the effect of proppant on filtration, the mass conservation equation of fracturing fluid flow is:

[0028]

[0029] In the above formula, w represents the width of hydraulic fracture; t represents the time step; v f Indicates the flow rate of fracturing fluid; Q inj represents the fracturing fluid pumping flow rate; c represents the proppant volume concentration; Q leak Indicates the filtration loss flow of fracturing fluid; μ f Indicates fracturing fluid viscosity; p frepresents the fluid pressure in the hydraulic fracture; γ represents the fluid static pressure gradient; z represents the gravity direction vector;

[0030] (3) The equivalent stress factor at the tip of the hydraulic fracture is introduced as the extension criterion of the hydraulic fracture in the reservoir matrix:

[0031] K e ≥K IC

[0032]

[0033] Predicting the propagation direction of hydraulic fractures by using the maximum circumferential stress criterion:

[0034]

[0035] In the above formula, K e Indicates the equivalent stress factor at the crack tip; K IC represents the fracture toughness of reservoir rock; α represents the crack extension angle in the local coordinate system of the crack tip; K I Represents the tensile stress intensity factor; K II represents the shear stress intensity factor;

[0036] When natural fractures develop in the reservoir, the intersection criterion of hydraulic fractures and natural fractures is introduced to describe the complex expansion behavior of hydraulic fractures after encountering natural fractures. The criteria for judging whether hydraulic fractures pass through natural fractures, hydraulic fractures extend along natural fractures, and shear slip occurs in natural fractures are as follows:

[0037] Through the cracks:

[0038] Steering extension:

[0039] Shear slip: τ rγ |≥τ 0 +k nf (σ γ -p nf ),γ=-θorπ-θ

[0040] In the above formula, σ γ represents the circumferential stress in the intersection area; σ r represents the radial stress in the intersection area; τ rγ represents the shear stress in the intersection area; σ n represents normal in-situ stress; T represents the tensile strength of reservoir rock; τ 0 represents the cohesion of natural fractures; p nf Fluid pressure in natural fractures; p f,int represents the fluid pressure at the intersection point; k nfrepresents the friction coefficient of the natural fracture wall; r represents the polar diameter in the polar coordinate system; γ represents the polar angle in the polar coordinate system; K IC,NF It indicates the fracture toughness of hydraulic fractures extending along natural fractures.

[0041] Furthermore, the step S03 establishes a flow distribution model between horizontal well fracturing clusters, which specifically includes:

[0042] (1) Calculate the flow friction of the fracturing fluid in the horizontal wellbore and perforation holes:

[0043]

[0044] In the above formula, Δp w represents the horizontal wellbore friction; Δp p represents the friction resistance of the perforation hole; μ f represents the viscosity of the fracturing fluid; ρ f Indicates the density of fracturing fluid; Q w Indicates the flow rate in the horizontal wellbore; Q p represents the flow rate through the perforation hole; L represents the length of the fracturing section; C d Indicates the hole flow coefficient; n p Indicates the number of perforation holes; d p represents the perforation aperture;

[0045] (2) Based on the pressure balance principle, the balance equation of the fluid pressure at the heel of the horizontal wellbore is established:

[0046] p h,1 =p h,2 =…=p h,m

[0047]

[0048] In the above formula, p h Indicates the fluid pressure at the heel of the horizontal wellbore; Δp w Indicates the flow friction of fracturing fluid in the horizontal wellbore; Δp p It indicates the friction resistance of fracturing fluid flowing through the perforation holes; p in Indicates the fluid pressure at the entrance of the hydraulic fracture;

[0049] (3) Based on the relationship between the perforation cluster flow rate and the fluid pressure balance equation, combined with the volume conservation and Newton iteration method, the flow rate of each perforation cluster entering the horizontal wellbore is calculated.

[0050] Furthermore, the step S04 constructs a liquid-solid two-phase flow model for proppant migration in hydraulic fractures, specifically including:

[0051] (1) Calculate the density and viscosity of the mixed sand solution formed by mixing the fracturing fluid and the proppant:

[0052] ρ s =c p ρ p +(1-c p )ρ f

[0053]

[0054] In the above formula, ρ s represents the density of the sand mixing liquid; ρ p represents the proppant density; μ s Indicates the viscosity of the sand mixing liquid; c p represents the volume concentration of proppant; c p,max Indicates the maximum allowable volume concentration of proppant; represents the dimensionless proppant concentration;

[0055] (2) Based on the Euler multiphase flow theory, the mass conservation equation of the flow of sand mixing fluid and proppant in the fracture is established:

[0056]

[0057] The flow rate calculation formula of sand mixing fluid and proppant in hydraulic fracture is:

[0058]

[0059] The proppant blocking function expression in hydraulic fracture is:

[0060]

[0061] The formula for calculating the sedimentation velocity of the proppant in the sand mixing fluid under the influence of gravity is:

[0062] v set =v stokes ·C iner ·C conc ·C frac

[0063]

[0064] In the above formula, v s Indicates the flow velocity of the sand mixing liquid; v p Indicates the proppant flow velocity; v set Indicates proppant settling velocity; v stokes It represents the Stokes settling velocity of the proppant particles; represents the dimensionless sand concentration; B represents the proppant plugging function; d p Indicates the proppant particle size; C iner The correction factor that indicates the effect of particle inertia on sedimentation velocity; C conc Indicates the correction coefficient of the effect of proppant concentration on sedimentation velocity; Cfrac It represents the correction coefficient of the effect of crack wall on settlement velocity.

[0065] Furthermore, the step S05 solves the fluid-solid coupling model and the proppant migration model to predict the proppant concentration distribution in the hydraulic fracture, specifically including:

[0066] (1) Combining the relationship between displacement and stress of hydraulic fracture units and the mass conservation equation of fluid flow in hydraulic fractures, a set of nonlinear fluid-solid coupling equations for three-dimensional hydraulic fracture expansion is established:

[0067]

[0068] In the above formula, A M represents the boundary influence coefficient moment; I represents the unit matrix; C f represents the fluid flow coefficient matrix; σ M represents the comprehensive stress matrix acting on the hydraulic fracture unit; Qs represents the fluid flow source matrix; w represents the hydraulic fracture width matrix; p f represents the fluid pressure matrix in the hydraulic fracture;

[0069] Before the proppant is pumped, the flow field in the above coupling equations adopts the mass conservation equation of the fracturing fluid. After the proppant is pumped, the flow field in the above coupling equations adopts the mass conservation equation of the sand mixing fluid. The nonlinear fluid-solid coupling equations are iteratively solved using Newton iteration to obtain the hydraulic fracture width and the fluid pressure in the fracture at the current moment. The iterative calculation format is:

[0070]

[0071] In the above formula, n represents the iteration step; Δw represents the iteration amount of hydraulic fracture width; Δp f represents the fluid pressure iteration amount;

[0072] (2) After solving the nonlinear fluid-solid coupling equations to obtain the hydraulic fracture width and the fluid pressure in the fracture, the flow velocity of the sand mixing fluid is calculated first, and then the sedimentation velocity of the proppant is calculated. The migration velocity of the proppant in the hydraulic fracture during the fracturing process is obtained by superposition calculation:

[0073] v p =B(w,d p )×(v f +v set )

[0074] (3) Substitute the proppant migration velocity into the proppant mass conservation equation to solve the distribution of proppant concentration in the fracture during the hydraulic fracture expansion process. First, the convection term and source term of the mass conservation equation are processed by the finite volume method:

[0075]

[0076]

[0077] The QUICK upwind format is further used for numerical discretization:

[0078]

[0079] The right side expression of the proppant flow mass conservation equation can be rearranged to obtain:

[0080]

[0081] In the above formula, dA represents the area of ​​the control unit; Δx and Δy represent the length and width of the control unit respectively; subscripts n, s, w, and e represent the upper, lower, left, and right boundaries of the control unit respectively; subscript C represents the center of the control unit; subscripts N, S, W, and E represent the centers of the adjacent units on the upper, lower, left, and right sides of the control unit respectively; NN, SS, WW, and EE represent the centers of the adjacent units on the upper side of N, the lower side of S, the left side of W, and the right side of E respectively; max() represents the maximum value function; A p represents the proppant flow coefficient matrix; represents the dimensionless concentration matrix of proppant; Q p,s represents the proppant source term matrix;

[0082] The fourth-order Runge-Kutta algorithm is introduced to discretize the unsteady-state term of the proppant flow mass conservation equation:

[0083]

[0084] In the above formula, t represents time; Δt represents the time step; F represents the right side term of the mass conservation equation; k1, k2, k3, k4 represent intermediate variables;

[0085] Under the condition of obtaining the hydraulic fracture width and fluid pressure, the concentration distribution of the proppant in the fracture during the hydraulic fracture expansion process can be calculated through the high-order discrete format of the above proppant flow mass conservation equation in the one-way coupling framework.

[0086] Furthermore, the step S06 calculates the global conductivity of the hydraulic fracturing fracture based on the fracture width and proppant distribution results, specifically including:

[0087] (1) According to the distribution results of proppant concentration in the hydraulic fracture unit after the fracture is closed after pumping, the hydraulic fracture unit is divided into two categories: unfilled fracture unit and sand-filled fracture unit:

[0088]

[0089] In the above formula, c prepresents the proppant volume concentration; represents the dimensionless proppant concentration; c p,max represents the maximum allowable volume fraction of proppant; c c represents the critical volume concentration;

[0090] (2) For unfilled fracture units, the fracture walls will gradually close and contact after the pump is stopped during construction. The critical closure width and contact permeability considering the roughness of the fracture wall are introduced to calculate the conductivity:

[0091] C 1 =w c ·k c

[0092] In the above formula, C 1 represents the conductivity of the unfilled fracture unit; w c represents the critical closing gap width; k c represents the contact permeability;

[0093] (3) For the sand-filled fracture unit filled with proppant, first calculate the porosity:

[0094]

[0095] The permeability of the sand-filled fracture unit is calculated according to the Kozeny-Carman model:

[0096]

[0097] According to the stress-displacement relationship of the hydraulic fracture unit and the fracture width distribution of the unfilled fracture unit, the fracture width of the sand-filled fracture unit after the pump is stopped and the fracture is closed is calculated:

[0098]

[0099] Combined with the permeability and width of the sand-filled fracture unit, the conductivity of the sand-filled fracture unit is calculated:

[0100] C 2 =w p ·k p

[0101] In the above formula, φ p represents the porosity of the sand-filled fracture unit; k p represents the permeability of the sand-filled fracture unit; V s Table proppant equivalent thickness; V b Indicates the proppant filling thickness; w 0 Indicates the width of the hydraulic fracture unit when the pump is stopped; w p k represents the width of the sand-filled fracture unit after the fracture is closed after pumping is stopped; K represents the Kozeny-Carman constant; SV Represents the surface area of ​​solid particles per unit volume; σ c represents the comprehensive stress acting on the crack unit; p f represents the fluid pressure in the hydraulic fracture; σ p represents the deformation stress of the proppant filling layer; σ f represents the contact stress between the proppant filling layer and the fracture wall; σ n represents the normal ground stress acting on the crack wall; C 2 It represents the conductivity of sand-filled fracture unit;

[0102] (4) Combining the conductivity results of unfilled fracture units with those of filled fracture units, the global conductivity distribution of hydraulic fracturing fractures is obtained, which can provide a reference for the design and optimization of reservoir hydraulic fracturing construction.

[0103] Compared with the prior art, the present invention has the following beneficial effects:

[0104] (1) The present invention can realize the coupled simulation of the dynamic expansion of the horizontal well hydraulic fracture and the proppant migration, and fully consider the mutual influence between the geometric morphology of the hydraulic fracture and the proppant transport state;

[0105] (2) The present invention introduces high-order numerical discrete formats such as the QUICK upwind format and the 4th-order Runge-Kutta algorithm to solve the nonlinear hyperbolic equation of proppant migration, which can accurately predict the migration and placement results of proppant in the fracturing fracture;

[0106] (3) The present invention can evaluate the global conductivity of complex fractures in horizontal well fracturing in fractured reservoirs, and provide a numerical simulation method for on-site fracturing construction design and parameter optimization. BRIEF DESCRIPTION OF THE DRAWINGS

[0107] In order to more clearly illustrate the technical solutions of the embodiments of the present invention, the drawings in the embodiments will be briefly introduced below. It should be understood that the following drawings only show certain embodiments of the present invention and therefore should not be regarded as limiting the scope. For ordinary technicians in this field, other related drawings can be obtained based on these drawings without creative work.

[0108] Figure 1 Schematic diagram of the implementation process of the present invention

[0109] Figure 2 Schematic diagram of fracture development interpretation results of reservoir imaging logging

[0110] Figure 3 Schematic diagram of the development trend of natural fractures in the plane

[0111] Figure 4 Schematic diagram of the modeling results of natural fractures in the reservoir

[0112] Figure 5 Schematic diagram of the pumping procedure for the embodiment

[0113] Figure 6 Schematic diagram of the solution process for coupling fracture propagation and proppant transport in horizontal wells

[0114] Figure 7 The three-dimensional fracture network extension trajectory and fracture width distribution results at the time of pump stop

[0115] Figure 8 The result diagram of proppant concentration distribution in the three-dimensional fracture network when the pump is stopped

[0116] Fig. 9 Schematic diagram of sand-filled fracture units and unfilled fracture units

[0117] Fig.10 The distribution result of the conductivity of the whole fracture area of ​​the horizontal well is shown in the figure Specific embodiments

[0118] In order to make the purpose, technical solution and advantages of the embodiments of the present invention clearer, the technical solution of the embodiments of the present invention will be clearly and completely described below in conjunction with the drawings of the embodiments of the present invention. Obviously, the described embodiments are part of the embodiments of the present invention, not all of the embodiments. Based on the described embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0119] Reference Figure 1 The present invention proposes a method for calculating proppant migration and global fracture conductivity in horizontal well hydraulic fractures, which mainly includes the following 6 steps:

[0120] S01. Establish reservoir geological model based on field logging exploration data;

[0121] S02. Establish a fluid-solid coupling model for hydraulic fracturing crack expansion:

[0122] S03, establishing a flow distribution model between horizontal well fracturing clusters;

[0123] S04. Construct a liquid-solid two-phase flow model for proppant migration in hydraulic fractures;

[0124] S05. Solve the fluid-solid coupling model and proppant migration model to predict the proppant concentration distribution in the hydraulic fracture;

[0125] S06. Calculate the global conductivity of hydraulic fracturing fractures based on the fracture width and proppant distribution results.

[0126] Specifically, step S01 of establishing a reservoir geological model based on on-site logging exploration data includes:

[0127] (1) Analyze the logging exploration data of the target block, determine the geological characteristics of the reservoir in the horizontal well fracturing section, and judge whether natural fractures are developed. Figure 2 , the natural fractures in the reservoir are developed, and the distribution of natural fractures in the plane is further obtained. Figure 3 .

[0128] (2) If natural fractures are not developed in the target block, the model is simplified to a homogeneous reservoir geological model. If natural fractures are developed in the target block, it is necessary to characterize and model the natural fractures. A natural fracture distribution model in the reservoir is established by using statistical methods and geological exploration data. Natural fractures are treated as vertical plane fractures. The Monte Carlo method is used to determine the geometric center distribution of natural fractures in the reservoir plane, and the surface density is introduced to calculate the number of natural fractures:

[0129]

[0130] The distribution of natural fracture lengths is described by the normal distribution probability density function:

[0131]

[0132] Fisher distribution is used to characterize the trend characteristics of natural fractures in reservoirs:

[0133]

[0134] In the above formula, ρ nf represents the surface density of natural fractures; N nf Indicates the number of natural fractures; A res represents the reservoir area; l nf represents the distribution of natural fracture length; θ nf represents the distribution of natural fracture strike angles; l ave represents the average length of natural fractures; σ l represents the standard deviation of natural fracture length; θ ave represents the average strike angle of natural fractures; K θ represents the Fisher distribution constant; K rand Represents a random number between 0 and 1.

[0135] In this example, the natural fractures in the selected reservoir are developed, and the natural fracture modeling in the reservoir is carried out by combining the above model with the natural fracture development trend data, where the natural fracture density is 0.18 / m 2 , the average length of natural fractures is 10m, the strike angles of natural fractures are 50° and 120°, and the distribution model of natural fractures in the reservoir is as follows: Figure 4 shown.

[0136] (3) Obtain the geological parameters, completion parameters, construction parameters and material parameters required for fracturing construction design.

[0137] In this embodiment, the required geological parameters, completion parameters, construction parameters and material parameters are shown in Table 1.

[0138] Table 1 Parameters of the embodiments

[0139] Serial number parameter Numeric Serial number parameter 1 Maximum horizontal principal stress / MPa 48.3 14 <![CDATA[Fracture toughness of natural fractures / MPa·m 0.5 > 2.0 2 Minimum horizontal principal stress / MPa 41.6 15 Wellbore inner diameter / mm 104.8 3 Vertical stress / MPa 46.1 16 Perforation aperture / mm 11 4 <![CDATA[Reservoir stress gradient / MPa·m -1 > 0.024 17 <![CDATA[Number of perforations / holes per cluster -1 > 8 5 Reservoir thickness / m 29.7 18 Perforation cluster number / cluster 4 6 Interlayer stress difference / MPa 3.10 19 Cluster spacing / m 15 7 Young's modulus / GPa 29.1 20 <![CDATA[Pump injection volume / m 3 > 576 8 Poisson's ratio / - 0.21 21 <![CDATA[Construction displacement / m 3 ·min -1 > 12 9 <![CDATA[Fracture toughness of rock / MPa·m 0.5 > 5.0 22 Fracturing fluid viscosity / mPa·s 5 10 Rock tensile strength / MPa 2.30 23 <![CDATA[Fracturing fluid density / kg·m -3 > 1000 11 Reservoir pore pressure / MPa 29.0 24 Proppant particle size / μm 220 12 Natural fracture permeability / mD 150 25 <![CDATA[Proppant density / kg·m -3 > 2560 13 Natural crack friction coefficient / - 0.70 26 <![CDATA[Maximum allowable concentration / m 3 ·m -3 > 0.585

[0140] Specifically, step S02 of establishing a fluid-solid coupling model for hydraulic fracturing crack extension includes:

[0141] (1) The relationship between hydraulic fracture unit displacement and stress is established based on the three-dimensional displacement discontinuity method:

[0142]

[0143] In the above formula, G represents the shear modulus; v represents the Poisson's ratio; σ represents the stress acting on the hydraulic fracture unit; D represents the displacement of the hydraulic fracture unit; I represents the kernel function; the subscripts N, H, and L represent the normal, dip, and strike directions along the hydraulic fracture unit, respectively; Represents the local coordinates of the hydraulic fracture element.

[0144] (2) The flow of fracturing fluid in hydraulic fractures is described by Poiseuille's law. Considering the effect of proppant on filtration, the mass conservation equation of fracturing fluid flow is:

[0145]

[0146] In the above formula, w represents the width of hydraulic fracture; t represents the time step; v f Indicates the flow rate of fracturing fluid; Q inj represents the fracturing fluid pumping flow rate; c represents the proppant volume concentration; Q leak Indicates the filtration loss flow of fracturing fluid; μ f Indicates fracturing fluid viscosity; p f represents the fluid pressure in the hydraulic fracture; γ represents the fluid static pressure gradient; z represents the gravity direction vector.

[0147] (3) The equivalent stress factor at the tip of the hydraulic fracture is introduced as the extension criterion of the hydraulic fracture in the reservoir matrix:

[0148] K e ≥K IC

[0149]

[0150] Predicting the propagation direction of hydraulic fractures by using the maximum circumferential stress criterion:

[0151]

[0152] In the above formula, K e Indicates the equivalent stress factor at the crack tip; K IC represents the fracture toughness of reservoir rock; α represents the crack extension angle in the local coordinate system of the crack tip; K I Represents the tensile stress intensity factor; K II represents the shear stress intensity factor.

[0153] When natural fractures develop in the reservoir, the intersection criterion of hydraulic fractures and natural fractures is introduced to describe the complex expansion behavior of hydraulic fractures after encountering natural fractures. The criteria for judging whether hydraulic fractures pass through natural fractures, hydraulic fractures extend along natural fractures, and shear slip occurs in natural fractures are as follows:

[0154] Through the cracks:

[0155] Steering extension:

[0156] Shear slip: |τ rγ |≥τ 0 +k nf (σ γ -p nf ),γ=-θorπ-θ

[0157] In the above formula, σ γ represents the circumferential stress in the intersection area; σ r represents the radial stress in the intersection area; τ rγ represents the shear stress in the intersection area; σ n represents normal in-situ stress; T represents the tensile strength of reservoir rock; τ 0 represents the cohesion of natural fractures; p nf Fluid pressure in natural fractures; p f,int represents the fluid pressure at the intersection point; k nf represents the friction coefficient of the natural fracture wall; r represents the polar diameter in the polar coordinate system; γ represents the polar angle in the polar coordinate system; K IC,NF It indicates the fracture toughness of hydraulic fractures extending along natural fractures.

[0158] In this embodiment, the horizontal wellbore is in the direction of the horizontal minimum principal stress, the number of perforation clusters is 4, and the initial fracture is along the direction of the horizontal maximum principal stress. Combining the displacement-stress relationship and the flow mass conservation equation, a fluid-solid coupling nonlinear equation group is constructed. Its Newton iteration format is shown below. The width of the hydraulic fracture and the fluid pressure are obtained through coupling solution.

[0159]

[0160] Specifically, step S03 of establishing a flow distribution model between horizontal well fracturing clusters includes:

[0161] (1) Calculate the flow friction of the fracturing fluid in the horizontal wellbore and perforation holes:

[0162]

[0163] In the above formula, Δp w represents the horizontal wellbore friction; Δp p represents the friction resistance of the perforation hole; μ f represents the viscosity of the fracturing fluid; ρ f Indicates the density of fracturing fluid; Q w Indicates the flow rate in the horizontal wellbore; Q p represents the flow rate through the perforation hole; L represents the length of the fracturing section; C d Indicates the hole flow coefficient; n p Indicates the number of perforation holes; d p Represents the perforation aperture.

[0164] (2) Based on the pressure balance principle, the balance equation of the fluid pressure at the heel of the horizontal wellbore is established:

[0165] p h,1 =p h,2 =…=p h,m

[0166]

[0167] In the above formula, p h Indicates the fluid pressure at the heel of the horizontal wellbore; Δp w Indicates the flow friction of fracturing fluid in the horizontal wellbore; Δp p It indicates the friction resistance of fracturing fluid flowing through the perforation holes; p in Represents the fluid pressure at the entrance of the hydraulic fracture.

[0168] p h,1 =p h,2 =…=p h,m

[0169]

[0170] (3) Based on the relationship between the perforation cluster flow rate and the fluid pressure balance equation, combined with the volume conservation and Newton iteration method, the flow rate of each perforation cluster entering the horizontal wellbore is calculated.

[0171] In this embodiment, the number of perforation clusters is 4, and the balance equation of the fluid pressure at the heel of the horizontal wellbore is:

[0172] p h,1 =p h,2 =p h,3 =ph,4

[0173]

[0174] The volume conservation formula of the flow rate into each perforation cluster and the total operation displacement of horizontal well fracturing is:

[0175] Q p,1 +Q p,2 +Q p,3 +Q p,4 =Q all

[0176] Construct Newton iteration calculation format:

[0177]

[0178] The dynamic flow rate entering each perforation cluster during the horizontal well fracturing process is obtained by iterative calculation based on the above formula.

[0179] Specifically, step S04 constructs a liquid-solid two-phase flow model for proppant migration in hydraulic fractures, including:

[0180] (1) Calculate the density and viscosity of the mixed sand solution formed by mixing the fracturing fluid and the proppant:

[0181] ρ s =c p ρ p +(1-c p )ρ f

[0182]

[0183] In the above formula, ρ s represents the density of the sand mixing liquid; ρ p represents the proppant density; μ s Indicates the viscosity of the sand mixing liquid; c p represents the volume concentration of proppant; c p,max Indicates the maximum allowable volume concentration of proppant; Represents the dimensionless proppant concentration.

[0184] (2) Based on the Euler multiphase flow theory, the mass conservation equation of the flow of sand mixing fluid and proppant in the fracture is established:

[0185]

[0186] The flow rate calculation formula of sand mixing fluid and proppant in hydraulic fracture is:

[0187]

[0188] The proppant blocking function expression in hydraulic fracture is:

[0189]

[0190] The formula for calculating the sedimentation velocity of the proppant in the sand mixing fluid under the influence of gravity is:

[0191] v set =v stokes ·C iner ·C conc ·C frac

[0192]

[0193] In the above formula, v s Indicates the flow velocity of the sand mixing liquid; v p Indicates the proppant flow velocity; v set Indicates proppant settling velocity; v stokes It represents the Stokes settling velocity of the proppant particles; represents the dimensionless sand concentration; B represents the proppant plugging function; d p Indicates the proppant particle size; C iner The correction factor that indicates the effect of particle inertia on sedimentation velocity; C conc Indicates the correction coefficient of the effect of proppant concentration on sedimentation velocity; C frac It represents the correction coefficient of the effect of crack wall on settlement velocity.

[0194] Specifically, step S05 solves the fluid-solid coupling model and the proppant migration model to predict the proppant concentration distribution in the hydraulic fracture, including:

[0195] Furthermore, the step S05 solves the fluid-solid coupling model and the proppant migration model to predict the proppant concentration distribution in the hydraulic fracture, specifically including:

[0196] (1) Combining the relationship between displacement and stress of hydraulic fracture units and the mass conservation equation of fluid flow in hydraulic fractures, a set of nonlinear fluid-solid coupling equations for three-dimensional hydraulic fracture expansion is established:

[0197]

[0198] In the above formula, A M represents the boundary influence coefficient moment; I represents the unit matrix; C f represents the fluid flow coefficient matrix; σ M represents the comprehensive stress matrix acting on the hydraulic fracture unit; Qs represents the fluid flow source matrix; w represents the hydraulic fracture width matrix; p f represents the fluid pressure matrix in the hydraulic fracture;

[0199] Before the proppant is pumped, the flow field in the above coupling equations adopts the mass conservation equation of the fracturing fluid. After the proppant is pumped, the flow field in the above coupling equations adopts the mass conservation equation of the sand mixing fluid. The nonlinear fluid-solid coupling equations are iteratively solved using Newton iteration to obtain the hydraulic fracture width and the fluid pressure in the fracture at the current moment. The iterative calculation format is:

[0200]

[0201] In the above formula, n represents the iteration step; Δw represents the iteration amount of hydraulic fracture width; Δp f represents the fluid pressure iteration amount;

[0202] (2) After solving the nonlinear fluid-solid coupling equations to obtain the hydraulic fracture width and the fluid pressure in the fracture, the flow velocity of the sand mixing fluid is calculated first, and then the sedimentation velocity of the proppant is calculated. The migration velocity of the proppant in the hydraulic fracture during the fracturing process is obtained by superposition calculation:

[0203] v p =B(w,d p )×(v f +v set )

[0204] (3) Substitute the proppant migration velocity into the proppant mass conservation equation to solve the distribution of proppant concentration in the fracture during the hydraulic fracture expansion process. First, the convection term and source term of the mass conservation equation are processed by the finite volume method:

[0205]

[0206]

[0207] The QUICK upwind format is further used for numerical discretization:

[0208]

[0209] The right side expression of the proppant flow mass conservation equation can be rearranged to obtain:

[0210]

[0211] In the above formula, dA represents the area of ​​the control unit; Δx and Δy represent the length and width of the control unit respectively; subscripts n, s, w, and e represent the upper, lower, left, and right boundaries of the control unit respectively; subscript C represents the center of the control unit; subscripts N, S, W, and E represent the centers of the adjacent units on the upper, lower, left, and right sides of the control unit respectively; NN, SS, WW, and EE represent the centers of the adjacent units on the upper side of N, the lower side of S, the left side of W, and the right side of E respectively; max() represents the maximum value function; A prepresents the proppant flow coefficient matrix; represents the dimensionless concentration matrix of proppant; Q p,s represents the proppant source term matrix;

[0212] The fourth-order Runge-Kutta algorithm is introduced to discretize the unsteady-state term of the proppant flow mass conservation equation:

[0213]

[0214] In the above formula, t represents time; Δt represents the time step; F represents the right side term of the mass conservation equation; k1, k2, k3, k4 represent intermediate variables;

[0215] Under the condition of obtaining the hydraulic fracture width and fluid pressure, the concentration distribution of the proppant in the fracture during the hydraulic fracture expansion process can be calculated through the high-order discrete format of the above proppant flow mass conservation equation in the one-way coupling framework.

[0216] In this embodiment, the schematic diagram of the horizontal well fracturing pumping program is as follows: Figure 5 As shown in the figure, the fracturing and proppant transport simulation is carried out by using a constant displacement and step-by-step sand addition method, including the main stages of pre-fluid fracturing, sand-carrying fluid sanding, and pump stop. The model solution process is as follows Figure 6 As shown in FIG. 1 , a sequential coupling method is adopted. First, the fracture width and the fluid pressure in the fracture are obtained by calculating the fracture extension fluid-solid coupling equations. Then, the migration velocity of the proppant in the hydraulic fracture is calculated. The proppant concentration distribution in the hydraulic fracture is obtained by introducing the numerical discrete format of proppant migration. After the pump is stopped in this embodiment, the three-dimensional extension trajectory and spatial distribution of the hydraulic fracture are shown in FIG. Figure 7 As shown in the figure, a complex hydraulic fracture network is formed under the interference conditions between natural fractures and interfracture stress. Figure 8 It is the distribution of proppant concentration in the fracturing fracture network after the pump is stopped. Due to the non-uniform distribution of fracture width, the distribution of proppant in the three-dimensional fracturing fracture is complex. Some fracturing fractures are not filled with proppant, which will cause differences in the conductivity and production contribution of each perforation cluster.

[0217] Specifically, step S06 calculates the global conductivity of hydraulic fractures based on the fracture width and proppant distribution results, including:

[0218] (1) According to the distribution results of proppant concentration in the hydraulic fracture unit after the fracture is closed after pumping, the hydraulic fracture unit is divided into two categories: unfilled fracture unit and sand-filled fracture unit:

[0219]

[0220] In the above formula, c p represents the proppant volume concentration; represents the dimensionless proppant concentration; cp,max represents the maximum allowable volume fraction of proppant; c c represents the critical volume concentration.

[0221] (2) For unfilled fracture units, the fracture walls will gradually close and contact after the pump is stopped during construction. The critical closure width and contact permeability considering the roughness of the fracture wall are introduced to calculate the conductivity:

[0222] C 1 =w c ·k c

[0223] In the above formula, C 1 represents the conductivity of the unfilled fracture unit; w c represents the critical closing gap width; k c Represents the contact permeability.

[0224] (3) For the sand-filled fracture unit filled with proppant, first calculate the porosity:

[0225]

[0226] The permeability of the sand-filled fracture unit is calculated according to the Kozeny-Carman model:

[0227]

[0228] According to the stress-displacement relationship of the hydraulic fracture unit and the fracture width distribution of the unfilled fracture unit, the fracture width of the sand-filled fracture unit after the pump is stopped and the fracture is closed is calculated:

[0229]

[0230] Combined with the permeability and width of the sand-filled fracture unit, the conductivity of the sand-filled fracture unit is calculated:

[0231] C 2 =w p ·k p

[0232] In the above formula, φ p represents the porosity of the sand-filled fracture unit; k p represents the permeability of the sand-filled fracture unit; V s Table proppant equivalent thickness; V b Indicates the proppant filling thickness; w 0 Indicates the width of the hydraulic fracture unit when the pump is stopped; w p k represents the width of the sand-filled fracture unit after the fracture is closed after pumping is stopped; K represents the Kozeny-Carman constant; S V Represents the surface area of ​​solid particles per unit volume; σ crepresents the comprehensive stress acting on the crack unit; p f represents the fluid pressure in the hydraulic fracture; σ p represents the deformation stress of the proppant filling layer; σ f represents the contact stress between the proppant filling layer and the fracture wall; σ n represents the normal ground stress acting on the crack wall; C 2 Represents the conductivity of the sand-filled fracture unit.

[0233] (4) Combining the conductivity results of unfilled fracture units with those of filled fracture units, the global conductivity distribution of hydraulic fracturing fractures is obtained, which can provide a reference for the design and optimization of reservoir hydraulic fracturing construction.

[0234] In this embodiment, according to Figure 8 The distribution of proppant in the fracturing fracture network at the time of pump stop, critical volume concentration c c The value is 0.01. First, the hydraulic fracture units are divided into two categories: unfilled fracture units and sand-filled fracture units. Fig. 9 As shown. For the unfilled fracture unit, the conductivity can be calculated based on the closed fracture width and contact permeability, where the closed fracture width and contact permeability can be measured based on experimental data. In this embodiment, the average values ​​are 0.57 mm and 3.46 D. For the sand-filled fracture unit filled with proppant, the fracture width and proppant concentration distribution of the fracturing fracture network after pumping is stopped, as shown Figure 7 and Figure 8 As shown in , the permeability and closed fracture width of the sand-filled fracture unit are first calculated, and then the conductivity of the proppant-filled area is calculated. The closed fracture width is iteratively calculated by the fluid pressure condition, stress condition and fracture width distribution in the fracture. The natural fractures, unfilled fractures and sand-filled fractures are superimposed to obtain the global conductivity distribution of hydraulic fractures, as shown in Fig.10 As shown in the figure, the simulation found that the accumulation of proppant at the front edge of the hydraulic fracture makes the conductivity of the front end of the hydraulic fracture larger, with a maximum value of about 11.30D·cm, and the conductivity of some areas not filled with proppant is extremely small. In order to improve the conductivity of the overall fracturing fracture to promote production increase, it is recommended to appropriately increase the construction displacement or fracturing fluid viscosity based on the existing pumping scheme to improve the fracture width to promote uniform expansion of the fracture and long-distance transport of proppant.

[0235] The above description is only a preferred embodiment of the present invention and does not limit the present invention in any form. Although the present invention has been disclosed as a preferred embodiment as above, it is not intended to limit the technical contents to make slight changes or modifications. Without departing from the scope of the technical solution of the present invention, the technical contents disclosed above can be used to make slight changes or modifications to equivalent embodiments of equivalent changes. 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 content of the technical solution of the present invention are still within the scope of the technical solution of the present invention.

Claims

1. A method for calculating proppant migration and global fracture conductivity in horizontal well hydraulic fractures, characterized in that: The following steps are involved: S01. Establish reservoir geological model based on field logging exploration data; S02. Establish a fluid-solid coupling model for hydraulic fracturing crack expansion: S03, establishing a flow distribution model between horizontal well fracturing clusters; S04. Construct a liquid-solid two-phase flow model for proppant migration in hydraulic fractures; S05. Solve the fluid-solid coupling model and proppant migration model to predict the proppant concentration distribution in the hydraulic fracture; S06. Calculate the global conductivity of hydraulic fracturing fractures based on the fracture width and proppant distribution results.

2. The method for calculating proppant migration and global fracture conductivity in horizontal well hydraulic fractures according to claim 1, characterized in that: The detailed steps of step S01 are: (1) Analyze the well logging exploration data of the target block, determine the reservoir geological characteristics in the horizontal well fracturing section, and judge whether natural fractures are developed; (2) If natural fractures are not developed in the target block, the model is simplified to a homogeneous reservoir geological model. If natural fractures are developed in the target block, it is necessary to characterize and model the natural fractures. A natural fracture distribution model in the reservoir is established by using statistical methods and geological exploration data. Natural fractures are treated as vertical plane fractures. The Monte Carlo method is used to determine the geometric center distribution of natural fractures in the reservoir plane, and the surface density is introduced to calculate the number of natural fractures: The distribution of natural fracture lengths is described by the normal distribution probability density function: Fisher distribution is used to characterize the trend characteristics of natural fractures in reservoirs: In the above formula, ρ nf represents the natural fracture surface density; N nf Indicates the number of natural fractures; A res Represents the reservoir area; l nf represents the distribution of natural fracture length; θ nf represents the distribution of natural fracture strike angles; l ave represents the average length of natural fractures; σ l represents the standard deviation of natural fracture length; θ ave represents the average strike angle of natural fractures; K θ represents the Fisher distribution constant; K rand Represents a random number between 0 and 1; (3) Obtain the geological parameters, completion parameters, construction parameters and material parameters required for fracturing construction design.

3. The method for calculating proppant migration and global fracture conductivity in horizontal well hydraulic fractures according to claim 1, characterized in that: The detailed steps of step S02 are: (1) The relationship between hydraulic fracture unit displacement and stress is established based on the three-dimensional displacement discontinuity method: In the above formula, G represents the shear modulus; v represents the Poisson's ratio; σ represents the stress acting on the hydraulic fracture unit; D represents the displacement of the hydraulic fracture unit; I represents the kernel function; the subscripts N, H, and L represent the normal, dip, and strike directions along the hydraulic fracture unit, respectively; represents the local coordinates of the hydraulic fracture unit; (2) The flow of fracturing fluid in hydraulic fractures is described by Poiseuille's law. Considering the effect of proppant on filtration, the mass conservation equation of fracturing fluid flow is: In the above formula, w represents the width of hydraulic fracture; t represents the time step; v f Indicates the flow rate of fracturing fluid; Q inj represents the fracturing fluid pumping flow rate; c represents the proppant volume concentration; Q leak Indicates the filtration loss flow of fracturing fluid; μ f Indicates fracturing fluid viscosity; p f represents the fluid pressure in the hydraulic fracture; γ represents the fluid static pressure gradient; z represents the gravity direction vector; (3) The equivalent stress factor at the tip of the hydraulic fracture is introduced as the extension criterion of the hydraulic fracture in the reservoir matrix: K e ≥K IC Predicting the hydraulic fracture propagation direction by the maximum circumferential stress criterion: In the above formula, K e Indicates the equivalent stress factor at the crack tip; K IC represents the fracture toughness of reservoir rock; α represents the crack extension angle in the local coordinate system of the crack tip; K I Represents the tensile stress intensity factor; K II represents the shear stress intensity factor; When natural fractures develop in the reservoir, the intersection criterion of hydraulic fractures and natural fractures is introduced to describe the complex expansion behavior of hydraulic fractures after encountering natural fractures. The criteria for judging whether hydraulic fractures pass through natural fractures, hydraulic fractures extend along natural fractures, and shear slip occurs in natural fractures are as follows: Through the cracks: Steering extension: Shearing slip: |τ rγ | ≥ τ0 + k nf (σ γ - p nf ), γ = -θ or π - θ In the above formula, σ γ represents the circumferential stress in the intersection area; σ r represents the radial stress in the intersection area; τ rγ represents the shear stress in the intersection area; σ n represents the normal in-situ stress; T represents the tensile strength of reservoir rock; τ0 represents the cohesive force of natural fractures; p nf Fluid pressure in natural fractures; p f,int represents the fluid pressure at the intersection point; k nf represents the friction coefficient of the natural fracture wall; r represents the polar diameter in the polar coordinate system; γ represents the polar angle in the polar coordinate system; K IC,NF It indicates the fracture toughness of hydraulic fractures extending along natural fractures.

4. The method for calculating proppant migration and global fracture conductivity in horizontal well hydraulic fractures according to claim 1, characterized in that: The detailed steps of step S03 are: (1) Calculate the flow friction of the fracturing fluid in the horizontal wellbore and perforation holes: In the above formula, Δp w represents the horizontal wellbore friction; Δp p represents the friction resistance of the perforation hole; μ f represents the viscosity of the fracturing fluid; ρ f Indicates the density of fracturing fluid; Q w Indicates the flow rate in the horizontal wellbore; Q p represents the flow rate through the perforation hole; L represents the length of the fracturing section; C d Indicates the hole flow coefficient; n p Indicates the number of perforation holes; d p represents the perforation aperture; (2) Based on the pressure balance principle, the balance equation of the fluid pressure at the heel of the horizontal wellbore is established: p h,1 =p h,2 =…=p h,m In the above formula, p h Indicates the fluid pressure at the heel of the horizontal wellbore; Δp w Indicates the flow friction of fracturing fluid in the horizontal wellbore; Δp p It indicates the friction resistance of fracturing fluid flowing through the perforation holes; p in Indicates the fluid pressure at the entrance of the hydraulic fracture; (3) Based on the relationship between the perforation cluster flow rate and the fluid pressure balance equation, combined with the volume conservation and Newton iteration method, the flow rate of each perforation cluster entering the horizontal wellbore is calculated.

5. The method for calculating proppant migration and global fracture conductivity in horizontal well hydraulic fractures according to claim 1, characterized in that: The detailed steps of step S04 are: (1) Calculate the density and viscosity of the mixed sand solution formed by mixing the fracturing fluid and the proppant: r s =c p r p +(1-c p )r f In the above formula, ρ s represents the density of the sand mixing liquid; ρ p Indicates proppant density; μ s Indicates the viscosity of the sand mixing liquid; c p represents the proppant volume concentration; c p,max Indicates the maximum allowable volume concentration of proppant; represents the dimensionless proppant concentration; (2) Based on the Euler multiphase flow theory, the mass conservation equation of the flow of sand mixing fluid and proppant in the fracture is established: The flow rate calculation formula of sand mixing fluid and proppant in hydraulic fracture is: The proppant blocking function expression in hydraulic fracture is: The formula for calculating the sedimentation velocity of the proppant in the sand mixing fluid under the influence of gravity is: v set =v stokes ·C iner ·C conc ·C frac In the above formula, v s Indicates the flow velocity of the sand mixing liquid; v p Indicates the proppant flow velocity; v set Indicates proppant settling velocity; v stokes It represents the Stokes settling velocity of the proppant particles; represents the dimensionless sand concentration; B represents the proppant plugging function; d p Indicates the proppant particle size; C iner The correction factor that indicates the effect of particle inertia on sedimentation velocity; C conc Indicates the correction coefficient of the effect of proppant concentration on sedimentation velocity; C frac It represents the correction coefficient of the effect of crack wall on settlement velocity.

6. The method for calculating proppant migration and global fracture conductivity in horizontal well hydraulic fractures according to claim 1, characterized in that: The detailed steps of step S05 are: (1) Combining the relationship between displacement and stress of hydraulic fracture units and the mass conservation equation of fluid flow in hydraulic fractures, a set of nonlinear fluid-solid coupling equations for three-dimensional hydraulic fracture expansion is established: In the above formula, A M represents the boundary influence coefficient moment; I represents the unit matrix; C f represents the fluid flow coefficient matrix; σ M represents the comprehensive stress matrix acting on the hydraulic fracture unit; Qs represents the fluid flow source matrix; w represents the hydraulic fracture width matrix; p f represents the fluid pressure matrix in the hydraulic fracture; Before the proppant is pumped, the flow field in the above coupling equations adopts the mass conservation equation of the fracturing fluid. After the proppant is pumped, the flow field in the above coupling equations adopts the mass conservation equation of the sand mixing fluid. The nonlinear fluid-solid coupling equations are iteratively solved using Newton iteration to obtain the hydraulic fracture width and the fluid pressure in the fracture at the current moment. The iterative calculation format is: In the above formula, n represents the iteration step; Δw represents the iteration amount of hydraulic fracture width; Δp f represents the fluid pressure iteration amount; (2) After solving the nonlinear fluid-solid coupling equations to obtain the hydraulic fracture width and the fluid pressure in the fracture, the flow velocity of the sand mixing fluid is calculated first, and then the sedimentation velocity of the proppant is calculated. The migration velocity of the proppant in the hydraulic fracture during the fracturing process is obtained by superposition calculation: v p =B(w,d p )×(v f +v set ) (3) Substitute the proppant migration velocity into the proppant mass conservation equation to solve the distribution of proppant concentration in the fracture during the hydraulic fracture expansion process. First, the convection term and source term of the mass conservation equation are processed by the finite volume method: The QUICK upwind format is further used for numerical discretization: The right side expression of the proppant flow mass conservation equation can be rearranged to obtain: In the above formula, dA represents the area of ​​the control unit; Δx and Δy represent the length and width of the control unit respectively; subscripts n, s, w, and e represent the upper, lower, left, and right boundaries of the control unit respectively; subscript C represents the center of the control unit; subscripts N, S, W, and E represent the centers of the adjacent units on the upper, lower, left, and right sides of the control unit respectively; NN, SS, WW, and EE represent the centers of the adjacent units on the upper side of N, the lower side of S, the left side of W, and the right side of E respectively; max() represents the maximum value function; A p represents the proppant flow coefficient matrix; represents the dimensionless concentration matrix of proppant; Q p,s represents the proppant source term matrix; The fourth-order Runge-Kutta algorithm is introduced to discretize the unsteady-state term of the proppant flow mass conservation equation: In the above formula, t represents time; Δt represents the time step; F represents the right side term of the mass conservation equation; k1, k2, k3, k4 represent intermediate variables; Under the condition of obtaining the hydraulic fracture width and fluid pressure, the concentration distribution of the proppant in the fracture during the hydraulic fracture expansion process can be calculated through the high-order discrete format of the above proppant flow mass conservation equation in the one-way coupling framework.

7. The method for calculating proppant migration and global fracture conductivity in horizontal well hydraulic fractures according to claim 1, characterized in that: The detailed steps of step S06 are: (1) According to the distribution results of proppant concentration in the hydraulic fracture unit after the fracture is closed after pumping, the hydraulic fracture unit is divided into two categories: unfilled fracture unit and sand-filled fracture unit: In the above formula, c p represents the proppant volume concentration; represents the dimensionless proppant concentration; c p,max represents the maximum allowable volume fraction of proppant; c c represents the critical volume concentration; (2) For unfilled fracture units, the fracture walls will gradually close and contact after the pump is stopped during construction. The critical closure width and contact permeability considering the roughness of the fracture wall are introduced to calculate the conductivity: C1=w c ·k c In the above formula, C1 represents the conductivity of the unfilled fracture unit; w c represents the critical closing gap width; k c represents the contact permeability; (3) For the sand-filled fracture unit filled with proppant, first calculate the porosity: The permeability of the sand-filled fracture unit is calculated according to the Kozeny-Carman model: According to the stress-displacement relationship of the hydraulic fracture unit and the fracture width distribution of the unfilled fracture unit, the fracture width of the sand-filled fracture unit after the pump is stopped and the fracture is closed is calculated: Combined with the permeability and width of the sand-filled fracture unit, the conductivity of the sand-filled fracture unit is calculated: C2=in p ·k p In the above formula, φ p represents the porosity of the sand-filled fracture unit; k p represents the permeability of the sand-filled fracture unit; V s Table proppant equivalent thickness; V b represents the proppant filling thickness; w0 represents the width of the hydraulic fracture unit when the pump is stopped; w p k represents the width of the sand-filled fracture unit after the fracture is closed after pumping is stopped; K represents the Kozeny-Carman constant; S V Represents the surface area of ​​solid particles per unit volume; σ c represents the comprehensive stress acting on the crack unit; p f represents the fluid pressure in the hydraulic fracture; σ p represents the deformation stress of the proppant filling layer; σ f represents the contact stress between the proppant filling layer and the fracture wall; σ n represents the normal in-situ stress acting on the fracture wall; C2 represents the conductivity of the sand-filled fracture unit; (4) Combining the conductivity results of unfilled fracture units with those of filled fracture units, the global conductivity distribution of hydraulic fracturing fractures is obtained, which can provide a reference for the design and optimization of reservoir hydraulic fracturing construction.

Images