Regular Article
A semiimplicit approach for the modeling of wells with inflow control completions
IFP Energies nouvelles, 1 et 4 avenue de BoisPréau, 92852 RueilMalmaison Cedex, France
^{*} Corresponding author: didieryu.ding@ifpen.fr
Received:
20
January
2020
Accepted:
29
April
2020
In the last two decades, new technologies have been introduced to equip wells with intelligent completions such as Inflow Control Device (ICD) or Inflow Control Valve (ICV) in order to optimize the oil recovery by reducing the undesirable production of gas and water. To optimally define the locations of the packers and the characteristics of the valves, efficient reservoir simulation models are required. This paper is aimed at presenting the specific developments introduced in a multipurpose industrial reservoir simulator to simulate such wells equipped with intelligent completions taking into account the pressure drop and multiphase flow. An explicit coupling or decoupling of a reservoir model and a well flow model with intelligent completion makes usually unstable and nonconvergent results, and a fully implicit coupling is CPU time consuming and difficult to be implemented. This paper presents therefore a semiimplicit approach, which links on one side to the reservoir simulation model and on the other side to the well flow model, to integrate ICD and ICV.
© E. Flauraud & D.Y. Ding, published by IFP Energies nouvelles, 2020
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
Nomenclature
${A}_{\mathrm{ICD}}$ : Surface area of ICD
${A}_{\mathrm{ICV}}$ : Surface area of ICV
${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
${\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}^{\mathrm{lim}}$ : Limit bottom hole pressure
${P}_{r,l}$ : Pressure in the reservoir grid block l
${\u2206P}_{k}^{a}$ : Pressure drop in the segment k of the annulus
${\u2206P}_{k}^{t}$ : Pressure drop in the segment k of the tubing
${\u2206P}_{k}^{t,\mathrm{hydro}}$ : Hydrostatic pressure drop in the segment k of the tubing
$\u2206{P}_{k}^{t,\mathrm{fric}}$ : Friction pressure drop in the segment k of the tubing
$\u2206{P}_{k}^{t,\mathrm{ICV}}$ : Pressure drop in IVC of the segment k of the tubing
$\u2206{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}^{\mathrm{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
${\u2206Z}_{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 : Total number of compartments or ICDs
${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
ncons : Total number of components
n : Refer to the discret time t^{n}
n + 1: Refer to the discret time t ^{n+1}
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., AlEnezi et al., 2012; AlKhelaiwi and Davies 2007; AlThuwaini et al., 2009; BroniBediako et al., 2019; Das and AlEnezi, 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., AlKhelaiwi, 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 prepacked 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 semianalytical type approach is used to describe the pressure drop through ICDs for various design configurations (see, e.g., Das and AlEnezi, 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 semianalytical approach (Das and AlEnezi, 2014; Das et al., 2012). This kind of approaches can be considered suitable for a rapid design, which focuses only on the immediate nearwell region and needs a fast response. However, to evaluate longterm 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 longterm 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 multisegment 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 multisegment 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 multisegment approach allows the well geometry to be represented more accurately and to better describe fluid flows in the wellbore. However, solving implicitly a multisegment 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 semiimplicit approach to integrate ICD/ICV model in a reservoir simulator with a simplified multisegment 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 semiimplicit 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 AlEnezi, 2014).
Fig. 1 Example of multilateral well discretization. 
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,N1$, 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 undercritical 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 steadystate momentum equations:

In the annulus, we consider only hydrostatic pressure loss:
$${P}_{k+1}^{a}{P}_{k}^{a}={\u2206P}_{k}^{a}\left({Q}_{w,k}^{a},{Q}_{o,k}^{a},{Q}_{g,k}^{a}\right)=\u2206{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, $\u2206{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}={\u2206P}_{k}^{t}\left({Q}_{w,k}^{t},{Q}_{o,k}^{t},{Q}_{g,k}^{t}\right)=\left(\u2206{P}_{k}^{t,\mathrm{hydro}}+\u2206{P}_{k}^{t,\mathrm{fric}}+\u2206{P}_{k}^{t,\mathrm{ICV}}\right),\hspace{1em}k=1,\cdots ,N,$$(2)
$$\mathrm{BHP}{P}_{N}^{t}={\u2206P}_{N}^{t}\left({Q}_{w,N}^{t},{Q}_{o,N}^{t},{Q}_{g,N}^{t}\right)=\left(\u2206{P}_{N}^{t,\mathrm{hydro}}+\u2206{P}_{N}^{t,\mathrm{fric}}+\u2206{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. $\u2206{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}={\u2206P}_{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:

where g is the gravity constant, ${\u2206Z}_{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}$.

$$\u2206{P}_{k}^{\mathrm{fric}}=\frac{8{l}_{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.

$$\u2206{P}_{k}^{\mathrm{ICV}}=\frac{{\overline{\rho}}_{k}^{\mathrm{mass}}{Q}_{T,k}^{2}}{2{A}_{\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.

$$\u2206{P}_{m}^{\mathrm{ICD}}=\frac{{\overline{\rho}}_{m}^{\mathrm{mass}}{Q}_{T,m}^{2}}{2{A}_{\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 m1}^{}}{\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}^{\mathrm{lim}}$ (or the sum of flow rates in several phases) at the surface or bottom conditions or in imposed bottom hole pressure ${P}^{\mathrm{lim}}$. Finally, we obtain the following system in the well:
$$\{\begin{array}{c}{P}_{k+1}^{a}{P}_{k}^{a}={\u2206P}_{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}={\u2206P}_{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}={\u2206P}_{k}^{t}\left({Q}_{w,k}^{t},{Q}_{o,k}^{t},{Q}_{g,k}^{t}\right)\hspace{1em}k=1,\hspace{1em}N1\\ \begin{array}{c}\mathrm{BHP}{P}_{N}^{t}={\u2206P}_{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},\u2206{P}_{bh,k}^{a}\right){Q}_{\alpha}^{\mathrm{lim}}=0\\ \begin{array}{c}\mathrm{or}\\ \mathrm{BHP}{P}^{\mathrm{lim}}=0,\end{array}\end{array}\end{array}\end{array}\end{array}$$(14)
where ${Q}_{\alpha}\left({X}_{r},\mathrm{BHP},\u2206{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 ${\u2206P}_{bh,k}^{a}$ (reps. ${\u2206P}_{bh,k}^{t}$) the pressure drop between ${P}_{j}^{a}$ (resp. ${P}_{j}^{t}$) and BHP:
$${\u2206P}_{bh,k}^{a}=\mathrm{BHP}{P}_{k}^{a}\hspace{1em}\mathrm{and}\hspace{1em}{\u2206P}_{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 ${\u2206P}_{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}{\u2206P}_{bh,k}^{a}\right).$$(18)
2.3.1 Explicit coupling method
As in the standard well model, the pressure differences ${\u2206P}_{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}{\u2206P}_{bh,k}^{a,n}\right),$$(19)
where n denotes the previous time step and n+1 the current time step.
${\u2206P}_{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,n1}\\ {P}^{t,0}={\mathrm{BHP}}^{n}+{\Delta P}_{bh}^{t,n1}\\ {\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}^{\mathrm{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},{\u2206P}_{\mathrm{bh}}^{a,n}\right){Q}_{\alpha}^{\mathrm{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}^{\mathrm{lim}}.$$(23)
The nonlinear system is solved using Newton’s method:
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}{\u2206P}_{\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 Semiimplicit coupling method
In order to keep the effectiveness of the previous method, we propose the semiimplicit coupling method which consists in calculating the pressure losses ${\u2206P}_{\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 nonlinear system (15) using the latest iterative solution in the reservoir gridblock as boundary conditions.
The semiimplicit algorithm can be written as follows:
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.
Fig. 2 Synthetic case with three permeability zones. 
In the simulation, this reservoir is uniformly discretized by 41 grid cells in the xdirection, 41 grid cells in the ydirection and 10 cells in the zdirection. 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 bottomhole pressure of 2900 psi.
In this synthetic case, the number and the size of the nozzles are designed from a singlephase flow problem. Figure 3 shows the flow rate of these four compartments with optimized ICDs in a singlephase 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 twophase flow simulations with an aquifer below.
Fig. 3 Oil rate in the four compartments. 
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 openhole 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 semiimplicit model.
Fig. 4 Oil and water production rate with and without ICDs. 
Fig. 5 Compartment oil (left) and water (right) productions without ICDs. 
Fig. 6 Compartment oil (left) and water (right) productions with ICDs. 
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 bottomhole pressure is 3000 psi for the producer and 6800 psi for the injector.
Fig. 7 Permeability distribution in Layer 3 and in a yz cross section. 
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 openhole 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 crosssection at the same time. Figure 11 presents the water repartition in the injection well. For the openhole 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.
Fig. 8 Oil and water rates on the production well. 
Fig. 9 Water saturation map in Layer 3 after 5 years water injection. 
Fig. 10 Water saturation map in a yz crosssection. 
Fig. 11 Water repartition in the injection well. 
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 bottomhole 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.
Fig. 12 Well bottom hole pressure. 
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 semiimplicit 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 openhole cases, while the ICV completion increases very slightly the CPU time with an average around 1.3%.
CPU time comparison (second).
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 semiimplicit 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.
Fig. 13 Reservoir permeability and well locations. 
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 watercuts controlled in a lower level (Fig. 16).
Fig. 14 Comparison of oil and water production at the four production wells. 
Fig. 15 Oil and water productions in the whole reservoir. 
Fig. 16 Water cut comparison. 
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 commonlyused ICD design, the proposed semiimplicit 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 semiimplicit 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 semiimplicit 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 semiimplicit 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 openhole completion.
References
 AlEnezi K., Das O.P., Aslam M., Ziyad K., Fipke S.R. (2012) Successful case histories of smart multilateral well with inflow control device and inflow control valve for lifecycle proactive reservoir management in high mobility reservoir, Minagish Field West Kuwait, in: SPE 161632, Abu Dhabi International Petroleum Conference and Exhibition, 11–14 November. [Google Scholar]
 AlKhelaiwi F.T.M. (2013) A comprehensive approach to the design of advanced well completions, PhD Dissertation, HeriotWatt University, Edinburgh. [Google Scholar]
 AlKhelaiwi F.T., Davies D.R. (2007) Inflow control devices: Application and value quantification of a developing technology, in: Paper SPE 108700 Presented at the 2007 International Oil Conference and Exhibition, Veracruz, Mexico, 27–30 June. [Google Scholar]
 AlThuwaini J., Shenawi S., Shenawi S., Yuen B. (2009) Simulation and optimisation of complex architecture wells with smart completions, in: Paper OTC 20131 Presented at the 2009 Offshore Technology Conference held in Houston, Texas, USA, 4–7 May. [Google Scholar]
 Brekke K., Lien S. (1992) New and simple completion methods for horizontal wells improved production performance in highpermeability thin oil zones, in: SPE 24762 Presented at the 67th SPE Annual Technical Conference and Exhibition, Washington, DC, USA, 4–7 October. [Google Scholar]
 BroniBediako E., Issaka Fuseini N., Ayitey Akoto R.N., Thompson Brantson E. (2019) Application of intelligent well completion in optimising oil production from oil rim reservoirs, Adv. GeoEnergy Res. 3, 4, 343–354. doi: 10.26804/ager.2019.04.01. [CrossRef] [Google Scholar]
 Christie M.A., Blunt M.J. (2001) Tenth SPE comparative solution project: A comparison of upscaling techniques, SPE Reserv. Evalu. Eng., pp. 308–317. [Google Scholar]
 Das O.P., AlEnezi K. (2014) A novel workflow for intelligent well inflow control valve design by integrating reservoir dynamics to facilitate proactive reservoir management in Minagish Field, West Kuwait, in: SPE 170803 Presented at the SPE Annual Technical Conference and Exhibition Held in Amsterdam, The Netherlands, 27–29 October. [Google Scholar]
 Das O.P., AlEnezi K., Aslam M., ElGezeeri T., Ziyab K., Fipke S.R., Ewens S. (2012) Novel design and implementation of Kuwait’s first smart multilateral well with inflow control device and inflow control valve for lifecycle reservoir management in high mobility reservoir, West Kuwait, in: SPE 159261 Presented at SPE Annual Technical Conference and Exhibition, 8–10 October, San Antonio, USA. [Google Scholar]
 Dimitios K., Hembling D., AlDawood N., AlQatari S., Simonian S., Salerno G. (2009) Optimizing horizontal well performance in nonuniform pressure environments using passive inflow control devices, in: Paper OTC 20129 Presented at the 2009 Offshore Technology Conference held in Houston, Texas, USA, 4–7 May. [Google Scholar]
 Gurses S., Chochua G., Rudic A., Kumar A. (2019) Dynamic Modeling and Design Optimization of Cyclonic Autonomous Inflow Control Devices, in: SPE 193824 presented at the SPE Reservoir Simulation Conference, Galveston, TX, USA, 1011 April. doi: 10.2118/193824MS [Google Scholar]
 HoJeen S., Dogru A. (2007) Modeling of equalizer production system and smart well applications in fullfield studies, in: Paper SPE 111288 Presented at the 2007 SPE/EAGE Reservoir Characterization and Simulation Conference, Abu Dhabi, UAE, 28–31 October. [Google Scholar]
 Holmes J.A., Barkve T., Lund O. (1998) Application of a multisegment well model to simulate flow in advanced wells, in: SPE 50646 Presented at the European Petroleum Conference, 20–22 October, The Hague, Netherlands. [Google Scholar]
 Holmes J.A., Byer T.J., Edwards D.A., Stone T.W., Pallister I., Shaw G.J., Walsh D. (2010) A unified wellbore model for reservoir simulation, in: SPE 134928 Presented at the SPE Annual Technical Conference and Exhibition, 19–22 September, Florence, Italy. [Google Scholar]
 Javid K., Mustafa M., Chitre S., Anurag A., Sayed M., Kuliyev M., Mishra A., Khalili Al Hosany K., Saeed Y. (2018) Comprehensive ICD/ICV completion design workflow practiced in Green oilfield offshore, Abu Dhabi, in: Paper SPE 192645 Presented at the Abu Dhabi International Petroleum Exhibition & Conference, Abu Dhabi, UAE, 12–15 November. [Google Scholar]
 Lauritzen J.E., Martiniussen I.B. (2011) Single and multiphase flow loop testing results for industry standard inflow control devices, in: SPE 146347 Presented at the Offshore Europe, 6–8 September, Aberdeen, UK. [Google Scholar]
 Raffn A.G., Hundsnes S., Kvernstuen S., Moen T. (2007) ICD screen technology used to optimize waterflooding in injector well, in: SPE 106018, Presented at the SPE Production and Operations Symposium, Oklahoma City, Oklahoma, USA, 31 March–3 April. [Google Scholar]
 Stone T., Moen T., Edwards D., Shadchnev A., Rashid K., Kvilaas G., Christoffersen K. (2015) Optimized design of autonomous inflow control devices for gas and water coning, in: SPE 173203 Presented at the SPE Reservoir Simulation Symposium, Houston, 23–25 February. [Google Scholar]
 Wan J., Dale B.A., Ellison T.K., Benish T.G., Grubert M. (2008) Coupled well and reservoir simulation models to optimize completion design and operations for subsurface control, in: Paper SPE 113635 Presented at the 2008 Europec/EAGE Annual Technical Conference and Exhibition, Rome, Italy, 9–12 June. [Google Scholar]
 Youngs B., Neylon K., Holmes J. (2009) Recent advances in modeling well inflow control devices in reservoir simulation, in: IPCT 13925, Doha, Qatar, 7–9 December. [Google Scholar]
All Tables
All Figures
Fig. 1 Example of multilateral well discretization. 

In the text 
Fig. 2 Synthetic case with three permeability zones. 

In the text 
Fig. 3 Oil rate in the four compartments. 

In the text 
Fig. 4 Oil and water production rate with and without ICDs. 

In the text 
Fig. 5 Compartment oil (left) and water (right) productions without ICDs. 

In the text 
Fig. 6 Compartment oil (left) and water (right) productions with ICDs. 

In the text 
Fig. 7 Permeability distribution in Layer 3 and in a yz cross section. 

In the text 
Fig. 8 Oil and water rates on the production well. 

In the text 
Fig. 9 Water saturation map in Layer 3 after 5 years water injection. 

In the text 
Fig. 10 Water saturation map in a yz crosssection. 

In the text 
Fig. 11 Water repartition in the injection well. 

In the text 
Fig. 12 Well bottom hole pressure. 

In the text 
Fig. 13 Reservoir permeability and well locations. 

In the text 
Fig. 14 Comparison of oil and water production at the four production wells. 

In the text 
Fig. 15 Oil and water productions in the whole reservoir. 

In the text 
Fig. 16 Water cut comparison. 

In the text 