Multiple Solutions for Stagnation-Point Flow of Unsteady Carreau Fluid along a Permeable Stretching/Shrinking Sheet with Non-Uniform Heat Generation

The aim of the present study was to explore the effect of a non-uniform heat source/sink on the unsteady stagnation point flow of Carreau fluid past a permeable stretching/shrinking sheet. The novelty of the flow model was enhanced with additional effects of magnetohydrodynamics, joule heating, and viscous dissipation. The nonlinear partial differential equations were converted into ordinary differential equations with the assistance of appropriate similarity relations and were then tackled by employing the Runge-Kutta-Fehlberg technique with the shooting method. The impacts of pertinent parameters on the dimensionless velocity and temperature profiles along with the friction factor and local Nusselt number were extensively discussed by means of graphical depictions and tables. The current results were compared to the previous findings under certain conditions to determine the precision and validity of the present study. The fluid flow velocity of Carreau fluid increased with the value of the magnetic parameter in the case of the first solution, and the opposite behavior was noticed for the second solution. It was seen that temperature of the Carreau fluid expanded with the higher values of unsteadiness and magnetic parameters. It was visualized from multiple branches that the local Nusselt number declined with the Eckert number parameter for both the upper and lower branch.


Introduction
The numerous application of non-Newtonian fluids in industry and commerce has prompted researchers to do study in this area. Important applications of these types of fluids include the chemical industry, such as paint manufacture, palm oil production, and shampoo production, as well as the food sector, such as mayonnaise production. The highly driven authors are therefore interested in the study of the rheology of non-Newtonian liquids. As the complicated numerical and analytical relationship between the shear rate and stress is represented by the non-Newtonian liquid substance, they are graded into dilatant, shear thinning, and shear thickening properties. While different fluid simulations are used in this respect to analyze the inherent advantages from the above components, there is no single scheme in this respect. Finally, the organized effort by Carreau (1972) [1] moved into motion and suggested the Carreau fluid scheme. He identified that the Carreau fluid scheme is a combination of the Newtonian and power law scheme that can express shear thinning characteristics at a reduced shear rate and thickening characteristics at a large shear rate. Good numbers of research articles are reported in literature which deals the characteristics of Carreau fluid flow induced by different stretched surfaces. Some of these are as follows. Olajuwan [2] stated that by enhancing the deformation rate, the constitutive expression of Carreau fluid decreased to the non-Newtonian fluid. Carreau's scheme was Coatings 2021, 11, 1012 2 of 14 described by Pan-tokratoras [3], utilizing the control parameter n. He explained that the characteristics of a shear thinning substance for 0 < n < 1 are expressed by a Carreau fluid, and a shear thickening for n > 1. Shadid and Eckert [4] investigated the dampers impact of a stretching cylinder on Carreau fluid. Khellaf and Lauriat [5] provided thermophysical characteristics of Carreau fluid in an inertial space between two concentric cylinders. Raju and Sandeep [6] were questioned about the effect of the nonlinear exposure on the flow of Carreau fluid. Khan and Hashim [7] illustrated the MHD boundary layer flow of the Carreau fluid owing to a stretched sheet. Some recent development in the heat transfer is found in Refs. [8][9][10][11][12][13][14][15][16][17][18][19].
MHD is the study of electrically conducting fluids that combine electromagnetism and fluid dynamics principles. The field is illustrated by the process of mixing nuclear reactor cooling and metals in an electrical boiler by introducing a magnetic field. Magnetic fields are used to track the heat transfers and momenta of various fluids in the boundary layer flow as they pass through a stretched surface. Some of the uses of MHD include crude oil purification, MHD generators, the petroleum sector, polymer technology, thermal protection, space vehicle propulsion, and MHD pumps. Such recent studies are as follows. Mukhopadhyay et al. [20] studied the slip mechanism of MHD flow of viscous fluid along a stretching cylinder. The MHD flow across the stretched surface with the condition of a convective boundary was examined by Ibrahim [21]. An experimental and numerical analysis was carried out by Hashizume [22] to determine the MHD flow in the liquid metal sheet process by identifying Li as a liquid metal that works. Currently, various researchers [23][24][25][26][27][28] have examined the MHD flow behavior of the lithium material in various duct structure applications.
Many researchers have looked at boundary layer flows caused by stretching/shrinking surfaces because of their ubiquitous use. Stretching/shrinking surfaces have flow and heat transfer characteristics that are commonly used in engineering processes, including lamination and melt-spinning, polymer industries, and continuous casting. In view of this, Zaimi et al. [29] obtained the dual branch solutions of the 2D flow of a nanofluid along a permeable shrinking surface. Freidoonimehr et al. [30] analytically discussed the combined effects of heat generation and chemical reaction on a MHD flow past a stretching/shrinking surface. In another article, Uddin et al. [31] numerically studied the aspects of solar radiation on a nanofluid flow past a shrinking sheet in the presence of a slip effect. The impact of a heat source and chemical reactive species on a stagnation point flow of MHD fluid due to the shrinking surface was considered by Dash et al. [32]. Merkin and Pop discussed the exothermic surface reaction on a stagnation point flow over a shrinking surface. Further, Mahabaleshwar et al. [33] discussed the mechanisms of the heat and mass transport of a Casson fluid along a stretching/shrinking sheet. The stagnation point flow of a nanofluid past a quadratically shrinking surface in the presence of nanoparticles was investigated by Anuar et al. [34]. Recently, Khan et al. [35] analyzed the heat transport features of Carreau fluid through a shrinking sheet in the presence of Soret and Dufour effects. Mousavi et al. [36] theoretically and experimentally analyzed the behavior of a Casson fluid flow over a shrinking sheet in the presence of nanoparticles.
The effect of joule heating, viscous dissipation, and a nonuniform heat source/sink on magneto-Carreau fluid past a permeable stretching was of fundamental importance in our study, which was motivated by the aforementioned reference work and numerous prospective industrial applications of the problem. As a result, the goal of this research was to apply the findings of Akbar et al. [37] to a broader problem that incorporates the effects of viscous dissipation and Ohmic heating. Dual branch solutions were captured numerically [38] in this study to evaluate the impact of various flow factors appearing in the governing equations, which were also illustrated with the help of graphs. The proposed physical issue was numerically solved utilizing the shooting scheme, which is more computationally efficient. It was found that the solutions of the similarity equations have a dual branch of solutions in a certain range of shrinking parameters.

Mathematical Formulation
Let us consider the 2D boundary layer flow of magneto-Carreau fluid near a stagnation point driven by a permeable stretching/shrinking sheet, as shown in Figure 1. The effects of viscous and Joule dissipation, as well as a non-uniform heat source/sink, are considered in the energy equation. It is presumed that the velocity of shrinking/stretching is u w = λu w , while λ is a constant which represents the stretching or shrinking surface. The stretching of the stretching/shrinking sheet is vertical to the y-axis. The velocity of the shrinking sheet is u w = λu w , where λ is a constant with (0 < λ > 0), referring to the stretched/shrunk surface, respectively, and u e (x, t) is the stagnation point velocity. It is presumed that T w , T ∞ , v w (x, t), and B 0 represent the surface temperature, ambient temperature of fluid, mass transport velocity, and transverse magnetic field, respectively.
Coatings 2021, 11, x FOR PEER REVIEW 3 of 15 computationally efficient. It was found that the solutions of the similarity equations have a dual branch of solutions in a certain range of shrinking parameters.

Mathematical Formulation
Let us consider the 2D boundary layer flow of magneto-Carreau fluid near a stagnation point driven by a permeable stretching/shrinking sheet, as shown in Figure 1. The effects of viscous and Joule dissipation, as well as a non-uniform heat source/sink, are considered in the energy equation. It is presumed that the velocity of shrinking/stretching is , while is a constant which represents the stretching or shrinking surface. The stretching of the stretching/shrinking sheet is vertical to the y-axis. The velocity of the shrinking sheet is , where is a constant with (0 0), referring to the stretched/shrunk surface, respectively, and (x, t) is the stagnation point velocity. It is presumed that , , (x, t), and represent the surface temperature, ambient temperature of fluid, mass transport velocity, and transverse magnetic field, respectively.  Under the above assumptions, including the viscous dissipation and boundary layer concept, the governing equations of the proposed problem can be stated as (Pop et al. [28]; Bhattacharrya [39]): ∂u ∂x The associated boundary condition is: We imagine that v w , u w , u e , T w , B 2 (t), and q m (x, t) have the subsequent forms: where parameter c represents the unsteadiness of the problem; a is a positive constant; S denotes the mass transport parameter with (0< S >0), for suction and injection, respectively; A and B are parameters of space-dependent and temperature-dependent heat generation/absorption. It is to be illustrated that both A and B are positive to internal heat source and negative to internal heat sink. Considering Equation (5), we assume that it is possible to consider the following similarity factors: where ψ is the stream function and can be defined as: v= − ∂ψ ∂x , u = ∂ψ ∂y . Replacing (6) with Equations (2) and (3), we obtain the following ODEs Pr with associated boundary conditions, In the above equation, β, S, A, n, B, We, Pr, M, and Ec represent the unsteadiness parameter, suction parameter, space-dependent parameter, power law index, temperaturedependent local Weissenberg number, Prandtl number, Hartman number, and Eckert number, respectively. These parameters are defined as follows: The physical properties of importance including C f x and Nu x are characterized by local skin friction and the local Nusselt number, respectively. where τ w and q w represent the shear stress of wall and heat flux, respectively, via In dimensionless form: where x is the local Reynolds number.

Solution Method
The set of self-similar Equations (7) and (8) along with assisting boundary conditions (9) has been tackled numerically through the Runge-Kutta-Fehlberg method. The final system is reduced into a set of first-order ordinary differential equations and is altered into the initial value problem as: with initial conditions: 

Results and Discussion
The simulation results of Equations (7) and (8) conducted at the boundary condition (9) were analyzed for numerous values of the pertinent parameters β, We, S, n, A, Pr, B, M, Ec, and λ and the developmental condition in which the dual (first and second branch) solution can occur in the unsteady flow across a shrunk surface. In Figure 2, dual nature of solutions are observed for local skin friction Re 1/2 C f for distinct values of n. The values n = 2.5, 2, 1.5 represent the shear thickening behavior of Carreau fluid found in dispersions of highly condensed colloid particles. Its viscosity increases the shear loading, making it useful in protective and impact resistance applications. As the critical values of λ changes from λ c = (−5.1323, −5.4061, −5.942), it is observed that surface drag force increases for first solution and opposite pattern is seen for second solution for higher values of n. Figure 3 displays the impact of Re 1/2 C f with λ for S when the dual solution λ c < λ < 0, where critical values of λ are λ c = (−4.5731, −4.9906, −5.4061). From Figure 3, it can be illustrated that first solution is significantly increasing, whereas second solution diminishes for the higher values of S. Similar phenomena are noted, although we examined Figure 4, with Figure 3 for We. The influence of Re −1/2 Nu compared to λ for various values of Pr is described in Figure 5. In this figure, we have computed that the λ c is −4.99. In the heat transfer mechanism, the Prandtl number controls the relative thickness of the momentum and thermal boundary layers. It is noticed that when the values of Pr are small, the heat diffuses quickly in the case of the second solution. Thus, the heat transfer rate at the surface enhances with higher values of Pr in the case of the first solution. As can be shown in Figure 6 for both solutions, the rising Ec results are reduced in the heat transport rate Re −1/2 Nu. Further, the Eckert number indicates if the transition of momentum energy to heat energy has a major impact on the fluid flow and heat transfer. Therefore, higher values of Ec reduce the heat transfer rate for both solutions. The variation of velocity profile f (η) for different values of β is revealed in Figure 7. From this Figure, it is noted that the velocity profile reduces for both branches for improving the values of β. The impact of M on velocity profile f (η) is distributed in Figure 8. From this plot, it is stated that higher values of M reduces the velocity profile. This is due to the fact that M provides the Lorentz force, which slows down the movement of fluid, thereby decaying the thickness of the boundary layer. Figure 9 illustrates that the temperature of the non-dimensional characteristic, thus the rate of heat transport, also improves for the first branch and decays for the second branch with the enhancement of B. On the other hand, an inverse behavior of A on the heat transport is detected, as can be noticed from Figure 10. The variation in θ(η) is displayed in Figures 11 and 12 for higher values of β and M. These figures illustrate that the heat transport improves for both branches for large values of β and M. These outcomes are well developed and have not been replicated here, for simplicity.                          Table 1 displays the scheme of changing parameters on Re 1/2 C f . Now, the Re 1/2 C f for the first solutions improve with the increase in S and decays for the second solution; in addition, Re 1/2 C f for both solutions improve with the increases in We and n. Table 2 depicts the characteristics of numerous parameters on Re −1/2 Nu. The Nusselt number decreases for both solutions with the improvement of Ec and also increases with Pr. It is observed that increasing S, We and n have positive and incremental effect on skin friction coefficient.