www.biomathforum.org/biomath/index.php/biomath ORIGINAL ARTICLE Peristaltic transport of viscoelastic bio-fluids with fractional derivative models Emilia Bazhlekova, Ivan Bazhlekov Institute of Mathematics and Informatics Bulgarian Acad. Sci., Sofia, Bulgaria e.bazhlekova@math.bas.bg, i.bazhlekov@math.bas.bg Received: 28 October 2015, accepted: 16 May 2016, published: 26 May 2016 Abstract—Peristaltic flow of viscoelastic fluid through a uniform channel is considered under the assumptions of long wavelength and low Reynolds number. The fractional Oldroyd-B constitutive vis- coelastic law is employed. Based on models for peristaltic viscoelastic flows given in a series of papers by Tripathi et al. (e.g. Appl Math Comput. 215 (2010) 3645–3654; Math Biosci. 233 (2011) 90– 97) we present a detailed analytical and numerical study of the evolution in time of the pressure gradient across one wavelength. An analytical ex- pression for the pressure gradient is obtained in terms of Mittag-Leffler functions and its behavior is analyzed. For numerical computation the frac- tional Adams method is used. The influence of the different material parameters is discussed, as well as constraints on the parameters under which the model is physically meaningful. Keywords-Riemann-Liouville fractional deriva- tive; Mittag-Leffler function; viscoelastic flow; frac- tional Oldroyd-B constitutive model; peristalsis. I. INTRODUCTION Recently, Fractional Calculus has gained con- siderable popularity mainly due to its numerous applications in diverse fields of science and en- gineering. Fractional Calculus allows integration and differentiation of arbitrary order, not necessar- ily integer. More precisely, it deals with integro- differential operators, where the integrals are of convolution type with weakly singular power-law kernels. Extensive applications of Fractional Calculus can be found in the constitutive modeling of viscoelasticity, see [3], [4], [10], [18], [20] and the references cited there. The fractional order constitutive models (proposed in the beginning in an implicit way, see for a historical overview [19], [22], [29]) appear to be a valuable tool for describing viscoelastic properties. Unlike the clas- sical models which exhibit exponential relaxation, the models of fractional order show power-law behavior which is widely observed in a variety of experiments. They provide a higher level of ade- quacy preserving linearity and give the possibility for relatively simple description of the complex behavior of non-Newtonian viscoelastic fluids. The generalized fractional Oldroyd-B constitu- tive law belongs to the class of linear fractional models of viscoelastic fluids. It is obtained by replacing the first order derivatives in the classi- cal Oldroyd-B model by derivatives of fractional order. The corresponding constitutive equation in the one-dimensional case is given by (1 + λα1 D α t ) τ(t) = η ( 1 + λ β 2D β t ) ε̇(t), (1) where τ(t) is the shear stress, ε(t) - shear strain, Citation: Emilia Bazhlekova, Ivan Bazhlekov, Peristaltic transport of viscoelastic bio-fluids with fractional derivative models, Biomath 5 (2016), 1605151, http://dx.doi.org/10.11145/j.biomath.2016.05.161 Page 1 of 12 http://www.biomathforum.org/biomath/index.php/biomath http://dx.doi.org/10.11145/j.biomath.2016.05.161 E. Bazhlekova et al., Peristaltic transport of viscoelastic bio-fluids ... the over-dot denotes the first time derivative, η > 0 is the dynamic viscosity of the fluid, λ1,λ2 ≥ 0 are parameters related to the relaxation and retardation times, respectively, and Dαt and D β t are fractional Riemann-Liouville derivatives in time of orders α and β, where 0 < α ≤ 1, 0 < β ≤ 1. The gener- alized Oldroyd-B model (1) encompasses a large class of fluids: Newtonian fluid (λ1 = λ2 = 0), fractional second grade fluid (λ1 = 0, λ2 > 0), fractional Maxwell fluid (λ2 = 0, λ1 > 0). In [23] a very good fit with experimental data is achieved for the fractional Oldroyd-B model. Unidirectional flows of viscoelastic fluids with the fractional Oldroyd-B constitutive law are studied in [5], [11], [13], [17], [21], to mention only few of many recent publications. The transportation of many biophysical fluids is controlled by a special mechanism called peristal- sis. The mechanism includes involuntary periodic contraction followed by relaxation and expansion of the ducts through which the fluids pass. It is inducted by the propagation of electrochemically generated waves along the vessels containing flu- ids. Examples from physiology where the peri- staltic transport is prevalent are the movement of chyme in the small intestine, transport of bile in bile ducts, transport of lymph in the lymphatic vessels, etc. The complex physical nature of peri- staltic flows of non-Newtonian fluids stimulated significant attention in the applied mathematics and engineering sciences research communities. For recent research on this topic we refer to [1], [2], [8], [9], [15], [28]. Fractional derivative models for peristaltic trans- port of viscoelastic fluids are derived in a series of papers by Tripathi et al. (e.g. [24], [25], [26], [27]) As noted in [26] such models are appropriate for describing the chyme movement in the small intestine, by considering the gastric chyme as a viscoelastic fluid. In [26] and [27] peristaltic transport through a cylindrical tube of fractional Oldroyd-B fluid is studied, with 0 < α ≤ β ≤ 1. In [26] inclined tube is considered and in [27] wall slip conditions are assumed. In [24], [25] the particular case of a fractional Maxwell model (λ2 = 0) is considered. For the numerical compu- tations the Adomian decomposition and homotopy analysis methods are used in the aforementioned articles. In the present work, peristaltic flow of viscoelas- tic fluid through a uniform channel is considered under the assumptions of long wavelength and low Reynolds number. The viscoelastic properties of the fluid are governed by the fractional Oldroyd-B constitutive equation. We employ the model pro- posed in [24] for the particular case of fractional Maxwell fluid (λ2 = 0) and generalize it in a straightforward way to cover also the case λ2 6= 0. Since the considered model is non-stationary in nature and contains parameters which are time- related such as the orders of the fractional time derivatives α and β, the relaxation and retardation times λ1, λ2, it is natural to study the time evo- lution of the physical quantities described by the model and the influence of the different parameters on this evolution. Our main contribution is a detailed analytical and numerical analysis of the time evolution of the pressure gradient across one wavelength in the peristaltic flow. To the best of our knowledge, this issue has not been discussed before in the general case λ2 6= 0. An explicit expression for the pressure gradient in terms of the Mittag-Leffler functions is derived and its behavior is studied. For the numerical computations a technique based on the fractional Adams method [6], [7] is used. Re- sults of several numerical examples are given and the influence of the different material parameters is discussed as well as constraints on the parameters under which the model is physically meaningful. The rest of the paper is organized as follows. In Section II the equation for the pressure gradient is derived. In Section III an analytical representation for the pressure gradient is obtained in terms of Mittag-Leffler functions and its behavior is ana- lyzed. In Section IV the numerical method used for the computations is described. The obtained numerical results are given in Section V and the influence of the different material parameters is discussed. Section VI contains conclusions. Some Biomath 5 (2016), 1605151, http://dx.doi.org/10.11145/j.biomath.2016.05.161 Page 2 of 12 http://dx.doi.org/10.11145/j.biomath.2016.05.161 E. Bazhlekova et al., Peristaltic transport of viscoelastic bio-fluids ... basic definitions and results from Fractional Cal- culus which are used in this work are summarized in an Appendix. II. MATHEMATICAL MODEL In this section we derive the equation for the pressure gradient in a peristaltic flow of a vis- coelastic fluid with fractional Oldroyd-B constitu- tive model. This equation was originally proposed in [24] for the case λ2 = 0. The necessary gener- alizations to the case λ2 6= 0 are straightforward. Here we give for completeness the main steps in the derivation. For further details we refer to [24] and the related works [17], [21], [25], [26], [27]. The fundamental equations governing the un- steady motion of an incompressible viscoelastic fluid are the continuity equation: ∇·V = 0 and the general Cauchy momentum equation: ρ ( ∂V ∂t + (V ·∇)V ) = ∇·σ (2) where V is the velocity vector, ρ - fluid density, σ - Cauchy stress tensor: σ = −pI + τ, (3) where p is the pressure and τ is the shear stress tensor, which for a viscoelastic fluid with the generalized fractional Oldroyd-B model satisfies the equation( 1 + λα1 Dα Dtα ) τ = η ( 1 + λ β 2 Dβ Dtβ ) A1. (4) Here η > 0 is the dynamic viscosity of the fluid, λ1 and λ2 are parameters related to relaxation and retardation times, respectively, satisfying (see [23]) λ1 ≥ λ2 ≥ 0, (5) α and β are fractional parameters, 0 < α ≤ 1, 0 < β ≤ 1, (6) A1 is the first Rivlin-Ericksen tensor given by A1 = ∇V + (∇V)T (7) and Dγ Dtγ denotes the upper convected time deriva- tive Dγτ Dtγ = D γ t τ +(V ·∇)τ − (∇V) ·τ −τ · (∇V) T , where Dγt is the Riemann-Liouville fractional derivative, see (51) in the Appendix for the defi- nition. It is assumed that the relevant Reynolds number is small enough for inertial effects to be negligible and the wavelength to diameter ratio is large enough for the pressure to be considered uniform over the cross-section of the channel. As in [24] we consider a uniform horizontal two-dimensional channel with h being the transverse displacement of the walls. Denote by x the axis along the channel and by y the transverse coordinate. Let u be the velocity of the flow in the direction of the channel. First, the above equations are rewritten in di- mensionless form (for details see [24]). Here, for simplicity we keep the same notations. According to the assumption of low Reynolds number, one obtains inserting (3) in the momentum equation (2): ∂p ∂x = ∂τxy ∂y , (8) ∂p ∂y = 0. (9) On the other hand, the constitutive equation (4) gives: (1 + λα1 D α t ) τxy = ( 1 + λ β 2D β t ) ∂u ∂y . (10) Applying the operator (1 + λα1 D α t ) to both sides of (8) the following identity is deduced: (1 + λα1 D α t ) ∂p ∂x = (1 + λα1 D α t ) ∂τxy ∂y (11) Differentiating with respect to y both sides of equation (10) one gets (1 + λα1 D α t ) ∂τxy ∂y = ( 1 + λ β 2D β t ) ∂2u ∂y2 . (12) Biomath 5 (2016), 1605151, http://dx.doi.org/10.11145/j.biomath.2016.05.161 Page 3 of 12 http://dx.doi.org/10.11145/j.biomath.2016.05.161 E. Bazhlekova et al., Peristaltic transport of viscoelastic bio-fluids ... From (11) and (12) one deduces the following equation for the pressure gradient ∂p ∂x (1 + λα1 D α t ) ∂p ∂x = ( 1 + λ β 2D β t ) ∂2u ∂y2 . (13) Boundary conditions for the velocity are given by ∂u ∂y ∣∣∣∣ y=0 = 0, u|y=h = 0. (14) In the following integrations we essentially use the fact that the pressure does not depend on the transverse coordinate y, see (9), i.e. p = p(x,t). Integrating Eq. (13) with respect to y twice and using boundary conditions (14) one obtains suc- cessively (1 + λα1 D α t ) ∂p ∂x y = ( 1 + λ β 2D β t ) ∂u ∂y . (15) y2 −h2 2 (1 + λα1 D α t ) ∂p ∂x = ( 1 + λ β 2D β t ) u. (16) Denote by Q the volumetric flow rate Q =∫h 0 u dy. Then after one more integration (16) implies − h3 3 (1 + λα1 D α t ) ∂p ∂x = ( 1 + λ β 2D β t ) Q. (17) Following [24] it is assumed that the wall of the channel undergoes contraction and relaxation given by h = 1 −φ cos2(πx), (18) where φ is the amplitude of the wave. The trans- formations between the wave and the laboratory frames (an established procedure in peristaltic fluid dynamics, see [14]) are given in dimensionless form by X = x−t, Y = y, U = u−1, θ = Q−h. (19) Let Q denotes the averaged volumetric flow rate Q = ∫ 1 0 Q dt. (20) Using the following relation from [14] Q = θ + 1 − φ 2 = Q−h + 1 − φ 2 , (21) Eq. (17) gives (1 + λα1 D α t ) ∂p ∂x = ( 1 + λ β 2D β t ) A, (22) where A = − 3 h3 ( Q + h− 1 + φ/2 ) . (23) For further details on this derivation we refer to [24]. Let us emphasize that the function A defined by (23) does not depend on time, but it depends on the spatial variable x via the peristaltic wave parameters h,φ and Q. Therefore we write A = A(x). III. ANALYTICAL PROPERTIES OF THE PRESSURE GRADIENT FUNCTION Since A(x) is independent of t equation (22) can be rewritten in the form (1 + λα1 D α t ) ∂p ∂x = ( 1 + λ β 2 t−β Γ(1 −β) ) A(x). (24) Here we have used the identity D β t 1 = t−β Γ(1 −β) , 0 < β < 1, (25) obtained by applying the definition (51) of the Riemann-Liouville fractional derivative and iden- tity (50). Eq. (24) implies that the pressure gradient can be expressed as follows ∂p ∂x = A(x)y(t), (26) where the function y(t) is a solution of the equa- tion (1 + λα1 D α t ) y(t) = 1 + λ β 2 t−β Γ(1 −β) . (27) Therefore, the time evolution of the pressure gra- dient is determined by the behavior of the function y(t). In the present work our study is limited to the properties of this function. In what follows we assume λ1 6= 0. Let us rewrite (27) in the form of the following fractional order equation Dαt y(t) = − 1 λα1 y(t) + F(t), (28) Biomath 5 (2016), 1605151, http://dx.doi.org/10.11145/j.biomath.2016.05.161 Page 4 of 12 http://dx.doi.org/10.11145/j.biomath.2016.05.161 E. Bazhlekova et al., Peristaltic transport of viscoelastic bio-fluids ... where F(t) = 1 λα1 ( 1 + λ β 2 t−β Γ(1 −β) ) . (29) Suppose the physically reasonable initial condition ∂p ∂x ∣∣∣∣ t=0 < ∞. (30) Therefore, from (26), y(0) < ∞ and thus lim t→0 (J1−αt y)(t) = 0, (31) where J1−αt is a Riemann-Liouville fractional in- tegral, see (47) in the Appendix. According to (63) the solution of equation (28) is given by y(t) = ∫ t 0 τα−1Eα,α ( − τα λα1 ) F(t− τ) dτ, (32) where Eα,α(·) denotes the Mittag-Leffler function (see (55) for the definition). From (32), (29) and the definition (47) of the fractional integration operator we deduce the representation y(t) = 1 λα1 ( J1t + λ β 2J 1−β t )( tα−1Eα,α ( − tα λα1 )) , which by the use of identities (61) and (59) implies that the function y(t) can be expressed in terms of the Mittag-Leffler functions as follows y(t) = 1 −Eα,1 ( − tα λα1 ) (33) + λ β 2 λα1 tα−βEα,α+1−β ( − tα λα1 ) . Inserting (55) into (33) one obtains the following series expansion y(t) = 1 − ∞∑ k=0 (−1)ktαk λαk1 Γ(αk + 1) (34) + λ β 2 ∞∑ k=0 (−1)ktαk+α−β λ α(k+1) 1 Γ(αk + α−β + 1) . In the particular case α = β (33) reduces to y(t) = 1 − ( 1 − ( λ2 λ1 )α) Eα,1 ( − tα λα1 ) , (35) while for λ2 = 0 it gives y(t) = 1 −Eα,1 ( − tα λα1 ) . (36) Based on the obtained expressions we study the behavior of the function y(t). Recall the restric- tions on the parameters λ1 ≥ λ2 ≥ 0, 0 < α ≤ 1 and 0 < β ≤ 1. In the simplest case λ2 = 0, based on rep- resentation (36) and the properties of Mittag- Leffler function, we easily infer that y(t) is a monotonically increasing function with y(0) = 0 and y(+∞) = 1. Moreover, for small times t the function y(t) grows faster when α is smaller, while for large t it grows faster (and approaches the value 1) when α is larger. This behavior can be seen also on Fig. 1. Qualitatively similar behavior can be deduced from the representation (35) for the case α = β, taking into account that λ1 ≥ λ2 (see also Fig. 5 and Fig. 6). However, there is one essential differ- ence: in this case y(t) does not vanish at t = 0, more precisely, (35) implies y(0) = ( λ2 λ1 )α , α = β. (37) To find the asymptotic behavior of y(t) for t → 0 we take the first terms in the series representation (34) and obtain for λ2 6= 0: y(t) = λ β 2 λα1 tα−β Γ(1 + α−β) + O ( tmin{α,2α−β} ) , (38) and for λ2 = 0: y(t) = 1 λα1 tα Γ(1 + α) + O ( t2α ) . (39) Therefore, for λ2 = 0 as well as for λ2 6= 0 and α > β the function y(t) vanishes for t → 0. If λ2 6= 0 and α = β then the initial value y(0) is as in (37). However, if λ2 6= 0 and α < β, then the asymptotic expansion (38) implies that the func- tion y(t) has a weak singularity as t → 0 (see also Fig. 3). This contradicts initial condition (30) and raises the question whether the model is physically correct in this case. Biomath 5 (2016), 1605151, http://dx.doi.org/10.11145/j.biomath.2016.05.161 Page 5 of 12 http://dx.doi.org/10.11145/j.biomath.2016.05.161 E. Bazhlekova et al., Peristaltic transport of viscoelastic bio-fluids ... Concerning large time behavior, asymptotic ex- pansion (56) in the Appendix implies y(t) = 1 −λα1 t−α Γ(1 −α) + λα2 t−β Γ(1 −β) + O ( t−min{α+β,2α} ) , t → +∞. (40) This asymptotic expansion is valid in all of the considered cases. Therefore, in all cases lim t→+∞ y(t) = 1. (41) This can also be observed on the figures. IV. NUMERICAL METHOD Let us note first that the explicit representation (34) derived from the series expansions of the Mittag-Leffler functions is appropriate for numer- ical computation of the function y(t) only for sufficiently small times. In [24], [25], [26], [27] two semi-numerical techniques are used for the solution of Eq. (27): Adomian decomposition method (ADM) and ho- motopy analysis method (HAM). These two meth- ods give a series of functions, which first terms are used for the numerical computation of y(t). It appears that the obtained approximations of y(t) by these two methods (for the chosen parameters in HAM ~ = −1 and p0 = 0) are the same as if we take the first terms of the series in (34). Therefore, it can be expected that the numerical techniques proposed in these studies retain the aforementioned disadvantage of using the series expansion (34) for numerical computation: they work only for sufficiently small times. In contrast, the numerical technique used in the present work is appropriate for all times. For numerical computation of the function y(t) we use an algorithm based on its representation as a solution of an integral equation. Applying the operator Jαt to both sides of equation (28) and using (31), (52), and the semi-group property (49), we obtain that y(t) satisfies the following integral equation y(t) = − 1 λα1 ∫ t 0 (t− τ)α−1 Γ(α) y(τ) dτ +H(t), (42) where H(t) = 1 λα1 ( tα Γ(α + 1) + λ β 2 tα−β Γ(α−β + 1) ) . (43) Equation (42) is used here for the numerical computation of the function y(t), applying the so- called fractional Adams method, originally pro- posed and analyzed by Diethelm et al. [6], [7]. This is a predictor-corrector method in which as predictor the fractional Adams-Bashforth method is used and as corrector the fractional Adams- Moulton method. For completeness, here we give the numerical scheme. To find a numerical solution of Eq. (42) in the time interval t ∈ [0,T] consider a uniform grid {tj = jh,j = 0, 1, ...,N} with some integer N and h = T/N. Denote by yj the approximation for y(tj). For the sake of brevity the notation λ = −λ−α1 is used. The predictor yPk+1 is determined by the formula yPk+1 = λ Γ(α) k∑ j=0 bj,k+1yj + H(tk), (44) where bj,k+1 = hα α ((k + 1 − j)α − (k − j)α). (45) The corrector formula is given by yk+1 = λ Γ(α)   k∑ j=0 aj,k+1yj + ak+1,k+1y P k+1   + H(tk+1), where aj,k+1 = hα α(α + 1) Aj,k+1 (46) and Aj,k+1 are defined by Aj,k+1=   kα+1−(k−α)(k+ 1)α if j = 0, (k − j + 2)α+1 + (k − j)α+1 −2(k − j + 1)α+1 if 1 ≤j ≤k, 1 if j = k + 1. Using this numerical algorithm, the function y(t) is computed for several values of the parameters. The performed numerical experiments indicate Biomath 5 (2016), 1605151, http://dx.doi.org/10.11145/j.biomath.2016.05.161 Page 6 of 12 http://dx.doi.org/10.11145/j.biomath.2016.05.161 E. Bazhlekova et al., Peristaltic transport of viscoelastic bio-fluids ... that this method is fast and stable. Although, due to the calculation of the integral, it is more time consuming for larger T , the numerical experiments with T = 10 and T = 100 indicate that the method works sufficiently well also for large time intervals. V. NUMERICAL RESULTS AND DISCUSSION In this section we discuss some results for the function y(t), obtained by the numerical technique described above. Recall that the graphs of y(t) rep- resent the time evolution of the pressure gradient. On Fig. 1 and Fig. 2 plots of the function y(t) are presented for λ2 = 0, which corresponds to the case of fractional Maxwell model, considered in [24]. Comparing these two figures to [24], Fig. 1 and Fig. 2, we observe the same behavior (for better comparison we have chosen the same values for the parameter α as in [24]). The time profiles on Fig. 1 exhibit increasing pressure gradient with time as for smaller α it increases faster for small t and slower for large t, whereas for larger α the situation is opposite: it increases slower for small t and faster for large t. This is in agreement with the theoretical observations in Section III based on the analytical representation (36). Unlike the figures in [24], where only the time interval t ∈ [0, 1] is considered, on Fig. 1 we also give plots for t > 1, which reveal that the pressure gradient does not grow infinitely with time and approaches a certain value (A(x)). This confirms the theoretical observations in Section III. On Fig. 2, where the influence of the relaxation time λ1 is illustrated, we see that the pressure gradient is smaller for larger values of λ1. Therefore this parameter resists the movement of the flow. Figures 3–6 correspond to the general case λ2 6= 0. On Fig. 3 the behavior of the pressure gradient function for α < β is illustrated. It is seen that the pressure gradient has a singularity at t = 0 (y(0) = +∞). This was also observed in Section III. To the best of our knowledge, this feature of the model has not been discussed before and raises the question of its physical adequacy. In the works 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 t y λ 1 =1, λ 2 =0 α=1/3 α=1/2 α=2/3 α=1 0 1 2 3 4 5 6 7 8 9 10 0 0.2 0.4 0.6 0.8 1 t y λ 1 =1, λ 2 =0 α=1/3 α=1/2 α=2/3 α=1 Fig. 1. Time profile of the pressure gradient for A = 1, λ1 = 1, λ2 = 0, and various values of α. Time interval [0,1] (above) and [0,10] (below). [26] and [27], where the case α ≤ β is considered, this issue has not been addressed. On Fig. 4 the case α > β is illustrated. Com- paring Fig. 3 and Fig. 4 it is seen that the behavior for α > β is qualitatively different from those for α < β. For small times the pressure gradient is monotonically decreasing for α < β (Fig. 3) and monotonically increasing for α > β (Fig. 4). However, this monotonic behavior is not retained for all t. The influence of the fractional parameters α and β observed on both figures is as follows. The effect of the fractional parameter α for small times is opposite to that for large times. The same holds for the fractional parameter β. In addition, the effects of the parameters α and β are found to be opposite to each other. Plots for the case α = β are given on Fig. 5 and Fig. 6. The influence of the relaxation time Biomath 5 (2016), 1605151, http://dx.doi.org/10.11145/j.biomath.2016.05.161 Page 7 of 12 http://dx.doi.org/10.11145/j.biomath.2016.05.161 E. Bazhlekova et al., Peristaltic transport of viscoelastic bio-fluids ... 0 1 2 3 4 5 6 7 8 9 10 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 t y λ 2 =0, α=0.5 λ 1 =0.4 λ 1 =1.0 λ 1 =3.0 λ 1 =10. Fig. 2. Time profile of the pressure gradient for A = 1, λ2 = 0, α = 0.5 and various values of λ1. 0 1 2 3 4 5 6 7 8 9 10 0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 t y λ 1 =λ 2 =1, α=0.2 β=0.4 β=0.6 β=0.8 β=1.0 0 1 2 3 4 5 6 7 8 9 10 0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 t y λ 1 =λ 2 =1, β=0.8 α=0.1 α=0.3 α=0.5 α=0.7 Fig. 3. Time profiles of the pressure gradient for α < β, A = 1 and λ1 = λ2 = 1. λ1 and retardation time λ2 observed on Fig. 5 is as follows. The pressure gradient increases with the retardation time λ2 whereas it decreases with the relaxation time λ1. On Fig. 6 the effect of 0 1 2 3 4 5 6 7 8 9 10 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 t y λ 1 =λ 2 =1, α=0.8 β=0.1 β=0.3 β=0.5 β=0.7 0 1 2 3 4 5 6 7 8 9 10 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 t y λ 1 =λ 2 =1, β=0.2 α=0.4 α=0.6 α=0.8 α=1.0 Fig. 4. Time profiles of the pressure gradient for α > β, A = 1 and λ1 = λ2 = 1. the fractional parameter α is examined. In order to capture the peculiarities of the function y(t) a larger time interval is considered t ∈ [0, 100]. The influence of the fractional parameter resemble those observed on Fig. 1 in the case of fractional Maxwell model. This is in agreement with the similarity in the explicit expressions (35) and (36). VI. CONCLUSIONS Employing the mathematical tools of Fractional Calculus we study in this work the time evolution of the pressure gradient in a viscoelastic peristaltic flow with fractional Oldroyd-B constitutive model. The analysis of the effect of different parameters shows that for α < β there is an unphysical singularity. This means that from the previously considered in [26] and [27] range 0 < α ≤ β ≤ 1 only for α = β the model is physically meaningful. This is also in agreement with the Biomath 5 (2016), 1605151, http://dx.doi.org/10.11145/j.biomath.2016.05.161 Page 8 of 12 http://dx.doi.org/10.11145/j.biomath.2016.05.161 E. Bazhlekova et al., Peristaltic transport of viscoelastic bio-fluids ... 0 1 2 3 4 5 6 7 8 9 10 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 t y α=β=0.5, λ 2 =0.2 λ 1 =0.4 λ 1 =1.0 λ 1 =3.0 λ 1 =10 0 1 2 3 4 5 6 7 8 9 10 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 t y α=β=0.5, λ 1 =20 λ 2 =0.4 λ 2 =1.0 λ 2 =3.0 λ 2 =10 Fig. 5. Time profiles of the pressure gradient for α = β = 0.5, A = 1 and various values of λ1 and λ2, λ1 > λ2. statement in [30], that the Oldroyd-B constitutive law is thermodynamically compatible only if the fractional orders α and β coincide and λ1 ≥ λ2. In both cases of physical interest: λ2 = 0, 0 < α ≤ 1 (fractional Maxwell model) and λ1 ≥ λ2 > 0, 0 < α = β ≤ 1 (thermodynam- ically compatible Oldroyd-B model) the pressure gradient across one wavelength is monotonically increasing with time and approaches a certain stationary value. The same qualitative behavior will hold for the pressure rise and friction force. The technique used in this work can be applied to more general fractional derivative viscoelastic models of peristaltic transport, such as models with more complicated geometry (non-uniform, cylindrical, inclined channels), flows with slip ef- fects, flows in porous media, etc. 0 1 2 3 4 5 6 7 8 9 10 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 t y λ 1 =10, λ 2 =1.0 α=β=0.2 α=β=0.4 α=β=0.6 α=β=0.8 0 10 20 30 40 50 60 70 80 90 100 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 t y λ 1 =10, λ 2 =1.0 α=β=0.2 α=β=0.4 α=β=0.6 α=β=0.8 Fig. 6. Time profiles of the pressure gradient for α = β, A = 1, λ1 = 10, λ2 = 1. Time interval [0,10] (above) and [0,100] (below). ACKNOWLEDGMENTS The work is partially supported by Grant DFNI-I02/9 from the Bulgarian National Science Fund; and the bilateral research project between Bulgarian and Serbian academies of sciences (2014-2016) ”Mathematical modeling via integral- transform methods, partial differential equations, special and generalized functions, numerical anal- ysis”. APPENDIX Here we summarize some facts from Fractional Calculus, for details see [12], [16]. The fractional order Riemann-Liouville integral Jαt is defined by Jαt f(t) = ∫ t 0 ωα(t− τ)f(τ) dτ, α > 0, (47) Biomath 5 (2016), 1605151, http://dx.doi.org/10.11145/j.biomath.2016.05.161 Page 9 of 12 http://dx.doi.org/10.11145/j.biomath.2016.05.161 E. Bazhlekova et al., Peristaltic transport of viscoelastic bio-fluids ... where ωα(t) = tα−1 Γ(α) , α > 0, t > 0. (48) Here Γ(·) is the Gamma function. Basic properties of this function are Γ(1) = 1, Γ(α + 1) = αΓ(α). Therefore ω1(t) ≡ 1. The operators of fractional integration satisfy the semi-group property: Jαt J β t = J α+β t , α,β > 0, (49) or, equivalently, Jαt ωβ = ωα+β, α,β > 0. (50) The Riemann-Liouville fractional derivative Dαt of order α ∈ (0, 1] is defined by D1t = d/dt and Dαt = D 1 t J 1−α t . (51) The Riemann-Liouville fractional derivatives and integrals are related via the identities: Dαt J α t f = f, α > 0, and Jαt D α t f = f − (J 1−α t f)(0)ωα(t), α ∈ (0, 1). (52) Application of the Laplace transform L{f(t)}(s) = f̂(s) = ∫ ∞ 0 e−stf(t) dt to the operator of fractional integration gives L{Jαt f}(s) = s −αf̂(s), α > 0, (53) which implies the following identity for the Riemann-Liouville fractional derivative of order α ∈ (0, 1): L{Dαt f}(s) = s αf̂(s) − (J1−αt f)(0). (54) Denote as usual by Eα,β(·) the two-parameter Mittag-Leffler function Eα,β(z) = ∞∑ k=0 zk Γ(αk + β) . (55) For α ∈ (0, 2),β > 0, the Mittag-Leffler function has the following asymptotic expansion as t → +∞ Eα,β(−t) = − N−1∑ k=1 (−t)−k Γ(β −αk) + O(t−N ). (56) An important particular case is E1,1(−t) = exp(−t) and some properties of the function Eα,1(−t) for 0 < α < 1 resemble the behavior of the exponential function: Eα,1(−t) is monotonically decreasing with Eα,1(0) = 1, Eα,1(−∞) = 0. However, unlike the fast exponential decay of exp(−t) for large t, the Mittag-Leffler function admits a slow algebraic decay, which is slower for smaller α. At t = 0 the opposite picture is observed: the Mittag-Leffler function admits a fast decay ( ddtEα,1(−t) →∞ for t → 0, see (60)), and this decay is faster for smaller α. Recall the Laplace transform pairs L{ωα(t)}(s) = s−α, (57) L { tβ−1Eα,β(λt α) } (s) = sα−β sα −λ . (58) The following identity is often useful J γ t ( tβ−1Eα,β(λt α) ) = tβ+γ−1Eα,β+γ(λt α), (59) where α,β,γ,t > 0. It can be proven by applying Laplace transform and using (53) and (58). Another useful property is the following d dt Eα,1(λt α) = λtα−1Eα,α(λt α). (60) It can be deduced again by applying Laplace transform and using (58) or directly from the series representation (55) of the Mittag-Leffler function. Integrating (60) and using that Eα,1(0) = 1 we obtain for α > 0 and t > 0 J1t ( tα−1Eα,α(λt α) ) = − 1 λ (1 −Eα,1(λtα)) (61) Let f(t) be an integrable on (0,T) function and let α ∈ (0, 1). Then the differential equation of fractional order Dαt y(t) = λy(t) + f(t), t > 0, (62) Biomath 5 (2016), 1605151, http://dx.doi.org/10.11145/j.biomath.2016.05.161 Page 10 of 12 http://dx.doi.org/10.11145/j.biomath.2016.05.161 E. Bazhlekova et al., Peristaltic transport of viscoelastic bio-fluids ... has a unique solution given by y(t) = y0t α−1Eα,α(λt α) (63) + ∫ t 0 τα−1Eα,α(λτ α)f(t− τ) dτ, where y0 = limt→0 J 1−α t y. This result can be found in [16], p. 137. The easiest way to prove it is by applying Laplace transform. REFERENCES [1] N.S. Akbar, M. Raza, R. Ellahi, Peristaltic flow with thermal conductivity of H2O + Cu nanofluid and entropy generation, Results Phys. 5 (2015) 115–124. http://dx.doi.org/10.1016/j.rinp.2015.04.003 [2] N.S. Akbar, M. Raza, R. Ellahi, Influence of induced magnetic field and heat flux with the suspension of carbon nanotubes for the peristaltic flow in a permeable channel, J Magn Magn Mater. 38 (2015) 405–415. http://dx.doi.org/10.1016/j.jmmm.2014.12.087 [3] T.M. Atanacković, S. Pilipović, B. Stanković, D. Zorica, Fractional Calculus with Applications in Mechanics: Vibrations and Diffusion Processes, John Wiley & Sons, 2014. [4] R.L. Bagley, P.J. Torvik, On the fractional calculus model of viscoelastic behavior, J Rheol. 30 1 (1986) 133-155. http://dx.doi.org/10.1122/1.549887 [5] E. Bazhlekova, I. Bazhlekov, Viscoelastic flows with fractional derivative models: computational approach via convolutional calculus of Dimovski, Fract Calc Appl Anal. 17 4 (2014) 954–976. http://dx.doi.org/10.2478/s13540-014-0209-x [6] K. Diethelm, N.J. Ford, A.D. Freed, Detailed error analysis for a fractional Adams method, Numer Algorithms 36 1 (2004) 31–52. http://dx.doi.org/10.1023/B:NUMA.0000027736.85078. be [7] K. Diethelm, N.J. Ford, A.D. Freed, A predictor-corrector approach for the numerical solution of fractional differ- ential equations, Nonlinear Dynam. 29 1 (2002) 3–22. http://dx.doi.org/10.1023/A:1016592219341 [8] R. Ellahi, F. Hussain, Simultaneous effects of MHD and partial slip on peristaltic flow of Jeffery fluid in a rectangular duct, J Magn Magn Mater. 393 (2015) 284– 292. http://dx.doi.org/10.1016/j.jmmm.2015.05.071 [9] R. Ellahi, S.U. Rahman, S. Nadeem, K. Vafai, The blood flow of Prandtl fluid through a tapered stenosed arteries in permeable walls with magnetic field, Commun Theor Phys. 63 03 (2015) 353–358. http://dx.doi.org/10.1088/0253-6102/63/3/353 [10] M. Fabrizio, Fractional rheological models for thermo- mechanical systems. Dissipation and free energies, Fract Calc Appl Anal. 17 1 (2014) 206–223. http://dx.doi.org/10.2478/s13540-014-0163-7 [11] C. Fetecau, M. Jamil, C. Fetecau, D. Vieru, The Rayleigh-Stokes problem for an edge in a generalized Oldroyd-B fluid, Z Angew Math Phys. 60 5 (2009), 921- 933. [12] R. Gorenflo, F. Mainardi, Fractional calculus: inte- gral and differential equations of fractional order. In: A. Carpinteri, F. Mainardi (Eds.), Fractals and Frac- tional Calculus in Continuum Mechanics, Springer- Verlag, Wien/New York (1997) 223–276. [13] J. Hristov, Integral-balance solution to the Stokes first problem of a viscoelastic generalized second grade fluid, Therm Sci. 16 (2012), 395–410. http://dx.doi.org/10.2298/TSCI110401077H [14] M. Jaffrin, A. Shapiro, Peristaltic pumping, Annu Rev Fluid Mech. 3 (1971) 13-36. http://dx.doi.org/10.1146/annurev.fl.03.010171.000305 [15] M.S. Kandelousi, R. Ellahi, Simulation of ferrofluid flow for magnetic drug targeting using Lattice Boltzmann method, Z Naturforsch. A 70 (2015) 115–124. http://dx.doi.org/10.1515/zna-2014-0258 [16] A.A. Kilbas, H.M. Srivastava, J.J. Trujillo, Theory and applications of fractional differential equations, North-Holland Mathematics studies, Amsterdam: Else- vier, 2006. [17] Y. Liu, L. Zheng, X. Zhang, Unsteady MHD Couette flow of a generalized Oldroyd-B fluid with fractional derivative, Comput Math Appl. 61 2 (2011) 443–450. http://dx.doi.org/10.1016/j.camwa.2010.11.021 [18] F. Mainardi, Fractional Calculus and Waves in Linear Viscoelasticity, London: Imperial College Press, 2010. [19] F. Mainardi, An historical perspective on fractional calculus in linear viscoelasticity, Fract Calc Appl Anal. 15 4 (2012), 712–717. http://dx.doi.org/10.2478/s13540-012-0048-6 [20] S.P. Näsholm, S. Holm, On a fractional Zener elastic wave equation, Fract Calc Appl Anal. 16 1 (2013) 26– 50. http://dx.doi.org/10.2478/s13540-013--0003-1 [21] H. Qi, M. Xu, Some unsteady unidirectional flows of a generalized Oldroyd-B fluid with fractional derivative, Appl Math Model. 33 (2009) 4184–4191. http://dx.doi.org/10.1016/j.apm.2009.03.002 [22] S. Rogosin, F. Mainardi, George William Scott Blair – the pioneer of fractional calculus in rheology, Commun Appl Ind Math. 6 1 (2014) e481. http://dx.doi.org/10.1685/journal.caim.481 [23] D.Y. Song, T.Q. Jiang, Study on the constitutive equa- tion with fractional derivative for the viscoelastic fluids – Modified Jeffreys model and its application, Rheol Acta 37 (1998) 512–517. http://dx.doi.org/10.1007/s003970050138 [24] D. Tripathi, S.K. Pandey, S. Das, Peristaltic flow of viscoelastic fluid with fractional Maxwell model through a channel, Appl Math Comput. 215 (2010) 3645–3654. http://dx.doi.org/10.1016/j.amc.2009.11.002 [25] D. Tripathi, Peristaltic transport of fractional Maxwell fluids in uniform tubes: Applications in endoscopy, Com- Biomath 5 (2016), 1605151, http://dx.doi.org/10.11145/j.biomath.2016.05.161 Page 11 of 12 http://dx.doi.org/10.1016/j.rinp.2015.04.003 http://dx.doi.org/10.1016/j.jmmm.2014.12.087 http://dx.doi.org/10.1122/1.549887 http://dx.doi.org/10.2478/s13540-014-0209-x http://dx.doi.org/10.1023/B:NUMA.0000027736.85078.be http://dx.doi.org/10.1023/B:NUMA.0000027736.85078.be http://dx.doi.org/10.1023/A:1016592219341 http://dx.doi.org/10.1016/j.jmmm.2015.05.071 http://dx.doi.org/10.1088/0253-6102/63/3/353 http://dx.doi.org/10.2478/s13540-014-0163-7 http://dx.doi.org/10.2298/TSCI110401077H http://dx.doi.org/10.1146/annurev.fl.03.010171.000305 http://dx.doi.org/10.1515/zna-2014-0258 http://dx.doi.org/10.1016/j.camwa.2010.11.021 http://dx.doi.org/10.2478/s13540-012-0048-6 http://dx.doi.org/10.2478/s13540-013--0003-1 http://dx.doi.org/10.1016/j.apm.2009.03.002 http://dx.doi.org/10.1685/journal.caim.481 http://dx.doi.org/10.1007/s003970050138 http://dx.doi.org/10.1016/j.amc.2009.11.002 http://dx.doi.org/10.11145/j.biomath.2016.05.161 E. Bazhlekova et al., Peristaltic transport of viscoelastic bio-fluids ... put Math Appl. 62 3 (2011) 1116–1126. http://dx.doi.org/10.1016/j.camwa.2011.03.038 [26] D. Tripathi, A mathematical model for the peristaltic flow of chyme movement in small intestine, Math Biosci. 233 (2011) 90–97. http://dx.doi.org/10.1016/j.mbs.2011.06.007 [27] D. Tripathi, O. Anwar Bég, J. L. Curiel-Sosa, Ho- motopy semi-numerical simulation of peristaltic flow of generalised Oldroyd-B fluids with slip effects, Comput Methods Biomech Biomed Engin. 17 4 (2014) 433–442. http://dx.doi.org/10.1080/10255842.2012.688109 [28] K. Vafai, A.A. Khan, S. Sajjad, R. Ellahi, The study of peristaltic motion of third grade fluid under the effects of Hall current and heat transfer, Z Naturforsch. A 70 4 (2015) 281–293 http://dx.doi.org/10.1515/zna-2014-0330 [29] D. Valério, J. Tenreiro Machado, V. Kiryakova, Some pioneers of the applications of fractional calculus, Fract Calc Appl Anal. 17 2 (2014) 552–578. http://dx.doi.org/10.2478/s13540-014-0185-1 [30] P. Yang, K.Q. Zhu, Thermodynamic compatibility and mechanical analogue of the generalized Jeffreys and generalized Oldroyd-B fluids with fractional derivatives, Sci China-Phys Mech Astron. 54 4 (2011) 737–742. http://dx.doi.org/10.1007/s11433-011-4271-7 Biomath 5 (2016), 1605151, http://dx.doi.org/10.11145/j.biomath.2016.05.161 Page 12 of 12 http://dx.doi.org/10.1016/j.camwa.2011.03.038 http://dx.doi.org/10.1016/j.mbs.2011.06.007 http://dx.doi.org/10.1080/10255842.2012.688109 http://dx.doi.org/10.1515/zna-2014-0330 http://dx.doi.org/10.2478/s13540-014-0185-1 http://dx.doi.org/10.1007/s11433-011-4271-7 http://dx.doi.org/10.11145/j.biomath.2016.05.161 Introduction Mathematical model Analytical properties of the pressure gradient function Numerical method Numerical results and discussion Conclusions References