Abstract
Oil recovery, including conventional and viscous oil, can be improved significantly by flooding with polymer solutions. This chemical flooding method can increase oil production, and it can improve the macrodisplacement efficiency and microsweep efficiencies. In this study, we establish physical models that include the dead-end and complex models based on the pore-network pattern etched into glass, using the snappyHexMesh solver in OpenFOAM. These models capture the complexity and topology of porous media geometry. We establish a mathematical model for transient flows of viscoelastic polymers using computational fluid dynamics simulations, and we study the distributions of pressure and velocity for different elasticity scenarios and different flooding process. The results demonstrate that the pressure difference increases as the relaxation time decreases, before the flow reaches its steady state. For a steady flow, elasticity can give rise to an additional pressure difference, which increases with increasing elasticity. Thus, the characteristics of pressure difference vary before and after the flow becomes steady; this phenomenon is very important. Velocity contours become more widely spaced with elasticity increase. This suggests that elasticity of the polymer solutions contributes to the microsweep efficiency. The results of the study provide the necessary theoretical foundation for laboratory experiments and development of methods for polymer flooding and can be helpful for the design and selection of polymers for polymer flooding.
1. Introduction
With continuing development of energy infrastructure, the supply and consumption of energy have come to play an important role in the industrialization and urbanization of world countries. Oil plays the most important role, because as of today it cannot be fully replaced by other energy sources. To ensure the sufficient production of oil, many measures have been adopted, including exploration for discovering new oil reserves, creating national oil strategic reserve systems, and avoiding import risk based on the market fluctuations of oil futures. Enhanced oil recovery from developed reservoirs is an important strategy to face the above-described challenges [1, 2]. Among various chemical flooding techniques for enhanced oil recovery, polymer flooding is especially promising. Of late, polymer flooding research has matured enough for commercial applications. Taking the Daqing oilfield as an example, 54 industrial blocks have been displaced by polymer flooding. Compared with water flooding, polymer flooding increased oil recovery by 13%. Polymer flooding-based recovery of second-type oil layers enhanced oil recovery by 10% compared with water flooding [3]. In addition, polymer flooding has been used in heavy-oil reservoirs [4], especially for offshore heavy oil. Polymer flooding is suitable for operation on offshore platforms because it does not require complex and additional surface facilities given a relatively limited platform space. Bohai Bay is the most important offshore oilfield in China, which produces heavy oil, and polymer flooding has significantly increased the oil recovery there. In three years, in the pilot area of four injection wells and six production wells, incremental oil production increased to 7.6 × 104 m3 [5, 6]. The main factor contributing to the success of this technique is that the viscosity of water can be increased by polymers, resulting in a better macroscale sweep efficiency [7, 8]. In terms of microscopic mechanisms, the polymer elasticity can contribute to pulling out residual oil droplets. Thus, on the microscale, polymer elasticity in porous media is important for increasing the oil recovery [9–12].
To investigate the pore-scale flow mechanisms, researchers have employed microscale visual experiments and theoretical computational fluid dynamics (CFD) methods and studied the effects of viscosity and elasticity on flow characteristics. A two-dimensional (2D) etched-silicon microscale model having the geometrical and topological characteristics of stone has been prepared. Experimental results for heavy oil displaced by water and polymer show that, at zero and low polymer concentrations, relatively long and wide fingers of injectant develop and slowly recover. As polymer concentration increases, many more relatively fine fingers are formed [13]. Wang and Xia compared viscoelastic polymers and pure viscous glycerin solutions having the same viscosity to displace oil, by conducting experiments with glass etched system [14–16]. By comparing images after flooding by two fluids, polymer elasticity was shown to be responsible for pulling out residual oil droplets. The oil recovery associated with polymer solutions was enhanced by 7%–16% compared with glycerin having the same viscosity. Vermolen et al. [17] designed and performed core flooding experiments using glycerol and polyacrylamide, in which they validated that polymer elasticity can enhance the displacement efficiency. The results of these experiments showed that, for low viscosity crude oil, extra oil was recovered when using polymers. This confirms that high-elasticity polymers can reduce the amount of residual oil compared with water flooding.
The polymer elasticity can enhance the displacement effect, as has been validated by the above experimental studies. Theoretical investigations about the flow characteristics of polymer solutions in the context of oil industry are lacking. Some researchers have studied steady-flow polymer solutions for 4 : 1 contraction and dead-end pores using theoretical method. The finite difference method was used to study the stream function, velocity, and pressure contours, for different viscoelastic polymer flows in microscale pores. The results of these studies show that, with increasing polymer elasticity, the flow area increases while the residual oil area decreases. The velocity and pressure also increase [18–20]. Static numerical simulations (water-based fluid flowing around an immobile oil droplet) showed that additional forces are present in the case of elastic fluids. Dynamic simulations revealed striking differences in the behavior of elastic and inelastic fluids; oil droplets underwent deformation over time when inelastic fluids were injected, but the droplet did not deform in the presence of elastic fluids [21]. Owing to the nonlinearity associated with elasticity, the numerical simulation of viscoelastic flows remains challenging. The stabilized method of finite elements was used in some studies [22]. This method yielded high elasticity (Weissenberg number > 2), and the results suggest that viscoelasticity may play a significant role in determining the behavior of polymer flows. Earlier studies of viscoelastic behavior focused on steady flows. In 2011, Jafari [23] developed a spectral element method to obtain an accurate numerical solution for viscoelastic fluid flows. The objective was to introduce a new algorithm for overcoming the drawbacks associated with the direct reformulation of the classical constitutive equation as a logarithmic one. The research effort culminated in a comprehensive study on temporal growth of spurious modes. Prasuna et al. [24] solved the problem of viscoelastic unsteady flow through two impermeable parallel plates using Laplace transformation technique. The equations were solved for stage-steady and unsteady motion, respectively. The researchers obtained expressions for the flow rate and shear stress on the walls. In most theoretical studies, simplified models such as the dead-end model, contraction structures, and expansion-contraction structures were adopted [25]. However, the structure of a real porous medium is very complex and so is likely to affect the flow characteristics. Hence, it is very important to consider the complexity of pores in studies of polymer elasticity. In this case, the finite volume method [26–29] is likely to be more stable. Also, the existing studies mainly focus on viscoelastic steady flows, and although some research has been done on unsteady flow, those studies address the algorithms rather than the flow characteristics of viscoelastic fluids in porous media. Transient flow characteristics (before a flow reaches a steady state) are important, leading to some potentially interesting phenomena that are ignored in experiments. In this study, we investigated the complex pore geometries based on etched-in glasses model. This 2D model can capture the geometrical and topological characteristics of pores in a realistic manner. Furthermore, we studied the microscopic characteristics of transient flows using the viscoelastic flow solver in OpenFOAM. These results provide theoretical support and guidance to studies that aim to enhance oil recovery using polymer flooding.
2. Physical Model
Tow physical models that represent the pores have been established in this study. The first one is the dead-end model, which captures a simple pore. The second one that models more complex pore structures is represented using the etched-in glass model. The dead-end model is schematically shown in Figure 1. The residual oil in the dead end often cannot be displaced by water. The real etched-in glass model is 40 × 40 mm in size and was used in microscale visual experiments with pure viscous glycerin and elastic polymer flooding. Its porosity is 35% and its permeability was 2200 × 10−3 μm. Figure 2(a) shows the saturated oil, imaged using a high-resolution camera. Figure 2(b) shows a region from Figure 1 as a digitized model to demonstrate the computational speed and convergence, and area under consideration is 1 × 1 mm. The diameter of the smallest pore was ~0.02 mm. In conclusion, this model provides a sufficiently detailed picture of a porous medium to allow studies of the polymer flow behavior.
(a) Original saturated oil micromodel
(b) The digitized model
(c) The stl file exported from Solidworks
(d) Meshes exported from snappyHexMesh
First, the image was digitized using Solidworks to obtain a black and white presentation, and, then, the black pixels corresponding to the pore boundary were extracted (Figure 1(c)). The skeleton of the pores must be imported using snappyHexMesh, a solver in OpenFOAM, and is shown in Figure 1(c). The nonorthogonal grid was obtained using the snappyHexMesh solver in OpenFOAM, and it is shown in Figure 1(d). The number of mesh elements in the model was ~399220. The minimum radius was 0.02 mm. The mesh file was designed so it could be accepted by any solver in OpenFOAM.
3. Mathematical Model
The following simulation addressed a viscoelastic transient flow in an incompressible and isothermal 2D system. The governing equations included the continuity equation, the momentum equation, and the constitutive equation.
The continuity equation isThe momentum equation iswhere is the velocity vector, in m/s, is the velocity in the direction, and is the velocity in the direction. Other parameters were as follows: flow time , in ; polymer density , in kg/m3; stress tension , in Pa, which had three components, with and corresponding to the first and second normal stresses, respectively, and corresponding to the shear stress; pressure , in Pa.
The polymer molecules with large relaxation time do not have sufficient time to stretch, recoil, and adjust themselves to the elongation and contraction of flow through pores and throats in the reservoir rock. Therefore, the resultant elastic strain produces high apparent viscosity. Viscoelastic polymer exhibits nonlinear rheological responses, which are captured using the constitutive equation. Strong experimental evidence indicates that the first normal stress difference is dominant and is larger than the second normal stress difference. The Oldroyd B equation can describe this property of polymer solutions [30], as follows:where ; is the polymer viscosity, in Pas; is the deformation rate tensor, in s−1; is the relaxation time, in s, which captures the elasticity of the polymer solutions. It has been suggested that the inverse of the frequency at which and intersect corresponds to the characteristic relaxation time of the polymer solution.
To ensure the convergence of numerical solutions, we used the DEVSS approach for solving the system of equations. The stress tensor was divided into two parts, including the solvent contribution part (viscous effect) and polymer contribution part (elastic effects); together, these allow representing clearly the nonlinear viscoelasticity phenomena [31]. The corresponding equation was where is the solvent contribution to stress, in Pa, and is the polymeric contribution, in Pa: The parameter could be calculated from the Oldroyd B constitutive differential equation, and the momentum equation was rewritten as follows:where is the solvent viscosity, in Pas, and is polymeric viscosity, in Pas.
Equations (1), (2), and (6) constitute the governing equations for solving the problem. The finite volume method was used to solve the above equations. The steps used for solving the system are given below.
(1) For a given initial velocity field, pressure field, and stress tensor field, momentum equation is solved, obtaining a new velocity field, where the pressure gradient and stress divergence are calculated by explicit method.
(2) Using the new velocity field, according to the RhieChow interpolation idea, construct the Poisson equation for pressure, solve for the pressure field, and correct the velocity to satisfy the continuity equation. For the transient problem we use the PISO method, and for the steady-state problem we use the SIMPLE method.
(3) Using the velocity field that satisfies the continuity equation, use the constitutive equation of the viscoelastic fluid to calculate the stress tensor.
(4) Repeat steps (1)–(3).
The equations were solved using the OpenFOAM’s viscoelasticFluidFoam solver.
In these models, the density of the polymer flow was 803 kg/m3 and viscosity was 8 mPas. The relaxation time values were 0.06 s, 0.1 s, and 0.15 s.
4. Flow Characteristics in the Dead-End Model
In these flow cases, the flow velocity has components in both the and directions. The flow velocity was defined:
Figure 3 shows the velocity contours at different times. It can be seen that the maximal velocity decreases as time increases. Because the process is unsteady, the polymer flows into the main flow pore from the inlet. This is manifested as an oscillation owing to the elasticity of the polymer solutions. The flow becomes steady at = 1 s. Thus, the flow velocity becomes constant at = 1 s. The maximal velocity does not change with time.
(a) 0.1 s
(b) 0.4 s
(c) 0.7 s
(d) 1 s
The pressure differences between the inlet and the outlet, for different relaxation times, are shown in Figure 4. For the same relaxation time, the pressure difference increases with the time, and it becomes constant for a steady flow. At the same time, for transient flows, the pressure difference increases with relaxation time decrease. For steady flows, the difference increases with relaxation time increase. The kinetic energy is transformed into elastic potential energy; thus, the pressure difference is larger for smaller relaxation times, before the flow becomes steady. Elastic potential energy is transformed into kinetic energy when the flow is steady; thus, in this regime the difference is larger for larger relaxation times.
5. Flow Characteristics in the Complex Model
Based on the dead-end model, we simulated the flow of viscoelastic polymers in a complex microscale model. Figures 5, 6, and 7 show the results of these simulations, for different solutions of elastic polymers.
(a) 0.1 s
(b) 0.5 s
(c) 2 s
(a) 0.1 s
(b) 0.5 s
(c) 2 s
(a) 0.1 s
(b) 0.5 s
(c) 2 s
Figures 5, 6, and 7 demonstrate that, in the porous media, the inlet-outlet pressure difference increases with time, and, for the same relaxation time, the pressure is much higher in pores compared with throats. The figures show that the pressure difference increases as the relaxation time decreases for the same injection time, when the flow is transient. This occurs because the kinetic energy is transformed into elastic energy, the elasticity increasing and the kinetic energy decreasing, and then the pressure difference reducing.
Figure 8 shows the pressure difference between the inlet and outlet of the polymer solution at different times. For the same injection time, before the flow becomes steady, the figure shows that the pressure difference increases with the elasticity decrease. For steady flows, the pressure difference increases as the relaxation time increases; thus, the elasticity of the polymer solution can produce an additional pressure drop. The greater the elasticity, the stronger the pressure drop. This result is consistent with the flow in the dead-end model.
Figure 9 shows the velocity contour with the flow time at the same relaxation time ( = 1.5 s). For flow times under 1 s, the flow velocity increases; however, for a steady flow, the velocity decreases and reaches a constant. Figure 10 shows the velocity contour with the different relaxation time for steady flows. As shown in Figure 10, the flow velocity increases with the elasticity increasing. Using the same maximal velocity, the obtained contour is shown in Figure 11.
(a) s
(b) s
(c) s
(a) s
(b) s
(c) s
(a) s
(b) s
(c) 5 s
Figure 11 also shows that when the maximal and minimal velocities are the same, for different relaxation times, the area swept by the same contours become larger as the relaxation time increases in three local areas. This phenomenon demonstrates that the elasticity of a polymer solution can increase the microsweep area and sweep efficiency in porous media.
6. Conclusions
In this paper, a mathematical model to describe the flow viscoelastic polymer solution in a porous medium was established and solved using the OpenFOAM solver. The microscale model, based on the etched-glass, can capture the complexity and topology of porous media.
For both models, the simulations of transient viscoelastic polymer flows revealed that the flow pressure difference increase as the elasticity decrease before the flow became steady. For a steady flow, the pressure difference increased as the elasticity increased. The flow velocity increased with flow time, and, in the flow steady state, the flow velocity was constant.
The simulations of velocity contours in the micromodel showed that the microsweep area increased as the relaxation time increased. The microsweep efficiency increased as the elasticity increased.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This work presented in this paper was financially supported by the National Natural Science Funds for Young Scholars of China (Grant no. 51604079), the Natural Science Foundation of Heilongjiang Province (Grant no. E2016015), and the University Nursing Program for Young Scholars with Creative Talents in Heilongjiang Province (Grant no. UNPYSCT-2016126), which are gratefully acknowledged.