Non-Unique Solutions of Magnetohydrodynamic Stagnation Flow of a Nanofluid towards a Shrinking Sheet Using the Solar Radiation Effect

This study aims to investigate the magnetohydrodynamic flow induced by a moving surface in a nanofluid and the occurrence of suction and solar radiation effects using the Buongiorno model. The numerical findings are obtained using MATLAB software. The effects of various governing parameters on the rates of heat and mass transfer along with the nanoparticles concentration and temperature profiles are elucidated graphically. Non-unique solutions are discovered for a specific variation of the shrinking strength. The temporal stability analysis shows that only one of them is stable as time passes. Furthermore, raising the Brownian motion parameter reduces both the local Sherwood number and the local Nusselt number for both solutions. It is also observed that increasing the thermophoresis parameter reduces the rate of heat transfer, whereas the opposite trend is observed for the rate of mass transfer.


Introduction
Renewable energy is becoming more popular as a way to minimize reliance on finite fossil fuel supplies and alleviate the effects of climate change. Since it transforms solar energy directly into heat and electricity without emitting greenhouse gases, solar energy has been established as the greatest and most inexpensive renewable energy source [1]. Studies during the past few years on solar energy have indicated that nanofluids play a major role in enhancing efficiency in solar systems due to their outstanding ability to transfer heat [2][3][4][5]. Choi and Eastman [6] pioneered the discovery of nanofluids by combining nanometer-sized particles with conventional fluids. Nanofluids are superior coolants in transportation, nuclear reactors, lubricants, thermal storage, domestic refrigerators, optical devices, and cancer therapeutics since nanoparticles have more thermal conductivity compared to base fluid [7][8][9][10]. Two models have been proposed by Tiwari and Das [11] and Buongiorno [12] for investigating the nanofluid properties. It is notable to mention that the model by Buongiorno is a non-homogeneous two-phase model with a slip velocity for the base fluid as well as nanoparticles that are not equal to zero. This model has several slip mechanisms: Brownian diffusion, gravity, thermophoresis, diffusiophoresis, fluid drainage, the Magnus effect, and inertia. Meanwhile, the Tiwari and Das model is a homogeneous single-phase model that investigates the influence of the solid volume fraction of nanoparticles on nanofluid behavior. Multiple scholars have utilized these two models to examine the characteristics of flow and heat transmission of a nanofluid under different physical situations [13][14][15][16].
Magnetohydrodynamic (MHD) is a field of study that investigates the magnetic field effect on electrically conductive fluids. In recent times, the study of magnetohydrodynamic has attracted numerous scholars given its wide implementations in engineering and industry, such as nuclear reactors, crystal growth, optical fiber filters, optical grafting, liquid metals, petroleum industries, and metallurgical processes [17]. Lund et al. [18] investigated the MHD mixed convection flow over a shrinking/stretching sheet in a nanofluid with a viscous dissipation effect. The results revealed that various solutions were produced; nevertheless, the first solution was determined to be the most stable. With the occurrence of injection/suction parameters, Naramgari and Sulochana [19] investigated the chemical reaction and thermal radiation effects on the constant magnetohydrodynamic flow of a nanofluid driven by a permeable shrinking/stretching surface. They discovered that as thermophoresis parameters and Brownian motion enhance, the nanoparticle concentration decreases while the mass transfer rate increases. The magnetic field influence on the free convection flow in a corrugated container loaded with nanofluids was investigated by Ahmed and Rashed [20]. The topic of hydromagnetic flow of third-grade nanofluid caused by a nonlinear deformable sheet under the impact of activation energy and chemical reactions was considered by Hayat et al. [21]. Afterward, Ayub et al. [22] used the Lorentz force impact to study the Carreau nanofluid flow, taking into account the influence of infinite shear rate viscosity. Ayub et al. [23] researched the infinite shearing rate influence on MHD Carreau nanofluid flow induced by a cylindrical channel. They discovered that increasing the impact of an inclined magnetic dipole tends to raise the energy while decreasing the velocity.
The stagnation flow demonstrates the behavior of fluid movement close to the stagnation point. Considering its wide applications in both scientific and industrial areas, the research of stagnation point flow is of great significance. In addition, Hiemenz [24] was the earliest to investigate the flow past a flat sheet at a two-dimensional stagnation point. Moreover, Dzulkifli et al. [25] examine the heat transfer rate and unstable stagnation-point flow towards an exponentially stretching surface in nanofluids, considering the velocity slip effect. Yashkun et al. [26] investigated the features of the stagnation point flow and heat transfer across a shrinking/stretching sheet with a suction effect. Subsequently, Kamal et al. [27] investigated the mixed effects of injection/suction, chemical reaction, and magnetic field on the flow in a nanofluid caused by a permeable shrinking/stretching surface. Finally, the constant two-dimensional magnetohydrodynamic stagnation flow of a nanofluid via a shrinking sheet was considered by Aladdin et al. [28]. It was inferred during the examination that when the magnetic field was enhanced, the fluid velocity, heat, and mass transmission rates increased. Still, the rate of heat and mass transfer is reduced when the thermophoresis parameter is enhanced. Since then, a great number of studies on the stagnation point flow have been completed in many respects in fluid dynamics, for instance [29][30][31][32][33][34].
Existing literature confirmed that no investigation has been done before on the magnetohydrodynamic flow of a nanofluid towards a shrinking sheet in the presence of solar radiation and suction using the Buongiorno model. To fill the above knowledge gap, a mathematical formulation of the current study is created with reference to the study conducted by Ghasemi and Hatami [35]. This research is different from that considered by Ghasemi and Hatami [35], where we consider a shrinking sheet with suction effects. The impact of the various physical parameters on the Sherwood number and the Nusselt number, along with the nanoparticles' concentration and temperature profiles, are depicted in the tabular form and graph. Besides, the emergence of non-unique solutions as well as their stability are discussed.

Mathematical Formulation
The two-dimensional magnetohydrodynamic flow towards a stretching/shrinking surface with solar radiation and suction effects is considered, as illustrated in Figure 1. It is assumed that u w (x) = ax is the surface velocity and u ∞ (x) = bx is the free stream It is noted that a > 0 is for stretching and a < 0 is for shrinking. Moreover, the ambient fluid temperature and the convective surface temperature are signified as T ∞ and T f , respectively. A magnetic field of uniform strength B 0 is implemented parallel to the direction of the flow. Given the previous assumptions, the governing equations can be written as [35]: subject to the boundary conditions:

Mathematical Formulation
The two-dimensional magnetohydrodynamic flow towards a stretching/shrinking surface with solar radiation and suction effects is considered, as illustrated in Figure 1. It is assumed that ( ) = is the surface velocity and ∞ ( ) = is the free stream velocity, where and are the positive constants. It is noted that > 0 is for stretching and < 0 is for shrinking. Moreover, the ambient fluid temperature and the convective surface temperature are signified as ∞ and , respectively. A magnetic field of uniform strength 0 is implemented parallel to the direction of the flow. Given the previous assumptions, the governing equations can be written as [35]: subject to the boundary conditions: Here and denote the component of velocity in the and axes, respectively, ( ) implies the mass flux velocity, represents the kinematic viscosity, signifies the electrical conductivity, and 0 indicates constant magnetic field in the direction. The following similarity transformation is employed to solve Equations (1)-(3): By utilizing the similarity transformation (4), the continuity Equation (1) is satisfied, while Equations (2) and (3) reduce to where = 0 2 ⁄ indicates the magnetic parameter, = ⁄ symbolizes the stretching/shrinking parameter, and = − √ ⁄ signifies the suction/injection parameter. On the other hand, the energy and concentration equations are respectively given by Here u and v denote the component of velocity in the x and y axes, respectively, v w (x) implies the mass flux velocity, ν f represents the kinematic viscosity, σ e signifies the electrical conductivity, and B 0 indicates constant magnetic field in the y direction. The following similarity transformation is employed to solve Equations (1)- (3): By utilizing the similarity transformation (4), the continuity Equation (1) is satisfied, while Equations (2) and (3) reduce to where M = σ e B 2 0 /ρ f b indicates the magnetic parameter, ε = a/b symbolizes the stretching/shrinking parameter, and s = −v w / √ aν f signifies the suction/injection parameter.
On the other hand, the energy and concentration equations are respectively given by where T represents the temperature, C f the allusion to the specific heat capacity, α implies the thermal diffusivity, and C denotes the nanoparticle concentration. The thermophoretic diffusion coefficient and the Brownian motion coefficient are signified by D T and D B , respectively. Further, q r addresses the radiative heat flux and τ = (ρc) p /(ρc) f refers to the proportion between the nanoparticle effective heat capacity and the heat capacity of the fluid. By utilizing the Rosseland approximations, we get the radiative heat flux in the following form (see Raptis [36], Brewster [37], and Sparrow and Cess [38]).
In Equation (9), σ * represents Stefan-Boltzman constant, while k * denotes the mean absorption coefficient. By utilizing Equation (9), Equation (7) becomes The boundary constraints are expressed as The non-dimensional temperature profile can be defined as 3 where R d = 16σ * T 3 ∞ /3kk * stands for the radiation parameter, and there is no thermal radiation effect when R d = 0. The last expression is then reduced to the following form: where Pr = ν f /α denotes the Prandtl number.
Equations (8) and (10) are transformed to the following form: subject to the transformed boundary conditions: Here ϕ(η) = c − c ∞ /c w − c ∞ and Bi iis the Biot number, N t is the thermophoresis parameter; and N b is the Brownian motion parameter, which are expressed as The surface heat and mass fluxes are expressed as: By using (16), the local Nusselt number Nu x = xq w /k(T f − T ∞ ) and the local Sherwood number where Re x = b/ν f x 2 is the local Reynolds number.

Stability Analysis
The numerical results demonstrate that there are two solutions for specific parameter values. As a result, it is necessary to look into the temporal stability of these solutions, to see which one is stable as time evolves. To perform a stability analysis of the dual solutions, as per Weidman et al. [39], this issue must be regarded in an unsteady form by introducing the new dimensionless time parameter τ. The unsteady case for the governing Equations (2), (8), and (10) is: where t refers to the time, the new similarity transformation is expressed as follows: Applying Equation (21) into Equations (18)- (20), the following equations are obtained: subjected to The following perturbations are introduced to investigate the flow stability in the long run, with f = f 0 (η), θ = θ 0 (η), and ϕ = ϕ 0 (η) which satisfy Equations (22)-(24) (see Weidman et al. [39]), ϕ(η.τ) = ϕ 0 (η) + e −γτ H(η) where γ is an unknown eigenvalue. The linearized problems are obtained by incorporating Equation (26) into Equations (22)- (25): 1 Pr and are subjected to the boundary conditions Here, the least eigenvalue γ determines the stability of the steady flow solutions f 0 (η), θ 0 (η) and ϕ 0 (η). Based on the study by Harris et al. [40], relaxation of the boundary condition on F(η), G(η) or H(η) is required to obtain the possible eigenvalues. Thus, the condition F 0 (η) → 0 as η → ∞ is relaxed, where the system (26)-(29) is solved along with the new boundary condition F 0 = 1.

Findings and Discussion
MATLAB software (MATLAB R2022a, MathWorks, Inc., Natick, MA, USA) is employed to obtain the numerical solutions to the ordinary differential equations (ODEs) (5), (13), and (14) aligned with the boundary conditions (6) and (15). The impacts of the governing parameters on the local Nusselt number Nu x Re −1/2 x (refers to the heat transfer rate) and local Sherwood number Sh x Re −1/2 x (refers to the mass transfer rate), along with the temperature profile θ(η), and concentration profile ϕ(η) are addressed graphically in Figures 2-9. It is prominent to note that the dual solutions are discovered to exist at a specified level of the parameters towards the shrinking region (ε > ε c ) and a unique solution is obtained when ε = ε c . However, no solutions have been discovered for ε < ε c , where ε c signifies the critical value of ε for which solutions exist. Moreover, a comparison of the heat transfer rate as well as the mass transfer rate is done with those recorded by Ghasemi and Hatami [35] and Khan and Pop [41], as shown in Tables 1 and 2. The comparison shows excellent agreement. As a result, the accuracy of the findings reported in this study is assured.      thermal boundary layer and therefore reducing the heat transfer rate at the surface. Contrary to this, Figure 5 displays reverse behavior when the increment of the thermophoresis parameter augments the local Sherwood number. The influence of thermophoresis will help the nanoparticles with high thermal conductivity penetrate further into the fluid, resulting in a thinner concentration boundary layer. Hence, the concentration gradient increases, which enhances the rate of mass transfer at the surface.   trary to this, Figure 5 displays reverse behavior when the increment of the thermophoresis parameter augments the local Sherwood number. The influence of thermophoresis will help the nanoparticles with high thermal conductivity penetrate further into the fluid, resulting in a thinner concentration boundary layer. Hence, the concentration gradient increases, which enhances the rate of mass transfer at the surface.   perature gradient and an elevation in the concentration gradient, resulting in a reduction in the heat transfer rate and an elevation of the mass transfer rate, as observed in Figures  4 and 5. Note that the thermophoresis force acts opposite the gradient of the temperature in nanofluids as well as helping nanoparticles migrate from the hot surface to the cold ambient fluid.  is positive, this signifies that the solution appears stable and indicates that the flow has just slight disturbances that have no effect on the flow characteristics or physical appearance. Otherwise, a negative value for the smallest eigenvalues indicates that the solution is unstable, implying that the disturbance impacting the flow system is growing. Moreover, we may deduce from Figure 10 that the first solution is stable, whereas the other solution is not, as time evolves. this signifies that the solution appears stable and indicates that the flow has just slight disturbances that have no effect on the flow characteristics or physical appearance. Otherwise, a negative value for the smallest eigenvalues indicates that the solution is unstable, implying that the disturbance impacting the flow system is growing. Moreover, we may deduce from Figure 10 that the first solution is stable, whereas the other solution is not, as time evolves.   Nt rises. This is due to the fact that the thermophoresis effect assists nanoparticles in transitioning from hot to cold regions, leading to a thicker thermal boundary layer and therefore reducing the heat transfer rate at the surface. Contrary to this, displays reverse behavior when the increment of the thermophoresis parameter Nt augments the local Sherwood number. The influence of thermophoresis will help the nanoparticles with high thermal conductivity penetrate further into the fluid, resulting in a thinner concentration boundary layer. Hence, the concentration gradient increases, which enhances the rate of mass transfer at the surface.
Moreover, the behavior of the temperature profile θ(η) as well as the concentration profile ϕ(η) with a variety of Brownian motion parameter Nb, where Pr = 6.2, Nt = 0.1,  Figures 8 and 9, respectively. In both solutions, raising the thermophoresis parameter Nt boosts the temperature θ(η) while decreasing the concentration ϕ(η). Additionally, the thermal boundary layer thickness rises as the thermophoresis parameter Nt rises, but the concentration boundary layer thickness drops. This culminates in a reduction in the temperature gradient and an elevation in the concentration gradient, resulting in a reduction in the heat transfer rate and an elevation of the mass transfer rate, as observed in Figures 4 and 5. Note that the thermophoresis force acts opposite the gradient of the temperature in nanofluids as well as helping nanoparticles migrate from the hot surface to the cold ambient fluid.
Lastly, by utilizing the bvp4c function in Matlab to solve Equations (27)- (30), an analysis of stability is executed. In Figure 10, the smallest eigenvalues γ for dual solutions when Pr = 6.2, Nb = 0.3, Nt = 0.1, Rd = 0.1, θ w = 0.3, Le = 3, Bi = 1, M = 0.1, and s = 0.1 are presented for several values of ε. The eigenvalues play a crucial role in determining the dual solutions stability. Provided that the smallest eigenvalue γ is positive, this signifies that the solution appears stable and indicates that the flow has just slight disturbances that have no effect on the flow characteristics or physical appearance. Otherwise, a negative value for the smallest eigenvalues indicates that the solution is unstable, implying that the disturbance impacting the flow system is growing. Moreover, we may deduce from Figure 10 that the first solution is stable, whereas the other solution is not, as time evolves.

Conclusions
The mathematical model of the MHD flow of a nanofluid over a shrinking sheet under influenced by suction and solar radiation has been investigated. By employing the proper similarity variables, the governing PDEs are transformed to ODEs, which are subsequently numerically solved by employing the bvp4c solver in Matlab software. The primary outcomes obtained are summarized below:

•
The outcomes deduce that the existence of dual (non-unique) solutions is provable for a given shrinking strength range ( < −1).

•
According to the temporal stability study, only the first solution is stable and, hence, physically significant.

•
The Sherwood number as well as the Nusselt number decreases when the Brownian motion parameter upsurges.

Conclusions
The mathematical model of the MHD flow of a nanofluid over a shrinking sheet under influenced by suction and solar radiation has been investigated. By employing the proper similarity variables, the governing PDEs are transformed to ODEs, which are subsequently numerically solved by employing the bvp4c solver in Matlab software. The primary outcomes obtained are summarized below:

•
The outcomes deduce that the existence of dual (non-unique) solutions is provable for a given shrinking strength range (ε < −1). • According to the temporal stability study, only the first solution is stable and, hence, physically significant.

•
The Sherwood number as well as the Nusselt number decreases when the Brownian motion parameter upsurges.

•
The rate of heat transfer reduces when the thermophoresis parameter Nt is elevated; however, the rate of mass transfer is found to increase. The concentration augments by incrementing the Brownian motion parameter Nb but reduces by elevating the thermophoresis parameter Nt.

•
The influence of the Brownian motion parameter Nb and the thermophoresis parameter Nt on the temperature profile reveals that the thermal boundary layer thicknesses as well as the temperature increase for both solutions. Data Availability Statement: Not applicable.

Conflicts of Interest:
The authors declare no conflict of interest.
Subscripts w condition at the wall ∞ ambient condition Superscript differentiation with respect to η