# Effect of adsorbed polymers on relative permeability and capillary pressure: a pore scale numerical study

Patrick Barreau, Didier Lasseux, Henri Bertin, Alain Zaitoun

► To cite this version:

Patrick Barreau, Didier Lasseux, Henri Bertin, Alain Zaitoun. Effect of adsorbed polymers on relative permeability and capillary pressure: a pore scale numerical study. Carlos A. Brebbia; Henry Power; Tim A. Osswald; S. Kim. Transactions on Modelling and Simulation, 10, WIT Press, pp.549-556, 1995, Transactions on Modelling and Simulation, 978-1-85312-324-5. 10.2495/BE950611. hal-03850316

HAL Id: hal-03850316

https://hal.science/hal-03850316

Submitted on 13 Nov 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Effect of adsorbed polymers on relative permeability and capillary pressure: a pore scale numerical study

P. Barreau,ª D. Lasseux,ª H. Bertin,ª A. Zaitounᵇ

ªLaboratoire "Energétique et Phénomènes de Transfert", ENSAM, 33405 Talence Cedex, France

ᵇInstitut Français du Pétrole, 92506 Rueil-Malmaison, France

## 1 Introduction

Waterflooding is one of the most commonly employed method for oil recovery from underground reservoirs. Due to different factors such as rock heterogeneities and unfavorable viscosity ratio between oil and water, this technique becomes however sometimes inefficient. To remedy the excessive water mobilization, injection of hydrosoluble polymers has been used. The aim of this technique is to place a selective barrier which remains permeable to oil but lowers water mobility (Zaitoun & Kohler¹, Broseta & al.²).

During a polymer-solution injection, macromolecules are adsorbed on pore walls and three tendencies result (Zaitoun & Kohler³) : i) water mobility (ratio of relative permeability to viscosity) is reduced since water preferentially flows under the form of a film on polymer-coated pore-walls ; ii) oil mobility is almost not modified since oil flows in the center of the pores while water films induce a lubrication effect ; iii) due to pore-size diminution rather than a change of oil to water interfacial tension, the capillary pressure is significantly increased. Up to date, these tendencies have been appreciated from a qualitative point of view at a macroscopic scale of observation and for further investigations, microscopic pore scale studies are required to quantify water-mobility reduction as well as capillary pressure increase as a function of the lubricating layer thickness for instance

In this work, we initiated a two-dimensional numerical study of a two-phase flow in a cell of a periodic pore-like structure. Starting from a guess on the free surface separating the two phases, our goal is first to compute the evolution of this surface towards the stationary solution for the given wetting-fluid saturation. Due to the presence of the fluid-fluid interface, we found it adequate to base our numerical scheme on a boundary element technique as reported in Wrobel & Brebbia⁴ using constant approximations. Subsequently, cellular quantities such as permeabilities and capillary pressure, will be computed.
<figure><img src="image_3.png" /></figure>

## 2 Governing equations

### 2.1 Description of the model - Hypotheses

We consider a two-phase flow in a periodic two-dimensional divergent-convergent channel. The wetting fluid (β-phase) flows near the pore walls while the non-wetting one (γ-phase) flows in the center of the channel. We assume that the channel diameter is small enough for gravity forces to be negligible in comparison to viscous and capillary effects.

<figure><img src="image_4.png" /><figcaption>Figure 1. Unit cell.</figcaption></figure>

Both fluids are assumed to be newtonian, non-compressible and immiscible and the flow is supposed to be slow enough to obey Stokes approximation.

Once periodicity and horizontal symmetry are taken into account, we can limit our study to the unit cell depicted in Figure 1.

### 2.2 Dimensionless boundary value problem

According to our hypotheses, the physical process under consideration can be described by the following dimensionless boundary value problem :

$$ (P_{\alpha}) \begin{cases} \nabla p_{\alpha} = \nabla^2 v_{\alpha} \\ \nabla \cdot v_{\alpha} = 0 \end{cases} \qquad (1) $$

with the corresponding boundary conditions :

$$ (B.C.1) \qquad v_{\beta} = 0 \qquad \text{on } A_{\beta\alpha} \qquad (2) $$

$$ (B.C.2) \qquad v_{\beta} = \frac{\mu_{\beta}}{\mu_{\gamma}} v_{\gamma} \qquad \text{on } A_{\beta\gamma} \qquad (3) $$

$$ (B.C.3) \qquad \Sigma_{\beta} \cdot n_{\beta\gamma} = \Sigma_{\gamma} \cdot n_{\beta\gamma} + \frac{C\sigma}{h_{\gamma}l} n_{\beta\gamma} \qquad \text{on } A_{\beta\gamma} \qquad (4) $$

Since we consider the flow in an infinite succession of identical cells, two extra conditions -one on the velocity and one on the pressure- are derived from periodicity considerations to close the problem. These two conditions are :

$$ v_{\alpha}(\mathbf{r}+\mathbf{l}) = v_{\alpha}(\mathbf{r}) \qquad (5) $$

$$ p_{\alpha}(\mathbf{r}+\mathbf{l}) = p_{\alpha}(\mathbf{r})+\mathbf{l} \cdot \mathbf{h}_{\alpha} \qquad (6) $$

where **r** is the position vector relative to a fixed origin, **l** represents the characteristic vector for periodicity and, as in Eq. (4), **h**_α is a uniform pressure gradient prescribed in the α-phase. Here we have :
<figure><img src="image_5.png" /></figure>

$$l = l e_1 \qquad (7)$$

$$\mathbf{h}_{\alpha} = h_{\alpha} \mathbf{e}_{1} \qquad (8)$$

A close attention to the consequence of periodicity conditions indicates that the two pressure gradients must be equal ($\mathbf{h}_{\beta} = \mathbf{h}_{\gamma} = \mathbf{h}$), and this can be easily justified by considering the boundary condition in Eq. (4) at points P and P' (see Figure 1) along with the fact that curvatures are equal at these two points.

Once the periodicity condition is taken into account between entrances and exits of our unit cell in Figure 1, this readily gives:

$$\mathbf{v}_{\alpha}|_{A_{2\alpha}} = \mathbf{v}_{\alpha}|_{A_{1\alpha}} \qquad (9)$$

$$(\mathbf{n}_{\alpha} \cdot \mathbf{\Sigma}_{\alpha})|_{A_{2\alpha}} = \mathbf{n}_{\alpha} \cdot (\mathbf{I} - \mathbf{\Sigma}_{\alpha})|_{A_{1\alpha}} \qquad (10)$$

In these relationships, $\mathbf{v}_\alpha$, $p_\alpha$, and $\mu_\alpha$ respectively represent the velocity, pressure and dynamic viscosity in the $\alpha$-phase, $\alpha$ referring to either $\beta$ or $\gamma$; C and $\sigma$ are the curvature of $A_{\beta\gamma}$ and the interfacial tension. We have used the notation $\mathbf{n}_{\beta\gamma}$ for the unit normal vector to $A_{\beta\gamma}$, directed from the $\beta$-phase towards the $\gamma$-phase, $\mathbf{n}_\alpha$ for the outwardly directed unit normal vector to $\Omega_\alpha$ on $A_{k\alpha}$ ($k=1,2$) and $\mathbf{\Sigma}_\alpha$ for the stress tensor in the $\alpha$-phase defined by:

$$\mathbf{\Sigma}_{\alpha} = -p_{\alpha} \mathbf{I} + (\nabla \mathbf{v}_{\alpha} + {}^{1}\nabla \mathbf{v}_{\alpha}) \qquad (\alpha=\beta \text{ or } \gamma) \qquad (11)$$

The first boundary condition in Eq. (2) expresses the classical no-slip condition on the solid wall. The second one in Eq. (3) indicates that no mass transfer and no slip occur on the fluid-fluid interface while the third one in Eq. (4) follows from the exact balance between the normal stress jump and capillary forces on $A_{\beta\gamma}$.

In the above equations, we have made use of dimensionless quantities defined from the corresponding dimensional ones -with superscript *- by:

$$\mathbf{v}_{\alpha} = \frac{\mu_{\alpha}}{h l^{2}} \mathbf{v}_{\alpha}^{*} \qquad (12)$$

$$p_{\alpha} = \frac{1}{h l} p_{\alpha}^{*} \qquad (13)$$

$$C = C^{*} l \qquad (14)$$

while lengths have been made dimensionless by $l$.

At this point of the analysis, the problem is not time dependent although the initial condition on the free-surface is not the stationary solution of the flow. Assuming the flow to be very slow and made of a succession of equilibrium states, we make use of the classical kinematic condition which can be expressed as:
<figure><img src="image_6.png" /></figure>

$$ \frac{D\Phi}{Dt} = \frac{\partial\Phi}{\partial t_{\alpha}} + v_{\alpha} \cdot \nabla\Phi = 0 \quad (15) $$

to displace the free surface $A_{\beta\gamma}$ of equation $\Phi(x, y, t) = 0$ at time t.

## 3 Boundary element formulation

According to the physical boundary conditions, we chose to use a boundary integral formulation based on velocities and normal stress variables (Ramachandran⁵). Since the flow under consideration involves two different fluids separated by an interface, a sub-region partition technique has been employed (Huyakorn & Pinder⁶). A boundary integral form of the initial Stokes problem is first written in each sub-region $\Omega_{\beta}$ and $\Omega_{\gamma}$ and compatibility conditions (Eqs. (3) and (4)) between the two formulations are used in a second step along with the boundary and periodicity conditions to yield the final global system that provides velocities where normal stresses are prescribed and vice-versa.

### 3.1 Continuous integral formulation

Since entrance and exit surfaces $A_{1\alpha}$ and $A_{2\alpha}$ are artificial boundaries on which periodicity is applied, we only consider the interior Stokes problem for each domain $\Omega_{\beta}$ and $\Omega_{\gamma}$. Starting from the formulation derived elsewhere with both interior and exterior Stokes problems (Da Costa Sequeira⁷), it is easy to obtain the integral form, interior to $\Omega_{\alpha}$ (Lasseux⁸):

$$ \int_{\partial\Omega_{\alpha}} \{ \sum_{\alpha} (v_{\alpha}(\mathbf{x}), p_{\alpha}(\mathbf{x})) \cdot \mathbf{u}_{\mathbf{k}}(\mathbf{x}-\mathbf{x}_{0}) - \sum_{\alpha} (\mathbf{u}_{\mathbf{k}}(\mathbf{x}-\mathbf{x}_{0}), t_{\mathbf{k}}(\mathbf{x}-\mathbf{x}_{0})) \cdot \mathbf{v}_{\alpha}(\mathbf{x}) \} \cdot \mathbf{n} \, d\Gamma = 0 \quad k=1, 2 \quad (16) $$

where $\partial\Omega_{\alpha}$ is the boundary of $\Omega_{\alpha}$, $\mathbf{n}$ is the outwardly directed unit normal to $\partial\Omega_{\alpha}$, $\mathbf{x}$ is a field point and $\mathbf{x}_0$ the source point; $t_{\mathbf{k}}$ and $\mathbf{u}_{\mathbf{k}}$ are the solution of the fundamental Stokes problem :

$$ \left\{ \begin{array}{l} \nabla t_{\mathbf{k}} - \nabla^2 \mathbf{u}_{\mathbf{k}} = \mathbf{k} \delta(\mathbf{x} - \mathbf{x}_0) \\ \nabla \cdot \mathbf{u}_{\mathbf{k}} = 0 \end{array} \right. \quad (17) $$

and are given by (Ladyzhenskaya⁹) :

$$ u_{ij}(\mathbf{x}-\mathbf{x}_0) = \frac{1}{4\pi} \left[ -\delta_{ij} \ln(r) + \frac{(x_{0i} - x_i)(x_{0j} - x_j)}{r^2} \right] \quad i, j=1, 2 \quad (18) $$

$$ t_j(\mathbf{x}-\mathbf{x}_0) = \frac{1}{2\pi} \frac{\partial \ln(r)}{\partial x_j} \quad i, j=1, 2 \quad (19) $$

In Eq. (17) $\mathbf{k}$ refers to one of the base vectors $\mathbf{e}_1$ or $\mathbf{e}_2$ and $\delta(\mathbf{x}-\mathbf{x}_0)$ is the Dirac delta function, while, in Eq. (18), $\delta_{ij}$ is the Kronecker delta and :
<figure><img src="image_7.png" /></figure>

$$r^2 = |\mathbf{x} - \mathbf{x}_0| \qquad (20)$$

From a physical point of view, $u_{ij}$ represents the ith component of velocity at point $\mathbf{x}$ due to a point force (Stokeslet) in the j-direction applied at point $\mathbf{x}_0$ (Pozrikidis¹⁰).

## 3.2 Discrete form

Boundaries $\partial\Omega_\beta$ and $\partial\Omega_\gamma$ are discretized with constant elements and a special care is taken to obtain the same discretization on $A_{\beta\gamma}$ for connection purpose. The choice of constant elements was suggested by its fast implementation, its good accuracy, as pointed out by Sugino & Tosaka¹¹, and by the fact that tricky problems to precisely evaluate normal derivatives at corner nodes (where the normal is discontinuous) are avoided (Mitra & Ingber¹²). Boundaries of the physical domain are replaced by a succession of straight line segments of constant length with nodes placed in the middle. On each element, velocity and normal stress are assumed to be constant functions defined by the nodal value. The two integral equations (16) are discretized accordingly leading to a set of elementary integrals which are evaluated analytically (Lasseux⁸). At node i of the interface, the curvature is calculated from the mean value of the rate of change of the norm of the tangential vector between elements i-1 and i and elements i and i+1.

Once both domains are connected together and boundary, symmetry as well as periodicity conditions are taken into account, one clearly sees that the system is still singular and this sounds natural since no pressure reference has been specified. To remove the undetermination, we chose to set a null x-normal stress component on the first element of $A_{1\gamma}$ above the horizontal axis. This choice remains of course arbitrary.

## 3.3 Algorithm

Numerical experiments begin with a guess on the free surface $A_{\beta\gamma}$ that satisfies periodicity conditions. Elements and nodes are placed along all the boundaries and the following algorithm is used:

1. each elementary integral is computed for both $\beta$- and $\gamma$-phases, and the linear system is formed. To do so, horizontal symmetry, periodicity and boundary conditions are taken into account,

2. the system is solved using a Gauss elimination method,

3. the interface is updated according to the computed velocity field on the free-surface. This step is performed using the following discrete form:

$$\mathbf{r}_i^{n+1} = \mathbf{r}_i^n + \delta t_\alpha \mathbf{v}_{\alpha i}^n \qquad (21)$$

which represents an explicit Euler scheme of the first order of Eq. (15). In Eqs. (21) $\mathbf{r}_i$ is the position vector of interfacial points and $\delta t_\alpha$ represents the time step associated to the $\alpha$-phase. Superscripts n and n+1 denote variables at time steps n and n+1.
<figure><img src="image_8.png" /></figure>

4. The free surface is re-gridded keeping the initial element size (Sugino & Tosaka¹¹).

Steps 1. through 4. are repeated until the stationary solution is found for the interface. This solution corresponds to a configuration for which velocities are everywhere tangential to the interface and we check this by forcing the maximum difference, $\delta$, between the derivative at each point of the interface and the ratio $v_{\alpha} \cdot e_2 / v_{\alpha} \cdot e_1$ at the corresponding point to be smaller than a given value (typically, 1%). Precision of computations is checked through the divergence integral over each fluid domain. This quantity is less than $10^{-4}$ in our simulations.

# 4 Results

## 4.1 Validation

Our numerical tool was first experienced on a problem which corresponds to a pair of plane parallel infinite plates (see Figure 1 with $\alpha=\pi$) distant from 2R with a wetting coating film of thickness $e$. Assuming that this configuration is stable, it clearly represents an attractive test since the solution is known a priori. At entrances and exits of each phase, it is characterized by constant and linear x- and y-normal stress components respectively and a Poiseuille-like horizontal flow-pattern given by:

$$v_{\beta} = v_{\beta} e_2 = 1/2 (R^2 - y^2) \quad (22)$$

$$v_{\gamma} = v_{\gamma} e_2 = 1/2 ((R-c)^2 - y^2) - \mu_{\gamma}/2\mu_{\beta} ((R-c)^2 - R^2) \quad (23)$$

On the free surface, the velocity is constant and horizontal while on the solid wall, the x- and y- normal stress components are constant and linear respectively.

<figure><img src="image_9.png" /><figcaption>Figure 2. Velocities on entrances and exits on the test problem.</figcaption></figure>

simulate this kind of problem.

In Figure 2, we represented the two above theoretical data along with the corresponding numerical results obtained with $\mu_{\gamma}/\mu_{\beta} = 11$ and $e=0.7R$. The excellent agreement between the two approaches allow to conclude that our numerical scheme is well adapted to
<figure><img src="image_10.png" /></figure>

## 4.2 Example

<figure><img src="image_11.png" /><figcaption>Figure 3. Domain configuration.</figcaption></figure>

In this section, we present an example of our results obtained for the stationary interface in a periodic convergent-divergent channel as represented in Figure 1.

The experiment was performed with the following parameters : $\mu_{\beta}=10^{-3}$ Pa.s ; $\mu_{\gamma}=11 \times 10^{-3}$ Pa.s ; $\sigma=30 \times 10^{-3}$ N/m ; $h=1$ Pa ($l$ is the period) and $\delta t_{\beta}=0.1$ in the geometrical configuration depicted in Figure 3. The evolution of the interface

from the initial guess towards the stationary solution is depicted in Figure 4 for $t_{\beta}=0, t_{\beta}=1.4$ and finally $t_{\beta}=16.4$ for which $\delta \le 0.01$. In fact, as shown in Figure 5, the interface derivative and the ratio $v_{\beta} \cdot e_2 / v_{\beta} \cdot e_1$ match (to 1%) for this last computation time. Figure 4 clearly shows how the two fluid domains are reshaped to satisfy the prescribed forces.

<figure><img src="image_12.png" /><figcaption>Figure 4. Interface evolution from the initial guess towards the stationary solution.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 5. Interface derivative at ratio $v_{\beta} \cdot e_2 / v_{\beta} \cdot e_1$ at $t_{\beta}=16.4$.</figcaption></figure>
<figure><img src="image_14.png" /></figure>

**5 Conclusion**

In this work, we performed numerical experiments on a two-phase two-
dimensional Stokes flow in periodic convergent-divergent channels
representative of a pore structure of a porous medium. Taking into account
symmetry and periodicity, we showed how, from a guess on the initial interface,
we obtained the stationary solution with a numerical scheme based on a
constant boundary integral element method. This study represents an important
step towards subsequent computations of cellular permeabilities and capillary
pressure during multiphase flow in porous material and especially for the
comprehension of the role of polymer injection in porous media during oil
recovery.

**References**

1. Zaitoun, A. & Kohler, N., The role of adsorption in polymer propagation through reservoir rocks. SPE 16274, SPE Int. Symp.on Oil Field Chemistry, San Antonio, Texas, February 4-6, 1987.

2. Broseta, D., Medjahed, F. Lecourtier, J. & Robin, M., Polymer adsorption/retention in porous media: effects of core wettability and residual oil. SPE/DOE 24149, SPE/DOE 8th Symp. on Enhanced Oil Recovery, Tulsa, Oklahoma, April 22-24, 1992.

3. Zaitoun, A. & Kohler, N., Two-phase flow through porous media : effect of an adsorbed polymer layer. SPE 18085, 63rd Annual Technical Conference and Exhibition of the Society of Petroleum Engineers, Houston, Texas, October 2-5, 1988.

4. Wrobel, L.C. & Brebbia, C.A. (Ed.), 'Fluid flow', Section 4, *Free Surface Flow*, Pro. 1st Int. Conf. on Computational Modelling of Free and Moving Boundary Problems, Southampton, CML Publications, Southampton, Boston, 1991.

5. Ramachandran, P.A., *Boundary element methods in transport phenomena*, Computational Mechanics Publications, Southampton, Boston, 1994.

6. Huyakorn, P.S. & Pinder, G.F., *Computational methods in subsurface flow*, Academic Press, 1983.

7. DaCosta Sequeira, A., Couplage entre la méthode des éléments finis et la méthode des équations intégrales : application au problème de Stokes dans le plan, *Thèse Doc. Université Paris VI*, 1981.

8. Lasseux, D., Caractérisations expérimentale, analytique et numérique d'un film dynamique lors du drainage d'un capillaire, *Thèse Doc. Université Bordeaux I*, 1990.

9. Ladyzhenskaya, O.A., *The mathematical theory of viscous incompressible flow*, Gordon & Breach, New York, 1969.

10. Pozrikidis, C., *Boundary integral and singularity methods for linearized viscous flow*, Cambridge University Press, 1992.

11. Sugino, R. & Tosaka, N., Boundary element analysis of unsteady nonlinear surface wave on water, Boundary Element XII, Springer Verlag, New York, Pro. 12th Int. Conf. on Boundary Elements in Engineering, Sapporo, 1990.

12. Mitra, A.K. & Ingber, M.S., Resolving difficulties in the BIEM caused by geometric corners and discontinuous boundary conditions, Pro. 9th Int. Boundary Element Conf., Stuttgart, Computational Mechanics Pubs., Southampton, 1987.

*Acknowledgement*: Authors gratefully acknowledge financial support from the 'Groupement Scientifique ARTEP - CNRS' including 'IFP, TOTAL, ELF and Gaz de France', without which this work could not have been performed.