A simulation method for differentiating mass transfer of multiphase multi-component fluid based on LBM

Abstract

This invention belongs to the technical field of numerical simulation of multiphase and multicomponent flow in porous media, specifically relating to a simulation method for differentiated mass transfer in multiphase and multicomponent fluids based on the LBM (Liquidity-Based Model). The method includes constructing a multicomponent thermodynamic system; correcting the critical parameters of each component using nanoconfining effects; and iteratively calculating and solving the multicomponent phase equilibrium using the PR (Proportional Relationship) equation of state. K The method includes: estimating the minimum miscibility pressure; calculating the miscibility factor; establishing a time-varying linkage mechanism between reservoir temperature and fluid properties; systematically calibrating LBM mesoscopic simulation parameters; constructing a multiphase, multicomponent flow simulation framework and embedding a differentiated mass transfer model; and quantitatively decomposing and verifying the closure error of four mass transfer mechanisms: expansion, immiscible interface extraction, miscible extraction, and displacement-carryover. This method can be used for the simulation and evaluation of multiphase, multicomponent mass transfer processes in porous media such as gas injection displacement and carbon dioxide geological storage in oil and gas reservoirs.

Description

Technical Field

[0001] This invention belongs to the technical field of numerical simulation of multiphase and multicomponent flow in porous media, and specifically relates to a simulation method for differentiated mass transfer of multiphase and multicomponent fluids based on LBM. Background Technology

[0002] Multiphase and multicomponent flow and mass transfer behavior in porous media are widely observed in engineering fields such as energy, environment, and chemical engineering. Taking CO2 displacement exploitation in oil and gas reservoirs as an example, after CO2 is injected into the reservoir, complex mass transfer processes occur between CO2 and various hydrocarbon components in the formation crude oil—including the dissolution of CO2 in the crude oil, crude oil volume expansion, selective extraction of light hydrocarbon components into the gas phase by CO2, and limited extraction of heavy hydrocarbon components at the interface. Accurately simulating and predicting the differentiated mass transfer behavior of each component and the relative contribution of each mass transfer mechanism during this process is of great guiding significance for optimizing injection parameters and improving recovery rate.

[0003] Lattice Boltzmann Method (LBM), as a mesoscale numerical simulation method, has become an important tool for simulating multiphase flow in porous media due to its advantages such as natural adaptability to complex pore boundaries, clear physical picture for easy handling of multiphase and multicomponent interactions, and high parallel computing efficiency. Within the LBM framework, multiphase and multicomponent flow is typically simulated using the Shan-Chen pseudopotential model to achieve phase separation or mixing between different fluid components. Its core idea is to characterize the repulsion or attraction between components at the mesoscale level by defining the pseudopotential function (effective density function) of each fluid component and the interaction strength parameter (G parameter) between component pairs, thereby naturally evolving the generation, deformation, and disappearance of phase interfaces (miscibility).

[0004] However, existing LBM-based multi-component mass transfer simulation methods suffer from the following three key technical shortcomings: First, there is a lack of a systematic calibration method for experimentally measurable LBM simulation parameters. Existing LBM multiphase flow simulations typically rely on trial and error to set key parameters such as G-parameters and pseudopotential characteristic densities, lacking a systematic calibration process to establish a quantitative correlation between these parameters and macroscopic physical properties that can be obtained experimentally, such as actual reservoir temperature, pressure, and fluid composition. Although the Peng-Robinson (PR) equation of state can accurately predict the phase equilibrium behavior of multi-component mixtures (including the phase equilibrium constant K-values ​​of each component and the gas-liquid phase composition distribution) based on basic physical properties such as the critical temperature, critical pressure, and eccentricity factor of the fluid, existing technologies have not yet provided a complete and traceable parameter transfer and calibration method for systematically mapping the thermodynamic calculation results of PR-EOS, especially the K-value gradient reflecting component volatility and the miscibility factor reflecting the degree of miscibility, to the G-parameters controlling phase separation / mixing behavior and the characteristic density controlling the pseudopotential shape in LBM.

[0005] Second, the thermodynamic effects of nanoscale pore confinement are not reflected in LBM parameter calibration. In tight reservoirs, the characteristic pore scale can be as low as 10 nm. At this scale, the confinement of fluid molecules by the pore walls causes a significant shift in the critical temperature and critical pressure compared to their bulk values ​​(i.e., the nanoscale confinement effect). Existing LBM models typically use bulk critical parameters directly when performing equation of state calculations, without correcting for the shift in critical parameters caused by nanoscale confinement. This results in LBM parameters calibrated based on the equation of state failing to accurately reflect the actual phase characteristics of the fluid in nanoporous systems.

[0006] Third, there is insufficient characterization of the differentiated mass transfer behavior of multiple components, and a lack of quantitative decomposition methods for mass transfer mechanisms. In actual crude oil systems, hydrocarbon components with different carbon numbers exhibit drastically different mass transfer behaviors during CO2 displacement due to significant differences in molecular weight, molecular size, volatility, and interaction strength with CO2: methane (C1) and light hydrocarbons (C2-C6) are readily extracted into the gas phase by CO2, while heavy hydrocarbons (C4-C6) are more readily extracted into the gas phase. 13+ Almost all components remain in the liquid phase. However, existing LBM mass transfer models typically apply the same mass transfer rules to all components, failing to reflect this differentiated mass transfer characteristic determined by thermodynamics. Furthermore, existing methods can only output overall indicators such as total oil recovery, and cannot quantitatively decompose total oil recovery into contributions from independent physical mechanisms such as displacement and carryover, dissolution and swelling, immiscible interface extraction, and miscible extraction.

[0007] Fourth, the simulation parameters lack a time-varying linkage mechanism that reflects changes in reservoir conditions (temperature, pressure). In actual reservoir development, reservoir temperature T is the core operating parameter affecting fluid properties. Increased temperature will simultaneously lead to a decrease in CO2 density, a decrease in CO2 viscosity, a decrease in the Henry's solubility constant (decreased solubility), an increase in mass transfer rate (Arrhenius effect), a decrease in the viscosity of light components of crude oil, a significant decrease in the viscosity of heavy components, an increase in the diffusion coefficient, and changes in binary interaction parameters. Temperature-driven correction is necessary. However, in existing LBM multiphase mass transfer simulation methods, the aforementioned temperature-driven fluid properties are typically fixed to calibration values ​​under a single operating condition, failing to automatically update in response to changes in operating conditions. When evaluating mass transfer behavior under different temperature conditions (e.g., 60°C to 150°C), researchers must recalibrate all LBM parameters for each temperature condition, lacking both physical consistency and the ability to support comparative analysis of multi-condition systems. Furthermore, during the simulation of the LBM main cycle, phase equilibrium... K Value field, Shan-Chen interaction strength parameter G, and miscibility factor All parameters evolve dynamically with the CO2 displacement process, but existing methods have failed to establish a complete mechanism for systematically updating these parameters over time with the simulation time step, resulting in simulation results that cannot accurately reflect the dynamic time-varying characteristics of phase evolution during CO2 displacement.

[0008] Therefore, there is an urgent need for a complete simulation framework that can start from the basic physical property parameters of reservoir fluids obtained from microfluidic experiments, systematically determine the LBM mesoscopic simulation parameters through nano-confined correction of equation of state thermodynamic calculations, construct a mass transfer model that reflects the differentiated characteristics of components, establish a time-varying linkage mechanism between reservoir temperature and other operating parameters and fluid properties, support multi-operating-condition system comparison, and realize quantitative decomposition and self-verification of mass transfer mechanism. Summary of the Invention

[0009] The purpose of this invention is to provide a simulation method for differential mass transfer in multiphase and multicomponent fluids based on Boltzmann Model (LBM) to address the four deficiencies mentioned above in existing technologies. This simulation method starts with reservoir fluid properties obtained from microfluidic experiments, calibrates lattice Boltzmann mesoscopic simulation parameters through thermodynamic calculations of the equation of state corrected by nanoconfining, and then constructs a multicomponent differential mass transfer model, achieving a complete simulation and analysis method that quantitatively decomposes the mass transfer mechanism and characterizes the time-varying properties of the parameters. It can be used for the simulation and evaluation of multiphase and multicomponent mass transfer processes in porous media such as gas injection displacement and carbon dioxide geological storage in oil and gas reservoirs.

[0010] This method is based on experimentally obtained basic physical property data such as reservoir temperature, crude oil carbon number distribution, density, viscosity and pore size. The critical parameters of each component phase are differentiated according to the nanoconfining ratio, and then substituted into the Peng-Robinson equation of state to solve for the multi-component phase equilibrium K value and miscibility factor. A time-varying linkage mechanism for parameters such as CO2 physical properties, Henry's dissolution constant, mass transfer rate, and crude oil viscosity was established with reservoir temperature as the driving variable. By comparing the three calibration paths of temperature, miscibility factor, and viscosity, the thermodynamic calculation results are systematically mapped to LBM mesoscopic parameters such as pseudopotential characteristic density, Shan-Chen interaction strength parameters, and relaxation time. We constructed an MRT-MCMP lattice Boltzmann simulation framework, embedding four differentiated mass transfer mechanisms: dissolution-desorption, fractional extraction, limited extraction of heavy components, and solvent effect-driven desorption, and introduced partitioned phase equilibrium dynamic feedback coupling. The contributions of four mass transfer mechanisms—expansion, immiscible interface extraction, miscible extraction, and displacement-carryover—were quantitatively decomposed and self-verified using mass conservation closure error.

[0011] The technical solution of this invention is: a simulation method for differential mass transfer of multiphase and multicomponent fluids based on LBM, comprising the following steps: (1) Constructing a multi-component thermodynamic system based on microfluidic experimental data: S1. Obtain basic physical property data of the target reservoir through microfluidic experiments and crude oil component testing; including: reservoir temperature T, initial reservoir pressure, etc. P reservoir 、 Characteristic pore size of reservoir rocks Crude oil density Crude oil viscosity Microfluidic calibration of diffusion coefficient, dynamic expansion factor, and CO2-crude oil minimum miscibility pressure. And the carbon number distribution of crude oil in reservoir fluids.

[0012] S2. Based on the carbon number distribution data of the crude oil in the reservoir fluid obtained, the crude oil is divided into different pseudo-components, including methane, light hydrocarbons, medium hydrocarbons and heavy hydrocarbons.

[0013] Furthermore, the typical classification of the pseudo-components is as follows: methane (C1, dissolved gas), light hydrocarbons (C2-C6, represented by pentane), and medium hydrocarbons (C7-C6). 12 (Represented by octane) and heavy hydrocarbons C 13+ (Represented by hexadecane).

[0014] S3. First, the mass fraction of each pseudo-component is obtained by summing the carbon number distributions; Then, the mole fraction of each pseudo-component is determined based on its mass fraction. Conversion: pseudo-component mole fraction = ; Finally, the obtained mole fractions of each pseudo-component are... Normalization yields the mole fractions of each hydrocarbon.

[0015] S4. Set the feed mole fraction for each component: First, inject fluid CO2, and calculate the CO2 feed mole fraction based on the amount of CO2 injected. ; Then, multiply the mole fractions of each hydrocarbon obtained in step S3 by (1- The feed mole fraction of each component was obtained.

[0016] Meanwhile, a dissolved CO2 phase is added to the original component system to independently track the concentration of CO2 that has dissolved into the liquid phase, serving as the state variable for calculating the mass transfer driving force and mass transfer of subsequent mechanism M1 (dissolution and desorption mass transfer of CO2 in the liquid phase).

[0017] In the Lattice Boltzmann model, the dissolved CO2 phase consists of two synergistic distribution functions: a density field distribution function, which participates in the calculation of inter-component interaction forces and contributes to liquid phase volume expansion and viscosity changes; and a concentration field distribution function, which does not participate in the calculation of interaction forces but is specifically used to provide a clean, locally dissolved CO2 concentration signal unaffected by interfacial numerical noise, for use in determining mass transfer driving forces. The specific construction methods of these two functions, their introduction into the LBM model, and their synergistic effects are detailed in step (5), mechanism M1.

[0018] (2) Based on the nanoconfining effect, the critical parameters of each component were corrected, and the phase equilibrium calculation of the PR equation of state was performed: S1. Obtain the critical temperature of the injected fluid and each pseudo-component in step (1) under bulk conditions. Critical pressure Eccentricity factor and molecular dynamics diameter .

[0019] S2. Correction of critical parameters for nano-confined execution: First, based on the characteristic pore size of the target reservoir rock and the molecular dynamic diameter of each component Differential determination of the number of adsorption layers of each component on the pore wall It follows the physical law that the larger the molecular size, the more adsorption layers there are.

[0020] Then, the adsorption layer thickness of each component is calculated according to formula (1). t ads，i and effective aperture d eff,i : (1); In the formula, t ads，i —Components i The thickness of the adsorption layer; n ads，i —Components i The number of adsorption layers; s i —Components i The molecular dynamics diameter; d eff,i —Components i Effective aperture; —Characteristic pore size of the target reservoir rock.

[0021] And calculate the confinement ratio of each component. x i , .

[0022] Finally, the nanoconfined critical parameter is corrected according to formula (2): (2); In the formula, —Components i The offset of the critical temperature relative to the bulk value; —Components i The offset of the critical pressure relative to the bulk value; —Components i Critical temperature under bulk conditions; —Components i Critical pressure under bulk conditions; f ( x i —Regarding the limit ratio x i Critical parameter offset correction function; a 1. a 2—calibration coefficient; where, a The typical value range for 1 is between -0.95 and -0.90. a The typical value of 2 is between 0.20 and 0.25; the calibration coefficient can be obtained through nanofluidic experiments or molecular dynamics simulation data. x i —Components i The limit ratio.

[0023] Finally, the corrected nanocritical parameters , As shown in formula (3): (3); In the formula, —Components i Critical temperature under confined conditions; —Components i Critical pressure under confined conditions.

[0024] Eccentricity factor It remains unchanged under nanoscale confinement conditions.

[0025] Due to the molecular dynamic diameter of each component Different and due to the number of adsorption layers Differentiated settings result in effective aperture They are also different, therefore the confinement ratios of each component are different. The components are different from each other, thus naturally achieving differentiated correction of the critical parameter shifts of each component. The larger the molecular size, the more adsorption layers, and the smaller the effective pore size of the component, the greater the shift in its critical parameter.

[0026] S3, PR-EOS phase equilibrium calculation: First, the nano-corrected critical temperature obtained in step S2 is... and nano-corrected critical pressure Substitute the parameters into the Peng-Robinson (PR) equation of state to calculate the phase thermodynamic parameters of each component.

[0027] ; In the formula, P —System pressure, v —molar volume R —Universal gas constant, T —Reservoir temperature, , Components i The gravitational parameters (including temperature correction) and co-volume parameters.

[0028] For components i The component parameters of the PR equation of state are calculated according to formula (4): (4); In the formula, —Base value of gravitational parameter of component i (without temperature correction, calculated only from nano-corrected critical parameter); b i —Co-volume parameter of component i; T r，i —Components i The relative temperatures; T —Reservoir temperature; —Components i Nanoscale modified critical temperature; —Components i Nanoscale modified critical pressure; R —Universal gas constant, R =8.314 J / (mol·K).

[0029] Temperature correction factor (T) Calculate according to formula (5): (5); In the formula, —Dimensionless temperature correction factor for component i; k i —Components i Temperature correction factor for the PR equation of state; T r，i —Components i The comparison temperature.

[0030] in, and The result calculated according to formula (6): (6).

[0031] The two empirical correlations in the above formula (6) are applicable to components with eccentricity factors not exceeding and exceeding 0.49, respectively, to take into account the fitting accuracy when the eccentricity factor of the heavy component is large.

[0032] Components i The temperature-corrected gravitational parameters are .

[0033] For multi-component mixtures, the classical mixing rules of van der Waals are used, and the gravitational parameters of the mixture are calculated according to formula (7). a mix Harmony volume parameters b mix .

[0034] (7); In the formula, a mix —The gravitational parameters of the mixture; b mix —Co-volume parameters of the mixture; —Componentsi mole fraction in the mixture; x j —Components j mole fraction in the mixture; i , j These are two different components within the same phase.

[0035] a i ( T ) — Components i Temperature-corrected gravitational parameters; a j ( T ) — Components j Temperature-corrected gravitational parameters; —Components i With components j Binary interaction parameters between them; b i —Components i The co-volume parameters.

[0036] in, The value of follows this rule: the greater the difference in carbon number between components, the higher the value of . The larger the value, the better.

[0037] For CO2-hydrocarbon systems, Calculate according to formula (8): (8); In the formula, ,base — Components i With components j The baseline value of the binary interaction parameters between them; — Components i With components j Lower bound of the binary interaction parameters between them; — Components i With components j Upper limit of binary interaction parameters between them; k 0 — The calibration coefficient (intercept term) of the logarithmic linear correlation that determines the baseline value of the binary interaction parameter kij. k1—Determine the calibration coefficients (slope terms) of the logarithmic linear correlation of the baseline value of the binary interaction parameter kij; n c,j —Components j The number of carbon atoms in hydrocarbons; T —Reservoir temperature; T ref —The reference temperature for the weak temperature-dependent correction term (usually the baseline temperature used in the calibration experiment). —Components i With components j Binary interaction parameters between them; T corr —Binary interaction parameters k ij Temperature correction factor, used to characterize k ij With temperature T Deviation from reference temperature T ref Minor corrections at that time.

[0038] Using hydrocarbon carbon number Determine a baseline value for the logarithmic linear relationship of the independent variable and multiply it by the bivariate interaction parameter. k ij Temperature correction factor T corr .

[0039] For hydrocarbon-hydrocarbon systems (i.e., components) i , j All are hydrocarbon pseudo-components, excluding CO2), binary interaction parameter baseline value. k ij,base This can be simplified to a linear function of the carbon number difference between components: ; In the formula, , —An empirical constant derived from regression calibration of experimental phase equilibrium data (its physical meaning is the same as that in formula (8) for CO2-hydrocarbon systems) k 0、 k 1. Similar, only the applicable objects are different); , —Components i , j The carbon number of representative hydrocarbons.

[0040] The linear function reflects the rule that the greater the difference in carbon number between hydrocarbon components, the greater the binary interaction parameter, which is consistent with the rule described in equation (8). Only the regression coefficient needs to be calibrated separately for the hydrocarbon-hydrocarbon system.

[0041] Then, based on the phase thermodynamic parameters of each component in the calculated PR equation of state, the multi-component phase equilibrium is solved iteratively.

[0042] Phase equilibrium constant K Value defined as component i The ratio of the mole fraction in the gas phase to the mole fraction in the liquid phase: (9); in, K i — Components i The phase equilibrium constant; — Components i Fugacity coefficient in the liquid phase; — Components i Fugacity coefficient in the gas phase.

[0043] The fugacity coefficient is calculated using formula (10): (10); in, (11); In the formula, — Components i The fugacity coefficient; b i —Components i Co-volume parameters; b mix —Co-volume parameters of the mixture; Z— Compressibility factor; obtained by solving the cubic form of the PR equation of state, with the smallest positive real root for the liquid phase and the largest positive real root for the gas phase; B —by the co-volume parameter of the mixture b mix The dimensionless parameters obtained by conversion; A —From the gravitational parameters of the mixture a mix The dimensionless parameters obtained by conversion; x j —Components jmole fraction in the mixture; a i ( T ) — Components i Temperature-corrected gravitational parameters; a j ( T ) — Components j Temperature-corrected gravitational parameters; —Components i With components j Binary interaction parameters between them; a mix —The gravitational parameters of the mixture; P— System pressure; R —Universal gas constant, R =8.314 J / (mol·K); T —Reservoir temperature.

[0044] K The initial value for iteration is obtained by formula (12): (12); K i (0) —Components i Phase equilibrium constant K Initial value for iteration; —Components i Critical pressure under confined conditions, i.e., nano-corrected critical pressure; P —System pressure; —Eccentricity factor; —Components i The critical temperature under confined conditions, i.e., the nano-corrected critical temperature; T —Reservoir temperature.

[0045] The K value is updated iteratively using a successive replacement method, specifically including: (1) With the current K Value (taken in the first iteration) K i (0) Substituting into formula (13) of the Rachford-Rice equation, we can solve for the gas phase fraction. V ; (2) V and current K Substituting the values ​​into formula (14), the liquid phase composition is obtained. x i Composition of gas phase y i ; (3) x i , y i Substitute into formulas (10) and (11) respectively to calculate the liquid phase fugacity coefficient. and gas phase fugacity coefficient ; (4) Update using the latest calculated values K value; (5) Determine convergence: Calculate the convergence in two steps. K The sum of squares of the logarithmic differences of the values, if the result is less than a preset tolerance (e.g., 1×10⁻⁶), -10 If the convergence occurs, output the current value. K value, x i , y i Otherwise, return to step (1) and iterate again with the updated K value until convergence.

[0046] In each iteration, the Rachford-Rice flash equation is called to solve for the gas phase fraction that satisfies formula (13). V : (13); In the formula, V —Gas phase fraction, which is the proportion of gas phase to the total molar amount of feed; z i —Components i Total mole fraction of feed; K i —Components i The phase equilibrium constant of the current iteration step; f ( V —Rachford-Rice objective function; The Newton-Raphson iterative method is used to solve for V.

[0047] Obtain gas phase fraction V Afterwards, the liquid and gas phase compositions are as shown in formula (14): (14); —Components iMole fraction in the liquid phase; z i —Components i Total mole fraction of feed; V —Gas phase fraction; K i —Components i The phase equilibrium constant of the current iteration step; —Components i Mole fraction in the gas phase.

[0048] The rule for judging the flash evaporation result is: if This is the total liquid phase ( V =0); if Then it is a complete gas phase ( V =1); otherwise, it indicates the coexistence of gas and liquid phases (0 < 1). V <1).

[0049] Meanwhile, the minimum miscibility pressure of each component is estimated using the power-law correlation formula (15): (15); In the formula, MMP i —Components i The corresponding minimum miscibility pressure (the unit in this calculation formula is MPa); —Components i The representative carbon number; f ( T —Temperature correction function, used to correct the minimum miscibility pressure estimated by the carbon number power-law correlation for temperature.

[0050] f ( T The result is obtained by formula (16): (16); In the formula, —Fahrenheit temperature.

[0051] Determined by experiment For key components (such as C7-C) 12 The MMP estimates are anchored and corrected to ensure that the MMP calculation results are consistent with the experimental calibration values.

[0052] Key components determined experimentally (such as C7-C) 12 )of MMP Experimental values MMP exp, compared with the estimated value of formula (15) MMP calc Perform a comparison and calculate the correction coefficient λ= MMP exp / MMP calc Apply this correction factor to other components (such as C2-C6, C4, C6) within the same crude oil system. 13+ )of MMP The estimated value is proportionally corrected, and after correction MMP i,corrected =λ×MMPi,calc, thereby ensuring that each component MMP The calculation results are consistent with the experimental calibration values.

[0053] Finally, calculate the injection pressure. The ratio of the minimum miscibility pressure of the component Then, the miscibility factor of each component is obtained by mapping using the piecewise function (17). ; (17).

[0054] (3) Establish a time-varying linkage mechanism between reservoir temperature and fluid physical parameters: This step is a common submodule that runs through multiple steps. Whenever the temperature T changes, the phase equilibrium calculation in step (2) and the LBM parameter calibration and mass transfer calculation in steps (4) / (5) need to be retried.

[0055] Using the reservoir temperature T determined in step (1) as the driving variable, time-varying linkage calculations are performed on the following fluid physical property parameters to ensure that the LBM simulation parameters always maintain physical consistency under different operating conditions: ① Time-varying correction of CO2 density: The density of supercritical CO2 decreases with increasing temperature; the current temperature... T CO2 injection density r gc (T) was obtained directly from the experiment, and the density values ​​at different temperatures were measured and used respectively.

[0056] ② Time-varying correction for CO2 viscosity: CO2 viscosity decreases as temperature increases. The CO2 lattice kinematic viscosity at the current temperature is also obtained directly from the experiment, and the collision relaxation time of the CO2 phase is updated accordingly to ensure that the flow parameters of the CO2 phase are adjusted synchronously with the measured temperature data.

[0057] ③Time-varying correction of Henry's solubility constant: The solubility of CO2 in the liquid phase decreases with increasing temperature. A time-varying calculation was performed using a solution enthalpy correction model. (Henry) T =Henrybase ×exp(ΔH diss / R ×(1 / T-1 / T base )); In the formula, Henry T —Henry solubility constant of the component at temperature T; Henry base —The Henry solubility constant, experimentally calibrated at the reference temperature; ΔH diss — Enthalpy of dissolution; approximately -16500 J / mol for light hydrocarbons and approximately -20000 J / mol for heavy hydrocarbons. A negative value indicates that dissolution is an exothermic process. R —gas constant; T —Current reservoir temperature; T base —Choose a specific temperature as a reference temperature in order to establish a correlation regarding temperature changes.

[0058] The time-varying correction is applied to the Henry constants specific to each component and takes effect in real time in the calculation of the mass transfer driving force and mass transfer of the mass transfer mechanism M1 (dissolution and desorption mass transfer of CO2 in the liquid phase) in step (5).

[0059] ④ Mass transfer rate k mt Arrhenius time-varying correction: The mass transfer rate increases with increasing temperature, and time-varying calculations are performed using the Arrhenius equation: k mt ( T )=k mt ( T base )×exp( Ea mt / R×(1 / T base -1 / T )); Where, k mt ( T )-temperature T Mass transfer rate coefficient (basic mass transfer rate common to M1-M4 mechanisms). k mt ( T base —Mass transfer rate coefficient calibrated experimentally or by molecular simulation at a reference temperature; Ea mt —Mass transfer activation energy; typical value is 12000 J / mol; R —gas constant; T —Current reservoir temperature; T base —Reference temperature.

[0060] This correction applies simultaneously to the base rate coefficients of the mass transfer mechanisms M1, M2, M3, and M4, causing the mass transfer rate to automatically increase with increasing temperature.

[0061] ⑤ Time-varying viscosity correction for each component of crude oil: Light oil components (C2-C6, C7-C) 12 Viscosity decreases with increasing temperature, which is corrected using a power-law empirical formula: nu f ( T )=nu f ( T base )×( T base / T ) 3.5 ; In the formula, nu f ( T )-temperature T Kinematic viscosity of light oil components (C2-C6, C7-C12); nu f ( T base — The kinematic viscosity of light oil components calibrated experimentally at the reference temperature; T —Current reservoir temperature; T base —Reference temperature.

[0062] Heavy oil components (C 13+ Viscosity is more sensitive to temperature, so we use: nu d ( T )=nu d ( T base )×( T base / T ) 5 ; In the formula, nu d ( T )-temperature T Kinematic viscosity of the heavy oil component (C13+); nu d ( T base— The kinematic viscosity of heavy oil components calibrated experimentally at the reference temperature; T —Current reservoir temperature; T base —Reference temperature.

[0063] The diffusion coefficient is corrected using the Stokes-Einstein method. D 0( T )= D 0( T base )×(T / T base ) / nu f ( T ); In the formula, D 0( T )-temperature T The diffusion coefficient of CO2 in crude oil; D 0( T base — The diffusion coefficient, experimentally calibrated at the reference temperature; T —Current reservoir temperature; T base —Reference temperature; nu f ( T )-temperature T Kinematic viscosity of light oil components (C2-C6, C7-C12).

[0064] The above viscosity time-varying correction results are further updated to update the relaxation time of each component, ensuring that the viscosity ratio (and thus the flowability difference) is automatically adjusted with temperature.

[0065] ⑥ Binary interaction parameters k ij Temperature time-varying correction: k ij The number decreases slightly with increasing temperature. When the temperature changes, the baseline obtained by calculating the carbon number correlation formula using formula (8) should be used. k ij The correction is made so that the binary interaction parameters automatically adjust with temperature, thereby affecting the mixture attraction parameters in the PR-EOS phase equilibrium calculation in step (2). a mix and K Value calculation result.

[0066] Items ① to ⑥ above are established based on experimental calibration values ​​at a reference temperature, and the corresponding physical property parameters are updated in real time as the temperature changes, enabling the simulation to reflect the CO2 density, viscosity, Henry's constant, mass transfer rate coefficient, diffusion coefficient, and other properties at different temperatures. k ij The real changes.

[0067] (4) Systematically calibrate LBM mesoscopic simulation parameters based on phase equilibrium calculation results: This step is the core calibration process that connects thermodynamic calculations and LBM mesoscopic simulations. It establishes three independent and traceable parameter transfer paths, mapping the thermodynamic information output from PR-EOS phase equilibrium calculations to the key mesoscopic parameters required for LBM simulations.

[0068] Calibration Path 1: Using comparative temperature T r,σ Calibration of pseudopotential characteristic density ; Calculate according to formula (18): (18); In the formula, —Pseudopotential feature density; —The baseline density can be taken as 1.0; The 1.5 in the min function is the upper limit for truncation.

[0069] Comparison temperature T r,σ Nanoscale correction of critical temperature The calculations characterize the degree to which a component deviates from its critical point at the operating temperature.

[0070] T r,σ Components with a concentration ≥1 (supercritical or near-critical state; such as CO2) T r,σ Approximately 1.32, C1 T r,σ Approximately 2.12) has a large r 0,σ Its pseudopotential function increases slowly with density, reflecting the high compressibility and low cohesion of the gas phase.

[0071] T r,σ Components with a concentration less than 1 (liquid, such as C7-C) 12 and C 13+ ) has a smaller Its pseudopotential function increases rapidly with density, reflecting the low compressibility and strong cohesion of the liquid phase, and forming a continuous transition with the gas phase components.

[0072] The pseudopotential function (effective density function) adopts an exponential form of the Shan-Chen type: (19); In the formula, —The pseudopotential function (effective density function) of component σ; r 0,σ —Pseudopotential feature density; p— Local lattice density of the component.

[0073] Furthermore, considering the additional suppression of spurious potential characteristic density by nanoconfinement in the near-wall region, a wall correction term is introduced: (20); In the formula, β σ —Dimensionless confinement strength factor, the larger the value, the more significant the wall confinement effect; x σ —The confinement ratio of component σ; s σ —The molecular dynamic diameter of component σ; d eff,σ —Effective pore size of component σ; —The pseudopotential characteristic density corrected for the near-wall region.

[0074] Calibration Path Two: Interaction strength parameters were calibrated using the misc misconception factor. : Component σ and The strength parameter of the Shan-Chen interaction between them is determined by formula (21): (twenty one); In the formula, —Components and The baseline value of the interaction strength parameter between them corresponds to the completely miscible phase. G Lower limit; —The adjustable range of the interaction strength parameter corresponds to the immiscible phase. G The upper limit of the increment above the baseline value; —Components The miscibility factor.

[0075] The components to be paired G After parameter differentiation calibration, it is necessary to satisfy CO2-heavy hydrocarbon C 13+of G Values ​​related to CO2 and light hydrocarbons (C2-C6) G A constraint that the ratio of values ​​is not less than 5 is imposed to ensure thermodynamic consistency.

[0076] If the calibrated ratio does not meet the constraint of ≥5, then the CO2-heavy hydrocarbon C ratio should be increased proportionally. 13+ of Or adjust the CO2-light hydrocarbon C2-C6 ratio. Substitute the result into formula (21) and recalculate, and repeat the verification until the constraint is satisfied.

[0077] If the requirements are still not met, further check whether there are any anomalies in the calculation of the miscibility factor of each component (especially the MMP estimate).

[0078] When injection pressure When the MMP of a certain component is much higher (misc approximately 1.0), (Minimum value) The Shan-Chen pseudopotential interaction is extremely weak, the interface between the two phases tends to disappear, and CO2 and this component exhibit miscibility.

[0079] When injection pressure When it is much lower than the MMP of a certain component (misc is approximately equal to 0). (Larger value) Strong interaction forces maintain a clear gas-liquid interface and exhibit immiscibility.

[0080] The value of determines the intensity of interfacial tension under immiscible conditions, for heavy components (such as C). 13+ Take the larger one For light components, take a smaller value .

[0081] Calibration Path 3: viscosity Calibration relaxation time First, the viscosity of crude oil was determined experimentally. Given the densities of each component, calculate the lattice kinematic viscosity of each component. lattice kinematic viscosity Calculated from physical viscosity using a dimensionless reference viscosity: ; In the formula, — The kinematic viscosity of component σ (dimensionless); —The physical kinematic viscosity of the components, calculated from the experimentally determined crude oil viscosity and the density of each component (unit: m). 2 / s); —The reference viscosity, composed of the lattice space step size Δx and the time step size Δt used in LBM simulation, is used to convert the physical viscosity into lattice units without dimension.

[0082] Then, it is mapped to the relaxation time of the MRT collision operator using formula (22). : (twenty two).

[0083] Different components The differences reflect the differences in their liquidity; supercritical CO2 and C1 have shorter relaxation times and are therefore highly liquid; heavy CO2... 13+ It has a long relaxation time and low liquidity.

[0084] For the near-wall region, a nano-wall viscosity enhancement effect is introduced, based on the confinement ratio. Regulation: (twenty three); In the formula, — The corrected collision relaxation time of component σ at the fluid node near the wall; —The wall viscosity enhancement coefficient of component σ; — The confinement ratio of component σ.

[0085] Relative body relaxation time The enhancement factor depends on the characteristic aperture. The pre-calibrated gradation increases with decreasing pore size and increasing component molecular size, reflecting the additional inhibitory effect of nanoconfinement on fluid flow near the wall.

[0086] Through the above three calibration paths, a system was established for the calibration of microfluidic experimental data ( T , P Carbon number distribution , , , Starting from ), via nano-correction of critical parameters ( , ) and PR-EOS phase equilibrium calculation ( T r, ,a(T), , , V, , ), which is ultimately mapped to LBM simulation parameters ( , , This provides a complete physical calibration link. Each LBM parameter in this link can be traced back to an upstream experimentally measurable and thermodynamically calculated intermediate quantity through an explicit formula, without the need for manual experience-based adjustments.

[0087] (5) Construct the MRT-MCMP lattice Boltzmann multiphase multicomponent flow simulation framework and embed the differential mass transfer model: S1. Based on the LBM mesoscopic simulation parameters calibrated in step (4), construct the multi-relaxation time multi-component Shan-Chen pseudopotential (MRT-MCMP) lattice Boltzmann simulation framework.

[0088] The D2Q9 discrete velocity model is adopted, which includes 9 discrete velocity directions. ( =0, 1, ..., 8). The collision operator employs a multiple relaxation time (MRT) scheme to handle systems where the viscosity differences between components can reach tens of times.

[0089] Each component The MRT-LBM evolution equation in velocity space is in the form of: (twenty four); In the formula, —migration trailing edge α Density distribution function in the direction; —Migration Frontier α Density distribution function in the direction; M —A 9×9 constant transformation matrix, wherein each row vector is pairwise orthogonal, satisfying the following condition: , D It is a diagonal matrix. M 1 It can be pre-calculated and stored once before the simulation begins, avoiding the extra computational overhead caused by repeated inversion at each time step; Λ σ —Diagonal relaxation matrix; —Density distribution function along the β direction of the migration front; —Equilibrium distribution function along the β direction; Δt —Time step; —Includes corrections for contributions from external forces.

[0090] The Λ σ The diagonal elements are: (25); The equilibrium distribution function is given by the second-order Taylor expansion of the Maxwell-Boltzmann distribution: (26); In the formula, — α The weighting coefficients for the direction; where, w 0 = 4 / 9 w 1-4 =1 / 9, w 5-8 =1 / 36; (0-8 correspond to the origin, east, north, west, south, northeast, northwest, southwest, and southeast directions). r σ —Components The macroscopic density is obtained by summing the distribution functions of the components in each direction: ; e α —Discrete velocity vector in direction α; —Components The equilibrium velocity; —The square of the lattice sound speed, with a value of 1 / 3; —The fourth power of the speed of sound in a lattice is 1 / 9.

[0091] Phase separation or mixing between components is achieved through Shan-Chen pseudopotential interactions. Components Subject to the components The pseudo-potential interaction force is calculated using the following formula: (27); In the formula, —The interaction strength parameters calibrated in step (4); —Components The pseudopotential value at position x; ; in, r 0, σ —The pseudopotential characteristic density of component σ; r σ ( x — The local macroscopic lattice density of component σ at position x.

[0092] Components It is also subject to fluid-solid wall forces, which are used to control wettability: (28); In the formula, —Interaction force between component σ and the wall surface; —Directional weighting coefficient; —Solid indicator function; where solid nodes are 1 and fluid nodes are 0; e α —Discrete velocity vector in direction α; —Fluid-solid interaction strength parameters; —Components The pseudopotential value at position x.

[0093] In addition, the fluid is also subjected to bulk forces. (For example, the displacement pressure gradient). The total external force on each component is: (29); In the formula, —The total external force acting on component σ at position x; —The pseudo-potential interaction force between components is given by formula (27); —The fluid-solid wall interaction force given by formula (28); — Volume driving force (such as constant volume force generated by displacement pressure gradient).

[0094] Components The macroscopic velocity is determined by the conservation of momentum and the contribution of external forces: (30); In the formula, —The macroscopic velocity of component σ itself; —The distribution function of component σ along direction α; —The discrete velocity vector in direction α; —The macroscopic density of component σ (given by formula 26); —Time step; —The total external force on component σ (given by formula 29); —The mixing equilibrium rate shared by all components, based on the density of each component. The weighted average is obtained.

[0095] The solid boundary conditions employ a half-step bounce scheme. The inlet boundary conditions use a non-equilibrium extrapolation scheme, decomposing the distribution function of the boundary nodes into equilibrium and non-equilibrium parts. The equilibrium part is constructed from the specified density and velocity at the inlet; the non-equilibrium part is approximated by the non-equilibrium parts of adjacent fluid nodes. The outlet boundary conditions employ fully developed convection outlet conditions.

[0096] S2. In each computational step of the LBM main loop, the following differentiated multi-mechanism mass transfer model is embedded. Unlike existing technologies that apply the same mass transfer rules to all components, this invention sets differentiated mass transfer parameters for different components that follow a unified gradient law, covering four types of mass transfer mechanisms.

[0097] Mass transfer mechanism M1: CO2 dissolution and desorption mass transfer in the liquid phase, driven by Henry's Law: For each fluid grid point, the effective Henry's constant is calculated weighted according to the current composition of the local liquid phase. (31); in, —Effective Henry constant; —Components The proportion in the current liquid phase; —Components A specific Henry's law of solubility; r σ —The local macroscopic density of component σ at this grid point; —Crude oil density at this grid point.

[0098] The setting follows the physical law that the heavier the component, the weaker its solubility; methane C1... H C1 Minimum, approximately 0.020, because methane is almost insoluble in CO2 under supercritical conditions; C2-C6 H C2-C6 Approximately 0.18; C7-C 12 of H C7-C12 Approximately 0.25 (baseline value); C 13+ of H C13+ Approximately 0.12 (heavy macromolecules have small gaps between them, making it difficult for CO2 molecules to be inserted).

[0099] The equilibrium concentration is: (32); in, —Local CO2 partial pressure calculated from local CO2 gas phase density; —CO2 gas phase density; —The square of the lattice sound speed, with a value of 1 / 3; —Pressure conversion factor, used to convert the gas phase density of lattice units into physical pressure units; —Local equilibrium dissolution concentration, which is the maximum equilibrium concentration of CO2 that can dissolve in the liquid phase under the current local pressure and temperature conditions, according to Henry's Law; —Effective Henry constant; C max This is the preset maximum solubility concentration limit.

[0100] The mass transfer driving force is: (33); in, —The driving force of mass transfer, the sign of which determines the direction of mass transfer (ΔC>0 indicates that the local dissolved CO2 concentration is lower than the equilibrium concentration, driving dissolution; ΔC<0 indicates that it is supersaturated, driving desorption). —Local equilibrium dissolution concentration; —Current dissolved CO2 concentration.

[0101] If the absolute value of the driving force exceeds the preset lower limit (e.g.) If the mass transfer is not greater than a certain proportion of the total amount of CO2 gas, then the mass transfer is calculated and limited. The mass transfer due to dissolution shall not exceed a certain proportion of the total amount of CO2 gas before dissolution, and the mass transfer due to desorption shall not exceed a certain proportion of the total amount of dissolved CO2.

[0102] Mass transfer is achieved by transferring mass from the CO2 gas phase distribution function to the dissolved CO2 auxiliary distribution function (dissolution) or vice versa (desorption).

[0103] The density field of the dissolved CO2 phase is solved independently using a set of D2Q9 distribution functions (isomorphic to the other components and involved in the calculation of inter-component interaction forces) according to the general pseudopotential framework described in formulas (19)-(24) and (27), denoted as Its macroscopic density Unlike the "repulsion" between gaseous CO2 and oil phase components (to support the gas-liquid interface), the interaction strength parameter between dissolved CO2 and each pseudo-component in the oil phase results in a "weak attraction," indicating that it has participated in volume expansion and viscosity changes as an internal component of the liquid phase, and no longer exhibits independent gas phase flow behavior.

[0104] If directly using this density field As the mass transfer driving force of formula (33) The concentration input is affected by non-physical density fluctuations near the gas-liquid interface due to spurious force fields, making the determination of the dissolution / desorption direction sensitive to local numerical noise. Therefore, an independent concentration field distribution function is added. The D2Q5 discrete velocity model (not involved in the force field construction of formula (27)) is used specifically to provide the clean local concentration signal required by formula (33), and its evolution equation is: ; ; In the formula, —Concentration field relaxation time; —Liquid phase transport rates obtained by density-weighting of each pseudo-component in the oil phase with dissolved CO2; —The source term is provided by the mass transfer rate. The concentration field is taken as the local Henry equilibrium concentration at the gas-liquid interface monolayer grid. Equation 31 represents the dynamic boundary condition, updated step-by-step with the gaseous CO2 density field; the concentration field at the inlet boundary (where only pure CO2 gas is injected and has not yet dissolved) is set to zero. Therefore, the mass transfer driving force and mass transfer rate are: ; ; In the formula, —Driving force for component mass transfer; —Local equilibrium dissolution concentration (given by formula (32)); —The current local dissolved CO2 concentration is taken from the above concentration field distribution function. The macroscopic value, not the dissolved CO2 density field. ; —The mass transfer rate of mechanism M1 at this grid point and at this time; --temperature T The mass transfer rate coefficient under the given temperature time-varying linkage mechanism (already given by the temperature time-varying linkage mechanism above); —The matriculation enhancement coefficient reflects the physical trend that the higher the degree of matriculation, the smaller the mass transfer resistance and the faster the dissolution / desorption process; —The overall miscibility factor of the system (given by formula (41)). This mass transfer rate is simultaneously injected into three locations equally according to the directional weighting coefficient: gas phase CO2 distribution function (Loss of mass), dissolved CO2 density field (Gain mass and update its contribution to the interaction forces), concentration field (Gain quality, update the criterion signal for the next time step). All three share the same transition amount and the same time step, and... Using the same lattice mass density dimension as the dissolved CO2 density field, the source term can adapt to both types of fields simultaneously without unit conversion.

[0105] Mass transfer mechanism M2: CO2 fractional extraction mass transfer of each displaced component: targeting C1, C2-C6, C7-C 12 Each extraction is performed independently.

[0106] Extraction is triggered when two conditions are met simultaneously: the local CO2 gas phase density is greater than the CO2 concentration threshold corresponding to that component. Furthermore, the local density of this component is greater than its corresponding oil phase concentration threshold. .

[0107] and They all follow the gradient principle that the lighter the component, the lower the threshold, for example... Approximately 0.06 (the lowest threshold, as C1 supercritical fluids are almost impossible to extract). Approximately 0.10; Approximately 0.15 (baseline).

[0108] The extraction rate formula is: (34); In the formula, m extract,σ —The extraction mass transfer rate of component σ (the mass extracted per unit time). t —Current simulation time; —The extraction rate coefficient of this component; following a gradient law where the smaller the carbon number, the higher the rate; among which, Approximately 0.008-0.010, Approximately 0.006-0.008, Approximately 0.006); —This component is in local phase equilibrium in the partition to which the current grid point belongs. K The value is assigned by the local K-value field provided by the partitioned local phase equilibrium dynamic feedback coupling mechanism (Equation 39–41), reflecting the tendency of the component to be phase equilibriumly distributed to the gas phase under the current local thermodynamic conditions. The larger the value, the easier it is to be extracted into the gas phase. —The volatility factor of the component; used to reflect the amplifying effect of the difference in volatility of the component on the extraction rate. The larger the value, the stronger the volatility and the faster the extraction. See the following examples for specific values. —The miscibility adjustment factor of this component is used to adjust the extraction rate according to the current local miscibility of this component. It decreases as the miscibility increases, reflecting the inverse relationship between the immiscible interface mechanism of extraction and the miscibility mixing mechanism. r σ —The local macroscopic density of component σ; —CO2 gas phase density; —The square of the speed of sound in the lattice is 1 / 3.

[0109] Mass transfer mechanism M3: CO2 affects heavy hydrocarbons C 13+ Limited extraction mass transfer: isomorphic to M2, but with the volatility factor term removed: M3 is formally similar to M2, but has essential differences in parameter settings.

[0110] (35); In the formula, —C 13+ Extraction rate coefficient of the component; —C 13+ The local phase equilibrium K value of the component (approximately 0.0005-0.001, fully reflecting the physical nature of heavy components being difficult, slow, and requiring little extraction). —C 13+ The miscibility factor of the components; —C 13+ Local density of components; —CO2 gas phase density; —The square of the speed of sound in a lattice.

[0111] Mass transfer mechanism M4: CO2 solvent effect-enhanced joint desorption mass transfer: When the proportion of local total CO2 concentration (gas phase + dissolved state) exceeds the preset trigger threshold and the liquid phase still exists, CO2 acts like an organic solvent, weakening the cohesive force between hydrocarbon molecules in the liquid phase and enhancing the desorption of each component from the liquid phase.

[0112] The local CO2 percentage (enhancing factor) is: (36); In the formula, —Local CO2 percentage (enhancing factor); —Local density of CO2 in the gas phase; —Local density of dissolved CO2; —The sum of the densities of all oil phase pseudo-components at this grid point.

[0113] For C1, C2-C6, C7-C12 For light and medium components, the separation rate is: (37); In the formula, —Component σ M 4. Separation rate; —The fundamental separation coefficient due to solvent effect; —The local CO2 percentage is given by formula (36); —Selectivity factor of component σ; —Local density of component σ; — The miscibility adjustment factor of component σ.

[0114] For heavy component C 13+ The separation rate is: (38); In the formula, ——Component C 13+ of M 4. Separation rate; —The fundamental separation coefficient due to solvent effect; —The local CO2 percentage is given by formula (36); ——Component C 13+ The selectivity factor; —Local density of component σ; ——Component C 13+ The miscibility adjustment factor.

[0115] It follows the gradient principle that the lighter the component, the easier it is to be carried out by the CO2 solvent.

[0116] The difference between equations (37) and (38) reflects the differentiated design of the solvent effect detachment mechanism of the heavy component in this invention: even under immiscible conditions, the light and medium components ( →0) still retains basic separation ability (lower limit of enhancement coefficient is 1), while the heavy component C 13+ The solvent effect of CO2 increases linearly from zero as the degree of miscibility increases—that is, only after the system gradually becomes more miscible can the solvent effect of CO2 be sufficient to overcome the strong cohesive forces between heavy macromolecules and separate them from the liquid phase. This setting is consistent with the thermodynamic nature of heavy components having high MMP and being difficult to miscible.

[0117] All differential parameters among the four mass transfer mechanisms mentioned above Following a unified gradient design principle—from C1 to C 13+ The mass transfer threshold monotonically increases, the mass transfer rate coefficient monotonically decreases, and the volatility / selectivity factor monotonically decreases.

[0118] Simultaneously, a dynamic feedback coupling mechanism based on partitioned local phase balance is introduced: every preset number of simulation steps... (e.g., every 200 steps), divide the entire computational domain into N parts along the x and y directions respectively. x and N y The system consists of 15 rectangular partitions (e.g., 5 x 3). To avoid instability in flash evaporation calculations due to excessively low local mass, a minimum partition mass threshold (min) can be set. mass_per_zone If the total density of non-solid mesh points within a partition is lower than this threshold, the local flash calculation for that partition is skipped in this round. This mechanism corresponds to step (2) described above.

[0119] For each partition, calculate the current lattice density of each component in all non-solid grid points within that partition, and convert it into the local molar composition of the partition: ; In the formula, —The local molar composition of component σ in this partition; , —Components , Local density at grid point x; , —Components , The molar mass; σ′\sigma' σ′——A dummy variable summing all components (including σ itself) within the partition.

[0120] The specific calculation method is as follows: ; In the formula, —Local estimated pressure in this partition; —The average density of gaseous CO2 at all non-solid grid points within this partition; — The square of the speed of sound in a lattice (=1 / 3); —Pressure conversion factor (consistent with the definition in Formula 32).

[0121] Composed of local moles in partitions The local pressure of the zone estimated by the above formula For input, repeat step (2) completely. Substitute into the PR state equation formulas (4), (5), and (11) to calculate. , and , Calculate according to formula (12) K The initial value is iterated, and the solution is obtained iteratively using the successive substitution method. K The values ​​and fugacity coefficients are substituted into the Rachford-Rice equation (13) to solve for the gas phase fraction of the partition. To obtain the local partition K Based on the field and local gas-liquid composition, the local enrichment coefficient is finally calculated: (39); In the formula, —Local enrichment coefficient of component σ; —The mole fraction of the vapor phase obtained from the local flash evaporation calculation of this zone; —The liquid phase mole fraction calculated from the local flash evaporation of this zone. The larger the value, the more likely the component is to be enriched in the gas phase under the current local conditions.

[0122] This reflects the extent to which the component is extracted into the gas phase under the current local composition conditions.

[0123] Based on the local enrichment coefficient, the global phase equilibrium K-value field and the interphase interaction intensity parameters are analyzed. G Perform a smooth update: (40); In the formula, —— K Value smoothing update coefficient (typical value 0.15-0.30); —— G Parameter smoothing coefficient (typical value 0.15); —The local K value obtained from this partitioning calculation; —Determined by local enrichment coefficients G Target value for parameters.

[0124] If the enrichment coefficient of the middle-grade component is high, it should be appropriately reduced. G Parameters are used to promote miscibility, and a higher enrichment factor is maintained if the enrichment factor of the heavy component is low. G Parameters to maintain a clear interface.

[0125] The smoothing process avoids abrupt parameter changes and numerical oscillations caused by compositional differences between adjacent partitions. This dynamic feedback mechanism creates a complete two-way coupled closed loop between thermodynamic calculations and LBM fluid dynamics simulations.

[0126] Based on the aforementioned dynamic feedback mechanism for phase equilibrium, the following time-varying parameter update mechanism within the simulation time step is further established to accurately capture the dynamic evolution characteristics of each physical quantity over time during CO2 displacement: ① Miscibility factor The time-varying update of this factor consists of a weighted fusion of two parts: the first part is a thermodynamic pressure term based on the ratio of global injection pressure to MMP, reflecting the static miscibility potential under fixed injection conditions; the second part is based on the interfacial tension at the current moment. IFT With the initial IFT The ratio is calculated using a preset mapping relationship for the interface status items. The latter gradually increases as the CO2 displacement process progresses and the interfacial tension decreases, causing the miscibility factor to evolve dynamically over time rather than remain at a fixed value.

[0127] Furthermore, a miscibility factor is constructed. , It is obtained by weighted average of two parts: (41); In the formula, —Overall miscibility factor of the system (time variable); —Static thermodynamic potential term based on the ratio of global injection pressure P to minimum miscibility pressure MMP, P / MMP; —A dynamic interface state item based on the ratio of the current interface tension to the initial interface tension; , —The weighting coefficients of the two items satisfy + =1.

[0128] These two parts respectively reflect the thermodynamic miscibility under global pressure conditions and the interfacial miscibility state reflected in the actual evolution of the flow field; the secondary clamping avoids overestimating the miscibility in the early stage when CO2 diffusion is insufficient.

[0129] Part One Based on injection pressure Ratio to minimum miscibility pressure According to the miscibility factor of a single component The same piecewise function is used to calculate and reflects the degree of thermodynamic miscibility under global pressure conditions: (42); In the formula, —The ratio of injection pressure to the minimum miscibility pressure of the component; --based on The static thermodynamic potential term.

[0130] Part Two Based on the ratio of the current interfacial tension to the initial interfacial tension during the simulation process : (43); In the formula, —The ratio of current interface tension to initial interface tension; --based on Dynamic interface status items.

[0131] The results are obtained through calculations based on a pre-defined mapping relationship, reflecting the interfacial miscibility state embodied in the actual evolution of the flow field.

[0132] ②Shan-Chen interaction strength parameters G Time-varying adaptive: During each partition balancing update, the enrichment coefficient of each partition is used as a basis. E σ Adaptive adjustment of CO2-middle-component interaction strength G fg Interaction strength with CO2-heavy components G dg The target value is determined by the smoothing coefficient α. G Asymptotic updates toward the target value, making G The parameters evolve dynamically with the displacement process.

[0133] Calibration path two (formula 21) gives the static calibration value (baseline value) of the G parameter, and this step is a dynamic adaptive update superimposed on this.

[0134] The higher the degree of extraction and enrichment, the better. G fg The smaller the value, the more miscible the phase; the lower the enrichment coefficient of the heavy component, the more likely it is to be miscible. G dg Maintain a large value to preserve the interface.

[0135] ③ Time-varying decomposition of the relative contributions of miscible extraction and immiscible extraction: In the four-mechanism decomposition in step (5), the contribution ratio of each mechanism varies with... It is continuously updated over time (see step (6) for the specific decomposition algorithm), and presents a quantifiable time-varying decomposition curve from the early stage of displacement (dominated by non-miscible extraction) to the middle and late stages of displacement (increased proportion of miscible extraction).

[0136] The code will automatically output data; simply record the data.

[0137] The above-mentioned time-varying mechanism and the temperature time-varying linkage correction in step (3) together constitute the complete parameter time-varying system of the present invention, ensuring the physical accuracy of the simulation results throughout the entire displacement time history.

[0138] (6) Quantitative decomposition of mass transfer mechanism based on multiple independent physical sources and self-verification of mass conservation closure: During the simulation, the density field of each component, the cumulative mass transfer of each mass transfer mechanism, the expansion factor SF, and the sweep efficiency are recorded at preset output intervals (e.g., every 100 steps). For each output time step, execute the following four-mechanism quantitative decomposition algorithm: ① Calculate the contribution of the expansion mechanism The value is determined by the larger of two independent physical sources, with a monotonically non-decreasing constraint applied. Source A (Inflation Factor Method): (44); In the formula, SF— Inflation factor; E v — Sweep efficiency; α sw —Expansion efficiency coefficient, typically around 0.55.

[0139] Source B (Net Solubility Method): (45); In the formula, —M 1. Net cumulative amount of dissolution (dissolution is positive, desorption is negative, and the amounts are accumulated with signs).

[0140] (46); In the formula, EOR —The current contribution of the expansion mechanism to this step; EOR —The expansion mechanism contribution from the previous step.

[0141] ② Calculate the extraction contribution of immiscible interfaces and the contribution of miscible extraction .

[0142] First, the cumulative mass transfer of each mass transfer mechanism is classified according to the light and heavy components of the displaced fluid: (47); In the formula, m light,MT — Total mass transfer of light components; M 2C1 — C1 extraction amount; M 2C2C6 — Extraction amount of C2-C6; M 2 b C7C12 — C7-C 12 Extraction volume; M 4cumul · f light,init—M 4. Lightweight components; (48); In the formula, m heavy,MT — Total mass transfer of heavy components; M 3 C13+ — C 13+ Extraction volume; M 4 cumul · f heavy,init —M 4. Heavy components.

[0143] Total mass transfer of light components For C1 / C2-C6 / C7-C 12 Extraction amount added M 4. Light component; Total mass transfer of heavy component C 13+ Extraction amount added M 4. Heavy component. Mass transfer of the light component. Direct allocation. A nonlinear decay correction for miscibility of heavy components is introduced: (49).

[0144] The physical basis for this square correction is: C 13+ The minimum miscibility pressure (estimated at approximately 65 MPa by the carbon number power-law correlation) is much higher than that of C7-C. 12 At the MMP (approximately 21 MPa), the probability of heavy components actually reaching a miscible state and being efficiently extracted under the same injection pressure is much lower than that of medium and light components. Therefore, a quadratic nonlinear suppression is applied to its miscible distribution ratio to accurately reflect the physical reality that heavy components are difficult to truly miscible under conventional injection pressure.

[0145] Then, the extraction contribution of the immiscible interface is calculated according to the following formula. and miscible extraction contribution ; (50); (51); In the formula, EOR mis — Miscible extraction contributes to oil recovery; m light,MT — Total mass transfer of light components; — Mass transfer of light components; m heavy,MT — Total mass transfer of heavy components; EOR immis — Extraction from immiscible interfaces contributes to the recovery rate.

[0146] ③ Determine the contribution of the displacement and carryover mechanism: (52); In the formula, EOR disp — Displacement and carryover mechanisms contribute to oil recovery; EOR total — Total recovery rate; EOR swell — Expansion mechanism contributes to oil recovery; EOR mis — Miscible extraction contributes to oil recovery; EOR immis — Extraction from immiscible interfaces contributes to the recovery rate.

[0147] Taking non-negative values, the physical logic for determining the displacement contribution using the residual method is as follows: expansion, immiscible extraction, and miscibility each have their own independent and non-overlapping physical calculation sources, while displacement (directly pushing the displaced fluid out of the pores through viscous propulsion and pressure gradient) is physically the remaining part after subtracting the mass transfer and expansion-related contributions from the total amount.

[0148] ④ Self-verification of mass conservation closure error: (53); If | | Less than the preset tolerance (e.g.) If the result satisfies the law of conservation of matter, the calculation process of each sub-step is correct, and the physical rationality of the decomposition result is verified.

[0149] (7) Based on the complete model constructed in steps (1) to (6) (covering thermodynamic phase equilibrium, LBM mesoscopic parameter calibration, temperature time-varying correction, differentiated mass transfer mechanism and partitioned feedback coupling, mechanism decomposition and conservation verification), simulation calculations are performed to obtain the spatiotemporal distribution of density fields of each component, cumulative mass transfer time series of each mass transfer mechanism, mass transfer efficiency curves of the overall and component levels, spatiotemporal evolution of saturation of each phase and gas-liquid phase composition, decomposition curves of contributions of the four mass transfer mechanisms and closure error verification results, spatiotemporal evolution of phase equilibrium K value field, and dynamic change data of key evaluation indicators such as interfacial tension, capillary number, expansion factor and sweep efficiency. The above data are then graphically presented using visualization software.

[0150] The beneficial effects of this invention are as follows: The method described in this invention establishes a phase equilibrium calculation based on the PR equation of state, starting from macroscopic physical property parameters measurable in microfluidic experiments (temperature, pressure, crude oil carbon number distribution and density, reservoir pore size, etc.), and corrected by nanoconfinement. It systematically calibrates LBM mesoscopic simulation parameters (pseudopotential characteristic density, Shan-Chen interaction strength parameters, relaxation time, etc.), establishes a time-varying linkage mechanism between reservoir temperature and multiple physical property parameters such as CO2 density, viscosity, Henry's dissolution constant, mass transfer rate, crude oil viscosity, and diffusion coefficient, and achieves phase equilibrium in the main simulation cycle. K Value field, Shan-Chen interaction strength parameters G The model incorporates a multi-component mass transfer model that covers at least four different mass transfer mechanisms, and ultimately achieves a complete simulation and analysis method for quantitative decomposition and mass conservation self-verification of the contributions of four mass transfer mechanisms: displacement and carry-over, dissolution and swelling, immiscible interface extraction, and miscible extraction.

[0151] The specific advantages are as follows: (1) This invention establishes a complete physical parameter transfer link that starts from macroscopic physical property parameters measurable in microfluidic experiments (temperature, pressure, crude oil carbon number distribution and density, reservoir pore size and minimum miscibility pressure), and performs phase equilibrium calculations based on the PR equation of state with nanoconfinement correction. This link is then systematically mapped to LBM mesoscopic simulation parameters through three independent calibration paths (calibrating pseudopotential characteristic density by comparing temperature, calibrating Shan-Chen interaction strength parameters by miscibility factor, and calibrating relaxation time by viscosity). Each LBM parameter in this link can be traced back and calculated from upstream experimental data using explicit formulas. This fundamentally differs from the existing approach of relying on trial and error to set LBM parameters based on human experience, ensuring the physical consistency of simulation parameters and their portability under different reservoir conditions.

[0152] (2) In the differentiated multi-mechanism mass transfer model of the present invention, the settings of all differentiated parameters (Henry constant, extraction threshold, rate coefficient, volatility factor, selectivity factor, etc.) follow the gradient of K values ​​of each component derived from the phase equilibrium calculation of PR-EOS (K C1 Much larger than K C2C6 Greater than K C7C12 Much larger than K C13+ The self-consistent gradient law—components with higher thermodynamic volatility are assigned lower extraction thresholds and higher extraction rates in the mass transfer model. This thermodynamically self-consistent parameter design enables the present invention to accurately output component-level differentiated recovery curves. Simultaneously, a dynamic feedback coupling mechanism based on partitioned phase equilibrium achieves a two-way coupling closed loop between thermodynamic calculations and LBM fluid dynamics simulations, overcoming the limitation that globally fixed parameters cannot reflect local component differences.

[0153] (3) The quantitative decomposition algorithm for the four-mechanism mass transfer contribution proposed in this invention calculates the contribution of each mechanism through multiple independent physical sources (expansion comes from the dual sources of expansion factor and net solubility of M1, taking the larger value and adding monotonicity constraint; immiscible / miscible comes from the mass transfer accumulation allocated according to the degree of miscibility factor; displacement comes from the residual method), and introduces the squared decay correction of the miscibility accessibility of heavy components ( This method demonstrates the nonlinear threshold effect of its high MMP (Metal-to-Main-Metal-Potential) and ultimately ensures the rationality of the material conservation in the decomposition results through self-verification of closure error. It provides engineers with a quantitative basis for determining the priority mass transfer mechanisms to be controlled under specific reservoir conditions, overcoming the limitation of existing technologies that can only output the total recovery curve but cannot distinguish the contributions of each mechanism.

[0154] (4) The entire simulation framework of this invention uses experimental data such as microfluidic experiments as input sources, and the test results are accurate, with good engineering applicability and scalability. Although the specific embodiment uses CO2 displacement of crude oil as a scenario, the method framework of this invention is applicable to any engineering scenario involving the differentiated mass transfer behavior between multiple fluid components in porous media, including carbon dioxide geological storage, groundwater pollution remediation, and micro / nanofluidic separation.

[0155] (5) This invention establishes a complete time-varying parameter system, which systematically solves the limitations of existing technologies in terms of fixed simulation parameters and inability to reflect dynamic processes from two dimensions. The first dimension is the time-varying linkage mechanism of operating temperature (step 3): with reservoir temperature T as the driving variable, the CO2 density, CO2 viscosity, Henry's dissolution constant of each component, and mass transfer rate are respectively adjusted by table lookup interpolation, van't Hoff equation and Arrhenius equation. k mt Viscosities of each component of crude oil (using power-law corrections of exponents 3.5 and 5.0 to differentiate the temperature sensitivity differences between light and heavy oils), diffusion coefficients, and binary interaction parameters.k ij To achieve fully automatic calculation with linked parameters, only a single temperature input value needs to be modified in the configuration module. T factor This allows for the automatic recalculation of all downstream physical properties and LBM mesoscopic parameters, supporting comparative analysis of multi-condition systems across a wide temperature range of 60°C to 150°C. This fundamentally eliminates the operational burden and risk of parameter inconsistency associated with manual parameter modification in multi-condition studies. The second dimension is the time-varying dynamic evolution mechanism of the displacement process (step 5): within each time step of the LBM main cycle, the phase equilibrium K-value field is dynamically updated every 200 steps through partitioned local flash evaporation calculations, reflecting the impact of local component changes on phase equilibrium; the miscibility factor... By incorporating the time-varying term of interface tension ( IFT current / IFT initial The Shan-Chen interaction strength parameter evolves continuously with the displacement process. G The system dynamically adjusts over time using an adaptive adjustment mechanism based on the enrichment coefficient Eσ. This two-dimensional time-varying system enables the present invention to output continuous time-varying parameter evolution curves and mechanistic decomposition dynamic curves throughout the entire displacement process, providing researchers with a unique means to quantitatively evaluate the evolution of the main controlling mechanisms at each stage of CO2 displacement over time.

[0156] This invention can be widely applied to engineering fields involving multiphase and multicomponent flow and differentiated mass transfer in porous media, such as oil and gas reservoir gas injection displacement simulation, carbon dioxide geological storage evaluation, groundwater pollutant migration and remediation prediction, and multicomponent separation process design in micro-nano fluid control devices. Attached Figure Description

[0157] Figure 1 This is a distribution map of crude oil components.

[0158] Figure 2 This is the viscosity-temperature curve of crude oil.

[0159] Figure 3 This is a model framework diagram.

[0160] Figure 4 The fitting curve for the experimental-simulated dynamic expansion factor.

[0161] Figure 5 This is a map showing the distribution of oil and gas at different times during the CO2 flooding process at 25 MPa.

[0162] Figure 6 The graph shows the variation of microscopic oil displacement efficiency of each component in CO2 flooding at 25 MPa over time. Detailed Implementation

[0163] The technical solution of the present invention will be described in detail below with reference to the accompanying drawings and embodiments.

[0164] 1. This embodiment takes a deep oil reservoir supercritical CO2 injection displacement crude oil system containing multiple hydrocarbon components as an example to demonstrate in detail the complete implementation process of the method of the present invention.

[0165] Obtain or generate the porous media geometry model used in the implementation: If a scanning electron microscope image of a real rock core exists, the binarization process of ImageJ software (enhancing contrast, 8-bit conversion, grayscale threshold adjustment, and noise reduction) is combined with the imread, im2bw, and imresize functions of Matlab to obtain a two-dimensional porous medium array S (0 = pore / fluid nodes, 1 = solid skeleton nodes) of lx multiplied by ly (300 multiplied by 150 in this embodiment), and the porosity is statistically analyzed.

[0166] If no real core image is available, an artificial porous medium with a specified porosity (approximately 0.38-0.42) is generated using a random circular particle stacking algorithm.

[0167] 2. The specific values ​​used in the examples are derived from measured data of high-temperature and high-pressure microfluidic experiments or from standard conversions based on experimental data.

[0168] Example 1 The simulation method for differential mass transfer of multiphase and multicomponent fluids based on LBM includes the following specific steps: (1) Constructing a multi-component thermodynamic system based on microfluidic experimental data: S1. The following reservoir physical property data were obtained through high-temperature and high-pressure microfluidic experiments: Reservoir temperature T = 388.15 K (115 °C), initial reservoir pressure P reservoir =55.07MPa, injection pressure =25MPa, characteristic pore size of reservoir rock ≈10nm (obtained from statistical analysis of scanning electron microscope images), crude oil density =843.5kg / m 3 Crude oil viscosity ≈5.45 mPa·s, minimum miscibility pressure of CO2-crude oil =23.73MPa (directly measured by microfluidic experiments).

[0169] The carbon number distribution of crude oil was obtained by analyzing the components of crude oil using gas chromatography-mass spectrometry (GC-MS). The mass fractions of each monomer hydrocarbon from C1 to C10 were obtained by GC-MS: C2=1%, C3=5%, iC4=2%, nC4=5%, iC5=3%, nC5=4%, C6=7%, C7=13%, C8=12%, C9=9%, C 10 =7% (69.5% in total, excluding C1, as C1 is a dissolved gas and is not obtained by GC-MS, see GOR conversion later). Additionally, the heavy component (C) was independently determined by SARA family composition analysis. 11 The content of (including gum, asphalt, etc.) is 34%.

[0170] Since GC-MS light-medium end analysis and SARA heavy end analysis are two independent testing methods, their data cannot be directly added to 100%. Therefore, C1 (GOR conversion fixed value) and C... 13 + (SARA fixed value) is the anchor point, for C2-C6, C7-C 12 Normalize according to the original GC-MS ratio, see S3 for details.

[0171] S2. Based on the obtained carbon number distribution data of the reservoir fluid crude oil, the crude oil is divided into different pseudo-components: methane (C1, dissolved gas); light hydrocarbons (C2-C6, represented by pentane); and medium hydrocarbons (C7-C6, represented by octane). 12 Heavy hydrocarbons, represented by hexadecane (C 13+ .

[0172] S3. First, sum the carbon number distributions above item by item according to the pseudo-component boundaries to obtain the mass fraction of each pseudo-component: The C1 content of methane, converted from a gas-oil ratio (GOR) of approximately 12, is approximately 1.5%. C1 (dissolved gas components) is not obtained by summing the carbon number distributions using GC-MS. This is because GC-MS typically analyzes only degassed crude oil samples from the ground and does not include dissolved light gases; therefore, it must be calculated separately from the field-measured gas-oil ratio (GOR). C2-C6 / C7-C12 / C13+, on the other hand, are obtained by summing the carbon number distributions using GC-MS. The specific calculation is as follows: ; In the formula, GOR is the gas-oil ratio (dissolved gas volume under standard conditions / crude oil volume at the surface, taken as 12 m³ in this example). 3 / m 3 ), —Density of methane under standard conditions —The density of degassed crude oil under standard conditions, when substituted into the formula, is approximately 1.5%.

[0173] The mass fraction of light hydrocarbons C2-C6 = 1 + 5 + 2 + 5 + 3 + 4 + 7 = 27% (accounting for 27% of the total oil); Medium hydrocarbons C7-C 12 The mass fractions were determined by GC-MS as C7, C8, C9, and C6. 10 The sum of the four terms is 13 + 12 + 9 + 7 = 41%; Heavy hydrocarbons C 13+ The mass fraction of heavy components such as asphaltene (determined independently by SARA family composition analysis) is 34%.

[0174] C1 and C 13+ The proportions are fixed, C2-C6 and C7-C 12 It needs to be normalized according to the original proportion.

[0175] Ultimately, the mass fractions of the four pseudo-components were determined as follows: C1=1.5%, C2-C6=25.61%, C7-C 12 =38.89%, C 13+ =34.00%, checksum =100%.

[0176] Then, the mole fraction of each pseudo-component is determined based on its mass fraction. Conversion: Pseudo-component mole fraction = .

[0177] Each pseudo-component i The molecular weights are as follows: methane (C1) molecular weight = 16.04 g / mol; light hydrocarbons (C2-C6, represented by pentane) molecular weight = 72.15 g / mol; and medium hydrocarbons (C7-C6) molecular weight = 72.15 g / mol. 12 (Using octane as a representative) Molecular weight = 114.23 g / mol, heavy hydrocarbon C 13+ (Taking hexadecane as an example) Molecular weight = 226.44 g / mol.

[0178] Calculate item by item: mole fraction of pseudo-component C1 =0.015 / 16.04=0.000935; mole fraction of pseudo-components C2-C6 =0.2561 / 72.15=0.003550; pseudo-component C7-C 12 mole fraction =0.3889 / 114.23=0.003405; pseudo-component C 13+ mole fraction =0.34 / 226.44=0.001502.

[0179] Total = 0.009392.

[0180] Finally, the obtained pseudo-component mole fractions Normalization yields the hydrocarbon mole fraction: The methane mole fraction C1 = 0.000935 / 0.009392 = 0.0996; The mole fraction of light hydrocarbons C2-C6 = 0.003550 / 0.009392 = 0.3780; Medium hydrocarbon mole fraction C7-C 12 =0.003405 / 0.009392=0.3625; Heavy hydrocarbon mole fraction C 13+ =0.001502 / 0.009392=0.1599.

[0181] S4. Set the feed mole fraction for each component: First, the CO2 feed mole fraction is introduced. The value is 5%; multiply the mole fractions of each hydrocarbon obtained in step S3 by 0.95.

[0182] Finally, the molar composition z of the five components of the feed was obtained. feed Specifically as follows: CO2 feed mole fraction = 0.0500, C1 feed mole fraction = 0.0946, C2-C6 feed mole fraction = 0.3591, C7-C 12 Feed mole fraction = 0.3445, C 13+ Feed mole fraction = 0.1519; normalized checksum = 1.0001.

[0183] The above feed composition is used for all subsequent PR-EOS phase equilibrium calculations (K-value iteration and Rachford-Rice flash evaporation) and MMP molar weighted calculations.

[0184] (2) Based on the nanoconfining effect, the critical parameters of each component are corrected and the phase equilibrium calculation of the PR equation of state is performed: S1. Obtain the critical temperature of the injected fluid and each pseudo-component in step (1) under bulk conditions. Critical pressure Eccentricity factor and molecular dynamics diameter .

[0185] The relevant parameters of each component are shown in Table 1.

[0186] Table 1

[0187] S2. Perform nano-confined critical parameter correction: First, based on the characteristic pore size of the target reservoir rock and the molecular dynamic diameter of each component Differential determination of the number of adsorption layers of each component on the pore wall The number of adsorption layers for each component is set according to the molecular diameter: n ads(CO2) =2, n ads(C1) =1 (methane molecules have the smallest size and only require a single layer for adsorption). n ads(C2-C6) =2, n ads(C7-C12) =3 (octane molecules are large and require three layers of adsorption). n ads(C13+) =4 (a hexadecane molecule requires a maximum of four adsorption layers).

[0188] Then, the adsorption layer thickness of each component is calculated according to formula (1). t ads，i and effective aperture d eff,i : (1); Characteristic pore size of target reservoir rock =10nm=1× m.

[0189] And calculate the confinement ratio of each component. x i , .

[0190] Adsorption layer thickness of each component t ads，i Effective aperture d eff,i and the limit ratio x i The calculation results are shown in Table 2.

[0191] Table 2 Pore parameters of each component

[0192] Finally, the nanoconfined critical parameter is corrected according to formula (2): (2).

[0193] Corrected nanocritical parameters , As shown in formula (3): (3).

[0194] In this embodiment, the critical parameter offset correction function is used. f ( x i )=-0.9409 x i +0.2415 .

[0195] Calibration coefficients in the correction function a 1. a 2. The molecular simulation results of Lennard-Jones fluid based on the Zarragoicoechea-Kuz nanoconfined model were calibrated to obtain: a 1 = -0.9409 a 2 = 0.2415.

[0196] The specific process for correcting the nanoconfined critical parameters of each component is as follows: ① Component CO2: f (0.03802) = -0.9409 × 0.03802 + 0.2415 × 0.03802 2 =-0.03577+0.00035=-0.03542.

[0197] =304.13×(1-0.03542)=293.36K.

[0198] =7.377×(1-0.03542)=7.1157MPa.

[0199] ②Component C1: f (0.04066) = -0.9409 × 0.04066 + 0.2415 × 0.04066 2 =-0.03826+0.00040=-0.03786.

[0200] =190.56×(1-0.03786)=183.35K.

[0201] =4.599×(1-0.03786)=4.4249MPa.

[0202] ③ Components C2-C6: f (0.06250) = -0.9409 × 0.06250 + 0.2415 × 0.062502 =-0.05881+0.00094=-0.05787.

[0203] =469.70×(1-0.05787)=442.53K.

[0204] =3.370×(1-0.05787)=3.1750MPa.

[0205] ④ Component C7-C 12 : f (0.10129) = -0.9409 × 0.10129 + 0.2415 × 0.10129 2 =-0.09528+0.00248=-0.09280.

[0206] =568.70×(1-0.09280)=515.92K.

[0207] =2.490×(1-0.09280)=2.2589MPa; ⑤ Component C 13+ : f (0.26563) = -0.9409 × 0.26563 + 0.2415 × 0.26563 2 =-0.24994+0.01704=-0.23290.

[0208] =723.00×(1-0.23290)=554.62K.

[0209] =1.400×(1-0.23290)=1.0739MPa.

[0210] The percentage shift of the critical temperature for each component is as follows: The critical temperature shift for CO2 is -3.54%, for C1 it is -3.79%, for C2-C6 it is -5.79%, and for C7-C... 12 Critical temperature deviation -9.28%, C 13+ Critical temperature deviation -23.29%.

[0211] It is evident that the larger the molecule and the smaller the effective pore size, the more significant the confinement shift, accurately reflecting the physical laws of the nano-confinement effect.

[0212] Eccentricity factor oh It remains unchanged under nanoscale confinement.

[0213] S3, PR-EOS phase equilibrium calculation: First, the nano-corrected critical temperature obtained in step S2 is... and nano-corrected critical pressure Substitute the parameters into the PR equation of state to calculate the phase thermodynamic parameters of each component.

[0214] The following uses C7-C 12 Using octane as an example, we will fully demonstrate the calculation process of PR-EOS parameters.

[0215] The component C7-C obtained in step S2 12 Nanoscale modified critical temperature =515.92K, nano-corrected critical pressure =2.2589× Substituting Pa into the Peng-Robinson (PR) equation of state Calculate component C7-C 12 The phase thermodynamic parameters.

[0216] also, oh =0.3996 (unchanged) T =388.15K, R =8.314 J / (mol·K).

[0217] .

[0218]

[0219] =1.4773× m 3 / mol.

[0220] Component C7-C 12 Comparison temperature: T r = =0.7523.

[0221] Bulk relative temperature (bulk reference value without nanoscale confinement correction) T r,bulk = T / T c,bulk =388.15 / 568.70=0.6826. T c,bulk The original bulk critical temperature (568.70 K) before correction.

[0222] As can be seen, the nano-correction increased the comparison temperature from 0.6826 to 0.7523, an increase of about 10.2%. The physical meaning is that the confinement effect makes the component closer to its critical point (more prone to vaporization) at the same temperature.

[0223] oh =0.3996≤0.49, therefore the first set of empirical correlations in formula (6) is selected: k i =0.37464 + 1.54226 × 0.3996 - 0.26992 × 0.3996 2 =0.37464+0.61629-0.04310 =0.9478.

[0224] Component C7-C 12 Temperature correction factor The calculation is as follows: =[1+0.9478×(1-0.8674)] 2 =[1+0.9478×0.1326] 2 =[1+0.1257] 2 =1.1257 2 =1.2672.

[0225] Component C7-C 12 Temperature-corrected gravitational parameters The calculation is as follows: =3.7246×1.2672=4.7201.

[0226] Comparison of body phases : PR equation of state component C7-C 12 Gravitational parameters The calculation is as follows: =0.45724×8.314 2 ×568.70 2 / (2.490× =3.7869.

[0227] Bulk temperature correction factor The calculation is as follows: =[1+0.9478×0.1738] 2 =1.16472 =1.3565.

[0228] Bulk gravitational parameters The calculation is as follows: =3.7869×1.3565=5.1370.

[0229] Nanoscale corrections affect gravitational parameters The gravitational parameter decreased from 5.1370 to 4.7201, a drop of approximately 8.1%. This decrease in gravitational parameter implies a weakening of the intermolecular attraction in the liquid phase, which is entirely consistent with the physical expectation that confinement effects suppress liquid phase stability.

[0230] The PR-EOS parameters of each component are summarized in Table 3 below.

[0231] Table 3. Parameters of five-component nano-modified PR-EOS

[0232] Calculate binary interaction parameters k ij : The baseline values ​​for the CO2-hydrocarbon system are determined by the logarithmic correlation of carbon atoms: k ij,base =0.0289+0.0421×ln(n c ).

[0233] k ij Temperature correction factor T corr = .

[0234] exist T = T ref When the value is 388.15K, the correction factor is always 1.0.

[0235] Each CO2-hydrocarbon pair k ij The value is: k ij (CO2-C1)=max(0.02,min(0.0289+0.0421×ln(1),0.15)) =max(0.02,min(0.0289,0.15)) =0.0289; k ij(CO2-C2-C6)=max(0.02,min(0.0289+0.0421×ln(5),0.18)) =max(0.02, min(0.0289+0.0421×1.6094, 0.18)) =max(0.02,min(0.0967,0.18)) =0.0967; k ij (CO2-C7-C 12 )=max(0.02, min(0.0289+0.0421×ln(8), 0.25)) =max(0.02, min(0.0289+0.0421×2.0794, 0.25)) =max(0.02,min(0.1164,0.25)) =0.1164; k ij (CO2-C 13+ )=max(0.02, min(0.0289+0.0421×ln(16), 0.35)) =max(0.02, min(0.0289+0.0421×2.7726, 0.35)) =max(0.02,min(0.1456,0.35)) =0.1456.

[0236] Various hydrocarbon-hydrocarbon systems k ij Based on the carbon number difference calculation, typical values ​​are as follows: k ij (C1-C2-C6) approximately 0.005; k ij (C1-C7-C 12 Approximately 0.010; k ij (C1-C 13+ Approximately 0.020; k ij (C2-C6-C7-C 12 Approximately 0.005; k ij (C2-C6-C 13+ Approximately 0.015; k ij (C7-C 12 -C13+ Approximately 0.008.

[0237] Construct a complete 5x5 symmetry k ij The matrix will be substituted into formula (7) (van der Waals mixing rule) to calculate. a mix , b mix Used for subsequent fugacity coefficient and K Value iterative calculation.

[0238] For multi-component mixtures, the classical mixing rules of van der Waals are used, and the gravitational parameters of the mixture are calculated according to formula (7). a mix Harmony volume parameters b mix .

[0239] The feed composition z in this embodiment feed (CO2=0.0500, C1=0.0946, C2-C6=0.3591, C7-C 12 =0.3445, C 13+ =0.1519), T Taking K = 388.15 as an example, the gravitational parameters of the mixture are calculated according to formula (7). a mix Harmony volume parameters b mix : Table 4 Components a i (T), b i

[0240] Step Two: k ij Matrix (CO2-hydrocarbons are calculated according to formula (8), and the hydrocarbon-hydrocarbon relationships are approximately taken as 0): k ij (CO2,Cl) = 0.0289 k ij (CO2,C2-C6)=0.0967, k ij (CO2,C7-C 12 ) = 0.1164, k ij (CO2,C 13+ )=0.1457 Step 3: Substitute into formula (7) to calculate amix : Diagonal term (i=j): ; Off-diagonal term (i<j): Table 5

[0241] After summing off-diagonal terms (i<j)=1.22041, calculate according to formula (7): ; .

[0242] Step 4: Substitute into b mix =Σ x i b i Calculate b mix : ; .

[0243] Calculation results: a mix =3.5915 Pa·m 6 / mol 2 , b mix =1.3788×10 -4 m 3 / mol. Substituting both into the aforementioned formula (7) allows entering the cubic form of the PR equation of state to solve for the compression factor Z and subsequent calculation of the fugacity coefficient.

[0244] Then, based on the phase thermodynamic parameters of each component obtained from the calculation of the PR equation of state, iteratively solve for multi-component phase equilibrium.

[0245] K Calculation of (Wilson) initial values for value iteration ( P =25MPa=25× Pa):

[0246] The Wilson initial values clearly reflect the volatility gradient of each component: is approximately 3783 times that of . This value gradient generated by pure thermodynamic calculation is K the physical basis for the subsequent design of differential mass transfer parameters.

[0247] The molar composition of the five components obtained in step (1) S4 is z feed For the feed composition, a successive replacement method is used. K The PR-EOS fugacity coefficients are solved iteratively.

[0248] In each iteration, the Rachford-Rice flash equation is called to solve for the gas phase fraction. V .

[0249] Under global feeding (CO2 accounts for only 5%), P Under the condition of 25MPa, the discriminant function value is: =0.0500×1.415+0.0946×3.110+0.3591×0.0495+0.3445×0.00762+0.1519×0.000822=0.07075+0.29427+0.01778+0.00263+0.00012=0.3856<1, the system is biased towards the liquid phase overall.

[0250] Convergence of each component K The value tends to 1 after 7 iterations (because the CO2 mole fraction is small, the system is in a near-miscible state).

[0251] Meanwhile, the minimum miscibility pressure of each component is estimated using the power-law correlation formula for carbon number. MMP i : Fahrenheit temperature =239°F.

[0252] Temperature correction factor The calculation is as follows: Cut off to the upper limit of 1.30.

[0253] Each component The calculation is as follows: =0.52× ×1.30× Pa.

[0254] =0.52× ×1.30× Pa = 0.52 × ×1.30× Pa = 0.68 MPa.

[0255] =0.52× ×1.30× Pa = 0.52 × 14.23 × 1.30 × =9.62MPa.

[0256] =0.52× ×1.30× Pa = 0.52 × 30.91 × 1.30 × =20.90MPa.

[0257] =0.52× ×1.30× Pa = 0.52 × 97.01 × 1.30 × =65.58MPa.

[0258] by MMP exp =23.73MPa was used as the anchoring value to calculate the correction factor: λ= MMP exp / MMP calc (C7-C 12 =23.73 / 20.90 = 1.1354 The correction factor is applied proportionally to the estimated MMP values ​​of the remaining components to obtain the corrected minimum miscibility pressures for each component: MMP C1,corrected =0.68 × 1.1354 = 0.77 MPa MMP C2-C6,corrected =9.62 × 1.1354 = 10.92 MPa MMP C13+,corrected =65.58 × 1.1354 = 74.46 MPa MMP C7-C12,corrected =23.73MPa (anchor point itself, consistent with experimental calibration value) Subsequent calculations of the miscibility factor for each component (misc i All of these are based on the corrected MMP values ​​mentioned above, see replacement position B for details.

[0259] Miscibility factor of each component Depend on = / The result is obtained through the following piecewise function.

[0260] (17); Among them, P injCalculations were performed based on the corrected MMP value at replacement position A under a pressure of 25 MPa: r P(C1) =25 / 0.77=32.5, which is much greater than 1, misc C1 =1.0; r P(C2-C6) =25 / 10.92=2.29, greater than 1, misc C2-C6 =1.0; r P(C7-C12) =25 / 23.73=1.054, greater than 1, misc C7-C12 =1.0; r P(C13+) =25 / 74.46=0.336, which is less than 0.5, misc C13+ =0.1×0.336 / 0.5≈0.0672.

[0261] (3) Establish a time-varying linkage mechanism between reservoir temperature and fluid physical parameters: In this embodiment, the pseudo-component division and the various physical property parameters in Tables 1 and 2 described in steps S1-S4 are all at the reference temperature. T base Measured or calibrated at 388.15K. (When simulating operating temperature) T Deviation T base At that time, the relevant physical property parameters should be corrected for time-varying characteristics as follows. T =420K (relative) T base Taking an elevation of 31.85K (simulating near-wellbore heating conditions) as an example: Binary interaction parameters: CO2-C7-C 12 For example, according to T corr =1-ck ij ×( T - T ref ), take ck ij =0.0015、 T ref = T base =388.15K, then T corr (420K) = 1 - 0.0015 × (420 - 388.15) = 1 - 0.0015 × 31.85 = 0.9522; Correspondingly, k ij ( T=420K)= k ij ( T base )× T corr (420K) / T corr (base) = 0.11644 × 0.9522 / 1.0 = 0.11086, meaning that the increase in temperature causes CO2-C7-C 12 As the interaction between the two elements increases, the binary interaction parameter decreases accordingly.

[0262] Henry's constant: in C7-C 12 For example, take the enthalpy of dissolution Δ of the component. H diss =-15kJ / mol (exothermic dissolution), according to calculate: , ,Right now The decrease in the Henry constant indicates that the solubility of CO2 in this component is relatively enhanced after the temperature increases.

[0263] CO2 density and viscosity: experimentally measured temperature T CO2 injection density and lattice kinematic viscosity at 420K are used to replace the values ​​in the table. T base The corresponding values ​​are used for the mass transfer rate coefficient, kinematic viscosity, and diffusion coefficient, respectively, in the same manner. T The experimental data or time-varying correction formulas at 420K will be updated, but will not be elaborated on here. After updating the above physical property parameters, substitute them back into steps (2), (4), (5), and (6) to obtain the results. T Complete simulation results under a 420K operating condition are used for comparison. T base A comparative analysis of operating conditions reveals the impact of temperature on CO2 displacement.

[0264] (4) Systematically calibrate LBM mesoscopic simulation parameters based on phase equilibrium calculation results: Calibration Path 1 (based on the comparison temperature of each component) T r As input, the pseudopotential feature density is calibrated. r 0,σ ): Nano-correction of each component T r : CO2: =1.0×[0.8+0.2×min(1.3231,1.5)] =0.8 + 0.2 × 1.3231 =0.8+0.2646 =1.0646; C1: =1.0×[0.8+0.2×min(2.1170, 1.5)] =0.8 + 0.2 × 1.5000 =0.8+0.3000 =1.1000; C2-C6: =1.0×[0.8+0.2×0.8771] =0.8+0.1754 =0.9754; C7-C 12 : =1.0×[0.8+0.2×0.7523] =0.8+0.1505 =0.9505; C 13+ : =1.0×[0.8 + 0.2×0.6999] =0.8+0.1400 =0.9400.

[0265] The final pseudopotential characteristic density after correction by the nanometer wall using formula (20) ranges from supercritical CO2 and C1 (approximately 1.05-1.10, gas phase characteristic) to liquid C 13+ (Approximately 0.93, liquid phase characteristics) exhibits a continuous transition, naturally reflecting the phase differences of each component.

[0266] Calibration path two (misc to) G (Calibrating interaction strength parameters).

[0267] Substitute the miscibility factors for each component: =0.02 + 0.2 × (1 - 1.0) = 0.020; C1 is completely miscible. G At extremely low values, the interface between CO2 and Cl tends to disappear.

[0268] =0.06 + 0.3 × (1 - 1.0) = 0.060; C2C6 is completely miscible. G It's very small.

[0269] =0.35 + 0.4 × (1 - 1.0) = 0.35; C7C 12 Even at 25 MPa, it was already miscible, using a reference... .

[0270] =0.30+2.0×(1-0.0672)=0.30+2.0×0.9328=0.30+1.848=2.166; C 13+ Non-mixing Take a large value to maintain a clear gas-liquid interface.

[0271] / =2.166 / 0.35=6.19≥5, which meets the constraint requirement of the differentiated interface regulation target of CO2-light components being nearly miscible and CO2-heavy components being immiscible.

[0272] Calibration path three (μ to τ calibration relaxation time).

[0273] First, establish a dimensionless unit system: Reference length =1× m, reference density =843.5kg / m 3 (Experimental crude oil density, corrected for temperature expansion), reference pressure = P reservoir =55.07× Pa, reference speed Approximately 255.5 m / s, reference kinematic viscosity Approximately 2.555× m 2 / s.

[0274] Lattice kinematic viscosity of each component Converted from the ratio of physical viscosity to reference viscosity: v CO2 =μ CO2 / ( × Approximately 0.050; corresponding to τ CO2 =3×0.050+0.5=0.650.

[0275] v C1 =μ C1 / (rho C1 × Approximately 0.075; corresponding to τ C1 =3×0.075+0.5=0.725.

[0276] v C2-C6 =μ C2C6 / (rho _C2C6 × Approximately 0.048; corresponding to τ C2-C6 =3×0.048+0.5=0.644).

[0277] v C7-C12 =μ C7C12 / (rho _C7C12 × Approximately 0.080; corresponding to τ C7-C12 =3×0.080+0.5=0.740.

[0278] v C13+ =μ C13+ / (rho C13+ × Approximately 0.150; corresponding to τ C13+ =3×0.150+0.5=0.950.

[0279] The effective relaxation time of near-wall lattice points increases when multiplied by the wall reinforcement factor, and under 10 nm pore size conditions, C7-C 12 The wall viscosity enhancement coefficient is 1.62 (corresponding to a wall viscosity τ of approximately 0.86), C 13+ The wall viscosity enhancement coefficient is 1.78 (corresponding to a wall τ of approximately 1.40), reflecting the additional inhibition effect of nanoconfinement on fluid flow near the wall.

[0280]

[0281]

[0282]

[0283] The nanowall viscosity enhancement effect is caused by step (2). and wall interaction coefficient Regulation.

[0284] (5) Construct the MRT-MCMP lattice Boltzmann multiphase multicomponent flow simulation framework and embed the differential mass transfer model: First, establish the D2Q9-MRT basic LBM model. The two-dimensional mesh size is lx=300 multiplied by ly=150.

[0285] The nine discrete velocity directions of the D2Q9 model are: e0 = (0, 0); e1, e3 = ( 1, 0); e2, e4 = (0, 1); e5=(1,1); e6=(-1,1); e7=(-1,-1); e8=(1,-1).

[0286] Weighting coefficients: w0 = 4 / 9, w 1-4 =1 / 9, w 5-8 =1 / 36.

[0287] Grid speed of sound .

[0288] The transformation matrix M is a 9x9 orthonormal matrix, and the inverse matrix M -1 It can be obtained directly by taking its transpose matrix (i.e., M). -1 =M T Since M is an orthogonal matrix, this avoids the extra computational overhead and errors caused by numerical inversion.

[0289] Each component uses an independent distribution function array (σ corresponds to CO2 gas phase, C1, C2-C6, C7-C respectively) 12 C 13+ The relaxation time calibration path for each of the six components (including dissolved CO2) is determined.

[0290] diagonal relaxation matrix .

[0291] Macroscopic density of each component at each grid point Macroscopic velocity is determined by momentum and external force.

[0292] Mixing equilibrium speed Used for calculating the equilibrium distribution function of each component.

[0293] In each calculation step, the pseudopotential value of each component under the current density field is first calculated. ( x ).

[0294] ; in, It has been calibrated by step (3).

[0295] Shan-Chen interaction forces between each component pair According to the calibration G Parameters and Equation (27). Fluid-solid wall interaction forces According to the calibration G s Calculate the parameters and equation (28). Total external force for each component. .

[0296] Then, perform MRT collision calculations: transform the distribution functions of each component from velocity space to moment space, perform relaxation collisions (towards equilibrium state) in moment space, superimpose the external force terms in moment space, and then inversely transform back to velocity space.

[0297] Next, the migration step is performed—the distribution function after the collision is propagated to the adjacent cells along the discrete velocity directions.

[0298] Finally, the boundary condition is applied—CO2 is injected into the inlet at a constant lattice velocity ve0 (the inlet CO2 density is set to the injection density rhog). c Approximately 0.65, the inlet density of each oil phase component is set to a minimum value amin= (This indicates that only pure CO2 is injected at the inlet), the outlet is a convective boundary condition, and the upper and lower walls are half-step bounce solid boundaries.

[0299] The initial densities of each component must satisfy the conservation of total oil lattice density: characteristic densities are set as C1=0.25 (small characteristic density for supercritical methane), C2-C6=0.50 (light liquid hydrocarbons), and C7-C... 12 =0.62 (medium octane), C 13+ =0.80 (heavy hexadecane).

[0300] Initial lattice density of each component = mass fraction × characteristic density × density correction factor. Density correction factor = total base oil density of two components (0.66 × 0.60 + 0.34 × 0.80 = 0.668) / total original density of four components (0.015 × 0.25 + 0.2561 × 0.50 + 0.3889 × 0.62 + 0.34 × 0.80 = 0.6449) = 0.668 / 0.6449 = 1.0358. The initial lattice density of each component is obtained as follows: C1=0.015×0.25×1.0358=0.0039; C2-C6=0.2561×0.50×1.0358=0.1327; C7-C 12 =0.3889×0.62×1.0358=0.2499; C 13+ =0.34×0.80×1.0358=0.2817.

[0301] The total is 0.6682, which is precisely consistent with the two-component baseline, ensuring the continuity of the denominator in the recovery rate calculation and its comparability with experimental data.

[0302] In mass transfer calculations, for each non-solid grid point within the computational domain, mass transfer determination and operations are performed sequentially in the order of M1 to M2 to M3 to M4. Gradient numerical settings for differentiated mass transfer parameters are then applied. All parameters below follow a strict gradient from C1 to C13+: Mass transfer mechanism M1: Calculate the effective Henry's constant, equilibrium dissolved concentration, and mass transfer driving force; perform dissolution or desorption according to the direction of the driving force by adjusting the CO2 gas phase distribution function g. α and the dissolved CO2 distribution function h _α and concentration distribution function g k_α To achieve mass transfer. Henry's constant. : H C1 =0.020 (Methane hardly dissolves CO2 because supercritical methane is itself in the gas phase). H C2C6 =0.18 (moderate solubility in light hydrocarbons). H C7C12 =0.25 (based on medium hydrocarbons - this value is derived from P) asp =14 MPa experimental calibration data anchoring, Henry light =0.25×(P / 22)^0.85), H C13+ =0.12 (Heavy hydrocarbons have smaller gaps between macromolecules and are less capable of accepting CO2 molecules than medium hydrocarbons).

[0303] Mass transfer mechanism M2 (component-by-component): Check whether the CO2 concentration and the concentration of the component meet their respective trigger thresholds. If they do, calculate the extraction rate and apply double limiting. Then, subtract the extraction amount from the component distribution function and add an equal amount of CO2 gas phase distribution function.

[0304] Mass transfer mechanism M3 (C only) 13+ ): Execute at a trigger threshold higher than M2 and a rate much lower than M2. CO2 concentration extraction threshold: =0.06 (lowest threshold - C1 supercritical state, can be extracted at almost any CO2 concentration). =0.10, =0.15 (benchmark), =0.25 (highest threshold - C is only triggered in areas with high CO2 enrichment) 13+ extraction).

[0305] Extraction rate coefficient : k extractC1 =0.008 (C1 has extremely high volatility and is the fastest to extract). k extractC2C6 =0.006 (rapid extraction of light hydrocarbons) k _extractC7C12 =0.006 (benchmark) k extractC13+ =0.0005 (approximately 1 / 12 of the baseline, with extremely limited extraction).

[0306] Volatile factors : volume C1 =2.0, volume C2C6 =1.5, volume C7C12 =1.0 (baseline). This reflects the volatility of C1 compared to C7-C. 12 The physical reality is twice the benchmark value.

[0307] Mass transfer mechanism M4: Check whether the total CO2 concentration exceeds the solvent effect trigger threshold. If it does, calculate the separation rate of each component based on the differential selectivity factor and limit it, and then perform mass transfer.

[0308] All mass transfer operations achieve mass transfer and conservation by increasing, decreasing, or displacing the corresponding component's lattice distribution function.

[0309] M4 Solvent Effect Selectivity Factor : sel C1 =1.5, sel C2C6 =1.2, sel C7C12 =1.0 (baseline) sel C13+ =0.20. This demonstrates that CO2 solvent has a much stronger ability to dissolve and carry light, small molecules than heavy, large molecules.

[0310] Maximum transmission quality limit (max) _rate_σ M2 C1 Approximately 0.010, M2 C2C6 Approximately 0.008, M2b C7C12 Approximately 0.008, M3 C13+ Approximately 0.001.

[0311] Meanwhile, the execution parameters for the partitioned local phase balance dynamic feedback coupling are: once every 200 steps; the number of partitions Nx multiplied by Ny = 5 multiplied by 3 = 15 partitions.

[0312] K-value smoothing coefficient =0.30; G-parameter smoothing coefficient =0.15.

[0313] To avoid instability in flash evaporation calculations due to excessively low local mass, a minimum zone mass threshold (min) can be set. mass_per_zone =0.5, partitions with quality below this threshold skip the local flash evaporation calculation in this round to save computational costs.

[0314] Every 200 steps, 15 partitions are traversed, and the local composition and local pressure of each partition are statistically analyzed. PR-EOS-K value iteration and Rachford-Rice flash evaporation are then performed to obtain the local K value and local enrichment coefficient. .

[0315] based on Update target G parameters: If C7-C 12 A higher enrichment coefficient reduces The target value promotes miscibility; if C 13+ If the enrichment coefficient is low (much less than 1), it will remain or increase. The target value is maintained by the interface. The global K-value field and G parameters are updated with a smoothing coefficient to complete the feedback loop.

[0316] The total number of steps t in the simulated main loop max =20,000 steps. Data is recorded every 100 steps—including total pore mass of each component, cumulative mass transfer of each mass transfer mechanism, expansion factor SF, and sweep efficiency. Interfacial tension (IFT), K-value field mean, and four-mechanism decomposition results.

[0317] (6) Quantitative decomposition of mass transfer mechanism based on multiple independent physical sources and self-verification of mass conservation closure: During the simulation, a four-mechanism decomposition calculation was performed every 100 steps. The simulation was completed at a certain intermediate point (injected pore volume multiple PV approximately 0.5, total recovery rate...). Taking approximately 28.5% as an example, the complete calculation process is shown below: ① Calculate the miscibility factor .current r P = P inj / MMP exp =25 / 23.73=1.054>1.0, misc P =1.0 indicates that the system has met the thermodynamic miscibility pressure condition under the current injection pressure.

[0318] Current interface tension IFT current (Depend on G fg × c s 2 ×ρ g,avg ×ρ o,avg (Estimation) and Initial IFT ratio r IFT Approximately 0.6, corresponding to dynamic interface items. miscellaneous IFT≈0.12 indicates that although the pressure conditions meet the miscibility requirements, the CO2 sweep is not yet sufficient (current injected pore volume multiple). PV ≈0.5, ripple efficiency E v (≈0.35), the actual interface state of the flow field has not yet evolved to a fully miscible state.

[0319] In summary, miscellaneous overall = w 1× miscellaneous P + w 2× miscellaneous IFT =0.80×1.0+0.20×0.12=0.824, which reflects the upper limit of the static miscibility potential that the system can reach after full saturation.

[0320] Since CO2 diffusion is not yet sufficient at this moment, it is necessary to introduce a concept related to diffusion efficiency. E v The associated secondary reduction avoids overestimating the actual degree of miscibility achieved by the system in the early stages when the impact is insufficient. After reduction, miscibility... overall At this moment ( PV ≈0.5) takes a value of approximately 0.15, which is used for the specific calculation of each sub-item of mechanism decomposition in step (6).

[0321] ② Calculate the inflation contribution. SF current Approximately 1.15 (the liquid phase expanded by about 15% due to CO2 dissolution). Approximately 0.35. =0.35×(1.15-1) / 1.15×0.55=0.35×0.1304×0.55=0.0251, approximately 2.51%. =M1 net cumulative dissolution amount / Approximately 0.028 / 0.668, approximately 4.19%. Taking the larger value of 4.19%, and comparing it with the expansion contribution from the previous time step (e.g., 3.8%), we obtain the larger value (monotonically non-decreasing constraint). Approximately 4.2%.

[0322] ③ Calculate the contributions of immiscible extraction and miscibility. Cumulative mass transfer (up to the current time): Total mass transfer of the light component. =M2 C1cum +M2 C2C6cum +M3 C7C12cum +M4 cum ×f light_init Approximately 0.085 (lattice units). Total mass transfer of heavy components. =M3 C13p_cum +M4 cum×f heav nit Approximately 0.018. Currently... Take 0.15 (because it is in the transition region from immiscible to near-immiscible). =0.15 2 =0.0225. The contribution of light miscible phase to dissolution is approximately 1.91% (0.085 × 0.15 / 0.668); the contribution of light immiscible phase to extraction is approximately 10.81% (0.085 × (1 - 0.15) / 0.668).

[0323] The contribution of heavy miscible phase mixing is approximately 0.06% (0.018 × 0.0225 / 0.668); the contribution of heavy immiscible phase extraction is approximately 2.63% (0.018 × (1 - 0.0225) / 0.668). =1.91% + 0.06% = 1.97%. =10.81% + 2.63% = 13.44%. It can be seen that the mass transfer of the heavy component is almost completely (97.75%) attributed to the immiscible phase extraction, with only 2.25% belonging to the miscible phase, accurately reflecting C... 13+ The physical reality is that extraction is mainly achieved through interfacial diffusion (rather than miscible mixing) under immiscible conditions.

[0324] ④ Calculate the displacement contribution. .

[0325] ⑤ Self-check for closure error. =8.89% + 4.2% + 13.44% + 1.97% - 28.5% = 28.50% - 28.5% = 0.00%. The absolute value of the closure error is much smaller than... This verifies that the decomposition results satisfy the law of conservation of mass, and the physical rationality of the four-mechanism quantitative decomposition algorithm at this moment.

[0326] Microfluidic experiments were conducted using the same parameters to verify the dynamic expansion factor SF, and the fitting result R was obtained. 2 =0.98, such as Figure 4 As shown.

[0327] Simulations up to the final time step (t=20000 steps, injected PV approximately 3.83) showed a total recovery rate of approximately 77.12%. The recovery rates of each component level exhibited a significant decreasing gradient—C1 had the highest recovery rate (96.25%, due to the near-complete removal of supercritical C1), followed by C2-C6 (94.94%), and then C7-C... 12 In the middle (91.97%), C 13+ The lowest (55.59%, mainly from displacement and carryover, and extraction at the immiscible interface; miscibility contribution was only 9.70% – subject to squared decay correction). _overall )2 inhibition).

[0328] The contribution percentages of the four mechanisms are approximately as follows: displacement and carryover contribute the most (approximately 44.3%), followed by extraction from immiscible interfaces (approximately 28.0%, mainly from C). 13+ Limited extraction), miscible extraction contributes approximately 25.2% (at P=25 MPa, it exceeds C1 / C2-C6 / C7-C 12 Under the conditions of their respective MMP values, the three light and medium components are completely miscible, with swelling contributing approximately 2.4%. Extraction selectivity ratio SR = (M2) _C1 +M2 _C2C6 ) / M3 _C7C12 The simulation exhibited a dynamic evolution—the initial SR was approximately 130 (lighter components were preferentially and rapidly extracted), gradually decreasing to approximately 0.246 as lighter components were depleted, thus validating the differential mass transfer model's ability to finely resolve the timing of component selective extraction. Throughout the simulation, the absolute value of the closure error remained within 10 at each time point. -4 The mass conservation consistency of the four-mechanism decomposition algorithm was verified to be below the order of magnitude (0.00% at the final moment) throughout the entire simulation period.

[0329] (7) After completing all 20,000 simulation steps using the constructed model, the data recorded at each time step, including the density field of each component, cumulative mass transfer, recovery curve, phase saturation, gas-liquid phase composition, K-value mean evolution, contribution ratio of the four mechanisms, and closure error, are visualized using Matlab to generate multi-angle analysis charts of the simulation results. These are then saved as structured MAT data files (.mat format) and Excel data tables (.xlsx format, containing 8 worksheets including recovery and mass balance, evaluation indicators, mass transfer mechanism, phase K-value, miscibility, gas-liquid phase composition, mechanism decomposition evaluation, and component mechanism decomposition matrix), for subsequent analysis and engineering applications. The results are as follows: Figure 5 , Figure 6 As shown.

Claims

1. A simulation method for differential mass transfer in multiphase, multicomponent fluids based on LBM, characterized in that, Includes the following steps: (1) Constructing a multi-component thermodynamic system based on microfluidic experimental data: S1. Obtain basic physical property data of the target reservoir through microfluidic experiments and crude oil component testing; S2. Based on the obtained carbon number distribution data of the reservoir fluid crude oil, the crude oil is divided into different pseudo-components, including methane (C1), light hydrocarbons (C2-C6), and medium hydrocarbons (C7-C6). 12 and heavy hydrocarbons C 13+ ; S3. Obtain the mass fraction of each pseudo-component and calculate the mole fraction of the pseudo-components. Conversion; S4. Set the feed mole fraction for each component; At the same time, a dissolved CO2 phase is added to the original component system to independently track the concentration of CO2 that has dissolved into the liquid phase; (2) Based on the nanoconfining effect, the critical parameters of each component are corrected and substituted into the PR equation of state to perform phase equilibrium iterative calculations to solve the phase equilibrium of the multi-component components. K Value; estimate minimum miscibility pressure; calculate miscibility factor; (3) Establish a time-varying linkage mechanism between reservoir temperature and fluid physical parameters; including CO2 density, CO2 viscosity, Henry's dissolution constant, mass transfer rate, viscosity of each component of crude oil, and binary interaction parameters. k ij The physical properties are calculated to be corrected for specific time variations with temperature. (4) Based on the phase equilibrium calculation results, systematically calibrate the LBM mesoscopic simulation parameters; set three specific calibration paths, namely pseudo potential characteristic density, Shan-Chen interaction strength parameter G, and MRT relaxation time τ; (5) Construct the MRT-MCMP lattice Boltzmann multiphase multicomponent flow simulation framework and embed the differential mass transfer model: S1. Based on the LBM mesoscopic simulation parameters calibrated in step (4), construct a multi-relaxation time multi-component Shan-Chen pseudopotential lattice Boltzmann simulation framework; S2. In each calculation step of the LBM main cycle, a differentiated multi-mechanism mass transfer model consisting of four mass transfer mechanisms—dissolution-desorption, fractional extraction, limited extraction of heavy components, and combined solvent effect desorption—is embedded. The four mass transfer mechanisms are as follows: Mass transfer mechanism M1: CO2 dissolution and desorption mass transfer in the liquid phase, driven by Henry's law; Mass transfer mechanism M2: fractional extraction mass transfer of CO2 for methane, light hydrocarbons, and medium hydrocarbons; targeting C1, C2-C6, and C7-C6 hydrocarbons. 12 Perform extraction independently; Mass transfer mechanism M3: CO2 on heavy hydrocarbon C 13+ Limited extraction mass transfer; Mass transfer mechanism M4: Co-departure mass transfer enhanced by CO2 solvent effect; S3. Introduce a dynamic feedback coupling mechanism based on partitioned local phase balance: every preset number of simulation steps The entire computational domain is divided into N parts along the x and y directions, respectively. x and N y One rectangular partition; For each partition, the current lattice density of each component in all non-solid grid points within that partition is calculated and converted into the local molar composition of the partition. And estimate the local pressure of the zone. ; Composed of local moles in partitions and the estimated local pressure of the zone For input, repeat step (2) completely; Get the local partition K The local gas-liquid composition and local gas-liquid composition are analyzed, and the local enrichment coefficient is finally calculated. S4. Based on the local enrichment coefficient, the global phase equilibrium is... K Value field and interphase interaction intensity parameters G Perform a smooth update; S5. Based on the aforementioned partitioned local phase balance dynamic feedback coupling mechanism, the following time-varying parameter update mechanism within the simulation time step is further established. ①Time-varying update miscibility factor ; ②Time-varying adaptive Shan-Chen interaction strength parameters G During each partition balancing update, the enrichment coefficient of each partition is used as a basis. E σ Adaptive adjustment of CO2-middle-component interaction strength G fg Interaction strength with CO2-heavy components G dg The target value is determined by the smoothing coefficient α. G Asymptotic updates toward the target value, making G The parameters evolve dynamically with the displacement process; ③ The relative contributions of time-varying decomposition micturition extraction and immiscible extraction; The contribution ratios of four mass transfer mechanisms—expansion, immiscible interface extraction, miscible extraction, and displacement-carryover—were quantitatively decomposed and continuously updated with the time-varying miscibility factor. (6) Based on the quantitative decomposition and mass conservation of four mass transfer mechanisms, namely expansion, immiscible interface extraction, miscible extraction and displacement and carryover, the closure error is self-verified. (7) Based on the complete model constructed in steps (1) to (6), which covers thermodynamic phase equilibrium, LBM mesoscopic parameter calibration, temperature time-varying correction, differentiated mass transfer mechanism and partitioned feedback coupling, mechanism decomposition and conservation verification, simulation calculation is performed to obtain the spatiotemporal distribution of density fields of each component, cumulative mass transfer time series of each mass transfer mechanism, mass transfer efficiency curves of the overall and component levels, spatiotemporal evolution of saturation of each phase and gas-liquid phase composition, decomposition curves of contributions of four mass transfer mechanisms and closure error verification results, spatiotemporal evolution of phase equilibrium K value field and dynamic change data of key evaluation indicators such as interfacial tension, capillary number, expansion factor and sweep efficiency. The dynamic change data is then graphically presented using visualization software.

2. The simulation method for differential mass transfer of multiphase and multicomponent fluids based on LBM according to claim 1, characterized in that, The corrected critical temperature under confined conditions in step (2) Critical pressure under confined conditions As shown below: ； In the formula, —Components i Critical temperature under bulk conditions; f ( x i —Regarding the limit ratio x i Critical parameter offset correction function; —Components i Critical temperature under confined conditions; —Components i Critical pressure under bulk conditions; —Components i Critical pressure under confined conditions; in, ; In the formula, —Components i The offset of the critical temperature relative to the bulk value; —Components i The offset of the critical pressure relative to the bulk value; —Components i Critical temperature under bulk conditions; —Components i Critical pressure under bulk conditions; f ( x i —Regarding the limit ratio x i Critical parameter offset correction function; a 1. a 2—calibration coefficient; where, a The typical value range for 1 is between -0.95 and -0.

90. a The typical value range for 2 is between 0.20 and 0.25; x i —Components i The limit ratio; In step (2), the multi-component phase equilibrium K The initial value for iteration is obtained by the following formula: ； K i (0) —Components i Phase equilibrium constant K Initial value for iteration; —Components i Critical pressure under confined conditions; P —System pressure; —Eccentricity factor; —Components i Critical temperature under confined conditions; T —Reservoir temperature; The minimum miscibility pressures corresponding to each component in step (2) are estimated using the power-law correlation of carbon number: ； In the formula, MMP i —Components i The corresponding minimum miscibility pressure; the unit in this calculation formula is MPa; —Components i The representative carbon number; f ( T —Temperature correction function; in, f ( T The result is obtained through the following formula: ； In the formula, —Fahrenheit temperature; The calculation of the miscibility factor in step (2) is as follows: First, calculate the injection pressure. The ratio of the minimum miscibility pressure of the component ; Then, the miscibility factor of each component is obtained through the following piecewise functional mapping. ; 。 3. The simulation method for differential mass transfer of multiphase and multicomponent fluids based on LBM according to claim 1, characterized in that, In step (4), the pseudo-potential characteristic density is used as calibration path one, which adopts the comparison temperature. T r,σ Calibration of pseudopotential characteristic density Calculate according to the following formula: ； In the formula, —Pseudopotential feature density; —The baseline density can be taken as 1.0; The 1.5 in the min function is the upper limit for truncation; T r,σ —Comparison temperature; The pseudopotential function adopts an exponential form of the Shan-Chen type: ； In the formula, —The pseudopotential function of component σ; ρ 0,σ —Pseudopotential feature density; ρ— Local lattice density of the component; At the same time, a wall correction item is introduced: ； In the formula, β σ —Dimensionless confinement strength factor; x σ —The confinement ratio of component σ; σ σ —The molecular dynamic diameter of component σ; d eff,σ —Effective pore size of component σ; —The pseudopotential characteristic density corrected for the near-wall region.

4. The simulation method for differential mass transfer of multiphase and multicomponent fluids based on LBM according to claim 1, characterized in that, In step (4), the Shan-Chen interaction strength parameter G is used as calibration path two, which uses the misc misc factor to calibrate the interaction strength parameter. : Component σ and The strength parameter of the Shan-Chen interaction between them is determined by the following formula: ； In the formula, —Components and The baseline value of the interaction strength parameter between them corresponds to the completely miscible phase. G Lower limit; —The adjustable range of the interaction strength parameter corresponds to the immiscible phase. G The upper limit of the increment above the baseline value; —Components The miscibility factor.

5. The simulation method for differential mass transfer of multiphase and multicomponent fluids based on LBM according to claim 1, characterized in that, In step (4), the MRT relaxation time τ is used as calibration path three, which adopts viscosity. Calibration relaxation time ; First, the viscosity of crude oil was measured experimentally. Given the densities of each component, calculate the lattice kinematic viscosity of each component. lattice kinematic viscosity Calculated from physical viscosity using a dimensionless reference viscosity: ； In the formula, — The kinematic viscosity of component σ, dimensionless; —The kinematic viscosity of the components, calculated from the experimentally determined crude oil viscosity and the density of each component, in m. 2 / s); —The reference viscosity is formed by the lattice space step size Δx and the time step size Δt used in LBM simulation; Then, the relaxation time is mapped to the MRT collision operator using the following formula. : ； For the near-wall region, a nano-wall viscosity enhancement effect is introduced, based on the confinement ratio. Regulation: ； In the formula, — The corrected collision relaxation time of component σ at the fluid node near the wall; —The wall viscosity enhancement coefficient of component σ; — The confinement ratio of component σ.

6. The simulation method for differential mass transfer of multiphase and multicomponent fluids based on LBM according to claim 1, characterized in that, Mass transfer mechanism M1 in step (5) S2: For each fluid grid point, the effective Henry's constant is calculated weighted according to the current composition of the local liquid phase: ； in, —Effective Henry constant; —Components The proportion in the current liquid phase; —Components The specific Henry's law of solubility; ρ σ —The local macroscopic density of component σ at this grid point; —Crude oil density at this grid point; The equilibrium concentration is: ； in, —Local CO2 partial pressure calculated from local CO2 gas phase density; —CO2 gas phase density; —The square of the lattice sound speed, with a value of 1 / 3; —Pressure conversion factor; —Local equilibrium dissolution concentration; —Effective Henry constant; C max This is the preset maximum solubility concentration limit; The mass transfer driving force is: ; in, —The driving force of mass transfer, the sign of which determines the direction of mass transfer; ΔC>0 indicates that the local dissolved CO2 concentration is lower than the equilibrium concentration, driving dissolution; ΔC<0 indicates that it is supersaturated, driving desorption; —Local equilibrium dissolution concentration; —Current dissolved CO2 concentration; Mass transfer mechanism M2 in step (5) S2: The extraction rate formula is: ； In the formula, m extract,σ —Extraction mass transfer rate of component σ; t —Current simulation time; —The extraction rate coefficient of this component; —This component is in local phase equilibrium in the partition to which the current grid point belongs. K value; —The volatility factor of this component; —The miscibility factor of this component; ρ σ —The local macroscopic density of component σ; —CO2 gas phase density; —The square of the lattice sound speed, with a value of 1 / 3; Mass transfer mechanism M3 in step (5) S2: The extraction rate formula is: ； In the formula, —C 13+ Extraction rate coefficient of the component; —C 13+ Local phase equilibrium of components K value; —C 13+ The miscibility factor of the components; —C 13+ Local density of components; —CO2 gas phase density; —the square of the speed of sound in a lattice; Mass transfer mechanism M4 in step (5) S2: The local CO2 percentage is: ； In the formula, —Local CO2 percentage; —Local density of CO2 in the gas phase; —Local density of dissolved CO2; —The sum of the densities of all oil phase pseudo-components at this grid point; For C1, C2-C6, C7-C 12 For light and medium components, the separation rate is: ； In the formula, —Component σ M 4. Separation rate; —The fundamental separation coefficient due to solvent effect; —Local CO2 percentage; —Selectivity factor of component σ; —Local density of component σ; —The miscibility factor of component σ; For heavy component C 13+ The separation rate is: ； In the formula, —Component C 13+ of M 4. Separation rate; —The fundamental separation coefficient due to solvent effect; —Local CO2 percentage; —Component C 13+ The selectivity factor; —Local density of component σ; —Component C 13+ The miscibility adjustment factor.

7. The simulation method for differential mass transfer of multiphase and multicomponent fluids based on LBM according to claim 1, characterized in that, The local molar composition of the partition in step (5) S3 The conversion is performed according to the following formula: ； In the formula, —The local molar composition of component σ in this partition; , —Components , Local density at grid point x; , —Components , The molar mass; σ′\sigma' σ′—a dummy variable summing all components (including σ itself) within this partition; Local pressure in the partition Estimate using the following formula: ； In the formula, —The local estimated pressure in this partition; —The average gaseous CO2 density at all non-solid grid points within this partition; —The square of the lattice sound speed, taking a value of 1 / 3; —Pressure conversion factor; The local enrichment coefficient is calculated according to the following formula: ； In the formula, —Local enrichment coefficient of component σ; —The mole fraction of the vapor phase obtained from the local flash evaporation calculation of this zone; —The liquid phase mole fraction calculated from the local flash evaporation of this zone.

8. The simulation method for differential mass transfer of multiphase and multicomponent fluids based on LBM according to claim 1, characterized in that, In step (5) S4, the global phase equilibrium is determined based on the local enrichment coefficient. K Value field and interphase interaction intensity parameters G Perform a smooth update: ； In the formula, — K Value smoothing update coefficient, typical value 0.15-0.30; — G (Parameter smoothing coefficient, typical value 0.15). —The local K value obtained from this partition calculation; —Determined by local enrichment coefficients G Target value for parameters.

9. The simulation method for differential mass transfer of multiphase and multicomponent fluids based on LBM according to claim 1, characterized in that, The time-varying parameter update mechanism within the simulation time step in step (5) S5: ① Miscibility factor Time-varying update: Constructing the miscibility factor , It is obtained by weighted average of two parts: ； In the formula, —Overall miscibility factor of the system; —Static thermodynamic potential term based on the ratio of global injection pressure P to minimum miscibility pressure MMP, P / MMP. —A dynamic interface state item based on the ratio of the current interface tension to the initial interface tension; , —The weighting coefficients of the two items satisfy + =1; in, Based on injection pressure Ratio to minimum miscibility pressure According to the miscibility factor of a single component The same piecewise function is used to calculate: ； In the formula, —The ratio of injection pressure to the minimum miscibility pressure of the component; -based on The static thermodynamic potential term; Based on the ratio of the current interfacial tension to the initial interfacial tension during the simulation process : ； In the formula, —The ratio of the current interface tension to the initial interface tension; -based on Dynamic interface status items.

10. The simulation method for differential mass transfer of multiphase and multicomponent fluids based on LBM according to claim 1, characterized in that, The quantitative decomposition of the contribution ratios of the four mass transfer mechanisms—expansion, immiscible interface extraction, miscible extraction, and displacement-carryover—in step (5) and step (6) is as follows: ① Calculate the contribution of the expansion mechanism Determined by the larger of two independent physical sources, with a monotonically non-decreasing constraint applied: Source A: Inflation Factor Method ； In the formula, SF— Inflation factor; E v — Sweep efficiency; α sw —Expansion efficiency coefficient, typically around 0.55; Source B Net Solubility Method: ； In the formula, — M 1. Net cumulative amount dissolved; dissolution is positive, desorption is negative, and the amounts are accumulated with signs. ； In the formula, EOR —Currently, the expansion mechanism of this step contributes to the recovery rate; EOR —The expansion mechanism in the previous step contributes to the recovery rate; ② Calculate the extraction contribution of immiscible interfaces and the contribution of miscible extraction : First, the cumulative mass transfer of each mass transfer mechanism is classified according to the light and heavy components of the displaced fluid: ； In the formula, m light,MT — Total mass transfer of light components; M 2C1 — C1 extraction amount; M 2C2C6 — Extraction amount of C2-C6; M 2 b C7C12 — C7-C 12 Extraction volume; M 4cumul · f light,init —M 4. Lightweight components; ； In the formula, m heavy,MT — Total mass transfer of heavy components; M 3 C13+ — C 13+ Extraction volume; M 4 cumul · f heavy,init —M 4. Heavy components; Light component mass transfer Direct allocation; nonlinear decay correction for miscibility of heavy components: ; Then, the extraction contribution of the immiscible interface is calculated according to the following formula. and miscible extraction contribution ; ； ； In the formula, EOR mis — Miscible extraction contributes to oil recovery; m light,MT — Total mass transfer of light components; — Mass transfer of light components; m heavy,MT — Total mass transfer of heavy components; EOR immis — Immiscible interface extraction contributes to recovery rate; ③ Determine the contribution of displacement and carryover mechanisms to EOR disp : ； In the formula, EOR disp — Displacement and carryover mechanisms contribute to oil recovery; EOR total — Total recovery rate; EOR swell — Expansion mechanism contributes to oil recovery; EOR mis — Miscible extraction contributes to oil recovery; EOR immis — Immiscible interface extraction contributes to recovery rate; In step (6), the mass conservation closure error is self-verified according to the following formula: ； If | If the result is less than the preset tolerance, it indicates that the decomposition result satisfies the law of conservation of matter, the calculation process of each sub-step is correct, and the physical rationality of the decomposition result is verified.

Images