- Email: [email protected]
New approach for describing transient pressure response generated by horizontal wells of arbitrary geometry
Applied Numerical Mathematics 40 (2002) 433–449 www.elsevier.com/locate/apnum
New approach for describing transient pressure response generated by horizontal wells of arbitrary geometry Reinaldo J. González-Requena a,∗ , Juan M. Guevara-Jordan b a Escuela de Ingeniería de Petróleo, Facultad de Ingeniería, Universidad Central de Venezuela, Caracas, Venezuela b Escuela de Matemáticas, Facultad de Ciencias, Universidad Central de Venezuela, Caracas, Venezuela
Abstract This paper is aimed at developing a methodology for studying the transient pressure behavior of horizontal wells with any curvilinear trajectory in an isotropic/anisotropic arbitrarily shaped reservoir. This methodology employs generalized functions to represent the tortuous horizontal well. A particular way of removing the singularities involved in the partial differential equation is based on reducing the original problem to the conventional solution of the homogeneous diffusivity equation under any given initial and boundary conditions. The Green function method and any standard numerical technique are combined in a single computational strategy to obtain the transient pressure response generated by a curved and twisted horizontal well in reservoirs with irregular boundaries. Analytical methods can be also used, whenever possible, to solve the reduced problem. This proposal can be easily broadened to analyze the performance of the pressure transient of any kind of reservoir sources or sinks that can be modeled using generalized functions. Some models are presented. 2001 Published by Elsevier Science B.V. on behalf of IMACS. Keywords: Horizontal well; Sources; Mathematical modeling; Finite elements method; Diffusivity equation; Transient pressure; Generalized functions; Curved/twisted trajectory
1. Introduction The number of horizontal wells has increased considerably since the early 1980s within the petroleum industry. Over the past twenty years, the advances in drilling and production techniques have paved the way to the development, study, and adoption of new models for representing horizontal wells which are more appropriate for the analysis of the well performance. This work presents a mathematical model of the transient pressure response generated by a horizontal and/or tortuous well which explicitly considers the modeled well as a generalized function in the * Corresponding author.
E-mail address: [email protected] (R.J. González-Requena). 0168-9274/01/$22.00 2001 Published by Elsevier Science B.V. on behalf of IMACS. PII: S 0 1 6 8 - 9 2 7 4 ( 0 1 ) 0 0 0 9 4 - 0
434
R.J. González-Requena, J.M. Guevara-Jordan / Applied Numerical Mathematics 40 (2002) 433–449
non-homogeneous diffusivity equation [2,19,35]. This modeling procedure is not naturalsince a well should be represented as a boundary condition [12,13,35]. For instance, when fluids are subtracted from reservoirs using production wells, they are characterized by intensive convergent flow patterns and pressure field singularities. This has led to the representation of wells by means of generalized functions or distributions. One of the purposes of the present paper is the explicit use of this type of functions in the diffusivity equation when representing a well of arbitrary trajectory. Evidence of the implicit use of such generalized functions for representing vertical, horizontal, and fractured wells can be found in the petroleum literature. For instance, the pressure distribution solutions in infinite and rectangular reservoirs can be seen in [3–8,21,31–34,37]. Likewise, Ligget and Liu [28] have used conveniently redefined Dirac delta functions to represent vertical wells, at the expenses of certain level of rigorousness in the generalized function model. On the other hand, Raghavan and Ozkan [35] provide solutions to many reservoir cases (in rectangular coordinates) based on a similar distributional approach. A specially designed numerical method is necessary to solve the well model developed in the case of 2D or 3D reservoirs with arbitrary boundaries due to the inherent difficulties involved in representing generalized functions using traditional numerical methods like finite differences or finite elements [27]. Therefore we hereby propose an algorithm that combines the removal of singularities (singularity programming), Green functions, and the finite elements method in a single computational strategy to numerically solve the Initial and Boundary Value (IBV) problem. This numerical method, hereafter called the Singularity Programming Method (SPM), is not new in its applicability. Different authors have used it to solve elliptic equations [1,15,18,23,36,38]. For the purposes of the diffusivity equation, the singularity removal method is not as direct due to the additional time nature of the Dirac delta function used in the partial differential equation, and therefore it must be reformulated. Astudillo and Guédez [2] originally applied this reformulated methodology in the case of vertical wells in rectangular domains. They used the finite differences method as a numerical approach. In this sense, one of the major contributions of the present paper is to extend and validate the applicability of SPM to the mathematical model proposed for horizontal or tortuous wells in arbitrarily shaped reservoirs. The authors can assure that this extension of the SPM has not been previously recorded in the specialized literature. The pressure distribution created by a “fictitious” horizontal well with a trajectory given by a circular helix will be ascertained in this paper. Such potential use is another contribution of this methodology. So far the cases of the most tortuous trajectory well found by the authors in the specialized literature are that of a curved, torsionless horizontal well proposed by Koh and Tiab [26], and those, also torsionless, wells, reviewed by Azar-Nejad et al. [3]. In addition to its potential for really addressing arbitrary well geometries, this methodology can be extended to represent other reservoir anomalies such as fractures. Section 2 presents some basic mathematical and physical considerations necessary to use generalized functions when modeling the pressure response of horizontal wells. Section 3 considers the fundamental solution of the diffusivity operator together with the associated source term solutions. They will be used to program the singularity removal algorithm in Section 4. Section 5, probably the most important, presents numerical issues used to account for the usual wellbore conditions. Section 6 focuses on examples and their analysis. Section 7 presents the conclusions and final considerations.
R.J. González-Requena, J.M. Guevara-Jordan / Applied Numerical Mathematics 40 (2002) 433–449
435
2. Mathematical model According to the standard considerations [3] about rock/fluid conditions, and taking its dimensionless form, the basic equation describing the continuous pressure distribution generated by a well of arbitrary trajectory γ can be simplified as follows [19,35,40]: ∂p − p = ∂t
b
q γ (u), t δ x − γ (u) γ (u) du
(1)
a
where p is pressure, t is time, q is production rate at the wellbore per length unit and per unit of time [25], δ is the Dirac delta function [25,41–43], and the line integral is calculated along the curve γ used for modeling the well trajectory. The analytical representation of this curve is γ (u) = (γ1 (u), γ2 (u), γ3 (u)). Here, γ (u) is the Euclidean norm of the velocity vector, which is tangent to the curve at γ (u). 1 It is worth pointing out that for the dimensioning of (1), the gravity term was considered through the Hubbert potential [3–6]. For simplicity’s sake, the symbols used in Eq. (1) are also used to represent dimensionless variables, such as time, space and Hubbert potential. In rigor, p is the dimensionless form of the potential, but from this point on we will refer to it just as pressure. Modeling and simulating the behavior of pressure generated by a tortuous well in a reservoir Ω requires Eq. (1) to be solved subject to certain conditions at the reservoir boundary ∂Ω with a given initial condition. For the purposes of the present paper, we shall take into consideration only the absence of flow in ∂Ω, i.e., a zero Neumann condition. Nevertheless, the method is perfectly applicable to the cases of constant pressure condition (dirichlet condition) or mixed conditions. Therefore, the following problem with initial and boundary value (IBV) is assessed: ∂p − p = ∂t
b
q γ (u), t δ x − γ (u) γ (u) du
in Ω × (0, ∞),
(2a)
in ∂Ω × (0, ∞),
(2b)
in Ω,
(2c)
a
∂p (x, t) = 0, ∂n p(x, 0) = 0,
where the semi-infinite interval is associated to time and the zero initial condition is consistent with the dimensioning used for these studies. The source term in Eqs. (1) and (2a) is an idealization of the real condition of the wellbore, which can be proposed as an inner boundary condition. In this approach, the well, Dε , is conceived as a tube of diameter ε around the axis trajectory γ , and the boundary condition equals the flow rate of the fluid through the tube surface ∂Dε with the known production rate of the well. This representation of the well is very close to the physical reality, but is not used in reservoir simulation. Engineers prefer to neglect the well tubing because its diameter is very small in comparison with the field size and the grid sizes used for numerical computations. It can be shown [11], that the tubular well with a small diameter is approached, in its “limit shape”, as a line source. 1 γ (u) is a vector function with values in R3 defined over the interval a u b. It must fulfill certain conditions of differentiability in order to assure proper curve modeling and the existence of the trajectory torsion [14,30].
436
R.J. González-Requena, J.M. Guevara-Jordan / Applied Numerical Mathematics 40 (2002) 433–449
From the physical point of view, the integral in the right side of (2a) has been used to represent the density of a continuous and longitudinal source [35] and it can be interpreted as the limit of point densities discretely distributed along the curve γ when this is partitioned without any limit [41,42]. On the other hand, only the total production rate of the well, Q(t), is known for general applications. This type of condition for the production rate is known as the uniform flow condition [21] and it is obviously related to the production rate per length unit and per unit of time q through the formula [35] Q(t) −
b
q γ (u), t γ (u) du.
(3)
a
3. Source term solution A very concise mathematical framework is exposed here in order to second the use of generalized functions (distributions) when modeling source terms in differential equations, and their corresponding source term solutions. It is assumed that most of these concepts are known by the reader. From the mathematical viewpoint, the integral in the right side of (2a), or more briefly denoted by the integral
q(γ , t)δ(x − γ ) dγ ,
γ
symbolizes the generalized function that operates under the rule [43]
(γ , t)δ(x − γ ) dγ , ϕ(x) =
γ
q(γ , t)ϕ(γ ) dγ
(4)
γ
for all test or “good” functions ϕ(x) which, in essence, are infinitely differentiable functions on Rn with compact support. Now, the equation ∂p(x, t) − p(x, t) = q(γ , t)δ(x − γ ) dγ (5) ∂t γ
can be interpreted as a differential equation in a distributional sense and only requires the generalized functions ∂p − p and (γ , t)δ(x − γ ) dγ ∂t γ
to be equal in the reservoir domain Ω for each t > 0 [25,35,43]. This distributional approach allows a natural interpretation of the physical situation and the use of the Fourier transform to obtain the solution of Eq. (5). The solution of Eq. (5) will be referred to as a source term solution, without consideration of boundary conditions. The expression fundamental solution will be used to denote the function traditionally known as free Green’s function. Using Fourier transform, Vladimirov [41,42] showed that a fundamental solution of the diffusivity operator is given by (3D case)
γ − x2 1 exp − . FS(x, t; γ , τ ) = (4π(t − τ ))3/2 4(t − τ )
(6)
R.J. González-Requena, J.M. Guevara-Jordan / Applied Numerical Mathematics 40 (2002) 433–449
437
On the other hand, the Fourier transform of the function f (x) is defined by the integral (provided it exists for every α) FTx [f ] ≡
···
f (x) exp −i(α ◦ x) dx1 · · · dxn = F (α).
(7)
Rn
The subscript x indicates that the Fourier transform has been applied to the variable x in Rn . The Fourier transforms of test functions are also test functions, which are adequately used to define the Fourier transforms of generalized functions [25,41–43]. In addition, if g(x, γ ) is a generalized function of x for all γ along the curve γ , it can be shown [25] that
g(x, γ ) dγ =
FTx γ
FTx g(x, γ ) dγ .
(8)
γ
To solve Eq. (5), we apply the Fourier transform to both sides and use expression (8). This generates an easy differential equation on t which solution is given by P (α, t) =
t
exp −α2 (t − τ ) − i(α ◦ γ ) q(γ , τ ) dγ dτ.
(9)
0 γ
This last expression is really the Fourier transform of the solution p of Eq. (5), so the solution p can be obtained using the inverse Fourier transform and is given by p(x, t) ≡ ST(x, t) =
t 0 γ
γ − x2 q(γ , τ ) exp − dγ dτ. (4π(t − τ ))3/2 4(t − τ )
(10)
4. Singularity programming The solution of (2) for a regular domain such as a parallelepiped or a rectangle is easily found through analytical methods [19]. When the frontier Ω is arbitrarily shaped, it is inevitable to resort to numerical methods, which, in general, are unable to adequately handle the Dirac delta-type function involved in the diffusivity [27]. According to the method here exposed, the solution p is considered as the addition of the source term solution ST, associated to the generalized function, and a smooth function, Preg , that can be solved through conventional numerical methods: p = Preg + ST.
(11)
It can be shown that Preg satisfies the following IBV problem ∂Rreg − Preg = 0 ∂t ∂ST ∂Preg (x, t) = − (x, t, γ ) ∂n ∂n Preg (x, 0) = 0,
in Ω × (0, ∞),
(12a)
in ∂Ω × (0, ∞),
(12b)
in Ω.
(12c)
The major advantage of (12) is that it represents a very conventional IBV problem for the diffusivity equation. It can be solved through any numerical method selected, or, given the case, by means of a
438
R.J. González-Requena, J.M. Guevara-Jordan / Applied Numerical Mathematics 40 (2002) 433–449
proper analytical method or a combination of both. In the present work, the Finite Element Method (FEM) was used to solve the reduced problem given by (12). Once Preg has been calculated, the pressure p solution of (2) is uniquely determined by (11). Lastly it is worth pointing out that using this method it possible to accept a trajectory for the horizontal well as curved and/or as twisted as can be mathematically handled.
5. Wellbore condition features The most important parameters used for the system of wells in reservoir simulation are flow rate and wellbore pressure. Consequently, two different conditions for the wellbore can be considered to address the problem under study (2). First, the uniform flux condition that assumes that the total production rate of the well is known, implying that the factor q(γ (u), t) in (10) is a known constant. Although the production rate is probably the best known information about oil wells, these kind of studies must be completed by including the capability of handling the wellbore pressure condition [3–6,24]. Here, the well production is conceived under the condition of uniform wellbore pressure, implying that the strength q(γ (u), t) of the production rate of the well is unknown. This pressure is defined as the reservoir pressure measured at the well radius and it is used in all cases to compute the unknown production rate of the well. As mentioned before, the proper solution of the reduced problem and the contribution of the corresponding source term solution give the required solution. However, the source term solution also deserves numerical treatment. For computational purposes, the spatial and time integrations involved in the source term solution (10) are discretely approximated. The following expression is then established:
ST x, t; γ (ui ), τj =
NU NT
q γ (ui ), τj C x, t; γ (ui ), τj γ (ui )ui τj .
(13)
j =1 i=1 U Here, {ui }N i=1 is a convenient partition of the interval [a, b]. Parameter u takes its values within this interval for generating each point of the well trajectory and is associated with the chosen numerical T method of integration. The set of times {τj }N j =1 is the corresponding subdivision of [0, t] linked to the selected method of numerical integration, with a last time step τNT equal (or very close) to t. Finally, C(x, t; γ (ui ), τj ) denotes the fundamental solution (6) evaluated at (γ (ui ), τj ). The discrete form (13) is used in Eq. (11) to obtain the numerical approximation of the pressure solution at any reservoir location x and at any instant t. Evidently, a corresponding numerical approximation is also used to treat the normal derivative of the source term that is present in the boundary condition of the reduced problem (12). When the uniform flux condition is imposed, the numerical solution given by (13) can be directly calculated due to the knowledge of q(γ (ui ), τj ). When the wellbore pressure condition is used, the strength q(γ (u), t) of the production rate of the well is unknown and must be previously established in order to find the solution using (13). Consequently, a different scheme must be proposed. NU Let {x w n }n=1 denote the wellbore positions (at the well radius) associated with the well trajectory, then, the necessary system of NU equations is given by
w pwf x w n , t = Preg x n , t +
NU NT
q γ (ui ), τj C x w n , t; γ (ui ), τ γ (ui ) ui τj
j =1 i=1
(14)
R.J. González-Requena, J.M. Guevara-Jordan / Applied Numerical Mathematics 40 (2002) 433–449
439
w for n = 1, 2, . . . , NU ; at each instant t, and where pwf (x w n , t) denotes the wellbore pressure at x n . Since the production rates q(γ (ui ), τj ) are known for j = 1, 2, . . . , Nτ − 1, the double summation associated with these indexes is called Γ (x w n , τNT −1 ) and the pressure differential is denoted by w w , t) = p (x , t) − p (x , t). Under these new variables, the system of equations (14) can be p(x w wf reg n n n written as NU
w w q γ (ui ), tNT C x w n , τNT ; γ (ui ), τNT −1 γ (ui ) ui τNT = p x n , t − Γ x n , τNT −1
(15)
i=1
for n = 1, 2, . . . , NU and at each instant t. These equations represent a generalized wellbore model for the horizontal well γ , because they relate the production rates with the pressure differential along the well. If wellbore pressures are given as input data, then this model is a linear system for the rates q(γ (ui ), τNT ). Its solution provides the production rates needed by the diffusivity equation solver in order to compute the pressure at the following time step. In this way, the extension of the new method to handle the wellbore pressure condition is achieved. In addition, the wellbore model can be used in reverse direction allowing the calculation of the wellbore pressure based on the rates. The wellbore model developed in this section is explicit. It means that the flow rate is determined before invoking the diffusivity solver on each time step. This kind of model is very simple but it may pose stability problems if long time steps are taken or strong changes in pressure and rate are observed. In these cases, shorter times steps could avoid this problem. If the times steps needed are too small then a more complex or implicit wellbore model is required. However, such development is not needed to extend and validate the numerical method presented in this article. More complex models for wells can be always developed following similar arguments to those presented in this section. Additional considerations will be presented to end this section. The discrete form of parameter u U induces a corresponding subdivision of points, {γ (ui )}N i=1 , on the well trajectory. These points can be considered as “gravity centers” of small arc elements composing the whole well trajectory. Moreover, the arc elements are here conceived as arc sources with respective wellbore pressures (concentrated at points γ (ui )) and production rates (given by q(γ (ui ), t)). This interpretation can be compared with the concept of Discrete Flux Elements (DFE) introduced by Azar-Nejad et al. [3–6] for modeling the potential (pressure) distribution generated by 2D finite sources: straight lines and curved lines. However, when these authors derive their definitions, they do it from results obtained for cases with straight horizontal wells parallel to coordinate axes, and in rectilinear reservoirs. Azar-Nejad et al. [3,5] assure that the source configuration must be straight to do the integral over the wellbore length. However, integrating over straight sources parallel to a coordinate axe implies that factor γ (u) in Eq. (10) is constant. Consequently, when the DFE concept is extrapolated to wells with arbitrary geometry (not twisted), the factor γ (u) is omitted, and it should be present in those integrals representing the curvilinear line sources. This length factor does not necessarily have to be constant for curved and twisted trajectories where, in fact, integration can also be done under adequate postulates of mathematical modeling. As has been seen, the factor γ (u) depends on the trajectory shape, and it has a different value at each point γ (u). Its absence could distort the accuracy of results, particularly when curved and twisted well trajectories are modeled.
440
R.J. González-Requena, J.M. Guevara-Jordan / Applied Numerical Mathematics 40 (2002) 433–449
6. Results and analysis 6.1. Finite element method The finite element method (FEM) [27,39] was used as the chosen numerical technique for solving the diffusivity equation in the numerical implementation of the SPM addressed by this paper. The PDE2D finite element software [39] was used to this end. It allows the solution of general systems of partial differential equations using polynomial bases of up to fourth order into 2D domains of arbitrary shape. FEM was selected because it is better to handle boundary conditions, since this is a major factor in implementing the SPM. Moreover, the FEM could be accommodated to tortuous (twisted) trajectories of horizontal wells using large and complex grid refinements [29], but the computational cost inherent to such procedure could render it unusable, particularly in 3D runs, unsteady flow cases, and domains with an arbitrary boundary. The polynomial base in all 2D examples was restricted to be of first order, which, under very general conditions [18], is equivalent to a scheme of finite differences. In 3D cases, the PDE2D software can only allow polynomial bases of third order into orthoedric-shaped boxes. This reflects the limited number of allowed blocks for studying 3D problems mostly due to restrictions in the available hardware (Sun Sparc II). The discretization caused in any dimension by PDE2D is completely implicit, preventing restrictions to the length of time steps based on the Courant–Friedrich–Lewy condition [27]. A more detailed FEM exposition lies beyond the purposes of this paper, as it focuses on showing the potential and flexibility of this approach for the tortuous well modeling, and to solving the selected cases. 6.2. Examples Four specific examples are presented below to validate the mathematical model, as well as the SPM used for the simulation of the pressure distribution generated by horizontal wells. A trapezoidal rule was used for the numerical integrals involved in these examples. 6.2.1. Example 1 (Validation against an analytic solution) For validation purposes, a rectangular reservoir given by Ω = [0, 2] × [0, 1] was selected. It contains a straight horizontal well γ of unit length parallel to the x axis and its center coincides with that of Ω, i.e., the coordinates (1, 0.5). An unitary uniform flux was assumed. Under these conditions, an analytical solution [19] and several numerical solutions obtained with two different methodologies were developed for solving the IBV problem (2). The numerical methodologies are the direct application of the FEM 2 to an equivalent problem, and the Singular Programming Methodology (SPM) proposed in this article, which also includes the use of the FEM. When using the direct application of the FEM, an approximation to the distributional representation of the horizontal well is considered, and the original differential equation (2a) becomes the following associated equation ∂pa − pa = ∂t
256, 0,
1 if (x, y) ∈ [0.5, 1.5] × 0.5 − 512 , 0.5 + otherwise.
1 512
,
2 The methodology that applies the FEM directly without removing singularities will be denoted by dFEM.
(16)
R.J. González-Requena, J.M. Guevara-Jordan / Applied Numerical Mathematics 40 (2002) 433–449
441
In this equation, the generalized function has been approximated by the indicator function shown at the right side of (16). It is well known [11] that these kinds of functions can be used at the limit process to approximate the singular function associated to the horizontal well representation. The validity of this function was verified by means of grid refinements and comparison with the analytical solution. Special grids were fitted to take into account the presence of the horizontal well, because the integration formulas can bypass the small band representing the well on uniform grids. The best approximation to the analytical solution using dFEM was found on a grid of, approximately, five thousand elements, at the dimensionless times 0.001, 0.01, 0.1 and 1. SPM was implemented using the same grids designed for dFEM. Although SPM does not need this kind of grid—it would work fine on global uniform grids—it was used just to match the dFEM working conditions and facilitate the comparison of the corresponding solutions. The best solution for the SPM was produced on a grid of, approximately, four hundred elements when it was compared with the analytical solution at the same times considered for dFEM. To obtain the convergence rate of the numerical methods in this test problem, the L2 (Ω) error of the difference between the numerical and analytical solutions was determined at the dimensionless time t −1. The pseudo-steady condition has been achieved at this moment, so no further changes in the solution are expected. Table 1 summarizes the errors for both cases and the grid dimension. Table 1 Errors in L2 (Ω) norm SPM Log(h)
Loganalytic − numerical
dFEM Grid
Log(h)
Loganalytic − numerical
Grid
−1.0986
−8.0102
72
−2.3979
−9.0612
968
−1.6094
−9.2567
200
−2.7081
−9.3093
1800
−1.9459
−9.6800
392
−3.2189
−9.7180
5000
Fig. 1. Example 1: L2 error vs. grid length. Logarithm values.
442
R.J. González-Requena, J.M. Guevara-Jordan / Applied Numerical Mathematics 40 (2002) 433–449
Table 2 CPU times for SMP and dFEM Number of points where solution is required
Normalized CPU times: time_SPM/time_dFEM
10 points
0.60
50 points
0.62
150 points
0.74
315 points
0.98
The contents of the table are shown in Fig. 1, where the slopes of the straight lines allow to conjecture the order of convergence for each method. The magnitude h represents the maximum length of all elements in the grid. Notice that the dFEM and the SPM, respectively, present almost first and second order of convergence. For each method, the CPU times were computed on a Sun workstation (Sparc II) on the respective finer grid. CPU times were found to be almost equal for both methods, in spite of the fact that dFEM was applied in a grid twelve times bigger than the one used in the SPM. However, if the solution is needed in a reduced number of points then SPM may be twice as fast as FEM as observed in Table 2. In any case, CPU times should be carefully analyzed because this is really a comparison between an efficient commercial code [39], and a domestic mixed algorithm developed for research purposes. 6.2.2. Example 2 (Applying the SPM to 2D irregular reservoirs) The SPM will be used in this example to numerically simulate the pressure behavior created by a horizontal well in a reservoir of irregular boundary. To this end, region Ω in Fig. 2(a) was considered with the same initial and boundary conditions given in (2). As in the previous example, well length is unity, and its total rate production is of magnitude one. Following the same procedure as before, we will also solve the associated equation to (2) where the generalized function is approximated by an indicator function. The solution of the associated problem will be obtained by directly using the FEM. This associated equation is given by
1 1 ∂pa 256 if (x, y) ∈ [0.5, 1.5] × 0.602 − , 0.5 + , 512 512 − pa = (17) ∂t 0 otherwise. Here, the well center is at the point (1, 0.602). Figs. 2(a) and (b) show the grids used for the application of SPM and dFEM, respectively. Both solutions converge to the correct solution. In this case, a numerical method has been assumed to converge on a certain grid if the difference between numerical solutions obtained in successive refinements are not punctually bigger than 10e−6. Again, it is worth mentioning that this kind of grids, which takes into account the presence/location of the horizontal well, are required only when using the FEM to solve the associated problem (17), but they are unnecessary with the application of the SPM for solving (2). This feature will be highlighted in the 3D example. The following figures represent a sequence of the pressure distribution obtained by application of SPM. Fig. 3(a) shows the drainage area generated by the horizontal well at dimensionless time t = 0.01. It can be observed in this case that pressure close to the well behaves as if the reservoir was unlimited. Additionally, the elliptical pressure isolines unambiguously define the geometric location of the well in the reservoir. Fig. 3(b) shows the final drainage area, which reflects pressure distribution at dimensionless
R.J. González-Requena, J.M. Guevara-Jordan / Applied Numerical Mathematics 40 (2002) 433–449
443
(a)
(b) Fig. 2. Example 2: (a) Coarse grid system for numerical solutions by SPM. (b) Refined grid system for numerical solutions by dFEM.
time t = 1.0 when the pseudo-stationary condition is reached. At this time, the drainage area is the whole reservoir and it has thus ceased to expand. Also, it can be observed that the contour lines close to the well are not concentric ellipses around the horizontal well, but rather around its center, which is a consequence of the adopted modeling condition, i.e., uniform production rate. For this constant flux case, a pressure gradient along the well is present because the flow has to enter the well uniformly along its length [9]. Such pressure behavior is also typical for fractured vertical wells, where fluid extracted through the fracture tends to move towards the center, onto the position where the vertical well is supposed to be pumping fluid to the surface. This is the reason why many of the available results for fractured wells are frequently applicable to horizontal wells. It is worth noting that fluid is not always pumped out of horizontal wells through their centers, but rather through one of its ends or using producing sub-intervals. Clearly, these situations require different sets of hypotheses for modeling purposes, but they can be also addressed by the SPM. In conclusion, this example shows the SPM capability to simulate the unsteady behavior of the pressure generated by a horizontal well in an arbitrarily shaped 2D reservoir. In addition, using the SPM in less
444
R.J. González-Requena, J.M. Guevara-Jordan / Applied Numerical Mathematics 40 (2002) 433–449
(a)
(b) Fig. 3. Example 2: (a) Pressure isolines around the horizontal well at time t = 0.01. (b) Isopressure lines around the horizontal well at time t = 1.0.
refined grids solutions were obtained which were comparable to those obtained with the dFEM for the equivalent problem in a more refined grid. 6.2.3. Example 3 (Applying the SPM to a well producing by wellbore pressure in 2D irregular reservoirs) This example gives us the opportunity of applying the wellbore model developed here to include wellbore pressure condition as input data. The reservoir and horizontal well are the same used in Example 2. Indeed, all the conditions of that example will be the same up to the simulation time t = 1. At this point the wellbore model is used to keep the wellbore pressure along the horizontal well at −0.51 from instant t = 1 to the time t = 1.5. This restriction on the boundary condition produces a completely different pressure behavior around the well. The pressure field at time t = 1.5 is presented in Fig. 4(a). It shows that the isopressure lines do not intercept the horizontal well because all the wellbore pressures, at this instant, have the value −0.51 imposed as a restriction.
R.J. González-Requena, J.M. Guevara-Jordan / Applied Numerical Mathematics 40 (2002) 433–449
445
(a)
(b) Fig. 4. Example 3: (a) Pressure distribution around the horizontal well at time t = 1.5. (b) Effects of wellbore model on pressure and rate.
In addition, Fig. 4(b) shows the effects of the wellbore model on the rate and pressure evaluated at the center of the horizontal well. It can be observed that the flow rate is constant until the time t = 1; after this moment the well rate is reduced in time by the wellbore model in order to keep the wellbore pressure constant at −0.51. The small oscillations of pressure around this value are a consequence of the explicit coupling of the wellbore model. All these observations of the behavior of pressure and production for the horizontal well are consistent with the physics of the problem at hand, so it may be concluded that the wellbore model developed for the SPM gives a satisfactory extension of this method. This extension depends neither on the reservoir boundary shape nor on a tortuously curved/twisted well trajectory. 6.2.4. Example 4 (Applying the SPM to wells with curved and twisted trajectory) The objective of this relevant case is to establish the unsteady pressure distribution generated by a horizontal well whose trajectory is as curved and twisted as a circular helix. The problem to be solved is
446
R.J. González-Requena, J.M. Guevara-Jordan / Applied Numerical Mathematics 40 (2002) 433–449
the same as expressed in (2). Here, γ (u) is the analytical representation of this tortuous trajectory [20], which is obtained by rotating the helix of equation (0.8u, 0.5 + 0.25 cos(π u), 0.5 + 0.25 sin(π u)), until its own axis is aligned with the direction of one of the major diagonals of the domain Ω = [0, 2] × [0, 1] × [0, 1]. Parameter u takes values between a = 0 and b = 2.5. Figs. 5(a) and (b) show a transient sequence for the pressure response generated by this twisted well at times 0.001 and 0.1, respectively. Under the flow scheme known as early radial, Fig. 5(a) shows pressure distribution (isosurfaces) honoring the curved and twisted shape of the helix-shaped trajectory of this well. This is similar to the elliptic pressure isolines that make up the elliptic drainage area reserved for
(a)
(b) Fig. 5. Example 4: (a) Pressure isosurfaces around a tortuous well at dimensionless time t = 0.001. (b) Pressure isosurfaces around a tortuous well at dimensionless time t = 0.1.
R.J. González-Requena, J.M. Guevara-Jordan / Applied Numerical Mathematics 40 (2002) 433–449
447
these times when 2D studies of straight-well response are performed. This can be also seen in Example 2 and in [3–6]. Fig. 5(b) presents the shape of the pressure isosurfaces resulting from the spherical flow at both ends of the well. Equivalent results have been reported by Azar-Nejad et al. [6] in their research with straight horizontal wells in a 2D rectangular reservoir. This application illustrates both the potential for emulating well trajectories as curved and twisted as the relevant mathematical representation can handle, as well as the total independence between the geometry/location of the well and the proposed grid for obtaining the numerical solution of the reduced problem.
7. Conclusions In this paper, the SPM was introduced and validated like a new methodology for assessing the transient pressure distribution generated by tortuous horizontal wells in arbitrary shaped reservoirs. In the 2D examples considered, the SPM presented excellent results even on coarse grids, and it was numerically established that its order of convergence is quadratic. Additionally, the SPM proved to have faster numerical convergence than the convergence reached by the direct application of the FEM to equivalent problems. Where possible, solutions using the SPM in coarse grids compared favorably with those obtained applying the dFEM in more refined grid. This evidenced the superiority of this methodology. The SPM applies equally well to wells that are producing either by constant rate of production or by uniform wellbore pressure. The SPM was shown to have potential for emulating well trajectories as curved and twisted as the relevant mathematical representation can treat. Likewise, this methodology has shown that arbitrarily shaped reservoirs are not a limitation. Summarizing, the main SPM advantages are the convergence in coarse grids; its elegant capability for emulating horizontal well trajectories as curved and twisted as required, and the total independence between numerical grids and the geometry/location of the trajectory well.
References [1] W.F. Ames, Numerical Methods for Partial Differential Equations, Academic Press, New York, 1977. [2] C. Astudillo, E. Guedez, Tratamiento de singularidades en el modelaje de pozos de un yacimiento, petrolífero, Graduate Thesis, Mathematics Department, School of Sciences, Universidad Central de Venezuela, Caracas, Venezuela, 1993. [3] F. Azar-Nejad, W. Tortike, S. Farouq Ali, Potential distribution around the sources with finite length (horizontal, vertical partially penetrating wells and fractures). Part II: Transient flow, Paper 35269, Society of Petroleum Engineers, 1996. [4] F. Azar-Nejad, W. Tortike, M. Farouq Ali, Potential distribution around the sources with finite length (horizontal, vertical partially penetrating wells and fractures). Part I: Steady state fluid flow, Paper 35270, Society of Petroleum Engineers, 1996. [5] F. Azar-Nejad, W. Tortike, S. Farouq Ali, 3D analytical solution to the diffusivity equation for finite sources with applications to horizontal wells, Paper 35512, Society of Petroleum Engineers, 1996.
448
R.J. González-Requena, J.M. Guevara-Jordan / Applied Numerical Mathematics 40 (2002) 433–449
[6] F. Azar-Nejad, W. Tortike, S. Farouq Ali, Performance of horizontal wells with irregular geometry, Paper 36550, Society of Petroleum Engineers, 1996. [7] D.K. Babu, A. Odeh, Productivity of horizontal wells, Paper 18298, Society of Petroleum Engineers, 1985. [8] D.K. Babu, A. Odeh, Transient flow horizontal well, pressure drawdown, and buildup analysis, Paper 18802, Society of Petroleum Engineers, 1989. [9] R.M. Butler, Horizontal Wells for the Recovery of Oil, Gas and Bitumen, Monograph, Vol. 2, The Petroleum Society of the Canadian Institute of Mining, Metallurgy and Petroleum, Calgary, AB, 1994. [10] H.S. Carslaw, J. Jaeger, Conduction of Heat in Solids, Oxford University Press, Oxford, 1959. [11] G. Chavent, J. Jaffre, Mathematical models and finite elements for reservoir simulation, in: J.L. Lions, G. Papanicolaou, H.F. Fujita, H.B. Keller (Eds.), Studies in Mathematics and its Applications, Vol. 17, Elsevier, Amsterdam, 1986. [12] B.C. Craft, M. Hawkins, Applied Petroleum Reservoir Engineering, Prentice-Hall, Englewood Cliffs, NJ, 1957. [13] L.P. Dake, Fundamentals of Reservoir Engineering, Elsevier, Amsterdam, 1978. [14] M.P. Do Carmo, Differential Geometry of Curves and Surfaces, Prentice-Hall, Englewood Cliffs, NJ, 1976. [15] J. Douglas, R. Ewing, M. Wheeler, Approximation of the pressure by a mixed method in the simulation of miscible displacement, RAIRO Analyse Numerique 17, 1983. [16] P. Duchateau, D.W. Zachmann, Partial Differential Equations, McGraw-Hill, New York, 1988. [17] R. Ewing, Problems arising in the modeling of processes for hydrocarbon recovery, in: R.E. Ewing (Ed.), The Mathematics of Reservoir Simulation, SIAM, Philadelphia, PA, 1983, Chapter 1. [18] R. Ewing, M.F. Wheeler, Galerkin Methods for Miscible Displacement Problems with Point Sources and Sinks-Unit Mobility Ratio, Mathematical Methods for Energy Research, SIAM, Philadelphia, PA, 1983. [19] A. Fermin, Método de Green para la simulación numérica de pozos horizontales, Graduate Thesis, Mathematics Department, School of Sciences, Universidad Central de Venezuela, Caracas, Venezuela, 1998. [20] R. González, Metodologías alternas para la solución numérica de la ecuación de difusión bajo condiciones de borde pertinentes a las aplicaciones en simulación de yacimientos, Ph.D. Thesis, Mathematics Department, School of Sciences, Universidad Central de Venezuela, 2000. [21] A.C. Gringarten, H.J. Ramey, The use of source and Green functions in solving unsteady flow problems in reservoirs, Paper 3818, Society of and Petroleum Engineers, 1973. [22] J. Guevara, A. Fermín, R. González, A new approach for horizontal wells singularities in petroleum engineering, Paper 51924, Society of Petroleum Engineers, 1999. [23] L. Hayes, R. Kendall, M. Wheeler, Treatment of source and sinks in steady-state reservoir engineering simulations, in: R. Vichnevetsky (Ed.), Advances in Computer Methods for Partial Differential Equations II, IMACS (AICA), 1977. [24] J.K. Jasti, V.R. Penmatcha, D.K. Babu, Use of analytical solutions in improvement of simulator accuracy, Paper 38887, Society of Petroleum Engineers, 1997. [25] D.S. Jones, Generalized Functions, McGraw-Hill, New York, 1966. [26] L.S. Koh, D. Tiab, A boundary element algorithm for modeling 3D horizontal well problems using 2D grids, Paper 26228, Society of Petroleum Engineers, 1993. [27] L. Lapidus, G. Finder, Numerical Solution of Partial Differential Equations in Sciences and Engineering, Wiley, New York, 1988. [28] J. Liggett, P. Liu, Boundary Integral Equation Method for Porous Media Flow, George Allen & Unwin, London, 1983. [29] N. Morita, Transient finite element code: a versatile tool for well performance analysis, Paper 25878, Society of Petroleum Engineers, 1993. [30] B. O’Neill, Elementary Differential Geometry, Academic Press, New York, 1976. [31] E. Ozkan, R. Raghavan, New solution for well test analysis problems: Part 1—analytical considerations, Paper 18615, Society of Petroleum Engineers, 1991.
R.J. González-Requena, J.M. Guevara-Jordan / Applied Numerical Mathematics 40 (2002) 433–449
449
[32] E. Ozkan, R. Raghavan, New solution for well test analysis problems: Part 2—computational considerations and applications, Paper 18616, Society of Petroleum Engineers, 1991. [33] E. Ozkan, R. Raghavan, New solution for well test analysis problems: Part 3—additional algorithms, Paper 23593, Society of Petroleum Engineers, 1994. [34] R. Raghavan, Well Test Analysis, Penn Well, Tulsa, OK, 1993. [35] R. Raghavan, E. Ozkan, A Method for Computing Unsteady Flows in Porous Media, Longman Scientific and Technical, Harlow, UK, 1994. [36] W.A. Roper, J. Merchant, C. Duvall, Combination of numerical and analytical techniques to improve waterflood model efficiency, Paper 2031, Society of Petroleum Engineers, 1962. [37] A.J. Rosa, R. Carvalho, A mathematical model for pressure evaluation in an infinite conductivity horizontal well, Paper 15967, Society of Petroleum Engineers, 1989. [38] T.F. Russell, M.F. Wheeler, Finite elements and finite differences methods for continuous flow in porous media, in: R.E. Ewing (Ed.), The Mathematics of Reservoir Simulation, SIAM, Philadelphia, PA, 1983, Chapter 2. [39] G.V. Sewell, Analysis of Finite Element Method: PDE2D/Protran, Springer, Berlin, 1985. [40] V. Villamizar, Personal communication, Mathematics Department, School of Sciences, Universidad Central de Venezuela, 1998. [41] V.S. Vladimirov, Equations of Mathematical Physics, Marcel Dekker, New York, NY, 1971. [42] V.S. Vladimirov, Generalized Functions in Mathematical Physics, Mir Publishers, Moscow, 1979. [43] A.H. Zemanian, Distribution Theory and Transform Analysis, Dover, New York, 1987.
New approach for describing transient pressure response generated by horizontal wells of arbitrary geometry Reinaldo J. González-Requena a,∗ , Juan M. Guevara-Jordan b a Escuela de Ingeniería de Petróleo, Facultad de Ingeniería, Universidad Central de Venezuela, Caracas, Venezuela b Escuela de Matemáticas, Facultad de Ciencias, Universidad Central de Venezuela, Caracas, Venezuela
Abstract This paper is aimed at developing a methodology for studying the transient pressure behavior of horizontal wells with any curvilinear trajectory in an isotropic/anisotropic arbitrarily shaped reservoir. This methodology employs generalized functions to represent the tortuous horizontal well. A particular way of removing the singularities involved in the partial differential equation is based on reducing the original problem to the conventional solution of the homogeneous diffusivity equation under any given initial and boundary conditions. The Green function method and any standard numerical technique are combined in a single computational strategy to obtain the transient pressure response generated by a curved and twisted horizontal well in reservoirs with irregular boundaries. Analytical methods can be also used, whenever possible, to solve the reduced problem. This proposal can be easily broadened to analyze the performance of the pressure transient of any kind of reservoir sources or sinks that can be modeled using generalized functions. Some models are presented. 2001 Published by Elsevier Science B.V. on behalf of IMACS. Keywords: Horizontal well; Sources; Mathematical modeling; Finite elements method; Diffusivity equation; Transient pressure; Generalized functions; Curved/twisted trajectory
1. Introduction The number of horizontal wells has increased considerably since the early 1980s within the petroleum industry. Over the past twenty years, the advances in drilling and production techniques have paved the way to the development, study, and adoption of new models for representing horizontal wells which are more appropriate for the analysis of the well performance. This work presents a mathematical model of the transient pressure response generated by a horizontal and/or tortuous well which explicitly considers the modeled well as a generalized function in the * Corresponding author.
E-mail address: [email protected] (R.J. González-Requena). 0168-9274/01/$22.00 2001 Published by Elsevier Science B.V. on behalf of IMACS. PII: S 0 1 6 8 - 9 2 7 4 ( 0 1 ) 0 0 0 9 4 - 0
434
R.J. González-Requena, J.M. Guevara-Jordan / Applied Numerical Mathematics 40 (2002) 433–449
non-homogeneous diffusivity equation [2,19,35]. This modeling procedure is not naturalsince a well should be represented as a boundary condition [12,13,35]. For instance, when fluids are subtracted from reservoirs using production wells, they are characterized by intensive convergent flow patterns and pressure field singularities. This has led to the representation of wells by means of generalized functions or distributions. One of the purposes of the present paper is the explicit use of this type of functions in the diffusivity equation when representing a well of arbitrary trajectory. Evidence of the implicit use of such generalized functions for representing vertical, horizontal, and fractured wells can be found in the petroleum literature. For instance, the pressure distribution solutions in infinite and rectangular reservoirs can be seen in [3–8,21,31–34,37]. Likewise, Ligget and Liu [28] have used conveniently redefined Dirac delta functions to represent vertical wells, at the expenses of certain level of rigorousness in the generalized function model. On the other hand, Raghavan and Ozkan [35] provide solutions to many reservoir cases (in rectangular coordinates) based on a similar distributional approach. A specially designed numerical method is necessary to solve the well model developed in the case of 2D or 3D reservoirs with arbitrary boundaries due to the inherent difficulties involved in representing generalized functions using traditional numerical methods like finite differences or finite elements [27]. Therefore we hereby propose an algorithm that combines the removal of singularities (singularity programming), Green functions, and the finite elements method in a single computational strategy to numerically solve the Initial and Boundary Value (IBV) problem. This numerical method, hereafter called the Singularity Programming Method (SPM), is not new in its applicability. Different authors have used it to solve elliptic equations [1,15,18,23,36,38]. For the purposes of the diffusivity equation, the singularity removal method is not as direct due to the additional time nature of the Dirac delta function used in the partial differential equation, and therefore it must be reformulated. Astudillo and Guédez [2] originally applied this reformulated methodology in the case of vertical wells in rectangular domains. They used the finite differences method as a numerical approach. In this sense, one of the major contributions of the present paper is to extend and validate the applicability of SPM to the mathematical model proposed for horizontal or tortuous wells in arbitrarily shaped reservoirs. The authors can assure that this extension of the SPM has not been previously recorded in the specialized literature. The pressure distribution created by a “fictitious” horizontal well with a trajectory given by a circular helix will be ascertained in this paper. Such potential use is another contribution of this methodology. So far the cases of the most tortuous trajectory well found by the authors in the specialized literature are that of a curved, torsionless horizontal well proposed by Koh and Tiab [26], and those, also torsionless, wells, reviewed by Azar-Nejad et al. [3]. In addition to its potential for really addressing arbitrary well geometries, this methodology can be extended to represent other reservoir anomalies such as fractures. Section 2 presents some basic mathematical and physical considerations necessary to use generalized functions when modeling the pressure response of horizontal wells. Section 3 considers the fundamental solution of the diffusivity operator together with the associated source term solutions. They will be used to program the singularity removal algorithm in Section 4. Section 5, probably the most important, presents numerical issues used to account for the usual wellbore conditions. Section 6 focuses on examples and their analysis. Section 7 presents the conclusions and final considerations.
R.J. González-Requena, J.M. Guevara-Jordan / Applied Numerical Mathematics 40 (2002) 433–449
435
2. Mathematical model According to the standard considerations [3] about rock/fluid conditions, and taking its dimensionless form, the basic equation describing the continuous pressure distribution generated by a well of arbitrary trajectory γ can be simplified as follows [19,35,40]: ∂p − p = ∂t
b
q γ (u), t δ x − γ (u) γ (u) du
(1)
a
where p is pressure, t is time, q is production rate at the wellbore per length unit and per unit of time [25], δ is the Dirac delta function [25,41–43], and the line integral is calculated along the curve γ used for modeling the well trajectory. The analytical representation of this curve is γ (u) = (γ1 (u), γ2 (u), γ3 (u)). Here, γ (u) is the Euclidean norm of the velocity vector, which is tangent to the curve at γ (u). 1 It is worth pointing out that for the dimensioning of (1), the gravity term was considered through the Hubbert potential [3–6]. For simplicity’s sake, the symbols used in Eq. (1) are also used to represent dimensionless variables, such as time, space and Hubbert potential. In rigor, p is the dimensionless form of the potential, but from this point on we will refer to it just as pressure. Modeling and simulating the behavior of pressure generated by a tortuous well in a reservoir Ω requires Eq. (1) to be solved subject to certain conditions at the reservoir boundary ∂Ω with a given initial condition. For the purposes of the present paper, we shall take into consideration only the absence of flow in ∂Ω, i.e., a zero Neumann condition. Nevertheless, the method is perfectly applicable to the cases of constant pressure condition (dirichlet condition) or mixed conditions. Therefore, the following problem with initial and boundary value (IBV) is assessed: ∂p − p = ∂t
b
q γ (u), t δ x − γ (u) γ (u) du
in Ω × (0, ∞),
(2a)
in ∂Ω × (0, ∞),
(2b)
in Ω,
(2c)
a
∂p (x, t) = 0, ∂n p(x, 0) = 0,
where the semi-infinite interval is associated to time and the zero initial condition is consistent with the dimensioning used for these studies. The source term in Eqs. (1) and (2a) is an idealization of the real condition of the wellbore, which can be proposed as an inner boundary condition. In this approach, the well, Dε , is conceived as a tube of diameter ε around the axis trajectory γ , and the boundary condition equals the flow rate of the fluid through the tube surface ∂Dε with the known production rate of the well. This representation of the well is very close to the physical reality, but is not used in reservoir simulation. Engineers prefer to neglect the well tubing because its diameter is very small in comparison with the field size and the grid sizes used for numerical computations. It can be shown [11], that the tubular well with a small diameter is approached, in its “limit shape”, as a line source. 1 γ (u) is a vector function with values in R3 defined over the interval a u b. It must fulfill certain conditions of differentiability in order to assure proper curve modeling and the existence of the trajectory torsion [14,30].
436
R.J. González-Requena, J.M. Guevara-Jordan / Applied Numerical Mathematics 40 (2002) 433–449
From the physical point of view, the integral in the right side of (2a) has been used to represent the density of a continuous and longitudinal source [35] and it can be interpreted as the limit of point densities discretely distributed along the curve γ when this is partitioned without any limit [41,42]. On the other hand, only the total production rate of the well, Q(t), is known for general applications. This type of condition for the production rate is known as the uniform flow condition [21] and it is obviously related to the production rate per length unit and per unit of time q through the formula [35] Q(t) −
b
q γ (u), t γ (u) du.
(3)
a
3. Source term solution A very concise mathematical framework is exposed here in order to second the use of generalized functions (distributions) when modeling source terms in differential equations, and their corresponding source term solutions. It is assumed that most of these concepts are known by the reader. From the mathematical viewpoint, the integral in the right side of (2a), or more briefly denoted by the integral
q(γ , t)δ(x − γ ) dγ ,
γ
symbolizes the generalized function that operates under the rule [43]
(γ , t)δ(x − γ ) dγ , ϕ(x) =
γ
q(γ , t)ϕ(γ ) dγ
(4)
γ
for all test or “good” functions ϕ(x) which, in essence, are infinitely differentiable functions on Rn with compact support. Now, the equation ∂p(x, t) − p(x, t) = q(γ , t)δ(x − γ ) dγ (5) ∂t γ
can be interpreted as a differential equation in a distributional sense and only requires the generalized functions ∂p − p and (γ , t)δ(x − γ ) dγ ∂t γ
to be equal in the reservoir domain Ω for each t > 0 [25,35,43]. This distributional approach allows a natural interpretation of the physical situation and the use of the Fourier transform to obtain the solution of Eq. (5). The solution of Eq. (5) will be referred to as a source term solution, without consideration of boundary conditions. The expression fundamental solution will be used to denote the function traditionally known as free Green’s function. Using Fourier transform, Vladimirov [41,42] showed that a fundamental solution of the diffusivity operator is given by (3D case)
γ − x2 1 exp − . FS(x, t; γ , τ ) = (4π(t − τ ))3/2 4(t − τ )
(6)
R.J. González-Requena, J.M. Guevara-Jordan / Applied Numerical Mathematics 40 (2002) 433–449
437
On the other hand, the Fourier transform of the function f (x) is defined by the integral (provided it exists for every α) FTx [f ] ≡
···
f (x) exp −i(α ◦ x) dx1 · · · dxn = F (α).
(7)
Rn
The subscript x indicates that the Fourier transform has been applied to the variable x in Rn . The Fourier transforms of test functions are also test functions, which are adequately used to define the Fourier transforms of generalized functions [25,41–43]. In addition, if g(x, γ ) is a generalized function of x for all γ along the curve γ , it can be shown [25] that
g(x, γ ) dγ =
FTx γ
FTx g(x, γ ) dγ .
(8)
γ
To solve Eq. (5), we apply the Fourier transform to both sides and use expression (8). This generates an easy differential equation on t which solution is given by P (α, t) =
t
exp −α2 (t − τ ) − i(α ◦ γ ) q(γ , τ ) dγ dτ.
(9)
0 γ
This last expression is really the Fourier transform of the solution p of Eq. (5), so the solution p can be obtained using the inverse Fourier transform and is given by p(x, t) ≡ ST(x, t) =
t 0 γ
γ − x2 q(γ , τ ) exp − dγ dτ. (4π(t − τ ))3/2 4(t − τ )
(10)
4. Singularity programming The solution of (2) for a regular domain such as a parallelepiped or a rectangle is easily found through analytical methods [19]. When the frontier Ω is arbitrarily shaped, it is inevitable to resort to numerical methods, which, in general, are unable to adequately handle the Dirac delta-type function involved in the diffusivity [27]. According to the method here exposed, the solution p is considered as the addition of the source term solution ST, associated to the generalized function, and a smooth function, Preg , that can be solved through conventional numerical methods: p = Preg + ST.
(11)
It can be shown that Preg satisfies the following IBV problem ∂Rreg − Preg = 0 ∂t ∂ST ∂Preg (x, t) = − (x, t, γ ) ∂n ∂n Preg (x, 0) = 0,
in Ω × (0, ∞),
(12a)
in ∂Ω × (0, ∞),
(12b)
in Ω.
(12c)
The major advantage of (12) is that it represents a very conventional IBV problem for the diffusivity equation. It can be solved through any numerical method selected, or, given the case, by means of a
438
R.J. González-Requena, J.M. Guevara-Jordan / Applied Numerical Mathematics 40 (2002) 433–449
proper analytical method or a combination of both. In the present work, the Finite Element Method (FEM) was used to solve the reduced problem given by (12). Once Preg has been calculated, the pressure p solution of (2) is uniquely determined by (11). Lastly it is worth pointing out that using this method it possible to accept a trajectory for the horizontal well as curved and/or as twisted as can be mathematically handled.
5. Wellbore condition features The most important parameters used for the system of wells in reservoir simulation are flow rate and wellbore pressure. Consequently, two different conditions for the wellbore can be considered to address the problem under study (2). First, the uniform flux condition that assumes that the total production rate of the well is known, implying that the factor q(γ (u), t) in (10) is a known constant. Although the production rate is probably the best known information about oil wells, these kind of studies must be completed by including the capability of handling the wellbore pressure condition [3–6,24]. Here, the well production is conceived under the condition of uniform wellbore pressure, implying that the strength q(γ (u), t) of the production rate of the well is unknown. This pressure is defined as the reservoir pressure measured at the well radius and it is used in all cases to compute the unknown production rate of the well. As mentioned before, the proper solution of the reduced problem and the contribution of the corresponding source term solution give the required solution. However, the source term solution also deserves numerical treatment. For computational purposes, the spatial and time integrations involved in the source term solution (10) are discretely approximated. The following expression is then established:
ST x, t; γ (ui ), τj =
NU NT
q γ (ui ), τj C x, t; γ (ui ), τj γ (ui )ui τj .
(13)
j =1 i=1 U Here, {ui }N i=1 is a convenient partition of the interval [a, b]. Parameter u takes its values within this interval for generating each point of the well trajectory and is associated with the chosen numerical T method of integration. The set of times {τj }N j =1 is the corresponding subdivision of [0, t] linked to the selected method of numerical integration, with a last time step τNT equal (or very close) to t. Finally, C(x, t; γ (ui ), τj ) denotes the fundamental solution (6) evaluated at (γ (ui ), τj ). The discrete form (13) is used in Eq. (11) to obtain the numerical approximation of the pressure solution at any reservoir location x and at any instant t. Evidently, a corresponding numerical approximation is also used to treat the normal derivative of the source term that is present in the boundary condition of the reduced problem (12). When the uniform flux condition is imposed, the numerical solution given by (13) can be directly calculated due to the knowledge of q(γ (ui ), τj ). When the wellbore pressure condition is used, the strength q(γ (u), t) of the production rate of the well is unknown and must be previously established in order to find the solution using (13). Consequently, a different scheme must be proposed. NU Let {x w n }n=1 denote the wellbore positions (at the well radius) associated with the well trajectory, then, the necessary system of NU equations is given by
w pwf x w n , t = Preg x n , t +
NU NT
q γ (ui ), τj C x w n , t; γ (ui ), τ γ (ui ) ui τj
j =1 i=1
(14)
R.J. González-Requena, J.M. Guevara-Jordan / Applied Numerical Mathematics 40 (2002) 433–449
439
w for n = 1, 2, . . . , NU ; at each instant t, and where pwf (x w n , t) denotes the wellbore pressure at x n . Since the production rates q(γ (ui ), τj ) are known for j = 1, 2, . . . , Nτ − 1, the double summation associated with these indexes is called Γ (x w n , τNT −1 ) and the pressure differential is denoted by w w , t) = p (x , t) − p (x , t). Under these new variables, the system of equations (14) can be p(x w wf reg n n n written as NU
w w q γ (ui ), tNT C x w n , τNT ; γ (ui ), τNT −1 γ (ui ) ui τNT = p x n , t − Γ x n , τNT −1
(15)
i=1
for n = 1, 2, . . . , NU and at each instant t. These equations represent a generalized wellbore model for the horizontal well γ , because they relate the production rates with the pressure differential along the well. If wellbore pressures are given as input data, then this model is a linear system for the rates q(γ (ui ), τNT ). Its solution provides the production rates needed by the diffusivity equation solver in order to compute the pressure at the following time step. In this way, the extension of the new method to handle the wellbore pressure condition is achieved. In addition, the wellbore model can be used in reverse direction allowing the calculation of the wellbore pressure based on the rates. The wellbore model developed in this section is explicit. It means that the flow rate is determined before invoking the diffusivity solver on each time step. This kind of model is very simple but it may pose stability problems if long time steps are taken or strong changes in pressure and rate are observed. In these cases, shorter times steps could avoid this problem. If the times steps needed are too small then a more complex or implicit wellbore model is required. However, such development is not needed to extend and validate the numerical method presented in this article. More complex models for wells can be always developed following similar arguments to those presented in this section. Additional considerations will be presented to end this section. The discrete form of parameter u U induces a corresponding subdivision of points, {γ (ui )}N i=1 , on the well trajectory. These points can be considered as “gravity centers” of small arc elements composing the whole well trajectory. Moreover, the arc elements are here conceived as arc sources with respective wellbore pressures (concentrated at points γ (ui )) and production rates (given by q(γ (ui ), t)). This interpretation can be compared with the concept of Discrete Flux Elements (DFE) introduced by Azar-Nejad et al. [3–6] for modeling the potential (pressure) distribution generated by 2D finite sources: straight lines and curved lines. However, when these authors derive their definitions, they do it from results obtained for cases with straight horizontal wells parallel to coordinate axes, and in rectilinear reservoirs. Azar-Nejad et al. [3,5] assure that the source configuration must be straight to do the integral over the wellbore length. However, integrating over straight sources parallel to a coordinate axe implies that factor γ (u) in Eq. (10) is constant. Consequently, when the DFE concept is extrapolated to wells with arbitrary geometry (not twisted), the factor γ (u) is omitted, and it should be present in those integrals representing the curvilinear line sources. This length factor does not necessarily have to be constant for curved and twisted trajectories where, in fact, integration can also be done under adequate postulates of mathematical modeling. As has been seen, the factor γ (u) depends on the trajectory shape, and it has a different value at each point γ (u). Its absence could distort the accuracy of results, particularly when curved and twisted well trajectories are modeled.
440
R.J. González-Requena, J.M. Guevara-Jordan / Applied Numerical Mathematics 40 (2002) 433–449
6. Results and analysis 6.1. Finite element method The finite element method (FEM) [27,39] was used as the chosen numerical technique for solving the diffusivity equation in the numerical implementation of the SPM addressed by this paper. The PDE2D finite element software [39] was used to this end. It allows the solution of general systems of partial differential equations using polynomial bases of up to fourth order into 2D domains of arbitrary shape. FEM was selected because it is better to handle boundary conditions, since this is a major factor in implementing the SPM. Moreover, the FEM could be accommodated to tortuous (twisted) trajectories of horizontal wells using large and complex grid refinements [29], but the computational cost inherent to such procedure could render it unusable, particularly in 3D runs, unsteady flow cases, and domains with an arbitrary boundary. The polynomial base in all 2D examples was restricted to be of first order, which, under very general conditions [18], is equivalent to a scheme of finite differences. In 3D cases, the PDE2D software can only allow polynomial bases of third order into orthoedric-shaped boxes. This reflects the limited number of allowed blocks for studying 3D problems mostly due to restrictions in the available hardware (Sun Sparc II). The discretization caused in any dimension by PDE2D is completely implicit, preventing restrictions to the length of time steps based on the Courant–Friedrich–Lewy condition [27]. A more detailed FEM exposition lies beyond the purposes of this paper, as it focuses on showing the potential and flexibility of this approach for the tortuous well modeling, and to solving the selected cases. 6.2. Examples Four specific examples are presented below to validate the mathematical model, as well as the SPM used for the simulation of the pressure distribution generated by horizontal wells. A trapezoidal rule was used for the numerical integrals involved in these examples. 6.2.1. Example 1 (Validation against an analytic solution) For validation purposes, a rectangular reservoir given by Ω = [0, 2] × [0, 1] was selected. It contains a straight horizontal well γ of unit length parallel to the x axis and its center coincides with that of Ω, i.e., the coordinates (1, 0.5). An unitary uniform flux was assumed. Under these conditions, an analytical solution [19] and several numerical solutions obtained with two different methodologies were developed for solving the IBV problem (2). The numerical methodologies are the direct application of the FEM 2 to an equivalent problem, and the Singular Programming Methodology (SPM) proposed in this article, which also includes the use of the FEM. When using the direct application of the FEM, an approximation to the distributional representation of the horizontal well is considered, and the original differential equation (2a) becomes the following associated equation ∂pa − pa = ∂t
256, 0,
1 if (x, y) ∈ [0.5, 1.5] × 0.5 − 512 , 0.5 + otherwise.
1 512
,
2 The methodology that applies the FEM directly without removing singularities will be denoted by dFEM.
(16)
R.J. González-Requena, J.M. Guevara-Jordan / Applied Numerical Mathematics 40 (2002) 433–449
441
In this equation, the generalized function has been approximated by the indicator function shown at the right side of (16). It is well known [11] that these kinds of functions can be used at the limit process to approximate the singular function associated to the horizontal well representation. The validity of this function was verified by means of grid refinements and comparison with the analytical solution. Special grids were fitted to take into account the presence of the horizontal well, because the integration formulas can bypass the small band representing the well on uniform grids. The best approximation to the analytical solution using dFEM was found on a grid of, approximately, five thousand elements, at the dimensionless times 0.001, 0.01, 0.1 and 1. SPM was implemented using the same grids designed for dFEM. Although SPM does not need this kind of grid—it would work fine on global uniform grids—it was used just to match the dFEM working conditions and facilitate the comparison of the corresponding solutions. The best solution for the SPM was produced on a grid of, approximately, four hundred elements when it was compared with the analytical solution at the same times considered for dFEM. To obtain the convergence rate of the numerical methods in this test problem, the L2 (Ω) error of the difference between the numerical and analytical solutions was determined at the dimensionless time t −1. The pseudo-steady condition has been achieved at this moment, so no further changes in the solution are expected. Table 1 summarizes the errors for both cases and the grid dimension. Table 1 Errors in L2 (Ω) norm SPM Log(h)
Loganalytic − numerical
dFEM Grid
Log(h)
Loganalytic − numerical
Grid
−1.0986
−8.0102
72
−2.3979
−9.0612
968
−1.6094
−9.2567
200
−2.7081
−9.3093
1800
−1.9459
−9.6800
392
−3.2189
−9.7180
5000
Fig. 1. Example 1: L2 error vs. grid length. Logarithm values.
442
R.J. González-Requena, J.M. Guevara-Jordan / Applied Numerical Mathematics 40 (2002) 433–449
Table 2 CPU times for SMP and dFEM Number of points where solution is required
Normalized CPU times: time_SPM/time_dFEM
10 points
0.60
50 points
0.62
150 points
0.74
315 points
0.98
The contents of the table are shown in Fig. 1, where the slopes of the straight lines allow to conjecture the order of convergence for each method. The magnitude h represents the maximum length of all elements in the grid. Notice that the dFEM and the SPM, respectively, present almost first and second order of convergence. For each method, the CPU times were computed on a Sun workstation (Sparc II) on the respective finer grid. CPU times were found to be almost equal for both methods, in spite of the fact that dFEM was applied in a grid twelve times bigger than the one used in the SPM. However, if the solution is needed in a reduced number of points then SPM may be twice as fast as FEM as observed in Table 2. In any case, CPU times should be carefully analyzed because this is really a comparison between an efficient commercial code [39], and a domestic mixed algorithm developed for research purposes. 6.2.2. Example 2 (Applying the SPM to 2D irregular reservoirs) The SPM will be used in this example to numerically simulate the pressure behavior created by a horizontal well in a reservoir of irregular boundary. To this end, region Ω in Fig. 2(a) was considered with the same initial and boundary conditions given in (2). As in the previous example, well length is unity, and its total rate production is of magnitude one. Following the same procedure as before, we will also solve the associated equation to (2) where the generalized function is approximated by an indicator function. The solution of the associated problem will be obtained by directly using the FEM. This associated equation is given by
1 1 ∂pa 256 if (x, y) ∈ [0.5, 1.5] × 0.602 − , 0.5 + , 512 512 − pa = (17) ∂t 0 otherwise. Here, the well center is at the point (1, 0.602). Figs. 2(a) and (b) show the grids used for the application of SPM and dFEM, respectively. Both solutions converge to the correct solution. In this case, a numerical method has been assumed to converge on a certain grid if the difference between numerical solutions obtained in successive refinements are not punctually bigger than 10e−6. Again, it is worth mentioning that this kind of grids, which takes into account the presence/location of the horizontal well, are required only when using the FEM to solve the associated problem (17), but they are unnecessary with the application of the SPM for solving (2). This feature will be highlighted in the 3D example. The following figures represent a sequence of the pressure distribution obtained by application of SPM. Fig. 3(a) shows the drainage area generated by the horizontal well at dimensionless time t = 0.01. It can be observed in this case that pressure close to the well behaves as if the reservoir was unlimited. Additionally, the elliptical pressure isolines unambiguously define the geometric location of the well in the reservoir. Fig. 3(b) shows the final drainage area, which reflects pressure distribution at dimensionless
R.J. González-Requena, J.M. Guevara-Jordan / Applied Numerical Mathematics 40 (2002) 433–449
443
(a)
(b) Fig. 2. Example 2: (a) Coarse grid system for numerical solutions by SPM. (b) Refined grid system for numerical solutions by dFEM.
time t = 1.0 when the pseudo-stationary condition is reached. At this time, the drainage area is the whole reservoir and it has thus ceased to expand. Also, it can be observed that the contour lines close to the well are not concentric ellipses around the horizontal well, but rather around its center, which is a consequence of the adopted modeling condition, i.e., uniform production rate. For this constant flux case, a pressure gradient along the well is present because the flow has to enter the well uniformly along its length [9]. Such pressure behavior is also typical for fractured vertical wells, where fluid extracted through the fracture tends to move towards the center, onto the position where the vertical well is supposed to be pumping fluid to the surface. This is the reason why many of the available results for fractured wells are frequently applicable to horizontal wells. It is worth noting that fluid is not always pumped out of horizontal wells through their centers, but rather through one of its ends or using producing sub-intervals. Clearly, these situations require different sets of hypotheses for modeling purposes, but they can be also addressed by the SPM. In conclusion, this example shows the SPM capability to simulate the unsteady behavior of the pressure generated by a horizontal well in an arbitrarily shaped 2D reservoir. In addition, using the SPM in less
444
R.J. González-Requena, J.M. Guevara-Jordan / Applied Numerical Mathematics 40 (2002) 433–449
(a)
(b) Fig. 3. Example 2: (a) Pressure isolines around the horizontal well at time t = 0.01. (b) Isopressure lines around the horizontal well at time t = 1.0.
refined grids solutions were obtained which were comparable to those obtained with the dFEM for the equivalent problem in a more refined grid. 6.2.3. Example 3 (Applying the SPM to a well producing by wellbore pressure in 2D irregular reservoirs) This example gives us the opportunity of applying the wellbore model developed here to include wellbore pressure condition as input data. The reservoir and horizontal well are the same used in Example 2. Indeed, all the conditions of that example will be the same up to the simulation time t = 1. At this point the wellbore model is used to keep the wellbore pressure along the horizontal well at −0.51 from instant t = 1 to the time t = 1.5. This restriction on the boundary condition produces a completely different pressure behavior around the well. The pressure field at time t = 1.5 is presented in Fig. 4(a). It shows that the isopressure lines do not intercept the horizontal well because all the wellbore pressures, at this instant, have the value −0.51 imposed as a restriction.
R.J. González-Requena, J.M. Guevara-Jordan / Applied Numerical Mathematics 40 (2002) 433–449
445
(a)
(b) Fig. 4. Example 3: (a) Pressure distribution around the horizontal well at time t = 1.5. (b) Effects of wellbore model on pressure and rate.
In addition, Fig. 4(b) shows the effects of the wellbore model on the rate and pressure evaluated at the center of the horizontal well. It can be observed that the flow rate is constant until the time t = 1; after this moment the well rate is reduced in time by the wellbore model in order to keep the wellbore pressure constant at −0.51. The small oscillations of pressure around this value are a consequence of the explicit coupling of the wellbore model. All these observations of the behavior of pressure and production for the horizontal well are consistent with the physics of the problem at hand, so it may be concluded that the wellbore model developed for the SPM gives a satisfactory extension of this method. This extension depends neither on the reservoir boundary shape nor on a tortuously curved/twisted well trajectory. 6.2.4. Example 4 (Applying the SPM to wells with curved and twisted trajectory) The objective of this relevant case is to establish the unsteady pressure distribution generated by a horizontal well whose trajectory is as curved and twisted as a circular helix. The problem to be solved is
446
R.J. González-Requena, J.M. Guevara-Jordan / Applied Numerical Mathematics 40 (2002) 433–449
the same as expressed in (2). Here, γ (u) is the analytical representation of this tortuous trajectory [20], which is obtained by rotating the helix of equation (0.8u, 0.5 + 0.25 cos(π u), 0.5 + 0.25 sin(π u)), until its own axis is aligned with the direction of one of the major diagonals of the domain Ω = [0, 2] × [0, 1] × [0, 1]. Parameter u takes values between a = 0 and b = 2.5. Figs. 5(a) and (b) show a transient sequence for the pressure response generated by this twisted well at times 0.001 and 0.1, respectively. Under the flow scheme known as early radial, Fig. 5(a) shows pressure distribution (isosurfaces) honoring the curved and twisted shape of the helix-shaped trajectory of this well. This is similar to the elliptic pressure isolines that make up the elliptic drainage area reserved for
(a)
(b) Fig. 5. Example 4: (a) Pressure isosurfaces around a tortuous well at dimensionless time t = 0.001. (b) Pressure isosurfaces around a tortuous well at dimensionless time t = 0.1.
R.J. González-Requena, J.M. Guevara-Jordan / Applied Numerical Mathematics 40 (2002) 433–449
447
these times when 2D studies of straight-well response are performed. This can be also seen in Example 2 and in [3–6]. Fig. 5(b) presents the shape of the pressure isosurfaces resulting from the spherical flow at both ends of the well. Equivalent results have been reported by Azar-Nejad et al. [6] in their research with straight horizontal wells in a 2D rectangular reservoir. This application illustrates both the potential for emulating well trajectories as curved and twisted as the relevant mathematical representation can handle, as well as the total independence between the geometry/location of the well and the proposed grid for obtaining the numerical solution of the reduced problem.
7. Conclusions In this paper, the SPM was introduced and validated like a new methodology for assessing the transient pressure distribution generated by tortuous horizontal wells in arbitrary shaped reservoirs. In the 2D examples considered, the SPM presented excellent results even on coarse grids, and it was numerically established that its order of convergence is quadratic. Additionally, the SPM proved to have faster numerical convergence than the convergence reached by the direct application of the FEM to equivalent problems. Where possible, solutions using the SPM in coarse grids compared favorably with those obtained applying the dFEM in more refined grid. This evidenced the superiority of this methodology. The SPM applies equally well to wells that are producing either by constant rate of production or by uniform wellbore pressure. The SPM was shown to have potential for emulating well trajectories as curved and twisted as the relevant mathematical representation can treat. Likewise, this methodology has shown that arbitrarily shaped reservoirs are not a limitation. Summarizing, the main SPM advantages are the convergence in coarse grids; its elegant capability for emulating horizontal well trajectories as curved and twisted as required, and the total independence between numerical grids and the geometry/location of the trajectory well.
References [1] W.F. Ames, Numerical Methods for Partial Differential Equations, Academic Press, New York, 1977. [2] C. Astudillo, E. Guedez, Tratamiento de singularidades en el modelaje de pozos de un yacimiento, petrolífero, Graduate Thesis, Mathematics Department, School of Sciences, Universidad Central de Venezuela, Caracas, Venezuela, 1993. [3] F. Azar-Nejad, W. Tortike, S. Farouq Ali, Potential distribution around the sources with finite length (horizontal, vertical partially penetrating wells and fractures). Part II: Transient flow, Paper 35269, Society of Petroleum Engineers, 1996. [4] F. Azar-Nejad, W. Tortike, M. Farouq Ali, Potential distribution around the sources with finite length (horizontal, vertical partially penetrating wells and fractures). Part I: Steady state fluid flow, Paper 35270, Society of Petroleum Engineers, 1996. [5] F. Azar-Nejad, W. Tortike, S. Farouq Ali, 3D analytical solution to the diffusivity equation for finite sources with applications to horizontal wells, Paper 35512, Society of Petroleum Engineers, 1996.
448
R.J. González-Requena, J.M. Guevara-Jordan / Applied Numerical Mathematics 40 (2002) 433–449
[6] F. Azar-Nejad, W. Tortike, S. Farouq Ali, Performance of horizontal wells with irregular geometry, Paper 36550, Society of Petroleum Engineers, 1996. [7] D.K. Babu, A. Odeh, Productivity of horizontal wells, Paper 18298, Society of Petroleum Engineers, 1985. [8] D.K. Babu, A. Odeh, Transient flow horizontal well, pressure drawdown, and buildup analysis, Paper 18802, Society of Petroleum Engineers, 1989. [9] R.M. Butler, Horizontal Wells for the Recovery of Oil, Gas and Bitumen, Monograph, Vol. 2, The Petroleum Society of the Canadian Institute of Mining, Metallurgy and Petroleum, Calgary, AB, 1994. [10] H.S. Carslaw, J. Jaeger, Conduction of Heat in Solids, Oxford University Press, Oxford, 1959. [11] G. Chavent, J. Jaffre, Mathematical models and finite elements for reservoir simulation, in: J.L. Lions, G. Papanicolaou, H.F. Fujita, H.B. Keller (Eds.), Studies in Mathematics and its Applications, Vol. 17, Elsevier, Amsterdam, 1986. [12] B.C. Craft, M. Hawkins, Applied Petroleum Reservoir Engineering, Prentice-Hall, Englewood Cliffs, NJ, 1957. [13] L.P. Dake, Fundamentals of Reservoir Engineering, Elsevier, Amsterdam, 1978. [14] M.P. Do Carmo, Differential Geometry of Curves and Surfaces, Prentice-Hall, Englewood Cliffs, NJ, 1976. [15] J. Douglas, R. Ewing, M. Wheeler, Approximation of the pressure by a mixed method in the simulation of miscible displacement, RAIRO Analyse Numerique 17, 1983. [16] P. Duchateau, D.W. Zachmann, Partial Differential Equations, McGraw-Hill, New York, 1988. [17] R. Ewing, Problems arising in the modeling of processes for hydrocarbon recovery, in: R.E. Ewing (Ed.), The Mathematics of Reservoir Simulation, SIAM, Philadelphia, PA, 1983, Chapter 1. [18] R. Ewing, M.F. Wheeler, Galerkin Methods for Miscible Displacement Problems with Point Sources and Sinks-Unit Mobility Ratio, Mathematical Methods for Energy Research, SIAM, Philadelphia, PA, 1983. [19] A. Fermin, Método de Green para la simulación numérica de pozos horizontales, Graduate Thesis, Mathematics Department, School of Sciences, Universidad Central de Venezuela, Caracas, Venezuela, 1998. [20] R. González, Metodologías alternas para la solución numérica de la ecuación de difusión bajo condiciones de borde pertinentes a las aplicaciones en simulación de yacimientos, Ph.D. Thesis, Mathematics Department, School of Sciences, Universidad Central de Venezuela, 2000. [21] A.C. Gringarten, H.J. Ramey, The use of source and Green functions in solving unsteady flow problems in reservoirs, Paper 3818, Society of and Petroleum Engineers, 1973. [22] J. Guevara, A. Fermín, R. González, A new approach for horizontal wells singularities in petroleum engineering, Paper 51924, Society of Petroleum Engineers, 1999. [23] L. Hayes, R. Kendall, M. Wheeler, Treatment of source and sinks in steady-state reservoir engineering simulations, in: R. Vichnevetsky (Ed.), Advances in Computer Methods for Partial Differential Equations II, IMACS (AICA), 1977. [24] J.K. Jasti, V.R. Penmatcha, D.K. Babu, Use of analytical solutions in improvement of simulator accuracy, Paper 38887, Society of Petroleum Engineers, 1997. [25] D.S. Jones, Generalized Functions, McGraw-Hill, New York, 1966. [26] L.S. Koh, D. Tiab, A boundary element algorithm for modeling 3D horizontal well problems using 2D grids, Paper 26228, Society of Petroleum Engineers, 1993. [27] L. Lapidus, G. Finder, Numerical Solution of Partial Differential Equations in Sciences and Engineering, Wiley, New York, 1988. [28] J. Liggett, P. Liu, Boundary Integral Equation Method for Porous Media Flow, George Allen & Unwin, London, 1983. [29] N. Morita, Transient finite element code: a versatile tool for well performance analysis, Paper 25878, Society of Petroleum Engineers, 1993. [30] B. O’Neill, Elementary Differential Geometry, Academic Press, New York, 1976. [31] E. Ozkan, R. Raghavan, New solution for well test analysis problems: Part 1—analytical considerations, Paper 18615, Society of Petroleum Engineers, 1991.
R.J. González-Requena, J.M. Guevara-Jordan / Applied Numerical Mathematics 40 (2002) 433–449
449
[32] E. Ozkan, R. Raghavan, New solution for well test analysis problems: Part 2—computational considerations and applications, Paper 18616, Society of Petroleum Engineers, 1991. [33] E. Ozkan, R. Raghavan, New solution for well test analysis problems: Part 3—additional algorithms, Paper 23593, Society of Petroleum Engineers, 1994. [34] R. Raghavan, Well Test Analysis, Penn Well, Tulsa, OK, 1993. [35] R. Raghavan, E. Ozkan, A Method for Computing Unsteady Flows in Porous Media, Longman Scientific and Technical, Harlow, UK, 1994. [36] W.A. Roper, J. Merchant, C. Duvall, Combination of numerical and analytical techniques to improve waterflood model efficiency, Paper 2031, Society of Petroleum Engineers, 1962. [37] A.J. Rosa, R. Carvalho, A mathematical model for pressure evaluation in an infinite conductivity horizontal well, Paper 15967, Society of Petroleum Engineers, 1989. [38] T.F. Russell, M.F. Wheeler, Finite elements and finite differences methods for continuous flow in porous media, in: R.E. Ewing (Ed.), The Mathematics of Reservoir Simulation, SIAM, Philadelphia, PA, 1983, Chapter 2. [39] G.V. Sewell, Analysis of Finite Element Method: PDE2D/Protran, Springer, Berlin, 1985. [40] V. Villamizar, Personal communication, Mathematics Department, School of Sciences, Universidad Central de Venezuela, 1998. [41] V.S. Vladimirov, Equations of Mathematical Physics, Marcel Dekker, New York, NY, 1971. [42] V.S. Vladimirov, Generalized Functions in Mathematical Physics, Mir Publishers, Moscow, 1979. [43] A.H. Zemanian, Distribution Theory and Transform Analysis, Dover, New York, 1987.