1 Introduction

The behavior of nanofluids^[1-4] in various heat transfer processes of different engineering, industrial^[5], and biomedical^[6] applications led us to rethink the development of appropriate mathematical models. The two-component nanofluid model based on seven slip mechanisms suggested by Buongiorno^[7] has attracted attention among various researchers and highlighted the importance of Brownian motion of nanoparticle and thermophoretic velocity due to temperature gradient. The problems of flow and heat transfer over stretching^[8] and shrinking^[9-11] sheet due to their various practical applications in plastic and wire industries, paper technology, polymer engineering, and polyethylene industries have been studied by many researchers. Khan and Pop^[12] and Rana and Bhargava^[13] extended these problems to nanofluid flow and heat transfer with linear and nonlinear stretching models, respectively. Kuznetsov and Nield^[14] studied natural convective boundary layer flow of a nanofluid past a vertical plate. The steady state problems on nanofluid have also been extended to magneto-hydrodynamics^[15-17], unsteady solution^{[16, 18-20]}, and non-Newtonian nanofluids^{[17, 21-24]} by various scientists and researchers.

In recent studies, many authors have found multiple solutions for different flow models. Earlier, Merkin^[25], Weidman et al.^[26], and Harris et al.^[27] found the two solutions in their respective studies. Stability analysis has been used on the basis of eigenvalues which predict physically realizable solutions. If the minimum eigenvalue is positive, then the solution is stable or physically possible, otherwise it is unstable. Moreover, the critical values calculated in these analyses tell us the existence of second solution. In recent studies of nanofluids, many authors have already found multiple solutions. Bachok et al.^[28] studied stagnation point nanofluid flow over the stretching/shrinking sheet using a single component based nanofluid transport model. They concluded that the critical values for the shrinking sheet ($\varepsilon \approx -1.246 5$) remain invariant with the inclusion of Cu, Al$_{2}$O$_{3}$, and TiO$_{2}$ nanoparticles. Later, the same authors extended their work for the exponentially stretching/shrinking sheet^[29]. Zaimi et al.^[30-31] studied the above problem in nanofluid using two-component Buongiorno' model both for linear and nonlinear stretching/shrinking sheets. The problem of bio-convection with gyrotactic micro-organisms in nanofluids was also studied by the same authors^[32]. Dhanai et al.^[17] extended the problem of Zaimi et al.^[31] using magnetohydrodynamic boundary layer flow and viscous dissipation effect. Later, Dhanai et al.^[23] extended the work to non-Newtonian power-law fluid under slip conditions and heat source/sink and found the critical value of power-law \linebreak index.

Recently, several authors introduced the existence of spin of nanoparticle which led them to extend the two-component nanofluid model on the basis of Eringen theory^[33-34] of micro-rotation which could potentially explain theoretical and experimental results. Gorla and Kuma- ri^[35] investigated the Eringen theory based micropolar nanofluid model for nonlinear stretching sheet. Hussain et al.^[36] implemented the different correlations for silver and copper nanoparticles considering the homogenous model using the Eringen micropolar model. Under these assumpti-ons^[37-39], we extend the paper of micropolar fluid by Yacob and Ishak^[40] to the nanofluid model. The passively controlled boundary conditions on nanoparticle which is coupled with temperature proposed by Kuznetsov and Nield^[41] are also taken in consideration so that we can have a more realistic model. The main focus is on finding the critical values for suction and shrinking parameters by taking multiple-slip effects^[42-43] and, further, the stability analysis is also performed. To the authors' best knowledge, this analysis has not been done before.

2 Mathematical formulation

Consider the flow of viscous incompressible micropolar nanofluid in the region $\overline{y}>0$ driven by a shrinking sheet of uniform temperature $T_{\mathrm w} $. The sheet velocity is assumed to vary linearly with the distance $\overline{x}$ from the fixed point on the sheet, i.e., $U_{\mathrm w} (\overline{x})=a_{1}\chi\overline{x}$, where $a_{1} $ is a positive constant. The flow model for the present problem is described in Fig. 1. The simplified two-dimensional equations governing the flow in the boundary layer of a steady, laminar, and incompressible micropolar nanofluid are^{[7, 44]}

The boundary conditions for the model obtained from Karniadakis et al.^[45] are given as

Here, $\overline{u}$ and $\overline{v}$ are the velocity components along the $\overline{x}$- and $\overline{y}$-axes, $\overline{N}$ is the micro-rotation/ angular velocity, $\rho $ is the fluid density, $U_{\mathrm w}$ and $T_{\mathrm w}$ are the velocity and the temperature of the sheet, respectively, $T_\infty$ is the ambient temperature (less than $T_{\mathrm w}$), $\nu $ is the kinematic viscosity, $T$ is the temperature, $C$ is the nanoparticle volume fraction concentration, $\mu $ is the dynamic viscosity, $j$ is the microinertia density, $\alpha $ is the thermal diffusivity, m is a constant, $\gamma $ is the micropolar spin gradient viscosity, $D_{1} $ is the thermal slip factor, $N_{1} $ is the hydrodynamic slip factor, $\kappa $ is the vortex viscosity, $D_{\mathrm B} $ is the Brownian diffusion coefficient, $D_{\mathrm T} $ is the thermophoretic diffusion coefficient, $\tau =\frac{(\rho c)_{\mathrm p} }{(\rho c)_{\mathrm f} } $ is the ratio of heat capacity of the nanoparticle and fluid, $(\rho c)_{\mathrm f} $ is the heat effective heat capacity of the fluid, and $(\rho c)_{\mathrm p} $ is the effective heat capacity of the nanoparticle material. Further, $V_{\mathrm w} $ is the mass flux through the permeable sheet with $V_{\mathrm w} ＜0$ for suction and $V_{\mathrm w}>0$ for injection. We notice that $m$ is a constant such that $0\le m\le 1$. Following Ahmadi^[46], we assume that $\gamma $ is given by

2.1 Boundary layer equations in dimensionless form

We consider the steady state flow and introducing the following dimensionless variables:

and hence introducing the dimensionless stream function $\psi$ defined by $u=\frac{\partial \psi }{\partial y}$ and $v=-\frac{\partial \psi }{\partial x} $ into Eqs. (2)-(5) , we derive

The boundary conditions become

The dimensionless parameters defined in Eqs. (9) -(13) are as follows. $\chi$ is the shrinking parameter, $K=\frac{\kappa }{\mu } $ is the micropolar parameter, $Pr =\frac{\nu }{\alpha } $ is the Prandtl number, $a=N_{1} \sqrt{a_{1} \nu } $ is the velocity slip parameter, $b=D_{1} \sqrt{\frac{a_{1} }{\nu } } $ is the thermal slip, $N_{\mathrm B}=\frac{\tau D_{\mathrm B} \Delta C}{\alpha } $ is the Brownian motion parameter, $N_{\mathrm T}=\frac{\tau D_{\rm T} \Delta T}{\alpha } $ is the thermophoresis parameter, $Sc=\frac{\nu }{D_{\mathrm B} } $ is the Schmidt number, and $s=-\frac{V_{\mathrm w} }{\sqrt{a_{1} \nu } } $ is the suction/injection parameter.

2.2 Search for similarity transformations using Lie group analysis

The solution of the system of partial differential equations (9) -(12) with the boundary conditions (13) is extremely difficult and computationally expensive. Hence, we transform them into a corresponding ordinary differential equations system using the scaling group of transformations. We define the scaling group of transformations as follows^[47]:

Here, $\varepsilon$ is the parameter of the group $\Gamma $, and $c_{i} (i=1, 2, \cdots, 6) $ are arbitrary real numbers not all zero, simultaneously. Equations (9) -(12) will remain invariant under the group of transformations in Eq. (14) if $c_i$ for different values of $i$ can be related as

Note that $y^{*} =y, ~\theta ^{*} =\theta , ~\phi ^{*} =\phi, $ i.e., $y, ~\theta$, and $\phi$ are invariants. Let us denote them by $\eta =y, ~ \theta =\theta (\eta )$, and $\phi =\phi(\eta)$. The characteristic equation is

Solving Eq. (16) , we obtain

Hence, the similarity transformations are

where $\eta $ is the similarity independent variable, $f(\eta )$ is the dimensionless velocity function, and $h(\eta)$ is the dimensionless angular velocity.

2.3 Similarity equations

Using Eq. (17) , from Eqs. (9) -(12) , we have

The boundary conditions (13) now become

where primes denote ordinary differentiation with respect to $\eta $.

2.4 Quantities of engineering interest

The physical quantities for this study are the dimensionless friction factor and the rate of heat and mass transfer. The non-dimensional friction factor is defined as

where

Using Eqs. (8) and (17) in Eq. (23) , we deduce

We finally obtain

The non-dimensional Nusselt number $Nu_{\overline{x}} $ and the Sherwood number $Sh_{\overline{x}} $ are defined as

where $q_{\mathrm w} $ and $q_{\mathrm m} $ are the wall heat and mass flux due to combined effect of Brownian motion and thermophoresis, respectively. Thus, $q_{\mathrm w} $ and $q_{\mathrm m} $ are described as

where

Using Eqs. (8) and (17) in Eqs. (25) and (26) , we have

In the above expressions, $Re_{\overline{x}} =\frac{U_{\mathrm w} \overline{x}}{\nu } $ is the local Reynolds number. The Sherwood number is calculated as zero due to the revised boundary condition

3 Flow stability analysis

To investigate the stability of the both branches (upper and lower), we consider the unsteady governing equations (1) -(5) . We introduce the following transformations:

Using Eq. (28) , the system of Eqs. (1) -(5) is reduced into the following dimensionless partial differential equations:

subject to the boundary conditions

To test the stability of the steady state solutions $f(\eta )=f_{0} (\eta)$, $h(\eta)=h_{0} (\eta)$, $\theta(\eta)=\theta _{0} (\eta )$, and $\phi (\eta )=\phi _{0} (\eta)$ satisfying the given boundary value problem, we have to consider^{[23, 25-27]}

where $\alpha $ is an unknown eigenvalue parameter, and $F(\eta , \tau )$, $H(\eta , \tau )$, $G(\eta , \tau )$, and $I(\eta , \tau )$ are small relative to $f_{0} (\eta )$, $h_{0} (\eta )$, $\theta _{0} (\eta )$, and $\phi _{0} (\eta )$, respectively. We obtain the following linearized problem:

along with the boundary conditions

The stability of the steady state solutions $f(\eta)=f_{0} (\eta )$, $h(\eta)=h_{0} (\eta)$, $\theta(\eta)=\theta _{0}(\eta)$, and $\phi (\eta)=\phi _{0} (\eta)$ can be done by setting $\tau =0$ and $F=F_{0} (\eta )$, $G=G_{0} (\eta )$, $H=H_{0} (\eta )$, and $I=I_{0} (\eta )$. The following linear eigenvalue problems are obtained:

with the boundary conditions

The stability of physical realizable solution either first solution (solid line) or second solution (dashed line) for the steady flow $f_{0}(\eta )$, $h_{0}(\eta)$, $\theta _{0}(\eta)$, and $\phi _{0} (\eta)$ is determined by the smallest eigenvalue $\alpha $ for fixed values of governing parameters. The range of possible eigenvalues can be obtained by relaxing a suitable boundary condition on $F_{0} (\eta)$, $G_{0} (\eta )$, $H_{0}(\eta)$, or $I_{0} (\eta)$. For the present problem, we relax the condition that $F'_{0}(\eta)\to 0$ as $\eta \to \infty $. Hence, for a fixed value of $\alpha $, we solve the final system (40) -(44) with the new boundary condition $F"_{0}=1$.

4 Results and discussion

In the present analysis, the system of ordinary differential equations (18) -(21) with the boundary conditions (22) is solved using the RKF45 method^{[17, 19, 23, 43, 48-49]}, and the results are validated with the ode solver in MATLAB. The calculations throughout the paper are carried out with the default set of governing parameter as $K$=1, $m$ = 0.5, $Pr$ = 6.2, $Sc$ = 10, $s$ = 3.0, $\chi =-1$, $N_{\mathrm B}=N_{\mathrm T}$ = 0.5, $a=b$ = 0.5, unless otherwise stated. The two solution branches are captured for different sets of controlling parameters, and the critical values ($s_{\mathrm c}, \chi _{\mathrm c} $) are determined for the micropolar nanofluid. The critical value for any parameter tells us the existence of no solution ($-\infty ＜s_{\mathrm c} , \chi _{\mathrm c} $), unique solution ($s_{\mathrm c} , \chi _{\mathrm c} $), and multiple solutions ($s_{\mathrm c}, \chi _{\mathrm c}>\infty $) in its domain. The linear regression analysis is also carried out for finding the overall impact of dominant parameters on the reduced Nusselt number for both no-slip and multiple-slip conditions. Under the no-flux boundary condition as suggested by Kuznetsov and Nield^[41] in their revised publications, we develop the correlations with the set of parametric values $K$=[0.1, 0.6, 1.3, 1.9], [$N_{\mathrm T}$, $N_{\mathrm B}$]=[0.1, 0.2, 0.3, 0.4, 0.5] as follows:

(i) For $a=b=0$ (no-slip)

(ii) For $a=b=0.5$ (multiple-slip)

The maximum errors in the above correlations for no-slip and multiple-slip are calculated to be 1.5% and 0.1%, respectively. Owing to this acceptable small error, the quadratic regression analysis is not required for the sake of brevity. It is clear from the correlations that the thermophoresis parameter is the key parameter to control the heat transfer. The behavior of the micropolar parameter is no longer the same as we go on increasing the slip parameters. With the slip parameters, the heat transfer drastically falls from 18.237 9 to 1.802 1 keeping other parameters zero. Both the Brownian motion and thermophoresis parameters are decreasing functions of the reduced Nusselt number. The contribution of Brownian motion to thermal boundary layer (which is proportional to nanoparticle concentration gradient) now tends to zero as the wall is approached, due to a passively controlled no-flux boundary condition. Moreover, the thermophoresis controls the migration of nanoparticle arises due to the temperature difference which has significant impact on the rate of heat transfer. Some of the results are presented in Table 1.

The multiple solutions (two branches) for the present shrinking sheet problem containing micropolar nanofluids are evaluated numerically. The first solution (upper branch) is easy to capture, while the second solution (lower branch) requires proper choice of the maximum boundary layer thickness during numerical computations. Out of these solutions, the physical realizable solution can only be determined by stability analysis. From the solution of linear eigenvalue problem (30) -(43) with the condition (44) , we construct a table for finding the minimum eigenvalues ($\alpha $) for different values of suction parameter and slip conditions. The positive eigenvalue reveals that the solution is physically possible, whereas the negative eigenvalue corresponds to unrealizable solution. Moreover, the eigenvalues are tending to zero as we are approaching to the critical value of suction in all selected cases (see Tables 2 and 3). Thus, it strongly validates the accuracy of the solution of eigenvalue problem.

In Fig. 2, the skin friction is plotted against the suction parameter for $a=b=0$ and $a=b=0.5$. The solid and dashed lines correspond to upper (stable) and lower branch (unstable), respectively. According to Miklavcic and Wang^[9] and Fang^[11], the flow induced by the shrinking sheet can be maintained by an adequate value of suction. From the solution of flow equations, they reported that only one solution exits at s=2, but the flow will generate into two solutions when $s>2$ and no solution for $s＜2$. We obtain similar results. The result obtained for $K$=0 (no-slip) strongly validates our numerical solution. As we increase the value of K from 0 to 2, the critical value for suction increases from 2.000 0 to 2.827 4. This is because of the negative contribution of rotation of nanoparticles to the velocity of fluid. There is an increase in the skin friction. However, with the consideration of both momentum and thermal slips, the critical values of suction decrease for each case of micropolar parameter ($K$=0, 1, 2) . The momentum slip at the boundary provides extra velocity which in turn decreases the skin friction.

Figures 3 and 4 depict the behavior of the critical value of the shrinking parameter for the skin friction and the Nusselt number, respectively. For the fixed value of suction which is taken as 3, the values of $\chi _{\mathrm c}＜\chi $ are $-$2.252 0, $-$1.150 2, and $-$1.125 5 for $K$=0, 1, and 2, respectively. Beyond that, there is no solution for the given boundary value problem. With the incorporation of slip conditions, the critical values rise to $-$2.971 5, $-$1.812 5, and $-$1.300 6 which reveal the opposite behavior for the shrinking parameter as compared with the suction parameter. At $\chi =-1$ and $s=3$, we can easily interpret that for the upper branch, the skin friction first decreases and then increases with $K$, but for the lower branch, the skin friction decreases with $K$. With the implication of slip conditions, the skin friction continuously increases from $K$ = 0 to 2 for first solution. The Nusselt number is found to decrease for both no-slip and multiple-slip conditions with the increase of the particle spin (micropolar parameter).

In Fig. 5, the velocity profile is established for different values of micropolar parameter. Both the branches upper and lower show opposite trend with the enhancement of particle rotation. The physically realizable upper branch for the velocity profile is detected to be diminishing with the increment of $K$. The slip conditions decrease the shrinking velocity at the boundary, but the pattern of velocity profile is the same as for the no-slip flow and thermal conditions. In Figs. 6-8, the velocity, temperature, and concentration profiles are analyzed under the effect of suction parameter. The upper branch slows down with the increase of suction, but the reverse effect is observed for lower branch. The increase in the suction parameter decreases the thermal boundary layer thickness for both cases. There is negligible difference for both the branches for temperature and concentration profiles. With the no-flux and no-slip conditions, the concentration first boosts up near the sheet which is higher for larger suction parameter, and later after $\eta \approx 0.05$, the concentration also tends to ambient concentration much faster for higher suction values. It should be kept in mind that the pore size is sufficiently large as compared with the nanoparticle diameter to prevent clogging of these particles and necessary surfactant should be added to avoid settling. Finally in Figs. 9(a) and 9(b), the streamlines are plotted for both solutions (upper and lower) with the default set of parameters. For both cases, the flows are symmetric about the $y$-axis. Due to high value of suction, the streamlines slightly bend towards the center for upper branch, and this impact is prominent for the case of lower branch.

5 Concluding remarks

In the present framework, the multiple solutions are captured for the problem of flow over a shrinking sheet. The non-Newtonian micropolar nanofluids are utilized for the present analysis which considers the particle spin of base fluid. With the implementation of physically realistic passively controlled boundary conditions, the Brownian motion parameter becomes insignificant to heat transfer. Multiple regression analysis shows that thermophoresis is now a key parameter to heat transfer. It is found to be a decreasing function of the reduced Nusselt number. Evaluation of critical values and stability analysis for realizable solution is performed.