On the Analysis and Numerical Formulation of Miscible Fluid Flow in Porous Media Using Chebyshev Wavelets Collocation Method ()
Received 16 May 2016; accepted 8 July 2016; published 11 July 2016
1. Introduction
Understanding the dynamics of fluid flow in porous media is undoubtedly a subject interesting to many scientists and engineers due to its applications in many areas of research [1] . In many occasions, research efforts have been made both experimentally and theoretically to better explain the dynamics of single-phase flow as well as heat transfer through varying porous media [2] - [6] , mostly encountered in diverse fields of science and engineering. In order to determine the appropriate conditions of operations as well as the machinery to use for the various operations and jobs, a number of parameters and the governing equations have to be predicted accurately. Several of the cases assumed that the flow follows a Darcy’s law [6] - [8] which was the fundamental principle used to describe the flow of fluids in a reservoir; however, some researchers have extended their work to non- Darcy flow [2] - [4] [9] . The flow in porous media has been studied by several authors using methods like the Implicit Pressure Explicit Saturation (IMPES) method [10] - [13] , fully implicit method [14] [15] , the finite volume method [16] , cell-centered finite difference method [17] , discontinuous Galerkin Method [18] , and sequential methods [15] [19] [20] .
Wavelet theory is a recent methodology which has a wide range of applications in many research areas in mathematics, physics, and engineering. The applicability of wavelets in different areas of research is their ability to represent a wider class of functions and operators efficiently and accurately. Unlike most of the methods used in solving fluid flow problems, a carefully constructed wavelet based method will be capable of capturing aspects of the function which other functional analysis methods may miss such as trends, breakdown points, discontinuities and self similarities [21] as well as integrating different data types [22] .
The use of wavelets methods for solving partial differential equations goes as far back as early 1990s [23] . Different wavelet families are applied in various papers for solving differential equations, in which the wavelet coefficients are computed based on either the Galerkin or Collocation methods [18] [24] . Application of wavelets methods reduces the problems to systems of algebraic equations.
In this study, the Chebyshev wavelets method constructed using the first kind of Chebyshev polynomial is proposed for numerically simulating the single-phase flow of fluids in a reservoir. The solution to the flow equations is obtained in one-dimension. The governing models of single-phase flow in reservoir are presented in Section 2. The Chebyshev wavelets are discussed in Section 3. In Section 5, the Chebyshev wavelets formulation of the single-phase flow model is presented. In Section 6, results from the numerical simulation are presented, discussed and finally we conclude in Section 7.
2. Single-Phase Flow Model
The single-phase flow in a porous medium considers the only one fluid phase or several completely miscible fluids in the reservoir at any given time. The single-phase flow in a porous medium is governed by partial differential equation resulting from the principle of mass conservation. This equation is given as [6]
(1)
where 
is density of fluid; 
is porosity; and u represents the velocity of the fluid which is related to the pressure gradient through the Darcy’s law. The Darcy’s law gives the velocity of the fluid to be
(2)
where K is the intrinsic permeability; 
the fluid viscosity and P represents the pressure at position x at time t. The porosity of the reservoir and the density of the fluid phase are both functions of pressure, and for an isothermal system
(3)
Introducing the compressibility relationship for the reservoir and fluid given as
(4)
then
(5)
Consequently, assuming that the permeability of the reservoir and the viscosity of the fluid are constant we have
(6)
We arrive at Equation (7) by substituting Equation (5) and Equation (6) into Equation (1).
(7)
Since the compressibility of the fluid is mostly small and for a low velocity flow the pressure gradient is also small, which is to say that 
and
, we obtain the single-phase flow model as
(8)
Initial condition is set for the flow problem at 
and the boundary conditions are set at 
and ![]()
for![]()
.
![]()
(9)
![]()
(10)
![]()
(11)
3. Principles of Wavelet Transform
Wavelets are basis functions obtained from dilating and translating a single function known as the mother wavelet. The wavelet transform is an integral operator obtained by taking the inner product of a function wavelets. The transform of a given function ![]()
is [25]
![]()
(12)
where a is the scaling parameter and b is the translation parameter. Given that the dilation parameter a and translation parameter b vary continuously then the wavelet family is said to be continuously [26] [27]
![]()
(13)
A wavelet function ![]()
satisfies the wavelet admissibility condition given as
![]()
(14)
where ![]()
is the Fourier transform of the wavelet function. To ensure![]()
, then
![]()
(15)
Restricting the a and b to assume discrete values as ![]()
and![]()
, where![]()
, gives a family of discrete wavelets [28] :
![]()
(16)
The family of wavelets form a wavelet basis of the ![]()
and these functions form an orthonormal basis if ![]()
and ![]()
[25] .
3.1. Chebyshev Polynomial
The Chebyshev polynomials are the eigenfunctions of singular Sturm-Liouville problem with many advantages. The Chebyshev polynomial ![]()
in this research is that of the first kind, with m been the degree of the polynomial. This polynomial can generally be represented by the recurrence relation
![]()
(17)
where M is a fixed positive integer greater than 2.
A set of Chebyshev Polynomials, ![]()
, are orthogonal with respect to the weight function ![]()
on the interval ![]()
[29] .
3.2. Chebyshev Wavelets
The Chebyshev wavelets are constructed based on the set of Chebyshev Polynomials. The Chebyshev wavelets ![]()
is a function of four arguments where![]()
, and m represents the order of the Chebyshev polynomial, is defined on the interval ![]()
as
![]()
(18)
where
![]()
(19)
To ensure orthogonality in dealing with the Chebyshev wavelets, the weight function ![]()
has to be dilated and translated as [30]
![]()
Given that ![]()
and![]()
; ![]()
and![]()
, the expansions are obtained for the Chebyshev wavelets as
![]()
![]()
![]()
![]()
![]()
![]()
3.3. Approximating Function
Any square intgrable function ![]()
may be representated by the Chebyshev wavelets as
![]()
(20)
where ![]()
in which ![]()
represents the inner product. The infinite series in Equation (20) can be truncated for finite values of n and m. The wavelet decomposition of the function ![]()
can then be written as
![]()
(21)
where C and ![]()
are ![]()
matrices given by
![]()
![]()
(22)
Likewise, a two variable function ![]()
defined on the square ![]()
which is square integrable can as well be expanded using the Chebyshev wavelets basis as:
![]()
(23)
where ![]()
is a ![]()
matrix.
3.4. Chebyshev Operational Matrix of Integration
The operational matrices are used to integrate or different a function (set of functions). This matrix was introduced by Chen and Hsiao in 1975 [26] . Given that Q is defined as ![]()
operational matrix for integration [31] , the integral of the vector ![]()
defined in Equation (22) can be obtained as
![]()
(24)
In this research, the operational matrix of integration Q is derived based on the Chebyshev wavelets. This is demonstrated for ![]()
and![]()
, the six Chebyshev basis functions can be integrated and the represented in terms of the wavelets function using the definition of the inner product which gives rise to the following results,
![]()
![]()
![]()
![]()
![]()
![]()
Based on Equation (24), the operational matrix of integration is obtained as
![]()
(25)
4. Convergence Analysis
The convergence of the Chebyshev wavelets basis is indicated in this section.
Theorem 1. If a function![]()
, with bounded second order derivative ![]()
can be expan- ded as a sum of infinite Chebyshev wavelets
![]()
(26)
then
![]()
(27)
which means the Chebyshev wavelets expansion converges uniformly to![]()
.
Proof. For the proof of this theorem you are referred to [32] .
Theorem 2. Let ![]()
be a continuous function defined on the interval ![]()
with second derivatives ![]()
bounded by![]()
, then the accuracy estimation
![]()
(28)
where
![]()
(29)
Proof. For the proof of this theorem you are referred to [33] .
5. Model Decomposition
The flow dynamics of fluids in the porous medium requires adequately solving the equations governing the flow process. In this section, the Chebyshev wavelets collocation method is used to analyze the governing equation of the single-phase flow in a porous medium. The unknown function, ![]()
in the problem is defined and decomposed at this stage using the Chebyshev wavelets in Equation (30).
![]()
(30)
Integrate with respect to t, making use of the Chebyshev operational matrix of integration (Q),
![]()
![]()
(31)
Substituting the initial condition![]()
, we obtain
![]()
(32)
Integrating Equation (30) twice with respect to x together with the boundary conditions gave
![]()
![]()
(33)
Substituting the boundary conditions, for ![]()
![]()
(34)
Therefore
![]()
(35)
Equation (32) and Equation (35) are substituted into Equation (8) which results to the wavelet representation of the single phase model governing the fluid flow in the media given in Equation (36).
![]()
(36)
Taking collocation points
![]()
(37)
we obtain a set of linear algebraic equations from Equation (36) based on the collocation points. These set of linear equations are solved for D. In order to reconstruct the pressure, ![]()
from its wavelet coefficient we obtain the expression for ![]()
by integrating Equation (35) with respect to t which gives
![]()
(38)
6. Simulation Results
In this section, we present the numerical solution of the single-phase flow model given in Equation (1) using the solution (Equation (38)) from Chebyshev wavelet collocation method discussed above and compare with the exact solution given in Equation (39) below:
![]()
(39)
This will allows us to estimate the error associated with the Chebyshev wavelet collocation method. The various parameter values necessary for describing the flow dynamics can be measured or estimated for the medium of choice and the occupying fluids. In this study, parameters for the numerical simulation of the flow problem were taken from Unsal et al. [34] . Chosen parameter values are for oil as the fluid phase flowing through a porous medium. The efficiency of the Chebyshev wavelets method was measured using [28] [34] the root mean square error (RMSE) and the maximum error (![]()
)
![]()
(40)
![]()
(41)
The wavelet equations resulting from the decomposition of the single phase flow model was solved algebraically for the wavelet coefficients and used to reconstruct the solution to the single phase flow model. The approximate solution of the pressure distribution in the reservoir obtained from the use of the Chebyshev wavelet collocation method and the corresponding exact solution are present in Figure 1.
The simulation provided estimates of the evolution of pressure of oil through the reservoir on which the pressure can either be maintained, increased or decrease to achieve the required quantity of oil to be produced. This will inform management on the necessary action to take to achieve target production. The absolute errors in the approximation of the pressure evolution are shown in Figure 2. Absolute error between analytic and numerical solutions in different time period of the flow is shown in Figure 3.
The simulation results based on the Chebyshev wavelet method, compared to the exact solution have fairly small errors measured which makes the Chebyshev wavelet method very efficient and accurate in approximating the pressure distribution in the reservoir from the flow model. Figure 4 shows a plot of the numerical solution compared to the exact at selected times in the reservoir showing a close approximation of the exact solution. In Figure 5, the plot shows the numerical solution at different depths of the reservoir.
The pressure in the medium was observed to drop gradually over the time period as seen in Figure 5. The pressure is greater for higher values of the spatial variable (depth). Figure 6 is a time plot of the root mean square errors and the maximum absolute error demonstrating the efficiency of the method. Some of the error values are tabulated in Table 1.
![]()
Figure 1. Pressure distribution in the reservoir.
![]()
Figure 2. Absolute error in different flow times t.
![]()
Figure 3. Absolute error in Chebyshev wavelets approximation of pressure in the medium for single-phase flow.
![]()
Figure 4. Comparing the numerical solution to the exact solution of pressure in the porous medium.
Root Mean Square Error Estimate and Maximum Absolute Error values for pressure distribution through the reservoir were calculated at each time within the simulation period. These were calculated using Equation (41). Table 1 shows respective values for the two error estimators at sample time intervals.
![]()
Figure 5. Pressure approximation by Chebyshev wavelets method over different time periods.
![]()
Figure 6. Time series of root mean square error (RMSE) and maximum absolute error (![]()
).
![]()
Table 1. Root mean square error estimate (RSME) and maximum absolute error (MAE) at time t.
7. Conclusion
This paper presented the Chebyshev wavelet method together with the operational matrix of integration as a numerical scheme for approximating the pressure distribution of the single phase flow of fluid in a porous medium. Both Root Mean Square Error Estimate ranging from 0.02 to 0.11 and the maximum absolute error between 0.04 and 0.21 indicate that the method is efficient for simulating the flow process. Aside the capacity to capture trends, breakdown points, discontinuities, self similarities in functions, the use of the Chebyshev wavelet method incorporates the boundary condition of the problem automatically making the method very convenient for solving boundary value problems. The problem is subsequently reduced to a set of algebraic equations. The resulting system of equations was solved to obtain the wavelet coefficient of the unknown functions from which the solution to the flow problem was reconstructed.