CN114880895A - Numerical simulation method for high-clay-content shale oil reservoir fracturing horizontal well - Google Patents

Abstract

The invention provides a numerical simulation method for a high-clay-content shale oil reservoir fracturing horizontal well, which comprises the following steps of: step S1, acquiring reservoir related parameters according to actual shale oil reservoir characteristics and rock core, well logging and fluid test data; step S2, establishing a shale oil reservoir fracturing well reservoir development dynamic parameter prediction model based on the discrete fracture model; step S3, reservoir development dynamic parameter prediction simulation of the fracturing process is carried out; step S4, performing reservoir development dynamic parameter prediction simulation in the well closing process; and step S5, performing reservoir development dynamic parameter prediction simulation in the fracturing flow-back process. According to the method, the osmotic pressure and heat flow solidification coupling effect are comprehensively considered, the fracturing and development dynamics of the shale oil reservoir with high clay content are predicted, the main influence factors of the fracturing and development dynamics are determined, and guidance is provided for shale oil reservoir fracturing optimization.

Description

Numerical simulation method for high-clay-content shale oil reservoir fracturing horizontal well

Technical Field

The invention belongs to the field of oil and gas field development, and particularly relates to a numerical simulation method for a shale oil reservoir fracturing horizontal well with high clay content.

Background

For shale oil reservoirs containing high clay minerals, the mineralization degree of formation water can reach hundreds of thousands of ppm, the osmotic pressure effect formed by the formation water and injected low-mineralization fracturing fluid in the fracturing process is obvious, and great influence is caused on fracturing and flowback. The shale oil reservoir fracturing process has a temperature (T), seepage (H), stress (M) and chemical (C) mutual coupling process, fluid in the reservoir flows in a porous medium and a crack, distribution of a temperature field can be influenced, and the change of the temperature field causes physical property change of the fluid and counteracts the seepage field. Also cause variations in osmotic pressure; the fluid movement can cause the convection action of salt ions, and the convection diffusion of the salt ions can influence the magnitude of osmotic pressure due to the action of the shale semipermeable membrane, so that the fluid movement speed is changed. Due to the constraint action of the ground stress, the volume expands and contracts to cause thermal stress, so that the stress field changes, the energy of the reservoir rock body changes due to the change of the stress field, and the temperature field changes according to the energy conservation principle. In the process of injection or exploitation, the change of pore pressure can cause the deformation of an oil reservoir rock framework, so that pore permeability is increased or reduced, and a seepage field is influenced. Therefore, it is necessary to comprehensively consider the osmotic pressure and heat flow solidification coupling effect, predict the fracturing and development dynamics of the shale oil reservoir with high clay content, determine the main influence factors of the fracturing and development dynamics, and further provide guidance for shale oil reservoir fracturing optimization.

Disclosure of Invention

According to the method, the osmotic pressure and heat flow solidification coupling effect are comprehensively considered, the fracturing and development dynamics of the shale oil reservoir with high clay content are predicted, the main influence factors of the fracturing and development dynamics are determined, and guidance is provided for shale oil reservoir fracturing optimization.

In a first aspect, an embodiment of the present application provides a numerical simulation method for a high clay content shale oil reservoir fractured horizontal well, including:

step S1, obtaining physical property parameters of the reservoir according to the actual shale oil reservoir characteristics and the rock core, well logging and fluid testing data; acquiring shale oil reservoir fracture geometric information, performing dimensionality reduction treatment on fractures in an oil reservoir as an inner boundary of the oil reservoir, establishing an oil reservoir geometric model, and then performing geometric subdivision on the oil reservoir geometric model by adopting triangular meshes to form discrete units;

step S2, establishing a shale oil reservoir fracturing well reservoir development dynamic parameter prediction model based on a discrete fracture model;

step S3, reservoir development dynamic parameter prediction simulation of the fracturing process is carried out;

step S4, performing reservoir development dynamic parameter prediction simulation in the well closing process;

and step S5, performing reservoir development dynamic parameter prediction simulation in the fracturing flow-back process.

Wherein, step S1 includes:

collected reservoir geological parameters including maximum and minimum horizontal principal stress distribution, rock Young's modulus, Poisson's ratio, formation temperature and pressure, porosity, permeability, semi-permeable membrane efficiency; collecting reservoir fluid parameters including viscosity, mineralization degree, heat conductivity coefficient, specific heat capacity, phase permeation curve and capillary force; collected completion information for fractured wells, including cluster spacing; collecting fracturing construction parameters including construction discharge capacity, viscosity of fracturing fluid, mineralization degree, heat conductivity coefficient and specific heat capacity; collecting the distribution of hydraulic fractures and natural fractures, and the length, width and seepage of the fractures; acquiring the geometrical information of the shale oil reservoir fracture according to the actual geological data or the existing geological model data of the reservoir: including crack position, size, density, strike, opening.

Wherein, step S2 includes:

step S21: respectively establishing an oil-water two-phase control equation of a matrix and a crack:

the matrix water phase control equation is:

the matrix oil phase equation is:

the fracture water phase equation is:

the fracture oil phase equation is:

step S22: establishing a salt ion migration control equation;

step S23: establishing a stress field control equation, wherein the stress field control is mainly controlled by a balance equation, a constitutive equation and a displacement equation, and the constitutive equation of the isotropic linear thermoelastic material is as follows:

σ′ _ij ＝2Gε _ij +λε _kk δ _ij -K′α _T Tδ _ij -α _B pδ _ij

in the formula, σ _ij Is a stress tensor, representing the stress state of a point, Pa; epsilon _ij Is a strain tensor, dimensionless; u. of _i Is the displacement tensor, m; λ is Lame constant, Pa,; g is shear modulus, Pa,; k' is the bulk modulus of elasticity, Pa,; alpha is alpha _B Biot coefficient, dimensionless; (ii) a F is the volume force, N/m, p' is the average pore pressure, 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

the Skempton corrected effective stress calculation formula proposed by Terzaghi is:

σ′＝σ-α _B p

the quasi-static equilibrium differential equation is:

σ′ _ij,j +F _i ＝0

the relationship between strain and displacement is:

in the formula of _ij Is a strain tensor; u. of _i Is the displacement tensor, m; f is the volume force, Pa;

the stress field equation uses a solid mechanics interface and corresponding parameters are input;

step S24: establishing a temperature field control equation, and obtaining the temperature field control equation according to the law of conservation of energy as follows:

wherein (ρ C) _p ) _eff Is the effective specific heat capacity, lambda, 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 Respectively the specific heat capacity, lambda, of the rock skeleton, oil and water _s ，λ _o ，λ _w The thermal conductivity coefficients of the rock skeleton, oil and water are respectively;

step S25: the solution was performed using COMSOL Multiphysic software.

Wherein, step S25 includes:

the COMSOL Multiphysic software is used and is based on a finite element method, a full-coupling solving method is provided for multiple physical fields, four solidified fields of heat flow can be combined to form a uniform coupling equation set for solving in coupling analysis, and dependent variables of the independent fields are obtained through calculation; the method is characterized in that a PDE module is used for solving the problem of oil-water two-phase flow, a porous medium heat transfer module is used for solving the problem of temperature field, solid mechanics is used for solving the problem of stress field, COMSOL Multiphysic combines the three to form a differential equation set expressed by a general formula, and THMC four-field full-coupling solution is realized.

Wherein, step S3 includes:

step S31: performing reservoir development dynamic parameter prediction simulation of the fracturing process according to the mathematical model and the boundary conditions;

step S32: and outputting numerical simulation results including reservoir temperature, pressure, water saturation and salt concentration.

Wherein, step S31 includes:

assuming that the initial time is 0 time, the injection time is t1 time, the shut-in time is t2 time, and the flow-back time is t3 time;

(1) initial conditions

The initial conditions of temperature and pressure in the fracturing process are the initial formation temperature and formation pressure of a reservoir, and are realized through the Dielder boundary condition, namely,

T＝T ₀ (t＝0)

p＝p ₀ (t＝0)

in the formula: t is ₀ Initial formation temperature, K; p is a radical of ₀ Initial formation pressure, Pa;

for the stress field, the initial state is in a stressed and unstrained state due to the reservoir being constrained by the geostress, which is achieved by adding boundary loads, i.e.,

σ _ij (x,y,z,t＝0)＝[σ _v σ _H σ _h ]

in the formula, σ _v Is vertical ground stress, Pa; sigma _H Pa for horizontal maximum ground stress; sigma _h Is the horizontal minimum ground stress, Pa;

(2) boundary condition

The boundary condition within the fracturing process is the flow boundary condition, which is mentioned in weak form, i.e.,

wherein A is the cross-sectional area, m ² ；q _w For the flow rate of the injected fracturing fluid, kg.m ³ /s，p _w Is the bottom hole flowing pressure;

the boundary conditions of the temperature field and the seepage field are similar,

wherein q is a heat source term due to an injection well, W/m ² ；C _pw The specific heat capacity of the injected water; t is _inj Is the injection temperature, K.

Wherein, step S4 includes:

step S41: performing reservoir development dynamic parameter prediction simulation of the well shut-in process according to the mathematical model and the boundary conditions;

step S42: and outputting numerical simulation results including reservoir temperature, pressure, water saturation and salt concentration.

Wherein, step S41 includes:

(1) initial conditions

The initial conditions of temperature and pressure in the shut-in process are the reservoir temperature and pressure at the end of fracturing, which are realized by the Dielder boundary conditions, that is,

T＝T ₁ (t＝t ₁ )

p＝p ₁ (t＝t ₁ )

in the formula: t is ₁ Is the reservoir temperature at the end of fracturing, K; p is a radical of ₁ Reservoir pressure at the end of fracturing, Pa;

for the stress field, the initial state is in a stressed and unstrained state due to the reservoir being constrained by the geostress, which is achieved by adding boundary loads, i.e.,

σ _ij (x,y,z,t＝0)＝[σ _v σ _H σ _h ]

in the formula, σ _v Is vertical ground stress, Pa; sigma _H Pa for horizontal maximum ground stress; sigma _h Is the horizontal minimum stress, Pa.

Wherein, step S5 includes:

step S51: performing reservoir development dynamic parameter prediction simulation of the fracturing flow-back process according to the mathematical model and the boundary conditions;

step S52: and outputting numerical simulation results including reservoir temperature, pressure, water saturation and salt concentration.

Wherein, step S51 includes:

assuming that the initial time is 0 time, the injection time is t1 time, the shut-in time is t2 time, and the flow-back time is t3 time;

(1) initial conditions

The initial conditions of the temperature and the pressure in the fracturing flow-back process are the reservoir temperature and the pressure at the end of well shut-in, and are realized through the Dielder boundary condition, namely,

T＝T ₂ (t＝t ₂ )

p＝p ₂ (t＝t ₂ )

in the formula: t is ₂ The reservoir temperature after the shut-in is finished, K; p is a radical of ₂ Reservoir pressure, Pa, after shut-in is finished;

for the stress field, the initial state is in a stressed and unstrained state due to the reservoir being constrained by the geostress, which is achieved by adding boundary loads, i.e.,

σ _ij (x,y,z,t＝0)＝[σ _v σ _H σ _h ]

in the formula, σ _v Is vertical ground stress, Pa; sigma _H Pa for horizontal maximum ground stress; sigma _h Is the horizontal minimum ground stress, Pa;

(2) boundary condition

The boundary condition in the flow-back process is a pressure boundary condition, the pressure boundary condition is a constant pressure boundary, that is,

wherein A is the cross-sectional area, m ² ；q _w For the flow rate of the injected fracturing fluid, kg.m ³ /s，p _w Is the bottom hole flowing pressure.

The numerical simulation method for the high-clay-content shale oil reservoir fracturing horizontal well has the following beneficial effects:

the application provides a numerical simulation method for a high-clay-content shale oil reservoir fracturing horizontal well, which comprises the following steps of: step S1, obtaining reservoir physical property parameters according to the actual shale oil reservoir characteristics and the rock core, well logging and fluid test data; acquiring shale oil reservoir fracture geometric information, performing dimensionality reduction treatment on fractures in an oil reservoir as an inner boundary of the oil reservoir, establishing an oil reservoir geometric model, and then performing geometric subdivision on the oil reservoir geometric model by adopting triangular meshes to form discrete units; step S2, establishing a shale oil reservoir fracturing well reservoir development dynamic parameter prediction model based on the discrete fracture model; step S3, reservoir development dynamic parameter prediction simulation of the fracturing process is carried out; step S4, performing reservoir development dynamic parameter prediction simulation in the well closing process; and step S5, performing reservoir development dynamic parameter prediction simulation in the fracturing flow-back process. According to the method, the osmotic pressure and heat flow solidification coupling effect are comprehensively considered, the fracturing and development dynamics of the shale oil reservoir with high clay content are predicted, the main influence factors of the fracturing and development dynamics are determined, and guidance is provided for shale oil reservoir fracturing optimization.

Drawings

The drawings of the present application are illustrative.

FIG. 1 is a flow diagram of a numerical simulation method for fracturing a horizontal well in a shale oil reservoir with high clay content according to an embodiment of the present application;

FIG. 2 is a coupling relation diagram in the numerical simulation method for fracturing a horizontal well in a shale oil reservoir with high clay content in the embodiment of the application;

FIG. 3 is a schematic diagram of multi-cluster intimate-cutting staged fracturing of a horizontal well;

FIG. 4 is a schematic diagram of the water saturation distribution of the solution; FIG. 5 is a schematic diagram of the temperature distribution of the solution; FIG. 6 is a schematic diagram of the pressure distribution of the solution; fig. 7 is a schematic view of the concentration distribution obtained by the solution.

Detailed Description

The present application is further described with reference to the following figures and examples.

In the following description, the terms "first" and "second" are used for descriptive purposes only and are not intended to indicate or imply relative importance. The following description provides embodiments of the invention, which may be combined or substituted for various embodiments, and this application is therefore intended to cover all possible combinations of the same and/or different embodiments described. Thus, if one embodiment includes feature A, B, C and another embodiment includes feature B, D, then this application should also be construed to include embodiments that include A, B, C, D in all other possible combinations, even though such embodiments may not be explicitly recited in the text that follows.

The following description provides examples, and does not limit the scope, applicability, or examples set forth in the claims. Changes may be made in the function and arrangement of elements described without departing from the scope of the disclosure. Various examples may omit, substitute, or add various procedures or components as appropriate. For example, the described methods may be performed in an order different than the order described, and various steps may be added, omitted, or combined. Furthermore, features described with respect to some examples may be combined into other examples.

Example one

As shown in fig. 1, the numerical simulation method for the fractured horizontal well of the high-clay-content shale oil reservoir comprises the following steps: step S1, obtaining reservoir physical property parameters according to the actual shale oil reservoir characteristics and the rock core, well logging and fluid test data; acquiring shale oil reservoir fracture geometric information, performing dimensionality reduction treatment on fractures in an oil reservoir as an inner boundary of the oil reservoir, establishing an oil reservoir geometric model, and then performing geometric subdivision on the oil reservoir geometric model by adopting triangular meshes to form discrete units; step S2, establishing a shale oil reservoir fracturing well reservoir development dynamic parameter prediction model based on the discrete fracture model; step S3, reservoir development dynamic parameter prediction simulation of the fracturing process is carried out; step S4, performing reservoir development dynamic parameter prediction simulation in the well closing process; and step S5, performing reservoir development dynamic parameter prediction simulation in the fracturing flow-back process.

According to the method, the osmotic pressure and heat flow solidification coupling effect are comprehensively considered, the fracturing and development dynamics of the shale oil reservoir with high clay content are predicted, the main influence factors of the fracturing and development dynamics are determined, and guidance is provided for shale oil reservoir fracturing optimization.

Example two

Step S1, obtaining reservoir physical property parameters according to the actual shale oil reservoir characteristics and the rock core, well logging and fluid test data comprises: collected reservoir geological parameters including maximum and minimum horizontal principal stress distribution, rock Young's modulus, Poisson's ratio, formation temperature and pressure, porosity, permeability, semi-permeable membrane efficiency; collecting reservoir fluid parameters including viscosity, mineralization degree, heat conductivity coefficient, specific heat capacity, phase permeation curve and capillary force; collected completion information for fractured wells, including cluster spacing; collecting fracturing construction parameters including construction discharge capacity, viscosity of fracturing fluid, mineralization degree, heat conductivity coefficient and specific heat capacity; collecting the distribution of hydraulic fractures and natural fractures, and the length, width and seepage of the fractures; acquiring the geometrical information of the shale oil reservoir fracture according to the actual geological data or the existing geological model data of the reservoir: including crack position, size, density, strike, opening.

Step S2, the establishment of a shale oil reservoir fractured well reservoir development dynamic parameter prediction model based on the discrete fracture model comprises the following steps:

step S21: respectively establishing an oil-water two-phase control equation of a matrix and a crack:

the conservation of mass equation for the fluid in the reservoir is:

(1)

in the formula, S is the saturation of the fluid and is dimensionless, and phi is the porosity and is dimensionless; q is a source; a is fluid type (oil/water). Porosity and density can both be written as a function of pressure:

ρ _a ＝ρ _a0 [1+C _la (p-p ₀ )] (2)

φ＝φ ₀ [C _f (p-p ₀ )] (3)

ρ _a φ＝ρ _a0 φ ₀ +ρ _a C _a (p-p ₀ ) (4)

in the formula, C _la Is the fluid compressibility in the reservoir, 1/Pa; c _f Is the fluid compression coefficient, 1/Pa; c _a To synthesize a compression factor, C _a ＝C _f +C _la ρ _a0 。

Seepage is seepage under the action of two pressures of pressure difference and osmotic pressure difference, and the motion equation of fluid seepage under the action of osmotic pressure difference is considered:

wherein u is the fluid seepage velocity; k is the fluid permeability, m ² μ is the fluid viscosity, pas; p is the formation pressure, Pa; c is the salt concentration, mol/m ³ ，E _op Is the efficiency of the semipermeable membrane, and has no dimension.

The theoretical calculation formula of the osmotic pressure is as follows:

wherein V is the molar volume of water and is 1.8X 10-5m3/mol, R is the gas constant and is 8.31X 103 Pa.L (mol. K) ^-1 (ii) a T is the temperature, K.

Fritz teaches that for electrolytes with a 1:1 ratio of positive to negative ions (assuming only Na + and Cl-' are present in the reservoir), the formula for osmolarity can be simplified as:

π≈vRTC (7)

wherein v is the number of ions constituting the solution, dimensionless; c is the concentration of the solution, mol/m ³ . Then, the equations can be solved by two auxiliary equations

S _wa +S _oa ＝1 (8)

p _ca ＝p _oa -p _wa (9)

Combining the above equations, the matrix water phase control equation is:

the matrix oil phase equation is:

the fracture water phase equation is:

the fracture oil phase equation is:

wherein the subscript m represents the matrix, f represents the fracture, w represents the water phase, o represents the oil phase;

u is the fluid seepage velocity;

k is the fluid permeability, m ² ；

k _m As the permeability of the matrix, m ² ；

k _f As the permeability of the matrix, m ² ；

k _rw,m Relative permeability of matrix water phase, dimensionless;

k _ro,m relative permeability of the matrix oil phase, dimensionless;

k _rw,f relative permeability of water phase of the crack and no dimension;

k _ro,f relative permeability of the fracture oil phase and no dimension;

S _w,m is the water saturation of the substrate, and has no dimension;

S _w,f the fracture water saturation is dimensionless;

P _w water phase pressure, Pa;

P _o,m is the matrix oil phase pressure, Pa;

P _o,f is the pressure of the fractured oil phase, Pa;

P _c,m capillary force in the matrix, Pa;

P _c,f capillary force in the crack, Pa;

ρ _o0 oil phase density at initial pressure of the reservoir, kg/m ³ ；

ρ _w0 The density of the water phase at the initial pressure of the reservoir, kg/m ³ ；

ρ _o,m Is the density of the oil phase in the matrix under pressure Po, kg/m ³ ；

ρ _w,m Is the density of the aqueous phase in the matrix under pressure Po, kg/m ³ ；

ρ _f,m Is the density of the aqueous phase in the fracture under pressure Po, kg/m ³ ；

ρ _o,f Is the density of the oil phase in the fracture under pressure Po, kg/m ³ ；

φ _m Is the porosity of the matrix, and has no dimension;

φ _f is crack porosity, dimensionless;

d _f is the seam width, mm;

C _w for the overall compressibility in the aqueous phase, it can be written as Cw ═ C _f +ΦC _l,w ；

C _o Is the overall compression factor of the oil phase and can be written as Co ═ C _f +ΦC _l,o ；

C _l,w The water phase compressibility in the reservoir is 1/Pa;

C _l,o is the compression coefficient of oil phase in reservoir, 1/Pa

C _f The compression coefficient of the rock is 1/Pa;

phi f is the porosity of the crack, and is dimensionless;

Φ _m is the porosity of the matrix, and has no dimension;

μ _w is the fluid viscosity, pas;

μ _o is the fluid viscosity, pas;

S _w,m the water saturation of the matrix is dimensionless;

S _w,f the fracture water saturation is dimensionless;

pi is osmotic pressure, Pa;

c is the salt concentration, mol/m ³ ；

C _m Is the salt ion concentration in the matrix, mol/m ³ ；

C _f Is the salt ion concentration in the fracture, mol/m ³ ；

V is the molar volume of water, taken at 1.8X 10 ^-5 m ³ /mol；

R is gas constant, and is 8.31X 10 ³ Pa·L(mol·K) ^-1 ；

T is the temperature, K;

a _Ⅰ represents the activity of the low salinity water without dimension;

a _Ⅱ represents the activity of the high salinity water without dimension;

v is the number of ions making up the solution, dimensionless;

F _diff flux generated for diffusion, mol/(m) ² ·s)；

E _op The efficiency of a semipermeable membrane is dimensionless, and the efficiency of an ideal semipermeable membrane is 1, i.e. no substance is allowed to pass through;

d is the diffusion coefficient, m ² /s；

F _adv Flux generated for convection, mol/(m) ² ·s)；

With S _w 、P _o Establishing a dependent variable through a coefficient differential equation interface, writing (1) and (2) into corresponding weak form equations, adding through weak contribution boundaries, and regarding the fracture as an inner boundary of the oil reservoir for simulation. The injection boundary conditions are just embodied in the weak form equation and are therefore also added by the weak contribution boundary conditions.

Coupling of stress field and seepage fieldMainly reflects the permeability change caused by pore deformation, the change rule of the porosity along with the stress is,

while the permeability and the porosity are satisfied,

k＝k ₀ (φ/φ ₀ ) ³ (15)

in the formula ₀ The porosity is in a stress-free state and is dimensionless; phi is a _r Residual porosity, dimensionless; k is a radical of ₀ Permeability in the unstressed state, m ² . Wherein

α _φ ＝5×10 ^-8 Pa ^-1 。

Step S22: and establishing a salt ion migration control equation.

Salt ion transport involves both convection and diffusion. Salt ion diffusion is related to concentration gradients, generally describing the constitutive equation for diffusion:

in the formula, F _diff Flux generated for diffusion, mol/(m) ² ·s)；E _op The efficiency of a semipermeable membrane is dimensionless, and the efficiency of an ideal semipermeable membrane is 1, i.e. no substance is allowed to pass through; d is the diffusion coefficient, m ² /s。

The ion transport due to the flow is:

F _adv ＝Cu (17)

in the formula, F _adv For diffusion flux, mol/(m) ² S); u is the velocity of the fluid flow, m/s.

Conservation equation of salt ion in matrix:

conservation equation of salt ion in the fracture:

the above equation is described by a stable convection-diffusion equation interface, and the expression of the corresponding term is input in the equation.

Step S23: establishing a stress field control equation, wherein the stress field control is mainly controlled by a balance equation, a constitutive equation and a displacement equation, and the constitutive equation of the isotropic linear thermoelastic material is as follows:

σ′ _ij ＝2Gε _ij +λε _kk δ _ij -K′α _T Tδ _ij -α _B pδ _ij (20)

in the formula, σ _ij Is a stress tensor, representing the stress state of a point, Pa; epsilon _ij Is a strain tensor, dimensionless; u. of _i Is the displacement tensor, m; λ is Lame constant, Pa,; g is shear modulus, Pa,; k' is the bulk modulus of elasticity, Pa,; alpha is alpha _B Biot coefficient, dimensionless; (ii) a F is the volume force, N/m, p' is the average pore pressure, 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 (21)

the Skempton corrected effective stress calculation formula proposed by Terzaghi is:

σ′＝σ-α _B p (22)

the quasi-static equilibrium differential equation is:

σ′ _ij，j +F _i ＝0 (23)

the relationship between strain and displacement is:

in the formula of _ij Is the strain tensor; u. of _i Is the displacement tensor, m; f is the volume force, Pa;

the stress field equation uses a solid mechanics interface and corresponding parameters are input;

step S24: establishing a temperature field control equation, and obtaining the temperature field control equation according to the law of conservation of energy as follows:

wherein (ρ C) _p ) _eff Is the effective specific heat capacity, lambda, 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 Respectively the specific heat capacity, lambda, of the rock skeleton, oil and water _s ，λ _o ，λ _w The thermal conductivity coefficients of the rock skeleton, oil and water are respectively.

For isotropic materials, pure shear deformation does not generate heat effect, and only volume change is carried out

A coupling effect between the strain field and the temperature field is caused. Compression

Release heat and expand

Absorbing heat. Corresponding to the addition amount of a heat source,

the above equation is described by a porous medium heat transfer interface, parameters such as effective specific heat capacity, effective thermal conductivity and effective thermal convection coefficient are defined firstly, and then corresponding parameters are input into the porous medium heat transfer interface for simulation.

Step S25: the solution was performed using COMSOL Multiphysic software.

The COMSOL Multiphysic software is used and is based on a finite element method, a full-coupling solving method is provided for multiple physical fields, four solidified fields of heat flow can be combined to form a uniform coupling equation set for solving in coupling analysis, and dependent variables of the independent fields are obtained through calculation; the method is characterized in that a PDE module is used for solving the problem of oil-water two-phase flow, a porous medium heat transfer module is used for solving the problem of temperature field, solid mechanics is used for solving the problem of stress field, COMSOL Multiphysic combines the three to form a differential equation set expressed by a general formula, and THMC four-field full-coupling solution is realized.

Step S3 includes: step S31: performing reservoir development dynamic parameter prediction simulation of the fracturing process according to the mathematical model and the boundary conditions; step S32: and outputting numerical simulation results including reservoir temperature, pressure, water saturation and salt concentration.

Step S31 includes: assuming that the initial time is 0, injection time is t1, shut-in time is t2, and flowback time is t 3;

(1) initial conditions

The initial conditions of temperature and pressure in the fracturing process are the initial formation temperature and formation pressure of a reservoir, and are realized through the Dielder boundary condition, namely,

T＝T ₀ (t＝0) (28)

p＝p ₀ (t＝0) (29)

in the formula: t is a unit of ₀ Initial formation temperature, K; p is a radical of ₀ Initial formation pressure, Pa;

for the stress field, the initial state is in a stressed and unstrained state due to the reservoir being constrained by the geostress, which is achieved by adding boundary loads, i.e.,

σ _ij (x,y,z,t＝0)＝[σ _v σ _H σ _h ] (30)

in the formula, σ _v Is vertical ground stress, Pa; sigma _H In order to level out the maximum ground stress,Pa；σ _h is the horizontal minimum ground stress, Pa;

(2) boundary condition

The boundary condition within the fracturing process is the flow boundary condition, which is mentioned in weak form, i.e.,

wherein A is the cross-sectional area, m ² ；q _w For the flow rate of the injected fracturing fluid, kg.m ³ /s，p _w Is bottom hole flowing pressure; the boundary conditions of the temperature field and the seepage field are similar,

wherein q is a heat source term due to an injection well, W/m ² ；C _pw The specific heat capacity of the injected water; t is _inj Is the injection temperature, K.

Step S4 includes: step S41: performing reservoir development dynamic parameter prediction simulation of the well shut-in process according to the mathematical model and the boundary conditions; step S42: and outputting numerical simulation results including reservoir temperature, pressure, water saturation and salt concentration.

Step S41 includes: (1) initial conditions

The initial conditions of temperature and pressure in the shut-in process are the reservoir temperature and pressure at the end of fracturing, which are realized by the Dielder boundary conditions, that is,

T＝T ₁ (t＝t ₁ ) (33)

p＝p ₁ (t＝t ₁ ) (34)

in the formula: t is ₁ Is the reservoir temperature at the end of fracturing, K; p is a radical of ₁ Reservoir pressure at the end of fracturing, Pa;

for the stress field, the initial state is in a stressed and unstrained state due to the reservoir being constrained by the geostress, which is achieved by adding boundary loads, i.e.,

σ _ij (x,y,z,t＝0)＝[σ _v σ _H σ _h ] (35)

in the formula, σ _v Is vertical ground stress, Pa; sigma _H Pa, horizontal maximum ground stress; sigma _h Is the horizontal minimum stress, Pa.

Step S5 includes: step S51: performing reservoir development dynamic parameter prediction simulation of the fracturing flow-back process according to the mathematical model and the boundary conditions; step S52: and outputting numerical simulation results including reservoir temperature, pressure, water saturation and salt concentration.

Step S51 includes: assuming that the initial time is 0 time, the injection time is t1 time, the shut-in time is t2 time, and the flow-back time is t3 time; (1) initial conditions

The initial conditions of the temperature and the pressure in the fracturing flow-back process are the reservoir temperature and the pressure at the end of well shut-in, and are realized through the Dielder boundary condition, namely,

T＝T ₂ (t＝t ₂ ) (36)

p＝p ₂ (t＝t ₂ ) (37)

in the formula: t is a unit of ₂ The reservoir temperature after the shut-in is finished, K; p is a radical of ₂ Reservoir pressure, Pa, after shut-in is finished;

for the stress field, the initial state is in a stressed and unstrained state due to the reservoir being constrained by the geostress, which is achieved by adding boundary loads, i.e.,

σ _ij (x,y,z,t＝0)＝[σ _v σ _H σ _h ] (38)

in the formula, σ _v Is vertical ground stress, Pa; sigma _H Pa for horizontal maximum ground stress; sigma _h Is the horizontal minimum stress, Pa;

(2) boundary condition

The boundary condition in the flow-back process is a pressure boundary condition, the pressure boundary condition is a constant pressure boundary, that is,

wherein A is the cross-sectional area, m ² ；q _w For the flow rate of the injected fracturing fluid, kg.m ³ /s，p _w Is the bottom hole flowing pressure.

According to the method, the osmotic pressure and heat flow solidification coupling effect are comprehensively considered, the fracturing and development dynamics of the shale oil reservoir with high clay content are predicted, the main influence factors of the fracturing and development dynamics are determined, and guidance is provided for shale oil reservoir fracturing optimization.

EXAMPLE III

Taking an X shale oil well as an example, specific reservoir geology, engineering parameters and fluid parameters are shown in a table 1, a horizontal well is subjected to multi-cluster close cutting staged fracturing, one cluster of fractures is simulated, the process of pumping for 100min, closing the well for 30 days and flowback (production) for 5 years is simulated,

TABLE 1 model parameters

The capillary force is generally obtained through direct measurement of experiments or through combination of experiments and empirical formulas, wherein the empirical formulas are shown in the following, and the capillary force can be obtained through measurement of interfacial tension, porosity and permeability of rocks.

Under different pressure and temperature, the density, viscosity and specific heat conductivity coefficient of oil and water are different, the related data can be obtained by researching related documents, the formula is fitted by experiments and related software, an empirical formula is given here for reference,

μ _w ＝1.3799-0.0212T+1.3604×10 ^-4 T ² -4.6454×10 ^-7 T ³ +8.9043×10 ^-10 T ⁴ -9.0791×10 ^-13 T ⁵ +3.8457×10 ^-16 T ⁶ ,T∈[273,413]

μ _w ＝0.0040-2.1075×10 ^-5 T+3.8577×10 ^-8 T ² -2.3973×10 ^-11 T ³ ,T∈[413,553]

C _pw ＝12010-80.4×T+0.3×T ² -5.4×10 ^-4 T ³ +3.6×10 ^-7 T ⁴ ,T∈[273,553]

ρ _w ＝838.4661+1.4005T-0.0030T ² +3.7182×10 ^-7 T ³ ,T∈[273,553]

λ _w ＝-0.8691+0.0089T-1.5837×10 ^-5 T ² +7.9754×10 ^-9 T ³ ,T∈[273,553]

μ _w ＝91.45245-1.33227T+0.00778T ² -2.27278×10 ^-5 T ³ +3.32420×10 ^-8 T ⁴ -1.94631×10 ^-11 T ⁵ ,T∈[273,453]

C _pw ＝-13408.1491+123.04415T-0.33540T ² +3.125×10 ^-4 T ³ ,T∈[273,453]

ρ _w ＝1055.04607-0.58175T-6.40532×10 ^-5 T ² ,T∈[273,453]

λ _w ＝0.134299-8.049738×10 ^-5 T,T∈[273,453]

based on the step 1, according to the parameters and the geological model in the table 1, an X shale oil well single cluster fracturing physical model is established.

Based on the step 2, oil-water two-phase flow is realized through a coefficient differential equation of a PDE (partial differential equation) module, the crack is regarded as the inner boundary of the whole domain, and a weak form of a corresponding control equation is added in the crack region to represent the crack equation before solving, so that the oil-water two-phase flow containing discrete cracks can be realized; solving a salt ion transport equation through a convection-diffusion equation of a PDE (partial differential equation) module; and solving the temperature field through the porous medium heat transfer module. And solving the stress field and the displacement field through a solid mechanics module. The coupling relation is realized by parameter variable transmission and adding source.

Based on step 3, as shown in fig. 4-7, the initial conditions in the fracturing process are the initial pressure, temperature and salt concentration of the reservoir, the seepage field uses the flow boundary condition, the temperature field uses the heat flow boundary condition, the stress field uses the boundary load to simulate the constraint of the ground stress on the reservoir, and the concentration field uses the dirichlet boundary condition. The results of the solution are shown in fig. 4(a), 5(a), 6(a), and 7 (a).

Based on step 4, as shown in fig. 4-7, the initial conditions during shut-in were reservoir temperature, pressure and salt concentration at the end of fracturing, with the end of the first study as the initial time of the study, by adding a study step, a simulation was run. The results of the solution are shown in fig. 4(b), 5(b), 6(b), and 7 (b).

Based on step 5, as shown in fig. 4-7, the initial conditions during flowback were reservoir temperature, pressure and salt concentration at the end of shut-in, and simulations were run with the end of the second study as the initial time of the study by adding the study steps. The seepage field uses a pressure boundary condition, and the stress field uses boundary load to simulate the restraint of the ground stress on the reservoir. The results of the solution are shown in fig. 4(c), 5(c), 6(c), and 7 (c).

The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention, and various modifications and changes may be made by those skilled in the art. Any modification, equivalent replacement, or improvement made within the spirit and principle of the present invention should be included in the protection scope of the present invention.

Claims (10)

1. A numerical simulation method for a high-clay-content shale oil reservoir fracturing horizontal well is characterized by comprising the following steps of:

step S1, obtaining reservoir physical property parameters according to the actual shale oil reservoir characteristics and the rock core, well logging and fluid test data; acquiring shale oil reservoir fracture geometric information, performing dimensionality reduction treatment on fractures in an oil reservoir as an inner boundary of the oil reservoir, establishing an oil reservoir geometric model, and then performing geometric subdivision on the oil reservoir geometric model by adopting triangular meshes to form discrete units;

step S2, establishing a shale oil reservoir fracturing well reservoir development dynamic parameter prediction model based on the discrete fracture model;

step S3, reservoir development dynamic parameter prediction simulation of the fracturing process is carried out;

step S4, performing reservoir development dynamic parameter prediction simulation in the well closing process;

and step S5, performing reservoir development dynamic parameter prediction simulation in the fracturing flow-back process.

2. The numerical simulation method for horizontal fractured wells of shale oil reservoirs with high clay content according to claim 1, wherein the step S1 comprises the following steps:

collected reservoir geological parameters including maximum and minimum horizontal principal stress distribution, rock Young's modulus, Poisson's ratio, formation temperature and pressure, porosity, permeability, semi-permeable membrane efficiency; collecting reservoir fluid parameters including viscosity, mineralization degree, heat conductivity coefficient, specific heat capacity, phase permeation curve and capillary force; collected completion information for fractured wells, including cluster spacing; collecting fracturing construction parameters including construction discharge capacity, viscosity of fracturing fluid, mineralization degree, heat conductivity coefficient and specific heat capacity; collecting the distribution of hydraulic fractures and natural fractures, and the length, width and seepage of fractures; acquiring the geometrical information of the shale oil reservoir fracture according to the actual geological data or the existing geological model data of the reservoir: including crack position, size, density, strike, opening.

3. The numerical simulation method for horizontal fractured wells of shale oil reservoirs with high clay content according to claim 1, wherein the step S2 comprises the following steps:

step S21: respectively establishing an oil-water two-phase control equation of a matrix and a crack:

the matrix water phase control equation is:

the matrix oil phase equation is:

the fracture water phase equation is:

the fracture oil phase equation is:

step S22: establishing a salt ion transport control equation;

step S23: establishing a stress field control equation, wherein the stress field control is mainly controlled by a balance equation, a constitutive equation and a displacement equation, and the constitutive equation of the isotropic linear thermoelastic material is as follows:

σ′ _ij ＝2Gε _ij +λε _kk δ _ij -K′α _T Tδ _ij -α _B pδ _ij

in the formula, σ _ij Is a stress tensor, representing the stress state of a point, Pa; epsilon _ij Is a strain tensor, dimensionless; u. of _i Is the displacement tensor, m; λ is Lame constant, Pa,; g is shear modulus, Pa,; k' is the bulk modulus of elasticity, Pa,; alpha is alpha _B Biot coefficient, dimensionless; (ii) a F is the volume force, N/m, p' is the average pore pressure, 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

the Skempton corrected effective stress calculation formula proposed by Terzaghi is:

σ′＝σ-α _B p

the quasi-static equilibrium differential equation is:

σ′ _ij,j +F _i ＝0

the relationship between strain and displacement is:

in the formula of _ij Is the strain tensor; u. of _i Is the displacement tensor, m; f is volume force, Pa;

the stress field equation uses a solid mechanics interface and corresponding parameters are input;

step S24: establishing a temperature field control equation, and obtaining the temperature field control equation according to the law of conservation of energy as follows:

wherein (ρ C) _p ) _eff Is the effective specific heat capacity, lambda, 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 Respectively the specific heat capacity, lambda, of the rock skeleton, oil and water _s ，λ _o ，λ _w The thermal conductivity coefficients of the rock skeleton, oil and water are respectively;

step S25: the solution was performed using COMSOL Multiphysic software.

4. The numerical simulation method for horizontal fractured wells of high clay content shale oil reservoirs of claim 3, wherein the step S25 comprises:

the COMSOL Multiphysic software is used and is based on a finite element method, a full-coupling solving method is provided for multiple physical fields, four solidified fields of heat flow can be combined to form a uniform coupling equation set for solving in coupling analysis, and dependent variables of the independent fields are obtained through calculation; the method is characterized in that a PDE module is used for solving the problem of oil-water two-phase flow, a porous medium heat transfer module is used for solving the problem of temperature field, solid mechanics is used for solving the problem of stress field, COMSOL Multiphysic combines the three to form a differential equation set expressed by a general formula, and THMC four-field full-coupling solution is realized.

5. The numerical simulation method for horizontal fracturing wells of high clay content shale oil reservoirs according to any one of claims 1-4, wherein the step S3 comprises:

step S31: performing reservoir development dynamic parameter prediction simulation of the fracturing process according to the mathematical model and the boundary conditions;

step S32: and outputting numerical simulation results including reservoir temperature, pressure, water saturation and salt concentration.

6. The numerical simulation method for horizontal fractured wells of shale oil reservoirs with high clay content according to claim 5, wherein the step S31 comprises the following steps:

assuming that the initial time is 0 time, the injection time is t1 time, the shut-in time is t2 time, and the flow-back time is t3 time;

(1) initial conditions

The initial conditions of temperature and pressure in the fracturing process are the initial formation temperature and formation pressure of a reservoir, and are realized through the Dielder boundary condition, namely,

T＝T ₀ (t＝0)

p＝p ₀ (t＝0)

in the formula: t is ₀ Initial formation temperature, K; p is a radical of ₀ Initial formation pressure, Pa;

for the stress field, the initial state is in a stressed and unstrained state due to the reservoir being constrained by the geostress, which is achieved by adding boundary loads, i.e.,

σ _ij (x,y,z,t＝0)＝[σ _v σ _H σ _h ]

in the formula, σ _v Is vertical ground stress, Pa; sigma _H Pa for horizontal maximum ground stress; sigma _h Is the horizontal minimum stress, Pa;

(2) boundary condition

The boundary condition within the fracturing process is the flow boundary condition, which is mentioned in weak form, i.e.,

wherein A is the cross-sectional area, m ² ；q _w For the flow rate of the injected fracturing fluid, kg.m ³ /s，p _w Is bottom hole flowing pressure;

the boundary conditions of the temperature field and the seepage field are similar,

wherein q is a heat source term due to an injection well, W/m ² ；C _pw The specific heat capacity of the injected water; t is a unit of _inj Is the injection temperature, K.

7. The numerical simulation method for horizontal fracturing wells of high clay content shale oil reservoirs according to any one of claims 1-4, wherein the step S4 comprises:

step S41: performing reservoir development dynamic parameter prediction simulation of the well shut-in process according to the mathematical model and the boundary conditions;

step S42: and outputting numerical simulation results including reservoir temperature, pressure, water saturation and salt concentration.

8. The numerical simulation method for horizontal fractured wells of high clay content shale oil reservoirs of claim 7, wherein the step S41 comprises:

(1) initial conditions

The initial conditions of temperature and pressure in the shut-in process are the reservoir temperature and pressure at the end of fracturing, which are realized by the Dielder boundary conditions, that is,

T＝T ₁ (t＝t ₁ )

p＝p ₁ (t＝t ₁ )

in the formula: t is ₁ Is the reservoir temperature at the end of the fracturing,K；p ₁ reservoir pressure at the end of fracturing, Pa;

for the stress field, the initial state is in a stressed and unstrained state due to the reservoir being constrained by the geostress, which is achieved by adding boundary loads, i.e.,

σ _ij (x,y,z,t＝0)＝[σ _v σ _H σ _h ]

in the formula, σ _v Is vertical ground stress, Pa; sigma _H Pa for horizontal maximum ground stress; sigma _h Is the horizontal minimum stress, Pa.

9. The numerical simulation method for horizontal fracturing wells of high clay content shale oil reservoirs according to any one of claims 1-4, wherein the step S5 comprises:

step S51: performing reservoir development dynamic parameter prediction simulation of the fracturing flow-back process according to the mathematical model and the boundary conditions;

step S52: and outputting numerical simulation results including reservoir temperature, pressure, water saturation and salt concentration.

10. The numerical simulation method for horizontal fractured wells of high clay content shale oil reservoirs of claim 9, wherein the step S51 comprises:

assuming that the initial time is 0 time, the injection time is t1 time, the shut-in time is t2 time, and the flow-back time is t3 time;

(1) initial conditions

The initial conditions of the temperature and the pressure in the fracturing flow-back process are the reservoir temperature and the pressure at the end of well shut-in, and are realized through the Dielder boundary condition, namely,

T＝T ₂ (t＝t ₂ )

p＝p ₂ (t＝t ₂ )

in the formula: t is ₂ The reservoir temperature after the shut-in is finished, K; p is a radical of ₂ Reservoir pressure, Pa, after shut-in is finished;

for the stress field, the initial state is in a stressed and unstrained state due to the reservoir being constrained by the geostress, which is achieved by adding boundary loads, i.e.,

σ _ij (x,y,z,t＝0)＝[σ _v σ _H σ _h ]

in the formula, σ _v Is vertical ground stress, Pa; sigma _H Pa for horizontal maximum ground stress; sigma _h Is the horizontal minimum ground stress, Pa;

(2) boundary condition

The boundary condition in the flow-back process is a pressure boundary condition, the pressure boundary condition is a constant pressure boundary, that is,

wherein A is the cross-sectional area, m ² ；q _w For the flow rate of the injected fracturing fluid, kg.m ³ /s，p _w Is the bottom hole flowing pressure.

Images