Coupling numerical simulation method for temporary plugging micro-transformation process of loose sandstone reservoir

Abstract

The invention discloses a coupling numerical simulation method for a temporary plugging micro-transformation process of a loose sandstone reservoir, and the method comprises the following steps: constructing a coupling control model; carrying out shear-induced anisotropic permeability evolution; modeling non-Newtonian rheology of the temporary plugging agent; judging a strain-driven permeability evolution path; carrying out three-field coupling iterative solution; and outputting a result and optimizing a scheme. The method is based on an independently developed multi-field joint solution framework, can be coupled with a common numerical platform through an intermediate interface for operation, and does not depend on a fixed software combination. The method is particularly suitable for offshore heterogeneous reservoirs and can be stably applied under different working conditions such as temporary plugging, capacity expansion, thermal recovery and profile control. Different from an existing one-way driving model, close linkage of the concentration field, the rheological effect and the mechanical response is achieved for the first time, and the simulation precision and decision-making efficiency of reservoir microtransformation are improved.

Description

Technical Field

[0001] This invention belongs to the field of oil and gas field production enhancement and transformation engineering technology, and in particular relates to a coupled numerical simulation method for the temporary plugging micro-reinforcement process of loose sandstone reservoirs. Background Technology

[0002] Loose sandstone-type heterogeneous reservoirs are widely distributed in offshore oil reservoirs such as the Bohai Bay Basin and the East China Sea. During the mid-stage of development, they are prone to problems such as seepage channel solidification and reduced injection-production efficiency. These reservoirs are generally characterized by loose structure, large differences in reservoir properties, easy interconnection of microfractures, and high sensitivity to injection pressure and shear disturbance.

[0003] During long-term water or steam injection, high-permeability zones tend to break through and form strong conductive channels, while low-permeability zones are difficult to utilize effectively, resulting in a continuous decline in injection and production efficiency. To overcome this bottleneck of "short-circuiting high-permeability channels – stagnation in low-permeability zones," engineering projects have gradually introduced temporary plugging agent modification systems with shear expansion capabilities. These systems aim to temporarily seal high-permeability channels and enhance the volumetric modification capacity of shear zones, thereby enabling the utilization of inefficient strata and the reconstruction of heterogeneity.

[0004] Currently, numerical simulation methods for reservoir stimulation applied in the oil and gas field development field generally have many problems. Some applications are limited to single scenarios and do not involve the simulation of reservoir micro-stimulation with temporary plugging-shear expansion-rheological coupling, thus falling outside the scope of reservoir production enhancement. Others have limited coupling, mostly relying on seepage or strain as the sole dominant factor, or using empirical factors (such as fracture coordinate transformation and thermal field multipliers), failing to achieve explicit simulation of plugging agent concentration migration and non-Newtonian rheological behavior. Some output indicators lack engineering relevance and are not suitable for widespread application. Summary of the Invention

[0005] The problem this invention aims to solve is to provide a coupled numerical simulation method for the temporary plugging and micro-remodeling process in loose sandstone reservoirs. This method integrates the shear expansion rheological properties of the plugging agent with a multi-field coupling mechanism, aiming to achieve channel identification, shear enhancement, and permeability tensor evolution simulation under complex geological conditions. This method features high nonlinearity, strong multi-physics coupling, strong shear sensitivity, and strong structural reconstruction capabilities, making it suitable for the evaluation and design of plugging-connection-micro-remodeling processes in heterogeneous offshore reservoirs.

[0006] To solve the above-mentioned technical problems, the technical solution adopted by this invention is: a coupled numerical simulation method for the temporary plugging micro-remodeling process of loose sandstone reservoirs, comprising the following steps:

[0007] S1: Construct a set of governing equations that couple the mechanical field, seepage field, and temporary plugging agent concentration field to dynamically characterize rock mass deformation, seepage field redistribution, and concentration field transport under injection disturbance, thereby driving the response evolution of the entire process of temporary plugging-micro-modification.

[0008] S2: Construct an anisotropic permeability tensor evolution model induced by shear strain, drive the dynamic change of the permeability tensor through the rock plastic strain tensor, and introduce the concentration of temporary plugging agent to correct the permeability;

[0009] S3: Establish a non-Newtonian shear-thinning rheological model for the temporary plugging agent, define the relationship between the shear viscosity of the temporary plugging agent and the shear rate, and couple the effect of the concentration field on the permeability regulation in the seepage calculation;

[0010] S4: Set up permeability evolution path branches driven by the strain state of rock mass units. During the simulation, dynamically identify the dominant response mechanism of the unit according to the stress-strain evolution state of the unit: when the unit is dominated by the temporary plugging agent retention effect, it is classified as a temporary plugging path; when the unit is dominated by the shear deformation expansion effect, it is classified as a shear modification path; select the corresponding permeability evolution model parameter set for different paths and adjust the unit permeability evolution in real time.

[0011] S5: Construct a three-field strongly coupled numerical model of pore pressure field, plastic strain field and permeability tensor field, and use a step-by-step loading and sub-step iteration control strategy to achieve the solution: Under each loading increment, the mechanical field is updated, the permeability field is updated and the seepage field is solved in sequence, and the pore pressure, permeability tensor and concentration variables are exchanged between the mechanical solver and the seepage solver through an intermediate interface, and the iteration is repeated until the solutions of each field converge.

[0012] S6: Outputs spatiotemporal evolution data of multi-dimensional key variables, generates injection pressure difference-injection volume response curves, rupture pressure maps of temporary plugging zones, shear modification radius distribution maps, and shear channel skeleton structure maps as visual evaluation indicators, and optimizes the temporary plugging agent layout scheme and injection construction parameters.

[0013] Furthermore, S1 includes the following steps:

[0014] S11: Establish the stress-strain relationship between the elastic and plastic stages;

[0015] S12: Plastic deformation adopts the Drucker-Prager improved yield criterion and the non-associated flow rule, defining the yield function F and the plastic potential function G, as follows:

[0016] F=(d0-P t tanβ) 2 +q 2 -p′tanβ-d=0

[0017] in: S=σ+p′I

[0018] G=ωd0tanψ+q 2 -p′tanψ=0

[0019] Where: d0 is the equivalent shear strength parameter of the yield surface, obtained by conversion from the rock cohesion-friction angle relationship, and the unit is Pa; P t β is the equivalent compressive yield parameter, used to characterize the compressive strength level of loose sandstone under low confining pressure, in Pa; β is the Drucker-Prager yield surface cone angle, dimensionless; q is the second deviatoric stress invariant, reflecting the shear strength, in Pa; p′ is the effective mean stress, reflecting the confining pressure level, in Pa; d is the yield function constant, used to control the overall position of the yield surface on the p′-q plane; ψ is the expansion angle, the potential surface inclination angle of the non-associated flow method, describing the degree of volume expansion caused by shear deformation, dimensionless; ω is the potential function weighting coefficient, ensuring that the overall dimensions of G are consistent with the yield function, usually determined by material test inversion; I is the unit second tensor.

[0020] S13: Construct the seepage-stress coupling control equation set.

[0021] Furthermore, step S2 includes the following steps:

[0022] S21: Construct the principal strain response model of the anisotropic permeability tensor for absolute permeability. The principal direction response models of the absolute permeability tensor are defined as follows:

[0023]

[0024] in: Initial absolute permeability in all directions, i.e., the inherent permeability before stress disturbance, in meters. 2 In the above formula, i = 1, 2, 3 correspond to the three principal directions in the principal stress coordinate system; a is the shear-induced permeability enhancement coefficient, i.e., the shear dilatation sensitivity coefficient, which reflects the degree of simultaneous increase in permeability in all directions under the action of strain in the principal direction, and its unit is consistent with the change in permeability, i.e., m. 2 b. Cross-coupling enhancement coefficient, describing the off-diagonal coupling permeability enhancement effect caused by the rearrangement of asymmetric particle structure and pore connectivity, in units of m. 2 The tensor structure reflects the enhanced lateral conduction capability under anisotropic perturbations; ε i The principal strain in each direction is the equivalent axial strain in the corresponding principal direction, and is dimensionless; in the above formula, i = 1, 2, 3 are consistent with the direction of the principal stress.

[0025] S22: The shear-induced enhancement coefficient λ is evolved, and its evolution form is as follows:

[0026] a=λa0, b=λb0

[0027] λ=0.5(1-h)λ min+0.5(1+h)λ max

[0028]

[0029] Where: a0 and b0 are the initial enhancement coefficients, corresponding to the shear-induced enhancement level before plastic disturbance, and are dimensionless; λ min , λ max The lower and upper limits of shear enhancement control the minimum and maximum values ​​of shear-induced permeability enhancement, respectively, and are dimensionless. The equivalent plastic strain is defined as: The cumulative inelastic deformation of the material during loading is dimensionless; ξ is the shear-induced reinforcement initiation point (e.g., 1%), dimensionless; m is a parameter that adjusts the reinforcement steepness (empirical value 40–60); h is an intermediate control variable, dimensionless.

[0030] S23: Construct a relative penetration rate model, which is expressed using the Corey empirical model:

[0031]

[0032] Where: k rw The relative permeability of the aqueous phase is a dimensionless parameter (0–1), reflecting the effective seepage capacity of the aqueous phase at a given water saturation level; k ro The relative permeability of the oil phase is a dimensionless parameter (0–1), reflecting the effective flow capacity of the oil phase under given water cut conditions; S we Effective water saturation, dimensionless; S wr S or These are the unusable water content and oil saturation, respectively, dimensionless; S w n represents water saturation, dimensionless; w n o These are the Corey indices for the aqueous and oil phases, respectively. They are dimensionless, determined experimentally, and used to describe the steepness of the relative permeability-saturation curve.

[0033] Furthermore, in S3, the generalized Carreau model is used to model its viscosity-shear rate relationship, where the viscosity μ of the temporary plugging agent is related to the shear rate. The non-Newtonian shear thinning relation is satisfied, as shown in the following formula:

[0034]

[0035] Where μ0 is the low shear viscosity, μ ∞ λ represents the high shear viscosity limit, and λ, n, and a are the fitting parameters.

[0036] Furthermore, in S4, when a rock mass unit is determined to be a temporarily blocked path, the permeability is controlled by the cumulative effect of the temporary plugging agent concentration. The permeability decreases according to a predetermined function as the concentration increases. The formula for the permeability evolution function in the temporarily blocked zone is as follows:

[0037] k TR =k0·(1-α·C p );

[0038] When the path is determined to be shear-modified, permeability is controlled by plastic shear strain. Permeability increases with increasing shear strain according to a predetermined function. The formula for the permeability evolution function in the TA zone is as follows:

[0039]

[0040] Where: k0 is the initial permeability, C p This refers to the concentration of the temporary plugging agent. This represents the maximum plastic strain.

[0041] Furthermore, S4 adopts a dynamic path discrimination mechanism without a preset fixed modification area. Based on the real-time stress-strain state of the unit and the confining pressure threshold law extracted from the experimental data, it automatically judges its response type and switches the corresponding permeability evolution model parameter set.

[0042] Furthermore, in the S5 coupling iteration process, several sub-loading increments are divided within a single time step. Within each increment, rock mechanical field calculation, permeability and porosity update, and seepage field calculation are performed sequentially. The process is iterated until the changes in the mechanical field and seepage field are both less than the convergence criterion.

[0043] Furthermore, this invention provides a coupled numerical simulation system for the temporary plugging micro-stimulation process of loose sandstone reservoirs, which runs the aforementioned coupled numerical simulation method for the temporary plugging micro-stimulation process of loose sandstone reservoirs, including:

[0044] The nonlinear elastoplastic constitutive modeling module is used to establish a nonlinear elastoplastic constitutive model of reservoir rocks and calculate the mechanical response field during the evolution of geostress under injected perturbation.

[0045] The temporary plugging agent rheological properties module is used to define the shear-thinning non-Newtonian viscosity model of the temporary plugging agent and couple the temporary plugging agent concentration field with the rheological model to dynamically control the permeability of the reservoir seepage field.

[0046] The anisotropic permeability tensor evolution module is used to calculate the anisotropic changes of the permeability tensor based on the rock plastic strain tensor, update the three-dimensional permeability field in real time, and realize the dynamic characterization of the pore channel directionality.

[0047] The response mechanism path branch selection module is used to determine the dominant response mechanism based on the strain state of the unit during the simulation process, and automatically select the corresponding permeability evolution path and parameter set to implement differentiated permeability update strategies for different types of regions.

[0048] The three-field strongly coupled solution module is used to realize the dynamic iterative coupled solution between the pore pressure field, rock stress-strain field and permeability field. It exchanges data synchronously with the external mechanical solver and seepage solver through a preset interface, and completes the update of multiple field variables in each sub-step iteration until convergence.

[0049] The response index extraction and visualization module is used to extract the evolution data of each key variable after the simulation is completed and generate multi-dimensional visualization results. At the same time, it automatically generates a micro-modification plan report based on the simulation output to assist in engineering design decisions.

[0050] Furthermore, the present invention provides an apparatus including a memory, a processor, and an algorithm stored in the memory and executable on the processor, wherein the processor, when executing the computer program, implements the data processing method as described above.

[0051] Furthermore, the present invention provides a computer-readable storage medium storing a computer algorithm, which, when executed by a processor, performs the aforementioned data processing.

[0052] The advantages and positive effects of this invention are:

[0053] 1. Construct a new path discrimination mechanism to improve transformation efficiency.

[0054] This invention proposes a "TR / TA path automaton," in which the unit can automatically select a temporary plugging or shearing evolution path based on strain and concentration conditions during simulation, and can switch between them when conditions are met. This mechanism overcomes the limitations of existing patents that rely on a single flow-dominated system, making the expansion process more consistent with the heterogeneous response of actual reservoirs, thereby significantly improving reservoir connectivity and stimulation efficiency.

[0055] 2. Achieve strong coupling of the three fields of mechanics, seepage, and concentration to reduce the uncertainty of the scheme.

[0056] Most existing technologies use weakly coupled models, failing to simultaneously consider rock mechanical response, plugging agent rheology, and permeability evolution. This invention establishes a three-field strongly coupled framework, introducing plugging agent concentration migration and non-Newtonian rheological properties, enabling simulation results to realistically reflect the dynamic alternation process of plugging control and expansion. This improvement effectively reduces the uncertainty in construction scheme design, ensuring that numerical predictions are closer to actual field conditions.

[0057] 3. Provide an innovative indicator system to optimize design and evaluation.

[0058] This invention not only outputs conventional pressure, strain, and permeability distributions, but also proposes for the first time novel evaluation indicators such as the "μ–κ consistency index," "path residence time distribution," and "shear channel skeleton diagram." These indicators can quantify the plugging range, conduction morphology, and expansion effect, providing a more operational basis for plugging agent placement, injection regime, and segmented design.

[0059] 4. Adaptable to various processes and well types.

[0060] This method can be widely applied to various micro-expansion techniques such as temporary plugging control, shear expansion, thermal recovery connection, and stratified profile modification, and can be flexibly configured according to different conditions of horizontal, vertical, or deviated wells. Compared with existing simulations of single operating conditions, this invention has greater advantages in terms of application scope and applicability.

[0061] 5. Ensures computational stability and is suitable for long-term evaluation.

[0062] By employing adaptive iteration and error control strategies, this method addresses the highly nonlinear problem under strong three-field coupling, ensuring computational stability even under complex boundary conditions and long-term loading. Particularly in the scenario of old thermal recovery areas, this method can continuously track the reservoir's plugging-expansion evolution during long-term injection and production processes, providing reliable support for the entire lifecycle development of oilfields. Attached Figure Description

[0063] Figure 1 This is a schematic diagram of the overall process of an embodiment of the present invention.

[0064] Figure 2 This is a flowchart of the parameter solving process of the coupled numerical simulation method in an embodiment of the present invention.

[0065] Figure 3 This is a schematic diagram of the anisotropic evolution of permeability under shear-induced conditions according to an embodiment of the present invention.

[0066] Figure 4 This is a pressure-volume response curve of the entire process of temporary plugging – shearing – reservoir micro-modification (expansion) using temporary plugging agent according to an embodiment of the present invention.

[0067] Figure 5-6 This is the modeling, simulation, and real-time monitoring and analysis interface for the temporary plugging-reservoir micro-modification process in this embodiment of the invention.

[0068] Figure 7-8 This is a fitting diagram of the second production run of well B11H according to an embodiment of the present invention.

[0069] Figure 9 This is a porosity evolution diagram of the B11H well reservoir after micro-modification, according to an embodiment of the present invention.

[0070] Figure 10This is a permeability evolution diagram of the B11H well reservoir after micro-modification, according to an embodiment of the present invention. Detailed Implementation

[0071] The technical solution of the present invention will now be clearly and completely described with reference to the accompanying drawings. Obviously, the described embodiments are only some, not all, of the embodiments of the present invention. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0072] The embodiments of the present invention will be further described below with reference to the accompanying drawings:

[0073] like Figure 1 As shown, a coupled numerical simulation method for the micro-modification process of temporary plugging in loose sandstone reservoirs is proposed. Based on the strong coupling mechanism of shear-induced permeability tensor evolution and concentration field regulation, a joint solution system of three fields—mechanical field, seepage field, and temporary plugging agent concentration field—is constructed, which includes the following steps.

[0074] S1: Coupled control model construction: Establish a three-field coupled control equation set covering the rock elastoplastic mechanical response, seepage migration behavior and temporary plugging agent concentration migration, forming a control system that can dynamically characterize rock mass deformation, seepage field redistribution and concentration field transport under injection disturbance, driving the response evolution of the entire process of temporary plugging-micro-modification.

[0075] This step aims to establish a dynamic geostress field that can respond to injection disturbances and rock mass deformation throughout the entire process, providing a mechanically driving boundary for subsequent seepage, shear expansion, and other behaviors. The model is based on nonlinear elastic-plastic theory and is divided into two state representations: an elastic stage and a plastic stage, which are respectively applicable to the formation response analysis during the initial disturbance and shear instability stages of micro-remodeling.

[0076] Specifically, S1 includes the following steps.

[0077] S11: Establish the stress-strain relationship between the elastic and plastic stages.

[0078] The elastic stage employs anisotropic nonlinear elastic constitutive relations, as shown in the following formula:

[0079] dσ ij =D ijkl dε kl

[0080] in:

[0081]

[0082] Where: dσ ij D is the increment of the stress tensor, in Pa. ijklIt is a fourth-order tensor, called the tangent modulus tensor or material stiffness tensor, with units of Pa; dε kl σ is the increment of the strain tensor, dimensionless; Pa is the reference pressure, in Pa; n is the exponent, controlling the sensitivity of the modulus to pressure, dimensionless; k is the volume modulus correlation coefficient, dimensionless; g is the shear modulus correlation coefficient, in Pa; P0 is the effective pressure measure, in Pa; σ ij Stress components, in Pa; S mn δ is the deviatoric stress tensor, with units of Pa; ij Let be the Kronecker δ tensor, which is dimensionless.

[0083] This formula fully considers the asymmetry of Poisson's ratio, the variation of bulk modulus, and the influence of shear stress on the modulus. It is a nonlinear generalization of the traditional isotropic elastic model.

[0084] During the plastic stage, when the shear disturbance exceeds the yield condition, the elastic-plastic tangential stiffness model is used to express the following:

[0085] dσ=D ep :dε

[0086] Where: dσ is the stress increment tensor, representing the minute change in stress during loading, with units of Pa; D ep dε is the elastoplastic tangential stiffness matrix, composed of the elastic modulus D, yield function F, plastic potential function G, and softening parameter, which has stronger nonlinear adaptability and is measured in Pa; dε is the strain increment tensor, which represents the minute change in strain of the material during loading and is dimensionless.

[0087] S12: Analyze the definition and advantages of yield criterion and plastic potential function.

[0088] Plastic deformation is treated using the Drucker-Prager improved yield criterion and the non-associated flow rule. The yield function F and plastic potential function G are defined as follows:

[0089] The yield function formula is as follows:

[0090] F=(d0-P t tanβ) 2 +q 2 -p′tanβ-d=0

[0091] in: S=σ+p′I

[0092] The formula for the plastic potential function (non-associated flow) is as follows:

[0093] G=ωd0tanψ+q 2 -p′tanψ=0

[0094] Where: d0 is the equivalent shear strength parameter of the yield surface, obtained by conversion from the rock cohesion-friction angle relationship, and the unit is Pa; P t β is the equivalent compressive yield parameter, used to characterize the compressive strength level of loose sandstone under low confining pressure, in Pa; β is the Drucker-Prager yield surface cone angle, dimensionless; q is the second deviatoric stress invariant, reflecting the shear strength, in Pa; p′ is the effective mean stress, reflecting the confining pressure level, in Pa; d is the yield function constant, used to control the overall position of the yield surface on the p′-q plane; ψ is the expansion angle, the potential surface inclination angle of the non-associated flow method, describing the degree of volume expansion caused by shear deformation, dimensionless; ω is the potential function weighting coefficient, ensuring that the overall dimensions of G are consistent with the yield function, usually determined by material test inversion; I is the unit second tensor.

[0095] This coupling method has the following advantages:

[0096] Strong nonlinear expressive power: The Drucker–Prager criterion is applicable to friction-dominated loose sandstone and can more accurately describe shear-induced yielding behavior.

[0097] Adjustable dilatation angle: ψ controls the volume expansion trend, which helps to fit the volume evolution during the formation of micro-modified shear channels.

[0098] Supports multi-type instability analysis: It can cover the entire process of expansion-rupture-closure, which is different from the single-stage description of classic models such as M-C and Tresca.

[0099] S13: Construct the seepage-stress coupling control equation set.

[0100] To simulate the dynamic response of the reservoir throughout the entire process of temporary plugging and micro-stimulation, a strongly coupled set of governing equations for the in-situ stress field and the seepage field needs to be constructed. Based on the theory of continuum mechanics, the following control system is constructed by jointly using nonlinear elastoplastic mechanical governing equations and non-Newtonian fluid seepage equations:

[0101] The mechanical equilibrium equations are as follows:

[0102]

[0103] in: The divergence operator (taking the divergence of a tensor) calculates and sums the first derivatives of a tensor or vector in spatial coordinates, increasing the dimension of the physical quantity by 1 / m compared to before the operation (e.g., taking the divergence of a stress tensor in Pa yields a quantity in Pa / m; σ′ is the effective stress tensor, in Pa; α is the Biot coefficient, dimensionless, describing the effect of pore pressure on skeletal stress; p is the pore pressure, in Pa; I is the unit tensor, a unit second tensor, dimensionless; b is the body force vector (e.g., gravity), in N / m). 3 Or Pa / m.

[0104] The mass conservation equation for seepage is as follows (considering the effects of consolidation and non-Newtonian fluids):

[0105]

[0106] in: Divergence operator (takes the divergence of a tensor); k a Absolute permeability tensor, m 2 ;k r Relative permeability (functional form as in the Corey model), dimensionless; μ is the shear-dependent viscosity of the non-Newtonian temporary plugging agent solution, in Pa·s; γ w The unit weight of water is N / m³. 3 p represents pore pressure, measured in Pa; ρ represents pore pressure. w This is the density of water, expressed in kg / m³. 3 g is the vector of gravitational acceleration, with units of m / s². 2 m T Equivalent compressibility coefficient, in units of 1 / Pa; ε v Volumetric strain is dimensionless; d is the infinitesimal component of volume change or strain, where d is not a specific physical quantity but a "differential operator"; dt is the infinitesimal element of time, with the unit being seconds.

[0107] Effective penetration rate k eff The expression for effective mobility Λ is as follows:

[0108] k eff =k a ·k r

[0109]

[0110] Where: k a This is the absolute permeability tensor, in meters. 2 Controlled by the principal strain (specifically modeled in S2); k r The relative permeability is dimensionless (specifically introduced in S2 using the Corey empirical model); μ is the effective fluid viscosity, in Pa·s; the aqueous phase is taken according to its actual viscosity μw (generally approximately 1.0 × 10⁻⁶).-3 P a The viscosity of HPAM temporary plugging agent is assigned according to the Carreau model based on the shear rate (·s).

[0111] S2: Shear-induced anisotropic permeability evolution: A shear strain-induced anisotropic permeability tensor evolution model is constructed. The permeability tensor is dynamically changed with tectonic orientation by driving the rock plastic strain tensor. The temporary plugging agent concentration correction is introduced to characterize the effect of plugging agent retention on pore conductivity, thereby directionally enhancing the permeability anisotropy of potential high-permeability channels.

[0112] Specifically, the permeability tensor k in S2 ij The evolution follows an empirical relationship of anisotropic permeability induced by strain:

[0113]

[0114] in, k represents the components of the plastic strain tensor. a Let be the intermediate permeability after concentration control, and b be the sensitivity coefficient of permeability to plastic shear strain. This relationship reflects the superimposed control effect of shear plastic deformation and plugging agent concentration on anisotropic permeability.

[0115] This step aims to realize the absolute permeability tensor k. a Relative penetration rate k r Evolutionary modeling along different principal strain directions reflects the reservoir conduction channel structure reconstruction mechanism induced by shear disturbance. Under the background of strong coupling of seepage, stress, and shear, the model integrates three sub-modules: plastic strain, anisotropic evolution, and shear enhancement coefficient, to construct the permeability tensor response relationship under the shear expansion control mechanism.

[0116] Specifically, S2 includes the following steps.

[0117] S21: Construct the principal strain response model of the anisotropic permeability tensor - absolute permeability.

[0118] Based on experimental observations and the shear expansion mechanism, the response models of the absolute permeability tensor in each principal direction are defined as follows:

[0119]

[0120] in: Initial absolute permeability in each direction (i.e., inherent permeability before stress disturbance), in meters. 2 In the above formula, i = 1, 2, 3 correspond to the three principal directions in the principal stress coordinate system; a is the shear-induced permeability enhancement coefficient (dilatation sensitivity coefficient), which reflects the degree of simultaneous increase in permeability in all directions under the action of strain in the principal direction, and its unit is consistent with the change in permeability, i.e., m.2 b is the cross-coupling enhancement factor, which describes the off-diagonal coupling permeability enhancement effect caused by the rearrangement of asymmetric particle structure and pore connectivity. The unit is also m. 2 The tensor structure reflects the enhanced transverse conduction capability under anisotropic perturbations; εi represents the principal strain in each direction, which is the equivalent axial strain in the corresponding principal direction and is dimensionless. In the above equation, i = 1, 2, 3 are consistent with the principal stress directions.

[0121] The model exhibits a significant directional enhancement characteristic after the plastic strain reaches the shear disturbance threshold, which is helpful for identifying the main seepage control channels during the micro-modification process.

[0122] S22: Evolution of the shear-induced enhancement coefficient λ.

[0123] Coefficients a and b are driven by a uniform shear-induced enhancement factor λ, and their evolution is as follows:

[0124] a=λa0, b=λb0

[0125] λ=0.5(1-h)λ min +0.5(1+h)λ max

[0126]

[0127] Where: a0 and b0 are the initial enhancement coefficients, corresponding to the shear-induced enhancement level before plastic disturbance, and are dimensionless; λ min , λ max The lower and upper limits of shear enhancement control the minimum and maximum values ​​of shear-induced permeability enhancement, respectively, and are dimensionless. Equivalent plastic strain is defined as: The cumulative inelastic deformation of the material during loading is dimensionless; ξ is the shear-induced reinforcement initiation point (e.g., 1%), dimensionless; m is a parameter that adjusts the reinforcement steepness (empirical value 40–60); h is an intermediate control variable, dimensionless.

[0128] This model can dynamically identify the degree of shear perturbation during the expansion process and control the evolution rate and amplitude of the permeability tensor through λ, which is a key path for three-field coupling modeling.

[0129] S23: Construct a relative penetration rate model.

[0130] During the seepage process, the relative permeability k r The assignment method for fluid viscosity μ is as follows:

[0131] Relative penetration rate is expressed using the Corey empirical model:

[0132]

[0133] Where: k rw The relative permeability of the aqueous phase is a dimensionless parameter (0–1), reflecting the effective seepage capacity of the aqueous phase at a given water saturation level; k ro The relative permeability of the oil phase is a dimensionless parameter (0–1), reflecting the effective flow capacity of the oil phase under given water cut conditions; S we Effective water saturation, dimensionless; S wr S or These are the unusable water content and oil saturation, respectively, dimensionless; S w n represents water saturation, dimensionless; w n o These are the Corey indices for the aqueous and oil phases, respectively. They are dimensionless, determined experimentally, and used to describe the steepness of the relative permeability-saturation curve.

[0134] Fluid viscosity assignment method:

[0135] For water: μ = 1;

[0136] For an aqueous solution of HPAM temporary plugging agent, the viscosity is a function of the shear rate, as shown in the following formula:

[0137]

[0138] A power-law model or a Carreau model can be used for fitting; where λ, n, and μ0 are determined by rheological experiments.

[0139] After completing the above model construction, k can be... a k r The effective permeability tensor k was constructed by combining μ and μ. eff , the specific formula is the same as above.

[0140] This tensor is directly used to solve the seepage mass conservation control equation in S13, realizing the simulation of the flow structure response driven by geostress disturbance.

[0141] S3: Non-Newtonian Rheological Modeling of Temporary Plugging Agent: A non-Newtonian shear-thinning rheological model of the temporary plugging agent is established, and the relationship between the shear viscosity of the temporary plugging agent and the shear rate is defined. The effect of concentration field on permeability regulation is introduced into the seepage calculation to simulate the viscosity reduction of the temporary plugging agent in the high shear stress zone and its dynamic plugging control effect on pore seepage.

[0142] Specifically, the temporary plugging agent in S3 is a partially hydrolyzed polyacrylamide (HPAM) system or a functionally equivalent biodegradable polymer colloid, whose viscosity μ is related to the shear rate. Satisfies the non-Newtonian shear thinning relation:

[0143]

[0144] Where μ0 is the low shear viscosity, μ ∞ The high-shear viscosity limit is defined by λ, n, and a, which are fitting parameters. This model is used to characterize the rheological behavior of temporary plugging agents where the viscosity decreases significantly in the high-shear region.

[0145] The construction and module independence of the HPAM temporary plugging agent non-Newtonian rheological model are explained below:

[0146] To accurately simulate the injection-shear-response behavior of HPAM temporary plugging agent in loose sandstone reservoirs, this step focuses on modeling the non-Newtonian rheological properties of the plugging agent fluid, constructing its dynamic viscosity evolution model under different shear conditions, and clarifying its connection with the seepage control equation.

[0147] Specifically, it includes the following aspects:

[0148] Selection and applicability of non-Newtonian viscosity models.

[0149] Considering the shear-thinning characteristics of HPAM temporary plugging agent aqueous solution during reservoir injection, the generalized Carreau model is selected to model its viscosity-shear rate relationship. This model has been used for viscosity assignment in the previous text, and its complete expression is as follows:

[0150]

[0151] in: Shear rate The apparent viscosity is expressed in Pa·s; μ0 is the zero-shear rate viscosity; λ is the shear response time parameter, expressed in s; and n is the shear exponent, dimensionless, reflecting the decreasing viscosity trend.

[0152] If a simplified approach is needed, a power-law model can also be used:

[0153]

[0154] in: Shear rate The apparent viscosity at the specified value is expressed in Pa·s. ρ is the shear rate, in units of 1 / s; n is the shear exponent, dimensionless; K is the consistency index, in units of Pa·s. n , used in power-law models.

[0155] The parameters required for the above model were determined by rheological experiment fitting, and it has good engineering adaptability and computational stability.

[0156] Explanation of the independence of modular design of the model.

[0157] Although this viscosity model has been introduced as a parameter in the seepage equation of S13 above, S3 still presents its modeling process independently for the following reasons:

[0158] The physical boundaries of the modules are clear: S13 only uses viscosity as the input variable, while this step establishes the physical relationship between the shear rate and viscosity functions, providing a basis for solving the governing equations in S5.

[0159] Facilitates unit testing and parameter tuning: The S3 module defines a separate viscosity function, which facilitates sensitivity analysis under different operating conditions (temperature, salinity, etc.).

[0160] It helps to optimize the design of control targets: In the response optimization of S6, the concentration of temporary plugging agent and the injection rate can be adjusted according to the shear response law, thereby improving the plugging stability and micro-modification efficiency.

[0161] The driving effect of the rheology module on the seepage response.

[0162] After completing the viscosity model, its output This will be passed as input to the seepage mass conservation equation, along with the absolute permeability tensor k. a Relative penetration rate k r Jointly construct the effective penetration tensor k eff This coupling mechanism forms a complete chain from shear rate control – non-Newtonian viscosity response – seepage path evolution, which is the key physical basis of the three-field coupling solution logic of this invention.

[0163] S4: Strain-Driven Permeability Evolution Path Identification: A permeability evolution path branch is defined based on the strain state of the rock mass elements. During the numerical simulation, the stress-strain evolution state of each grid element is monitored in real time, and its dominant response mechanism is dynamically identified using the confining pressure response threshold criterion: when the element is dominated by the temporary plugging agent retention effect, it is classified as a temporary plugging path (TR); when it is dominated by the shear deformation expansion effect, it is classified as a shear modification path (TA). For different paths, corresponding permeability evolution model parameter sets are selected, and the element permeability evolution is adjusted in real time to adapt to plugging or expansion requirements.

[0164] Specifically, when a rock mass unit is identified as a temporary plugging path (TR), its permeability is mainly controlled by the cumulative effect of the plugging agent concentration, decreasing according to a predetermined function as the concentration increases. When identified as a shear-modified path (TA), its permeability is mainly controlled by plastic shear strain, increasing according to a predetermined function as the shear strain increases. The permeability evolution functions corresponding to the two paths are as follows. The permeability evolution path in the TR zone is concentration-dominated, satisfying: k TR =k0·(1-α·C p )

[0165] The permeability evolution path in the TA zone is strain-controlled, satisfying the following:

[0166] Where: k0 is the initial permeability, C p This refers to the concentration of the temporary plugging agent. The maximum plastic strain is represented by the TR path, which controls the contraction of pore channels through the concentration field, while the TA path controls the expansion of pore channels through the strain field, thus enabling a quantitative characterization of the dynamic changes in permeability under different mechanisms.

[0167] A dynamic path discrimination mechanism without pre-defined fixed stimulation zones is adopted: during numerical simulation, the range of temporary plugging zones or shear zones in the reservoir is not pre-defined. Instead, the response type is automatically determined based on the real-time stress-strain state of the elements, and the corresponding permeability evolution model parameter set is switched accordingly. This path discrimination algorithm utilizes the segmented characteristics of the stress-strain curves extracted from triaxial cyclic loading tests and the confining pressure threshold law to construct criteria, enabling the model to adaptively identify the locations of temporary plugging failure points and shear channel formation points. Compared with the fixed zone assumption, this dynamic mechanism can more realistically simulate the gradual evolution of plugging paths during reservoir stimulation.

[0168] This step focuses on addressing "how to apply the geostress-seepage-rheology three-field model to specific reservoir structures," achieving spatial partitioning, parameter assignment, boundary setting, and physical state initialization before model solving. This method combines measured data, geological modeling platforms (such as Petrel), modeling tools (such as Geomis), and numerical analysis requirements to complete the entire modeling and initialization process.

[0169] Specifically, it includes the following steps.

[0170] The construction logic of the geostress field. The geostress field is constructed using a method of "measured control + tectonic constraints + local inversion". The core process is as follows: The data source comes from engineering measurements, including results from small-scale fracturing well tests and backflow method geostress tests. The initial values ​​of the three principal stresses are set based on the measured results and the extrapolation function of the depth relationship, and are fitted in layers. The tectonic stress influence zone is jointly corrected by tectonic interpretation and fault slip direction, especially by joint control with the fracture attributes (such as slip direction, fault displacement, and extension range) in the Petrel model. The boundary stress and constraint stress fields are tensor interpolated and boundary projected through the finite element subroutine in Geomis, forming a numerical input field that conforms to the matching law of tectonic distribution and stress field. Region division and mesh control method. The model space discretization adopts a structured mesh and a dual-control meshing method of structure and stress: the reservoir partitioning is consistent with the Petrel model, preserving the original boundary properties of fault planes, interlayer boundaries, and lithological abrupt change surfaces; the mesh is locally refined for injection wells, dominant fault zones, and expected shear zones to meet the accuracy requirements of microscale stress response; the modeling module supports the use of self-developed programs to call Petrel output meshes (.GRDECL or .EGRID files) to directly construct the finite element computation domain; if the Geomis platform is running, the computational FEM mesh (supporting triangular and tetrahedral element partitioning) can be generated with one click, and can be exported in formats such as ABAQUS, COMSOL, and FLAC3D, adapting to multiple numerical platforms.

[0171] The assignment strategy for physical property parameters and initial states. Parameter initialization is performed according to the principle of "regional uniformity + node adjustability": Rock mechanical parameters: taken from triaxial core experiments and in-situ acoustic wave inversion results, including nonlinear elastic modulus, Poisson's ratio, yield parameter, etc.; Fluid rheological parameters: the viscosity-shear rate relationship of HPAM temporary plugging agent is derived from shear rate scanning experiments, and parameters μ0, λ, and n are obtained by indoor fitting, which are described in detail in S3. Initial pore pressure: the fitted value of the formation pressure curve is used, and the fault zone is corrected according to the degree of tectonic stress concentration; Porosity-permeability tensor initialization: based on the shear-induced evolution model in S2, different directional strain fields are input to generate anisotropic initial tensor k. ij All parameters can be imported through Geomis or Petrel attribute tables and automatically mapped to numerical grid cells, ensuring that the input data is consistent with the geological platform.

[0172] Boundary control strategy and dynamic evolution characteristics. Unlike the traditional fixed-boundary-value condition (Dirichlet / Neumann) static assignment method, this invention adopts a boundary tensor response mechanism to support the dynamic boundary response throughout the entire process of formation deformation, stress redistribution, seepage disturbance, and shear coupling: fixed and slip boundaries are automatically classified based on seismic interpretation fault planes; the expansion stress generated by injected disturbances can be projected along the fault boundary to adjacent elements, realizing the nonlinear response of fault plane opening, shear slip, and closure rebound; it supports loading history time series (such as staged injection – steady injection – closure process) to achieve step-by-step loading and coupled feedback. This mechanism avoids the weakness of static models in failing to capture expansion-induced boundary stress changes, effectively enhancing the model's adaptability to complex disturbance conditions.

[0173] S5: Three-Field Coupled Iterative Solution: A strongly coupled numerical model of the pore pressure field, plastic strain field, and permeability tensor field is constructed. A step-by-step loading and sub-step iterative control strategy is employed for the solution: Under each loading increment, the mechanical field is updated sequentially, the permeability field is updated, and the seepage field is solved. A middle interface is used to exchange variables such as pore pressure, permeability tensor, and concentration between the mechanical solver and the seepage solver, repeating the iteration until the solutions for each field converge. This process enables dynamic response prediction for the entire "temporary closure-conduction-expansion" stage.

[0174] Specifically, the three-field coupled calculation is achieved through cross-platform co-simulation. ABAQUS software is used as the stress-strain field solver, and CMG's STARS module is used as the flow-concentration field solver. The two exchange data and perform synchronous iterative control through a pre-developed interface layer. This interface is responsible for transferring the pore compaction strain and plastic zone distribution calculated by ABAQUS to the seepage solution, while simultaneously feeding back the updated pore pressure, fluid saturation, and plugging agent concentration from STARS to ABAQUS, achieving seamless integration of the mechanical-seepage solution in each iteration step.

[0175] The coupled iterative process employs a sub-step loading and nested loop control strategy: within a single time step, several sub-loading increments are defined. Within each increment, rock mechanical field calculations are performed sequentially, the permeability tensor and effective porosity of each element are updated based on the mechanical results, and then the seepage field (including concentration field) is calculated. After completing the above sequence, the convergence of multi-field decoupling is checked. If convergence is not achieved, the number of iterations is increased in the current time step until the changes in both the mechanical field and the seepage field are less than the convergence criterion, before proceeding to the next time step. Through this sub-step iterative loop, strong coupling and stable convergence of multiple fields are achieved, improving computational accuracy.

[0176] This step, based on the aforementioned model construction and parameter initialization, proceeds to the numerical solution stage of the main governing equations. Unlike traditional weakly coupled or explicit sequential solutions, this method employs a nonlinear strongly coupled strategy to jointly solve the three fields of ground stress, seepage, and shear rheology, realistically simulating the multi-physics feedback and structural response throughout the entire process of temporary closure, shearing, and micro-modification.

[0177] Specifically, S5 includes the following steps.

[0178] S51: Solution process and solution path.

[0179] The joint solution of the three field equations adopts a double-nested structure: the outer layer controls the coordinated convergence of the total field quantities (including pressure field, displacement field, and shear state). The inner layer is dominated by Newton-Raphson nonlinear iteration. In each iteration: the principal strain is calculated using the displacement field from the previous step. The permeability tensor and relative permeability are updated in the shear expansion region. The shear rate field is updated and the viscosity is calculated using the rheological model. The updated seepage equation is solved iteratively to obtain the pressure field. The pressure is fed back to the effective stress term to update the plastic region. Volume expansion, seepage velocity, and shear enhancement factor are corrected until convergence. This process can achieve a closed-loop solution of the coupled loop of "stress disturbance – permeability evolution – flow channel adjustment – ​​mechanical feedback".

[0180] S52: Numerical control mechanism and stability assurance.

[0181] To address the numerical oscillations and dissolution issues caused by strong nonlinear coupling, this method introduces multiple control mechanisms:

[0182] Residual control + damping iteration mechanism: Introducing damping coefficients to smooth variable updates in the plastic evolution region improves convergence in the shear expansion stage.

[0183] Variable time-step control algorithm: Adaptively reduce the time step during drastic response phases such as expansion, penetration, and contraction to prevent local non-convergence.

[0184] Local reconstruction of shear-sensitive elements: Local element tensor reconstruction is performed in regions where tensile strain exceeds the threshold to improve the computational stability of boundary slip and local displacement regions.

[0185] Anisotropic tensor symmetry correction: Tensor symmetry checks are automatically performed after each principal stress / permeability update to avoid ill-conditioned elements interfering with convergence.

[0186] S53: Engineering adaptation and platform integration capabilities.

[0187] The coupled solution module boasts excellent openness and engineering compatibility: it can be embedded into the Geomis system to achieve "one-click start of solution – one-click visualization of results." Solution results can be exported to the Petrel / CMG platform for result feedback and fault reinterpretation. Multi-threading supports FLAC3D / ABAQUS steady-state module co-solution, adapting to different owner modeling systems. All parameter fields can generate standard VTK files for visualization of multiple physical quantities such as structural settlement and shear failure zones.

[0188] S6: Results Output and Scheme Optimization: Output multi-dimensional response results, including spatiotemporal evolution data of key variables such as shear rate distribution field, permeability tensor distribution field, pore pressure field, plastic strain field, and temporary plugging agent concentration field. Based on the above simulation results, generate visualization evaluation indicators such as injection pressure difference-injection volume response curve, rupture pressure map of the temporary plugging zone (TR zone), shear modification radius distribution map, and shear channel skeleton structure map. These are used to quantitatively analyze the channel blocking range, shear conduction channel morphology, modification influence radius, and effect area during the temporary plugging micro-modification process, and optimize the temporary plugging agent layout scheme and injection construction parameters accordingly.

[0189] The outputs of the coupled simulation include, but are not limited to, the temporal and spatial evolution data of key physical fields such as the shear rate field, anisotropic permeability tensor field, pore pressure field, plastic strain field, and temporary plugging agent concentration field. Based on the above spatiotemporal data, the following characteristic curves and maps are further generated: injection pressure difference – cumulative injection volume response curve (reflecting the overall pressure response of plugging and diversion), crack initiation pressure distribution map of the temporary plugging zone (TR region) (reflecting the temporary plugging failure point), prediction map of the modification influence radius (reflecting the shear expansion range), and shear channel skeleton structure map (reflecting the network morphology of high-permeability channels), etc. The above output results are used to quantitatively evaluate the plugging control path, response area, modification boundary, and parameter optimization in the temporary plugging micro-modification process, serving as the basis for micro-modification scheme design and on-site construction decisions.

[0190] The method supports multi-condition comparative simulations with multiple parameter combinations as inputs, including but not limited to the following adjustable parameters: initial concentration of the temporary plugging agent and concentration gradient during injection, injection rate and total injection volume, fluid viscosity model parameters, shear yield coefficient (rock shear strength parameter), rock elastic modulus and plastic hardening modulus (affecting rock mass strain response), original confining pressure conditions of the formation, initial distribution of porosity and permeability, etc. By changing the above variable combinations and conducting repeated simulations, the influence of various factors on the effect of temporary plugging-micro-remodeling can be analyzed, and the optimal temporary plugging agent formulation and injection scheme can be selected accordingly. The main controlling factors of the response mechanism under different reservoir conditions can be summarized, and targeted construction parameter sensitivity strategies and optimization guidelines can be formulated.

[0191] After solving the coupled equations, this step focuses on the system output of the simulation results, response field reconstruction, and multidimensional evaluation, providing quantitative basis and engineering optimization suggestions for the entire process of temporary blocking-shearing-micro-modification. This module constructs a unified output system of data extraction, index analysis, and engineering correlation, with multiple functions such as automatic decoupling, clear engineering criteria, and platform-based data transmission.

[0192] Specifically, S6 includes the following steps.

[0193] S61: Response variable extraction and field output.

[0194] After the coupled model calculation is completed, the following core response variables can be output: In-situ stress response field (principal stress tensor, shear stress concentration zone); Shear expansion index (equivalent plastic strain field). Critical shear zone identification). Seepage structure field (equivalent permeability tensor k). eff Streamline diagrams; plugging strength evolution diagrams (fluid pressure drop zone, temporary plugging agent migration path); channel opening behavior identification (shear penetration zone, local pressure drop zone); volume expansion / settlement curves (micro-deformation response caused by reservoir disturbance). All response variables support hourly output and have a three-dimensional visualization structure, which can be used for post-processing tools for profile analysis, profile export, and time series comparison.

[0195] S62: Design of Engineering Objectives and Evaluation Indicators.

[0196] To scientifically evaluate the response effect and technical value of the method of this invention in the temporary plugging-micro-modification process, a number of quantifiable engineering evaluation indicators are set, covering aspects such as plugging strength, shear induction effect, seepage structure evolution and parameter optimization, and a complete evaluation framework is constructed.

[0197] First, regarding the sealing performance, the blocking pressure difference index (denoted as Δp) is used. block This indicator reflects the ability of HPAM temporary plugging agent to form a high retardation pressure zone within the pore structure. This indicator can be directly used to compare the plugging efficiency under different injection concentrations or shear conditions.

[0198] Secondly, regarding the shear response characteristics, a shear enhancement coefficient difference λ is introduced. max -λ min As a quantitative measure, it is used to characterize the magnitude and distribution of permeability enhancement within the main strain-triggered region, and this index can identify the shear sensitivity of the micro-modified region.

[0199] To reflect the degree of effective reconstruction of reservoir spatial structure under shear-induced conditions, two key volume indices are defined: the expansion trigger zone (TR zone, denoted as V). TR ) and the shear conduction region (TA region, volume denoted as V) TAAmong them, the TR region represents the active region of plastic expansion, and the TA region represents the permeable main control channel formed after shearing. The two constitute the main physical path of micro-modification.

[0200] Meanwhile, to determine the critical moment of reservoir fracture or shear channel opening, a breakthrough critical time t is introduced. break As an important dividing point in the evolutionary process, this value is extracted based on the abrupt change of equivalent variability and local pressure drop, and can be used in engineering to identify sensitive nodes of injection operations.

[0201] In terms of evaluating seepage structures, a structural bending index is proposed. Where L real L represents the actual path length of the fluid in the simulation. geo This represents the geometric shortest path. This index reflects the degree of flow around the structure and structural disturbance of the temporary plugging agent fluid during the plugging-opening-breaking process, and is an important indicator for judging the heterogeneous evolution of seepage.

[0202] Furthermore, the model supports parameter sensitivity inversion and optimal solution recommendation. The injection concentration C can be adjusted based on the simulation process. p Injection rate q inj With injection duration T inj The response trend is used to output a set of optimal parameter combinations (C). p ,q inj ,T inj ) opt It is used to guide subsequent on-site construction design and parameter selection, significantly improving the operability and economy of the project.

[0203] This invention also provides a coupled numerical simulation system for the temporary plugging micro-stimulation process of loose sandstone reservoirs, which runs the aforementioned coupled numerical simulation method for the temporary plugging micro-stimulation process of loose sandstone reservoirs, including:

[0204] The nonlinear elastoplastic constitutive modeling module is used to establish a nonlinear elastoplastic constitutive model of reservoir rocks and calculate the mechanical response field during the stress evolution process under injected disturbances.

[0205] The temporary plugging agent rheological properties module is used to define the shear-thinning non-Newtonian viscosity model of the temporary plugging agent and couple the temporary plugging agent concentration field with the rheological model to dynamically control the permeability of the reservoir seepage field.

[0206] The anisotropic permeability tensor evolution module is used to calculate the anisotropic changes of the permeability tensor based on the rock plastic strain tensor, update the three-dimensional permeability field in real time, and realize the dynamic characterization of the pore channel directionality.

[0207] The response mechanism path branch selection module is used to determine the dominant response mechanism based on the strain state of the unit during the simulation process, and automatically select the corresponding permeability evolution path and parameter set to implement differentiated permeability update strategies for different types of regions.

[0208] The three-field strongly coupled solution module is used to realize the dynamic iterative coupled solution between the pore pressure field, rock stress-strain field and permeability field. It exchanges data synchronously with the external mechanical solver and seepage solver through a preset interface, and completes the update of multiple field variables in each sub-step iteration until convergence.

[0209] The three-field strongly coupled solution module includes a built-in CMG–ABAQUS linkage interface program, supporting coupled calculations with CMG reservoir numerical simulation software and ABAQUS finite element software. This interface module can automatically complete the data mapping between the seepage grid and the mechanical grid, and connect the solution step size of the two software programs to achieve synchronous updates and iterative control of permeability field, pressure field, and stress-strain field variables, thereby ensuring the stability and accuracy of multi-field coupled calculations.

[0210] The response index extraction and visualization module is used to extract the evolution data of key variables after the simulation is completed and generate multi-dimensional visualization results, including injection pressure difference-volume curves, shear rate cloud maps, permeability structure profiles, and impact range maps of temporary plugging modification. At the same time, it automatically generates micro-modification scheme reports based on simulation outputs, including suggestions for the layout of temporary plugging areas, suggestions for optimizing injection construction systems, and modification radius control strategies, to assist in engineering design decisions.

[0211] Among them, the response index extraction and visualization module has intelligent analysis and report generation functions. It can automatically identify the shear hyperpermeability channel structure formed based on the evolution trajectory of the permeability tensor, extract the optimal injection parameter combination of the temporary plugging agent by combining the pressure difference-injection volume response, and automatically generate a micro-modification scheme design report that includes the optimized layout of the temporary plugging zone, injection cycle and discharge suggestions, expected modification radius and production increase effect, providing quantitative basis for on-site implementation.

[0212] The coupled simulation mechanism described herein is applicable to offshore loose sandstone reservoirs with strong heterogeneity and sensitivity, and can be used in various production enhancement and stimulation processes such as temporary plugging and displacement, shear expansion, thermal recovery and conduction, and selective profile modification. Through the simulation of this invention, it is possible to identify and locate high-permeability channels within the reservoir, optimize the design of combined temporary plugging and shear stimulation paths, and quantitatively assist in the formulation of construction parameters (such as injection pressure and injection volume), thereby meeting the requirements of on-site engineering for the predictive accuracy and design efficiency of micro-stimulation schemes.

[0213] The experimental data used for model establishment and parameter initialization include, but are not limited to: sequence stratigraphy of the geological model, spatial heterogeneous distribution of reservoir porosity and permeability (obtained from core analysis and well logging interpretation); mechanical-flow coupling data such as rock elastic modulus, yield criterion parameters, permeability-effective stress response curves, and shear expansion thresholds obtained from triaxial and true triaxial water injection experiments; and chemical experimental data related to temporary plugging agents / expansion materials, such as viscosity-shear rate curves and concentration-permeability influence relationships. The method utilizes the aforementioned initial parameters for model calibration, and during the simulation process, it continuously inverts key uncertain parameters through multiple rounds of "simulation-verification-parameter tuning" cycles, ultimately achieving a closed-loop optimization design of the reservoir micro-stimulation scheme from parameter selection to effect prediction.

[0214] The following uses the second round of reservoir micro-stimulation in well B11H of a loose sandstone oil reservoir in the Bohai Sea as an example to illustrate the invention in detail:

[0215] This well area is characterized by loose pore-throat structure, significant elastoplastic response, complex polymer injection history, and strong reservoir heterogeneity. To address issues such as fixed dominant channels and poor response of inefficient layers, the project team established a complete coupled numerical simulation system of temporary plugging-shearing-micro-modification and achieved seamless integration with the Petrel-CMG integrated modeling platform.

[0216] S1: As Figure 2 As shown, the entire modeling process begins with the parameter solution logic. First, rock mechanics and seepage parameters such as elastic-plastic modulus, Biot coefficient, initial porosity, and permeability are obtained through triaxial fatigue tests and true triaxial water injection tests. Then, combined with backflow-assisted in-situ stress test data, an elastic-plastic-yield-seepage coupled control equation set is constructed. Test results show that this region exhibits a significant nonlinear permeability tensor evolution trend under shear disturbance conditions, requiring close consideration of anisotropic enhancement paths.

[0217] In the core modeling module, the shear-induced permeability evolution model proposed in this invention is adopted, such as... Figure 3 As shown, the direction of the principal axis of the permeability tensor is adjusted by driving the plastic strain tensor, thereby achieving directional linkage between the seepage structure and rock mass deformation. This mechanism effectively reflects the formation of new conductive regions along the "weak expansion axis" by the temporary plugging agent under shear disturbance, and dynamically updates the seepage field distribution. Unlike traditional isotropic expansion simulations, this model can gradually adjust the direction of the seepage channels during the pressurization stage, exhibiting a stronger capacity for expansion in low-permeability zones.

[0218] Response simulation of the entire process of temporary plugging agent injection – temporary plugging – shear expansion, as follows: Figure 4As shown, the simulated pressure-volume relationship curve closely matches the field monitoring data, successfully capturing the stress reconstruction behavior in the three stages of plugging, unblocking, and rupture. The shear-induced λ enhancement coefficient in the model varies nonlinearly with the principal strain, significantly improving the directional selectivity of the channel after plugging, preventing the continued conduction of high-permeability channels, and increasing the overall effective expansion ratio.

[0219] In terms of simulation and monitoring platform construction, a real-time feedback-parameter tuning system was built, such as... Figure 5 , Figure 6 As shown, the model output data is compared with the real-time reservoir monitoring curves to dynamically update shear zone parameters, plugging agent concentration field, and anisotropic permeability tensor evolution path, achieving closed-loop control of "production, modeling, and optimization simultaneously." This platform integrates the Petrel static model, the CMG-STARS expansion module, and self-compiled coupled control equations, realizing unified dual-scale prediction of reservoir macro-response and micro-evolution.

[0220] like Figure 7-10 The figure shows the production fitting curve and porosity-permeability evolution results of the second round of reservoir micro-reform in well B11H. After considering non-Newtonian fluid transport and shear-induced tensor evolution, the model achieves a production fitting error better than ±5%, with a 10–15% increase in porosity in the expanded region and a maximum permeability increase of over 3 times along the principal shear axis, demonstrating the model's accurate inversion capability for engineering control responses. During the fitting process, the model uses the actual injected concentration field and incorporates iterative correction based on the inter-well pressure response in the micro-reformed area to ensure consistency between the numerical simulation and the field response.

[0221] In summary, this invention successfully achieved reservoir shear expansion directionality control and quantitative assessment of micro-stirring in this well area, proving that the model not only possesses the integrity of coupled control but also has field-verifiable engineering applicability, providing strong technical support for subsequent reservoir-like implementation.

[0222] The advantages and positive effects of this invention are:

[0223] 1. Construct a new path discrimination mechanism to improve transformation efficiency.

[0224] This invention proposes a "TR / TA path automaton," in which the unit can automatically select a temporary plugging or shearing evolution path based on strain and concentration conditions during simulation, and can switch between them when conditions are met. This mechanism overcomes the limitations of existing patents that rely on a single flow-dominated system, making the expansion process more consistent with the heterogeneous response of actual reservoirs, thereby significantly improving reservoir connectivity and stimulation efficiency.

[0225] 2. Achieve strong coupling of the three fields of mechanics, seepage, and concentration to reduce the uncertainty of the scheme.

[0226] Most existing technologies use weakly coupled models, failing to simultaneously consider rock mechanical response, plugging agent rheology, and permeability evolution. This invention establishes a three-field strongly coupled framework, introducing plugging agent concentration migration and non-Newtonian rheological properties, enabling simulation results to realistically reflect the dynamic alternation process of plugging control and expansion. This improvement effectively reduces the uncertainty in construction scheme design, ensuring that numerical predictions are closer to actual field conditions.

[0227] 3. Provide an innovative indicator system to optimize design and evaluation.

[0228] This invention not only outputs conventional pressure, strain, and permeability distributions, but also proposes for the first time novel evaluation indicators such as the "μ–κ consistency index," "path residence time distribution," and "shear channel skeleton diagram." These indicators can quantify the plugging range, conduction morphology, and expansion effect, providing a more operational basis for plugging agent placement, injection regime, and segmented design.

[0229] 4. Adaptable to various processes and well types.

[0230] This method can be widely applied to various micro-expansion techniques such as temporary plugging control, shear expansion, thermal recovery connection, and stratified profile modification, and can be flexibly configured according to different conditions of horizontal, vertical, or deviated wells. Compared with existing simulations of single operating conditions, this invention has greater advantages in terms of application scope and applicability.

[0231] 5. Ensures computational stability and is suitable for long-term evaluation.

[0232] By employing adaptive iteration and error control strategies, this method addresses the highly nonlinear problem under strong three-field coupling, ensuring computational stability even under complex boundary conditions and long-term loading. Particularly in the scenario of old thermal recovery areas, this method can continuously track the reservoir's plugging-expansion evolution during long-term injection and production processes, providing reliable support for the entire lifecycle development of oilfields.

[0233] The foregoing has provided a detailed description of one embodiment of the present invention, but this description is merely a preferred embodiment and should not be construed as limiting the scope of the invention. All equivalent variations and modifications made within the scope of the claims of this invention should still fall within the patent coverage of this invention.

Claims

1. A coupling numerical simulation method of a loose sandstone reservoir temporary plugging micro-reformation process, characterized in that: The method comprises the following steps: S1: constructing a control equation set of three-field coupling of a mechanical field, a seepage field and a temporary plugging agent concentration field, to dynamically represent rock mass deformation, seepage field redistribution and concentration field transmission under injection disturbance, and drive response evolution of the whole process of temporary plugging-micro-reformation; S2: constructing an anisotropic permeability tensor evolution model induced by shear strain, driving the permeability tensor to dynamically change through a rock plastic strain tensor, and introducing temporary plugging agent concentration to correct the permeability; S3: establishing a non-Newtonian shear thinning rheological model of the temporary plugging agent, defining the relationship between shear viscosity of the temporary plugging agent and shear rate, and coupling the concentration field to the regulation effect of the permeability in seepage calculation; S4: setting a permeability evolution path branch driven by a rock mass unit strain state, in the simulation process, according to the stress-strain evolution state of the unit, dynamically identifying the dominant response mechanism thereof: when the unit is mainly affected by the temporary plugging agent retention effect, it is classified as a temporary plugging type path; when the unit is mainly affected by the shear deformation dilatancy effect, it is classified as a shear reformation type path; corresponding permeability evolution model parameter groups are selected for different paths, and the unit permeability evolution is adjusted in real time; S5: constructing a three-field strong coupling numerical model of a pore pressure field, a plastic strain field and a permeability tensor field, and realizing solving by using a step-by-step loading and sub-step iteration control strategy: under each loading increment, the mechanical field is updated, the permeability field is updated and the seepage field is solved in sequence, and the pore pressure, the permeability tensor and the concentration variable are exchanged between the mechanical solver and the seepage solver through an intermediate interface, and iteration is repeated until the solutions of the fields converge; S6: outputting time-space evolution data of multi-dimensional key variables, generating injection pressure difference-injection volume response curves, temporary plugging zone fracture pressure maps, shear reformation radius distribution maps and shear channel skeleton structure maps, visual evaluation indexes, optimizing temporary plugging agent arrangement schemes and injection construction parameters.

2. The method according to claim 1, wherein the method is characterized by: The S1 comprises the following steps: S11: establishing a stress-strain relationship in an elastic stage and a plastic stage; S12: using a Drucker-Prager improved yield criterion and a non-associated flow rule for plastic deformation, defining a yield function F and a plastic potential function G, and the formulas are as follows: F = (d0 - P t tan β) 2 + q 2 - p'tan β - d = 0 where: S = σ + p' I G = ωd0 tan ψ + q 2 - p' tan ψ = 0 where d0 is the equivalent shear strength parameter of the yield surface, which is converted from the cohesion-friction angle relationship of rock, with the unit of Pa; P t is the equivalent compressive yield parameter, which is used to represent the compression strength level of loose sandstone under low confining pressure, with the unit of Pa; β is the Drucker-Prager yield surface cone angle, dimensionless; q is the second deviatoric stress invariant, which reflects the shear action strength, with the unit of Pa; p' is the effective mean stress, which reflects the confining pressure level, with the unit of Pa; d is the yield function constant, which is used to control the overall position of the yield surface on the p'-q plane; ψ is the dilatancy angle, which is the potential surface inclination of the non-associated flow method, and describes the degree of volume expansion caused by shear deformation, dimensionless; ω is the potential function weight coefficient, which ensures that the overall dimension of G is consistent with the yield function, and is usually determined according to material test inversion; I is the unit second order tensor; S13: constructing a seepage-stress coupling control equation set.

3. The coupled numerical simulation method of the temporary plugging and micro-reformation process of the unconsolidated sandstone reservoir according to claim 1 or 2, characterized in that: The S2 comprises the following steps: S21: constructing a principal strain response model of an anisotropic permeability tensor, and the response model of each principal direction of the absolute permeability tensor is defined as follows: wherein: Absolute permeability in each direction at the initial state, i.e. the intrinsic permeability before stress disturbance, unit: m 2 ; in the above formula, i = 1, 2, 3 respectively correspond to the three principal directions under the principal stress coordinate system; a is a shear-induced permeability enhancement coefficient, i.e. a shear dilation sensitivity coefficient, reflecting the simultaneous promotion degree of the permeability in each direction under the action of the principal direction strain, the unit is consistent with the permeability change, i.e. m 2 ; b is a cross-coupling enhancement coefficient, describing the non-diagonal coupling permeability enhancement effect caused by the rearrangement of the asymmetric particle structure and pore connectivity; the unit is m 2 ; the tensor structure reflects the response enhancement of the cross-coupling under the non-isotropic disturbance; ε i Principal strain in each direction, which is the equivalent axial strain in the corresponding principal direction, dimensionless; in the above formula, i = 1, 2, 3 are consistent with the principal stress direction; S22: evolving a shear-induced enhancement coefficient λ, and the evolution form is as follows: a = λa0, b = λb0 λ = 0.5(1 - h)λ min + 0.5(1 + h)λ max where: a0, b0are initial enhancement coefficients, dimensionless, corresponding to the shear-induced enhancement level before plastic disturbance occurs; λ min , λ max Shear enhancement lower and upper limit values, respectively control the minimum and maximum values of shear-induced permeability enhancement, dimensionless; is the equivalent plastic strain, defined as: characterizes the accumulated inelastic deformation of the material during loading, dimensionless; ξ is the shear-induced enhancement starting point (e.g., 1%), dimensionless; m is a parameter that adjusts the steepness of the enhancement (empirically valued between 40-60); h is an intermediate control variable, dimensionless; S23: constructing a relative permeability model, and the relative permeability model is expressed by using a Corey empirical model: where: k rw Relative permeability of the water phase, dimensionless parameter (0-1) reflecting the effective flow capacity of the water phase at a given water saturation; k ro Relative permeability of the oil phase, dimensionless parameter (0-1) reflecting the effective flow capacity of the oil phase at a given water saturation; S we Effective water saturation, dimensionless; S wr , S or are the irreducible water and oil saturations, respectively, dimensionless; S w is the water saturation, dimensionless; n w , n o are the Corey exponents for the water and oil phases, respectively, dimensionless, determined experimentally, and used to describe the steepness of the relative permeability-saturation curve.

4. The method according to claim 1 or 2, characterized in that: In the S3, the generalized Carreau model is selected to model the viscosity-shear rate relationship, the viscosity μ of the temporary plugging agent and the shear rate satisfies the non-Newtonian shear thinning relationship, the formula is as follows: where μ0is the low shear viscosity, μ ∞ is the high shear viscosity limit, and λ, n, a are fitting parameters.

5. The method according to claim 1 or 2, wherein the method is characterized by: In the S4, when the rock mass unit is determined to be a temporary plugging type path, the permeability is controlled by the concentration accumulation effect of the temporary plugging agent, the permeability decreases with the increase of the concentration according to a predetermined function, and the permeability evolution function formula of the temporary plugging type zone is as follows, k TR = k0·(1 - a·C p ); When it is determined to be a shear reformation type path, the permeability is controlled by the plastic shear strain, the permeability increases with the increase of the shear strain according to a predetermined function, and the permeability evolution function formula of the TA zone is as follows: where: k0 is the initial permeability, C p is the temporary plugging agent concentration, is the maximum plastic strain.

6. The method according to claim 1 or 2, wherein the method is characterized by: The S4 adopts a dynamic path discrimination mechanism without preset fixed reconstruction area, automatically discriminates the response type based on the real-time stress-strain state of the unit and the confining pressure threshold law extracted from experimental data, and switches the corresponding permeability evolution model parameter set.

7. The method according to claim 1 or 2, wherein the method is characterized by: In the S5, the three-field coupling calculation is simulated by ABAQUS software as a stress-strain field solver and CMG STARS module as a flow-concentration field solver, data exchange and synchronous iteration control are realized through a pre-developed interface middle layer, and in the coupling iteration process, a plurality of sub-loading increments are divided in a single time step, and the rock mechanics field calculation, permeability and porosity update, and seepage field calculation are sequentially performed in each increment, and the iteration is repeated until the changes of the mechanics field and the seepage field are less than the convergence criterion.

8. A coupled numerical simulation system of a temporary plugging micro-reformation process for a loose sandstone reservoir, characterized in that: The loose sandstone reservoir temporary plugging micro-reformation process coupling numerical simulation method of any one of claims 1-7 is run, comprising, a nonlinear elastoplasticity constitutive modeling module for establishing a nonlinear elastoplasticity constitutive model of the reservoir rock, and calculating the mechanical response field in the evolution process of the ground stress under the injected disturbance; a temporary plugging agent rheological property module for defining a shear thinning non-Newtonian viscosity model of the temporary plugging agent, and coupling the temporary plugging agent concentration field with the rheological model to dynamically regulate the reservoir seepage field permeability; an anisotropic permeability tensor evolution module for calculating the anisotropic change of the permeability tensor according to the rock plastic strain tensor, and real-time updating the three-dimensional permeability field to realize dynamic characterization of the pore channel directionality; a response mechanism path branch selection module for determining the dominant response mechanism of the unit based on the unit strain state during the simulation, and automatically selecting the corresponding permeability evolution path and parameter set to implement differentiated permeability updating strategies for different types of regions; a three-field strong coupling solving module for realizing dynamic iterative coupling solving among the pore pressure field, the rock stress-strain field, and the permeability field, and realizing data synchronous exchange with the external mechanics solver and the seepage solver through a preset interface, and completing multi-field variable updating in each sub-step iteration until convergence; a response index extraction and visualization module for extracting the evolution data of each key variable after the simulation is completed, and generating multi-dimensional visualization results, and automatically generating a micro-reformation scheme report according to the simulation output to assist engineering design decisions.

9. An apparatus comprising a memory, a processor, and an algorithm stored in the memory and executable on the processor, wherein: The processor executes the computer program to realize the data processing method of any one of claims 1-7.

10. A computer readable storage medium having stored thereon a computer algorithm, the computer algorithm comprising: The computer algorithm is executed by the processor to realize the data processing of any one of claims 1-7.

Images