Convective E ﬀ ect on Magnetohydrodynamic (MHD) Stagnation Point Flow of Casson Fluid over a Vertical Exponentially Stretching / Shrinking Surface: Triple Solutions

: In the current study, the characteristics of heat transfer of a steady, two-dimensional, stagnation point, and magnetohydrodynamic (MHD) ﬂow of shear thickening Casson ﬂuid on an exponentially vertical shrinking / stretching surface are examined in attendance of convective boundary conditions. The impact of the suction parameter is also considered. The system of governing partial di ﬀ erential equations (PDEs) and boundary conditions is converted into ordinary di ﬀ erential equations (ODEs) with the suitable exponential similarity variables of transformations and then solved using the shooting method with the fourth order Runge–Kutta method. Similarity transformation is an important class of phenomena in which scale symmetry allows one to reduce the number of independent variables of the problem. It should be noted that solutions of the ODEs show the symmetrical behavior of the PDES for the proﬁles of velocity and temperature. Similarity solutions are obtained for the case of stretching / shrinking and suction parameters. It is revealed that there exist two ranges of the solutions in the speciﬁc ranges of the physical parameters, three solutions depend on the opposing ﬂow case where stagnation point ( A ) should be equal to 0.1, two solutions exist when λ 1 = 0 where λ 1 is a mixed convection parameter and A > 0.1 , and a single solution exists when λ 1 > 0 . Moreover, the e ﬀ ects of numerous applied parameters on velocity, temperature distributions, skin friction, and local Nusselt number are examined and given through tables and graphs for both shrinking and stretching surfaces.


Introduction
The topic of boundary layer fluid flow on the stretching surfaces is widely investigated for single solution cases. Stretching surface has many applications in many processes of manufacturing, such as the extrusion of molten polymers due to the slit die in the production of plastic sheets, paper production, wires as well as in the fiber coating process of foodstuffs. The qualities of final products in such processes largely depend upon the cooling rate in the process of heat exchange. A magnetohydrodynamic parameter is an important parameter to use, such that the rate of the cooling might be controlled so that the desired products with quality can be obtained. Crane [1] gave the method of solution for the incompressible steady-state 2D boundary layer flow in the case of viscous fluid on the stretching surface. Furthermore, the problem examined by Crane has been studied more by extending it into many diverse aspects-some of its recent important directions over the stretching flows were considered by Mabood et al., [2], Rana et al., [3], Hamad [4], Hassan et al. [5] and Haq et al., [6]. Hamad and Ferdows [7] examined metallic nanofluid over a stretching surface. Vajravelu [8] and Cortell [9] considered viscous fluid on the non-linear stretching sheet. In this paper, we also considered the stretching surface because of its extensive usages in different areas.
During the past few decades, the stagnation-point flow has engrossed the intention of many scholars because of its growing industrial and scientific utilization-for example, electronic devices cooled by fans, cooling of nuclear reactors during an emergency shutdown, solar central receivers exposed to wind currents, hydrodynamic processes, etc. In fluid mechanics problems, the point in the flow field where the fluid local velocity becomes zero is known as a stagnation point. This point occurs at an object facade where the fluid is at rest because of a force exerted by the object. Hayat et al. [10] considered the stagnation point of Casson fluid flow on the vertical linear surface by including the viscous dissipation effect. The governing equations were resolved by homotopic analysis and the velocity field and fluid temperature were calculated to decrease Prandtl count functions. Casson fluid can be defined as a shear thinning fluid which is presumed to include an unlimited viscosity at zero shear rate, a yield stress at which no flow exists, and zero viscosity at an infinite shear rate. Abbas et al. [11] inspected Casson fluid on the stretching/shrinking surface where a fluid is the stagnation point flow. It was concluded during the examination that the velocity field rises by increasing the Casson parameter, but the thickness of the momentum boundary layer declines. In an analysis of stagnation point flow, Bhattacharyya et al. [12] obtained a dual solution. Recently, Nadeem et al. [13] investigated the three-dimensional stagnation point flow of nanofluids and obtained dual solutions. They concluded that a stable solution is a primary solution. The significant research is listed on such a phenomenon can be seen in the papers of Bhattacharyya [14]; Hayat et al. [15]; Ramesh et al. [16]; Haldar et al. [17]; Shafiq et al. [18]; and Jafarimoghaddam [19].
In the last couple of years, the non-uniqueness of solutions of fluid flow problems has been considered on the moving, shrinking and stretching sheets in the presence and absence of a bouncy effect by many researchers [20,21]. These multiple solutions exist due to several physical parameters' impacts, such as the suction parameter, mixed convection parameter, and so forth. Most past investigations demonstrated that the probabilities of occurrence of numerous solutions of fluid flow on a shrinking surface are higher, as a contrast with stretching surfaces [12]. The possibilities of non-uniqueness solutions of the flow field on a stretching surface are probable when a flow is opposing or when there is stagnation point flow. Similarly, multiple solutions for Newtonian fluids can be obtained easily as compared to non-Newtonian fluids. Several researchers have proposed that because of non-linearity in the fluid model, numerous solutions exist [20,22]. However, the models of non-Newtonian fluids contain many non-linear terms and because of these terms, multiple solutions are not going to exist. The value of several solutions for fluid boundary layers has significantly increased due to the occurrence of these fluids in many industrial applications [23]. Numerous solutions in actual situations of fluid flow problems cannot be visualized, because in reality many researchers have failed to notice multiple solutions [24]. The multiple solutions of magnetohydrodynamic (MHD) fluid flow problems have been examined experimentally as well as numerically by several scholars. Ishak et al. [25] considered the mixed convective MHD flow next to incompressible viscous fluid flow and obtained dual solutions. Subhashini et al. [26] considered the MHD convection mixed flow over a vertical surface and found the dual solution using the shooting technique. Ridha and Curie [27] were possibly the first to indicate that there are dual solutions in the opposing flow system and that the flow regime persists. The goal of this research is to examine the stagnation point flow of Casson fluid on a permeable exponentially vertical stretching and shrinking surfaces unanimously with the effect of convective boundary condition numerically.
The governing PDEs are first transformed into ODEs. The obtained system of equations is converted from boundary value problems (BVPs) to initial value problems (IVPs) with the help of the shooting method. Finally, the system of equations of IVPs is solved by the Runge Kutta (RK) method with the aid of maple software in the shootlib function. Furthermore, the three-stage Labatto three-A method is applied to performed stability analysis with the help of a bvp4c solver in MATLAB.

Mathematical Description of the Problem
We considered an incompressible, steady, 2-D stagnation point flow of Casson electrically conducting fluid on an exponentially vertical shrinking/stretching sheet with the convective boundary condition effect, as revealed in Figure 1. It is also supposed that the state of the isotropic rheological equation and the incompressible flow of the Casson fluid is defined as (see Lund et al. [28]): where µ B is the plastic dynamic viscosity, P y denotes the yield stress of the fluid, and the product of deformation rate component is denoted by π where π = e ij e ij is the (i, j)-th deformation rate component and the critical value of π is π c . The uniform magnetic field of the strength is B, which was applied perpendicularly to the stretching/shrinking surface (refer to Figure 1). It is supposed that a constant convective sheet temperature is T f , while T ∞ is the ambient fluid temperature where T f > T ∞ . Due to the very low value of the magnetic Reynolds number, the induced magnetic field vanished. The governing PDEs are first transformed into ODEs. The obtained system of equations is converted from boundary value problems (BVPs) to initial value problems (IVPs) with the help of the shooting method. Finally, the system of equations of IVPs is solved by the Runge Kutta (RK )method with the aid of maple software in the shootlib function. Furthermore, the three-stage Labatto three-A method is applied to performed stability analysis with the help of a bvp4c solver in MATLAB.

Mathematical Description of the Problem
We considered an incompressible, steady, 2-D stagnation point flow of Casson electrically conducting fluid on an exponentially vertical shrinking/stretching sheet with the convective boundary condition effect, as revealed in Figure 1. It is also supposed that the state of the isotropic rheological equation and the incompressible flow of the Casson fluid is defined as (see Lund et al. [28]): where is the plastic dynamic viscosity, denotes the yield stress of the fluid, and the product of deformation rate component is denoted by where π = is the (i, j)-th deformation rate component and the critical value of is . The uniform magnetic field of the strength is B, which was applied perpendicularly to the stretching/shrinking surface (refer to Figure 1). It is supposed that a constant convective sheet temperature is , while is the ambient fluid temperature where . Due to the very low value of the magnetic Reynolds number, the induced magnetic field vanished.   The equations of steady Casson fluid flow over the vertical surface with all assumptions can be expressed as follows: Along with boundary conditions In these equations, u and v are the velocity components along the direction of xand y-axis, respectively. Moreover, ρ is the fluid density, β is the Casson fluid parameter, ϑ is the kinematic viscosity, σ is the electrical conductivity of fluid, B = B 0 e x 2l is the magnetic field where B 0 is the magnetic constant, T is the fluid temperature, and α is the thermal diffusivity. Furthermore, T w = T ∞ + T 0 e 2x is the temperature of wall, and T ∞ is the ambient temperature. In addition, the thermal conductivity of Casson fluid is k f and the convective heat transfer coefficient is S where S is the suction/injunction parameter, and u w = a e x l is the velocity of surface. We look for the similarity solutions of Equations (2)-(4) by introducing the following transformations and the stream function ψ can be expressed in velocity components as, By applying Equations (6) and (7) in Equations (2)-(5), we have 1 Pr Along with these conditions where M = 2lσ(B 0 ) 2 ρU w is the magnetic parameter, λ is the stretching/shrinking parameter, Re 2 is the mixed convection parameter where Gr is the heat Grashof number. It should be noted that the flow in opposite (assisting) flow is indicated by λ 1 > 0 (λ 1 < 0). The physical quantities of interest are the skin friction coefficient and the local Nusselt number, which are written in the following terms: where Re x is the local Reynolds number.

Stability Analysis
In the current study, we found triple solutions of Equations (8) and (9) along with boundary conditions (10) for both surfaces. Thus, it is needed to perform stability analysis to notice a stable solution which can be physically reliable after the time passes. In order to accomplish the solutions' temporal stability, the unsteady form of Equations (3) and (4) must be considered by proposing the new dimensionless time variable τ where τ is related to the initial solutions of Equations (8) and (9) (without perturbation functions). Equations (3) and (4) can be expressed for the unsteady state flow as follows The new dimensionless time-dependent variable is τ = a 2l e x l .t where t is time. Therefore, Equation (6) can be written as follows: By putting Equation (15) in Equations (13) and (14), it is obtained 1 Pr With corresponding boundary conditions According to Lund et al. [29,30], Hamid et al. [31], and Mustafa et al. [32], the stability analysis of the dual solutions is verified 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 the unknown eigenvalue is γ. It is worth mentioning that perturbated function is considered in the form of an exponential as compared to power function, as these functions increase and decrease more rapidly as compared to the power functions. By putting Equation (19) into Equations (16) and (17), we get the linearized eigenvalue problem as follows: 1 with boundary conditions System of linearized eigenvalue problem Equations (20) and (21) along with boundary conditions (22) is solved and obtained the infinite set of eigenvalues (γ 1 < γ 2 < γ 3 < . . .). The solution is said to be a 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.

Three-Stage Lobatto IIIA Formula
The three-stage Lobatto IIIa formula is a well-known numerical method. This method can easily solve all kinds of non-linear and linear ODEs. It is developed with the support of finite difference code in bvp4c. After that, the study of stability is conducted by employing of the bvp4c solver function. According to Rehman et al. [33], "this collocation formula and the collocation polynomial provides a C 1 continuous solution that is fourth-order accurate uniformly in [a , b]. Mesh selection and error control are based on the residual of the continuous solution". Furthermore, the tolerance of relative error is fixed at 10 −5 . The suitable mesh determination is created and returned in the field sol.x. The bvp4c returns solution, named sol.y., as a construction. In any case, values of the solution are gotten from the array named sol.y, relating to the field sol.x. The general procedure of this technique, along with the stability analysis, is shown in Figure 2.
Symmetry 2020, 12, x FOR PEER REVIEW 6 of 16 with boundary conditions System of linearized eigenvalue problem Equations (20) and (21)

Three-Stage Lobatto IIIA Formula
The three-stage Lobatto IIIa formula is a well-known numerical method. This method can easily solve all kinds of non-linear and linear ODEs. It is developed with the support of finite difference code in bvp4c. After that, the study of stability is conducted by employing of the bvp4c solver function. According to Rehman et al. [33], "this collocation formula and the collocation polynomial provides a continuous solution that is fourth-order accurate uniformly in , . Mesh selection and error control are based on the residual of the continuous solution". Furthermore, the tolerance of relative error is fixed at 10 . The suitable mesh determination is created and returned in the field sol.x. The bvp4c returns solution, named sol.y., as a construction. In any case, values of the solution are gotten from the array named sol.y, relating to the field sol.x. The general procedure of this technique, along with the stability analysis, is shown in Figure 2.

Discussion
In the following section, we discuss the numerical results of Equations (8) and (9) with their boundary conditions (10). These equations are highly non-linear third order ODEs, therefore, there is the possibility of the presence of multiple solutions. In this regard, we solved these equations

Discussion
In the following section, we discuss the numerical results of Equations (8) and (9) with their boundary conditions (10). These equations are highly non-linear third order ODEs, therefore, there is the possibility of the presence of multiple solutions. In this regard, we solved these equations numerically by using of shooting method with the fourth order RK technique in Maple software with the help of the shootlib function. It is noticed that there exist triple solutions in the ranges of the different physical parameters. It is worth indicating that triple solutions depend on the opposing flow case where A should be equal to 0.1, dual solutions exist when λ 1 = 0 and A > 0.1, and a single solution exists when λ 1 > 0. In this investigation, we tried to obtain all possible solutions of the considered problem over the stretching and shrinking surfaces; therefore, we kept A = 0.1 and λ 1 = −0.2. It should be noted 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.
We compared − 1 + 1 β f (0) findings with those of Hussain et al. [34] in order to verify the numerical coding of the analysis for various values of β and M in Table 1. It can be concluded from the comparison that our numerical method and coding are working properly and can be used effectively.  Table 2 is given for the smallest eigenvalue values of the solutions. The three-stage Labatto III-A method is adopted to get the smallest eigenvalues by solving Equations (20) and (21) with a relaxed boundary ( F 0 (η) → 0 as η → ∞ is converted to F 0 (0) = 1) condition (22). It is observed that only the first solution is stable, since the values of γ 1 are in the sign of plus, which shows the initial decay of the disturbance, while two remaining solutions are the unstable solutions as the values of γ 1 are in the sign of minus, which indicate the initial growth of the disturbance. The effect of A on the field of velocity f (η) is presented in Figure 3. Figure 3 shows that velocity f (η) is a decreasing (growing) function of A for the shrinking (stretching) case in the first solution. It is detected for the second solution that initially velocity field f (η) rises and then declines as A rises in the case of the shrinking surface, while opposing behavior is noticed in the case of the stretching surface. The third solution exits only when A = 0, 0.1-the dual nature of the velocity field on both surfaces can easily be noticed. Symmetry 2020, 12, x FOR PEER REVIEW 8 of 16 Figures 4 and 5 are presented to notice an effect of the Casson parameter ( ) on the velocity field and temperature field, respectively. From Figure 5, it is noted that the increments in decrease the velocity field and the corresponding thickness of the boundary layer on both surfaces in the first solution. Moreover, the velocity filed increases at first and then begins to decrease in the second and the third solutions when 1, while decreasing behavior is noticed when 1. For the third solution. It seems that there is no big variation in the temperature field for higher values of for first solution for both surfaces (see Figure 5). The temperature of the fluid increases (decreases) over the shrinking (stretching) surface in the second solution when increases. This might happen since the rise in upsurges (decreases) the forces of viscous, and forces produce (impede) some heat energy in the flow. Hence, the temperature distribution with the corresponding thickness of the thermal boundary layer enhances (reduces) with advance values of . For the third solution, the temperature reduces (initially decreases and then increases) over the shrinking (stretching) surface.   Figure 5, it is noted that the increments in β decrease the velocity field and the corresponding thickness of the boundary layer on both surfaces in the first solution. Moreover, the velocity filed increases at first and then begins to decrease in the second and the third solutions when λ = −1, while decreasing behavior is noticed when λ = 1. For the third solution. It seems that there is no big variation in the temperature field for higher values of β for first solution for both surfaces (see Figure 5). The temperature of the fluid increases (decreases) over the shrinking (stretching) surface in the second solution when β increases. This might happen since the rise in β upsurges (decreases) the forces of viscous, and forces produce (impede) some heat energy in the flow. Hence, the temperature distribution with the corresponding thickness of the thermal boundary layer enhances (reduces) with advance values of β. For the third solution, the temperature reduces (initially decreases and then increases) over the shrinking (stretching) surface.
The impact of a magnetic parameter (M) on the fields of velocity and temperature is revealed in Figures 6 and 7. It is perceived that when the effect of M increases, as a result a decrement in velocity profile over both surfaces in the first solution occur. An increase in M yields a resistive kind of force known as Lorentz force, which is created in the flow-as an effect, velocity field curves decrease. The same behaviors of declination of the momentum boundary layer thickness are noticed in the second and third solutions for the case of stretching surface. When λ = −1, the velocity field increases in unstable solutions. From Figure 7, it can be seen that expansion in M raises the temperature distribution over both surfaces in all solutions. Physically, it can be justified as "some extra heat is produced in the flow due to the Lorentz force". Henceforth, an increment in the field of the magnetic parameter amplifies the thickness of the thermal boundary layer. the third solutions when 1, while decreasing behavior is noticed when 1. For the third solution. It seems that there is no big variation in the temperature field for higher values of for first solution for both surfaces (see Figure 5). The temperature of the fluid increases (decreases) over the shrinking (stretching) surface in the second solution when increases. This might happen since the rise in upsurges (decreases) the forces of viscous, and forces produce (impede) some heat energy in the flow. Hence, the temperature distribution with the corresponding thickness of the thermal boundary layer enhances (reduces) with advance values of . For the third solution, the temperature reduces (initially decreases and then increases) over the shrinking (stretching) surface.  The impact of a magnetic parameter ( ) on the fields of velocity and temperature is revealed in Figures 6 and 7. It is perceived that when the effect of increases, as a result a decrement in velocity profile over both surfaces in the first solution occur. An increase in yields a resistive kind of force known as Lorentz force, which is created in the flow-as an effect, velocity field curves decrease. The same behaviors of declination of the momentum boundary layer thickness are noticed in the second and third solutions for the case of stretching surface. When 1, the velocity field increases in unstable solutions. From Figure 7, it can be seen that expansion in raises the temperature same behaviors of declination of the momentum boundary layer thickness are noticed in the second and third solutions for the case of stretching surface. When 1, the velocity field increases in unstable solutions. From Figure 7, it can be seen that expansion in raises the temperature distribution over both surfaces in all solutions. Physically, it can be justified as "some extra heat is produced in the flow due to the Lorentz force". Henceforth, an increment in the field of the magnetic parameter amplifies the thickness of the thermal boundary layer.   Figures 8 and 9 are demonstrated to check the impact of Prandtl ( ) on profiles of velocity and temperature field, respectively. In general, the Prandtl number does not affect the velocity field, except when the flow is flowing over the vertical surface. For the shrinking (stretching) surface, it is found that the velocity field decreases in the first (second) solution for the higher magnitudes of . Moreover, it is determined that the field of velocity is the rising function of in the second and third solutions for the case of a shrinking surface and the third solution for the case of a stretching surface. From Figure 9, the temperature of the fluid drops concerning the higher effect of in all solutions and both surfaces, as is expected. Physically, this happens since the Prandtl number ( ) is the ratio of the momentum diffusion to thermal, which implies a lower thermal effect for a higher effect of .  Figures 8 and 9 are demonstrated to check the impact of Prandtl (Pr) on profiles of velocity and temperature field, respectively. In general, the Prandtl number does not affect the velocity field, except when the flow is flowing over the vertical surface. For the shrinking (stretching) surface, it is found that the velocity field decreases in the first (second) solution for the higher magnitudes of Pr. Moreover, it is determined that the field of velocity is the rising function of Pr in the second and third solutions for the case of a shrinking surface and the third solution for the case of a stretching surface. From Figure 9, the temperature of the fluid drops concerning the higher effect of Pr in all solutions and both surfaces, as is expected. Physically, this happens since the Prandtl number (Pr) is the ratio of the momentum diffusion to thermal, which implies a lower thermal effect for a higher effect of Pr.
Moreover, it is determined that the field of velocity is the rising function of in the second and third solutions for the case of a shrinking surface and the third solution for the case of a stretching surface. From Figure 9, the temperature of the fluid drops concerning the higher effect of in all solutions and both surfaces, as is expected. Physically, this happens since the Prandtl number ( ) is the ratio of the momentum diffusion to thermal, which implies a lower thermal effect for a higher effect of .    where i = 1, 2, 3 and S ci is the critical value where solutions also exist. It has been perceived that the f (0) increases (declines) with increasing suction and Casson parameters in the first (second and third) solution. It is worth mentioning that the shear stress is negative (positive) in the third (first and second) solution. A negative sign suggests that the sheet applies a dragging strength on fluid, and a positive sign suggests the opposite. Moreover, the rise in β leads to an increase in −θ (0) for the third and first solutions, while it decreases in the second solution (see Figure 11).  There exist triple solutions and a single solution when  and  , respectively,  where 1,2,3 and is the critical value where solutions also exist. It has been perceived that the ′′ 0 increases (declines) with increasing suction and Casson parameters in the first (second and third) solution. It is worth mentioning that the shear stress is negative (positive) in the third (first and second) solution. A negative sign suggests that the sheet applies a dragging strength on fluid, and a positive sign suggests the opposite. Moreover, the rise in leads to an increase in ′ 0 for the third and first solutions, while it decreases in the second solution (see Figure 11).             Meanwhile, it decreases in the second solution (see Figure 14). Physically, this is due to the fact that Equation (10) contains a negative sign in front of ; therefore, a higher value of Biot number causes opposite convective heat to be advanced in the thermal equilibrium. Biot number can be explained as "it is a measure of the ratio of surface internal and external thermal resistances". It should be noted that when 0, the sheet's lower side contains hot fluid which is entirely isolated, and therefore the inner thermal resistance of the sheet is tremendously high and no convective heat transfer to the cold fluid over the surface happens. Meanwhile, it decreases in the second solution (see Figure 14). Physically, this is due to the fact that Equation (10) contains a negative sign in front of Bi; therefore, a higher value of Biot number causes opposite convective heat to be advanced in the thermal equilibrium. Biot number can be explained as "it is a measure of the ratio of surface internal and external thermal resistances". It should be noted that when Bi = 0, the sheet's lower side contains hot fluid which is entirely isolated, and therefore the inner thermal resistance of the sheet is tremendously high and no convective heat transfer to the cold fluid over the surface happens.

Conclusions
The flow of steady, stagnation point, MHD of a Casson fluid on a vertical exponentially stretching/shrinking sheet is examined. The effect of the Biot number is also taken into account. The triple solutions are noticed. The significant outcomes of the current study are briefed as follows: 1.
Triple solutions for the coefficient of skin friction, the gradient of temperature, velocity, and temperature profiles occur for specific values of the applied quantity examined in the current examination.

2.
The critical value and the range of the first and second solutions for the coefficient of skin friction rise with a higher magnitude of the Casson parameter.

3.
For the stable solution, the velocity of the Casson fluid reduces (as expected) over both surfaces for the strong field of the Lorentz force.

4.
For the shrinking surface, additional mass suction is required for the occurrence of single and multiple solutions for non-Newtonian Casson fluid (S c1 , S c2 corresponding to β = 1.5, 5) compared to Newtonian fluid (S c3 when β = ∞).
Fluid temperature reduces in all solutions and both surfaces when the effect of Pr increases. 7.
The only first solution is stable from triple solutions.