A numerical simulation method for fracturing horizontal wells in shale oil reservoirs with high clay content

Abstract

The present invention provides a numerical simulation method for fractured horizontal wells in shale oil reservoirs with high clay content. The method includes: Step S1, obtaining reservoir-related parameters according to the actual characteristics of the shale oil reservoir and core, logging, and fluid test data; Step S2, establishing a prediction model for reservoir development dynamic parameters of fractured wells in shale oil reservoirs based on the discrete fracture model; Step S3, performing a prediction simulation of reservoir development dynamic parameters during the fracturing process; Step S4, performing a prediction simulation of reservoir development dynamic parameters during the shut-in process; Step S5, performing a prediction simulation of reservoir development dynamic parameters during the fracturing flowback process. The present invention comprehensively considers the coupling effects of osmotic pressure and heat flow solidification to predict the fracturing and development dynamics of shale oil reservoirs with high clay content, determine its main influencing factors, and thus provide guidance for the optimization of shale oil reservoir fracturing.

Description

Technical Field

[0001] The present invention belongs to the field of oil and gas field development, and specifically relates to a numerical simulation method for fracturing horizontal wells in shale oil reservoirs with high clay content. Background Art

[0002] For shale oil reservoirs with relatively high clay minerals, the salinity of formation water can be as high as hundreds of thousands of ppm. During the fracturing process, the osmotic pressure effect formed with the injected low-salinity fracturing fluid is significant, which has a great impact on fracturing and flowback. During the fracturing process of shale oil reservoirs, there is a coupled process of temperature (T), seepage (H), stress (M), and chemistry (C). The fluid in the reservoir flows in porous media and fractures, which will affect the distribution of the temperature field. The change of the temperature field causes the change of fluid physical properties, which reacts on the seepage field. It will also cause the change of osmotic pressure; the fluid movement will cause the convective action of salt ions. Due to the shale semi-permeable membrane effect, the convective diffusion of salt ions will affect the magnitude of osmotic pressure, resulting in the change of fluid movement speed. Due to the constraint of in-situ stress, the volume expands and contracts, causing thermal stress, resulting in the change of the stress field. The change of the stress field causes the change of the energy of reservoir rock mass. According to the principle of energy conservation, it will lead to the change of the temperature field. During the injection or production process, the change of pore pressure will cause the deformation of the oil reservoir rock skeleton, resulting in the increase or decrease of porosity and permeability, thus affecting the seepage field. Therefore, it is necessary to comprehensively consider the osmotic pressure and the coupling effect of heat flow solidification, predict the fracturing and development dynamics of shale oil reservoirs with high clay content, determine its main influencing factors, and then provide guidance for the optimization of shale oil reservoir fracturing. Summary of the Invention

[0003] The present invention comprehensively considers the osmotic pressure and the coupling effect of heat flow solidification, predicts the fracturing and development dynamics of shale oil reservoirs with high clay content, determines its main influencing factors, and then provides guidance for the optimization of shale oil reservoir fracturing.

[0004] In a first aspect, an embodiment of the present application provides a numerical simulation method for fracturing horizontal wells in shale oil reservoirs with high clay content, including:

[0005] Step S1: Obtain reservoir physical property parameters according to the actual characteristics of the shale oil reservoir and core, logging, and fluid test data; obtain the fracture geometric information of the shale oil reservoir, take the fractures in the reservoir as the internal boundaries of the reservoir for dimension reduction processing, establish a reservoir geometric model, and then use triangular grids to perform geometric meshing on the reservoir geometric model to form discrete elements;

[0006] Step S2: Establish a prediction model for reservoir development dynamic parameters of a fractured well in a shale oil reservoir based on a discrete fracture model;

[0007] Step S3: Perform a prediction simulation of reservoir development dynamic parameters during the fracturing process;

[0008] Step S4, perform predictive simulation of reservoir development dynamic parameters during the well shut-in process;

[0009] Step S5, perform predictive simulation of reservoir development dynamic parameters during the fracturing flowback process.

[0010] Among them, Step S1 includes:

[0011] Collect reservoir geological parameters, including the maximum and minimum horizontal principal stress distributions, Young's modulus of the rock, Poisson's ratio, formation temperature and pressure, porosity, permeability, semipermeable membrane efficiency; collect reservoir fluid parameters, including viscosity, salinity, thermal conductivity, specific heat capacity, relative permeability curve, capillary force; collect completion information of the fractured well, including cluster spacing; collect fracturing construction parameters, including construction displacement, fracturing fluid viscosity, salinity, thermal conductivity, specific heat capacity; collect the distributions of hydraulic fractures and natural fractures, and the fracture length, width, porosity and permeability; obtain shale oil reservoir fracture geometric information based on the actual geological data of the reservoir or existing geological model data: including fracture location, size, density, orientation, aperture.

[0012] Among them, Step S2 includes:

[0013] Step S21: Establish two-phase oil-water control equations for the matrix and fractures respectively:

[0014] The matrix water phase control equation is:

[0015]

[0016] The matrix oil phase equation is:

[0017]

[0018] The fracture water phase equation is:

[0019]

[0020] The fracture oil phase equation is:

[0021]

[0022] Step S22: Establish a salt ion migration control equation;

[0023] Step S23: Establish a stress field control equation. The stress field control is mainly controlled by the equilibrium equation, constitutive equation and displacement equation. The constitutive equation of an isotropic linear thermoelastic material is:

[0024] σ′ ij =2Gε ij +λε kk δ ij -K′α T Tδ ij -αB pδ ij

[0025] In the formula, σ ij is the stress tensor, representing the stress state at a point, in Pa; ε ij is the strain tensor, dimensionless; u i is the displacement tensor, in m; λ is the Lame constant, in Pa,; G is the shear modulus, in Pa,; K′ is the bulk modulus of elasticity, in Pa,; α B is the Biot coefficient, dimensionless;; F is the body force, in N / m, p′ is the mean pore pressure, in Pa;

[0026] p′ = S w p w + S o p o = S w (p o - p c ) + S o p o = p o - S w p c

[0027] The effective stress calculation formula proposed by Terzaghi and modified by Skempton is:

[0028] σ′ = σ - α B p

[0029] The quasi-static equilibrium differential equation is:

[0030] σ′ ij,j + F i = 0

[0031] The relationship between strain and displacement is:

[0032]

[0033] In the formula, ε ij is the strain tensor; u i is the displacement tensor, in m; F is the body force, in Pa;

[0034] The stress field equation uses the solid mechanics interface, and corresponding parameters can be input;

[0035] Step S24: Establish the temperature field control equation. According to the law of conservation of energy, the temperature field control equation is:

[0036]

[0037]

[0038] In the formula, (ρCp ) eff is the effective specific heat capacity of the rock mass, λ eff is the effective thermal conductivity of the rock mass, η eff is the effective heat convection coefficient of the fluid, C s , C o , C w are the specific heat capacities of rock skeleton, oil and water, respectively, s ,λ o ,λ w are the thermal conductivity of rock skeleton, oil and water respectively;

[0039] Step S25: Use COMSOL Multiphysic software to solve.

[0040] Wherein, step S25 comprises:

[0041] COMSOL Multiphysic software is used. This software is based on the finite element method and provides a fully coupled solution method for multiple physical fields. In the coupled analysis, the four fields of heat flow and solidification can be combined to form a unified set of coupled equations for solution, and the dependent variables of each independent field can be calculated at the same time. The PDE module is used to solve the oil-water two-phase flow problem, the porous medium heat transfer module is used to solve the temperature field problem, and the solid mechanics is used to solve the stress field problem. COMSOL Multiphysic combines the three to form a set of differential equations expressed in a general formula, realizing the fully coupled solution of the four fields of THMC.

[0042] Wherein, step S3 comprises:

[0043] Step S31: performing a prediction simulation of reservoir development dynamic parameters during the fracturing process according to the mathematical model and boundary conditions;

[0044] Step S32: Output numerical simulation results, including reservoir temperature, pressure, water saturation, and salt concentration.

[0045] Wherein, step S31 includes:

[0046] Assume that the initial time is 0, the injection time is t1, the shut-in time is t2, and the flowback time is t3;

[0047] (1) Initial conditions

[0048] The initial conditions of temperature and pressure in the fracturing process are the initial formation temperature and pressure of the reservoir, which are realized by the Dirichlet boundary condition, that is,

[0049] T=T0(t=0)

[0050] p=p0(t=0)

[0051] where: T0 is the initial formation temperature, K; p0 is the initial formation pressure, Pa;

[0052] For the stress field, since the reservoir is constrained by in-situ stresses, it is in a state of having stress but no strain in the initial state, which is achieved by applying boundary loads, i.e.,

[0053] σ ij (x,y,z,t = 0) = [σ v σ H σ h

[0054] In the formula, σ v is the vertical in-situ stress, Pa; σ H is the maximum horizontal in-situ stress, Pa; σ h is the minimum horizontal in-situ stress, Pa;

[0055] (2) Boundary conditions

[0056] During the fracturing process, the inner boundary condition is the flow rate boundary condition, which is mentioned in the weak form, i.e.,

[0057]

[0058] In the formula, A is the cross-sectional area, m 2 ; q w is the flow rate of the injected fracturing fluid, kg·m 3 / s, p w is the bottom-hole flowing pressure;

[0059] The boundary conditions of the temperature field and the seepage field are similar,

[0060]

[0061] In the formula, q is the heat source term caused by the injection well, W / m 2 ; C pw is the specific heat capacity of the injected water; T inj is the injection temperature, K.

[0062] Among them, step S4 includes:

[0063] Step S41: Predict and simulate the reservoir development dynamic parameters during the shut-in process according to the mathematical model and boundary conditions;

[0064] Step S42: Output the numerical simulation results, including reservoir temperature, pressure, water saturation, and salt concentration.

[0065] Among them, step S41 includes:

[0066] (1) Initial conditions

[0067] ​The initial conditions of temperature and pressure during the well shut-in process are the reservoir temperature and pressure at the end of fracturing, which are achieved through Dirichlet boundary conditions, i.e.,

[0068] T = T1(t = t1)

[0069] p = p1(t = t1)

[0070] where: T1 is the reservoir temperature at the end of fracturing, K; p1 is the reservoir pressure at the end of fracturing, Pa;

[0071] For the stress field, since the reservoir is constrained by in-situ stresses, the initial state is in a state of having stress but no strain, which is achieved by applying boundary loads, i.e.,

[0072] σ ij (x,y,z,t = 0) = [σ v σ H σ h

[0073] where σ v is the vertical in-situ stress, Pa; σ H is the maximum horizontal in-situ stress, Pa; σ h is the minimum horizontal in-situ stress, Pa.

[0074] Among them, step S5 includes:

[0075] Step S51: Predict and simulate the reservoir development dynamic parameters during the fracturing flowback process according to the mathematical model and boundary conditions;

[0076] Step S52: Output the numerical simulation results, including reservoir temperature, pressure, water saturation, and salt concentration.

[0077] Among them, step S51 includes:

[0078] Assume that the initial moment is the 0 moment, inject for t1 time, shut in for t2 time, and flow back for t3 time;

[0079] (1) Initial conditions

[0080] The initial conditions of temperature and pressure during the fracturing flowback process are the reservoir temperature and pressure at the end of well shut-in, which are achieved through Dirichlet boundary conditions, i.e.,

[0081] T = T2(t = t2)

[0082] p = p2(t = t2)

[0083] where: T2 is the reservoir temperature after well shut-in, K; p2 is the reservoir pressure after well shut-in, Pa;

[0084] ​For the stress field, since the reservoir is constrained by in-situ stresses and is in a state of having stresses but no strains in the initial state, this is achieved by applying boundary loads, that is,

[0085] σ ij (x, y, z, t = 0) = [σ v σ H σ h

[0086] In the formula, σ v is the vertical in-situ stress, Pa; σ H is the maximum horizontal in-situ stress, Pa; σ h is the minimum horizontal in-situ stress, Pa;

[0087] (2) Boundary conditions

[0088] During the flowback process, the inner boundary condition is the pressure boundary condition, and the pressure boundary condition is a constant-pressure boundary, that is,

[0089]

[0090] In the formula, A is the cross-sectional area, m 2 ; q w is the flow rate of the injected fracturing fluid, kg·m 3 / s, and p w is the bottom-hole flowing pressure.

[0091] The numerical simulation method for a fractured horizontal well in a high clay content shale oil reservoir in the embodiments of the present application has the following beneficial effects:

[0092] The present application provides a numerical simulation method for a fractured horizontal well in a high clay content shale oil reservoir, including: Step S1, obtaining reservoir physical property parameters according to the actual shale oil reservoir characteristics and core, logging, and fluid test data; obtaining the fracture geometric information of the shale oil reservoir, treating the fractures in the reservoir as the inner boundaries of the reservoir for dimensionality reduction, establishing a reservoir geometric model, and then performing geometric subdivision on the reservoir geometric model using triangular meshes to form discrete elements; Step S2, establishing a prediction model for reservoir development dynamic parameters of a fractured well in a shale oil reservoir based on the discrete fracture model; Step S3, performing prediction simulation of reservoir development dynamic parameters during the fracturing process; Step S4, performing prediction simulation of reservoir development dynamic parameters during the shut-in process; Step S5, performing prediction simulation of reservoir development dynamic parameters during the fracturing flowback process. The present invention comprehensively considers the coupling effects of osmotic pressure and heat flow solidification, predicts the fracturing and development dynamics of a high clay content shale oil reservoir, determines its main influencing factors, and further provides guidance for fracturing optimization of shale oil reservoirs. Description of the Drawings

[0093] The drawings of the present application are exemplary.

[0094] Figure 1 ​Schematic flow chart of the numerical simulation method for fractured horizontal wells in high-clay-content shale oil reservoirs in the embodiments of the present application;

[0095] Figure 2 Coupling relationship diagram in the numerical simulation method for fractured horizontal wells in high-clay-content shale oil reservoirs in the embodiments of the present application;

[0096] Figure 3 Schematic diagram of multi-cluster intensive segmented fracturing of horizontal wells;

[0097] Figure 4 Schematic diagram of the water saturation distribution of the solution result; Figure 5 Schematic diagram of the temperature distribution of the solution result; Figure 6 Schematic diagram of the pressure distribution of the solution result; Figure 7 Schematic diagram of the concentration distribution of the solution result. Specific implementation manners

[0098] The present application will be further introduced below in conjunction with the drawings and embodiments.

[0099] In the following description, the terms "first" and "second" are only for the purpose of description and cannot be construed as indicating or implying relative importance. The following description provides multiple embodiments of the present invention, and different embodiments can be replaced or combined. Therefore, the present application can also be considered to include all possible combinations of the same and / or different embodiments described. Thus, if one embodiment includes features A, B, and C, and another embodiment includes features B and D, then the present application should also be considered to include embodiments containing all other possible combinations of A, B, C, and D, even though such embodiments may not be explicitly described in the following content.

[0100] The following description provides examples and does not limit the scope, applicability, or examples set forth in the claims. Changes can be made to the functions and arrangements of the described elements without departing from the scope of the content of the present application. Various processes or components can be appropriately omitted, substituted, or added to each example. For example, the described method can be executed in a different order from the described order, and various steps can be added, omitted, or combined. In addition, the features described in some examples can be combined into other examples.

[0101] Embodiment 1

[0102] As Figure 1As shown in the figure, the numerical simulation method for fracturing horizontal wells in high clay content shale oil reservoirs of this application includes: Step S1, obtaining reservoir physical property parameters according to actual shale reservoir characteristics and core, logging, and fluid test data; obtaining the geometric information of fractures in the shale oil reservoir, taking the fractures in the reservoir as internal boundaries of the reservoir for dimensionality reduction processing, establishing a reservoir geometric model, and then using triangular meshes to perform geometric dissection on the reservoir geometric model to form discrete elements; Step S2, establishing a prediction model for reservoir development dynamic parameters of a fractured well in a shale oil reservoir based on a discrete fracture model; Step S3, performing prediction simulation of reservoir development dynamic parameters during the fracturing process; Step S4, performing prediction simulation of reservoir development dynamic parameters during the shut-in process; Step S5, performing prediction simulation of reservoir development dynamic parameters during the fracturing flowback process.

[0103] The present invention comprehensively considers the osmotic pressure and the coupled effect of heat flow solidification, predicts the fracturing and development dynamics of high clay content shale oil reservoirs, determines its main influencing factors, and further provides guidance for the optimization of shale reservoir fracturing.

[0104] Example 2

[0105] Step S1 obtaining reservoir physical property parameters according to actual shale reservoir characteristics and core, logging, and fluid test data includes: collecting reservoir geological parameters, including the distribution of maximum and minimum horizontal principal stresses, Young's modulus of rock, Poisson's ratio, formation temperature and pressure, porosity, permeability, and semipermeable membrane efficiency; collecting reservoir fluid parameters, including viscosity, salinity, thermal conductivity, specific heat capacity, relative permeability curve, and capillary force; collecting the completion information of the fractured well, including cluster spacing; collecting fracturing construction parameters, including construction displacement, fracturing fluid viscosity, salinity, thermal conductivity, and specific heat capacity; collecting the distribution of hydraulic fractures and natural fractures and the fracture length, width, porosity, and permeability; obtaining the geometric information of fractures in the shale oil reservoir according to the actual geological data of the reservoir or existing geological model data: including fracture location, size, density, strike, and aperture.

[0106] Step S2 establishing a prediction model for reservoir development dynamic parameters of a fractured well in a shale oil reservoir based on a discrete fracture model includes:

[0107] Step S21: Establishing the oil-water two-phase control equations for the matrix and fractures respectively:

[0108] The mass conservation equation of the fluid in the reservoir is: (1)

[0110] In the formula, S is the saturation of the fluid, dimensionless, φ is the porosity, dimensionless; Q is the source phase; a is the fluid type (oil / water). Both the porosity and the density can be written as functions of pressure:

[0111] ρ a =ρ a0 [1 + Cla (p - p0)] (2)

[0112] φ = φ0[C f (p - p0)] (3)

[0113] ρ a φ = ρ a0 φ0 + ρ a C a (p - p0) (4)

[0114] In the formula, C la is the fluid compressibility in the reservoir, 1 / Pa; C f is the fluid compressibility, 1 / Pa; C a is the comprehensive compressibility, C a = C f + C la ρ a0 .

[0115] Seepage is the seepage under the action of two pressures, namely the pressure difference and the osmotic pressure difference. The motion equation of fluid seepage considering the action of the osmotic pressure difference:

[0116]

[0117] In the formula, u is the fluid seepage velocity; k is the fluid permeability, m 2 , μ is the fluid viscosity, Pa·s; p is the formation pressure, Pa; C is the salt concentration, mol / m 3 , E op is the semi - permeable membrane efficiency, dimensionless.

[0118] Among them, the theoretical calculation formula of the osmotic pressure is:

[0119]

[0120] In the formula, V is the molar volume of water, taking 1.8×10 - 5 m3 / mol, R is the gas constant, taking 8.31×103 Pa·L(mol·K) -1 ; T is the temperature, K.

[0121] Fritz proposed that for an electrolyte with a cation - anion ratio of 1:1 (this application assumes that only Na+ and Cl - exist in the reservoir), the formula for the osmotic pressure can be simplified as:

[0122] π≈vRTC (7)

[0123] In the formula, v is the number of ions forming the solution, dimensionless; C is the concentration of the solution, mol / m 3 . Then, with the help of two auxiliary equations, the equation

[0124] Swa +S oa = 1 (8)

[0125] p ca = p oa -p wa (9)

[0126] Combining the above equations, the matrix aqueous phase control equation is as follows:

[0127]

[0128] The matrix oil phase equation is as follows:

[0129]

[0130] The fracture aqueous phase equation is as follows:

[0131]

[0132] The fracture oil phase equation is as follows:

[0133]

[0134] Among them, the subscript m represents the matrix, f represents the fracture, w represents the aqueous phase, and o represents the oil phase;

[0135] u is the fluid seepage velocity;

[0136] k is the fluid permeability, m 2 ;

[0137] k m is the matrix permeability, m 2 ;

[0138] k f is the matrix permeability, m 2 ;

[0139] k rw,m is the relative permeability of the matrix aqueous phase, dimensionless;

[0140] k ro,m is the relative permeability of the matrix oil phase, dimensionless;

[0141] k rw,f is the relative permeability of the fracture aqueous phase, dimensionless;

[0142] k ro,f is the relative permeability of the fracture oil phase, dimensionless;

[0143] S w,m is the matrix water saturation, dimensionless;

[0144] S w,f is the fracture water saturation, dimensionless;

[0145] P w is the aqueous phase pressure, Pa;

[0146] P o,m is the matrix oil phase pressure, Pa;

[0147] P o,f is the fracture oil phase pressure, Pa;

[0148] P c,m is the capillary force in the matrix, Pa;

[0149] P c,f is the capillary force in the fracture, Pa;

[0150] ρ o0 is the oil phase density at the initial reservoir pressure, kg / m 3 ;

[0151] ρ w0 is the aqueous phase density at the initial reservoir pressure, kg / m 3 ;

[0152] ρ o,m is the oil phase density in the matrix at pressure Po, kg / m 3 ;

[0153] ρ w,m is the aqueous phase density in the matrix at pressure Po, kg / m 3 ;

[0154] ρ f,m is the aqueous phase density in the fracture at pressure Po, kg / m 3 ;

[0155] ρ o,f is the oil phase density in the fracture at pressure Po, kg / m 3 ;

[0156] φ m is the matrix porosity, dimensionless;

[0157] φ f is the fracture porosity, dimensionless;

[0158] d f is the fracture width, mm;

[0159] C w is the comprehensive compressibility of the aqueous phase, and can be written as Cw = C f +ΦC l,w ;

[0160] C o is the comprehensive compressibility of the oil phase, and can be written as Co = C f +ΦCl,o ;

[0161] C l,w is the compressibility of the aqueous phase in the reservoir, 1 / Pa;

[0162] C l,o is the compressibility of the oil phase in the reservoir, 1 / Pa

[0163] C f is the rock compressibility, 1 / Pa;

[0164] Φf is the fracture porosity, dimensionless;

[0165] Φ m is the matrix porosity, dimensionless;

[0166] μ w is the fluid viscosity, Pa·s;

[0167] μ o is the fluid viscosity, Pa·s;

[0168] S w,m is the matrix water saturation, dimensionless;

[0169] S w,f is the fracture water saturation, dimensionless;

[0170] π is the osmotic pressure, Pa;

[0171] C is the salt concentration, mol / m 3 ;

[0172] C m is the salt ion concentration in the matrix, mol / m 3 ;

[0173] C f is the salt ion concentration in the fracture, mol / m 3 ;

[0174] V is the molar volume of water, taking 1.8×10 -5 m 3 / mol;

[0175] R is the gas constant, taking 8.31×10 3 Pa·L(mol·K) -1 ;

[0176] T is the temperature, K;

[0177] a Ⅰ represents the activity of low-salinity water, dimensionless;

[0178] a Ⅱ represents the activity of high-salinity water, dimensionless;

[0179] v is the number of ions constituting the solution, dimensionless;

[0180] F diff is the flux generated by diffusion, mol / (m 2 ·s);

[0181] E op is the efficiency of the semi-permeable membrane, dimensionless. The efficiency of an ideal semi-permeable membrane is 1, that is, it does not allow any substance to pass through;

[0182] D is the diffusion coefficient, m 2 / s;

[0183] F adv is the flux generated by convection, mol / (m 2 ·s);

[0184] Taking S w and P o as the dependent variables, established through the coefficient differential equation interface. Write (1) and (2) as the corresponding weak form equations, add through the weak contribution boundary, and regard the fracture as the internal boundary of the reservoir for simulation. The injection boundary condition is exactly reflected in the weak form equation, so it is also added through the weak contribution boundary condition.

[0185] For the coupling effect of the stress field and the seepage field, it is mainly reflected in the change of permeability caused by pore deformation. The variation law of porosity with stress is

[0186] And the permeability and porosity satisfy

[0187] k = k0(φ / φ0) 3 (15)

[0188] where φ0 is the porosity in the stress-free state, dimensionless; φ r is the residual porosity, dimensionless; k0 is the permeability in the stress-free state, m 2 . Among them α φ = 5×10 -8 Pa -1 .

[0189] Step S22: Establish the salt ion migration control equation.

[0190] Salt ion migration includes two parts: convection and diffusion. Salt ion diffusion is related to the concentration gradient. The general constitutive equation describing diffusion:

[0191]

[0192] In the formula, F diffis the flux generated by diffusion, mol / (m 2 ·s); E op is the efficiency of the semi-permeable membrane, dimensionless. The efficiency of an ideal semi-permeable membrane is 1, i.e., no substance is allowed to pass through; D is the diffusion coefficient, m 2 / s.

[0193] The ion migration caused by flow is:

[0194] F adv = Cu (17)

[0195] In the formula, F adv is the diffusion flux, mol / (m 2 ·s); u is the velocity of fluid flow, m / s.

[0196] The salt ion conservation equation in the matrix:

[0197]

[0198] The salt ion conservation equation in the crack:

[0199]

[0200] The above equations are described through the steady convection-diffusion equation interface, and the expressions of the corresponding terms are input into the equations.

[0201] Step S23: Establish the stress field control equation. The stress field control is mainly governed by the equilibrium equation, constitutive equation, and displacement equation. The constitutive equation of an isotropic linear thermoelastic material is:

[0202] σ′ ij = 2Gε ij + λε kk δ ij - K′α T Tδ ij - α B pδ ij (20)

[0203] In the formula, σ ij is the stress tensor, representing the stress state at a point, Pa; ε ij is the strain tensor, dimensionless; u i is the displacement tensor, m; λ is the Lame constant, Pa,; G is the shear modulus, Pa,; K′ is the bulk modulus of elasticity, Pa,; α B is the Biot coefficient, dimensionless;; F is the body force, N / m, p′ is the mean pore pressure, Pa;

[0204] p′ = S w p w + S op o = S w (p o - p c ) + S o p o = p o - S w p c (21)

[0205] The effective stress calculation formula proposed by Terzaghi and modified by Skempton is:

[0206] σ′ = σ - α B p (22)

[0207] The quasi-static equilibrium differential equation is:

[0208] σ′ ij，j + F i = 0 (23)

[0209] The relationship between strain and displacement is:

[0210]

[0211] where ε ij is the strain tensor; u i is the displacement tensor, m; F is the body force, Pa;

[0212] The stress field equation uses the solid mechanics interface, and corresponding parameters can be input;

[0213] Step S24: Establish the temperature field control equation. According to the law of conservation of energy, the temperature field control equation is:

[0214]

[0215]

[0216] where, (ρC p ) eff is the effective specific heat capacity of the rock mass, λ eff is the effective thermal conductivity of the rock mass, η eff is the effective heat convection coefficient of the fluid, C s , C o , C w are the specific heat capacities of the rock skeleton, oil, and water respectively, λ s , λ o , λ w are the thermal conductivities of the rock skeleton, oil, and water respectively.

[0217] For isotropic materials, pure shear deformation does not produce thermal effects, only volume changes Only then does it lead to the coupling effect between the strain field and the temperature field. Compression Exothermic, expansion Endothermic. It is equivalent to an additional heat source amount,

[0218]

[0219] The above equations are described through the porous medium heat transfer interface. First, parameters such as the effective specific heat capacity, effective thermal conductivity, and effective heat convection coefficient are defined, and then the corresponding parameters are input into the porous medium heat transfer interface for simulation.

[0220] Step S25: Solve using COMSOL Multiphysic software.

[0221] Use COMSOL Multiphysic software. This software is based on the finite element method and provides a fully coupled solution method for multi-physical fields. It can combine the heat flow solidification four fields into a unified coupled equation system for solution in the coupled analysis, and at the same time calculate the dependent variables of each independent field; use the PDE module to solve the oil-water two-phase flow problem, the porous medium heat transfer module to solve the temperature field problem, and solid mechanics to solve the stress field problem. COMSOL Multiphysic combines the three to form a differential equation system expressed in a general form to achieve the full coupling solution of the THMC four fields.

[0222] Step S3 includes: Step S31: Predict and simulate the reservoir development dynamic parameters during the fracturing process according to the mathematical model and boundary conditions; Step S32: Output the numerical simulation results, including reservoir temperature, pressure, water saturation, and salt concentration.

[0223] Step S31 includes: Assume that the initial moment is time 0, inject for time t1, shut in for time t2, and flowback for time t3;

[0224] (1) Initial conditions

[0225] The initial conditions of temperature and pressure during the fracturing process are the initial formation temperature and formation pressure of the reservoir, which are realized through the Dirichlet boundary condition, that is,

[0226] T = T0 (t = 0) (28)

[0227] p = p0 (t = 0) (29)

[0228] In the formula: T0 is the initial formation temperature, K; p0 is the initial formation pressure, Pa;

[0229] For the stress field, since the reservoir is constrained by in-situ stress and is in a state of having stress but no strain in the initial state, it is realized by adding boundary loads, that is,

[0230] σij (x, y, z, t = 0) = [σ v σ H σ h (30)

[0231] In the formula, σ v is the vertical in-situ stress, Pa; σ H is the maximum horizontal in-situ stress, Pa; σ h is the minimum horizontal in-situ stress, Pa;

[0232] (2) Boundary conditions

[0233] During the fracturing process, the inner boundary condition is the flow rate boundary condition, which is mentioned in the weak form, that is,

[0234]

[0235] In the formula, A is the cross-sectional area, m 2 ; q w is the flow rate of the injected fracturing fluid, kg·m 3 / s, p w is the bottom-hole flowing pressure; The boundary conditions of the temperature field and the seepage field are similar,

[0236]

[0237] In the formula, q is the heat source term caused by the injection well, W / m 2 ; C pw is the specific heat capacity of the injected water; T inj is the injection temperature, K.

[0238] Step S4 includes: Step S41: Predict and simulate the reservoir development dynamic parameters during the shut-in process according to the mathematical model and boundary conditions; Step S42: Output the numerical simulation results, including reservoir temperature, pressure, water saturation, and salt concentration.

[0239] Step S41 includes: (1) Initial conditions

[0240] The initial conditions of temperature and pressure during the shut-in process are the reservoir temperature and pressure at the end of fracturing, which are realized through the Dirichlet boundary condition, that is,

[0241] T = T1(t = t1) (33)

[0242] p = p1(t = t1) (34)

[0243] In the formula: T1 is the reservoir temperature at the end of fracturing, K; p1 is the reservoir pressure at the end of fracturing, Pa;

[0244] For the stress field, since the reservoir is constrained by in-situ stresses, its initial state is one of stress without strain, which is achieved by applying boundary loads, i.e.,

[0245] σ ij (x,y,z,t = 0) = [σ v σ H σ h (35)

[0246] In the formula, σ v is the vertical in-situ stress, Pa; σ H is the maximum horizontal in-situ stress, Pa; σ h is the minimum horizontal in-situ stress, Pa.

[0247] Step S5 includes: Step S51: Predict and simulate the reservoir development dynamic parameters during the fracturing flowback process according to the mathematical model and boundary conditions; Step S52: Output the numerical simulation results, including reservoir temperature, pressure, water saturation, and salt concentration.

[0248] Step S51 includes: Assume that the initial moment is moment 0, inject for time t1, shut in for time t2, and flow back for time t3; (1) Initial conditions

[0249] The initial conditions of temperature and pressure during the fracturing flowback process are the reservoir temperature and pressure at the end of shut-in, which are achieved through Dirichlet boundary conditions, i.e.,

[0250] T = T2(t = t2) (36)

[0251] p = p2(t = t2) (37)

[0252] In the formula: T2 is the reservoir temperature after shut-in, K; p2 is the reservoir pressure after shut-in, Pa;

[0253] For the stress field, since the reservoir is constrained by in-situ stresses, its initial state is one of stress without strain, which is achieved by applying boundary loads, i.e.,

[0254] σ ij (x,y,z,t = 0) = [σ v σ H σ h (38)

[0255] In the formula, σ v is the vertical in-situ stress, Pa; σ H is the maximum horizontal in-situ stress, Pa; σ h is the minimum horizontal in-situ stress, Pa;

[0256] (2) Boundary conditions

[0257] During the flowback process, the inner boundary condition is the pressure boundary condition, and the pressure boundary condition is the constant pressure boundary, that is,

[0258]

[0259] In the formula, A is the cross-sectional area, m 2 ; q w is the flow rate of the injected fracturing fluid, kg·m 3 / s, p w is the bottom-hole flowing pressure.

[0260] The present invention comprehensively considers the osmotic pressure and the coupled effect of heat flow curing, predicts the fracturing and development dynamics of shale oil reservoirs with high clay content, determines its main influencing factors, and further provides guidance for the optimization of shale oil reservoir fracturing.

[0261] Example 3

[0262] Taking Well X in the shale oil reservoir as an example, the specific reservoir geology, engineering parameters, and fluid parameters are shown in Table 1. Horizontal well multi-cluster dense segmented fracturing is adopted, and one cluster of fractures is selected for simulation. The process of simulating 100 min of pumping injection, shutting in for 30 days, and flowing back (producing) for 5 years is carried out.

[0263] Table 1 Model parameters

[0264]

[0265] The capillary force is generally obtained directly through experiments or by combining experiments with empirical formulas. The empirical formula is as follows and can be obtained by measuring the interfacial tension and the porosity and permeability of the rock.

[0266]

[0267]

[0268] At different pressures and temperatures, the density, viscosity, specific heat capacity, and thermal conductivity of oil and water will also be different. The relevant data can be obtained by: ① researching relevant literature, ② fitting formulas through experiments and relevant software. Here, empirical formulas are given for reference.

[0269] μ w = 1.3799 - 0.0212T + 1.3604×10 -4 T 2 - 4.6454×10 -7 T 3 + 8.9043×10 -10 T 4 - 9.0791×10 -13 T 5 + 3.8457×10 -16 T 6, T ∈ [273, 413]

[0270] μ w = 0.0040 - 2.1075×10 -5 T + 3.8577×10 -8 T 2 - 2.3973×10 -11 T 3 , T ∈ [413, 553]

[0271] C pw = 12010 - 80.4×T + 0.3×T 2 - 5.4×10 -4 T 3 + 3.6×10 -7 T 4 , T ∈ [273, 553]

[0272] ρ w = 838.4661 + 1.4005T - 0.0030T 2 + 3.7182×10 -7 T 3 , T ∈ [273, 553]

[0273] λ w = - 0.8691 + 0.0089T - 1.5837×10 -5 T 2 + 7.9754×10 -9 T 3 , T ∈ [273, 553]

[0274] μ w = 91.45245 - 1.33227T + 0.00778T 2 - 2.27278×10 -5 T 3 + 3.32420×10 -8 T 4 - 1.94631×10 -11 T 5 , T ∈ [273, 453]

[0275] C pw = - 13408.1491 + 123.04415T - 0.33540T 2 + 3.125×10 -4 T 3 , T ∈ [273, 453]

[0276] ρ w = 1055.04607 - 0.58175T - 6.40532×10-5 T 2 , T ∈ [273, 453]

[0277] λ w = 0.134299 - 8.049738×10 -5 T, T ∈ [273, 453]

[0278] Based on Step 1, establish a physical model of single-cluster fracturing for X shale oil wells according to the parameters in Table 1 and the geological model.

[0279] Based on Step 2, realize the two-phase flow of oil and water through the coefficient differential equation of the PDE (partial differential equation) module. Regard the fracture as the internal boundary of the whole domain. Before solving, characterize the fracture equation by adding the weak form of the corresponding control equation in the fracture area, then the two-phase flow of oil and water with discrete fractures can be realized; solve the salt ion transport equation through the convection-diffusion equation of the PDE (partial differential equation) module; solve the temperature field through the heat transfer module of porous media. Solve the stress field and displacement field through the solid mechanics module. The coupling relationship is realized through parameter variable transfer and adding source phases.

[0280] Based on Step 3, as Figures 4 - 7 shown, the initial conditions during the fracturing process are the initial pressure, temperature and salt concentration of the reservoir. The flow field uses the flow boundary condition, the temperature field uses the heat flux boundary condition, the stress field uses the boundary load to simulate the constraint of in-situ stress on the reservoir, and the concentration field uses the Dirichlet boundary condition. The solution results are as Figure 4 (a), 5(a), 6(a), 7(a).

[0281] Based on Step 4, as Figures 4 - 7 shown, the initial conditions during the well shut-in process are the reservoir temperature, pressure and salt concentration at the end of fracturing. By adding research steps, use the end time of the first research as the initial time of this research for simulation. The solution results are as Figure 4 (b), 5(b), 6(b), 7(b).

[0282] Based on Step 5, as Figures 4 - 7 shown, the initial conditions during the flowback process are the reservoir temperature, pressure and salt concentration at the end of well shut-in. By adding research steps, use the end time of the second research as the initial time of this research for simulation. The flow field uses the pressure boundary condition, and the stress field uses the boundary load to simulate the constraint of in-situ stress on the reservoir. The solution results are as Figure 4 (c), 5(c), 6(c), 7(c).

[0283] The above introduction is only the preferred embodiment of the present invention and is not intended to limit the present invention. For those skilled in the art, various modifications and variations can be made to the present invention. Any modification, equivalent replacement, improvement, etc. made within the spirit and principle of the present invention shall be included within the protection scope of the present invention.

Claims

1. A numerical simulation method for fracturing horizontal wells in shale oil reservoirs with high clay content, characterized in that include: Step S1, obtaining reservoir physical property parameters according to actual shale reservoir characteristics and core, well logging and fluid testing data; Obtain the geometric information of shale oil reservoir fractures, use the fractures in the reservoir as the internal boundaries of the reservoir for dimensionality reduction, establish a reservoir geometric model, and then use triangular meshes to geometrically divide the reservoir geometric model into discrete units; Step S2, establishing a shale oil reservoir fracturing well reservoir development dynamic parameter prediction model based on a discrete fracture model; comprising: Step S21: Establish the matrix and fracture oil-water two-phase control equations respectively: The governing equation for the matrix water phase is: The matrix oil phase equation is: The fracture water phase equation is: The fracture oil phase equation is: where k m is the matrix permeability; k f is the fracture permeability; k rwm Matrix aqueous relative permeability; k rom Relative permeability of matrix oil phase; k rwf Relative permeability of the fracture water phase; k rof Relative permeability of fissure oil phase; S wm is the matrix water saturation; S wf is the fracture water saturation; P om is the matrix oil phase pressure; P of is the oil phase pressure of the fracture; ρ om is the density of the oil phase in the matrix at pressure Po; ρ wm is the density of the aqueous phase in the matrix at pressure Po; ρ of is the density of the oil phase in the fracture under the pressure Po; d f is the seam width; C w is the comprehensive compressibility coefficient of the aqueous phase; C o is the comprehensive compressibility coefficient of the oil phase; φ f is the fracture porosity; φ m is the matrix porosity; μ w is the viscosity of the aqueous fluid; μ o is the viscosity of the oil-phase fluid; π is the osmotic pressure; C is the salt concentration; E op is the efficiency of the semi-permeable membrane; Step S22: establishing a salt ion transport control equation; Step S23: Establish the stress field control equation. The stress field control is mainly controlled by the equilibrium equation, the constitutive equation, and the displacement equation. The constitutive equation of the isotropic linear thermoelastic material is: σ′ ij = 2Gε ij + λε kk δ ij - K′α T Tδ ij - α B pδ ij where, σ′ ij is the stress tensor, representing the stress state at a point, in Pa; ε ij is the strain tensor, dimensionless; λ is the Lame constant, in Pa; G is the shear modulus, in Pa; K′ is the bulk modulus of elasticity, in Pa; α B is the Biot coefficient, dimensionless; p is the mean pore pressure, in Pa; p′ = S w p w + S o p o = S w (p o - p c )+ S o p o = p o - S w p c Proposed by Terzaghi, the Skempton modified effective stress calculation formula is: σ′ = σ - α B p The quasi-static equilibrium differential equation is: σ′ ij,j +F i =0 The relationship between strain and displacement is: where ε ij is the strain tensor; The stress field equation uses the solid mechanics interface and only requires inputting the corresponding parameters. Step S24: Establish the temperature field control equation. According to the law of conservation of energy, the temperature field control equation is: Where, (ρC p ) eff is the effective specific heat capacity of the rock mass, λ eff is the effective thermal conductivity of the rock mass, η eff is the effective heat convection coefficient of the fluid, C s , C o , C w are the specific heat capacities of the rock skeleton, oil, and water respectively, λ s , λ o , λ w are the thermal conductivities of the rock skeleton, oil, and water respectively; Step S25: solving using COMSOL Multiphysic software; Step S3, performing a reservoir development dynamic parameter prediction simulation during the fracturing process; Step S4, performing a prediction simulation of reservoir development dynamic parameters during the well shut-in process; Step S5, performing a prediction simulation of reservoir development dynamic parameters during the fracturing flowback process.

2. The numerical simulation method for fracturing horizontal wells in high clay content shale oil reservoirs according to claim 1, wherein Step S1 includes: The collected reservoir geological parameters include the maximum and minimum horizontal principal stress distribution, rock Young's modulus, Poisson's ratio, formation temperature and pressure, porosity, permeability and semi-permeable membrane efficiency; the collected reservoir fluid parameters include viscosity, salinity, thermal conductivity, specific heat capacity, phase permeability curve and capillary force; the collected completion information of fracturing wells, including cluster spacing; the collected fracturing construction parameters, including construction displacement, fracturing fluid viscosity, salinity, thermal conductivity and specific heat capacity; the collected distribution of hydraulic fractures and natural fractures, as well as the fracture length, fracture width and pore permeability; the obtained shale oil reservoir fracture geometry information based on the actual geological data of the reservoir or the existing geological model data: including fracture location, size, density, direction and aperture.

3. The numerical simulation method for fractured horizontal wells in high clay content shale oil reservoirs according to claim 1, wherein Step S25 includes: COMSOL Multiphysic software is used. This software is based on the finite element method and provides a fully coupled solution method for multiple physical fields. In the coupled analysis, the four fields of heat flow and solidification can be combined to form a unified set of coupled equations for solution, and the dependent variables of each independent field can be calculated at the same time. The PDE module is used to solve the oil-water two-phase flow problem, the porous medium heat transfer module is used to solve the temperature field problem, and the solid mechanics is used to solve the stress field problem. COMSOL Multiphysic combines the three to form a set of differential equations expressed in a general formula, realizing the fully coupled solution of the four fields of THMC.

4. The numerical simulation method for fractured horizontal wells in high clay content shale oil reservoirs according to any one of claims 1-3, characterized in that, Step S3 includes: Step S31: Predict and simulate the reservoir development dynamic parameters during the fracturing process according to the mathematical model and boundary conditions; Step S32: Output the numerical simulation results, including reservoir temperature, pressure, water saturation, and salt concentration.

5. The numerical simulation method for fracturing horizontal wells in high clay content shale oil reservoirs according to claim 4, characterized in that Step S31 includes: Assume that the initial time is time 0, the injection time is t1, the shut-in time is t2, and the flowback time is t3; (1) Initial conditions The initial conditions of temperature and pressure during the fracturing process are the initial formation temperature and formation pressure of the reservoir, which are realized through the Dirichlet boundary conditions, that is, T = T0 (t = 0) p = p0 (t = 0) where: T0 is the initial formation temperature, K; p0 is the initial formation pressure, Pa; For the stress field, since the reservoir is constrained by in-situ stresses and is in a state of having stress but no strain in the initial state, it is realized by adding boundary loads, that is, σ ij (x, y, z, t = 0) = [σ v σ H σ h ​ where σ v is the vertical in-situ stress, Pa; σ H is the maximum horizontal in-situ stress, Pa; σ h is the minimum horizontal in-situ stress, Pa; (2) Boundary conditions The inner boundary condition during the fracturing process is the flow rate boundary condition, which is mentioned in the weak form, that is, where A is the cross-sectional area, m 2 ; p w is the bottom-hole flowing pressure; The boundary conditions of the temperature field and the seepage field are similar. where q w，T is the heat source term caused by the injection well, W / m 2 ; C pw is the specific heat capacity of the injected water; T inj is the injection temperature, K.

6. The numerical simulation method for a fractured horizontal well in a high clay content shale oil reservoir according to any one of claims 1-3, characterized in that, Step S4 includes: Step S41: Predict and simulate the reservoir development dynamic parameters during the shut-in process according to the mathematical model and boundary conditions; Step S42: Output the numerical simulation results, including reservoir temperature, pressure, water saturation, and salt concentration.

7. The numerical simulation method for fracturing horizontal wells in high clay content shale oil reservoirs according to claim 6, characterized in that, Step S41 includes: (1) Initial conditions The initial conditions of temperature and pressure during the shut-in process are the reservoir temperature and pressure at the end of fracturing, which are realized through the Dirichlet boundary conditions, that is, T = T1 (t = t1) p = p1 (t = t1) where: T1 is the reservoir temperature at the end of fracturing, K; p1 is the reservoir pressure at the end of fracturing, Pa; For the stress field, since the reservoir is constrained by in-situ stresses and is in a state of having stress but no strain in the initial state, it is realized by adding boundary loads, that is, σ ij (x, y, z, t = 0) = [σ v σ H σ h ​ where, σ v is the vertical in-situ stress, Pa; σ H is the maximum horizontal in-situ stress, Pa; σ h is the minimum horizontal in-situ stress, Pa.

8. The numerical simulation method for fractured horizontal wells in high clay content shale oil reservoirs according to any one of claims 1 to 3, characterized in that, Step S5 includes: Step S51: Predict and simulate the reservoir development dynamic parameters during the fracturing flowback process according to the mathematical model and boundary conditions; Step S52: Output the numerical simulation results, including reservoir temperature, pressure, water saturation, and salt concentration.

9. The numerical simulation method for fracturing horizontal wells in high clay content shale oil reservoirs according to claim 8, wherein Step S51 includes: Assume that the initial time is time 0, the injection time is t1, the shut-in time is t2, and the flowback time is t3; (1) Initial conditions The initial conditions of temperature and pressure during the fracturing flowback process are the reservoir temperature and pressure at the end of shut-in, which are realized through the Dirichlet boundary conditions, that is, T = T2 (t = t2) p = p2 (t = t2) where: T2 is the reservoir temperature after shut-in, K; p2 is the reservoir pressure after shut-in, Pa; For the stress field, since the reservoir is constrained by in-situ stresses and is in a state of having stress but no strain in the initial state, it is realized by adding boundary loads, that is, σ ij (x,y,z,t=0)=[σ v σ H σ h ​ where, σ v is the vertical in-situ stress, Pa; σ H is the maximum horizontal in-situ stress, Pa; σ h is the minimum horizontal in-situ stress, Pa; (2) Boundary conditions The inner boundary condition during the flowback process is the pressure boundary condition, and the pressure boundary condition is a constant pressure boundary, that is, where p w is the bottom-hole flowing pressure.

Images