A fluid-solid coupling reservoir numerical simulation method considering the dynamic evolution of fractures
By constructing a fractured reservoir model and an embedded discrete fracture grid, real-time coupling of the seepage field and the stress field is achieved, which solves the problems of dynamic updating of fracture geometric topology and insufficient coupling between the seepage field and the stress field in the existing technology, and improves the simulation accuracy of the oil and gas reservoir production process.
Patent Information
- Application Number
- CN202511015724.8
- Authority / Receiving Office
- CN · China
- Patent Type
- Patents(China)
- Current Assignee / Owner
- Filing Date
- 2025-07-23
- Publication Date
- 2025-09-30
- Estimated Expiration
- 2045-07-23
AI Technical Summary
Existing technologies make it difficult to simulate the dynamic updating of fracture geometric topology and the real-time coupling of seepage and stress fields at the reservoir scale, resulting in inaccurate simulation of the production process of unconventional oil and gas reservoirs.
A fluid-solid coupling reservoir numerical simulation method considering the dynamic evolution of fractures is adopted. By constructing a fractured reservoir model, a fracture grid is constructed based on the embedded discrete fracture model, and the mass conservation and momentum conservation equations are solved at each time step to determine whether the fracture is open, expanding or closing, thereby realizing real-time two-way interaction between the seepage field and the stress field.
Accurately describing the opening, expansion, and closure processes of fractures improves the accuracy of predictions of fluid displacement and effects in unconventional oil and gas reservoirs, providing a reliable basis for production prediction and program design.
Smart Images

Figure CN120524871B_ABST
Abstract
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 reservoir numerical simulation method considering the dynamic evolution of fractures. Background Art
[0002] Unconventional reservoirs (such as shales and ultra-low permeability sandstones) have low porosity and permeability, resulting in complex matrix-fracture dual-media flow characteristics. Conventional single-porosity media models struggle to accurately capture the flow coupling between fractures and the matrix. Early studies often employed dual-porosity / dual-permeability models, treating the reservoir as two parallel media: matrix and fractures. The mass transfer between the fractures and matrix is described by the seepage exchange coefficient between the wellbore and the simulation block. Such models, which account for the seepage characteristics of both microscopic matrix pores and macroscopic fracture networks at multiple scales, have been widely used in fracture history matching and production prediction. However, because they often assume that the fracture network geometry, permeability, and exchange coefficients remain constant throughout the simulation, they struggle to capture the real-time propagation, closure, and connectivity changes of induced fractures. To more accurately simulate fracture evolution, some researchers have proposed the discrete fracture network (DFN) model. This approach explicitly characterizes the fracture geometry in the form of surfaces or line elements within a numerical grid, enabling detailed simulation of fracture conduction pathways. The DFN model is highly flexible in terms of fracture distribution, strike direction, and pore volume distribution. However, it still lacks consideration of the formation stress field and fracture mechanical properties. Fractures are often regarded as static diversion structures, making it difficult to simulate the fracture generation, extension, and closure processes.
[0003] Based on the above technical status, there is an urgent need for a technical solution that can realize dynamic updating of fracture geometry topology and bidirectional real-time coupling of seepage field and stress field in reservoir-scale simulation, while taking into account both computational efficiency and accuracy, so as to improve the reliability and guidance of unconventional oil and gas reservoir and production process simulation. Summary of the Invention
[0004] To solve the above technical problems, the present invention discloses a fluid-solid coupling reservoir numerical simulation method that takes into account the dynamic evolution of fractures. This method can accurately describe the opening, expansion and closure processes of fractures during reservoir development, thereby providing technical support for production prediction.
[0005] To achieve the above object, the present invention adopts the following technical solutions:
[0006] A fluid-solid coupled reservoir numerical simulation method considering the dynamic evolution of fractures includes the following steps:
[0007] s1. Considering the conservation of mass and momentum between fractures and matrix, a fractured reservoir model is constructed.
[0008] s2. Construct a fracture grid based on the embedded discrete fracture model. All fracture grids are initially in a closed state.
[0009] s3. At each time step, solve the mass conservation and momentum conservation equations and update the seepage and stress field parameters of the model;
[0010] s4. At each time step, determine whether the crack is open, expanding, or closing based on the crack properties and the seepage field and stress field parameters updated in step s3;
[0011] s5. Construct the mass and momentum coupling relationship between the newly opened fracture mesh and the matrix and other fracture meshes, and disconnect the mass and momentum coupling relationship between the closed fracture mesh and the matrix and other fracture meshes; when constructing the mass and momentum coupling relationship between the newly opened fracture mesh and the matrix and other fracture meshes, comply with the mass conservation equation and momentum conservation equation in step s1.
[0012] s6. Iterate steps s3 to s5 until the simulation ends.
[0013] Optionally, in step s1, considering the conservation of mass and momentum between the crack and the matrix, the mass conservation equation is:
[0014] ;
[0015] Where, 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 time step, is the total number of matrix and fracture grids adjacent to the current grid; is the nabla operator, which is used to obtain the gradient of the parameters in each coordinate direction; 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, is the j-phase pressure of the current grid, is the j-phase flow rate between the current grid and other matrix and fracture grids, is the source and sink term of phase j in the current grid.
[0016] The momentum conservation equation is:
[0017] ;
[0018] 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:
[0019] ;
[0020] Where, 、 and are the maximum principal stress, the intermediate principal stress, and the minimum principal stress, respectively.
[0021] Optionally, in step s2, the method for constructing the fracture grid based on the embedded discrete fracture model is:
[0022] (1) Crack segment extraction and mesh creation: Perform geometric intersection detection on each crack and matrix unit, and intercept the line segment or surface patch of the crack in each matrix mesh that intersects with the crack to generate an independent "crack mesh" for it;
[0023] (2) Constructing fracture-matrix and fracture-fracture coupling connections: Based on the mass conservation equation and momentum conservation equation, for each fracture grid, its corresponding matrix grid is identified, and a mass and momentum coupling relationship is established between the two. At the same time, the intersection between different fracture segments is detected, and a mass and momentum coupling relationship is also established between the intersecting fracture grid units to achieve a bidirectional exchange of mass and momentum between the fracture and the matrix, and between the fractures.
[0024] (3) Synchronous allocation of source / sink terms: Scan all matrix grids. If a source / sink term such as injection or production is found, the value is assigned to the fracture grid contained in the matrix grid to ensure the consistency of the source / sink terms between the fracture and the matrix.
[0025] When the crack segment (or surface) is in a closed state, the crack mesh corresponding to the crack segment is constructed only according to step (1); when the crack segment (or surface) is in an open state, steps (1) to (3) are implemented for the crack segment.
[0026] Optionally, in step s4, whether the crack is open, extended or closed is determined by:
[0027] (1) If all the crack grids of the current crack are in a closed state before this time step, the pressure and stress relationship of all the crack grids of the current crack is calculated to determine whether there is an open crack grid; if there is an open crack grid, the currently opened crack grid is marked as the tip of the crack; if there are multiple open crack grids, the current time step is shortened, the simulation is repeated, and the crack grid that opens first is found and marked as the tip of the crack;
[0028] (2) If a crack already contains a crack tip, the pressure and stress relationship of the crack grid corresponding to the crack tip is calculated to determine whether the crack tip is closed; if the crack tip is closed, the end grid of all the crack grids currently open is marked as the crack tip, and the non-end grid is marked as the main body of the crack; if after the crack tip is closed, there is only the end grid among all the crack grids open, then the main body of the crack is not marked; if after the crack tip is closed, there is no open crack grid, then the crack tip and main body are not marked;
[0029] (3) If it is determined in step (2) that the tip of the current crack is still in an open state, the pressure and stress relationship of the crack grid in front of the crack tip along the crack expansion direction is calculated to determine whether there is an open crack grid (i.e., whether the crack is expanding). If there is an open crack grid, the currently open crack grid is marked as the tip of the crack, and the original crack tip grid is marked as the main body of the crack.
[0030] Optionally, in step s4, the crack grid is judged to be open or closed by:
[0031] If the pressure and stress conditions of the current fracture grid meet the following formula, the fracture grid is opened:
[0032] ;
[0033] Where, is the pressure of the current crack grid, is the normal stress of the current fracture grid.
[0034] If the pressure and stress conditions of the current fracture grid satisfy the following equation, the fracture grid is closed:
[0035] ;
[0036] Where, is the crack closure coefficient, when When the value is 0, the opened crack mesh will no longer be closed.
[0037] Optionally, in step s4, the opening of the newly opened crack grid is:
[0038] ;
[0039] Where, is the crack width, is the volume flow rate of the fluid in the current fracture grid, is the viscosity of the fluid in the crack, is the crack length, E is the elastic modulus of rock, is the crack height.
[0040] The permeability of the newly opened fracture grid is:
[0041] ;
[0042] Where, is the absolute permeability of the fracture.
[0043] The beneficial effect of the present invention is that the method of the present invention is applicable to both finite difference and finite element simulators, dynamically capturing the entire process of fracture initiation, expansion and closure, and accurately reflecting the topological evolution of the fracture network. The method incorporates the update of fracture aperture and permeability into the coupling system, realizing real-time two-way interaction between the seepage field and the stress field; it can improve the prediction accuracy of fluid displacement and effects in unconventional oil and gas reservoirs, and provide a reliable basis for production prediction and opening plan design. BRIEF DESCRIPTION OF THE DRAWINGS
[0044] Figure 1 This is a flow chart of a fluid-solid coupling reservoir numerical simulation method considering the dynamic evolution of fractures according to the present invention;
[0045] Figure 2 Schematic diagram of the structure of the five-point well pattern water flooding black oil model I and II shown in an embodiment of the present invention;
[0046] Figure 3A A schematic diagram of the fracture shape in the middle of the reservoir (the 10th layer of Model I) after 10 years of production in Model I according to an embodiment of the present invention;
[0047] Figure 3B A schematic diagram of the fracture shape in the middle of the reservoir after 20 years of production in Model I according to an embodiment of the present invention;
[0048] Figure 4A A schematic diagram of the fracture width in the middle of the reservoir after 10 years of production in Model I according to an embodiment of the present invention;
[0049] Figure 4B A schematic diagram of the fracture width in the middle of the reservoir after 20 years of production in Model I according to an embodiment of the present invention;
[0050] Figure 5A A schematic diagram of oil saturation in the middle of the reservoir after 20 years of production in Model I according to an embodiment of the present invention;
[0051] Figure 5B A schematic diagram of oil saturation in the middle of the reservoir after 20 years of production in Model II according to an embodiment of the present invention;
[0052] Figure 6A A schematic diagram of water saturation in the middle of the reservoir after 20 years of production in Model I according to an embodiment of the present invention;
[0053] Figure 6B A schematic diagram of water saturation in the middle of the reservoir after 20 years of production in Model II according to an embodiment of the present invention;
[0054] Figure 7 The daily oil production variation curves of the entire area for Model I and Model II shown in the embodiment of the present invention;
[0055] Figure 8 The daily gas production variation curves of the entire area for Model I and Model II shown in the embodiment of the present invention;
[0056] Figure 9 The daily water production variation curve of the entire area of Model I and Model II shown in the embodiment of the present invention;
[0057] Figure 10 The bottom hole pressure difference change curves of the production wells of Model I and Model II shown in the embodiments of the present invention. DETAILED DESCRIPTION
[0058] 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.
[0059] A fluid-solid coupling reservoir numerical simulation method considering the dynamic evolution of fractures is used to simulate and predict the black oil reservoir model with five-point well pattern water flooding. Figure 1 The specific steps are as follows:
[0060] s1. Considering the conservation of mass and momentum between fractures and matrix, a five-point well pattern water flooding black oil reservoir model is constructed and recorded as Model I. Figure 2As shown in Figure 1, 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 50×50×20 (X×Y×Z) grids, and each grid is 10×10×25 meters in size. The model contains four production wells and one water injection well, arranged in a five-point well pattern. Both the production well and the water injection well are perforated throughout the reservoir. In Model I, there is a natural fracture surrounding the water injection well, as shown in Figure 1. Figure 2 As shown by the black box surrounding the injection well, the natural fracture is initially closed. The reservoir model's geological and production parameters are shown in Table 1, and the formation crude oil physical properties are shown in Table 2. For comparison, another model, designated Model II, was constructed with the same physical dimensions and parameters as Model I but without considering the opening and propagation of natural fractures.
[0061] Table 1 Reservoir geology and production parameters
[0062]
[0063] Table 2 Physical properties of formation crude oil
[0064]
[0065] Considering the conservation of mass and momentum between the crack and the matrix, the mass conservation equation is:
[0066] ;
[0067] Where, 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 time step, is the total number of matrix and fracture grids adjacent to the current grid; is the nabla operator, which is used to obtain the gradient of the parameters in each coordinate direction; 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, is the j-phase pressure of the current grid, is the j-phase flow rate between the current grid and other matrix and fracture grids, is the source and sink term of phase j in the current grid.
[0068] The momentum conservation equation is:
[0069] ;
[0070] 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:
[0071] ;
[0072] Where, 、 and are the maximum principal stress, the intermediate principal stress, and the minimum principal stress, respectively.
[0073] s2. Construct a fracture grid based on the embedded discrete fracture model. All fracture grids are initially in a closed state. Specifically:
[0074] (1) Crack segment extraction and mesh creation: Perform geometric intersection detection on each crack and matrix unit, and intercept the line segment or surface patch of the crack in each matrix mesh that intersects with the crack to generate an independent "crack mesh" for it;
[0075] (2) Constructing fracture-matrix and fracture-fracture coupling connections: Based on the mass conservation equation and momentum conservation equation, for each fracture grid, its corresponding matrix grid is identified, and a mass and momentum coupling relationship is established between the two. At the same time, the intersection between different fracture segments is detected, and a mass and momentum coupling relationship is also established between the intersecting fracture grid units to achieve a bidirectional exchange of mass and momentum between the fracture and the matrix, and between the fractures.
[0076] (3) Synchronous allocation of source / sink terms: Scan all matrix grids. If a source / sink term such as injection or production is found, the value is assigned to the fracture grid contained in the matrix grid to ensure the consistency of the source / sink terms between the fracture and the matrix.
[0077] When the crack segment (or surface) is in a closed state, the crack mesh corresponding to the crack segment is constructed only according to step (1); when the crack segment (or surface) is in an open state, steps (1) to (3) are implemented for the crack segment.
[0078] s3. At each time step, solve the mass conservation and momentum conservation equations and update the seepage and stress field parameters of the model;
[0079] s4. At each time step, determine whether the crack is open, expanding, or closing based on the crack properties and the seepage field and stress field parameters updated in step s3;
[0080] The methods for determining whether a crack is open, expanding or closing are:
[0081] (1) If all the crack grids of the current crack are in a closed state before this time step, the pressure and stress relationship of all the crack grids of the current crack is calculated to determine whether there is an open crack grid; if there is an open crack grid, the currently opened crack grid is marked as the tip of the crack; if there are multiple open crack grids, the current time step is shortened, the simulation is repeated, and the crack grid that opens first is found and marked as the tip of the crack;
[0082] (2) If a crack already contains a crack tip, the pressure and stress relationship of the crack grid corresponding to the crack tip is calculated to determine whether the crack tip is closed; if the crack tip is closed, the end grid of all the crack grids currently open is marked as the crack tip, and the non-end grid is marked as the main body of the crack; if after the crack tip is closed, there is only the end grid among all the crack grids open, then the main body of the crack is not marked; if after the crack tip is closed, there is no open crack grid, then the crack tip and main body are not marked;
[0083] (3) If it is determined in step (2) that the tip of the current crack is still in an open state, the pressure and stress relationship of the crack grid in front of the crack tip along the crack expansion direction is calculated to determine whether there is an open crack grid (i.e., whether the crack is expanding). If there is an open crack grid, the currently open crack grid is marked as the tip of the crack, and the original crack tip grid is marked as the main body of the crack.
[0084] The method to determine whether the crack grid is open or closed is:
[0085] If the pressure and stress conditions of the current fracture grid meet the following formula, the fracture grid is opened:
[0086] ;
[0087] Where, is the pressure of the current crack grid, is the normal stress of the current fracture grid.
[0088] If the pressure and stress conditions of the current fracture grid satisfy the following equation, the fracture grid is closed:
[0089] ;
[0090] Where, is the crack closure coefficient, when When the value is 0, the opened crack mesh will no longer be closed.
[0091] The opening of the newly opened crack grid is:
[0092] ;
[0093] Where, is the crack width, is the volume flow rate of the fluid in the current fracture grid, is the viscosity of the fluid in the crack, is the crack length, E is the elastic modulus of rock, is the crack height.
[0094] The permeability of the newly opened fracture grid is:
[0095] ;
[0096] Where, is the absolute permeability of the fracture.
[0097] s5. Construct the mass and momentum coupling relationship between the newly opened fracture grid and the matrix and other fracture grids, and disconnect the mass and momentum coupling relationship between the closed fracture grid and the matrix and other fracture grids;
[0098] s6. Iterate steps s3 to s5 until the simulation ends.
[0099] Figure 3A 、 Figure 3B and Figure 4A 、 Figure 4B The shapes and widths of natural fractures in Model I at different production times are shown. Figure 3A 、 Figure 3B The grid with color value 2 represents the tip of natural fracture, the grid with color value 1 represents the main body of natural fracture, and the grid with color value 0 does not contain natural fracture. Figure 3A and Figure 3B As shown in , as the simulation progresses, the natural fracture extends toward both ends along the X-axis, reaching 100 m and 300 m in the 10th and 20th years, respectively. Figure 4A and Figure 4B As shown in Figure 2, the width distribution of natural fractures presents a pattern of high in the center and low at both ends. This is because the center of the natural fracture is the injection well, where the pressure is the highest, and therefore the fracture width is the largest.
[0100] like Figure 5A and Figure 5B As shown in Figure 2, compared with Model II which does not consider the opening and expansion of natural fractures, the water flooding range of Model I which considers the dynamic evolution of natural fractures is significantly larger than that of Model II, and the injected water in Model I is centered on the natural fracture and diffuses along the Y axis to both ends; similarly, as shown in Figure 2, Figure 6A and Figure 6BAs shown in the figure, the oil saturation in most of the central area of Model I is significantly lower than that in Model II, indicating that natural fractures have a significant impact on the simulation results of water and oil saturation.
[0101] Figure 7 、 Figure 8 and Figure 9 The daily oil production, daily gas production and daily water production curves of the entire area for Model I and Model II are shown respectively. Figure 7 and Figure 8 It can be seen that the daily oil and gas production of Model I is significantly higher than that of Model II. This is because the natural fractures in Model I increase the permeability around the injection well after opening, so the water injection volume of the injection well is larger, and the reservoir pressure can be replenished faster, thus maintaining a larger production pressure difference between the production well and the reservoir (e.g. Figure 10 (as shown in the figure), which results in higher daily oil and gas production in Model I. At the same time, higher water injection rates and reservoir pressures can displace reservoir crude oil to production wells more quickly, so the daily water production in Model I is close to zero. However, because the natural fractures in Model II are closed, the reservoir permeability is low, and the flow of reservoir crude oil to the production wells is slow. As a result, the crude oil around the production wells cannot be replenished in a timely manner. Therefore, the daily water production in Model II is higher than that in Model I. This shows that accurately describing natural fractures has a significant impact on the simulation results of oil, gas, and water production.
[0102] In summary, the dynamic evolution of fractures, such as their opening and expansion, significantly impacts simulation results such as reservoir water saturation, oil saturation, and production. The proposed method accurately describes the geometric characteristics and evolution of fractures, providing reliable technical support for production performance analysis and optimization strategies in unconventional reservoir development.
[0103] 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 reservoir numerical simulation method considering the dynamic evolution of fractures, characterized by: The steps include: s1. Considering the conservation of mass and momentum between fractures and matrix, a fractured reservoir model is constructed. s2. Detect cracks contained in the matrix and construct a crack mesh based on the embedded discrete crack model. All crack meshes are initially in a closed state. 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, determine whether the crack is open, expanding, or closing based on the crack properties and the seepage field and stress field parameters updated in step s3; s5. Construct the mass and momentum coupling relationship between the newly opened fracture grid and the matrix and other fracture grids, and disconnect the mass and momentum coupling relationship between the closed fracture grid and the matrix and other fracture grids; s6. Iterate steps s3 to s5 until the simulation ends.
2. The method for numerical simulation of fluid-solid coupled reservoirs considering the dynamic evolution of fractures according to claim 1, characterized in that: In step s1, considering the conservation of mass and momentum between the crack and the matrix, the mass conservation equation is: Where, is the fluid porosity, ρ j is the density of phase j, where j is w, o or g, representing the water phase, oil phase or gas phase respectively; S j is the saturation of phase j, t is the time step, and n is the total number of matrix and fracture grids adjacent to the current grid; is the nabla operator, which is used to obtain the gradient of the parameters in each coordinate direction; k is the absolute permeability of the current grid, k rj is the relative permeability of phase j in the current grid, μ j is the j-phase viscosity of the current grid, P j is the j-phase pressure of the current grid, q j is the j-phase flow rate between the current grid and other matrix and fracture grids, and ψ is the source and sink term of the j-phase in the current grid.
3. The method for numerical simulation of fluid-solid coupled reservoirs considering the dynamic evolution of fractures according to claim 1, characterized in that: In step s4, it is determined whether the crack is open, expanding or closing. The determination method is: (1) If all the crack grids of the current crack are in a closed state before this time step, the pressure and stress relationship of all the crack grids of the current crack is calculated to determine whether there is an open crack grid; if there is an open crack grid, the currently open crack grid is marked as the tip of the crack; If there are multiple crack meshes open, reduce the current time step, repeat the simulation, find the crack mesh that opens first, and mark it as the tip of the crack; (2) If a crack already contains a crack tip, calculate the pressure and stress relationship of the crack grid corresponding to the crack tip to determine whether the crack tip is closed; if the crack tip is closed, mark the end grid of all crack grids currently open as the crack tip, and the non-end grid as the main body of the crack; If after the crack tip is closed, the only open crack grid for that crack is the end grid, the main body of the crack will not be marked; If there is no open crack grid for the crack after the crack tip is closed, the crack tip and body are not marked; (3) If it is determined in step (2) that the tip of the current crack is still in an open state, the pressure and stress relationship of the crack grid in front of the crack tip along the crack expansion direction is calculated to determine whether there is an open crack grid, that is, whether the crack is expanding. If there is an open crack grid, the currently open crack grid is marked as the tip of the crack, and the original crack tip grid is marked as the main body of the crack.
4. The method for numerical simulation of fluid-solid coupled reservoirs considering the dynamic evolution of fractures according to claim 3, characterized in that: In step s4, the method for determining whether the crack grid is open or closed is: If the pressure and stress conditions of the current fracture grid meet the following formula, the fracture grid is opened: (P f -s n )>0; Where, P f is the pressure of the current crack grid; σ n is the normal stress of the current crack grid; If the pressure and stress conditions of the current fracture grid satisfy the following equation, the fracture grid is closed: (P f -s n )·λ<0; Where λ is the crack closure coefficient. When the value of λ is 0, the opened crack grid is no longer closed.
5. The method for numerical simulation of fluid-solid coupled reservoirs considering the dynamic evolution of fractures according to claim 1, characterized in that: In step s4, the opening of the newly opened crack grid is: Where w f is the fracture width, Q is the volume flow rate of the fluid in the current fracture grid, μ is the viscosity of the fluid in the fracture, L is the fracture length, E is the elastic modulus of the rock, and h is the fracture height; The permeability of the newly opened fracture grid is: Where k f is the absolute permeability of the fracture.