Effect of Viscous Dissipation in Heat Transfer of MHD Flow of Micropolar Fluid Partial Slip Conditions: Dual Solutions and Stability Analysis

Abstract

In this study, first-order slip effect with viscous dissipation and thermal radiation in micropolar fluid on a linear shrinking sheet is considered. Mathematical formulations of the governing equations of the problem have been derived by employing the fundamental laws of conservations which then converted into highly non-linear coupled partial differential equations (PDEs) of boundary layers. Linear transformations are employed to change PDEs into non-dimensional ordinary differential equations (ODEs). The solutions of the resultant ODEs have been obtained by using of numerical method which is presented in the form of shootlib package in MAPLE 2018. The results reveal that there is more than one solution depending upon the values of suction and material parameters. The ranges of dual solutions are $S\ge S_{ci}$, $i=0,1,2$ and no solution is $S<S_{ci}$ where $S_{ci}$ is the critical values of $S$. Critical values have been obtained in the presence of dual solutions and the stability analysis is carried out to identify more stable solutions. Variations of numerous parameters have been also examined by giving tables and graphs. The numerical values have been obtained for the skin friction and local Nusselt number and presented graphically. Further, it is observed that the temperature and thickness of the thermal boundary layer increase when thermal radiation parameter is increased in both solutions. In addition, it is also noticed that the fluid velocity increases in the case of strong magnetic field effect in the second solution.

1. Introduction

The study of the magnetohydrodynamic flow of non-Newtonian fluids has fascinated many scientists, researchers and mathematicians because there are numerous technical applications of non-Newtonian material such as greases, processed food, paints, ketchup, dyes, etc., [1]. The relation of shear stress and strain is non-linear the fluid is said to be a non-Newtonian fluid. Moreover, the dealing with the non-Newtonian fluid is more complicated as compared to Newtonian as the equation of Navier Stokes is related to many constraints, other physical parameters, and numbers [2]. It is not possible that one fluid model keeps all characteristics of the non-Newtonian fluid. Therefore, researchers introduced many non-Newtonian models such as model of Casson [3,4], Micropolar [5,6], Maxwell [7,8]; Jeffrey [9,10]; Burgers [11,12], etc. The micropolar model was introduced by Eringen in 1966 [13] in which he added the microrotation effect in the equation of Navier Stokes as the equation of Navier Stokes fails to interpret the behavior of such fluids. The micropolar fluid was examined by Bhattacharjee et al. [14] in which they also considered the porous effects and found only a single solution. Koriko et al. [15] investigated the micropolar fluid in the presence of thermal effects. Lakshmi et al. [16] studied the micropolar fluid and its heat characteristics by consideration of the second slip effect. Further, the following advanced studies have been carried out for micropolar fluids [17,18,19,20,21,22]. It is worth to note that the focus of all the above-mentioned studies is only in a single solution. Few related studies of micropolar fluid can be seen in these articles in which authors found multiple solutions successfully [23,24,25,26,27]. It has been concluded from the vast survey of published literature that the multiple solutions of micropolar fluid with the effect of viscous dissipation, thermal radiation and slip effect has not been studied so far.

Nowadays, nonlinear boundary values problems (BVPs) of fluid flows and their multiple solutions are important due to their wide applications in engineering and scientific research. In light of the nonlinear problems and their solutions, there may exist multiple solutions in numerous nonlinear fluid flow problems whose all possible solutions are hard to find [28]. Since all possible solutions satisfy the solution conditions of the considered fluid flow problems, it is essential not neglect any solution. Unfortunately, most of the current numerical techniques fail to recognize the multiple solutions because they are near to each other and may oscillate between two adjacent solutions [28]. These solutions exist in boundary layer fluid flow problems due to non-linearity in equations [29,30] as encountered in material science, electrochemistry, viscoelasticity, electro-magnetics, and acoustics. Soid et al. [31] considered steady MHD flow over a radially stretching or shrinking disk and noticed that dual solutions only exist on shrinking surface with small values of suction and the magnetic parameters. Khan et al. [32] investigated Carreau nanofluids and found the dual solutions by using the BVP4C method in MATLAB 2017 version. Raju et al. [33] discovered that dual solution exists only in a certain range of Power-law index values and found that reduction in the coefficient of skin friction and heat transfer rate due to increments in a field of a magnetic parameter. Lok et al. [34] examined the stagnation point flow of micropolar fluid and found dual solutions. Dual solutions of micropolar fluid in the presence of chemical reaction were investigated by Sandeep et al. [35]. Yacob et al. [36] found two solutions for the shrinking case when they dealt with micropolar fluid with the effect of stagnation. Turkyilmazoglu [37] extended the work of Bhattacharyya et al. [38] and found dual solutions of micropolar fluid in exact form. From an extensive review of published literature, it is noticed that micropolar fluid with the effect of partial slip conditions and viscous dissipation has not been investigated for the case of multiple solutions. Therefore, the main aim of this study is to find all possible solutions of MHD flow of the micropolar fluid with slip effect over a linear shrinking surface.

In the present decade, multiple solutions are considered in the fluid flow problem due to their wide range of applications in different areas of industries, engineering and so on. In order to determine which solution is physically realizable and stable, the stability analysis is employed by several researchers in their works. Based on the published works, Wilks and Bramley [39] were the first who found dual solutions in mixed convection flow and performed stability analysis of the solutions. After that, Merkin [40] also found the dual solutions in the fluid flow problem of mixed convection in a porous medium and carried out stability analysis to indicate that the first solution is more stable. The general definition of the theory of stability is regarded with measuring the behavior of perturbations of infinitesimal around the base state of finite-amplitude. If the perturbed flow returns to the base state, the flow is considered stable. On the other hand, if the small-amplitude perturbations diverge from the base state, the flow is considered unstable. Recently, Ali et al., [41] performed the stability analysis for MHD mixed convection flow of viscous fluid on the vertical surface and found that the first solution is a stable one. Lund et al. [42] examined the Williamson fluid, found dual solutions, and performed stability analysis. Khan et al. [43] investigated the non-Newtonian fluid and did stability analysis in order to indicate stable solutions. With existence of generation and absorption of heat, Jusoh et al. [44] found multiple solutions and examined the stability of solutions. Triple solutions have been found for the micropolar nanofluid over shrinking surfaces by Lund et al. [45]. The results of stability analysis exposed that only stable solution is the first one and the remaining two solutions are unstable. However, many researchers who found multiple solutions in fluid flow problems did not conduct this analysis [35,46,47,48]. It is worth to highlight that stability analysis should be considered in fluid flow problems when multiple solutions occur. Thus, this analysis is considered in this study to determine a more stable solution.

2. Problem Formulation

We considered MHD flow of micropolar fluid over the vertical shrinking surface at $y=0$. In addition, it is supposed that the velocity of a linear surface is $U_w=-ax$ and the temperature distribution of the sheet is $T_w\left(x\right)=T_{\infty }+T_0{x}^2$ where $T_{\infty }$ is the outside temperature of thermal boundary layer. A non-uniform magnetic field $B\equiv \left[0,B_0,0\right]$ is applied normal to the surface. We neglect the induced magnetic field by considering low magnetic Reynolds number. Moreover, the effect of the electric filed is not present (refer to Figure 1 for details). Along with all mentioned conditions the governing problems can be expressed in vector form as follows:

Now, Equations (1)–(4) can be written in the boundary layer form for steady flow as follows:

$$\mathit{\nabla }·\mathit{V}=\mathbf{0}$$

(1)

$$\rho \frac{d\mathit{V}}{dt}=-\mathit{\nabla }\mathit{P}+\left(\mu +\kappa \right){\mathit{\nabla }}^2\mathit{V}+\kappa \left(\mathit{\nabla }\times \mathit{N}\right)+\left(\mathit{J}\times \mathit{B}\right)$$

(2)

$$\rho j\frac{d\mathit{N}}{dt}=\gamma {\mathit{\nabla }}^2\mathit{N}-\kappa \left(2\mathit{N}-\mathit{\nabla }\times \mathit{V}\right)$$

(3)

$$\rho c_p\frac{d\mathit{T}}{dt}=k_1{\mathit{\nabla }}^2\mathit{T}-\mathit{\nabla }·{\mathit{q}}_{\mathit{r}}+\left(\mu +\kappa \right)\mathsf{\Phi }$$

(4)

where $t$ is the time, $\mathit{V}=\left[\overline{u}\left(x,y\right),\overline{v}\left(x,y\right)\right]$ is the velocity vector, $\rho =\rho \left(x,t\right)$ is the fluid density, $\mu$ is the dynamic viscosity of fluid, $\mathit{P}$ is the pressure, $=\frac{\partial }{\partial \overline{x}}\mathit{i}+\frac{\partial }{\partial \overline{y}}\mathit{j}$ is del or gradient operator, ${\mathit{\nabla }}^2=\frac{{\partial }^2}{\partial \overline{x^2}}+\frac{{\partial }^2}{\partial \overline{y^2}}$ is the Laplacian operator, $\mathit{T}=\overline{T}\left(x,y\right)$ is the temperature, $\mathit{J}\times \mathit{B}=-\mathit{\sigma }{\mathit{B}}_{\mathbf{0}}^{\mathbf{2}}\mathit{V}$ is the Lorentz force where $\mathit{B}$ denotes the total magnetic field which is the sum of applied magnetic field ${\mathit{B}}_{\mathbf{0}}$ and the induced magnetic field $\mathit{b}$, $\mathit{J}$ represents the current density, $\mathit{N}$ the micromotional or angular velocity, $\gamma =\left(\mu +\frac{\kappa }{2}\right)j$ is the spin gradient viscosity, $j$ is the gyration parameter of fluid, $\kappa$ is the vortex viscosity. Moreover, ${\mathit{q}}_{\mathit{r}}$ shows the Rosseland approximation which can be defined as ${\mathit{q}}_{\mathit{r}}=-\frac{4{\sigma }_1}{3k^{*}}\left[\frac{\partial T^4}{\partial x},\frac{\partial T^4}{\partial y}\right]$. Further, $\mathsf{\Phi }=2\left[{\left(\frac{\partial \overline{u}}{\partial \overline{x}}\right)}^2+{\left(\frac{\partial \overline{v}}{\partial \overline{y}}\right)}^2\right]+{\left(\frac{\partial \overline{v}}{\partial \overline{x}}+\frac{\partial \overline{u}}{\partial \overline{y}}\right)}^2$ is the dissipation function:

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

(5)

$$u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}=\left(\vartheta +\frac{\kappa }{\rho }\right)\frac{{\partial }^{2}u}{\partial y^2}+\frac{\kappa }{\rho }\frac{\partial N}{\partial y}-\frac{{\sigma }^{*}{B}^{2}u}{\rho }$$

(6)

$$u\frac{\partial N}{\partial x}+v\frac{\partial N}{\partial y}=\frac{1}{\rho j}\left[\gamma \frac{{\partial }^{2}N}{\partial y^2}-\kappa \left(2N+\frac{\partial u}{\partial y}\right)\right]$$

(7)

$$u\frac{\partial T}{\partial x}+v\frac{\partial T}{\partial y}=\left(\frac{k_1}{\rho c_p}+\frac{16{\sigma }_1{T}_{\infty }^3}{3k^{*}\rho c_p}\right)\frac{{\partial }^{2}T}{\partial y^2}+\left(\frac{\mu +\kappa }{\rho c_p}\right){\left(\frac{\partial u}{\partial y}\right)}^2$$

(8)

whose associated boundary conditions are:

$$\begin{gathered} v=v_0;u=U_w+D_1\frac{\partial u}{\partial y};N=-m\frac{\partial u}{\partial y};T=T_w+D_2\frac{\partial T}{\partial y}aty=0 \\ u\to 0;N\to 0;T\to T_{\infty };\mathrm{as}y\to \infty \end{gathered}$$

(9)

Here respective velocities of x-axis and y-axis are $u$ and $v$, $\vartheta$ is the kinematic viscosity, $k_1$ is the thermal conductivity, ${\sigma }^{*}$ is the electrical conductivity, ${\sigma }_1$ is the Stefan-Boltzmann constant, $k^{*}$ is the mean absorption coefficient, $m$ is the boundary parameter, $c_p$ is the specific heat at constant pressure, $v_0$ is the mass transfer velocity of the porous sheet, while $D_1$ and $D_2$ are the respective slip factors of velocity and temperature.

Following similarity variables are used to convert the PDEs into ODEs:

$$(=y\sqrt{\frac{a}{\vartheta }},u=ax{f}^{\prime }\left(\eta \right),v=-\sqrt{a\vartheta }f\left(\eta \right),N=\sqrt{\frac{a}{\vartheta }}axg\left(\eta \right),\theta \left(\eta \right)=\frac{\left(T-T_{\infty }\right)}{\left(T_w-T_{\infty }\right)}$$

(10)

By using of Equation (10) in Equations (5)–(9), we obtain:

$$\left(1+K\right)f^{‴}+f{f}^{″}-{f\prime }^2+K{g}^{\prime }-M{f}^{\prime }=0$$

(11)

$$\left(1+\frac{K}{2}\right)g^{″}+fg\prime -g{f}^{\prime }-2Kg-K{f}^{″}=0$$

(12)

$$\frac{1}{Pr}\left(1+\frac{4}{3}Rd\right){\theta }^{″}+f{\theta }^{\prime }-2f^{\prime }\theta +Ec\left(1+K\right)f^{″2}=0$$

(13)

with boundary conditions:

$$\begin{array}{c}f\left(0\right)=S,f^{\prime }\left(0\right)=-1+\alpha f^{″}\left(0\right),g\left(0\right)=-m{f}^{″}\left(0\right),\theta \left(0\right)=1+\beta \theta \prime \left(0\right)\\ f^{\prime }\left(\eta \right)\to 0;g\left(\eta \right)\to 0;\theta \left(\eta \right)\to 0;as\eta \to \infty \end{array}$$

(14)

where differentiation with respect to $\eta$ is denoted by prime, $K=\frac{\kappa }{\mu },M=\frac{{\sigma }^{*}{B}_0^2}{\rho a},Pr=\frac{\mu c_p}{k_1},Rd=\frac{4T_{\infty }^3{\sigma }_1}{k_1{k}^{*}},Ec=\frac{a^2}{T_0{c}_p},S=-\frac{v_w}{\sqrt{a\vartheta }},\alpha =D_1\sqrt{\frac{a}{\vartheta }},\beta =D_2\sqrt{\frac{a}{\vartheta }}$ are micropolar parameter, Hartmann number, Prandtl number, thermal radiation parameter, Eckert number, suction/blowing parameter, velocity slip parameter, and thermal slip parameter, respectively. It is worth noting that strong concentration can be observed when $m=0$ which means that microelements are unable to rotate. Further, weak concentration can be examined when $m=0$.5 which indicates that antisymmetric part of the stress tensor is vanishing. Similarly, $m=1$ shows the turbulent flow model of boundary layer.

The physical quantities of interest are skin friction coefficient $C_f$ with surface shear stress ${\tau }_w$ and local Nusselt number $N_u$ with surface heat flux $q_w$ are given by:

$$\begin{array}{c}{C}_f=\frac{2{\tau }_w}{\rho u_w^2},\mathrm{where}{\tau }_w={\left[\left(\mu +\kappa \right)\frac{\partial u}{\partial y}+\kappa N\right]}_{y=0}\\ N_u=\frac{x{q}_w}{k_1\left(T_w-T_{\infty }\right)},q_w=-k_1{\left[\left(1+\frac{16{\sigma }^{*}{T}_{\infty }^3}{3k_1{k}^{*}}\right)\frac{\partial T}{\partial y}\right]}_{y=0}\end{array}$$

(15)

Using Equation (10) in Equation (15), yields:

$$\begin{array}{c}{C}_f{\left(R{e}_x\right)}^{\frac{1}{2}}=\left(1+\left(1-m\right)K\right)f^{″}\left(0\right)\\ N_u{\left(R{e}_x\right)}^{-\frac{1}{2}}=-\left(1+\frac{4}{3}Rd\right){\theta }^{\prime }\left(0\right)\end{array}$$

(16)

where $R{e}_x=a{x}^2/\vartheta$ is local Reynolds number.

3. Stability Analysis

According to Weidman et al. [49] and Lund et al. [50], the existence of multiple solutions is possible in fluid flow problems in some ranges of physical parameters. It is necessary to carry out stability analysis in order to determine which one solution is stable. For this purpose, a new dimensionless time variable $\tau =at$ is presented to model the problem. Dimensionless similarity variables in terms of $\tau$ can be expressed as follows:

$$\begin{gathered} \eta =y\sqrt{\frac{a}{\vartheta }},\tau =at,u=ax{f}^{\prime }\left(\eta ,\tau \right),v=-\sqrt{a\vartheta }f\left(\eta ,\tau \right), \\ N=\sqrt{\frac{a}{\vartheta }}axg\left(\eta ,\tau \right),\theta \left(\eta ,\tau \right)=\frac{\left(T-T_{\infty }\right)}{\left(T_w-T_{\infty }\right)} \end{gathered}$$

(17)

As a result, Equations (6)–(8) are written as follows:

$$\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}=\left(\vartheta +\frac{\kappa }{\rho }\right)\frac{{\partial }^{2}u}{\partial y^2}+\frac{\kappa }{\rho }\frac{\partial N}{\partial y}-\frac{{\sigma }^{*}{B}^{2}u}{\rho }$$

(18)

$$\frac{\partial N}{\partial t}+u\frac{\partial N}{\partial x}+v\frac{\partial N}{\partial y}=\frac{1}{\rho j}\left[\gamma \frac{{\partial }^{2}N}{\partial y^2}-\kappa \left(2N+\frac{\partial u}{\partial y}\right)\right]$$

(19)

$$\frac{\partial T}{\partial t}+u\frac{\partial T}{\partial x}+v\frac{\partial T}{\partial y}=\left(\frac{k_1}{\rho c_p}+\frac{16{\sigma }^{*}{T}_{\infty }^3}{3k^{*}\rho c_p}\right)\frac{{\partial }^{2}T}{\partial y^2}+\left(\frac{\mu +\kappa }{\rho c_p}\right){\left(\frac{\partial u}{\partial y}\right)}^2$$

(20)

Substituting Equation (17) in Equations (18)–(20), we have:

$$\left(1+K\right)\frac{{\partial }^{3}f}{\partial {\eta }^3}+f\frac{{\partial }^{2}f}{\partial {\eta }^2}-{\left(\frac{\partial f}{\partial \eta }\right)}^2+K\frac{\partial g}{\partial \eta }-M\frac{\partial f}{\partial \eta }-\frac{{\partial }^{2}f}{\partial \tau \partial \eta }=0$$

(21)

$$\left(1+\frac{K}{2}\right)\frac{{\partial }^{2}g}{\partial {\eta }^2}+f\frac{\partial g}{\partial \eta }-g\frac{\partial f}{\partial \eta }-2Kg-K\frac{{\partial }^{2}f}{\partial {\eta }^2}-\frac{\partial g}{\partial \tau }=0$$

(22)

$$\frac{1}{Pr}\left(1+\frac{4}{3}Rd\right)\frac{{\partial }^2\theta }{\partial {\eta }^2}+f\frac{\partial \theta }{\partial \eta }-2\theta \frac{\partial f}{\partial \eta }+Ec\left(1+K\right)\frac{{\partial }^{2}f}{\partial {\eta }^2}-\frac{\partial \theta }{\partial \tau }=0$$

(23)

subject to boundary conditions:

$$\begin{array}{c}f\left(0,\tau \right)=S;\frac{\partial f\left(0,\tau \right)}{\partial \eta }=-1+\alpha \frac{{\partial }^{2}f\left(0,\tau \right)}{\partial {\eta }^2}\mathrm{and}\frac{\partial f\left(\eta ,\tau \right)}{\partial \eta }\to 0\mathrm{as}\eta \to \infty \\ h\left(0,\tau \right)=-m\frac{{\partial }^{2}f\left(0,\tau \right)}{\partial {\eta }^2}\mathrm{and}h\left(\eta ,\tau \right)\to 0\mathrm{as}\eta \to \infty \\ \theta \left(0,\tau \right)=1+\beta \frac{\partial \theta \left(0,\tau \right)}{\partial \eta }\mathrm{and}\theta \left(\eta ,\tau \right)\to 0\mathrm{as}\eta \to \infty \end{array}$$

(24)

The steady flow solutions ($f\left(\eta \right)=f_0\left(\eta \right)$, $g\left(\eta \right)=g_0\left(\eta \right)$, $\theta \left(\eta \right)={\theta }_0\left(\eta \right)$) must satisfy the corresponding boundary conditions in order to perform the solution’s stability. The complete functions of solution are expressed as:

$$\begin{array}{c}f\left(\eta ,\tau \right)=f_0\left(\eta \right)+e^{-\epsilon \tau }F\left(\eta ,\tau \right)\\ g\left(\eta ,\tau \right)=g_0\left(\eta \right)+e^{-\epsilon \tau }H\left(\eta ,\tau \right)\\ \theta \left(\eta ,\tau \right)={\theta }_0\left(\eta \right)+e^{-\epsilon \tau }G\left(\eta ,\tau \right)\end{array}$$

(25)

where $F\left(\eta ,\tau \right)$, $H\left(\eta ,\tau \right)$, and $G\left(\eta ,\tau \right)$ are small relative to $f_0\left(\eta \right)$, $g_0\left(\eta \right)$, and ${\theta }_0\left(\eta \right)$ respectively and $\epsilon$ is the unknow eigenvalue. By using Equation (23) in Equations (21)–(24) and by keeping $\tau =0$ (the steady state), we get the linearized form of the equations:

$$\left(1+K\right)F_0^{‴}+f_0{F}_0^{″}+F_0{f}_0^{″}-2f_0^{\prime }{F}_0^{\prime }+K{H}_0^{\prime }-M{F}_0^{\prime }+\epsilon F_0^{\prime }=0$$

(26)

$$\left(1+\frac{K}{2}\right)H_0^{″}+f_0{H}_0^{\prime }+F_0{g}_0^{\prime }-g_0{F}_0^{\prime }-g_0{F}_0^{\prime }-2K{H}_0-K{F}_0^{″}+\epsilon H_0=0$$

(27)

$$\frac{1}{Pr}\left(1+\frac{4}{3}Rd\right)G_0^{″}+f_0{G}_0^{\prime }+F_0{\theta }_0^{\prime }-2f_0^{\prime }{G}_0-2F_0^{\prime }{\theta }_0+2Ec.\left(1+K\right)f_0^{″}{F}_0^{″}+\epsilon G_0=0$$

(28)

with the boundary conditions:

$$\begin{array}{c}{F}_0\left(0\right)=0,F_0^{\prime }\left(0\right)=\alpha F_0^{″}\left(0\right),H_0\left(0\right)=-m{F}_0^{″}\left(0\right),G_0\left(0\right)=\beta G_0^{\prime }\left(0\right)\\ F_0^{\prime }\left(\eta \right)\to 0,G_0\left(\eta \right)\to 0,H_0\left(\eta \right)\to 0\mathrm{as}\eta \to \infty \end{array}$$

(29)

The stability of the steady solution has been studied by finding the values of smallest eigenvalue. According to Haris et al. [51] and Ali et al. [52], we need to relax the condition in order to determine the range of possible eigenvalues. For current problem, $F_0^{\prime }\left(\eta \right)\to 0$ as $\eta \to \infty$ is relaxed to $F_0^{″}\left(0\right)=1$.

4. Numerical Method

When the system of ODEs is non-linear then it is not possible to find the solution of the system in an exact form of solution. Henceforth, we look to the numerical approach. In this paper, the shooting method has been used to solve the system of non-linear ODEs with boundary conditions. The shooting method converts the boundary conditions into initial conditions and a higher order of ODEs to first-order ODEs. The resultant system of first-order ODEs has been solved by the Runge-Kutta 4th order method. The detail of the present method is as follows:

$$p=f\prime ,q=p\prime ,p_1=g\prime ;\left(1+K\right)q\prime =p^2-fq-K{p}_1+Mp$$

(30)

$$p_1=g\prime ;\left(1+\frac{K}{2}\right)p_1^{\prime }=gp-f{p}_1-2Kg-Kq$$

(31)

$$r=\theta \prime ;\left(1+\frac{4}{3}Rd\right)r\prime =Pr\left\{2p\theta -fr-Ec\left(1+K\right)q^2\right\}$$

(32)

with conditions:

$$\begin{array}{c}f\left(0\right)=S;p\left(0\right)=-1+\alpha q\left(0\right);q\left(0\right)={\alpha }_1\\ g\left(0\right)=-m{\alpha }_1;p_1\left(0\right)={\alpha }_2\\ \theta \left(0\right)=1+\beta r\left(0\right);r\left(0\right)={\alpha }_3\end{array}$$

(33)

where unknown initial conditions are denoted by ${\alpha }_1$, ${\alpha }_2$ and ${\alpha }_3$. It is worth mentioning that the missing initial values of ${\alpha }_1$, ${\alpha }_2$ and ${\alpha }_3$ are needed to shoot such that profiles of the solutions satisfy the original boundary conditions which are $f^{\prime }\left(\eta \right)\to 0;g\left(\eta \right)\to 0;\mathrm{and}\theta \left(\eta \right)\to 0\mathrm{as}\eta \to \infty$.

5. Results and Discussion

Steady MHD flow of electrically conducted micropolar fluid in the presence of viscous dissipation and thermal radiation have been examined. Moreover, the effect of velocity and thermal slip condition have also been taken into account. Shooting method with 4th order RK method has been employed with help of shootlib function in Maple 2018 to get the solutions of the governing non-linear ODEs as stated in Equations (11)–(13). The numerical results revealed that dual solutions exist for different ranges of applied parameters. In addition, the system of governing Equations (26)–(28) with boundary condition (29) has been solved by employing the BVP4C package presented in MATLAB in order to perform the stability analysis of solutions. BVP4C package is created on the finite difference scheme by applying three-stage Lobatto IIIa formula possessing fourth-order accuracy. A good initial guess is important to predict the existence of the dual solution. Normally, fix numbers are kept for the first solution (it is always not necessary), while the second solution depends upon the good initial guesses. These good guesses should satisfy the boundary conditions and fulfill boundary conditions asymptotically. Further, if any guess fails to generate the convergent solution then an appropriate set of values of parameters are needed to be selected using trial approach in order to make appropriate guesses for finding the second solution. The results of the smallest eigenvalue are given in

Table 1. The values of smallest eigenvalue for different values of S and K where M = 0.1, m = 0.5 and $\alpha =0.1$ are kept fixed.

K	S	$\mathit{\epsilon }$
		1st Solution	2nd Solution
0	1.85	0.43203	−0.02392
	1.95	0. 83673	−0.2750
	2.35	1.0874	−0.92698
1	2.35	0.3292	−0.2685
	2.5	0.7038	−0.4163

. It can be seen that the signs of the smallest eigenvalue for the first (second) solution are positive (negative) which mean that first (second) solution is stable (unstable). In order to validate our numerical results, we compared the values of skin friction coefficient with the exact solutions of Equation (11). As shown in

Table 2. The comparison of numerical solutions and exact solutions for Skin friction co-efficient.

K	S	${\mathit{f}}^{″}\left(0\right)$			${\mathit{f}}^{″}\left(0\right)$
		1st Solution			2nd Solution
		Exact	Numerical	Absolute Error	Exact	Numerical	Absolute Error
0	2	1.00	1.00009998	$9.0\times {10}^{-5}$	--	--	--
	2.5	2.00	1.99999999	$1.0\times {10}^{-8}$	0.50	0.49710246	$8.9754\times {10}^{-3}$
0.1	3	2.47185249	2.47185249	0	0.38529036	0.38486773	$4.2263\times {10}^{-4}$
	3.5	3.01773913	3.01773913	0	0.31559419	0.31386568	$1.72851\times {10}^{-3}$
0.5	3	2.00	2.00000000	0	0.40	0.39990306	$9.694\times {10}^{-5}$
	3.5	2.47703296	2.47703306	$1.0\times {10}^{-7}$	0.32296703	0.31959198	$3.3751\times {10}^{-3}$

, both results were found in excellent agreement with maximum absolute errors of the stable (unstable) solution is ${10}^{-5}$ (${10}^{-3}$). It worth to mention that if $\alpha =M=0$ and $m=0.5$, Equation (11) is reduced to $f\left(\eta \right)=S-\frac{1-e^{-\lambda \eta }}{\lambda }$ which implies that $f″\left(0\right)=\frac{S\pm \sqrt{S^2-2K-4}}{\left(2+K\right)}$ (refer to [37]). The dimensionless coupled linear boundary layer ODEs have been solved for velocity, micro-rotation, and temperature profiles with the impact of the emerging micropolar parameter, Hartmann number, Prandtl number, thermal radiation parameter, Eckert number, Suction/Blowing parameter, velocity slip parameter, and thermal slip parameter.

Figure 2 exhibits the velocity profile ($f\prime \left(\eta \right)$) with the effect of material parameter ($K$). The fluid velocity increases as values of $K$ get higher and hence the boundary layer thickness becomes thicker in the first solution. On the other hand, reverse behavior of velocity of fluid has been noticed in the second solution. This decreasing behavior is caused by the created resistance of micropolar parameter on the velocity field which reduces the fluid velocity. The velocity distribution, displayed in Figure 2, reveals that the motion rate is significantly increases for the no-slip condition. It is also observed that the thickness of momentum boundary layer reduces as the increase of velocity slip parameter as compared to $\alpha =0$ in both solutions. Figure 3 depicts the profile of velocity ($f\prime \left(\eta \right)$) for several values of $M$. The fluid velocity increases for the strong effect of the magnetic field in the second solution. However, the velocity and thickness of the momentum boundary layer tend to decrease in the first solution as expected. Physically, the magnetic force called a Lorentz force creates the intense drag force which resists the flow of fluid. Thus, the fluid velocity and hydrodynamic boundary layer thickness are reduced. The variation of material parameter $\left(K\right)$ on the micro-rotation profile ($g\left(\eta \right)$) is revealed in Figure 4. For the first solution when $m=0,1$ the micro-rotation profile rises in the range $0\le K\le 0.5$ and decreases in the range $0.5<K\le 1$. On the other hand, same behavior is noted for the second solution as noticed in the first solution.

The graph of the temperature profile ($\theta \left(\eta \right)$) in the presence and absence of Eckert number ($Ec$) for various values of Prandtl number ($Pr$) is drawn in Figure 5. Temperature of the fluid is noticed to reduce when $Pr$ increases in both solutions. Physically, the enhancement in the effect of the Prandtl number reduces the thermal diffusivity which causes the reduction in the thermal boundary layer thickness. Therefore, it can be concluded that the Prandtl number controls the thickness of thermal boundary layer and cooling in conducting flows. In the case of small Prandtl number, heat spreads quickly from the surface to the fluid flow which implies that the thickness of thermal boundary layer is much higher for the liquid metals, i.e., range of Prandtl number for liquid metals is 0.01 to 0.001. The interesting behavior is noticed for the second solution when $Ec=0$ where the thickness of thermal boundary layer increases when $1.5\le Pr\le 7$ and decreases when $7<Pr\le 10$. Figure 6 demonstrates the behavior of temperature profiles ($\theta \left(\eta \right)$) for increasing values of the thermal radiation ($Rd$). It is detected that the thickness of thermal boundary layer and temperature increase when thermal radiation parameter is increased in both solutions. This agrees with the physical fact that thermal radiation spreads heat from surface to fluid flow and consequently the thermal boundary layer becomes thicker. The situation can also be further explained by relationship between the Prandtl number and thermal radiation parameter given by $\frac{1}{Pr}\left(1+\frac{4}{3}Rd\right)$. It can be clearly seen that if the values of $Rd$ increase then the effect of thermal diffusivity increases. On the other hand, when the value of $Pr$ decreases, the thickness of thermal boundary and temperature of fluid increase.

The temperature distributions ($\theta \left(\eta \right)$) for various values of $K$ in presence and absence of thermal slip effect is shown in Figure 7. It is reckoned that the thickness of the thermal boundary layer tends to increase for the higher effect of the material parameter in the first solution. On the other hand, the reverse trend of temperature is found in the second solution. It is worth to notice that all plots (2–8) satisfy asymptotically the far-field boundary conditions for both solutions which show the validity of the shooting methods and generated numerical results. The variation of Eckert number on the non-dimensional fluid temperature profiles is revealed in Figure 8, indicating that both temperature and thickness of thermal boundary layer increase in the first solution. This increment of temperature is caused by smaller effect of enthalpy difference of boundary layer since Eckert number is the ratio of the kinetic energy flow and an enthalpy difference of boundary layer. Eckert number helps the converting of kinetic energy into internal energy by opposing fluid stresses. It is worth mentioning that the loss of temperature from the surface to the fluid flow is possible for $0<Ec\ll 1$. On the other hand, the opposite behavior of temperature profiles is examined.

It can be easily seen from Figure 9, Figure 10 and Figure 11 that multiple solutions exist for different ranges of $S$ which depend upon the values of the material parameter. From those figures, we observed that dual solutions exist for $S\ge S_{ci}$, $i=0,1,2$ and no solution exists when $S<S_{ci}$ where $S_{ci}$ is the critical values of $S$ when dual solutions occur. The variation of critical values totally depends upon the values of $K$. Variation of skin friction coefficient ($f″\left(0\right)$) with respect to $K$ for different values of $S$ is demonstrated in Figure 9. It can be examined easily that increments in suction produce more drag forces as compared to the material parameter ($K$) in the first solution. As a result, skin friction increases in the first solution. On contrary, the opposite trend is noticed for the second solution. The same behavior of couple stress coefficient ($g\prime \left(0\right)$) is examined in Figure 11 which we noticed for the skin friction coefficient in both solutions. The effect of suction and material parameter ($K$) on the heat transfer rate ($-\theta \prime \left(0\right)$) was drawn in Figure 10. There exists discontinuity in the second solution which indicate instability of the solution. Finally, it is also noticed that the heat transfer rate is lower for the higher values of material parameters in the first solutions. Moreover, the increasing of the suction would reason the rate of heat transfer to accelerate because the temperature is simply diffused over the permeable sheet relatively than an impenetrable sheet.

6. Conclusions

Studies on a micropolar fluid considering the effects of thermal radiation, velocity slip, thermal slip and viscous dissipation have been carried out to investigate the existence of all possible multiple solutions. Mathematical formulations of the micropolar problem has been gotten by employing the fundamental laws of conservations. The governing equations are transformed to equations of boundary layers. The highly non-linear coupled partial differential equations (PDEs) are then converted into non-dimensional ordinary differential equations (ODEs) by using of linear transformations. Some important outcomes of this study can be summarized as follows:

(1)

Dual solutions exist over the permeable shrinking sheet and the critical values $S_c$ of suction parameter depends on material parameter $K$.

(2)

Velocity field and the thickness of hydrodynamic boundary layer increase (decrease) in the first solution and decrease (increases) in the second solution when $K$ $\left(M\right)$ is increased.

(3)

Both solutions increase (decrease) for the strong effect of $Ec$ and $Rd$ ($Pr$).

(4)

The material parameter and mass suction parameter are inversely proportional to each other’s.

(5)

The results of numerical solutions are closed to the exact solutions.

(6)

The results of the smallest eigenvalues show that second solution is not stable as compared to the first solution.

(7)

The realizable flow situation is only possible for the first solution.