# Simultaneous Flow of n-Immiscible Fractional Maxwell Fluids with Generalized Thermal Flux and Robin Boundary Conditions

Authors: Abdul Rauf; Qammar Rubbab; Nehad Ali Shah; Kaleem Razzaq Malik
Journal: Advances in Mathematical Physics
 (2021)

---

## Abstract

In a rectangular region, the multilayered laminar unsteady flow and temperature distribution of the immiscible Maxwell fractional fluids by two parallel moving walls are studied. The flow of the fluid occurs in the presence of Robin’s boundaries and linear fluid-fluid interface conditions due to the motion of the parallel walls on its planes and the time-dependent pressure gradient. The problem is defined as a mathematical model which focuses on the fluid memory, which is represented by a constituent equation with the Caputo time-fractional derivative. The integral transformations approach (the Laplace transform and the finite sine-Fourier transform) is used to determine analytical solutions for velocity, shear stress, and the temperature fields with fluid interface, initial, and boundary conditions. For semianalytical solutions, the algorithms of Talbot are used to calculate the Laplace inverse transformation. We used the Mathcad software for graphical illustration and numerical computation. It has been observed that the memory effect is significant on both fluid motion and temperature flow.


---

## Body

## 1. Introduction

In nature, there often exist flows of immiscible materials. Due to its broad application in research, medicine, geophysics, industry, petroleum engineering, and hydrogeology, the study of simultaneous flow of two or more immiscible fluids is significant [1–4]. The numerous applications include the recovery of petroleum oil, blood flow through the veins of a capillary vessel, the treatment of machinery, the processing of organic film and mucus in living cells, the removal of carbon dioxide from the environment, the control of groundwater, crude oil pipeline flows, and the formation of blisters in microfluid and bubble trains.In some industrial problems, fluid flow is multicomponent, and therefore, there are layers of fluids having different densities and viscosities. The interface of these layers creates moving boundaries in between the walls of the channel in which fluid is flowing. This causes flow phenomenon to be not only nonlinear but also very complex and its study challenging.A long-wave technique was used for the first study of the linear stability of the viscoelastic two-layered simultaneous Poiseuille and Couette flow by Yih [5]. It has been seen that the Kelvin-Helmholtz instability can occur due to viscosity and density stratification. Several scientists subsequently researched the stability/instability of the immiscible fluid flow in two or multifaceted layers [6–8]. The existence and uniqueness of the simultaneous multilayered Couette/Poiseuille fluid motions in channel/pipes were investigated by Le Meur [9], and the approximated Oldroyd differential component and the viscosity proportions were found important to a unique result. The Couette-Poiseuille motion of the two-layer fluids was taken into account by Kalogirou and Blyth [10] to examine stability. The fluid in the bottom layer is surrounded by surfactants which adsorbed on the interface. The thickness ratio of the fluid viscosity is considerably higher, and when the flow is relatively stable, the surfactant is soluble enough.Immiscible fluid flows are frequently encountered in the design of industrial processes and equipment. Indeed, several flow patterns of interest exist for different flow conditions in liquid-liquid flows. In [11], the two-phase flow of immiscible fluids in porous media is described in continuous models as the momentum exchanges between the two steps by simply making generalizations of the Darcy law for single phases of development and adding saturation-dependent porosity and permeability. By extending the principle of Buckley-Leverett, authors examine the effect of the cross-cutting words on concentration profiles and pressure losses for various fluids. It has been shown that the outcomes on dual-phase flow might increase when the impact of the fluid-fluid interfaces appears close to that of the solid-fluid interfaces with the permeability of the porous medium. Hisham et al. [12] reported a two-layer analysis in the presence of time-dependent pressure gradient of the immiscible Maxwell fluids between two simultaneous moving plates. With the support of integral transformation, Laplace and finite Fourier Sine transformation, analytical solutions for velocity, and shear stress are retrieved. The increase in kinematic viscosity decreases the maximum speed value. In a rectangular channel with two parallel translating plates in the presence of a time-dependent pressure gradient, Rauf et al. [13] proposed analytical and semianalytical solutions for the velocity fields and temperature fields for the simultaneous flow of the n-immiscible Maxwell fractional fluid. It has been found that the heat transfer in ordinary fluids is higher than those of thermal memory fluids; however, the memory parameters affect the fluids’ velocities as accelerating factors. Some other important multilayer flow problems are studied in [14–16]. The parallel flow of the fractional Maxwell fluid inside a cylindrical domain has been studied by Rauf et al. [17] in the presence of pressure gradient in the flow direction. The Laplace transformations, in combination with the finite Weber and Hankel transformation of zeroth order provide analytical solutions to flow velocities and shear stresses. The fluid velocity was shown to decrease as the values of the fractional parameters have been increased. The study of multilayer flow of generalized immiscible Maxwell fluids between two parallel plates with Robin boundary conditions is still lacking in the literature. The aim of this article is to fill this gap in the literature.The fractional operators have been intensively incorporated in recent decades for mathematical modeling of the numerous topics in real life. Many dynamic memory processes can be investigated using time-fractional derivative operators. The fractional calculus thus became an essential component in biology, chemistry, physics, and many areas of engineering. Hristov [18, 19], Povstenko and Kyrylych [20, 21], Zheng et al. [22], Baleanu et al. [23, 24], and Hilfer [25] have a comprehensive description of the theory of fractional operators with their properties and applications. References [26–38] include mathematical observations and potential implementations and applications of the fractional differential Calculus.In this paper, we have studied the multilayer flow between the two parallel plates of immiscible fractional Maxwell fluids. In the vicinity of the fractional heat flux in fluid layers, we assumed an unsteady, incompressible, and fully established one-dimensional fluid motion induced by the motion of the boundary plates and by the applied pressure gradient as a function of time. Moreover, we took into account the Robin boundary conditions on the boundary plates and the linear interfacial fluid-fluid conditions between two consecutive layers. We have used finite Fourier sine-cosine transformation, which is ideal for the Robin-type boundary conditions, along with the Laplace transformation, to explore analytical solutions for velocities and shear stresses. Using the Laplace transform along with Tablot’s techniques for the numerical inversion of Laplace, a semianalytical solution is recovered for thermal profile.


## 2. Mathematical Modeling

The flow region is(1)D′=x′,y′,z′x′,z′∈−∞,∞,y′∈0,h′,with the boundary plates positioned at y′=0 and y′=h′>0. To begin with at time t′=0, the two walls and the fluids inside are at rest with the atmospheric temperature T2o. After this instance, the boundary plate positioned at y′=0 move with velocity u′1o=U0f′1t′ along the x′-axis and the wall temperature Tog′1t′, while the channel plate y′=h′ moves with the wall velocity u′2o=U0f′2t′ analogous with the x′-axis and the wall temperature Tog′2t′ (Figure 1). It has been assumed that the functions f1′t′,f′2t′, g′1t′, and g′2t′ are piece-wise continuous functions, and f′10=g′10=f′20=g′20=0. We consider that the velocity vector is of the form V′→=u′y′,t′,0,0. We assume the simultaneous n-immiscible Maxwell fluids between two parallel boundary planes. It divides the domain of flow 0,h′ into n subdomains 0=h′0<h′1<h′2⋯<h′n=h′. In the region y′∈hi−1′,hi′,h′i−1<h′i, Maxwell fluid flow has the viscosity μi, relaxation time λ′i=μi/Gi, density ρi, Gi, the elastic modulus, temperature T′iy′,t′, velocity u′iy′,t′, the shear stress τ′iy′,t′, and the thermal flux q′iy′,t′, where i=1,2,⋯,n.Figure 1
Geometry of the problem.The Maxwell fluids are supposed to be incompressible and immiscible, and the stream is linear, unsteady, and completely established. The fluid flow is produced by the applied pressure gradient as a function of time towards flow and by the motion of the boundary walls. In view of the assumptions on the velocity vectors, the continuity equation is indistinguishably fulfilled by all velocitiesu′iy′,t′, i∈In1, where In1≔1,2,⋯,n, and we assume that In−11≔1,⋯,n−1. The system of constitutive equations of motion and the initial, boundary, and the fluid-fluid interface conditions are given, respectively, as follows:
(i)
The n-linear momentum equations:(2)ρi∂ui′∂t′=∂τi′∂y′−∂p′∂x′,i∈In1.(ii)
The n-constitutive equations:(3)τi′+λ′i∂τi′∂t′=μi∂ui′∂y′,i∈In1.(iii)
The corresponding n-initial conditions:(4)ui′y′,0=0,τi′y′,0=0,i∈In1.(iv)
The boundary condition:(5)atthelowerwally′=0,u′10,t′−l′0∂u1′y,t′∂y′y′=0=U0f′1t′,attheupperwally′=hn′=h′,un′h′,t′+ln′∂un′y′,t′∂y′y′=h′=U0f2′t′.(v)
At the interface, we assume the continuity of Robin conditions:(6)ui′hi′,t′+li′∂ui′y′,t′∂y′y′=hi′=ui+1′hi′,t′−li′∂ui+1′y′,t′∂y′y′=hi′,τi′hi′,t′=τi+1′hi′,t′,i∈In−11.The system of constitutive relations for temperature and heat flux is given by the following:(i)
The n-thermal transport balance equations:(7)ρicip∂Ti′y′,t′∂t′=−∂qi′y′,t′∂y′,i∈In1.(ii)
The n-thermal fluxes with Fourier’s law:(8)qi′y′,t′=−ki∂Ti′y′,t′∂y′,i∈In1,where ki and cip are the thermal conductivity and the specific heat at constant pressure, respectively. We consider the following:
(i)
n-Initial conditions:(9)Ti′y′,0=T2o,qi′y′,0=q0,i∈In1.(ii)
Boundary conditions:(10)T1′0,t′−T2o−δ′0∂T1′y,t′∂y′y′=0=Tog′1t′,Tn′h′,t′−T2o+δ′n∂T′ny′,t′∂y′y′=h′=Tog2′t′.(iii)
Interface conditions:(11)Ti′hi′,t′−T2o+δ′i∂Ti′y′,t′∂y′y′=hi′=Ti+1′hi′,t′−T2o−δ′i∂Ti+1′y′,t′∂y′y′=hi′,q′ih′i,t′=q′i+1h′i,t′,i=1,2,⋯,n−1.Consider the nondimensional variables(12)x=x′h′,y=y′h′,t=ν1t′h′2,ui=ui′U0,τi=h′τi′μ1U0,νi=μiρi,p=h′p′μ1U0,λi=λ′iν1h′2,ai=ρiρ1,bi=μiμ1,f1t=f1′h′2tν1,g1t=g1′h′2tν1,Ti=Ti′−T2oTo,qi=qi′q0,σi=νiν1,Pri=μicipki,f2t=f2′h′2tν1,hi=hi′h′,g2t=g2′h′2tν1,lj=l′jh′,δj=δ′jh′i∈In1,j=0,1,2,⋯,n.The dimensionless form of the governing Equations (2)–(11) are reduced into the following form:
(13)ai∂ui∂t=−∂p∂x+∂τi∂y,i∈In1,(14)τi+λi∂τi∂t=bi∂ui∂y,i∈In1,(15)ν1γ0ρicip∂Ti∂t=−∂qi∂y,i∈In1,(16)qi=−γ0ki∂Ti∂y,i∈In1,along with the following:
(i)
The dimensionless n-initial conditions:(17)Tiy,0=0,qiy,0=1,uiy,0=0,τiy,0=0,i∈In1.(ii)
The dimensionless boundary condition:(18)onthebottomplatey=0,u10,t−l0∂u1y,t∂yy=0=f1t,T10,t−δ0∂T1y,t∂yy=0=g1t,(19)ontheupperplatey′=hn=1,un1,t+ln∂uny,t∂yy=1=f2t,Tn1,t+δn∂Tny,t∂yy=1=g2t.(iii)
The nondimensional interface conditions:(20)uihi,t+li∂uiy,t∂yy=hi=ui+1hi,t−li∂ui+1y,t∂yy=hi,τihi,t=τi+1hi,t,(21)Tihi,t+δi∂Tiy,t∂yy=hi=Ti+1hi,t−δi∂Ti+1y,t∂yy=hi,qihi,t=qi+1hi,t,i∈In−11,where γ0=To/q0h′ and To is the characteristic scale.


### 2.1. Generalized Constitutive Mathematical Model

Consider the following generalized constitutive mathematical relations:(22)τi+λiDtαiτi=bi∂ui∂y,αi∈0,1,i∈In1,(23)qi=−γ0kiDt1−βi∂Ti∂y,βi∈0,1,i∈In1,where the Caputo derivative Dtσ is defined as
(24)DtσFy,t=1Γ1−σ∫0tt−ρ−σ∂Fy,ρ∂ρdρ,0≤σ<1.For the special case=1,Dt1Fy,t=∂Fy,t/∂t. It is significant to note that the fractional models (22) and (23) have the following equivalent formulae as
(25)τiy,t=biλitαi−1Eαi,αi−tαiλi∗∂uiy,t∂y=∫0tbiλit−ταi−1Eαi,αi−t−ταiλi∂uiy,τ∂ydτ,i∈In1,(26)qiy,t=tβi−1Γβi∗∂2Tiy,t∂t∂y=∫0tt−τβi−1Γβi∂2Tiy,τ∂τ∂ydτ,i∈In1.The relations (25) and (26) depict that the histories of the thermal and the velocity gradients impact the time-variation of heat flux and the shear stress, separately. Moreover, it has been seen from the above relations that the nonlocality kernel of the heat flux observes the power-law tβi−1/Γβi, while the nonlocality kernel for the shear stress is the function bi/λitαi−1Eαi,αi−tαi/λi,Eαi,αi. being the Mittag-Leffler function. In this study, we consider the pressure gradient in the flow direction as a known function, namely
(27)−∂p∂x=Pt,where Pt is a peice-wise continuous function on 0,∞.


## 2.1. Generalized Constitutive Mathematical Model

Consider the following generalized constitutive mathematical relations:(22)τi+λiDtαiτi=bi∂ui∂y,αi∈0,1,i∈In1,(23)qi=−γ0kiDt1−βi∂Ti∂y,βi∈0,1,i∈In1,where the Caputo derivative Dtσ is defined as
(24)DtσFy,t=1Γ1−σ∫0tt−ρ−σ∂Fy,ρ∂ρdρ,0≤σ<1.For the special case=1,Dt1Fy,t=∂Fy,t/∂t. It is significant to note that the fractional models (22) and (23) have the following equivalent formulae as
(25)τiy,t=biλitαi−1Eαi,αi−tαiλi∗∂uiy,t∂y=∫0tbiλit−ταi−1Eαi,αi−t−ταiλi∂uiy,τ∂ydτ,i∈In1,(26)qiy,t=tβi−1Γβi∗∂2Tiy,t∂t∂y=∫0tt−τβi−1Γβi∂2Tiy,τ∂τ∂ydτ,i∈In1.The relations (25) and (26) depict that the histories of the thermal and the velocity gradients impact the time-variation of heat flux and the shear stress, separately. Moreover, it has been seen from the above relations that the nonlocality kernel of the heat flux observes the power-law tβi−1/Γβi, while the nonlocality kernel for the shear stress is the function bi/λitαi−1Eαi,αi−tαi/λi,Eαi,αi. being the Mittag-Leffler function. In this study, we consider the pressure gradient in the flow direction as a known function, namely
(27)−∂p∂x=Pt,where Pt is a peice-wise continuous function on 0,∞.


## 3. Solution of the Problem

We have applied the finite sine-Fourier transform coupled with the Laplace transform, to obtain explicit solutions of Equations (13), (15), (22), and (23) with conditions (17)–(20). We apply the Laplace transform to Equations (13), (15), and (18)–(23) along with the initial condition (17); we have
(28)aisu¯iy,s=∂τ¯iy,s∂y+P¯s,i∈In1,(29)τ¯iy,s=bi1+λisαi∂u¯iy,s∂y,i∈In1,(30)ν1γ0ρicipsT¯iy,s=−∂q¯iy,s∂y,i∈In1,(31)q¯i=−γ0kis1−βi∂T¯iy,s∂y,i∈In1,(32)u¯10,s−l0∂u¯1y,s∂yy=0=f¯1s,u¯n1,s+ln∂u¯ny,s∂yy=1=f¯2s,(33)T¯10,s−δ0∂T¯1y,s∂yy=0=g¯1s,T¯n1,s+δn∂T¯ny,s∂yy=1=g¯2s,(34)u¯ihi,s+li∂u¯iy,s∂yy=hi=u¯i+1hi,s−li∂u¯i+1y,s∂yy=hi,(35)τ¯ihi,s=τ¯i+1hi,s,i=1,2,⋯,n−1,(36)T¯ihi,s+δi∂T¯iy,s∂yy=hi=T¯i+1hi,s−δi∂T¯i+1y,s∂yy=hi,(37)q¯ihi,s=q¯i+1hi,s,i=1,2,⋯,n−1.


### 3.1. Analytical Solutions for Velocity and Shear Stress

Combining Equations (28) and (29), we can write
(38)ais1+λisαiu¯iy,s=1+λisαiP¯s+bi∂2u¯iy,s∂y2,i∈In1,where φ¯y,s=∫0∞φy,texp−stdt represents the Laplace transform of the function φy,t. In case of nonhomogeneous Robin boundary and fluid-fluid interface conditions (32) and (34), the finite Fourier sine-cosine transform of the function u¯iy,s, i∈In1, can be defined with the help of the Fourier series theory and the Sturm-Liouville theory as [39, 40].
(39)u¯^iηmi,s=1hi−hi−1∫hi−1hiu¯iy,sΦmiydy,where Φmiy=sinηmiy−hi−1/hi−hi−1+li−1/hi−hi−1ηmicosηmiy−hi−1/hi−hi−1 and ηmi are the positive roots of the transcendental equation tanηmi=hi−hi−1li−1+liηmi/li−1liηmi2−hi−hi−12. The inverse Fourier sine-cosine transform of u¯^iηmi,s is defined by
(40)u¯iy,s=∑m=1∞u¯^iηmi,sJmiΦmiy,where
(41)Jmi=1hi−hi−1∫hi−1hiΦmi2ydy=li−12ηmi2−hi−hi−124ηmihi−hi−12sin2ηmi+hi−hi−12+li−1hi−hi−1+li−12ηmi2−li−1hi−hi−1cos2ηmi2hi−hi−12.By direct computations, using the robin boundary condition (32) and interface fluid-fluid condition (34), we can write
(42)∂2u¯^i∂y2ηmi,s=ηmihi−hi−12cosηmihi−hi−12+li−12ηmi2li−1liηmi2−hi−hi−12w¯ihi,li,s+v¯ihi−1,li−1,s−ηmi2hi−hi−12u¯^iηmi,s,i∈In1,where we rename v¯ihi−1,li−1,s=u¯ihi−1,s−li−1∂u¯ihi−1,s/∂y, w¯ihi,li,s=u¯ihi,s+li∂u¯ihi,s/∂y, i∈In−11,v¯0h0,l0,s=f¯1s, and w¯nhn,ln,s=f¯2s. With the application of finite Fourier sine-cosine transform (39) to Equation (38) corresponding to the robin boundary conditions (32) and interface fluid-fluid conditions (34), and using Equation (42), the transformed velocities take the form
(43)u¯^iηmi,s=hi−hi−121+λisαiP¯saihi−hi−12s1+λisαi+biηmi2+biηmiaihi−hi−12s1+λisαi+biηmi2×cosηmihi−hi−12+li−12ηmi2li−1liηmi2−hi−hi−12w¯ihi,li,s+v¯ihi−1,li−1,s,i∈In1.To apply the inverse Fourier sine-cosine transform, we rewrite Equation (44) in the following suitable form:
(44)u¯^iηmi,s=li−1ηmisinηmi+hi−hi−11−cosηmiηmihi−hi−1hi−hi−1+liv¯ihi−1,li−1,s+li−1w¯ihi,li,shi−hi−1+li−1+li+−li−1ηmi−hi−hi−1−li−1ηmicosηmi+hi−hi−1+li−1ηmi2sinηmihi−hi−1ηmi2w¯ihi,li,s−v¯ihi−1,li−1,shi−hi−1hi−hi−1+li−1+li+biηmiaihi−hi−12s1+λisαi+biηmi2+C1miv¯ihi−1,li−1,s+biηmicosηmiaihi−hi−12s1+λisαi+biηmi2hi−hi−12+li−12ηmi2li−1liηmi2−hi−hi−12+C2miw¯ihi,li,s+hi−hi−121+λisαiP¯saihi−hi−12s1+λisαi+biηmi2,i∈In1,where
(45)C1mi=−1ηmi+hi−hi−1li−1+liηmicosηmi+hi−hi−12sinηmi−lili−1ηmi2sinηmihi−hi−1hi−hi−1+li−1+liηmi2,C2mi=hi−hi−12ηmicosηmi−hi−hi−12+li−1hi−hi−1+li−1ηmi2sinηmihi−hi−1hi−hi−1+li−1+liηmi2.Considering the auxiliary relations along with their Fourier sine-cosine transformv(46)φ1iy=1,y∈hi−1,hi,φ^1iηmi=li−1ηmisinηmi+hi−hi−11−cosηmiηmihi−hi−1,(47)φ2iy=y−hi−1hi−hi−1,y∈hi−1,hi,(48)φ^2iηmi=−li−1ηmi−hi−hi−1−li−1ηmicosηmi+hi−hi−1+li−1ηmi2sinηmihi−hi−1ηmi2,(49)m=1,2,⋯,i∈In1,h0=0,hn=1,the inverse Fourier sine-cosine transform of Equation (46) takes the form
(50)u¯iy,s=hi−hi−1+liv¯ihi−1,li−1,s+li−1w¯ihi,li,shi−hi−1+li−1+li+w¯ihi,li,s−v¯ihi−1,li−1,shi−hi−1+li−1+liy−hi−1+∑m=1∞biηmiaihi−hi−12s1+λisαi+biηmi2+C1miΦmiyv¯ihi−1,li−1,sJmi+∑m=1∞biηmicosηmiaihi−hi−12s1+λisαi+biηmi2hi−hi−12+li−12ηmi2li−1liηmi2−hi−hi−12+C2miΦmiyw¯ihi,li,sJmi+∑m=1∞hi−hi−121+λisαiaihi−hi−12s1+λisαi+biηmi2ΦmiyP¯sJmi,i∈In1.Now, from (29) and (50), we obtain
(51)τ¯iy,s=T¯i,1y,sw¯ihi,li,s−T¯i,2y,sv¯ihi−1,li−1,s+T¯i,3y,sP¯s,i∈In1,where
(52)T¯i,1y,s=bi1+λisαihi−hi−1+li−1+li+bihi−hi−11+λisαi∑m=1∞ηmibiηmicosηmiaihi−hi−12s1+λisαi+biηmi2hi−hi−12+li−12ηmi2li−1liηmi2−hi−hi−12+C2mi×cosηmiy−hi−1/hi−hi−1−li−1/hi−hi−1ηmisinηmiy−hi−1/hi−hi−1Jmi,T¯i,2y,s=bi1+λisαihi−hi−1+li−1+li−bihi−hi−11+λisαi∑m=1∞ηmibiηmiaihi−hi−12s1+λisαi+biηmi2+C1mi×cosηmiy−hi−1/hi−hi−1−li−1/hi−hi−1ηmisinηmiy−hi−1/hi−hi−1Jmi,T¯i,3y,s=bihi−hi−1P¯s∑m=1∞ηmicosηmiy−hi−1/hi−hi−1−li−1/hi−hi−1ηmisinηmiy−hi−1/hi−hi−1Jmiaihi−hi−12s1+λisαi+biηmi2.Using Equations (32) and (34) in Equation (51), we get the following linear system (for detail, see “Appendix”):
(53)M¯sW¯s=N¯s,where
(54)M¯s=Y¯1,1−Y¯1,200⋯0Y¯2,1Y¯2,2−Y¯2,30⋯00Y¯3,2Y¯3,3−Y¯3,4⋯0⋮⋮⋱⋱⋱⋮00⋯Y¯n−3,n−4Y¯n−3,n−3−Y¯n−3,n−200⋯Y¯n−2,n−3Y¯n−2,n−2−Y¯n−2,n−100⋯0Y¯n−1,n−2Y¯n−1,n−1,W¯s=w¯1h1,l1,sw¯2h2,l2,s⋮w¯n−1hn−1,ln−1,s,N¯s=b¯1b¯2⋮b¯n−1.Finally, we obtain(55)W¯s=M¯−1sN¯s.Now,w¯ihi,li,s=v¯i+1hi,li,s,i∈In−11, are known functions; therefore, the velocities u¯1y,s,⋯,u¯ny,s are known. In order to obtain the inverse Laplace transforms of the functions u¯iy,s,i∈In1, we consider the following auxiliary functions:
(56)H¯i0m,s=1aihi−hi−12s1+λisαi+ηmi2=s−1aihi−hi−12λisαi+λi−1·11+ηmi2s−1/aihi−hi−12λisαi+λi−1=1aiλihi−hi−12∑k=0∞−1kηmi2aihi−hi−12λiks−k−1sαi+λi−1k+1,H¯i1m,s=1+λisαihi−hi−12ais1+λisαi+biϑmi2=1aihi−hi−12∑k=0∞−1kηmi2hi−hi−12aiλiks−k−1sαi+λi−1k,H¯i4m,s=11+λisαiaihi−hi−12s1+λisαi+biϑmi2=1aiλi2hi−hi−12∑k=0∞−1kηmi2aihi−hi−12λiks−k−1sαi+λi−1k+2.Since the generalized G-Lorenzo-Hartley function is defined by [41]
(57)Gσ1,σ2,σ3t,σ=L−1sσ2sσ1−σσ3=∑k=0∞Γk+σ3σktk+σ3σ1−σ2−1k!Γk+σ3σ1−σ2Γσ3,Res>0,Reσ1σ3−σ2>0,σsσ1<1,and for αi,βi>0,(58)L−1sαi−βisαi−d=tβi−1Eαi,βidtαi,where Eαi,βi. is the Mittag-Leffler function [42]. The inverse Laplace transform of H¯i0m,s, H¯i1m,s, and H¯i4m,s takes the form
(59)Hi0m,t=1aiλihi−hi−12∑k=0∞−1kηmi2aihi−hi−12λikGαi,−k−1,k+1t,−λi−1,(60)Hi1m,t=1hi−hi−12ai∑k=0∞−1kηmi2khi−hi−12aiλik+1Gαi,−k−1,kt,−λi−1,(61)Hi4m,t=1aiλi2hi−hi−12∑k=0∞−1kηmi2hi−hi−12aiλikGαi,−k−1,k+2t,−λi−1,where δt is the Dirac delta function. Using Equations (59), (60), and (50), we obtain for velocities uiy,t,i∈In1, with the following expressions:
(62)uiy,t=hi−hi−1+liv¯ihi−1,li−1,t+li−1w¯ihi,li,thi−hi−1+li−1+li+w¯ihi,li,t−v¯ihi−1,li−1,thi−hi−1+li−1+liy−hi−1+∑m=1∞biηmiHi0m,t∗vihi−1,li−1,t+C1mivihi−1,li−1,tΦmiyJmi+∑m=1∞hi−hi−12+li−12ηmi2biηmicosηmili−1liηmi2−hi−hi−12Hi0m,t∗wihi,li,t+C2miwihi,li,tΦmiyJmi+∑m=1∞hi−hi−12Hi1m,t∗PtΦmiyJmi,i∈In1,where h1t∗h2t=∫0th1t−τh2τdτ is the convolution product of the functions h1t and h2t. The system of shear stresses τiy,t, i∈In1, can be determined by applying inverse Laplace transform to Equation (51) and using Equations (59) and (61).
(63)τiy,t=Ti1y,t∗wihi,li,t−Ti2y,t∗vihi−1,li−1,t+Ti3y,t∗Pt,where
(64)Ti,1y,t=biλi−1tαi−1Eαi,αi−λi−1tαihi−hi−1+li−1+li+bihi−hi−1∑m=1∞ηmibiηmicosηmihi−hi−12+li−12ηmi2li−1liηmi2−hi−hi−12Hi4m,t+C2miλi−1tαi−1Eαi,αi−λi−1tαicosηmiy−hi−1/hi−hi−1−li−1/hi−hi−1ηmisinηmiy−hi−1/hi−hi−1Jmi,T¯i,2y,s=biλi−1tαi−1Eαi,αi−λi−1tαihi−hi−1+li−1+li−bihi−hi−1∑m=1∞ηmibiηmiHi4m,t+C1miλi−1tαi−1Eαi,αi−λi−1tαicosηmiy−hi−1/hi−hi−1−li−1/hi−hi−1ηmisinηmiy−hi−1/hi−hi−1Jmi,T¯i,3y,s=bi∑m=1∞ηmihi−hi−1Hi0m,t×cosηmiy−hi−1/hi−hi−1−li−1/hi−hi−1ηmisinηmiy−hi−1/hi−hi−1Jmi.At the end of this section, we mention that, in the case of a single fluid, there are interesting studies in the literature on the flow of Maxwell fluids with slip on the channel walls. Thus, in the particular caseα=1,n=1,Pt=f2t=0, our results are equivalent with those obtained in ([43], Equations (34) and (55)). Moreover for the Newtonian case with n=1,Pt=0,λ1=0, the velocity profile given by Equation (62) is equivalent to the solution obtained in [44], Equation (65) with φ=0 and ρβf=0.


### 3.2. Solution for the Thermal Transport

Using Equations (30) and (31), we can write the system of n-equation describing the Laplace transformed temperature profile:
(65)∂2T¯i∂y2=QisβiT¯i,i∈In1,where Qi=ai/biPri. The general solution of Equation (65) is
(66)T¯iy,s=Ci1se−yQisβi+Ci2seyQisβi,i∈In1,where the unknown parameters C11,C21,⋯,Cn1,C12,C22,⋯,Cn2 are to be computed by using Equations (33) and (36) and are given by the following linear system:
(67)AsYs=Bs,where the As matrix is
(68)As=a000⋯0a1−c10⋯00a2−c2⋯0⋮⋮⋱⋱⋮00⋯an−1−cn−1c000⋯0c1−d10⋯00c2−d2⋯0⋮⋮⋱⋱⋮00⋯cn−1−dn−1−p1q10⋯00−p2q2⋯0⋮⋮⋱⋱⋮00⋯−pn−1qn−100⋯0qnr1−s10⋯00r2−s2⋯0⋮⋮⋱⋱⋮00⋯rn−1−sn−100⋯0sn,and the vectors Ys and Bs are
(69)Ys=C11C21⋮Cn1C12C22⋮Cn2,Bs=g¯1s00⋮0g¯ns1×2n.Here(70)a0=1+δ0Q1sβ1,ai=1−δiQisβie−hiQisβi,c0=1−δ0Q1sβ1,ci=1+δiQisβiehiQisβi,ci=1+δiQi+1sβi+1e−hiQi+1sβi+1,di=1−δiQi+1sβi+1ehiQi+1sβi+1,pi=liQisβie−hiQisβi,ri=liQisβiehiQisβi,li=kiki+1,qi=Qi+1sβi+1e−hiQi+1sβi+1,si=Qi+1sβi+1ehiQi+1sβi+1,i∈In−11,qn=1−δnQnsβne−Qnsβn,sn=1+δnQnsβneQnsβn.We incorporate the following notations to describe the linear system (67) in an appropriate format:
(71)M=a000⋯0a1−c10⋯00a2−c2⋯0⋮⋮⋱⋱⋮00⋯an−1−cn−1n×n;i.e.M=Miji,j∈In−10,M0j=a0δ0,j,j∈In−10,Mij=aiδi,j+1−biδi+1,j+1,i∈In−11,j∈In−10,N=c000⋯0c1−d10⋯00c2−d2⋯0⋮⋮⋱⋱⋮00⋯cn−1−dn−1n×n;i.e.N=Niji,j∈In−10,N0j=c0δ0,j,j∈In−10,Nij=ciδi,j+1−diδi+1,j+1,i∈In−11,j∈In−10,Q=0−p2q2⋯0⋮⋮⋱⋱⋮00⋯−pn−1qn−100⋯0qnn×n;i.e.Q=Qiji,j∈In−10,Qij=−pi+1δi+1,j+1+qi+1δi+2,j+1,i∈In−20,j∈In−10,Qn−1j=qnδn−1,j,j∈In−10,R=r1−s10⋯00r2−s2⋯0⋮⋮⋱⋱⋮00⋯rn−1−sn−100⋯0snn×n;i.e.R=Riji,j∈In−10,Rij=ri+1δi+1,j+1−si+1δi+2,j+1,i∈In−20,j∈In−10,Rn−1j=snδn−1,j,j∈In−10,C1=g¯1s0⋮0n×1,C2=00⋮0g¯2sn×1,A=C11C21⋮Cn1n×1,B=C12C22⋮Cn2n×1,where δij is the Kronecker tensor. The linear system (67) is equivalent to
(72)MA+NB=C1,QA+RB=C2.MatricesM,N,Q,R are invertible and triangular. We can rewrite Equation (72) as
(73)A=S−1N−1C1−R−1C2,B=R−1C2−R−1QS−1N−1C1−R−1C2.,where we suppose that S=N−1M−R−1Q is invertible. An easy computation shows that
(74)R0=R−1C2=g¯2ss1s2⋯sn−1r1r2⋯rn−1sng¯2ss2⋯sn−1r2⋯rn−1sn⋮g¯2ssn−1rn−1sng¯2ssn,C0≔N−1C1−R−1C2=g¯1sa0−g¯2ss1s2⋯sn−1r1r2⋯rn−1sng¯1sb1a0a1−g¯2ss2⋯sn−1r2⋯rn−1sng¯1sb1b2a0a1a2−g¯2ss3⋯sn−1r3⋯rn−1sn⋮g¯1sb1b2⋯bn−2a0a1a2⋯an−2−g¯2ssn−1rn−1sng¯1sb1b2⋯bn−1a0a1a2⋯an−1−g¯2ssn,where
(75)N−1=1a000⋯0c1a0a1−c10⋯0c1c2a0a1a2−c1c2a2−c2⋯0⋮⋮⋱⋱⋮c1c2⋯cn−1a0a1a2⋯an−1−c1c2⋯cn−1a2⋯an−1⋯−cn−2cn−1an−1−cn−1n×n,that is
(76)N−1=niji,j∈In−10,n0j=1a0δ0,j,j∈In−10,ni0=∏ℓ=1icℓa0aℓ,i∈In−11,nij=−∏m=1icmam∏ℓ=0j−1aℓ+1cℓ∑k=0n−1δi,j+k,i,j∈In−11,R−1=1r1s1r1r2s1s2r1r2r3⋯s1s2⋯sn−2r1r2r3⋯rn−1s1s2⋯sn−1r1r2r3⋯rn−1sn01r2s2r2r3⋯s2⋯sn−2r2r3⋯rn−1s2⋯sn−1r2r3⋯rn−1sn⋮⋮⋮⋱⋮⋮000⋯01snn×n,that is
(77)R−1=ρiji,j∈In−10,ρij=∑k=0n−1∏m=0irmsm∏ℓ=0jsℓrℓ+1δi=k,j,r0=s0=rn=1,i,j∈In−10.Now, the matrixS=Siji,j∈In−10 is defined by the elements
(78)Sij=∑k=0n−1nikMkj−ρikQkj.This reduces the system (73) to
(79)A=S−1C0,B=R0−R−1QS−1C0.Now, we know the auxiliary functionsCj1,Cj2, j∈In1; from the linear system (79), the analytical expressions for the solution of temperature profiles Tiy,t can be obtained from Equation (66), by using the inverse Laplace transform as
(80)Tiy,t=12πι∫σ−i∞σ+i∞estCi1se−yQisβi+Ci2seyQisβids,i∈In1.Since the auxiliary functionsCj1,Cj2, j∈In1, involved in Equation (79) are intricate, we therefore have used the following numerical Talbot’s algorithms [45, 46] for the computation of the inverse Laplace transform.Consider the functiongy,t has the Laplace transform G¯y,s. The Talbot algorithm [45] approximates the function gy,t as
(81)gy,t≅rMexprt2G¯y,r+∑k=1M−1ReexptzεkG¯y,zεk1+iζεk,where
(82)r=2M5t,zε=rεcotε+i,ε∈−π,π,ζε=ε+εcotε−1cotε,εk=kπM.The functiongy,t can be approximated by another method, the improved Talbot algorithm [46].
(83)gy,t≅1t∑k=1Mexptz1σkG¯y,z1σkν+iζ1σk,where
(84)z1ε=Mtνiε+μεcotαε−ξ,ε∈−π,π,ζ1ε=αμε+μαεcotαε−1cotαε,σk=2k−1πM−π.Here,α,M,ν,μ,ξ are variables the user must define.


## 3.1. Analytical Solutions for Velocity and Shear Stress

Combining Equations (28) and (29), we can write
(38)ais1+λisαiu¯iy,s=1+λisαiP¯s+bi∂2u¯iy,s∂y2,i∈In1,where φ¯y,s=∫0∞φy,texp−stdt represents the Laplace transform of the function φy,t. In case of nonhomogeneous Robin boundary and fluid-fluid interface conditions (32) and (34), the finite Fourier sine-cosine transform of the function u¯iy,s, i∈In1, can be defined with the help of the Fourier series theory and the Sturm-Liouville theory as [39, 40].
(39)u¯^iηmi,s=1hi−hi−1∫hi−1hiu¯iy,sΦmiydy,where Φmiy=sinηmiy−hi−1/hi−hi−1+li−1/hi−hi−1ηmicosηmiy−hi−1/hi−hi−1 and ηmi are the positive roots of the transcendental equation tanηmi=hi−hi−1li−1+liηmi/li−1liηmi2−hi−hi−12. The inverse Fourier sine-cosine transform of u¯^iηmi,s is defined by
(40)u¯iy,s=∑m=1∞u¯^iηmi,sJmiΦmiy,where
(41)Jmi=1hi−hi−1∫hi−1hiΦmi2ydy=li−12ηmi2−hi−hi−124ηmihi−hi−12sin2ηmi+hi−hi−12+li−1hi−hi−1+li−12ηmi2−li−1hi−hi−1cos2ηmi2hi−hi−12.By direct computations, using the robin boundary condition (32) and interface fluid-fluid condition (34), we can write
(42)∂2u¯^i∂y2ηmi,s=ηmihi−hi−12cosηmihi−hi−12+li−12ηmi2li−1liηmi2−hi−hi−12w¯ihi,li,s+v¯ihi−1,li−1,s−ηmi2hi−hi−12u¯^iηmi,s,i∈In1,where we rename v¯ihi−1,li−1,s=u¯ihi−1,s−li−1∂u¯ihi−1,s/∂y, w¯ihi,li,s=u¯ihi,s+li∂u¯ihi,s/∂y, i∈In−11,v¯0h0,l0,s=f¯1s, and w¯nhn,ln,s=f¯2s. With the application of finite Fourier sine-cosine transform (39) to Equation (38) corresponding to the robin boundary conditions (32) and interface fluid-fluid conditions (34), and using Equation (42), the transformed velocities take the form
(43)u¯^iηmi,s=hi−hi−121+λisαiP¯saihi−hi−12s1+λisαi+biηmi2+biηmiaihi−hi−12s1+λisαi+biηmi2×cosηmihi−hi−12+li−12ηmi2li−1liηmi2−hi−hi−12w¯ihi,li,s+v¯ihi−1,li−1,s,i∈In1.To apply the inverse Fourier sine-cosine transform, we rewrite Equation (44) in the following suitable form:
(44)u¯^iηmi,s=li−1ηmisinηmi+hi−hi−11−cosηmiηmihi−hi−1hi−hi−1+liv¯ihi−1,li−1,s+li−1w¯ihi,li,shi−hi−1+li−1+li+−li−1ηmi−hi−hi−1−li−1ηmicosηmi+hi−hi−1+li−1ηmi2sinηmihi−hi−1ηmi2w¯ihi,li,s−v¯ihi−1,li−1,shi−hi−1hi−hi−1+li−1+li+biηmiaihi−hi−12s1+λisαi+biηmi2+C1miv¯ihi−1,li−1,s+biηmicosηmiaihi−hi−12s1+λisαi+biηmi2hi−hi−12+li−12ηmi2li−1liηmi2−hi−hi−12+C2miw¯ihi,li,s+hi−hi−121+λisαiP¯saihi−hi−12s1+λisαi+biηmi2,i∈In1,where
(45)C1mi=−1ηmi+hi−hi−1li−1+liηmicosηmi+hi−hi−12sinηmi−lili−1ηmi2sinηmihi−hi−1hi−hi−1+li−1+liηmi2,C2mi=hi−hi−12ηmicosηmi−hi−hi−12+li−1hi−hi−1+li−1ηmi2sinηmihi−hi−1hi−hi−1+li−1+liηmi2.Considering the auxiliary relations along with their Fourier sine-cosine transformv(46)φ1iy=1,y∈hi−1,hi,φ^1iηmi=li−1ηmisinηmi+hi−hi−11−cosηmiηmihi−hi−1,(47)φ2iy=y−hi−1hi−hi−1,y∈hi−1,hi,(48)φ^2iηmi=−li−1ηmi−hi−hi−1−li−1ηmicosηmi+hi−hi−1+li−1ηmi2sinηmihi−hi−1ηmi2,(49)m=1,2,⋯,i∈In1,h0=0,hn=1,the inverse Fourier sine-cosine transform of Equation (46) takes the form
(50)u¯iy,s=hi−hi−1+liv¯ihi−1,li−1,s+li−1w¯ihi,li,shi−hi−1+li−1+li+w¯ihi,li,s−v¯ihi−1,li−1,shi−hi−1+li−1+liy−hi−1+∑m=1∞biηmiaihi−hi−12s1+λisαi+biηmi2+C1miΦmiyv¯ihi−1,li−1,sJmi+∑m=1∞biηmicosηmiaihi−hi−12s1+λisαi+biηmi2hi−hi−12+li−12ηmi2li−1liηmi2−hi−hi−12+C2miΦmiyw¯ihi,li,sJmi+∑m=1∞hi−hi−121+λisαiaihi−hi−12s1+λisαi+biηmi2ΦmiyP¯sJmi,i∈In1.Now, from (29) and (50), we obtain
(51)τ¯iy,s=T¯i,1y,sw¯ihi,li,s−T¯i,2y,sv¯ihi−1,li−1,s+T¯i,3y,sP¯s,i∈In1,where
(52)T¯i,1y,s=bi1+λisαihi−hi−1+li−1+li+bihi−hi−11+λisαi∑m=1∞ηmibiηmicosηmiaihi−hi−12s1+λisαi+biηmi2hi−hi−12+li−12ηmi2li−1liηmi2−hi−hi−12+C2mi×cosηmiy−hi−1/hi−hi−1−li−1/hi−hi−1ηmisinηmiy−hi−1/hi−hi−1Jmi,T¯i,2y,s=bi1+λisαihi−hi−1+li−1+li−bihi−hi−11+λisαi∑m=1∞ηmibiηmiaihi−hi−12s1+λisαi+biηmi2+C1mi×cosηmiy−hi−1/hi−hi−1−li−1/hi−hi−1ηmisinηmiy−hi−1/hi−hi−1Jmi,T¯i,3y,s=bihi−hi−1P¯s∑m=1∞ηmicosηmiy−hi−1/hi−hi−1−li−1/hi−hi−1ηmisinηmiy−hi−1/hi−hi−1Jmiaihi−hi−12s1+λisαi+biηmi2.Using Equations (32) and (34) in Equation (51), we get the following linear system (for detail, see “Appendix”):
(53)M¯sW¯s=N¯s,where
(54)M¯s=Y¯1,1−Y¯1,200⋯0Y¯2,1Y¯2,2−Y¯2,30⋯00Y¯3,2Y¯3,3−Y¯3,4⋯0⋮⋮⋱⋱⋱⋮00⋯Y¯n−3,n−4Y¯n−3,n−3−Y¯n−3,n−200⋯Y¯n−2,n−3Y¯n−2,n−2−Y¯n−2,n−100⋯0Y¯n−1,n−2Y¯n−1,n−1,W¯s=w¯1h1,l1,sw¯2h2,l2,s⋮w¯n−1hn−1,ln−1,s,N¯s=b¯1b¯2⋮b¯n−1.Finally, we obtain(55)W¯s=M¯−1sN¯s.Now,w¯ihi,li,s=v¯i+1hi,li,s,i∈In−11, are known functions; therefore, the velocities u¯1y,s,⋯,u¯ny,s are known. In order to obtain the inverse Laplace transforms of the functions u¯iy,s,i∈In1, we consider the following auxiliary functions:
(56)H¯i0m,s=1aihi−hi−12s1+λisαi+ηmi2=s−1aihi−hi−12λisαi+λi−1·11+ηmi2s−1/aihi−hi−12λisαi+λi−1=1aiλihi−hi−12∑k=0∞−1kηmi2aihi−hi−12λiks−k−1sαi+λi−1k+1,H¯i1m,s=1+λisαihi−hi−12ais1+λisαi+biϑmi2=1aihi−hi−12∑k=0∞−1kηmi2hi−hi−12aiλiks−k−1sαi+λi−1k,H¯i4m,s=11+λisαiaihi−hi−12s1+λisαi+biϑmi2=1aiλi2hi−hi−12∑k=0∞−1kηmi2aihi−hi−12λiks−k−1sαi+λi−1k+2.Since the generalized G-Lorenzo-Hartley function is defined by [41]
(57)Gσ1,σ2,σ3t,σ=L−1sσ2sσ1−σσ3=∑k=0∞Γk+σ3σktk+σ3σ1−σ2−1k!Γk+σ3σ1−σ2Γσ3,Res>0,Reσ1σ3−σ2>0,σsσ1<1,and for αi,βi>0,(58)L−1sαi−βisαi−d=tβi−1Eαi,βidtαi,where Eαi,βi. is the Mittag-Leffler function [42]. The inverse Laplace transform of H¯i0m,s, H¯i1m,s, and H¯i4m,s takes the form
(59)Hi0m,t=1aiλihi−hi−12∑k=0∞−1kηmi2aihi−hi−12λikGαi,−k−1,k+1t,−λi−1,(60)Hi1m,t=1hi−hi−12ai∑k=0∞−1kηmi2khi−hi−12aiλik+1Gαi,−k−1,kt,−λi−1,(61)Hi4m,t=1aiλi2hi−hi−12∑k=0∞−1kηmi2hi−hi−12aiλikGαi,−k−1,k+2t,−λi−1,where δt is the Dirac delta function. Using Equations (59), (60), and (50), we obtain for velocities uiy,t,i∈In1, with the following expressions:
(62)uiy,t=hi−hi−1+liv¯ihi−1,li−1,t+li−1w¯ihi,li,thi−hi−1+li−1+li+w¯ihi,li,t−v¯ihi−1,li−1,thi−hi−1+li−1+liy−hi−1+∑m=1∞biηmiHi0m,t∗vihi−1,li−1,t+C1mivihi−1,li−1,tΦmiyJmi+∑m=1∞hi−hi−12+li−12ηmi2biηmicosηmili−1liηmi2−hi−hi−12Hi0m,t∗wihi,li,t+C2miwihi,li,tΦmiyJmi+∑m=1∞hi−hi−12Hi1m,t∗PtΦmiyJmi,i∈In1,where h1t∗h2t=∫0th1t−τh2τdτ is the convolution product of the functions h1t and h2t. The system of shear stresses τiy,t, i∈In1, can be determined by applying inverse Laplace transform to Equation (51) and using Equations (59) and (61).
(63)τiy,t=Ti1y,t∗wihi,li,t−Ti2y,t∗vihi−1,li−1,t+Ti3y,t∗Pt,where
(64)Ti,1y,t=biλi−1tαi−1Eαi,αi−λi−1tαihi−hi−1+li−1+li+bihi−hi−1∑m=1∞ηmibiηmicosηmihi−hi−12+li−12ηmi2li−1liηmi2−hi−hi−12Hi4m,t+C2miλi−1tαi−1Eαi,αi−λi−1tαicosηmiy−hi−1/hi−hi−1−li−1/hi−hi−1ηmisinηmiy−hi−1/hi−hi−1Jmi,T¯i,2y,s=biλi−1tαi−1Eαi,αi−λi−1tαihi−hi−1+li−1+li−bihi−hi−1∑m=1∞ηmibiηmiHi4m,t+C1miλi−1tαi−1Eαi,αi−λi−1tαicosηmiy−hi−1/hi−hi−1−li−1/hi−hi−1ηmisinηmiy−hi−1/hi−hi−1Jmi,T¯i,3y,s=bi∑m=1∞ηmihi−hi−1Hi0m,t×cosηmiy−hi−1/hi−hi−1−li−1/hi−hi−1ηmisinηmiy−hi−1/hi−hi−1Jmi.At the end of this section, we mention that, in the case of a single fluid, there are interesting studies in the literature on the flow of Maxwell fluids with slip on the channel walls. Thus, in the particular caseα=1,n=1,Pt=f2t=0, our results are equivalent with those obtained in ([43], Equations (34) and (55)). Moreover for the Newtonian case with n=1,Pt=0,λ1=0, the velocity profile given by Equation (62) is equivalent to the solution obtained in [44], Equation (65) with φ=0 and ρβf=0.


## 3.2. Solution for the Thermal Transport

Using Equations (30) and (31), we can write the system of n-equation describing the Laplace transformed temperature profile:
(65)∂2T¯i∂y2=QisβiT¯i,i∈In1,where Qi=ai/biPri. The general solution of Equation (65) is
(66)T¯iy,s=Ci1se−yQisβi+Ci2seyQisβi,i∈In1,where the unknown parameters C11,C21,⋯,Cn1,C12,C22,⋯,Cn2 are to be computed by using Equations (33) and (36) and are given by the following linear system:
(67)AsYs=Bs,where the As matrix is
(68)As=a000⋯0a1−c10⋯00a2−c2⋯0⋮⋮⋱⋱⋮00⋯an−1−cn−1c000⋯0c1−d10⋯00c2−d2⋯0⋮⋮⋱⋱⋮00⋯cn−1−dn−1−p1q10⋯00−p2q2⋯0⋮⋮⋱⋱⋮00⋯−pn−1qn−100⋯0qnr1−s10⋯00r2−s2⋯0⋮⋮⋱⋱⋮00⋯rn−1−sn−100⋯0sn,and the vectors Ys and Bs are
(69)Ys=C11C21⋮Cn1C12C22⋮Cn2,Bs=g¯1s00⋮0g¯ns1×2n.Here(70)a0=1+δ0Q1sβ1,ai=1−δiQisβie−hiQisβi,c0=1−δ0Q1sβ1,ci=1+δiQisβiehiQisβi,ci=1+δiQi+1sβi+1e−hiQi+1sβi+1,di=1−δiQi+1sβi+1ehiQi+1sβi+1,pi=liQisβie−hiQisβi,ri=liQisβiehiQisβi,li=kiki+1,qi=Qi+1sβi+1e−hiQi+1sβi+1,si=Qi+1sβi+1ehiQi+1sβi+1,i∈In−11,qn=1−δnQnsβne−Qnsβn,sn=1+δnQnsβneQnsβn.We incorporate the following notations to describe the linear system (67) in an appropriate format:
(71)M=a000⋯0a1−c10⋯00a2−c2⋯0⋮⋮⋱⋱⋮00⋯an−1−cn−1n×n;i.e.M=Miji,j∈In−10,M0j=a0δ0,j,j∈In−10,Mij=aiδi,j+1−biδi+1,j+1,i∈In−11,j∈In−10,N=c000⋯0c1−d10⋯00c2−d2⋯0⋮⋮⋱⋱⋮00⋯cn−1−dn−1n×n;i.e.N=Niji,j∈In−10,N0j=c0δ0,j,j∈In−10,Nij=ciδi,j+1−diδi+1,j+1,i∈In−11,j∈In−10,Q=0−p2q2⋯0⋮⋮⋱⋱⋮00⋯−pn−1qn−100⋯0qnn×n;i.e.Q=Qiji,j∈In−10,Qij=−pi+1δi+1,j+1+qi+1δi+2,j+1,i∈In−20,j∈In−10,Qn−1j=qnδn−1,j,j∈In−10,R=r1−s10⋯00r2−s2⋯0⋮⋮⋱⋱⋮00⋯rn−1−sn−100⋯0snn×n;i.e.R=Riji,j∈In−10,Rij=ri+1δi+1,j+1−si+1δi+2,j+1,i∈In−20,j∈In−10,Rn−1j=snδn−1,j,j∈In−10,C1=g¯1s0⋮0n×1,C2=00⋮0g¯2sn×1,A=C11C21⋮Cn1n×1,B=C12C22⋮Cn2n×1,where δij is the Kronecker tensor. The linear system (67) is equivalent to
(72)MA+NB=C1,QA+RB=C2.MatricesM,N,Q,R are invertible and triangular. We can rewrite Equation (72) as
(73)A=S−1N−1C1−R−1C2,B=R−1C2−R−1QS−1N−1C1−R−1C2.,where we suppose that S=N−1M−R−1Q is invertible. An easy computation shows that
(74)R0=R−1C2=g¯2ss1s2⋯sn−1r1r2⋯rn−1sng¯2ss2⋯sn−1r2⋯rn−1sn⋮g¯2ssn−1rn−1sng¯2ssn,C0≔N−1C1−R−1C2=g¯1sa0−g¯2ss1s2⋯sn−1r1r2⋯rn−1sng¯1sb1a0a1−g¯2ss2⋯sn−1r2⋯rn−1sng¯1sb1b2a0a1a2−g¯2ss3⋯sn−1r3⋯rn−1sn⋮g¯1sb1b2⋯bn−2a0a1a2⋯an−2−g¯2ssn−1rn−1sng¯1sb1b2⋯bn−1a0a1a2⋯an−1−g¯2ssn,where
(75)N−1=1a000⋯0c1a0a1−c10⋯0c1c2a0a1a2−c1c2a2−c2⋯0⋮⋮⋱⋱⋮c1c2⋯cn−1a0a1a2⋯an−1−c1c2⋯cn−1a2⋯an−1⋯−cn−2cn−1an−1−cn−1n×n,that is
(76)N−1=niji,j∈In−10,n0j=1a0δ0,j,j∈In−10,ni0=∏ℓ=1icℓa0aℓ,i∈In−11,nij=−∏m=1icmam∏ℓ=0j−1aℓ+1cℓ∑k=0n−1δi,j+k,i,j∈In−11,R−1=1r1s1r1r2s1s2r1r2r3⋯s1s2⋯sn−2r1r2r3⋯rn−1s1s2⋯sn−1r1r2r3⋯rn−1sn01r2s2r2r3⋯s2⋯sn−2r2r3⋯rn−1s2⋯sn−1r2r3⋯rn−1sn⋮⋮⋮⋱⋮⋮000⋯01snn×n,that is
(77)R−1=ρiji,j∈In−10,ρij=∑k=0n−1∏m=0irmsm∏ℓ=0jsℓrℓ+1δi=k,j,r0=s0=rn=1,i,j∈In−10.Now, the matrixS=Siji,j∈In−10 is defined by the elements
(78)Sij=∑k=0n−1nikMkj−ρikQkj.This reduces the system (73) to
(79)A=S−1C0,B=R0−R−1QS−1C0.Now, we know the auxiliary functionsCj1,Cj2, j∈In1; from the linear system (79), the analytical expressions for the solution of temperature profiles Tiy,t can be obtained from Equation (66), by using the inverse Laplace transform as
(80)Tiy,t=12πι∫σ−i∞σ+i∞estCi1se−yQisβi+Ci2seyQisβids,i∈In1.Since the auxiliary functionsCj1,Cj2, j∈In1, involved in Equation (79) are intricate, we therefore have used the following numerical Talbot’s algorithms [45, 46] for the computation of the inverse Laplace transform.Consider the functiongy,t has the Laplace transform G¯y,s. The Talbot algorithm [45] approximates the function gy,t as
(81)gy,t≅rMexprt2G¯y,r+∑k=1M−1ReexptzεkG¯y,zεk1+iζεk,where
(82)r=2M5t,zε=rεcotε+i,ε∈−π,π,ζε=ε+εcotε−1cotε,εk=kπM.The functiongy,t can be approximated by another method, the improved Talbot algorithm [46].
(83)gy,t≅1t∑k=1Mexptz1σkG¯y,z1σkν+iζ1σk,where
(84)z1ε=Mtνiε+μεcotαε−ξ,ε∈−π,π,ζ1ε=αμε+μαεcotαε−1cotαε,σk=2k−1πM−π.Here,α,M,ν,μ,ξ are variables the user must define.


## 4. Numerical Results and Discussions

The unsteady, laminar flow with thermal conductivity of the simultaneous n-layer fractional immiscible Maxwell fluids in a rectangular channel has been examined. The motion of these fluids is produced by the time-based pressure gradient in the flow direction and by the displacement of the channel boundaries with the time-based fluid-fluid interfacial conditions.In this problem, the generalization puts into consideration the fractional constitutive equation of the Maxwell fluids based on the Caputo time-fractional derivative; therefore, the velocity gradient histories influence the fluid behavior. Such type of flow is so-called the flow with memory.On the solid boundaries, the Robin boundary conditions are taken into account, whereas the velocity and shear stress are assumed continuous at the fluid-fluid interfacey=d0. Semianalytical results of the velocities, shear stresses, and temperature profiles are determined with the help of the Laplace transformation and the Talbot algorithms used for the numerical inverse Laplace transforms. Moreover, an analytical solution for the flow is recovered by using the Laplace transform in conjunction with the finite sine-Fourier transform.The findings obtained are generic; therefore, several special cases can be considered. Multilayer flows of ordinary/fractional Maxwell fluids can be analyzed as special cases by allowing certain fractional parameters to be equal to one.With the help of the Mathcad software, numerical results have been illustrated for the obtained solutions of fluid velocities and temperatures. These results are shown in Figures2–14. For graphical illustration of the fluid velocities, we have used the material parameters
(85)ρ1=1000,μ1=0.05,G1=1.2,cp1=0.25,k1=2.5,l1=0.1,σ0=0.05,ρ2=1300,μ2=0.2,G2=2.4,cp2=0.27,k2=2.7,l2=0.15,σ1=0.1,ρ3=1500,μ3=0.6,G3=5.0,cp3=0.29,k3=2.9,l3=0.2,σ2=0.15,σ3=0.2.Figure 2
Temperature profilesTiy,t, i=1,2,3 versus y for fractional multilayer Maxwell fluid at t=0.002,β2=0.6, β3=0.8, and different values of β1.Figure 3
Temperature profilesTiy,t, i=1,2,3 versus y for fractional multilayer Maxwell fluid at t=0.002,β1=0.2, β3=0.8, and different values of β2.Figure 4
Temperature profilesTiy,t, i=1,2,3 versus y for fractional multilayer Maxwell fluid at t=0.002,β1=0.2, β2=0.4, and different values of β3.Figure 5
Temperature profilesTiy,t, i=1,2,3 versus t for fractional multilayer Maxwell fluid at β1=0.0.4, β2=0.6, and different values of y.Figure 6
Temperature profilesTiy,t, i=1,2,3 for small values of the time t (β1=0.4,β2=0.6,andβ3=0.8).Figure 7
Temperature profilesTiy,t, i=1,2,3 for large values of the time t (β1=0.4,β2=0.6,andβ3=0.8).Figure 8
Variation of the kernelstβi/Γβi with the fractional parameters βi, i=1,2,3.Figure 9
Velocity profilesuiy,t, i=1,2,3 versus y for fractional multilayer Maxwell fluid at α1=0.2,α2=0.4,α3=0.8,λ1=0.021,λ2=0.042,λ3=0.06, and different values of time t.Figure 10
Velocity profilesuiy,t, i=1,2,3 versus y for fractional multilayer Maxwell fluid at t=0.19,λ1=0.021,λ2=0.042,λ3=0.06, and different values of fractional parameters α1=α2.Figure 11
Velocity profilesuiy,t, i=1,2,3 versus y for fractional multilayer Maxwell fluid at t=0.19,λ1=0.021,λ2=0.042,λ3=0.06, and different values of fractional parameters α2=α3.Figure 12
Velocity profilesuiy,t, i=1,2,3 versus y for fractional multilayer Maxwell fluid at t=0.19,λ1=0.021,λ2=0.042,λ3=0.06, and different values of fractional parameters α1=α3.Figure 13
The profiles of velocityuiy,t, i=1,2,3, for small values of the time t (α1=0.4,α2=0.6,α3=0.8a1=1,a2=1.3,a3=1.5,b1=1,b2=4,b3=12, h1=0.2, h2=0.6, and λ1=0.021,λ2=0.042,λ3=0.060).Figure 14
The profiles of velocityuiy,t, i=1,2,3 for large values of the time t (α1=0.4,α2=0.6,α3=0.8a1=1,a2=1.3,a3=1.5,b1=1,b2=4,b3=12, h1=0.2, h2=0.6, and λ1=0.021,λ2=0.042,λ3=0.060).To analyze the temperature profiles of the three-layer fluids withh0=0,h1=0.3,h2=0.7,h3=1, we take the special instance when the functions g1t,g2t appearing in the Robin boundary conditions are constant, i.e., g1t=0.6Ht, g2t=1.2Ht, where Ht is the Heaviside function. It is shown in Figure 2 that the variance of the fractional temperature variable β1 of the fluid located in the first layer has a noticeable effect on the thermal profile in the first two layers; however, the effect on the temperature from the topmost layer is minor. This has been caused by the decrease in thermal profile in the first two layers with the increase in fractional variable β1. It is revealed in Figure 2 that the variance of the fractional temperature variable β2 of the fluid located in the second layer has a noticeable effect on the thermal profile in the last two layers; however, the effect on the temperature from the first layer is negligible. This has been caused by the decrease in thermal profile in the last two layers with the increase in fractional variable β1. It is exposed in Figure 4 that the variance of the fractional temperature variable β3 of the fluid located in the third layer has a noticeable effect on the thermal profile in the first and the last layers; however, the effect on the temperature from the middle layer is minor. This has been caused by the decrease in thermal profile in the first and the last layers with the increase in fractional variable β3.It is recognized by Equation (26) that the heat flux memory kernel is the relation hβit=tβi−1/Γβi whose plot is shown in Figure 8. It is described in Figure 8 that, for βi≥0.2 and t=0.002, the kernel hβit, i=1,2,3, declines with βi; therefore, the temperature gradient distribution relation declines (memory impacts are softer). Time-evolution of the thermal layers, in various channel locations and for differing values of the temperature fractional variables, is shown in Figure 5. As shown in Figure 5, the thermal profile as a function of the fractional variables 5 has different characteristics for low, versus high, time values. This attitude reflects prior studies shown in Figures 2–4 and also with the progression of the heat memory kernel, hβit, i=1,2,3. We have to mention that the temperatures of all layers vary for small time fluctuations, but after a low time valuation for t (t=1 in examined scenarios), the temperatures are just about unchanged (see Figures 6 and 7).In order to examine the three layers, fluid velocities withh0=0,h1=0.2,h2=0.6,h3=1, we considered the particular case when the functions in the slip boundary conditions are constant, i.e., f1t=0.5Ht, f2t=0.8Ht, where Ht is the Heaviside function and the considered pressure gradient is Pt=sint. Figure 9 illustrates the time effect on the velocity profile uiy,t, i=1,2,3. The velocity profile is observed to increase with the increase in time. Figures 10–12 are plotted to study the influence of the velocity fractional parameters αi on the velocity fields. It is observed that the velocity profile is increasing with the increase in the fractional parameters αi1=αi2 and for fixed αi3,i1,i2,i3=1,2,3. The fractional variables have braking impacts. The first two layers’ flows are accelerated; however, the last layer’s fluid is slowed down.In Figures13 and 14, the profiles of velocity uiy,t, i=1,2,3, for different values of the time t, are presented. It is observed that, for t≥40, the profiles of velocity are unchanged; therefore, the velocity is given by the permanent solution. It can be observed from Figures 6 and 13 that at a very short time, the profiles for velocities and temperatures are similar to the initial conditions, that is, zero everywhere.


## 5. Conclusions

Time-dependent simultaneous n-layer fluid flow in a rectangular channel was examined through heat exchange of Maxwell immiscible fluids with generic constitutive equations for the shear stress and heat flux.The Caputo time-fractional derivative defines the generic constitutive relations; thus, the behavior of the fluid is determined by the histories of the temperature and velocity gradient.We have used the finite Fourier and Laplace transform coupled with numerical Laplace inversions for the analytical and semianalytical results for the velocity, temperature, and shear stress profiles with the assumption that for the adjacent layers, the interfacial heat fluxes and shear stresses are equal and in the presence of the interface Robin-type interfacial conditions.The results of this study attained are of a general nature; however, many specific situations can be produced. Multilayer flows of ordinary/fractional Maxwell fluids can be analyzed as special cases by allowing certain fractional parameters to be equal to one.With the help of the Mathcad software, numerical results have been illustrated for the obtained solutions of fluid velocities and temperatures.For small and, respectively, large values of timet, the fluid flow and the heat transfer differ. Such specific characteristics are due to the differences in time and fractional parameters of the thermal/velocity kernels, so the memory influences have a tremendous impact on the fluids.



---

*Source: 5572823-2021-04-14.xml*