A thermal flow solidification multi-field coupling crack evolution adaptive upscaling method

Abstract

The application discloses a thermal flow solidification multi-field coupling crack evolution self-adaptive upscaling method, relates to the technical field of unconventional oil and gas reservoir development and reservoir numerical simulation, and comprises the following steps: establishing a fine scale model and a coarse scale model; establishing a simplified model and obtaining a temperature field by running; running a full-physical coarse scale model and correcting a predicted fine scale temperature; iteratively calculating a fine scale reaction rate and upscaling; self-adaptive dynamic correction of a coarse scale model reaction frequency factor; coarse scale multi-field coupling solution and thermal induced crack identification; thermal induced crack merging evolution and multi-field physical property parameter updating; time stepping and simulation termination.The application has the beneficial effect that, under the premise of ensuring fine scale calculation accuracy, the method can realize significant acceleration of multi-field coupling crack dynamic evolution simulation, and can be widely applied to in-situ conversion of oil shale, underground coal gasification and in-situ upgrading of oil sand and various underground energy engineering mine scale numerical simulation with solid phase reaction as the main control mechanism.

Description

Technical Field

[0001] This invention relates to the field of unconventional oil and gas reservoir development and reservoir numerical simulation technology, and in particular to an adaptive upscaling method for thermal flow solidification multi-field coupled fracture evolution. Background Technology

[0002] As conventional oil and gas resources gradually enter the mid-to-late stages of development, global underground energy development is accelerating its expansion towards unconventional energy sources and new underground energy systems controlled primarily by solid-phase reactions or strong solid-fluid coupling. Examples include in-situ conversion processes (ICP) of oil shale, underground coal gasification, in-situ upgrading of oil sands, underground hydrogen storage, and CO2 plume geothermal energy. These projects generally involve multi-field coupling behaviors related to heat conduction, multiphase and multi-component seepage, rock mechanical failure, and multi-stage chemical reactions. Furthermore, the fracture system, as the main transport channel for fluids and heat, undergoes continuous dynamic evolution throughout the entire development lifecycle, significantly impacting development indicators. Therefore, accurately characterizing these multi-field coupling processes and rationally formulating development technology policies necessitates the establishment of accurate numerical models for thermal flow solidification (THMC) multi-field coupling.

[0003] In-situ conversion of oil shale is one of the representative solid-phase reaction-driven scenarios in the aforementioned projects. The basic principle of ICP is as follows: heat is supplied to the reservoir through densely arranged electrically heated wells, raising the temperature of solid kerogen to the thermal decomposition temperature range of 300–400°C, causing it to decompose into flowable oil and gas, which is then extracted through production wells. This process simultaneously exhibits strongly nonlinear solid-phase pyrolysis reactions, drastic pore pressure changes, and the resulting multi-field coupled behavior of thermally induced fracture initiation and propagation. It comprehensively tests the accuracy, efficiency, and stability of multi-field coupled fracture evolution numerical simulation methods, and therefore was selected as a typical application and implementation verification case in this invention.

[0004] However, existing multi-field coupled numerical simulation methods for solid-phase reaction-driven processes suffer from the following problems: First, traditional thermo-fluid-chemical (THC) coupled models mostly ignore geostress changes and fracture dynamic evolution, treating fractures as static features or using equivalent porosity-permeability approximations, making it difficult to accurately characterize the thermally induced fracture initiation and propagation behavior induced by solid phase decomposition and a sharp increase in pore pressure. Second, while some scholars have proposed all-physics THMC four-field coupled models that can describe multi-field interactions and fracture dynamic evolution with relatively fine detail, their simulations are usually only possible on laboratory-scale models smaller than 1 meter. Due to the large number of meshes, the computational cost at the field scale is extremely high, making it difficult to meet the needs of practical engineering development applications. Third, most existing reservoir upscaling methods are designed for single-phase or multi-phase seepage problems without temperature changes and chemical reactions. A few upscaling methods for solid-phase reaction processes often ignore the fracture system or treat fractures only as static features, thus making them unsuitable for multi-field coupled simulations of thermally induced fracture dynamic evolution. Fourth, some existing upscaling methods require modification of the internal source code of the numerical simulator, while commercial numerical simulators (such as Eclipse and CMG) are limited in scope. STARS are generally closed-source software, and their users cannot access the internal source code, so the above-mentioned scaling methods cannot be applied directly.

[0005] In summary, existing multi-field coupled numerical simulation methods for solid-state reaction-driven THMC suffer from bottlenecks, such as difficulty in balancing accuracy and computational efficiency, and difficulty in using them with commercial numerical simulators. There is an urgent need to study a method for upscaling multi-field coupled numerical simulation of thermally induced crack evolution that can accurately characterize the dynamic evolution of thermally induced cracks, significantly reduce computational costs, and be directly applied to closed-source commercial simulation software. Summary of the Invention

[0006] To address the aforementioned technical problems, this invention discloses an adaptive upscaling method for the evolution of multi-field coupled cracks in thermal flow solidification. This method focuses on multi-field coupled systems driven by solid-phase reactions or strong solid-fluid coupling. It constructs a two-stage upscaling framework combining a fine-scale reaction dynamics prediction mechanism based on temperature field decoupling with a thermally induced fracture dynamic evolution mechanism based on an embedded discrete fracture model (EDFM). In the first stage, the fluid flow processes of the coarse and fine-scale models are decoupled and simplified to construct a simplified model considering only heat conduction. The simplified fine-scale temperature is then corrected and predicted based on the coarse-scale temperature difference, and the reaction rate of each fine-scale grid is iteratively calculated and upscaled to the coarse scale. In the second stage, at each time step, the reaction frequency factor of the coarse-scale model is adaptively and dynamically corrected based on the upscaled reaction rate. Furthermore, the initiation, propagation, and merging of thermally induced fractures are identified and reconstructed in real time based on the embedded discrete fracture model. This invention significantly accelerates the simulation of multi-field coupled fracture dynamic evolution while ensuring the accuracy of fine-scale calculations. It can be widely applied to field-scale numerical simulations of various underground energy engineering projects, such as in-situ oil shale conversion, underground coal gasification, and in-situ oil sands upgrading, where solid-phase reactions are the primary controlling mechanism.

[0007] To achieve the above objectives, the present invention adopts the following technical solution:

[0008] An adaptive upscaling method for multi-field coupled crack evolution in thermal flow solidification includes the following steps:

[0009] S1. Considering the coupling relationship between natural fractures and heat transfer field, seepage field, stress field and chemical reaction field, a multi-field coupled fine-scale model of in-situ conversion heat flow solidification of oil shale is established. The mesh of the fine-scale model is uniformly coarsened to obtain the corresponding coarse-scale model, and the spatial mapping relationship between the coarse-scale mesh and its internal fine-scale sub-mesh is established.

[0010] S2. Simplify the coarse-scale model and the fine-scale model respectively to obtain the simplified coarse-scale model and the simplified fine-scale model. Run the simplified coarse-scale model and the simplified fine-scale model respectively to obtain the grid temperature values ​​of each model at each time step.

[0011] S3. Run the unsimplified coarse-scale model to obtain the temperature values ​​of each coarse-scale grid at each time step. Based on the temperature difference between the simplified coarse-scale model and the unsimplified coarse-scale model at each time step, correct the temperature of each grid in the simplified fine-scale model at each time step to obtain the predicted full physical fine-scale grid temperature.

[0012] S4. Based on the corrected fine-scale grid temperature, iteratively calculate the reaction rate and kerogen concentration of each fine-scale grid at each time step, and then weight the reaction rate of all fine-scale sub-grids in each coarse-scale grid according to the pore volume to obtain the average reaction rate of the coarse-scale grid after scaling up.

[0013] S5. Initialize the multiphysics parameters of the unsimplified coarse-scale model. At each time step, adaptively and dynamically correct the reaction frequency factor of the coarse-scale model based on the average reaction rate of the coarse-scale grid obtained in step S4.

[0014] S6. Using the corrected coarse-scale model, solve the mass, energy and momentum conservation equations, and detect the key stress data of each rock matrix grid at each time step. Based on the key stress data, determine the grids that have undergone tensile or shear failure, identify the location of cracks, and construct new crack grids based on the embedded discrete crack model.

[0015] S7. When a new fracture is generated at a time step, detect the intersection characteristics between the new fracture and the existing fracture, merge fractures with the same opening direction and adjacent positions, and update multi-field physical property parameters (including non-adjacent connections between matrix and fracture, as well as fracture aperture, permeability, and rock matrix porosity).

[0016] S8, from the first The time step advances to the 1st Repeat steps S5 through S7 for each time step until the coarse-scale model simulation ends.

[0017] Furthermore, in step S1, considering the coupling relationship between the cracks and the heat transfer field, seepage field, stress field, and chemical reaction field, the coupling method for establishing the multi-field coupling model of in-situ heat flow solidification of oil shale is as follows:

[0018] For the mass conservation equation:

[0019] (1);

[0020] In the formula, This is for the simulated duration; For fluid pore volume; The volume of a single matrix or crack mesh block; express The density of the phase, where, , , , These represent the aqueous phase, oil phase, gas phase, and solid phase, respectively. for Phase saturation; The components in the current grid exist The mass fraction in the phase; Components exist The diffusion coefficient in the phase; For the current grid and the first Components of adjacent grids exist The difference in mass fraction between the phases; For the current grid and the first The contact area between adjacent cracks or matrix grids; The geometric center of the current mesh and the first The projection length of the line connecting the geometric centers of two adjacent cracks or matrix meshes in the direction of the normal vector of the two mesh contact surfaces; This represents the number of matrix grids adjacent to the current grid. This represents the number of crack meshes adjacent to the current mesh. The total number of chemical reactions; For fluid pore volume; and Components In the Stoichiometric coefficients of products and reactants in a chemical reaction; For the first The volumetric reaction rate of a chemical reaction; for The source and sink of phases; For the current grid and the first Between adjacent cracks or matrix meshes The conductivity coefficient of the phase; Indicates the current grid and its neighboring grids The phase potential energy difference is calculated using the following formula:

[0021] (2);

[0022] In the formula, For the current grid Phase pore pressure; For grid of Phase pore pressure;

[0023] For the energy conservation equation:

[0024] (3);

[0025] In the formula, It represents the total pore volume containing both fluid and chemically reactive solid components. The volume of an inert rock matrix that does not undergo chemical reactions; express The mass and internal energy of a phase; The internal energy of the mass of a solid phase that can undergo chemical reactions; The internal energy of the inert rock matrix is ​​mentioned; The mass concentration of the solid component that can undergo a chemical reaction; for Enthalpy of phase mass; The overall thermal conductivity of the saturated rock matrix; This represents the temperature difference between the current grid and its neighboring grids. For the first Mass enthalpy change of a chemical reaction; Energy input provided for heating the well;

[0026] For the momentum conservation equation:

[0027] (4);

[0028] In the formula, Poisson's ratio; Biot coefficient; This represents the total pore pressure of the current grid. It is the linear thermal expansion coefficient; It is the bulk modulus; It is a volume force; For temperature; For normal mean stress, The calculation method is as follows:

[0029] (5);

[0030] in, , and These are the maximum principal stress, intermediate principal stress, and minimum principal stress, respectively.

[0031] Furthermore, in step S2, the process of simplifying the coarse-scale model and the fine-scale model is as follows:

[0032] The porosity, permeability, and initial saturation of each fluid component in both the coarse-scale and fine-scale models were set to zero, and the production wells were shut down, leaving only the electrically heated wells supplying heat to the reservoir.

[0033] The coupling relationships between the crack, heat transfer field, seepage field, chemical reaction field, and stress field contained in the mass conservation equation, energy conservation equation, and momentum conservation equation include:

[0034] (1) Coupling between fractures and 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. For example, fractures can act as channels for heat flow, causing heat to be transferred rapidly along the fractures. Conversely, high temperature can cause thermal expansion or thermal stress concentration in the rocks around the fractures, thereby causing the generation of new fractures or the expansion of existing fractures.

[0035] (2) Coupling between fractures and 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. Conversely, fluid migration causes changes in reservoir pressure, which in turn affects the stress distribution of the reservoir, leading to the generation of new fractures or the expansion of existing fractures.

[0036] (3) Coupling between fractures and stress field: The generation and propagation of fractures will redistribute the stress field in the reservoir, leading to local stress concentration or release; conversely, changes in geostress will cause the generation of new fractures or the propagation of existing fractures.

[0037] (4) Coupling between fractures and chemical reaction fields: fractures provide channels for the rapid transport of fluids and reactants, which accelerates the transport of reactants and thus affects the rate of chemical reaction; conversely, chemical reaction products (such as gases generated by thermal cracking) increase the pore pressure inside the reservoir, which may trigger new fractures or extend existing fractures.

[0038] (5) Coupling of heat transfer field and seepage field: Temperature change will change the viscosity and density of fluid, thus affecting the fluid transport characteristics in the seepage field; conversely, fluid transport will carry away heat or replenish cold fluid, thus changing the heat distribution in the reservoir.

[0039] (6) Coupling of heat transfer field and stress field: Temperature change will cause thermal expansion or thermal contraction, which will change the stress distribution in the reservoir; conversely, the change of stress state will lead to the generation of new cracks, thereby affecting the heat transfer path and efficiency.

[0040] (7) Coupling of heat transfer field and chemical reaction field: Temperature is an important controlling factor for chemical reaction. For example, the thermal cracking rate of solid kerogen in oil shale increases significantly with increasing temperature; conversely, chemical reaction releases or absorbs heat, thereby changing the local temperature distribution and heat conduction characteristics.

[0041] (8) Coupling of seepage field and stress field: Fluid migration causes changes in reservoir pressure, which in turn affects the stress distribution of the reservoir; conversely, stress changes lead to the generation, propagation and pore morphology changes, which in turn affect permeability and fluid flow path.

[0042] (9) Coupling of seepage field and chemical reaction field: fluid migration leads to changes in reactant concentration, thereby affecting the rate and range of chemical reaction; conversely, chemical reaction changes the pore structure of reservoir or generates gas, affecting fluid flow characteristics and pressure distribution.

[0043] (10) Coupling of stress field and chemical reaction field: Changes in stress state can 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 reaction, and indirectly regulating the chemical reaction rate; conversely, the phase change caused by chemical reaction and the oil, gas and water products generated will cause changes in pore pressure, which in turn affect the reservoir stress distribution.

[0044] In step S2, the method for simplifying the coarse-scale and fine-scale models is as follows: the porosity, permeability, and initial saturation of each fluid component in both models are set to zero, and the production wells are shut down, leaving only the electrically heated wells supplying heat to the reservoir. This simplification decouples the heat transfer process from the multiphase flow process in the simplified model, allowing it to simulate only the temperature field evolution based on heat conduction. Since oil shale reservoirs inherently possess low porosity and low permeability, the contribution of multiphase flow to the overall heat transfer process is extremely limited. Therefore, the impact of this simplification on the accuracy of temperature field prediction is negligible, but it significantly reduces the computational cost of the numerical simulator.

[0045] Furthermore, in step S3, the process of correcting the temperature prediction for the simplified fine-scale grid is as follows:

[0046] Step 3.1: Based on the temperatures of each mesh in the simplified coarse-scale model obtained in step S2 at each time step, and the temperatures of each mesh in the unsimplified coarse-scale model obtained in step S3 at each time step, calculate the temperature difference between the two meshes at each time step. For the... Coarse-scale grid:

[0047] (6);

[0048] In the formula, For the unsimplified coarse-scale model and the simplified coarse-scale model in the 1st... The grid, the first Temperature difference at each time step; For the unsimplified coarse-scale model in the first The grid, the first Temperature at each time step; To simplify the coarse-scale model in the first The grid, the first Temperature at each time step;

[0049] Step 3.2: Add the temperature difference obtained in Step 3.1 to the temperature of the corresponding mesh in the simplified fine-scale model to obtain the predicted temperature of the full physical fine-scale mesh.

[0050] (7);

[0051] In the formula, For the prediction of the full physical fine scale The grid in the first Temperature at each time step; To simplify the fine-scale model The grid in the first The temperature at time step n, where the nth time step n is the temperature at ... The fine-scale grid is located at the... Within a coarse-scale grid, the two correspond to each other through the spatial mapping relationship established in step S1.

[0052] The temperature correction mechanism constructed by equations (6) and (7) accurately predicts the full physical temperature field of the fine-scale grid without performing full physical fine-scale simulation, which is the key to achieving efficient prediction of fine-scale reaction dynamics.

[0053] Further, in step S4, the process of iteratively calculating the reaction rate of each fine-scale grid and scaling up is performed:

[0054] Step 4.1: Initialize the reaction kinetic parameters of each fine-scale mesh, including the reaction frequency factor, activation energy, initial kerogen concentration, and pore volume.

[0055] Step 4.2, for the first At the nth time step, based on the corrected fine-scale mesh temperature, the Arrhenius reaction kinetic model is used to calculate the nth time step. The fine-scale grid in the first... Reaction rate at each time step:

[0056] (8);

[0057] In the formula, For the first In the nth fine-scale grid The reaction rate at each time step; The frequency factor for fine-scale reactions remains constant throughout the simulation. Activation energy; It is the ideal gas constant; The proportion of fluids and chemically reactive solid components Total pore volume of a fine-scale grid; For the first The fine-scale grid in the first... The mass concentration of the solid-phase reactants at each time step;

[0058] Step 4.3, update the following formula: Kerogen concentration at each time step in the fine-scale grid:

[0059] (9);

[0060] In the formula, For the first The fine-scale grid in the first... The mass concentration of the solid-phase reactants at each time step; For the first The time step size of each time step;

[0061] Step 4.4, from the first The time step advances to the 1st Repeat steps 4.2 to 4.3 for each time step until the simulation ends;

[0062] Step 4.5, for each coarse-scale grid The reaction rates of all fine-scale subgrids in the middle are weighted by pore volume to obtain the first... All fine-scale grids corresponding to the coarse-scale grid in the first coarse-scale grid are in the second coarse-scale grid. Average reaction rate at each time step This is the reaction rate of the coarse-scale grid after scaling up.

[0063] Furthermore, in step S5, the process of adaptively and dynamically correcting the response frequency factor of the coarse-scale model is as follows:

[0064] Step 5.1, calculate the first coarse-scale model. The grid in the first Reaction rate at each time step:

[0065] (10);

[0066] In the formula, For the first The coarse-scale grid in the first... The reaction rate at each time step; For the first The coarse-scale grid in the first... The uncorrected response frequency factor at each time step is equal to the value at the previous time step. Corrected response frequency factor (If the current time step is the initial time step at the start of the simulation, i.e., there is no previous time step, then the value is the same as the fine-scale reaction frequency factor.) Consistent); For the first The coarse-scale grid in the first... The mass concentration of the solid-phase reactants at each time step;

[0067] Step 5.2: Based on the average reaction rate of the coarse-scale mesh obtained in step S4, calculate the total pore volume of the fluid and the chemically reactive solid components in the first coarse-scale mesh in the second step. Correction factor for response frequency factor at each time step:

[0068] (11);

[0069] In the formula, For the first The coarse-scale grid in the first... The reaction frequency factor correction factor for each time step; The first one obtained in step S4 All fine-scale grids corresponding to the coarse-scale grid in the first coarse-scale grid are in the second coarse-scale grid. The average reaction rate at each time step;

[0070] Step 5.3, correct the first according to the following formula. The coarse-scale grid in the first... Response frequency factor at each time step:

[0071] (12);

[0072] In the formula, For the first The coarse-scale grid in the first... Corrected response frequency factor for each time step;

[0073] Step 5.4: Substitute the corrected reaction frequency factor into the unsimplified coarse-scale model, and proceed to step S6 for the next step. Multi-field coupled crack evolution simulation at multiple time steps. Through this adaptive dynamic correction mechanism, the local reaction rate of the coarse-scale model at each time step is strictly converged to the reaction rate of the fine-scale model after scaling up. This suppresses the simulation deviations caused by the artificially high reaction rate, temperature overshoot, and excessively rapid accumulation of pore pressure caused by coarse-scale discretization. This is the core mechanism of this invention to achieve high-precision coarse-scale multi-field coupled crack evolution simulation.

[0074] Further, in step S6, the process of identifying the rock matrix mesh that generates thermally induced fractures is as follows:

[0075] If the stress condition of the rock matrix mesh satisfies the following formula, then the mesh is considered to have undergone tensile failure:

[0076] (13);

[0077] In the formula, The tensile strength of the rock matrix, The minimum effective principal stress of the rock matrix;

[0078] If the stress condition of the rock matrix mesh satisfies the following equation, then the mesh is considered to have undergone shear failure:

[0079] (14);

[0080] In the formula, For shear stress, For cohesion, It is normal stress. It is the internal friction angle;

[0081] Shear stress and normal stress The calculation formulas are as follows:

[0082] (15);

[0083] (16);

[0084] In the formula, Rock failure surface and minimum total principal stress The angle between them.

[0085] The formation mechanism of thermally induced cracks is as follows: Under high temperature, kerogen undergoes thermal decomposition, generating oil, gas, and water products, leading to a significant increase in pore pressure, which weakens the effective stress of the rock. This change reduces the shear strength of the rock and provides conditions for tensile fracture. As the pore pressure further increases, the shear stress may exceed the shear strength of the rock, or act in conjunction with tensile stress, prompting cracks to open in tensile or shear form.

[0086] In step S7, the method for merging and evolving newly generated thermally induced fractures with existing fractures is as follows: when a new thermally induced fracture is generated at a given time step, the existing natural fractures and thermally induced fractures are traversed to determine whether the opening direction of the newly generated fracture is consistent with that of the existing fractures and whether their spatial positions are adjacent. When the newly generated fracture has the same opening direction as the existing fracture and is adjacent in position, it is merged with the existing fracture, which is equivalent to the existing fracture extending along its opening direction to the newly generated fracture. At the same time, the non-adjacent connection relationship between the matrix mesh, the fracture mesh and the wellbore is reconstructed based on the embedded discrete fracture model, and the fracture aperture, fracture permeability, rock matrix porosity, rock matrix permeability, comprehensive thermal conductivity and enthalpy values ​​of each phase are updated according to the merged fracture geometric parameters.

[0087] The beneficial effects of this invention are that, compared with the prior art, this method has the following advantages:

[0088] (1) For the first time, a two-stage upscaling numerical simulation method for solid-phase reaction-driven thermal flow solidification multi-field coupled fracture evolution mechanism based on temperature field decoupling is coupled with the thermally induced fracture dynamic evolution mechanism based on embedded discrete fracture model. By adaptively and dynamically correcting the coarse-scale reaction frequency factor at each time step, the problem of excessive reaction rate, excessive pore pressure and premature fracture initiation caused by coarse-scale discretization is effectively suppressed. The coarse-scale model maintains a high degree of consistency with the fine-scale reference results in terms of reaction dynamics, thermally induced fracture initiation timing, rock physical property evolution and final oil and gas production.

[0089] (2) By using a simplified model to predict fine-scale temperature fields through decoupling from the fluid flow process, the high computational cost of directly running a full physical fine-scale model is avoided, and the computational speed can be significantly improved while ensuring simulation accuracy. Validation results using Shell's Mahogany Demonstration Project-South (MDP-S) oil shale in-situ conversion pilot test as the benchmark model show that the method proposed in this invention can achieve a computational speedup of more than 33 times compared to full physical fine-scale simulation; at the same time, it significantly improves the convergence and robustness of the numerical model.

[0090] (3) The precise characterization of thermally induced fracture initiation, propagation, and merging enables the upscale model to accurately reproduce the dynamic control effect of the fracture network on the heat-fluid coupling mass transfer path. This invention is not only applicable to the representative underground energy development process of in-situ oil shale conversion, which is mainly controlled by solid-phase pyrolysis reaction, but can also be extended to multi-field coupled engineering problems involving solid-phase reactions, such as underground coal gasification and in-situ oil sands upgrading. It has important theoretical and engineering application value for the design of mine-scale development schemes, prediction of development indicators, and optimization of production dynamics in the above-mentioned fields. Attached Figure Description

[0091] Figure 1 This is a schematic diagram of the process of the present invention;

[0092] Figure 2 This is a three-dimensional schematic diagram of the fine-scale and coarse-scale in-situ transformation heat flow solidification multi-field coupling model of oil shale established in an embodiment of the present invention.

[0093] Figure 3 The above are comparative curves showing the change of kerogen thermal cracking reaction rate over time in the grid near the heated well using the fine-scale model, coarse-scale model, and upscale model in the embodiments of the present invention.

[0094] Figure 4 This is a schematic diagram of the pore pressure distribution of the second rock matrix on day 141 in the fine-scale, coarse-scale, and up-scale models of this invention.

[0095] Figure 5 This is a schematic diagram of the crack network distribution on day 48 of the coarse-scale model, day 141 of the fine-scale model, and day 140 of the up-scale model in an embodiment of the present invention.

[0096] Figure 6 This is a comparison chart of the daily oil production, cumulative oil production, daily gas production, and cumulative gas production curves of the coarse-scale, fine-scale, and up-scale models in this embodiment of the invention. Detailed Implementation

[0097] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only a part of the embodiments of the present invention, and not all of them. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0098] This embodiment uses the Mahogany Demonstration Project-South (MDP-S) oil shale in-situ conversion pilot test of Shell as a reference case to illustrate the method proposed in this invention.

[0099] An adaptive upscaling method for multi-field coupled crack evolution in thermal flow solidification is illustrated in the following schematic diagram: Figure 1 As shown, the specific steps are as follows:

[0100] Step S1: Establish fine-scale and coarse-scale models.

[0101] Considering the coupling relationship between natural fractures and the heat transfer field, seepage field, stress field, and chemical reaction field, a fine-scale model of multi-field coupling for in-situ heat flow solidification of oil shale is established, denoted as the fine-scale model. X and Y represent two mutually perpendicular horizontal directions, and Z represents the vertical direction perpendicular to both X and Y. The fine-scale model has 51×51×3 grids in the X, Y, and Z directions, with each grid cell measuring 0.22×0.22×8.25m. The fine-scale model uses a hexagonal well pattern and includes 6 electrically heated wells, such as... Figure 2 As shown in the cylindrical shape; one central production well, such as Figure 2 The cylinder with the top arrow is shown in the middle. All heating and production wells are vertical wells, and the heating well power is 700W / m.

[0102] The mesh of the fine-scale model is uniformly coarsened to obtain the corresponding coarse-scale model. The coarse-scale model has 17×17×3 meshes in the X, Y, and Z directions, and the size of a single coarse-scale mesh is 0.67×0.67×8.25m. Figure 2The cylinder with the top arrow is shown in the image. Establish the spatial mapping relationship between the coarse-scale grid and its internal fine-scale subgrids: each coarse-scale grid corresponds to 3×3=9 fine-scale subgrids in the X and Y directions, and a one-to-one correspondence in the Z direction.

[0103] Both the fine-scale and coarse-scale models employ adiabatic boundary conditions and include a kerogen thermal cracking reaction and multiple stages of hydrocarbon component cracking reactions. The main hydraulic, thermodynamic, and rock mechanics parameters involved in the models are shown in Table 1, and the chemical equations, reaction frequency factors, and activation energy parameters for the kerogen thermal cracking reaction and the multiple stages of hydrocarbon component cracking reactions are shown in Table 2.

[0104] Table 1. Main parameters of the multi-field coupled model for in-situ thermal solidification of oil shale.

[0105]

[0106] Table 2. Equations and kinetic parameters for kerogen thermal cracking and multi-stage oil and gas cracking.

[0107]

[0108] Step S2: Build a simplified model and run it to obtain the temperature field.

[0109] The fine-scale and coarse-scale models were simplified by setting matrix porosity, fracture porosity, matrix permeability, fracture permeability, and initial saturation of each fluid phase to zero, and shutting down production wells, leaving only electrically heated wells supplying heat to the reservoir. This resulted in simplified fine-scale and coarse-scale models. Running the simplified coarse-scale and fine-scale models yielded the grid temperature values ​​at each time step. Because the fluid flow process is completely decoupled after simplification, the model only simulates the temperature field evolution based on heat conduction, thus significantly shortening the simulation time compared to a full-physics simulation.

[0110] Step S3: Run the full physical coarse-scale model and correct the predicted fine-scale temperature.

[0111] Run the unsimplified, fully physical coarse-scale model to obtain the temperature values ​​of each coarse-scale grid at each time step. For each coarse-scale grid The temperature difference between the full physical coarse-scale model and the simplified coarse-scale model is calculated according to equation (6). Furthermore, for each fine-scale subgrid falling into this coarse-scale grid... The temperature is simplified to a fine scale according to equation (7). Subtracting the temperature difference mentioned above yields the predicted temperature of the full physical fine-scale mesh. .

[0112] Because oil shale reservoirs are characterized by low porosity and low permeability, the contribution of multiphase flow to the overall heat transfer process is extremely limited. Therefore, the temperature difference between the simplified model and the full physical model on the same coarse-scale grid is minimal. It can accurately characterize the local correction of the temperature field by the flow effect, so that Equation (7) can accurately predict the full physical temperature field of the fine-scale grid without performing full physical fine-scale simulation.

[0113] Step S4: Iteratively calculate the fine-scale reaction rate and scale up.

[0114] Initialize the reaction kinetic parameters for each fine-scale mesh. Begin, for each fine-scale grid Based on the corrected temperature The first reaction was calculated according to the Arrhenius nonequilibrium reaction kinetic formula (8). The reaction rate at time step 1; update the reaction rate at time step 2 according to equation (9). The concentration of kerogen at each time step. Iterate the above process until the simulation ends.

[0115] For each coarse-scale grid The reaction rate of all its fine-scale subgrids The weighted average based on pore volume yields the... All fine-scale grids corresponding to the coarse-scale grid in the first coarse-scale grid are in the second coarse-scale grid. Average reaction rate at each time step This is the reaction rate of the coarse-scale grid after scaling up.

[0116] Step S5: Adaptive dynamic correction of the frequency factor of the coarse-scale model response.

[0117] Initialize the unsimplified coarse-scale model multiphysics parameters at each time step. :

[0118] (5.1) Calculate the coarse-scale model according to equation (10). The grid in the first The reaction rate at each time step Among them, the uncorrected reaction frequency factor The value is equal to the previous time step. Corrected response frequency factor If the current time step is the initial time step at the start of the simulation, i.e., there is no previous time step, then Values ​​and fine-scale response frequency factors Consistent;

[0119] (5.2) Calculate the reaction frequency factor correction factor according to formula (11). ;

[0120] (5.3) Calculate the corrected response frequency factor according to equation (12). ;

[0121] (5.4) Substitute the corrected reaction frequency factor into the unsimplified coarse-scale model and proceed to step S6 for the next step. Multi-field coupled crack evolution simulation at each time step.

[0122] This adaptive dynamic correction mechanism ensures that the reaction rate of the coarse-scale model at each time step strictly converges to the reaction rate of the fine-scale model after scaling up, thereby suppressing the problems of artificially high reaction rates, temperature overshoot, non-physical accumulation of pore pressure, and premature crack initiation caused by coarse-scale discretization.

[0123] Step S6: Coarse-scale multi-field coupling solution and thermal crack identification.

[0124] Solving the coarse-scale model after reaction frequency factor correction on the 1st The mass conservation equation (1), energy conservation equation (3), and momentum conservation equation (4) are obtained for each time step; the coarse-scale rock matrix grid is obtained in the 1st time step. Key stress data at each time step, including the maximum total principal stress. Intermediate total principal stress Minimum total principal stress Minimum effective principal stress and pore pressure P .

[0125] For each coarse-scale rock matrix grid, tensile failure is determined according to Equation (13); shear failure is determined according to Equations (14) to (16). When the tensile failure criterion is met, the crack initiation direction is related to the minimum effective principal stress. The direction is perpendicular; when the shear failure criterion is met, the crack initiation direction is perpendicular to the minimum total principal stress. Angle Based on the embedded discrete crack model, new crack meshes are constructed within the meshes identified as having failed.

[0126] Step S7: Thermally induced crack merging and evolution and multi-field property parameter update.

[0127] For the For each time step that generates a new crack, iterate through the existing cracks and determine whether the crack's opening direction is consistent with the existing cracks and whether the cracks are adjacent in space. If they are consistent and adjacent, merge the newly generated crack with the existing crack, which is equivalent to the existing crack extending along its opening direction to the newly generated crack. Otherwise, retain it as an independent new crack.

[0128] Based on the merged and evolved fracture topology, the non-adjacent connections between the matrix mesh, fracture mesh, and wellbore are reconstructed. Furthermore, based on the solid mass loss, pore pressure changes, and fracture geometric parameter changes caused by kerogen pyrolysis, multi-field physical properties (including matrix porosity, matrix permeability, fracture aperture, fracture permeability, overall thermal conductivity, and enthalpy of each phase) are updated, thus completing the first... Multi-field coupled crack evolution simulation at each time step.

[0129] Step S8: Time stepping and simulation termination.

[0130] From the The time step advances to the 1st Repeat steps S5 through S7 for each time step until the coarse-scale model simulation ends.

[0131] Application examples

[0132] To fully verify the superiority of the proposed method in terms of simulation accuracy and computational efficiency, this embodiment simultaneously establishes a fine-scale model, an uncorrected coarse-scale model, and the upscaling model proposed in this invention to simulate the Shell MDP-S oil shale in-situ conversion pilot test. The comparison results of the three models in terms of reaction kinetics, pore pressure, fracture initiation timing, rock matrix porosity-permeability evolution, and key indicators of oil and gas production are as follows:

[0133] (1) Reaction kinetics: such as Figure 3 As shown, in the grid near the heating well (e.g., the H1A grid), the fine-scale model reveals that the kerogen pyrolysis reaction rate exhibits multiple pulse-like fluctuations, reaching 6.37, 5.36, and 2.87 kg / (m³) on days 14, 79, and 103, respectively. 3 The local peak value (on day 53) was observed; however, the uncorrected coarse-scale model failed to capture the aforementioned transient response kinetics, exhibiting only a single artificial peak value of 14.69 kg / (m³) on day 53. 3 The upscaling model proposed in this invention accurately predicts the pulse timing and peak size of the fine-scale reaction rate through dynamic correction of the reaction frequency factor, maintaining a high degree of consistency with the fine-scale reference solution.

[0134] (2) Regarding the relationship between pore pressure and crack initiation timing: such as Figure 4 As shown, the uncorrected coarse-scale model, due to an artificially inflated reaction rate, caused the pore pressure near the heated well to rise to 5000–6000 kPa on day 48, resulting in premature tensile failure of the rock matrix and the premature initiation of 340 fractures. In contrast, the fine-scale model showed a pressure of only about 4000 kPa in the same area on day 48, with no fractures forming until day 141, when a fracture network of 5103 fractures developed. Figure 5As shown, the upscaling model proposed in this invention accurately reproduces the evolution of the fine-scale reaction rate and pore pressure, and the timing of thermally induced crack initiation is precisely delayed to day 140, which strictly matches the fine-scale reference timing, effectively eliminating the non-physical premature crack initiation phenomenon of the coarse-scale model.

[0135] (3) Regarding oil and gas production: such as Figure 6 As shown, the uncorrected coarse-scale model exhibits a non-physical premature rise in the daily oil production curve on day 48 due to premature fracturing and an artificially inflated reaction rate, with a peak value as high as 3.48m. 3 / day, with a peak daily gas production of 229.44m³. 3 / day, approximately the fine-scale result (86.22m) 3 The daily oil and gas production curves of the upscaling model proposed in this invention are in high agreement with the fine-scale reference results, and the cumulative oil production of the upscaling model is 115.98 m³ / day. 3 , compared with the fine-scale reference solution (116.53m) 3 The results were almost identical, with a cumulative gas production of 2.08 × 10⁻⁶. 4 m 3 , compared with the fine-scale reference solution (2.03×10 4 m 3 The error is less than 2.5%.

[0136] (4) Computational efficiency: On a 6-core Intel Xeon Platinum 8260 CPU (2.4 GHz) platform, the full physical fine-scale simulation takes approximately 1.16 × 10⁻⁶ seconds. 7 The uncorrected coarse-scale simulation takes approximately 7.58 × 10⁻⁶ seconds (about 134.65 days). 5 The simulation time for the scale improvement method described in this invention is approximately 3.46 × 10 seconds. 5 The upscaling simulation achieved a computational speedup of 33.50 times compared to the fine-scale simulation (approximately 4.01 days). At the same time, the upscaling simulation was faster than the uncorrected coarse-scale simulation. The main reason is that the excessively fast reaction rate and excessively high pore pressure in the coarse-scale model were effectively improved, which significantly improved the convergence and robustness of the numerical model.

[0137] In summary, the adaptive upscaling method for thermal flow solidification multi-field coupled fracture evolution proposed in this invention maintains a high degree of consistency with the fine-scale reference results in terms of reaction kinetics, timing of thermally induced fracture initiation, evolution of rock matrix porosity and permeability, and key indicators for oil and gas production prediction. At the same time, it achieves a computational speedup of approximately 33 times compared to full physical fine-scale simulation. It has significant theoretical and engineering application value for the design of field-scale development schemes, prediction of development indicators, and optimization of production dynamics in in-situ conversion of oil shale.

[0138] Of course, the above description is not intended to limit the present invention, and the present invention is not limited to the examples given above. Any changes, modifications, additions or substitutions made by those skilled in the art within the scope of the present invention should also fall within the protection scope of the present invention.

Claims

1. An adaptive upscaling method for multi-field coupled crack evolution in thermal flow solidification, characterized in that, Includes the following steps: S1. Considering the coupling relationship between natural fractures and heat transfer field, seepage field, stress field and chemical reaction field, a multi-field coupled fine-scale model of in-situ conversion heat flow solidification of oil shale is established. The mesh of the fine-scale model is uniformly coarsened to obtain the corresponding coarse-scale model, and the spatial mapping relationship between the coarse-scale mesh and its internal fine-scale sub-mesh is established. S2. Simplify the coarse-scale model and the fine-scale model respectively to obtain the simplified coarse-scale model and the simplified fine-scale model. Run the simplified coarse-scale model and the simplified fine-scale model respectively to obtain the grid temperature values ​​of each model at each time step. S3. Run the unsimplified coarse-scale model to obtain the temperature values ​​of each coarse-scale grid at each time step. Based on the temperature difference between the simplified coarse-scale model and the unsimplified coarse-scale model at each time step, correct the temperature of each grid in the simplified fine-scale model at each time step to obtain the predicted full physical fine-scale grid temperature. S4. Based on the corrected fine-scale grid temperature, iteratively calculate the reaction rate and kerogen concentration of each fine-scale grid at each time step, and then weight the reaction rate of all fine-scale sub-grids in each coarse-scale grid according to the pore volume to obtain the average reaction rate of the coarse-scale grid after scaling up. S5. Initialize the multiphysics parameters of the unsimplified coarse-scale model. At each time step, adaptively and dynamically correct the reaction frequency factor of the coarse-scale model based on the average reaction rate of the coarse-scale grid obtained in step S4. S6. Using the corrected coarse-scale model, solve the mass, energy and momentum conservation equations, and detect the key stress data of each rock matrix grid at each time step; Based on key stress data, determine the meshes that have experienced tensile or shear failure, identify the location of crack initiation, and construct new crack meshes based on an embedded discrete crack model. S7. When a new crack is generated at a time step, detect the intersection characteristics between the new crack and the existing crack, merge cracks with the same opening direction and adjacent positions, and update the multi-field physical property parameters. S8, from the first The time step advances to the 1st Repeat steps S5 through S7 for each time step until the coarse-scale model simulation ends.

2. The adaptive upscaling method for multi-field coupled crack evolution in thermal flow solidification as described in claim 1, characterized in that, In step S1, the coupling method for establishing the multi-field coupled model of in-situ thermal flow solidification of oil shale is as follows: For the mass conservation equation: ； In the formula, This is for the simulated duration; For fluid pore volume; The volume of a single matrix or crack mesh block; express Phase density; for Phase saturation; The components in the current grid exist The mass fraction in the phase; Components exist The diffusion coefficient in the phase; For the current grid and the first Components of adjacent grids exist The difference in mass fraction between the phases; For the current grid and the first The contact area between adjacent cracks or matrix grids; The geometric center of the current mesh and the first The projection length of the line connecting the geometric centers of two adjacent cracks or matrix meshes in the direction of the normal vector of the two mesh contact surfaces; This represents the number of matrix grids adjacent to the current grid. This represents the number of crack meshes adjacent to the current mesh. The total number of chemical reactions; For fluid pore volume; and Components In the Stoichiometric coefficients of products and reactants in a chemical reaction; For the first The volumetric reaction rate of a chemical reaction; for The source and sink of phases; For the current grid and the first Between adjacent cracks or matrix meshes The conductivity coefficient of the phase; Indicates the current grid and its neighboring grids The phase potential energy difference is calculated using the following formula: ； In the formula, For the current grid Phase pore pressure; For grid of Phase pore pressure; For the energy conservation equation: ； In the formula, It represents the total pore volume containing both fluid and chemically reactive solid components. The volume of an inert rock matrix that does not undergo chemical reactions; express The mass and internal energy of a phase; The internal energy of the mass of a solid phase that can undergo chemical reactions; The internal energy of the inert rock matrix is ​​mentioned; The mass concentration of the solid component that can undergo a chemical reaction; for Enthalpy of phase mass; The overall thermal conductivity of the saturated rock matrix; This represents the temperature difference between the current grid and its neighboring grids. For the first Mass enthalpy change of a chemical reaction; Energy input provided for heating the well; For the momentum conservation equation: ； In the formula, Poisson's ratio; Biot coefficient; This represents the total pore pressure of the current grid. It is the linear thermal expansion coefficient; It is the bulk modulus; It is a volume force; For temperature; For normal mean stress, The calculation method is as follows: ； in, , and These are the maximum principal stress, intermediate principal stress, and minimum principal stress, respectively.

3. The adaptive upscaling method for multi-field coupled crack evolution in thermal flow solidification as described in claim 1, characterized in that, Step S2 involves simplifying the coarse-scale and fine-scale models: The porosity, permeability, and initial saturation of each fluid component in both the coarse-scale and fine-scale models were set to zero, and the production wells were shut down, leaving only the electrically heated wells supplying heat to the reservoir.

4. The adaptive upscaling method for multi-field coupled crack evolution in thermal flow solidification as described in claim 1, characterized in that, Step S3 involves correcting and predicting the temperature of the simplified fine-scale grid. Step 3.1: Based on the temperatures of each mesh in the simplified coarse-scale model obtained in step S2 at each time step, and the temperatures of each mesh in the unsimplified coarse-scale model obtained in step S3 at each time step, calculate the temperature difference between the two meshes at each time step. For the... Coarse-scale grid: ； In the formula, For the unsimplified coarse-scale model and the simplified coarse-scale model in the 1st... The grid, the first Temperature difference at each time step; For the unsimplified coarse-scale model in the first The grid, the first Temperature at each time step; To simplify the coarse-scale model in the first The grid, the first Temperature at each time step; Step 3.2: Add the temperature difference obtained in Step 3.1 to the temperature of the corresponding mesh in the simplified fine-scale model to obtain the predicted temperature of the full physical fine-scale mesh. ； In the formula, For the prediction of the full physical fine scale The grid in the first Temperature at each time step; To simplify the fine-scale model The grid in the first Temperature at each time step.

5. The adaptive upscaling method for multi-field coupled crack evolution in thermal flow solidification as described in claim 1, characterized in that, In step S4, the process of iteratively calculating the reaction rate of each fine-scale grid and scaling up is performed: Step 4.1: Initialize the reaction kinetic parameters of each fine-scale mesh, including the reaction frequency factor, activation energy, initial kerogen concentration, and pore volume. Step 4.2, for the first At the nth time step, based on the corrected fine-scale mesh temperature, the Arrhenius reaction kinetic model is used to calculate the nth time step. The fine-scale grid in the first... Reaction rate at each time step: ； In the formula, For the first In the nth fine-scale grid The reaction rate at each time step; The frequency factor for fine-scale reactions remains constant throughout the simulation. Activation energy; It is the ideal gas constant; The proportion of fluids and chemically reactive solid components Total pore volume of a fine-scale grid; For the first The fine-scale grid in the first... The mass concentration of the solid-phase reactants at each time step; Step 4.3, update the following formula: Kerogen concentration at each time step in the fine-scale grid: ； In the formula, For the first The fine-scale grid in the first... The mass concentration of the solid-phase reactants at each time step; For the first The time step size of each time step; Step 4.4, from the first The time step advances to the 1st Repeat steps 4.2 to 4.3 for each time step until the simulation ends; Step 4.5, for each coarse-scale grid The reaction rates of all fine-scale subgrids in the middle are weighted by pore volume to obtain the first... All fine-scale grids corresponding to the coarse-scale grid in the first coarse-scale grid are in the second coarse-scale grid. Average reaction rate at each time step This is the reaction rate of the coarse-scale grid after scaling up.

6. The adaptive upscaling method for multi-field coupled crack evolution in thermal flow solidification as described in claim 1, characterized in that, Step S5 involves the adaptive dynamic correction of the response frequency factor of the coarse-scale model. Step 5.1, calculate the first coarse-scale model. The grid in the first Reaction rate at each time step: ； In the formula, For the first The coarse-scale grid in the first... The reaction rate at each time step; For the first The coarse-scale grid in the first... The uncorrected response frequency factor at each time step is equal to the value at the previous time step. Corrected response frequency factor ; For the first The coarse-scale grid in the first... The mass concentration of the solid-phase reactants at each time step; Step 5.2: Based on the coarse-scale grid average reaction rate obtained in step S4, calculate the... The coarse-scale grid in the first... Correction factor for response frequency factor at each time step: ； In the formula, For the first The coarse-scale grid in the first... The reaction frequency factor correction factor for each time step; The first one obtained in step S4 All fine-scale grids corresponding to the coarse-scale grid in the first coarse-scale grid are in the second coarse-scale grid. The average reaction rate at each time step; Step 5.3, correct the first according to the following formula. The coarse-scale grid in the first... Response frequency factor at each time step: ； In the formula, For the first The coarse-scale grid in the first... Corrected response frequency factor for each time step; Step 5.4: Substitute the corrected reaction frequency factor into the unsimplified coarse-scale model, and proceed to step S6 for the next step. Multi-field coupled crack evolution simulation at each time step.

7. The adaptive upscaling method for multi-field coupled crack evolution in thermal flow solidification as described in claim 1, characterized in that, In step S6, the process of identifying the rock matrix mesh that generates thermally induced fractures is as follows: If the stress condition of the rock matrix mesh satisfies the following formula, then the mesh is considered to have undergone tensile failure: ； In the formula, The tensile strength of the rock matrix, The minimum effective principal stress of the rock matrix; If the stress condition of the rock matrix mesh satisfies the following equation, then the mesh is considered to have undergone shear failure: ； In the formula, For shear stress, For cohesion, It is normal stress. It is the internal friction angle; Shear stress and normal stress The calculation formulas are as follows: ； ； In the formula, Rock failure surface and minimum total principal stress The angle between them.

Images