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.

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 [14]. 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 [68]. 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 [1416]. 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 [2638] 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 with the boundary plates positioned at and . To begin with at time , the two walls and the fluids inside are at rest with the atmospheric temperature After this instance, the boundary plate positioned at move with velocity along the -axis and the wall temperature , while the channel plate moves with the wall velocity analogous with the -axis and the wall temperature (Figure 1). It has been assumed that the functions , , and are piece-wise continuous functions, and . We consider that the velocity vector is of the form . We assume the simultaneous n-immiscible Maxwell fluids between two parallel boundary planes. It divides the domain of flow into subdomains . In the region , Maxwell fluid flow has the viscosity , relaxation time , density , , the elastic modulus, temperature , velocity , the shear stress , and the thermal flux where .

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 velocities , , where , and we assume that 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: (ii)The n-constitutive equations: (iii)The corresponding n-initial conditions: (iv)The boundary condition: (v)At the interface, we assume the continuity of Robin conditions:

The system of constitutive relations for temperature and heat flux is given by the following: (i)The n-thermal transport balance equations:(ii)The n-thermal fluxes with Fourier’s law:where and are the thermal conductivity and the specific heat at constant pressure, respectively. We consider the following: (i)n-Initial conditions: (ii)Boundary conditions: (iii)Interface conditions:

Consider the nondimensional variables

The dimensionless form of the governing Equations (2)–(11) are reduced into the following form: along with the following: (i)The dimensionless n-initial conditions: (ii)The dimensionless boundary condition: (iii)The nondimensional interface conditions: where and is the characteristic scale.

2.1. Generalized Constitutive Mathematical Model

Consider the following generalized constitutive mathematical relations: where the Caputo derivative is defined as

For the special case . It is significant to note that the fractional models (22) and (23) have the following equivalent formulae as

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 , while the nonlocality kernel for the shear stress is the function being the Mittag-Leffler function. In this study, we consider the pressure gradient in the flow direction as a known function, namely where is a peice-wise continuous function on .

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

3.1. Analytical Solutions for Velocity and Shear Stress

Combining Equations (28) and (29), we can write where represents the Laplace transform of the function . In case of nonhomogeneous Robin boundary and fluid-fluid interface conditions (32) and (34), the finite Fourier sine-cosine transform of the function , , can be defined with the help of the Fourier series theory and the Sturm-Liouville theory as [39, 40]. where and are the positive roots of the transcendental equation . The inverse Fourier sine-cosine transform of is defined by where

By direct computations, using the robin boundary condition (32) and interface fluid-fluid condition (34), we can write where we rename , , , and 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

To apply the inverse Fourier sine-cosine transform, we rewrite Equation (44) in the following suitable form: where

Considering the auxiliary relations along with their Fourier sine-cosine transformv the inverse Fourier sine-cosine transform of Equation (46) takes the form

Now, from (29) and (50), we obtain where

Using Equations (32) and (34) in Equation (51), we get the following linear system (for detail, see “Appendix”): where

Finally, we obtain

Now, , are known functions; therefore, the velocities are known. In order to obtain the inverse Laplace transforms of the functions , we consider the following auxiliary functions:

Since the generalized G-Lorenzo-Hartley function is defined by [41] and for where is the Mittag-Leffler function [42]. The inverse Laplace transform of , , and takes the form where is the Dirac delta function. Using Equations (59), (60), and (50), we obtain for velocities with the following expressions: where is the convolution product of the functions and . The system of shear stresses , can be determined by applying inverse Laplace transform to Equation (51) and using Equations (59) and (61). where

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 , our results are equivalent with those obtained in ([43], Equations (34) and (55)). Moreover for the Newtonian case with , the velocity profile given by Equation (62) is equivalent to the solution obtained in [44], Equation (65) with and