Dual Solutions and Stability Analysis of Micropolar Nanoﬂuid Flow with Slip E ﬀ ect on Stretching / Shrinking Surfaces

: The purpose of the present paper is to investigate the micropolar nanoﬂuid ﬂow on permeable stretching and shrinking surfaces with the velocity, thermal and concentration slip e ﬀ ects. Furthermore, the thermal radiation e ﬀ ect has also been considered. Boundary layer momentum, angular velocity, heat and mass transfer equations are converted to non-linear ordinary di ﬀ erential equations (ODEs). Then, the obtained ODEs are solved by applying the shooting method and in the results, the dual solutions are obtained in the certain ranges of pertinent parameters in both cases of shrinking and stretching surfaces. Due to the presence of the dual solutions, stability analysis is done and it was found that the ﬁrst solution is stable and physically feasible. The results are also compared with previously published literature and found to be in excellent agreement. Moreover, the obtained results reveal the angular velocity increases in the ﬁrst solution when the value of micropolar parameter increases. The velocity of nanoﬂuid ﬂow decreases in the ﬁrst solution as the velocity slip parameter increases, whereas the temperature proﬁles increase in both solutions when thermal radiation, Brownian motion and the thermophoresis parameters are increased. Concentration proﬁle increases by increasing N t and decreases by increasing N b .


Introduction
The boundary flow and heat transfer phenomena of all types of fluids have remained of great interest for many researchers, especially for purpose of the practical implementation. Originally, fluids are classified into Newtonian and non-Newtonian fluids. There are many industrial applications of non-Newtonian fluids, such as drilling mud, polycrystal melts, oils, certain paints, volcanic lava, cosmetic product, fluid suspensions, molten polymers and food product, etc. There are many mathematical models for various constitutive equations that deal with flow phenomena with different parameters. Among such models, the micropolar fluid model is one of the non-Newtonian models introduced by Eringen [1]. This model deals with the flow of the micropolar fluids, whereas some examples of the micropolar fluids are liquid crystals with rigid molecules, suspensions or colloidal solutions, some biological fluids, exocytic lubricants and the blood of the animals, etc. In many industries, fluids are used for heat transport purposes, but common fluids are weak source of the heat transport. So, with development of nanotechnology, in the field of fluid mechanics, Choi and Eastman [2] introduced a new kind of modern fluid that could suspend the solid nanoparticles into common base fluids, later named as nanofluid. Nanofluids possess high thermal conductivity as compared to all concerned host fluids. Nanofluids is a sub-class of the modern rapidly growing field of nanoscience and technology.
Several practical problems have been discussed by many researchers in order to enhance thermal conductivity for the last many years. There are two well-known models that support transport phenomena of nanofluids, which are the Tiwari and Das model [3] and Buongiorno's model [4]. It should be noted that Buongiorno's model is considered a two-phase non-homogeneous model in which the slip velocity of the base fluid and nanoparticles are not equal to zero. This model contains seven slip parameters, including diffusiophoresis, inertia, fluid drainage, gravity, Brownian motion, Magnus and thermophoresis [5]. On the other hand, the Tiwari-Das model is one of the single-phase (homogeneous) models in which the thermophysical characteristics of the base fluid are improved with the impact of nanoparticles allied with viscosity and thermal conductivity. Many researchers considered these two models to study all types of nanofluids as well as micropolar nanofluids, such as Hsiao [6], Lund et al. [7,8], Hashemi et al. [9], Hussain et al. [10], Dero et al. [11][12][13], Alarifi et al. [14] Pourfattah et al. [15], Alsarraf et al. [16], Jafarimoghaddam et al. [17,18], and Pal and Mandal [19]. Further work on nanofluid is referenced in [20][21][22]. In the present research, Buongiorno's model is considered in order to develop the mathematical model for boundary layer micropolar nanofluid flow with slip parameters and radiation effect. It is observed that a little attention is given to this model for the case of existing multiple solutions in flow of the micropolar nanofluid.
The velocity, temperature and concentration slip factors may exist in case of rising the relative difference among velocity, temperature and the concentration entities. Awais et al. [23] studied the different slip conditions on a stretching sheet. They stated that if the temperature in the fluid flow is controlled by the thermal slip parameter then there is the possibility the concentration may also be controlled through transport phenomena of mass transfer. Das [24] examined the velocity slip on micropolar viscous fluid and found the increasing value of the velocity slip parameter will decrease the velocity as well as the thickness of the boundary layer. Ramya et al. [25] considered the nanofluid with thermal and velocity slip parameters. Some studies of micropolar nanofluid with different slip conditions can be found in [26][27][28][29].
It has been observed in previous literature that the work on the non-Newtonian nanofluid flow is very limited. Therefore, in present research work, a non-Newtonian micropolar nanofluid is considered by using a two-phase model. There study of Hayat et al. [30] examined the two-dimensional incompressible micropolar nanofluid flow over a linearly stretching/shrinking sheet with velocity, thermal and concentration slip effects. To the best of the author's knowledge, the multiple similarity solutions with stability analysis of micropolar nanofluid flow on a permeable stretching and shrinking surface along the velocity, thermal and the concentration slip effects have not been studied by any other researcher in the past years. Thus, the purpose of the present work is to determine the multiple similarity solutions with stability analysis of the micropolar nanofluid flow with thermal radiation and slip effects. Buongiorno's model is used due to its novelty of thermophoresis and Brownian motion of small particles. Furthermore, numerical solutions have been obtained by applying the shooting technique with shootlib function in the Maple software. There has been stability analysis performed in order to determine the stable solution with the help of BVP4C solver. The graphical results of different applied parameters are illustrated in the discussions. Furthermore, the physical behaviors of the skin friction coefficient, and the local Nusselt and Sherwood number are discussed. The main focus of this paper is to determine the multiple similarity solutions with stability analysis of micropolar boundary layer nanofluid and to understand the behavior of different parameters in the flow and heat transfer characteristics. It is expected that the present study will help those researchers interested in multiple similarity solutions with stability analysis in micropolar nanofluid flows.

Flow Analysis
Two-dimensional incompressible micropolar nanofluid flow through a linearly shrinking and stretching surfaces with velocity, thermal and concentration slip effects are shown in Figure 1 Furthermore, the effects of Brownian motion and the thermophoresis have been considered in Newtonian heat and concentration equations with thermal radiation effect. By using the above-mentioned assumptions as well as boundary layer approximations, the system of governing equations of micropolar boundary layer nanofluid flow will be written as: By using similarity transformations (7) in Equations (1)-(5), the following ODEs are obtained, subject to boundary conditions

Subject to boundary conditions
Where u and v represent velocity components in the directions of the x-axis and y-axis, respectively. While ϑ, ρ, κ, N,γ and j are kinematic energy, fluid density, vortex viscosity, angular velocity, spin gradient viscosity and micro-inertia density, respectively. Furthermore, T, α, K * , σ * , τ w , D B , D T , T ∞ and C are temperature, thermal diffusivity of nanofluid, mean absorption coefficient, Stefan-Boltzmann constant, ratio of effective heat capacity of nanoparticles of solid material to the heat capacity of fluid, thermophoretic diffusion coefficient, coefficient of Brownian diffusion, free stream temperature and concentration field, respectively. v 0 is mass flux velocity, u w (x) = ax is shrinking/stretching velocity of surface and λ is shrinking/stretching parameter where λ < 0 indicates shrinking surface while λ > 0 indicates stretching surface and β, β 1 and β 2 are the velocity, heat transfer and the concentration slip parameters, respectively.
The following similarity variables are applied to determine the similarity solutions of the problem.
By using similarity transformations (7) in Equations (1)-(5), the following ODEs are obtained, 1 1 subject to boundary conditions where, prime indicates differentiation with respect to the variable η, K = κ µ is the micropolar material parameter, Pr = ϑ α is a Prandtl number, Rd = is a thermophoresis parameter and Sc = ϑ D B is the Schmidt number. Moreover, f w = − v 0 √ aϑ is the suction parameter ( f w > 0), δ = β a ϑ is a velocity slip, ϑ is a thermal slip and δ C = β 2 a ϑ is concentration slip parameters. It should be noted that the values of m = 0 show that the microelements are closed to surface with a strong concentration, which means the particles are not rotatable. In the same manner, m = 0.5 shows that the concentration is week due to vanishing of the anti-symmetric part of the stress tensor. In addition, m = 1 specifies the turbulent boundary layer flow model of Ramzan et al. [31].
Quantities of the physical interests are the coefficient of the skin friction C f , local Nusselt number N u and Sherwood number S h that are stated as, Using Equation (7) in (13), it is obtained where Re x = ax 2 /ϑ is a local Reynolds number.

Stability Analysis
The numerical computation of the Equations (8)-(11) with boundary conditions (12) shows the occurrence of the dual solutions. Therefore, it is necessary to perform the stability analysis to find the stable and physically feasible solution. For stability of solution, the first step is to convert governing Equations (2)-(5) into unsteady form by introducing a new time dependent variable τ as mentioned by [32][33][34].
Now, the new similarity variables have been introduced as follows By using Equation (19) in Equations (15)-(18), it is obtained Subjected to boundary conditions Now, there has been perturbed the basic solutions where smallest eigenvalue is ε and F(η), G(η), H(η) and S(η) are small relative to f 0 (η), g 0 (η), θ 0 (η) and ∅ 0 (η) respectively. In the last one, the following system of linearized eigenvalue problems was obtained: with boundary conditions The above-mentioned linearized Equations (26)- (29) with boundary conditions (30) are solved by using BVP4c solver function in MATLAB software and found values of smallest eigenvalue ε. In order to find the smallest eigenvalues, there is to be relaxed one of the boundary conditions into initial condition as recommended by Harris et al. [35] and Ali et al. [36]. In this particular problem, there is relaxed the G 0 (η) → 0 as η → ∞ into G 0 (0) = 1. It should be noted that negative smallest eigenvalues indicate the initial growth of disturbance, henceforth the solution of fluid flow becomes unstable. Whereas, if smallest eigenvalue value is positive, that flow of fluid is stable and physically realizable.

Result and Discussion
The governing equations of boundary layer micropolar nanofluid flow [8][9][10][11] subject to boundary conditions (12) are solved by shooting method with help of Maple software. The results of the momentum, angular velocity, temperature and concentration profiles are demonstrated graphically for different values of various physical parameters such as micropolar material, Brownian motion, thermal radiation, thermophoresis, and slip parameters. The obtained results are compared to those obtained by Hayat et al. [30] using homotopy analysis method (HAM) for the values of skin friction coefficients presented in Table 1. It is found that the present results show good agreement with those obtained by Hayat et al. [30]. So, it can be said that the obtained results are correct and reliable. Furthermore, the numerical solutions indicate the occurrence of dual solutions in the micropolar nanofluid flow problem for all profiles. Therefore, another technique, the three-stage Lobatto IIIa formula is created in BVP4c with the help of finite difference technique. Afterward, stability analysis is performed by using of the BVP4c solver function. According to Lund et al. [37,38], "this collocation formula and the collocation polynomial provides a C 1 continuous solution that is fourth order accurate uniformly in [a , b]. Mesh selection and error control are based on the residual of the continuous solution." Values of the smallest eigenvalue are given in Table 2. With the help of these values, it is expressed that the first solution (upper branch solution) is stable and physically realizable because the smallest eigenvalue is positive. Whereas the second solution (lower branch solution) is unstable because of the smallest eigenvalue is negative, so it is not physically realizable.  Numerical values of skin friction coefficient f (0), couple stress coefficient g (0), local Nusselt number −θ (0) and local Sherwood number −∅ (0) with variation of the stretching/shrinking parameter (λ) for different values of suction parameter f w are illustrated in Figures 2-5, where values of       In this work, the obtained critical points for variation of the λ are λc 0 = −1.4312, λc 1 − 1.6535 and λc 2 − 1.8933 for the different values of f w = 2.8, 3, 3.2, respectively. Thus, |λc| increases as the f w is increased. It is interesting to mention that dual solutions for this problem exist for the linearly stretching (λ > 0) as well as for the shrinking (λ < 0) surfaces. Figure 2 illustrates that in the second solution, the skin friction coefficient f (0) reduces as the value of f w is increased. Whereas in the first solution, the skin friction coefficient increases with an increment in the f w for λ ≤ 0 and for λ > 0 an opposite result is noticed. The effect of f w with variation of λ on the coefficient of the couple stress g (0) is presented in Figure 3. It is observed that the couple stress coefficient g (0) decreases in first solution for λ > 0 but at λ ≤ 0 opposite trend of the result is observed. In the second solution, it can be seen that the coefficient of the couple stress g (0) decreases in both cases of the stretching and shrinking surfaces throughout the flow. Moreover, Figure 4 shows the rate of heat transfer −θ (0) increases in both first and second solutions when λc ≤ λ for a different values of f w . The absolute value of critical point λc increases as the value of the suction parameter is increased. The rate of concentration (local Sherwood number) decreases in both first and second solutions with increment in the values of f w at variation of the λ when λ > −0.4 but for λ < −0.4 the rate of concentration is increased that is mentioned in Figure 5.   It is interesting to mention that dual solutions for this problem exist for the linearly stretching 0) as well as for the shrinking 0) surfaces. Figure 2 illustrates that in the second solution, the skin friction coefficient 0 reduces as the value of is increased. Whereas in the first solution, the skin friction coefficient increases with an increment in the for 0 and for 0 an opposite result is noticed. The effect of with variation of on the coefficient of the couple stress ′ 0 is presented in Figure 3. It is observed that the couple stress coefficient ′ 0 decreases in first solution for 0 but at 0 opposite trend of the result is observed. In the second solution, it can be seen that the coefficient of the couple stress ′ 0 decreases in both cases of the stretching and shrinking surfaces throughout the flow. Moreover, Figure 4 shows the rate of heat transfer ′ 0  Figure 6 represents the effect of material parameter (K) on skin friction coefficient f (0). It is analyzed when K = 0 (Newtonian case), dual similarity solutions exist in the range for f w ≥ 1.90906 and no solution at f w < 1.90906. It is worthwhile to specify the increasing values of K reduces the range of multiple solutions versus suction parameter f w . For K = 1 (K = 2), dual solutions occur when f w ≥ 2.35601 ( f w ≥ 2.7314) and no solution at f w < 2.35601 ( f w < 2.7314). Furthermore, any increment in the value of the suction parameter increases the rate of skin friction in first solution and decline in the second solution. The effect of K with variation of f w on couple stress coefficient, rate of heat and concentration transfers are shown in Figures 7-9, respectively. The couple stress is increased in first solution and diminished in second solution when K is increased along variation of f w . Heat transfer rate increases in both solutions by increasing the value of K along f w . Figure 9 shows that the rate of concentration decreases in first solution and against it, it increases in second solution as value of K is enhanced with variation of f w . Another point that is important is that in this figure, the increasing value of K decreases the range of variation of f w .
Energies 2019, 12, x; doi: FOR PEER REVIEW www.mdpi.com/journal/energies stress coefficient, rate of heat and concentration transfers are shown in Figures 7-9, respectively. The couple stress is increased in first solution and diminished in second solution when is increased along variation of . Heat transfer rate increases in both solutions by increasing the value of along . Figure 9 shows that the rate of concentration decreases in first solution and against it, it increases in second solution as value of is enhanced with variation of . Another point that is important is that in this figure, the increasing value of decreases the range of variation of .     Figure 10 represents the variations of velocity profile for different values of K. It is analyzed that the thickness of momentum boundary layer and velocity of micropolar nanofluid flow increase in the first solution as K is increased. On the other hand, velocity decreases in the second solution as the effect of micropolar material parameter (K) is enhanced, this is because of the increasing the rate of material increases the viscosity that makes resistance during fluid flow. The effect of the velocity slip parameter at the velocity profile is presented in Figure 11. The velocity boundary layer thickness and nanofluid velocity decrease as the velocity slip parameter δ is increased in the first solution. Whereas fluctuation is seen in second solution, the velocity of the nanofluid decreases at the start, but after a point it increases as δ is enhanced.     is enhanced, this is because of the increasing the rate of material increases the viscosity that makes resistance during fluid flow. The effect of the velocity slip parameter at the velocity profile is presented in Figure 11. The velocity boundary layer thickness and nanofluid velocity decrease as the velocity slip parameter is increased in the first solution. Whereas fluctuation is seen in second solution, the velocity of the nanofluid decreases at the start, but after a point it increases as is enhanced.  Figures 12 and 13 show the impact of m and K on microrotation profiles. In both, the microrotation boundary layer thicknesses and microrotation profiles are decreasing in the first and second solutions with increment in the value of the constants m and K. The behavior of temperature profile on account of thermal radiation is illustrated in Figure 14. The thermal radiation parameter Rd defines the relative contribution of conduction heat transfer to thermal radiation transfer. The expanding rate of radiation increase the thermal boundary layer and the rate of the heat transfer in both solutions. The impact of N t and N b on temperature profiles is shown in Figures 15 and 16, respectively. It is examined that any increment in thermophoresis parameter N t is caused to produce the high thermal conductivity of the nanoparticles in the base fluid, where the thickness of thermal boundary layer and temperature profile are enhanced in both solutions (see Figure 15). On the other hand, enhancement in Brownian motion parameter N b expands the temperature of the micropolar nanofluid (mentioned in Figure 16). This is due to "gradual development in nanoparticles percentage with Brownian motion parameter". Figure 17 demonstrates the effect of thermal slip parameter δ T on the temperature profile.
Both the thermal boundary layer thickness and temperature profile decline when thermal slip is enhanced in both solutions. Originally, enhancement in thermal slip parameter retards the fluid motion, which consequently shows a decline in net molecular movement. Consequently, less molecular momentum decreases the thermal boundary layer and the temperature profile subsequently diminishes.
Energies 2019, 12, x; doi: FOR PEER REVIEW www.mdpi.com/journal/energies Figure 9. The concentration transfer rate with different values with variation of . Figure 10 represents the variations of velocity profile for different values of . It is analyzed that the thickness of momentum boundary layer and velocity of micropolar nanofluid flow increase in the first solution as is increased. On the other hand, velocity decreases in the second solution as the effect of micropolar material parameter is enhanced, this is because of the increasing the rate of material increases the viscosity that makes resistance during fluid flow. The effect of the velocity slip parameter at the velocity profile is presented in Figure 11. The velocity boundary layer thickness and nanofluid velocity decrease as the velocity slip parameter is increased in the first solution. Whereas fluctuation is seen in second solution, the velocity of the nanofluid decreases at the start, but after a point it increases as is enhanced.    Figures 12 and 13 show the impact of and on microrotation profiles. In both, the microrotation boundary layer thicknesses and microrotation profiles are decreasing in the first and second solutions with increment in the value of the constants and . The behavior of temperature profile on account of thermal radiation is illustrated in Figure 14. The thermal radiation parameter          The effect of the thermophoresis parameter on the concentration profile is revealed in Figure 18. It can be clearly observed that any increment in increases the concentration profile and its boundary thickness in both solutions. Contrary to this, Figure 19 shows reverse behavior in both solutions when the Brownian motion parameter is increased. This is due to the fact that the increment in the magnitude of Brownian motion raises the rate of time at which nanoparticles  The effect of the thermophoresis parameter on the concentration profile is revealed in Figure 18. It can be clearly observed that any increment in increases the concentration profile and its boundary thickness in both solutions. Contrary to this, Figure 19 shows reverse behavior in both solutions when the Brownian motion parameter is increased. This is due to the fact that the increment in the magnitude of Brownian motion raises the rate of time at which nanoparticles The effect of the thermophoresis parameter on the concentration profile is revealed in Figure 18. It can be clearly observed that any increment in N t increases the concentration profile and its boundary thickness in both solutions. Contrary to this, Figure 19 shows reverse behavior in both solutions when the Brownian motion parameter N b is increased. This is due to the fact that the increment in the magnitude of Brownian motion raises the rate of time at which nanoparticles transport with distinct velocities in arbitrary direction. Figure 20 shows a small enhancement in concentration profile when δ C is enlarged in both solutions.
Energies 2019, 12, x; doi: FOR PEER REVIEW www.mdpi.com/journal/energies transport with distinct velocities in arbitrary direction. Figure 20 shows a small enhancement in concentration profile when is enlarged in both solutions.    Figure 20 shows a small enhancement in concentration profile when is enlarged in both solutions.

Conclusion
The steady two-dimensional laminar boundary layer flow of micropolar nanofluid over a linearly shrinking/stretching surface has been studied with the effects of thermal radiation and slip parameters. The governing equations of micropolar nanofluid flow are transformed to ordinary differential equations by applying similarity transformations. The equations are solved by the shooting method in the Maple software. The effects of the different physical parameters involved in the system of the equations were suction, micropolar material, radiation, Brownian motion, thermophoresis parameters, and the slip parameters. The main findings of this study are as follows:

Conclusions
The steady two-dimensional laminar boundary layer flow of micropolar nanofluid over a linearly shrinking/stretching surface has been studied with the effects of thermal radiation and slip parameters. The governing equations of micropolar nanofluid flow are transformed to ordinary differential equations by applying similarity transformations. The equations are solved by the shooting method in the Maple software. The effects of the different physical parameters involved in the system of the equations were suction, micropolar material, radiation, Brownian motion, thermophoresis parameters, and the slip parameters. The main findings of this study are as follows:

1.
Dual solutions exist for skin friction coefficient, couple stress coefficient, local Nusselt number and local Sherwood number for certain parameters.

2.
From stability analysis, it is examined that the first solution is stable and physically realizable, while the second solution is not stable so it is not physical realizable.

3.
There exist different ranges of dual similarity solutions. For K = 1 (K = 2), dual solutions occur when f w ≥ 2.35601 ( f w ≥ 2.7314) and no solution when f w < 2.35601 ( f w < 2.7314).

4.
The velocity profile declines in the first solution by increasing the values of the velocity slip parameter.

5.
The impact of m and K on microrotation profiles show that the microrotation boundary layer thicknesses and microrotation profiles decrease in both solutions. 6.
Temperature profiles increase in both solutions when thermal radiation, thermophoresis and the Brownian motion parameters are enhanced. 7.
The concentration of nanoparticles increases by increasing N t and decreases by increasing N b .