Introduction

At present, the main oil fields in China have entered the late stage of development after a long period of exploitation. Most of them have entered the stage of high and extremely high water cut, among which the water cut of the major reservoirs is higher. Under such circumstances, it is urgent to carry out depth development of the old oil field and take corresponding measures to enhance the recovery (Nan et al. 2015). At the same time, with the increase in water cut, the distribution of underground oil and water presents a pattern of "overall high dispersion and local relative enrichment" (Dakuang 2007). Accurate prediction of the distribution state of the remaining oil in the reservoir has become the basis and key for improving oil recovery. From 2002 to 2004, Shengli oilfield used petroleum geology, geophysics, reservoir engineering, computer technology and other knowledge to predict and infuse the well pattern in the remaining oil accumulation area. It brings a cumulative oil increase of 875 × 104t and a profit of 28.9 × 108yuan. The large-scale systematic technology was developed and brought high economic benefits at the same time. The comprehensive application of the knowledge and technology of various disciplines to guide the exploitation of the remaining oil will achieve good practical results (Yang 2003; Longxin 1999).

The flow unit study method is the main method to guide the exploration of residual oil. The flow unit study method can simultaneously reflect both the physical properties of the reservoir and the seepage characteristics of the fluid (Shi-Mi 2001a). It is of great significance to the prediction and adjustment of the remaining oil accumulation area in the oil field. However, the research method of reservoir flow units is not completely consistent due to the limitations of specific geological conditions and actual data (Hearn et al. 1984). Past research method of flow unit, which is mainly composed of geological research, can be divided into the flow unit division method (method of sedimentary facies, the analytic hierarchy process (AHP), the heterogeneous synthetic index method) and is given priority to with the mathematical method of reservoir parameter analysis method (flow layer index FZI method (Ebanks 1987), pore throat geometric shape method, production dynamic parameter method, multi-parameter synthesis) two kinds (Amaefule and Altunbay 1993). The key problems existing in current flow units are as follows: the hierarchical structure analysis of flow units (Liu and Xu 2003). The division method of flow units based on geological research is mainly realized by statistical analysis, and there is no clear theoretical formula and division standard (Mingzhen et al. 2015). Due to the inherent particularity of the flow unit, there is no good distinction between the flow unit and the physical phase of the rock (Miall 1996). The dynamic concept of flow does not well reflect the distribution law of flow and change in reservoir fluid at different development stages (Gao et al. 2005). To sum up, the existing flow unit can only be used to evaluate the quality of the reservoir. It cannot reflect the oil and water distribution rules at different development stages (Jin-Ying et al. 2006). And it is difficult to objectively analyze the scope and influence of the reservoir invalid flooding at the later stage of development. It is also an urgent problem to be solved in the current flow unit research.

In 2016, Wu et al. put forward the theory of EWU and defined the concepts and internal relations of flow unit, drive unit, and EWU (Yan-Xin et al. 2004). And based on hydrodynamics principle, two-phase seepage flow theory and mathematical statistic method, the “flow unit” is divided into the criteria of high speed and invalid drive, to establish the basic theory and division method of EWU (Shi-Mi 2001b). The advantage of EWU mainly reflects in: EWU each flow unit is defined based on the theory of fluid mechanics and seepage mechanics and has clear physical meaning and classification standard; the science of flow units is classified and divided into the flow speed; is described in addition to the effective drive and invalid drive in the process of displacement; makes it match with the development effect, good combining geology, and reservoir; describes the late development of invalid drive problem; and provides a basis for the oilfield development adjustment. Moreover, the theory of EWU is proposed based on the displacement process of crude oil in a porous medium of reservoir, which solves the problem that the existing flow unit cannot reflect the flow regularity of fluid in reservoir. The classification criterion is related to the parameters such as water injection time, which solves the problem that the distribution law of oil and water dynamics cannot be reflected in the development stage (Bandilla et al. 2009).

In this paper, based on the two-dimensional EWU, a three-position EWU model is established using shape function and the superposition principle of hydrodynamic potential. Then, the streamline distribution of reservoir is obtained, and the corresponding water cut and saturation changes are obtained through the saturation distribution model before water breakthrough. The model is applied to the reservoirs under rhythmic conditions to obtain the corresponding streamline and saturation distribution. The structure of this paper is as follows: The second section introduces the content and division principle of two-dimensional EWU theory and the establishment of three-dimensional EWU MODEL; in Sect. 3, the flow line and saturation distribution of the reservoir are obtained according to the positive rhythm, anti-rhythm and compound rhythm, and the corresponding mining suggestions are given according to the obtained results. Section 4 is the conclusion.

Model Introduction

Introduction to the theory of two-dimensional EWU

Definition and classification of flow unit, drive unit, displacement unit and EWU:

Flow unit (Jin-Ying et al. 2006): It is the seepage area with similar flow laws and reservoir characteristics within the range of water flooding in the reservoir. Based on the theory of hydrodynamics, the “flow unit” is divided into “high-speed flow area” and “low-speed flow area" through the "non-uniform flow distribution curve.”

EWU: Based on the definition of the drive unit, according to the two-phase seepage theory, the effective and invalid flooding in the flow unit is divided. And the area with a water production rate higher than the limit water cut of 98% is defined as invalid flooding. Based on the hydrodynamic theory and the flow tube model, through "flow contribution rate, flow non-uniform distribution curve, flow intensity difference coefficient and limit water cut," the pressure gradient control that will overcome the starting pressure gradient and contribute to production ‘Area” is divided into Class I (high-speed flow invalid flooding), Class II (high-speed flow effective flooding), Class III (low-speed flow invalid flooding), Class IV (low-speed flow effective flooding).

Establishment of three-dimensional EWU model

In this study, the flow function method (Fen and Kiuchi 2000) is used to study the fluid flow characteristics of thick oil layers under different well pattern types under the three-dimensional condition of gravity. Since the flow function is defined by the incompressible fluid under two-dimensional plane conditions, a single flow function cannot describe the three-dimensional flow problem. Moreover, for three-dimensional heterogeneous reservoirs, heterogeneous boundary shapes are diverse, and it is difficult to characterize and calculate in rectangular coordinate systems. Therefore, in order to solve the three-dimensional flow problem of heterogeneous reservoir, the orthogonal coordinate system is used to characterize the heterogeneity distribution characteristics. Then, the three-dimensional streamline is represented by the intersection of two surfaces in the coordinate system. Basic assumptions are made:

  1. 1.

    The fluid is incompressible;

  2. 2.

    Rock is slightly compressible;

  3. 3.

    Regardless of changes in reservoir temperature;

  4. 4.

    The fluid component is oil–water two-phase;

  5. 5.

    The oil and water phases are not miscible;

  6. 6.

    Consider the influence of gravity;

For the incompressible deformation of the volume, the divergence in the velocity field is 0, and the interior of the flow field is a passive field; then, the velocity vector at a certain point can be expressed as (Zhenfan and Zhiye 2000):

$$\vec{v} = {\text{grad}} \varpi = \frac{\partial \varpi }{x}i + \frac{\partial \varpi }{y}j + \frac{\partial \varpi }{z}k = v_{x} \vec{i} + v_{t} \vec{j} + v_{z} \vec{k}$$
(1)

where \(\varpi = \varpi \left( {x,y,z} \right)\) is the speed potential function; \(v_{x} ,v_{y} ,v_{z}\) is the component of velocity in three directions of three-dimensional space.

According to the principle that the streamline is orthogonal to the equipotential surface, the direction of the streamline is always the same as the direction of the potential gradient. This principle is applicable in all coordinate systems. Therefore, for a three-dimensional flow in a curvilinear coordinate system, the velocity vector can be expressed by two flow surface functions \(\psi\) and \(\phi\); then, the velocity vector can be expressed as:

$$\vec{v} = \nabla \psi \left( {x,y,z} \right) \times \nabla \phi \left( {x,y,z} \right) = v_{x} \vec{i} + v_{y} \vec{j} + v_{z} \vec{k}$$
(2)

where \(\psi\) and \(\phi\) are two stream functions in the plane and in the vertical direction, \(v_{x} , \, v_{y} , \, v_{z}\) are the velocities in the x, y, z directions.

Among them,

$$\nabla \psi \left( {x,y,z} \right) \times \nabla \phi \left( {x,y,z} \right) = \left| {\begin{array}{*{20}c} i & j & k \\ {\frac{\partial \psi }{x}} & {\frac{\partial \psi }{y}} & {\frac{\partial \psi }{z}} \\ {\frac{\partial \phi }{x}} & {\frac{\partial \phi }{y}} & {\frac{\partial \phi }{z}} \\ \end{array} } \right|$$
(3)

So for rectangular coordinates, the speed can be expressed as:

$$\left\{ {\begin{array}{*{20}l} {v_{x} = \frac{\partial \psi }{{\partial y}}\frac{\partial \phi }{{\partial z}} - \frac{\partial \psi }{{\partial z}}\frac{\partial \phi }{{\partial y}}} \hfill \\ {v_{y} = \frac{\partial \psi }{{\partial z}}\frac{\partial \phi }{{\partial x}} - \frac{\partial \psi }{{\partial x}}\frac{\partial \phi }{{\partial z}}} \hfill \\ {v_{z} = \frac{\partial \psi }{{\partial x}}\frac{\partial \phi }{{\partial y}} - \frac{\partial \psi }{{\partial y}}\frac{\partial \phi }{{\partial x}}} \hfill \\ \end{array} } \right.$$
(4)

Determination of shape function of rhythm reservoir

For the flow rate in the injection-production well pattern, q is the flow rate and A is the actual flow area perpendicular to the main permeability direction. The XZ plane is shown in Fig. 1:

Fig. 1
figure 1

Schematic diagram of flow velocity and section

Considering the relationship between sections A, Ax, v, vx, we can get:

$$A_{x} = A \cdot \cos \theta ,v_{x} = v \cdot \cos \theta$$
(5)

\(\theta\) is the angle between the velocity direction and the horizontal direction, so:

$$v_{x} = \frac{q}{{A_{x} }} \cdot \cos^{2} \theta$$
(6)

Among them,\(q_{x} = q \cdot \cos^{2} \theta\).

In the injection-production well network, a source–sink injection–production unit is selected for streamline calculation. The cross-sectional schematic diagram of Ax is as Fig. 2:

Fig. 2
figure 2

Schematic diagram of Ax section

The blue curve in Fig. 2 is the characteristic streamline, which is used to represent the general direction of the streamline. The figure shows:

$$A_{x} = \, h\left( x \right) \cdot d\left( x \right)$$
(7)
$${\text{d(}}x{) = }\left\{ {\begin{array}{*{20}c} {{2}x \cdot {\text{tan}}\alpha {\text{ x}} \le \frac{l}{2}} \\ {2(l{ - }x{)} \cdot {\text{tan}}\alpha {\text{ x}} > \frac{l}{{2}} \, } \\ \end{array} } \right.$$
(8)

where α is the angle determined by the interwell position of the unit injection–production unit of the well pattern. l is the horizontal well spacing of one injection and one production, and when \(\, x{ = }\frac{l}{2}\), the potential of injection-production wells is equal, h(x) is the distance between the corresponding point and the bottom of the reservoir.

In order to solve h(x), XZ section is selected for analysis, as shown in Fig. 3:

Fig. 3
figure 3

Schematic diagram of XZ section

The blue curve in Fig. 3 is the characteristic streamline, which is used to represent the general direction of the streamline. According to the above picture, we can see:

$$h\left( x \right) = H - h_{0} \left( x \right)$$
(9)

where H is the overall reservoir thickness and h0(x) is the distance from the corresponding point to the top. Let Δh = H/N be the thickness of each layer and N be the total number of layers.

Then,

$$\theta_{{\text{i}}} = \frac{\pi }{2} - \arctan \left(\frac{{k_{i} - \frac{{\sum\limits_{i}^{N} {k_{i} } }}{N}}}{{\frac{{\sum\limits_{i}^{N} {k_{i} } }}{N}}}\right) - \arctan (\frac{l}{H})$$
(10)

Among them, ki is the permeability of the layer i reservoir.

$$x_{i} = \left\{ {\begin{array}{*{20}l} 0 & {i = 0} \\ {x_{i - 1} + \frac{\Delta h}{{\tan \theta_{i} }}} & {i \ne 0} \\ \end{array} } \right.$$
(11)
$$\begin{array}{*{20}c} {h(x) = H - (x - x_{i - 1} ) \cdot \tan \theta_{i} - \Delta h \cdot (i - 1)} & {x_{i - 1} } \\ \end{array} \le x < x_{i + 1}$$
(12)

For rhythm reservoirs, the shape function of an injection–production unit is represented by two shape functions h(x) and d(x). The determination of the shape function depends on the interface characteristics of the reservoir plane and longitudinal flow.

For the flow in the x direction,

$$v_{x} = \frac{\partial \psi }{{\partial y}}\frac{\partial \phi }{{\partial z}} - \frac{\partial \psi }{{\partial z}}\frac{\partial \phi }{{\partial y}} = \frac{q}{{A_{x} }}{\text{*cos}}^{2} \theta { = }\frac{{q_{{\text{x}}} }}{h(x)d(x)}$$
(13)

In order to simplify formula (13), decomposition equation right, let \(\psi ,\phi\) satisfy:

$$\frac{\partial \psi }{{\partial z}} = \frac{{q_{{\text{x}}} }}{d(x)},\frac{\partial \phi }{{\partial y}} = - \frac{1}{h(x)}$$
(14)

After integration:

$$\psi (x,z) = \frac{{q_{{\text{x}}} z}}{d(x)},\phi (x,y) = \frac{ - y}{{h(x)}}$$
(15)

The formula is only related to two parameters; we can get:

$$\frac{\partial \psi }{{\partial y}}{ = 0},\frac{\partial \phi }{{\partial z}}{ = 0}$$
(16)

According to formula (16) and seepage mechanics, simplify formula (4). Then, the velocity function satisfies:

Under \(x \le \frac{l}{2}\)

$$\left\{ {\begin{array}{*{20}l} {v_{x} = - \frac{\partial \psi }{{\partial z}}\frac{\partial \phi }{{\partial y}} = \frac{{q_{x} }}{{{2}x{\text{tan}}\alpha \left[ {H - (x - x_{i - 1} ) \cdot \tan \theta_{i} - \Delta h \cdot (i - 1)} \right]}}} \hfill \\ {v_{y} \left( {x,y} \right) = \frac{\partial \psi }{{\partial z}}\frac{\partial \phi }{{\partial x}} = \frac{{yq_{x} \tan \theta_{i} }}{{{2}x{\text{tan}}\alpha \left[ {H - (x - x_{i - 1} ) \cdot \tan \theta_{i} - \Delta h \cdot (i - 1)} \right]^{2} }}} \hfill \\ {v_{z} \left( {x,z} \right) = \frac{\partial \psi }{{\partial x}}\frac{\partial \phi }{{\partial y}} = \frac{{ - {\text{2tan}}\alpha zq_{x} }}{{\left[ {{2}x \cdot {\text{tan}}\alpha } \right]^{2} \left[ {H - (x - x_{i - 1} ) \cdot \tan \theta_{i} - \Delta h \cdot (i - 1)} \right]}}} \hfill \\ \end{array} } \right.$$
(17)

Under \({\text{x > }}\frac{l}{{2}} \,\)

$$\left\{ {\begin{array}{*{20}l} {v_{x} = - \frac{\partial \psi }{{\partial z}}\frac{\partial \phi }{{\partial y}} = \frac{{q_{x} }}{{2(l{\text{ - x)tan}}\alpha \left[ {H - (x - x_{i - 1} ) \cdot \tan \theta_{i} - \Delta h \cdot (i - 1)} \right]}}} \hfill \\ {v_{y} \left( {x,y} \right) = \frac{\partial \psi }{{\partial z}}\frac{\partial \phi }{{\partial x}} = \frac{{yq_{x} \tan \theta_{i} }}{{2(l{\text{ - x)}} \cdot {\text{tan}}\alpha \left[ {H - (x - x_{i - 1} ) \cdot \tan \theta_{i} - \Delta h \cdot (i - 1)} \right]^{2} }}} \hfill \\ {v_{z} \left( {x,z} \right) = \frac{\partial \psi }{{\partial x}}\frac{\partial \phi }{{\partial y}} = \frac{{{\text{2tan}}\alpha zq_{x} }}{{\left[ {2(l{\text{ - x)}} \cdot {\text{tan}}\alpha } \right]^{2} \left[ {H - (x - x_{i - 1} ) \cdot \tan \theta_{i} - \Delta h \cdot (i - 1)} \right]}}} \hfill \\ \end{array} } \right.$$
(18)

From Darcy's law, we can obtain the pressure distribution equation as follows.

Under \({\text{x}} \le \frac{l}{2}\)

$$\left\{ {\begin{array}{*{20}l} {\frac{\partial p}{{\partial x}} = - \frac{{\mu v_{x} }}{{K_{x} }} = - \frac{{\mu q_{x} }}{{{2}x \cdot {\text{tan}}\alpha \cdot K_{x} \left[ {H - (x - x_{i - 1} ) \cdot \tan \theta_{i} - \Delta h \cdot (i - 1)} \right]}}} \hfill \\ {\frac{\partial p}{{\partial y}} = - \frac{{\mu v_{y} \left( {x,y} \right)}}{{K_{y} }} = \frac{{yq_{x} \tan \theta_{i} \mu }}{{{2}x{\text{tan}}\alpha K_{y} \left[ {H - (x - x_{i - 1} ) \cdot \tan \theta_{i} - \Delta h \cdot (i - 1)} \right]^{2} }}} \hfill \\ {\frac{\partial p}{{\partial z}} = - \frac{{\mu v_{z} \left( {x,z} \right)}}{{K_{z} }} = \frac{{ - {\text{2tan}}\alpha zq_{x} \mu }}{{K_{z} \left[ {{2}x \cdot {\text{tan}}\alpha } \right]^{2} \left[ {H - (x - x_{i - 1} ) \cdot \tan \theta_{i} - \Delta h \cdot (i - 1)} \right]}} \pm \rho g} \hfill \\ \end{array} } \right.$$
(19)

Under \(x{ > }\frac{l}{{2}} \,\)

$$\left\{ {\begin{array}{*{20}l} {\frac{\partial p}{{\partial x}} = - \frac{{\mu v_{x} }}{{K_{x} }} = - \frac{{\mu q_{x} }}{{{2}(l{\text{ - x)tan}}\alpha K_{x} \left[ {H - (x - x_{i - 1} ) \cdot \tan \theta_{i} - \Delta h \cdot (i - 1)} \right]}}} \hfill \\ {\frac{\partial p}{{\partial y}} = - \frac{{\mu v_{y} \left( {x,y} \right)}}{{K_{y} }} = \frac{{yq_{x} \tan \theta_{i} \mu }}{{2(l{\text{ - x)}} \cdot {\text{tan}}\alpha K_{y} \left[ {H - (x - x_{i - 1} ) \cdot \tan \theta_{i} - \Delta h \cdot (i - 1)} \right]^{2} }}} \hfill \\ {\frac{\partial p}{{\partial z}} = - \frac{{\mu v_{z} \left( {x,z} \right)}}{{K_{z} }} = \frac{{{\text{2tan}}\alpha zq_{x} \mu }}{{K_{z} \left[ {2(l{\text{ - x)}} \cdot {\text{tan}}\alpha } \right]^{2} \left[ {H - (x - x_{i - 1} ) \cdot \tan \theta_{i} - \Delta h \cdot (i - 1)} \right]}} \pm \rho g} \hfill \\ \end{array} } \right.$$
(20)

For different heterogeneous conditions, the three-dimensional flow characteristics are controlled by the setting of boundary conditions in order to realize the streamline characterization of three-dimensional space.

Reservoir streamlines and saturation distribution under rhythmic conditions

According to formula, given the parameters of the shape function, combined with the given q, the three-dimensional flow function model can be established using MATLAB software to construct the flow equation of the rhythm reservoir. Combined with the three-dimensional potential distribution of the reservoir constructed by the three-dimensional potential function model, the velocity and streamline distribution of the reservoir are obtained. Because the main feature of the prosody condition is reflected in the Z direction, the distribution of the reservoir streamlines in the XZ direction is drawn according to the velocity field distribution: follow the water flooding front, determine the specific position of the front through seepage velocity and integrate the saturation. Then, the flow units in the reservoir are divided according to the EWU theory through the streamline distribution and velocity distribution and then obtain the saturation change at different times on the streamline. Finally, the overall water cut and recovery ratio of the entire flow unit are obtained, and the corresponding recommended measures are given in combination with the above information.

Positive rhythm reservoir

The design reservoir is a three-layer positive rhythm reservoir. The permeability ratios are set to 2, 4, 8 and 16, respectively. The permeability of the three layers is 100mD, 50mD and 25mD from bottom to top. We calculate the reservoir streamline and saturation distribution and then calculate the variation law of reservoir water cut with time and production degree under different conditions. Relative permeability of reservoir is shown in Table 2:

In order to study the influence of different heterogeneous conditions on the water cut and saturation of the reservoir, the entire development process is divided into four stages as shown in Fig. 4 according to the difference in the water cut of the reservoir:

Fig. 4
figure 4

Water cut growth curve

Stage 1-initial water flooding stage: This stage is the most efficient stage of oilfield development. The water flooding front has not reached the oil well. The water cut of the reservoir is low at this stage, generally below 10%.

Stage 2-rapid rise of water cut stage: This stage is a breakthrough in the front of water flooding. The oil well begins to contain water in a large area. At this stage, the displacement efficiency begins to decrease, and the water cut rises rapidly, and the water cut is < 80%.

Stage 3-the water cut slowly rising stage: 95% > water cut > 80% in this stage, the oil displacement efficiency is poor at this stage, the reservoir begins to appear high water cut, and the water cut increase rate decreases.

Stage 4-High water cut stage: In this stage, the water cut is > 95%, most of the reservoir is flooded, the oil displacement efficiency is < 5%, and the remaining oil saturation of the reservoir is basically unchanged.

The division and characteristics of effective flow units in the corresponding stages are shown in Fig. 5.

Fig. 5
figure 5

Schematic diagram of effective drive unit characterization at 1–4 stages

The I–IV areas in Fig. 5 correspond to the main existing areas of I-IV EWU in Table 1. In stages 1 and 2, the invalid water injection cycle has not yet formed, and there only exist classes II (red streamline—the lowermost streamline), III (an area not covered by a streamline) and IV (red streamline—the uppermost streamline) unit. There are only high-speed flow effective areas and low-speed flow effective areas in the flow area between injection and production units; in the stage 3, when the main flow line is flooded, the class I EWU is the high-speed flow ineffective area (yellow flow line—the lowermost streamline). Before the comprehensive water cut rises to 98%, four class of EWUs exist simultaneously and have obvious characteristics. The ineffective water injection cycle area (class I) will expand significantly over time. In stage 4, when the comprehensive water cut reaches 98%, the class II flow unit (the effective area of high-speed flow) disappears. Currently, there only exist class I (yellow streamline—the lowermost streamline), III (an area not covered by a streamline) and IV (yellow streamline—the uppermost streamline) unit. At this time, the remaining oil in the flow area is mainly concentrated in the class IV EWU.

Table 1 EWU division and characteristics (Du et al. 2017)
Table 2 Relative permeability of reservoir

The blue curve in Fig. 6 is the reservoir streamline. In order to better analyze, I and II areas are marked in Fig. 6. The density of streamlines is used as a standard for dividing high-speed and low-speed units, and the distribution of water cut is used as a simple standard for distinguishing between effective and ineffective units (Du et al. 2017). Observation area I can be seen that the streamline is dense in the lower part of the reservoir, indicating that the flow velocity in this part is faster. Combined with the distribution of water cut, it can be found that the early stage is high-speed flow effective flooding and the later stage is high-speed flow invalid flooding. The area II reservoir upper streamline is relatively sparse, indicating that the flow velocity is slow and the area is small, and the oilfield development is low. It belongs to low-speed flow effective flooding, and later turns into low-speed flow invalid flooding. The rhythm reservoir is affected by the uneven permeability rate, and the lower permeability is large. At the same time, due to the gravity effect, this distribution is formed, and the lower part is easy to form a dominant channel, which ultimately leads to a significant reaction in development efficiency.

Fig. 6
figure 6

Streamline distribution of positive rhythm reservoir.

It can be seen from Fig. 7 that for positive rhythm reservoirs, the ultimate recovery of the reservoir is higher than 51%, the greater the permeability difference, the smaller the reservoir recovery rate, because the greater the permeability difference, the greater the storage capacity. The stronger the vertical heterogeneity of the layer is, the easier it is to form dominant permeability channels in the lower part of the reservoir. The greater the difference, the greater the permeability in the lower part, which is small, ultimately resulting in a low recovery factor.

Fig. 7
figure 7

Distribution of recovery ratios of reservoirs with positive rhythm at different permeability ratios

Anti-rhythm reservoir

The designed anti-rhythm reservoir has the same permeability ratios as those of the positive rhythm. The permeability of the three layers is 25mD, 50mD and 100mD from bottom to top. We calculate the reservoir streamline and saturation distribution and then calculate the variation law of reservoir water cut with time and production degree under different conditions.

The division and characteristics of EWUs in the corresponding stages are shown in Fig. 8.

Fig. 8
figure 8

Schematic diagram of effective drive unit characterization at 1–4 stages

The I-IV areas in Fig. 8 correspond to the main existing areas of I-IV EWU in Table 1. In stages 1 and 2, ineffective injection cycles have not yet been established, with only class II (red streamlines—top streamlines), III (areas not covered by streamlines) and IV (red streamlines—bottom streamlines) elements. In the flow area between the injection and the production units, there are only high-speed flow effective area and low-speed flow effective area; in stage 3, when the main flow line is flooded, the class I EWU is the high-speed flow ineffective area (yellow flow line—the uppermost streamline). In stage 4, when the comprehensive water cut reaches 98%, the class II flow unit (the effective area of high-speed flow) disappears. The remaining oil in the flow area is mainly concentrated in the class IV EWU.

Compared with positive rhythm reservoirs, the remaining oil is mainly located in the class IV EWU. It can be seen that determining the location of the class IV EWU is a feasible and efficient method when determining the main areas of remaining oil.

The blue curve in Fig. 9 is the reservoir streamline. In order to better analyze, I and II areas are marked in Fig. 6. From observation area I, we can see that the upper stream line of the reservoir is denser, indicating that the flow velocity in this part is faster. Combining with the distribution of water cut, it was class II EWU in the early stage and become class I EWU in the later stage. The lower streamlines of the oil reservoirs in area II are relatively sparse, indicating that the flow velocity is slow and the degree of oilfield development is low. Even under the influence of gravity, the degree of bottom development is still very low. It can be seen that permeability is a more important factor, which determines the degree of reservoir exploitation.

Fig. 9
figure 9

Distribution of anti-prosody reservoir streamlines.

It can be seen from Fig. 10 that compared with the positive rhythm reservoir, the overall water flooding effect of the reservoir is better due to gravity, and the reservoir time is later.

Fig. 10
figure 10

Comparison of water cut curves of positive and anti-rhythmic reservoirs (permeability ratio = 2)

It can be seen from Fig. 11 that for reverse rhythm reservoirs, the ultimate recovery rate of the reservoir is higher than 56%. The greater the difference in permeability, the smaller the recovery rate of the reservoir. Compared with positive rhythm reservoirs, because of gravity, the formation of dominant channels in the upper part of the reservoir takes longer and the flooding displacement has a wider range, so the overall recovery rate is higher than that of positive rhythm reservoirs.

Fig. 11
figure 11

Distribution of recovery ratios of extremely reverse anti-rhythmic reservoirs with different permeability ratios

Influence of injection-production conditions

In order to understand and analyze the form and distribution of residual oil in rhythmic reservoir more accurately, we study the occurrence state of remaining oil by designing different injection production combinations and provide guidance for further design of potential tapping method of remaining oil in rhythmic reservoir.

Zonal injection-commingled production

In this paper, three layers of positive rhythm reservoir are designed, then each layer is injected with water and the middle layer is produced. Figure 12a shows the distribution of reservoir saturation after 1500 days of production, and Fig. 12b, c, d shows the distribution of streamline and remaining oil saturation at different times.

Fig. 12
figure 12

a Distribution of recovery ratios of reservoirs with positive rhythm at different permeability ratios. b 300 days, c 800 days, d 1500 days

As can be seen from Fig. 12a, b, c, d, type II EWU of positive rhythm reservoir is firstly generated in the lowest layer, and the type IV EWU is formed in the middle layer and upper layer. With the development, water is first seen in the bottom layer of the reservoir, which means the form a dominant channel and the loss of type II EWU. With further water injection, water appears in the middle layer gradually, and the type IV EWU decreases gradually. With the joint influence of gravity and other functions, the remaining oil range of the last time gradually tends to be stable, and the overall water content of the reservoir is gradually stable.

Commingled injection-zonal production

In this paper, three layers of positive rhythm reservoir are designed, and then, the middle layer is injected, all three layers are produced. Figure 13a shows the distribution of reservoir saturation after 1500 days of production, and Fig. 13b, c, d shows the distribution of streamline and remaining oil saturation at different times.

Fig. 13
figure 13

a Distribution of recovery ratios of reservoirs with positive rhythm at different permeability ratios. b 300 days, c 800 days, d 1500 days

As shown in Fig. 13a, c, d, type II EWU of positive rhythm reservoir is firstly generated in the lowest layer, and the type IV EWU is formed in the middle layer, and the upper reservoir is basically not driven. With the development, water is first seen in the bottom layer of the reservoir, which means the form of I EWU and the loss of type II EWU. With further water injection, water appears in the middle layer gradually, and the type IV EWU decreases gradually.

Figure 14 is the water cut of reservoir under different injection production conditions, which can show that the increase in water cut in zonal injection-commingled production model is slower than that in commingled injection-zonal production model, and the effective displacement time is longer. Comparing the two mining methods, the final recovery of zonal injection-commingled production mode is higher and the interlayer interference is relatively small, which makes the upper reservoir in the effective driving range. When the reservoir reaches the stage of high water cut, the remaining oil can be extracted by turning off the high permeability layer and increasing the water injection rate of the low permeability layer. The reservoir produced by the commingled injection-zonal production can increase the sweep area of water injection by injecting plugging agent into the high permeability layer.

Fig. 14
figure 14

Water cut of reservoir under different injection production conditions

Recommended measures

For the positive rhythm reservoir, the key to improve the effect of water injection development is to seal the high permeability section at the bottom of sandstone and to fracture the low permeability zone at the top, so as to promote the displacement of injected water along the low permeability zone at the top and increase the sweep volume. The saturation distribution before and after adjustment is shown in Fig. 15.

Fig. 15
figure 15

Saturation distribution of reservoirs with positive rhythm

For the anti-rhythm reservoir, the main idea of improving oil recovery is flow field transposition. Two main measures are taken. On the one hand, perforation is used to increase the water injection pressure in the low-permeability area at the bottom to expand the water flooding area; on the other hand, the top area is properly sealed to slow the formation of dominant seepage channels. The saturation distribution before and after adjustment is shown in Fig. 16.

Fig. 16
figure 16

Saturation distribution of reservoirs with anti-rhythm

Conclusion

In this paper, a theoretical model of the three-dimensional EWU is established through the heterogeneous shape function and the principle of superposition of hydrodynamic potential; it is applied to the reservoir under the pros and cons of rhythm. The corresponding streamline and saturation distribution are obtained and then divide the EWU according to the distribution. Then, the location of the effective displacement unit of low-speed flow is determined, which is the main enrichment area of remaining oil, and corresponding recommended measurements are given for the subsequent remaining oil production in the oil field.

The three-dimensional EWU theoretical model is based on the two-dimensional EWU theory (Nan et al. 2016), because the flow function is defined by the incompressible fluid under two-dimensional plane conditions. In order to solve the three-dimensional flow problem of heterogeneous reservoirs, an orthogonal coordinate system is used in this paper to characterize the heterogeneous distribution characteristics. The characteristic of the three-dimensional streamline is the spatial streamline obtained by the intersection of two surfaces in the coordinate system. According to the principle of hydrodynamic potential superposition, the reservoir geological conditions are expressed in the form of potential, superimposed with the seepage potential, and then, the boundary conditions are determined according to the heterogeneous shape function to complete the expression of the three-dimensional flow function model. The corresponding water cut and saturation can be predicted by the presented model. And when the remaining oil enrichment area through the obtained streamline distribution is determined, the production measures can be recommended correspondingly.

For the positive rhythm reservoir, the lower streamline in the flow unit is dense, indicating that the area has a fast waterflood speed and is a high-speed flow unit. It is the main area where the dominant channel is formed. The front edge of the waterflood breaks through the lower part of the reservoir as the rate difference increases, and the reservoir seepage time becomes longer. The remaining oil is mainly distributed in the class IV EWU, that is, the upper part of the reservoir. The ultimate recovery rate of the reservoir under different permeability ratios is about 55%. Compared with positive rhythm reservoirs, due to the gravity effect, the formation of dominant channels in the upper part of the reservoir takes longer time and the flooding displacement has a wider range. The remaining oil is mainly distributed in the class IV EWU that is, the bottom part of the reservoir. The ultimate recovery rate of the reservoir under different permeability ratios is about 62%. The overall recovery rate is 7% higher than that of positive rhythm reservoirs.

For rhythmic reservoirs, the dominant channel of seepage flow is located in the class I EWU, and the remaining oil is located in the class IV EWU. After the reservoir reaches the high water-cut stage, we can adopt the method of plugging the dominant channel to expand the effective displacement range of the reservoir. This will increase the overall sweep efficiency of water injection development, reduce the damage to the reservoir caused by ineffective cycles and submit the final recovery factor.