Oil shale pyrolysis reaction analysis method and system

Abstract

The present application provides a kind of oil shale pyrolysis reaction analysis method and system, it is related to oil shale pyrolysis simulation technical field, the present application is first based on the fractal porous medium transport theory calculation fractal tortuosity factor, to accurately depict micro-pore structure;Using weighted least squares inversion first-order reaction parameter vector, guarantee the accuracy of kinetic parameters;Combined with semi-coke cumulative mass to construct dynamic permeability model, quantifies the pore blockage effect caused by solid deposition;Established with heavy oil residence time as the core secondary cracking rate model, restores the secondary loss process of heavy oil in low-permeability environment;Using explicit time step algorithm step by step to solve mass conservation equation, realize heat-flow-hydrolysis multi-field coupling and the fidelity simulation of component evolution;Finally, based on product deterioration index determines the optimal fracturing time, provides scientific engineering decision basis for in-situ exploitation, reduces resource waste.

Description

Technical Field

[0001] This invention relates to the field of oil shale pyrolysis simulation technology, specifically to a method and system for analyzing oil shale pyrolysis reactions. Background Technology

[0002] Oil shale, as an important unconventional oil and gas resource, is of strategic significance for ensuring energy security through its efficient development. Currently, in-situ oil shale extraction technology has become the mainstream approach in the industry due to its minimal environmental disturbance and relatively controllable energy consumption. This technology involves heating underground rock formations and using thermochemical reactions to break down solid kerogen into oil and gas fluids. With the deep integration of industrial software and traditional mining, numerical simulation technology has become a core tool for assessing mine reserves, predicting production capacity, and optimizing heating schemes. However, the underground in-situ environment is extremely complex, involving strong coupling of thermal, fluid, and chemical fields. How to accurately recreate the evolution of microscopic pores and the laws governing phase transitions of matter using simulation methods, and achieve a leap from extensive heating to refined management, is a key direction for current industrial technology upgrading.

[0003] In the prior art, CN119227297A discloses a method and system for analyzing the pyrolysis reaction of oil shale based on kinetic model calibration. This technology focuses on solving the problem of accurately obtaining chemical kinetic parameters. Its core approach involves collecting thermogravimetric experimental data of oil shale and using optimization algorithms to optimize and calibrate parameters of multiple parallel reaction models, thereby retrieving the activation energy and pre-exponential factor that characterize the kerogen pyrolysis properties. This method effectively uses laboratory measured data to correct the deviations of the theoretical model, achieving numerical simulation of the primary pyrolysis rate of kerogen and the product formation law, representing the current typical technical status of using experimental data to drive numerical calculations in the industry.

[0004] However, the aforementioned existing technologies still have significant limitations in handling the coupling of microscopic transport and complex reactions. First, existing models often neglect the fractal characteristics of oil shale pore structure and the dynamic blockage effect of semi-coke deposition on fluid channels, leading to large deviations in the calculation of real-time permeability. Second, traditional methods mainly focus on the generation of primary products, failing to fully consider the retention effect of heavy oil in low-permeability environments, and cannot quantify the secondary cracking reaction caused by high-temperature retention, resulting in simulation results that often overestimate oil production. Furthermore, due to the lack of a real-time monitoring mechanism for the evolution trend of product components, existing technologies struggle to capture the critical point of product quality degradation and effectively reduce secondary cracking losses caused by heavy oil retention due to semi-coke blockage.

[0005] The information disclosed in the background section is only intended to enhance the understanding of the background of this disclosure, and therefore may include information that does not constitute prior art known to those skilled in the art. Summary of the Invention

[0006] The purpose of this invention is to provide a method and system for analyzing the pyrolysis reaction of oil shale, so as to solve the problems mentioned in the background art.

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

[0008] A method for analyzing the pyrolysis reaction of oil shale, comprising the following steps:

[0009] In the oil shale mining area, the area used for pyrolysis reaction analysis was selected as the target area. The physical property data of oil shale in the target area were collected, and the thermogravimetric time series data of oil shale in the target area were obtained. The fractal tortuosity factor was obtained by analyzing the physical property data using the fractal porous media transport theory.

[0010] Based on thermogravimetric time-series data, multiple parallel first-order reaction kinetic models corresponding to different kerogen components were constructed. Using thermogravimetric time-series data as observation constraints, the weighted least squares method was used to perform parameter inversion on the multiple parallel first-order reaction kinetic models to obtain the first-order reaction parameter vector.

[0011] The initial permeability is determined based on the fractal tortuosity factor combined with the fractal porous media fluid dynamics theory, and a dynamic permeability model is established for calculating the real-time permeability. The dynamic permeability model uses the cumulative mass of semi-coke as the pore blockage variable, and establishes a secondary cracking reaction rate model with heavy oil residence time as the cracking inducing factor and heavy oil concentration as the concentration correction variable.

[0012] During the pyrolysis reaction, a set of mass conservation differential equations containing kerogen molar concentration, heavy oil concentration and semi-coke cumulative mass is established based on the first-order reaction parameter vector. The set of mass conservation differential equations is solved step by step over time. At each time step, the real-time permeability is updated by calling the dynamic permeability model based on the current time step's semi-coke cumulative mass. The heavy oil residence time is calculated based on the updated real-time permeability.

[0013] Based on the heavy oil concentration and residence time of the current time step, the secondary cracking reaction rate model is called to update the secondary cracking rate, and the updated secondary cracking rate is introduced into the mass conservation equation to solve, thus obtaining the component evolution dataset updated by secondary cracking.

[0014] The product degradation index at each time step is calculated based on the component evolution dataset and compared with the preset degradation index threshold. Based on the comparison results, the time step when the product degradation index meets the preset degradation judgment condition is taken as the optimal fracturing time.

[0015] Furthermore, the physical properties of oil shale include: total organic carbon content, initial porosity, specific surface area, pore fractal dimension, pore connectivity coefficient, pore shape factor, and stress sensitivity coefficient.

[0016] The oil shale sample was heated at a preset constant heating rate until it reached the preset target temperature. The heating was then stopped. Temperature, time, and the real-time remaining mass of the sample were recorded at equal time intervals to obtain thermogravimetric time-series data and generate thermogravimetric time-series data.

[0017] Furthermore, the calculation principle of the fractal tortuosity factor of oil shale samples is as follows:

[0018]

[0019] in, The fractal tortuosity factor represents the oil shale sample. This represents the pore connectivity coefficient of an oil shale sample. This indicates the initial porosity of the oil shale sample. This represents the pore fractal dimension of the oil shale sample.

[0020] Furthermore, the central difference method was used to perform numerical differentiation on the thermogravimetric time series dataset to calculate the experimental conversion rate at the i-th time step.

[0021] The principle of calculating the experimental conversion rate at time i:

[0022]

[0023] in, Indicates the first The experimental conversion rate at each time point. This indicates the preset time interval. This represents the real-time remaining mass of the sample at time 1. Indicates the first The real-time remaining mass of the sample at each moment. Indicates the first The real-time remaining mass of the sample at each moment. Indicates the first The real-time remaining mass of the sample at each moment. Indicates the total number of recorded moments;

[0024] Based on chemical kinetics theory, a mathematical model is constructed containing a predetermined total number of parallel first-order reactions. A deviation objective function is established in conjunction with the experimental conversion rate. The optimal primary reaction parameters for each parallel first-order reaction are iteratively found using a global optimization algorithm.

[0025] The mathematical model for parallel first-order reactions is as follows:

[0026]

[0027] in, This indicates the total number of parallel reactions pre-set by relevant staff. Represents the pre-exponential factor of the j-th component. This represents the activation energy of the j-th component. This indicates the number of points calculated using the model. The overall cumulative conversion rate of parallel first-order reactions at each time point. This represents the natural exponential function. Represents the ideal gas constant. Indicates the first The absolute temperature at any given moment;

[0028] The expression for calculating the deviation objective function is as follows:

[0029]

[0030] in, Indicates the first The initial mass fraction of each component, and satisfying the condition from arrive The sum of the initial mass fractions of each group is 1. This represents the vector of first-order reaction parameters to be solved, specifically... ={ , ..., , ..., }, This represents the deviation function value calculated by substituting the first-order reaction parameter vector P into the equation. This indicates that the first one was preset by relevant staff. The weighting coefficient at each time point.

[0031] Furthermore, the initial permeability of the oil shale sample was calculated based on the fractal tortuosity factor and the fractal porous media fluid dynamics theory.

[0032] The principle of initial permeability calculation for oil shale samples:

[0033]

[0034] in, This indicates the initial permeability of the oil shale sample. This represents the pore shape factor of an oil shale sample. This indicates the specific surface area of ​​the oil shale sample.

[0035] Based on the initial permeability of the oil shale sample and combined with the accumulated mass of semi-coke, the real-time permeability of the current time step caused by semi-coke deposition is calculated using a dynamic permeability model.

[0036] The real-time penetration rate model for the current time step is as follows:

[0037]

[0038] in, This indicates that during the in-situ pyrolysis simulation of the target region, the first... Real-time penetration rate at any given moment. This indicates that during the in-situ pyrolysis simulation of the target region, the first... The real-time porosity at each moment was obtained through solid-phase volume equilibrium calculations based on the principle of conservation of mass during the pyrolysis simulation. Indicates the stress sensitivity coefficient. This indicates the pressure of the overlying strata in the preset target area. This indicates that during the in-situ pyrolysis simulation of the target region, the first... The real-time pore pressure at each moment was obtained by calculating using the law of conservation of fluid mass during the pyrolysis simulation. This represents the preset semi-coke clogging factor. This indicates that the target region to be input is in the in-situ pyrolysis simulation process, the first... The cumulative mass of half-focal length at a given moment. This represents the preset initial pore volume, which is obtained by calculating the product of the initial porosity of the oil shale sample and the volume of the target region.

[0039] Based on the heavy oil residence time at the current time step, the cracking induction probability factor is calculated using threshold activation theory and the Sigmoid activation function.

[0040] Principle of calculating the cleavage induction probability factor:

[0041]

[0042] in, This represents the probability factor for cleavage induction. This indicates the residence time of heavy oil at the current time step to be input. This represents the preset time-temperature threshold coefficient. This indicates that during the in-situ pyrolysis simulation of the target region, the first... The temperature at that moment This represents the preset probability smoothing factor;

[0043] Based on the cracking induction probability factor and the heavy oil concentration at the current time step, the secondary cracking rate at the current time step is calculated using the law of mass action and a probability correction to the Aronne equation.

[0044] The secondary fragmentation rate model for the current time step is as follows:

[0045]

[0046] in, This indicates that during the in-situ pyrolysis simulation of the target region, the first... The rate of secondary cleavage reaction at time 1. This represents the pre-split exponential factor of the second quadratic split. This represents the preset activation energy for secondary pyrolysis. This indicates that the target region to be input is in the in-situ pyrolysis simulation process, the first... The molar concentration of heavy oil at each moment.

[0047] Furthermore, at the initial moment of the simulation, the molar concentration of kerogen in each parallel first-order reaction is initialized based on the total organic carbon content and the oil shale density obtained through density logging.

[0048] The principle of initializing the molar concentration of kerogen in each parallel first-order reaction is as follows:

[0049]

[0050] in, This represents the initial molar concentration of kerogen in the first-order reaction of group j. Indicates the density of oil shale. Indicates the total organic carbon content. This indicates the preset molar mass of kerogen;

[0051] Using the first-order reaction parameter vector, the total rate of the first-order parallel reaction at the current time step is calculated by summation. Combined with the secondary cracking reaction rate, the generation and consumption terms of kerogen, heavy oil, gas and semi-coke components are constructed. Using a time-step-based discrete iterative algorithm, based on the concentration state and the total rate of the first-order parallel reaction at the current time step, the molar concentrations of kerogen, heavy oil and semi-coke after the pyrolysis reaction of oil shale at the current time are calculated according to the law of conservation of mass.

[0052] The principle for calculating the rate of first-order parallel reactions in group j:

[0053]

[0054] in, This indicates that during the in-situ pyrolysis simulation of the target region, the first... The first-order parallel reactions in group 1 were in the 1st... The reaction rate at each moment This indicates that during the in-situ pyrolysis simulation of the target region, the first... The first-order parallel reactions in group 1 were in the 1st... The molar concentration of kerogen at each time point;

[0055] Principle of updating molar concentration of kerogen components:

[0056]

[0057] in, This indicates that during the in-situ pyrolysis simulation of the target region, the first... The first-order parallel reactions in group 1 were in the 1st... The molar concentration of kerogen at each time point, This indicates the preset simulation time interval step size;

[0058] The principle of heavy oil component molar concentration update after first-order parallel reaction:

[0059]

[0060] in, This indicates that during the in-situ pyrolysis simulation of the target region, at the [missing information] th [missing information] The heavy oil molar concentration obtained after a first-order parallel reaction at a given time point. This indicates that during the in-situ pyrolysis simulation of the target region, at the [missing information] th [missing information] The molar concentration of heavy oil obtained after secondary cracking at a given time point, At that time, the first Molar concentration of heavy oil at time 1 , This represents the preset stoichiometric coefficient for kerogen oil production;

[0061] The principle of half-coke component molar concentration update after first-order parallel reaction:

[0062]

[0063] in, This indicates that during the in-situ pyrolysis simulation of the target region, at the [missing information] th [missing information] The half-joule molar concentration obtained at time point after a first-order parallel reaction. This indicates that during the in-situ pyrolysis simulation of the target region, at the [missing information] th [missing information] The half-joule molar concentration obtained after secondary pyrolysis at a given time point This represents the pre-defined stoichiometric coefficient for kerogen semi-coke production;

[0064] Semi-coke cumulative quality update principle:

[0065]

[0066] in, This indicates that during the in-situ pyrolysis simulation of the target region, at the [missing information] th [missing information] The cumulative mass of half-focal length at a given moment. This indicates the volume of the preset target area. This indicates the molar mass of a semi-joule;

[0067] Based on the in-situ pyrolysis simulation process of the target region, in the first... The cumulative mass of semi-coke at each time step is used to update the real-time permeability using the dynamic permeability model. Based on the updated real-time permeability, the residence time of heavy oil at the current time step is calculated.

[0068]

[0069] in, This indicates that during the in-situ pyrolysis simulation of the target region, the first... The residence time of heavy oil at a given moment Indicates the preset temperature The dynamic viscosity of the oil. This indicates the length of the preset target area.

[0070] Furthermore, based on the heavy oil concentration and residence time of the current time step, the secondary cracking reaction rate model is called to update the secondary cracking rate, and the updated secondary cracking rate is introduced into the mass conservation equation to solve, thus obtaining the component evolution dataset updated by secondary cracking.

[0071] The principle of updating the molar concentration of heavy oil components after secondary cracking:

[0072]

[0073] in, This indicates that during the in-situ pyrolysis simulation of the target region, at the [missing information] th [missing information] The heavy oil molar concentration obtained after secondary cracking at a certain time point. This indicates that during the in-situ pyrolysis simulation of the target region, at the [missing information] th [missing information] Heavy oil convection transport at a given moment This indicates the preset Darcy flow rate in the oil phase;

[0074] The principle of updating the molar concentration of half-coke components after secondary pyrolysis:

[0075]

[0076] in, This indicates that during the in-situ pyrolysis simulation of the target region, at the [missing information] th [missing information] The half-joule molar concentration obtained after secondary pyrolysis at time +1. This represents the preset stoichiometric coefficient for semi-coke produced by heavy oil cracking;

[0077] Principle of updating molar concentration of gas components:

[0078]

[0079] in, This indicates that during the in-situ pyrolysis simulation of the target region, at the [missing information] th [missing information] The gas molar concentration at time +1, This indicates that during the in-situ pyrolysis simulation of the target region, at the [missing information] th [missing information] The gas molar concentration at a given time. This represents the preset stoichiometric coefficient for kerogen gas production. This represents the preset stoichiometric coefficient for gas production from heavy oil cracking. This indicates that during the in-situ pyrolysis simulation of the target region, at the [missing information] th [missing information] Gas convection transport at a given moment, This indicates the preset gas-phase Darcy velocity.

[0080] Furthermore, based on the component evolution dataset, and combined with the preset target region volume and preset molar masses of heavy oil and gas, spatial integral statistics are performed on all target regions of the entire oil shale mining area to calculate the values ​​at time t and time t during the in-situ pyrolysis simulation of the entire oil shale mining area. The total cumulative mass of heavy oil and the total cumulative mass of gas at time +1;

[0081] Based on the calculations obtained during the in-situ pyrolysis simulation of the entire oil shale mining area, at time t and at... The total cumulative mass of heavy oil and the total cumulative mass of gas at time +1 are used to construct a formula for calculating the product degradation index and calculate the product degradation index.

[0082] The formula for calculating the product degradation index is as follows:

[0083]

[0084] in, During the in-situ pyrolysis simulation of the entire oil shale mining area, the product degradation index at time t is... This represents the total accumulated mass of gas in the entire oil shale mining area at time t during the in-situ pyrolysis simulation. This represents the total accumulated mass of heavy oil in the entire oil shale mining area at time t during the in-situ pyrolysis simulation. During the in-situ pyrolysis simulation of the entire oil shale mining area, the total mass of the accumulated gas at time t+1 is... This represents the total accumulated mass of heavy oil in the entire oil shale mining area at time t+1 during the in-situ pyrolysis simulation. This represents the preset regularization coefficient to prevent the denominator from being zero;

[0085] The product degradation index is compared with the preset degradation index threshold. If the product degradation index fails to exceed the preset degradation index threshold for three consecutive moments during the in-situ pyrolysis simulation of the entire oil shale mining area, it is considered that the current oil shale pyrolysis reaction does not require fracturing. If the product degradation index exceeds the preset degradation index threshold for three consecutive moments, it is considered that the current oil shale pyrolysis reaction requires fracturing, and the first time it exceeds the threshold among the three consecutive moments is taken as the optimal fracturing time.

[0086] An analysis system for the pyrolysis reaction of oil shale, the analysis system being used to implement the above-mentioned analysis method, comprising:

[0087] The data acquisition module is used to select the area for pyrolysis reaction analysis in the oil shale mining area as the target area, collect the physical property data of oil shale in the target area, obtain the thermogravimetric time series data of oil shale in the target area, and use the fractal porous media transport theory to analyze the physical property data to obtain the fractal tortuosity factor.

[0088] The first-order reaction inversion module is used to construct multiple parallel first-order reaction kinetic models corresponding to different kerogen components based on thermogravimetric time-series data. Using thermogravimetric time-series data as observation constraints, the weighted least squares method is used to perform parameter inversion on multiple parallel first-order reaction kinetic models to obtain the first-order reaction parameter vector.

[0089] The secondary cracking parameter construction module is used to determine the initial permeability based on the fractal tortuosity factor combined with the fractal porous media fluid dynamics theory, and to establish a dynamic permeability model for calculating the real-time permeability. The dynamic permeability model uses the cumulative mass of semi-coke as the pore blockage variable, and establishes a secondary cracking reaction rate model with heavy oil residence time as the cracking inducing factor and heavy oil concentration as the concentration correction variable.

[0090] The first-order reaction simulation module is used to establish a set of mass conservation differential equations based on the first-order reaction parameter vector during the pyrolysis reaction process. This set of equations includes the molar concentration of kerogen, the concentration of heavy oil, and the cumulative mass of semi-coke. The module solves the set of mass conservation differential equations step by step. At each time step, the module calls the dynamic permeability model to update the real-time permeability based on the cumulative mass of semi-coke at the current time step. The module then calculates the residence time of heavy oil based on the updated real-time permeability.

[0091] The secondary cracking update module is used to update the secondary cracking rate by calling the secondary cracking reaction rate model based on the heavy oil concentration and residence time of the current time step, and then introduce the updated secondary cracking rate into the mass conservation equation to solve, thereby obtaining the component evolution dataset after the secondary cracking update.

[0092] The optimal time determination module is used to calculate the product degradation index at each time step based on the component evolution dataset and compare it with the preset degradation index threshold. Based on the comparison results, the time step when the product degradation index meets the preset degradation determination condition is taken as the optimal fracturing time.

[0093] Compared with the prior art, the beneficial effects of the present invention are:

[0094] This invention calculates the fractal tortuosity factor based on the fractal porous media transport theory and uses a dynamic permeability model to calculate the real-time permeability caused by semi-coke deposition. By introducing the fractal tortuosity factor, the complex microscopic pore topology of oil shale can be more realistically depicted, overcoming the shortcomings of traditional methods that simplify the reservoir into a homogeneous model. More importantly, this invention uses the accumulated mass of semi-coke generated in real time during the chemical reaction as a key physical variable and feeds it back into the fluid dynamics calculation. This approach effectively quantifies the physical phenomena of decreased porosity and increased fluid flow resistance caused by semi-coke deposition as the pyrolysis reaction proceeds, ensuring the accuracy of the real-time permeability of the target area dynamically changing with semi-coke accumulation during the simulation. It also solves the problem of distortion in the velocity and pressure field calculations caused by the failure to consider the solid phase deposition blockage effect, significantly improving the simulation accuracy of the fluid transport capacity of complex underground porous media.

[0095] This invention also calculates the heavy oil residence time and the cracking induction probability factor to modify the Arunnes equation and obtain the secondary cracking rate, and determines the optimal fracturing time based on the product degradation index. This invention uses the heavy oil residence time to calculate the cracking induction probability factor, reconstructing the reverse loss process of heavy oil cracking into gas and semi-coke due to slow flow in a low-permeability environment, and correcting the production prediction bias caused by a single first-order reaction model. Based on this, this invention further calculates the product degradation index and compares it with a preset degradation index threshold. By monitoring the trend of this index, the critical time window for secondary cracking, where the product composition transitions from high oil production to high gas production, can be sensitively captured. This data-driven decision-making logic transforms complex component evolution data into optimal fracturing time instructions that directly guide engineering operations, providing a scientific basis for on-site intervention before losses escalate, and effectively reducing resource waste and reduced returns caused by overheating or excessively long heavy oil residence time. Attached Figure Description

[0096] Figure 1 This is a schematic diagram of the overall method flow of the present invention.

[0097] Figure 2 This is a schematic diagram of the overall system structure of the present invention. Detailed Implementation

[0098] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to specific embodiments.

[0099] It should be noted that, unless otherwise defined, the technical or scientific terms used in this invention should have the ordinary meaning understood by one of ordinary skill in the art to which this invention pertains. The terms "first," "second," and similar terms used in this invention do not indicate any order, quantity, or importance, but are merely used to distinguish different components. Terms such as "comprising" or "including" mean that the element or object preceding the word encompasses the elements or objects listed following the word and their equivalents, without excluding other elements or objects. Terms such as "connected" or "linked" are not limited to physical or mechanical connections, but can include electrical connections, whether direct or indirect. Terms such as "upper," "lower," "left," and "right" are used only to indicate relative positional relationships; when the absolute position of the described object changes, the relative positional relationship may also change accordingly.

[0100] Example:

[0101] Please see Figure 1 The present invention provides a technical solution:

[0102] A method for analyzing the pyrolysis reaction of oil shale, comprising the following steps:

[0103] Step 1: Select the area in the oil shale mining area for pyrolysis reaction analysis as the target area, collect the physical property data of oil shale in the target area, obtain the thermogravimetric time series data of oil shale in the target area, and use the fractal porous media transport theory to analyze the physical property data to obtain the fractal tortuosity factor.

[0104] In this embodiment, a region in the oil shale mining area selected for pyrolysis reaction analysis is designated as the target region. Oil shale samples are collected from each target region, and multiphysics experiments are performed on the collected oil shale samples. The multiphysics experiments include:

[0105] The total organic carbon content of the oil shale sample was determined using an organic carbon analyzer.

[0106] Thermogravimetric analysis was used to test the pyrolysis characteristics. The oil shale sample was heated at a preset constant heating rate. The temperature, time and real-time remaining mass of the sample were recorded according to the preset thermogravimetric analysis time window step to obtain thermogravimetric time series data and generate a thermogravimetric time series dataset.

[0107] Microscopic porosity tests were conducted using mercury intrusion porosimetry to determine the initial porosity and specific surface area of ​​oil shale samples. Simultaneously, the fractal dimension of the pores was calculated based on the slope of the double logarithmic curve of mercury intrusion pressure versus mercury intrusion saturation. The pore connectivity coefficient of the oil shale samples was determined using nuclear magnetic resonance technology.

[0108] Pyrolysis experiments were conducted on mined oil shale samples. The pore shape factor was obtained using image analysis. SEM images of the cross-section of the oil shale samples were taken. The pore edges were identified using image processing software. The perimeter and area of ​​the pores were calculated. The pore shape factor of the oil shale samples was calculated according to the shape factor calculation formula.

[0109] Mechanical properties were tested using a triaxial rock pressure experimental setup. The permeability of oil shale under different confining pressures was measured, and the stress sensitivity coefficient was obtained by fitting.

[0110] In-situ pyrolysis of oil shale is a dynamic process accompanied by phase transitions and property evolution. Total organic carbon (TOC) content is the material basis of the pyrolysis reaction, directly determining the initial reserves and potential oil production of kerogen per unit volume of rock. Thermogravimetric time-series data is the core characterization of reaction kinetics, recording the mass change trajectory of the sample over time at a specific heating rate. It is the sole data source for subsequent inversion of first-order reaction kinetic parameters, determining the accuracy of chemical field calculations. Initial porosity, specific surface area, and stress sensitivity coefficient constitute the physical boundaries of fluid dynamics, especially the stress sensitivity coefficient, which reflects the influence of high confining pressure at deep underground levels on pore closure. Obtaining these key parameters using an organic carbon analyzer, thermogravimetric analyzer, mercury intrusion porosimetry, and triaxial rock mechanics experimental setup can eliminate errors introduced by theoretical assumptions, ensuring a high degree of consistency between the simulation model and the physical reality of the target mining area.

[0111] The principle for calculating the fractal tortuosity factor of oil shale samples is as follows:

[0112]

[0113] in, The fractal tortuosity factor represents the oil shale sample. This represents the pore connectivity coefficient of an oil shale sample. This indicates the initial porosity of the oil shale sample. This represents the pore fractal dimension of the oil shale sample.

[0114] The fractal tortuosity factor combines several key microscopic physical factors, including pore connectivity coefficient, initial porosity, and pore fractal dimension. The resulting fractal tortuosity factor is a dimensionless physical quantity that characterizes the effect of microstructure on flow obstruction. It reflects the ratio of the actual flow path length of fluid molecules in the complex porous medium of oil shale to the macroscopic straight-line distance. By introducing the fractal tortuosity factor, the correction effect of the irregularity of the microscopic pore topology on macroscopic seepage capacity is quantified, overcoming the shortcomings of traditional models that simplify pores into regular channels, and significantly improving the physical accuracy of subsequent permeability calculations. Initial porosity directly affects the calculation result of the fractal tortuosity factor. A smaller initial porosity means a higher proportion of solid-phase skeleton within the rock, a narrower and more severely obstructed space for fluid flow, and the fluid must bypass more obstacles to pass through, resulting in a more tortuous and longer actual flow path. Therefore, a larger fractal tortuosity factor corresponds to a larger initial porosity; that is, initial porosity and the fractal tortuosity factor are negatively correlated. This formula also fully considers the topological complexity of the pore network. The pore connectivity coefficient reflects the complexity of the connections between pore channels. The larger the pore connectivity coefficient, the higher the probability that the fluid will detour or take a complex path during transmission, and therefore the larger the fractal tortuosity factor. In addition, the pore fractal dimension, as part of the exponential term, describes the roughness and space-filling characteristics of the pore surface. The larger the pore fractal dimension, the more complex the pore structure and the stronger the self-similarity, leading to a non-linear increase in the fractal tortuosity factor, thus more sensitively capturing the influence of microstructural changes on flow resistance.

[0115] Step 2: Construct multiple parallel first-order reaction kinetic models corresponding to different kerogen components based on thermogravimetric time-series data. Using thermogravimetric time-series data as observation constraints, the weighted least squares method is used to perform parameter inversion on multiple parallel first-order reaction kinetic models to obtain the first-order reaction parameter vector.

[0116] In this embodiment, the central difference method is used to perform numerical differentiation on the thermogravimetric time series dataset to calculate the experimental conversion rate at the i-th time step.

[0117] The principle of calculating the experimental conversion rate at time i:

[0118]

[0119] in, Indicates the first The experimental conversion rate at each time point. This indicates the preset time interval. This represents the real-time remaining mass of the sample at time 1. Indicates the first The real-time remaining mass of the sample at each moment. Indicates the first The real-time remaining mass of the sample at each moment. Indicates the first The real-time remaining mass of the sample at each moment. Indicates the total number of recorded moments;

[0120] No. The experimental conversion rate at time step 1 is a kinetic indicator reflecting the intensity of the kerogen pyrolysis reaction at that instant. Specifically, it reflects the rate at which the sample completes normalized mass conversion per unit time. Its technical advantage lies in smoothing experimental noise through central difference, providing a reaction rate curve with a high signal-to-noise ratio. This calculation is correlated with the sample mass change at adjacent time steps; the numerator in the formula essentially characterizes the mass change from time step 1 to time step 2. Time to the The proportion of mass loss at any given time relative to the total volatile mass. The experimental conversion rate at time step 1 is positively correlated with the mass difference between adjacent time steps. This means that within the same time step, the faster the sample mass decreases, the more rapidly the pyrolysis reaction proceeds, and the higher the calculated conversion rate. Meanwhile, the... The experimental conversion rate at each moment is negatively correlated with the preset thermogravimetric analysis time window step. With the time step as the denominator, the shorter the sampling interval, the faster the instantaneous reaction rate, given a constant mass change.

[0121] Based on chemical kinetics theory, a mathematical model is constructed containing a predetermined total number of parallel first-order reactions. A deviation objective function is established in conjunction with the experimental conversion rate. The optimal primary reaction parameters for each parallel first-order reaction are iteratively found using a global optimal algorithm.

[0122] The mathematical model for parallel first-order reactions is as follows:

[0123]

[0124] in, This indicates the total number of parallel reactions pre-set by relevant staff. Represents the pre-exponential factor of the j-th component. This represents the activation energy of the j-th component. This indicates the first time the model is calculated. The overall cumulative conversion rate of parallel first-order reactions at each time point. This represents the natural exponential function. Represents the ideal gas constant. Indicates the first The absolute temperature at any given moment;

[0125] The expression for calculating the deviation objective function is as follows:

[0126]

[0127] in, Indicates the first The initial mass fraction of each component, and satisfying the condition from arrive The sum of the initial mass fractions of each group is 1. This represents the vector of first-order reaction parameters to be solved, specifically... ={ , ..., , ..., }, This represents the deviation function value calculated by substituting the first-order reaction parameter vector P into the equation. This indicates that the first one was preset by relevant staff. Weighting coefficients at each time point;

[0128] Kerogen in oil shale is not a single component, but a mixture of various organic macromolecules with different bond energies. Traditional single-reaction models cannot simultaneously fit the main peak and shoulder peak characteristics in the pyrolysis curve. Therefore, this method constructs a mathematical model containing multiple parallel reactions, assuming that kerogen undergoes independent first-order decomposition by several virtual components with different activation energies. To find the optimal parameter combination for these virtual components, a bias objective function must be constructed. The significance of the bias objective function lies in quantifying the total error between theoretical predictions and experimental facts. By minimizing this error, the kinetic parameters that best match the characteristics of the oil shale in this mining area can be derived. The bias objective function value is a statistical index used to evaluate the accuracy of the kinetic parameters, specifically reflecting the overall deviation between the theoretical reaction rate curve generated by the currently set first-order reaction parameter vector and the rate curve measured in actual experiments. The technical advantage of the bias objective function value is that it transforms the complex curve fitting problem into a mathematical extremum problem, making it possible to automatically find the optimal parameters using a global optimum algorithm. The deviation objective function value is positively correlated with the difference between the experimental conversion rate and the theoretically calculated rate. This means that the closer the rate calculated by the theoretical model is to the experimentally measured value, the smaller the difference within the parentheses, and the smaller the final accumulated deviation objective function value, indicating that the parameters are more accurate. Simultaneously, this function is also affected by the... The adjustment of the weight coefficient at time step n, the first time step n The weighting coefficients at each time point are positively correlated with the bias contribution value. By increasing the weighting coefficients of the most intense reaction phases, the optimization algorithm can be forced to prioritize ensuring the fitting accuracy of the main reaction phase, thereby improving the reliability of the inversion results in engineering applications. Finally, by minimizing this function value, the system outputs a first-order reaction parameter vector containing the activation energies, pre-exponential factors, and initial mass fractions of each component, completing the transformation from macroscopic experimental data to microscopic simulation parameters.

[0129] Step 3: Determine the initial permeability based on the fractal tortuosity factor and the fractal porous media fluid dynamics theory, and establish a dynamic permeability model for calculating the real-time permeability. The dynamic permeability model uses the cumulative mass of semi-coke as the pore blockage variable, and establishes a secondary cracking reaction rate model with heavy oil residence time as the cracking inducing factor and heavy oil concentration as the concentration correction variable.

[0130] In this embodiment, the initial permeability of the oil shale sample is calculated based on the fractal tortuosity factor and the fractal porous media fluid dynamics theory.

[0131] The principle of initial permeability calculation for oil shale samples:

[0132]

[0133] in, This indicates the initial permeability of the oil shale sample. This represents the pore shape factor of an oil shale sample. This indicates the specific surface area of ​​the oil shale sample.

[0134] Initial permeability is a macroscopic physical quantity characterizing the fluid flow capacity of oil shale in its original state. Specifically, it reflects the ease with which fluid can pass through the rock medium before pyrolysis. Oil shale is dense and heterogeneous, and its initial flow capacity is constrained by its microscopic pore structure. Traditional empirical formulas often neglect the irregularity of pore shape and the tortuousness of the path, leading to distorted estimations of baseline permeability. This embodiment introduces the fractal tortuosity factor calculated in step 1 to derive the macroscopic initial permeability from the microscopic mechanism, thereby establishing a flow benchmark that incorporates microscopic topological features. This method fully considers the geometric resistance of fluid flow. The denominator in the formula includes a comprehensive consideration of specific surface area, pore shape factor, and fractal tortuosity factor, which also determines the influence relationship between the parameters: initial porosity is positively correlated with initial permeability. The larger the initial porosity, the wider the effective channel cross section for fluid flow and the smaller the flow resistance. On the other hand, the specific surface area of ​​oil shale samples is negatively correlated with initial permeability. The larger the specific surface area, the larger the friction surface in contact between the fluid and the solid skeleton, and the significantly increased viscous resistance. At the same time, the fractal tortuosity factor is negatively correlated with initial permeability. The greater the tortuosity, the longer the actual path of fluid flow and the greater the pressure loss along the flow path, thus leading to a decrease in permeability.

[0135] Based on the initial permeability of the oil shale sample, the accumulated mass of semi-coke at the current time moment is retrieved from step five, and the real-time permeability at the current time step caused by semi-coke deposition is calculated using the dynamic permeability model.

[0136] The principle of real-time penetration rate calculation at the current time step:

[0137]

[0138] in, This indicates that during the in-situ pyrolysis simulation of the target region, the first... Real-time penetration rate at any given moment. This indicates that during the in-situ pyrolysis simulation of the target region, the first... The real-time porosity at each moment was obtained through solid-phase volume equilibrium calculations based on the principle of conservation of mass during the pyrolysis simulation. This represents the stress sensitivity coefficient, with units of 1000 kJ / m². , The overlying strata pressure representing the preset target area is calculated using density logging and logging integration methods. This indicates that during the in-situ pyrolysis simulation of the target region, the first... The real-time pore pressure at each moment was obtained by calculating using the law of conservation of fluid mass during the pyrolysis simulation. This represents the preset semi-coke clogging factor, in units of... The semi-coke plugging factor is used to characterize the equivalent plugging intensity caused by a unit mass of semi-coke in a unit initial pore volume, leading to a decrease in the permeability of oil shale. The preset semi-coke plugging factor value is related to the pore throat radius distribution of oil shale and is calibrated through core flow experiments. Under fixed temperature and pressure conditions, semi-coke-containing fluid is injected into the core of collected oil shale samples, and the relationship between the decrease in permeability and the mass of sediment is measured. The larger the semi-coke plugging factor, the more sensitive the permeability of the oil shale in the target area is to semi-coke plugging. This indicates that the target region to be input is in the in-situ pyrolysis simulation process, the first... The cumulative mass of half-focal length at a given moment. This represents the preset initial pore volume, which is obtained by calculating the product of the initial porosity of the oil shale sample and the volume of the target region.

[0139] Real-time permeability is a dynamic variable coupled with geomechanical and chemical depositional effects, specifically reflecting the remaining fluid transport capacity at the current moment after being subjected to both effective stress compression and semi-coke blockage. During in-situ heating, the permeability of underground oil shale is not constant. The pressure of the overlying strata causes pore closure, and the semi-coke generated by the pyrolysis reaction deposits in the pore throats, causing physical blockage. Without calculating real-time permeability, the model cannot simulate the real physical phenomena of more intense reactions, more severe blockage, and slower flow rates. This formula achieves real-time feedback of the cumulative mass of semi-coke on the permeability of underground oil shale. In the formula's variable relationships, real-time porosity is positively correlated with real-time permeability; the expansion of pore space is beneficial to fluid flow. However, the difference between the overlying strata pressure and the real-time pore pressure is negatively correlated with real-time permeability. As the difference between the overlying strata pressure and the real-time pore pressure increases, the rock skeleton is compressed, pore throats close, leading to an exponential decrease in the permeability of underground oil shale. More importantly, the formula introduces a semi-coke blockage correction term. The cumulative mass of semi-coke at each moment is negatively correlated with the real-time permeability. As the pyrolysis reaction proceeds, the more semi-coke is generated, the larger the pore volume it occupies, and the stronger the physical blocking effect on the fluid channel, which directly leads to a linear decrease in the real-time permeability.

[0140] Based on the heavy oil residence time at the current time step, the cracking induction probability factor is calculated using threshold activation theory and the Sigmoid activation function.

[0141] Principle of calculating the cleavage induction probability factor:

[0142]

[0143] in, This represents the probability factor for cleavage induction. This indicates the residence time of heavy oil at the current time step to be input. This represents the preset time-temperature threshold coefficient. This indicates that during the in-situ pyrolysis simulation of the target region, the first... The temperature at that moment This represents the preset probability smoothing factor;

[0144] The cracking induction probability factor is a dimensionless correction coefficient based on threshold activation theory, ranging from 0 to 1. It specifically reflects the statistical probability of a secondary cracking reaction occurring in a given fluid element under specific temperature and flow rate conditions. Traditional chemical kinetic models typically assume that the reaction will occur at a fixed rate once the temperature reaches the reaction conditions, neglecting the inhibitory or promoting effect of fluid flow on the reaction. In in-situ oil shale extraction, if heavy oil has a short residence time, it may flow out of the high-temperature zone before cracking even at high temperatures. Conversely, if blockage leads to retention, cracking may occur even at low temperatures. This method introduces the cracking induction probability factor, essentially adding a flow-reaction coupling switch to the Aronunis equation, solving the problem of traditional models overestimating or underestimating the amount of secondary cracking in heterogeneous flow fields. This formula uses a modified form of the Sigmoid activation function, which profoundly reflects the nonlinear influence of various physical quantities on the possibility of cracking: the residence time of heavy oil is positively correlated with the cracking induction probability factor. The longer the residence time, the longer the heavy oil molecules are heated, and the higher the probability of bond breaking when crossing the energy barrier. When the residence time exceeds a certain threshold, the probability factor approaches 1. At the same time, temperature is also positively correlated with the cracking induction probability factor. Physically, this means that the higher the temperature, the shorter the residence time threshold required for secondary cracking. High temperature environment will significantly reduce the sensitivity of cracking reaction to time, making the reaction easier to induce. The probability smoothing factor controls the steepness of the probability change from 0 to 1, and is used to simulate the transition range characteristics of the reaction.

[0145] Based on the cracking induction probability factor and the heavy oil concentration at the current time step, the secondary cracking rate at the current time step is calculated using the law of mass action and a probability correction to the Aronne equation.

[0146] The secondary fragmentation rate model for the current time step is as follows:

[0147]

[0148] in, This indicates that during the in-situ pyrolysis simulation of the target region, the first... The rate of secondary cleavage reaction at time 1. The pre-split factor for secondary cracking was obtained by consulting literature on heavy oil pyrolysis kinetics. The preset activation energy for secondary cracking was obtained by consulting literature on heavy oil pyrolysis kinetics. This indicates that the target region to be input is in the in-situ pyrolysis simulation process, the first... The molar concentration of heavy oil at each moment;

[0149] The secondary cracking reaction rate is a kinetic physical quantity characterizing the speed of conversion of heavy oil into gas and semi-coke. Specifically, it reflects the rate of conversion during in-situ pyrolysis simulation in the target region. At a given moment, this represents the amount of material lost per unit volume of heavy oil per unit time due to degradation. The core objective of oil shale extraction is to obtain liquid oil, while secondary cracking is a reverse reaction that leads to decreased oil yield, increased gas production, and exacerbated semi-coke blockage. By modifying the standard kinetic equation using the probability factors calculated in the previous step, this formula can reconstruct the dynamic loss process of crude oil cracking into gas and coke, providing the most crucial source term data for subsequent calculations of product degradation indices and determination of fracturing timing. In the variable relationships of the formula, the molar concentration of heavy oil is positively correlated with the secondary cracking reaction rate, following the law of mass action: the higher the reactant concentration, the higher the molecular collision frequency, and the faster the reaction rate. Temperature is strongly positively correlated with the secondary cracking reaction rate. As the temperature increases, the exponential term increases sharply, indicating that high temperature is the dominant factor leading to the intensification of secondary cracking. Most importantly, the cracking induction probability factor, as a multiplicative correction term, is positively correlated with the secondary cracking reaction rate. It plays a role in regulating the actual effectiveness. Even if the temperature and the molar concentration of heavy oil are both very high, if the fluid flow is extremely fast, causing the cracking induction probability factor to approach 0, the final calculated effective cracking rate will be very small, thus conforming to the physical fact of rapid displacement to protect oil products.

[0150] Step 4: During the pyrolysis reaction, a set of mass conservation differential equations containing kerogen molar concentration, heavy oil concentration and semi-coke cumulative mass is established based on the first-order reaction parameter vector. The set of mass conservation differential equations is solved step by step over time. At each time step, the real-time permeability is updated by calling the dynamic permeability model based on the current time step's semi-coke cumulative mass. The heavy oil residence time is calculated based on the updated real-time permeability.

[0151] In this embodiment, at the initial moment of the simulation, the molar concentration of kerogen in each parallel first-order reaction is initialized based on the total organic carbon content and the oil shale density obtained through density logging.

[0152] The principle of initializing the molar concentration of kerogen in each parallel first-order reaction is as follows:

[0153]

[0154] in, This represents the initial molar concentration of kerogen in the first-order reaction of group j. Indicates the density of oil shale. Indicates the total organic carbon content. This indicates the preset molar mass of kerogen;

[0155] This calculation realizes the mapping from geological parameters to chemical parameters. Its calculation principle is based on the definition of the amount of substance. By introducing the density of oil shale and the total organic carbon content, the total mass of organic matter contained in a unit volume of rock is calculated. Then, combined with the initial mass fraction of each component obtained from the inversion in step 2, the total mass is allocated to each parallel virtual reaction channel, thereby determining the molar concentration of kerogen of each component reaction at the initial time of the simulation.

[0156] Using the first-order reaction parameter vector, the total rate of the first-order parallel reaction at the current time step is calculated by summation. Combined with the secondary cracking reaction rate, the generation and consumption terms of kerogen, heavy oil, gas and semi-coke components are constructed. Using a time-step-based discrete iterative algorithm, based on the concentration state and the total rate of the first-order parallel reaction at the current time step, the molar concentrations of kerogen, heavy oil and semi-coke after the pyrolysis reaction of oil shale at the current time are calculated according to the law of conservation of mass.

[0157] The principle for calculating the rate of first-order parallel reactions in group j:

[0158]

[0159] in, This indicates that during the in-situ pyrolysis simulation of the target region, the first... The first-order parallel reactions in group 1 were in the 1st... The reaction rate at each moment This indicates that during the in-situ pyrolysis simulation of the target region, the first... The first-order parallel reactions in group 1 were in the 1st... The molar concentration of kerogen at each time point;

[0160] Principle of updating molar concentration of kerogen components:

[0161]

[0162] in, This indicates that during the in-situ pyrolysis simulation of the target region, the first... The first-order parallel reactions in group 1 were in the 1st... The molar concentration of kerogen at each time point, This indicates the preset simulation time interval step size;

[0163] The updated calculations of the first-order parallel reaction rate and kerogen component molar concentration in group j are based on the Arrhenius law and the explicit Euler integral method. The calculation principle lies in utilizing the in-situ pyrolysis simulation of the target region, the... The temperature at the moment and the moment The first-order parallel reactions in group 1 were in the 1st... The real-time reaction rate is calculated from the molar concentration of kerogen at each time step, and then discretely integrated over time. Since kerogen, as a solid reactant, is consumed only by unidirectional pyrolysis and does not involve flow migration, its molar concentration at the next time step depends only on the current concentration minus the amount consumed during that time step, reflecting the source-sink term handling mechanism in mass conservation.

[0164] The principle of heavy oil component molar concentration update after first-order parallel reaction:

[0165]

[0166] in, This indicates that during the in-situ pyrolysis simulation of the target region, at the [missing information] th [missing information] The heavy oil molar concentration obtained after a first-order parallel reaction at a given time point. This indicates that during the in-situ pyrolysis simulation of the target region, at the [missing information] th [missing information] The molar concentration of heavy oil obtained after secondary cracking at a given time point, At that time, the first Molar concentration of heavy oil at time 1 , The preset stoichiometric coefficients for kerogen oil production were obtained by consulting existing literature in petroleum geochemistry.

[0167] The principle of half-coke component molar concentration update after first-order parallel reaction:

[0168]

[0169] in, This indicates that during the in-situ pyrolysis simulation of the target region, at the [missing information] th [missing information] The half-joule molar concentration obtained at time point after a first-order parallel reaction. This indicates that during the in-situ pyrolysis simulation of the target region, at the [missing information] th [missing information] The half-joule molar concentration obtained after secondary pyrolysis at a given time point The stoichiometric coefficients for semi-coke production from kerogen are obtained by consulting existing literature in petroleum geochemistry.

[0170] Semi-coke cumulative quality update principle:

[0171]

[0172] in, This indicates that during the in-situ pyrolysis simulation of the target region, at the [missing information] th [missing information] The cumulative mass of half-focal length at a given moment. This indicates the volume of the preset target area. This indicates the molar mass of a semi-joule;

[0173] As a solid product, semi-coke's evolution logic focuses on in-situ accumulation. The calculation principle summarizes all solid residues. More importantly, by converting the updated semi-coke molar concentration into the cumulative semi-coke mass, this cumulative mass is fed back into the dynamic permeability model in real time, ensuring that the real-time permeability calculation at the next moment accurately reflects the pore blockage effect caused by the formation of semi-coke at the current moment.

[0174] Based on the in-situ pyrolysis simulation process of the target region, in the first... The cumulative mass of semi-coke at each time step is used to update the real-time permeability using the dynamic permeability model. Based on the updated real-time permeability, the residence time of heavy oil at the current time step is calculated.

[0175]

[0176] in, This indicates that during the in-situ pyrolysis simulation of the target region, the first... The residence time of heavy oil at a given moment Indicates the preset temperature The dynamic viscosity of the oil was obtained by referring to the standard heavy oil viscosity-temperature curve in the Petroleum Engineering Handbook. This indicates the length of the preset target area. This indicates that during the in-situ pyrolysis simulation of the neighboring regions of the preset target area, the first... The real-time pore pressure at a given moment was obtained through fluid mass conservation calculations during the pyrolysis simulation, based on the law of fluid mass conservation. This value was taken from the value obtained during the in-situ pyrolysis simulation of the target region at the given moment. Real-time pore pressure at a given moment The neighboring region with the largest difference in the in-situ pyrolysis simulation process is the first... Real-time pore pressure at any given moment;

[0177] Heavy oil residence time is a bridging parameter between fluid dynamics and chemical reaction kinetics, specifically reflecting the average physical time required for heavy oil particles to traverse a target region. If the generated heavy oil cannot be discharged in time, it will undergo secondary cracking at high temperatures, resulting in losses. Traditional Darcy's law can only calculate flow velocity, not directly quantify how long the oil remains in the target region. This method calculates heavy oil residence time by converting fluid flow rate into a chemical reaction time window, providing a direct temporal basis for subsequent calculations of the probability of secondary cracking. The formula reflects the competitive relationship between driving force and resistance. It presupposes a positive correlation between the dynamic viscosity of heavy oil at different temperatures and its residence time; the more viscous the heavy oil, the slower the flow, the longer the residence time, and the higher the risk of secondary cracking. Conversely, the dynamic viscosity of heavy oil at different temperatures is positively correlated with its residence time. The real-time permeability at any given moment is negatively correlated with the residence time of heavy oil. The higher the permeability, the smoother the fluid discharge and the shorter the residence time. At the same time, the real-time pore pressure difference between the neighboring region and the target region is also negatively correlated with the residence time of heavy oil. The greater the pressure difference, the stronger the driving force, and the faster the heavy oil can be discharged from the current target region, thereby reducing the residence time and the probability of being cracked and lost.

[0178] Step 5: Based on the heavy oil concentration and residence time of the current time step, call the secondary cracking reaction rate model to update the secondary cracking rate, and introduce the updated secondary cracking rate into the mass conservation equation to solve, so as to obtain the component evolution dataset updated by secondary cracking.

[0179] In this embodiment, based on the heavy oil concentration and heavy oil residence time at the current time step, the secondary cracking reaction rate model is called to update the secondary cracking rate, and the updated secondary cracking rate is introduced into the mass conservation equation to solve, thereby obtaining the component evolution dataset updated by secondary cracking.

[0180] The principle of updating the molar concentration of heavy oil components after secondary cracking:

[0181]

[0182] in, This indicates that during the in-situ pyrolysis simulation of the target region, at the [missing information] th [missing information] The heavy oil molar concentration obtained after secondary cracking at a certain time point. This indicates that during the in-situ pyrolysis simulation of the target region, at the [missing information] th [missing information] Heavy oil convection transport at a given moment This represents the preset Darcy velocity in the oil phase, calculated using Darcy's law for multiphase flow.

[0183] The principle of updating the molar concentration of half-coke components after secondary pyrolysis:

[0184]

[0185] in, This indicates that during the in-situ pyrolysis simulation of the target region, at the [missing information] th [missing information] The half-joule molar concentration obtained after secondary pyrolysis at time +1. This represents the preset stoichiometric coefficient for semi-coke produced by heavy oil cracking, which is obtained by consulting existing chemical handbooks related to heavy oil pyrolysis or petroleum hydrocarbon cracking reactions.

[0186] Principle of updating molar concentration of gas components:

[0187]

[0188] in, This indicates that during the in-situ pyrolysis simulation of the target region, at the [missing information] th [missing information] The gas molar concentration at time +1, This indicates that during the in-situ pyrolysis simulation of the target region, at the [missing information] th [missing information] The gas molar concentration at a given time. This represents the pre-defined stoichiometric coefficient for kerogen gas production, obtained by consulting existing literature in petroleum geochemistry. This represents the preset stoichiometric coefficient for gas production from heavy oil cracking, obtained by consulting existing chemical engineering handbooks related to heavy oil pyrolysis or petroleum hydrocarbon cracking reactions. This indicates that during the in-situ pyrolysis simulation of the target region, at the [missing information] th [missing information] Gas convection transport at a given moment, This represents the preset gas-phase Darcy velocity, calculated using Darcy's law for multiphase flow.

[0189] The calculation of heavy oil component molar concentration and gas component molar concentration is the core of this invention for achieving multiphysics coupling. Its calculation principle strictly follows the convection-reaction equation: for heavy oil, the increase in heavy oil component molar concentration originates from the generation term of primary kerogen cracking, minus the consumption term of secondary cracking, and is superimposed with the net convection flux driven by Darcy flow. For gas, the increase in gas component molar concentration originates simultaneously from primary kerogen cracking and secondary cracking of heavy oil, clearly quantifying the contribution of oil cracking to gas conversion. Through this iterative process, the complex behavior of heavy oil in porous media during pyrolysis—characterized by simultaneous generation, flow, and cracking—can be dynamically simulated. By using the latest heavy oil concentration after calculating convection transport and secondary cracking consumption as the heavy oil concentration for the next first-order parallel reaction, it ensures that the heavy oil molar concentration obtained after the first-order parallel reaction in the target region at the next moment during the in-situ pyrolysis simulation is based on the actual remaining oil volume within the target region. This avoids erroneously calculating high cracking amounts in areas where heavy oil has already flowed away, thereby improving the realism of the secondary cracking rate and enhancing the physical fidelity of the simulation model's prediction of cracking losses.

[0190] Step 6: Calculate the product degradation index based on the component evolution dataset and compare it with the preset degradation index threshold. According to the comparison results, the time when the product degradation index meets the preset degradation judgment condition is taken as the optimal fracturing time.

[0191] In this embodiment, based on the component evolution dataset, combined with the preset target region volume and preset molar masses of heavy oil and gas, spatial integral statistics are performed on all target regions of the entire oil shale mining area to calculate the values ​​at time t and time t during the in-situ pyrolysis simulation of the entire oil shale mining area. The total cumulative mass of heavy oil and the total cumulative mass of gas at time +1 are obtained. The molar mass of the gas is calculated by performing component analysis of the generated gas using a gas chromatograph during thermogravimetric analysis in step one, and by weighting the molar fraction of the gas components and their known molar masses.

[0192] Based on the calculations obtained during the in-situ pyrolysis simulation of the entire oil shale mining area, at time t and at... The total cumulative mass of heavy oil and the total cumulative mass of gas at time +1 are used to construct a formula for calculating the product degradation index and calculate the product degradation index.

[0193] The formula for calculating the product degradation index is as follows:

[0194]

[0195] in, During the in-situ pyrolysis simulation of the entire oil shale mining area, the product degradation index at time t is... This represents the total accumulated mass of gas in the entire oil shale mining area at time t during the in-situ pyrolysis simulation. This represents the total accumulated mass of heavy oil in the entire oil shale mining area at time t during the in-situ pyrolysis simulation. During the in-situ pyrolysis simulation of the entire oil shale mining area, the total mass of the accumulated gas at time t+1 is... This represents the total accumulated mass of heavy oil in the entire oil shale mining area at time t+1 during the in-situ pyrolysis simulation. This represents the preset regularization coefficient to prevent the denominator from being zero;

[0196] The product degradation index is a characteristic physical quantity reflecting the dynamic evolution trend of the quality of oil shale pyrolysis products over time. Specifically, it reflects the severity of the transition from liquid-dominated to gas-dominated product structure at the current time step. The core economic value of in-situ oil shale mining lies in obtaining liquid heavy oil. However, in the later stages of high-temperature extraction, the retained heavy oil undergoes secondary cracking, transforming into gas and semi-coke, leading to a decrease in oil production and resource waste. Simply observing cumulative production cannot accurately detect the point of qualitative change. This method, by constructing a product degradation index, can quantify the degradation rate of product components, thereby accurately capturing the critical time window for large-scale secondary cracking and providing a scientific basis for engineering intervention. In the formula's construction logic, the ratio of the total cumulative mass of gas to the total cumulative mass of heavy oil is used to characterize the relative abundance of gas and oil. The total accumulated mass of gas in the entire field is positively correlated with the product degradation index. As the pyrolysis reaction of oil shale proceeds, if heavy oil decreases due to cracking, this ratio increases, leading to a higher degradation index, indicating that product quality is rapidly deteriorating. Conversely, the total accumulated mass of heavy oil in the entire field, as the denominator of the ratio, is negatively correlated with the product degradation index. The more heavy oil accumulates, the lower the gas-oil ratio, and the smaller the product degradation index, indicating that the system is still in a healthy peak oil production period. Furthermore, the formula calculates the change in this ratio between adjacent time points and divides it by the time step, converting the difference in accumulation into a degradation rate. This makes the index more sensitive to sudden secondary cracking reactions. Simultaneously, by dividing by the sum of the total mass of gas and heavy oil, the dimensional differences caused by different total reserves in the mining area are eliminated, ensuring the universality of this index across different mining blocks.

[0197] The product degradation index is compared with the preset degradation index threshold. If the product degradation index fails to exceed the preset degradation index threshold for three consecutive moments during the in-situ pyrolysis simulation of the entire oil shale mining area, it is considered that the current oil shale pyrolysis reaction does not require fracturing. If the product degradation index exceeds the preset degradation index threshold for three consecutive moments, it is considered that the current oil shale pyrolysis reaction requires fracturing, and the first time it exceeds the threshold among the three consecutive moments is taken as the optimal fracturing time.

[0198] Through the above steps, this invention transforms complex thermo-fluidochemical coupling simulation data into a single decision command that directly guides field operations. Compared to the traditional, crude method of determining fracturing time based solely on experience or temperature field distribution, this invention's product degradation index-based judgment mechanism possesses extremely high sensitivity and scientific rigor. It not only identifies the inflection point where oil production efficiency begins to decline but also effectively filters out accidental fluctuations in numerical calculations through a continuous three-time-period judgment mechanism, ensuring the stability of the decision. This refined control method allows engineers to implement fracturing at the optimal time before significant secondary cracking losses occur in heavy oil, maximizing the preservation of liquid oil resources and reducing energy waste and economic losses caused by overheating or delayed extraction.

[0199] Please see Figure 2 The present invention also provides an oil shale pyrolysis reaction analysis system, the analysis system being used to implement the above-mentioned analysis method, comprising:

[0200] The data acquisition module is used to select the area for pyrolysis reaction analysis in the oil shale mining area as the target area, collect the physical property data of oil shale in the target area, obtain the thermogravimetric time series data of oil shale in the target area, and use the fractal porous media transport theory to analyze the physical property data to obtain the fractal tortuosity factor.

[0201] The first-order reaction inversion module is used to construct multiple parallel first-order reaction kinetic models corresponding to different kerogen components based on thermogravimetric time-series data. Using thermogravimetric time-series data as observation constraints, the weighted least squares method is used to perform parameter inversion on multiple parallel first-order reaction kinetic models to obtain the first-order reaction parameter vector.

[0202] The secondary cracking parameter construction module is used to determine the initial permeability based on the fractal tortuosity factor combined with the fractal porous media fluid dynamics theory, and to establish a dynamic permeability model for calculating the real-time permeability. The dynamic permeability model uses the cumulative mass of semi-coke as the pore blockage variable, and establishes a secondary cracking reaction rate model with heavy oil residence time as the cracking inducing factor and heavy oil concentration as the concentration correction variable.

[0203] The first-order reaction simulation module is used to establish a set of mass conservation differential equations based on the first-order reaction parameter vector during the pyrolysis reaction process. This set of equations includes the molar concentration of kerogen, the concentration of heavy oil, and the cumulative mass of semi-coke. The module solves the set of mass conservation differential equations step by step. At each time step, the module calls the dynamic permeability model to update the real-time permeability based on the cumulative mass of semi-coke at the current time step. The module then calculates the residence time of heavy oil based on the updated real-time permeability.

[0204] The secondary cracking update module is used to update the secondary cracking rate by calling the secondary cracking reaction rate model based on the heavy oil concentration and residence time of the current time step, and then introduce the updated secondary cracking rate into the mass conservation equation to solve, thereby obtaining the component evolution dataset after the secondary cracking update.

[0205] The optimal time determination module is used to calculate the product degradation index at each time step based on the component evolution dataset and compare it with the preset degradation index threshold. Based on the comparison results, the time step when the product degradation index meets the preset degradation determination condition is taken as the optimal fracturing time.

[0206] The above formulas are all dimensionless calculations. The formulas are derived from software simulations based on a large amount of collected data to obtain the most recent real-world results. The preset parameters in the formulas are set by those skilled in the art according to the actual situation.

[0207] The above embodiments can be implemented, in whole or in part, by software, hardware, firmware, or any other combination thereof. When implemented in software, the above embodiments can be implemented, in whole or in part, as a computer program product. Those skilled in the art will recognize that the units and algorithm steps of the various examples described in conjunction with the embodiments disclosed herein can be implemented by electronic hardware, or a combination of computer software and electronic hardware. Whether these functions are implemented in hardware or software depends on the specific application and design constraints of the technical solution.

[0208] The units described as separate components may or may not be physically separate. The components shown as units may or may not be physical units; they may be located in one place or distributed across multiple network units. Some or all of the units can be selected to achieve the purpose of this embodiment according to actual needs.

[0209] The above description is merely a specific embodiment of this application, but the scope of protection of this application is not limited thereto. Any changes or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in this application should be included within the scope of protection of this application.

Claims

1. A method for analyzing the pyrolysis reaction of oil shale, characterized in that, The specific steps include: In the oil shale mining area, the area used for pyrolysis reaction analysis was selected as the target area. The physical property data of oil shale in the target area were collected, and the thermogravimetric time series data of oil shale in the target area were obtained. The fractal tortuosity factor was obtained by analyzing the physical property data using the fractal porous media transport theory. Based on thermogravimetric time-series data, multiple parallel first-order reaction kinetic models corresponding to different kerogen components were constructed. Using thermogravimetric time-series data as observation constraints, the weighted least squares method was used to perform parameter inversion on the multiple parallel first-order reaction kinetic models to obtain the first-order reaction parameter vector. The initial permeability is determined based on the fractal tortuosity factor combined with the fractal porous media fluid dynamics theory, and a dynamic permeability model is established for calculating the real-time permeability. The dynamic permeability model uses the cumulative mass of semi-coke as the pore blockage variable, and establishes a secondary cracking reaction rate model with heavy oil residence time as the cracking inducing factor and heavy oil concentration as the concentration correction variable. During the pyrolysis reaction, a set of mass conservation differential equations containing kerogen molar concentration, heavy oil concentration and semi-coke cumulative mass is established based on the first-order reaction parameter vector. The set of mass conservation differential equations is solved step by step over time. At each time step, the real-time permeability is updated by calling the dynamic permeability model based on the current time step's semi-coke cumulative mass. The heavy oil residence time is calculated based on the updated real-time permeability. Based on the heavy oil concentration and residence time of the current time step, the secondary cracking reaction rate model is called to update the secondary cracking rate, and the updated secondary cracking rate is introduced into the mass conservation equation to solve, thus obtaining the component evolution dataset updated by secondary cracking. The product degradation index at each time step is calculated based on the component evolution dataset and compared with the preset degradation index threshold. Based on the comparison results, the time step when the product degradation index meets the preset degradation judgment condition is taken as the optimal fracturing time. The initial permeability of the oil shale sample was calculated based on the fractal tortuosity factor and the fractal porous media fluid dynamics theory. The principle of initial permeability calculation for oil shale samples: in, This indicates the initial permeability of the oil shale sample. This represents the pore shape factor of an oil shale sample. This indicates the specific surface area of ​​the oil shale sample. Based on the initial permeability of the oil shale sample and combined with the accumulated mass of semi-coke, the real-time permeability of the current time step caused by semi-coke deposition is calculated using a dynamic permeability model. The dynamic penetration rate model for the current time step is as follows: in, This indicates that during the in-situ pyrolysis simulation of the target region, the first... Real-time penetration rate at any given moment. This indicates that during the in-situ pyrolysis simulation of the target region, the first... The real-time porosity at each moment was obtained through solid-phase volume equilibrium calculations based on the principle of conservation of mass during the pyrolysis simulation. Indicates the stress sensitivity coefficient. This indicates the pressure of the overlying strata in the preset target area. This indicates that during the in-situ pyrolysis simulation of the target region, the first... The real-time pore pressure at each moment was obtained by calculating using the law of conservation of fluid mass during the pyrolysis simulation. This represents the preset semi-coke clogging factor. This indicates that the target region to be input is in the in-situ pyrolysis simulation process, the first... The cumulative mass of half-focal length at a given moment. This represents the preset initial pore volume, which is obtained by calculating the product of the initial porosity of the oil shale sample and the volume of the target region. Based on the heavy oil residence time at the current time step, the cracking induction probability factor is calculated using threshold activation theory and the Sigmoid activation function. Principle of calculating the cleavage induction probability factor: in, This represents the probability factor for cleavage induction. This indicates the residence time of heavy oil at the current time step to be input. This represents the preset time-temperature threshold coefficient. This indicates that during the in-situ pyrolysis simulation of the target region, the first... The temperature at that moment This represents the preset probability smoothing factor; Based on the cracking induction probability factor and the heavy oil concentration at the current time step, the secondary cracking reaction rate at the current time step is calculated using the law of mass action and a probability correction to the Aronne equation. The rate model for the secondary cleavage reaction at the current time step is as follows: in, This indicates that during the in-situ pyrolysis simulation of the target region, the first... The rate of secondary cleavage reaction at time 1. This represents the pre-split exponential factor of the second quadratic split. This represents the preset activation energy for secondary pyrolysis. This indicates that the target region to be input is in the in-situ pyrolysis simulation process, the first... The molar concentration of heavy oil at each moment.

2. The method for analyzing the pyrolysis reaction of oil shale according to claim 1, characterized in that: The physical properties of oil shale include: total organic carbon content, initial porosity, specific surface area, pore fractal dimension, pore connectivity coefficient, pore shape factor, and stress sensitivity coefficient. The oil shale sample was heated at a preset constant heating rate until it reached the preset target temperature. The heating was then stopped. Temperature, time, and the real-time remaining mass of the sample were recorded at equal time intervals to obtain thermogravimetric time-series data and generate thermogravimetric time-series data.

3. The method for analyzing the pyrolysis reaction of oil shale according to claim 2, characterized in that: The principle for calculating the fractal tortuosity factor of oil shale samples is as follows: in, The fractal tortuosity factor represents the oil shale sample. This represents the pore connectivity coefficient of an oil shale sample. This indicates the initial porosity of the oil shale sample. This represents the pore fractal dimension of the oil shale sample.

4. The method for analyzing the pyrolysis reaction of oil shale according to claim 3, characterized in that: The experimental conversion rate at time i was calculated by numerical differentiation using the central difference method on the thermogravimetric time series dataset. The principle of calculating the experimental conversion rate at time i: in, Indicates the first The experimental conversion rate at each time point. This indicates the preset time interval. This represents the real-time remaining mass of the sample at time 1. Indicates the first The real-time remaining mass of the sample at each moment. Indicates the first The real-time remaining mass of the sample at each moment. Indicates the first The real-time remaining mass of the sample at each moment. Indicates the total number of recorded moments; Based on chemical kinetics theory, a mathematical model is constructed containing a predetermined total number of parallel first-order reactions. A deviation objective function is established in conjunction with the experimental conversion rate. The optimal primary reaction parameters for each parallel first-order reaction are iteratively found using a global optimization algorithm. The mathematical model for parallel first-order reactions is as follows: in, This indicates the total number of parallel reactions pre-set by relevant staff. Represents the pre-exponential factor of the j-th component. This represents the activation energy of the j-th component. This indicates the first time the model is calculated. The overall cumulative conversion rate of parallel first-order reactions at each time point. This represents the natural exponential function. Represents the ideal gas constant. Indicates the first The absolute temperature at any given moment; The expression for calculating the deviation objective function is as follows: in, Indicates the first The initial mass fraction of each component, and satisfying the condition from arrive The sum of the initial mass fractions of each group is 1. This represents the vector of first-order reaction parameters to be solved, specifically... ={ , ..., , ..., }, This represents the deviation function value calculated by substituting the first-order reaction parameter vector P into the equation. This indicates that the first one was preset by relevant staff. The weighting coefficient at each time point.

5. The method for analyzing the pyrolysis reaction of oil shale according to claim 4, characterized in that: At the initial moment of the simulation, the molar concentration of kerogen in each parallel first-order reaction is initialized based on the total organic carbon content and the oil shale density obtained through density logging. The principle of initializing the molar concentration of kerogen in each parallel first-order reaction is as follows: in, This represents the initial molar concentration of kerogen in the first-order reaction of group j. Indicates the density of oil shale. Indicates the total organic carbon content. This indicates the preset molar mass of kerogen; Using the first-order reaction parameter vector, the total rate of the first-order parallel reaction at the current time step is calculated by summation. Combined with the secondary cracking reaction rate, the generation and consumption terms of kerogen, heavy oil, gas and semi-coke components are constructed. Using a time-step-based discrete iterative algorithm, based on the concentration state and the total rate of the first-order parallel reaction at the current time step, the molar concentrations of kerogen, heavy oil and semi-coke after the pyrolysis reaction of oil shale at the current time are calculated according to the law of conservation of mass. The principle for calculating the rate of first-order parallel reactions in group j: in, This indicates that during the in-situ pyrolysis simulation of the target region, the first... The first-order parallel reactions in group 1 were in the 1st... The reaction rate at each moment This indicates that during the in-situ pyrolysis simulation of the target region, the first... The first-order parallel reactions in group 1 were in the 1st... The molar concentration of kerogen at each time point; Principle of updating molar concentration of kerogen components: in, This indicates that during the in-situ pyrolysis simulation of the target region, the first... The first-order parallel reactions in group 1 were in the 1st... The molar concentration of kerogen at each time point, This indicates the preset simulation time interval step size; The principle of heavy oil component molar concentration update after first-order parallel reaction: in, This indicates that during the in-situ pyrolysis simulation of the target region, at the [missing information] th [missing information] The heavy oil molar concentration obtained after a first-order parallel reaction at a given time point. This indicates that during the in-situ pyrolysis simulation of the target region, at the [missing information] th [missing information] The molar concentration of heavy oil obtained after secondary cracking at a given time point, At that time, the first Molar concentration of heavy oil at time 1 , This represents the preset stoichiometric coefficient for kerogen oil production; The principle of half-coke component molar concentration update after first-order parallel reaction: in, This indicates that during the in-situ pyrolysis simulation of the target region, at the [missing information] th [missing information] The half-joule molar concentration obtained at time point after a first-order parallel reaction. This indicates that during the in-situ pyrolysis simulation of the target region, at the [missing information] th [missing information] The half-joule molar concentration obtained after secondary pyrolysis at a given time point This represents the pre-defined stoichiometric coefficient for kerogen semi-coke production; Semi-coke cumulative quality update principle: in, This indicates that during the in-situ pyrolysis simulation of the target region, at the [missing information] th [missing information] The cumulative mass of half-focal length at a given moment. This indicates the volume of the preset target area. This indicates the molar mass of a semi-joule; Based on the in-situ pyrolysis simulation process of the target region, in the first... The cumulative mass of semi-coke at each time step is used to update the real-time permeability using the dynamic permeability model. Based on the updated real-time permeability, the residence time of heavy oil at the current time step is calculated. in, This indicates that during the in-situ pyrolysis simulation of the target region, the first... The residence time of heavy oil at a given moment Indicates the preset temperature The dynamic viscosity of the oil. This indicates the length of the preset target area.

6. The method for analyzing the pyrolysis reaction of oil shale according to claim 5, characterized in that: Based on the heavy oil concentration and residence time of the current time step, the secondary cracking reaction rate model is called to update the secondary cracking rate, and the updated secondary cracking rate is introduced into the mass conservation equation to solve, thus obtaining the component evolution dataset updated by secondary cracking. The principle of updating the molar concentration of heavy oil components after secondary cracking: in, This indicates that during the in-situ pyrolysis simulation of the target region, at the [missing information] th [missing information] The heavy oil molar concentration obtained after secondary cracking at a certain time point. This indicates that during the in-situ pyrolysis simulation of the target region, at the [missing information] th [missing information] Heavy oil convection transport at a given moment This indicates the preset Darcy flow rate in the oil phase; The principle of updating the molar concentration of half-coke components after secondary pyrolysis: in, This indicates that during the in-situ pyrolysis simulation of the target region, at the [missing information] th [missing information] The half-joule molar concentration obtained after secondary pyrolysis at a certain time point This represents the preset stoichiometric coefficient for semi-coke produced by heavy oil cracking; Principle of updating molar concentration of gas components: in, This indicates that during the in-situ pyrolysis simulation of the target region, at the [missing information] th [missing information] The gas molar concentration at time +1, This indicates that during the in-situ pyrolysis simulation of the target region, at the [missing information] th [missing information] The gas molar concentration at a given time. This represents the preset stoichiometric coefficient for kerogen gas production. This represents the preset stoichiometric coefficient for gas production from heavy oil cracking. This indicates that during the in-situ pyrolysis simulation of the target region, at the [missing information] th [missing information] Gas convection transport at a given moment, This indicates the preset gas-phase Darcy velocity.

7. The method for analyzing the pyrolysis reaction of oil shale according to claim 6, characterized in that: Based on the component evolution dataset, and combined with the preset target region volume and preset molar masses of heavy oil and gas, spatial integral statistics are performed on all target regions of the entire oil shale mining area to calculate the values ​​at time t and time t during the in-situ pyrolysis simulation of the entire oil shale mining area. The total cumulative mass of heavy oil and the total cumulative mass of gas at time +1; Based on the calculations obtained during the in-situ pyrolysis simulation of the entire oil shale mining area, at time t and at... The total cumulative mass of heavy oil and the total cumulative mass of gas at time +1 are used to construct a formula for calculating the product degradation index and calculate the product degradation index. The formula for calculating the product degradation index is as follows: in, During the in-situ pyrolysis simulation of the entire oil shale mining area, the product degradation index at time t is... This represents the total accumulated mass of gas in the entire oil shale mining area at time t during the in-situ pyrolysis simulation. This represents the total accumulated mass of heavy oil in the entire oil shale mining area at time t during the in-situ pyrolysis simulation. During the in-situ pyrolysis simulation of the entire oil shale mining area, the total mass of the accumulated gas at time t+1 is... This represents the total accumulated mass of heavy oil in the entire oil shale mining area at time t+1 during the in-situ pyrolysis simulation. This represents the preset regularization coefficient to prevent the denominator from being zero; The product degradation index is compared with the preset degradation index threshold. If the product degradation index fails to exceed the preset degradation index threshold for three consecutive moments during the in-situ pyrolysis simulation of the entire oil shale mining area, it is considered that the current oil shale pyrolysis reaction does not require fracturing. If the product degradation index exceeds the preset degradation index threshold for three consecutive moments, it is considered that the current oil shale pyrolysis reaction requires fracturing, and the first time it exceeds the threshold among the three consecutive moments is taken as the optimal fracturing time.

8. A system for analyzing the pyrolysis reaction of oil shale, characterized in that: The analysis system is used to implement the analysis method according to any one of claims 1-7, including: The data acquisition module is used to select the area for pyrolysis reaction analysis in the oil shale mining area as the target area, collect the physical property data of oil shale in the target area, obtain the thermogravimetric time series data of oil shale in the target area, and use the fractal porous media transport theory to analyze the physical property data to obtain the fractal tortuosity factor. The first-order reaction inversion module is used to construct multiple parallel first-order reaction kinetic models corresponding to different kerogen components based on thermogravimetric time-series data. Using thermogravimetric time-series data as observation constraints, the weighted least squares method is used to perform parameter inversion on multiple parallel first-order reaction kinetic models to obtain the first-order reaction parameter vector. The secondary cracking parameter construction module is used to determine the initial permeability based on the fractal tortuosity factor combined with the fractal porous media fluid dynamics theory, and to establish a dynamic permeability model for calculating the real-time permeability. The dynamic permeability model uses the cumulative mass of semi-coke as the pore blockage variable, and establishes a secondary cracking reaction rate model with heavy oil residence time as the cracking inducing factor and heavy oil concentration as the concentration correction variable. The first-order reaction simulation module is used to establish a set of mass conservation differential equations based on the first-order reaction parameter vector during the pyrolysis reaction process. This set of equations includes the molar concentration of kerogen, the concentration of heavy oil, and the cumulative mass of semi-coke. The module solves the set of mass conservation differential equations step by step. At each time step, the module calls the dynamic permeability model to update the real-time permeability based on the cumulative mass of semi-coke at the current time step. The module then calculates the residence time of heavy oil based on the updated real-time permeability. The secondary cracking update module is used to update the secondary cracking rate by calling the secondary cracking reaction rate model based on the heavy oil concentration and residence time of the current time step, and then introduce the updated secondary cracking rate into the mass conservation equation to solve, thereby obtaining the component evolution dataset after the secondary cracking update. The optimal time determination module is used to calculate the product degradation index at each time step based on the component evolution dataset and compare it with the preset degradation index threshold. Based on the comparison results, the time step when the product degradation index meets the preset degradation determination condition is taken as the optimal fracturing time.

Images