The first geometry is used to establish a simple model for a network of three caves, with water from the below aquifer and oil is displaced upwards, as shown in the following pictures.

Real parameters:

$$\rho_{water} /\rho_{oil}$$ (kg/m^{3}) |
$$\nu_{water} /\nu_{oil}$$ (Pa.s) | $$ \sigma_{oil/water} $$ (N/m) |

1000; 800 | 0.001;0.01 | 0.025 |

It is observed a lot of numerical diffusion on the interface with the real parameters, therefore it is difficult to distinguish the two phases,. A large area of the simulation would have the volume fraction of oil around 0.5, which is impossible because oil and water are not miscible.

This phenomenon can be explained by the fact that in a same cell small portion of oil and water are present, therefore the oil volume fraction can not be either zero or one.

To avoid this problem, several modifications can be done, including :

- refine the mesh or increase the time step, but this would increase the time of calculation

- increase the surface tension to avoid the emulsion of the fluid

-increase the viscosity to avoid the effect of the separation of the fluid in others smaller structures.

Firstly, the viscosity and the surface tension are increased in order to analyze if the results are the same as the one with the real parameter with the help of another software fluent.

The parameters of the simulation are:

$$\rho_{water} /\rho_{oil}$$ (kg/m^{3}) |
$$\nu_{water} /\nu_{oil}$$ (Pa.s) | $$ \sigma_{oil/water} $$ (N/m) | $$ dt_{min}/dt_{max}$$ (s) | inlet/outlet lenght (m) | inlet velocity (m/s) | mesh size (m) | mesh size (cells) | inlet Re |

1000; 800 | 1;10 | 0.487 | 1e-7; 1e-5 | 0.056; 0.025 | 0.2 | 0.76x0.44 | 122x88 | 10 |

This simulation shows no numerical diffusion, with the post treatment, the outlet flow rate and the oil remaining and recovery ratio will be analyzed versus time.

For the outlet flow rate:

$$ Q_{outlet}=\frac{1}{n}\sum \limits_{i=1}^n \tau_i v_i L_{outlet}$$

$$ratio_{oil-remaining}=\frac{oil_{initial}-Q_{outlet}t}{oil_{initial}}$$

$$ratio_{oil-recovery}=1-ratio_{oil-remaining}$$

with n: number of outlet cells

$\tau_i$ : oil volume fraction of cell i

v_{i}: outlet velocity of cell i

$$oil_{initial}=\sum \limits_{i=1}^m \sum \limits_{j=1}^l \tau_{ij} S_{cell}$$

m x l: numbers of cells of the domain

$S_{cell}$: area of a cell

To plot the data, the dimensionless number is applied, the flow rate reference is $Q_{inlet}$ and the time reference is $t=\frac{S_{domain}-S_{obstacles}} {L_{inlet} V_{inlet}}$, which is the time for water to fill the domain.

We have flow rate conservation due to the incompressibility of the fluids, therefore when only oil flow out the inlet and outlet flow are equal, after the dimensionless time around 0.3, water begins to flow out, so the oil outflow rate decreases and a change of the gradient in the oil remaining ratio plot which increases because less oil flows out in time.

Same parameters as the simulation with Jadim

$$\rho_{water} /\rho_{oil}$$ (kg/m |
$$\nu_{water} /\nu_{oil}$$ (Pa.s) | $$ \sigma_{oil/water} $$ (N/m) | dt (s) | inlet/outlet lenght (m) | inlet velocity (m/s) | mesh size (m) | mesh size (cells) | inlet Re |

1000; 800 | 1;10 | 0.487 | 1e-3; | 0.056; 0.025 | 0.2 | 0.76x0.44 | 122x88 |
10 |

The results are almost the same, the differences can be explained by the difference of the two codes the way to solve the equations are not the same.

The fluctuations observed on the outflow rate of oil are due to the bubbles of oil crossing the outlet.

Real parameters

$$\rho_{water} /\rho_{oil}$$ (kg/m |
$$\nu_{water} /\nu_{oil}$$ (Pa.s) | $$ \sigma_{oil/water} $$ (N/m) | dt (s) | inlet/outlet lenght (m) | inlet velocity (m/s) | mesh size (m) | mesh size (cells) | inlet Re |

1000; 800 | 0.001;0.01 | 0.025 | 1e-3; | 0.056; 0.025 | 0.024 | 0.76x0.44 | 244x176 |
1200 |

Decreasing the surface tension would increase the emulsion of the fluid and more bubbles would be created.

Compare with the last simulation, the viscosity has been divided by 1000, the velocity by 10 and therefore the Reynolds number obtained is 12000. If the same velocity was kept before the Re number would have been 120000, which give a turbulent flow. That is why the velocity has been multiply by ten.

Compare with the viscosity, the decrease of the velocity is negligible.

Therefore it is shown the influence of the viscosity comparing the two simulations.

Comparison of the two fluent simulations: effects of the viscosity.

The oil recovery is larger when the viscosity increases, because it creates big bubbles which are flowed out by water.

In oil industry, the chemical products are used to decrease viscosity in the case of porous medium reservoirs, which decreases the capillarity pressure and therefore helps oil to go up.

It is estimated in the case of the caves network reservoir, increase the viscosity with chemical products can be considered as a method of oil enhanced recovery ratio.

For the second geometry,the case of two caves of different sizes connected by a fracture is considered, with water injection from the bottom, thus oil is displaced towards the upside, as shown in the figure below.

The parameters such as Reynolds numbers, viscosity of the fluids, interfaces tension, as shown in the table below, will be varied to see the effects on the oil recovery, sweeping efficiency, distribution of the remaining oil,flow rate, etc. The CFD softwares Jadim and fluent, are used separately to simulate this case. The sensitivity of mesh density is also tested.

Simulation parameters

$$\rho_{water} /\rho_{oil}$$ (kg/m^{3}) |
$$\nu_{water} /\nu_{oil}$$ (Pa.s) | $$ \sigma_{oil/water} $$ (N/m) | $$ dt_{min}/dt_{max}$$(s) | inlet/outlet lenght (m) | inlet velocity (m/s) | mesh size (m) | mesh size (cells) | inlet Re |

1000;800 | 2;20 | 0.487 | 0.1e-7;0.1e-4 | 0.02;0.05 | 2 | 0.45*0.35 | 90*70 | 20 |

1000;800 | 2;20 | 0.487 | 0.1e-7;0.1e-4 | 0.02;0.05 | 5 | 0.45*0.35 | 90*70 | 50 |

1000;800 | 0.1;1 | 0.1;1 | 0.1e-7;0.1e-4 | 0.02;0.025 | 5 | 0.45*0.35 | 90*70 | 200 |

1000;800 | 2;20 | 2;20 | 0.1e-7;0.1e-4 | 0.02;0.025 | 2 | 0.45*0.35 | 360*280 | 20 |

1000;800 | 0.01;0.1 | 0.01;0.1 | 0.1e-7;0.1e-4 | 0.02;0.025 | 1 | 0.45*0.35 | 360*280 | 2000 |

Firstly, a coarse mesh and real parameters of fluids are applied to simulate the flow with VOF method by jadim software.

However, the simulation is not reasonable because the water and oil is mixed. This could be caused by the high Reynolds number or the density of mesh. Therefore, the inlet velocity is decrease and the viscosity is raised for the simulation.

$$\rho_{water} /\rho_{oil}$$ (kg/m^{3}) |
$$\nu_{water} /\nu_{oil}$$ (Pa.s) | $$ \sigma_{oil/water} $$ (N/m) | $$ dt_{min}/dt_{max}$$ (s) | inlet/outlet lenght (m) | inlet velocity (m/s) | mesh size (m) | mesh size (cells) | inlet Re |

1000;800 | 0.1;1 | 0.025 | 0.1e-7;0.1e-4 | 0.02;0.05 | 5 | 0.45*0.35 | 90*70 | 200 |

This time, the phenomenon of diffusion is not as severe. But the mixture of the two phase still happens.

Thus, to avoid the diffusion, we need to use:

-a higher value of viscosity (100 times) to avoid this situation.

-a refine mesh.

$$\rho_{water} /\rho_{oil}$$ (kg/m^{3}) |
$$\nu_{water} /\nu_{oil}$$ (Pa.s) | $$ \sigma_{oil/water} $$ (N/m) | $$ dt_{min}/dt_{max}$$ (s) | inlet/outlet lenght (m) | inlet velocity (m/s) | mesh size (m) | mesh size (cells) | inlet Re |

1000;800 | 2;20 | 0.487 | 0.1e-7;0.1e-4 | 0.02;0.05 | 2 | 0.45*0.35 | 90*70 | 20 |

1000;800 | 2;20 | 0.487 | 0.1e-7;0.1e-4 | 0.02;0.05 | 5 | 0.45*0.35 | 90*70 | 50 |

With the viscosity of the fluids 1000 times of the initial values and Reynolds within 100, the calculation finally reaches convergence. Below is the initial and final state for the Reynolds 50.

As mentioned before, to solve the problem of convergence during calculation and test the sensitivity of results to the mesh density. A mesh of density 16 times compared with the former one is applied.

To achieve the above objectives, two group of parameters are used for the refine mesh:

A. A group of parameters same with one of those used for the coarse mesh, which is:

$$\rho_{water} /\rho_{oil}$$ (kg/m^{3}) |
$$\nu_{water} /\nu_{oil}$$ (Pa.s) | $$ \sigma_{oil/water} $$ (N/m) | $$ dt_{min}/dt_{max}$$ (s) | inlet/outlet lenght (m) | inlet velocity (m/s) | mesh size (m) | mesh size (cells) | inlet Re |

1000;800 | 2;20 | 0.487 | 0.45*0.35 | 0.2;0.025 | 2/5 | 0.45*0.35 | 360*280 | 20/50 |

B. The interface tension as initial values. Viscosity 10 times and 100 times to the initial values.

$$\rho_{water} /\rho_{oil}$$ (kg/m^{3}) |
$$\nu_{water} /\nu_{oil}$$ (Pa.s) | $$ \sigma_{oil/water} $$ (N/m) | $$ dt_{min}/dt_{max}$$ (s) | inlet/outlet lenght (m) | inlet velocity (m/s) | mesh size (m) | mesh size (cells) | inlet Re |

1000;800 | 0.01;0.1 | 0.025 | 0.1e-7;0.1e-4 | 0.02;0.025 | 1 | 0.45*0.35 | 360*280 | 2000 |

1000;800 | 0.1;1 | 0.025 | 0.1e-7;0.1e-4 | 0.02;0.025 | 1 | 0.45*0.35 | 360*280 | 200 |

When the viscosity is 10 times, the problem of convergence still exists. However, for the case where viscosity is 100 times, the calculation achieves convergence. For the group of parameters same with one of those adopted with the coarse mesh, the calculation converges; however it take almost 14 times to finish the simulation.

For case B, the movement of the fluids is similar with what we have obtained from the coarse mesh. Also, the interface between the two fluids is clearer.

To compare the differences results given by fluent and jadim. Firstly, same parameters with one group which is adopted by jadim in utilized, as shown in the table below.

$$\rho_{water} /\rho_{oil}$$ (kg/m^{3}) |
$$\nu_{water} /\nu_{oil}$$ (Pa.s) | $$ \sigma_{oil/water} $$ (N/m) | inlet/outlet lenght (m) | inlet velocity (m/s) | mesh size (m) | inlet Re |

1000;800 | 2;20 | 0.487 | 0.02;0.05 | 5 | 0.45*0.35 |
50 |

With fluent, we have more dispersed phase due to the high viscosity. The differences for the results of jadim of fluent will be discussed later.

Then the real parameters of the fluids is used to see the differences.

$$\rho_{water} /\rho_{oil}$$ (kg/m^{3}) |
$$\nu_{water} /\nu_{oil}$$ (Pa.s) | $$ \sigma_{oil/water} $$ (N/m) | inlet/outlet lenght (m) | inlet velocity (m/s) | mesh size (m) | inlet Re |

1000;800 | 0.001;0.01 | 0.025 | 0.02;0.05 | 0.066 | 0.45*0.35 |
1320 |

For the real parameters of the fluids, the viscosity is much smaller, the velocity is small, The weber number is very small, so the interface tension plays an important role in the process, so the interface does not change.

**1. Evolution of flow rate and oil recovery**

In the very beginning before the dimensionless time reaches 0.3, the outlet oil flow rate equals with the inlet flow flow, because of mass conservation and that water has not reached the outlet. Afterwards, water reaches the outlet, the outlet oil flow rate starts to slow down, and the gradient for the oil recovery rate starts to slow down.

**2. Parameters sensitivity**

**A. Inlet Re**

Re=20 Re=50

To evaluate the effects of Reynolds number on the oil recovery, all other parameters are set the same value except that the Re number are different at the inlet. For a high Re, water reaches the outlet earlier, and the value of oil recovery is higher. This could lead by the recirculation caused by a higher Re.

**B. mesh density (0.00125m/cell & 0.005m/cell)**

*Coarse mesh (Final state) Refine mesh(Final state)*

From the figures of oil recovery ratio and flow rate, for the simulation with the refine mesh, the water reaches the outlet later, because the water diffuses less with the refine mesh. Therefore, we have a lower oil recovery ratio. And the interface in clearer with the refine mesh. However, it take about 14 times to finish the calculation with the refine mesh. Thus, both calculation duration and the precise requirements should be taken into consideration when choosing the mesh density.

**3. Comparison for fluent and jadim**

To compare the the differences of the softwares jadim and fluent, all parameters are set as the same value. From the figures, by using software fluent, water reaches the outlet with a shorter time. One possible reason is that there is more diffusion with fluent. Another reason for the difference could be the method used to calculate the flow rate at the outlet. In jadim, the oil outlet oil outlet velocity is calculated by the average of the multiplication of the velocity and oil ratio in every cell by every time step. However, for fluent, this velocity is calculated by the average velocity and average oil ratio by every time step.

The third geometry is used to established a more complex network of three caves with several connectivity using ibm method (immersed boundary method) with jadim. The aim of this geometry is to discover whether oil would be blocked or water would be produced early .

$$\rho_{water} /\rho_{oil}$$ (kg/m^{3}) |
$$\nu_{water} /\nu_{oil}$$ (Pa.s) | $$ \sigma_{oil/water} $$ (N/m) | $$ dt_{min}/dt_{max}$$ (s) | inlet/outlet lenght (m) | inlet velocity (m/s) | mesh size (m) | mesh size (cells) | inlet Re |

1000; 800 | 0.1;1 | 0.487 | 1e-7; 1e-3 | 0.5; 0.13 | 0.01 | 1.5x0.6 | 350x200 | 50 |

No oil has been blocked, but water begins to be produced at dimensionless time 1.3, where the curve of oil remaining and recovery ratio's gradients change.

Water production is not profitable for the industry, therefore "only" 80% of oil will be produced which is really important in comparison of the standard reservoir in porous medium in the world where the mean recovery ratio is around 35%.