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

Abstract

In this paper, the steady three-dimensional magnetohydrodynamic (MHD) mixed convection flow of nanofluids over a permeable vertical stretching/shrinking sheet with slip conditions is investigated in a numerical manner. As such, two types of nanofluids, Cu-water and Ag-water, had been considered. A similarity transformation was employed to reduce the governing equations to ordinary differential equations, which were then solved numerically using the MATLAB (Matlab R2015a, MathWorks, Natick, MA, USA, 1984) programme bvp4c. The numerical solutions derived from the ordinary differential equations subjected to the associated boundary conditions, were obtained to represent the values of the mixed convection parameter. Dual (upper and lower branch) solutions were discovered in the opposing flow region of the mixed convection parameter. A stability analysis was carried out to prove that the upper branch solution was indeed stable, while the lower branch solution was unstable. The significant effects of the governing parameters on the reduced skin friction coefficients, the reduced local Nusselt number, as well as the velocity and temperature profiles, were presented graphically and discussed in detail.

1. 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₂O₃, CuO, and TiO₂) 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.

2. 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. Thermophysical properties of base fluid and nanoparticles.

Physical Properties	Water	Cu	Ag
$C_p$(J kg−1 K−1)	4179	385	235
$\rho$ (k gm−1)	997.1	8933	10,500
$k$ (Wm−1 K−1)	0.613	400	429
$\beta \times {10}^{-5}$ (K−1)	21	1.67	1.89
$\sigma$ (Sm−1)	5.5 $\times$ 10−6	59.6 $\times$ 106	6.3 $\times$ 107

[40,41,42]. Moreover, the sheet is extended nonlinearly in two lateral x- and y- directions, 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, $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_{\infty }$, which are determined later. Under these conditions, the governing boundary layer equations for both flow and heat transfer are given in the following:

$$\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}+\frac{\partial w}{\partial z}=0,$$

(1)

$$u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}+w\frac{\partial u}{\partial z}={\nu }_{nf}\frac{{\partial }^{2}u}{\partial z^2}-\frac{{\sigma }_{nf}{B}^{2}u}{{\rho }_{nf}}+\frac{g{\left(\rho \beta \right)}_{nf}}{{\rho }_{nf}}\left(T-T_{\infty }\right),$$

(2)

$$u\frac{\partial v}{\partial x}+v\frac{\partial v}{\partial y}+w\frac{\partial v}{\partial z}={\nu }_{nf}\frac{{\partial }^{2}v}{\partial z^2}-\frac{{\sigma }_{nf}{B}^{2}v}{{\rho }_{nf}},$$

(3)

$$u\frac{\partial T}{\partial x}+v\frac{\partial T}{\partial y}+w\frac{\partial T}{\partial z}={\alpha }_{nf}\frac{{\partial }^{2}T}{\partial z^2},$$

(4)

along with the following boundary conditions:

$$\begin{array}{l}u=c{u}_w+N_1\frac{\partial u}{\partial z},v=c{v}_w+N_1\frac{\partial v}{\partial z},w=w_w,T=T_w+K_1\frac{\partial T}{\partial z}\mathrm{at}z=0,\\ u\to 0,v\to 0,T\to T_{\infty }\mathrm{as}z\to \infty ,\end{array}$$

(5)

where u, v and w are the velocity components along the x-, y- and z- axes, 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, ${\nu }_{nf}$ denotes the kinematic viscosity of the nanofluid, ${\left(\rho \beta \right)}_{nf}$ denotes the thermal expansion coefficient of the nanofluid, ${\rho }_{nf}$ refers to the density of the nanofluid, ${\sigma }_{nf}$ is the electrical conductivity of the nanofluid and ${\alpha }_{nf}$ 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]):

$$\begin{array}{l}{\mu }_{nf}=\frac{{\mu }_f}{{(1-\varphi )}^{2.5}},{\alpha }_{nf}=\frac{k_{nf}}{{\left(\rho C_p\right)}_{nf}},{\rho }_{nf}=\left(1-\varphi \right){\rho }_f+\varphi {\rho }_s,\\ {\left(\rho \beta \right)}_{nf}=\left(1-\varphi \right){\left(\rho \beta \right)}_f+\varphi {\left(\rho \beta \right)}_s,{\left(\rho C_p\right)}_{nf}=\left(1-\varphi \right){\left(\rho C_p\right)}_f+\varphi {\left(\rho C_p\right)}_s,\\ \frac{{\sigma }_{nf}}{{\sigma }_f}=1+\frac{3\left({\sigma }_s/{\sigma }_f-1\right)\varphi }{{\sigma }_s/{\sigma }_f+2-\left({\sigma }_s/{\sigma }_f-1\right)\varphi },\frac{k_{nf}}{k_f}=\frac{k_s+2k_f-2\varphi \left(k_f-k_s\right)}{k_s+2k_f+\varphi \left(k_f-k_s\right)},\end{array}$$

(6)

where ${\mu }_{nf}$ is the dynamic viscosity of the nanofluid, as described by the Brinkman model, ${\mu }_f$ denotes the dynamic viscosity of the base fluid, $\varphi$ denotes the solid volume fraction of nanoparticles, ${\rho }_f$ and ${\rho }_s$ are the density of the base fluid and the solid nanoparticle, respectively, ${\left(\rho C_p\right)}_{nf}$ is the heat capacitance of the nanofluid, $k_{nf}$ 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:

$$\begin{array}{l}{u}_w=v_w=a{\left(x+y\right)}^n,w_w=-\frac{n+1}{2}\sqrt{a{\nu }_f}{\left(x+y\right)}^{\left(n-1\right)/2}s,T_w=T_{\infty }+T_0{\left(x+y\right)}^{2n-1},\\ B=B_0{\left(x+y\right)}^{\left(n-1\right)/2},N_1=N_0{\left(x+y\right)}^{\left(1-n\right)/2},K_1=K_0{\left(x+y\right)}^{\left(1-n\right)/2}.\end{array}$$

(7)

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 $\theta ,$ which are related to the similarity variable $\eta ,$ as follows:

$$\begin{array}{l}u=a{\left(x+y\right)}^n{f}^{\prime }\left(\eta \right),v=a{\left(x+y\right)}^n{h}^{\prime }\left(\eta \right),\\ w=-\sqrt{a{\nu }_f}{\left(x+y\right)}^{\left(n-1\right)/2}\left[\frac{n+1}{2}\left(f+h\right)+\frac{n-1}{2}\eta \left(f^{\prime }+h^{\prime }\right)\right],\\ \theta \left(\eta \right)=\frac{T-T_{\infty }}{T_w-T_{\infty }},\eta =\sqrt{\frac{a}{{\nu }_f}}{\left(x+y\right)}^{\left(n-1\right)/2}z.\end{array}$$

(8)

Note that f and h are the dimensionless quantities, $f^{\prime }$ and $h^{\prime }$ are the dimensionless velocity, $\theta$ is the dimensionless temperature and prime denotes differentiation with respect to $\eta .$ Substituting similarity transformation (8) into governing Equations (1)–(4), the following equations are obtained:

$$\begin{array}{l}\frac{1}{{(1-\varphi )}^{2.5}}{f}^{‴}+\left(1-\varphi +\varphi \frac{{\rho }_s}{{\rho }_f}\right)\left[\frac{n+1}{2}\left(f+h\right)f^{″}-n\left(f^{\prime }+h^{\prime }\right)f^{\prime }\right]\\ -\left[1+\frac{3\left({\sigma }_s/{\sigma }_f-1\right)\varphi }{{\sigma }_s/{\sigma }_f+2-\left({\sigma }_s/{\sigma }_f-1\right)\varphi }\right]M{f}^{\prime }+\left[1-\varphi +\varphi \frac{{\left(\rho \beta \right)}_s}{{\left(\rho \beta \right)}_f}\right]\lambda \theta =0,\end{array}$$

(9)

$$\begin{array}{l}\frac{1}{{(1-\varphi )}^{2.5}}{h}^{‴}+\left(1-\varphi +\varphi \frac{{\rho }_s}{{\rho }_f}\right)\left[\frac{n+1}{2}\left(f+h\right)h^{″}-n\left(f^{\prime }+h^{\prime }\right)h^{\prime }\right]\\ -\left[1+\frac{3\left({\sigma }_s/{\sigma }_f-1\right)\varphi }{{\sigma }_s/{\sigma }_f+2-\left({\sigma }_s/{\sigma }_f-1\right)\varphi }\right]M{h}^{\prime }=0,\end{array}$$

(10)

$$\frac{k_{nf}}{k_f}\left(\frac{1}{\mathrm{Pr}}\right){\theta }^{″}+\left[1-\varphi +\varphi \frac{{\left(\rho C_p\right)}_s}{{\left(\rho C_p\right)}_f}\right]\left[\frac{n+1}{2}\left(f+h\right){\theta }^{\prime }-\left(2n-1\right)\left(f^{\prime }+h^{\prime }\right)\theta \right]=0,$$

(11)

with the following boundary conditions:

$$\begin{array}{l}f\left(0\right)=s,h\left(0\right)=0,\theta \left(0\right)=1+K{\theta }^{\prime }\left(0\right),\\ f^{\prime }\left(0\right)=c+N{f}^{″}\left(0\right),h^{\prime }\left(0\right)=c+N{h}^{″}\left(0\right),\\ f^{\prime }\left(\eta \right)=0,h^{\prime }\left(\eta \right)=0,\theta \left(\eta \right)=0\mathrm{as}\eta \to \infty .\end{array}$$

(12)

In fact, throughout this paper, the governing parameters are as follows: Pr refers to the Prandtl number, $\lambda$ denotes the mixed convection parameter with $\lambda <0$ corresponding to opposing flow and $\lambda >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:

$$\mathrm{Pr}=\frac{{\nu }_f}{{\alpha }_f},\lambda =\frac{G{r}_x}{{\mathrm{Re}}_x^2},M=\frac{{\sigma }_f{B}_0^2}{a{\rho }_f},N=N_0\sqrt{\frac{a}{{\nu }_f}},K=K_0\sqrt{\frac{a}{{\nu }_f}}.$$

(13)

Here, the local Grashof number, $G{r}_x$ and the local Reynolds number, ${\mathrm{Re}}_x$ along x-direction, given by the following:

$$G{r}_x=\frac{g{\beta }_f\left(T_w-T_{\infty }\right){\left(x+y\right)}^3}{{\nu }_f^2},{\mathrm{Re}}_x=\frac{u_w\left(x+y\right)}{{\nu }_f}.$$

(14)

The important physical quantities in this study are the skin friction coefficients, $C_{fx}$ and $C_{fy}$ and the local Nusselt number, $N{u}_x$ which are defined as follows:

$$C_{fx}=\frac{{\tau }_{zx}}{{\rho }_f{u}_w^2},C_{fy}=\frac{{\tau }_{zy}}{{\rho }_f{v}_w^2},N{u}_x=\frac{\left(x+y\right)q_w}{k_f\left(T_w-T_{\infty }\right)},$$

(15)

where ${\tau }_{zx}$ and ${\tau }_{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:

$${\tau }_{zx}={\mu }_{nf}{\left(\frac{\partial u}{\partial z}+\frac{\partial w}{\partial x}\right)}_{z=0},{\tau }_{zy}={\mu }_{nf}{\left(\frac{\partial v}{\partial z}+\frac{\partial w}{\partial y}\right)}_{z=0},q_w=-k_{nf}{\left(\frac{\partial T}{\partial z}\right)}_{z=0}.$$

(16)

Using (7) in (15) and (16), the following are obtained:

$${\mathrm{Re}}_x^{1/2}{C}_{fx}=\frac{1}{{\left(1-\varphi \right)}^{2.5}}{f}^{″}\left(0\right),{\mathrm{Re}}_y^{1/2}{C}_{fy}=\frac{1}{{\left(1-\varphi \right)}^{2.5}}{h}^{″}\left(0\right),{\mathrm{Re}}_x^{-1/2}N{u}_x=-\frac{k_{nf}}{k_f}{\theta }^{\prime }\left(0\right),$$

(17)

where ${\mathrm{Re}}_x=u_w\left(x+y\right)/{\nu }_f$ is the local Reynolds number along x-direction and ${\mathrm{Re}}_y=v_w\left(x+y\right)/{\nu }_f$ is the local Reynolds number along y-direction.

3. Stability Analysis

The present numerical results depict the existence of dual solutions in a certain range of $\lambda$, 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:

$$\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}+w\frac{\partial u}{\partial z}={\nu }_{nf}\frac{{\partial }^{2}u}{\partial z^2}-\frac{{\sigma }_{nf}{B}^{2}u}{{\rho }_{nf}}+\frac{g{\left(\rho \beta \right)}_{nf}}{{\rho }_{nf}}\left(T-T_{\infty }\right),$$

(18)

$$\frac{\partial v}{\partial t}+u\frac{\partial v}{\partial x}+v\frac{\partial v}{\partial y}+w\frac{\partial v}{\partial z}={\nu }_{nf}\frac{{\partial }^{2}v}{\partial z^2}-\frac{{\sigma }_{nf}{B}^{2}v}{{\rho }_{nf}},$$

(19)

$$\frac{\partial T}{\partial t}+u\frac{\partial T}{\partial x}+v\frac{\partial T}{\partial y}+w\frac{\partial T}{\partial z}={\alpha }_{nf}\frac{{\partial }^{2}T}{\partial z^2},$$

(20)

where t denotes the time. Additionally, the new dimensionless functions of u, v, w, and $\theta ,$ are initiated in relation to the similarity variable, $\eta$ and the second similarity variable, $\tau$ for the unsteady-state problem (18)–(20) as follows:

$$\begin{array}{l}u=a{\left(x+y\right)}^n\frac{\partial f}{\partial \eta }\left(\eta ,\tau \right),v=a{\left(x+y\right)}^n\frac{\partial h}{\partial \eta }\left(\eta ,\tau \right),\\ w=-\sqrt{a{\nu }_f}{\left(x+y\right)}^{(n-1)/2}\left[\frac{n+1}{2}\left(f+h\right)+\frac{n-1}{2}\eta \left(\frac{\partial f}{\partial \eta }+\frac{\partial h}{\partial \eta }\right)+\left(n-1\right)\tau \left(\frac{\partial f}{\partial \tau }+\frac{\partial h}{\partial \tau }\right)\right],\\ \theta \left(\eta ,\tau \right)=\frac{T-T_{\infty }}{T_w-T_{\infty }},\eta =\sqrt{\frac{a}{{\nu }_f}}{\left(x+y\right)}^{\left(n-1\right)/2}z,\tau =a{\left(x+y\right)}^{n-1}t.\end{array}$$

(21)

Thus, the unsteady-state problem (18)–(20) can be written as:

$$\begin{array}{l}\frac{1}{{\left(1-\varphi \right)}^{2.5}}\frac{{\partial }^{3}f}{\partial {\eta }^3}+\left(1-\varphi +\varphi \frac{{\rho }_s}{{\rho }_f}\right)[\frac{n+1}{2}\left(f+h\right)\frac{{\partial }^{2}f}{\partial {\eta }^2}-n\left(\frac{\partial f}{\partial \eta }+\frac{\partial h}{\partial \eta }\right)\frac{\partial f}{\partial \eta }\\ -\left(n-1\right)\tau \left(\frac{\partial f}{\partial \eta }+\frac{\partial h}{\partial \eta }\right)\frac{{\partial }^{2}f}{\partial \eta \partial \tau }+\left(n-1\right)\tau \left(\frac{\partial f}{\partial \tau }+\frac{\partial h}{\partial \tau }\right)\frac{{\partial }^{2}f}{\partial {\eta }^2}-\frac{{\partial }^{2}f}{\partial \eta \partial \tau }]\\ -\left[1+\frac{3\left({\sigma }_s/{\sigma }_f-1\right)\varphi }{{\sigma }_s/{\sigma }_f+2-\left({\sigma }_s/{\sigma }_f-1\right)\varphi }\right]M\frac{\partial f}{\partial \eta }+\left(1-\varphi +\varphi \frac{{\left(\rho \beta \right)}_s}{{\left(\rho \beta \right)}_f}\right)\lambda \theta =0,\end{array}$$

(22)

$$\begin{array}{l}\frac{1}{{\left(1-\varphi \right)}^{2.5}}\frac{{\partial }^{3}h}{\partial {\eta }^3}+\left(1-\varphi +\varphi \frac{{\rho }_s}{{\rho }_f}\right)[\frac{n+1}{2}\left(f+h\right)\frac{{\partial }^{2}h}{\partial {\eta }^2}-n\left(\frac{\partial f}{\partial \eta }+\frac{\partial h}{\partial \eta }\right)\frac{\partial h}{\partial \eta }\\ -\left(n-1\right)\tau \left(\frac{\partial f}{\partial \eta }+\frac{\partial h}{\partial \eta }\right)\frac{{\partial }^{2}h}{\partial \eta \partial \tau }+\left(n-1\right)\tau \left(\frac{\partial f}{\partial \tau }+\frac{\partial h}{\partial \tau }\right)\frac{{\partial }^{2}h}{\partial {\eta }^2}-\frac{{\partial }^{2}h}{\partial \eta \partial \tau }]\\ -\left[1+\frac{3\left({\sigma }_s/{\sigma }_f-1\right)\varphi }{{\sigma }_s/{\sigma }_f+2-\left({\sigma }_s/{\sigma }_f-1\right)\varphi }\right]M\frac{\partial h}{\partial \eta }=0,\end{array}$$

(23)

$$\begin{array}{l}\frac{1}{\mathrm{Pr}}\frac{k_{nf}}{k_f}\frac{{\partial }^2\theta }{\partial {\eta }^2}+\left(1-\varphi +\varphi \frac{{\left(\rho C_p\right)}_{nf}}{{\left(\rho C_p\right)}_f}\right)[\frac{n+1}{2}\left(f+h\right)\frac{\partial \theta }{\partial \eta }-\left(2n-1\right)\left(\frac{\partial f}{\partial \eta }+\frac{\partial h}{\partial \eta }\right)\theta \\ -\left(n-1\right)\tau \left(\frac{\partial f}{\partial \eta }+\frac{\partial h}{\partial \eta }\right)\frac{\partial \theta }{\partial \tau }+\left(n-1\right)\tau \left(\frac{\partial f}{\partial \tau }+\frac{\partial h}{\partial \tau }\right)\frac{\partial \theta }{\partial \eta }-\frac{\partial \theta }{\partial \tau }]=0,\end{array}$$

(24)

subject to the boundary conditions:

$$\begin{array}{l}f\left(0,\tau \right)=s,h\left(0,\tau \right)=0,\theta \left(0,\tau \right)=1+K\frac{\partial \theta }{\partial \eta }\left(0,\tau \right),\\ \frac{\partial f}{\partial \eta }\left(0,\tau \right)=c+N\frac{{\partial }^{2}f}{\partial {\eta }^2}\left(0,\tau \right),\frac{\partial h}{\partial \eta }\left(0,\tau \right)=c+N\frac{{\partial }^{2}h}{\partial {\eta }^2}\left(0,\tau \right),\\ \frac{\partial f}{\partial \eta }\left(\infty ,\tau \right)\to 0,\frac{\partial h}{\partial \eta }\left(\infty ,\tau \right)\to 0,\theta \left(\infty ,\tau \right)\to 0.\end{array}$$

(25)

Furthermore, in order to test the stability of solutions $f\left(\eta \right)=f_0\left(\eta \right),$ $h\left(\eta \right)=h_0\left(\eta \right)$, and $\theta \left(\eta \right)={\theta }_0\left(\eta \right)$ in satisfying the steady-state problem (9)–(12), the following is written:

$$\begin{array}{l}f\left(\eta ,\tau \right)=f_0\left(\eta \right)+e^{-\gamma \tau }F\left(\eta ,\tau \right),\\ h\left(\eta ,\tau \right)=h_0\left(\eta \right)+e^{-\gamma \tau }H\left(\eta ,\tau \right),\\ \theta \left(\eta ,\tau \right)={\theta }_0\left(\eta \right)+e^{-\gamma \tau }G\left(\eta ,\tau \right),\end{array}$$

(26)

where $\gamma$ is an unknown eigenvalue parameter. Besides, note that $F\left(\eta ,\tau \right)$, $H\left(\eta ,\tau \right)$, and $G\left(\eta ,\tau \right)$ are smaller relative to $f_0\left(\eta \right),$ $h_0\left(\eta \right)$, and ${\theta }_0\left(\eta \right).$ Furthermore, by substituting (26) into Equations (22)–(24), along with the boundary conditions (25), the following equations are obtained:

$$\begin{array}{l}\frac{1}{{\left(1-\varphi \right)}^{2.5}}\frac{{\partial }^{3}F}{\partial {\eta }^3}+\left(1-\varphi +\varphi \frac{{\rho }_s}{{\rho }_f}\right)[\frac{n+1}{2}\left(f_0\frac{{\partial }^{2}F}{\partial {\eta }^2}+h_0\frac{{\partial }^{2}F}{\partial {\eta }^2}+f_0^{″}F+f_0^{″}H\right)\\ -n\left(f_0^{\prime }\frac{\partial H}{\partial \eta }+2f_0^{\prime }\frac{\partial F}{\partial \eta }+h_0^{\prime }\frac{\partial F}{\partial \eta }\right)+\left(n-1\right)\tau \left(\gamma f_0^{\prime }\frac{\partial F}{\partial \eta }+\gamma h_0^{\prime }\frac{\partial F}{\partial \eta }-f_0^{\prime }\frac{{\partial }^{2}F}{\partial \eta \partial \tau }-h_0^{\prime }\frac{{\partial }^{2}F}{\partial \eta \partial \tau }\right)\\ -\left(n-1\right)\tau \left(\gamma f_0^{″}F+\gamma f_0^{″}H-f_0^{″}\frac{\partial F}{\partial \tau }-f_0^{″}\frac{\partial H}{\partial \tau }\right)+\gamma \frac{\partial F}{\partial \eta }-\frac{{\partial }^{2}F}{\partial \eta \partial \tau }]\\ -\left[1+\frac{3\left({\sigma }_s/{\sigma }_f-1\right)\varphi }{{\sigma }_s/{\sigma }_f+2-\left({\sigma }_s/{\sigma }_f-1\right)\varphi }\right]M\frac{\partial F}{\partial \eta }+\left(1-\varphi +\varphi \frac{{\left(\rho \beta \right)}_s}{{\left(\rho \beta \right)}_f}\right)\lambda G=0,\end{array}$$

(27)

$$\begin{array}{l}\frac{1}{{\left(1-\varphi \right)}^{2.5}}\frac{{\partial }^{3}H}{\partial {\eta }^3}+\left(1-\varphi +\varphi \frac{{\rho }_s}{{\rho }_f}\right)[\frac{n+1}{2}\left(f_0\frac{{\partial }^{2}H}{\partial {\eta }^2}+h_0\frac{{\partial }^{2}H}{\partial {\eta }^2}+h_0^{″}F+h_0^{″}H\right)\\ -n\left(f_0^{\prime }\frac{\partial H}{\partial \eta }+2h_0^{\prime }\frac{\partial H}{\partial \eta }+h_0^{\prime }\frac{\partial F}{\partial \eta }\right)+\left(n-1\right)\tau \left(\gamma f_0^{\prime }\frac{\partial H}{\partial \eta }+\gamma h_0^{\prime }\frac{\partial H}{\partial \eta }-f_0^{\prime }\frac{{\partial }^{2}H}{\partial \eta \partial \tau }-h_0^{\prime }\frac{{\partial }^{2}H}{\partial \eta \partial \tau }\right)\\ -\left(n-1\right)\tau \left(\gamma h_0^{″}F+\gamma h_0^{″}H-h_0^{″}\frac{\partial F}{\partial \tau }-h_0^{″}\frac{\partial H}{\partial \tau }\right)+\gamma \frac{\partial H}{\partial \eta }-\frac{{\partial }^{2}H}{\partial \eta \partial \tau }]\\ -\left[1+\frac{3\left({\sigma }_s/{\sigma }_f-1\right)\varphi }{{\sigma }_s/{\sigma }_f+2-\left({\sigma }_s/{\sigma }_f-1\right)\varphi }\right]M\frac{\partial H}{\partial \eta }=0,\end{array}$$

(28)

$$\begin{array}{l}\frac{1}{\mathrm{Pr}}\frac{k_{nf}}{k_f}\frac{{\partial }^{2}G}{\partial {\eta }^2}+\left(1-\varphi +\varphi \frac{{\left(\rho C_p\right)}_{nf}}{{\left(\rho C_p\right)}_f}\right)[\frac{n+1}{2}\left(f_0\frac{\partial G}{\partial \eta }+h_0\frac{\partial G}{\partial \eta }+{\theta }_0^{\prime }F+{\theta }_0^{\prime }H\right)\\ -\left(2n-1\right)\left(f_0^{\prime }G+h_0^{\prime }G+{\theta }_0\frac{\partial F}{\partial \eta }+{\theta }_0\frac{\partial H}{\partial \eta }\right)+\left(n-1\right)\tau \left(\gamma f_0^{\prime }G+\gamma h_0^{\prime }G-f_0^{\prime }\frac{\partial G}{\partial \tau }-h_0^{\prime }\frac{\partial G}{\partial \tau }\right)\\ -\left(n-1\right)\tau \left(\gamma {\theta }_0^{\prime }F+\gamma {\theta }_0^{\prime }H-{\theta }_0^{\prime }\frac{\partial F}{\partial \tau }-{\theta }_0^{\prime }\frac{\partial H}{\partial \tau }\right)+\gamma G-\frac{\partial G}{\partial \tau }]=0,\end{array}$$

(29)

subject to the boundary conditions:

$$\begin{array}{l}F\left(0,\tau \right)=0,H\left(0,\tau \right)=0,G\left(0,\tau \right)=K\frac{\partial G}{\partial \eta }\left(0,\tau \right),\\ \frac{\partial F}{\partial \eta }\left(0,\tau \right)=N\frac{{\partial }^{2}F}{\partial {\eta }^2}\left(0,\tau \right),\frac{\partial H}{\partial \eta }\left(0,\tau \right)=N\frac{{\partial }^{2}h}{\partial {\eta }^2}\left(0,\tau \right),\\ \frac{\partial F}{\partial \eta }\left(\infty ,\tau \right)\to 0,\frac{\partial H}{\partial \eta }\left(\infty ,\tau \right)\to 0,G\left(\infty ,\tau \right)\to 0.\end{array}$$

(30)

Furthermore, as depicted by Weidman et al. [46], in order to identify the initial growth or the decay of solution (26), $\tau =0,$ is applied and hence, the following are obtained: $F=F_0\left(\eta \right),$ $H=H_0\left(\eta \right)$, and $G=G_0\left(\eta \right).$ Furthermore, in order to test the numerical procedure, the following linear eigenvalue problem has to be solved:

$$\begin{array}{l}\frac{1}{{\left(1-\varphi \right)}^{2.5}}{F}_0^{‴}+\left(1-\varphi +\varphi \frac{{\rho }_s}{{\rho }_f}\right)[\frac{n+1}{2}\left(f_0{F}_0^{″}+h_0{F}_0^{″}+f_0^{″}{F}_0+f_0^{″}{H}_0\right)\\ -n\left(2f_0^{\prime }{F}_0^{\prime }+f_0^{\prime }{H}_0^{\prime }+h_0^{\prime }{F}_0^{\prime }\right)+\gamma F_0^{\prime }]+\left(1-\varphi +\varphi \frac{{\left(\rho \beta \right)}_s}{{\left(\rho \beta \right)}_f}\right)\lambda G_0\\ -\left[1+\frac{3\left({\sigma }_s/{\sigma }_f-1\right)\varphi }{{\sigma }_s/{\sigma }_f+2-\left({\sigma }_s/{\sigma }_f-1\right)\varphi }\right]M{F}_0^{\prime }=0,\end{array}$$

(31)

$$\begin{array}{l}\frac{1}{{\left(1-\varphi \right)}^{2.5}}{H}_0^{‴}+\left(1-\varphi +\varphi \frac{{\rho }_s}{{\rho }_f}\right)[\frac{n+1}{2}\left(f_0{H}_0^{″}+h_0{H}_0^{″}+h_0^{″}{F}_0+h_0^{″}{H}_0\right)\\ -n\left(2f_0^{\prime }{H}_0^{\prime }+f_0^{\prime }{F}_0^{\prime }+h_0^{\prime }{H}_0^{\prime }\right)+\gamma H_0^{\prime }]-\left[1+\frac{3\left({\sigma }_s/{\sigma }_f-1\right)\varphi }{{\sigma }_s/{\sigma }_f+2-\left({\sigma }_s/{\sigma }_f-1\right)\varphi }\right]M{H}_0^{\prime }=0,\end{array}$$

(32)

$$\begin{array}{l}\frac{1}{\mathrm{Pr}}\frac{k_{nf}}{k_f}{G}_0^{″}+\left(1-\varphi +\varphi \frac{{\left(\rho C_p\right)}_{nf}}{{\left(\rho C_p\right)}_f}\right)[\frac{n+1}{2}\left(f_0{G}_0^{\prime }+h_0{G}_0^{\prime }+{\theta }_0^{\prime }{F}_0+{\theta }_0^{\prime }{H}_0\right)\\ -\left(2n-1\right)\left(f_0^{\prime }{G}_0+h_0^{\prime }{G}_0+{\theta }_0{F}_0^{\prime }+{\theta }_0{H}_0^{\prime }\right)+\gamma G_0]=0,\end{array}$$

(33)

subject to the boundary conditions:

$$\begin{array}{l}{F}_0\left(0\right)=0,H_0\left(0\right)=0,F_0^{\prime }\left(0\right)=N{F}_0^{″}\left(0\right),H_0^{\prime }\left(0\right)=N{H}_0^{″}\left(0\right),G_0\left(0\right)=K{G}_0^{\prime }\left(0\right),\\ F_0^{\prime }\left(\eta \right)\to 0,H_0^{\prime }\left(\eta \right)\to 0,G_0\left(\eta \right)\to 0\mathrm{as}\eta \to \infty .\end{array}$$

(34)

The solutions $f_0\left(\eta \right),$ $h_0\left(\eta \right)$, and ${\theta }_0\left(\eta \right)$ 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\left(\eta \right),$ $H_0\left(\eta \right)$, or $G_0\left(\eta \right)$ to better find the range of $\gamma .$ Thus, the condition $F_0^{\prime }\left(\eta \right)\to 0$ as $\eta \to \infty$ is relaxed and for a fixed value of $\gamma$, the linear eigenvalue problem (31)–(34) are solved, along with the new boundary condition; $F_0^{″}\left(0\right)=1.$ The obtained solutions offers an infinite set of possible eigenvalues ${\gamma }_1<{\gamma }_2<{\gamma }_3<\dots$, where ${\gamma }_1$ is the minimum eigenvalue. Furthermore, if ${\gamma }_1>0$, then there is an initial decay of disturbances and the flow becomes stable. Other than that, if ${\gamma }_1>0$, then there is an initial growth of disturbances and the flow becomes unstable.

4. 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⁻⁵. A suitable finite value of $\eta \to \infty$ had been opted, namely $\eta ={\eta }_{\infty }=10$ for the upper and lower branch solutions. The numerical results are retrieved for varied values that involved dimensionless parameters, namely mixed convection parameter $\lambda$, 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, $\varphi ,$ c and s are n = 3, $\varphi$ = 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 $\lambda$, 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. Numerical values of $f^{″}\left(0\right)$ for several values of M and n.

M	n	Mahanthesh et al. [35]	Jain and Choudhary [49]	Present
0	1	−1.41421	−1.419135111	−1.414214
	3	−2.29719	−2.301331346	−2.297186
2	1	–	−2.000000000	−2.000000
	3	–	−2.701221616	−2.701216

presents the comparative studies of the present results for reduced skin friction coefficient $f^{″}\left(0\right)$ for several values of M and n, when c = 1 and $\lambda$ = s = $\varphi$ = 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 ${\gamma }_1$, for several values of $\lambda$, N, K, and M, are given in

Table 3. Minimum eigenvalue, ${\gamma }_1$ at several values of N, K, M, and $\lambda$.

N	K	M	$\mathit{\lambda }$	${\mathit{\gamma }}_1$ (Upper Branch)	${\mathit{\gamma }}_1$ (Lower Branch)
1	1	1	−9	3.4985	−0.2172
			−9.2	1.8120	−0.1051
1	1	0.8	−6	5.0945	−0.1339
			−6.1	3.3007	−0.0133
1	0.5	1	−5	5.2992	−0.4028
			−5.5	1.8051	−0.1048
0.9	1	1	−5	8.9037	−0.4313
			−5.9	5.1995	−0.1139

Note: M is the magnetic parameter, n is the positive constant, N is the velocity slip parameter, K is the temperature slip parameter, $\lambda$ is the mixed convection parameter.

. The results in Table 3 clearly show that ${\gamma }_1>0$ is for upper branch solution, while ${\gamma }_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.

Additionally, Figure 2, Figure 3, Figure 4, Figure 5, Figure 6, Figure 7, Figure 8, Figure 9, Figure 10, Figure 11, Figure 12, Figure 13, Figure 14, Figure 15, Figure 16, Figure 17, Figure 18 and Figure 19 display the variations of the reduced skin friction coefficients $f^{″}\left(0\right),h^{″}\left(0\right)$ and the reduced local Nusselt number $-{\theta }^{\prime }\left(0\right)$ with $\lambda$ 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 $\lambda c<\lambda ,$ where $\lambda c$ is the critical value of parameter $\lambda$, while no solution exists for $\lambda <\lambda c$. Besides, interestingly, the unique (upper branch) solution is attained when $\lambda >0$ (assisting flow), while the dual solutions are found only for $\lambda <0$ (opposing flow). Furthermore, for $\lambda <\lambda c,$ both full Navier-Stokes and energy equations have to be solved. In addition, it is clear that all $\left|\lambda c\right|$ 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, Figure 2, Figure 3, Figure 4, Figure 5, Figure 6 and Figure 7 illustrate the variations of $f^{″}\left(0\right),h^{″}\left(0\right)$, and $-{\theta }^{\prime }\left(0\right)$ against $\lambda$ 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 $\left|\lambda c\right|$ 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 $\lambda$ until it reaches a critical value at $\lambda c$, for both opposing and assisting flows. Besides, the corresponding lower branch solution decreases with the increase of $\lambda$ for an opposing flow. Nonetheless, the values of $-{\theta }^{\prime }\left(0\right)$ (see Figure 6 and Figure 7) for the lower branch solution decreased rapidly with the increase of $\lambda$. Moreover, the values of $-{\theta }^{\prime }\left(0\right)$ for the upper branch solution are extremely large, when compared to the corresponding lower branch solution in the neighborhood of $\lambda =0.$ Furthermore, Figure 2, Figure 3, Figure 4 and Figure 5 demonstrate that at any $\lambda$ station within the range of $\lambda c<\lambda ,$ the values of $f^{″}\left(0\right)$ and $h^{″}\left(0\right)$ for the dual solutions decrease as N increases. In contrast, Figure 6 and Figure 7 display that at any $\lambda$ station in the range of $\lambda c<\lambda <0,$ the values of $-{\theta }^{\prime }\left(0\right)$ for the dual solutions escalate with the increase in N.

Meanwhile, in Figure 8, Figure 9, Figure 10, Figure 11, Figure 12 and Figure 13, the variations of $f^{″}\left(0\right)$, $h^{″}\left(0\right)$, and $-{\theta }^{\prime }\left(0\right)$ are plotted with $\lambda$ for several values of temperature slip parameter K for both nanoparticles in water–based nanofluids. As noted from the figures, the critical value $\left|\lambda c\right|$ 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, Figure 8, Figure 9, Figure 10 and Figure 11 point out that the values of $f^{″}\left(0\right)$ and $h^{″}\left(0\right)$ for the upper branch solution increase with increment in K within the domain of $\lambda c<\lambda <0.$ On the other hand, the opposite trend is observed when $\lambda >0.$ Additionally, Figure 12 and Figure 13 portray that at any $\lambda$ within the range of $\lambda c<\lambda <0,$ the values of $-{\theta }^{\prime }\left(0\right)$ for both solution branches increase as K increases.

Furthermore, Figure 14, Figure 15, Figure 16, Figure 17, Figure 18 and Figure 19 explain the impact of magnetic parameter M on $f^{″}\left(0\right)$, $h^{″}\left(0\right)$ and $-{\theta }^{\prime }\left(0\right)$. In addition, these figures indicate that the strength of $\left|\lambda c\right|$ 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 $\lambda$ within the domain of $\lambda <\lambda c.$ The corresponding lower branch solution decreases as $\lambda$ is increased when $\lambda c<\lambda <0.$ Nonetheless, Figure 18 and Figure 19 reveal that the values of $-{\theta }^{\prime }\left(0\right)$ for the lower branch solution display rapid decrease with the increase of $\lambda$. Moreover, the values of $-{\theta }^{\prime }\left(0\right)$ for the lower branch solution are large, in comparison to the corresponding upper branch solution in a region close to $\lambda =0.$ Besides, Figure 14, Figure 15, Figure 18 and Figure 19 illustrate that at any $\lambda$ station within the range of $\lambda c<\lambda ,$ for increment of M, the values of $f^{″}\left(0\right)$ (within the domain of $\lambda c<\lambda$) and $-{\theta }^{\prime }\left(0\right)$ (within the domain of $\lambda c<\lambda <0$) portray an increment in the upper branch solution, as well as a decrease in the lower branch solution. Furthermore, Figure 16 and Figure 17 point out that the values of $h^{″}\left(0\right)$ for the upper branch solution display a similar trend with that of $f^{″}\left(0\right)$. Nevertheless, the trend for the lower branch solution differs, when compared to the trend observed for the upper branch solution.

Finally, Figure 20, Figure 21 and Figure 22 offer sample of dimensionless velocity $f^{\prime }\left(\eta \right)$, $h^{\prime }\left(\eta \right)$, and temperature $\theta \left(\eta \right)$ profiles against $\eta$, 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 Figure 2, Figure 3, Figure 4, Figure 5, Figure 6, Figure 7, Figure 8, Figure 9, Figure 10, Figure 11, Figure 12, Figure 13, Figure 14, Figure 15, Figure 16, Figure 17, Figure 18 and Figure 19. Furthermore, the upper branch solution for these profiles portray a thinner boundary layer, in comparison to that of the lower branch solution. Additionally, Figure 20 and Figure 21 present for a fixed values of n, $\varphi$, c, s, $\lambda$, N, K, and M, the upper branch solution for $f^{\prime }\left(\eta \right)$ and $h^{\prime }\left(\eta \right)$ are higher for Ag-water nanofluid. Other than that, Figure 22 shows that the upper branch solution of $\theta \left(\eta \right)$ 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.

5. 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^{″}\left(0\right)$, $h^{″}\left(0\right)$, and $-{\theta }^{\prime }\left(0\right)$ 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 $\lambda$ in the opposing flow region with the solution curves bifurcating at the critical values of $\lambda c$, while the lower branch solution could not be continued further to the point where $\lambda$ 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 $\left|\lambda c\right|$ 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.
• Ag-water nanofluid is more capable of delaying the boundary layer separation, in comparison to Cu-water nanofluid.
• Ag-water nanofluid displayed better enhancement for heat transfer, when compared to that for Cu-water nanofluid.