Stability Analysis of the Magnetized Casson Nanoﬂuid Propagating through an Exponentially Shrinking / Stretching Plate: Dual Solutions

: In this research, we intend to develop a dynamical system for the magnetohydrodynamic (MHD) ﬂow of an electrically conducting Casson nanoﬂuid on exponentially shrinking and stretching surfaces, in the presence of a velocity and concertation slip e ﬀ ect, with convective boundary conditions. There are three main objectives of this article, speciﬁcally, to discuss the heat characteristics of ﬂow, to ﬁnd multiple solutions on both surfaces, and to do stability analyses. The main equations of ﬂow are governed by the Brownian motion, the Prandtl number, and the thermophoresis parameters, the Schmid and Biot numbers. The shooting method and three-stage Lobatto IIIa formula have been employed to solve the equations. The ranges of the dual solutions are f wc 1 ≤ f w and λ c ≤ λ , while the no solution ranges are f wc 1 > f w and λ c > λ . The results reveal that the temperature of the ﬂuid increases with the extended values of the thermophoresis parameter, the Brownian motion parameter, and the Hartmann and Biot numbers, for both solutions. The presence of dual solutions depends on the suction parameter. In order to indicate that the ﬁrst solution is physically relevant and stable, a stability analysis has been performed.


Introduction
There are many kinds of non-Newtonian fluids; namely, Casson fluid, Carreau fluid, Williamson fluid, Maxwell fluid, Micropolar fluid, Jeffrey fluid, Second and third-grade fluids, Burgers' fluid, Eyring-Powell fluid, and so forth [1][2][3][4][5][6][7]. Non-Newtonian fluids are those which are composed of blends. Examples of this type of fluid are pastes, blood, slurries, ketchup, polymer solutions and gels. In simple words, fluids that don't comply with Newton's law of viscosity are non-Newtonian fluids, where shear stress is not directly proportional to the gradient of velocity. However, most previous attempts have been made to derive single solutions of non-Newtonian fluids over different surfaces. In the current research, the MHD flow of a non-Newtonian nanofluid in the presence of velocity and concentration slip parameters has been investigated, along with the effect of the convective boundary condition on exponentially stretching and shrinking sheets. It is worth mentioning that the conduct

Mathematical Formulation
The flow of an electrically conducting two-dimensional (2D) Casson nanofluid on exponentially shrinking and stretching surfaces is assumed. The state rheological equations for the isotropics of Casson fluid are expressed as [8][9][10] where P y , µ B are the yield stress and plastic dynamic viscosity of Casson fluid, respectively. Moreover, π c is a critical value of π based on the non-Newtonian model, and e ij is the (i, j)th component of the deformation rate, expressed as π = e ij e ij . The surface velocity is presumed as u w = λae x l , where a is the characteristic velocity of the surface and λ < 0 and λ > 0 indicate the shrinking and stretching surface, respectively. Further, we supposed that the variable magnetic field has a strength of B = B 0 e x 2l , where B 0 is the constant magnetic field. The induced magnetic field is neglected by assuming the number of low magnetic Reynolds.
The field of the velocity of the boundary layers of the Casson nanofluid's flow, with temperature and concentration equations, can be expressed as The related boundary conditions (2)(3)(4)(5) are where u is the velocity component along the x-axis, v is the velocity component along the y-axis, P y , ρ, σ * , α, D B , D T are the viscosity of the fluid, the Casson parameter, the density of fluid, the electrical conductivity, the thermal diffusivity, the Brownian motion and the thermophoresis, respectively. τ w is the heat capacity of the nanofluid and the effective heat capacity of the nanoparticle material, and T ∞ and C ∞ are ambient temperature and concentration correspondingly, such that is the velocity condition, where A * 1 is the slip factor of velocity, and N = N 1 e −x 2l is the concentration condition where N 1 is the slip factor of concentration. We will employ similar variables (7) in Equations (2)-(6) in order to obtain similar solutions.
where the stream function has the following relations with the velocity components u = ∂ψ ∂y , v = − ∂ψ ∂x . By substitution of the stream function relations with Equation (7) in (2)-(6), we obtain 1 Pr along with the boundary conditions Here, prime denotes the derivative with respect to η, is the magnetic field, Pr = ϑ α is the Prandtl number, and N t = are the thermophoresis and Brownian motion parameters, respectively. Furthermore, Sc = ϑ D B denotes the Schmidt number, λ is the stretching/shrinking parameter, v w = − ϑa 2l e x 2l f w is the suction/blowing parameter, δ = A * 1 ϑa 2l is the velocity slip parameter, δ C = N 1 a 2ϑl is the concentration slip parameter, and the Biot number As regards the coefficient of skin friction, the local Nusselt and Sherwood numbers are physical quantities of interest, and can be defined as follows where Re x = ax ϑ is the Reynold number.

Stability Analysis
Weidman et al. [33] proposed stability solutions, and Haris and Pop [34] proceeded with them in their studies, wherein multiple solutions for the boundary layer have been considered. From their work, it is determined that the additional solution is unrealizable, and said to be an unstable solution. On the other hand, the first solution is the physically realizable solution, which is said to be a stable solution. The initial step of obtaining the stable solution is to change Equations (2)-(5) to the unsteady governing equations, by presenting a new dimensionless variable, τ = a 2l e x l .t, as suggested by Merkin [35]. Thus, we have ∂u ∂t ∂C ∂t It is interesting to highlight that the presence of τ is associated with initial value problems, which correspond to the stable solution. After the integration of τ = a 2l e x l .t with η = y a 2ϑl e x 2l , we have following new similarity transformation variables Replacing Equations (13)- (15) with Equation (16), we obtain 1 Pr subject to the boundary conditions Subsequently, to classify the stability of the solution f (η) = f 0 (η), θ(η) = θ 0 (η) and ∅(η) = ∅ 0 (η), which solved the boundary value problems of Equations (17)-(20), we will follow the steps of Rosca and Pop [36]; where F(η, τ), G(η, τ) and H(η, τ) are corresponding small relatives of f 0 (η), θ 0 (η) and ∅ 0 (η), and ε is the unknown eigenvalue parameter which must be found. By applying Equation (21) in Equations (17)- (19), including boundary condition (20) and keeping the value of τ = 0, as recommended in paper [33] to compute the initial decay and growth of the solutions of (21), the functions F(η, τ), G(η, τ) and H(η, τ) are expressed as F 0 (η), G 0 (η) and H 0 (η), respectively. It should be noted that these functions' ranges must be 0 < F(η, τ) < 1, 0 < G(η, τ) < 1 and 0 < H(η, τ) < 1. Henceforth, we have a system for the linearized eigenvalue problem. 1 with the boundary condition According to the previous study of Haris et al. [34], in order to find the smallest possible eigenvalues, we need to relax one boundary condition of the functions F 0 (η), G 0 (η) and H 0 (η). In the current problem, we relaxed the F 0 (η) → 0 as η → ∞ , and then solved the Equations (22)-(25) with the new relaxed initial conditions F 0 (0) = 1. It is worth mentioning that if the value of ε > 0, then the initial decay of disturbances exists, and the flow becomes stable and physically realizable. Further, if the value of ε < 0, this indicates the initial growth of disturbances in the system, and the flow becomes unstable.

Numerical Methods
In this section, the numerical procedures of the shooting method and the three-stage Lobatto IIIa formula are described. These methods are very simple and highly accurate in computing results for the boundary layer flow problems, and many researchers use them to obtain multiple solutions and the stability analyses of the solutions.

Shooting Method
The shooting method is an iterative method that changes the boundary value problems (BVPs) into equivalent the initial value problems (IVPs) by using a sequence of initial condition guesses, obtained using an appropriate iterative method, until the required accuracy is reached. The resultant IVPs are then solved by employing the suitable IVP method presented in the literature. In this problem, we have used the well-known Runge-Kutta method of the fourth order, which is known because of its high accuracy. The descriptions of the method are as follows where the unknown initial conditions are denoted by α 1 , α 2 and α 3 . It is worth mentioning that the missing initial values of α 1 , α 2 and α 3 are needed to shoot such that the profiles of the solutions satisfy the original boundary conditions, which are f (η) → 0, θ(η) → 0, and ∅(η) → 0 as η → ∞ .

Three-Stage Lobatto IIIA Formula
This is a well-known numerical method. All types of linear and non-linear differential equations are solved by this method easily. The three-stage Lobatto IIIa formula is created in bvp4c with the help of a finite difference code. Afterward, stability analysis is performed using the bvp4c solver function. According to Rehman et al., [37] "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". Further, the tolerance of the relative error is fixed at 10 −5 . The suitable mesh determination is created and returned in the field sol.x. The bvp4c returns the solution, called 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 method, along with stability analysis, is presented in Figure 1.
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". Further, the tolerance of the relative error is fixed at 10 . The suitable mesh determination is created and returned in the field sol.x. The bvp4c returns the solution, called 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 method, along with stability analysis, is presented in Figure 1.

Result and Discussion
Providing clear comprehension of the present problem, we performed numerical computation using two numerical methods; specifically, the shooting technique along with the Runge-Kutta method of the fourth order in MAPLE software, and the three-stage Lobatto IIIa formula in bvp4c by employing the finite difference code in MATLAB software. We attempted to plot all the graphs of skin friction ( (0)), heat (− (0)) and concentration (−∅ (0)) transfer rates, and the temperature, velocity and concentration profiles for numerous values of directly involved parameters, such as the Casson ( ) and magnetic parameters ( ) , the Schmidt ( ) and Prandtl ( ) numbers, the thermophoresis ( ) and Brownian motion ( ) parameters, the suction ( ) and velocity slip ( ) parameters, and the convection parameter ( ) and concentration slip parameter ( ). To validate the correctness of the shooting method, the values of (0), − (0) and −∅ (0) are compared with those derived by Rehman et al. [37] and Mustafa et al. [38] in Table 1, and we establish excellent concurrences with them; henceforth, we can use the present code confidently. The system of highly non-linear ordinary differential equations (ODEs) (8-10) subject to boundary condition (11) is solved by employing the shooting method with the Runge Kutta (RK)-4th order method in MAPLE software. The shooting technique is an iterative method which changes BVPs into equivalent IVPs by using a sequence of initial condition guesses, obtained using an appropriate iterative method until the required accuracy is reached, and then solved by the shooting technique. This is the accepted method

Result and Discussion
Providing clear comprehension of the present problem, we performed numerical computation using two numerical methods; specifically, the shooting technique along with the Runge-Kutta method of the fourth order in MAPLE software, and the three-stage Lobatto IIIa formula in bvp4c by employing the finite difference code in MATLAB software. We attempted to plot all the graphs of skin friction ( f (0)), heat (−θ (0)) and concentration (−∅ (0)) transfer rates, and the temperature, velocity and concentration profiles for numerous values of directly involved parameters, such as the Casson (β) and magnetic parameters (M), the Schmidt (Sc) and Prandtl (Pr) numbers, the thermophoresis (N t ) and Brownian motion (N b ) parameters, the suction ( f w ) and velocity slip (δ) parameters, and the convection parameter (A) and concentration slip parameter (δ C ). To validate the correctness of the shooting method, the values of f (0), −θ (0) and −∅ (0) are compared with those derived by Rehman et al. [37] and Mustafa et al. [38] in Table 1, and we establish excellent concurrences with them; henceforth, we can use the present code confidently. The system of highly non-linear ordinary differential equations (ODEs) (8-10) subject to boundary condition (11) is solved by employing the shooting method with the Runge Kutta (RK)-4th order method in MAPLE software. The shooting technique is an iterative method which changes BVPs into equivalent IVPs by using a sequence of initial condition guesses, obtained using an appropriate iterative method until the required accuracy is reached, and then solved by the shooting technique. This is the accepted method for solving the boundary value problems of ODEs. We kept the minimum and maximum range of η ∞ = 6 and η ∞ = 10, respectively, which satisfied the hydrodynamic, temperature and concentration boundary layers asymptotically.
It is worth mentioning that the solution to this problem is the local similar solution. Figure 2 is drawn to show the presence of two solutions in both surfaces, namely the shrinking and stretching surfaces. Figures 3-5 show the effect of λ on the shear stress f (0), the local heat flux −θ (0) and the local mass flux −∅ (0). From these figures, it can be concluded that two sorts of solutions exist; explicitly dual solutions (λ c ≤ λ) and no solution (λ c > λ). It has been observed that the skin friction coefficient is higher for the shrinking surface as compared to the stretching surface in both solutions (refer Figure 3). Further, the behaviors of the heat and concentration transfer seem the same, where −θ (0) and −∅ (0) increase over the stretching surface, as compared to the shrinking surface in the first solution. On the other hand, dual behaviors are noticed in the second solution. The impact of the Casson parameter (β) on the shear stress with various values of the suction parameter ( f w ) is drawn in Figure 6. From studies of Miklavčič and Wang [39], Fang and Zhang [40] and Bhattacharyya [41], it can be concluded that "flow due to the exponentially shrinking surface needs more suction when contrasted with the flow over a linearly shrinking surface so as to keep the larger amount of vorticity produced inside the boundary layer". For the case of a Newtonian fluid, β = ∞, the range of dual solution is f wc1 ≤ f w , and the no solution range is f wc1 > f w . For Casson fluid cases, multiple solutions rely upon the values of the Casson parameter; for example, high suction is expected to keep up the flow inside the boundary layer for small values of the Casson parameter, as shown in Figure 6. On the other hand, the Casson parameter is relative to the shear stress for the first solution, where the coefficient of skin friction rises when β diminishes. For the second solution, the reverse trend was noticed. Figures 7  and 8 display the impact of β on local Nusselt and Sherwood numbers separately. The local Nusselt and Sherwood numbers are higher for the greater value of β, contrasted with a smaller one, which is because we have 1 β in the momentum equation. Table 1.                    The impact of β on the velocity profile appears in Figure 9. The momentum of the boundary layer improves the first solution when 1 β changes in the absence or the presence of the magnetic parameter. On the other hand, the double nature of the velocity distribution and its boundary layer thickness is noticed in the second solution case when β rots; physically, this tends to be explained by the fact that decreases in β are caused in order to produce greater viscosity, and as a result resistance is produced in the flow of fluid. Consequently, the hydrodynamic boundary layer becomes thinner. Figure 10 exhibits the impact of M on the velocity distribution with and without the velocity slip effect. The slip effect reduces the velocity boundary layer monotonically in both solutions, when contrasted with δ = 0. The momentum boundary layer ends up slenderer in the first solution, and the physical magnetic parameter creates progressively more resistance in the flow of the fluid due to the moving electric charges in nanofluid, and thus the velocity of the fluid deaccelerated. On the other hand, the opposite pattern is perceived for the second solution. Figure 11 demonstrates the effect of the Casson parameter (β) on temperature distribution. The reduction in β produces a contrary force in the flow of fluid due to the high viscosity in the fluid flow, thus the temperature of the fluid and the thermal boundary layer thickness decrease continually when β declines in both solutions. Figure 12 presents the effect of different Prandtl numbers on the distribution of temperature. The higher convectional number (A) improves the temperature of the surface and the boundary layer becomes thicker in both solutions, as expected. It is worth noting that A → ∞ indicates the constant temperature of the wall. Further, we know the Prandtl number (Pr) is the ratio of momentum diffusivity to thermal diffusivity, which implies an enhancement in Pr is the cause of the flimsier thermal diffusivity, and therefore the thickness of the thermal boundary layer becomes thinner, as found in Figure 12. The impact of the Brownian motion (N b ) and thermophoresis (N t ) parameters on the temperature profile is exhibited in Figures 13 and 14, respectively. The thicker thermal boundary layer is present as the Brownian motion and thermophoresis parameters rise in all the solutions. Practically, thermophoresis and the Brownian motion help the nanoparticles to change their places from hot to cold regions; in this manner, the thickness of the thermal boundary layer is enhanced monotonically. Figures 15 and 16 reveal the impact of N t and N b on concentration distributions. The concentration profile causes a decline in the Brownian motion parameter (N b ), which is improved in both solutions (refer to Figure 15). On the contrary, the thermophoresis parameter (N t ) is provably relative to the concentration profile. It is additionally seen that the effect of the concentration slip (δ C ) is the reason for the development of the concentration boundary layer in both solutions (alluded to in Figures 15 and 16). The concentration profile for numerous values of the Schmidt number (Sc) is drawn in Figure 17. Lower molecular diffusivity is found for a higher Schmidt number (Sc). Henceforth, a thicker concentration boundary layer in both solutions is possible with the lower values of the Schmidt number (Sc).
In this study, to measure the stability of the solution, stability analyses are performed. The fact behind this investigation is that if there occurs more than one solution in any problem. The main objective of performing this examination is to determine what solution is the first solution, which is linearly stable and physically related, by utilizing the bvp4c solver in MATLAB programming. The solution's stability is based on the sign of the minimum value of ε. To find the unknown eigenvalue (ε) from Equation (21), we have to solve Equations (22)-(24) with relaxed F 0 (0) = 1 boundary conditions (25). Table 2 displays the values of ε 1 for different values of f w and λ. From Table 2, it is obvious to us that the negative values of ε 1 are demonstrated in the second solution, which shows the initial growth of disturbance in the flow, and this is supposed to be an unstable solution. On the other hand, positive values of ε 1 are shown in the first solution, which signifies the initial decay of disturbance; thus the flow becomes stable.

Conclusion Remarks
In the present research paper, the MHD flow of an electrically conducting Casson nanofluid on an exponentially shrinking and stretching sheets, with the effects of velocity and concertation slips along with the convective condition, is investigated numerically. By using exponential similarity variables, governing equations are reduced into highly non-linear ODEs, which are solved by the shooting method and three-stage Lobatto IIIa formula. In order to find a stable solution, stability analysis is also conducted. These are the main results of the problem:

1.
For both surfaces, there are dual solutions; 2.
For the stable solution, the velocity profile decreases for higher values of Casson parameter β and the magnetic field M; 3.
With a rise in the intensity of the convection parameter A, the thermal boundary layer is enhanced. Yet, the Prandtl number Pr has the inverse relationship with the temperature profile in both solutions.

4.
Stability analysis reveals that an unstable (/stable) solution is a second (/first) solution.

5.
Dual solutions vary, and no solution depends on the parameters involved.