Fracture-vug type reservoir well wall stability evaluation method based on damage heat flow solidification coupling

Abstract

The invention belongs to the technical field of well drilling, and relates to a fracture-vug type reservoir well wall stability evaluation method based on damage heat flow solidification coupling, which comprises the following steps: constructing a discrete fracture-vug network model based on fracture-vug geometric parameters; constructing a fracture-vug type reservoir heat flow solidification coupling mathematical model in which the damage variables are introduced; based on the discrete fracture-vug network model and the heat flow solidification coupling mathematical model, performing fracture-vug type reservoir heat flow solidification coupling numerical simulation; based on the numerical simulation result and the heat flow solidification coupling mathematical model, a damage criterion is introduced to evaluate the stability of the well wall; according to the method, by integrating discrete fracture-cavity network modeling, heat flow solidification multi-field coupling simulation introducing damage variables and dynamic stability evaluation, a technical scheme closer to actual geological conditions and physical processes is provided for well wall stability evaluation in the fracture-cavity type reservoir drilling and completion process, and the accuracy of evaluation results is effectively improved.

Description

Technical Field

[0001] This invention relates to the field of drilling technology, and more specifically, to a method for evaluating wellbore stability in fractured-vuggy reservoirs based on damage-thermal-solidification coupling. Background Technology

[0002] Currently, numerical simulation technology has become a key tool for studying multiphysics processes in reservoirs and evaluating wellbore stability. Among existing technologies, numerical models for fluid-structure interaction and thermo-fluid-structure interaction are relatively mature and can describe the basic physical processes of fractured-vuggy reservoirs to a certain extent. However, existing models have significant shortcomings:

[0003] First, regarding model completeness, most existing models fail to systematically integrate chemical solute transport processes and damage evolution effects caused by progressive rock failure. The absence of a chemical field ignores the potential impact on rock mechanical properties; the lack of damage evolution means the models cannot describe the time-varying behavior of material properties as deformation deteriorates. Second, existing instability criteria are overly simplistic. Current wellbore stability evaluation systems typically rely on classical strength theories such as the Mohr-Coulomb criterion and tensile criteria. These criteria are essentially static, based on critical states, and can only identify the instantaneous threshold of failure, failing to characterize the dynamic, progressive process of wellbore instability from microcrack initiation and propagation to macroscopic penetration. During this process, the evolution of key parameters such as permeability and elastic modulus caused by damage accumulation exacerbates the development of failure, a damage mechanism that existing methods cannot reflect.

[0004] Therefore, existing evaluation methods suffer from insufficient accuracy and lagging predictive ability in assessing wellbore stability of fractured-vuggy reservoirs due to incomplete models and overly simplified instability criteria. Summary of the Invention

[0005] This invention provides a wellbore stability evaluation method for fractured-vuggy reservoirs based on damage-thermal-solidification coupling. Its purpose is to solve the problems of insufficient evaluation accuracy and lagging prediction ability caused by incomplete models and oversimplified instability criteria in the existing technology.

[0006] The technical solution of the present invention is as follows:

[0007] A wellbore stability evaluation method for fractured-vuggy reservoirs based on damage-thermal-solidification coupling is characterized by the following steps:

[0008] S1. Construct a discrete slot network model based on the geometric parameters of the slots;

[0009] S2. Construct a mathematical model of thermal flux-solidification coupling in fractured reservoirs that incorporates damage variables;

[0010] S3. Numerical simulation of thermal-fluid-solidification coupling in fractured reservoirs is carried out based on discrete fractured-vuggy network model and thermal-fluid-solidification coupling mathematical model;

[0011] S4. Based on numerical simulation results and a coupled thermal flow-solidification mathematical model, a failure criterion is introduced to evaluate wellbore stability.

[0012] Furthermore, the geometric parameters include: crack geometric parameters and cave geometric parameters, wherein the crack geometric parameters include the location, size, orientation, aperture, and distribution density of the crack; and the cave geometric parameters include the shape, location, size, orientation, and distribution density of the cave.

[0013] Furthermore, the evolution equations for the damage variables include

[0014] Shear damage evolution equation:

[0015]

[0016] Tensile damage evolution equation:

[0017]

[0018] In the formula:

[0019] This represents shear damage and is dimensionless. Represents the shear damage evolution coefficient, dimensionless;

[0020] The uniaxial compressive strength is expressed in MPa. The strain at the peak of the maximum compressive principal stress is dimensionless.

[0021] The first effective principal stress is MPa; The second effective principal stress, MPa;

[0022] The angle between rock particles resisting shear slip, in degrees;

[0023] Poisson's ratio is dimensionless. The initial elastic modulus is GPa;

[0024] Indicates tensile damage; dimensionless. The initial tensile strength is given in MPa.

[0025] Represents the residual tensile strength, where The residual strength coefficient is dimensionless.

[0026] It represents the equivalent elastic strain and is dimensionless.

[0027] It represents the elastic principal strain and is dimensionless. The first elastic principal strain is dimensionless.

[0028] The second elastic principal strain is dimensionless. The third elastic principal strain is dimensionless.

[0029] This represents the strain under tensile strength and is dimensionless. It represents the ultimate tensile strain and is dimensionless.

[0030] Furthermore, S1 includes the following steps:

[0031] S1.1. Based on fractal geometry theory, determine the quantity and distribution of cracks and karst caves;

[0032] S1.2. Based on Fisher distribution and Weibull distribution, determine the directional distribution and center location distribution of cracks and caves.

[0033] Furthermore, the thermal-fluid-solidification coupling mathematical model in S2 includes: solid field equations, fluid field equations, temperature field equations, chemical field equations, and crack deformation model that incorporate damage variables.

[0034] Furthermore, the fluid field equations include:

[0035] The seepage equations for rock matrix, fractures, caverns, and the NS equations within caverns, along with the boundary of BJS-coupled seepage and free flow regions, are introduced with damage variables.

[0036] Furthermore, the temperature field equations include:

[0037] Energy conservation equations in rock matrix, fractures, and caverns with damage variables.

[0038] Furthermore, the chemical field equations include:

[0039] Solute transport equations in porous media with damage variables, solute transport equations in cracks, and solute transport equations in karst caves.

[0040] Furthermore, S3 includes the following steps:

[0041] S3.1. Based on the physical parameters of reservoir rock matrix, fractures, caverns and fluids, the discrete fracture-cavity network model is combined with the damage-thermal-fluid-solidification coupled mathematical model to carry out numerical simulation of thermal-fluid-solidification coupled fracture-cavity reservoirs;

[0042] S3.2. Perform numerical solutions in numerical simulation software, output simulation results including temperature, pressure, stress, displacement, flow rate and solute transfer, and plot the corresponding images and curves of temperature, pressure, flow rate and solute transfer.

[0043] Furthermore, S4 includes the following steps:

[0044] S4.1. The damage variables of each part in the thermal flow-solidification coupling mathematical model are corrected by the damage evolution equation;

[0045] S4.2. Determine the formation collapse pressure by using the single weak surface criterion, the Mohr-Coulomb criterion, and the shear damage evolution criterion;

[0046] S4.3. Determine the formation fracture pressure by using the tensile fracture criterion and tensile damage evolution criterion of the rock.

[0047] S4.4. Output the cloud map of the perimeter failure area and the safe mud density window curve to complete the comprehensive evaluation of wellbore stability.

[0048] The beneficial effects of this invention are as follows:

[0049] This invention effectively solves the problems of insufficient evaluation accuracy and lagging prediction ability in existing technologies by integrating discrete fractured-vuggy network modeling, introducing multi-field coupled simulation of thermal flow solidification with damage variables, and dynamic stability evaluation. It provides a technical solution that better reflects actual geological conditions and physical processes for wellbore stability evaluation during drilling and completion of fractured-vuggy reservoirs, effectively improving the accuracy of evaluation results. Specifically:

[0050] First, this invention employs fractal geometry theory, Fisher distribution, and Weibull distribution to quantify the scale distribution, spatial orientation characteristics, and location probability of fractures and cavities, respectively, and constructs a three-dimensional discrete fracture-cavity network model. This effectively solves the defects of "fuzzy geometric features of fractures and cavities and difficulty in quantifying distribution patterns" in traditional reservoir modeling. The model constructed by this invention can realistically reproduce the heterogeneous structural characteristics of fracture-cavity reservoirs, providing a geometric carrier that is highly matched with the actual reservoir for subsequent multi-field coupled simulations. This avoids simulation errors caused by distortion in reservoir structure characterization and lays a geological foundation for the accurate calculation of subsequent multi-physics interaction processes.

[0051] Secondly, this invention constructs a multi-field coupled model of thermal flow solidification with damage correction, linking the stress field, seepage field, temperature field, and chemical field with the fracture deformation model through damage variables. It uses the damage evolution equation to dynamically correct key parameters such as elastic modulus, permeability, thermal conductivity, and solute diffusion coefficient, effectively making up for the shortcomings of the prior art in "ignoring the chemical solute transfer effect and the time-varying characteristics of damage-induced parameters". The coupled model of this invention can accurately depict the dynamic cyclic process of "temperature change, thermal stress, damage evolution, seepage channel reshaping, fluid intrusion and stress redistribution", effectively restoring the coupling mechanism of multi-field interaction and damage deterioration during drilling and completion, making the simulation results more consistent with the real physical response of the reservoir, and avoiding the physical distortion caused by the traditional simplified coupled model.

[0052] Furthermore, this invention employs a dynamic judgment method combining a single weak surface criterion, the Mohr-Coulomb criterion, and the tensile fracture criterion with the damage evolution equation, replacing the traditional static evaluation logic that only characterizes the critical state of failure. Simultaneously, it comprehensively determines the safe mud density window by integrating collapse pressure and fracture pressure and outputs the distribution of wellbore instability risk areas, effectively solving the problem of "difficulty in adapting to the gradual instability process of the wellbore" in the background technology. The evaluation method constructed by this invention can quantify the dynamic impact of damage deterioration on rock strength, effectively capturing the instability evolution path of "micro-crack initiation-expansion-through," not only providing a clear mud density operating range for drilling and completion operations but also locating high-risk areas to ensure on-site construction safety. Attached Figure Description

[0053] To more clearly illustrate the technical solutions of the embodiments of the present invention, the accompanying drawings used in the embodiments will be briefly introduced below. It should be understood that the following drawings only show some embodiments of the present invention and should not be regarded as a limitation of the scope. For those skilled in the art, other related drawings can be obtained from these drawings without creative effort.

[0054] Figure 1 This is a flowchart of the method provided by the present invention;

[0055] Figure 2 A schematic diagram of the discrete crack network geometric model for the method provided by this invention;

[0056] Figure 3 A schematic diagram of the geometric model of the karst network for the method provided by this invention;

[0057] Figure 4 A schematic diagram of the discrete slot network geometric model for the method provided by this invention;

[0058] Figure 5 A schematic diagram of the density window calculated using only the thermal-fluid-curing coupling model for the method provided by the present invention;

[0059] Figure 6 A schematic diagram of the density window calculated using the thermal flow-curing coupling model for introducing damage in the method provided by this invention. Detailed Implementation

[0060] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only a part of the embodiments of the present invention, not all of them. Therefore, the following detailed description of the embodiments of the present invention provided in the accompanying drawings is not intended to limit the scope of the claimed invention, but merely to represent selected embodiments of the invention. All other embodiments obtained by those skilled in the art based on the embodiments of the present invention without inventive effort are within the scope of protection of the present invention.

[0061] Example

[0062] This embodiment provides a wellbore stability evaluation method for fractured-vuggy reservoirs based on damage-thermal-solidification coupling, including the following steps:

[0063] S1. Construct a discrete slot network model based on the geometric parameters of the slots;

[0064] S2. Construct a mathematical model of thermal flux-solidification coupling in fractured reservoirs that incorporates damage variables;

[0065] S3. Numerical simulation of thermal-fluid-solidification coupling in fractured reservoirs is carried out based on discrete fractured-vuggy network model and thermal-fluid-solidification coupling mathematical model;

[0066] S4. Based on numerical simulation results and a thermal-fluid-solidification coupled mathematical model, a failure criterion is introduced to evaluate wellbore stability.

[0067] Furthermore, the geometric parameters include: crack geometric parameters and cave geometric parameters, wherein the crack geometric parameters include the location, size, orientation, aperture, and distribution density of the crack; and the cave geometric parameters include the shape, location, size, orientation, and distribution density of the cave.

[0068] Furthermore, the evolution equations for the damage variables include:

[0069] Shear damage evolution equation:

[0070]

[0071] Tensile damage evolution equation:

[0072]

[0073] In the formula:

[0074] This represents shear damage and is dimensionless. Represents the shear damage evolution coefficient, dimensionless;

[0075] The uniaxial compressive strength is expressed in MPa. The strain at the peak of the maximum compressive principal stress is dimensionless.

[0076] The first effective principal stress is MPa; The second effective principal stress, MPa;

[0077] The angle between rock particles resisting shear slip, in degrees;

[0078] Poisson's ratio is dimensionless. The initial elastic modulus is GPa;

[0079] Indicates tensile damage; dimensionless. The initial tensile strength is given in MPa.

[0080] Represents the residual tensile strength, where The residual strength coefficient is dimensionless.

[0081] It represents the equivalent elastic strain and is dimensionless.

[0082] It represents the elastic principal strain and is dimensionless. The first elastic principal strain is dimensionless.

[0083] The second elastic principal strain is dimensionless. The third elastic principal strain is dimensionless.

[0084] This represents the strain under tensile strength and is dimensionless. It represents the ultimate tensile strain and is dimensionless.

[0085] Furthermore, S1 includes the following steps:

[0086] S1.1. Based on fractal geometry theory, determine the quantity and distribution of cracks and karst caves;

[0087] S1.2. Based on Fisher distribution and Weibull distribution, determine the directional distribution and center location distribution of cracks and caves.

[0088] Specifically as follows:

[0089] Based on the quantity, direction, and center location distribution of cracks and caves, a discrete crack-cavity network model is constructed.

[0090] S1.1. Based on fractal geometry theory and power-law distribution, determine the quantity distribution of cracks and karst caves;

[0091] In this embodiment, considering the complexity and statistical self-similarity of crack distribution, fractal geometry theory is used as a method to describe the crack network, quantifying the power-law distribution of "crack / cavity length-number"; where the number of cracks with a length equal to or greater than l can be expressed as:

[0092] (1)

[0093] In the formula:

[0094] N is the number of cracks with a length equal to or greater than l, and is dimensionless;

[0095] C is a proportionality constant, which is dimensionless;

[0096] D f It is the fractal dimension of the fracture, which can be obtained through core observation and measurement.

[0097] The differential of equation (1) with respect to l yields:

[0098] (2)

[0099] Let the minimum and maximum lengths of the cracks in the target area be l, respectively. min and l max The length is in l min and l max Number of cracks between It can be estimated as follows:

[0100] (3)

[0101] From equations (2) and (3), the probability density function for the number of cracks can be defined as:

[0102] (4)

[0103] Define the distribution function F(l) for crack length, representing l min The total number of cracks between l and l and from l min Count up to l max The ratio of the total number of cracks; the function F(l) is as follows:

[0104] (5)

[0105] Let l∈(l min , lmax If so, then it can be inferred that F(l)∈(0,1); it is assumed that the function F(l) conforms to a certain distribution pattern; the maximum and minimum fracture lengths l are specified respectively. max l min Then the random fracture length can be determined as:

[0106] (6)

[0107] The probability density of the major and minor axes of the caves and the number of caves can be obtained using the same method, by replacing "crack length l" with "major axis length of the cave ellipse" and introducing shape parameters. Since a is the major axis and b is the minor axis, the distribution pattern, probability density, and random generation method of the major axis of the cave can be obtained similarly.

[0108] S1.2. Based on Fisher's distribution and Weibull's distribution, determine the directional distribution and central location distribution of cracks and caves;

[0109] The orientation and center location of fissures and caves also need to follow a certain distribution pattern, usually Fisher's distribution is used to describe the orientation of fissures:

[0110] (7)

[0111] In the formula:

[0112] θ represents the azimuth angle of the opening, in degrees; k is the dispersion coefficient, dimensionless.

[0113] The center location of cracks and caverns is determined using the Weibull distribution:

[0114] (8)

[0115] In the formula:

[0116] a is the major axis of the ellipse of the cave; b is the minor axis of the ellipse;

[0117] This completes the construction of the discrete slot network model.

[0118] Furthermore, the thermal-fluid-solidification coupling mathematical model in S2 includes: solid field equations, fluid field equations, temperature field equations, chemical field equations, and crack deformation model that incorporate damage variables.

[0119] When establishing the solid field equation (damage-stress field equation), the stress generated by deformation stress and pore elasticity effect of the rock matrix, the thermal stress caused by thermal expansion, and the expansion stress during the hygroscopic expansion process should be considered simultaneously, and the coefficients should be corrected to take into account the influence of damage.

[0120] The solid field equations include constitutive equations, as follows:

[0121] (9)

[0122] In the formula:

[0123] Represents total stress, in MPa; Represents total strain, dimensionless;

[0124] represents pore pressure, MPa; D represents damage variable, dimensionless.

[0125] It represents volumetric strain and is dimensionless. Reservoir temperature, expressed in °C;

[0126] This represents the mass fraction of the solute and is dimensionless.

[0127] and Represents reservoir temperature and solute mass fraction under stress-free conditions; dimensionless.

[0128] in, , , , and The definition is as follows:

[0129] (10)

[0130] In the formula:

[0131] This represents the damage-corrected shear modulus, in MPa. This represents Poisson's ratio, which is dimensionless.

[0132] E is the initial Young's modulus. Young's modulus after damage correction, MPa;

[0133] This represents the damaged rock mass bulk modulus, in MPa.

[0134] Represents the biot coefficient, which is dimensionless; The biot coefficient, dimensionless, represents the correction for chemical expansion and damage.

[0135] The bulk modulus of the material constituting the porous framework, expressed in MPa;

[0136] Expresses the molar mass of the solute, in kg / mol; The coefficient of water absorption expansion is expressed in MPa.

[0137] Represents the universal gas constant. ; This represents the solvent mass fraction and is dimensionless.

[0138] The stress coupling term, representing the effects of temperature change and chemical field, is dimensionless.

[0139] The coefficient of thermal expansion of the matrix after damage correction, in °C. -1 ;

[0140] This represents a reference value for the specific fluid entropy at the average system temperature and solute mass fraction. .

[0141] in The definition is as follows:

[0142] (11)

[0143] This completes the construction of the damage-solid field (stress field) equations.

[0144] The fluid field equations include:

[0145] The rock matrix seepage equations, fracture seepage equations, seepage equations within caverns, NS equations within caverns, and BJS condition-coupled seepage and free flow region boundaries are as follows: (Details omitted)

[0146] Rock matrix seepage equation:

[0147] (12)

[0148] In the formula:

[0149] Represents total strain, dimensionless; The fluid density is expressed in kg⋅m⁻³; t represents time in seconds.

[0150] Reservoir temperature, expressed in °C; Represents radial strain, dimensionless;

[0151] This represents the mass fraction of the solute and is dimensionless. The fluid dynamic viscosity is expressed in Pa·s.

[0152] Represents the reservoir temperature under stress-free conditions, in °C; The chemical diffusion-percolation coupling coefficient is dimensionless.

[0153] This represents the mass fraction of solute under stress-free conditions and is dimensionless.

[0154] The thermal conductivity of the matrix is ​​expressed in m² / s / K.

[0155] , , , , respectively represent the pressure-permeability coupling coefficient, temperature-permeability coupling coefficient, strain-permeability coupling coefficient, and solute-permeability coupling coefficient after damage correction;

[0156] in, , , , , and The expression is as follows:

[0157] (13)

[0158] In the formula:

[0159] Fluid density, kg⋅m -3 ; The bulk modulus of a fluid is expressed in MPa.

[0160] Represents the coefficient of thermal expansion of a fluid, in °C. -1 ; The thermal conductivity in the matrix is ​​expressed in m. 2 / s / K;

[0161] The permeability of a porous medium after damage correction is expressed in m. 2 ;

[0162] The permeability of porous media is expressed in m. 2 ; It is the permeability coefficient related to damage, and is dimensionless;

[0163] It is the width of the damaged area, in meters (m). It is the thermal conductivity related to damage, W / m / K;

[0164] Indicates matrix porosity; b0 represents the matrix porosity after damage correction, dimensionless; b0 represents the initial crack aperture under stress-free conditions, in mm.

[0165] in:

[0166] (14)

[0167] In the formula:

[0168] Represents residual porosity, dimensionless; Indicates the first principal stress, in MPa;

[0169] Indicates the second principal stress, in MPa; Indicates the third principal stress, in MPa;

[0170] The expression for the permeability coefficient related to effective stress is as follows:

[0171] (15)

[0172] In the formula:

[0173] The matrix drainage compressibility coefficient is dimensionless.

[0174] b0 represents the initial crack aperture, in mm; b r Indicates the residual crack aperture, in mm;

[0175] K n Indicates normal stiffness, GPa / m.

[0176] Crack seepage equation:

[0177] (16)

[0178] Among them, coefficient , , , and The expression is as follows:

[0179] (17)

[0180] In the formula:

[0181] Represents the dimensionless porosity of cracks after damage correction.

[0182] Represents the unit vector normal to the crack; Expressed as the flow exchange between the matrix and the fluid at the fracture surface, kg·m2 ꞏs)-1;

[0183] The permeability of the natural fissure follows the cubic law, m 2 ;

[0184] Crack aperture, mm; The permeability of the fracture after damage correction is expressed in m. 2 ;

[0185] The thermal conductivity in the crack is represented by m. 2 / s / K;

[0186] The biot coefficient of the crack after damage correction is dimensionless.

[0187] This represents the coefficient of thermal expansion of the crack after damage correction, in °C. -1 ;

[0188] The bulk modulus of the crack is expressed in MPa.

[0189] The mass conservation equation for caves:

[0190] (18)

[0191] in, It is the fluid velocity in the free-flow region;

[0192] Navier-Stokes equations in caves:

[0193] (19)

[0194] in, , It is the volumetric strain tensor;

[0195] g represents gravitational acceleration; It is the fluid pressure inside the cave, in MPa;

[0196] BJS boundary conditions for seepage and free flow coupling:

[0197] (20)

[0198] In the formula:

[0199] P p This represents the pore pressure in a porous medium, expressed in MPa.

[0200] n and τ represent the unit normal vector and unit tangential vector at the interface between the porous medium and the cavern, respectively;

[0201] u p The velocity of seepage in porous media is expressed in m·s. -1 K is the permeability in the porous medium, m 2 ;

[0202] β is the tangential drag coefficient, which is dimensionless.

[0203] The Darcy flow equations for the rock matrix region, the Darcy flow equations for the fracture region, and the free flow equations for the cave region are coupled to simulate the effect of approximate free flow in the cave within the strata.

[0204] This completes the construction of the fluid field equations (damage-seepage field equations);

[0205] In this embodiment, energy conservation equations are established for the rock matrix, cracks, and caverns, and the thermal conductivity coefficient is corrected by damage variables.

[0206] The temperature field equations include: energy conservation equations in the rock matrix, fractures, and caverns, modified according to damage evolution; specifically as follows:

[0207] The energy conservation equation in the rock matrix is:

[0208] (twenty one)

[0209] In the formula:

[0210] Indicates the temperature of the substrate, in °C; Represents the equivalent volumetric specific heat capacity, J·(m³). 3 ·K) -1 ;

[0211] Represents the porosity of the matrix after damage correction, dimensionless;

[0212] The density of a fluid is expressed in kg / m³. 3 ; The heat capacity of a fluid under constant pressure, kcal·kg -1 ·℃ -1 ;

[0213] This represents the Darcy velocity in the porous medium after damage correction, in m·s. -1 ;

[0214] The equivalent thermal conductivity after damage correction, expressed in W·m. -1 ·K-1 ;

[0215] in , The expression is as follows:

[0216] (twenty two)

[0217] Let J be the heat capacity of the matrix under constant pressure, in J·(kg·K). -1 ;

[0218] This indicates the density of the porous medium, in kg / m³. 3 ;

[0219] The thermal conductivity of a fluid is expressed in W·m⁻¹·K⁻¹.

[0220] This represents the thermal conductivity of the damaged matrix, expressed in W·m. -1 ·K -1 ;

[0221] The energy conservation equation in the crack is:

[0222] (twenty three)

[0223] In the formula:

[0224] Indicates the hydraulic aperture of the crack after damage correction, in mm;

[0225] Represents the dimensionless porosity of cracks after damage correction.

[0226] The velocity of fluid seepage in the crack is expressed in m·s. -1 ;

[0227] This represents the heat transfer between the porous medium and the crack, expressed in W·m. -3 ;

[0228] in, , The expression is as follows:

[0229]

[0230] The energy conservation equation in a cave is:

[0231]

[0232] In the formula:

[0233] The heat capacity of a fluid under constant pressure, expressed in kcal·kg⁻¹ -1 ·℃ -1 ;

[0234] This represents the heat exchange rate between porous media and fluids in caverns, expressed in W·m. -3 ;

[0235] A is the area of ​​the contact zone between the cave and the matrix, in m² 2 ;

[0236] The volume of the contact area between the cavern and the matrix is ​​m. 3 ;

[0237] The temperature of the fluid is expressed in °C.

[0238] Chemical field equations include:

[0239] Solute transport equations in porous media modified by damage evolution, solute transport equations in fractures, and solute transport equations in caverns are detailed below.

[0240] Solute transport equations in porous media:

[0241] (26)

[0242] In the formula:

[0243] This indicates the solute concentration, in mol / L. Indicates the concentration adsorbed on the pore surface, in mol / L;

[0244] The diffusion coefficient of a solute in a porous medium is represented by m. 2 / s;

[0245] The diffusion coefficient of heat in a porous medium is represented by m. 2 ·s·K -1 ;

[0246] The transport equation for solute in the crack is:

[0247] (27)

[0248] In the formula:

[0249] Indicates the concentration adsorbed on the surface of the crack, in mol / L;

[0250] The diffusion coefficient of the solute in the crack is represented by m. 2 / s;

[0251] The coefficient of heat diffusion in the crack is m. 2 ·s·K -1 ;

[0252] This represents the amount of solute exchanged in porous media and cracks, expressed in kg / (m²·s).

[0253] Equations for solute transport in caves:

[0254] (28)

[0255] In the formula:

[0256] This indicates the concentration adsorbed on the surface of the cave wall, in mol / L;

[0257] The diffusion coefficient of a solute in a cave is represented by m. 2 / s;

[0258] The coefficient of heat diffusion in a cave is m. 2 ·s·K -1 ;

[0259] This represents the amount of solute exchanged between porous media and karst caves, expressed in kg / (m²·s).

[0260] S2.2. Establish a crack deformation model that incorporates damage variables.

[0261] The aperture of a natural crack is affected by factors such as initial aperture, directional stress, and shear slip. A deformation model of the crack is established to determine the crack aperture.

[0262] The crack deformation model is as follows:

[0263]

[0264]

[0265]

[0266] In the formula:

[0267] This represents the tangential stress on the crack surface, in MPa. This represents the normal stress on the crack surface, in MPa.

[0268] and These represent the tangential and normal stiffness of the crack after damage correction, respectively, in GPa / m;

[0269] and These represent the tangential and normal displacements of the crack, respectively.

[0270] Indicates the initial aperture of the crack under stress-free conditions, in mm;

[0271] This indicates the crack opening caused by normal stress, in mm;

[0272] Indicates the additional aperture caused by crack shear dilatation, in mm;

[0273] This represents the effective normal stress on the crack surface after damage correction, in MPa.

[0274] Indicates the maximum opening of the crack under normal stress after damage correction, in mm;

[0275] Indicates the initial normal stiffness of the crack; Indicates the shear displacement of the crack;

[0276] This represents the displacement corresponding to the occurrence of shear slip; This indicates the maximum permissible slip displacement of the crack;

[0277] This represents the shear stress corresponding to the occurrence of shear slip in the crack;

[0278] , This indicates the shear dilatation angle of the crack.

[0279] The specific setup process is as follows:

[0280] The stress and displacement on the crack surface satisfy the following:

[0281] (29)

[0282] In the formula:

[0283] This indicates the tangential stress on the crack surface; This represents the normal stress on the crack surface;

[0284] and These represent the tangential and normal stiffness of the crack after damage correction, respectively.

[0285] and These represent the tangential and normal displacements of the crack, respectively.

[0286] The aperture of a natural crack consists of the initial aperture, the closing aperture caused by normal stress, and the shear dilatation aperture caused by shear slip of the natural crack.

[0287] (30)

[0288] In the formula:

[0289] Indicates the initial aperture of the crack under stress-free conditions, in mm;

[0290] This indicates the crack opening caused by normal stress, in mm;

[0291] The additional aperture caused by crack shear dilatation is expressed as follows:

[0292] (31)

[0293] (32)

[0294] In the formula:

[0295] This represents the tangential stress on the crack surface, in MPa. This represents the normal stress on the crack surface, in MPa.

[0296] and These represent the tangential and normal stiffness of the crack after damage correction, respectively, in GPa / m;

[0297] and These represent the tangential and normal displacements of the crack, respectively, in mm;

[0298] Indicates the initial aperture of the crack under stress-free conditions, in mm;

[0299] This indicates the crack opening caused by normal stress, in mm;

[0300] Indicates the additional aperture caused by crack shear dilatation, in mm;

[0301] This represents the effective normal stress on the crack surface after damage correction, in MPa.

[0302] Indicates the maximum opening of the crack under normal stress after damage correction, in mm;

[0303] The initial normal stiffness of the crack is expressed in GPa / m. Represents the shear displacement of the crack, in mm;

[0304] The displacement corresponding to the occurrence of shear slip is indicated in mm; Indicates the maximum permissible slip displacement of the crack, in mm;

[0305] This represents the shear stress, expressed in MPa, corresponding to the occurrence of shear slip in the crack.

[0306] , The shear dilatation angle of the crack is expressed in degrees.

[0307] This completes the establishment of the thermal-fluid-curing coupling mathematical model that incorporates damage variables.

[0308] S3. Numerical simulation of discrete crack models based on thermal flow solidification coupling mathematical model;

[0309] Furthermore, S3 includes the following steps:

[0310] S3.1. Based on the physical parameters of reservoir rock matrix, fractures, caverns and fluids, the discrete fracture-cavity network model is combined with the damage-thermal-fluid-solidification coupled mathematical model to carry out numerical simulation of thermal-fluid-solidification coupled fracture-cavity reservoirs;

[0311] The coupling principle is as follows:

[0312] Temperature changes first alter pore volume and fluid viscosity through thermal expansion and the heat-permeability coupling coefficient, thus affecting pore pressure and seepage. Changes in fluid pressure and chemical concentration, in turn, drive the deformation and damage evolution of the solid framework through Biot's effective stress, chemical expansion coefficient, and permeability evolution model. Strain and damage in the solid, in turn, regulate permeability, porosity, and thermal conductivity, causing changes in flow and heat transfer characteristics. Therefore, temperature, pore pressure, chemical potential, and stress form a closed-loop feedback: heat changes flow, flow affects chemistry, chemistry changes stress state, and stress, in turn, controls seepage and heat transfer paths, achieving dynamic coupling between the material's internal microstructure, physical parameters, and multiple physical boundary conditions.

[0313] Damage variables are obtained by inputting the formula for damage evolution into the expression. By adding the damage variables to the definitions of other relevant parameters, dynamic damage evolution can be achieved during iterative calculations.

[0314] S3.2. Perform numerical solutions in numerical simulation software, output simulation results including temperature, pressure, stress, displacement, flow rate, and solute transfer, and plot corresponding graphs and curves for temperature, pressure, flow rate, and solute transfer.

[0315] In this embodiment, COMSOL simulation software is selected. First, based on the model coupling logic of S3.1, the corresponding multiphysics module (solid mechanics, Darcy's law, heat transfer, rare matter transfer interface, etc.) is selected in the COMSOL simulation software. The key physical parameters of the reservoir rock matrix, fractures, caverns, and fluids (including permeability, elastic modulus, thermal conductivity, solute diffusion coefficient, etc. after damage correction) are input, and the boundary conditions of each field are configured (including BJS boundary conditions and other coupling boundary settings for free flow and seepage fields). After the configuration is completed, the COMSOL software will automatically identify the coupling relationship of each field, complete the boundary interaction calculation of free flow and seepage fields, and perform multi-field coupling numerical solution in combination with the preset damage variable iteration logic.

[0316] After the solution is completed, the output includes full simulation data including temperature, pressure, stress, displacement, flow rate, and solute transport. Based on this data, corresponding images and curves such as temperature field distribution cloud map, pressure gradient curve, flow velocity vector map, and solute concentration diffusion curve are plotted, providing a quantitative analysis basis for the subsequent wellbore stability evaluation of S4.

[0317] S4. Dynamically evaluate wellbore stability based on simulation results and damage evolution criteria.

[0318] Through damage-thermal-solidification coupled numerical simulation, the distribution of physical quantities such as temperature, pressure, stress, displacement, and seepage field of the formation around the well during the drilling process was obtained. These results reflect the dynamic response of the surrounding rock mass under the coupling effect of multiple physics fields. Based on the physical field data obtained from the above simulation, failure criteria were applied to determine the stability of the wellbore.

[0319] In numerical software (such as COMSOL), the single weak surface criterion, Mohr-Coulomb criterion, tensile fracture criterion and related damage evolution equations involved in S4 can be written as custom expressions, and the stress, strain and other field variables output by S3 can be called for calculation, thereby realizing a coherent analysis from coupled simulation to stability evaluation at the numerical level.

[0320] Relevant failure criteria are introduced to determine wellbore instability. Wellbore stability ultimately depends on the stress state of the surrounding rock and the rock's failure strength. Therefore, comparing the stress on the wellbore with the strength criteria allows for the determination of wellbore instability. Insufficient mud density leads to shear failure, resulting in well collapse or borehole narrowing; excessive mud density causes tensile failure, leading to leakage. Therefore, there exists a reasonable density range for drilling fluid use, with the upper limit being the fracturing pressure. The lower limit is the collapse pressure. .

[0321] Furthermore, S4 includes the following steps:

[0322] S4.1. Combine damage variables to determine the formation collapse pressure by using the single weak surface criterion, the Mohr-Coulomb criterion, and the shear damage evolution criterion;

[0323] S4.2. Determine the formation fracture pressure by using the tensile fracture criterion and tensile damage evolution criterion of the rock.

[0324] S4.3. Output the cloud map of the well perimeter failure area and the safe mud density window curve to complete the comprehensive evaluation of wellbore stability.

[0325] Specifically as follows:

[0326] S4.1. Combine damage variables to determine the formation collapse pressure by using the single weak surface criterion, the Mohr-Coulomb criterion, and the shear damage evolution criterion;

[0327] First, the wellbore stress distribution in the wellbore coordinate system is obtained:

[0328] (33)

[0329] In the formula:

[0330] , Represents the wellbore stress in the x-axis and y-axis directions of stratified formations, in MPa;

[0331] The shear stress around the well, expressed in MPa, represents the stress of a well-behaved stratum tangent to the xy plane.

[0332] Indicates tangent to in wellbore column coordinates Shear stress in a plane, MPa.

[0333] Based on the stress distribution around the well, the triaxial principal stress at any location on the well wall can be determined, thereby determining the maximum and minimum principal stresses.

[0334] (34)

[0335] Based on the principal stress distribution characteristics of the wellbore, determining the direction of the principal stresses and the spatial location of the bedding planes is sufficient to determine whether the wellbore is unstable. The angle between the maximum principal stress and the axial direction z is:

[0336] (35)

[0337] The angle between the maximum principal stress and the normal to the bedding plane is:

[0338] (36)

[0339] Wherein, the bedding normal vector n is determined by the dip angle of the bedding plane. and orientation (azimuth of bedding planes) )Sure:

[0340] (37)

[0341] The direction vector N of the maximum principal stress in the wellbore is determined by the wellbore angle. , well inclination angle Wellbore azimuth The angle between the maximum principal stress and the axis z Joint control, specifically in the form of:

[0342] (38)

[0343] In the formula:

[0344] i,j,k represent the three coordinate basis vectors of the wellbore coordinate system.

[0345] S4.1.2. The single weak surface criterion with damage correction should be used first to determine the shear failure of the bedding plane;

[0346] Based on the damage coefficient and the single weak surface criterion, the critical failure judgment equation for the bedding plane is obtained:

[0347] (39)

[0348] In the formula:

[0349] The cohesive force of the weak surface is MPa; The angle of friction on the weaker surface.

[0350] S4.1.3. If the weak surface failure is not satisfied, the damage-corrected Mohr-Coulomb criterion shall be used to determine the rock shear failure.

[0351] When shear failure occurs but the weak-plane failure condition is not met, the damage criterion can be described by the Mohr-Coulomb criterion:

[0352] (40)

[0353] In the formula:

[0354] The third effective principal stress, It represents uniaxial compressive strength.

[0355] S4.1.4. Introduce the shear damage evolution equation to describe the strength degradation process of rock after peak strength;

[0356] When shear damage conditions are met, the elastobrittle dual-damage constitutive relation is used to describe the failure process of intact rocks. The evolution equation of shear damage is as follows:

[0357] (41)

[0358] In the formula:

[0359] This represents the strain at the peak of the maximum compressive principal stress;

[0360] This represents the third elastic principal strain (maximum compressive principal strain).

[0361] The first effective principal stress, Indicates the second effective principal stress;

[0362] The residual strength coefficient, Indicates residual compressive strength;

[0363] For Young's modulus, Poisson's ratio, It is a cohesive force.

[0364] in, It is the shear failure initiation point obtained by deriving the formulas from Hooke's Law and the yield criterion when the failure criterion is satisfied (for example, by adding a "Safety" module to the Solid Mechanics - Linear Elasticity Materials module in COMSOL to introduce the failure criterion); its calculation method is existing technology, such as "it can be solved by COMSOL simulation software, You can add custom formulas to variables in COMSOL software. This can be obtained immediately, and will not be explained further here;

[0365] It is the strain calculated based on the current stress state. Its value is the "current actual strain value" calculated based on the actual stress state of the current engineering scenario (such as real-time load, confining pressure, etc.); its calculation method is based on existing technology (…). (It can be directly called in COMSOL), so no further explanation is needed here.

[0366] It is worth noting that in this embodiment, tensile strain is defined as positive and compressive strain as negative (consistent with the default strain sign in COMSOL simulation software). Therefore, when At that time, in actual numerical terms The absolute value is greater than The absolute value of , at which point no damage occurs.

[0367] Damage evolution correction , , The revised shear failure criterion is obtained as follows:

[0368] (42)

[0369] Indicates compressive strength. This is shear damage.

[0370] S4.2. Determine the formation fracture pressure;

[0371] S4.2.1. A preliminary judgment is made using the tensile fracture criterion;

[0372] Fracturing pressure analysis model: Formation fracturing is the result of tensile stress, using the tensile fracture criterion:

[0373] (43)

[0374] In the formula:

[0375] The first effective principal stress, This refers to tensile strength.

[0376] S4.2.2. Introduce the tensile damage evolution equation to describe the process of rock tensile strength decreasing with increasing tensile strain;

[0377] Evolution equation of tensile damage state:

[0378] (44)

[0379] In the formula:

[0380] Represents equivalent elastic strain; Indicates the elastic principal strain;

[0381] This represents the strain under tensile strength.

[0382] Represents the ultimate tensile strain, where The limiting strain coefficient;

[0383] Represents the residual tensile strength, where This is the residual strength coefficient.

[0384] Similarly, it can be corrected through damage evolution. , , The revised tensile failure criterion is obtained as follows:

[0385] (45)

[0386] For tensile damage, The first effective principal stress is MPa;

[0387] The initial tensile strength, This is the equivalent tensile stress.

[0388] S4.3. Output the cloud map of the well perimeter failure area and the safe mud density window curve to complete the comprehensive evaluation of wellbore stability.

[0389] To better illustrate this, the evaluation method of the present invention is used to evaluate the stability of a fractured-vuggy reservoir:

[0390] (1) Parameter settings;

[0391] The basic parameters required for the evaluation are set as follows: rock cohesion ranges from 4 to 35 MPa, and tensile strength ranges from 0.5 to 14 MPa; and a set of specific rock mechanical parameters are given, including Young's modulus E of 18 GPa, Poisson's ratio μ of 0.24, and formation pore pressure P0 of 27 MPa.

[0392] (2) Numerical simulation and dynamic stability assessment;

[0393] The above parameters were input into the COMSOL numerical simulation software, and the drilling process was simulated based on the established thermal-fluid-solidification coupling model that incorporates damage variables. During the simulation, the failure area and damage state of the formation around the well were determined in real time according to the failure criteria dynamically corrected by the damage variables (such as the modified Mohr-Coulomb criterion).

[0394] (3) Iterative calculations are performed to determine the critical pressure;

[0395] Adjust the drilling fluid density and other relevant engineering parameters, and repeat the simulation calculation until the critical state of wellbore instability is accurately determined, thereby obtaining the formation collapse pressure (lower limit of the safety window) and fracture pressure (upper limit of the safety window).

[0396] (4) Introduce the analysis of the safe density window pattern of damage;

[0397] Under conditions of low rock tensile strength and low cohesion, simultaneous tensile and shear failure is highly likely. The calculated collapse pressure equivalent density surface is larger than the fracture pressure equivalent density surface, therefore, there is no safe density window under these conditions. As rock strength increases, the collapse pressure equivalent density surface decreases, while the fracture pressure equivalent density surface increases, thus increasing the safe density window.

[0398] refer to Figure 6 Simulation results show that the range of the formation safety density window is significantly affected by rock damage. Specifically, the safety density window is first calculated using a thermal-fluid-solidification coupling model without damage correction. After introducing damage correction, both the fracture pressure equivalent density surface and the collapse pressure equivalent density surface decrease. Since damage has a greater impact on tensile failure, the decrease in the fracture pressure equivalent density surface is also greater, thus narrowing the safety density window.

[0399] (5) Comparison of typical working conditions and engineering guidance.

[0400] To illustrate this further, let's compare two typical operating conditions:

[0401] Operating Condition 1:

[0402] refer to Figure 5 , Figure 5 Density window calculated using only the model with thermal flow-curing coupling;

[0403] When the temperature difference between the drilling fluid and the reservoir is -60℃, the rock cohesion is 14MPa, and the rock tensile strength S t0 At 8 MPa, without introducing damage, the equivalent density of collapse pressure calculated using only the thermal-fluid-solidification coupling model is 1.28 g / cm³, and the equivalent density of rupture pressure is 1.35 g / cm³. The calculated safe density window is [1.28 g / cm³]. 3 1.35g / cm 3 ].

[0404] Operating Condition 2:

[0405] refer to Figure 6 , Figure 6 Density window calculated for a thermal-fluid-solidification coupling model incorporating damage;

[0406] Under the same negative temperature difference (-60℃), with cohesion (14MPa) and tensile strength (8MPa), the collapse pressure equivalent density calculated using the damage-corrected thermo-fluid-solidification coupling model decreased to 1.02 g / cm³, while the rupture pressure equivalent density decreased to 1.07 g / cm³. However, because damage has a more significant impact on rupture pressure, the rupture pressure equivalent density decreased more, and the safe density window decreased to [1.02 g / cm³]. 3 1.07 g / cm³ 3 ].

[0407] Damage can cause nonlinearities in mechanical parameters, leading to oscillations in the equivalent density. In summary, the damage introduced in this invention improves the accuracy of actual simulations and can provide precise guidance for drilling fluid density design.

[0408] The above description is not intended to limit the present invention in any way. Although the present invention has been disclosed above through embodiments, it is not intended to limit the present invention. Any person skilled in the art can make some modifications or alterations to the above-disclosed technical content to create equivalent embodiments without departing from the scope of the present invention. Any simple modifications, equivalent changes and alterations made to the above embodiments based on the technical essence of the present invention without departing from the scope of the present invention shall still fall within the scope of the present invention.

Claims

1. A method for evaluating wellbore stability in fractured-vuggy reservoirs based on damage-thermal-solidification coupling, characterized in that, Includes the following steps: S1. Construct a discrete slot network model based on the geometric parameters of the slots; S2. Construct a mathematical model of thermal flux-solidification coupling in fractured reservoirs that incorporates damage variables; S3. Numerical simulation of thermal-fluid-solidification coupling in fractured reservoirs is carried out based on discrete fractured-vuggy network model and thermal-fluid-solidification coupling mathematical model; S4. Based on numerical simulation results and a coupled thermal flow-solidification mathematical model, a failure criterion is introduced to evaluate wellbore stability.

2. The method according to claim 1, characterized in that, The geometric parameters include: crack geometric parameters and cave geometric parameters. Crack geometric parameters include the location, size, orientation, aperture, and distribution density of cracks; cave geometric parameters include the shape, location, size, orientation, and distribution density of caves.

3. The method according to claim 1, characterized in that, The evolution equations of the damage variables include Shear damage evolution equation: , Tensile damage evolution equation: , In the formula: This represents shear damage and is dimensionless. Represents the shear damage evolution coefficient, dimensionless; The uniaxial compressive strength is expressed in MPa. The strain at the peak of the maximum compressive principal stress is dimensionless. The first effective principal stress is MPa; The second effective principal stress, MPa; The angle between rock particles resisting shear slip, in degrees; Poisson's ratio is dimensionless. The initial elastic modulus is GPa; Indicates tensile damage; dimensionless. The initial tensile strength is given in MPa. Represents the residual tensile strength, where The residual strength coefficient is dimensionless. It represents the equivalent elastic strain and is dimensionless. It represents the elastic principal strain and is dimensionless. The first elastic principal strain is dimensionless. The second elastic principal strain is dimensionless. The third elastic principal strain is dimensionless. This represents the strain under tensile strength and is dimensionless. It represents the ultimate tensile strain and is dimensionless.

4. The method according to claim 1, characterized in that, S1 includes the following steps: S1.

1. Based on fractal geometry theory, determine the quantity and distribution of cracks and karst caves; S1.

2. Based on Fisher distribution and Weibull distribution, determine the directional distribution and center location distribution of cracks and caves.

5. The method according to claim 1, characterized in that, The thermal-fluid-solidification coupling mathematical model in S2 includes: solid field equations, fluid field equations, temperature field equations, chemical field equations, and crack deformation model, all incorporating damage variables.

6. The method according to claim 5, characterized in that, The fluid field equations include: The seepage equations for rock matrix, fractures, caverns, and the NS equations within caverns, along with the boundary of BJS-coupled seepage and free flow regions, are introduced with damage variables.

7. The method according to claim 5, characterized in that, The temperature field equations include: Energy conservation equations in rock matrix, fractures, and caverns with damage variables.

8. The method according to claim 5, characterized in that, Chemical field equations include: Solute transport equations in porous media with damage variables, solute transport equations in cracks, and solute transport equations in karst caves.

9. The method according to claim 1, characterized in that, S3 includes the following steps: S3.

1. Based on the physical parameters of reservoir rock matrix, fractures, caverns and fluids, the discrete fracture-cavity network model is combined with the damage-thermal-fluid-solidification coupled mathematical model to carry out numerical simulation of thermal-fluid-solidification coupled fracture-cavity reservoirs; S3.

2. Perform numerical solutions in numerical simulation software, output simulation results including temperature, pressure, stress, displacement, flow rate and solute transfer, and plot the corresponding images and curves of temperature, pressure, flow rate and solute transfer.

10. The method according to claim 1, characterized in that, S4 includes the following steps: S4.

1. The damage variables of each part in the thermal flow-solidification coupling mathematical model are corrected by the damage evolution equation; S4.

2. Determine the formation collapse pressure by using the single weak surface criterion, the Mohr-Coulomb criterion, and the shear damage evolution criterion; S4.

3. Determine the formation fracture pressure by using the tensile fracture criterion and tensile damage evolution criterion of the rock. S4.

4. Output the cloud map of the perimeter failure area and the safe mud density window curve to complete the comprehensive evaluation of wellbore stability.

Images