Application of CFD technique to simulate enhanced oil recovery processes: current status and future opportunities

Nowadays, because of the reduction in oil resources and the passage of the first and second life period of current reservoirs, using enhanced oil recovery (EOR) methods is of great importance. In recent years, due to the developments in technology and the advent of powerful computers, using simulation methods in enhanced oil recovery processes is on the rise. The computational fluid dynamics (CFD) method, as a branch of fluid mechanics, is a suitable method for studying and simulating EOR methods. In this study, a review was done on the application of CFD studies for simulating EOR methods. Also, potentials for future studies and the challenges researchers may face in this method were mentioned. Although using this method in enhanced oil recovery processes has recently started, different areas for more studies still exist. To optimize the usage of this method in future studies, the necessity of multiphase models and solution methods development, as well as considering all microscopic parameters such as interfacial tension and viscosity in investigating oil recovery factor is of great importance.


Introduction
Most of the petroleum resources are in the third phase of their lives, and the number of unexplored resources for fossil fuels is limited. Therefore, it is necessary to increase the recovery of petroleum from current reservoirs. Extracting oil from petroleum reservoirs is done in three steps. The first step is performed on its own by natural forces in the reservoir. In the second step, water or gas will be injected into the reservoir to be studied. Nevertheless, more than 60% of the petroleum storage will remain in the reservoir (Zhang and Morrow 2006). Thus, in the third stage, called enhanced oil recovery, the oil recovery will be improved using different methods. It is believed that enhanced oil recovery technologies will play an essential role in the coming years due to the decrease in reservoir exploration and the heavy oil which remains in the reservoirs.
Enhanced oil recovery methods are technologies for extracting oil which is not producible with conventional methods (Babadagli 2003;Manrique et al. 2010;Pal et al. 2018;Park et al. 2018). These methods are categorized as follows: thermal methods (Bondarenko et al. 2017), chemical methods (Bera et al. 2017;Elyaderani and Jafari 2019;Fortenberry et al. 2016;Singh and Mahto 2017;Yousefvand and Jafari 2018), gas injection (Alagorni et al. 2015;Ren and Nguyen 2017;Salehi et al. 2014), and others (Bera and Babadagli 2015;Geetha et al. 2018;Kazemzadeh et al. 2019). Figure 1 shows different oil recovery methods (Kokal and Al-Kaabi 2010). Field studies and the investigation of enhanced oil recovery methods for optimizing these methods are almost impossible. Likewise, laboratory studies of these methods can be costly and dangerous. In addition, they would not provide an accurate insight into the problems which are hard to be investigated in a laboratory. Due to the reasons mentioned above, in recent years, researchers have investigated enhanced oil recovery methods by numerical methods and simulations.
One of the advantages of simulation methods is that all of the reservoir conditions such as temperature, pressure and different forces could be modeled. Another advantage is that they use results obtained from laboratory tests, whose scales Edited by Yan-Hua Sun * Arezou Jafari ajafari@modares.ac.ir 1 are smaller than those of field tests. As a result, when these simulations are scaled up to real conditions of the reservoirs, operational risks will decrease (Jafari 2008).
In this research, a comprehensive review of the studies using computational fluid dynamics (CFD) methods in different oil recovery enhancements has been carried out. Through categorizing each study done in various methods of enhancing oil recovery, governing equations, potentials, and the opportunities for using this technique were discussed. Finally, by critically reviewing previous studies, the challenges which researchers may face were discussed and some suggestions for future studies were offered.

Fundamentals of computational fluid dynamics
Computational fluid dynamics is one of the powerful computational methods that has been considered by many scientists and engineers in recent years (Fayal et al. 2010;Moradi et al. 2014). In the past, the CFD method was mostly used by mechanical and aerospace engineers; however, at the moment, chemical and petroleum engineers are exploring the benefits of this technique (Abdi et al. 2010;Raynal et al. 2016). CFD is a common method for modeling flow systems in many fields of chemical and petroleum engineering. It also allows engineers to better understand fluid flow calculations. The use of CFD as an appropriate tool for analyzing systems as well as studying fluid flow in hydrocarbon reservoirs is increasing. In fact, CFD simulation is a great tool for researchers to better understand the laws of hydrodynamics.
By examining specific factors, such as shear stress and fluid velocity inside the system, and their changes, it is possible to better predict their effects and achieve optimal values.
Despite experimental methods, CFD simulation enables visual observation of fluid phenomena like pressure and temperature distributions, velocity profiles, and phase displacement behavior in porous media from the pore scale to the field and reservoir scale. There are plenty of software packages available for CFD to reduce the time of data processing, but the response time is dependent on the complexity of the problem and the type of meshing.
In the modeling and analysis of hydrocarbon reservoirs, the study begins with the formulation of mass and momentum equations. To obtain more accurate information about the system, the momentum, mass, and energy conservation equations, if needed, must be considered simultaneously. Generally, the results of a CFD simulation include local information on how fluid moves in the reservoir and the calculation of recovery efficiency, to name but two.

CFD advantages
Three different methods, including experimental studies, theoretical approaches, and numerical simulations, are used for solving fluid flow fields. Empirical methods are based on experimental measurements and are usually based on Buckingham's Pi Theorem (Zohuri 2017). In these methods, by making equations dimensionless, a smaller sample can be used instead of a true model with real sizes and specific conditions. Then, by using the π-Buckingham theory, the results are generalized to the original model.
Theoretical or analytical methods have been developed based on solving governing equations in fluid mechanics and heat transfer. In most cases, the formulation of fundamental fluid mechanics and heat transfer rules are treated as second-order partial differential equations, which have only analytical and precise solutions in some exceptional states. This is because the governing equations in fluid mechanics create a set of nonlinear and dependent differential equations that must be solved with different initial and boundary conditions. Therefore, in most cases, the analytical solution to the equations of fluid mechanics is very limited, and these restrictions increase by applying boundary conditions. The method employed in recent years to solve these equations is the computational fluid dynamics, which is based on numerical computations. CFD advantages over experimental methods as follows: • Using CFD analysis makes it possible to examine the conditions that are impossible or dangerous to achieve under laboratory conditions. In other words, the study of special and critical conditions is possible with this 1 3 method. Also, in scaling up experimental results, CFD simulations could be very helpful. • The cost of performing CFD calculations is lower than that of an experimental test. In many cases, simulation has a faster rate than experimental work. For this reason, quite a bit of cost and time can be saved. • Using CFD, it would be possible to calculate all the parameters in each point of the current, but in the experimental method, only a few points can be obtained. To explain it more clearly, by using the CFD technique, complete information and very precise details can be obtained for solving the problem. • As new hybrid EOR methods have been developed, in order to apply these methods on a reservoir scale, it is important for companies to have a comprehensive knowledge of different aspects that affect the process. In such conditions, CFD simulations help to reduce risks, costs, and damages to the reservoir.
On the other hand, the CFD method has been used by many researchers in recent years and has been mentioned as a powerful method to study different issues in the oil industry. This method has many advantages in comparison with traditional methods. Usually, conventional software based on Black oil/Composition models and so on considered the porous medium as a black box and cannot solve the pore network modeling (PNM) problems well. Also, conventional reservoir software cannot analyze the physic of the system well. But CFD technique can examine all events in the porous medium like sediment formation and so on every moment and in details. This technique can solve microscopic problems like fluid flow in a glass micromodel, can model the connection between the pores appropriately, and can simulate the complex systems like nanofluid flooding and many other concerns. Scientific software based on CFD technique has good flexibility in using different equations like Darcy's law and Brinkman equation or new developed equation as well as different EOS to study the PVT behavior of phases. Also, in comparison with traditional software in the oil industry, CFD technique is able to simulate the mechanical equipment and instruments in order to design and make accurate equipment in oil production processes.

Governing equations
In every EOR process, different fluids with distinctive exist in the reservoir, which in the simplest state, can be oil and the injected fluid. In this condition, there is a two-phase flow regime in the porous medium. Also, in some conditionssuch as the existence of gas bubbles, solid particles, precipitations, and formation water-the complexity of the system will increase (Wu 2015). Therefore, to study this type of process, an appropriate multiphase model must be used. There are various models for investigating the multiphase flow regime (Chen et al. 2006). Eulerian-Eulerian and Eulerian-Lagrangian are two basic approaches for modeling multiphase flows. In the first approach, Euler basis is applied for all phases (without explicit calculation of the boundary layer between two phases), and in the other method, Euler basis is employed for continuous phase and Lagrangian basis is used for dispersed phases. It should be noted that, in previous studies on the simulation of enhanced oil recovery processes, Euler-Euler equations have been widely used. Later, according to the physics and assumptions of the problem, some common models applied in studies were presented (Ferziger and Peric 2012;Safaei et al. 2016;Shahmohammadi and Jafari 2014).

Mixture model
The mixture model solves the continuity, momentum and energy equations for the mixture. It also solves a volume fraction equation for the secondary phases. Microscopic governing equations of continuity, momentum, volume fraction and energy are as follows: where n is the number of phases, t is time, ⃗ F is body force, k is the volume fraction of phase k, S E includes any volumetric heat sources and the mixture velocity, density, and viscosity are as follows: The drift velocity of kth phase is: The slip velocity (relative velocity) is defined as the velocity of a secondary phase (p) relative to the primary phase (q) velocity: Also, for an incompressible phase E k = h k , ( h k is enthalpy for phase k) and for a compressible phase: In this study, oil (the main phase in the porous medium) and the injected fluid were assumed as the primary and secondary phases, respectively. The drift velocity which is related to the relative velocity is:

VOF model
The volume of fluid (VOF) model is applicable to immiscible multiphase flow regimes. This model can simulate a liquid-liquid multicomponent system by solving a set of momentum equations and using the volume fraction of each fluid ). The VOF model can also incorporate the effect of surface tension in the cross sections of both phases. Microscopic equations of continuity, momentum, and volume fraction are presented hereunder (Chung 2010).
The volume fraction equation is: where q is the volume fraction of phase q. The continuity equation for phase q is: where q , �� ⃗ v q , and ṁ qp represent density of phase q, velocity of phase q, mass transfer from phase q to phase p, respectively.
The momentum equation is: where and P represent viscosity and pressure, respectively. The VOF model is an appropriate multiphase model for studying wettability and surface tension. In Eq. 14, the effect of the IFT between oil and the injected fluid is included through F term proposed by Brackbill et al. (Brackbill et al. 1992): where is surface tension coefficient, is the curvature of the interface, n is the interface normal vector, and is the Dirac delta function. Also, to calculate wall adhesion, the static contact angle ( ) of the primary phase can be provided by the following equation: where n W and t W are the unit vectors normal and tangential to the wall, respectively. The energy equation is: where the mass-averaged energy is:

Eulerian-Eulerian model
Eulerian-Eulerian multiphase model is a complex but accurate method that allows for mixing and separation of phases.
It is an appropriate model to simulate droplets or bubbles of secondary phase(s) dispersed in the continuous fluid phase (as primary phase). Also, the Eulerian model can be used to define granular particles and their diameters. This model uses a single-pressure field for all phases and solves momentum, enthalpy, continuity, and turbulence equations for each phase, and additionally, it tracks volume fractions. In this model, heat and mass transfer between phases can happen (Jayanti 2018). The volume fraction of phase q is defined by: where n ∑ q=1 q = 1 The general conservation equations of this model are as follows: The continuity equation for phase q is: and conservation of momentum is where q is the qth phase stress-strain tensor, ⃗ R pq is interphase force, and ⃗ F lift,q is the lift force acting on a secondary phase.
Also, the conservation of energy can be where h q is the specific enthalpy of the qth phase, ⃗ q q is the heat flux, S q is a source term that includes sources of enthalpy (e.g., due to chemical reactions or radiations), Q pq is the intensity of heat exchange between the pth and qth phases, and h pq is the interphase enthalpy (e.g., the enthalpy of the vapor at the temperature of the droplets in the case of evaporation).

Average density and viscosity equations
The following equations are used to calculate the average properties of mixture fluids in a computational cell. These equations are based on volume fraction weighted average (Lv and Wang 2015).
where , and represent density, viscosity and volume fraction, respectively. Also, subscribes w and o refer to water and oil.

Boundary conditions
CFD simulations are based on numerical calculations. For this reason, adjusting a good initial condition is very necessary and can help the convergence of the solution. Almost all CFD problems are based on Navier-Stokes equations. Besides, the Navier-Stokes equations are usually a boundary value problem (Ferziger and Peric 2012). Thus, solving all CFD problems and fluid dynamics simulations depends on the correct application of the values of the variables in the boundary nodes, and boundary conditions define how the system operates.
The most common boundary conditions are categorized into inlet boundary conditions, outlet boundary conditions, wall boundary conditions, constant pressure boundary conditions, axisymmetric boundary conditions, symmetric boundary conditions, and periodic or cyclic boundary conditions. Choosing the appropriate boundary condition depends on problem perspective. In other words, depending on the dynamic or static conditions of the system, there are different boundary conditions. In each simulation of EOR process, we deal with the injection of a fluid from an injection well and production of oil from another well. Therefore, it is very important to set an appropriate boundary condition which models the entering/exiting fluid from a surface well. In fluid dynamics, a wide range of boundary condition types exists which permit a fluid to enter and exit the system. Table 1 shows different types of the boundary conditions for inlet and outlet of the solution domain.
Generally, at other surfaces or boundaries, no fluid enters nor goes out. Hence, the wall boundary condition, which is the most commonly used boundary condition in a CFD analysis, can be assigned. However, it should be noted that when the energy equation is solved (like thermal EOR methods), thermal boundary conditions must be defined at wall boundaries. Thus, the wall boundaries can be assumed as fixed heat flux, fixed temperature, convective heat transfer, external radiation heat transfer, as well as combined external radiation and convection heat transfer. Also, in other conditions, such as chemical or physical reactions, sediment formation, gas diffusion, mass transfer, and so on, different boundary conditions can be selected. Finally, it is worth mentioning that choosing any of the specified boundary conditions depends on the problem, solution domain and the purpose of simulation.

Application of CFD in enhanced oil recovery
So far, various studies have been conducted on the application of CFD technique in the simulation of EOR processes. Different software and multiphase models have been developed and presented to solve the governing equations of EOR processes. Also, these studies have been done on a different scale from microscopic to macroscopic scales. With the advent of modern methods of oil extraction such as nanofluid, polymer, thermal, and so on, some researchers have turned to use this technique in EOR investigations.  Highlights: Using a direct numerical simulation method to investigate the drainage process of CO 2 Zhao and Wen (2017) VOF Pore Fluid dynamics Fluent Nanofluid flooding Highlights: Pore-scale simulation of wettability and interfacial tension effects on flooding processes 1 3 shows a review of the studies made on numerical studies and simulations by the CFD method in EOR processes.
From Table 2, it can be found that different types of multiphase models like mixture, VOF, and Eulerian have been used to solve the multiphase flow in a porous medium, as well as a number of researches coupled the Navier-Stokes with Darcy's law equation. So far few studies have been conducted on a field scale and some of them did not consider heat and mass transfer phenomena in their simulations. In other words, most studies have investigated the process in terms of fluid dynamics physical point of view. Finally, it should be mentioned that different types of commercial software such as Fluent, COMSOL Multiphysics, ANSYS Polyflow, FOAMpro, and ANSYS CFX have been used to simulate the EOR processes.

Water injection
Water, as an inexpensive and most available fluid, is a good choice and the first candidate to be used in EOR processes. The water is mainly taken from seawater, surface water, as well as produced water, and then injected into the reservoir. The most important reason for injecting water into a reservoir is to maintain the reservoir's pressure so that it can continue the production of oil. This process generally has a low recovery factor and should be formed with high precision to prevent water production in the production well (Setiawan et al. 2014). Many parameters such as the injection rate, injection pressure, permeability, and porosity of porous media are important in this process. Therefore, these factors should be optimized so that we can have a high oil recovery factor in this process. Using CFD techniques, some researchers studied the oil recovery factor, wettability alteration, and so on in the water injection.
Silva et al. (Silva et al. 2017) conducted a numerical simulation of water injection in a homogeneous oil reservoir using ANSYS CFX software. Its geometry is a fraction of an oil reservoir (270 m × 180 m × 15 m) created by ICEM CFD software. The results showed that, by increasing the water injection rate from 45.8 to 137.5 m 3 /day, the recovery factor increased from approximately 12.64% to 19.11%. Also, they considered two points as the injection well: inside and at the top of the oil zone. The oil recovery of these two points showed a small difference. In comparison with injection at the top of the reservoir, the internal injection had a better recovery because, in this state, the injected fluid reached the production wells faster. It should be mentioned that, in order to have a high oil recovery factor in the water injection process, the point of injection is less important than the wettability of porous medium, the viscosity difference between oil and water, and the mobility ratio of the injected fluid. In addition, it was realized that the injected water tends to create a path toward the production well, which is an important issue in terms of the difference between water and oil. This could lead to coning phenomena that could cause problems in oil production (Silva et al. 2017).

Gas injection
One of the common methods to increase oil recovery is injecting gas. Some gases such as nitrogen and flue gas have more complicated injection conditions-including high pressure and injection in deeper levels-than CO 2 , is used more in EOR processes. Another problem with nitrogen injection is in re-separation reservoirs, preventing the remaining gas being used as fuel (Taber et al. 1997). It has been revealed by experimental results that CO 2 injection has more recovery than does N 2 (Tang et al. 2016). Therefore, the focus of researchers is on the simulation of CO 2 injection, using the CFD method. CO 2 injection is an accessible and low-cost method with environmental advantages. This method is one of the most common ways for increasing oil and gas recovery (Rafiee et al. 2018). As it is known, the behavior of CO 2 at 31 °C and 72.8 atm is supercritical, while at the normal temperature and pressure, it behaves like other gases (Agarwal et al. 2017). Therefore, injecting supercritical CO 2 is a promising method, which has been applied by many researchers. Tang et al. (2016) studied a new method for CO 2 -EOR, which focuses on huff and puff supercritical CO 2 injection. They made a comparison between continuous injection at high pressures for oil recovery from fractured tight reservoirs (Tang et al. 2016). They studied the effect of oil composition, permeability, CO 2 injection patterns, gas properties, and rock mineralogy on ultimate oil recovery factor.
Later, carbon dioxide injection tests on core scales were simulated by the CFD method. A trial and error procedure was employed for estimating the efficiency factor of carbon dioxide diffusion coefficient, which demonstrated the best agreement with oil recovery factor and oil saturation in the core samples. The simulated distribution of oil saturation versus time approved that carbon dioxide emissions happen from the fracture into the matrix or no-flow boundaries. As a result, it was deducted that the huff and puff cyclic method with the pressure drop process causes a high cumulative oil recovery and increases oil production rate about 10%. It also consumes less CO 2 than does continuous injection method (Tang et al. 2016).
Zhao et al. (Zhao et al. 2015) studied the effect of supercritical CO 2 injection on pressure and temperature alterations by a simulation through Fluent software. Also, for phase change calculations of boiling liquid expanding vapor explosion (BLEVE), which Fluent cannot model, user-defined functions (UDF) were imported to Fluent. It was observed that the pressure at the throat opening moment drops dramatically, and after that, the process of phase change became weak gradually. In other words, there were two stages in the process of changing pressure. The pressure increasing stage occurred at the moment of throat opening, and a small amount of CO 2 was produced in the phase change stage. When the pressure reached its maximum value, liquid boiling began and became weaker gradually. However, the temperature changes were reported as stable and uniformly decreased throughout the process (Zhao et al. 2015).
In Safi et al. (2016) and Agarwal et al. (2017) studies, numerical simulations of subsurface flow in the EOR process were performed, using COZSim/COZVeiw as a multiphase flow solver. Moreover, for Enhanced Gas Recovery (EGR), TOUGH2, which is a CFD solver, was utilized. By employing a genetic algorithm (GA) code, COZSim/COZ-Veiw and TOUGH2 were modified. The flowchart of the GA-integrated COZSim solver is presented in Fig. 2. In order to optimize the CO 2 injection rate, GA-COZSim/GA-COZView and GA-TOUGH2 were utilized for the two cases of constant mass and pressure injections. In the simulated reservoir, such parameters as porosity, reservoir temperature and pressure, oil gravity, and boundary conditions were considered (Agarwal et al. 2017). Furthermore, the effect of different injection rates on the EOR was investigated.
It is shown in Fig. 3 that, with increasing the injection rate, the recovery factor sharply increased initially, and after the injection rate of 4600 Mscf/day, it became almost constant. Also, in the EGR optimizing process, the CO 2 injection rates led to a shorter well shut-in time, and the methane recycling rate increased dramatically.
Understanding the phenomenon of CO 2 displacement on a microscopic scale is an important topic in examining its  (Safi et al. 2016) impact on the EOR. Zhu et al. (2017) studied this phenomenon, using the direct numerical simulation (DNS) method in an oil-wet porous medium. The DNS method is generalizable to different models, which is an important feature. In their study, the DNS method was coupled with COMSOL Multiphysics for the numerical solution of Navier-Stokes and phase field equations. The phase field method is an approach for studying diffuse interfaces in a multiphase system (Badalassi et al. 2003). The results showed that, after the carbon dioxide broke through in the outlet, the pressure in the mainstream of CO 2 gas flow was considerably reduced. Then, the oil began to flow in the larger pores, and increasing the viscous forces was the most important mechanism for improving oil recovery. Also, gravity fingers increased the CO 2 sweep area when the viscose forces were low, and, as shown in Fig. 4, as the contact angle decreases, oil recovery declines as well. Xing et al. (2013Xing et al. ( , 2014 examined the consequences of CO 2 environmental impacts and damages which would be done by possible blowout in the CO 2 -EOR process. Several CO 2 dispersion tests were performed on an experimental scale since it is not impossible to do them on a larger scale. After that, by using a series of scale-up rules, effective parameters-such as length, velocity, and time, achieved from the experimental tests-were generalized to real conditions. In order to validate the scale-up processes, numerical simulations were conducted by the CFD method. The k − ε, RNG k − ε, and SST k − ε turbulence models were utilized for modeling heavy gas dispersion. By Fluent software and transient species transport modeling, a simulation was performed to calculate the CO 2 concentration as a function of time. Figure 5 shows the comparison between the experimental and simulation results at different CO 2 flow rates. The results showed that k − ε and SST k − ε were in agreement with the observed results from experiments, and the RNG k − ε results, because of its poor statistical performance, were undesirable.
Finally, the CO 2 concentration at the centerline of the porous medium was measured. The results of the simulation represented an appropriate agreement with the experimental results, except for near the jet nozzle as the source of CO 2 emission at the bottom of the prism. It was concluded that the scale-up process (as much as 10 times larger) offered an acceptable result at low flow rates (Xing et al. 2013(Xing et al. , 2014.

Polymer flooding
Polymer flooding increases the sweep efficiency by adjusting the mobility ratio of injected fluid at an optimum level and, as a result, increases the oil recovery factor up to 30%. One of the important elements that make polymers useful for EOR processes is their ability in increasing the viscosity of flooding liquid. Also, choosing the best type of polymer with considering several factors is the most important step to apply in EOR applications (Zerkalov 2015).
The viscoelastic polymer is one of the new methods used for enhanced oil recovery processes. A fluid that includes viscoelastic materials such as hydrolyzed polyacrylamide (HPAM) is a non-Newtonian fluid that has both the characterizations of viscous liquid and those of elastic solids. It can be said that the viscoelastic fluid has a memory, because, after the deformation, it partially returns to its original state (Clarke et al. 2016).
CFD technique is a general method to simulate different chemical systems. CFD software like COMSOL Multiphysics is a powerful tool to design, predict the performance and analyze multiphysics systems. This technique has a great potential to simulate non-Newtonian and viscoelastic fluid flow in porous media as well (Zhong et al. 2017). Using computational fluid dynamics leads to the simulation of the viscoelastic flow around the oil droplets, resulting in a better understanding of sustainable reduction where Q is the volumetric flow rate, is the relaxation time of polymer molecules, and R t is the tightest radius of pore throats. They found that for De ≪ 1, elastic effects were negligible. In their study, as shown in Fig. 6, two geometries were simulated. The first one was pore center, and the second one was constriction. In the pore center geometry, the viscoelastic fluid pressure profile was not symmetric. This asymmetry increased with the increase in Deborah (25) De = Q πR 3 t dimensionless number because, unlike the Newtonian fluids, viscoelastic fluids required a relaxing time for returning to their original shape. Therefore, as it is presented in Fig. 7, the velocity profile for a viscoelastic fluid is not symmetric, and the velocity near the center of the geometry is higher than that of other places. Despite the difficulty of simulating the constriction model due to the fluid contraction at the tightest radius and high elasticity in pore throats, they conducted a constriction model for oil droplets whose size was larger than the diameter of the pore. They found that, by increasing the Deborah dimensionless number, the overall pressure drop increased and that a small pressure drop occurred in the front of oil droplets. Finally, they considered constant pressure boundaries instead of fixed rates in the constriction geometry. They concluded that, by increasing the elasticity of the fluid, normal stress as well as flow resistance increased, and the velocity profile in the throat was lower than that in the center. Also, they observed that, by increasing the Deborah dimensionless number, pressure forces reduced. Additionally, the improvement in oil mobility, which was due to the presence of viscoelastic polymer, resulted from combined effects, two of which were normal forces and additional pressure drop. These two occurred due to the increased viscosity. Another impact of combined effects was pulling the oil droplet mass, which resulted from the asymmetry in the flow lines. By using computational fluid dynamics, it was concluded that normal forces were very important for viscoelastic fluids rather than for Newtonian fluids. They also found out that the effect of normal forces increased significantly by increasing Deborah dimensionless number (Afsharpoor and Balhoff 2013;Afsharpoor et al. 2012).
High viscosity non-Newtonian polymer injection can cause excessive pressure drop through the well. Accurate models are needed to study and investigate this pressure drop. Jackson et al. (Jackson et al. 2011) examined a new approximation model by using computational fluid dynamics, which was developed for the pressure drop in horizontal wells. A real horizontal well and the modeled boxed section are presented in Fig. 8. The correction to pressure ( P * ) relationship equation of the developed model is presented in Eq. 26, which is based on some previous studies done by Marshall and Trowbridge (Marshall and Trowbridge 1974); this equation has some corrections for non-Newtonian fluids and high rates.
(26) where * new = k w R w t res is the modified conductivity ratio, R w is the well radius, z is the length of the tube, and which includes pressure drawdown ( P 0 − p T , P 0 is entry pressure, p T is outlet pressure), viscosity of fluid at the well wall ( well ) , conductivity ratio ( * new ), well radius ( R w ), and mean inlet velocity ( V 0 m ). In their study, the results of numerical solutions by COM-SOL Multiphysics software showed a noticeable agreement with the EA Marshall's analytical results (Marshall and Trowbridge 1974). Then, in order to assess the impact of the pressure drop in well sweep efficiency using UTCHEM software, a series of tests were performed based on a CFD model to determine the effect of the axial pressure drop of fluid injection with high viscosity. These tests also were conducted for homogeneous and heterogeneous wells. Finally, it was concluded that a significant pressure reduction (over 10%) happened in efficiency at early periods, i.e., 10 days, in homogeneous reservoirs. Also, with the passage of time, the reduction in oil recovery factor became slower, and the axial pressure drop was more in heterogeneous reservoirs than that in homogeneous reservoirs (Jackson et al. 2011).
The pore shapes and the connection between pore throats are crucial parameters in the investigation of the process on microscopic scales, one of which is pore network modeling (PNM). In these cases, microscopic geometry software such as CFD software is required for the creation of porous media. To do so, a precise method that can be applied is to use the SEM images obtained from thin sections of rock samples. In the study of Clemens et al. (Clemens et al. 2013), the polymer was injected into a micromodel, showing that the presence of polymer improves the mobility ratio. Then, a CFD simulation was performed using SEM images from a Berea sandstone thin section. Figure 9 shows the procedure of taking SEM images for digitizing and simulating a unit block with CFD technique. In this process, a quarter of the thin section was simulated, using CorelDraw and CATIA software. The final porous medium was constructed from putting together a series of previous simulated sections. Although it was not an accurate procedure, it was a simple way to create a porous medium based on a real section of the reservoir rock. The displacement processes were simulated on a microscopic scale, and the results showed a substantial agreement with laboratory data. CFD simulations were also employed to inject a polymer (FLOPAAM 3630) solution, and the results showed a considerable variation in viscosity of the polymer solution in a porous medium. It was observed that the displacement efficiency for assuming non-Newtonian polymer shear-thinning behavior was about 5% better than that for assuming Newtonian solution. This was because of the shear-thinning behavior of non-Newtonian polymer solution, which, as compared to Newtonian polymers, led to a higher viscosity in the pores.

Microbial enhanced oil recovery
Microbial enhanced oil recovery (MEOR) is a method which makes beneficial changes such as stable emulsions of oil and water. MEOR methods use microorganisms and their products to increase the oil recovery via different mechanisms, namely reducing IFT, moving oil to the top of the reservoir through blocking the permeable zones, and increasing the oil mobility (Laskin et al. 2003). Jafari et al. (2016) studied the impact of biosurfactant flooding on additional oil recovery (AOR), using the CFD method. They utilized COMSOL Multiphysics, as their CFD tool, for the creation of their micromodel. By investigating pressure drop profile, they found out that: at the inlet of the porous medium, the pressure was at its highest value; by moving through the micromodel, the pressure decreased; and, at the outlet, the pressure was at its lowest value.
They concluded that the biosurfactant injection increased the oil recovery due to reducing the pressure drop rate and lowering the fingering effect. Additionally, they observed that the lower shear force made improvements in the biosurfactant injection process and the oil recovery. Also, if the inlet velocity of the injection fluid increased, the AOR would reduce. Moreover, a reduction in the surface tension value between oil and surfactant (from 0.4 to 0.2 dyne/cm) enhanced AOR. Figure 10 presents the profiles of the studied parameters.  Gharibshahi et al. (2015) developed a method of computational fluid dynamics to simulate the effects of pore morphology and pore distribution in a two-dimensional micromodel on the enhanced oil recovery factor during the nanofluid injection. In Fig. 11, seven different types of micromodels with different patterns and shapes are illustrated. The effects of pore heterogeneity and pore conductivity with the presence or absence of the throat line and those of pore shapes on the enhanced oil recovery were investigated. The mixture model, as a simple multiphase model, was used to solve the governing equations. This model, which was able to select and calculate the properties of granular phases, was appropriate for the simulation of liquid-solid flows.

Nanofluid injection
The results showed that, in comparison with real conditions at field scale, in a model in which the throat lines connected pores, the oil recovery factor was more than that in the non-connected pores. To investigate the effect of pore shapes, two scenarios of the presence or absence of throats were considered.
In the first scenario, three micromodels with different forms of pores, spherical, quadratic, and triangular were selected. The oil recovery factor in the pores for the quadratic model quite matched the oil recovery factor in experimental data (Maghzi et al. 2012). In other words, the relative error in the quadratic model was smaller than that in other patterns. The second scenario considered models without throat lines. They chose two micromodels with circular and quadratic pore shapes. The results showed that the spherical pore shape model had a more accurate prediction than did the quadratic model .
Also, by simulating the homogeneous and heterogeneous models, the flow in the model which has distributed pores heterogeneously is similar to the flow in the reservoir rocks, although the trap effect in these models was not adequately investigated. In addition, because of the connection between pores, the fingering effect was reduced. This phenomenon was due to the similarity of the front flow. At the microscopic level with quadratic and triangular models, the pores which had corners were suitable for the investigation of the effect of the trapping fluid.  Dezfully et al. (2015) studied the effect of nanoparticles in the supercritical CO 2 injection for improving the oil recovery through the CFD method. For validating the obtained results, the simulation results were compared with experimental data (Maghzi et al. 2012), and the relative error of 5.17% was obtained, which was an acceptable value. It was observed that the presence of nanoparticles improved the oil recovery factor, and the oil recovery increased through the enhancement of nanoparticle concentration. It was found that, by increasing the volume fraction of nanoparticles in the base fluid (supercritical CO 2 ), the properties of the base fluid such as viscosity and density would change. Then, the mobility ratio of the injected fluid became adjusted at an optimum ratio; hence, the fingering effect phenomenon decreased. In conclusion, more surface area of the porous medium would be in contact with the injected fluid, and much oil would be produced. Figure 12 shows the result of the oil recovery factor at different concentrations of nanosilica. Zhao and Wen (2017) investigated the effects of wettability and interfacial tension in the scale of a pore on the enhanced oil recovery, using the VOF method. The finite volume method, based on ICEM CFD code, in Fluent software was utilized for the numerical simulation of Navier-Stokes equations. The effect of wettability on the flooding process was modeled in an unsaturated oil model, and their simulation presented a substantial agreement with the previous experimental data (Zhang et al. 2014). Based on previous studies (Mousavi et al. 2013), it was observed that nanoparticles could alter rock wettability. Different nanoparticles with different concentrations resulted in different ultimate wettability. Figure 13 illustrates six cases with different ultimate wettability in the process of oil/nanofluid flooding. In addition, they studied the effect of IFT changes, which were caused due to the use of surfactants, on the water flooding process. The results showed that the reduction in IFT between oil and water could be effective in better extraction of the oil phase. Furthermore, both wettability and capillary effects had significant influences on the EOR. Nandwani et al. (2019) conducted coreflood tests for surfactant flooding and applied the CFD technique to simulate the performed process. The used surfactant solution was a high salinity solution containing the active surface ionic liquid (C16mimBr) and the non-ionic surfactant (TERGITOL 15-S-9). ANSYS Workbench software was utilized to design the studied geometry, which was a calcite powder cuboidal packed bed, and the mixture multiphase model was used to solve the governing equations in ANSYS Fluent software. A comparison was made between the results of additional oil recovery gained from the experiment and those gained from simulation, and a substantial agreement between the results was observed. It was also observed that the surfactant mixture showed a higher efficiency due to the reduction of the fingering phenomenon and the decreased IFT.

Thermal methods
Thermal methods are used for increasing reservoir temperature, leading to a decrease in the viscosity of heavy oil or its evaporation part, which causes an increase in the oil recovery (Stahl et al. 1987).

SAGD
SAGD is a useful method for increasing the recovery of heavy oil bitumen. A correct and efficient injection of steam is crucial since producing steam is very difficult and costly. For optimizing this process, flow control devices (FCDs) are used. In designing FCDs, the safest output velocity for decreasing the damage to the casing is chosen. McChesney and Edlebeck (2015) showed how to determine the lowest bulk velocity by using CFD. To determine which geometry of closing sleeve could provide the best conditions for entering steam flow by FCDs, the two-phase CFD simulation for each case was done. In Design 1, it was concluded that recirculation, which occurs in the steam chamber, causes an adverse effect on the velocity rate. From Design 2, it was understood that full chamber length for the flow path in the closing sleeve was not the best condition for steam flow. Designs 3 to 5 showed less recirculation, since an offset slot was used in them. Simulation results showed that the casing would not have a collision in the steam flow velocities over 100 ft/s such that corrosion would not happen. Also, the results showed that, even if the closing sleeve was not installed, a 50% equilibrium between the left and right nozzles would be made. Figure 14 shows designs from 1 to 5.
In a study done by Quiroga (2018), numerical simulations were performed by ANSYS software, and the effects of oil saturation and contact angle changes in a square capillary were investigated with increasing the temperature. According to the distribution of steam in the porous medium, inflow control device (ICP) was also evaluated for steam injection. According to the results, it was found that, with increasing the size of the square capillary, the saturation of the residual oil increased. However, as the contact angle increased, the residual oil saturation decreased. Finally, it was observed that, as the temperature increased, the residual oil saturation decreased, which was in substantial agreement with the laboratory results (Argüelles-Vivas and Babadagli 2015). Gharibshahi et al. (2019) developed a CFD model, as a thermal EOR method, for simulating the simultaneous injection of metal nanoparticles and vapor in a two-dimensional glass micromodel. In this research, in order to construct the geometry of porous medium and solve the governing equations of the system, Gambit and Fluent software were used, respectively.
By using Taguchi test design method, four parametersincluding the type of nanoparticles, weight percent of nanoparticles in the base fluid, nanoparticle diameter, and input temperature of injecting fluid-were investigated. In addition, to study the displacement of the injected fluid in the porous medium, the mixture multiphase model was utilized, and the following equations were used to find the thermal properties of the oil: where T, K, and C are referred to temperature (in °C), thermal conductivity (in BTU °C −1 h −1 ft −1 ), and heat capacity (in BTU lbm −1 °F −1 ), respectively. Since the main mechanism of the oil recovery in the thermal EOR methods is to reduce the viscosity of the produced oil, the following equation was introduced into the software as a user-defined function to study the variations in the viscosity of the produced oil relative to the temperature.
ln (ln ( )) = 0.07547 + 5.76588 ln(API) − 0.00101(1.8T + 32) ig. 13 Final phase distribution at different wettability conditions (Zhao and Wen 2017) 1 3 where µ and T represent crude oil viscosity in cP and temperature in °C, respectively. In Fig. 15, the flow of the injected fluid in the model is illustrated. As it is evident in this model, the fingering phenomenon and the trapping of oil are appropriately modeled.
The obtained results are shown in Fig. 16, indicating that Al 2 O 3 nanoparticles, compared to CuO and Fe 3 O 4 nanoparticles, will have the greatest impact on the final oil recovery factor. Also, by reducing the diameter of the nanoparticles, increasing the weight percentage of the nanoparticles in the base fluid, and increasing the inlet temperature of the injected fluid, the oil recovery factor will increase. Finally, the temperature variations in the porous medium and the reduction of the produced oil viscosity under optimal injection conditions were obtained.
Argüelles-Vivas and Babadagli (2016) investigated heavy oil displacement by gas injection at high temperatures. Moreover, the residual oil saturation (S or ) on the capillary porous medium surface was evaluated, using CFD. Simulation results of the displacement at the temperatures of 55 °C and 85 °C were in reasonable agreement with laboratory data (Argüelles-Vivas and Babadagli 2015).
In their study, ANSYS CFX software was used to evaluate the homogenous multiphase model. In order to discretize the momentum transfer as well as the continuity and volume fractions equations, the VOF method was used. The residual oil saturation behavior at 200 °C was investigated, and, by increasing the temperature, a decrease in the amount of S or was observed at a constant air injection rate. In the last step, the 3D behavior of S or in the oil displacement by air in a square capillary was also investigated with a CFD method.   1 3 The results presented in Fig. 17 show that S or reduced with increasing the temperature and contact angle. Also, it can be seen that, in the thermal EOR method, in comparison with contact angle, temperature has more impact on the oil recovery factor.

Hot water flooding
Hot water flooding decreases the viscosity of oil due to the raising of the oil temperature such that it increases the oil mobility and improves the oil production (Kermen 2011). Hot water flooding is a thermal oil recovery method, which can be implemented for the special formations which are sensitive to freshwater. In comparison with other thermal methods, hot water flooding is a low-cost method (Torabi et al. 2012). Also, during hot water flooding, the pressure drop is low because the oil viscosity decreases due to the warming caused by hot water (Lv and Wang 2015). Lv and Wang (2015) tried to investigate the hot water injection process in a two-dimensional porous medium, using Fluent software as a CFD tool. They compared hot water with 100 °C temperature injection to the conventional water injection and observed that hot water increased the oil recovery factor about 8% more than does regular water. Figure 18 illustrates the difference in residual oil and water distribution between conventional and hot water flooding after water breaks through. Additionally, in their study, the mentioned process was investigated for oil-wet, water-wet, and intermediate wettability conditions. They found that, in the oil-wet system, the initial water and then a part of the oil were removed, and after that, the water breakthrough occurred. However, in systems with middle wettability, oil droplets and adjacent water were mixed, and then oil and water were removed periodically. In addition, it was shown that, in an oil-wet system, the oil recovery and the original oil saturation increased with raising the hot water temperature and decreasing the oil viscosity. For a system with intermediate wettability, three cases with three different characteristics concerning the droplets and symmetry of them were investigated, and it was concluded that the pattern of oil/water played an important role in the two-phase flow of water/oil system (Lv and Wang 2015).

Challenges and future prospects
The application of computational fluid dynamics techniques in the study of EOR methods is a considerable challenge. Currently, many researchers are inclined to use this method. Although an acceptable progress has already been made in this area, to fix the shortcomings of this technique, more comprehensive and profound studies are necessary. It is expected that, in the near future, there will be significant advances in this field.
Although CFD technique has been used for simulating several EOR processes, further studies on this subject should be done in the future. For example, several studies have already been done on the application of CFD method to simulate water flooding and gas injection, individually. However, unfortunately, there is no study of the simulation of WAG and SWAG operations with CFD technique. Also, due to the potential of CFD technique in simulating multiphysics problems, investigations on different EOR methods like combustion, electromagnetic heating, MEOR, and so on are strongly needed in future studies.
EOR methods require a complete understanding of the reservoir conditions and forces affecting fluid flow in porous media on the macroscopic and microscopic scales. In order to consider all the conditions of the system, most of the forces and processes, which determine their capabilities, such as phase change materials and deposition processes of nanoparticles in nanofluid injection process, do not exist in common simulation software. These conditions should be considered by entering user-defined functions (UDF) into simulators.
The processes of injecting fluid into porous media are generally time-consuming and complicated. It is acceptable that in the large-scale simulations, the microscopic effects and phenomena are less necessary. For example, in the field scale EOR simulations, the fluid flow on the pore scale and microscopic recovery are less important than sweep efficiency and macroscopic recovery on the reservoir scale. Then, in the large-scale simulations, the size of mesh becomes larger and the number of grids and nodes decrease, consequently. But, due to the dimensions of reservoirs and the number of grids, simulating the whole reservoir is complex. Therefore, because of memory and computational limitations, the entire well should be simulated by modeling the well in segments separately (Jackson et al. 2011).
In order to reduce computational costs and solve the governing equations of the system, appropriate multiphase models and methods should be developed for solving equations with respect to the limitations of computer memories and their processing power.
The wettability and interfacial tension alteration inside porous media are not yet fully modelable. For this reason, comprehensive equations should be developed in the future for considering these parameters in processes.
The investigation of microscopic effects such as the impact of ions on the low salinity and smart water injection processes as well as understanding the mechanism of oil recovery processes require precise examinations on molecular scales. Therefore, for investigating the interactions between atoms and molecules (liquid-liquid and solid-fluid interactions) based on the laws of physics, CFD software can be coupled with other techniques such as molecular dynamics.
So far, no specific and accurate relationship for considering the effect of all nanoparticles on the oil recovery factor has been presented. Thus, CFD technique allows Hot water flooding Conventional water flooding Fig. 18 Difference in residual oil and water distribution between conventional and hot water flooding after water broke through (Lv and Wang 2015) researchers to investigate nanoparticle injection processes without much cost to examine different effects such as the changes of diameter and shape of nanoparticles, concentration of nanoparticles in the base fluid, and so on. As mentioned above, the CFD technique has the potential to simulate the processes of the oil recovery on a reservoir scale. However, little research has been done so far. In large-scale simulation studies, minimizing disk space, memory and CPU usage, optimizing configuration for a parallel modeling, choosing the optimal computational domain, performing optimal meshing and reducing number of nodes and grids are necessary to increase computational speed (Gorobets 2016). Although, some options like finite volume method (FVM) is an appropriate approach in CFD codes to optimize memory usage and solution speed, especially for large-scale problems, it is hoped that these issues will be solved in the near future. On the other hand, it is necessary to use this technique to investigate the parameters involved before applying new EOR methods such as nanofluid injection on the reservoir scale. Hence, such parameters as the fracture distribution in porous media as well as the thickness and number of fractures can be studied.
One of the most important features of CFD is its ability to simulate mechanical equipment and instruments. Although it is possible that downhole equipment interacts with its surrounding area in the reservoir, studying this issue have not been considered in any previous studies, so far. Therefore, this technique can be used to design and make accurate equipment in oil extraction processes.
Optimization of a process is a method to make a system work better or more efficient by adjusting the effective parameters on the system at an optimal level without violating some constraint. It can minimize the cost and time while maximizing throughput and/or efficiency. By optimizing an EOR process, the amount of injecting materials will be minimized, all of devices and equipment will work optimally and the power consumption will be reduced, the risk of operation and formation damages will be reduced and finally more amount of oil can be recovered from the reservoir. Due to the advantages of CFD simulations such as reducing costs and time, obtaining complete information and very precise details as well as examining all conditions that are impossible or dangerous to achieve under laboratory conditions, this technique can be used wieldy in optimization of a process (Gharibshahi et al. 2019;Thévenin and Janiga 2008). Therefore, coupling of optimization and CFD method is a good candidate for screening the best operational conditions in each EOR process.
Finally, it should be noted that most studies have been conducted at ambient temperature and pressure. Therefore, in future studies, the effects of temperature and pressure during flooding operations should be investigated.

Conclusion
Given the advancement of computer equipment, the use of simulation techniques such as computational fluid dynamics techniques has become possible for the researchers. Using these methods leads to the reduction of the costs and risks related to experimental studies. Furthermore, more accurate and complete information could be gained from enhanced oil recovery processes.
In this research, different studies in the field of enhanced oil recovery processes, using CFD methods such as water injection; gas injection; chemical methods such as polymer flooding, MEOR, and nanofluid injection; as well as thermal processes, including SAGD and hot water flooding were reviewed.
Choosing an appropriate multiphase model to solve the governing equations is a critical step. Several multiphase models such as the mixture model, VOF, and the Eulerian model were used to obtain the most accurate results. In every CFD simulation problem, the grid independency of the results must be checked. Then, the validity of numerical results with experimental data should be examined. Also, different pieces of software based on CFD technique such as ANSYS Fluent, ANSYS CFX, and COMSOL Multiphysics, to name but three, should be applied to study the EOR processes.
Over the past decade, several researchers used CFD techniques to solve EOR problems; however, we are in the early stages of CFD simulations of enhanced oil recovery processes. With further developments of numerical models and considering more factors in future studies, better and more reliable results will be achieved. In conclusion, the obtained results can be used on a field scale and can provide more accurate predictions of enhanced oil recovery performances.
Finally, it should be mentioned that using CFD technique to simulate EOR processes is very young and this technique can be improved in the future. With further studies and the development of various equations based on this method, a very clear future for this method can be imagined in comparison with conventional models for using in oil industry problems.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creat iveco mmons .org/licen ses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.