Patent Grant
11686190
The present invention relates to the technical field of oilfield exploration, in particular to a 4D quantitative and intelligent diagnosis method and system for spatio-temporal evolution of oil-gas reservoir damage types and extent.
In various stages of oilfield exploration and development, due to the influence of a variety of internal and external factors, the original physical, chemical, thermodynamic and hydrodynamic equilibrium of a reservoir change, which inevitably causes a reduction in reservoir internal permeability in regions close to a well wall and even regions far from the well wall of the reservoir to clog flowing of fluid so as to cause reservoir damage and a well production decline, or even “kill” the reservoir, such that an oil-gas well produces no oil or gas. The causes of reservoir damage are diverse and complex. Especially during production, a storage-permeability space, surface wettability, a hydrodynamic field, a temperature field, rock types, etc. of reservoir rock continuously change, such that the damage mechanism changes over time, and the damage lasts for a long period and covers a wide range, and is more complex and superimposed. Once reservoir damage occurs, corresponding declogging measures must be taken according to reservoir damage to restore fluid flow channels, so as to increase oil-gas well production and water well injection capacity.
Usually, there are various causes and types of oil-gas reservoir damage, and multiple types of damage coexist, each type causing a different extent of oil-gas reservoir damage. It is very difficult but very important to determine which types (several or more than a dozen types) on earth lead to the oil-gas reservoir damage, the damage extent caused by each damage type or a contribution rate of each damage type to a total damage extent, a damaged zone radius and a total damaged zone radius due to each damage type, a spatio-temporal evolution law of oil-gas reservoir damage by each damage type, etc. It is indispensable not only for avoiding the occurrence of oil-gas reservoir damage for wells in normal production to increase production (or injection capacity), taking unclogging measures for wells that have stopped production to restore production, and taking protective measures for wells facing stop of production to extend the production life, but also for improving the precision of numerical simulation of oil pools and accurately predicting “sweet spots”, etc. It is of great significance not only for strong theoretical significance, but also for a high practical value.
Currently, methods for diagnosing reservoir damage can be divided into field diagnostic methods and indoor evaluation methods. The field diagnostic methods include a well testing method. Although the well testing method can quantitatively provide important parameters such as skin factor, clogging ratio, and additional pressure drop that characterize a damage extent of a reservoir in a preset region of a well to be diagnosed, the skin factor characterized thereby is associated with other parameters. In other words, the skin factor obtained by the well testing method not only reflects a true reservoir damage characteristic, but is also an overall manifestation of various aspects and factors (i.e., it is the sum of a true damage skin factor and a pseudo-skin factor composed of a well deviation skin factor, a reservoir shape skin factor, a partial penetration skin factor, a non-Darcy flow skin factor, a perforation skin factor and the like), and the skin factor must be decomposed to obtain the true damage skin factor. The indoor evaluation method includes a core flow test method. The core flow test method uses a permeability change before and after core displacement to get a damage extent, and is suitable for studying reservoir damage by single-factor, but it is difficult to reflect reservoir damage laws on a larger scale. Moreover, as indoor core experiment conditions are idealized, cores used for evaluation are in an original state, and dynamic changes of reservoir properties cannot be taken into account, experimental results differ greatly from the real damage of the downhole reservoir. That is to say, so far, researches in China and other countries for decades have not achieved accurate quantitative and rapid diagnosis on reservoir damage causes, types and extent, not to mention spatio-temporal evolution diagnosis, which has become a key and major difficult problem at core of the strategy of “increasing reserves and production” in the international arena.
An objective of the present invention is to provided a 4D quantitative and intelligent diagnosis method and system for spatio-temporal evolution of oil-gas reservoir damage types and extent, which can quantitatively simulate characteristic parameters of reservoir damage caused by relevant factors and a total characteristic parameter of the reservoir damage, so that the oil-gas reservoir damage of an oil-gas production well or an injection well in various stages or aspects of exploration and development is “transparent, digitized, dynamic, visualized and intelligentized” in the spatial and temporal domains, and an experience of being personally on the scene is achieved, thereby, for a well without reservoir damage, performing quantitative prediction of reservoir damage and spatio-temporal deduction of damage laws, which is of scientific guidance significance for preventing or avoiding reservoir damage, and formulating development plans for oil pools and subsequent well stimulation measures, and for a well with reservoir damage, also performing quantitative diagnosis of reservoir damage and spatio-temporal deduction of damage laws, which is of very great significance for achieving optimal design of a declogging measure and improvement or restoration of oil-gas well production and water well injection capacity. Therefore, this makes the present invention be an indispensable technology in exploration and development of oil-gas pools.
To achieve the above objective, in a first aspect, the present invention provides a method for determining a damage extent of a reservoir, including: based on a spatio-temporal evolution simulation equation of reservoir damage by each of a plurality of factors, determining a characteristic parameter characterizing reservoir damage by each of the plurality of factors, wherein the reservoir is located in a preset region of a well to be diagnosed; and determining an effective characteristic parameter characterizing the damage extent of the reservoir based on the characteristic parameter characterizing reservoir damage by each of the plurality of factors.
By adopting the above technical solution, according to the present invention based on the spatio-temporal evolution simulation equation of reservoir damage by each of the plurality of factors, a characteristic parameter characterizing reservoir damage by each of the plurality of factors are creatively determined; and the effective characteristic parameter characterizing the damage extent of the reservoir is determined based on the characteristic parameter characterizing reservoir damage by each of the plurality of factors. Thus, by using the spatio-temporal evolution simulation equations of reservoir damage by the relevant factors, the characteristic parameters (such as permeability) of reservoir damage caused by the factors respectively and a total characteristic parameter (such as total permeability or effective permeability) of reservoir damage caused by the plurality of relevant factors can be quantitatively simulated. Therefore performing quantitative prediction of reservoir damage and spatio-temporal deduction of damage laws of scientific guidance significance for preventing or avoiding reservoir damage, and formulating development plans for oil pools and subsequent well stimulation measures for a well without reservoir damage, and is of very great significance for optimal design of a declogging measure and improvement or restoration of oil well production and water well injection capacity for damaged wells, and improvement of numerical simulation precision of oil pools.
In a second aspect, the present invention further provides a system for determining a damage extent of a reservoir, including: a first parameter determination device configured to, based on a spatio-temporal evolution simulation equation of reservoir damage by each of a plurality of factors, determine a characteristic parameter characterizing reservoir damage by each of the plurality of factors, wherein the reservoir is located in a preset region of a well to be diagnosed; and a second parameter determination device configured to determine an effective characteristic parameter characterizing the damage extent of the reservoir based on the characteristic parameter characterizing reservoir damage by each of the plurality of factors.
The system for determining the damage extent of the reservoir has the same advantages as the above method for determining the damage extent of the reservoir with respect to the prior art, which will not be described in detail here.
Correspondingly, in a third aspect, the present invention further provides a machine-readable storage medium that stores instructions which are configured to cause a machine to execute the method for determining a damage extent of a reservoir.
Other features and advantages of embodiments of the present invention will be described in detail in the subsequent section of detailed description of the embodiments.
The accompanying drawings are used to provide further understanding of the embodiments of the present invention and form part of the specification, and are used, together with the following detailed description of the embodiments, for explaining the embodiments of the present invention, but do not limit the embodiments of the present invention. In the drawings:
The specific implementations of the present invention will be described in detail below with reference to the accompanying drawings. It should be understood that the specific implementations described here are only used for illustrating and explaining the present invention, instead of limiting the present invention.
During development of an oilfield, different types of reservoir damage are produced by various reasons, resulting in lower production or even shutdown of oil wells and reduction of water injection capacity of water wells to cause huge economic losses to the oilfield. Implementing declogging on a damaged reservoir is an important measure to increase the production and improve the recovery ratio. However, an optimized decision of the declogging measure must be based on quantitative diagnosis of a reservoir damage type and extent. Due to a long period and a wide range of the reservoir damage during oilfield development, and more complexity and superimposition of the damage, the diagnosis of the reservoir damage is more difficult, so far, there is a lack of diagnosis methods and technology with high accuracy and practicality at home and abroad, so it is impossible to quantitatively diagnose the proportion of each damage type in a total damage extent of a well to be diagnosed, and even more impossible to achieve quantitative simulation of spatio-temporal evolution.
For the complex mechanism and relevant parameters that influence characteristic parameters such as permeability of a reservoir in a region around the well to be diagnosed, embodiments (e.g., embodiments 1-13) of the present invention creatively provide spatio-temporal evolution simulation equations that quantitatively simulate reservoir damage by relevant factors (factor 1—extraneous solid-phase particles, factor 2—clay swelling, factor 3—inorganic precipitation, factor 4—the migration of fine particle within a reservoir, factor 5—water lock effect, factor 6—stress sensitivity, factor 7—sand production, factor 8—wettability reversal, factor 9—emulsification, factor 10—organic scale, factor 11—Jamin effect, factor 12—bacteria, and factor 13—polymer), then the parameters such as the permeability of the reservoir damaged by the corresponding factors can be determined quantitatively by using the various spatio-temporal evolution simulation equations, and finally the total damage extent of the well to be diagnosed in different time and space and contributions of each relevant factor to the total damage extent can be determined in conjunction with the principle of superposition.
Specifically, for different types of wells to be diagnosed, factors that influence the permeability and other characteristic parameters of reservoirs in regions around the different types of wells to be diagnosed are obviously different, as shown in Table 1; and for the same type of wells to be diagnosed, relevant factors that influence the permeability and other characteristic parameters of reservoirs in regions around the same type of wells to be diagnosed in different development stages (e.g., a drilling stage, a water injection/a polymer injection stage, and an oil production stage) are also obviously different, as shown in Table 1. That is, the plurality of factors are different for different types of wells to be diagnosed; and the plurality of factors are also different for the same type of wells to be diagnosed in different development stages.
Table 1 relevant factors influencing the permeability of a reservoir in a region around a well to be diagnosed
For the step S102, the determining an effective characteristic parameter characterizing a damage extent of the reservoir may include: determining an effective characteristic parameter F({right arrow over (r)}, t) characterizing the damage extent of the reservoir based on a characteristic parameter Fi({right arrow over (r)}, t) characterizing reservoir damage by an ith factor in the plurality of factors and the following formula,
where Li is a weight of Fi({right arrow over (r)}, t) (which is reasonably set according to the actual situation; for example, it may be set to 1); {right arrow over (r)} is any location in the reservoir; t is any time; and n is the number of the plurality of factors.
Furthermore, in the case where the characteristic parameter Fi({right arrow over (r)}, t) characterizing reservoir damage by the ith factor in the plurality of factors is obtained, based on the characteristic parameter Fi({right arrow over (r)}, t) characterizing reservoir damage by the ith factor in the plurality of factors and the following formula, the proportion δi({right arrow over (r)}, t) or Fi({right arrow over (r)}, t) in the the effective characteristic parameter characterizing the damage extent of the reservoir may also be determined,
wherein the characteristic parameter may be permeability, skin factor and/or permeability damage rate, etc, which may be reasonably set according to actual needs in practical application, and some or all of the above three characteristic parameters may be calculated.
In the case where the well to be diagnosed is a water injection well, a polymer injection well or an oil production well and is in a drilling stage, the plurality of factors include at least two of: extraneous solid-phase particles, clay swelling, the migration of fine particle within the reservoir (which may be referred to as “fine particle migration”), inorganic precipitation, and water lock effect, which may be reasonable set according to actual needs in practical application.
In an embodiment, in the case where the well to be diagnosed is a water injection well, a polymer injection well or an oil production well and is in a drilling stage, the plurality of factors may include extraneous solid-phase particles, clay swelling, and fine particle migration.
For the step S102, the effective characteristic parameter F({right arrow over (r)}, t) characterizing the damage extent of the reservoir is determined based on characteristic parameters F1({right arrow over (r)}, t), F2({right arrow over (r)}, t) and F4({right arrow over (r)}, t) (subscripts thereof correspond to the factor 1—extraneous solid-phase particles, the factor 2—clay swelling and the factor 4—fine particle migration in the drilling stage in Table 1, and the three characteristic parameters may be determined by spatio-temporal evolution simulation equations shown in the following Embodiments 1, 2 and 4, respectively) characterizing reservoir damage respectively by the extraneous solid-phase particles, the clay swelling, and the fine particle migration, and the following formula:
F({right arrow over (r)},t)=F1({right arrow over (r)},t)+F2({right arrow over (r)},t)+F4({right arrow over (r)},t).
In another embodiment, the plurality of factors may further include at least one of:
inorganic precipitation and water lock effect. For example, the plurality of factors may further include both inorganic precipitation and water lock effect.
For the step S102, the effective characteristic parameter F({right arrow over (r)}, t) characterizing the damage extent of the reservoir is determined based on characteristic parameters F1({right arrow over (r)}, t), F2({right arrow over (r)}, t), F4({right arrow over (r)}, t), F3({right arrow over (r)}, t) and F5({right arrow over (r)}, t) (which may be determined by spatio-temporal evolution simulation equations shown in the following Embodiments 1, 2, 4, 3 and 5 respectively) characterizing reservoir damage respectively by the extraneous solid-phase particles, the clay swelling, the fine particle migration, the inorganic precipitation, and the water lock effect, and the following formula:
F({right arrow over (r)}, t)=Σi=15Fi({right arrow over (r)}, t).
In the case where the well to be diagnosed is an oil production well and is in an oil production stage, the plurality of factors include at least two of: fine particle migration, sand production, emulsification, Jamin effect, stress sensitivity, wettability reversal, and organic scale, which may be reasonably set according to actual needs in practical application.
In an embodiment, in the case where the well to be diagnosed is an oil production well and is in an oil production stage, the plurality of factors may include fine particle migration, sand production, emulsification, and organic scale.
For the step S102, the effective characteristic parameter F({right arrow over (r)}, t) characterizing the damage extent of the reservoir is determined based on characteristic parameters F1({right arrow over (r)}, t), F3({right arrow over (r)}, t), F5({right arrow over (r)}, t) and F6({right arrow over (r)}, t) (subscripts thereof correspond to the factor 1—fine particle migration, the factor 3—sand production, the factor 5—emulsification, and the factor 6—Jamin effect in the oil production stage in Table 1, and the four characteristic parameters may be determined by spatio-temporal evolution simulation equations shown in the following Embodiments 4, 7, 9 and 10 respectively) characterizing reservoir damage respectively by the fine particle migration, the sand production, the emulsification, and the organic scale, and the following formula:
F({right arrow over (r)}, t)=F1({right arrow over (r)}, t)+F3({right arrow over (r)}, t)+F5({right arrow over (r)}, t)+F6({right arrow over (r)}, t).
In another embodiment, the plurality of factors may further include at least one of: stress sensitivity, wettability reversal, and organic scale. For example, the plurality of factors may further include three of stress sensitivity, wettability reversal, and Jamin effect.
For the step S102, the effective characteristic parameter F({right arrow over (r)}, t) characterizing the damage extent of the reservoir is determined based on characteristic parameters F1({right arrow over (r)}, t), F3({right arrow over (r)}, t), F5({right arrow over (r)}, t), F6({right arrow over (r)}, t), F2({right arrow over (r)}, t), F4({right arrow over (r)}, t) and F7({right arrow over (r)}, t) (which may be determined by spatio-temporal evolution simulation equations shown in the following Embodiments 4, 7, 9, 10, 6, 8 and 11 respectively) characterizing reservoir damage respectively by the fine particle migration, the sand production, the emulsification, the Jamin effect, the stress sensitivity, the wettability reversal, and the organic scale, and the following formula:
F({right arrow over (r)}, t)=Σi=17Fi({right arrow over (r)}, t).
In the case where the well to be diagnosed is a water injection well and is in a water injection stage, the plurality of factors include at least two of: clay swelling, bacteria, water lock effect, extraneous solid-phase particles, fine particle migration, and inorganic precipitation, which may be reasonably set according to actual needs in practical application.
In an embodiment, in the case where the well to be diagnosed is a water injection well and is in a water injection stage, the plurality of factors may include clay swelling, bacteria, and water lock effect.
For the step S102, the effective characteristic parameter F({right arrow over (r)}, t) characterizing the damage extent of the reservoir is determined based on characteristic parameters F2({right arrow over (r)}, t), F5({right arrow over (r)}, t) and F6({right arrow over (r)}, t) (subscripts thereof correspond to the factor 2—clay swelling, the factor 5—water lock effect, and the factor 6—bacteria in the water injection stage in Table 1, and the three characteristic parameters may be determined by spatio-temporal evolution simulation equations shown in the following Embodiments 2, 12 and 5, respectively) characterizing reservoir damage respectively by the clay swelling, the bacteria, and the water lock effect, and the following formula:
F({right arrow over (r)}, t)=F2({right arrow over (r)}, t)+F5({right arrow over (r)}, t)+F6({right arrow over (r)}, t).
In another embodiment, the plurality of factors may further include at least one of: extraneous solid-phase particles, fine particle migration, and inorganic precipitation. For example, the plurality of factors may further include three of extraneous solid-phase particles, fine particle migration, and inorganic precipitation.
For the step S102, the effective characteristic parameter F({right arrow over (r)}, t) characterizing the damage extent of the reservoir is determined based on characteristic parameters F2({right arrow over (r)}, t), F5({right arrow over (r)}, t), F6({right arrow over (r)}, t), F1({right arrow over (r)}, t), F4({right arrow over (r)}, t) and F3({right arrow over (r)}, t) (which may be determined by spatio-temporal evolution simulation equations shown in the following Embodiments 2, 12, 5, 1, 4 and 3 respectively) characterizing reservoir damage respectively by the clay swelling, the bacteria, the water lock effect, the extraneous solid-phase particles, the fine particle migration, and the inorganic precipitation, and the following formula:
F({right arrow over (r)}, t)=Σi=16Fi({right arrow over (r)}, t).
In the case where the well to be diagnosed is a polymer injection well and is in a polymer injection stage, the plurality of factors include at least two of: polymer, clay swelling, extraneous solid-phase particles, fine particle migration, and inorganic precipitation, which may be reasonably set according to actual needs in practical application.
In an embodiment, in the case where the well to be diagnosed is a polymer injection well and is in a polymer injection stage, the plurality of factors may include polymer and clay swelling.
For the step S102, the effective characteristic parameter F({right arrow over (r)}, t) characterizing the damage extent of the reservoir is determined based on characteristic parameters F2({right arrow over (r)}, t) and F5({right arrow over (r)}, t) (subscripts thereof correspond to the factor 2—clay swelling, and the factor 5—polymer in the polymer injection stage in Table 1, and the two characteristic parameters may be determined by spatio-temporal evolution simulation equations shown in the following Embodiments 13 and 2, respectively) characterizing reservoir damage respectively by the polymer, and the clay swelling, and the following formula:
F({right arrow over (r)}, t)=F2({right arrow over (r)}, t)+F5({right arrow over (r)}, t).
In another embodiment, the plurality of factors may further include at least one of:
extraneous solid-phase particles, fine particle migration, and inorganic precipitation. For example, the plurality of factors may further include three of extraneous solid-phase particles, fine particle migration, and inorganic precipitation.
For the step S102, the effective characteristic parameter F({right arrow over (r)}, t) characterizing the damage extent of the reservoir is determined based on characteristic parameters F2({right arrow over (r)}, t), F5({right arrow over (r)}, t), F1({right arrow over (r)}, t), F4({right arrow over (r)}, t), and F3({right arrow over (r)}, t) (which may be determined by spatio-temporal evolution simulation equations shown in the following Embodiments 2, 13, 1, 4 and 3 respectively) characterizing reservoir damage respectively by the polymer, the clay swelling, the extraneous solid-phase particles, the fine particle migration, and the inorganic precipitation, and the following formula:
F({right arrow over (r)}, t)=Σi=15Fi({right arrow over (r)}, t).
In practical application, the plurality of factors may be reasonably set according to actual needs, i.e., not limited by the combinations of above embodiments.
Processes of establishing the spatio-temporal evolution simulation equation of reservoir damage by each of the 13 relevant factors in Table 1 described above and determining the related characteristic parameters by using the spatio-temporal evolution simulation equations is described below respectively (see Embodiments 1-13 for details, the 13 embodiments respectively corresponding to the 13 factors in Table 1). For the same factor involved in different stages, parameters such as permeability of the reservoir damaged by the same factor may be determined by referring to corresponding embodiments; specific values of some parameters in the corresponding embodiments described in different stages may be different.
It is to be noted that for the sake of simple description, for physical and chemical quantities that evolve with time and space in embodiments of the present invention, the variable ({right arrow over (r)}, t) may be omitted, for example, K({right arrow over (r)}, t) may be shortened to K.
The essence of extraneous solid-phase particles invading a reservoir and causing clogging is migration and deposition after invasion of the exogenous solid-phase particles into a medium. Thus, the core of each embodiment of the present invention is to establish a kinetic model of migration and deposition of the solid-phase particles. Specifically, based on mass conservation, a diffusion relationship, and the like, a spatio-temporal evolution control phenomenological model (containing a concentration C of flowing particles and a concentration Cd of deposited particles) for concentration distribution of the extraneous solid-phase particles in a reservoir around a well to be diagnosed is established, and in conjunction with a relationship between a deposition concentration and a characteristic parameter characterizing the damage extent of the reservoir such as permeability, spatio-temporal field distribution of the characteristic parameter such as permeability can be diagnosed.
Step S1101: determining a velocity of a fluid containing flowing particles in a reservoir.
The reservoir is located in a preset region of a well to be diagnosed (e.g., a water injection well, an oil production well, or the like).
For the step S1101, the determining a velocity of a fluid containing flowing particles in a reservoir may include: establishing a pressure conduction equation for the fluid entering the reservoir; and determining the velocity of the fluid according to the pressure conduction equation and a Darcy formula.
Specifically, a pressure is a force that drives a solid-liquid mixture (i.e., a fluid containing flowing particles) to continuously intrude from a wellbore of the water injection well into the reservoir around the well to be diagnosed, whereby the pressure conduction equation for the fluid entering the reservoir as expressed in formula (1-1) can be established:
and then the velocity of the fluid can be determined according to formula (1-1) and the Darcy formula (e.g., formula (1-2) below),
where P({right arrow over (r)}, t) is a pressure of the fluid; ϕ is a porosity of the reservoir; μ is fluid viscosity; ct is a fluid-rock integrated compression coefficient; K({right arrow over (r)}, t) is a permeability of the reservoir; and τ is a tortuosity of the reservoir.
Step S1102: establishing a mass balance equation between the fluid and deposited particles on the reservoir, based on a convection parameter and a diffusion parameter of the fluid.
According to conservation of mass, a change in the mass of the fluid is equal to a negative change in the mass of the deposited particles. For the step S1102, the establishing a mass balance equation between the fluid and the deposited particles on the reservoir may include: establishing the mass balance equation expressed in the following formula (1-3), based on the convection parameter and the diffusion parameter of the fluid,
where ρ is density of the fluid; ϕ is the porosity of the reservoir; w({right arrow over (r)}, t) is the mass fraction (which may also be called a mass concentration) of the flowing particles; u is a Darcy apparent velocity; j({right arrow over (r)}, t) is a diffusion flow rate, j({right arrow over (r)}, t)=−ϕρLD({right arrow over (r)}, t)∇w({right arrow over (r)}, t), where ρL is a density of the fluid, D({right arrow over (r)}, t) is a diffusion coefficient, D({right arrow over (r)}, t)=αν({right arrow over (r)}, t), α is a vertical diffusivity, and ν({right arrow over (r)}, t) is a velocity of the fluid;
where {dot over (m)} ({right arrow over (r)}, t) is an accumulated mass of the deposited particles per unit time, and t is time.
Step S1103: establishing a connection condition equation between a volume concentration of the deposited particles and a volume concentration of the fluid, based on the convection parameter and the diffusion parameter of the fluid.
For the step S1103, the establishing a connection condition equation between a volume concentration of the deposited particles and to volume concentration of the fluid may include: establishing the connection condition equation expressed in the following formula (1-4), based on the convection parameter and the diffusion parameter of the fluid,
where ρP is particle density; Cd({right arrow over (r)}, t) is the volume concentration of the deposited particles; and k({right arrow over (r)}, t)=k0({right arrow over (r)})G1(Cd)F1(T), where k0 is an original fluid loss coefficient,
Since the correlation between F1(T) and temperature is measured by exp(1/T), and in a common temperature range (e.g. 300 K to 400 K), the change of this function is actually very slow and actually close to an isothermal process, thus
where Cd({right arrow over (r)}, t) is the volume concentration of the deposited particles, Cd max is a maximum volume concentration of the deposited particles, and mk is a first empirical value. All of the above parameters can be either constants, or parameters that vary with space, i.e., in a non-homogeneous situation.
Step S1104: determining a spatio-temporal evolution simulation equation of reservoir damage by the particles according to a relationship between a mass fraction of the flowing particles and a volume concentration of the flowing particles, the velocity of the fluid, the mass balance equation and the connection condition equation.
Wherein the spatio-temporal evolution simulation equation is used to simulate a four-dimensional (4D) spatio-temporal evolution process of characteristics of reservoir damage caused by the extraneous solid-phase particles.
Wherein the relationship between the mass fraction of the flowing particles and the volume concentration of the flowing particles may be
where ρP is the density of the deposited particles; ρL is the density of the fluid; and C({right arrow over (r)}, t) is the volume concentration of the flowing particles. The spatio-temporal evolution simulation equation of reservoir damage by the particles may include: a spatio-temporal evolution simulation equation of reservoir damage by particle migration expressed by formula (1-5), and a spatio-temporal evolution simulation equation of reservoir damage by particle deposition expressed by formula (1-6).
For the step S1104, the determining a spatio-temporal evolution simulation equation of reservoir damage by the particles may include: determining the spatio-temporal evolution simulation equation of reservoir damage by particle migration expressed by formula (1-5) according to the relationship between the mass fraction of the flowing particles and the volume concentration of the flowing particles, the velocity of the fluid, and the mass balance equation expressed by formula (1-3):
and determining the spatio-temporal evolution simulation equation of reservoir damage by particle deposition expressed by formula (1-6) according to the relationship between the mass fraction of the flowing particles and the volume concentration of the flowing particles, the velocity of the fluid, and the connection condition equation expressed by formula (1-4):
where C({right arrow over (r)}, t) is the volume concentration of the flowing particles; ν({right arrow over (r)}, t) is the velocity of the fluid; τ is the tortuosity of the reservoir; ρP is the density of the deposited particles; ρL is the density of the fluid;
where k0({right arrow over (r)}) is an initial value of the fluid loss coefficient of the reservoir; Cd({right arrow over (r)}, ti+1) is the volume concentration of the deposited particles; Cd max is a maximum volume concentration of the deposited particles; mk is a first empirical value; α is a vertical diffusivity; and ϕ is the porosity of the reservoir. k0({right arrow over (r)})=ƒ(NR, NPe, NA, NDL, NE1NE2, NG, NLo, NvdW, ζp(g)), where NR, NPe, NA, NDL, NE1, NE2, NG, NLo, NvdW, ζp(g) are a radius number, a Peclet number, an attraction number, an electrical double layer number, a first electric potential force number, a second electric potential force number, a gravity number, a London force number, a van der Waals force number, and potentials of flowing particles and matrix particles, respectively (for details of relevant expressions of the parameters, see Table 2).
In summary, according to the present invention, the mass balance equation between the fluid containing the flowing particles and the deposited particles on the reservoir is creatively established; the connection condition equation between the volume concentration of the deposited particles and the volume concentration of the fluid is established; and the spatio-temporal evolution simulation equation of reservoir damage by the particles is determined according to the relationship between the mass fraction of the flowing particles and the volume concentration of the flowing particles, the velocity of the fluid, the mass balance equation and the connection condition equation. Thus, by using the determined spatio-temporal evolution simulation equation, a four-dimensional spatio-temporal evolution process of characteristics of reservoir damage caused by extraneous solid-phase particles can be quantitatively simulated, which is of very great significance for optimal design of a declogging measure and improvement or restoration of oil well production and water well injection capacity for damaged wells, and improvement of numerical simulation precision of oil pools.
Step S1201: determining the volume concentration of the deposited particles based on the spatio-temporal evolution simulation equation established by the modeling method for reservoir damage by the extraneous solid-phase particles.
For the spatio-temporal evolution simulation equation of reservoir damage by particle migration expressed by formula (1-5) described above, in a one-dimensional situation, such an equation can be arranged into the following general form:
where aa, bb, cc may be a constant (e.g., the diffusion coefficient) or a function (e.g., the velocity of the fluid); and ƒ may be a pressure, substance concentration, stress, or the like. A backward difference is used for time and a central difference is used for space. Then the above equation may use the following difference format:
where i=1, 2, 3 . . . Ni,
n=1, 2, 3, . . . , t=nΔt, Ni being the number of discrete space points.
A solution interval is x ∈ (0, xmax), and Δx and Δt are space and time steps. Simultaneously, an initial condition ƒin|n=0=ƒi0, i=1, 2, 3 . . . , Ni and a boundary condition (ƒin|i=1=ƒ0, n=1, 2, 3 . . . (at a well wall) and
(at a boundary of a preset range or several meters from the well wall, with a virtual grid i+1 being constructed) are considered.
First, for i=2, 3, . . . , Ni−1, the above difference format is arranged as follows:
where A1i, A2i and A3i are respectively:
Meanwhile, according to formula (1-5), we obtain:
Substituting formula (1-11) into formula (1-10) yields a specific expression of the iterative relational expression (1-9), which is not listed here as the specific expression of the iterative relational expression (1-9) is complicated. Then, iterative calculation is performed by using the initial condition and the boundary condition to obtain a value of the field ƒ.
Next, a difference solution procedure for illustrating the boundary condition is described.
The above iterative relational expression (1-9) is applicable to a non-boundary grid. For i=1 (at the well wall), as a point-centered grid is used and it is a Dirichlet boundary condition, the following relational expression can be obtained directly:
ƒ1nƒ0(constant), i=1. (1-12)
For i=N (several meters from the well wall, at a boundary of a preset range), which is a Neumann or second-type (Neumann) boundary condition, a virtual grid i=Ni+1 is added, and ƒi+1n=ƒi−1n is obtained from
and ƒi+1n=ƒi−1n is substituted into formula (1-9) to yield:
ƒin=(A1i+A2i)ƒi−1n+AQ3iƒin−1(i=N). (1-13)
The spatio-temporal variation of the field function ƒ can be solved according to the above process. Since the above numerical model is built for a reservoir near the wellbore of the well to be diagnosed, a cylindrical coordinate system needs to be used when the distribution of a physical quantity ƒ around the well is solved. Thus, formula
needs to be transformed into
This form is not conducive to an equally spaced difference, so that a coordinate transformation can be introduced: r=rwex′, where rw is a radius of the wellbore, and x′ is a dimensionless spatial coordinate. Substituting the transformation into a general equation yields an equation about x′:
serve as new equation coefficients, the above formula is essentially same as
Thus, an equally spaced difference may be applied to the x′ coordinate, and the iterative format described before is used. After the value of ƒ is calculated, the spatial coordinate may be mapped back from x′ to r to obtain ƒ(r, t).
The spatio-temporal evolution simulation equation established by the above-mentioned modeling method for reservoir damage by the extraneous solid-phase particles comprehensively considers the influence of various physicochemical factors on the reservoir damage when the solid-phase particles invade the reservoir, and thus the volume concentration of the deposited particles obtained by the solution of this step S1201 is very precise.
Step S1202: determining a characteristic parameter characterizing the damage extent of the reservoir in a preset region of a well to be diagnosed, based on the volume concentration of the deposited particles.
Wherein the characteristic parameter may be permeability of the reservoir and/or a fluid loss coefficient of the reservoir.
In an embodiment, the characteristic parameter may be the permeability of the reservoir.
For the step S1202, the determining a characteristic parameter characterizing the damage extent of the reservoir in a preset region of a well to be diagnosed may include: determining the permeability K({right arrow over (r)}, t) the reservoir based on the volume concentration Cd({right arrow over (r)}, t) of the deposited particles and formula (1-15),
In an embodiment, the characteristic parameter may be the fluid loss coefficient of the reservoir.
For the step S1202, the determining a characteristic parameter characterizing the damage extent of the reservoir in a preset region of a well to be diagnosed may include: determining the fluid loss coefficient k({right arrow over (r)}, t) of the reservoir based on the volume concentration Cd({right arrow over (r)}, t) of the deposited particles and formula (1-16),
where ϕ0 is an initial value of the porosity; Cd max is a maximum volume concentration of the deposited particles; mk and mK are a first empirical value and a second empirical value, respectively; K0 ({right arrow over (r)}) is an initial value of the permeability of the reservoir; and k0 ({right arrow over (r)}) is an initial value of the fluid loss coefficient of the reservoir.
wherein the characteristic parameter may be a skin factor of the reservoir.
For the step S1202, the determining a characteristic parameter characterizing the damage extent of the reservoir in a preset region of a well to be diagnosed may include: determining the permeability K({right arrow over (r)}, t) of the reservoir based on the volume concentration Cd({right arrow over (r)}, t) of the deposited particles and formula
and determining the skin factor s({right arrow over (r)}, t) of the reservoir based on the permeability K({right arrow over (r)}, t) of the reservoir and formula (1-17),
where K0 ({right arrow over (r)}) is the initial value of the permeability of the reservoir; and
The characteristic parameter obtained by the step S1202 (e.g., the permeability K({right arrow over (r)}, t) and the skin factor S({right arrow over (r)}, t) of the reservoir) is a result of 4D quantitative simulation of spatio-temporal evolution. Quantitative prediction of reservoir damage and spatio-temporal deduction of damage laws can be performed according to its evolution characteristics, which is of scientific guidance significance for preventing or avoiding reservoir damage, and formulating development plans for oil pools and subsequent well stimulation measures.
In summary, the volume concentration of the deposited particles can be determined by using the determined spatio-temporal evolution simulation equation, and then a characteristic parameter (e.g., the permeability and/or the skin factor of the reservoir) characterizing the damage extent of the reservoir in the preset region of the well to be diagnosed is determined based on the volume concentration of the deposited particles, whereby the four-dimensional spatio-temporal evolution process of the characteristics of reservoir damage caused by the extraneous solid-phase particles can be simulated quantitatively. Therefore, for a well with reservoir damage, quantitative simulation and spatial-temporal evolution of reservoir damage are achieved by using historical parameters, which is of great significance for optimal design of a declogging measure and improvement of numerical simulation precision of oil pools; and for a well without reservoir damage, quantitative prediction of reservoir damage and spatio-temporal deduction of damage laws can be performed according to physical property parameters and engineering parameters to be implemented, which is of scientific guidance significance for preventing or avoiding reservoir damage, and formulating development plans for oil pools and subsequent well stimulation measures.
A process of diffusion of water molecules of an extraneous fluid (e.g., injected water) through a solid-liquid interface (between rock and the fluid in a reservoir) to the interior of a solid-phase medium (the rock in the reservoir) can be considered as a phenomenon in any local region of the solid-liquid interface, and for any sufficiently small-scale local region on the interface, a diffusion direction of the water molecules can be considered to be perpendicular to a tangential direction of a point (e.g., point O) in the region (i.e., the diffusion direction is perpendicular to a plane where the region is located), as shown in
Step S2101: determining a Darcy apparent velocity of a fluid in a reservoir in a preset region of a well to be diagnosed.
wherein the well to be diagnosed may be, for example, a water injection well or oil production well.
For the step S2101, determining the velocity of the fluid in the reservoir may include: establishing a pressure conduction equation for the fluid entering the reservoir; and determining the Darcy apparent velocity of the fluid according to the pressure conduction equation and a Darcy formula.
Specifically, a pressure is a force that drives a solid-liquid mixture to continuously intrude from a wellbore of the water injection well into the reservoir around the well, whereby the pressure conduction equation for the fluid entering the reservoir as expressed in formula (2-1) can be established:
and then the Darcy apparent velocity of the fluid can be determined according to formula (2-1) and the Darcy formula (e.g., formula (2-2) below),
where P({right arrow over (r)}, t) is the pressure of the fluid; ϕ is the porosity of the reservoir; μ is fluid viscosity; ct is a fluid-rock integrated compression coefficient; and K({right arrow over (r)}, t) is the permeability of the reservoir.
Step S2102: establishing a mass balance equation for water molecules in the fluid according to the Darcy apparent velocity of the fluid and a diffusion coefficient of the water molecules in the fluid.
Under reservoir conditions, water contents at different locations within pores in the reservoir satisfy a mass conservation equation. The variation of the water contents within the reservoir is mainly determined by two processes: convection and diffusion. Specifically, for the step S2102, the establishing a mass balance equation for water molecules in the fluid may include: establishing the mass balance equation expressed in the following formula (2-3) according to the Darcy apparent velocity u of the fluid and the diffusion coefficient Dw of the water molecules:
where ϕ0 is an initial value of the porosity of the reservoir; c1({right arrow over (r)}, t) is a water volume fraction of the pores in the reservoir; and {right arrow over (r)} is a spatial location of any point in the reservoir.
An initial condition for the mass balance equation for the water molecules in the fluid is c1({right arrow over (r)}, t=0)=ϕ0·Swc, for the mass balance equation for the water molecules in the fluid is c1(|{right arrow over (r)}=rw, t)=ϕ0 (that is, the pores in the reservoir in a well wall of the water injection well is completely filled with water, i.e., the water saturation in the pores is 1), where ϕ0 is the initial value of the porosity of the reservoir; rw is a wellbore radius of the well to be diagnosed; and Swc is an irreducible water saturation in the reservoir.
Step S2103: establishing a diffusion equation for diffusion of the water molecules in the fluid to the interior of rock in the reservoir according to Fick's law of diffusion.
It should be noted that c(n, t) is a water volume fraction of rock in the reservoir at time t, at a coordinate n in a one-dimensional coordinate system (the direction of the coordinate axis is a normal direction of the solid-liquid interface at the location {right arrow over (r)}) established with a location {right arrow over (r)} (e.g., point O in
A one-dimensional coordinate system n that is perpendicular to the solid-liquid interface and points to the interior of the solid phase can be established, wherein n=0 at the interface, and n>0 inside the solid phase, as shown in
where n is a coordinate in the one-dimensional coordinate system established with the location {right arrow over (r)} in the reservoir as the origin and with the normal direction of the interface between the fluid and the rock at the location {right arrow over (r)} as the coordinate axis; t is time; Dw is a diffusion coefficient of the water molecules; and c(n, t) is a water volume fraction of the rock in the reservoir.
Wherein an initial condition of the diffusion equation is c(n, t=0)=c0; and a boundary condition of the diffusion equation is
Dw is the diffusion coefficient of the water molecules; c(n, t) is the water volume fraction of the rock in the reservoir; and kƒ is a membrane exchange coefficient.
Step S2104: determining a spatio-temporal evolution simulation equation of reservoir damage by the clay swelling according to the diffusion equation and the mass balance equation.
For the step S2104, as shown in
Step S2301: determining a water volume fraction of pores in the reservoir according to the mass balance equation, and the boundary condition and the initial condition of the mass balance equation.
The water volume fraction c1({right arrow over (r)}, t) of the pores in the reservoir can be obtained according to the above formula (2-3) and the initial condition and the boundary condition of the mass balance equation for the water molecules in the fluid.
Step S2302: determining a water absorption rate of the rock in the reservoir according to the diffusion equation, the boundary condition of the diffusion equation and the initial condition of the diffusion equation, and the water volume fraction of the pores in the reservoir.
First, the water volume fraction c(n, t) of the rock in the reservoir can be solved according to the above equations (2-4) and the boundary condition of the diffusion equation and the initial condition of the diffusion equation,
where
is a complementary error function, and L−1{·} represents an inverse Laplace transform.
Then, {dot over (A)} ({right arrow over (r)}, t) can be obtained according to c1 ({right arrow over (r)}, t), c(n, t) and the definition of the water absorption rate of the rock in the reservoir
Specifically, the relative magnitudes of C1 and C0 determine whether {dot over (A)} is positive or negative. If c1>c0, then {dot over (A)}>0, indicating that the water content in the pores is greater than the water content in the solid phase and water will diffuse into the solid phase; conversely, if c1<c0, {dot over (A)}<0 according to the two formulas, which means that the rock solid phase loses water. Finally, under reservoir conditions, the water content in the pores is always greater than or equal to the water content in the solid phase, and a limiting condition should be added for {dot over (A)} such that it is equal to 0 when c1<c0. Thus {dot over (A)} is expressed as:
Step S2303: determining a spatio-temporal evolution simulation equation of reservoir damage by the clay swelling according to the water absorption rate of the rock in the reservoir.
For the step S2303, the determining a spatio-temporal evolution simulation equation of reservoir damage by the clay swelling may include: determining the spatio-temporal evolution simulation equation of reservoir damage by the clay swelling expressed by the following formula (2-6) according to the water absorption rate {dot over (A)} ({right arrow over (r)}, t) of the reservoir:
where ϕ ({right arrow over (r)}, t) is the porosity of the reservoir; and λ is a clay swelling coefficient.
Specifically, the clay expansion coefficient
where Cc is a mass percentage of clay in the rock; PI is a plasticity coefficient of the rock (dimensionless), where if PI<1˜2, the rock is brittle rock; if 2<PI<6, the rock is plastic-brittle rock; and if PI>6, the rock is plastic rock; and k′ is an empirical parameter.
{dot over (A)}({right arrow over (r)}, t) can be obtained according to the above formula (2-6). If {dot over (A)}>0 (i.e., the water absorption rate is positive), the clay swells so that
the porosity decreases.
In summary, according to the present invention, the Darcy apparent velocity of the fluid in the reservoir in a preset region of a well to be diagnosed is creatively determined; the mass balance equation for water molecules in the fluid is established according to the Darcy apparent velocity of the fluid and the diffusion coefficient of the water molecules in the fluid; the diffusion equation for diffusion of the water molecules in the fluid to the interior of rock in the reservoir is established; and the spatio-temporal evolution simulation equation of reservoir damage by the clay swelling is determined according to the diffusion equation and the mass balance equation. Thus, using the determined spatio-temporal evolution simulation equation, a four-dimensional spatio-temporal evolution process of characteristics of reservoir damage caused by the clay swelling can be quantitatively simulated. Therefore, performing quantitative prediction of reservoir damage and spatio-temporal deduction of damage laws is of scientific guidance significance for preventing or avoiding reservoir damage, and formulating development plans for oil pools and subsequent well stimulation measures for a well without reservoir damage, and is of very great significance for optimal design of a declogging measure and improvement or restoration of oil well production and water well injection capacity for damaged wells, and improvement of numerical simulation precision of oil pools.
Step S2401: determining porosity of the reservoir based on the spatio-temporal evolution simulation equation established by the modeling method for reservoir damage by clay swelling.
For the solution of the spatio-temporal evolution simulation equation for reservoir damage by the clay swelling expressed by the above formula (2-6), c1 ({right arrow over (r)}, t) needs to be calculated according to formula (2-3). For the specific solving process, reference can be made to the solving process of the volume concentration of the deposited particles in the above Embodiment 1, which will not be described here.
After the water volume fraction c1 (r, t) of the pores in the reservoir is calculated by the above method, the water absorption rate {dot over (A)} of the reservoir can be calculated according to the above formula (2-5) (
Step S2402: determining a characteristic parameter characterizing the damage extent of the reservoir in a preset region of a well to be diagnosed, based on the determined porosity of the reservoir.
Wherein the characteristic parameter may be the permeability of the reservoir.
For the step S2402, the determining a characteristic parameter characterizing the damage extent of the reservoir in a preset region of a well to be diagnosed may include: determining the permeability K({right arrow over (r)}, t) of the reservoir based on the porosity ϕ({right arrow over (r)}, t) of the reservoir and formula (2-7):
where ϕ0 is an initial value of the porosity of the reservoir; mK is a second empirical value; and K0 ({right arrow over (r)}) is an initial value of the permeability of the reservoir.
The characteristic parameter may be a fluid loss coefficient of the reservoir.
For the step S2402, the determining a characteristic parameter characterizing the damage extent of the reservoir in a preset region of a well to be diagnosed may include: determining the fluid loss coefficient k({right arrow over (r)}, t) of the reservoir based on the porosity ϕ({right arrow over (r)}, t) of the reservoir and formula (2-8):
where ϕ0 is the initial value of the porosity of the reservoir; ϕd max is maximum porosity of the reservoir, mk is a first empirical value; and k0 ({right arrow over (r)}) is an initial value of the fluid loss coefficient of the reservoir.
In the case where the characteristic parameter is the permeability of the reservoir and the fluid loss coefficient of the reservoir, the permeability of the reservoir can be determined by formula (2-7), and the fluid loss coefficient of the reservoir is determined by formula (2-8).
Wherein the characteristic parameter may be a skin factor of the reservoir.
For the step S2402, the determining a characteristic parameter characterizing the damage extent of the reservoir in a preset region of a well to be diagnosed may include: determining the permeability K({right arrow over (r)}, t) of the reservoir based on the porosity ϕ(r, t) of the reservoir and formula
and determining the skin factor S ({right arrow over (r)}, t) of the reservoir based on the permeability K ({right arrow over (r)}, t) of the reservoir and formula (2-9):
where K0 ({right arrow over (r)}) is the initial value of the permeability of the reservoir; and
The characteristic parameter obtained by the step S2402 (e.g., the permeability K({right arrow over (r)}, t) and the skin factor S({right arrow over (r)}, t) of the reservoir) is a result of 4D quantitative simulation of spatio-temporal evolution (as shown in
where Kmax({right arrow over (r)}, t) is a maximum value of K({right arrow over (r)}, t)), and a working person concerned can visually confirm the damage extent of the reservoir from
In summary, according to the present invention, the porosity of the reservoir can be creatively determined by using the determined spatio-temporal evolution simulation equation, then the characteristic parameter (e.g., the permeability and/or the skin factor of the reservoir) characterizing the damage extent of the reservoir in the preset region of the well to be diagnosed can be determined according to the porosity of the reservoir, and thus, a four-dimensional spatio-temporal evolution process of characteristics of reservoir damage caused by the clay swelling can be quantitatively simulated. Therefore, performing quantitative prediction of reservoir damage and spatio-temporal deduction of damage laws is of scientific guidance significance for preventing or avoiding reservoir damage, and formulating development plans for oil pools and subsequent well stimulation measures for a well without reservoir damage, and is of very great significance for optimal design of a declogging measure and improvement or restoration of oil well production and water well injection capacity for damaged wells, and improvement of numerical simulation precision of oil pools.
Inorganic precipitation may occur when an extraneous fluid is incompatible with a fluid in a reservoir. The inorganic precipitation can clog fluid flow channels, thereby causing reservoir damage. General inorganic precipitation may include: calcium carbonate (CaCO3), calcium sulfate (CaSO4), strontium sulfate (SrSO4), barium sulfate (BaSO4) and other inorganic precipitates.
Step S3101: determining a Darcy apparent velocity of a fluid in a reservoir in a preset region of a well to be diagnosed (e.g., a water injection well or an oil production well).
For the step S3101, determining the velocity of the fluid in the reservoir may include: establishing a pressure conduction equation for the fluid entering the reservoir; and determining the Darcy apparent velocity of the fluid according to the pressure conduction equation and a Darcy formula.
For the specific determination process, reference can be made to the process of determining a Darcy apparent velocity in the above Embodiment 2 (i.e., the above formulas (2-1) and (2-2) and related description thereof).
Step S3102: determining an ion concentration loss corresponding to each ion of a plurality of ions in an extraneous fluid.
Wherein the ion concentration loss is caused by a precipitation reaction between the each ion and a corresponding ion in the fluid in the reservoir.
For the step S3102, the determining an ion concentration loss corresponding to each ion of a plurality of ions in an extraneous fluid may include: determining the ion concentration loss Li({right arrow over (r)}, t) corresponding to the each ion i according to a concentration Ci({right arrow over (r)}, t) of the each ion i, a concentration Cj({right arrow over (r)}, t) of an ion j of at least one ion that undergoes the precipitation reaction with the each ion i, and the following equation (3-1),
where kij is a reaction rate of the each ion i and the ion j; Ni is the number of at least one ion that undergoes the precipitation reaction with the each ion i; and βi and βj are reaction coefficients of the ion i and the ion j, respectively.
Using SO42− in the extraneous fluid as an example, the SO42− can react with Ca2+, Ba2+ and Sr2+ in the fluid in the reservoir to produce precipitates. At time t, at any spatial location {right arrow over (r)} within the reservoir, an ion concentration loss resulting from the precipitation reaction between SO42− in the extraneous fluid and Ca2+, Ba2+ and Sr2+ in the fluid in the reservoir is
where
are reaction rates of the precipitation reactions between SOr2− and Ca2+, Ba2+ and Sr2+, respectively;
are concentrations of the ions SO42−, Ca2+, Ba2+ and Sr2+ at time t, at the reservoir space {right arrow over (r)}, respectively; and
are a reaction coefficient of the ion SO42− and the ion Ca2+, a reaction coefficient of the ion SO42− and the ion Ba2+, and a reaction coefficient of the ion SO42− and the ion or Sr2+, respectively. Similarly, an ion concentration loss in the reservoir resulting from a precipitation reaction between each of other ions in the extraneous fluid and a corresponding ion in the reservoir can be determined.
Specifically, the reaction rate kij of the each ion i and the ion j is determined by a scaling index of a corresponding precipitate produced by the each ion i and the ion j. For example, the reaction rate kij of the each ion i and the ion j may satisfy the following relational expression (3-2):
where kij0 is a constant; and ISij is the scaling index of the corresponding precipitate produced by the each ion i and the ion j. More specifically, the scaling index ISij may be determined by the concentration of the each ion i, the ionic strength of the fluid, the concentration of the ion j, the temperature of the fluid and the pressure of the fluid.
Specifically, at a certain temperature T, pressure P and ion concentration ([Me] is a concentration of a (free) cation, wherein the cation may be a calcium ion, a strontium ion, a barium ion, or the like, and [An] is a concentration of a (free) anion, wherein the anion may be a bicarbonate ion, a sulfate ion, or the like), whether a precipitation reaction in a reservoir solution occurs is usually determined by the scaling index, which is expressed as:
where Si is the ionic strength of the fluid; and KC is a solubility product coefficient of the precipitation reaction.
If IS≤0, the solution is in an undersaturated or saturated state, an inorganic precipitate is generated; if IS>0, the solution is in a supersaturated state, there is a tendency to generate an inorganic precipitate. Is varies with the reservoir location, ion concentration, temperature and pressure, and is a function related to time and space.
According to a Tomson-Oddo calculation method, scaling indices of four inorganic precipitates under reservoir conditions are respectively as follows:
(1) Calcium carbonate CaCO3:
Is(P,T,Si)=log([Ca2+][CO3−])+pH−2.42+0.02T−1.53×10−5T2−6.33×1031 3P−2.02Si1/2+0.727Si
(2) Calcium sulfate CaSO4:
(3) Barium sulfate BaSO4:
IS(P,T,Si)[BaSO4]=log([Ba2+][SO42−])+10.03−4.8×10−3T+11.4×10−6T2−4.8×10−5P−2.62Si1/2+0.89Si−2.0×10−3Si1/2T,
(4) Strontium sulfate SrSO4:
IS(P,T,Si)[SrSO4]=log([Sr2+][SO42−])+3.11+2.0×10−3T+6.4×10−6T2−4.6×10−5P−1.89Si1/2+0.67Si−1.9×10−3Si1/2T.
The pH in the above formulas is a pH value of a liquid in the reservoir (a liquid formed by mixing an original liquid in the reservoir with an extraneous liquid).
Step S3103: establishing a migration equation for the each ion according to a Darcy apparent velocity of the fluid, the ion concentration loss corresponding to the each ion, and a diffusion coefficient of the each ion.
In an ion concentration control equation, a flow J of ion migration includes two parts: convection and diffusion, using the ion i as an example:
Ji=Jid+Jic=−Di∇Ci+uCi,
where Ji is a migration flow of the ion i; Jid is a diffusion flow of the ion i; Jic is a convection flow of the ion i; u is the Darcy apparent velocity of the fluid; Ci is the concentration of the each ion i; and Di is the diffusion coefficient of the each ion i.
For the step S3103, the precipitation reaction results in an ion concentration loss Li({right arrow over (r)}, t) at a reservoir space {right arrow over (r)} at time t, and a migration equation for the each ion i that can be established according to the law of mass conservation may include:
where u is the Darcy apparent velocity of the fluid; ϕ({right arrow over (r)}, t) is the porosity of the reservoir; and Di is the diffusion coefficient of the each ion i.
Step S3104: determining a spatio-temporal evolution simulation equation of reservoir damage by the inorganic precipitation according to the migration equation for the each ion and a reaction coefficient of a precipitate produced by the each ion.
Wherein the spatio-temporal evolution simulation equation is used to simulate a four-dimensional spatio-temporal evolution process of characteristics of reservoir damage caused by the corresponding precipitates produced by the plurality of ions.
For the step S3104, the determining a spatio-temporal evolution simulation equation of reservoir damage by the inorganic precipitation may include: determining the concentration Ci({right arrow over (r)}, t) of the each ion i according to the mass balance equation (3-4) of the each ion; and determining the spatio-temporal evolution simulation equation of reservoir damage by the inorganic precipitation expressed by the following formula (3-5) according to the concentration Ci({right arrow over (r)}, t) of the each ion i and a reaction coefficient coij of a precipitate produced by the each ion:
where dij({right arrow over (r)}, t) is an accumulated concentration of the precipitate produced by the precipitation reaction between the ion i and the ion j at time t and at the reservoir space {right arrow over (r)}; and dij({right arrow over (r)}, t+dt) is an accumulated concentration of the precipitate produced by the precipitation reaction between the ion i and the ion j at time t+dt and at the reservoir space {right arrow over (r)}.
In summary, according to the present invention, the migration equation for the each ion of the plurality of ions in the extraneous fluid is creatively established according to the Darcy apparent velocity of the fluid in the reservoir in the preset region of the well to be diagnosed, the ion concentration loss corresponding to the each ion, and the diffusion coefficient of the each ion; and the spatio-temporal evolution simulation equation of reservoir damage by the inorganic precipitation is determined according to the migration equation for the each ion and the reaction coefficient of the precipitate produced by the each ion. Thus, by using the determined spatio-temporal evolution simulation equation, a four-dimensional spatio-temporal evolution process of characteristics of reservoir damage caused by the inorganic precipitation can be quantitatively simulated. Therefore, performing quantitative prediction of reservoir damage and spatio-temporal deduction of damage laws is of scientific guidance significance for preventing or avoiding reservoir damage, and formulating development plans for oil pools and subsequent well stimulation measures for a well without reservoir damage, and is of very great significance for optimal design of a declogging measure and improvement or restoration of oil well production and water well injection capacity for damaged wells, and improvement of numerical simulation precision of oil pools.
Step S3201: determining a volume concentration of inorganic precipitation based on the spatio-temporal evolution simulation equation established by the modeling method for reservoir damage by the inorganic precipitation, and a molar mass and density of each precipitate produced in the reservoir.
Wherein the each precipitate is produced by a precipitation reaction of an each ion of a plurality of ions in an extraneous fluid and a corresponding ion in a fluid in the reservoir, and the volume concentration of inorganic precipitation is a total volume concentration of the each precipitate.
For the migration equation of the each ion in reservoir damage by the inorganic precipitation expressed by the above formula (3-4), a water volume fraction c1 ({right arrow over (r)}, t) of pores in the reservoir can be solved by referring to the process of solving a volume concentration of deposited particles in the above Embodiment 1, which will not be described here.
After the concentration Ci({right arrow over (r)}, t) of the each ion i is calculated by the above method, an accumulated concentration dij({right arrow over (r)}, t) of each precipitate can be calculated according to the above formula (3-5); and then a volume concentration Cd
of all precipitates can be determined. where Ni is the number of ions which are located in the reservoir and suffer from a precipitation reaction with the ion i; and M is the number of a plurality of ions in the extraneous fluid. The spatio-temporal evolution simulation equation established by the above modeling method for reservoir damage by the inorganic precipitation comprehensively considers the influence of various physicochemical factors on the reservoir damage when inorganic precipitation occurs in the reservoir, and thus the volume concentration of the inorganic precipitation obtained by the solution of the step S3201 is very precise.
Step S3202: determining a characteristic parameter characterizing the damage extent of the reservoir in a preset region of a well to be diagnosed, based on the determined volume concentration of inorganic precipitation.
Wherein the characteristic parameter may be permeability of the reservoir and/or a fluid loss coefficient of the reservoir.
In an embodiment, the characteristic parameter may be the permeability of the reservoir.
For the step S3202, the determining a characteristic parameter characterizing the damage extent of the reservoir in a preset region of a well to be diagnosed includes: determining the permeability K({right arrow over (r)}, t) of the reservoir based on the volume concentration Cd({right arrow over (r)}, t) of the inorganic precipitation and formula (3-6):
In an embodiment, the characteristic parameter may be the permeability of the reservoir.
For step S3202, the determining a characteristic parameter characterizing the damage extent of the reservoir in a preset region of a well to be diagnosed includes: determining the fluid loss coefficient k({right arrow over (r)}, t) of the reservoir based on the volume concentration Cd({right arrow over (r)}, t) of the inorganic precipitation and formula (3-7):
where ϕ0 is an initial value of the porosity of the reservoir; Cd max is a maximum volume concentration of the inorganic deposition; mk and mK are a first empirical value and a second empirical value, respectively; Ko ({right arrow over (r)}) is an initial value of the permeability of the reservoir; and k0 ({right arrow over (r)}) is an initial value of the fluid loss coefficient of the reservoir.
Wherein the characteristic parameter is a skin factor of the reservoir.
For the step S3202, the determining a characteristic parameter characterizing the damage extent of the reservoir in a preset region of a well to be diagnosed may include: determining the permeability K({right arrow over (r)}, t) of the reservoir based on the volume concentration Cd({right arrow over (r)}, t) of the inorganic precipitation and formula
and determining the skin factor S({right arrow over (r)}, t) of the reservoir based on the permeability K({right arrow over (r)}, t) of the reservoir and formula (3-8):
where Ko ({right arrow over (r)}) is the initial value of the permeability of the reservoir; and
The characteristic parameter (e.g., the permeability K({right arrow over (r)}, t) and the skin factor S({right arrow over (r)}, t) of the reservoir) obtained by the step S3202 is a result of 4D quantitative simulation of spatio-temporal evolution (not shown). Therefore, quantitative prediction of reservoir damage and spatio-temporal deduction of damage laws can be performed according to evolution characteristics of the permeability or the skin factor, which is of scientific guidance significance for preventing or avoiding reservoir damage, and formulating development plans for oil pools and subsequent well stimulation measures.
In summary, the volume concentration of the inorganic precipitation can be determined by using the determined spatio-temporal evolution simulation equation, and then the characteristic parameter (e.g., the permeability and/or the skin factor of the reservoir) characterizing the damage extent of the reservoir in the preset region of the well to be diagnosed can be determined based on the determined volume concentration of the inorganic precipitation, whereby a four-dimensional spatio-temporal evolution process of the characteristics of reservoir damage caused by the inorganic precipitation can be simulated quantitatively. Therefore, performing quantitative prediction of reservoir damage and spatio-temporal deduction of damage laws is of scientific guidance significance for preventing or avoiding reservoir damage, and formulating development plans for oil pools and subsequent well stimulation measures for a well without reservoir damage, and is of very great significance for optimal design of a declogging measure and improvement or restoration of oil well production and water well injection capacity for damaged wells, and improvement of numerical simulation precision of oil pools.
The essence of clogging by fine particles (i.e., solid-phase particles with particle sizes smaller than a preset size (e.g., 37 microns)) within a reservoir is migration and deposition of the fine particles within the reservoir. Thus, the core of embodiments of the present invention is to establish a kinetic model of migration and deposition of the fine particles within the reservoir. Specifically, based on mass conservation, a diffusion relationship, and the like, a spatio-temporal evolution control phenomenological model (containing a concentration C of migrating fine particles and a concentration Cd of deposited fine particles) of concentration distribution of the fine particles within the reservoir around the well to be diagnosed is established, and in conjunction with a relationship between a deposition concentration and a characteristic parameter characterizing the damage extent of the reservoir such as permeability, spatio-temporal field distribution of the characteristic parameter such as permeability can be diagnosed.
It should be noted that for the sake of simple description, for physical and chemical quantities that evolve with time and space in embodiments of the present invention, a variable ({right arrow over (r)}, t) may be omitted, for example, ϕw({right arrow over (r)}, t) may be shortened to ϕw; and K({right arrow over (r)}, t) may be shortened to K.
Step S4101: determining a velocity of a fluid in a reservoir.
Wherein the reservoir is located in a preset region of a well to be diagnosed (e.g., a water injection well, or an oil production well).
For the step S4101, the determining a velocity of a fluid in a reservoir may include:
establishing a pressure conduction equation for the fluid entering the reservoir; and determining the velocity of the fluid according to the pressure conduction equation and a Darcy formula.
Specifically, a pressure is a force that drives a solid-liquid mixture (i.e., a fluid containing the migrating fine particles) to continuously intrude from a wellbore of the water injection well into the reservoir around the well, whereby the pressure conduction equation for the fluid entering the reservoir as expressed in formula (4-1) can be established:
and then the velocity of the fluid can be determined according to formula (4-1) and the Darcy formula (e.g., formula (4-2) below):
where P({right arrow over (r)}, t) is the pressure of the fluid; ϕ is the porosity of the reservoir; μ is fluid viscosity; ct is a fluid-rock integrated compression coefficient; K({right arrow over (r)}, t) is the permeability of the reservoir; and τ is the tortuosity of the reservoir.
Step S4102: establishing a mass balance equation between the fluid and deposited fine particles on rock in the reservoir, based on a convection parameter and a diffusion parameter of the fluid and a mass change rate of the migrating fine particles in the fluid.
There is a correlation between the mass change rate of the migrating fine particles and the velocity of the fluid. A process of obtaining the mass change rate of the migrating fine particles is described in detail below.
In establishment of a fine particle migration damage model, first a critical velocity of the fluid when the deposited fine particles start to migrate is considered, and then how the migrating fine particles change a solid-liquid flow deposition equation is considered.
According to a fine particle starting model, forces on elastic solid fine particles with a radius rs on a rough inner surface of a rock pore is shown in
Specifically, the deposited fine particles are subjected to a dragging force Fd, a gravity Fg, an electrostatic force Fe, and a lifting force Fl in the same direction as a flow velocity ν. Since the magnitudes of both the dragging force Fd and the lifting force Fl are functions of the flow velocity ν, the corresponding flow velocity when their moments are balanced is the critical velocity νcr. According to the forces on the deposited fine particles, the moment balance equation for the deposited fine particles expressed by the following equation (4-3) is established:
Fd·ld=(Fe−Fl+Fg)·ln, (4-3)
The dragging force Fd generated by a flow field near a rough inner surface of a rock pore, on the deposited fine particles attached to this surface can be obtained from an asymptotic solution of a Navier-Stokes equation, which is expressed as follows:
Fd=ωπμrsνcr, (4-4)
where μ is the viscosity of the fluid, rs is the radius of the deposited fine particles, νcr is a flow velocity at a distance rs from the surface, and ω is a drag coefficient (e.g., ω=6×1.7 and in particular, if the value of ω is 6, the dragging force in this case corresponds to a dragging force on the solid-phase fine particle in free borderless flow).
The lifting force Fl on the deposited fine particles by a shear flow field can be expressed as follows:
Fl=χ[ρμ(rsνcr)3]1/2, (4-5)
where ρ and μ are the density and viscosity of the fluid, respectively, and χ is a lift coefficient.
The gravity Fg of the deposited fine particles can be expressed as follows:
where g is a gravitational acceleration, ρs is the density of the deposited fine particle, and ρ is the density of the fluid.
Generally, the magnitude of the electrostatic force Fe on the deposited fine particles is calculated from a derivative of an electrostatic potential with respect to space:
where V(h) is total electrostatic potential energy, which includes three parts: VLVA (London-Van der Waals potential), VDLR (electrical double layer potential) and VBR (Bonn potential). That is, V(h) can be expressed as:
V=VLVA+VDLR+VBR,
wherein VLVA and VDLR can be derived from the well-known DLVO theory:
wherein, VBR can be expressed as:
H is a Hamaker constant, h is a distance between the surface of the deposited fine particles and the surface of a medium (e.g., rock), ε0 is a dielectric constant of the deposited fine particles, De is e an electrical double layer constant of the deposited fine particles, ψ01 and ψ02 are surface potential energy of the deposited fine particles and rock framework fine particles, respectively, σI,J is a Lennard-Jones potential constant of atomic or molecular interaction, and κ is an inverse Debye length (i.e., its magnitude is 1 divided by the length).
νcr and the flow velocity ν of the fluid (the true Darcy velocity, ν in the fine particle migration model) have a relationship as follows:
in the formula, rc is an average radius of reservoir pore throats (i.e., an average radius of pores); and a relationship between the flow velocity ν of the fluid and a dragging force on the surface fine particles by the fluid can be determined by solving the equations (4-4) and (4-8) simultaneously.
A force arm In (of a normal force Fn) is approximately a radius of a contact deformation surface of the deposited fine particle with a matrix (e.g., rock) under the action of the normal force Fn(Fn=Fe−Fl+Fg):
where rs is a radius of a starting fine particle (under the scouring action of the fluid, part of the deposited fine particles start to migrate and enter the fluid to become migrating fine particles, and the part of the fine particles become starting fine particles), K is a complex Young's modulus, and
where E1 and E2 are Young's moduli of the starting fine particle and the matrix, respectively, and ν1 and ν2 are the Poisson's ratios of the starting fine particle and the matrix, respectively.
Once ln is obtained, ld can be determined from a simple geometric relationship:
ld=√{square root over (rs2−ln2)}, (4-10)
and by solving simultaneous equations of the above formulas (4-3)-(4-10), the critical velocity νcr can be obtained, which is expressed as follows:
According to formula (4-11), the critical velocity νcr is related to the mechanical, physical, and chemical properties of the deposited fine particles and the medium (e.g., rock). Only when the actual velocity of the fluid in the reservoir exceeds the critical velocity νcr, can the deposited fine particles migrate under the action of the fluid to become migrating fine particles (or a migrating material source). In general, the flow velocity of the fluid closer to the center of the wellbore of the well to be diagnosed is higher, so a fine particle migration region should be an annular band near the wellbore.
According to a mass equation, assuming that the mass change rate of the migrating fine particles (i.e., the quantity of fine particles released) is qs, qs has the following properties:
In other words, only when the velocity of the fluid (which can also be called fluid flow velocity) exceeds the critical velocity, can the fine particles incipiently move and enter the fluid to participate in migration, thereby increasing the mass of the fluid-solid mixture. Therefore, for the step S4102, the establishing a mass balance equation between the fluid and deposited fine particles on rock in the reservoir may include: establishing the mass balance equation expressed in the following formula, based on a convection parameter and a diffusion parameter of the fluid,
where ρ is the density of the fluid; ϕ is the porosity of the reservoir; w({right arrow over (r)}, t) is the mass fraction (which may also be called a mass concentration) of the deposited fine particles; u is a Darcy apparent velocity; j({right arrow over (r)}, t) is a diffusion flow rate, j({right arrow over (r)}, t)=−ϕρLD∇w({right arrow over (r)}, t), where ρL is the density of the fluid, D ({right arrow over (r)}, t) is a diffusion coefficient of the migrating fine particles, D({right arrow over (r)}, t)=αν({right arrow over (r)}, t), α is a vertical diffusivity, and ν ({right arrow over (r)}, t) is the velocity of the fluid;
is an accumulated mass of the deposited fine particles per unit time; t is time; and qs is the mass change rate of the migrating fine particles.
Wherein the mass change rate qs of the migrating fine particles is obtained by: determining the intensity Q(r) of a release field of the deposited fine particles; determining a decay function Y(t) of the intensity of the release field; and determining the mass change rate qs=Q(r)Y(t) of the migrating fine particles according to the intensity Q(r) of the release field and the decay function Y(t) of the intensity of the release field. Specifically, the intensity Q(r) of the release field may be a constant (q0), and the decay function Y(t) may be an exponential decay function that can vary with time (e.g., e−λt, where λ is a decay constant).
Step S4103: establishing a connection condition equation between a volume concentration of the deposited fine particles and a volume concentration of the fluid, based on the convection parameter and the diffusion parameter of the fluid.
For the step S4103, the establishing a connection condition equation between a volume concentration of the deposited fine particles and a volume concentration of the fluid may include: establishing the connection condition equation expressed in the following formula (4-13), based on the convection parameter and the diffusion parameter of the fluid:
where ρp is the density of the deposited fine particles; Cd({right arrow over (r)}, t) is the volume concentration of the deposited fine particles; and k({right arrow over (r)}, t)=k0({right arrow over (r)})Gl(Cd)F1(T), where k0 is an original fluid loss coefficient,
Since the correlation between F1(T) (F1(T) is an exponential function related to temperature) and temperature is measured by exp(1/T) and in a common temperature range (e.g. 300 K to 400 K), the change of this function is actually very slow and actually close to an isothermal process, thus
where Cd({right arrow over (r)}, t) is the volume concentration of the deposited fine particles, Cd max is a maximum volume concentration of the deposited fine particles, and mk is a first empirical value. All of the above parameters can be either constants, or parameters that vary with space, i.e., in a non-homogeneous situation.
Step S4104: determining a spatio-temporal evolution simulation equation of reservoir damage by the fine particles within the reservoir according to a relationship between a mass fraction of the migrating fine particles and a volume concentration of the migrating fine particles, the velocity of the fluid, the mass balance equation and the connection condition equation.
Wherein the spatio-temporal evolution simulation equation is used to simulate a four-dimensional spatio-temporal evolution process of characteristics of reservoir damage caused by the fine particles.
Wherein the relationship between the mass fraction of the migrating fine particles and the volume concentration of the migrating fine particles may be
where ρp is the density of the deposited fine particles; ρL is the density of the fluid; w({right arrow over (r)}, t) is the mass fraction of the migrating fine particles; and C({right arrow over (r)}, t) is the volume concentration of the migrating fine particles. The spatio-temporal evolution simulation equation of reservoir damage by the fine particles may include: a spatio-temporal evolution simulation equation of reservoir damage by fine particle migration expressed by formula (4-14), and a spatio-temporal evolution simulation equation of reservoir damage by fine particle deposition expressed by formula (4-15).
For the step S4104, the determining a spatio-temporal evolution simulation equation of reservoir damage by the fine particles may include: determining the spatio-temporal evolution simulation equation of reservoir damage by fine particle migration expressed by formula (4-14) according to the relationship between the mass fraction of the migrating fine particles and the volume concentration of the migrating fine particles, the velocity of the fluid, and the mass balance equation expressed by formula (4-12):
and determining the spatio-temporal evolution simulation equation of reservoir damage by fine particle deposition expressed by formula (4-15) according to the relationship between the mass fraction of the migrating fine particles and the volume concentration of the migrating fine particles, the velocity of the fluid, and the connection condition equation expressed by formula (4-13):
where C({right arrow over (r)}, t) is the volume concentration of the migrating fine particles; ν({right arrow over (r)}, t) is the velocity of the fluid; τ is the tortuosity of the reservoir; ρp is the density of the deposited fine particles; ρL is the density of the fluid;
and k0({right arrow over (r)}) is an initial value of the fluid loss coefficient of the reservoir; Cd({right arrow over (r)}, t) is the volume concentration of the deposited fine particles; Cd max is the maximum volume concentration of the deposited fine particles; mk is the first empirical value; α is the vertical diffusivity; ϕ is the porosity of the reservoir; and qs is the mass change rate of the migrating fine particles. k0({right arrow over (r)})=ƒ(NR, NPe, NA, NDL, NE1, NE2, NG, NLo, NvdW, ζp(g)), where NR, NPe, NA, NDL, NE1, NE2, NG, NLo, NvdW, ζp(g) are a radius number, a Peclet number, an attraction number, an electrical double layer number, a first electric potential force number, a second electric potential force number, a gravity number, a London force number, a van der Waals force number, and potentials of migrating fine particles and matrix particles (i.e., particles deposited on the rock), respectively (for details of relevant expressions of the parameters, see Table 2).
In summary, according to the present invention, the mass balance equation between the fluid and the deposited fine particles on rock in the reservoir is creatively established according to the convection parameter and the diffusion parameter of the fluid in the reservoir and the mass change rate of migrating fine particles; the connection condition equation between the volume concentration of the deposited fine particles and the volume concentration of the fluid is established according to the convection parameter and the diffusion parameter of the fluid; and the spatio-temporal evolution simulation equation of reservoir damage by the fine particles within the reservoir is determined according to the relationship between the mass fraction of the migrating fine particles and the volume concentration of the migrating fine particles, the velocity of the fluid, the mass balance equation and the connection condition equation. Thus, by using the determined spatio-temporal evolution simulation equation, a four-dimensional spatio-temporal evolution process of characteristics of reservoir damage caused by the fine particles within the reservoir can be quantitatively simulated. Therefore thereby, performing quantitative prediction of reservoir damage and spatio-temporal deduction of damage laws is of scientific guidance significance for preventing or avoiding reservoir damage, and formulating development plans for oil pools and subsequent well stimulation measures for a well without reservoir damage, and is of very great significance for optimal design of a declogging measure and improvement or restoration of oil well production and water well injection capacity for damaged wells, and improvement of numerical simulation precision of oil pools.
Step S4301: determining a volume concentration of the deposited fine particles based on the spatio-temporal evolution simulation equation established by the modeling method for reservoir damage by the fine particles within the reservoir.
For the spatio-temporal evolution simulation equation of reservoir damage by fine particle migration expressed by the above formula (4-14), the volume concentration C({right arrow over (r)}, t) of the deposited fine particles can be determined by referring to the process of solving a volume concentration of deposited particles in the above Embodiment 1, which will not be described here.
After the volume concentration C({right arrow over (r)}, t) of the migrating fine particles is calculated by the above method, the volume concentration Cd({right arrow over (r)}, t) of the deposited fine particles can be calculated according to the above formula (4-15), and thus the spatio-temporal evolution simulation equation established by the above modeling method for reservoir damage by the fine particles comprehensively considers the influence of various physical and chemical factors on reservoir damage during fine particle migration, so the volume concentration of the deposited fine particles obtained by the step S4301 is very precise.
Step S4302: determining a characteristic parameter characterizing the damage extent of the reservoir in a preset region of a well to be diagnosed, based on the volume concentration of the deposited fine particles.
Wherein the characteristic parameter may be permeability of the reservoir and/or a fluid loss coefficient of the reservoir.
In an embodiment, the characteristic parameter may be the permeability of the reservoir.
For the step S4302, the determining a characteristic parameter characterizing the damage extent of the reservoir in a preset region of a well to be diagnosed may include: determining the permeability K({right arrow over (r)}, t) of the reservoir based on the volume concentration Cd({right arrow over (r)}, t) of the deposited fine particles and formula (1-15).
In an embodiment, the characteristic parameter may be the fluid loss coefficient of the reservoir.
For the step S4302, the determining a characteristic parameter characterizing the damage extent of the reservoir in a preset region of a well to be diagnosed may include: determining the fluid loss coefficient k({right arrow over (r)}, t) of the reservoir based on the volume concentration Cd({right arrow over (r)}, t) of the deposited fine particles and formula (1-16),
The characteristic parameter may be a skin factor of the reservoir.
For the step S4302, determining the characteristic parameter characterizing the damage extent of the reservoir in the preset region of the well to be diagnosed may include: determining the permeability K({right arrow over (r)}, t) of the reservoir based on the volume concentration Cd({right arrow over (r)}, t) of the deposited fine particles and formula
and determining the skin factor S({right arrow over (r)}, t) of the reservoir based on the permeability K({right arrow over (r)}, t) of the reservoir and formula (1-17),
The characteristic parameter obtained by the step S4302 (e.g., the permeability K({right arrow over (r)}, t) and the skin factor S({right arrow over (r)}, t) of the reservoir) is a result of 4D quantitative simulation of spatio-temporal evolution (as shown in
where Kmax({right arrow over (r)}, t) is a maximum value of K({right arrow over (r)}, t)), a working person concerned can visually confirm the damage extent of the reservoir from
In summary, the volume concentration of the deposited fine particles can be determined by using the determined spatio-temporal evolution simulation equation, and then the characteristic parameter (e.g., the permeability and/or the skin factor of the reservoir) characterizing the damage extent of the reservoir in the preset region of the well to be diagnosed can be determined according to the volume concentration of the deposited fine particles, whereby a four-dimensional spatio-temporal evolution process of characteristics of reservoir damage caused by the fine particles within the reservoir can be simulated quantitatively. Therefore, performing quantitative prediction of reservoir damage and spatio-temporal deduction of damage laws is of scientific guidance significance for preventing or avoiding reservoir damage, and formulating development plans for oil pools and subsequent well stimulation measures for a well without reservoir damage, and is of very great significance for optimal design of a declogging measure and improvement or restoration of oil well production and water well injection capacity for damaged wells, and improvement of numerical simulation precision of oil pools.
A water lock effect mainly occurs in water-wettable reservoir rock. As water is a wetting phase, water always first occupies small pores and then medium and large pores to split oil gas into a dispersed phase, thus significantly reducing the permeation of oil gas in the reservoir (e.g., reducing the permeability of the reservoir). The water lock effect is controlled by a variety of factors such as reservoir lithology, physical properties, pore structures, and an invading fluid. Especially, geometrical characteristics of a reservoir medium have a great influence on reservoir damage by the water lock effect. Different pore throat structure distribution modes and complexity can lead to significant changes in the distribution mode of the water wetting phase in the rock, thus affecting the permeability of the reservoir.
Therefore, the core of the embodiments of the present invention is to establish a kinetic model (i.e., a diffusion equation for the diffusion of the water molecules through a solid-liquid interface from a liquid phase in the pores to the interior of a solid phase and a convection diffusion equation for the fluid in the pores) of diffusion of water molecules into the interior of the rock and water content variations within the pores in the reservoir. Specifically, based on Fick's law of diffusion, a convection diffusion relationship of the fluid in the pores in the reservoir, and the like, a spatio-temporal evolution control phenomenological model (containing a water volume fraction c1 of the pores in the reservoir and an initial value c0 of a water volume fraction of the rock in the reservoir) of porosity distribution in the reservoir around the well to be diagnosed influenced by the water lock effect is established, and in conjunction with a relationship between porosity of the reservoir and a characteristic parameter characterizing the damage extent of the reservoir such as permeability, spatio-temporal field distribution of the characteristic parameter such as permeability can be diagnosed.
Step S5101: determining a Darcy apparent velocity of a fluid in a reservoir in a preset region of a well to be diagnosed.
Wherein the well to be diagnosed may be, for example, a water injection well.
For the step S5101, determining a velocity of a fluid in a reservoir may include: establishing a pressure conduction equation for the fluid entering the reservoir; and determining the Darcy apparent velocity of the fluid according to the pressure conduction equation and a Darcy formula.
For the specific determination process, reference can be made to the process of determining the Darcy apparent velocity in the above Embodiment 2 (i.e., the above formulas (2-1) and (2-2) and related description thereof).
Step S5102: establishing an aqueous phase motion equation of the reservoir according to the Darcy apparent velocity of the fluid and a diffusion coefficient of water molecules in the fluid.
Under reservoir conditions, water contents at different locations within pores in the reservoir satisfy a mass conservation equation. Wherein the motion of an extraneous aqueous phase within the reservoir is mainly determined by two processes: convection and diffusion. Specifically, for the step S5102, the establishing an aqueous phase motion equation of the reservoir may include: establishing a mass balance equation expressed in the following formula according to the Darcy apparent velocity u of the fluid and the diffusion coefficient Dw of the water molecules:
where ϕ0 is an initial value of porosity of the reservoir; ϕw({right arrow over (r)}, t) is absolute porosity with pores in the reservoir being occupied by the aqueous phase; and {right arrow over (r)} is a spatial location of any point in the reservoir (e.g., using the center of the well to be diagnosed as an origin).
The aqueous phase motion equation expressed by the following formula (5-1) is established according to the mass balance equation and a spatio-temporal distribution function
of an aqueous phase saturation of the reservoir:
An initial condition for the aqueous phase motion equation is Sw({right arrow over (r)}, t=0)=Swc, and a boundary condition for the aqueous phase motion equation is Sw(|{right arrow over (r)}|=rw, t)=1 (that is, reservoir pores at a well wall of the water injection well is completely filled with water, i.e., the aqueous phase saturation in the pores is 1), where ϕ0 is the initial value of the porosity of the reservoir; rw is a wellbore radius of the well to be diagnosed; and Swc is an irreducible water saturation in the reservoir.
Step S5103: establishing a permeability distribution equation of the reservoir according to a pore size distribution characteristic of the pores in the reservoir and a preset permeability model of the reservoir.
For the step S5103, as shown in
Step S5201: determining a volume density function of pores with a pore size λ and a pore size distribution equation of the aqueous phase saturation of the reservoir according to the pore size distribution characteristic of the pores in the reservoir.
To quantitatively describe the pore structures of the reservoir, a fractal theory is used to study the water lock effect in the pore structures. According to the geometric principle of fractal, if pore size distribution of the reservoir has a fractal characteristic, the number N (>λ) of pores with a pore size larger than λ in the reservoir has the following power function relationship with λ:
where D is a fractal dimension of the pores (0<D<3).
In the case where the pore size distribution characteristic of the pores in the reservoir is that the number N(>λ) of the pores with the pore size larger than λ in the reservoir satisfies the above formula (5-2), as shown in
Step S5301: determining a total volume of the pores in the reservoir to be Φmax=A(λmax3-D−λmin3-D) according to the number N(>λ) of the pores with the pore size larger than λ in the reservoir.
Specifically, the total number N(>λmin) of the pores in the reservoir can then be obtained according to the above formula (5-2):
the following formula (5-4) can be obtained according to formulas (5-2) and (5-3):
then a relationship between the number N(>λ) of the pores with the pore size greater than λ and λ is a power function relationship expressed by the following formula (5-5):
N(>λ)=∫λλmaxƒ(λ)dλ=aλ31 D, (5-5)
λ, λmin and λmax in the formula (5-5) are a pore size, a minimum pore size and a maximum pore size of the pores, respectively (λmin and λmax can be obtained from an average pore size and a standard deviation of the pore size distribution; generally
and a is a proportional constant.
Next, from formula (5-5), a pore size distribution density function ƒ(λ) of the reservoir can be obtained, which satisfies the following formula (5-6):
in the formula, α′=−Da is a proportional constant.
A fractal expression of the total volume of the pores in the reservoir can be obtained from the pore size distribution density function expressed by the above formula (5-6):
Φmax=∫λ
where α is a constant related to the shape of the pores (α=1 if the shape of the pores is cube, or α=π/6 if the shape of the pores is sphere), and by integration, we can obtain:
Φmax=A(λmax3-D−λmin3-D), (5-8)
similarly, the volume of the pores with the pore size smaller than λ in the reservoir is Φλ=∫λ
Step S5302: according to the total volume Φmax of the pores in the reservoir and the volume Φλ=A(λ3-D−λmin3-D) of the pores with the pore size smaller than λ in the reservoir, determining the volume density function of the pores with the pore size λ as:
where D is the fractal dimension of the pores; and λ, λmin and λmax are the pore size, minimum pore size and maximum pore size of the pores, respectively; and A=αα′/(3-D) (a constant).
Step S5202: establishing the permeability distribution equation of the reservoir according to the preset permeability model, the volume density function of the pores with the pore size λ and the pore size distribution equation of the aqueous phase saturation.
As shown in
Step S5401: determining a volume of pores occupied by a non-aqueous phase to be Φnw(λ)=A(λmax3-D−λ3-D) according to the number N(>λ) of the pores with the pore size larger than λ.
Assuming that the pores with the pore size smaller than λ are completely occupied by the aqueous phase, and the pores with the pore size larger than λ are completely occupied by the non-aqueous phase (i.e., the rock in the reservoir is water-wettable (i.e., hydrophilic)), the volume Φnw(λ) of the pores occupied by the non-aqueous phase can be obtained in conjunction with the above formula (5-8),
Φnw(λ)=A(λmax3-D−λ3-D). (5-10)
Step S5402: determining the pore size distribution equation
of the aqueous phase saturation expressed by the following formula according to the total volume Φmax of the pores in the reservoir and the volume Φnw(λ) of the pores occupied by the non-aqueous phase.
Wherein D is the fractal dimension of the pores; and λ, λmin and λmax are the pore size, minimum pore size and maximum pore size of the pores, respectively; and A=αα′/(3-D).
Specifically, a pore size distribution equation of a non-aqueous phase saturation can be determined according to formula (5-8) and formula (5-10) to be
and then the pore size distribution equation of the aqueous phase saturation can be determined from the pore size distribution equation Snw(λ) of the non-aqueous phase saturation,
According to a linear Hagen-Poiseuille viscous flow, permeability of a capillary bundle model can be expressed as
In the embodiments of the present invention, a permeable channel of the reservoir can be regarded as an accumulation of multiple capillary bundles. Due to the continuity of the pore size distribution, the expression of the permeability of the capillary bundle model can be written in integral form as:
In the case where the preset permeability model of the reservoir satisfies
establishing the permeability distribution equation of the reservoir may include: establishing a permeability distribution equation (not shown) of the reservoir according to the preset permeability model
of the reservoir, the volume density function dξ of the pores with the pore size λ and the pore size distribution equation of the aqueous phase saturation.
According to the established permeability distribution equation, a distribution equation of a permeability damage rate expressed by the following formula can be further established,
Specifically, first, substituting formula (5-9) into formula (5-12) can yield a pore size distribution function of the permeability of the reservoir:
and then, replacing a variable λ in the above pore size distribution function K(λ) of the permeability by Sw according to a relational expression between λ and the water saturation Sw(λ) in formula (5-11) can yield formula (5-13).
Step S5104: determining the spatio-temporal evolution simulation equation of reservoir damage by the water lock effect according to the permeability distribution equation and the aqueous phase motion equation, wherein the spatio-temporal evolution simulation equation is used to simulate a four-dimensional spatio-temporal evolution process of characteristics of reservoir damage caused by the water lock effect.
Specifically, the spatio-temporal distribution function Sw({right arrow over (r)}, t) of the aqueous phase saturation of the reservoir can be obtained according to formula (5-1), and Sw({right arrow over (r)}, t) is substituted into a four-dimensional spatio-temporal distribution form of the permeability of the reservoir, i.e., the spatio-temporal evolution simulation equation of reservoir damage by the water lock effect is obtained.
In summary, according to the present invention, the Darcy apparent velocity of the fluid in the reservoir in the preset region of the well to be diagnosed is creatively determined; the aqueous phase motion equation of the reservoir is established according to the Darcy apparent velocity of the fluid and the diffusion coefficient of the water molecules in the fluid; the permeability distribution equation of the reservoir is established; and the spatio-temporal evolution simulation equation of reservoir damage by the water lock effect is determined according to the permeability distribution equation and the aqueous phase motion equation. Thus, by using the determined spatio-temporal evolution simulation equation, a four-dimensional spatio-temporal evolution process of characteristics of reservoir damage caused by the water lock effect can be quantitatively simulated. Therefore, performing quantitative prediction of reservoir damage and spatio-temporal deduction of damage laws is of scientific guidance significance for preventing or avoiding reservoir damage, and formulating development plans for oil pools and subsequent well stimulation measures for a well without reservoir damage, and is of very great significance for optimal design of a declogging measure and improvement or restoration of oil well production and water well injection capacity for damaged wells, and improvement of numerical simulation precision of oil pools.
Correspondingly, another embodiment of the present invention further provides a method for determining a damage extent of a reservoir. The method may include: determining a characteristic parameter characterizing the damage extent of the reservoir in a preset region of a well to be diagnosed, based on the spatio-temporal evolution simulation equation established by the modeling method for reservoir damage by the water lock effect.
For the solution of the spatio-temporal evolution simulation equation for reservoir damage by the water lock effect described above, Sw({right arrow over (r)}, t) needs to be calculated according to formula (5-1). For the specific solving process, reference can be made to the solving process of the volume concentration of the deposited particles in the above Embodiment 1, which will not be described here.
After the aqueous phase saturation Sw({right arrow over (r)}, t) of the reservoir is calculated by the above method, the permeability K({right arrow over (r)}, t) of the reservoir can be calculated according to the above formula (5-13) (of course, variations of the permeability K({right arrow over (r)}, t) at a particular location r in the reservoir with time may also be obtained, as shown in
The characteristic parameter characterizing the damage extent of the reservoir in the preset region of the well to be diagnosed can be calculated based on the permeability of the reservoir.
In an embodiment, the characteristic parameter may be a permeability damage rate of the reservoir.
Correspondingly, the determining a characteristic parameter characterizing the damage extent of the reservoir in a preset region of a well to be diagnosed may include: determining the permeability K({right arrow over (r)}, t) of the reservoir based on the spatio-temporal evolution simulation equation; and determining the permeability damage rate I({right arrow over (r)}, t) of the reservoir based on the permeability K({right arrow over (r)}, t) of the reservoir and formula (5-14):
where Kmax ({right arrow over (r)}, t) is a maximum value of K ({right arrow over (r)}, t).
In another embodiment, the characteristic parameter may be a skin factor of the reservoir.
The determining a characteristic parameter characterizing the damage extent of the reservoir in a preset region of a well to be diagnosed may include: determining the permeability K({right arrow over (r)}, t) of the reservoir based on the spatio-temporal evolution simulation equation; and determining the skin factor S({right arrow over (r)}, t) of the reservoir based on the permeability K({right arrow over (r)},t) of the reservoir and formula (5-15):
where Ko ({right arrow over (r)}) is an initial value of the permeability of the reservoir; and
The characteristic parameters obtained by the above embodiments (e.g., the permeability K({right arrow over (r)}, t) and the skin factor S({right arrow over (r)}, t) of the reservoir) are a result of 4D quantitative simulation of spatio-temporal evolution (
The distribution of pore structures and solid matrices (i.e., rock matrix or framework particles) within a reservoir is very complex. A microscopic solid matrix (i.e., a cluster) forms a seepage channel, and multiple pores (equivalent to hydraulic bundles) are formed between multiple solid matrices (i.e., clusters), and the multiple solid matrices and the multiple pores form a bundle of capillary tubes (also called a capillary bundle) as shown in
Moreover, the reservoir is subjected to a downward pressure due to the action of an overlying rock layer of the reservoir on the one hand, and the reservoir is subjected to an upward pressure due to the support of a fluid within the reservoir on the other hand. As a result, under the action of the two pressures (i.e., effective stress) on the reservoir, the capillary bundle shrinks (as shown in
Step S6201: determining an effective stress on a reservoir in a preset region of a well to be diagnosed.
For the step S6201, the determining an effective stress on a reservoir in a preset region of a well to be diagnosed may include: establishing a pressure conduction equation for a fluid entering the reservoir; and determining the effective stress on the reservoir according to the pressure conduction equation and a pressure from an overlying rock layer of the reservoir.
Specifically, a pressure is a force that drives a solid-liquid mixture to continuously intrude from a wellbore of a water injection well into the reservoir around the well to be diagnosed, whereby the pressure conduction equation for the fluid entering the reservoir as expressed in formula (6-1) can be established:
An initial condition for the pressure conduction equation is P({right arrow over (r)}, t=0)=P0, and a boundary condition for the pressure conduction equation is P(|{right arrow over (r)}|=rw, t)=Pw (that is, the pressure at a well wall of the water injection well is Pw).
Then, a pressure P({right arrow over (r)}, t) (i.e., a pore pressure) generated by the fluid on the reservoir can be obtained according to formula (6-1), and an effective stress σ({right arrow over (r)}, t) on the reservoir can then be determined according to P({right arrow over (r)}, t) and a pressure PA({right arrow over (r)}, t) from the overlying rock layer:
σ({right arrow over (r)}, t)=PA({right arrow over (r)}, t)−P({right arrow over (r)}, t) (6-2)
Step S6202: determining a flow rate of the fluid in the reservoir according to pore size distribution characteristics of pores in the reservoir, a diameter and length of each capillary bundle of the reservoir under the effective stress, and a fluid flow formula.
Wherein the capillary bundle is composed of a plurality of solid matrices and pores between the plurality of solid matrices.
For the step S6202, as shown in
Before performing step S6301, the determining a flow rate of the fluid in the reservoir may further include: determining the diameter of each capillary bundle of the reservoir under the effective stress according to an elastic modulus of the reservoir; and determining the length of each capillary bundle of the reservoir under the effective stress according to the elastic modulus and a Poisson's ratio of the reservoir.
Processes of determining the diameter and length of each capillary bundle of the reservoir under the effective stress are described below respectively.
A transverse wave velocity of the well to be diagnosed is a key feature for determining important data such as the Poisson's ratio and elastic modulus of the reservoir. During field logging, shear wave logging is generally not performed for the well to be diagnosed, but rather the Poisson's ratio is determined from a silt content and a linear formula about the silt content. Wherein the linear formula is fitted by particular rock samples, so the precision of the Poisson's ratio determined by the linear formula is very low for different reservoirs. In the embodiment, the transverse wave velocity of the well to be diagnosed can be indirectly calculated according to acquired transverse wave velocities and longitudinal wave velocities of a plurality of adjacent wells around the well to be diagnosed, and a longitudinal wave velocity acquired by field logging of the well to be diagnosed.
Before performing the step of determining the length of each capillary bundle of the reservoir under the effective stress, the determining a flow rate of the fluid in the reservoir may further include: determining a transverse wave velocity of the reservoir according to a transverse wave velocity and a longitudinal wave velocity of each of a plurality of particular adjacent wells located in the preset region; and determining the elastic modulus and the Poisson's ratio of the reservoir according to the transverse wave velocity of the reservoir.
Wherein the particular adjacent wells may be wells for which transverse and longitudinal wave velocities have been acquired by means of logging.
Specifically, the transverse wave velocity νt of the well to be diagnosed may be determined by a longitudinal wave velocity νl and a relational expression νt=ανl−b. The coefficients a and b in the relational expression may be calculated based on the transverse and longitudinal wave velocities of the adjacent wells. For example, the coefficients a and b are obtained by linear regression using a least square method:
where n is the number of adjacent wells in the preset area for which the transverse and wave velocities are obtained by logging; νti is a transverse wave velocity of an ith adjacent well; νli is a longitudinal wave velocity of the ith adjacent well;
Compared with an existing method of obtaining a Poisson's ratio described above, the result of the transverse wave velocity of the well to be diagnosed obtained from logging data such as the transverse and longitudinal wave velocities of the adjacent wells in the embodiment is very precise because the logging data is very precise, and thus the Poisson's ratio and other data obtained subsequently are also very accurate. Next, how to acquire the Poisson's ratio and the elastic modulus of the reservoir by using the transverse wave velocity of the well to be diagnosed will be introduced below.
Specifically, the Poisson's ratio v and the elastic modulus E of the reservoir can be obtained from the following formula (6-3) and the transverse wave velocity of the well to be diagnosed:
where ρ is a volume density of the reservoir; Δtt is a transverse wave time difference (which can be measured by a single-transmitter dual-receiver sound velocity measurement device, i.e., a time difference between reception, of transverse waves of glide waves formed by sound waves, by two receivers spaced by a distance L in a vertical direction); and Δt is a sound wave time difference (which can be measured by a single-transmitter dual-receiver sound velocity measurement device, i.e., a time difference between reception, of glide waves formed by sound waves, by two receivers spaced by a distance L in the vertical direction).
Based on the acquired elastic modulus and Poisson's ratio, the diameter and length of any one capillary bundle in the reservoir under a certain stress are determined below.
First, the determining the diameter of each capillary bundle of the reservoir under the effective stress may include: determining the diameter λp of the particular capillary bundle of the reservoir under the effective stress σ({right arrow over (r)}, t) according to the elastic modulus E of the reservoir and the following equation:
λp=F(1+σ({right arrow over (r)},t)/E)λc0,
where F is a shape factor of the particular capillary bundle; and λc0 is an initial diameter of the particular capillary bundle of the reservoir under no effective stress.
Second, the determining the length of each capillary bundle of the reservoir under the effective stress may include: determining an initial length Lc0 of the particular capillary bundle of the reservoir under no effective stress according to an apparent length of the particular capillary bundle of the reservoir under no effective stress, the initial diameter of the particular capillary bundle of the reservoir under no effective stress, and a fractal scale law; and determining a length Lp of the particular capillary bundle of the reservoir under the effective stress σ({right arrow over (r)}, t) according to the elastic modulus E of the reservoir, the Poisson's ratio v of the reservoir, the initial length Lc0 of the particular capillary bundle of the reservoir under no effective stress, and the following formula: Lp=[1−σ({right arrow over (r)}, t)/(vE)]Lc0.
Wherein the determining the initial length Lc0 of the particular capillary bundle of the reservoir under no effective stress may include: determining the initial length Lc0 of the particular capillary bundle of the reservoir under no effective stress according to the apparent length L0 of the particular capillary bundle of the reservoir under no effective stress, the initial diameter λc0 of the particular capillary bundle of the reservoir under no effective stress, and the fractal scale law Lc0=λc01-D
Step S6301: determining a flow ratio of the fluid according to the diameter and length of each capillary bundle of the reservoir under the effective stress, and a fluid flow ratio formula.
The flow ratio q(λc0, σ({right arrow over (r)},t) of the fluid is determined according to the obtained diameter λp and length Lp of the particular capillary bundle of the reservoir under the effective stress σ({right arrow over (r)}, t) described above, and the fluid flow ratio formula
wherein G is a geometric factor of a transverse section of the particular capillary bundle; Δp is a displacement pressure at both ends of the transverse section of the reservoir; and F is a shape factor of the particular capillary bundle.
Step S6302: determining a flow rate of the fluid according to the pore size distribution characteristics of pores of the reservoir, the flow ratio of the fluid in the reservoir, and the fluid flow formula.
For the step S6302, the determining a flow rate of the fluid includes: determining the flow rate Q(σ({right arrow over (r)}, t)) of the fluid according to the pore size distribution characteristics
of pores of the reservoir, the flow ratio q(λc0, σ({right arrow over (r)}, t)) of the fluid in the reservoir, and the fluid flow formula Q=∫λ
Specifically, based on the theory of rock fractal geometry, number density Nc of pores of a solid matrix of with a diameter λc0 (i.e., number cumulative distribution of solid matrices of a certain diameter per unit volume, or pore size distribution characteristics of the pores) satisfies the relation:
For a two-dimensional (reservoir) model, 0<Dcf<2; and for a three-dimensional (reservoir) model, 0<Dcf<3.
According to the pore size distribution characteristics of the pores in the reservoir described above, dNc=Dcfλc0maxD
Step S6203: determining a spatio-temporal evolution simulation equation of reservoir damage by stress sensitivity according to a permeability model of the reservoir and the flow rate of the fluid in the reservoir.
Wherein the spatio-temporal evolution simulation equation is used to simulate a four-dimensional spatio-temporal evolution process of characteristics of reservoir damage caused by stress sensitivity.
For the step S6203, the determining a spatio-temporal evolution simulation equation of reservoir damage by stress sensitivity may include: determining the spatio-temporal evolution simulation equation of reservoir damage by the stress sensitivity according to the permeability model
of the reservoir and the flow rate Q(σ({right arrow over (r)}, t)) of the fluid in the reservoir of the reservoir. Where μ is the viscosity of the fluid; L0 is the apparent length of the particular capillary bundle of the reservoir under no effective stress; A is an area of a transverse section of the reservoir; and Δp is the displacement pressure at both ends of the transverse section of the reservoir. Specifically, A=L02=∫λ
Specifically, the spatio-temporal evolution simulation equation of reservoir damage by the stress sensitivity expressed by the following formula (6-6) can be obtained according to the permeability model
of the reservoir and the flow rate Q(σ({right arrow over (r)}, t)) expressed by the above formula (6-5),
In summary, according to the present invention, the effective stress on the reservoir in the preset region of the well to be diagnosed is creatively determined; the flow rate of the fluid in the reservoir is determined according to the pore size distribution characteristics of the pores of the reservoir, the diameter and length of each capillary bundle of the reservoir under the effective stress, and the fluid flow formula; and the spatio-temporal evolution simulation equation of reservoir damage by the stress sensitivity is determined according to the permeability model of the reservoir and the flow rate of the fluid in the reservoir. Thus, by using the determined spatio-temporal evolution simulation equation, a four-dimensional spatio-temporal evolution process of the characteristics of reservoir damage caused by the stress sensitivity can be quantitatively simulated. Therefore, performing quantitative prediction of reservoir damage and spatio-temporal deduction of damage laws is of scientific guidance significance for preventing or avoiding reservoir damage, and formulating development plans for oil pools and subsequent well stimulation measures for a well without reservoir damage, and is of very great significance for optimal design of a declogging measure and improvement or restoration of oil well production and water well injection capacity for damaged wells, and improvement of numerical simulation precision of oil pools.
Correspondingly, another embodiment of the present invention further provides a method for determining a damage extent of a reservoir. The method includes: determining a characteristic parameter characterizing the damage extent of the reservoir in a preset region of a well to be diagnosed, based on the spatio-temporal evolution simulation equation by the modeling method for reservoir damage by stress sensitivity.
For the solution of the spatio-temporal evolution simulation equation for reservoir damage by pressure sensitivity expressed by the above formula (6-6), P({right arrow over (r)}, t) needs to be calculated according to formula (6-1). For the specific solving process, reference can be made to the solving process of the volume concentration of the deposited particles in the above Embodiment 1, which will not be described here. Then an aqueous phase saturation Sw({right arrow over (r)}, t) of the reservoir is determined according to the P({right arrow over (r)}, t) obtained by calculation.
After P({right arrow over (r)}, t) and the aqueous phase saturation Sw({right arrow over (r)}, t) of the reservoir are calculated by the above method, the permeability K({right arrow over (r)}, t) of the reservoir can be calculated according to the above formula (6-6), and thus the spatio-temporal evolution simulation equation established by the above modeling method for reservoir damage by pressure sensitivity comprehensively considers the influence of various physical and chemical factors on reservoir damage during damage by pressure sensitivity, so the permeability of the reservoir obtained by the embodiment is very precise.
A characteristic parameter characterizing the damage extent of the reservoir in the preset region of the well to be diagnosed can be calculated based on the permeability of the reservoir.
In an embodiment, the characteristic parameter may be a fluid loss coefficient of the reservoir.
Correspondingly, the determining a characteristic parameter characterizing the damage extent of the reservoir in a preset region of a well to be diagnosed may include: determining permeability K({right arrow over (r)}, t) of the reservoir based on the spatio-temporal evolution simulation equation; and determining a permeability damage rate I({right arrow over (r)}, t) of the reservoir based on the permeability K({right arrow over (r)}, t) of the reservoir and formula (5-14).
In another embodiment, the characteristic parameter may be a skin factor of the reservoir.
The determining a characteristic parameter characterizing the damage extent of the reservoir in a preset region of a well to be diagnosed may include: determining the permeability K({right arrow over (r)}, t) of the reservoir based on the spatio-temporal evolution simulation equation; and determining the skin factor S of the reservoir based on the permeability K({right arrow over (r)}, t) of the reservoir and formula (5-15).
The characteristic parameters (e.g., the permeability K({right arrow over (r)}, t) the skin factor S({right arrow over (r)}, t) and the permeability damage rate of the reservoir) obtained by the above embodiments are a result of 4D quantitative simulation of spatio-temporal evolution (as shown in
According to a sand production mechanism, key parameters of reservoir damage by sand production are a threshold flow velocity and a critical flow velocity, so the core of quantitative simulation of spatio-temporal evolution of reservoir damage by sand production is to solve the threshold flow velocity and the critical flow velocity. The threshold flow velocity is a flow velocity for incipient motion of sand grains. When a fluid flow velocity is greater than the threshold flow velocity, part of the sand grains start to move, and sand produced at the moment is attached sand; and if the fluid flow velocity exceeds the critical flow velocity, shear failure occurs on a rock framework, a large amount of sand starts to be produced, and the sand produced at the moment includes attached sand and framework sand. When the fluid flow velocity exceeds the threshold flow velocity and is lower than the critical flow velocity, the reservoir starts to produce sand partly, which has almost no influence on the reservoir (e.g., permeability), and proper sand production is helpful to the permeability instead; and when the fluid flow velocity is greater than the critical flow velocity such that the reservoir produces a large amount of sand, relatively great damage is caused to the reservoir (e.g., permeability), so only the situation where the fluid flow velocity is greater than the critical flow velocity such that the reservoir is influenced is considered in the embodiments of the present invention. The critical flow velocity is related to a critical production differential pressure (or critical production), so a specific process of determining the critical production differential pressure (or critical production) is involved herein. On this basis, based on mass conservation, a diffusion relationship, and the like, a spatio-temporal evolution control phenomenological model (containing a concentration C of sand grains and a concentration Cd of deposited sand grains) of concentration distribution of the sand grains in a reservoir around a well to be diagnosed is established, and in conjunction with a relationship between a deposition concentration and a characteristic parameter characterizing the damage extent of the reservoir such as permeability, spatio-temporal field distribution of the characteristic parameter such as permeability can be diagnosed.
Step S7101: determining a velocity of a fluid in a reservoir.
Wherein the reservoir is located in a preset region of a well to be diagnosed (e.g., a water injection well).
For the step S7101, the determining a velocity of a fluid in a reservoir may include: establishing a pressure conduction equation for the fluid entering the reservoir; and determining the velocity of the fluid according to the pressure conduction equation and a Darcy formula.
For the specific determination process, reference can be made to the process of determining a velocity of a fluid in the above Embodiment 1 (i.e., the above formulas (1-1) and (1-2) and related description thereof).
Step S7102: establishing a mass balance equation between the fluid and deposited sand grains on rock in the reservoir, based on a convection parameter and a diffusion parameter of the fluid and a mass change rate of sand grains in the fluid.
There is a correlation between the mass change rate of the sand grains and a crude oil production of the reservoir. Specifically, the correlation between the mass change rate of the sand grains and the crude oil production of the reservoir includes: the mass change rate of the sand grains is greater than 0 in the case where the crude oil production of the reservoir is greater than critical production.
In an embodiment, the critical production is obtained by: determining a pressure of the fluid according to the pressure conduction equation; determining a critical flowing bottom hole pressure of the fluid at the time the reservoir starts to produce sand according to an effective radial stress and an effective circumferential stress of the reservoir and a Mohr-Coulomb criterion; and determining the critical production according to the critical flowing bottom hole pressure of the fluid, the pressure of the fluid and a Dupuit formula.
In establishment of a sand production damage model, first a critical velocity of the fluid when the sand grains incipiently move and migrate is considered, and then how the sand grains change a solid-liquid flow deposition equation is considered.
When a stratum is drilled, stress distribution around the well changes and stress concentration occurs at a well wall.
Under the combined effect of a downhole liquid column pressure (a first term on the right side of the following formula), a crustal stress (second and third terms on the right side of the following formula) and fluid seepage (a fourth term on the right side of the following formula), a radial stress σr ({right arrow over (r)}, t) and a circumferential stress σθ({right arrow over (r)}, t) of the reservoir (e.g., the rock of the reservoir is an isotropic uniform elastomer) are respectively:
where rw is a wellbore radius of the well to be diagnosed; P ({right arrow over (r)}, t) is the pressure of the fluid; σH and σh are a maximum horizontal crustal stress and a minimum horizontal crustal stress, respectively; ν({right arrow over (r)}, t) is a flow velocity of the fluid; ϕ is porosity of the reservoir; β is a pore Biot coefficient
where Cr is a rock basement compression coefficient; and Cb is a rock volume compression coefficient); θ is an included angle between radial and horizontal maximum crustal stress directions at {right arrow over (r)} in the reservoir; and δ is 1 when permeation occurs at the well wall, and δ is 0 when permeation does not occur at the well wall. A situation where permeation occurs at the well wall (i.e., δ is 1) is considered here.
For a porous permeable reservoir, a relationship between the radial stress and an effective radial stress σ′r; of the reservoir satisfies the following formula (7-3), and a relationship between the circumferential stress and an effective circumferential stress σ′θ of the reservoir satisfies the following formula (7-4):
σθ=σ′θ+βP({right arrow over (r)}, t); (7-3)
σr=σ′r+βP({right arrow over (r)}, t). (7-4)
A tectonic stress is neglected, and according to an Anderson uniaxial strain model (σH=σh), an oilfield crustal stress is:
where σv is a pressure from an overlying rock layer of the reservoir. According to density logging data, the overlying rock layer pressure σv can be obtained by the following formula:
If a theoretical value of an overlying rock layer gradient is adopted, σv=22.7 H, where H is a depth; and if the overlying rock layer pressure is assumed to increase uniformly with the depth, σv=[ρS(1−ϕ)+ρLϕ]gH, where ρs is the average density of the rock framework; ρL is the density of the fluid; and H is the depth.
According to the above formulas (7-1)-(7-5), it can be determined that the effective radial stress and an effective axial stress of the reservoir satisfy the following equations, respectively:
When r=rw and the depth is the depth where a bottom hole is located, P({right arrow over (r)}, t)=Pwf. As a differential pressure is the largest at the well wall, only after the reservoir at the well wall produces sand, can other locations produce sand. In the embodiment, the effective radial stress and the effective circumferential stress on the surface of the well wall are considered as:
Under the combined effect of the crustal stress, pore pressures in the stratum and friction during fluid seepage on stratum rock, the effective radial stress increases and exceeds a yield condition in an unconsolidated sandstone reservoir, which will cause instability and plastic flow of the reservoir rock to produce sand, and reservoir rock failure follows the Mohr-Coulomb criterion. Specifically, when a maximum principal stress σmax and a minimum principal stress σmin are used, fluid pressures in the pores of the reservoir are considered, the Mohr-Coulomb criterion can be expressed as:
where C is a sand-mudstone cohesive force; φ is a stratum internal frictional angle; νpo is a Poisson's ratio of the rock; νp is a longitudinal wave velocity; and νcl is a shale content (%).
For the sand-mudstone cohesive force C, the sand-mudstone cohesive force CSo and a sound wave propagation velocity νp satisfy the following relationship:
where ρb(H)ρb(H) is rock volume density in the reservoir at the depth H. The shale content Vmcl can be calculated according to spontaneous potential logging data and using an empirical formula:
where PSP is a spontaneous potential of argillaceous sandstone; and SSP is a static spontaneous potential of water-bearing clean sandstone in the preset region. The stratum internal frictional angle φ can be calculated by the following equation:
ϕ=2.564 log10[M+(M2+1)1/2]+20, where M=58.93−1.785C.
σmax=σ′θσmin=σ′r, and,
and equations (7-8)-(7-9) are substituted into the above formula (7-10) to obtain the critical flowing bottom hole pressure at the time when the well to be diagnosed starts to produce sand:
and hence, according to the above formula (7-11), a critical production differential pressure ΔPcr can be obtained:
ΔPcr=P({right arrow over (r)}, 0)−Pcr. (7-12)
After the critical production differential pressure ΔPcr is obtained, the critical production Qcr can be determined from the Dupuit formula (7-13):
where Bo is an oil phase volume factor; re an oil pool radius; rw is the wellbore radius of the well to be diagnosed; μo is crude oil viscosity; and K is the permeability of the reservoir.
According to formula (7-13), it can be known that the critical production Qcr is closely related to the velocity of the fluid. Only when the actual velocity of the fluid in the reservoir exceeds the critical velocity (or the actual production of crude oil exceeds the critical production), can the reservoir produce sand under the action of the fluid.
According to a mass equation, assuming that the mass change rate of the sand grains (i.e., the amount of sand grains released) is qs, qs has the following properties:
In other words, only when the production of crude oil exceeds the critical production, can the reservoir produce sand and the sand grains formed by sand production enter the fluid to participate in migration, thereby increasing the mass of the fluid-solid mixture. Therefore, for step S7102, establishing a mass balance equation between the fluid and the deposited sand grains on rock in the reservoir may include: establishing the mass balance equation expressed in the following formula, based on a convection parameter and a diffusion parameter of the fluid,
where ρ is the density of the fluid; ϕ is the porosity of the reservoir; w ({right arrow over (r)}, t) is the mass fraction (which may also be called a mass concentration) of the deposited sand grains; u is a Darcy apparent velocity; j ({right arrow over (r)}, t) is a diffusion flow rate, j({right arrow over (r)}, t)=−ϕρLD∇w({right arrow over (r)}, t), where ρL is the density of the fluid (i.e. ρL=ρ), D({right arrow over (r)}, t) is a diffusion coefficient of the sand grains, D({right arrow over (r)}, t)=αν({right arrow over (r)}, t), α is a vertical diffusivity, and ν({right arrow over (r)}, t) is the velocity of the fluid;
{dot over (m)}({right arrow over (r)}, t) is an accumulated mass of the deposited sand grains per unit time; t is time; and qs is the mass change rate of the sand grains.
The mass change rate qs of the sand grains is obtained by: determining the intensity Q(r) of a release field of the deposited sand grains; determining a decay function Y(t) of the intensity of the release field; and determining the mass change rate qs=Q(r)Y(t) of the sand grains according to the intensity Q(r) of the release field and the decay function Y(t) of the intensity of the release field. Specifically, the intensity Q(r) of the release field may be a constant (q0), and the decay function Y(t) may be an exponential decay function (e.g. e−λt, where λ is a decay constant) that can vary with time.
Step S7103: establishing a connection condition equation between a volume concentration of the deposited sand grains and a volume concentration of the fluid, based on the convection parameter and the diffusion parameter of the fluid.
For the step S7103, the establishing a connection condition equation between a volume concentration of the deposited sand grains and a volume concentration of the fluid may include: establishing the connection condition equation expressed in the following formula (7-15), based on the convection parameter and the diffusion parameter of the fluid,
where ρp is the density of the deposited sand grains; Cd is the volume concentration of the deposited sand grains; and k({right arrow over (r)}, t)=k0({right arrow over (r)})Gl(Cd)Fl(T), where k0 is an original fluid loss coefficient,
Since the correlation between F1(T) and temperature is measured by exp(1/T) and in a common temperature range (e.g., 300 K to 400 K), the change of this function is actually very slow and actually close to an isothermal process, thus
where Cd({right arrow over (r)}, t) is the volume concentration of the deposited sand grains, Cd max is a maximum volume concentration of the deposited sand grains, and mk is a first empirical value. All of the above parameters can be either constants, or parameters that vary with space, i.e., in a non-homogeneous situation.
Step S7104: determining a spatio-temporal evolution simulation equation of reservoir damage by the sand production according to a relationship between the mass fraction of the sand grains and the volume concentration of the sand grains, the velocity of the fluid, the mass balance equation and the connection condition equation.
Wherein the relationship between the mass fraction of the migrating sand grains and the volume concentration of the migrating sand grains may be
where ρp is the density of the deposited sand grains; ρL is the density of the fluid; w({right arrow over (r)}, t) is the mass fraction of the sand grains; and C({right arrow over (r)}, t) is the volume concentration of the sand grains. The spatio-temporal evolution simulation equation of reservoir damage by the sand production may include: a spatio-temporal evolution simulation equation of reservoir damage by the sand production expressed by formula (7-16), and a spatio-temporal evolution simulation equation of reservoir damage by the sand grain deposition expressed by formula (7-17).
For the step S7104, the determining a spatio-temporal evolution simulation equation of reservoir damage by the sand production may include: determining the spatio-temporal evolution simulation equation of reservoir damage by the sand production expressed by formula (7-16) according to the relationship between the mass fraction of the sand grains and the volume concentration of the sand grains, the velocity of the fluid, and the mass balance equation expressed by formula (7-14):
and determining the spatio-temporal evolution simulation equation of reservoir damage by the sand grain deposition expressed by formula (7-17) according to the relationship between the mass fraction of the sand grains and the volume concentration of the sand grains, the velocity of the fluid, and the connection condition equation expressed by formula (7-15):
where C({right arrow over (r)}, t) where the volume concentration of the sand grains; ν({right arrow over (r)}, t) is the velocity of the fluid; τ is the tortuosity of the reservoir; ρp is the density of the deposited sand grains; ρL is the density of the fluid;
and k0({right arrow over (r)}) is an initial value of the fluid loss coefficient of the reservoir; Cd({right arrow over (r)}, t) is the volume concentration of the deposited sand grains; Cd max is the maximum volume concentration of the deposited sand grains; mk is a first empirical value; α is vertical diffusivity; ϕ is porosity of the reservoir; and qs is a mass change rate of the sand grains.
k0({right arrow over (r)})=ƒ(NR, NPe, NA, NDL, NE1, NE2, NG, NLo, NvdW, ζp(g)), where NR, NPe, NA, NDL, NE1, NE2, NG, NLo, NvdW, ζp(g) are a radius number, a Peclet number, an attraction number, an electrical double layer number, a first electric potential force number, a second electric potential force number, a gravity number, a London force number, a van der Waals force number, and potentials of sand grains and matrix particles (i.e., particles deposited on the rock), respectively (for details of relevant expressions of the parameters, see Table 2).
In summary, according to the present invention, the mass balance equation between the fluid and the deposited sand grains on the rock in the reservoir is creatively established based on the convection parameter and the diffusion parameter of the fluid and the mass change rate of sand grains in the fluid, wherein there is a correlation between the mass change rate of the sand grains and a crude oil production of the reservoir; the connection condition equation between the volume concentration of the deposited sand grains and the volume concentration of the fluid is established based on the convection parameter and the diffusion parameter of the fluid; and the spatio-temporal evolution simulation equation of reservoir damage by the sand production is determined according to the relationship between the mass fraction of the sand grains and the volume concentration of the sand grains, the velocity of the fluid, the mass balance equation and the connection condition equation. Thus, by using the determined spatio-temporal evolution simulation equation, a four-dimensional spatio-temporal evolution process of characteristics of reservoir damage caused by the sand production can be quantitatively simulated. Therefore, performing quantitative prediction of reservoir damage and spatio-temporal deduction of damage laws is of scientific guidance significance for preventing or avoiding reservoir damage, and formulating development plans for oil pools and subsequent well stimulation measures for a well without reservoir damage, and is of very great significance for optimal design of a declogging measure and improvement or restoration of oil well production and water well injection capacity for damaged wells, and improvement of numerical simulation precision of oil pools.
Step S7201: determining the volume concentration of the deposited sand grains based on the spatio-temporal evolution simulation equation established by the modeling method for reservoir damage by the sand production.
For the spatio-temporal evolution simulation equation of reservoir damage by the sand production expressed by the above formula (7-16), the volume concentration C({right arrow over (r)}, t) of the sand grains can be calculated by referring to the process of solving a volume concentration of deposited particles in the above Embodiment 1.
After the volume concentration C({right arrow over (r)}, t) of the sand grains is calculated by the above method, the volume concentration Cd({right arrow over (r)}, t) of the deposited sand grains can be calculated according to the above formula (7-17), and thus the spatio-temporal evolution simulation equation established by the above modeling method of reservoir damage by the sand grains comprehensively considers the influence of various physical and chemical factors on reservoir damage during sand production, so the volume concentration of the deposited sand grains obtained by the step S7201 is very precise.
Step S7202: determining a characteristic parameter characterizing the damage extent of the reservoir in a preset region of a well to be diagnosed, based on the volume concentration of the deposited sand grains.
In an embodiment, the characteristic parameter may be permeability of the reservoir.
For the step S7202, the determining a characteristic parameter characterizing the damage extent of the reservoir in a preset region of a well to be diagnosed may include: determining the permeability K({right arrow over (r)}, t) of the reservoir based on the volume concentration Cd({right arrow over (r)}, t) of the deposited sand grains and formula (1-15).
In an embodiment, the characteristic parameter may be a fluid loss coefficient of the reservoir.
For the step S7202, the determining a characteristic parameter characterizing the damage extent of the reservoir in a preset region of a well to be diagnosed may include: determining the permeability K({right arrow over (r)}, t) of the reservoir based on the volume concentration Cd({right arrow over (r)}, t) of the deposited sand grains and formula (1-15); and determining the fluid loss coefficient k({right arrow over (r)}, t) of the reservoir based on the volume concentration Cd({right arrow over (r)}, t) of the deposited sand grains and formula (1-16).
Wherein the characteristic parameter may be a skin factor of the reservoir.
For the step S7202, the determining a characteristic parameter characterizing the damage extent of the reservoir in a preset region of a well to be diagnosed may include: determining the permeability K({right arrow over (r)}, t) of the reservoir based on the volume concentration of the deposited sand grains and formula
and determining the skin factor S({right arrow over (r)}, t) of the reservoir based on the permeability K({right arrow over (r)}, t) of the reservoir and formula (1-17).
The characteristic parameter (e.g., the permeability K({right arrow over (r)}, t) and the skin factor S({right arrow over (r)}, t) of the reservoir) obtained by the step S7202 is a result of 4D quantitative simulation of spatio-temporal evolution (as shown in
where Kmax({right arrow over (r)}, t) is a maximum value of K({right arrow over (r)}, t)), and a working person concerned can visually confirm the damage extent of the reservoir from
In summary, the volume concentration of the deposited sand grains can be determined by using the determined spatio-temporal evolution simulation equation, and then the characteristic parameter (e.g., the permeability and/or the skin factor of the reservoir) characterizing the damage extent of the reservoir in the preset region of the well to be diagnosed can be determined according to the volume concentration of the deposited sand grains, whereby a four-dimensional spatio-temporal evolution process of the characteristics of reservoir damage caused by the sand production can be simulated quantitatively. Therefore, performing quantitative prediction of reservoir damage and spatio-temporal deduction of damage laws is of scientific guidance significance for preventing or avoiding reservoir damage, and formulating development plans for oil pools and subsequent well stimulation measures for a well without reservoir damage, and is of very great significance for optimal design of a declogging measure and improvement or restoration of oil well production and water well injection capacity for damaged wells, and improvement of numerical simulation precision of oil pools.
Wettability reversal is a phenomenon in which surfaces of pores in a reservoir change from hydrophilic to lipophilic, which weakens circulation of an oil phase in the pores, resulting in poor permeability of the reservoir. When an oil-phase saturation is high, the oil phase occupies large pores and is good in connectivity, and its flow has characteristics similar to a capillary flow, and the permeability is in a linear relation with the square of the saturation; and when the oil-phase saturation is low, the oil phase is mainly dispersed and attached to wall surfaces of small pores, and is poor in connectivity, and as the saturation decreases, the permeability of the oil phase decreases at a faster rate (e.g., changing based on the 4th power law).
On the one hand, a wettability damage extent of the reservoir is determined by a relationship between relative permeability of the oil phase and the oil (or water) saturation; on the other hand, pressure distribution in the reservoir also influences a fluid flow velocity and permeability in the pores. Therefore, the core of the embodiments of the present invention is to establish a pressure distribution field of the oil phase (considering pressure distribution of the oil phase and an aqueous phase in the fluid, respectively) and a convection diffusion law of the oil phase. Specifically, the pressure distribution field of the oil phase is determined according to a pressure distribution equation of the reservoir in a preset region around a well to be diagnosed and a force balance condition of capillaries; and spatio-temporal field distribution of a characteristic parameter characterizing the damage extent of the reservoir such as permeability can be diagnosed in conjunction with a Darcy formula and the convection diffusion law of the oil phase.
Step S8101: determining a relationship between a pressure distribution field of an aqueous phase and a pressure distribution field of capillaries in a reservoir according to a pressure distribution equation of the reservoir in a preset region of a well to be diagnosed.
Wherein the capillaries are formed by the wettability reversal of a contact surface between the aqueous phase and an oil phase in the reservoir.
Usually, a pressure field of a liquid (mixture of an oil phase and an aqueous phase) is used to describe a pressure distribution field of the reservoir as a whole, but in the embodiment, the oil phase is separated from the aqueous phase, and pressure distribution of the reservoir is considered by using pressure distribution of the oil phase and pressure distribution of the aqueous phase, respectively, so that the pressure distribution of the reservoir can be simulated more closely to the actual situation in the reservoir, and thus a spatio-temporal evolution simulation equation of reservoir damage by wettability reversal can be simulated more precisely (i.e., a very precise reservoir permeability result is obtained) by means of the pressure distribution field of the oil phase.
For the step S8101, the relationship between the pressure distribution field Pw({right arrow over (r)}, t) of the aqueous phase and the pressure distribution field Pc({right arrow over (r)}, t) of the capillaries in the reservoir is determined according to the pressure distribution equation of the reservoir in the preset region of the well to be diagnosed expressed in the following formula (8-1):
where ϕ is porosity (a constant) of the reservoir; ct is an integrated compression coefficient (a constant) of the reservoir; k, krw and kro are permeability of the reservoir, relative permeability of the oil phase and relative permeability of the aqueous phase of the reservoir; and μw and μo are the viscosity of the aqueous phase and the viscosity of the oil phase of the reservoir.
Step S8102: determining a pressure distribution field of the oil phase according to the relationship between the pressure distribution field of the aqueous phase and the pressure distribution field of the capillaries and a force balance condition of the capillaries.
The force balance condition of the capillaries may be a three-force balance condition expressed by the following formula (8-2),
Pc({right arrow over (r)}, t)=Po({right arrow over (r)}, t)−Pw({right arrow over (r)}, t), (8-2)
where Pc({right arrow over (r)}, t) is a pressure of the capillaries and Pc({right arrow over (r)}, t) is determined by an effective water saturation in the capillaries; Po({right arrow over (r)}, t) is the pressure distribution field of the oil phase; and Pw({right arrow over (r)}, t) is the pressure distribution field of the aqueous phase.
Determining Pc({right arrow over (r)}, t) by an effective water saturation in the capillaries may include: determining Pc({right arrow over (r)}, t) according to the effective water saturation and the following formula (8-3):
where S*w is the effective water saturation in the capillaries; m is a Corey constant; and Pce is a pressure threshold of the capillaries.
Specifically, the effective water saturation S*w may be determined by: determining, according to a saturation of the oil phase and a saturation of the aqueous phase, the effective water saturation expressed by the following formula (8-4):
S*w({right arrow over (r)}, t)=[(Sw({right arrow over (r)}, t)−Swir)/(1−So({right arrow over (r)}, t)−Swir)], (8-4)
wherein So({right arrow over (r)}, t) is the saturation of the oil phase; Sw({right arrow over (r)}, t) is the saturation of the aqueous phase, and So({right arrow over (r)}, t)+Sw({right arrow over (r)}, t)=1; and Swir is an irreducible water saturation in the capillaries.
In other words, the pressure distribution field Po({right arrow over (r)}, t) of the oil phase can be obtained according to the above formulas (8-1)-(8-4).
Step S8103: determining a velocity distribution field of the oil phase according to the pressure distribution field of the oil phase and a Darcy formula.
Specifically, the velocity distribution field of the oil phase can be determined by substituting the pressure distribution field Po({right arrow over (r)}, t) of the oil phase into the Darcy formula expressed by the following formula (8-5):
where μo is viscosity of the oil phase; and K({right arrow over (r)}, t) is permeability of the reservoir.
Step S8104: determining a spatio-temporal evolution simulation equation of reservoir damage by wettability reversal according to a convection diffusion law of the oil phase, the velocity distribution field and a dispersion coefficient of the oil phase.
Wherein the spatio-temporal evolution simulation equation is used to simulate a four-dimensional spatio-temporal evolution process of characteristics of reservoir damage caused by wettability reversal.
For the step S8104, the determining a spatio-temporal evolution simulation equation of reservoir damage by wettability reversal may include: determining, according to the convection diffusion law of the oil phase, the velocity distribution field uo({right arrow over (r)}, t) and the dispersion coefficient Do of the oil phase, the spatio-temporal evolution simulation equation of reservoir damage by the wettability reversal expressed by the following formula:
where ϕ({right arrow over (r)}, t) is the porosity of the reservoir; and So({right arrow over (r)}, t) is the saturation of the oil phase.
That is, the saturation So({right arrow over (r)}, t) of the oil phase (i.e., the spatio-temporal evolution simulation equation of reservoir damage by the wettability reversal) can be determined according to formulas (8-1)-(8-6).
In summary, according to the present invention, the relationship between the pressure distribution field of the aqueous phase and the pressure distribution field of the capillaries in the reservoir is creatively determined according to a pressure distribution equation of the reservoir in the preset region around the well to be diagnosed; then the pressure distribution field of the oil phase is determined according to the relationship between the pressure distribution field of the aqueous phase and the pressure distribution field of the capillaries and the force balance condition of the capillaries; subsequently the velocity distribution field of the oil phase is determined according to the pressure distribution field of the oil phase and the Darcy formula; and finally the spatio-temporal evolution simulation equation of reservoir damage by wettability reversal is determined according to the convection diffusion law of the oil phase, the velocity distribution field and the dispersion coefficient of the oil phase. Thus, by using the determined spatio-temporal evolution simulation equation, a four-dimensional spatio-temporal evolution process of characteristics of reservoir damage caused by wettability reversal can be quantitatively simulated. Therefore, performing quantitative prediction of reservoir damage and spatio-temporal deduction of damage laws is of scientific guidance significance for preventing or avoiding reservoir damage, and formulating development plans for oil pools and subsequent well stimulation measures for a well without reservoir damage, and is of very great significance for optimal design of a declogging measure and improvement or restoration of oil well production and water well injection capacity for damaged wells, and improvement of numerical simulation precision of oil pools.
Step S8201: determining a saturation of an oil phase in a reservoir in a preset region of a well to be diagnosed based on the spatio-temporal evolution simulation equation established by the modeling method for reservoir damage by wettability reversal.
For the solution of the spatio-temporal evolution simulation equation of reservoir damage by wettability reversal expressed by the above formula (8-6), reference can be made to the process of solving the volume concentration of the deposited particles in the above Embodiment 1, so that the saturation So({right arrow over (r)}, t) of the oil phase can be calculated.
The saturation So({right arrow over (r)}, t) of the oil phase can be calculated by the above method. As the spatio-temporal evolution simulation equation established by the above modeling method for reservoir damage by wettability reversal comprehensively considers the influence of various physical and chemical factors on reservoir damage during wettability reversal, the saturation of the oil phase obtained by the step S8201 is very precise.
Step S8202: determining a characteristic parameter characterizing the damage extent of the reservoir based on the determined saturation of the oil phase.
For the step S8202, the characteristic parameter may be relative permeability of the oil phase. Correspondingly, the determining a characteristic parameter characterizing the damage extent of the reservoir may include: determining the relative permeability Kro({right arrow over (r)}, t) of the oil phase based on the saturation So({right arrow over (r)}, t) of the oil phase and a relationship between the relative permeability and the saturation of the oil phase expressed by the following formula (8-7):
Kro({right arrow over (r)}, t)=α0α1So({right arrow over (r)}, t)+αSSo2({right arrow over (r)}, t)+α3So3({right arrow over (r)}, t)+α4So4({right arrow over (r)}, t), (8-7)
where α1, α2, α3, α4, α5 are constants, and 0≤Kro({right arrow over (r)}, t)≤1.
Further, the permeability K({right arrow over (r)}, t) of the oil phase can be determined according to the relative permeability Kro({right arrow over (r)}, t).
In an embodiment, the characteristic parameter may be a permeability damage rate of the reservoir.
Correspondingly, the determining a characteristic parameter characterizing the damage extent of the reservoir may include: calculating the permeability damage rate I({right arrow over (r)}, t) of the reservoir based on the permeability K({right arrow over (r)}, t) of the reservoir and formula (8-8),
where Kmax({right arrow over (r)}, t) is a maximum value of K({right arrow over (r)}, t).
In another embodiment, the characteristic parameter may be a skin factor of the reservoir. Correspondingly, the determining a characteristic parameter characterizing the damage extent of the reservoir may include: calculating the skin factor S({right arrow over (r)}, t) of the reservoir based on the permeability K({right arrow over (r)}, t) of the reservoir and formula (8-9):
where Ko ({right arrow over (r)}) is an initial value of the permeability of the reservoir; and
The characteristic parameter (e.g., the permeability K({right arrow over (r)}, t) and the skin factor S({right arrow over (r)}, t) of the reservoir) obtained by the step S8202 is a result of 4D quantitative simulation of spatio-temporal evolution (as shown in
In summary, according to the present invention, the saturation of the oil phase can be creatively calculated by using the determined spatio-temporal evolution simulation equation, then the characteristic parameter (e.g., the permeability and/or the skin factor of the reservoir) characterizing the damage extent of the reservoir in the preset region of the well to be diagnosed can be determined based on the determined saturation of the oil phase, and thus, a four-dimensional spatio-temporal evolution process of characteristics of reservoir damage caused by wettability reversal can be quantitatively simulated. Therefore, performing quantitative prediction of reservoir damage and spatio-temporal deduction of damage laws is of scientific guidance significance for preventing or avoiding reservoir damage, and formulating development plans for oil pools and subsequent well stimulation measures for a well without reservoir damage, and is of very great significance for optimal design of a declogging measure and improvement or restoration of oil well production and water well injection capacity for damaged wells, and improvement of numerical simulation precision of oil pools.
In a porous medium in a reservoir, a low interfacial tension and a high mechanical shear between an oil phase and an aqueous phase are main factors that cause emulsion formation. Therefore, a radius of emulsified droplets formed is determined based on the influence of a temperature field of the reservoir on the viscosity of the oil phase and an emulsification condition of a fluid; then, a spatio-temporal evolution control phenomenological model of a clogging probability is determined based on a pore size distribution function of the reservoir and the radius of the emulsified droplets; and subsequently, in conjunction with a relationship between the clogging probability and a characteristic parameter characterizing the damage extent of the reservoir such as permeability, spatio-temporal field distribution of the characteristic parameter such as permeability can be diagnosed.
Step S9101: determining a Darcy apparent velocity of a fluid in a reservoir in a preset region of a well to be diagnosed.
Wherein the well to be diagnosed may be, for example, an oil production well.
For the step S9101, determining the velocity of the fluid in the reservoir may include:
establishing a pressure conduction equation for the fluid entering the reservoir; and determining the Darcy apparent velocity of the fluid according to the pressure conduction equation and a Darcy formula.
For the specific determination process, reference can be made to the process of determining the Darcy apparent velocity in the above Embodiment 2 (i.e., the above formulas (2-1) and (2-2) and related description thereof).
Step S9102: determining a viscosity of an oil phase in the reservoir according to a temperature field of the reservoir and a function relationship between the viscosity of the oil phase and a temperature.
Before performing the step S9102, the modeling method may further include: determining the temperature field of the reservoir according to a thermal conductivity coefficient of the fluid, a thermal diffusion coefficient of the reservoir, an average flow velocity of an extraneous fluid, and a heat balance equation of the reservoir.
Specifically, the temperature field T({right arrow over (r)}, t) of the reservoir is determined according to the thermal conductivity coefficient Dƒ of the fluid, the thermal diffusion coefficient Dl of the reservoir, the average flow velocity um of the extraneous fluid, and the heat balance equation of the reservoir expressed by the following formula (9-1):
In conjunction with an initial condition and a boundary condition, temperature distribution (i.e., the temperature field) of the reservoir under different spatio-temporal conditions can be calculated by this formula (9-1). See the following description for details of the specific calculation procedure.
For the step S9102, the viscosity μ0({right arrow over (r)}, t) of the oil phase is determined according to the temperature field T({right arrow over (r)}, t) of the reservoir and the function relationship between the viscosity and the temperature of the oil phase expressed by the following formula (9-2):
log[log(μ0({right arrow over (r)}, t)+1)]=a−bγAPI−c log(T({right arrow over (r)}, t)), (9-2)
where γAPI is a gravitational parameter of the oil phase; a and b are both constants; and {right arrow over (r)} is a spatial location of any point in the reservoir. Viscosity distribution of the oil phase under different spatio-temporal conditions of the reservoir can be calculated according to formulas (9-1) and (9-2). Since the viscosity distribution of the oil phase is greatly influenced by the temperature distribution of the reservoir, a calculation result of the viscosity of the oil phase obtained by considering the temperature field of the reservoir is more accurate, and reservoir damage by the emulsification clogging can be simulated more precisely according to the viscosity of the oil phase.
Step S9103: determining a radius of an emulsified droplet formed by an emulsification of the fluid according to the Darcy apparent velocity of the fluid, the viscosity of the oil phase, and an emulsification condition of the fluid.
Wherein the emulsification condition of the fluid may be a critical condition expressed by the following formula (9-3):
wherein μw is viscosity of the aqueous phase in the fluid; μo is viscosity of the oil phase; Kw is permeability of the aqueous phase; Ko is permeability of the oil phase; ρw is density of the aqueous phase; ρo is density of the oil phase; g is a gravitational acceleration; μ0({right arrow over (r)}, t) is viscosity of the oil phase; u({right arrow over (r)}, t) is a Darcy apparent velocity of the fluid; and r is a spatial location of any point in the reservoir.
That is, if inequality (9-3) is satisfied, it indicates that oil-water emulsification occurs in the reservoir.
For the step S9103, the determining a radius of an emulsified droplet formed by an emulsification of the fluid may include: determining, according to the Darcy apparent velocity u({right arrow over (r)}, iΔt) of the fluid, μ0({right arrow over (r)}, iΔt) the viscosity of the oil phase, and the emulsification condition of the fluid, the radius of the emulsified droplet expressed by the following formula (9-4):
(9-4) where σ is an oil-water interfacial tension;
Step S9104: determining a spatio-temporal evolution simulation equation of reservoir damage by emulsification clogging according to a pore size distribution function of pores of the reservoir and the radius of the emulsified droplet.
Wherein the spatio-temporal evolution simulation equation is used to simulate a four-dimensional spatio-temporal evolution process of reservoir damage characteristics caused by the emulsification clogging.
For the step S9104, the determining a spatio-temporal evolution simulation equation of reservoir damage by emulsification clogging may include: determining, according to the pore size distribution function N(λ, μs, σs) of the pores of the reservoir and the radius λo({right arrow over (r)}, iΔt) of the emulsified droplet, the spatio-temporal evolution simulation equation of reservoir damage by emulsification clogging expressed by the following formula (9-5):
(9-5) where β(λo({right arrow over (r)}, iΔt)) is a clogging probability of the reservoir; μs and σs are a first pore size distribution characteristic parameter and a second pore size distribution characteristic parameter, respectively; iΔt is the ith time increment, i being a non-negative integer; and {right arrow over (r)} is a spatial location of any point in the reservoir.
In an embodiment, the pore size distribution function N(λ) of the pores in the reservoir may be approximated as a log-normal function expressed by the following formula (9-6):
μs and σs in the above formula can be calculated by the following process: a pore size average value and standard deviation can be calculated as
(i.e.,
according to field data, and then expressions of μs and σs shown in the following formula (9-7) can be obtained according to the pore size average value and standard deviation:
Since the pore size average value and standard deviation are known quantities, the corresponding μs and σs, can be obtained, and then μs and σs are substituted into the above formula (9-6) to obtain a specific form of the pore size distribution function.
For an emulsified droplet with a radius λo, only part of the pores smaller than λo in the distribution function is clogged by the emulsified droplet, so the clogging probability β (i.e., the cumulative distribution from 0 to λo) in the above formula (9-5) can be specifically expressed as the following formula (9-8):
where erfc( ) is a complementary error function:
In summary, according to the present invention, the Darcy apparent velocity of the fluid in the reservoir in the preset region of the well to be diagnosed is creatively determined; the viscosity of the oil phase in the reservoir is determined according to the temperature field of the reservoir and the function relationship between the viscosity of the oil phase and the temperature; the radius of the emulsified droplet formed by an emulsification of the fluid is determined according to the emulsification condition of the fluid; and the spatio-temporal evolution simulation equation of reservoir damage by emulsification clogging is determined according to the pore size distribution function of the pores of the reservoir and the radius of the emulsified droplet. Thus, by using the determined spatio-temporal evolution simulation equation, a four-dimensional spatio-temporal evolution process of characteristics of reservoir damage caused by emulsification clogging can be quantitatively simulated. Therefore performing quantitative prediction of reservoir damage and spatio-temporal deduction of damage laws is of scientific guidance significance for preventing or avoiding reservoir damage, and formulating development plans for oil pools and subsequent well stimulation measures for a well without reservoir damage, and is of very great significance for optimal design of a declogging measure and improvement or restoration of oil well production and water well injection capacity for damaged wells, and improvement of numerical simulation precision of oil pools.
Step S9201: determining a clogging probability of the reservoir based on the spatio-temporal evolution simulation equation established by the modeling method for reservoir damage by emulsification clogging.
For the solution of the spatio-temporal evolution simulation equation for reservoir damage by emulsification clogging expressed by the above formula (9-8), T({right arrow over (r)}, t) needs to be calculated according to formula (9-1). For the specific solving process, reference can be made to the solving process of the volume concentration of the deposited particles in the above Embodiment 1, which will not be described here.
After the temperature field T({right arrow over (r)}, t) of the reservoir is calculated by the above method, the clogging probability of the reservoir can be calculated according to the above formulas (9-2), (9-4) and (9-8), and thus the spatio-temporal evolution simulation equation established by the above modeling method for reservoir damage by emulsification clogging comprehensively considers the influence of various physical and chemical factors on reservoir damage during emulsification clogging, so the clogging probability of the reservoir obtained by the step S9201 is very precise.
Step S9202: determining a characteristic parameter characterizing the damage extent of the reservoir in a preset region of a well to be diagnosed, based on the determined clogging probability of the reservoir.
Wherein the characteristic parameter may be permeability of the reservoir.
For the step S9202, the determining a characteristic parameter characterizing the damage extent of the reservoir in a preset region of a well to be diagnosed may include: determining the permeability Kd({right arrow over (r)}, t) of the reservoir based on the determined clogging probability β(λo({right arrow over (r)}, iΔt)) of the reservoir and the following formula (9-9):
where μs is a first pore size distribution characteristic parameter; λo({right arrow over (r)}, iΔt) is a radius of an emulsified droplet; mK is a second empirical value; and n is a total number of time increments Δt.
Wherein the characteristic parameter may be a skin factor of the reservoir.
For the step S9202, the determining a characteristic parameter characterizing the damage extent of the reservoir in a preset region of a well to be diagnosed may include: determining the permeability Kd({right arrow over (r)}, t) of the reservoir based on the determined clogging probability) β(λo({right arrow over (r)}, iΔt)) of the reservoir and the following formula:
and determining the skin factor S({right arrow over (r)}, t) of the reservoir based on the permeability Kd({right arrow over (r)}, t) of the reservoir and formula (9-10):
where μs is the first pore size distribution characteristic parameter; λo({right arrow over (r)}, iΔt) is the radius of the emulsified droplet; mK is the second empirical value; n is the total number of time increments Δt; rw is a wellbore radius of the well to be diagnosed, and rsw is a damage radius of the reservoir.
The characteristic parameter (e.g., the permeability Kd({right arrow over (r)}, t) and the skin factor S({right arrow over (r)}, t) of the reservoir) obtained by the step S9202 is a result of 4D quantitative simulation of spatio-temporal evolution (as shown in
where Kmax({right arrow over (r)}, t) is a maximum value of K({right arrow over (r)}, t)), and a working person concerned can visually confirm the damage extent of the reservoir from
In summary, according to the present invention, the clogging probability of the reservoir can be creatively calculated by using the determined spatio-temporal evolution simulation equation, then the characteristic parameter (e.g., the permeability and/or the skin factor of the reservoir) characterizing the damage extent of the reservoir in the preset region of the well to be diagnosed can be determined based on the determined clogging probability, and thus, a four-dimensional spatio-temporal evolution process of characteristics of reservoir damage caused by emulsification clogging can be quantitatively simulated. Therefore, performing quantitative prediction of reservoir damage and spatio-temporal deduction of damage laws is of scientific guidance significance for preventing or avoiding reservoir damage, and formulating development plans for oil pools and subsequent well stimulation measures for a well without reservoir damage, and is of very great significance for optimal design of a declogging measure and improvement or restoration of oil well production and water well injection capacity for damaged wells, and improvement of numerical simulation precision of oil pools.
During oilfield development, a pressure equilibrium state in an oil reservoir is destroyed, resulting in overflow of some light components and precipitation of some heavy components in crude oil, and the precipitation forms organic scale, which clogs oil gas flow channels to cause serious damage to the reservoir. Generally speaking, the organic scale such as asphaltene is partially dissolved and partially suspended as a colloid in the crude oil, and the colloidal precipitation of the organic scale such as asphaltene is related to its dissolving capacity. Thus, the core of the embodiments of the present invention is to establish a kinetic model of variations of the dissolving capacity of the organic scale with a reservoir pressure. Specifically, based on the influence of the pressure on the dissolving capacity of the organic scale, and the like, a spatio-temporal evolution control phenomenological model of organic scale particle distribution in a reservoir around a well to be diagnosed influenced by the organic scale is established, and in conjunction with a relationship between organic scale particle distribution and both the porosity of the reservoir and a characteristic parameter characterizing the damage extent of the reservoir (such as permeability), spatio-temporal field distribution of the characteristic parameter such as permeability can be diagnosed.
Step S10101: determining a pressure of a reservoir in a preset region of a well to be diagnosed.
First, a cylindrical model of the reservoir (or an oil pool) as shown in
Due to flow symmetry of a fluid, the pressure inside each cylindrical shell has equal magnitude (e.g., the pressure inside an ith cylindrical shell is P(ri)). According to a continuity equation of the pressure of the reservoir expressed by the following formula (10-1), an inner boundary condition P(rw)=Pw, (i.e., a horizontal distance from a wall of the well to be diagnosed to the reference central axis is rw, and the pressure at the wall of the well to be diagnosed is Pw) and an outer boundary condition P(re)=Pe (i.e., a horizontal distance from the cylindrical shell where the outermost part of the reservoir is located to the reference central axis is re, and the pressure at the cylindrical shell where the outermost part of the reservoir is located is Pe), expression (10-2) of P(ri) can be obtained,
In conjunction with the above formula (10-2) and a Dupuit formula expressed in the following formula (10-3), a reservoir pressure P(ri, t) at the ith cylindrical shell expressed by the following (10-4) can be obtained:
where Qo(t) is crude oil production of the well to be diagnosed; K is permeability of the reservoir; H is the thickness of the reservoir; and μ is viscosity of the fluid within the reservoir.
For the step S10101, the determining a pressure of a reservoir in a preset region of a well to be diagnosed may include: in the case where the reservoir is divided into a plurality of cylindrical shells with a central axis of the well to be diagnosed as a reference central axis and having a preset thickness, determining, according to the continuity equation of the pressure of the reservoir and the Dupuit formula, the pressure P(ri) of the reservoir expressed by the following formula (10-4):
where ri is an average horizontal distance from the ith cylindrical shell of the plurality of cylindrical shells to the reference central axis; re is an oil pool radius of the reservoir; Qo(t) is the crude oil production of the well to be diagnosed; K is permeability of the reservoir; H is the thickness of the reservoir; and μ is the viscosity of the fluid within the reservoir.
Step S10102: determining a first relational expression in which a maximum dissolved quantity of the organic scale in a crude oil produced from the reservoir varies with the pressure of the reservoir according to a bubble point pressure of the reservoir, a molar volume of the crude oil at the bubble point pressure, a solubility parameter of the crude oil, a solubility parameter of the organic scale in the crude oil, and a molar volume of the organic scale.
Wherein the solubility parameter of the crude oil is acquired by: determining solubility parameters of a plurality of preset components according to boiling point temperatures, critical temperatures and molar volumes of the plurality of preset components and the temperature of the reservoir, wherein the plurality of preset components are asphaltenes (e.g. C7, C8, and other C7+ asphaltenes) with a plurality of preset carbon contents; and determining the solubility parameter of the crude oil according to the solubility parameters and volume fractions of the plurality of preset components in the crude oil.
Specifically, the behavior of cohesive energy of crude oil molecules per unit volume (i.e., the solubility parameter SL of the crude oil) is the most complex, and the solubility parameters δi(P) of the components (e.g., asphaltenes with a plurality of preset carbon contents (e.g., C7, C8, and other C7+ asphaltenes)) need to be solved separately first. For the component i,
where Tci and Tbi are a critical temperature and a boiling point temperature of the component i, respectively; Vi is a molar volume of the component i (Vi(P)=xiV (P), where xi is a molar fraction of the component i; V(P) can be calculated from
where a and b are a first empirical coefficient and a second empirical coefficient, respectively); T is the temperature of the reservoir; and R is a gas constant.
Then, in conjunction with the above formula (10-5) and by using δL=Σinϕiδi(P), δL is calculated, where ϕi is the volume fraction (which can be obtained from analytical data of oil physical properties) of the component i; and n is the number of the components.
For the step S10102, the determining a first relational expression in which a maximum dissolved quantity of the organic scale in a crude oil produced from the reservoir varies with the pressure of the reservoir includes the following steps S10301-S10302, as shown in
Step S10301: determining a molar volume of the crude oil under the pressure of the reservoir according to the bubble point pressure of the reservoir and the molar volume of the crude oil produced from the reservoir at the bubble point pressure.
For the step S10301, the determining a molar volume of the crude oil under the pressure of the reservoir may include: determining, according to the bubble point pressure Pb of the reservoir and the molar volume VLb of the crude oil produced from the reservoir at the bubble point pressure, the molar volume VL(P(ri, t)) of the crude oil under the pressure P(ri, t) of the reservoir expressed by the following formula (10-6):
where Cƒ is a compression coefficient (e.g., (10˜440)×10−4 MPa−1) of the crude oil; and Bo is a compression coefficient (usually 1.0-1.2) of the crude oil.
When the pressure P is higher than the bubble point pressure, a pure compression process occurs, and the volume decreases as the pressure increases; and when the pressure P is lower than the bubble point pressure, on the one hand, the volume increases as the pressure decreases, and on the other hand, the crude oil decomposes into a gas phase and the volume decreases, forming an extreme value at a certain pressure point.
Step S10302: determining the first relational expression according to the molar volume of the crude oil under the pressure of the reservoir, the solubility parameter of the crude oil, the solubility parameter of the organic scale in the crude oil, and the molar volume of the organic scale.
For the step S10302, the determining the first relational expression may include: determining, according to the molar volume VL(P(ri, t)) of the crude oil under the pressure of the reservoir, the solubility parameter δL(P(ri, t)) of the crude oil, the solubility parameter δa of the organic scale in the crude oil, and the molar volume Va of the organic scale, the first relational expression expressed by the following formula (10-7):
where exp{ } is an exponential function with a natural constant e as its base, T is the temperature of the reservoir; R is the gas constant; and (ϕa)max may be in %. In an embodiment, the value of the parameter may be calculated by: δa=9.99×(1−5.94×10−4T).
Step S10103: determining a second relational expression in which a mole number of organic scale particles in the crude oil varies with both the pressure of the reservoir and the maximum dissolved quantity of the organic scale in the crude oil, according to a distribution function of organic scale particles in the organic scale and a mole number of the crude oil.
Wherein the distribution function is a proportional function of a mole number of organic scale particles with a particle size greater than a preset particle size to a total mole number of the organic scale particles.
For the step S10103, determining the second relational expression may include:
determining, according to the distribution function ƒtrap(Rp) of organic scale particles in the organic scale and the mole number ηo of the crude oil, the second relational expression in which the mole number η(P(ri, t), (ϕa)max(P(ri, t))) of organic scale particles in the crude oil varies with both the pressure P(ri, t) of the reservoir and the maximum dissolved quantity (ϕa)max(P(ri, t)) of the organic scale in the crude oil, which is expressed by the following formula (10-8):
η(P(ri, t), (ϕa)max(P(ri, t)))=∫0tηo[ϕa−(ϕa)max(P(ri, t))]ƒtrap(Rp)dt, (10-8)
where ηo is the mole number of the crude oil; ϕa is the total content of the organic scale in the crude oil; and Rp is the preset particle size.
The parameters in the above formula are explained and described below. For the mole number
of the crude oil, ρo is density of the crude oil; Qo(t) is the crude oil production of the well to be diagnosed; and MWo is an average molar mass of the crude oil. ηo[ϕa−(ϕa)max(P(ri, t))] represents a mole number change rate of organic scale particles in the ith cylindrical shell at time t; and for the distribution function ƒtrap(Rp), ƒtrap(Rp)=∫R
Step S10104: determining a spatio-temporal evolution simulation equation of reservoir damage by organic scale according to the first relational expression, the second relational expression and the pressure of the reservoir.
Wherein the spatio-temporal evolution simulation equation is used to simulate a four-dimensional spatio-temporal evolution process of reservoir damage characteristics caused by organic scale.
Specifically, formulas (10-4) and (10-8) are substituted into formula (10-7) to determine the spatio-temporal evolution simulation equation of reservoir damage by organic scale. As a result, the specific form of the spatio-temporal evolution simulation equation is very complex and is not illustrated here. That is, the spatio-temporal evolution simulation equation of reservoir damage by organic scale is equivalent to an equation set composed of formulas (10-4), and (10-7)-(10-8).
The embodiments of the present application mainly discuss the specific case of asphaltene, that is, the two concepts of organic scale and asphaltene are interchangeable.
In summary, according to the present invention, the first relational expression in which the maximum dissolved quantity of the organic scale in the crude oil varies with the pressure of the reservoir is determined creatively according to the bubble point pressure of the reservoir, the molar volume of crude oil produced from the reservoir at the bubble point pressure, the solubility parameter of the crude oil, the solubility parameter of the organic scale in the crude oil, and the molar volume of the organic scale; the second relational expression in which the mole number of the organic scale particles in the crude oil varies with both the pressure of the reservoir and the maximum dissolved quantity of the organic scale in the crude oil is determined according to the distribution function of organic scale particles in the organic scale and the mole number of the crude oil; and the spatio-temporal evolution simulation equation of reservoir damage by organic scale is determined according to the first relational expression, the second relational expression and the pressure of the reservoir. Thus, by using the determined spatio-temporal evolution simulation equation, a four-dimensional spatio-temporal evolution process of reservoir damage characteristics caused by organic scale can be quantitatively simulated. Therefore, performing quantitative prediction of reservoir damage and spatio-temporal deduction of damage laws is of scientific guidance significance for preventing or avoiding reservoir damage, and formulating development plans for oil pools and subsequent well stimulation measures for a well without reservoir damage, and is of very great significance for optimal design of a declogging measure and improvement or restoration of oil well production and water well injection capacity for damaged wells, and improvement of numerical simulation precision of oil pools.
Step S10401: determining a mole number of organic scale particles in the crude oil based on the spatio-temporal evolution simulation equation established by the modeling method for reservoir damage by organic scale.
For the pressure equation expressed by the above formula (10-1), the pressure P(ri, t) of the reservoir and the mole number η(P(ri, t), (ϕa)max(P (ri, t))) of the organic scale particles in the crude oil can be calculated by referring to the process of solving the volume concentration of the deposited particles in the above Embodiment 1.
The pressure P(ri, t) of the reservoir and the mole number η(P(ri, t), (ϕa)max(P(ri, t))) of the organic scale particles in the crude oil may be calculated by the above method, as the spatio-temporal evolution simulation equation established by the above modeling method for reservoir damage by organic scale comprehensively considers the influence of various physical and chemical factors on reservoir damage when the reservoir is clogged by organic scale particles, the mole number of the organic scale particles in the crude oil obtained by the step S10401 is very precise.
Step S10402: determining a characteristic parameter characterizing the damage extent of the reservoir in the preset region of the well to be diagnosed, based on the mole number of the organic scale particles in the crude oil.
In an embodiment, the characteristic parameter is a permeability of the reservoir.
For the step S10402, the determining a characteristic parameter characterizing the damage extent of the reservoir in the preset region of the well to be diagnosed may include: determining porosity ϕ(ri, t) of the reservoir based on the mole number η(P(ri, t), (ϕa)max(P(ri, t))) of the organic scale particles in the crude oil and a formula)
Where ϕ0 is an initial value of the porosity; mK is a second empirical value; and K0(ri) is an initial value of the permeability of the reservoir.
In an embodiment, the characteristic parameter is a fluid loss coefficient of the reservoir.
For the step S10402, the determining a characteristic parameter characterizing the damage extent of the reservoir in the preset region of the well to be diagnosed d may include: determining the porosity ϕ(ri, t) of the reservoir based on the mole number η(P(ri, t), (ϕa)max(P(ri, t))) of the organic scale particles in the crude oil and the formula
where dr is the preset thickness of the cylindrical shell; and determining the fluid loss coefficient k(ri, t) of the reservoir based on the porosity ϕ(ri, t) of the reservoir and the formula
Where ϕ0 is the initial value of the porosity; mk is a first empirical value; and k0(ri) is an initial value of the fluid loss coefficient of the reservoir.
In an embodiment, the characteristic parameter is a skin factor of the reservoir.
For the step S10402, the determining a characteristic parameter characterizing the damage extent of the reservoir in the preset region of the well to be diagnosed may include: determining the porosity ϕ(ri, t) of the reservoir based on the mole number η(P(ri, t), (ϕa)max(P(ri, t))) of the organic scale particles in the crude oil and the formula
determining the permeability K(ri, t) reservoir based on the porosity ϕ(ri, t) of the reservoir and the formula
and determining the skin factor S(ri, t) of the reservoir based on the permeability K(ri, t) of the reservoir and the formula
Where ϕ0 is the initial value of the porosity; mK is the second empirical value; K0(ri) is the initial value of the permeability of the reservoir;
The characteristic parameter (e.g., the permeability K(ri, t) and the skin factor S(ri, t) of the reservoir) obtained by the step S10402 is a result of 4D quantitative simulation of spatio-temporal evolution (as shown in
where Kmax(ri, t) is a maximum value of K(ri, t)), and a working person concerned can visually confirm the damage extent of the reservoir from
In summary, according to the present invention, the mole number of the organic scale particles in the crude oil is creatively determined based on the spatio-temporal evolution simulation equation established by the modeling method for reservoir damage by organic scale; and the characteristic parameter characterizing the damage extent of the reservoir in the preset region of the well to be diagnosed is determined based on the mole number of organic scale particles in the crude oil. Thus, a four-dimensional spatio-temporal evolution process of characteristics of reservoir damage caused by organic scale can be quantitatively simulated. Therefore, performing quantitative prediction of reservoir damage and spatio-temporal deduction of damage laws, which is of scientific guidance significance for preventing or avoiding reservoir damage, and formulating development plans for oil pools and subsequent well stimulation measures for a well without reservoir damage, and is of very great significance for optimal design of a declogging measure and improvement or restoration of oil well production and water well injection capacity for damaged wells, and improvement of numerical simulation precision of oil pools.
During flowing of an oil phase and an aqueous phase within a low-permeability reservoir, there are usually a large quantity of dispersed oil droplets and air bubbles. When these individuals move to narrow pore throats (for example, pores with a pore size of 2.5 μm) in the reservoir, as the diameters of these individuals are larger than the diameters of the pore throats, the flowing of these individuals is obstructed, forcing the oil droplets, air bubbles or the like to deform. Therefore, the phenomenon of additional capillary resistance on the oil droplets, bubbles or the like during movement in the reservoir with uneven diameter distribution is called a Jamin effect. Compared with the water lock effect, it is the oil droplets, air bubbles and the like that are forced to deform by the resistance in the small pore throats, and a wetting phase is water and a non-wetting phase is the oil droplets or bubbles in this case.
Usually, the Jamin effect occurs when the following conditions are satisfied at the same time: firstly, crude oil should be in a dispersed state and cannot move in the form of a continuous oil flow; secondly, diameters of channels of the geological reservoir should be relatively small. When a dispersed oil droplet moves into a narrow channel of the reservoir, as its diameter is larger than the diameter of the channel, the oil droplet cannot continue to move due to resistance, or when an oil droplet flows in a variable diameter channel, a capillary effect occurs, which leads to Jamin damage (a situation where permeability of the reservoir decreases due to the Jamin effect).
Specifically, the droplet moves from left to right under the action of a displacement pressure P on the left and right sides of the reservoir. When the oil droplet moves to a pore 1, due to the obstruction of the pore 1, a front end of the oil droplet (2R1=λ, where R1 is a radius of curvature of a liquid surface at the front end, and λ is a pore size of the pore 1) is subjected to a pre-bubble pressure P1; a rear end of the oil droplet (its diameter 2R2=
As a result, the additional resistance on the oil droplet due to the Jamin effect is:
and due to the front end
of the droplet, the back end
of the droplet, and
where σ is a surface tension of the liquid droplet. When the displacement pressure (i.e., a pressure difference outside the orifice throat) reaches at least Pc, the air bubble or droplet can pass through the pore throat. Otherwise, it is clogged.
The Jamin effect is influenced by various factors such as pore structures, lithology, physical properties and an invading fluid of the reservoir. Jiamin damage is closely related to geometric properties of a pore medium of rock. Different pore throat structure distribution modes and complexity can lead to significant changes in the distribution pattern of a water wetting phase in the rock, thus affecting the permeability of the reservoir. Therefore, the core of the embodiments of the present invention is to establish a kinetic model (i.e., an aqueous phase motion equation and a permeability distribution equation in the reservoir) of aqueous phase saturation variations within the pores in the reservoir. Specifically, the aqueous phase motion equation in the reservoir is established based on a convection diffusion relation (i.e. mass balance equation) for a fluid within the pores in the reservoir, and the like; the permeability distribution equation is established based on pore size distribution characteristics of the pores and a preset permeability model; then a spatio-temporal evolution control phenomenological model of permeability distribution in the reservoir around a well to be diagnosed influenced by the Jamin effect is determined according to the aqueous phase motion equation and the permeability distribution equation, and thus spatio-temporal field distribution of characteristic parameter characterizing the damage extent of the reservoirs such as permeability can be diagnosed.
Step S11101: determining a Darcy apparent velocity of a fluid in a reservoir in a preset region of a well to be diagnosed.
Wherein the permeability of the reservoir is lower than preset permeability; and the well to be diagnosed may be, for example, an oil production well.
For the step S11101, determining the velocity of the fluid in the reservoir may include: establishing a pressure conduction equation for the fluid entering the reservoir; and determining the Darcy apparent velocity of the fluid according to the pressure conduction equation and a Darcy formula.
For the specific determination process, reference can be made to the process of determining the Darcy apparent velocity in the above Embodiment 2 (i.e., the above formulas (2-1) and (2-2) and related description thereof).
Step S11102: establishing an aqueous phase motion equation of the reservoir according to the Darcy apparent velocity of the fluid and a diffusion coefficient of water molecules in the fluid.
Under reservoir conditions, water contents at different locations within pores in the reservoir satisfy a mass conservation equation. Motion of an extraneous aqueous phase within the reservoir is mainly determined by two processes: convection and diffusion. Specifically, for the step S11102, the establishing an aqueous phase motion equation of the reservoir may include: establishing, according to the Darcy apparent velocity u of the fluid and the diffusion coefficient Dw of the water molecules, a mass balance equation expressed in the following formula:
where ϕ0 is an initial value of porosity of the reservoir; ϕw({right arrow over (r)}, t) is absolute porosity with the pores in the reservoir being occupied by the aqueous phase; and {right arrow over (r)} is a spatial location of any point in the reservoir (e.g., using the center of the well to be diagnosed as an origin).
The aqueous phase motion equation expressed by the following formula (11-1) is established according to the mass balance equation and a spatio-temporal distribution function
of an aqueous phase saturation of the reservoir:
An initial condition for the aqueous phase motion equation is Sw({right arrow over (r)}, t=0)=Swc, and a boundary condition for the aqueous phase motion equation is Sw(|{right arrow over (r)}|=rw, t)=1 (that is, reservoir pores in a well wall of the water injection well are completely filled with water, i.e., the aqueous phase saturation in the pores is 1), where ϕ0 is an initial value of the porosity of the reservoir; rw is a wellbore radius of the well to be diagnosed; and Swc is an irreducible water saturation in the reservoir.
Step S11103: establishing a permeability distribution equation of the reservoir according to a pore size distribution characteristic of the pores of the reservoir and a preset permeability model of the reservoir.
For the step S11103, as shown in
Step S11201: determining a volume density function of pores with a pore size λ and a pore size distribution equation of the aqueous phase saturation of the reservoir according to the pore size distribution characteristic of the pores of the reservoir.
Pore structures of the rock of reservoir are strongly irregular and controlled by a variety of factors, such that pores in a pore structure model based on classical geometry are greatly different from real pores. To quantitatively describe the pore structures of the reservoir, a fractal theory is used to study the Jamin effect of two-phase flow in the pore structures, and a wetting film effect of the aqueous phase on the inner surfaces of pore channels is ignored. According to the geometric principle of fractal, if pore size distribution of the reservoir has a fractal characteristic, a number N (>λ) of pores with a pore size larger than λ in the reservoir has the following power function relationship with λ:
where D is a fractal dimension (2<D<3) of the pores.
In the case where the pore size distribution characteristic of the pores of the reservoir is that the number N(>λ) of the pores with the pore size larger than λ in the reservoir satisfies the above formula
as shown in
Step S11301: determining a total volume of the pores in the reservoir to be Φmax=A(λmax3-D−λmin3-D) according to the number N(>λ) of the pores with the pore size larger than λ in the reservoir.
Specifically, the total number N(>λmin) of the pores in the reservoir and the number λ of pores with the pore size larger than λ can be obtained according to the above formula
according to formulas (11-2) and (11-3), the following formula (11-4) can be obtained:
ƒ(λ) in the above formula (11-4) is a pore size distribution density function of the reservoir, and a relationship between the number N(<λpc) of the pores with the pore size smaller than λpc and λ is a power function relationship expressed by the following formula (11-5):
N(<λpc)=∫λλ
λ, λmin, λmax and λpc in the above related formula are a pore size, a minimum pore size, a maximum pore size of the pores (λmin and λmax can be obtained from an average pore size and a standard deviation of the pore size distribution; generally
and a maximum diameter with the Jamin effect, respectively (i.e., a particular pore size of a pore where a non-aqueous phase is subjected to minimum capillary resistance); and α is a proportional constant.
Next, from formula (11-5), the pore size distribution density function ƒ(λ) of the reservoir can be obtained, which satisfies the following formula (11-6):
in the formula, α′=−Da is a proportional constant.
A fractal expression of the total volume of the pores in the reservoir can be obtained from the pore size distribution density function expressed by the above formula (11-6):
Φmax=∫λminλmaxƒ(λ)α3dλ, (11-7)
where α is a constant related to the shape of the pores (α=1 if the shape of the pores is a cube, or α=π/6 if the shape of the pores is a sphere), and by integration, we can obtain:
ΦmaxA(λmax3-D−λmin3-D), (11-8)
similarly, the volume of the pores with the pore size smaller than λ in the reservoir is Φλ=∫λ
Step S11302: according to the total volume Φmax of the pores in the reservoir and the volume Φλ=A(λ3-D−λmin3-D) of the pores with the pore size smaller than λ in the reservoir, determining the volume density function of the pores with the pore size λ as:
wherein D is the fractal dimension of the pores; and λ, λmin and λmax are the pore size, minimum pore size and maximum pore size of the pores, respectively; and A=αα′/(3-D) (a constant).
Step S11202: establishing a permeability distribution equation of the reservoir according to the preset permeability model, the volume density function of the pores with the pore size λ and the pore size distribution equation of the aqueous phase saturation.
As shown in
Step S11401: determining a volume of pores occupied by a non-aqueous phase to be Φnw(λ)=A(λpc3-D−λ3-D) according to the number N(>λ) of the pores with the pore size larger than λ.
Wherein λpc is a particular pore size of a pore where the non-aqueous phase is subjected to minimum capillary resistance (i.e., a maximum pore throat diameter with the Jamin effect).
In the case where the reservoir is water-wettable (i.e., hydrophilic), a displacement pressure P (gradually increasing from zero to P) is applied to both sides of the reservoir, and oil droplets or air bubbles can smoothly pass through pores with a pore size in the range of λmax to λpc (i.e., communicating pores are opened successively from the maximum pore size λmax to λpc) according to the conditions for producing the Jamin effect. When the pressure P does not reach Pc, pore throats with λ larger than λpc are completely occupied by water, while pore throats with λ smaller than λpc are still clogged by the oil droplets or air bubbles.
Since the pore size distribution density function ƒ(λ) determined by the above formula (11-6) can be determined according to the number N(>λ) of pores with a pore size greater than λ in the reservoir, thus a volume Φnw(λ) of pores occupied by the non-aqueous phase when the pressure P does not reach Pc can be obtained according to ƒ(λ),
Φnw(λ)=A(λpc3-D−λ3-D). (11-10)
Step S11402: determining, according to the total volume Φmax of the pores in the reservoir and the volume Φnw(λ) of the pores occupied by the non-aqueous phase, the pore size distribution equation
of the aqueous phase saturation expressed by the following formula, where D is the fractal dimension of the pores; and λ, λmin and λmax are the pore size, minimum pore size and maximum pore size of the pores, respectively; and A=αα′/(3-D).
Specifically, the volume of the pores occupied by the aqueous phase can be determined according to formulas (11-8) and (11-10) to be Φw(λ)=Φmax−Φnw(λ)=A(λmax3-D−λmin3-D−λpc3-D+λ3-D); and the following formula can be obtained in conjunction with formula (11-8):
Since λmin<<λmax, and
the above equation can be written as the following formula (11-11), i.e., the pore size distribution equation (11-11) of the aqueous phase saturation can be determined:
According to the theory of fluid mechanics, under the effect of a differential pressure ΔP, a total flow rate of the fluid with viscosity μ passing through a capillary bundle can be described by a Hagen-poiseuille equation:
in the above formula, λi is a diameter of a pore channel i (i.e., a pore size of a pore i), Vi is a volume of the pore channel i, Vi=Φ0MLmξi, Φ0 is the porosity of the reservoir, ξi is a pore volume fraction of the pore channel i, Lc is a bending length of the pore channel i, Lm is a straight length of the pore channel i, M is an average pore throat cross-sectional area,
is tortuosity of the pore channel i (it may also be calculated by empirical formula
0.35). Substituting the above quantities into the above formula yields:
and according to a Darcy's law, permeability of a capillary bundle model can be expressed as
In the embodiments of the invention, a permeable channel of the reservoir can be regarded as an accumulation of multiple capillary bundles. Due to the continuity of the size distribution of the pores, the expression of the permeability of the capillary bundle model can be written in an integral form as:
In the case where the preset permeability model of the reservoir satisfies
the establishing apermeability distribution equation of the reservoir may include: establishing, according to the preset permeability model
of the reservoir, the volume density function dξ of the pores with the pore size λ and the pore size distribution equation of the aqueous phase, the permeability distribution equation of the reservoir:
Wherein an initial condition for the pressure conduction equation is K(t=0)=1. K(Sw) is a function of permeability with respect to water saturation (as shown in
(a larger D indicates a larger proportion of small pore throats; for a conventional reservoir λmax>20 μm, D<2.6; and for a low permeability reservoir λmax<0.2 μm, D>2.8), τ is the tortuosity, λmax is the maximum pore throat diameter, and λmin is a minimum pore throat diameter).
Specifically, first, the displacement pressure is gradually increased until P is achieved, at the moment, pore throats with a diameter larger than λ are completely occupied by water, wherein the absolute porosity with pores being occupied by the aqueous phase is ϕw, and formula (11-9) is substituted into formula (11-12) to obtain the pore size distribution function of the permeability of the reservoir:
and then, the variable λ in the above pore size distribution function K(λ) of the permeability is replaced by Sw according to the relational expression between λ and the water saturation Sw(λ) in formula (11-11), and in conjunction with ϕw=ϕ0Sw, to obtain formula (11-13) (i.e., obtaining an expression of the dimensionless permeability, which is a function of the water saturation).
Step S11104: determining a spatio-temporal evolution simulation equation of reservoir damage by the Jamin effect according to the permeability distribution equation and the aqueous phase motion equation, wherein the spatio-temporal evolution simulation equation is used to simulate a four-dimensional spatio-temporal evolution process of characteristics of reservoir damage caused by the Jamin effect.
Specifically, the spatio-temporal distribution function Sw({right arrow over (r)}, t) of the aqueous phase saturation of the reservoir (as shown in
In summary, according to the present invention, the Darcy apparent velocity of the fluid in the reservoir in the preset region in the well to be diagnosed is creatively determined, wherein the permeability of the reservoir is lower than preset permeability; the aqueous phase motion equation of the reservoir is established according to the Darcy apparent velocity of the fluid and the diffusion coefficient of water molecules in the fluid; the permeability distribution equation of the reservoir is established; and the spatio-temporal evolution simulation equation of reservoir damage by the Jamin effect is determined according to the permeability distribution equation and the aqueous phase motion equation. Thus, by using the determined spatio-temporal evolution simulation equation, a four-dimensional spatio-temporal evolution process of reservoir damage characteristics caused by the Jamin effect can be quantitatively simulated. Therefore, performing quantitative prediction of reservoir damage and spatio-temporal deduction of damage laws is of scientific guidance significance for preventing or avoiding reservoir damage, and formulating development plans for oil pools and subsequent well stimulation measures for a well without reservoir damage, and is of very great significance for optimal design of a declogging measure and improvement or restoration of oil well production and water well injection capacity for damaged wells, and improvement of numerical simulation precision of oil pools.
Correspondingly, another embodiment of the present invention further provides a method for determining a damage extent of a reservoir. The method may include: determining a characteristic parameter characterizing the damage extent of the reservoir in the preset region of the well to be diagnosed, based on the spatio-temporal evolution simulation equation by the modeling method for reservoir damage by the Jamin effect.
For the solution of the spatio-temporal evolution simulation equation for reservoir damage by the Jamin effect expressed by the above formula (11-13), Sw({right arrow over (r)}, t) needs to be calculated according to formula (11-1). For the specific solving process, reference can be made to the solving process of the volume concentration of the deposited particles in the above Embodiment 1, which will not be described here.
After the aqueous phase saturation Sw({right arrow over (r)}, t) of the reservoir is calculated by the above method, the permeability K({right arrow over (r)}, t) of the reservoir can be calculated according to the above formula (11-13), and thus the spatio-temporal evolution simulation equation established by the above modeling method of reservoir damage by the Jamin effect comprehensively considers the influence of various physical and chemical factors on reservoir damage during damage by the Jamin effect, so the permeability of the reservoir obtained by the embodiment is very precise.
A characteristic parameter characterizing the damage extent of the reservoir in the preset region of the well to be diagnosed is calculated based on the permeability of the reservoir.
In an embodiment, the characteristic parameter may be a permeability damage rate of the reservoir.
Correspondingly, the determining a characteristic parameter characterizing the damage extent of the reservoir in the preset region of the well to be diagnosed may include: determining the permeability K({right arrow over (r)}, t) of the reservoir based on the spatio-temporal evolution simulation equation; and determining the permeability damage rate I({right arrow over (r)}, t) of the reservoir based on the permeability K({right arrow over (r)}, t) of the reservoir and formula (5-14).
In another embodiment, the characteristic parameter may be a skin factor of the reservoir.
The determining a characteristic parameter characterizing the damage extent of the reservoir in the preset region of the well to be diagnosed may include: determining the permeability K({right arrow over (e)}, t) of the reservoir based on the spatio-temporal evolution simulation equation; and determining the skin factor S({right arrow over (r)}, t) of the reservoir based on the permeability K({right arrow over (r)}, t) of the reservoir and formula (5-15).
The characteristic parameters (e.g., the permeability damage rate I(Sw) of the reservoir) obtained by the above embodiments are a result of 4D quantitative simulation of spatio-temporal evolution (as shown in
The metabolism of bacteria is controlled by catalysis of various enzymes in their cells, and the activity of the enzymes is extremely sensitive to temperature; for example, a too high or too low temperature can inactivate the enzymes in the cells. During water injection, a temperature change interval is formed between the interior of a well and a reservoir, and seriously affects the rate of bacterial metabolism, such that bacteria are attached to rock surfaces and form biofilms, thereby reducing the porosity of the reservoir. Thus, the core of the embodiments of the present invention is to establish a variation relationship between an apparent concentration distribution equation of nutrients in a fluid within the reservoir and a temperature and a variation relationship between an apparent concentration distribution equation of bacteria in the fluid and a temperature. Specifically, based on energy conservation, mass conservation, a diffusion relationship, and the like, a spatio-temporal evolution control phenomenological model of concentration distribution of the nutrients and the bacteria in the fluid within the reservoir around a well to be diagnosed is established, and in conjunction with a relationship between concentration distribution and a characteristic parameter characterizing the damage extent of the reservoir such as permeability, spatio-temporal field distribution of the characteristic parameter such as permeability can be diagnosed.
Before performing step S12101, the modeling method further includes: determining a Darcy apparent velocity of a fluid in a reservoir in a preset region of a well to be diagnosed.
Determining the velocity of the fluid in the reservoir may include: establishing a pressure conduction equation for the fluid entering the reservoir; and determining the Darcy apparent velocity of the fluid according to the pressure conduction equation and a Darcy formula.
For the specific determination process, reference can be made to the process of determining the Darcy apparent velocity in the above Embodiment 2 (i.e., the above formulas (2-1) and (2-2) and related description thereof).
Before performing step S12101, the modeling method further includes: determining a temperature distribution field of the reservoir according to the Darcy apparent velocity of the fluid, a thermal conductivity coefficient and thermal diffusivity of the fluid, and an energy conservation law.
During water injection, energy is transferred between reservoir rock and the fluid by means of temperature changes due to a temperature difference between the temperature of the injected water and the temperature of the reservoir. In this case, a mathematical model of field-scale reservoir temperature distribution can be established according to the Darcy apparent velocity u({right arrow over (r)}, t) of the fluid in the reservoir in the preset region of the well to be diagnosed, the thermal conductivity coefficient Dcon and the thermal diffusivity Ddis of the fluid, and the energy conservation law, to obtain an expression of a reservoir temperature distribution control equation.
where T({right arrow over (r)}, t) is the reservoir temperature distribution; and the thermal diffusivity may be expressed by thermal conductivity.
Step S12101: determining a growth rate of the bacteria according to the temperature distribution field of the reservoir and an actual concentration of nutrients in the fluid.
For the step S12101, determining the growth rate of the bacteria (including bacteria attached to rock surfaces and bacteria in the fluid) may include the following steps S12201-S12202, as shown in
Step S12201: determining a maximum growth rate of the bacteria according to the temperature distribution field of the reservoir and a maximum growth rate formula of the bacteria.
Currently, a square root model containing parameters such as activation energy and frequency factor is usually used to describe the maximum growth rate of the bacteria, but the field applicability of this model is poor (mainly embodied in inaccurate prediction results). However, after continuous research, the inventor found that the temperature of the reservoir is a main factor that influences the maximum growth rate of the bacteria and the distribution of the bacteria in the reservoir. Therefore, under the condition of a sufficient injected water source and bacterial nutrient source, the embodiment uses the maximum growth rate formula of the bacteria (including bacteria attached to rock surfaces and bacteria in the fluid) with the temperature as a main variable to simulate the spatial distribution of the maximum growth rate of the bacteria in the reservoir.
For the step S12201, the determining a maximum growth rate of the bacteria includes:
determining the maximum growth rate gmax({right arrow over (r)}, t) of the bacteria according to the temperature distribution field T({right arrow over (r)}, t) of the reservoir and the maximum growth rate formula of the bacteria expressed by the following formula (12-2):
gmax({right arrow over (r)}, t)=[b1(T({right arrow over (r)}, t)−Tmin)]2*{1−exp[c1(T({right arrow over (r)}, t)−Tmax)]}, (12-2)
where b1 and c1 are a first bacterial growth empirical parameter and a second bacterial growth empirical parameter, respectively; Tmax and Tmin are a maximum temperature and a minimum temperature for bacterial growth, respectively.
Step S12202: determining the growth rate of the bacteria according to the maximum growth rate of the bacteria and the actual concentration of the nutrients in the fluid.
Considering that bacterial growth is an irreversible first-order reaction or that bacterial metabolism consumes the nutrients, for the step S12202, the determining the growth rate of the bacteria includes: determining the growth rate gactual bacteria according to the maximum growth rate gmax({right arrow over (r)}, t) of the bacteria, the actual concentration cnu of the nutrients in the fluid, a Monod half growth coefficient ks, and the following formula (12-3):
Step S12102: determining the total amount of the bacteria on rock surfaces according to an amount of the bacteria attached to the rock surfaces in the reservoir associated with both an apparent concentration of the bacteria in the fluid and a total amount of the bacteria on the rock surfaces, and the growth rate and a decay rate of the bacteria.
Specifically, main factors that influence the amount of the bacteria attached to the rock surfaces in the reservoir are bacterial adsorption and desorption rates, and the bacteria attached to the rock surfaces of the reservoir form biofilms, so the amount Cdeposition of the bacteria attached to the rock surfaces in the reservoir can be determined according to a clogging rate kclogging (a constant which may be in 1/day), a declogging rate kdeclogging (a constant which may be in 1/day), the apparent concentration Cbacteriantran of the bacteria in the fluid, the total amount Vbacteriatran of the bacteria on the rock surfaces, and the following formula (12-4):
Cdeposition=kcloggingCbacterian−kdecloggingVbacteriatran. (12-4)
The bacteria attached to the rock surfaces form biofilms, thereby reducing the porosity of the reservoir. In the embodiment, the amount of the bacteria on the rock surfaces is determined mainly by two factors: a net value added by bacterial growth and decay, and the attached amount.
For the step S12102: the determining the total amount of the bacteria on the rock surfaces may include: determining the total amount Vbacteriatran of the bacteria on the rock surfaces according to the amount Cdeposition of the bacteria attached to the rock surfaces in the reservoir, and the growth rate gactual and the decay rate kdecay of the bacteria, and the following formula:
In the above formula (12-5), the variation of the total amount of the bacteria on the rock with time is illustrated on the left side, and a total growth value of metabolic decay of the bacteria on the rock surfaces and a quantity of deposited bacteria are illustrated on the right side respectively.
Step S12103: establishing an apparent concentration distribution equation of the bacteria in the fluid according to the Darcy apparent velocity of the fluid, a dispersion coefficient of the bacteria, the growth rate and the decay rate of the bacteria, and the amount of the bacteria attached to the rock surfaces in the reservoir.
In the embodiment, main considerations are diffusion of the bacteria due to a concentration gradient, convection migration of the bacteria due to water injection, and growth and decay of the bacteria, and the first four factors cause changes in bacterial concentration (including a decrease in concentration due to the formation of biofilms by the bacteria attached to rock of the reservoir). Since the macroscopic effect of irregular motion of bacterial flagellar oscillation in the reservoir is not obvious, the irregular motion effect generated by flagellar oscillation is merged into a convection term in the model; and the effect of bacterial Brownian motion is covered by the effect of bacterial diffusion. In addition, since bacterial chemotaxis has a small effect on distribution rules of the bacteria in the reservoir, this effect is covered by the convective effect, and is combined with the convection term into one term.
For the step S12103, the establishing the apparent concentration distribution equation of the bacteria in the fluid may include: establishing, according to the Darcy apparent velocity u of the fluid, the dispersion coefficient Dsum of the bacteria, the growth rate gactual and the decay rate kdecay of the bacteria, and the amount Cdeposition of the bacteria attached to the rock surfaces in the reservoir, the apparent concentration distribution equation of the bacteria in the fluid expressed by the following formula (12-6):
where cnufrol is apparent concentration distribution of the bacteria in the fluid.
Step S12104: establishing an apparent concentration distribution equation of the nutrients according to the Darcy velocity of the fluid, the dispersion coefficient of the nutrients, the total amount of the bacteria on the rock surfaces, and the apparent concentration of the bacteria.
In the embodiment, main considerations are diffusion of the nutrients due to a concentration gradient, and convection migration of the nutrients due to water injection, and the two factors cause changes in nutrient concentration.
For the step S12104, the establishing the apparent concentration distribution equation of the nutrients may include: establishing, according to the Darcy velocity u({right arrow over (r)}, t) of the fluid, the dispersion coefficient Dnusum of the nutrients, the total amount Vbacteriatran of the bacteria on the rock surfaces in the reservoir, and the apparent concentration Cbacteriatran of the bacteria, the apparent concentration distribution equation of the nutrients expressed by the following formula (12-7):
where gactual is the growth rate of the bacteria; Y is a yield coefficient of the bacteria; and Cnutran is apparent concentration distribution of the nutrients. A one-dimensional form of the above formula (12-7) may be written as
Step S12105: determining a spatio-temporal evolution simulation equation of reservoir damage by bacteria according to the apparent concentration distribution equation of the nutrients and the apparent concentration distribution equation of the bacteria in the fluid.
Wherein the spatio-temporal evolution simulation equation is used to simulate a four-dimensional spatio-temporal evolution process of characteristics of reservoir damage caused by the bacteria.
Specifically, the spatio-temporal evolution simulation equation of reservoir damage by bacteria can be obtained according to the apparent concentration distribution equation of the nutrients expressed by the above formula (12-7) and the apparent concentration distribution equation of the bacteria in the fluid expressed by the above formula (12-6), and in conjunction with other equations (12-1)-(12-5). That is, the spatio-temporal evolution simulation equation of reservoir damage by bacteria is equivalent to an equation set composed of formulas (12-1)-(12-7).
In summary, according to the present invention, the growth rate of the bacteria is creatively determined according to the temperature distribution field of the reservoir and the actual concentration of the nutrients in the fluid; the total amount of the bacteria on rock surfaces is determined according to the amount of the bacteria attached to the rock surfaces in the reservoir, and the growth rate and the decay rate of the bacteria; the apparent concentration distribution equation of the bacteria in the fluid is established according to the Darcy apparent velocity of the fluid, the dispersion coefficient of the bacteria, the growth rate and the decay rate of the bacteria, and the amount of the bacteria attached to the rock surfaces in the reservoir; the apparent concentration distribution equation of the nutrients is established according to the Darcy velocity of the fluid, the dispersion coefficient of the nutrients, the total amount of the bacteria on the rock surfaces, and the apparent concentration of the bacteria; and the spatio-temporal evolution simulation equation of reservoir damage by bacteria is determined according to the apparent concentration distribution equation of the nutrients and the apparent concentration distribution equation of the bacteria in the fluid. Thus, by using the determined spatio-temporal evolution simulation equation, a four-dimensional spatio-temporal evolution process of characteristics of reservoir damage caused by bacteria can be quantitatively simulated. Therefore, performing quantitative prediction of reservoir damage and spatio-temporal deduction of damage laws is of scientific guidance significance for preventing or avoiding reservoir damage, and formulating development plans for oil pools and subsequent well stimulation measures for a well without reservoir damage, and is of very great significance for optimal design of a declogging measure and improvement or restoration of oil well production and water well injection capacity for damaged wells, and improvement of numerical simulation precision of oil pools.
Step S12301: determining an amount of bacteria attached to rock surfaces, based on the spatio-temporal evolution simulation equation established by the modeling method for reservoir damage by bacteria.
For the spatio-temporal evolution simulation equation of reservoir damage by bacteria expressed by the above formulas (12-6)-(12-7), the amount of the attached bacteriadeposition can be solved by referring to the process of solving the volume concentration of deposited particles in the above Embodiment 1.
The amount of the attached bacteria Cdeposition is calculated by the above method. As the spatio-temporal evolution simulation equation established by the above modeling method for reservoir damage by bacteria comprehensively considers the influence of various physical and chemical factors on reservoir damage during migration of the bacteria and nutrients in the fluid, the amount of the attached (deposited) bacteria obtained by the step S12301 is very precise.
Step S12302: determining a characteristic parameter characterizing the damage extent of the reservoir in the preset region of the well to be diagnosed, based on the amount of the bacteria attached to the rock surfaces.
Wherein the characteristic parameter is permeability of the reservoir.
For the step S12302, the determining a characteristic parameter characterizing the damage extent of the reservoir in the preset region of the well to be diagnosed may include: determining the permeability K({right arrow over (r)}, t) of the reservoir based on the amount Cdeposition of the bacteria attached to the rock surfaces and the density ρ of the bacteria and the following formula (12-8),
where ϕ0 is an initial value of porosity of the reservoir; and Ko ({right arrow over (r)}) is an initial value of the permeability of the reservoir. Wherein the characteristic parameter is a skin factor of the reservoir.
For the step S12302, the determining a characteristic parameter characterizing the damage extent of the reservoir in the preset region of the well to be diagnosed includes: determining the permeability K({right arrow over (r)}, t) of the reservoir based on the apparent concentration Cnutran({right arrow over (r)}, t) and the actual concentration cnu({right arrow over (r)}, t) of the nutrients in the fluid and formula
and determining the skin factor S({right arrow over (r)}, t) of the reservoir based on the permeability K({right arrow over (r)}, t) of the reservoir and the following formula (12-9):
where Ko ({right arrow over (r)}) is an initial value of the permeability of the reservoir; and
The characteristic parameter (e.g., the permeability K({right arrow over (r)}, t) and the skin factor S({right arrow over (r)}, t) of the reservoir) obtained by the step S12302 is a result of 4D quantitative simulation of spatio-temporal evolution (not shown). Therefore, quantitative prediction of reservoir damage and spatio-temporal deduction of damage laws can be performed according to evolution characteristics of the permeability or the skin factor, which is of scientific guidance significance for preventing or avoiding reservoir damage, and formulating development plans for oil pools and subsequent well stimulation measures.
In summary, the apparent concentration of the nutrients in the fluid can be determined by using the determined spatio-temporal evolution simulation equation, and then the characteristic parameter (e.g., the permeability and/or the skin factor of the reservoir) characterizing the damage extent of the reservoir in the preset region of the well to be diagnosed can be determined based on the apparent concentration and the actual concentration of the nutrients in the fluid, whereby a four-dimensional spatio-temporal evolution process of the reservoir damage characteristic caused by bacteria can be simulated quantitatively, Therefore, performing quantitative prediction of reservoir damage and spatio-temporal deduction of damage laws is of scientific guidance significance for preventing or avoiding reservoir damage, and formulating development plans for oil pools and subsequent well stimulation measures for a well without reservoir damage, and is of very great significance for optimal design of a declogging measure and improvement or restoration of oil well production and water well injection capacity for damaged wells, and improvement of numerical simulation precision of oil pools.
According to the adsorption mechanism of polymers on the surfaces of porous media, adsorption of polymers on a reservoir includes layer adsorption (as shown in
The essence of clogging by polymer adsorption is migration and adsorption of polymers in a fluid within the reservoir. Thus, the core of the embodiments of the present invention is to establish a kinetic model of migration and adsorption of a polymer according to the law of mass conservation. Specifically, based on mass conservation, a diffusion relationship, and the like, a spatio-temporal evolution control phenomenological model (containing a concentration C of the polymer in a fluid) of concentration distribution of the polymer in a reservoir around a well to be diagnosed is established, and in conjunction with a relationship between the concentration C and a characteristic parameter characterizing the damage extent of the reservoir such as permeability, spatio-temporal field distribution of the characteristic parameter such as permeability can be diagnosed.
Step S13201: determining a velocity of a fluid in a reservoir.
Wherein the reservoir is located in a preset region of a well to be diagnosed (e.g., a polymer injection well).
For the step S13201, determining the velocity of the fluid containing the polymer in the reservoir includes: establishing a pressure conduction equation for the fluid entering the reservoir; and determining the velocity of the fluid according to the pressure conduction equation and a Darcy formula.
For the specific determination process, reference can be made to the process of determining the velocity of the fluid in the above Embodiment 1 (i.e., the above formulas (1-1) and (1-2) and related description thereof).
At the time t, a total number of of molecular chains of the polymer contained in a unit control volume of the fluid at a location {right arrow over (r)} in a porous medium of the reservoir is N ({right arrow over (r)}, t), and the unit control volume flows at a velocity ν in a pore for time Δt with a displacement r. Layer adsorption and bridging adsorption cause the number of molecular chains of the polymer adsorbed on the reservoir to change, so the number of molecular chains of the polymer contained in the unit control volume becomes N ({right arrow over (r)}+Δ{right arrow over (r)}, t+Δt). During the above adsorption process, the total number of polymer molecular chains in the entire reservoir remains unchanged.
Before performing step S13202, the modeling method may further include: establishing, based on the number of molecular chains of the polymer subjected to layer adsorption on the reservoir and the number of molecular chains of the polymer subject to bridging adsorption on the reservoir, a number distribution equation of molecular chains of the polymer in the fluid expressed by the following formula (13-1):
where χ is a mass proportion of the polymer subjected to layer adsorption at any location {right arrow over (r)} of the reservoir at time t (i.e., the proportion of the mass of the polymer subject to layer adsorption at any location {right arrow over (r)} in the reservoir in a total mass of all adsorbed polymer at that location at time t);
where m is an average mass of molecular chains of the adsorbed polymer; ρs is density of the adsorbed polymer; ASl and ASb are a specific surface area of the reservoir where layer adsorption occurs and a specific surface area of the reservoir where bridging adsorption occurs, respectively; ϕ0 is porosity of the reservoir (the porosity of the reservoir can be considered to be unchanged in various embodiments of the present invention); Γl, Γb are layer adsorption density and bridging adsorption density, respectively; and N ({right arrow over (r)}, t) is the number of molecular chains of the polymer in the fluid at any location {right arrow over (r)} in the reservoir at time t,
Step S13202: establishing a balance equation of a mass concentration of the polymer in the fluid and a proportion distribution equation of molecular chains of adsorbed polymer adsorbed on the reservoir in the polymer according to the velocity of the fluid, the number distribution equation of the molecular chains of the polymer in the fluid and a diffusion coefficient of the polymer.
Wherein the proportion of the molecular chains of the adsorbed polymer is the proportion of a number of the molecular chains of the adsorbed polymer in an initial number of molecular chains of the polymer in the fluid.
For the step S13202, the establishing a balance equation of a mass concentration of the polymer in the fluid and a proportion distribution equation of molecular chains of adsorbed polymer adsorbed on the reservoir in the polymer may include the following steps S13301-S13303, as shown in
Step S13301: determining a function relationship in which a mass concentration of the adsorbed polymer varies with layer adsorption density and bridging adsorption density, based on the number distribution equation of the molecular chains of the polymer in the fluid.
Both the left and right sides of the above formula (13-1) represent the number of molecular chains subjected to polymer adsorption per unit control volume. Multiplying the two sides of formula (13-1) by m (the mass of the polymer subject to adsorption) and taking a limit value Δt→0 yields a formula with two sides that can represent a molecular mass (i.e., mass concentration Cd(Γl, Γb)) of the polymer subject to polymer adsorption per unit control volume,
represents a mass concentration of the polymer subjected to layer adsorption; and
represents a mass concentration of the polymer subject to bridging adsorption.
Step S13302: establishing the balance equation of the mass concentration of the polymer in the fluid according to the velocity of the fluid, the diffusion coefficient of the polymer, and the function relationship in which the mass concentration of the adsorbed polymer varies with the layer adsorption density and the bridging adsorption density.
For the adsorption process of the polymer, the polymer molecular chains follow the law of mass conservation, i.e., the change in the mass of flowing polymer molecules is equal to the negative change in the mass of polymer molecules adsorbed on surfaces of the porous medium of the reservoir:
where j ({right arrow over (r)}, t) is a diffusion flow rate of the polymer, and
where D represents the diffusion coefficient of the polymer; md ({right arrow over (r)}, t) represents the mass of the polymer subjected to adsorption; u ({right arrow over (r)}, t) is the velocity of the fluid; ρ is density of the fluid; w ({right arrow over (r)}, t) is a mass fraction of the flowing polymer molecules in the fluid; and ϕ0 is an initial value of the porosity of the reservoir. If the change in the mass per unit control volume (i.e., the change in the mass concentration) is considered, due to microscopic changes within the pores (not considering the porosity ϕ0), the above formula can be changed to the mass concentration shown in following formula (13-3).
For the step S13302, the establishing the balance equation of the mass concentration of the polymer in the fluid may include: establishing, according to the velocity ν({right arrow over (r)}, t) of the fluid, the diffusion coefficient D of the polymer, and the function relationship Cd (Γl, Γb) in which the mass concentration of the adsorbed polymer varies with the layer adsorption density and the bridging adsorption density, the balance equation of the mass concentration of the polymer in the fluid expressed by the following formula (13-3):
where C({right arrow over (r)}, t) is the mass concentration of the polymer in the fluid (C=ρw({right arrow over (r)}, t)); and Γl and Γb are the layer adsorption density and the bridging adsorption density, respectively.
Step S13303: establishing the proportion distribution equation of the molecular chains of the adsorbed polymer according to the number distribution equation of the molecular chains of the polymer in the fluid.
For the step S13303, the establishing the proportion distribution equation of the molecular chains of the adsorbed polymer may include: establishing, according to the number distribution equation of the molecular chains of the polymer in the fluid and
the proportion distribution equation of the molecular chains of the adsorbed polymer expressed by the following formula (13-4):
where C0 is an initial mass concentration of the polymer in the fluid; N0 is an initial number of the molecular chains of the polymer in the fluid; Ppoly is the proportion of the number of the molecular chains of the adsorbed polymer in the initial number of the molecular chains of the polymer in the fluid; and ν({right arrow over (r)}, t) is the velocity of the fluid.
Specifically, Taylor first-order expansion (ignoring the second-order term) is carried out for the first term at the left side of the number distribution equation of the molecular chains of the polymer in the fluid expressed by the equation (13-1), and then both sides of the equation are simultaneously divided by Δt, where
and when Δt→0, the above formula (13-1) becomes:
represent an adsorption rate for layer adsorption and an adsorption rate for bridging adsorption (i.e., the number of polymer molecular chains adsorbed per unit time), respectively:
where ρs represents the density of the adsorbed polymer; Asl and Asb represent a specific surface area of the medium where layer adsorption occurs and a specific surface area of the medium where bridging adsorption occurs, respectively (there may be Asl=Asb=As); ϕ0 represents the porosity of the reservoir; Γl represents the layer adsorption density of the polymer; and Γb represents the bridging adsorption density of the polymer. Substituting the above two formulas into formula (7) yields the following formula:
and then dividing both sides of the above formula (13-6) by N0 (i.e., the initial number of the molecular chains of the polymer in the fluid) yields the proportion distribution equation of the molecular chains of the adsorbed polymer expressed by formula (13-4).
Step S13203: determining a layer adsorption density and a bridging adsorption density of the polymer when an adsorption of the polymer does not reach a saturation state, according to the layer adsorption rate and the layer desorption rate, the bridging adsorption rate and the bridging desorption rate, and the law of conservation of kinetic energy.
The adsorption process of the polymer also follows the law of conservation of kinetic energy, i.e., adsorption and desorption of polymer molecules per unit time meets a dynamic equilibrium and finally reach an adsorption saturation state. The desorption and adsorption dynamic equilibrium of layer adsorption and bridging adsorption can be expressed as formulas (13-7) and (13-8).
For part of the step S13203, determining the layer adsorption density of the polymer when the adsorption of the polymer does not reach the saturation state may include: determining the layer adsorption density of the polymer according to the layer adsorption rate and the layer desorption rate and the law of conservation of kinetic energy expressed by the following formula:
where κal is the layer adsorption rate; κdl is the layer desorption rate; Γl∞ is a saturated layer adsorption density of the polymer when the adsorption of the polymer reaches the saturation state; and Γl is the layer adsorption density.
For part of the step S13203, determining the bridging adsorption density of the polymer when the adsorption of the polymer does not reach the saturation state may include: determining the bridging adsorption density of the polymer when the adsorption of the polymer does not reach saturation according to the bridging adsorption rate and the bridging desorption rate and the law of conservation of kinetic energy expressed by the following formula:
where κab is the bridging adsorption rate; κdb is the bridging desorption rate; Γb∞ is a saturated bridging adsorption density of the polymer when the adsorption of the polymer reaches the saturation state; Γb is the bridging adsorption density; and Ppoly is the proportion of the number of the molecular chains of the adsorbed polymer in the initial number of the molecular chains of the polymer in the fluid.
In the above two formulas (13-7)-(13-8), on the left side, the term
represents the layer adsorption density of the polymer per unit time when the adsorption of the polymer does not reach the saturation state, and on the right side, κai(Γi∞−Γi) represents dynamic adsorption density of the polymer per unit time, and κdiΓi represents dynamic desorption density per unit time, where i=l (layer adsorption), b (bridging adsorption).
Step S13204: determining a spatio-temporal evolution simulation equation of reservoir damage by the polymer according to the balance equation of the mass concentration of the polymer in the fluid, the proportion distribution equation of the molecular chains of the adsorbed polymer, the layer adsorption density and the bridging adsorption density.
Substituting formula (13-2) to the right side of formula (13-3) yields the following formula (13-9):
and then substituting the above formulas (13-7)-(13-8) into formulas (13-4) and (13-9) yields the spatio-temporal evolution simulation equation of reservoir damage by the polymer. That is, the spatio-temporal evolution simulation equation of reservoir damage by the polymer is equivalent to an equation set composed of formulas (13-4), and (13-7)-(13-9).
In summary, according to the present invention, the balance equation of the mass concentration of the polymer in the fluid and the proportion distribution equation of the molecular chains of adsorbed polymer adsorbed on the reservoir in the polymer is creatively established according to the velocity of the fluid, the number distribution equation of the molecular chains of the polymer in the fluid and the diffusion coefficient of the polymer; the layer adsorption density and the bridging adsorption density of the polymer when the adsorption of the polymer does not reach the saturation state is determined according to the layer adsorption rate and the layer desorption rate, the bridging adsorption rate and the bridging desorption rate, and the law of conservation of kinetic energy; and the spatio-temporal evolution simulation equation of reservoir damage by the polymer is determined according to the balance equation of the mass concentration of the polymer in the fluid, the proportion distribution equation of the molecular chains of the adsorbed polymer, the layer adsorption density and the bridging adsorption density. Thus, by using the determined spatio-temporal evolution simulation equation, a four-dimensional spatio-temporal evolution process of reservoir damage characteristics caused by the polymer can be quantitatively simulated. Therefore, performing quantitative prediction of reservoir damage and spatio-temporal deduction of damage laws is of scientific guidance significance for preventing or avoiding reservoir damage, and formulating development plans for oil pools and subsequent well stimulation measures for a well without reservoir damage, and is of very great significance for optimal design of a declogging measure and improvement or restoration of oil well production and water well injection capacity for damaged wells, and improvement of numerical simulation precision of oil pools.
Step S13401: determining a mass concentration of the polymer in the fluid based on the spatio-temporal evolution simulation equation established by the modeling method for reservoir damage by the polymer.
For the polymer migration equation of reservoir damage by polymer adsorption expressed by the above formula (13-7), the volume concentration C({right arrow over (r)}, t) of the polymer can be calculated by referring to the process of solving the volume concentration of the deposited particles in the above Embodiment 1.
The mass concentration C({right arrow over (r)}, t) of the polymer in the fluid can be calculated by the above method. As the spatio-temporal evolution simulation equation established by the above modeling method of reservoir damage by polymer adsorption comprehensively considers the influence of various physical and chemical factors on reservoir damage during polymer adsorption within the reservoir, the mass concentration of the polymer in the fluid obtained by the step S13401 is very precise.
Step S13402: determining a characteristic parameter characterizing the damage extent of the reservoir in the preset region of the well to be diagnosed, based on the mass concentration of the polymer in the fluid.
Wherein the characteristic parameter is permeability of the reservoir.
For the step S13402, the determining a characteristic parameter characterizing the damage extent of the reservoir in the preset region of the well to be diagnosed may include: determining the permeability K({right arrow over (r)}, t) of the reservoir based on the mass concentration C({right arrow over (r)}, t) of the polymer in the fluid and formula (13-10):
where ϕ0 is an initial value of porosity of the reservoir; C0({right arrow over (r)}, t) is an initial mass concentration of the polymer in the fluid; C∞ is a maximum adsorption mass concentration of the polymer; mK is a second empirical value; and Ko({right arrow over (r)}) is an initial value of the permeability of the reservoir.
Wherein the characteristic parameter is a fluid loss coefficient of the reservoir.
For the step S13402, the determining a characteristic parameter characterizing the damage extent of the reservoir in the preset region of the well to be diagnosed may include: determining the fluid loss coefficient k({right arrow over (r)}, t) of the reservoir based on the mass concentration of the polymer in the fluid and formula (13-11):
where ϕ0 is the initial value of the porosity of the reservoir; C0({right arrow over (r)}, t) is the initial mass concentration of the polymer in the fluid; C∞ is the maximum adsorption mass concentration of the polymer; mk is a first empirical value; and k0({right arrow over (r)}) is an initial value of the fluid loss coefficient of the reservoir.
Wherein the characteristic parameter is a skin factor of the reservoir.
For the step S13402, the determining a characteristic parameter characterizing the damage extent of the reservoir in the preset region of the well to be diagnosed includes: determining the permeability K({right arrow over (r)}, t) of the reservoir based on the mass concentration C({right arrow over (r)}, t) of the polymer in the fluid and formula
and determining the skin factor S({right arrow over (r)}, t) of the reservoir based on the permeability K({right arrow over (r)}, t) of the reservoir and formula (13-12):
where Ko ({right arrow over (r)}) is the initial value of the permeability of the reservoir; ϕ0 is the initial value of the porosity; C0({right arrow over (r)}, t) is the initial mass concentration of the polymer in the fluid; mK is the second empirical value;
The characteristic parameter (e.g., the permeability K({right arrow over (r)}, t) and the skin factor S({right arrow over (r)}, t) of the reservoir) obtained by the step S13402 is a result of 4D quantitative simulation of spatio-temporal evolution (as shown in
where Kmax({right arrow over (r)}, t) is a maximum value of K({right arrow over (r)}, t)), and a working person concerned can visually confirm the damage extent of the reservoir from
In the case where the permeability of the reservoir when each factor damages the reservoir is obtained in each of the above embodiments, the skin factor (or permeability damage rate) of the reservoir when each factor damages the reservoir can be obtained according to a relational expression (e.g., formula (3-8) or (5-14)) between the permeability and the skin factor (or permeability damage rate). Then, for various stages of the same diagnostic well, weighted summation is performed on the permeability (or skin factor or permeability damage rate) of the reservoir when each of the plurality of factors damages the reservoir respectively to determine the permeability (or skin factor or permeability damage rate) of the reservoir when the plurality of factors damage the reservoir simultaneously. Further, a contribution proportion of the permeability of the reservoir when each factor damages the reservoir to the total permeability (or total skin factor or total permeability damage rate) can also be determined.
In summary, according to the present invention, the mass concentration of the polymer in the fluid can be calculated creatively by using the determined spatio-temporal evolution simulation equation, then the characteristic parameter (e.g., the permeability and/or the skin factor of the reservoir) characterizing the damage extent of the reservoir in the preset region of the well to be diagnosed can be determined based on the determined mass concentration of the polymer, and thus, a four-dimensional spatio-temporal evolution process of characteristics of reservoir damage caused by polymer adsorption can be quantitatively simulated. Therefore, performing quantitative prediction of reservoir damage and spatio-temporal deduction of damage laws is of scientific guidance significance for preventing or avoiding reservoir damage, and formulating development plans for oil pools and subsequent well stimulation measures for a well without reservoir damage, and is of very great significance for optimal design of a declogging measure and improvement or restoration of oil well production and water well injection capacity for damaged wells, and improvement of numerical simulation precision of oil pools.
Field Verification and Application
Using relevant parameters provided on site, a total skin factor (as shown in
The following is an example of well G01s1 in Suizhong oilfield 36-1. Relevant parameters and measured skin factors of the well are shown in Table 3.
Table 3 Measured results of basic parameters and skin factors of well G01s1 in Suizhong oilfield 36-1
The parameters in Table 3 are input into software loaded with the above methods to obtain variations of damaged zone radii and damage extents (e.g., permeability damage rates) caused by six damage types: “sand production, stress sensitivity, emulsification, fine particle migration, wettability reversal, and organic scale” with time (as shown in
The simulation results show that at day 40, a total skin factor is 35.65, a measured skin factor is 36.5, and a pseudo skin factor due to well deviation and partial completion is 0.49, with a relative error of only 0.99%, which indicates very high conformity. Analysis based on the damage factors shows that the major damage factors are emulsification, fine particle migration, organic scale, and sand production for reservoir damage, while the damage caused by wettability reversal and stress sensitivity is very small. Of course, it is also possible to simulate the distribution and evolution of any reservoir damage parameter in a 4D space, such as at different depth points and in different directions, as needed.
In addition, the 4D quantitative simulation technology of spatio-temporal evolution established herein has been widely verified and applied in different types of oilfields in China, such as those developed by China National Offshore Oil Corporation, with an average accuracy of 95% or above, which fully proves the reliability, accuracy and practicality of this simulation technology. A comparison of simulation results for some of the wells is shown in Table 4.
Table 4 Comparison of simulation calculated values and measured values of skin factors of some wells
Thus, the embodiments described above can achieve the following advantages.
In summary, according to the present invention, based on a spatio-temporal evolution simulation equation of reservoir damage by each of a plurality of factors, a characteristic parameter characterizing reservoir damage by each of the plurality of factors is creatively determined; and the effective characteristic parameter characterizing the damage extent of the reservoir is determined based on the characteristic parameter characterizing reservoir damage by each of the plurality of factors. Thus, by using the spatio-temporal evolution simulation equations of reservoir damage by the relevant factors, the characteristic parameters (such as permeability) of reservoir damage caused by the factors respectively and a total characteristic parameter (such as total permeability or effective permeability) of reservoir damage caused by the plurality of relevant factors can be quantitatively simulated. Therefore, performing quantitative prediction of reservoir damage and spatio-temporal deduction of damage laws is of scientific guidance significance for preventing or avoiding reservoir damage, and formulating development plans for oil pools and subsequent well stimulation measures for a well without reservoir damage, and is of very great significance for optimal design of a declogging measure and improvement or restoration of oil well production and water well injection capacity for damaged wells, and improvement of numerical simulation precision of oil pools.
Correspondingly, an embodiment of the present invention further provides a system for determining a damage extent of a reservoir, the system including: a first parameter determination device (not shown) configured to, based on a spatio-temporal evolution simulation equation of reservoir damage by each of a plurality of factors, determine a characteristic parameter characterizing reservoir damage by each of the plurality of factors, wherein the reservoir is located in a preset region of a well to be diagnosed; and a second parameter determination device (not shown) configured to determine an effective characteristic parameter characterizing the damage extent of the reservoir based on the characteristic parameter characterizing reservoir damage by each of the plurality of factors.
The system for determining the damage extent of the reservoir has the same advantages as the above-mentioned method for determining the damage extent of the reservoir with respect to the prior art, which will not be described here.
Correspondingly, an embodiment of the present invention further provides a machine-readable storage medium that stores instructions which are configured to cause a machine to execute the method for determining the damage extent of the reservoir.
The machine-readable storage medium includes, but is not limited to, phase-change memory (short for phase change random access memory, PRAM, also called RCM/PCRAM), static random access memory (SRAM), dynamic random access memory (DRAM), other types of random access memory (RAM), read-only memory (ROM), electrically erasable programmable read-only memory (EEPROM), flash memory or other memory technologies, compact disc read-only memory (CD-ROM), digital versatile disk (DVD) or other optical storage, magnetic cartridge tape, magnetic tape disk storage or other magnetic storage devices, and various other media that can store program codes.
The steps in the above embodiments can be performed by a computer, and the processing procedures of various physicochemical quantities involved in certain steps (e.g., steps S1101-S1104) achieve simulation of a spatio-temporal evolution field of reservoir damage by the factors, and the processing procedures of various physicochemical quantities involved in certain steps (e.g., steps S1201-S1202) achieve prediction of spatio-temporal evolution of reservoir damage by the factors. The steps in the above embodiments can be performed by a processor.
Preferred implementations of the present invention are described above in detail with reference to the accompanying drawings. However, the present invention is not limited to the specific details in the above implementations. Within the scope of the technical concept of the present invention, various simple modifications can be made to the technical solutions of the present invention, and these simple modifications are all encompassed within the protection scope of the present invention.
In addition, it should be noted that the specific technical features described in the above specific implementations may be combined in any suitable manner without contradiction. To avoid unnecessary repetition, various possible combinations will not be described separately in the present invention.
In addition, various different implementations of the present invention may also be combined optionally, and the combinations should also be regarded as contents disclosed in the present invention so long as they do not depart from the idea of the present invention.
| Number | Date | Country | Kind |
|---|---|---|---|
| 202010873143.9 | Aug 2020 | CN | national |
This application is a bypass continuation of PCT/CN2021/114851, filed Aug. 26, 2021, which claims benefits of Chinese patent application 202010873143.9 filed on Aug. 26, 2020, which is hereby incorporated by reference in its entirety.
| Number | Date | Country |
|---|---|---|
| 103806904 | May 2014 | CN |
| 105181927 | Dec 2015 | CN |
| 108240955 | Jul 2018 | CN |
| 108804779 | Nov 2018 | CN |
| 109594982 | Apr 2019 | CN |
| 109946436 | Jun 2019 | CN |
| 110376112 | Oct 2019 | CN |
| Entry |
|---|
| Zou et al., Investigating Reservoir Pressure Transmission for Three Types of Coalbed Methane Reservoirs in the Qinshui Basin in Shan'xi Province, China, 2013, Petroleum Geoscience, vol. 19, pp. 375-383 (Year: 2013). |
| English Translation of International Search Report from PCT/CN2021/114851 dated Oct. 28, 2021 (3 pages). |
| Jiang et al. Quantitative Prediction Method of Water Blocking Damage in Low and Extra-Low Permeability Reservoir) Petroleum Drilling Techniques, vol. 40, No. 1, Jan. 31, 2012, p. 69-73 (English Abstract). |
| Wang et al. Mathematical Models Forecasting the Formation Impairment Induced FR-OM Drilling Fluid Drilling & Production Technology Dec. 31, 2002 , pp. 6 and 80-83 (English Abstract). |
| Yi et al. Numerical Simulation of Reservoir Damage by Oil/Water Emulsification Based on Hydrodynamic Instability[J]. Drilling Fluid & Completion Fluid, 2017, 34(6): pp. 22-128. doi: 10.3969/j.issn.1001-5620.2017.06.02 with English abstract. |
| Li et al., A review of the current status of induced seismicity monitoring for hydraulic fracturing in unconventional tight oil and gas reservoirs, Fuel, vol. 242, 2019, pp. 195-210. |
| Number | Date | Country | |
|---|---|---|---|
| 20230026538 A1 | Jan 2023 | US |
| Number | Date | Country | |
|---|---|---|---|
| Parent | PCT/CN2021/114851 | Aug 2021 | US |
| Child | 17864224 | US |
