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.

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 .

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: where The general solution of Equation (65) is where the unknown parameters are to be computed by using Equations (33) and (36) and are given by the following linear system: where the matrix is and the vectors and are

Here

We incorporate the following notations to describe the linear system (67) in an appropriate format: where is the Kronecker tensor. The linear system (67) is equivalent to

Matrices are invertible and triangular. We can rewrite Equation (72) as where we suppose that is invertible. An easy computation shows that where that is that is

Now, the matrix is defined by the elements

This reduces the system (73) to

Now, we know the auxiliary functions , ; from the linear system (79), the analytical expressions for the solution of temperature profiles can be obtained from Equation (66), by using the inverse Laplace transform as

Since the auxiliary functions , , 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 function has the Laplace transform . The Talbot algorithm [45] approximates the function as where

The function can be approximated by another method, the improved Talbot algorithm [46]. where

Here, 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 interface . 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 Figures 214. For graphical illustration of the fluid velocities, we have used the material parameters

To analyze the temperature profiles of the three-layer fluids with , we take the special instance when the functions appearing in the Robin boundary conditions are constant, i.e., , , where is the Heaviside function. It is shown in Figure 2 that the variance of the fractional temperature variable 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 . It is revealed in Figure 2 that the variance of the fractional temperature variable 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 . It is exposed in Figure 4 that the variance of the fractional temperature variable 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 .

It is recognized by Equation (26) that the heat flux memory kernel is the relation whose plot is shown in Figure 8. It is described in Figure 8 that, for and , the kernel , declines with ; 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 24 and also with the progression of the heat memory kernel, , . We have to mention that the temperatures of all layers vary for small time fluctuations, but after a low time valuation for ( in examined scenarios), the temperatures are just about unchanged (see Figures 6 and 7).

In order to examine the three layers, fluid velocities with , we considered the particular case when the functions in the slip boundary conditions are constant, i.e., , , where is the Heaviside function and the considered pressure gradient is . Figure 9 illustrates the time effect on the velocity profile , The velocity profile is observed to increase with the increase in time. Figures 1012 are plotted to study the influence of the velocity fractional parameters on the velocity fields. It is observed that the velocity profile is increasing with the increase in the fractional parameters and for fixed The fractional variables have braking impacts. The first two layers’ flows are accelerated; however, the last layer’s fluid is slowed down.

In Figures 13 and 14, the profiles of velocity , , for different values of the time , are presented. It is observed that, for , 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 time , 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.

Appendix

Using the boundary conditions (32) and the interface conditions (34) in Equation. (63), we obtained the following algebraic system: where

is the Kronecker delta. The system (A.1) can be written in the following equivalent form:

Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest

The authors declare that they have no conflicts of interest.