Peristaltic Transport of Viscoelastic Bio-ﬂuids with Fractional Derivative Models

a mathematical study of the peristaltic ﬂow of a viscoelastic ﬂuid through a non-uniform channel is presented under the assumptions of long wavelength and low Reynolds number. The viscoelastic properties of the ﬂuid are modeled by using the fractional Burgers’ constitutive equation, containing four fractional Riemann-Liouville derivatives of orders α, 2 α, β and 2 β , where 0 < α ≤ β < 1. This model covers the fractional Oldroyd-B, Maxwell and second-grade ﬂuid as particular cases. of peristaltic According to the can be obtained as a solution of a diﬀerential equation of fractional order. Our main contri-bution is a detailed study of this equation. Applying Laplace transform, its analytical is represented in terms of Mittag-Leﬄer functions. To the some numerical techniques are compared. The results of several of numerical are and the eﬀects of the diﬀerent parameters on the pressure diﬀerence across one wavelength are


I. INTRODUCTION
Recently, Fractional Calculus has gained considerable popularity mainly due to its numerous applications in diverse fields of science and engineering.Fractional Calculus allows integration and differentiation of arbitrary order, not necessarily 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 classical 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 adequacy preserving linearity and give the possibility for relatively simple description of the complex behavior of non-Newtonian viscoelastic fluids.
The generalized fractional Oldroyd-B constitutive law belongs to the class of linear fractional models of viscoelastic fluids.It is obtained by replacing the first order derivatives in the classical Oldroyd-B model by derivatives of fractional order.The corresponding constitutive equation in the one-dimensional case is given by 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 generalized 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 peristalsis.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 fluids.Examples from physiology where the peristaltic 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 peristaltic 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 transport 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 computations the Adomian decomposition and homotopy analysis methods are used in the aforementioned articles.
In the present work, peristaltic flow of viscoelastic 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 proposed 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 = 0. Since the considered model is non-stationary in nature and contains parameters which are timerelated such as the orders of the fractional time derivatives α and β, the relaxation and retardation times λ 1 , λ 2 , it is natural to study the time evolution 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 = 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.Results 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 analyzed.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 basic definitions and results from Fractional Calculus 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 viscoelastic fluid with fractional Oldroyd-B constitutive model.This equation was originally proposed in [24] for the case λ 2 = 0.The necessary generalizations to the case λ 2 = 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 unsteady motion of an incompressible viscoelastic fluid are the continuity equation: and the general Cauchy momentum equation: where V is the velocity vector, ρ -fluid density, σ -Cauchy stress tensor: 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 Here η > 0 is the dynamic viscosity of the fluid, λ 1 and λ 2 are parameters related to relaxation and retardation times, respectively, satisfying (see [23]) α and β are fractional parameters, A 1 is the first Rivlin-Ericksen tensor given by and D γ Dt γ denotes the upper convected time derivative where D γ t is the Riemann-Liouville fractional derivative, see (51) in the Appendix for the definition.
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 dimensionless 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): On the other hand, the constitutive equation (4) gives: Applying the operator (1 + λ α 1 D α t ) to both sides of (8) the following identity is deduced: Differentiating with respect to y both sides of equation ( 10) one gets From ( 11) and ( 12) one deduces the following equation for the pressure gradient ∂p ∂x Boundary conditions for the velocity are given by ∂u ∂y y=0 = 0, u| y=h = 0.
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 successively Denote by Q the volumetric flow rate Q = h 0 u dy.Then after one more integration (16) implies Following [24] it is assumed that the wall of the channel undergoes contraction and relaxation given by where φ is the amplitude of the wave.The transformations between the wave and the laboratory frames (an established procedure in peristaltic fluid dynamics, see [14]) are given in dimensionless form by Let Q denotes the averaged volumetric flow rate Using the following relation from [14] Eq. ( 17) gives where 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 (24) Here we have used the identity obtained by applying the definition (51) of the Riemann-Liouville fractional derivative and identity (50).Eq. ( 24) implies that the pressure gradient can be expressed as follows where the function y(t) is a solution of the equation Therefore, the time evolution of the pressure gradient 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 = 0. Let us rewrite (27) in the form of the following fractional order equation where Suppose the physically reasonable initial condition Therefore, from (26), y(0) < ∞ and thus lim t→0 where J 1−α t is a Riemann-Liouville fractional integral, see (47) in the Appendix.According to (63) the solution of equation ( 28) is given by 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 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 Inserting (55) into (33) one obtains the following series expansion .
In the simplest case λ 2 = 0, based on representation (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 difference: in this case y(t) does not vanish at t = 0, more precisely, (35) implies 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 = 0: and for λ 2 = 0: Therefore, for λ 2 = 0 as well as for λ 2 = 0 and α > β the function y(t) vanishes for t → 0. If λ 2 = 0 and α = β then the initial value y(0) is as in (37).However, if λ 2 = 0 and α < β, then the asymptotic expansion (38) implies that the function 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.
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 numerical 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 homotopy analysis method (HAM).These two methods 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 p 0 = 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 where .
(43) Equation ( 42) is used here for the numerical computation of the function y(t), applying the socalled fractional Adams method, originally proposed 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 {t j = jh, j = 0, 1, ..., N } with some integer N and h = T /N .Denote by y j the approximation for y(t j ).For the sake of brevity the notation λ = −λ −α 1 is used.The predictor y P k+1 is determined by the formula where The corrector formula is given by and A j,k+1 are defined by Using this numerical algorithm, the function y(t) is computed for several values of the parameters.The performed numerical experiments indicate 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) represent 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.
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 [26] and [27], where the case α ≤ β is considered, this issue has not been addressed.
On Fig. 4 the case α > β is illustrated.Comparing 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    λ 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 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 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 (thermodynamically 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 effects, flows in porous media, etc.

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 integraltransform methods, partial differential equations, special and generalized functions, numerical analysis".

APPENDIX
Here we summarize some facts from Fractional Calculus, for details see [12], [16].
The operators of fractional integration satisfy the semi-group property: or, equivalently, The Riemann-Liouville fractional derivative D α t of order α ∈ (0, 1] is defined by D 1 t = d/dt and The Riemann-Liouville fractional derivatives and integrals are related via the identities: which implies the following identity for the Riemann-Liouville fractional derivative of order α ∈ (0, 1): 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 ( d dt E α,1 (−t) → ∞ for t → 0, see (60)), and this decay is faster for smaller α.