Multi-field coupling method of thermal flow solidification for in-situ transformation of oil shale considering thermally induced fractures

Abstract

The invention discloses a multi-field coupling method for thermal flow solidification of oil shale in-situ conversion considering thermally induced fractures, which relates to the technical field of complex oil reservoir development, including: considering the coupling relationship between fractures and heat transfer field, seepage field, stress field and chemical reaction field, establishing a multi-field coupling model for thermal flow solidification; constructing a natural fracture network in the conversion model at the initial time step; simulating the multi-field coupling model for thermal flow solidification, and detecting the key stress data of each grid in the model at each time step; identifying the rock matrix grid for generating thermally induced fractures according to the key stress data of each grid at each time step, and constructing a thermally induced fracture grid; when a new thermally induced fracture is generated in a time step, detecting the intersection characteristics of the new fracture and the existing fracture, and merging the fractures with the same opening direction and adjacent positions; continuing the numerical simulation of the next time step until the simulation ends. The invention has a wide range of application value in the prediction of thermally induced fractures during the in-situ conversion of oil shale.

Description

Technical Field

[0001] The invention relates to the technical field of complex oil and gas reservoir development, and in particular to a multi-field coupling method for in-situ conversion of oil shale with thermal flow solidification taking thermally induced fractures into consideration. Background Art

[0002] Underground in-situ heating conversion is a key technology to achieve effective resource utilization. However, since in-situ conversion requires huge amounts of electricity to heat the reservoir, the current development cost is still relatively high. In order to accurately characterize the in-situ conversion process and rationally formulate development technology policies, it is necessary to establish an accurate in-situ conversion development multi-field coupling numerical model.

[0003] However, a key issue that affects the accuracy of numerical modeling of oil shale in-situ conversion is that during the entire development life cycle, the fracture system, which is the main channel for seepage, continues to evolve dynamically due to the complex multi-field coupling effects of heat conduction, mechanical damage, seepage and chemical reaction involved in the in-situ conversion of oil shale, and its mechanism has not yet been clarified. In addition, the current numerical simulation methods for thermal fractures are mostly focused on a single field effect, making it difficult to achieve a detailed description of the evolution law of thermal fractures under multi-field coupling. This results in a low accuracy in the characterization of oil and gas seepage during the development process, which restricts the design and optimization of the development plan.

[0004] Therefore, it is urgent to study a multi-field coupling method for in-situ conversion of oil shale with thermal flow and solidification considering thermally induced fractures. Summary of the invention

[0005] To solve the above technical problems, the present invention discloses a thermal flow solidification multi-field coupling method for in-situ conversion of oil shale taking into account thermally induced fractures. The method is based on embedded discrete fracture simulation EDFM technology. By constructing a thermal flow solidification multi-field coupling model and combining it with dynamic fracture generation and merging mechanisms, the method can accurately simulate the evolution behavior of fractures during the in-situ conversion of oil shale.

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

[0007] The multi-field coupling method of heat flow solidification for in-situ transformation of oil shale considering thermally induced fractures includes the following steps:

[0008] s1. Considering the coupling relationship between fractures and heat transfer field, seepage field, stress field and chemical reaction field, a three-dimensional oil shale in-situ conversion heat flow solidification multi-field coupling model is established;

[0009] s2. At the initial time step, a natural fracture network in the thermal flow-solidification multi-field coupling model is constructed based on the embedded discrete fracture model EDFM;

[0010] s3. Simulate the multi-field coupling model of thermal flow and solidification, and detect the key stress data of each grid in the model at each time step;

[0011] s4. In each time step, according to the stress key data of each grid obtained in step s3, the rock matrix grid where thermally induced fractures are generated is identified, and the thermally induced fracture grid is constructed based on the embedded discrete fracture model EDFM;

[0012] s5. When a new thermally induced fracture is generated in a time step, the intersection features of the new fracture and the existing fracture are detected, and the fractures with the same opening direction and adjacent positions (including natural fractures and thermally induced fractures) are merged; the numerical simulation of the next time step is continued until the simulation is completed.

[0013] Optionally, in step s1, the coupling relationship between the crack and the heat transfer field, the seepage field, the stress field and the chemical reaction field is considered, and the coupling method is:

[0014] For the mass conservation equation:

[0015] ;

[0016] Where, t is the simulation time, in seconds; V f is the fluid pore volume, in m 3 ; V is the volume of a single matrix or fracture grid block, in m 3 ρ j Indicates the density of phase j in kg / m 3 , where j is w, o, g, s, representing water phase, oil phase, gas phase, and solid phase, respectively; S j is the saturation of phase j, dimensionless; M is the difference in mass fraction of component c in phase j between the current grid and the adjacent grid, dimensionless; cj is the mass fraction of component c in phase j, dimensionless; n m is the number of matrix grids adjacent to the current grid, n f is the number of crack grids adjacent to the current grid; n r is the total number of chemical reactions; is the fluid pore volume, dimensionless; and kc are the stoichiometric numbers of component c as product and reactant in the kth chemical reaction, dimensionless; r k is the volumetric reaction rate of the kth chemical reaction, in kg / (m 3 ·s);q j is the source and sink term of phase j, in m 3 / s;λ j is the conductivity coefficient of phase j between two adjacent fractures or matrix grid blocks; Represents the potential energy difference between the current grid and its adjacent grid, in kPa, where Φ j The calculation formula is:

[0017] ;

[0018] In the formula, p j is the pore pressure of the current grid j phase, in kPa; g is the gravitational acceleration, in m / s 2 ; h is the depth of the current grid from the reference plane, in meters;

[0019] For the energy conservation equation:

[0020] ;

[0021] Where V v is the total pore volume containing fluid and solid components that can undergo chemical reactions, in m 3 ; V r is the volume of the chemically inert rock matrix, in m 3 ; U j represents the mass internal energy of phase j, in J / kg; U s Indicates the mass internal energy of the solid phase that can undergo chemical reactions, in J / kg; U r is the internal energy of the inert rock matrix, in J / m 3 ; C s is the total mass concentration of solid components that can undergo chemical reactions, in kg / m 3 ;H j is the mass enthalpy of phase j, in J / kg; κ is the comprehensive thermal conductivity of the saturated rock matrix, in W / (m·K); is the temperature difference between the current grid and its adjacent grid, in °C or K; H rk is the mass enthalpy change of the kth chemical reaction, in J / kg; q H The energy input provided to the heating well, in J / s;

[0022] For the momentum (stress) conservation equation:

[0023] ;

[0024] In the formula, is Poisson's ratio, dimensionless; α is Biot's coefficient, dimensionless; P is the total pore pressure of the current grid, in kPa; β is the linear thermal expansion coefficient, in 1 / ℃ or 1 / K; K is the bulk modulus, in kPa; F is the body force, in N / m 3 ; T is temperature, in °C or K; σ m is the normal average stress, in kPa, σ m The calculation method is:

[0025] ;

[0026] Where, σ1, σ2 and σ3 are the maximum principal stress, intermediate principal stress and minimum principal stress, respectively, and the unit is kPa;

[0027] The coupling relationship between the cracks, heat transfer field, seepage field, chemical reaction field and stress field contained in the mass conservation equation, energy conservation equation and momentum (stress) conservation equation includes:

[0028] (1) Coupling of cracks and heat transfer fields:

[0029] Fractures affect the heat transfer field: Fractures affect the distribution of heat in the reservoir by changing the thermal conductivity (such as thermal conductivity) and local thermal convection characteristics of the reservoir. For example, fractures can serve as channels for heat flow, resulting in rapid heat transfer along the fractures.

[0030] Heat transfer field affects cracks: high temperature can cause thermal expansion or thermal stress concentration in the rocks around the cracks, thereby causing the formation of new cracks or the expansion of existing cracks.

[0031] (2) Coupling of fractures and seepage fields

[0032] Fractures affect the seepage field: Fractures provide preferential channels for fluid flow, thereby changing the seepage path and overall permeability distribution of the reservoir. The morphology of the fracture network directly determines the efficiency of fluid flow.

[0033] Seepage field affects fractures: fluid migration causes changes in reservoir pressure, which in turn affects the stress distribution of the reservoir, leading to the formation of new fractures or the expansion of existing fractures.

[0034] (3) Coupling between cracks and stress fields

[0035] Fractures affect the stress field: The generation and expansion of fractures redistribute the stress field in the reservoir, resulting in local stress concentration or release.

[0036] Stress field affects cracks: changes in ground stress cause the formation of new cracks or the expansion of existing cracks.

[0037] (4) Coupling of cracks and chemical reaction fields

[0038] Cracks affect the chemical reaction field: Cracks provide channels for the rapid migration of fluids and reactants, which speeds up the transport of reactants and thus affects the rate of chemical reactions.

[0039] Chemical reaction fields affect fractures: Chemical reaction products (such as gases produced by thermal cracking) increase the pore pressure inside the reservoir, which may induce new fractures or extend existing fractures.

[0040] (5) Coupling of heat transfer field and seepage field

[0041] The heat transfer field affects the seepage field: temperature changes will change the viscosity and density of the fluid, thereby affecting the fluid migration characteristics in the seepage field.

[0042] The seepage field affects the heat transfer field: the migration of fluid takes away heat or supplements cold fluid, thereby changing the heat distribution in the reservoir.

[0043] (6) Coupling of heat transfer field and stress field

[0044] Heat transfer field affects stress field: temperature changes will cause thermal expansion or thermal contraction, thereby changing the stress distribution in the reservoir.

[0045] The stress field affects the heat transfer field: changes in the stress state will cause new cracks to form, thereby affecting the heat transfer path and efficiency.

[0046] (7) Coupling of heat transfer field and chemical reaction field

[0047] Heat transfer field affects chemical reaction field: Temperature is an important controlling factor of chemical reactions. For example, the thermal cracking rate of solid kerogen in oil shale increases significantly with increasing temperature.

[0048] The chemical reaction field affects the heat transfer field: chemical reactions release or absorb heat, thereby changing the local temperature distribution and heat conduction characteristics.

[0049] (8) Coupling of seepage field and stress field

[0050] Seepage field affects stress field: fluid migration causes changes in reservoir pressure, which in turn affects the stress distribution of the reservoir.

[0051] The stress field affects the seepage field: stress changes lead to crack generation, expansion and pore morphology changes, which in turn affect permeability and fluid flow paths.

[0052] (9) Coupling of seepage field and chemical reaction field

[0053] The seepage field affects the chemical reaction field: fluid migration drives changes in reactant concentrations, thereby affecting the rate and scope of chemical reactions.

[0054] Chemical reaction field affects the seepage field: chemical reactions change the pore structure of the reservoir or generate gas, affecting fluid flow characteristics and pressure distribution.

[0055] (10) Coupling of stress field and chemical reaction field

[0056] The stress field affects the chemical reaction field: changes in the stress state will lead to the formation of new cracks, thereby affecting the heat transfer path and efficiency as well as the porosity and permeability of the reservoir, changing the reservoir temperature and the distribution of reactants and products of chemical reactions, and indirectly regulating the chemical reaction rate.

[0057] Chemical reaction field affects stress field: Phase change caused by chemical reaction and the products such as oil, gas and water generated by it will cause changes in pore pressure, thus affecting reservoir stress distribution.

[0058] Optionally, in step s4, the principle of thermal fracture generation is: under high temperature, kerogen undergoes thermal cracking to generate products such as oil, gas, and water, resulting in a significant increase in pore pressure, thereby weakening the effective stress of the rock. This change reduces the shear strength of the rock and provides conditions for the occurrence of tensile fractures. As the pore pressure further increases, the shear stress may exceed the shear strength of the rock, or work together with the tensile stress to cause the fracture to open in the form of tension or shear.

[0059] The method to identify the rock matrix grid that generates thermally induced fractures is:

[0060] If the stress condition of the mesh satisfies the following equation, the mesh will be tensilely damaged:

[0061] ;

[0062] In the formula, σ T is the tensile strength of the rock matrix, is the minimum effective principal stress of the rock matrix;

[0063] If the stress condition of the mesh satisfies the following equation, the mesh will be sheared:

[0064] ;

[0065] In the formula, is the shear stress, c is the cohesive force, σ n is the normal stress, is the internal friction angle, shear stress and normal stress σ n The calculation formulas are:

[0066] ;

[0067] ;

[0068] Where σ1 and σ3 are the maximum total principal stress and the minimum total principal stress, respectively, and θ is the angle between the rock failure surface and the minimum total principal stress σ3.

[0069] The beneficial effects of the present invention are:

[0070] (1) The present invention fully considers the mechanical failure characteristics of oil shale at high temperature and the dynamic influence of cracks on multi-field coupling effects. It can accurately simulate the dynamic evolution process of thermally induced cracks and the mutual influence between heat conduction, mechanical failure, seepage field and chemical reaction field, which is of great significance for improving the numerical simulation accuracy of the oil heat flow solidification multi-field coupling model;

[0071] (2) The present invention provides an efficient numerical simulation tool that can accurately describe the mechanism of the evolution of thermally induced fractures on heat transfer, seepage and chemical reactions. This method has broad application value in the prediction of thermally induced fractures during the in-situ conversion of oil shale, and can provide accurate theoretical support for the production dynamic analysis and production plan design of the in-situ conversion of oil shale. BRIEF DESCRIPTION OF THE DRAWINGS

[0072] Figure 1 This is a schematic flow chart of a multi-field coupling method for in-situ conversion of oil shale with thermal flow solidification taking thermally induced fractures into consideration;

[0073] Figure 2A It is a three-dimensional schematic diagram of the multi-field coupling model I and model II of the in-situ conversion heat flow solidification of oil shale of the present invention;

[0074] Figure 2B It is a top view of the multi-field coupling model I and model II of the in-situ conversion heat flow solidification of oil shale of the present invention;

[0075] Figure 3A It is a three-dimensional schematic diagram of the multi-field coupling model III of the in-situ conversion heat flow solidification of oil shale of the present invention;

[0076] Figure 3B It is a top view of the multi-field coupling model III of the in-situ conversion heat flow solidification of oil shale of the present invention;

[0077] Figure 4 The daily oil production curves of Model I, Model II and Model III shown in the embodiment of the present invention;

[0078] Figure 5A A top view of crack distribution of Model I and Model II at the initial moment shown in an embodiment of the present invention;

[0079] Figure 5B A top view of crack distribution of Model I shown in an embodiment of the present invention after 2.5 years of production;

[0080] Figure 5C A top view of crack distribution of Model I shown in an embodiment of the present invention after 5 years of production;

[0081] Figure 5D A top view of crack distribution of Model I shown in an embodiment of the present invention after 7.5 years of production;

[0082] Figure 5E A top view of crack distribution of Model I shown in an embodiment of the present invention after 10 years of production;

[0083] Fig. 6A The three-dimensional crack distribution diagram of Model I and Model II at the initial moment shown in the embodiment of the present invention;

[0084] Figure 6B The three-dimensional crack distribution diagram of Model I shown in the embodiment of the present invention after 2.5 years of production;

[0085] Figure 6C The three-dimensional crack distribution diagram of Model I shown in the embodiment of the present invention after 5 years of production;

[0086] Fig.6D The three-dimensional crack distribution diagram of Model I shown in the embodiment of the present invention after 7.5 years of production;

[0087] Fig. 6E The three-dimensional crack distribution diagram of Model I shown in the embodiment of the present invention after 10 years of production. DETAILED DESCRIPTION

[0088] 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, rather than 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 claimed for protection, but merely represents the 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.

[0089] The present invention discloses a multi-field coupling method for in-situ conversion of oil shale with heat flow solidification considering thermal fractures, and performs upscaling on a three-dimensional in-situ conversion of oil shale with heat flow solidification multi-field coupling model, such as Figure 1 As shown, the specific steps are as follows:

[0090] s1. Considering the coupling relationship between fractures and heat transfer field, seepage field, stress field and chemical reaction field, a three-dimensional oil shale in-situ conversion heat flow solidification coupling model is established and recorded as Model I, as shown in Figure 2A and Figure 2B shown.

[0091] X and Y represent two mutually perpendicular horizontal directions, and Z represents the vertical direction perpendicular to the X and Y directions. The number of grids in the fine-scale model I is 15×15×15 (X×Y×Z), and the size of each grid is 3×3×5 meters. The model contains 7 hexagonal heating well patterns. Figure 2A The cylinder with a triangle on the top is a production well (a total of 6 wells), and the cylinder without a triangle on the top is an electric heating well (a total of 24 wells). Both the production wells and the electric heating wells are vertical wells, and the power distribution of the electric heating wells ranges from 600 to 900 W / m. For comparison, two other models with the same physical dimensions and parameters as Model I were established, denoted as Model II and Model III; Model II only contains natural fractures in the initial time step, and does not contain thermal fractures generated subsequently, such as Figure 2A and Figure 2B As shown; Model III does not contain any cracks, such as Figure 3A and Figure 3B Models I, II and III all contain a natural fracture and a kerogen thermal cracking reaction. The chemical equation of the kerogen thermal cracking reaction is shown in Table 1. The key parameters of models I, II and III are shown in Table 2.

[0092] Table 1

[0093]

[0094] Table 2

[0095]

[0096] Considering the coupling relationship between cracks and heat transfer field, seepage field, stress field and chemical reaction field, the coupling method is:

[0097] For the mass conservation equation:

[0098] ;

[0099] Where, t is the simulation time, in seconds; V f is the fluid pore volume, in m 3 ; V is the volume of a single matrix or fracture grid block, in m 3 ρ j Indicates the density of phase j in kg / m 3 , where j is w, o, g, s, representing water phase, oil phase, gas phase, and solid phase, respectively; S j is the saturation of phase j, dimensionless; M is the difference in mass fraction of component c in phase j between the current grid and the adjacent grid, dimensionless; cj is the mass fraction of component c in phase j, dimensionless; n m is the number of matrix grids adjacent to the current grid, n fis the number of crack grids adjacent to the current grid; n r is the total number of chemical reactions; is the fluid pore volume, dimensionless; and kc are the stoichiometric numbers of component c as product and reactant in the kth chemical reaction, dimensionless; r k is the volumetric reaction rate of the kth chemical reaction, in kg / (m 3 ·s);q j is the source and sink term of phase j, in m 3 / s;λ j is the conductivity coefficient of phase j between two adjacent fractures or matrix grid blocks; Represents the potential energy difference between the current grid and its adjacent grid, in kPa, where Φ j The calculation formula is:

[0100] ;

[0101] In the formula, p j is the pore pressure of the current grid j phase, in kPa; g is the gravitational acceleration, in m / s 2 ; h is the depth of the current grid from the reference plane, in meters;

[0102] For the energy conservation equation:

[0103] ;

[0104] Where V v is the total pore volume including fluid and kerogen components, in m 3 ; V r is the volume of the chemically inert rock matrix, in m 3 ; U j represents the mass internal energy of phase j, in J / kg; U s Indicates the mass internal energy of the solid phase that can undergo chemical reactions, in J / kg; U r is the internal energy of the inert rock matrix, in J / m 3 ; C s is the total mass concentration of the solid components of the kerogen component, in kg / m 3 ;H j is the mass enthalpy of phase j, in J / kg; κ is the comprehensive thermal conductivity of the saturated rock matrix, in W / (m·K); is the temperature difference between the current grid and its adjacent grid, in °C or K; H rk is the mass enthalpy change of the kth chemical reaction, in J / kg; q HThe energy input provided to the heating well, in J / s;

[0105] For the momentum (stress) conservation equation:

[0106] ;

[0107] In the formula, is Poisson's ratio, dimensionless; α is Biot's coefficient, dimensionless; P is the total pore pressure of the current grid, in kPa; β is the linear thermal expansion coefficient, in 1 / ℃ or 1 / K; K is the bulk modulus, in kPa; F is the body force, in N / m 3 ; T is temperature, in °C or K; σ m is the normal average stress, in kPa, σ m The calculation method is:

[0108] ;

[0109] Where, σ1, σ2 and σ3 are the maximum principal stress, intermediate principal stress and minimum principal stress, respectively, and the unit is kPa;

[0110] The coupling relationships among the cracks, heat transfer fields, seepage fields, chemical reaction fields and stress fields included in the mass conservation equation, energy conservation equation and momentum (stress) conservation equation include: coupling between cracks and heat transfer fields; coupling between cracks and seepage fields; coupling between cracks and stress fields; coupling between cracks and chemical reaction fields; coupling between heat transfer fields and seepage fields; coupling between heat transfer fields and stress fields; coupling between heat transfer fields and chemical reaction fields; coupling between seepage fields and stress fields; coupling between seepage fields and chemical reaction fields; coupling between stress fields and chemical reaction fields.

[0111] s2. At the initial time step, the natural fracture networks in Model I and Model II are constructed based on the embedded discrete fracture model EDFM.

[0112] s3. Run Model I and check the key stress data of each grid in Model I at each time step; similarly, run Model II and Model III directly to the end of the simulation without checking Model II and Model III at each time step.

[0113] s4. In each time step, according to the stress key data of each grid in model I obtained in step s3, the rock matrix grid that generates thermally induced fractures is identified, and the thermally induced fracture grid is constructed based on EDFM.

[0114] The method to identify the rock matrix grid that generates thermally induced fractures is:

[0115] If the stress condition of the mesh satisfies the following equation, the mesh will be tensilely damaged:

[0116] ;

[0117] In the formula, σ T is the tensile strength of the rock matrix, is the minimum effective principal stress of the rock matrix;

[0118] If the stress condition of the mesh satisfies the following equation, the mesh will be sheared:

[0119] ;

[0120] In the formula, is the shear stress, c is the cohesive force, σ n is the normal stress, is the internal friction angle, shear stress and normal stress σ n The calculation formulas are:

[0121] ;

[0122] ;

[0123] Where σ1 and σ3 are the maximum total principal stress and the minimum total principal stress, respectively, and θ is the angle between the rock failure surface and the minimum total principal stress σ3.

[0124] s5. When model I generates new cracks at a certain time step, detect the intersection features of the new cracks and the existing cracks, and merge the cracks with the same opening direction and adjacent positions (including natural cracks and thermal cracks), which is equivalent to the extension of the existing cracks to the new cracks; continue the numerical simulation of model I at the next time step until the simulation ends.

[0125] like Figure 4 As shown in the figure, there are significant differences between the daily oil production curves of model I, model II and model III, indicating that natural fractures and the dynamic evolution of fractures have a significant impact on the simulation results of in-situ conversion of oil shale. This further verifies the effectiveness of the method of the present invention in accurately describing the dynamic evolution process of fractures, which is of great significance for mastering the production dynamics of in-situ conversion of oil shale; on the other hand, by comparing the daily oil production curves of model II and model III, it can be seen that the daily oil production peak of model II containing natural fractures is higher than that of model III without natural fractures, which indicates that the presence of natural fractures is conducive to increasing the daily oil production; by comparing the daily oil production curves of model I and model II, it can be seen that the daily oil production curve of model I is significantly higher than that of model II, and due to the generation of new fractures during the simulation process, the daily oil production curve of model I on the 800th day rises faster, which indicates that the generation of new fractures is also conducive to increasing the daily oil production.

[0126] like Figures 5A to 5E , Figures 6A to 6EAs shown in the figure, with the continuous production, the number of fractures in model I gradually increases, and covers 45.33% of the area in model I after 10 years of production. The fractures are mainly concentrated in the well network area. Therefore, the impact of the dynamic evolution of fractures on the in-situ transformation process of oil shale cannot be ignored. The method of the present invention is of great significance for the production dynamic analysis and development scheme design of in-situ transformation development of oil shale. In summary, the dynamic evolution of fractures has an important influence on the in-situ transformation process of oil shale. The method of the present invention can accurately simulate the dynamic evolution process of fractures, and has important application value for the production dynamic analysis and development scheme design in the in-situ transformation development of oil shale.

[0127] 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 protection scope of the present invention.

Claims

1. A multi-field coupling method for in-situ conversion of oil shale with thermal flow solidification considering thermally induced fractures, characterized in that: The steps include: s1. Considering the coupling relationship between fractures and heat transfer field, seepage field, stress field and chemical reaction field, a three-dimensional oil shale in-situ conversion heat flow solidification multi-field coupling model is established; s2. At the initial time step, a natural fracture network in the thermal flow-solidification multi-field coupling model is constructed based on the embedded discrete fracture model; s3. Simulate the multi-field coupling model of thermal flow and solidification, and detect the key stress data of each grid in the model at each time step; s4. In each time step, according to the stress key data of each grid obtained in step s3, the rock matrix grid where thermally induced fractures are generated is identified, and the thermally induced fracture grid is constructed based on the embedded discrete fracture model; s5. When a new thermally induced crack is generated in a time step, the intersection characteristics of the new crack and the existing crack are detected, and the cracks with the same opening direction and adjacent positions are merged; the numerical simulation of the next time step is continued until the simulation is completed; In step s1, the coupling relationship between the crack and the heat transfer field, seepage field, stress field and chemical reaction field is considered, and the coupling method is: For the mass conservation equation: Where, t is the simulation time, in seconds; V f is the fluid pore volume, in m 3 ; V is the volume of a single matrix or fracture grid block, in m 3 ; ρ j Indicates the density of phase j in kg / m 3 , where j is w, o, g, s, representing water phase, oil phase, gas phase, and solid phase, respectively; S j is the saturation of phase j, dimensionless; ΔM cj M is the difference in mass fraction of component c in phase j between the current grid and the adjacent grid, dimensionless; cj is the mass fraction of component c in phase j, dimensionless; n m is the number of matrix grids adjacent to the current grid, n f is the number of crack grids adjacent to the current grid; n r is the total number of chemical reactions; is the fluid pore volume, dimensionless; s′ kc and kc are the stoichiometric numbers of component c as product and reactant in the kth chemical reaction, dimensionless; r k is the volumetric reaction rate of the kth chemical reaction, in kg / (m 3 ·s);q j is the source and sink term of phase j, in m 3 / s;λ j is the conductivity coefficient of phase j between two adjacent cracks or matrix grid blocks; ΔΦ j Represents the potential energy difference between the current grid and its adjacent grid, in kPa, where Φ j The calculation formula is: Φ j =p j -ρ j gh; In the formula, p j is the pore pressure of the current grid j phase, in kPa; g is the gravitational acceleration, in m / s 2 ; h is the depth of the current grid from the reference plane, in meters; For the energy conservation equation: Where V v is the total pore volume containing fluid and solid components that can undergo chemical reactions, in m 3 ; V r is the volume of the chemically inert rock matrix, in m 3 ; U j represents the mass internal energy of phase j, in J / kg; U s Indicates the mass internal energy of the solid phase that can undergo chemical reactions, in J / kg; U r is the internal energy of the inert rock matrix, in J / m 3 ; C s is the total mass concentration of solid components that can undergo chemical reactions, in kg / m 3 ;H j is the mass enthalpy of phase j, in J / kg; κ is the comprehensive thermal conductivity of the saturated rock matrix, in W / (m·K); ΔT is the temperature difference between the current grid and its adjacent grid, in °C or K; H rk is the mass enthalpy change of the kth chemical reaction, in J / kg; q H The energy input provided to the heating well, in J / s; For the momentum conservation equation: Where υ is Poisson's ratio, dimensionless; α is Biot's coefficient, dimensionless; P is the total pore pressure of the current grid, in kPa; β is the linear thermal expansion coefficient, in 1 / ℃ or 1 / K; K is the bulk modulus, in kPa; F is the body force, in N / m 3 ; T is temperature, in °C or K; σ m is the normal average stress, in kPa, σ m The calculation method is: s m =σ1+σ2+σ3; Where, σ1, σ2 and σ3 are the maximum principal stress, intermediate principal stress and minimum principal stress, respectively, and the unit is kPa; The coupling relationships among the cracks, heat transfer fields, seepage fields, chemical reaction fields and stress fields included in the mass conservation equation, energy conservation equation and momentum conservation equation include: coupling between cracks and heat transfer fields, coupling between cracks and seepage fields, coupling between cracks and stress fields, coupling between cracks and chemical reaction fields, coupling between heat transfer fields and seepage fields, coupling between heat transfer fields and stress fields, coupling between heat transfer fields and chemical reaction fields, coupling between seepage fields and stress fields, coupling between seepage fields and chemical reaction fields, and coupling between stress fields and chemical reaction fields.

2. The multi-field coupling method for in-situ conversion of oil shale considering thermal fractures according to claim 1, characterized in that: In step s4, the method for identifying the rock matrix grid generating thermally induced fractures is: If the stress condition of the mesh satisfies the following equation, the mesh will be tensilely damaged: s T +σ3'≤0; In the formula, σ T is the tensile strength of the rock matrix, σ3' is the minimum effective principal stress of the rock matrix; If the stress condition of the mesh satisfies the following equation, the mesh will be sheared: τ≥c+(σ n -P)tan(γ); In the formula, τ is the shear stress, c is the cohesive force, σ n is the normal stress, γ is the internal friction angle, shear stress τ and normal stress σ n The calculation formulas are: Where σ1 and σ3 are the maximum total principal stress and the minimum total principal stress, respectively, and θ is the angle between the rock failure surface and the minimum total principal stress σ3.

Images