/compile/output.dvi EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS Vol. 7, No. 2, 2014, 210-229 ISSN 1307-5543 – www.ejpam.com Approximate Analytical Study of Fingero-Imbibition Phenomena of Time-Fractional Type in Double Phase Flow through Porous Media Olaniyi S. Iyiola∗, Samson Babatunde Folarin Department of Mathematics and Statistics, King Fahd University of Petroleum and Minerals, KFUPM, Dhahran, Saudi Arabia Abstract. We consider the non-linear partial differential equation of time-fractional type describing the spontaneous imbibition of water by an oil-saturated rock (double phase flow through porous media). The fact that oil and water form two immiscible liquid phases and water represents preferentially wet- ting phase are the basic assumption of this work. The Homotopy Analysis Method is used to obtain the saturation of injected water. We obtain the graphical representation of solution using MATLAB R2007b and Microsoft Excel 2010 with different fractional order (α > 0) and the comparison is made with the solution obtained in [19] using Adomian Decomposition Method when α = 1 including numerical values. 2010 Mathematics Subject Classifications: 65M99,76S05 Key Words and Phrases: Fractional derivative, Fingero-imbibition, Double phase flow in porous media, Immiscible fluid, Homotopy Analysis Method 1. Introduction In Caputo [6] and He [10], the approach used to account for the effects of changing flux is to embody the effects of memory which has to do with posing problem in terms of fractional calculus. Levy-flight type of transport is a well known diffusion process which is described by a fractional system. Motivated by this idea, we propose a fractional type Fingero-Imbibition phenomena equation in double phase flow through porous media and obtain analytical ap- proximate solution using Homotopy Analysis Method, HAM. We consider equation of form C DαT S(X , T ) = � ∂ S(X , T ) ∂ X �2 + S(X , T ) ∂ 2S(X , T ) ∂ X 2 (1) ∗Corresponding author. Email address: samuel@kfupm.edu.sa (O. S. Iyiola) http://www.ejpam.com 210 c© 2014 EJPAM All rights reserved. O. Iyiola, S. Folarin / Eur. J. Pure Appl. Math, 7 (2014), 210-229 211 where C DαT is the Caputo fractional derivative with appropriate initial condition. When there is difference in the viscosity of two flowing phases due to wetting difference, then we have fingering phenomena. The importance of this phenomenon has gained attention by various concerned fields such as geophysics, geo-hydrology, reservoir engineering etc, with little or no attention to the fractional type. Generally, for the past three decades, fractional calculus has been considered with great importance due to its various applications in fluid flow, control theory of dynamical systems, chemical physics, electrical networks, and so on. The quest of getting accurate methods for solving resulted non-linear model involving fractional order is of almost concern of many re- searchers in this field today. Various methods have been put to use successfully to obtain analytical solutions such as Adomian Decomposition Method (ADM) [2, 19, 20], Variational Iteration Method (VIM) [11, 20], Homotopy Perturbation Method (HPM) [9], and EXP-function Method [24] see also [12– 14, 23]. One of the powerful analytical approach to solving non-linear differential equations is Homotopy Analysis Method (HAM) [1, 4]. Recent works have been done using this method (HAM) to obtain analytical solutions of some differential equations given improvement on the method [3, 5, 7]. Xu et al. [25] recently applied the HAM to linear, homogeneous one and two dimensional fractional heat-like PDEs subject to the Neumann boundary conditions. Very recently, HAM was shown to be capable of solving both linear and non-linear systems of fractional partial differential equations [15]. This paper considers equation (1) subject to some appropriate initial where C DαT (· ) = ∂ α(.) ∂ Tα is a Caputo fractional differential operator. We compare the result obtained by HAM with that of ADM [19], when α = 1 to affirm the reliability of the method including numerical values. We obtain numerical results at the end with different values of α and h (h is auxiliary parameter introduced in HAM) given the effect of both parameters on solution of fingero- imbibition phenomena of fractional type in double phase flow through porous media. 2. Preliminaries Here we state necessary tools to the actualization of the aim of this paper including def- initions and some known results. This work adopts Caputo’s definition to some concepts of fractional derivatives which is a modification of the Riemann-Liouville’s definition and has the advantage of dealing properly with initial value problems. The initial conditions are given in terms of the field variables and their integer order which is the case in many physical processes. Definition 1. A real function l is said to be in the space Cµ, µ ∈ R, x > 0, if there exists a real number p(> µ) such that l(x) = x p l1(x), where l1 ∈ C[0,∞) and it is said to be in the space Cm µ if and only if l(m) ∈ Cµ, m ∈ N. Definition 2. The Riemann-Liouville’s (RL) fractional integral operator of order 0< α < 1, of a O. Iyiola, S. Folarin / Eur. J. Pure Appl. Math, 7 (2014), 210-229 212 function f ∈ L1(a, b) is given as Jα f (t) = 1 Γ(α) ∫ t 0 (t −τ)α−1 f (τ)dτ, t > 0, (2) where Γ is the Gamma function and J0 f (t) = f (t). Definition 3. The Riemann-Liouville’s (RL) fractional derivative of order 0< α < 1, of a function f is Dα0+ f (t) = DJ1−α 0+ f (t). (3) provided the right-hand side exists where D = d/d t. Definition 4. The fractional derivative in the Caputo’s sense is defined as [21], C Dα f (t) = Jn−αDn f (t) = 1 Γ(n−α) ∫ t 0 (t −τ)n−α−1 f (n)(τ)dτ, (4) where n− 1< α≤ n, n ∈ N, t > 0. Caputo’s fractional derivative also has a useful property [8] JαC Dα f (t) = f (t)− n−1 ∑ k=0 f (k)(0+) tk k! , (5) where n− 1< α≤ n. Lemma 1. Let α≥ 0, β ≥ 0 and f ∈ C L(a, b). Then Jαa Jβ f (t) = Jα+βa f (t), (6) for all t ∈ (a, b]. Lemma 2. Let t ∈ (a, b]. Then � Jαa (t − a)β � (t) = Γ(β + 1) Γ(β +α+ 1) (t − a)β+α, α¾ 0, β > 0. (7) Remark 1. From the definitions given above, we observed that the Riemann-Liouville fractional derivative of a constant function is not equal to zero while that of Caputo fractional derivative of constant function is zero. O. Iyiola, S. Folarin / Eur. J. Pure Appl. Math, 7 (2014), 210-229 213 3. Method of Solution Consider a non-linear fractional partial differential equation of the form ¨ Dαt � f (x , t) � = A � f , fx , fx x � + B � f , fx , fx x � + C(x , t) n− 1< α≤ n, t > 0 f (k)(x , 0) = gk(x) k = 0,1,2,3, . . . , n− 1, (8) where A is a linear operator and B is a non-linear operator both of which might include other fractional derivatives of order less than α and C is a known analytic function. Using (5), we obtain f (x , t) = n−1 ∑ k=0 gk(x) tk k! + JαC(x , t) + JαA � f , fx , fx x � + JαB � f , fx , fx x � , n− 1< α≤ n, t > 0. (9) 3.1. The Zeroth-Order Deformation Equation Let L denotes an auxiliary linear operator, f0(x , t) is an initial approximation of f (x , t), satisfied by the initial condition in (8) In (8), we have linear operator Dαt which in this work is different from linear operator A and we can choose it to be L(ϕ) = Dαt � ϕ � (10) with corresponding initial approximation f0(x , t) = m−1 ∑ k=0 gk(x) tk k! + JαC(x , t). (11) The non-linear operator by (8) can be defined for simplicity sake as N(ϕ) = Dαt (ϕ)− A � ϕ,ϕx ,ϕϕx x � + B � ϕ,ϕx ,ϕx x � − C(x , t). (12) We can now construct the zeroth-order deformation in the frame of Homotopy Analysis Method (HAM) [17] as (1− r)L � F(x , t; r)− f0(x , t) � = rhN (F(x , t; r)) , (13) with initial conditions: F (k)(x , 0; r) = gk(x), k = 0,1,2,3, . . . , m− 1. (14) where h 6= 0 is an auxiliary parameter, r ∈ [0,1] is the embedding parameter and F(x , 0; r) is an unknown function on the independent variables x , t, and r. We observe that ϕ = 0 is a solution of L(ϕ) = 0 where r = 0 since f0(x , t) satisfies all the initial conditions (8). O. Iyiola, S. Folarin / Eur. J. Pure Appl. Math, 7 (2014), 210-229 214 So, F(x , t; 0) = f0(x , t). (15) Also, the zeroth-order deformation equations (13) and (14) are equivalent to the original equations (8) provided r = 1 and F(x , t; 1) = f (x , t). (16) Using r, we expand, in Taylor series, F as F(x , t; r) = f0(x , t) + ∞ ∑ m=1 fm(x , t)rm. (17) where fm(x , t) = 1 m! ∂ mF(x , t; r) ∂ rm � � � � r=0 . (18) If we assume that the auxiliary linear operator L, the initial guess f0 and the auxiliary pa- rameter h are properly chosen such that the series (18) converges at r = 1, then by (16) we have f (x , t) = f0(x , t) + ∞ ∑ m=1 fm(x , t). (19) 3.2. The mth-Order Deformation Equation We consider the vector ~fm = { f0(x , t), f1(x , t), . . . , fm(x , t)}. (20) Differentiating (13) m times with respect to the (embedding) parameter r, then evaluating at r = 0 and finally dividing them by m!, we have the so called mth-order deformation equation (Lioa [16, 17]) as L � fm(x , t)−χm fm−1(x , t) � = hRm � ~fm−1 � , (21) with initial conditions f (k)m (x , 0) = 0, k = 0,1,2, ..., m− 1. (22) where Rm � ~fm−1 � = 1 (m− 1)! ∂ m−1N(F(x , t)) ∂ rm−1 � � � � r=0 (23) and χm = ¨ 0 m¶ 1 1 m> 1, (24) When we substitute (12) into (23), since A is a linear operator we get Rm � ~fm−1 � =Dαt − A � f(m−1), f(m−1)x , f(m−1)x x � O. Iyiola, S. Folarin / Eur. J. Pure Appl. Math, 7 (2014), 210-229 215 − 1 (m− 1)! ∂ m−1B(F, Fx , Fx x) ∂ rm−1 � � � � r=0 − (1−χm)C(x , t) (25) According to (10), Jα can be applied to both sides of (21) to get JαDα � fm(x , t)−χm fm−1(x , t) � = hJα � Rm � ~fm−1 �� . (26) Combining property (5) and the initial conditions in (22), we have fm(x , t) = χm fm−1(x , t) + hJα � Rm � ~fm−1 �� . (27) And, finally we will approximate the HAM solution (19) for the purpose of computation by truncated series ϕm(x , t) = m−1 ∑ k=0 fk(x , t). (28) Lemma 3 ([3]). Suppose the series f (x , t) = f0(x , t) + ∞ ∑ m=1 fm(x , t) converges, where fm(x , t) is governed by (21) with definitions (23) and (24). Then f must be a solution of (8). 3.3. Mathematical Analysis Considering a finite cylindrical piece of homogeneous porous matrix which is saturated with native liquid A surrounded completely by an impermeable surface except for an end of the cylinder labelled as the imbibition face x = 0. This end is opened to an adjacent formation of injected liquid B, see Diagram 1. The phenomenon of fingering will occur simultaneously with imbibition for a less viscous and preferentially wetting phase of liquid B which describes a one-dimensional fingero-imbibition phenomena for which the injection is started by imbibition and resulting displacement produce instabilities. We assume that the validity of Darcy’s law for the double phase flow system [22] the seepage velocities of wetting phase (vw) and the non-wetting (v0) as: ( vw = − kw µw K ∂ Pw ∂ x v0 = − k0 µ0 K ∂ P0 ∂ x , (29) where kw, k0 are relative permeability, Pw, P0 are pressure and µw, µ0 are kinematic viscosi- ties (constant) of wetting phase and non-wetting phase respectively and the permeability of homogeneous medium is K . The coordinate x is measured along the axis of the cylindrical medium, the origin being located at the imbibition face x = 0. We have that vw = −v0 (30) O. Iyiola, S. Folarin / Eur. J. Pure Appl. Math, 7 (2014), 210-229 216 Figure 1: The Fingero-Imbibition (counter-current) Phenomena in Fractured Reservoir. for a counter current flow. Hence, (29) give kw µw K ∂ Pw ∂ x + k0 µ0 K ∂ P0 ∂ x = 0 (31) Mehta [18], gives the definition of capillary pressure Pc as Pc = P0 − Pw (32) That is ∂ Pc ∂ x = ∂ P0 ∂ x − ∂ Pw ∂ x (33) Then (31) and (33) give � kw µw + k0 µ0 � ∂ Pw ∂ x + k0 µ0 ∂ Pc ∂ x = 0 (34) From here, using (34), we can write (29) explicitly as vw = K � kw µw �� k0 µ0 � ∂ Pc ∂ x � kw µw + ∂ k0 ∂ µ0 �−1 (35) For wetting phase, equation of continuity is given by φ ∂ Sw ∂ t + ∂ vw ∂ x = 0 (36) where Sw is the saturation of the wetting phase and φ is the porosity of the medium. Substituting the value of vw of (35) into (36), we obtain φ ∂ Sw ∂ t + ∂ ∂ x � K kwk0 kwµ0 + k0µw ∂ Pc ∂ x � = 0 (37) O. Iyiola, S. Folarin / Eur. J. Pure Appl. Math, 7 (2014), 210-229 217 Equation (37) is a non-linear partial differential equation that describes that the fingero- imbibition phenomenon of two immiscible fluids flow through homogeneous porous cylindrical medium with impervious bounding surface on three sides. We assume standard forms of (Scheidagger and Johnson [22]) for the analytical relation- ship between the relative permeability, phase saturation and capillary pressure phase satura- tion knowing that fictitious relative permeability is the function of displacing fluid saturation. k0 = 1− ξSw, (38) where ξ= 1.11 kw = Sw (39) Pc = βSw (40) The model we are considering involves water and viscous oil and so according to (Scheidegger [22]), we have k0kw kwµ0 + k0µw ≈ k0 µ0 = 1− ξSw µ0 = S µ0 , where S = 1− ξSw. (41) Hence, by substituting (41), (40), (39), and (38) into (37), we arrived at ∂ Sw ∂ t = kβ µ0φ ∂ ∂ x � (1− ξSw) ∂ Sw ∂ x � (42) For a dimensionless form of (42), we choose the following new variables X = x L and T = kβ φµ0 L2 t (43) to get ( ∂ Sw ∂ T = ∂ ∂ X � (1− ξSw) ∂ Sw ∂ X � Sw(X , 0) = f (X ) = e−X . (44) The choice of initial condition is due to the fact that the saturation of injected water decreases exponentially when x increases (Mehta [18]). Equation (44) with S = 1− ξSw, gives ( ∂ S ∂ T = ∂ ∂ X � S ∂ S ∂ X � S(X , 0) = f (X ) = 1− ξe−X . (45) 3.4. Main Result This present paper considers the fingero-imbibition phenomena of time-fractional type (α > 0 order) through porous media using HAM to obtain analytical solution ( DαT S = ∂ ∂ X � S ∂ S ∂ X � S(X , 0) = f (X ) = 1− ξe−X . (46) O. Iyiola, S. Folarin / Eur. J. Pure Appl. Math, 7 (2014), 210-229 218 When α= 1, we obtain (45) the usual fingero-imbibition phenomena through porous media. Equation (46) can be written as in (1) DαT S(X , T ) = � ∂ S(X , T ) ∂ X �2 + S(X , T ) ∂ 2S(X , T ) ∂ X 2 (47) We construct zeroth-order deformation from (13) ¨ (1− r)L � S̄(X , T ; r)− S̄0(X , T ) � = rhN � S̄(X , T ; r) � S̄(X , 0) = f (X ) = 1− ξe−X . (48) where N(ψ) = DαTψ−ψ 2 x −ψψx x (49) The auxiliary linear operator can be chosen as L(ψ) = DαTψ (50) with property L(ψ) = 0 for ψ = 0. (51) While the initial guess is S0(X , T ) = 1− ξe−X . (52) The high-order deformation equation from (21) as ¨ L � Sm(X , T )−χmS(m−1)(X , T ) � = hRm � ~S(m−1)(X , T ) � S(k)m (X , 0) = 0, k = 0,1, (53) where Rm � ~S(m−1) � = 1 (m− 1)! ∂ m−1N [S(X , T ; r)] ∂ rm−1 � � � � r=0 (54) Then, Rm � ~S(m−1) � can be given by Rm � ~S(m−1) � = DαT S(m−1) − m−1 ∑ i=0 Si xS(m−1−i)x − m−1 ∑ i=0 SiS(m−1−i)x x (55) Accordingly, the governing equation is as follows: Sm = χmSm−1 + hJα  DαT S(m−1) − m−1 ∑ i=0 Si xS(m−1−i)x − m−1 ∑ i=0 SiS(m−1−i)x x   (56) where m¾ 1. Using (56), we obtain the first few terms of HAM series solutions as follows S0(X , T ) = 1− ξe−X , (57) O. Iyiola, S. Folarin / Eur. J. Pure Appl. Math, 7 (2014), 210-229 219 S1(X , T ) =χ1S0 + hJα � C DαT S0 − S2 0x − S0S0x x � =− hJα � 2ξ2e−2X − ξe−X � =− 2hξ2e−2X Tα Γ(α+ 1) + hξe−X Tα Γ(α+ 1) , (58) S2(X , T ) =χ2S1 + hJα � C DαT S1 − 2S0xS1x − S0S1x x − S1S0x x � =− 2hξ2e−2X Tα Γ(α+ 1) + hξe−X Tα Γ(α+ 1) − 2h2ξ2e−2X Tα Γ(α+ 1) + h2ξe−X Tα Γ(α+ 1) + 12h2ξ2e−2X T2α Γ(2α+ 1) − 18h2ξ3e−3X T2α Γ(2α+ 1) − h2ξe−X T2α Γ(2α+ 1) . (59) The remaining terms for m= 3,4, . . . can be obtained using Mathematica or MATLAB Hence, the HAM series solution is obtained as S(X , T ) =S0(X , T ) + S1(X , T ) + S2(X , T ) + S3(X , T ) + . . . =1− ξe−X − 4hξ2e−2X Tα Γ(α+ 1) + 2hξe−X Tα Γ(α+ 1) − 2h2ξ2e−2X Tα Γ(α+ 1) + h2ξe−X Tα Γ(α+ 1) + 12h2ξ2e−2X T2α Γ(2α+ 1) − 18h2ξ3e−3X T2α Γ(2α+ 1) − h2ξe−X T2α Γ(2α+ 1) + . . . . (60) Table 1 shows the few terms HAM approximation of (1) when α= 0.10,0.25,0.50,0.75,0.90 and 1.00 with h = −0.6 and −1.0 and their respective solution plots by Matlab are given in Figures 2-3. Curves representing the effect of different values of h on the saturation S for fixed T and α are shown in Figures 4-5; Figure 6 represents the effect of different values of α on the saturation for fixed T and h. O . Iy io la , S . Fo la rin / E u r. J. P u re A p p l. M a th , 7 (2 0 1 4 ), 2 1 0 -2 2 9 2 2 0 Table 1: HAM Approximation of (1) Varying α and h α= 0.10 α= 0.25 α = 0.50 α= 0.75 α= 0.90 α= 1.00 T X h= −0.6 h= −1.0 h= −0.6 h= −1.0 h= −0.6 h= −1.0 h= −0.6 h= −1.0 h= −0.6 h= −1.0 h= −0.6 h= −1.0 0.02 0.1000 -0.7577 -3.0573 -0.1072 -0.8505 0.0624 -0.0345 0.0330 0.0213 0.0171 0.0136 0.0099 0.0083 0.2000 -0.4092 -2.0031 0.0119 -0.5400 0.1279 0.0350 0.1128 0.0939 0.1037 0.0954 0.0995 0.0946 0.3000 -0.1685 -1.2864 0.1039 -0.3206 0.1898 0.0984 0.1871 0.1629 0.1833 0.1715 0.1815 0.1740 0.4000 -0.0010 -0.8020 0.1775 -0.1629 0.2483 0.1576 0.2560 0.2280 0.2563 0.2420 0.2563 0.2471 0.5000 0.1174 -0.4758 0.2389 -0.0464 0.3038 0.2133 0.3198 0.2893 0.3232 0.3073 0.3245 0.3142 0.6000 0.2034 -0.2559 0.2921 0.0429 0.3562 0.2663 0.3787 0.3468 0.3844 0.3676 0.3867 0.3757 0.7000 0.2682 -0.1063 0.3396 0.1143 0.4056 0.3168 0.4330 0.4905 0.4404 0.4231 0.4433 0.4321 0.8000 0.3196 -0.0022 0.3833 0.1742 0.4521 0.3650 0.4830 0.4505 0.4915 0.4741 0.4949 0.4836 0.9000 0.3624 0.0730 0.4240 0.2267 0.4957 0.4108 0.5289 0.4970 0.5389 0.5210 0.5418 0.5306 1.0000 0.4010 0.1304 0.4625 0.2743 0.5365 0.4544 0.5710 0.5400 0.5806 0.5640 0.5844 0.5736 0.04 0.1000 -0.9102 -3.5498 -0.2189 -1.2668 0.0599 -0.1306 0.0522 0.0213 0.0334 0.0231 0.0232 0.0181 0.2000 -0.5073 -2.3261 -0.0613 -0.8213 0.1240 -0.0412 0.1236 0.0848 0.1131 0.0952 0.1072 0.0961 0.3000 -0.2308 -1.4956 0.0548 -0.5125 0.1833 0.0338 0.1911 0.1461 0.1873 0.1638 0.1848 0.1693 0.4000 -0.0404 -0.9360 0.1431 -0.2965 0.2388 0.0990 0.2547 0.2053 0.2561 0.2287 0.2563 0.2377 0.5000 0.0924 -0.5613 0.2129 -0.1431 0.2912 0.1578 0.3144 0.2621 0.3199 0.2899 0.3221 0.3014 0.6000 0.1868 -0.3109 0.2705 -0.0310 0.3409 0.2121 0.3702 0.3163 0.3788 0.3473 0.3825 0.3606 0.7000 0.2564 -0.1430 0.3200 0.0539 0.3881 0.2631 0.4223 0.3678 0.4331 0.4009 0.4378 0.4154 0.8000 0.3100 -0.0287 0.3642 0.1214 0.4328 0.3115 0.4707 0.4166 0.4831 0.4509 0.4885 0.4660 0.9000 0.3537 0.0518 0.4048 0.1779 0.4752 0.3577 0.5156 0.4625 0.5290 0.4873 0.5349 0.5127 1.0000 0.3914 0.1114 0.4428 0.2275 0.5152 0.4019 0.5572 0.5057 0.5711 0.5403 0.5773 0.5557 0.06 0.1000 -1.0111 -3.8724 -0.3102 -1.5916 0.0461 -0.2376 0.0649 0.0098 0.0471 0.0270 0.0354 0.0251 0.2000 -0.5721 -2.5373 -0.1209 -1.0388 0.1136 -0.1205 0.1304 0.0694 0.1209 0.0918 0.1143 0.0958 0.3000 -0.2719 -1.6318 0.0156 -0.6585 0.1738 -0.7285 0.1930 0.1271 0.1904 0.1545 0.1878 0.1634 0.4000 -0.0661 -1.0227 0.1165 -0.3958 0.2289 0.0471 0.2526 0.1829 0.2556 0.2150 0.2562 0.2276 0.5000 0.0762 -0.6159 0.1937 -0.2126 0.2803 0.1118 0.3092 0.2370 0.3166 0.2728 0.3196 0.2883 0.6000 0.1764 -0.3455 0.2555 -0.0820 0.3288 0.1695 0.3626 0.2891 0.3734 0.3278 0.3783 0.3453 0.7000 -0.2491 -0.1655 0.3070 0.0139 0.3747 0.2222 0.4129 0.3393 0.4262 0.3799 0.4324 0.3987 0.8000 0.3044 -0.0443 0.3520 0.0876 0.4184 0.2717 0.4601 0.3872 0.4752 0.4290 0.4822 0.4486 0.9000 0.3487 0.0395 0.3927 0.1473 0.4600 0.3185 0.5042 0.4329 0.5205 0.4751 0.5281 0.4949 1.0000 0.3864 0.1006 0.4305 0.1984 0.4995 0.3633 0.5454 0.4762 0.5622 0.5182 0.5702 0.5380 O. Iyiola, S. Folarin / Eur. J. Pure Appl. Math, 7 (2014), 210-229 221 (a) α= 0.1 (b) α= 0.25 (c) α= 0.5 (d) α= 0.75 (e) α= 0.9 (f) α= 1.0 Figure 2: HAM Solution for h= −0.6 and Various α O. Iyiola, S. Folarin / Eur. J. Pure Appl. Math, 7 (2014), 210-229 222 (a) α= 0.1 (b) α= 0.25 (c) α= 0.5 (d) α= 0.75 (e) α= 0.9 (f) α= 1.0 Figure 3: HAM Solution for h= −1.0 and Various α O. Iyiola, S. Folarin / Eur. J. Pure Appl. Math, 7 (2014), 210-229 223 (a) α= 0.1 (b) α= 0.25 (c) α= 0.5 (d) α= 0.75 (e) α= 0.9 (f) α= 1.0 Figure 4: Effect of h on the HAM Solution for Fixed T = 0.02 and Various α O. Iyiola, S. Folarin / Eur. J. Pure Appl. Math, 7 (2014), 210-229 224 (a) α= 0.1 (b) α= 0.25 (c) α= 0.5 (d) α= 0.75 (e) α= 0.9 (f) α= 1.0 Figure 5: Effect of h on the HAM Solution for Fixed T = 0.04 and Various α O. Iyiola, S. Folarin / Eur. J. Pure Appl. Math, 7 (2014), 210-229 225 (a) T = 0.02, h= −0.6 (b) T = 0.02, h= −1.0 (c) T = 0.04, h= −0.6 (d) T = 0.04, h= −1.0 (e) T = 0.06, h= −0.6 (f) T = 0.06, h= −1.0 Figure 6: Effect of α on the HAM Solution for Various T and h O. Iyiola, S. Folarin / Eur. J. Pure Appl. Math, 7 (2014), 210-229 226 When α = 1 and h = −0.03, the problem solved by [19], using Adomian Decomposition Method, ADM , coincides with our problem and so the solutions were compared in Table 2 and is represented in Figure 7. Table 2: Comparison of Solutions Using HAM vs. ADM T = 0.02 T = 0.04 T = 0.06 X SHAM SADM X SHAM SADM X SHAM SADM 0.1000 0.0099 0.0146 0.1000 0.0232 0.0315 0.1000 0.0354 0.0471 0.2000 0.0995 0.1052 0.1000 0.1072 0.1178 0.1000 0.1143 0.1294 0.3000 0.1815 0.1878 0.1000 0.1848 0.1969 0.1000 0.1878 0.2053 0.4000 0.2563 0.2629 0.1000 0.2563 0.2692 0.1000 0.2562 0.2751 0.5000 0.3245 0.3313 0.1000 0.3221 0.3354 0.1000 0.3196 0.3392 0.6000 0.3867 0.3834 0.1000 0.3825 0.3958 0.1000 0.3783 0.3981 0.7000 0.4433 0.4499 0.1000 0.4378 0.4510 0.1000 0.4324 0.4520 0.8000 0.4949 0.5012 0.1000 0.4885 0.5013 0.1000 0.4822 0.5012 0.9000 0.5418 0.5479 0.1000 0.5349 0.5471 0.1000 0.5281 0.5464 1.0000 0.5844 0.5902 0.1000 0.5773 0.5888 0.1000 0.5702 0.5876 (a) T = 0.02 (b) T = 0.04 (c) T = 0.06 Figure 7: Comparison of HAM solution and ADM solution for Various T REFERENCES 227 4. Discussion and conclusion The basis for many scientific and engineering applications is now depending largely on the study of physics of flow through porous media. The applications ranges from hydrologist in his study of the migration of underground water, the petroleum engineer in his study of movement of oil, gas and water through the reservoir of oil or gas field, the chemical engineer in connection with filtration processes. We also have applications in soil mechanics, ceramic engineering, oil recovery process, water purification and powder metallurgy. In the recent years, Caputo [6] and He [10], fractional order have been incorporated into existing differen- tial equations which have deep meaning and applications inexhaustible giving hope for future research work in the field of sciences and other related fields. In this work, we propose the fractional Fingero-imbibition equation to model the spon- taneous imbibition of water by an oil-saturated rock in a double phase flow through porous media. Homotopy Analysis Method (HAM) was implemented to obtain the approximate ana- lytical solution of (1). The convergence region of the series solution obtained by HAM can be adjusted and controlled by the auxiliary parameter h such adjustments are demonstrated rep- resented in Figures 4-5. It is interesting to note that different values α have significant effect on the saturation level and so it is worth investigating. The results are shown in Figure 6. Our results extend the work of Meher et al. in [19] by moving away from the classical Fingero-imbibition equation to fractional order and also using HAM (instead of ADM) which has been proven to be more accurate due to the presence of control parameter h. Several numerical results are presented including graphs to illustrates different scenario. Table 2 and Figures 7a-7c give the analysis of the comparison made between HAM and that of ADM in [19] when α = 1 and for some fixed T to show the efficiency and accuracy of the suggested method. The numerical simulations are obtained using MATLAB R2007b and Microsoft Excel 2010. ACKNOWLEDGEMENTS The financial support received from King Fahd University of Petroleum and Minerals (KFUPM) is gratefully acknowledged. References [1] S. Abbasbandy. An approximation solution of a nonlinear equation with Riemann-Liouville’s fractional derivatives by He’s variational iteration method, Journal of Computational and Applied Mathematics, vol. 207, no. 1, pp. 53-58, 2007. [2] O. Abdulaziz, I. Hashim, M. S. H. Chowdhury, and A. K. Zulkifle. Assessment of decompo- sition method for linear and nonlinear fractional differential equations, Far East Journal of Applied Mathematics, vol. 28, no. 1, pp. 95-112, 2007. [3] O. Abdulaziz, I. Hashim, and A. Saif. Series Solution of Time-Fractional PDEs by Homotopy Analysis Method, Hindawi Publishing Corporation, Differential Equations and Nonlinear Mechanics, pps. 16, 2008. REFERENCES 228 [4] A. S. Bataineh, M. S. M. Noorani, and I. Hashim. Approximate analytical solutions of systems of PDEs by homotopy analysis method, Computers and Mathematics with Appli- cations, vol. 55, no. 12, pp. 2913-2923, 2008. [5] A. S. Bataineh, M. S. M. Noorani, and I. Hashim. Approximate solutions of singular two- point BVPs by modified homotopy analysis method, Physics Letters A, vol. 372, no. 22, pp. 4062-4066, 2008. [6] M. Caputo. Diffusion of fluids in porous media with memory, Geothermics 28, 113-130, 1999. [7] G. Domairry and M. Fazeli. Homotopy analysis method to determine the fin efficiency of convective straight fins with temperature-dependent thermal conductivity, Communications in Nonlinear Science and Numerical Simulation, vol. 14, no. 2, pp. 489-499, 2009. [8] R. Gorenflo and F. Mainardi. Fractional calculus: integral and differential equations of fractional order, Fractals and Fractional Calculus in Continuum Mechanics (Udine, 1996), vol. 378 of CISM Courses and Lectures, pp. 223-276, Springer, Vienna, Austria, 1997. [9] J.-H. He. Application of homotopy perturbation method to nonlinear wave equations, Chaos, Solitons and Fractals, vol. 26, no. 3, pp. 695-700, 2005. [10] J.-H. He. Approximate analytical solution for seepage flow with fractional derivatives in porous media, Computer Methods in Applied Mechanics and Engineering. 167, 57-68, 1998. [11] J.-H. He and X.-H. Wu. Construction of solitary solution and compacton-like solution by variational iteration method, Chaos, Solitons and Fractals, vol. 29, no. 1, pp. 108-113, 2006. [12] C. H. C. Hussin and A. Kılıcman. On the Solution of Fractional Order Nonlinear Boundary Value Problems By Using Differential Transformation Method, European Journal of Pure and Applied Mathematics, vol. 4, No. 2, 174-185. 2011. [13] B. Ibis, M. Bayram and A. G. Agargün. Applications of Fractional Differential Transform Method to Fractional Differential-Algebraic Equations, European Journal of Pure and Ap- plied Mathematics, vol. 4, No. 2, 129-141. 2011. [14] O.S. Iyiola. Solving k-Fractional Hilfer Differential Equations via Combined Fractional In- tegral Transform Methods, British Journal of Mathematics and Computer Science, vol. 4, No. 10, 1427-1436, 2014. [15] H. Jafari and S. Seifi. Solving a system of nonlinear fractional partial differential equations using homotopy analysis method, Communications in Nonlinear Science and Numerical Simulation, vol. 14, no. 5, pp. 1962-1969, 2009. REFERENCES 229 [16] S.-J. Liao. An approximate solution technique not depending on small parameters: a special example, International Journal of Non-linear Mechanics, vol. 30, no. 3, pp. 371-380, 1995. [17] S.-J. Liao. Beyond Perturbation: Introduction to Homotopy Analysis Method, vol. 2 of CRC Series: Modern Mechanics and Mathematics, Chapman and Hall/CRC, Boca Raton, Fla, USA, 2004. [18] M.N. Mehta. An asymtotic expansion in fluid flow through porous media; PhD Thesis: S.G. University, Surat (India), 1978. [19] R. Meher, M. N. Mehta, and S. K. Meher. Adomian Decomposition Approach to Fingero- Imbibition Phenomena in Double Phase Flow through Porous Media, International Journal of Applied Mathematics and Mechanics. 6 (9): 34-46, 2012. [20] S. Momani and Z. Odibat. Analytical approach to linear fractional partial differential equa- tions arising in fluid mechanics, Physics Letters A, vol. 355, no. 4-5, pp. 271-279, 2006. [21] I. Podlubny. Fractional Differential Equations, Vol. 198 of Mathematics in Science and Engineering, Academic Press, San Diego, Calif, USA, 1999. [22] A. E. Scheidegger. The Physics of Flow through Porous Media, The Macmillan Co., New York, 1960. [23] V. Turut1 and N. Güzel. On solving Partial Differential Equations of Fractional Order by Using the Variational Iteration Method and Multivariate Padé Approximations, European Journal of Pure and Applied Mathematics, Vol. 6, No. 2, 147-171. 2013. [24] X.-H. Wu and J.-H. He. EXP-function method and its application to nonlinear equations, Chaos, Solitons and Fractals, vol. 38, no. 3, pp. 903-910, 2008. [25] H. Xu, S.-J. Liao, and X.-C. You. Analysis of nonlinear fractional partial differential equa- tions with the homotopy analysis method, Communications in Nonlinear Science and Nu- merical Simulation, vol. 14, no. 4, pp. 1152-1156, 2009.