A method and system for predicting induced sand production and sand production area in deep water shallow gas reservoirs

Abstract

This invention relates to the field of rock mechanics engineering, and more particularly to a method and system for predicting sand production and sand-producing areas induced by depressurization mining in deep-water shallow gas reservoirs. The method includes: characterizing the physical properties of rock samples from weakly cemented deep-water reservoirs to obtain an initial set of pore structure parameters; conducting high-pressure triaxial shear softening experiments to extract strain softening constitutive parameters and establish dynamic mechanical parameter evolution equations; using thick-walled cylindrical depressurization displacement experiments to invert erosion parameters of the fluid-sand grain interaction behavior under the failure state described by the dynamic mechanical parameter evolution equations, obtaining the critical initiation velocity and erosion coefficient; establishing a mechanical yield criterion and a mass erosion rate equation, and coupling them with the porosity-permeability dynamic evolution equation to construct a multi-field coupled control equation set for fluid-solid sand; and numerically solving the multi-field coupled control equation set for fluid-solid sand to obtain the dynamic sand-producing area and cumulative sand production. This invention improves the accuracy of quantitative prediction of sand production.

Description

Technical Field

[0001] This invention relates to the field of rock mechanics engineering technology, and in particular to a method and system for predicting the amount of sand produced and the sand-producing area induced by depressurization mining in deep water shallow gas reservoirs. Background Technology

[0002] Deep-sea shallow gas reservoirs are rich in reserves and represent an important successor area for future offshore natural gas exploration and development. These reservoirs are generally characterized by large water depths, shallow burial depths, unconsolidated or weakly cemented formations, high porosity, high permeability, and low strength. During depressurization extraction, as pore fluids are continuously extracted, the effective formation stress increases dramatically, easily leading to instability in the mechanical state of the surrounding rocks and changes in the hydrodynamic field, resulting in severe sand production problems. Sand production not only wears downhole and subsea manifold equipment and clogs gas production devices, but in severe cases, it can also cause large-scale wellbore collapse and seafloor subsidence, resulting in significant economic losses and marine safety and environmental accidents. Therefore, accurately predicting the amount and area of ​​sand production during depressurization extraction in deep-water shallow gas reservoirs is a prerequisite for development strategies for this type of gas field.

[0003] Existing sand production prediction technologies typically employ shear or tensile failure models based on continuum mechanics, estimating the sand production range by determining whether the surrounding rock has entered the plastic yield zone. However, these methods only consider the rock mass mechanical yield failure caused by the increase in effective stress, neglecting the complete process of weakly cemented rock transforming into loose sand particles after yield failure, and subsequently undergoing dynamic erosion and migration under the drag force of high-speed fluids. This results in the predicted mechanical failure area generally being larger than the actual sand production area, and it is impossible to quantitatively calculate the cumulative sand production, making it difficult to meet the needs of safe depressurization mining and sand control completion engineering optimization in deep water and shallow gas formations. Summary of the Invention

[0004] This invention provides a method and system for predicting the amount of sand produced and the sand-producing area induced by depressurization mining in deep water shallow gas reservoirs, in order to overcome the deficiencies of existing technologies.

[0005] This invention provides a method for predicting the amount of sand produced and the sand-producing area induced by depressurization mining in deep-water shallow gas reservoirs, including:

[0006] S1. Physical property characterization of rock samples from deep-water weakly cemented reservoirs to obtain an initial set of pore structure parameters; S2. Based on the initial pore structure parameter set, conduct high-pressure triaxial shear softening experiments, extract strain softening constitutive parameters, and establish dynamic mechanical parameter evolution equations. S3. Using a thick-walled cylindrical pressure-reducing displacement experiment, the erosion parameters of the fluid-sand grain interaction behavior under the failure state described by the dynamic mechanical parameter evolution equation are inverted to obtain the critical start-up velocity and erosion coefficient. S4. Based on the dynamic mechanical parameter evolution equation and the critical start-up speed, establish the mechanical yield criterion and the mass erosion rate equation, and couple the porosity dynamic evolution equation to construct a multi-field coupled control equation set for fluidized solid sand. S5. Input the initial pore structure parameter set, the dynamic mechanical parameter evolution equation and the erosion coefficient into the finite volume method solver to numerically solve the multi-field coupled control equation set of the fluidized solid sand to obtain the dynamic sand discharge area and the cumulative sand discharge.

[0007] According to the method for predicting the amount of sand produced and the sand-producing area induced by depressurization mining in deep water shallow gas reservoirs provided by the present invention, step S1 further includes: S11. Using in-situ compaction data of seabed sediments, rock samples were compacted and consolidated under equal stress to obtain reconstructed core samples. S12. The pore distribution of the reconstructed core sample is scanned by nuclear magnetic resonance to obtain pore throat radius distribution data; S13. The three-dimensional spatial structure of the reconstructed core sample is imaged by CT scanning, and combined with the pore throat radius distribution data, an initial pore structure parameter set including initial porosity, initial permeability and mineral composition is obtained.

[0008] According to the method for predicting the amount of sand produced and the sand-producing area induced by depressurization mining in deep water shallow gas reservoirs provided by the present invention, step S2 further includes: S21. Under the conditions of simulated in-situ effective overlying rock pressure and pore pressure, consolidated drained and consolidated undrained triaxial shear tests were carried out on the reconstructed core samples, and the stress-strain curves were recorded to extract the initial elastic modulus, Poisson's ratio and peak strength. S22. For the strain softening stage after the peak of the stress-strain curve, the evolution relationship between cohesion and internal friction angle with respect to equivalent plastic strain is fitted using multiple sets of confining pressure data to establish a dynamic mechanical parameter evolution equation; the expression of the dynamic mechanical parameter evolution equation is:

[0009] in, For equivalent skew plastic strain, The dynamic cohesion varies with equivalent plastic strain. The initial peak cohesion, For residual cohesion, The critical plastic strain threshold for achieving the residual state.

[0010] According to the method for predicting the amount of sand produced and the sand-producing area induced by depressurization mining in deep water shallow gas reservoirs provided by the present invention, step S3 further includes: S31. Using a thick-walled cylindrical core with a central hole to simulate the wellbore shape, apply an equivalent radial confining pressure to the outer wall, reduce the internal central hole pressure at a set rate, and inject gas or liquid into it for displacement. Record the apparent flow rate of the fluid and the quality of the produced sand body in real time under each production pressure differential condition. S32. The erosion parameters of the produced sand body are fitted and inverted to obtain the critical starting velocity and erosion coefficient; the expression for the erosion coefficient is:

[0011] in, The erosion coefficient is... For time intervals, Time interval The total mass of sand actually produced was collected from the thick-walled cylindrical experiment. The effective volume of the rock that has undergone erosion. The density of the rock skeleton particles. Real-time porosity of the rock Let Darcy's apparent velocity vector be the fluid velocity vector. This is the critical start-up speed.

[0012] According to the method for predicting the amount of sand produced and the sand-producing area induced by depressurization mining in deep water shallow gas reservoirs provided by the present invention, step S4 further includes: S41. Based on the effective stress principle, the evolution equation of the dynamic mechanical parameters is substituted into the modified Mohr-Coulomb criterion to construct a mechanical yield criterion; the expression of the mechanical yield criterion is:

[0013] in, The yield function value, For the maximum principal effective stress, For the minimum principal effective stress, For equivalent skew plastic strain, The dynamic cohesion varies with equivalent plastic strain. The dynamic internal friction angle varies with equivalent plastic strain; S42. When the yield function value output by the mechanical yield criterion is greater than or equal to 0 and the local flow velocity is greater than the critical initiation velocity, calculate the mass erosion rate of rock per unit volume; the expression for the mass erosion rate is:

[0014] in, The mass erosion rate of rock per unit volume. The erosion coefficient is... The density of the rock skeleton particles. Real-time porosity of the rock Let Darcy's apparent velocity vector be the fluid velocity vector. The critical starting speed; S43. Couple the mass erosion rate to the dynamic evolution equations of porosity and permeability to construct a set of multi-field coupled control equations for fluid-solid sand; the expression of the set of multi-field coupled control equations for fluid-solid sand is as follows:

[0015]

[0016] in, Indicates time, Real-time permeability of the rock Initial penetration rate, This represents the initial porosity.

[0017] According to the method for predicting the amount of sand produced and the sand-producing area induced by depressurization mining in deep water shallow gas reservoirs provided by the present invention, step S5 further includes: S51. Establish a three-dimensional computational domain for the near-wellbore zone, densify the grid in the sand-sensitive area near the wellbore wall, and assign the initial porosity and initial permeability of the initial pore structure parameter set to a three-dimensional heterogeneous space according to the well logging interpretation results. S52. The fixed stress splitting method is used to solve the fluid-solid-sand multi-field coupled control equations in a sequential iterative manner. The solution order of the sequential iterative solution is as follows: implicitly solve the fluid pressure field, and substitute the pore pressure gradient as a volume force into the mechanics module to solve the stress field and equivalent plastic strain. S53. Determine the yield state based on the mechanical yield criterion, calculate the mass erosion rate, and explicitly update the porosity and permeability of multiple grid cells to obtain the single time step progression process. S54. Continuously optimize the single time step process until the set total mining time, extract the yield zone distribution, porosity field and velocity field at each moment, and obtain the prediction parameters; S55. Output the dynamic sand discharge area and cumulative sand discharge volume based on the predicted parameters.

[0018] According to the method for predicting the amount of sand produced and the sand-producing area induced by depressurization mining in deep water shallow gas reservoirs provided by the present invention, step S55 further includes: S551. Define the set of spatial grid cells that satisfy the mass erosion rate being greater than 0 and the real-time porosity being greater than the macroscopic collapse limit porosity of the formation as the dynamic sand production area. S552. The cumulative sand output is obtained by performing a full-time-space volume integral calculation on the mass erosion rate within the dynamic sand output area; the expression for the cumulative sand output is:

[0019] in, To calculate the cumulative sand output, The density of the rock skeleton particles. Indicates time, This is a dynamic sand-producing area. The mass erosion rate of rock per unit volume. Let be the volume of the mesh element.

[0020] The method for predicting sand production and sand-producing areas induced by depressurization mining in deep water shallow gas reservoirs according to the present invention further includes: S6. Based on the dynamic sand discharge area and the cumulative sand discharge volume, construct a production pressure difference capacity sand discharge risk response chart and determine the safe critical production pressure difference window; According to the method for predicting the amount of sand produced and the sand-producing area induced by depressurization mining in deep water shallow gas reservoirs provided by the present invention, step S6 further includes: S61. Input the logging physical property parameters and in-situ geostress parameters of the target gas well into the finite volume method solver, set multiple sets of different bottom hole flowing pressure boundary conditions, run the coupled model for each working condition, and extract the corresponding gas production curve and the response data of cumulative sand production as a function of pressure difference. S62. Using production pressure difference, gas production capacity, and sand production risk as three coordinate dimensions, the response data is plotted as a three-dimensional response chart. S63. Based on the maximum allowable sand holding capacity limit of the formation or the upper limit of the surface sand control equipment, determine the safe critical production pressure differential window by inversion from the three-dimensional response chart.

[0021] This invention also provides a prediction system for sand production and sand-producing areas induced by depressurization exploitation in deep-water shallow gas reservoirs, used to execute a prediction method for sand production and sand-producing areas induced by depressurization exploitation in deep-water shallow gas reservoirs as described in any of the above claims, comprising: Acquisition module: Used to acquire the initial pore structure parameter set for physical property characterization of rock samples from deep-water weakly cemented reservoirs; Extraction module: used for high-pressure triaxial shear softening experiments on the initial pore structure parameter set, extracting strain softening constitutive parameters, and establishing dynamic mechanical parameter evolution equations; Inversion module: used to perform erosion parameter inversion on the fluid-sand grain interaction behavior under the failure state described by the dynamic mechanical parameter evolution equation using a thick-walled cylindrical pressure reduction displacement experiment, to obtain the critical start-up velocity and erosion coefficient; Establishment Module: Based on the evolution equation of the dynamic mechanical parameters and the critical start-up speed, establish the mechanical yield criterion and the mass erosion rate equation, and couple the porosity and seepage dynamic evolution equation to construct a multi-field coupled control equation set for fluidized solid sand; The solution module is used to input the initial pore structure parameter set, the dynamic mechanical parameter evolution equation, and the erosion coefficient into the finite volume method solver to numerically solve the multi-field coupled control equation set of the fluid-solid sand, and obtain the dynamic sand discharge area and the cumulative sand discharge.

[0022] This invention provides a method and system for predicting sand production and sand-producing areas induced by depressurization mining in deep-water shallow gas reservoirs. By combining high-pressure triaxial shear softening experiments with thick-walled cylindrical depressurization displacement experiments, the method simultaneously obtains dynamic mechanical parameter evolution equations and key parameters such as erosion coefficient and critical initiation rate from an experimental perspective. This ensures that the constitutive equations and erosion equations relied upon for subsequent numerical simulations are supported by sufficient experimental data, fundamentally avoiding prediction biases caused by relying on empirical estimations for parameter values. This significantly improves the reliability of describing the mechanical response of deep-water weakly cemented reservoirs under in-situ high-pressure conditions. Secondly, by introducing a strain softening constitutive relation with equivalent plastic strain as an internal variable and coupling it with the modified Mohr-Coulomb yield criterion, the method can continuously track the complete evolution path of rocks from elastic deformation, peak failure to residual loose state, overcoming the limitations of traditional single plastic zone methods. The discriminant method has an inherent flaw: overestimating the sand-producing area due to neglecting the softening process. Further, by superimposing fluid-hydraulic erosion conditions on the mechanical yield criterion, the calculation of the mass erosion rate is triggered only when both yield failure and supercritical flow velocity conditions are simultaneously met. The erosion process is then fed back in real-time to the nonlinear dynamic evolution equations of porosity and permeability, achieving a strong coupling dynamic response between sand particle transport and reservoir property changes. This makes the predicted sand-producing area boundary and cumulative sand production more closely resemble the actual formation response. Finally, the sequential iterative solution strategy based on the finite volume method strictly guarantees the local conservation of mass and momentum during multi-field coupling. The constructed three-dimensional response chart of production pressure differential, production capacity, and sand production risk transforms complex numerical simulation results into intuitive engineering decision-making basis, providing quantitative support for determining the critical safe production pressure differential window and selecting the precision of sand control screens. Attached Figure Description

[0023] To more clearly illustrate the technical solutions in this invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are some embodiments of this invention. For those skilled in the art, other drawings can be obtained from these drawings without creative effort.

[0024] Figure 1 A schematic diagram of a method for predicting the amount of sand produced and the sand-producing area induced by depressurization mining in deep water shallow gas reservoirs, provided by the present invention; Figure 2A schematic diagram of a prediction system for predicting sand production and sand production area induced by depressurization mining in deep water shallow gas reservoirs, provided by the present invention; Figure 3 This is a schematic diagram illustrating the dynamic evolution of the radius of the sand-producing area around the well during depressurization mining of deep-water weakly cemented reservoirs as a function of mining time, provided in an embodiment of the present invention. Figure 4 This is a schematic diagram illustrating the dynamic evolution of the cumulative sand production volume of a single well and a single layer during the depressurization and exploitation process of a deep-water weakly cemented reservoir, as provided in an embodiment of the present invention. Detailed Implementation

[0025] To make the objectives, technical solutions, and advantages of this invention clearer, the technical solutions of this invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some, not all, embodiments of this invention, and should not be construed as limiting the invention. All other embodiments obtained by those skilled in the art based on the embodiments of this invention without creative effort are within the scope of protection of this invention. In the description of this invention, it should be understood that the terminology used is for descriptive purposes only and should not be construed as indicating or implying relative importance.

[0026] To better understand this invention, the research background of this invention will be explained in detail below.

[0027] While existing sand production prediction technologies and theories have been widely applied in conventional medium-deep consolidated sandstone reservoirs, they face a series of pressing problems when dealing with the special conditions of deep-water, shallow-gas, and weakly cemented reservoirs. Specifically: Traditional empirical prediction methods (such as well logging sonic transit time method, critical porosity method, etc.) mainly rely on historical statistical patterns of specific blocks, lack clear rock mechanics and seepage mechanics mechanisms, and cannot be accurately applied to deep-sea shallow gas exploration areas that lack long-term mining data.

[0028] Shear or tensile failure models based on continuum mechanics only consider rock mass yielding failure caused by an increase in effective stress. These methods neglect the dynamic erosion and migration process of weakly cemented rocks after yielding, which transforms them into loose sand grains and are subsequently dragged by high-speed fluids. This often leads to a predicted mechanical failure zone larger than the actual sand-producing zone, and makes it impossible to quantitatively calculate the cumulative sand production.

[0029] Existing fluid-structure interaction models mostly focus on macroscopic formation deformation analysis or single-phase and multiphase fluid seepage calculations, lacking a precise description of the strain constitutive relationship for deep-water shallow sediments. More importantly, existing models rarely consider the large-scale, nonlinear dynamic evolution of local reservoir porosity and permeability after rock skeleton destruction and massive sand migration and loss, failing to achieve true dynamic coupling modeling of sand migration.

[0030] In terms of experimental testing methods, conventional rock mechanics tests are mostly aimed at diagenetic and consolidated rock masses, making it difficult to accurately characterize the consolidation and drainage characteristics of shallow, weak sediments under in-situ high-pressure environments. Furthermore, traditional mechanical failure tests and sand production experiments lack effective cross-validation, resulting in a lack of systematic experimental calibration of the core parameters required for numerical simulations.

[0031] Some current technologies have attempted to study sand production problems in weakly cemented formations from specific perspectives (such as using a single continuous medium damage constitutive correction, a simplified local fluid drag force formula, or a simple sand control and displacement experiment only for low-pressure shallow sand). However, there is still a lack of a predictive method that can integrate physical simulation experiments, multi-field coupling algorithms, and rigorously perform wellbore-scale quantitative evaluation of sand production. This method cannot meet the needs of safe and efficient depressurization mining and sand control completion engineering optimization in deep water and shallow gas formations.

[0032] like Figure 1 As shown, this invention provides a method for predicting the amount of sand produced and the sand-producing area induced by depressurization mining in deep-water shallow gas reservoirs, including: S1. Physical properties of rock samples from deep-water weakly cemented reservoirs were characterized to obtain an initial set of pore structure parameters.

[0033] Step S1 further includes: S11. Using in-situ compaction data of seabed sediments, rock samples were prepared by isostress compaction and consolidation to obtain reconstructed core samples.

[0034] To address the characteristics of deep-water shallow gas layer rock samples being easily disturbed and poorly cemented, this invention first collects in-situ compaction data of seafloor sediments, including parameters such as in-situ effective overburden pressure, pore pressure, and sediment density. Then, based on the aforementioned in-situ compaction data, isostress compaction and consolidation treatment, consistent with the in-situ strata conditions, is applied to the undisturbed rock samples in the laboratory. Specifically, axial pressure and confining pressure are simultaneously applied to the rock samples in a triaxial pressure chamber, allowing the rock samples to complete compaction and consolidation under conditions equivalent to the in-situ stress state of the seafloor, ultimately yielding reconstructed core samples.

[0035] S12. The pore distribution of the reconstructed core sample is scanned by nuclear magnetic resonance to obtain pore throat radius distribution data.

[0036] In step S12, nuclear magnetic resonance (NMR) scanning is performed on the reconstructed core sample. The basic principle of NMR technology is that hydrogen protons in the rock pores generate relaxation signals under the excitation of an external magnetic field. The relaxation time is directly related to the pore size; that is, the larger the pore, the longer the relaxation time, and the smaller the pore, the shorter the relaxation time. Specifically, this invention processes the collected relaxation time distribution data through T2 spectrum inversion to convert it into pore throat radius distribution data, that is, the volume proportion distribution of pore throats of different sizes, thereby quantitatively describing the size composition of the pore channels inside the reconstructed core sample.

[0037] S13. The three-dimensional spatial structure of the reconstructed core sample is imaged by CT scanning, and combined with the pore throat radius distribution data, an initial pore structure parameter set including initial porosity, initial permeability and mineral composition is obtained.

[0038] After obtaining the pore throat radius distribution data, in step S13, the present invention further performs micron-level CT scanning on the same reconstructed core sample. The CT scan uses X-rays to perform layer-by-layer cross-sectional scanning of the core, generating a series of two-dimensional grayscale images. These two-dimensional cross-sectional images are then stacked and fused by a three-dimensional reconstruction algorithm to obtain a three-dimensional spatial geometric structure image of the pores and skeleton inside the core.

[0039] Subsequently, the present invention performs image segmentation processing on the three-dimensional reconstructed image, that is, distinguishes the porous phase and the solid framework phase in the image according to the gray level threshold, and calculates the ratio of the volume of the porous phase to the total volume to directly obtain the initial porosity; at the same time, the initial permeability is calculated by applying virtual seepage simulation to the three-dimensional pore network; in addition, mineral identification analysis is performed on the gray level characteristics of the solid framework phase, and the mineral composition is determined by combining energy dispersive spectroscopy or diffraction data.

[0040] Finally, the initial porosity, initial permeability, mineral composition and pore throat radius distribution data are integrated to form an initial pore structure parameter set, which serves as the input benchmark for subsequent mechanical experiments and numerical simulations.

[0041] S2. Based on the initial pore structure parameter set, a high-pressure triaxial shear softening experiment is conducted to extract strain softening constitutive parameters and establish dynamic mechanical parameter evolution equations.

[0042] Step S2 further includes: S21. Under the conditions of simulated in-situ effective overlying strata pressure and pore pressure, consolidated drained and consolidated undrained triaxial shear tests were carried out on the reconstructed core samples, and the stress-strain curves were recorded to extract the initial elastic modulus, Poisson's ratio and peak strength.

[0043] In step S21, after obtaining the initial pore structure parameter set, this invention uses the porosity, permeability, and mineral composition data therein to guide the setting of confining pressure conditions and the selection of drainage schemes for triaxial experiments. Specifically, in the triaxial pressure chamber, based on the in-situ effective overlying strata pressure and pore pressure of the target reservoir, consolidated drained triaxial shear experiments and consolidated undrained triaxial shear experiments are conducted on the reconstructed core samples. The difference between the two types of experiments lies in whether the drainage valve is open: in the consolidated drained experiment, pore water is freely discharged during the shearing process, and the pore pressure remains constant; in the consolidated undrained experiment, the drainage valve is closed, pore water cannot be discharged during the shearing process, and the pore pressure changes with shear deformation. Subsequently, at least three different confining pressure gradients are selected for both types of experiments, covering the in-situ effective stress range of the target reservoir. The axial stress and axial strain under each confining pressure condition are synchronously acquired throughout the entire process, and the complete stress-strain curve is recorded. Subsequently, the initial elastic modulus is extracted from the slope of the linear elastic segment of the stress-strain curve, the Poisson's ratio is extracted from the ratio of transverse strain to axial strain, and the peak strength is extracted from the highest point of the curve, forming a set of elastic mechanical parameters.

[0044] S22. For the strain softening stage after the peak of the stress-strain curve, the evolution relationship between cohesion and internal friction angle with respect to equivalent plastic strain is fitted using multiple sets of confining pressure data to establish a dynamic mechanical parameter evolution equation; the expression of the dynamic mechanical parameter evolution equation is:

[0045] in, For equivalent skew plastic strain, The dynamic cohesion varies with equivalent plastic strain. The initial peak cohesion, For residual cohesion, The critical plastic strain threshold for achieving the residual state.

[0046] In step S22, this invention extracts the post-peak data under each confining pressure condition for the data segment following the peak strength point in the stress-strain curve, i.e., the strain softening stage. For the post-peak data under each confining pressure condition, the portion of the axial strain exceeding the peak strain is converted into equivalent plastic strain, and the corresponding axial stress value at that plastic strain level is recorded. Then, based on the linear regression method of the Mohr-Coulomb strength criterion, the cohesion and internal friction angle values ​​at that plastic strain level are calculated by solving simultaneous equations for multiple sets of confining pressure data. The above calculation process is repeated for the entire plastic strain interval from the peak to the residual state to obtain a discrete data point sequence showing the variation of cohesion and internal friction angle with equivalent plastic strain.

[0047] Finally, piecewise linear fitting was performed on the above discrete data point sequence. Within the range from zero to the critical plastic strain threshold, the cohesion decreased linearly from the initial peak cohesion to the residual cohesion. After the plastic strain exceeded the critical plastic strain threshold, the cohesion remained unchanged at the residual cohesion. Thus, the evolution equation of dynamic mechanical parameters was established, and the cohesion and internal friction angle were expressed as explicit functions of the equivalent plastic strain.

[0048] S3. Using a thick-walled cylindrical pressure-reducing displacement experiment, the erosion parameters of the fluid-sand grain interaction behavior under the failure state described by the dynamic mechanical parameter evolution equation are inverted to obtain the critical start-up velocity and erosion coefficient.

[0049] Step S3 further includes: S31. Using a thick-walled cylindrical core with a central hole to simulate the wellbore shape, apply an equivalent radial confining pressure to the outer wall, reduce the internal central hole pressure at a set rate, and inject gas or liquid into it for displacement. Record the apparent flow rate of the fluid and the quality of the produced sand body in real time under each production pressure differential condition.

[0050] In step S31, after obtaining the evolution equation of dynamic mechanical parameters, the fluid-sand grain interaction behavior under mechanical failure state is further quantitatively characterized through a thick-walled cylindrical core depressurization displacement experiment. The experimental sample is a cylindrical rock core with a central hole, where the central hole simulates the wellbore channel, and the outer wall simulates the far-field formation boundary. The experimental setup uniformly applies radial confining pressure to the outer wall, the value of which is set based on the in-situ horizontal stress of the target reservoir, thus simulating the lateral constraint of the formation on the surrounding rock of the wellbore. Subsequently, a mixture of gas and formation water is injected into the central hole, while the fluid pressure inside the central hole is gradually reduced from its initial value at a set depressurization rate. Under the condition that the external confining pressure remains constant, a gradually increasing production pressure difference is formed between the central hole and the outer wall, simulating the depressurization process during trial production. Under each differential pressure condition, after the experiment stabilizes, the fluid flowing out of the central hole is collected at a certain time interval. The fluid volume during this time interval is recorded by a flow meter and divided by the time interval and the cross-sectional area of ​​the central hole to calculate the Darcy apparent velocity of the fluid. At the same time, the solid particles in the produced fluid are filtered, dried, and weighed to obtain the total mass of the produced sand body during this time interval.

[0051] S32. The erosion parameters of the produced sand body are fitted and inverted to obtain the critical starting velocity and erosion coefficient; the expression for the erosion coefficient is:

[0052] in, The erosion coefficient is... For time intervals, Time interval The total mass of sand actually produced was collected from the thick-walled cylindrical experiment. The effective volume of the rock that has undergone erosion. The density of the rock skeleton particles. Real-time porosity of the rock Let Darcy's apparent velocity vector be the fluid velocity vector. This is the critical start-up speed.

[0053] Furthermore, after obtaining the apparent flow velocity and mass of produced sand under different pressure differential conditions, for each pressure differential condition, based on the expression for the erosion coefficient, the candidate value of the erosion coefficient under the corresponding condition is calculated with the total mass of produced sand as the numerator and the product of the effective volume of eroded rock, the density of skeleton particles, the real-time solid volume fraction, the supercritical flow velocity and the time interval as the denominator.

[0054] For all data points under differential pressure conditions, this invention uses the least squares method to jointly fit the erosion coefficient and the critical start-up speed. The minimum sum of squared fitting residuals is used as the convergence criterion to determine the optimal values ​​of the erosion coefficient and the critical start-up speed that minimize the overall error of all data points, thus completing the erosion parameter inversion.

[0055] S4. Based on the dynamic mechanical parameter evolution equation and the critical start-up speed, establish the mechanical yield criterion and the mass erosion rate equation, and couple the porosity and seepage dynamic evolution equation to construct a multi-field coupled control equation set for fluidized solid sand.

[0056] Step S4 further includes: S41. Based on the effective stress principle, the evolution equation of the dynamic mechanical parameters is substituted into the modified Mohr-Coulomb criterion to construct a mechanical yield criterion; the expression of the mechanical yield criterion is:

[0057] in, The yield function value, For the maximum principal effective stress, For the minimum principal effective stress, For equivalent skew plastic strain, The dynamic cohesion varies with equivalent plastic strain. The dynamic internal friction angle varies with equivalent plastic strain.

[0058] In step S41, the effective stress principle refers to the fact that the actual stress borne by the formation rock skeleton is equal to the total stress minus the product of the pore fluid pressure and the Biot coefficient. That is, changes in pore pressure directly change the stress state of the rock skeleton. When constructing the mechanical yield criterion, this invention substitutes the dynamic cohesion and dynamic internal friction angle established in step S22, both of which are functions of equivalent plastic strain, into the modified Mohr-Coulomb criterion to obtain the mechanical yield criterion shown above. This criterion uses the maximum principal effective stress, the minimum principal effective stress, the dynamic cohesion, and the dynamic internal friction angle as inputs to calculate the yield function value: when the yield function value is less than zero, the mesh element is in an elastic state; when the yield function value is greater than or equal to zero, the mesh element is determined to have undergone shear yielding and entered a mechanical failure state.

[0059] S42. When the yield function value output by the mechanical yield criterion is greater than or equal to 0 and the local flow velocity is greater than the critical initiation velocity, calculate the mass erosion rate of rock per unit volume; the expression for the mass erosion rate is:

[0060] in, The mass erosion rate of rock per unit volume. The erosion coefficient is... The density of the rock skeleton particles. Real-time porosity of the rock Let Darcy's apparent velocity vector be the fluid velocity vector. This is the critical start-up speed.

[0061] In step S42, at the yield function value In the grid cell, it is further determined whether the fluid Darcy apparent velocity modulus in the cell exceeds the critical start-up velocity obtained by inversion in step S32. The fluid Darcy apparent velocity is the macroscopic seepage velocity calculated from the pressure gradient field according to Darcy's law, and the critical start-up velocity is the minimum velocity threshold that triggers the loose sand particles to begin migration.

[0062] Only when both of the above conditions are met simultaneously, this invention calculates the unit volume mass erosion rate of the grid cell by multiplying the erosion coefficient, the skeleton particle density, the real-time solid volume fraction, and the supercritical flow velocity (i.e., the difference between the actual flow velocity and the critical start-up velocity), thus obtaining the mass of sand body carried away and eroded by the fluid per unit volume and per unit time in each grid cell that meets the dual conditions within the current time step.

[0063] S43. Couple the mass erosion rate to the dynamic evolution equations of porosity and permeability to construct a set of multi-field coupled control equations for fluid-solid sand; the expression of the set of multi-field coupled control equations for fluid-solid sand is as follows:

[0064]

[0065] in, Indicates time, Real-time permeability of rocks, Initial penetration rate, This represents the initial porosity.

[0066] After obtaining the mass erosion rate, step S43 of this invention couples it to the porosity dynamic evolution equation. Specifically, the mass erosion rate is divided by the density of the skeleton particles to obtain the increase in porosity per unit time, i.e., the time partial derivative of porosity, thereby directly converting the mass loss of sand particles into an increase in pore volume.

[0067] After the porosity is updated, this invention calculates the real-time permeability using the Kozeny-Carman permeability model. This model corrects the initial permeability by multiplying the cube of the ratio of real-time porosity to initial porosity by the square of the ratio of initial solid volume fraction to real-time solid volume fraction, and outputs the real-time permeability for the current time step. At this point, both porosity and permeability are updated, forming a complete fluid-solid-sand multi-field coupled control equation set with the fluid pressure field and stress field.

[0068] S5. Input the initial pore structure parameter set, the dynamic mechanical parameter evolution equation and the erosion coefficient into the finite volume method solver to numerically solve the multi-field coupled control equation set of the fluidized solid sand to obtain the dynamic sand discharge area and the cumulative sand discharge.

[0069] Step S5 further includes: S51. Establish a three-dimensional computational domain for the near-wellbore zone, densify the grid in the sand-sensitive area near the wellbore wall, and assign the initial porosity and initial permeability of the initial pore structure parameter set to a three-dimensional heterogeneous space according to the well logging interpretation results.

[0070] In step S51, the present invention establishes a three-dimensional computational domain near the wellbore, centered on the target gas well. The domain size covers the wellbore to the far-field formation boundary. The mesh is refined in the sand-sensitive area near the wellbore wall. The mesh size in the refined area is set according to the actual wellbore size and the expected sand-affected radius (in one specific embodiment, the near-wellbore mesh accuracy is 2 mm). Subsequently, the initial porosity and initial permeability obtained in step S13 are assigned values ​​one by one in each mesh cell of the three-dimensional computational domain according to the well logging interpretation results, forming a three-dimensional initial physical property field reflecting the heterogeneity of the reservoir space, which serves as the initial input for subsequent solutions.

[0071] S52. The fixed stress splitting method is used to solve the fluid-solid-sand multi-field coupled control equations in a sequential iterative manner. The solution order of the sequential iterative solution is as follows: implicitly solve the fluid pressure field, and substitute the pore pressure gradient as a volume force into the mechanics module to solve the stress field and equivalent plastic strain.

[0072] In step S52, the fixed stress splitting method is a sequential iterative strategy for handling fluid-structure interaction nonlinear equations. Its core is to separate the solution of the fluid pressure field from the solution of the solid mechanical field within each time step. In each time step, this invention first implicitly discretizes the fluid mass conservation equation, using the porosity and permeability fields from the previous time step as coefficients, and the initial pore pressure at the outer boundary and the bottom-hole flowing pressure at the inner boundary as boundary conditions, to obtain the full-field fluid pressure distribution for the current time step. Then, the gradient of this pressure field is multiplied by the Biot coefficient and substituted into the solid momentum balance equation as an equivalent volume force to solve for the stress tensor of each grid element. Finally, the equivalent plastic strain field is iteratively calculated based on the elastoplastic constitutive relation.

[0073] S53. Determine the yield state based on the mechanical yield criterion, calculate the mass erosion rate, and explicitly update the porosity and permeability of multiple grid cells to obtain the single time step progression process.

[0074] In step S53, the present invention substitutes the equivalent plastic strain field and stress field output in step S52 into the mechanical yield criterion in step S41, and judges whether the yield function value is greater than or equal to zero for each grid cell; for the yield grid cell, the fluid pressure gradient field output in step S52 is further read, and the fluid Darcy apparent velocity of each grid cell is calculated according to Darcy's law and compared with the critical start-up velocity; for the grid cell that simultaneously satisfies the yield condition and the supercritical velocity condition, the erosion rate is calculated by substituting it into the mass erosion rate equation, and then the porosity and permeability of the grid cell are updated in an explicit format to complete the update of all field quantities in a single time step.

[0075] S54. Continuously optimize the single time step process until the set total mining time, extract the yield zone distribution, porosity field and velocity field at each moment, and obtain the prediction parameters.

[0076] Furthermore, in step S54, the updated porosity field and permeability field from step S53 are used as the initial input for the next time step, and the complete solution process from steps S52 to S53 is re-executed. This process continues step by step until the cumulative time reaches the set total mining time. At the end of each time step, the yield zone distribution, porosity field, and velocity field of the entire field at the current moment are extracted and stored to form a time series of prediction parameters.

[0077] S55. Output the dynamic sand discharge area and cumulative sand discharge volume based on the predicted parameters.

[0078] Step S55 further includes: S551. The set of spatial grid cells that satisfy the condition of mass erosion rate being greater than 0 and real-time porosity being greater than the macroscopic collapse limit porosity of the formation is defined as the dynamic sand production area.

[0079] For dynamic sand-producing regions, this invention performs dual-condition screening on the porosity field stored at each time step, grid by grid: first, the mass erosion rate of the grid cell is greater than 0; second, the real-time porosity of the grid cell exceeds the formation's macroscopic collapse limit porosity (this limit value is determined by rock mechanics experiments or well logging interpretation, characterizing the critical value when porosity increases to the point of overall structural instability and collapse). The set of grid cells that simultaneously meet the above two conditions is defined as a dynamic sand-producing region, whose spatial range dynamically expands as the time step progresses.

[0080] S552. The cumulative sand output is obtained by performing a full-time-space volume integral calculation on the mass erosion rate within the dynamic sand output area; the expression for the cumulative sand output is:

[0081] in, To calculate the cumulative sand output, The density of the rock skeleton particles. Indicates time, This is a dynamic sand-producing area. The mass erosion rate of rock per unit volume. Let be the volume of the mesh element.

[0082] In step S552, the present invention multiplies the mass erosion rate of all grid cells in the dynamic sand discharge area by the volume of each grid cell in the spatial domain and sums them up over all time steps in the time domain. Finally, it divides the result by the density of the skeleton particles and converts the mass integral result into a volume quantity to obtain the cumulative sand discharge amount. The calculation formula is as above.

[0083] This also includes: S6. Based on the dynamic sand discharge area and the cumulative sand discharge volume, construct a production pressure difference capacity sand discharge risk response chart and determine the safe critical production pressure difference window; Step S6 further includes: S61. Input the logging physical property parameters and in-situ geostress parameters of the target gas well into the finite volume method solver, set multiple sets of different bottom hole flowing pressure boundary conditions, run the coupled model for each working condition, and extract the corresponding gas production curve and the response data of cumulative sand production as a function of pressure difference.

[0084] In step S61, the present invention inputs the logging physical properties (including porosity, permeability, reservoir thickness, etc.) and in-situ stress parameters (including overburden pressure, maximum and minimum horizontal principal stresses, and initial pore pressure) of the target gas well into the finite volume method solver. The initial pore pressure of the far-field formation is set at the outer boundary of the model, and a series of different bottom-hole flowing pressures are sequentially set at the inner boundary of the wellbore. The difference between each bottom-hole flowing pressure and the initial pore pressure corresponds to a production pressure differential condition. For each pressure differential condition, the complete solution process from S52 to S55 is run, extracting the gas production curve (obtained by integrating the fluid pressure gradient and permeability at the inner boundary of the wellbore) and the cumulative sand production time curve under that condition, forming the response data for each condition.

[0085] S62. Using production pressure difference, gas production capacity, and sand production risk as three coordinate dimensions, the response data is plotted as a three-dimensional response chart.

[0086] In step S62, the present invention uses production pressure difference as the first coordinate axis, gas production capacity as the second coordinate axis, and sand production risk (characterized by cumulative sand production) as the third coordinate axis. The response data points extracted under each pressure difference condition in step S61 are plotted as a three-dimensional response chart. Each data point in the chart corresponds to a combination of production capacity and sand production under a pressure difference condition.

[0087] S63. Based on the maximum allowable sand holding capacity limit of the formation or the upper limit of the surface sand control equipment, determine the safe critical production pressure differential window by inversion from the three-dimensional response chart.

[0088] In step S63, the present invention uses the maximum allowable sand holding limit of the formation or the maximum sand output of the ground sand control equipment as the constraint boundary. The constraint plane is marked in the three-dimensional response chart. The intersection of the constraint plane with the response data surface corresponds to a set of critical production pressure difference values. The pressure difference value that maximizes the production capacity while the sand output does not exceed the upper limit of the constraint is read from it and determined as the safe critical production pressure difference window.

[0089] like Figure 2 As shown, the present invention also provides a prediction system for sand production and sand-producing areas induced by depressurization mining in deep-water shallow gas reservoirs, comprising: Acquisition Module 100: Used to acquire the initial pore structure parameter set for physical property characterization of rock samples from deep-water weakly cemented reservoirs; Extraction module 200: used for high-pressure triaxial shear softening experiments on the initial pore structure parameter set, extracting strain softening constitutive parameters, and establishing dynamic mechanical parameter evolution equations; Inversion module 300: used to perform erosion parameter inversion on the fluid-sand grain interaction behavior under the failure state described by the dynamic mechanical parameter evolution equation using a thick-walled cylindrical pressure reduction displacement experiment, to obtain the critical start-up velocity and erosion coefficient; Establishment Module 400: Based on the dynamic mechanical parameter evolution equation and the critical start-up speed, establish the mechanical yield criterion and the mass erosion rate equation, and couple the porosity and seepage dynamic evolution equation to construct a multi-field coupled control equation set for fluidized solid sand. Solver module 500: Used to input the initial pore structure parameter set, the dynamic mechanical parameter evolution equation and the erosion coefficient into the finite volume method solver, and numerically solve the multi-field coupled control equation set of the fluid-solid sand to obtain the dynamic sand discharge area and the cumulative sand discharge.

[0090] The device embodiments described above are merely illustrative. The units described as separate components may or may not be physically separate. The components shown as units may or may not be physical units; that is, they may be located in one place or distributed across multiple network units. Some or all of the modules can be selected to achieve the purpose of this embodiment according to actual needs. Those skilled in the art can understand and implement this without any creative effort.

[0091] Through the above description of the embodiments, those skilled in the art can clearly understand that each embodiment can be implemented by means of software plus necessary general-purpose hardware platforms, and of course, it can also be implemented by hardware. Based on this understanding, the above technical solutions, in essence or the part that contributes to the prior art, can be embodied in the form of a software product. This computer software product can be stored in a computer-readable storage medium, such as ROM / RAM, magnetic disk, optical disk, etc., including several instructions to cause a computer device (which may be a personal computer, server, or network device, etc.) to execute the method for predicting the amount of sand induced and the sand-producing area in deep-water shallow gas layer depressurization mining as described in various embodiments or some parts of embodiments.

[0092] The following section uses a deep-water shallow sedimentary layer as an example to introduce the experimental, modeling techniques, parameter analysis, and implementation methods of the present invention, which includes a method and system for predicting sand production and sand-producing areas induced by depressurization mining in deep-water shallow gas layers.

[0093] Step S101: Taking a test well as an example, predict sand production during pressure reduction and optimize sand control strategies. The formation of this gas field is in an incompletely consolidated state, belonging to a typical deep-water, weakly cemented, high-porosity, and high-permeability gas reservoir.

[0094] Step S102: Extract authentic core samples from the target stratigraphic layer, with the lithology mainly consisting of argillaceous siltstone. Physical property tests and microstructure scanning were conducted in the laboratory, determining the initial porosity of the reservoir to be 38%, the average initial permeability to be 800 mD, and the clay mineral content to be approximately 12%. The undisturbed rock samples were then re-stress-compacted and consolidated in the laboratory to prepare multiple sets of standard cylindrical test cores.

[0095] Step S103: Simulate the effective overlying strata pressure of 3.5 MPa in the target reservoir and conduct a consolidation-drained triaxial experiment. The experimental records show that the initial Young's modulus of the rock was 1.2 GPa and the Poisson's ratio was 0.35. The rock sample underwent mechanical shear yielding when the axial strain reached 2.5%, and the initial peak cohesion was measured to be 1.5 MPa, with an internal friction angle of 28 degrees.

[0096] Step S104: Completely record the data of the entire strain softening process after peak strength. As the equivalent plastic strain continues to increase, the rock sample exhibits a significant nonlinear strength decay characteristic. When the plastic strain reaches 5%, the rock skeleton is completely destroyed, the residual cohesion drops sharply to 0.1 MPa, and the residual internal friction angle drops to 22 degrees, completing the testing and calibration of the parameters for the entire stage of mechanical yielding and softening.

[0097] Step S105: A thick-walled cylindrical core with a central hole is used to simulate the microstructure of the surrounding rock in a real wellbore. A radial confining pressure of 5 MPa is uniformly applied externally to simulate in-situ horizontal in-situ stress. Formation water and gas are injected into the central hole, while the fluid pressure in the central hole is slowly reduced from 5 MPa at a rate of 0.5 MPa / d to simulate the depressurization process during trial production.

[0098] Step S106: Real-time monitoring of the sand content and apparent velocity of the produced fluid under different production pressure differentials. Experiments showed that when the pressure differential at the central orifice reached 1.5 MPa, the gas-water mixture began to produce a large amount of silt. Through inversion calculations derived from erosion theory, the critical initiation velocity for fluid-triggered loose sand particle migration was found to be 0.03 m / s, and the formation erosion coefficient was fitted to a value of 0.015.

[0099] Step S107: Extract the macroscopic physical properties of the formation and the key parameters of micromechanical softening and hydraulic erosion obtained from experiments. Establish a three-dimensional hydrodynamic solid-phase sand grain multi-field coupling prediction model in the near-wellbore zone based on the finite volume method. The model mesh is divided with an accuracy of 2 mm in the sand-producing sensitive area near the wellbore. The initial porosity field and permeability field are assigned three-dimensional spatial heterogeneous values ​​based on the well logging interpretation results.

[0100] Step S108: Set the outer boundary of the model to the initial pore pressure of the far-field formation at 16.5 MPa, and the inner boundary to the fluid production boundary at a constant bottom hole pressure. Set the initial production test pressure difference to 3 MPa, solve the 30-day transient pressure reduction production dynamic simulation, and alternately calculate and update the pressure field, stress field, and porosity field of the grid cells within each time step.

[0101] Step S109: On the third day of mining, shear yielding occurs in the direction of the minimum horizontal principal stress around the well, resulting in local loss of cohesion and transformation into loose, skeletonless sand body. The high-speed degassing airflow generates a strong hydraulic drag erosion effect in this damaged area, causing the local porosity to rapidly increase to over 50%.

[0102] Step S1010: The boundaries of the dynamic grid regions with porosity exceeding the formation collapse limit and a mass erosion rate greater than 0 are determined and updated in real time. After 30 days, the maximum radial intrusion depth of the well-perimeter erosion sand-producing area expands to 0.45 m. Through dual integration calculations over the entire spatial domain and the simulation time step, the predicted total volume of sand produced in the first month for a single well and single layer is 1.8 m³, forming a clear dominant sand-producing channel.

[0103] Step S1011: Continuously change the bottom-hole flowing pressure boundary conditions in the numerical model to simulate the dynamic response of reservoir sand production under different production pressure differential gradients from 1 MPa to 5 MPa. Extract the relative growth rate of gas production and the inflection point of the cumulative sand production surge under different pressure differential conditions, and draw a three-dimensional assessment and decision-making chart of production pressure differential, production capacity benefits, and sand production risk specific to the target shallow gas reservoir block.

[0104] Step S1012: Optimize the field test production parameters based on the actual operating conditions of the deepwater project. The forecast chart shows that the critical differential pressure window for maintaining trace sand production in this weakly cemented deepwater gas layer is 1.8 MPa. To balance the single-well test production capacity benefits with the safe service life of the downhole sand control string, the field test production engineering scheme design strictly limits the maximum production differential pressure in the initial stage of production to within 1.5 MPa, and explicitly recommends using 200 mm precision metal mesh sand control screens for sand control completion operations.

[0105] Figure 3 This is a curve showing the dynamic evolution of the radius of the sand-producing area around the well as the extraction time during the depressurization and extraction process of deep-water weakly cemented reservoirs using the method of this invention. Figure 3 The horizontal axis represents the extraction time in days (d), and the vertical axis represents the radius of the sand-producing area in meters (m). The curve shows that in the initial extraction phase, approximately 3 days, shear yielding first occurs in the direction of the minimum horizontal principal stress around the well, and the radius of the sand-producing area begins to increase from zero. Subsequently, with the continuous effect of the production pressure differential, the radius of the sand-producing area continues to expand outward, reaching its maximum expansion rate around 20 days of extraction. By the 30th day of extraction, the maximum radial penetration depth of the sand-producing area has expanded to 0.45 m. The above evolution pattern indicates that, under the set production pressure differential conditions, the sand-producing area around the well is not fully formed at the initial extraction time, but rather gradually expands radially from near the wellbore to the outer formation as the effective stress continues to redistribute and porosity accumulates.

[0106] Figure 4This is a dynamic evolution curve of the cumulative sand production volume of a single well and a single layer during the depressurization and exploitation process of deep-water weakly cemented reservoirs using the method of this invention, as a function of exploitation time. Figure 4 The horizontal axis represents mining time in days (d), and the vertical axis represents cumulative sand production volume in cubic meters (m³). The curve shows that the cumulative sand production volume increases slowly in the early stages of mining, corresponding to the sand-producing area being in its initial expansion phase with a small erosion area. As mining progresses to approximately 10 days, the growth rate of the cumulative sand production volume accelerates significantly, corresponding to a rapid expansion of the sand-producing area radius and more grid cells simultaneously meeting erosion conditions. By the 30th day of mining, the total cumulative sand volume produced in a single well and single layer reaches 1.8 m³. This curve is consistent with... Figure 3 The evolution curves of the sand-producing area radius correspond to each other on the time axis. The time period when the expansion rate of the sand-producing area radius accelerates is consistent with the time period when the slope of the cumulative sand-producing volume increases, reflecting a direct quantitative correlation between the spatial expansion of the sand-producing area and the growth of the cumulative sand-producing volume.

[0107] This invention constructs a quantitative prediction model for sand production that simultaneously considers both mechanical yielding and hydraulic erosion mechanisms. It overcomes the limitations of traditional single-mechanical plasticity zone discrimination methods, which tend to overestimate the sand production area, and empirical statistical methods, which lack physical mechanisms and are difficult to apply widely. It reconstructs the complete dynamic evolution process of deep-water weakly cemented gas layers, from matrix strain softening to sand formation and then to fluid drag migration. Simultaneously, through a combined experimental method of in-situ high-pressure triaxial shear softening and thick-walled cylindrical pressure-reducing displacement, it achieves experimental measurement of the mechanical and erosion parameters of deep-sea shallow soft sediments under in-situ high-pressure conditions, providing reliable support for nonlinear constitutive equations and fluid erosion equations. Furthermore, the multi-field coupled dynamic evolution simulation of fluid-solid sand particles based on the finite volume method strictly ensures the local conservation of mass and momentum during multiphase transport, overcoming the technical bottleneck of traditional continuous medium numerical simulations' inability to uniformly quantify the cumulative sand production volume and predict the dynamic sand production boundary region. The graphical production pressure differential, capacity response, and sand production risk optimization method transforms complex three-dimensional numerical simulation results into intuitive engineering application charts, providing visualized and quantitative decision support for setting critical safe production pressure differentials, formulating reasonable production allocation schemes, and selecting the accuracy of sand control completion tubing in deep-sea gas field development projects.

[0108] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, and not to limit them; although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some of the technical features; and these modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the spirit and scope of the technical solutions of the embodiments of the present invention.

Claims

1. A method for predicting sand production and sand-producing areas induced by depressurization mining in deep-water shallow gas reservoirs, characterized in that, include: S1. Physical property characterization of rock samples from deep-water weakly cemented reservoirs to obtain an initial set of pore structure parameters; S2. Based on the initial pore structure parameter set, conduct high-pressure triaxial shear softening experiments, extract strain softening constitutive parameters, and establish dynamic mechanical parameter evolution equations. S3. Using a thick-walled cylindrical pressure-reducing displacement experiment, the erosion parameters of the fluid-sand grain interaction behavior under the failure state described by the dynamic mechanical parameter evolution equation are inverted to obtain the critical start-up velocity and erosion coefficient. S4. Based on the dynamic mechanical parameter evolution equation and the critical start-up speed, establish the mechanical yield criterion and the mass erosion rate equation, and couple the porosity dynamic evolution equation to construct a multi-field coupled control equation set for fluidized solid sand. S5. Input the initial pore structure parameter set, the dynamic mechanical parameter evolution equation and the erosion coefficient into the finite volume method solver to numerically solve the multi-field coupled control equation set of the fluidized solid sand to obtain the dynamic sand discharge area and the cumulative sand discharge.

2. The method for predicting sand production and sand-producing areas induced by depressurization mining in deep-water shallow gas reservoirs according to claim 1, characterized in that, Step S1 further includes: S11. Using in-situ compaction data of seabed sediments, rock samples were compacted and consolidated under equal stress to obtain reconstructed core samples. S12. The pore distribution of the reconstructed core sample is scanned by nuclear magnetic resonance to obtain pore throat radius distribution data; S13. The three-dimensional spatial structure of the reconstructed core sample is imaged by CT scanning, and combined with the pore throat radius distribution data, an initial pore structure parameter set including initial porosity, initial permeability and mineral composition is obtained.

3. The method for predicting the amount of sand produced and the sand-producing area induced by depressurization mining in deep water shallow gas reservoirs according to claim 1, characterized in that, Step S2 further includes: S21. Under the conditions of simulated in-situ effective overlying rock pressure and pore pressure, consolidated drained and consolidated undrained triaxial shear tests were carried out on the reconstructed core samples, and the stress-strain curves were recorded to extract the initial elastic modulus, Poisson's ratio and peak strength. S22. For the strain softening stage after the peak of the stress-strain curve, the evolution relationship between cohesion and internal friction angle with respect to equivalent plastic strain is fitted using multiple sets of confining pressure data to establish a dynamic mechanical parameter evolution equation; the expression of the dynamic mechanical parameter evolution equation is: in, For equivalent skew plastic strain, The dynamic cohesion varies with equivalent plastic strain. The initial peak cohesion, For residual cohesion, The critical plastic strain threshold for achieving the residual state.

4. The method for predicting the amount of sand produced and the sand-producing area induced by depressurization mining in deep water shallow gas reservoirs according to claim 1, characterized in that, Step S3 further includes: S31. Using a thick-walled cylindrical core with a central hole to simulate the wellbore shape, apply an equivalent radial confining pressure to the outer wall, reduce the internal central hole pressure at a set rate, and inject gas or liquid into it for displacement. Record the apparent flow rate of the fluid and the quality of the produced sand body in real time under each production pressure differential condition. S32. The erosion parameters of the produced sand body are fitted and inverted to obtain the critical starting velocity and erosion coefficient; the expression for the erosion coefficient is: in, The erosion coefficient is... For time intervals, Time interval The total mass of sand actually produced was collected from the thick-walled cylindrical experiment. The effective volume of the rock that has undergone erosion. The density of the rock skeleton particles. Real-time porosity of the rock Let Darcy's apparent velocity vector be the fluid velocity vector. This is the critical start-up speed.

5. The method for predicting sand production and sand-producing areas induced by depressurization mining in deep-water shallow gas reservoirs according to claim 1, characterized in that, Step S4 further includes: S41. Based on the effective stress principle, the evolution equation of the dynamic mechanical parameters is substituted into the modified Mohr-Coulomb criterion to construct a mechanical yield criterion; the expression of the mechanical yield criterion is: in, The yield function value, For the maximum principal effective stress, For the minimum principal effective stress, For equivalent skew plastic strain, The dynamic cohesion varies with equivalent plastic strain. The dynamic internal friction angle varies with equivalent plastic strain; S42. When the yield function value output by the mechanical yield criterion is greater than or equal to 0 and the local flow velocity is greater than the critical initiation velocity, calculate the mass erosion rate of rock per unit volume; the expression for the mass erosion rate is: in, The mass erosion rate of rock per unit volume. The erosion coefficient is... The density of the rock skeleton particles. Real-time porosity of the rock Let Darcy's apparent velocity vector be the fluid velocity vector. The critical starting speed; S43. Couple the mass erosion rate to the dynamic evolution equations of porosity and permeability to construct a set of multi-field coupled control equations for fluid-solid sand; the expression of the set of multi-field coupled control equations for fluid-solid sand is as follows: in, Indicates time, Real-time permeability of the rock Initial penetration rate, This represents the initial porosity.

6. The method for predicting sand production and sand-producing areas induced by depressurization mining in deep water shallow gas reservoirs according to claim 1, characterized in that, Step S5 further includes: S51. Establish a three-dimensional computational domain for the near-wellbore zone, densify the grid in the sand-sensitive area near the wellbore wall, and assign the initial porosity and initial permeability of the initial pore structure parameter set to a three-dimensional heterogeneous space according to the well logging interpretation results. S52. The fixed stress splitting method is used to solve the fluid-solid-sand multi-field coupled control equations in a sequential iterative manner. The solution order of the sequential iterative solution is as follows: implicitly solve the fluid pressure field, and substitute the pore pressure gradient as a volume force into the mechanics module to solve the stress field and equivalent plastic strain. S53. Determine the yield state based on the mechanical yield criterion, calculate the mass erosion rate, and explicitly update the porosity and permeability of multiple grid cells to obtain the single time step progression process. S54. Continuously optimize the single time step process until the set total mining time, extract the yield zone distribution, porosity field and velocity field at each moment, and obtain the prediction parameters; S55. Output the dynamic sand discharge area and cumulative sand discharge volume based on the predicted parameters.

7. The method for predicting sand production and sand-producing areas induced by depressurization mining in deep-water shallow gas reservoirs according to claim 6, characterized in that, Step S55 further includes: S551. Define the set of spatial grid cells that satisfy the mass erosion rate being greater than 0 and the real-time porosity being greater than the macroscopic collapse limit porosity of the formation as the dynamic sand production area. S552. The cumulative sand output is obtained by performing a full-time-space volume integral calculation on the mass erosion rate within the dynamic sand output area; the expression for the cumulative sand output is: in, To calculate the cumulative sand output, The density of the rock skeleton particles. Indicates time, This is a dynamic sand-producing area. The mass erosion rate of rock per unit volume. Let be the volume of the mesh element.

8. The method for predicting the amount of sand produced and the sand-producing area induced by depressurization mining in deep water shallow gas reservoirs according to claim 1, characterized in that, Also includes: S6. Based on the dynamic sand discharge area and the cumulative sand discharge volume, construct a production pressure difference capacity sand discharge risk response chart and determine the safe critical production pressure difference window.

9. The method for predicting sand production and sand-producing areas induced by depressurization mining in deep water shallow gas reservoirs according to claim 8, characterized in that, Step S6 further includes: S61. Input the logging physical property parameters and in-situ geostress parameters of the target gas well into the finite volume method solver, set multiple sets of different bottom hole flowing pressure boundary conditions, run the coupled model for each working condition, and extract the corresponding gas production curve and the response data of cumulative sand production as a function of pressure difference. S62. Using production pressure difference, gas production capacity, and sand production risk as three coordinate dimensions, the response data is plotted as a three-dimensional response chart. S63. Based on the maximum allowable sand holding capacity limit of the formation or the upper limit of the surface sand control equipment, determine the safe critical production pressure differential window by inversion from the three-dimensional response chart.

10. A prediction system for sand production and sand-producing areas induced by depressurization exploitation in deep-water shallow gas reservoirs, used to execute the prediction method for sand production and sand-producing areas induced by depressurization exploitation in deep-water shallow gas reservoirs as described in any one of claims 1-9, characterized in that, include: Acquisition module: Used to acquire the initial pore structure parameter set for physical property characterization of rock samples from deep-water weakly cemented reservoirs; Extraction module: used for high-pressure triaxial shear softening experiments on the initial pore structure parameter set, extracting strain softening constitutive parameters, and establishing dynamic mechanical parameter evolution equations; Inversion module: used to perform erosion parameter inversion on the fluid-sand grain interaction behavior under the failure state described by the dynamic mechanical parameter evolution equation using a thick-walled cylindrical pressure reduction displacement experiment, to obtain the critical start-up velocity and erosion coefficient; Establishment Module: Based on the evolution equation of the dynamic mechanical parameters and the critical start-up speed, establish the mechanical yield criterion and the mass erosion rate equation, and couple the porosity and seepage dynamic evolution equation to construct a multi-field coupled control equation set for fluidized solid sand; The solution module is used to input the initial pore structure parameter set, the dynamic mechanical parameter evolution equation, and the erosion coefficient into the finite volume method solver to numerically solve the multi-field coupled control equation set of the fluid-solid sand, and obtain the dynamic sand discharge area and the cumulative sand discharge.

Images