Three-Dimensional Magnetohydrodynamic Mixed Convection Flow of Nanofluids over a Nonlinearly Permeable Stretching / Shrinking Sheet with Velocity and Thermal Slip

Anuar Jamaludin 1, Roslinda Nazar 2,* and Ioan Pop 3 1 Department of Mathematics, Universiti Pertahanan Nasional Malaysia, Kuala Lumpur 57000, Malaysia; mohdanuar@upnm.edu.my 2 School of Mathematical Sciences, Faculty of Science and Technology, Universiti Kebangsaan Malaysia, UKM Bangi, Selangor 43600, Malaysia 3 Department of Mathematics, Babeş-Bolyai University, R-400084 Cluj-Napoca, Romania; popm.ioan@yahoo.co.uk * Correspondence: rmn@ukm.edu.my


Introduction
At present, mixed convection or combined forced and free (natural) convection has emerged as an essential aspect in studies pertaining to heat transfer that revolves around both nature and engineering applications, for example, electronic cooling, nuclear reactors technology, and nanotechnology, to name a few.Besides, since the past decade, most studies concerning mixed convection flow analysis have emphasized the existence of dual (upper and lower branch) solutions to oppose and assist flows in a certain range of the buoyancy (or mixed convection) parameter.For instance, Ramachandran et al. [1] looked into the laminar mixed convection in two-dimensional stagnation flow adjacent to vertical heated surfaces with specified wall temperature and specified wall heat flux cases.They discovered several dual solutions for some ranges of buoyancy parameter in the opposing flow region.Meanwhile, Devi et al. [2] extended the work of Ramachandran et al. [1] for unsteady flow.Dual solutions were found to exist for a certain range of the buoyancy parameter in the opposing flow region.In contrast to the studies conducted by Ramachandran et al. [1] and Devi et al. [2], Rhida and Curie [3] revealed that dual solutions exist in both opposing and assisting flow regions.Other than that, a broader perspective was investigated by Merkin [4] who discovered that the upper branch of solutions was stable, but unstable for the lower branch.On top of that, he also looked into the issue of mixed convection flow in a saturated porous medium.Later, Merkin and Mahmmod [5] analysed two-dimensional mixed convection flow with prescribed wall heat flux, hence describing a comprehensive report of how both lower and upper branch solutions behave over all possible ranges of the buoyancy parameter deduced.In addition, Ishak et al. [6] extended the work of Merkin and Mahmmod [5] by determining the effects of suction and injection.They also examined the existence of dual solutions in both opposing and assisting flows for two cases: permeable and impermeable.Moreover, several other published papers have reported the existence of dual solutions in mixed convection flow (see [7][8][9][10] for more details).
Furthermore, studies concerning momentum and heat transfer characteristics of the boundary layer flow subject to magnetic field have gained great interest due to their benefits in controlling the cooling rate.In fact, this cooling rate has a significant function in enhancing the desired attributes of the final product to be applied in industrial fields, such as polymer technology.In view of this, numerous researchers like Ishak et al. [11], Makinde and Chinyoka [12], as well as Gireesha et al. [13], have examined the magnetohydrodynamic (MHD) effect on boundary layer flow and heat transfer issues in both Newtonian and non-Newtonian fluids with various aspects.Nevertheless, these studies dismissed the notion of mixed convection.In fact, researches about MHD mixed convection flow have been carried out by many authors.For example, Ishak et al. [14] extended the work of Ramachandran et al. [1] by investigating the steady MHD flow, where they discovered that magnetic parameter increased the range of solutions, which in turn, decelerated the separation between the boundary layers.Other than that, Ali et al. [15] studied the effects of MHD on mixed convection stagnation point flow over a vertical surface.As a result, they discovered that the velocity, temperature, and induced magnetic field profiles were affected by the buoyancy and magnetic parameters for both opposing and assisting flows.Besides, studies on unsteady MHD mixed convection flow in the presence of suction/injection and radiation, which is applied to a vertical stretching/shrinking surface, had been studied by Sandeep et al. [16].Meanwhile, in an incompressible Casson fluid, Sharada and Shankar [17] looked into the steady three-dimensional mixed convection flow past an exponential stretching sheet, together with the effects of magnetic field and heat generation.Moreover, recently, Sivasankaran et al. [18] included the combined effects of MHD, chemical reactions, radiation, and velocity slip on mixed convection flow near a stagnation point over a vertical plate in a porous medium, along with convective boundary condition.
Additionally, the field of nanofluids has emerged as a significant topic in the heat transfer research at present due to its various applications, such as in domestic refrigerator, cancer therapy, solar water heating, cooling of electronics, coolants for nuclear reactors, vehicle engine cooling, transformer cooling and the list is endless (see Wong and Leon [19] and Saidur et al. [20]).Furthermore, nanofluid is a dilute suspension of ultra-fine nano-sized solid particles (nanoparticles with size smaller than ~100 nm, made from various materials, for instance, Cu, Ag, Al 2 O 3 , CuO, and TiO 2 ) in the original base fluids, for instance, water, oil, glycerol, and ethylene glycol.Nonetheless, these base fluids have low thermal conductivity.Therefore, an innovative approach has been developed to enhance the thermal conductivity of these base fluids.Moreover, several experimental studies have indicated increment in thermal conductivity of these base fluids when the nanoparticles are suspended in the base fluids (Eastman et al. [21], Liu et al. [22] and Hwang et al. [23]).In fact, the thermal conductivity of nanofluid is substantially higher than that of the base fluid, which in turn, boosts the performance of heat transfer.Besides, the two common models that have always been applied by researchers in studying boundary layer flow and heat transfer of nanofluid are the Tiwari-Das model [24] and the Buongiorno model [25].While the Buongiorno model highlights the thermophoresis and Brownian motion effects, the Tiwari-Das model focuses on the volumetric fraction of nanoparticles.
Meanwhile, Yasin et al. [26] studied the steady mixed convection flow towards a vertical permeable surface in a porous medium saturated by nanofluids, whereas Pal et al. [27] looked into the effects of internal heat generation/absorption, suction/injection, and thermal radiation on steady mixed convection stagnation point flow of nanofluids.Other than that, Pal and Mandal [28] determined the influences of thermal radiation and viscous dissipation on a steady mixed convection flow of nanofluids near the stagnation point region past a stretching/shrinking sheet that is embedded in the porous medium in the presence of non-uniform heat source/sink.Furthermore, the study of unsteady mixed convection flow of nanofluids was also carried out by Mahdy [29], Abdullah et al. [30] and Oyelakin et al. [31].However, flow of MHD was disregarded in their studies.Nevertheless, the study of mixed convection flow of nanofluids under the influence of magnetic field was indeed investigated by Yazdi et al. [32], Haroun et al. [33], Mustafa et al. [34], as well as Mahanthesh et al. [35], to name a few.
On the other hand, increasing interest to study the effect of slip conditions on both flow and heat transfer of nanofluids with various physical conditions has been noted [36][37][38][39].Recently, Mahanthesh et al. [35] examined the combined effects of MHD, thermal radiation, velocity slip, and thermal slip on three-dimensional mixed convection flow of nanofluids over a nonlinear stretching sheet.Motivated by the above mentioned studies, the objective of this present study is to extend the work of Mahanthesh et al. [35] on the steady three-dimensional MHD mixed convection flow over a nonlinear stretching/shrinking sheet with suction and slip conditions using two different types of water-based nanofluids, namely Cu-water and Ag-water nanofluids, through the use of mathematical nanofluid model, which has been proposed by Tiwari and Das [24].Nonetheless, the effect of thermal radiation has been excluded from this study.Furthermore, the case of shrinking sheet in the presence of suction with stability analysis has been determined for dual (upper and lower branch) solutions which have not been considered by Mahanthesh et al. [35].Therefore, it is believed that the results retrieved from this study are new and differ from those reported by Mahanthesh et al. [35].Moreover, it must be noted here that the stability analysis of the dual solutions had been performed to prove the stability of upper branch solutions, as well as the instability of the lower branch solutions.In addition, the behaviour of the reduced skin friction coefficients and the reduced local Nusselt number on selected governing parameters, namely mixed convection, velocity slip, temperature slip, and magnetic parameters, appear to be the main focus of the present study.These numerical results are presented graphically to display the effects of these parameters on the reduced skin friction coefficients and the reduced local Nusselt number.

Problem Formulation
Steady, three-dimensional, mixed convection, MHD laminar boundary layer flow and heat transfer with electrically conducting fluid over a nonlinear, vertical, permeable, and stretching/shrinking sheet has been considered for investigation.Besides, this study has employed water-based nanofluid with varied types of metallic solid nanoparticles, namely copper (Cu) and silver (Ag).The thermophysical properties of these base fluids and nanoparticles are presented in Table 1 [40][41][42].Moreover, the sheet is extended nonlinearly in two lateral xand ydirections, as illustrated in Figure 1.It further portrays that coordinate z-is measured in the direction of normal to the sheet, while coordinates x and y are in the plane of the sheet with x-axis measured in vertical direction.Besides, it has been assumed that the flow subjected to the variable of transverse magnetic field B, is applied parallel to z-axis.The induced magnetic field is also assumed to be small compared to the applied magnetic field, hence it is neglected.Other than that, the vertical sheet is assumed to move continuously in both x-and y-directions with velocities, w u and w v , respectively, while the mass flux velocity is , w w the temperature of the sheet is w T and the constant temperature of the ambient fluid is T  , which are determined later.Under these conditions, the governing boundary layer equations for both flow and heat transfer are given in the following: along with the following boundary conditions: where u, v and w are the velocity components along the x-, y-and z-axes, respectively, c is the stretching/shrinking parameter with 0 c  for a shrinking sheet and 0 c  for a stretching sheet, g refers to the acceleration due to gravity, B is the constant magnetic field,  It further portrays that coordinate zis measured in the direction of normal to the sheet, while coordinates x and y are in the plane of the sheet with x-axis measured in vertical direction.Besides, it has been assumed that the flow subjected to the variable of transverse magnetic field B, is applied parallel to z-axis.The induced magnetic field is also assumed to be small compared to the applied magnetic field, hence it is neglected.Other than that, the vertical sheet is assumed to move continuously in both xand y-directions with velocities, u w and v w , respectively, while the mass flux velocity is w w , the temperature of the sheet is T w and the constant temperature of the ambient fluid is T ∞ , which are determined later.Under these conditions, the governing boundary layer equations for both flow and heat transfer are given in the following: along with the following boundary conditions: where u, v and w are the velocity components along the x-, yand zaxes, respectively, c is the stretching/shrinking parameter with c < 0 for a shrinking sheet and c > 0 for a stretching sheet, g refers to the acceleration due to gravity, B is the constant magnetic field, ν n f denotes the kinematic viscosity of the nanofluid, (ρβ) n f denotes the thermal expansion coefficient of the nanofluid, ρ n f refers to the density of the nanofluid, σ n f is the electrical conductivity of the nanofluid and α n f is the thermal diffusivity of the nanofluid.Furthermore, the effective properties of the nanofluid are described in the following (see [40,41,[43][44][45]): where µ n f is the dynamic viscosity of the nanofluid, as described by the Brinkman model, µ f denotes the dynamic viscosity of the base fluid, φ denotes the solid volume fraction of nanoparticles, ρ f and ρ s are the density of the base fluid and the solid nanoparticle, respectively, ρC p n f is the heat capacitance of the nanofluid, k n f is the thermal conductivity of the nanofluid, as approximated by Maxwell-Garnett model, in which the subscript 'f ' refers to the base fluid while 's' refers to the solid nanoparticle.Furthermore, N 1 and K 1 are the velocity slip factor and the thermal slip factor, respectively (see Ibrahim and Shankar [37]), and it has been assumed that u w , v w , w w , T w , B, N 1 , and K 1 take the following forms: Here, s refers to the suction/injection parameter with s < 0 for injection and s > 0 for suction, B 0 , N 0 and K 0 are the constants, T 0 is the constant characteristic temperature with T 0 < 0 that implies cooled surface (opposing flow) while T 0 > 0 corresponds to heated surface (assisting flow), N and K are the dimensionless slip parameters, and lastly, a and n are the positive constants.
Moreover, the governing Equations ( 1)-( 4) along with boundary conditions (5), can be converted into ordinary differential equations by the dimensionless functions u, v, w and θ, which are related to the similarity variable η, as follows: Note that f and h are the dimensionless quantities, f and h are the dimensionless velocity, θ is the dimensionless temperature and prime denotes differentiation with respect to η. Substituting similarity transformation (8) into governing Equations ( 1)-( 4), the following equations are obtained: Appl.Sci.2018, 8, 1128 with the following boundary conditions: In fact, throughout this paper, the governing parameters are as follows: Pr refers to the Prandtl number, λ denotes the mixed convection parameter with λ < 0 corresponding to opposing flow and λ > 0 corresponding to assisting flow, M denotes the magnetic parameter, N is the velocity slip parameter and K is the temperature slip parameter, which are defined as: Here, the local Grashof number, Gr x and the local Reynolds number, Re x along x-direction, given by the following: The important physical quantities in this study are the skin friction coefficients, C f x and C f y and the local Nusselt number, Nu x which are defined as follows: where τ zx and τ zy are shear stresses at the walls denoted as zx and zy, respectively, while q w is the heat flux at surface z = 0.The shear stresses are defined as: Using ( 7) in ( 15) and ( 16), the following are obtained: where Re x = u w (x + y)/ν f is the local Reynolds number along x-direction and Re y = v w (x + y)/ν f is the local Reynolds number along y-direction.

Stability Analysis
The present numerical results depict the existence of dual solutions in a certain range of λ, N, K, and M. Therefore, the stability of dual solutions has to be determined.In fact, prior studies have established that the stability of dual solutions can be determined by performing stability analysis, as suggested by Merkin [4].Thus, for the purpose of conducting the stability analysis, an unsteady-state problem is considered.Equation (1) holds, while Equations ( 2)-( 4) are substituted as follows: where t denotes the time.Additionally, the new dimensionless functions of u, v, w, and θ, are initiated in relation to the similarity variable, η and the second similarity variable, τ for the unsteady-state problem ( 18)-( 20) as follows: Thus, the unsteady-state problem ( 18)-( 20) can be written as: subject to the boundary conditions: Furthermore, in order to test the stability of solutions f (η) = f 0 (η), h(η) = h 0 (η), and θ(η) = θ 0 (η) in satisfying the steady-state problem ( 9)-( 12), the following is written: where γ is an unknown eigenvalue parameter.Besides, note that F(η, τ), H(η, τ), and G(η, τ) are smaller relative to f 0 (η), h 0 (η), and θ 0 (η).Furthermore, by substituting (26) into Equations ( 22)- (24), along with the boundary conditions (25), the following equations are obtained: Appl.Sci.2018, 8, 1128 1 Pr subject to the boundary conditions: Furthermore, as depicted by Weidman et al. [46], in order to identify the initial growth or the decay of solution ( 26), τ = 0, is applied and hence, the following are obtained: F = F 0 (η), H = H 0 (η), and G = G 0 (η).Furthermore, in order to test the numerical procedure, the following linear eigenvalue problem has to be solved: 1 1 Pr (33) subject to the boundary conditions: The solutions f 0 (η), h 0 (η), and θ 0 (η) are determined from the steady-state problem ( 9)-( 12) before plugging back into Equations ( 31)- (33), in which the linear eigenvalue problem (31)-( 34) could be solved.Moreover, Harris et al. [47] suggested relaxing an appropriate boundary condition on F 0 (η), H 0 (η), or G 0 (η) to better find the range of γ.Thus, the condition F 0 (η) → 0 as η → ∞ is relaxed and for a fixed value of γ, the linear eigenvalue problem (31)- (34) are solved, along with the new boundary condition; F 0 (0) = 1.The obtained solutions offers an infinite set of possible eigenvalues γ 1 < γ 2 < γ 3 < . .., where γ 1 is the minimum eigenvalue.Furthermore, if γ 1 > 0, then there is an initial decay of disturbances and the flow becomes stable.Other than that, if γ 1 > 0, then there is an initial growth of disturbances and the flow becomes unstable.

Results and Discussion
The steady-state problem ( 9)-( 12) is solved in a numerical manner by using the bvp4c programme from MATLAB, as described by Shampine et al. [48].The bvp4c further requires an initial guess, which should reveal the attributes of the dual solutions and satisfy boundary conditions (12).Besides, the relative error tolerance on residuals has been set to 10 −5 .A suitable finite value of η → ∞ had been opted, namely η = η ∞ = 10 for the upper and lower branch solutions.The numerical results are retrieved for varied values that involved dimensionless parameters, namely mixed convection parameter λ, velocity slip parameter N, temperature slip parameter K, and magnetic parameter M, in the presence of nanoparticles.On top of that, nanoparticles Cu and Ag have been considered with water as base fluids.In fact, in order to cover the limitation of experimental data, the choice of values is as determined by prior researchers.The value of Prandtl number is considered as Pr = 6.2 (water-based nanofluid), while the values of n, φ, c and s are n = 3, φ = 0.2, c = −1 and s = 1, respectively.These values are fixed throughout this study (except for comparison purpose).Besides, the values of dimensionless parameters λ, N, K, and M are given in the respective figures and tables.
Next, in order to validate the present numerical method employed, the present results of this study are compared with findings obtained from prior works of Mahanthesh et al. [35] and Jain and Choudhary [49].Table 2 presents the comparative studies of the present results for reduced skin friction coefficient f (0) for several values of M and n, when c = 1 and λ = s = φ = N = K = 0.Moreover, as noted in Table 2, it is interesting to observe that the present results are in very good agreement with the above mentioned works.Hence, one can claim that the present results are indeed correct and accurate.
Nonetheless, interestingly, the present numerical results do indicate the existence of dual (upper and lower branch) solutions.As discussed earlier, the stability analysis is performed in this study in order to test the stability of the dual solutions.Therefore, the linear eigenvalue problem ( 31)-( 34) is solved numerically using the MATLAB programme bvp4c.In addition, the minimum eigenvalues γ 1 , for several values of λ, N, K, and M, are given in Table 3.The results in Table 3 clearly show that γ 1 > 0 is for upper branch solution, while γ 1 < 0 reflects the lower branch solution.These results are consistent with those of other studies (see [4,8,9,39,47]) and suggest that the upper branch solution is indeed stable, but unstable for the lower branch solution.Note: M is the magnetic parameter, n is the positive constant, N is the velocity slip parameter, K is the temperature slip parameter, λ is the mixed convection parameter.
Additionally, Figures 2-19 display the variations of the reduced skin friction coefficients f (0), h (0) and the reduced local Nusselt number −θ (0) with λ for varied values of parameters N, K, and M for both Cu-water and Ag-water nanofluid.The solid lines are used throughout this paper in order to refer to the upper branch solution, whereas the dotted lines represent the lower branch solution.One can observe that the dual (upper and lower branch) solutions do exist for Equations ( 9)-( 11), subjected to boundary conditions (12) within the range of λc < λ, where λc is the critical value of parameter λ, while no solution exists for λ < λc.Besides, interestingly, the unique (upper branch) solution is attained when λ > 0 (assisting flow), while the dual solutions are found only for λ < 0 (opposing flow).Furthermore, for λ < λc, both full Navier-Stokes and energy equations have to be solved.In addition, it is clear that all |λc| values for Ag-water nanofluid are higher than those for Cu-water nanofluid.Hence, one can assume that the presence of Ag nanoparticles in the water-based nanofluid is more capable of delaying the boundary layer separation, in comparison to that for Cu nanoparticles.
Other than that, Figures 2-7 illustrate the variations of f (0), h (0), and −θ (0) against λ for varied values of velocity slip parameter N for both Cu-water and Ag-water nanofluids.As a result, it was observed that the value of |λc| increased as N increased.This indicates that the impact of N escalates the range of existence of the solutions for the boundary value problem ( 9)-( 12).On top of that, this hints that the increment of velocity slip could decelerate the separation of boundary layer in the opposing flow region.Meanwhile, for a fixed value of N, the upper branch solution decreases with the decrease of λ until it reaches a critical value at λc, for both opposing and assisting flows.Besides, the corresponding lower branch solution decreases with the increase of λ for an opposing flow.Nonetheless, the values of −θ (0) (see Figures 6 and 7) for the lower branch solution decreased rapidly with the increase of λ.Moreover, the values of −θ (0) for the upper branch solution are extremely large, when compared to the corresponding lower branch solution in the neighborhood of λ = 0. Furthermore, Figures 2-5 demonstrate that at any λ station within the range of λc < λ, the values of f (0) and h (0) for the dual solutions decrease as N increases.In contrast, Figures 6 and 7 display that at any λ station in the range of λc < λ < 0, the values of −θ (0) for the dual solutions escalate with the increase in N.          increases, hinting that the magnetic effect delays the boundary layer separation for the opposing flow region, as well as increment within the range of existence for solutions of Equations ( 9)-( 11    increases, hinting that the magnetic effect delays the boundary layer separation for the opposing flow region, as well as increment within the range of existence for solutions of Equations ( 9)-( 11   portray an increment in the upper branch solution, as well as a decrease in the lower branch solution.Meanwhile, in Figures 8-13, the variations of f (0), h (0), and −θ (0) are plotted with λ for several values of temperature slip parameter K for both nanoparticles in water-based nanofluids.As noted from the figures, the critical value |λc| escalates as the parameter K is increased.Thus, it is possible that the increment of temperature slip can delay the separation of the boundary layers in the opposing flow region.When seen together, these results indicate that the range of existence of the solutions for the boundary value problem ( 9)-( 12) increases with increment in K.In addition, Figures 8-11 point out that the values of f (0) and h (0) for the upper branch solution increase with increment in K within the domain of λc < λ < 0. On the other hand, the opposite trend is observed when λ > 0. Additionally, Figures 12 and 13 portray that at any λ within the range of λc < λ < 0, the values of −θ (0) for both solution branches increase as K increases.

 
 profiles against  , for both Cu-water and Ag-water nanofluids.In addition, it is vivid from these figures that the aforementioned profiles do satisfy the far field boundary conditions (12) asymptotically.Hence, it validates and supports the numerical results retrieved for Equations ( 9)- (11), which are subject to boundary conditions (12), and admits the existence of dual nature of the solutions illustrated in Figures 2-19.Furthermore, the upper branch solution for these profiles portray a thinner boundary layer, in comparison to that of the lower branch solution.Additionally,

 
 profiles against  , for both Cu-water and Ag-water nanofluids.In addition, it is vivid from these figures that the aforementioned profiles do satisfy the far field boundary conditions (12) asymptotically.Hence, it validates and supports the numerical results retrieved for Equations ( 9)- (11), which are subject to boundary conditions (12), and admits the existence of dual nature of the solutions illustrated in Figures 2-19.Furthermore, the upper branch solution for these profiles portray a thinner boundary layer, in comparison to that of the lower branch solution.Additionally,   confirms the notion that nanoparticle with lower thermal conductivity, i.e., Ag, displays better improvement for heat transfer, when compared to Cu.Furthermore, Figures 14-19 explain the impact of magnetic parameter M on f (0), h (0) and −θ (0).In addition, these figures indicate that the strength of |λc| escalates as parameter M increases, hinting that the magnetic effect delays the boundary layer separation for the opposing flow region, as well as increment within the range of existence for solutions of Equations ( 9)- (11) with boundary conditions (12).Simultaneously, for a fixed value of M, the upper branch solution is found to decrease with the decrease of λ within the domain of λ < λc.The corresponding lower branch solution decreases as λ is increased when λc < λ < 0. Nonetheless, Figures 18 and 19 reveal that the values of −θ (0) for the lower branch solution display rapid decrease with the increase of λ.Moreover, the values of −θ (0) for the lower branch solution are large, in comparison to the corresponding upper branch solution in a region close to λ = 0. Besides, Figures 14,15,18 and 19 illustrate that at any λ station within the range of λc < λ, for increment of M, the values of f (0) (within the domain of λc < λ) and −θ (0) (within the domain of λc < λ < 0) portray an increment in the upper branch solution, as well as a decrease in the lower branch solution.Furthermore, Figures 16 and 17 point out that the values of h (0) for the upper branch solution display a similar trend with that of f (0).Nevertheless, the trend for the lower branch solution differs, when compared to the trend observed for the upper branch solution.Finally, Figures 20-22 offer sample of dimensionless velocity f (η), h (η), and temperature θ(η) profiles against η, for Cu-water and Ag-water nanofluids.In addition, it is vivid from these figures that the aforementioned profiles do satisfy the far field boundary conditions (12) asymptotically.Hence, it validates and supports the numerical results retrieved for Equations ( 9)- (11), which are subject to boundary conditions (12), and admits the existence of dual nature of the solutions illustrated in Figures 2-19.Furthermore, the upper branch solution for these profiles portray a thinner boundary layer, in comparison to that of the lower branch solution.Additionally, Figures 20 and 21 present for a fixed values of n, φ, c, s, λ, N, K, and M, the upper branch solution for f (η) and h (η) are higher for Ag-water nanofluid.Other than that, Figure 22 shows that the upper branch solution of θ(η) for Cu-water and Ag-water nanofluids are almost identical.Nonetheless, the temperature for Ag-water nanofluid is found to be lower.Therefore, this study confirms the notion that nanoparticle with lower thermal conductivity, i.e., Ag, displays better improvement for heat transfer, when compared to Cu.

Conclusions
The problem related to steady three-dimensional MHD mixed convection flow over a permeable vertical stretching/shrinking sheet with slip conditions in a nanofluid has been studied numerically in order to demonstrate the combined effects of velocity slip, temperature slip, and magnetic field on Cu-water and Ag-water nanofluids.Besides, with the aid of similarity transformations, the governing equations have been reduced to ordinary differential equations.Furthermore, the resulting equations subjected to the associated boundary conditions have been solved numerically for both the assisting and opposing flow cases, especially for the representative values of the selected governing parameters.In addition, a detailed discussion pertaining to the aforementioned effects of


A stability analysis was performed to confirm that the upper branch solution is indeed stable, whereas unstable for the lower branch solution.

Conclusions
The problem related to steady three-dimensional MHD mixed convection flow over a permeable vertical stretching/shrinking sheet with slip conditions in a nanofluid has been studied numerically in order to demonstrate the combined effects of velocity slip, temperature slip, and magnetic field on Cu-water and Ag-water nanofluids.Besides, with the aid of similarity transformations, the governing equations have been reduced to ordinary differential equations.Furthermore, the resulting equations subjected to the associated boundary conditions have been solved numerically for both the assisting and opposing flow cases, especially for the representative values of the selected governing parameters.In addition, a detailed discussion pertaining to the aforementioned effects of f (0), h (0), and −θ (0) is given.Other than that, several significant observations attained from this study are given in the following:

•
Dual (upper and lower branch) solutions are found for each value of λ in the opposing flow region with the solution curves bifurcating at the critical values of λc, while the lower branch solution could not be continued further to the point where λ reflected zero.

•
A stability analysis was performed to confirm that the upper branch solution is indeed stable, whereas unstable for the lower branch solution.

•
The value of |λc| increases with the increase of N, K, and M, thus, the boundary layer separation can be delayed by increasing the values of the aforementioned parameters.

Figure 1 .
Figure 1.Physical model and coordinate system.
nf  denotes the kinematic viscosity of the nanofluid,   nf  denotes the thermal expansion coefficient of the

Figure 1 .
Figure 1.Physical model and coordinate system.
Appl.Sci.2018, 8, x FOR PEER REVIEW 12 of 26 decreased rapidly with the increase of  .Moreover, the values of large, when compared to the corresponding lower branch solution in the neighborhood of 0.  Furthermore, Figures 2-5 demonstrate that at any  station within the range of , solutions decrease as N increases.In contrast, Figures6 and 7display that at any  station in the range of solutions escalate with the increase in N.

Figure 2 .
Figure 2. Variation of f (0) with λ for several values of N for Cu-water nanofluid.

Figure 3 .
Figure 3. Variation of f (0) with λ for several values of N for Ag-water nanofluid.

Figures 14 -.
Figures 14-19 explain the impact of magnetic parameter M on

Figure 4 .
Figure 4. Variation of h (0) with λ for several values of N for Cu-water nanofluid.

Figure 5 .
Figure 5. Variation of h (0) with λ for several values of N for Ag-water nanofluid.

Furthermore, Figures 16
and 17 point out that the values of   0 h for the upper branch solution display a similar trend with that of   0 f  .Nevertheless, the trend for the lower branch solution differs, when compared to the trend observed for the upper branch solution.

Figures 20 and 21
Figures 20 and 21 present for a fixed values of n,  , c, s,  , N, K, and M, the upper branch solution for

Figure 8 .
Figure 8. Variation of f (0) with λ for several values of K for Cu-water nanofluid.

Figures 20 and 21
Figures 20 and 21 present for a fixed values of n,  , c, s,  , N, K, and M, the upper branch solution for

Figure 9 .
Figure 9. Variation of f (0) with λ for several values of K for Ag-water nanofluid.

Figure 15 .
Figure 15.Variation of f (0) with λ for several values of M for Ag-water nanofluid.

Figure 18
Figure 18.Variation of
given.Other than that, several significant observations attained from this study are given in the following:  Dual (upper and lower branch) solutions are found for each value of  in the opposing flow region with the solution curves bifurcating at the critical values of c  , while the lower branch solution could not be continued further to the point where  reflected zero.

Table 1 .
Thermophysical properties of base fluid and nanoparticles.

Table 2 .
Numerical values of f (0) for several values of M and n.
• Ag-water nanofluid is more capable of delaying the boundary layer separation, in comparison to Cu-water nanofluid.•Ag-waternanofluid displayed better enhancement for heat transfer, when compared to that for Cu-water nanofluid.