Triple Local Similarity Solutions of Darcy-Forchheimer Magnetohydrodynamic (MHD) Flow of Micropolar Nanoﬂuid Over an Exponential Shrinking Surface: Stability Analysis

: In this paper, the MHD ﬂow of a micropolar nanoﬂuid on an exponential sheet in an Extended-Darcy-Forchheimer porous medium have been considered. Buongiorno’s model is considered in order to formulate a mathematical model with di ﬀ erent boundary conditions. The governing partial di ﬀ erential equations (PDEs) of the nanoﬂuid ﬂow are changed into a third order non-linear quasi-ordinary di ﬀ erential equation (ODE), using the pseudo-similarity variable. The resultant ODEs of the boundary value problems (BVPs) are renewed into initial value problems (IVPs) using a shooting method, and then the IVPs are solved by a fourth order Runge-Kutta (RK) method. The e ﬀ ects of various physical parameters on the proﬁles of velocity, temperature, microrotation velocity, concentration, skin friction, couple stress coe ﬃ cients, heat, and concentration transfer are demonstrated graphically. The results reveal that triple solutions appear when S ≥ 2.0337 for K = 0.1 and S ≥ 2.7148 for K = 0.2. A stability analysis has been performed to show the stability of the solutions; only the ﬁrst solution is stable and physically possible, whereas the remaining two solutions are not stable.


Introduction
Micropolar fluid is a polar fluid which contains rigid randomly oriented or spherical particles. It can be defined as a fluid with micro structures and belongs to the nonsymmetric stress tensor [1]. Furthermore, this fluid model is employed to analyze the behavior of liquid crystals and exotic polymeric fluid or lubricant colloidal suspensions. Ariman et al. [2,3], Eringen [4][5][6], and Lukaszewicz [7] discussed the properties and applications of the micropolar fluid in details. The concept of the electrically conducting fluids motion in the presence of a magnetic field is called magnetohydrodynamics, or MHD for short. The word MHD is the combination of the words magneto, hydro, and dynamics, which mean magnetic, fluid and motion, respectively. MHD is also known as magnetofluid dynamics and hydromagnetic, which can be defined as the study of the dynamics of the electromagnetic field and the electrically conducting fluids. Recently, Kumar et al. [8] examined the MHD flow of micropolar fluid with a porous medium. Micropolar fluid with an MHD effect on a shrinking sheet along a weak concentration has been considered by Gupta et al. [9]. Turkyilmazoglu [10] found the exact solution of micropolar fluid within the existence of the MHD effect. The MHD flow of micropolar fluids with a porous medium had been considered by many researchers, such as Sheikh et al. [11], Akhter et al. [12], Siddiq et al. [13], Dero et al. [14], Hayat et al. [15,16], Ahmed et al. [17], and Waqas et al. [18].
In the last couple of years, the use of nanofluid as a convectional fluid, in order to increase the heat transfer rate, has pulled in extensive consideration among researchers. Research demonstrated that dissolving different sorts of nanoparticles, such as nonmetal, polymeric and metal mixed in the base fluids, provides good thermal properties [19,20]. The term nanofluid, which was introduced by Choi and Eastman in 1995 [21], can be defined as a fluid that is a mixture of regular (base) fluids with nano-meter sized particles (less than 100 nm). These particles may contain oxides, carbon nanotubes, and metals. On the other hand, oil, ethylene glycol, and water are generally considered to be the base fluids. These fluids have different physical and chemical properties from regular fluids [22]. There are two approaches to study nanofluids, namely the experimental and numerical one. Many researchers considered the numerical approach to understand the behavior of nanofluids and introduced new concepts to understand them. Khanafer et al. [23] built up a model to contemplate the heat transfer improvement of Cu-water nanofluid in a two-dimensional enclosure. Meanwhile, Buongiorno [24] constructed a new non-homogeneous model in which velocity of base fluids and nanoparticles are not equal to zero. This model consists of seven slip parameters, which are Brownian diffusion, diffusiophoresis, gravity settling, fluid drainage, inertia, thermophoresis, and the Magnus effect. The references of the development of nanofluids can be found in the book by Nield and Bejan [25] and also in the published review articles on nanofluid, such as Mahian et al. [26][27][28] and Wong and Leon [29]. Recently, a few researchers have considered nanoparticles with non-Newtonian base fluid in the presence of MHD effects, such as Mahdy [30], Rehman et al. [31], Hamid et al. [32], Eid et al. [33], and Prabhakar et al. [34].
It can be observed from previously published literature that not much work has been done on the Extended-Darcy-Forchheimer porous medium, due to the fact that the governing equations cannot be reduced to self-similarity solutions through the use of a similarity transformation, particularly when using exponential similarity variables. Similarly, the MHD flow of micropolar nanofluid over an exponential shrinking surface has also not been considered because the equation of the angular velocity cannot be transformed into a self-similarity solution. Keeping in view these drawbacks, we attempt to employ a new approach which is a pseudo-similarity variable in the governing equations of fluid flow in order to obtain a local similar solution, as adopted by a few researchers in their studies [35][36][37][38].
The key objective of the present work is to consider the MHD flow of micropolar nanofluid over an exponential shrinking surface in an Extended-Darcy-Forchheimer porous medium. The resultant equations, after performing the pseudo-similarity variable in the form of a third-order non-linear quasi-ordinary differential equation, have been solved using the shooting method with the RK-method; we found triple solutions. When multiple solutions exist in any problem, it is necessary to conduct a stability analysis in order to determine the stable solutions. Consequently, this analysis is also taken into account in this research.

Problem Description and Formulation
The steady incompressible two-dimensional MHD flow of a micropolar nanofluid on an exponentially shrinking surface in an Extended-Darcy-Forchheimer porous medium is considered by + bu u in the Navier Stokes equation. The velocity of the shrinking surface in the form of exponential terms is given by U w (x) = U 0 e 2x , while the uniform magnetic field of the strength B 0 has been normally applied to it ( Figure 1). Due to a small value of the magnetic Reynolds number, the induced magnetic field is ignored. Under the consideration of the mentioned assumptions, the boundary layer equations of motion for the micropolar nanofluid, heat and concentration equations are expressed as: subject to the following boundary conditions: ∅ + ∅ + = 0 (15) subject to the boundary conditions below: Here, prime stands for the differentiation with respect to the new independent variable , = is the non-Newtonian parameter, is the permeability parameter, where = ⁄ is the local Reynolds number. We considered the following similarity transformations, as adopted by Sanjayanand and Khan [39], to solve Equations (1)-(5), with boundary condition (6): where u = ∂ψ ∂y and v = − ∂ψ ∂x are components of the velocity along the directions x and y respectively,ρ is the fluid density, ϑ is the kinematic viscosity, σ is the electrical conductivity of the fluid, B = B 0 e x 2l is the magnetic field with a constant magnetic strength B 0 , K 1 is the permeability of the porous medium, b is the local inertia coefficient, κ is the vortex viscosity, N is the microrotation, γ indicates the spin gradient viscosity, j is the ratio of the micro inertia and unit mass, T is the fluid temperature, and α is the thermal diffusivity of the micropolar nanofluid. Furthermore, τ 1 = (ρc) p (ρc) f is the ratio between the effective heat capacity of the nanoparticle material and the capacity of the fluid, D B is the Brownian diffusion coefficient, D T is the thermophoretic diffusion coefficient, T w is the temperature of the wall, T ∞ is the ambient temperature, C w is the concentration of the wall C ∞ is the ambient concentration, and 2l is the velocity slip factor. It might be mentioned that the range of m is 0 ≤ m ≤ 1; however, m is constant. In the case of m = 0, we have N = 0, which indicates that the strong concentration and micro-elements are near to the wall and are not rotatable. Furthermore, m = 0.5 shows a weak concentration, which causes the anti-symmetric part of the stress tensor to vanish. On the other hand, m = 1 indicates the turbulent boundary layer flows modeling (see [40,41]).
Using Equation (7) in Equations (2)-(5), we get the following partial differential equations: 1 Pr Furthermore, many authors considered γ = µ + κ 2 j = µ 1 + K 2 j, where κ = µK is the material parameter [41] in their work. In our problem, the terms of the Extended-Darcy-Forchheimer porous U 0 do not allow it to have self-similar solutions. For this reason, by using the pseudo-similarity variable, a local similarity solution can be obtained by equating the derivative of the functions of f , g, θ and ∅ with respect to x being equal to zero. This implies [39]. As a result, all the terms on the right-hand side become zero, and we get the following third-order non-linear quasi-ordinary differential equation: 1 Pr subject to the boundary conditions below: Here, prime stands for the differentiation with respect to the new independent variable η, K = κ µ is the non-Newtonian parameter, K 1 is the permeability parameter, 2l is the velocity slip, and S < 0 and S > 0 are the mass injunction and suction parameter, respectively.
The physical quantities of interest are the coefficient of the skin friction, the local Nusselt number and local Sherwood number, which are given by: where Re x = lu w /ϑ is the local Reynolds number.

Stability Analysis
Weidman et al. [42] initiated a study of the stability analysis of multiple solutions. Since then, some researchers, such as Rosca and Pop [43] and Lund et al. [44,45], performed stability analyses in their studies on multiple solutions of fluid flow problems. They found that only the first or upper solution has a physical meaning, while all of the remaining solutions (second or third) are not physically relevant or, in other words, are said to be unstable solutions. The first step in finding the stability of the solutions is to change the momentum, heat, and concentration equations into an unsteady form by considering a new variable τ. In our case, we have τ = U 0 2l e x l ·t, as defined in the paper of Rehman et al. [46]: The presence of τ is associated with initial value problems of the stable solution. Equating the derivative of functions with respect to x being equal to zero leads to the new similarity transfer variables in the presence of τ and η, which can be expressed as: By applying Equation (22) in Equations (18)-(21), we get: 1 Pr subject to the boundary conditions: In order to indicate the solution stability of f (η) = f 0 (η), g(η) = g 0 (η), θ(η) = θ 0 (η) and ∅(η) = ∅ 0 (η), which satisfy the equation of the boundary value problem (23)-(26) with boundary condition (27), we follow the suggestion of Rehman et al. [46] by introducing the following functions: where F(η, τ), G(η, τ), H(η, τ) and S (η, τ) are small relative to f 0 (η), g 0 (η), θ 0 (η) and ∅ 0 (η) of the steady state solutions. It should be noted that the range of these functions are 0 < F(η, τ) < 1, 0 < G(η, τ) < 1, 0 < H(η, τ) < 1 and 0 < S(η, τ) < 1. Furthermore, ε is an unknown eigenvalue parameter, which needs to be found. Substituting the values of the functions and their derivatives from Equation (28) in Equations (23)-(26) with the boundary condition (27), we have: 1 Pr subject to the boundary conditions: We assumed τ = 0 for Equations (23)- (26) in order to calculate the initial growth and decay of the solution of Equation (28), as recommended by Alarifi et al. [47]. Under these circumstances, F(η, τ), G(η, τ), H(η, τ) and S (η, τ) can be written as F 0 (η), G 0 (η), H 0 (η) and S 0 (η), respectively.
It is stated in the studies of Lund et al. [44,45] and Haris et al. [48] that eigenvalues can be determined if and only if the boundary condition of any one function of the following functions F 0 (η), G 0 (η), H 0 (η) and S 0 (η) can be relaxed into the initial condition by differentiating that function one more time. In this study, we relaxed F 0 (η) → 0 as η → ∞ and then solved the system of Equations (29)-(32) subject to Equation (33) along with the new relaxed boundary condition F 0 (0) = 1. It is worth mentioning that the sign of the smallest eigenvalues (ε) determines the stability of the solutions. The smallest eigenvalue is negative (positive), which indicates that the solution of the flow is unstable (stable) and that there is an initial growth (decay) of disturbances.

Results and Discussion
In order to fully understand the considered fluid flow model, the numerical study has been carried out for various important physical parameters, such as the magnetic parameter M, permeability parameter K 1 , Forchheimmer parameter F S , non-Newtonian parameter K, thermophoresis parameter N t , Brownian motion parameter N b , etc. The highly non-linear system of the quasi-ordinary differential Equations (12)- (15), along with the boundary conditions (16), have been solved by using the shooting method, and triple solutions were found. The value of η ∞ is chosen from 4 to 8, and it is worth noting that the value of η ∞ increases until the profiles of the velocity, temperature, and concentration converge asymptotically to the momentum, temperature, and concentration boundary layers, respectively. Figure 2 was drawn to analyze the effect of the material parameter K on the velocity profiles. The thickness of the velocity boundary layer increases when the micropolar parameter K is increased in the first and the second solutions, due to the fact that the material parameter reduced the drag force and that the hydrodynamic boundary layer therefore rose. We noted that when K = 0 (Newtonian fluid), the second solution did not exist. On the other hand, the dual nature of the velocity profile was observed in the third solution. Figure 3 depicts the variation of the velocity profiles for the different values of the Forchheimmer parameter F S . We observed that, due to increments in F S , the resistant force occurred when the fluid flow was flowing on the porous surface, and hence the velocity of the flow declined in the third solutions. The dual behavior of the velocity profile was noticed in the second solution. However, no change could be seen in the first solution when F S was increased. The velocity and thickness of the momentum boundary layers are inversely (directly) proportional to K 1 and M in the first (third) solution. On the other hand, the dual behavior of the velocity profile was seen in the second solution, as illustrated in Figures 4 and 5. The effect of the slip parameters λ and m on the velocity profiles are shown in Figures 6 and 7. In the first solution, the velocity boundary thickness was reduced as λ and m increased; this was due to the fact that the velocity of the fluid and surface have a big difference when the velocity slip factor is enhanced. For the second (third) solution, the velocity profiles decreased (increased) initially after they inclined (declined) when λ and m improved. The dual nature of the flow in some sense indicates that there is an initial growth of disturbance. Initially, the microrotation profile was reduced and then started to rise with increasing values of the micropolar parameter in the first solution, as shown in Figure 8. The thickness of the microrotation boundary layer was enhanced (reduced) as K was enhanced in the third (second) solution. The dual behavior of the microrotation profile was been observed in all solutions except the first solution when the value of m increased. The thickness of the microrotation boundary layer increased with increasing values of m in the first solution, as demonstrated in Figure 9. Figure 10 was drawn to analyze the variation of the temperature profiles for different values of the Prandtl number Pr. The thermal boundary layer thickness and temperature were incrementally reduced in the values of the Prandtl number Pr for all of the solutions, as expected. This is due to the fact that a high Prandtl number causes the thermal conductivity of nanofluid to diminish, and as a result the temperature is reduced. The variation of the temperature profiles for different values of the Brownian motion parameter N b is shown in Figure 11. It was observed in all solutions of the temperature profiles that the temperature and thermal layer thickness were enhanced with increasing values of N b . This is due to fact that the Brownian motion N b increases the kinetic energy of the nanofluid; thus, the temperature of the nanofluid increases. The temperature profile, with an effect resulting from the thermophoresis parameter N t , is illustrated in Figure 12. In all three solutions of the nanofluid flow problem, we noted that as the thermophoresis parameter N t rose, the temperature and thickness of the thermal layer increased. This is because the thermophoretic force is generated by N t and the temperature gradient, which pushes the flow of the nanofluid far from the boundary layer as a resulting thickness of the thermal boundary layer increases. The thickness of the concentration boundary layer declined with increasing values of the Brownian motion parameter N b in all solutions in Figure 13, which was expected. This was physically justified by the fact that Brownian motion is generated when nanoparticle and base fluid are mixed together in a nanofluid system. Since Brownian diffusion shows the conduction of heat under those circumstances, the thickness of the concentration boundary layer decreases. Figure 14 was sketched to examine the effect of N t on the concentration profile of nanoparticles. In all three solutions, the thickness of the concentration boundary layer was enhanced when thermophoresis increased.                                          Wang [49] and many researchers also stated in their studies that the similarity solution of the fluid flow problems over the shrinking surface is possible to obtain when sufficient wall mass suction is applied. The flow of Newtonian fluid is different from that of non-Newtonian fluid; it is observed for micropolar nanofluid that when the value of the micropolar parameter increases, a strong mass suction is required to obtain the solution. In this study, we discovered that there exist two regions for similarity solutions, namely multiple solutions and single solution, depending on the mass suction parameter. For = 0.1( = 0.2) there is a range of triple solutions when S ≥ 2.0337(S ≥ 2.7148), and a single similarity solution exists, S < 2.0337(S < 2.7148), as shown in Figure 15. When the suction is enhanced, the skin friction increases in the first solution and decreases in the second and the third solutions. Figures 16-18 were drawn to examine the effect of the suction S and micropolar parameter on the couple stress coefficient, and the heat and concentration transfer rate, respectively. In all graphs, when the suction is increased, the couple stress coefficient, heat, and concentration transfer rate are enhanced for all of the solutions. Wang [49] and many researchers also stated in their studies that the similarity solution of the fluid flow problems over the shrinking surface is possible to obtain when sufficient wall mass suction is applied. The flow of Newtonian fluid is different from that of non-Newtonian fluid; it is observed for micropolar nanofluid that when the value of the micropolar parameter K increases, a strong mass suction is required to obtain the solution. In this study, we discovered that there exist two regions for similarity solutions, namely multiple solutions and single solution, depending on the mass suction parameter. For K = 0.1(K = 0.2) there is a range of triple solutions when S ≥ 2.0337(S ≥ 2.7148), and a single similarity solution exists, S < 2.0337(S < 2.7148), as shown in Figure 15. When the suction is enhanced, the skin friction increases in the first solution and decreases in the second and the third solutions. Figures 16-18 were drawn to examine the effect of the suction S and micropolar parameter K on the couple stress coefficient, and the heat and concentration transfer rate, respectively. In all graphs, when the suction is increased, the couple stress coefficient, heat, and concentration transfer rate are enhanced for all of the solutions.  Wang [49] and many researchers also stated in their studies that the similarity solution of the fluid flow problems over the shrinking surface is possible to obtain when sufficient wall mass suction is applied. The flow of Newtonian fluid is different from that of non-Newtonian fluid; it is observed for micropolar nanofluid that when the value of the micropolar parameter increases, a strong mass suction is required to obtain the solution. In this study, we discovered that there exist two regions for similarity solutions, namely multiple solutions and single solution, depending on the mass suction parameter. For = 0.1( = 0.2) there is a range of triple solutions when S ≥ 2.0337(S ≥ 2.7148), and a single similarity solution exists, S < 2.0337(S < 2.7148), as shown in Figure 15. When the suction is enhanced, the skin friction increases in the first solution and decreases in the second and   By performing a stability analysis, the stability of the fluid flow solutions is achieved in this research. We need to perform a stability analysis when more than one solution exists in the flow problem. The main focus of this analysis is to determine which solution is stable and physically possible and which one is unstable. It is worth noting that the stability of the solution depends on the sign of the smallest eigenvalue. The value of the smallest eigenvalue is determined through Equation (24), for which we have to solve Equations (25)- (28), along with the boundary conditions (29). The smallest eigenvalues are demonstrated in Table 1 for different values of the suction and non-Newtonian parameters. It is clear from Table 1 that the negative (positive) values of the smallest   By performing a stability analysis, the stability of the fluid flow solutions is achieved in this research. We need to perform a stability analysis when more than one solution exists in the flow problem. The main focus of this analysis is to determine which solution is stable and physically possible and which one is unstable. It is worth noting that the stability of the solution depends on the sign of the smallest eigenvalue. The value of the smallest eigenvalue is determined through Equation (24), for which we have to solve Equations (25)- (28), along with the boundary conditions (29). The smallest eigenvalues are demonstrated in Table 1 for different values of the suction and non-Newtonian parameters. It is clear from Table 1 that the negative (positive) values of the smallest   By performing a stability analysis, the stability of the fluid flow solutions is achieved in this research. We need to perform a stability analysis when more than one solution exists in the flow problem. The main focus of this analysis is to determine which solution is stable and physically possible and which one is unstable. It is worth noting that the stability of the solution depends on the sign of the smallest eigenvalue. The value of the smallest eigenvalue is determined through Equation (24), for which we have to solve Equations (25)- (28), along with the boundary conditions (29). The smallest eigenvalues are demonstrated in Table 1 for different values of the suction and non-Newtonian parameters. It is clear from Table 1 that the negative (positive) values of the smallest By performing a stability analysis, the stability of the fluid flow solutions is achieved in this research. We need to perform a stability analysis when more than one solution exists in the flow problem. The main focus of this analysis is to determine which solution is stable and physically possible and which one is unstable. It is worth noting that the stability of the solution depends on the sign of the smallest eigenvalue. The value of the smallest eigenvalue is determined through Equation (24), for which we have to solve Equations (25)- (28), along with the boundary conditions (29). The smallest eigenvalues ε are demonstrated in Table 1 for different values of the suction and non-Newtonian parameters. It is clear from Table 1 that the negative (positive) values of the smallest eigenvalue ε indicate an initial growth (decay) of the disturbance that will interrupt (resume) the boundary layer separation and flow from becoming unstable (stable). It is worth mentioning that the stable solution always provides a good physical meaning which can be realized.

Conclusions
In this research, the MHD flow of micropolar nanofluid over an exponentially shrinking surface was considered with the effect of the porous and velocity slip. Exponential similarity variables were used to convert the partial differential equations into quasi-ordinary differential equations. The resultant equations were converted from BVPs to IVPs using a shooting method, after which the IVPs were solved by an RK-4th order method. After the findings of multiple solutions of nanofluid flow, a stability analysis was performed in order to indicate the stable solution by using the BVP4C solver in MATLAB software. The main summary findings of our research are as follow: • Triple solutions exist when S ≥ 2.0337 for K = 0.1 and when S ≥ 2.7148 for K = 0.2.

•
Dual solutions exist for the Newtonian case K = 0.

•
The study of critical points acknowledges the range of multiple solutions and single solutions.

•
The study of the stability analysis indicates that only the first solution is stable and that the remaining two solutions are unstable. Acknowledgments: Authors would like to thank YUTP 015LCO-078 for the financial support. The authors would also like to thank Universiti Utara Malaysia (UUM) for the moral and financial support in conducting this research.