Numerical Modeling of Nanotechnology-Boosted Chemical Enhanced Oil Recovery Methods

Abstract

Since it was theorized more than 50 years ago, nanotechnology has become the perfect boost for existing old technologies. The unique properties exhibited by materials at these scales have a potential to improve the performance of mature oil fields along with enhanced oil recovery (EOR) processes. Regarding polymer flooding, the influence of the (macro) molecules’ architecture on the fluid properties has been lately stressed. This chapter presents the numerical simulation of the combination of both agents in a single, combined recovery process. The presence of the nanoparticles affects the rheological behavior and the rock’s wettability, increasing the organic phase mobility. Undesirable effects such as (nano) particle aggregation and sedimentation are also considered. The polymer’s architecture has a major influence on the recovery process, improving the rheological and viscoelastic properties. On the other hand, although nanoparticles improve the viscosity as well, its main mechanism is their adsorption onto the rock and wettability modification. This chapter shows the importance of a good polymer characterization for EOR, the potential of nanoparticles acting as a boost of traditional EOR processes, and the vital role CFD techniques play to assess the potential of these agents and the optimization of the recovery strategies.

1. Introduction

The so-called era of discovery and exploitation of the denominated “easy oil” is since some years reaching to an end [1, 2, 3, 4, 5]. The exploitation of a conventional oil field can be mainly divided in three stages, which depend on the physical mechanisms acting during the oil recovery [6, 7, 8]. The first step consists in taking advantage of the natural-driven mechanisms present in the oil field, without the injection of fluids or specific agents. This stage, known as primary recovery, finished when the pressure in the reservoir or the amount of oil produced is no longer enough to render a profitable production. Subsequently, water or gas starts being injected with a dual goal: repressurize the rock formation and sweep the remaining oil toward the producer wells (Figure 1). Both primary and secondary recoveries account for approximately a 50% of the original oil in place (OOIP). During the last 50 years, oil companies began applying more advanced processes after secondary recovery, which are known as enhanced oil recovery (EOR) or tertiary processes, which involve the injection of different fluids in order to modify the physical properties of the different fluids and/or rock formation. This renders an increase in the oil production and lifetime of conventional oil sources, which is necessary while newer and greener energy sources are developed and optimized (Figure 2).

Figure 1.

Anticlinal type petroleum trap [9].

Figure 2.

Schematic representation as a function of time of the different oil recovery stages and their respective productivities [12].

Among EOR techniques, chemical EOR comprises the addition of certain agents (e.g., polymers, surfactants, alkali) and presents a great potential in low- and medium-viscosity mature oil fields to increase their productivity. Recently, polymers were also applied in high-viscosity reservoirs showing promising results when specific geological conditions are met. Nevertheless, there are certain problems associated with these products which hinder their efficiency and limit their applicability. Thus, during the last 20 years, researchers have begun making use of the nanotechnology in order to boost the efficiency or EOR processes, based on the novel properties of materials exhibited at these scales. The combination of chemical EOR agents with different nanomaterials has shown increased recovery efficiencies in rock formations which had otherwise reached their operational limit [10, 11].

However, the potential advantages of the application of nanomaterials in oil recovery should also be carefully considered from an environmental point of view. Some concerns have been raised since the same features that make nanotechnology so attractive to oil recovery processes might also have a negative impact on the environment and human health. These include the potential long-term side effects associated with medical applications as well as with the biodegradability of nanomaterials being used [13, 14, 15, 16, 17, 18, 19]. Even though during the last years many nanomaterials have been inserted into the market, the amount of information over their impact on the environment is minimal. For instance, there is almost no information about the associated risks in the manufacturing, usage, and final disposal of nanomaterials [20]. Focusing on EOR techniques, it has been shown that a percentage of the nanoparticles will remain underground and deposited in the rock formation, remaining for many years, and this could cause the contamination of groundwater sources. Thus, one of the desirable properties of these nanoparticles should be high durability and recyclability in a cost-effective process to decrease their impact on the environment [20].

An important point to mention is that the field test of these new techniques involves a significant use of resources and time in order to assess their efficiency. Therefore, for the last 40 years, scientists started using computers to perform this assessment and save considerable time and physical resources, which is known as reservoir simulation. This is a branch of computational fluid dynamics (CFD) which involves solving the balance equations present in porous media, which renders a number of coupled, highly nonlinear systems of equations dealing with temporal and spatial variations of pressure and mass concentrations. Different numerical and physical techniques have been applied during the years in order to simulate different recovery processes as well as to increase the numerical accuracy of the simulation.

The objective of this chapter is to present the potential of nanomaterials as a boost of traditional EOR techniques, focusing on the development of numerical models for reservoir simulation, especially in the combination of chemical EOR agents with nanoparticles, studying their advantages and synergy in order to increase the productivity of conventional oil sources. Reservoir simulation consists broadly speaking of three parts: geological, fluid, and well models which describe the main parts of the extraction system [21, 22, 23, 24]. The accurate mathematical representation of the whole system is still a topic in which further research is needed. The risks associated with the uncertainties in the numerical model might lead to the failure of exploration and production (E&P) projects. This rendered a multiphase, multicomponent model, considering all the effects that chemical agents and nanoparticles provoke both in the fluid and rock formation. Regarding the nanoparticles, this includes the aggregation, retention, rheology, and changes in permeability and porosity. The chemical agent studied in this chapter is a polymer, which is modeled considering the influence of the (macro)molecules’ architecture on the fluid properties. The salt present in the reservoir is also considered in the water phase. This model rendered a novel simulator, combining the benefits of nanotechnology with chemical EOR processes.

2. Model description

2.1 Physical model

The study of this combined model of polymer-boosted EOR flooding is presented in a 2D domain, based on a well configuration used in the oil industry. Indeed, the five-spot scheme consists in a square domain, in which the injection well is placed in the center and four production wells in each of the vertices (Figure 3a). From this, a mathematic simplification consists in dividing this domain using its symmetry, which is known as quarter five-spot (Figure 3b).

Figure 3.

Five-spot scheme (a) indicating the well’s location and the possible presence of faults and its simplification to the quarter five-spot used during this chapter (b) [12, 32].

The physical model ($\mathrm{\Omega }$) represents an oil field of known geometric and physical properties, i.e., absolute permeabilities ($K$) and porosity ($\varphi$) and a rock compressibility ($c_f$). Furthermore, the fluid flow is considered isothermal and incompressible. Darcy’s law is valid, and the gravitational forces are negligible compared to the viscous ones [25, 26, 27]. Thus, numerically speaking, the domain is discretized in a number of $n_x\times n_y$ blocks. The grid size is chosen based on the analysis of the representative elementary volume (REV), which is determined by the minimum size in which the rock properties remain approximately the same.

The recovery process consists in a two-phase (aqueous and organic), multicomponent (water, salt, polymer nanoparticles, and petroleum) flow. These components may be also mixtures of a number of pure ones, e.g., petroleum is a mixture of many hydrocarbons, water contains dissolved minerals (other than salt itself), and finally polymer is composed of different molecules of different lengths and architectures [27]. The polymer properties are determined by the average molecular weight, which in this model is assumed to be identical for all the molecules, which renders a polydispersity index ($\mathit{PDI}$) equal to unity. The nanoparticles affect the water phase and the rheology, using a function of their concentration and size. With respect to this, aggregation mechanisms present in the system tend to increment the nanoparticles’ average size (Figure 4) [28, 29, 30, 31].

Figure 4.

Schematic representation as a function of time of the different oil recovery stages and their respective productivities [33].

The mathematical description of the system is represented by a number of strongly nonlinear partial differential equations complemented by a set of algebraic relationships describing the physical properties of fluids and rock, which are aggregation of nanoparticles, degradation of polymer molecules, interfacial tension, residual phase saturations, relative permeabilities, rock wettability, phase viscosities, capillary pressure, adsorption and retention of both polymer and nanoparticles onto the formation, inaccessible pore volume (IAPV), disproportionate permeability reduction (DPR), nanoparticles-polymer interactions, and dispersion.

2.2 Mathematical model

To study the flow of multiphase, multicomponent system in porous media, the mass, momentum, and energy balance equations are applied. Therefore, the equations used to describe the process are Darcy’s equation applied to each phase and the mass conservation valid for each component [33]. Since in most of chemical EOR processes the exchange of energy (i.e., temperature changes) is not significant, the energy balance equation is not considered in the simulation. The compositional approach is chosen because of its versatility to model the different physical properties according to the components’ concentrations. These equations are then applied on a REV of the porous medium. Considering Darcy’s equation for each phase first,

$${\overrightarrow{u}}^j=-\underset{¯}{\underset{¯}{K}}\cdot \frac{k_r^j}{{\mu }^j}\cdot \overrightarrow{\nabla }{p}^j;j=o,a$$E1

$$\frac{\partial \left(\varphi z_i\right)}{\partial t}+\nabla \cdot \sum _j{V}_i^j\cdot {\overrightarrow{u}}^j-\nabla \cdot \sum _j\underset{¯}{\underset{¯}{D_i^j}}\cdot \nabla \cdot V_i^j=-\frac{\partial \left(\varphi {\mathit{Ad}}_i\right)}{\partial t}+q_i;i=p,\mathit{np},w,s,\mathit{pol}$$E2

$$\underset{¯}{\underset{¯}{D_i^j}}={\mathit{dm}}_i^j\cdot \varphi \cdot S^j\cdot {\delta }_{\mathit{ij}}+\parallel {\overrightarrow{u}}^j\parallel \cdot \left[\frac{{\mathit{dl}}^j}{\parallel {\overrightarrow{u}}^j{\parallel }^2}\cdot \left|\begin{array}{cc}{\left(u_x^j\right)}^2& u_x^j\cdot u_y^j\\ u_y^j\cdot u_x^j& {\left(u_y^j\right)}^2\end{array}\right|+{\mathit{dt}}^j\cdot \left|\begin{array}{cc}1-\frac{{\left(u_x^j\right)}^2}{\parallel {\overrightarrow{u}}^j{\parallel }^2}& -\frac{u_x^j\cdot u_y^j}{\parallel {\overrightarrow{u}}^j{\parallel }^2}\\ -\frac{u_y^j\cdot u_x^j}{\parallel {\overrightarrow{u}}^j{\parallel }^2}& 1-\frac{{\left(u_y^j\right)}^2}{\parallel {\overrightarrow{u}}^j{\parallel }^2}\end{array}\right|\right]$$E3

$$\varphi c_r\frac{\partial p^a}{\partial t}+\overrightarrow{\nabla }\cdot \left(\lambda \cdot \nabla p^a\right)=\frac{\partial }{\partial t}\left(\varphi \cdot \sum _i{\mathit{Ad}}_i\right)-\overrightarrow{\nabla }\cdot \left({\lambda }^o\cdot \nabla p_c\right)+q_t$$E4

As in every compositional simulation, the number of equations generated by the conservation laws is not sufficient to mathematically determine the system. Thus, a number of auxiliary relationships are needed in order to find the solution, which for a five-component, two-phase system this renders a total of $N_{\mathit{\text{comp}}}\cdot \left(N_{\mathit{\text{phases}}}-1\right)=5$ equations, determined by system’s phase behavior [24, 34, 35].

2.2.1 Discretization of the partial differential equations

The above system of equations are discretized and then solved in this chapter by a finite difference method, which is one of the most well-known techniques in CFD. The first equation is the aqueous phase pressure (Eq. (4)), which is implicitly discretized using a centered scheme for the pressure terms and a second-order Taylor approximation for the time derivatives. This scheme chosen for the simulator is often used in systems with derived second-order with coefficients that are not constants. Besides the Darcy momentum equation, the discretization of the total and aqueous velocities is also done using a centered difference scheme. Therefore, Eqs. (1) and (4) are discretized as

$$\begin{array}{c}{c}_r{\left(\varphi +\frac{\mathrm{\Delta }t}{2}\frac{\partial \varphi }{\partial t}\right)}_{m,n}^{<n+1>,\left[k\right]}{\left(\frac{p_{m,n}^{a,<n+1>}-p_{m,n}^{a,<n>}}{\mathrm{\Delta }t}\right)}^{\left[k+1\right]}+\frac{{\lambda }_{x,m+12,n}^{<n+1>,\left[k\right]}}{\mathrm{\Delta }{x}^2}\cdot {\left(p_{m+1,n}^a-p_{m,n}^a\right)}^{<n+1>,\left[k+1\right]}-\dots \\ -\frac{{\lambda }_{x,m-12,n}^{<n+1>,\left[k\right]}}{\mathrm{\Delta }{x}^2}\cdot {\left(p_{m,n}^a-p_{m-1,n}^a\right)}^{<n+1>,\left[k+1\right]}+\frac{{\lambda }_{y,m+12,n}^{<n+1>,\left[k\right]}}{\mathrm{\Delta }{y}^2}\cdot {\left(p_{m,n+1}^a-p_{m,n}^a\right)}^{<n+1>,\left[k+1\right]}-\dots \\ -\frac{{\lambda }_{y,m-12,n}^{<n+1>,\left[k\right]}}{\mathrm{\Delta }{y}^2}\cdot {\left(p_{m,n}^a-p_{m,n-1}^a\right)}^{<n+1>,\left[k+1\right]}={\left(\varphi +\frac{\mathrm{\Delta }t}{2}\frac{\partial \varphi }{\partial t}\right)}_{m,n}^{<n+1>,\left[k\right]}{\left(\frac{{\mathit{Ad}}_{m,n}^{<n+1>}-{\mathit{Ad}}_{m,n}^{<n>}}{\mathrm{\Delta }t}\right)}^{\left[k+1\right]}+\dots \\ +{\left(\mathit{Ad}+\frac{\mathrm{\Delta }t}{2}\frac{\partial \mathit{Ad}}{\partial t}\right)}_{m,n}^{<n+1>,\left[k\right]}{\left(\frac{{\varphi }_{m,n}^{<n+1>}-{\varphi }_{m,n}^{<n>}}{\mathrm{\Delta }t}\right)}^{\left[k+1\right]}+\frac{{\lambda }_{x,m+12,n}^{o,<n+1>,\left[k\right]}}{\mathrm{\Delta }{x}^2}\cdot {\left(p_{m+1,n}^c-p_{m,n}^c\right)}^{<n+1>,\left[k+1\right]}-\dots \\ -\frac{{\lambda }_{x,m-12,n}^{o,<n+1>,\left[k\right]}}{\mathrm{\Delta }{x}^2}\cdot {\left(p_{m,n}^c-p_{m-1,n}^c\right)}^{<n+1>,\left[k+1\right]}+\frac{{\lambda }_{y,m+12,n}^{o,<n+1>,\left[k\right]}}{\mathrm{\Delta }{y}^2}\cdot {\left(p_{m,n+1}^c-p_{m,n}^c\right)}^{<n+1>,\left[k+1\right]}-\dots \\ -\frac{{\lambda }_{y,m-12,n}^{o,<n+1>,\left[k\right]}}{\mathrm{\Delta }{y}^2}\cdot {\left(p_{m,n}^c-p_{m,n-1}^c\right)}^{<n+1>,\left[k+1\right]}+q_{m,n}^{t,<n+1>,\left[k\right]}\end{array}$$E5

$$\begin{array}{c}{\overrightarrow{u}}_{m,n}^{<n+1>,\left[k+1\right]}=\left[-\frac{{\lambda }_{x,m,n}^{\left[k\right]}}{2\cdot \mathrm{\Delta }x}\cdot {\left(p_{m+1,n}^a-p_{m-1,n}^a\right)}^{\left[k+1\right]}-\dots \right.\\ {\left.-\frac{{\lambda }_{x,m,n}^{o,\left[k\right]}}{2\cdot \mathrm{\Delta }x}\cdot {\left(p_{c,m+1,n}^a-p_{c,m-1,n}^a\right)}^{\left[k+1\right]}\right]}^{<n+1>}\cdot \widehat{i}+\left[-\frac{{\lambda }_{y,m,n}^{\left[k\right]}}{2\cdot \mathrm{\Delta }y}\cdot {\left(p_{m,n+1}^a-p_{m,n-1}^a\right)}^{\left[k+1\right]}-\dots \right.\\ {\left.-\frac{{\lambda }_{y,m,n}^{o,\left[k\right]}}{2\cdot \mathrm{\Delta }y}\cdot {\left(p_{c,m,n+1}^a-p_{c,m,n-1}^a\right)}^{\left[k+1\right]}\right]}^{<n+1>}\cdot \widehat{j}\\ \end{array}$$E6

$$\begin{array}{c}{\overrightarrow{u}}_{m,n}^{a,<n+1>,\left[k+1\right]}={\left[-\frac{{\lambda }_{x,m,n}^{a,\left[k\right]}}{2\cdot \mathrm{\Delta }x}\cdot {\left(p_{m+1,n}^a-p_{m-1,n}^a\right)}^{\left[k+1\right]}\right]}^{<n+1>}\cdot \widehat{i}+\dots \\ +{\left[-\frac{{\lambda }_{y,m,n}^{a,\left[k\right]}}{2\cdot \mathrm{\Delta }y}\cdot {\left(p_{m,n+1}^a-p_{m,n-1}^a\right)}^{\left[k+1\right]}\right]}^{<n+1>}\cdot \widehat{j}\\ \end{array}$$E7

where $m,n$ represent the cells of the numerical domain $\left(x,y\right)=\left(m\cdot \mathrm{\Delta }x,n\cdot \mathrm{\Delta }y\right)$, respectively, $<n>$ is the temporal step $\left(\text{time}=<n>\cdot \mathrm{\Delta }t\right)$, and $\left[k\right],\forall k\in {\mathbb{N}}^{+},$ is the iteration number within each time-step. Finally, mass conservation equation is discretized using a second-order approach. Eq. (2) is the typical advection-diffusion PDE used to describe the mass transport in porous media. The advective terms are of hyperbolic nature, and first-order numerical schemes cause an artificial diffusion in the solution. There are different approaches in order to overcome this, and one of them is the use of higher-order schemes. The proposed simulator in this chapter uses a full second-order explicit discretization scheme in time and space, based on total variation diminishing (TVD) and flux-limiting techniques. This increases the numerical accuracy of the simulator as well as decreases the influence of numerical diffusion and dispersion. Diffusive terms are discretized using a centered second-order scheme. The second order in time is achieved using a Taylor expansion up to the second order [12]. The flux-limiting techniques require to establish a functional relationship between the gradient of the volumetric concentration and the limiting function $\psi$. Different second-order methods have been studied, with the functions utilized in the development of this simulator presented in Table 1 (together with the standard upwind method). These functions depend on the ratio of the concentrations’ consecutive gradients in the numerical grid $\left(r_{x,i}=\left(V_{i,m,n}^{j,\left[k\right]}-V_{i,m-1,n}^{j,\left[k\right]}\right)/\left(V_{i,m+1,n}^{j,\left[k\right]}-V_{i,m,n}^{j,\left[k\right]}\right)\right)$.

Flux limiter	Function
Upwind	0
Superbee	$max\left[0,min\left(2r,1\right),min\left(r,2\right)\right]$
Minmod	$max\left[0,min\left(r,1\right)\right]$
MUSCL	$max\left[0,min\left(2r,\frac{1+r}{2},2\right)\right]$

Table 1.

General parameters used for the simulations.

Thus, the discretized mass conservation equation is

$$\begin{array}{c}\frac{{\mathcal{C}}_1}{\mathrm{\Delta }t}{z}_i^{<n+1>}={\mathcal{C}}_2{z}_i^{<n>}+\frac{{\mathcal{C}}_3}{\mathrm{\Delta }x}\cdot {\sum }_j{F}_{\mathit{LIM},x}^{j,<n+1>,\left[k+1\right]}{\left(u_{x,m,n}^{j,\left[k+1\right]}\cdot V_{i,m,n}^{j,\left[k\right]}-u_{x,m-1,n}^{j,\left[k+1\right]}\cdot V_{i,m-1,n}^{j,\left[k\right]}\right)}^{<n+1>}+\dots \\ +\frac{{\mathcal{C}}_3}{\mathrm{\Delta }y}\cdot {\sum }_j{F}_{\mathit{LIM},y}^{j,<n+1>,\left[k+1\right]}{\left(u_{y,m,n}^{j,\left[k+1\right]}\cdot V_{i,m,n}^{j,\left[k\right]}-u_{y,m,n-1}^{j,\left[k+1\right]}\cdot V_{i,m,n-1}^{j,\left[k\right]}\right)}^{<n+1>}+\dots \\ +\frac{1}{\mathrm{\Delta }{x}^2}\cdot {\sum }_j{\left[{\left(S^j\varphi {\mathit{dm}}_i^j\right)}_{m+1/2,n}\cdot \left(V_{i,m+1,n}^j-V_{i,m,n}^j\right)-{\left(S^j\varphi {\mathit{dm}}_i^j\right)}_{m-1/2,n}\cdot \left(V_{i,m,n}^j-V_{i,m-1,n}^j\right)\right]}^{<n+1>,\left[k\right]}+\dots \\ +\frac{1}{\mathrm{\Delta }{y}^2}\cdot {\sum }_j{\left[{\left(S^j\varphi {\mathit{dm}}_i^j\right)}_{m,n+1/2}\cdot \left(V_{i,m,n+1}^j-V_{i,m,n}^j\right)-{\left(S^j\varphi {\mathit{dm}}_i^j\right)}_{m,n-1/2}\cdot \left(V_{i,m,n}^j-V_{i,m,n-1}^j\right)\right]}^{<n+1>,\left[k\right]}-\dots \\ -\frac{1}{\mathrm{\Delta }t}{\left(\varphi +\mathrm{\Delta }t\frac{\partial \varphi }{\partial t}\right)}_{m,n}^{<n+1>,\left[k+1\right]}\cdot {\left({\mathit{Ad}}_i^{<n+1>}-{\mathit{Ad}}_i^{<n>}\right)}_{m,n}^{\left[k\right]}-\frac{{\mathit{Ad}}_{i,m,n}^{<n+1>,\left[k\right]}}{\mathrm{\Delta }t}\cdot {\left({\varphi }^{<n+1>}-{\varphi }^{<n>}\right)}_{m,n}^{\left[k\right]}+\dots \\ +q_{i,m,n}^{<n+1>,\left[k+1\right]}+\sum _j\frac{u_{x,m,n}^{j,\left[k+1\right]}\mathrm{\Delta }t}{2{\varphi }_{m,n}^{<n+1>,\left[k+1\right]}}\cdot {\left(\frac{\partial \mathit{Ad}}{\partial x}\frac{\partial \varphi }{\partial t}+\mathit{Ad}\frac{{\partial }^2\varphi }{\partial t\partial x}+\frac{\partial \varphi }{\partial x}\frac{\partial \mathit{Ad}}{\partial t}\right)}_{i,m,n}^{<n+1>,\left[k\right]}+\dots \\ +{\sum }_j\frac{u_{y,m,n}^{j,\left[k+1\right]}\mathrm{\Delta }t}{2{\varphi }_{m,n}^{<n+1>,\left[k+1\right]}}\cdot {\left(\frac{\partial \mathit{Ad}}{\partial y}\frac{\partial \varphi }{\partial t}+\mathit{Ad}\frac{{\partial }^2\varphi }{\partial t\partial y}+\frac{\partial \varphi }{\partial y}\frac{\partial \mathit{Ad}}{\partial t}\right)}_{i,m,n}^{<n+1>,\left[k\right]}-\dots \\ -{\sum }_j{\left(\frac{u_{x,m,n}^{j,\left[k+1\right]}\mathrm{\Delta }t}{2}\frac{{\partial }^2\mathit{Ad}}{\partial t\partial x}+\frac{u_{y,m,n}^{j,\left[k+1\right]}\mathrm{\Delta }t}{2}\frac{{\partial }^2\mathit{Ad}}{\partial t\partial y}\right)}_{i,m,n}^{<n+1>,\left[k\right]}+\dots \\ +{\sum }_j\frac{\mathrm{\Delta }t}{2{\varphi }_{m,n}^{<n+1>,\left[k+1\right]}}{\left[{\left(u_{x,m,n}^{j,\left[k+1\right]}\right)}^2\frac{{\partial }^2{V}_{i,m,n}^{j,\left[k\right]}}{\partial x^2}+{\left(u_{y,m,n}^{j,\left[k+1\right]}\right)}^2\frac{{\partial }^2{V}_{i,m,n}^{j,\left[k\right]}}{\partial y^2}+2u_{x,m,n}^{j,\left[k+1\right]}{u}_{y,m,n}^{j,\left[k+1\right]}\frac{{\partial }^2{V}_{i,m,n}^{j,\left[k\right]}}{\partial x\partial y}\right]}^{<n+1>}\end{array}$$E8

and the coefficients ${\mathcal{C}}_{\mathrm{1,2,3}}$ are calculated as follows:

$$\begin{array}{c}{\mathcal{C}}_1={\left({\varphi }_{m,n}+\mathrm{\Delta }t\frac{\partial \varphi }{\partial t}-\frac{u_{\mathit{tx},m,n}\mathrm{\Delta }t}{2{\varphi }_{m,n}}\frac{\partial \varphi }{\partial x}-\frac{u_{\mathit{ty},m,n}\mathrm{\Delta }t}{2{\varphi }_{m,n}}\frac{\partial \varphi }{\partial y}\right)}^{<n+1>,\left[k+1\right]}\\ {\mathcal{C}}_2={\left(\frac{{\varphi }_{m,n}}{\mathrm{\Delta }t}-\frac{u_{\mathit{tx},m,n}}{2{\varphi }_{m,n}}\frac{\partial \varphi }{\partial x}-\frac{u_{\mathit{ty},m,n}}{2{\varphi }_{m,n}}\frac{\partial \varphi }{\partial y}+\frac{u_{\mathit{tx},m,n}\mathrm{\Delta }t}{2{\varphi }_{m,n}}\frac{{\partial }^2\varphi }{\partial t\partial x}+\frac{u_{\mathit{ty},m,n}\mathrm{\Delta }t}{2{\varphi }_{m,n}}\frac{{\partial }^2\varphi }{\partial t\partial y}\right)}^{<n+1>,\left[k+1\right]}\\ {\mathcal{C}}_3={\left(1-\frac{\mathrm{\Delta }t}{2{\varphi }_{m,n}}\frac{\partial \varphi }{\partial t}\right)}^{<n+1>,\left[k+1\right]}\end{array}$$E9

2.2.2 Boundary conditions

At the beginning of the EOR process, the oil saturation in the reservoir is assumed to be equal to values after the application of primary schemes or the saturation after a waterflooding which reached the economical threshold at the producing well. Therefore, there is no chemical components present, and the initial pressure is constant throughout the reservoir. Thus:

$$t=0;\forall \left(x,y\right)\in \mathrm{\Omega }:z_c=0;z_p=S^{\mathit{\text{or}}H};p^a=p_i$$E10

The combined EOR process begins with the injection for a certain period of time of a polymer and/or nanoparticles at a constant concentration. After this period, the chemical slug is followed by a water-bank in order to sweep the remaining oil. As the boundary conditions in the domain, a “no flow” is assumed on the oil field contour ($\mathrm{\Gamma }$), since it is assumed that the porous rock is surrounded by an impermeable rock layer. Regarding the advection mechanisms, this is satisfied imposing a zero mobility on the boundaries. As far as the diffusive mechanisms are concerned, Fick’s law is applied rendering

$$\mathbf{\text{Injecting}}\mathbf{well}\Rightarrow \left\{\begin{array}{l}0\le t\le t_{\mathit{\text{in}}}:z_c=z_{\mathit{\text{in}}}\\ t>t_{\mathit{\text{in}}}:z_{\mathit{\text{in}}}=z_w,z_c=0\end{array}\right.$$E11

$$\mathbf{\text{Boundaries}}\Rightarrow {\lambda }_{m,n}^j=0\wedge \frac{\partial z_i}{\partial {\stackrel{⌣}{n}}_{\mathrm{\Gamma }}}=0;i=p,c;\forall t\wedge \forall \left(m,n\right)\in \mathrm{\Gamma }$$E12

3. Physical properties

The properties for a polymer flooding are well described in the literature [12]. These are considered in this model, but during this chapter only the novel approaches will be discussed and presented, including the polymer architecture and the nanoparticles’ modifications. Hence, the following phenomena are not included in the scope of this chapter: residual saturations, disproportionate permeability reduction (DPR), inaccessible pore volume (IAPV), polymer degradation, and capillary pressure.

3.1 Chemical component partition

As mentioned in the previous section, the conservation laws must be complemented with a number of extra relationships to mathematically find the solution of the system of equations. These come from what is known as the system’s phase behavior. This is considered essential in order to understand how the components distribute among the phases. In this model it is assumed that both oil and water remain in the organic and aqueous phases, respectively. Moreover, polymer and salt remain only in the aqueous phase. As far as the phase behavior of the nanoparticles is concerned, a similar concept to the one commonly used for surfactants is employed in order to calculate how the (nano)particles distribute between the phases. Depending on the particles’ wettability, namely: hydrophobic and lipophilic (HLPN); neutral wet (NWPN); or lipophobic and hydrophilic polysilicon nanoparticles (LHPN), these will be present in the organic, both, or the aqueous phase, respectively. Thus, three dimensionless relationships are used to account for this partition:

$$\text{Solubilization}\text{coefficient}=L_{\mathit{pnp}}^a=\frac{V_p^a}{V_{\mathit{np}}^a}$$E13

$$\text{Swelling}\text{coefficient}=L_{\mathit{wnp}}^o=\frac{V_w^o}{V_{\mathit{np}}^o}$$E14

$$\text{Partition}\text{coefficient}=k_{\mathit{np}}=\frac{V_{\mathit{np}}^o}{V_{\mathit{np}}^a}$$E15

Since in this model petroleum and water remain in their respective phases, the swelling and solubilization coefficients are zero $\left(L_{\mathit{pnp}}^a=L_{\mathit{wnp}}^o=0\right)$. The value of the nanoparticle’s partition coefficient is then determined by the physical characteristics of the nanomaterial used in the flooding, namely, HLPN (for $k_{\mathit{np}}\gg 1$), LHPN (for $k_{\mathit{np}}\approx 0$), and NWPN (for $k_{\mathit{np}}\approx 1$). Finally, the polymer and salt are assumed to remain completely in the aqueous phase, hence $V_s^o=V_{\mathit{pol}}^o=0$.

3.2 Phase viscosities

The main functionality of the polymer addition is to increase the rheological and viscoelastic properties of the aqueous phase, improving the recovery efficiency and avoiding the water fingering phenomenon [6, 26, 27]. Although the nanoparticles also affect the viscosity, its influence is not as important as with the polymer. Several correlations have been proposed to account for the influence of nanoparticles on the viscosity, since Einstein’s original work, used for low concentrations, to expanded studies considering, among others, the size and type of the nanoparticles, temperature, and the characteristics of the carrier fluid [36, 37, 38, 39, 40]. In order to consider these components in the calculation, a stepwise procedure is adopted. Since the polymer’s architecture is also considered in this model, it is necessary to determine the viscosity ratio between a linear polymer used as reference and a (hyper)branched one which has the same total molecular weight [41, 42, 43, 44, 45, 46]. Eq. (16) is used in the proposed model to calculate this relationship (Figure 5).

Figure 5.

Viscosity ratio $g_{\mathit{\text{viscosity}}}$ as a function of the numbers of arms and the molecular weight relationship [12].

$$g_{\mathit{\text{viscosity}}}=\frac{1}{{\left(1+\rho \cdot f\right)}^3}\left[1+2f\cdot \rho +\left(2f+f^2\right){\rho }^2+\left(3f^2-2f\right){\rho }^3\right]$$E16

where $f$ is the polymer’s number of arms and $\rho$ is the relationship between molecular weights of arms and backbone. The number of arms depends on the polymer used, with some authors reporting polymers for EOR up to 17 arms [47]. Subsequently, the aqueous phase viscosity at zero shear rate can be calculated using the Mark-Houwink relationship based on the polymer’s molecular weight:

$${\mu }_{0\mathit{sr}}={\mu }_w\left\{1+\left[k_1\left(g_v,f\right)\cdot \left[\eta \right]\cdot V_{\mathit{pol}}^a+k_2(g_v,f)\cdot \left[\eta \right]\cdot V_{\mathit{pol}}^{a^2}+k_3(g_v,f)\cdot \left[\eta \right]\cdot V_{\mathit{pol}}^{a^3}\right]C_{\mathit{SEP}}^{\mathit{Sp}}\right\}$$E17

The influence of nanoparticles is considered using a relationship for the nanofluid relative viscosity, expressed as a function of the particles’ diameter, which is a function of time (see Section 3.3), the carrier fluid molecules diameter, and the particles’ volumetric concentration [48, 49, 50, 51]:

$${\mu }_{\mathit{np}}={\mu }_{\mathit{cf}}\cdot \left[1+\left(2.5+{\eta }_1{e}^{-{\mathit{Dia}}_{\mathit{np}}{d}_w}\right)\cdot V_{\mathit{np}}^a+\left(6.2+{\eta }_2{e}^{-{\mathit{Dia}}_{\mathit{np}}{d}_w}\right)\cdot V_{\mathit{np}}^{a^2}\right]$$E18

where ${\eta }_{1,2}$ are correlation constants. It is important to mention that in the present model, the particular effects of associating polymers are not considered [40, 52, 53]. It is therefore recommended that future studies should be performed in order to establish a relationship between the rheological properties and the presence of both nanoparticles and these types of polymers.

3.3 Aggregation of nanoparticles

The aggregation of nanoparticles is a well-known and studied phenomenon, reported by many experiments which show that nanoparticles may aggregate until a critical size is reached, followed by their sedimentation [54, 55, 56, 57, 58]. This can be defined as the formation of clusters by (nano)particles or when small clusters aggregate to form bigger ones due to the forces present in the system, which are explained, for instance, by the Derjaguin-Landau-Verwey-Overbeek (DLVO) theory (Figure 6). As the flooding process evolves in the rock formation, the gradual increase in the clusters’ size may provoke the sedimentation, affecting negatively the aqueous phase properties regarding the oil recovery efficiency. On the other hand, these clusters may also be split into smaller ones if their size reaches a critical limit (Figure 7) [59, 60, 61, 62, 63, 64, 65].

Figure 6.

Different forces at the molecular scale as a function of the particles’ distance [65].

Figure 7.

Scheme of the aggregation mechanisms of nanoparticles and their subsequent splitting [33].

The mathematical modeling of this phenomenon was done in this chapter using a simple ODE with exponential growing, using the particle size and an aggregation constant to take into account the change in the particles’ diameter as a function of time. This size is calculated for every time-step of the simulation since this affects also other related phenomena during the flooding process (see Section 3.5). Furthermore, a maximum possible size for the clusters is also assumed, above which it splits into two identical, smaller clusters [59]:

$$\frac{{\mathit{\text{dDia}}}_{\mathit{np}}}{\mathit{dt}}=K_{\mathit{\text{aggreg}}}\cdot {\mathit{Dia}}_{\mathit{np}}$$E19

3.4 Diffusion of nanoparticles

The diffusion process can modify significantly the particles’ transport and therefore their influence on the recovery process. In this chapter the influence of the polymers’ molecules in this process will be also considered, which has not been reported previously [66]. In the proposed model in this chapter, the starting point is the Brownian diffusion model calculated by the Einstein-Stokes equation (Eq. (20)), which is valid for low concentrations, and it is therefore corrected for different values of the particles’ volumetric concentration (Eq. (21)):

$$D_{\mathit{np}}^a=\frac{k_{B}T}{3{\mu }_{\mathit{cf}}\pi {\mathit{Dia}}_{\mathit{np}}{f}_{\mathit{\text{corr}}}}$$E20

$$f_{\mathit{\text{corr}}}={\left(1-V_{\mathit{np}}^a\right)}^{-6.55}$$E21

Then, this diffusion coefficient must be corrected in order to take into account the presence of the polymer molecules and their characteristics (architecture, composition, radius of gyration, and molecular weight). The radius of gyration depends on its structure, its chemical composition, and molecular weight. The calculation for a linear molecule is presented in Eq. (22) [67]:

$$R_{g,\mathit{\text{linear}}}=b_{\mathit{rg}}\cdot {\left(\frac{M_{w,\mathit{bb}}+f\cdot M_{w,\mathit{arm}}}{M_{w,\mathit{mon}}}\right)}^{0.588}$$E22

where $b_{\mathit{rg}}$ is the segment length and $M_{w,\mathit{mon}}$ is the monomer’s molecular weight. Using this value, the radius of gyration of the branched polymer can be calculated using a similar procedure to the one used in the viscosity calculation, expressing this relationship as a function of the number of arms (Figure 8):

Figure 8.

Radius of gyration ratio $g_{\mathit{rg}}$ as a function of the numbers of arms and the molecular weight relationship [12].

$$g_{\mathit{rg}}=\frac{3f-2}{f^2}$$E23

The nanoparticles’ diffusion coefficient is then calculated using a stepwise function, expressed as a function of the particles’ size, radius of gyration, and overlapping concentration, based on the study made by Flory [68, 69, 70]:

$${\varphi }^{\ast }=\frac{3\left(M_{w,\mathit{bb}}+f\cdot M_{w,\mathit{arm}}\right)}{4{\mathit{\pi \rho }}_{\mathit{pol}}{N}_{\mathit{Av}}{R}_g^3}$$E24

$$\chi =R_{g,\mathit{pol}}{\left(\frac{V_{\mathit{pol}}^a}{{\varphi }^{\ast }}\right)}^{-0.76}$$E25

$$D_{\mathit{NPp}}^a\Rightarrow \left\{\begin{array}{l}{D}_{\mathit{NP}}^a\cdot e^{{\alpha }_D{\left(\frac{{\mathit{Dia}}_{\mathit{NP}}}{2\chi }\right)}^{{\delta }_D}}\mathbf{\text{if}}\left(\frac{{\mathit{Dia}}_{\mathit{NP}}}{2}\right)\le R_g\\ D_{\mathit{NP}}^a\cdot e^{{\alpha }_D{\left(\frac{R_g}{\chi }\right)}^{{\delta }_D}}\mathbf{\text{if}}\left(\frac{{\mathit{Dia}}_{\mathit{NP}}}{2}\right)\ge R_g\end{array}\right.$$E26

where ${\rho }_{\mathit{pol}}$ is the polymer’s density, $N_{\mathit{Av}}$ is the Avogadro number, and ${\alpha }_D$ and ${\delta }_D$ are empirical constants.

3.5 Retention and adsorption

The adsorption takes place when nanoparticles or polymer molecules are deposited onto the surface of the rock formation. This phenomenon causes a loss of the EOR agents in the porous media, decreasing the efficiency and eventually making the whole process unprofitable, in case of high adsorption rates [48]. Similar to what is reported in combined chemical EOR processes, the presence of two or more agents will cause a competitive adsorption process, modifying their deposition rates, since the polymer molecules will cover part of the rock’s surface, rendering a smaller available area for the adsorption of nanoparticles. This is modeled using two factors, one affecting the polymer’s adsorption, and the other the nanoparticles, which are function of the injection scheme. This also alters the porous media properties (i.e., porosity and permeability) if the particles’ size is larger or of the same order than the pore size or if large volumes of particles accumulate [71, 72, 73, 74, 75]. The loss is quantified in Eq. (4) and it can be calculated as [66, 72].

$${\mathit{Ad}}_{\mathit{np}}=\left(F_{\mathit{SP}}\cdot v_1+v_2+v_3\right)$$E27

where $v_1$ is the volume of nanoparticles available on the pore surfaces, $v_2$ is the volume of nanoparticles entrapped in the pores of the porous medium due to plugging and bridging, and $v_3$ is the release rate from pore walls by colloidal forces, considering the salinity and the nanoparticles’ possible charge [66]. The first term is usually expressed as a function of the critical velocity for surface deposition. Below this, only the retention phenomenon occurs, and above this value a combination retention and entrainment takes place:

$$v_1^{<n>}\Rightarrow \left\{\begin{array}{l}{\alpha }_1\cdot u^j\cdot V_{\mathit{np}}^a\mathbf{\text{if}}{u}^j\le u_{\mathit{\text{crit}}}^a\\ {\alpha }_1\cdot u^j\cdot V_{\mathit{np}}^a-{\alpha }_2\cdot v_1^{<n-1>}\cdot \left(u^j-u_{\mathit{\text{crit}}}^j\right)\mathbf{\text{if}}{u}^j\ge u_{\mathit{\text{crit}}}^a\\ v_1^{<0>}=0\forall \mathrm{\Omega }\end{array}\right.$$E28

where ${\alpha }_1$ is the coefficient for surface retention and ${\alpha }_2$ is the coefficient for entrainment. The critical velocity can be expressed as a function of the particles’/clusters’ size, which is a function of time (Figure 9) [71, 72, 76]. It is important to note here that this velocity is also a function of the particles’ shape, which can vary considerable in the case of clusters:

Figure 9.

Critical velocity (a) and retention parameters (b) as a function of the particle size [33].

$$u_{\mathit{\text{crit}}}^j\left[\raisebox{1ex}{$m$}\!\left/ \!\raisebox{-1ex}{$\mathit{day}$}\right.\right]=0.00992736\cdot {\mathit{Dia}}_{\mathit{NP}}\left[\mathit{nm}\right]+0.0009936$$E29

Regarding the term $v_2$, a formulation similar to the retention term in $v_1$ is adopted, expressed according to the following stepwise function:

$$v_2^{<n>}\Rightarrow \left\{\begin{array}{l}{\alpha }_3\cdot u^j\cdot V_{\mathit{np}}^a\\ v_2^{<0>}=0\forall \mathrm{\Omega }\end{array}\right.$$E30

where ${\alpha }_3$ is the pore blocking empirical constant. Lastly, $v_3$ is calculated based on the salinity, expressed as the TDS in the rock formation [66]:

$$v_3^{<n>}\Rightarrow \left\{\begin{array}{l}-{\alpha }_4\cdot v_1^{<n>}\cdot \left(C_{\mathit{sc}}-V_s^a\right)2\mathit{mm}\\ v_3^{<0>}=0\forall \mathrm{\Omega }\end{array}\right.$$E31

where ${\alpha }_4$ is the rate of colloidally induced mobilization. Eq. (31) implies that the colloidal release of particles from any phase $j$ is limited by a critical salinity $C_{\mathit{sc}}$, which is a function of the nanoparticles’ type and the mineralogy of the rock formation.

3.6 Change in absolute porosity and porosity

The processes studied in the previous sections alter, to a lesser or greater extent, the physical properties of the rock formation, i.e., porosity and absolute permeability. The nanoparticles’ sedimentation, adsorption, and retention (by blocking and bridging) affect the EOR process and are one of the main causes of formation damage. The porosity decreases with the nanoparticles’ deposition caused, among others, by aggregation mechanisms. Ju [72, 77, 78, 79] developed a numerical model to calculate the instantaneous porosity, which is independent of the alteration by compressibility:

$${\varphi }_{i,j}^{<n>}={\varphi }_{i,j}^{<n-1>}-\sum _{i=1}^3\left(v_i^{<n>}-v_i^{<n-1>}\right)$$E32

The absolute permeability is also affected by these processes, and its variation is related to the parameters calculated in the previous section, according to Eq. (33):

$$K_{i,j}^{<n>}=K_{i,j}^{<0>}{\left[\left(1-f_K\right)k_f+f_K\frac{{\varphi }_{i,j}^{<n>}}{{\varphi }_{i,j}^{<0>}}\right]}^{n_K}$$E33

where $K_{i,j}^{<0>}$ is the initial permeability, $k_f$ is a constant for fluid seepage allowed by plugged pores, and $n$ is an experimental coefficient which varies from 2.5 to 3.5. The term $f_K$ is the fraction of original cross-sectional area or unplugged pores open to the flow:

$$f_K=1-\sum _j{\gamma }_{\mathit{fi}}\cdot v_2^j$$E34

where ${\gamma }_{\mathit{fi}}$ is the coefficient of flow efficiency for the nanoparticles. The last part consists in describing how the nanoparticles also affect the rock’s wettability by the retention/adsorption processes described previously in this chapter and, hence, the relative permeabilities. Firstly, how much area of the rock is covered by nanoparticles and how much area is available should be calculated. This is based on experimental observations which evidenced a functionality of the type of nanoparticles, their shape and size, the rock formation, its origin and initial porosity, permeability, and wettability. It is assumed for this calculation that the nanoparticles are spherical and touching each other [13]. Thus, the specific area of the particles is

$$a_{\mathit{np}}=\frac{\mathit{\text{Area}}}{\mathit{\text{Volume}}}=\frac{6}{{\mathit{Dia}}_{\mathit{np}}}$$E35

Therefore, with the total amount of nanoparticles adsorbed or entrapped and the specific area, the total surface of the porous medium covered by the particles is calculated according to Eq. (36):

$$a_{t,i,j}=\frac{6\cdot {\beta }_{\mathit{np}}}{{\mathit{Dia}}_{\mathit{np}}}{\left(v_1+v_2+v_3\right)}_{i,j}^{<n>}$$E36

where ${\beta }_{\mathit{np}}$ is the surface area coefficient. To calculate the total specific area, the Kozeny-Carman equation is used. This correlation originally evidenced a problem in its applicability to all rock formations and, moreover, not valid for complex poral geometries [80, 81]. Thus, Carman modified the expression adding a variable to Kozeny’s equation, known as hydraulic tortuosity [82, 83]:

$$S_v=\frac{\mathit{Ko}\cdot {\varphi }^3}{T^2\cdot \sqrt{K_x^2+K_y^2}}$$E37

where $S_v$ is the total specific area, $\mathit{Ko}$ is the Kozeny constant, and $T$ is the tortuosity. This new variable was initially thought to be independent of the rock properties, but later it was proven to be a function of the porosity, with a minimum of 1 when the latter tended to unity, and then increasing monotonically as the porosity decreases (Figure 10). There are several relationships between tortuosity and porosity, for example, Eq. (38) based on the work presented by Matyka et al. [84]:

Figure 10.

Tortuosity-porosity relationship [12].

$$T=1-0.77\cdot ln\varphi$$E38

Using the results from previous equations, the fraction of the rock covered by the nanoparticles can be calculated, hence the wettability alteration (Figure 11). The reference condition is the relative permeability curves without nanoparticles in the system ($a_{t,i,j}=0$). On the contrary, when the entire surface of the rock is covered by the particles, the other extreme case is obtained $a_{t,i,j}/S_v\ge 1$, and then the maximum wettability change is reached, obtaining the new relative permeabilities $k_r^{a^{\ast }}$ and $k_r^{o^{\ast }}$. Between these two conditions means a partial coverage of the rock’s surface ($0<a_{t,i,j}/S_v<1$), adopting a linear interpolation between the extreme cases to calculate the relative permeability modification (Eq. (39)) [71, 72, 78, 79, 85]:

Figure 11.

Relative permeability modification based on the adsorption of nanoparticles [12].

$$k_{r,\mathit{NP}}^{a,o}=min\left(1,k_r^{a,o}\cdot \left[1+\frac{a_{t,i,j}}{S_v}\left({\theta }_{a,o}-1\right)\right]\right)$$E39

where ${\theta }_{a,o}$ represent the maximum achievable wettability modifications due to the nanoparticles’ adsorption (${\theta }_{a,o}=k_r^{a,o^{\ast }}/k_r^{a,o}$).

4. Results

The objective of this chapter, besides the development of a new CFD model for the simulation of a nanoparticle/polymer flooding process, is to present the potential advantages of combined techniques in EOR as well as the benefits of using CFD tools in order to optimize the recovery strategies in oil fields. Thus, a series of simulations were performed in a 2D oil field with a random permeability field (Figure 12), starting with a standard waterflooding process to be used as benchmark; followed by a traditional EOR polymer flooding using commercial, linear polymers to increase the viscosity of the sweeping phase; and finally a combined flooding with nanoparticles + branched polymers.

Figure 12.

Absolute permeability fields in the X (left) and Y (right) directions for the refined mesh, expressed in mD [12].

The recovery process was simulated during a period of 3000 days for the three different mechanisms. The results of both the recovery factor and the fractional flow as a function of time are presented in Figure 13.

Figure 13.

Oil recovery and fractional flow as a function of time for the reference cases and the nanoparticles and polymer scheme in the refined mesh [12].

At the beginning of the process, all the techniques render the same results, but the influence of the EOR agents is noticeable once the oil slug reaches the producing well. In the case of the linear polymer, there is a reduction of the oil in the field due to the increase in the viscosity. Furthermore, when the combined process with nanoparticles is used, these modify not only the rheology and viscoelastic properties of the polymer solution but also the wettability of the rock formation, allowing a further recovery of oil. This can be appreciated in Figures 14 and 15. The final oil saturation in the field decreases significantly when the nanoparticles are used, especially alongside the diagonal line connecting both wells, where the velocities reach the highest values.

Figure 14.

Oil saturation after 3000 days for the waterflooding (a), linear polymer (b), and the nanoparticle and polymer (c) EOR flooding schemes [12].

Figure 15.

Oil saturation after 3000 days for a linear polymer flooding (a) and a combined nanoparticles + polymer flooding (b) [12].

The influence of the nanoparticles is more evident in the areas with low permeabilities (Figures 14 and 15) in which both waterflooding and linear polymer could not desaturate completely. On the other hand, the polymer + nanoparticle flooding modified the wettability of the formation, increasing the mobility of the oil phase, rendering lower residual saturations.

Numerically speaking, the solution of a chemical EOR process implies solving a set of highly nonlinear and coupled balance equations. In the proposed simulator, an iterative approach for each time-step of the simulation has been proposed. This rendered a pseudo-implicit scheme, guaranteeing the numerical stability of the system. This behavior can be proved using a matricial stability analysis of the system [12].

5. Conclusions

The goal of this chapter was to present the current strategies in oil recovery and how traditional techniques can be boosted by means of the nanotechnology, introducing the development of a reservoir simulation using CFD techniques. There is a need of optimizing the production of conventional oil sources while more sustainable energy resources are developed and a smooth transition between these can be carried out. One of the techniques used to evaluate the performance of these methods is reservoir simulation, a branch of engineering that emerged in recent years, used to justify and analyze the execution of E&P investments. Among EOR processes, chemical agents show a great potential in different oil fields, being mostly used in low- and medium-viscosity fields. A way to improve their efficiency is to use the nanotechnology in order boost the advantages of these chemical agents.

A novel mathematical model of porous media flow for a combined EOR/nanotechnology process is presented during this chapter, using a (hyper)branched polymer with several possible architectures, coupled with nanoparticles of different wettabilities. The mathematical model is represented by the momentum (Darcy) and mass conservation laws, using a compositional approach due to its versatility to model multiphase, multicomponent systems. There are several physical phenomena present in EOR flooding, and the combination of chemicals and nanoparticles affects some of them, studied in this chapter, presenting a set of formulas to implement these in a reservoir simulator. The polymer architecture is key factor in the oil recovery, with branched (e.g., comb/star) polymers yielding better recovery factors than linear ones. On the other hand, nanoparticle flooding increases the oil recovered by altering the rock wettability, allowing the organic phase to flow more easily. Thus, the synergy between both agents presents a great potential for its application in field tests.

All in all, nanotechnology-enhanced chemical EOR flooding could represent a novel and improved technique, considering the advantages and synergy of the agents being injected. Nanotechnology represents a breakthrough in EOR processes, and it is a perfect example of how well-developed, standard techniques can be enhanced by using the advantages of materials exhibited at the nanoscale.

Nomenclature

$\mathit{Ad}$

component adsorption $\left[1/\mathrm{day}\right]$

$c_r$

rock compressibility $\left[1/\mathrm{Pa}\right]$

$\underset{¯}{\underset{¯}{\mathrm{D}}}$

dispersion tensor

$f$

number of arms (polymer)

$K$

absolute permeability $\left[\mathrm{mD}\right]$

$k_r$

relative permeability

$p$

reservoir pressure $\left[\mathrm{Pa}\right]$

$p_{\mathit{wf}}$

bottom-hole pressure $\left[\mathrm{Pa}\right]$

$q$

flowrate $\left[{\mathrm{m}}^3/\mathrm{day}\right]$

$r_w$

well radius $\left[\mathrm{m}\right]$

$S$

phase saturation

$s$

well skin factor

$u$

Darcy velocity $\left[\mathrm{m}/\mathrm{day}\right]$

$v$

nanoparticle adsorption

$V$

volumetric concentration

$z$

overall concentration

Greek letters$\mathrm{\Gamma }$

domain boundary

$\lambda$

phase mobility

$\mu$

absolute viscosity $\left[\mathrm{Pa}\mathrm{s}\right]$

$\sigma$

interfacial tension $\left[\mathrm{mN}/\mathrm{m}\right]$

$\varphi$

formation porosity

$\mathrm{\Omega }$

reservoir domain

Superscripts$a$

aqueous phase

$c$

capillary

$H$

water-oil system (no chemical)

$j$

phase

$<n>$

time-step

$o$

oil phase

$r$

residual

Subscripts$i$

component

$\mathit{\text{in}}$

injection

$m,n$

spatial grid blocks

$\mathit{np}$

nanoparticle component

$p$

hydrocarbon component

$\mathit{pol}$

polymer component

$s$

salt component

$t$

total

$w$

water component