Dual similarity solutions of MHD stagnation point flow of Casson fluid with effect of thermal radiation and viscous dissipation: stability analysis

In this paper, the rate of heat transfer of the steady MHD stagnation point flow of Casson fluid on the shrinking/stretching surface has been investigated with the effect of thermal radiation and viscous dissipation. The governing partial differential equations are first transformed into the ordinary (similarity) differential equations. The obtained system of equations is converted from boundary value problems (BVPs) to initial value problems (IVPs) with the help of the shooting method which then solved by the RK method with help of maple software. Furthermore, the three-stage Labatto III-A method is applied to perform stability analysis with the help of a bvp4c solver in MATLAB. Current outcomes contradict numerically with published results and found inastounding agreements. The results reveal that there exist dual solutions in both shrinking and stretching surfaces. Furthermore, the temperature increases when thermal radiation, Eckert number, and magnetic number are increased. Signs of the smallest eigenvalue reveal that only the first solution is stable and can be realizable physically.

www.nature.com/scientificreports/ applications. Khan and Husain 13 explored Casson liquid on the circle disk where the system of governing ordinary differential equations was solved by utilizing homotopy method. Meanwhile, Hamid et al. 14 considered Casson liquid encased in the cavity. Further, double solutions were found by Hamid et al. 15 within the sight of thermal radiation. Stability analysis was additionally performed by utilizing of BVP4C in the MATLAB programming. Some ongoing improvements on Casson fluid can be found in these articles [16][17][18][19][20] .
The study of boundary layer Casson fluid flow on stretching and shrinking surfaces was widely investigated for single solution cases. Stretching surface has many applications in many manufacturing such as the extrusion of the molten polymers due to the slit die in productions of plastic sheets, paper productions, wires as well as in the fibers coating process of the food stuff. The qualities of final products in such processes depend heavily upon the cooling rate in the process of the heat exchange. Therefore, a MHD parameter is an important element to be considered so that the rate of the cooling can be controlled in order to obtain the desired quality products. Crane 21 proposed a method for solving incompressible steady-state 2-D viscous fluid on which later extended to many diverse aspects. Some of its recent important directions over-stretching flows can be seen in these articles [22][23][24][25][26] . Vajravelu 27 and Cortell 28 considered viscous fluid on a non-linear stretching sheet. In this paper, we also considered a stretching surface with the shrinking surface due to its extensive utilizations in different fields.
From the previous couple of years, multiple similarity solutions of fluid flow problems have been considered on the stretching and shrinking sheets in the presence and absence of bouncy effect by numerous scholars [29][30][31][32] . These multiple solutions exist due to several physical parameters impacts such as suction parameter, mixed convection parameter and so forth. Furthermore, the past researches demonstrated the probabilities of the occurrence of multiple similarity solutions of fluid flow over a shrinking sheet are more than over a stretching surface 31 . The possibilities of non-uniqueness solution of fluid flow on a stretching surface are probable when the flow is stagnation point flow or opposing flow. Similarly, multiple solutions for Newtonian fluids can be gotten easily as compared to non-Newtonian fluids. It is stated in the previously published literature that non-uniqueness of the solutions occurs because of the existence of non-linearity in the fluid equations 29,30 . However, the models of non-Newtonian fluids contain many non-linear terms and thus leads to non-existence of multiple solutions.
The boundary layer flows of fluid and their similarity solutions are gotten much attention due to their vast applications in many industrial fields 32 . In real situations, multiple solutions cannot be visualized in the boundary layer problems and difficult to be detected. Hence, many researchers fail to notice multiple solutions 33,34 . Multiple solutions of MHD fluid flow problems have been examined theoretically as well as numerically by numerous researchers. Dero et al. 35 obtained triple solutions during the investigation of micropolar fluid with thermal radiation effect. Further, Dero et al. 36 inspected the unsteady flow of nanofluid on the shrinking sheet and found double solutions in deaccelerated case. It seems that Ridha and Curie 37 are the pioneers who found dual solutions in the opposing flow. The motivation behind this study is to consider multiple similarity solutions of the MHD stagnation point flow of Casson fluid on a permeable exponentially stretching and shrinking surfaces unanimously with viscous dissipation and thermal radiation effect.

Mathematical formulation
There has been studied as steady incompressible 2-D stagnation point flow of Casson electrically leading fluid on an exponentially shrinking and stretching surfaces with the impact of thick viscous dissipations and thermal radiation. There has additionally been supposed that rheological equation of the state for the isotropic and the incompressible progression of the Casson fluid which are reported as (allude 10 ): where the plastic dynamic viscosity of the non-Newtonian fluid is meant by µ B , π signifies the result of deformation rate segment, that is, , π = e ij e ij is (i, j)th deformation rate part and π c is critical value of the π which depends on the non-Newtonian model, and P y indicates the yield stress of the fluid. Moreover, a system of cartesian coordinate is taken into account, where the x-axis is supposed alongside the shrinking/Stretching surface and the y-axis is normal to it. Further, u w = ae x/l is the shrinking and stretching velocity of surface. The uniform magnetic field is applied to the normal of the fluid flow B=B 0 e x/2l where B 0 is the constant magnetic strength (Fig. 1). The field of induced magnetic is disregarded as a result of the small estimation of the magnetic Reynolds number. With the above assumptions, we get following equations along the following boundary conditions S where S is the suction/injunction parameter, and u e = be x/l is the stagnation point.
Following similarity variables are used to get the similarity solutions and ψ is the stream function and can be written in velocity component as follows, By employing Eqs. (6) and (7) in Eqs.
(2)-(5), we obtained along boundary conditions is the Hartmann number, Pr = ϑ α denotes the Prandtl number, R d = (Tw−T∞)cp is Eckert number, is the stretching and shrinking parameter, and A = b a is the velocity ratio parameter of the stagnation point.
Physical quantities of interests are coefficient of skin friction, and the local Nusselt number, which are reported as:

Stability analysis
In this study, we found dual solutions of Eqs. (8)-(9) along with the boundary conditions (10) for both surfaces. Thus, it is needed to do the stability analysis to detect a stable solution which can be physically reliable after the time passes. According to Hamid et al. 15 , the upper branch solutions always show stability and physically reliability. On the other hand, the second solutions are unstable and consequently physically unreliable. The same remarks were reported by many researchers [refer to [38][39][40][41][42] ].
In order to accomplish the solutions' temporal stability, the unsteady form of Eqs. (3)-(4) must be considered by proposing the new dimensionless time variable τ where τ is related to solutions of Eqs. (8)- (9). Equations (3)-(4) can be expressed for the unsteady state flow as follows The new dimensionless time-dependent variable τ = a 2l e x/l .t is presented. Therefore, Eq. (6) can be written as follows: By substituting Eq. (15) in Eqs. (13)-(14), we get With corresponding boundary conditions According to Lund et al. 43 , Ismail et al. 44 , and Naganthran et al. 45 , the stability of the dual solutions is tested by perturbing the steady solution by using following functions where corresponding small relatives of f 0 (η) and θ 0 (η) are F 0 (η) and G 0 (η) and unknow eigenvalue is γ. It is worth to mention that perturbated function has been considered in the form of exponential as compared to power function as these functions increase and decrease more rapidly as compared to the power functions. By putting Eq. (19) into Eqs. (16)(17), we get linearized eigenvalue problem as follows: Subject to boundary conditions System of linearized eigenvalue problem Eqs. (20)-(21) along boundary conditions (22) is solved and obtained the infinite set of eigenvalues (γ 1 < γ 2 < γ 3 < · · · ).
The solution is said to be stable flow if and only if the sign of the γ 1 is positive which shows the initial decay, as time passes. On the other hand, if the sign of the γ 1 is negative, at that point the flow solution shows the initial growth of development and the solution is said to be an unstable solution, as time passes.

Result and discussion
In this segment of the article, we discuss about the effect of numerous physical rising parameters on temperature, velocity, rate of heat transfer, and coefficient of skin friction profiles. In order to validate the results of our numerical technique, the numerical results of − 1 + 1 β f ′′ (0) were compared with results obtained by Hussain et al. 46  Figure 2 exhibits the profile of velocity for numerous estimations of the stagnation point A in both shrinking and stretching surfaces. It is worth to note that A < 1(A > 1) demonstrates that surface velocity is more(less) than the free stream velocity, while A = 1 implies that the surface and free stream velocities are equivalent. The velocity profile decreases (rises) when the values of A are expanded in the first solution on the shrinking (stretching) surface. On the other hand, increasing and decreasing behaviors are seen in the second solution on both   www.nature.com/scientificreports/ surfaces. Velocity profile for various estimations of the Casson parameter β on the shrinking/stretching surface is illustrated in Fig. 3. The velocity boundary layer thickness declines at the higher values of β in the first solution over both surfaces due to the resistance created by β in the fluid flow. In addition, it is worth to mention that when β → ∞ , the fluid flow behaves like a simple viscous fluid. In other words, it becomes Newtonian fluid. In the second solution, the fluid velocity is expanded initially and gradually decreased in both surfaces. The impact of Hartmann number M on the profile of velocity is delineated in Fig. 4. Hydro boundary layer becomes thinner and the velocity of the fluid is also deaccelerated in both solutions for the maximum intensity of the magnetic parameter. Physically, this is caused by the expansion of Lorentz force which creates the resistance in the fluid flow inside the boundary. Therefore, the velocity of fluid declines. Figure 5 reveals the relationship between temperature profile and Casson parameter. It has been detected that singularity exists in the temperature distribution in the second solution for the case of stretching surface. Moreover, it is noticed that when Casson parameter is increased then thermal boundary layer and temperature  www.nature.com/scientificreports/ of fluid decrease in both solutions and surfaces. Physically, this reduction is caused by lower values of β which is associated to increment (declination) in the yield (shear) stress. The behavior of temperature profile for higher values of M is illustrated in Fig. 6. It is noticed that singularity exists in the range 0 < M ≤ 0.25 for the case of shrinking surface. This singularity assures us the instability of the second solution and explained in detail in the stability analysis section. It has been found that the fluid's temperature and thickness of thermal boundary layer rise when magnetic effect enhances for stable and unstable solutions on both surfaces. Physically, it happens due to the fact that higher values of magnetic parameter create a strong force of Lorentz which is an agent to produce more heat from the surface to the fluid. The effect of temperature profile for the thermal radiation was plotted in Fig. 7. The temperature of fluid enhances in both solutions and surfaces for the higher values of Rd . Figure 8 exposes the temperature profile for the various values of Ec . It is also found that the temperature and thickness of thermal boundary layer increase in both solutions and surfaces as well when viscous dissipation impact is  www.nature.com/scientificreports/ increased in the form of the Eckert number. It is worth mentioning that Ec << 1 indicates the balanced between convection and conduction in the energy equation. Practically, the increments in the temperature profiles can be explained as "the reduction of the boundary layer enthalpy difference for advanced values of Eckert number" 47 . Figure 9 displays the nature of temperature profile for the increasing values of the Prandtl number Pr. It is noticed that thinner thermal boundary is for the larger values of the Prandtl number in both solutions and surfaces. In addition, there exists singularity in the unstable solution for the stretching case. Moreover, this reduction in temperature and thickness of the thermal boundary layer is affected by the lower thermal diffusivity since the relationship between the Pr with thermal diffusivity is reciprocal to each other. Henceforth, the temperature of fluid decreases for the higher values of the Prandtl number. Figure 10 displays the correlation between the skin friction coefficients with λ and S . It is analyzed that the reduction in the suction lowers the skin friction in both solutions which infers that the contact of surface with molecules of fluid is decreasing for the lower effect of the suction. Moreover, an interesting behavior is observed  It is also noticed that dual solutions exist on both surfaces. Further, there are the ranges of the dual solutions ( ≥ ci ) and no-solution ( < ci ) which depend upon the values of ci where i = 1, 2, 3 . The same behavior is also noticed in Fig. 11. It is worth mentioning that stagnation point A has direct relationship with the skin friction coefficient. Figure 12 demonstrates the impact of higher values of β on f ′′(0) . It is examined that f ′′(0) is lower for the higher values of the non-Newtonian parameter in the first solution. This decreasing behavior of f ′′(0) is due to the inverse relation of shear stress and the yield stress in the fluid equations. On contrary, the reverse behavior is noted for the higher values of β in the second solution. Further, increments in the suction produce more (less) drag force in the first (second) solution. It is also found that no solution exists when S < S ci while dual solution obtained when for i = 1, 2, 3 . Figure 13 describes the behavior of heat transfer rate for different values of velocity ratio parameter . There are multiple singularities occur in the unstable solution which demonstrate the instability of the flow for the second solution. Moreover, the rate of heat transfer enhances for the high effect of the suction in the first solution. Finally, the temperature gradient was plotted in Fig. 14. The rate of heat transfer is advanced for the more noteworthy estimations of mass suction in the first solution. Then again, the turnaround pattern is perceived for the second solution.    www.nature.com/scientificreports/

Conclusion
The rate of heat transfer of the steady MHD stagnation point flow of Casson fluid on the shrinking/stretching surface in the existence of thermal radiation and viscous dissipation has been examined. By using similarity transformation, the self-similar nonlinear ODEs have been gotten and then solved by using the shooting method in MAPLE software. Double solutions occur for the different ranges of the velocity ratio parameter and mass suction parameter. Stability analysis is done for the solutions by using BVP4C solver in the MATLAB software and the results suggest that only the first solution is the stable. It is found that the velocity and its boundary layer thickness decrease for the greater A in the first solution on the shrinking surface. The thickness of the thermal boundary enhances for the advanced values of the Eckert number and thermal radiation parameter while opposite behavior of temperature profile is noticed for the higher values of the Prandtl number in both solutions. The coefficient of skin friction is reduced by the velocity ratio parameter for both solutions. The coefficient of skin friction decreases for the wall mass transfer parameter in the second solution. In the first solution, however, the behavior of the skin friction contradicts with the second solution.