Fluid flow in a reservoir drained by a multiple fractured horizontal well

S.V. Golovin¹¹1Corresponding author. E-mail: golovin@hydro.nsc.ru, Address: pr. Lavrentyeva 15, Novosibirsk, 630090, Russia K.A. Gadylshina Lavrentyrev Institute of Hydrodynamics, Novosibirsk, Russia Novosibirsk State University, Novosibirsk, Russia

Abstract

A mathematical model for computation of the fluid pressure in a reservoir drained by a horizontal multiple fractured well is proposed. The model is applicable for an arbitrary network of fractures with different finite conductivities of each segment, for variable in space and time physical parameters of the reservoir and for different field development plans. The variational formulation of the model allows effective numerical simulation using the finite element method. Case studies demonstrate how the main flow characteristics (well productivity, pressure distribution) depend on the geometrical and physical characteristics of the reservoir and of the fracture network. The presented model is suitable for estimation of the productivity of a multiple fractured well and as an optimization tool for efficient reservoir development.

keywords:

multiple fracturing, horizontal well , productivity optimizaiton, numerical modeling , finite element method

1 Introduction

Technologies of horizontal drilling and multi-stage fracturing play a key role in a field development for low-permeable reservoirs. Modern advances in engineering allows creation of exact design protocols that consider peculiar properties of the reservoir under development. The appropriate mathematical modelling of future functioning of the multiple fractured horizontal well, estimation of fluid inflow and prediction of the dynamics of the depression zone are important for the rational planning of the field development. The model should take into account geometrical characteristics of the reservoir, variable permeability of the rock, finite hydraulic resistance of fractures and the wellbore, and also be capable to compute the result within a reasonable time on a PC. The prognosis of the inflow and of the geometry of the depression zone can be used as a quick estimation of the quality of the planned wellbore stimulation, or for determining initial data for more advanced industry reservoir simulators.

There are a number of analytical solutions for estimation of productivity of a multifractured horizontal well [1, 2, 3, 4, 5]. Although requiring very limited computational resources, these solutions are obtained under restrictive assumptions and simplifications, which limits the area of their application. Hybrid [6, 7, 8, 9, 10, 11, 12] and numerical models [17, 18, 19, 20, 21, 22] are mostly use not obvious hypotheses regarding the character of the flow near the fractures and the wellbore. Besides, none of the cited works provide a complete time-dependent pressure field over the reservoir, which might be useful as initial data in commercial simulators.

In present work we develop the hydraulic model proposed in [13, 14] by presenting more efficient algorithm of conjugation of flows in the domains of different dimensions and by taking into account variability of physical properties of the reservoir and fractures. The model describes a 3D filtration Darcy flow in the reservoir conjugated with the 2D flow in fractures and with the 1D flow in the wellbore. The conjugation is achieved by a proper choice of conditions over common boundaries of the domains of different dimensions. Our approach is similar to the one presented in [15] for description of fractures as interfaces in porous medium. For efficient numerical computations in case of a simplified fractures geometry where all fractures extend from the bottom to the top of the reservoir, we derive a simplified 2D model by averaging pressure along the vertical coordinate. The 2D model has an advantage of faster computation although still reflecting the main geometrical properties of the reservoir.

The model is reformulated in a weak form such that equations for all components of the fluid flow are incorporated into one weak problem suitable for numerical solution by the finite element method. The advantage of the proposed algorithm is that it computes the flow in all segments simultaneously in one time step by the implicit numerically stable scheme. The algorithm also allows taking into account variability of physical properties of the reservoir, finite conductivity of fractures and wellbore. The case studies given in the last section of the paper reveals the dependence of the productivity of the fractured wellbore on the geometrical and physical characteristics of the reservoir and of the fracture network. This analysis demonstrate that the model can be used as an optimization tool for the proper planning of the reservoir development.

Nomenclature
$p$ , $p_{0}$ , $p_{w}$ , $p_{\infty}$	pressure, MPa	$\rho$	density, g/cm³
$m$ , $m_{0}$	porosity, $\sim$	$k$	permeability, D
$\mu$	dynamic viscosity, Pa $\cdot$ s	$t$	time, s
$L=L_{x}$ , $L_{y}$ , $H=L_{z}$	sizes of the reservoir, m	$\Omega$	the reservoir
$h_{j}$ , $j=1,\ldots,N$	fracture heights, cm	$\Sigma$	outer boundary of $\Omega$
$d_{j}$ , $j=1,\ldots,N$	fracture aperture, cm	$\Gamma=\bigcup\limits_{i=0}^{N}\Gamma_{i}$	inner boundary of $\Omega$
$k_{j}$ , $j=1,\ldots,N$	fracture permeability, D	$\theta_{j}$	flow rate coefficient, m³s^-1MPa^-1
$\gamma_{j}$ , $j=0,\ldots,N$	imperfection coefficient, $\sim$	$R$	wellbore radius, m
$\varepsilon$	the elastic capacity coefficient, MPa^-1	$L_{w}$	wellbore length, m
$\nu$	outer normal vector	$\psi$	test function
$Q_{0}$	flow rate at wellbore, g/s	$U$	pressure within the wellbore, MPa
$\mathbf{q}_{j}$ , $j=1,\ldots,N$	total discharge of fluid in fractures, cm²/s	$s$	local coordinate along wellbore
$q_{j}$ , $j=0,\ldots,N$	velocity of fluid inflow to the wellbore ( $j=0$ ) and fractures ( $j>0$ ), cm/s	$\alpha_{0}$	proportionality coefficient, m²s^-1MPa^-1

2 Formulation of the mathematical model

We observe a single-fluid model of fluid filtration in a rectangular reservoir $\Omega$ with dimensions $L_{x}\times L_{y}\times L_{z}$ (see Figure 1). For brevity, we also use the notations $L=L_{x}$ and $H=L_{z}$ for the reservoir’s horizontal and vertical sizes. The reservoir is exploited by a horizontal well with multiple hydraulic fractures placed arbitrarily along the wellbore. The equation of the one-phase filtration follows from the mass continuity equation and the Darcy law in the form

\frac{\partial(m\rho)}{\partial t}-\mathrm{div\,}\left(\rho\frac{k}{\mu}\nabla p\right)=0,\quad\mathbf{x}\in\Omega\setminus\Gamma.

(1)

Here $\Omega\subset\mathbb{R}^{3}$ is the modelling domain (the reservoir), $p$ is the pressure, $\rho$ is the fluid density, $m$ is the porosity, $k$ is the given permeability of the reservoir, $\mu$ is the fluid viscosity. We denote by $\Sigma$ is the outer boundary of $\Omega$ , and by $\Gamma=\bigcup\limits_{i=0}^{N}\Gamma_{i}$ the inner boundaries of $\Omega$ , comprising the wellbore $(i=0)$ and fractures $(1\leqslant i\leqslant N)$ . In present work we neglect the gravity force.

Refer to caption — Figure 1: Schematics of the rectangular reservoir $\Omega$ with the outer boundary $\Sigma$ drained by a horizontal wellbore $\Gamma_{0}$ with multiple fractures $\Gamma_{i}$ $i=1,\ldots,N$ . The borehole is located at point $O$ .

In what follows at the modelling of the oil reservoir we neglect the compressibility of the fluid phase $(\rho=\mathrm{const})$ while taking into account the compressibility of the rock by assuming that the porosity depends linearly on the pressure [23]:

m=m_{0}\bigl{(}1+\varepsilon(p-p_{0})\bigr{)}

(2)

where $m_{0}$ and $p_{0}$ are the reference porosity and pressure, and $\varepsilon$ is the elastic capacity coefficient.

It is supposed that the horizontal wellbore $\Gamma_{0}$ is straight and is parallel to the $x$ -axis, has a fixed radius $R$ and given coordinates of the origin $O=(X_{w},Y_{w},Z_{w})$ (the bottom hole) and the end $E$ . Hydraulic fractures $\Gamma_{j}$ $(j=1,\ldots,N)$ are modelled by slots of fixed width $d_{j}$ as a planar simply connected figures (mostly, rectangles or ellipses), intersecting the horizontal wellbore $\Gamma_{0}$ at given angles. Position of each fracture $\Gamma_{j}$ with respect to the wellbore is individual, at that, fractures heights $h_{i}$ are less or equal to the reservoir thickness $H$ .

Interaction of the filtration flow with the well and the fractures is modelled by the boundary condition of the second kind:

q_{j}=-\left(\frac{k}{\mu}\frac{\partial p}{\partial\nu}\right)_{\Gamma_{j}},\quad j=0,\ldots,N.

(3)

Here $q_{j}$ is the velocity of the inflow to the wellbore or the fracture, $\nu$ is the outer (with respect to the domain $\Omega$ ) normal to the boundary $\Gamma_{j}$ .

The fluid flow along the wellbore is governed by the continuity equation

\frac{\partial(S\rho)}{\partial t}+\frac{\partial Q_{w}}{\partial s}=2\pi R\rho q_{0}.\quad 0<s<L_{w}.

(4)

Here $S=\pi R^{2}$ is the sectional area of the wellbore, $s$ is the local coordinate along the wellbore, $Q_{w}=S\rho v_{w}$ is the flow rate through the section of the wellbore, $v_{w}$ is the average fluid velocity along the wellbore, $L_{w}$ is the length of the wellbore. Assuming laminar flow in the wellbore, the fluid velocity $v_{w}$ is given by the Hagen-Poiseuille formula

v_{w}=-\gamma_{0}\frac{R^{2}}{8\mu}\frac{\partial p}{\partial s}.

(5)

Here $0<\gamma_{0}\leqslant 1$ is the imperfection coefficient that takes into account possible additional hydraulic resistance of wellbore due to imperfections. In the general 3-D model we do not distinguish between the reservoir and the wellbore pressure. The model with the individual wellbore pressure that would take into account the wellbore casing or formation of a filter cake on the wall will be observed in a separate paper. Under the assumption of invariability of the wellbore radius and incompressibility of fluid we obtain the relation between the pressure derivatives by the normal and along the wellbore.

\left.\frac{\partial p}{\partial\nu}\right|_{\Gamma_{0}}=\frac{\gamma_{0}R^{3}}{16k}\frac{\partial^{2}p}{\partial s^{2}}

(6)

For the flow in hydraulic fractures we use similar conservation laws:

\frac{\partial(d_{j}\rho)}{\partial t}+\mathrm{div\,}_{2}\mathbf{Q_{j}}=\rho q_{j},

(7)

The fluid flow rate is calculated by the formula $\mathbf{Q}_{j}=\rho\mathbf{q}_{j}$ with the total discharge $\mathbf{q}_{j}$ given by one of the following relations depending on the problem setup. For clean fracture of width $d_{j}$ we use the Poiseuille law for the flow in a slot between two plates [24]:

\mathbf{q}_{j}=-\frac{\gamma_{j}d_{j}^{3}}{12\mu}\nabla_{2}p.

(8)

This formula suits well for the injection wells. The case of “dirty” fracture filled with proppant is treated by the Darcy law

\mathbf{q}_{j}=-\frac{\gamma_{j}k_{j}d_{j}}{\mu}\nabla_{2}p.

(9)

Equation (9) can be used in case when proppant concentration in the fracture is such that the notion of fracture permeability $k_{j}$ becomes reasonable. The 2-D differential operators $\mathrm{div\,}_{2}$ and $\nabla_{2}$ are calculated in local coordinates for each fracture. In both cases $0<\gamma_{j}\leq 1$ is a coefficient characterizing imperfection of fracture walls. Both formulae (8) and (9) are unified by introduction of the flow rate coefficient

\theta_{j}=\left\{\begin{array}[]{ll}\frac{\gamma_{j}d_{j}^{3}}{12\mu}&\mbox{ --- clean fracture,}\\[5.69054pt] \frac{\gamma_{j}k_{j}d_{j}}{\mu}&\mbox{ --- dirty fracture},\end{array}\right.\quad j=1,\ldots,N.

(10)

By combining equations (7) and (10) we obtain the expression of the normal (to the fracture’s wall) derivative of pressure in terms of derivatives of pressure along the fracture:

\left.\frac{\partial p}{\partial\nu}\right|_{\Gamma_{j}}=\theta_{j}\Delta_{2}p.

(11)

Here $\Delta_{2}$ is a 2-D Laplace operator taken in fracture’s local coordinates. In the use of formula (11) one should pay attention to the fact that the fluid inflow should be calculated over both walls of the fracture.

2.1 Weak formulation of the problem

Let us write the weak statement of the problem. Due to the incompressibility of fluid $\rho=\mathrm{const}$ and using equation (2) we obtain

\frac{\partial m}{\partial t}=\frac{\partial}{\partial t}\Bigl{(}m_{0}(1+\varepsilon(p-p_{0}))\Bigr{)}=m_{0}\varepsilon\frac{\partial p}{\partial t}.

Equation (1) gives

m_{0}\varepsilon\frac{\partial p}{\partial t}-\mathrm{div\,}\left(\frac{k}{\mu}\nabla p\right)=0.

By multiplication of this equation to an arbitrary function $\psi(\mathbf{x})$ and integration over the domain $\Omega$ we find

0=\iiint\limits_{\Omega}\psi\left(m_{0}\varepsilon\frac{\partial p}{\partial t}-\text{div}\left(\frac{k}{\mu}\nabla p\right)\right)d\Omega=\\ =\iiint\limits_{\Omega}\left(m_{0}\varepsilon\psi\frac{\partial p}{\partial t}-\text{div}\left(\psi\frac{k}{\mu}\nabla p\right)+\frac{k}{\mu}\nabla\psi\nabla p\right)d\Omega=\\ =\iiint\limits_{\Omega}\left(m_{0}\varepsilon\psi\frac{\partial p}{\partial t}+\frac{k}{\mu}\nabla\psi\nabla p\right)d\Omega-\frac{k}{\mu}\iint\limits_{\partial\Omega}\psi\frac{\partial p}{\partial\nu}dS

(12)

Let us transform the boundary integral. Over the outer boundary $\Sigma$ depending on the problem’s formulation we can state either Neumann or Dirichlet condition, i.e. either the zero flow rate:

\frac{\partial p}{\partial\nu}\Big{|}_{\Sigma}=0,

or a given pressure:

p\Big{|}_{\Sigma}=p_{\infty},\mbox{ at that }\psi\Big{|}_{\Sigma}=0.

By using expressions for the normal derivative over the inner boundaries (6), (11), we find

\iint\limits_{\partial\Omega}\psi\frac{\partial p}{\partial\nu}dS=\iint\limits_{\Gamma_{0}\setminus\cup(\Gamma_{0}\cap\Gamma_{j})}\psi\frac{\partial p}{\partial\nu}dx+\sum\limits_{j=1}^{N}\iint\limits_{\Gamma_{j}}\psi\frac{\partial p}{\partial\nu}dS=\\ =\frac{\gamma_{0}R^{3}}{16k}\iint\limits_{\Gamma_{0}\setminus\cup(\Gamma_{0}\cap\Gamma_{j})}\psi\frac{\partial^{2}p}{\partial x^{2}}dS+\sum\limits_{j=1}^{N}\theta_{j}\iint\limits_{\Gamma_{j}}\psi\Delta_{2}p\,dS=\\ =\frac{\gamma_{0}\pi R^{4}}{8k}\psi\frac{\partial p}{\partial x}\Big{|}_{O}^{E}+\frac{\gamma_{0}\pi R^{4}}{8k}\sum\limits_{j=1}^{N}\psi\frac{\partial p}{\partial x}\Big{|}_{\Gamma_{j}^{+}\cap\Gamma_{0}}^{\Gamma_{j}^{-}\cap\Gamma_{0}}-\frac{\gamma_{0}R^{3}}{16k}\iint\limits_{\Gamma_{0}}\frac{\partial\psi}{\partial x}\frac{\partial p}{\partial x}dS+\\ +\sum\limits_{j=1}^{N}\theta_{j}\Big{(}\int\limits_{\Gamma_{j}\cap\Gamma_{0}}\psi\frac{\partial p}{\partial\nu}dl-\iint\limits_{\Gamma_{j}}\nabla_{2}\psi\cdot\nabla_{2}p\,dS\Big{)}.

(13)

Here by $\Gamma_{j}^{\pm}\cap\Gamma_{0}$ we denote circles obtained at the intersection of the side surface of $j$ -th fracture with the cylindrical wellbore $\Gamma_{0}$ . The terms calculated over these intersections, have the meaning of inflows form the fractures to the wellbore, and of the outflows into the wellbore. All together these inflows and outflows compensate each other, hence, the corresponding terms cancel out. At that, integrals over $\Gamma_{0}$ should be calculated without taking into account the gaps on the intersections of the fractures and the wellbore.

The fluid flow at the end of the fracture $E$ is assumed to vanish, hence,

\iint\limits_{\partial\Omega}\psi\frac{\partial p}{\partial\nu}dS=\frac{\mu}{k}Q_{0}(t)\psi\Big{|}_{O}-\frac{\gamma_{0}R^{3}}{16k}\iint\limits_{\Gamma_{0}}\frac{\partial\psi}{\partial x}\frac{\partial p}{\partial x}dS-\sum\limits_{j=1}^{N}\theta_{j}\iint\limits_{\Gamma_{j}}\nabla_{2}\psi\cdot\nabla_{2}p\,dS.

(14)

By integration over $\Gamma_{j}$ it is necessary to compute integral over both walls of the fracture.

The variational statement of the problem for the multiple fractured reservoir with the horizontal wellbore reads as follows: Find function $p\in W^{1,2}(\Omega)$ , satisfying the equation

\iiint\limits_{\Omega}\left(m_{0}\varepsilon\psi\frac{\partial p}{\partial t}+\frac{k}{\mu}\nabla\psi\nabla p\right)d\Omega-Q_{0}(t)\psi\Big{|}_{O}+\frac{\gamma_{0}R^{3}}{16\mu}\iint\limits_{\Gamma_{0}}\frac{\partial\psi}{\partial x}\frac{\partial p}{\partial x}dS\\ +\sum\limits_{j=1}^{N}\theta_{j}\iint\limits_{\Gamma_{j}}\nabla_{2}\psi\cdot\nabla_{2}p\,dS=0,\quad p|_{t=0}=p_{0}(\mathbf{x})

(15)

for every $\psi(\mathbf{x})\in W^{1,2}(\Omega)$ at every time $t\in[0,T]$ .

2.2 Dimensionless formulation

It is convenient to introduce the dimensionless variables

t^{\prime}=t/T,\quad p^{\prime}=p/p_{\ast},\quad q^{\prime}_{0}=TQ_{0}/(LR^{2}),\quad x^{\prime}=x/L,\quad y^{\prime}=y/L,\quad z^{\prime}=z/L.

Here $T$ and $p_{\ast}$ are some characteristic scales of time and fluid pressure. Writing the problem (15) in dimensionless variables leads to

\iiint\limits_{\Omega^{\prime}}\left(\psi\frac{\partial p^{\prime}}{\partial t^{\prime}}+a_{1}\nabla^{\prime}\psi\nabla^{\prime}p^{\prime}\right)d\Omega^{\prime}+a_{2}\iint\limits_{\Gamma_{0}^{\prime}}\frac{\partial\psi}{\partial x^{\prime}}\frac{\partial p^{\prime}}{\partial x^{\prime}}dS^{\prime}-a_{3}q^{\prime}_{0}(t)\psi\Big{|}_{O^{\prime}}\\ +\sum\limits_{j=1}^{N}a_{j+3}\iint\limits_{\Gamma^{\prime}_{j}}\nabla_{2}^{\prime}\psi\cdot\nabla_{2}^{\prime}p^{\prime}\,dS^{\prime}=0,\quad p^{\prime}|_{t=0}=p^{\prime}_{0}(\mathbf{x}),

(16)

where

a_{1}=\frac{Tk}{\mu m_{0}\varepsilon L^{2}},\quad\quad a_{2}=\frac{\gamma_{0}R^{3}T}{16\mu L^{3}m_{0}\varepsilon},\quad\quad a_{3}=\frac{R^{2}}{L^{2}m_{0}\varepsilon p_{\ast}},\quad a_{j+3}=\frac{\theta_{j}T}{L^{3}m_{0}\varepsilon}

(17)

In case of non-constant reservoir parameters there will be dimensionless multipliers under integrals in the corresponding terms of equation (16).

The problem under consideration has a small parameter: a ratio $H/L$ of reservoir’s width to its length. This observation allows us to construct a simplified 2-D model of the process as described in the subsequent sections.

3 Two-dimensional model

For the construction of a simplified 2-D model we will assume that all the hydraulic fractures extend from the bottom to the top of the reservoir: $h_{j}=H$ , $j=1,\ldots,N$ . In such a case the reservoir fluid flow is mainly $z$ -independent and two-dimensional. Let us introduce the operation of averaging along the vertical axis as

\bar{f}(t,x,y)=\frac{1}{H}\int\limits_{0}^{H}f(t,x,y,z)dz.

By substitution of a $z^{\prime}$ -independent test function $\psi(t^{\prime},x^{\prime},y^{\prime})$ in equation (16), one can perform integration along $Oz^{\prime}$ -axis in integrals over the domain $\Omega^{\prime}$ as well as in integrals over the fractures $\Gamma^{\prime}_{j}$ , $j=1,\ldots,N$ . As a result of integration, pressure $p^{\prime}$ transforms to the average pressure $\bar{p}^{\prime}$ with the multiplier $H/L$ . The source term at point $O^{\prime}$ in the averaged model has the same meaning and does not change. It is remained to average correctly the term, containing the integration over the wellbore $\Gamma_{0}^{\prime}$ .

In the main volume of the reservoir and in the vicinity of the hydraulic fractures the average pressure $\bar{p}$ does not differ much from the original pressure $p$ , but this is not true near the wellbore. At the averaging along the reservoir’s hight, the wellbore is substituted by a fictitious fracture extended from the bottom to the top of the reservoir. Parameters of the fictitious fracture should be selected in order to satisfy two conditions: (a) Equality of inflows of the reservoir’s fluid to the fictitious fracture and to the original wellbore; (b) Equality of hydrodynamical resistance for the flow along the fictitious fracture and along the wellbore.

Condition (b) is easily satisfied by a proper choice of the width of the fictitious fracture by equating the flow rates for flow in the fracture and in the wellbore with the same pressure gradient using formulae (5) and (8):

\frac{\pi R^{4}}{8\mu}=\frac{Hd_{0}^{3}}{12\mu}\quad\Rightarrow\quad d_{0}=\left(\frac{3\pi R^{4}}{2H}\right)^{1/3}.

On the contrary, condition (a) in case of non-zero inflow can not be satisfied in terms of the present one-pressure model. Indeed, the inflow from the reservoir is determined by the difference of the pressure in the wellbore (fracture) and at some distant point. In case of the wellbore there is a logarithmic singularity at the fracture’s axis, whereas in case of the inflow to the fracture the pressure is continuous. This implies, that fluid inflow to the fracture requires much smaller pressure gradient between the fracture and the reservoir than the same inflow to the wellbore. However, due to the coincidence of the hydrodynamical resistance, pressure gradients providing the flow along the fracture and the wellbore are the same. Since the pressure at the origin of the wellbore is known, it is impossible to choose the required pressure along the fictitious fracture that guarantees implementation of both conditions (a) and (b).

For the solution of the declared problem we will consider pressure in the wellbore separately from the reservoir’s pressure.

3.1 Computation of the wellbore pressure

Let us denote the fluid pressure in the wellbore by symbol $U$ . Instead of equation (3) for $j=0$ we will use the relation

2\pi Rq_{0}=-2\pi R\left(\frac{k}{\mu}\frac{\partial\bar{p}}{\partial\nu}\right)_{\Gamma_{j}}=\alpha_{0}(\bar{p}-U)\Bigr{|}_{\Gamma_{j}}

(18)

that declares that the fluid inflow to the wellbore is proportional to the difference between the averaged reservoir’s pressure $\bar{p}$ and the wellbore pressure $U$ . Coefficient $\alpha_{0}$ is chosen from the condition of equality of inflows to the planar fracture under the action of the averaged pressure $\bar{p}$ , and to the wellbore due to the original 3-D pressure:

\alpha_{0}=-\frac{2\pi k}{\mu\ln\left(\frac{2R\pi}{H}\sin\bigl{(}\frac{Z_{w}\pi}{H}\bigr{)}\right)}.

Details of computation of coefficient $\alpha_{0}$ are given in A. Equation (4) remain unchanged, whereas in equation (5) symbol $p$ should be substituted by $U$ .

In the new definition of the problem, the fluid flow along the wellbore is driven by the gradient of pressure $U$ , which satisfies the continuity equation in the following form

\frac{\partial}{\partial s}\left(-\frac{\gamma_{0}\,\pi\,R^{4}}{8\mu}\frac{\partial U}{\partial s}\right)=\alpha_{0}\bigl{(}\bar{p}-U\bigr{)}+\sum\limits_{j=1}^{N}\alpha_{j}d_{j}(\bar{p}-U)\delta(s-s_{j}).

(19)

Here by $\alpha_{j}$ we denote the proportionality coefficients for the wellbore ( $j=0$ ) and fractures ( $j=1,\ldots,N$ ) computed according to A as follows:

\alpha_{j}=-\frac{2\pi\theta_{j}}{d_{j}p\ln\left(\frac{2R\pi}{H}\sin\bigl{(}\frac{Z_{w}\pi}{H}\bigr{)}\right)},\quad j=1,\ldots,N.

Pressure $\bar{p}$ in equation (19) should be taken at the corresponding point of the wellbore. Values $s_{j}$ in the argument of the Dirac delta-function $\delta$ are coordinates of the hydraulic fractures in the local coordinates along the wellbore.

Boundary conditions to equation (19) can be chosen either as a given flow rate or a given pressure at the origin $O$ . At the end $E$ of the wellbore we assume zero flow rate:

\left.\frac{\partial U}{\partial s}\right|_{s=0}=-\frac{8\mu Q(t)}{\pi R^{4}\gamma_{0}},\mbox{ or }U|_{s=0}=p_{w}(t);\quad\left.\frac{\partial U}{\partial s}\right|_{s=L_{w}}=0.

(20)

Let us choose the following dimensionless variables:

U^{\prime}=U/p_{\ast},\quad\delta^{\prime}=L\delta.

Equation (19) and boundary conditions (20) in the dimensionless variables take the form (primes are omitted):

\begin{array}[]{l}\displaystyle-\frac{\partial^{2}U}{\partial s^{2}}=b_{1}(\bar{p}-U)+\sum\limits_{j=1}^{N}b_{j+1}(\bar{p}-U)\delta(s-s_{j}),\\[11.38109pt] -\displaystyle\left.\frac{\partial U}{\partial s}\right|_{s=0}=b_{0}q_{0}(t),\quad U|_{s=L_{w}}=\bar{p}\end{array}

(21)

Dimensionless parameters $b_{j}$ have the following expressions in terms of the dimensional variables:

b_{0}=\frac{8\,\mu\,L^{2}}{\gamma_{0}\,\pi\,R^{2}\,T\,p_{\ast}},\quad b_{1}=\frac{8\,\mu\,L^{2}\,\alpha_{0}}{\gamma_{0}\,\pi\,R^{4}},\quad b_{j+1}=\frac{8\,\mu\,d_{j}\,L\,\alpha_{j}}{\gamma_{0}\,\pi\,R^{4}}.

Equation (21) describes the pressure distribution in the wellbore and takes into account the fluid inflow from the reservoir and from the fractures.

3.2 Weak form of the 2-D model

For the numerical calculations it is convenient to use the weak form of the problem. Formula (14) for the fluid inflow to the wellbore and to the fractures is modified as follows:

\iint\limits_{\partial\Omega}\psi\frac{\partial\bar{p}}{\partial\nu}dS=-\frac{\mu}{k}\int\limits_{O}^{E}\alpha_{0}\psi(\bar{p}-U)ds-\frac{\mu}{k}\sum\limits_{j=1}^{N}\alpha_{j}d_{j}\psi(\bar{p}-U)|_{\Gamma_{j}\cap\Gamma_{0}}-\sum\limits_{j=1}^{N}\theta_{j}H\int\limits_{\Gamma_{j}}\frac{\partial\psi}{\partial s}\frac{\partial\bar{p}}{\partial s}\,ds.

(22)

Thus, we obtain the following equation for the averaged pressure $\bar{p}$ in the dimensionless variables (primes are omitted):

\iint\limits_{\Omega_{2}}\left(\psi\frac{\partial\bar{p}}{\partial t}+a_{1}\nabla\psi\cdot\nabla\bar{p}\right)d\Omega+a_{2}\int\limits_{O}^{E}\psi(\bar{p}-U)ds+\sum\limits_{j=1}^{N}\left(a_{2j+1}\psi(\bar{p}-U)|_{\Gamma_{j}\cap\Gamma_{0}}+a_{2j+2}\int\limits_{\Gamma_{j}}\frac{\partial\psi}{\partial s}\frac{\partial\bar{p}}{\partial s}\,ds\right)=0,

(23)

where

a_{1}=\frac{Tk}{\mu m_{0}\varepsilon L^{2}},\quad a_{2}=\frac{\alpha_{0}T}{m_{0}\varepsilon HL},\quad a_{2j+1}=\frac{\alpha_{j}d_{j}T}{m_{0}\varepsilon HL^{2}},\quad a_{2j+2}=\frac{\theta_{j}T}{L^{3}m_{0}\varepsilon}.

(24)

Here by $\Omega_{2}$ and $\Gamma_{j}$ , $j=0,\ldots,N$ we denote projections of domain $\Omega$ , of the wellbore and of the fractures to $Oxy$ plane, represented by a rectangle and lines respectively. Equation (23) can be modified with the help of equation (21) as

\iint\limits_{\Omega_{2}}\left(\psi\frac{\partial\bar{p}}{\partial t}+a_{1}\nabla\psi\nabla\bar{p}\right)d\Omega+c_{1}\int\limits_{O}^{E}\frac{\partial\psi}{\partial s}\frac{\partial U}{\partial s}ds-c_{0}q(t)\psi(O)+\sum\limits_{j=1}^{N}a_{2j+2}\int\limits_{\Gamma_{j}}\frac{\partial\psi}{\partial s}\frac{\partial\bar{p}}{\partial s}\,ds=0,

(25)

where

c_{0}=\frac{R^{2}}{m_{0}\varepsilon HLp_{\ast}},\quad c_{1}=\frac{\gamma_{0}\pi R^{4}T}{8\mu m_{0}\varepsilon HL^{3}}.

(26)

Multiplication of equation (21) to the test function $\eta$ and integration on $s$ from $O$ to $E$ with the use of the boundary conditions (20) leads to

\begin{array}[]{l}\displaystyle\int\limits_{O}^{E}\frac{\partial\eta}{\partial s}\frac{\partial U}{\partial s}ds-b_{1}\int\limits_{O}^{E}\eta(\bar{p}-U)ds-\sum\limits_{j=1}^{N}b_{j+1}(\bar{p}-U)\eta|_{s=s_{j}}-b_{0}q_{0}(t)\eta|_{s=0}=0.\end{array}

(27)

The weak formulation of the 2-D problem is now states as follows: Find function $p\in W^{1,2}(\Omega_{2})$ and function $U\in W^{1,2}(\Gamma_{0})$ satisfying equations (23) (or (25)) and (27) for every functions $\psi\in W^{1,2}(\Omega_{2})$ and $\eta\in W^{1,2}(\Gamma_{0})$ at every time $t\in[0,T]$ .

4 Case studies

The purpose of the following numerical experiments is to demonstrate the agreement of our numerical results with known data and adequacy of the behaviour of the main flow parameters as functions of the geometry and physical properties of the problem.

4.1 Numerical realization

Numerical calculations are performed in a freely accessible finite element solver FreeFEM++ [16]. For the calculations we use weak dimensionless formulation (23)–(27). Symmetry of the problem with respect to the unknown functions $(p,U)$ and test functions $(\psi,\eta)$ assures the symmetry of the stiffness matrix which is important for the correctness of the numerical algorithm.

Time derivative is approximated by a finite difference: $\partial p/\partial t\approx(p^{n+1}-p^{n})/\tau$ where $\tau$ is the time step. The upper index denotes the time instant: $p^{n}=p(t_{n},\mathbf{x})$ , $t_{n}=n\tau$ . Time iterations are performed using the first-order backward (implicit) Euler method. Due to unconditional stability of the method, the value of the time step $\tau$ does not depend on the characteristic diameter of the mesh cells.

For the construction of a numerical mesh we take 200 mesh vertices on the outer border $\Sigma$ , and a proportional number of vertices over the wellbore and fractures. The total number of mesh vertices varies from 6000 to 13000 depending on the complexity of the geometry. Dimensionless time step is equal to 0.05 which is equivalent to 1 hour 12 minutes. In all cases time of calculation on 4000 time steps was about 10–15 of minutes on a usual PC.

4.2 Common parameters

The main reservoir and fluid parameters used in numerical tests are listed in Table 1. In all tests we prescribe the borehole pressure $p|_{O}=p_{w}$ and calculate the volume of produced fluid using formula (27) with $\eta=1$ . Unless otherwise mentioned, we use the Neumann’s (no-flow) outer boundary condition: $\partial p/\partial n|_{\Sigma}=0$ . Initial pressure in undisturbed reservoir is assumed to be zero: $p(0,\mathbf{x})=0$ . In our model there is no difference between production and injection wells, therefore for the sake of convenience we assume positive pressures by taking $p_{w}>0$ .

Table 1: Parameters set for the reservoir

Parameter	Symbol	Value
reservoir size	$L_{x}\times L_{y}\times L_{z}$	$2800\times 2600\times 20$ m
wellbore origin	$O=(X_{w},Y_{w},Z_{w})$	$(800,1300,10)$ m
wellbore radius	$R$	10 cm
wellbore length	$L_{w}$	1100 m
fracture aperture	$d_{j}$	1 cm, $j=1,\ldots,N$
fracture permeability	$k_{j}$	1000 D $\approx 0.987\cdot 10^{-9}\mathrm{m}^{2}$ , $j=1,\ldots,N$
borehole pressure	$p_{w}$	10 MPa
pressure at infinity	$p_{\infty}$	0 MPa
fluid dynamic viscosity	$\mu$	1 cP $\approx 10^{-3}\mathrm{Pa}\cdot\mathrm{sec}$
porosity	$m_{0}$	0.1
elastic capacity coefficient	$\varepsilon$	1 GPa^-1
imperfection coefficient	$\gamma_{j}$	$\in(0,1]$ , $j=0,\ldots,N$

Table 2: Set of instrumental parameters for Case 1

index	rock permeability $k$ , mD	imperfection coefficient $\gamma_{0}$	number of fractures $N$
0	1	1	0
1	10	$3\cdot 10^{-5}$	2

4.3 Case 1: Influence of rock permeability, well conductivity and presence of fractures

We begin with the test that demonstrate the influence of the rock permeability, wellbore conductivity and presence of fractures to the well production. We number the observed cases by a 3-digits binary number as indexed in Table 2 where the instrument parameters are listed. For example, the case #101 corresponds to $k=10$ mD, $\gamma_{0}=1$ , $N=4$ . In cases of fractured wellbore we consider two linear fractures intersecting the wellbore at right angle as shown in Figure 2. The first fracture has symmetrical wings of length 260 m each, the second fracture has non-symmetrical wings of 260 an 305 m. The imperfection coefficients are $\gamma_{1}=0.5$ for the first fracture and $\gamma_{2}=1$ for the second fracture.

The calculated volume of produced fluid in all 8 cases during the period of 180 days is shown in Figure 3. Higher productivity is expectably obtained in cases with hydraulic fracture #xx1 in contrast to cases of a single wellbore #xx0.

Lower values of imperfection coefficient $\gamma_{0}$ correspond to higher hydraulic resistance of the wellbore which limits the volume of produced fluid. Difference in the wellbore conductivity affects the well productivity weakly in cases of the lower rock permeability (pairs 000 — 010, 000 — 011) and strongly in cases of higher permeability (pairs 100 — 110, 101 — 111). The influence of the rock permeability is also expectable: higher volumes of fluid are obtained for the higher permeabilities.

Figure 3 contains graphs that are drawn with the same data as some cases observed in [13], namely, cases #110 and #111 correspond to graphs number 2 in Figures 2 and 3 in [13] respectively. One can see a good agreement between these pairs of graphs which confirms the reliability of our numerical algorithm.

4.4 Case 2: Influence of the relative positions of fractures

In this set of cases we vary the number of fractures, lengthes and distances between fractures. We observe three cases: (a) four symmetrical fractures on the equal intervals of 220 m, each of the total length (both wings) of 1040 m; (b) two fractures separated by the intervals of 367 m, each of the total length 2080 m; (c) four symmetrical fractures separated by the shorter interval of 110 m of the total length 1040 m each. In all cases the rock permeability is $k=1$ mD and imperfection coefficients are $\gamma_{j}=1$ , $j=0,\ldots,4$ , The idea is to compare the productivity of fractures of the same total length 4160 m in different geometrical setups.

Figure 4 demonstrates the volume of produced fluid versus time in each of the observed cases. One can see that cases (a) and (c) do not distinguish within first 10 days of production until the drainage zones of fractures start to influence each other. In case (c) of denser fractures the mutual influence of the drainage zones decreases the productivity in comparison to case (a) of fractures placed more rare. In contrast, two long fractures in case (b) of the same total length as in case (a) are less productive during first 40 days due to the finite conductivity of fractures. However, for longer time two long fractures are more productive than 4 shorter ones due to the larger area of the drainage zone. Pressure distribution for all three cases at $t=60$ days is shown in Figure 5 where the difference of the drainage zones is clearly seen.

Thus, this example demonstrate the potential of the model for optimization of the fractures network for given reservoir conditions.

4.5 Case 3: An “arbitrary” fracture net

The model allows specifying any 2D net of fractures (also with self-intersections) with different conductivities of its segments. In this case we observe a set of fractures shown in Figure 6. The produced volume of fluid is compared for different rock permeabilities, reservoir heights and outer boundary conditions. The data set enumerated by a binary 3-digits number as given in Table 3. We use two types of outer boundary conditions: Neumann (no-flow): $\partial p/\partial n|_{\Sigma}=0$ and Dirichlet (given pressure): $p|_{\Sigma}=0$ .

The produced volume of fluid is shown in Figure 7. The difference in outer boundary conditions is clearly seen for high permeability (compare cases 100–110, 101–111). For the Neumann outer boundary condition the reservoir contains only a finite volume of fluid, which implies that for large time the produced volume is limited by the total volume in contrast with the case of the Dirichlet boundary condition where the total volume is unlimited and the produced volume tends to a linear function of time. Difference in the reservoir’s heights gives the proportional increase of the productivity as is seen from the comparison of cases #xx0 and #xx1. The dependence on the permeability is also expectable: higher volumes of produced fluid correspond to higher permeability.

Table 3: Set of instrumental parameters for Case 3

index	rock permeability $k$ , mD	outer boundary conditions	reservoir’s hight $L_{z}$
0	1	no-flow	20 m
1	10	given pressure	10 m

Pressure distribution for cases #000 and #111 at $t=5$ , 30, 50 and 180 days is shown in Figure 8. One can see that the depression zone spreads faster in case of the higher rock permeability. The Dirichlet outer boundary condition in case #111 allows a fluid inflow to the reservoir from the surrounding media. Together with the high permeability this brings the flow to the stationary regime starting at $t\approx 40$ h as it follows from the comparison of the bottom-right column of pictures in Figure 8 and from the corresponding productivity curve in Figure 7. In contrast, the Neumann boundary condition in case #000 leads to the leveling of pressure for higher $t$ .

5 Conclusion

The mathematical model of fluid flow in a rectangular reservoir drained by a horizontal wellbore with an arbitrary network of vertical hydraulic fractures is proposed. The model allows computation of the time-dependent fluid pressure distribution within all segments of the computational domain (reservoir, fractures and wellbore) as well as integral characteristics of the flow such as the wellbore productivity. The model is applicable for reservoirs with variable in time and space physical properties (i.e. rock permeability, porosity, imperfection coefficients, etc.) and various scenarios of production (prescribed borehole pressure or flow rate as a function of time). The model does not use restrictive assumptions regarding the structure of the fluid flow and involve only a limited set of empirical parameters (wellbore and fractures imperfection coefficients). The weak formulation of the model given in the paper is suitable for the numerical solution of the model using the finite element method.

In case when all fractures extend from the bottom to the top of the reservoir, the model is reduced to a two-dimensional one by averaging along the vertical coordinate. The 2D model is implemented in a numerical code using the finite element solver FreeFEM++. Case studies for the 2D model demonstrate the dependence of the wellbore productivity on the physical and geometrical characteristics of the reservoir, of the outer boundary conditions and of the characteristics of the fracture network. This analysis suggests further applications of the model as an optimization tool for estimation of wellbore productivity under different reservoir development plans. Ability of the model to produce the fluid pressure field in the reservoir for various production scenarios allows one to compute initial data for industrial reservoir simulators. The model can also be used for the direct examination of validity of hypotheses regarding the character of fluid flow in different stages of the reservoir development, used in analytical and semi-analytical models that are cited in Introduction.

6 Acknowledgements

The work was partially supported by President Grant for Leading Scientific Schools of the Russian Federation (grant N.Sh.-2133.2014.1) and by RFBR (grant 16-01-00610).

Appendix A Computation of coefficients $\alpha_{j}$ for the 2-D model

For the computation of coefficient $\alpha_{j}$ let us assume a problem of filtration of incompressible fluid in an infinite stripe of width $H$ between two impermeable planar walls, generated by either a point or a slit sink of a given intensity $Q$ . The geometry of the problems is described in Figure 9.

The slit sink coincide with the interval $OT$ over $Oz$ -axis. The point sink $S$ is located on the hight $a$ from the bottom of the stripe.

The flow is governed by the relation $\mathbf{v}=-\sigma\nabla p$ ( $\sigma=k/\mu$ for the fracture with proppant or $\sigma=d^{2}/(12\mu)$ for the clean fracture of width $d$ ) and the continuity equation $\mathrm{div\,}\mathbf{v}=0$ . Thus, pressure $p$ satisfies the Laplace equation with no-flow conditions over the boundaries of the stripe. The flow rate at the infinity is given as

\frac{\partial p}{\partial z}\Bigr{|}_{z=0,H}=0,\quad-H\sigma\frac{\partial p}{\partial y}\Bigr{|}_{|y|\to\infty}=\frac{Q}{2}.

The exact solution of the problem on a slit sink in a stripe has the form

p_{f}=\frac{Q}{\sigma}\left(\frac{|y|}{2H}-\frac{\ln 2}{\pi}\right).

The solution of the problem with a point sink reads as follows:

p_{s}=\frac{Q}{4\pi`\sigma}\left[\ln\left(\sinh^{2}\Bigl{(}\frac{\pi y}{2H}\Bigr{)}+\sin^{2}\Bigl{(}\frac{\pi(z-a)}{2H}\Bigr{)}\right)+\ln\left(\sinh^{2}\Bigl{(}\frac{\pi y}{2H}\Bigr{)}+\sin^{2}\Bigl{(}\frac{\pi(z+a)}{2H}\Bigr{)}\right)\right]

The constant arbitrariness in the definition of function $p$ is exploited to satisfy the matching condition: $\lim\limits_{|y|\to\infty}(p_{f}(y,z)-p_{s}(y,z))=0$ . Let us calculate the asymptotical behaviour of function $p_{s}$ at $r=\sqrt{y^{2}+(z-a)^{2}}\to 0$ :

p_{s}\bigr{|}_{r\to 0}=\frac{Q}{2\pi\sigma}\left[\ln r+\ln\left(\frac{\pi}{2H}\sin\Bigl{(}\frac{a\pi}{H}\Bigr{)}\right)\right]+\frac{Q(z-a)}{4\sigma H}\cot\Bigl{(}\frac{a\pi}{H}\Bigr{)}+O(r^{2})

We require the following condition to be satisfied:

Q=\alpha\bigl{(}p_{f}|_{y=0}-{p_{s}}|_{r=R}\bigr{)}

accurate to the small terms of order $R/H$ . After cancellation of $Q$ accurate to the small terms we obtain

1=\alpha\sigma\left[-\frac{\ln 2}{\pi}-\frac{\ln R}{2\pi}-\frac{1}{2\pi}\ln\left(\frac{\pi}{2H}\sin\Bigl{(}\frac{a\pi}{H}\Bigr{)}\right)\right].

This finally gives the formula for $\alpha$ as

\alpha=-\frac{2\pi\sigma}{\ln\left(\frac{2R\pi}{H}\sin\bigl{(}\frac{a\pi}{H}\bigr{)}\right)}.

References

[1] H. Sadrpanah, T. Charles, J. Fulton Explicit Simulation of Multiple Hydraulec Fractures in Horizontal Wells. 2006, SPE 99575
[2] B. Cvetkovic, G. Halvorsen, J. Sagen, E.N. Rigatos, Modelling the productivity of a multifractured horizontal well. SPE 71076, 2001.
[3] B. Guo, X. Yu, M. Khoshgahdam, A simple analytical model for predicting productivity of multifractured horizontal wells. SPE 114452-PA, 2009.
[4] T. M. Hegre, L. Larsen, Productivity of multifractured horizontal wells. SPE 28845, 1994.
[5] H. Shojaei, E. S. Tajer, Analytical solution of transient multiphase flow to a horizontal well with multiple hydraulic fractures. SPE 165703, 2013.
[6] J. Ren, P. Guo Performance of vertical fractured wells with multiple finite-conductivity fractures. J. Geophys. Eng. 2015. V. 12. No. 6. Pp. 978 -987.
[7] M. Al-Kobaisi, E. Ozkan, and H. Kazemi, A hybrid numerical-analytical model of finite-conductivity vertical fractures intercepted by a horizontal well. SPE 92040, 2004.
[8] J. Wan, K. Aziz, Semi-analytical well model of horizontal wells with multiple hydraulic fractures. SPE 81190-PA, 2002.
[9] W. Zhou, R. Banerjee, B.D. Poe, J. Spath, M. Thambynayagam, Semianalytical production simulation of complex hydraulic-fracture networks. SPE 157367-PA, 2014.
[10] Sh. Yao, F. Zeng, H. Liu, G. Zhao. A semi-analytical model for multi-stage fractured horizontal wells. Journal of Hydrology. 2013. V. 507. Pp. 201- 212.
[11] H.-T. Wang. Performance of multiple fractured horizontal wells in shale gas reservoirs with consideration of multiple mechanisms. Journal of Hydrology. 2014. V. 510. Pp. 299 -312.
[12] Yu-L. Zhao, L.-H. Zhang, J.-X. Luo, B.-N. Zhang. Performance of fractured horizontal well with stimulated reservoir volume in unconventional gas reservoir. Journal of Hydrology. 2014. V. 512. Pp. 447- 456.
[13] Kashevarov A. A., Hydraulic model for oil and gas reservoirs, Journal of Applied Mechanics and Technical Physics, Vol. 51, No. 6, pp. 868 – 876, 2010.
[14] Kuznetsov, D. S., Cheremisin, A. N., Chesnokov, A., Golovin, S. Advanced horizontal well model. 2010. In North Africa Technical Conference and Exhibition. Society of Petroleum Engineers. SPE-127742-MS.
[15] V. Martin, J. Jaffre, J.E. Roberts. Modeling fractures and barriers as interfaces for flow in porous media. SIAM J. Sci. Copmput. 2005. Vol. 26, No. 5, Pp. 1667 -1691
[16] Hecht F. New development in FreeFem++. J. Numer. Math. 2012. V. 20. No. 3-4, 251 – 265. 65Y15
[17] S. Kocberber, Finite-element black oil simulation system for heterogeneous reservoirs with horizontal wells having vertical hydraulic fractures. SPE 25269, 1993.
[18] A. Zerzar and Y. Bettam, Interpretation of multiple hydraulically fractured horizontal wells in closed systems. In proceedings of the SPE International Improved Oil Recovery Conference in Asia Pacific (IIORC 03), Pp. 295 307, 2003.
[19] Y. Jiang, J. Zhao, Y. Li, H. Jia, L. Zhang, Extended finite element method for predicting productivity of multifractured horizontal wells. Mathematical Problems in Engineering. V. 2014, ID 810493, 9 pages.
[20] Z. Lei, S. Cheng, X. Li, H. Xiao, A new method for prediction of productivity of fractured horizontal wells based on non-steady flow. Journal of Hydrodynamics, Ser. B. V. 19, Iss. 4, Pp. 494 500.
[21] R. Vicente, T. Ertekin, Modeling of coupled reservoir and multifractured horizontal well flow dynamics. SPE 101929, 2006.
[22] Guo J., Zhang L., Wang H., Feng G. Pressure Transient Analysis for Multi-stage fractured horizontal wells in shale gas reservoirs. Transport in Porous Media. 2012. V. 93. Pp.635 653.
[23] Biot, M.A., Theory of elasticity and consolidation for a porous anisotropic solid. J. Appl. Phys. 1955. 26 (2), 182 185.
[24] Batchelor, G. K. An Introduction to Fluid Dynamics. 1967 Cambridge University Press.

Fluid flow in a reservoir drained by a multiple fractured horizontal well

Abstract

keywords:

1 Introduction

2 Formulation of the mathematical model

2.1 Weak formulation of the problem

2.2 Dimensionless formulation

3 Two-dimensional model

3.1 Computation of the wellbore pressure

3.2 Weak form of the 2-D model

4 Case studies

4.1 Numerical realization

4.2 Common parameters

4.3 Case 1: Influence of rock permeability, well conductivity and presence of fractures

4.4 Case 2: Influence of the relative positions of fractures

4.5 Case 3: An “arbitrary” fracture net

5 Conclusion

6 Acknowledgements

Appendix A Computation of coefficients αj\alpha_{j} for the 2-D model

References

Appendix A Computation of coefficients $\alpha_{j}$ for the 2-D model