Multi-layer flows of immiscible fractional Maxwell fluids in a cylindrical domain.

Multi-layer flows of immiscible fractional Maxwell fluids in a cylindrical domain.

Journal Pre-proof

Multi-layer flows of immiscible fractional Maxwell fluids in a cylindrical domain. Abdul Rauf, Yasir Mahsud, Imran Siddique PII: DOI: Reference:

S0577-9073(19)30921-9 https://doi.org/10.1016/j.cjph.2019.09.015 CJPH 944

To appear in:

Chinese Journal of Physics

Received date: Revised date: Accepted date:

22 May 2019 11 September 2019 20 September 2019

Please cite this article as: Abdul Rauf, Yasir Mahsud, Imran Siddique, Multi-layer flows of immiscible fractional Maxwell fluids in a cylindrical domain., Chinese Journal of Physics (2019), doi: https://doi.org/10.1016/j.cjph.2019.09.015

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 The Physical Society of the Republic of China (Taiwan). Published by Elsevier B.V. All rights reserved.

Highlights • Stratified flows of Maxwell fluids in cylindrical domain with damped shear stresses are studied. • The memory effects have a significant influence on the motion of the fluids • It is found that for increasing values of the fractional parameter the fluid velocity is decreasing

1

Multi-layer flows of immiscible fractional Maxwell fluids in a cylindrical domain.

Abdul Rauf 1 ∗ , Yasir Mahsud 2 , Imran Siddique 3 1 Dept. of Computer Science and Engineerings, Air University Multan, Pakistan 2 Abdus Salam School of Mathematical Sciences, GC University Lahore, Pakistan 3 Dept. of Mathematics, University of Management & Technology, Lahore, Pakistan

Abstract Laminar unsteady multilayer axial flows of fractional immiscible Maxwell fluids in a circular cylinder are investigated. The flow of fluids is generated by a time-dependent pressure gradient in the axial direction and by the translational motion of a cylinder along his axis. The considered mathematical model is based on the fractional constitutive equation of Maxwell fluids with Caputo time-fractional derivatives. Analytical solutions for the fractional differential equations of the velocity fields with boundary and interfaces conditions have been determined by using the Laplace transform coupled with the Hankel transform of order zero and the Weber transform of order zero. The influence of the memory effects on the motion of the fluid has been investigated for the particular case of three fractional Maxwell fluids. It is found that for increasing values of the fractional parameter the fluid velocity is decreasing. The memory effects have a stronger influence on the velocity of the second layer. Keywords: Maxwell fluid; Fractional constitutive equation; Multi-layer flow; Integral transforms; Velocity field.

1

Introduction

Flows of layered immiscible fluids in a pipe/channel are widely observed in nature. These fluid flows are common in many physiological and industrial applications, such as in the gas-assisted injected moulding, food processing, oil recovery, flows of crude oil in pipelines, equipment cleaning, biomedical contexts (biofilms and mucus flow in human cells), etc [1–4]. Good motivation for the study of multi-layered fluid in a cylindrical domain comes from the common multi-layered flow processes involved in the drilling and the construction of gas and oil wells. Moreover, the drilling and hydraulic fracturing, CO2 sequestration, geothermal and the domestic wells for water distribution are cemented with similar techniques as used in the gas and oil industry. All through the essential establishing process, a progression of ∗

Address correspondence to Abdul Rauf.

E-mail: abdul.rauf @aumc.edu.pk

2

fluids is siphoned down the casing, to remove the boring mud and additionally other in-situ fluids. Failure in such processes is highly costly and catastrophic. It is conventional often use a water-based mud miscible with cement slurry for drilling. But for the nonconventional rapid outputs, the use of water-based mud causes several problems with differential sticking and swelling shales. However, this problem can be resolved if we use an immiscible oil-based mud instead of water-based mud [5]. A wide range of study has been performed to analyze the stability of immiscible fluids [6–9]. Linear stability of flows of viscoelastic fluids with interfaces has been analyzed by Bogaerds et al. [10] by using the time-dependent finite element algorithm. One-dimensional flows of sandwiched viscoelastic fluids which are defined by the Johnson-Segalman and ordinary/modified Phan-Thien Tanner constitutive laws have been investigated by Herve Le Meur [11]. In this paper, the existence and uniqueness for the multilayer Couette/ Poiseuille flows in pipes/channels have been investigated. It is found that the viscosity ratios and the parameter of interpolated Oldroyd derivative play an important role in the uniqueness of the solution. A large volume of analytical, experimental and computational research is available in the literature on the study of the displacement flow of the miscible flow [12–14]. However, the study on the immiscible fluid systems is limited due to the complexity in the fluids flow caused by the fluids interfacial tension. Another important constraint to our study on immiscible systems is the negligible interfacial tension between the two fluids. This reduces the immiscible fluid system to miscible fluid system with infinite Peclet number, or equivalently, zero molecular diffusion [15]. Hasnain et al., [16] studied the two-layer immiscible iso-viscous displacement flow of Newtonian fluids in a cylindrical domain inclined at an angle. It is observed that when the exchange flow is not present in a capillary, the blockage may occur hindering the Rayleigh-Taylor instabilities. In [17], Klubertanz et al. have made a comparative study for the multiphase flow of the immiscible and the miscible fluids in the porous medium. Even if the microscopic basis for both types of the immiscible and the miscible fluids may be rather similar, the macroscopic formulation is quite different as expected. Barrios-Pia et al. [18] studied the effects of vegetation on the multilayered fluid flow pattern. In an open flow channel, emergent and submerged vegetation are studied and compared. It is observed that when the flow is covered by the vegetation, a poor agreement between the experimental measurements and predictions is observed with both emergent and submerged vegetation. Other interesting results for the study of the flow pattern of immiscible fluids can be found in [19, 20]. It is worth point out that many complex systems in chemistry, mechanics or biology are modelled by differential and integral operators of integer order. In recent years, very precise experiments have shown the existence of anomalous processes which cannot be satisfactorily described by differential equations of integer order; these complex systems being described more accurate by the fractional differential equations. There are many fractional differential operators studied by the fractional calculus. Riemann-Liouville fractional integral/derivative [21], Caputo fractional derivative [22], Caputo-Fabrizio fractional derivative [23], AtagnanaBaleanu fractional derivative [24], Yang-Srivastava-Machado fractional derivative [25], YangGao-Machado-Baleanu fractional derivative [26] and the Yang-Abdel-Aty-Cattani fractional derivative [27] are only some from the fractional differential operators used in viscoelasticity, heat and mass transport processes. A Caputo type fractional derivative of variable order has been developed in [28], which has been observed efficient in modelling the concentrations in the complex transport process. A recent study on the fractional calculus operators with 3

local, non-local, singular and non-singular kernels can be seen in [29–33] and the references therein. Hristov [34] has applied the Riemann-Liouville fractional derivative to study the transient space-fractional diffusion with power-law superdiffusivity. Mahsud et al. [35] have investigated boundary layer flows of an upper-convected Maxwell fluid modelled by the time-fractional Caputo-Fabrizio derivatives. They found that fractional fluids flow slower than ordinary fluids. Wang et al. [36] have investigated the unsteady electro-osmotic flows of generalized Maxwell fluids in cylindrical domains and made a comparison between ordinary and generalized Maxwell fluids. Theoretical results and practical applications of the fractional calculus can be found in the references [37–39]. In this paper, we studied laminar unsteady multilayer axial flows of fractional immiscible Maxwell fluids in a moving circular cylinder. The flow of fluids is due to a time-dependent pressure gradient in the axial direction and by the translational motion of the cylinder along his axis. The mathematical model is based on the generalized fractional constitutive equation of Maxwell fluids with Caputo time-fractional derivatives. Analytical solutions for the fractional differential equations of the velocity fields with boundary and interfaces conditions have been determined by using the Laplace transform coupled with the Hankel transform of order zero and the Weber transform of order zero. To carry out numerical calculations, the positive roots of the transcendental equations with Bessel functions J0 (h1 x) = 0 and B0 (hk , x) = J0 (hk x)Y0 (hk−1 x) − J0 (hk−1 x)Y0 (hk x) = 0, k = 2, 3 are determined by using the subroutine “root (f (x), x, a, b)” of the Mathcad 15 package. Obtained numerical results for the velocity field are presented in graphical illustrations for the particular case of three fractional Maxwell fluids. The influence of the memory effects on the fluids motion has been investigated. It is found that for increasing values of the fractional parameter the fluid velocity is decreasing. The memory effects have stronger influence on the velocity of the second layer.

2

Mathematical formulation of the problem

The flow domain is circular cylinder ˜ z˜) ∈ [0, R] × [0, 2π] × (−∞, ∞)}, ˜ = {(˜ D r, θ, with axis of the flow along z˜-axis. Initially at time t˜ = 0, the cylinder and the enclosed fluids are at rest. After this moment, the cylinder is moving along the z˜-axis with the velocity u˜m (R, t˜) = U0 f1 (t˜) (Fig.1). The function f1 (t˜) is a differentiable function of exponential order at infinity and f1 (0) = 0. In this paper the velocity field is considered to be of the form ∼ →
(1) (2) r, t˜), 0, 0 . Here we fix the sets Jm := {1, 2, · · · , m}, Jm := {2, 3, · · · , m}, m ∈ N. V = u˜(˜ (2) We consider here m-layers Maxwell fluids flow in domains r˜ ∈ [0, R1 ], r˜ ∈ [Rk−1 , Rk ], k ∈ Jm inside the cylinder of radius R := Rm . In the domain r˜ ∈ [0, R1 ], flows a Maxwell fluid with ˜ 1 = µ1 , G1 the elastic modulus, velocity u˜1 (˜ r, t˜) the density ρ1 , viscosity µ1 , relaxation time λ G1 (2) and the shear stress τ˜1 (˜ r, t˜). In the annulur domains r˜ ∈ [Rk−1 , Rk ] : Rk−1 < Rk k ∈ Jm ˜ k = µk , Gk the flows a Maxwell fluid with the density ρk , viscosity µk , relaxation time λ Gk elastic modulus, velocity u˜k (˜ r, t˜) and the shear stress τ˜k (˜ r, t˜). The fluids are assumed to be incompressible and immiscible and the flow is unsteady, onedimensional and fully developed. The flows of the fluids are generated by the movement of the 4

Figure 1: Geometry of the problem boundary of the cylinder and by the time-dependent pressure gradient in the flow direction. Based on the assumptions about velocity fields, the continuity equation is identically satisfied (1) by all velocities u˜k (˜ r, t˜), k ∈ Jm . The governing equations of motion are [40] • the system of linear momentum equations ρk

∂ u˜k ∂ p˜ 1 ∂ (1) =− + (˜ rτ˜k ), k ∈ Jm , ˜ ∂˜ r r˜ ∂˜ r ∂t

(1)

• the system of constitutive equations ˜ k ∂ τ˜k = µk ∂ u˜k , k ∈ J (1) . τ˜k + λ m ∂˜ r ∂ t˜

(2)

Along with the above equations, we consider the following initial, boundary and the interface fluid-fluid conditions: • initial conditions, (1) , u˜k (˜ r, 0) = 0, τ˜k (˜ r, 0) = 0, k ∈ Jm

(3)

∂ u˜1 (˜ r, t˜) |r˜=0 = 0, u˜m (R, t˜) = U0 f1 (t˜), ∂˜ r

(4)

• boundary condition,

5

• the interfaces conditions, (1) ˜ j , t˜) = τ˜j+1 (Rj , t˜), j ∈ Jm−1 u˜j (Rj , t˜) = u˜j+1 (Rj , t˜), τ˜j (R .

(5)

Introducing the non-dimensional variables ˜ τk r˜ , t = UR0 t , uk = Uu˜k0 , τk = µR˜ , p = ρ1p˜U 2 , R k U0 0  ˜ Rt k , f (t) = f λk = λkRU0 , ak = ρρk1 , bk = ρ1µRU 1 U0 , hj 0 (1) (1) k ∈ Jm ; j ∈ Jm−1 ,

r=

=

Rj , R

into Eqs. (1)-(5), the mathematical model for dimensionless governing equations takes the form, ak

∂p 1 ∂ ∂uk (1) , = − + bk (rτk ), k ∈ Jm ∂t ∂r r ∂r τk + λk

(6)

∂τk ∂uk (1) = , k ∈ Jm , ∂t ∂r

(7)

along with • the non-dimensional initial conditions (1) uk (r, 0) = 0, τk (r, 0) = 0, k ∈ Jm ,

(8)

• the non-dimensional boundary condition ∂u1 (r, t) |r=0 = 0, um (1, t) = f (t), ∂r

(9)

• the non-dimensional interface conditions (1)

uj (hj , t) = uj+1 (hj , t), τj (hj , t) = τj+1 (hj , t), j ∈ Jm−1 .

2.1

(10)

Mathematical model with fractional constitutive equations

In the following , we consider the generalized constitutive equations [35, 39], τk + λk Dαt k τk =

∂uk (1) , αk ∈ (0, 1], k ∈ Jm , ∂r

(11)

where Dγt F (y, τ )

1 = Γ(1 − γ)

Z

t

0

(t − τ )−γ

∂F (y, τ ) dτ, 0 ≤ γ < 1 ∂τ

is the Caputo derivative [23]. For particular case γ = 1, D1t F (y, t) = that the pressure gradient is a known function of t, therefore −

∂F (y,t) . ∂t

(12) We suppose

∂p = P (t), ∂r

where P (t) is a local integrable function on [0, ∞) of the exponential order at infinity. 6

(13)

3

Solution of the problem

3.1

The case k = 1

For the case k = 1, the fractional model of the governing equations describing the fluid flow in the cylindrical domain r ∈ [0, h1 ] can be written as a1

1 ∂ ∂u1 = P (t) + b1 (rτ1 ), ∂t r ∂r

τ1 + λ1 Dαt 1 τ1 =

∂u1 , α1 ∈ (0, 1], ∂r

(14)

(15)

along with • the non-dimensional initial conditions u1 (r, 0) = 0, τ1 (r, 0) = 0,

(16)

• the non-dimensional boundary condition ∂u1 (r, t) |r=0 = 0. ∂r

(17)

In order to determine analytical solutions of the Eqs. (14), (15) with conditions (16), (17) we use the Laplace transform with respect to the variable t, coupled with the finite Hankel transform with respect to the variable r, namely, Z h1 u˜1 (r1n , t) = u1 (r, t)rJ0 (rr1n )dr, n = 1, 2, · · · , (18) 0

with the inverse Hankel transform ∞ 2 X J0 (rr1n ) u˜1 (r1n , t), u1 (r, t) = 2 h1 n=1 J12 (h1 r1n )

(19)

where r1n , n = 1, 2, · · · , are the positive roots of the transcendental equation J0 (h1 x) = 0. Applying the Laplace transform to Eqs. (14), (15), (17) and using the initial conditions (16), we obtain the problem in transformed domain 1 ∂ a1 q¯ u1 (r, q) = P¯ (q) + b1 (r¯ τ1 (r, q)), r ∂r

τ¯1 (r, q) =

∂ u¯1 (r, q) 1 , α 1 + λ1 q 1 ∂r

∂ u¯1 (r, q) |r=0 = 0. ∂r 7

(20)

(21)

(22)

Substituting Eq. (21) in Eq. (20) gives a1 q¯ u1 (r, q) = P¯ (q) +

1 ∂ ∂ u¯1 (r, q) b1 (r ), 1 + λ1 q α1 r ∂r ∂r

Eq. (23) with the application of finite Hankel transform (18) takes the form   J (h r ) b ' ' 1 1 1n 1 2 a1 q u1 = h1 P¯ (q) + h1 u¯1 (h1 , q)r1n J1 (h1 r1n ) − r1n u1 r1n 1 + λ1 q α1

(23)

(24)

'

Rearranging Eq. (24) we get the transformed form of the velocity u1 '

u1 (r1n , q) =

(1 + λ1 q α1 )h1 P¯ (q) J1 (h1 r1n ) b1 h1 u¯1 (h1 , q) r1n J1 (h1 r1n ) (25) + α 2 1 a1 q(1 + λ1 q ) + b1 r1n r1n a1 q(1 + λ1 q α1 ) + b1 r1n 2

In order to apply the inverse Hankel transform (19) to Eq. (25) we rewrite Eq. (25) in the following equivalent form '

u1 (r1n , q) =

u¯1 (h1 , q) 2 2 (h1 r1n − 4)J1 (h1 r1n ) 3 h1 r1n

(26)

+P¯ (q)A1n (q)J1 (h1 r1n ) + B1n (q)J1 (h1 r1n )¯ u1 (h1 , q), where A1n (q) =

B1n (q) =

(1 + λ1 q α1 )h1 , r1n (a1 q(1 + λ1 q α1 ) + b1 r1n 2 )

2 2 4b1 r1n − (h21 r1n − 4)(1 + λ1 q α1 )a1 q . 3 (a1 q(1 + λ1 q α1 ) + b1 r1n 2 ) h1 r1n

Applying inverse Hankel transform (19) to Eq. (26) we get ∞ ∞ ¯1n (q) u¯1 (h1 , q) 2 2P¯ (q) X A¯1n (q) 2¯ u1 (h1 , q) X B u¯1 (r, q) = r + J (rr ) + J0 (rr1n ) (27) 0 1n h21 h21 n=1 J1 (h1 r1n ) h21 J (h r ) n=1 1 1 1n

To obtain the inverse Laplace transforms of the Eq. (27) we consider the following (1) auxiliary functions for k ∈ Jm , ¯ k0 (n, q) = L

=

1 2 ak q(1 + λk q αk ) + bk rkn

q −1 · ak λk (q αk + λ−1 k ) 1+

1 2 q −1 bk rkn ak λk (q αk +λ−1 k )

 2 i   ∞ 1 X q −i−1 i bk rkn , = (−1) i+1 ak λk i=0 ak λk (q αk + λ−1 k ) 8

(28)

¯ k4 (n, q) = L

1 ¯ k0 (n, q) L 1 + λk q αk

 2 i   ∞ q −i−1 1 X i bk rkn (−1) , = i+2 ak λ2k i=0 ak λk (q αk + λ−1 k )

(29)

¯ 13 (n, q) = q L ¯ 10 (n, q) L  2 i   ∞ 1 X q −i i b1 r1n = (−1) , i+1 a1 λ1 i=0 a1 λ 1 (q α1 + λ−1 1 ) ¯ k1 (n, q) = L

(30)

(1 + λk q αk ) ¯ k0 (n, q) = λk (q αk + λk −1 )L 2 ak q(1 + λk q αk ) + bk rkn

 2 i   ∞ 1X q −i−1 i bk rkn = , (−1) i ak i=0 ak λk (q αk + λ−1 k )

(31)

and ¯ k2 (n, q) := q L ¯ k1 (n, q) L  2 i   ∞ 1 1X q −i i bk rkn = + (−1) , i ak ak i=1 ak λk (q αk + λ−1 k ) Since the generalized G-Lorenzo-Hartley function is defined by [41]  X  ∞ qb Γ(k + c)dk t(k+c)a−b−1 −1 Ga,b,c (t, d) = L = , (q a − d)c k!Γ((k + c)a − b)Γ(c) k=0 Re(q) > 0, Re(ac − b) > 0, |

(32)

(33)

d | < 1, qa

and for αi , βi > 0, the Mittag-Leffler function are defined as [42]  αi −βi  −1 q L = tβi −1 Eαi ,βi (dtαi ), q αi − d

(34)

we have:   2 i (bk rkn ) −1 Lk0 (n, t) = (−1) Gαk ,−i−1,i+1 (t, −λk ) , i+1 (a λ ) k k i=0

(35)

  ∞ 2 i 1X −1 i (bk rkn ) L14 (n, t) = (−1) Gαk ,−i−1,i+2 (t, −λk ) , λk i=0 (ak λk )i+1

(36)

∞ X

i

9

  2 i ) (b1 r1n −1 L13 (n, t) = (−1) Gα1 ,−i,i+1 (t, −λ1 ) , i+1 (a λ ) 1 1 i=0

(37)

  ∞ 2 i 1X −1 i (bk rkn ) Lk1 (n, t) = Gαk ,−i−1,i (t, −λk ) , (−1) ak i=0 (ak λk )i

(38)

  ∞ 2 i 1X 1 −1 i (bk rkn ) (−1) Gαk ,−i,i (t, −λk ) , Lk2 (n, t) = δ(t) + ak ak i=1 (ak λk )i

(39)

∞ X

i

and

Applying inverse Laplace transform to Eq. (27) we get the velocity distribution profile u1 (r, t)  ∞ 2 X J0 (rr1n ) h1 u1 (h1 , t) 2 r + 2 P (t) ∗ L11 (n, t) (40) u1 (r, t) = h21 h1 n=1 J1 (h1 r1n ) r1n  2 − 4) 4b1 a1 (h21 r1n + u1 (h1 , t) ∗ L10 (n, t) − u1 (h1 , t) ∗ L12 (n, t) , 3 h1 r1n h1 r1n Rt where h1 (t) ∗ h2 (t) = 0 h1 (t − τ )h2 (τ )dτ is the convolution product of the functions h1 (t), h2 (t). With the help of Eqs. (21) and (27), the transformed shear stress τ¯1 (r, q) takes the form τ¯1 (r, q) = T¯21 (r, q)¯ u1 (h1 , q) + T¯31 (r, q),

(41)

where T¯21 (r, q) =

2 h21 (1+λ1 q α1 )

T¯31 (r, q) =



r−

2P¯ (q) − h2 (1+λ α 1q 1 ) 1

∞ P

n=1

∞ P

n=1



¯1n (q)r1n B J (rr1n ) J1 (h1 r1n ) 1

, (42)

¯1n (q)r1n A J (rr1n ). J1 (h1 r1n ) 1

Applying inverse Laplace transform to Eq. (41) we get the following expression of shear stress τ1 (r, t) τ1 (r, q) = T21 (r, t) ∗ u1 (h1 , t) + T31 (r, t),

(43)

where T21 (r, t) = 1 − 8b h3 1

P∞

2r α1 −1 t Eα1 ,α1 (− λ11 tα1 ) h21 λ1

J1 (rr1n ) n=1 J1 (h1 r1n ) L14 (n, t) ¯

T31 (r, t) = − 2Ph(q) 1

+

2a1 h31

P∞

P∞

2 −4) J1 (rr1n ) (h21 r1n L13 (n, t), 2 n=1 J1 (h1 r1n ) r1n

J1 (rr1n ) n=1 J1 (h1 r1n ) P (t)

10

∗ L10 (n, t).

(44)

3.2

(2)

Solution for the annular domain r ∈ [hk−1 , hk ], k ∈ Jm : (2)

For the case k ∈ Jm , the fractional model of the governing equations describing the fluid flow in the cylindrical domain r ∈ [hk−1 , hk ] can be written as ak

1 ∂ ∂uk (2) , = P (t) + bk (rτk ), k ∈ Jm ∂t r ∂r

τk + λk Dαt k τk =

(45)

∂uk (2) , , αk ∈ (0, 1], k ∈ Jm ∂r

(46)

along with • the non-dimensional initial conditions (2) uk (r, 0) = 0, τk (r, 0) = 0, k ∈ Jm ,

(47)

• the non-dimensional boundary condition um (1, t) = f (t),

(48)

• the non-dimensional interface conditions (1)

uj (hj , t) = uj+1 (hj , t), τj (hj , t) = τj+1 (hj , t), j ∈ Jm−1 .

(49)

In order to determine analytical solutions of the Eqs. (45), (46) with conditions (47), (48) and (49) we use the Laplace transform with respect to the variable t, coupled with the finite Weber transform [43] with respect to the variable r namely, Z hk ∗ uk (r, t)rB0 (r, rkn )dr, n = 1, 2, · · · , (50) uk (rkn , t) = hk−1

with the inverse transform uk (r, t) =

∞ 2 π 2 X rkn J02 (hk rkn )B0 (r, rkn ) ∗ uk (rkn , t), 2 n=1 J02 (hk−1 rkn ) − J02 (hk rkn )

(51)

where B0 (r, rkn ) = J0 (rrkn )Y0 (hk−1 rkn ) − J0 (hk−1 rkn )Y0 (rrkn ) (2)

and rkn , n = 1, 2, · · · , k ∈ Jm , being the positive roots of the transcendental equation B0 (hk , x) = 0. Applying the Laplace transform to Eqs. (45), (46), (48)and (49) and using the initial conditions (47), we obtain the problem in the transformed domain ak q¯ uk (r, q) = P¯ (q) + bk

τ¯k (r, q) =

1 ∂ (2) (r¯ τk (r, q)), k ∈ Jm , r ∂r

1 ∂ u¯k (r, q) (2) , k ∈ Jm , 1 + λk q αk ∂r 11

(52)

(53)

u¯m (1, q) = f¯(q),

(54) (1)

u¯j (hj , q) = u¯j+1 (hj , q), τ¯j (hj , q) = τ¯j+1 (hj , q), j ∈ Jm−1 .

(55)

Substituting Eq. (53) into the Eq. (52) gives ak q¯ uk (r, q) = P¯ (q) +

1 ∂ ∂ u¯k (r, q) bk (2) (r ), k ∈ Jm , α k 1 + λk q r ∂r ∂r

(56)

Eq. (56) with the application of finite Weber transform (50) takes the form ak q

∗

− uk (rkn , q)

∗ bk 2P¯ (q) J0 (hk−1 rkn ) 2P¯ (q) 2 − − − r uk (rkn , q)+ = 2 2 πrkn J0 (hk rkn ) πrkn 1 + λk q αk kn

(57)

bk 2 J0 (hk−1 rkn ) bk 2 (2) , u¯k (hk , q) − u¯k (hk−1 , q), k ∈ Jm α α 1 + λk q k π J0 (hk rkn ) 1 + λk q k π 2 where with the help of Wronskian’s Bessel relation J0 (z)Y1 (z) − J1 (z)Y0 (z) = − πz , and other basic properties of Bessel functions we have, ∗

1= Z

hk

hk−1

2 2 J0 (hk−1 rkn ) − 2 , 2 πrkn J0 (hk rkn ) πrkn

(58)

∗ ∂ ∂ u¯k (r, q) 2 J0 (hk−1 rkn ) 2 2 − (r ) B0 (r, rkn )dr = −rkn uk (rkn , q) + u¯k (hk , q) − u¯k (hk−1 , q), ∂r ∂r π J0 (hk rkn ) π

n = 1, 2, · · · .

(59)

Rearranging Eq. (57) we get the transformed form of the velocity ∗

−

uk (rkn , q) =

+

∗

− uk (rkn , q)

2P¯ (q) 1 + λk q αk 1 + λk q αk J0 (hk−1 rkn ) 2P¯ (q) − (60) 2 2 2 2 πrkn ak q(1 + λk q αk ) + bk rkn J0 (hk rkn ) πrkn ak q(1 + λk q αk ) + bk rkn

bk 2 J0 (hk−1 rkn ) bk 2 (2) u¯k (hk , q) − u¯k (hk−1 , q), k ∈ Jm . 2 2 α α k k ak q(1 + λk q ) + bk rkn π J0 (hk rkn ) ak q(1 + λk q ) + bk rkn π

To apply the inverse Weber transform, we rewrite the Eq. (60) into the following equivalent form ∗

− uk (rkn , q)

∗

∗

= F1 (rkn )¯ uk (hk−1 , q) + F2 (rkn )¯ uk (hk , q)

(2) +Skn (q)¯ uk (hk−1 , q) + Qkn (q)¯ uk (hk , q) + Rkn (q)P (q), k ∈ Jm ,

where   2 2rkn (h2k − h2k−1 ) + 8 1 8J0 (hk−1 rkn ) F1 (rkn ) = 2 , − 4 4 hk − h2k−1 πrkn J0 (hk rkn ) πrkn ∗

12

(61)

 2 2  (2rkn (hk − h2k−1 ) − 8)J0 (hk−1 rkn ) 8 1 + 4 , F2 (rkn ) = 2 4 J0 (hk rkn ) hk − h2k−1 πrkn πrkn ∗

∗

k Skn (q) = − π[ak q(1+λ2b α 2 ] − F1 (rkn ) k q k )+bk r kn

∗

¯ k0 (n, q) − F1 (rkn ), = − 2bπk L Qkn (q) = = Rkn (q) = =

bk 2 J0 (hk−1 rkn ) 2 π J (h r ak q(1+λk q αk )+bk rkn 0 k kn ) 2bk J0 (hk−1 rkn ) ¯ Lk0 (n, q) π J0 (hk rkn )

∗

− F2 (rkn )

∗

− F2 (rkn ),

J0 (hk−1 rkn )−J0 (hk rkn ) 1+λk q αk 2 2 a q(1+λ q αk )+b r 2 J0 (hk rkn ) πrkn k kn k k 2 J0 (hk−1 rkn )−J0 (hk rkn ) ¯ Lk1 (n, q). 2 J0 (hk rkn ) πrkn

Applying the inverse Weber transform (51) to Eq. (61) we get −

uk (r, q) =

r2 − h2k−1 h2k − r2 u ¯ (h , q) + u¯k (hk , q) k k−1 h2k − h2k−1 h2k − h2k−1

(62)

∞ 2 X π2 rkn J02 (hk rkn )Skn (q) B0 (r, rkn ) + u¯k (hk−1 , q) 2 J 2 (hk−1 rkn ) − J02 (hk rkn ) n=1 0 ∞ 2 X π2 rkn J02 (hk rkn )Qkn (q) + u¯k (hk , q) B0 (r, rkn ) 2 J 2 (hk−1 rkn ) − J02 (hk rkn ) n=1 0

+πP (q)

∞ X n=1

¯ k1 (n, q) J0 (hk rkn )L (2) B0 (r, rkn ), k ∈ Jm . J0 (hk−1 rkn ) + J0 (hk rkn )

Notice that dB0 (r, rkn ) = −rkn (J1 (rrkn )Y0 (hk−1 rkn ) − J0 (hk−1 rkn )Y1 (rrrrkn )) := −rkn B10 (r, rkn ). (63) dr Using Eq. (53) in Eq. (62) τ¯k (r, q) = −

+

h2k

2r u¯k (hk−1 , q) 2r u¯k (hk , q) + 2 2 2 α − hk−1 1 + λk q k hk − hk−1 1 + λk q αk ∞

3 π 2 u¯k (hk−1 , q) X rkn J02 (hk rkn )Skn (q) B10 (r, rkn ) 2 1 + λk q αk n=1 J02 (hk−1 rkn ) − J02 (hk rkn )

13

(64)

−

∞ 3 J02 (hk rkn )Qkn (q) π 2 u¯k (hk , q) X rkn B10 (r, rkn ) 2 1 + λk q αk n=1 J02 (hk−1 rkn ) − J02 (hk rkn )

−π P¯ (q)

∞ X ¯ k0 (n, q) rkn J0 (hk rkn )L (2) B10 (r, rkn ), k ∈ Jm . J (h r ) + J0 (hk rkn ) n=1 0 k−1 kn

With the notations, T¯1k (r, q) = − 1+λ1k qαk T¯2k (r, q) =

1 1+λk q αk





2r h2k −h2k−1

2r h2k −h2k−1

P T¯3k (r, q) = −π P¯ (q) ∞ n=1

+

−

π2 2

π2 2



P∞

3 J 2 (h r rkn k kn )Skn (q) 0 n=1 J02 (hk−1 rkn )−J02 (hk rkn ) B10 (r, rkn )



P∞

3 J 2 (h r rkn k kn )Qkn (q) 0 2 n=1 J0 (hk−1 rkn )−J02 (hk rkn ) B10 (r, rkn )

¯ k0 (n,q) rkn J0 (hk rkn )L B (r, rkn ), J0 (hk−1 rkn )+J0 (hk rkn ) 10

,

,

(65)

(2)

k ∈ Jm ,

Eq. (64) can be written in equivalent form (2) τ¯k (r, q) = T¯1k (r, q)¯ uk (hk−1 , q) + T¯2k (r, q)¯ uk (hk , q) + T¯3k (r, q), k ∈ Jm .

(66)

Using the interface conditions (55) in Eqs. (27) and (66) we obtained the following algebraic system (see Appendix 1) ¯ (q)U¯ (q) = N ¯ (q), M

(67)

where ¯ (q) = (M ¯ kj (q)) M (1) ; k,j∈Jm−1 ¯ kj (q) = T¯1k (hk , q)δk−1,j + (T¯2k (hk , q) − T¯1(k+1) (hk , q))δk−1,j−1 M (1) −T¯2(k+1) (hk , q)δk−1,j−2 , k, j ∈ Jm−1 , ¯ (q) = (N ¯k (q)) (1) ; N k∈Jm−2 (1) ¯ ¯ ¯ Nk (q) = T3(k+1) (hk , q) − T3k (hk , q), k ∈ Jm−2 , ¯m−1 (q) = T¯3m (hm−1 , q) − T¯3(m−1) (hm−1 , q) − T¯2m (hm−1 , q)¯ N um (hm = 1, q),

(68)

U¯ (q) = (¯ uk (hk , q))k∈J (1)

m−1 ;

(1) U¯k (q) = u¯k (hk , q), k ∈ Jm−1 ,

where u¯n (hn = 1, q) = f¯(q) is known. Finally, we obtain ¯ −1 (q)N ¯ (q). U¯ (q) = M

(69)

(1)

Now, uk (hk , q), k ∈ Jm are known functions, therefore the velocities u¯1 (r, s), · · · , u¯m (r, s) are known. Applying inverse Laplace transform to Eq. (62) we get the velocities distribution (2) profiles uk (r, t), k ∈ Jm , as uk (r, t) =

r2 − h2k−1 h2k − r2 u (h , t) + uk (hk , t) k k−1 h2k − h2k−1 h2k − h2k−1 14

(70)

+

∞ 2 J02 (hk rkn )B0 (r, rkn ) π 2 X rkn uk (hk−1 , t) ∗ Skn (t) 2 n=1 J02 (hk−1 rkn ) − J02 (hk rkn )

∞ 2 π 2 X rkn J02 (hk rkn )B0 (r, rkn ) + uk (hk , t) ∗ Qkn (t) 2 n=1 J02 (hk−1 rkn ) − J02 (hk rkn ) ∞ 2 X rkn J0 (hk rkn )B0 (r, rkn ) (2) +π P (t) ∗ Lk1 (n, t), k ∈ Jm , J (h r ) + J (h r ) 0 k−1 kn 0 k kn n=1

where Skn (t) = − Qkn (t) =

∗ 2bk Lk0 (n, t) − F1 (rkn )δ(t), π

∗ 2bk J0 (hk−1 rkn ) Lk0 (n, t) − F2 (rkn )δ(t). π J0 (hk rkn ) (2)

Applying inverse Laplace transform to Eq. (66) we get the shear stresses τk (r, t), k ∈ Jm , as (2) τk (r, t) = T1k (r, t) ∗ uk (hk−1 , t) + T2k (r, t) ∗ uk (hk , t) + T3k (r, t), k ∈ Jm ,

(71)

where T1k (r, t) = 

2bk Lk4 (n, t) π

T2k (r, t) = 

tαk −1 2r Eαk ,αk − h2 −h 2 λk k k−1

+

4

−

tαk λk



∗

+

α −1 F1 (rkn ) t λkk Eαk ,αk

2r tαk −1 − h2 −h Eαk ,αk 2 λk k k−1

2bk J0 (hk−1 rkn ) Lk4 (n, t) π J0 (hk rkn )

T3k (r, t) = −π



−

P∞



−

tαk λk



∗



−

π2 2

− π2 2

α −1 F1 (rkn ) t λkk Eαk ,αk

P∞

3 J 2 (h r rkn k kn )B10 (r,rkn ) 0 n=1 J02 (hk−1 rkn )−J02 (hk rkn ) ×

tαk λk



∗ uk (hk−1 , t),

P∞

3 J 2 (h r rkn k kn )B10 (r,rkn ) 0 n=1 J02 (hk−1 rkn )−J02 (hk rkn ) ×



−

tαk λk



rkn J0 (hk rkn )B10 (r,rkn ) n=1 J0 (hk−1 rkn )+J0 (hk rkn ) (Lk4 (n, t)

(72)

∗ uk (hk−1 , t), ∗ P (t)),

(2)

k ∈ Jm .

Numerical results and discussions

Unsteady multilayer axial flows of fractional immiscible Maxwell fluids in a circular cylinder of radius R have been investigated. The flow of fluids is influenced by a time-dependent pressure gradient in the axial direction and an external magnetic field in the radial direction. The mathematical model is based on the generalized fractional constitutive equation for the 15

Maxwell fluids. The time-fractional Caputo derivative operator is used for the generalization of the constitutive equation. It is important to note that the constitutive equation (11) can be written in the equivalent form τk (r, t) =

Zt 0

(t − σ)αk −1 αk
∂u(r, σ) (t − σ) Eαk ,αk −λ−1 dσ k λk ∂r

(73)

where Eα,β (·) is the two parameters Mittag-Leffler function. It is observed from Eq. (73) that the histories of the velocity derivative with respect to the variable r influence the shear stress of the fluid, therefore the memory effects are considered in the mathematical model. In order to investigate physical aspects of the fluids motion, we have considered the case of three Maxwell fluids in a circular cylinder of radius R = 0.3 [m]. In the inner core r ∈ [0, h1 ], h1 = 0.1 flows a fractional Maxwell fluid with ρ1 = 1000 [kg/m3 ], µ1 = 0.05 [N s/m2 ], G1 = 1.2 [N/m2 ], in the annular region r ∈ [h1 , h2 ], h2 = 0.4 flows a fractional Maxwell fluid characterized by ρ2 = 1300 [kg/m3 ], µ2 = 0.2 [N s/m2 ], G2 = 2.4 [N/m2 ] and, in the external layer r ∈ [h2 , h3 ], h3 = 1 flows a fractional Maxwell fluid with ρ3 = 1500 [kg/m3 ], µ3 = 0.6 [N s/m2 ], G3 = 5.0 [N/m2 ]. To the characteristic velocity we considered the value U0 = 0.1 [m/s], for the time= P (t) = 2 cos(t) and for the boundary velocity u(1, t) = f (t) = pressure gradient −∂p ∂r 5(1 − cos(t)). Figures 2-4 are plotted to investigate the influence of the fractional parameters of Caputo time-fractional derivatives on the fluid velocity. In Fig. 2 the fractional parameter α1 of the fractional Maxwell fluid in the inner core is kept constant while, the fractional parameters α2 and α3 have different values. In Fig. 3 the fractional parameter α2 is constant; respectively in Fig. 4 the fractional parameter α3 is constant. It is observed from these figures that the influence of memory effects on the motion of the fluids in the inner core and the external layer is insignificant. The velocity of the fluid from the second layer is strongly affected by the memory effects and due to the interfaces interactions. Fluid velocity decreases with the fractional parameters, because, the weigh
α −1 αk t in Eq. (73) is an increasing function with respect function Fαk (t) = t λkk Eαk ,αk −λ−1 k to the parameter αk ∈ [0, 0.7] (See Fig. 8), therefore, when the values of the fractional parameter increase the values of the shear stress increases and the velocity decreases. Figures 5-7 have been plotted to analyze the influence of the non-dimensional parameters bk , k = 1, 2, 3 on the fluid motion. In Figs. 5-7 we used the following values of the material coefficients: ρ1 = 900, µ1 = 0.05, G1 = 1.2, ρ2 = 1400, µ1 = 0.2, G2 = 2.4, ρ3 = 2000, µ3 = 0.6, G3 = 5.0 For these figures we have used the characteristic velocities U0 ∈ {0.100, 0.105, 0.110} which correspond to the following parameters: b1 = 1.852 × 10−3 , b2 = 7.407 × 10−3 , b3 = 0.022 to the Fig. 5, b1 = 1.764 × 10−3 , b2 = 7.055 × 10−3 , b3 = 0.021 to the Fig. 6, respectively b1 = 1.684 × 10−3 , b2 = 6.734 × 10−3 , b3 = 0.020 to the Fig. 7. It is observed from these figures that the influence of the memory effects is stronger when the values of the parameters bk , k = 1, 2, 3 are small. This property is more visible on the graphs of the velocity u2 (r, t) in Figs 5-7.

16

Finally, we give a comparison between the fractional model with the Caputo timefractional derivative and the fractional model described by the Yang-Srivastava-Machado time-fractional derivative [25]. By considering now the constitutive equation of the shear stress in the form τk + λk Y SM Dtαk τk =

∂uk (1) , αk ∈ (0, 1], k ∈ Jm , ∂r

(74)

with Y SM

Dtαk τk (r, t)

1 d = 1 − αk dt

Z

t

exp

0



 αk (t − σ) − τk (r, σ)dσ, 1 − αk

the Yang-Srivastava-Machado derivative, we can write that Z t ∂u(r, σ) ∂u(r, t) Gαk (t − σ) = dσ, τk (r, t) = Gαk (t) ∗ ∂r ∂r 0

(75)

(76)

where, the kernel Gαk (t) is given by Gαk (t) = 1 −

2mk exp(−mk t) sinh(nk t), nk

(77)

and αk αk2 − 4λk (1 − αk ) 2 mk = , nk = . 2(1 − αk ) 4(1 − αk )2

(78)

Profiles of kernel Gα (t), α ∈ (0, 1) and t ∈ {0.15, 0.20, 0.30} are plotted in Fig. 9. Comparing Figs. 8 and 9 it can be observed that, in the two models the memory effect is opposite, namely; in the model based on Caputo time-fractional derivative the damping of velocity gradient is bigger for larger values of the fractional parameter, while in the model based on the Yang-Srivastava-Machado time-fractional derivative, the damping is bigger for small values of the fractional parameter.

Appendix 1 Using the interface conditions (55) in Eqs. (27) and (66) we obtained the following algebraic system a ¯1,1 u¯1 (h1 , q) − a ¯1,2 u¯2 (h2 , q) = ¯b1 , (2) a ¯k,k−1 u¯k−1 (hk−1 , q) + a ¯k,k u¯k (hk , q) − a ¯k,k+1 u¯k+1 (hk+1 , q) = ¯bk , k ∈ Jm−2 ,

(79)

a ¯m−1,m−2 u¯m−2 (hm−2 , q) + a ¯m−1,m−1 u¯m−1 (hm−1 , q) = ¯bm−1 , where a ¯1,1 = T¯2,1 (h1 , q) − T¯1,2 (h1 , q), a ¯1,2 = T¯2,2 (h1 , q), ¯ ¯ a ¯k,k−1 = T1,k (hk , q), a ¯k,k = T2,k (hk , q) − T¯1,k+1 (hk , q), (2) a ¯k,k+1 = T¯2,k+1 (hk , q), k ∈ Jm−1 , 17

(80)

and (1) ¯bj = T¯3,j+1 (hj , q) − T¯3,j (hj , q) + δj,m−1 a ¯m−1,m , j ∈ Jm−1 ,

(81)

δj,m−1 is Kronecker delta. The system (79) can be written in the following equivalent form ¯ (q)U¯ (q) = N ¯ (q), M

(82)

where 

     ¯ M (q) =     

a ¯1,1 −¯ a1,2 0 0 ··· 0 a ¯2,1 a ¯2,2 −¯ a2,3 0 ··· 0 0 a ¯3,2 a ¯3,3 −¯ a3,4 ··· 0 .. .. .. .. .. .. . . . . . . 0 0 ··· a ¯m−3,m−4 a ¯m−3,m−3 −¯ am−3,m−2 0 0 ··· a ¯m−2,m−3 a ¯m−2,m−2 −¯ am−2,m−1 0 0 ··· 0 a ¯m−1,m−2 a ¯m−1,m−1 

and

  U¯ (q) =  

u¯m−1 (hm−1 , q)

  ¯ N (q) =  

¯b1 ¯b2 .. . ¯bm−1

     ,    



u¯1 (h1 , q) u¯2 (h2 , q) .. .





  , 



  . 

(83)

(84)

(85)

References [1] Nelson, E. B. & Gulliot, D. 2006 Well Cementing, 2nd edn. Schlumberger Educational Services. [2] Cole, P.A., Asteriadou, K., Robbins, P.T., Owen, E.G., Montague, G.A. and Fryer, P.J., 2010. Comparison of cleaning of toothpaste from surfaces and pilot scale pipework. Food and bioproducts processing, 88(4), pp.392-400. [3] Burfoot, D., Middleton, K.E. and Holah, J.T., 2009. Removal of biofilms and stubborn soil by pressure washing. Trends in Food Science & Technology, 20, pp.S45-S47. [4] Joseph, D.D., Bai, R., Chen, K.P. and Renardy, Y.Y., 1997. Core-annular flows. Annual Review of Fluid Mechanics, 29(1), pp.65-90. [5] Al Khayyat, B., Morris, J.A., Salahudin, S., Ravi, K. and Faria, J., 1999, January. Successes in production-liner cementing in oil-based mud: A case study. In SPE/IADC Middle East Drilling Technology Conference. Society of Petroleum Engineers. [6] Yih, C.S., 1967. Instability due to viscosity stratification. Journal of Fluid Mechanics, 27(2), pp.337352.

18

Fig. 2. Profiles of velocities u1(r, t ) u2 (r, t ) u3 (r, t ) for different values of the fractional parameters  2 ,  3 with  constant and the time t = 0.2. 1

19

Fig. 3. Profiles of velocities u1(r, t ) u2 (r, t ) u3 (r, t ) for different values of the fractional parameters 1,  3 with  constant and the time t = 0.2. 2

20

Fig. 4. Profiles of velocities u1(r, t ) u2 (r, t ) u3 (r, t ) for different values of the fractional parameters 1,  2 with  3 constant and the time t = 0.2.

21

Fig. 5. Profiles of velocities u1(r, t ) u2 (r, t ) u3 (r, t ) for different values of the fractional parameters 1,  2 ,  3 and the time t = 0.25.

22

Fig. 6. Profiles of velocities u1(r, t ) u2 (r, t ) u3 (r, t ) for different values of the fractional parameters 1,  2 ,  3 and the time t = 0.25.

23

Fig. 7. Profiles of velocities u1(r, t ) u2 (r, t ) u3 (r, t ) for different values of the fractional parameters 1,  2 ,  3 and the time t = 0.25.

24

25

[7] Joseph, D.D., Renardy, M. and Renardy, Y., 1984. Instability of the flow of two immiscible liquids with different viscosities in a pipe. Journal of Fluid Mechanics, 141, pp.309-317. [8] Kouris, C. and Tsamopoulos, J., 2001. Dynamics of axisymmetric core-annular flow in a straight tube. I. The more viscous fluid in the core, bamboo waves. Physics of Fluids, 13(4), pp.841-858. [9] C. Kouris and J. Tsamopoulos, Dynamics of axisymmetric core-annular flow in a straight tube. II. The less viscous fluid in the core, saw tooth waves, Phys. Fluids 14, 1011 (2002) [10] Bogaerds, A.C., Hulsen, M.A., Peters, G.W. and Baaijens, F.P., 2003. Time dependent finite element analysis of the linear stability of viscoelastic flows with interfaces. Journal of non-newtonian fluid mechanics, 116(1), pp.33-54. [11] Le Meur, H., 1997. Non-uniqueness and linear stability of the one-dimensional flow of multiple viscoelastic fluids. ESAIM: Mathematical Modelling and Numerical Analysis, 31(2), pp.185-211. [12] Taghavi, S.M., Alba, K. and Frigaard, I.A., 2012. Buoyant miscible displacement flows at moderate viscosity ratios and low Atwood numbers in near-horizontal ducts. Chemical engineering science, 69(1), pp.404-418. [13] Alba, K., Taghavi, S.M. and Frigaard, I.A., 2014. Miscible heavy-light displacement flows in an inclined two-dimensional channel: a numerical approach. Physics of Fluids, 26(12), p.122104. [14] Taghavi, S.M., Alba, K., Son, T., Wielage-Burchard, K., Martinez, D.M. and Frigaard, I.A., 2012. Miscible displacement flows in near-horizontal ducts at low Atwood number. Journal of Fluid Mechanics, 696, pp.175-214. [15] Petitjeans, P. and Maxworthy, T., 1996. Miscible displacements in capillary tubes. Part 1. Experiments. Journal of Fluid Mechanics, 326, pp.37-56. [16] Hasnain, A., Segura, E. and Alba, K., 2017. Buoyant displacement flow of immiscible fluids in inclined pipes. Journal of Fluid Mechanics, 824, pp.661-687. [17] Klubertanz, G., Bouchelaghem, F., Laloui, L. and Vulliet, L., 2003. Miscible and immiscible multiphase flow in deformable porous media. Mathematical and Computer Modelling, 37(5-6), pp.571-582. [18] Barrios-Pia, H., Ramrez-Len, H., Rodrguez-Cuevas, C. and Couder-Castaeda, C., 2014. Multilayer numerical modeling of flows through vegetation using a mixing-length turbulence model. Water, 6(7), pp.2084-2103. [19] Li, J. and Renardy, Y., 2000. Numerical study of flows of two immiscible liquids at low Reynolds number. SIAM review, 42(3), pp.417-439. [20] Redapangu, P.R., Chandra Sahu, K. and Vanka, S.P., 2012. A study of pressure-driven displacement flow of two immiscible liquids using a multiphase lattice Boltzmann approach. Physics of Fluids, 24(10), p.102110. [21] I. Podlubny, Fractional Differential Equations, Vol. 198 of Mathematics in Science and Engineering, Academic Press, San Diego, USA, 1999. [22] M. Caputo, Linear model of dissipation whose Q is almost frequency independent II, The Geophys. J. R. A. Soc., 13 (1967) 529-539. [23] M. Caputo, M. Fbrizio, A new definition of fractional derivative without singular kernel, Progr. Fract. Differ. Appl., 1 (2) (2015) 73-85. [24] A. Atangana, D. Baleanu, New fractionalderivatives with non-local and non-singular kernel: Theory and application to heat transfer model, Therm. Sci., 20 (2) (2016) 763-769.

26

[25] Yang, X.J., Srivastava, H.M. and Machado, J.A., 2015. A new fractional derivative without singular kernel: application to the modelling of the steady heat flow. arXiv preprint arXiv:1601.01623. [26] Yang, X.J., Gao, F., Machado, J.T. and Baleanu, D., 2017. A new fractional derivative involving the normalized sinc function without singular kernel. The European Physical Journal Special Topics, 226(16-18), pp.3567-3575. [27] Xiao-Jun, Yang., Abdel-Aty, M. and Cattani, C., 2019. A new general fractional-order derivataive with Rabotnov fractional-exponential kernel applied to model the anomalous heat transfer. Thermal Science, 23(3). [28] Yang, X.J. and Machado, J.T., 2017. A new fractional operator of variable order: application in the description of anomalous diffusion. Physica A: Statistical Mechanics and its Applications, 481, pp.276283. [29] Yang, X.J., 2019. General Fractional Derivatives: Theory, Methods and Applications. Chapman and Hall/CRC. [30] Yang, X.J., Gao, F., Ju, Y. and Zhou, H.W., 2018. Fundamental solutions of the general fractionalorder diffusion equations. Mathematical Methods in the Applied Sciences, 41(18), pp.9312-9320. [31] Yang, X.J., Feng, Y.Y., Cattani, C. and Inc, M., 2019. Fundamental solutions of anomalous diffusion equations with the decay exponential kernel. Mathematical Methods in the Applied Sciences, 42(11), 4054-4060. [32] Yang, X.J., 2018. New rheological problems involving general fractional derivatives with nonsingular power-law kernels. Proceedings of the Romanian Academy Series A-Mathematics Physics Technical Sciences Information Science, 19(1), pp. 45-52. [33] Yang, X.J., Gao, F. and Srivastava, H.M., 2018. A new computational approach for solving nonlinear local fractional PDEs. Journal of Computational and Applied Mathematics, 339, pp.285-296. [34] J. Hristov, Transient space-fractional diffusion with a poer-law superdiffusivity: Approximate integralbalance approach, Fundamenta Informaticae, 151 (2017) 371-388, doi: 10.3233/FI-2017-1498. [35] Y. Mahsud, N. A. Shah, D. Vieru, Influence of time-fractional derivatives on the boundary layer flow of Maxwell fluids, Chinese J. Phys., 55 (20170, 1340-1351. [36] S. Wang, M. Zhao, X. Li, Transient electro-osmotic flow of generalized Maxwell fluids in a straight pipe of circular cross section, Cent. Eur. J. Phys. 12 (6) (2014) 445-451. [37] J. Hristov, Integral-balance approach to 1-D space-fractional diffusion models, Mathematical Methods in Engineering Applications in Dynamics of Complex Systems, Springer, 2018. [38] J. Hristov, Integral-Balance Solution to Nonlinear Subdiffusion Equation, Ch. 3, in Frontiers in Fractional Calculus, 1 (2017), 71-106. [39] J. Hristov, A transient flow of a non-Newtonian fluid modeled by a mixed time-space derivative: An improved integral-balance approach, pp. 1-20, in Mathematical Methods in Engineering-Theory, Edited by KenanTas, DumitruBaleanu and J. A. T. Machado, Springer, 2018. [40] Y. Nakayama, R.F. Boucher, Introduction to fluid Mechanics, Butterworth Heinemann,1999. [41] C. F. Lorenzo, T. T. Hartley, Generalized functions for the fractional calculus, NASA/TP-1999209424/REV1. [42] M. Arshad, J.Choi, S. Mubeen, K. S. Nisar, G. Rahman, A new extension of the Mittag-Lefler function, Commun. Korean Math. Soc. 33(2) 2018 549-560, doi; 104134/CKMSC170216. [43] L Debnath, D. Bhatta, Integral Transforms and Their Applications, Second Edition, Chapman & Hall/CRC, Boca Raton, London, New York, 2007.

27