A fluid-solid coupling simulation method for CO2 huff-and-puff in shale reservoirs considering fracture evolution

Abstract

The present invention discloses a fluid-solid coupling simulation method for CO2 injection in shale reservoirs that considers fracture evolution. The method relates to the technical field of complex oil and gas reservoir development, including: determining shale geological parameters and CO2 injection parameters, and constructing a numerical model considering the mutual coupling relationship between fractures, seepage fields, and geological stress fields; constructing hydraulic fractures around the well; solving mass conservation and momentum conservation equations at each time step, and updating the model's seepage field and stress field parameters; at each time step, based on the tensile strength criterion and the Mohr-Coulomb criterion, if any form of mechanical failure of the matrix generates fractures, constructing new fractures based on a discrete fracture network model; updating the fracture aperture and permeability; and iterating repeatedly until the simulation is complete. Compared with existing static fracture models, the present invention combines fracture initiation with multi-field coupling effects to more accurately characterize the CO2 injection process, providing a reliable basis for optimizing CO2 displacement development strategies in shale reservoirs.

Description

Technical Field

[0001] The present invention relates to the technical field of complex oil and gas reservoir development, and in particular to a fluid-solid coupling simulation method for CO2 huff-and-puff in shale oil reservoirs taking fracture evolution into consideration. Background Art

[0002] Shale reservoirs are typically characterized by low porosity and low permeability: pores are primarily nanoscale micropores, often rich in organic matter, with effective porosity typically below 10% and extremely low permeability (generally well below 1 mD). In this context, the rock matrix permeability is almost negligible, and the pore-fracture system becomes the primary reservoir space and seepage pathway for oil and gas. CO2 huff-and-puff, as an effective technique for enhancing unconventional oil and gas recovery, has recently garnered widespread attention in shale reservoir development. Studies have shown that CO2 injection into the reservoir can increase crude oil mobility through mechanisms such as dissolution, volume expansion, and solution gas drive, significantly improving recovery efficiency. Furthermore, studies have found that the presence of fractures during CO2 huff-and-puff can significantly increase initial well production rates and ultimate recovery. Fractures increase the contact area between CO2 and crude oil, promoting the recovery of more oil. Therefore, CO2 huff-and-puff is particularly suitable for the development of ultra-low permeability shale reservoirs and offers significant potential for enhanced oil recovery.

[0003] However, shale reservoirs have complex geological conditions, high and uneven geostress, and strong rock brittleness. Fracture properties are controlled by factors such as fluid pressure within the fractures and geostress. Existing multi-physics coupling simulation methods primarily focus on the seepage-geostress coupling in reservoirs after CO2 injection. These methods often use the fracture network as a static input, failing to capture the generation and evolution of fractures driven by high pressure. This leads to an inadequate description of fluid-structure coupling effects, compromising simulation accuracy and engineering guidance.

[0004] In summary, there is an urgent need to develop a numerical simulation method that couples seepage and geological stress and can accurately describe the process of fracture generation and dynamic evolution, so as to accurately reproduce the changes in fracture topology and the evolution of seepage channels and provide reliable theoretical support for shale development. Summary of the Invention

[0005] To solve the above technical problems, the present invention discloses a fluid-solid coupling simulation method for CO2 huff-and-puff in shale oil reservoirs that takes into account fracture evolution. Based on discrete fracture grid technology, this method constructs a shale CO2 huff-and-puff model that takes into account fracture evolution and seepage-in-situ stress coupling. It can accurately describe the fracture evolution process during CO2 huff-and-puff, thereby providing technical support for production prediction.

[0006] To achieve the above object, the present invention adopts the following technical solutions:

[0007] A fluid-solid coupling simulation method for CO2 huff-and-puff in shale oil reservoirs considering fracture evolution includes the following steps:

[0008] s1. Determine the geological parameters of shale and CO2 huff-and-puff parameters, consider the coupling relationship between fractures, seepage fields, and geological stress fields, and construct a numerical model;

[0009] s2. Construct hydraulic fractures around the well based on the discrete fracture network model;

[0010] s3. At each time step, solve the mass conservation and momentum conservation equations and update the seepage and stress field parameters of the model;

[0011] s4. At each time step, based on the model stress field parameters updated in step s3, determine whether the matrix has undergone tensile failure or shear failure according to the tensile strength criterion and the Mohr-Coulomb criterion, respectively. If the matrix undergoes any form of mechanical failure and generates cracks, construct new cracks based on the discrete fracture network model. The aperture of the new crack is the initial crack aperture, and the permeability of the crack is determined by the crack aperture.

[0012] s5. At each time step, update the fracture aperture and permeability according to the model stress field and seepage field parameters updated in step s3;

[0013] s6. Iterate steps s3 to s5 until the simulation ends.

[0014] Optionally, in step s1, the coupling relationship between the fracture, seepage field and geological stress field is considered, and the coupling method is:

[0015] For the mass conservation equation:

[0016] ;

[0017] Where, is the time step, V is the grid volume, is the fluid porosity, is the density of phase j, where j is w, o, or g, representing the water phase, oil phase, or gas phase, respectively; is the j-phase saturation, is the mass fraction of the current grid component c in phase j, is the total number of matrix grids adjacent to the current grid, is the total number of crack grids adjacent to the current grid, is the contact area between the current grid and the lth adjacent grid, d is the straight-line distance between the center of the current grid and the center of the lth adjacent grid, is the pressure potential energy difference between the current grid and the adjacent grid, is the diffusion coefficient of phase j, is the difference in mass fraction of component c in phase j between the current grid and the lth adjacent grid, is the source and sink term of phase j, is the j-phase fluidity of the current grid, and the calculation method is:

[0018] ;

[0019] Where, is the absolute permeability of the current grid, is the relative permeability of phase j in the current grid, is the j-phase viscosity of the current grid;

[0020] For the momentum conservation equation:

[0021] ;

[0022] Where, is Poisson's ratio, is the Biot coefficient, is the total pore pressure of the current grid, is the linear thermal expansion coefficient, is the bulk modulus, is the body force, is the temperature, is the normal mean stress of the grid, The calculation method is:

[0023] ;

[0024] Where, 、 and are the maximum principal stress, the intermediate principal stress, and the minimum principal stress, respectively.

[0025] Optionally, in step s2, the method for constructing the cracks using the discrete crack network model is:

[0026] (1) For each fracture, a certain number of fracture grids are added to the model. The number of newly added fracture grids is consistent with the number of matrix grids that the fracture passes through. Each fracture grid represents a fracture segment distributed in a matrix grid. The permeability, porosity, and volume of the newly added fracture grids are consistent with the permeability, porosity, and volume of the fracture.

[0027] (2) Based on the mass conservation equation and momentum conservation equation, for each newly added fracture segment represented by the fracture grid, a mass and momentum coupling relationship is established between the fracture grid representing the fracture segment and the matrix grid where the fracture segment is located; if a fracture segment of the current fracture intersects with a fracture segment of another fracture, a mass and momentum coupling relationship is established between the fracture grids corresponding to the intersecting fracture segments;

[0028] (3) If the matrix grid where a certain fracture segment of the current fracture is located contains a source-sink term, in order to ensure the mass exchange between the source-sink term and the fracture, the same source-sink term is also added to the fracture grid representing the fracture segment.

[0029] Optionally, in step s4, whether the matrix undergoes tensile failure or shear failure is determined according to the tensile strength criterion and the Mohr-Coulomb criterion, respectively. The tensile strength criterion is:

[0030] If the stress condition of the matrix grid satisfies the following formula, the matrix grid will be tensilely damaged:

[0031] ;

[0032] Where, is the tensile strength of the rock matrix, is the minimum effective principal stress of the rock matrix; the tensile strength criterion shows that when the tensile strength of the rock is fixed, the smaller the effective stress of the rock, the more likely it is to undergo tensile failure.

[0033] The Mohr-Coulomb criterion is:

[0034] If the stress condition of the matrix grid satisfies the following formula, the matrix grid will undergo shear failure:

[0035] ;

[0036] Where, is the shear stress, For cohesion, is the normal stress, is the internal friction angle, shear stress and normal stress The calculation formulas are:

[0037] ;

[0038] ;

[0039] Where, and are the maximum total principal stress and the minimum total principal stress, respectively. is the rock failure surface and the minimum total principal stress The angle between them.

[0040] Optionally, in step s4, the permeability of the fracture is determined by the fracture aperture, and the calculation method is:

[0041] ;

[0042] Where, is the absolute permeability when the fracture is first opened, The opening of the crack when it first opens.

[0043] Optionally, in step s5, the aperture and permeability of all fractures are updated according to the model stress field and seepage field parameters updated in step s3, and the updating method is:

[0044] ;

[0045] ;

[0046] Where, is the current absolute permeability of the fracture; is the fluid-structure coupling coefficient, is the normal stress in the crack plane.

[0047] The beneficial effects of the present invention are as follows: (1) based on the tensile fracture and Mohr-Coulomb shear failure criteria, matrix failure is determined in real time and a fracture grid is automatically generated. At the same time, the mass conservation and momentum balance equations are solved, the seepage field and stress field parameters are dynamically updated, the fracture aperture and permeability are updated and incorporated into the fluid-solid coupling system, the entire process of fracture initiation and evolution under high-pressure CO2 displacement conditions is simulated, and real-time interaction between the fracture network and the fluid-solid field is achieved.

[0048] (2) By accurately characterizing the evolution of fracture topology and changes in seepage channels, it can provide reliable data support for production prediction and injection scheme optimization of CO2 huff-and-puff or other gas oil recovery processes, and has strong promotion value. BRIEF DESCRIPTION OF THE DRAWINGS

[0049] Figure 1 This is a flow chart of a fluid-solid coupling simulation method for CO2 huff-and-puff in shale oil reservoirs taking into account fracture evolution according to the present invention;

[0050] Figure 2 This is a schematic diagram of a shale CO2 huff-and-puff fluid-solid coupling model according to an embodiment of the present invention;

[0051] Figure 3 A schematic diagram of a hydraulic fracturing crack at the initial moment of a model shown in an embodiment of the present invention;

[0052] Figure 4A This is a crack distribution diagram of the model shown in an embodiment of the present invention after one throughput cycle (180 days);

[0053] Figure 4B This is a crack distribution diagram of the model shown in an embodiment of the present invention after two throughput rounds (360 days);

[0054] Figure 4C This is a crack distribution diagram of the model shown in an embodiment of the present invention after three throughput rounds (540 days);

[0055] Figure 5 This is a curve showing the average permeability variation of hydraulic fractures in the model according to an embodiment of the present invention. DETAILED DESCRIPTION

[0056] In order to make the purpose, technical solutions and advantages of the embodiments of the present invention clearer, the technical solutions in the embodiments of the present invention will be clearly and completely described below in conjunction with the drawings in the embodiments of the present invention. Obviously, the described embodiments are part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative work are within the scope of protection of the present invention. Therefore, the following detailed description of the embodiments of the present invention provided in the drawings is not intended to limit the scope of the invention for which protection is sought, but merely represents selected embodiments of the present invention. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative work are within the scope of protection of the present invention.

[0057] A CO2 huff-and-puff fluid-solid coupling simulation method for shale oil reservoirs considering fracture evolution is proposed to simulate and predict the CO2 huff-and-puff fluid-solid coupling model. Figure 1 The specific steps are as follows:

[0058] s1. Based on the geological and throughput parameters of the CO2 huff-and-puff project in Canada's Bakken shale oil reservoir (as shown in Table 1), a fluid-solid coupling numerical model for CO2 huff-and-puff was established, taking into account the mutual coupling relationship between fractures, seepage fields, and geological stress fields. Due to the lack of some mechanical parameters of Canada's Bakken shale oil reservoir, the application example of this invention uses typical shale mechanical parameters (rock tensile strength and cohesion are 0, and internal friction angle is 30°). Figure 2 As shown in the figure, X and Y represent two mutually perpendicular horizontal directions, and Z represents the vertical direction perpendicular to the X and Y directions. The model has a grid size of 41 × 41 × 3 (X × Y × Z), and each grid is 5 × 5 × 10 m. The model contains one well for oil production and CO2 injection.

[0059] Table 1 Model geological and throughput parameters

[0060]

[0061] Considering the mutual coupling relationship between fractures, seepage field and geological stress field, a numerical model is constructed; the coupling method is:

[0062] For the mass conservation equation:

[0063] ;

[0064] Where, is the time step, V is the grid volume, is the fluid porosity, is the density of phase j, where j is w, o, or g, representing the water phase, oil phase, or gas phase, respectively; is the j-phase saturation, is the mass fraction of the current grid component c in phase j, is the total number of matrix grids adjacent to the current grid, is the total number of crack grids adjacent to the current grid, is the contact area between the current grid and the lth adjacent grid, d is the straight-line distance between the center of the current grid and the center of the lth adjacent grid, is the pressure potential energy difference between the current grid and the adjacent grid, is the diffusion coefficient of phase j, is the difference in mass fraction of component c in phase j between the current grid and the lth adjacent grid, is the source and sink term of phase j, is the j-phase fluidity of the current grid, and the calculation method is:

[0065] ;

[0066] Where, is the absolute permeability of the current grid, is the relative permeability of phase j in the current grid, is the j-phase viscosity of the current grid;

[0067] For the momentum conservation equation:

[0068] ;

[0069] Where, is Poisson's ratio, is the Biot coefficient, is the total pore pressure of the current grid, is the linear thermal expansion coefficient, is the bulk modulus, is the body force, is the temperature, is the normal mean stress of the grid, The calculation method is:

[0070] ;

[0071] Where, 、 and are the maximum principal stress, the intermediate principal stress, and the minimum principal stress, respectively.

[0072] s2. Based on the discrete fracture network model, eight hydraulic fractures are constructed around the well, such as Figure 3 The method of constructing the discrete fracture network model is as follows:

[0073] (1) For each fracture, a certain number of fracture grids are added to the model. The number of newly added fracture grids is consistent with the number of matrix grids that the fracture passes through. Each fracture grid represents a fracture segment distributed in a matrix grid. The permeability, porosity, and volume of the newly added fracture grids are consistent with the permeability, porosity, and volume of the fracture.

[0074] (2) Based on the mass conservation equation and momentum conservation equation, for each newly added fracture segment represented by the fracture grid, a mass and momentum coupling relationship is established between the fracture grid representing the fracture segment and the matrix grid where the fracture segment is located; if a fracture segment of the current fracture intersects with a fracture segment of another fracture, a mass and momentum coupling relationship is established between the fracture grids corresponding to the intersecting fracture segments;

[0075] (3) If the matrix grid where a certain fracture segment of the current fracture is located contains a source-sink term, in order to ensure the mass exchange between the source-sink term and the fracture, the same source-sink term is also added to the fracture grid representing the fracture segment.

[0076] s3. At each time step, solve the mass conservation and momentum conservation equations and update the seepage field and stress field parameters of the model.

[0077] s4. At each time step, based on the model stress field parameters updated in step s3, determine whether the matrix has undergone tensile failure or shear failure according to the tensile strength criterion and the Mohr-Coulomb criterion, respectively. If the matrix undergoes any form of mechanical failure and generates cracks, construct new cracks based on the discrete fracture network model. The aperture of the new crack is the initial crack aperture, and the permeability of the crack is determined by the crack aperture.

[0078] The tensile strength criteria are:

[0079] If the stress condition of the matrix grid satisfies the following formula, the matrix grid will be tensilely damaged:

[0080] ;

[0081] Where, is the tensile strength of the rock matrix, is the minimum effective principal stress of the rock matrix; the tensile strength criterion shows that when the tensile strength of the rock is fixed, the smaller the effective stress of the rock, the more likely it is to undergo tensile failure.

[0082] The Mohr-Coulomb criterion is:

[0083] If the stress condition of the matrix grid satisfies the following formula, the matrix grid will undergo shear failure:

[0084] ;

[0085] Where, is the shear stress, For cohesion, is the normal stress, is the internal friction angle, shear stress and normal stress The calculation formulas are:

[0086] ;

[0087] ;

[0088] Where, and are the maximum total principal stress and the minimum total principal stress, respectively. is the rock failure surface and the minimum total principal stress The angle between them.

[0089] The permeability of a fracture is determined by the fracture aperture and is calculated as:

[0090] ;

[0091] Where, is the absolute permeability when the fracture is first opened, The opening of the crack when it first opens.

[0092] s5. At each time step, update the fracture aperture and permeability based on the model stress field and seepage field parameters updated in step s3. The updating method is:

[0093] ;

[0094] ;

[0095] Where, is the current absolute permeability of the fracture; is the fluid-structure coupling coefficient, is the normal stress in the crack plane.

[0096] s6. Iterate steps s3 to s5 until the simulation ends.

[0097] Figures 4A to 4C The crack distribution diagram after different throughput rounds is shown. As the simulation progresses, after one throughput round, Figure 4AAs shown in Figure 2, cracks are generated at x = 5-15 m and 185-200 m, and y = 90-120 m. After two throughput rounds, as shown in Figure 2, cracks are generated at x = 5-15 m and 185-200 m, and y = 90-120 m. Figure 4B As shown in Figure 2, the number of cracks increases significantly, and the crack range extends to x = 0 ~ 65 m and 150 ~ 205 m, and y = 80 ~ 120 m. After three throughput rounds, as shown in Figure 2, the number of cracks increases significantly, and the crack range extends to x = 0 ~ 65 m and 150 ~ 205 m, and y = 80 ~ 120 m. Figure 4C As shown in the figure, the cracks almost penetrate the entire simulation section in the x-direction, and the crack density at z = 3050-3060 m is significantly higher than that in the z = 3060-3080 m area. This is because the effective stress in the z-direction increases with depth. According to the tensile failure criterion, the shallower the formation, the more likely the rock will undergo tensile failure.

[0098] Figure 5 The average permeability of the eight hydraulic fractures around the well changes with the number of huff and puff cycles. The initial permeability of the hydraulic fractures is 8.33×10 -7 m 2 During the three gas injection cycles, the average permeability of the hydraulic fractures increased to 1.03×10 -6 m 2 , 1.43×10 -6 m 2 and 5.33×10 -6 m 2 , showing an accelerating growth trend. In summary, under fluid-solid coupling conditions, fracture formation and evolution have a significant impact on CO2 huff and puff in shale reservoirs. The method presented in this paper can accurately simulate the dynamic evolution of fracture initiation and seepage parameters under fluid-solid coupling conditions, providing important technical support for production prediction and optimized design of CO2 huff and puff in shale reservoirs.

[0099] The proposed method updates the pressure, stress, and seepage parameters of each grid in the model in real time by solving the mass and momentum conservation equations at each time step. Simultaneously, it determines whether fractures have formed based on the tensile strength criterion and the Mohr-Coulomb shear failure criterion, constructs new fractures based on a discrete fracture network model, and updates the fracture aperture and permeability in real time, thereby establishing a fracture network that truly reflects reservoir conditions. This method's advantage lies in coupling the fracture, seepage, and stress fields during CO2 stimulation in shale, as well as the dynamic evolution of fracture properties. This allows for a more accurate description of the rock and fluid dynamics within the reservoir, providing a reliable scientific basis for production forecasting and planning for shale reservoir development during CO2 stimulation.

[0100] Of course, the above description is not a limitation of the present invention, and the present invention is not limited to the above examples. Changes, modifications, additions or substitutions made by technicians in this technical field within the essential scope of the present invention should also fall within the scope of protection of the present invention.

Claims

1. A fluid-solid coupling simulation method for CO2 huff-and-puff in shale reservoirs considering fracture evolution, characterized in that: The steps include: s1. Determine the geological parameters of shale and CO2 huff-and-puff parameters, consider the coupling relationship between fractures, seepage fields, and geological stress fields, and construct a numerical model; s2. Construct hydraulic fractures around the well based on the discrete fracture network model; s3. At each time step, solve the mass conservation and momentum conservation equations and update the seepage and stress field parameters of the model; s4. At each time step, based on the model stress field parameters updated in step s3, determine whether the matrix has undergone tensile failure or shear failure according to the tensile strength criterion and the Mohr-Coulomb criterion, respectively. If the matrix undergoes any form of mechanical failure and generates cracks, construct new cracks based on the discrete fracture network model. The aperture of the new crack is the initial crack aperture, and the permeability of the crack is determined by the crack aperture. s5. At each time step, update the fracture aperture and permeability according to the model stress field and seepage field parameters updated in step s3; s6. Iterate steps s3 to s5 until the simulation ends. In step s1, the mutual coupling relationship among the fracture, seepage field and geological stress field is considered, and the coupling method is as follows: the numerical model includes the mass conservation equation and the momentum conservation equation; For the mass conservation equation: ； Where, is the time step, V is the grid volume, is the fluid porosity, is the density of phase j, where j is w, o, or g, representing the water phase, oil phase, or gas phase, respectively; is the j-phase saturation, is the mass fraction of the current grid component c in phase j, is the total number of matrix grids adjacent to the current grid, is the total number of crack grids adjacent to the current grid, is the contact area between the current grid and the lth adjacent grid, d is the straight-line distance between the center of the current grid and the center of the lth adjacent grid, is the pressure potential energy difference between the current grid and the adjacent grid, is the diffusion coefficient of phase j, is the difference in mass fraction of component c in phase j between the current grid and the lth adjacent grid, is the source and sink term of phase j, is the j-phase fluidity of the current grid, and the calculation method is: ； Where, is the absolute permeability of the current grid, is the relative permeability of phase j in the current grid, is the j-phase viscosity of the current grid; In step s2, the method of constructing the cracks in the discrete crack network model is: (1) For each crack, a certain number of crack grids are added to the model. The number of newly added crack grids is consistent with the number of matrix grids that the crack passes through. Each crack grid represents the crack segment distributed in a matrix grid. (2) According to the mass conservation equation and momentum conservation equation, for each fracture segment represented by the newly added fracture grid, a mass and momentum coupling relationship is established between the fracture grid representing the fracture segment and the matrix grid where the fracture segment is located.

2. The method for fluid-solid coupling simulation of CO2 huff-and-puff in shale oil reservoirs considering fracture evolution according to claim 1, characterized in that: In step s4, the permeability of the fracture is determined by the fracture aperture, which is calculated as follows: ； Where, is the absolute permeability when the fracture is first opened, The opening of the crack when it first opens.

3. The method for fluid-solid coupling simulation of CO2 huff-and-puff in shale oil reservoirs considering fracture evolution according to claim 1, characterized in that: In step s5, the aperture and permeability of all fractures are updated according to the model stress field and seepage field parameters updated in step s3. The updating method is: ； ； Where, is the current absolute permeability of the fracture; is the fluid-structure coupling coefficient, is the normal stress in the crack plane.

Images