- Email: [email protected]
Multi-layer flows of immiscible fractional Maxwell fluids with generalized thermal flux
Multi-layer flows of immiscible fractional Maxwell fluids with generalized thermal flux.
Journal Pre-proof
Multi-layer flows of immiscible fractional Maxwell fluids with generalized thermal flux. Abdul Rauf, Yasir Mahsud, Itrat Abbas Mirza, Qammar Rubbab PII: DOI: Reference:
S0577-9073(19)30946-3 https://doi.org/10.1016/j.cjph.2019.10.006 CJPH 969
To appear in:
Chinese Journal of Physics
Received date: Revised date: Accepted date:
20 December 2018 14 June 2019 9 October 2019
Please cite this article as: Abdul Rauf, Yasir Mahsud, Itrat Abbas Mirza, Qammar Rubbab, Multi-layer flows of immiscible fractional Maxwell fluids with generalized thermal flux., Chinese Journal of Physics (2019), doi: https://doi.org/10.1016/j.cjph.2019.10.006
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 with damped shear stresses are studied. • Memory effects are considered in the thermal transport process. • For large values of the time t, the fluid temperature is smaller in fractional fluids. • Fluids characterized by the fractional constitutive equations flow slower than the ordinary fluids.
1
Multi-layer flows of immiscible fractional Maxwell fluids with generalized thermal flux.
Abdul Rauf 1 ∗ Yasir Mahsud 2 Itrat Abbas Mirza 3 Qammar Rubbab 4 1 Dept. of Computer Science and Engineering, Air University Multan, Pakistan 2 Abdus Salam School of Mathematical Sciences, GCU Lahore, Pakistan 3 Dept. of Mathematics, Khwaja Fareed University of Engineering & information Technology, Rahim Yar Khan, Pakistan 4 Dept. of Mathematics, The Woman University Multan, Pakistan
Abstract Unsteady laminar flows and heat transfer of n-immiscible fractional Maxwell fluids in a channel are investigated under influence of time-dependent pressure gradient. The isothermal channel walls have translational motions in their planes with time-dependent velocities. Governing equations of the mathematical model are based on the generalized constitutive equations for shear stress and thermal flux described by the time-fractional Caputo derivative. Analytical and semi-analytical solutions for velocity, shear stress, and temperature fields are obtained by using finite sine-Fourier and Laplace transforms. In the case of semi-analytical solutions, the inverse Laplace transforms are obtained numerically by employing the Talbots algorithms. Using the software Mathcad, numerical calculations have carried out and results are presented in graphical illustrations in order to analyze the memory effects on the fluid temperature and motion. It is found that in fluids with thermal memory the heat transfer is slower compared with the ordinary fluid, while the fractional velocity parameters act as braking/accelerating factors of the fluids. Keywords: n-layered immiscible fluids; generalized thermal flux; fractional Maxwell fluids; analytical and semi-analytical solutions
1
Introduction
Flows of layered immiscible fluids in a pipe/channel are important for many industrial and physiological applications. For example, the study of two-layer flows of immiscible fluids ∗
Address correspondence to Abdul Rauf. E-mail: attari [email protected]; abdul.rauf @aumc.edu.pk
2
could be useful to the optimal transport of crude oils, polymer solutions or mucus transport in the respiratory tract [1]. To obtain polymer films or plaques, the co-extrusion process of different polymers is used. Sometimes, due to technological process imperfections, the finished product may have some wavy interfaces (interfacial instability). Theoretical studies of the stability of multi-layer plane Poiseuille flows give important information about parameters that control the stability of multi-layer flows such as layer thickness ratios (flow-rate ratios), viscosity ratios or, density ratios [2, 3]. Anturkar et al. [4] have investigated the linear interfacial stability of multilayer plane Poiseuille flows of Oldroyd-B fluids with shear rate dependent viscosity. They determined asymptotic solutions at long wavelengths and found that multilayer flows of viscoelastic fluids can be stable for all wavelengths. Also, it is obtained that, symmetric and nearly symmetric configurations with three layers are unstable if the inner fluid is more viscous than the outer fluids [4]. Linear stability of flows of viscoelastic fluids with interfaces has been analyzed by Bogaerds et al. [5] 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 [6]. In this interesting work, the existence and uniqueness for the multilayer Couette/ Poiseuille flows in channels/pipes 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 solution. Vasquez et al. [7] have modeled and simulated the mucus flow in human bronchial epithelial cell cultures. They studied idealized axisymmetric swirling flows considering four kinds of fluids Newtonian fluid, upper-convected Maxwell fluid, Giesekus fluid and 5-mode Giesekus fluid. Barannyk et al [8] have studied nonlinear dynamic of stratified immiscible, perfect dielectric fluids in horizontal channels in the presence of an electric field perpendicular on the gravitational field. They found that a strong electric field could stabilize the Rayleigh -Taylor interfacial instability. Flows of three-layer of different fluids driven by the gravity and a pressure gradient in a slightly inclined channel have been investigated by Papaefthymiou and Papageorgiou [9]. They observed that nonlinear stability depends on the existence of bounded invariant regions in the hyperbolic state space of the interfacial. Zhao et al. [10] have studied unsteady boundary layer natural convection heat transfer of the fractional Maxwell fluids. Mainly, they found that these coefficients are several times larger than those in non-swirling flows. The interfacial coefficient is well correlated in terms of the liquid volume fraction and the gas Reynolds number, while both gas Reynolds number and liquid Reynolds number are required for correlating the wall friction coefficient. Other interesting problems can be found in [11–14]. It is worth point out that many complex systems in chemistry, mechanics or biology are modeled 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 are described more accurate by the fractional differential equations. There are many fractional differential operators studied by the fractional calculus. Riemann-Liouville fractional integral/derivative [15], Caputo fractional derivative [16], Caputo-Fabrizio fractional derivative [17], AtagnanaBaleanu fractional derivative [18] or Prabhakar fractional derivative [19] are only some from the fractional differential operators used in viscoelasticty, thermal and mass transport. Hristov [20] has applied the Riemann-Liouville fractional derivative to study the transient space-fractional diffusion with power-law superdiffusivity. Mahsud et al. [21] have 3
investigated boundary layer flows of an upper-convected Maxwell fluid modeled by the timefractional Caputo-Fabrizio derivatives. They found that fractional fluids flow slower than the ordinary fluids. Using Laplace and Hankel transforms, Jiang et al. [22] have studied electroosmotic flows of Oldroyd-B fluids with a fractional constitutive equation based on the time-fractional Riemann-Liouville derivative. Wang et al. [23] have investigated the unsteady electroosmotic flows of generalized Maxwell fluids in cylindrical domains and made comparison between ordinary and generalized Maxwell fluids. Ying et al. [24] have studied the electroosmotic flows of the fractional Maxwell fluids in a rectangular microchannel using a fully discrete spectral method for the temporal variable and a Legendre spectral method in the spatial direction. Theoretical results and practical applications of the fractional calculus can be found in the references [25–28]. In this paper we investigate multi-layer flows with heat transfer of Maxwell fluids with generalized constitutive equations for the shear stress and thermal flux. The flow domain is a rectangular channel bounded of two parallel isothermal plates which have translational motion in their planes with time-dependent velocities. The generalized constitutive equations, used in this work, are described by the Caputo time-fractional derivative, therefore the histories of the velocity gradient and histories of temperature gradient influence the fluid shear stress and temperature. Analytical solutions for the Laplace transforms of temperature, velocity and shear stress fields are obtained by employing the Laplace and finite sine-Fourier transforms and hypothesis that interfacial velocities, shear stresses, thermal fluxes and temperatures are equal for adjacent fluids. The obtained results have a general character; therefore many particular cases can be obtained. Multi-layer flows of ordinary/fractional Maxwell fluids are studied as particular cases by letting some fractional parameters equal with one. Numerical values of the temperature and velocity fields are obtained by using the Mathcad software and the numerical Talbots algorithms for the inverse Laplace transforms. Numerical results are used for graphical illustrations of the fluid temperature and velocity. According with the generalized constitutive equations for the shear stress and thermal flux, the variations of the thermal/velocity kernels with the fractional parameter are different for small, respectively large time. It is found that at small values of time t, the heat transfer is smaller in fluids with generalized thermal flux than in the ordinary. At large values of time t, the behavior is opposite. The velocity memory fractional parameters act as braking/accelerating factors of the fluids.
2
Mathematical Modeling
The flow domain is D0 = {(x0 , y 0 , z 0 )| x0 ∈ (−∞, ∞), y 0 ∈ [0, h], z 0 ∈ (−∞, ∞)}, with the boundary walls situated in planes y 0 = 0 and y 0 = h > 0 . Initially at time t0 = 0, the two plates and the enclosed fluids are at rest with the ambient temperature is T2o . After this moment, the channel flat wall situated in the plane y 0 = 0 is moving along the x0 -axis with the velocity u01o = U0 f1 (t0 ) and the temperature T1o , while the plane wall y 0 = h translates parallel with the x0 -axis with the velocity u02o = U0 g1 (t0 ) and the temperature T2o (Fig.1). Functions f1 (t0 ), g1 (t0 ) are differentiable functions of exponential order at infinity and f1 (0) = g1 (0) = 0. In this paper the velocity field is considered to be of the form 4
Figure 1: Geometry of the problem. V~ 0 = (u0 (y 0 , t0 ), 0, 0). We consider here n-layers Maxwell fluids between two parallel plates. Let h0 = 0 and hn = h. In the region y 0 ∈ [hi−1 , hi ], hi−1 < hi flows a Maxwell fluid with the density ρi , viscosity µi , relaxation time λ0i = Gµii , Gi the elastic modulus, velocity u0i (y 0 , t0 ), temperature Ti0 (y 0 , t0 ), thermal flux qi0 (y 0 , t0 ) and the shear stress τi0 (y 0 , t0 ), where i = 1, 2, . . . , n. The fluids are assumed to be immiscible and incompressible and the flow is one-dimensional, unsteady and fully developed. The flows of the fluids are generated by the time-dependent pressure gradient in the flow direction and by the movements of the two walls of the channel. Based on the assumptions about velocity fields, the continuity equation is identically 0 := satisfied by all velocities ui (y 0 , t0 ), i = 1, 2, · · · , n. We fix here In1 := {1, 2, · · · , n}, In−1 {0, 1, · · · , n − 1}. The governing equations of motion along with initial, boundary and the interface fluid-fluid conditions are • the system of linear momentum equations ρi
∂u0 i ∂τ 0 i ∂p0 = − 0 , i ∈ In1 , ∂t0 ∂y 0 ∂x
(1)
• the system of constitutive equations τi0 + λ0i
∂τ 0 i ∂u0 i = µ , i ∈ In1 , i ∂t0 ∂y 0
(2)
• the initial conditions u0i (y 0 , 0) = 0, τi0 (y 0 , 0) = 0, i ∈ In1 , 5
(3)
• the boundary condition on the bottom plate y 0 = 0, u01 (0, t0 ) = U0 f1 (t0 ),
(4)
on the upper plate y 0 = hn = h, u0n (h, t0 ) = U0 g1 (t0 ), • the interfaces conditions 0 1 u0i (hi , t0 ) = u0i+1 (hi , t0 ), τi0 (hi , t0 ) = τi+1 (hi , t0 ), i ∈ In−1 .
(5)
The system of governing equations for thermal transport and heat flux are given by • the balance equations of thermal transport ρi cip
∂Ti0 (y 0 , t0 ) ∂qi0 (y 0 , t0 ) = − , ∂t0 ∂y 0
• the Fourier’s law of thermal flux qi0 (x0 , t0 ) = −ki
(6)
∂Ti0 (y 0 , t0 ) , ∂y 0
(7)
where i ∈ In1 , cip and ki are the specific heat at constant pressure and thermal conductivity, respectively. To the thermal transport we consider the following initial, boundary and the interface fluid-fluid conditions: • initial conditions Ti0 (y 0 , 0) = T2o , qi0 (y 0 , 0) = q0 , i ∈ In1 ,
(8)
T 0 (0, t0 ) = T1o , Tn0 (h, t) = T2o ,
(9)
• boundary conditions
• interface conditions 0 0 1 Ti0 (hi , t0 ) = Ti+1 (hi , t0 ), qi0 (hi , t0 ) = qi+1 (hi , t0 ), i ∈ In−1 .
(10)
Introducing the non-dimensional variables 0 0 0 0 0 x0 , y = yh , t = νh12t , ui = uU0i , τi = µhτ1 Ui0 , νi = µρii , p = µhp , h U 1 0 2 2 0 λi = λhi2ν1 , ai = ρρ1i , bi = µµ1i , f (t) = f1 hν1t , g(t) = g1 hν1t , di Ti0 −T2o q0 µc Ti = T1o−T , qi = q0i , σi = νν1i , P ri = ikiip , i ∈ In1 , 2o
x=
=
hi , h
into Eqs. (1)-(10), the mathematical model for dimensionless governing equations takes the form, ai
∂ui ∂p ∂τi =− + , i ∈ In1 , ∂t ∂x ∂y 6
(11)
τi + λi
∂τi ∂ui = bi , i ∈ In1 , ∂t ∂y
ν1 γ0 ρi cip
∂Ti ∂qi = − , i ∈ In1 , ∂t ∂y
qi = −γ0 ki
∂Ti , i ∈ In1 , ∂y
(12)
(13)
(14)
along with • the non-dimensional initial conditions ui (y, 0) = 0, τi (y, 0) = 0, Ti (y, 0) = 0, qi (y, 0) = 1, i ∈ In1 ,
(15)
• the non-dimensional boundary condition u1 (0, t) = f (t), un (1, t) = g(t), T1 (0, t) = 1, Tn (1, t) = 0,
(16)
• the non-dimensional interface conditions ui (di , t) = ui+1 (di , t), τi (di , t) = τi+1 (di , t),
(17)
1 Ti (di , t) = Ti+1 (di , t), qi (di , t) = qi+1 (di , t), i ∈ In−1 ,
where γ0 =
2.1
T1o −T2o . q0 h
Mathematical model with fractional constitutive equations
In the following , we consider the generalized constitutive equations, τi + λi Dαt i τi = bi
qi =
−γ0 ki Dt1−βi
Dγt F (y, τ )
1 = Γ(1 − γ)
∂ui , αi ∈ (0, 1], i ∈ In1 , ∂y
(18)
∂Ti , βi ∈ (0, 1], i ∈ In1 , ∂y
(19)
where Z
0
t
(t − τ )−γ
∂F (y, τ ) dτ, 0 ≤ γ < 1 ∂τ
(20)
(y,t) is the Caputo derivative. For particular case γ = 1, D1t F (y, t) = ∂F∂t . It is important to note that the fractional constitutive equations (18) and (19) can be written in the equivalent forms bi αi −1 tαi ∂ui (y, t) τi (y, t) = t Eαi ,αi − ∗ (21) λi λi ∂y
7
=
Z
0
t
bi (t − τ )αi −1 Eαi ,αi λi
(t − τ )αi − λi
tβi −1 ∂ 2 Ti (y, t) qi (y, t) = ∗ = Γ(βi ) ∂t∂y
Z
t
0
∂ui (y, τ ) dτ, i ∈ In1 , ∂y
(t − τ )βi −1 ∂ 2 Ti (y, t) dτ, i ∈ In1 . Γ(βi ) ∂t∂y
(22)
The Eqs. (21) and (22) highlight the fact that the histories of the velocity gradient and temperature gradient influence the time-evaluation of shear stress and thermal flux, respectively. Also it is observed from theaboveequations that, the nonlocality kernel for the shear stress is the function
bi αi −1 t Eαi ,αi λi
α
− tλii , Eαi ,αi (.) being the Mittag-Leffler function; while
the nonlocality kernel of the thermal flux is the power-law function the pressure gradient is a known function of t, therefore −
tβi −1 . Γ(βi )
We suppose that
∂p = P (t), ∂x
(23)
where P (t) is a local integrable function on [0, ∞) of the exponential order at infinity.
3
Solution of the problem
In order to determine analytical solutions of the Eqs. (11), (13, (18), (19) with conditions (15)-(17) we use the Laplace transform coupled with the finite sine-Fourier transform [29]. Applying the Laplace transform to Eqs. (11), (13), (16)-(19) and using the initial conditions (15), we obtain the problem in transformed domain ∂ τ¯i (y, s) , i ∈ In1 , ∂y
(24)
bi ∂ u¯i (y, s) , i ∈ In1 , α 1 + λi s i ∂y
(25)
ai s¯ ui (y, s) = P¯ (s) +
τ¯i (y, s) =
∂ q¯i (y, s) , i ∈ In1 , ν1 γ0 ρi cip sT¯i (y, s) = − ∂y
q¯i = −γ0 ki s1−βi
∂ T¯i (y, s) , i ∈ In1 , ∂y
(26)
(27)
1 u¯1 (0, s) = f¯(s), u¯n (1, s) = g¯(s), T¯1 (0, s) = , T¯n (1, s) = 0, s
(28)
u¯j (dj , s) = u¯j+1 (dj , s), τ¯j (dj , s) = τ¯j+1 (dj , s),
(29)
1 T¯j (dj , s) = T¯j+1 (dj , s), q¯j (dj , s) = q¯j+1 (dj , s), j ∈ In−1 .
8
3.1
Solution for the thermal transport
Using Eqs. (26) and (27) we can write the Laplace transformed system of equations for thermal transports T¯i (y, s), i ∈ In1 , in the form ∂ 2 T¯i = ζi sβi T¯i , i ∈ In1 , ∂y 2 where ζi =
P ri . σi
(30)
The general solution of Eq. (30) is T¯i (y, s) = Ai (s)e−y
√
ζ i sβ i
+ Bi (s)ey
√
ζ i sβ i
, i ∈ In1 ,
(31)
where the unknown parameters A1 , A2 . . . , An , B1 , B2 , . . . , Bn are to determined by applying the boundary conditions (28) and interface conditions (29) to Eq. (31) and are given by the following linear system M(s)X(s) = C(s), where the matrix M(s) is given by 1 0 0 ··· 0 a −b 0 ··· 0 1 1 0 a −b · · · 0 2 2 .. . .. . . .. .. .. . . 0 0 · · · an−1 −bn−1 M(s) = p −q 0 ··· 0 1 1 0 p −q · · · 0 2 2 .. . .. ... ... .. . . 0 0 · · · pn−1 −qn−1 0 0 ··· 0 pn
and the columns X(s) and C(s) are
Here
−dj
aj (s) = e
√
−2
pn (s) = e `j (s) =
kj+1 ,j kj
X(s) =
√
ζ j sβ j
ζn sβn
An B1 B2 .. . Bn
−dj
, bj (s) = e
, qj (s) = −`j s
1 ∈ In−1 .
A1 A2 .. .
√
1 1 a1
0 .. .
(32)
0 − b11
βj −βj+1
1 a2
··· ··· ··· .. .
1 s
0 0 .. . 0
,
.
√ β p ζ s j β j , pj (s) = − ζj s aj (s), rj (s) = aj j(s) ,
p bj (s) ζj+1 sβj+1 , sj (s) = 9
0 0 0 .. .
.. . 1 1 − bn−1 0 0 · · · an−1 r1 −s1 0 ··· 0 0 r2 −s2 · · · 0 .. .. .. .. .. . . . . . 0 0 · · · rn−1 −sn−1 0 0 ··· 0 1
, C(s) =
ζj+1 sβj+1
0 0 − b12 .. .
√
`j sβj −βj+1 ζj+1 sβj+1 , bj (s)
In order to write the system (32) in a suitable form, we introduce the notations
M =
N =
Q=
R=
1 0 0 ··· 0 a1 −b1 0 ··· 0 0 a2 −b2 · · · 0 .. .. .. ... ... . . . 0 0 · · · an−1 −bn−1
1
0 − b11
1 a1
1 a2
0 .. .
.. . 0
0
0 0 − b12 .. .
··· ··· ··· .. .
···
1
− bn−1
r1 −s1 0 ··· 0 0 r2 −s2 · · · 0 .. .. .. ... ... . . . 0 0 · · · rn−1 −sn−1 0 0 ··· 0 1
C1 =
1 s
0 .. . 0
n×1
, C2 =
1
p1 −q1 0 ··· 0 0 p2 −q2 · · · 0 .. .. .. ... ... . . . 0 0 · · · pn−1 −qn−1 0 0 ··· 0 pn
0 0 .. . 0
; i.e.
Mij = ai δi,j+1 − bi δi+1,j+1 , n×n
0 1 , ,j ∈ In−1 i ∈ In−1
0 , N0j = δ0,j , j ∈ In−1
; i.e.
Nij = n×n
1 δ ai i,j+1
−
(34)
1 δ , bi i+1,j+1
1 0 i ∈ In−1 ,j ∈ In−1 ,
Q = (Qij )i,j∈I 0 , n−1
; i.e.
Qij = pi+1 δi+1,j+1 − qi+1 δi+2,j+1 , i∈
n×n
0 ,j In−2
∈
(35)
0 , In−1
0 Q(n−1)j = pn δn−1,j , j ∈ In−1 ,
R = (Rij )i,j∈I 0 ,
n−1
; i.e.
Rij = ri+1 δi+1,j+1 − si+1 δi+2,j+1 , i∈
n×n
0 In−2 ,j
,A =
A1 A2 .. . An
n×1
,B =
where δij is the Kronecker tensor. The system (32) is written as M A + N B = C1 , QA + RB = C2 . 10
∈
(36)
0 In−1 ,
0 R(n−1)j = δn−1,j , j ∈ In−1 ,
n×1
(33)
n−1
0 , M0j = δ0,j , j ∈ In−1
N = (Nij )i,j∈I 0 ,
n−1
0 0 0 .. .
an−1
M = (Mij )i,j∈I 0 ,
B1 B2 .. . Bn
,
(37)
n×1
(38)
Matrices M, N, Q, R are triangular and invertible matrices, Eqs. (38) can be written as (N −1 M − R−1 Q)A = N −1 C1 , (39) B = −R−1 QA. Suppose that matrix S = N −1 M − R−1 Q is invertible, we have A = S −1 N −1 C1 , B = −R−1 QS −1 N −1 C1 .
(40)
An easy computation shows that
1 C0 := N −1 C1 = s
where,
N −1
that is,
=
1
1 b1 a1 b1 b2 a1 a2
.. .
b1 b2 ···bn−1 a1 a2 ···an−1
b1 a1 b1 b2 a1 a2
.. .
0 −b1 − b1a2b2 .. .
0 0 −b2 .. .
b1 b2 ···bn−1 a1 a2 ···an−1
···bn−1 − b1a2b2···a n−1
···
··· ··· ··· .. .
(41)
0 0 0 .. .
bn−1 − bn−2 −bn−1 an−1
,
n×n
0 N −1 = (nij )i,j∈In−1 ,
0 , ni0 = n0j = δ0,j , j ∈ In−1
nij = −
i Q
m=1
bm am
This reduces the system (40) to
j−1 Q `=0
a`+1 b`
n−1 P k=0
i Q
`=1
b` ,i a`
1 , ∈ In−1
(42)
1 δi,j+k , b0 = 1, i, j ∈ In−1 .
A = S −1 C0 , B = −R−1 QS −1 C0 .
(43)
The inverse matrices of M and R are 0 M −1 = (mij )i,j∈In−1
0 m0j = δ0,j , j ∈ In−1 , mi0 =
mij = −
i Q
m=1
am bm
j−1 Q `=0
b` a`+1
n−1 P k=0
i Q
`=1
a` ,i b`
1 ∈ In−1 ,
1 δi,j+k , b0 = 1, i, j ∈ In−1 ;
11
(44)
0 R−1 = (ρij )i,j∈In−1 ,
ρij =
n−1 P k=0
i Q
m=0
rm sm
j Q
`=0
s` r`+1
(45) δi−k,j , r0 = s0 = rn = 1, i, j ∈
0 In−1 .
0 Now, the matrix S = (Sij )i,j∈In−1 is defined by the elements
Sij =
n−1 X k=0
(nik Mkj − ρik Qkj ).
(46)
After having determined the parameters Aj , Bj , j ∈ In1 , from the system (43) the analytic solution for the system of thermal transport Ti (y, t) can be determined by applying the inverse Laplce transform to Eq. (31) 1 Ti (y, t) = 2πι
Z
σ+i∞
est (Ai (s)e−y
√
ζi sβi
+ Bi (s)ey
√
σ−i∞
ζi sβi
)ds, i ∈ In1 .
(47)
In this paper, the inverse Laplace transforms are obtained by using Talbot’s algorithms [31, 32]. ¯ s) be the Laplace transform of the function g(y, t). The function g(y, t) can be Let G(y, approximated with the help of Talbot algorithm [31] for the Laplace transform inversion, ( ) M −1 X r exp(rt) ¯ r) + ¯ z(θk )) (1 + iζ(θk )) , g(y, t) ∼ G(y, Re exp(tz(θk ))G(y, (48) = M 2 k=1 where r=
2M , 5t
z(θ) = rθ(cot θ + i), θ ∈ (−π, π),
ζ(θ) = θ + (θcotθ − 1)cotθ, θk =
(49)
kπ . M
Another method to approximate the function g(y, t) is the improved Talbot algorithm for the inverse Laplace transform and is given by [32] M
g(y, t) ∼ =
1X ¯ z1 (σk )) (ν + iζ 1 (σk )), exp(tz1 (σk ))G(y, t k=1
(50)
where, z1 (θ) =
M t
[νiθ + µθ cot(αθ) − ξ] , θ ∈ [−π, π],
ζ1 (θ) = αµθ + µ(αθcot(αθ) − 1)cot(αθ), σk =
(2k−1)π M
Here α, M, ν, µ, ξ are parameters to be specified by the user.
12
(51) − π.
3.2
Analytical solutions for velocity and shear stress
Multiplying Eq. (24) by 1 + λi sαi both sides and using Eq. (25) we obtained the system of equations for Laplace transformed velocities, ∂ 2 u¯i (y, s) , i ∈ In1 , ai s(1 + λi sαi )¯ ui (y, s) = (1 + λi sαi )P¯ (s) + bi ∂y 2 ¯ s) = where φ(y,
R∞
(52)
φ(y, t) exp(−st)dt is the Laplace transform of the function φ(y, t) [29].
0
In order to find the exact solution of velocities ui (y, t) we will use the finite sine-Fourier transform along with the Laplace transform. The finite sine-Fourier transform of the function ¯ s), y ∈ [a, b] is defined as [30] φ(y, e (s) = φ m
Zb a
¯ s) sin(ϑm (y − a))dy, ϑm = mπ , b > a, m = 1, 2, · · · , φ(y, b−a
(53)
along with the inverse Fourier transform defined by ¯ s) = φ(y,
∞ 2 X e (s). sin(ϑm (y − a))φ m b − a m=1
(54)
With the application of finite Fourier sine transform (53) to Eq. (52) along with the boundary conditions (28) and (29), the transformed velocities take the form eim (s) = u
1−(−1)m
(i)
(1+λi sαi )P¯ (s)
(i) ϑm
ϑm =
sαi )+b
ai s(1+λi mπ ,m = di −di−1
(i) 2 i (ϑm )
+
1, 2, · · · , i
(i)
bi ϑm u ¯i (di−1 ,s) ai s(1+λi ∈ In1 ,
sαi )+b
(i) 2 i (ϑm )
−
(i)
(−1)m bi ϑm u ¯i (di ,s) (i) ai s(1+λi sαi )+bi (ϑm )2
,
(55)
where for i = 1, u¯1 (d0 , s) = f¯(s) and for i = n, u¯n (dn , s) = g¯(s). In order to apply the inverse Fourier sine transform, we rewrite the Eq. (55) in the following suitable form (−1)m+1 u¯i (di , s) u¯i (di−1 , s) 1 − (−1)m (1 + λi sαi )P¯ (s) e − uim (s) = + + (i) (i) (i) (i) ϑm ϑm ϑm ai s(1 + λi sαi ) + bi (ϑm )2 ai s(1 + λi sαi )¯ ui (di−1 , s) (i)
(i)
ϑ(i) m
(−1)m ai s(1 + λi sαi )¯ ui (di , s)
, (i) (i) ϑm (ai s(1 + λi sαi ) + bi (ϑm )2 ) mπ = , m = 1, 2, · · · , i ∈ In1 . di − di−1
ϑm (ai s(1 + λi sαi ) + bi (ϑm )2 )
+
(56)
Using the following auxiliary functions along with their inverse Fourier sine transform ϕ1i (y) =
di −y , di −di−1 y−di−1 , di −di−1
y ∈ [di−1 , di ], ϕ˜1im =
ϕ2i (y) = y ∈ [di−1 , di ], ϕ˜2im = m = 1, 2, · · · , i ∈ In1 , d0 = 0, dn = 1,
13
1 (i) , ϑm (−1)m+1 (i)
ϑm
,
(57)
the inverse Fourier sine transform of Eq. (56) takes the form i−1 i −y u¯i (y, s) = did−d u¯i (di−1 , s) + dy−d u¯i (di , s)+ i−1 i −di−1 ∞ m P (i) (1+λi sαi )P¯ (s) 2 sin(ϑm (y − di−1 )) 1−(−1) (i) 2 − (i) αi di −di−1
2 di −di−1 2 di −di−1
m=1 ∞ P m=1 ∞ P
ϑm
(i)
sin(ϑm (y − di−1 )) (i) sin(ϑm (y
m=1 (i)
ϑm =
− di−1 ))
mπ ,m di −di−1
ai s(1+λi s
)+bi (ϑm )
ai s(1+λi sαi )¯ ui (di−1 ,s) (i) (i) ϑm (ai s(1+λi sαi )+bi (ϑm )2 ) ui (di ,s) (−1)m ai s(1+λi sαi )¯ (i)
(i)
ϑm (ai s(1+λi sαi )+bi (ϑm )2 )
+
(58)
,
= 1, 2, · · · , i ∈ In1 .
Now, from (25) and (58), we obtain τ¯i (y, s) =
bi ∂ u¯i (y, s) = T¯i1 (y, s)¯ ui (di , s) − T¯i2 (y, s)¯ ui (di−1 , s) + T¯i3 (y, s)P¯ (s), (59) 1 + λi sαi ∂y
where ∞ (i) m αi X b (−1) a (1 + λ s )s cos(ϑ (y − d )) m i i i i−1 T¯i1 (y, s) = 1+2 , αi ) + b (ϑ(i) )2 ) (1 + λi sαi )(di − di−1 ) (a s(1 + λ s m i i i m=1 T¯i2 (y, s) =
∞ (i) X bi ai (1 + λi sαi )s cos(ϑm (y − di−1 )) 1+2 , αi ) + b (ϑ(i) )2 ) (1 + λi sαi )(di − di−1 ) (a s(1 + λ s m i i i m=1
T¯i3 (y, s) =
2bi di − di−1
X ∞
(i)
(1 − (−1)m ) cos(ϑm (y − di−1 ))
m=1
ai s(1 + λi
sαi )
+
(i) bi (ϑm )2
, i ∈ In1 .
(60)
With the help of the interface liquid-liquid conditions (29a), the Laplace transformed velocities on the interfaces y = di , i = 1, 2, . . . , n − 1, can be written as ¯ i1 (s) L 1 , i ∈ In−1 , u¯i (di , s) = u¯i+1 (di , s) = ¯ Li0 (s)
(61)
where, ¯ i0 (s) = T¯i1 (di , s) − T¯(i+1)1 (di , s), L ¯ i1 (s) = (T¯i2 (di , s) − T¯(i+1)2 (di , s)))¯ L ui (di−1 , s) + (T¯(i+1)3 (di , s) − T¯i3 (di , s))P¯ (s). In order to obtained the inverse Laplace transforms of the functions u¯i (y, s), i ∈ In1 and 1 u¯j (dj , s) = u¯j+1 (dj , s), j ∈ In−1 , we consider the following auxiliary functions, ¯ i0 (m, s) = H
=
1 (i)
ai s(1 + λi sαi ) + bi (ϑm )2
s−1 · ai λi (sαi + λ−1 i ) 1+ 14
1 (i)
bi (ϑm )2 s−1 ai λi (sαi +λ−1 i )
∞ (i) 2 k s−k−1 1 X k bi (ϑm ) = (−1) , k+1 ai λi k=0 ai λ i (sαi + λ−1 i )
(62)
¯ i3 (m, s) = sH ¯ i0 (m, s) H ∞ (i) 2 k s−k 1 X k bi (ϑm ) = (−1) , k+1 ai λi k=0 ai λ i (sαi + λ−1 i ) ¯ i1 (m, s) = H
(63)
(1 + λi sαi ) (i)
ai s(1 + λi sαi ) + bi (ϑm )2
∞ (i) 2 k s−k−1 1X k bi (ϑm ) (−1) , = k ai k=0 ai λi (sαi + λ−1 i )
(64)
and ¯ i2 (m, s) := sH ¯ i1 (m, s) H ∞ (i) 2 k 1X 1 s−k k bi (ϑm ) , = + (−1) k ai ai k=1 ai λ i (sαi + λ−1 i )
(65)
Since the generalized G-Lorenzo-Hartley function is defined by [33] −1
Gα,β,γ (t, σ) = L
X ∞ sβ Γ(k + γ)σ k t(k+γ)α−β−1 , = (sα − σ)γ k!Γ((k + γ)α − β)Γ(γ) k=0
Re(s) > 0, Re(αγ − β) > 0, |
(66)
σ | < 1, sα
and for αi , βi > 0, L
−1
sαi −βi = tβi −1 Eαi ,βi (dtαi ), sαi − d
(67)
¯ i0 (m, s), where Eαi ,βi (.) is the Mittag-Leffler function [34]. The inverse Laplace transform of H ¯ i3 (m, s), H ¯ i1 (m, s) and H ¯ i2 (m, s) takes the form H Hi0 (m, t) =
∞ X k=0
(i) 2 k −1 k (bi (ϑm ) ) (−1) Gαi ,−k−1,k+1 (t, −λi ) , (ai λi )k+1
∞ (i) 2 k X −1 k (bi (ϑm ) ) Hi3 (m, t) = (−1) Gαi ,−k,k+1 (t, −λi ) , k+1 (a i λi ) k=0 15
(68)
(69)
∞ (i) 2 k 1X k bi (ϑm ) Gαi ,−k−1,k (t, −λ−1 Hi1 (m, t) = (−1) i ), ai k=0 ai λi
(70)
∞ (i) 2 k δ(t) 1X k bi (ϑm ) Hi2 (m, t) = + Gαi ,−k,k (t, −λ−1 (−1) i ), ai ai k=1 ai λi
(71)
and
where δ(t) is dirac delta function. Using Eqs. (58), (70) and (71), we obtain for velocities ui (y, t), i ∈ In1 , the following expressions: i−1 i −y ui (y, t) = did−d ui (di−1 , t) + dy−d ui (di , t)+ i−1 i −di−1 ∞ m P (i) 2 Hi1 (m, t) ∗ P (t)− sin(ϑm (y − di−1 )) 1−(−1) (i) di −di−1
2 di −di−1 2 di −di−1
m=1 ∞ P m=1 ∞ P
ϑm
(i)
sin(ϑm (y − di−1 ))
ai (i) Hi2 (m, t) ϑm
(i) sin(ϑm (y
(−1)m ai
m=1 (i)
ϑm =
− di−1 ))
mπ ,m di −di−1
(i)
ϑm
∗ ui (di−1 , t)+
(72)
Hi2 (m, t) ∗ ui (di , t),
= 1, 2, · · · , i ∈ In1 .
Rt where h1 (t) ∗ h2 (t) = 0 h1 (t − τ )h2 (τ )dτ is the convolution product of the functions h1 (t), h2 (t) and ui (di , t) = ui+1 (di , t) is given by ui (di , t) = ui+1 (di , t) =
=
Li1 (t) Li0 (t)
(Ti2 (di , t) − T(i+1)2 (di , t)) ∗ ui (di−1 , t) + (T(i+1)3 (di , t) − Ti3 (di , t)) ∗ P (t) , Ti1 (di , t) − T(i+1)1 (di , t)
where, ∞ X bi −1 αi −1 αi −1 Ti1 (di , t) = ai Hi3 (m, t) , λ t Eαi ,αi (−λi t ) + 2 (di − di−1 ) i m=1 ∞ X bi+1 m −1 αi −1 −1 αi T(i+1)1 (y, t) = λ t Eαi ,αi (−λi+1 t ) + 2 (−1) ai+1 H(i+1)3 (m, t) , (di+1 − di ) i+1 m=1 ∞ X bi m −1 αi −1 −1 αi Ti2 (di , t) = λ t Eαi ,αi (−λi t ) + 2 (−1) ai Hi3 (m, t) , (di − di−1 ) i m=1 ∞ X bi+1 −1 αi −1 −1 αi Eαi ,αi (−λi+1 t ) + 2 ai+1 H(i+1)3 (m, t) , T(i+1)2 (di , t) = λ t (di+1 − di ) i+1 m=1 16
(73)
2bi Ti3 (di , t) = di − di−1 2bi+1 T(i+1)3 (di , t) = di+1 − di
X ∞
X ∞
m=1
m=1
((−1) − 1)Hi0 (m, t) , m
1 . (1 − (−1) )H(i+1)0 (m, t) , i ∈ In−1 m
The system of shear stresses τi (y, t), i ∈ In1 , can be determined by applying inverse Laplace transform to Eq. (59) and using Eqs. (68) and (69) τ¯i (y, s) =
bi ∂ u¯i (y, s) = T¯i1 (y, s)¯ ui (di , s) − T¯i2 (y, s)¯ ui (di−1 , s) + T¯i3 (y, s)P¯ (s), (74) α i 1 + λi s ∂y
where ∞ X bi m −1 αi −1 −1 αi (i) λ t Eαi ,αi (−λi t ) + 2 (−1) ai Hi3 (m, t) cos(ϑm (y − di−1 )) , Ti1 (y, t) = di − di−1 i m=1 ∞ X bi −1 αi −1 αi −1 (i) ai Hi3 (m, t) cos(ϑm (y − di−1 )) , λ t Eαi ,αi (−λi t ) + 2 Ti2 (y, t) = di − di−1 i m=1 ∞ X 2bi (i) Ti3 (y, t) = (1 − (−1)m )Hi0 (m, t) cos(ϑm (y − di−1 )), i ∈ In1 . di − di−1 m=1
3.3
Semi-analytical solution for velocity and shear stress
This subsection is devoted to determine a new form of solutions to the problem given by Eqs. (24), (25), (28a), (29a) with the help of the Laplace transform coupled with the classical method for the ordinary differential equations. The general solutions of the system of the equations (52) are p p 1 ¯ P (s), (75) u¯i (y, s) = Di (s) exp −y wi (s) + Ei (s) exp y wi (s) + ai s where
wi (s) =
ai s(1 + λi sαi ), bi
i ∈ In1 .
(76)
Now with Eq. (25) and (75), the Laplace transformed form of the shear stresses τi (y, t), i ∈ In1 , can be written as τ¯i (y, s) =
h i p p p p bi −D (s) w (s) exp(−y w (s)) + E (s) w (s) exp(y w (s)) . (77) i i i i i i 1 + λi sαi
Using the boundary conditions (28) and the interface conditions (29), we have M1 (s)X1 (s) = C1 (s), 17
(78)
where the matrix M(s) is given by 1 0 0 ··· 0 e −f 0 · · · 0 1 1 0 e −f2 · · · 0 2 .. . .. . . .. .. .. . . 0 0 · · · e −f n−1 n−1 M1 (s) = c −d 0 · · · 0 1 1 0 c −d2 · · · 0 2 . . .. . . .. .. .. .. . 0 0 · · · cn−1 −dn−1 0 0 ··· 0 cn
1 1 e1
0 .. .
0 − f11 1 e2
0 0 − f12 .. .
··· ··· ··· .. .
0 0 0 .. .
.. . 1 1 0 0 · · · en−1 − fn−1 g1 −h1 0 ··· 0 0 g2 −h2 · · · 0 .. .. .. .. .. . . . . . 0 0 · · · gn−1 −hn−1 0 0 ··· 0 hn
and the columns X(s) and C(s) are
−dj
ej (s) = e Here gj (s) = hj (s) =
p
√
D1 D2 .. .
X1 (s) = wj (s)
Dn E1 E2 .. . En −dj
, fj (s) = e
wj (s)edj
√
wj (s)
P0 (s) P1 (s) .. .
Pn−1 (s) , C1 (s) = 0 0 .. . 0
√
wj+1 (s)
b
√
,
.
√ p −dj wj (s) , cj (s) = − wj (s)e , (1+λ sαj )
j+1 j , dj (s) = − bj (1+λ αj+1 ) j+1 s
bj+1 (1+λj sαj ) p wj+1 (s)edj bj (1+λj+1 sαj+1
wj+1 (s)
√ p wj+1 (s)e−dj wj+1 (s) ,
,
¯ ¯ 1 − a1j ) P (s) P0 (s) = f¯ − Pa(s) , Pj (s) = ( aj+1 , j = 1, 2, ..., n − 1. s 1s Similar to the case of thermal transport, the system (78) can be written in equivalent form as, M1 D + N1 E = C3 , (79) Q1 D + R1 E = C4 .
18
where .. . N1 = Nij1 i,j∈I 0 , n−1 .. . 1 0 1 0 .. N0j = δ0,j , j ∈ In−1 , M0j = δ0,j , j ∈ In−1 , . .. Mij1 = ei δi,j+1 − fi δi+1,j+1 , . Nij1 = e1i δi,j+1 − f1i δi+1,j+1 , .. . 0 1 0 1 .. , , j ∈ In−1 i ∈ In−1 , , j ∈ In−1 i ∈ In−1 . .. 1 . R1 = Rij , Q1 = Q1ij i,j∈I 0 , 0 i,j∈In−1 n−1 .. . 1 = gi+1 δi+1,j+1 − hi+1 δi+2,j+1 , Q1ij = ci+1 δi+1,j+1 − di+1 δi+2,j+1 , ... Rij .. 0 0 0 0 . , , j ∈ In−1 , i ∈ In−2 , j ∈ In−1 i ∈ In−2 .. . 0 1 0 .. , = δn−1,j , j ∈ In−1 R(n−1)j , Q1(n−1)j = cn δn−1,j , j ∈ In−1 . M1 = Mij1
C3 =
P0 (s) P1 (s) .. . Pn−1 (s)
n×1
, C4 =
0 i,j∈In−1
0 0 .. . 0
,
n×1
,D =
D1 D2 .. . Dn
n×1
,E =
E1 E2 .. . En
, n×1
where δij is the Kronecker tensor. Matrices M1 , N1 , Q1 , R1 are triangular invertible matrices, therefore Eqs. (79) can be written as (N1−1 M1 − R1−1 Q1 )D = N1−1 C3 , (80) E = −R1−1 Q1 D. Suppose that matrix S1 = N1−1 M1 − R1−1 Q1 is invertible, we have D = S1−1 N1−1 C3 , E = −R1−1 Q1 S1−1 N1−1 C3 ,
(81)
where, 0 M1 −1 = (m1ij )i,j∈In−1 ,
0 m10j = δ0,j , j ∈ In−1 , m1i0 =
m1ij = −
i Q
m=1
em fm
j−1 Q `=0
f` e`+1
n−1 P k=0
i Q
`=1
e` ,i f`
1 ∈ In−1 ,
1 δi,j+k , f0 = 1, i, j ∈ In−1 ;
19
(82)
0 N1−1 = (n1ij )i,j∈In−1 ,
0 n10j = δ0,j , j ∈ In−1 , n1i0 =
nij = −
i Q
m=1
fm em
j−1 Q `=0
e`+1 f`
n−1 P k=0
i Q
`=1
f` ,i e`
1 ∈ In−1 ,
1 ; δi,j+k , f0 = 1, i, j ∈ In−1
0 R1−1 = (ρ1ij )i,j∈In−1 ,
ρ1ij
=
n−1 P k=0
i Q
m=0
gm hm
j Q
`=0
h` g`+1
0 δi−k,j , g0 = h0 = gn = 1, i, j ∈ In−1 ,
and the matrix S1 = (Sij1 )i,j=0,1,...,n−1 is defined by the elements Sij1
=
n−1 X k=0
1 n1ik Mkj − ρ1ik Q1kj .
The Laplace transforms given by Eqs. (75), (77) and (81) are complicated, therefore analytical methods are difficult to use for the inverse Laplace transforms of these functions. Here we used two accuracy numerical algorithms, namely the fixed Talbot algorithm [31] and the improved Talbot algorithm [32] for numerical Laplace transform inversion defined in Eqs. (48)-(51).
4
Numerical results and discussions
Unsteady multi-layer laminar flows with heat transfer of Maxwell fluids within a rectangular channel have been investigated by using fractional constitutive equations for the shear stress and thermal flux. The motion of fluids is generated by the channel walls, which are moving with the time-dependent velocities in their planes, and by the time-dependent pressure gradient in the flow direction. The generalized constitutive equations, used in this work, are described by the Caputo time-fractional derivative, therefore, the histories of the velocity gradient and histories of temperature gradient influence the behavior of fluids. Analytical solutions for the Laplace transforms of temperature and velocity fields have been obtained by employing the Laplace transform and Talbots algorithms for the numerical inversion of the Laplace transform. Closed form for fluids velocity has been determined by coupling the Laplace transform with finite sine-Fourier transforms. Analytical and semianalytical solutions have been obtained along with the hypotheses that interfacial velocities, shear stresses, thermal fluxes and temperatures are equal for adjacent fluids. The obtained results have a general character; therefore many particular cases can be obtained. Multi-layer flows of ordinary/fractional Maxwell fluids can be studied as particular cases by letting some fractional parameters equal with one. Using the software Mathcad, numerical calculations were carried out for the obtained solutions. Numerical results are presented in Figs. 2- 13. 20
Figs. 2 - 5 have been plotted in order to analyze the influence of the thermal fractional parameters βi , i = 1, 2, 3 on the layers temperature. The considered fluids are characterized by the following material parameters ρ1 = 1000, µ1 = 0.05, G1 = 1.2, cp1 = 0.25, k1 = 2.5, d1 = 0.3, ρ2 = 1300, µ2 = 0.2, G2 = 2.4, cp2 = 0.27, k2 = 2.7, d2 = 0.7, ρ3 = 1500, µ3 = 0.6, G3 = 5.0, cp3 = 0.29, k1 = 2.9, d1 = 1.0. Profiles of layers temperature, given by Eq. (31), versus the spatial coordinate y, for different values of the thermal fractional parameters are sketched in Fig. 2 for small time, t = 0.02, and Fig. 3 for large time, t = 2. It is observed from Fig. 2 that the variation of the thermal fractional parameter β1 of the fluid situated in the first layer, has significant influence on the temperature in the first and second layers but, the influence on the temperature from the third layer is insignificant. This behavior is due to the decreasing of temperature in the first and second layer with the fractional parameter β1 . It is known from Eq. (22) that the memory kernel of thermal flux is β−1 the function hβ (t) = tΓ(β) whose graph is given in Fig. 6. It is observed from Fig 6 that, for t = 0.02 and β1 ≥ 0.4 the kernel hβ (t) is decreasing with β1 , therefore the weight function of the thermal gradient decreases (the memory effects are weaker). An opposite behavior has been found for the large time when the memory kernel increases with the fractional parameter (See Fig 6, for t = 2). In this case, memory effects are stronger and the layers temperature increase with the fractional parameter β1 (See Fig. 3). Time-evolution of layers temperature, in different positions of the channel and for different values of the thermal fractional parameters are presented in Figs. 4 and 5. As seen in Figs. 4 and 5, the temperature as function of the fractional parameter, has different behavior for small, respectively for large values of the time t. This behavior is in accordance with the previous results presented in Figs. 2 and 3 and with the evolution of the thermal memory kernel hβ (t). We must note that the layers temperatures have small variations in time and after a low value of the time t (t = 1 in analyzed cases), the temperatures are almost constant. Also, from numerical simulations, was observed that the influence of the second/third thermal fractional parameter on the fluid temperature is smaller than the influence of the first fractional parameter. This behavior is due to the thermal boundary conditions because, only the wall y = 0 is heated, so the interfacial temperatures decreases away from this wall. To investigate the velocity of fluids of the three layers, we considered the particular case when the velocities of channel walls are constant, namely, u1 (0, t) = f (t) = 0.5H(t), u3 (1, t) = g(t) = 0.4H(t), H(t) being the unit step Heaviside function and the pressure gradient is P (t) = 1 + sin(t). Figs. 7-9 were plotted in order to study the influence of the velocity fractional parameters αi on the velocity fields. The fractional parameters have effects of braking. The flows of the first two layers are slowed down but, the fluid in the third layer is accelerated. Even if, the wall y = 0 moves with the constant velocity, fluid velocity decreases due viscosities and the tα b α−1 damping of the velocity gradient with the kernel χα (t) = λ t Eα,α − λ . The profiles of the function χα (t) are drawn in Fig. 13. For t = 1, the values of stress kernel decrease with α ≥ 0.5, therefore the memory effects are weakened. The influence of the viscosities ratios bi = µµ1i on the fluid velocity is analyzed in Figs. 1012. It is observed that, for decreasing values of the ratio bi , fluid velocity increases. Indeed, in this case the viscosities of fluids in second and third layer are smaller than viscosity of the 21
first fluid; therefore after the first interface the viscous forces are smaller than the viscous force in the first layer. For increasing values of the ratio bi , fluid flows slower. In order to verify the accuracy of numerical results, we used the analytical solution of velocity given by Eq. (72) and the semi-analytical solution given by Eq. (75) for following values of parameters: d0 = 0.3, d1 = 0.7, d2 = 1, t = 2.5, α1 = 0.8, α2 = 0.7, α3 = 0.4, µ1 = 0.6, µ2 = 0.3, µ3 = 0.6 Numerical results are given in Table 1 and Fig. 14. There is a good agreement between the results obtained with the two formulas.
5
Conclusions
Unsteady multi-layer flows with heat transfer of Maxwell immiscible fluids with generalized constitutive equations for the shear stress and thermal flux in a rectangular channel have been investigated. The generalized constitutive equations are described by the Caputo time-fractional derivative, therefore the fluids behavior is influenced by the histories of the velocity gradient and of temperature gradient. Analytical and semi-analytical solutions for the temperature, velocity and shear stress fields are obtained by employing the Laplace and finite sine-Fourier transforms and hypothesis that interfacial velocities, shear stresses, thermal fluxes and temperatures are equal for adjacent fluids. The obtained results have a general character; therefore many particular cases can be obtained. Multi-layer flows of ordinary/fractional Maxwell fluids are studied as particular cases by letting some fractional parameters equal with one. Numerical values of the temperature and velocity fields are obtained by using the Mathcad software and the numerical Talbots algorithms for the inverse Laplace transforms. Numerical results are used for graphical illustrations of the fluid temperature and velocity. The heat transfer and fluid motions are different for small, respectively large values of the time t. These different behaviors are due to the thermal/ velocity kernels variations with respect to time and fractional parameters, therefore the memory effects have significant influence on the fluids. This behavior is in agreement with the variation of the thermal kernel hβ (t) (See Fig. 6). For the values of the thermal fractional parameter β ≤ 0.35 and small values of the time t, the thermal kernel hβ (t) increases with β therefore, the values of the temperature gradient are strongly damped and the heat transfer is reduced. For the thermal fractional parameter β > 0.35, the values of thermal kernel decreases, the damping of the thermal gradient is more weakly and the heat transfer increases. The stress damping kernel χα (t) has the variation similar to the thermal damping kernel (See Fig. 13), therefore, the effect on the fluid velocity is stronger for small values of the stress fractional parameter α.
References [1] D. K. Satpathi, B. V. Rathish, P. Chandra, Unsteady-state laminar flow of viscoelastic gel and air in a channel : Application to mucus transport in a cough machine simulating trachea, math. Comput. Modelling, 38 (2003) 63-75, DOI: 10.1016/S0893-7177(03)00192-4.
22
Figure 2: Temperature profile Ti (y, t) versus y for small values of time and different values of fractional parameters.
Table 1: Comparison between results obtained with analytical and semi-analytical solutions
23
Figure 3: Temperature profile Ti (y, t) versus y for large values of time and different values of fractional parameters.
24
Figure 4: Time-variation of temperature profile Ti (y, t) with the fractional parameter β1 .
25
Figure 5: Time-variation of temperature profile Ti (y, t) with the fractional parameter β2 .
26
Figure 6: Variation of the kernel hβ (t) = tβ−1 /Γ(β) with the fractional parameter β.
Figure 7: Profiles of velocities ui (y, t) for the fractional parameter α1 variable.
27
Figure 8: Profiles of velocities ui (y, t) for the fractional parameter α2 variable.
Figure 9: Profiles of velocities ui (y, t) for the fractional parameter α3 variable.
28
Figure 10: Profiles of velocities ui (y, t) for µ1 variable.
Figure 11: Profiles of velocities ui (y, t) for µ2 variable.
29
Figure 12: Profiles of velocities ui (y, t) for µ3 variable.
Figure 13: Stress kernel χα (t).
30
Figure 14: Comparison between results obtained with analytical and semi-analytical solutions. [2] D. D. Joseph, Y. Y. Renardy, Fundamentals of two-fluid dynamics, Part I: Mathematical Theory and Applications, Springer-Verlag, 1992. [3] Y. Y. Su, B. Khomami, Interfacial stability of multilayer viscoelastic fluids in slit and converging channel dies geometries, J. Rheol., 36 (1992) 357-387. [4] N. R. Anturkar, T. C. Papanastasiou, J. O. Wilkes, Stability of multilayer extrusion of viscoelastic liquids, AIChE Journal, 36 (5) (1990), 710-724. [5] A. C. B. Bogaerds, M. A. Hulsen, G. W. M. Peters, F. P. T. Baaijens, Time dependent finite element analysis of the linear stability of viscoelastic flows with interfaces, J. Non-Newtonian Fluid Mech. 116 (2003), 33-54. [6] H. Le Meur, Non-uniqueness and linear stability of the one-dimensional flow of multiple viscoelastic fluids, Math. Modell. Numer. Anal. 31 (2) (1997) 185-212. [7] P. A. Vasquez, Y. Jin, E. Palmer, D. Hill, M. G. Forest, Modeling and simulation of mucus flow in human bronchial epithelial cell cultures-Part I: Idealized Axysymmetric swirling flow, PLoS Comput. Biol. 12 (8) (2016) e1004872, doi: 10.1371/journal.pcbi.1004872. [8] L. L. Barannyk, D. T. Papageorgiou, P. G. Petropoulos, J. M. Vanden-Broeck, Nonlinear dynamics and wall touch-up in unstability stratified multilayer flows in horizontal channel, SIAM J. Appl. Math. 75 (1) (2015) 92-113, doi: 10.1137/140968070. [9] E. S. Papaefthymiou, D. T. Papageorgiou, Nonlinear stability in three-layer channel flows, J. Fluid Mech., 829 R2 (2017), doi: 10. 1017/jfm.2017.605. [10] J. Zhao, L. Zheng, X. Zhang, F. Liu, Unsteady boundary layer natural convection heat transfer of fractional Maxwell viscoelastic fluid over a vertical plate, Int. J. Heat Mass Transfer, 97 (2016)760-766. [11] X. Chen, W. Yang, X. Zhang, F. Liu, Unsteady boundary layer flow of viscoelastic MHD fluid with a double fractional Maxwell model, Appl. Math. Letters, 95 (2019) 143-149. [12] W. Chen, J. Zhang, J, Zhang, A variable-order time-fractional derivative model for chloride ions subdiffusion in concrete structures, Fractional Calculus and Appl. Ananl., 16(1) (2013) 7692.
31
[13] X. H. Ying, J. X. Yun, Time fractional dual-phase lag heat conduction equation, Chin. Phys. B, 24 (3), (2015) 034401. [14] Muhammad Danial Hisham , Abdul Rauf , Dumitru Vieru , Aziz Ullah Awan , Analytical and semianalytical solutions to flows of two immiscible Maxwell fluids between moving plates, Chinese Journal of Physics (2018), doi: https://doi.org/10.1016/j.cjph.2018.10.009 [15] I. Podlubny, Fractional Differential Equations, Vol. 198 of Mathematics in Science and Engineering, Academic Press, San Diego, USA, 1999. [16] M. Caputo, Linear model of dissipation whose Q is almost frequency independent II, The Geophys. J. R. A. Soc., 13 (1967) 529-539. [17] M. Caputo, M. Fbrizio, A new definition of fractional derivative without singular kernel, Progr. Fract. Differ. Appl., 1 (2) (2015) 73-85. [18] 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. [19] T. R. Prabhakar, A singular integral equation with a generalized Mittag-Leffler function in the kernel, Yokohama Math. J., 19 (1971) 7-15. [20] J. Hristov, Transient space-fractional diffusion with a power-law superdiffusivity: Approximate integralbalance approach, Fundamenta Informaticae, 151 (2017) 371-388, doi: 10.3233/FI-2017-1498. [21] 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. [22] Y. Jiang, H. Qi, H. Xu, X. Jiang, Transient electroosmotic slip flow of fractional Oldroyd-B fluids, Microfluid nanofluid, 21 : 7 (2017) doi : 10.1007/s10404-016-1843-x. [23] 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. [24] X. Yang, H. Qi, X. Jiang, Numerical analysis for electroosmotic flow of fractional Maxwell fluid, Appl. Math. Lett., 78 (2018) 1- 8. [25] Y. Z. Povstenko, Fractional heat conduction equation and associated thermal stress, J. Thermal Stresses, 28 (1) (2004) 83-102. [26] J. Hristov, Integral-balance approach to 1-D space-fractional diffusion models, Mathematical Methods in Engineering Applications in Dynamics of Complex Systems, Springer, 2018. [27] J. Hristov, Integral-Balance Solution to Nonlinear Subdiffusion Equation, Ch. 3, in Frontiers in Fractional Calculus, 1 (2017), 71-106. [28] 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 Kenan Tas, Dumitru Baleanu and J. A. T. Machado, Springer, 2018. [29] D. Brian, Integral Transforms and Their Applications (third ed.) New York: Springer, 2002. [30] R. N. Bracewell, The Fourier Transform and its Applications (Third ed.) Boston: Mc Graw-Hill, 2000. [31] J. Abate, P. P. Valko, Multi-precision Laplace transform inversion, Int. J. Numer. Meth. Engng., 60 (2004) 979-993, doi:10.1002/nme.995. [32] B. Dingfelder, J. A. C. Weideman, An improved Talbot method for numerical Laplace transform inversion, Numer. Algor., 68 (2015) 167-183, doi: 10.1007/s11075-014-9895-z.
32
[33] C. F. Lorenzo, T. T. Hartley, Generalized functions for the fractional calculus, NASA/TP-1999209424/REV1. [34] 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.
33
Journal Pre-proof
Multi-layer flows of immiscible fractional Maxwell fluids with generalized thermal flux. Abdul Rauf, Yasir Mahsud, Itrat Abbas Mirza, Qammar Rubbab PII: DOI: Reference:
S0577-9073(19)30946-3 https://doi.org/10.1016/j.cjph.2019.10.006 CJPH 969
To appear in:
Chinese Journal of Physics
Received date: Revised date: Accepted date:
20 December 2018 14 June 2019 9 October 2019
Please cite this article as: Abdul Rauf, Yasir Mahsud, Itrat Abbas Mirza, Qammar Rubbab, Multi-layer flows of immiscible fractional Maxwell fluids with generalized thermal flux., Chinese Journal of Physics (2019), doi: https://doi.org/10.1016/j.cjph.2019.10.006
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 with damped shear stresses are studied. • Memory effects are considered in the thermal transport process. • For large values of the time t, the fluid temperature is smaller in fractional fluids. • Fluids characterized by the fractional constitutive equations flow slower than the ordinary fluids.
1
Multi-layer flows of immiscible fractional Maxwell fluids with generalized thermal flux.
Abdul Rauf 1 ∗ Yasir Mahsud 2 Itrat Abbas Mirza 3 Qammar Rubbab 4 1 Dept. of Computer Science and Engineering, Air University Multan, Pakistan 2 Abdus Salam School of Mathematical Sciences, GCU Lahore, Pakistan 3 Dept. of Mathematics, Khwaja Fareed University of Engineering & information Technology, Rahim Yar Khan, Pakistan 4 Dept. of Mathematics, The Woman University Multan, Pakistan
Abstract Unsteady laminar flows and heat transfer of n-immiscible fractional Maxwell fluids in a channel are investigated under influence of time-dependent pressure gradient. The isothermal channel walls have translational motions in their planes with time-dependent velocities. Governing equations of the mathematical model are based on the generalized constitutive equations for shear stress and thermal flux described by the time-fractional Caputo derivative. Analytical and semi-analytical solutions for velocity, shear stress, and temperature fields are obtained by using finite sine-Fourier and Laplace transforms. In the case of semi-analytical solutions, the inverse Laplace transforms are obtained numerically by employing the Talbots algorithms. Using the software Mathcad, numerical calculations have carried out and results are presented in graphical illustrations in order to analyze the memory effects on the fluid temperature and motion. It is found that in fluids with thermal memory the heat transfer is slower compared with the ordinary fluid, while the fractional velocity parameters act as braking/accelerating factors of the fluids. Keywords: n-layered immiscible fluids; generalized thermal flux; fractional Maxwell fluids; analytical and semi-analytical solutions
1
Introduction
Flows of layered immiscible fluids in a pipe/channel are important for many industrial and physiological applications. For example, the study of two-layer flows of immiscible fluids ∗
Address correspondence to Abdul Rauf. E-mail: attari [email protected]; abdul.rauf @aumc.edu.pk
2
could be useful to the optimal transport of crude oils, polymer solutions or mucus transport in the respiratory tract [1]. To obtain polymer films or plaques, the co-extrusion process of different polymers is used. Sometimes, due to technological process imperfections, the finished product may have some wavy interfaces (interfacial instability). Theoretical studies of the stability of multi-layer plane Poiseuille flows give important information about parameters that control the stability of multi-layer flows such as layer thickness ratios (flow-rate ratios), viscosity ratios or, density ratios [2, 3]. Anturkar et al. [4] have investigated the linear interfacial stability of multilayer plane Poiseuille flows of Oldroyd-B fluids with shear rate dependent viscosity. They determined asymptotic solutions at long wavelengths and found that multilayer flows of viscoelastic fluids can be stable for all wavelengths. Also, it is obtained that, symmetric and nearly symmetric configurations with three layers are unstable if the inner fluid is more viscous than the outer fluids [4]. Linear stability of flows of viscoelastic fluids with interfaces has been analyzed by Bogaerds et al. [5] 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 [6]. In this interesting work, the existence and uniqueness for the multilayer Couette/ Poiseuille flows in channels/pipes 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 solution. Vasquez et al. [7] have modeled and simulated the mucus flow in human bronchial epithelial cell cultures. They studied idealized axisymmetric swirling flows considering four kinds of fluids Newtonian fluid, upper-convected Maxwell fluid, Giesekus fluid and 5-mode Giesekus fluid. Barannyk et al [8] have studied nonlinear dynamic of stratified immiscible, perfect dielectric fluids in horizontal channels in the presence of an electric field perpendicular on the gravitational field. They found that a strong electric field could stabilize the Rayleigh -Taylor interfacial instability. Flows of three-layer of different fluids driven by the gravity and a pressure gradient in a slightly inclined channel have been investigated by Papaefthymiou and Papageorgiou [9]. They observed that nonlinear stability depends on the existence of bounded invariant regions in the hyperbolic state space of the interfacial. Zhao et al. [10] have studied unsteady boundary layer natural convection heat transfer of the fractional Maxwell fluids. Mainly, they found that these coefficients are several times larger than those in non-swirling flows. The interfacial coefficient is well correlated in terms of the liquid volume fraction and the gas Reynolds number, while both gas Reynolds number and liquid Reynolds number are required for correlating the wall friction coefficient. Other interesting problems can be found in [11–14]. It is worth point out that many complex systems in chemistry, mechanics or biology are modeled 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 are described more accurate by the fractional differential equations. There are many fractional differential operators studied by the fractional calculus. Riemann-Liouville fractional integral/derivative [15], Caputo fractional derivative [16], Caputo-Fabrizio fractional derivative [17], AtagnanaBaleanu fractional derivative [18] or Prabhakar fractional derivative [19] are only some from the fractional differential operators used in viscoelasticty, thermal and mass transport. Hristov [20] has applied the Riemann-Liouville fractional derivative to study the transient space-fractional diffusion with power-law superdiffusivity. Mahsud et al. [21] have 3
investigated boundary layer flows of an upper-convected Maxwell fluid modeled by the timefractional Caputo-Fabrizio derivatives. They found that fractional fluids flow slower than the ordinary fluids. Using Laplace and Hankel transforms, Jiang et al. [22] have studied electroosmotic flows of Oldroyd-B fluids with a fractional constitutive equation based on the time-fractional Riemann-Liouville derivative. Wang et al. [23] have investigated the unsteady electroosmotic flows of generalized Maxwell fluids in cylindrical domains and made comparison between ordinary and generalized Maxwell fluids. Ying et al. [24] have studied the electroosmotic flows of the fractional Maxwell fluids in a rectangular microchannel using a fully discrete spectral method for the temporal variable and a Legendre spectral method in the spatial direction. Theoretical results and practical applications of the fractional calculus can be found in the references [25–28]. In this paper we investigate multi-layer flows with heat transfer of Maxwell fluids with generalized constitutive equations for the shear stress and thermal flux. The flow domain is a rectangular channel bounded of two parallel isothermal plates which have translational motion in their planes with time-dependent velocities. The generalized constitutive equations, used in this work, are described by the Caputo time-fractional derivative, therefore the histories of the velocity gradient and histories of temperature gradient influence the fluid shear stress and temperature. Analytical solutions for the Laplace transforms of temperature, velocity and shear stress fields are obtained by employing the Laplace and finite sine-Fourier transforms and hypothesis that interfacial velocities, shear stresses, thermal fluxes and temperatures are equal for adjacent fluids. The obtained results have a general character; therefore many particular cases can be obtained. Multi-layer flows of ordinary/fractional Maxwell fluids are studied as particular cases by letting some fractional parameters equal with one. Numerical values of the temperature and velocity fields are obtained by using the Mathcad software and the numerical Talbots algorithms for the inverse Laplace transforms. Numerical results are used for graphical illustrations of the fluid temperature and velocity. According with the generalized constitutive equations for the shear stress and thermal flux, the variations of the thermal/velocity kernels with the fractional parameter are different for small, respectively large time. It is found that at small values of time t, the heat transfer is smaller in fluids with generalized thermal flux than in the ordinary. At large values of time t, the behavior is opposite. The velocity memory fractional parameters act as braking/accelerating factors of the fluids.
2
Mathematical Modeling
The flow domain is D0 = {(x0 , y 0 , z 0 )| x0 ∈ (−∞, ∞), y 0 ∈ [0, h], z 0 ∈ (−∞, ∞)}, with the boundary walls situated in planes y 0 = 0 and y 0 = h > 0 . Initially at time t0 = 0, the two plates and the enclosed fluids are at rest with the ambient temperature is T2o . After this moment, the channel flat wall situated in the plane y 0 = 0 is moving along the x0 -axis with the velocity u01o = U0 f1 (t0 ) and the temperature T1o , while the plane wall y 0 = h translates parallel with the x0 -axis with the velocity u02o = U0 g1 (t0 ) and the temperature T2o (Fig.1). Functions f1 (t0 ), g1 (t0 ) are differentiable functions of exponential order at infinity and f1 (0) = g1 (0) = 0. In this paper the velocity field is considered to be of the form 4
Figure 1: Geometry of the problem. V~ 0 = (u0 (y 0 , t0 ), 0, 0). We consider here n-layers Maxwell fluids between two parallel plates. Let h0 = 0 and hn = h. In the region y 0 ∈ [hi−1 , hi ], hi−1 < hi flows a Maxwell fluid with the density ρi , viscosity µi , relaxation time λ0i = Gµii , Gi the elastic modulus, velocity u0i (y 0 , t0 ), temperature Ti0 (y 0 , t0 ), thermal flux qi0 (y 0 , t0 ) and the shear stress τi0 (y 0 , t0 ), where i = 1, 2, . . . , n. The fluids are assumed to be immiscible and incompressible and the flow is one-dimensional, unsteady and fully developed. The flows of the fluids are generated by the time-dependent pressure gradient in the flow direction and by the movements of the two walls of the channel. Based on the assumptions about velocity fields, the continuity equation is identically 0 := satisfied by all velocities ui (y 0 , t0 ), i = 1, 2, · · · , n. We fix here In1 := {1, 2, · · · , n}, In−1 {0, 1, · · · , n − 1}. The governing equations of motion along with initial, boundary and the interface fluid-fluid conditions are • the system of linear momentum equations ρi
∂u0 i ∂τ 0 i ∂p0 = − 0 , i ∈ In1 , ∂t0 ∂y 0 ∂x
(1)
• the system of constitutive equations τi0 + λ0i
∂τ 0 i ∂u0 i = µ , i ∈ In1 , i ∂t0 ∂y 0
(2)
• the initial conditions u0i (y 0 , 0) = 0, τi0 (y 0 , 0) = 0, i ∈ In1 , 5
(3)
• the boundary condition on the bottom plate y 0 = 0, u01 (0, t0 ) = U0 f1 (t0 ),
(4)
on the upper plate y 0 = hn = h, u0n (h, t0 ) = U0 g1 (t0 ), • the interfaces conditions 0 1 u0i (hi , t0 ) = u0i+1 (hi , t0 ), τi0 (hi , t0 ) = τi+1 (hi , t0 ), i ∈ In−1 .
(5)
The system of governing equations for thermal transport and heat flux are given by • the balance equations of thermal transport ρi cip
∂Ti0 (y 0 , t0 ) ∂qi0 (y 0 , t0 ) = − , ∂t0 ∂y 0
• the Fourier’s law of thermal flux qi0 (x0 , t0 ) = −ki
(6)
∂Ti0 (y 0 , t0 ) , ∂y 0
(7)
where i ∈ In1 , cip and ki are the specific heat at constant pressure and thermal conductivity, respectively. To the thermal transport we consider the following initial, boundary and the interface fluid-fluid conditions: • initial conditions Ti0 (y 0 , 0) = T2o , qi0 (y 0 , 0) = q0 , i ∈ In1 ,
(8)
T 0 (0, t0 ) = T1o , Tn0 (h, t) = T2o ,
(9)
• boundary conditions
• interface conditions 0 0 1 Ti0 (hi , t0 ) = Ti+1 (hi , t0 ), qi0 (hi , t0 ) = qi+1 (hi , t0 ), i ∈ In−1 .
(10)
Introducing the non-dimensional variables 0 0 0 0 0 x0 , y = yh , t = νh12t , ui = uU0i , τi = µhτ1 Ui0 , νi = µρii , p = µhp , h U 1 0 2 2 0 λi = λhi2ν1 , ai = ρρ1i , bi = µµ1i , f (t) = f1 hν1t , g(t) = g1 hν1t , di Ti0 −T2o q0 µc Ti = T1o−T , qi = q0i , σi = νν1i , P ri = ikiip , i ∈ In1 , 2o
x=
=
hi , h
into Eqs. (1)-(10), the mathematical model for dimensionless governing equations takes the form, ai
∂ui ∂p ∂τi =− + , i ∈ In1 , ∂t ∂x ∂y 6
(11)
τi + λi
∂τi ∂ui = bi , i ∈ In1 , ∂t ∂y
ν1 γ0 ρi cip
∂Ti ∂qi = − , i ∈ In1 , ∂t ∂y
qi = −γ0 ki
∂Ti , i ∈ In1 , ∂y
(12)
(13)
(14)
along with • the non-dimensional initial conditions ui (y, 0) = 0, τi (y, 0) = 0, Ti (y, 0) = 0, qi (y, 0) = 1, i ∈ In1 ,
(15)
• the non-dimensional boundary condition u1 (0, t) = f (t), un (1, t) = g(t), T1 (0, t) = 1, Tn (1, t) = 0,
(16)
• the non-dimensional interface conditions ui (di , t) = ui+1 (di , t), τi (di , t) = τi+1 (di , t),
(17)
1 Ti (di , t) = Ti+1 (di , t), qi (di , t) = qi+1 (di , t), i ∈ In−1 ,
where γ0 =
2.1
T1o −T2o . q0 h
Mathematical model with fractional constitutive equations
In the following , we consider the generalized constitutive equations, τi + λi Dαt i τi = bi
qi =
−γ0 ki Dt1−βi
Dγt F (y, τ )
1 = Γ(1 − γ)
∂ui , αi ∈ (0, 1], i ∈ In1 , ∂y
(18)
∂Ti , βi ∈ (0, 1], i ∈ In1 , ∂y
(19)
where Z
0
t
(t − τ )−γ
∂F (y, τ ) dτ, 0 ≤ γ < 1 ∂τ
(20)
(y,t) is the Caputo derivative. For particular case γ = 1, D1t F (y, t) = ∂F∂t . It is important to note that the fractional constitutive equations (18) and (19) can be written in the equivalent forms bi αi −1 tαi ∂ui (y, t) τi (y, t) = t Eαi ,αi − ∗ (21) λi λi ∂y
7
=
Z
0
t
bi (t − τ )αi −1 Eαi ,αi λi
(t − τ )αi − λi
tβi −1 ∂ 2 Ti (y, t) qi (y, t) = ∗ = Γ(βi ) ∂t∂y
Z
t
0
∂ui (y, τ ) dτ, i ∈ In1 , ∂y
(t − τ )βi −1 ∂ 2 Ti (y, t) dτ, i ∈ In1 . Γ(βi ) ∂t∂y
(22)
The Eqs. (21) and (22) highlight the fact that the histories of the velocity gradient and temperature gradient influence the time-evaluation of shear stress and thermal flux, respectively. Also it is observed from theaboveequations that, the nonlocality kernel for the shear stress is the function
bi αi −1 t Eαi ,αi λi
α
− tλii , Eαi ,αi (.) being the Mittag-Leffler function; while
the nonlocality kernel of the thermal flux is the power-law function the pressure gradient is a known function of t, therefore −
tβi −1 . Γ(βi )
We suppose that
∂p = P (t), ∂x
(23)
where P (t) is a local integrable function on [0, ∞) of the exponential order at infinity.
3
Solution of the problem
In order to determine analytical solutions of the Eqs. (11), (13, (18), (19) with conditions (15)-(17) we use the Laplace transform coupled with the finite sine-Fourier transform [29]. Applying the Laplace transform to Eqs. (11), (13), (16)-(19) and using the initial conditions (15), we obtain the problem in transformed domain ∂ τ¯i (y, s) , i ∈ In1 , ∂y
(24)
bi ∂ u¯i (y, s) , i ∈ In1 , α 1 + λi s i ∂y
(25)
ai s¯ ui (y, s) = P¯ (s) +
τ¯i (y, s) =
∂ q¯i (y, s) , i ∈ In1 , ν1 γ0 ρi cip sT¯i (y, s) = − ∂y
q¯i = −γ0 ki s1−βi
∂ T¯i (y, s) , i ∈ In1 , ∂y
(26)
(27)
1 u¯1 (0, s) = f¯(s), u¯n (1, s) = g¯(s), T¯1 (0, s) = , T¯n (1, s) = 0, s
(28)
u¯j (dj , s) = u¯j+1 (dj , s), τ¯j (dj , s) = τ¯j+1 (dj , s),
(29)
1 T¯j (dj , s) = T¯j+1 (dj , s), q¯j (dj , s) = q¯j+1 (dj , s), j ∈ In−1 .
8
3.1
Solution for the thermal transport
Using Eqs. (26) and (27) we can write the Laplace transformed system of equations for thermal transports T¯i (y, s), i ∈ In1 , in the form ∂ 2 T¯i = ζi sβi T¯i , i ∈ In1 , ∂y 2 where ζi =
P ri . σi
(30)
The general solution of Eq. (30) is T¯i (y, s) = Ai (s)e−y
√
ζ i sβ i
+ Bi (s)ey
√
ζ i sβ i
, i ∈ In1 ,
(31)
where the unknown parameters A1 , A2 . . . , An , B1 , B2 , . . . , Bn are to determined by applying the boundary conditions (28) and interface conditions (29) to Eq. (31) and are given by the following linear system M(s)X(s) = C(s), where the matrix M(s) is given by 1 0 0 ··· 0 a −b 0 ··· 0 1 1 0 a −b · · · 0 2 2 .. . .. . . .. .. .. . . 0 0 · · · an−1 −bn−1 M(s) = p −q 0 ··· 0 1 1 0 p −q · · · 0 2 2 .. . .. ... ... .. . . 0 0 · · · pn−1 −qn−1 0 0 ··· 0 pn
and the columns X(s) and C(s) are
Here
−dj
aj (s) = e
√
−2
pn (s) = e `j (s) =
kj+1 ,j kj
X(s) =
√
ζ j sβ j
ζn sβn
An B1 B2 .. . Bn
−dj
, bj (s) = e
, qj (s) = −`j s
1 ∈ In−1 .
A1 A2 .. .
√
1 1 a1
0 .. .
(32)
0 − b11
βj −βj+1
1 a2
··· ··· ··· .. .
1 s
0 0 .. . 0
,
.
√ β p ζ s j β j , pj (s) = − ζj s aj (s), rj (s) = aj j(s) ,
p bj (s) ζj+1 sβj+1 , sj (s) = 9
0 0 0 .. .
.. . 1 1 − bn−1 0 0 · · · an−1 r1 −s1 0 ··· 0 0 r2 −s2 · · · 0 .. .. .. .. .. . . . . . 0 0 · · · rn−1 −sn−1 0 0 ··· 0 1
, C(s) =
ζj+1 sβj+1
0 0 − b12 .. .
√
`j sβj −βj+1 ζj+1 sβj+1 , bj (s)
In order to write the system (32) in a suitable form, we introduce the notations
M =
N =
Q=
R=
1 0 0 ··· 0 a1 −b1 0 ··· 0 0 a2 −b2 · · · 0 .. .. .. ... ... . . . 0 0 · · · an−1 −bn−1
1
0 − b11
1 a1
1 a2
0 .. .
.. . 0
0
0 0 − b12 .. .
··· ··· ··· .. .
···
1
− bn−1
r1 −s1 0 ··· 0 0 r2 −s2 · · · 0 .. .. .. ... ... . . . 0 0 · · · rn−1 −sn−1 0 0 ··· 0 1
C1 =
1 s
0 .. . 0
n×1
, C2 =
1
p1 −q1 0 ··· 0 0 p2 −q2 · · · 0 .. .. .. ... ... . . . 0 0 · · · pn−1 −qn−1 0 0 ··· 0 pn
0 0 .. . 0
; i.e.
Mij = ai δi,j+1 − bi δi+1,j+1 , n×n
0 1 , ,j ∈ In−1 i ∈ In−1
0 , N0j = δ0,j , j ∈ In−1
; i.e.
Nij = n×n
1 δ ai i,j+1
−
(34)
1 δ , bi i+1,j+1
1 0 i ∈ In−1 ,j ∈ In−1 ,
Q = (Qij )i,j∈I 0 , n−1
; i.e.
Qij = pi+1 δi+1,j+1 − qi+1 δi+2,j+1 , i∈
n×n
0 ,j In−2
∈
(35)
0 , In−1
0 Q(n−1)j = pn δn−1,j , j ∈ In−1 ,
R = (Rij )i,j∈I 0 ,
n−1
; i.e.
Rij = ri+1 δi+1,j+1 − si+1 δi+2,j+1 , i∈
n×n
0 In−2 ,j
,A =
A1 A2 .. . An
n×1
,B =
where δij is the Kronecker tensor. The system (32) is written as M A + N B = C1 , QA + RB = C2 . 10
∈
(36)
0 In−1 ,
0 R(n−1)j = δn−1,j , j ∈ In−1 ,
n×1
(33)
n−1
0 , M0j = δ0,j , j ∈ In−1
N = (Nij )i,j∈I 0 ,
n−1
0 0 0 .. .
an−1
M = (Mij )i,j∈I 0 ,
B1 B2 .. . Bn
,
(37)
n×1
(38)
Matrices M, N, Q, R are triangular and invertible matrices, Eqs. (38) can be written as (N −1 M − R−1 Q)A = N −1 C1 , (39) B = −R−1 QA. Suppose that matrix S = N −1 M − R−1 Q is invertible, we have A = S −1 N −1 C1 , B = −R−1 QS −1 N −1 C1 .
(40)
An easy computation shows that
1 C0 := N −1 C1 = s
where,
N −1
that is,
=
1
1 b1 a1 b1 b2 a1 a2
.. .
b1 b2 ···bn−1 a1 a2 ···an−1
b1 a1 b1 b2 a1 a2
.. .
0 −b1 − b1a2b2 .. .
0 0 −b2 .. .
b1 b2 ···bn−1 a1 a2 ···an−1
···bn−1 − b1a2b2···a n−1
···
··· ··· ··· .. .
(41)
0 0 0 .. .
bn−1 − bn−2 −bn−1 an−1
,
n×n
0 N −1 = (nij )i,j∈In−1 ,
0 , ni0 = n0j = δ0,j , j ∈ In−1
nij = −
i Q
m=1
bm am
This reduces the system (40) to
j−1 Q `=0
a`+1 b`
n−1 P k=0
i Q
`=1
b` ,i a`
1 , ∈ In−1
(42)
1 δi,j+k , b0 = 1, i, j ∈ In−1 .
A = S −1 C0 , B = −R−1 QS −1 C0 .
(43)
The inverse matrices of M and R are 0 M −1 = (mij )i,j∈In−1
0 m0j = δ0,j , j ∈ In−1 , mi0 =
mij = −
i Q
m=1
am bm
j−1 Q `=0
b` a`+1
n−1 P k=0
i Q
`=1
a` ,i b`
1 ∈ In−1 ,
1 δi,j+k , b0 = 1, i, j ∈ In−1 ;
11
(44)
0 R−1 = (ρij )i,j∈In−1 ,
ρij =
n−1 P k=0
i Q
m=0
rm sm
j Q
`=0
s` r`+1
(45) δi−k,j , r0 = s0 = rn = 1, i, j ∈
0 In−1 .
0 Now, the matrix S = (Sij )i,j∈In−1 is defined by the elements
Sij =
n−1 X k=0
(nik Mkj − ρik Qkj ).
(46)
After having determined the parameters Aj , Bj , j ∈ In1 , from the system (43) the analytic solution for the system of thermal transport Ti (y, t) can be determined by applying the inverse Laplce transform to Eq. (31) 1 Ti (y, t) = 2πι
Z
σ+i∞
est (Ai (s)e−y
√
ζi sβi
+ Bi (s)ey
√
σ−i∞
ζi sβi
)ds, i ∈ In1 .
(47)
In this paper, the inverse Laplace transforms are obtained by using Talbot’s algorithms [31, 32]. ¯ s) be the Laplace transform of the function g(y, t). The function g(y, t) can be Let G(y, approximated with the help of Talbot algorithm [31] for the Laplace transform inversion, ( ) M −1 X r exp(rt) ¯ r) + ¯ z(θk )) (1 + iζ(θk )) , g(y, t) ∼ G(y, Re exp(tz(θk ))G(y, (48) = M 2 k=1 where r=
2M , 5t
z(θ) = rθ(cot θ + i), θ ∈ (−π, π),
ζ(θ) = θ + (θcotθ − 1)cotθ, θk =
(49)
kπ . M
Another method to approximate the function g(y, t) is the improved Talbot algorithm for the inverse Laplace transform and is given by [32] M
g(y, t) ∼ =
1X ¯ z1 (σk )) (ν + iζ 1 (σk )), exp(tz1 (σk ))G(y, t k=1
(50)
where, z1 (θ) =
M t
[νiθ + µθ cot(αθ) − ξ] , θ ∈ [−π, π],
ζ1 (θ) = αµθ + µ(αθcot(αθ) − 1)cot(αθ), σk =
(2k−1)π M
Here α, M, ν, µ, ξ are parameters to be specified by the user.
12
(51) − π.
3.2
Analytical solutions for velocity and shear stress
Multiplying Eq. (24) by 1 + λi sαi both sides and using Eq. (25) we obtained the system of equations for Laplace transformed velocities, ∂ 2 u¯i (y, s) , i ∈ In1 , ai s(1 + λi sαi )¯ ui (y, s) = (1 + λi sαi )P¯ (s) + bi ∂y 2 ¯ s) = where φ(y,
R∞
(52)
φ(y, t) exp(−st)dt is the Laplace transform of the function φ(y, t) [29].
0
In order to find the exact solution of velocities ui (y, t) we will use the finite sine-Fourier transform along with the Laplace transform. The finite sine-Fourier transform of the function ¯ s), y ∈ [a, b] is defined as [30] φ(y, e (s) = φ m
Zb a
¯ s) sin(ϑm (y − a))dy, ϑm = mπ , b > a, m = 1, 2, · · · , φ(y, b−a
(53)
along with the inverse Fourier transform defined by ¯ s) = φ(y,
∞ 2 X e (s). sin(ϑm (y − a))φ m b − a m=1
(54)
With the application of finite Fourier sine transform (53) to Eq. (52) along with the boundary conditions (28) and (29), the transformed velocities take the form eim (s) = u
1−(−1)m
(i)
(1+λi sαi )P¯ (s)
(i) ϑm
ϑm =
sαi )+b
ai s(1+λi mπ ,m = di −di−1
(i) 2 i (ϑm )
+
1, 2, · · · , i
(i)
bi ϑm u ¯i (di−1 ,s) ai s(1+λi ∈ In1 ,
sαi )+b
(i) 2 i (ϑm )
−
(i)
(−1)m bi ϑm u ¯i (di ,s) (i) ai s(1+λi sαi )+bi (ϑm )2
,
(55)
where for i = 1, u¯1 (d0 , s) = f¯(s) and for i = n, u¯n (dn , s) = g¯(s). In order to apply the inverse Fourier sine transform, we rewrite the Eq. (55) in the following suitable form (−1)m+1 u¯i (di , s) u¯i (di−1 , s) 1 − (−1)m (1 + λi sαi )P¯ (s) e − uim (s) = + + (i) (i) (i) (i) ϑm ϑm ϑm ai s(1 + λi sαi ) + bi (ϑm )2 ai s(1 + λi sαi )¯ ui (di−1 , s) (i)
(i)
ϑ(i) m
(−1)m ai s(1 + λi sαi )¯ ui (di , s)
, (i) (i) ϑm (ai s(1 + λi sαi ) + bi (ϑm )2 ) mπ = , m = 1, 2, · · · , i ∈ In1 . di − di−1
ϑm (ai s(1 + λi sαi ) + bi (ϑm )2 )
+
(56)
Using the following auxiliary functions along with their inverse Fourier sine transform ϕ1i (y) =
di −y , di −di−1 y−di−1 , di −di−1
y ∈ [di−1 , di ], ϕ˜1im =
ϕ2i (y) = y ∈ [di−1 , di ], ϕ˜2im = m = 1, 2, · · · , i ∈ In1 , d0 = 0, dn = 1,
13
1 (i) , ϑm (−1)m+1 (i)
ϑm
,
(57)
the inverse Fourier sine transform of Eq. (56) takes the form i−1 i −y u¯i (y, s) = did−d u¯i (di−1 , s) + dy−d u¯i (di , s)+ i−1 i −di−1 ∞ m P (i) (1+λi sαi )P¯ (s) 2 sin(ϑm (y − di−1 )) 1−(−1) (i) 2 − (i) αi di −di−1
2 di −di−1 2 di −di−1
m=1 ∞ P m=1 ∞ P
ϑm
(i)
sin(ϑm (y − di−1 )) (i) sin(ϑm (y
m=1 (i)
ϑm =
− di−1 ))
mπ ,m di −di−1
ai s(1+λi s
)+bi (ϑm )
ai s(1+λi sαi )¯ ui (di−1 ,s) (i) (i) ϑm (ai s(1+λi sαi )+bi (ϑm )2 ) ui (di ,s) (−1)m ai s(1+λi sαi )¯ (i)
(i)
ϑm (ai s(1+λi sαi )+bi (ϑm )2 )
+
(58)
,
= 1, 2, · · · , i ∈ In1 .
Now, from (25) and (58), we obtain τ¯i (y, s) =
bi ∂ u¯i (y, s) = T¯i1 (y, s)¯ ui (di , s) − T¯i2 (y, s)¯ ui (di−1 , s) + T¯i3 (y, s)P¯ (s), (59) 1 + λi sαi ∂y
where ∞ (i) m αi X b (−1) a (1 + λ s )s cos(ϑ (y − d )) m i i i i−1 T¯i1 (y, s) = 1+2 , αi ) + b (ϑ(i) )2 ) (1 + λi sαi )(di − di−1 ) (a s(1 + λ s m i i i m=1 T¯i2 (y, s) =
∞ (i) X bi ai (1 + λi sαi )s cos(ϑm (y − di−1 )) 1+2 , αi ) + b (ϑ(i) )2 ) (1 + λi sαi )(di − di−1 ) (a s(1 + λ s m i i i m=1
T¯i3 (y, s) =
2bi di − di−1
X ∞
(i)
(1 − (−1)m ) cos(ϑm (y − di−1 ))
m=1
ai s(1 + λi
sαi )
+
(i) bi (ϑm )2
, i ∈ In1 .
(60)
With the help of the interface liquid-liquid conditions (29a), the Laplace transformed velocities on the interfaces y = di , i = 1, 2, . . . , n − 1, can be written as ¯ i1 (s) L 1 , i ∈ In−1 , u¯i (di , s) = u¯i+1 (di , s) = ¯ Li0 (s)
(61)
where, ¯ i0 (s) = T¯i1 (di , s) − T¯(i+1)1 (di , s), L ¯ i1 (s) = (T¯i2 (di , s) − T¯(i+1)2 (di , s)))¯ L ui (di−1 , s) + (T¯(i+1)3 (di , s) − T¯i3 (di , s))P¯ (s). In order to obtained the inverse Laplace transforms of the functions u¯i (y, s), i ∈ In1 and 1 u¯j (dj , s) = u¯j+1 (dj , s), j ∈ In−1 , we consider the following auxiliary functions, ¯ i0 (m, s) = H
=
1 (i)
ai s(1 + λi sαi ) + bi (ϑm )2
s−1 · ai λi (sαi + λ−1 i ) 1+ 14
1 (i)
bi (ϑm )2 s−1 ai λi (sαi +λ−1 i )
∞ (i) 2 k s−k−1 1 X k bi (ϑm ) = (−1) , k+1 ai λi k=0 ai λ i (sαi + λ−1 i )
(62)
¯ i3 (m, s) = sH ¯ i0 (m, s) H ∞ (i) 2 k s−k 1 X k bi (ϑm ) = (−1) , k+1 ai λi k=0 ai λ i (sαi + λ−1 i ) ¯ i1 (m, s) = H
(63)
(1 + λi sαi ) (i)
ai s(1 + λi sαi ) + bi (ϑm )2
∞ (i) 2 k s−k−1 1X k bi (ϑm ) (−1) , = k ai k=0 ai λi (sαi + λ−1 i )
(64)
and ¯ i2 (m, s) := sH ¯ i1 (m, s) H ∞ (i) 2 k 1X 1 s−k k bi (ϑm ) , = + (−1) k ai ai k=1 ai λ i (sαi + λ−1 i )
(65)
Since the generalized G-Lorenzo-Hartley function is defined by [33] −1
Gα,β,γ (t, σ) = L
X ∞ sβ Γ(k + γ)σ k t(k+γ)α−β−1 , = (sα − σ)γ k!Γ((k + γ)α − β)Γ(γ) k=0
Re(s) > 0, Re(αγ − β) > 0, |
(66)
σ | < 1, sα
and for αi , βi > 0, L
−1
sαi −βi = tβi −1 Eαi ,βi (dtαi ), sαi − d
(67)
¯ i0 (m, s), where Eαi ,βi (.) is the Mittag-Leffler function [34]. The inverse Laplace transform of H ¯ i3 (m, s), H ¯ i1 (m, s) and H ¯ i2 (m, s) takes the form H Hi0 (m, t) =
∞ X k=0
(i) 2 k −1 k (bi (ϑm ) ) (−1) Gαi ,−k−1,k+1 (t, −λi ) , (ai λi )k+1
∞ (i) 2 k X −1 k (bi (ϑm ) ) Hi3 (m, t) = (−1) Gαi ,−k,k+1 (t, −λi ) , k+1 (a i λi ) k=0 15
(68)
(69)
∞ (i) 2 k 1X k bi (ϑm ) Gαi ,−k−1,k (t, −λ−1 Hi1 (m, t) = (−1) i ), ai k=0 ai λi
(70)
∞ (i) 2 k δ(t) 1X k bi (ϑm ) Hi2 (m, t) = + Gαi ,−k,k (t, −λ−1 (−1) i ), ai ai k=1 ai λi
(71)
and
where δ(t) is dirac delta function. Using Eqs. (58), (70) and (71), we obtain for velocities ui (y, t), i ∈ In1 , the following expressions: i−1 i −y ui (y, t) = did−d ui (di−1 , t) + dy−d ui (di , t)+ i−1 i −di−1 ∞ m P (i) 2 Hi1 (m, t) ∗ P (t)− sin(ϑm (y − di−1 )) 1−(−1) (i) di −di−1
2 di −di−1 2 di −di−1
m=1 ∞ P m=1 ∞ P
ϑm
(i)
sin(ϑm (y − di−1 ))
ai (i) Hi2 (m, t) ϑm
(i) sin(ϑm (y
(−1)m ai
m=1 (i)
ϑm =
− di−1 ))
mπ ,m di −di−1
(i)
ϑm
∗ ui (di−1 , t)+
(72)
Hi2 (m, t) ∗ ui (di , t),
= 1, 2, · · · , i ∈ In1 .
Rt where h1 (t) ∗ h2 (t) = 0 h1 (t − τ )h2 (τ )dτ is the convolution product of the functions h1 (t), h2 (t) and ui (di , t) = ui+1 (di , t) is given by ui (di , t) = ui+1 (di , t) =
=
Li1 (t) Li0 (t)
(Ti2 (di , t) − T(i+1)2 (di , t)) ∗ ui (di−1 , t) + (T(i+1)3 (di , t) − Ti3 (di , t)) ∗ P (t) , Ti1 (di , t) − T(i+1)1 (di , t)
where, ∞ X bi −1 αi −1 αi −1 Ti1 (di , t) = ai Hi3 (m, t) , λ t Eαi ,αi (−λi t ) + 2 (di − di−1 ) i m=1 ∞ X bi+1 m −1 αi −1 −1 αi T(i+1)1 (y, t) = λ t Eαi ,αi (−λi+1 t ) + 2 (−1) ai+1 H(i+1)3 (m, t) , (di+1 − di ) i+1 m=1 ∞ X bi m −1 αi −1 −1 αi Ti2 (di , t) = λ t Eαi ,αi (−λi t ) + 2 (−1) ai Hi3 (m, t) , (di − di−1 ) i m=1 ∞ X bi+1 −1 αi −1 −1 αi Eαi ,αi (−λi+1 t ) + 2 ai+1 H(i+1)3 (m, t) , T(i+1)2 (di , t) = λ t (di+1 − di ) i+1 m=1 16
(73)
2bi Ti3 (di , t) = di − di−1 2bi+1 T(i+1)3 (di , t) = di+1 − di
X ∞
X ∞
m=1
m=1
((−1) − 1)Hi0 (m, t) , m
1 . (1 − (−1) )H(i+1)0 (m, t) , i ∈ In−1 m
The system of shear stresses τi (y, t), i ∈ In1 , can be determined by applying inverse Laplace transform to Eq. (59) and using Eqs. (68) and (69) τ¯i (y, s) =
bi ∂ u¯i (y, s) = T¯i1 (y, s)¯ ui (di , s) − T¯i2 (y, s)¯ ui (di−1 , s) + T¯i3 (y, s)P¯ (s), (74) α i 1 + λi s ∂y
where ∞ X bi m −1 αi −1 −1 αi (i) λ t Eαi ,αi (−λi t ) + 2 (−1) ai Hi3 (m, t) cos(ϑm (y − di−1 )) , Ti1 (y, t) = di − di−1 i m=1 ∞ X bi −1 αi −1 αi −1 (i) ai Hi3 (m, t) cos(ϑm (y − di−1 )) , λ t Eαi ,αi (−λi t ) + 2 Ti2 (y, t) = di − di−1 i m=1 ∞ X 2bi (i) Ti3 (y, t) = (1 − (−1)m )Hi0 (m, t) cos(ϑm (y − di−1 )), i ∈ In1 . di − di−1 m=1
3.3
Semi-analytical solution for velocity and shear stress
This subsection is devoted to determine a new form of solutions to the problem given by Eqs. (24), (25), (28a), (29a) with the help of the Laplace transform coupled with the classical method for the ordinary differential equations. The general solutions of the system of the equations (52) are p p 1 ¯ P (s), (75) u¯i (y, s) = Di (s) exp −y wi (s) + Ei (s) exp y wi (s) + ai s where
wi (s) =
ai s(1 + λi sαi ), bi
i ∈ In1 .
(76)
Now with Eq. (25) and (75), the Laplace transformed form of the shear stresses τi (y, t), i ∈ In1 , can be written as τ¯i (y, s) =
h i p p p p bi −D (s) w (s) exp(−y w (s)) + E (s) w (s) exp(y w (s)) . (77) i i i i i i 1 + λi sαi
Using the boundary conditions (28) and the interface conditions (29), we have M1 (s)X1 (s) = C1 (s), 17
(78)
where the matrix M(s) is given by 1 0 0 ··· 0 e −f 0 · · · 0 1 1 0 e −f2 · · · 0 2 .. . .. . . .. .. .. . . 0 0 · · · e −f n−1 n−1 M1 (s) = c −d 0 · · · 0 1 1 0 c −d2 · · · 0 2 . . .. . . .. .. .. .. . 0 0 · · · cn−1 −dn−1 0 0 ··· 0 cn
1 1 e1
0 .. .
0 − f11 1 e2
0 0 − f12 .. .
··· ··· ··· .. .
0 0 0 .. .
.. . 1 1 0 0 · · · en−1 − fn−1 g1 −h1 0 ··· 0 0 g2 −h2 · · · 0 .. .. .. .. .. . . . . . 0 0 · · · gn−1 −hn−1 0 0 ··· 0 hn
and the columns X(s) and C(s) are
−dj
ej (s) = e Here gj (s) = hj (s) =
p
√
D1 D2 .. .
X1 (s) = wj (s)
Dn E1 E2 .. . En −dj
, fj (s) = e
wj (s)edj
√
wj (s)
P0 (s) P1 (s) .. .
Pn−1 (s) , C1 (s) = 0 0 .. . 0
√
wj+1 (s)
b
√
,
.
√ p −dj wj (s) , cj (s) = − wj (s)e , (1+λ sαj )
j+1 j , dj (s) = − bj (1+λ αj+1 ) j+1 s
bj+1 (1+λj sαj ) p wj+1 (s)edj bj (1+λj+1 sαj+1
wj+1 (s)
√ p wj+1 (s)e−dj wj+1 (s) ,
,
¯ ¯ 1 − a1j ) P (s) P0 (s) = f¯ − Pa(s) , Pj (s) = ( aj+1 , j = 1, 2, ..., n − 1. s 1s Similar to the case of thermal transport, the system (78) can be written in equivalent form as, M1 D + N1 E = C3 , (79) Q1 D + R1 E = C4 .
18
where .. . N1 = Nij1 i,j∈I 0 , n−1 .. . 1 0 1 0 .. N0j = δ0,j , j ∈ In−1 , M0j = δ0,j , j ∈ In−1 , . .. Mij1 = ei δi,j+1 − fi δi+1,j+1 , . Nij1 = e1i δi,j+1 − f1i δi+1,j+1 , .. . 0 1 0 1 .. , , j ∈ In−1 i ∈ In−1 , , j ∈ In−1 i ∈ In−1 . .. 1 . R1 = Rij , Q1 = Q1ij i,j∈I 0 , 0 i,j∈In−1 n−1 .. . 1 = gi+1 δi+1,j+1 − hi+1 δi+2,j+1 , Q1ij = ci+1 δi+1,j+1 − di+1 δi+2,j+1 , ... Rij .. 0 0 0 0 . , , j ∈ In−1 , i ∈ In−2 , j ∈ In−1 i ∈ In−2 .. . 0 1 0 .. , = δn−1,j , j ∈ In−1 R(n−1)j , Q1(n−1)j = cn δn−1,j , j ∈ In−1 . M1 = Mij1
C3 =
P0 (s) P1 (s) .. . Pn−1 (s)
n×1
, C4 =
0 i,j∈In−1
0 0 .. . 0
,
n×1
,D =
D1 D2 .. . Dn
n×1
,E =
E1 E2 .. . En
, n×1
where δij is the Kronecker tensor. Matrices M1 , N1 , Q1 , R1 are triangular invertible matrices, therefore Eqs. (79) can be written as (N1−1 M1 − R1−1 Q1 )D = N1−1 C3 , (80) E = −R1−1 Q1 D. Suppose that matrix S1 = N1−1 M1 − R1−1 Q1 is invertible, we have D = S1−1 N1−1 C3 , E = −R1−1 Q1 S1−1 N1−1 C3 ,
(81)
where, 0 M1 −1 = (m1ij )i,j∈In−1 ,
0 m10j = δ0,j , j ∈ In−1 , m1i0 =
m1ij = −
i Q
m=1
em fm
j−1 Q `=0
f` e`+1
n−1 P k=0
i Q
`=1
e` ,i f`
1 ∈ In−1 ,
1 δi,j+k , f0 = 1, i, j ∈ In−1 ;
19
(82)
0 N1−1 = (n1ij )i,j∈In−1 ,
0 n10j = δ0,j , j ∈ In−1 , n1i0 =
nij = −
i Q
m=1
fm em
j−1 Q `=0
e`+1 f`
n−1 P k=0
i Q
`=1
f` ,i e`
1 ∈ In−1 ,
1 ; δi,j+k , f0 = 1, i, j ∈ In−1
0 R1−1 = (ρ1ij )i,j∈In−1 ,
ρ1ij
=
n−1 P k=0
i Q
m=0
gm hm
j Q
`=0
h` g`+1
0 δi−k,j , g0 = h0 = gn = 1, i, j ∈ In−1 ,
and the matrix S1 = (Sij1 )i,j=0,1,...,n−1 is defined by the elements Sij1
=
n−1 X k=0
1 n1ik Mkj − ρ1ik Q1kj .
The Laplace transforms given by Eqs. (75), (77) and (81) are complicated, therefore analytical methods are difficult to use for the inverse Laplace transforms of these functions. Here we used two accuracy numerical algorithms, namely the fixed Talbot algorithm [31] and the improved Talbot algorithm [32] for numerical Laplace transform inversion defined in Eqs. (48)-(51).
4
Numerical results and discussions
Unsteady multi-layer laminar flows with heat transfer of Maxwell fluids within a rectangular channel have been investigated by using fractional constitutive equations for the shear stress and thermal flux. The motion of fluids is generated by the channel walls, which are moving with the time-dependent velocities in their planes, and by the time-dependent pressure gradient in the flow direction. The generalized constitutive equations, used in this work, are described by the Caputo time-fractional derivative, therefore, the histories of the velocity gradient and histories of temperature gradient influence the behavior of fluids. Analytical solutions for the Laplace transforms of temperature and velocity fields have been obtained by employing the Laplace transform and Talbots algorithms for the numerical inversion of the Laplace transform. Closed form for fluids velocity has been determined by coupling the Laplace transform with finite sine-Fourier transforms. Analytical and semianalytical solutions have been obtained along with the hypotheses that interfacial velocities, shear stresses, thermal fluxes and temperatures are equal for adjacent fluids. The obtained results have a general character; therefore many particular cases can be obtained. Multi-layer flows of ordinary/fractional Maxwell fluids can be studied as particular cases by letting some fractional parameters equal with one. Using the software Mathcad, numerical calculations were carried out for the obtained solutions. Numerical results are presented in Figs. 2- 13. 20
Figs. 2 - 5 have been plotted in order to analyze the influence of the thermal fractional parameters βi , i = 1, 2, 3 on the layers temperature. The considered fluids are characterized by the following material parameters ρ1 = 1000, µ1 = 0.05, G1 = 1.2, cp1 = 0.25, k1 = 2.5, d1 = 0.3, ρ2 = 1300, µ2 = 0.2, G2 = 2.4, cp2 = 0.27, k2 = 2.7, d2 = 0.7, ρ3 = 1500, µ3 = 0.6, G3 = 5.0, cp3 = 0.29, k1 = 2.9, d1 = 1.0. Profiles of layers temperature, given by Eq. (31), versus the spatial coordinate y, for different values of the thermal fractional parameters are sketched in Fig. 2 for small time, t = 0.02, and Fig. 3 for large time, t = 2. It is observed from Fig. 2 that the variation of the thermal fractional parameter β1 of the fluid situated in the first layer, has significant influence on the temperature in the first and second layers but, the influence on the temperature from the third layer is insignificant. This behavior is due to the decreasing of temperature in the first and second layer with the fractional parameter β1 . It is known from Eq. (22) that the memory kernel of thermal flux is β−1 the function hβ (t) = tΓ(β) whose graph is given in Fig. 6. It is observed from Fig 6 that, for t = 0.02 and β1 ≥ 0.4 the kernel hβ (t) is decreasing with β1 , therefore the weight function of the thermal gradient decreases (the memory effects are weaker). An opposite behavior has been found for the large time when the memory kernel increases with the fractional parameter (See Fig 6, for t = 2). In this case, memory effects are stronger and the layers temperature increase with the fractional parameter β1 (See Fig. 3). Time-evolution of layers temperature, in different positions of the channel and for different values of the thermal fractional parameters are presented in Figs. 4 and 5. As seen in Figs. 4 and 5, the temperature as function of the fractional parameter, has different behavior for small, respectively for large values of the time t. This behavior is in accordance with the previous results presented in Figs. 2 and 3 and with the evolution of the thermal memory kernel hβ (t). We must note that the layers temperatures have small variations in time and after a low value of the time t (t = 1 in analyzed cases), the temperatures are almost constant. Also, from numerical simulations, was observed that the influence of the second/third thermal fractional parameter on the fluid temperature is smaller than the influence of the first fractional parameter. This behavior is due to the thermal boundary conditions because, only the wall y = 0 is heated, so the interfacial temperatures decreases away from this wall. To investigate the velocity of fluids of the three layers, we considered the particular case when the velocities of channel walls are constant, namely, u1 (0, t) = f (t) = 0.5H(t), u3 (1, t) = g(t) = 0.4H(t), H(t) being the unit step Heaviside function and the pressure gradient is P (t) = 1 + sin(t). Figs. 7-9 were plotted in order to study the influence of the velocity fractional parameters αi on the velocity fields. The fractional parameters have effects of braking. The flows of the first two layers are slowed down but, the fluid in the third layer is accelerated. Even if, the wall y = 0 moves with the constant velocity, fluid velocity decreases due viscosities and the tα b α−1 damping of the velocity gradient with the kernel χα (t) = λ t Eα,α − λ . The profiles of the function χα (t) are drawn in Fig. 13. For t = 1, the values of stress kernel decrease with α ≥ 0.5, therefore the memory effects are weakened. The influence of the viscosities ratios bi = µµ1i on the fluid velocity is analyzed in Figs. 1012. It is observed that, for decreasing values of the ratio bi , fluid velocity increases. Indeed, in this case the viscosities of fluids in second and third layer are smaller than viscosity of the 21
first fluid; therefore after the first interface the viscous forces are smaller than the viscous force in the first layer. For increasing values of the ratio bi , fluid flows slower. In order to verify the accuracy of numerical results, we used the analytical solution of velocity given by Eq. (72) and the semi-analytical solution given by Eq. (75) for following values of parameters: d0 = 0.3, d1 = 0.7, d2 = 1, t = 2.5, α1 = 0.8, α2 = 0.7, α3 = 0.4, µ1 = 0.6, µ2 = 0.3, µ3 = 0.6 Numerical results are given in Table 1 and Fig. 14. There is a good agreement between the results obtained with the two formulas.
5
Conclusions
Unsteady multi-layer flows with heat transfer of Maxwell immiscible fluids with generalized constitutive equations for the shear stress and thermal flux in a rectangular channel have been investigated. The generalized constitutive equations are described by the Caputo time-fractional derivative, therefore the fluids behavior is influenced by the histories of the velocity gradient and of temperature gradient. Analytical and semi-analytical solutions for the temperature, velocity and shear stress fields are obtained by employing the Laplace and finite sine-Fourier transforms and hypothesis that interfacial velocities, shear stresses, thermal fluxes and temperatures are equal for adjacent fluids. The obtained results have a general character; therefore many particular cases can be obtained. Multi-layer flows of ordinary/fractional Maxwell fluids are studied as particular cases by letting some fractional parameters equal with one. Numerical values of the temperature and velocity fields are obtained by using the Mathcad software and the numerical Talbots algorithms for the inverse Laplace transforms. Numerical results are used for graphical illustrations of the fluid temperature and velocity. The heat transfer and fluid motions are different for small, respectively large values of the time t. These different behaviors are due to the thermal/ velocity kernels variations with respect to time and fractional parameters, therefore the memory effects have significant influence on the fluids. This behavior is in agreement with the variation of the thermal kernel hβ (t) (See Fig. 6). For the values of the thermal fractional parameter β ≤ 0.35 and small values of the time t, the thermal kernel hβ (t) increases with β therefore, the values of the temperature gradient are strongly damped and the heat transfer is reduced. For the thermal fractional parameter β > 0.35, the values of thermal kernel decreases, the damping of the thermal gradient is more weakly and the heat transfer increases. The stress damping kernel χα (t) has the variation similar to the thermal damping kernel (See Fig. 13), therefore, the effect on the fluid velocity is stronger for small values of the stress fractional parameter α.
References [1] D. K. Satpathi, B. V. Rathish, P. Chandra, Unsteady-state laminar flow of viscoelastic gel and air in a channel : Application to mucus transport in a cough machine simulating trachea, math. Comput. Modelling, 38 (2003) 63-75, DOI: 10.1016/S0893-7177(03)00192-4.
22
Figure 2: Temperature profile Ti (y, t) versus y for small values of time and different values of fractional parameters.
Table 1: Comparison between results obtained with analytical and semi-analytical solutions
23
Figure 3: Temperature profile Ti (y, t) versus y for large values of time and different values of fractional parameters.
24
Figure 4: Time-variation of temperature profile Ti (y, t) with the fractional parameter β1 .
25
Figure 5: Time-variation of temperature profile Ti (y, t) with the fractional parameter β2 .
26
Figure 6: Variation of the kernel hβ (t) = tβ−1 /Γ(β) with the fractional parameter β.
Figure 7: Profiles of velocities ui (y, t) for the fractional parameter α1 variable.
27
Figure 8: Profiles of velocities ui (y, t) for the fractional parameter α2 variable.
Figure 9: Profiles of velocities ui (y, t) for the fractional parameter α3 variable.
28
Figure 10: Profiles of velocities ui (y, t) for µ1 variable.
Figure 11: Profiles of velocities ui (y, t) for µ2 variable.
29
Figure 12: Profiles of velocities ui (y, t) for µ3 variable.
Figure 13: Stress kernel χα (t).
30
Figure 14: Comparison between results obtained with analytical and semi-analytical solutions. [2] D. D. Joseph, Y. Y. Renardy, Fundamentals of two-fluid dynamics, Part I: Mathematical Theory and Applications, Springer-Verlag, 1992. [3] Y. Y. Su, B. Khomami, Interfacial stability of multilayer viscoelastic fluids in slit and converging channel dies geometries, J. Rheol., 36 (1992) 357-387. [4] N. R. Anturkar, T. C. Papanastasiou, J. O. Wilkes, Stability of multilayer extrusion of viscoelastic liquids, AIChE Journal, 36 (5) (1990), 710-724. [5] A. C. B. Bogaerds, M. A. Hulsen, G. W. M. Peters, F. P. T. Baaijens, Time dependent finite element analysis of the linear stability of viscoelastic flows with interfaces, J. Non-Newtonian Fluid Mech. 116 (2003), 33-54. [6] H. Le Meur, Non-uniqueness and linear stability of the one-dimensional flow of multiple viscoelastic fluids, Math. Modell. Numer. Anal. 31 (2) (1997) 185-212. [7] P. A. Vasquez, Y. Jin, E. Palmer, D. Hill, M. G. Forest, Modeling and simulation of mucus flow in human bronchial epithelial cell cultures-Part I: Idealized Axysymmetric swirling flow, PLoS Comput. Biol. 12 (8) (2016) e1004872, doi: 10.1371/journal.pcbi.1004872. [8] L. L. Barannyk, D. T. Papageorgiou, P. G. Petropoulos, J. M. Vanden-Broeck, Nonlinear dynamics and wall touch-up in unstability stratified multilayer flows in horizontal channel, SIAM J. Appl. Math. 75 (1) (2015) 92-113, doi: 10.1137/140968070. [9] E. S. Papaefthymiou, D. T. Papageorgiou, Nonlinear stability in three-layer channel flows, J. Fluid Mech., 829 R2 (2017), doi: 10. 1017/jfm.2017.605. [10] J. Zhao, L. Zheng, X. Zhang, F. Liu, Unsteady boundary layer natural convection heat transfer of fractional Maxwell viscoelastic fluid over a vertical plate, Int. J. Heat Mass Transfer, 97 (2016)760-766. [11] X. Chen, W. Yang, X. Zhang, F. Liu, Unsteady boundary layer flow of viscoelastic MHD fluid with a double fractional Maxwell model, Appl. Math. Letters, 95 (2019) 143-149. [12] W. Chen, J. Zhang, J, Zhang, A variable-order time-fractional derivative model for chloride ions subdiffusion in concrete structures, Fractional Calculus and Appl. Ananl., 16(1) (2013) 7692.
31
[13] X. H. Ying, J. X. Yun, Time fractional dual-phase lag heat conduction equation, Chin. Phys. B, 24 (3), (2015) 034401. [14] Muhammad Danial Hisham , Abdul Rauf , Dumitru Vieru , Aziz Ullah Awan , Analytical and semianalytical solutions to flows of two immiscible Maxwell fluids between moving plates, Chinese Journal of Physics (2018), doi: https://doi.org/10.1016/j.cjph.2018.10.009 [15] I. Podlubny, Fractional Differential Equations, Vol. 198 of Mathematics in Science and Engineering, Academic Press, San Diego, USA, 1999. [16] M. Caputo, Linear model of dissipation whose Q is almost frequency independent II, The Geophys. J. R. A. Soc., 13 (1967) 529-539. [17] M. Caputo, M. Fbrizio, A new definition of fractional derivative without singular kernel, Progr. Fract. Differ. Appl., 1 (2) (2015) 73-85. [18] 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. [19] T. R. Prabhakar, A singular integral equation with a generalized Mittag-Leffler function in the kernel, Yokohama Math. J., 19 (1971) 7-15. [20] J. Hristov, Transient space-fractional diffusion with a power-law superdiffusivity: Approximate integralbalance approach, Fundamenta Informaticae, 151 (2017) 371-388, doi: 10.3233/FI-2017-1498. [21] 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. [22] Y. Jiang, H. Qi, H. Xu, X. Jiang, Transient electroosmotic slip flow of fractional Oldroyd-B fluids, Microfluid nanofluid, 21 : 7 (2017) doi : 10.1007/s10404-016-1843-x. [23] 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. [24] X. Yang, H. Qi, X. Jiang, Numerical analysis for electroosmotic flow of fractional Maxwell fluid, Appl. Math. Lett., 78 (2018) 1- 8. [25] Y. Z. Povstenko, Fractional heat conduction equation and associated thermal stress, J. Thermal Stresses, 28 (1) (2004) 83-102. [26] J. Hristov, Integral-balance approach to 1-D space-fractional diffusion models, Mathematical Methods in Engineering Applications in Dynamics of Complex Systems, Springer, 2018. [27] J. Hristov, Integral-Balance Solution to Nonlinear Subdiffusion Equation, Ch. 3, in Frontiers in Fractional Calculus, 1 (2017), 71-106. [28] 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 Kenan Tas, Dumitru Baleanu and J. A. T. Machado, Springer, 2018. [29] D. Brian, Integral Transforms and Their Applications (third ed.) New York: Springer, 2002. [30] R. N. Bracewell, The Fourier Transform and its Applications (Third ed.) Boston: Mc Graw-Hill, 2000. [31] J. Abate, P. P. Valko, Multi-precision Laplace transform inversion, Int. J. Numer. Meth. Engng., 60 (2004) 979-993, doi:10.1002/nme.995. [32] B. Dingfelder, J. A. C. Weideman, An improved Talbot method for numerical Laplace transform inversion, Numer. Algor., 68 (2015) 167-183, doi: 10.1007/s11075-014-9895-z.
32
[33] C. F. Lorenzo, T. T. Hartley, Generalized functions for the fractional calculus, NASA/TP-1999209424/REV1. [34] 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.
33