© E. Flauraud & D.Y. Ding, published by IFP Energies nouvelles, 2020

Nomenclature

$A_{\mathrm{ICD}}$ : Surface area of ICD

$A_{\mathrm{ICV}}$ : Surface area of ICV

BHP: Bottom Hole Pressure

$C_{\mathrm{ic},l}^{\alpha }$ : Molar fraction of the component ic in the phase α in the reservoir grid block l

$C_{\mathrm{ICD}}$ : Constriction factor of ICD

$C_{\mathrm{ICV}}$ : Constant constriction factor of ICV

d_k : Diameter of the segment k

f_k : Friction factor of the segment k

$F_{\alpha ,k}^{\mathrm{mol}}$ : Molar fraction of the phase α in the segment k

$F_{\alpha ,m}^{\mathrm{mol}}$ : Molar fraction of the phase α in the ICD m

g : Gravity constant

${\mathrm{kr}}_{\alpha ,l}$ : Relatives permeability of the phase α in the reservoir grid block l

l_k : Length of the segment k

$P_k^a$ : Pressure in the node k of the annulus

$P_k^t$ : Pressure in the node k of the tubing

$P^{\lim }$ : Limit bottom hole pressure

$P_{r,l}$ : Pressure in the reservoir grid block l

${∆P}_k^a$ : Pressure drop in the segment k of the annulus

${∆P}_k^t$ : Pressure drop in the segment k of the tubing

${∆P}_k^{t,\mathrm{hydro}}$ : Hydrostatic pressure drop in the segment k of the tubing

$∆P_k^{t,\mathrm{fric}}$ : Friction pressure drop in the segment k of the tubing

$∆P_k^{t,\mathrm{ICV}}$ : Pressure drop in IVC of the segment k of the tubing

$∆P_m^{\mathrm{ICD}}$ : Pressure drop in the ICD m

$q_{\mathrm{ic},j}^{\mathrm{mol}}$ : Molar flow rate of the component ic between the reservoir and the well at the perforation j

$Q_{\alpha }^{\lim }$ : Limit flow rate of the phase α

$Q_{T,k}$ : Total volumetric flow rate in the segment k

$Q_{T,m}$ : Total volumetric flow rate in the ICD m

$Q_{w,k}^a,Q_{o,k}^a,Q_{g,k}^a$ : Volumetric flow rates of the water, oil, and gas phase in the annulus segment k

$Q_{w,k}^t,Q_{o,k}^t,Q_{g,k}^t$ : Volumetric flow rates of the water, oil, and gas phase in the tubing segment k

$Q_{w,m}^{\mathrm{ICD}},Q_{o,m}^{\mathrm{ICD}},Q_{g,m}^{\mathrm{ICD}}$ : Volumetric flow rates of the water, oil, and gas phase in the ICD m

${\mathrm{WI}}_j$ : Well index at the perforation j

$X_r$ : The set of unknowns in the reservoir grid blocks

$Z_{\mathrm{ref}}$ : Bottom hole reference depth of the well

${∆Z}_{k,k+1}$ : Height difference between the nodes k and k + 1

${\overline{\rho }}_k^{\mathrm{mass}}$ : Mass density of the fluid mixture in the segment k

${\overline{\rho }}_m^{\mathrm{mass}}$ : Mass density of the fluid mixture in the ICD m

${\rho }_{\alpha ,k}^{\mathrm{mass}}$ : Mass density of the phase α in the segment k

${\rho }_{\alpha ,k}^{\mathrm{mol}}$ : Molar density of the phase α in the segment k

${\rho }_{\alpha ,m}^{\mathrm{mol}}$ : Molar density of the phase α in the ICD m

${\mu }_{\alpha ,l}$ : Viscosity of the phase α in the reservoir grid block l

Subscripts and superscripts

m : Compartment or ICD index

M : Total number of compartments or ICDs

K : Node or segment index

$k_1^m,k_2^m$ : The first and the last node of the compartment m

N : Total number of nodes in the annulus and in the tubing

ic : Component index

ncons : Total number of components

n : Refer to the discret time tⁿ

n + 1: Refer to the discret time t ⁿ⁺¹

1 Introduction

Intelligent completions have become increasingly popular, especially for horizontal and multilateral wells, as the companies are striving to maximize the oil production and minimize the number of wells. Subsequently, Inflow Control Devices (ICDs) are widely equipped in the well in order to equalize wellbore pressure drop to achieve an evenly distributed flow profile along the well or to control the production in some particular branches. Active control devices such as Autonomous Inflow Control Device (AICD) or ICV are also more and more used to optimize the oil recovery by reducing the undesirable production of gas and water.

The technique of inflow control devices was introduced in the 90s (Brekke and Lien, 1992) to adjust the pressure drop along the well and to balance the inflow along the length of the completion. This technique has also the ability to improve sweep efficiency, delay water or gas breakthrough by reducing the localized drawdown and redistributing the influx over a long wellbore length, by taking into consideration reservoir heterogeneities and pressure drops across the completion. Later on, active control devices such as AICD were developed to optimally choke back unwanted fluids. ICVs, which can be controlled from the surface in order to reduce undesired fluid production, are also greatly used. In recent years both ICDs and ICVs have gained popularity and are being applied to a wide range of fields. Their efficiency has been studied through a variety of field cases (see, e.g., Al-Enezi et al., 2012; Al-Khelaiwi and Davies 2007; Al-Thuwaini et al., 2009; Broni-Bediako et al., 2019; Das and Al-Enezi, 2014; Dimitios et al., 2009; Gurses et al., 2019; Javid et al., 2018; Raffn et al., 2007).

Various types of intelligent well completions have been proposed to improve well performances and ultimate recoveries (see, e.g., Al-Khelaiwi, 2013). These completions are capable of managing the fluid flow into or out of the length of the wellbore in order to optimize the well production. Among the advanced completions, an ICD is a sandface completion technology specifically developed to help balance the flow contribution along the length of the wellbore. It is a device that restricts or chokes the inflow of fluids from the screen section. The nozzle type of ICD is the most commonly used. The size and the number of nozzles are designed to balance the inflow profile along the wellbore. The original ICD concept had a number of labyrinth channels installed within a pre-packed screen mounted on a solid base pipe. Other types of ICDs such as helical channels, slots, tubes, nozzles, and orifices, etc. have also been developed by various suppliers. Autonomous ICDs generally involve a more complex flow path to generate more flow restriction with higher pressure drops when flowing fluids are lighter or less viscous than oil.

Appropriate ICD/ICV completion design is the main factor for successful completion performance over the entire production life of wells. In many cases, a simplified semi-analytical type approach is used to describe the pressure drop through ICDs for various design configurations (see, e.g., Das and Al-Enezi, 2014; Lauritzen and Martiniussen, 2011). During the design phase, the fluid flow along the well is taken into consideration in the modeling, and sometimes, annular flow is also considered. However, the inflow computation from the reservoir is generally simplified, using an analytical or semi-analytical approach (Das and Al-Enezi, 2014; Das et al., 2012). This kind of approaches can be considered suitable for a rapid design, which focuses only on the immediate near-well region and needs a fast response. However, to evaluate long-term production, the model should be coupled with a numerical reservoir simulator, which takes into account reservoir heterogeneities, aquifers and gas cap, as well as at the presence of other injection/production wells.

From reservoir study point of view, there is also a need to model ICDs with a reservoir simulator to predict long-term well productions. The integration of the reservoir simulator with a wellbore model and intelligent completions forms an essential part of optimal reservoir management. In most reservoir simulators, this kind of intelligent wells is modeled using the multi-segment well model (see, e.g., HoJeen and Dogru, 2007; Holmes et al., 1998; Holmes et al., 2010; Stone et al., 2015; Wan et al., 2008; Youngs et al., 2009). In the multi-segment well model, the wellbore is discretized by several segments that have their own characteristics (length, diameter, inclination, roughness, etc.) and may or may not be connected to one or more reservoir grid blocks. Besides, additional segments can easily be tuned to represent special devices such as ICDs. Compared to the standard well model, the multi-segment approach allows the well geometry to be represented more accurately and to better describe fluid flows in the wellbore. However, solving implicitly a multi-segment well model coupled with the reservoir model is often complex and penalizes the simulation time, and an explicit decoupled modeling is usually not stable or convergent.

In this paper, we present a semi-implicit approach to integrate ICD/ICV model in a reservoir simulator with a simplified multi-segment concept. The annulus is discretized with a node in each perforated grid cell. A Well Index (WI) is used to link the reservoir and the annulus. A wellbore flow model is used to simulate various pressure losses, including friction, acceleration, and hydrostatic pressure drops, inside the tubing along the well. ICDs are used to link the annulus and the tubing, and ICVs are generally placed inside the tubing. The reservoir model and the well flow model can be two standalone modules. The proposed approach is based on a locally implicit scheme to make the connection between these two models. In most cases, the pressure drop across an ICD/ICV is generally much higher than that in the annulus or tubing. So the modeling in the annulus can be simplified. We will first present the numerical approach of this semi-implicit scheme, and then present the numerical examples, which illustrate the effectiveness of this coupling strategy in terms of production efficiency and computational performance.

2 The well model

2.1 The well discretization

The configuration of the multilateral well with ICDs considered in this paper is illustrated by the Figure 1. The annulus part of well branches is subdivided into isolated compartments with packers. Each compartment is equipped with ICDs. We consider that in an annulus, there is only one flow path for the fluid to go into the tubing through ICDs. Besides, for the multilateral wells, it is convenient to add in the tubing Inflow Control Valves (ICVs) to balance the flow rates from laterals and then optimize the production (Das and Al-Enezi, 2014).

At the discrete level, the annulus and the tubing are discretized into several segments. Each segment is delimited by two neighboring nodes which are symbolically represented by the center of the perforated grid blocks. In the following, we denote by N the number of well perforations and by M the number of compartments. The nodes in each compartment are numbered from $k_1^m$ to $k_2^m$ (from toe to heel) where $k_1^m$ and $k_2^m$ are respectively the index of the first and last node of the compartment m. We still denote by k, for $k=1,N-1$, the segment delimited by the nodes k and k+1.

An ICD is represented through a connection between a node of each compartment to its immediate neighbor in the tubing. Concerning the ICVs, they are associated to an existing tubing segment and add an additional pressure drop in the segment. Finally, an additional top segment connects the last node of the tubing to the well’s bottom hole reference depth $Z_{\mathrm{ref}}$ where we define the Bottom Hole Pressure (BHP). In the following, we take the example of nozzle type ICD in the modeling, and only under-critical conditions are considered for the ICV.

2.2 The pressure drop calculation

At the time scale of the reservoir, we suppose that the flow is stationary in the wells. Then the pressure losses in each segment in the well are given by the following steady-state momentum equations:

• In the annulus, we consider only hydrostatic pressure loss:

$$P_{k+1}^a-P_k^a={-∆P}_k^a\left(Q_{w,k}^a,Q_{o,k}^a,Q_{g,k}^a\right)=-∆P_k^{a,\mathrm{hydro}},\hspace{1em}k=k_1^m,\cdots ,k_2^m\hspace{1em}\mathrm{and}\hspace{1em}m=1,\cdots ,M,$$(1)

• where $P_k^a$ and $P_{k+1}^a$ are the pressures at the perforations k and k+1 of the annulus, $∆P_k^{a,\mathrm{hydro}}$ the hydrostatic pressure difference on the segment k and $Q_{w,k}^a,Q_{o,k}^a,Q_{g,k}^a$ are the volumetric flow rates in the segment k.

• In the tubing, we take into account the hydrostatic pressure loss, the frictional pressure loss and eventually the ICV pressure loss:

$$P_{k+1}^t-P_k^t=-{∆P}_k^t\left(Q_{w,k}^t,Q_{o,k}^t,Q_{g,k}^t\right)=-\left(∆P_k^{t,\mathrm{hydro}}+∆P_k^{t,\mathrm{fric}}+∆P_k^{t,\mathrm{ICV}}\right),\hspace{1em}k=1,\cdots ,N,$$(2)

$$\mathrm{BHP}-P_N^t=-{∆P}_N^t\left(Q_{w,N}^t,Q_{o,N}^t,Q_{g,N}^t\right)=-\left(∆P_N^{t,\mathrm{hydro}}+∆P_N^{t,\mathrm{fric}}+∆P_N^{t,\mathrm{ICV}}\right),$$(3)

• where $P_k^t$ and $P_{k+1}^t$ are the pressures at the nodes k and k+1 of the tubing. $∆P_k^{t,\mathrm{ICV}}$ is the pressure loss through the ICV placed in the segment k of the tubing.

• In the ICD, the momentum equation is given by:

$$P_k^t-P_k^a=-{∆P}_m^{\mathrm{ICD}}\left(Q_{w,m}^{\mathrm{ICD}},Q_{o,m}^{\mathrm{ICD}},Q_{g,m}^{\mathrm{ICD}}\right),\hspace{1em}k=k_2^m\hspace{1em}\mathrm{and}\hspace{1em}m=1,\cdots ,M.$$(4)

The different pressure losses are defined as follows:

• $$∆P_k^{\mathrm{hydro}}={\overline{\rho }}_k^{\mathrm{mass}}g{∆Z}_{k,k+1},$$(5)

• where g is the gravity constant, ${∆Z}_{k,k+1}$ is the height difference between the nodes k and k+1, and ${\overline{\rho }}_k^{\mathrm{mass}}$ denotes the mass density of the fluid mixture in the segment k defined by,

$${\overline{\rho }}_k^{\mathrm{mass}}={\displaystyle \sum _{\alpha \in \left\{w,o,g\right\}}}{\beta }_{\alpha ,k}{\rho }_{\alpha ,k}^{\mathrm{mass}},$$(6)

with ${\rho }_{\alpha ,k}^{\mathrm{mass}}$ the mass density of the phase α in the segment k and ${\beta }_{\alpha ,k}$ the volume fraction of the phase α: ${\beta }_{\alpha ,k}=\frac{Q_{\alpha ,k}}{Q_{T,k}}$ and $Q_{T,k}=Q_{w,k}+Q_{o,k}+Q_{g,k}$.

• $$∆P_k^{\mathrm{fric}}=\frac{8l_k{f}_k{\overline{\rho }}_k^{\mathrm{mass}}{Q}_{T,k}^2}{{\pi }^2{d}_k^5},$$(7)

where $l_k$, d _k and $f_k$ are respectively the length, the diameter and the friction factor of the segment k.

• $$∆P_k^{\mathrm{ICV}}=\frac{{\overline{\rho }}_k^{\mathrm{mass}}{Q}_{T,k}^2}{2A_{\mathrm{ICV}}^2{C}_{\mathrm{ICV}}^2},$$(8)

where $A_{\mathrm{ICV}}$ is the surface area of the ICV and $C_{\mathrm{ICV}}$ is the constant constriction factor of the valve.

• $$∆P_m^{\mathrm{ICD}}=\frac{{\overline{\rho }}_m^{\mathrm{mass}}{Q}_{T,m}^2}{2A_{\mathrm{ICD}}^2{C}_{\mathrm{ICD}}^2\left(\mathrm{Re}\right)},$$(9)

where $A_{\mathrm{ICD}}$ is the surface area of the ICD and $C_{\mathrm{ICD}}$ is the constriction factor of the device which depends of the Reynolds number (Re) (Lauritzen and Martiniussen, 2011).

All the pressure drops defined above are functions of the volume flow rates in the segments or in the devices. Considering a compositional fluid flow model and according to the principle that the total molar flow of each component is conservative, the phase volume flow rates are computed as follows:

• In the annulus segments:

$$Q_{\alpha ,k}^{a,m}=\frac{F_{\alpha ,k}^{\mathrm{mol}}}{{\rho }_{\alpha ,k}^{\mathrm{mol}}}{\displaystyle \sum _{j=k_1^m}^k}{\displaystyle \sum _{\mathrm{ic}=1}^{\mathrm{ncons}}}{q}_{\mathrm{ic},j}^{\mathrm{mol}},\hspace{1em}k=k_1^m,\dots ,k_2^m-1\hspace{1em}\mathrm{and}\hspace{1em}m=1,\cdots ,M.$$(10)

• In the ICDs:

$$Q_{\alpha ,m}^{\mathrm{ICD}}=\frac{F_{\alpha ,m}^{\mathrm{mol}}}{{\rho }_{\alpha ,m}^{\mathrm{mol}}}{\displaystyle \sum _{j=k_1^m}^{k_2^m}}{\displaystyle \sum _{\mathrm{ic}=1}^{\mathrm{ncons}}}{q}_{\mathrm{ic},j}^{\mathrm{mol}},\hspace{1em}m=1,\cdots ,M.$$(11)

• In the tubing’s segments:

$$Q_{\alpha ,k}^t=\frac{F_{\alpha ,k}^{\mathrm{mol}}}{{\rho }_{\alpha ,k}^{\mathrm{mol}}}{\displaystyle \sum _{m\mathrm{\prime }\le m-1}^{}}{\displaystyle \sum _{j=k_1^{m\mathrm{\prime }}}^{k_2^{m\mathrm{\prime }}}}{\displaystyle \sum _{\mathrm{ic}=1}^{\mathrm{ncons}}}{q}_{\mathrm{ic},j}^{\mathrm{mol}},$$(12)

• where $F_{\alpha ,k}^{\mathrm{mol}}$ and ${\rho }_{\alpha ,k}^{\mathrm{mol}}$ (resp. $F_{\alpha ,m}^{\mathrm{mol}}$ and ${\rho }_{\alpha ,m}^{\mathrm{mol}}$) are respectively the molar fraction and the molar density of the phase α in the segment k (resp. in the ICDs) and $q_{\mathrm{ic},j}^{\mathrm{mol}}$ is the total molar flow rate of the component ic between the reservoir and the well at the perforation j,

$$q_{\mathrm{ic},j}^{\mathrm{mol}}={\displaystyle \sum _{\alpha \in \left\{w,o,g\right\}}}{\mathrm{WI}}_j\left(\frac{{\mathrm{kr}}_{\alpha ,l}{\rho }_{\alpha ,l}^{\mathrm{mol}}}{{\mu }_{\alpha ,l}}\right)C_{\mathrm{ic},l}^{\alpha }\left(P_{r,l}-P_j^a\right),$$(13)

• where ${\mathrm{WI}}_j$ is the well index, ${\mathrm{kr}}_{\alpha ,l}$ and ${\mu }_{\alpha ,l}$ are the relative permeability and viscosity of the phase α and $C_{\mathrm{ic},l}^{\alpha }$ is the molar fraction of the component ic in the phase α. All these quantities are defined in the reservoir grid block l. $P_{r,l}$ is the reservoir pressure in the cell l. In the following, we denote by X _r the set of unknowns in the reservoir grid block. The phase molar fractions $F_{\alpha ,k}^{\mathrm{mol}}$are obtained from a thermodynamic flash calculation in each segment or device.

To close the system (1)–(4), we add an equation that defines the well’s control mode. A well can operate either in an imposed phase flow rate $Q_{\alpha }^{\lim }$ (or the sum of flow rates in several phases) at the surface or bottom conditions or in imposed bottom hole pressure $P^{\lim }$. Finally, we obtain the following system in the well:

$$\{\begin{array}{c}{P}_{k+1}^a-P_k^a=-{∆P}_k^a\left(Q_{w,k}^a,Q_{o,k}^a,Q_{g,k}^a\right)\hspace{1em}k=k_1^m,k_2^m-1\hspace{1em}\mathrm{and}\hspace{1em}m=1,\dots ,M\\ P_k^t-P_k^a=-{∆P}_m^{\mathrm{ICD}}\left(Q_{w,m}^{\mathrm{ICD}},Q_{o,m}^{\mathrm{ICD}},Q_{g,m}^{\mathrm{ICD}}\right)\hspace{1em}k=k_2^m\hspace{1em}\mathrm{and}\hspace{1em}m=1,\dots ,M\\ \begin{array}{c}{P}_{k+1}^t-P_k^t=-{∆P}_k^t\left(Q_{w,k}^t,Q_{o,k}^t,Q_{g,k}^t\right)\hspace{1em}k=1,\hspace{1em}N-1\\ \begin{array}{c}\mathrm{BHP}-P_N^t=-{∆P}_k^t\left(Q_{w,N}^t,Q_{o,N}^t,Q_{g,N}^t\right)\\ \{\begin{array}{c}{Q}_{\alpha }\left(X_r,\mathrm{BHP},∆P_{bh,k}^a\right)-Q_{\alpha }^{\lim }=0\\ \begin{array}{c}\mathrm{or}\\ \mathrm{BHP}-P^{\lim }=0,\end{array}\end{array}\end{array}\end{array}\end{array}$$(14)

where $Q_{\alpha }\left(X_r,\mathrm{BHP},∆P_{bh,k}^a\right)$ is the α phase volume flow rate which is computed from the sum of the molar inflow rate at the perforations $q_{\mathrm{ic},j}^{\mathrm{mol}}$.

The system (14) is formed of 2N + 1 nonlinear equations that are rewritten in the following condensed form:

$$G_w\left(P^a,P^t,\mathrm{BHP},X_r\right)=0.$$(15)

If we suppose that the reservoir’s unknowns Xr are fixed as boundary conditions of the system (15), then this nonlinear system can be solved using the Newton’s method to find the 2N+1 unknowns Pkak=1N, Pktk=1N and BHP.

• Initial guess: $P^{a,0}$, $P^{t,0}$ and, ${\mathrm{BHP}}^0$.

• While $\Vert G_w^{\lambda }\Vert >\epsilon$ do.

• Solve the linear system:

$$\left(\begin{array}{c}\begin{array}{ccc}& & \end{array}\\ \begin{array}{ccc}{\left(\frac{\partial G_w}{\partial P^a}\right)}^{\lambda }& {\left(\frac{\partial G_w}{\partial P^t}\right)}^{\lambda }& {\left(\frac{\partial G_w}{\partial \mathrm{BHP}}\right)}^{\lambda }\end{array}\\ \begin{array}{ccc}& & \end{array}\end{array}\right)\left(\begin{array}{c}\delta P^a\\ \delta P^t\\ \delta \mathrm{BHP}\end{array}\right)=-G_w\left(P^{a,\lambda },P^{t,\lambda },{\mathrm{BHP}}^{\lambda },X_r\right),$$

• Update the iterates: $P^{a,\lambda +1}=P^{a,\lambda }+\delta P^a$, $P^{t,\lambda +1}=P^{t,\lambda }+\delta P^t$ and ${\mathrm{BHP}}^{\lambda +1}={\mathrm{BHP}}^{\lambda }+\delta \mathrm{BHP}$.

• At convergence, set $P^a=P^{a,\lambda +1}$, $P^t=P^{t,\lambda +1}$ and $\mathrm{BHP}={\mathrm{BHP}}^{\lambda +1}$.

Once the pressures in the tubing and the annulus are computed, we can deduce ${∆P}_{bh,k}^a$ (reps. ${∆P}_{bh,k}^t$) the pressure drop between $P_j^a$ (resp. $P_j^t$) and BHP:

$${∆P}_{bh,k}^a=\mathrm{BHP}-P_k^a\hspace{1em}\mathrm{and}\hspace{1em}{∆P}_{bh,k}^t=\mathrm{BHP}-P_k^t\hspace{1em}k=1,\cdots ,N.$$(16)

2.3 Wells and reservoir coupling strategy

The coupling between the reservoir’s unknowns and the well’s unknowns comes from the molar flow rate of each component at the perforations:

$$q_{\mathrm{ic},k}^{\mathrm{mol}}={\displaystyle \sum _{\alpha \in \left\{w,o,g\right\}}}{\mathrm{WI}}_k\left(\frac{{\mathrm{kr}}_{\alpha ,l}{\rho }_{\alpha ,l}^{\mathrm{mol}}}{{\mu }_{\alpha ,l}}\right)C_{\mathrm{ic},l}^{\alpha }\left(P_{r,l}-P_k^a\right).$$(17)

Using the definition of ${∆P}_{bh,k}^a$ (16), we can rewrite these fluxes as:

$$q_{\mathrm{ic},k}^{\mathrm{mol}}={\displaystyle \sum _{\alpha \in \left\{w,o,g\right\}}}{\mathrm{WI}}_k\left(\frac{{\mathrm{kr}}_{\alpha ,l}{\rho }_{\alpha ,l}^{\mathrm{mol}}}{{\mu }_{\alpha ,l}}\right)C_{\mathrm{ic},l}^{\alpha }\left(P_{r,l}-\mathrm{BHP}-{∆P}_{bh,k}^a\right).$$(18)

2.3.1 Explicit coupling method

As in the standard well model, the pressure differences ${∆P}_{bh,k}^a$ are computed at the beginning of each time step and kept constant during the time step:

$$q_{\mathrm{ic},k}^{\mathrm{mol},n+1}={\displaystyle \sum _{\alpha \in \left\{w,o,g\right\}}}{\mathrm{WI}}_k{\left(\frac{{\mathrm{kr}}_{\alpha ,l}{\rho }_{\alpha ,l}^{\mathrm{mol}}}{{\mu }_{\alpha ,l}}\right)}^{n+1}{C}_{\mathrm{ic},l}^{\alpha ,n+1}\left(P_{r,l}^{n+1}-{\mathrm{BHP}}^{n+1}-{∆P}_{bh,k}^{a,n}\right),$$(19)

where n denotes the previous time step and n+1 the current time step.

${∆P}_{bh,k}^{a,n}$ are computed by solving the nonlinear system (15) with the Algorithm 1 in which the reservoir unknowns are fixed at the previous time step $X_r=X_r^n$and the initial guess for the Newton Algorithm is given by,

$$\begin{array}{c}{P}^{a,0}={\mathrm{BHP}}^n+{\Delta P}_{bh}^{a,n-1}\\ P^{t,0}={\mathrm{BHP}}^n+{\Delta P}_{bh}^{t,n-1}\\ {\mathrm{BHP}}^0={\mathrm{BHP}}^n.\end{array}$$(20)

Thus, the only well variable coupled with the reservoir is the Bottom Hole Pressure (BHP) and the global system to be solved at each time step is written as follows:

$$\{\begin{array}{c}{F}_r\left(X_r^{n+1},{\mathrm{BHP}}^{n+1}\right)=0\\ F_w\left(X_r^{n+1},{\mathrm{BHP}}^{n+1}\right)=0,\end{array}$$(21)

$F_r$ represents the set of conservation equations for the number of moles of the components in the reservoir grid blocks and $F_w$ defines the well constraint equation.

If the well is under volume flow rate control ($Q_{\alpha }^{\lim }$ with $\alpha \in \left\{w,o,g,\mathrm{wo}\right\}$), $F_w$ is,

$${F_w\left(X_r^{n+1},{\mathrm{BHP}}^{n+1}\right)=Q}_{\alpha }\left(X_r^{n+1},{\mathrm{BHP}}^{n+1},{∆P}_{\mathrm{bh}}^{a,n}\right)-Q_{\alpha }^{\lim }.$$(22)

If the well is under bottom hole pressure control $F_w$ is,

$$F_w\left(X_r^{n+1},{\mathrm{BHP}}^{n+1}\right)={\mathrm{BHP}}^{n+1}-P^{\lim }.$$(23)

The nonlinear system is solved using Newton’s method:

Explicit coupling method

Compute ${∆P}_{\mathrm{bh},k}^{a,n}$ by solving $G_w\left(P^{a,n},P^{t,n},{\mathrm{BHP}}^n,X_r^n\right)=0$ Define $q_{\mathrm{ic},k}^{\mathrm{mol},n+1}={\displaystyle \sum _{\alpha \in \left\{w,o,g\right\}}}{\mathrm{WI}}_k{\left(\frac{{\mathrm{kr}}_{\alpha ,l}{\rho }_{\alpha ,l}^{\mathrm{mol}}}{{\mu }_{\alpha ,l}}\right)}^{n+1}{C}_{\mathrm{ic},l}^{\alpha ,n+1}\left(P_{r,l}^{n+1}-{\mathrm{BHP}}^{n+1}-{∆P}_{\mathrm{bh},k}^{a,n}\right)$ Initialization: $\left(X_r^0,{\mathrm{BHP}}^0\right)=\left(X_r^n,{\mathrm{BHP}}^n\right)$ While $\Vert F_r^{\lambda }\Vert >{\epsilon }_r$ and $\Vert F_w^{\lambda }\Vert >{\epsilon }_w$ do  Solve the linear system: $$\left(\begin{array}{cc}\frac{\mathrm{\partial }{F}_r^{\lambda }}{\mathrm{\partial }{X}_r}& \frac{\mathrm{\partial }{F}_r^{\lambda }}{\mathrm{\partial }\mathrm{BHP}}\\ \frac{\mathrm{\partial }{F}_w^{\lambda }}{\mathrm{\partial }{X}_r}& \frac{\mathrm{\partial }{F}_w^{\lambda }}{\mathrm{\partial }\mathrm{BHP}}\end{array}\right)\left(\begin{array}{c}\\ \begin{array}{c}\begin{array}{c}\delta X_r\\ \end{array}\\ \delta \mathrm{BHP}\end{array}\\ \end{array}\right)=-\left(\begin{array}{c}\\ \begin{array}{c}{F}_r\left(X_r^{\lambda },{\mathrm{BHP}}^{\lambda }\right)\\ \begin{array}{c}\\ F_w\left(X_r^{\lambda },{\mathrm{BHP}}^{\lambda }\right)\end{array}\end{array}\\ \end{array}\right)$$(24)  Update the iterates: $X_r^{\lambda +1}=X_r^{\lambda }+\delta X_r$ and ${\mathrm{BHP}}^{\lambda +1}={\mathrm{BHP}}^{\lambda }+\delta \mathrm{BHP}$. At convergence, set $X_r^{n+1}=X_r^{\lambda +1}$ and ${\mathrm{BHP}}^{n+1}={\mathrm{BHP}}^{\lambda +1}$.

This method is efficient and works well as long as the well is operating at a prescribed flow rate. However, as soon as the well switches to fixed bottom hole pressure constraint, we observe numerical oscillations in the production curves due to the crossflow phenomena.

To deal with these instabilities, the pressure drop in equation (19) must be treated implicitly

$$q_{\mathrm{ic},k}^{\mathrm{mol},n+1}={\displaystyle \sum _{\alpha \in \left\{w,o,g\right\}}}{\mathrm{WI}}_k{\left(\frac{{\mathrm{kr}}_{\alpha ,l}{\rho }_{\alpha ,l}^{\mathrm{mol}}}{{\mu }_{\alpha ,l}}\right)}^{n+1}{C}_{\mathrm{ic},l}^{\alpha ,n+1}\left(P_{r,l}^{n+1}-{\mathrm{BHP}}^{n+1}-{∆P}_{\mathrm{bh},k}^{a,n+1}\right).$$(25)

The main drawback of this approach is that all the well’s unknowns and the perforated grid block’s unknowns are coupled and the equations of the system (15) and of the system (21) must be solved simultaneously.

2.3.2 Semi-implicit coupling method

In order to keep the effectiveness of the previous method, we propose the semi-implicit coupling method which consists in calculating the pressure losses ${∆P}_{\mathrm{bh},k}^{a,n+1}$ at the beginning of each global Newton’s iteration and to keep them constant during the iteration. These pressure drops in the well are then calculated by solving the non-linear system (15) using the latest iterative solution in the reservoir grid-block as boundary conditions.

The semi-implicit algorithm can be written as follows:

Semi-implicit coupling method

Compute ${∆P}_{\mathrm{bh},k}^{a,n}$ by solving $G_w\left(P^{a,n},P^{t,n},{\mathrm{BHP}}^n,X_r^n\right)=0$ Set $q_{\mathrm{ic},k}^{\mathrm{mol},n+1}={\displaystyle \sum _{\alpha \in \left\{w,o,g\right\}}}{\mathrm{WI}}_k{\left(\frac{{\mathrm{kr}}_{\alpha ,l}{\rho }_{\alpha ,l}^{\mathrm{mol}}}{{\mu }_{\alpha ,l}}\right)}^{n+1}{C}_{\mathrm{ic},l}^{\alpha ,n+1}\left(P_{r,l}^{n+1}-{\mathrm{BHP}}^{n+1}-{∆P}_{\mathrm{bh},k}^{a,n}\right)$ Initialization: $\left(X_r^0,{\mathrm{BHP}}^0\right)=\left(X_r^n,{\mathrm{BHP}}^n\right)$ While $\Vert F_r^{\lambda }\Vert >{\epsilon }_r$ and $\Vert F_w^{\lambda }\Vert >{\epsilon }_w$ do  If ${\mathrm{BHP}}^n=P^{\lim }$ then   Compute ${∆P}_{\mathrm{bh},k}^{a,\lambda }$ by solving $G_w\left(P^{a,\lambda },P^{t,\lambda },{\mathrm{BHP}}^{\lambda },X_r^{\lambda }\right)=0$   Set $q_{\mathrm{ic},k}^{\mathrm{mol},\lambda +1}={\displaystyle \sum _{\alpha \in \left\{w,o,g\right\}}}{\mathrm{WI}}_k{\left(\frac{{\mathrm{kr}}_{\alpha ,l}{\rho }_{\alpha ,l}^{\mathrm{mol}}}{{\mu }_{\alpha ,l}}\right)}^{\lambda +1}{C}_{\mathrm{ic},l}^{\alpha ,\lambda +1}\left(P_{r,l}^{\lambda +1}-{\mathrm{BHP}}^{\lambda +1}-{∆P}_{\mathrm{bh},k}^{a,\lambda }\right)$  Solve the linear system: $$\left(\begin{array}{cc}\frac{\mathrm{\partial }{F}_r^{\lambda }}{\mathrm{\partial }{X}_r}& \frac{\mathrm{\partial }{F}_r^{\lambda }}{\mathrm{\partial }\mathrm{BHP}}\\ \frac{\mathrm{\partial }{F}_w^{\lambda }}{\mathrm{\partial }{X}_r}& \frac{\mathrm{\partial }{F}_w^{\lambda }}{\mathrm{\partial }\mathrm{BHP}}\end{array}\right)\left(\begin{array}{c}\delta X_r\\ \delta \mathrm{BHP}\end{array}\right)=-\left(\begin{array}{c}{F}_r\left(X_r^{\lambda },{\mathrm{BHP}}^{\lambda }\right)\\ F_w\left(X_r^{\lambda },{\mathrm{BHP}}^{\lambda }\right)\end{array}\right).$$(26)  Update the iterates: $X_r^{\lambda +1}=X_r^{\lambda }+ \delta X_r$ and ${\mathrm{BHP}}^{\lambda +1}={\mathrm{BHP}}^{\lambda }+\delta \mathrm{BHP}$. At convergence, set $X_r^{n+1}=X_r^{\lambda +1}$ and ${\mathrm{BHP}}^{n+1}={\mathrm{BHP}}^{\lambda +1}$.

The advantage of this method, in addition to being more stable, is that we do not have to change the global linear solver because the systems (24) and (26) have exactly the same structure. This method is not intrusive and therefore easy to implement.

3 Numerical examples

3.1 Example 1 – a synthetic case.

Consider a sugar box reservoir of size 4100 × 4100 × 700 ft with three different permeability regions: 50 mD, 20 mD, and 800 mD, as shown in Figure 2. A horizontal well of length 2000 ft traverses these three zones, and the perforated drain is subdivided into four compartments equipped with nozzle type ICDs. Each compartment has a length of 500 ft. The first compartment is in the zone of 50 mD permeability, the second and the third compartments are in the zone of 20 mD permeability and the fourth compartment is in the zone of 800 mD permeability. The top of the reservoir is 9000 ft, and the horizontal drain is located at a depth of 9315 ft. An aquifer is found below the horizontal well, at the depth of 9350 ft.

In the simulation, this reservoir is uniformly discretized by 41 grid cells in the x-direction, 41 grid cells in the y-direction and 10 cells in the z-direction. The cell sizes are 100 × 100 × 70 ft. The horizontal drain contains 20 active well cells, and each compartment connects to five well cells. The initial reservoir pressure is 3700 psi. The control mode adopted in this study is a maximum total liquid rate of 2000 bbl/day and a minimum bottom-hole pressure of 2900 psi.

In this synthetic case, the number and the size of the nozzles are designed from a single-phase flow problem. Figure 3 shows the flow rate of these four compartments with optimized ICDs in a single-phase flow simulation. At the beginning, Compartment 4 has a very high flow rate due to the perforation in the high permeability zone, and the other compartments show low flow rates. Due to ICD controls, the flow rate in Compartment 4 quickly drops to 500 bbl/day, and the rates in other compartments increase to around 500 bbl/day. The productions in the four compartments are quite close. These ICDs are then used to oil/water two-phase flow simulations with an aquifer below.

Figure 4 presents the oil and water production of the horizontal well with and without ICD equipment. It is clear that with ICDs the oil rate is significantly increased while the water production is reduced. Figure 5 presents the oil and water rates in each compartment for the base case without ICDs, and Figure 6 shows the results for the case equipped with ICDs. For the open-hole horizontal well without ICD equipment, very high flow rates are observed in the high permeability zone (Compartment 4) for both oil and water productions. The production in this zone can be greatly reduced with ICD completions. As a consequence, the water production is globally reduced. The advantage of ICD completions is simulated using the proposed semi-implicit model.

We now consider the case where the reservoir permeability is generated using a geostatistical model as shown by Figure 7. The vertical to horizontal permeability ratio is 0.1. The top six layers is the oil zone, and the bottom four layers is the aquifer. One producer (PRO) with four branches is drilled in Layer 3, while the injector (WINJ) having two branches is drilled with one branch at Layer 3 in the oil zone and the other branch at Layer 9 in the aquifer (Fig. 7). A total flow rate of 12 000 bbl/day is imposed for both the injector and the producer. The initial reservoir pressure is 3700 psi, and the limited bottom-hole pressure is 3000 psi for the producer and 6800 psi for the injector.

As the upper branch (Branch 1) of the injector is located in the same layer as the production well, early water breakthrough may occur due to the water flooding from the upper injection branch. To get a better injection repartition, one solution is to equip the injection well with ICDs. Another solution is to install ICV in the top branch (Branch 1) to reduce the injection volume in that branch. For the ICD completion, each branch is completed with four compartments equipped with nozzle type ICDs. The nozzle sizes are selected according to the reservoir heterogeneity, petrophysical and fluid properties, as well as the sweep efficiency. Small nozzles are installed in Branch 1 and in the last compartment of Branch 2.

Figure 8 compares the oil and water productions simulated with open-hole completion, ICD completion and ICV completion in the injection well. Equipped with intelligent completions in the injection well can control the sweep efficiency and hence delay the water breakthrough. Figure 9 presents the water saturation maps for different scenarios in Layer 3 after 5 years’ water injection, and Figure 10 shows the water saturation in a cross-section at the same time. Figure 11 presents the water repartition in the injection well. For the open-hole completion, the water injection is almost equally repartitioned in these two branches, while using ICD or ICV the volume of injected water in the upper branch (Br1 in Fig. 11) is greatly reduced, and the reservoir pressure support is maintained through the lower branch injection. The water sweeping front in the oil zone is amortized, and therefore the water production is delayed.

It has to be mentioned that for these well configurations, it is difficult to improve the oil production by other strategies. For example, implementation of ICD or ICV in the four branched production well cannot help to reduce the undesired water production.

Figure 12 shows the well bottom hole pressure for different completions. At the beginning, the bottom-hole pressures do not reach their limit conditions, and the total flow rate is imposed at the wells. The production well flowing pressure is constrained by the limit pressure at around 7–9 years, and the injection well pressure is constrained by the upper pressure limit at around 8 years with the ICD completion.

For these two cases (the case of three permeability zones and the heterogeneous case), each model is repeatedly simulated 10 times on a linux machine to evaluate the CPU time. Table 1 shows the mean CPU time and its variance. The semi-implicit approach for ICD and ICV modeling does not much impact on the CPU time in this example. With the ICD completion, CPU time is almost unchanged comparing to the open-hole cases, while the ICV completion increases very slightly the CPU time with an average around 1.3%.

3.2 Example 2. Case SPE 10

SPE 10 is a highly heterogeneous case initially developed for upscaling comparison (Christie and Blunt, 2001). We use this fine scale heterogeneity field to simulate ICD completions with the proposed semi-implicit approach. The fine grid contains 60 blocks in the x direction, 220 blocks in the y direction and 85 layers. The reservoir contains oil and immobile water in the initial state. One fully penetrated injector is located at the center of this field and four fully penetrated vertical producers are drilled near the four corners (Fig. 13). A water rate of 5000 bbl/day is imposed at the injector, and the four production wells are constrained by a pressure limit of 4000 psi.

As the reservoir is very heterogeneous, ICDs are installed in the production wells to control the water breakthrough coming from connected channel layers. Each production well is completed with five compartments and equipped with nozzle type ICDs.

Figure 14 compares the oil and water productions at the four producers for the case with and without ICD completions, and Figure 15 shows the field oil and water productions. It is clear that the field oil production is increased, while the whole water production is almost unchanged or even slightly reduced with ICD installations. The water sweeping is slowed down in some permeable layers due to the implementation of ICDs, and its sweep efficiencies are increased in other layers. The equipment of ICD shows its great efficiency in Wells P1 and P4 as well as P3, where the water productions are reduced and oil productions increased. Wells P2 is also a successful candidate. Although the water production is increased by 69.7 Mbbl in P2, a higher volume of oil of 85.1 Mbbl is produced. With ICD installations, a higher oil production is found in all the wells with the water-cuts controlled in a lower level (Fig. 16).

The simulation without ICD equipment needs 2037 Newton iterations with 15 637 s in CPU time, while the simulation with ICDs requires 2044 Newton iterations and 15 737 s in CPU time. The total CPU time is increased only around 0.6% in this very heterogeneous field. This augmentation in CPU time is not significant. The simulation time and the number of Newton iterations with ICD completions depend on the ICD design. If extremely small nozzle sizes are used, a very high pressure drop may be created across the ICD, and this may impact on the convergence issue and CPU time. For the commonly-used ICD design, the proposed semi-implicit approach for the reservoir and well model coupling can simulate the flow through the intelligent completions with few impacts on the CPU time.

4 Conclusion

A semi-implicit approach is presented to couple a standard reservoir simulator and a well flow model with ICD/ICV completions. Further details have been given for this semi-implicit coupling. The presence of ICD/ICV equipment, which generates generally very high pressure drop across the device, can penalize CPU time in the coupled modeling. The proposed semi-implicit approach provides an easy way for the coupling implementation with very limited impact on computational time. Numerical examples show that the proposed model can simulate conveniently ICD/ICV completions with almost the same CPU time as for the case of open-hole completion.

References

All Tables

All Figures