Capillary action plays an important role in oil recovery by water flooding. As the pore channel radius decreases, the capillary action increases, which seriously affects reservoir development, especially in a low-permeability sandstone reservoir. The Ordos Basin is a typical low-permeability sandstone reservoir in China. Studying how variations in the capillary force affect the remaining oil production on the pore scale helps in understanding how the capillary action improves the development of unconventional reservoirs. In this study, the core of the Chang 6 Formation in the Ansai Oilfield, Ordos Basin was scanned by computed tomography. Then, the digital core model was established. The oil–water two-phase flow in pores was described using the method based on the Navier–Stokes equation coupled with the method of the volume of fluid simulation. The water flooding process was simulated on the pore scale. The results show that in the process of pore scale water flooding, the oil–water interface stays at the position between the throat channel and the pore area, where the oil–water interface reverses and the capillary force presents resistance, forming the capillary barrier or capillary valve. Affected by the capillary barrier, the oil–water two-phase flow in the process of water flooding is described by a “step-by-step” model. The pore structure characteristics at the junction of the pore area and the throat channel control the movement of the oil–water interface and affect the water flooding production and the ultimate recovery factor. As the liquid injection rate increases, the oil on both sides of the main channel is produced. While the oil recovery rate reaches 66%, the remaining oil on the edges becomes increasingly difficult to be produced. This difficulty is closely related to the viscosity of the injection fluid, interfacial tension, injection rate, pore radius, and pore wall wettability.

Waterflood development is the most widely used secondary oil recovery method worldwide, which can keep the reservoir pressure stable and improve oil displacement efficiency. However, the water cut increases continuously with waterflood development and the efficient development of the remaining oil in the stratum becomes a key scientific issue [1]. Although measures such as profile modification, water plugging, and well pattern adjustment can be taken to slow down the rise of water cut [24], the water cut still increases. Additionally, most reservoirs in China are continental clastic deposits, accounting for about 92%. They are highly heterogeneous and demonstrate high crude oil viscosity in the longitudinal and horizontal directions, resulting in low oil recovery by water flooding [1]. Currently, most reservoirs in China have entered the mature development stage, with water cuts higher than 80% and recovery rates lower than 30% [5, 6]. Therefore, using existing scientific and technological methods, the industry has focused on how to maximize the efficiency of oil displacement in a high water cut period. The occurrence state and distribution characteristics of the remaining oil in the high water cut stage are closely related to the pore properties of the reservoir. The influence of microscopic forces such as surface tension and capillary force on the remaining oil production becomes particularly important [710]. Particularly in unconventional reservoirs represented by low-permeability reservoirs, coalbed methane reservoirs, and shale reservoirs, the scale of pore and throat channel is small, and the capillary force has a significant impact on the occurrence and production of reservoir fluids [11, 12]. A deep understanding of the capillary action could greatly help improve the development effect of unconventional reservoirs [13, 14].

Currently, the microflow process of water flooding is described by physical and numerical simulations [1518]. The physical simulations can only derive the distribution law of oil and water in pore space and the characteristics of macroscopic flow and pressure. Owing to measurement limitations, physical simulations are limited in the ability to study the mechanism behind the macrocharacteristics of water flooding. Numerical simulations can quantitatively describe the multiphase flow process in pores and reproduce the movement disciplinarian of oil and water at the pore scale. Over the past few years, numerical simulations have gradually become the new method for studying multiphase flows at the pore scale. To further investigate how the oil–water interface affects the oil–water two-phase flow on the pore scale, it is crucial to accurately determine the position and tension of the oil–water interface. The research methods for interface position tracking include the front tracking method, level set method, and volume of the fluid method [1921]. The volume of fluid method can be extended from the two to the three dimensions, the volume is conserved in the calculation process, and the interface movement can be accurately tracked by a high-precision solution algorithm or accurate interface reconstruction technology. Open-source Field Operation and Manipulation (OpenFOAM) is a C++ code package that is open source, readable, and expansible [22]. There are many two-phase solvers in OpenFOAM. interFoam is used to solve the insoluble two-phase movement process with a free surface. The Navier–Stokes equation is used to describe the fluid movement, the volume of fluid method is used to solve the position of the two-phase free surface, and the finite volume method is the most perfect and basic two-phase solver in OpenFOAM.

The Ordos Basin, an important oil-bearing basin in the west of the North China Craton, is a typical representative of continental low-permeability reservoirs in China [2326]. Low porosity and low permeability are the significant characteristics of sandstone reservoirs in the Yanchang Formation. The porosity ranges mainly between 6% and 16%, while the permeability ranges mainly between 0.1 and 3.0 mD [27, 28]. Understanding the intrinsic relationship between capillary action and the characteristics of water flooding for low-permeability reservoirs and providing a theoretical basis for analyzing the formation mechanism and distribution law of the remaining oil is very important.

In this study, the images of the core of the T3Ch6 Formation in the Ansai Oilfield, Ordos Basin from computed tomography (CT) scans were used to construct the digital core model to simulate the core structure. Based on the interFoam solver in OpenFOAM, the calculation process uses the Navier–Stokes equation to describe the two-phase medium movement and the volume of fluid method to track the position of the two-phase interface. The pressure implicit with splitting of operators (PISO) algorithm [29] was used to decouple the pressure and velocity equations and simulate the two-phase medium movement. The dynamic positions of oil and water during water flooding with pore size are simulated, and the changes in the capillary force and the influence of the capillary barrier on oil–water movement are analyzed. The analysis of the remaining oil production by increasing the liquid injection rate is also discussed.

The Ordos Basin is the second largest sedimentary basin in China, which is located in the Yinshan Mountain in the north, the Qinling Orogenic Belt in the south, the Lvliangshan Orogenic Belt in the east, and the Liupanshan Thrust Belt in the west (Figure 1). It is currently composed of six first-order tectonic units, including the Yimeng Uplift, the Weibei Uplift, the Jinxi Tectonic Belt, the Yishan Slope, the Tianhuan Depression, and the Western Thrust Zone [30]. It is adjacent to the Hetao basin, the Yinchuan Graben, and the Weihe Basin. The development of the Ordos Basin in the Paleo-Mesozoic can be divided into three evolutionary stages: the divergent boundary craton basin from Cambrian to Early Ordovician, the convergent boundary craton basin from Middle Ordovician to Middle Triassic, and the residual craton basin in the plate from Late Triassic to Early Cretaceous [28, 31]. The present basin is a residual basin that has undergone multiple transformation stages. The basin basement is the Precambrian crystalline metamorphic rock. The Upper Proterozoic–Lower Paleozoic marine carbonate deposition, the Upper Paleozoic marine–continental coal-bearing clastic deposition, the Mesozoic continental fluvial–delta–lacustrine facies deposition, the Cenozoic aeolian loess, and the fluvial facies deposition are deposited on the basement [32]. The total thickness is about 5000–6000 m. The oil-bearing strata in the Ordos Basin are the Yanchang Formation of the Middle–Upper Triassic and the Middle Jurassic Formation, and the gas-producing strata are the Carboniferous–Permian and Ordovician Formations [31]. The Triassic Yanchang Formation is rich in low-permeability reservoirs, and the proven tight oil reserves are about 2.0 billion tons [28]. It mainly occurs in tight sandstone reservoirs of the T3Ch6-8 oil-bearing intervals nearby or is interbedded with oil shale layers, without migrating to large distances.

The Ansai Oilfield is located in the middle east of the Yishan Slope in the Ordos Basin, and its tectonic setting is a gentle west-dipping monocline [33]. The Ansai Oilfield is a typical low-permeability and low-yield oilfield. Its reservoir is the Yanchang Formation of the Triassic, buried at a depth of 1000–1300 m, and it is a sedimentary system dominated by the delta of an inland freshwater lake [34]. The T3Ch6 oil-bearing formation is an important oil-bearing formation in the Wangyao area, which can be divided into T3Ch61, T3Ch62, and T3Ch63 layers. T3Ch61 is the main production layer. The average porosity of the T3Ch6 reservoir is 13.30%, the average permeability is 2.3 mD, and the average thickness of the oil layer is about 12 m. These data suggest that the reservoir demonstrates low porosity, low permeability, and low production [33].

3.1. Pore Scale Numerical Model

In this study, the volume of fluid method is used to track the position of the oil–water interface. The principle is to define a volume fraction α. When α is equal to 0 or 1, it means that the grid cell is occupied by the specified phase fluid or no specified fluid. When α ranges between 0 and 1, the grid cell has a free interface, which is called the interface element. The free surface was determined by tracking the interface position and direction with the volume fraction. In this study, the grid cell is filled with water when α is equal to 1, and the grid cell is filled with oil when α is equal to 0. The pore scale dynamic model of the oil–water two-phase flow in porous media can be constructed by the following equations:
  • (1)

    Continuity equation:

(1)u=0,
where u is the velocity, m s-1
  • (2)

    Volume fraction equation:

(2)αt+αu=0,
where α is the volume fraction
  • (3)

    Momentum equation:

(3)ρut+ρuuμτ=p+ρg,
where ρ is the fluid density, kg m-3; μ is the dynamic viscosity coefficient, Pa s; P is the dynamic pressure, Pa; and τ is the rate of the strain tensor, s−1
In Equation (3), ρ is the average density of the two-phase mixed fluid, μ is the average dynamic viscosity coefficient of the two-phase fluid, and their expressions are the following:
(4)ρ=αρ1+1αρ2,(5)μ=αμ1+1αμ2,
where ρ1 is the density of phase 1, kg m-3;ρ2 is the density of phase 2, kg m-3; μ1 is the dynamic viscosity coefficient of phase 1, Pa S;μ2 is the dynamic viscosity coefficient of phase 2, Pa S; α is the volume fraction of phase 1, which is explicitly tracked by Equation (2) in the calculation process.

3.2. Solving Algorithm of the Pore Scale Numerical Model

3.2.1. Finite Volume Discretization Method

The basic idea of computational fluid dynamics is to discretize continuous partial differential equations in grid cells and then solve the simplified algebraic equations to obtain the relevant variables. The commonly used numerical discretization methods include the finite volume method, finite difference method, and finite element method. The finite volume method satisfies the law of conservation of mass, is accurate, and can handle the meshing problem of the complex boundary area. Therefore, using the finite volume discretization method, this study discretized the partial differential equations composed of Equations (1), (2), and (3) and obtained the linear algebraic equations.

In the discretization process of the equation using the finite volume method, for any partial differential equation,
(6)φt+uφDφφ=Sφ,
where φ is relevant physical quantity; Dφis the diffusion coefficient; t is the time, s; u is the velocity, m s-1; and Sφ is a source item.
By integrating Equation (6), we can get
(7)VφtdV+VuφdVVDφφdV=VSφdV.
According to Gauss’s theorem,
(8)φtV+SuφdSSDφφ=VSpφ+Su,
where V is the volume of the control unit, SP is the linearization coefficient of the source item Sφ, and Su is the linearization constant of the source item Sφ.
The finite volume discretization method divides the calculation area into multiple control units containing multiple surfaces. Each surface is composed of a sequence of points. The physical and geometric properties of the surfaces are obtained through a series of point calculations. After discretizing the control units, the following results are obtained:
(9)φpn+1φpnΔtV+fufSfφffDφφfSf=VSpφ+Su,
where φPn+1 is the value of a physical quantity at time n+1, φPn is the value of a physical quantity at time n, uf is the speed on the grid cell surface, φf is the value of a physical quantity in the element surface, and Sfis the area vector of the element face.
The convection term in Equation (9), namely, the second term on the left, is interpolated:
(10)fufSfφf=ψffWfφpn+1+1Wfφneighbour,fn+1,(11)ψf=ufSf,(12)φf=Wfφpn+1+1Wfφneighbour,fn+1,
where ψf is the flow rate on the element surface, Wfis the weight coefficient of linear interpolation, and φneighbour,fn+1 is the value of a physical quantity on the adjacent element at time n+1.
The diffusion term in Equation (9), namely, the third term on the left, is interpolated as follows:
(13)fDφφfSf=fDφΔφn+1rneighbour,frprneighbour,frp2Sf=fΔφn+1βf,(14)Δφn+1=φneighbour,fn+1φpn+1,(15)βf=SfDφrneighbour,frpΔφn+1,
where rneighbour,f is the position of the face center of the adjacent unit, rp is the location of the cell center, Δφn+1 is the difference between a physical quantity on the adjacent element surface and the core of the element body at time n+1; and βf is a physical quantity obtained by linear interpolation, which can be obtained by linear interpolation.
Substituting Equations (10) and (13) into Equation (9), the formula is obtained as follows:
(16)φpn+1φpnVΔt+ψffφffΔφn+1βf=VSpφ+Su.
Further, simplification Equation (16) results in the following:
(17)apφpn+1=faN,fφneighbour,fn+1+VΔtφpn+VSu,(18)ap=VΔt+fψfWffβfVSp,(19)aN,f=ψf1Wf+βf,
where αp is the linearization coefficient and αN,f is the linearization coefficient.
Write the organized equation in the form of a system of linear equations:
(20)apφpn+1=faN,fφneighbour,fn+1+b,(21)b=VΔtφpn+VSu,
where b is the linearization constant.

3.2.2. Solver Algorithm and Process

There are currently two types of numerical solutions: coupled and separated. The separated method offers a simple and convenient way to calculate solutions and is widely used. The most widely used method for solving fluid problems is the pressure correction method, which includes SIMPLE, SIMPLEC, PISO, and other algorithms [29, 35]. The PISO algorithm can greatly reduce the number of iterations of the unsteady state problem and further improve calculation efficiency. In this study, the PISO algorithm is used to decouple the rate and pressure.

According to the finite volume discretization process of the polyhedral mesh, the momentum Equation (3) is discretized, and the derived semidiscretization scheme is the following:
(22)apu=AHpdpc,
whereαp is the upper diagonal value of the coefficient matrix obtained by discretizing Equation (3), Pd is dynamic pressure, Pc is the pressure gradient owing to surface tension, and AH is presented as
(23)AH=NaNu+b,
where αN is the implicit contribution coefficient of the node next to the current node in the discretization process of Equation (3) and b represents the explicit discrete contributions other than pressure.
Divide both sides of Equation (23) by ap and arrange the equation as follows:
(24)u=AHappdappcap.
The velocity field obtained from Equation (24) should satisfy the continuity equation. Substitute Equation (24) into the continuity Equation (1) and get
(25)1appd=1apAH1appc.
The pressure Equation (25) is derived by the PISO algorithm, and the new pressure obtained by solving this equation is substituted into Equation (24) to update the velocity. Additionally, to facilitate the diffusion of convective terms the next time, the flow rate at the interface of the element should be updated at the same time as the velocity at the center of the body. The surface flow rate is updated by the interpolation of Equation (24) on the surface:
(26)φf=AHapfSfpdapfSfpcapfSf,
where φf is surface flow rate.

In the numerical simulation process of this study, the flow of calculating the solutions of equations is as follows: First, the finite volume method is used for the discrete solution of the momentum Equation (3), and the velocity u is derived. Then, the PISO algorithm is used to derive the pressure Poisson Equation (25), and the volume fraction Equation (2) is solved to obtain the volume fraction at the next moment. Finally, Equations (4) and (5) are used to track the fluid properties. Then, repeat the first step until the end of the cycle is reached or convergence accuracy is achieved.

3.3. Digital Core Model and Parameters

The core sample was drilled in the Ansai Oilfield of the Ordos Basin, whose lithology is sandstone. A cylinder with a diameter of 2.00 mm is cut in the selected area on the core column, and the sample is scanned in a micron CT scanner with a scanning resolution of 0.70 μm.

The micron CT used in this study is the Zeiss Xradia 510 Versa micron CT scanner produced by the Zeiss company. The 3D spatial resolution is less than 0.70 μm, the pixel size of the 40x detector (objective lens) at the maximum magnification is 0.34 μm, and the pixel size of the 20x detector (objective lens) ranges from 0.65 to 0.68 μm. The moving range of the X-ray source is 215.00 mm, and the moving range of the detector is 290.00 mm. Simultaneously, built-in in situ tensile or compression test bench, equipment model DEBEN-CT5000-RT-SS, maximum support 5 kN loading force, and loading speed from 0.03 to 0.30 mm/min. The sample was scanned by a micron CT scanner. The scan results are shown in Figure 2(a). The cross-section obtained by scanning indicates that the clastic particles have medium-poor sorting, poor roundness, and general pore connectivity, belong to a grain-support structure. The local quartz particles are sutured contacts. According to the scanned data, image segmentation is performed to obtain the image binary map (Figure 2(b)). Based on this, the pore network model can be established to calculate pore throat parameters and perform seepage flow simulations.

In this study, the digital core model used in the test was selected based on the CT scan data (Figure 2(c)). The physical size of the digital core is 1.44mm×1.44mm, with particles in the white area and pore channels in the gray area. Initially, the core is filled with oil. During the simulation, water is injected into the pore from A on the left, displacing oil in the core, and liquid is drained from B. All other positions are closed. The boundary conditions used in the calculation are shown in Table 1, and the physical conditions of the fluid are shown in Table 2.

Water was injected into the core from inlet A at a constant injection rate, and the distribution characteristics of oil and water during water flooding at this injection rate were obtained (Figure 3). The injection rate used in this core model is 0.005 m/s, the capillary number is 7.15×105, the oil–water viscosity ratio is 5 : 1, and the wetting angle is 45°.

Figure 3 shows the oil–water distribution and the positions of the oil–water interface at (a) 0.180 s, (b) 0.230 s, (c) 0.325 s, (d) 0.400 s, (e) 0.465 s, (f) 0.585 s, (g) 0.775 s, (h) 0.820 s, and (i) final. The oil–water interface stays between the throat channel and the pore area, where the throat channel or pore radius changes abruptly. For example, at 0.180 s, the oil–water interface stays at the positions of a1, b1, and c1 (Figure 3(a)); at 0.230 s, it stays at position b2 (Figure 3(b)); at 0.325 s, it stays at positions b3 and e1 (Figure 3(c)); at 0.4 s, it stays at positions e2 and e3 (Figure 3(d)); at 0.465 s, it stays at positions b4 and b5 (Figure 3(e)); at 0.585 s, it stays at position b6 (Figure 3(f)); at 0.775 s, it stays at positions e5 and b7 (Figure 3(g)); at 0.820 s, it stays at positions b8 and b9 (Figure 3(h)). In this case, the capillary tension may appear as resistance even under wet conditions.

5.1. Variation of Capillary Force in Pore Scale during Water Flooding

During water flooding, the capillary force greatly affects oil and water distribution [15, 18, 36, 37]. The oil–water interface is easy to stop in the area where the radius of the throat channel and pore area changes abruptly (Figure 3).

To explain this phenomenon, we present the morphology of the oil–water interface at different positions and the variation in the capillary force at the corresponding positions (Figure 4). When the oil–water interface is in the throat channel, the direction of the capillary force (pointing to the side of the concave liquid surface) under water-wet conditions is consistent with the direction of the oil–water movement and the capillary force is the power (A ⟶ B). When water flows out of the throat channel and θ+β<π/2 (θ is the three-phase contact antenna of the oil–water–pore, β is the opening angle of the pore), the direction of the oil–water movement remains the same (B ⟶ D), but the capillary force gradually decreases from Pcmax (Figures 4(a) and 4(b)). When θ+β>π/2, the oil–water two-phase interface reverses, and the capillary force points to the waterside, preventing oil and water movement (B ⟵ D) (Figures 4(c) and 4(d)). After the oil–water interface moves to position B, the oil–water interface only deforms but does not move forward (the three-phase contact line stays at position B) until the angle between the deformed interface and the waterside angle of pore wall equals the equilibrium wetting angle θ (the oil–water morphology b is shown in Figure 4(c)). During this process, the capillary force formed on the pore wall of oil–water changes from Pcmax to 0 (the oil–water morphology d is shown in Figure 4(c)); thereafter, the oil–water interface reverses and the capillary force becomes negative until the negative maximum Pcmin. When the fluid power is insufficient to overcome this maximum resistance, the oil–water interface stops moving. Therefore, the driving pressure must be higher than the negative maximum Pcmin to promote the fluid flow. This threshold is called the capillary valve.

Figure 3 shows that water flows into the pore from the left inlet. At 0.180 s, the oil–water interface encounters the capillary barrier at position b and stops moving (Figure 3(a)). The effect of the capillary valve is formed here, which causes the lateral spread of oil and water. At 0.230 s, the effect of the capillary valve is formed at positions a and c, preventing further movement of oil and water (Figure 3(b)). After the oil–water interface at positions a and c is blocked, the upstream pressure increases, causing the capillary valve at b to be broken. As a result, the oil–water interface continues to migrate to d and stops at positions f and e at 0.325 s (Figure 3(c)). The blockage at position f causes the oil–water interface to migrate downward and stop at positions g and h at 0.4 s (Figure 3(d)). Simultaneously, water moves further along the center channel and blocks at positions i and h at 0.465 s (Figure 3(e)). This blockage will cause the capillary valves at g and j to be broken and the oil in the fj channel to be displaced (Figure 3(f)). Then, the water spreads forward along the center at position j and blocks at positions k, g, and j (Figure 3(g)). This lateral blockage will cause the capillary valve at position j of the main channel to break and block at positions l and m at 0.820 s (Figure 3(h)). Finally, the oil–water interface moves forward along the channel where m is located and flows out along the outlet.

5.2. Capillary Barrier and the “Step-by-Step” Model of Oil–Water Movement

During water flooding, the oil–water movement is not stable and continuous, but a “step-by-step” movement, which is related to the capillary barrier (Figures 3, 5 and 6). If a capillary barrier is formed in the flow direction and prevents further movement of the oil–water interface, the oil and water will move laterally. If a capillary barrier forms laterally and prevents further sweep, the water drive front will advance along the main flow direction. Then, the process is repeated, that is, the “step-by-step” oil and water movement. Therefore, the capillary barrier controls the movement model of oil and water and then affects their spatial sweep range.

Figure 5 shows the oil and water “step-by-step” movement and the pressure change after the capillary barrier is formed at positions a, b, and c. After water injection, capillary barriers are formed at points a, b, and c, and the pressure of the oil phase is greater than that of the water phase at the corresponding points (Figures 5(a), 5(d) and 6). If water injection is continuous, the capillary pressure breaks at point c and forms a capillary barrier at point d (Figure 5(b)). With the break of c, the pressure state changes, and the downstream pressure of point a decreases, while the downstream pressures of points b and d increase (Figures 5(e) and 6). The downstream pressure rise will impede fluid movement in the parallel channel. After the pressure barrier forms at point d, the oil–water interface breaks again from point a and forms a capillary barrier at point e (Figure 5(c)). The downstream pressure of the corresponding point a increases again, while the downstream pressures of b and d decrease (Figures 5(f) and 6).

5.3. Analysis of Remaining Oil Production after Increasing Liquid Injection Rate

To investigate how different injection rates affect the degree of the remaining oil production, we simulated the distribution of the remaining oil under five injection rates. First, the core was filled with oil, and then, water was injected into the core from inlet A at rate of 0.005 m/s, with B as the extraction outlet. The oil–water viscosity ratio was set to 5 : 1, while the wetting angle was set to 45°. Other flow parameters are listed in Tables 1 and 2. The water flooding process continues until no oil is produced at outlet B. At this time, the oil and water distribution is shown in Figure 7(a), and the process is shown in Figure 3. Then, water was injected from inlet A at a rate of 0.010 m/s, 0.050 m/s, 0.100 m/s, and 0.500 m/s, respectively. The inlet and outlet positions and the fluid properties remain unchanged when the injection rate increases.

Figure 7 shows the final oil–water distribution driven by water at different injection rates. Figure 7(a) shows the initial state of oil and water distribution with an injection rate (0.005 m/s). Oil and water migrate along Channel 1, and the remaining oil is distributed on both sides of the main channel. The oil–water interface at af is at the junction of the pore and throat channel, which forms the capillary barrier. The channel pressure is balanced with the capillary resistance at these interfaces so that the oil–water interface at the upper and lower sides of the main channel is stalled. When water flooding was carried out at a low injection rate (0.010 m/s), the remaining oil was spatially distributed in the lower area of the main channel (Channel 1), and point e was broke, causing a part of the remaining oil to move down along Channel 2, while the upper area did not change (Figure 7(b)). When the injection rate is 0.050 m/s, the pressure increases, exceeding the capillary force at interfaces a and g, so that the capillary barrier at interfaces a and g breaks, and the remaining oil moves to both sides of the main channel along the Channel 3 and Channel 4 (Figure 7(c)). This allows the remaining oil to be produced. Additionally, the existence of an oil block at c blocks the main Channel 1. The pressure on the main channel is balanced with the capillary force at the interface on both sides of the oil block, which is in a stagnant state and drives out the remaining oil at the lower side from the right outlet along Channel 4. When the injection rate increases to 0.100 m/s, the barrier at b is breached, and the remaining oil is produced and moves along Channel 5 (Figure 7(d)). When the injection rate increases to 0.500 m/s, the capillary pressures at c and d on the upper side break, and the capillary valve is formed again at b. The remaining oil on the upper side is produced and moves along Channel 7 and Channel 8, and the barrier at f on the lower side is broken and Channel 6 is formed (Figure 7(e)).

Overall, as the injection rate gradually increases, the oil on both sides of the main channel will be produced successively, and the remaining oil on the edge is becoming increasingly difficult to be produced. The final recovery rate shows that an increase in the injection rate from 0.005 m/s to 0.070 m/s causes the rapid production of the remaining oil and increases the recovery rate to 66% with a linear change (Figure 7(f)). Subsequently, the recovery rate slows down and the remaining oil was increasingly difficult to produce. The recovery rate is only 77% when the extraction rate is 0.500 m/s.

5.4. Mechanical Mechanism of Remaining Oil Production

To quantitatively describe the conditions of the remaining oil production, we analyze the distribution characteristics of oil and water under pore conditions. Figure 8 shows the final oil–water distribution characteristics of water flooding when the injection rate is 0.010 m/s. The main water drive channel (Channel 1) is formed between A and B, and part of the remaining oil moves down along Channel 2 (A ⟶ C), while the remaining oil in the upper area is not produced.

The critical condition of the remaining oil production is related to the pressure drop between A and B, and the capillary barrier at a and b. Water in Channel 1 is moving, and there is a pressure drop p, which can be expressed by formula (27):
(27)p=8μwuLR12,
where μw is water viscosity, u is the injection rate, L is the length of Channel 1, and R1 is the radius of Channel 1.
A can be connected to B through Channel 2 and Channel 3, but there is no movement of oil and water in these two channels because there are capillary barriers pa,cap,pb,cap at points a and b, respectively, which are expressed by the following formula:
(28)pa,cap=σcosθ+βR2,(29)pb,cap=σcosθR3,
where σ is the interfacial tension coefficient, R2 is the radius of Channel 2, R3 is the radius of Channel 3, θ is the wetting angle, and β is the opening angle.
If remaining oil from paths A ⟶ C ⟶ B is to be produced, the pressure drop between A and B should be greater than the sum of the capillary barriers at a and b:
(30)p>pa,cap+pb,cap.
By integrating equations (27)–(30), we can obtain
(31)u>σcosθ8R3μwLR12σcosθ+β8R2μwLR12.

It can be concluded that if the remaining oil is to be produced, it is necessary to increase the viscosity of water, reduce the interfacial tension, and increase the injection rate. Simultaneously, it is also closely related to the pore radius and the wettability of the pore wall. Changing the wettability may not have an obvious effect on the production of the remaining oil.

In this study, a digital core model is established based on CT scan data, and the numerical method used in the simulation process of pore scale water flooding is analyzed. Then, the model of oil–water movement using a constant injection rate is simulated, and the influence of capillary force on oil–water movement during water flooding is discussed.

  • (1)

    The oil–water interface stays at the position between the throat channel and the pore area, the oil–water interface reverses, and the capillary force presents resistance, forming the capillary barrier or capillary valve

  • (2)

    The oil and water movement during water flooding is not stable and continuous, but a “step-by-step” movement model, which is closely related to the capillary barrier

  • (3)

    In the low-permeability sandstone reservoir of the Ansai Oilfield in the Ordos Basin, as the liquid injection rate gradually increases, the oil on both sides of the main channel is successively produced. After the recovery rate reaches 66%, the remaining oil on the edge becomes increasingly difficult to be used

  • (4)

    By increasing the viscosity of water, reducing the interfacial tension, and increasing the injection rate to a certain value, the degree of the remaining oil production can be improved. The difficulty of the remaining oil production is closely related to the viscosity of the injection fluid, interfacial tension, influence of injection velocity, pore radius, and wettability of the pore wall

The data that support the findings of this study are available from the corresponding author (Junwei Su) upon reasonable request.

The authors declare that there is no conflict of interest regarding the publication of this paper.

This work is funded by the National Natural Science Foundation of China (No. 12072256) and the Key Research and Development Program of Shaanxi (Program No. 2021GXLH-Z-071).

Exclusive Licensee GeoScienceWorld. Distributed under a Creative Commons Attribution License (CC BY 4.0).