Quantitative prediction method and device for sand production pressure difference and sand production amount of high-temperature and high-pressure sandstone reservoir

Abstract

The invention belongs to the technical field of oil and gas production and development, and discloses a sand production pressure difference and sand production quantity quantitative prediction method suitable for a high-temperature and high-pressure sandstone reservoir, which comprises the following steps: calculating static rock mechanical parameters of a target well; predicting initial crustal stress; establishing a porosity and permeability dynamic evolution model; establishing a three-dimensional finite element model; establishing a seepage stress coupling transient nonlinear finite element equation; applying well bottom pressure drop to simulate production, and determining sand production pressure difference; and establishing a sand production rate prediction model, and quantitatively predicting the sand production rate. According to the method, factors such as pore permeability stress, temperature change and water content influencing rock strength are comprehensively considered, a three-dimensional finite element model is established, a seepage stress coupling finite element equation and solving steps are established based on a reservoir elastic-plastic consolidation theory, critical plastic strain is introduced to quantitatively characterize the sand production risk, a sand production prediction model is established, and the rock strength is predicted. And finally, quantitative and accurate calculation and scientific pressure control of the critical sand production pressure difference and the sand production amount are realized.

Description

Technical Field

[0001] The invention belongs to the technical field of oil and gas production and development, and in particular relates to a method and a device for quantitatively predicting the sand production pressure difference and sand production amount applicable to high-temperature and high-pressure sandstone reservoirs. Background Art

[0002] High-temperature and high-pressure fractured sandstone oil and gas reservoirs are characterized by large structural stress and high risk of sand production. Sand production in the formation will lead to reservoir collapse, sand burial in the wellbore, and ground leakage. In the worst case, it will cause a well repair or even cause a blowout and other disastrous consequences. For example, the sand content of a certain well A in the trial production was as high as 6%, which was equivalent to 86.4m of sand production. 3 , resulting in the puncture of the downhole packer; during the second oil test of Well B, a large-scale trial production caused sand production in the reservoir. 484 days after it was put into production, sand production from the formation blocked the wellbore, and a second completion test and well repair were necessary to restore its production channel; during the oil test of Well C, signs of sand production appeared when the production pressure difference was 7MPa, and continuous tubing was used to flush sand to solve the problem of sand burial in the pipe column.

[0003] Therefore, accurate quantitative prediction of reservoir sand production risk and scientific formulation of pressure control system are the key to ensure safe oil and production test of high temperature and high pressure sandstone reservoirs. However, the existing sand production prediction models are mostly derived from production experience, with large errors and little reference value. It is urgent to establish a new method for accurately predicting the sand production of high temperature and high pressure sandstone reservoirs.

[0004] Patent application CN109826622A discloses a simulation system for simulating sand production from sandstone reservoirs, wherein a sand production device is arranged in a confining pressure loading device, and a rock sample is arranged in the sand production device; an axial pressure loading device is arranged above the sand production device and contacts the upper end of the sand production device, and a fluid providing device is connected to the sand production device; the sand production device is provided with a fluid outlet, and a sand collecting device is connected to the sand production device through the fluid outlet; the axial pressure loading device is used to apply axial pressure to the rock sample; the confining pressure loading device is used to apply confining pressure to the rock sample; the fluid providing device is used to introduce fluid into the sand production device when the confining pressure is greater than a preset threshold value, so as to carry the discrete sand generated by the rock sample to the sand collecting device through the fluid; the sand collecting device is used to calculate the total mass of the discrete sand, record the total time used by the sand production device to discharge the discrete sand, and determine the sand production rate based on the total mass and the total time. The sand production simulation system of this patent can simulate the sand production situation of sandstone reservoirs and provide a basis for the actual mining process; however, the simulation system of this patent simulates the sand production situation of sandstone reservoirs through indoor experiments, which is only applicable to sandstone reservoirs and does not simulate high-temperature and high-pressure sandstone reservoirs. At the same time, the equipment investment cost is high, and it consumes manpower and time. The simulation efficiency is low and the simulation results cannot support the actual construction site. Summary of the invention

[0005] The present invention aims to address the technical problems existing in the prior art and provide a method and device for quantitatively predicting the sand production pressure difference and sand production amount suitable for high-temperature and high-pressure sandstone reservoirs. The method comprehensively considers factors such as pore-permeability stress sensitivity, temperature changes, and the effect of water content on rock strength. A three-dimensional finite element model is established based on reservoir static rock mechanics, geostress, pore pressure, porosity, perforation and completion parameters. Based on the reservoir elastoplastic consolidation theory, a seepage stress coupling finite element equation and solution steps are established, and the critical plastic strain is introduced to quantitatively characterize the sand production risk, and a sand production prediction model is established. Finally, the critical sand production pressure difference and the sand production amount are quantitatively and accurately calculated, which provides an important theoretical basis for scientifically formulating a pressure control system and ensuring safe oil testing and trial production of high-temperature and high-pressure sandstone reservoirs.

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

[0007] A method for quantitatively predicting sand production pressure difference and sand production volume applicable to high-temperature and high-pressure sandstone reservoirs, the method comprising the following steps:

[0008] Step S1: Calculate the static rock mechanical parameters of the target well;

[0009] Step S2: predicting initial geostress based on static rock mechanics parameters, combined with image logging data and drilling reports;

[0010] Step S3: Combined with the pore-permeability stress sensitivity and temperature change, a porosity dynamic evolution model and a permeability dynamic evolution model are established to achieve full coupling calculation of fluid-solid-temperature three fields and predict reservoir physical property parameters;

[0011] Step S4: Establishing a three-dimensional finite element model based on the parameters obtained in steps S1-S3, combined with the wellbore size and perforation hole parameters;

[0012] Meshing the three-dimensional finite element model, assigning reservoir static rock mechanical parameters, reservoir physical property parameters, produced fluid density and viscosity parameters, and applying load and displacement boundary conditions;

[0013] Step S5: establishing a transient nonlinear finite element equation for seepage stress coupling;

[0014] Step S6: applying a bottom hole pressure drop to the three-dimensional finite element model in step S4 to simulate production. At the same time, using the seepage stress coupled transient nonlinear finite element equation in step S5, numerically iteratively solve the elastic-plastic deformation behavior of the model pores caused by the pressure drop, and determine the sand production pressure difference according to the critical plastic strain value of sand production and the distribution range of the plastic strain;

[0015] Step S7: Establish a sand production prediction model, and perform quantitative prediction of sand production based on the reservoir fluid seepage velocity and the critical plastic strain value of sand production.

[0016] Furthermore, the step S1 specifically includes:

[0017] Step S101: For the target reservoir, static rock mechanics parameters are experimentally obtained, and dynamic rock mechanics parameters are calculated based on logging data;

[0018] Step S102: obtaining a conversion coefficient between the static parameters and dynamic parameters of different rock mechanics parameters;

[0019] Step S103: Analyze the logging data of the target well, calculate the rock mechanics material characteristic parameters at various locations in the target well formation, and correct the rock mechanics material characteristic parameters with the help of conversion coefficients to obtain the static rock mechanics parameters of the target well.

[0020] Furthermore, in step S1, well logging data analysis is performed on the target well to calculate the rock mechanical material characteristic parameters at various locations in the target well formation, where:

[0021] The logging data at least include: Gama ray, resistivity, acoustic wave and density logging data; the acoustic wave includes compression wave and shear wave;

[0022] The rock mechanics material characteristic parameters include at least: elastic modulus, Poisson's ratio, shear strength and internal friction angle.

[0023] Further, in step S2, the initial geostress includes: the magnitude and direction of vertical geostress, maximum horizontal principal stress and minimum horizontal principal stress;

[0024] The calculation formula of the initial geostress is as follows:

[0025]

[0026] In the above formula, σ H is the maximum horizontal principal stress; σ h is the minimum horizontal principal stress; σ v is the vertical stress; α is the Biot coefficient; p w is the formation pore pressure; ν is the Poisson's ratio; β1 is the tectonic stress coefficient in the direction of the maximum horizontal principal stress; β2 is the tectonic stress coefficient in the direction of the minimum horizontal principal stress; TVD is the vertical depth of the reservoir; ρ is the reservoir rock density; g is the gravity acceleration vector.

[0027] Further, in step S3, the porosity dynamic evolution model is:

[0028]

[0029] In the above formula, E is the elastic modulus of sandstone; n is the porosity; n0 is the initial porosity; ν is the Poisson's ratio; γ is the thermal expansion coefficient; ε vis the product strain; dp w is the change value of pore pressure; (T-T0) is the change value of reservoir temperature;

[0030] Furthermore, in step S3, the permeability dynamic evolution model is:

[0031]

[0032] In the above formula, K is the permeability; K0 is the initial permeability; ε v is the product strain; n is the porosity; n0 is the initial porosity; E is the elastic modulus of sandstone; ν is the Poisson's ratio; γ is the thermal expansion coefficient; Δp w is the change in pore pressure; ΔT is the change in reservoir temperature.

[0033] Furthermore, in step S4, the applied load conditions include: gravity, initial geostress, initial pore pressure, and zero displacement constraint on the symmetry plane.

[0034] Furthermore, in step S4, the displacement boundary conditions applied include: the perforation surface pressure boundary is the bottom hole pressure, and the wellbore surface is the casing boundary.

[0035] Furthermore, the step S5 specifically includes:

[0036] The sandstone skeleton equilibrium equation, geometric equation, constitutive equation, fluid continuity equation and fluid seepage equation are combined to establish a seepage stress coupled transient nonlinear finite element equation with node displacement and pore pressure as basic unknowns.

[0037] Furthermore, the expression of the seepage stress coupled transient nonlinear finite element equation is:

[0038]

[0039] in,

[0040]

[0041] In the above formula, K is the coupling coefficient matrix of node pore pressure and rock skeleton bulk modulus; L is the coefficient matrix of elastic strain and displacement relationship; H is the coefficient matrix of plastic strain and displacement relationship; S is the coefficient matrix of stress and displacement differential relationship; a is the time increment; M is the coupling coefficient matrix of node pore pressure and shape function; N is the nonlinear finite element shape function; F is the external force vector; V is the seepage velocity vector; Δt is the time increment from the nth increment to the n+1th increment; Δ(u) n+1 is the displacement increment of the discrete unit node in the n+1th incremental step; Δ(p w ) n+1 is the pore pressure increment of the discrete unit node in the n+1th incremental step; (u)n is the node displacement of the discrete element at the nth incremental step; (p w ) n is the pore pressure of the discrete unit node at the nth incremental step; (u) n+1 is the node displacement of the discrete element at the n+1th incremental step; (p w ) n+1 is the pore pressure of the discrete element node in the n+1th incremental step.

[0042] Furthermore, the sand production prediction model established in step S7 is:

[0043] V 砂 =η(1-n)C1v w

[0044]

[0045] In the above formula, V 砂 is the sand production rate; η is the sand production coefficient; η1 is the empirical coefficient; n is the porosity; C1 is the density of sand and gravel; ν w is the reservoir fluid seepage velocity; ε pl is the equivalent plastic strain; is the critical plastic strain value for sand production.

[0046] Furthermore, the critical plastic strain value for sand production is 2.5%-3.5%.

[0047] Furthermore, the method further comprises the following steps:

[0048] Step S8: Establish a parameterized calculation module, construct a parameterized simulation method and parameter modification interface for predicting the sand production pressure difference and sand production volume, and calculate in real time the changes in pore pressure and plastic strain around the wellbore caused by the increase in bottom hole flow pressure.

[0049] In addition, the present invention also provides an electronic device suitable for quantitative prediction of sand production pressure difference and sand production volume of high-temperature and high-pressure sandstone reservoirs. The electronic device includes a processor, a memory, and a computer program stored in the memory and run on the processor. When the processor executes the program, it implements the method described in any one of the above items.

[0050] Compared with the prior art, the beneficial effects of the present invention are:

[0051] (1) The present invention provides a quantitative prediction method for the sand production pressure difference and sand production volume applicable to high-temperature and high-pressure sandstone reservoirs. The method comprehensively considers factors such as pore-permeability stress sensitivity, temperature change, and the effect of water content on rock strength. A three-dimensional finite element model is established based on reservoir static rock mechanics, geostress, pore pressure, porosity, perforation completion and other parameters. Based on the reservoir elastoplastic consolidation theory, a seepage stress coupling finite element equation and solution steps are established. The critical plastic strain is introduced to quantitatively characterize the sand production risk, and a sand production prediction model is established. Finally, the critical sand production pressure difference and sand production volume are quantitatively and accurately calculated, thereby scientifically guiding the formulation of a pressure control system, which is of great significance for ensuring safe oil testing and trial production of high-temperature and high-pressure sandstone reservoirs.

[0052] (2) The method of the present invention is used to quantitatively predict the sand production of high-temperature and high-pressure sandstone reservoirs, and pressure control is performed based on the quantitative prediction of sand production. The method has been applied to 5 wells in engineering practice, of which 4 wells did not produce sand under normal production pressure differential conditions, and 1 well produced sand under the predicted sand production pressure differential, with a compliance rate of 100%, verifying the effectiveness of the method of the present invention. At the same time, the method of the present invention was used to create a sand production pressure differential control chart considering the coupling effect of water content and formation pressure depletion, to guide oil testing and achieve precise control of pressure differential at oil and gas production sites. The method was successfully applied to the resumption of production of a well B in the target reservoir, with a maximum daily oil production of 92.8 tons and a controlled pressure differential of 48 MPa for safe production. BRIEF DESCRIPTION OF THE DRAWINGS

[0053] Figure 1 is a flow chart of a method according to an embodiment of the present invention;

[0054] Figure 2 is a rock mechanical parameter profile of a target well in an embodiment of the present invention;

[0055] Figure 3 A three-dimensional finite element mesh model of a target well according to an embodiment of the present invention;

[0056] Figure 4 A finite element mesh model for predicting sand production in a target well according to an embodiment of the present invention;

[0057] Figure 5 This is a flow chart of a numerical solution process for the sand pressure difference according to an embodiment of the present invention;

[0058] Figure 6 1 is a graph showing the pore elastic-plastic deformation results of a target well according to an embodiment of the present invention, wherein (a) is a step-by-step cloud graph of pore pressure, (b) is a step-by-step cloud graph of equivalent stress, and (c) is a step-by-step cloud graph of equivalent plastic strain;

[0059] Figure 7 A user subroutine interface of an embodiment of the present invention;

[0060] Figure 8is a sand production curve diagram of a target well in an embodiment of the present invention;

[0061] Fig. 9 This is a schematic diagram of the interface of the parameterized calculation module according to an embodiment of the present invention;

[0062] Fig.10 A pore pressure distribution cloud diagram obtained by simulating a target well in an embodiment of the present invention;

[0063] Fig.11 A cloud diagram of plastic strain distribution around a borehole obtained by simulating a target well in an embodiment of the present invention;

[0064] Note: Unclear coordinates or text in the figure does not affect the disclosure and understanding of the technical solution of the present invention. DETAILED DESCRIPTION

[0065] The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the 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.

[0066] Generally speaking, sand production in oil and gas wells is caused by high geostress and low rock strength. If the rock strength is strong enough or the geostress is small enough, there will be no sand production problem. When the geostress is large and the rock strength is small, one of the mechanical behavior characteristics of the formation is the plastic deformation of the rock, and many indoor tests and engineering experience have proved that sand production is closely related to the plastic deformation of the rock.

[0067] Example 1

[0068] The embodiment of the present invention provides a method for quantitatively predicting the sand production pressure difference and sand production amount suitable for high-temperature and high-pressure sandstone reservoirs. It comprehensively considers factors such as pore-seepage stress sensitivity, temperature changes, and the effect of water content on rock strength, and establishes a three-dimensional finite element model based on reservoir static rock mechanics, geostress, pore pressure, porosity, perforation and completion parameters. Based on the reservoir elastoplastic consolidation theory, the seepage stress coupling finite element equation and solution steps are established, and the critical plastic strain is introduced to quantitatively characterize the sand production risk, so as to achieve quantitative and accurate calculation of the critical sand production pressure difference and sand production amount.

[0069] Taking a well B in the target reservoir as an example, the method is described in detail.

[0070] Combination Figure 1 As shown, the method mainly comprises the following steps:

[0071] Step S1: Calculate the static rock mechanical parameters of the target well;

[0072] Step S1 specifically includes the following steps:

[0073] Step S101: For the target reservoir, static rock mechanics parameters are experimentally obtained, and dynamic rock mechanics parameters are calculated based on logging data;

[0074] Step S102: obtaining a conversion coefficient between the static parameters and dynamic parameters of different rock mechanics parameters;

[0075] Step S103: Analyze the logging data of the target well, calculate the rock mechanics material characteristic parameters at various locations in the target well formation, and correct the rock mechanics material characteristic parameters with the help of conversion coefficients to obtain the static rock mechanics parameters of the target well.

[0076] The static rock mechanics parameters obtained from the experiment are accurate and reliable, and can restore the mechanical properties of rocks under formation conditions, but the data points are limited; the dynamic rock mechanics parameters calculated from logging data are continuous, easy to obtain, economical and fast, but cannot be directly used for engineering calculations.

[0077] To this end, the embodiment of the present invention conducts core experiments to establish the conversion relationship between the dynamic and static rock mechanics parameters of the target reservoir, as shown in Table 1 below:

[0078] Table 1 Conversion relationship between dynamic and static rock mechanical parameters of target reservoir

[0079]

[0080]

[0081] Well logging data analysis is performed on the target well B to obtain the rock mechanical material characteristic parameters at various locations in the formation. The required logging data include: Gama ray, resistivity, acoustic wave (including compression wave and shear wave), density logging data, etc.; the obtained rock mechanical material characteristic parameters will include elastic modulus, Poisson's ratio, uniaxial compressive strength, cohesion, internal friction angle, etc.

[0082] The rock mechanics material characteristic parameters calculated at this time are dynamic parameters, which cannot be directly used for engineering calculations. They need to be corrected according to the conversion coefficients in Table 1 to obtain static rock mechanics parameters. Taking the elastic modulus as an example, the elastic modulus parameter profile before and after correction is as follows: Figure 2 shown.

[0083] In the embodiment of the present invention, the corrected static rock mechanics parameters are: elastic modulus of the perforated well section = 12 GPa, Poisson's ratio = 0.27, bonding strength = 10 MPa, and internal friction angle = 30 degrees.

[0084] Step S2: predicting initial geostress;

[0085] According to the static rock mechanics parameters obtained in step S1, combined with image logging data and drilling reports and other information, the magnitude and direction of the three components of initial geostress are obtained. The three components of initial geostress include vertical geostress, maximum horizontal principal stress and minimum horizontal principal stress. The calculation formula is as follows:

[0086]

[0087] In the above formula (1), σ H is the maximum horizontal principal stress; σ h is the minimum horizontal principal stress; σ v is the vertical stress; α is the Biot coefficient; p w is the formation pore pressure; ν is the Poisson's ratio; β1 is the tectonic stress coefficient in the direction of the maximum horizontal principal stress; β2 is the tectonic stress coefficient in the direction of the minimum horizontal principal stress; TVD is the vertical depth of the reservoir; ρ is the reservoir rock density; g is the gravity acceleration vector.

[0088] For the target well B, the parameter value of 5770 meters in depth was taken to obtain the magnitude and direction of the three components of the initial geostress, among which Sv=146.26MPa, Sh=141.89MPa, and SH=144.19MPa.

[0089] Step S3: Establishing a porosity and permeability dynamic evolution model;

[0090] Taking into account the influence of pore-permeability stress sensitivity and temperature change, the porosity dynamic evolution model and permeability dynamic evolution model are established to realize the full coupling calculation of fluid-solid-temperature three fields, improve the calculation accuracy, and realize accurate prediction of reservoir physical parameters, including:

[0091] (1) Porosity dynamic evolution model

[0092] Before optimization:

[0093]

[0094] The porosity is defined as:

[0095]

[0096] According to the effective stress principle, when fluid is extracted from the reservoir, the pore pressure of the sandstone will decrease, and the external stress of the sandstone, i.e., the confining pressure, will remain unchanged. Therefore, the effective stress of the sandstone will increase accordingly, compressing the rock skeleton, and the skeleton compression coefficient is:

[0097]

[0098] According to the A Bane and BA Maximov formulas, taking the differentials on both sides yields:

[0099]

[0100] Substituting formula (5) into formula (4), we get:

[0101]

[0102] In addition, the compression coefficient of the rock skeleton can be calculated using the following formula:

[0103]

[0104] Combined with the elastic modulus E of sandstone skeleton s Relationship with sandstone elastic modulus E and porosity n:

[0105]

[0106] Substituting equations (6) and (7) into equation (4), we can obtain the volume change of rock skeleton caused by the change of pore pressure:

[0107]

[0108] After considering the changes in pore pressure and temperature of the tight sandstone reservoir, the deformation of the tight sandstone skeleton is:

[0109]

[0110] Substituting equation (10) into equation (3), we get the optimized porosity dynamic equation:

[0111]

[0112] In the above formulas (2)-(11), n ​​is the porosity; n0 is the initial porosity; ε v is the product strain; V c V is the rock unit control volume; s is the volume of the dense sandstone skeleton; ΔV s is the volume deformation of the dense sandstone skeleton; c s is the compression coefficient of dense sandstone skeleton, MPa -1 ; σ s is the effective stress borne by the rock skeleton; dp w is the pore pressure change value, MPa; K g is the bulk modulus of the dense sandstone skeleton matrix; ν is the Poisson's ratio of the dense sandstone, dimensionless; E s is the elastic modulus of the dense sandstone skeleton, MPa; E is the elastic modulus of sandstone; ν is the Poisson's ratio; γ is the thermal expansion coefficient; (T-T0) is the reservoir temperature change value.

[0113] (2) Permeability dynamic evolution model

[0114] Before optimization:

[0115]

[0116] Based on the Kozeny-Carman permeability equation, the permeability ratio of the current state to the initial state is:

[0117]

[0118] According to the definitions of surface ratio, volume strain and porosity, we can get:

[0119]

[0120] in,

[0121]

[0122] Substituting the above formula (15) into formula (14), we get the optimized permeability dynamic equation:

[0123]

[0124] In the above formulas (12)-(16), K is permeability; K0 is initial permeability; n is porosity; n0 is initial porosity; ε v is the product strain; S v is the specific surface area; k z is a dimensionless parameter; S v0 is the initial specific surface area; β is the experimental coefficient, which is taken as 0.1; V v is the total pore volume of reservoir rock; ΔV v V is the change in pore volume of reservoir rock; c V is the rock unit control volume; s is the volume of the dense sandstone skeleton; E is the elastic modulus of sandstone; ν is the Poisson's ratio; γ is the thermal expansion coefficient; p w is the formation pore pressure; dp w is the change in pore pressure; ΔT is the change in reservoir temperature.

[0125] Based on Python language, Equation (11) and Equation (16) are edited into user subroutines to realize the full coupling calculation of fluid-solid-temperature three fields and improve the calculation accuracy. The edited user subroutines are as follows:

[0126]

[0127]

[0128] Step S4: Establishing a three-dimensional finite element model, specifically including:

[0129] According to the parameters obtained in steps S1-S3, combined with the wellbore size and perforation hole parameters, a three-dimensional finite element model is established, such as Figure 3 As shown, the three-dimensional finite element model includes a casing, a cement ring and a reservoir are arranged outside the casing in sequence, and perforations are arranged on the side wall of the casing.

[0130] The three-dimensional finite element model is meshed, and the reservoir static rock mechanical parameters (obtained in step S1), reservoir physical property parameters (obtained in step S3, including porosity and permeability), produced fluid density and viscosity parameters are assigned, and load and displacement boundary conditions are applied.

[0131] The load conditions include: gravity, initial geostress (obtained in step S2), initial pore pressure, and zero displacement constraint on the symmetry plane; the displacement boundary conditions include: the perforation surface pressure boundary is the bottom hole pressure, and the well wall surface is the casing boundary.

[0132] The finite element mesh model for sand production prediction of the target well according to the embodiment of the present invention is shown in Figure 4 shown.

[0133] Step S5: establishing a transient nonlinear finite element equation for seepage stress coupling;

[0134] Based on the effective stress principle, the total stress borne by the reservoir is jointly borne by the effective stress borne by the sandstone skeleton and the reservoir pore pressure; during the production process, the change in pore pressure will cause the change in effective stress, which will lead to changes in reservoir porosity and permeability. At the same time, these changes in turn affect the distribution of reservoir fluid seepage velocity and pore pressure, that is, there is an interaction between the strain generated in the reservoir and the seepage of the reservoir fluid. This interaction is the dynamic coupling of reservoir seepage and stress, also known as fluid-solid coupling.

[0135] The sandstone skeleton equilibrium equation, geometric equation, constitutive equation, fluid continuity equation and fluid seepage equation are combined to form a seepage stress coupled transient nonlinear finite element equation with node displacement and pore pressure as basic unknowns.

[0136] (1) Sandstone skeleton equilibrium equation and tensor form

[0137] The sandstone skeleton equilibrium equation is:

[0138]

[0139] In the above formula (17), σ x , σ y , σ z , τ xy , τ yx , τ xz , τ zx , τ zy, τ yz are the stress components of the six faces of the unit cell;

[0140] The tensor expression for the balance is:

[0141]

[0142] In the above formula (18), σ ij,j is the tensor form of the rock skeleton stress matrix; is the tensor form of the effective stress matrix of the rock skeleton; p w is the pore pressure; δ ij is the Planck coefficient; (δ ij p w ) ,j is the tensor form of the pore pressure matrix; f i is the external force of the unit.

[0143] (2) Geometric equations

[0144]

[0145] In the above formula (19), ε x , ε y , ε z , γ xy , γ yz , γ zx are the strain components of the six surfaces of the unit body; u, v, w are the displacement components in three directions.

[0146] (3) Reservoir rock constitutive equation

[0147]

[0148] In the above formula (20), {dσ} is the element stress differential matrix; D e is the elastic constitutive matrix of rock; D p is the plastic constitutive matrix; {dε} is the element strain differential matrix; F 准则 is the Mohr-Coulomb yield equation; E is the elastic modulus; ν is the Poisson's ratio; A is the hardening function, which is determined by the experimental stress-plastic strain curve of rock mechanics; T is the symbol for transposed matrix calculation; σ is stress; I1 is the first invariant of the stress tensor; I2 is the second invariant of the stress deviator; θ σ is the Lode stress angle; is the internal friction angle; is the gradient vector of the yield surface.

[0149] (4) Reservoir fluid continuity equation

[0150]

[0151] In the above formula (21), ρ w is the density of the wet fluid; v w is the reservoir fluid seepage velocity; x is the partial differential in the seepage direction; Among them, n w is the free fluid volume ratio of the reservoir; n t is the volume ratio of the bound fluid in the reservoir; J is the volume change ratio of the reservoir fluid; V w is the volume of free fluid in the pores of reservoir rock, V' is the total volume of reservoir fluid, V t is the bound fluid volume in the pores of reservoir rocks, V0 is the volume of reservoir fluid in the initial state;

[0152] (5) Reservoir fluid flow equation

[0153]

[0154] In the above formula (22), v w is the reservoir fluid seepage velocity (the velocity of the fluid relative to the solid skeleton); K is the permeability vector, in Darcy or millidarcy, 1D = 1000mD = 10 -3 μm 2 ;n w =ns, n is porosity, s is saturation. This paper assumes that the tight sandstone reservoir is full of saturated fluid, so n w =ns (with n in formula (21) w are two calculation methods); μ is the dynamic viscosity of the reservoir fluid; g is the gravity acceleration vector; p w is the pore pressure; is the gradient; x is the partial differential in the seepage direction.

[0155] Combining equations (17)-(22), the transient nonlinear finite element equation of seepage stress coupling is established:

[0156]

[0157] in:

[0158]

[0159] In the above formulas (23)-(24), K is the coupling coefficient matrix of node pore pressure and rock skeleton bulk modulus; L is the coefficient matrix of elastic strain and displacement relationship; H is the coefficient matrix of plastic strain and displacement relationship; S is the coefficient matrix of stress and displacement differential relationship; a is the time increment; M is the coupling coefficient matrix of node pore pressure and shape function; N is the nonlinear finite element shape function; F is the external force vector; V is the seepage velocity vector; Δt is the time increment from the nth increment to the n+1th increment; Δ(u) n+1 is the displacement increment of the discrete unit node in the n+1th incremental step; Δ(pw ) n+1 is the pore pressure increment of the discrete unit node in the n+1th incremental step; (u) n is the node displacement of the discrete element at the nth incremental step; (p w ) n is the pore pressure of the discrete unit node at the nth incremental step; (u) n+1 is the node displacement of the discrete element at the n+1th incremental step; (p w ) n+1 is the pore pressure of the discrete element node in the n+1th incremental step.

[0160] Step S6: applying bottom hole pressure drop to simulate production and determine the sand production pressure difference, specifically including:

[0161] According to the three-dimensional finite element model of step S4, the bottom hole pressure drop is applied to simulate production. The seepage stress coupled transient nonlinear finite element equation in step S5 is used to numerically iteratively solve the elastic-plastic deformation behavior of the model pore caused by the pressure drop, and the sand production pressure difference is determined based on the critical plastic strain value (3%) of sand production and the distribution range of plastic strain.

[0162] The seepage stress coupled transient nonlinear finite element equation (23) established in step S5 is numerically solved by using the iterative formula (24) in combination with the critical plastic strain value of sand production. The solution process is as follows: Figure 5 As shown, the specific steps include:

[0163] Reservoir production Finally, in the entire solution calculation process ending at the n+1th analysis step, the Δt time period is divided into n+1 analysis steps, and each analysis step can be divided into multiple incremental steps. The incremental size of this incremental step can be adjusted according to the calculation convergence situation.

[0164] Specifically, at the initial moment, the initial geological parameters of the reservoir are known. For rock and rock engineering problems, the initial geostress balance must be performed first, so as to obtain the reservoir state with initial geostress but no initial formation strain after hundreds of millions of years of geological deposition, that is, the unit node displacement (u) at the initial moment. 0 and the unit node pore pressure (p w ) 0 As a known condition, taking this time as the base point, the reservoir is first produced for a period of time Δt1, and the analysis step calculation is performed to form the finite element equation coefficient matrix at the i=0 analysis step, that is, the initial moment. If it does not converge at this time, the incremental step size of the analysis step will be reduced, and the coefficient matrix will be adjusted to obtain better convergence; if the calculation converges, the unit node displacement increment Δ(u) at the i=1 analysis step will be output 1 and the increment of pore pressure at the unit node Δ(p w ) 1, and thus the node displacement (u) is obtained 1 , node pore pressure (p w ) 1 The corresponding reservoir strain distribution and stress distribution are further calculated by the reservoir rock geometry equation and the constitutive equation. The redistribution of reservoir stress, strain and formation pressure will affect the changes in reservoir porosity, permeability and production dynamics. This change is controlled by the dynamic change model of reservoir physical property parameters derived in the previous chapter. The subroutine is secondary developed based on this model on the ABAQUS finite element platform. The data calculated by the subroutine will be passed to the main program, thus forming the coefficient matrix of the fluid-solid coupling finite element equation for the next analysis step.

[0165] Similarly, by repeatedly iterating and solving, we can find Δ(u) in the n+1th analysis step. n+1 and Δ(p w ) n+1 , thus obtaining (u) at the n+1th analysis step n+1 and (p w ) n+1 , the strain is further calculated by the geometric equation, and the calculated stress and pore pressure distribution are obtained through the constitutive equation and the effective stress principle, until the plastic strain reaches the critical plastic strain value of sand production (plastic strain ≥ 3%), the whole calculation process is completed, and the corresponding sand production pressure difference is obtained.

[0166] The pore elastic-plastic deformation result of the target well finally calculated by the embodiment of the present invention is as follows: Figure 6 As shown, Figure (a) is the step-by-step cloud diagram of pore pressure, Figure (b) is the step-by-step cloud diagram of equivalent stress, and Figure (c) is the step-by-step cloud diagram of equivalent plastic strain.

[0167] Step S7: Establishing a sand production prediction model;

[0168] Considering the fluid drag force, a sand production prediction model is established according to the reservoir fluid seepage velocity and the critical plastic strain value of sand production:

[0169] V 砂 =η(1-n)C1v w (25)

[0170]

[0171] In the above formulas (25)-(26), V 砂 is the sand production rate; η is the sand production coefficient; η1 is the empirical coefficient, which is taken as 4 here; n is the porosity; C1 is the density of sand and gravel; ν w is the reservoir fluid seepage velocity; ε pl is the equivalent plastic strain; is the critical plastic strain value for sand production, that is, the maximum plastic strain that can be tolerated. Here, it is taken as 3%. This value is the sand production judgment value obtained through a large number of indoor experimental studies, numerical simulations and on-site evaluation verification.

[0172] Based on Python language, edit formula (25) and formula (26) as follows Figure 7 The user subroutine shown can realize the quantitative prediction of sand production.

[0173] The sand production curve finally calculated by the embodiment of the present invention is as follows: Figure 8 shown.

[0174] After the simulation and calculation of the above steps, the results show that when the production pressure difference of the target well B in the secondary oil test stage is 49.48MPa, the plastic strain reaches the critical value of 3%, and there is a risk of sand production; when the pressure difference is 59.5MPa, sand production begins, and the sand production per unit area of ​​the hole is 0.14kg / s, and the converted sand production rate of the whole wellbore is <2.05L / min; the actual sand production pressure difference is 61.6MPa, and the error is only 3.4%.

[0175] The present invention was applied to a total of 5 wells in the target reservoir, of which 4 wells did not produce sand under normal production pressure differential conditions, and 1 well produced sand after the pressure differential for sand production was higher than the predicted one, with a compliance rate of 100%, verifying the effectiveness of the method of the present invention; at the same time, the method of the present invention was applied to create a sand production pressure differential control chart considering the coupling effect of water content and formation pressure depletion (specifically, by changing the water content and formation pressure depletion, a large number of calculations were carried out, and sand production pressure differentials under different water contents and different formation pressures were formed, and a cloud map was drawn), to guide oil testing and achieve precise control of pressure differentials in oil and gas production sites, and it was successfully applied to the resumption of production of a well B in the target reservoir, with a maximum daily oil production of 92.8 tons and a controlled pressure differential of 48 MPa for safe production.

[0176] Example 2

[0177] Based on the method of Example 1, the embodiment of the present invention further adds the following steps:

[0178] Step S8: Establishing a parameterized calculation module;

[0179] In order to reduce modeling time and improve prediction efficiency, a Python language program was compiled (see Figure 7 ), a parametric simulation method for sand pressure difference and sand volume prediction and a parameter modification interface were established (see Fig. 9 ), which can calculate the changes in pore pressure and plastic strain around the wellbore caused by the increase in bottom hole flow pressure in real time (see Figure 10-11 ), computing efficiency increased by 80%.

[0180] Example 3

[0181] An embodiment of the present invention provides an electronic device suitable for quantitatively predicting the sand production pressure difference and sand production amount of high-temperature and high-pressure sandstone reservoirs. The electronic device includes a processor, a memory, and a computer program stored in the memory and executable on the processor. When the processor executes the program, it implements the method for quantitatively predicting the sand production pressure difference and sand production amount suitable for high-temperature and high-pressure sandstone reservoirs as described in Example 1 or Example 2.

[0182] The above description is only an embodiment of the present invention and is not intended to limit the present invention. Any modification, equivalent replacement and improvement made within the application scope of the present invention shall be included in the protection scope of the present invention.

Claims

1. A quantitative prediction method for sand production pressure difference and sand production volume applicable to high temperature and high pressure sandstone reservoirs, characterized in that: The method comprises the following steps: Step S1: Calculate the static rock mechanical parameters of the target well; Step S2: predicting initial geostress based on static rock mechanics parameters, combined with image logging data and drilling reports; Step S3: Combined with the pore-permeability stress sensitivity and temperature change, a porosity dynamic evolution model and a permeability dynamic evolution model are established to achieve full coupling calculation of fluid-solid-temperature three fields and predict reservoir physical property parameters; Step S4: Establishing a three-dimensional finite element model based on the parameters obtained in steps S1-S3, combined with the wellbore size and perforation hole parameters; Meshing the three-dimensional finite element model, assigning reservoir static rock mechanical parameters, reservoir physical property parameters, produced fluid density and viscosity parameters, and applying load and displacement boundary conditions; Step S5: establishing a transient nonlinear finite element equation for seepage stress coupling; Step S6: applying a bottom hole pressure drop to the three-dimensional finite element model in step S4 to simulate production. At the same time, using the seepage stress coupled transient nonlinear finite element equation in step S5, numerically iteratively solve the elastic-plastic deformation behavior of the model pores caused by the pressure drop, and determine the sand production pressure difference according to the critical plastic strain value of sand production and the distribution range of the plastic strain; Step S7: Establish a sand production prediction model, and perform quantitative prediction of sand production based on the reservoir fluid seepage velocity and the critical plastic strain value of sand production.

2. The method according to claim 1, characterized in that The step S1 specifically includes: Step S101: For the target reservoir, static rock mechanics parameters are experimentally obtained, and dynamic rock mechanics parameters are calculated based on logging data; Step S102: obtaining a conversion coefficient between the static parameters and dynamic parameters of different rock mechanics parameters; Step S103: Analyze the logging data of the target well, calculate the rock mechanics material characteristic parameters at various locations in the target well formation, and correct the rock mechanics material characteristic parameters with the help of conversion coefficients to obtain the static rock mechanics parameters of the target well.

3. The method according to claim 2, characterized in that In step S1, well logging data analysis is performed on the target well to calculate the rock mechanical material characteristic parameters at various locations in the target well formation, where: The logging data at least include: Gama ray, resistivity, acoustic wave and density logging data; the acoustic wave includes compression wave and shear wave; The rock mechanics material characteristic parameters include at least: elastic modulus, Poisson's ratio, shear strength and internal friction angle.

4. The method according to claim 1, characterized in that: In step S2, the initial geostress includes: the magnitude and direction of vertical geostress, maximum horizontal principal stress and minimum horizontal principal stress; The calculation formula of the initial geostress is as follows: In the above formula, σ H is the maximum horizontal principal stress; σ h is the minimum horizontal principal stress; σ v is the vertical stress; α is the Biot coefficient; p w is the formation pore pressure; ν is the Poisson's ratio; β1 is the tectonic stress coefficient in the direction of the maximum horizontal principal stress; β2 is the tectonic stress coefficient in the direction of the minimum horizontal principal stress; TVD is the vertical depth of the reservoir; ρ is the reservoir rock density; g is the gravity acceleration vector.

5. The method according to claim 1, characterized in that In step S3, the porosity dynamic evolution model is: In the above formula, E is the elastic modulus of sandstone; n is the porosity; n0 is the initial porosity; ν is the Poisson's ratio; γ is the thermal expansion coefficient; ε v is the product strain; dp w is the change of pore pressure; (T-T0) is the change of reservoir temperature.

6. The method according to claim 5, characterized in that In step S3, the permeability dynamic evolution model is: In the above formula, K is the permeability; K0 is the initial permeability; ε v is the product strain; n is the porosity; n0 is the initial porosity; E is the elastic modulus of sandstone; ν is the Poisson's ratio; γ is the thermal expansion coefficient; Δp w is the change in pore pressure; ΔT is the change in reservoir temperature.

7. The method according to claim 1, characterized in that In step S4, the applied load conditions include: gravity, initial geostress, initial pore pressure, and zero displacement constraint on the symmetry plane.

8. The method according to claim 1, characterized in that In step S4, the displacement boundary conditions applied include: the perforation surface pressure boundary is the bottom hole pressure, and the well wall surface is the casing.

9. The method according to claim 1, characterized in that: The step S5 specifically includes: The sandstone skeleton equilibrium equation, geometric equation, constitutive equation, fluid continuity equation and fluid seepage equation are combined to establish a seepage stress coupled transient nonlinear finite element equation with node displacement and pore pressure as basic unknowns.

10. The method according to claim 9, characterized in that The expression of the seepage stress coupled transient nonlinear finite element equation is: in, In the above formula, K is the coupling coefficient matrix of node pore pressure and rock skeleton bulk modulus; L is the coefficient matrix of elastic strain and displacement relationship; H is the coefficient matrix of plastic strain and displacement relationship; S is the coefficient matrix of stress and displacement differential relationship; a is the time increment; M is the coupling coefficient matrix of node pore pressure and shape function; N is the nonlinear finite element shape function; F is the external force vector; V is the seepage velocity vector; Δt is the time increment from the nth increment to the n+1th increment; Δ(u) n+1 is the displacement increment of the discrete unit node in the n+1th incremental step; Δ(p w ) n+1 is the pore pressure increment of the discrete unit node in the n+1th incremental step; (u) n is the node displacement of the discrete element at the nth incremental step; (p w ) n is the pore pressure of the discrete unit node at the nth incremental step; (u) n+1 is the node displacement of the discrete element at the n+1th incremental step; (p w ) n+1 is the pore pressure of the discrete element node in the n+1th incremental step.

11. The method according to claim 1, characterized in that: The sand production prediction model established in step S7 is: V 砂 =η(1-n)C1v w In the above formula, V 砂 is the sand production rate; η is the sand production coefficient; η1 is the empirical coefficient; n is the porosity; C1 is the density of sand and gravel; ν w is the reservoir fluid seepage velocity; ε pl is the equivalent plastic strain; is the critical plastic strain value for sand production.

12. The method according to claim 1, characterized in that The critical plastic strain value for sand production is 2.5%-3.5%.

13. The method according to any one of claims 1 to 12, characterized in that: The method further comprises the steps of: Step S8: Establish a parameterized calculation module, construct a parameterized simulation method and parameter modification interface for predicting the sand production pressure difference and sand production volume, and calculate in real time the changes in pore pressure and plastic strain around the wellbore caused by the increase in bottom hole flow pressure.

14. An electronic device suitable for quantitatively predicting the sand production pressure difference and sand production volume of a high-temperature and high-pressure sandstone reservoir, the electronic device comprising a processor, a memory, and a computer program stored in the memory and executable on the processor, characterized in that: When the processor executes the program, the method according to any one of claims 1 to 13 is implemented.

Images