Dual Stratiﬁed Nanoﬂuid Flow Past a Permeable Shrinking/Stretching Sheet Using a Non-Fourier Energy Model

: The present study emphasizes the combined effects of double stratiﬁcation and buoyancy forces on nanoﬂuid ﬂow past a shrinking/stretching surface. A permeable sheet is used to give way for possible wall ﬂuid suction while the magnetic ﬁeld is imposed normal to the sheet. The governing boundary layer with non-Fourier energy equations (partial differential equations (PDEs)) are converted into a set of nonlinear ordinary differential equations (ODEs) using similarity transformations. The approximate relative error between present results (using the boundary value problem with fourth order accuracy (bvp4c) function) and previous studies in few limiting cases is sufﬁciently small (0% to 0.3694%). Numerical solutions are graphically displayed for several physical parameters namely suction, magnetic, thermal relaxation, thermal and solutal stratiﬁcations on the velocity, temperature and nanoparticles volume fraction proﬁles. The non-Fourier energy equation gives a different estimation of heat and mass transfer rates as compared to the classical energy equation. The heat transfer rate approximately elevates 5.83% to 12.13% when the thermal relaxation parameter is added for both shrinking and stretching cases. Adversely, the mass transfer rate declines within the range of 1.02% to 2.42%. It is also evident in the present work that the augmentation of suitable wall mass suction will generate dual solutions. The existence of two solutions (ﬁrst and second) are noticed in all the proﬁles as well as the local skin friction, Nusselt number and Sherwood number graphs within the considerable range of parameters. The implementation of stability analysis asserts that the ﬁrst solution is the real solution.


Introduction
The boundary layer theory has been a subject of interest for many researchers because of its significance in many industrial and engineering applications, i.e., airfoil design of of airplanes, the automobile industry and friction drag of a ship. The analytical solution on the boundary layer flow of a viscous fluid was introduced by Crane [1] and Miklavčič and Wang [2] for stretching and shrinking cases, respectively. A stretching sheet solution would generate a far field suction towards the sheet but opposite flow behaviors were observed for the shrinking case. The shrinking sheet seemed to distract the velocity away from the sheet. An adequate wall mass suction is needed to induce the flow due to the unconfined vorticity within the boundary layer. Wang [3] analyzed the stagnation point flow over a shrinking sheet and discovered that the solutions were non-unique and did not exist for larger shrinking rates. Fang et al. [4] analytically investigated viscous flow over shrinking rate with a slip condition and also concluded that the imposition of wall mass suction would induce two solutions. Jusoh et al. [5] also concluded that the application of suction at certain values of the involving parameters contributed to the appearance of the second solution. Very recently, Soomro et al. [6] obtained dual solutions for mixed convective flow of a viscous fluid with copper nanoparticles over a permeable shrinking cylinder. The findings showed that the mixed effects of opposing flow and wall mass suction might generate two solutions. However, energy generation becomes a major issue in the industrial requirements. Attention was focused on enhancing the heat transfer rate of the systems, which include power stations, chemical plants, air conditioning and petrochemical industry. Nanofluid is created by dispersing the nanoparticles (size < 100 nm) into a base fluid such as water, ethylene or propylene glycol to inflate the thermal conductivity. Traditional Newtonian fluids such as water and air have been utilized as cooling fluids in many industrial and manufacturing processes but due to the final inspection of the product quality, it does not seem that the viscous fluid is a good choice for the cooling fluid. Al-Waeli et al. [7] reported that the intensification of the working fluid's thermal conductivity can increase the thermal and electrical as well as the heat transfer efficiencies. The invention of nanofluids with great thermophysical properties can improve the performance for massive engineering operations such as in nuclear cooling systems, solar thermal energy, biomedical applications, lubrication, coolant in automobile radiators and others [8][9][10][11][12][13][14][15]. Two remarkable nanofluid models were developed by Buongiorno [16] and Tiwari and Das [17]. Buongiorno [16] highlighted seven slip mechanisms which can contribute relative velocity between the nanoparticles and the working fluids; however, only Brownian motion and thermophoresis are significant to model the nanofluids. Tiwari and Das [17] proposed a model that investigated the effect of nanoparticles volume fraction on the enhancement of the thermal properties. Khan and Pop [18] intiated the investigation on the flow of a nanofluid towards a stretching sheet by using Buongiorno's model. Latterly, Alsarraf et al. [19] explored the nanoparticles shape effect (brick, blade, cylindrical, platelet and spherical) on the boehmite alumina nanofluid flow characteristics. There are also a few boundary layer studies that modeled the non-Newtonian fluid as the operating fluid with the presence of nanoparticles as reported by Mahmood et al. [20], Aziz and Jamshed [21], Jamshed and Aziz [22] and Aziz et al. [23].
Most of the boundary layer flow and heat transfer problems considered the classical Fourier model which represents the thermal conduction heat transfer. On the other hand, the model is not precise for specific situations in real engineering phenomena since it produces a parabolic energy equation [24]. In 1948, Cattaneo [25] modified the Fourier model known as Maxwell-Cattaneo (MC) law by including a thermal relaxation time term to present the thermal inertia. Christov [26] re-examined the MC law by considering the time derivative of heat flux which resulting in a hyperbolic energy equation. Some studies were conducted to justify the uniqueness and stability of the solutions for the energy equation using the Cattaneo-Christov heat flux model (see [27][28][29]). In the past few years, many researchers have investigated the boundary layer problem of various fluids due to a stretching sheet using the non-Fourier model. Mustafa [30] used this model for Maxwell fluid and noticed that there was a relation between the thermal boundary layer thickness and the thermal relaxation time. Salahuddin et al. [31] emphasized the magnetohydrodynamics (MHD) flow of Williamson fluid with variable thickness. They discovered that the temperature profile decreased using the heat flux model compared to the classical Fourier model. Hayat et al. [32] considered stagnation flow of Jeffrey fluid towards a nonlinear stretching surface which resulted in the fluid temperature declining for higher thermal relaxation parameters. Malik et al. [33] analyzed the MHD flow of Casson fluid and also found that the temperature decreased using the Cattaneo Christov heat flux model. Very recently, a similar study conducted for nanofluids [24,34] and diverse non-Newtonian fluid models such as Oldroyd B-fluid [35], Eyring Powell [36,37], Carreau fluid [38], viscoelastic fluid [39,40] and second grade fluid [41,42] resulted in the non-Fourier model leading to the reduction of the temperature distribution. Mahmood et al. [43] designed the boundary layer flow of Casson nanofluid using Cattaneo-Christov heat flux model while Jamshed and Aziz [44] utilized the Casson hybrid TiO 2 -CuO/EG nanofluid on the stretched flow problem with entropy generation.
Double stratification of fluids appears due to the temperature and concentration variations, respectively. From the engineering perspective, many applications utilize the concept of double stratification, such as heat rejection from environment, thermal energy storage systems and solar energy because the better stratification corresponds to the higher energy performance. There are a few studies focused on mathematical analysis for the convective flow of nanofluids with the presence of stratification [45][46][47][48][49][50][51][52][53]. Ibrahim and Makinde [45] analyzed the natural convection and stable doubly stratified flow with the presence of nanoparticles. Mat Yasin et al. [46] dealt with the mixed convective flow in a thermally stratified porous medium saturated by a nanofluid using Tiwari and Das model. The impact of double stratification and mixed convective flow of the non-Newtonian fluids with the appearance of nanoparticles have been studied by Hussain et al. [47] and Abbasi et al. [48] using Maxwell and Jeffrey nanofluid models. Besthapu et al. [49] examined the influence of both thermal and solutal stratification on MHD nanofluid due to an exponentially stretching sheet. Daniel et al. [50,51] considered the effects of viscous dissipation and Joule heating on both steady and unsteady MHD mixed convection of doubly stratified nanofluid. Kandasamy et al. [52] investigated the effect of double stratification on MHD nanofluid over a porous vertical plate. The previously published results showed that the temperature and concentration distribution are decreased with the intensity of thermal and solutal stratification, respectively.
The preceding literatures give an inspiration for the authors to examine the behavior of MHD nanofluid flow with the double stratification phenomenon over a permeable shrinking/stretching sheet. The non-Fourier energy equation is utilized in the governing model whereas the water-based nanofluid is represented by the Buongiorno s model. The partial differential equations (PDEs) are simplified into a set of similarity differential equations by employing a suitable transformations and then numerically computed using the remarkable boundary value problem with fourth order accuracy (bvp4c) code in the MATLAB software. The results are graphically manifested for the nanofluid velocity, temperature, and concentration within the applicable range of the pertinent parameters. The study is also concerned with whether the non-Fourier energy model has relevant impact on the heat and mass transfer rates as compared to the classical (Fourier) energy model. Since the present work contemplates both shrinking flow and suction effect, the existence of non-unique solutions is expected . The implementation of stability analysis may certify the physical or real solution. To the best of the authors' knowledge, the results are new and have not been published before.

Problem Formulation
An incompressible and steady flow of a nanofluid towards a shrinking/stretching sheet with linear velocity U w (x) = ax is considered in the present study as illustrated in Figure 1. The wall temperature and concentration is represented by T w (x) = T 0 + Ax and C w (x) = C 0 + Ex, correspondingly where T w > T 0 and C w > C 0 . The ambient temperature and concentration are in the linear stratified form of T ∞ (x) = T 0 + Bx and C ∞ (x) = C 0 + Fx such that T 0 and C 0 are the beginning ambient temperature and concentration of the nanofluid. The following assumptions for the physical model are also examined in the present work: • The nanoparticles and the base fluid are assumed in a thermal equilibrium state.

•
The Buongiorno's model of nanofluid is used to incorporate the mixed effects of Brownian motion and thermophoresis.

•
The induced magnetic field is negligible as a result of the insignificant value of the magnetic Reynolds number

•
The Hall current effect is excluded with an assumption of no external electric field is applied to the physical model.

•
Thermal and solutal stratifications exist due to the appearance of thermal and solutal buoyancy forces.

Cattaneo Christov Heat Flux Model
The Cattaneo Christov heat flux model is in the vector form of: or can be written as given that q is the heat flux and For incompressible flow, Equation (3) will reduce to and eliminating q from Equations (1) and (4) where It is noticed that for λ 2 = 0, the Cattaneo Christov heat flux model is reduced to the classical Fourier energy equation.

Steady Flow Formulation
Under the boundary layer approximation, the governing model in the present study which consists of continuity, momentum, energy and concentration equations are given by (see Akbar et al. [24], Besthapu et al. [49] and Anwar et al. [54]): in conjunction with the initial and boundary conditions where (u, v) are the velocity components along the (x, y) directions, respectively, T is the nanofluid temperature, C is the nanoparticles volume fraction, ν is the kinematic viscosity, σ M is the electrical conductivity of the fluid, ρ f is the density of base fluid, ρ p is the density of nanoparticles, λ 2 is the thermal relaxation time, α = k ρC p is the thermal diffusivity of the fluid, τ 1 = (ρc) p (ρc) f is the ratio of heat capacity of the nanoparticles to the base fluid, D B is the Brownian diffusion coefficient, D T is the thermophoretic diffusion coefficient, V 0 = − √ aνS is the mass flux velocity, ε is the multiplier of shrinking/stretching parameter such that ε < 0 for shrinking sheet and ε > 0 for stretching sheet, accordingly whereas ε = 0 for static sheet.
The governing boundary layer problem in Equations (8)-(10) with the initial and boundary conditions (see Equations (11) and (12)) are transformed into a system of ordinary (similarity) differential equations by utilizing a set of transformations as listed below: where η is the similarity variable and ψ is the stream function such that which identically satisfies Equation (7). The resulting transformed ordinary differential equations (ODEs) align with the initial and boundary conditions are: 1 Pr where and prime implies the derivative with respect to the similarity variable η. Furthermore, f is the dimensionless stream function along the x−direction, θ is the dimensionless temperature, φ is the dimensionless nanoparticle volume fraction, M is the magnetic parameter, λ is the thermal buoyancy parameter, N is the solutal buoyancy parameter, Pr is the Prandtl number, Le is the Lewis number, Nb is the Brownian motion parameter, Nt is the thermophoresis parameter, δ 1 is the thermal stratification parameter, δ 2 is the solutal stratification parameter, γ T is the thermal relaxation parameter and S is the suction parameter, which are defined as It is worthwhile to mention here that λ > 0 and λ < 0 corresponds to the assisting and opposing flow, respectively. The system of ODEs in Equations (15)-(17) is reduced to the model by Akbar et al. [24] for the problem of unstratified nanofluid.
The physical quantities of interest in the present study are the dimensionless skin friction coefficient, local Nusselt number (heat transfer rate) and local Sherwood number (mass transfer rate) which are in the form of: respectively. Equation (21) is obtained by substituting Equations (13), (14) and (22) (surface heat and mass fluxes): into the dimensional form of local skin friction C f x , Nusselt number N u x and Sherwood number S h x :

Stability Analysis
A boundary layer problem may produce zero, unique, dual or multiple solutions depending on the imposition of appropriate physical parameters. Generally, for non-unique solutions, the first solution which satisfies the boundary condition i.e., (see Equation (19)) is denoted as the physical or real solution as compared to the other solutions. Hence, it is important to identify all the possible solutions that may arise in the boundary layer problem to avoid misinterpretation of the fluid flow characteristics. In certain cases, the second solution may exhibitthe same pattern of the real flow characteristics based on the velocity, temperature or concentration profiles. Therefore, it is necessary to validate the real solution through a proper analysis. The execution of the stability analysis is mathematically performed to verify the physical or real solution among all the solutions. There has been much current literature that discussed the importance, formulation and execution of the stability analysis (see Ismail et al. [55][56][57], Bakar et al. [58,59], Anuar et al. [60], Salleh et al. [61,62], Najib et al. [63], Jamaludin et al. [64] and Yahaya et al. [65]).
The solution is considered unstable or not real if there is an initial growth of disturbance in the solution. Theoretically, the disturbance may exponentially decay or grow with time and as a consequence, the stability formulation is initiated by considering an unsteady (time dependent) problem. Hence, for the present analysis, an unsteady form of Equations (8)-(10) need to be examined: New dimensionless time-based transformations with τ = at are adapted to the unsteady Equations (24)-(26): so that 1 Pr 1 PrLe where restrict to the conditions The representations in Equation (34) are used in the stability process where f 0 (η) = f (η), θ 0 (η) = θ(η) and φ 0 (η) = φ(η) are the similarity solutions of the steady problem (see Equations (15)-(19)), σ is an unidentified eigenvalue, F(η), G(η) and H(η) are small relative to f 0 (η), θ 0 (η) and The linearized eigenvalue problem is attained by substituting Equation (34) into Equations (28)- (33). 1 Pr 1 PrLe where align with the conditions The stability of the similarity (steady) solutions f 0 (η), θ 0 (η) and φ 0 (η) depends on the smallest eigenvalue, σ 1 by solving the governing linearized eigenvalue model in Equations (35)- (40). Relaxation of a boundary condition is required to possess a possible range of the computed eigenvalues. Hence, in the present work, the boundary condition F (η) → 0 as η → ∞ (see Equation (40)) is relaxed and replaced with the normalizing boundary condition F (0) = 1.

Results and Discussion
The present steady flow problem represented by the system of nonlinear ODEs in Equations (15)-(17) with the conditions (see Equation (18) and (19)) is carried out using the bvp4c (boundary value problem with fourth order accuracy) function in MATLAB software. To solve the suitable boundary value problems, 3-stage Lobatto IIIa (finite difference method) is used in the built-in-bvp4c function [66][67][68]. The main interest in the present study is to examine the effect of the pertinent parameters namely suction S, thermal relaxation γ T , thermal stratification δ 1 and solutal stratification δ 2 on the dimensionless velocity, temperature and concentration profiles. The authors are also concerned with the existence of dual solutions since suction and shrinking parameters are taken into account. In the present work, M = Le = 1, Nb = 0.3, Nt = 0.1, λ = −0.5, N = 0.5, δ 1 = 0.1, δ 2 = 0.2, ε = −1, S = 2.5 and γ T = 0.01 are considered as the main value for the numerical computations whereas Pr = 6.2 is selected to symbolize the water as the base fluid. It is worth to mention here that the results are new, hence, the velocity, temperature and concentration profiles with the inflation of control parameters may be different from the existing studies. However, a few of comparisons is made with the previously published results for limiting cases to validate the bvp4c method used in the current work as tabulated in Tables 1 and 2. The approximate relative error ε a , is also calculated and it shows that ε a between present and previous results are adequately small (0-0.3694%). It is clear from Table 2 that the thermal relaxation parameter γ T which resulting from the non-Fourier energy model can intensify the rate of heat transfer while reduce the rate of mass transfer.    Table 3 compares the numerical values of the heat and mass transfer rates for the problem of dual stratified nanofluid over a shrinking and stretching sheet. For the shrinking flow, −θ (0) approximately increases 5.83% and 11.59% when γ T = 0.005 and 0.01, respectively. In contrast, −φ (0) decreases 1.13% when γ T = 0.005 and 2.42% when γ T = 0.01. However, the percentage of increment for −θ (0) in the case of stretching flow is slightly higher as compared to the shrinking flow. In addition, the assisting flow case (λ = 0.3) and opposing flow case (λ = −0.5) in Tables 2 and 3, correspondingly have same pattern (increment/decrement) of the heat and mass transfer rates. Table 3. Numerical values of −θ (0) and −φ (0) for the case of dual stratified nanofluid over shrinking sheet (ε = −1) and stretching sheet (ε = 1) when M = Le = 1, λ = −0.5, N = 0.5, S = 2.3, Pr = 6.2, The existence of dual solutions for S > 2 in the present study is in accordance to the result by Miklavčič and Wang [2], Fang et al. [4], Jusoh et al. [5] and Soomro et al. [6]. It is apparent from these previous studies that the flow will produce two solutions by employing S > 2 for a suitable combination of the control parameters. It also can be seen from Figure 2 that the skin friction coefficient inclines as ε → −ε while declines when the stretching sheet (ε → +ε) is considered. The sheet is at static motion when ε = 0 for all values of S. Different characteristics of f (0) are noticed for the second solution which can justify the reliability of first solution as the physical solution. The utilization of suction for both shrinking and stretching flow also boosts the heat and mass transfer rates as illustrated in Figures 3 and 4. Suction may assist the movement of heated fluid particles towards the wall and subsequently, increases both heat and mass transfer rates. The nanofluid velocity profile, as depicted in Figure 5, elevates while the temperature and concentration (see Figures 6 and 7) deteriorate with an upsurge of the suction parameter. This is due to the suction which can diminish the momentum boundary layer thickness and therefore enhance the flow near to the sheet [5]. Besides, suction also allows the heated fluid movement towards the wall surface and develops both thermal and concentration boundary layer thickenesses.      Figures 12-14. Surprisingly, the combined effects of suction S, thermal buoyancy λ and shrinking (ε = −1) parameters gave the reverse effect of the magnetic parameter M on all the profiles. Theoretically, a magnetic field may induce a drag reduction known as Lorentz force which gives resistance to the fluid flow and hence decelerates the velocity as reported by Rashidi et al. [69]. However, the shrinking flow with a high value of suction parameter (S > 2) may contribute to the increment of velocity profile as plotted in Figure 12 and depreciate both thermal and concentration boundary layer thickness (see Figures 13 and 14).      The present boundary layer problem produced two solutions within a certain range of the pertinent parameters. Hence, it is necessary to prove the real or physical solution by a valid method rather than conclude based on the velocity, temperature or concentration profiles. The stability analysis is successfully conducted by solving Equations (35)-(40) using the bvp4c algorithm in MATLAB software. Figure 15 illustrated the smallest eigenvalue σ 1 of both solutions towards ε where σ 1 act as a determinant of the stability solutions. Positive σ 1 implies that the flow is stable whereas negative σ 1 indicates an initial growth of disturbances which resulting that the flow is unstable [55][56][57][58][59][60][61][62][63][64][65]. It is validated from Figure 15 that the first and second solutions have positive and negative σ 1 , respectively which indicates that the first solution is the real solution. Furthermore, the smallest eigenvalue σ 1 for both solutions is σ 1 → 0 as ε → ε c which proves the formulation of the present stability analysis.

Conclusions
The present work scrutinizes the boundary layer flow coupled with heat and mass transfer of a doubly stratified nanofluid using a non-Fourier energy model, or known as the Cattaneo-Christov heat flux model. The governing model in PDEs are reduced to a system of similarity differential equations and then numerically solved using bvp4c code which is programmed in the MATLAB software. The conclusions are as follows: • Thermal relaxation parameter γ T which is resulting from the application of the non-Fourier energy model, gives a different estimation of the heat and mass transfer rates as compared to the classical energy equation. The rate of heat transfer increases whereas the mass transfer rate decreases with the improvement of the thermal relaxation parameter for both stretching and shrinking cases.

•
The skin friction coefficient inflates for the shrinking case (ε < 0), deflates for the stretching case (ε > 0) and has no change for the static plate (ε = 0) when the suction parameter S is enhanced. Higher values of the suction parameter i.e., S > 2 also contribute to the appearance of the dual solutions in the present work.

•
The stability analysis is conducted and it is mathematically proved that the first solution is the physical/real solution.

•
The combined effects of buoyancy parameter (λ and N) with higher values of the suction parameter S in the shrinked flow gave the reverse effect of the Lorentz force to the velocity, temperature and concentration profiles.

•
Both temperature and concentration slightly inclines with the imposition of the thermal δ 1 and solutal δ 2 stratification parameters, respectively.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: