Magnetohydrodynamic (MHD) Flow of Micropolar Fluid with Effects of Viscous Dissipation and Joule Heating Over an Exponential Shrinking Sheet: Triple Solutions and Stability Analysis

A numerical study was carried out to examine the magnetohydrodynamic (MHD) flow of micropolar fluid on a shrinking surface in the presence of both Joule heating and viscous dissipation effects. The governing system of non-linear ordinary differential equations (ODEs) was obtained from the system of partial differential equations (PDEs) by employing exponential transformations. The resultant equations were transformed into initial value problems (IVPs) by shooting technique and then solved by the Runge–Kutta (RK) method. The effects of different parameters on velocity, angular velocity, temperature profiles, skin friction coefficient, and Nusselt number were obtained and demonstrated graphically. We observed that multiple solutions occurred in certain assortments of the parameters for suction on a surface. The stability analysis of solutions was performed, and we noted that the first solution was stable while the remaining two solutions were not. The results also showed that the velocity of the fluid increased as the non-Newtonian parameter rose in all solutions. Furthermore, it was detected that the temperature of fluid rose at higher values of the Eckert number in all solutions.


Introduction
The boundary layer flow and heat transfer over a stretching/shrinking sheet have raised extensive interest worldwide, as seen throughout engineering literature, since it has many applications including glass-fiber production, plastic pieces aero-dynamic extrusion, hot rolling, paper production, etc. [1,2]. The micropolar fluid model has attracted more attention from researchers as compared to other non-Newtonian fluid models due to its ability to determine the impacts of local configuration and micro-motions of the fluid components that are disregarded by conventional models [3]. Several operative models have been carefully derived in the steady-state regime. Many researchers considered various non-Newtonian fluid models in their studies. Khan et al. [4] analyzed theoretically 1.
To extend the problem of Aurangzaib et al. [15] by considering the effect of magnetic, viscous, and Joule heating functions.

2.
To find the maximum number of the multiple solutions. 3.
To perform stability analysis of multiple solutions in order to determine a stable solution.
To the best of our knowledge, no such model has been considered for multiple solutions. Therefore, these results are new and would provide a good understanding of micropolar fluid.

Mathematical Formulation
A two-dimensional, laminar boundary layer electrically-conducting incompressible micropolar fluid flow with viscous dissipation along an exponentially shrinking surface is considered. The u, v are velocity components corresponding the directions x and y, respectively. The x-axis is with the surface of the sheet while the y-axis is perpendicular to it. Physical model of the problem is shown in Figure 1. According to the above-mentioned assumptions, along with boundary layer approximations, the concerned boundary layer equations of motion along with angular momentum and heat transfer can be written as [15,35]: Subject to these boundary conditions where η is the similarity variable, f (η) is non-dimensional form of the stream function, and θ(η) denotes the non-dimensional temperature whereas, prime stands for differentiation with respect to variable η.
The following exponential form of the similarity transformations are used [15], By using Equation (6) in order to transform Equations (1)-(4) in the system of non-dimensional ordinary differential equations, we have 1 Pr and boundary conditions will be as Here, prime indicates differentiation due to η, K = κ µ is the micropolar parameter (non-dimensional material parameter), we assume that γ is the given by ρU w is the magnetic field parameter, Pr = ϑ α indicates Prandtl number, Ec = U 2 w c p T 0 is the Eckert number, and S is the suction/injection parameter. The coefficient of skin friction, local couple stress, and Nusselt numbers is obtained in dimensionless form as where Re x is the Reynolds number.
Here, ′ prime indicates differentiation due to , = is the micropolar parameter (nondimensional material parameter), we assume that γ is the given by where = [15]. Further, = ( ) is the magnetic field parameter, = indicates Prandtl number, = is the Eckert number, and S is the suction/injection parameter. The coefficient of skin friction, local couple stress, and Nusselt numbers is obtained in dimensionless form as where is the Reynolds number.

Stability Analysis
In this paper triple solutions were found. When there exists more than one solution, it is necessary to conduct the analysis of stability of solutions in order to indicate the stable solution. According to Usama et al. [36] and Lund et al. [37], the first step of performing stability analysis is to change the governing equations in unsteady flow form by introducing a new similarity variable, = . .

The new governing equations become
where is the time. The similarity variables in Equation (6) with can be written as

Stability Analysis
In this paper triple solutions were found. When there exists more than one solution, it is necessary to conduct the analysis of stability of solutions in order to indicate the stable solution. According to Usama et al. [36] and Lund et al. [37], the first step of performing stability analysis is to change the governing equations in unsteady flow form by introducing a new similarity variable, τ = U w 2l e x l .t. The new governing equations become where t is the time. The similarity variables in Equation (6) with τ can be written as By substituting Equation (15) into Equations (12)-(14), we get In order to indicate the solution stability of f (η) = f 0 (η), g(η) = g 0 (η), and θ(η) = θ 0 (η) which satisfies the equation of boundary value problem in Equations (16)- (19). For this purpose, we follow the suggestion of Rehman et al. [38] by introducing the following functions: where ε is the unknown eigenvalue which is needed to determine, F(η, τ), G(η, τ), and H(η, τ) which are the small relatives of f 0 (η), g 0 (η), and θ 0 (η), respectively. Further, Equation (20) is submitted to Equations (16)- (19) by keeping τ = 0, we get 1 Subject to boundary conditions It is necessary to determine the initial growth of decay and the initial growth of disturbance of the values of the smallest eigenvalue. In order to find the initial growth of decay and disturbance of the eigenvalue, we have to relax one boundary condition as suggested by Haris et al. [39] and Lund et al. [40]. Physically, this relaxation in boundary condition provides the results of the smallest eigenvalue without loss of generality [41]. Mathematically, it is the step of the algorithm of the stability analysis. In this study, we relaxed the F 0 (η) → 0 as η → ∞ into F 0 (0) = 1.

Results and Discussion
In this study, the MHD flow of micropolar fluid on an exponentially shrinking surface was investigated. The governing dimensionless equations are solved numerically using the shooting method and multiple solutions are obtained. The comparison of skin friction with the existing work is presented in Figure 2 for the same values of K and m. Both solutions show the same critical values at parameter S. These dual solutions only exist for certain values of S c . At S = S c , both solutions meet each other. Beyond this critical value, no similarity solution exists as the boundary layer separates from the sheet. The variation of the shear stress with the parameter S is depicted in Figure 3 for various values of micropolar parameter K. As anticipated, in the first solution, the shear stress rises with all values of S and declines with K. This increment in the shear stress is due to the high values of suction since suction produces more resistance in the fluid flow. Yet, for the second and third solutions, the shear stress reduces with increasing S and has dual solutions at S c1 = 2.1255 and S c2 = 2.67528, respectively corresponding to K = 0.1 and K = 0.2. The variation of the couple stress coefficient with S is depicted in Figure 4 for different values of K. This variation is found to be negligible in the first solution, whereas the couple stress coefficient decreases (increases) with an increase in S (K) in the second solution. Physically, this is due to the fact that the material parameter supports the particles of the skew-symmetric of the fluid and therefore the couple stress coefficient increases in the second solution.
In the third solution, the variation of the couple stress coefficient is totally different. It increases with S and decreases with K. Like shear stress, the dimensionless heat transfer coefficient shows the dual nature at the same critical points ( Figure 5). All other parameters are kept constant. The first solution for the dimensionless heat transfer coefficient also shows the same nature as the first solution for shear stress. Moreover, heat transfer enhances in the first and second solutions for the high values of mass suctions.

Results and Discussion
In this study, the MHD flow of micropolar fluid on an exponentially shrinking surface was investigated. The governing dimensionless equations are solved numerically using the shooting method and multiple solutions are obtained. The comparison of skin friction with the existing work is presented in Figure 2 for the same values of K and m. Both solutions show the same critical values at parameter S. These dual solutions only exist for certain values of Sc. At S = Sc, both solutions meet each other. Beyond this critical value, no similarity solution exists as the boundary layer separates from the sheet. The variation of the shear stress with the parameter S is depicted in Figure 3 for various values of micropolar parameter K. As anticipated, in the first solution, the shear stress rises with all values of S and declines with K. This increment in the shear stress is due to the high values of suction since suction produces more resistance in the fluid flow. Yet, for the second and third solutions, the shear stress reduces with increasing S and has dual solutions at Sc1 = 2.1255 and Sc2 = 2.67528, respectively corresponding to K = 0.1 and K = 0.2. The variation of the couple stress coefficient with S is depicted in Figure 4 for different values of K. This variation is found to be negligible in the first solution, whereas the couple stress coefficient decreases (increases) with an increase in S (K) in the second solution. Physically, this is due to the fact that the material parameter supports the particles of the skew-symmetric of the fluid and therefore the couple stress coefficient increases in the second solution. In the third solution, the variation of the couple stress coefficient is totally different. It increases with S and decreases with K. Like shear stress, the dimensionless heat transfer coefficient shows the dual nature at the same critical points ( Figure 5). All other parameters are kept constant. The first solution for the dimensionless heat transfer coefficient also shows the same nature as the first solution for shear stress. Moreover, heat transfer enhances in the first and second solutions for the high values of mass suctions.             The velocity profiles of the micropolar fluid for the pertinent parameters are presented in Figures 6-8. Figure 6 illustrates the multiple solutions of the dimensionless velocity for several values of M. The remaining parameters are kept constant. It is important to note that these solutions satisfy both boundary conditions. When the intensity of the magnetic field is enlarged, the Lorentz force increases and thus, the velocity decreases in the first solution. Physically, this reduction in the thickness of the momentum boundary layer is due to the high resistant which is created by the Lorentz force. The dimensionless velocity increases when the effect of the magnetic field is enhanced in both the second and third solutions. The figure also displays that the hydrodynamic boundary layer thickness increases in the second and third solutions. Similarly, the effects of micropolar and slip parameters on the variation of the dimensionless velocity are depicted in Figures 7 and 8, respectively. In each case, three solutions are obtained. It is observed that the velocity of the fluid and its thickness of boundary increases for the higher values of micropolar parameter in all solutions (Figure 7). Physically, this enhancement in the thickness of the momentum boundary layer indicates that the micropolar parameter produces less drag force in the boundary layer separation, while the velocity profiles reduce as the slip parameter enhances (Figure 8), as expected. In each case, the boundary layer thickness is higher in the second and third solutions. This reduction is due to the existence of the singularities in the profiles of the velocity. The velocity profiles of the micropolar fluid for the pertinent parameters are presented in Figures  6-8. Figure 6 illustrates the multiple solutions of the dimensionless velocity for several values of M. The remaining parameters are kept constant. It is important to note that these solutions satisfy both boundary conditions. When the intensity of the magnetic field is enlarged, the Lorentz force increases and thus, the velocity decreases in the first solution. Physically, this reduction in the thickness of the momentum boundary layer is due to the high resistant which is created by the Lorentz force. The dimensionless velocity increases when the effect of the magnetic field is enhanced in both the second and third solutions. The figure also displays that the hydrodynamic boundary layer thickness increases in the second and third solutions. Similarly, the effects of micropolar and slip parameters on the variation of the dimensionless velocity are depicted in Figures 7 and 8, respectively. In each case, three solutions are obtained. It is observed that the velocity of the fluid and its thickness of boundary increases for the higher values of micropolar parameter in all solutions (Figure 7). Physically, this enhancement in the thickness of the momentum boundary layer indicates that the micropolar parameter produces less drag force in the boundary layer separation, while the velocity profiles reduce as the slip parameter enhances (Figure 8), as expected. In each case, the boundary layer thickness is higher in the second and third solutions. This reduction is due to the existence of the singularities in the profiles of the velocity.     The variation in the dimensionless angular velocity with the corresponding micropolar and slip parameters are depicted in Figures 9 and 10. In each case, three solutions are obtained, except = 0. In the first solution, the dimensionless angular velocity increases asymptotically and satisfies the farfield boundary conditions when and are increased. In the second solution, the angular velocity rises initially and then declines. The angular velocity decreases with an increase in the micropolar The variation in the dimensionless angular velocity with the corresponding micropolar and slip parameters are depicted in Figures 9 and 10. In each case, three solutions are obtained, except K = 0. In the first solution, the dimensionless angular velocity increases asymptotically and satisfies the far-field boundary conditions when K and m are increased. In the second solution, the angular velocity rises initially and then declines. The angular velocity decreases with an increase in the micropolar parameter. In the third solution, the angular velocity decreases to a minimum value and then increases to attain the far-field boundary conditions. The variation of the temperature profiles inside the thickness of the thermal layer for different values of magnetic field M, micropolar fluid parameter K, Eckert number Ec, and Prandtl numbers Pr are revealed in Figures 11-14, respectively. In each case, multiple solutions are investigated. Figure 11 reveals that the first solution is the real possible solution and there is small variation in the profile of velocity. In the second (third) solution, the dimensionless temperature decreases (increases) with the magnetic field. The increment in the temperature of the fluid is expected when the magnetic field increases but the interesting scenario happens in the second solution. The reduction in the profile of temperature in the second solution is due to the fact that the Lorentz force does not produce the resistance which is expected for the higher value of M. Similar behavior for the dimensionless temperature against different values of K can be observed in Figure 12. The Eckert number describes heat dissipation within the thermal boundary layer. It is also noticed that when K = 0 there are two solutions only. Further, as the Eckert number increases, the heat dissipation increases and consequently, the dimensionless temperature increases ( Figure 13). Physically, it can be explained that the thermal boundary layer thickness increases for the intense impact of kinetic energy which is created for higher values of the Eckert number, as a result, profiles of temperature and thickness of the thermal boundary layer enhance. Figure 14 explains the variation of the dimensionless temperature with the Prandtl number. As the Prandtl number increases, the thermal resistance increases and consequently, the dimensionless temperature decreases, and the momentum diffusivity dominates the behavior. parameter. In the third solution, the angular velocity decreases to a minimum value and then increases to attain the far-field boundary conditions. . The variation of the temperature profiles inside the thickness of the thermal layer for different values of magnetic field M, micropolar fluid parameter K, Eckert number Ec, and Prandtl numbers Pr are revealed in Figures 11-14, respectively. In each case, multiple solutions are investigated. Figure   Figure 9. Angular velocity profile for different values of K.
Symmetry 2020, 12, x FOR PEER REVIEW 10 of 16 parameter. In the third solution, the angular velocity decreases to a minimum value and then increases to attain the far-field boundary conditions. . The variation of the temperature profiles inside the thickness of the thermal layer for different values of magnetic field M, micropolar fluid parameter K, Eckert number Ec, and Prandtl numbers Pr are revealed in Figures 11-14, respectively. In each case, multiple solutions are investigated. Figure   Figure 10. Angular velocity profile for different values of m. created for higher values of the Eckert number, as a result, profiles of temperature and thickness of the thermal boundary layer enhance. Figure 14 explains the variation of the dimensionless temperature with the Prandtl number. As the Prandtl number increases, the thermal resistance increases and consequently, the dimensionless temperature decreases, and the momentum diffusivity dominates the behavior.      Finally, Table 1 shows the values of the smallest eigenvalue. It is worth highlighting that the stability of solutions depends on the sign of the values of smallest eigenvalue . If the value of is positive, then the initial growth of decay exists, and the flow becomes stable. Hence, the first solution Finally, Table 1 shows the values of the smallest eigenvalue. It is worth highlighting that the stability of solutions depends on the sign of the values of smallest eigenvalue ε. If the value of ε is positive, then the initial growth of decay exists, and the flow becomes stable. Hence, the first solution is stable because the sign of ε is positive. On the other hand, if the value of ε is negative, then an initial growth of disturbance exists, and therefore the second and third solutions are said to be unstable.

Conclusions
MHD flow of micropolar fluid over an exponential shrinking surface was addressed. Hartmann number, viscous dissipation, and Joule heating have been accounted for. Multiple solutions and stability analysis were the main objectives of this study. The shooting method was implemented to get the solutions and the three-stage Lobatto IIIA formula, developed in bvp4c with the help of finite difference code, was used to determine the stability of the solutions. The main points of this study are as follows:

1.
Two solutions exist for the case of Newtonian fluid.

2.
Three solutions exist for the case of no-Newtonian fluid in the specific values of the suction parameter.

3.
Ranges of single and multiple solutions are dependent on the suction parameter.

4.
Results of the stability analysis of solutions indicate that only the first solution is stable.

5.
The thickness of the momentum boundary layer enhances in all solutions when the slip parameter is increased. 6.
The thickness of the thermal boundary layer is directly proportional to the values of the Eckert number. As the Eckert number increases, the temperature of fluid also rises due to the high impact of the kinetic energy. 7.
The occurrence of a higher velocity of the fluid is possible for unstable solutions when the magnetic parameter increases. 8.
The thermal field has been noted as lower corresponding to a larger Prandtl number in all solutions. 9.
The angular velocity of the micropolar fluid increased in the first and third solutions for the higher values of material and slip parameters.