On the Analysis and Numerical Formulation of Miscible Fluid Flow in Porous Media Using Chebyshev Wavelets Collocation Method


In this paper, the Chebyshev wavelet method, constructed from the Chebyshev polynomial of the first kind is proposed to numerically simulate the single-phase flow of fluid in a reservoir. The method was used together with the operational matrices of integration which resulted in an algebraic system of equations. The system of equation was solved for the wavelet coefficient and used to construct the solutions. The efficiency and accuracy of the method were demonstrated through error measurements. Both the root mean square and the maximum absolute error analysis used in the study were within significantly close range. The Chebyshev wavelet collocation method subsequently was observed to closely approximate the analytic solution to the single phase flow model quite well.

Share and Cite:

Amoako-Yirenkyi, P. , Awashie, G. and Dontwi, I. (2016) On the Analysis and Numerical Formulation of Miscible Fluid Flow in Porous Media Using Chebyshev Wavelets Collocation Method. Journal of Applied Mathematics and Physics, 4, 1210-1221. doi: 10.4236/jamp.2016.47126.

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]


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


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


Introducing the compressibility relationship for the reservoir and fluid given as




Consequently, assuming that the permeability of the reservoir and the viscosity of the fluid are constant we have


We arrive at Equation (7) by substituting Equation (5) and Equation (6) into Equation (1).


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


Initial condition is set for the flow problem at and the boundary conditions are set at and for.




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]


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]


A wavelet function satisfies the wavelet admissibility condition given as


where is the Fourier transform of the wavelet function. To ensure, then


Restricting the a and b to assume discrete values as and, where, gives a family of discrete wavelets [28] :


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


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




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


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


where C and are matrices given by


Likewise, a two variable function defined on the square which is square integrable can as well be expanded using the Chebyshev wavelets basis as:


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


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


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




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




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).


Integrate with respect to t, making use of the Chebyshev operational matrix of integration (Q),


Substituting the initial condition, we obtain


Integrating Equation (30) twice with respect to x together with the boundary conditions gave


Substituting the boundary conditions, for




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).


Taking collocation points


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


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:


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 ()



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.

Conflicts of Interest

The authors declare no conflicts of interest.


[1] Scheidegger, A.E. (1974) The Physics of Flow through Porous Media. University of Toronto Press, Toronto.
[2] Chai, Z., Shi, B., Lu, J. and Guo, Z. (2010) Non-Darcy Flow in Disordered Porous Media: A Lattice Boltzmann Study. Computers and Fluids, 39, 2069-2077.
[3] Bernsdorf, J., Hasert, M. and Roller, S. (2011) Lattice Boltzmann Simulation of Non-Darcy Flow in Porous Media. International Conference on Computational Science, Procedia Computer Science, 4, 1048-1057.
[4] Holditch, S.A. and Morse, R.A. (1976) The Effects of Non-Darcy Flow on the Behavior of Hyfraulically Fractured Gas Wells.
[5] Zeng, Z. and Grigg, R. (2006) A Criterion for Non-Darcy Flow in Porous Media. Transport in Porous Medium, 63, 57-69.
[6] Soulaine, C., Davit, Y. and Quintard, M. (2013) A Two-Pressure Model for Slightly Compressible Single Phase Flow in Bi-Structed Porous Media. Oxford Centre for Collaborative Applied Mathematics, 96, 55-70.
[7] Hassanizadeh, S.M. and Gray, W.G. (1993) Toward an Improved Description of the Physics of Two-Phase Flow. Advances in Water Resources, 16, 53-67.
[8] Unsal, E., Matthai, S.K. and Blunt, M.J. (2009) Simulation of Multiphase Flow in Fractured Reservoirs Using a Fracture-Only Model with Transfer Functions. Computer Geosciences, 14, 527-538.
[9] Al-Otaibi, A. and Wu, Y.S. (2010) Transient Behavior and Analysis of Non-Darcy Flow in Porous and Fractured Reservoirs according to the Barree and Conway Model. SPE International.
[10] Chen, Z., Li, B. and Huan, G. (2004) An Improved Impes Method for Two-Phase Flow in Porous Media. Transport in Porous Medium, 3, 361-376.
[11] Foroozesh, J., Barzegari, D., Ayatollahi, Sh. and Jahanmiri, A. (2008) Simulation of Water Comning in Oil Reservoir Using a Corrected Impes Method. Iranian Journal of Chemical Engineering, 5, 4.
[12] Sheldon, J.W. and Cardwell, W.T. (1959) One-Dimensional, Incompressible, Noncapillary, Two-Phase Fluid Flow in a Porous Medium. Society of Petroleum Engineers, 216, 290-296.
[13] Young, L.C. and Russell, T.F. (1993) Implementation of an Adaptive Implicit Method. SPE Symposium on Reservoir Simulation, New Orleans, 28 February-3 March 1993, 14 p.
[14] Douglas, Jr., J., Peaceman, D.W. and Rachford Jr., H.H. (1959) A Method for Calculating Multi-Dimensional Immiscible Displacement. Society of Petroleum Engineers, 216, 297-308.
[15] Li, B., Chen, A. and Huan, G. (2003) The Sequential Method for the Black-Oil Reservoir Simulation on Unstructured Grids. Journal of Computational Physics, 192, 36-72.
[16] Eymard, R. and Gallouet, T. (2002) Finite Volume Schemes for Two-Phase Flow in Porous Media. Proceedings of ALGORITMY, Vysoke Tatry, 8-13 September 2002, 31-40.
[17] El-Amin, M.F., Meftah, R., Salama, A. and Sun, S. (2015) Numerical Treatment of Two-Phase Flow in Porous Media Including Specific Interfacial Area. International Conference on Computational Science, 51, 1249-1258.
[18] Epshteyn, Y. and Riviere, B. (2005) On the Solution of Imcompressible Two-Phase Flow by a P-Version Disconti-Nous Galerkin Method. Communications in Numerical Methods in Engineering, 22, 741-751.
[19] Chen, Z. (2001) Formulations and Numerical Methods of the Black Oil Model in Porous Media. Society for Industrial and Applied Mathematics Journal on Numerical Analysis, 38, 489-514.
[20] Chen, Z., Qin, G. and Ewing, R.E. (2006) Analysis of a Compositional Model for Fluid Flow in Porous Media. Society for Industrial and Applied Mathematics, 60, 747-777.
[21] Yue, W. and Guo, T. (2003) Reservoir Fluid Identification by the Wavelet Transform. Chinese Journal of Geophysics, 46, 1241-1250.
[22] Lu, P. and Horne, R.N. (2000) A Multiresolution Approach to Reservoir Parameter Estimation Using Wavelet Analysis. Society of Petroleum Engineers.
[23] Lepik, U. (2011) Solving Pdes with the aid of Two-Dimensional Haar Wavelets. Computers and Mathematics with Applications, 61, 1873-1879.
[24] Kou, J. and Sun, S. (2010) A New Treatment of Capillarity to Improve the Stability of Impes Two-Phase Flow Formulation. Elsevier Computers and Fluids, 39, 1923-1931.
[25] Calderbank, A.R., Daubechies, I., Sweldens, W. and Yeo, B. (1998) Wavelet Transforms That Map Integers to Integers. Applied and Computational Harmonic Analysis, 5, 332-369.
[26] Chih-Fan, C. and Chi-Huang, H. (1975) Design of Piecewise Constant Gains for Optimal Control via Walsh Functions. IEEE Transactions on Automatic Control, 20, 596-603.
[27] Guf, J.S. and Jiang, W.S. (1996) The Haar Wavelets Operational Matrix of Integration. International Journal of Systems Science, 27, 623-628.
[28] Hooshmandasl, M.R., Heydari, M.H. and Maalek Ghaini, F.M. (2012) Numerical Solution of the One-Dimensional Heat Equation by Using Chebyshev Wavelets Method. Applied and Computational Mathematics, 1, 2168-9679.
[29] Heydari, M.H, Hooshmandasl, M.R. and Maalek Ghaini, R.M. (2014) A New Approah of Chebyshev Wavelets Method for Partial Differential Equations with Boundary Conditions of the Telegraph Type. Applied Mathematical Modelling, 38, 1597-1606.
[30] Fariborzi Araghi, M.A., Daliri, S. and Bahmanpour, M. (2012) Numerical Solution of Integro-Differential Equation by Using Chebyshev Wavelet Operational Matrix of Integration. International Journal of Mathematical Modelling & Computations, 2, 127-136.
[31] Chen, C.F. and Hsiao, C.H. (1975) Design of Piecewise Constant Gains for Optimal Control via Walsh Functions. IEEE Transactions on Automatic Control, 20, 596-603.
[32] Adibi, H. and Assari, P. (2010) Chebyshev Wavelet Method for Numerical Solution of Fredholm Integral Equations of the First Kind. Mathematical Problems in Engineering, 2010, Article ID: 138408.
[33] Mohammadi, F. (2015) A Chebyshev Wavelet Operational Method for Solving Stochastic Volterra-Fredholm integral Equations. International Journal of Applied Mathematical Research, 4, 217-227.
[34] Unsal, E., Matthai, S.K. and Blunt, M.J. (2009) Simulation of Multiphase Flow in Fractured Reservoirs Using a Fracture-Only Model with Transfer Functions. Computational Geoscience, 14, 527-538.

Copyright © 2024 by authors and Scientific Research Publishing Inc.

Creative Commons License

This work and the related PDF file are licensed under a Creative Commons Attribution 4.0 International License.